地震资料去噪是地震资料处理中重要环节之一, 而变换域去噪则是地震资料去噪的主要方法之一^[1]。变换域去噪是通过某一数学变换方法将地震资料变换至对应的变换域中, 并根据有效信号和噪声在变换域中的特点, 去掉噪声并保留有效信号。频率域去噪^[2]是最常用的去噪方法之一, 但该方法无法压制与有效信号频带范围重合的噪声。小波变换去噪^[3]能够将信号进行多尺度分解后再去噪, 但该方法只能沿单一方向进行信号变换, 无法适应同时在多个方向上变化的信号。脊波变换、曲波变换^[4-5]、Shearlet变换^[6]、S变换^[7]等都是小波变换的一种衍生方式, 且都能够处理在多个方向上变化的信号, 因此在地震资料去噪中得到了广泛应用。基于时频分析的算法, 如结合经验模态分解和独立分量分析等^[8-10], 在实际地震资料处理中也取得了不错的应用效果, 但这些算法都是基于图像处理发展起来的, 并不针对地震资料。

为了得到更加适用于地震数据的变换方法, FOMEL等^[11-12]提出了一种全新的数学变换, 即Seislet变换。Seislet变换也属于小波变换类, 但Seislet变换是沿着地震资料同相轴的方向进行变换处理, 而沿测线方向的小波变换可以看作是零倾角的Seislet变换。小波变换、曲波变换等主要以信号采样点为基本变换单位, 而Seislet变换则以地震道为变换单位, 因此Seislet变换能更好地适应地震数据, 且能进一步应用于地震资料去噪。刘洋等^[13]将高阶Seislet变换应用于地震数据随机噪声处理, 使随机噪声得到了有效压制。勾福岩等^[14]利用Seislet变换压制海洋地震资料的涌浪噪声, 取得了较好的应用效果。为获得更好的去噪效果, 刘财等^[15]首次将阈值函数引入到Seislet变换压制地震资料噪声, 但阈值函数本身会影响地震资料去噪效果。

1995年, DONOHO^[16]最早提出了小波阈值去噪算法, 介绍了硬阈值函数和软阈值函数。硬阈值函数在阈值处存在断点, 软阈值函数则存在恒定的偏差^[17]。针对这一缺陷, 一些学者提出了新的阈值函数。GAO等^[18]利用半软阈值函数改进软阈值函数, 并取得了一定的效果。许丽群^[19]进一步提出了硬软阈值函数, 能够一定程度地改善硬、软阈值函数的缺陷。何柯等^[20]在常规阈值函数的基础上提出一种补偿阈值函数算法, 对之前受到损伤的与噪声曲波系数重叠的弱信号进行补偿。余江奇等^[21]根据数据噪声在曲波域的特点, 采用了一种新的阈值函数, 并将其引入到曲波变换稀疏反褶积中。张华等^[22]则在曲波域中的每一个尺度选取不同的阈值因子, 达到曲波域中局部阈值去噪的目的。陈浩等^[23]针对Shearlet变换阈值去噪方法不能随尺度和方向变化的不足, 提出了随尺度和方向变化的自适应阈值。

前人研究表明, 传统阈值去噪算法不能有效地分离差异较小的有效信号与噪声, 且传统阈值确定准则也不能很好地适用于Seislet域。本文在前人研究的基础上, 应用Riemann-Liouville分数阶积分理论推导出分数阶阈值函数; 再根据地震数据在Seislet域不同尺度的分布特点, 提出基于尺度加权的阈值确定方法; 最后利用合成加噪模型数据和实际地震数据验证了本文方法的可行性和有效性。

1 方法原理 1.1 Seislet变换基本原理

Seislet变换主要是利用小波提升算法^[24](见附录A), 并应用平面波解构滤波器(PWD)^[25], 沿着地震同相轴的方向对地震数据进行分解。

PWD-Seislet变换定义的预测算子P和更新算子U如下:

$ P[e]_{i}=\frac{S_{i}^{+}\left[e_{j, i-1}\right]+S_{i}^{-}\left[e_{j, i}\right]}{2} $

(1)

$ U[r]_{i}=\frac{S_{i}^{+}\left[r_{j, i-1}\right]+S_{i}^{-}\left[r_{j, i}\right]}{4} $

(2)

