叠前重加权L1范数稀疏约束的地震反演方法

2023-12-12赵云文晓涛尹川韩文明李陈龙

石油地球物理勘探 2023年6期

赵云，文晓涛*，尹川，韩文明，李陈龙

（1.成都理工大学油气藏地质及开发工程重点实验室，四川成都 610059；2.成都理工大学地球物理学院，四川成都 610059；3.中海石油国际能源服务（北京）有限公司，北京 100010）

0 引言

与叠后地震反演方法相比，叠前地震反演方法可获得更多的弹性参数，能更有效地进行储层预测和流体识别[1-3]。为了提高叠前地震反演结果的准确性，可通过约束方法[4-6]将先验信息引入反演。Berkhout[7]利用反射系数的稀疏信息和初始模型的趋势信息，应用L2范数约束添加先验约束信息，提高了反演结果的准确性。然而，Zhang 等[8]发现采用L2范数约束先验信息会出现低稀疏性，导致反演结果分辨率变低，相邻岩层边界位置呈现非聚焦性平滑过渡。由于地层反射系数向量通常具有稀疏性，向量中存在稀疏非零值，因此在稀疏表示过程中如何提高先验约束信息的稀疏性非常重要。

Taylor等[9]提出了利用L1范数稀疏约束地震反演方法，随后Thueune 等[10]、Kong 等[11]、Zhang 等[12]均验证了L1范数的稀疏性优于L2范数。Guitton 等[13]优化了L1范数，提出了结合柯西约束的混合范数。李昕洁等[14]将该方法用于全波形反演，取得了较好的效果。Zhang 等[8]提出了基于L1范数优化方法的基追踪反演(BPI)算法，传统叠前AVO 反演边界尖锐问题得到一定程度的解决，同时证实了L1范数约束的稀疏性更好。张雨强等[15]、Sun 等[16]在地球物理反演中引入了自适应Lp范数稀疏表示方法，可反演出“块状”地层，反演结果比L1范数纵向分辨率更高。Pérez 等[17]提出了一种加权混合L2,1范数的叠前AVO 反演方法，优化了反演函数约束项的稀疏性，可获得更多的块状解。王治强等[18]采用了全变分(TV)约束波阻抗的方法，既可以有效压制随机噪声，还能较好地刻画地层边界。Li等[19]基于分裂Bregman 迭代算法，将L1范数与全变分正则化结合，实现了更高横向分辨率的叠前地震多道同时反演。Huang 等[20]、Wang 等[21]、耿伟恒等[22]将非凸L1-2正则化方法拓展到叠前AVO 反演，证明了L1-2范数比L1范数更稀疏。然而，上述稀疏约束忽略了岩层边界中的地震反射振幅信息，导致反演结果出现速度伪层，为后续参数预测积累了误差。

为此，本文提出了一种叠前重加权L1范数稀疏约束的地震反演方法(简称“本文方法”)，在反演目标函数稀疏项的构建中引入了基于反射系数表示的重加权矩阵，即“重加权L1范数”。相比L1范数稀疏约束(简称“传统方法”)，重加权L1范数的优越性主要体现在：其稀疏约束表示可以通过重加权矩阵补充地层边界信息。加权地层反射边界信息的引入能够提高L1范数约束的稀疏性，本文将重加权L1范数稀疏约束和初始低频模型约束同时用于构造反演目标函数，采用交替方向乘子法(ADMM)和软阈值收缩算法(ISTA)进行迭代求解[23-25]。理论模型测试和实际工区应用结果均表明，本文方法明显提高了叠前多参数同时反演的分辨率，对层速度边界刻画更清晰，极大减弱了速度伪层和密度伪层现象。

1 方法原理

1.1 正演模拟

地震信号一般可以表示为

式中：S为地震记录；ω为地震子波；R为地层反射系数；N为噪声；“*”表示褶积运算。

引入子波卷积矩阵W，将式(1)中的褶积运算转化为简单线性运算，其中W定义为

式中m为地震子波的长度。式(1)可简化为

基于Aki和Richard反射系数近似公式[26]

其中

针对式(5)，定义变量

式中：RVP、RVS、Rρ分别为纵波速度、横波速度、密度反射系数；a(θ)、b(θ)、c(θ)皆为角度系数。

近似地将RVP表示为

将MVP记作VP的自然对数，即用矩阵形式表示

根据式(7)、式(8)，将纵波的反射系数表示成矩阵形式

