基于ATV正则化与初始模型约束的波阻抗反演*

李 曙，李 冉，金鹏程，曹继华，陈泽平，周元玲，杨 喜

(吉首大学信息科学与工程学院，湖南 吉首 416000)

1 波阻抗反演方法研究现状

地震波阻抗是一种重要的岩石属性，它与岩性、孔隙度等岩石物理参数密切相关，因此如何获取一直是勘探地震学领域的研究热点.叠后地震反演是得到波阻抗的有效方法，传统的基于Tikhonov正则化的波阻抗反演方法，通过求解L2范数最优化问题来得到波阻抗.这类方法具有简单和理论成熟的优点，但是存在反演结果易受噪声影响、垂直分界面过于光滑等问题，从而导致反演得到的波阻抗剖面的垂向分辨率不高、多解性较大等.为了克服这些问题，地球物理学家们尝试引入新的理论.近年来，以压缩感知和稀疏最优化为代表的新理论引起了他们的注意.基于压缩感知的算法，在只有少数观测数据的前提下，能够以极高的精度对稀疏(或变换域稀疏)信号进行重构.[1]基于此，大量地球物理研究人员提出了基于压缩感知理论的地震信号噪声衰减、地震反褶积、叠前弹性参数反演和波阻抗反演等方法[2-6].在波阻抗反演方面，2010年，Wang Y F[7]提出了一种基于L1范数正则化和梯度下降的方法，该方法验证了L1范数正则化在波阻抗反演方面的可行性.L1范数正则化可以提高反演结果的抗噪声性能，但是没有充分利用地下岩层的先验信息，如一般认为岩层是呈块状分布的，即同一岩层内部的岩石属性近似.为了将这一信息引入波阻抗反演以得到更符合地质先验的结果，2014年，Zhang F C等[8]提出了基于L1范数误差函数和全变分(Total Variation,TV)正则化的波阻抗反演方法，将数据保真项用L1范数表征的同时将波阻抗的全变分作为约束项，利用迭代重加权最小二乘方法求出反射系数后，借助波阻抗和反射系数之间的递推关系求出波阻抗.实验结果表明：在已知0时刻的波阻抗值情况下，该方法对高斯噪声和脉冲噪声具有较好的鲁棒性；全变分正则化的加入使得该方法能更清晰地刻画地层的边界.与Zhang F C等用L1范数表征数据保真项的做法不同，2015年，Liu C等[9]提出了用L2范数表征数据保真项并对反射系数项施加L1范数约束的波阻抗反演方法.与传统的阻尼最小二乘反演方法相比，基于L1范数正则化的波阻抗反演方法具有明显的抗噪声优势.

上述波阻抗反演方法只考虑了单道情况，没有考虑地震数据各道之间的联系.H Hamid等[10]发现了这一问题，提出了具有横向约束的多道波阻抗同时反演方法.研究结果显示，多道同时反演能充分利用地震数据的横向连续性信息，提高反演结果的横向连续性和抗噪声性能，从而揭示更多的地质结构特征.A Gholami[11]提出了一种基于全变分正则化的非线性多道波阻抗同时反演方法，在进行多道同时反演时，利用地震数据在时空信息方面的特点和全变分正则化在获得块状波阻抗方面的优势，并采用DCT变换提升计算速度.Yuan S Y等[12]提出了具有变换域稀疏提升的多道波阻抗同时反演方法，利用信号和噪声在变换域中系数值不同的特点来抑制噪声.实验结果表明，利用该方法可以稳定反演过程，提高反演结果的横向连续性，对高波数噪声有较强的抑制能力.

