数据驱动的块排列正则化多道叠前地震反演

王宁， 周辉， 王玲谦， 赵海波

1 东北石油大学地球科学学院， 大庆 163318 2 中国石油大学(北京)， 北京 102249 3 大庆油田勘探开发研究院， 大庆 163318

0 引言

地震反演作为地震勘探中一种重要的解释方法，将地震观测数据转化为弹性参数，可以有效去除地震观测数据中带限子波的影响，补充低频和高频信息(Tarantola, 1987; Shuey, 1985; 高刚等, 2013; 曹俊兴等, 2011).根据反演得到的弹性参数和不同的岩石物理模型，可以获得孔隙度、泥质含量等各种岩石物理参数(Karimi et al., 2010; Zong et al., 2012, 2013; 印兴耀等, 2014；白俊雨，2020).传统的叠后地震反演只能提供波阻抗剖面，不足以描述流体类型和地下精细构造(Sharma and Chopra, 2015；张岩等，2021；周单等，2022).相比之下，叠前地震数据包含了振幅随偏移距变化(AVO)的信息，且与纵波速度、横波速度和密度的反射系数有关(Aki and Richards, 1980).因此，叠前地震反演为流体预测和岩性识别提供了更丰富、可靠的信息(宗兆云等, 2012; Chen et al., 2020; Liu et al., 2020a, b; Chen et al., 2021；吴奎等，2021).

由于地震观测数据的带限性和噪声，地震反演存在严重的不适定性(Varela et al., 2006).为了获得稳定、准确的反演结果，正则化方法可以对反演结果施加特定的约束，从而获得具有期望特征的反演结果.目前有两种典型的正则化地震反演方法，包括平滑约束和块约束(VanDecar and Snieder, 1994; Ulrych and Sacchi, 2005；曹丹平等, 2010；印兴耀等, 2020a；杨俊等，2020).Tikhonov(1963)提出了L2范数来约束反演参数，可以根据初始模型的低频趋势得到平滑的反演结果.Buland和Omre(2003)将贝叶斯框架与地震反演相结合，假设反演参数服从高斯先验分布，使反演结果平滑.然而，这些平滑约束会模糊边界和断层，并不能提高反演结果的分辨率.为了加强反演结果的不连续性，目前提出了多种形式的块约束.Alemie和Sacchi(2011)引入贝叶斯框架下的三元柯西概率分布，基于稀疏的假设条件，重构地下反射系数序列，提高反演结果的分辨率.Zhang等(2013)将L1范数引入叠前反射系数重构中，通过提供准确的地下反射系数序列，增强了反演弹性参数的块状特征，更清晰地刻画了反演结果中的边缘和断层.块状约束的主要挑战和局限性是地下薄层和弱反射层的预测(Wang et al., 2019).此外，这些方法都是逐道进行的，忽略了地下构造的空间连续性.单道处理方法降低了反演剖面的信噪比，不可避免地干扰了进一步的解释和储层预测(Ma et al., 2017, 2019).

目前，有许多研究致力于多道同时反演和反射系数重构(Gholami and Sacchi, 2013; Hamid et al., 2018; Huang et al., 2021; 霍国栋等, 2017; 印兴耀等, 2020b).Gholami(2015, 2016)利用全变分约束进行多道非线性波阻抗反演，该方法可以使反演结果具有块状特征，增强其空间连续性.Hamid和Pidlisecky(2015)将所有道连接成一道，并对相邻道反演结果施加水平平滑约束.当真实情况违背事先的假设条件，这些多道反演方案不能提供准确的反演结果.因此，后续发展出了数据驱动的多道地震反演方法.Karimi(2015)利用平面波分解方法提取地下构造倾角，然后利用光滑算子沿局部斜率约束反演结果.Hamid和Pidlisecky(2016)利用观测数据计算地震记录同相轴的局部倾角，并将倾角信息和测井数据引入到地震反演中，提高反演结果的横向连续性.Huang等(2022)从地震数据中估算了地震斜率属性，并将其纳入全变分正则化中，获得高分辨率反演结果.

