地震反射系数反演在推断地下构造和展现地震剖面中更多细节等方面扮演着重要角色, 也得到了广泛研究^[1-4]。传统的提高分辨率手段, 如反褶积和谱白化方法, 只能通过增强或恢复地震主频带范围内的信息来获得带限的反射系数, 因而并不能从根本上拓宽地震数据的频带宽度以获得脉冲反射系数, 难以准确描述地层边界和一些特殊地质体, 如薄储层、透镜体、单砂体以及地层岩性尖灭位置等^[5-6]。基于反演的反射系数估计方法则可以从带限地震数据中重构地震频带以外的反射系数的频率成分, 获得宽频带的脉冲反射系数来提高地震数据分辨率和地震数据解释的精度^[7]。

随机稀疏反射系数反演是一种反射系数反演方法, 其思路是在地层反射系数稀疏的假设下, 将反射系数的估计过程看作是两个部分的综合体:第一部分是一个非线性问题, 即非零反射系数位置的确定问题; 第二部分是一个线性问题, 即确定位置后的非零反射系数的幅值求取问题^[8]。这种反射系数反演思路可以被看成是最严格意义上的稀疏反射系数反演^[9], 并且, 采用该思路的反演方法具有可以使用不同的初始种子进行多次运算以减少反演的不确定性的优点^[8]。获得的脉冲反射系数可以有效拓展地震数据的频谱, 对噪声和子波的敏感性较小, 能够较好地保护弱反射信号, 结合宽频带子波等手段可获得高分辨率地震数据^[8-10]。但是, 此类方法均在单道模型下实施, 且存在以下两个问题。①反演的预设参数, 包括正则化因子和非零反射系数个数, 通常根据经验人工选定, 不合理的参数选择通常会引发过拟合问题^[11]。VELIS^[9]详细分析了非零反射系数个数的选取方法和非零反射系数个数与正则化因子的相互关系, 发现非零反射系数的个数设定在很大程度上会影响反演的最终结果, 不仅会影响非零反射系数的位置, 也会影响非零反射系数的振幅估计, 建议谨慎选取非零反射系数的个数预设值, 并且在反演前进行大量实验。WANG等^[10]在频率域引入有限更新率理论, 改善了选取预设非零反射系数的稳定性, 但未进行统计性分析, 证据略显不足。SEN等^[12]将预设非零反射系数个数当作一个未知数, 在估计反射系数的同时估计非零反射系数个数, 但这种策略使得反演的未知数维度增加, 成本增大。②在反演过程中, 算法的收敛取决于地震数据和反演反射系数的正演模拟数据之间的匹配程度, 具有较大幅值的反射系数会对反演的趋势产生巨大影响。因此, 振幅较大的异常值对此类反演的结果危害很大^{[11, 13]}, 而基于单道模型的随机稀疏反射系数反演方法对异常值影响问题未给出合适的解决办法。

为了解决上述两个问题, 假设反射系数稀疏和横向连续, 引入相邻道数据作为约束条件, 提出了多道随机稀疏反射系数反演模型, 进而构建了一种多道随机稀疏反射系数反演方法。该方法通过对目标地震道数据的非零反射系数的位置和倾角信息进行同时非线性搜索, 在一定程度上解决了前文分析的第一个问题, 有效提升了反演预设参数的选择稳定性。与此同时, 由于采用了多道模型, 可以有效缓解地震数据中存在的少数异常值对反演的影响, 解决了第二个问题。作为单道随机稀疏反射系数反演方法的改进, 该方法保留了进行多次反演实现的能力和对估计出的反射系数的不确定性进行分析的能力。

1 方法原理

叠后地震数据可以描述为反射系数与地震子波的褶积^[14]。假设子波在某一时间和深度范围内是稳态的, 那么褶积模型可以写成如下形式:

$ \begin{array}{*{20}{c}} {{s^{(l)}}(t) = w(t) * {r^{(l)}}(t) + {n^{(l)}}(t)}\\ {l = 1,2, \cdots ,L\quad t = 1,2, \cdots ,N} \end{array} $

