随着油气勘探开发的不断深入, 小尺度(薄互层)岩性油气藏的识别已成为愈发重要的问题^[1], 在褶积模型的假设下, 要求能够从地震信号中反褶积出准确的反射系数信息, 包括它的位置和振幅^[2]。有限频地震信号的反褶积常常是一个多解问题, 而在实际应用中, 地质模型的成层性决定了我们需要的是具有稀疏性质的解^[2]。压缩感知理论的信号恢复问题就是求解信号在某个基底下的稀疏表示系数, 这与地震反褶积有类似之处。VELIS^[2]将模拟退火算法和线性最小二乘法相结合, 求取反射系数序列中脉冲的时间位置和振幅, 然而模拟退火算法非常耗时, 应用于大规模数据的计算成本较高。匹配追踪(MP)算法可用于稀疏问题的求解, YANG等^[3]利用MP算法实现了压缩感知框架下的重力信号恢复, 但MP算法通常在信号超稀疏且噪声不显著时才有优势^[4]。目前广泛使用的稀疏问题求解办法是L₁正则化约束。ZHANG等^[5]认为地震反射系数可以用一系列代表楔状模型的奇偶字典线性表示, 求取表示系数时使用L₁范数进行稀疏性约束, 实现了模型约束意义下的反演。OLIVEIRA等^[6]使用L₁正则化约束实现了衰减介质中的稀疏反褶积。余慧敏等^[7]利用L₁正则化约束实现了压缩感知框架下的探地雷达成像。蒋川东等^[8]实现了基于L₁范数的低场核磁共振T₂谱稀疏反演。然而, 压缩感知理论的一些实践表明, 在稀疏恢复能力和适应噪声方面, L₁正则化不是最佳的稀疏正则化^[9], 这也启示我们必须在地震稀疏反褶积中寻找更好的稀疏约束方法。

XU等^[10-12]对非凸的L_q(0 < q < 1)正则化进行了系统研究, 揭示了L_1/2正则化的特殊重要性, 发现它有很强的稀疏约束能力, 同时具有特定求解算法。ZENG等^[13]将建立的“L_1/2正则化理论”应用于稀疏雷达成像。CHEN等^[14]将L_1/2正则化用于生物荧光层析成像, 增强了解的稀疏性和对噪声的鲁棒性。SUN等^[15]将L_1/2正则化用于遥感数据高光谱解混, 得到了优于L₁正则化的结果。姚兰^[16]将L_1/2正则化用于支持向量机分类研究中, 获得了更好的分类效果。本文将L_1/2正则化理论引入到地震稀疏反褶积, 并和L₁正则化作对比, 采用合成数据和实际数据验证它的有效性。

1 方法原理 1.1 稀疏反褶积

在自激自收条件下, ROBINSON^[17]提出利用褶积模型描述地震响应, 地震记录可近似看作是地震子波与地层反射系数的褶积结果, 即:

$ s\left( t \right) = w\left( t \right) * r\left( t \right) + n\left( t \right) $

(1)

式中:s(t)表示合成地震记录; w(t)表示地震子波; r(t)表示反射系数序列; n(t)表示随机噪声。在矩阵格式下(1)式可以表达为:

$ \mathit{\boldsymbol{s}} = \mathit{\boldsymbol{Wr}} + \mathit{\boldsymbol{n}} $

(2)

式中:W为卷积核矩阵。求取反射系数问题通常存在多解性, 然而层状地层的假设决定了反射系数r是一个稀疏脉冲序列, 于是反褶积问题可以描述为在约束条件下的稀疏求解问题。

$ \min {\left\| \mathit{\boldsymbol{r}} \right\|_0}\;\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\mathit{\boldsymbol{s}} = \mathit{\boldsymbol{Wr}} + \mathit{\boldsymbol{n}} $

(3)

式中:‖r‖₀是L₀范数, 代表向量的非零项系数的个数。类似(3)式的稀疏问题也可以从压缩感知理论导出。由压缩感知理论可知, 当信号在某个基底下能稀疏表示时(如余弦变换基、小波变换基), 可以用低于奈奎斯特采样率的方式感知信号, 且仍能高概率地恢复出信号^[18]。压缩感知的信号恢复问题可以表示为:

$ \min {\left\| \mathit{\boldsymbol{x}} \right\|_0}\;\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\mathit{\boldsymbol{y}} = \mathit{\boldsymbol{\varphi \psi x}} + \mathit{\boldsymbol{\varepsilon }} $

(4)

式中:ψ表示基底; φ表示测量矩阵; y表示测量值(维度小于原始信号维度); x表示原始信号的稀疏表示系数; ε表示噪声。可令A=φψ。(3)式和(4)式都可以直接建模为L₀正则化问题^[12], 以反褶积为例:

$ \mathop {\min }\limits_\mathit{\boldsymbol{r}} \left\{ {{{\left\| {\mathit{\boldsymbol{s}} - \mathit{\boldsymbol{Wr}}} \right\|}_2} + \lambda {{\left\| \mathit{\boldsymbol{r}} \right\|}_0}} \right\} $

(5)

式中:λ是正则化参数。我们面临的问题是, L₀正则化问题是一个NP-hard问题, 无法直接求解。已有的许多关于压缩感知和稀疏反演的文献中, 问题((5)式)通常凸松弛为L₁范数问题^[12]:

$ \mathop {\min }\limits_\mathit{\boldsymbol{r}} \left\{ {{{\left\| {\mathit{\boldsymbol{s}} - \mathit{\boldsymbol{Wr}}} \right\|}_2} + \lambda {{\left\| \mathit{\boldsymbol{r}} \right\|}_1}} \right\} $

(6)

L₁范数问题是凸优化问题, 常见算法包括梯度投影法、同伦法、迭代收缩阈值法、近似梯度法和增广拉格朗日乘子法^[19]。本文使用迭代收缩阈值算法(iterative shrinkage-thresholding algorithm, ISTA), BECK等^[20]给出了该算法的详细推证。ISTA算法源自一般化的最小化问题:

$ \mathop {\min }\limits_\mathit{\boldsymbol{x}} \left\{ {f\left( \mathit{\boldsymbol{x}} \right) + \lambda g\left( \mathit{\boldsymbol{x}} \right)} \right\} $

(7)

式中:f(x), g(x)是凸函数, 且对任意的x, y∈Rⁿ, 有‖Δf(x)-Δf(y)‖≤L(f)‖x-y‖, 其中, L(f)是Lipschitz常数, ‖·‖是Euclid范数。ISTA算法用迭代算法求解该问题, 算法的迭代格式为:

$ \begin{array}{l} {\mathit{\boldsymbol{x}}_k} = \mathop {\rm{argmin}}\limits_x \left\{ {f\left( {{\mathit{\boldsymbol{x}}_{k - 1}}} \right) + \left\langle {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{k - 1}},\nabla f\left( {{\mathit{\boldsymbol{x}}_{k - 1}}} \right)} \right\rangle + } \right.\\ \;\;\;\;\;\;\;\;\left. {\frac{1}{{2{t_k}}}{{\left\| {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{k - 1}}} \right\|}^2} + \lambda g\left( \mathit{\boldsymbol{x}} \right)} \right\}\\ \;\;\;\;\;\; = \mathop {\rm{argmin}}\limits_x \left\{ {\frac{1}{{2{t_k}}}\left\| {\mathit{\boldsymbol{x}} - \left[ {\left( {{\mathit{\boldsymbol{x}}_{k - 1}} - } \right.} \right.} \right.} \right.\\ \;\;\;\;\;\;\;\;\;\left. {{{\left. {\left. {{t_k}\nabla f\left( {{\mathit{\boldsymbol{x}}_{k - 1}}} \right)} \right]} \right\|}^2} + \lambda g\left( \mathit{\boldsymbol{x}} \right)} \right\} \end{array} $

(8)

式中:t_k为步长。对于L₁范数问题((6)式), f(x)=‖s-Wr‖₂, g(x)=‖r‖₁, (8)式进一步转化为:

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{r}}_k} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{r}} \left\{ {\frac{1}{{2{t_k}}}\left\| {\mathit{\boldsymbol{r}} - \left[ {\left( {{\mathit{\boldsymbol{r}}_{k - 1}} - 2{t_k}{\mathit{\boldsymbol{W}}^{\rm{T}}}\left( {\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{r}}_{k - 1}} - } \right.} \right.} \right.} \right.} \right.}\\ {\left. {{{\left. {\left. {\left. \mathit{\boldsymbol{s}} \right)} \right]} \right\|}^2} + \lambda {{\left\| \mathit{\boldsymbol{r}} \right\|}_1}} \right\}} \end{array} $

(9)

问题(9)由于变量可分, r_k的每一个分量可以通过求解一个简单的软阈值函数T_soft得到:

$ {\mathit{\boldsymbol{r}}_k} = {T_{{\rm{soft}}}}\left[ {\left( {{\mathit{\boldsymbol{r}}_{k - 1}} - 2{t_k}{\mathit{\boldsymbol{W}}^{\rm{T}}}\left( {\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{r}}_{k - 1}} - \mathit{\boldsymbol{s}}} \right)} \right.} \right] $

(10)

