地震勘探是一种重要的物探方法(付燕，2002).通过人工激发地震波，分析地震子波在地层中的传播规律，从而获得不同地层间的反射强度，了解地层结构，达到勘探的目的.由于真实的地质结构非常复杂，因此通过地层剖面来分析地质结构是非常困难的.稀疏分解方法一经提出，许多研究学者就将其应用到地震数据解释中.

在传统的人工地震记录(徐伯勋等，2001)中，是把大地地层看成一个线性系统，即假设激发震源产生一个尖脉冲，经过地层传播后仍然保持脉冲的波形不变.然而，在实际的地震记录中情况是非常复杂的.激发震源所产生的尖脉冲经地层介质的吸收后，尖脉冲的波形将会被拉长，这种具有时延的单个反射波被称为地震子波.从而实际的地震记录 x(t)由反射系数γ(t)和地震子波w(t)褶积组成(陈玉东，2006)，公式为

由于震源信号经过地层传播时，信号在地层间的反射是有所差别的，从而地层反射系数γ(t)可以抽象为一有限长序列γ(n).不同地层之间对震源信号的反射系数可以简化为γ(n)上的非零数值，而同一地层内相对较小的反射系数则简化为0(Charles等，2005).因此地震记录的离散化公式为

信号稀疏分解的思想最早是由Mallat和Zhang(Mallat et al.，1993)在1993年提出来的，信号稀疏分解是信号在过完备原子库上根据信号特点选取最佳匹配原子，用选取的最佳原子的线性叠加得到一个非常简洁的稀疏表达式的过程.稀疏分解以其灵活、简洁、稀疏这一优良特性已经被广泛应用到地震资料的解释中，如地震信号的去噪(赵天资等，2008)、谱分解(张显文等，2010；Wang et al.，2007)、反演(Ismael et al.，2012)、存储预测(武国宁等，2012)、快速分解(陈发宇等，2007； 张繁昌等，2012)、滤波(李学良等，2013；石颖等，2013)等方面.

目前，在地震信号处理中常用的稀疏分解算法是匹配追踪(Matching Pursuit，MP)算法(王建英等，2006).Chakraboty等(1995)最早将MP算法用于地震信号时频属性分析；Liu等(2004)提出基于Ricker小波匹配算法的地震信号分解；王成梅等(王成梅等，2010)提出了改进的Ricker子波的地震信号MP分解，降低信号分解后的均方误差.以上方法应用MP对地震信号进行处理，取得了比较好的效果，但MP选取的原子间存在相关性，因而信号的分解不够稀疏.为了克服上述问题，本文采用正交匹配追踪(Orthogonal Matching Pursuit，OMP)方法(杨福生，2006)对信号进行稀疏分解，要求在每步选取的最佳原子必须与已选取的原子正交，来提高信号分解的稀疏度.与MP比较， OMP中添加了正交化过程，但同时增加了算法的计算复杂度.为了降低计算复杂度，采用一种全局搜索算法—重复加权提升搜索(Repeated Weighted Boosting Search，RWBS)算法，来快速搜索最佳匹配原子.

目前稀疏分解的算法很多，大多是先建立一个过完备原子库，再通过一定的准则从过完备原子库中选取与观测信号最佳匹配的原子，最后利用选取的最佳原子进行线性叠加得到与原始信号最佳逼近的稀疏表达式.

雷克子波(Ricker Wavelet)是人工合成地震记录时常用的零相位子波，并在地震信号滤波处理中得到广泛应用(邓小英等，2011).Ricker子波在时间域上的表达式为

图 1

Fig. 1

图 1 一个Ricker原子 Fig. 1 A Ricker atom

设 G={g _ri} r _i∈Γ为建立的过完备原字库， g _ri为过完备原字库中的原子， r _i这个指标与构建原子的参数及参数个数有关，Γ为由r_i这个指标组成的集合.一般将情况下需对原字库中的所有原子进行归一化处理，即=1.假设研究的信号是 y ，信号的长度为n，构建过完备原子库的原子为零相位Ricker子波.

1.1 基于MP的信号稀疏分解

首先，在过完备原字库 G 中寻找与待分解信号y最佳匹配的原子，记作 g _n1，其满足条件

然后，从 y中减去y沿着g _n1方向的投影，就得到残余信号 R₁y ，同时在原字库中剔除选取的最佳原子，即有

按照上面的方法依次选择与残余信号最配的原子，经过k次匹配后，信号被分解为

其中 R₀y=y ，R_ky为原始信号分解后的残余，剩余原字库为 G _k= G _k-1－{ g _nk}.当残余能量小于或等于某个给定阈值ξ₁(在对加性噪声信号进行分解时，ξ₁的值一般要比加性噪声方差稍大一些；在无噪声情况下，ξ₁的取值要尽量小)时停止迭代.选取的原子数目愈少时，信号就分解的越稀疏.

由于在MP算法过程中，只能保证 R _k y与g _nk正交，但不能确保 R _k y与g _n1, g _n2,…, g _nk-1正交，因此所得到的结果是次最优的，即信号分解不是最稀疏的.研究MP的改进算法来减少分解所需要的原子个数已成为信号稀疏分解的热点.

1.2 基于OMP的信号稀疏分解

与MP算法相比较，OMP算法在MP算法每选取一个最佳匹配原子之后，即寻找下一个最佳匹配原子之前，将原字库中剩余原子与已选取原子所张成的子空间的一组基做Gram-Schmidt正交归一化处理.第一次提取的原子 g _n1就是第一个正交基 e ₁，它所张成的空间记作V₁.则y在e₁方向上的投影为 P _V1y=〈 y,e ₁〉 e ₁，剩余原字库为 G₁=G－{g _n1}.

