基于自监督网络学习与自适应中值滤波的地震数据去噪方法研究

Research on Seismic Data Denoising Method Based on Self Supervised Network Learning and Adaptive Median Filtering

摘要: 在地震图像采集过程中，由于存在各种影响因素，采集到的地震资料常常夹杂着各种噪声，严重影响了后续研究的准确性。因此，图像去噪成为至关重要的一步。自适应中值滤波作为一种经典的去噪方法，在去除椒盐噪声方面表现优异。近来，一种结合传统全变差与自监督学习网络模型的方法——S2S-WTV被提出，该方法不仅拥有神经网络的强大表示能力，还有传统模型的泛化能力，在针对高斯噪声去除方面展现出了出色的性能。基于此，本文提出了一种结合传统方法与神经网络的多种噪声去噪模型，显著提升了去噪后的视觉效果与图像结构相似性。

Abstract: During seismic image acquisition, various factors often result in the data being noisy, severely impacting the accuracy of subsequent studies. Therefore, image denoising is critical. Adaptive median filtering, a classic denoising technique, excels at removing salt and pepper noise. Recently, S2S-WTV, a method that integrates traditional total variation with self-supervised learning network models, has been introduced. This approach combines the strong representational capabilities of neural networks with the generalization abilities of traditional models, showing outstanding performance in reducing Gaussian noise. Building on this, the paper introduces a multi-noise denoising model that combines traditional methods with neural networks, significantly enhancing visual quality and image structure similarity post-denoising.

1. 引言

现如今，地震资料在军事，建筑、遥感以及农业等与我们息息相关的各个领域中都有极其重要的作用。然而，由于地震资料在采集的过程中，受环境、设备以及人工等影响，导致采集到的图像中往往含有一定程度的噪声。为了后续能充分利用地震资料的有效信息，去噪这一步骤就显得尤为重要。

在图像去噪领域，相关研究一直是学术界关注的热点。随着人工智能技术的不断进步，去噪方法逐渐分为两大类：一类是基于传统模型的改进方法，另一类是基于深度学习的神经网络方法。传统去噪方法包括滤波器[1] [2]——通过设计合适的滤波器，将图像中的噪声成分去除只保留有效信息；频率域滤波[3] [4]——将图像转换到频率域，去除噪声频率后再将图像恢复；小波变换[5]——将图像在特定尺度上进行分解，将图像的不同频率分离开方便对噪声进行去除；曲波变换[6] [7]——作为改进版的小波变换，可以更好地捕捉图像的几何特征；模态分解[8]——将图像信号分解为多个固有的模态函数以此来分离信号等。经典算法如BM3D [9]、WNNM [10]以及张量分解方法[11]在去噪效果上取得了不错的成果。然而，这些方法在复杂噪声场景中的表现仍有局限。

与传统方法相比，深度学习技术在图像去噪领域的发展更加迅速。其中，卷积神经网络(CNN) [12]及其改进网络[13]是研究的重点。CNN通过在含噪图像中提取特征，与干净图像进行损失函数计算，经过不断优化迭代，实现噪声的去除。随后，残差网络[14]、FFDnet [15]和Adnet [16]等进一步提升了去噪效果。然而，深度学习方法也面临挑战，尤其是地震数据中，获取足够的训练数据对较为困难。因此，一些研究提出了通过同一图像的不同含噪版本进行训练的网络[17]，如S2S网络[18]，通过伯努利采样生成多个训练数据，以提升泛化能力。尽管深度学习方法在图像去噪效果上表现优异，但其泛化能力较弱，无法完全取代传统方法。结合自监督学习与加权全变差的S2S-WTV网络在此方面提供了一个折中方案[19]，既保留了深度学习的表达能力，也增强了泛化性能。

然而我们知道，野外采集到的地震数据常常含有不同类型的噪声，其中最常见的就是高斯噪声与椒盐噪声。虽然S2S-WTV方法在去除高斯噪声方面效果显著，但对于椒盐噪声的去除效果有限。在去除椒盐噪声方面，作为传统方法的均值滤波有较好的效果，但是会导致图像过度平滑，损失掉一定的细节，影响后续实验。相较而言自适应中值滤波能达到较好的效果。因此，将上述两种方法进行结合可以有效去除地震资料中的不同噪声。最后通过实验验证我们的方法能获得较好的去噪效果。

2. 相关研究

1) 噪声

高斯噪声是一种加性噪声，其概率密度服从高斯分布。而椒盐噪声是一种逻辑噪声，又被称为脉冲噪声，该噪声在图像上显示的情况是在图像上随机出现黑白点。

2) 自适应中值滤波

