基于低秩近似与自监督学习的地震数据去噪方法研究

Research on Seismic Data Denoising Method Based on Low-Rank Approximation and Self-Supervised Learning

摘要: 地震图像由不同频段的图像组成，提供丰富的地层和纹理信息。然而，地震图像常常受到环境、设备因素干扰，伴随着复杂的随机噪声。传统的地震图像去噪方法依赖约束优化方法，对于先验知识的选择也十分重要。由于受限于手工先验，传统方法的去噪性能急需提升。近来，深度学习被广泛应用于图像处理领域，但是需要训练大量的成对样本，对未经过训练的噪声分布泛化能力差。因此，一种基于优化的方法和深度学习的深度低秩分解的自监督网络的方法——Flex-DLD被提出，通过用一个可训练的神经网络模块(DLD)替代多个复杂的手工求解器，并将其他先验转化为损失函数项，极大地简化了算法的设计和实现流程。本文成功地将基于模型的优化方法与数据驱动的深度学习方法(深度低秩分解)相结合。既充分利用了神经网络强大的非线性表示能力替代传统的低秩分解(如SVD)，又通过WTV来引导网络学习，以更好在去噪的同时保护地震剖面中的结构和边缘信息。

Abstract: Seismic images consist of multi-frequency bands, providing rich stratigraphic and textural information. However, seismic images are often contaminated by environmental and equipment factors, accompanied by complex random noise. Traditional seismic image denoising methods rely on constrained optimization approaches, where the selection of prior knowledge is also crucial. Due to limitations imposed by manual priors, the denoising performance of traditional methods urgently requires improvement. Recently, deep learning has been widely applied in image processing, but it requires training on large amounts of paired samples and exhibits poor generalization ability for untrained noise distributions. Therefore, Flex-DLD—a method combining optimization-based approaches with deep low-rank decomposition via self-supervised networks—has been proposed. By replacing multiple complex manual solvers with a single trainable neural network module (DLD) and converting other priors into loss function terms, the design and implementation of the algorithm are greatly simplified. This paper successfully combines model-based optimization methods with data-driven deep learning methods (deep low-rank decomposition). It fully utilizes the powerful nonlinear representation capabilities of neural networks to replace traditional low-rank decompositions (such as SVD), and uses WTV to guide network learning, so as to better preserve structural and edge information in seismic profiles while denoising.

1. 引言

在地震资料采集过程中由于环境所产生的随机噪声，也被称为非相干噪声，会影响与信号相关的成像的准确性。所以抑制随机噪声对地震数据处理和解释具有重要意义。地震数据去噪领域的传统方法有基于滤波器、变换、降阶等。基于滤波器的方法将信号与噪声分离，比如均值滤波[1]，[2]、双边滤波[3]等。基于变换的方法将原始数据转换到另一个域，产生稀疏信号，从而可以轻松地从噪声中分离出来，包括傅里叶变换[4]、curvelet、小波变换[5]等。基于降阶的方法是指地震信号经过一定的预变换后达到较低的秩，而噪声提高了秩，从而去噪问题可以看作是矩阵降维的问题。

传统的地震数据去噪方法在信噪比高时十分有效，但不适合信噪比低的地震信号。相比这些传统方法，低秩近似在模型方法中的应用最为广泛。低秩近似是通过找到原始矩阵的低秩矩阵来降低矩阵的维数。而张量分解则是扩展了低秩近似，将高维张量分解为一组低维因子，充分利用其各个维度的信息，从而提高了去噪的效果。比如，常见的Tucker分解[6]和CP分解[7] [8]，Tucker分解为几个矩阵与核心张量[9]的乘积，CP分解可以将输入的张量分解为秩一张量[9]的加权组合形式。低秩近似和张量分解技术可以减少模型大小和计算复杂性，但是基于低秩的方法需要手动调整模型的参数和先验，容易陷入局部最优，使得去噪效果不佳。

相比之下，深度学习可以通过学习数据集来调整参数，具有较强的学习能力。比如卷积神经网络(CNN) [10]能通过结构优化提升特征表达能力，后续出现的VGGNet [11]、DnCNN [12]、ResNet [13]，提升了去噪的效果。虽然深度学习方法较传统方法具有更强的泛化与表达能力，但是常用的监督学习方法需要大量的配对的地震数据，时间和成本花费太大。因此，出现了仅需单幅噪声图像的自监督方法，自监督方法如DIP [14]、S2S [15]和S2S-WTV [16]，无需成对样本，结构轻量、内存占用低，有效提升了去噪性能与泛化能力。

