地球物理学报  2015, Vol. 58 Issue (10): 3758-3770   PDF    
引用本文
杨继东, 黄建平, 王欣等. 2015. 复杂地表条件下叠前菲涅尔束偏移方法. 地球物理学报, 58(10): 3758-3770, doi: 10.6038/cjg20151026   
Yang J D,Huang J P,Wang X,et al. 2015. Prestack Fresnel beam migration method under complex topographic conditions. Chinese J. Geophys. (in Chinese), 58(10): 3758-3770, doi: 10.6038/cjg20151026   
复杂地表条件下叠前菲涅尔束偏移方法
杨继东, 黄建平, 王欣, 李振春, 段心意    
中国石油大学(华东)地球物理系, 青岛 266580
投稿时间：2014-12-27    最后修改时间：2015-07-26 .
基金项目：国家973课题(2014CB239006,2011CB202402),国家自然科学基金(41104069,41274124),山东省自然科学基金(ZR2011DQ016),中央高校科研业务费专项基金(R1401005A)联合资助.
作者简介：杨继东,男,1990年生,在读硕士,主要从事地震波束正演与成像研究.E-mail:yangjidong_china@163.com
通讯作者：黄建平,男,1982年生,理学博士,副教授,博士生导师(院聘).长期从事地震波正演模拟、高斯束偏移及最小二乘偏移方法研究.E-mail:jphuang_china@gmail.com
摘要：高斯束偏移虽然克服了Kirchhoff偏移不能处理多波至和单程波动方程偏移不能对陡倾构造准确成像的问题,但在复杂地表条件下其偏移精度取决于所选择的初始束宽度,即当初始宽度较小时,近地表成像精度较高,但此时中深层成像质量较差;反之当初始宽度较大时,中深层成像质量提高,但近地表成像精度降低.针对高斯束偏移中深层和浅层成像精度的矛盾,本文发展了一种适用于陆地复杂地表条件的叠前菲涅尔束偏移方法.基于惠更斯-菲涅尔原理,本文首先给出了菲涅尔束的概念及其表征的格林函数,并采用有效邻域波场近似理论和反褶积成像条件,导出了复杂地表条件下叠前保幅深度偏移公式.最后,针对常规旁轴射线追踪中的数值噪音,给出了一种压制策略.同高斯束偏移相比,本文方法不仅解决了中深层和浅层成像精度的矛盾,而且提高了复杂地表条件下平面波的分解精度,使得偏移结果更加准确可靠.典型的模型算例验证了本文方法的有效性和稳健性.
关键词： 复杂地表     高斯束     菲涅尔束     保幅成像     叠前深度偏移    
Prestack Fresnel beam migration method under complex topographic conditions
YANG Ji-Dong, HUANG Jian-Ping, WANG Xin, LI Zhen-Chun, DUAN Xin-Yi    
Department of Geophysics, China University of Petroleum (East China), Qingdao 266580, China
Abstract: Gaussian-beam migration overcomes the limitation of Kirchhoff migration in imaging the multiple arrivals and of one-way wave equation migration in imaging the steep-dip structures, but its imaging accuracy depends on the initial beam width under complex topographic conditions, i.e., that with small initial beam width, the imaging accuracy is high at near-surface but low at depth, whereas with large initial beam width, the imaging accuracy is improved in the deep but declines in the near-surface. Aiming at the problem in Gaussian-beam migration between near-surface and deep areas, we have developed a method of prestack Fresnel beam migration method under complex topographic conditions in this paper.
Based on the principle of Huygens-Fresnel, we use the first Fresnel zone radius to limit the beams' effective width and derive of the formula of Fresnel beam. And then, we adopt the steepest descent method to obtain the approximate Green function constructed with Fresnel beams. Further, with the theory of wave field approximation in effective vicinity and deconvolution imaging condition, we obtain the prestack amplitude-preserved migration formula under complex topography conditions. Finally, aiming at the numerical noise of conventional paraxial ray tracing algorithm, we proposed an optimized strategy.
Compared with Gaussian beam migration, the proposed method not only solves the contradiction of imaging accuracy in near-surface and at depth, but also improves the accuracy of plane wave decomposition, making the migration results more accurate and reliable. Typical numerical examples proved the validity and robust of our method.
Key words: Complex topography     Gaussian beam     Fresnel beam     Amplitude-preserved imaging     Prestack depth migration    
1 引言

