几种时频分析方法比较

Dec. 2006

几种时频分析方法比较

陈雨红1, 杨长春1, 曹齐放2, 李波涛1, 尚永生1

(1.中国科学院地质与地球物理研究所,北京100029; 2.胜利油田物探公司地震勘探研究所开发地震室,东营257100)

摘 要 地震信号由于各种原因,往往是一种非线性、非平稳信号,基于平稳信号理论的常规傅立叶变化方法不能

刻画任一时刻的频率成分.时频分析能同时保留时间与频率信息,目前已经出现了很多时频分析方法.本文介绍了Hilbert 变换、H ilbert -Huang 变换、正弦曲线拟合、雷克子波匹配、短时傅立叶变换、小波变换、S 变换以及Co hen 类这八种方法,并从时间分辨率、频率分辨率,以及对多频率成份信号适应能力等各方面阐述了各种方法的优缺点,对其中的一些方法结合了理论记录进行了试算,进一步阐述了这些方法的长处和不足之处.关键词 时频分析,Hilbert 变换,正弦曲线拟合,S 变换,H ilber t -H uang 变换中图分类号 P631 文献标识码 A

文章编号 1004-2903(2006)04-1180-06

The comparison of some time -frequency analysis methods

CH EN Yu -hong 1, YANG Chang -chun 1, CAO Q-i fang 2, LI Bo -tao 1, SH ANG Yong -sheng 1

(1.Institu te of Geolog y and Ge op hy sics ,Chinese A cad emy of S ciences ,Beij in g 100029,China ;

2.S heng li exp lor ation r esearc h institute ,Dong ying 257100,China )

Abstract T he seismic data usually is a non -linear non -statio nar y ser ies fo r many reasons,so the F our ier transfor m based on stat ionary signal can no t sho w the new s of fr equency at sometime.T he time -frequency analy sis can give the new s of time and f requency at the same time,no w there have many methods of time -frequency analy sis.T his paper intro duces eig ht metho ds o f t ime -frequency analysis,such as H ilber t t ransfo rm,H ilbert -H uang tr ansfo rm,sine curv e fitting ,R icker w avelet matching ,sho rt time Fo urier t ransfo rm,w avelet t ransform,S tr ansfo rm and Co hen metho ds.T his pa per ex plains the merits and defects o f those met ho ds fro m the t ime r eso lutio n,frequency r eso lutio n and the adaptiv e abilit y of mult-i frequency signals.M o reo ver,this paper explains the advantages and short co mings of so me metho ds deeply by some test cases.Keywords

time -f requency analysis,hilbert t ransfo rm,sine curv e fitting ,S tr ansfor m,hilbert -huang transfor m

收稿日期 2006-02-10; 修回日期 2006-05-20.

基金项目 中科院院知识创新工程重大项目(kzcx1-s w -18)资助.

作者简介 陈雨红,男,1979年4月生,2004年获石油大学(北京)硕士学位,现为中科院地质与地球物理研究所博士生.(E -mail:yuhong.

chen @sohu.com )

0 引 言

在地震勘探过程中,由于地层吸收,孔隙流体等各种原因,地震信号往往具有非线性、非平稳特征,即实际地震资料的频率成分是随时间变化而变化的,因而基于平稳信号分析理论的常规傅立叶变换

和功率谱估计方法就失去了物理意义,也就不适合于表征某一时间点的频率成分分布情况.因此需要把整体谱推广到局部谱中来,目前已经发展了很多时频分析方法.

时频分析方法在地震勘探中有很多应用,如:瞬

时属性的提取、描述地震相[1]、判定P 波到达时[2]、计算泥岩体积[3]

、地层吸收补偿[4,5]

、信号识别

[6]

、压

制随机噪音[7]、相干体[8]、薄层识别[9]

等.

本文分析对比了八种时频分析方法:H ilbert 变换、H ilbert -H uang 变换、正弦曲线拟合、雷克子波匹配、短时傅立叶变换、小波变换、S 变换以及Cohen 类双线性变换方法.主要从时间分辨率、空间分辨率以及多频率信号识别能力等方面比较详细地阐述了各种方法的优缺点.

4期陈雨红,等:几种时频分析方法比较

1 Hilbert 变换

Tarner [10]等应用复地震道分析方法计算瞬时振幅、瞬时相位和瞬时频率,该方法中设地震道为,通过H ilber t 变换得到复地震道