(1)

$ {r^{(l)}}(t) = \sum\limits_{k = 1}^K {a_k^{(l)}} \delta (t - \tau _k^{(l)}) $

(2)

将公式(2)代入公式(1)中, 可以获得:

$ \begin{array}{*{20}{l}} {{s^{(l)}}(t) = w(t) * {r^{(l)}}(t) + {n^{(l)}}(t)}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{k = 1}^K {a_k^{(l)}} w(t - \tau _k^{(l)} + 1) + {n^{(l)}}(t)} \end{array} $

(3)

式中:s(t)为地震数据; w(t)为地震子波; n(t)为噪声; r(t)为反射系数; K为非零反射系数个数; t为时间采样点; N为时间采样点总数; L为拟反演地震道总数; 上标l表示当前地震道的标识; a_k^(l)为第l道第k个非零反射系数的振幅; τ_k^(l)为时间延迟。多道地震数据的形式可以表示为:

$ \begin{array}{l} \begin{array}{*{20}{l}} {{s^{(l + \Delta l)}}(t) = w(t) * {r^{(l + \Delta l)}}(t) + {n^{(l + \Delta l)}}(t)}\\ {\quad = \sum\limits_{k = 1}^K {a_k^{(l + \Delta l)}} w(t - \tau _k^{(l + \Delta l)} + 1) + {n^{(l + \Delta l)}}(t)} \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \Delta l = \cdots , - 1,0,1, \cdots \end{array} $

(4)

式中:Δl为目标地震道数据和相邻地震道数据的道间隔。为了加入相邻道数据的约束, 假设a_k^(l+Δl)≈a_k^(l), 即相邻道数据的振幅变化小, 并引入斜率p_k^l, p=Δt/Δl, 其中, Δt为时间采样间隔。令

$ \tau _k^{(l + \Delta l)} = \tau _k^{(l)} + {\rm{floor}} (p_k^l \times \Delta l) $

(5)

其中, floor(x)代表向下取整函数, 即取不超过实数x的最大整数。公式(5)的作用是利用目标地震道数据中非零反射系数的位置计算出相邻地震道数据中非零反射系数的位置。公式(4)用矩阵形式表示为:

$ {\mathit{\boldsymbol{s}}^{(l + \Delta l)}} = \mathit{\boldsymbol{A}}(\tau _k^{l + \Delta l},p_k^l){\mathit{\boldsymbol{a}}^{(l)}} + {\mathit{\boldsymbol{n}}^{(l + \Delta l)}} $

(6)

式中:A为N×K矩阵, 其中的元素A_ik=w(i-τ_k^(l+Δl)+1)=w{i-[τ_k^(l)+floor(p_k^l×Δl)]+1}。

通过最小化地震数据与正演计算获得数据的残差进行最优化估计, 构建反演目标函数:

$ J = \sum\limits_{\Delta l} {{{\left\| {\mathit{\boldsymbol{A}}(\tau _k^{l + \Delta l},p_k^l){\mathit{\boldsymbol{a}}^{(l)}} - {\mathit{\boldsymbol{s}}^{(l + \Delta l)}}} \right\|}^2}} $

(7)

$ \Delta l = \cdots , - 1,0,1 $

求解过程中通过非线性搜索非零反射系数的位置和斜率, 并利用最小二乘方法获得反射系数的估计, 公式(7)对应的是常规最小二乘形式:

$ \frac{{\partial {{\left\| {\mathit{\boldsymbol{A}}(\tau _k^l,p_k^l){{\mathit{\boldsymbol{\hat a}}}^{(l)}} - {\mathit{\boldsymbol{s}}^{(l)}}} \right\|}^2}}}{{\partial \mathit{\boldsymbol{\hat a}}}} = 0 $

(8)

由(8)式可以获得:

$ \mathit{\boldsymbol{\hat a}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{s}} $

(9)

