三维自适应最小二乘逆时偏移技术

匡 斌 唐祥功 张 猛 赵庆国 单联瑜

(中国石化胜利油田分公司物探研究院，山东东营 257022)

1 引言

随着中国油气勘探、开发工作的不断深入，勘探重点逐步由常规油气储层向深部碳酸盐岩及非常规油气储层转移，其中保幅处理以及高精度成像是关键环节，与野外地震数据采集紧密相关，可为反演以及岩性解释提供扎实的基础数据。

双程波动方程逆时偏移技术是实现复杂构造成像的较好方法。该技术对波动方程的近似较少，具有无倾角限制的优点，可以对转换波、多次波、棱柱波等特殊波场进行成像，而且能够处理纵、横向变速问题。但利用逆时偏移技术进行复杂构造成像时还存在如下问题：①成像结果受低频噪声影响较大；②成像结果分辨率有限，受采集孔径、地震波带宽等因素的限制，无法满足复杂油气藏精细成像的技术要求；③无法提供准确的成像振幅信息，深部成像振幅弱且不均衡。上述问题制约了逆时偏移技术的广泛应用，为此需要开发反演偏移算法，以提高常规逆时偏移技术的成像精度。

传统成像技术基于互易原理对数据进行共轭转置，并不是正演算法的逆运算，可以解释为反演成像的初始迭代结果。最小二乘偏移是利用Born近似对非线性问题进行局部寻优迭代求解的一种反演方法。在最小二乘偏移理论框架下，发展了不同的实现方法：包括最小二乘Kirchhoff偏移方法、最小二乘Beam偏移方法、最小二乘单程波偏移方法和逆时偏移方法等。Kirchhoff类偏移方法已经不能满足高精度成像的要求。单程波偏移是兼顾成像质量和效率的中间方法，成像精度较双程波偏移低。随着计算机技术的发展，最小二乘逆时偏移方法将摆脱计算效率的限制，在实际数据处理中发挥重要作用。

自从Bamberger等[1]将最小二乘法引入地震反演，人们认识到反演和偏移计算存在内在联系，偏移仅仅是反演的第一步。Dong等[2]首次将逆时偏移作为最小二乘正演和偏移算子实现了最小二乘逆时偏移，证实最小二乘逆时偏移的分辨率更高，可以消除低频噪声，得到真振幅成像结果。Yao等[3]考虑了数据与模型的非线性关系，实现了非线性最小二乘逆时偏移。为了克服最小二乘逆时偏移耗时太大的问题，Dai等[4]、Godwin等[5]将相位编码技术引入最小二乘逆时偏移，通过编码超道集的形成，提高了计算效率和最小二乘偏移技术的实用能力。

随着中国深部储层的勘探、开发市场的不断增大，人们也开展了基于反演且成像精度高的最小二乘偏移方法研究[6-12]。杨其强等[13]实现了基于傅里叶有限差分算子的最小二乘偏移算法，消除了成像噪声，提高了分辨率。王彦飞等[14]、张洪宙等[15]提出一种新的偏移反演成像迭代正则化模型，并论述了求解该模型的混合型共轭梯度算法。黄建平等[16]研究了针对碳酸盐岩裂缝型储层的最小二乘偏移成像方法，取得了较好的反演效果。李振春等[17]进行了黏声介质时空域最小平方逆时偏移，在校正黏滞性吸收的同时，解决了频域方法的计算稳定性难题。胡江涛[18]提出了最小二乘逆时偏移及角度道集提取方法，重点分析了正演算子精度、梯度质量、梯度计算效率以及逆时偏移角度道集提取方法。

最小二乘逆时偏移在反演理论框架下建立目标函数，并借助伴随状态法推导梯度计算、步长搜索等基础步骤的一系列迭代反演理论算法。该技术可以克服常规逆时偏移对深部储层成像分辨率不高且振幅不均衡等问题，对非常规油气勘探及中、深部碳酸盐岩储层勘探具有重要意义。

2 最小二乘偏移原理

在常密度声介质中，地震波传播受如下声波方程控制

(1)