前面介绍的反演方法都只是联合使用1～2种正则化方法，没有充分利用先验信息.在这种背景下，Li S等[4]提出了具有多种正则化约束的地震波阻抗反演方法，综合了反射系数的L1范数稀疏约束、波阻抗的各向同性全变分约束、初始模型约束和多道同时反演等手段.实验结果显示，与基于L1范数正则化或全变分正则化的多道同时反演方法相比，该方法具有更优的抗噪声性能和更小的反演误差，且能更清晰地刻画地层的分界面.研究结果同时表明：与基于L1范数正则化多道同时反演方法相比，基于全变分正则化多道同时反演的反演效果有明显提升；与基于全变分正则化多道同时反演方法相比，基于L1范数正则化和全变分正则化的多道同时反演方法的反演效果只是略有提升，却面临更复杂的正则化参数选择问题.

上述反演方法中用到的全变分方法都是各向同性全变分(Isotropic Total Variation,ITV)方法，与ITV正则化方法相比，基于L1范数最小化的各向异性全变分正则化(Anisotropic Total Variation,ATV)方法对地震信号的稀疏性描述能力更强且实现简单.由于地球的滤波效应将导致地表接收到的地震信号的低频成分缺失，而在反演算法中加入初始模型约束可有效缓解这一问题，因此笔者将ATV正则化和初始模型约束相结合并采用多道同时反演的方式，实现基于ATV正则化和初始模型约束的多道波阻抗同时反演.

2 波阻抗反演方法

2.1 正演模型

2.2 反演方法

笔者将ATV正则化与初始模型约束相结合，提出了基于ATV正则化与初始模型约束的多道波阻抗同时反演方法(简称新方法).新方法的目标函数为数据保真项、ATV正则化项和初始模型约束项，即

(1)

等式右边分别为数据保真项、ATV正则化项和初始模型约束项.为了求解该目标函数，笔者采用split Bregman迭代方法[13]，(1)式可进一步表示为

(2)

(3)

(4)

(5)

(6)

至此，目标函数的求解已完成.为了更清晰地说明新方法的执行过程，下面给出算法的伪代码：

Input：μ，λ，α，tol，L，S，L′

While

k=k+1

End

End

3 仿真实验

为了测试新方法的效果，笔者利用Marmousi2模型[14]中的部分数据进行仿真实验，真实的波阻抗模型如图1所示.从图1可以看出，深度方向有450个采样点，距离方向有500道.利用图1所示的真实波阻抗剖面先求得反射系数剖面(图2)，然后用反射系数剖面和一个主频为40 Hz的雷克子波褶积产生合成地震记录，再向该合成记录中添加20%的高斯随机噪声，产生反演所需的合成地震剖面(图3)，最后对真实波阻抗模型进行高斯低通滤波，并将滤波后的结果作为初始波阻抗模型(图4).

图1 真实的波阻抗模型Fig. 1 True Acoustic Impedance Model

图2 反射系数剖面Fig. 2 Reflectivity Section

图3 合成地震剖面Fig. 3 Synthetic Seismic Section

图4 初始波阻抗模型Fig. 4 Initial Acoustic Impedance Model

设置参数μ=4×10-5,λ=4.5×10-3,α=0.6×10-5，反演结果如图5所示.为便于观察单道的反演情况，笔者抽取第200道的反演结果，并与相应位置的真实波阻抗曲线、初始模型波阻抗曲线作对比，结果如图6所示.

图5 反演结果Fig. 5 Inversion Result

图6 第200道的反演结果Fig. 6 Inversion Result of Trace 200

将图5的反演结果剖面与图1的真实波阻抗剖面进行对比，就可发现即使在较强噪声干扰的情况下，新方法仍然能得到较好的反演效果：不仅能清晰地勾勒出大多数层位，而且反演结果与真实值之间的相对误差也较小.图6所示的单道反演结果，从微观层面表明新方法的反演结果与真实波阻抗值之间的误差较小.

4 结语

提出了一种基于ATV正则化与初始模型约束的多道波阻抗同时反演方法，利用该方法既能得到块状的反演结果，又能缓解地震信号缺失低频信息的问题.此外，为了求解目标函数，笔者提出了一种基于split Bregman的迭代算法，并选用Marmousi 2模型中的部分数据进行测试.测试结果表明，即使在噪声干扰较大的情况下仍然能获得较好的反演结果，反演结果与真实值之间的相对误差较小.