0 引言

随着勘探的不断深入，我国油气战略重点转向了具有复杂地表和地下条件的双复杂地区.像中国西部山区这样的起伏地表条件限制了采集孔径范围，造成了深层构造照明不足、成像振幅不均衡，从而影响了成像质量和后续的属性分析.因此，开展起伏地表条件下的照明补偿方法研究对于提高起伏地表地区的成像品质有着重要的现实意义.

传统地震资料处理方法基于水平地表假设，不能直接用于地表起伏条件下的地震成像.起伏地表面的影响和深层弱照明问题是地表起伏地区成像面临的两个主要问题.目前工业生产用的各种静校正方法(高程、折射波、层析等)都是基于地表一致性的假设条件，地震射线垂直入射和出射，但是实际受高程和近地表速度的变化射线并不是垂直入射和出射，从而引起静校不“静”，波场延拓的方法能够解决这种静校不“静”的问题.针对高程基准面校正所带来的问题，Berryhill(1979)提出了波动方程基准面校正的方法，将已知任意形状观测面上的波场延拓到一个特定的基准面上.零速层法(Beasley and Lynn, 1992)是在地表面与基准面之间插入一个速度值很小的虚拟层，利用高程基准面静校正将野外数据校正到水平基准面上，地震波在该层中横向传播可以忽略不计，即用波动方程的方式抵消了高程校正的时移量.直接下延法可以避免偏移前的高程静校正(Reshef，1991；Yang et al., 1999；叶月明等，2008).Al-Saleh等(2009)直接从起伏面开始构建横向可变深度步长的显式波场延拓算子，相对于零速层技术有明显的成像优势.Shragge和Sava(2005)提出了基于保角变换的方法建立黎曼坐标系统以适应地表起伏面，再按照常规偏移方法成像.张新彦等(2017)利用数学变换手段将笛卡尔坐标系的不规则模型映射到曲线坐标系的规则模型，提出了起伏地形下的高精度反射波走时层析成像方法.随着计算机技术的发展，相对耗时的起伏地表逆时偏移和最小二乘偏移方法也在被逐步应用(曲英铭等，2015).侯爵等(2018)引入一种精度无损的处理起伏边界的模型参数化方法，采用零延迟归一化互相关成像条件实现了起伏地表条件下的弹性波场逆时偏移成像.周华敏等(2014)以地震波场在地下的照明强度作为迭代反演的预条件算子，加快了最小二乘偏移的收敛速度.近年来，基于高斯束的复杂地表偏移方法以其较高的偏移效率而被应用，岳玉波等(2012)提出了复杂地表条件下保幅高斯束偏移，黄建平等(2016)研发了基于有效邻域波场近似的起伏地表保幅高斯束偏移.有限的采集孔径造成地下照明不足，为提高地下照明度，Wu等(2004)提出了水平地表情况下的局部角度域照明补偿方法.基于倾斜叠加的局部角度域道集求取非常耗时(Xie and Wu, 2002; Xie et al., 2004)，而基于小波束分解的方法可以较快地得到局部角度域成像和照明补偿因子(Cao and Wu, 2009;毛剑等，2010).陈生昌等(2007)、陈生昌和王汉闯(2010)提出基于平面波照明的补偿方法，相对于局部角度域的照明补偿方法具有较高的计算效率.

有限偏移孔径造成的照明不足问题在起伏地表条件下更为严重，为了能够改善起伏地表地区的成像品质，本文提出了一种基于小波束偏移算子的起伏地表区照明补偿方法.在基准面与起伏面之间插入虚拟层，通过波场逐步累加的延拓方法解决了起伏面对波场延拓的影响，并引入空间滤波函数压制了虚拟层内偏移噪音.应用了小波束波场延拓算子，其局部窗化的参考速度减小了背景速度扰动量，提高了波场延拓算子的精度.局部指数标架小波束波场分解求得了局部波数域成像矩阵和照明补偿因子，最后在局部倾角域实现起伏地表条件下的照明补偿.SEG起伏地表模型补偿测试和振幅分析验证了本方法的有效性，深层照明度得到加强的同时也提高了成像振幅均衡度.

