近年来，小波凭借其时频域上良好的局部化特性，在地震数据随机噪声压制方面应用非常广泛.常用的二维小波是由一维小波张量积构成的可分离小波，只能对水平、垂直、对角三个方向进行滤波，其基无方向性，对各个方向进行相同程度滤波，滤波系数难以稀疏表示同相轴的取向信息，对方向边缘细节具有模糊作用，一维小波精确的时频定位性在二维中难以体现(Zhang and Yin，2004；Yu and Whitcombe，2008).为了克服小波无方向性的不足，引进了Ridgelet变换(Candès，1998)，先进行Radon变换，把面上的直线映射到面上的点，再进行小波变换，一定程度提高小波对线的检测.然而实际中，同相轴曲线居多，不能理想的转化为集中的点，在地震信号中较少使用.后来，引进了块剖分的构造思想，在局部中，曲率足够小，曲线可近似作直线来处理，即Curvelet变换.Curvelet解决了小波方向性不足的缺点，进行了子带划分，具有明确方向的基，能够稀疏的表示图像平滑区域，也能较为稀疏的表示边缘部分，有较好的适用性.但Curvelet的方向性是靠波数增加、方向滤波器个数增多的方式来增强，冗余度增加，数据存储空间庞大，同时分块过程必然引入边缘效应(Candès and Donoho，1999；Candès et al.，2006；林春和王绪本，2009；Oliveira et al.，2012；袁艳华等，2013).Contourlet变换在低冗余度的情况下，对Curvelet进行离散化实现，采用拉普拉斯金字塔实现多尺度、多方向的分析，虽然在对地震资料处理中有很好的效果，但是依旧无法克服Curvelet变换本身的不足(Do and Vetterli，2002；Li et al.，2008).

我们考虑是否有一种滤波器，在方向上连续，能够不进行插值计算得到任意方向的滤波器.Freeman和Adelson(1991)等提出了方向可控滤波器的理论框架，其基本思想是沿任意方向的滤波器都可以表示成一组基滤波器函数的线性组合形式，这样一来好似滤波器可以沿任意方向转动，因此得名.按照方向可控滤波器的思想，输入图像与固定的基滤波器卷积即可高效的估计出一副图像与任一旋转角度滤波器的卷积结果，极大的减少了计算量；在对方向特征(边缘、脊等)分析上，方向可控滤波器有很多应用.结合Hilbert变换构成正交滤波器组，郑林等(2002)提出基于方向可控滤波器的图像融合新方法；施俊等(2009)提出一种封闭几何图形的自动检测算法；张大明等(2005)则将方向可控滤波器与小波分析结合，也达到了很好的图像融合的目的，并证明了该方法效果优于传统的小波融合图像方法.根据定向滤波器的思想，方向可控滤波器经常应用于金字塔算法中，对图像进行多尺度、多方向分析.董立娟和练秋生(2008)根据此思路提出了一种改进的视觉信息处理模型，有效的克服了人脸识别中的一个重要难题；林春等(2011，2012)将方向可控金字塔用于对三维地震切片的分解中，不仅能够分析不同尺度方向上断层的特征及走向，还能获取断层的信息.通过方向可控滤波器获取纹理走向，被应用于掌纹(王艳霞和阮秋琦，2008)、虹膜(李志慧，2007)等模式识别.可见，方向可控滤波器主要集中于图像处理方面，如果将其成功引入到地震信号保幅去噪中，将非常有意义.本文提出了一类新型方向可控滤波器用于地震勘探信号去噪处理，其高方向选择性使得可以通过分析输入图像每个像素的局部模式角度信号来分离不同方向特征以及噪声，从而达到轴的增强，噪声的压制.