随着油气勘探的逐渐发展，陆上地震勘探的工作重心逐渐向山地、戈壁滩、黄土塬、沙漠等复杂地表地区转移(贾承造，2004；赵文智等，2007).在这些探区，近地表高程和速度的急剧变化不再满足地表一致性假设，使用基于常规静校正的数据处理方法难以获得高质量的地震剖面(李振春，2014).然而，针对这些复杂地表探区的特殊偏移方法可以消除地表高程和近地表速度变化对地震波走时和振幅的影响，对地下地质构造准确成像，为后续的地震资料解释、属性提取以及储层预测提供可靠的地震资料，因此也一直是国内外地球物理学者研究的热点.

基于零阶射线近似的Kirchhoff偏移可以直接从起伏地表进行波场的积分延拓，且无需规则观测的数据，是解决复杂地表条件下偏移成像问题有效方法之一.Wiggins(1984)首先提出起伏地表条件下基于kirchhoff积分的波场外延和偏移公式.Gray和Marfurt(1995)进一步证明了直接从起伏地表进行Kirchhoff偏移要优于静校正+Kirchhoff偏移的处理结果.Jäger(2003)通过考虑实际观测面的倾角和局部观测排列的关系，提出了一种直接从起伏地表进行的真振幅Kirchhoff偏移方法.Kirchhoff偏移方法虽然对复杂地表条件具有天然的适应性，但是由于零阶射线理论本身的不足，不能对焦散区、阴影区和多波至区域准确成像，因此在复杂地表和复杂构造探区该方法往往不能满足精细勘探的需求.

基于波动方程的“基准面校正”、“零速层”法、“逐步累加”法和“波场上延”法也是解决复杂地表条件偏移成像问题的有效方法.Berryhill(1979)首先提出了“波动方程基准面校正”的思想，之后Yilmaz和Lucas(1986)，Schneider等(1995)先后将其应用于起伏海底和起伏陆表的“层替换”中，在一定程度上改善了起伏海底探区和陆地逆冲带的成像质量.Beasley和Lynn(1992)把“零速层”的思想引入到起伏地表成像中，使得常规偏移算法在起伏地表探区的应用成为现实，之后Lynn等(1993)、MacKay(1994)和Gray(1997)对其应用中存在的问题做了深入剖析.Reshef(1991)提出了直接在深度偏移时进行“逐步累加”的波场外推策略，在一定程度上消除了起伏地表对地震波场的影响.何英等(2002)结合“零速度层”和“逐步累加”延拓的思想提出了“波场上延”法，也能较好地消除起伏地表对地下构造成像的影响.以上四种方法均基于地震波场的整体外推算法，一般需要规则的地震数据，然而在复杂地表探区往往很难满足这一条件，因此上述方法的实用性有所限制.

基于旁轴射线理论的高斯束偏移方法不但克服了Kirchhoff偏移不能处理多波至和单程波动方程偏移不能对陡倾构造准确成像的缺点，而且保留了射类偏移方法对复杂地表条件的天然适应性.Hill(1990，2001)首先提出了高斯束偏移的基本框架，Gray(2005)、Gray和Bleistein(2009)进一步将其拓展到了起伏地表条件和真振幅偏移中.岳玉波等(2010)通过考虑起伏地表对地震波场走时的线性影响，提出了一种适用于复杂地表条件的“保幅延拓”高斯束偏移方法，Yang等(2014)将起伏地表引起的二次走时校正项和振幅校正项引入到高斯束偏移中，改善了起伏地表探区的近地表成像质量.

根据Zhu(2013)的文献资料可知，在地表高程和近地表速度场急速变化的探区，高斯束偏移的成像精度取决于所选择的初始束宽度，即当初始宽度较小时，近地表成像精度较高，但中深层成像质量较差；反之，当初始宽度较大时，中深层成像质量有所提高，但近地表的成像质量明显下降.