式中:e为偶数序列; r为奇数序列; S_i⁺表示沿着地震同相轴向左进行预测, 数值上等于PWD向左求得的局部倾角; S_i^-表示沿着地震同相轴向右进行预测, 数值上等于PWD向右求得的局部倾角; j是分解层数; i为序列序号。

Seislet变换不再像小波变换沿着单一角度进行分解, 而是根据地震同相轴的方向沿着整个地震数据空间进行分解, 能更为精细地刻画地下地质特征。因此有效信号在Seislet域中能被很好地压缩, 其Seislet系数集中在较低的尺度。

1.2 分数阶阈值算法

阈值去噪算法最早应用于小波变换去噪, 现在广泛应用于曲波变换去噪、Shearlet变换去噪等方法中。稀疏变换将信号投影到对应的变换域, 有效信号在合适的变换域下具有稀疏特征, 即有效信号的能量集中在有限的变换域系数范围, 而噪声则遍布整个变换空间。一般认为代表有效信号的变换系数幅值大于代表噪声的变换系数幅值, 在确定了阈值之后, 大于阈值的变换系数被认为是由有效信号变换而来的, 小于阈值的变换系数则被认为是由噪声变换而来的, 这时可以采用阈值函数对变换系数进行处理, 最后将处理后的系数进行反变换, 达到压制噪声、保留有效信号的目的。

传统的阈值函数主要有硬阈值函数和软阈值函数。

硬阈值函数是将系数绝对值小于阈值的系数直接置为0, 将系数绝对值大于等于阈值的系数保持原值, 但在阈值处存在断点, 如图 1a所示。硬阈值函数如下:

$ T(s)=\left\{\begin{array}{ll} {s} & {|s| \geqslant \lambda} \\ {0} & {|s|<\lambda} \end{array}\right. $

(3)

式中:s为Seislet稀疏变换后的系数; λ为阈值。

软阈值函数将系数绝对值大于阈值的系数等于阈值与该系数的差值, 这样就能确保阈值函数的连续性, 但软阈值函数处理后的系数存在较大的偏差, 重构信号与原信号差异较大, 如图 1b所示。软阈值函数如下:

$ T(s)=\left\{\begin{array}{cc} {\operatorname{sgn}(s)(|s|-\lambda)} & {|s| \geqslant \lambda} \\ {0} & {|s|<\lambda} \end{array}\right. $

(4)

式中:sgn(·)为符号函数。

硬阈值函数在阈值附近存在断点, 重构信号会产生伪吉布斯现象, 一般较少采用。常用的主要是软阈值函数, 但软阈值函数处理后的系数与原系数相比具有恒定的偏差。近几年来, 分数阶微积分理论在科学和工程学的很多领域得到了应用。本文引入Riemann-Liouville分数阶积分理论^[25-27], 其公式如下:

$ J^{\alpha} f(t)=\frac{1}{\Gamma(\alpha)} \int_{0}^{t} \frac{f(\tau)}{(t-\tau)^{1-\alpha}} \mathrm{d} \tau $

(5)

式中:J^αf(t)表示对函数f(t)进行分数阶积分处理; $\Gamma(\alpha)=\int_{0}^{\infty} y^{\alpha-1} e^{-y} \mathrm{d} y$为Gamma函数; α为分数阶阶数。

将分数阶积分理论应用到阈值函数中, 推导得到分数阶阈值函数:

$ \left\{\begin{array}{l} {T(s)=J^{\alpha} f(s)=\frac{1}{\Gamma(\alpha)} \int_{0}^{s} \frac{f(t)}{(s-t)^{1-\alpha}} \mathrm{d} t} \\ {f(s)=\left\{\begin{array}{ll} {s} & {|s| \geqslant \lambda} \\ {0} & {|s|<\lambda} \end{array}\right.} \end{array}\right. $

(6)

式中:f(s)为某一载体函数, 此处用硬阈值函数作为载体函数; s为Seislet变换系数。经过多次测试, 本文选用的分数阶阶数α为0.5阶。被该分数阶积分公式积分后的阈值函数T(s)既能增强函数本身的连续性, 又能保留原函数的特征。公式(6)的数学意义相当于对硬阈值函数做0.5阶的分数阶积分处理, 在保留硬阈值函数原有特征的基础上, 进一步增强其在阈值处的连续性, 进而得到分数阶阈值函数。与传统的硬阈值函数和软阈值函数相比, 分数阶阈值函数很好地控制了系数偏差, 在端点处也得到了一定的平滑, 如图 1c所示。

阈值是有效信号与噪声存在区别的临界点。阈值越小, 噪声压制的程度越小, 阈值越大, 噪声压制的程度越大, 所以阈值的选择会直接影响阈值去噪算法的效果^[16]。噪声的Seislet系数与有效信号的Seislet系数在Seislet域的不同尺度不同, 所以需要对每一个尺度的噪声Seislet系数和有效信号Seislet系数进行分析, 以此计算阈值。本文根据地震数据在Seislet域的分布特点, 并在DONOHO给出的全局阈值^[28]基础上提出了一种尺度加权阈值, 尺度越小, 权重越大, 尺度越大, 权重越小, 且得到的阈值在每一个尺度上都不同, 公式如下:

$ \left\{\begin{array}{l} {\lambda_{j}=\frac{\sigma \sqrt{2 \ln \left(N_{j} N_{i}\right)} w}{\max (j)-j}} \\ {\sigma=\frac{\operatorname{median}\left(\left|s_{j, i}\right|\right)}{0.6745}} \end{array}\right. $

(7)

式中:σ是噪声方差估计公式^[29]; max(·)表示取最大值; max(j)-j为权重; N_j表示尺度的层数; N_i表示每一层系数的个数; j为尺度; i为尺度上每个系数的序号; w为调节参数, 最大不超过分解的最大尺度, w越大对噪声压制的程度越大; s_{j, i}为不同尺度上的Seislet系数; λ_j为不同尺度上的阈值; median(·)表示取中值。

1.3 Seislet分数阶阈值算法流程

Seislet分数阶阈值算法流程见图 2。该流程使用平面波解构滤波器(PWD)来估计地震数据的局部倾角, 也可以使用线性Randon变换、预测误差滤波器、有限差分法等方法来估计地震数据的局部倾角^[30]。在该流程中, 尺度加权阈值的设置非常重要, 第一次设置的阈值w为最大尺度的一半, 再根据去噪效果的好坏进行调节, 直到达到最好的去噪效果。尺度加权阈值利用Seislet域低尺度中有效信号分量远多于高尺度中有效信号分量的特点, 在Seislet域中设置更为合理的阈值, 也可以根据去噪效果人为调控参数w, 得到更好的去噪效果。去噪效果的好坏由人为抽取的单道振幅频率、去噪结果剖面以及压制的噪声进行综合分析得出。

2 数值测试 2.1 反射波双曲线模型测试

为了验证本文提出的Seislet分数阶阈值算法的可行性和有效性, 首先利用共炮点反射波时距曲线公式直接模拟得到反射波双曲线模型。地震反射波时距曲线公式如下:

$ t=\sqrt{\frac{(2 h)^{2}+x^{2}}{v^{2}}} $

(8)

式中:x是炮检距; v为地下介质速度; h为深度。

该模型共有2条反射曲线, 如图 3a所示, 道间距为10m, 共计512道, 时间采样间隔为2 ms, 采样时间为1.024s。

利用公式(9)计算信噪比, 将添加噪声的反射波双曲线模型(图 3b)信噪比R_SN控制在0.25dB。

$ R_{\mathrm{sN}}=\frac{\sum\limits_{i} \sum\limits_{j} S_{i, j}}{\sum\limits_{i} \sum\limits_{j} N_{i, j}} $

(9)

式中:S_{i, j}表示有效信号第i道第j个采样点的振幅值; N_{i, j}表示添加的随机噪声第i道第j个采样点的振幅值, 计算结果的单位为dB。

进行Seislet变换之前, 需要获得局部倾角估计数据。本文使用平面波解构滤波器(PWD)对含噪模型(图 3b)进行局部倾角估计, 得到的局部倾角估计值如图 3c所示, 由图 3可知, 因为噪声的幅值较大, 对局部倾角的估计造成一定的影响, 能量团的分布较为散乱, 但在局部倾角估计图中也能够识别出反射波曲线。图 3d为含噪模型的Seislet系数分布, 在尺度小于50的范围内系数值比较集中, 主要为有效信号的Seislet系数, 但在整个尺度范围都有Seislet系数的分布, 说明随机噪声的Seislet系数在Seislet域中仍然分布杂乱, 难以集中。在模型测试中, 分数阶阈值算法和传统硬、软阈值算法都采用全局阈值。

对含噪模型进行Seislet变换后, 分别采用硬阈值算法、软阈值算法、分数阶阈值算法去噪, 得到不同的去噪结果(图 4)。图 4a为经过硬阈值去噪后的结果, 噪声得到部分衰减, 但还存在部分噪声。经过软阈值算法去噪后(图 4b), 大部分噪声得到了压制, 整体信噪比高于使用硬阈值算法去噪的结果。采用分数阶阈值算法的去噪结果中(图 4c), 噪声几乎全被压制, 去噪效果明显优于前2种算法。因为局部倾角估计算法受到噪声影响导致估计结果不准, 所以图 4c的去噪结果存在部分扭曲现象。采用公式(9)计算去噪后数据的信噪比, 并以其衡量3种算法的去噪效果发现, 硬阈值和软阈值算法去噪结果的信噪比分别为1.72dB和2.13dB, 而分数阶阈值算法去噪结果的信噪比为3.82dB, 去噪后的信噪比明显高于传统硬、软阈值算法。

为了进一步验证方法的有效性, 对3种算法去噪后的结果、加噪模型数据以及未加噪模型数据分别抽取第256道进行频谱分析, 结果如图 5所示, 可以看出, 未加噪数据(图 5a)的主频在50Hz左右, 幅值由中间向两边降低, 呈现近似的正态分布; 加噪数据(图 5b)的频率在整个频带范围都呈现出一种无规律的分布, 高频率分量较多; 虽然图 5c、图 5d、图 5e的高频分量都有所减少, 但图 5e的高频分量最少, 且其主频率成分与图 5a的最为接近。综合分析图 4和图 5可知, 传统硬、软阈值算法对随机噪声的压制效果不甚理想。

2.2 实际地震资料测试

使用某陆上工区的实际单炮地震记录(图 6a)进行测试。该单炮记录的道间距为10m, 共128道, 采样间隔为1 ms, 采样时间1.4s。由图 6a可知, 该单炮记录浅部能量和初至能量较强, 底部能量较弱但存在异常振幅, 整个单炮记录上存在一定的随机噪声和线性噪声。图 6b为该单炮记录的局部倾角估计, 由图可见, 能量团的分布层次分明, 如果某块区域能量团颜色较为一致, 则说明这块区域的局部倾角值大致相同。利用图 6b所示局部倾角估计结果进行Seislet变换得到Seislet系数(图 6c)。在实际地震资料测试中, 首先在分数阶阈值算法、硬阈值算法、软阈值算法中采用尺度加权阈值来进行地震资料去噪, 然后再在分数阶阈值算法中分别采用尺度加权阈值和全局阈值进行去噪效果对比。

对实际单炮记录进行Seislet变换后, 再进行硬阈值算法、软阈值算法、分数阶阈值算法进行去噪, 结果见图 7。经过硬阈值算法去噪后的结果(图 7a)能够在一定程度上压制随机噪声, 但整体去噪程度不如图 7b和图 7c。软阈值算法(图 7b)和分数阶阈值算法(图 7c)都能较好地压制噪声, 在局部区域的同相轴连续性和噪声压制方面, 分数阶阈值去噪算法的效果优于软阈值函数。将未去噪的单炮数据分别和3种不同方法去噪后的结果相减, 得到的噪声记录如图 8所示。硬阈值算法去除的噪声相对较少, 且未能压制一些振幅较强的噪声(黄色箭头); 软阈值算法压制的噪声中却去掉了一定的有效信号(红色圈和红色箭头), 且软阈值算法去掉了部分初至信息; 分数阶阈值算法的去噪效果比前两种阈值算法更好, 对有效信号的伤害最轻。在分数阶阈值算法去噪中分别采用尺度加权阈值和全局阈值, 并对去噪结果和压制的噪声作对比分析(图 9)发现, 采用全局阈值的分数阶阈值去噪算法未能压制底部的异常振幅(黄色箭头), 且同相轴附近仍存在部分噪声(红色箭头), 而采用尺度加权阈值的分数阶阈值去噪算法能更有效地压制噪声。

为进一步验证分数阶阈值去噪算法的有效性, 分别对实际单炮记录、硬阈值算法、软阈值算法、分数阶阈值算法的去噪结果抽取第64道数据进行频谱分析, 结果如图 10所示。图 10a为实际单炮记录的第64道频谱, 可见存在一定的高频分量, 但幅值不高, 频率抖动次数较多, 这些高频分量为对应的背景噪声, 主频应在50~80Hz。硬阈值算法处理后(图 10b), 频谱的高频分量(大于50Hz)并没有太过明显的改观。相较于硬阈值算法, 软阈值算法处理结果的单道频谱在50Hz和80Hz处的幅值均有所下降, 且趋于平缓(图 10c)。结合硬、软阈值算法去噪的结果(图 7a和图 7b)可见, 软阈值算法的去噪效果优于硬阈值算法。经过分数阶阈值算法去噪后的频谱(图 10d)的高频分量幅值低于软阈值算法去噪后的结果, 且频率变化也较为缓慢(见红色箭头), 再结合不同阈值算法的去噪结果(图 7c)可见, 分数阶阈值去噪算法能够有效地压制背景噪声, 且效果优于硬阈值去噪算法和软阈值去噪算法。

3 结论

本文在Riemann-Liouville分数阶积分公式的基础上得到分数阶阈值函数, 增强了阈值处的连续性, 进而提出一种用于Seislet域去噪的尺度加权阈值函数, 实现了阈值计算的局部精细化, 应用了尺度加权阈值的分数阶阈值算法能有效压制噪声并降低有效信号的损失。分别采用硬阈值算法、软阈值算法、分数阶阈值算法对理论模型数据和实际地震资料进行处理, 通过对去噪效果、信噪比数据以及频谱对比分析, 结果表明:传统硬、软阈值算法噪声压制不彻底, 且不能有效分离差异较小的有效信号和噪声; 相较于传统阈值算法, 分数阶阈值算法能更好地压制噪声, 减少对有效信号的损伤; 与全局阈值相比, 本文提出的尺度加权阈值在Seislet域中具有更好的适应性, 在阈值去噪中效果更佳。

虽然Seislet变换能够对较为复杂的地震数据进行变换处理, 但Seislet变换较为依赖局部倾角的估计结果, 一旦信噪比过低, 便会影响Seislet变换的准确性。如何提高局部倾角的估计精度或者减少Seislet变换对局部倾角数据的依赖是今后的改进方向。

附录 A 小波提升算法

小波提升算法是将一个小波通过一系列步骤构建出一个新的小波, 主要分为分裂、更新、预测等环节, 具体步骤如下。

1) 将该序列y_j划分为奇序列和偶序列。

$ \left\{\begin{array}{l} {y_{j}=o_{j+1}+e_{j+1}} \\ {o_{j+1}=y_{j \cdot 2 m+1}} \\ {e_{j+1}=y_{j \cdot 2 m}} \end{array}\right. $

(A1)

式中:y_j是初始序列; o_j+1是奇序列; e_j+1是偶序列; o_j+1和e_j+1为初始序列采样点个数的一半。

2) 用预测算子对偶序列进行预测, 再用奇序列减去偶序列的预测结果, 得到残差序列, 也就是小波系数r_{j+1, m}。此时的残差序列可以代表序列的高频信息。

$ r_{j+1, m}=o_{j+1}-p\left(e_{j+1}\right) $

(A2)

式中:r_{j+1, m}是残差序列; p是预测算子。

预测算子p是求偶数序列e_j+1中两个相邻数的均值, 其定义如下:

$ p\left(e_{j+1, k}\right)=\frac{e_{j , k}+e_{j \cdot k-1}}{2} $

(A3)

3) 再对偶序列进行更新, 能够得到一个序列中的相关信息分量, 也就是小波变换尺度系数。

$ c_{j+1, m}=e_{j+1}+u\left(r_{j+1, m}\right) $

(A4)

式中:c_{j+1, m}为近似序列(尺度系数), 反映序列的低频信息; u为更新算子, 定义如下:

$ u\left(r_{j+1, k}\right)=\frac{r_{j+1, k}+r_{j+1, k+1}}{4} $

(A5)

4) 将近似序列c作为新的输入数据, 在下一个尺度j+2上重复步骤1)~步骤3)。