为了充分利用神经网络的强大学习能力，同时解决传统低秩方法仅利用局部信息和配对样本不足的问题，本文利用DLD模型的框架，以自监督学习的方式进行完全训练，不利用SVD技术和现有的深度的降噪器。基于地震去噪的低秩模型来解释模型的网络设置，根据地层的特点选择相应的卷积层，并且将传统的先验知识相结合对地震数据进行去噪。

2. 相关方法

2.1. 低秩近似

低秩近似寻求秩为$R=\min \left(M,N\right)$ 的矩阵$X\in {ℝ}^{M\times N}$ 表示为两个较小的矩阵$U\in {ℝ}^{M\times K}$ 、$V\in {ℝ}^{N\times K}$ 的乘积，其中$K$ 通常被认为小于$R$ 。对$U$ 、$V$ 和所需的目标的特定约束产生了不同类型的低秩近似算法。低秩近似常用于去噪、压缩，常用的有奇异值分解(SVD)，可以表示为

$\underset{U,V}{\min }{\Vert X-U{V}^{\text{T}}\Vert }_F\text{s}\text{.t}\text{.}U=V^{\text{T}}$

2.2. 通道注意力

后续的读者对SE模块进行改进，提出通道注意力(Channel Attention Mechanism, CAM) [17]通过显式地建模卷积特征通道之间的相互依赖关系来改善网络产生的表示的质量。CAM通过分析通道间的关系，动态调整各通道权重，强调重要特征通道，抑制无关通道。联合使用两种池化，分别利用平均池化捕获全局统计特征最大池化提取局部显著特征，通过共享MLP生成通道权重，表示如下：

$\begin{array}{c}{M}_c\left(F\right)=\sigma \left(\text{MLP}\left(\text{AvgPool}\left(F\right)\right)+\text{MLP}\left(\text{MaxPool}\left(F\right)\right)\right)\\ =\sigma \left(W_1\left(W_0\left(F_{\text{avg}}^c\right)\right)+W_1\left(W_0\left(F_{\max }^c\right)\right)\right)\end{array}$

其中，$\sigma$ 代表Sigmoid激活函数，$W_0\in {ℝ}^{C/r\times C}$ 、$W_1\in {ℝ}^{C/r\times C}$ 是MLP的权重，$M_C\in {ℝ}^{C\times 1\times 1}$ 是通道注意力权重。$F_{\text{avg}}^C$ 、$F_{\max }^C$ 分别通道注意力阶段表示经过全局平均池化和最大池化后的特征。

2.3. Flex-DLD去噪框架

Flex-DLD [18]将传统先验灵活结合深度低秩模型，其优势在于以自监督学习的方式进行完全训练，没有利用SVD技术和现有的深度降噪器，能够将HIS的空间和光谱平滑度等传统先验知识纳入优化框架，并且解决了DHP [19]/DeepTensor [20]面临的挑战。同时还解决了传统模型需要训练大量样本的问题。

对于HIS去噪问题，可以表示为

$Y=X+N+S$

$Y=\left[y_1,y_2,\cdots ,y_B\right]$

其中，$Y\in {ℝ}^{N\times B}$ ($N=H\times W$ 表示图像大小)是$\mathcal{Y}$ 的向量化表示，$N\in {ℝ}^{N\times B}$ 表示高斯噪声，$S\in {ℝ}^{N\times B}$ 表示非高斯(稀疏)噪声。$y_i$ 表示第$i$ 个波段图像的矢量化形式。

去噪的优化问题为

$X=\underset{X}{\arg \min }\frac{1}{2}{\Vert Y-X-S\Vert }_F^2+\lambda g\left(X\right)$

低秩分解的目的是将矩阵$Y$ 近似表示为两个低秩矩阵的乘积

$Y\approx L=U{V}^{\text{T}}$

其中，$L\in {ℝ}^{N\times B}$ 是$Y$ 的低秩表示，$U\in {ℝ}^{N\times P}$ 和$V\in {ℝ}^{S\times B}$ 可以看作奇异值分解的左奇异值向量和右奇异值向量，且$P\ll \min \left(N,B\right)$ 。所以$L$ 的秩远低于$Y$ 。