针对于此问题，Liu等(2014)利用第一菲涅尔带约束各向异性介质走时表的成像范围，在一定程度上改善了偏移质量.本文从初始束参数动态选择的角度，对中心射线附近波场走时和振幅的延拓范围加以控制，提出了一种新型地震波束——菲涅尔束，并将其应用于叠前保幅深度成像中，实现了一种直接从起伏地表进行偏移处理的新型地震波束成像方法.此外，针对传统旁轴射线追踪中的数值噪音，本文基于有效半宽度的概念给出了一种较好的压制策略.在实现上述偏移方法的基础上，本文通过两个典型的模型算例，验证本文方法的有效性和稳健性.

2 理论方法 2.1 菲涅尔束

首先，本小节简单介绍高斯束有效半宽度和第一菲涅耳带半径的求取方法.然后，根据惠更斯-菲涅耳原理，利用第一菲涅耳带半径限制高斯束的有效半宽度，导出菲涅耳束的表达式.最后，简单讨论菲涅耳束在传播过程中的基本性质.

根据Červený等(1982)和Müler(1984)可知，高斯束有效半宽度的表达式为

其中s为中心射线的弧长，ωref为参考频率，ε(s)为复值初始束参数，[p1(s)，q1(s)]和[p2(s)，q2(s)]为动力学射线追踪方程的两组基础解，“Im”为取虚部运算符.

由菲涅尔体追踪理论(Červený，1992)可知，第一菲涅尔带半径可近似表示为

由于这种方法在震源处失效，本文结合均匀介质中菲涅尔带半径的解析解，在此给出一种修正表达式

其中λave=2πvave/ωref，vave为模型的平均速度.

根据惠更斯-菲涅耳原理可知，中心射线附近的波场能量主要分布在第一菲涅尔带内，为使高斯束 的能量分布符合这一原则，我们使用如下的约束条件：

进一步将(2)式和(3)式带入(4)式，并考虑到 ε(s)=εr(s)－iεi(s)、εi(s)=0和q1p2－q2p1=1，可得到

注意到与高斯束方法不同，上式中ε(s)不再是一常数，而是随射线弧长动态变化的函数，这也是本文方法与高斯束方法的本质区别.

为了方便讨论，在此将(5)式中初始参数决定的地震波束命名为菲涅尔束，进一步参考高斯束的表达式，可知其在频率域具有如下形式：

其中ω为圆频率，s和n分别为二维射线中心坐标系的弧长分量和横向分量，τ(s)为中心射线的走时.

为了了解菲涅尔束传播的基本性质，在此考虑速度为2000 m·s-1均匀介质.根据高斯束和菲涅尔束的定义，可计算出如图 1所示的地震波束传播形态和有效半宽度.可以看出，菲涅尔束的有效半宽度在震源处约为λ/4，并且在后续的传播路径上以平方根规律增大(见图 1a和1c)；相应的，高斯束的有效半宽度在震源处约为 2 λ，并且在后续的射线路径上以双曲规律增大(见图 1b和1c).因此，相对于高斯束来说，菲涅尔束在整条射线路径上都具有较小的有效束宽度，这使得中心射线附近波场的走时和振幅计算更加准确，从而提高了波动方程解的精度.

图 1 均匀介质中高斯束(a)和菲涅尔束(b)的传播形态及它们的有效半宽度(c)(介质速度为2000 m·s-1， 波束参考频率为10 Hz) Fig. 1 Propagation (a, b) and the effective half-width (c) of Gaussian beam and Fresnel beam in homogenous medium where velocity is 2000 m·s-1, referenced frequency is 10 Hz
2.2 由菲涅尔束表征的格林函数

根据Červený等(1982)和Popov(1982)的相关文献可知，地下某点M的格林函数可近似表示为

式中UΦ(s，n)为以角度Φ出射的菲涅尔束，Φ(Φ，s)为其积分权系数.

如果上式中的积分权系数已知，便可利用(7)式求取地下任意一点的地震波场，为此我们将通过对比均匀介质中格林函数的解析解和由方程(7)导出的近似解求得Φ(Φ，s).

考虑速度为v0的均匀介质，并令

则(7)式可化简为

根据如图 2所示的几何关系，可知

图 2 由菲涅尔束表征的格林函数 Fig. 2 Schematic of representation of Green′s function with Fresnel beam summation
将它们带入(8)式可得

其中

对(11)式中f(Φ)两侧微分可得

其中