1 方法原理 1.1 起伏地表条件下的波场延拓

常规偏移方法都是基于地表面水平的假设，无法在地表面起伏的条件下直接应用.本文采用虚拟层内波场零速传播的方法来解决起伏面对波场延拓的影响.在起伏面与基准面间插入虚拟层(图 1)，通过对虚拟层内波场延拓的控制，使得虚拟层内无波场传播.为避免波场在起伏面处发生透射或反射，在虚拟层内填充了近地表速度.以下行波场u_d(x, z, w)为例，从基准面开始波场沿深度方向正向延拓，而在整个虚拟层内是以零速传播，当遇到起伏地表面时，将对应的震源波场累加到延拓波场中并继续正向延拓.对于上行波场u_u(x, z, w)，从基准面开始波场沿深度方向逆向延拓，在整个虚拟层内同样也是以零速传播，当遇到起伏面时，将对应的地震记录累加到延拓波场中，并继续沿深度方向逆向延拓，如式(1)和式(2)所示：

(1)

(2)

其中，ω是频率，(x, z)是空间位置，u_dc(x, z, ω)和u_uc(x, z, ω)分别上一层延拓至(x, z)处的下行波场和上行波场，u_s(x_s, z_s, ω)和u_r(x_r, z_r, ω)别是起伏地表面上对应震源位置(x_s, z_s)与检波点位置(x_r, z_r)的震源波场和地震记录波场，在每个延拓深度层都会加权空间滤波因子f(x, z), 以压制虚拟层内的噪音，如(3)式所示：

(3)

其中，u_df(x, z, ω)和u_uf(x, z, ω)分别是加权滤波因子后的下行和上行波场，空间滤波因子f(x, z)如式(4)所示：

(4)

1.2 基于局部余弦基的小波束波场延拓

传统单程波偏移算子中的背景速度一般为偏移孔径内的平均速度或最小速度，在地表高程起伏较大的地区，虚拟层内的速度与实际地层间的速度可能会存在较大差异，导致真实速度与背景速度间扰动量偏大，影响算子精度.小波束波场延拓算子在偏移孔径范围内划分多个窗，背景速度为每个窗内的平均速度，局部窗化的参考速度降低了速度扰动量，也提高了波场延拓算子精度.本文应用了基于局部余弦基(Local Cosine Basis，简称LCB)的小波束波场延拓算子(Wu et al., 2008)，该算子得益于快速LCB变换而具有较高的计算效率.

在频率-空间域，标量声波方程可以表示为

(5)

其中，v(x, z)是速度，u(x, z, ω)代表频率-空间域的波场，为了方便起见，后文用u(x, z)来代替频率-空间域的波场.

将深度z上的波场沿空间x方向窗化后，按照式(6)分解波场如下：

(6)

其中，是在空间窗x_n和波数下的小波束分解系数.〈, 〉表示内积，其中，〈u(x, z), b_mn(x)〉=∫u(x, z)b_mn^*(x)dx，*代表复共轭.b_mn(x)是分解的LCB小波束，表示为

(7)

其中，基原子由空间位置x_n和间距L_n=x_n+1-x_n决定，B_n(x)是相互重叠的平滑钟函数，将式(6)代入式(5)得到标量声波方程式(8)：

(8)

在上式中，是一系列关于深度z的常系数，传播的影响包含在小波束的演化过程中.对于单程波近似，可以将小波束演化表示为式(9)：

(9)

其中，a_mn(x)是小波束b_mn(x)经过介质传播后的波束，A_n是平方根算子.经过介质的传播后a_mn(x)不再是一个小波束，应用同样的基函数对a_mn(x)再进行小波束分解表示为

(10)

其中，传播算子矩阵在小波束域可以表示为

(11)

那么，延拓至深度z+Δz的波场在小波束域可以表示为

(12)