由于地震数据受噪声影响, 为使反演结果更加稳定, 反演时需要增加正则化约束。本文采用了二次型正则化约束^[8], 则公式(9)改写为:

$ \mathit{\boldsymbol{\hat a}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} + \lambda \mathit{\boldsymbol{I}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{s}} $

(10)

与公式(9)对应的、加入了正则化约束后的最终目标函数可以写为:

$ \begin{array}{*{20}{c}} {{J_a} = \sum\limits_{\Delta l} {{{\left\| {\mathit{\boldsymbol{A}}(\tau _k^{l + \Delta l},p_k^l){\mathit{\boldsymbol{a}}^{(l)}} - {\mathit{\boldsymbol{s}}^{(l + \Delta l)}}} \right\|}^2} + \lambda {{\left\| {{\mathit{\boldsymbol{a}}^{(l)}}} \right\|}^2}} }\\ {\Delta l = \cdots , - 1,0,1, \cdots } \end{array} $

(11)

求解目标函数是一个迭代的过程, 在每一次迭代中, 利用模拟退火方法^[10]获得非零反射系数的位置和斜率, 进而通过最小二乘方法计算获得反射系数, 当地震数据残差收敛或迭代次数达到预期后停止迭代^{[8, 11]}, 其中淬火温度等参数的选择, 本文参照了王万里等^[8]、VELIS^[9]和韩璇颖等^[15]给出的策略。相邻道数据的介入作为一种横向约束, 对斜率的搜索可视为是对目标地震道数据和相邻道数据之间相关性的运用。同时搜索位置和斜率会影响计算效率, 但可以通过给出斜率的先验估计减小斜率的搜索范围, 使得计算效率达到可以接受的范围, 斜率的先验估计可以通过地层倾角估计方法获得^[16-17]。同时, 虽然采用反射波同相轴横向变化不剧烈的假设(a_k^(l+Δl)≈a_k^(l)), 但该假设对反演的反射系数幅值不会产生过多影响, 这是因为每一个目标地震道数据的反射系数振幅主要由其对应的地震数据的振幅来决定, 并且, 目标地震道数据的反演是逐道进行的。在利用模拟退火方法获取非零反射系数的位置和斜率的过程中, 可以进行多次独立实验, 获取的多次实现结果可以进行不确定性分析或利用多次实现的均值降低反演结果的不确定性^{[5, 8-9]}。多道反射系数反演方法的技术流程具体如下。

1) 选定预设非零反射系数的数值、同时反演道数、最大迭代次数和终止迭代阈值。

2) 初始化时间延迟和斜率, 计算当前时间延迟和斜率下的相邻地震道数据的时间延迟。

3) 在模拟退火算法的每一次迭代中进行以下计算:①更新目标地震道数据的时间延迟和斜率, 计算更新的时间延迟和斜率下的相邻地震道数据的时间延迟; ②利用公式(10)更新A和振幅; ③利用公式(11)计算J_a; ④通过模拟退火算法中的接受准则, 决定是否接受更新后的时间延迟和斜率; ⑤利用公式(7)计算J, 如果J大于终止迭代阈值或迭代次数还未达到最大迭代次数, 回到步骤①。

4) 在步骤3)收敛后, 利用时间延迟和振幅, 根据公式(1)重构反射系数序列。

2 模拟数据测试

