石油物探  2019, Vol. 58 Issue (5): 625-644  DOI: 10.3969/j.issn.1000-1441.2019.05.001
0
文章快速检索     高级检索

引用本文 

曲英铭. 起伏地表直接成像技术研究进展[J]. 石油物探, 2019, 58(5): 625-644. DOI: 10.3969/j.issn.1000-1441.2019.05.001.
QU Yingming. Research progress of topographic imaging methods[J]. Geophysical Prospecting for Petroleum, 2019, 58(5): 625-644. DOI: 10.3969/j.issn.1000-1441.2019.05.001.

基金项目

国家自然科学基金(41774133, 41904101)、国家科技重大专项(2016ZX05024-003-011, 2016ZX05006-002-003, 2016ZX05002-005-007HZ, 2016ZX05014-001-008HZ, 2016ZX05026-002-002)、国家重点研发计划(2016YFC060110501)、中央高校基本科研业务费专项资金(19CX2010A)、中国石化地球物理重点实验室开放基金(wtyjy-wx2017-01-04, wtyjy-wx2018-01-06)、山东自然科学基金(ZR2019QD004)、中国石油大学(华东)人才引进费(20180041)共同资助

作者简介

曲英铭(1990—), 男, 博士, 副教授, 硕士生导师, 主要从事地震波传播理论、成像与反演方法等方面的研究工作。Email:quyingming@upc.edu.cn

文章历史

收稿日期:2019-06-26
改回日期:2019-07-19
起伏地表直接成像技术研究进展
曲英铭1,2     
1. 中国石油大学(华东)地球科学与技术学院地球物理系, 山东 青岛 266580;
2. 中国石油化工股份有限公司地球物理重点实验室, 江苏 南京 211103
摘要:常用的高程静校正方法需要满足地表一致性假设的基本条件, 但在我国西部与南方山前带油气勘探区内, 地表起伏剧烈, 横向变速明显, 地表一致性假设不再成立。为了对复杂山前带进行准确的成像, 起伏地表直接成像方法得到了广泛研究与快速发展。对起伏地表直接成像技术的研究进展进行了综述, 包括基于射线理论的起伏地表偏移, 基于单程波方程的起伏地表偏移及基于双程波方程的起伏地表变网格、曲网格、三角网格逆时偏移, 起伏地表粘声逆时偏移, 起伏地表各向异性逆时偏移, 起伏地表最小二乘偏移等, 并对起伏地表速度反演方法进行了总结。从简单的起伏地表构造成像到同时存在剧烈起伏地表及复杂地下结构的双复杂构造直接成像回顾了起伏地表成像的总体发展趋势, 指出基于双程波方程的起伏地表直接成像方法是未来的主要研究方向, 该方法在理论上可以对任意复杂构造进行精确成像, 对速度的依赖性强和计算量庞大是其主要的发展瓶颈。
关键词起伏地表    直接成像    叠前深度偏移    射线类    单程波    逆时偏移(RTM)    最小二乘逆时偏移    速度建模    各向异性    
Research progress of topographic imaging methods
QU Yingming1,2     
1. Department of Geophysics, School of Geosciences, China University of Petroleum (East China), Qingdao 266580, China;
2. SINOPEC Key Laboratory of Geophysics, Nanjing 211103, China
Foundation item: This research is financially supported by the National Natural Science Foundation of China (Grant Nos.41774133, 41904101), the National Science and Technology Major Project of China (Grant Nos.2016ZX05024-003-011, 2016ZX05006-002-003, 2016ZX05002-005-007HZ, 2016ZX05014-001-008HZ, 2016ZX05026-002-002), the National Key Research and Development Program(Grant No.2016YFC060110501), the Fundamental Research Funds for the Central Universities (Grant No.19CX2010A), the Open Funds of SINOPEC Key Laboratory of Geophysics (Grant Nos.wtyjy-wx2017-01-04, wtyjy-wx2018-01-06), the Natural Science Foundation of Shandong Province (Grant No.ZR2019QD004), and the Talent Introduction Fund of China University of Petroleum (East China) (Grant No.20180041)
Abstract: The assumption of surface regularity is the basis for the commonly used elevation-static correction method.In some areas such as piedmont area of western and southern China, the assumption is hardly satisfied if severe surface undulations are present, or evident lateral variations of near-surface velocity exist.Several methods for topographic images were developed in recent years, with the aim of imaging piedmont areas accurately.In this study, various topographic pre-stack depth migration methods have been reviewed.These methods, by which the imaging is performed on the rugged surface directly, include the topographic pre-stack depth migration based on ray theory, the one-way wave equation, and the two-way wave equation (reverse time migration, RTM).The RTM can be implemented on the basis of dual-variable grids, curvilinear grids, or triangular grids.In addition to these methods, the visco-acoustic RTM, anisotropic RTM, and least-squares RTM (LSRTM) for rugged topography have been reviewed, and the methods for topographic velocity inversion also have been summarized.The general trend in the technologies for imaging undulated surfaces shows a transition from the direct imaging of undulated surfaces to the direct imaging of dual-complexity structures, featuring both undulated surfaces and complex subsurface structures.Methods for topographic imaging based on the two-way wave equation represent the main focus of future research, for which the main challenges are the heavy dependence on high-precision velocity measurements and the large computational costs.
Keywords: rugged topography    topographic imaging    prestack depth migration    ray theory    one wave equation    reverse time migration    least-squares reverse time migration (LSRTM)    velocity modeling    anisotropy    

随着我国油气资源的日益短缺以及东部油气勘探程度的不断提高, 地震勘探的重点逐步转向了勘探难度更大的西部与南方。这些地区起伏剧烈的地表及复杂的近地表结构给地震勘探工作带来了极大的困难与挑战[1]:由于地表剧烈起伏、近地表横向变速剧烈, 地表一致性假设不再满足, 高程静校正后的成像结果会发生畸变, 成像精度无法满足地质解释的需求。

起伏地表的地震资料处理经历了以下发展过程:常规静校正、改进的高精度静校正、基准面校正[2-4]以及起伏地表直接成像方法。其中基于起伏地表构造直接进行成像是解决复杂近地表地区成像问题的根本方法。

起伏地表直接成像方法总体上包括以下3种:射线类方法[5-6]、单程波方法[7-9]及双程波方法[10-11]。射线类成像方法对复杂近地表条件具有较好的适应性, 可直接在起伏地表上进行波场延拓与偏移成像[12-14]。该方法主要分Kirchhoff偏移法和束偏移法两种。Kirchhoff偏移法由绕射扫描叠加方法衍生而来, 基于地震记录的加权绕射叠加, 应用波动方程的积分解进行地震波场的反向传播与成像[12]。为了压制Kirchhoff偏移方法产生的偏移噪声, GRAY[13]将束偏移引入到Kirchhoff方法中, 既保留了Kirchhoff方法计算速度快、对复杂地表适应性强的优点, 又提高了偏移剖面的信噪比。岳玉波等[15]基于局部平面波分解, 提出了一种起伏地表保幅高斯束成像方法, 不仅能够克服起伏地表的影响, 还能得到较为保幅的地下构造成像结果。杨继东等[16]将改进的高斯束-菲涅尔束应用到起伏地表束偏移方法中, 提高了起伏地表平面的分解精度, 均衡了浅、中、深层的成像精度。黄建平等[17]将岳玉波等[15]提出的起伏地表保幅高斯束成像方法发展到起伏地表弹性波。起伏地表单程波成像方法主要包括零速层法[18]、逐步累积法[19]、直接下延法[20]和波场上延法[21]等。直接下延法和波场上延法是由零速层法与逐步累积法相结合发展而来的。学者们对波场延拓算子做了许多改进:程玖兵等[22]提出了基于优化系数傍轴近似方程的频率-空间域偏移算子, 该算子能较好地适用于横向变速情况下的复杂近地表构造成像; 王成祥等[23]提出了一种混合法起伏地表单程波偏移方法, 提高了起伏地表条件下的成像精度; 叶月明等[24-25]在频率-空间域实现了带误差补偿的起伏地表有限差分单程波深度偏移方法, 压制了常规频率-空间域有限差分单程波偏移方法产生的偏移噪声, 并引入保幅算子, 提出了双复杂构造单程波保幅偏移方法。基于双程波的逆时偏移方法[26]对于高陡构造等复杂地质体成像精度高, 近年来随着计算机技术的发展已成为学者们研究的热点。关于逆时偏移的研究主要针对提高计算效率[27-28]、优化成像条件[29]、压制成像噪声[30-33]以及起伏地表条件下的改善成像质量[34-36]等。

本文对起伏地表直接成像技术的研究进展进行综述, 主要针对射线类、单程波及双程波起伏地表叠前深度偏移方法。

1 基于射线理论的起伏地表偏移成像

1984年, WIGGINS[12]提出了一种基于起伏地表地震数据的Kirchhoff积分公式及偏移方法, GRAY[13]验证了Kirchhoff起伏地表直接成像方法的成像优势, JAGER等[14]优化并发展了起伏地表Kirchhoff成像方法, 提出了保幅Kirchhoff起伏地表直接成像方法。Kirchhoff成像方法计算效率高, 能够较好地处理起伏地表, 因此也是目前工业界常用的偏移方法。但该方法对复杂构造成像精度较低, 且存在严重的偏移噪声。

束偏移方法实际上是一种改进的Kirchhoff成像方法。该方法通过倾斜叠加将局部区域地震数据分解为不同方向的波数并进行成像[37-41]。2005年, GRAY[13]提出了一种基于局部静校正的起伏地表高斯束成像方法(图 1), 但当地表起伏剧烈时, 静校正会导致地震波场产生畸变, 影响地下构造成像的精度。针对这个问题, 岳玉波等[15]提出了针对起伏地表的保幅高斯束成像方法, 将高斯束格林函数(图 2a)与起伏地表条件下利用瑞雷积分计算的逆时传播的地震波场(图 2b)结合起来, 推导并实现了起伏地表条件下的高斯束逆时延拓算子, 利用反褶积成像条件实现了起伏地表保幅高斯束成像。该方法的优点是在复杂的近地表处直接进行局部平面波分解, 消除了几何扩散对地震波振幅的影响, 从而消除了剧烈起伏地表的影响, 且对成像结果进行了保幅。图 3对比了采用岳玉波等[15]提出的起伏地表保幅高斯束成像方法和GRAY[13]提出的基于局部静校正的起伏地表高斯束成像方法对Canadian Foothills模型成像的结果。Canadian Foothills模型存在起伏剧烈的地表和复杂的近地表速度变化, 地震成像极为困难。由图 3可见, 相对基于局部静校正的起伏地表高斯束成像结果(图 3c), 基于保幅延拓法的起伏地表高斯束成像结果(图 3b)成像精度和信噪比更高。