P_{jl, mn}是小波束传播算子矩阵，决定着小波束的传播和交互偶合，通过式(12)可以使波场在小波域从深度z延拓到z+Δz层.波场重建是波场分解的反过程，小波域的波场在频率空间域重建后可以表示为式(13)：

(13)

小波束波场延拓算子的特点在于局部窗化的参考速度，将每一炮的偏移孔径范围分成若干个窗，在每个窗内计算出该窗内的背景速度和每个网格点上的扰动速度，能够确定波场在空间局部位置处的波数，每个小窗内的波场就是一个小波束.实际资料处理中的每个小波束宽度取值在120 m左右，就能够满足小波束波场延拓算子的精度需要.

1.3 局部倾角域成像与照明补偿

本文的照明补偿是在角度域完成，通过波场分解将频率空间域波场变换到局部倾角域.在波场延拓过程中采用的是基于局部余弦基的小波束波场延拓算子，其中局部余弦基函数(式(7))还可以表示为式(14)，两个指数项表明该类小波束有两个对称的波瓣，即沿深度轴对称的两个方向传播，缺乏单一定义的方向性，不利于照明补偿.

(14)

事实上，如式(14)所示的基函数都是由局部指数标架线性组合而来的，局部指数标架定义为

(15)

基于式(15)局部指数标架分解的波场具有单一方向的角度，能够用于照明补偿.利用式(15)进行频率波数域波场重构(毛剑等，2010)，得到局部波数域上行波场和下行波场ω)，局部波数域成像条件表示为(Cao and Wu, 2009)

(16)

对于某一波数对应平面波传播角度如式(17)所示，根据该式可以将局部波数域成像转换到局部角度域成像L(x, z, θ_s, θ_g)：

(17)

其中，θ_s和θ_g分别表示地下局部位置入射角和散射角，地层的局部倾角θ_n更能反映出地下的局部真实情况，各角度之间的关系如图 2所示，倾角与地层反射角之间的关系如式(18)所示：

(18)

通过式(18)将局部角度域成像L(x, z, θ_s, θ_g)转换到局部倾角域L(x, z, θ_n)，局部角度域成像能够反映出地下局部反射点在不同方向的入射波场和反射波场对成像的贡献，称之为局部成像矩阵(Local Image Matrix, 简称LIM).局部散射矩阵(Local Scattering Matrix，简称LSM)则是反映了地下固有特征，描述振幅与散射能量间的一定比例关系.由于采集孔径限制以及传播路径等影响，LIM是被扭曲了的LSM，照明补偿的目的就是为了让LIM更加逼近LSM，从而最终反映出地下的真实情况(Wu et al., 2004).与水平地表情况类似，局部波数域照明补偿因子可以表示为(Cao and Wu, 2009)

(19)

其中，G_F是格林函数正演过程，A(x_g, x_s)是空间采集孔径.根据式(17)和式(18)，将局部波数域照明补偿因子变换到局部倾角域F(x, z, θ_n)，于是，局部倾角域照明补偿计算如式(20)：

(20)

其中，θ_min和θ_max分别是最小和最大地层倾角，ε是阻尼因子，以保证计算的稳定性，目前本方法的阻尼因子是基于经验值和测试后选定，选取的基本原则是比倾角域补偿因子小两个数量级，这样既保证了计算的稳定性，也不改变相对关系.

2 模型测试

为了验证本文方法的有效性，对加拿大逆掩断层模型(Gray and Marfurt, 1995)进行试算分析.该模型地表高程较大(图 3a)，最大高程近2 km，而且地表起伏剧烈，局部范围2.5 km内高程差能够达到1.5 km，并且发育多套高陡断层和逆掩构造(图 3b)，起伏地表面和上覆复杂构造影响了深层两套倾斜地层的照明度，图 3c是原始单炮记录，地表起伏造成的同相轴扭曲现象严重.模拟的激发方式是中间放炮，共277炮、480道接收、采样率2000、4 ms采样、道间距15 m、最大炮检距3600 m、最小炮检局为15 m、炮间距为90 m.

