2011, Vol. 54
目前为止，Kirchhoff叠前时间偏移仍然是我国石油工业界广泛采用的叠前时间偏移方法.计算效率高以及适应观测系统能力强是Kirchhoff偏移方法的优点.然而，由于使用均方根速度场，传统的Kirchhoff偏移方法适应速度横向变化的能力较弱.刘洪等提出了适用于横向变速情况下的非对称走时计算方法[1， 2]，提高了走时计算的精度从而改进了Kirchhoff偏移的聚焦性能.然而在横向变速介质中，基于常速介质中的格林函数计算得到的振幅加权因子通常是不精确的.相比而言波动方程波场外推算子对振幅的处理更加自然.Dubrulle[3]推广了Gazdag[4]的零偏移距相移偏移方法至共偏移距道集，提出了共偏移距道集的频率-波数域波动方程叠前时间偏移方法.其本质就是对不同方向的平面波分量进行不同的相移.但该方法在频率-波数域实现，无法适应速度的横向变化，Ekren 和Ursin[5]简化了Dubrulle[3]的共偏移距偏移方法，采用叠加速度作为偏移速度场.通过对偏移前后的数据分别实施几何扩散校正，实现了振幅保真的二维频率-波数域共偏移距偏移.Mosher[6]提出了中点-偏移距域内的二维偏移距平面波叠前偏移方法，并讨论了叠前偏移对AVO 的影响.Mosher[7]给出了偏移距平面波深度偏移方法，对每个共射线参数剖面的偏移类似于叠后零偏移距偏移，只是偏移方程中多了一项与地层倾角有关的偏移距倾角项.Pestana[8]提出了一个二维平面波叠前偏移方法，首先对炮集做倾斜叠加，然后分选至共射线参数剖面.对于共射线参数剖面的偏移类似于叠后偏移，只是多了一项有关于射线参数的校正.王华忠等[9]发展了二维偏移距平面波有限差分法叠前时间偏移，通过入射角与偏移距波数的关系，得到每个偏移距波数对应的平面波方程，进而差分求解该方程，并把该方法用于叠前实际数据的处理.
本文把偏移距平面波叠前时间偏移推广到三维空间中，导出了三维偏移距平面波方程，并从射线理论的角度给出了其几何解释，同时指出该方法的假设前提，并用局部倾斜叠加方法来消除对CMP 道集进行离散Radon变换(τ-p变换)引入的噪音.针对实际三维数据的不规则性(中点-偏移距域内方位角展布不均匀及偏移距采样不规则)，提出了与方位角无关的三维倾斜叠加方法，可以得到高信噪比的偏移距平面波数据体；最后我们给出了三维偏移距平面波方程的有限差分解法，并用实际数据的处理结果证明了本文方法的有效性.
2 三维偏移距平面波方程在垂向双程走时(τ)域，双平方根方程(Claerbout[10])有如下频散关系:
|
(1) |
其中，km= (kmx，kmy)为中点波数，kh= (khx，khy)为偏移距波数，kτ 为截距时间波数.ω 为角频率.公式(1)可以整理为
|
(2) |
偏移距波数kh与偏移距射线参数ph有如下关系:
|
(3) |
将(3)式代入(2)式，有
|
(4) |
公式(4)右端的最后一项称偏移距倾角项.在小倾角或者小的射线参数时该项很小(Mosher[7]).忽略偏移距倾角项可得偏移距平面波的频散方程
|
(5) |
在频率-空间域，上述频散方程可以写为
|
(6) |
公式(6)描述了在横向不变速的介质中偏移距平面波的传播(王华忠等[9]).
我们可以从射线理论的角度来解释偏移距平面波的几何含义.在地下局部反射点处，入射角γ 与偏移距波数kh有如下关系(Weglein[11]):
|
(7) |
(1) 式可等价地写为
|
(8) |
由(7)式，可定义如下关系
|
(9) |
将式(9)代入(8)式，整理得到
|
(10) |
在频率-空间域，上述频散方程可以写为
|
(11) |
但是，上式含有平面波在地下局部反射点处的入射角，因而不能直接用来进行实际计算.在地下局部反射点处，入射角γ、反射界面倾角α 与偏移距射线参数ph满足一定的几何关系.如图 1所示，在二维情况下，炮点、检波点与反射点共面.偏移距射线参数ph可以表示为
|
图 1 地下局部反射点处的人射、反射射线示意图 Fig. 1 A sketch of incident and reflected ray on local reflector |
|
(12) |
在三维情形下，上式中的倾角α 变为视倾角αa, 它与反射界面的真倾角θ 和方位角φ (地下反射界面的倾向与炮检连线确定的视倾向之间的夹角)的关系为tanαa =tanθcosφ.为讨论方便，假定炮点、反射点和检波点连线共面.在此面内，(12)式依然成立.偏移距射线参数ph用真倾角θ 和方位角φ 来可以表示为
|
(13) |
显然，当视倾角方向和真倾角方向的夹角等于零时，(13)式退化为(12)式.
进一步地，若忽略反射界面的倾角，(12)式中偏移距射线参数ph可以近似表示为
|
(14) |
结合(11)、(14)式，可以得到偏移距平面波方程(6).显然，只有在横向不变速的介质中偏移距平面波方程才能精确成立.
3 三维偏移距平面波分解 3.1 局部倾斜叠加偏移距平面波分解可以通过倾斜叠加CMP 道集来实现.倾斜叠加(τ-p变换或线性Radon变换)可以在炮检域或中点-偏移距域中实施.Tieman[12]指出，由于CMP道集中同相轴的可预测性，对CMP道集做倾斜叠加可以获得更好的效果.由于受到实际数据的有限采集孔径以及离散采样等因素影响，离散线性Radon变换结果受噪音干扰严重.许多学者研究过如何压制倾斜叠加过程中引入的噪音(Schultz和Claerbout[13]；Stoffa[14]；Kostov和Biondi[15]；Mitchell和Kelamis[16]).其中，较为简单有效的是Mitchell 和Kelamis[16] 提出的双曲速度滤波(Hyperbolicvelocityfiltering)，只要仔细选择叠加速度就可以较好地压制倾斜叠加引入的噪音.Mosher[6]指出，可以在最小二乘意义下求解τ-p域的平面波数据体，但求解过程较为复杂.
为了压制倾斜叠加引入的噪音，一个比较直观的方式就是在局部范围内应用倾斜叠加.类似的思想在高斯射线束偏移(Gaussian Beam Migration)中被广泛应用.在GBM 中，通过局部τ-p变换将高斯窗内的地震道变换到射线参数-双程走时域.本文提出了局部倾斜叠加(Local Slant Stacking)方法.对CMP道集的偏移距平面波分解，可以通过局部倾斜叠加方法来实现.它可以显著减少倾斜叠加引入的噪音.与GBM 不同的是，偏移距平面波分解生成的是偏移距射线参数-截距时间域的平面波数据.
下面以二维情况为例，讨论局部倾斜叠加的思想.CMP道集的时距关系为
|
(15) |
其中，t0 为垂向双程走时，vt0 为均方根速度.
偏移距射线参数定义为
|
(16) |
地震波走时t、截距时间τ 和偏移距射线参数ph有如下关系:
|
(17) |
由(15)，(16)式，h可以表示为偏移距射线参数ph、截距时间τ 和均方根速度vt0 的函数:
|
(18) |
公式(18)事实上指明了中点-偏移距域内的CMP道集与τ-p域中的平面波的对应关系.如图 2所示，我们可以将CMP道集中以(t，h)为中心的一段反射同相轴映射到τ-p域中、以(τ=t-ph，p=dt/dh)为中心的一个邻域内.在实际应用中，考虑到均方根速度总是有一定误差，因而公式(18)计算出的偏移距h必然在一个局部范围内.我们可以在h周围选取一个矩形窗.矩形窗的半宽度可以表示为
|
图 2 局部倾斜叠加方法示意图（Burnet[17]).局部倾斜 叠加将CMP道集中的双曲时距关系（a)映射到了 τ-p 域的椭圆时距关系（b) Fig. 2 A sketch of local slant-stacking approacl (Burnet[17]). The hyperbolic event in a CMP gather (a^ is mapped into an elliptic event (b) in theτ-p domain |
|
W与截距时间τ、射线参数ph和速度误差∣Δv∣有关.显然，局部倾斜叠加的实施与均方根速度有关，但对速度并不敏感.均方根速度影响矩形窗的位置及宽度，因而控制着偏移距平面波数据的信噪比.若取矩形窗为整个空间孔径，则局部倾斜叠加退化为离散线性Radon变换.
3.2 与方位角无关的三维倾斜叠加理论上来讲，二维局部倾斜叠加可直接推广到三维情况.此时，CMP 道集中的绕射双曲面被映射到τ-p域中的旋转椭球体.但对于实际三维地震数据，即使采用了局部倾斜叠加，在τ-p域仍然受假频干扰严重.对于野外炮集记录，炮点的空间采样不足可以导致合成炮集平面波受假频干扰严重.即使对于陆上宽方位角采集的数据，在中点-偏移距域内同样表现为方位角展布及偏移距采样不均匀.中点-偏移距域中数据的不规则性必然导致三维倾斜叠加后假频严重，在射线参数方向严重欠采样.这样得到的平面波数据是不能用来偏移的.
考虑到在偏移距平面波方程的推导过程中假设地下介质是水平层状的，此时方位角的影响可忽略.本文提出了与方位角无关的三维倾斜叠加(azimuth-independent 3D slant-stacking)，将矢量平面波数据沿方位角叠加(公式(19))，生成标量平面波数据体.
其中，Ψ(ph；x)为矢量偏移距平面波数据体，Φ(ph；x)为与方位角无关的偏移距平面波数据体.
|
(19) |
事实上，公式(19)仅是与方位角无关的三维倾斜叠加的一个形式表达.在实施过程中，并没有生成Ψ(ph；x).考虑在局部倾斜叠加过程中，若地震道在以公式(18)为中心的空间窗W内，可以直接映射为由射线参数的模ph所标识的偏移距平面波，而不考虑∣ph∣的方向.因而可以在局部倾斜叠加的过程中直接实现与方位角无关的三维倾斜叠加.这样做可以产生与二维平面波分解相似的数据体.同时，在方位角叠加过程中大大减少了偏移距不均匀导致的空间假频，从而提高了偏移距平面波剖面的信噪比.由于偏移距平面波方程的推导忽略了地层倾角的影响，这样的近似是合理的.当然，三维平面波传播是与方位角有关的，方程(6)中的偏移距射线参数矢量也表明理论上如此.实际过程中，也可以考虑ph的方向，即在不同的方位角范围内进行(19)式定义的方位角叠加运算.
4 三维偏移距平面波方程的有限差分解对于速度横向缓变速的介质，偏移距平面波方程(6)可以近似成立.对其取单平方根近似，有
|
(20) |
对方程(20)的右端项做连分式展开(Claerbout[10])，可以得到如下方程组
|
(21) |
|
(22) |
方程(22)中的系数可以用各种优化方法来确定.用隐式有限差分法求解(22)式，可得如下矩阵方程
|
其中，
|
Ui，jn= $\tilde{U}$(iΔx，jΔy，nΔτ；ω).关于ε 的取值，Claerbout[10]有较详细的讨论.一般地，取ε =$\frac{1}{6}$或ε=$\frac{1}{4}-\frac{1}{{{\pi }^{2}}}$.
平方根算子连分式展开以及拉普拉斯算子双向分裂引入的误差可以用类似Li[18]定义的误差校正算子(附录A)来校正.通过串联求解方程(21)和(23)，可以实现三维偏移距平面波有限差分叠前时间偏移.
对于每个射线参数标识的偏移距平面波的偏移，其计算量与频率-空间域波动方程叠后有限差分偏移相当.总的计算量取决于平面波剖面的个数.三维情况下，对约100~200个的偏移距平面波剖面的偏移就可以得到较好的结果.因此，该方法的计算效率非常高.
5 数值试验 5.1 三维偏移距平面波分解图 3为陆上窄方位角采集的某CMP 的炮检点坐标分布图.很明显，即使是在炮集内数据采样规则，在中点-偏移距域内，CMP 道集仍然表现为方位角展布不均匀及偏移距采样不规则.图 4(a, b)分别为该CMP道集和对应的平面波数据体.偏移距平面波射线参数范围为0~500μs/m.由于CMP道集中缺失某些偏移距的地震道，在τ-p域中相应的位置也出现了空白.由于采用了局部倾斜叠加，平面波数据的信噪比有了显著提高.图 5为陆上窄方位角采集的某CMP 的炮检点坐标分布图，其方位角展布比图 4中略宽.图 6(a, b)分别为该CMP道集和对应的平面波数据体.偏移距平面波射线参数范围为0~400μs/m.可以明显看出，在方位角叠加的过程中克服了τ-p域假频问题.图 7 为陆上宽方位角采集的某CMP的炮检点坐标分布图.图 8(a, b)分别为该CMP 道集和对应的平面波数据体.偏移距平面波射线参数范围为0~400μs/m.该地区地下构造较为平缓，用本文中的平面波分解方法得到了高质量的平面波数据体.
|
图 3 窄方位角观测的CMP面元的炮检坐标分布 Fig. 3 Shot-receiver coordinates of a CMP gather from a narrow-azimuth geometry |
|
图 4 CMP道集（a)及对应的偏移距平面波数据（b) Fig. 4 CMP gather (a) and the corresponding offset plane-wave dataset (b) |
|
图 5 某窄方位角观测的CMP面元的炮检坐标分布 Fig. 5 Shot-receiver coordinates of a CMP gather from a narrow-azimuth geometry |
|
图 6 CMP道集（a)及对应的平面波数据（b) Fig. 6 CMP gather (a) and the corresponding offset plane-wave dataset (b) |
|
图 7 某宽方位角观测的CMP面元的炮检坐标分布 Fig. 7 Shot-receiver coordinates of a CMP gather from a wide-azimuth geometry |
|
图 8 CMP道集（a)及对应的平面波数据（b) Fig. 8 CMP gather (a) and the corresponding offset plane-wave dataset (b) |
用本文中的方法对不同探区的CMP 道集进行平面波分解，均可得到较好的平面波数据体.充分证明了文中方法的有效性.
5.2 三维偏移距平面波偏移实际三维数据来自中国南方某探区.该工区的满覆盖面积为54.17km2，数据大小为66.4G，总共594万道，记录时间为6s, 采样间隔为2ms.该工区的特点是有负向背斜构造，包含比较大的陡倾角构造.首先对该数据进行平面波分解.平面波分解的参数为:射线参数总数Np=160，最小射线参数Pmin=2.5μs/m, 最大射线参数Pmax=400μs/m.图 9 为某条线对应的偏移距平面波剖面，平面波参数为25μs/m.平面波分解后的共射线参数剖面与共偏移距剖面十分相似.分别用三维Kirchhoff PSTM(OMEGA软件)和三维偏移距平面波PSTM进行偏移(使用160个进程)，Omega处理系统的Kirchhoff积分偏移耗时23h, 本文方法偏移耗时18h.所用机群环境软硬件参数如表 1和表 2所示.平面波偏移的效率相对于Kirchhoff偏移有22% 的提升.值的注意的是，平面波偏移的计算量并不是随着射线参数个数的增加而线性增长的.事实上，较大的射线参数对应的平面波数据体只对浅部的成像有贡献，因此偏移外推深度可以随着射线参数的增大而递减.
|
图 9 三维偏移距平面波剖面 Fig. 9 An offset plane-wave section |
|
表 1 计算机群系统软硬件环境 Table 1 Hardware and software environments |
|
|
表 2 两种偏移方法处理相同数据的效率对比 Table 2 A comparison of the computational time on the same field data set between Kirchhoff and offset plane-wave PSTM |
KirchhoffPSTM (OMEGA 软件)和平面波PSTM 偏移剖面分别如图 10和图 11所示.对比两图可以看出，平面波偏移结果较Kirchhoff积分法的结果在局部区域有明显改善.平面波偏移对绕射波的收敛更加干净(图 11 左侧圆圈区域内)，对断层的刻画更加清晰.
|
图 10 Khxhhoff积分法偏移剖面 Fig. 10 A migrated section produced by Kirchhoff PSTM |
|
图 11 三维偏移距平面波偏移剖面 Fig. 11 A migrated section produced by offset plane-wave PSTM |
图 12a 为Kirchhoff偏移输出的共成像点道集，图 12b为偏移距平面波偏移输出的偏移距射线参数域共成像点道集(射线参数范围是0~200μs/m, 经过切除处理).对比两图可以看出，平面波输出的共成像点道集的波形特征更好，且信噪比更高.图 11由40个共射线参数剖面的叠加(射线参数个数Np=40，最小射线参数Pmin=5μs/m, 最大射线参数Pmax=200μs/m)而得到.仅采用40个平面波偏移剖面进行切除叠加(约相当于表 2 中平面波偏移计算时间的50%)，就能达到与Kirchhoff偏移质量相仿的剖面，进一步表明了平面波偏移在计算效率上的优势.在此基础之上，基于成像道集剩余曲率的剩余速度分析方法(Feng[19]，Al-Yahya[20])可以用来更新偏移所用的时间域层速度场.
|
图 12 Khxhhoff偏移（a)与平面波偏移输出的成像道集（b) Fig. 12 Common-image gathers produced by Kirchhoff PSTM (a) and OPW PSTM (b) |
本文提出的局部斜叠加方法可以有效压制离散线性Radon变换过程所引入的线性噪音，生成高质量的偏移距平面波数据体.与方位角无关的三维偏移距平面波分解，是处理三维实际数据的一种折衷的做法，在横向变速介质中会对成像精度产生一定的影响.由于偏移距平面波方程的推导过程中已经忽略了地层倾角的影响，这种折衷在时间偏移的框架下是可以接受的.在数值实验中，对不同偏移距/方位角分布的CMP道集的平面波分解结果也证明了文中方法的稳健性.
偏移距平面波时间偏移采用有限差分法实现波场外推，使用时间域层速度场，可以适应缓横向变速介质.基于波动方程的波场外推算子对振幅的处理也更加可靠.通过与方位角无关的三维偏移距平面波分解，仅需百余个左右的平面波数据体就可以实现偏移，保证了偏移距平面波时间偏移具有很高的计算效率.同时，三维偏移距平面波偏移可以输出偏移距射线参数域共成像点道集.偏移速度场可以通过基于共成像点道集的剩余速度分析方法来更新.偏移距平面波偏移可作为水平地表三维叠前时间偏移的一个很好的解决方案.此外，偏移距平面波的频散方程还可以推广到黏声介质中，从而实现黏声介质中的偏移距平面波叠前时间偏移.
附录A我们可以定义类似Li[18]的误差校正算子:
|
(A1) |
其中，Sx=
|
(A2) |
将(A2)式中的二阶微分算子
|
(A3) |
在频率-波数域上式可以写为
|
(A4) |
其中，
经过误差校正后的波场为
|
(A5) |
在频率-波数域中，对每一层的外推波场
感谢南京石油物探技术研究院提供实际数据及软硬件测试平台.
| [1] | 刘洪, 刘国锋, 李博, 等. 基于横向导数的走时计算方法及其在叠前时间偏移中的应用. 石油物探 , 2009, 48(1): 3–10. Liu H, Liu G F, Li B, et al. The travel time calculation method via lateral derivative to velocity and its application in pre-stack time migration. Geophysical Prospecting for Petroleum (in Chinese) (in Chinese) , 2009, 48(1): 3-10. |
| [2] | 刘国峰, 刘洪, 李博, 等. 山地地震资料叠前时间偏移方法及其GPU实现. 地球物理学报 , 2009, 52(12): 3101–3108. Liu G F, Liu H, Li B, et al. Method of prestack time migration of seismic data of mountainous regions and its GPU implementation. Chinese J. Geophys. (in Chinese) , 2009, 52(12): 3101-3108. |
| [3] | Dubrulle A A. Numerical methods for the migration of constant-offset sections in homogeneous and horizontally layered media. Geophysics , 1983, 48(9): 1195-1203. DOI:10.1190/1.1441542 |
| [4] | Gazdag J. Wave-equation migration with the phase-shift method. Geophysics , 1978, 43(7): 1342-1351. DOI:10.1190/1.1440899 |
| [5] | Ekren B O, Ursin B. True-amplitude frequency-wavenumber constant-offset migration. Geophysics , 1999, 64(3): 915-924. DOI:10.1190/1.1444599 |
| [6] | Mosher C C, Keho T H, Weglein A B, et al. The impact of migration on AVO. Geophysics , 1996, 61(6): 1603-1615. DOI:10.1190/1.1444079 |
| [7] | Mosher C C, Foster D J, Hassanzadeh S. Seismic imaging with offset plane waves. In: Mathematical Methods in Geophysical Imaging IV. SPIE Proceedings Series,Vol. 1996, 2822: 52~63 |
| [8] | Pestana R, Stoffa P L. Plane wave pre-stack time migration. 70th Annual International Meeting, SEG, Expanded Abstracts , 2000: 810-813. |
| [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) (in Chinese) , 2009, 48(1): 11-19. |
| [10] | Claerbout J F. Imaging the Earth's Interior. London: Blackwell Scientific Publications , 1985. |
| [11] | Weglein A B, Stolt R H. Migration-inversion revisited. The Leading-Edge , 1999, 18(8): 950-975. DOI:10.1190/1.1438416 |
| [12] | Tieman H J. Improving plane-wave decomposition and migration. Geophysics , 1997, 62(1): 195-205. DOI:10.1190/1.1444119 |
| [13] | Schultz P S, Claerbout J F. Velocity estimation and downward continuation by wavefront synthesis. Geophysics , 1978, 43(4): 691-714. DOI:10.1190/1.1440847 |
| [14] | Stoffa P L, Buhl P, Diebold J B, et al. Direct mapping of seismic data to the domain of intercept time and ray parameter—A plane-wave decomposition. Geophysics , 1981, 46(3): 255-267. DOI:10.1190/1.1441197 |
| [15] | Kostov C, Biondi B. Improved resolution of slant stacks using beam stacks. 57th Annual International Meeting, SEG, Expanded Abstracts , 1987: 792-794. |
| [16] | Mitchell A R, Kelamis P G. Efficient τ-p hyperbolic velocity filtering. Geophysics , 1990, 55(5): 619-625. DOI:10.1190/1.1442873 |
| [17] | Burnett W, Fomel S. A Gaussian beam analysis of the Radon transform. 78th Annual International Meeting, SEG, Expanded Abstracts , 2008: 2993-2997. |
| [18] | Li Z M. Compensating finite difference errors in 3-D migration and modeling. Geophysics , 1991, 56(10): 1650-1660. DOI:10.1190/1.1442975 |
| [19] | Feng B, Wang H Z, Liu S Y. 3D offset plane-wave finite-difference pre-stack time migration. 79th Annual International Meeting, SEG, Expanded Abstracts , 2009: 2934-2938. |
| [20] | Al-Yahya K. Velocity analysis by iterative profile migration. Geophysics , 1989, 54(6): 718-729. DOI:10.1190/1.1442699 |