VTI介质角度域叠前深度偏移

李 江 李庆春

(①中国煤炭科工集团西安研究院有限公司，陕西西安 710077； ②长安大学地质工程与测绘学院，陕西西安 710054)

地球介质一般具有各向异性特征，其中VTI介质(具有垂直对称轴的横向各向同性介质)是一种常见的各向异性介质。在地震勘探中一般将各向异性视为速度各向异性，在研究区域构造时是可行的。随着勘探目标体尺度的缩小，忽略地层各向异性影响的地震偏移成像导致目标区空间反射波归位不准确，绕射波不能完全收敛，甚至产生构造扭曲或构造假象，严重降低了地震成像精度[1-3]。近年来，人们研究了VTI介质的偏移方法。例如：Alkhalifah[4]实现了VTI介质的相移偏移方法，但该方法无法适应介质速度的横向变化；吴国忱等[5]实现了VTI介质广义屏叠前深度偏移方法，但其中对各向异性频散关系的近似使该方法仅局限于弱各向异性介质；韩建光等[6-7]提出了VTI介质高斯束叠前深度偏移方法，数值实验表明，该方法不受各向异性程度的限制，能适应强各向异性介质。虽然高斯束偏移方法具有高效、灵活的优点，但成像效果明显受速度横向变化剧烈的影响，因此人们利用波动方程叠前深度偏移对复杂区精确成像。刘礼农等[8-9]发展了VTI介质波动方程最优分裂傅里叶偏移方法；程玖兵等[10-11]研究了qP波在各向异性介质中的传播模式； Sun等[12]等发展了VTI介质逆时偏移方法；Zhang 等[13]提出了一种三维TTI介质纯P波方程求解方法； Hua 等[14]、黄建平等[15]、朱峰等[16]发展了VTI介质和黏声VTI介质的有限差分叠前深度偏移方法。但以上波动方程类偏移方法大多是在炮域实现叠前偏移，研究表明，在横向变速剧烈情况下，共炮偏移生成的炮域共成像点道集在运动学和动力学上存在假象。为此， Hoop等[17]、Xu等[18]研究了角度域共成像点道集对射线多路径问题的适应性，提出通过角道集的叠加成像可以避免构造假象，有效提高偏移剖面的成像效果。Stork等[19]进一步证明，基于波动方程偏移提取的角度域共成像点道集(ADCIG)不存在假象。Xie等[20]通过Gabor-Daubechies框架小波束域波场分解进行叠前角度域偏移； James等[21]通过波动方程偏移提取炮检距域共成像点道集和ADCIG； Save等[22]提出一种采用波动方程偏移提取ADCIG的方法，将偏移过程和ADCIG计算过程分离，即在成像空间通过傅里叶域的径向道变换计算出以入射角为参数的角道集，并对众多提取ADCIG的方法进行了分类和总结。Biondi[23]讨论了各向异性偏移的ADCIG；李振春等[24]、张敏等[25]对比、分析了单平方根和双平方根偏移提取ADCIG的方法，并发展了保幅的角度域偏移成像方法。Tang等[26]发展了角度域平面波逆时偏移方法。由于ADCIG是目前唯一不存在假象的道集，因此在地震数据处理中具有优势[27]。

本文研究了VTI介质角度域偏移方法，以各向同性双平方根方程的角度域偏移方法为基础，从VTI介质qP波频散关系出发，推导出VTI介质角度域偏移的波场延拓算子，在频率—波数域处理以横向均匀速度传播的波场，在空间域处理具有速度扰动特征的波场以提高波场延拓精度。断层模型和Hess VTI模型测试结果表明：该方法不受介质各向异性程度的影响，可使地震波场准确归位、绕射波收敛，从而对VTI介质精确成像；与各向同性偏移方法相比，该方法消除了各向异性效应对波场传播的影响，明显提高了成像精度；偏移生成的ADCIG避免了各向同性偏移方法使成像道集无法校平的缺陷，且具有更丰富的角度信息。实际资料的应用效果证实了本文方法的实用性。

VTI介质描述了由周期性的薄互层和平行排列的微裂隙引起的各向异性特征。在各向同性介质中，纵、横波以各自的射线速度(群速度)独立传播。在VTI介质中，纵、横波相互耦合，并且在一般情况下波的极化方向与传播方向既不垂直也不平行，因此VTI介质中的纵、横波又被称为qP波、qSV波和qSH波，相互耦合的纵、横波传播波动方程为[28]

(1)

式中:U=(ux,uy,uz)T为波场位移，ux、uy、uz分别为U在空间坐标x、y、z的分量；ρ为介质密度；F=(fx,fy,fz)T为单位质量元素的体力，fx,fy,fz分别为F在x、y、z的分量；t为时间；C为VTI介质刚度矩阵；L为偏导数算子矩阵。