传统的中值滤波方法是在一定的窗口范围内，将所有的像素点进行排列，选出其中的中值来作为这个范围内中心点的值。这个方式对于椒盐噪声尤为有效，因为椒盐噪声的像素值要不就是极大或者极小，选取中值来替代中心点的值可以有效去除噪声点，保留有效信息。但是，这个方法受窗口大小影响较大，如果窗口较大，确实能获得较好的去噪效果，但是同时会一定程度地损失一部分有效信息，如一些精细的纹理等。但是如果窗口选取过小，则会导致图像的去噪效果欠佳。

因此，在上述基础上，一种自适应中值滤波方法被提出，该方法可以动态地选择滤波窗口的大小，在去除噪声和保留细节两个方面达到了较好的平衡。其原理为：首先假设噪声图像中某一点的像素值为$x_{i,j}$ ，当前选取的窗口为G，窗口中像素最大值为$g_{max}$ ，像素中值为$g_{med}$ ，像素最小值为$g_{min}$ ，M表示窗口的最大尺寸，则有：

(1) 如果$g_{min}<g_{med}<g_{max}$ ，则转到步骤(2)，否则就增大窗口尺寸，如果窗口尺寸 ≤ M，则重复步骤(1)，否则就输出$g_{med}$ 。

(2) 如果$g_{min}<x_{i,j}<g_{max}$ ，则输出$x_{i,j}$ ，否则就输出$g_{med}$ 。

3) S2S-WTV

S2S-WTV方法的优势之一就在于可以通过一张噪声图片恢复出干净的图片，其训练数据集是通过将观测到的地震Y进行伯努利采样来获取的：

${\widehat{Y}}_n=M_n\odot Y,n=1,2,\cdots ,N$

由于地震信号相邻列之间相关性较强，因此我们对数据的一些伯努利采样列施加掩码$M_n$ (伯努利矩阵)，然后使用未掩码的列来预测被掩码的列。

其网络结构是基于U-Net网络结构，包含编码阶段、解码阶段和跳过连接。在编码阶段利用了掩膜引导的残差卷积(MGRConv)而不是部分卷积来对当前卷积层的输出进行加权，它可以在随机掩码下更好地提取有用信息。在每个解码块中使用dropout机制来防止过拟合问题和方差暴增问题。

为了利用地震资料列相关的特性，S2S-WTV还引入了加权TV来保存去噪图像的光滑性，增强其泛化能力。得到的优化方程如下：

$\arg \min {\displaystyle \sum }_{n-1}^N\left({\Vert \left(Y-{\widehat{X}}_n\right)\odot \left(1-M_n\right)\Vert }_F^2+\gamma {\Vert {\widehat{X}}_n\Vert }_{TV,\omega }\right),\text{where}{\widehat{X}}_n=f_{\theta }\left({\widehat{Y}}_n\right)$

其中${\widehat{X}}_n$ 表示干净的信号，$\odot$ 表示逐元素相乘，$M_n$ 表示用于生成样本的掩码，$\gamma$ 表示权衡参数，$f_{\theta }(\cdot )$ 表示具有可学习参数$\theta$ 的卷积神经网络。在上述优化方程中，对于每个实例${\widehat{Y}}_n$ 保真度部分仅在$M_n$ 遮挡的元素上计算，但是由于有许多具有不同掩码的实例，所以可以确保在训练过程中包含Y的所有元素。

最后，通过解码阶段的dropout将训练好的CNN生成多个CNN，这些生成的CNN可能具有一定程度的独立性，因此可以通过对所有输出进行平均来降低去噪结果的噪声。

3. 本文方法

首先，对于含噪的地震图像，我们可以先将其表示为Y，由于原始图像中可能包含不止一种噪声，因此，我们考虑首先去除椒盐噪声。根据上一节中的内容，我们先将原始的地震资料进行自适应中值滤波，以此来去除椒盐噪声，而后对于去噪后的地震资料，我们将其表示为$Y_{med}$ 。对于$Y_{med}$ ，我们只考虑里面的高斯噪声，由于高斯噪声是一种加性噪声，所以我们可以将地震资料表示为如下模型：

$Y_{med}=X+N$

其中X表示干净图像，N表示高斯噪声。而去噪就是在已知$Y_{med}$ 的情况下还原出干净图像X。我们常用的方法是建立正则化模型，通过在模型中引入保真度项来保证噪声图像与干净图像之间偏差不会太大，而正则化项如低秩正则化等，则是通过对数据的某些先验知识来进行约束，以达到更好的去噪效果。通常情况下，正则化模型可以表示为如下形式：

$\min {\Vert Y-X\Vert }_2^2+F\left(X\right)$

其中加号的前部分表示保真度项，而$F\left(X\right)$ 则表示正则化项。

对于本文而言，我们采用S2S-WTV的方法，首先对原始图像进行伯努利采样，获得训练图像，并用其对网络进行训练，得到训练好的网络。将网络训练部分作为正则化项的一部分引入到公式中，得到的模型如下：