该框架引入一个即插即用的深度低秩分解模型。通过两个生成式神经网络$f_U(\cdot )$ 和$f_V(\cdot )$ 来获取$U$ 和$V$ 。$U$ 是网络$f_U(\cdot )$ 的输出，为了更好地捕获空间信息，$f_U(\cdot )$ 输入和输出是形状为$H\times W\times P$ 的三维张量。$V$ 是网络$f_V(\cdot )$ 的输出，捕获频谱上的相关性。$f_U(\cdot )$ 按照编码器–解码器的框架，且其中的编码器和解码器都包含$S_U$ 。每个阶段包含两个二维卷积块和两个通道注意力模块。卷积块由2D卷积层、2D BN和ReLU激活函数。

$f_V(\cdot )$ 的编码器和解码器都由$S_V$ 阶段组成。每个编码器阶段包括两个1D卷积层、两个BN层和两个ReLU激活函数。每个解码器阶段由一个上采样、1D卷积层、BN和激活函数组成。生成$U$ 和$V$ 以获得$Y$ 的低秩表示可以通过训练网络来实现。设$W_U$ 和$W_V$ 为$f_U(\cdot )$ 和$f_V(\cdot )$ 的可学习参数。网络参数通过最小化以下目标函数进行联合优化。

${ℒ}_{\text{DLD}}={\Vert Y-U{V}^{\text{T}}\Vert }_F^2$

其中，$U=f_U\left(E_U\right)$ ，$V=f_V\left(E_V\right)$ 。因此，DLD是以自监督方式训练的，无需干净地震数据作为标签。所以，优化问题可以写为

$\underset{W_u,W_V}{\min }{\Vert Y-f_U\left(E_U\right)f_V{\left(E_V\right)}^{\text{T}}\Vert }_F^2$

通过深度学习和基于模型的框架联合作用实现去噪，并且结合传统先验知识TV [21]、SSTV [22]，进一步提高去噪和恢复的性能，并且DLD框架可以与其他轻量网络和先验知识灵活结合。

3. 本文方法

地震资料具有独特的性质，常见的有稀疏性和低秩性，通常通过平滑先验和低秩近似等方法去进行图像去噪。通常采集的地震数据包含随机噪声，复杂而多样，本文只考虑高斯白噪声，利用$T$ 和$S$ 分别表示时间采样点和地震道数。$X\in {ℝ}^{T\times S}$ 表示待恢复的干净地震数据。其退化模型可以表示为

$Y=X+N$

其中，$N\in {ℝ}^{T\times S}$ 表示加性噪声。其目的是从观测到的$Y$ 中还原出干净的$X$ 。

从最大后验概率估计出发，该问题可形式化为

$\widehat{X}=\arg \underset{X}{\max }\log P\left(Y|X\right)+\log P\left(X\right)$

图像去噪问题通常建立正则化模型，包含两项，前一项是数据保真项，确保观测到的图像与还原出来的干净图像不会有太大的误差，后一项是正则化项，用于平滑或者约束图像，提高局部或者全局的去噪性能，并且保留图像的细节信息。

假设噪声N服从独立同分布的高斯分布，则上述问题等价于以下最小化问题

$\widehat{X}=\arg \underset{X}{\min }\frac{1}{2}{\Vert Y-X\Vert }_F^2+\lambda g\left(X\right)$

其中，$\lambda$ 表示平衡因子，$g\left(X\right)$ 为正则化项，表示对干净地震数据的先验知识。

受Flex-DLD启发，将低秩矩阵分解与深度学习结合，提出DLD-WTV的地震数据去噪模型。将加权全变差先验加入DLD深度学习框架。WTV对地震数据施加自适应平滑约束。

$\text{WTV}\left(L\right)={\Vert W_t\odot \left(D_{t}L\right)\Vert }_1+{\Vert W_s\odot \left(D_{s}L\right)\Vert }_1$

其中，$D_t$ 和$D_s$ 分别代表沿时间方向和道方向的差分算子，${\Vert \cdot \Vert }_1$ 表示$L_1$ 范数，$\odot$ 表示逐元素相乘，$W_t$ ，$W_s\in {ℝ}^{T\times S}$ 为权重矩阵。

所以带有WTV先验的DLD的目标函数为

${ℒ}_{\text{DLD-WTV}}={ℒ}_{\text{DLD}}+{\lambda }_R\cdot \text{WTV}\left(L\right)$

其中${\lambda }_R$ 表示平衡因子。优化过程：最小化损失函数来联合训练网络参数$f_U(\cdot )$ 和$f_V(\cdot )$ 。利用随机梯度下降算法，如ADAM作为优化算法。