受数据驱动多道反演方法优势的启发，本文将Wang等(2021)提出的多道叠后地震波阻抗反演方法扩展到叠前地震反演.与以往的工作不同的是，该方法没有直接从叠前地震数据中提取构造信息.考虑到叠后地震数据具有更高的信噪比，本文利用叠后地震记录的局部相似性，提取了记录构造延展的块排列矩阵.准确的构造信息可以有效降低叠前反演的病态性，获得稳定可靠的反演结果.此外，与多道叠后地震反演不同，内存和计算效率是多道叠前地震反演的主要限制因素.本文将逐道地震记录不匹配项与多道构造约束相结合，对原目标函数进行优化，可以节省大型正演矩阵的内存，提高计算效率.

1 基本原理

1.1 块排列正则化

从叠后地震数据中提取一个块排列矩阵P，期望将反演参数重新排序为分段序列.与传统的直接平滑方法不同，该方法首先根据观测数据的局部相似性，将特征相近的采样点重新排列在一起，然后再采用约束的方法去除噪声，最后重排回原始位置.为了得到块排序矩阵，构造目标函数为：

(1)

其中D表示一阶差分矩阵，x表示待反演参数.该目标函数的求解等价于传统的TSM(Travelling Salesman)最优化问题(Cormen et al., 2001).每个地震记录块中各个采样点的振幅可以等价于TSM问题里的城市坐标，商人需要一次性通过各个城市并选取最短的路径.随机最近邻启发式算法是求解该优化问题的一种常用方法，本文采用该方法获取块排列矩阵.图2展示了块排列矩阵应用到推覆体模型后的结果.可以观察到，二维速度剖面经过块排列矩阵处理后可以转换为一维块状序列.

基于提取的块排列矩阵，可以构建具有多道构造约束的块排列正则项：

r(m)=‖QLPm‖1,(2)

(3)

可以看到，当相邻采样点存在边界时，地震记录振幅的二阶差分相对较大，即分母比较大，使相应的权重系数变小，达到减弱边界平滑的作用.

图1 地震记录块的提取(a) 二维实际资料； (b) 采用滑动窗从地震记录剖面提取的地震记录块.Fig.1 The seismic patch extraction from seismic profile(a) The 2D field seismic data; (b) The seismic patches decomposed from the seismic profile with a sliding window.

1.2 多道叠前地震反演

Zoeppritz方程可以描述振幅随偏移距变化的现象，但它具有严重的非线性，在正演模拟和反演时计算复杂(Zhi et al., 2016).本文采用Fatti近似方程(Fatti et al., 1994)进行叠前地震反演:

(4)

图2 块排列矩阵的作用(a) 合成叠后地震记录； (b) 真实纵波速度模型； (c) 真实模型所有道连接的一维速度序列； (d) 块排列矩阵作用后的一维速度序列.Fig.2 The function of patch-ordering matrix(a) The synthetic seismic data; (b) True P-wave velocity model; (c) The contaminated 1D sequence of the true model; (d) The reordered 1D velocity sequence via patch-ordering scheme.

图3 (a)Fatti近似三项系数随入射角度变化；(b)Fatti近似三项随入射角度变化Fig.3 (a) The amplitude variation with incident angle of three coefficients inFatti approximation; (b) The amplitude variation with incident angle of three terms in Fatti approximation

基于褶积模型，正演模拟可以写成矩阵的形式：

dj=WADmj,(5)

其中，dj表示第j道不同入射角的观测记录，W表示地震子波矩阵，A表示方程(4)中Fatti近似公式的系数，mj表示第j道反演的纵波阻抗和横波阻抗，其可以展开为：

(6)

考虑到横向连续性，将叠前地震反演目标函数写成多道的形式：

(7)

为了直接使用L-BFGS迭代类算法求解目标函数(Nocedal and Wright, 2006)，将目标函数中L1范数写成平滑的形式：

(8)

基于块排列正则化的多道叠前地震反演，其工作流程总结如下：

(1)利用高信噪比叠后地震记录，提取块排列矩阵，自适应地记录地下构造信息；

(2)采用井震标定的方法，将测井数据进行时深转换，并估计用于反演的地震子波；

