地震数据的缺失重建是地震数据处理的主要任务之一。目前, 地震数据缺失重建的方法有三大类:预测滤波法、波动方程法以及数学变换法。预测滤波法有F-X域预测误差滤波法^[1]以及基于加权最小二乘回归的F-X域预测滤波插值法等^[2], 该类方法不仅计算量偏大, 而且较容易引入误差。波动方程法有AMO局部叠加法^[3]以及NMO与逆DMO联合法^[4]等, 该类方法需要地下介质速度参数且计算量大。数学变换法有基于曲波变换的地震数据重建方法^[5]、基于傅里叶变换的地震数据重建方法^[6]以及字典学习地震数据重建方法^[7-8], 该类方法理论公式直观, 计算结果偏差较小, 比较适用于实际生产。

由压缩感知理论可知:倘若数据具有稀疏性, 那么只要采集少量的数据, 就可以利用合适的重建算法准确地重建原始数据^[9-10]。而地震信号可以通过数学变换进行稀疏表示, 满足压缩感知应用于地震数据重建的先决条件, 因此, 压缩感知技术在基于傅里叶变换^[11]、小波变换^[12]、曲波变换^[13]和Shearlet变换^[14]的地震数据重建中得到了较大发展。其中, 选择一个合适的稀疏基十分重要, 它直接影响到地震数据的重建效果^[15]。傅里叶变换缺乏描述地震信号局部时频特征的能力; 短时傅里叶变换能够解决局部时频问题, 但是窗口难以完全自适应; 小波变换虽然弥补了窗口难以自适应的缺陷, 但是对于二维图像中线性特征的描述能力较弱。DO等^[16]提出了Contourlet变换, 首先利用拉普拉斯金字塔结构对信号进行多尺度分解并捕获奇异点, 然后利用方向滤波器合成Contourlet域中的稀疏系数。相对小波变换和短时傅里叶变换, 它能够更好地描述地震信号中的线性特征。目前, Contourlet变换已被广泛应用于地震数据处理, 如:在Contourlet域中对地震数据的稀疏系数进行阈值处理, 从而完成地震数据重建^[17];采用Contourlet和Curvelet相结合的方法对地震数据进行稀疏表示, 再通过正则化方法对地震数据进行去噪及重建^[18];采用Contourlet变换来增强地震剖面的横向连续性^[19]。本文提出一种基于压缩感知和Contourlet变换的地震数据重建方法, 通过合成数据测试分析和实际地震数据重建验证了方法的有效性。

1 基本原理 1.1 压缩感知

根据压缩感知理论, 地震数据的缺失重建是将M维不完整的地震数据f(f∈R^M)恢复为N维完整的地震数据y(y∈R^N, N>M):

$ \mathit{\boldsymbol{f}} = \mathit{\phi} \mathit{\boldsymbol{y}} $

(1)

式中, ϕ(ϕ∈R^M×N)为测量矩阵算子。压缩感知技术针对稀疏信号, 因此将y在可变域上变稀疏:

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

(2)

式中:φ为稀疏变换算子, s为y在可变域上的稀疏系数。将(1)式和(2)式合写为:

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

(3)

式中:φ^H为稀疏逆变换算子。引入重构算子A=ϕφ^H, 将其代入(3)式中, 有如下表示:

$ \mathit{\boldsymbol{f}} = \mathit{\boldsymbol{As}} $

(4)

公式(4)中的未知量个数大于方程个数(N>M), 所以为病态方程, 无法进行求解。CANDES等^[20]指出, 在s具有稀疏性和A满足有限等距性质(restricted isometry property, RIP)的情况下, 可以对欠定方程(4)进行求解。即对于稀疏度为k的任意稀疏信号s, 存在常数ν_k∈(0, 1), 使A满足:

$ \left( {1 - {\nu _k}} \right){\left\| \mathit{\boldsymbol{s}} \right\|_2} \le {\left\| {\mathit{\boldsymbol{As}}} \right\|_2} \le \left( {1 - {\nu _k}} \right){\left\| \mathit{\boldsymbol{s}} \right\|_2} $

