在地震勘探资料中，有效信号往往被湮没在强随机噪声中，阻碍有效信息的恢复.因此，压制低SNR地震资料中随机噪声并最大限度保留有效信号，是地震勘探学界研究的难点和热点.

时频峰值滤波(TFPF)是消减地震勘探随机噪声的有效方法之一(金雷等，2005;Lin et al.，2007;林红波等，2008)，具有约束条件少，对低SNR地震资料适应性强等优势.为了有效压制噪声且恢复有效信号，则要求信号是线性的(Boashash and Mesbah，2004).因此，在应用TFPF处理实际地震资料时，采用加窗的Wigner-Ville分布(WVD)，即PWVD来实现地震信号的局部线性化(李月等，2009).现有研究表明，有效的压制随机噪声和信号幅度保持对窗长的要求是矛盾的，为了得到理想的去噪效果需要大窗长，但在大窗长内信号是非线性的，从而导致信号幅度衰减(Tian and Li，2014).此外，由于地震勘探资料具有非平稳的特征，单一窗长无法同时满足复杂地震信号的局部线性化条件，从而导致该方法处理含噪信号时存在有效信号损失及噪声压制不足的问题(Lin et al.，2014;Zhuang et al.，2015).为此，Zeng等(2015)提出了基于凸集和Viterbi算法的自适应TFPF(ATFPF).该方法将不同窗长的TFPF滤波结果组成一个凸集，然后在这个凸集上通过Viterbi算法寻找最小化给定的全局二次泛函所在位置的滤波结果作为有效信号的最优估计.该方法在去噪的基础上，能够快速地跟踪信号变化，无畸变地恢复出有效信号.但此方法也存在一定的问题：首先，该方法是时域一维滤波，并没有考虑到信号的空间信息；其次，该方法使用的全局二次泛函是能量函数，所以当处理较低SNR的含噪信号时，不能很好地估计出有效信号.

为了解决上述问题，本文提出了空间局部加权回归自适应TFPF(SLWR-ATFPF)的二维滤波算法.首先，对含噪信号进行空间局部加权回归(SLWR)，使处理之后的信号包含空间信息(Deng et al.，2010，2011).然后，再对包含空间信息的含噪信号进行基于凸集和Viterbi算法的自适应TFPF，从而完成时空二维自适应滤波.

本文所考虑的地震勘探资料可以用如下方程建模：

其中，x∈表示有效信号，g∈表示随机噪声.利用空间加权回归对含噪信号s进行处理，得到空间加权信号S可以表示为：

式中，C_rk为道间回归系数，C_tk为时间回归系数.

由于地震勘探信号具有非平稳的性质，我们引入局部的思想.对含噪信号s进行SLWR(玄海燕等，2013).选用的局部空间的大小为m×n(m，n均为奇数).则局部道间回归系数c_rk可以通过局部加权最小二乘得到.具体的过程如下：

设对应于采样点t_(n-1)/2上r₁～r_m地震道的观测值为s_(n-1)/2，j=[s_(n-1)/2，1，…，s_(n-1)/2，m]，进行加权最小二乘拟合，即使得

-…-β_q(r_(n-1)/2，j)r^q_(n-1)/2，j]²最小.

令R为Vandermonde矩阵，β_r为拟合系数矩阵，W_r为加权矩阵，且需要满足文献(Cleveland，1979)中的4个要求.利用加权最小二乘理论，得到的道间加权回归值为(虞乐和肖基毅，2012)：

则在r_(n-1)/2，j处的拟合值为：

其中，(t_(n-1)/2，j)为拟合系数的估计值，c_rk(r_(n-1)/2，j)为局部道间回归系数.由于采用的是滑动窗局部加权回归，因此只需要求得局部道间回归系数c_rk(r_{(n-1)/2，(m-1)/2})，即为c_rk，求解完毕.

上述求解c_rk的过程可以通过相应的变换进行简化.如下所示：

由于加权矩阵W_r为对角阵，因此可以对其进行Cholesky分解，得到：

把(5)式代入(3)式中，可以得到：

对公式(6)的证明见附录A.令B=S_rR，P_r=B(B^TB)^-1B^T，则得到S_(n-1)/2，j=P_rs_(n-1)/2，j.不难发现，P_r为正交投影矩阵(庞善起，2005).对矩阵B进行正交三角分解，得到：

则

综上可知，局部道间回归系数：c_rk=P_r(:，(m-1)/2).

从而，完成对局部道间回归系数求解的简化.同理可以得到沿采样点方向的局部时间回归系数c_tk=P_t((n-1)/2，:)，其中P_t为时间方向上的正交投影矩阵.则