图 1 局部静校正法[13]
图 2 高斯束格林函数(a)和复杂地表条件下的反向延拓波场(b)[15]
图 3 Canadian Foothills模型及叠前深度偏移结果对比 a Canadian Foothills模型; b基于保幅延拓法的起伏地表高斯束成像; c基于局部静校正的起伏地表高斯束成像[13]

采用起伏地表保幅高斯束成像方法对A探区起伏地表实际资料(图 4)进行了偏移处理。图 4a是该探区的速度场, 图 4b是某一炮的单炮记录, 图 4c是起伏地表的高程变化, 该地区最大高程为58m, 局部变化剧烈。实际资料共计329炮, 最大道数为320道, 采样点数为1501, 采样率为4ms, 道间距为50m, 炮间距为100m, CDP间距为25m;速度场横向1297个CDP, 纵向深度7560m, 深度采样间隔是4m。图 5为A探区起伏地表实际资料保幅高斯束偏移成像结果, 可以看出起伏地表的影响得到了很好的消除。

图 4 A探区起伏地表实际资料 a速度场; b单炮记录; c起伏地表高程
图 5 A探区起伏地表实际资料保幅高斯束偏移成像结果

与波动方程类偏移方法相比, 基于射线理论的偏移方法本身就能够很好地适应起伏地表构造, 因此, 各种新的射线类偏移成像方法都可以发展为起伏地表直接成像方法。随着弹性波多分量起伏地表高斯束成像方法的发展, 射线类起伏地表直接成像方法已经较为成熟。此外, 射线类成像方法具有计算速度快的优势, 因此, 射线类起伏地表直接成像方法在工业界得到了广泛应用。但是, 该方法基于射线理论, 主要反映地下构造的运动学特征, 难以准确地对地下复杂构造的动力学特征进行刻画, 且该方法无法对高陡构造等复杂构造进行成像。

2 基于单程波方程的起伏地表保幅叠前深度偏移

随着波动方程叠前深度偏移技术的逐渐完善以及在复杂构造地区的广泛应用, 波动方程偏移已经被证实是复杂构造成像较为理想的途径。相对Kirchhoff偏移方法, 波动方程叠前深度偏移方法在构造成像方面有着明显的优势。基于波动方程的保幅偏移成像分两步:第一步是波场传播, 第二步是建立共成像点道集的成像原则。SAVA等[42]基于波动方程偏移建立了角度域共成像点道集(CIGs), 并建立了一种保幅成像原则。朱绪峰[43]对保幅单程波传播理论进行了研究。

基准面校正是地震数据处理必不可少的一步, 在近地表横向速度变化剧烈以及地表起伏地区, 基准面校正更加重要。鉴于常规高程基准面校正无法适应起伏地表和近地表速度横向变化剧烈地区, BERRYHILL[44]提出了针对叠后数据的波动方程基准面校正方法, 并将该方法发展到叠前数据[45]。波动方程波场外推技术可以将野外地震数据从地面延拓到任一平面, 从而消除了地形剧烈变化对构造成像的影响。BEASLEY等[18]提出“零速层”方法, 将水平基准面挪到了地表最高点或最高点以上, 在该水平基准面和实际地表之间插入虚拟层, 虚拟层的速度用常速度(零速度或非常小的速度)填充, 然后从水平基准面开始进行偏移, 这期间只考虑地震波的垂向传播, 遇到实际地层再恢复正常。该方法的优势在于无需做高程静校正, 只对速度模型做微小改变, 巧妙地化解了起伏地表的影响。1991年, RESHEF[19]提出了“逐步-累加”的波场延拓思想。近年来, 学者们基于BEASLEY等[18]和RESHEF[19]的方法提出了“波场上延”法和“直接下延”法, 有效地压制了地表剧烈变化的影响。如2006年, 田文辉等[2]提出了复杂地表和地下地质条件下的直接下延波动方程叠前深度偏移方法; 叶月明等[24-25]应用直接下延法研究双复杂介质条件下的叠前深度偏移, 在保持较高计算效率的同时, 有效地改善了双复杂介质的成像效果。

常规傅里叶有限差分法(FFD)在速度变化很大的情况下不稳定。为了克服这一缺点, MILLINAZZO等[3]对基于Pade逼近得到的平方根算子做复Pade逼近, 推导出基于复Pade逼近的傅里叶有限差分算子, 通过旋转平方根的分支截断得到较好的成像效果。该方法的优点在于减少了偏移噪声, 同时对延拓步长的要求不高, 提高了计算效率。图 6对比了该方法和常规FFD方法用于2D-SEG起伏地表模型(图 6a)的成像结果[20], 相对常规FFD(图 6b), 基于复Pade逼近的FFD方法(图 6c)不仅压制了偏移噪声, 而且很好地改善了深层成像效果。

图 6 2D-SEG起伏地表模型叠前深度偏移 a 2D-SEG速度模型; b常规FFD偏移方法; c基于复Pade逼近的FFD偏移方法[20]

采用基于复Pade逼近的FFD方法对A探区起伏地表实际资料(图 4)进行了偏移处理测试。图 7a图 7c分别是基于常规FFD算子、旋转角度α=5°的复Pade逼近FFD算子和复Pade逼近的保幅FFD算子, 采用直接下延法进行叠前深度偏移的成像结果。可以看出:基于复Pade逼近的FFD算子的成像效果要好于常规FFD算子的成像效果, 偏移噪声得到了明显压制; 基于复Pade逼近的保幅FFD算子的成像效果要好于常规FFD算子的成像效果和基于复Pade逼近的FFD算子的成像效果, 在压制偏移噪声的同时, 有效保持了深层反射能量, 浅中层和深层的构造形态都得到了较好的成像。

图 7 A探区起伏地表实际资料叠前深度偏移 a常规FFD偏移方法; b基于复Pade逼近的FFD方法; c基于复Pade逼近的保幅FFD方法

波动方程成像方法相对射线理论成像方法能够更好地对复杂地下构造进行成像。在计算机技术无法承受逆时偏移巨大的计算量和内存压力之前, 单程波成像方法得到了广泛应用。但基于单程波方程的成像方法对起伏地表的适应性差, 需要对起伏地表进行特殊处理。当学者们认识到起伏地表直接成像的必要性之后, 基于单程波方程的起伏地表直接成像方法发展迅速, 出现了很多用于消除剧烈起伏地表影响的单程波成像方法。随着计算机技术的快速发展及逆时偏移技术的推广应用, 近些年基于单程波方程的起伏地表直接成像方法发展较为缓慢, 这是因为单程波成像方法在成像精度上不如逆时偏移方法, 在计算效率和对起伏地表的适应性上不如射线类方法。

3 基于双程波方程的起伏地表叠前逆时偏移 3.1 双程波波场延拓算子

逆时偏移跟正演密不可分, 正演技术的发展一定程度上带动了逆时偏移的发展。关于起伏地表对地震波场的影响研究[46]可以追溯到20世纪40年代, 只是因为当时技术不发达, 这项工作一直没有得到展开。到了20世纪70年代, 随着计算机技术的迅速发展, 地震模拟技术也随之迅速发展, 20世纪八、九十年代掀起了起伏地表正演模拟研究的热潮。正演模拟总体上分为两类:射线类和波动方程类。射线类方法最早由WIGGINS[12]提出, 证明了起伏地表情况下Kirchhoff积分偏移的适应性, 但该方法终因基于高频近似、在速度横向变化剧烈地区适应性差等缺陷发展缓慢。波动方程类方法是对建立的地质模型进行网格划分, 划分后的模型由有限个离散点组成。该方法对介质没有横向变化的限制, 如果网格足够小, 得到的解会非常精确, 而且综合考虑了地震波的运动学和动力学特征, 因此得到了广泛应用。

几十年来, 解决起伏地表问题的波动方程方法基本可以归为四大类:有限差分法、有限元法、谱元法和边界元法。其中有限差分法因其良好的模拟精度和计算效率应用最为广泛, 该方法用差分算子代替微分算子, 将波动方程和计算区域离散化, 通过不断更新迭代获得各个离散点上的波场值。20世纪60年代, ALTERMAN等[47]开创性提出了层状介质中的弹性波有限差分方法, 自此有限差分方法在地震勘探领域不断发展。1986年, VIRIEUX[48]提出了一阶速度-应力波动方程交错网格有限差分法, 随后又实现了SH波、P-SV波的波场模拟, 极大地提高了波场模拟精度和计算效率, 且没有增加存储空间。针对起伏地表频散严重的问题, JASTRAM等[49]采用可变网格的方法进行网格剖分, TESSMER等[50]将坐标变换的思想应用于起伏地表。董良国[51]用纵向坐标变换实现了起伏地表条件下的弹性波数值模拟, 之后又基于泰勒展开发展了非规则网格差分方法。褚春雷[52]在前人研究的基础上采用non-Sibson插值方法进行声波方程正演模拟以及逆时偏移, 取得了不错的效果。有限元法是变分和剖分插值, 该方法的核心部分是有限元网格的生成, 目前有各向异性网格、贴体网格、并行网格和自适应网格等。1974年, THOMPSON等[53]系统研究了贴体网格, 提出了偏微分数值求解方法, 证明了Laplace方程和Poisson方程可以实现保角变换。孙建国等[54]采用自适应网格技术生成网格, 并将其应用于复杂地表情况下的地震波数值模拟中。总体来说, 有限元法的优点是能够很好地处理起伏地表问题, 缺点是内存要求高, 计算量大。谱元法是结合有限元和伪谱法发展起来的一种数值模拟方法。KOMATITSCH等[55]用谱元法实现了起伏地表条件下的二维和三维弹性波数值模拟, 并且分析了起伏地表引起的各种干扰波, 但是未研究地表剧烈变化条件下的有效反射问题。CHE等[56]通过修改Chebyshev谱元算法实现了起伏地表弹性波数值模拟, 证明了该方法计算效率高、收敛速度快, 并对产生的干扰波进行了讨论, 但未研究如何消除干扰波。总体来说, 谱元法具有精度高、稳定性好、耗时少等优点, 但因其计算量仍然较大未得到广泛应用。边界元法的计算量较大, 而且其本身的半解析性质难以适应地表剧烈变化的情况, 因此其应用受到了限制。

3.2 起伏地表逆时偏移 3.2.1 逆时偏移的基本原理

