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

(完整word版)时频分析在地震数据处理中的应用时频分析在地震数据处理中的应用摘要：经典的傅立叶分析只适用于分析平稳信号，而不适用于非平稳信号。为了分析非平稳信号，我们采用时频分析方法。时频分析能够清楚的揭示信号的时变频谱特征，是对时变、非平稳信号进行分析与处理的有力工具。本文介绍了时频分析中常用的STFT和WignerVille分布和S变换的特点，并利用Matlab对一个地震信号进行时频分析的实现。关键词：时频分析；短时傅立叶变换（STFT）；S变换；WignerVille分布引言在传统的信号处理领域，基于 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、WVDWigner-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.2 & t=0 & 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);