时频分析在地震数据处理中的应用.doc

(完整word版)时频分析在地震数据处理中的应用
时频分析在地震数据处理中的应用
摘要：经典的傅立叶分析只适用于分析平稳信号，而不适用于非平稳信号。为了分析非平稳信号，我们采用时频分析方法。时频分析能够清楚的揭示信号的时变频谱特征，是对时变、非平稳信号进行分析与处理的有力工具。本文介绍了时频分析中常用的STFT和Wigner－Ville分布和S变换的特点，并利用Matlab对一个地震信号进行时频分析的实现。
关键词：时频分析；短时傅立叶变换（STFT）；S变换；Wigner－Ville分布
引言
在传统的信号处理领域，基于 Fourier 变换的信号频域表示及其能量的频域分布揭示了信号在频域的特征，它们在传统的信号分析与处理的发展史上发挥了极其重要的作用。但是，Fourier 变换是一种整体变换，即对信号的表征要么完全在时域，要么完全在频域，作为频域表示的功率谱并不能告诉我们其中某种频率分量出现在什么时候及其变化情况。然而，在许多实际应用场合，信号是非平稳的，其统计量（如相关函数、功率谱等）是时变函数。这时，只了解信号在时域或频域的全局特性是远远不够的，最希望得到的乃是信号频谱随时间变化的情况。为此，需要使用时间和频率的联合函数来表示信号，这种表示简称为信号的时频表示。
时频分析方法旨在通过构造一种时间和频率的密度函数，将一个一维的时间信号以二维的时间一频率函数形式表示出来，以揭示信号中所包含的频率分量及其随时间的变化特性。这使我们不但能够同时掌握非平稳信号的时域及频域信息，而且可以清楚地了解非平稳信号的频率是如何随时间变化的。
通过时频分析方法技术对地震信号进行分析处理可在获得地震信号的瞬时频率、瞬时相位、瞬时振幅等瞬时参数的同时获得时频谱图等重要时频域信息，实现对地震信号的边缘检测、属性提取等。
时频分析是非平稳信号分析处理领域的重要方法，时频分布的基本任务是建立一个函数，要求这个函数能够同时用时间和频率来描述信号的能量密度。如果
有了这样的一个分布，就可以计算某一确定的频率和时间范围内能量的百分率、
计算某一特定时刻的频率密度、计算该分布的整体和局部的各阶矩。即寻找一个联合密度函数P(t，f)，使p(t，f)=在时间t和频率f的强度，或者P(t，f)在时间t和频率f，在时一频单元内的部分能量。
地震信号的时频分析
一般将时频分析方法分为线性和非线性两种。典型的线性时频表示有短时傅叶变换(简记为STFT)、Gabor展开和小波变换(Wavelet Transformation，简记为WT)等。非线性时频方法是一种二次时频表示方法(也称为双线性)，最典型的是WVD(Wigner-Ville Distribution)和Cohen类。
采用的地震信号如下图所示：

图1 时间域的信号

图2 用FFT方法求取的能量谱
1、短时傅里叶变换（STFT）
传统的傅立叶变换只在频率域具备了局部分析的能力，而在时间域不具有这种能力。要得到满足同时获得时间和频率的局部分析能力的要求，一种最基本的方法就是:取出信号在所关心时刻附近的一小段，而忽略信号的其它部分，对其作傅立叶变换，即可得到这一特定时刻的频率分量。因为所取的时间长度与整个信号相比很短，所以将这种方法称为短时傅立叶变换(STFT)，它是时频分析中最简单的形式。短时傅立叶变换(STFT)的基本思想:用窗函数来截取信号，假定信号在窗内是平稳的，采用傅立叶变换来分析窗内信号，以便确定在那个时间存在的频率，然后沿着信号移动窗函数，得到信号频率随时间的变化关系，这就得到了我们所需要的时频分布。
STFT的物理意义在于，对于一定的分析时刻t, 可以视为信号s(t)在该时刻的“局部频谱”，从而整个变换的结果也就能揭示信号频谱的变化特性。短时Fourier变换的时频分辨率受制于窗函数的形状和宽度短时傅立叶变换的时间分辨率与分析窗函数的时间域宽度成正比，而其频率分辨率与分析窗的频宽成正比。从而，一个好的时间分辨率需要一个短的窗函数，而一个好的频率分辨率需要一个长的窗函数。因此，短时傅立叶变换不能同时兼顾时间分辨率和频率分辨率。当选用的窗函数为Gaussian函数时，该变换为Gabor变换。