(3)将井数据进行低通滤波，沿层位进行插值，获得低频的纵波阻抗，横波阻抗和密度的低频初始模型；

(4)使用L-BFGS算法求解多道叠前地震反演的目标函数.

2 合成数据测试

为了说明基于块排列正则化的多道叠前地震反演方法的可靠性，将该方法应用于Marmousi模型的叠前反演中，并与常规的基于模型约束的单道反演方法对比.在合成数据测试中，采用均方根误差衡量反演结果的精度：

(9)

图4a为Marmousi模型合成的叠后地震记录，使用主频为30 Hz，采样间隔为1 ms的雷克子波.图4b—d展示了添加20%随机噪声的叠前地震角道集.真实的Marmousi模型纵波阻抗和横波阻抗如图5a、b所示.初始模型通过采用81×81的平滑窗得到，如图5c、d所示.

为了便于比较，本文也提供了传统的基于模型约束的反演结果，如图6a、b所示.可以看到，传统的反演结果存在严重的噪声干扰，很难刻画地下构造信息.相比之下，提出的方法可以提供更干净准确的反演剖面.为了更直观地比较传统方法和提出方法的反演结果，对反演结果进行局部放大，如图7所示.传统反演方法得到的纵波阻抗和横波阻抗的均方根误差(Root Mean Square Error，RMSE)分别为747.0和565.4.提出方法的反演结果，其对应的RMSE分别为568.56和506.10.

在基于块排列正则化的多道叠前地震反演中，目标函数包含超参数λP和λS，本文采用试错法确定超参数，得到最优的反演结果，虽然目前有很多超参数确定的方法，但对不同的目标函数，其方法的适用性依然存在问题(张宏兵和杨长春, 2003; 张宏兵等, 2005).RMSE随超参数λP和λS的变化如图8a所示，其中叉号为纵波阻抗的RMSE，圆圈为横波阻抗的RMSE.理论上，较大的超参数会导致反演结果平滑甚至模糊，很难刻画边界和断层等细节信息.过小的超参数不能有效压制反演结果中的噪声.为了进一步验证L-BFGS算法在求解目标函数过程中的收敛性，记录了目标函数值随迭代次数的变化，如图8b所示，可以看出目标函数随着迭代次数的增加，逐渐收敛到一个稳定的值，因此求解多道反演目标函数的策略是有效的.

图4 观测地震记录(a) 叠后地震记录；角度为(b)0°； (c) 15°； (d) 30°的含噪角道集.Fig.4 The observed seismic record(a) The post-stack seismic data; The noisy angle gathers of (b) 0°, (c) 15°, and (d) 30°.

图5 真实模型和初始模型真实(a)纵波阻抗和(b)横波阻抗模型； (c) 纵波阻抗和(d)横波阻抗的初始模型.Fig.5 The true and initial modelsThe true (a) IP and (b) IS models; The initial (c) IP and (d) IS models.

图6 反演结果传统基于模型约束的(a)纵波阻抗和(b)横波阻抗反演结果；块排列正则化的(c)纵波阻抗和(d)横波阻抗反演结果.Fig.6 The inversion resultsThe conventional model-based (a) IP and (b) IS inversion results; The patch-ordering regularized (c) IP and (d) IS inversion results.

图7 局部放大的传统(a)纵波阻抗和(b)横波阻抗反演结果，块排列正则化(c)纵波阻抗和(d)横波阻抗反演结果Fig.7 Zoomed-in view of conventional (a) IP and (b) IS inversion results,patch-ordering regularized (c) IP and (d) IS inversion results

图8 (a) RMSE随超参数λP和λS的变化； (b) 目标函数值随迭代次数的变化Fig.8 (a) The RMSE variation with hyper-parameters λP and λS; (b) The value of objective function versus iterations

图9 (a) 叠后地震记录； (b) 小角度、(c)中角度和(d)大角度部分叠加角道集Fig.9 (a) The post-stack seismic data; Partially stacked angle gathers of (b) small-angle range,(c) middle-angle range, and (d) large-angle range