z (t)=s r (t)+j s i (t)=s(t)+j s( )

t - d .

(1)

则瞬时振幅为

a(t)=s 2r (t)+s 2

i (t).

(2)

瞬时相位为

(t)=tan -1

s i (t)s r (t)

.(3)

瞬时频率为

f (t)=

12 d (t)

d t

.

(4)

该方法假定信号为单频率成分信号,通过H ilbert 变换能够很好的给出信号瞬时特征,有很高的时间分辨率,但是信号两边界处难以处理,往往误差较大,而且也会传播到其它地方,影响频率计算精度.如图1b,该方法在频率间断点处也能很好的反映瞬时特征,但是在两端频率值抖动较多.

图1a 单频信号,0.1~0.2s 是频率为40Hz 的正

弦信号,其余为频率为20H z 的正弦信号

Fig.1a Mono -f requency signal,the sig nal betw een 0.1s to 0.2s is s ine -cur ve of 40H z

,a nd other s is sine -cur ve of 20H z

图1b 图1a 信号利用Hilbert 变换得到的瞬时频率Fig.1b The instantaneous f r equency of the signal

(Fig.1a)by H ilber t tr ansf or m

该方法对于多频率成分信号情况,由于H ilbert

变换没有实际物理意义,因而会出现频率尖脉冲情况和负频率情况.如图2b,通过该方法计算的瞬时频率没有任何物理意义.

图2a 多频信号,该信号是三个信号叠加而成

第一个信号:振幅为2,频率为80H z (0~0.3s);第二个信号:振幅为1,频率为20H z (0~0.15s);第三个信号:振幅为1,频率为40H z (0.15s ~0.3s).Fig.2a Multi -f r eq uency sig nal ,the s ignal comp osed by thr ee sig nals.the f ir st sine -curv e sig nal (0~0.3s)is amp litude of 2and f r equency of 80H z ;the second (0~0.15s)is amp litude of 1and f requency of 20H z ;the last(0.15~0.3s)is amp litud e of 1and f r equency of 40H z

图2b 图2a 信号利用H ilber t 变换得到的瞬时频率Fig.2b T he instantaneous f r equency of the si gnal

(Fig.2a)by H ilber t trans f or m

2 H il bert -H uang 变换

N.E.H uang (1998)等人首先介绍了H ilbert –H uang 变换,它是为了解决H ilbert 变换不能处理多值频率而发展起来的,在国内石春香[11]等首先

用于地震资料时频分析中,该方法的关键是经验模态分离(EM D),它认为任何复杂的时间序列都是由一些相互不同的、简正的固有模态函数(IMF )组

1181

地 球 物 理 学 进 展21卷

成,因此,可以从复杂的时间序列直接分离出从高频到低频的若干阶基本时间序列,即IM F ,然后对每条IMF 利用H ilber t 变换计算其瞬时属性

.

图3a

多频信号,该曲线由多频率成分组成

Fig.3a Mul ti -f r equency signal ,the cur v e comp osed by mul ti -f r equency comp onent sine -cur

ve

图3b 图3a 中信号用常规H ilber t -H uang

变换得到的时频分布情况

Fig.3b The time -f r equency d is tr ibution of the signal (Fig.3a)by H ilber t -H uang tr ansf or

m

图3c 图3a 中信号用改进后的H ilber t -H uang

变换得到的时频分布情况

ig.3c T he time -f requency distr ibution of the s ignal (Fig.3a)by imp r oved H il ber t -H uang tr ansf or m

该方法计算的瞬时频率有很高的时间分辨率和较高的频率分辨率,然而,这里有三个主要缺点:在低频区会产生不合理的频谱特征;首先得到的几个IM F 会包含太宽的频率范围(前面几个IMF 往往反映高频信息),因而不能满足H ilber t 变换要求的信号为单频率;不能分离出信号为低能量成分.

Peng

[12]

等对该方法进行了改进,它利用小波包

来进行经验模态分离,更好的保证了每个IMF 具有窄的频带范围,同时对每个IMF 与原始信号做互相关,筛选掉互相关程度弱的信号,可以筛选掉不合理的低频成分,做互相关时进行了归一化处理,这也可以保留了能量低但相关程度高的频率成分.如图3b,