式中：v为准确速度场；f为震源子波；t为时间；X=X(x,y,z)为空间位置点；utotal为地震全波场，包含声介质中的所有波现象，如直达波、回转波、折射波、一次和多次反射波等。

地震勘探过程可以表达为经典的反演问题，即已知观测数据，利用波动方程投影关系估计地下介质参数场信息。如果该反问题同时用到了地震波中的所有波现象，就是全波形反演(FWI)，是一个非线性问题。为了降低非线性反演的难度，通常可以将非线性正问题进行逐步线性近似。

采用Born近似实现非线性问题的逐步线性化，假设速度场分解为背景速度v0和速度扰动δv，全波场utotal也分解为背景波场u0和散射波场us

v(X)=v0(X)+δv(X)

(2)

utotal=u0+us

(3)

将式(2)、式(3)代入式(1)得

(4)

在速度扰动较小的假设下，存在如下关系

(5)

将式(5)代入式(4)，分别得到背景波场和散射波场的传播方程[18]

(6)

(7)

背景波场传播方程(式(6))描述由炮点直接传播到检波点的波场，近似代表直达波和回转波。在散射波场传播方程(式(7))中，假设速度扰动较小，则方程右端的高阶散射项可以省略，散射波场退化为一次散射波场，近似代表一次反射波和折射波，令

(8)

m可以近似看作反射系数，则一次散射波场的传播方程为

(9)

定义L1(m)为一阶Born正演模拟算子，代表同时求解式(6)和式(9)。一次反射波数据与反射系数呈线性关系，利用反射波数据dobs进行反射系数估计，其泛函可以写为

J(m)=(dobs-L1m)T(dobs-L1m)

(10)

求解式(10)实现最小二乘偏移，以估计反射系数。可以采用梯度导引方法(如最速下降法、共轭梯度法等)构建最小二乘偏移的迭代反演算法，通过迭代不断修正偏移成像结果，达到对深部储层保幅成像的目的。

3 自适应最小二乘逆时偏移技术

基于反演算法的最小二乘逆时偏移技术近年来受到广泛关注，但是在实际应用中面临的主要问题有：①实际资料子波的不确定性；②不够准确的深度偏移速度场；③庞大的计算量。因此需要进行子波估计、数据预处理、正演和实际数据一致性处理、成像质量控制(QC)等实用化措施改善反演算法的收敛性和计算的稳定性，提高偏移成像精度并减少迭代次数。

最小二乘偏移技术的技术流程与FWI相同，在最小平方意义下使正演数据与实际采集信号的差值波场最小化，反演目标为地震反射系数。在理论上正演数据和实际采集信号的差值波场应与深度偏移结果和真实地下反射系数的差值波场对应。但是深度域速度误差不仅在正演和偏移两个环节影响计算精度，而且还造成正演模拟数据和实际采集信号之间的时移误差，从而影响差值波场品质，进而影响最小二乘反演偏移效果。

Zeng等[19]提出的自适应最小二乘逆时偏移技术是一种动态成像聚焦技术，采用自适应策略动态计算正演模拟记录和实际采集信号在波形和能量之间的最佳匹配。因为深度域速度误差对时深转换的影响经过正演、偏移两个相反的计算过程后会部分消除，因此只要增强合成数据和实际数据的一致性，就可以减少速度误差对最小二乘偏移结果的影响，从而改善偏移效果。该技术只需要在第一次迭代时计算一次时移修正量，基本不增加最小二乘偏移的总体计算量。主要技术措施包括：动态时间调整、RMS能量均衡、自适应加权因子、计算孔径自适应变化等。该技术解决了复杂构造地区聚焦收敛速度慢的技术问题，可减少梯度计算中的算子和数据噪声。