图3 用STFT做时频分析
2、S变换
主要对S变换的定义、推导和特性进行详细阐述，S变换综合短时傅立叶变换和小波变换的优点，又避免了它们的不足:它与傅立叶变换有着直接的联系，具有无损可逆性;与短时傅立叶变换和小波变换一样，也是一种线性时频表示，因此不存在交义项的干扰;S变换具有多种分辨率，克服了短时傅立叶变换固定分辩率的不足;S变换中含有相位因子，这是小波变换所不具备的特性。总体来说，S变换是近几年发展起来的一种新的时频分析方法.S变换结合了短时傅里叶变换和小波变换的优点,具有相位信息,同时该变换与小波变换一样,其时频窗可以调节大小以适应非平稳信号的特点.S变换的这些优点,使得它在地球物理方面得到广泛的应用。
S变换首先是由Stockwell等人提出的，是以Morlet小波为基本小波的连续小波变换的延展。在S变换中，简谐波与高斯函数的乘机构成了基本小波，因为简谐波在时间域可以作伸缩变换，而高斯函数则进行伸缩和平移。ST也可以认为是CWT的“相位校正”。函数的S变换表示为：

其中：

可以用表示的傅里叶变换，其中和是相同的意义。

图4 S变换时频分析
3、WVD
Wigner-Ville变换是1932年由Wigner首次提出的，并应用于量子力学领域，后来Ville等人将其引入到信号分析处理领域。1966年，Cohen发现各种发现各种时频分析只是Wigner-Ville变换的不同形式，可以统一起来，成为Cohen类双线性时频分析。
信号的Wigner-Ville变换用公式表示为

式中是的解析信号，即

是实信号的Hilbert变换。
Wigner-Ville变换也可以用解析信号的频谱来表示

