各向异性地层深度域速度建模

李汉青,李向阳,侯 嵩

(1.中国石油大学(北京) CNPC物探重点实验室，北京 102249；2.CGG公司英国技术中心,英国 雷德希尔 RH11PP)

偏移成像是常规地震数据处理的最后一步，也是十分关键的步骤。但在实际处理中地下构造复杂且目标储层深度较大，叠后时间偏移成像方法效果往往不佳。叠前深度偏移(Pre-stack Depth Migration，PSDM) 是解决速度横向变化剧烈的复杂地区的地震资料成像的关键技术。它可以使散射波收敛，并将地震反射波移动到其在深度域中的正确位置。它可以提供比其他成像方法更好的复杂地下结构图像[1-2]。

当前对于纵波和转换波的多分量偏移处理，是通过波场分离将数据分离为纵波和横波场，然后基于声学假设对分离后的波场进行纵波和横波标量偏移成像[3-4]。这种处理方法没有考虑弹性矢量波场的特性，在波场分离的过程中，不同波形的能量不能完全分离。剩余的非纯波能会在偏移结果中产生大量的串扰噪声，成像效果严重降低[5]。

基于矢量波场延拓的弹性波偏移成像技术为解决此类问题提供了有效方法。Pao和Varatharajulu[6]推导出了弹性波Kirchhoff-Helmholtz积分形式。秦福浩等[7]在此基础上推出了弹性波Kirchhoff积分偏移公式。荣骏召等[8]从多分量地震波矢量场的特征出发，基于均匀各向同性弹性波波动方程，利用弹性波波动方程的矢量特性，在偏移过程中将泄漏的地震波信息还原到PP波与PS波分量中，实现了偏移过程中的波场分离。

与其他偏移方法相比，PSDM需要更精确的速度模型。传统的速度模型更新方法可能不够精确。叠前深度偏移本身对速度模型的敏感性高，是建立速度模型的有效工具。偏移速度分析有两个主要原则：深度聚焦分析(Depth-Focusing Analysis，DFA)和剩余曲率分析(Residual-Curvature Analysis，RCA)。深度聚焦分析是在利用叠加功率测量速度误差的基础上进行的，剩余曲率分析是基于剩余时差来测量速度误差的。在DFA中，如果偏移深度和聚焦深度的差为零，偏移速度是可以接受的。在RCA中，如果不同偏移量的成像深度差为零，偏移速度是可以接受的[9-10]。

再者，由于实际地下介质具有弹性各向异性，因此在深度偏移过程中需要对介质中的各向异性进行补偿。Crampin等[11-14]的研究发现，地下介质中存在广泛的各向异性。岩石的层状结构、定向排列的裂隙都会引起各向异性。横波对于各向异性的敏感，使其在多波地震资料的处理和解释中有了重要的价值。转换波速度建模比纵波速度建模难度更大。因为，对于P-S转换波，下行波为P波，上行波为S波，这使得速度分析可靠性降低。在各向异性介质中，波的传播速度不仅是空间坐标的函数，而且与入线方向有关。如果不考虑各向异性的走时，共成像点(Common Image Point，CIP)道集在大偏移量下将不会被拉平，从而影响成像质量。为此提出了一种估计地下各向异性信息并将其应用于叠前深度偏移的方法。

针对这些不足，提出了基于矢量波波场延拓的弹性波叠前深度偏移速度分析和各向异性参数的估计方法，更新了速度结构，建立了纵波和转换波速度模型。该方法基于RCA的最基本原理，避免了数学逼近，可用于复杂速度结构的建立。

1 弹性波Kirchhoff积分法

鉴于传统多分量处理方法没有考虑弹性矢量波场的特征，剩余能量会造成偏移结果中的噪声干扰。本文采用矢量波波场延拓的弹性波偏移成像技术。弹性波Kirchhoff-Helmholtz积分方程，可以同时向下延拓P波和PS波的波场。Pao等[6]推导出了弹性波Kirchhoff-Helmholtz积分形式[7]。Kuo等[15]将此方法应用于单界面模型数据的弹性波Kirchhoff叠前深度偏移，同时输入X、Y、Z三分量的单炮数据，分别得到X、Y、Z三个方向的偏移剖面。在自由表面边界条件下，地表应力为零，可以得到简化的偏移方程：

