自适应结构方向复扩散压制沙漠地震随机噪声

林红波, 张丹丹

(吉林大学 通信工程学院, 长春 130012)

0 引 言

在地震勘探中, 所获得的地震勘探数据往往伴随着大量的随机噪声, 严重影响地震数据处理的精度, 给进一步的地质解释工作带来很多困难。因此, 压制地震勘探随机噪声, 提高地震勘探数据质量是地震勘探数据处理中重要的基础环节。近些年来, 学者们提出很多地震勘探随机噪声压制方法, 如维纳滤波[1]、经验模式分解[2]、奇异值分解[3]、小波变换[4]和曲波变换[5]等。这些方法在地震勘探随机噪声压制方面取得了很好的效果。但在低信噪比的情况下, 有效信号的恢复不够理想, 尤其在沙漠地震勘探数据中, 沙漠随机噪声频率低, 能量强, 与有效信号具有相似的时空特性和频率特性, 使辨识和提取微弱信号变得更加困难。

基于PDE(Partial Differential Equation)的多尺度分析方法在图像去噪方面效果显著, 最早的PDE模型是实数域线性扩散方程[6], 扩散方程的解等同于含噪图像和时变标准偏差的高斯核的卷积, 扩散系数为常数。该方法能有效平滑随机噪声, 同时也模糊了一些边缘和细节。在此基础上, 非线性扩散模型构建了基于梯度变化的扩散函数, 有效改善了图像结构信息的保留效果[7-10]。非线性扩散滤波方法不依赖先验知识, 如频率、线性、低秩以及噪声统计特性等, 因此在复杂随机噪声压制方面具有明显的优势。然而, 这种实数域扩散方法利用梯度辨识边缘并调整扩散强度, 在低信噪比情况下易混淆边缘和噪声。基于这个问题, Gilboa等[11]提出复数域扩散方法。非线性复扩散(NCD: Nonlinear Complex Diffusion)的虚部为图像平滑的二阶导, 能抑制噪声对边缘检测的影响, 扩散系数是基于虚部的单调递减函数, 相比于基于梯度变化的扩散系数能更好适应复杂的结构特征。基于NCD的算法适应低信噪比数据和复杂结构, 成功应用于地震勘探数据去噪[12]。然而, 该扩散方法仍然利用梯度描述同相轴的局部特征, 无法准确刻画沙漠地震勘探数据中具有复杂结构同相轴的方向特征, 易导致同相轴衰减。笔者提出了自适应结构方向复扩散算法(ASOCD：Adaptive Structure-Oriented Complex Diffusion), 通过在复扩散方程中引入基于结构方向的方向导数, 使扩散过程沿同相轴方向演化, 更好地保留同相轴；此外, 构建了基于局部协方差矩阵的自适应扩散阈值, 有效地解决了沙漠随机噪声压制不彻底问题。合成数据和实际数据处理验证了ASOCD方法的有效性。

1 线性复扩散与非线性复扩散

1.1 线性复扩散

Gilboa等[11]受到量子力学启发, 将实数域的线性扩散和自由薛定谔方程结合得到线性复扩散。自由薛定谔方程定义为

(1)

线性扩散方程定义为

(2)

其中c是常扩散系数, 初始条件为u(x)|t=0=u0(x),u0(x)为初始的含噪数据。基于自由薛定谔方程, 将扩散方程推广到复数域, 得到线性复扩散方程

(3)

其中c=rejθ,r是c的模。当θ趋近于0时, 式(3)的基本解可以表示为

(4)

1.2 非线性复扩散

线性复扩散的常扩散系数会导致在去除噪声的同时平滑掉边缘。为解决这个问题, Gilboa等[11]引入以虚部为变量的扩散系数, 提出了非线性复扩散NCD。对于含噪数据u∈R2, 非线性复扩散定义为

(5)

其中扩散系数c是关于虚部单调递减的函数, 扩散阈值k和虚部共同作用调节扩散强度。为了应用复扩散算法, 需要在离散形式下求解, 用u(p,q)表示地震勘探数据u的第q道的第p个采样点。对式(5)的离散化, 笔者采用有限差分离散化方法, 得NCD方程离散化公式为