观察(13)式可知，当Φ=Φ0时，df(Φ)/dΦ=0，所以Φ0为函数f(Φ)的驻点，又f(Φ)的二阶导数为

因此根据积分的稳相近似原理，(9)式可近似表示为

对比上式和零阶近似射线理论(ART)中格林函数的解析表达式

可得到

因此由菲涅尔束表征的格林函数的最终表达式为

注意到上式与高斯束方法不同的是：格林函数的积分权系数不仅与射线出射角有关，也与中心射线的弧长有关，是根据每一射线步菲涅耳带的大小自适应计算得到的.

为了讨论角度间隔dΦ、数据频率ω和初始束宽度l(s0)对高斯束和菲涅耳束表征格林函数的影响，在此仍考虑速度为2000 m·s-1的常速介质.如图 3 所示，红色和蓝色线条分别为根据(19)式和高斯束 方法计算的格林函数，黑色线条为零阶渐进射线理论(ART)的解析解.可以看出：(1)在近源处，由菲涅尔束计算的格林函数比使用高斯束方法计算的格林函数精度高，这是因为高斯束较大的初始束宽使其在近源处主要表现为平面波，而菲涅尔束合理的能量分布范围使其在近源处主要表现为球面波，因此由菲涅尔束构造的点源波场更加符合理论值；(2)角度间隔对菲涅尔束方法和高斯束方法计算的格林函数影响不大(见图 3a—3c)，这是因为在一定射线密度下二者均能使计算区域得到充分照明；(3)数据频率越大，菲涅尔束方法和高斯束方法计算的格林函数精度越高，这是由于它们都是波动方程的高频渐进解，频率越高近似解的精度越高(见图 3a、3d、3e和3f)；(4)当初始宽度较大时，由高斯束表征的格林函数的精度在近源和远源处均较低，当 初始宽度较小时，近源处格林函数的精度有所提高，但远源处格林函数的精度明显下降(见图 3a、3g和3h).

图 3 零阶射线近似方法(ART)、高斯束方法和菲涅尔束方法表征的格林函数
介质速度为2000 m·s-1，参考频率为10 Hz. (a) dΦ=1°, ω=20 Hz, l(s0)=λ; (b) dΦ=5°, ω=20 Hz, l(s0)=λ; (c) dΦ=10°, ω=20 Hz, l(s0)=λ; (d) dΦ=1°, ω=5 Hz, l(s0)=λ; (e) dΦ=1°, ω=10 Hz, l(s0)=λ; (f) dΦ=1°, ω=50 Hz, l(s0)=λ; (g) dΦ=1°, ω=20 Hz, l(s0)=4λ； (h) dΦ=1°, ω=20 Hz, l(s0)=0.25λ. Fig. 3 Green′s function constructed by ART, Gaussian beam summation and Fresnel beam summation, wherein velocity is 2000 m·s-1, referenced frequency is 10 Hz
2.3 基于有效邻域波场近似的复杂地表菲涅尔束偏移公式

根据Červený和Peník(1984)提出的有效邻域波场近似理论，从起伏地表接收点R到地下成像点P的格林函数可近似表达为(具体推导见Yang等(2014))：

式中 x P=(xP，zP)、 x R=(xR，zR)和x L=(xL，zL)分别为P、R和L点的笛卡尔坐标，i= －1，pLx和pLz分别为束中心L的水平和垂直慢度，εpxL和vL分 别为束中心L的菲涅尔束初始参数和速度，A(x P，L)和T(x P，L)分别为从L传播至P的菲涅尔束振幅和走时，x =(xR-xL，zR-zL)T，上标‘T’代表矩阵转置，，αL为束中心处中心射线的出射角，

根据Gray和Bleistein(2009)所提出的波场保幅延拓公式，可以得到复杂地表条件下上、下行波场的表达式：

其中θR=αR－βR，αR和βR分别为检波点处射线的出射角和起伏地表的局部倾角，由于在初始点菲涅尔束为平面波处，所以有αR=αL，θR=αL－βR(见图 4)，θS具有类似的意义；dr为起伏地表上检波点距 离，PU(x R，x S，ω)为 x S点激发 x R接收的地震记录，G*(x P，x R，ω)为格林函数G(x P，x R，ω)的复值共轭.