式(1)描述了VTI介质中各质点在不同时刻的位移情况和弹性波传播规律。虽然式(1)的每一项都有明确的物理意义，但并没有体现波在VTI介质中的传播特性，如相速度、群速度或频散，因此有必要推导体现波在VTI介质中传播的特性方程。去掉式(1)的体力项，其平面波的解为

U=Aexp[ik(n·x-vt)]

(2)

(3)

其中

(4)

式(3)、式(4)中：Ax、Ay、Az分别为A在x、y、z的分量；kx、ky、kz分别为k在x、y、z的分量；c11、c66、c44、c12、c13、c33分别为C中不同位置的元素。

式(3)即为Kelvin-Christoffel方程，在数学上描述了本征值问题，是关于ρω2的一元三次方程。在VTI介质中，给定波场传播方向，由式(3)得到3个相速度解，分别对应qP波、qSV波、qSH波。地震波在各向异性介质中传播时纵、横波是相互耦合的，在计算qP波垂直波数时，qSV波速度影响较小，可以忽略。为了从式(3)中分离qP波，令横波速度为零，进行拟声波近似[4]。将Thomsen各向异性参数VP0、ε、δ[29]代入式(4)，有

(5)

将式(5)代入式(3)求解，可得VTI介质qP波频散方程，即

(6)

如果仅考虑xoz平面内的情况，求解式(6)可得二维VTI介质qP波垂直波数表达式

(7)

当ε=δ=0时，对应各向同性介质；当ε=δ≠0时，对应椭圆各向异性介质；当ε≠0,δ≠0且ε≠δ时，对应一般各向异性介质。

本文的角度域偏移采用双平方根(DSR)方程定义波场U的传播特征，在频率—波数域可表示为

(8)

式中:km、kh分别为中点、半炮检距分量的波数；kz为同时包含炮点和检波点的垂直波数，即

kz=kzs+kzg

(9)

由VTI介质的频散方程得

(10)

式中kmh=(km∓kh)/2，当kz=kzg时取正号，当kz=kzs时取负号。将式(9)、式(10)代入式(7)，可得VTI介质角度域偏移的波场传播方程。在同一延拓层内假设介质纵波速度和各向异性参数不变，则可得对应的相移偏移波场延拓算子

U(km,kh,z+Δz,ω,ε,δ)=U(km,kh,z,ω,ε,δ)×

exp[iΔz(kzs0+kzg0)]

(11)

式中kzs0、kzg0为v=v0、ε=ε0、δ=δ0时的背景波数，可表示为

(12)

相移法适用于横向均匀层状介质。为了解决复杂构造的成像问题，本文借用双域法的研究思路处理横向非均匀VTI介质。将横向非均匀VTI介质分解为横向均匀的速度、各向异性参数及其扰动量，在频率—波数域处理横向均匀的速度及各向异性参数背景场，然后在空间域基于速度及各向异性扰动对波数域计算结果进行时移校正。

将速度及各向异性参数分解为背景场(·)和扰动场Δ(·)，其中背景参数由模型平滑而得，各参数的扰动量同样包括炮点和检波点两部分。为方便推导，统一讨论各参数炮点和检波点分量，在实际计算过程中将其分开。各参数分量表示为

(13)

则VTI介质的垂直波数kz可分解为各向异性背景介质垂直波数kz0和扰动介质散射波场垂直波数kr

kz(x,z,km,kh)=kz0(z,km,kh)+kr(x,z,km,kh)

(14)

则

kr(x,z,km,kh)=kz(x,z,km,kh)-kz0(z,km,kh)

(15)

将kz、kz0在各向异性背景参数处展开，在小扰动假设下，忽略各向异性参数扰动的二阶及以上项，整理后可得

kr(x,z,km,kh)=as(z,kmh)Δs(x,z)+

aε(z,kmh)Δε(x,z)+aδ(z,kmh)Δδ(x,z)

(16)

其中

而

根据各向异性相关理论可知：在各向异性介质波数展开式中，慢度扰动项Δs(x,z)主要影响低阶次项，表现为对小角度传播波场的影响，各向异性参数扰动项主要影响高阶次项，表现为对高角度传播波场的影响[5]。因此总的扰动介质垂直波数kr可由三部分组成

kr(x,z,km,kh)=kΔs(x,z,km,kh)+

kΔε(x,z,km,kh)+kΔδ(x,z,km,kh)

(17)

结合式(13)，可将VTI介质波场传播算子P分解为背景波场的相移传播算子PPS、慢度扰动传播算子PΔs、ε扰动传播算子PΔε和δ扰动传播算子PΔδ。即

P=PPS+PΔs+PΔε+PΔδ

(18)

其中

PPS=exp(ikz0Δz)exp(iωΔsΔz)

