广角地震反射特征及反演研究

2017-11-01武泗海

石油地球物理勘探 2017年5期

关键词：平面波广角波场

武泗海赵虎尹成贾鹏

(①西南石油大学地球科学与技术学院，四川成都 610500； ②西南石油大学天然气地质四川省重点实验室，四川成都 610500； ③川庆钻探地球物理勘探公司，四川成都 610213)

·综合研究·

广角地震反射特征及反演研究

武泗海*①②赵虎①②尹成①②贾鹏③

武泗海，赵虎，尹成，贾鹏﹒广角地震反射特征及反演研究.石油地球物理勘探，2017,52(5):1005-1015.

在小炮检距情况下，地震反射波振幅随着入射角的变化而变化，但相位基本保持不变。当入射角接近或大于临界角时，地震反射波相位也随之发生变化。振幅与相位都是波形信息的重要特征，因此，本文对广角反射波特征进行了研究。首先以平层介质为例，结合广角地震记录合成方法，研究了t-x域、τ-p域的反射波振幅/相位随入射角变化的特征，并希望据此进行广角的振幅/相位反演。但由于τ-p变换受采样间隔等因素的影响，本文又利用逆时偏移成像提取共成像点道集的方法，获取广角反射信息，并将其应用于弹性介质参数的反演试验。最后以平层介质和局部Marmousi-2模型为例做了相关实验，反演结果也验证了广角振幅与相位及联合反演的有效性。

广角反演τ-p变换共成像点道集逆时偏移

1 引言

在实际勘探工区中，如逆掩断层区、海洋深水区、火山岩、碳酸盐岩上覆层发育区，采集的资料经常遇到超临界角反射的问题。当入射角大于临界角时，由于反射系数变为复数产生附加相角，反射波相位产生变化，波形产生畸变。王志等[1]通过高斯射线正演模拟广角反射的方法，讨论了广角反射特征及影响波形识别的主要因素；孙成禹等[2]通过分析广角反射数据的特征，研究广角数据的校正方法及存在的问题；孟宪军等[3]通过研究地震反射系数相角变化特征，指出了临界角前后的相角变化特征。

AVO/AVA技术根据振幅随炮检距变化的特征进行地下介质的岩性信息识别[4,5]，但这仅仅考虑振幅大小及正负极性对反演的影响，实际上，反射波相位的变化也同样重要。对于子波而言，振幅和相位对波形都具有非常重要的影响。此外，由于振幅在传播过程中存在透射、波前扩散及反射临界角等问题，振幅信息的提取精度也影响了反演效果，而相位的影响较小，其精度相对较高。因此，进一步研究相位变化特征对广角反射资料的应用具有一定实际意义。

在广角反射存在的情况下，由于Zoeppritz方程是基于平面波入射的假设，这就导致平面波反射系数(PRCs)与实际地震反射特征并不完全一致。针对这一问题，Haase[6]和Ursenbach等[7]研究了在大炮检距情况下的球面波AVO响应，并提出了一种简洁的平缓界面球面波反射系数(SRCs)的计算方法；Ayzenberg等[8]提出了有效反射系数(ERCs)的概念，将其应用于弯曲界面的AVO分析。另外，前人就球面波场的平面波分解问题进行了研究，使Zoeppritz方程在广角范围依然适用，van der Baan等[9]研究了球面波的τ-p域分析方法，使平面波反射系数继续适用于振幅随慢度变化的反演。Zhu等[10]研究了广角反射在τ-p域的响应，通过反射率法模拟结果证实τ-p域反射特征与平面波反射系数基本一致。但τ-p变换也存在不足之处，为了避免假频，对空间采样提出了严格的要求。