图 4 起伏地表成像原理图 Fig. 4 Schematic of imaging principle

在波场反向延拓中，为了使检波点接收的地震记录与菲涅尔束表征的波场相匹配，需要对地震记录进行加窗处理，在此采用Hill(1990)给出的高斯窗函数

式中ΔL为束间隔.

将(20)式和(22)式带入(21)式可得到上、下行波场的最终表达式：

其中

进一步利用反褶积成像条件

可以得到最终的菲涅耳束成像公式：

为了提高计算效率，在此采用Gray(2005)提出的稳相近似法，可得到如下的快速菲涅耳束偏移公式：

其中

(27)式中，ΔL，l0和dpx是菲涅尔束偏移的关键参数，它们决定了计算效率和成像精度，在此给出它们的经验表达式:

其中ωh、Zm和Φmax分别为地震资料的最高频率、速度模型的最大深度和射线追踪的最大初始角度.

2.4 旁轴射线追踪中数值噪音压制方法

在地震波束偏移的射线追踪中，首先利用初值射线追踪方法计算中心射线路径的坐标、实值走时和振幅，然后在中心射线附近利用递归算法在粗网格上计算旁轴射线的复值走时和振幅，并把振幅的大小作为是否进行下一个递归运算的判定条件.但是此算法在中心射线最后步长附近往往出现如图 5a箭头所示的偏移噪音和图 5b箭头所示的旁轴射线追踪误区，其产生的本质原因在于递归算法终止条件的不完备.

图 5 常规旁轴射线追踪算法
(a) Marmousi数据集单炮偏移结果； (b)单条波束传播形态. Fig. 5 Conventional paraxial ray tracing algorithm
(a) One-shot migration result of Marmousi dataset; (b) The propagation of single beam.

为了消除这种由于计算方法导致的数值噪音，在此我们在计算最后射线步振幅时增添另一个判定条件，即如果当前粗网格位于最后射线步的有效邻域内，则计算其振幅和走时，并返回逻辑值“真”，进行下一个递归运算，否则返回逻辑值“假”，并终止当前递归运算.有效邻域的半径定义为

其中send为最后射线步的弧长.

将上述优化算法应用到与图 5相同的数据中，可得到如图 6所示的偏移结果和波束传播状态，可以看出，单炮偏移结果中的数值噪音和旁轴射线追踪误区均已消失，验证了上述优化算法的有效性.

图 6 优化后的旁轴射线追踪算法
(a) Marmousi数据集单炮偏移结果； (b)单条波束传播形态. Fig. 6 Optimized paraxial ray tracing algorithm
(a) One-shot migration result of Marmousi dataset; (b) Propagation of single beam.
3 数值算例

在实现上述偏移方法的基础上，本小节首先通过常速起伏地表模型对方法的正确性和保幅性进行了测试.

如图 7a所示，此常速起伏地表模型网格大小为501×501，纵横向采样间隔分别为10 m和15 m，速度为2000 m·s-1，该模型中的反射是由介质密度差引起的，并且三个水平界面的垂直反射系数相同.有限差分正演的单炮记录如图 7b所示，该炮集共501道，道间距为15 m，采样间隔为4 ms，记录长度为5 s.图 8为采用高斯束方法和菲涅尔束方法得到的偏移结果，图 9为相应的单道归一化振幅提取值.可以看出，虽然高斯束方法能够消除起伏地表和偏移距对同相轴的走时影响，得到基本正确的成像结果，但较大的束宽度会使其在反射界面附近产生如图 8a箭头和图 9a—9c红色箭头所示的相干噪音，这些偏移噪音又进一步影响了该方法的保幅性(见图 9a—9c绿色椭圆)，使其不能获得准确的振幅信息.相对来说，本文方法使用了窄而稳定的菲涅尔束，使得偏移结果更加清晰、干净，并且可以使有效孔径内(图 8b箭头中间)中深层的反射能量得到一定恢复(见图 9d—9f，三个反射界面振幅基本一致)，验证了本文方法具有更好的振幅保持性.

图 7 常速起伏地表模型
(a)密度模型；(b)单炮记录. Fig. 7 Constant-velocity irregular topographic model
(a) Density model; (b) Single shot record.