设在第k次选取最佳匹配子时，需要将剩余原字库 G _k-1中所有的原子与已得到的一组正交基 e₁,e₂,…,e _k-1 做Gram-Schmidt正交归一化处理，即有

在正交归一化处里后的原子中选取e_k的目标函数为

从而由 e₁,e₂,…,e _k 张成的子空间记为V_k，则 y减去y 在V_k上的投影，就得到新的残余和剩余原字库，即为

按照上面的(7-11)式进行迭代，直到小于或等于某个给定的阈值ξ₁(ξ₁的取法同上)时，就停止迭代.

由上面MP和OMP算法的介绍可知，MP和OMP算法的区别在于OMP算法选取的所有原子都是两两相互正交的，而MP算法没有正交过程.正因为原子得到了正交化处理，所以若达到相同的重构精度，OMP算法的分解的步数应该比MP算法的少，但同时也增加算法的计算复杂度.此外，MP和OMP的分解效果都与建立的过完被原子库的大小有关，当原子库越大时，分解的结果越稀疏但速度较慢；当原子库越小时，分解的速度提高了但不够稀疏.为了克服这些问题，本文提出了基于RWBS的正交匹配追踪算法.

由于基于匹配追踪算法的稀疏分解的每一步，实际上是解一个最优化问题，所以可以将寻找最佳匹配原子的过程转化为利用智能算法求函数最优化的问题.选取最佳匹配原子的搜索方法有很多，如人工鱼群算法(舒伟杰等，2009)、模拟退火算法(徐帅等，2008)、遗传算法(张静等，2008)等等.本文将应用RWBS算法，它是一种新的全局搜索算法，易于实现.在使用该算法时，需要设定一些参数.为了方便下面算法的描述，设定了如下参数： Ps为种群初始个数，NG为外部循环迭代次数，NB为内部循环迭代最大次数，ξ为内部循环终止条件之一.

RWBS的内部循环是一个加权提升搜索过程.它先结合信号的先验信息，随机的产生Ps个种群：{ u ′ _i|i=1,2,…,Ps}，根据(3，7-9)式计算J(u).再寻找当前最优和最差种群，将其分别记作 u _best=argminJ(u), u _worst=argmaxJ( u ).对种群{ u ′ _i|i=1,2,…,Ps}用凸组合的方式产生一个新的种群 u ′ _Ps+1，同时利用镜像原理产生一个 u ′ _Ps+1关于 u _best的对称点 u ′ _Ps+2.最后用 u ′ _Ps+1和 u ′ _Ps+2中使得J( u)最小种群来替换u _worst.该过程重复进行直到达到最大迭代步数NB，或 u ′ _Ps+2－ u ′ _{Ps+1 2}≤ξ时停止迭代.

图 2是一个简单的说明加权提升搜索算法的二维例子，在这个例子中Ps=3， u₁,u₂,u₃是三个初始种群，且由图可知u _best= u₁,u _worst= u₂.u₄是u₁,u₂,u₃加权组合得到的种群，u₅是u ₄的镜像，且J( u ₄)>J( u₅)，因此用u₅来代替u₂.在下一步循环中的三个种群就是u₁,u₃,u ₅，按照上面的过程循环，直到达到终止条件.在这个过程中加权组合的方式很关键，关于它的具体说明可以参考文献Ron等(2003).

图 2

Fig. 2

图 2 简单加权搜索最优化过程 Fig. 2 The optimization process of a simple weighted search

加权提升搜索属于局部最优化过程，对初始选择的种群依赖性很强.为了克服局部最优的特性，将加权提升搜索过程重复多次来得到一个全局最优解.这种得到全局最优解的过程就是RWBS算法，关于RWBS算法的更多细节可以参考文献Sheng等(2005).当用RWBS搜索第k个匹配原子时，算法流程如图 3.

图 3

Fig. 3

图 3 用RWBS寻找第k个匹配原子的算法流程图 Fig. 3 The algorithm flow chart of looking for the kth matching pursuit atom by RWBS

RWBS采用加权提升搜索最优过程来寻找一个较优的种群 u ^k_t,t=0,1,…,NG－1，然后将其融入下一步循环中.保留全局最优值 u ^k_NG作为新选择的匹配原子参数组，删除其余的种群.

本文提出的基于RWBS的快速分解算法步骤如下：

步骤1 给定参数Ps,NG，NB以及ξ的值，以及最终的迭代终止阈值ξ₁的值，置初始残余 R _o y=y,k =1；

步骤2 利用RWBS 寻找与 R _k-1 y最佳匹配的原子，记作u _k，再将 R _k-1 y沿着u _k的方向投影，就得到残余信号 R _ky；

步骤3 当 R _k y 小于阈值ξ₁时，停止迭代；否则令k=k+1，返回到步骤2.

针对现有的基于MP算法的地震信号稀疏分解方法稀疏度不够高的局限性，提出了结合OMP和RWBS的改进算法.该算法不仅有效地提高信号分解的稀疏度，而且降低了OMP算法的运行速度.无论是对人工合成的地震信号还是对实际的地震信号，改进算法对地震信号的分解都比MP更加稀疏，而且对受噪声污染的信号具有较好的鲁棒性.

匹配追踪方法的思想是选择与待分解信号最佳匹配的原子，进行信号的稀疏分解.但匹配追踪的难点在于如何选取核模型，如果与构成原始信号的子波波型相近，则不仅会反应出信号的一些有用的特征，而且还会提高分解的稀疏度.因此，如何构建更好的核模型来进行信号的稀疏分解是接下来研究的重点.