逆时偏移(RTM)作为目前成像精度最高的偏移技术，在提高复杂介质的成像质量方面表现出优越的性能，而且RTM输出的共成像点道集还能提供振幅和相位信息，为后续属性解释提供依据[11-13]。目前，RTM技术已从二维到三维，从声波方程逐步发展到弹性波方程。许多学者针对RTM中的成像条件、边界条件和成像噪声压制方面的问题做了大量的研究工作。薛东川[14]对比了两大类成像条件：入射波到达时间类与入射波/反射波的相关类成像条件的特点、计算量及存储量；张智等[15]提出了弹性波逆时偏移的稳定激发振幅成像条件，与归一化成像条件相比，该成像条件具有更好的成像能力和空间分辨率；陈康等[16]提出了频率域修正的四阶拉普拉斯滤波方法，该方法在改善成像噪声压制效果的同时，还保持了子波的振幅与相位。此外，为解决RTM波场存储和计算量大的问题，王保利等[17]分析了目前常用的几类存储方案，并结合存储边界策略和检查点技术，提出了有效存储边界策略，在保证波场延拓精度的同时，大幅降低了RTM的存储需求；王维红等[18]提出了针对VSP RTM的PML边界存储方法，有效地节约了存储空间；李博等[19]利用CUDA编程架构实现CPU/GPU异构并行技术，解决了逆时偏移的计算量大的瓶颈，并给出了偏移算法的实施策略。

本文首先在研究广角反射波波形变化特征的基础上，利用τ-p变换进行平面波分解后提取振幅与相位信息，对比分析t-x域、τ-p域的反射波振幅/相位随入射角变化的特征。由于τ-p变换的限制，在此基础上，本文利用RTM提取角度域共成像点道集，并以平层和局部Marmousi-2模型为例说明和实现了广角AVA/PVA反演，反演结果也进一步验证了广角反演的有效性。

2 广角反射波变化特征研究

Zoeppritz方程精确地描述了平面波在弹性分界面的反射及透射变化特征，由于方程繁杂，其中不仅包含反射、透射变化，还存在横波及转换波等不同波型的变化。因此，本文以PP波反射情况为研究对象，为方便起见，采用以下PP波反射系数公式[20]

Rpp=

(1)

其中

A=(r4-Q)2T1T3

B=(r4-Q-1)2T0T1T2T3

由于广角反射波仅出现在临界角之后，根据Castagna等[21]提出的四类AVO响应，本文从中选取Ⅰ类和Ⅲ类为例分析广角反射波的变化特征,其对应模型的弹性参数如表1所示。

表1 Ⅰ类和Ⅲ类AVO模型的弹性参数

根据式(1)计算上述两类模型的反射系数—振幅及相位随入射角变化曲线，如图1所示，通过对比可以得出以下结论。

(1)Ⅰ类模型的AVO响应曲线：在小入射角情况下，振幅随着入射角增大而变化，相位基本保持不变。当入射角大于第一临界角时，振幅存在突变；相位则随着入射角增加而变为负值; 在第二临界角位置，由于相位角度范围限制在-180°～180°,相位产生突变。

(2)Ⅲ类模型的AVO相应曲线：由于该模型中不存在临界角，相位基本保持不变，只有振幅随着入射角的增加而变化。

为了直观地分析广角地震反射特征，下面给出根据Ⅰ类模型生成的广角反射记录，如图2所示。

图1 Ⅰ类和Ⅲ类AVO模型的反射振幅(a)与相位(b)随入射角变化曲线

图2 Ⅰ类AVO模型的合成广角反射记录

从图中可以看出：广角记录所表现的特征与平面波反射系数特征一致，广角的反射振幅比临界前强，且存在明显的相位变化特征。

3 广角反射的τ -p域分析

在小入射角情况下，平面波反射系数与球面波反射系数差异较小；但当入射角大于临界角时，平面波反射波系数(PRC)不再适用于远炮检距t-x域球面波记录。正如前文所述，通过将平面波分解，可以为广角信息在τ-p域利用平面波反射系数进行反演提供可能。本文通过波动方程正演模拟得到Ⅰ类模型的平层介质共炮点t-x域记录(图3)，并利用τ-p变换实现了平面波分解，得到对应的τ-p域记录(图4)。

通过提取τ-p域不同入射角的反射波记录，将提取波形按角度排列得到如图5所示的角道集，并结合前文合成的理论角道集，可以看出，两者无论是临界角所处的位置，还是振幅和相位变化特征都基本吻合。

