三维地震数据频域无监督随机噪声压制方法

薛亚茹，苏军利*，冯璐瑜，张程，梁琪

（1.中国石油大学（北京）信息科学与工程学院，北京 102249；2.中国石油大学（北京）油气资源与工程全国重点实验室，北京 102249）

0 引言

地震数据往往受自然环境、采集仪器等因素的影响而被不同程度噪声干扰，从而降低了信噪比和保真度。因此，对采集的地震数据去噪是不可或缺的重要环节。

依据地震数据特征和不同先验假设提出了各种去噪理论和实现方法，主要可分为以下几类。①基于预测的滤波方法，如t-x域预测滤波[1]、f-x域反褶积滤波[2]、非平稳性预测滤波[3]等。该类方法基于地震道之间的可预测性，适用于线性同相轴。②基于稀疏变换的滤波方法，如傅里叶变换[4]、小波变换[5]、曲波变换[6]等。此类方法假设地震数据在变换域能被稀疏表征，通过抑制较小的参数值达到去噪目的，在一定程度上可以获得良好的去噪效果，但变换域噪声压制的阈值等超参数难以确定。 ③矩阵降秩方法。该类方法假设地震数据可以表征为一个低秩矩阵，随机噪声的存在会增加矩阵的秩，因此可以将去噪问题转化为低秩矩阵逼近问题。如Cadzow 滤波[7]、奇异值分解[8]、多通道奇异谱分析(Multichannel Singular Spectrum Analysis，MSSA)[9]等方法。字典学习本质也属于降秩方法，它能够从地震数据中自适应地学习基函数，从而对数据进行稀疏表示，以便将信号与噪声进行分离，达到好的去噪效果[10]，然而其计算成本较高，而且算法在字典更新步骤中会破坏稀疏系数的结构。模态分解也是降秩方法之一，通过提取有限的模态分量达到信噪分离的目的，目前使用最为广泛的有经验模态分解[11]和变分模态分解[12]。

近年来，随着计算机性能的迅速提高及大数据时代的来临，深度学习研究备受青睐。不少学者将深度学习引入到地球物理领域实现对地震数据去噪处理[13-15]。最初的研究多是采用有监督深度学习模型，相较于依赖先验信息的传统算法，有监督模型在地震数据去噪和插值方面表现出了优异的性能[16-18]，但是海量高质量标签数据的获取，无疑增加了人工和计算成本。因此，无监督学习受到了广泛关注。基于无监督学习的自编码网络模型通过对输入数据编码压缩和解码重构达到去噪的目的[19-20]。Zhang 等[21]通过对网络中间层加入稀疏约束提高网络去噪性能；Yang等[22]提出了一种全连接型深度去噪无监督网络用于三维地震数据噪声压制；Saad等[23]将无监督网络用于增强微地震信号和压制随机噪声，通过对原始数据进行切片操作，初步减小输入数据的维度，之后利用自编码网络实现信噪分离。为了强化网络对地震数据高阶特征的学习能力，在编码器和解码器之间加入跳跃连接，使得解码器在获得全局特征的同时可以融入编码器提取的局部特征[24]。考虑到叠前地震资料经动校正(NMO)后具有较高的自相似性，Liu 等[25]提出了一种基于无监督深度学习压制面波的方法，将经过动校正后的地震数据作为神经网络的输入，利用UNet 型卷积网络提取数据特征，从而达到信噪分离的目的。

上述基于神经网络的地震数据去噪方法大多在时空域完成。将地震数据剖分成小的时空切片[24]，利用切片之间的相似性[25]，通过神经网络提取有效信号特征，实现信噪分离。此类方法将数据切分，未曾考虑切片之间的连续性，也未考虑三维数据体之间的相关性。地震数据在频域是可预测的，MSSA 法利用该特性实现噪声压制和数据重构[9]，并能够较好地保持横向连续性，但该方法需要构建大型Hankel 矩阵，计算量大；同时基于奇异值分解理论的压缩去噪方法提取的地震数据特征是相互正交的，对信号进行稀疏表征的程度有限，而且有效信号秩的优选也比较困难。

