0 引言

地震数据重建是地震信号分析领域的一个重要研究方向。地震数据采集时由于采集环境以及采集成本等因素^[1]，采集到的数据通常呈现不规则或稀疏分布，不完整的地震数据缺失部分地球物理信息，据其不能准确地预测地下岩层的性质。因此，为保证满足后续地震资料处理对数据完整性的要求，必须恢复缺失的地震数据，以提高储层预测的准确性，促进油气的勘探与开发。显然，重建缺失的地震数据具有非常重要的现实意义。

传统的地震数据重建方法主要有基于滤波的重建方法^[2-4]、基于波动方程的重建方法^[5-7]以及基于变换函数的重建方法^[8-11]等三种类型。滤波重建法是通过褶积插值滤波器实现地震数据重建，主要用于重建规则采样地震数据，且计算量大。波动方程法是通过Kirchhoff积分算子实现地震数据重建^[12]，该方法需地下信息作为先验条件，在地下信息未知或精度低时对重建结果影响很大。变换函数法是利用最小二乘法将信号变换到某一新的变换域，再通过反变换重构地震数据。常用的变换函数有傅里叶变换、曲波变换、离散余弦变换等。变换函数法结构简单，计算复杂度低，可实现规则和非规则采样地震数据的重建，且不需将地下地层信息作为先验条件，但该方法对采样率要求较高，造成采集成本高昂，且要求地震数据的采集规模较大。

压缩感知理论^[13-16](Compressive sensing，CS)的应用，克服了传统地震数据重建方法的缺陷，解决了对采样率要求较高的问题^[17]，使稀疏信号可在低于奈奎斯特采样频率的条件下用稀疏优化算法重建完整信号。压缩感知地震数据重建主要分为稀疏表示阶段和稀疏重建阶段，重建阶段相当于稀疏优化问题。白兰淑等^[18]结合稀疏促进的曲波恢复迭代算法和Bregman迭代阈值算法提出改进的联合迭代算法，在压缩感知理论框架下，用Curvelet变换稀疏表示地震数据，该改进算法能稀疏优化重建地震数据，具有收敛快、重建信噪比高的优点，但其计算量大、计算复杂度高。周亚同等^[19]提出压缩感知下K-SVD字典学习自适应稀疏表示地震数据、正则化正交匹配追踪(Regularized orthogonal matching pursuit，ROMP)稀疏优化重建地震数据的方法，能加快复杂实际地震数据重建速度。但该稀疏优化算法得到的不是全局最优解，故重建精度较低。

交替乘子方向算法(Alternating direction method of multipliers，ADMM)^[20-22]具有简化复杂性、高重建速度和良好重建质量的优点，在地震数据领域引起了广泛关注。Sternfels等^[23]将ADMM算法与压缩感知理论结合实现含噪地震数据的降噪。Zhang等^[24]将ADMM算法应用于缺失地震数据的恢复，实现压缩感知地震数据的精确重构。李慧等^[25]在此基础上提出压缩感知框架下K-SVD字典学习稀疏表示地震数据，利用ADMM算法优化重建地震数据，达到提高地震数据重建精度的目的。

尽管ADMM算法在地震数据重建中效果显著，但由于该算法中平衡因子无法根据地震数据自适应地进行调整，易损害重建精度。另外，在基于ADMM算法的压缩感知地震数据重建模型中，由于输入训练样本较少，算法参数较多以及模型复杂度较高等问题，容易出现过拟合现象。针对上述问题，本文对ADMM算法进行改进，提出平方正则交替乘子方向算法(Square regular-alternating direction method of multipliers，SR-ADMM)。一方面在ADMM算法迭代中加入平方正则项；另一方面自适应选取平衡因子。最终将SR-ADMM算法与压缩感知理论相结合，实现基于压缩感知的SR-ADMM算法地震数据重建。通过模拟地震数据、实际单炮数据以及大庆油田实际地震数据对本文所提算法进行了验证。

1 压缩感知基本理论

压缩感知理论具有三个重要的前提条件：地震数据的稀疏性、测量矩阵的不相干性以及合适的优化算法。压缩感知的数据重建模型为

$ \mathit{\boldsymbol{y = \varPhi x}} $

(1)

式中：x∈R^N，为原始信号；Φ∈R^M×N，为测量矩阵；y∈R^M，为测量矩阵下原始信号的观测值。该式表示原始信号从N维信号降维到M维获得观测值，信号降维过程可用图 1所示图像表示。

