石油地球物理勘探  2019, Vol. 54 Issue (5): 1075-1083  DOI: 10.13810/j.cnki.issn.1000-7210.2019.05.015
0
文章快速检索     高级检索

引用本文 

李庆洋. T分布起伏地表最小二乘逆时偏移. 石油地球物理勘探, 2019, 54(5): 1075-1083. DOI: 10.13810/j.cnki.issn.1000-7210.2019.05.015.
LI Qingyang. Least-squares reverse time migration on undulated surface based on T-distribution. Oil Geophysical Prospecting, 2019, 54(5): 1075-1083. DOI: 10.13810/j.cnki.issn.1000-7210.2019.05.015.

本项研究受国家科技重大专项“东濮凹陷油气富集规律与!储领域”(2016ZX05006-004)、中国博士后科学基金面上项目“变换域最小二乘逆时偏移成像理论与实用化研究”(2019M652601)和中国石化中原油田分公司科技攻关项目“基于GPU加速的伪深度域最小二乘逆时偏移成像方法研究”(2019111)联合资助

作者简介

李庆洋  在站博士后, 助理研究员, 1988年生; 2011年获中国石油大学(华东)地球物理学专业学士学位, 2017年获该校地质资源与地质工程专业博士学位; 现为SEG和EAGE协会会员, 就职于中国石化中原油田物探研究院, 主要从事复杂介质地震波正演模拟与反演方法研究

李庆洋, 河南省濮阳市中原路庆山街3号中原油田物探研究院, 457001。Email: liqingyang.1988@163.com

文章历史

本文于2019年1月30日收到,最终修改稿于同年6月13日收到
T分布起伏地表最小二乘逆时偏移
李庆洋     
中国石化中原油田物探研究院, 河南濮阳 457001
摘要:在复杂构造、复杂岩性、复杂地表条件等地区进行油气地震勘探面临诸多难题,其中剧烈起伏的地表已成为制约上述地区地震勘探发展的瓶颈之一。最小二乘逆时偏移(Least-Squares Reverse Time Migration,LSRTM)相对于常规偏移具有更高的成像分辨率、振幅保真性及均衡性等优势,但基于矩形网格的LSRTM在面对复杂地表时无法很好地适用山前带等剧烈起伏地形。此外,山前带地震资料中包含较多的噪声,T分布相比Huber范数和混合模,在缺失数据的条件下更稳健,且没有多余参数,因而简单实用,而Huber范数和混合模的结果严重依赖参数选取,需要大量的尝试。为此,将全交错网格引入贴体网格,将T分布推广到起伏地表LSRTM,进一步推导了贴体网格线性Born正演方程,在此基础上提出了基于贴体全交错网格的起伏地表LSRTM算法,较好地克服了起伏地形的影响。模型试算验证了算法的有效性和对复杂模型的适应性。
关键词起伏地表    全交错网格    最小二乘逆时偏移    T分布    目标泛函    
Least-squares reverse time migration on undulated surface based on T-distribution
LI Qingyang     
Geophysical Exploration Research Institute, Zhongyuan Oilfield Branch Co., SINOPEC, Puyang, Henan 457001, China
Abstract: The hydrocarbon seismic exploration in areas with complex structure, complex lithology, and complex surface conditions has many difficulties, among which the severely undulating surface is one of the bottlenecks. The least-squares reverse time migration (LSRTM) has higher imaging resolution, amplitude fidelity, and balance than conventional migration methods. However, the LSRTM based on rectangular grid cannot deal with complex surface and overcome severe undulating terrains such as piedmont zones. In addition, seismic data in these areas often contain a lot of noise. Compared with Huber norm and mixed norm, the T distribution has a better performance in incomplete data, and it has no redundant parameters. So it is simple and practical. The results of Huber norm and mixed norm methods depend heavily on the parameter selection, which requires a lot of attempts. Therefore, a fully-staggered grid is introduced into the body-fitted grid, and the T distribution is extended to LSRTM with undulating surfaces. The linear Born forward equation of body-fitted is further deduced. On this basis, the LSRTM algorithm for undulating surface based on body-fitted fully staggered grid is proposed, which can better solve influences of undulating surfaces. Model tests verify the effectiveness and complex-model suitability of the proposed algorithm.
Keywords: irregular topography    fully-staggeredgrid    least-squares reverse time migration (LSRTM)    T distribution    object function    
0 引言