$\min {\Vert \left(Y-X\right)\odot \left(1-M_n\right)\Vert }_2^2\text{s}\text{.t}\text{.}X=f_{\theta }\left({\widehat{Y}}_{med,n}\right)$

其中$\odot$ 表示逐元素相乘，$M_n$ 表示用于生成样本的掩码，$f_{\theta }(\cdot )$ 表示具有可学习参数$\theta$ 的S2S网络，${\widehat{Y}}_{med,n}$ 表示原始图像中值滤波后再进行了伯努利采用得到的数据。

在上述模型的基础上，我们引入加权全变差来作为正则化项，以约束地震图像的平滑性，由于地震图像之间相邻列更相似，所以我们只考虑水平方向：

${\Vert X\Vert }_{WTV}={\Vert W\odot \nabla X\Vert }_1$

即梯度仅在水平方向上进行计算：

$\nabla X=X_{i,j+1}-X_{i,j}$

于是我们可以得到我们的模型如下：

$\arg \min {\Vert \left(Y-f_{\theta }\left({\widehat{Y}}_{med,n}\right)\right)\odot \left(1-M_n\right)\Vert }_2^2+\lambda {\Vert {\widehat{Y}}_{med,n}\Vert }_{WTV}$

其中$\lambda$ 表示权衡参数。对于这个模型，我们采用ADMM算法对其进行求解，首先我们在模型中引入新变量，然后通过拉格朗日乘子法将引入变量之后的模型转化为无约束优化的形式：

$\begin{array}{c}L\left(\theta ,P,Q\right)=\arg \min {\displaystyle \sum }_{n=1}^N{\Vert \left(Y-f_{\theta }\left({\widehat{Y}}_{med,n}\right)\right)\odot \left(1-M_n\right)\Vert }_2^2+\lambda {\Vert P_n\Vert }_{WTV}\\ +\frac{\mu }{2}{\Vert P_n-\nabla f_{\theta }\left({\widehat{Y}}_{med,n}\right)-Q_n\Vert }_F^2\end{array}$

其中$\mu$ 是惩罚参数，而之所以要求和是由于在网络训练之前我们通过伯努利采样生成了N个数据，然而在实际计算中，由于同时计算N个数据计算量过大，因此在每次迭代时我们将随机选取一个数据，这里我们假设在第$k+1$ 次迭代时选取第m个数据。对于上述无约束优化式子，我们通过ADMM将其转化为多个子问题进行求解，依次对几个变量进行不断轮换的迭代优化。

求解关于变量$\theta$ 的子问题：

对于第k次迭代，我们使用自适应矩估计算法(Adam)来对CNN的参数$\theta$ 进行迭代，即在ADMM每次迭代的过程中使用一步Adam，并将以下式子作为损失函数来获得更新后的参数${\theta }^{k+1}$ 。

$\underset{\theta }{\arg \min }{\Vert \left(Y-f_{\theta }\left({\widehat{Y}}_{med,m}\right)\right)\odot \left(1-M_m\right)\Vert }_2^2+\frac{\mu }{2}{\Vert P_m-\nabla f_{\theta }\left({\widehat{Y}}_{med,m}\right)-Q_m\Vert }_F^2$

求解关于变量P的子问题：

$P_m^{k+1}=\arg \underset{P}{\min }\lambda {\Vert P_m\Vert }_{WTV}+\frac{\mu }{2}{\Vert P_m-\nabla f_{{\theta }^{k+1}}\left({\widehat{Y}}_{med,m}\right)-Q_m\Vert }_F^2$

对于这个子问题，我们采用软阈值法进行求解。

求解关于变量Q的子问题：

$Q_m^{k+1}=Q_m^k+\mu \left(\nabla f_{{\theta }^k}\left({\widehat{Y}}_{med,m}\right)-Q_m^k\right)$

对于权重矩阵而言，每一百次更新一次，更新公式如下所示：