逆时偏移的研究历史可以追溯到20世纪80年代。BAYSAL等[57]和WHITMORE[26]率先提出了RTM的思想, 证明了RTM的优越性, DAN[58]和MUFTI[59]将这一思想应用到地震勘探领域, 然而早期的RTM研究都是针对叠后资料的。后期针对叠前数据的RTM研究可以分为声波方程RTM和弹性波方程RTM两类。LESAGE等[60]和JONES等[61]将声波RTM应用于复杂构造的成像, 证明了声波RTM对于复杂构造成像的优势。随着多分量地震采集技术和计算机技术的发展, 弹性波RTM也得到了广泛的研究[62-65]。弹性波RTM主要包括两种类型:一是标量RTM, 该方法对纵横波地震记录分别进行RTM成像[66], 由于准确分离纵横波记录还存在一定的难度, 因此该方法的应用受到限制。二是矢量弹性波RTM, 将纵横波波场的分离放在应用成像条件之后[67]。在各向同性介质中, 对波场应用亥姆霍兹分解定理可以得到纵横波分离的波场[68]。另外, 在纵横波成像(PS成像和SP成像)中, 需要校正横波的极性[69-70]。虽然RTM对于复杂构造成像具有很大的优势, 但关于如何提高RTM的成像精度和计算效率问题, 还需要进一步研究。近年来, 学者们通过改造波动方程或进行波场分解[71-74]、修改成像条件[75-76]和成像后滤波[34]等提高了RTM的成像精度和计算效率。

对于二维各向同性介质, 声波方程可以通过降阶表示为:

$ \left\{\begin{array}{l}{\frac{\partial u}{\partial t}=-\rho v_{\mathrm{P}}^{2}\left(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{z}}{\partial z}\right)} \\ {\frac{\partial v_{x}}{\partial t}=\frac{1}{\rho} \frac{\partial u}{\partial x}} \\ {\frac{\partial v_{z}}{\partial t}=\frac{1}{\rho} \frac{\partial u}{\partial z}}\end{array}\right. $ (1)

式中:u表示应力; vxvz分别表示x方向和z方向质点的速度; ρ表示密度; vP表示纵波速度。其时间二阶、空间十阶的交错网格有限差分格式可以表示为:

$ \begin{array}{l} U_{i,j}^{k + \frac{1}{2}} = U_{i,j}^{k - \frac{1}{2}} - \Delta t\rho v_{\rm{P}}^2\left[ {\frac{1}{{\Delta x}}\sum\limits_{n = 1}^5 {C_n^{(5)}} \left( {P_{i + \frac{{2n - 1}}{2},j}^k - } \right.} \right.\\ \;\;\;P_{i - \frac{{2n - 1}}{2},j}^k\left. { + \frac{1}{{\Delta z}}\sum\limits_{n = 1}^5 {C_n^{(5)}} \left( {Q_{i,j + \frac{{2n - 1}}{2}}^k - Q_{i,j - \frac{{2n - 1}}{2}}^k} \right)} \right]\\ P_{i + \frac{1}{2},j}^k = P_{i + \frac{1}{2},j}^{k - 1} + \frac{{\Delta t}}{{\Delta x \cdot \rho }}\sum\limits_{n = 1}^5 {C_n^{(5)}} \left[ {U_{i + n,j}^{k - \frac{1}{2}} - U_{i - (n - 1),j}^{k - \frac{1}{2}}} \right]\\ \;\;\;\;\;\;\;\;\;\;\;Q_{i,j + \frac{1}{2}}^{k - 1} = Q_{i,j + \frac{1}{2}}^{k - 1} + \frac{{\Delta t}}{{\Delta z \cdot \rho }}\sum\limits_{n = 1}^5 {C_n^{(5)}} \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {U_{i,j + n}^{k - \frac{1}{2}} - U_{i,j - (n - 1)}^{k - \frac{1}{2}}} \right] \end{array} $ (2)

式中:Δx和Δz分别代表水平方向x、垂直方向z的网格间距; k表示时间上的离散值; UPQ分别代表uvxvz; Δt表示时间步长; ij分别代表水平方向x和垂直方向z的离散值; C为有限差分系数。由此可以进行声波方程交错网格正演, 同时对检波点处波场进行反传, 根据成像条件完成成像。叠前RTM的成像条件包括三类:一是激发时间成像条件; 二是互相关成像条件; 三是震源波场与检波点波场振幅比成像条件。其中互相关成像条件可表示为:

$ I(x, y, z)=\int S(x, y, z, t) R(x, y, z, t) \mathrm{d} t $ (3)

式中:xy是水平坐标; z为深度坐标; t为时间。弹性波动方程逆时偏移的基本原理与声波相似, 这里不再赘述。

3.2.2 变网格逆时偏移

在起伏地表情况下, 采用常规有限差分地震波模拟方法必须取很小的网格间距才能保证计算精度, 而这样会带来巨大的计算量以及局部过采样问题。为了解决这一问题, JASTRAM等[77]基于二维声波方程提出了变网格步长的算法。其后, 马光克等[78]实现了可变网格高阶有限差分算法, 并将其应用于逆时偏移成像。当空间网格变化时, 时间采样需要满足精细网格对应的稳定性条件, 这会增加偏移计算的时间。为此, TAE-SEOB[79]、HUANG等[80]、李振春等[81]先后实现了基于交错网格的时空双变算法(图 8), 并进行了虚假反射误差分析。

图 8 双变网格原理[81]

在进行变网格RTM时, 需要注意的是波场外推过程中局部精细网格和粗糙网格之间的转移。图 9为变网格RTM流程图[81]。利用近地表起伏模型进行了声波RTM测试, 分别采用粗网格和双变网格进行了逆时偏移, 得到的成像结果如图 10所示[81]。可见双变网格RTM的成像效果(图 10c)比粗网格RTM的成像效果(图 10b)好。

图 9 变网格逆时偏移流程[81]
图 10 逆时偏移成像结果 a近地表起伏模型; b粗网格; c双变网格[81]
3.2.3 坐标变换的有限差分逆时偏移

相对于变网格方法, 曲网格对于起伏地表更具优势。曲网格地震正演模拟最初只适用于内部界面[82], 随后有学者根据HESTHOLM等[83]推导的自由边界条件公式, 推导出适用于起伏地表的曲网格显式精确条件, 提高了计算精度。如鲁雁翔等[84]推导了弹性波正向传播和逆时传播曲线网格差分格式, 实现了曲线交错网格有限差分法逆时偏移成像。

近年来, 贴体网格因其对起伏地表良好的适应性而受到关注, 学者们将贴体网格应用于起伏地表正演模拟, 取得了不错的效果。贴体网格生成的原则是使离散后的网格边界与地表形态吻合, 以避免人为产生的阶梯边界引起的虚假散射[53]。ZHANG等[85]采用DRP/opt MacCormack格式计算空间导数, 实现了贴体同位网格上的一阶速度-应力方程有限差分正演模拟, 但该方法计算量较大。随后, APPELO等[86]、兰海强等[87]推导并实现了贴体网格上的二阶波动方程正演模拟算法。然而二阶位移方程在泊松比较大的介质条件下容易出现不稳定, 且该算法不易推广到高阶。CHENG等[88]采用贴体旋转交错网格实现了起伏地表条件下的地震波正演模拟, 李庆洋等[89]发展了曲线坐标系下基于仿真型有限差分的贴体全交错网格正演模拟算法。如图 11所示, 贴体网格可以通过计算空间到物理空间的坐标变换来获得[90]

图 11 计算空间到物理空间的坐标变换[90]

利用含有高斯山峰、山谷的洼陷模型进行了逆时偏移成像测试(图 12)。图 12a为含有高斯山峰、山谷的洼陷模型, 该模型共4层, 速度分别为2500、3500、4000、4500m/s, 模型大小为301×201个网格点, 网格间距为10m, 最大高程差可达400m。图 12b为贴体网格剖分图, 可以看出, 贴体网格对复杂界面具有强适应性及正交性。观测系统:301个检波器均匀分布在地表以下10m处, 道间距10m, 第一炮位于(0, 20m), 炮间距30m, 共101炮。计算参数:时间间隔0.6ms, 最大记录时间1.8s, 震源为主频20Hz的雷克子波。图 12c为声波逆时偏移成像结果, 可以看出, 起伏地表的影响得到了很好的消除, 洼陷构造得到了较准确的成像。

图 12 洼陷模型逆时偏移成像 a含高斯山峰、山谷的洼陷模型; b贴体网格剖分; c贴体网格逆时偏移成像结果[91]
3.2.4 三角网格

非结构性网格适用于复杂区域, 可以控制网格的疏密, 较好地模拟曲界面。目前, 非结构性网格有很多种形式, 其中不规则三角网格作为一种很典型的非结构性网格形式(图 13)[36]受到了国内外学者的广泛关注。如GUO等[92]使用三角网格法处理起伏边界, 其本质就是对任意形状的自由表面边界进行坐标变换, 通过引入不规则网格差分算子, 使应力-速度方程能够适用于任意形状的边界, 将波动方程的差分计算转换到新坐标系下进行。孙小东等[36]利用三角网格进行差分离散, 并将其应用于复杂地表叠前逆时偏移。

图 13 起伏地表三角网格剖分 a凹形地表剖分; b凸形地表剖分[36]

对起伏地表模型(图 14a)进行三角网格剖分, 然后做RTM, 图 14b图 14c分别是模拟记录和逆时偏移的结果。可以看出, 模型深部同相轴得到了很好的成像, 成像剖面的信噪比也比较高。

图 14 起伏地表模型三角网格逆时偏移成像 a起伏地表模型; b模拟记录; c逆时偏移结果[36]
3.2.5 起伏地表粘声逆时偏移