由于式(1)为欠定方程，有无穷多的解，很难直接通过观测值y恢复出原始信号x。因此可用一个与测量矩阵不相关的变换矩阵对原始信号做稀疏表示，得到稀疏矩阵。其数学表达式为

$ \mathit{\boldsymbol{x = D\varTheta }} $

(2)

式中：D∈R^N×N，为稀疏变换矩阵；Θ为稀疏系数，是长度为N的一维向量。式(2)的原始信号稀疏表示过程可用图 2所示图像表示。

因此观测值y可以表示为

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{\varPhi x}} = \mathit{\boldsymbol{\varPhi \varPhi \varTheta }} = \mathit{\boldsymbol{A\varTheta }} $

(3)

式中A=ΦD，为M×N阶感知矩阵。A满足约束等距性质，即可通过求解L₀范数优化问题很好地从观测值y中恢复出原始信号x

$ \min {\left\| \mathit{\boldsymbol{\varTheta }} \right\|_0}\quad {\rm{ s}}{\rm{.t}}{\rm{. }}\quad \mathit{\boldsymbol{y = A\varTheta }} $

(4)

式中稀疏系数Θ的L₀范数表示该向量中所有非零元素的个数。但L₀范数的最优化问题是一个非确定性多项式问题，在稀疏变换矩阵D与测量矩阵Φ不相干时，求解Θ的L₀范数问题等价于求解Θ的L₁范数，故可将式(4)转化为求解L₁范数优化问题

$ \min {\left\| \mathit{\boldsymbol{\varTheta }} \right\|_1}\quad {\rm{ s}}{\rm{.t}}{\rm{. }}\quad \mathit{\boldsymbol{y = A\varTheta }} $

(5)

式(5)与式(4)等价，即求解满足约束条件的最稀疏Θ转换为求解满足约束条件的元素绝对值和最小Θ。通过求解式(5)得到最稀疏系数$\mathit{\boldsymbol{ \boldsymbol{\tilde \varTheta} }}$，最后根据$\mathit{\boldsymbol{ \boldsymbol{\tilde \varTheta} }}$恢复重建原始信号

$ \tilde{\boldsymbol{x}}=\tilde{\boldsymbol{D}} \tilde{\boldsymbol{\varTheta}} $

(6)

2 SR-ADMM算法

适用的重构算法是压缩感知精确重构的重要前提。ADMM算法的优点在于将对偶上升法的可分解性与乘子法的收敛性质结合起来，通过对偶上升法使线性稀疏项不断接近最优解，保证在惩罚参数较小的情况下满足精度要求；克服乘子法中变量二次项无法拆分的问题，将原始问题拆分求解，可实现并行运算且提高运算效率。使用ADMM算法求解压缩感知地震数据重建问题的数学模型为

$ \widetilde{\boldsymbol{\varTheta}}=\operatorname{argmin} \frac{1}{2}\|\boldsymbol{A} \boldsymbol{\varTheta}-\boldsymbol{y}\|_{2}^{2}+\boldsymbol{\lambda}\|\boldsymbol{\varTheta}\|_{1} $

(7)

式中λ表示平衡因子，用于控制前、后两项之间的权重。

引入辅助变量Z

$ \begin{gathered} \widetilde{\boldsymbol{\varTheta}}=\operatorname{argmin} \frac{1}{2}\|\boldsymbol{A} \boldsymbol{\varTheta}-\boldsymbol{y}\|_{2}^{2}+\boldsymbol{\lambda}\|\boldsymbol{Z}\|_{1} \\ \text { s.t. } \quad \boldsymbol{\varTheta}-\boldsymbol{Z}=0 \end{gathered} $

(8)

该式的增广拉格朗日函数为

$ \begin{gathered} L_{\rho}(\boldsymbol{\varTheta}, \boldsymbol{Z}, \boldsymbol{y})=\frac{1}{2}\|\boldsymbol{A} \boldsymbol{\varTheta}-\boldsymbol{y}\|_{2}^{2}+\boldsymbol{\lambda}(\boldsymbol{\varTheta}-\boldsymbol{Z})+ \\ \frac{\rho}{2}\|\boldsymbol{\varTheta}-\boldsymbol{Z}\|_{2}^{2} \end{gathered} $

(9)

根据上式可知，Θ的更新就是在固定参数Z和λ的情况下，求解增广拉格朗日函数的最小值问题，即