由于随机噪声在各个位置的方向随机性，以及有效信号在各个位置的方向确定性，可以得到：

因此，经过SLWR处理之后得到的空间加权信号中，有效信号部分得到了增强，随机噪声得到了衰减.

TFPF是一种基于时频分析理论的噪声消减算法(Boashash and Mesbah，2004).以(1)式含噪信号的模型为例，简述一下TFPF的实现过程.

首先，对含噪信号s进行调频，使其成为调频信号的瞬时频率(IF)，

其中，μ>0为调频系数.

然后，计算关于z的PWVD：

其中，w∈为窗函数，且其长度为2L+1.

最后，从(W)_i，f中得到峰值处的瞬时频率作为x的估计.

简化上述三个步骤，TFPF最终可以等效成一个低通滤波(Zeng et al.，2015)，表示成如下形式：

其中，*表示卷积运算，h∈表示TFPF的冲激响应，且和窗函数w具有相同的长度.

由上面可知TFPF等效成一个低通滤波器h，则对于不同窗长L_v的TFPF可以得到一个冲激响应集合H={h_Lv:L_v≥1，v=1，2，…}，则H的凸包为：

显然，conv(H)为凸集，h_c为不同窗长TFPF冲激响应的凸组合(Zeng et al.，2015).根据卷积运算的线性特点，空间加权信号S的估计集合Y_c为：

也为凸集.且Y_c=conv(Y)， Y={y_Lv=h_Lv*S，h_Lv∈H}.从Y_c中选择y等效于从conv(H)中选择h_c.再根据Hilbert投影定理，可知某定义在的距离泛函J(y)一定可以在Y_c上取得全局最小值(Boyd and Vandenberghe，2004).

令定义在Y_c上的距离泛函J(y)为：

其中，A，B∈^n×n为对角阵，l，u∈Y_c为Y_c的边界，D∈^(n-1)×n为差分矩阵(Lax，2002).

J(y)的表达式中，第一项表示对观测信号的跟踪，第二项表示对观测信号的平滑.A，B控制着跟踪与平滑程度之间的平衡.为了使y能够自适应于信号的变化，一个合理的想法是选取参数为(Zeng et al.，2015)：

其中，A体现了对Y_c内信号最大衰减幅度的跟踪，B体现了对噪声最大衰减的跟踪.h_ui和h_li分别为上下边界处等效成的滤波器，σ为噪声的标准差的估计值.由于该最优化问题往往是病态的，故选择Viterbi算法作为求解工具(Forney，1973).SLWR-ATFPF的具体步骤如下：

(1)对含噪信号s进行SLWR处理，得到空间加权信号S.

(2)给定滤波器序列H，求解空间加权信号S的估计集合Y_c，以及Y_c的上下边界u，l.

(3)在每个时间采样点处取区间[l_i，u_i]的分割点集P_i，使得P_i={y_i，k:y_i，k=l_i+k(u_i-l_i)/N，k=0，1，2，…，N}，其中N为层数，以y_i，k为顶点，可以组成网格图 1.

图 1

Fig. 1

图 1 SLWR-ATFPF的网格图Fig. 1 Trellis diagram of SLWR-ATFPF

(4)根据泛函J(y)定义相邻顶点的距离如下表达式所示：

其中，y_i，a，y_i-1，b分别为P_i和P_i-1中的顶点，且其方向为P_i-1指向P_i.S_i∈S，为对应空间加权信号中，第i个采样点处的值.

(5)根据Viterbi算法，在网格图中寻找距离最短路径y，也就是使得J(y)最小的最优解，则y即为有效信号x的最优估计.

应用SLWR-ATFPF处理单炮反射地震资料模型.采用Ricker子波构建模拟地震记录，包含两个反射轴，主频分别为45 Hz，30 Hz，对应的速度为2100 m·s^-1，2600 m·s^-1，采样频率为500 Hz.加入高斯白噪声使SNR为-5 dB，如图 2a所示.分别采用TFPF(窗长为7)，ATFPF(H={h₁，h₁₃，h₅₃，h₁₀₂₃})以及本文提出的SLWR-ATFPF(H={h₁，h₁₃，h₅₃，h₁₀₂₃})对含噪信号进行处理，得到的结果如图 2所示.对比图 2a—2d，不难发现图 2d轴向更连续.为了使效果更加明显，我们取出第23道进行幅度和频谱的对比，如图 3所示.通过对比不难发现，SLWR-ATFPF在波峰、波谷处幅度的保持，以及频谱成分的跟踪均优于前两种方法.这主要是由于TFPF仅采用单一窗长，而且没有包含道间相关信息；ATFPF虽然在时间方向上可以得到最优解，但是当SNR较低时，基于能量泛函J(y)的最优解判断，就会把有效信号当成噪声去掉；而SLWR-ATFPF克服了上述问题，通过时空二维滤波，先将信号进行空间加权，再处理，从而得到较为理想的结果.表 1列举了不同方法处理之后的SNR，不难发现SLWR-ATFPF处理之后得到的SNR最高.

