第6卷第6期2009年12月 CHINESEJOURNALOFENGINEERINGGEOPHYSICSVol16,No16Dec1,2009
文章编号:1672)7940(2009)06)0740)06doi:1013969/j1issn11672-7940120091061015
Riker子波匹配追踪算法及其改进屈念念,刘江平,李家斌(中国地质大学地球物理与空间信息学院,武汉430074)
作者简介:屈念念(1986-),女,硕士研究生,主要从事地震数据处理方法研究。E-mail:q120080265@yahoo.cn刘江平(1957-),男,教授,博士生导师,主要从事地震勘探的科研与教学工作。E-mail:liujp@cug.edu.cn
摘 要:匹配追踪(MatchingPursuit)算法的基本思想是基于信号的可分解和重构,是在一个确定的函数集
合中自适应地选择一些函数来表示一个信号的计算过程,函数集合中的每个函数都称为原子。本文利用奇异值分解对传统的匹配追踪算法进行了改进,提高收敛速度、计算速度以及重构精度,并将得到的时频分布与其他方法进行对比,测验结果证明了改进算法的高效性和有效性。关键词:匹配追踪算法;最小二乘;奇异值分解;时频分布
中图分类号:P631文献标识码:A收稿日期:2009-11-11
RevisedMatchingPursuitAlgorithmUsingRickerWavelets
QuNiannian,LiuJiangping,LiJiabin(InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China)
Abstract:ThebasicideaofMatchingPursuitalgorithmisbasedonthepossibilityofdecom-positionandreconstructionofseismicdata.Itisaprocessinwhichaseriesoffunctionscanbechosenfromacertainfunctionsettorepresentagivensignal,andtheeachandeveryfunctionwithinthefunctionsetiscalledasanatom.Thispaperusessingularvaluedecom-positionmethodtorevisematchingpursuitalgorithm,andimprovesthespeedofconver-genceandcomputation,andthencomparesthetime-frequencydistributionoftherevisedmatchingpursuitwithothermethods.Theresultprovestheeffectivenessoftherevisedmatchingpursuitalgorithm.Keywords:matchingpursuitalgorithm;least-square;singularvaluedecomposition;time-frequencydistribution
匹配追踪(MatchingPursuit)算法是在一个确定的函数集合中自适应地选择一些函数来表示一个信号的计算过程,函数集合中的每个函数都称为原子。其核心思想是将信号表示为一系列与信号局部结构特征最佳匹配的时频原子的线性组合,然后求各时频原子的时频分布并将其叠加,得到信号的时频分布。MP算法虽然计算复杂度较高,但有很多优良的性质,如对信号自适应的灵活表达,这是传统的傅立叶变化或小波变换[1~4]所无法比拟的,另外还具有较高的时间)频率分辨率,暂态结构的局部自适应性,信号结构的参数表示更加灵活等优点。MP算法一经提出,便很快被应用于地球物理的地震信号处理领域。在1996年应用MP分解算法对压缩的地震信号进行Kirchhoff偏移计算;在2003年独立多分辨率分析和MP算法对计算量和数据进行统计与分析;在2004年以Ricker子波为原子对地震信号进行时频分解;在2005年采用Morlet小波为原子对地震信号进行MP分解等。MP算法是一种重复迭代逼近的贪婪算法,因此其核心问题就是如何建立有效的原子库,并快速地检索出匹配的原子,提高检索效率,加快计算速度。DurkaPJ和QianS提出了在原子生成的层面,通过借助于前一次原子的生成结果来便捷地进行当前原子的计算生成,而原子字典的索引方式则是先依照离散化尺度参数,再根据尺度参数来确定其他原子参数的计算方式。尹忠科等将最相关原子的匹配与遗传算法的整体寻优过程联系起来以实现快速定位最佳匹配原子的目的,从而提高匹配逼近过程的计算效率。WangYH从瞬时频率和瞬时相位的角度出发,通过计算原始信号的瞬时频率和瞬时相位[5~9]来粗略地确定时频原子的频率参数和相位参数,然后在相对集中的字典子集中进行最佳原子的匹配搜索,同样也提高了计算速度。本文在上述学者研究工作的基础上对匹配追踪算法进行了改进,提高收敛的速度和计算速度,并将得到的时频分布与其他方法进行对比,测验结果证明了改进算法的高效性和有效性。
1 匹配追踪基本原理MP算法的核心思想就在一个确定的函数集合中挑选最能体现信号特征的一系列函数,其中每个函数称为原子[10],MP算法原子库的类别有许多种,如:Gabor原子,ChirpLet原子,FMmlet原子,Ricker子波原子,阻尼正弦函数原子等。假设D为进行信号分解的函数集合,原始信号为u,长度为N。MP算法是通过把函数u垂直投影到函数集合D的元素上来进行重复的迭代估算,在进行一次迭代后,函数u可以表示为:u=3u,gC4gC
+Ru(1)
其中gCID,Ru是把u函数在函数集合D上进行垂直投影后的残余矢量,为了使残余矢量尽可能的小,就必须使内积计算项|&尽可能大。很显然,gC与Ru是正交的,因此:+u+2=3u,gC42++Ru+2(2)假设R0u=u,且进行了n轮的迭代(n\0)得到残余矢量Rnu,此时在函数集合D选择一个原子,使其匹配逼近这个残余矢量Rnu,因此:Rnu=3Rnu,grn4grn+Rn+1u(3)其中Rn+1u就是进行了n+1次迭代得到的残余矢量,gCn是第n次迭代时选择的原子。MP算法就是这样的一个重复迭代的过程,如果重复这样的迭代过程m次,即可将u表示为如下的形式:u=Emn=13Rn-1u,grn4grn+Rnu(4)因此经过m次迭代分解计算后,原始信号u可用m个原子的合成来近似表示,其误差为第m次迭代计算后的残余矢量。匹配追踪算法的每次迭代时都按照一定的索引扫描整个函数集合,搜索出与当前信号最大相关的原子参数,然后剔除最相关原子的能量形成残余信号,残余信号再进入下一轮迭代计算,直到迭代循环终结或残值信号满足设定的阀值。2 Ricker子波匹配追踪算法原理 Ricker子波的时间域表示为:XR(t,fj)=(1-2P2f2jt2)exp(-P2f2jt2)(5)利用Ricker子波原子进行匹配追踪分解时,每次迭代都剔除与地震道中最大相关的Ricker子波[12~16]。迭代过程一直持续到剩余的地震道幅度小于给定的阈值为止。此过程便把地震道分解为Ricker子波原子的线性组合,其结果是对应原始地震道信号的一系列具有不同到达时刻和振幅的Ricker子波。每一道有限带宽的地震信号u(t)都可以表示为Ricker子波的线性组合:u(t)=Ejaj*w(t-tj,fj,Hj)+Noise(6)其中Aj,tj,fj,Hj为分解出的第j次Ricker子波w的振幅,中心时间,波峰频率和相位。tj是通过u(t)的瞬时包络局部极大值来估计的,即通过对信号做Hilbert变换求信号的瞬时包络,其中tj可粗略估计为包络峰值处对应的时间t,瞬时频率favg则为tj处的频率,可由Hilbert变换求出,则对于Ricker子波的主频为:
741 第6期 屈念念等:Riker子波匹配追踪算法及其改进fj=P2*favg(7)为了有效的求得第j次Ricker子波的振幅Aj和相位Hj,对信号u(t)和Ricker子波w做Hilbert变换可得复地震道U(t)和复子波W(t,fj):U(t)=u(t)+iuH(t)(8)W(t)=w(t)+iwH(t)(9)则(6)式可表示为:
U(t)=EjAj*W(t-tj,fj)+Noise,
(10)其中Aj=ajeiUj,(6)式中的幅度Aj表示为的Aj
模,相位Hj表示为的Aj相位。
为了使信号u(t)与匹配的Ricker子波集之间的相对误差:
R(t)=3U(t)-EJ1[Aj*W(t-tj,fj)]42(11)为最小,利用最小二乘法来求的Aj,则A=[WTW+EI]-1WTU(12)因为在实际资料中地震信号初始相位均归零,所以在分解的过程中可以不用求取相位这个参数,默认相位Hj为零。因此在利用匹配追踪法分解时求取振幅Aj,主频fj和中心时间tj。为了有效的提高时频分辨率,避免对过于庞大的原子库进行搜索,可以在由上述步骤中所求得的tj和fj加上自定义的范围,如对于所得的第j次Ric-ker子波的中心时间tj,可以在以tj为中心的20ms范围内搜索,fj亦如此,从而可以得到与原信号最大相关的Ricker子波,就可以利用这个最佳的原子进行下一步最小二乘法反演,这样就能更好地表示原始信号。
3 算法改进从理论角度来讲,求取的Aj要使搜索出来的一系列子波与原始信号尽可能接近,这样能够减少迭代次数,加快收敛速度。从理论上讲,通过最小二乘法就能够得到J个Ricker子波。对于不存在接近于零的特征值和无穷大特征值的行列式WTW,最小二乘法足以利用适当的迭代次数解出最佳的子波集合。但是在实际的地震信号中,行列式WTW的特征值有一部分会趋近于零,利用最小二乘法所计算出的振幅误差较大,然后减去
子波后的信号就要在原来的已经计算的包络处再次计算一个子波,这样会增加迭代次数,降低效率和重构精度。因此,采用新的方法来计算振幅。笔者尝试了奇异值分解[19]的方法来计算振幅,当矩阵W为奇异矩阵时,可进行奇异值分解:W=SrCrVTr(13)其中:r为矩阵的秩,Cr是由WTW之r个非零特征值之正根组成的对角线矩阵,Sr,Vr均为WTW的特征向量矩阵。则A=VrC-1rSTrU(14)此种方法不受行列式特征值的影响,可以有效地减少计算振幅的误差,而且可以减少迭代次数,提高运算速度。另外,在计算过程中,除了设定剩余的地震道幅度为阈值外,还设定了另外一个迭代终止准则,即收敛准则,当迭代N次的剩余信号幅值与迭代N+1次的剩余信号相对幅值差值接近于零时,迭代结果已经不再收敛,停止迭代。如果不设置这个收敛准则,则因为计算误差所剩余的相同延时和频率的子波将会不断分解,只会延长计算时间,而对于重构信号毫无用处。