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

Vol. 37 No. 5

2019年9月

Journal of Jilin University (Information Science Edition)

Sept. 2019

文章编号：1671-56 (2019) 05 ~0463 ~07

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

林红波，张丹丹

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

摘要：非线性复扩散（NCD: Nonlinear Complex Diffusion)能针对不同结构特征调节扩散强度，有效压制地震勘 探随机噪声。但沙漠地震勘探数据结构复杂，NCD方法不能有效估计出同相轴方向，会导致同相轴衰减。 为此，基于同相轴结构特征提出自适应结构方向复扩散（ASOCD: Adaptive Structure-Oriented Complex Diffusion) 滤波算法，该算法沿同相轴方向定向复扩散，在压制随机噪声同时增强有效信号。此外，该算法利用地震勘探 数据的区域特性自适应调节扩散阈值，在信号区域减小扩散阈值以更好地恢复地震勘探信号。合成和实际沙 漠地震数据结果验证了所提方法的有效性，与NCD方法相比，滤波后信噪比提高了 6 dB左右。关键词：非线性复扩散；局部协方差矩阵；定向复扩散；沙漠随机噪声；地震勘探 中图分类号：TN911.7

文献标识码：A

Adaptive Structure-Oriented Complex Diffusion Filtering for

Desert Seismic Random Noise Attenuation

LIN Hongbo, ZHANG Dandan

(College of Communication Engineering, Jilin University, Changchun 130012, China)

Abstract： NCD (Nonlinear Complex Diffusion) can adjust the diffusion intensity according to different structural

features and effectively suppress the random noise in seismic data. However, the morphology of desert seismic data is complex, so that the NCD method is unable to effectively estimate the direction of the reflection event, leading to the attenuation of the reflection event. Aiming at this problem, an ASOCD ( Adaptive Structure- Oriented Complex Diffusion) algorithm is proposed based on the structural features of the reflection event. The oriented complex diffusion is performed along the direction of refection event to suppress the random noise and enhance the effective signal. And the proposed algorithm adaptively adjusts the diffusion threshold according to the regional characteristics of seismic data, thus better recovers the seismic signal using a small diffusion threshold for the signal region. The results of the synthetic and real desert seismic data verify the effectiveness of the proposed method, and the proposed method can increase the SNR ( Signal-To-Noise Ratio) of filtered data by

6 dB more than the NCD method.

Key words ： nonlinear complex diffusion ( NCD) ； local covariance matrix ； oriented complex diffusion ； desert

random nose； seismic exploration

0引言

在地震勘探中，所获得的地震勘探数据往往伴随着大量的随机噪声，严重影响地震数据处理的精

度，给进一步的地质解释工作带来很多困难。因此，压制地震勘探随机噪声，提高地震勘探数据质量是

收稿日期：20丨94741

基金项目：国家自然科学基金资助项目（41774117);吉林省科技发展计划基金资助项目（20170101167JC)作者简介：林红波（1973— ），女，长春人，吉林大学教授，博士生导师，主要从事信号与信息处理研究，（Tel)86-13844167429

(E-mail) hblin@jlu. edu. cn〇

4吉林大学学报(信息科学版)第37卷

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

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

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

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

1.1线性复扩散薛定谔方程定义为

其中幻是量子力学中的波动函数，爪是质量，ft是普朗克常量，4是拉普拉斯算子，j=

线性扩散方程定义为

—(«,*) =cAu(t,x), dt

其中C是常扩散系数，初

始

条

件

为

«

将扩散方程推广到复数域，得到线性复扩散方程

= cAu(t,x) , 0 < c e C