$ \boldsymbol{L}^{\prime}{}_{\boldsymbol{\varTheta}}=\boldsymbol{A}^{\mathrm{T}}\left(\boldsymbol{A} \boldsymbol{\theta}^{k+1}-\boldsymbol{y}\right)+\boldsymbol{\lambda}^{k}+\rho\left(\boldsymbol{\theta}^{k+1}-\boldsymbol{z}^{k}\right) $

(10)

求解可知ADMM算法第k+1次的迭代表达式为

$ \left\{\begin{array}{l} \boldsymbol{\theta}^{k+1}=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}+\rho \boldsymbol{I}\right)\left[\boldsymbol{A}^{\mathrm{T}} \boldsymbol{y}+\rho\left(\boldsymbol{z}^{k}-\boldsymbol{\lambda}^{k}\right)\right] \\ \boldsymbol{z}^{k+1}=\boldsymbol{S}_{\boldsymbol{\lambda} / \rho}\left(\boldsymbol{\theta}^{k+1}+\frac{\boldsymbol{\lambda}^{k}}{\rho}\right) \\ \boldsymbol{u}^{k+1}=\boldsymbol{u}^{k}+\left(\boldsymbol{\theta}^{k+1}-\boldsymbol{z}^{k+1}\right) \end{array}\right. $

(11)

式中：u=λ/ρ，其中ρ表示惩罚参数；S_λ/ρ为软阈值函数，定义为

$ \boldsymbol{S}_{\boldsymbol{\lambda} / \rho}(x)= \begin{cases}x-k & x>k \\ 0 & |x| \leqslant k \\ x+k & x<-k\end{cases} $

(12)

以K-SVD字典学习对地震数据进行稀疏表示，并用ADMM算法解决稀疏优化问题，构成基于ADMM算法和K-SVD字典学习的压缩感知地震数据重建模型。在该模型中，由于存在输入训练样本较少、ADMM算法迭代参数较多及模型复杂度较高等问题，容易出现过拟合现象，因此增加平方邻近项正则化式中的子问题，防止过拟合。正则化交替乘子方向算法公式如下

$ \left\{\begin{array}{l} \boldsymbol{\theta}^{k+1}=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}+\rho \boldsymbol{I}\right)\left[\boldsymbol{A}^{\mathrm{T}} \boldsymbol{y}+\rho\left(\boldsymbol{z}^{k}-\boldsymbol{\lambda}^{k}\right)\right] \\ \boldsymbol{z}^{k+1}=\boldsymbol{S}_{\boldsymbol{\lambda} / \rho}\left(\boldsymbol{\theta}^{k+1}+\frac{\boldsymbol{\lambda}^{k}}{\rho}\right)+\frac{1}{2}\left\|\boldsymbol{z}-\boldsymbol{z}^{k}\right\|^{2} \\ \boldsymbol{u}^{k+1}=\boldsymbol{u}^{k}+\left(\boldsymbol{\theta}^{k+1}-\boldsymbol{z}^{k+1}\right)+\frac{1}{2}\left\|\boldsymbol{u}-\boldsymbol{u}^{k}\right\|^{2} \end{array}\right. $

(13)

另外，由于ADMM算法中平衡因子的取值只能根据经验选取，若平衡因子λ较大，则说明稀疏项权重更大，所求值更稀疏而不会过多考虑约束；反之，重建数据更接近于原始信号。为了使λ更精确，从而获得更高重建精度，本文将λ的选取改进为可自适应方式。对平衡因子的改进方案如下

$ \boldsymbol{\lambda}_{1, a}=\boldsymbol{\lambda} \times \alpha^{l}, \boldsymbol{\lambda}_{2, a}=\boldsymbol{\lambda} \times \alpha^{l-1}, \cdots $

(14)

$ \boldsymbol{\lambda}_{1, b}=\boldsymbol{\lambda} / \alpha^{l}, \boldsymbol{\lambda}_{2, b}=\boldsymbol{\lambda} / \alpha^{l-1}, \cdots $

(15)

首先根据经验设定初始平衡因子λ。两式中α的作用是选取λ的邻近值，可选择略大于1或略小于1的数，如α=0.95。式(14)和式(15)分别表示对初始平衡因子做乘法和除法，获得多个近似值。接着分别用所获近似值以式(13)进行迭代实现重建，选取使地震数据重建误差最小的平衡因子为最优平衡因子，通过一轮迭代就可得到平衡因子最优值；之后，继续完成后续迭代过程。两式中常数l表征近似解的个数，若l取值越大，就能获得数量更多的平衡因子近似值，理论上重建效果更好。但l过大，会造成时间空间负荷，因此需根据地震数据实情选取。