同时算法收敛要求步长t_k满足t_k∈(0, 1/‖W^T·W‖)。

尽管L₁范数约束已被广泛使用, 但压缩感知领域的研究结果表明, 使用L₁范数约束不能保证有最稀疏的结果, 也不能确保在较小采样率的情况下精确恢复信号^[12], 此外, 从地震稀疏反褶积得到的结果也能发现, L₁范数约束得到的结果在幅值上与真实反射系数的拟合度也有不足^[21]。

1.2 L_1/2正则化理论

针对L₁范数存在的问题, XU等^[10-12]系统地研究了形如(11)式的L_q正则化问题, 提出了“L_1/2正则化理论”。

$ \mathop {\rm{min}}\limits_\mathit{\boldsymbol{x}} \left\{ {{{\left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{Ax}}} \right\|}_2} + \lambda \left\| \mathit{\boldsymbol{x}} \right\|_q^q} \right\} $

(11)

XU等^[10-12]通过压缩感知研究方面的大量计算证实, q∈(0, 1/2]时, L_q正则化性能没有显著差异, q∈(1/2, 1]时, L_q正则化性能递减, 他们的工作揭示了L_1/2正则化在各种稀疏约束中的突出的代表性现象。L_1/2正则化是非凸的, 常规的凸优化工具难以求解, XU等^[12]基于严格的数学分析给出了求解L_1/2正则化问题的特定阈值迭代算法, 使得该问题能快速高效地求解。该算法的迭代格式为:

$ {\mathit{\boldsymbol{x}}_{n + 1}} = {H_\lambda }\left( {B\left( {{\mathit{\boldsymbol{x}}_n}} \right)} \right) $

(12)

其中, $B(\mathit{\boldsymbol{x}}) = \mathit{\boldsymbol{x}} + {\left({\frac{\mathit{\boldsymbol{A}}}{{\left\| \mathit{\boldsymbol{A}} \right\|}}} \right)^{\rm{T}}}\left({\frac{\mathit{\boldsymbol{y}}}{{\left\| \mathit{\boldsymbol{A}} \right\|}} - \frac{\mathit{\boldsymbol{A}}}{{\left\| \mathit{\boldsymbol{A}} \right\|}}\mathit{\boldsymbol{x}}} \right)$, H_λ(·)为阈值算子, 定义为:

$ {H_\lambda }\left( \mathit{\boldsymbol{x}} \right) = {\left( {{h_\lambda }\left( {{\mathit{\boldsymbol{x}}_1}} \right),{h_\lambda }\left( {{\mathit{\boldsymbol{x}}_2}} \right), \cdots ,{h_\lambda }\left( {{\mathit{\boldsymbol{x}}_n}} \right)} \right)^{\rm{T}}} $

(13)

其中, h_λ(·)定义为:

$ {h_\lambda }\left( z \right) = \left\{ {\begin{array}{*{20}{c}} {\frac{4}{3}z{{\cos }^2}\left\{ {\frac{{\rm{ \mathsf{ π} }}}{3} - \frac{1}{3}{{\cos }^{ - 1}}\left[ {\frac{\lambda }{8}{{\left( {\frac{{\left| z \right|}}{3}} \right)}^{ - 1.5}}} \right]} \right\}}&{\left| z \right| \ge \frac{1}{4}\sqrt[3]{{54}}{\lambda ^{2/3}}}\\ 0&{其它} \end{array}} \right. $

(14)

XU等^[12]有针对性地对比了L_1/2正则化和L₁正则化在压缩感知信号恢复问题研究方面的性能差异, 在有噪声和无噪声条件下, L_1/2正则化的恢复精度都更高, 且对欠采样的稳健性更强。这些工作为我们提供了另一种地震稀疏反褶积的思路, 即将L_1/2范数作为反射系数的稀疏约束项, 并用该特定阈值迭代算法进行求解。

2 数值实验

本文假设子波已知, 实验中选取主频为30Hz的雷克子波(图 1)。在实际处理中, 子波可以从地震记录中通过盲提取获得, 如利用高阶统计量的方法^[22]。

2.1 简单模型

在各类正则化问题中, 正则化参数λ对反演效果有显著的影响。通常依赖经验来设定λ, 为了对比在常规取值范围内、不同λ条件下L_1/2范数约束和L₁范数约束的反演效果, 构建了一个简单的反射系数序列(图 2a), 时间采样间隔为2ms, 合成对应的地震数据, 并加入15%的噪声, 结果见图 2b所示。本文使用SRE(signal-to-reconstruction error)来评价反演效果, SRE定义为:

$ {\rm{SRE}} = 10\lg \frac{{E\left( {\left\| \mathit{\boldsymbol{r}} \right\|_2^2} \right)}}{{E\left( {\left\| {\mathit{\boldsymbol{r}} - \mathit{\boldsymbol{\hat r}}} \right\|_2^2} \right)}} $

(15)

式中:$\mathit{\boldsymbol{\hat r}}$代表反演得到的反射系数序列; E(·)代表期望。SRE值越大, 反演效果越好。通过遍历不同的正则化参数并在两种范数约束下进行反褶积, 得到的评价结果如图 3所示。从图 3中可以看出, 随着正则化参数的变化, L_1/2正则化在大多数情况下优于L₁正则化。图 4展示了不同噪声水平下这两种正则化方法的性能差异(它们使用各自优选后的正则化参数), 可以看出, 噪声水平较大时, 两者的性能都有损失, 但L_1/2正则化对噪声的适应能力更强。

图 5a是构建的另一个反射系数模型, 图 5b是对应的加入了15%噪声的合成地震记录。分别用L₁和L_1/2范数约束进行反褶积, 结果如图 6a和图 6b所示, 图 7a和图 7b是对应的反演结果与真实反射系数的残差。综合图 6和图 7可以看出, 这两种稀疏约束均能得到较准确的反射系数模型形态特征, 但L_1/2范数约束的残差更小, 能更好地保护反射系数幅值。

2.2 复杂模型

利用Marmousi2 P波速度模型的局部(图 8红色方框中)构建了一个反射系数模型, 如图 9a所示。图 9b和图 9c分别为对应的不含噪声的和含20%噪声的地震记录, 共150道, 时间采样间隔为2ms。

图 10是采用不含噪声的数据反演得到的结果, 其中图 10a是采用L₁范数约束得到的结果, 图 10b是采用L_1/2范数约束得到的结果。大体来看, 采用这两种方法都能得到反射系数模型主要格架, 但从标注红色方框的部分可以看出, 采用L₁范数约束得到的结果中存在反射系数缺失的情况, 不利于横向连续性的保持和弱反射系数的保护, 而采用L_1/2范数约束能较好地保存反射系数信息, 横向连续性好, 刻画地层尖灭的能力更强。

图 11是采用含噪声的数据反演得到的结果, 其中图 11a是采用L₁范数约束得到的结果, 图 11b是采用L_1/2范数约束得到的结果。由于噪声的存在, 此时采用这两种方法反演的质量相比无噪声情形下的结果都有所下降。从红色方框指示的位置可以看到, 采用L₁范数约束得到的反射系数模型存在明显的层位缺失的现象, 有的弱反射层被整体破坏; 采用L_1/2范数约束得到的反射系数模型虽然也存在缺失或振幅削弱的地方, 但尚能保存大致轮廓, 对地层格架特征的描述依然优于L₁范数约束。

3 实际数据测试

为了验证上述数值实验结果, 对某实际二维地震数据进行了应用测试, 其叠加剖面如图 12所示。图 13是用统计性方法从实际地震记录中提取的子波。该区域薄层发育, 纵向分辨率的提升十分必要。

图 14是L_1/2正则化约束反演结果, 从标注红色方框的部分和箭头指示位置来看, L_1/2正则化约束能够较好地揭示某些薄层结构的稀疏反射系数, 有较强的刻画薄层和微小透镜体(箭头指示的部分)的能力, 也能较好地保持原始地层模型的纵向展布特征。图 15展示了反褶积前、后地震数据频谱的变化, 可以看出, 采用本文方法反褶积之后, 子波影响得到了有效消除, 地震资料有效频带宽度显著提高, 剖面的地层结构分辨率得到了明显提升。

4 结束语

本文将用于压缩感知领域的L_1/2正则化理论引入到地震稀疏反褶积中, 采用了它的特定阈值迭代算法进行求解, 形成了一种新的非凸正则化约束的稀疏反褶积方法。理论合成数据测试结果表明, L_1/2正则化的表现优于L₁正则化, 适应噪声的能力更强, 能更好地保护弱反射系数振幅。实际资料应用结果表明, L_1/2正则化突出的稀疏表达能力, 使它能够很好地刻画薄层结构和透镜体, 为地震资料高分辨率处理提供了有力的工具。当然, 作为一种数据处理手段, 无论采用何种算法, 反褶积的结果都需要其它手段验证, 如测井信息的标定, 这是实际工作中不能忽视的。

稀疏问题在地震信号处理中普遍存在, 如基于原子库的信号稀疏分解, 将信号用一系列具有特定时频特征的原子稀疏表示, 可以达到去噪的目的, 用L_1/2正则化理论求取稀疏表示系数, 有助于分离有效信息和噪声, 这也是下一步研究方向之一。