Freeman和Adelson(1991)采用的滤波核为厄米特-高斯函数(HGFs：Hermite-Gauss Functions)，由厄米特多项式和无向性的高斯函数之积来表示.HGFs的一个固有局限是由于高斯函数的无方向性，滤波时对各方向都进行了平滑作用，细节保护不好.于是在高斯谐波函数中引入一个参数对HGFs进行拉伸，形成如图 1d的拉伸厄米特高斯函数(E-HGFs：Elongated Hermite-Gauss Functions).E-HGFs具有明确方向性的特征，对轴敏感度增强的同时，对梯度方向的平滑作用也减弱(Perona，1995).图 1给出了HGFs和E-HGFs的简单对比，从图 1e和图 1f的滤波结果来看，E-HGFs滤波去噪的同时能更好定位轴的位置.E-HGFs有个重大的缺陷，即不具备可控性，需要设计出可控的E-HGFs的拟合函数，使其具有方向可控的特性，在2.3部分将给出拟合过程.利用E-HGFs的高方向选择性和高效性，在局部进行滤波处理，根据噪声无方向性，考虑构造和层序的方向，结合统计特性，我们就可以沿着地震数据同相轴方向进行滤波.在第3部给出的理论模型和共炮点资料处理结果表明，本文算法可以很好的达到压制随机噪声、保持幅度的目的.

图 1

Fig. 1

图 1 HGFs和E-HGFs滤波的简单对比 (a)(b) 合成原始图像及其加噪图像; (c) 无方向性的滤波核即HGFs; (d) 拉伸之后的滤波核; (e)为(b)与(c)的卷积结果; (f)为(b)与(d)的卷积结果.Fig. 1 Comparison of filters HGFs and E-HGFs (a) and (b) Synthetic image and noisy image; (c) Isotropic kernel; (d) Elongated convolution kernel; (e) Output of convolution of the image in (b) with (c); (f) Output of convolution of image in (b) with (d).

如果一组滤波器f(r，λ)(r为笛卡尔坐标向量，λ为这组滤波器所依赖的参数)能够展开为一组基滤波器V_n(r)(n=1，2，…N)的线性组合，且线性系数只与λ有关，那么就说该滤波器具有可控性.可用公式(1)表示为

其中a(λ)为λ的函数，V(r)称为基滤波器.

如果我们用这组滤波器卷积处理一副图像I(r)时，可用公式(2)表达为

由(2)式可看出只要I(r)*V_n(r)能够得到，a(λ)确定，那么任意λ值的滤波结果就能被估计出来.通过线性卷积，对λ无需离散化取值，可以认为Y(r，λ)是关于λ的连续函数，具有精确可导可积的性质，克服了传统通过离散化参数λ得到滤波器产生的离散插值误差，同时能减少大量的卷积操作.

，对一个滤波核的旋转操作为f(r，θ)=R_θ(f).根据李群法(Teo and Hel-Or，1998)和奇异值分解(Perona，1995)可以得到，当V_n(r)为圆谐函数CHFs(circular harmonic functions)时，此时具有可控性，我们把r换到极坐标(ρ，Φ)上，V_n(ρ，Φ)=U_n(ρ)e^inΦ，则：

如果一个滤波核T(r)可以展开成一系列适当的CHF之和，公式为

那么对于T(r)的旋转变换可化为

对比(1)式和(5)式，可以得到滤波核T(r)具有方向可控性，我们称之为方向可控滤波器.处理图像结构(Perona，1995)如图 2所示.

图 2

Fig. 2

图 2 对图像任意方向处理结构Fig. 2 Structure of processing an image in any orientation

若二维函数H_k，l(x，y;μ)满足式(6)，H_k，l(x，y;μ)即为拉伸厄米特高斯函数(E-HGFs)，公式为

其中H_p(x)=(-1)^pe^x2(d/dx)^pe^-x2是p阶Hermite多项式，通过参数μ来控制拉伸程度.当μ=1时，高斯项退化为无方向性，即经典的HGFs.HGFs函数简单，用起来较为方便，然而由于其无方向性的局限，滤波时对各个方向均进行平滑作用，降低了检测精度.E-HGFs凭借其高方向性，一方面克服了折中检测效果和定位精度的Canny(1986)局限，另一方面可以通过分析不同方向的响应来区分噪声和感兴趣的特征，应用前景更为广泛.图 3为圆盘的滤波效果(此图引用Papari(2012)文献中Fig. 2)，(c)中响应只在垂直方向的中间位置有较大值，很好地展示了E-HGFs对方向定位的高精确度.