基于上式，我们可以得到其优化模型

$\widehat{X}=\arg \underset{X}{\min }\frac{1}{2}{\Vert Y-X\Vert }_F^2+\text{DLD-WTV}\left(X\right)$

首先在模型中引入辅助变量$Z$ ，然后通过拉格朗日乘子法将引入变量之后的模型转换为无约束优化的形式。

$\left(\widehat{X},\widehat{Z}\right)=\arg \underset{\left(X,Z\right)}{\min }\frac{1}{2}{\Vert Y-X\Vert }_F^2+{\Vert Z-U{V}^{\text{T}}\Vert }_F^2+{\lambda }_R{\displaystyle \sum _{i,j}{\omega }_{i,j}}\left|\nabla X\left(i,j\right)\right|$

$\text{s}.\text{t}.X=Z$

构造增广拉格朗日函数

${ℒ}_{\rho }\left(X,Z,u\right)=\frac{1}{2}{\Vert Y-X\Vert }_F^2+{\Vert Z-U{V}^{\text{T}}\Vert }_F^2+{\lambda }_R{\displaystyle \sum _{i,j}{\omega }_{i,j}\left|\nabla X\left(i,j\right)\right|}+\frac{\rho }{2}{\Vert X-Z+u\Vert }_F^2$

其中$u$ 为拉格朗日乘子，$\rho >0$ 为惩罚函数。通过最小化${ℒ}_{\rho }$ ，依次更新$X$ 、$Z$ 和$u$ 。

故$X$ 、$Z$ 和$u$ 的最小化问题如下：