PΔs= i[as(z,kmh)-as(z,0)]ΔsΔz×

exp(ikz0Δz)exp(iωΔsΔz)

PΔε=iaεΔεΔzexp(ikz0Δz)exp(iωΔsΔz)

PΔδ=iaδΔδΔzexp(ikz0Δz)exp(iωΔsΔz)

可以看出P为双域传播算子。首先基于速度扰动校正相移背景波场，还考虑了各向异性参数扰动对高角度入射波场传播的影响，并对相移结果进行进一步补偿，较好地提高了非均匀各向异性介质中波场延拓精度。

综上所述， 将VTI介质角度域偏移的波场延拓过程分述如下。

(1)频率—波数域相移。

U1(km,kh,z+Δz,ω,ε,δ)

=U(km,kh,z,ω,ε,δ)exp(iΔzkz0)

(19)

其中

相移处理实现了背景波场在VTI介质中的传播。

(2)频率—空间域时移。

U2(km,kh,z+Δz,ω,ε,δ)=U1exp[iΔzωΔs]

(20)

(3)频率—空间域各向异性参数校正。

U(km,kh,z+Δz,ω,ε,δ)

=U2exp{iΔz[as(z,kmh)-as(z,0)]Δs+

aεΔε+aδΔδ}

≈U2{1+iΔz[(as(z,kmh)-as(z,0))Δs+

aεΔε+aδΔδ]}

(21)

校正由各向异性参数扰动引起的大角度入射波场传播误差。

波场延拓后，对成像空间的地震波场进行数据映射获得角度域共成像点道集[22-25,30]。对于常速介质，反射波旅行时t为炮点、检波点位置和成像深度的函数(图1)，即t=t(z,s,g)。t对z求导，得

(22)

成像点局部半炮检距射线参数可表示为

(23)

由以上两式得

(24)

在频率—波数域，有

(25)

式中kh和kz分别为局部半炮检距波数和局部深度波数。式(25)即为成像空间炮检距域向角度域的映射公式，将波场由kh—kz域变换到β—kz域。

图1 常速介质中的射线传播路径示意图

在完成每步波场延拓后，首先提取零时间非零炮检距成像值，得到炮检距域共成像点道集(ODCIG)，局部成像公式为

(26)

I(m,h,z)⟹I(m,β,z)

(27)

式中：“⟹”表示由炮检距域向角度域的转换；I(m,β,z)表示坐标(m,z)处的ADCIG，它是空间位置和入射角的函数，因此可以解决常规成像道集遇到的多路径问题。通过上述方法获得不同CMP点的ADCIG，将某一入射角范围内的所有道集叠加便得到一定角度的叠前深度偏移剖面。

图2展示了水平层状模型ADCIG提取过程。首先使用局部成像条件获得波动方程ODCIG(图2b)，ODCIG在成像点深度聚焦为一个点；然后在频率—波数域利用局部炮检距和入射角的映射关系，将局部炮检距转化为角度信息，最终获得ADCIG(图2c)，ADCIG在成像点深度为拉平的直线，同时包含了丰富的角度信息，具有更高的信噪比和分辨率。

图2 水平层状模型ADCIG

为了验证VTI介质角度域偏移方法，设计了VTI介质模型(图3)进行正演，图4为正演记录。由图可见，炮集覆盖范围分别对应高速异常体和各向异性异常体位置，反射波较明显，较好地反映了地下构造，但由于存在异常体及断层，炮集记录中同时存在严重的绕射波干扰。

图3 VTI介质模型

图4 正演记录

以此数据为基础，分别进行常规角度域偏移和考虑各向异性影响的VTI介质角度域偏移。图5为各向同性、各向异性角度域叠前深度偏移剖面。由图可见：①由于未考虑各向异性参数的影响，各向同性偏移成像效果相对较差，浅层高速异常体边界成像模糊，分辨率较低；反映各向异性异常体的反射同相轴能量很弱，几乎没有反映下边界；对于右侧的断层成像，断点绕射波没有完全收敛(图5a)。②各向异性偏移提高了高速异常体的成像分辨率，边界刻画清晰；各向异性异常体成像能量增强，且明显反映了下边界；提高了断层区域成像信噪比，断点绕射波收敛较好，构造特征更清晰；深层反射波能量增强，提高了信噪比、分辨率，明显改善了总体成像效果(图5b)。

为了更直观地展示各向异性偏移的优势，分别提取具有典型构造特征位置处的ADCIG进行对比。图6为由各向同性、各向异性偏移方法提取的ADCIG。由图可见：①在两个异常体位置处，各向同性偏移ADCIG同相轴无法校平，存在不同程度的向上弯曲(图6a左、图6a中)；各向异性偏移ADCIG同相轴较为平直，而且可观测到更宽的角度范围(图6b左、图6b中)。②在断层处各向同性偏移ADCIG同相轴短小、弯曲、连续性较差，在局部区域还产生多组同相轴，出现构造假象(图6a右)；各向异性偏移ADCIG同相轴平直，消除了同相轴假象，剖面信噪比大大提高，较准确地反映了真实的构造特征(图6b右)。