常规方法在高频区很混乱,实际上就是IMF 频带太宽,不能保证是单一频率成分造成的,而且很多弱能量信号也被掩盖掉了;如图3c,改进后在高频区频率分辨率也很高,而且弱能量信号也能体现出来.

3 正弦曲线拟合

Choi [13]

等首先使用正弦曲线拟合的方法来分

析语音信号,H ar dy 等[3]将该方法应用于地震信号瞬时频率的估计.该方法是利用一组不同频率的正弦曲线来拟合一时窗内信号,取拟合误差最小的正弦曲线为该信号的最佳拟合曲线,同时可以给出频率、振幅及相位信息,取拟合的剩余曲线重新拟合,

直到拟合误差小到一定程度为止,这样就可以得到一组频率、振幅及相位,因此该方法可以用于多频率信号.但是该方法分析信号必须在一个时窗内进行,时窗内有滤波效果,因而该方法时间分辨率没有H ilber t 变换高;而且不同频率成分的信号会互相干扰,所以在提取一个频率成分的信号时会破坏其它频率信号,只有窗长选择半个周期及其整数倍附近才对其它频率影响较小,因而在不同频率进行拟合时最好采用选取不同的窗长的方法进行拟合,而H ar dy 采用的算法是固定窗长,也就造成了时频曲线振荡情形,另外该方法中采用固定时窗,在频率间断点处会出现不合适的频率成分,因而建议在分析每个点频率时滑动时窗,避免固定时窗方法中时窗内前后采样点频率不同相互影响.但是这加大了不

少计算量.

如图4、图5,该方法在频率间断点周围存在一个斜

坡,这就是受窗口的影响,时间分辨率没有Hilber t 变换方法高,而且在间断点出现很多散点,这就是窗口内存在两种频率成分,它们相互影响造成虚假频率成分的出现,但是该方法能够分辨多值频率.

1182

4期陈雨红,等:

几种时频分析方法比较

图4 图1a 中信号利用正弦曲线拟合

方法得到的时频分布情况

Fig.4 T he time -f r equency distr ibution of the s ignal

(Fig.1a)by sine -cur ve f itting metho

d

图5 图2a 中信号利用正弦曲线拟合

方法得到的时频分布情况

Fig.5 The time -f r eq uency dis tr ibution of the sig nal (Fig.2a)by s ine -cur ve f itting method

4 雷克子波匹配

Liu

[14]

首先给出了基于雷克子波分解方法进行

时频分解,由于地震记录比较符合雷克子波情况而不是正弦曲线,因而该方法有可能比正弦曲线拟合方法能得到更高时间精度和频率精度的时频分布特征.然而在这个算法中沿用H ilbert 变换得到的瞬时频率、瞬时相位和瞬时振幅来计算雷克子波的主频、延时及振幅,这对于在薄互层情况下出现的多个子波叠加情况,也就是出现多值频率情况,用H ilber t 变换得到的频率没有任何物理意义,因而它的频率、振幅及相位均不反映实际情况,所以用该算法计算的瞬时属性在多值频率下可信性很低.但是如果对频率及延时完全按照扫描拟合得到,那么它的时间分辨率和频率分辨率都比较高,但计算量太大,目前也没发现有人这么做.不过该思路也许会有很大的

发展前景.

5 短时傅立叶变换

短时傅立叶变换是通过滑动时窗来计算其频谱,因而它的时间分辨率和频率分辨率受H eisenber g 测不准原理约束,也就是利用短窗口有较高的时间分辨率,但是频率分辨能力差;当利用长窗口时有较高的频率分辨率,但时间分辨能力就弱.如图6a,利用的是窗长较大的窗口,频率分辨能力还可以,但时间方向上就分辨不出来了;如图6b,利用的是窗长较小的窗口,时间方向上分辨能力

还可以,但频率方向上就难以分辨了.

图6a 利用长窗长ST FT 变换得到的时频分布Fig.6a The time -f r equency distribution by short

t ime f r equency analy sis of long windth window

图6b 利用短窗长ST FT 变换得到的时频分布Fig.6b The time -f requency distribution by shor t t ime f r equency analysis of short windth window

6 小波变换

Chakr aborty [15]等阐述了利用小波变换的方法进行时频分析,它是一种多尺度(MRA )方法,对不同频率用不同尺度进行分析,它能给出比较好的时间精度,在低频区有很好的频率精度,而在高频区频率分辨能力较弱,但是在实际应用中还可以满足要