l-u(a：) |(_〇— Uq{x)0 < c e R (2)

。（幻为初始的含噪数据。基于自由薛定谔方程，

⑶

其中c = re'「是〇的模。当0趋近于0时，式(3)的基本解可以表示为

0 —♦ 0

lim Re(u) = g-a*u0,

d —*〇 lim

tf

= tAg^. * u0 (4)

其中G表示标准差为&的高斯核，*表示卷积。式(4)表明，实部Re(u)为含噪数据与高斯核的卷积， 虚部Im(u)是含噪数据与高斯核卷积的拉普拉斯变换（平滑的二阶导数）。由此可看出，复扩散方程解的 实部相当于高斯平滑的滤波结果，复扩散方程解的虚部能表征结构特征，并克服了拉普拉斯算子对噪声 敏感的问题。

第5期

林红波，等：自适应结构方向复扩散压制沙漠地震随机噪声

465

1.2非线性复扩散

线性复扩散的常扩散系数会导致在去除噪声的同时平滑掉边缘。为解决这个问题，

引人

以虚部为变量的扩散系数，提出了非线性复扩散NCD。对于含噪数据ue/?2,非线性复扩散定义为

^ = V • (c(Im(u)) V u)

c(Im(it))

= eie/(l + (lm(u)/k6)2)其中扩散系数c是关于虚部单调递减的函数，扩散阈值A和虚部共同作用调节扩散强度。为了应用复扩 散算法，需要在离散形式下求解，用《(/»，？)表示地震勘探数据u的第9道的第个采样点。对式（5)的 离散化，笔者采用有限差分离散化方法，得NCD方程离散化公式为

u\"+l

= u\" + A[c( Im(u1^) )DN(un) + c(Im(ul) )Ds(un) +c(lm(unz))DE(un) + c(lm(unw))Dw(u\")]

(6)

其中表示第

n

次迭代的去噪结果，A表示时间步长，是由前向差分对如/扣进行离散化得到的，下角

标 N、S、E、W 分别是北南东西 4 个方向的缩写，i^=

u\"(/5-l，g)，i4

=wn(p + l,g)，i4=

u\"(p，g + l

)，

u\"w = u\"(P- 1)，％ (〇、£>s(〇、(f)、(〇是利用中心差分对▽ u进行离散化得到的，分别为

N、S、E、W

4方向的方向导数估计，可表示为

DN(un) = ul - un(p,q) ,