(1)

下标标注方式满足爱因斯坦求和约定。

根据射线出射角将矢量位移分别投影到PP分量和PS分量：

(2)

其中，θ为位移水平分量与x轴夹角；φ为位移与轴夹角。考虑二维多分量勘探，将uy略去，得到二维弹性波Kirchhoff积分法矢量偏移公式：

(3)

(4)

实际介质是各向异性的，地下绕射点到地表接收点的PP波与PS波的射线路径不再是直线，两者的射线路径、走时和出射角均不相同。因此需要对均匀各向同性介质偏移公式进行改进。走时t，tP，tS，路径出射角θP、θS和等效纵横波速度VP、VS通过各向异性介质的射线追踪和速度分析得到。求得方程(3)、(4)写为：

(5)

(6)

式中，tp为上行P波走时，ts为上行SV波走时；SP、SS分别为上行P、SV波的射线路径长度；θP、θS分别为P、SV波射线出射角；VP、VS分别为P、SV波速度。

与基于声波方程的Kirchhoff积分偏移算法比较，弹性波Kirchhoff积分矢量偏移算法具有很多优点：具有波场分离的特性，偏移时同时分离纵波和转换波；偏移算法具有保幅性；各向异性射线追踪和速度分析使算法适用于各向异性介质。

2 深度偏移速度更新方法

利用叠前深度偏移对速度模型的敏感性，将它作为一种有力的速度分析工具。在进行速度建模时，可以通过叠前时间偏移结果和Dix公式得到一个初始的速度模型。在这个模型基础上给定一个初始速度范围和步长，通过叠前深度偏移得到不同速度下的CIP道集。

剩余曲率分析是根据剩余旅行时差测量速度误差，将剩余曲率分析原理应用于共成像点(Common Imaging Point，CIP)道集，得到共成像点道集拉平准则，即当偏移速度准确时共成像点道集上不同炮检距的成像深度相同；当偏移速度偏低时，同相轴上翘；反之，偏移速度偏高时，同相轴下垂。

对于各向异性地层，Sarkar和Tsvankin[16-18]给出各向异性近似下纵波成像深度与炮检距关系

(7)

同理，上述将方程拓展到SV波，方程形式与纵波方程相同[2]，可以表示为

(8)

首先按照介质为各向同性时进行速度分析，此时r2为0。本文采用的速度更新方法就是在此基础上提出的。速度更新步骤如下：

(1)在一定的速度范围内，通过叠前深度偏移构建CIP道集。

(2)通过参数r1拉平CIP道集，并拾取每个速度值所对应的r1的值(图1)。

图1 不同偏移速度得到的CIP道集的r1参数谱Fig.1 r1parametric spectrum of CIP gathers obtained at different migration velocities

(3)通过线性回归的方法计算r1=0时的速度值，用这个速度作为估计速度更新速度模型。

3 深度偏移各向异性参数反演

在各向同性理论下，地震波速度不随地下介质空间坐标变化；而在各向异性介质中，地震波速度随着射线角度的变化而改变。

Thomsen[13]用5个各向异性参数VP0、VS0、ε、δ、γ描述了速度的弱各向异性，VP0、VS0是垂直纵波和横波速度，横向各向同性介质(Vertical Transverse Isotropy，VTI)中地震波速度与这5个参数的关系式：

VP(θ)=VP0(1+δsin2θcos2θ+εsin4θ)

(9)

(10)

VSH(θ)=VS0(1+γsin2θ)

(11)

上一节方程(7)给出了各向异性地层偏移速度、偏移深度和真实速度的关系。由方程可以看出CIP道集上的剩余时差是非双曲函数，其大小由参数Vnm0和η决定。其中Vnm0决定近偏移距时差(二次)，而η对时差的影响在远偏移距变得更大(4次)。如果偏移值和真实值一致，则偏移深度z(h)与半跑间距h无关，即同相轴变平。虽然方程(7)是假设反射层水平时推导而来，但Sarkar证明即使在地层倾角较大时也成立[18]。