实际地下介质存在明显的粘弹性, 对地震波有强烈的吸收衰减作用。消除粘弹性对地震波传播的影响一直是学者们研究的重要课题。到目前为止, 主要发展了两大类方法:反Q滤波和反Q偏移。其中反Q滤波技术在对地震波的吸收衰减进行补偿的同时改善了相位特征, 早期应用较为广泛[93-94]。随着偏移技术的发展, 反Q偏移技术逐渐成为主流。根据所依据的理论不同, 反Q偏移可分为基于射线理论、基于单程波理论和基于双程波理论3种。基于射线理论的反Q偏移主要有带衰减补偿的Kirchhoff偏移方法[95]、变换域衰减补偿方法[96]以及叠前深度偏移衰减补偿方法[97]。1994年, DAI等[98]率先提出了基于单程波方程的反Q偏移方法, 在进行偏移的同时补偿吸收衰减的影响。随后, 学者们使用单程波方法进行了深度偏移补偿。陈雪等[99]依据单程波方程对粘弹介质进行了地震波场数值模拟, 并对正演及照明结果进行了对比分析。随着逆时偏移的发展, 学者们又提出了基于双程波方程的逆时偏移Q补偿算法, 提高了成像精度。如LIU等[100]推导了一种适用于起伏地表的粘弹介质Q补偿逆时偏移方法。简单说来, 粘声RTM的基本原理就是在常规声波RTM正、反向延拓过程中同时进行吸收衰减补偿。对于起伏地表粘声RTM, 需要同时校正起伏地表和粘弹性对成像的影响。QU等[91]提出了一种粘声起伏地表逆时偏移方法, 推导了适应起伏地表贴体网格的辅助坐标系粘声拟微分方程, 该方法既能够校正地下粘弹性引起的振幅衰减及相位畸变, 又能够克服剧烈起伏地表对成像的影响。通过实际资料试处理对起伏地表粘声逆时偏移方法进行了测试。该资料共1180炮, 每炮906个检波点。偏移速度模型和Q模型如图 15所示, 速度范围为1910~6480m/s, Q的范围为36~110。震源函数为25Hz的雷克子波, 时间采样间隔为2ms, 记录时间为8s。图 16a图 16b分别为粘声和声波起伏地表RTM成像结果, 可见声波RTM的深部构造成像能量弱, 且信噪比较低, 而起伏地表粘声RTM成像结果经过Q补偿, 断层、不整合面和其它构造都得到了很好的成像。图 17为从图 16中水平位置4000m处抽取的波数谱曲线[91], 证明了粘声RTM经过衰减补偿, 成像结果的能量和分辨率都得到了提高。

图 15 实际资料的速度场(a)和Q模型(b)[91]
图 16 起伏地表逆时偏移成像结果对比 a粘声; b声波[91]
图 17图 16水平位置4000m处抽取的波数谱曲线[91]
3.2.6 起伏地表各向异性逆时偏移

不同于粘声逆时偏移在波场补偿过程中需要引入特殊处理压制指数级增长的高频噪声, 各向异性逆时偏移的核心是各向异性介质的正演模拟, 因此这里只介绍各向异性介质数值模拟技术。各向异性介质也是当前学者们研究的热点, 其中最受关注的主要有水平横向各向同性介质(HTI)、垂直横向各向同性介质(VTI)和倾斜横向各向同性介质(TTI)三类, 统称为TI介质。在各向同性介质中, RTM、LSRTM以及全波形反演(FWI)使用的都是标量纵波方程, 但是在TI介质中, 由于存在纵横波耦合问题, 简单的标量纵波方程难以表述地震波的运动学和动力学特征。针对这一问题, 国内外学者提出了很多解决方法:TSVANKIN[101]推导了VTI介质qP-qSV波精确相速度公式, ALKHALIFAH[102]对该公式进行了简化, 提出了著名的“声学近似”公式。简单说来, 各向异性逆时偏移的发展历程, 就是qP波方程的简化历程。根据前人的研究成果, 其大致可分为两类:第一类是qP-qSV波耦合波动方程, 该方程在倾角急剧变化时极不稳定。为了解决这一难题, FLETCHER等[103]引入非零横波速度来改进qP-qSV波耦合波动方程, 从而保证成像的稳定性, DUVENECK等[104]提出一种新的稳定波动方程。第二类是完全不含伪横波干扰的TI介质纯qP波控制方程。DU等[105]在平方根近似和弱各向异性近似条件下, 推导出一种TTI介质时间-波数域纯qP波解耦方程, 该方程在运动学上较为准确, 但是在时空域为拟微分形式, 导致数值计算存在困难。为此, SONG等[106]、LIU等[107]提出了时空域最佳分离近似(OSA)方法, 该方法虽然计算精度高, 但成本也非常高; PESTANA等[108]提出了快速展开法(REM), 该方法在VTI介质逆时偏移成像中得到了广泛应用; ZHAN等[109]提出了基于REM的TTI介质混合有限差分-伪谱法, 该方法需要进行多次傅里叶变换和插值, 计算效率很低; CHU等[110]推导了一种不含分数阶算子的时空域纯qP波近似方程。针对起伏地表条件下的TI介质RTM成像, 需要在校正各向异性的同时应用起伏地表成像方法。兰海强等[87]实现了VTI介质起伏地表地震波数值模拟, 李庆洋[111]基于贴体网格实现了起伏地表TTI介质地震波正演模拟。

利用如图 18所示的高斯山峰、山谷模型进行了各向异性介质地震波数值模拟测试。高斯山峰、山谷模型的高程表达式为:

$ \begin{array}{l} y = 0 - 200\exp \left[ {\frac{{ - {{(x - 800)}^2}}}{{{{110}^2}}}} \right] + \\ 200\exp \left[ {\frac{{ - {{(x - 2200)}^2}}}{{{{110}^2}}}} \right]\;\;x \in [0,3000] \end{array} $ (4)
图 18 高斯山峰、山谷模型速度场(a)及其网格剖分(b)

模型从上至下分为两层, 模型参数分别为C110=10.0GPa, C130=2.5GPa, C330=6.0GPa, C550=2.0GPa, C110=16.0GPa, C130=4.48GPa, C330=16.0GPa, C550=5.76GPa。模型网格大小为301×201, 网格间距为10m。震源为主频20Hz的雷克子波, 时间采样间隔0.6ms, 震源位于(1500m, 10m)处。图 19a图 19b分别为模拟得到的VTI介质和TTI介质0.45s时刻的水平分量波场快照。从图 19可以看出, 在qSV波三叉区, 山峰、山谷相当于二次震源, 所有经过山峰、山谷的地震波都会激发出反射波和转换波, 极大地丰富了地震波场, 增加了波场复杂性。总之, 起伏地表的存在使得地震波在地下各向异性介质中的传播模式发生了很大变化, 使波场变得更加复杂。

图 19 高斯山峰、山谷模型0.45s时刻的水平分量波场快照 a VTI介质; b TTI介质
3.2.7 起伏地表最小二乘逆时偏移

常规起伏地表逆时偏移虽然可以得到较正确的地下成像结果, 满足构造成像的需求, 但仍然存在如下问题:①偏移噪声大。Laplacian滤波虽然能有效去除低频噪声, 但去除不彻底, 且引入了高频噪声。②偏移剖面上反射同相轴中间能量强、两侧能量较弱, 即振幅均衡性不佳。③由于地下照明强度随深度的增加而减弱, 因而RTM结果深部能量较弱, 振幅保真性差。采用起伏地表LSRTM算法可以有效解决常规RTM存在的问题。

近些年来, 最小二乘偏移因其具有信噪比高、保真性好、分辨率高等优点受到国内外学者的广泛关注, 其中最小二乘逆时偏移的研究最为广泛, 多震源最小二乘逆时偏移可以有效提高计算效率。为了压制多震源串扰噪声, 国内外学者提出了多种编码技术, 主要包括极性编码[112-113]、振幅编码[113]、分频编码[114]、平面波编码[115-116]等。2016年, HOU等[117]提出了基于扩展模型空间的最小二乘逆时偏移, 降低了最小二乘偏移对偏移速度的依赖性。随后, 徐凯等[118]和ZHU等[119]分别提出了校正粘弹性和各向异性的最小二乘逆时偏移方法。起伏地表最小二乘逆时偏移需要将声波方程转化到曲坐标系下, 根据Born近似理论得到迭代公式, 通过不断迭代得到成像结果。

利用3.2.3节图 12a所示含高斯山峰、山谷的洼陷模型进行了起伏地表LSRTM测试。观测系统和计算参数与3.2.3节相同。采用震源极性编码方式将101炮地震数据组合成一个超道集, 使其计算量相当于单炮情形, 从而大大缓解了计算的需求。图 20a为反射系数模型, 图 20b为基于相位编码的起伏地表LSRTM得到的第30次迭代结果。可以看出, 地下构造清晰可见, 相对起伏地表RTM结果(图 12c)在振幅保真性、均衡性、压制低频噪声等方面都有了较大改善, 但由编码引入的高频串扰噪声也清晰可见。图 20c为起伏地表LSRTM得到的第80次迭代结果, 该结果与理论反射系数模型(图 20a)非常接近, 相对第30次迭代结果(图 20b)有效压制了串扰噪声。

图 20 起伏地表LSRTM测试 a反射系数模型; b LSRTM第30次迭代结果; c LSRTM第80次迭代结果

近年来, 弹性波最小二乘逆时偏移也是一个热点方向。为了实现多分量数据成像, JIN等[120]建立了能够对纵横波进行同时反演的弹性反演理论。QU等[121]提出了基于波场分离技术的弹性波最小二乘逆时偏移, 对纵横波速度及密度分量进行分别成像, 并提出了基于波场分离的起伏地表弹性波最小二乘逆时偏移[122]

4 起伏地表速度反演方法

速度分析的准确性对于地震成像非常重要, 叠前深度偏移对速度分析的要求更高, 常规的速度分析方法对于起伏地表条件下的地震成像显得力不从心, 为此, 国内外学者研究了针对叠前深度偏移的速度反演方法。主要包括两类:第一类是起伏地表层析速度反演方法, 第二类是起伏地表全波形反演方法。两种方法均可以独立地获取叠前偏移速度场, 也可以进行联合, 实现多级优化反演、层析反演用以获取较高质量的低波数速度信息, 波形反演用以恢复速度场中的高波数成分, 从而得到优化的叠前偏移速度场。

4.1 起伏地表层析反演

针对起伏地表条件, 基于角道集的成像域走时层析是行之有效的速度反演方法。ZHANG等[123]和秦宁等[124-125]研究了基于角道集的层析速度反演方法, 给出了起伏地表成像域层析速度反演流程(图 21)[124]。首先利用高斯束偏移方法提取角度域共成像点道集, 由SAVA等[42]提出的成像条件将不同偏移张角对应的成像值累加到共成像点道集所对应的角度范围内即可; 然后根据实际的层析目标地质体, 选择最合适的参数化形式, 以最低的计算成本得到精确的层析速度。走时残差的拾取是成像域走时层析速度反演的关键步骤, 基于角道集能够获得精度较高的走时残差。走时残差和深度残差转换关系式为:

$ \Delta \mathit{\boldsymbol{t}} = 2s\Delta \mathit{\boldsymbol{z}}\cos \alpha \cos \beta $ (5)
图 21 起伏地表成像域层析速度反演流程[124]

式中:Δt为走时残差向量; Δz为深度残差向量; s为成像点处的局部慢度值; α为反射层倾角; β为射线入射角, 对应角度域共成像点道集的角度。图 22为走时残差与深度残差转化关系示意图[124]。另外, 在反演过程中, 可以通过加入正则化和井约束来提高层析反演的稳定性和精度。图 23a图 23b为起伏地表模型及其炮集记录, 图 23c图 23d为层析速度场及深度偏移剖面[126], 可以看出, 起伏地表模型得到了较为清晰的成像。