(5)

一般情况下, A是否满足RIP很难由(5)式判定。BARANIUK^[21]提出, 当φ^H与ϕ不相干时, A很大概率上满足RIP。因为A=ϕφ^H, φ^H与ϕ的相干性可以定义为^[22]:

$ \mu \left( \mathit{\boldsymbol{A}} \right) = \mathop {\max }\limits_{i,j,i \ne j} \frac{{\left\langle {\mathit{\boldsymbol{a}}_i^{\rm{T}},{\mathit{\boldsymbol{a}}_j}} \right\rangle }}{{\left\| {{\mathit{\boldsymbol{a}}_i}} \right\| \cdot \left\| {{\mathit{\boldsymbol{a}}_j}} \right\|}} $

(6)

式中, a_i为A的第i列向量。在地震数据采集中, 因为各个方位的检波器采集的是点震源发出的脉冲, 所以ϕ是一个单位矩阵, 它与φ^H不相干。当该地震道未采集到数据时, ϕ的a_i列元素为0, 采集到数据时, ϕ的a_i列元素为1^[23]。

此时, 求解方程(4)等价于求解下述模型:

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

(7)

式中: ${\mathit{\boldsymbol{\tilde s}}}$为近似稀疏系数, ε为重构误差。求解s的最小L₀范数是一个NP(non-deterministic polynomial)难题^[24], DONOHO^[25]将压缩感知中求解最小L₀范数问题转化为求解最小L₁范数问题:

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

(8)

(8) 式是一个凸优化问题, 可改写为:

$ \mathit{\boldsymbol{\tilde s}} = \arg \min \frac{1}{2}{\left\| {\mathit{\boldsymbol{As}} - \mathit{\boldsymbol{f}}} \right\|_2} + \lambda {\left\| \mathit{\boldsymbol{s}} \right\|_1}\;\;\;\;\lambda > 0 $

(9)

式中:λ为正则化参数。本文采用快速迭代收缩阈值算法(FISTA)求解(9)式, 得到${\mathit{\boldsymbol{\tilde s}}}$后再由φ^H得到重构的地震数据${\mathit{\boldsymbol{\tilde y}}}$, 从而完成地震数据的重建:

$ \mathit{\boldsymbol{\tilde y}} = {\mathit{\boldsymbol{\varphi }}^{\rm{H}}}\mathit{\boldsymbol{\tilde s}} $

(10)

在压缩感知技术中, 稀疏变换阶段及重建阶段的算法选择十分重要, 直接影响到重建效果。

1.2 Contourlet变换

本文在压缩感知的稀疏变换阶段选择Contourlet变换^[26]。该变换主要由拉普拉斯金字塔分解和方向滤波器组采集两部分组成。拉普拉斯金字塔分解即将待分解的地震数据y与低通滤波算子d卷积, 第l层的低频信息y_l可以表示为:

$ \begin{array}{l} {\mathit{\boldsymbol{y}}_l} = \left( {g,h} \right) = \sum\limits_{m = - 2}^2 {\sum\limits_{n = - 2}^2 {d\left( {m,n} \right)y\left( {2g + m} \right)\left( {2h + n} \right)} } \\ l \in \left( {0,F} \right]\;\;\;g \in \left( {0,{A_l}} \right]\;\;\;h \in \left( {0,{B_l}} \right] \end{array} $

(11)

式中:F为金字塔分解总层数, A_l为地震数据的行, B_l为地震数据的列。第l层的高频信息y_lh可以表示为:

$ {\mathit{\boldsymbol{y}}_{l{\rm{h}}}} = \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{y}}_l} $

(12)

方向滤波器组的采集矩阵E_p^l可以表示为:

$ \mathit{\boldsymbol{E}}_p^l = \left\{ \begin{array}{l} {\rm{diag}}\left( {{2^{l - 1}},2} \right)\;\;\;\;0 \le p \le {2^{l - 1}}\\ {\rm{diag}}\left( {2,{2^{l - 1}}} \right)\;\;\;\;{2^{l - 1}} \le p \le {2^l} \end{array} \right. $

(13)

式中:p为方向角度。方向滤波器组对每一层的y_lh进行分解, 得到相应的稀疏系数s_l, 以内积形式表示为:

$ {\mathit{\boldsymbol{s}}_l} = \left\langle {{\mathit{\boldsymbol{y}}_{l{\rm{h}}}},\mathit{\boldsymbol{\gamma }}} \right\rangle $

(14)

式中:γ为Contourlet系数空间R^Z(Z=0, 1, …, 2^l-1)的基函数。

1.3 快速迭代收缩阈值算法

本文在压缩感知的重建阶段采用快速迭代收缩阈值算法FISTA)^[27-28], 它是ISTA^[29]的一种改进算法。ISTA的迭代形式为:

$ {{\mathit{\boldsymbol{\tilde s}}}_{n + 1}} = {\rm{soft}}\left[ {{{\mathit{\boldsymbol{\tilde s}}}_n} + \frac{1}{\alpha }{\mathit{\boldsymbol{A}}^{\rm{H}}}\left( {\mathit{\boldsymbol{f}} - \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\tilde s}}}_n}} \right),\frac{\lambda }{{2\alpha }}} \right] $

(15)

式中: ${{\mathit{\boldsymbol{\tilde s}}}_n}$为迭代n次后得到的近似稀疏系数; α为控制步长的参数, 且α大于或者等于A^HA的最大特征值; A^H为A的逆; soft为阈值算子。

$ {\rm{soft}}\left( {\mathit{\boldsymbol{\tilde s}},T} \right) = \left\{ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tilde s}} - T}&{\mathit{\boldsymbol{\tilde s}} > T}\\ 0&{ - T \le \mathit{\boldsymbol{\tilde s}} \le T}\\ {\mathit{\boldsymbol{\tilde s}} + T}&{\mathit{\boldsymbol{\tilde s}} < T} \end{array}} \right. $

(16)

式中:T为阈值。FISTA和ISTA的区别在于每一步迭代时初始点的选择不同, 当迭代到n+1次时, ISTA选择${{\mathit{\boldsymbol{\tilde s}}}_n}$作为初始点, 而FISTA定义了一个新参数Ω_n+1作为初始点, 它表示前两次迭代的线性组合:

$ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{n + 1}} = {{\mathit{\boldsymbol{\tilde s}}}_n} + \frac{{{\beta _n} - 1}}{{{\beta _{n + 1}}}}\left( {{{\mathit{\boldsymbol{\tilde s}}}_n} - {{\mathit{\boldsymbol{\tilde s}}}_{n - 1}}} \right) $

(17)

式中:β_n+1是一个迭代更新的参量, 定义为:

$ {\beta _{n + 1}} = \frac{{1 + \sqrt {1 + 4\beta _n^2} }}{2} $

(18)

β_n+1的初始值为1。将(15)式中的${{\mathit{\boldsymbol{\tilde s}}}_n}$替换成Ω_n+1, 得到FISTA的迭代形式:

$ {{\mathit{\boldsymbol{\tilde s}}}_{n + 1}} = {\rm{soft}}\left[ {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{n + 1}} + \frac{1}{\alpha }{\mathit{\boldsymbol{A}}^{\rm{H}}}\left( {\mathit{\boldsymbol{f}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{n + 1}}} \right),\frac{\lambda }{{2\alpha }}} \right] $

(19)

2 合成数据测试及分析

本文采用信噪比R_SN来衡量地震数据的重建效果:

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

(20)

式中: ${\mathit{\boldsymbol{\tilde y}}}$为重构的地震数据; y为原始未缺失的地震数据。采用如图 1a所示合成炮集数据测试了本文方法的重建性能, 该炮集有90个地震道, 每道有900个采样点, 时间采样率为0.001s。图 1b为随机缺失50%地震道的合成炮集数据。

采用压缩感知技术对图 1b中随机缺失50%地震道的合成炮集数据进行重建。图 2a至图 2c分别为采用傅里叶稀疏基、小波稀疏基和Contourlet稀疏基重建的结果。重建过程中90%较小的稀疏系数被阈值滤除, 前后两次迭代的信噪比相对误差小于1%时停止迭代。

