2017, Vol. 60
大量岩石物理实验以及野外测量表明，各向异性广泛存在于各种岩石类型中.地震各向异性集中体现为速度各向异性，即地震波波速随传播方向变化，如果在地震资料处理中忽略各向异性效应，特别是走时误差，必然给走时层析、偏移成像以及参数建模等处理方法带来一些不可估量的错误，从而误导地质解释工作(Thomsen，1986；李振春等, 2016a, 2016b；印兴耀和刘倩，2016；刘文卿等，2016；Shi and Wang, 2016；Wu et al., 2016).为了提升各向异性地震资料处理质量，近年来，随着地震勘探的不断深入与计算机技术的迅速发展，针对各向异性的高精度地震成像和参数反演技术早已成为业界研究焦点并逐步应用于实际生产当中.
目前，最受学者关注的各向异性介质有三类：垂直横向各向同性介质(VTI，Vertical Transversely Isotropic meida)、水平横向各向同性介质(HTI，Horizontal Transversely Isotropic meida)和倾斜横向各向同性介质(TTI，Tilted Transverse Isotropic media)，统称为TI介质(Transversely Isotropic meida)，尽管各类介质的地质成因差异较大，但在等效理论框架下，其介质参数仅仅归结为倾角不同.与各向同性介质类似，在TI介质中主要存在两类波型，即纵波与横波，但偏振特征有所区别，除了特定方向之外，纵波的极化方向与传播方向不平行，横波的极化方向与传播方向也不垂直，因此称纵波为准纵波(即qP波)，称横波为准横波，此外，在各向异性介质中，两类横波的传播速度并不一致，因此又分别称为快、慢横波(也即qSV波和SH波)，这一现象称为横波分裂.正演模拟是认识地震波传播规律的有效手段，也是实现逆时偏移(RTM，Reverse Time Migration)、波动方程层析以及全波形反演(FWI，Full Waveform Inversion)的核心，各向异性全弹性波方程可准确且全面地描述各类地震波型在地层中的传播情况，但是全弹性波正演、成像与反演还面临诸多困难：①野外多分量地震数据采集十分昂贵；②需要输入更多的背景参数，如横波速度、横波各向异性强度参数等等，而以当前的技术，这些参数尚难完全准确地获取；③计算成本过高，I/O及内存需求较大，算法设计复杂，增加了成像和反演的难度；④波型分离与分解非常复杂，存在分离不彻底或计算量大等问题(Yong et al., 2016).因此，在TI介质正演模拟、逆时偏移以及参数反演中仍然采用简化的qP波方程为主.
众所周知，简单实用的标量纵波方程和高效稳健的计算方法是编制地震波正演与逆时偏移成像软件的基础，但在各向异性介质中始终存在纵横波耦合，很难导出仅仅刻画qP波运动学与动力学特征的精确纯qP波方程.为此，多年来许多学者提出了一系列理论：Thomsen (1986)首先提出了弱各向异性近似理论，随后，Tsvankin (1996)推导了用Thomsen参数表征的VTI介质qP-qSV波精确相速度公式，为了简化该式，Alkhalifah (1998)提出了著名的“声学近似”，这些成果成为后续各向异性研究的理论基础，在构建qP波控制方程、推动逆时偏移技术应用等方面发挥了重要作用.追溯qP波方程的简化历程，可大致归为两类：第一类是qP-qSV波耦合波动方程，Alkhalifah (2000)在声学近似下，基于qP波精确频散方程首次导出了VTI介质四阶波动方程，但四阶方程存在两组不合物理规律的人为假解，而且高阶时-空混合偏导数求解困难，因此，Zhou等(2006a)通过引入辅助函数将其简化为二阶耦合波动方程组，随后采用局部旋转坐标变换将其扩展到TTI介质并应用于逆时偏移成像中(Zhou et al., 2006b)，该方程在倾角急剧变化的区域中存在严重的不稳定现象，且不稳定性随传播时间迅速加剧.针对不稳定问题，一种解决方案是重新引入非零横波速度项来修正上述方程组，以减缓由声学近似造成的不稳定，但增强了伪横波噪音强度，同时计算效率降低，另一种方案是将模型不稳定区域设置为椭圆各向异性，但改变了模型的局部特征(Fletcher et al., 2009；Yoon et al., 2010)，对此，很多专家认为仅从运动学频散关系出发构建qP波正演算子，不足以准确描述qP波动力学特征，难以克服方程固有的不稳定性(Duveneck et al., 2008；Zhang and Zhang, 2009a；徐文才等，2016)，Zhang等(2011)与Duveneck和Bakker (2011)分别从虎克定律出发，并结合声学近似，推导出一种新的稳定的波动方程，并证明该类方程在正演和成像过程中能量始终保持稳定，程玖兵等(2013)则借助相似变换，推导出在动力学上qP波占优的伪纯模式波动方程.总之，沿用声学近似并不能完全消除非对称轴方向的剪切波能量，方程中必然存在纵横波耦合及由此造成的数值频散和不稳定现象，很难准确成像，而且只适用于η≥0的各向异性区域.
第二类是完全不含伪横波干扰的TI介质纯qP波控制方程，要实现qP与qSV波的完全分离，关键在于推导简单、灵活的纯qP控制方程和设计有效的拟微分求解算法.Klié和Toro (2001)最早开始讨论qP与qSV波的解耦问题，其目的是压制Alkhalifah (2000)提出的四阶波动方程的人为假解，并缓解由于解的指数增长所导致的不稳定问题，随后，Du等(2005)和Zhang等(2005)基于弱各向异性近似和平方根近似，推导出一种在运动学上较为准确的TTI介质时间-波数域纯qP波解耦方程，如果变换到时空域，即表现为拟微分形式，因此，主流的有限差分解法不再适用，为了处理该拟微分算子，先后发展了多种算法：第一，时空域最佳分离近似(OSA，Optimized Separable Approximation)方案(Song，2001；Chen and Liu, 2004；Liu et al., 2009)，此方法计算精度较高但分离成本十分昂贵；第二，快速展开法(REM，Rapid Expansion Method)，正广泛应用于VTI介质逆时偏移成像中(Pestana and Stoffa, 2010；Pestana et al., 2011)；第三，基于REM的TTI介质混合有限差分与伪谱法，该方法需进行多次傅氏变换(二维7次，三维22次)和插值处理，效率很低(Zhan et al., 2012, 2013；黄建平等，2016；王璐琛等，2016)；第四，拟微分算子分解策略，即将传播算子分解成易于处理的微分算子和标量算子，然后分别求解，该策略计算效率较高，并且被成功推广至TTI和TOA介质偏移成像中(Sheng and Zhou, 2014a，2014b).近年来，Alkhalifah等(2013)另辟蹊径，通过qP波走时场建立等效的各向同性模型，进而计算纯qP波场，此外，Chu等(2013)导出一种不含分数阶算子的时空域纯qP波近似方程，并构建了经典的有限差分法计算方案.综上所述，该类方程能够较为准确地描述qP波走时特征、消除伪横波噪音并且避免不稳定现象，但是方程推导过程中需要进行两次或多次近似，动力学特征难以保证，而且求解过程中计算成本与内存需求也较高.
为了克服上述纯qP波模拟存在的诸多缺陷，兼顾计算效率与精度，本文通过借鉴当前较为新颖的波场延拓工具--Low-rank分解，来设计纯qP波正、反向延拓算法.目前，对Low-rank近似的研究主要集中在TCCS (Texas Consortium for Computational Seismology)小组，Fomel等(2010, 2013)首次将Low-rank分解引入到地震波场延拓中，并指出该算法在构建各类介质正演与成像算子、降低计算成本等方面拥有巨大的潜力.紧接着，相继提出了适用于各向同性与TTI介质的高精度Low-rank有限差分及Low-rank傅里叶有限差分正演模拟方法，与传统的有限差分和傅里叶有限差分相比，这两类方法的波场传播更加稳定，且不存在数值频散，还能直接处理变密度情形(Song et al., 2011, 2013).随后，Song和Alkhalifah (2013)又将Low-rank分解扩展到正交各向异性介质正演模拟中，获得了理想的纯qP波场.为了继承交错网格的特点，Fang等(2013, 2014, 2015)及方刚等(2014)则提出将Low-rank分解策略与交错网格高阶有限差分相结合，先后实现了交错网格Low-rank (SGL，Staggered-Grid Lowrank Method)正演和交错网格Low-rank有限差分(SGLFD，Staggered-Grid Lowrank Finite-Differences Method)正演模拟方法，目前，该算法只用于各向同性介质中.此外，鉴于Low-rank分解在处理空间-波数混合域问题中的天然优势，Cheng和Fomel (2014)将其应用于各向异性弹性波波型分离中，取得良好的分离效果.近年来，Sun等(2015)又发展了一种基于Low-rank近似的黏介质正演与成像方法，较为准确地刻画了地震波振幅衰减特征.从已出版的文献中可见，Low-rank分解在复杂TTI介质成像处理中应用较少.
本文借鉴前人经验，首先构建了复杂VTI与TTI介质纯qP波传播算子，然后借助Low-rank分解策略求取该算子中的空间-波数混合域矩阵，同时引入Cerjan衰减边界条件压制人工边界反射，最终建立了复杂TTI介质高精度正演模拟与逆时偏移成像流程，在实现方法的基础上，通过开展均匀模型正演试验及复杂TTI介质成像测试，展示了本方法在刻画纯qP波运动学特征方面的优势及其在成像应用中的潜力.
2 方法原理 2.1 纯qP波传播算子在各向异性介质中，构建一种既简单又精确的纯qP波传播算子是实现TI介质高精度正演模拟、逆时偏移以及全波形反演的首要任务.仿照经典的各向同性伪解析法原理(Zhang and Zhang, 2009b；Etgen and Brandsberg-Dahl, 2009；Crawley et al., 2010)，当介质速度与各向异性参数随着空间位置变化时，地震波场时间外推公式可由空间-波数混合域相移算子来近似表征：
|
(1) |
其中，t是时间坐标，x是空间坐标，k是空间波矢量；Φ代表角频率函数，由介质参数及传播方向共同决定；P(k, t)与P(k, t+Δt)分别是相邻两时刻的波数域地震波场，Δt表示时间步长；在点震源的情况下，±号分别表示会聚波和发散波.本文中，空间傅氏变换对满足如下定义：
|
(2) |
公式(2)中，p和P分别表示傅氏变换前后的空间域和波数域地震波场.
从公式(1)可知，角频率函数是决定延拓算子性质的关键，为了构建数值模拟可行的纯纵波延拓算子，将公式(1)代入TI介质波动方程，在高频近似下可导出角频率关系式(Fomel，2004)，即
|
(3) |
公式(3)中，符号v表示地震波相速度，对于各向同性介质，v仅为空间位置坐标的函数，而在TI介质中，由于地震波传播速度随传播方向变化，因此，v不仅与空间位置有关，还是波矢量k的函数，可见TI介质相速度本身即为空间-波数混合域形式；▽表示梯度算子，如果相速度梯度▽ v或时间步长Δt足够小，可以进一步忽略高阶项，于是方程(3)可简化为v(x, k)| k.
为了进一步推导基于Thomsen参数表征的角频率公式，可借助于qP波频散关系，在空间-波数域，VTI介质qP波精确频散关系可表示为(Tsvankin，1996；Fowler et al., 2003)
|
(4) |
公式(4)中，vv表示沿着对称轴方向的qP波相速度，vh表示在各向同性面内的qP波相速度，
|
(5) |
公式(5)中，
|
(6) |
最后，将公式(6)代入公式(1)中即可得到VTI介质纯qP波传播算子.值得一提的是，本文采用Low-rank分解地震波延拓方案可以直接处理上述空间-波数混合域算子，无需推导显式的时空域波动方程，且在横波速度已知的前提下，可进一步采用精确频散方程(即公式(4))来构建纯qP波传播算子以提高计算精度.
接下来，考察TTI介质.受构造运动的影响，TI介质对称轴并不总是竖直向下，相对于观测坐标系，经常存在倾斜与旋转，即极化各向异性和方位各向异性两种情况.为方便数学描述，首先定义如图 1所示的观测坐标系，倾角θ为介质对称轴方向与z轴的夹角，并规定绕y轴以逆时针方向旋转为正，方位角φ为对称轴方向在xoy面内的投影与x轴的夹角，且在图中以顺时针方向旋转为正，根据上述规定，TTI介质对称轴方向矢量可表示为
|
(7) |
|
图 1 观测坐标系 Fig. 1 Survey coordinates |
借助局部旋转坐标变换可以将VTI介质纯qP波传播算子扩展到TTI的情形，即采用新的空间波数分量
|
(8) |
在二维(2D)情况下，只存在极化倾角，令方位角为0，且不考虑y方向，公式(8)可简化为
|
(9) |
联立公式(1)、(6)、(8)，即可得到TTI介质纯qP波传播算子.如果直接采用一步法延拓求解上述纯qP波算子，虽然波场延拓过程稳定，且人工边界易于处理，但是该方法涉及复矩阵分解，在时间步长一定的情况下，一步法延拓需耗费更多的计算量和存储空间(Du et al., 2014).为了避免复数处理，本文采用简化的两步法时间延拓方案，下面讨论纯qP波两步法延拓公式.
所谓两步法时间延拓，即为欲求下一时刻的地震波场，必须预先知道前两个时刻的波场.首先基于公式(1)可得一步法逆向时间延拓公式为
|
(10) |
其次，将公式(1)与公式(10)左右两端对应相加，再将其反变换到空间域，可得两步法时间延拓公式(Tal-Ezer，1986；Etgen, 1986, 1989；Soubaras and Zhang, 2008)，即
|
(11) |
对于均匀介质，公式(11)无条件稳定，即取任意时间步长，都能保证地震波稳定地传播；而对于非均匀或各向异性介质，只要时间步长Δt足够小，也能提供较为合理的数值近似解(Fomel et al., 2013).
观察公式(11)可知，方程右端被积函数中的空间-波数混合域算子刻画了纯qP波所有传播特征，本文称之为纯qP波传播矩阵，记为
|
(12) |
方程(11)可直接采用高维空间傅氏变换进行求解，但对复杂介质，在时间方向每延拓一步，需进行一次傅氏正变换(FFT，Fast Fourier Transforms)和多次(Nx×Nz次，Nx、Nz表示模型网格点数)傅氏反变换(IFFT，Inverse Fast Fourier Transforms)，计算成本极其昂贵，不利于大规模复杂模型地震波计算；另一种是采取基于Low-rank分解的波场延拓方案，可将多次IFFT减少为一次或几次IFFT，从而降低计算成本.下一节将探讨Low-rank分解延拓算法.
2.2 Low-rank延拓算法高效稳健的计算方法是设计地震成像处理软件的核心，也是提高计算效率的关键.Low-rank分解的基本思想是从矩阵W中分别选择一组具有代表性的空间位置x和一组具有代表性的空间波数k，去近似表征全矩阵W，从而降低计算成本，它本质上是一种广义的插值算法，在插值中通过Low-rank近似来选择最佳的参考速度和权值(Fomel et al., 2010；Song et al., 2013；Wu and Alkhalifah, 2014；Yan and Liu, 2016).
为实现空间-波数混合域算子的分离与求解，通过Low-rank分解将传播矩阵W近似分解为三个子矩阵，即
|
(13) |
其中，W 1、W 2都是全矩阵W的子矩阵，W 1由矩阵W中M个线性无关的列向量组成，W 2由矩阵W中N个线性无关的行向量组成，amn是连接W 1与W 2的系数矩阵，也是实现矩阵最佳分离的权系数.
根据Fomel等(2010)建立的Low-rank矩阵分解流程，在模型给定的情况下，通过预先设置容许误差和时间步长Δt，即可求出上述两个子矩阵和系数矩阵，以及对应子矩阵的秩M和N.本文在这里不再详细讨论其解法，具体实现过程与数学证明可查阅(Fomel et al., 2013).最后，将公式(13)带入公式(11)，即可导出基于Low-rank分解的两步法时间延拓公式：
|
(14) |
|
图 2 基于Low-rank分解的纯qP波场延拓算法流程图 Fig. 2 The flowchart of using Low-rank decomposition algorithm to compute pure qP wavefield |
从公式(14)可见：Low-rank分解是一种解决空间-波数混合域算子的有效工具，为伪解析法波场延拓提供了新的思路.需要强调的是，秩N的大小决定时间层每延拓一步需执行的IFFT次数，也即控制着计算成本的高低，秩越低延拓速度越快.事实上，M表示选取的具有代表性的特定波数的个数，N表示选取的具有特定介质参数大小的空间位置个数，其大小都受控于介质结构的复杂程度、纵横向参数变化情况和预定误差大小等因素，介质结构越复杂、参数变化越剧烈、容许误差越小，值就越大，这是因为模型越复杂，必将导致更多线性无关的空间位置和波数位置，要在给定误差范围内准确表征原始混合域矩阵W，必然需要更大的M和N值(Song and Alkhalifah, 2013).
根据公式(14)可设计如下基于Low-rank分解的纯qP波场延拓算法流程：
从上述算法流程可以看出：①在时间方向每递推一步，需要执行一次FFT和N次IFFT运算，以及必要的乘法与求和运算，这些运算构成了波场延拓算法的主要计算量；②正是因为Low-rank分解方法兼顾了所有的波数成分和空间位置，特别是高波数段、高速岩体、参数剧变区域，因此在计算过程中不会产生数值频散和不稳定现象；此外，③对频散关系始终未进行平方运算，也不必推导显式的波动方程，所以方程不会产生增根，也就避免了伪横波干扰的产生以及对模型参数的限制.
2.3 逆时偏移流程逆时偏移是一种基于双程波动方程的精确成像方式，不受反射界面倾角的限制，在偏移速度较为可靠的前提下，能够同时使回转波和多次波准确成像.实现基于Low-rank分解的各向异性逆时偏移成像主要包括三个步骤，即震源波场正向延拓、检波点波场反向延拓以及零延迟互相关提取成像剖面.
首先，以震源函数作为初始条件，通过地震波沿时间正向延拓构建震源波场，其定解问题可描述为
|
(15) |
式中，pf(x, t)表示t时刻的震源波场，即正传波场，f(t)代表震源子波，x s为震源坐标.
其次，以炮记录作为边界条件，通过地震波沿时间逆向外推构建检波点波场，其定解问题可表示为
|
(16) |
其中，pb(x, t)表示t时刻的检波点波场，也即反传波场，r(x g, t)代表地震炮记录，x g为接收点坐标.
最后，对震源波场和检波点波场进行零延迟互相关，即可提取RTM成像剖面m，即
|
(17) |
由上述流程可知：在两步法时间延拓方案中，地震波正、反向延拓过程共享同一传播矩阵，即只需进行一次Low-rank分解，就能实现地震波正、逆双向外推，然而在一步法延拓中需要分别对正、反向传播矩阵进行Low-rank分解，增加了额外计算量与存储空间.
此外，对于正演模拟与逆时偏移成像，都面临人工边界反射问题.完全匹配层(PML)边界条件具有最好的边界吸收效果(王永刚等，2007)，但由于边界层(特指含有衰减因子的扩展层)的引入，传播矩阵变得更加复杂，导致秩N变大，从而成倍地增加计算量.为此，本文统一采用Cerjan指数衰减边界条件解决边界反射，即直接对含有扩展层的速度模型进行Low-rank分解，并在波场延拓过程中对扩展层内的波场值施加指数衰减.
3 模型测试/数值模拟本文分别采用三类模型测试来检验基于Low-rank分解的纯qP波正演模拟及逆时偏移方法的正确性和适应性.第一类是四组均匀模型对比试验，第二类是简单TTI模型分析试验，第三类是BP2007复杂模型成像试处理.本文观测系统统一采取成本较低的地面反射波接收方式，通过Cerjan衰减边界条件处理人工边界反射.在成像处理中，使用震源归一化互相关成像条件来补偿中深部成像能量.
3.1 均匀模型 3.1.1 模型参数首先，为了验证纯qP波传播算子及Low-rank延拓算法的正确性和有效性，本文根据控制变量法设计了四组典型模型正演模拟对比试验.模型参数如表 1所示，前1~3组模型具有相同的速度和各向异性参数，但倾角不同，依次为0°、90°、30°，分别代表VTI、HTI和TTI介质，即在等效介质理论框架下，不同类型的TI介质只存在倾角的差异，第4组与第3组类似，只是调换了各向异性参数ε与δ，用于模拟ε < δ的情况.模型网格点数为701×701，纵横向采样间距是10 m×10 m，模型大小总计7.0 km×7.0 km，震源函数选用主频为30 Hz的雷克子波，震源位置处于模型中心(3.5 km，3.5 km)处，时间采样间距是2.0 ms.
|
表 1 典型参数 Table 1 Parameters |
图 3展示了不同方法不同介质在1.0 s时刻的波场快照.其中，图 3a-3d分别对应上述四组介质模拟结果，正演方程采用声学近似下的耦合波动方程(具体形式可见文献Zhou et al., 2006b中公式(7)和(8))，图中可见：①由于各模型倾角不同，波前面出现相应旋转，且旋转角度等于极化倾角(如图 3c所示)；②在qP波正演快照中，还存在冗余的低速伪横波假象，其响应曲线表现为能量较弱的菱形状，这与声学近似(令对称轴方向的横波速度等于零)下横波波前形态一致；③当各向异性Thomsen参数ε < δ时，波场计算不稳定，该类方法不再适用.类似地，图 3e-3h是选取带非零横波速度项的耦合波动方程(具体形式可见文献Fletcher et al., 2009中公式(2))计算的结果，对照可知：①与图 3a-3d有相似的结论，不同之处在于低速横波形态由菱形变为方形，这是由于引入非零横波所致，其目的是为了缓解复杂TI介质波场延拓不稳定；②该类方程仍不能克服Thomsen参数限制.同样地，Low-rank纯qP波延拓方法的正演结果如图 3i-3l所示，经观察对比可得：① qP波走时特征与上述两种方法的计算结果类似，然而残余伪横波假象被彻底消除；②在ε < δ的情况下，波场延拓依然稳定，也即克服了对模型参数的限制.为了进一步考察上述纯qP波正演结果的正确性，图 3m-3p分别给出了对应四种参数的群速度理论曲线(图中，白色实线代表qP波，黑色实线代表qSV波)，可以发现：计算结果与理论曲线基本吻合，表明本方法能够较为准确地描述qP波运动学特征.
|
图 3 均匀TI介质在t=1 s时刻的波场快照及理论波前曲线 (a)-(d)采用声学近似(vs=0)下的耦合波动方程计算结果；(e)-(h)采用带非零横波速度项(vs=0.577vv)的耦合波动方程计算结果；(i)-(l)本文方法计算结果；(m)-(p) qP波与qSV波的理论波前面曲线. Fig. 3 Wave-field snapshots at time t=1 s in homogeneous TI medium and theoretical wavefront curves (a)-(d) Wave-field snapshots are generated from coupled equation with vs=0; (e)-(h) Wave-field snapshots are also generated from coupled equation but with nonzero SV-wave velocity vs=0.577 vv; (i)-(l) Wave-field snapshots are computed using our method; (m)-(p) Theoretical wavefront of a qP-wave & a qSV-wave. |
通过对比试验可以看到：耦合型波动方程正演模拟始终存在低速伪横波假象，且不能适用于ε < δ的情形，而本文Low-rank纯qP波延拓方法既能准确表征qP波传播规律，又能消除伪横波干扰且克服模型参数限制.
3.1.3 偏移脉冲响应测试分别选取表 1中第Ⅰ、Ⅲ、Ⅳ三组代表性模型进行逆时偏移脉冲响应测试，用以分析说明伪横波假象对成像结果的影响以及各成像算子的优缺点.测试结果如图 4所示，其中第1~3各行所采取的正演方程依次为声学近似下的耦合波动方程、带非零横波速度项的耦合波动方程和Low-rank纯qP波延拓算子.经仔细对比可以得出：①与各向同性RTM类似，脉冲响应在大倾角位置几乎不存在能量衰减和形变(如图 4a虚线箭头所示)，表明RTM没有倾角限制，能够准确成像陡倾构造，然而不同之处在于TI介质偏移脉冲呈椭圆或方形特征，其形态取决于各向异性参数和倾角，在成像中该形态决定着同相轴成像位置；②在耦合型方程的qP波脉冲响应图像中(如图 4a, b, d, e)，除qP波有效能量以外，还含有较强的串扰噪音及SV波频散假象(频散主要见于图 4a, 4b)，这是由于震源波场与接收点波场中的残余qSV波能量所致，进一步可将串扰噪音分为qSV与qSV自相关和qSV与qP互相关两种类型(如图 4(a, d, b, e)实线剪头所指)，因此，伪横波干扰势必降低成像质量；③通过开展本文纯qP波逆时偏移脉冲响应测试，串扰噪音被彻底消除，信噪比显著提高(如图 4g, h所示)；④对于同一类介质，不同偏移方法所求取的qP波脉冲响应基本一致(如图 4a, d, g或图 4b, e, h)，这也从另一个侧面证明了本方法的正确性；⑤然而，当模型参数ε < δ时，耦合型方程不再适用(如图 4c, f)，但本方法仍能获得清晰的纯qP波脉冲响应(如图 4i)，体现了Low-rank纯qP波偏移方法的优势.
|
图 4 TI介质qP波逆时偏移脉冲响应测试 (a)-(c)声学近似(vs=0)下的耦合波动方程；(d)-(f)带非零横波速度项(vs=0.577vv)的耦合波动方程；(g)-(i)本文方法. Fig. 4 Impulse response of qP-wave RTM in TI media (a)-(c) Coupled equation by setting SV-wave velocity vs=0; (d)-(f) Coupled equation by setting nonzero SV-wave velocity vs=0.577vv; (g)-(i) Our approach. |
综上所述，地震各向异性造成偏移脉冲形态畸变，因此，在成像中如果不考虑各向异性效应，可能会产生偏移不足或偏移过量，形成偏移噪音与偏移假象，甚至歪曲同相轴成像位置，造成成像紊乱.对于各向异性耦合型方程逆偏移，存在伪横波干扰及模型参数限制，而Low-rank纯qP波方法克服了上述缺陷.
3.2 简单TTI模型 3.2.1 模型参数其次，通过简单TTI模型正演与成像测试，分析探讨本文算法的稳定性及特点.测试模型选用Overthrust模型，模型参数如图 5所示，图 5a-5d依次为沿着对称轴方向的qP波速度、极化倾角模型、Thomsen参数ε及δ.网格大小301×301，网格间距10 m×6 m，模型范围总计3.0 km×1.8 km，浅层纵波速度为2.0 km·s-1，高速地层可达4.0 km·s-1.该模型整体构造较为简单，有利于观察地震波传播规律，特别是运动学特性；此外，模型中部区域各向异性特征明显，且倾角变化剧烈，是测试分析成像方法稳定性的标准模型.
|
图 5 二维TTI Overthrust模型参数 (a)纵波速度场；(b)倾角模型；(c) ε模型；(d) δ模型. Fig. 5 2D TTI Overthrust model (a) Velocity field; (b) Tilt angle θ; (c) Thomsen parameter ε and (d) δ. |
根据公式(14)，基于Low-rank分解的纯qP波波场延拓可通过如下三个步骤实现：首先，利用公式(12)合成Overthrust模型的真实传播矩阵W，抽取空间位置为(1.5 km，0.9 km)处(vv=3.0 km·s-1，θ=50°，ε=0.15，δ=0.08)的真实传播矩阵，如图 6a所示，图中可见：矩阵元素大小随波矢量模的增加而减小.随后，按照Fomel等(2010)建立的Low-rank矩阵分解流程对上述传播矩阵进行分解，形成三个规模较小的子矩阵，本文预先设定容许误差为1.0×10-4，时间步长为0.5 ms，经计算，可求出三个子矩阵及对应子矩阵的秩M=N=7，这表明在地震波传播过程中，时间方向每延拓一步，需耗费1次傅氏正变换(FFT)和7次傅氏反变换(IFFTs)；矩阵分解精度决定着地震波计算精度，为了直观显示Low-rank矩阵分解的准确程度，本文通过公式(13)重构传播矩阵.图 6b展示了对应空间位置处的重构结果，对比可知：利用子矩阵重构的近似结果接近于真实值；另外，图 6c给出了两者差剖面，可见矩阵元素的误差水平小于6.0×10-8，证实了Low-rank矩阵分解的正确性.需要指出的是，在具体实现过程中，上述传播矩阵的合成与分解是交替进行的，其目的在于避免大型传播矩阵的存储.
|
图 6 空间位置(1.5 km，0.9 km)处的传播矩阵W (x，k)=2cos[((x，k)(t] (a)真实矩阵；(b)重构矩阵；(c) Low-rank近似误差.该位置处模型参数分别是：vp=3000 m·s-1，eps=0.15，dlt=0.08，θ=50°. Fig. 6 The propagation matrix at location (1.5 km, 0.9 km) (a) True matrix; (b) Reconstructed matrix; (c) Error for the Low-rank approximation. |
此外，由于秩N的大小决定计算成本的高低，且受控于速度和各向异性参数结构的复杂程度、预定误差水平及时间步长，为此，本文详细给出了N随时间步长Δt及误差水平的变化情况如表 2所示，仔细分析表 2可以发现：同一模型中，在误差水平一定的条件下，随着时间步长的增加，秩N的大小逐步增大，计算成本也随之加大；同一模型中，如果给定时间步长，秩N的大小随容许误差的减小而增大，计算效率随之降低；因此，在实际生产中，选择合适的时间步长与容许误差应当合理权衡计算成本与模拟精度.
|
|
表 2 秩N随时间步长或误差水平变化情况 Table 2 Rank N variation with time step size Δt or error level |
最后，选定主频为25 Hz的雷克子波作为震源子波，并根据图 2所示的Low-rank延拓算法流程进行纯qP波正演模拟.
3.2.3 波场快照与炮记录上述模拟所得0.7 s时刻的正演快照如图 7a所示，从图中可以看到：地震波场清晰稳定，直达波、反射波、透射波、折射波以及绕射波等各种波型丰富齐全，且无数值频散与伪横波假象.为了暴露常规耦合型波动方程高阶有限差分正演模拟的不足，同时展现本文方法在波场延拓中的优越性，图 7b给出了声学近似下耦合型波动方程正演结果，仔细观察可见：qP波运动学特征得到较为准确的刻画，但在强各向异性及倾角急剧变化的区域，地震波出现明显的不稳定现象(如图 7b实线剪头所指)；为了缓解上述不稳定问题，学者提出引入非零横波速度项，其改进方程的模拟结果如图 7c所示，对比可知：虽然其不稳定性得到一定地改善，但qSV波噪音能量也显著增强(如图 7c虚线箭头所指)，在成像中，此类噪音表现为错位的同相轴干扰，并穿插于真实同相轴之间形成偏移假象或偏移噪音，最终降低成像质量，因此，要实现各向异性高精度成像，对qSV波必须予以消除.
|
图 7 不同方法波场快照对比 (a) Low-rank纯qP波正演；(b)声学近似下耦合型波动方程；(c)带非零横波速度项的耦合型波动方程. Fig. 7 Comparison of wave-field snapshots calculated with different methods (a) Using pure qP-wave modeling method based on Low-rank decomposition; (b) Using coupled equation by setting vs=0; (c) Using coupled equation by setting vs=0.577 vv. |
为了进一步测试本文方法对大时间步长及高频震源的适应性，图 8a-8c分别展示了不同参数下的正演结果，其主频从a到c依次为25、60、60 Hz，而对应时间步长分别为1、1、2 ms.对比图 8a与8b可以发现：随着地震子波主频增大，空间子波的分辨率明显提高，图中表现为子波形态更为精细，也即空间子波的波长显著减小，这符合频率与波长之间的反比例关系；经计算，当子波主频为60 Hz，网格间距为10 m时，近地表地震波长减小为33.3 m，因此每个波长内所采的网格点数为3.3个，如果采用常规高阶有限差分进行正演模拟，势必产生严重的数值频散假象，从而降低模拟精度；然而本文方法计算结果清楚干净，没有出现频散噪音，表明本方法能够在不减小网格间距的前提下进行高频震源的地震波高精度正演模拟，该特点非常适用于海上高频地震资料处理.
|
图 8 不同参数Low-rank纯qP波波场快照对比 (a) f=25 Hz, Δt=1 ms；(b) f=60 Hz, Δt=1 ms；(c) f=60 Hz, Δt=2 ms. Fig. 8 Comparison of wave-field snapshots from Low-rank pure qP-wave modeling methods with |
对比图 8b和8c可以看出：在子波主频一定的情况下，采用2 ms时间步长与采用1 ms时间步长所计算的地震波场精度基本一致，而计算效率提高了1倍；根据CFL (Courant-Friedrichs-Lewy)稳定性条件，代入最大速度与网格间距可得出时间步长应满足Δt≤0.75 ms，即有限差分法正演模拟的时间步长不能超过0.75 ms，这说明本文方法优于有限差分，能够在较大时间步长的情况下保证地震波稳定传播.基于该特点，在正演及成像处理中，可以适当增大时间步长来显著提高计算效率，减少计算成本，且不损失计算精度.
上面分析了某一时刻的地震波波形特征，接下来讨论近地表质点的振动规律.本次试验采用301道地表满接收观测方式，共合成51炮地震记录，道间距10 m，炮间距60 m，炮点井深12 m，总记录时间为2.4 s，时间步长0.5 ms，每道4801个采样点.其中，第1炮地震记录如图 9a所示，从该记录中可以清楚地看到：记录波场的同相轴清晰稳定，反射折射等物理特征明显，没有出现频散假象，而不足之处是记录中伴有轻微的边界反射干扰(如图 9a箭头标注)，这是由于Cerjan衰减边界条件吸收衰减不完全造成；图 9b是对应声学近似下耦合型波动方程的正演记录，很显然，随着传播时间的增加，耦合型方程正演记录逐渐出现不稳定现象，且不稳定性随时间增大而迅速加剧，在1.2 s以后，真实反射波几乎被噪声完全掩盖；为了抑制波场传播过程发生畸变，通过增加非零横波速度项来缓和方程的不稳定性，图 8c展示了改进方程的模拟记录，从图中可见：加入非零横波速度项抑制了振幅畸变，波场中没有出现异常振幅，各波型特征与图 8a接近，但该记录总体上模糊不清，信噪比低，分辨率差.
|
图 9 不同方法正演记录对比 (a)-(c)的计算方法分别与图 7(a)-(c)一致. Fig. 9 Comparison of shot gathers calculated with different methods The methods of (a)-(f) are consistent with Fig.7 (a)-(f), respectively. |
通过以上对比分析可知：本文Low-rank纯qP波正演模拟方法能较准确地模拟强TTI介质中地震波传播过程，且能在一定程度上改善常规二阶耦合型方程正演模拟的不稳定问题.下面进行偏移成像测试，为了更好地成像，本文切除了能量较强的直达波.
3.2.4 成像测试各种方法成像结果分别如图 10a-10d所示.图 10a为各向同性方法成像结果，可以看出：由于忽略强各向异性造成的地震波走时畸变，导致成像剖面中TTI异常体下覆地层的同相轴连续性变差(如图 10a虚线框标示)，而TTI区域的倾斜同相轴出现扭曲(如图 10a箭头所指)，总体而言，模型的构造位置成像欠佳且噪声干扰严重.基于耦合型方程的成像结果如图 10b所示，对比可知：通过各向异性逆时偏移矫正了各向同性方法的不足，使得反射波正确归位，地下构造准确成像，特别是遭受各向异性影响的倾斜地层及深部反射层都恢复到真实的空间位置；但剖面中出现了虚假同相轴及严重的偏移噪音，这是因为残余qSV波能量弥散于真实剖面之中互相干涉，造成成像紊乱.
|
图 10 各种偏移方法成像结果对比 (a)传统各向同性逆时偏移结果；(b)耦合型方程各向异性逆时偏移结果；(c)本方法逆时偏移结果；(d)各向异性高斯束偏移结果. Fig. 10 Comparison of migration results obtained using different approaches (a) Conventional acoustic RTM image; (b) TTI RTM image by coupled qP wave equations; (c) TTI RTM image by our methods; (d) Anisotropic Gaussian beam migration image. |
然而，通过采用本文Low-rank纯qP波逆时偏移方法，基本消除了伪横波造成的偏移假象与偏移噪音，成像结果如图 10c所示，观察可见：剖面同相轴更加连续，振幅更加均衡，信噪比明显提高，成像质量显著改善.此外，计算效率是制约方法实用化的一个重要因素，在此，我们将本方法与TTI介质高斯束成像策略进行了对比，高斯束成像剖面如图 10d所示，从中可以发现：该方法成像效果与图 10c相近，但深部同相轴聚集性略显不足，其单炮计算时间约为77.281 s，而本方法单炮成像时间约为754.249 s，其计算效率较低，原因在于逆时偏移成像算子的计算时间比射线类长，此外，Low-rank波场延拓方式较有限差分法实现需要更多的计算量.上述测试结果表明：各向异性主要引起地震波相位畸变(即走时误差)，忽略各向异性效应势必造成偏移不足；TTI介质高斯束偏移兼顾了计算效率与成像精度，而本文方法具有最高的成像精度，是复杂各向异性地质体准确成像的有力工具.
3.3 复杂BP2007模型 3.3.1 BP2007模型最后，在验证方法正确性的基础上，将上述纯qP波叠前逆时偏移成像流程应用于复杂BP2007模型成像试处理，以进一步验证本方法对更复杂模型的适用性.
图 11a-11d分别展示了TTI介质BP2007模型的纵波速度场、倾角模型、各向异性Thomsen参数ε及δ.该模型在沉积地层的背景下受高速岩体侵入，高陡构造及断裂系统发育，反射结构复杂；在侵入岩体侧翼，由于遭受强烈的构造挤压，沉积产生的横向各向同性地层(VTI)被改造，其对称轴方向在局部挤压应力的作用下发生倾斜，最终形成极化各向异性，宏观上呈现出TTI介质特征.模型网格数为6298×901，网格间距Δx=Δz=12.5 m，因此，模型横向距离是78712.5 m，垂向深度为11250.0 m，此外，浅水层的速度是1492 m·s-1.
|
图 11 2D TTI BP2007模型参数 (a)速度模型；(b)倾角模型；(c) ε模型；(d) δ模型. Fig. 11 2D TTI BP2007 model (a) Model velocity field vv; (b) Model dip parameter θ; (c) Model Thomsen parameter ε and (d) δ. |
为了使模型边界充分照明，模型左右两端分别扩充了802道和263道，也即模型范围扩展成：-10025~82000 m.观测系统采用单边接收方式：BP2007模型共有1641炮地震数据，第一炮炮点位置是x=0 m，激发井深6 m，炮点间距50 m；每炮均为800道接收，接收点深度8 m，近偏移距是-37.5 m，远偏移距是-10025 m，道间距为12.5 m；记录长度9.2 s，采样间隔8 ms，所以每道1151个采样点；震源主频是40 Hz，子波延迟时间48 ms.图 12a-12c分别展示了BP2007模型第1炮、第821炮和第1641炮地震记录，对应激发点坐标分别是x=0 m、x=41.0 km和x=82.0 km，从中可以看出炮记录信噪比高，反射波等接收信息丰富.
|
图 12 地震炮记录.激发点坐标分别是(a) x=0 km；(b) x=41.0 km；(c) x=82.0 km Fig. 12 Common-shot gathers used for RTM at shot location (a) x=0 km, (b) x=41.0 km and (c) x=82.0 km, respectively |
图 12为BP2007模型Low-rank纯qP波逆时偏移成像剖面，图中可见：模型中的主要反射体都归位到正确的空间位置，两套高速岩体得到很好的刻画，岩体轮廓及侧翼成像清晰，高陡断层的边界展布清楚分明，由浅层至深层的各反射同相轴连续均衡，该剖面总体上信噪比高，成像效果较好.上述成像特点表明：本方法计算稳定且不存在伪横波干扰，对复杂模型成像具有较好的适应性，在改善深部构造的成像质量方面也具有一定的优势，此外，各向异性对构造成像具有较大的影响，针对强各向异性区域，特别是倾角急剧变化的区域，通过考虑各向异性偏移成像算法，有利于显著改善复杂构造的成像质量.
|
图 13 基于Low-rank分解的纯qP波叠前逆时偏移成像结果 Fig. 13 Prestack pure-qP wave RTM image of the 2D BP2007 TTI model using Low-rank decomposition algorithm |
有限差分一直是地震波正演模拟与逆时偏移成像领域最为流行的数值计算方法.但为了更好地处理TTI介质纯qP波正演及成像问题，本文首先构建了一种适用于任意TI介质的纯qP波传播算子，然后借助Low-rank分解求取该算子中的空间-波数域矩阵，同时引入Cerjan衰减边界条件来压制边界反射干扰，最终发展了一种间接的纯qP波正演模拟方案，并将其成功应用于复杂TI介质逆时偏移成像中.通过开展三组典型模型正演与复杂TTI介质成像测试，得出如下几点认识.
① TI介质Low-rank纯qP波正演模拟方法能够较为准确且稳定地刻画qP波运动学传播规律，既避免了纯qP波方程复杂的数学推导，又克服了耦合型波动方程对模型参数的限制，同时，还能适用于较大时间或空间步长与高频震源正演模拟；
②本文方法能够彻底消除残余伪横波干扰及数值频散假象，从而避免该类噪音对成像结果的污染，最终提高成像剖面的信噪比，改善成像质量；
③波场延拓的精度可控，这是由于地震波计算精度受控于传播矩阵的分解精度，而矩阵分解好坏决定于容许误差及时间步长，利用该项特点，即可建立相应的评价机制来监控地震波正演模拟精度，这有利于正演及成像软件编制和技术实用化；
④各向异性逆时偏移是解决复杂地质体准确成像的有力工具，本文成像方法对地下复杂构造成像拥有良好的适应性，在改善深部构造的成像质量方面也具有一定的优势；
⑤该类Low-rank正演与成像处理模式能够直接扩展到其他各类介质中，且只需利用相应介质中对应波型的频散方程来修改传播矩阵，即可实现该类波型的正、反向波场延拓，进而实现正演模拟或逆时偏移.
尽管基于Low-rank分解的各向异性逆时偏移成像方法能够有效提高TTI介质复杂构造的成像质量，获得较为理想的成像结果，但是该方法也存在明显的不足之处.
①计算效率较低，相对于TTI介质高斯束偏移，其计算时间增加了约9.8倍；此外，当模型更为复杂时，Low-rank分解的秩N会明显增加，必然导致计算成本随之成倍增加.
②采用Cerjan衰减边界条件很难完全消除边界反射干扰，势必造成传播波场中深部有效信号被边界反射污染，甚至淹没，进而使正演精度变低，成像质量变差，不利于高精度正演与成像处理.
③本方法主要解决构造成像问题，离保幅成像还有一定的差距，因为本文所构建的传播矩阵主要利用地震波相位信息，暂时还未能考虑密度变化及多波多分量的情况，因此地震波振幅存在一定程度的失真.
为此，在今后的研究过程中，首先应当引入完全匹配层(PML)边界处理技术、Low-rank有限差分或Low-rank傅里叶有限差分以及CPU/GPU协同并行加速等优化措施，来进一步提升计算精度与计算效率；其次，应重点考虑如何将优化的Low-rank纯qP波延拓算法应用于多分量成像、层析速度反演、最小二乘偏移、全波形反演以及实际资料处理中；总之，发展一种兼具效率与精度的Low-rank纯qP波速度分析与偏移成像技术系列是下一步的主要研究方向.
致谢作者首先感谢两位匿名审稿专家提出宝贵意见；同时，感谢美国石油公司(Amoco)及BP公司提供二维TTI Overthrust模型和二维TTI BP2007模型及对应炮数据；最后感谢中国石油大学(华东)研究生创新工程项目资助.
| 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 | |
| Alkhalifah T, Ma X X, bin Waheed U, et al. 2013. Efficient anisotropic wavefield extrapolation using effective isotropic models.//13th International Congress of the Brazilian Geophysical Society & EXPOGEF. Rio de Janeiro, Brazil, 1540-1543. | |
| Chen J B., Liu H. 2004. Optimization approximation with separable variables for the one-way wave operator. Geophysical Research Letters, 31(6): L06613. | |
| Cheng J B, Fomel S. 2014. Fast algorithms for elastic-wave-mode separation and vector decomposition using low-rank approximation for anisotropic media. Geophysics, 79(4): C97-C110. DOI:10.1190/geo2014-0032.1 | |
| 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 | |
| Crawley S, Brandsberg-Dahl S, McClean J, et al. 2010. TTI reverse time migration using the pseudo-analytic method. The Leading Edge, 29(11): 1378-1384. DOI:10.1190/1.3517310 | |
| Du X, Bancroft J C, Lines L R. 2005. Reverse-time migration for tilted TI media.//SEG Technical Program Expanded Abstracts. SEG, 1930-1934. | |
| Du X, Fowler P J, Fletcher R P. 2014. Recursive integral time-extrapolation methods for waves:A comparative review. Geophysics, 79(1): T9-T26. DOI:10.1190/geo2013-0115.1 | |
| Duveneck E, Milcik P, Bakker P M, et al. 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration.//SEG Technical Program Expanded Abstracts. SEG, 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 | |
| Etgen J. 1986. High order finite-difference reverse time migration with the two way nonreflecting wave equation. Stanford Exploration Project report 48, 133-146. | |
| Etgen J. 1989. Accurate wave equation modeling. Stanford Exploration Project report 60, 131-147. | |
| Etgen J T, Brandsberg-Dahl S. 2009. The pseudo-analytical method:Application of pseudo-Laplacians to acoustic and acoustic anisotropic wave propagation.//SEG Technical Program Expanded Abstracts. SEG, 2552-2556. | |
| Fang G, Du Q Z, Fomel S. 2013. Seismic wave extrapolation on a staggered grid using low-rank decomposition and low-rank finite differences.//SEG Technical Program Expanded Abstracts. SEG, 3433-3438. | |
| Fang G, Fomel S, Du Q Z, et al. 2014. Lowrank seismic-wave extrapolation on a staggered grid. Geophysics, 79(3): T157-T168. DOI:10.1190/geo2013-0290.1 | |
| Fang G, Fomel S, Du Q Z. 2014. Lowrank finite difference on a staggered grid and its application on reverse time migration. Journal of China University of Petroleum (in Chinese), 38(2): 44-51. | |
| Fang G, Hu J W, Fomel S. 2015. Weighted least square based lowrank finite difference for seismic wave extrapolation.//SEG Technical Program Expanded Abstracts. SEG, 3554-3559. | |
| 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 | |
| Fomel S. 2004. On anelliptic approximations for qP velocities in VTI media. Geophysical Prospecting, 52(3): 247-259. DOI:10.1111/gpr.2004.52.issue-3 | |
| Fomel S, Ying L X, Song X L. 2010. Seismic wave extrapolation using lowrank symbol approximation.//SEG Technical Program Expanded Abstracts. SEG, 3092-3096. | |
| Fomel S, Ying L X, Song X L. 2013. Seismic wave extrapolation using lowrank symbol approximation. Geophysical Prospecting, 61(3): 526-536. DOI:10.1111/gpr.2013.61.issue-3 | |
| Fowler P J. 2003. Practical VTI approximations:a systematic anatomy. Journal of Applied Geophysics, 54(3-4): 347-367. DOI:10.1016/j.jappgeo.2002.12.002 | |
| Huang J P, Huang J Q, Li Z C, et al. 2016. Staggered grid pseudo-spectral numerical simulation of first-order quasi P-wave equation for TTI media. OGP (in Chinese), 51(3): 487-496. | |
| Klié H, Toro W. 2001. A new acoustic wave equation for modeling in anisotropic media.//SEG Technical Program Expanded Abstracts. SEG, 1171-1174. | |
| Li Z C, Yang F S, Wang X D. 2016a. High-precision numerical simulation of first-order qP-waves in the transversely isotropic (TI) medium optimized by the LS-RSGFD method. Chinese J. Geophys. (in Chinese), 59(4): 1477-1490. DOI:10.6038/cjg20160428 | |
| Li Z C, Yong P, Huang J P, et al. 2016b. Elastic wave reverse time migration based on vector wavefield seperation. Journal of China University of Petroleum (in Chinese), 40(1): 42-48. | |
| Liu F Q, Morton S A, Jiang S S, et al. 2009. Decoupled wave equations for P and SV waves in an acoustic VTI media.//SEG Technical Program Expanded Abstracts. SEG, 2844-2848. | |
| Liu W Q, Wang X W, Wang Y C, et al. 2016. A regularized qP-wave equation for TTI media and its application to reverse time migration. Chinese J. Geophys. (in Chinese), 59(3): 1059-1069. DOI:10.6038/cjg20160326 | |
| Pestana R C, Stoffa P L. 2010. Time evolution of the wave equation using rapid expansion method. Geophysics, 75(4): T121-T131. DOI:10.1190/1.3449091 | |
| Pestana R C, Ursin B, Stoffa P L. 2011. Separate P-and SV-wave equations for VTI media.//SEG Technical Program Expanded Abstracts. SEG, 163-167. | |
| Sheng X, Zhou H B. 2014a. Accurate simulations of pure quasi-P-waves in complex anisotropic media. Geophysics, 79(6): T341-T348. DOI:10.1190/geo2014-0242.1 | |
| Sheng X, Zhou H B. 2014b. Pure quasi-P wave calculation in anisotropic media.//SEG Technical Program Expanded Abstracts. SEG, 3887-3891. | |
| Shi Y, Wang Y. 2016. Reverse time migration of 3D vertical seismic profile data. Geophysics, 81: S31-S38. DOI:10.1190/geo2015-0277.1 | |
| Song J. 2001. Optimal representation of multivariate functions or data in visualizable low-dimensional spaces. Chinese Science Bulletin, 46(16): 1337-1345. DOI:10.1007/BF03183384 | |
| Song X L, Fomel S, Ying L X, et al. 2011. Lowrank finite-differences for wave extrapolation.//SEG Technical Program Expanded Abstracts. SEG, 3372-3376. | |
| Song X L, Alkhalifah T. 2013. Modeling of pseudoacoustic P-waves in orthorhombic media with a low-rank approximation. Geophysics, 78(4): C33-C40. DOI:10.1190/geo2012-0144.1 | |
| Song X L, Fomel S, Ying L X. 2013. Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation. Geophysical Journal International, 193(2): 960-969. DOI:10.1093/gji/ggt017 | |
| Soubaras R, Zhang Y. 2008. Two-step explicit marching method for reverse time migration.//7SEG SEG Technical Program Expanded Abstracts. SEG, 2272-2276. | |
| Sun J Z, Zhu T Y, Fomel S. 2015. Viscoacoustic modeling and imaging using low-rank approximation. Geophysics, 80(5): A103-A108. DOI:10.1190/geo2015-0083.1 | |
| Tal-Ezer H. 1986. Spectral methods in time for hyperbolic equations. SIAM J. Numer. Anal., 23(1): 11-26. DOI:10.1137/0723002 | |
| Thomsen L. 1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966. DOI:10.1190/1.1442051 | |
| Tsvankin I. P-wave signatures and notation for transversely isotropic media:An overview.1996: 467-483. | |
| Wang L C, Chang X, Wang Y B. 2016. Forward modeling of pseudo P waves in TTI medium using staggered grid. Chinese J. Geophys. (in Chinese), 59(3): 1046-1058. DOI:10.6038/cjg20160325 | |
| Wang Y G, Xing W J, Xie W X, et al. 2007. Study of absorbing boundary condition by perfectly matched layer. Journal of China University of Petroleum (in Chinese), 31(1): 19-24. | |
| Wu D, Yao G, Cao J, et al. 2016. Least-squares RTM with L1 norm regularization. J. Geophys. Eng., 13: 666-673. DOI:10.1088/1742-2132/13/5/666 | |
| Wu Z D, Alkhalifah T. 2014. The optimized expansion based low-rank method for wavefield extrapolation. Geophysics, 79(2): T51-T60. DOI:10.1190/geo2013-0174.1 | |
| Xu W C, Yang G Q, Li Z C, et al. 2016. Pseudo acoustic equation for TI medium attenuation based on the GSLS model. Chinese J. Geophys. (in Chinese), 59(6): 2232-2244. DOI:10.6038/cjg20160626 | |
| Yan J, Liu H. 2016. Modeling of pure acoustic wave in tilted transversely isotropic media using optimized pseudo-differential operators. Geophysics, 81(3): T91-T106. DOI:10.1190/geo2015-0111.1 | |
| Ying X Y, Liu Q. 2016. Anisotropic rock physics modeling of tight sandstone and applications. Journal of China University of Petroleum (in Chinese), 40(2): 52-58. | |
| Yong P, Huang J, Li Z, et al. 2016. Elastic-wave reverse-time migration based on decoupled elastic-wave equations and inner-product imaging condition. J. Geophys. Eng., 13: 953-963. DOI:10.1088/1742-2132/13/6/953 | |
| Yoon K, Suh S, Ji J, et al. 2010. Stability and speedup issues in TTI RTM implementation.//SEG Technical Program Expanded Abstracts. SEG, 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 pseudospectral/finite-difference scheme for solving TTI pure P-wave equation.//13th International Congress of the Brazilian Geophysical Society & EXPOGEF. Rio de Janeiro, Brazil, 1322-1327. | |
| Zhang L B, Rector J W Ⅲ, Hoversten G M. 2005. Finite-difference modelling of wave propagation in acoustic tilted TI media. Geophysical Prospecting, 53(6): 843-852. DOI:10.1111/gpr.2005.53.issue-6 | |
| Zhang Y, Zhang H. 2009a. A stable TTI reverse time migration and its implementation.//SEG Technical Program Expanded Abstracts. SEG, 2794-2798. | |
| Zhang Y, Zhang G Q. 2009b. One-step extrapolation method for reverse time migration. Geophysics, 74(4): A29-A33. DOI:10.1190/1.3123476 | |
| 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 | |
| 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.//SEG Technical Program Expanded Abstracts. SEG, 194-198. | |
| 程玖兵, 康玮, 王腾飞. 2013. 各向异性介质qP波传播描述Ⅰ:伪纯模式波动方程. 地球物理学报, 56(10): 3474–3486. DOI:10.6038/cjg20131022 | |
| 方刚, FomelS, 杜启振. 2014. 交错网格Lowrank有限差分及其在逆时偏移中的应用. 中国石油大学学报(自然科学版), 38(2): 44–51. | |
| 黄建平, 黄金强, 李振春, 等. 2016. TTI介质一阶qP波方程交错网格伪谱法正演模拟. 石油地球物理勘探, 51(3): 487–496. | |
| 李振春, 杨富森, 王小丹. 2016a. 基于LS-RSGFD方法优化的横向各向同性(TI)介质一阶qP波高精度数值模拟. 地球物理学报, 59(4): 1477–1490. DOI:10.6038/cjg20160428 | |
| 李振春, 雍鹏, 黄建平, 等. 2016b. 基于矢量波场分离弹性波逆时偏移成像. 中国石油大学学报(自然科学版), 40(1): 42–48. | |
| 刘文卿, 王西文, 王宇超, 等. 2016. 带正则化校正的TTI介质qP波方程及其逆时偏移方法和应用. 地球物理学报, 59(3): 1059–1069. DOI:10.6038/cjg20160326 | |
| 王璐琛, 常旭, 王一博. 2016. TTI介质的交错网格伪P波正演方法. 地球物理学报, 59(3): 1046–1058. DOI:10.6038/cjg20160325 | |
| 王永刚, 邢文军, 谢万学, 等. 2007. 完全匹配层吸收边界条件的研究. 中国石油大学学报(自然科学版), 31(1): 19–24. | |
| 徐文才, 杨国权, 李振春, 等. 2016. 基于GSLS模型TI介质衰减拟声波方程. 地球物理学报, 59(6): 2232–2244. DOI:10.6038/cjg20160626 | |
| 印兴耀, 刘倩. 2016. 致密储层各向异性地震岩石物理建模及应用. 中国石油大学学报(自然科学版), 40(2): 52–58. | |