图 22 走时残差与深度残差转化关系[124]
图 23 起伏地表模型层析反演 a起伏地表速度模型; b炮集记录; c层析速度场; d层析后的深度偏移剖面[126]
4.2 起伏地表全波形反演

在全波形反演过程中, 通常采用从低频到高频逐步反演的多尺度反演策略。Wiener滤波器能够将源信号转化为十分接近目标信号的形式, BOONYASIRIWAT等[126]提出了基于Wiener滤波器的高效时间域全波形反演方法, 可以将地震数据的频率转化为较低频率, 以实现多尺度反演。考虑到计算量和空间采样问题, 采用变网格可以很好地适应起伏地表的情况, 因此, 时间域多尺度双变网格全波形反演可以较好地解决起伏地表的速度反演问题。采用如图 24a所示的起伏地表模型对时空双变网格全波形反演方法进行了测试, 图 24b图 24c分别为时空双变网格FWI和常规粗网格FWI经过50次迭代得到的反演速度[127]。可以看出, 常规粗网格FWI的结果在起伏地表处出现了较大误差, 而时空双变网格FWI的结果较为准确, 深部速度界面和起伏地表附近的速度值都得到了正确的反演。与时间域多尺度双变网格全波形反演类似, 在频率域同样可以实现多尺度双变网格全波形反演。

图 24 50次迭代的反演速度 a起伏地表模型; b双变网格FWI; c常规粗网格FWI[127]

常规的矩形网格模拟起伏地表时因其阶梯状离散会引起虚假散射和绕射, 因此难以准确模拟起伏地表的波场。针对这一问题, 我们采用曲网格来进行离散化。垂向曲网格-矩形网格耦合方法是将起伏地表附近的区域网格化为曲网格, 而深层区域则网格化为矩形网格, 可以较好地模拟起伏地表的波场。然而, 垂向曲网格坐标变换技术只能在垂向上进行网格映射[49], 使该方法受到严格的稳定性条件限制, 因此无法稳定地模拟地震波在剧烈起伏地表中的波场。QU等[128]根据地表附近正交曲网格的坐标变换方程[82], 提出了曲线系统下的正交曲网格-矩形网格耦合方法, 用以反演起伏地表条件下的速度。采用SEG逆掩断层模型(图 25)对正交曲网络-矩形网格方法进行了测试, 图 26为最终速度反演结果[128]。该反演结果表明, 正交曲网格-矩形网格耦合方法可以在复杂近地表情况下较好地反演出速度值。

图 25 SEG逆掩断层模型 a vP; b vS[128]
图 26 SEG逆掩断层模型反演结果 a vP; b vS[128]
5 应用中的瓶颈问题及技术方案

起伏地表射线类方法能够很好地适应起伏地表, 但主要反映地下构造的运动学信息, 难以对地下复杂构造的动力学信息进行准确的刻画, 且该方法难以对高陡构造等复杂构造进行准确成像。因此, 在研究起伏地表射线类方法时, 要尽可能考虑地下构造的复杂性。①围绕高斯束格林函数构建方式, 考虑模型介质水平及纵向速度变化对束形态的影响, 发展基于速度驱动的起伏地表自适应高斯束偏移方法。同时考虑浅部及深部构造成像分辨率, 引入声学中的“菲涅尔带”构建思想, 发展起伏地表投影菲涅尔束偏移方法。②考虑复杂地下构造的粘弹性和各向异性, 发展起伏地表粘弹性及各向异性射线偏移方法。③引入反演思想, 发展起伏地表最小二乘Kirchhoff偏移或高斯束偏移。

起伏地表单程波成像方法在成像精度上不如逆时偏移方法, 在计算效率和对起伏地表的适应性上不如射线类方法, 因此后续的研究重点是突出两者的优点, 弱化两者的缺点。

基于双程波的逆时偏移在近十年来得到了飞速发展。该方法从理论上可以对任意复杂构造进行精确成像, 近年来大量的应用实例也证明该方法比射线类和单程波类成像方法得到更好的成像结果, 但该方法对剧烈起伏地表的适应性较差, 需要进行特殊处理。起伏地表逆时偏移的关键是起伏地表波场延拓算子的构建方法, 从求解波动方程的角度可将其分为有限差分法、有限元法、伪谱法、谱元法等。其中, 有限差分法的突出优势是计算速度快、占用内存少, 但其对剧烈起伏地表的适应性最差。制约起伏地表逆时偏移方法的瓶颈问题主要有两个:①对速度的依赖性强; ②计算量庞大, 特别是近些年基于反演的最小二乘逆时偏移, 对地下构造粘弹性和各向异性的引入更是给该方法带来了巨大的计算负担。起伏地表逆时偏移方法对速度的依赖性强, 但在山前带等剧烈起伏地表、复杂构造区域, 精确的速度建模非常困难, 为此, 地球物理学家们研究了多种对起伏地表进行准确速度反演的方法。然而, 任意一种速度建模方法都不能获得准确的速度, 需要将多种速度建模方法配合使用, 充分发挥速度分析、层析反演、波形反演等多种速度建模方法的优势, 并采用多信息融合的浅—中—深层、低—中—高频逐级优化的反演策略, 提高起伏地表速度反演的精度和效率。

在计算效率提升方面, 可以采取以下策略:①利用各种近似, 在保证计算精度的前提下最大限度地降低计算量。比如校正各向异性时, 使用弱各向异性假设, 校正粘弹性时, 只补偿对成像影响相对大的振幅衰减的影响, 忽略相位频散的影响。②对起伏地表波场模拟算子进行优化, 提高波场模拟的效率。③充分利用GPU/CPU协同加速、MPI/Openmp并行加速及波场区域分解等技术。④在起伏地表最小二乘逆时偏移中, 利用编码技术压缩地震数据, 并发展提高最小二乘逆时偏移收敛速度的反演算法等。

6 结论与展望

地震成像是一门不断发展、不断深入的学科。成像技术不断发展的最终目的是使地下复杂构造精确成像。起伏地表直接偏移技术不仅具有处理灵活、计算效率高的特点, 而且能对复杂构造进行精确成像, 消除常规偏移方法中的噪声影响, 对于我国西南地区复杂山地、丘陵构造的成像有着重要意义。然而, 目前起伏地表直接偏移技术尚处于实验探索阶段, 存在很多问题需要解决。地下岩性及储层精细描述的需求必然会带动起伏地表直接偏移技术的发展, 笔者认为起伏地表直接偏移技术的发展方向有以下两个方面。

在提高偏移成像精度方面, 此前的成像方法都是基于各种假设使模型简单化, 但地下介质受构造运动的影响其性质多种多样, 需要对复杂介质进行不断探索才能实现起伏地表的精确成像。偏移速度的准确性也会对成像质量有很大的影响, 如何通过反演得到准确的速度也是目前需要深入探索的方向。

在提高成像精度的同时, 计算量的增加降低了工作效率也是一个需要考虑的实际问题。为了消除起伏地表的影响, 对起伏地表进行网格剖分的方式会影响成像质量。双变网格计算速度快、效率高, 但无法消除起伏地表噪声的影响。曲网格成像质量高, 但进行坐标变换无疑会增加计算量, 降低工作效率。

综上所述, 如何在提高复杂构造成像效果的同时, 保证算法的灵活简便, 降低计算成本, 是广大地球物理学家需要继续探索的方向。