从上面两种不同形式的Wigner-Ville变换表达式中可以看出，式中不包含任何窗函数，从而避免了线性时频变换中实践分辨率和频率分辨率的相互牵制，难以兼顾的问题。因为Wigner-Ville变换的时间-带宽积可达到Heisenberg测不准原理(不确定原理给出的下界，故可以证明，没有任何一种时频变换方法的时间-频率分辨率及聚集性能出其右。
在Wigner-Ville变换表达式中，信号出现了两次，故称之为双线性时频变换。Wigner-Ville变换不是线性的，即两信号之和的Wigner-Ville变换不等于每一个信号的Wigner-Ville变换之和，其中多出了一个附加项。
令
则有

其中

式中前两项是自由项(auto terms)，第三项是交叉项(cross terms)。交叉项常常导致时频平面上出现伪影现象：交叉项是实的，混杂于自由项成分当中，并且其幅度是自由项成分的两倍；交叉项是震荡型的，每两个信号分量就会产生一个交叉项。交叉项的存在严重地干扰着对Wigner-Ville分布的解释，当信号变得非常复杂时，Wigner-Ville分布甚至变得毫无意义。

图5 WVD时频分析
4、小波变换
在短时傅立叶变换和Gabor展开中我们都使用了固定的时间窗函数，这就引出了时间分辨率和频率分辨率的概念，时间分辨率和频率分辨率是一对矛盾。根据海森堡的测不准原理，即时间窗函数的长度越长，频率分辨率就越高，而对于时间分辨率则越差。为了平衡时间分辨率和频率分辨率这个矛盾，可以采取对存在高频分量的部分采用高的时间分辨率和低的频率分辨率，而对于低频分量则采用高的频率分辨率和低的时间分辨率的方法，这就是多分辨分析的思想。 小波变换是一种在时间-尺度平面内，利用多分辨率分析思想分析非平稳信号的方法。所谓小波，就是一个满足容许条件的一个函数族

可以看出函数族是由窗函数在时间上平移b，在尺度上伸缩a，再乘上归一化因子后的结果，所以非平稳信号的连续小波变换定义为 :

其中是小波基函数的共轭。
将小波变换和短时傅立叶变换两者的基函数相比较，可以看出，小波变换基函数的尺度参数决定了小波变换的多分辨分析特性，即利用时间-尺度联合函数来分析非平稳信号的“变焦距”法，以达到分析信号局部特性的目的。
小波变换由于其本身分辨力的优良性能，因此一经提出，很快就成了非平稳信号分析和处理的一大热点，经过近 20 年的发展，小波变换取得了突破性的发展，形成了多分辨分析，框架和滤波器组三大完整丰富的小波变换理论体系。现在小波变换己经被广泛地应用在信号的奇异性检测、计算机视觉、图像处理、语音分析与合成等等诸多领域、在分形和混沌理论中也有了很多的应用。
声音信号分析：
1. 女声版江南style：

图6 女声版声音信号
图7 FFT所得振幅信息

图8 短时傅里叶变换时频分析

图9 S变换时频分析

图10 WVD时频分析

图11 小波变换时频分析

图12 Hilbert变换后小波变换时频分析

图13 Hilbert变换后WVD时频分析
2. 男声版江南style：

图12 男声版声音信号

图13 FFT所得振幅信息

图14 短时傅里叶变换时频分析

图15 S变换时频分析

图16 WVD时频分析

图17 小波变换时频分析

图18 加入Hilbert变换小波变换时频分析

图19 加入Hilbert变换WVD时频分析
对比男生与女生版两段相似的声音信号发现女生的频率段基本在0.35kHz左右，而男生版的声音信号则有一段明显的低频段。再者由于声音信号内存在的其他干扰，可能是的信号的频率与能量等有差异，但是最基本的频率信息还是可以看出来的。
另外，交叉项是震荡型的，每两个信号分量就会产生一个交叉项。交叉项的存在严重地干扰着对Wigner-Ville分布的解释，当信号变得非常复杂时，Wigner-Ville分布甚至变得毫无意义。就像上面分析的声音信号，除了看出简单的对称性之外，其他的分析意义基本没有。
几种不同的Wigner-Ville变换
原始信号： x1=cos(2*pi*200*t).*(t>=0 & t<=0.8);
x2=cos(2*pi*50*t).*(t>=0.2 & t<=1);
x=x1+x2;

图20 原始信号

图21 三种WVD时频分析

图22 WVD时频分析三维显示

结论
由于地下介质的复杂性，地震信号是一种非平稳随机信号，时频分析是分析处理地震信号的一种非常有效的方法。前人在此领域已做了大量工作，将很多时频分析的方法应用在地震数据处理解释中。本文对各种时频分析方法的比较研究和分析其各自特点。并将时频分析的方法用于地震信号分析处理中，结果表明:通过时频分析方法对地震剖面进行处理，可以得到瞬时地震属性参数(包括频率类瞬时属性参数、瞬时振幅谱)，这些新的地震剖面对地质解释有一定的辅助作用。另外后面声音信号的分析说明，时频分析对于频率有很好的敏感性，对于区分不同频带的信息有着非常好的分辨率。像不同的乐器，相同物体不同的激发声音方式等，都可以通过时频分析进行解释。
时频分析以联合时频分布的形式来表示信号的特性，具有很多的优点:
1）它克服了傅里叶分析时域和频域完全分离的缺陷，将时频两域联合起来对信号进行分析，能同时考虑到两个方面的性能。
2）弥补了信号的时间能量密度和频谱能量密度不能充分描述信号的物理特性的缺陷。
3）在时频相平面上，可以精确地定位在某一时刻出现了哪些频率分量，以及某一分量出现在哪些时刻。
4）对于不同情况的信号，通过对核函数施加一些约束条件，就可以设计符合期望性能的时频分布，来满足处理信号的要求。
5）为非平稳信号的分析提供了有效的工具，为信号的分析开辟了新的途径。
时频分析虽然具有很多优点，但同时也具有不少缺点:
1）由于双线性形式的时频分布是非线性的，使得两个信号和的时频分布已不再是两个信号各自分布的和，即存在交叉项。
2）在时频域进行去噪时计算复杂，不易实现。
3）时频分布虽然反映了信号的能量分布，但不能用信号的“瞬时能量”来解释某一时刻或某一频率处的时频分布。
现实中很多信号，比如语音信号，都是时变非平稳的，时变非平稳特性是现实信号的普遍规律，联合时频分析技术正是应现实的科学和工程应用需求而产生和发展起来的。对于许多信号，仅用时域或频域里的各种方法去分析往往不能揭示信号内部的局部特征和信息，而时频分析作为一种能将频谱随时间的演变关系明确表现出来的新手段，自然更符合实际应用的需要。
参考文献
1. 邹红星,周小波,李衍达. 时频分析:回溯与前瞻[J]. 电子学报,2000,09:78-84.
2. 刘葵,刘招君,朱建伟,张峰. 时频分析在石油地球物理勘探中的应用[J]. 世界地质,2000,03:282-285.
3. 刘丽娟. 时频分析技术及其应用[D].成都理工大学,2008.
4. 陈斌. 时频分析及在地震信号分析中的应用研究[D].成都理工大学,2007.
5. 于振江. 一种基于Matlab的语音信号采集与分析系统设计[J]. 科技情报开发与经济,2012,12:112-114.
6. 陈家焱,陈冬娇,张达响. 基于Matlab的声音信号采集与分析处理[J]. 计算机与现代化,2005,06:91-92+96.
7. 题原,张劲松. 基于MATLAB的语音信号采集和分析系统的可视化设计[J]. 齐齐哈尔大学学报,2006,06:43-46.
8. 蒋济同,唐世振. 基于MATLAB的振动信号采集与分析系统的研究[A]. 中国力学学会、中国振动工程学会、中国航空学会、中国机械工程学会、中国宇航学会.第九届全国振动理论及应用学术会议论文集[C].中国力学学会、中国振动工程学会、中国航空学会、中国机械工程学会、中国宇航学会:,2007:5.
9. 宋杨洁. 基于LabVIEW与MATLAB的语言信号的采集与分析[D].武汉理工大学,2012.
10. 周渊,王炳和,刘斌胜. 基于MATLAB的噪声信号采集与分析系统的设计[J]. 电声技术,2004,07:52-54.
附：
1、时频分析基于Matlab的实现
clc;
clear;
nfft=512;
fid=fopen('xinhao1.txt','r');
fclose(fid);
%x=u(:,200);
figure(1);
tfrstft(u);
figure(2);
tfrsp(u);
figure(3);
tfrwv(u);
2、 不同的Wigner-Ville变换
clear all;
N=2048;
fs=2048;
n=0:N;
t=n*1/2000;
x1=cos(2*pi*200*t).*(t>=0 & t<=0.8);
x2=cos(2*pi*50*t).*(t>=0.2 & t<=1);
x=x1+x2;
figure(1);
plot(t,x);
title('原始信号');
xlabel('时间 t');
ylabel('幅值 A');
grid on;
figure(2);
subplot(311);
[tfr,t,f]=tfrwv(x',1:N,N);
contour(t/fs,f(1:length(f)/2)*fs,abs(tfr(1:length(f)/2,:)));
title('Wigner-Ville分布时频图');
xlabel('时间 t');
ylabel('频率 f');
grid on;
subplot(312);
[tfr,t,f]=tfrpwv(x',1:N,N);
contour(t/fs,f(1:length(f)/2)*fs,abs(tfr(1:length(f)/2,:)));
title('伪Wigner-Ville分布时频图');
xlabel('时间 t');
ylabel('频率 f');
grid on;
subplot(313);
[tfr,t,f]=tfrspwv(x',1:N,N);
contour(t/fs,f(1:length(f)/2)*fs,abs(tfr(1:length(f)/2,:)));
title('平滑伪Wigner-Ville分布时频图');
xlabel('时间 t');
ylabel('频率 f');
grid on;
figure(3);
mesh(t/fs,f(1:length(f)/2)*fs,abs(tfr(1:length(f)/2,:)));
title('三维时频图');
xlabel('时间 t');
ylabel('频率 f');
zlabel('幅值 A');