1 引言

随着我国油气勘探的不断发展，陆上地震勘探的工作重心逐渐向西部和南部的山地、戈壁滩、黄土塬和沙漠等复杂起伏地表地区转移(赵文智等，2007；贾承造，2004).在这些起伏地表探区，近地表高程和速度的急剧变化使得地表一致性假设不再成立，因此基于地表一致性假设的常规数据处理方法在这些探区难以满足勘探精度的要求.然而，直接从这些复杂地表进行偏移的成像方法可以较好地消除地表高程和近地表速度变化引起的走时、振幅误差，对地下地质结构准确成像，为后续的地震资料解释、属性提取以及储层预测提供可靠的地震资料，所以一直是国内外地球物理学者研究的热点和重点.

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

基于渐进射线理论的高斯束偏移是一种较为准确、灵活、高效的深度域地震成像方法，并可以直接从起伏地表进行波场的延拓成像，无需规则的观测数据，对复杂地表条件有天然的适应性.Hill(1990，2001)首先提出了共偏移距域的高斯束偏移方法，之后Gray(2005)将其推广到共激发点道集内，并探讨了几种实现方法的优劣.Gray和Bleistein(2009)基于单程波保幅延拓公式，提出了真振幅高斯束偏移方法，并建议用“局部静校正”处理陆地探区的起伏地表问题.郭朝斌等(2011)对高斯束偏移中的角度域成像、起伏地表成像及束参数的选取做了详细阐述.岳玉波等(2010)针对于复杂地表条件提出了一种“保幅延拓”方法，可较好地消除起伏地表对不同角度出射高斯束的走时影响.“局部静校正”法和“保幅延拓”法均只考虑了束中心附近高斯窗内地表高程和近地表速度场对走时的线性影响，忽略它们对二次时差和振幅的影响，而当地表高程和近地表速度变化剧烈时，这必然导致成像质量的降低.

本文基于Červený和Pšenčík(1984)所提出的高斯束中心射线附近有效邻域波场近似理论，推导了起伏地表高程和近地表速度场对地震波场的一次和二次时差校正项以及振幅校正项，在验证所导出公式正确性的基础上，将其应用于保幅高斯束深度偏移中，实现了一种精度更高的起伏地表保幅高斯束偏移方法.此外，考虑到高斯束偏移是在角度域实现，本文进一步给出了一种更加准确的基于有效邻域走时近似的旁轴射线角度计算方法.在实现上述成像方法的基础上，本文通过两个典型的模型算例，验证方法的有效性和适应性.

2 理论方法 2.1 起伏地表条件下有效邻域波场近似公式

假设高斯束中心射线Ω在起伏地表L点出射，接收点R在中心射线的投影为R₀(如图 1所示，蓝色为射线中心坐标系，绿色为笛卡儿坐标系)，则在射线中心坐标系下该条高斯束对接收点波场贡献为：

(1)

其中ω为圆频率，A(R₀)和τ(R₀)分别为R₀点的复值振幅和实值走时，M(R₀)=P(R₀)/Q(R₀)，P、Q和M为高斯束的复值动力学参数，n为接收点在射线中心坐标系中的横向坐标.根据Červený和Pšenčík(1984)所提出的有效邻域波场近似方法，可得到：

(2)

式中

(3)

其中α^L为中心射线在L点与笛卡儿坐标系z轴正方向的夹角，(x^R，z^R)和(x^L，z^L)分别为接收R和束中心L的笛卡儿坐标，v为中心射线的速度，s为射线中心坐标系中的弧长，“T”代表矩阵转置.将以上各参数代入(2)式，经化简可得到(具体推导见附录A)：

(4)

式中ω₀和l₀分别为动力学射线追踪初始参数中的参考频率和水平初始宽度，p_x^L和p_z^L分别为束中心的水平和垂直慢度.

考虑(4)式各项的物理含义可知：第一个指数项为接收点波场相对于束中心波场的振幅校正项，第二个指数项为相应的相位校正项.振幅校正项中的是由接收点和束中心水平坐标不一致引起的，在常规高斯束偏移中表现为局部倾斜叠加的指数衰减项；和是由于地表高程差异引起的振幅衰减项.目前在国内外起伏地表高斯束偏移方法中还未曾考虑这两项的作用.相位校正项中τ(L)为束中心的实值走时，在偏移处理中作为时间的初始值，一般取为零；p_x^L(x^R－x^L)是由接收点和束中心水平坐标不一致引起的走时差，在常规高斯束偏移中体现为局部倾斜叠加的相位校正项，最早由Hill(1990)提出；p_z^L(z^R－z^L)是由于地表高程差异引起的相位校正项，岳玉波等(2010)在“保幅延拓”法中考虑了此项的影响，并在一定程度上改善了复杂地表高斯束成像的质量；为由起伏地表和近地表速度变化引起的二次时差校正项，目前国内外学者未曾考虑该项对偏移结果的影响，但当地表高程和近地表速度剧烈变化时此项的影响不可忽略.