波场在介质中是以小波束的形式传播，图 4a和4b分别是15 Hz和40 Hz情况下的小波束传播，二者都是下行波场，以垂直地表的方向传播.对比可以看到，较高频率的小波束在地下传播的角度扩散相对更大一些.图 5a、图 5b和图 5c分别是横向位置相同、不同深度三个位置的局部角度域成像矩阵，横、纵轴分别代表了入射角和散射角.从图中可以看出，深度8 km位置处的局部成像角度贡献更为均匀(图 5a)，局部成像矩阵中心在入射角和散射角都是0°附近.深度6 km位置处，局部成像矩阵中心在入射角-30°，散射角0°附近，而在深度4 km处，局部成像矩阵成椭圆形分布.局部倾角域能够更加客观地反演局部地层倾向，因此在局部倾角域进行照明补偿能够根据不同角度的照明情况来针对性地补偿.图 6a、6b、6c和图 6d分别是局部倾角为-30°、0°、30°和传统自动增益控制(Auto Gain Control，简称AGC)的补偿因子，不同局部倾角域补偿因子也表现出不同补偿的方向性.图 7a、7b和7c分别是局部倾角为-30°、0°、30°情况下的局部倾角域成像，图 7d、7e和7f分别是局部倾角为-30°、0°、30°情况下的局部倾角域成像照明补偿后结果，从图中可以明显看出，照明补偿后，不同角度的成像都有改善，尤其是解决了由于深层照明不足对构造成像带来的影响.

相对于传统的分步傅里叶叠前深度偏移(图 8a)，LCB小波束叠前深度偏移算子局部窗化的参考速度降低了背景速度扰动量，尤其是在起伏地表面附近的浅层成像效果有了显著的改善(图 8b).将所有倾角域补偿后的像叠加得到最终局部倾角域照明补偿结果，如图 8d所示，与补偿前(图 8b)相比，深层构造成像振幅的能量有了显著的提升，而且补偿效果上更优于传统AGC增益处理(图 8c).为了进一步分析照明补偿后的振幅，沿图 9a箭头所指层位，拾取了照明补偿前(图 8b)、AGC补偿后(8c)和照明补偿后(图 8d)的成像振幅如图 9b所示，AGC增益处理只是增强了深层振幅值，在层位上拾取的振幅值和照明补偿前的LCB偏移结果相近，不同反射角对成像贡献的相对关系并没有变，并没有达到照明补偿的目的，应用本文方法照明补偿后，不同反射角的成像贡献更加的均衡，如图 9b矩形框内振幅曲线可以看出，照明补偿后的振幅值局部范围内最高值和最低值的差最小，相对更加均衡.在图 9a中选了6个点进行局部倾角域照明补偿前后的振幅随角度(AVA)分析，图 10a和图 10b分别是这6个点补偿前后的局部倾角域成像振幅值，照明补偿后成像振幅值在各个角度上的均衡性相对于补偿前有着显著的改善.

3 结论

针对起伏地表条件下的有限采集孔径和深层弱照明问题，本文提出了一种基于小波束的起伏地表照明补偿方法，解决了起伏面对波场延拓的影响，也补偿由于采集孔径不足造成的深层成像能量不均衡问题，为后续的AVA分析提供更为可靠的偏移数据.在偏移过程中，基于波场逐步累加的思路，使得波场延拓能够适应起伏的地表条件，引入LCB小波束偏移算子，在每一层上将背景速度窗化，降低了背景速度扰动量，提高了成像精度，克服了传统单程波算子全局参考速度导致的起伏面附近强成像噪音问题.局部倾角域的照明补偿不但提升了深层照明度，而且成像振幅均衡度也得到了有效改善，有助于后期的属性分析和反演.对复杂的SEG起伏地表模型测试说明了该方法的正确性.在近地速度精度较高情况下，本方法能够适用于地表高程差较大和地表横向起伏变化剧烈地区的成像，并有望在起伏地表地区的实际资料处理中取得较好的效果，为山地、盆山边界等复杂地表地区地震勘探发挥重要作用.