图10 (a)纵波阻抗和(b)横波阻抗的初始模型；传统基于模型约束的(c)纵波阻抗和(d)横波阻抗反演结果；基于块排列正则化的(e)纵波阻抗和(f)横波阻抗的反演结果Fig.10 Initial model of (a) IP and (b) IS; Model-based inversion results of (c) IP and (d) IS; Patch-ordering regularized inversion results of (e) IP and (f) IS

图11 (a)和(b) 传统纵波阻抗反演结果局部放大； (c)和(d) 块排列正则化纵波阻抗反演结果局部放大Fig.11 Zoomed-in view of (a) and (b) conventional IP inversion results,(c) and (d) patch-ordering regularized IP inversion results

图12 (a)和(b) 传统横波阻抗反演结果局部放大； (c)和(d) 块排列正则化横波阻抗反演结果局部放大Fig.12 Zoomed-in view of (a) and (b) conventional IS inversion results,(c) and (d) patch-ordering regularized IS inversion results

图13 (a)纵波阻抗和(b)横波阻抗的反演结果与测井数据比较Fig.13 Comparison between borehole data and inversion results of (a) IP and (b) IS

3 实际数据测试

将提出的方法应用于实际的地震数据中，部分叠加角道集和叠后地震记录如图9所示.小角度，中角度和大角度部分叠加角道集角度范围为0°～15°、10°～25°和20°～35°.纵波阻抗和横波阻抗的初始模型通过测井数据10～15 Hz低通滤波和插值得到，如图10a、b所示.块排列矩阵通过叠后地震记录提取，然后，利用提取的块排列矩阵构建多道反演的目标函数.基于块排列正则化的反演结果如图10e、f所示.可以看到在反演剖面上有一个大的阻抗异常体，反演结果清晰地刻画出研究区域潜山构造.为了进一步说明提出方法的优势，传统的基于模型约束的反演结果如图10c、d所示.可以看到反演的剖面信噪比较低，干扰了地下层位和构造的刻画.为了更直观地比较，将图10中反演结果的浅层和深层进行局部放大.图11a、c展示了浅层基于模型和块排列正则化反演的纵波阻抗，图11b、d展示了深部纵波阻抗反演结果.可以看到，提出的多道反演策略可以在提高分辨率的同时加强反演结果的横向连续性.此外，基于模型约束和块排列正则化的横波阻抗反演结果如图12所示，可以得到相同的结论.

为了定量评价块排列正则化反演策略的可靠性，反演结果与位于CDP 34处的测井数据进行比较，如图13所示.块排列正则化反演得到的纵波阻抗和横波阻抗(蓝色曲线)和测井曲线(红色曲线)匹配较好.它们的相关系数为0.85和0.84.黑色曲线表示基于模型约束的反演结果，可以发现，基于模型约束的反演结果虽然低频趋势与井数据基本一致，但信噪比较低，出现很多不规则抖动.基于模型约束的反演结果与测井曲线的相关系数分别为0.81和0.79.因此，提出的方法可以提供可靠的反演结果，并且能更清晰地刻画地下构造的空间展布.

4 结论

本文提出了一种基于块排列正则化的多道叠前地震反演方法.该方法利用叠后地震数据提取出一个可以记录构造方向的块排列矩阵.求取块排列矩阵的目标函数等价于传统的TSM问题，该目标函数可以用随机最近邻启发式方法有效地求解.然后利用提取的块排列矩阵构造正则项，考虑到重新排序的弹性参数序列的局部空间连续性，采用变权重系数的策略，在断层或分界面处采用较小的权重系数，在连续构造处采用较大权重系数，增强连续性.因此，这种正则化方法可以在保护边缘的同时降低反演结果中的噪声.目标函数由观测数据不匹配项、纵波阻抗和横波阻抗块排列正则项构成.在合成数据测试中，与传统的基于模型约束的反演结果相比，块排列正则化反演结果中的纵波阻抗和横波阻抗的变化更清晰、更准确.在实际资料测试中，提出方法的反演结果可以稳定、准确地描述地下地层的延展.因此，该方法可以作为今后储层预测的有效工具.此外，该方法有望应用于地震数据处理和全波形反演中，在增强信号或反演结果空间连续性的同时，保护边缘.