本文提出一种频域降秩的深度学习去噪方法。首先分析三维地震数据体的低秩特性，利用奇异值分解理论指导自编码网络的建立；针对频域噪声分布特征，引入K-L(Kullback-Leibler)散度优化损失函数。合成记录和实际地震资料处理结果表明，本文方法在随机噪声压制效果和处理效率等方面均优于MSSA和K-SVD(K-奇异值分解)方法。

1 方法原理

1.1 三维地震数据的低秩表征

地震勘探通过多炮激发实现地下构造反射成像。为提高成像精度，震源在探区多次激发，构成了炮点坐标(xs，ys)、检波点坐标(xr，yr)及时间t的五维地震数据，可表示为d(xs，ys，xr，yr，t)。为方便解释，本文以一条测线的炮集为例进行分析，但同样适用其他三维数据体。假设地层是水平层状介质，第k炮数据可表示为dk(xs，xr，t)，则第k+1炮记录可表示为

式中：Δxs为炮间距；p为射线参数。变换到频率域，有

由上式可知，相邻炮记录之间的频谱仅有线性相位差s=exp(-i2πfpΔxs)。对某个频率而言，各炮数据构成低秩矩阵，可以表示为

式中：n为检波点个数；m为炮数。由此可以看到该矩阵的秩是1，也就是说各炮记录是相似的。但是考虑到实际地层并非水平，也并非各向同性，上述矩阵秩不一定是1，但仍满足低秩假设，因此对于三维地震数据, 其频率切片满足低秩特征。

M的奇异值分解可表示为

式中：r为矩阵M的秩；U是矩阵MMT的特征向量矩阵，由各炮f-xr谱特征向量组成，且相互正交；V是矩阵MTM的特征向量矩阵，由各检波点的f-xs谱特征向量组成；S表示各特征值构成的对角矩阵；σk为第k个特征值；uk、vk为第k个特征值对应的特征向量。

根据上述分析，频率切片会呈现一定的低秩特征(或稀疏性)，对任一炮的某一频率分量可以表示为特征向量的稀疏表示，即

式中：mj表示第j炮频率分量；cj=S(VT)j为稀疏表征参数。该表示方式也可以理解为稀疏字典概念，其中U也称为字典。通过适当的降秩，可以达到压制噪声的目的。

奇异值分解方法获得的各特征向量是相互正交的，这限制了数据稀疏表征程度。为了解决该问题，提出了独立分量分析、基于数据驱动的冗余字典等方法。本文将采用数据驱动的神经网络框架，在频域实现地震数据压缩去噪。

1.2 无监督自编码网络模型

自编码器是Lecun[26]提出的一种具有对称结构的无监督型神经网络。对于输入数据为X，编码器函数可表示为

该过程期望获得数据X的稀疏编码。式中：h为隐层神经元的输出；w为网络输入的权重；b为网络输入的偏置。与式(5)相比，实质就是寻找稀疏编码系数c；编码器各神经元的权重系数w近似为式(5)中的特征向量矩阵U的转置，但没有U各特征向量的正交性限制。自编码网络可以通过神经网络的训练过程优化各特征向量，打破正交约束，大大提高数据特征的稀疏性。

然后，网络通过将来自隐层的低维空间特征进行解压缩达到数据重构和还原的目的，这一过程也被称为解码。解码器函数可以表示为

该过程即为式(6)的逆过程，利用编码器的稀疏系数重构数据。式中：Y为神经网络的输出；w´为网络输出的权重；b´为网络输出的偏置。在整个过程中由于压缩作用，会将不相干的噪声压制，而只保留相干性较强的有效信号。