动态时间调整技术(Dynamic Time Warping，DTW)来自语音识别技术，该算法基于动态规划的思想，将测试信号与参考信号模板进行比较以获取相似度，按照某种距离测度得出两信号的相似程度并选择最佳对齐路径。DTW技术在地震勘探邻域也有广泛应用，也称为动态波形匹配(Dynamic Waveform Matching)技术[20]，可以用于对齐理论合成数据和实际数据的波形、测井合成记录的标定、P波和S波成像结果标定、NMO偏移量校正、不同炮检距传播角成像结果的一致化校正等方面。Hale[21]应用改进DTW技术实现正演数据和实际炮数据的动态波形匹配，可减少短波长速度误差对偏移结果的影响。改进DTW技术与基于局部互相关的常规DTW技术相比，具有受噪声干扰影响小、可以适应位移量剧烈波动的技术特点。定义时间域正演数据f(i)是由实际数据g(i)偏移后得到的，其中时移量s(i)=l(整数，单位为样点数)，定义误差函数

e(i,l)=[f(i)-g(i+l)]2

(11)

DTW技术采用递归方法把全局误差函数e(i,l)分解为局部子问题进行求解，使用嵌套算法简化问题，以避免搜索所有可能的路径。具体分为叠加和回溯两个步骤实现，并估算动态偏移量s(i)。叠加算法首先沿正时间顺序计算地震波传播过程中的累加误差函数d(距离)，具体实现方法由如下伪代码表示

(12)

如果把误差函数e(i,l)看成是非起伏地表面的高程，该算法的物理意义可以理解为尽量沿平缓地形寻找1，2，…，N-1之间距离最短的路径。回溯算法沿反时间顺序从前一步计算的误差函数d中提取最优动态偏移量s， 具体实现方法由如下伪代码表示

(13)

式(12)、式(13)的约束条件是相邻样点时移量的差值不能大于1， 即

|s(i)-s(i-1)|≤1

(14)

这里的1是可选择参数，可以根据需要而动态调整。使用DTW方法估算时空域多维不同数据之间的时间(或深度)偏移量是一个NP完全问题[22]，受数据信噪比、构造复杂性等因素制约，实现难度很大。

使用RMS均衡方法对正演数据与实际数据的振幅做一致性校正，以消除Born近似算法对反偏移正演数据振幅的影响。振幅校正不需要实现完全匹配，需要保护振幅的变化趋势以免影响偏移结果的保幅性。

在x-t域二维滑动时窗(h，t)下计算正演数据和实际数据的归一化互相关评估数据匹配的可信度(式(15))。自适应加权因子

(15)

图1为单炮记录的正演数据和实际数据匹配可信度加权系数，图2为使用可信度加权系数加权前、后数据的最小二乘逆时偏移效果对比。可见，使用可信度加权系数加权的最小二乘逆时偏移数据的噪声被削弱，同相轴连续性得到改善。

图1 单炮记录正演数据和实际数据匹配可信度加权系数

因为大倾角构造或者其他类似构造往往会带来较强反射能量，叠加求和后会压制周边照明弱、成像差的砂砾岩体等微幅度构造—岩性油气藏的成像效果。作为目的层的岩性油气藏构造往往倾角较小，成像不需要很大的偏移孔径。根据构造复杂性，针对目标区域使用可变的计算孔径。具体实现策略为：第一次迭代使用最大孔径以得到一个较好的初始结果并完成对大倾角构造的成像；第二次迭代先从相对较小孔径开始，然后在每次迭代、更新时逐渐加大孔径，并根据目标区的偏移效果确定迭代孔径变化量。该策略相当于对目标区进行动态控制照明，可实现针对目标区聚焦保幅成像，并减少迭代次数，节省计算量。

图2 使用可信度加权系数加权前(上)、后(下) 数据的最小二乘逆时偏移效果

4 Y地区三维实际数据测试

胜利油田通过20多年的三维地震采集工作，积累了丰富的常规及高精度三维地震资料。随着勘探目标不断向中深层、薄互层及特殊岩性体等油气藏转移，圈闭的隐蔽性不断增强。储层预测也面临着从定性到定量、从储层岩性到物性、含流体识别的挑战，精度要求越来越高。常规深度偏移成果资料成像精度、分辨率和保幅性较低，一直是制约精细储层预测技术发展的瓶颈。

最小二乘逆时偏移技术实用化所面临的难题与FWI类似，研究难度不亚于算法研发本身。相关参考文献较少，且基本局限于二维资料、海上资料。本文基于GPU加速技术，使用三维自适应最小二乘叠前逆时深度偏移技术在胜利油田Y探区等多块三维资料获得了较好的应用效果。