(6)

(7)

从式(7)中可看出, NCD方法利用相邻点表征式(5)的中的u, 从而实现扩散沿着含噪数据的边缘方向和梯度方向, 并结合扩散强度实现噪声的抑制以及边缘的保留。

2 自适应结构方向复扩散原理

由于接收沙漠地震数据的检波器间隔较远, 相邻检波器接收到反射子波的时间延迟比较长, 因此同相轴往往具有较大的倾角。而NCD方法仅利用相邻点表征结构信息, 所描述的边缘方向和沙漠地震数据的同相轴方向不能吻合, 会导致同相轴的失真。另外, 部分沙漠地震随机噪声和有效信号具有相似的虚部, 对信号和噪声采用相同的扩散阈值, 无法在保留有效信号的同时抑制随机噪声。针对这些问题, 笔者提出自适应结构方向复扩散(ASOCD)方法, 构造基于结构方向的方向导数, 使扩散能沿同相轴方向演化, 并根据区域特性自适应调节扩散阈值, 从而在抑制噪声的同时增强有效信号。

在离散NCD模型的基础上, 引入基于结构方向的方向导数, 得到自适应结构方向复扩散模型

(8)

图1 扫描示意图

其中cT表示结构方向的扩散系数,cN表示为垂直于结构方向的扩散系数,DR(un)表示结构方向的方向导数估计,DL(un)表示与结构方向相反方向的方向导数估计。如图1所示, 方向导数估计为沿同相轴方向求差分, 定义为

(9)

其中sp为结构方向的斜率。式(8)表明, 自适应结构方向复扩散能沿着同相轴方向定向扩散, 扩散强度由扩散系数调节。为了实现对随机噪声的抑制, 需要增强噪声区域的扩散系数cT和cN；而为了更好地恢复同相轴, 则信号区域的扩散系数cT要比噪声区域小, 并尽可能减少扩散系数cN。为此, 笔者构建了自适应复扩散系数。

2.1 自适应复扩散系数

基于地震勘探数据的区域特征, 自适应地调整ASOCD模型的扩散系数。结构方向的扩散系数cT定义为

cT(Im(u(p,q)))=ejθ/(1+(Im(u(p,q))/k′(p,q)θ)2)

(10)

其中阈值k′(p,q)能在信号区域Ω和噪声区域φ自适应调整, 定义为

(11)

满足Kn>Ks, 垂直与结构方向的扩散系数cN定义为

(12)

由式(10)和式(12)可知, 结构方向复扩散模型在噪声区域和信号区域具有不同的扩散强度。沙漠随机噪声具有弱相似性, 和有效信号具有相似虚部, 此时利用较大的阈值Kn可以增强扩散强度, 能解决常阈值扩散滤波对噪声压制不彻底的问题。在信号区域Ω, 结构方向复扩散只进行结构方向的定向扩散, 在同相轴方向采用较小阈值控制扩散强度, 在垂直同相轴方向的扩散系数为0, 避免同相轴模糊。

笔者利用局部协方差矩阵辨识信号区域Ω和噪声区域φ。由于有效信号存在较强的空间相关性, 沙漠地震勘探随机噪声没有或存在较弱的空间相关性, 而协方差矩阵能反映矩阵中任意两个变量之间的相关性[13], 则能反映在一定窗长内的任意两道地震数据之间的相关性, 因此能识别信号区域和噪声区域。对于点u(p,q)的局部协方差矩阵Fp,q的计算, 以该点为中心, 计算大小为A×B的邻域L的协方差矩阵, 则局部协方差矩阵Fp,q中任意一点Fp,q(p′,q′)可表示为

(13)

其中up′表示矩阵L中的第p′列,μp′表示矩阵L中p′列的均值, T代表转置。局部协方差矩阵最大的值在对角线上取得, 对于噪声区域, 局部协方差矩阵对角线上的值即是噪声方差, 对于信号区域, 局部协方差矩阵对角线上的值即是含噪信号的方差。可以设置区分信号和噪声的阈值τ, 表示为

τ=σuσn

(14)