与t-x域相比，τ-p域广角反射振幅与相位更接近于平面波反射系数。图5中道集的振幅与相位随入射角变化的特征曲线分别如图6、图7所示。由于在计算反正切值时的角度范围限制，导致了图1中相位在第二临界角附近产生突变。为了降低该突变对反演结果的影响，图7中进行了相位展开处理，使得相位变化更加连续，并给出了与t-x域提取振幅与相位及理论反射波振幅和相位曲线对比。

图3 Ⅰ类AVO模型共炮点记录的垂向速度分量(局部)

图4 Ⅰ类AVO模型共炮点记录的τ -p域波场

图5 τ -p域选排后的PP波广角反射记录(0°～70°)

图6 Ⅰ类AVO模型的广角振幅随入射角变化曲线

图7 Ⅰ类AVO模型的广角相位随入射角变化曲线

从图中不难发现：在入射角较小时，t-x域反射波振幅和相位与理论曲线基本一致，但随着入射角的增大，提取的振幅和相位曲线与平面波理论反射曲线存在明显误差。与之相反，τ-p域提取的振幅和相位信息在广角则能保持较高的精度。由于大角度的反射的振幅更强，相位变化特征也更加明显，这在一定程度上也有利于提高介质参数反演的稳定性和精度。

4 RTM角道集AVA/PVA反演

通过广角反射的τ-p域分析，在大入射角情况下，τ-p域的反射特征与理论结果基本一致。但τ-p变换要求规则的空间采样间隔和横向为均匀介质[10]，这与实际情况并不完全符合。因此，本文利用RTM方法获取角度域共成像点道集(ADCIGs)，并将其应用于广角AVA/PVA反演。

由于RTM是基于双程波波动方程，与射线类和单程波类成像方法相比，RTM成像结果更为精确。因此，RTM是目前主流的提取ADCIGs的方法，基于RTM提取ADCIGs主要分为基于方向矢量类方法、局部平面波分解类方法和局部移动成像条件法的三大类[22]。

为了获取可靠的角度域共成像点道集，并将其应用于广角AVA/PVA反演，本文采用Zhang等[23]提出的ADCIGs提取方法。ADCIGs提取方法的主要难点在于共成像点位置的局部入射角的计算，即波场传播方向角度与反射界面角度的法向之差。

4.1 局部入射角的计算

由于平面P波的位移、速度及加速度的极化方向平行于传播方向，因此波场传播方向可以通过P波波场的梯度得到，对于弹性波波动方程，P波波场可以通过对波场求散度进行分离，即

(2)

式中：U表示地震波场的位移；λ和μ表示弹性参数；ρ表示密度。

在二维情况下，其传播方向余弦(γ1,γ2)可表示为

(3)

(4)

一般情况下，波场传播方向角度α为

(5)

对于反射界面法向，可以将其RTM成像结果进行Hilbert变换，再利用变换结果的瞬时相位的空间梯度(瞬时波数)得到[24]。假定s(x,z)为偏移图像，与其Hilbert变换q(x,z)构成的复信号c(x,z)可表示为

c(x,z)=s(x,z)+iq(x,z)

(6)

对于复数图像信息c(x,z)，对其瞬时相位φ求取空间梯度可以得到瞬时波数，空间梯度计算公式为

(7)

式中w1和w2分别表示反射体瞬时波数的水平和垂直分量。由此可估算反射体的法向角度β。在二维情况下，反射体法向角度β为

(8)

通过式(5)和式(8)可求得波场传播角度α及反射体法向角度β。如图8所示，在二维情况下，MM′表示反射界面，则成像点处的入射角θ为

θ=α-β

(9)

图8 局部入射角计算示意图

4.2 RTM角道集的提取

为了获取较精确的反射波振幅与相位信息,成像条件的选择至关重要。根据Sandip等[25]对成像条件的保幅性研究，振幅比成像条件具有更高的成像分辨率，能够真实地反映反射系数的大小。因此，采用振幅比成像条件，其成像时刻为成像点的激发时刻。

利用弹性波一阶应力—速度方程式的高阶差分计算进行震源波场的正推和检波点波场的逆时延拓。在震源/检波点的波场延拓过程中，利用散度算子提取弹性波场中P波分量用于成像[26]，与常规的水平/垂直分量成像相比，其成像振幅更加精确。在每一个共炮点记录的波场反传的同时，一方面根据成像时刻及振幅比成像条件计算成像点的成像值；另一方面，需要保存当前炮在该成像点成像时刻的传播角度α,为下一步提取角道集做准备。