为了定量说明起伏地表引起的振幅校正和二次时差校正的作用，在此对束中心速度函数为v=(2000+5x+5z)m·s^-1，x^R－x^L和z^R－z^L均为-50 m的模型进行理论分析.根据局部静校正法、保幅延拓法和公式(4)可以得到如图 2所示的试算结果：由图 2a可知，当入射角较小时(±10°左右)，由于p_z≈1/v₀，p_x≈0，常规局部静校正方法就成了岳玉波保幅延拓法的特例，此时两者走时校正量基本一致；但是当入射角较大时，由于p_z≠1/v₀，p_x≠0，两者将存在明显差异；在此基础之上，公式(4)进一步考虑了束中心附近速度梯度的影响，当入射角较小时，由于p_x和p_z对时差校正量影响较小，所以速度梯度对时差的影响就凸显出来(如图 2a所示，绿色曲线在小角度处和另外两条相差较大)，而当入射角较大时，p_x和p_z对时差校正量的影响又起主要作用，所以此时本文方法和保幅延拓法一致(如图 2a所示，绿色与蓝色曲线在大角度处基本重合).观察图 2b可知，常规方法的振幅校正量在此模型中仅为本文方法在入射角为0°和66.5°时的特例，这是因为常规起伏地表高斯束偏移只考虑了x^R－x^L对振幅的影响，而忽略了地表高程和入射角的作用.

由以上分析可以看出，为了在地表高程和近地表速度剧烈变化的地区获得更准确的成像结果，应将本文推导的振幅校正项和二次时差校正项引入起伏地表高斯束偏移中.

2.2 保幅高斯束偏移公式

高斯束偏移的关键是使用不同出射方向的高斯束来表征格林函数，如果高斯束从L点出发传播至成像点P，那么基于2.1节推导的相位和振幅校正项，从束中心L有效邻域内接收点R传播至P点的格林函数可近似表示为：

(5)

其中A(x^P，x^L)和T(x^P，x^L)分别为从L传播至P点高斯束的复值振幅和复值走时，W和X的表达式见(3)式，上标代表所对应的位置.

式(5)中的格林函数既包含了Hill(1990)所提出的由检波点和束中心水平坐标差异所引起的相位校正项(iωp_x^L(x^R－x^L))，也包含了岳玉波等(2010)所提出的由检波点和束中心高程差异所引起的相位校正项(iωp_z^L(z^R－z^L))，以及2.1节推导的由近地表速度梯度所引起的二次相位校正项和高程差异引起的振幅校正项，其近似精度更高.在具体实现时，如果傅氏变换核的指数为正，则相位校正的符号也为正(本文中采用物理学家所选择的符号约定，即取正傅里叶核)(Claerbout，1991)，否则为负.

根据Gray和Bleistein(2009)真振幅高斯束偏移中记录波场的反向延拓公式，我们可得到起伏地表条件下的上行波场表达式:

(6)

式中θ^R=α^R－β^R，α^R、β^R分别为接收点处中心射线的出射角和界面的倾角，由于高斯束在初始点处波前为平面，故α^R=α^L，所以有θ^R=α^L－β^R(角度关系见图 3)，dr为检波点之间的实际距离，P_U(x^R，x^S，ω)为记录波场，G^*(x^P，x^R，ω)为格林函数的复值共轭.

考虑到高斯束振幅随横向距离呈指数衰减的性质，以及为了确保由高斯束所构造的波场在地表与记录波场相匹配，需要对记录波场沿水平方向做加窗处理，在此选择Hill(1990)所给出的窗函数：

(7)

其中ΔL为水平束间隔.把(5)和(7)式代入(6)式中可以得到起伏地表条件下由高斯束表征的上行波场：

(8)

其中

(9)

为基于有效邻域波场近似的局部倾斜叠加，此特殊的局部倾斜叠加不仅考虑到了倾角和高程对地震波走时的线性影响，还考虑到它们对振幅以及走时的二次影响.

同理可得到由高斯束表征的下行波场：

(10)

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

(11)

可得到起伏地表条件下保幅高斯束偏移公式：

(12)