其中σn是利用Chen等[14]提出的噪声估计方法估计出的噪声标准差,σu表示含噪数据的标准差。由于含噪信号方差大于噪声方差, 则对于信号区域, 局部协方差矩阵对角线上的值要比阈值大, 对于噪声区域, 局部协方差矩阵对角线上的值要比阈值小。则识别的信号区域Ω和噪声区域φ表示为

(15)

2.2 结构方向估计

结构方向复扩散模型基于方向导数沿地震勘探数据的结构方向自适应扩散, 其中结构方向的估计至关重要。基于地震勘探同相轴的空间相关性, 笔者采用非局部统计特性扫描方法估计地震勘探数据的结构方向。由于一部分沙漠随机噪声具有弱空间相似性, 若对随机噪声也进行方向扩散, 易产生具有一定空间相似性的噪声残留。为更好地压制随机噪声, 笔者在估计结构方向时增加约束条件, 使噪声区域的结构方向为水平方向, 此时结构方向扩散退化为非线性复扩散, 并通过自适应调节扩散系数增强对沙漠随机噪声抑制能力。

(16)

(17)

3 仿真实验

3.1 模拟记录的处理结果

为了验证ASOCD方法的有效性, 对合成的沙漠地震记录进行处理。合成的沙漠地震记录包含由Ricker构成的3条同相轴, 主频率从上到下依次为20 Hz, 18 Hz和16 Hz。采样频率是1 000 Hz, 如图2a所示。添加真实的沙漠随机噪声, 使信噪比达到-5.46 dB(见图2b)。将ASOCD滤波方法应用于含噪数据处理并和NCD方法进行对比, 去噪结果以及滤除的噪声如图2c～图2f所示。可以看出, 这两种方法都能抑制噪声, 但ASOCD滤波方法能获得更干净的背景。在信号保留方面, NCD方法滤除的噪声中可以看出同相轴, 则NCD方法在抑制噪声的同时损失了有效信号, 而ASOCD算法去除的噪声中几乎不存在有效信号, 因此ASOCD算法能够抑制噪声的同时保留有效信号。

a 纯净信号 b 含噪信号 c NCD去噪结果

d NCD滤除的噪声 e ASOCD去噪结果 f ASOCD滤除的噪声

对去噪结果进行量化分析, 表1给出了不同沙漠噪声水平下, 去噪前后数据的信噪比比较, 其中使用信噪比公式为

(18)

其中us和ud分别为原始无噪数据和去噪后的数据, 数据结果如表1所示。从表1中可看出, 不同噪声水平下, ASOCD算法去噪后的信噪比相比于传统NCD算法平均提高了6 dB。

表1 不同沙漠噪声水平下的信噪比对比

3.2 实际记录的处理结果

为了进一步验证提出方法的有效性, 对真实的沙漠地震数据进行处理。图3a为中国沙漠地区包含107道和1 400个采样点的真实沙漠地震数据。显然有效信号被噪声严重损坏了, 某些噪声区域的幅值甚至超过了有效信号。对该地震数据分别应用NCD方法和ASOCD方法去噪, 对比去噪后的结果, ASOCD方法(见图3c)相比NCD方法(见图3b)可以获得更加干净的背景。尤其在黑色框所示区域, ASOCD方法能获得更加清晰更加连续的同相轴, 噪声也能得到更好的抑制。

a 原始地震数据 b ASOCD去噪结果 c NCD去噪结果

4 结 语

笔者提出ASOCD算法对沙漠地震数据中的随机噪声进行压制, 利用非局部特性扫描方法计算同相轴的方向, 并构造基于结构方向的方向导数, 使扩散能沿着同相轴方向进行。此外, 利用局部协方差矩阵提高了信号区域和噪声区域的识别能力, 能更有效区分信号和噪声在低信噪比情况下, 使得扩散阈值能自适应区域特性进行调整。ASOCD算法实现了在噪声区域进行扩散阈值较大的非线性复扩散, 有效地抑制沙漠随机噪声, 在信号区域沿着同相轴方向进行扩散阈值较小的定向复扩散, 有效地保留同相轴。合成的沙漠地震数据以及真实的沙漠地震数据验证了所提算法的有效性。同NCD算法相比, ASOCD算法具有更好的信号恢复能力, 以及更优越的噪声抑制能力。