随着油气勘探、开发的深入,勘探难度逐渐加大,在复杂构造、复杂岩性、复杂地表条件等地区进行油气地震勘探面临诸多难题,其中剧烈起伏的地表已成为制约上述地区地震勘探发展的瓶颈之一[1]。传统的偏移成像方法大多基于水平地表假设,对起伏地形的适应性不强,容易导致复杂地表条件下的地下成像质量不高[2]

相比于其他偏移方法,逆时偏移(RTM)应用双程波动方程延拓波场,避免了对波动方程的近似,无倾角限制,是目前较为精确的偏移成像方法,但仍属于常规偏移的范畴,其偏移算子是线性正演算子的共轭转置,而不是它的逆[3]。当地震数据采集不足或不规则、地下构造复杂、波场带宽有限时,常规偏移方法只能对地下构造模糊成像,无法满足岩性油气藏勘探、开发的需求。最小二乘偏移(Least-Squares Migration, LSM)将成像看作最小二乘意义下的反演问题,通过不断拟合误差泛函,得到振幅保真性更好、分辨率更高、偏移噪声更少的成像结果[4-6]。早期的LSM主要基于射线理论[7-8]和单程波理论[9-13],近年来基于双程波理论的LSRTM得到广泛关注[14-19],但目前的LSRTM算法大多基于水平地表假设,没有考虑起伏地形的影响,在一定程度上限制了其推广应用。

正演模拟是LSRTM的组成单元,有限差分法计算效率高、占用内存小且易简单实现,是当前应用最广泛的求解双程波动方程方法,但在模拟复杂地表及地下复杂界面时,在应用笛卡尔坐标系下的规则网格剖分时,易出现阶梯状边界而形成虚假绕射波。贴体网格能很好地拟合任意起伏地形,宜于处理起伏地表问题[20]。但贴体网格需要坐标变换,而求解变换后的曲坐标系波动方程时,标准交错网格会引入插值误差而降低模拟精度。为此,本文将全交错网格[21]引入贴体网格中,实现了起伏地表有限差分正演模拟,进一步推导了贴体网格线性Born正演方程,在此基础上提出了基于贴体全交错网格的起伏地表LSRTM算法,较好地克服了起伏地形的影响。

此外,起伏地形条件的地震数据常含有较多噪声,常用的基于L2模拟合的LSRTM算法对噪声非常敏感,尤其是当地震数据中含有奇异值时,常规LSRTM结果被噪声严重干扰。L1模对强噪声的容忍性较L2模更好,但由于L1模在零值处不可导,因而常用Huber范数或L1、L2混合模替代[22]。T分布是另一种较L2模更稳健的方法,已被成功用于噪声数据全波形反演[23]。Aravkin等[24]认为,T分布相比Huber范数和混合模,在缺失数据的条件下更稳健,且没有多余参数,因而简单实用,而Huber范数和混合模的结果严重依赖参数选取,需要大量的尝试。本文将T分布推广到起伏地表LSRTM,通过处理含噪数据,证实了方法的有效性。同时考虑到LSRTM的计算量,引入动态相位编码技术,将多炮数据组合成一个超道集,明显节省了计算量、提高了计算效率。

1 贴体网格正演模拟

众所周知,RTM及LSRTM算法的核心为波场延拓,波场模拟方法的优劣直接决定LSRTM的成败,在起伏地形条件下更是如此。因此,在起伏地形条件下开展LSRTM,首先要研究起伏地形条件的正演模拟。传统有限差分法在面对起伏地形时存在困难,而贴体网格则是一种很好的解决方案,在山前带具有广阔的应用前景。

贴体网格是一种适合复杂地表介质的网格离散方法,网格生成的原则是使离散后的网格边界与地表形态吻合,以避免人为产生的阶梯边界引起的虚假散射,贴体网格可以由计算空间到物理空间的坐标变换获得(图 1)。

图 1 贴体网格映射示意图 左侧笛卡尔坐标系下的弯曲网格为物理空间,右侧曲线坐标系下的规则方形网格为计算空间

贴体网格生成之后,笛卡尔坐标系(x, z)与曲线坐标系(ξ, η)的网格点建立了一一对应关系,从而可求得两个坐标系间的映射关系(即雅克比矩阵)ξx, ξz, ηx, ηz,且有$ {f_x} \equiv \frac{{\partial f}}{{\partial x}}, {f_z} \equiv \frac{{\partial f}}{{\partial z}}$f代表ξη