需要注意的是，在成像过程中，几何效应导致波场在垂向上存在拉伸畸变[27,28]，严重降低了广角信息的分辨率和精度。因此在进行成像之前首先对子波进行了拉伸校正，经过校正后的广角信息有了明显改善。最后，通过对所有炮的偏移结果进行叠加，得到最终偏移成像结果。

图9 RTM角道集计算的流程

根据式(8)，结合最终成像剖面可以求得每个成像点位置的法向角度β，进而求得局部入射角θ。在每一炮内，根据入射角进行角度分选，提取成像值，按照入射角度进行选排，即可得到ADCIGs。ADCIGs提取的总体流程如图9所示。一般情况下，由于角道集的离散角度与真实反射角存在差异，因此需要对抽取的角道集进行平滑处理，其平滑公式为

(10)

式中：E表示抽取的角度域共成像点道集；kθ表示角度道中的离散入射角； Image表示成像值；θ表示当前时刻成像点的入射角；σ表示高斯函数的方差；Num_shot表示偏移总炮数。

4.3 AVA/PVA反演策略

本文利用模拟退火法的反演方法进行广角AVA/PVA反演试验工作，模拟退火法反演有不依赖于初始模型、不需要求解雅克比矩阵、且不易陷入局部最优解等优点，在地球物理资料反演中有许多成功实例。

根据式(1),PP波反射系数RPP可以表示为入射角θ及r1、r2、r3、r4的函数

RPP=RPP(θ,r1,r2,r3,r4)

(11)

由于r1、r2、r3分别表示利用vP1规范化后的vP2、vS1、vS2；r4是利用ρ1规范化后的ρ2，因此式(11)中只有入射角θ为独立变量。通常情况下，可通过层析成像获得上层介质纵波速度vP1，而ρ1则可根据加德纳(Gardner)公式计算，或从实验室测定得到。

当已知vP1和ρ1参数时，对其他四个参数进行反演，在L2范数意义下，目标函数φ(x)可以表示为

(12)

5 模型试算

5.1 层状介质模型

为了说明RTM角道集抽取和广角反射振幅与相位信息提取方法的有效性，以层状介质模型为例(图10)，模型尺寸为8km×2km，上/下介质纵、横波速度(km/s)及密度(g/cm3)分别为4.0、1.54、2.3和2.0、0.88、2.0。常规的逆时偏移包括震源波场的正传和检波点的逆时延拓两个波场模拟过程，采用高阶交错网格有限差分法实现波场的模拟，其中炮点与检波器均匀布设在水平地表。图11为炮点位于4km处的单炮记录的速度水平分量与垂直分量。在切除直达波后，该记录作为波场逆时延拓的边界条件。

图10 单一层状介质模型

图11 切除直达波后的单炮记录水平分量(a)及垂直分量(b)

通过前文所述方法得到层状介质的单炮偏移结果(图12a)及经过子波拉伸校正[29]的结果(图12b)。然后根据入射角信息，从多炮偏移结果中提取共成像点的不同角度反射信息获取ADCIGs(图13)。

图12 PP波偏移结果拉伸校正前(a)、后(b)对比

图13 水平位置在4km处的ADCIGs

由于广角子波波形存在畸变，如果提取峰值振幅作为反射系数，会存在较大的误差。此处采用求取子波包络峰值的方法，降低波形变化带来的误差。另外，相位则利用互相关法计算，通常垂直入射/反射的波形常常也并不是零相位，因此，通常以零度入射的相位作为参照子波进行相位的相关计算。通过对ADCIGs中提取的振幅信息和相位随入射角变化曲线(图14)可以看出：实际振幅在小角度与理论值基本吻合，大角度由于球面波效应，两者的差异逐渐增加；相位则基本保持不变。这与前文所述Ⅰ类AVO曲线特征基本一致。

图14 ADCIGs提取的振幅和相位信息与理论值对比

5.2 Marmousi-2模型