式中D0为单参数一阶差分矩阵，可定义为

同理可得

式中：MVS为横波速度的自然对数矩阵；Mρ为密度的自然对数矩阵。

最终推导出正演模型的矩阵表达式为

式中：B为角度系数矩阵；D为差分矩阵；M和D分别定义为

1.2 反演目标函数构建

用于反演纵波速度、横波速度和密度的自然对数与地震记录之间的关系为

式中||·||2为L2范数。

基于传统L1范数稀疏约束包含低频背景约束项，传统方法的反演目标函数为

式中：M0为低频模型矩阵；λ为低频模型约束项系数；α为L1范数稀疏约束项系数；||·||1为L1范数。

重加权矩阵数学表示形式[27]为

式中：Q为自适应重加权矩阵；qi为重加权矩阵对角线第i个元素；ri为反射系数向量R的第i个元素；N为单道采样点。矩阵Q的具体形式为

为了将重加权矩阵元素大小与反射系数建立联系，利用边界反射振幅信息将qi定义为

式中ξ为一个很小的常量，本文称为加权因子。qi与反射系数的绝对值呈负相关。当上、下两层参数变化大时，即为岩层边界时，反射系数较大，此时权重减小，使稀疏约束项的影响减小，故可以在岩层边界处产生较大反射系数；而当上、下两层参数变化小或不变时，即代表同一层内，此时反射系数小或无，将导致权重增大，间接增大了稀疏约束的比重。通过这种动态调节约束项权重，本文方法可以得到比传统方法(不重加权)更稀疏的反演结果。

参照张雨强等[15]对比范数稀疏性的方法，将本文重加权L1范数和传统L1范数的稀疏问题分别表示为

式中：H为预测值；HR是实际值；η为极小正数。

由图1 可见，交点位置即为问题的最优解，图1b 交点(蓝点)位置比图1a 更趋近坐标纵轴，这说明重加权L1范数比传统L1范数更具稀疏性，可获得更稀疏的反演结果。因此，本文构建的反演目标函数为

图1 L1范数（a）与重加权L1范数（b）的稀疏性对比

1.3 反演算法

对于传统方法的目标函数(式(17))和本文方法的目标函数(式(23))，本文采用交替方向乘子方法(ADMM)将含多个未知参数的目标函数分解成多个单一参数的子问题进行求解。

以本文方法为例，首先采用拉格朗日乘子P代替L1范数中的QBDM矩阵进行运算，将上述非线性问题转化为有约束的线性问题，得到目标函数

进一步引入二次惩罚项和对偶项[28]，将式(24)有约束的目标函数转化为无约束的目标函数，得到

于是，将式(25)拆解成分别关于M、P和C的三个子优化问题，即

①与M相关的目标函数为

②与P相关的目标函数为

③与C相关的目标函数为

针对式(26)采用函数梯度方法求解，可得

式中：“·”为标量乘法运算符号；k为迭代次数；I为单位矩阵。

针对式(27)，引入软阈值收缩算法求解同时含有L1范数和L2范数的目标函数，可得

式中sign(·)为符号函数，其定义为

针对式(28)，可根据式(26)相同的求解方法，得到

基于式(20)，本文引入的重加权矩阵对角线元素的迭代表达式为

通过对目标函数的详细求解，本文方法的地震反演流程框架如表1所示。

表1 叠前重加权L1范数的反演流程

2 模型测试

选择经典Marmousi Ⅱ模型测试本文方法和传统方法。该模型为横向500道、纵向500个采样点，其纵波速度、横波速度和密度的剖面分别如图2所示。该模型中包含多种不同厚度、不同形状的砂岩地层，其中多组薄层将作为验证本文算法的重要地层。该模型还发育边界清晰的断层和背斜构造，适合验证本文方法在刻画地层边界、减弱伪层现象方面的优越性。

图2 不同参数的理论模型

2.1 无噪测试

首先基于式(4)，利用理论模型数据计算入射角为10°、20°和30°的地层反射系数；然后分别与频率为30 Hz的雷克子波褶积合成对应的近、中、远角度的地震记录(图3)。

图3 无噪情况下不同角度数据的合成地震记录

对反演目标函数添加可以补充背景信息的低频模型，对理论模型的纵波速度、横波速度和密度数据进行高斯低通滤波，即可获得叠前纵波速度、横波速度和密度的低频模型(图4)。