利用链式法则和常规笛卡尔坐标系波动方程,可导出曲线坐标系二维声波波动方程[25]

$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial {v_x}}}{{\partial t}} = {\xi _x}\frac{{\partial p}}{{\partial \xi }} + {\eta _x}\frac{{\partial p}}{{\partial \eta }}}\\ {\rho \frac{{\partial {v_z}}}{{\partial t}} = {\xi _z}\frac{{\partial p}}{{\partial \xi }} + {\eta _z}\frac{{\partial p}}{{\partial \eta }}}\\ {\frac{{\partial p}}{{\partial t}} = \rho {v^2}\left( {{\xi _x}\frac{{\partial {v_x}}}{{\partial \xi }} + {\xi _z}\frac{{\partial {v_z}}}{{\partial \xi }} + {\eta _x}\frac{{\partial {v_x}}}{{\partial \eta }} + {\eta _z}\frac{{\partial {v_z}}}{{\partial \eta }}} \right)} \end{array}} \right. $ (1)

式中: vxvz为质点振动速度;p为声压;ρ为介质密度;v为声波速度。

由式(1)可知,曲坐标系下每个变量都要在同一网格点的两个方向计算空间导数,常规的标准交错网格已不满足需求(图 2a),计算vx的垂向导数与vz的横向导数均需复杂的波场插值,降低了模拟精度。由于同位网格的一阶中心差分缺少耗散特性,引起高频振荡现象,因此需要人为滤波,增加了实现的复杂性、降低了模拟精度。全交错网格机制的主要思想是速度和声压的不同分量交错定义在同一网格点(图 2b)。目前全交错网格主要用于水平地表和矩形规则网格,本文将全交错网格引入曲线坐标系[25]。由全交错网格的网格定义机制可知,其满足式(1)的交错分布特性,避免了插值和滤波运算,提高了模拟精度、降低了实现复杂性。由于同一变量定义在同一网格的两个不同位置,所以需要分别更新、计算,具体差分格式与常规标准交错网格相同,因而简单、易实现。

图 2 网格剖分示意图 (a)标准交错网格;(b)全交错网格

边界条件是正演模拟的关键,受计算机内存和计算速度等因素限制,数值模拟的范围必须是有限的,除上边界外的其他三个边界都会产生人工边界问题。为此,本文采用完全匹配层吸收边界[26],使波能够自由地穿过边界,而不产生反射;上边界为自由边界,本文采用牵引力镜像法实施[20]

2 起伏地表LSRTM基本原理 2.1 起伏地表线性化波动方程

为了推导的简便,仿照一阶速度—应力方程与二阶弹性波方程的相互转换方式,可将式(1)重新整理为

$ \begin{array}{l} {s^2}\frac{{{\partial ^2}p}}{{\partial {t^2}}} = \left( {\xi _x^2 + \xi _z^2} \right)\frac{{{\partial ^2}p}}{{\partial {\xi ^2}}} + 2\left( {{\xi _x}{\eta _x} + {\xi _z}{\eta _z}} \right)\frac{{{\partial ^2}p}}{{\partial \xi \partial \eta }} + \\ \;\;\;\;\;\;\left( {\eta _x^2 + \eta _z^2} \right)\frac{{{\partial ^2}p}}{{\partial {\eta ^2}}} + \left( {{\xi _{xx}} + {\xi _{zz}}} \right)\frac{{\partial p}}{{\partial \xi }} + \left( {{\eta _{xx}} + {\eta _{zz}}} \right)\frac{{\partial p}}{{\partial \eta }} \end{array} $ (2)

式中:s为慢度场;$ {f_{xx}} \equiv \frac{{{\partial ^2}f}}{{\partial {x^2}}}, {f_{zz}} \equiv \frac{{{\partial ^2}f}}{{\partial {z^2}}}$f代表ξη

将慢度场的平方s2分解为背景慢度场平方s02与扰动慢度场平方Δs2的叠加,即

$ {s^2} = s_0^2 + \Delta {s^2} $ (3)

由场的叠加原理可知,总波场p可理解为由背景介质产生的背景波场p0和由扰动介质产生的扰动波场ps叠加而成,即

$ p = {p_0} + {p_s} $ (4)

p0p都满足波动方程,分别将其代入式(2)并相减,然后应用Born近似,由p0替代p0+ps,可得ps的控制方程

