2012, Vol. 55
在地震资料处理与解释领域，信号的稀疏表示[1-2]近年来越来越引起人们的关注.理想的信号分解方式应该根据信号的特点，自适应选择基函数来分解信号[3]，这对于具有时域和频域局部化特征的地震信号来说尤其重要.Mallat等人提出的匹配追踪算法(MP)正是具有此优点的信号分解方法，该算法通过创建超完备时频原子库，将信号在时频原子库中进行展开，从而实现信号的自适应分解[4].
尽管匹配追踪能够给出较好的信号表示，但该方法为了寻找最佳匹配原子，在信号分解过程中需要大量内积运算，导致计算量过于庞大.此外，为了保持匹配追踪精度，时频原子库往往规模很大，这更增加了匹配追踪的计算量.因此，在保持匹配追踪算法精度的同时，降低其计算量是亟待解决的问题之一.
原匹配追踪方法中，由于时频原子库是超完备的，匹配追踪在某一次迭代中得到时频原子不会与之前迭代得到的时频原子正交，所以在新一轮迭代中又引入了新的时频原子分量，这是由匹配追踪算法的“贪婪"策略决定的，称之为贪婪算法[5].同时，在算法执行过程中，很可能因为某次迭代的“短视"，选择了一个不合适的时频原子，而算法本身并没有提供在后续迭代中重新选择该时频原子的机制，只能通过选择更多的时频原子进行迭代来纠正错误，从而产生“过匹配"现象[6].如何解决匹配追踪算法的过匹配现象是第二个需要解决的问题.
目前随着大规模高密度三维地震勘探的普及，地震数据越来越庞大，因而对匹配追踪来讲，在保持其算法精度的前提下提高计算效率，降低计算量，对于该算法在地震信号处理和解释中的推广是非常重要的[7].本文在分析匹配追踪算法的基础上，利用地震信号的时频局部性特征作为先验信息，采用动态最优搜索策略寻找时频原子，大大降低了匹配追踪分解的计算量.同时，通过时频原子的正交变换，使每次迭代中不再包含已选的时频原子，消除了时频原子库中的冗余分量.
2 匹配追踪算法剖析 2.1 时频原子库的构建时频原子库D是一系列时频原子的集合，常用的时频原子是通过窗函数g(t)的伸缩、平移等操作得到的.窗函数g(t)∈L2(R)满足以下条件：
(1) g(t)是实的连续可微函数；
(2) ‖g(t)‖ =1；
(3) ∫g(t)dt≠0，且g(0)≠0[4, 8].令γ= (σ，u，ξ)，
则时频原子定义为
|
(1) |
其中σ(σ>0)为尺度，1/$\sqrt \sigma $ 用来对时频原子进行标准化，使‖gγ(t)‖ =1；u表示中心时间；ξ 为频率因子.从(1)式看出，时频原子主要能量集中在u附近，大小与σ 成比例.(1)式的频域表示为
|
(2) |
可见|ĝγ(f)| 是偶对称的，它以频率f=ξ 为对称中心，能量主要集中在ξ 附近，大小与1/σ 成正比.σ将窗函数的能量集中在点(u，ξ)附近，从而获得紧凑表示.
为了适应地震信号不同地震子波波形稀疏分布[9]的特征，以高斯函数gσ(t)∈ L2(R)作为母函数，以φ 表示相位，γ = (φ，u，ξ)，则有：
|
(3) |
(3) 式即为地震信号匹配追踪分解中常用的相位时频原子库.
2.2 匹配追踪分解机理设信号f可以分解为gγ0方向的分量及垂直于
|
(4) |
其中〈，〉表示内积运算，R(f)是信号f在gγ0方向投影后的残差信号，如图 1所示.
|
图 1 匹配追踪分解机理 Fig. 1 Principle of matching pursuit decompositon |
由于gγ0和R(f)是两个正交方向，因此|f|2 =|〈f，gγ0〉|2 + |R(f)|2.为了使残差能量|R(f)| 最小，须选择gγ0 ∈ D使|〈f，gγ0〉|最大.
匹配追踪采用迭代算法，开始时设R0(f)=f，迭代到第n次时，残差信号为Rn(f).选择gγn∈ D使gγn与Rn(f)最匹配，即Rn(f)被分解为
|
(5) |
重复此分解过程，直到残差能量小于容许值，信号最终被分解为
|
(6) |
(6) 式表示信号f最终被分解成一系列时频原子的线性叠加.尽管时频原子{gγk}0≤k≤n未必两两正交，但能量保持守恒[10-11].
3 动态搜索快速匹配追踪 3.1 基本原理在匹配追踪每次迭代过程中，为了寻找最佳时频原子，贪婪算法不惜将原子库中所有的原子搜索一遍，这种搜索方式为全局最优搜索(GOS)，它使信号的分解拥有最高的自由度，使“短视搜索"现象出现的可能性降到最低，但该方式是以庞大的计算量为代价的.
由于时频原子具有能量集中性的特点，根据这个特点，在时频原子搜索过程中，可以将参数u限制在以u0 为中心的范围 [u0-nΔu，u0+nΔu](n为正整数，Δu为u的采样间隔)内进行搜索.在该范围内，匹配追踪分解仍然是完全自由的，这样会大大减小计算量，这种搜索策略称为局部最优搜索(LOS)方式.
在匹配追踪分解之前事先得到信号的频率、相位等先验信息，在迭代分解过程中，以u0 处的先验信息来约束时频原子库中相应频率、相位等参数的搜索范围，会在更大程度上减少匹配追踪的计算量.由于先验信息点的位置随着每次迭代的变化而变化，即搜索范围是动态变化的，将这种搜索方式称为动态最优搜索(DOS).
在地震信号匹配追踪中，利用复地震道分析[12]技术，将瞬时属性引入到信号分析中，其中瞬时频率和瞬时相位给出了信号在任意时间点处比较准确的频率和相位信息，因而在匹配追踪分解时，可以引入地震信号的瞬时频率和瞬时相位作为动态最优搜索的先验信息[13-15]，实现快速高效分解.
3.2 模型测试为检验动态搜索匹配追踪的分解效果和计算效率，利用4个高斯函数合成了图 2 左边所示的合成记录，表 1是组成该合成记录的每个高斯函数的参数.分别采用全局最优搜索、局部最优搜索和动态最优搜索三种搜索方式对图 2左边的合成记录进行匹配追踪分解.
|
图 2 合成记录及构成合成记录的高斯函数 Fig. 2 Synthe tictraceand Gaussion functions |
|
表 1 高斯函数的参数列表 Table 1 List of Gaussian function parameters |
匹配追踪分解结果如图 3所示，不论何种搜索方式，迭代都只进行了4次，且每次迭代得到的时频原子都与图 2中的一个高斯函数对应.三种搜索方式得到的第4个高斯函数的参数如表 2 所示.对比表 1中第4个高斯函数的实际参数可以看出，三种搜索方式都能够保持匹配追踪的分解精度.
|
图 3 不同搜索方法匹配追踪得到的时频原子 (a)全局最优搜索；(b)局部最优搜索；(c)动态最优搜索. Fig. 3 Time-frequency atoms by matching pursuit decomposition of different searching strategies (a) Global optimal searching (GOS) ； (b) Local optimal searching (LOS) ； (c) Dynamic optimal searching |
|
|
表 2 三种搜索方式结果对比(以第4个高斯函数为例) Table 2 Comparison of the three searching strategies(the 4th Gaussian function) |
图 4给出了三种搜索方式的计算效率，图中y轴为搜索方式，x轴为各种搜索方式下的不同迭代次数，z轴为计算时间.由图 4a可见不同搜索方式的计算效率差别非常显著，局部最优搜索和动态最优搜索的计算时间远远小于全局最优搜索.将局部最优搜索和动态最优搜索的计算时间放大，如图 4b所示，可以看出动态最优搜索的计算效率是局部搜索方式的14倍，是全局搜索的1000倍.
|
图 4 不同搜索方式计算效率对比 (a)全局最优搜索、局部最优搜索和动态最优搜索的计算效率对比；(b)局部最优搜索和动态最优搜索的计算效率对比. Fig. 4 Comparison of computational efficiency between different searching strategies (a) Comparison of computational efficiency between GOS, LOS and DOS; (b) Comparison of computational efficiency between LOS and DOS. |
由于时频原子库是超完备的，匹配追踪在第i次迭代中得到时频原子gγi并不与之前的时频原子{gγk}0≤k<i正交，所以又引入了gγi分量，这是由匹配追踪算法的“贪婪"策略决定的.为了解决这一问题，运用Schmidt方法对时频原子{gγk}0≤k<i进行正交变换，以消除时频原子库中的冗余分量，方法如下：设R0(f)=f，gγ0所在的空间为V0，PV0(f)为f在V0 上的投影，则迭代一次后的残差信号为R1(f)=f-PV0(f)=f- 〈f，gγ0〉gγ0 . (7)定义Schmidt正交基的第一个分量为u0 =gγ0 .若已有n个线性独立的时频原子{gγp}0≤p<n并有相应的正交基(up)0≤p<n，于是：
|
(8) |
选择gγn∈D使gγn与Rn(f)最匹配，减去信号在空间Vn-1 上的投影，得到正交基的第n个向量un：
|
(9) |
(up)0≤p<n+1 构成了空间Vn的正交基，于是信号f的正交时频原子匹配追踪分解方程表示为
|
(10) |
由(10)式可见，正交时频原子匹配追踪算法将信号f分解成一系列正交时频原子的线性叠加.通过时频原子的正交变换，使每次迭代过程中不仅减去了信号在gγn上的投影，而且减去了信号在所有已选时频原子上的投影，残差信号中不再包含已选时频原子分量.随着迭代分解次数的增加，式(10)右端时频原子的线性叠加可以任意逼近原始信号.
4.2 算法测试图 5左边是一道实际地震信号.以Morlet子波[7, 16]为基本原子创建时频原子库，并采用动态最优搜索策略，首先对该信号进行非正交匹配追踪分解，经过231次迭代，信号被完全分解.将各时频原子进行两两内积运算，其内积值分布如图 6所示，可见非零内积值广泛分布，说明时频原子之间并不两两正交，正是这些不正交的时频原子导致了非正交匹配追踪算法出现“过匹配"现象.
|
图 5 地震信号及重构信号的对比 Fig. 5 Comparison between seismic signal and reconstructed signal |
|
图 6 非正交匹配追踪各时频原子的内积值分布 Fig. 6 Inner product distribution of time-frequency atoms by non-orthogonal matching pursuit |
将图 5左边的地震信号采用动态搜索的正交时频原子匹配追踪算法进行分解，迭代次数为162次.图 5中间的信号是正交时频原子重构的地震信号，右边是实际信号与重构信号的残差，可见重构信号与实际地震信号非常相似，残差极小，说明正交匹配追踪能够实现地震信号的分解.为了对比，同样计算了所有时频原子的内积值，其分布如图 7所示，可见所有内积值均为零，说明时频原子两两正交，不存在冗余分量.
|
图 7 正交匹配追踪时频原子的内积值分布 Fig. 7 Inner product distribution of time-frequency atoms by orthogonal matching pursuit |
图 8给出了两种算法残差能量随迭代次数的变化，其中虚线是非正交匹配追踪算法的，实线是正交匹配追踪(OMP)算法的.从图中可以看出，在开始的几次迭代中，两种算法的收敛速度几乎相同，大约10次迭代之后，正交匹配追踪算法的收敛速度开始占优势.
|
图 8 残差能量随迭代次数的变化 Fig. 8 Residual energy variation with iterative number |
在计算用时方面，由于非正交算法的独立性，每次迭代所消耗的时间是相同的.而在正交匹配追踪算法中，对已选时频原子进行正交变换需要额外的计算时间，从(9)式可以看出，随着迭代次数的增加，这种额外的消耗也随之增加.图 9a中红点和蓝点分别表示每次迭代非正交匹配追踪和正交匹配追踪的计算用时，实线为拟合出的用时趋势，其中红线为非正交算法的用时趋势，蓝线为正交算法的用时趋势，拟合线的走势很好地说明了两种算法的用时特点.
|
图 9 不同匹配追踪方法的计算用时对比 (a)每次迭代用时对比；(b)累计用时对比. Fig. 9 Comparison of computing tme between different methods (a) Comparison of computing time in each iteration； (b) Comparison of total computing time. |
图 9b是累计用时对比，可见虽然每次迭代正交匹配追踪用时较多，但由于正交算法减少了迭代次数，因而总体用时较少.
5 结 论匹配追踪算法能够根据地震信号的时频局部化特征实现自适应分解.尽管匹配追踪有很多优点，但原贪婪算法计算效率极低，满足不了巨量地震数据实时处理的要求，大大限制了该方法在实际数据处理中的应用.
通过对匹配追踪算法的分析发现，计算效率低下的主要原因有两个：一是时频原子库规模庞大，时频原子超完备，存在大量冗余分量；二是在每次迭代过程中，为了寻找最佳时频原子，算法不惜将时频原子库中的所有原子都搜索一遍.针对时频原子超完备问题，本文采用正交变换方法，使信号在每次迭代时不再包含已选的时频原子分量，最终分解为一系列正交时频原子的线性叠加.为了避免搜索过程的“贪婪性"，充分利用地震信号的瞬时信息，采用动态最优搜索策略寻找合适的时频原子，大大降低了匹配追踪分解的计算量.数据分析结果表明，通过动态最优搜索和时频原子的正交变换，不仅保持了匹配追踪的分解精度，而且大幅度提高了算法的计算效率，从而快速实现地震数据的分解.
| [1] | Friedman J H, Stuezle W. Projection pursuit regression. Journal of the American Statistical Association , 1981, 76(376): 817-823. DOI:10.1080/01621459.1981.10477729 |
| [2] | LaMotte L R, Hocking R R. Computational efficiency in the selection of regression variables. Technometrics , 1970, 12(1): 83-93. |
| [3] | Cohen L. Time-Frequency Analysis. Prentice Hall , 1995. |
| [4] | Mallat S G, Zhang Z F. Matching pursuits with time-frequency dictionaries. IEEE Transactions on Signal Processing , 1993, 41(12): 3397-3415. DOI:10.1109/78.258082 |
| [5] | Davis G M, Mallat S G, Zhang Z F. Adaptive time-frequency decompositions. Optical Engineering , 1994, 33(7): 2183-2191. DOI:10.1117/12.173207 |
| [6] | Chen S S, Donoho D L, Saunders M A. Atomic decomposition by basis pursuit. SIAM Journal of Sci. Comput. , 1999, 20(1): 33-61. |
| [7] | Liu J L, Marfurt K J. Matching pursuit decomposition using Morlet wavelets. //Expanded Abstracts of 75th SEG Annual International Meeting. 2005:786-789. |
| [8] | Chakraborty A, Okaya D. Frequency-time decomposition of seismic data using wavelet-based methods. Geophysics , 1995, 60(6): 1906-1916. DOI:10.1190/1.1443922 |
| [9] | Castagna J P, Sun S, Siegfried R W. Instantaneous spectral analysis: detection of low-frequency shadows associated with hydrocarbons. The Leading Edge , 2003, 22(2): 120-127. DOI:10.1190/1.1559038 |
| [10] | 陈发宇, 杨长春. 基于MP方法的地震信号快速分解算法. 地球物理学进展 , 2007, 22(6): 1692–1697. Chen F Y, Yang C C. Seismic signal's decomposition based on matching pursuits method. Progress in Geophysics (in Chinese) , 2007, 22(6): 1692-1697. |
| [11] | 陈发宇, 杨长春. 对依据频率的匹配追踪快速算法的改进. 石油物探 , 2009, 48(1): 80–83. Chen F Y, Yang C C. An improvement for frequency-dominated fast matching pursuits algorithm. Geophysical Prospecting for Petroleum (in Chinese) , 2009, 48(1): 80-83. |
| [12] | Taner M T, Koehler F, Sheriff R E. Complex seismic trace analysis. Geophysics , 1979, 44(6): 1041-1063. DOI:10.1190/1.1440994 |
| [13] | Wang Y H. Seismic time-frequency spectral decomposition by matching pursuit. Geophysics , 2007, 72(1): 13-20. |
| [14] | 陈林, 宋海斌. 基于Morlet小波匹配追踪算法的地震时频属性提取. 石油地球物理勘探 , 2008, 43(6): 673–679. Chen L, Song H B. Extraction of seismic time-frequency attribute based on Morlet wavelet match tracing algorithm. Oil Geophysical Prospecting (in Chinese) , 2008, 43(6): 673-679. |
| [15] | Liu J, Marfurt K J. Instantaneous spectral attributes to detect channels. Geophysics , 2007, 72(2): 23-31. |
| [16] | Morlet G, Arens G, Fourgeau E, et al. Wave propagation and sampling theory-Part II: sampling theory and complex waves. Geophysics , 1982, 47(2): 222-236. DOI:10.1190/1.1441329 |