在地震勘探中, 为了采集到正确的地震数据, 采样数据必须要满足Nyquist采样定理, 即采样频率必须大于最大频率的两倍。在地震数据的野外采集过程中, 由于建筑物、湖泊和冰川等地形条件的限制, 采集到的地震数据通常在空间上不完整且不规则, 致使数据存在缺失或混叠现象, 对后续处理产生严重的干扰。所以如何用不完整、不规则的数据重建较为完整、规则的数据在实际应用中有巨大作用。

传统的地震数据重建方法主要分为3类, 分别为预测滤波器方法^[1]、波动方程方法^[2-3]及信号分析和统计的方法^[4-6]。压缩感知理论是在基于信号分析和统计基础上提出的一种以低采样率数据获得完整信号的方法, 突破了Nyquist采样定理的限制, 为“两宽一高”地震勘探技术提供了新思路。王常波^[7]在压缩感知的理论框架下, 利用Shearlet稀疏基重构信噪比更高的地震数据。张良等^[8]对比离散余弦稀疏基、傅里叶稀疏基、小波稀疏基、Curvelet稀疏基以及Shearlet稀疏基在合成记录和实际地震记录的应用效果, 指出地震信号在Shearlet变换域内的稀疏程度高于其他稀疏基, 基于Shearlet稀疏基重构出的效果最好。杨冠雨等^[9]在Shearlet变换基础上, 利用L₁范数和广义全变分双正则化方法对连续缺失数据进行更有效地重建。但是, 传统的压缩感知方法需要人为设定阈值, 在实际问题中需要耗费大量的人力和经验, 并且该方法在数据大量缺失的情况下效果不佳。王华忠^[10]指出地震数据插值与地下介质的复杂程度息息相关, 由于地下介质的未知、不可预测性, 直接限制了压缩感知方法的应用, 甚至止步不前。

近年来随着计算机硬件性能的提升以及人工智能技术的发展, 使用人工智能手段处理地震数据已成为新的研究热点。深度学习是一种基于数据驱动的方法。不依赖于物理方程与公式理论, 仅通过大量数据训练, 便可以学习输入数据和标签数据间的映射关系, 并利用这种映射关系解决相似的问题。在地球物理领域中, 深度学习已经在噪声压制、初至波拾取、油气储层分布预测、波形分类与识别等方面取得了进展。MANDELLI等^[11]在随机缺失的地震道中引入CNN实现地震道插值。WANG等^[12]进一步证明CNN能够有效地应用于规则数据和不规则数据重建, 并从时间域、空间域、f-k域与常规方法的结果对比, 指出残差网络在地震数据重建中能够表现出很好的效果。王峰^[13]利用CNN同时对地震数据进行去噪和重建处理, 更高效地得到了高精度的结果。在进行网络训练前, 加入预处理过程使输入数据和输出数据之间的映射简单, 进而得到更好的重建结果。GAO等^[14]采用Shepard插值方法对数据进行预处理以降低输入数据和标签数据之间的不确定性, 通过CNN获得信噪比较高的重建数据。WANG等^[15]对规则数据进行3次样条插值, 随后利用残差网络完成规则地震数据的重建。JIA等^[16]在训练前对数据进行双3次插值, 并利用支持向量机方法进行重建, 重建效果优于传统的f-x插值方法。利用深度学习方法对地震数据重建的过程中, 对数据进行初步处理能够提高输入数据和输出数据的质量, 加快网络的收敛速度, 同时也能获得高精度的结果。

基于模型驱动的压缩感知重建基础上, 结合数据驱动方法, 提出一种基于压缩感知和深度学习的地震数据重建方法。首先在先验知识的约束下, 利用地震信号在Shearlet域内的稀疏性对其进行压缩感知初步重建, 随后利用深度学习方法进行监督学习实现数据重建的复杂迭代, 在训练过程中采用数据块处理的方式构建训练数据以学习地震信号微小特征, 最后利用训练好的网络模型实现不规则地震数据的重建。合成地震数据和实际地震数据重建结果验证了这种方法的有效性。