$ \begin{array}{l} s_0^2\frac{{{\partial ^2}{p_s}}}{{\partial {t^2}}} = \left( {\xi _x^2 + \xi _z^2} \right)\frac{{{\partial ^2}{p_s}}}{{\partial {\xi ^2}}} + 2\left( {{\xi _x}{\eta _x} + {\xi _z}{\eta _z}} \right)\frac{{{\partial ^2}{p_s}}}{{\partial \xi \partial \eta }} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {\eta _x^2 + \eta _z^2} \right)\frac{{{\partial ^2}{p_s}}}{{\partial {\eta ^2}}} + \left( {{\xi _{xx}} + {\xi _{zz}}} \right)\frac{{\partial {p_s}}}{{\partial \xi }} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {{\eta _{xx}} + {\eta _{zz}}} \right)\frac{{\partial {p_s}}}{{\partial \eta }} - \Delta {s^2}\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}} \end{array} $ (5)

式(5)即为曲线坐标系线性化Born正演(反偏移)方程。可见,ps是由p0与Δs2的相互作用作为二次震源在背景介质中传播的波场,与Δs2呈线性关系,具有明确的物理含义。

定义背景介质中的格林函数G满足

$ \begin{array}{l} s_0^2\frac{{{\partial ^2}G}}{{\partial {t^2}}} = \left( {\xi _x^2 + \xi _z^2} \right)\frac{{{\partial ^2}G}}{{\partial {\xi ^2}}} + 2\left( {{\xi _x}{\eta _x} + {\xi _z}{\eta _z}} \right)\frac{{{\partial ^2}G}}{{\partial \xi \partial \eta }} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left( {\eta _x^2 + \eta _z^2} \right)\frac{{{\partial ^2}G}}{{\partial {\eta ^2}}} + \left( {{\xi _{xx}} + {\xi _{zz}}} \right)\frac{{\partial G}}{{\partial \xi }} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left( {{\eta _{xx}} + {\eta _x}} \right)\frac{{\partial G}}{{\partial \eta }} + \delta \left( {x - {x_s}} \right) \end{array} $ (6)

利用格林函数可将式(5)中的ps表示为

$ {p_s} = - \int_\Omega {\rm{d}} x\int_t {\rm{d}} tG\Delta {s^2}\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}} = \int_\Omega {\rm{d}} x\int_t {\rm{d}} tG\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}}m $ (7)

式中:m=-Δs2,将其定义为模型参数;Ω为积分空间范围。为方便后续的推导,式(7)可写成算子的形式

$ {p_s} = \mathit{\boldsymbol{Lm}} $ (8)

式中$ \mathit{\boldsymbol{L}}{\rm{ = }}\int_\Omega {{\rm{d}}\mathit{x}} \int_t {{\rm{d}}\mathit{tG}\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}}} $为起伏地表扰动波场的线性正演算子,形式上与水平地表线性正演算子相同,不同之处在于格林函数G是由起伏地表正演方程(式(6))计算获得,考虑了起伏地形的波场传播特征。

2.2 基于T分布的起伏地表LSRTM

基于反演的成像方法寻求最优的地下介质模型,以使正演波场与观测波场残差的模最小,是一个最小范数问题。常规基于L2模的LSRTM目标泛函为

$ {J_{{{\rm{L}}_2}}}(\mathit{\boldsymbol{m}}) = \frac{1}{2}{\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right)^2} $ (9)

式中pobs为观测记录。

L1模相比L2模更稳健,但当数据误差接近于零时,应用L1范数准则求取目标函数梯度会出现不稳定,因而常用Huber范数或混合模代替,但由于最后的反演结果严重依赖于参数的选取,一般需要多次尝试,大大增加了计算量[27]

与混合模不同的是,T分布对参数的依赖性较弱,因而更简单、高效。基于T分布的LSRTM目标泛函为

$ {J_{\rm{T}}}(\mathit{\boldsymbol{m}}) = \frac{{\alpha + 1}}{2}\lg \left[ {1 + \frac{1}{{\alpha {\sigma ^2}}}{{\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right)}^2}} \right] $ (10)

式中:α为自由度参数;σ为尺度参数。特别地,当α=1时,式(10)为柯西分布。大量试算表明,α=2、σ=1是一组较好的参数组合,适用于绝大多数情况。

采用梯度导引类算法(最速下降或共轭梯度)求解,需要计算目标泛函关于模型参数的梯度,L2模和T分布目标泛函的梯度公式可分别写为