通过上述分析可见，自编码网络工作原理与奇异值分解相似，这为自编码网络提供了工作机理，克服了神经网络的黑箱特性。图1解释了自编码网络与奇异值分解在结构和原理上的对应关系。

图1 自编码网络与奇异值分解原理映射图

基于上述自编码网络去噪原理，设计了无监督网络模型(图2)。首先，应用快速傅里叶变换(FFT)将三维数据变换到频率域，由于频率域数据是复数，将某一切片数据实部和虚部串接作为神经网络的输入。然后，采用自编码网络实现噪声压制，将处理过的所有频率切片重构回三维数据频谱，进行快速傅里叶逆变换(IFFT)得到去噪后数据。

图2 无监督去噪模型框架

在该处理流程中，自编码网络结构是影响去噪的关键因素。一是网络不同层之间的连接方式。本文在隐层添加跳跃连接块以增强网络对地震数据特征的学习，如此，解码层的输入包含两部分，即来自隐层的输出和来自跳跃连接块的输出。跳跃连接块的使用可以加速无监督模型收敛，增强网络提取高阶特征的能力。解码器重建抽象特征，然后将这些抽象特征映射回原始地震数据空间，以实现随机噪声衰减。通过迭代学习抽象特征表示，尽可能减少编码器和解码器的差异。

二是网络中间层的维度决定了数据压缩程度。维度太小，将导致信号压缩过度，损失有效信号；维度太大，压缩不到位，噪声去除效果不佳。采用深度学习的目的是希望获得比SVD 方法更好的压缩编码效果，因此利用SVD 分解后特征值能量占比为90%时对应的阶数，作为中间层维度的参考。此策略保证自编码网络的稀疏度低于SVD 方法，能学习到更稀疏的数据特征。

1.3 基于K-L 散度的优化模型

均方误差函数是回归问题中较为常见的损失函数。根据Bayes 反演理论，最小化均方误差函数是以零均值高斯白噪声为假设条件。但是由于本文方法是在频率域完成，白噪声的存在导致频率切片的均值偏移，因此对去噪结果进行均值约束，有利于数据归位于其真值。K-L 散度可以度量两个分布之间的距离，定义为

式中：Pde=mean(Y)，表示对去噪后数据Y取均值；Pref为参考均值，理论上应该等于干净信号的均值。式（8）度量了去噪数据均值与参考均值的距离，当Pde=Pref时，该值达到最小。

将上述K-L 散度约束与均方误差函数相结合，定义新的损失函数为

正则化项的加入可以驱动去噪结果均值接近参考值。

2 实验结果分析

使用一组合成数据和两组海洋实际数据验证本文方法的有效性，并根据信噪比(Signal-Noise Ratio，SNR)、局部相似性对方法进行评价。

区域(包括本矿区)开展过1∶20万、1∶5万水系沉积物地球化学测量以及1∶1万土壤地球化学测量工作。1∶5万水系沉积物测量成果显示，本区元素异常以Mo、W为主，二者强度高、规模大,且与1∶20万水系沉积物地球化学钼异常总体形态相似，套合较好。1∶1万土壤地球化学测量Mo、W、Ag、Cu异常强度高、规模较大，是矿区的主要成矿元素和伴生元素。Mo异常集中分布在上棋盘北侧和南侧，以及大尖山一带。上棋盘北侧和南侧的Mo异常，总体展布方向呈北北西—南南东向，大尖山一带的Mo异常，总体展布方向呈北西—南东向，与区域上Mo异常成矿带基本一致。经普查地质工作，化探钼异常多为矿致异常。

SNR作为全局衡量标准，定义为

2.1 合成数据