困扰单道随机稀疏反演方法的问题之一是在反演前需要预设非零反射系数的个数, 而预设的个数只能通过人为估计和预测来确定, 并且需要进行反复试验才可获得该个数的大概估计。而这个参数的选取会直接影响反演结果, 并且会影响正则化因子的作用效果。一个相对大的预设参数会降低反演的分辨率, 并出现高频的毛刺状的伪反射系数。本文提出的多道随机稀疏反射系数反演可以有效缓解这个问题。为了对多次独立反演结果进行统计分析, 排除其它干扰, 本文设计了一个简单的反射系数模型, 该模型包含10道地震数据, 每一道中含有3个非零反射系数, 如图 1a所示。将反射系数与采样间隔为1ms、主频为45Hz的雷克子波褶积, 并加入随机噪声, 生成了信噪比为5的模拟数据, 如图 1b所示。随机设定模拟退火算法中的初始种子值, 进行了50次独立的反演试验, 并统计分析每一次的反演结果中的非零反射系数总数。很显然, 在理想情况下, 每一次的反演结果中非零反射系数的理论值应该是3×10=30个。然而, 由于地震数据中存在随机噪声干扰, 且预设的非零反射系数个数大于真实值, 使得反演结果中的非零反射系数个数必定大于理论值。图 2为50次独立反演获得的非零反射系数个数的统计结果。在不同的预设K条件下, 对比单道随机稀疏反射系数反演和多道随机稀疏反射系数反演的结果发现, 多道随机稀疏反射系数反演估计的反射系数个数更加接近理论值, 证明多道随机稀疏反射系数反演的压制噪声和抑制高频的毛刺状反射系数的能力比单道随机稀疏反射系数反演的更强。

在图 1所示的模拟含噪数据中加入异常值干扰, 得到另一种模拟数据(图 3a)。常规单道随机稀疏反演方法力求用最少的反射系数个数达到最好的数据匹配效果, 因此振幅较大的异常值往往会被视为有效值而保留下来, 对反演结果造成负面影响, 如图 3c所示。而图 3b中的多道随机稀疏反演获得的结果没有受奇异值影响, 较好地还原了模型反射系数特征, 展示了本文方法可以降低异常值干扰的特点。图 4给出了一个稍复杂并含有界面尖灭的模拟数据试验结果, 选取的同时反演道个数为3, 即Δl∈[-1, 0, 1], 预设非零反射系数个数参数K=12, 正则化参数λ=0.1。对比单道随机稀疏反射系数反演和多道随机稀疏反射系数反演结果中红框区域可以看出, 多道随机稀疏反射系数反演结果保留了识别较薄界面的能力, 且对噪声的抵抗能力更强, 反演结果更加稳定。

3 实际数据应用

实际数据来自川东北某工区, 该区地下构造较为复杂, 地层横向变化快, 地层构造和起伏的描述对后续解释工作十分重要。反演所用子波来自井震标定结果。选取同时反演道个数为3, 即Δl∈[-1, 0, 1], 预设非零反射系数个数参数K=20, 正则化参数λ=0.25。地震数据如图 5a所示。图 5b和图 5c分别为多道随机稀疏反射系数反演和单道随机稀疏反射系数反演得到的高分辨率地震剖面, 获取过程是利用宽频带子波与反演反射系数结果进行褶积。对比图 5b和图 5c可以看出, 多道随机稀疏反射系数反演结果中断层和构造的变化展现得更加清晰, 振幅较强同相轴间的信息得以展示, 且剖面背景干净, 同相轴叠置关系清晰, 证实了多道随机稀疏反射系数反演在体现界面细节、断层和反演结果稳定程度上均具有一定的优势。实际数据的应用体现了本文方法在实际数据中的应用潜力和优势及其对改善地震解释精度的意义。

4 结论

本文提出了一种多道随机稀疏反射系数反演方法, 通过同时非线性搜索反射系数位置和斜率, 将多道反射系数褶积模型和单道随机反射系数反演方法相结合, 该方法可以用于地震数据的高分辨率处理和反演。多道随机稀疏反射系数反演方法保留了单道随机稀疏反演方法识别薄层的能力和进行不确定性分析的优势, 并提高了选取预设参数的稳定性和抵抗奇异值噪声的能力。模拟数据测试结果证明了本文方法的有效性和对单道随机稀疏反射系数反演的改进, 其中的统计性试验证明了方法对参数选择稳定性的提高。实际地震数据的应用揭示了本文方法提高地震资料分辨率和解释效果的潜力。本文方法仍存在进一步研究的空间, 如何加入测井资料约束是下一步研究的方向。