$ \frac{{\partial {J_{{{\rm{L}}_2}}}}}{{\partial \mathit{\boldsymbol{m}}}} = {\mathit{\boldsymbol{L}}^*}\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right) $ (11)
$ \frac{{\partial {J_{\rm{T}}}}}{{\partial \mathit{\boldsymbol{m}}}} = {\mathit{\boldsymbol{L}}^*}\frac{{\alpha + 1}}{{\alpha {\sigma ^2} + {{\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right)}^2}}}\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right) $ (12)

式中上角“*”表示共轭转置。由式(11)、式(12)可以看出,两种范数的不同之处在于:L2模直接利用波场残差计算梯度,T分布则利用加权后的残差记录计算梯度,其他更新流程完全相同。因此,基于常规L2模的各种加速方法都可直接应用于T分布,如相位编码技术。

LSRTM的计算量过于庞大,从而限制了其推广应用。考虑到LSRTM的计算量与炮数成线性关系,因而通过相位编码技术[28]将多个炮集组合成一个超道集,可有效减小计算量,在此基础上应用共享存储并行编程(OpenMP)技术,可进一步提高计算效率,本文选用动态震源极性编码[29]

3 模型试算

本文给出三个模型算例,以验证本文算法的有效性和适用性。

3.1 倾斜界面模型

以倾斜界面的网格化离散为例,说明贴体网格在刻画不规则界面时的优势。图 3为倾斜界面模型剖分示意图。由图可见:常规矩形网格不能很好地离散倾斜界面,在地表界面处产生很多阶梯状的毛刺(图 3a);贴体网格则可很好地逼近复杂界面(图 3b),不仅可以完全拟合真实界面,且在界面处满足很好的正交性,有利于实施边界条件。

图 3 倾斜界面模型剖分示意图 (a)常规矩形网格;(b)贴体网格
速度为3000m/s,地表坡度为6°

图 4为0.3s时刻波场快照。由图可见:常规有限差分算法对倾斜界面的离散作用较差,产生较强的虚假绕射和散射,严重干扰了有效波场(图 4a);贴体全交错网格算法很好地拟合了倾斜界面(图 4b)。

图 4 0.3s时刻波场快照 (a)常规有限差分算法;(b)贴体全交错网格算法
纵、横向网格数目均为201,纵、横向网格间距约为10m。震源为主频20Hz的雷克子波,位于(1000m,400m)处,时间采样间隔为0.6ms
3.2 起伏地表洼陷模型

图 5为起伏地表洼陷模型,图 6为贴体网格剖分及其山峰处的局部放大图。可见,贴体网格对复杂界面的适应性较强。图 7为本文算法得到的RTM结果及其Laplace滤波结果。由图可见:首先,在本文算法得到的RTM结果中地下构造基本被低频噪声掩盖(图 7a)。其次,RTM结果的Laplace滤波结果虽然可正确成像地下构造,但仍然存在如下问题(图 7b):①偏移噪声大。Laplace滤波不能彻底去除低频噪声,且还引入了高频噪声;②反射同相轴中间能量强、两侧能量较弱,即振幅均衡性不佳;③由于地下照明强度随深度的增大而减弱,因而RTM结果深部能量较弱,振幅保真性差。

图 5 起伏地表洼陷模型 模型网格数为301×201,纵、横向网格间距均为10m,最大高程差达400m

图 6 贴体网格剖分(a)及其山峰处的局部放大图(b)

图 7 本文算法得到的RTM结果(a)及其Laplace滤波(b) 301个检波器均匀分布在地表以下10m处,道间距为10m,第1炮位于(0, 20m),炮间距为30m,共101炮。时间间隔为0.6ms,最大记录时间为1.8s,震源为主频20Hz的雷克子波

采用本文提出的起伏地表LSRTM算法可以有效地解决常规RTM存在的问题。然而LSRTM的计算量过于庞大,目前的计算机资源无法满足运算需要。多震源技术可有效缓解计算量问题,但会引入串扰噪声,采用动态相位编码技术可很好地压制串扰噪声,在大幅降低计算量的同时,可得到与常规LSRTM算法相当的结果。将101炮地震数据利用震源极性编码方式组合成一个超道集,使计算量相当于单炮情形,从而大大缓解了计算需求。图 8为基于相位编码的起伏地表LSRTM成像结果。由迭代30次的LSRTM成像结果可见,地下构造清晰,在振幅保真性、均衡性、压制低频噪声等方面较常规RTM结果明显改善,但存在由编码引入的较强高频串扰噪声(图 8a);迭代80次的LSRTM成像结果与理论反射率模型非常接近,有效地压制了串扰噪声(图 8b)。

