2010, Vol. 53
2. 中国科学院地质与地球物理研究所, 北京 100029
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
炮点和检波点不在同一水平面上，是陆上地震勘探经常遇到的问题.针对这一类起伏地表采集的地震数据已发展了许多处理方法.这些方法可归纳为两类：一是静校正技术[1]，它假设地震波是沿铅垂线传播的；另一类是基准面重建技术[2]，它正确考虑了地震波在近地表传播的实际路径.静校正技术仅适用于近地表存在明显低速层情况（多数情况这一假设是成立的），而基准面重建方法需要分别对共炮点道集和共检波点道集进行波场延拓，计算成本较高.这两类方法都需要事先确定近地表速度.由于起伏地表采集的地震资料信噪比较低，因此获得准确的近地表速度模型是一个相当困难的任务[3, 4].
为了避免以上方法存在的问题，也发展了直接基于起伏地表采集数据的叠前偏移成像方法，一类是基于浮动基准面的叠前时间偏移方法[1]，另一类是直接考虑起伏地表形状的叠前深度偏移技术[5].浮动基准面叠前时间偏移方法是静校正技术与叠前时间偏移的结合，当有高速层（老地层）出露时，垂直入射和出射的静校正的应用条件已不成立，这一方法将产生较大的误差.起伏地表叠前深度偏移技术需要准确的近地表速度和地下构造的层速度模型，而两类速度又相互影响，与海上或水平地表情况下确定深度偏移的层速度模型相比，速度建模的难度增加了许多；起伏地表叠前深度偏移实际上将问题的难点前推至速度建模，目前尚没有取得好的实际应用效果.
叠前时间偏移可对一类倾角、断层较为复杂，但速度横向变化不是很剧烈的构造进行较好成像.叠前时间偏移的一个主要优点是只需使用均方根速度；这样可简单地通过速度扫描等方式得到恰当的速度模型，回避了使用叠前深度偏移方法面临的速度建模困难.而将叠前时间偏移的这一特点推广、应用于解决起伏地表问题，这是本文研究的出发点；不同于浮动基准面叠前时间偏移方法，本文将在正确考虑地震波在近地表传播的实际路径的基础上展开研究.作者已在文献[6]中发展了一个基于双速度参数的叠前时间偏移方法，较好地将叠前时间偏移技术推广应用于起伏地表问题；但双速度参数将速度扫描扩大到二维，在实际应用中确定速度模型还是不够简单.本文方法则将两个速度参数解耦，只需分别扫描确定两个独立的速度参数，使得速度建模变得更简单.
本文方法是通过将叠前时间偏移与基准面重建技术结合而发展的.不同于基于波场延拓的基准面重建方法，本文在起伏表面上填充了虚拟介质，形成均匀的虚拟层，通过结合瑞利积分和相移法，将炮点和检波点反向延拓重建到虚拟层的顶界面上.通过扫描均匀虚拟介质的速度，获得了最适用于下一步叠前时间偏移处理的基准面重建效果.对重建后的数据采用常规的叠前时间偏移技术，即可较好地解决起伏地表采集数据的偏移成像问题.为减少叠前时间偏移的计算量，本文发展了基于最大成像角度自适应确定偏移孔径的叠前时间偏移算法.为验证本文方法对近地表有较复杂的速度横向变化及高速层出露情况的处理能力，文中给出了应用该方法处理SEG复杂地表模型理论数据[7]的例子；为验证本文方法的广泛适用性，文中也将这一方法应用于近地表存在明显低速层的理论模型数据.
2 混合法虚拟基准面波场重建本文的关键是通过建立虚拟基准面，实现基于均匀介质波场延拓的基准面重建方法；这既简化重建算法，也可以通过速度扫描，获取最佳的虚拟介质速度.针对虚拟的均匀介质中波场反向基准面重建问题，本文发展了结合瑞利积分和相移法的混合法.对虚拟基准面上的炮点，通过瑞利积分获得其在真实炮点上的走时和幅值，形成对应基准面上炮点的合成炮记录；对合成炮记录的所有处于起伏表面的检波点，应用基于波动方程的相移法，将检波点正传到虚拟基准面上，从而实现虚拟基准面的波场重建.这一混合算法所需计算量较小，使得扫描最佳虚拟介质的速度成为可能.下面将分别讨论这两个过程.
2.1 合成炮记录在起伏表面上填充了虚拟介质，可将虚拟介质顶面作为虚拟基准面.对在虚拟基准面上的炮点（偶极源），若将其看作空间脉冲函数，由瑞利积分Ⅱ [8]可得，其在真实炮点上的波场可表达为
|
(1) |
式中（xs，ys，zs）为真实炮点坐标，（x0，y0，z0）为基准面上炮点的坐标，j为单位虚数，f（ω）为基准面上震源的频谱，Δ为空间采样的特征尺寸，c是虚拟介质的速度，
二维情况下，式（1）表达式可变为
|
(2) |
式中H1(2)（ωr/c）为一阶二型汉克尔函数，r=
|
(3) |
在积分法偏移中，一般均应用远场解计算地震波的走时与幅值.由于成像目标一般远离炮点和检波点，这一近似不影响成像效果.但在基准面重建时，就低频分量而言，基准面上的炮点和真实炮点的距离并不够远，因此这种情况下忽略近场解将带来明显的误差.图 1给出了二维情况下远场解与准确解在不同频率下的相对误差，从图中可看出，为保证相对误差小于1%，25、50 Hz和100 Hz对应的距离分别需要大于740、370 m和185 m；频率越低，应用远场解的起始距离就越大.为保证基准面重建的精度，建议直接采用（1）或（2）式计算走时与幅值，而不是像偏移算法那样采用远场解.
|
图 1 远场解与准确解在不同频率下的相对误差 Fig. 1 Relative errors of the far field and exact solutions for different frequencies |
将原始炮记录转换到频率域，记为G（xs，ys，zs；x，y，z，ω）.利用（1）或（2）式的结果，可得基准面上炮点（x0，y0，z0）对应的合成炮记录：
|
(4) |
式中S代表叠加孔径，它是综合考虑成像角度、z0 -zs值以及频率ω决定的；选择恰当的孔径，既可减少合成炮的计算量，也可压制噪音.我们将在第四节讨论偏移孔径时一并讨论叠加孔径的选择.为保持实际信号的波阻特征，在计算式（4）中的P（xs，ys，zs，ω）时，令对应的f（ω）为一常实数.
合成炮记录对应着炮点在虚拟基准面，但检波点（x，y，z）在实际起伏地表上的单炮记录.对这一炮记录进行反向正传，即可得到检波点同样在虚拟基准面上的基准面重建数据.下节将讨论如何实现这一过程.
2.2 基准面波场重建假设在起伏表面上填充了虚拟介质后，起伏地表上的检波点均处于均匀的虚拟介质中.因此可应用相移法实现基准面波场重建.对每一合成炮记录，统计检波点的z坐标，按其大小及容许误差，形成深度系列zn，zn-1，…，z1，z0，其中z0既是基准面的深度坐标.基于相移法[9]和离散的傅氏变换，可得波场延拓的递推公式如下：
|
(5) |
|
(6) |
式中叠加是针对所有属于深度zi-1的检波点.若起伏地表上检波点的平均水平采样间距分别为Δx和Δy，合成炮的最大水平孔径分别是xmax和ymax，则式（5）和（6）中的水平波数分别为
|
(7) |
|
(8) |
为提高计算效率，当式（5）和（6）中右端的离散傅氏变换涉及较多的项叠加时，可应用非规则采样快速傅里叶变换[10]（NFFT）进行计算.
由式（5）可得到
式（5）和（6）的波场延拓算法不仅具有较高的计算效率，还可以在处理检波点处于起伏地表的同时，灵活处理检波点的水平分布不均匀问题；而这是实际起伏地表采集时经常遇到的.因此，本文方法也避免了现行基于波场延拓的基准面重建方法存在的炮点和检波点水平分布不均匀时的困难.
3 速度建模本文算法利用虚拟介质的速度和基于虚拟基准面的均方根速度，实现起伏地表采集数据的叠前成像.虚拟介质的速度用于基准面数据重建；均方根速度用于基准面重建数据的叠前时间偏移.这样，两个速度参数是解耦，可单独确定各个速度参数，从而简化了速度建模过程.
合理的虚拟介质速度应与实际的近地表速度匹配，使得由两者结合形成的虚拟层更趋近均匀介质层.这样，上节基于均匀介质的基准面重建更加合理，而包括虚拟层的实际速度结构也更适宜应用叠前时间偏移技术.如果虚拟介质的速度选择不合理，上节的基准面重建就不能得到合理的结果；这个不合理的结果将表现在合成炮上的同相轴连续性较差.因此，本文通过速度扫描确定虚拟介质的最佳速度：选择几个典型点，对不同速度求得这些点为炮点的合成炮，比较合成炮，使得同相轴光滑连续的速度即为虚拟介质的最佳速度.
下面以SEG逆掩断层模型（图 2）为例说明这一方法.在图 2上选取不同的炮点位置，分别用3000、3500、4000 m/s和4500 m/s进行基准面重建.图 3给出了炮点在15030 m处的实际单炮和不同虚拟介质速度的合成炮记录比较，这里对合成炮仅取了局部.对比图 3中的各个合成炮可知，4000 m/s速度下的合成炮其同相轴连续性更好，且深层反射的同相轴更清晰（与4500 m/s速度对比），因此确定虚拟介质的最佳速度为4000 m/s.后面的数值算例将证明，这一速度确实可取得较好的成像结果.
|
图 2 SEG逆掩断层构造速度模型 Fig. 2 SEG overthrust velocity model including surface topography |
|
图 3 SEG模型的原始炮集与不同虚拟速度下的合成炮对比(炮点:15030 m) (a)原始; (b)速度:3000 m/s (c)速度:3500 m/s (d)速度:4000 m/s (e)速度:4500 m/s. Fig. 3 Comparison of the original l shot gather with the artificial shot gathers under different virtual velocities (shot location:15030 m) |
当完成基准面数据重建，炮点和检波点已被校正到同一水平面上，这时可采用常规的叠前时间偏移技术确定基于这一新的基准面的均方根速度.本文采用如下流程：（1）对几个典型CMP道集作动校正，拾取动校速度，对几个道集的结果作平均，求得一横向均匀的均方根速度场；（2）基于横向均匀的均方根速度进行叠前时间偏移，形成各个CDP位置的CRP道集；（3）对所有CRP道集按横向均匀的均方根速度作反动校，然后作动校分别拾取各CDP位置处的均方根速度，对拾取的均方根速度场作平滑，即得到本文用于叠前时间偏移的均方根速度场.
通过这样一个简单过程，我们就得到了本文方法需要的两类（虚拟介质的速度和基于虚拟基准面的均方根速度）速度模型.
4 算法实现流程为减少叠前时间偏移的计算量并压制偏移噪音，本文发展了基于最大成像角度自适应确定偏移孔径的叠前时间偏移算法.本文成像算法采用输入道成像方式，最大成像角度是由单道数据成像在空间的展开孔径决定的，给定最大成像角度，我们就可确定不同垂直走时对应的成像孔径，以二维问题为例，成像孔径x可由下式决定：
|
(9) |
式中vrms为偏移用的均方根速度，τ为垂直走时，h为半偏移距，θ为定义的最大成像角度.计算合成炮的叠加孔径的选择将与（9）式的偏移孔径相匹配，即落在孔径中的实际地表上的全部炮点，就是合成虚拟地表上该虚拟炮点的叠加孔径.基于（9）式的偏移孔径，将与计算合成炮的孔径相匹配，这既保证了有效信号正确成像，又较好地压制了噪音.由于减少了每个输入道需要计算的区域，也极大地减少了偏移的计算量.
本文方法的计算流程如下：（1）利用第2节和第3节所述的方法确定最适合叠前时间偏移的虚拟层速度，形成炮点和检波点位于同一虚拟界面的合成炮记录；（2）由第3节所述过程对新合成数据进行速度分析，获得叠前时间偏移所需的均方根速度场；（3）对合成炮记录进行时间偏移成像，形成CRP道集，将CRP道集校平叠加即实现了对整个数据的叠前时间偏移.
5 数值算例为验证本文方法处理近地表有较复杂的速度横向变化和高速层出露的能力，文中给出了应用这一方法处理SEG复杂地表模型理论数据的例子；实际上，这一理论模型是为验证叠前深度偏移方法发展的，它假设已获得了准确的速度模型，仅验证叠前深度偏移方法能否在速度模型已知情况下解决复杂起伏地表采集数据的成像问题.从图 2速度模型，特别是近地表速度分布来看，直接由反射数据获得这样的速度模型是很困难的.
为验证本文方法的广泛适用性，文中也将这一方法应用于近地表存在明显低速层的理论模型数据.对这类模型，浮动基准面叠前时间偏移方法也是可以得到较好成像效果的.
对这两类极端问题的实验，将表明本文方法可处理较为复杂的起伏地表成像问题.
5.1 二维断陷盆地模型该断陷盆地模型近地表速度变化较为平缓（图 4），存在一层最大厚度20 m左右，速度为1400 m/s的风化层，风化层以下地层的速度为2200 m/s，近地表构造简单并且速度较低.利用声波方程数值模拟方法[11]获得的数据集共151炮，观测系统不规则，数据最大偏移距6000 m，记录长度4s，时间采样率为2 ms.直达波和面波已被剔除，但表面和层间多次波包含在炮记录中.定义新的基准面在1500 m.
|
图 4 断陷盆地速度和构造模型 Fig. 4 Velocity model of rift-subsidence basin |
首先通过扫描，确定了虚拟介质的最佳速度为2000 m/s.图 5给出了由合成炮抽取的最小偏移距剖面与原始数据的最小偏移距剖面的对比，可以看出合成数据的同相轴连续性较好，较好地反映了速度模型中的主要构造；从同相轴的连续性来看，这一最小偏移距剖面应该优于直接采用静校正方法获得的统一基准面数据的对应剖面.图 6是由基准面合成数据求得的均方根速度模型，对比速度模型图 4可知，叠前时间偏移使用的均方根速度较实际的层速度变化更为平缓.图 7给出了偏移叠加剖面，对比图 4可知，构造得到了正确成像，断点归位正确.由于基准面定义为水平面，避免了采用浮动基准面时的转换流程.图 7中深部成像的一些频散是由于数据正演过程中深部网格剖分较大，炮记录本身存在的频散所致；剖面中的干扰波为原始数据中存在的多次波.
|
图 5 断陷盆地模型最小偏移距剖面对比 (a)原始; (b)合成. Fig. 5 Comparison of the minium-offset sections (a) original; (b) artificial. |
|
图 6 断陷盆地模型的均方根速度场 Fig. 6 Model of rms velocity for the ritt-subsidence basin |
|
图 7 断陷盆地模型叠前时间偏移剖面 Fig. 7 Prestack time migration stacked section of the rift-subsidence basin |
从图 2速度模型中可以看出，地表起伏较为剧烈，最大高差达1237 m，构造复杂，速度的横向及垂向梯度均较大.该模型高速层出露地表，静校正的假设条件已不能满足.数据集共278炮，炮间距90 m，道间距15 m，最大偏移距3600 m，采用中间放炮两边接收的观测系统，记录总长度8 s，采样率4 ms.定义新的基准面在0 m高程.
如第3节中图 3所示那样，首先确定4000 m/s为最佳的虚拟介质速度.然后由合成炮记录求得了如图 8所示的均方根速度.图 9给出了偏移叠加剖面，对比图 2可知，主要构造得到正确成像，复杂表面和逆掩推覆构造下的单倾界面被正确成像；直界面变弯曲是因为成像是在垂直旅行时坐标下的缘故.与参考文献[12, 13]中的结果进行对比，可以看出本文偏移剖面中的同相轴更加连续，信噪比得到提高，这表明本文方法较已有方法更有效；这是因为本文方法在描述地震波的实际传播路径方面更进了一步.由于叠前时间偏移方法本身是不适宜处理横向速度剧烈变化问题的，因此偏移剖面中存在较强的偏移噪音，且部分大倾角地层没能正确成像.实际上，本文方法并不是针对解决这类复杂问题而发展的.
|
图 8 SEG逆掩断层模型的均方根速度场 Fig. 8 Model of rms velocity for SEG overthrust model |
|
图 9 SEG逆掩断层模型叠前时间偏移剖面 Fig. 9 Prestack time migration stacked section for SEG overthrust model |
本文发展了一个针对起伏地表采集数据的叠前时间偏移方法和新流程.该方法用基准面重建取代了基于浮动基准面的静校正技术，因此克服了高速层出露时静校正方法的误差.文中提出了一种结合虚拟界面、瑞利积分和相移法的混合的基准面重建方法，有效地避免了确定近地表速度的困难；本文的基准面重建方法可灵活地处理炮点和检波点水平分布不均匀，具有很高的计算效率.本文方法的基础是叠前时间偏移，因此这一方法仅适用于速度横向中等变化问题.就我国勘探的地质目标而言，这一方法是适宜的.两类较极端问题的理论模型数据成像结果表明，本文方法具有较好的适应性，既可处理近地表存在明显低速层的问题，也可较好地处理近地表速度变化剧烈、存在高速层出露的问题.与各类野外静校正方法一样，这一方法仅能考虑长波长静校正，对实际存在的短波长静校正量，还需配合应用剩余静校正方法，这一点对各类考虑起伏地表的成像方法都是共同的.尽管所给的数值算例是二维，但文中讨论同时考虑了二维和三维问题，将这一方法直接应用于三维问题不存在任何障碍.
| [1] | Yilmaz O. Seismic Data analysis:Processing, Inversion, and Interpretation of seismic data. Society of Exploration Geophysicists, Tulsa, 2001 |
| [2] | Berryhill J R. Wave-equation datuming before stack. Geophysics , 1984, 49: 2064-2066. DOI:10.1190/1.1441620 |
| [3] | Chang X, Liu Y K, Wang H, et al. 3-D tomographic static corrections. Geophysics , 2002, 67: 1275-1285. DOI:10.1190/1.1500390 |
| [4] | Kelamis P G, Erickson K E, Verschuur D J, et al. Velocity-independent redatuming:A new approach to the near-surface problem in land seismic data processing. The Leading Edge , 2002, 21: 730-735. DOI:10.1190/1.1503185 |
| [5] | Reshef M. Depth migration from irregular surface with the depth extrapolation methods. Geophysics , 1991, 56: 119-122. DOI:10.1190/1.1442947 |
| [6] | 董春晖, 张剑锋. 起伏地表下的直接叠前时间偏移. 地球物理学报 , 2009, 52(1): 239–244. Dong C H, Zhang J F. Prestack time migration including surface topography. Chinese J. Geophys. (in Chinese) , 2009, 52(1): 239-244. |
| [7] | Gray S H, Marfurt K J. Migration from topography:Improving the near-surface image. Can.J.Expl.Geophys. , 1995, 31: 18-24. |
| [8] | Berkhout A J. Seismic migration:Imaging of acoustic energy by wavefield extrapolation. Elsevier Science Publication Company, 1980 |
| [9] | Gazdag J. Wave equation migration with the phase-shift method. Geophysics , 1978, 43: 1342-1351. DOI:10.1190/1.1440899 |
| [10] | 熊登, 张剑锋. 表驱动的二维非规则采样快速傅里叶变换. 地球物理学报 , 2008, 51(6): 1860–1867. Xiong D, Zhang J F. Tabel-Driven 2-D non-uniform sampling fast Fourier transform. Chinese J. Geophys. (in Chinese) , 2008, 51(6): 1860-1867. |
| [11] | 徐义, 张剑锋. 地震波数值模拟的非规则网格PML吸收边界. 地球物理学报 , 2008, 51(5): 1520–1526. Xu Y, Zhang J F. An irregular-grid perfectly matched layer absorbing boundary for seismic wave modeling. Chinese J. Geophys. (in Chinese) , 2008, 51(5): 1520-1526. |
| [12] | Sun C W, Jing P G, Pu Y. Seismic imaging in the mountainous area of Southern China:Statics correction compared with migration from topography. Expanded Abstracts of 79th Annual Internat SEG Meeting , 2009: 2965-2969. |
| [13] | 方伍宝, 朱海波, 王汝珍, 等. 基于SEG起伏地表模型数据的偏移方法比较. 勘探地球物理进展 , 2007, 30(5): 382–387. Fang W B, Zhu H B, Wang R Z, et al. Comparison of different migration methods based on SEG irregular topography model data. PEG (in Chinese) , 2007, 30(5): 382-387. |