从图 2可以看出, 基于压缩感知技术的各稀疏基都可以完成合成数据的重建。其中傅里叶稀疏基的重建用时最短, 为10.003s, 但是重建结果稍差, 含有部分分布在同相轴周围的噪声, 能量较弱的同相轴湮没在噪声之中, 勉强能识别出来。小波稀疏基的重建用时为11.193s, 重建结果比傅里叶稀疏基好, 只含有微量噪声, 同相轴几乎都能识别。而Contourlet稀疏基重建的结果最好, 同相轴清晰可见, 几乎不含噪声, 重建用时为14.751s。图 3对比了本文方法不同迭代次数重建的合成炮集数据Contourlet域中最外层稀疏系数; 图 4为不同稀疏基重建的合成炮集数据采样率和信噪比关系曲线; 表 1对比了本文方法不同迭代次数重建随机缺失50%地震道的合成炮集数据信噪比。

从图 3可以看出, 随着迭代次数的增加, Contourlet域中最外层稀疏系数显示的信息逐渐增加, 对数据细节的描述越来越清晰。当迭代次数达到15次后, 最外层稀疏系数不再有明显变化。

利用压缩感知技术对随机缺失20%~70%地震道的合成炮集数据进行重建的结果均表明, 采用Contourlet稀疏基比采用傅里叶稀疏基和小波稀疏基重建的结果信噪比高。

3 实际地震数据重建及分析

采用如图 5a所示实际炮集数据测试了本文方法对实际地震数据的重建性能。该炮集有100个地震道, 每道有1000个采样点, 时间采样率为0.001s。图 5b为随机缺失50%地震道的实际炮集数据。

采用压缩感知技术对图 5b中随机缺失50%地震道的实际炮集数据进行重建。图 6a至图 6c分别为采用傅里叶稀疏基、小波稀疏基和Contourlet稀疏基重建的结果, 重建过程中90%较小的稀疏系数被阈值滤除, 前后两次迭代的信噪比相对误差小于1%时停止迭代。

从图 6可以看出, 基于压缩感知技术的各稀疏基都能完成实际地震数据的重建。其中傅里叶稀疏基的重建用时为11.848s, 重建结果一般, 含有微量噪声, 重建同相轴能量较弱, 模糊不清。小波稀疏基的重建用时为13.786s, 重建结果稍好, 所含噪声较少, 能识别出能量较强的同相轴。Contourlet稀疏基的重建用时为18.044s, 重建结果比前两者好, 几乎不含噪声, 能量很弱的同相轴也可以识别出来。图 7对比了本文方法不同迭代次数重建的实际炮集数据Contourlet域中最外层稀疏系数, 图 8为不同稀疏基重建的实际炮集数据采样率和信噪比关系曲线, 表 2对比了本文方法不同迭代次数重建随机缺失50%地震道的实际炮集数据信噪比。

从图 7可以看出, 随着迭代次数的增加, Contourlet域中最外层稀疏系数对实际地震数据的细节描述越来越丰富, 当迭代次数达到18次后, Contour-let域中最外层稀疏系数不再有明显增加。

利用压缩感知技术对随机缺失20%~70%地震道的实际地震数据进行重建的结果均表明, 采用Contourlet稀疏基比采用傅里叶稀疏基和小波稀疏基重建的结果信噪比高。

4 结束语

本文在压缩感知技术框架下引入Contourlet稀疏基, 提出了一种基于压缩感知和Contourlet变换的地震数据重建方法。

合成数据和实际地震数据测试结果表明, 在随机缺失20%~70%地震道时, 利用压缩感知技术能够完成地震数据的重建, Contourlet稀疏基相比傅里叶稀疏基和小波稀疏基具有更好的重建效果, 所增加的耗时有限, 在可接受的范围之内。

虽然本文方法重建的地震数据信噪比要比压缩感知技术常用的稀疏基高, 但还不能令人完全满意, 因此, 我们下一步的研究方向是构造重建效果更好的稀疏基。