图 8 常速起伏地表模型偏移结果
(a) 高斯束方法; (b)菲涅尔束方法. Fig. 8 Migration results for constant-velocity irregular topography model
(a) Gaussian beam method; (b) Fresnel beam method.

图 9 常速起伏地表模型偏移结果单道提取值
(a) 高斯束方法，CDP=150；(b) 高斯束方法，CDP=250； (c) 高斯束方法，CDP=350； (d) 菲涅尔束，CDP=150； (e) 菲涅尔束，CDP=250； (d) 菲涅尔束，CDP=35. Fig. 9 Traces extracted from migration results of constant-velocity irregular topographic model
(a) Gaussian beam method, CDP=150; (b) Gaussian beam method, CDP=250; (c) Gaussian beam method, CDP=350; (d) Fresnel beam method, CDP=150; (e) Fresnel beam method, CDP=250; (f) Fresnel beam method, CDP=350.

在计算效率方面，由于本文方法的初始束宽度选为λave/4，相同炮集内分解的高斯窗约为高斯束方法的4倍.因此理论上来说菲涅尔束偏移的处理时间也应是高斯束方法的4倍，但是由于其在每个高斯窗内地震道数减少，单次计算的效率有所提高，所以实际计算时间约为高斯束方法的2~3倍(在Dell-Optiplex390微机上，图 8a计算时间为86 s，图 8b计算时间为203 s).

为了进一步验证本文方法的适用性以及对偏移速度场的敏感性，在此对加拿大逆掩断层数据集进行测试.该数据集共278炮，单炮最大道数为480道，道间距15 m，中间激发，两侧接收，时间采样间 隔4 ms，记录长度8 s，最大高程差达1537 m.图 10a为其速度模型，网格大小为1668×1000，纵横向间隔分别为10 m和15 m.该模型不但地表高程和近地表速度变化剧烈，还具有大曲率褶皱、逆冲断层等复杂地质构造，是测试深度域起伏地表偏移方法的经典模型.

图 10 加拿大逆冲断层模型及偏移结果
(a)速度误差为5%的偏移速度场； (b)高斯束偏移(旁轴射线追踪优化前)，l0=λave； (c)高斯束偏移(旁轴射线追踪优化后)， l0=λave； (c)高斯束偏移，l0=4λave； (d)高斯束偏移，l0=0.25λave； (d)菲涅尔束偏移结果. Fig. 10 Canada over-thrust faults model and corresponding migrated results
(a) Velocity model with 5% error; (b) Gaussian beam migration (before paraxial ray tracing optimization) with l0=λave; (c) Gaussian beam migration (after paraxial ray tracing optimization) with l0=λave; (d) Gaussian beam migration with l0=4λave; (e) Gaussian beam migration with l0=0.25λave; (f) Fresnel beam migration.

图 10b和10c分别为旁轴射线追踪优化前后的高斯束偏移结果.可以看出，除了图 10b深部偏移噪声稍弱外(箭头处)，两者差异较小.究其原因可知，该数据集的高覆盖次数使得单炮偏移中由于传统递归射线追踪造成的数值噪音(如图 5a所示)淹没在了多次叠加后的强反射同相轴中.因此本文所给出的旁轴射线追踪改进方法只是偏移算法的一种优化策略，对单炮成像结果的改进效果显著，对多次叠加后的成像结果影响不大.

为了解初始束宽度对高斯束偏移质量的影响，在此进一步选用初始半宽度为 l0=4λave和l0=0.25λave对该数据集进行试算，得到如图 10d和10e所示的偏移结果，图 10f为菲涅耳束偏移结果，图 11为它们在CDP=800~1201、Depth=0~3000 m范围内的局部放大结果.观察对比图 10和11可知：(1)当高斯束初始宽度为l0=4λave时，不仅浅层成像噪音较大(见图 10d矩形框和图 11b椭圆)，中深层偏移精度也较低(见图 10d椭圆)，这是由于较大的初始束宽度使整条射线路径附近延拓的走时和振幅精度较低造成的；(2)当l0=λave时(Hill建议的初始束参数)，虽然中深层成像质量有所改善，但浅部仍存在少量噪音(见图 10c矩形框和图 11a椭圆)，究其原因可知，虽然初始束宽度减小为λave，但是由其对应高斯束构造的地震波场在近源处仍表现为平面波性质，与实际的球面波场不一致，产生了少量偏移噪音；(3)当l0=0.25λave时，近地表的偏移噪音己基本消除，但中深层反射界面变得模糊不清(见图 10e椭圆和图 11c箭头)，这主要是因为较小的初始束宽使高斯束的有效半宽度在远源处变得很大，导致中深层中心射线附近计算的走时和振幅精度较低，从而降低了偏移精度；(4)观察图 10f可知，菲涅耳束偏移不仅可以消除浅层的偏移噪音(见图 10f矩形框和图 11d)，也可对中深层的反射界面清晰成像，这主要取决于其合理的能量分布范围，使其构造的波场不仅在近源处与实际的球面波场一致，而且在远源处也满足地震波的能量分布规律.