参考文献
[1]
刘定进, 刘志成, 蒋波. 面向复杂山前带的深度域地震成像处理研究[J]. 石油物探, 2016, 55(1): 49-59.
LIU D J, LIU Z C, JIANG B. The processing workflow of depth domain imaging facing the complex piedmont belt[J]. Geophysical Prospecting for Petroleum, 2016, 55(1): 49-59. DOI:10.3969/j.issn.1000-1441.2016.01.007
[2]
田文辉, 李振春, 张辉. 直接下延法波动方程叠前深度偏移[J]. 石油物探, 2006, 45(5): 447-451.
TIAN W H, LI Z C, ZHANG H. Direct decentralized wave equation prestack depth migration[J]. Geophysical Prospecting for Petroleum, 2006, 45(5): 447-451. DOI:10.3969/j.issn.1000-1441.2006.05.003
[3]
MILLINAZZO F A, ZALA C A, BROOKE G H. Square-root approximations for parabolic equation algorithms[J]. The Journal of the Acoustical Society of America, 1997, 101(2): 760-766. DOI:10.1121/1.418038
[4]
李录明, 罗省贤. 波场延拓表层模型校正[J]. 石油地球物理勘探, 2001, 36(5): 572-583.
LI L M, LUO X X. Surface model correction by wave-field continuation[J]. Oil Geophysical Prospecting, 2001, 36(5): 572-583. DOI:10.3321/j.issn:1000-7210.2001.05.008
[5]
秦宁, 王延光, 杨晓东, 等. 非水平地表高斯束叠前深度偏移及山前带应用实例[J]. 石油地球物理勘探, 2017, 52(1): 81-86.
QIN N, WANG Y G, YANG X D, et al. Gaussian beam prestack depth migration for undulating-surface area in piedmont zone[J]. Oil Geophysical Prospecting, 2017, 52(1): 81-86.
[6]
黄建平, 杨继东, 李振春, 等. 基于有效邻域波场近似的起伏地表保幅高斯束偏移[J]. 地球物理学报, 2016, 59(6): 2245-2256.
HUANG J P, YANG J D, LI Z C, et al. An amplitude-preserved Gaussian beam migration based on wave field approximation in effective vicinity under irregular topographical conditions[J]. Chinese Journal of Geophysics, 2016, 59(6): 2245-2256.
[7]
陈林谦.基于起伏地表单程波波动方程的叠前深度偏移成像技术应用与研究[D].北京: 中国石油大学, 2017
CHEN L Q.Application and method study of pre-stack one-wave equation depth migration based on irregular topography[D].Beijing: China University of Petroleum, 2017 http://cdmd.cnki.com.cn/Article/CDMD-11414-1019805584.htm
[8]
段心标, 王华忠, 白英哲, 等. 基于GPU的三维起伏地表单程波叠前深度偏移[J]. 石油物探, 2016, 55(2): 223-230.
DUAN X B, WANG H Z, BAI Y Z, et al. 3D one-way wave equation pre-stack depth migration from topography based on the acceleration of GPU[J]. Geophysical Prospecting for Petroleum, 2016, 55(2): 223-230. DOI:10.3969/j.issn.1000-1441.2016.02.008
[9]
张宇, 徐升, 张关泉, 等. 真振幅全倾角单程波方程偏移方法[J]. 石油物探, 2007, 46(6): 582-587, 642.
ZHANG Y, XU S, ZHANG G Q, et al. True amplitude full tilt one-way wave equation migration method[J]. Geophysical Prospecting for Petroleum, 2007, 46(6): 582-587, 642. DOI:10.3969/j.issn.1000-1441.2007.06.010
[10]
刘红伟, 刘洪, 李博, 等. 起伏地表叠前逆时偏移理论及GPU加速技术[J]. 地球物理学报, 2011, 54(7): 1883-1892.
LIU H W, LIU H, LI B, et al. Pre-stack reverse time migration for rugged topography and GPU acceleration technology[J]. Chinese Journal of Geophysics, 2011, 54(7): 1883-1892. DOI:10.3969/j.issn.0001-5733.2011.07.022
[11]
王延光, 匡斌. 起伏地表叠前逆时深度偏移与并行实现[J]. 石油地球物理勘探, 2012, 47(2): 266-271.
WANG Y G, KUANG B. Pre-stack reverse-time depth migration on rugged topography and parallel computation realization[J]. Oil Geophysical Prospecting, 2012, 47(2): 266-271.
[12]
WIGGINS J W. Kirchhof fintegral extrapolation and migration of nonplanar data[J]. Geophysics, 1984, 49(8): 1239-1248. DOI:10.1190/1.1441752
[13]
GRAY S H. Gaussian beam migration of common-shot records[J]. Geophysics, 2005, 70(4): S71-S77. DOI:10.1190/1.1988186
[14]
JAGER C, HERTWECK T, SPINER M. True-amplitude Kirchhoff migration from topography[J]. Expanded Abstracts of 73rd Annual Internat SEG Mtg, 2003, 909-913.
[15]
岳玉波, 李振春, 钱忠平, 等. 复杂地表条件下保幅高斯束偏移[J]. 地球物理学报, 2012, 55(4): 1376-1383.
YUE Y B, LI Z C, QIAN Z P, et al. Amplitude-preserved Gaussian beam migration under complex topographic conditions[J]. Chinese Journal of Geophysics, 2012, 55(4): 1376-1383. DOI:10.6038/j.issn.0001-5733.2012.04.033
[16]
杨继东, 黄建平, 王欣, 等. 复杂地表条件下叠前菲涅尔束偏移方法[J]. 地球物理学报, 2015, 58(10): 3758-3770.
YANG J D, HUANG J P, WANG X, et al. Pre-stack Fresnel beam migration method under complex topographic conditions[J]. Chinese Journal of Geophysics, 2015, 58(10): 3758-3770. DOI:10.6038/cjg20151026
[17]
黄建平, 袁茂林, 段新意, 等. 一种解耦的起伏地表弹性波高斯束偏移方法[J]. 石油地球物理勘探, 2015, 50(3): 460-468.
HUANG J P, YUAN M L, DUAN X Y, et al. Decoupled elastic Gaussian beam migration for rugged topography[J]. Oil Geophysical Prospecting, 2015, 50(3): 460-468.
[18]
BEASLEY C J, LYNN W. The zero velocity layer:migration from irregular surfaces[J]. Geophysics, 1992, 57(11): 1435-1443. DOI:10.1190/1.1443211
[19]
RESHEF M. Depth migration from irregular surfaces with the depth extrapolation methods[J]. Geophysics, 1991, 56(1): 119-122. DOI:10.1190/1.1442947
[20]
雷秀丽, 李振春, 孔雪, 等. 起伏地表条件下基于复Pade逼近的叠前深度偏移[J]. 地球物理学进展, 2012, 27(5): 2100-2106.
LEI X L, LI Z C, KONG X, et al. Pre-stack depth migration using complex Pade approximations for irregular surfaces[J]. Progress in Geophysics, 2012, 27(5): 2100-2106.
[21]
何英, 王华忠, 马在田, 等. 复杂地形条件下波动方程叠前深度成像[J]. 勘探地球物理进展, 2002, 25(3): 14-19.
HE Y, WANG H Z, MA Z T, et al. Pre-stack wave equation depth migration for irregular topography[J]. Progress in Exploration Geophysics, 2002, 25(3): 14-19.
[22]
程玖兵, 王华忠, 马在田. 频率-空间域有限差分法叠前深度偏移[J]. 地球物理学报, 2001, 44(3): 390-395.
CHENG J B, WANG H Z, MA Z T. Pre-stack depth migration with finite difference method in frequency-space domain[J]. Chinese Journal of Geophysics, 2001, 44(3): 390-395.
[23]
王成祥, 赵波, 张关泉. 基于起伏地表的混合法叠前深度偏移[J]. 石油地球物理勘探, 2002, 37(3): 219-223.
WANG C X, ZHAO B, ZHANG G Q. Rugged surface oriented hybrid pre-stack depth migration[J]. Oil Geophysical Prospecting, 2002, 37(3): 219-223. DOI:10.3321/j.issn:1000-7210.2002.03.004
[24]
叶月明, 李振春, 仝兆岐, 等. 双复杂条件下带误差补偿的频率-空间域有限差分法叠前深度偏移[J]. 地球物理学进展, 2008, 23(1): 136-145.
YE Y M, LI Z C, TONG Z Q, et al. Xwfd pre-stack depth migration based on dual-complexity with error compensation correction[J]. Progress in Geophysics, 2008, 23(1): 136-145.
[25]
叶月明, 李振春, 仝兆岐, 等. 双复杂介质条件下频率-空间域有限差分法保幅偏移[J]. 地球物理学报, 2008, 51(5): 1511-1519.
YE Y M, LI Z C, TONG Z Q, et al. Preserved-amplitude migration with dual-complexity[J]. Chinese Journal of Geophysics, 2008, 51(5): 1511-1519. DOI:10.3321/j.issn:0001-5733.2008.05.025
[26]
WHITMORE N D. Iterative depth migration by backward time propagation[J]. Expanded Abstracts of 53rd Annual Internat SEG Mtg, 1983, 382-385.
[27]
刘红伟, 李博, 刘洪, 等. 地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报, 2010, 53(7): 1725-1733.
LIU H W, LI B, LIU H, et al. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J]. Chinese Journal of Geophysics, 2010, 53(7): 1725-1733. DOI:10.3969/j.issn.0001-5733.2010.07.024
[28]
XU S, ZHANG Y, TANG B. 3D common image gathers from reverse time migration[J]. Expanded Abstracts of 80th Annual Internat SEG Mtg, 2010, 3257-3261.
[29]
KAELIN B, GUITTON A. Imaging condition for reverse time migration[J]. Expanded Abstracts of 76th Annual Internat SEG Mtg, 2006, 2594-2598.
[30]
刘红伟, 刘洪, 邹振. 地震叠前逆时偏移中的去噪与存储[J]. 地球物理学报, 2010, 53(9): 2171-2180.
LIU H W, LIU H, ZOU Z. The problems of denoise and storage in seismic reverse time migration[J]. Chinese Journal of Geophysics, 2010, 53(9): 2171-2180. DOI:10.3969/j.issn.0001-5733.2010.09.017
[31]
BULCÃO A, FILHO D M S, MANSUR W J. Improved quality of depth images using reverse time migration[J]. Expanded Abstracts of 77th Annual Internat SEG Mtg, 2007, 2407-2411.
[32]
GUITTON A, KAELIN B, BIONDI B. Least-squares attenuation of reverse-time-migration artifacts[J]. Geophysics, 2008, 73(1): 19-23.
[33]
杜启振, 朱钇同, 张明强, 等. 叠前逆时深度偏移低频噪声压制策略研究[J]. 地球物理学报, 2013, 56(7): 2391-2401.
DU Q Z, ZHU Y T, ZHANG M Q, et al. A study on the strategy of low wave number noise suppression for pre-stack reverse-time depth migration[J]. Chinese Journal of Geophysics, 2013, 56(7): 2391-2401.
[34]
刘阳, 赵虎, 尹成, 等. 起伏地表条件下的弹性波叠前逆时偏移[J]. 煤田地质与勘探, 2019, 47(1): 181-186.
LIU Y, ZHAO H, YIN C, et al. Elastic wave pre-stack reverse time migration on undulating surface[J]. Coal Geology & Exploration, 2019, 47(1): 181-186. DOI:10.3969/j.issn.1001-1986.2019.01.028
[35]
曲英铭, 黄建平, 李振春, 等. 分层坐标变换法起伏自由地表弹性波叠前逆时偏移[J]. 地球物理学报, 2015, 58(8): 2896-2911.
QU Y M, HUANG J P, LI Z C, et al. Elastic wave modeling and pre-stack reverse time migration of irregular free-surface based on layered mapping method[J]. Chinese Journal of Geophysics, 2015, 58(8): 2896-2911.
[36]
孙小东, 李振春, 王小六. 三角网格有限差分法叠前逆时偏移方法研究[J]. 地球物理学进展, 2012, 27(5): 2077-2083.
SUN X D, LI Z C, WANG X L. Pre-stack reverse-time migration using a finite difference method based on triangular grids[J]. Progress in Geophysics, 2012, 27(5): 2077-2083.
[37]
HILL N R. Gaussian beam migration[J]. Geophysics, 1990, 55(11): 1416-1428. DOI:10.1190/1.1442788
[38]
HILL N R. Pre-stack Gaussian-beam depth migration[J]. Geophysics, 2001, 66(4): 1240-1250. DOI:10.1190/1.1487071
[39]
SUN Y H, QIN F H, CHECKLES S, et al. 3-D pre-stack Kirchhoff beam migration for depth imaging[J]. Geophysics, 2000, 65(5): 1592-1603. DOI:10.1190/1.1444847
[40]
COCKSHOTT I. Specular beam migration:a low cost 3-D pre-stack depth migration[J]. Expanded Abstracts of 76th Annual Internat SEG Mtg, 2006, 2634-2638.
[41]
TING C O, WANG D L. Controlled beam migration applications in Gulf of Mexico[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 368-372.
[42]
SAVA P C, FOMEL S. Angle-domain common-image gathers by wavefield continuation methods[J]. Geophysics, 2003, 68(3): 1065-1074. DOI:10.1190/1.1581078
[43]
朱绪峰.单程波方程的保幅偏移及其角度域成像[D].北京: 中国石油大学, 2008
ZHU X F.One-way wave equation preserved-amplitude migration and imaging in angles domain[D].Beijing: China University of Petroleum, 2008 http://cdmd.cnki.com.cn/Article/CDMD-10425-2008199525.htm
[44]
BERRYHILL J R. Wave equation datuming[J]. Geophysics, 1979, 44(8): 1329-1344. DOI:10.1190/1.1441010
[45]
BERRYHILL J J R. Wave equation datuming before stack[J]. Geophysics, 1984, 49(11): 2064-2066. DOI:10.1190/1.1441620
[46]
WIDESS M B. Effect of surface topography on seismic mapping[J]. Geophysics, 1946, 11(3): 362-372. DOI:10.1190/1.1437260
[47]
ALTERMAN Z, KARAL F C. Propagation of elastic waves in layered media by finite difference methods[J]. Bulletin of the Seismological Society of America, 1968, 58(1): 367-398.
[48]
VIRIEUX J. P-SV wave propagation in heterogeneous media:velocity-stress finite-difference method[J]. Geophysics, 1986, 51(11): 1933-1942.
[49]
JASTRAM C, TESSMER E. Elastic modelling on a grid with vertically varying spacing[J]. Geophysical Prospecting, 2010, 42(4): 357-370.
[50]
TESSMER E, KOSLOFF D. 3-D elastic modeling with surface topography by a Chebychev spectral method[J]. Geophysics, 1994, 59(3): 464-473. DOI:10.1190/1.1443608
[51]
董良国. 复杂地表条件下地震波传播数值模拟[J]. 勘探地球物理进展, 2005, 28(3): 187-194.
DONG L G. Numerical simulation of seismic wave propagation under complex surface conditions[J]. Progress in Exploration Geophysics, 2005, 28(3): 187-194.
[52]
褚春雷.非规则网格有限差分法声波方程地震正演模拟及逆时偏移[D].青岛: 中国海洋大学, 2003
CHU C L.Seismic modeling and reserse-time migration with a finite-difference method on irregular grids[D].Qingdao: Ocean University of China, 2003 http://cdmd.cnki.com.cn/Article/CDMD-10423-2003093453.htm
[53]
THOMPSON J F, THAMES F C, MASTIN C W. Automatic numerical generation of body-fitted curvilinear coordinates for a field containing any number of arbitrary two-dimensional bodies[J]. Journal of Computational Physics, 1974, 15(3): 299-319.
[54]
孙建国, 蒋丽丽. 用于起伏地表条件下地球物理场数值模拟的正交曲网格生成技术[J]. 石油地球物理勘探, 2009, 44(4): 494-500, 528.
SUN J G, JIANG L L. Orthogonal curvilinear grid generation technique used for numeric simulation of geophysical fields in undulating surface condition[J]. Oil Geophysical Prospecting, 2009, 44(4): 494-500, 528. DOI:10.3321/j.issn:1000-7210.2009.04.020
[55]
KOMATITSCH D, VILOTTE J P. The spectral element method:an efficient tool to simulate the seismic response of 2D and 3D geological structures[J]. Bulletin of the Seismological Society of America, 1998, 88(2): 368-392.
[56]
CHE C X, WANG X M, LIN W J. The Chebyshev spectral element method using staggered predictor and corrector for elastic wave simulations[J]. Applied Geophysics, 2010, 7(2): 174-184, 195. DOI:10.1007/s11770-010-0242-9
[57]
BAYSAL E, DAN D K, SHERWOOD J W C. Reverse time migration[J]. Geophysics, 1983, 48(11): 1514-1524. DOI:10.1190/1.1441434
[58]
DAN I R M. Reversed time migration in spatial frequency domain[J]. Geophysics, 1983, 48(5): 627-635. DOI:10.1190/1.1441493
[59]
MUFTI I R. Large-scale three-dimensional seismic models and their interpretive significance[J]. Geophysics, 1990, 55(9): 1166-1182. DOI:10.1190/1.1442933
[60]
LESAGE A C, ZHOU H B, CELA J M, et al. 3D reverse-time migration with hybrid finite difference:pseudo spectral method[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 2258-2261.
[61]
JONES I F, GOODWIN M C, BERRANGER I D, et al. Application of anisotropic 3D reverse time migration to complex North Sea imaging[J]. Expanded Abstracts of 77th Annual Internat SEG Mtg, 2007, 2140-2144.
[62]
CHANG W, MCMECHAN G. 3-D elastic pre-stack reverse-time depth migration[J]. Geophysics, 1994, 59(4): 597-609. DOI:10.1190/1.1443620
[63]
YAN J, SAVA P. Isotropic angle-domain elastic reverse-time migration[J]. Geophysics, 2008, 73(6): S229-S239. DOI:10.1190/1.2981241
[64]
李振春, 雍鹏, 黄建平, 等. 基于矢量波场分离弹性波逆时偏移成像[J]. 中国石油大学学报(自然科学版), 2016, 40(1): 42-48.
LI Z C, YONG P, HUANG J P, et al. Elastic wave reverse time migration based on vector wavefield seperation[J]. Journal of China University of Petroleum(Edition of Natural Science), 2016, 40(1): 42-48. DOI:10.3969/j.issn.1673-5005.2016.01.006
[65]
周熙焱, 常旭, 王一博, 等. 基于纵横波解耦的三维弹性波逆时偏移[J]. 地球物理学报, 2018, 61(3): 1038-1052.
ZHOU X Y, CHANG X, WANG Y B, et al. 3D elastic reverse time migration based on P-and S-wave decoupling[J]. Chinese Journal of Geophysics, 2018, 61(3): 1038-1052.
[66]
SUN R, CHOW J, CHEN K J. Phase correction in separating P-and S-waves in elastic data[J]. Geophysics, 2001, 66(5): 1515-1518. DOI:10.1190/1.1487097
[67]
王维红, 张伟, 石颖, 等. 基于波场分离的弹性波逆时偏移[J]. 地球物理学报, 2017, 60(7): 2813-2824.
WANG W H, ZHANG W, SHI Y, et al. Elastic reverse time migration based on wavefield separation[J]. Chinese Journal of Geophysics, 2017, 60(7): 2813-2824.
[68]
DELLINGER J, ETGEN J. Wave-field separation in two-dimensional anisotropic media[J]. Geophysics, 1990, 55(7): 914-919. DOI:10.1190/1.1442906
[69]
BALCH A H, ERDEMIR C. Sign-change correction for pre-stack migration of P-S converted wave reflections[J]. Geophysical Prospecting, 1994, 42(6): 637-663. DOI:10.1111/j.1365-2478.1994.tb00233.x
[70]
DU Q, ZHU Y, BA J. Polarity reversal correction for elastic reverse time migration[J]. Geophysics, 2012, 77(2): S31-S41. DOI:10.1190/geo2011-0348.1
[71]
ROBIN P, FLETCHER E A. Suppressing unwanted internal reflections in pre-stack reverse-time migration[J]. Geophysics, 2006, 71(6): E79-E82. DOI:10.1190/1.2356319
[72]
LIU F Q. An anti-dispersion wave equation for modeling and reverse-time migration[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 2277-2281.
[73]
SUH S Y, CAI J. Reverse-time migration by fan filtering plus wavefield decomposition[J]. Expanded Abstracts of 79th Annual Internat SEG Mtg, 2009, 2804-2808.
[74]
李博, 刘志成, 李小爱, 等. 基于复数域波场分解的保幅逆时偏移成像方法[J]. 石油物探, 2019, 58(2): 237-244.
LI B, LIU Z C, LI X A, et al. Wavefield decomposition in complex domain-based amplitude-preserved reverse time migration[J]. Geophysical Prospecting for Petroleum, 2019, 58(2): 237-244. DOI:10.3969/j.issn.1000-1441.2019.02.009
[75]
YOON K, MARFURT K J. Reverse-time migration using the Poynting vector[J]. Exploration Geophysics, 2006, 37(1): 102-107.
[76]
DU Q, GUO C F, ZHAO Q, et al. Vector-based elastic reverse time migration based on scalar imaging condition[J]. Geophysics, 2017, 82(2): S111-S127. DOI:10.1190/geo2016-0146.1
[77]
JASTRAM C, BEHLE A. Acoustic modeling on a vertically varying grid[J]. Geophysical Prospecting, 1992, 40(2): 157-169. DOI:10.1111/j.1365-2478.1992.tb00369.x
[78]
马光克, 李洋森, 孙万元, 等. 可变网格高阶有限差分法逆时偏移研究[J]. 石油物探, 2016, 55(5): 728-736.
MA G K, LI Y S, SUN W Y, et al. Acoustic pre-stack reverse time migration using variable grid finite-difference method[J]. Geophysical Prospecting for Petroleum, 2016, 55(5): 728-736. DOI:10.3969/j.issn.1000-1441.2016.05.012
[79]
TAE-SEOB K. Finite-difference seismic simulation combining discontinuous grids with locally variable timesteps[J]. Bulletin of the Seismological Society of America, 2004, 94(1): 207-219. DOI:10.1785/0120030080
[80]
HUANG J P, QU Y M, LI Q Y, et al. Variable-coordinate forward modeling of irregular surface based on dual-variable grid[J]. Applied Geophysics, 2015, 12(1): 101-110, 123. DOI:10.1007/s11770-014-0476-2
[81]
李振春, 李庆洋, 黄建平, 等. 一种稳定的高精度双变网格正演模拟与逆时偏移方法[J]. 石油物探, 2014, 53(2): 127-136.
LI Z C, LI Q Y, HUANG J P, et al. A stable and high-precision dual-variable grid forward modeling and reverse time migration method[J]. Geophysical Prospecting for Petroleum, 2014, 53(2): 127-136. DOI:10.3969/j.issn.1000-1441.2014.02.001
[82]
FORNBERG B. The pseudo spectral method:accurate representation of interfaces in elastic wave calculations[J]. Geophysics, 1988, 53(5): 625-637. DOI:10.1190/1.1442497
[83]
HESTHOLM ST O, RUUD B O. 2D finite-difference viscoelastic wave modelling including surface topography[J]. Geophysical Prospecting, 2000, 48(2): 341-373. DOI:10.1046/j.1365-2478.2000.00185.x
[84]
鲁雁翔, 白超英. 起伏层状介质中曲线交错网格有限差分弹性波逆时偏移成像[J]. 地球物理学进展, 2019, 34(2): 605-614.
LU Y X, BAI C Y. Time reverse migration based on curvilinear stagger-grid finite difference method[J]. Progress in Geophysics, 2019, 34(2): 605-614.
[85]
ZHANG W, SHEN Y, ZHAO L. Three-dimensional anisotropic seismic wave modelling in spherical coordinates by a collocated-grid finite-difference method[J]. Geophysical Journal International, 2012, 188(3): 1359-1381. DOI:10.1111/j.1365-246X.2011.05331.x
[86]
APPELO D, PETERSSON N A. A stable finite difference method for the elastic wave equation on complex geometries with free surfaces[J]. Communications in Computational Physics, 2009, 5(1): 84-107.
[87]
兰海强, 刘佳, 白志明. VTI介质起伏地表地震波场模拟[J]. 地球物理学报, 2011, 54(8): 2072-2084.
LAN H Q, LIU J, BAI Z M. Wave-field simulation in VTI media with irregular free surface[J]. Chinese Journal of Geophysics, 2011, 54(8): 2072-2084. DOI:10.3969/j.issn.0001-5733.2011.08.014
[88]
CHENG J W, FAN N, ZHANG Y Y, et al. Irregular surface seismic forward modeling by a body-fitted rotated-staggered-grid finite difference method[J]. Applied Geophysics, 2018, 15(Z1): 420-431, 567.
[89]
李庆洋, 黄建平, 李振春. 起伏地表贴体全交错网格仿真型有限差分正演模拟[J]. 石油地球物理勘探, 2015, 50(4): 633-642.
LI Q Y, HUANG J P, LI Z C. Undulating surface body-fitted grid seismic modeling based on fully staggered-grid mimetic finite difference[J]. Oil Geophysical Prospecting, 2015, 50(4): 633-642.
[90]
THOMPSON J F, WARSI Z U A, MASTIN C W. Numerical grid generation:foundations and applications[M]. Amsterdam: North-Holland, 1985: 8-17.
[91]
QU Y, LI J. Q-compensated reverse time migration in viscoacoustic media including surface topography[J]. Geophysics, 2019, 84(4): S201-S217. DOI:10.1190/geo2018-0313.1
[92]
GUO S J, LI Z C, SUN X D, et al. Post-stack reverse-time migration using a finite difference method based on triangular grids[J]. Applied Geophysics, 2008, 5(2): 115-120. DOI:10.1007/s11770-008-0001-y
[93]
ZHANG X W, HAN L G, ZHANG F J, et al. An inverse Q-filter algorithm based on stable wavefeild continuation[J]. Applied Geophysicis, 2007, 4(4): 263-270. DOI:10.1007/s11770-007-0040-9
[94]
HARGREAVES N D, CALVERT A J. Inverse Q filtering by Fourier transform[J]. Geophysics, 1991, 56(4): 519-527. DOI:10.1190/1.1443067
[95]
TRAYNIN P, LIU J, REILLY J M. Amplitude and bandwidth recovery beneath gas zones using Kirchhoff pre-stack depth Q-migration[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 2412-2416.
[96]
刘喜武, 年静波, 刘洪. 基于广义S变换的吸收衰减补偿方法[J]. 石油物探, 2006, 45(1): 9-14.
LIU X W, NIAN J B, LIU H. Absorption attenuation compensation method based on generalized S transform[J]. Geophysical Prospecting for Petroleum, 2006, 45(1): 9-14. DOI:10.3969/j.issn.1000-1441.2006.01.002
[97]
ZHANG J, WU J, LI X. Compensation for absorption and dispersion in pre-stack migration:an effective Q approach[J]. Geophysics, 2013, 78(1): S1-S14.
[98]
[99]
陈雪, 韩立国, 杨贺龙, 等. 粘弹VTI介质单程波正演模拟与照明分析[J]. 石油物探, 2015, 54(6): 643-651.
CHEN X, HAN L G, YANG H L, et al. One-way wave equation forward modeling and illumination analysis in viscoelastic VTI media[J]. Geophysical Prospecting for Petroleum, 2015, 54(6): 643-651. DOI:10.3969/j.issn.1000-1441.2015.06.002
[100]
LIU X, CHEN J, ZHAO Z, et al. Q-compensated reverse time migration in viscoelastic media with irregular free surface[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 4473-4477.
[101]
TSVANKIN I. P-wave signatures and notation for transversely isotropic media[J]. Geophysics, 1996, 61(1): 467-483.
[102]
ALKHALIFAH T. Acoustic approximations for processing in transversely isotropic media[J]. Geophysics, 1998, 63(2): 623-631. DOI:10.1190/1.1444361
[103]
FLETCHER R P, DU X, FOWLER P J. Reverse time migration in tilted transversely isotropic (TTI) media[J]. Geophysics, 2009, 74(6): WCA179-WCA187. DOI:10.1190/1.3269902
[104]
DUVENECK E, MILCIK P, BAKKER P M, et al. Acoustic VTI wave equations and their application for anisotropic reverse-time migration[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 2186-2190.
[105]
DU X, BANCROFT J C, LINES L R. Reverse-time migration for tilted TI media[J]. Expanded Abstracts of 75th Annual Internat SEG Mtg, 2005, 1930-1934.
[106]
SONG X, FOMEL S, YING L X. Low rank finite-differences and low rank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation[J]. Geophysical Journal International, 2013, 193(2): 960-969. DOI:10.1093/gji/ggt017
[107]
LIU F Q, MORTON S A, JIANG S S, et al. Decoupled wave equations for P and SV waves in an acoustic VTI media[J]. Expanded Abstracts of 79th Annual Internat SEG Mtg, 2009, 2844-2848.
[108]
PESTANA R C, STOFFA P L. Time evolution of the wave equation using rapid expansion method[J]. Geophysics, 2010, 75(4): T121-T131. DOI:10.1190/1.3449091
[109]
ZHAN G, PESTANA R C, STOFFA P L. An efficient hybrid pseudo spectral/finite-difference scheme for solving TTI pure P-wave equation[J]. Journal of Geophysics and Engineering, 2013, 10(2): 1322-1327.
[110]
CHU C L, MACY B K, ANNO P D. Pure acoustic wave propagation in transversely isotropic media by the pseudospectral method[J]. Geophysical Prospecting, 2013, 61(3): 556-567. DOI:10.1111/j.1365-2478.2012.01077.x
[111]
李庆洋.复杂介质地震波正演模拟与最小二乘偏移研究[D].青岛: 中国石油大学, 2017
LI Q Y.Study on seismic forward modeling and least-squares migration in complex media[D].Qingdao: China University of Petroleum, 2017 http://cdmd.cnki.com.cn/Article/CDMD-10425-1019836689.htm
[112]
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
[113]
HU J, WANG H, FANG Z, et al. Efficient amplitude encoding least-squares reverse time migration using cosine basis[J]. Geophysical Prospecting, 2016, 64(6): 1483-1497. DOI:10.1111/1365-2478.12356
[114]
HUANG Y, SCHUSTER G T. Multisource least-squares migration of marine streamer and land data with frequency-division encoding[J]. Geophysical Prospecting, 2012, 60(4): 663-680. DOI:10.1111/j.1365-2478.2012.01086.x
[115]
李闯, 黄建平, 李振春, 等. 平面波最小二乘逆时偏移编码策略分析[J]. 石油物探, 2015, 54(5): 592-601.
LI C, HUANG J P, LI Z C, et al. Analysis on the encoding strategies of plane-wave least-square reverse time migration[J]. Geophysical Prospecting for Petroleum, 2015, 54(5): 592-601. DOI:10.3969/j.issn.1000-1441.2015.05.012
[116]
田文彬, 张凯, 李振春. 平面波最小二乘逆时偏移方法的优化[J]. 石油物探, 2019, 58(2): 245-252.
TIAN W B, ZHANG K, LI Z C. Optimization of plane wave least squares reverse time migration[J]. Geophysical Prospecting for Petroleum, 2019, 58(2): 245-252. DOI:10.3969/j.issn.1000-1441.2019.02.010
[117]
HOU J, SYMES W W. Accelerating extended least-squares migration with weighted conjugate gradient iteration[J]. Geophysics, 2016, 81(4): S165-S179. DOI:10.1190/geo2015-0499.1
[118]
徐凯, 孙赞东. 基于粘声衰减补偿的最小二乘逆时偏移[J]. 石油物探, 2018, 57(3): 419-427.
XU K, SUN Z D. Least-squares reverse time migration based on visco-acoustic attenuation compensation[J]. Geophysical Prospecting for Petroleum, 2018, 57(3): 419-427. DOI:10.3969/j.issn.1000-1441.2018.03.011
[119]
ZHU F, HUANG J, YU H. Least-squares Fourier finite-difference pre-stack depth migration for VTI media[J]. Journal of Geophysics and Engineering, 2018, 15(2): 421-437. DOI:10.1088/1742-2140/aa9a0a
[120]
JIN S, MADARIAGA R, VIRIEUX J, et al. Two-dimensional asymptotic iterative elastic inversion[J]. Geophysical Journal International, 2010, 108(2): 575-588.
[121]
QU Y M, LI J L, HUANG J P, et al. Elastic least-squares reverse time migration with velocities and density perturbation[J]. Geophysical Journal International, 2018, 212(2): 1033-1056. DOI:10.1093/gji/ggx468
[122]
QU Y, GUAN Z, LI Z. Topographic elastic least:squares reverse time migration based on vector P‐and S‐wave equations in the curvilinear coordinates[J]. Geophysical Prospecting, 2019, 67(5): 1271-1295. DOI:10.1111/1365-2478.12775
[123]
ZHANG K, LI Z C, ZENG T S, et al. Residual curvature migration velocity analysis for angle domain common imaging gathers[J]. Applied Geophysics, 2010, 7(1): 49-56. DOI:10.1007/s11770-010-0006-1
[124]
秦宁, 李振春, 杨晓东. 基于角道集的井约束层析速度反演[J]. 石油地球物理勘探, 2011, 46(5): 725-731.
QIN N, LI Z C, YANG X D. Tomography velocity inversion by well constraint based on the angle domain common imaging gathers[J]. Oil Geophysical Prospecting, 2011, 46(5): 725-731.
[125]
秦宁.地震走时层析与波形反演方法研究[D].青岛: 中国石油大学, 2013
QIN N.Research on seismic traveltime tomography and waveform inversion[D].Qingdao: China University of Petroleum, 2013 http://cdmd.cnki.com.cn/Article/CDMD-10425-1015024320.htm
[126]
BOONYASIRIWAT C, VALASEK P, ROUTH P, et al. An efficient multiscale for time-domain waveform tomography[J]. Geophysics, 2009, 74(6): WCC59-WCC68. DOI:10.1190/1.3151869
[127]
曲英铭, 李振春, 黄建平, 等. 基于多尺度双变网格的时间域全波形反演[J]. 石油物探, 2016, 55(2): 241-250.
QU Y M, LI Z C, HUANG J P, et al. Full waveform inversion based on multi-scale dual-variable grid in time domain[J]. Geophysical Prospecting for Petroleum, 2016, 55(2): 241-250. DOI:10.3969/j.issn.1000-1441.2016.02.010
[128]
QU Y M, LI Z C, HUANG J P, et al. Elastic full waveform inversion for surface topography[J]. Geophysics, 2017, 82(5): R269-R285. DOI:10.1190/geo2016-0349.1