使用褶积模型合成三维数据体。该模型共128×128 道，每道256 个采样点，时间采样间隔为2 ms。对合成数据添加方差为0.3 的高斯噪声，信噪比低至-5.96 dB(图3a)，分别采用MSSA、K-SVD 和本文方法去噪。K-SVD 方法去噪效果最差，信噪比仅提升至2.72dB，同相轴附近残留噪声较多，残差道集也可以看出同相轴附近痕迹明显，说明该方法去除噪声并不彻底(图3b)；MSSA 方法较K-SVD 方法去噪效果有所改善，信噪比提高至4.52 dB，但残差道集有同相轴残留(图3c)，说明该方法对有效信号造成了一定程度的损伤；本文方法去噪较为彻底、效果最好，信噪比提升至6.55 dB，同相轴刻画较为清晰，残差数据也基本看不到同相轴痕迹(图3d)，说明本文方法在压制随机噪声的同时不会对有效信号造成损伤。

图3 合成数据三种方法去噪效果三维对比

为方便观察各方法的数据处理细节，抽取该三维数据Inline 方向第11 个剖面进行对比(图4)。与原始合成数据(图4a 左)相比，加噪数据的同相轴被淹没在强噪声中，信噪比低至-5.98 dB(图4a 右)。K-SVD方法去噪结果的信噪比仅提高到2.47 dB(图4b 左)，残差剖面中同相轴痕迹明显(图4b右)。MSSA 方法去噪结果的信噪比达到4.37 dB，但同相轴处仍有明显的噪声残留(图4c 左)，残差剖面中也留有同相轴痕迹(图4c右)。本文方法去噪后，同相轴更连续，信噪比提高到6.37 dB(图4d 左)，残差剖面中同相轴痕迹不明显(图4d 右)。由此可以看出，本文方法压制噪声效果优于其他两种方法。

图4 合成数据三种方法去噪效果剖面对比

进一步抽取主频为30 Hz 的频率剖面，测试重构阶数或中间层维度对SVD 算法和无监督网络去噪结果的影响，即比较编码稀疏度对去噪结果的影响(图5)。当中间层维度取19 时，无监督网络去噪结果的SNR 最高、效果最好(图5中绿色虚直线②所示)，恰是该频率切片奇异值能量占比为90%所对应的维度。从图中不难发现，无监督网络去噪结果的SNR 一直高于SVD 方法，即当两者有相同的稀疏度时，无监督网络性能更优(图5 中绿色虚直线①所示)。即使在更小的稀疏参数下，无监督网络也可以获得比SVD 方法更佳的去噪效果。

图5 中间层维度对去噪性能的影响

提取出无监督网络第二编码层的特征向量，并计算其协方差矩阵，即各特征向量之间的相关谱(图6)，除了主对角线上，其余位置仍然有不少非零值，可见其特征向量并未完全正交，即采用无监督网络打破了传统SVD 方法特征向量正交这一约束，有利于提高数据表征的稀疏性。

图6 第二编码层特征向量相关谱

去噪方法的时间成本同样是一个值得关注的问题。对于合成三维数据体，MSSA、K-SVD 和本文方法处理时间分别为792、1358、276 s。由于K-SVD 方法每次迭代都需要奇异值分解运算，耗时最大；MSSA 也需要大规模矩阵的奇异值分解，耗时也较大；而本文方法采用的数据结构简单，自编码网络通过反向传播优化网络，计算简单，可以快速实现。

2.2 实际数据一

为了定量评价本文方法的实际数据处理能力，对实际海洋地震数据进行去噪实验。该数据共有400炮，每一炮包含400道数据，每道数据有800个时间采样点，采样间隔为2 ms。图7a是原始地震道集及其加入高斯白噪声的数据，信噪比降至-7.30 dB，原始数据中弱信号同相轴被噪声淹没。图7b～图7d 分别是K-SVD、MSSA 和本文方法的处理结果及其残差道集，可以看出本文方法的残差道集中有效信号能量最弱。对图7 红框区域进行放大显示，如图8所示。可以明显看出，MSSA 和K-SVD 方法去噪后的数据同相轴附近仍残留大量噪声，有效信号能量较弱，而本文方法去噪后同相轴最清晰。K-SVD、MSSA 和本文方法处理结果的SNR 分别为2.76、4.73 和5.91 dB，显然本文方法的去噪效果最好。