${\left(W^{k+1}\right)}_{\left(i,j\right)}=\{\begin{array}{l}\frac{{\Vert Y-f_{{\theta }^k}\left({\widehat{Y}}_m\right)\Vert }_F^2}{2HW\left|{\left({\nabla }_h{f}_{{\theta }^k}\left({\widehat{Y}}_m\right)\right)}_{\left(i,j\right)}\right|},k\mathrm{mod}100=0\\ {\left(W^k\right)}_{\left(i,j\right)},k\mathrm{mod}100\ne 0\end{array}$

对于上述迭代公式，当迭代次数达到我们设置的最大次数时则停止迭代，在本文中迭代次数设置为5000次。迭代完成之后我们将得到训练好的神经网络$f_{\theta }$ ，然后通过解码阶段的dropout生成多个神经网络，对这些网络得到的去噪结果进行平均则能得到我们最终的去噪结果。

4. 实验与分析

为了验证所提出方法的有效性，我们选取一张合成资料对其进行加噪处理，然后对噪声图像采用我们提出的方法进行去噪实验，随后通过定量和定性两种方式进行分析。除了合成资料以外，由于椒盐噪声会在运输与储存的过程中产生，我们还进一步选用了有椒盐噪声的实际地震资料来验证所提方法的去噪效果。对于上述两个资料，我们均截取大小为256 × 256的图像，并将图像全部转化为灰度图像，方便后续处理。以下为我们选取的干净与添加了一定噪声(Gaussian/Salt Pepper = 30/0.05)的合成资料以及加了椒盐噪声的实际资料(图1)。

Figure 1. Original data graph

图1. 原始数据图

除了用本文方法进行实验以外，我们还将本文方法与一些经典方法，如自适应中值滤波、BM3D、S2S-WTV进行对比验证。实验结果将在后续进行展示。

以表1为对合成资料添加不同程度的高斯和椒盐噪声进行去噪得到的峰值信噪比(PSNR)以及结构相似性(SSIM)对比。通过观察表格中的数据可以发现，较好的数值结果主要是出现在本文所提出的方法中，尤其是在结构相似度方面本文方法有显著提升，说明我们的方法在去除地震资料噪声方面有较为优秀的效果。而自适应中值滤波方法由于其对椒盐噪声有较好的去除效果，但对高斯噪声效果欠佳，所以其整体去噪效果表现较差。而BM3D主要是用于去除高斯噪声，所以对含有椒盐噪声的图像处理效果较差。

Table 1. Comparative analysis

表1. 对比分析

Noise (Gaussian/Salt Pepper)	$\sigma =30/0.05$		$\sigma =50/0.07$		$\sigma =70/0.09$
method	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
observd	14.88	0.1382	12.92	0.0857	11.65	0.0603
自适应中值滤波	20.69	0.3411	15.59	0.2268	14.00	0.1605
BM3D	21.38	0.6195	17.48	0.3073	16.18	0.2296
S2S-WTV	23.81	08469	18.12	0.6256	15.73	0.4737
mine	24.01	0.8609	18.93	0.7549	15.84	0.6432

除了上述定量分析外，我们还对五种方法去噪后的结果进行定性分析，通过观察不同方法去噪后得到的结果图我们可以发现，自适应中值滤波对于椒盐噪声去噪效果较好。BM3D对于高斯噪声效果较好，去噪后的图像还能看到明显椒盐噪声。S2S-WTV方法虽然数值效果较好，但是视觉效果较差。而本文方法处理后的图像可以看到不管是高斯噪声还是椒盐噪声都有较好的去除效果，恢复出的图像同相轴也较为清晰，表现出的视觉效果较好(图2)。

Figure 2. From left to right, the denoising images of adaptive median filtering, BM3D, S2S-WTV and the method in this paper

图2. 从左往右依次为自适应中值滤波、BM3D、S2S-WTV以及本文方法去噪图

除了上述合成资料以外，我们还对实际的地震资料进行实验。由于实际资料我们无法获得其干净的图像，所以难以通过定量实验来分析去噪效果的优劣，但是我们同样可以从定性分析入手，观察不同方法去噪后得到的效果，同样还可以观察残差图来辅助判断去噪效果。以下为几种方法的去噪结果图：我们可以发现自适应中值滤波去掉了一部分有效信号和椒盐噪声，但是大部分噪声未被去除。BM3D去除了大部分椒盐噪声但是还有些不规则噪声未去除，并且也有一定程度的信号泄漏。S2S-WTV虽然没有去除有效信号，但对噪声的去除有限。本文方法虽然有一定程度的噪声泄漏，但是去除了图像中的绝大部分噪声(图3)。

5. 结论

本文提出了一种结合自适应中值滤波和S2S-WTV的去噪方法，旨在去除地震图像中的高斯噪声

Figure 3. The first row from left to right shows the denoising images of adaptive median filtering, BM3D, S2S-WTV and the method proposed in this paper; the second row from left to right shows the residual images of adaptive median filtering, BM3D, S2S-WTV and the method proposed in this paper

图3. 第一行从左往右依次为自适应中值滤波、BM3D、S2S-WTV以及本文方法去噪图；第二行从左往右依次为自适应中值滤波、BM3D、S2S-WTV以及本文方法的残差图

与椒盐噪声的同时保留更多的细节信息，以便后续的处理。通过合成资料与实际资料实验验证得知本文方法在去除高斯噪声与椒盐噪声方面取得了良好的效果，与几种经典的方法进行对比，在PSNR、SSIM尤其是SSIM数值方面都有一定的提升，其视觉效果相较而言也较为优秀。