图 3

Fig. 3

图 3 HGFs和E-HGFs一阶导对合成圆盘的滤波响应 (a) 合成圆盘图片; (b)(c) HGFs、E-HGFs一阶导分别与(a)的卷积结果(Papari et al.，2012). Fig.3 Filter response to the first-order derivative of HGFs and E-HGFs for synthetic image (a) Synthetic image; (b) and (c) Results of convolution with first-order derivative of HGFs and E-HGFs with (a)(Papari et al.，2012).Fig. 3

然而，E-HGFs是不可控函数，即在式(5)中N取值无限，因此必须找出可以拟合E-HGFs的可控函数.在方向可控滤波器发展史上，主要有基于李群理论的方法以及基于奇异值分解的方法.李群法结果为频谱连续的函数，无法进行有限因子分解.奇异值分解法得到的是一系列正交可控性函数，解决了李群法的无限问题，但结果经常是无法解析的，难以解出数值解.为克服以上方法的不足，本文采用Giuseppe Papari(2012)等提出的在奇异值分解理论框架上的通过引入特定矩阵算子来构建紧凑可控的小波函数的方法，以下为其推导：

首先应用Papari(2012)定义的移位矩阵算子S，对于一个任意向量r，{r}ⁿ表示向量r中的第n个元素，那么{S^kr}ⁿ={r}^n+k，即移位了k个元素，然后又定义三个方便使用的算子为

现在令向量 为基滤波器H_k，l(x，y;μ)的径向成分，那么根据厄米特多项式与高斯函数的关系可以得到公式为(Papari et al.，2012)

其中，为第一类型修正k阶贝塞尔函数，即.由(7)(8)式可以得到公式为

根据(9)式以及厄米特多项式与高斯函数的关系，可以得到U_k，l(ρ，μ)用贝塞尔函数展开的式子，Papari(2012)等给出了几个详细的结果.本文采用k=0，l=0的滤波器对地震记录进行方向滤波，其表达式为

图 4列出了n={0，1，2，3，4，5}的基滤波器.

图 4

Fig. 4

图 4 k=0,l=0的基滤波器(a,b,c,d,e,f分别对应n={0,1,2,3,4,5})Fig. 4 Basic filters when k=0,l=0 and n={0,1,2,3,4,5}

本文采用E-HGFs对地震数据进行方向滤波，首先根据每一点各个方向的滤波响应方差初步区分噪声点与有效信号点，然后再与局部小块均值进行对比，进一步确认，最后按照方差大小比例对数据进行非线性滤波，这样处理的实质是对信号引入了方向维度，滤波器沿着地震数据同相轴方向进行滤波.算法流程总结如下：

(1)采用Papari等(2012)提出的方法拟合出方向可控的E-HGFs滤波器.

(2)输入共炮点记录，进行方向滤波，构建“时间-空间-方向”三维数据.

(3)根据噪声的无向性，在方向维度进行方差计算，方差小于一定门限时，初步判定为噪声点.

(4)计算邻域局部均值，结合第(3)步的方差，再次判定，防止误判.

(5)对判定的噪声点进行压制，其余点取其最大的方向滤波器结果，最后用带宽足够大的低通滤波器去除毛刺.

模拟一幅200道，采样时间0.002 s的人工合成地震记录.浅层和深层同相轴的主频分别为30 Hz和23 Hz，速度分别为1800 m·s^-1和2300 m·s^-1，在记录中加入-5 db高斯白噪声.图 5a和b分别为原始记录和加噪记录，(c)(d)(e)分别为小波变换、Curvelet变换及本文提出方法的滤波结果，(f)(g)(h)分别为小波变换、Curvelet变换及本文提出方法的差图像.Wavelet去噪阈值算法经典中有4种选择：(1)Stein无偏似然估计阈值；(2)全局阈值 ，其中σ为标准差，N为对应层像素点数，在含噪记录中σ未知，一般通过第一层分解的高频子带HH₁的中值来获得，即σ=median(HH₁)/0.6745；(3)混合型阈值(‘heursure’阈值)，这种阈值算法是前两者算法的一种折中方式；(4)最小最大准则阈值.