为了验证本文提出的广角AVA/PVA联合反演的有效性，以局部Marmousi-2模型(图15)为例,将模型的尺寸设为2.0km(宽)×0.88km(高)。由于相位变化仅存在于临界角后，因此将目的层上、下层介质的纵波速度(km/s)、横波速度(km/s)和密度(g/cm3)分别填充为2.0、0.88、2.0和4.0、1.54、2.3。逆时偏移的模拟参数分别为：震源为40Hz主频Ricker子波;采样间隔为0.4ms;空间采样间隔均为4m；记录时间为1400ms；炮数为60个;检波点个数为501个，均匀分布于地表。

在震源正向延拓过程中，根据每个成像点的激发时间，结合式(5)获得震源纵波传播方向角度。图16为震源位于1.2km处的P波传播方向角度。

在检波点逆时延拓过程中，根据成像点的激发时刻，应用振幅比成像条件进行成像，完成单炮偏移，通过叠加所有单炮结果获得最终的偏移剖面，再由式(7)获得其复图像信号相位得到瞬时波数，进而计算反射体法向角度(图17)。已知震源传播方向角度和反射体法向角度，即可计算在当前炮中该模型的入射角分布情况(图18)。

图15 Marmouis-2模型(红色框区为示例模型)

图16 水平距离1km处震源的P波传播方向角度

图17 局部Marmousi-2模型的反射体法向角度

在多炮偏移中，不同位置的炮点可以提供目标点的不同角度反射信息。在确定该角度后，结合式(10)可以从不同位置炮点的偏移结果中抽取角度域共成像点道集。图19所示分别为原始模型与修改模型的5组角道集，在原始模型中(图19a)，由于速度差异较小，目标深度上层反射层次清晰。为了凸显广角反射特征，修改后的模型(图19b)在目标深度速度差异较大，对上层角道集产生了一定的影响。

图18 局部Marmousi-2模型的入射角

从目标层位的角道集中可以提取广角反射波的振幅和相位信息，如图19所示。本节将其应用于中广角反演中，通常情况下，广角畸变的存在导致难以获取精确的大角度反射信息。本文假定在上层纵波速度vP1及密度ρ1已知的条件下，利用广角振幅和相位进行其他弹性参数的反演。

将从角道集中提取的广角振幅与相位信息带入式(12)进行迭代反演。本文将AVA、PVA及联合反演的结果进行对比，如图20所示。由于受到所提取ADCIGs质量的影响，本文仅采用入射角为0°～50°的振幅与相位数据。

图19 模型修改前(a)、后(b)的ADCISs，角度范围为-90°～90°

图20 模型的AVA/PVA反演结果

表2 三种反演结果相对误差对比

表2给出了上述三类方法的反演结果与模型真实参数的相对误差值。通过对比分析反演结果可知：反射波振幅主要取决于上、下界面的波阻抗比值，而AVA反演又依赖于振幅信息，因此vP2和ρ2的相对误差呈符号相反的特征，但是波阻抗的反演结果较为准确。由于振幅信息受透射及波前扩散影响较大，AVA反演结果并不十分理想；相比较而言，相位信息较为可靠，其对vS2反演的效果非常好。同时，研究还发现，由于相位对密度变化极不敏感，因而导致ρ2的反演误差较大；在AVA/PVA联合反演中，观测数据与反演结果拟合较好，通过表2也可以看出，联合反演的结果较两者都更接近于模型参数。

6 结论和认识

通过对广角地震反射波振幅/相位的分析和研究，得到以下结论和认识。

(1)在反射波临界角存在的情况下，随入射角的增大，反射波振幅随之变化的同时，相位也并不总是零相位，在大入射角时，甚至会出现极性反转的情况。

(2)地震波形信息包含了振幅与相位，特别是在广角条件下，相位特征更明显。在波形反演中隐含了反射波的振幅和相位信息，本文则以AVA/PVA反演的方式利用了广角反射的波形信息。

(3)通过τ-p域的分析结果可以看出，由于Zoeppritz方程以平面波入射为假设条件，其并不适用于广角地震反射波特征的分析。本文通过模型说明了广角τ-p域分析的可行性，但由于τ-p变换尚存不足之处，τ-p域的广角反演还有待于进一步研究。

