2011, Vol. 54
由于各向异性情况下声波方程的渐进解已变的很复杂，基于渐进解的积分类解法在计算效率上的优势已不很明显.基于单程波方程的波动方程偏移在各向异性介质偏移上有了更大的应用前景[1~3].各向同性介质波动方程深度偏移算子研究中，双域法是一类有效的单程波求解思路，已发展的偏移方法如分裂步傅里叶法[4]、傅里叶有限差分法[5]、广义屏法[6， 7]，以及最优分裂傅里叶偏移(OSP)方法[8~11]等均为双域法.已有学者尝试将这类方法应用至各向异性介质的波动方程深度偏移算子研究中，如广义屏方法[7]、宽角屏方法[12]及最优分裂傅里叶方法等[10， 11].但广义屏方法中对各向异性频散关系的近似使得这一偏移算子的应用只能局限于弱各向异性介质，宽角屏方法引入的有限差分导致这一算子在扩展到三维时产生与差分类方法[5， 13]同样的问题，即数值各向异性.
一般而言，具有垂直对称轴的横向各向同性(VTI)介质是对地质构造各向异性的一个很好的近似.本文一方面进一步应用双域法的研究思路来处理横向非均匀的三维VTI各向异性介质.即将横向非均匀各向异性介质分解为横向均匀的速度及Thomsen 各向异性参数[1]场和速度及Thomsen 各向异性参数扰动，在波数域中考虑横向均匀的速度及Thomsen 各向异性参数场，而在空间域基于速度及Thomsen 各向异性参数扰动对波数域的结果进行修正.另一方面，本文着重三维VTI介质中基于单程波算法的波动方程偏移方法在实际工业应用中遇到的两个实际问题开展研究.一是各向异性波动方程深度偏移方法中，除对单程波算子进行研究，发展具有更高效率和精度的单程波算法外，在深度偏移实现上如何提高其计算效率.二是当地下构造沿inline及crossline方向均存在较陡角度时，波动方程偏移方法如何依据现行观测的crossline方向采样较稀疏的三维地震数据对其准确成像.针对这两个问题，本文提出了双域法各向异性波动方程深度偏移方法工业应用中如何提升计算效率及实现陡倾角构造成像的策略.这也是作者以往各向同性波动方程深度偏移研究成果[14， 15]针对各向异性介质的发展.
2 三维VTI各向异性介质中OSP 波动方程叠前深度偏移方法各向异性介质偏移方法的研究首先必须确定使用何种宏观速度模型.VTI各向异性总共有5个参数ε，δ，γ 和Vp0、VS0，较弱的各向异性可忽略准S波速度对准P波传播的影响，因此对本文研究的基于P波成像的VTI各向异性介质偏移，需1个垂直方向的P波速度Vp0，2个各向异性参数ε，δ.相应所应用的宏观速度模型是介质沿垂直方向的射线速度和Thomsen 各向异性参数.若各向异性参数为零，这一方法即退化为各向同性介质中的偏移.
Tsvankin定义了VTI介质中由Thomsen参数表示的相速度表达式V2(θ)[2]:
|
(1) |
其中f=1-Vs02/VP02，θ 为相角，有
|
(2) |
式中ω 为频率，kx，ky，kz分别为x，y，z方向波数.式(1)可以适应于强各向异性情况，但此时必须考虑S波速度的影响，由于一般很难获得VS0，我们在VS0未知的情况下假定f=0.5.由式(2)可得
|
(3) |
将式(2)和(3)代入式(1)得:
|
(4) |
式中
|
与准P波有关的两个根是
|
(5) |
它们分别对应上行和下行的准P波.
式(5)中具有正或零的实部的解kz即对应下行的准P波，由kz可得到对应的相移算子[16]
|
(6) |
式中，j是虚数单位，Δz为波场传播层的厚度.波场反传(即深度延拓)算子可近似得到，为
|
(7) |
利用式(6)和(7)的相移算子，即可得到VTI各向异性介质偏移的相移法，但这一方法仅能适应于横向均匀的层状介质.这里，我们借用双域法的研究思路来处理横向非均匀的各向异性介质.即将横向非均匀各向异性介质分解为横向均匀的速度及Thomsen各向异性参数场和速度及Thomsen 各向异性参数扰动，在波数域中考虑横向均匀的速度及Thomsen各向异性参数场，而在空间域基于速度及Thomsen各向异性参数扰动对波数域的结果进行修正.
定义各向异性介质中分裂的双域单程波算子如下:
|
(8) |
|
(9) |
式中z 是该层介质的中点深度，VP00(z)是该层介质横向均匀的背景速度，VP0(x，y，z)每点的实际速度，kz0(ω，kx，ky)是由背景速度决定的垂直波数；u=1/VP02(x，y，z)-1/VP00(z)是速度扰动，横向均匀时它为零；r1 和r2 是针对各向异性介质引入的各向异性参数扰动，定义为r1 =ε(x，y，z)-ε0(z)，r2 =δ(x，y，z)-δ0(z)，ε0(z)和δ0(z)是该层介质横向均匀的背景各向异性参数；横向均匀的背景速度、背景各向异性参数可由速度和各向异性参数模型平均而得(速度项需采用慢度平均)，ai就是待定的算子系数，我们通过求取ai(i=1，2，3，4)，实现式(8)和(9)对实际单程波算子的最佳逼近.
系数ai(i=1，2，3，4)可通过令
|
(10) |
最小.式中s(ω，kx，ky)是人为定义的权因子，xn，m和yn，m是离散的水平坐标值，(kz)n，m是根据(xn，m，yn，m，z)处的速度和各向异性参数由式(5)解得的垂直波数(实部为正)；区间Ω 被定义如下:
|
(11) |
这里θmaxx 和θmaxx 是事先定义的沿inline(x-)和crossline(y-)方向的最大成像角度，V(θmaxx)和V(θmaxx)是将θmaxx 和θmaxx 代入式(1)得到的相速度.
令
图 1给出了在文献[17]中给出的典型的弱(ε=0.11，δ =0.07)与强(ε=0.36，δ=0.229)各向异性介质情况下准P 波的频散关系，也即式(4)所示方程的解，图中kr= 
|
图 1 （a)与强（b)各向异性介质情况下准P波的频散关系 图中实线与虚线分别为kz的实部与虚部. Fig. 1 Dispersion relations of qP-wave for the weak (a) and strong (b) anisotropy The solid and dashed lines denote the real and imaginary parts of kz , respectively. |
|
图 2 弱与强各向异性介质情况下相移算子与OSP实部（a)与虚部（b)的比较 Fig. 2 Comparisons between the optimum split-step Fourier one-way operator (OSP) and the exact anisotropic phase-shift operator for the weak and strong anisotropy (a) The real parts of the operators； (b) The imaginary parts of the operators. |
求得系数ai(i=1，2，3，4)后，沿用我们在各向同性介质情况下发展的OSP方法[8]，即可得到VTI各向异性介质中波场延拓算法:
|
(12) |
式中，P± (x，y，z+Δz)为单程波波场，± 分别代表向下和向上深度延拓，R2 和R-2分别代表二维正反傅里叶变换.利用5个二维傅里叶正变换和1个二维傅里叶反变换，即可实现三维VTI各向异性介质中波场的深度延拓.
为保证波场延拓的稳定性，需对(12)式的延拓结果做稳定化处理，即
|
(13) |
式中，σ 是一个小实数.
3 三维VTI介质中波动方程叠前深度偏移的变步长深度延拓方法实际工业应用中，上节发展的VTI各向异性介质最优分裂傅里叶偏移方法是基于炮域偏移进行的.炮域波动方程叠前深度偏移，其计算量与ns×nω×nz×NxNylgNxNy成正比.其中，ns为叠前数据炮集数目，nω 为参与波场深度延拓的频率个数，nz为波场深度延拓次数，NxNy为共中心点个数.因而，已发展的多项提高计算效率的措施均是针对上式中各个环节开展的.如为减少炮集个数而发展的相位编码及面炮技术；在共炮集及频率循环中开展并行计算以及针对波动方程深度偏移流程的优化频率采样及在波场深度延拓中采用较大步长[18， 19]以减少深度延拓次数等.
当采用较大步长进行波场深度延拓时，受两个因素的直接影响，一是所用单程波算子在介质横向非均匀时能够适应多大步长；二是大步长深度延拓不能直接给出所有等间距成像点上的地震波场；而深度偏移成像(尤其是陡倾角构造的成像)希望在深度方向每一采样点均进行成像，因而需同时发展相应的波场插值方法以得到等间距成像点上的波场，进而应用成像条件成像.本节将我们已发展的各向同性介质中频率相关的变步长波场深度延拓方法及相应单点插值技术拓展到各向异性介质.
引入新的常量Z，令Δz=Z/ω，将Δz代入(8)式，可得
|
(14) |
修改(4)式可得
|
(15) |
式中
|
上式中的
|
(16) |
此时求解的ai(i=1，2，3，4)将使得频率相关的波场深度延拓步长取得最好的精度；而为使频率相关的波场深度延拓步长等于当前频率地震波场波长的1/4，通常取常量Z=0.5π(Vp0)min, 其中(Vp0)min 是最小的背景速度.
若记速度模型中的最小和最大速度等间距选择的一组参考速度为VP0l，其中l= 1，2，…，m.则得到一组频率无关的变量Zl=ωΔz=0.5VP0lπ，频率相关变步长波场深度延拓算法可表示为:
|
(17) |
式中，
|
为使深度延拓对所有频率的波场延拓都不会出现空间假频，选取最靠当前延拓层最小速度的VP0l对应的Zl，kz0=kz0/ω 可由(15)式解得.
变步长深度延拓不能直接给出等间距成像点上所有频率分量的地震波场；而深度偏移成像希望在深度方向每一采样点均进行成像，因此必须应用波场插值方法得到等间距成像点上的所有频率分量的波场，进而应用成像条件成像.对VTI各向异性介质而言，由于波场插值方向是沿对称轴方向进行的，因而插值计算可沿用各向同性介质的插值方法[15].相对于工业界使用的大步长延拓加插值的偏移成像策略[4， 18]和直接在单程波算子中考虑介质纵向变化的大步长单程波算子[19]，这一方法既减少了偏移的计算量(由于较大减少了中、低频分量延拓时的延拓层数，总计算量减少近三分之二)，又保证了成像质量.
4 三维VTI介质中波动方程叠前深度偏移的反假频深度延拓方法现行野外观测的地震数据在crossline方向采样通常是稀疏的.而对基于单程波方程的深度偏移算法而言，直接应用稀疏采样数据偏移将导致大角度传播的波场不能正确归位，从而使得陡角度构造不能正确成像且产生空间假频.针对这一实际应用中的问题，本节建议一个三维VTI各向异性介质中不需插值的陡倾角构造成像策略.即首先通过填充空道，得到满足陡角度构造成像所需的、正确模拟大角度波场传播的(形式上的)理想采样数据，再构造反假频波场延拓算子，进行反假频和压制(偏移)噪音的波场深度延拓；这样即可不通过叠前数据的插值重建而实现陡倾角地层的正确成像.这里的“形式上"是指并不需要真的在偏移处理中输入的地震数据中添加许多零值，而仅在波场延拓时考虑这些空道即可，从而避免了对叠前数据填充空道所产生的巨大存储需要.
在单程波方程的三维波场深度延拓算法中，两个方向的空间采样Δx和Δy决定了波数域的可能取值范围，即
|
(18) |
超出这一范围的波数将在空间域产生假频；因此，空间采样完全决定了波场深度延拓算法可正确模拟的地震波场的最大传播角度(与垂直方向夹角).
实际勘探中，空间采样往往不够密，且在inline方向和crossline方向采样不同，此时，最大传播角度θmax 往往小于90o.对VTI各向异性介质，在xz平面中传播的波场，其最大传播角度θmaxx 与inline方向空间采样Δx的关系如下:
|
(19) |
而对yz平面中传播的波场有
|
(20) |
这里ωmax 是地震数据的最大有效频率.
式(18)所定义的是最大的波数范围，超出这一范围的波数将在波场深度延拓中产生空间假频；同时，我们也可根据成像的需要(为压制噪音)，人为定义波场在xz平面和yz平面中传播的最大传播角度，并据此定义有效的波数域范围.因此，我们修正有效的波数范围Ω 如下:
|
(21) |
由于射线参数在波场的传播过程中是不变的，因此式中用最大射线参数
|
来定义地震波场的最大传播角度.衰减式(21)所定义的波数范围以外的波场，即可实现三维VTI各向异性介质中反假频及控制波场最大传播角度的成像.定义反假频因子D(pmaxx, pmaxy):
|
(22) |
即反假频因子D(pmaxx, pmaxy)在式(21)所定义的波数范围内为1，而在所定义的范围外按指数衰减；α 为定义的衰减系数.在上节发展的三维VTI各向异性介质波动方程深度延拓算法中的波数域计算中引入这一反假频因子，即可避免空间假频且控制参与成像的波场的最大传播角度，得到反假频波场延拓算子.据此，式(17)所示的波场深度延拓公式可修改为:
|
(23) |
为验证本文方法的正确性和工业应用的可能性，分别将其应用至SEG 二维Hess VTI各向异性介质理论模型数据及我国东部某工区野外数据集观测系统下的三维脉冲响应计算，前者证明了本文理论和方法的正确性，后者表明了本文方法初步具有实际工业应用价值.
5.1 SEG二维Hess各向异性VTI模型成像SEG 二维Hess各向异性VTI模型数据共720炮，炮间距100ft, 最大道数656道，道间距40ft, 单边激发，记录长度为7.992s, 采样间隔6 ms, 最大偏移距26200ft, 最小偏移距0ft.速度模型及各向异性参数模型横向和纵向采样间隔均为20ft.该模型设计的成像考察点主要有位于工区左部的盐丘A、位于工区中部的自上而下的三个透镜体B，C，D以及位于工区右部的贯穿整个深度的大断层E.
图 3展示了SEG 二维HessVTI各向异性介质模型宏观模型参数VP0，ε 及δ 剖面，需要指出的是透镜体B 在VP0剖面上不存在，仅存在于各向异性参数剖面上.图 4a给出了SEG二维HessVTI各向异性介质模型波动方程各向异性OSP 方法偏移剖面，作为对比，图 4b 给出了SEG 二维HessVTI各向异性介质模型波动方程各向同性OSP 方法偏移剖面.仔细观察可见，当忽略各向异性时，盐丘A 边界成像模糊，其顶边界左侧翼成像位置不正确(参见图 5)，右侧翼与上覆地层接触关系不清楚、底边界成像丢失.透镜体B 未能成像，导致其相邻的水平地层成像亦受到干扰；透镜体C，D 成像杂乱.工区右侧断层E 的断面波未能收敛，其成像位置在水平方向及深度上均存在误差；整个剖面信噪比较低.在考虑各向异性的剖面上，盐丘顶底边界更为清楚，与上覆地层的关系得到正确刻画；透镜体B，C，D 均得到准确成像；工区右侧断层E 断面波收敛，断点清晰，成像位置正确；整个剖面信噪有了较大幅度的提高.从而表明，本文发展的各向异性VTI介质波动方程叠前深度偏移方法可较好地处理强非均匀的VTI各向异性介质的情况.
|
图 3 SEG二维Hess VTI各向异性介质模型宏观模型参数VP0(a)，$\varepsilon $(b)及$\delta $(c) Fig. 3 SEG 2D Hess VTI anisotropic model： Vp0 (a)$\varepsilon $(b) and $\delta $(c) , respectively |
|
图 4 SEG 2D HeSS VTI各向异性模型波动方程各向异性（a)和各向同性（b)OSP方法偏移剖面 Fig. 4 The anisotropic (a) and isotropic (b) OSP prestack depth migration result of SEG 2D Hess VTI anisotropic datasets |
|
图 5 SEG 2D Hess VTI各向异性模型盐丘左侧翼各向异性OSP(a)与各向同性OSP(b)偏移结果对比 Fig. 5 Comparisons of top left salt boundary migration images of SEG 2D Hess VTI model between the anisotropic (a) and isotropic (b) OSP prestack depth migration results |
受各向异性参数建模技术制约，本文方法的实际资料应用测试存在困难.实际应用中基于炮域的波动方程偏移计算结果是由炮集内所有接收地震道叠前脉冲响应的叠加，因此本节应用本文方法进行实际数据规模的叠前脉冲响应计算以测试其对工业应用的适应能力.将某野外炮集数据集的速度模型的速度取值为各向异性测试中速度VP0，其各向异性Thomsen参数则由一定的规则由VP0计算得到，速度及Thomsen 参数模型的空间采样分别为Δx=25m, Δy=50m, Δz=5m.输入震源子波为一道实际地震道.
图 6a给出这一测试的三维VTI各向异性介质叠前脉冲响应inline及crossline方向切片，图中竖直线标识了输入震源地震道的位置.横直线标识了用于与忽略各向异性比较的同相轴的位置.作为对比，图 6b给出了该测试忽略各向异性时对应位置的三维叠前脉冲响应inline及crossline方向切片，图中竖直线标识了震源子波的位置，曲线则标识了考虑各向异性时脉冲响应同相轴的位置.图 6c进一步展示了这一测试三维VTI各向异性介质叠前脉冲响应与忽略各向异性时三维叠前脉冲响应深度为2km的水平切片对比，图中“+"标识了震源的横向位置.
|
图 6 (a)三维VTI各向异性介质叠前脉冲响应inline方向（Une130)及cr〇ssine(CDP295)方向切片，图中竖直线标识了震源的位置.水平直线标识了各向异性叠前脉冲同相轴的位置；（b)三维各向同性介质叠前脉冲响应inline方向 (line130)及crosline方向（CDP295)切片，图中竖直线标识了震源的位置.曲线标识了图a中各向异性叠前脉冲同相轴的 位置；（)三维VTI各向异性介质叠前脉冲响应（左）与三维各向同性介质叠前脉冲响应（右）深度为2km的水平切片对比，图中“ + ”标识了震源的横向位置;输人震源为一道地震记录 Fig. 6 (a) 3-D VTI prestack migration impulse response of the anisotropic OSP method. The position of the sourced, e. a real seismic record) is indicated by the vertical line. The horizontal lines shown the events employed to comprise with the isotropic OSP method. (b) 3-D prestack migration impulse response of the isotropic OSP method. The position of the source is indicated by the vertical line. The curve lines shown the events obtained by the anisotropic OSP method. (c) Comparisons of horizontal slices at the depth of 2 km between the anisotropic and isotropic OSP 3 D prestack migration impulse response. The horizontal position of the source is indicated by the+ |
综合比较图 6 可见，各向异性叠前脉冲响应与各向同性叠前脉冲响应波形形态相似，在输入震源为实际地震道的情况下二者的信噪比基本处于同等水平.比较图中所标识的同相轴可见，在垂直方向二者走时相同，在水平方向则由于各向异性情况下受Thomsen参数影响，波传播较各向同性情况快.考虑到本文测试的观测系统及测试参数均是应用了某野外数据集的一炮数据的偏移体量，因而可表明本文方法具有较强的适应速度及各向异性Thomsen参数横向变化及实际数据规模偏移计算的能力.
6 结语本文将各向同性复杂构造介质中的最优分裂Fourier方法发展为各向异性介质中的波场延拓算法，并着重就其工业应用中遇到的提高计算效率及反假频成像策略等两个问题展开讨论.文中建议了三维VTI各向异性介质中可提高其计算效率的频率相关变步长波场深度延拓算法，其核心是对波场的高频分量使用较小延拓步长，而对中、低频分量使用较大延拓步长，进而结合波场插值实现理想深度采样的深度成像.文中基于对地层倾角的估计，通过(形式上)填充空道重建理想采样的地震数据并引入以射线参数为基准的反假频因子，可实现三维VTI各向异性介质波动方程深度偏移陡倾角地层的准确成像.SEG 二维HessVTI各向异性介质理论模型数据测试及野外数据集观测系统下的三维脉冲响应计算表明本文提出了一种具有工业应用潜力的三维各向异性波动方程叠前深度方法.
| [1] | Thomesen L. Weak elastic anisotropy. Geophysics , 1986, 51(10): 1954-1966. DOI:10.1190/1.1442051 |
| [2] | Tsvankin I. P-wave signatures and notation for transversely isotropic media: an overview. Geophysics , 1996, 61(2): 467-483. DOI:10.1190/1.1443974 |
| [3] | Tsvankin I, Gaiser J, Grechka V, et al. Seismic anisotropy in exploration and reservoir characterization: An overview. Geophysics , 2010, 73(5): 75A15-75A29. |
| [4] | Stoffa P L, Fokkema J T, de Luna Freire R M, et al. Split-step Fourier migration. Geophysics , 1990, 55(4): 410-421. DOI:10.1190/1.1442850 |
| [5] | Ristow D, Rühl T. Fourier finite-difference migration. Geophysics , 1994, 59(12): 1882-1893. DOI:10.1190/1.1443575 |
| [6] | Le Rousseau J H, De Hoop M V. Modeling and imaging with the scalar generalized-screen algorithms in isotropic media. Geophysics , 2001, 66(5): 1551-1568. DOI:10.1190/1.1487101 |
| [7] | Le Rousseau J H, de Hoop M V. Scalar generalized-screen algorithms in transversely isotropic media with a vertical symmetry axis. Geophysics , 2001, 66(5): 1538-1550. DOI:10.1190/1.1487100 |
| [8] | Liu L N, Zhang J F. 3D wavefield extrapolation with optimum split-step Fourier method. Geophysics , 2006, 71(3): T95-T108. DOI:10.1190/1.2197493 |
| [9] | Zhang J F, Liu L N. Optimum split-step Fourier 3D depth migration: Developments and practical aspects. Geophysics , 2007, 72(3): S167-S175. DOI:10.1190/1.2715658 |
| [10] | 刘礼农, 高红伟, 刘洪, 等. 三维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 3D VTI media. Chinese J. Geophys. (in Chinese) , 2005, 48(2): 406-414. DOI:10.1002/cjg2.v48.2 |
| [11] | Liu L N, Zhang J F. Optimum split-step Fourier one-way operators for seismic modeling and imaging in 3D VTI media. Geophysical Research Letters , 2006, 33: L09308. |
| [12] | Han Q Y, Wu R S. A one-way dual-domain propagator for scalar qP-waves in VTI media. Geophysics , 2005, 70(2): D9-D17. DOI:10.1190/1.1884826 |
| [13] | Fei T W, Liner C L. Hybrid Fourier finite-difference 3D depth migration for anisotropic media. Geophysics , 2008, 73(2): S27-S34. DOI:10.1190/1.2828704 |
| [14] | 刘礼农, 陈树民, 张剑锋. 稀疏采样下陡角度构造的波动方程深度偏移成像. 地球物理学报 , 2006, 49(5): 1452–1459. Liu L N, Chen S M, Zhang J F. Steep dipping structure depth imaging using wave equation based migration schemes on 3-D sparse dataset. Chinese J. Geophys. (in Chinese) , 2006, 49(5): 1452-1459. |
| [15] | 张剑锋, 卢宝坤, 刘礼农. 波动方程深度偏移的频率相关变步长延拓方法. 地球物理学报 , 2008, 51(1): 222–228. Zhang J F, Lu B K, Liu L N. Frequency-dependent varying-step depth extrapolation scheme for wave equation based migration. Chinese J. Geophys. (in Chinese) , 2008, 51(1): 222-228. |
| [16] | Zhang J F, Verschuur D J, Wapenaar C P A. Depth migration of shot records in heterogeneous, transversely isotropic media using optimum explicit operators. Geophys. Prosp., , 49(3): 287-299. DOI:10.1046/j.1365-2478.2001.00255.x |
| [17] | Wang Z J. Seismic anisotropy in sedimentary rocks, part 2: Laboratory data. Geophysics , 2003, 67(5): 1423-1440. |
| [18] | Fu L Y. Wavefield interpolation in the Fourier wavefield extrapolation. Geophysics , 2004, 69(1): 257-264. DOI:10.1190/1.1649393 |
| [19] | 刘洪, 袁江华, 陈景波, 等. 大步长波场深度延拓的理论. 地球物理学报 , 2006, 49(6): 1779–1793. Liu H, Yuan J H, Chen J B, et al. Theory of large step wavefield depth extrapolation. Chinese J. Geophys. (in Chinese) , 2006, 49(6): 1779-1793. |