图4 不同参数低频初始模型

由式(23)可知，正则化参数α用于控制稀疏约束项的权重，α越大，稀疏约束的力度越大；参数λ用于控制低频模型约束的权重，λ越大，低频模型补充的先验背景信息的影响越大。在反演迭代过程中，低频模型约束项的主要作用是保证迭代的稳定性，而不对反演效果起主导作用，因此本文引入L曲线[29]确定本文方法和传统方法的最优正则化参数α。经过大量启发式参数测试，在无噪情况下，传统方法的最优参数为λ1=5.1×10－4、λ2=5.0×10－4、λ3=5.0×10－3、α=1.5×10－8、μ=1.5×10－8；本文方法的最优参数为λ1=5.1×10－6、λ2=5.0×10－6、λ3=9.0×10－5、α=8.5×10－10、μ=1.5×10－9、ξ=1.0×10－5。

分别采用传统方法与本文方法对纵波速度、横波速度和密度进行同时反演。对比图5 与图6 可知，本文方法反演的地层边界更清晰，纵向分辨率明显更高(黑色箭头处)，尤其针对薄气砂岩储层刻画能力更强(黑色虚线圆圈)。

图5 无噪情况下传统方法不同参数的反演结果

图6 无噪情况下本文方法不同参数的反演结果

为了更准确地分析两种方法的反演效果，图7 展示了无噪情况下第100 道(单道)不同参数反演结果。由图可见，在传统方法的反演结果中，同一地层出现较强的“正弦形”扰动(红色箭头处)，导致同层中产生许多伪速度层和伪密度层(小波浪状)，反演结果不够准确；而在本文方法反演结果中，边界更清晰，很大程度地减弱了伪层现象，保证了同层反演的稳定性。传统方法反演整体呈脉冲式(黑色箭头处)，局部平滑了相邻地层边界，这是降低纵向分辨率的主要原因，而本文方法的反演结果与理论模型吻合度更高。传统方法识别薄层的分辨率低、精度低(蓝色箭头处)，而本文方法通过加权反射系数的稀疏性，对薄层的识别能力明显提高，同时也在一定程度上降低了密度反演结果的不稳定性。

2.2 抗噪性测试

为测试本文方法的抗噪性，对远、中、近角度的地震记录添加20%的高斯白噪声(图8)，再分别采用传统方法和本文方法反演。在调参数过程中，随着加噪程度的增大，建议适当增大加权因子ξ。同样，得出加噪情况下传统方法的最优参数为λ1=1.0×10－1、λ2=5.5×10－2、λ3=7.0×10－1、α=5.7×10－1、μ=2.0×10－4；本文方法的最优参数为λ1=1.9×10－2、λ2=1.5×10－2、λ3=2.4×10－1、α=1.5×10－1、μ=2.9×10－4、ξ=1.0×10－5。

对比图9与图10可知，本文方法反演的地层边界分辨率更高(黑色箭头处)；薄气砂岩储层边界更清晰，伪层现象大幅减少(黑色虚线圈处)。由此证明，在加噪情况下本文方法仍然能够获得比传统方法分辨率更高的纵波速度、横波速度和密度剖面。

图9 添加20%的高斯白噪声的传统方法不同参数反演结果

图10 添加20%的高斯白噪声的本文方法不同参数反演结果

图11展示了加噪情况下第100道(单道)反演效果对比。由图可见，相比传统方法，本文方法反演的纵波速度、横波速度和密度边界更清晰，且对薄层反演精度更高，识别更准确(黑色箭头处)。

图11 添加20%的高斯白噪声的第100 道不同参数反演结果

为了进一步测试本文方法的抗噪性，对地震记录加50%的高斯白噪声，得到的地震记录如图12所示。分别采用本文方法与传统方法反演，通过启发式确定传统方法的最优参数为λ1=1.0×10－1、λ2=5.5×10－2、λ3=7×10－1、α=5.7×10－1、μ=7.0×10－1；本文方法的最优参数为λ1=1.0×10－2、λ2=1.4×10－2、λ3=1×10－1、α=8.7×10－1、μ=5.9×10－4、ξ=1.0×10－2。

传统方法和本文方法得到的反演结果分别如图13 和图14所示。由图可见，本文方法对细薄气砂岩储层的形态刻画更清晰，识别能力更强(黑色虚线圈内)；同时，本文方法对地层边界刻画更清晰，而传统方法模糊了地层边界。