(4)本文利用逆时偏移获取角度域共成像点道集，并将提取的广角振幅与相位信息用于AVA/PVA反演，模型实例验证了该方法的可行性。

(5)反演结果对比发现：AVA/PVA的反演结果各有所长，AVA对纵波速度/密度反演效果较好，PVA则对横波速度反演较好，联合反演在一定程度上弥补了两者的不足。由此可见，在多波多分量采集、大炮检距和超临界角采集中，广角反射的资料具有一定的应用价值。

基于逆时偏移提取的角道集质量直接影响了广角反演的结果，尤其是在复杂模型中，激发时刻的计算及基于振幅比成像条件也存在不稳定性，而且广角地震反射特征复杂，存在严重的拉伸效应；另外，折射首波很可能与一次反射波相混淆。这就导致了基于方向矢量类的反射角计算方法并不能完全适用。因此，逆时偏移方法中入射角的计算和基于相关成像条件下的拉伸校正、角道集提取和保幅性的研究是下一步的工作方向。

[1] 王志,贺振华,黄德济等.广角反射波场特征研究及正演模拟分析.地球物理学进展,2003,18(1):116-121. Wang Zhi,He Zhenhua,Huang Deji et al.The wave-field features researching of wide-angle reflection and the analyzing of forward modeling.Progress in Geophysics,2003,18(1):116-121.

[2] 孙成禹,倪长宽,李胜军等.广角地震反射数据特征及校正方法研究.石油地球物理勘探,2007,42(1)：24-29. Sun Chengyu,Ni Changkuan,Li Shengjun et al.Feature of wide-angle reflection data and correction me-thod.OGP,2007,42(1):24-29.

[3] 孟宪军,刘福平,王玉梅.地震反射系数的相角变化.石油地球物理勘探,2010,45(增刊1):121-124. Meng Xianjun,Liu Fuping,Wang Yumei.The phase angle change for seismic reflection coefficient.OGP,2010,45 (S1):121-124.

[4] 代荣获,张繁昌,刘汉卿.叠前地震资料AVA变参数反演方法.地球物理学进展,2015,30(1):261-266. Dai Ronghuo,Zhang Fanchang,Liu Hanqing.AVA inversion using replacing parameters method for pre-stack seismic data.Progress in Geophysics,2015,30(1):261-266.

[5] 代荣获,张繁昌,刘汉卿等.基于贝叶斯理论的逐次迭代非线性AVA反演方法.吉林大学学报(地球科学版),2014,44(6):2026-2033. Dai Ronghuo,Zhang Fanchang,Liu Hanqing et al.Non-linear pre-stack seismic AVA inversion based on Bayesian theory using successive iteration method.Journal of Jilin University(Earth Science Edition),2014,44(6):2026-2033.

[6] Haase A B.Plane waves,spherical waves and angle-dependent P-wave reflectivity in elastic VTI-models.SEG Technical Program Expanded Abstracts,2002,21:304-307.

[7] Ursenbach C P,Haase A B.Efficient spherical-wave AVO modeling.The Leading Edge,2007,26(12):1584-1589.

[8] Ayzenberg M,Tsvankin I,Aizenberg A et al.Effective reflection coefficients for curved interfaces in transversely isotropic media.Geophysics,2009,74(5):33-53.

[9] van der Baan M and Smit D.Amplitude analysis of isotropic P-wave reflections.Geophysics,2006,71(6):93-103.

[10] Zhu Xinfa,McMechan G.Analysis of phase variation with ray parameter in theτ-pdomain.SEG Technical Program Expanded Abstracts,2012,31:1-5.

[11] 王保利,高静怀,陈文超等.逆时偏移中用Poynting矢量高效地提取角道集.地球物理学报,2013,56(1):262-268. Wang Baoli,Gao Jinghuai,Chen Wenchao et al.Extracting efficiently angle gathers using Poynting vector during reverse time migration.Chinese Journal of Geophysics,2013,56(1):262-268.

[12] 杨午阳,张厚柱,撒利明等.逆时偏移关键问题探讨.石油地球物理勘探，2016,51(6):1251-1262. Yang Wuyang,Zhang Houzhu,Sa Liming et al.Some issues about reverse time migration.OGP,2016,51(6):1251-1262.