图 11 偏移结果局部放大对比图(CDP范围901~1200，深度范围0~3000 m)
(a) 高斯束偏移，l0=λave； (b) 高斯束偏移，l0=4λave； (c) 高斯束偏移，l0=0.25λave； (d) 菲涅尔束偏移结果. Fig. 11 Comparison of local migration results CDP= [901 1200] and z= [0 3000]
(a) Gaussian beam migration with l0=λave; (b) Gaussian beam migration with l0=4λave; (c) Gaussian beam migration with l0=0.25λave; (d) Fresnel beam migration.

为了验证本文方法对速度误差的敏感性，在此进一步使用误差为10%和20%的速度场(见图 12a和12c)进行偏移测试，得到了如图 12b和12d的偏移结果.对比图 10f、12b和12d可以发现，随着速度误差的增大，成像结果的偏移噪音有所增加(见图 12c，12d箭头)，但整体的地质结构面貌仍清晰可见，大曲率褶皱、陡倾断面和底部倾斜反射层仍能够准确成像.以上分析说明，本文方法对深度域速度场不慎敏感，对速度误差具有较好的容忍度.

图 12 速度误差分别为10%(a)和20%(c)的偏移速度场及相应的菲涅尔束偏移结果(b，d) Fig. 12 Migration velocity models with 5% error (a), 10% error (c), and corresponding Fresnel beam migration results (b, d)
4 结论与讨论

基于惠更斯-菲涅耳原理和有效邻域波场近似理论，本文提出了一种适用于陆地复杂地表条件的菲涅耳束偏移方法，并通过典型模型算例，验证了本文方法的有效性和稳健性.通过使用菲涅耳带半径限制地震波束的有效半宽度，导出了一种新型的地震波束，并将其应用于叠前深度偏移中，解决了高斯束偏移中深层和浅层成像精度矛盾的问题；通过使用有效邻域波场近似理论导出的二次走时校正项和振幅校正项，提高了记录波场在复杂地表的平面波分解精度，进一步改善了成像质量.模型测试结果表明，本文方法虽然牺牲一定的计算效率，但是却大大 改善了复杂地表探区的成像质量，尤其是对浅部反射层.此外，试算结果也表明本文方法对偏移速度误差也具有较好的容忍度.由此可见，菲涅耳束偏移将会成为地震叠前深度成像的有效工具之一.

致谢 感谢两位外审专家对本文提出的宝贵意见和建议.

