2011, Vol. 54
波动方程叠前深度偏移能适应横向变速强烈情况，在复杂地质构造地区产生很好的成像.但是，波动方程叠前偏移计算量很大.为此，许多学者研究了平面波偏移在减少计算量方面的潜力.Berkhout[1]提出面炮偏移；张叔伦等[2]实现了傅里叶有限差分平面波偏移；陈生昌等[3]实现平面波深度偏移；Zhang等[4]实现3D 延迟放炮偏移；Stoffa等[5]对共检波点道集和共炮道集分别作平面波偏移再叠加以提高成像精度；崔兴福等[6]利用平面波偏移进行了分角度成像和AVA 分析；叶月明等[7]研究了起伏地表下的平面波偏移；Shan 等[8]提出倾斜坐标系下平面波偏移以提高角度成像准确度；王华忠等[9]实现Offset平面波有限差分法叠前时间偏移；陈生昌和马在田[10]从地震波场的线性叠加角度从广义上对这类地震数据合成作了说明.平面波偏移还可以自然的产生共入射角角度道集，用来进行速度更新[4].但之前关于平面波偏移的研究基本都是在各向同性介质下，没有考虑各向异性介质情况.
各向异性在石油勘探中普遍存在，且在地震成像中变得越来越重要[11， 12].VTI(VerticalTransverseIsotropy)介质是一种最简单也是很实用的各向异性介质近似形式.但VTI模型只对简单的水平沉积地层有效.当地层为倾斜薄层类型构造时，如加拿大落矶山脉丘陵地带以及靠近盐丘体位置[11]；或在各向同性介质中嵌入的并行的倾斜断裂情况[13]，介质的对称轴不再是竖直的，不能由VTI介质表征，而需要将其考虑为TTI(TiltTransverseIsotropy)介质.TTI介质层在各向同性或VTI偏移中会引起严重的畸变[11， 14].为了对TTI覆盖层下的目标储层准确成像，偏移算法必须能够处理倾斜对称轴.
与VTI介质相比，TTI介质的频散关系要复杂的多.各向同性的频散关系很简单，且有显式的表达；VTI介质和TTI介质频散关系都是一个四次方程，但VTI介质的频散关系是一个偶函数，假设S波速度为零，可以导出相应的频散关系显式表示；而TTI介质频散关系需要一个偶函数和一个奇函数共同表示，不能导出相应的频散关系显式表达，需要数值的方法解.传统的隐式有限差分[15-18]需要Taylor序列分析和Pade展开近似显式频散关系，但对于TTI介质很难导出其频散关系的Taylor序列.Shan[19]提出了优化的方法，将非线性问题降解为线性优化问题，再利用线性优化解出的系数求解差分系数.Shan 的方法不受弱各向异性假设的限制.本文沿用优化的思想，首先用解析方法解VTI介质频散关系，再用数值的方法将其旋转到TTI介质情况，然后直接建立非线性优化问题，使用非线性最小二乘求解差分系数，得到TTI介质中IFD 外推算子.
文章首先介绍了各向同性介质下的平面波偏移理论，接着阐述了TTI介质中IFD 算子系数设计，结合二者，将平面波偏移推广到TTI介质中.文章将设计的TTI波场外推算子和偏移脉冲响应分别与理论解作了对比，并用TTI弹性波方程拟合，并用极化滤波分离出的P 波数据验证TTI介质平面波偏移的准确性和高效性.
2 各向同性介质平面波偏移常规叠前偏移是通过对单个炮集进行独立的偏移得到单炮成像，并通过叠加所有的单炮成像得到整个地下成像.波动方程单炮偏移过程需要两步:首先，将震源波场和接收波场从地表外推到地下所有的深度；其次，在每个深度通过对震源波场和接收波场应用成像条件(如做互相关)得到成像.
炮集记录可以合成一个新的数据体，来表征一个现实中没有进行的物理实验.这种数据合成一个重要的例子就是组合炮集记录建立平面波震源道集.数学上可以通过对共检波点(或共炮点)道集进行倾斜叠加(线性Radon变换)实现:
|
(1) |
其中，R=R(sx ，x，z，ω)是接收波场，sx 是震源位置，psx 是相应于sx 的射线参数，rx 是地表检波器位置.相应的地表平面波震源波场为:
|
(2) |
和常规波动方程共炮法偏移一样，典型的平面波偏移也需要两步.首先，将震源波场Sp 和接收波场Rp 分别用单程波动方程独立的外推到地下所有深度点；其次，含有射线参数Psx 的平面波成像通过对震源波场和用圆频率ω 加权了的检波波场作互相关实现:
|
(3) |
其中Sp* 是震源波场Sp 的复共轭.整个图像通过叠加所有Psx 值的平面波成像得到:
|
(4) |
因为倾斜叠加和偏移都是线性算子，在连续情况下，平面波偏移与传统的炮剖面偏移效果一样[2].在实际离散形式下，需要有足够数量的P值使两种成像等同.Zhang等[4]讨论了在平面波偏移方法实际应用中需要的共P剖面数量Np 及分量采样步长Δp:
|
(5) |
其中，f=ω/2π，Ns为一个共检波点道集中的炮数，Δxs 为炮间距，所以NsΔxs 为一个共检波道集的长度，入射角α1 ≤αs ≤α2，vs 为地表速度.
3 TTI介质中IFD 算子系数设计 3.1 TTI介质中的频散关系在VTI介质中，相速度qP-和qSV-用Thomsen参数可以表示为[20]:
|
(6) |
其中，θ 为波的速度的传播角，即波传播方向与竖向轴之间的夹角；f=1-(vs0/vp0)2；vp0和vs0分别为qP-和qSV-波在竖向方向的速度.在式(6)中，当根号前符号为正时，v(θ)为qP-波的相速度，为负时，v(θ)为qSV-波的相速度.
根据式(6)，并假设f≈1可以导出VTI介质频散关系的显式表达[21]:
|
(7) |
将式(6)对称轴由竖直旋转到一个倾斜角度φ，可以得到TTI介质的相速度方程，对称轴与竖直方向成角度φ:
|
(8) |
这里与方程(6)式中不同的是，ε 和δ 变为与对倾斜对称轴方向相同的量.Vp0为qP-波平行于对称轴传播的速度.
通过式(8)及相位角θ 和波数的关系，可以得到TTI介质中的频散关系[18]:
|
(9) |
其中sz=kzvp0/ω，系数d0 ，d1 ，d2 ，d3 ，d4 分别为:
|
其中，${s_x} = \frac{{{k_x}{v_{{p_0}}}}}{\omega }$.TTI介质中的频散关系式(9)是一个有四个解的四次方程，其中两个解代表上行/下行qP-波，另外两解代表上行/下行qSV-波.我们只考虑qP-波的外推.
3.2 IFD算子设计虽然TTI 介质频散关系比各向同性介质和VTI介质要复杂的多，但设计IFD 算子模式基本相同.
式(9)是一个四次方程，但不像各向同性和VTI介质可以得到sz的显式表示，即使假设f=1，对于TTI介质，也很难得到sz以sx 为变量的直接表达式.传统的隐式有限差分需要Taylor 序列和Pade展开分析近似显式频散关系，因此，对于TTI介质很难应用Pade展开求解频散关系.本文采用解析的方法解VTI介质频散关系式(7)，并通过数值的方法将VTI解(sx ，sz)VTI 旋转一个角度φ 得到TTI介质频散关系解(sx ，sz)，再通过构建近似函数，用非线性最小平方匹配求解差分算子系数.
对于各向同性或VTI介质，sz是sx 的偶函数.可以用加权偶函数近似频散关系，如:
|
(10) |
表示.但对于TTI介质，sz不是sx 的对称函数.所以，除了偶函数之外，还需要奇函数近似频散关系，如sx ，sx3.本文使用如下一个加权函数近似TTI介质频散关系:
|
(11) |
其中，sz0 =sz(0)，系数ai，bi，ci，i=1，2，…，n，可以使用最小平方方法估计得到.式(11)整体保持着传统近似方法的形式，加入奇数项是为了突破各向同性或VTI介质中差分竖向对称的限制.
本文构建如下最小平方优化问题:
|
(12) |
式(12)是一个非线性优化问题，要用非线性最小二乘法解[22].系数ai，bi，ci，i=1，2，…，n是各向异性参数ε，δ 以及倾斜角φ 的函数，它们随水平位置变化.在波场外推过程中如果在每个网格都进行最小平方估计计算量太大，给定各向异性参数ε，δ 及倾斜角φ 范围，可以提前算出系数并存大表格中，在波场外推过程中直接从表格中查询代入有限差分算法即可.从表格中查找到系数后，再根据相应于式(11)的隐式有限差分模式式(18)(见附录)实现TTI介质有限差分，过程和计算量和各向同性下就相似了.
图 1是在各向异性参数为ε=0.2，δ=0.15，φ=45°时，对几种近似方法与真实的频散关系曲线进行了比较.“真实TTI曲线"是由式VTI频散关系式(7)计算并旋转φ=45°得到的TTI“真实"频散关系；“各向同性曲线"是不考虑各向异性影响时的频散关系；“VTI曲线"是由式(7)解析解出的频散关系；“优化的TTI曲线"为用本文建立的优化式(11)解出的频散关系，由图可以看出这种近似与真实频散关系几乎完全匹配，精度很高.几种近似方法对应的四阶(n=2)差分系数见表 1，其中优化的TTI除差分系数外还需要一个相移加权常数，这里sz0=0.9211.用同样的办法也可以拟合更高阶差分系数.
|
图 1 频散关系曲线比较 Fig. 1 Dispersion relation comparison |
|
表 1 几种方法求算子系数(ε=0.2，δ=0.15，φ=45°) Table 1 Coefficients of different approximately operators |
图 2 是在TTI参数为ε=0.2，δ=0.15，φ=45°，速度vp0 =2000m/s, 时间t=0.8s时几种方法的偏移脉冲响应，其中白线为根据相速度方程式(8)，及速度时间等参数算出的理论曲线(理论曲线单程按t=0.4s计算).图 2a为各向同性IFD80°差分算子对各向异性介质下的脉冲响应，从图中可以明显看出，忽略各向异性将会带来很大误差，尤其是在横向方向；图 2b为由VTI有限差分近似式(10)计算得到的差分算子的脉冲响应，从图中可以看出，忽略TTI介质对称轴倾角同样会带来较大误差；图2c为根据TTI频散关系近似式(11)，直接解本文构建的非线性优化问题式(12)得到的差分系数算子的偏移脉冲响应，从图中可以看出脉冲响应与理论曲线匹配的非常好，这就证明本文设计TTI隐式有限差分系数的方法非常准确.
|
图 2
几种方法的脉冲响应，vp0 =2000m/s, TTI介质参数:ε=0.2，δ=0.15，t=0.8s.白色线是理论曲线，由式(8)中相速度方程直接计算(近似f=1)
Fig. 2
Impulse responses of a tilted TTI data with ，vp0 =2000 m/s, ε = 0.2，δ=0.15 ， |
各向同性介质中，平面波偏移通过解上/下行波方程，再利用式(3)和式(4)的成像过程成像.在TTI介质中，平面波的合成过程和成像物理过程是不变的，不同的是应用的波传播方程不一样了，波场外推算子改为上面构建的TTI介质中IFD 外推算子.各向异性介质中波场外推算子的精度越高，成像越精确，但往往算子精度越高需要求解的优化问题计算量就越大，这就导致了各向异性偏移计算量一般要高于各向同性情况.应用文中的方法，使用四阶差分(n=2)，比各向同性计算量多出的部分是求取差分系数部分，而每个模型这种系数只需求一次，而且如果使用表格驱动模式，建立TTI隐式差分系数库后，所有模型只需查找系数时间即可，可以忽略.影响TTI介质差分系数的各向异性参数有ε 和δ，在建立或搜索系数库时需要遍历这两个参数即可.在求出各点差分系数后，计算量与各向同性介质情况下就基本相同了.另外，从应用平面波偏移这方面说，较常规各向异性偏移却又大大减少了偏移的计算量.
5 模拟算例考虑到TTI介质对称轴方向一般是横向变化的[23]，建立如图 3所示TTI介质模型.图 3a为竖向速度，图 3(b, c, d)分别为各向异性参数φ，δ，ε 的变化.TTI模型正演数据由全弹性波TTI波动方程正演程序正演[24]，然后利用极化滤波的方法将P 波与S波分离[25].图 4(a, b)分别为弹性波正演的z分量和x分量，图 4(c, d)分别为经极化滤波分离后的P波和P-SV 波场.我们只考虑P波偏移.拟合数据共200炮，道间距5 m, 炮间距20 m, 采样间隔1ms, 采样时长为2.4s.
|
图 3 TTI介质模型图 (a)沿对称轴方向的P波速度；(b)倾斜对称轴的倾斜角度φ；(c)各向异性参数ε；(d)各向异性参数δ. Fig. 3 TTI model (a) Velocity along symmetry axis； anisotropic parameters: (b)φ，(c)ε，(d)δ. |
|
图 4 极化滤波分离纵波和横波 (a)z分量数据；(b)x分量数据；(c)极化滤波后P波数据；(d)极化滤波得到的PSV 波数据. Fig. 4 P and S separation using polarization analysis method (a) z component； (b)x component； (c) P component and (d) PSV component. |
图 5为单个P剖面偏移过程.图 5(a, b)分别为p=0.2 m/km 时的震源波场和接收波场，图 5c为相应的偏移剖面.根据式(5)我们计算产生30 个共P面，其中入射角范围为-30°≤α≤30°，最终平面波偏移图像由所有30个单个共P 面偏移结果叠加而成，如图 6d所示.
|
图 5 p=0.2m/km 时平面波偏移 (a)平面波震源波场；(b)平面波记录；(c)相应的单个P偏移剖面. Fig. 5 Single plane wave migration process (a) Sourcewave-field； (b) Plane wave receiver wave-field；(c) Single p profile migration result. |
|
图 6 不同方法对拟合TTI介质数据偏移结果 (a)各向同性程序偏移结果；(b)VTI偏移程序偏移结果；(c)TTI偏移程序偏移结果；(d)TTI平面波偏移程序30个共P面偏移结果. Fig. 6 Migration results obtained by different programs (a) Isotropic media migration； (b) VTI media migration； (c) Conventional shot-profile TTI migration； (d) 30 numbers plane wave migration in TTI media. |
图 6为几种不同方法对TTI拟合数据的偏移结果.图 6a为各向同性隐式有限差分偏移结果；图6b为VTI隐式有限差分偏移结果；图 6c为用本文设计的TTI隐式有限差分算子计算的常规炮剖面法偏移结果，由201 炮偏移叠加得到；图 6d为本文中使用的TTI隐式有限差分平面波偏移结果，由30个共P面偏移结果叠加得到.将几种图与原始模型的对比，从图 6a和图 6b可以看出，各向同性和VTI偏移程序都不能很好地处理TTI介质产生的数据；从图 6d可以看出TTI平面波偏移是非常有效的，能准确偏移出原始模型构造.并且用少几倍的计算量达到了与炮剖面法(图 6c)质量相当的成像.在图6d中大倾角界面处能量稍弱，这是由于平面波偏移采集的参数P的范围都是-30°≤α≤30°，如果加入大角度的P值，大倾角界面能量可以得到加强.
6 结论TTI介质频散关系要比各向同性和VTI介质复杂，很难导出其频散关系的显式表示，因此也无法使用传统Taylor序列和Pade展开法进行隐式有限差分.本文通过解析解VTI介质频散关系，再用数学的方法旋转为TTI介质频散关系.传统的隐式有限差分频散关系近似函数为偶函数，而单独偶函数无法拟合TTI介质频散关系，文中在偶数基础之上加入一个奇函数，并建立非线性最优化问题求解TTI介质隐式差分算子.将设计的IFD 波场外推算子频散曲线与理论解对比，对比证明算子精度很高.对几种介质的偏移脉冲响应对比分析，证明了各向同性和VTI介质IFD 算子在TTI介质中都会带来很大误差，尤其是横向方向，而文中设计的TTI介质IFD 算子与TTI介质相速度方程解析解匹配的很好.从计算量方面考虑，TTI介质IFD 算子比各向同性算子多出的计算量部分为求取差分系数部分，在求取系数后，计算量和各向同性下基本相同.
文章将平面波理论与TTI介质IFD 算子结合实现TTI介质IFD 平面波偏移.将偏移算法用于拟合TTI数据，并与各向同性、VTI介质IFD 偏移及TTI介质IFD 常规炮剖面法偏移做了对比.偏移结果对比表明，各向同性和VTI介质IFD 偏移不能很好地对TTI数据成像，而应用设计的TTI介质IFD算子的常规炮剖面法偏移和平面波偏移都能很准确地偏移出原始模型构造.与炮剖面法相比，TTI介质平面波偏移的一个优势是在保证成像质量的同时，大大降低了计算量，这对各向异性偏移研究有重要意义.
附录 TTI介质有限差分模式在近似的频散关系方程(11)中，当四阶(n=2)时，根据逆Fourier变换将sz和sx 用偏微分算子$\frac{\partial }{{\partial z}}/\frac{\omega }{{{V_{P0}}}}$和$i\frac{\partial }{{\partial x}}/\frac{\omega }{{{V_{P0}}}}$代替，可以得到一个偏微分方程:
|
(A1) |
方程(A1)可以通过串级法(cascading)解:
|
(A2) |
|
(A3) |
|
(A4) |
方程(A2)可以通过在空间域相移实现，与各向同性不同的是它有个加权项Sz0.设Pin=P(ω，nΔz，jΔx)，其中Δx和Δz是有限差分网格大小.在方程(A3)中，离散形式为
|
(A5) |
用有限差分算子代替偏微分算子，有
|
导出如下有限差分方程:
|
(A6) |
其中，
式(A6)可以用Crank-Nicholson 法解[26]，当c1=0时，式(A6)与各向同性介质下隐式有限差分相同，所以式(A6)的计算量几乎和各向同性介质下相同.
致谢感谢张显文博士有意义的建议和讨论,郝奇博士提供的TTI介质全波正演程序,匿名审稿专 家细致且富有建设性的意见,在此一并致以诚挚的谢意．
| [1] | Berkhout A J. Areal shot record technology. J. Seis. Expl. , 1992, 1(3): 251-264. |
| [2] | 张叔伦, 孙沛勇. 基于平面波合成的傅里叶有限差分叠前深度偏移. 石油地球物理勘探 , 1999, 34(1): 1–7. Zhang S L, Sun P Y. Prestack depth migration based on both plane wave synthesizing and Fourier finite difference. Oil Geophysical Prospecting (in Chinese) , 1999, 34(1): 1-7. |
| [3] | 陈生昌, 曹景忠, 马在田. 叠前地震数据的平面波深度偏移法. 地球物理学报 , 2003, 46(6): 821–826. Chen S C, Cao J Z, Ma Z T. A method of plane wave depth migration for prestack seismic data. Chinese J.Geophys. (in Chinese) , 2003, 46(6): 821-826. |
| [4] | Zhang Y, Sun J, Notfors C, et al. Delayed-shot 3D depth migration. Geophysics , 2005, 70(5): 21-28. |
| [5] | Stoffa P L, Sen M K, Seifoullaev R K, et al. Plane-wave depth migration. Geophysics , 2006, 71(6): 261-272. DOI:10.1190/1.2357832 |
| [6] | 崔兴福, 刘卫东, 刘桂宝, 等. 平面波偏移、分角度成像与AVA道集生成. 石油物探 , 2007, 46(6): 615–620. Cui X F, Liu W D, Liu G B, et al. Plane wave migrating,separated angle imaging and AVA gather generating. Geophysical Prospecting for Petroleum (in Chinese) , 2007, 46(6): 615-620. |
| [7] | 叶月明, 李振春, 仝兆岐, 等. 起伏地表条件下的合成平面波偏移及其并行实现. 石油地球物理勘探 , 2007, 42(6): 622–628. Ye Y M, Li Z C, Tong Z Q, et al. Synthetic plane wave migration in relief surface condition and implementation in parallel way. Oil Geophysical Prospecting (in Chinese) , 2007, 42(6): 622-628. |
| [8] | Shan G J, Biondi B. Plane-wave migration in tilted coordinates. Geophysics , 2008, 73(5): 185-194. DOI:10.1190/1.2957891 |
| [9] | 王华忠, 冯波, 任浩然. 二维Offset平面波有限差分法叠前时间偏移. 石油物探 , 2009, 48(1): 11–19. Wang H Z, Feng B, Ren H R. Finite-difference pre-stack time migration of 2D offset plane wave. Geophysical Prospecting for Petroleum (in Chinese) , 2009, 48(1): 11-19. |
| [10] | 陈生昌, 马在田. 广义地震数据合成及其偏移成像. 地球物理学报 , 2006, 49(4): 1144–1149. Chen S C, Ma Z T. Generalized synthesis of seismic data and its migration. Chinese J. Geophys. (in Chinese) , 2006, 49(4): 1144-1149. |
| [11] | Vestrum R W, Lawton D C, Schmid R. Imaging structures below dipping TI media. Geophysics , 1999, 64(4): 1239-1246. DOI:10.1190/1.1444630 |
| [12] | Yan L, Lines L, Lawton D. Influence of seismic anisotropy on prestack depth migration. The Leading Edge , 2004, 23(1): 30-36. DOI:10.1190/1.1645453 |
| [13] | Grechka V, Tsvankin I. Characterization of dipping fractures in a transversely isotropic background. Geophysical Prospecting , 2004, 52(1): 1-10. DOI:10.1046/j.1365-2478.2004.00396.x |
| [14] | Isaac J H, Lawton D C. Image mispositioning due to dipping TI media: A physical seismic modeling study. Geophysics , 1999, 64(4): 1230-1238. DOI:10.1190/1.1444629 |
| [15] | Claerbout J F. Imaging the Earth's Interior. Oxford: Blackwell Scientific Publications, 1985 . |
| [16] | Ristow D. Migration in transversely isotropic media using implicit finite-difference operators. Journal of Seismic Exploration , 1999, 8(1): 39-55. |
| [17] | Ristow D, Rühl T. 3-D implicit finit-difference migration by multiway splitting. Geophysics , 1997, 62(2): 554-567. DOI:10.1190/1.1444165 |
| [18] | 刘礼农, 高红伟, 刘洪, 等. 三维VTI 介质中波动方程深度偏移的最优分裂Fourier方法. 地球物理学报 , 2005, 48(2): 406–414. Liu L N, Gao H W, Liu H, et al. Wave equation depth migration with optimum split-step Fourier method in 3-D VTI media. Chinese J. Geophys. (in Chinese) , 2005, 48(2): 406-414. DOI:10.1002/cjg2.v48.2 |
| [19] | Shan G. Optimized implicit finite-difference migration for TTI media. 77th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2007. 2290-2294 |
| [20] | Tsvankin I. P-wave signatures and notation for transversely isotropic media: An overview. Geophysics , 1996, 61(2): 467-483. DOI:10.1190/1.1443974 |
| [21] | Shan G. Optimized implicit finite-difference migration for VTI media. 76th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2006. 2367-2371 |
| [22] | Bjrck . Numerical Methods for Least Squares Problems. Philadelphia: SIAM, 1996 . |
| [23] | Kumar D, Sen M K, Ferguson R J. Traveltime calculation and prestack depth migration in tilted transversely isotropic media. Geophysics , 2004, 69(1): 37-44. DOI:10.1190/1.1649373 |
| [24] | 郝 奇. 三维TTI介质地震波场正演模拟. 长春: 吉林大学, 2007. Hao Q. Seismic wavefield modeling in 3D TTI media [Master's thesis] (in Chinese). Changchun: Jilin University, 2007 http://cdmd.cnki.com.cn/Article/CDMD-10183-2007095087.htm |
| [25] | Perelberg A, Hornbostel S. Applications of seismic polarization analysis. Geophysics , 1994, 59(1): 119-130. DOI:10.1190/1.1443522 |
| [26] | Ashyralyev A. On modified Crank-Nicholson difference schemes for Stochastic parabolic equation. Numerical Functional Analysis and Optimization , 2008, 29(3-4): 268-282. DOI:10.1080/01630560801998138 |