Ds(un) =ul~u{p,q)DE(r

un) =ul~un(p,q), Dw(u\") =uw-un(p,q)(7)

从式(7)中可看出，NCD方法利用相邻点表征式(5)的中的Vu，从而实现扩散沿着含噪数据的边缘方向

和梯度方向，并结合扩散强度实现噪声的抑制以及边缘的保留。

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

由于接收沙漠地震数据的检波器间隔较远，相邻检波器接收到反射子波的时间延迟比较长，因此同

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

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

un+l = un + A[cA,(Im(«^)

+ cN(lm(^u\"s) )Ds(u\") +cT(hn(u\"L) )D,(u\") + cT(lm(unR))DR(un) ]

(8)

其中•表示结构方向的扩散系数，cw表示为垂直于

结构方向的扩散系数

，

表示结构方向的方向导

数估计，久U\")表示与结构方向相反方向的方向导数 估计。如图1所示，方向导数估计为沿同相轴方向求 差分，定义为

Dl{u) = u\"L - u{p,q) uL = un(p + s q - \\)图1 扫描示意图

Dr(u\")

= uR - u\"(p,q)

(9)

Fig. 1 Scanning sketch map

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

466

吉林大学学报(信息科学版)第37卷

2.1自适应复扩散系数

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

cr(Im(u(j3,9) ) )

= e,e/(\\ + {\\m(u{p,q))/k'(p,q)d)1)

K , Ks,

u{p,q) e (p u(p,q) e {1

(10)

其中阈值F(/>,9)能在信号区域和噪声区域自适应调整，定义为

(

11)

满足心>K,垂直与结构方向的扩散系数&定义为

C/v(—

))) ={々(1+(I—

)’以)2)，心，9)

(12)

!■〇， u(p,q) e Q

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

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

= J^Y[(*V-/V)T(Mr/V)]，

其中表示矩阵

i

p' - 1,2,q' = (13)

中的第，列，/V表示矩阵L中y列的均值，T代表转置。局部协方差矩阵最大的值在

对角线上取得，对于噪声区域，局部协方差矩阵对角线上的值即是噪声方差，对于信号区域，局部协方 差矩阵对角线上的值即是含噪信号的方差。可以设置区分信号和噪声的阈值t，表示为

T = (Tu(Tn

(14)

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

|/2 = \\u{p,q)\\Fpq{p' ,q') > T, 3(p',q')\\

10 = \\u(p,q) \\Fpq(p' ,q') ^ r, \\/(p',q')\\⑴)

2.2 结构方向估计

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

随机噪声抑制能力。

对于地震勘探数据u(p,9)，非局部统计特性扫描方法首先设定一个非局部搜索区间，对该区间的每 个点，以当前点为中心，以tp( -M幻f率彳

~G~~

^ (u(p - itp,q + i)=-G

(16)

M

)为斜率的扫描轨线~进行扫描。以M(P，9)为例，如图1所

示，计算轨线与2C + 1个地震道交点的相似性M,p[l5]，并将取彳^最大相似性的斜率作为该点的初始斜

第5期

林红波，等：自适应结构方向复扩散压制沙漠地震随机噪声

467

其中爪>是斜率为*p的轨线和2C + 1个地震道相交的采样点的均值。统计搜素框内所有点的初始斜率， 将得到数目最多的斜率作为点》(P，？）的同相轴斜率l(P,9)。则含噪地震数据的结构方向的斜率\\ 可表示为

sP(p,q). «(p,g) e n

l〇, u(p,q) e (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所示。可以看出，这两种方

法都能抑制噪声，但AS0CD

滤波方法能获得更干净的背景。在信号保留方面，NCD方法滤除的噪声中

可以看出同相轴，则

NCD方法在抑制噪声的同时损失了有效信号，而AS0CD

算法去除的噪声中几乎不

存在有效信号，因此AS0CD算法能够抑制噪声的同时保留有效信号。

100

埯200

2〇〇

张300

张300

400

20

30

40

50

500

0

10

20

30

40

50

地震道地震道

b含噪信号NCD去噪结果

^ 200

味300

400

0

10

20

30

40

50

0

10

20

30

40

50

地震道地震道

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

f AS0CD滤除的噪声

图2

模拟沙漠地震数据的去噪结果对比

Fig. 2 Results of synthetic desert seismic data

对去噪结果进行量化分析，表i

给出了不同沙漠噪声水平下，去噪前后数据的信噪比比较，其中使

用信噪比公式为

Rsm = 101g( || us ||

|| us-ud\\\\2)

(18)

其中和心分别为原始无噪数据和去噪后的数据，数据结果如表1所示。从表1中可看出，不同噪声

水平下，ASOCD算法去噪后的信噪比相比于传统NCD算法平均提高了 6 dB。

表1

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

Tab. 1

Comparison of SNR with different desertnoise levels

/dB

初始尺SNR

NCDASOCD初始‘H

NCDASOCD- 0.834.009.65-5.46-0.325.79-3.33

1.68

7.44

-6.85

-1.98

4.04

468吉林大学学报（信息科学版)

第37卷

3.2 实际记录的处理结果

为了进一步验证提出方法的有效性，对真实的沙漠地震数据进行处理。图3a为中国沙漠地区包含 107道和1 400个采样点的真实沙漠地震数据。显然有效信号被噪声严重损坏了，某些噪声区域的幅值 甚至超过了有效信号。对该地震数据分别应用NCD方法和ASOCD方法去噪，对比去噪后的结果，

ASOCD方法（见图3c

)相比NCD方法（见图3b)可以获得更加干净的背景。尤其在黑色框所示区域，

00.3^ 0.6

I

00.3

-

ASOCD方法能获得更加清晰更加连续的同相轴，噪声也能得到更好的抑制。

0

0.3 ^ 0.60.9

1.2

0

麵

17

39

61

83

道数

0-6

0.9 1.20

17

39

61

道数

^ 0.9

1.2

83

17

39

61

道数

83 107

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

图3 实际沙漠地震数据的去噪结果对比

Fig. 3 Results of real desert seismic data

4 结语

笔者提出ASOCD算法对沙漠地震数据中的随机噪声进行压制，利用非局部特性扫描方法计算同相

轴的方向，并构造基于结构方向的方向导数，使扩散能沿着同相轴方向进行。此外，利用局部协方差矩 阵提高了信号区域和噪声区域的识别能力，能更有效区分信号和噪声在低信噪比情况下，使得扩散阈值

能自适应区域特性进行调整。ASOCD算法实现了在噪声区域进行扩散阈值较大的非线性复扩散，有效 地抑制沙漠随机噪声，在信号区域沿着同相轴方向进行扩散阈值较小的定向复扩散，有效地保留同相 轴。合成的沙漠地震数据以及真实的沙漠地震数据验证了所提算法的有效性。同NCD算法相比,

ASOCD算法具有更好的信号恢复能力，以及更优越的噪声抑制能力。

参考文献：

[1 ]李娟，孟可心，李月，等.基于相似匹配维纳滤波的地震资料噪声压制[J].吉林大学学报：工学版，2017, 47(6):

19-1968.

LI Juan, MENG Kexin, LI Yue, et al. Seismic Data Noise Suppression Based on Similar Matching Wiener Filtering [ J]. Journal of Jilin University： Engineering Edition, 2017, 47(6) ： 19-1968.

[2 ] HOU Wenlong, JIA Ruisheng, SUN Hongmei, et al. Random Noise Reduction in Seismic Data by Using Bidimensional

Empirical Mode Decomposition and Shearlet Transform [J]. IEEE Access, 2019, 7： 71734-71386.

[3 ] LIANG Xiaoqiang, LI Yue, ZHANG Chao. Noise Suppression for Microseismic Data by Non-Subsampled Shearlet Transform

Based on Singular Value Decomposition [J]. Geophysical Prospecting, 2018, 66： 4-903.

[4 ]LIU Shucong, GAO Ergen, CHEN Xun. Seismic Data Denoising Simulation Research Based on Wavelet Transform [J].

Applied Mechanics & Materials, 2014, 490/491 ： 1356-1360.

[5]张岩，任伟建，唐国维•基于曲波变换的地震数据压缩感知重构算法[J].吉林大学学报：信息科学版，2015,

33(5)： 570-577.

ZHANG Yan, REN Weijian, TANG Guowei. Seismic Data Compressed Sensing Reconstruction Algorithm Based on Curvelet Transform [J]. Journal of Jilin University： Information Science Edition, 2015, 33(5) ： 570-577.

[6 ] WITKIN A P. Scale-Space Filtering [C] /^Proceedings of the 7th International Joint Conference on Artificial Intelligence. San

Francisco, CA, USA： Morgan Kaufmann Publishers, 1983： 1019-1022.

[7 ]TEBINI S, MBARKI Z, SEDDIK H, et al. Rapid and Efficient Image Restoration Technique Based on New Adaptive

Anisotropic Diffusion Function [J]. Digital Signal Processing, 2016, 48： 201-215.

第5期

林红波，等：自适应结构方向复扩散压制沙漠地震随机噪声

469

[8 ] CHAO S M, TSAI D M. An Improved Anisotropic Diffusion Model for Detail-and Edge-Preserving Smoothing [ J]. Pattern

Recognition Letters, 2010, 31(13) ： 2012-2023.

[9 ]XU Jiangtao, JIA Yuanyuan, SHI Zaifeng, et al. An Improved Anisotropic Diffusion Filter with Semi-adaptive Threshold for

Edge Preservation [J]. Signal Processing, 2016, 119： 80-91.

[10] RAFSANJANI H K, SEDAAGHI M H, SARYAZDI S. An Adaptive Diffusion Coefficient Selection for Image Denoising [J].

Digital Signal Processing, 2017 , ： 71-82.

[11] GILBOA G, SOCHEN N, ZEEVI Y Y. Image Enhancement and Denoising by Complex Diffusion Processes [J]. IEEE

Transactions on Pattern Analysis and Machine Intelligence, 2004, 26(8) ： 1020-1036.

[12] ZHANG Yushu, LIN Hongbo. The Adaptive Complex Shock Diffusion for Seismic Random Noise Attenuation [C] //2017 IEEE

Global Conference on Signal and Information Processing. Montreal： IEEE, 2018： 333-337.[13] 罗鹏飞，张文明.随机信号分析与处理[M].北京：清华大学出版社，2006.

LUO Pengfei, ZHANG Wenming. Random Signal Analysis and Processing [M]. Beijing： Tsinghua University Press, 2006.[14] CHEN Guangyong, ZHU Fengyuan, HENG P A. An Efficient Statistical Method for Image Noise Level Estimation [C] /^IEEE

International Conference on Computer Vision. Santiago, Chile： IEEE, 2015： 477-485.

[15] 曹善舒，李月，吴宁.自适应径向道TFPF压制地震记录随机噪声[J].吉林大学学报：信息科学版，2016, 34(2):

166-173.

CAO Shanshu, LI Yue, WU Ning. Adaptive Radial TFPF Suppression of Random Noise from Seismic Recording [ J]. Journal of Jilin University： Information Science Edition, 2016, 34(2) ： 166-173.

(责任编辑：张洁）