表 1(Table 1)

表 1 信噪比定量对比(单位：dB) Table 1 Quantitative comparison of SNR (unit: dB) 滤波前 TFPF处理后 ATFPF处理后 SLWR-ATFPF 处理后 0 3.8038 11.2520 12.1613 -3 2.9149 8.0332 9.6366 -5 2.0628 4.9838 6.0511

表 1 信噪比定量对比(单位：dB) Table 1 Quantitative comparison of SNR (unit: dB)

图 2

Fig. 2

图 2 合成地震记录结果比较 (a) 含噪合成地震记录; (b) TFPF处理结果；(c) ATFPF处理结果; (d) SLWR-ATFPF处理结果.Fig. 2 Comparison of results of synthetic records (a) Noisy synthetic seismic records; (b) Result of TFPF; (c) Result of ATFPF; (d) Result of SLWR-ATFPF.

图 3

Fig. 3

图 3 单道合成地震记录结果比较 (a) 23道45 Hz的波形; (b) 23道45 Hz的频谱；(c) 23道30 Hz的波形; (d) 23道30 Hz的频谱.Fig. 3 Comparison of results of single channel synthetic records (a) Waveform of 45 Hz of 23rd channel; (b) Spectrum of 45 Hz of 23rd channel; (c) Waveform of 30 Hz of 23rd channel; (d) Spectrum of 30 Hz of 23rd channel.

将本文提出的SLWR-ATFPF应用于168道，4000采样点，采样频率为1000 Hz共炮点记录(见图 4).其中，图 4a为原始记录，图 4b为窗长为9的TFPF处理结果，图 4c是ATFPF(H={h₁，h₉，h₄₉，h₄₉₉})处理的结果，图 4d是本文提出方法(H={h₁，h₉，h₄₉，h₄₉₉})的处理结果.对比四幅图片，不难看出SLWR-ATFPF处理的结果中对随机噪声有更好的压制效果，为了体现改进算法的保幅效果，我们对方框里的结果进行了放大，对应的放大图如图 5所示.通过对比不难发现，SLWR-ATFPF处理之后，同相轴更加连续.

图 4

Fig. 4

图 4 共炮点记录处理结果对比 (a) 23道45 Hz的波形; (b) 23道45 Hz的频谱；(c) 23道30 Hz的波形; (d) 23道30 Hz的频谱.Fig. 4 Comparison of results of single channel synthetic records (a) Waveform of 45 Hz of 23rd channel; (b) Spectrum of 45 Hz of 23rd channel; (c) Waveform of 30 Hz of 23rd channel; (d) Spectrum of 30 Hz of 23rd channel.

图 5

Fig. 5

图 5 局部共炮点记录处理结果对比 (a) 局部实际地震勘探记录； (b) TFPF处理局部放大结果； (c) ATFPF处理局部放大结果； (d) SLWR-ATFPF处理局部放大结果.Fig. 5 Comparison of local common shot point seismic data processing results (a) Local real seismic record; (b) Local enlarged result of TFPF; (c) Local enlarged result of ATFPF; (d) Local enlarged result of SLWR-ATFPF.

本文针对传统的TFPF仅采用单一窗长且没用利用空间相关信息的问题，研究了SLWR对TFPF的影响，并提出了SLWR-ATFPF.通过分析TFPF的特点，将SLWR处理之后的空间加权信号作为输入，再进行ATFPF，从而完成时空二维滤波.仿真实验和共炮点记录处理结果表明，本文提出的算法降低了滤波后有效信号在波峰和波谷处的衰减，具有更好的保幅性，同时对噪声有良好的压制作用；在实际的地震资料处理中，该方法也得到了较好的处理效果.

已知S_(n-1)/2，j=R[R^TS_r^TS_rR]^-1R^TS_r^TS_rs_(n-1)/2，j.

对上面公式进行整理得到：S_(n-1)/2，j=R[(S_rR)^TS_rR]^-1(S_rR)^TS_rs_(n-1)/2，j.

已知S_r为对角矩阵，利用交换律，得到S_(n-1)/2，j=(S_rR)(S_rR)^T(S_rR)^-1(S_rR)^Ts_(n-1)/2，j.