[13] 吴国忱,秦海旭.各向同性弹性波随机边界逆时偏移.地球物理学进展,2014,29(4):1815-1821. Wu Guochen,Qin Haixu.Elastic reverse time migration in isotropic medium based on random boundary.Progress in Geophysics,2014,29(4):1815-1821.

[14] 薛东川.几种叠前逆时偏移成像条件的比较.石油地球地球物理勘探,2013,48(2):222-227. Xue Dongchuan.Imaging condition comparison of pre-stack reverse time migration.OGP,2013,48(2):222-227.

[15] 张智,刘有山,徐涛等.弹性波逆时偏移中的稳定激发振幅成像条件.地球物理学报,2013，56(10):3523-3533. Zhang Zhi,Liu Youshan,Xu Tao et al.A stable excitation amplitude imaging condition for reverse time migration in elastic wave equation.Chinese Journal of Geophysics,2013,56(10):3523-3533.

[16] 陈康,吴国忱.逆时偏移拉普拉斯算子滤波改进算法.石油地球物理勘探,2012,47(2):249-255. Chen Kang,Wu Guochen.An improved filter algorithm based on Laplace operator in reverse-time migration.OGP,2012,47(2):249-255.

[17] 王保利,高静怀,陈文超等.地震叠前逆时偏移的有效边界存储策略.地球物理学报,2012,55(7):2412-2421. Wang Baoli,Gao Jinghuai,Chen Wenchao et al.Efficient boundary storage strategies for seismic reverse time migration.Chinese Journal of Geophysics,2012,55(7):2412-2421.

[18] 王维红,郭雪豹,石颖等．VSP逆时偏移及其存储策略研究．石油地球物理勘探,2015,50(1)：77-83. Wang Weihong,Guo Xuebao,Shi Ying et al．VSP reverse time migration and its data storage strategy．OGP,2015,50(1)：77-83．

[19] 李博,刘红伟,刘国峰等.地震叠前逆时偏移算法的CPU/GPU实施对策.地球物理学报,2010,53(12):2938-2943. Li Bo,Liu Hongwei,Liu Guofeng et al.Computational strategy of seismic reverse time migration on CPU/GPU.Chinese Journal of Geophysics,2010,53(12):2938-2943.

[21] Castagna J P and Swan H W.Principles of AVO cross-plotting.The Leading Edge,1997,16(4):337-344.

[22] Jin H,McMechan G A and Guan H M.Comparison of methods for extracting ADCIGs from RTM.Geophysics,2014,79(3):89-103.

[23] Zhang Q S and McMechan G A.Direct vector-field method to obtain angle-domain common-image gathers from isotropic acoustic and elastic reverse time migration.Geophysics,2011,76(5):135-149.

[24] Havlicek J P,Havlicek J W,Mamuya N D et al.Skewed 2D Hilbert transforms and computed AM-FM models.IEEE International Conference on Image Proce-ssing,1998,602-606.

[25] Chattopadhyay S,McMechan G A.Imaging conditions for pre-stack reverse-time migration.Geophysics,2008,73(3):81-89.

[26] 李振春,张华,刘庆敏等.弹性波交错网格高阶有限差分法波场分离数值模拟.石油地球物理勘探,2007,42(5):249-255. Li Zhenchun,Zhang Hua,Liu Qingmin et al.Numeric simulation of elastic wavefield separation by staggering grid high-order finite difference algorithm.OGP,2007,42(5):249-255.

[27] Tygel M,Schleicher J and Hubral P.Pulse distortion in depth migration.Geophysics,1994,59(10):1561-1569.

[28] Zhu Xinfa and McMechan G A.Stretch-free migration imaging condition.SEG Technical Program Expanded Abstracts,2013,32:3882-3887.

[29] 杨佳佳，何兵寿，陈婷．逆时深度偏移中的子波拉伸校正．石油地球物理勘探，2016，51(1)：135-140. Yang Jiajia，He Bingshou，Chen Ting．Wavelet stretch correction for reverse-time migration．OGP，2016，51(1)：135-140．

(本文编辑：刘英)