图5 各向同性(a)、各向异性(b)角度域叠前深度偏移剖面

简单模型的成像测试证明了本文的VTI介质角度域偏移方法的可行性，可以处理VTI介质成像问题。

为了进一步验证本文的VTI介质角度域叠前深度偏移方法的成像效果，选取SEG 2D Hess VTI模型(图7)进行偏移测试。图8为各向同性偏移和VTI介质偏移剖面对比。由图可见：①各向同性偏移虽然可以对模型进行构造成像，并可以刻画各个界面，但成像分辨率较低，部分绕射波没有完全收敛，深层同相轴能量较弱，特别是高速盐丘体边界成像模糊，并在盐丘上部出现构造异常；无法对上部透镜体正确成像，中、下部的透镜体成像杂乱，分辨率较低；大断层处存在断面绕射波，成像位置出现偏差，整体成像效果较差(图8a)。②VTI介质角度域偏移对各个层位的成像更清晰，清晰地刻画了盐丘底界面、三个透镜体，大断层的断点清晰，断点绕射波收敛较好，断层位置更精确，明显提高了剖面信噪比(图8b)。图9为透镜体区域各向同性和各向异性偏移结果对比。由图可见，各向异性偏移消除了构造假象，绕射波完全收敛，成像精度更高(图9b)。

图6 由各向同性(a)、各向异性(b)偏移方法提取的ADCIG

图7 SEG 2D Hess VTI模型

图8 各向同性偏移(a)和VTI介质偏移(b)剖面对比

图9 透镜体区域各向同性(a)和各向异性(b)偏移结果对比

Hess VTI模型测试结果表明，本文的VTI介质角度域偏移方法可以较好地处理复杂地质构造的成像问题。

图10为各向同性和各向异性偏移提取的ADCIG。由图可见：各向同性偏移的ADCIG同相轴无法完全校平，还残留部分绕射波，信噪比和分辨率较低(图10a)；VTI介质偏移的ADCIG同相轴较为平直，信噪比和分辨率较高，进一步说明VTI介质角度域偏移的可行性(图10b)。

图10 各向同性(a)和各向异性(b)偏移提取的ADCIG

受各向异性参数建模技术的制约，VTI介质偏移方法在实际数据处理中存在困难。为了验证本文方法对实际数据处理的适应性，将偏移速度分析获得的速度模型作为各向异性偏移中的qP波速度，通过

(28)

求取各向异性参数ε、δ[15]。计算获得的ε、δ最大值分别为0.45和0.36，在地球介质各向异性参数取值范围内。将以上参数用于实际数据的VTI介质角度域偏移。图11为各向同性偏移和VTI介质偏移剖面。

图11 实际数据各向同性偏移(a)和VTI介质偏移(b)剖面

图12 CDP560处提取的ADCIG

由图可见，各向同性偏移与各向异性偏移结果的构造形态基本一致，但后者在成像细节方面有所改进(图11b)，表现为：剖面左侧断点绕射波进一步收敛，断层下盘反射波能量增强，断面刻画更清晰；剖面右侧向斜区域成像清晰，中深层反射波能量得到增强，主要层位的反射同相轴连续性更好，提高了剖面分辨率。图12为CDP560处提取的ADCIG。由图可见，VTI介质偏移的角道集同相轴更平直，连续性更好，可观测到更宽的角度范围(图12b)。

将VTI介质qP波的频散关系引入双平方根方程，以各向同性偏移的波场延拓为基础，推导了VTI介质双平方根偏移的波场延拓算子，考虑了各向异性因素对波场传播引起的扰动，在波场延拓过程中进行频率—空间域误差校正。通过对延拓波场进行局部炮检距成像和角度转化，最终实现VTI介质角度域叠前深度偏移。

模型试算和实际资料处理结果表明：各向同性偏移方法由于未考虑各向异性参数的影响，绕射波不能完全收敛，波场聚焦效果差，降低了成像剖面的分辨率和信噪比，不能对地质构造精确成像；VTI介质角度域偏移可对断层、盐丘、小尺度地质体精确成像。对于角度域偏移产生的ADCIG而言，各向同性偏移的ADCIG同相轴无法校平，残留断点绕射波，波场无法正确聚焦，不能正确反映局部地质特征；VTI介质偏移的ADCIG同相轴较平直，角度范围更宽，波场归位准确，精度较高。因此，利用VTI介质角度域偏移方法可对复杂构造精确成像。

石油地球物理勘探2019年2期