Y地区位于东营北部陡坡带东段、东营凹陷东部，在构造上位于坨胜永断裂带和中央隆起带东部，地下断裂系统发育，具有台地式边界断裂特征，普遍发育砂砾岩扇体，油气资源丰富。该区沙三段、沙四段时期广泛发育砂砾岩体，其油藏类型以构造—岩性油藏为主。选择Y地区6束线的单炮数据(223G)进行三维自适应最小二乘偏移测试，一共有6347炮数据，每炮1501个样点，采样率为4ms，偏移深度为6000m，深度采样间隔为10m。偏移迭代的成像面积约为100km2。使用16个K20 GPU节点并行计算，完成一次迭代使用的时间为123.7h， 共进行了5次迭代。

图3为原始单炮数据与正演数据。可见，原始单炮数据与正演数据的主要同相轴吻合较好。图4、图5分别展示了A、B测线的三维自适应最小二乘逆时偏移和常规逆时偏移效果。可见，三维自适应最小二乘逆时偏移剖面分辨率明显提高，对构造形态、断裂特征、地层接触关系及砂砾岩体的可识别能力更强。图6为常规逆时偏移与三维自适应最小二乘逆时偏移水平切片。可见，在三维自适应最小二乘逆时偏移切片上洼陷沉积区域砂砾岩边界更清楚、内幕信息更丰富。图7为常规逆时偏移、三维自适应最小二乘逆时偏移剖面(转换到时间域)。由图可见，三维自适应最小二乘逆时偏移剖面改善了复杂断块的成像质量，目的层的构造—岩性油气藏细节显示清楚。总体来说，三维自适应最小二乘逆时偏移的构造成像分辨率、精度和保幅性较常规逆时偏移有明显改善，取得了较好的地质效果。

图3 原始单炮数据(左)与正演数据(右)

图4 A测线常规逆时偏移(左)与三维自适应最小二乘逆时偏移(右)剖面 完全为纯波显示，未做修饰性处理

图5 B测线常规逆时偏移(左)与三维自适应最小二乘逆时偏移(右)剖面

图6 常规逆时偏移(左)与三维自适应最小二乘逆时偏移(右)水平切片(z=3000m，对应主要目的层)

图7 常规逆时偏移(左)、三维自适应最小二乘逆时偏移(右)过井剖面 Y551井钻遇鼻状构造，Y552、Y556井的目的层为砂砾岩体

5 结束语

作为反演偏移技术，最小二乘逆时偏移在背景速度较准确的情况下，能够反演地震波场的高波数成分，对精细储层描述有重要意义。但是也不能忽视最小二乘逆时偏移的实现难度、对信噪比和精确速度场的依赖和应用局限性。本文的研究成果证明了最小二乘逆时偏移技术在东部探区岩性勘探中具有较大应用价值。

本文得到同济大学王华忠教授、胡江涛博士，石油大学李振春教授、邓文志硕士等同仁的技术指导和帮助，在此表示感谢。

[1] Bamberger G,Chavent C,Hemon et al.Inversion of normal incidence seismograms.Geophysics,1982,47(5):757-770.

[2] Dong S,Cai J,Guo M et al.Least-squares reverse time migration: towards true amplitude imaging and improving the resolution.SEG Technical Program Expanded Abstracts,2012,31.

[3] Yao G and Jakubowicz H et al.Non-linear least-squares reverse-time migration.SEG Technical Program Expanded Abstracts,2012,31.

[4] Dai W and Schuster G T.Least-squares migration of multisource data with a deblurring filter.Geophy-sics,2011,76(5):R135-R146.

[5] Godwin J and Sava P.A comparison of shot-encoding schemes for wave-equation migration.Geophysical Prospecting,2013,61(5):391-408.

[6] 郭振波,李振春.最小平方逆时偏移真振幅成像.石油地球物理勘探,2014,49(1):113-120. Guo Zhenbo and Li Zhenchun.True amplitude imaging on least square reverse time migration.OGP,49(1):113-120．