1 方法原理 1.1 基于压缩感知的数据重建方法

地震数据重构可以归纳为:

$ d_{\mathrm{obs}}= \boldsymbol{Md} $

(1)

式中: d_obs代表不规则的观测地震数据, M代表采样矩阵, d代表需要重建的规则数据。

采样矩阵M是由0, 1元素构成的对角矩阵, 矩阵对角线上的元素与d_obs相对应, 数据缺失时对角线元素为0, 反之, 对角线元素为1^[17]。

压缩感知理论的提出为欠定问题的求解提供了一种新思路^[18], 该理论指出采样数据满足以下3个条件便可以重建出完整的原始信号, 分别为: ①信号d是稀疏的或者是可压缩的; ②采样矩阵是随机的且和信号本身不相干; ③需要某种促进稀疏性的重建算法。

假设原始地震数据是稀疏的或在某种变换域内是稀疏的, 即可以通过某种变换使得数据能够进行稀疏表示。相比于其它的变换, Shearlet变换结合了多尺度的方法^[19], 具有较好的方向性和稀疏性, 能够更好地刻画图像边缘信息。在二维情况下, 连续的Shearlet变换系统可以表示为:

$ \varphi_{\text {ast }}=a^{-\frac{3}{4}} \varphi\left[\boldsymbol{A}_a^{-1} \boldsymbol{B}_b^{-1}(x-t)\right] $

(2)

式中: $\boldsymbol{A}_a=\left[\begin{array}{cc}a & 0 \\ 0 & \sqrt{a}\end{array}\right]$代表各向异性膨胀矩阵; B_b=$\left[\begin{array}{ll}1 & b \\ 0 & 1\end{array}\right]$代表剪切矩阵; φ代表Shearlet母函数; a代表尺度参数; b代表剪切参数。

采用Shearlet变换对数据进行稀疏表示, 该过程可以表示为:

$ \boldsymbol{x}=\boldsymbol{S d} $

(3)

式中: x代表信号d在Shearlet域中的稀疏系数; S代表Shearlet变换。

令矩阵D=MS^－1, 其中S^－1代表Shearlet逆变换, 则(1)式可写为:

$ d_{\mathrm{obs}}=\boldsymbol{Dx} $

(4)

压缩感知理论指出(4)式可以用下列约束进行求解:

$ \left\{\begin{array}{l} \tilde{\boldsymbol{x}}=\operatorname{argmin} _\boldsymbol{x}\|\boldsymbol{x}\|_0 \quad \text { s.t. } \quad \boldsymbol{D} \boldsymbol{x}=\boldsymbol{d}_{\mathrm{obs}} \\ \tilde{\boldsymbol{d}}=\boldsymbol{S}^{-1} \tilde{\boldsymbol{x}} \end{array}\right. $

(5)

式中: ‖·‖₀代表L₀范数, 即表示向量中非0元素的个数。由于0范数的求解存在NP难问题, 在实际应用中, 经常用L₁范数代替L₀范数, 所以(5)式可以改写为下方的无约束优化问题:

$ \left\{\begin{array}{l} \tilde{\boldsymbol{x}}_\lambda=\arg \min _\boldsymbol{x}\left(\left\|\boldsymbol{d}_{\mathrm{obs}}-\boldsymbol{D} \boldsymbol{x}\right\|_2^2+\lambda\|\boldsymbol{x}\|_1\right) \\ \tilde{\boldsymbol{d}}=\boldsymbol{S}^{-1} \tilde{\boldsymbol{x}}_\lambda \end{array}\right. $

(6)

式中: ‖·‖₁代表L₁范数: ‖·‖₂代表L₂范数, λ>0代表正则化参数。

在求解该类问题时, 重建质量主要取决于正则化参数λ的取值, 可以通过不断调整参数的取值达到更高精度的需求, 首先选取比较大的λ求得稀疏逼近解, 然后利用指数下降的方式λ^k=0.9^kλ_max, λ_max=max(|D^－1d_obs|)减小λ的取值, 使得结果不断逼近真实解。

$ \lambda=\sup _\lambda\left\{\lambda:\left\|\boldsymbol{d}_{\mathrm{obs}}-\boldsymbol{D} \boldsymbol{x}\right\|_2 \leqslant \varepsilon\right\} $

(7)

式中: ε代表重建误差。当重建结果与原始数据差值小于该值时退出迭代循环。

采用Landweber迭代阈值算法^[20]结合冷却策略(cooling strategy)^[21]的方式求解(6)式, 通过迭代减小(5)式中二次方项的大小, 进而得到Shearlet稀疏系数x, 迭代公式为:

$ x^k=T_\lambda\left[x^{k-1}+\boldsymbol{D}^{-1}\left(\boldsymbol{d}_{\mathrm{obs}}-\boldsymbol{D} x^{k-1}\right)\right] $

(8)

式中: 阈值函数T_λ=: sgn(x)·max(0, |x|－|λ|); k代表迭代次数。表 1算法伪代码展示了Landweber迭代阈值算法的基本流程。

1.2 基于卷积神经网络的数据重建方法

卷积神经网络主要包括前向传播和反向传播两种传播算法。前向传播算法是将输入数据从前向后推导形成输出数据的算法, 由于需要给定初始权重矩阵W和偏置矩阵b, 所以网络的输出结果与真实结果存在偏差。反向传播算法则与前向传播算法方向相反, 采取从后向前进行推导的方式, 实现从输出层、隐藏层到输入层的计算, 通过对比进行权重矩阵的更新, 减小输出结果与真实结果的差值。

在利用深度学习方法进行地震数据重建的过程中, 通常将随机缺失后的地震道作为输入数据x_inp, 完整的地震数据作为标签数据y_lab, 通过大量数据的学习获得输入数据和标签数据的映射关系, 从而恢复缺失的地震数据。

本文利用DnCNN网络结构进行训练^[22], 该模型是一种基于残差学习的网络模型, 该网络不是直接学习输入数据和标签数据之间的映射, 而是学习输入数据和残差(输入数据和标签数据的差值)之间的映射关系。由于残差图像大部分值较小, 网络更容易优化, 能够很好地避免网络训练中的梯度弥散问题, 其基本结构如图 1所示。

该网络包括卷积层(Conv), 批规范化层(batch normalizing, BN)和非线性化层(ReLU), 利用CNN方法对地震数据进行重建的数学表达式定义为:

$ \boldsymbol{y}_{\mathrm{out}}=\operatorname{Net}\left(\boldsymbol{x}_{\mathrm{inp}}, \theta\right) $

(9)

式中: y_out代表CNN的重建结果; Net代表网络模型; θ代表网络参数。为了增强网络的非线性, 在每层处理前加入非线性激活函数ReLU(·), 其具体数学表达式为:

$ \operatorname{ReLU}(\boldsymbol{x})= \begin{cases}x & \boldsymbol{x} \geqslant 0 \\ 0 & \boldsymbol{x} <0\end{cases} $

(10)

同时该网络引入了BN层^[23], 用以减少内部协变量转移现象(internal covariate shift)。BN层能将每一层的输入变成均值为0, 标准差为1的呈正态分布的输入数据, 进而减少初始参数对模型训练的影响, 加快网络收敛速度, 减少训练时间, 一定程度上减少过拟合。BN层的处理过程如表 2所示。

1.3 压缩感知与深度学习相结合的地震数据重建方法

输入数据和标签数据之间的关联性越大或者两者之间的不确定性越小, 利用卷积神经网络重建出的地震数据效果越好。因此, 为增强输入数据和标签数据之间的相关性, 本文将基于模型驱动的压缩感知方法和基于数据驱动的深度学习方法相结合, 经过模型驱动获得深度学习的输入数据, 以及数据驱动完成地震信号重建的复杂迭代, 给出了求解地震数据重建问题的新策略。

模型驱动和数据驱动相结合的目标函数O(m, Net)可以写为:

$ \left\{\begin{array}{l} O(\boldsymbol{m}, \mathrm{Net})=\left\|\boldsymbol{d}_{\mathrm{obs}}-F(\boldsymbol{m})\right\|+P(\boldsymbol{m})+\| \boldsymbol{m}- \\ \quad \operatorname{Net}\left(\boldsymbol{d}_{\mathrm{inp}}\right)\|+\| \boldsymbol{m}_{\text {lab }}-\operatorname{Net}\left(\boldsymbol{d}_{\mathrm{inp}}\right) \|=O_1(\boldsymbol{m})+ \\ O_2(\boldsymbol{m})+O_3(\mathrm{Net}) \\ O_1(\boldsymbol{m})=\left\|\boldsymbol{d}_{\mathrm{obs}}-F(\boldsymbol{m})\right\|+P(\boldsymbol{m}) \\ O_2(\boldsymbol{m})=\left\|\boldsymbol{m}-\operatorname{Net}\left(\boldsymbol{d}_{\mathrm{inp}}\right)\right\| \\ O_3(\mathrm{Net})=\left\|\boldsymbol{m}_{\mathrm{lab}}-\operatorname{Net}\left(\boldsymbol{d}_{\mathrm{inp}}\right)\right\| \end{array}\right. $

(11)

式中: m代表待求的重建数据; P代表正则化约束; F代表对数据的采集; d_inp代表CNN的输入数据; m_lab代表CNN的标签数据。

对于(11)式而言, 同时得到最优的m和Net有些困难。因此, 采用分步求解的算法, 即首先通过O₁(m)对数据初步重构, 为O₃(m)提供合适的输入数据d_inp, 使得O₃(m)中的映射简单, 网络易优化。然后, 利用O₃(m)构建深度学习网络Net, 基于上述两步即可获得合适的输入数据d_inp和深度学习网络Net。最后, 利用O₂(m)得到最终的重构数据m, 图 2为该方法的网络结构示意。

对于多层卷积神经网络结果, 可以在(9)式的基础上展开为:

$ \left\{\begin{array}{l} \operatorname{Net}_1\left(\boldsymbol{d}_{\mathrm{inp}} ; \theta_1\right)=\operatorname{ReLU}\left(\boldsymbol{W}_1 * \boldsymbol{d}_{\mathrm{inp}}+\boldsymbol{b}_1\right) \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ldots \\ \operatorname{Net}_l\left(\boldsymbol{d}_{\mathrm{inp}} ; \theta_l\right)=\operatorname{ReLU}\left(\boldsymbol{W}_l * \operatorname{Net}_{l-1}\left(\boldsymbol{d}_{\mathrm{inp}} ; \theta_{l-1}\right)+\boldsymbol{b}_{l-1}\right) \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ldots \\ \operatorname{Net}_L\left(\boldsymbol{d}_{\mathrm{inp}} ; \theta_L\right)=\boldsymbol{W}_L * \operatorname{Net}_{L-1}\left(\boldsymbol{d}_{\mathrm{inp}} ; \theta_{L-1}\right)+\boldsymbol{b}_{L-1} \end{array}\right. $

(12)

式中: l=1, …, L; l代表网络层数; W_l代表第l层网络的权重矩阵; b_l代表第l层网络的偏置矩阵。

为了评估网络预测结果与真实值的误差, 引入损失函数。损失函数的值越小, 代表模型的鲁棒性越好, 采用均方差函数作为损失函数, 其数学表达式为:

$ l(\theta)=\frac{1}{2 N} \sum\limits_{i=1}^N\left\|\operatorname{Net}\left(\boldsymbol{d}_{\mathrm{inp}}^i, \theta\right)-\boldsymbol{m}_{\mathrm{lab}}^i\right\|_2^2 $

(13)

式中: N代表批处理数据的样本数; Net(d_inpⁱ, θ)代表输入数据d_inpⁱ经过网络后的预测值; m_labⁱ代表对应的完整地震数据。采用自适应矩估计Adam方法^[24]优化损失函数, 即使在内存较低的情况下也能进行高效计算。

2 模拟数据试验

对合成的叠前地震反射数据进行网络训练, 训练数据集包括4种模拟地震数据, 即1994 BP静态基准模型(1994 BP statics benchmark model)数据、1997 BP 2.5D偏移基准模型(1997 BP 2.5D migration benchmark model)数据、2007 BP各向异性速度基准模型(2007 BP anisotropic velocity benchmark)数据和Hess VTI偏移基准模型(Hess VTI migration benchmark)数据, 整合后共有11059个单炮记录。由于相邻单炮记录特征相似, 并且考虑到计算机的性能, 在11059个单炮记录中每10个单炮记录中选取1个, 共选取1106个单炮记录。再从中选取1086个单炮记录作为训练数据, 选取20个单炮记录作为测试数据。图 3为模型训练的收敛曲线。对各单炮数据进行随机缺失50%地震道的处理, 随后对其进行压缩感知重建, 得到CNN模型的输入数据, 将完整的单炮数据作为标签数据, 并将其对应裁剪为50×50的数据块(patch)学习地震数据局部特征。

利用DnCNN网络结构进行训练。该网络具有17个卷积层。第1层由卷积层和非线性化层组成, 第2层至第16层的每一层均由卷积层、批规范化层和非线性化层组成, 最后一层由一个卷积层组成。每层卷积核大小均为3×3, 第一层卷积核的通道数为1, 个数为64, 第2层至第16层卷积核的通道数为64, 个数为64, 最后一层卷积核的通道数为64, 个数为1, 共进行50轮训练。前30轮的学习率为0.001, 后20轮的学习率为0.0001。在进行卷积运算的时候, 利用填零方法使得每层输入数据与输出数据大小相同。训练环境为Windows 10系统的工作站, 处理器为至强W-2245 3.9GHz 64GB的内存和10GB的NVIDIA RTX3080 GPU, 网络训练总共需要10个小时。其中, 结合压缩感知和深度学习方法在CNN方法的基础上, 加入了压缩感知预处理环节, 训练时间增加40分钟。

为了评价地震数据的重建效果, 引入信噪比公式:

$ \mathrm{SNR}=20 \lg \left(\frac{\|\boldsymbol{f}\|_2}{\|\boldsymbol{f}-\hat{\boldsymbol{f}}\|_2}\right) $

(14)

式中: f代表完整数据; $\hat{\boldsymbol{f}}$代表重建后的数据。

为对比验证不同方法的效果, 对20个测试数据进行随机缺失70%地震道处理, 随后对该缺失数据分别进行压缩感知方法重建、CNN方法重建以及结合压缩感知和深度学习方法重建。其中, 压缩感知方法迭代次数为50次, 而结合压缩感知和深度学习方法在前期处理时, 迭代次数仅为20次。将3种方法得到的结果从信噪比和计算时间两个角度分别进行对比, 结果如表 3所示。

由表 3可看出, 当随机缺失70%地震道数据时, 压缩感知方法已经不能恢复缺失数据, 信噪比均值为10.61dB, 与完整数据差值较大。而CNN方法及结合压缩感知和深度学习方法(本文方法)都能对缺失数据进行有效恢复, 能够较好地重建出完整的地震数据, 信噪比均值分别为21.53dB和26.76dB。结合压缩感知和深度学习方法在压缩感知理论的基础上结合深度学习方法对部分缺失数据进行进一步处理, 具有更高的重建精度。同时在测试阶段, 虽然结合压缩感知和深度学习方法的计算时间长于CNN方法, 但是时间比传统的压缩感知方法少50%左右。

为了进一步验证结合压缩感知和深度学习方法的有效性, 选取另外一组测试数据分别进行随机缺失30%, 50%, 70%地震道处理, 随后分别利用压缩感知方法、CNN方法以及结合压缩感知和深度学习方法对数据进行重建。该测试数据大小为1000×499;时间采样间隔为6ms; 道间距为12.5m;如图 4所示。为了更进一步直观地展现3种方法的重建结果和残差剖面, 选取道号300~500, 时间为400~600ms的局部数据进行对比(图 4黑框区域)。图 5至图 7分别展现了随机缺失30%, 50%, 70%地震道后黑框区域内的重建结果和残差剖面。

从图 5至图 7重建结果和残差剖面可以看出, 在缺失少量数据的情况下, 3种方法都能对缺失信息进行有效重建, 但是随着缺失数据增多, 压缩感知重建的结果逐渐变差。在随机缺失70%地震道数据时, 出现同相轴间断、残差呈挂面条状分布的现象(图 7d), 已经不能很好地重建出完整的地震数据。然而, CNN方法及结合压缩感知和深度学习方法(本文方法)相较于压缩感知方法, 在数据大量缺失的情况下, 仍可以得到较完整的重建结果。本文方法在CNN方法的基础上结合模型驱动的先验知识, 对比图 7e和图 7f中的箭头部分可以看出, 本文方法能够进一步刻画微小特征, 能够基本重建出完整的地震数据。

表 4列出了在3种缺失情况下, 利用(14)式计算出的3种方法的信噪比结果以及重建所需的计算时间。可以看出, 无论是在低缺失还是高缺失的情况下, 本文方法的重建精度都比其它方法的重建精度高。

此外, 基于数据驱动的CNN方法及本文方法需要一定时间训练网络参数, 但在重建时能够快速地获得高精度的结果, 本文方法在CNN方法的基础上只需增加部分训练时间便可以进一步改善重建效果。

3 实际数据试验

选取实际海洋地震数据^[25]进行重建, 所选数据大小为512×430, 时间采样间隔为4ms, 道间距为12.5m, 实际地震数据与随机缺失50%地震道数据如图 8所示。

分别采用压缩感知方法、CNN方法及本文方法对图 8b中缺失的数据进行重建, 其结果以及残差剖面如图 9所示, 这3种方法重建所需计算时间分别为5.21s, 1.18s和3.52s。从图 9可以看出, 3种方法都能很好地重建出实际地震数据, 但是本文方法重建出的结果横向连续性最强, 能够恢复出大量微小的特征。从信噪比角度而言, 3种方法在随机缺失50%地震道的情况下的信噪比分别为12.62dB, 21.25dB和24.37dB。由此可见, 本文方法重建出的地震数据信噪比最高, 与实际数据的差值最小。

4 结论与认识

本文将基于模型驱动的压缩感知方法和基于数据驱动的深度学习方法相结合, 提出了基于压缩感知和深度学习的地震数据重建方法。在网络中融入模型的先验知识, 利用大量数据进一步学习模型先验知识, 进而有效降低了压缩感知方法对于先验知识的要求, 提高了数据在大量缺失情况下的重建精度。合成地震数据和实际地震数据应用结果表明, 与传统的压缩感知方法和深度学习方法相比, 结合压缩感知和深度学习方法的重建效果和重建精度都有很大提高, 并且在数据大量缺失的情况下, 也能得到很好的结果。在计算效率方面, 本文方法的计算时间主要集中于网络训练阶段, 当网络训练完成后, 重建效率比传统压缩感知方法提高了一倍。同时, 采用数据块处理的方式构建训练数据, 能够更好地学习地震数据的局部特征, 具有较好的泛化性。将模型驱动和数据驱动进行初步串联获得了很好的结果。如何将模型驱动和数据驱动进一步有机结合是下一步的研究方向。同时本文方法在训练过程中计算时间较长, 如何在保证高精度的条件下提高计算效率是另一个研究方向。