图7 实际海洋数据三种方法去噪效果对比

图8 图7 红框区域的放大显示

为了验证本文方法保持数据横向连续性的优势，抽取第10 个共炮检距道集0.4～1.0 s 部分(图9)。由图9 可以看出，本文方法在强噪声背景下仍较好地恢复了同相轴连续性，而且最清晰，尤其是在图中箭头处，效果明显优于其他两种方法。

图9 实际海洋数据三种方法去噪结果的共炮检距道集对比

对于实际数据一，MSSA、K-SVD 和本文方法处理时间分别为6524、1315、1160s。由于实际数据维度较大，因此MSSA 的计算效率低于K-SVD。可见，无论是合成数据还是实际数据，本文方法的计算效率最高。

本文方法中参数Pref的取值会影响各个子网络对频率切片的去噪性能。选取不同的Pref值对30 Hz含噪频率切片数据进行去噪处理并计算其SNR(图10)。含噪切片数据的SNR 为-6.59dB，当Pref取值0.03 时，去噪后SNR 最高，这也恰是无噪数据的均值。随着Pref逐渐偏离真值，去噪后SNR 也会随之下降。由此可以看出，K-L 散度约束对去噪结果有一定影响。

图10 参数Pref对去噪效果的影响

2.3 实际数据二

为了评价本文方法的实用性，对实际含噪数据二进行处理。实际数据一中的随机噪声是人为添加，而实际数据二中是真实噪声。图11a 是原始含噪道集，噪声导致同相轴的清晰度降低；图11b 为K-SVD 方法去噪结果及残差道集，可以看出，噪声去除并不彻底，尤其是远炮检距仍然可以看到大量噪声存在；图11c 为MSSA 方法去噪结果及残差剖面，噪声得到了有效压制，但是同相轴的连续性并未得到恢复，并且残差道集中有明显的有效信号痕迹；图11d是本文方法去噪结果及残差剖面，噪声压制比较彻底、同相轴更清晰，残差道集中有效信号不明显。需承认的是，本文方法在去噪的同时对弱信号会造成轻微的损伤，如图中4.8 s 以后弱信号难以恢复。这是由于在频率域，浅层强信号和深层弱信号耦合，经过自编码网络的压缩会压制弱信号。

图11 实际地震数据二三种方法去噪效果对比

去噪结果与去除的噪声之间的局部相似性能评价信号泄漏程度[27]，高相似异常就表明去噪结果在相应的位置有信号泄漏。图12 展示了K-SVD、MSSA和本文方法去噪结果与去除的噪声之间的相似性。由图可见：K-SVD 方法处理后在2.4～3.2 s 大炮检距处，残差信号与去噪结果存在较高的相似性，说明去噪结果中仍然保留了部分噪声，压制不彻底；经过MSSA 方法去噪后的数据在1.8～3.4 s除了大炮检距处噪声残留，对有效信号损伤也较大，这在图11c的残差道集上也很明显；而本文方法去除的噪声与保留的信号相干性较小，说明在去噪时没有损伤有效信号，具有较高的保真度。

图12 三种方法去噪结果与残差的局部相似性对比

3 结论

本文利用三维数据频域相似性，建立了一种奇异值分解理论指导的自编码无监督去噪网络模型。通过将奇异值分解映射为自编码网络，克服了奇异值分解各特征向量正交的约束，提高了数据表示的稀疏性。根据频率域噪声分布特点在损失函数中引入K-L散度约束，驱动数据恢复其均值分布。本文方法实现了一种运算速度快、又保留传统方法物理可解释性的无监督学习算法模型。合成和实际数据处理结果均表明，本文方法在去噪的同时能够保留信号细节，具有较高的保真度。