地震数据采集中, 常规采样方式得到的地震数据存在数据缺失的问题。根据香农采样定理, 为了不失真地重建地震数据, 地震数据采样频率应该不小于数据最高频率的两倍。传统的地震数据缺失重建方法有很多, SPITZ^[1]提出了在f-x域进行插值重构的方法, 但是这种方法的计算量大; RONEN^[2]根据波动方程, 利用NMO与逆DMO联合重构缺失地震数据, 但当参数未知时, 该方法需要手动调整参数, 结果不精确。ZWARTJES等^[3]提出了在频率-波数域进行插值的地震数据重建方法; HERRMANN等^[4]提出了基于曲波变换的地震数据恢复算法; NAGHIZADEH等^[5]提出了使用多维预测滤波的方法来重建地震数据, 该方法首先采用最小加权范数插值法将重要的低频信息正则化, 然后再利用重建的数据对高频信息进行插值重建。

与传统方法不同, 压缩感知理论下的信号采样方式突破了香农采样定理的限制, 即使采样率很低也可以重建地震数据。2011年, WANG等^[6]将基于最小化L1范数的重建算法用于压缩感知下地震数据的重建, 并采用小波变换对信号稀疏表示, 取得了很好的重建效果。

利用压缩感知方法进行地震数据重建, 必须满足3个条件:地震数据具有稀疏性、合理的测量矩阵以及非线性重建算法^[7-11]。也就是要求测量算子和稀疏变换算子高度不相干, 还需要选择合适的重建算法, 本文主要研究重建算法。压缩感知中常用的重建算法包括L0范数类算法(如OMP、迭代式硬阈值算法等)以及最小化L1范数类算法^[12](如SPGL1算法、FPC算法等)。其中, L0类算法大多都以稀疏度K作为信号成功重构的先决条件之一, 但在现实中很难得到稀疏度K的精确值。CANDES等^[13-14]将凸优化理论引入到压缩感知, 并编写了L1-magic软件, 但该软件运行耗时较长。

2008年, HALE等^[15]首次将FPC算法应用于大规模的稀疏优化问题, 并证明了其收敛性能优于传统的重建算法。随着压缩感知理论的发展, FPC算法也逐渐被应用到很多领域, MA等^[16]将不动点迭代与Bregman迭代结合求解矩阵最小化问题; 刘晓曼等^[17]在压缩感知框架下将FPC算法应用于核磁共振成像并取得良好的效果。FPC算法基于算子分裂的思想, 属于迭代收缩类算法, FFPC算法是FPC算法的改进版, 引入了软阈值步长参数和正则化参数的收缩参数。与FPC算法相比, FFPC算法的收敛速度更快, 耗时更短。本文将FFPC算法应用于缺失地震数据的重建, 研究发现, 与传统的缺失地震数据重建方法相比, 无论是在重建精度还是耗时问题上, 压缩感知理论下的FFPC算法都具有明显优势, 对实际缺失地震数据的重建效果也较好。

1 压缩感知框架下基于FFPC算法的地震数据重建 1.1 压缩感知基本原理

信号的稀疏性是压缩感知能够实现的关键要素, 常用的稀疏表示方法有很多, 如小波变换、脊波变换、字典训练等, 其中小波变换使用简单, 对地震数据的稀疏表示效果也较好。在众多小波基中, Symlet小波基具有良好的对称性, 在一定程度上能够减少地震数据分析和重构时的相位失真, 因此本文采用Symlet小波基对地震数据X进行稀疏表示

$ \mathit{\boldsymbol{s}} = \mathit{\boldsymbol{\varphi X}} $

(1)

式中:φ表示Symlet小波变换算子; X∈R^N为重建后的完整地震数据; s为X在小波域上的稀疏系数。压缩感知方法的优点是能够从缺失的地震数据Y中精确地恢复出原始地震X的全部信息。将其用于地震数据重建的表达式为:

$ \mathit{\boldsymbol{Y}} = \phi {\mathit{\boldsymbol{\varphi }}^{\rm{H}}}\mathit{\boldsymbol{s}} $

(2)