$\{\begin{array}{l}\widehat{X}=\arg \underset{X}{\min }\frac{1}{2}{\Vert Y-X\Vert }_F^2+\frac{\rho }{2}{\Vert X-Z+u\Vert }_F^2\\ \widehat{Z}=\arg \underset{Z}{\min }\text{DLD-WTV}\left(Z\right)+\frac{\rho }{2}{\Vert Z-X-u\Vert }_F^2\\ \widehat{u}=u+\left(\widehat{X}-\widehat{Z}\right)\end{array}$

WTV中自适应权重W的更新，在每次迭代$k+1$ 开始时，根据前一轮的干净估计$X^k$ 更新WTV权重W

$W_{i,j}^{k+1}=\frac{1}{{\left|\nabla X^k\right|}_{i,j}+ϵ}$

其中${\left|\nabla X^k\right|}_{i,j}$ 为$X^k$ 在$\left(i,j\right)$ 处的局部梯度模，$ϵ$ 为稳定常数。该权重$W^{k+1}$ 用于后续Z-子问题的计算。

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

$X^{k+1}=\frac{1}{\rho +1}\left[Y+\rho \left(Z^K-u^k\right)\right]$

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

$Z^{k+1}=U^{k+1}{V}^{k+1^{\text{T}}}$

其中，$U^{k+1}$ 和$V^{k+1}$ 分别作为每次迭代时网络$f_U(\cdot )$ 和$f_V(\cdot )$ 的最终输出，所以具有WTV先验的DLD模型，$U^{k+1}$ 和$V^{k+1}$ 通过以下方式获得

$\left(U^{k+1},V^{k+1}\right)=\underset{\left(U^{k+1},V^{k+1}\right)}{\arg \min }{ℒ}_{\text{DLD-WTV}}^{k+1}$

${ℒ}_{\text{DLD-WTV}}^{k+1}={\Vert \left(X^{k+1}+u^k\right)-U^{k+1}{V}^{k+1^{\text{T}}}\Vert }_F^2+{\lambda }_R\text{WTV}\left(U^{k+1}{V}^{k+1^{\text{T}}}\right)$

求解关于$u$ 的子问题：

$u^{k+1}=u^k+\left(X^{k+1}-Z^{k+1}\right)$

4. 实验与分析

为了验证所提方法的有效性，我们选取模拟资料对其进行添加加性高斯白噪声，然后对加噪的图像应用所提出的方法，获得实验结果，并采用定量与定性两种方法进行结果分析。此外，还选取了含有高斯噪声的实际地震图像验证所提出方法的去噪效果。对于以上两种资料，我们截取256 × 256大小的地震图像，灰度值范围为0~255，以便后续实验的开展。选取的干净的模拟地震资料、以及添加一定噪声的模拟资料和实际地震资料图像如图1所示。

Figure 1. Original data graph

图1. 原始数据图

为了评估所提出的方法DLD-WTV的去噪性能，为模拟地震资料分别加上方差为30、50、70的高斯白噪声，再将所提出的方法与BM3D、DLD、DLD-TV对不同的噪声进行实验，验证其去噪性能。由图2所见，所提出的DLD-WTV方法，相对于BM3D、DLD、DLD-TV，去噪后的视觉效果更好。相对于BM3D，DLD-WTV的去噪效果更好，轮廓更加清晰；相比于DLD、DLD-TV，DLD-WTV对边缘和细节的保留更加完整，并且相对于前两者对于平滑和边缘区域赋予相同权重，DLD-WTV更能平衡平坦区域和边缘区域的惩罚强度，有更好的去噪性能和更多的纹理信息。

Figure 2. From left to right, the denoising images of BM3D, DLD, DLD-TV and the method in this paper

图2. 从左往右依次为BM3D、DLD、DLD-TV以及本文方法去噪图

为了更加清楚的展示DLD-WTV的去噪效果，分别计算了去噪后图像的峰值信噪比(PSNR)和结构相似度(SSIM)，如表1所示。通常认为PSNR大于30时有较好的去噪效果，PSNR越大，去噪效果越好。SSIM则是越接近1越好，表示恢复的图像与原始图像越接近。所以可以发现DLD-WTV的去噪效果在PSNR、SSIM上有很大的提升，噪声水平越高时，去噪效果越好，并且有很好的稳定性。

Table 1. Comparative analysis

表1. 对比分析

Noise (Gaussian)	$\sigma =30$		$\sigma =50$		$\sigma =70$
Method	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
BM3D	36.77	0.95	33.52	0.93	31.01	0.92
DLD	37.25	0.96	33.67	0.95	32.89	0.93
DLD-TV	37.88	0.97	34.23	0.95	33.22	0.94
DLD-WTV	38.21	0.98	35.56	0.97	35.05	0.96

除了对模拟地震资料进行实验外，我们还对实际地震资料进行了实验。因为实际地震资料无法获得其干净的图像，所以难以通过定量指标PSNR、SSIM去判断去噪效果的好坏。因此，我们用定性的方法对实际地震资料去噪后的效果进行分析，并且可以通过残差图对其进一步补充分析。通过图3中的去噪结果图，我们发现BM3D和DLD在一定程度上去除了噪声，但是还有部分噪声未被去除，并且观察其残差图像，发现都存在着一定程度的信号泄漏。DLD-TV相比前两种方法，几乎不存在信号泄漏，去噪性能和平滑性能都表现更好，但是去噪效果有限，不能彻底去除噪声。DLD-WTV方法能够有区分性的对平坦区域和边缘进行去噪，从去噪的效果图和残差图来看，具有良好的去噪效果，保留信号的能力也比较出色。

Figure 3. The first row from left to right shows the denoising images of BM3D, DLD, DLD-TV and the method proposed in this paper; the second row from left to right shows the residual images of BM3D, DLD, DLD-TV and the method proposed in this paper

图3. 第一行从左往右依次为BM3D、DLD、DLD-TV以及本文方法的去噪图；第二行从左往右依次为BM3D、DLD、DLD-TV以及本文方法的残差图

5. 结论

考虑到低秩近似方法的局限性与深度学习的强大学习能力，本文提出一种基于优化的方法与深度学习结合的深度低秩分解的自监督网络，通过修改其子网络以适应地震资料去噪，并结合自适应的加权全变差，用一个可训练的神经网络模块(DLD)替代多个复杂的手工求解器，并将其他先验转化为损失函数项，极大地简化了算法的设计和实现流程，提高了模型的可解释性。本文旨在去除地震资料中的高斯噪声的同时保留更多的细节和纹理信息，方便后续的分析。实验结果表明，本文方法除了能够发挥低秩近似与自监督学习的优点，还能有效保留边缘细节，并且在模拟地震图像和实际地震图像上都有较好的去噪效果。但与常用的监督学习方法相比，该方法需要更多的时间来进行去噪。并且需要更强大的GPU硬件进行计算，训练期间的一些参数也是手动筛选。在后续的工作中，将考虑设计一个更加轻量的网络以节省运行的时间成本，并考虑增加自适应机制来自动筛选参数。进一步考虑张量分解与深度学习结合的可能性，通过实验深入验证去噪性能。