图 8 基于相位编码的起伏地表LSRTM成像结果 (a)迭代30次;(b)迭代80次

图 9为残差曲线与起伏地表洼陷模型上界面中点处振幅谱。由图可见:①随着迭代次数增加,数据残差和模型残差都逐渐减小,开始下降趋势明显,随着迭代次数的进一步增大,下降趋势减慢;由于目标泛函是数据空间的拟合,因此数据残差能收敛到更低值(图 9a)。②随着迭代次数增大,高频成分逐渐得到恢复,振幅趋近于真值,表明LSRTM可以反演高频成分(图 9b)。

图 9 残差曲线(a)与起伏地表洼陷模型上界面中点处振幅谱(b)

测试发现,三种范数对随机噪声的容忍度基本相同,噪声较弱时效果较好,噪声较强时效果较差。当数据中含有异常值时,混合模及T分布较L2模具有较大优势。图 10为含噪单炮记录,图 11为含脉冲噪声的LSRTM结果。由图可见:L2模LSRTM结果含有较多干扰(图 11a);混合模(图 11b)和T分布(图 11c)LSRTM结果显著消除了脉冲噪声的影响,但前者需要多次尝试以寻求最优参数,因此T分布LSRTM更简单、实用。

图 10 含噪单炮记录 (a)仅含随机噪声;(b)既含随机噪声又包含脉冲噪声

图 11 含脉冲噪声的LSRTM结果 (a)L2模;(b)混合模;(c)T分布
3.3 起伏地表中原模型

图 12为速度模型及其反射率模型。由图可见,起伏地表中原模型不仅地表起伏剧烈,且地下构造复杂,可检验偏移算法的成像效果。采用动态相位编码技术将218炮数据组合成一个超道集,然后利用本文算法迭代计算。图 13为起伏地表中原模型LSRTM结果。由图可见:在1次迭代的LSRTM结果中存在非常强的近地表低频噪声,不能识别地下构造,且由编码引入的高频串扰也很强(图 13a);20次迭代的LSRTM结果(图 13b)已明显压制了低频噪声,基本可识别地下构造,但仍存在较强的高频串扰,且横向振幅均衡性较差;随迭代次数进一步增大,LSRTM结果质量变好,基本消除了高频串扰和低频噪声,显著提高了振幅均衡性和保真性(图 13c图 13d)。

图 12 速度模型(a)及其反射率模型(b) 网格数为436×250,网格间距约为10m。时间采样间隔为0.8ms,最大记录时间为3.2s,爆炸震源为主频20Hz的雷克子波。436个检波器均匀分布于地下10m处,道间距为10m,共218炮均匀分布在地表以下30m处,炮间距为20m

图 13 起伏地表中原模型LSRTM结果 (a)迭代1次;(b)迭代20次;(c)迭代40次;(d)迭代80次
4 结束语

针对起伏地表地震偏移成像存在的问题,本文发展了基于贴体全交错网格的起伏地表LSRTM算法。通过理论分析及模型试算,得到以下认识:

(1) 贴体网格能很好地拟合起伏界面,避免了矩形网格中由于阶梯离散产生的虚假绕射波,且可适应任意复杂地表情形。在曲线坐标系下采用全交错网格,避免了标准交错网格的波场插值和同位网格的高频振荡问题,提高了模拟精度,减小了实现复杂度。

(2) 与常规RTM相比,LSRTM可实现真振幅成像,提高了分辨率、减小了偏移噪声,显著提高了振幅均衡性和保真性。本文测试的三种范数在随机噪声下的表现基本相同,但相比常规L2模,混合模与T分布能较好地压制脉冲噪声的影响,且T分布不依赖于参数选取,因此更简单、实用。

(3) 在起伏地表LSRTM算法中,加入动态相位编码技术不仅能极大地降低计算成本,且可有效压制串扰噪声,在此基础上应用OpenMP技术,可显著提高计算效率,从而将LSRTM的计算成本降低到同常规RTM相同的水平。

尚需指出,本文算法仅进行了理论模型测试,还未用于实际资料处理,且采用的模型近地表变化较为平缓,在横向变速剧烈的近地表条件下算法的适应性还需进一步探讨。此外,如何将其推广到GPU等快速计算设备上进一步提高计算效率,如何通过预条件和正则化算子等增加算法的稳定性并加快收敛速度等是下一步的研究方向。