式中:Y∈R^M为采集到的缺失地震数据, $N \gg M$; ϕ为测量矩阵; φ^H表示小波变换算子的共轭转置。压缩感知理论恢复地震数据的要求是稀疏矩阵和测量矩阵满足有限等距性质(restricted isometry property, RIP)^[18]。地震勘探中, 首先利用点震源激发尖脉冲, 然后利用各个方位的检波器接收地震记录, 故可认为测量过程是在狄拉克(Dirac)基中完成的, 即测量矩阵是单位矩阵。Symlet小波是近似对称的紧支集正交小波, HERRMANN等^[4]证明了Dirac基与所有的正交变换基都不相关, 所以Dirac基与Symlet小波不相关, 即地震数据满足压缩感知理论的有限等距性质。本文引入约束矩阵T, 将约束矩阵T作用于测量矩阵ϕ, 从而得到一个随机缺失的地震数据。

在公式(2)中, 我们引入感知矩阵A, 令A=Tϕφ^H, 则Y=As。这样将缺失地震数据的重建问题表示为最小L0范数问题:

$ \mathit{\boldsymbol{s}} = \arg \min {\left\| \mathit{\boldsymbol{s}} \right\|_0}\;\;\;\;{\rm{s}}.{\rm{t}}.\;\;\;\;\;{\left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{As}}} \right\|_2} \le \varepsilon $

(3)

式中:ε是误差估计参数。公式(3)中求解L0范数是一个NP(non-deterministic polynomial)难的问题^[19], DONOHO^[20]将其改为求解最小L1范数的形式:

$ \mathit{\boldsymbol{s}} = \arg \min {\left\| \mathit{\boldsymbol{s}} \right\|_1}\;\;\;\;{\rm{s}}.{\rm{t}}.\;\;\;\;\;{\left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{As}}} \right\|_2} \le \varepsilon $

(4)

1.2 FFPC算法

传统FPC方法主要用于求解下式所示的L1范数问题:

$ \mathop {\min }\limits_{\mathit{\boldsymbol{s}} \in {{\bf{R}}^N}} {\left\| \mathit{\boldsymbol{s}} \right\|_1} + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{As - Y}}} \right\|_M^2 = \theta \left( \mathit{\boldsymbol{s}} \right) $

(5)

其不动点迭代式为:

$ {\mathit{\boldsymbol{s}}^{k + 1}} = {R_v}h\left( {{\mathit{\boldsymbol{s}}^k}} \right)\;\;\;\;v = \frac{\tau }{\mu } $

(6)

式中:R_v(·)表示软阈值收缩算子; k是迭代次数; v是阈值。当公式(5)中M=2时, 我们用正则化参数μ来平衡L1范数和L2范数各自所占的比重, FPC方法的收敛速度由v=τ/μ和二次函数部分的Hessian矩阵决定。为了提高重建地震数据的质量和加快算法的收敛速率, 我们在此引入软阈值步长参数t和正则化收缩参数β。FFPC算法利用软阈值步长参数t来组合前两步迭代结果, 并通过引入正则化收缩参数β在迭代过程中来缩小正则化参数μ, 从而提高算法的收敛速度。具体步骤如下:对于不动点迭代公式(6), 我们参考了快速迭代收缩阈值(fast iterative shrinkage thresholding, FIST)算法, 利用前两次迭代结果对Z^k进行更新^[21]。

$ {\mathit{\boldsymbol{Z}}^k} = {\mathit{\boldsymbol{s}}^k} + \frac{{{t_{k - 1}}}}{{{t_{k + 1}}}}\left( {{\mathit{\boldsymbol{s}}^k} - {\mathit{\boldsymbol{s}}^{k - 1}}} \right) $

(7)

式中:Z^k是为了简化公式而引入的序列; s^k和s^k-1依次为第k次和第k-1次迭代的结果。利用最优化算法求得最优解s^*:

$ {\mathit{\boldsymbol{s}}^*} = \mathop {\arg \min }\limits_{\mathit{\boldsymbol{X}} \in {{\bf{R}}^N}} {\left\| \mathit{\boldsymbol{s}} \right\|_1} + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{As - Y}}} \right\|_2^2 $

(8)

对于公式(7), 我们采用类似于最速下降方法来更新s:

$ \theta \left( {{\mathit{\boldsymbol{s}}^{k + 1}}} \right) = \theta \left[ {{\mathit{\boldsymbol{s}}^k} - {t_k}\nabla \theta \left( {{\mathit{\boldsymbol{s}}^k}} \right)} \right] = \mathop {\min }\limits_t \theta \left[ {{\mathit{\boldsymbol{s}}^k} - t\nabla \theta \left( {{\mathit{\boldsymbol{s}}^k}} \right)} \right] $

(9)

在满足公式(9)时, 利用下式来寻找最优步长t_k:

$ {t_k} = \frac{{1 + \sqrt {4t_{k - 1}^2 + 1} }}{2} $

(10)

由此, 公式(6)可改写为:

$ \begin{array}{l} {\mathit{\boldsymbol{s}}^{k + 1}} = {R_v}h\left( {{\mathit{\boldsymbol{Z}}^k}} \right)\\ \;\;\;\;\;\; = {\mathop{\rm sgn}} \left[ {{\mathit{\boldsymbol{Z}}^k} - \tau g\left( {{\mathit{\boldsymbol{Z}}^k}} \right)} \right] \odot \max \left\{ {\left| {{\mathit{\boldsymbol{Z}}^k} - \tau g\left( {{\mathit{\boldsymbol{Z}}^k}} \right)} \right| - } \right.\\ \;\;\;\;\;\;\left. {\frac{\tau }{\mu },0} \right\} \end{array} $

(11)

式中:⊙表示按元素逐个作乘积的算子; 阈值v=τ/μ, 其中v值越大, 算法的收敛速率越快, 因此只要使τ值尽可能大, 或使μ值尽可能小, 都能达到加快算法收敛速率的效果。

本文中采用B-B(Barzilai-Borwein)步长来更新参数τ^[22]:

$ {\tau ^k} = \frac{{{{\left\| {{\mathit{\boldsymbol{s}}^k} - {\mathit{\boldsymbol{s}}^{k - 1}}} \right\|}^2}}}{{{{\left( {{\mathit{\boldsymbol{s}}^k} - {\mathit{\boldsymbol{s}}^{k - 1}}} \right)}^{\rm{T}}}\left( {{\mathit{\boldsymbol{g}}^k} - {\mathit{\boldsymbol{g}}^{k - 1}}} \right)}} $

(12)

为了使收敛速率加快, 我们缩小μ的值。通常将δmax(A^TY)(δ>0)作为μ的初值，迭代中令μ=β·μ(0 < β < 1), 使μ逐渐减小，同时保证μ仍可平衡二次函数和L1范数在迭代过程中的比重。

在FFPC算法中, 我们采用计算两次相邻迭代结果的误差来决定是否终止迭代:

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{G}}_k}\left( {{\mathit{\boldsymbol{s}}^k},{\mathit{\boldsymbol{s}}^{k - 1}}} \right) = \frac{{{{\left\| {{\mathit{\boldsymbol{s}}^k} - {\mathit{\boldsymbol{s}}^{k - 1}}} \right\|}_2}}}{{\max \left\{ {{{\left\| {{\mathit{\boldsymbol{s}}^k}} \right\|}_2},1} \right\}}} < {s_{{\rm{tol}}}}\\ {\left\| {\mathit{\boldsymbol{g}}\left( {{s^k}} \right)} \right\|_\infty } - 1 < {g_{{\rm{tol}}}} \end{array} \right. $

(13)

式中:s_tol和g_tol为迭代终止条件, 实验中s_tol=1×e^-4, g_tol=1×e^-2。当这两个条件同时满足时, FFPC算法停止迭代。

1.3 压缩感知框架下基于FFPC算法的缺失地震数据重建步骤

压缩感知框架下基于FFPC算法的缺失地震数据重建步骤见图 1。

2 数据实验和分析

为了验证算法的有效性和实用性, 我们采用Marmousi模型合成数据和实际地震数据对算法进行验证。利用信噪比值(R_SN)来衡量缺失地震数据的重建效果:

$ {R_{{\rm{SN}}}} = 10\lg \frac{{\left\| \mathit{\boldsymbol{X}} \right\|_2^2}}{{\left\| {\mathit{\boldsymbol{\bar X}} - \mathit{\boldsymbol{X}}} \right\|_2^2}} $