为了提高计算效率，将式(12)从炮点和接收点的水平慢度域(p_x^S，p_x^L)变换到中心点和偏移距的水平慢度域(p_x^m，p_x^h)中，并利用最速下降近似求取关于p_x^h的内层积分和震源照明项P_D(x^P，x^S，ω)P_D^*(x^P，x^S，ω)(Gray和Bleistein，2009)，可以得到高效的偏移公式：

(13)

其中：T^*是复值走时T的复值共轭，

(14)

此偏移公式充分考虑了起伏地表的高程、倾角及近地表速度梯度对地震波走时和振幅的影响，并做了相应的校正，具有更高的成像精度.

2.3 旁轴射线传播角度计算方法

根据Červený和Pšenčík(1984)有效邻域走时近似理论，可知旁轴射线上一点P的实值走时可表示为：

(15)

P、Q的关系如图 4所示，p_x^Q和p_z^Q分别为Q点的水平和垂直慢度，Re(W)为2×2的复值矩阵W的实部，具体表达见(3)式.

对式(15)沿x方向求导可得：

(16)

其中W₁₁和W₁₂分别为矩阵W左上角和右上角的元素.考虑到，可得：

(17)

同理可得到：

(18)

则P点的传播角度为：

(19)

式中arctan2为值域为[-π，π]的反余弦函数，在C和Fortran编程语言中均有此库函数.

由于在具体实现时，经搜索求得中心射线上与P点相距最近的时间步往往为其有效邻域内任意的一点，如图 4中的Q点，而不一定是S₀点，所以就具体实现来说，用(19)式求得的旁轴射线的传播角度比郭朝斌等(2011)和岳玉波等(2010)给出的方法更加准确.由于Červený和Pšenčík(1984)有效邻域波场近似理论中只给出了复值走时和振幅的计算方法，而没有给出旁轴射线传播角度的计算方法，所以式(18)和式(19)实质上也是对这一理论的完善和补充.

3 数值算例

在实现本文方法的基础上，本小节首先通过简单的常速起伏地表层状模型对方法的保幅性和正确性进行了测试.如图 5a所示，此起伏地表模型的网格大小为801×1001，纵横向采样间隔分别为4 m和10 m，速度为2000 m·s^-1，在此模型中反射是由介质的密度差引起，并且两个水平界面的垂直反射系数相同，起伏地表高程最大相差800多米.图 5b为该模型的声波正演单炮记录，该炮集共有801道，道间距为10 m，采样间隔为1 ms，记录长度为4 s.图 6为使用不同偏移方法和成像条件得到的成像结果以及在CDP=400处提取的单道归一化振幅，图 7为利用局部校正法和本文方法得到的CDP=400处的角度域共成像点道集(ADCIG)，仔细对比图 6和图 7可以得出以下结论：(1)相对于图 6a和6d的互相关成像结果，本文方法可以较好地补偿由几何扩散所引起的振幅衰减，准确地恢复了地下界面的垂直反射系数(见图 6f)，证明了本文方法具有一定的振幅保持性；(2)如图 6b和6e所示，基于局部静校正和反褶积成像条件的高斯束偏移结果虽然在反射界面处正确恢复了垂直反射系数，但是仔细观察图 6e可知，在反射界面之外的深度处存在因时差校正不准确造成的偏移噪声(见图 6e红色矩形框)，而本文方法由于采用了较为准确的走时校正项，可消除这些成像噪声(见图 6f)；(3)从提取的ADCIG也可看出，相对局部静校正方法，本文方法精确的时差校正使得同相轴连续性更好(见图 7箭头)，准确的振幅校正使得大角度的偏移能量得到了较好的控制(见图 7椭圆)，减少了ADCIG叠加时的成像噪声.

此简单模型的试算结果表明，本文方法在成像精度和噪声压制方面较常规高斯束偏移方面有一定的提高，这主要是因为：基于有效领域波场近似的精确时差校正和振幅校正，不仅使得有效反射能量在相干叠加时相位更加一致，而且对远偏移距处的能量有较好的控制，因此从根本上压制了偏移噪声的出现；另外，本文方法采用了保幅的波场延拓公式和反褶积成像条件，可以对中深部的反射能量进行一定补偿，较为准确地恢复界面的垂直反射系数.