图13 添加50%的高斯白噪声的传统方法不同参数反演结果

图14 添加50%的高斯白噪声的本文方法不同参数反演结果

通过对比第100 道(单道)不同参数反演结果(图15)，本文方法对边界刻画能力更强(黑色箭头处)，反演结果更准确(蓝色箭头处)。

图15 添加50%的高斯白噪声的第100 道不同参数反演结果

综上所述，通过添加不同的噪声，本文方法仍然能够获得比较合理、可靠的反演结果。

2.3 鲁棒性测试

为测试本文算法的鲁棒性，对地震记录添加大小为地震振幅一半的异常值，异常值数量为地震信号长度的3%。此时，传统方法的最优参数为λ1=8.0×10－2、λ2=5.0×10－2、λ3=4.0×10－1、α=1.5×10－6、μ=1.5×10－4；本文方法的最优参数为λ1=1.0×10－2、λ2=1.0×10－2、λ3=9.0×10－1、α=9.5×10－4、μ=1.0×10－4、ξ=1.0×10－5。

分别采用本文方法和传统方法对第100道进行单道反演，结果如图16所示。由图可见，两种方法均受到异常值的较大影响，传统方法产生大量不准确的脉冲解，伪层现象和非聚焦平滑过渡现象加重(黑色箭头)；而本文方法表现出更强的稳定性，边界刻画更准确，这说明本文方法鲁棒性更强，反演结果更加稳定、可靠。

2.4 反演效果评价

引入信噪比(SNR)和归一化均方根误差(NRMSE)[30]用于评价反演效果。与均方根误差(RMSE)相比，NRMSE 更能反映出反演结果与理论模型之间的平均差异(或偏差)和差异的可变性，以从横向和纵向的角度更直观、准确地对比两种方法的反演效果。它们表达式分别为

式中：Y为理论参数值；为理论参数的均值；X为反演得到的参数；Ymax为理论参数的最大值；Ymin为理论参数的最小值；Tr为剖面总道数。

由表2可知，本文方法反演的纵波速度、横波速度和密度的信噪比均高于传统方法；从纵向分析，本文方法的归一化均方根误差均低于传统方法；从横向分析，本文方法与传统方法获得的纵波速度、横波速度和密度中，纵波速度、横波速度的归一化均方根误差均小于密度，说明纵波速度、横波速度与模型差异的可变性更小，反演能力更强，反演结果质量更高。

表2 无噪情况下不同方法反演结果评价

表3和表4分别为加噪20%和50%的评价结果。与表2 相似，均证明了本文方法比传统方法反演精度更高，抗噪性更强。

表3 加噪20%情况下不同方法反演结果评价

表4 加噪50%情况下不同方法反演结果的定量评价

3 实际资料应用

实际资料来源于中国西部地区，目的层储层主要为页岩。选择A 井(直井)测井数据参与实际应用。

首先，对1.770 s 处目的层按优势响应角度筛选出三组中心入射角的地震数据(图17)，并用于纵波速度、横波速度和密度的叠前三参数反演，它们依次为3°(1°～5°)、15°(13°～17°)、27°(25°～29°)。然后，对A井的纵波速度、横波速度和密度井数据进行内插、外推和低通滤波，构建出各参数的低频模型。再分别采用传统方法和本文方法进行反演，此时传统方法的最优参数为λ1=4.0×109、λ2=1.5×1010、λ3=5.0×1011、α=8.7×1013、μ=2.5×1010；本文方法的最优参数为λ1=4.0×109、λ2=5.0×1010、λ3=8.0×1011、α=8.7×1012、μ=2.5×105、ξ=5.0×10－2。

图17 过A 井的共角度叠加地震剖面

图18为过A 井的单道反演结果，可见本文方法比传统方法与实际数据吻合得更好。

图18 过A 井的单道不同方法、不同参数反演结果

同样，由图19 可见，本文方法反演结果与实际数据吻合较好(黑色箭头处)。对于纵波速度反演结果来说，传统方法受伪层干扰，导致地层分辨率低(蓝色虚线圈内)，出现混层现象；本文方法对薄层识别能力明显增强(红色箭头处)，地层的横向连续性更好(黑色虚线圈内)，对地层边界刻画更清晰(蓝色虚线圈内)，纵向分辨率明显提高。类似地，横波速度和密度反演结果也能证明以上结论，这表明本文方法在实际应用中的可行性和优越性。