参考文献
[1]
张华, 李振春, 韩文功. 起伏地表条件下地震波数值模拟方法综述[J]. 勘探地球物理进展, 2007, 30(5): 334-339.
ZHANG Hua, LI Zhenchun, HAN Wen'gong. Review of seismic wave numerical simulation from irregular topography[J]. Progress in Exploration Geophysics, 2007, 30(5): 334-339.
[2]
孙建国. 复杂地表条件下地球物理场数值模拟方法评述[J]. 世界地质, 2007, 26(3): 345-362.
SUN Jianguo. Methods for numerical modeling of geo-physical fields under complex topographical conditions:a critical review[J]. Global Geology, 2007, 26(3): 345-362. DOI:10.3969/j.issn.1004-5589.2007.03.015
[3]
Claerbout J F. Earth Soundings Analysis:Processing versus Inversion[M]. Blackwell Scientific Publications, Boston, 1992.
[4]
Tarantola A. Linearized inversion of seismic reflection data[J]. Geophysical Prospecting, 1984, 32(6): 998-1015. DOI:10.1111/j.1365-2478.1984.tb00751.x
[5]
Nemeth T, Wu C J, Schuster G T. Least-squares mig-ration of incomplete reflection data[J]. Geophysics, 1999, 64(1): 208-221. DOI:10.1190/1.1444517
[6]
陈生昌. 基于反射波动方程的地震反射数据波形成像[J]. 石油地球物理勘探, 2018, 53(3): 487-501.
CHEN Shengchang. Seismic reflection data waveform imaging based on reflection wave equation[J]. Oil Geophysical Prospecting, 2018, 53(3): 487-501.
[7]
刘玉金, 李振春, 吴丹, 等. 局部倾角约束最小二乘偏移方法研究[J]. 地球物理学报, 2013, 56(3): 1003-1011.
LIU Yujin, LI Zhenchun, WU Dan, et al. The research on local slope constrained least-squares migration[J]. Chinese Journal of Geophysics, 2013, 56(3): 1003-1011.
[8]
黄建平, 李振春, 孔雪, 等. 碳酸盐岩裂缝型储层最小二乘偏移成像方法研究[J]. 地球物理学报, 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[J]. Chinese Journal of Geophysics, 2013, 56(5): 1716-1725.
[9]
Kuehl H, Sacchi M D. Least-squares wave-equation migration for AVP/AVA inversion[J]. Geophysics, 2003, 68(1): 262-273. DOI:10.1190/1.1543212
[10]
Clapp M L, Clapp R G, Biondi B L.Regularized least-squares inversion for 3-D subsalt imaging[C]. SEG Technical Program Expanded Abstracts, 2005, 24: 1814-1817. https://www.onepetro.org/conference-paper/SEG-2005-1814
[11]
杨其强, 张叔伦. 最小二乘傅里叶有限差分法偏移[J]. 地球物理学进展, 2008, 23(2): 433-437.
YANG Qiqiang, ZHANG Shulun. Least-squares Fourier finite-difference migration[J]. Progress in Geophysics, 2008, 23(2): 433-437.
[12]
Tang Y. Target-oriented wave-equation least-squares migration/inversion with phase-encoded Hessian[J]. Geophysics, 2009, 74(6): WCA95-WCA107. DOI:10.1190/1.3204768
[13]
周华敏, 陈生昌, 任浩然, 等. 基于照明补偿的单程波最小二乘偏移[J]. 地球物理学报, 2014, 57(8): 2644-2655.
ZHOU Huamin, CHEN Shengchang, REN Haoran, et al. One-way wave equation least-squares migration based on illumination compensation[J]. Chinese Journal of Geophysics, 2014, 57(8): 2644-2655.
[14]
Dai W, Fowler P, Schuster G T. Multi-source least-squares reverse time migration[J]. Geophysical Prospecting, 2012, 60(4): 681-695. DOI:10.1111/j.1365-2478.2012.01092.x
[15]
黄建平, 曹晓莉, 李振春, 等. 最小二乘逆时偏移在近地表高精度成像中的应用[J]. 石油地球物理勘探, 2014, 49(1): 107-112.
HUANG Jianping, CAO Xiaoli, LI Zhenchun, et al. Least square reverse time migration in high resolution imaging of near surface[J]. Oil Geophysical Prospecting, 2014, 49(1): 107-112.
[16]
郭振波, 李振春. 最小平方逆时偏移真振幅成像[J]. 石油地球物理勘探, 2014, 49(1): 113-120.
GUO Zhenbo, LI Zhenchun. True-amplitude imaging based on least-squares reverse time migration[J]. Oil Geophysical Prospecting, 2014, 49(1): 113-120.
[17]
李庆洋, 黄建平, 李振春, 等. 去均值归一化互相关最小二乘逆时偏移及其应用[J]. 地球物理学报, 2016, 59(8): 3006-3015.
LI Qingyang, HUANG Jianping, LI Zhenchun, et al. Mean-residual normalized cross-correlation least-squares reverse time migration and its application[J]. Chinese Journal of Geophysics, 2016, 59(8): 3006-3015.
[18]
李庆洋, 黄建平, 李振春, 等. 优化的多震源最小二乘逆时偏移[J]. 石油地球物理勘探, 2016, 51(2): 334-341.
LI Qingyang, HUANG Jianping, LI Zhenchun, et al. Optimized multi-source least-squares reverse time mig-ration[J]. Oil Geophysical Prospecting, 2016, 51(2): 334-341.
[19]
薛浩, 刘洋, 杨宗青. 基于优化时空域频散关系的声波方程有限差分最小二乘逆时偏移[J]. 石油地球物理勘探, 2018, 53(4): 745-753.
XUE Hao, LIU Yang, YANG Zongqing. Least-square reverse time migration of finite-difference acoustic wave equation based on an optimal time-space dispersion relation[J]. Oil Geophysical Prospecting, 2018, 53(4): 745-753.
[20]
Zhang W, Chen X F. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[J]. Geophysical Journal of the Royal Astronomical, 2010. DOI:10.1111/j.1365-246X.2006.03113.x
[21]
Lisitsa V, Vishnevskiy D. Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity double dagger[J]. Geophysical Prospecting, 2010, 58(4): 619-635. DOI:10.1111/j.1365-2478.2009.00862.x
[22]
Brossier R, Operto S, Virieux J. Which data residualnorm for robust elastic frequency-domain full waveform inversion[J]. Geophysics, 2010, 75(3): R37-R46. DOI:10.1190/1.3379323
[23]
Woodon J, Minji K, Shinwoong K, et al. Full waveform inversion using student's t distribution:anumerical study for elastic waveform inversion and simultaneous-source method[J]. Pure and Applied Geophysics, 2015, 172(6): 1491-1509. DOI:10.1007/s00024-014-1020-7
[24]
Aravkin A, Leeuwen T V, Herrmann F.Robust full-waveform inversion using the student's T-distribution[C]. SEG Technical Program Expanded Abstracts, 2011, 30: 2669-2673.
[25]
李庆洋, 黄建平, 李振春, 等. 起伏地表贴体全交错网格仿真型有限差分正演模拟[J]. 石油地球物理勘探, 2015, 50(4): 633-642.
LI Qingyang, HUANG Jianping, LI Zhenchun, et al. Undulating surface body-fitted grid seismic modeling based on fully staggered-grid mimetic finite difference[J]. Oil Geophysical Prospecting, 2015, 50(4): 633-642.
[26]
田坤, 黄建平, 李振春, 等. 多轴卷积完全匹配层吸收边界条件[J]. 石油地球物理勘探, 2014, 49(1): 143-152.
TIAN Kun, HUANG Jianping, LI Zhenchun, et al. Multi-axial convolution perfectly matched layer absorption boundary condition[J]. Oil Geophysical Prospecting, 2014, 49(1): 143-152.
[27]
李娜. 基于Huber范数的多震源最小二乘逆时偏移[J]. 石油地球物理勘探, 2017, 52(5): 941-947.
LI Na. Multi-source least-squares reverse time migration based on Huber norm[J]. Oil Geophysical Prospecting, 2017, 52(5): 941-947.
[28]
Romero L A, Ghiglia D C, Ober C C, et al. Phase encoding of shot records in prestack migration[J]. Geophysics, 2000, 65(2): 426-436. DOI:10.1190/1.1444737
[29]
Schuster G T, Wang X, Huang Y, et al. Theory of multisource crosstalk reduction by phase-encoded statics[J]. Geophysical Journal International, 2011, 184(3): 1289-1303. DOI:10.1111/j.1365-246X.2010.04906.x