参考文献
[1] Beasley C, Lynn W. 1992. The zero-velocity layer: Migration from irregular surfaces. Geophysics, 57(11): 1435-1443.
[2] Berryhill J R. 1979. Wave-equation datuming. Geophysics, 44(8): 1329-1344.
[3] Červený V, Popov M M, Pšeník L. 1982. Computation of wave fields in inhomogeneous media-Gaussian beam approach. Geophysical Journal International, 70(1): 109-128.
[4] Červený V, Soares J E P. 1992. Fresnel volume ray tracing. Geophysics, 57(7): 902-915.
[5] Červený V, Pšeník L. Gaussian beams in elastic 2-D laterally varying layered structures. Geophysical Journal International, 1984, 78(1):65-91.
[6] Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration. Geophysics, 74(2): S11-S23.
[7] Gray S H, Marfurt K J. 1995. Migration from topography: Improving the near-surface image. Canadian Journal of Exploration Geophysics, 31(1-2): 18-24.
[8] Gray S H. 1997. Where is the zero-velocity layer? Geophysics, 62(1): 266-269.
[9] Gray S H. 2005. Gaussian beam migration of common-shot records. Geophysics, 70(4): S71-S77.
[10] He Y, Wang H Z, Ma Z T, et al. 2002. Pre-stack wave equation depth migration for irregular topography. J. Progress in Exploration Geophysics (in Chinese), 25(3): 13-19.
[11] Hill N R. 1990. Gaussian beam migration. Geophysics, 55(11): 1416-1428.
[12] Hill N R. 2001. Prestack Gaussian-beam depth migration. Geophysics, 66(4): 1240-1250.
[13] Jäger C. 2003. True-amplitude Kirchhoff migration from topography.//SEG Technical Program Expanded Abstracts.
[14] Jia C Z. 2004. The great discovery and future strategy of Chinese petroleum exploration in recent years.//Chinese Petroleum Geological Conference Proceedings (in Chinese).
[15] Li Z C. 2014. Research status and development trends for seismic imaging technology. Oil Geophysical Prospecting (in Chinese), 49(1): 1-21.
[16] Liu S Y, Liu T C, Wang H Z. 2014. Fresnel Beam Migration in Anisotropic Media.//76th EAGE Conference and Exhibition.
[17] Lynn W, MacKay S, Beasley C J. 1993. Efficient migration through complex water-bottom topography. Geophysics, 58(3): 393-398.
[18] MacKay S. 1994. Efficient wavefield extrapolation to irregular surfaces using finite differences: Zero-velocity datuming.//Technical Program Expanded Abstracts, 1564-1567.
[19] Müller G. 1984. Efficient calculation of Gaussian-beam seismograms for two-dimensional inhomogeneous media. Geophysical Journal International, 79(1): 153-166.
[20] Popov M M. 1982. A new method of computation of wave fields using Gaussian beams. Wave Motion, 4(1): 85-97.
[21] Reshef M. 1991. Depth migration from irregular surfaces with depth extrapolation methods. Geophysics, 56(1): 119-122.
[22] Schneider Jr W A, Phillip L D, Paal E F. 1995. Wave-equation velocity replacement of the low-velocity layer for overthrust-belt data. Geophysics, 60(2): 573-579.
[23] Wiggins J W. 1984. Kirchhoff integral extrapolation and migration of nonplanar data. Geophysics, 49(8): 1239-1248.
[24] Yang J D, Huang J P, Wang X, et al. 2014. Amplitude-Preserved Gaussian Beam Migration Based on Wave Field Approximation in Effective Vicinity Under Rugged Topography Condition.//2014 SEG Annual Meeting. Society of Exploration Geophysicists.
[25] Yilmaz O, Lucas D. 1986. Prestack layer replacement. Geophysics, 51(7): 1355-1369.
[26] Yue Y B, Li Z C, Zhang P, et al. 2010. Prestack Gaussian beam depth migration under complex surface conditions. Applied Geophysics (in Chinese), 7(2): 143-148.
[27] Zhao W Z, Hu S Y, Dong D Z, et al. 2007. Petroleum exploration progresses during the l0th Five-Year Plan and key exploration domains for the future in China. Chinese J. Petroleum Exploration and Development (in Chinese), 34(5): 513-520.
[28] Zhu T F. Complex-Beam Migration: Formulation and Comparisons.//2013 SEG Annual Meeting. Society of Exploration Geophysicists.
[29] 何英, 王华忠, 马在田等. 2002. 复杂地形条件下波动方程叠前深度成像. 勘探地球物理进展, 25(3): 13-19.
[30] 贾承造. 2004. 中国石油近年油气重大发现与未来勘探战略. 中国石油地质年会论文集.
[31] 李振春. 2014. 地震偏移成像技术研究现状与发展趋势. 石油地球物理勘探, 49(1): 1-21.
[32] 岳玉波, 李振春, 张平等. 2010. 复杂地表条件下高斯波束叠前深度偏移. 应用地球物理, 7(2): 143-148.
[33] 赵文智, 胡素云, 董大忠等. 2007. “十五” 期间中国油气勘探进展及未来重点勘探领域. 石油勘探与开发, 34(5): 513-520.