图 5

Fig. 5

图 5 (a)(b)为原始信号和信噪比为-5dB含噪信号; (c)(d)(e)分别为小波变换、Curvelet变换以及 本文方法的滤波结果; (f)(g)(h)分别为小波变换、Curvelet变换以及本文方法的差图像Fig. 5 (a) Raw seismic record with two reflection events; (b) Record with the SNR of -5 dB; (c),(d)and(e) Results after the conventional Wavelet transform, Curvelet transform and steerable filters; (f),(g) and (h) Difference diagrams of conventional wavelet transform, curvelet transform and steerable filters

在实验中，我们将这4种算法均用于测试，由于第1种和第2种阈值算法选择规则相对保守，去噪效果较差，第4种阈值算法出来的波形产生了一定程度的畸变，相对而言，第2种全局阈值的算法有效的去除了噪声并保持了原始信号的特征.因此本文实验最终采用全局阈值算法，对信号进行三层分解，通过阈值进行噪声压制.Curvelet变换应用的是林春和王绪本(2009)提出的算法，采用硬阈值去噪，T=Cσ_λσ，其中σ为噪声的估计标准偏差，σ_λ是每个曲波系数标准偏差近似值，C是一个依赖于尺度和方向的常数，不同的尺度和方向下去不同的值.从滤波结果的视觉效果可看出，本文方法压制背景噪声效果更为理想，差图像可以看出信号保幅性也具有较大的优势.图 6a给出了单道时域波形的对比，图 6b给出了其频域波形的对比，无论从时域效果还是频域性能来看，经本文方法处理后的信号都更接近原始信号.

图 6

Fig. 6

图 6 (a)单道时域波形对比(b)单道频谱对比Fig. 6 (a) Waveforms of signal in the 102nd channels. (b) Spectrum of signal in 102 channel

图 7a绘制了在-15 db到5 db噪声下，数据间隔为1 db，采用三种方法对含噪记录进行处理的结果曲线，从图中可以看出，在不同信噪比记录中，本文方法较传统小波变换和Curvelet去噪后的记录的信噪比均有所提高.图 7b给出了三种方法在信噪比分别为-15 db到5 db的处理结果的平均幅值保持，从图中可以看出，本文方法保幅性明显存在极大优势.从图 7b中可以看出，噪声对保幅性影响不大，算法决定了幅值保持度.

图 7

Fig. 7

图 7 (a)不同噪声背景下的去噪信噪比(b)不同噪声背景下的平均幅值保持Fig. 7 (a) SNR after processing by different methods under different background noise， (b) Keeping average amplitudes after processing by different methods under different background noise

图 8

Fig. 8

图 8 (a)原始地震记录(b)本文方法处理结果Fig. 8 (a) Original seismic records， (b) Results of processing by proposed method

图 9

Fig. 9

图 9 (a) 局部原始地震记录; (b) 本文方法处理结果Fig. 9 (a) Local original seismic records; (b) Results of processing by proposed method

图 10

Fig. 10

图 10 实际频谱与处理后频谱对比Fig. 10 Comparison of spectra of original seismic records and processed signals.

本文首次将拉伸厄米特高斯核函数拟合的方向可控滤波器用于地震信号去噪处理，这种算法出发点在于根据同相轴具有空间几何特性，利用拉伸厄米特高斯核函数具有明确的方向性，使得构造出来的方向可控滤波器具有高方向分辨率并且能够准确定位轴的位置，在减少计算量的同时沿着轴向进行滤波处理，有效地压制随机噪声，保持信号的幅值.本文算法使处理后的地震记录信噪比明显提升，地震同相轴更清晰、连续性更好.方向滤波后的数据选择是根据统计特征进行局部阈值处理，今后，将进一步研究如何对数据进行自适应选择，取得更好的滤波效果.