1183

地 球 物 理 学 进 展21卷

求.然而小波变换有一个致命的弱点,小波变换是沿着时间方向绝对平移暗含了调制沿着包络方向传播,这样造成了相位信息仅仅是局部的,会失去其物理意义造成难以解释.

7 S 变换

Stockw ell [16]对S 变换进行了详细阐述,Pinnegar

[18]

将该方法推广,并应用于地震信号时频

分析中来,它是短时傅立叶变换和小波变换的组合,其一维连续正变换为:

S( ,f )= +

- h(t)

|f |2

e -

f 2(t - )22

-j 2 f t d t .(5)

该式中比连续小波变换乘以了一个相位项

e j 2 f

(6)它解决了小波变换相位局部化问题,

它和小波变换一样是线性算子,对多频率信号没有双线性类型变换那样的交叉项干扰问题.

图7a 图1a 信号用S 变换得到的时频分布Fig.7a T he time -f requency distr ibution of the s ignal (Fig.1a)by

s tr ansf o rm

图7b 图2a 信号用S 变换得到的时频分布Fig.7b The time -f r equency d is tr ibution of the s ignal (Fig.2a)by s tr ansf o rm

如图7a,7b 所示,利用s 变换能够分辨多频率信号,但由于它实际上隐含存在一个窗口,因而在频率间断的地方频率刻画不准确,而且它也难以准确刻画振幅值,在峰值(谷值)与过零点地方振幅值相差较大.

8 Cohen 类方法

Cohen 类方法本质上是一种双线性变换,它是对信号定义瞬时相关函数.Wig ner -Ville 分布(WVD )定义为:W VD x (t,f )= +

- x t +

2

x *

x -

2

e -j 2 f

d .(7)

相应的模糊函数(AF )定义为:

A F x ( ,v)= + - x t +

2

x *

x - 2

e j 2

f d .(8)

可以将W VD 改写为如下形式:

W VD x (t,f )= +

-

+

-

AF x ( ,v)e -j 2 (tv + f )

d v d t .

(9)

该方法对于多分量信号或者具有复杂调制规律的信号,信号间的相互作用产生交叉项干扰,它降低了信号时频分布的分辨率,模糊了信号的原始特征,使时频分布难以解释.人们对AF 或WVD 进行了各种平滑的改进方法,它们可以以统一的形式表示出:p (t,f )= +

- +

- A F x ( ,v) ( ,v)e -j 2 (tv + f )

d v d t .(10)

( ,v)为核函数,当

( ,v)=h( ).(11)

得到的是伪Wigner -Ville 分布(P W VD ),它是对变量 加窗函数来达到减小交叉项的目的,也可以对方向也可以加窗平滑.当

( ,v)=e -(2 v )2

.(12)

时得到Choi -Williams 分布(CWD );

当

( ,v)=

sin ( v )

.(13)

时得到Bor n -J or dan 分布(BJ D );当

( ,v)=co s ( v).

(14)

时得到Mar genau -H ill 分布(MH D );当

( ,v)==g ( )| |sin (

v) v

.(15)

1184

时得到锥形核分布(CK D).

虽然采用核函数进行了平滑,可以减小交叉项的影响,然而这也损害了W VD的许多有用的特性,信号响的时频聚焦性也会有所降低,而且也不能完全消除交叉项的影响.

9 结 论

通过分析对比这八种时频分析方法,可以得出以下结论:

(1)H ilber t变换对于单频率信号有很高的时间分辨率和频率分辨率,但对于多频率成分信号它失去了物理意义,因而其瞬时属性就会没有任何意义,造成解释的困难.

(2)H ilber t-H uang变换对于多频率成分信号按照理论也有很高的频率分辨率和时间分辨率,然而它的核心技术经验模态分离方法还很不成熟,它往往会在高频区频率分辨率不高,在低频区也会出现不合理的频率成分,容易掩盖低能量的频率成分.

(3)正弦曲线拟合方法有比较高的时间分辨率和频率分辨率,适应多频率成分,但是它在提取多频率成分时,各频率成分可能相互影响,而且它对于频率有间断地方受窗口的影响会出现虚假频率成分,此外该方法计算量很大.

(4)雷克子波匹配现有算法由于受Hilbert变换方法影响不适应多频率成分信号.