(15)

式中:X为未缺失的原始地震数据; X为重建后的地震数据。

2.1 Marmousi模型合成数据重建结果分析

图 2为Marmousi模型的某炮道集合成的地震数据, 一共96道, 每道750个采样点, 时间采样率为0.001s, 我们采用图 3中的测量矩阵来对图 2地震数据进行缺失, 为了使图 2地震数据缺失40%, 选择测量矩阵大小为58×96, 该测量矩阵为0-1矩阵, 图 3中黄色亮点处值为1, 其余部分都为0。图 4为图 2随机缺失40%地震道后的地震数据。图中缺失的地震数据复杂度高, 重建难度大。

采用压缩感知方法对图 4中缺失的地震数据进行重建, 重建算法分别采用OMP算法、SPGL1算法、FPC算法和FFPC算法, 得到的重建地震数据以及与原始数据的差异如图 5~图 8所示。从重建精度上来看, 基于FPC算法和FFPC算法的重建精度最高, 噪声极少, 同相轴清晰; 基于OMP算法的重建效果最差, 噪声很多; 基于SPGL1算法的重建效果居中。在算法的耗时上, OMP算法、SPGL1算法、FPC算法、FFPC算法耗时分别为1.6661, 2.0927, 3.2708, 2.1667s。可以看出, OMP算法耗时最少, 而FFPC算法相对FPC算法, 尽管重建精度上没太多提高, 但耗时更短, 在处理一些大型数据时, 可以节省很多的时间。

采用同样的方法使图 2合成地震数据随机缺失20%~80%的地震道, 在压缩感知框架下, 分别采用OMP算法、SPGL1算法、FPC算法和FFPC算法对缺失后的地震数据进行重建。基于上述4种算法缺失地震数据重建结果的R_SN值曲线如图 9所示, 可以看出, FFPC算法和FPC算法相对OMP算法重建数据的R_SN值更高, 效果更好; 相对SPGL1算法也有优势。

2.2 实际地震数据重建结果分析

图 10为实际地震数据, 共200道, 每道900个采样点, 时间采样率0.001s。图 11是本次实验中采用的测量矩阵, 为了使地震数据缺失40%, 测量矩阵选为120×200的0-1矩阵, 黄色部分为1, 其余为0。图 12为图 10随机缺失40%地震道后的地震数据。

采用压缩感知方法对图 12中缺失的地震数据进行重建。图 13~图 16分别为基于OMP算法、SPGL1算法、FPC算法和FFPC算法得到的重建后地震数据及其与原始数据的差异。可以看出, 基于OMP算法得到的地震记录噪声最多, SPGL1算法次之, FPC算法和FFPC算法几乎不含噪声。可见, 在实际地震数据缺失重建中, FPC算法和FFPC算法得到的结果最好。实验中OMP算法、SPGL1算法、FPC算法和FFPC算法耗时分别为2.7123, 2.6997, 4.8905, 3.2709s。由此可见, 改进后的FFPC算法相对FPC算法耗时更短。

图 17为对随机缺失20%~80%地震道的实际地震数据采用不同算法重建结果的R_SN值曲线, 可以看出, 在对实际缺失地震数据进行重建时, 压缩感知理论下基于FPC算法和FFPC算法的重建结果R_SN值总是最高, 重建效果比OMP算法和SPGL1算法要好。

3 结论

FFPC算法是在FPC算法的基础上引入了软阈值步长收缩参数和正则化收缩参数, 算法的运行效率得到很大提升, 本文将FFPC算法引用到缺失地震数据的重建中, 取得了较好的实验结果。不论是对复杂的Marmousi模型合成地震数据缺失, 还是对实际地震数据缺失, 压缩感知框架下的FFPC算法都能很好的完成重建, 相对FPC算法, FFPC算法的耗时更短, 效率更高; 相对传统的OMP和SPGL1等算法, 采用FFPC算法得到的重建结果R_SN值更高, 效果更好。压缩感知理论下的数据采样方式相比传统采样方式得到的数据量更小, 效率更高, 这为以后地震数据采样方法的发展提供了很好的参考。