为了定量分析方程(7)中各向异性参数ε、δ对偏移深度的影响，建立含有2个反射层的各向异性介质速度模型，第1层是倾斜层，第2层为水平层。第1层和第2层速度以及各向异性参数相同。VP0=2 000 m/s，ε=0.3，δ=0.1，则η=ε-δ/1+2δ=0.167。运用CWPSU程序susynlvfti正演合成数据。再运用第一节所述基于矢量波波场延拓的弹性波偏移技术进行偏移得到CIP道集。由于各向异性参数的存在，偏移的关键步骤是各向异性射线追踪。这里运用SU程序ray2dan进行各向异性射线追踪计算走时和出射角，该程序依据旁轴近似将每条射线计算的旅行时外推到相邻的网格点上。

图2给出了在VP0=2 000 m/s正确时，各向异性参数ε、δ的错误取值和正确取值时的CIP道集。当选用正确的参数ε=0.3,δ=0.1进行偏移时，CIP道集同相轴是平直的(图2(c))。图2(a)是假设地层为各向同性时(即ε=δ=η=0)的偏移结果，由于忽略了各向异性参数的影响，同相轴上翘。图2(b)是各向异性参数欠估计时的偏移结果(ε=0.3,δ=0.05，η=0.136)，与之对应图2(d)是各向异性参数过估计时的偏移结果(ε=0.4,δ=0.15，η=0.192)，同相轴下倾。

图2 各项异性参数对CIP道集同相轴的影响Fig.2 Influence of anisotropy parameters on CIP gathers

这种速度各向异性在CIP道集上表现为远偏移距的同相轴无法拉平。由方程(9)、(10)可知，P波和SV波在传播过程中是耦合在一起的，同时由Thompson参数ε、δ决定。因为P波剖面有较强的能量和较高的信噪比，所以本文通过纵波反演Thompson参数ε、δ。

如上分析，当地层存在各向异性时，即使VP0准确，在偏移获得的CIP道集中，成像深度仍然随炮检距变化。由于各向异性的存在，不同偏移距应当采用不同的偏移速度。因为成像点Zh和Z0的走时相同，有

(12)

这里Vθ即公式(9)中的VP(θ)，将方程(9)与(11)联立，得

(13)

方程(13)可改写为

Aδ+Bε=C

(14)

(15)

(16)

(17)

A、B、C的值与成像深度和炮检距有关，成像深度和相对应的炮检距可以从CIP道集上拾取，如图3、图4所示。

图3 成像深度为炮检距的函数Fig.3 Relation between imaging depths and offsets

图4 拾取成像深度与炮检距Fig.4 Picking imaging depths and offsets

然后，根据公式(17)计算得到列向量A、B、C。最后通过求解广义逆矩阵的方法，通过公式(16)求解Thompson参数ε和δ。

同理，对于Thompson参数γ，可以在纯横波CIP道集上进行各向异性参数分析。

(18)

γ=A-1·B

(19)

计算得到列向量A、B(图4)。再通过求解广义逆矩阵的方法，通过式(19)求解Thompson参数。

4 实际工区模型测试

以胜利油田罗家工区Ken71模型进行试验。Ken71模型为复杂的地质构造模型。该地区具有断裂系统复杂、微幅构造难以识别、薄互层发育等地质特征，是典型的河流相与三角洲相沉积类型油藏[19]。针对薄互层地质特点，建立了典型的Ken71地质模型并进行了正演模拟。各向异性速度建模和各向异性参数分析流程如图5所示，P波和S波速度场模型如图6所示。

图5 各向异性速度建模和各向异性参数分析流程Fig.5 Anisotropic velocity modeling and anisotropic parameter analysis process

图6 Ken71速度场模型Fig.6 Velocity model of Ken71

图中红色括号标示出目标薄互层位置，1 480～1 630 m，其各向异性参数为ε=0.20,δ=0.05，γ=0.10。网格大小为2 766×901，横向的采样间隔为10 m，而纵向的采样间隔为2 m。正演弹性波记录共有75炮，每炮202道接收，道间距为20 m，此模型时间采样点共有3 000个，采样点间时间的间隔为1 ms[20],从中可以看出地质构造的复杂性。