3 基于压缩感知的SR-ADMM地震数据重建

本文基于压缩感知的SR-ADMM算法地震数据重建的流程如图 3所示。

在压缩感知框架下，以K-SVD字典学习对地震数据做稀疏表示，用SR-ADMM算法稀疏优化实现缺失地震数据重建。缺失数据重建具体步骤为：

(1) 输入训练样本，对重建模型进行训练，得到初始字典D₀、初始稀疏系数Θ₀；

(2) 输入原始缺失地震数据x，用初始字典D₀对该地震数据进行稀疏表示；

(3) 用SR-ADMM算法求解最优化问题(式(7))，得到初步稀疏近似$\mathit{\boldsymbol{ \boldsymbol{\tilde \varTheta} }}$，求解$\mathit{\boldsymbol{\tilde x = D \boldsymbol{\tilde \varTheta} }}$得到地震数据的初步重建；

(4) 以K-SVD字典学习算法更新稀疏变换字典D，用更新后的字典D对步骤(3)初步重建地震数据做稀疏表示；

(5) 用SR-ADMM算法求解更新得到稀疏近似$\mathit{\boldsymbol{ \boldsymbol{\tilde \varTheta} }}$，更新$\mathit{\boldsymbol{\tilde x = D \boldsymbol{\tilde \varTheta} }}$；

(6) 将步骤(3)迭代得到的地震数据替换为下一次迭代的输入数据，重复步骤(4)~步骤(5)，直到收敛；输出迭代后的重建地震数据$\mathit{\boldsymbol{\tilde x}}$，实现缺失地震数据的重建。

4 地震数据仿真验证

将信噪比作为地震数据重建性能的评判指标

$ \mathrm{SNR}=10 \lg \frac{\|\boldsymbol{x}\|_{2}^{2}}{\|\boldsymbol{x}-\tilde{\boldsymbol{x}}\|_{2}^{2}} $

(16)

式中$\mathit{\boldsymbol{\tilde x}}$为由原始地震数据x得到的重建数据。

4.1 模拟地震数据重建仿真验证

验证模型的模拟地震数据(图 4)共76道，522个采样点，采样间隔为1ms。图 5为该炮点集随机采样数据。ADMM算法中平衡因子根据多次实验后得出，选取λ=0.15；SR-ADMM算法平衡因子则根据模拟地震数据自适应地选取。

OMP算法是压缩感知框架下常用于解决优化问题的一种稀疏优化算法。分别采用OMP(图 6a)、ADMM(图 6b)和SR-ADMM(图 6c)三种算法对图 5所示随机采样模拟地震数据进行重建。

上述三种算法的重建信噪比和均方误差如表 1所示。将这三种算法重建结果(图 6，深蓝色波形)与原始地震数据(黑色波形)进行对比，并将缺失较严重的第19~第33道数据进行局部放大。可见对于该随机缺失模拟地震数据，基于压缩感知理论的三种算法都能实现数据重建，OMP算法(图 6a)重建效果一般，从随机缺失的第21道~第26道数据，易见缺失部分重建振幅与原始振幅存在误差，且该算法对于未缺失部分(如第45道)数据产生影响，重建的均方误差相对其他两种方法较大，重建信噪比相对其他两种方法最低；ADMM算法(图 6b)重建效果较好，不会对未缺失部分数据产生影响，但缺失数据的振幅尚未得到完全恢复；SR-ADMM算法(图 6c)重建的效果最好，对数据细节的恢复更完备，同相轴清晰可见，重建数据振幅与原始地震数据振幅一致，且该算法重建信噪比最高，重建均方误差最小。

由于SR-ADMM算法中平衡因子是根据地震数据自适应地选取，所选平衡因子使数据重建误差最小，从而达到提高重建精度的目的。而且，基于该算法的压缩感知数据重建模型防止了数据的过拟合现象，一定程度上也提高了数据重建精度。因此，针对随机缺失模拟地震数据，采用SR-ADMM算法重建的信噪比高于OMP和ADMM算法，即SR- ADMM算法的重建精度高于经典ADMM算法。

4.2 实际地震数据仿真验证