(5)短时傅立叶变换方法在利用长窗口时,频率分辨率较高,但时间分辨率低;短窗口时,时间分辨率较高,但频率分辨率低.

(6)小波变换方法有较高的时间分辨率,低频区有较高频率分辨率,高频区频率分辨率较低,但其相位局部化,造成各频率相位基准不一,从而造成解释的困难.而S变换对相位进行了校正.

(7)W VD方法有很好的时频聚焦性,但受交叉项干扰的影响,它的各种平滑改进方法能一定程度上消除交叉项干扰影响,但又降低了时频聚焦性.

参 考 文 献(Ref erences):

[1] Philip p e S,Guy D.Seismic sequence analysis and attribu te

ex tr action using qua dra tic time-f r equ ency

rep re sentations[J].Geop hy sics2001,66(6):1947~1959.

[2] Rober t P C,La lu M.The S-transf or m w ith w indow s of

arbitray and v ar ying shap e[J].Geop hysics,68(1):381~

385.

[3] Har dy H H,Ric hard A,Beier,Jonathan D.Gaston.

Frequency e stimate s of seismic traces[J].Geoph ysic s,2003,

68(1):370~380.

[4] 白桦,李鲲鹏.基于时频分析的地层吸收补偿[J].石油地球物

理勘探,1999,34(6):2~8.

[5] 李鲲鹏,李衍达,张学工.基于小波包分解的地层吸收补偿

[J].地球物理学报,2000,43(4):542~549.

[6] 高静怀,满蔚仕,陈树民.广义S变换域有色噪音与信号识别方

法[J].地球物理学报,2004,47(5):869~875.

[7] 金雷,李月,杨宝俊.用时频峰值滤波方法消减地震勘探资料

中随机噪声的初步研究[J].地球物理学进展,2005,20(3):

724~728.

[8] 王西文,杨孔庆,周立宏,王娟,刘洪,李幼铭.基于小波变

换的地震相干体算法研究[J].地球物理学报,2002,45(6):

847~853.

[9] 高静怀,陈文超,李幼铭,田芳.广义S变换与薄互层地震响

应分析[J].地球物理学报,2003,46(4):526~532.

[10] Taner M T,K oehler F,Sh erif f R E.Comp le x se ismic

trace analysis[J].Geop hysics,1979,44(6):1041~1063.

[11] 石春香,罗奇峰.时程信号的H ilbert-H uang变换与小波分

析[J].地震学报,2003,25(4):398~405.

[12] Peng Z K,Pete r W,Tes F,Ch u L.An imp r ov ed

H ilber t-H uang transf or m and its app lic ation in v ibration

signal analysis[J].Journal of S ound and Vibration.2005,

286:187~205.

[13] Ch oi A.Real-time f unda me ntal f requency estimation by

least-square f itting[J].IE EE T rans.Sp eech A ud io

p rocessing.1997,5,201~205.

[14] Liu J L.Time-Fr equency d ecomp osition based on R ic ker

w av elet[J].SEG E xp and ed A bstarc ts,2004,23:1933. [15] Chakr aborty A,Okay a D.Frequency-time dec omp osition of

seismic data using w av elet-base d method s[J].Geop hysics,

1995,60(6):1906~1916.

[16] Stockw e ll R G,Mansinha L,L owe R P.Locatiz ation of the

c omp le x spe ctr um:the S-tr ansf orm[J].IEE E T ransactions

on S ignal Pr oce ssing,1996,44(4):998~1001.

[17] Por tniag uine O.Inv erse spe ctr al Decomp osition[J].SE G

E xp and ed A bstarc ts,2004,23:1786.

[18] Saba dell F J.Sp ec tr al localiz ation of seismic data w ith a

p hase c orr ected w av e let tr ansf orm[J].SE G Ex pand ed

A bstarc ts,2001,20:595.

[19] Arthur E,Bar nes.The calculation of I nstan taneous

f requenc y and instan taneous band w idth[J].Geop hysics,

1992,57(11):1520~1524.

[20] 高静怀,汪文秉,朱光明.小波变换与信号瞬时特征分析[J].

地球物理学报,1997,40(6):821~832.

[21] 殷文,印兴耀.基于MP I的时频分布的改进及应用[J].地球

物理学进展,2005,20(1):165~169.

1185