正演弹性波场的X分量和Z分量单炮记录如图7所示，可以看出由于地下地质构造的复杂性导致地震记录信息也非常复杂。

图7 Ken71模型正演记录Fig.7 Synthetic shot gathers of Ken71

运用图5提出的速度建模和各向异性参数分析流程，得到纵波和转换波速度模型，如图8所示。在速度更新时为简化层剥过程，只拾取了7个反射界面，逐层计算了纵波和转换波的速度。

图8 速度更新得到的速度模型Fig.8 Velocity model after velocity analysis

在假设速度模型准确的基础上对第六层深度为1 400～1 680 m处当作各向异性地层(VTI)处理，反演得到各向异性参数平均结果约为ε=0.22,δ=0.05，γ=0.12。深度偏移得到X=3000 m处纵波CIP道集如图9所示，其中图9(a)为忽略各向异性(即各向同性ε=δ=γ=0，偏移得到CIP道集，图9(b)为利用各向异性参数ε=0.22,δ=0.05、偏移得到CIP道集。对比二图看到各向异性建模偏移的CIP道集拉平效果较好。

图9 深度偏移得到的纵波CIP道集Fig.9 P-wave CIP gathers from depth migration

类似图9，利用各向同性和各向异性速度建模最终偏移的纵波成像结果和转换波成像结果如图10和图11所示。

图10 深度偏移得到的叠加剖面Fig.10 Stacked profile from depth migration

图11 转换波深度偏移得到的叠加剖面Fig.11 Stacked profiles obtained by depth migration of converted waves

深度偏移得到的转换波CIP道集如图12所示，其中图12(a)为各向同性偏移得到的CIP道集，图12(b)为利用各向异性参数γ=0.12偏移得到的CIP道集，在VTI层(上部红色箭头处)，二者拉平效果接近，但对深层区别较大(下部红色箭头处)。纵横波成像深度基本一致。各向异性参数对实际资料拉平效果没有前述模拟资料明显，主要原因可能是实际资料各向异性参数较小，另外实际资料的信噪比较低和Ken71地层复杂的结构有关。但总体拉平效果有所改善，估计的参数值可以接受。

图12 深度偏移得到的转换波CIP道集Fig.12 Converted wave CIP gathers obtained by depth migration

图12中圆圈处显示偏移结果得到了提高，对于纵波，同相轴连续性变得更好，对于转换波由于CIP道集拉平，叠加结果使叠加剖面同相轴更清晰。图13是各向异性速度建模最终偏移的纵波成像结果和转换波成像结果对比。

图13 Ken71地震数据纵波和转换波深度偏移剖面对比Fig.13 Comparison of P-wave and converted-wave depth migration profiles of Ken71 seismic data

从图13中可以看出，两者对地下复杂的地质构造例如起伏地表、断层等都具有不错的成像效果，结合地震解释资料，深度剖面上对含油透镜体、油水薄互层、油顶油水同层以及油气水同层等主要层位都进行了准确的刻画，转换波的分辨率高于纵波。可以看到在深度1 400～1 680 m处正好是油水薄互层，薄互层一般等效为VTI地层，与前文设定的各向异性层对应，层位偏移归位到真实位置。

5 结论

(1)本文利用Kirchhoff积分叠前深度偏移公式，考虑各向异性介质的射线追踪，实现了各向异性地层矢量偏移，输入多分量地震数据，同时分离得到纵波和转换波的成像剖面。

(2)之后分析了纵波和转换波深度域地层速度模型建模方法。本方法基于RCA原理，在CIP道集上进行速度更新，并将偏移速度与真实速度的差异参数化，通过拾取能量团的方式估计真实速度。然后对目标层进行叠前深度偏移，拾取速度结构，实现速度建模。在速度模型准确的基础上，本文提出了新的各向异性参数分析方法，推导出了各向异性参数、偏移距与成像深度之间的关系式，通过求解广义逆矩阵的方法反演各向异性参数，实际数据测试取得了较好的效果，验证了方法的可行性。