选取的M油田实际单炮地震数据(图 7)共有170道，采样点数为1300，时间采样率为1ms；其缺失率为40%的单炮地震数据如图 8a所示。分别使用OMP、ADMM及SR-ADMM三种算法对缺失数据(图 8a)进行重建(图 8b~图 8d)。ADMM算法中平衡因子据多次试验后设定为λ=0.17，SR-ADMM算法的平衡因子据实际地震数据自适应地选取。三种算法对应的重建信噪比和均方误差如表 2所示。

结合表 2和图 8b~图 8d可看出：OMP算法(图 8b)重建效果欠理想，同相轴连续性较差，重建信噪比较低，均方误差较大；ADMM算法(图 8c)重建效果比OMP算法稍好，大致实现了缺失数据的重建，但重建后的数据附近出现抖动，重建信噪比较高，均方误差中等；SR-ADMM算法(图 8d)重建效果很好，同相轴更光滑，重建振幅连续性良好，在三种算法中重建信噪比最高，均方误差最小。

为了更直观地对比三种优化算法的重建效果，将随机采样原始数据中缺失较多的第59~第109道进行局部放大(图 9a)，其中第59~第67道、第72~第74道、第76道、第78~第87道、第90~第95道、第100道及第104道为缺失数据。分别采用OMP(图 9b)、ADMM(图 9c)和SR-ADMM(图 9d)三种算法得到相应重建结果的局部放大图。

对比图 9a与图 9b可看出部分缺失数据实现了重建(图 9b蓝色箭头所指)，部分缺失数据未得到恢复(红色箭头所指)，且重建数据振幅与原始数据有差距；对比图 9a与图 9c可见，缺失数据虽得到了恢复(图 9c蓝色箭头)，但重建后的振幅与原始数据存在误差；对比图 9a与图 9d发现，所有缺失数据都得到有效恢复(图 9d蓝色箭头)，且重建数据振幅与原始数据一致。

因此，采用SR-ADMM算法可以很好地重建实际地震数据复杂的细节，重建性能较好。

4.3 三维实际地震数据应用

上述模拟和实际地震单炮记录验证了本文方法的有效性，进一步将本文方法应用于中国石油大庆油田实际三维数据。从该三维地震数据的一条测线中选取200道(图 10)，时间采样率为1ms，采样点为500个。图 11a为随机缺失30%的该三维数据实例。ADMM算法中平衡因子据多次实验后设定为λ=0.20，SR-ADMM算法的平衡因子基于实际地震数据自适应地选取。

同样分别采用OMP、ADMM及SR-ADMM三种算法对随机采样实际三维地震数据进行重建，得到相应的重建结果(图 11b~图 11d)。可见OMP算法的数据重建(图 11b)效果一般，未能很好地实现缺失数据的恢复重建，同相轴欠连续；ADMM算法重建数据(图 11c)效果相对较好，同相轴连续性提高，但存在同相轴抖动现象；SR-ADMM算法很好地实现了三维数据重建(图 11d)，缺失位置的数据与相邻位置同相轴的走向相符合，且很连续。

为便于比较三种算法重建效果，将三种重建结果的第50~第60道与原始数据进行放大对比(图 12)。三种算法的重建信噪比及均方误差如表 3所示。

综合分析表 3和图 12所示局部放大图，可以看出OMP算法恢复地震数据(图 12a)的振幅与原始地震数据(黑色)相差较大，且重建信噪比较低，均方误差较大；ADMM算法重建(图 12b)的振幅与原始地震数据更契合，其重建信噪比高于OMP算法，均方误差较小；而SR-ADMM算法重建后(图 12c)振幅与原始地震数据几乎一致，重建信噪比最高，均方误差最小。显然，SR-ADMM算法在实际三维地震数据重建中也有效，且可提高实际三维地震数据的重建精度。

5 结束语

本文提出基于压缩感知的SR-ADMM算法地震数据重建方法。该算法在原始ADMM算法基础上加入正则项正则化迭代中子问题，并且自适应地选取平衡因子，利用该算法解决最优化问题实现缺失地震数据的重建。仿真验证表明利用本文算法能有效地实现缺失地震数据，且其重建精度高于常规压缩感知地震数据重建优化算法。

需要指出的是，数据的稀疏性是压缩感知数据重建的前提，对重建效果有很大影响，数据越稀疏，重建效果越好。因此，如何更好的对地震数据进行稀疏表示是下一步的研究方向。