为了进一步验证本文方法的适用性，在此对国际标准SEG起伏地表模型(Amoc和BP公司设计的加拿大起伏地表逆掩断层模型)数据进行测试，该数据共278炮，单炮最大道数为480道，道间距15 m，中间激发，采样率为4 ms，记录长度8 s，最大高程差达1537 m.图 8a为其速度模型，网格大小为1668×1000，纵横向间隔分别为10 m和15 m，此模型不但地表高程、近地表速度变化剧烈，又具有地下复杂的褶皱和逆冲断层，是测试起伏地表偏移方法的经典模型.图 8b为该数据集的第250炮记录，可以看出该剖面的同相轴已不再呈双曲形态(见图 8b椭圆)，并且互相交错(见图 8b箭头)，从侧面反映了该逆冲断层模型的结构复杂性.

图 9为采用不同偏移方法和成像条件得到的偏移结果，图 10为CDP=800的归一化反射系数和偏移结果的单道提取值.观察图 9a和10b可知，基于有效邻域近似波场和互相关成像条件的高斯束偏移方法，虽然可以对地下复杂构造正确成像，但是深层能量较弱(见图 9a箭头和图 10b)；另外，近地表大量偏移噪声(见图 9a矩形框)，掩盖了近地表的有效反射层(见图 10b，红色箭头处振幅远大于相应深度处的反射系数)；观察图 9b和图 10c可知，利用局部静校正和反褶积成像条件的成像高斯束方法，虽然深层能量得到一定补偿(见图 10c，5~10 km的反射能量相对图 10b中相应深度处的反射能量较大)，但是浅层噪声仍然较大(见图 9b矩形框和图 10c箭头)，使得深层有效反射能量相对较弱(见图 9b箭头和图 10c中5~10 km区域).观察图 9c和图 10d可知，利用本文方法，浅层噪声明显得到压制(见图 9c矩形框和图 10d箭头)，深层能量得到明显补偿(见图 10d中5~10 km的反射能量相对图 10b和10c相应深度处明显增强)，并且整体能量分布较均衡(见图 10d和10a，偏移结果与反射系数更一致)，这主要取决于本文方法引入的振幅校正项和更准确的时差校正项.但是由于使用反褶积成像条件也相对地放大了炮记录中噪声的成像结果(见图 10d零轴附近的高频毛刺干扰)，关于这些高频噪声的压制还需要做进一步的研究.

为了验证本文方法在单一CDP上的成像精度，在此提取了如图 11所示的CDP=500、800和1200的ADCIG道集，从图上可以看出同一深度上的同相轴均已拉平，并且深浅层能量基本一致，从而证实了本文方法对复杂地表条件的适应性和有效性.

4 结论与讨论

本文在常规高斯束偏移和有效邻域波场近似理论的基础上，发展了一种适用于复杂起伏地表条件的保幅高斯束偏移方法，并通过典型模型算例，验证了本文方法的有效性和适应性，通过测试结果的对比分析，可以得到以下认识：

(1) 本文方法在起伏地表进行平面波分解时，不但考虑了起伏地表的高程、倾角信息，还注意到近地表速度梯度的影响，并相应地引入二次时差校正项，获得更高的平面波分解精度，从而使得成像结果中因相位校正不准产生的偏移噪声大大减少；

(2) 本文方法首次在起伏地表高斯束偏移中引入因地表高程差引起的振幅校正项，使得远离束中心道集的能量得到合理控制，从而进一步使近地表的偏移噪声得到相应地衰减；

(3) 由于采用了Gray所提出的保幅延拓公式和反褶积成像条件，本文方法可以较好地恢复地下反射层垂直反射系数，获得振幅相对保持的成像结果.

虽然本文实现了基于有效邻域波场近似的起伏地表高斯束偏移方法，但是下一步还需要在如下几个方面做进一步的深入研究：

(1) 由于在高斯束偏移中，要想提高深层复杂构造的成像精度必须引入深度聚焦的思想，但这使得局部倾斜叠加的邻道依赖性增大，而要想提高近地表的成像精度又必须得减小这种依赖性，因此为了提高深层和浅层总体的偏移精度还需要进一步研究适应不同工区地质条件的地震波束偏移方法；

(2) 由于反褶积成像条件放大了炮记录中噪声的成像效果，因此要进一步从成像条件和滤波两方面探索相应的噪声压制方法；

(3) 结合反演的思想，进一步研究基于高斯束理论的反问题，探索适应于复杂起伏地表条件的高斯束反演方法.

附录A

在地表L点处，以角度α^L出射高斯束的动力学参数P(s)和Q(s)分别为：

(A1)

其中ω₀为参考频率，l₀为水平初始宽度，v^L为L点的初始速度，i为虚数单位.

结合(3)式可求得矩阵A的实部与虚部：

(A2)

将(A2)式代入(2)式，可得到出射点L附近有效邻域内高斯束的近似表达式：

(A3)

进一步令，就可得到(4)式.