2017, Vol. 60
随着高精度地震勘探技术的不断发展，常规地震采集正在向低频、大偏移距以及宽方位地震采集方向转变，最小二乘偏移(LSM)、全波形反演(FWI)也逐步走向生产实际，特别是海上地震勘探已有不少成功实例.然而，地震各向异性集中表现为速度各向异性，对于各向异性探区，如果仅仅在各向同性假设条件下进行LSM或FWI，很可能会由于模拟数据与观测数据不匹配造成成像位置不准确、偏移噪声严重、迭代不收敛甚至发散等问题，严重影响反演成像质量.各向异性介质理论能够更加真实且全面地描述地震波在地层中的传播规律(Thomsen，1986；Tsvankin，1996；李娜等，2014)，因此，探讨并发展一种适用于各向异性介质的正演模拟、偏移成像及反演方法具有重要的理论和现实意义.
针对各向异性介质拟声波正演模拟及偏移成像，国内外学者进行了大量研究.第一类是基于波动方程的偏移成像方法：推导一种精确描述qP波传播特征的拟声波方程是各向异性数据准确成像的关键.Alkhalifah(1998)首次提出了TI介质声学近似，即通过设置沿着对称轴方向的剪切波速度为零，将TI介质全弹性波方程简化为四阶qP-qSV波耦合波动方程，进而在频率-空间域进行求解(Alkhalifah，2000；吴国忱和梁锴，2006)；Zhou等(2006a)进一步引入辅助函数使VTI介质四阶波动方程简化为二阶，在此基础上，很多学者将其推广到TTI介质并用于逆时偏移成像中，该类方法计算量较小、求解相对简单、但存在伪横波干扰及严重的不稳定现象(Grechka et al.，2004；Zhou et al.，2006b；李博等，2012；王娟等，2012；张岩和吴国忱，2013)；为了克服波场不稳定问题，当前有四种解决方案：第一，Flecter等(2009)提出的重新引入非零横波速度来修正上述波动方程以保证成像稳定；第二，从虎克定律出发，推导一种新的稳定波动方程用于成像处理中(Duveneck et al.，2008；Duveneck and Bakker 2011；Zhang and Zhang，2009；Zhang et al.，2011)；第三，采取Yoon等(2010)提出的各向异性参数匹配技术，将TTI介质极化倾角急剧变化的区域设置为椭圆各向异性；第四，采用基于纯qP波动方程的正演及成像算子，该类方法计算稳定，无伪横波干扰和由低速横波造成的频散假象，但方程中存在复杂的拟微分算子，计算效率低(Du et al.，2005；Chu et al.，2013；Pestana et al.，2011；Zhan et al.，2012，2013；程玖兵等，2013).此外，Hestholm(2009)借鉴一阶波动方程的计算优势，构建了VTI介质一阶波动方程正演模拟算子；吴国忱(2006)对各向异性qP波单程双域传播算子进行了深入研究(吴国忱等，2008；梁锴等，2009).
第二类是基于射线理论的偏移成像方法：相对于波动方程偏移，射线类成像方法兼具灵活、高效等优点.Kirchhoff偏移计算效率最高，仍是当今工业界最受欢迎的成像方法(刘立彬，2012；张洪宙和王彦飞，2014)，但是无法解决多波至、焦散区以及阴影区等成像问题；相比Kirchoff偏移，束偏移方法不仅计算效率高，而且在保幅成像、高陡构造成像、角道集直接提取等方面都具有明显的优势，具有较高的成像精度，近年来，各向异性介质叠前高斯束偏移成像方法成为当前业界的研究热点(Alkhalifah，1995；Zhu et al.，2007；刘鹏等2013；段鹏飞等，2013；韩建光等，2015；段新意等，2015).但总的来说，从反演的角度建立目标泛函通过迭代寻优的方式来处理各向异性成像问题较为少见.
为了对地下界面信息进行更加精确地成像，最小二乘逆时偏移(LSRTM)是当前最有潜力的成像方式之一.LSRTM旨在利用合成数据与观测数据之间的匹配程度来矫正常规RTM成像结果振幅失真、采样不规则、地震数据缺失等问题，通过迭代使成像结果逐渐趋于真实模型(Dai and Schuster，2009；Dai et al.，2011；刘玉金等，2013).Tarantola(1984)最早建立了最小二乘偏移的理论基础，随后，Nemeth等(1999)将最小二乘思想引入到各向同性声介质反演成像中，相继产生了最小二乘Kirchoff偏移、最小二乘单程波偏移以及LSRTM，该类方法具有均衡成像振幅，提高空间分辨率，减少偏移假象等优点(Duquet et al.，2000；Aoki and Schuster.，2009；Symes，2008；Dai et al.，2012；Dong et al.，2012；王彦飞，2013；黄建平等，2014a)；此外，在我国西部碳酸盐岩裂缝型储层成像、复杂介质保幅成像、近地表高精度成像等方面，LSM也具有一定的优势(黄建平等，2013a，2013b，2014b)；近年来，国内外学者还发展了黏声介质LSRTM，该方法能在成像过程中自动校正衰减所致的振幅误差和相位畸变(Dutta and Schuster，2014；李振春等，2014).
虽然LSRTM能够显著地改善偏移成像质量，但其计算成本也成倍增加(LSRTM迭代一次的计算量相当于标准RTM的2倍).鉴于LSRTM能很好地适应多震源数据成像，不少专家提出将多震源编码策略与LSRTM结合，既可提高计算效率又能保证成像质量，根据编码策略不同，主要分为线性编码与非线性编码，其中，平面波LSRTM就是一种经济高效的线性编码成像方式.平面波LSRTM通过合成平面波记录压缩地震数据，进而减少成像算法中炮循环次数和I/O成本(Tang，2009；Zhang et al.，2005；Huang and Schuster，2012；Dai and Schuster，2013；黄建平等，2014c；黄建平等，2015).然而，传统的平面波LSRTM忽略了介质的速度各向异性，导致模拟数据与观测数据不匹配，进而影响成像质量及收敛速度.
针对各向同性正演记录与各向异性观测记录不匹配的问题，本文首先实现了VTI介质LSRTM成像方法(VTI-LSRTM)，用以校正介质各向异性引起的相位畸变及振幅误差；其次，为了降低计算成本提高计算效率，进一步引入平面波编码策略，发展了基于平面波加速的VTI介质LSRTM方法(VTI-P-LSRTM)；最后通过简单模型及复杂Marmousi模型成像对比测试验证了本文方法在处理各向异性相位误差、提高成像信噪比及分辨率、补偿中深层能量、均衡反射波振幅等方面的优势，有利于解决小尺度缝洞、复杂小断块以及非均质体边界难以分辨的问题.
2 方法原理 2.1 LSRTM基本理论LSRTM旨在利用合成数据(或反偏移记录)与观测数据之间的匹配程度来矫正常规RTM成像结果中振幅失真及采样不规则等问题，可以通过如下最优化反演问题来描述(Dai and Schuster，2009；刘玉金等，2013；黄建平等，2015)：
|
(1) |
其中，J是目标泛函；dcal与dobs分别表示预测数据和观测数据，预测数据可以通过矩阵向量表示为dcal=Lm，L表示线性正演算子，m是与速度扰动有关的反射系数场，也即偏移成像结果.
反射系数场m可以通过最速下降法(Steepest Descent Method，SDM)或共轭梯度法(Conjugate Gradient Method，CGM)等进行迭代求解(Nemeth et al.，1999；Dai et al.，2012；王彦飞，2013)，其中CGM方法是在SDM方法的基础上，通过修正梯度方向和步长来加快收敛速度的一种改进方法.预条件CGM算法具体实现步骤如下：
|
(2) |
公式(2)中，k表示LSRTM迭代次数；gk和β分别表示最速下降法梯度方向和步长，其中梯度方向是通过对炮数据残差(即：Lm(k)－dobs)进行RTM成像得到；z和α分别表示在最速下降方向和步长的基础上进行修正后的共轭梯度方向及对应的步长；T为共轭转置或者伴随；C为预条件算子，引入合理的正则化约束条件可使反演问题寻优更稳定同时收敛速度更快，本文采用震源波场照明作为预条件算子(Plessix and Mulder，2004).
从上述算法可以得出，LSRTM方法主要计算包括五个步骤：通过①线性正演算子计算预测波场，进而求取炮数据残差；②利用偏移算子及梯度公式计算梯度方向，注意成像条件决定于梯度公式；③通过预条件算子修正梯度方向和步长得到共轭梯度方向和新的步长；④计算模型更新量，更新成像结果；⑤将数据残差或者模型残差与预先设定的阈值进行对比，判断是否需要进行下一次迭代.
2.2 VTI-LSRTM理论由上述介绍可知，VTI-LSRTM成像的关键在于推导适用于该类介质的线性正演算子、偏移算子以及梯度更新公式.下面将分别进行推导：
2.2.1 VTI介质反偏移算子构建对于各向异性介质，Duveneck等(2008)提出的二维VTI介质常密度拟声波方程组被广泛应用于qP波正演模拟与逆时偏移成像中，其时间-空间域表达式如下：
|
(3) |
其中，vpz，ε，δ都是描述介质各向异性特征的Thomsen参数，vpz表示沿着对称轴方向的qP波相速度，ε是度量qP波各向异性强度的无量纲因子，δ是连接纵横向qP波相速度的过渡性参数；p，q代表波场分量，且只是纵横向应力分量的线性组合，并无实际物理意义；fx，fz分别是x，z方向的体力源；x是位置坐标，t是时间.
根据扰动理论，相速度及各向异性参数都可以表示为背景介质与扰动介质之和的形式(Dutta and Schuster，2014)，即：
|
(4) |
公式(4)中，vpz0，ε0，δ0表示背景参数，dvpz，dε，dδ表示参数的微小扰动.为了简化，本文借鉴Dutta和Schuster(2014)提出的多参数扰动简化思路，只考虑影响程度最大的速度扰动，做如下两点假设：①密度是常数；②各向异性参数不存在扰动，即，
|
(5) |
此外，由于地震波满足线性叠加原理，在介质参数存在扰动时，地震波总场可由入射(或背景)波场与扰动(或散射)波场的叠加来表征：
|
(6) |
其中，p0，q0是通过背景介质利用公式(3)求得的背景波场；dp、dq是在介质存在扰动时产生的扰动波场.
将公式(4—6)代入公式(3)，并应用VTI介质Born弱散射假设条件(Clayton and Stolt，1981；吴国忱，2006；Plessix and Cao，2011)：即在vpz0≥dvpz的情况下，p0≥dp，q0≥dq.此外定义反射系数模型，即
|
(7) |
公式(7)即为VTI介质反偏移公式，该公式建立了反射系数模型与扰动波场之间的线性关系.从公式(7)可知，欲求扰动场，需进行两步正演模拟：第1步是根据公式(3)，以背景参数、震源函数作为输入求解入射波场p0和q0；第2步是通过公式(7)，以背景参数、入射波场以及每次迭代的近似成像结果作为输入求取扰动场dp和dq.与各向同性介质反偏移公式不同之处在于，①构建正演算子需同时输入背景速度和各向异性参数，而不仅仅是纵波速度；②输出的地震波场也不再是单一的声压波场，而是一对地震波场分量.
仿照各向同性LSRTM，公式(7)还可以通过矩阵形式表示为
|
(8) |
这里，L即为线性正演算子，反映模型参数扰动对数据扰动的影响；m是反射系数模型，代表LSRTM每次迭代的成像结果；d表示考虑各向异性的Born正演数据，通过公式(3)和(7)分两步进行计算.
2.2.2 VTI介质偏移算子及梯度更新公式在数学上，偏移算子是线性正演算子的伴随(Plessix，2006)，偏移算子可以通过伴随状态法求解.首先定义(p，q)的伴随变量为(r，w)，则伴随变量满足如下公式：
|
(9) |
公式(9)中，Δprecv和Δqrecv表示模拟数据与观测数据的残差；xs为震源坐标，xg为检波点坐标.
其次，VTI介质误差泛函可表示为p、q两个分量的残差在L2模的条件下同时达到最小，如下所示：
|
(10) |
其中，U、dobs分别表示模拟数据向量和观测数据向量，如公式(11)所示；对于VTI介质，模型参数m可以是速度或者各向异性参数等.
|
(11) |
LSRTM的另一个关键步骤就是推导泛函梯度.Tarantolla(1984)提出可以通过剩余波场反传结果与正传波场的互相关来计算负梯度方向，进而实现波形迭代反演.本文推导了VTI介质误差泛函的梯度更新公式：
|
(12) |
进一步结合照明补偿预条件算子(Plessix and Mulder，2004)，得到震源归一化梯度更新公式：
|
(13) |
在公式(13)中，通过公式(3)计算正传波场(p，q)，基于公式(9)以炮数据残差作为输入计算伴随波场(r，w)，最终求得梯度方向g.用矩阵表示为
|
(14) |
式中，LT即为偏移算子，反映模型参数扰动对成像扰动的影响；Δd表示数据向量残差.因为偏移算子只是正演算子的共轭，而不是其逆，所以第一次偏移结果仅仅是真实模型的一个近似.从公式(14)可以得出：泛函梯度是通过偏移算子来构建的，将偏移算子作用于数据向量残差即可求取泛函梯度，其成像条件决定于梯度公式.
至此，正演算子、偏移算子及梯度更新公式推导完毕.基于上述推导，可以得到VTI-LSRTM方法处理流程.
从流程图 1可知，VTI-LSRTM是将VTI-RTM的成像剖面作为初始模型，随后，通过迭代不断更新使成像结果逼近真实反射系数.与各向同性介质LSRTM方法不同，在VTI-LSRTM中，正演算子L和偏移算子LT的构建都需要同时输入背景参数vpz，ε，δ，并根据相应的波动方程数值求解算法来实现；在FWI中，背景参数vpz，ε，δ随着迭代次数不断更新，而在LSRTM中保持不变；此外，为了保证最小二乘迭代收敛，偏移算子与反偏移算子必须互为共轭，其共轭特性可以通过点乘测试来验证(Claerbout，1992).
|
图 1 VTI-LSRTM算法流程 Fig. 1 Flow chart of LSRTM for VTI media |
对于常规逐炮偏移，VTI-LSRTM方法占用内存大、计算成本高，而平面波偏移是一种采取多炮同时偏移来提高计算效率的加速手段(Dai et al.，2012；黄建平等，2015).其第1步是地震数据压缩，即将炮域的地震数据按照线性时移编码(也称平面波编码，plane-wave encoding)合成平面波记录，在x-t域，可用射线参数p来表示经过不同时移的平面波记录：
|
(15) |
其中，θ是地表激发倾角，v是在地表位置θ方向的纵波相速度.对于二维观测方式，平面波记录可表示为
|
(16) |
在公式(16)中，d(xg，t;xs)表示炮域的地震数据；d(xg，t;p)表示平面波域地震数据，即平面波记录；δ(t－p·｜xs｜)是时移函数，*表示卷积运算，xs是每一炮震源坐标，xg是接收点坐标，时移量为p·｜xs｜，它随震源坐标线性变化.因此，平面波编码过程可以表述为：首先将炮域的地震数据进行逐炮时移，然后进行线性叠加，最终将多炮地震数据压缩成1个平面波记录.
第2步即为平面波偏移，使用相同编码方式的平面波震源激发产生震源波场，将平面波记录作为边界进行反向传播获得检波点波场，互相关成像即获得平面波成像结果，可见其计算量仅相当于1炮地震数据的偏移成本.
然而单个平面波成像结果通常会遭受严重的串扰噪声，为了在提高计算效率的同时保证计算精度，通常将多个不同编码方式的平面波成像结果进行叠加，以此来增强有效信号，减弱串扰噪声.即平面波域最终成像结果可以表式为
|
(17) |
在公式(17)中，pi是合成第i个平面波记录所选用的射线参数；mi是第i个平面波记录作为输入获得的成像结果，m是最终成像结果；N是总的平面波记录个数，通常远小于炮数.
与处理流程1类似，通过输入编码后的平面波记录，采用上述平面波偏移策略进行成像处理和梯度计算，即可完成基于平面波加速的VTI介质逆时偏移(VTI-P-LSRTM).需要说明的是，在构建平面波反偏移算子时，其震源也需要按照合成平面波时的编码方式进行编码.
3 模型测试本文分别选用了两组模型来检验VTI-LSRTM及VTI-P-LSRTM方法的正确性和适应性.一组是包含强各向异性地质体的简单层状模型，另一组是根据Alkhalifah(1997)的观点进行改造的各向异性复杂Marmousi模型.本文采用时间二阶-空间十阶Lebedev网格高阶有限差分正演模拟方法合成观测数据，通过Cerjan衰减边界条件处理边界反射.在成像处理中，使用震源照明补偿作为预条件算子以保证收敛速度.
3.1 简单模型试算 3.1.1 模型参数与观测系统首先，通过对一个简单层状模型进行成像试算，观察分析地震各向异性对偏移成像的影响，并初步测试本文方法的正确性.图 2(a，b)分别是真实速度场和偏移速度场，利用真实速度和偏移速度根据反射系数模型的定义可求取相对慢度扰动(Relative slowness perturbation).图 2(c，d)分别为经过平滑处理后的Epsilon和Delta模型，对模型参数进行适当地平滑处理可以使得地震波正、反向延拓过程更加稳定.模型中各向异性异常体的Thomsen参数大小分别为ε=0.2，δ=-0.1，这表明①地震波传播至该异常体中时会表现出很强的非椭圆各向异性特征；且②地震波横向传播速度大约为垂向传播速度的1.18倍.本文设置了强各向异性倾斜地层，目的在于突出基于各向同性假设的偏移成像算法在TI介质成像中的不足，同时展示本文方法在成像处理中的优越性.模型网格的大小为451×301，纵、横向网格间距为10×10 m，因此，横向范围总计4500 m，最大深度3000 m.浅层近地表速度为2500 m·s-1，正演子波选取主频为30 Hz的雷克子波，经计算，每个空间波长内至少能采8个点，既避免了数值频散，又兼顾了计算效率和精度.
|
图 2 简单模型参数 (a)真实纵波速度;(b)偏移纵波速度;(c)Epsilon模型;(d)Delta模型. Fig. 2 A simple model (a)True P-velocity model;(b)Migration P-Velocity model;(c)Smoothed epsilon; (d)Delta model used for VTI-RTM and VTI-LSRTM. |
为了形成照明充分的观测数据，本文设计如下观测系统：设置226炮激发，451道地表满接收的观测方式，炮点井深10 m，炮间距20 m，道间距10 m，第一炮炮点位于CDP号等于1的位置.计算参数如下：总记录时间为2.0 s，时间采样步长为0.5 ms，每道4001个采样点，代入改造后的CFL(Courant-Friedrichs-Lewy)条件可知：上述参数能够保证计算过程稳定(Courant et al.，1928).利用图 2中的真实速度与各向异性参数建立VTI介质正演算子并合成相应的观测数据，再以偏移速度及正确的各向异性参数构建偏移算子，依次选取各向同性RTM/LSRTM、各向异性RTM/LSRTM以及平面波成像方法对观测数据进行成像处理，处理结果将分别在3.1.2和3.1.3节进行详细对比分析.
3.1.2 常规LSRTM成像剖面图 3(a—d)分别展示了传统的各向同性声波RTM、声波LSRTM、VTI-RTM、以及VTI-LSRTM迭代30次的成像结果.对比图 3(a，b)声波RTM及LSRTM成像结果可知：①对于浅层同相轴，与标准RTM相比，LSRTM成像结果有相对更高的成像分辨率和较少的成像噪声；然而，②在各向异性地层中，无论是RTM还是LSRTM成像剖面，都存在较为严重的伪横波成像干扰(如图中实线箭头所示)，且倾斜地层出现明显的分叉现象，多出一条上拉的虚假同相轴(如图中椭圆线圈所示)，这是由于速度各向异性(其他方向的波速大于垂向速度)造成观测记录中该反射层的反射波走时减小所致；此外，③各向异性异常体下方的深部同相轴出现了明显的相位畸变、轻微的振幅异常、绕射波不收敛等，且同相轴连续性较差(如图中虚线箭头所示).值得一提的是，随着声波LSRTM迭代次数的增加，上述三种现象不但没有得到改善，反而更加严重，这表明介质各向异性严重影响了经典的地震波传播规律，特别是地震波相位信息，因此，在迭代过程中合成数据与观测数据出现不匹配现象，导致成像效果不佳.
|
图 3 常规成像剖面 (a)各向同性声波RTM;(b)声波LSRTM迭代30次的成像结果;(c)VTI-RTM;(d)VTI-LSRTM迭代30次的成像结果. Fig. 3 The images with(a)the conventional acoustic RTM；(b)acoustic LSRTM with 30 iterations；(c)the VTI-RTM；(d)VTI-LSRTM with 30 iterations. |
采用VTI-RTM及本文提出的VTI-LSRTM成像方法对同样的观测数据进行成像试算，成像结果分别如图 3(c，d)所示.①首先对比各向同性与VTI介质成像结果可知，通过考虑介质各向异性特征，VTI-RTM与VTI-LSRTM方法均消除了由于各向异性异常体造成的倾斜地层成像分叉、伪横波成像干扰、下覆地层相位畸变等现象，绕射波准确归位，成像质量得到显著改善.②其次，从图 3(b，d)可见，与标准的各向同性声波LSRTM相比，VTI-LSRTM成像方法既有相似之处，又有其独特优势，相似之处体现在浅层各向同性地层成像中，两种方法均能获得较高分辨率且振幅相对保持的成像结果，其优势在于当处理各向异性体以及下覆地层时，随着迭代次数的增加，反射同相轴能够正确归位，成像振幅变得更加均衡，不会因数据不匹配导致成像紊乱或不收敛现象.③最后，通过比较VTI介质RTM与LSRTM成像结果(图 3(c，d))可以清楚地看到，VTI-LSRTM在偏移噪声压制、深部构造刻画、震源效应抑制等方面都明显优于VTI-RTM的成像结果，整体上成像剖面也变得更加清晰，在下文的复杂模型测试中将进一步验证这一观点.需要指出的是，VTI-RTM与VTI-LSRTM与常规各向同性成像方法类似，仍然存在一些当前难以克服的成像噪声，如图 3(c，d)虚线矩形框所示，VTI-RTM的计算时间大约是标准RTM的5倍.
在最小二乘反演方法研究中，收敛曲线经常作为判定最小二乘成像方法优劣的标准.为了更加直观地展示VTI-LSRTM在成像中的优越性，图 4分别给出了各向同性和VTI-LSRTM方法的数据残差随迭代次数变化的关系曲线.图中带圈的实线与带点的实线分别表示各向同性LSRTM与VTI-LSRTM的收敛曲线，通过仔细观察对比可以得出如下结论： ① 在处理遭受各向异性影响的地震数据时，VTI-LSRTM的收敛速度比声波LSRTM更快，这是因为考虑各向异性特征的正演算子与偏移算子互为共轭，合成数据与观测数据更加匹配； ② 两种方法的反演过程都较为稳定，总的趋势是随着迭代次数的增加，数据残差逐渐减小，合成数据逐步逼近于观测数据，成像结果趋于真实模型； ③ 随着迭代过程的进行，残差减小的速度逐渐变慢，因此，在实际生产当中，选择合适的迭代次数需要兼顾成像精度与计算成本.
|
图 4 数据残差归一化收敛曲线 Fig. 4 The normalized convergence curves of data residuals |
第1步需要按照公式(16)合成如下三组平面波记录，第I组，仅含1个p=0.0的平面波记录；第Ⅲ组，具有15个不同参数p的平面波记录，p的取值范围是-0.2~0.2 s·km-1(对应的倾角范围是-30°~30°)，并且按照线性关系变化；第Ⅲ组，30个平面波记录，p的取值范围及变化趋势与第Ⅲ组相同.图 5(a—c)分别展示了p=-0.2、0.0及0.2(对应的地表激发倾角是-30°、0°和30°)三种不同射线参数的平面波记录，图中可见： ① 通过平面波编码，将多炮数据压缩成几个类似单炮的数据体; ② 记录中主要的反射波形态由“双曲型”变为“直线型”，且斜率与射线参数及反射界面形态有关，这正是因为采取了线性时移编码方式，为此，从平面波记录中可以粗略地识别浅层构造形态(如图 5b); ③ 然而复杂构造及其下覆地层所对应的反射波同相轴仍然十分杂乱，这是由于球面波记录相互干涉叠加的结果.然后分别采用基于平面波加速的各向同性LSRTM(P-LSRTM)与VTI介质LSRTM(VTI-P-LSRTM)对上述三组平面波地震记录进行成像测试.
|
图 5 平面波记录射线参数分别是(a)p=-0.2 s·km-1；(b)p=0.0；(c)p=0.2 s·km-1. Fig. 5 Plane-wave gathers with(a)p=-0.2 s·km-1;(b)p=0.0;(c)p=0.2 s·km-1. |
成像结果分别如图 6(a—f)所示，第一行分别对应以上三组记录各向同性介质平面波LSRTM迭代30次的成像剖面，而图 6底部的三个剖面分别是相应的VTI介质平面波LSRTM成像结果.通过仔细观察对比可知：①通过增加平面波个数，无论是各向同性还是VTI介质平面波偏移结果，其串扰噪声都基本被压制，剖面的信噪比显著提高，这说明将多个不同编码方式的平面波成像结果进行叠加，有利于增强有效信号，压制成像噪声；②从算法方面来看，如果考虑介质的各向异性，其偏移剖面的成像质量明显优于各向同性介质平面波LSRTM成像结果，特别是当平面波个数大于1，串扰噪声得到较好压制时其算法对成像效果的改善尤为突出，这是因为VTI介质平面波LSRTM成像算法有效地校正了由于各向异性造成的相位畸变、同相轴分叉、振幅异常以及伪横波干扰等现象，提高了剖面信噪比及成像分辨率；③此外，采取平面波加速策略之后的VTI介质LSRTM成像效果与常规LSRTM方法相近，如图 6(d，f)所示，其成像质量基本接近于图 3d，但计算成本将成倍减少，以15个平面波偏移结果为例，其计算效率相对于常规LSRTM可提高约15倍，这表明基于平面波加速的VTI介质LSRTM成像方法能够在保证成像质量的同时，显著提高计算效率，减少计算成本，是一种极具潜力的成像策略.
|
图 6 平面波成像结果 (a)与(b)分别是1个平面波记录；(c)与(d)分别是15个平面波；(e)与(f)分别是30个平面波的ISO-LSRTM与VTI-LSRTM迭代30次的成像结果. Fig. 6 The images obtained with(a)plane-wave ISO-LSRTM and(b)plane-wave VTI-LSRTM of one plane-wave records and 30 iterations，(c)plane-wave ISO-LSRTM and(d)plane-wave VTI-LSRTM of fifteen plane-wave records and 30 iterations，(e)plane-wave ISO-LSRTM and(f)plane-wave VTI-LSRTM of thirty plane-wave records and 30 iterations. |
综上所述：对于遭受各向异性地质体影响的地震数据，如果采取本文提出的VTI介质平面波LSRTM成像算法处理相位畸变等问题，同时通过增加平面波个数压制串扰噪声，且不断迭代更新成像结果，那么三者结合即可有效地提高剖面信噪比和成像分辨率，改善成像质量，获得清晰可靠的地震剖面.
为了更全面地展示采用平面波加速策略的优势，本文给出了上述第Ⅲ组平面波LSRTM的数据残差及模型残差收敛曲线，分别如图 7(a，b)所示.从图 7可见：随着迭代次数的增加，无论是各向同性还是VTI介质平面波LSRTM方法，其反演过程都较为稳定，数据残差和模型残差都不断减小，成像结果逐渐逼近真实反射率模型；但是基于平面波加速的VTI-LSRTM方法具有更好的收敛效果，主要表现在①当迭代次数一定时，其成像结果更加接近于真实模型，如图 7中在第30次迭代结束时，成像结果已基本收敛，数据残差趋于0，模型残差降到40%以下；②随着迭代次数增加，其残差减小的速度明显大于各向同性的情况，如图 7所示两条曲线之间拉开的距离越来越大，模型残差更为明显；③当迭代过程进行到一定程度时，各向同性介质成像方法的残差几乎不再下降，而VTI介质成像方法将仍然保持平稳收敛.
|
图 7 归一化收敛曲线(30个平面波记录) (a)数据残差；(b)模型残差. Fig. 7 The normalized convergence curves of different LSRTM method with thirty plane-wave records (a)The data residuals；(b)The model residuals. |
通过简单模型成像对比试验，基本验证了本文VTI-LSRTM成像方法的正确性及平面波加速策略在VTI介质反演成像中的效率优势.
3.2 复杂Marmousi模型验证 3.2.1 Marmousi速度模型及各向异性参数进一步，将本方法应用于复杂Marmousi模型成像试处理，以检验方法对复杂模型的适应性.图 8(a，b)分别展示了经过抽稀处理后的纵波真实速度场及用于成像处理的偏移速度场，另外，作者借鉴Alkhalifah(1997)提出的观点，建立了各向异性Thomsen参数模型，分别如图 8(c，d)所示.该模型反射结构复杂、断裂及高陡构造发育、纵横向速度变化剧烈，是检验LSRTM成像方法保幅性及成像适应性的经典模型.模型网格大小为369×188，计算过程中网格间距为12.5 m，因此，模型水平方向尺度为4600.0 m，垂直方向尺度为2337.5 m，近地表速度为1500 m·s-1.
|
图 8 改造的Marmousi模型 (a)真实速度；(b)偏移速度；(c)Epsilon模型；(d)Delta模型. Fig. 8 The modified Marmousi models (a)True P-velocity model；(b)Migration P-Velocity model；(c)Smoothed epsilon；(d)Delta model used for VTI-RTM and VTI-LSRTM. |
在正演及成像过程中，观测系统如下：初始炮点位置在第一个网格点处，炮点深度12.5 m，炮点间隔25.0 m，共185炮激发，在地表每个网点都接收，总计369个接收道.计算参数设计为：震源选取主频为25 Hz的雷克子波，时间采样间隔0.5 ms，每道共记录7001个采样点，总记录时间为3.5 s，经计算，符合CFL稳定性条件.使用上述参数计算得到地震炮记录之后，根据公式(16)合成15个不同编码方式的平面波记录，其p的取值范围及变化趋势与上述3.1.3 节第Ⅲ组相同，基于公式(15)，经计算可知对应的倾角范围是-17.5°~17.5°.
图 9(a—c)分别展示了射线参数p=-0.2、0.0和0.2三种情况下的平面波记录，观察可见：道集中含有较为丰富的反射波信息，但由于模型结构呈现出严重的非均质性，导致平面波响应特征十分复杂，难以直接从平面波记录中有效地辨别出地下构造形态，需要进一步执行成像处理，从而使反射波归位和绕射波收敛.在成像处理之前，首先切除了能量较强的直达波信号，目的在于压制偏移噪声，提高浅层成像质量，最后将上述编码前后的地震记录用于成像测试中.
|
图 9 平面波记录射线参数依次为(a)p=-0.2 s·km-1；(b)p=0.0；(c)p=0.2 s·km-1. Fig. 9 Plane-wave gathers with(a)p=-0.2 s·km-1；(b)p=0.0；(c)p=0.2 s·km-1 |
各种成像方法对应的成像结果分别如图 10(b—f)所示，图 10a是相对慢度扰动，需要说明的是，它并不是真实的反射系数模型，只是速度扰动与背景速度的综合反映，因为VTI介质的反射系数是由介质界面两侧的垂向纵波速度、各向异性参数以及入射角度共同决定的，所以成像结果仅仅是反射系数的近似.
|
图 10 慢度扰动及成像剖面对比 (a)相对慢度扰动模型；(b)各向同性LSRTM迭代30次的成像结果；(c)和(d)分别是VTI-LSRTM迭代1次及30次的成像剖面； (e)和(f)分别是基于平面波加速的VTI-LSRTM迭代1次及30次的结果. Fig. 10 Relative slowness perturbation and the comparison among images from(b)acoustic LSRTM with 30 iterations，(c)VTI-RTM，(d)VTI-LSRTM with 30 iterations，(e)plane-wave VTI-RTM with thirty plane-wave records，(f)plane-wave VTI-LSRTM of thirty plane-wave records and 30 iterations |
图 10b是各向同性介质LSRTM成像结果，从图中可见，由于没有考虑各向异性效应的影响，致使同相轴发生错乱，断层边界及深部构造难以分辨，整个剖面被杂乱的噪声所掩盖；图 10(c，d)分别是VTI介质RTM及LSRTM成像结果，对比可知，两剖面的总体成像效果都明显优于图 10b，反射同相轴正确归位，成像位置较为清晰，由相位畸变等造成的同相轴紊乱现象基本被消除，但是VTI-RTM成像剖面中浅层同性轴受震源效应的影响较为严重、中深部同性轴能量较弱、采用Laplace滤波压制低频噪声导致振幅保持性较差，然而VTI-LSRTM通过迭代自动压制低频噪声，剖面中震源效应被抑制，断层、背斜、高陡构造以及不整合面等地下复杂构造准确成像，中深层反射能量也得到补偿，其成像质量最高；图 10(e，f)分别是基于平面波加速的VTI-RTM及VTI-LSRTM成像结果，其成像质量分别接近于VTI-RTM及VTI-LSRTM.
需要强调的是，VTI-LSRTM是一种逐炮偏移方法，计算效率低，而平面波偏移是一种多炮同时进行偏移的加速策略，本例中VTI介质平面波LSRTM的计算效率约为VTI-LSRTM的6.2倍，此外，通过平面波编码将地震记录压缩成平面波记录，有效地降低了I/O以及内存需求，因此，本文提出的VTI介质平面波LSRTM是一种兼具效果与效率的双赢成像策略.
3.3 速度模型敏感性测试最后，为了探索本方法对背景速度模型的精度要求，以进一步指导其在实际生产中的应用，本文还进行了方法的速度误差敏感性分析测试.
在测试过程中，假定各向异性参数是正确的，成像误差仅仅来源于偏移速度误差.真实模型如图 8(a，c，d)所示，并通过如下两种方式来引入偏移速度误差，第I种是对速度模型按照一定的比例同时增大或减小，模拟系统误差；第Ⅲ种是对速度模型进行平滑处理，模拟随机误差，选择两种方式来引入偏移速度误差的原因在于不同类型的误差对成像效果的影响程度不尽相同(曹文俊等，2007；Luo and Hale，2014).据此，本文按照相对误差大小依次设计了4组带有误差的偏移速度模型，分别如图 11A一列所示，依次命名为第1、2、3、4组.第1组与第2组系统误差相同，第2组与第3组随机误差一致，第4组设计较大的误差值，其相对误差大小分别为：3.14%、5.79%、7.49%、9.58%，从图中可见： ① 随着随机误差的增大，偏移速度场逐渐变得模糊，断层等构造形态变得难以分辨； ② 系统误差虽不改变速度结构的形态，但从颜色条可知，速度值被同时减小.分别将上述4组偏移速度场作为输入进行偏移成像测试，观测系统、计算参数以及平面波记录(平面波个数等于15)与3.2.1 节所述一致.
|
图 11 偏移速度模型及对应成像剖面 (A)含不同误差的偏移速度模型；(B)平面波VTI-LSRTM迭代30次的成像结果. Fig. 11 Sensitivity of plane-wave VTI-LSRTM to errors in the migration velocity model (A)The figures in the left column show the different velocity models used for plane-wave VTI-LSRTM； (B)The figures in the right column are the corresponding plane-wave VTI-LSRTM images after 30 iterations. |
VTI介质平面波LSRTM迭代30次的成像结果分别如图 11B一列所示，由图中仔细观察： ① 对比第1、2组成像结果可知，当系统误差一定时，偏移速度场中的随机误差几乎不会引入成像误差，成像剖面比较清晰； ②对比第2、3组成像剖面可得，当随机误差一定时，随着系统误差的增大，反射层变得复杂、错乱，特别是深部同相轴难以辨认，构造也未能准确成像，由于速度模型总体降低，因此成像剖面有整体抬升的趋势； ③ 随着相对误差的逐渐增大，成像效果随之变差，但根本影响因素在于系统误差增大.此外，本方法对各向异性参数误差的敏感性测试将另文发表.
为了更清楚地展示速度误差对反演迭代过程的影响，本文给出了上述4组测试的数据残差收敛曲线，如图 12所示，对比可知： ① 相对于采用真实速度模型的成像结果，由于速度误差的存在，收敛效果总体较慢； ② 随着相对误差的增大，特别是系统误差的增加，收敛效果逐渐变差； ③ 当误差增加到8.0%左右时，随着迭次次数增加，误差基本不收敛.经推测，本方法可容许的最大速度相对误差大小大约为8.0%.
|
图 12 数据残差归一化收敛曲线 Fig. 12 The normalized convergence curves of data residuals |
本文首先将最小二乘反演思想引入到各向异性介质的成像处理中，实现了VTI介质LSRTM方法(VTI-LSRTM)；进一步将其与平面波编码技术结合，发展了基于平面波加速的VTI介质LSRTM方法.通过简单模型以及改造后的复杂Marmousi模型成像试验可知，相对于各向同性声波LSRTM，该方法除了具备标准LSRTM的优点以外，还有如下几点优势：
① VTI-LSRTM方法弥补了各向同性正演算法的不足，消除了预测数据与观测数据之间由于模拟方法带来的成像误差，特别是走时(或者相位)误差，有利于解决各向异性探区成像效果不佳的问题；
② 基于平面波加速的VTI介质LSRTM通过将炮域的地震数据合成平面波记录，极大地减少了内存及I/O成本，是一种兼具成像精度与计算效率的成像方式；
③ 随着迭代次数增加，数据误差逐步降低，VTI-LSRTM成像结果始终保持稳定，深部构造成像质量明显改善.
因此，与标准LSRTM相比，VTI-LSRTM是一种较为有效的提高VTI介质成像质量的方法.尽管该类方法对VTI介质成像有一定的优势，但是也应该看到该类方法存在的不足：
① 受构造运动的影响，地球介质还表现为极化各向异性、方位各向异性以及正交各向异性特征，VTI-LSRTM对此类介质的精确成像尚有不足；
② 该类方法需要同时具备相对正确的偏移速度场和各向异性Thomsen参数.目前，尽管已经发展多种各向异性参数提取方法(刘玉柱等，2015)，但对各向异性Thomsen参数的估计通常是不正确的，其反演难度比速度估计大得多，这势必制约本方法的实用性；
③ 本文分别采用了震源波场照明补偿和线性编码技术来提高计算精度和效率，还可以进一步考虑其他预条件算子、非线性编码等技术来提高本方法的适应性；
④ 地震子波、背景参数的精度、对地球介质的理论假设以及地震数据品质等等都是影响最小二乘偏移效果的关键因素，要真正实现模拟记录和观测记录的最佳匹配，还必须考虑地震子波，数据品质等因素的影响.
综上所述，在今后的研究过程中，可以重点研究各向异性参数反演策略来提高参数估计的准确性，同时修正目标泛函降低方法对速度误差的敏感性，针对地震子波等因素的影响，可以借鉴各向同性介质中互相关目标泛函来解决子波估计困难的问题；其次，还应引入非线性编码、寻求可行性约束、考虑GPU加速等，尽可能在保证LSRTM成像精度的同时进一步提高其计算效率；此外，发展TTI介质纯纵波LSRTM方法也具有极其重要的意义.
附录A在该附录中，给出了本文VTI介质反偏移公式(7)的详细推导过程.首先，将公式(4)和公式(6)代入公式(3)可得：
|
(A1) |
利用假设条件(即简化公式(5))消除影响较小的各向异性参数扰动，此外，考虑到背景参数及由此产生的入射波场(或称为背景波场)应当满足如下波动方程：
|
(A2) |
为此，将式(A2)代入式(A1)并经过简单的移项处理，可将式(A1)重新写为
|
(A3) |
观察可知，该方程组左右两端都存在未知波场，扰动波场与速度扰动之间表现为非线性关系，如果进行直接求解，其计算成本极其昂贵(Clayton and Stolt，1981；吴国忱，2006；Plessix and Cao，2011).为了推导两者之间的线性关系，本文在公式(A3)的基础上，通过应用Born弱散射假设条件，即当vpz0≥dvpz时，p0≥dp、q0≥dq，对此，分别采用入射波场p0及q0代替方程组(A3)右端项中的总场p0+dp及q0+dq，并再次结合入射波场所满足的方程(A2)，可将式(A3)进一步简化为
|
(A4) |
忽略高阶速度扰动项(dvpz0)2，并代入反射系数模型m的定义式，最终可以得到：
|
(A5) |
公式(A5)建立了反射系数模型与扰动波场之间的线性关系.分析公式(A5)可知：扰动波场与入射波场具有相同的格林函数，都是通过背景参数借助有限差分等方式来构建的，而不同之处在于震源项，以入射波场的时间二阶导数与反射系数模型的乘积作为震源，即可激发产生扰动波场.
致谢作者感谢两位匿名审稿专家的宝贵意见.感谢国家重点基础研究发展计划(2014CB239006)，国家自然基金(41274124)，中央高校科研业务费专项基金(R1401005A)项目的联合资助.
| Alkhalifah T. 1995. Gaussian beam depth migration for anisotropic media. Geophysics, 60(5): 1474-1484. DOI:10.1190/1.1443881 | |
| Alkhalifah T. 1997. An anisotropic Marmousi model. Stanford Exploration Project, SEP-95:265-282. | |
| Alkhalifah T. 1998. Acoustic approximations for processing in transversely isotropic media. Geophysics, 63(2): 623-631. DOI:10.1190/1.1444361 | |
| Alkhalifah T. 2000. An acoustic wave equation for anisotropic media. Geophysics, 65(4): 1239-1250. DOI:10.1190/1.1444815 | |
| Aoki N, Schuster G T. 2009. Fast least-squares migration with a deblurring filter. Geophysics, 74(6): WCA83-WCA93. DOI:10.1190/1.3155162 | |
| Cao W J, Li Z C, Hang W G. 2007. Sensitivity analysis of wave equation pre-stack depth migration to errors in migration velocity model. Natural Gas Industry (in Chinese), 27(S1): 120-122. | |
| Cheng J B, Kang W, Wang T F. 2013. Description of qP-wave propagation in anisotropic media, Part I:Pseudo-pure-mode wave equations. Chinese J. Geophys. (in Chinese), 56(10): 3474-3486. DOI:10.6038/cjg20131022 | |
| Chu C L, Macy B K, Anno P D. 2013. Pure acoustic wave propagation in transversely isotropic media by the pseudospectral method. Geophysical Prospecting, 61(3): 556-567. DOI:10.1111/gpr.2013.61.issue-3 | |
| Claerbout J F. 1992. Earth Soundings Analysis:Processing Versus Inversion. London:Blackwell Scientific Publications. | |
| Clayton R W, Stolt R H. 1981. A Born-WKBJ inversion method for acoustic reflection data. Geophysics, 46(11): 1559-1567. DOI:10.1190/1.1441162 | |
| Courant R, Friedrichs K, Lewy H. 1928. über die partiellen differenzengleichungen der mathematischen physik. Mathematische Annalen, 100(1): 32-74. DOI:10.1007/BF01448839 | |
| Dai W, Schuster J. 2009. Least-squares migration of simultaneous sources data with a deblurring filter.//79th Annual International Meeting, SEG, Expanded Abstracts, 2990-2994. | |
| Dai W, Wang X, Schuster G T. 2011. Least-squares migration of multisource data with a deblurring filter. Geophysics, 76(5): R135-R146. DOI:10.1190/geo2010-0159.1 | |
| Dai W, Fowler P, Schuster G T. 2012. Multi-source least-squares reverse time migration. Geophysical Prospecting, 60(4): 681-695. DOI:10.1111/gpr.2012.60.issue-4 | |
| Dai W, Schuster G T. 2013. Plane-wave least-squares reverse-time migration. Geophysics, 78(4): S165-S177. DOI:10.1190/geo2012-0377.1 | |
| Dong S, Cai J, Guo M, et al. 2012. Least-squares reverse time migration:Towards true amplitude imaging and improving the resolution.//82nd Annual International Meeting, SEG, Expanded Abstracts, 10:1190-1194. | |
| Du X, Bancroft J C, Lines L R. 2005. Reverse-time migration for tilted TI media.//75th Annual International Meeting, SEG, Expanded Abstracts, 24:1930-1934. | |
| Duan P F, Cheng J B, Chen A P, et al. 2013. Local angle-domain Gaussian beam prestack depth migration in a TI medium. Chinese J. Geophys. (in Chinese), 56(12): 4206-4214. DOI:10.6038/cjg20131223 | |
| Duan X Y, Li Z C, Huang J P, et al. 2014. A prestack Gaussian beam depth migration in common-shot domain for anisotropic media. Geophysical Prospecting for Petroleum (in Chinese), 53(5): 579-586. | |
| Duquet B, Marfurt K J, Dellinger J A. 2000. Kirchhoff modeling, inversion for reflectivity, and subsurface illumination. Geophysics, 65(4): 1195-1209. DOI:10.1190/1.1444812 | |
| Dutta G, Schuster G T. 2014. Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation. Geophysics, 79(6): S251-S262. DOI:10.1190/geo2013-0414.1 | |
| Duveneck E, Milcik P, Bakker P M, et al. 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration.//78 Annual International Meeting, SEG, Technical Program Expanded Abstracts, 2186-2190. | |
| Duveneck E, Bakker P M. 2011. Stable P-wave modeling for reverse-time migration in tilted TI media. Geophysics, 76(2): S65-S75. DOI:10.1190/1.3533964 | |
| Fletcher R P, Du X, Fowler P J. 2009. Reverse time migration in tilted transversely isotropic (TTI) media. Geophysics, 74(6): WCA179-WCA187. DOI:10.1190/1.3269902 | |
| Grechka V, Zhang L B, Rector J W. 2004. Shear waves in acoustic anisotropic media. Geophysics, 69(2): 576-582. DOI:10.1190/1.1707077 | |
| Han J G, Wang Y, Zhang X B. 2015. Gaussian beam prestack depth migration in VTI media. Oil Geophysical Prospecting (in Chinese), 50(2): 267-273. | |
| Hestholm S. 2009. Acoustic VTI modeling using high-order finite differences. Geophysics, 74(5): T67-T73. DOI:10.1190/1.3157242 | |
| Huang J P, Li Z C, Kong X, et al. 2013a. A study of least-squares migration imaging method for fractured-type carbonate reservoir. Chinese J. Geophys. (in Chinese), 56(5): 1716-1725. DOI:10.6038/cjg20130529 | |
| Huang J P, Li Z C, Liu Y J, et al. 2013b. The least square pre-stack depth migration on complex media. Progress in Geophysics (in Chinese), 28(6): 2977-2983. DOI:10.6038/pg20130619 | |
| Huang J P, Xue Z G, Bu C C, et al. 2014a. The study of least-squares migration method based on split-step DSR. Journal of Jilin University:Earth Science Edition, 44(1): 369-374. DOI:10.13278/j.cnki.jjuese.201401306 | |
| Huang J P, Cao X L, Li Z C, et al. 2014b. Least square reverse time migration in high resolution imaging of near surface. Oil Geophysical Prospecting (in Chinese), 49(1): 107-112. | |
| Huang J P, Sun Y S, Li Z C, et al. 2014c. Least-squares split-step migration based on frequency-division encoding. Oil Geophysical Prospecting (in Chinese), 49(4): 702-707. | |
| Huang J P, Li C, Li Q Y, et al. 2015. Least-squares reverse time migration with static plane-wave encoding. Chinese J. Geophys. (in Chinese), 58(6): 2046-2056. DOI:10.6038/cjg20150619 | |
| Huang Y S, Schuster G T. 2012. Multisource least-squares migration of marine streamer and land data with frequency-division encoding. Geophysical Prospecting, 60(4): 663-680. DOI:10.1111/gpr.2012.60.issue-4 | |
| Li B, Li M, Liu H W, et al. 2012. Stability of reverse time migration in TTI media. Chinese J. Geophys. (in Chinese), 55(4): 1366-1375. DOI:10.6038/j.issn.0001-5733.2012.04.032 | |
| Li N, Huang J P, Li Z C, et al. 2014. The study on numerical simulation method of Lebedev Grid with dispersion improvement coefficients in TTI media. Chinese J. Geophys. (in Chinese), 57(1): 261-269. DOI:10.6038/cjg20140121 | |
| Li Z C, Guo Z B, Tian K. 2014. Least-squares reverse time migration in visco-acoustic medium. Chinese J. Geophys. (in Chinese), 57(1): 214-228. DOI:10.6038/cjg20140118 | |
| Liang K, Wu G C, Ying X Y, et al. 2009. Wave equation decomposition in 3-D TTI medium. Oil Geophysical Prospecting (in Chinese), 44(1): 19-27. | |
| Liu L B. 2012. Angle-domain Kirchhoff prestack time migration approach in VTI media. Journal of China University of Petroleum (in Chinese), 36(5): 51-55. | |
| Liu P, Yang C C, Wang Y F. 2013. Three-dimensional Green's function calculation in TTI media using the Gaussian Beam Method. Progress in Geophysics (in Chinese), 28(5): 2547-2553. DOI:10.6038/pg20130533 | |
| Liu Y J, Li Z C, Wu D, et al. 2013. The research on local slope constrained least-squares migration. Chinese J. Geophys. (in Chinese), 56(3): 1003-1011. DOI:10.6038/cjg20130328 | |
| Liu Y Z, Wang G Y, Yang J Z, et al. 2015. Multi-parameter full-waveform inversion for VTI media based on Born sensitivity kernels. Chinese J. Geophys. (in Chinese), 58(4): 1305-1316. DOI:10.6038/cjg20150418 | |
| Luo S, Hale D. 2014. Least-squares migration in the presence of velocity errors. Geophysics, 79(4): S153-S161. DOI:10.1190/geo2013-0374.1 | |
| Nemeth T, Wu C J, Schuster G T. 1999. Least-squares migration of incomplete reflection data. Geophysics, 64(1): 208-221. DOI:10.1190/1.1444517 | |
| Pestana R C, Ursin B, Stoffa P L. 2011. Separate P-and SV-wave equations for VTI media.//12th International Congress of the Brazilian Geophysical Society, 163-167. | |
| Plessix R E, Mulder W A. 2004. Frequency-domain finite-difference amplitude-preserving migration. Geophysical Journal International, 157(3): 975-987. DOI:10.1111/gji.2004.157.issue-3 | |
| Plessix R E. 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysical Journal International, 167(2): 495-503. DOI:10.1111/gji.2006.167.issue-2 | |
| Plessix R E, Cao Q. 2011. A parametrization study for surface seismic full waveform inversion in an acoustic vertical transversely isotropic medium. Geophysical Journal International, 185(1): 539-556. DOI:10.1111/gji.2011.185.issue-1 | |
| Symes W W. 2008. Approximate linearized inversion by optimal scaling of prestack depth migration. Geophysics, 73(2): R23-R35. DOI:10.1190/1.2836323 | |
| Tang Y X. 2009. Target-oriented wave-equation least-squares migration/inversion with phase-encoded Hessian. Geophysics, 74(6): WCA95-WCA107. DOI:10.1190/1.3204768 | |
| Tarantola A. 1984. Linearized inversion of seismic reflection data. Geophysical Prospecting, 32(6): 998-1015. DOI:10.1111/gpr.1984.32.issue-6 | |
| Thomsen L. 1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966. DOI:10.1190/1.1442051 | |
| Tsvankin I. 1996. P-wave signatures and notation for transversely isotropic media:An overview. Geophysics, 61(2): 467-483. DOI:10.1190/1.1443974 | |
| Wang J, Li Z C, Sun X D, et al. 2012. Reverse time migration in tilt transversely isotropic(TTI) media. Oil Geophysical Prospecting (in Chinese), 47(4): 573-577. | |
| Wang Y F. 2013. Comparison of interferometric migration and preconditioned regularizing least squares migration inversion methods in seismic imaging. Chinese J. Geophys. (in Chinese), 56(1): 230-238. DOI:10.6038/cjg20130123 | |
| Wu G C. 2006. Seismic Wave Propagation and Imaging in Anisotropic Media (in Chinese). Dongying:China University of Petroleum Press. | |
| Wu G C, Liang K. 2006. Quasi P-wave forward modeling in frequency-space domain in VTI media. Oil Geophysical Prospecting (in Chinese), 40(5): 535-545. | |
| Wu G C, Liang K, Wang H Z. 2008. High-order one-way generalized-screen propagation operator of qP-wave in VTI medium. Oil Geophysical Prospecting (in Chinese), 42(6): 640-650. | |
| Yoon K, Suh S, Ji J, et al. 2010. Stability and speedup issues in TTI RTM implementation.//80th Annual International Meeting, SEG, Expanded Abstracts, 3221-3225. | |
| Zhan G, Pestana R C, Stoffa P L. 2012. Decoupled equations for reverse time migration in tilted transversely isotropic media. Geophysics, 77(2): T37-T45. DOI:10.1190/geo2011-0175.1 | |
| Zhan G, Pestana R C, Stoffa P L. 2013. An efficient hybrid pseudo spectral/finite-difference scheme for solving the TTI pure P-wave equation. Journal of Geophysics and Engineering, 10(2): 025004. DOI:10.1088/1742-2132/10/2/025004 | |
| Zhang H Y, Wang Y F. 2014. Iterative regularization migration inversion in VTI media. Oil Geophysical Prospecting (in Chinese), 49(6): 1083-1090. | |
| Zhang Y, Sun J, Notfors C, et al. 2005. Delayed-shot 3D depth migration. Geophysics, 70(5): E21-E28. DOI:10.1190/1.2057980 | |
| Zhang Y, Zhang H Z. 2009. A stable TTI reverse time migration and its implementation.//79th Annual International Meeting, SEG, Expanded Abstracts. | |
| Zhang Y, Zhang H Z, Zhang G Q. 2011. A stable TTI reverse time migration and its implementation. Geophysics, 76(3): WA3-WA11. DOI:10.1190/1.3554411 | |
| Zhang Y, Wu G C. 2013. Methods of removing pseudo SV-wave artifacts in TTI media qP-wave reverse-time migration. Chinese J. Geophys. (in Chinese), 56(6): 2065-2076. DOI:10.6038/cjg20130627 | |
| Zhou H, Zhang G, Bloor R. 2006a. An anisotropic acoustic wave equation for VTI media.//68th EAGE Conference and Exhibition. SPE, EAGE. | |
| Zhou H B, Zhang G Q, Bloor R. 2006b. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media.//76th Annual International Meeting, SEG, Expanded Abstracts, 194-198. | |
| Zhu T F, Gray S H, Wang D L. 2007. Prestack Gaussian-beam depth migration in anisotropic media. Geophysics, 72(3): S133-S138. DOI:10.1190/1.2711423 | |
| 曹文俊, 李振春, 韩文功. 2007. 波动方程叠前深度偏移对速度模型的敏感性分析. 天然气工业, 27(增刊A): 120–122. | |
| 程玖兵, 康玮, 王腾飞. 2013. 各向异性介质qP波传播描述I:伪纯模式波动方程. 地球物理学报, 56(10): 3474–3486. DOI:10.6038/cjg20131022 | |
| 段鹏飞, 程玖兵, 陈爱萍, 等. 2013. TI介质局部角度域高斯束叠前深度偏移成像. 地球物理学报, 56(12): 4206–4214. DOI:10.6038/cjg20131223 | |
| 段新意, 李振春, 黄建平, 等. 2014. 各向异性介质共炮域高斯束叠前深度偏移. 石油物探, 53(5): 579–586. | |
| 韩建光, 王赟, 张晓波, 等. 2015. VTI介质高斯束叠前深度偏移. 石油地球物理勘探, 50(2): 267–273. | |
| 黄建平, 李振春, 孔雪, 等. 2013a. 碳酸盐岩裂缝型储层最小二乘偏移成像方法研究. 地球物理学报, 56(5): 1716–1725. DOI:10.6038/cjg20130529 | |
| 黄建平, 李振春, 刘玉金, 等. 2013b. 复杂介质最小二乘叠前深度偏移成像方法. 地球物理学进展, 28(6): 2977–2983. DOI:10.6038/pg20130619 | |
| 黄建平, 薛志广, 步长城, 等. 2014a. 基于裂步DSR的最小二乘偏移方法. 吉林大学学报(地球科学版), 44(1): 369–374. DOI:10.13278/j.cnki.jjuese.201401306 | |
| 黄建平, 曹晓莉, 李振春, 等. 2014b. 最小二乘逆时偏移在近地表高精度成像中的应用. 石油地球物理勘探, 49(1): 107–112. | |
| 黄建平, 孙郧松, 李振春, 等. 2014c. 一种基于分频编码的最小二乘裂步偏移方法. 石油地球物理勘探, 49(4): 702–707. | |
| 黄建平, 李闯, 李庆洋, 等. 2015. 一种基于平面波静态编码的最小二乘逆时偏移方法. 地球物理学报, 58(6): 2046–2056. DOI:10.6038/cjg20150619 | |
| 李博, 李敏, 刘红伟, 等. 2012. TTI介质有限差分逆时偏移的稳定性探讨. 地球物理学报, 55(4): 1366–1375. DOI:10.6038/j.issn.0001-5733.2012.04.032 | |
| 李娜, 黄建平, 李振春, 等. 2014. Lebedev网格改进差分系数TTI介质正演模拟方法研究. 地球物理学报, 57(1): 261–269. DOI:10.6038/cjg20140121 | |
| 李振春, 郭振波, 田坤. 2014. 黏声介质最小平方逆时偏移. 地球物理学报, 57(1): 214–228. DOI:10.6038/cjg20140118 | |
| 梁锴, 吴国忱, 印兴耀, 等. 2009. 三维TTI介质波动方程分解. 石油地球物理勘探, 44(1): 19–27. | |
| 刘立彬. 2012. VTI介质角度域Kirchhoff叠前时间偏移方法. 中国石油大学学报(自然科学版), 36(5): 51–55. | |
| 刘鹏, 杨长春, 王彦飞. 2013. 三维TTI介质格林函数的高斯束计算方法. 地球物理学进展, 28(5): 2547–2553. DOI:10.6038/pg20130533 | |
| 刘玉金, 李振春, 吴丹, 等. 2013. 局部倾角约束最小二乘偏移方法研究. 地球物理学报, 56(3): 1003–1011. DOI:10.6038/cjg20130328 | |
| 刘玉柱, 王光银, 杨积忠, 等. 2015. 基于Born敏感核函数的VTI介质多参数全波形反演. 地球物理学报, 58(4): 1305–1316. DOI:10.6038/cjg20150418 | |
| 王娟, 李振春, 孙小东, 等. 2012. TTI介质逆时偏移成像. 石油地球物理勘探, 47(4): 573–577. | |
| 王彦飞. 2013. 地震波干涉偏移及预条件正则化最小二乘偏移成像 方法对比. 地球物理学报, 56(1): 230–238. DOI:10.6038/cjg20130123 | |
| 吴国忱. 各向异性介质地震波传播与成像.东营: 中国石油大学出版社, 2006. | |
| 吴国忱, 梁锴. 2006. VTI介质频率-空间域准P波正演模拟. 石油地球物理勘探, 40(5): 535–545. | |
| 吴国忱, 梁锴, 王华忠. 2008. VTI介质qP波广义高阶屏单程传播算子. 石油地球物理勘探, 42(6): 640–650. | |
| 张洪宙, 王彦飞. 2014. VTI介质迭代正则化偏移反演算法. 石油地球物理勘探, 49(6): 1083–1090. | |
| 张岩, 吴国忱. 2013. TTI介质qP波逆时偏移中伪横波噪声压制方法. 地球物理学报, 56(6): 2065–2076. DOI:10.6038/cjg20130627 | |