[7] 陈云峰,王华忠,任浩然等.线性反演最小二乘叠前偏移的矩阵形式解析.石油地球物理勘探,2014,49(6):1091-1096. Chen Yunfeng,Wang Huazhong,Ren Haoran et al.Matrix representation of least square prestack migration based on linearization.OGP,49(6):1091-1096．

[8] 黄建平,孙陨松,李闯等.一种基于分频编码的最小二乘裂步偏移方法.石油地球物理勘探,2014,49(4):702-707. Huang Jianping,Sun Yunsong,Li Chuang et al.The least squares DSR migration method based on frequency selection encoding.OGP,2014,49(4):702-707.

[9] 黄建平,曹晓莉,李振春等.最小二乘逆时偏移在近地表高精度成像中的应用.石油地球物理勘探,2014,49(1):107-112. Huang Jianping,Cao Xiaoli,Li Zhenchun et al.Application of least square reverse time migration in high resolution imaging of near surface.OGP,49(1):107-112．

[10] 李振春,李闯,黄建平等.基于先验模型约束的最小二乘逆时偏移方法.石油地球物理勘探,2016,51(4):738-744. Li Zhenchun,Li Chuang,Huang Jianping et al.Regularized least square reverse time migration with prior modeling.OGP,51(4):738-744．

[11] Li Chuang,Huang Jianping,Li Zhenchun et al.Pre-conditioned prestack plane-wave least squares reverse time migration with singular spectrum constraint.

Applied Geophysics,2017,14(1):73-86.

[12] 王华忠,胡江涛,郭颂.最小二乘叠前深度偏移理论和方法.石油物探,2017,56(2):159-170. Wang Huazhong,Hu Jiangtao,Guo Song.Theory and method of least square prestack depth migration and imaging.GPP,2017,56(2):159-170.

[13] 杨其强,张叔伦.最小二乘傅立叶有限差分偏移.地球物理学进展,2008,23(2):433-437. Yang Qiqiang and Zhang Shulun.Least squares Fourier finite difference migration.Progress in Geophy-sics,2008,23(2):433-437.

[14] 王彦飞,杨长春,段秋梁.地震偏移反演成像的迭代正则化方法研究.地球物理学报,2009,52(6):1615-1624. Wang Yanfei,Yang Changchun and Duan Qiuliang.On iterative regularization methods for migration deconvolution and inversion in seismic imaging.Chinese Journal of Geophysics,52(6):1615-1624.

[15] 张洪宙,王彦飞.VTI介质迭代正则化偏移反演算法.石油地球物理勘探,2014,49(6):1083-1090. Zhang Hongzhou and Wang Yanfei.Iterative regularization migration inversion in VTI media.OGP,2014,49(6):1083-1090．

[16] 黄建平,李振春,孔雪等.碳酸盐岩裂缝型储层最小二乘偏移成像方法研究.地球物理学报,2013,56(5):1716-1725. Huang Jianping,Li Zhenchun,Kong Xue et al.A study of least-squares migration imaging method for fractured-type carbonate reservoir.Chinese Journal of Geophysics,2013,56(5):1716-1725.

[17] 李振春,郭振波,田坤.粘声介质最小平方逆时偏移.地球物理学报,2014,57(1):214-228. Li Zhenchun,Guo Zhenbo,Tian Kun et al.Least-squares reverse time migration in visco-acoustic medium.Chinese Journal of Geophysics,2014,57(1):214-228.

[18] 胡江涛.最小二乘逆时偏移及角度道集提取方法研究[学位论文].上海：同济大学，2015.

[19] Zeng Chong,Dong Shuqian,Wu Zhaohong et al.Adaptive least-squares RTM and application to Freedom WAZ subsalt imaging.SEG Technical Program Expanded Abstracts,2015,34:4059-4063.

[20] Anderson K and Gaby J.Dynamic waveform ma-tching.Information Sciences,1983,31(3):221-242.

[21] Hale D.Dynamic warping of seismic images.Geophy-sics,2013,78(2):S105-S115.

[22] Keysers D and Unger W.Elastic image matching is NP-complete.Pattern Recognition Letters,2003,24(1-3):445-453.