1 引言

反射地震勘探目的是由反射地震数据得到地下的地质构造图像，包括目标区域波场重建和应用成像条件两个步骤.波动方程描述了地震波在介质中传播所遵循的规律，为了提高效率，出现了若干对波动方程近似求解的偏移方法，主要可分为基于射线理论的克希霍夫、高斯束偏移方法以及基于波动方程波场外推的偏移方法.在复杂地质构造环境下，基于波动方程的算法会产生较精确的成像结果，因为利用该类算法对波动方程在一定频带上求解，能够自然解决多路径到达问题.按照成像条件不同，又可分为炮域成像和观测系统沉降法成像.

在反射地震勘探中，震源能量从地表到反射界面再返回到地表检波器，得到的是一个双程时间记录，而且每条地震记录都与激发源和检波器的位置有关.观测系统沉降法将地表接收到的所有炮数据同时向下延拓，每个深度得到的地震波场等效为把源和检波器布置在该深度平面上所能接收到的反射地震数据^[1].就单条时间记录来看，该偏移方法将每条地震记录放置在源处和检波器处分两次向下延拓，在到达反射界面并且源和检波器位置重合处，对应该反射的子波正好移动至零时刻，因此该偏移方法成像条件为提取零时刻零偏移距波场的值作为地层反射强度.从偏移过程可知，该方法建立在互易定理基础上.互易定理表明在检波器和源位置互换后将得到相同的地震图.物理上，互易定理成立的原因在于声波速度在同一路径上沿相反方向传播速度相同^[2].

基于观测系统沉降概念，Yilmaz与Claerbout^[3]提出了双平方根(Double Square Root)算子偏移并将其应用于横向速度缓慢变化介质情况下的成像.为了使该算法适应于横向强速度变化，Popovici^[4]和Jin等^[5]分别把裂步傅里叶和广义屏波场传播算子应用于观测系统沉降法偏移中，并分别在二维Marmousi模型和三维SEG盐丘模型上得到好的成像效果.Biondi等^[6-7]提出了共方位角(CommonAzimuth)和方位角校正(Azimuth Moveout)成像策略，为窄方位角数据深度域成像提供了有效方法.孙沛勇等^[8]通过增加屏的阶数，提高了横向剧烈变速条件下DSR广义屏叠前深度偏移的精度.程玖兵等^[9-10]对比分析了几种基于稳相近似的DSR偏移算法并讨论了窄方位地震数据DSR方程偏移.刘定进等^[11]推导了基于DSR方程的保幅地震偏移方法.上述DSR偏移算法都是基于单程波且在中心点-偏移距坐标系下完成的，但观测系统沉降过程也可以在源-检波器坐标系下实现，而且可以用双程波算子延拓^[12-14].

小波束偏移^[15-17](Beamlet migration)方法基于小波变换和波动方程的局部扰动理论，在横向速度剧烈变化的地质环境中(如涉及到盐下成像)，空间局域化算子对介质的近似程度高，可以大大提高波的传播精度，显著提高盐下成像质量.小波束变换将频率域波场沿空间方向做小波变换，小波束偏移是将小波束变换后的波场直接在相空间(小波束空间)沿深度方向延拓.小波束偏移不仅可以得到高精度的地质结构成像结果，而且在波场延拓过程中，在小波束域(局部空间-波数域)可以提取目标结构丰富的局部方向性信息，从而为各种局域性和方向性相关的分析研究提供有力工具^[18-19].

大部分的小波束偏移方法都在炮域实现，Luo等^[20]探索了偏移距域小波束叠前深度偏移算法.本文将局部余弦基小波束偏移^[21]与观测系统沉降法结合，在源-检波器坐标系下对频率域波场进行局部余弦小波束分解，分解后波场沿共源小波束和共检波器小波束分别向下层深度延拓.本文第二部分介绍了局部余弦基的定义及其对叠前数据的分解.第三部分推导了在源-检波器坐标系下的局部余弦基观测系统沉降法传播算子.第四部分通过二维SEG/EAGE和Marmousi模型算例验证该方法的正确性以及在盐丘模型和复杂介质中的成像效果.第五部分分析对比了观测系统沉降法和炮域偏移计算效率.第六部分是结论.

2 局部余弦基及其对频率域源-检波器坐标系下叠前地震数据的分解 2.1 局部余弦基

具有局域化特性的基可以用来分解在不同时间或空间区间上有不同行为的信号，Coifman和Meyer^[22]在20世纪90年代初提出的局部余弦变换(Local Cosine Transform，LCT)就是其中的一种.局部余弦变换中的基函数是由平滑、紧支撑的钟函数与一个余弦函数的乘积构成，这些基保持了正交性并且具有较小的Heisenberg乘积.局部余弦基与加窗傅里叶变换有很多相似处，但Balian-Low定理证明窗口调制的指数函数不能构成Heisenberg乘积为有限值的正交基，也即加窗傅里叶变换在原空间和映射空间不能同时获得好的局域化性质.

钟函数起到把一个长信号划分为若干与钟函数长度相同的不同区间短信号的作用，且有如下性质：

(1)

其中，x_n表示第n个窗的起始位置，ε和ε′为区间左、右重叠半径，重叠半径起到平滑相邻区间信号的作用，避免产生“块效应”，该性质同时保证了局部余弦基的正交性和提供快速算法.图 1给出了一种钟函数形状及其对空间的划分.

局部余弦基函数表达式定义如下：

(2)

其中，L_n=x_n+1-x_n为区间长度，m为波数指数，B_n(x)为窗中心在x_n处的钟函数.局部余弦基在波数域的表达式如下：

(3)

其中，B_n(ξ)为钟函数B_n(x)的傅里叶变换，且ξ_m=.

局部余弦变换及其反变换可以用快速算法实现^[23]，可以达到O(N_xlog₂N_w)的计算复杂度，其中N_x和N_w为信号和钟函数采样点数.在N_x较大情况下，局部余弦变换与快速傅里叶变换(计算复杂度为O(N_xlog₂N_x))计算效率将会有显著的差异.

2.2 频率域源-检波器坐标系下叠前地震数据分解与重构

在二维地震数据采集中，每条地震数据由产生该数据时源的位置x^s，检波器接收位置x^r和记录时间t决定.因此地表的反射地震记录是一个三维变量数据体u(x^s，x^r，t).定义从时间域到频率域的傅里叶变换如下：

(x^s，x^r，ω)表示频率域波场，ω表示时间角频率，为了简便，后文中在不引起混淆时用(x^s，x^r)表示频率域波场.

观测系统沉降法成像要求把采集到的数据一次全部读入内存，并按源和检波器的相对位置排列，构成共源道集和共检波器道集.先对共源道集进行小波束分解：

(4)

其中，U(x_q^r, ξ_q^r, x^s)为共源道集的小波束系数，x_q^r和ξ_q^r分别表示局部检波器位置和局部检波器波束.对每一个共小波束检波器(x_q^r，ξ_q^r)道集U(x_q^r，ξ_q^r，x^s)，再进行小波束分解，

(5)

经过两次分解后，地震数据从共源-共检波器道集(x^s，x^r)变换到共小波束源-共小波束检波器道集U(x_n^s，ξ_m^s，x_q^r，ξ_q^r)，小波束传播算子直接在共小波束源-共小波束检波器域内对波场延拓.

以二维SEG/EAGE盐丘模型为例，显示对地表接收到的地震数据分解后小波束系数分布.该模型共有325炮，每一炮左侧有176检波器接收，相邻炮之间有一个空间采样点间隔(50m).图 2a和图 2b分别显示的是该数据体主频(15Hz)的谱值和小波束分解后系数的绝对值.源和检波器局部余弦变换窗的长度均为16个采样点.为了清楚展示窗口内的系数分布情况，在图 2b中，我们选择两个窗口并把其系数放大显示.在偏移过程中，在同一个窗内的系数根据该窗口的局部速度向下一层延拓.每一个窗口内的系数都有一定的稀疏性，在传播时，可以只传播一些大的系数(例如10^-4乘以该层中系数绝对值的最大值)来加快偏移速度，而成像质量几乎不受影响.

频率域波场是复数，而局部余弦基为实数基，在变换时，对波场的实部和虚部分别局部余弦变换.小波束域波场可以通过下列公式恢复到空间域：

(6)

其中，为了简化表示，用μ表示小波束分解系数指数(m，n，p，q)，用b_μ^sr表示源和检波器小波束的张量积.

3 小波束域观测系统沉降法偏移成像

大多数传统的观测系统沉降法偏移在中心点-偏移距域完成，这样可以避免把数据从共源道集排列成共检波器道集，而且可以节省计算和内存使用量.然而，详细的数学分析表明，在中心点-偏移距道集实现的方法忽略了中心点和偏移距波数的交叉项，这些项在有强速度变化模型中对成像质量有巨大影响^[24].因此，考虑到横向速度变化剧烈的情况，我们选择在源-检波器坐标系下完成观测系统沉降过程.

3.1 双平方根与单平方根算子

Claerbout^[2]给出了频率域源-接收器坐标系下的观测沉降公式：

(7)

式中v^s和v^r分别代表源和检波器处的速度，当处理横向速度变化的介质时，两者的值将不相同.式(7)中的两个相加的平方根算子被称为源-检波器坐标系下的双平方根算子(DSR)，双平方根算子把源和检波器一起向下延拓.然而，把源和检波器同时延拓与把源和检波器分开依次延拓完全等价，(7)式可以等效为两个公式：

(9)

公式(8)和公式(9)都是由单平方根算子(Single Square Root)构成，(8)式将源向下延拓，而(9)式将检波器向下延拓，而且两式具有完全相同的形式，这意味着把源和检波器向下延拓的传播算子是完全相同的，而且完全由垂直波数(频率和速度)决定.利用这一性质，可以加速观测系统沉降过程.根据不同的垂直波数，可以提前计算单平方根传播算子，并存储在一个表格中.在偏移时，这些算子可以通过查表被反复多次调用，这样可以节省大量重复计算.由于小波束传播算子的稀疏性^[21]，传播算子在表格中存储时，可以利用稀疏矩阵的存储技巧存储，传播过程也变成稀疏矩阵的相乘运算.

3.2 利用局部余弦基分解单平方根算子

鉴于(8)式和(9)式具有相同的形式，下面以(9)式为例推导小波束域的观测系统沉降算子.式(9)的解为

(10)

式中，_z(x^s，x^r)为第z层波场，_z+Δz^r(x^s，x^r)是把所有的检波器向下延拓Δz后的波场.经过(10)式作用后得到的波场_z+Δz^r(x^s，x^r)(x^s，x^r)等价于源在深度z激发，而检波器在深度z+Δz接收到的波场.(10)式中算子

(11)

是(9)式的精确解.

根据在波束域传播的局部扰动理论^[21]，空间域波场被分解到波束域，而精确单平方根传播算子(11)式可以被分成局部背景传播算子P_SSR^S0和局部相位屏校正算子P_SSR^S1：

(12)

其中局部背景传播算子可以表示为

(13)

其中，v_p^r表示在第p个检波器窗口中的局部匀速参考速度，P_SSR^S1=表示局部相位校正.

小波束原子在空间中的传播必须满足波动方程，经过传播后，一个小波束原子会有不同程度的畸变，而且会扩散到其它窗口或其它小波束中.对传播后的小波束原子的重新分解就会得到小波束域传播算子：

(14)

式中，μ′=(m′，n′，p′，q′)表示经过传播后波场的小波束系数指数，〈，〉代表内积运算.小波束域背景传播算子可以表示为

(15)

由于局部余弦基的正交性，式(15)可以简化为

(16)

其中，.公式(16)表明，源-检波器坐标系下观测系统沉降法小波束传播算子与炮域单平方根传播算子完全相同^[21].只是观测系统沉降过程中从共点源变成了共小波束源，而且对检波器传播时，共小波束源不会扩散到其它小波束源.

下面用波动方程在频率波数域的解析解来推导背景介质中的小波束算子：

(17)

把波数域局部余弦基表达式(3)代入公式(17)得

(18)

如果采用归一化的局部余弦基而且保持所有的钟型函数长度和形状相同，那么位置在x_q^r处的钟函数B_q(ξ)和位置在原点处x_q^r=0钟函数B₀(ξ)有如下关系：

(19)

则背景介质中的局部余弦传播算子(18)式变为

(20)

因为小波束对局部波场的有效表示，通常背景传播算子P_SSR^S0(x_q'^r，ξ_q'^r；xrq，ξ_q^r)是高度稀疏的矩阵.而且，小波束传播矩阵与特征方程分解很相似，相对空间域算子，小波束域算子有很多对传播精度影响很小的值.在实际应用中，可以设置一个阈值对传播算子压缩，例如10^-7，或者更小，取决于对精度的要求.

3.3 频率空间域局部扰动校正

一个完整的小波束传播过程包括波数域局部参考速度的背景传播和在空间域的局部相位屏校正.鉴于这些小的相对局部参考速度的扰动很小，对于大角度在不同窗口之间传播的波，大部分的相位误差会相互抵消.所以只对窗内速度扰动进行校正就可以有很好的精度.下面给出局部相位校正过程：

(1)对每个频率，波场(x^s，x^r)乘以真实速度相位屏；

(2)沿共源道集和共检波器道集用钟型函数折叠波场(局部余弦变换的一个步骤)；

(3)再对波场乘以负的背景相位屏

(4)对波场进行反折叠(局部余弦反变换的一个步骤).

3.4 成像条件

在观测系统沉降过程中，炮点与接收点重合处零时刻的波场值就反应了该空间点的波阻抗特征^[9].叠前数据量通常很大，频率域偏移一般把单个频率波场从地表延拓到地下并得到单频像场的值，最后成像结果由所有频率结果相加(零时刻条件)得到.通过下式得到偏移成像结果：

(21)

I(x，z)表示空间域的像，real (·)表示取实部.公式(21)也可以看做在源和检波器相遇时从频率到时间域的逆变换，并提取零时刻的波场值^[2].

3.5 余弦基小波束源-检波器观测系统沉降法偏移算法

总结小波束域观测系统沉降法偏移算法如下：

(1)初始成像值I(x，z)=0；

(2)对所有的频率ω波场(x^s，x^r)，对其进行小波束分解，得到小波束系数U(x_n^s，ξ_m^s，x_q^r，ξ_q^r)；

(3)对所有的深度n=1，…，N，首先对每一个共小波束源波场向下延拓一层Δz，

其中_z+Δz^r (x_n^s，ξ_m^s，x_q'^r，ξ_q'^r)是一个中间结果，表示源在z而检波器在z+Δz所接收到的波场；

(4)对_z+Δz^r(x_n^s，ξ_m^s，x_q'^r，ξ_q'^r)沿共检波器小波束继续延拓Δz，

得到_z+Δz^r(x_n'^s，ξ_m'^s，x_q'^r，ξ_q'^r)即表示源与检波器都在深度z+Δz的波场；

(5)把波场变换到空间域进行局部相位校正；

(6)应用零时刻，零偏移距成像条件，I(x，z+Δz)=I(x，z+Δz)+real (_z+Δz^r(x^s，x^r)|_x^s=x^r=x).

4 数值算例

首先，在二维SEG/EAGE盐丘模型上测试该算法在横向强变速介质模型中的成像精度.该模型在水平方向上有1200个采样点，深度方向有150个采样点，采样间隔均为25m.其中最小速度为1524m/s，最大速度为4480m/s，速度模型如图 3a所示.图 3b和图 3c分别给出了炮域和观测系统沉降法局部余弦基小波束偏移成像结果.总体来说，在该盐丘模型上两者结果非常类似.Biondi^[25]在理论上也证明了炮域和源-检波器观测系统沉降法偏移在一定条件下是完全等价的.

接下来在二维Marmousi模型上测试源-检波器局部余弦小波束偏移算法的表现.该模型水平方向和深度方向分别有737和750个采样点，采样间隔分别为12.5m和4m.一共有240炮，每炮有96个检波器，同样为左边接收.产生数据时使用Ricker子波主频为25Hz，每条地震记录时间长度为3s，时间采样率4ms.图 4a给出了偏移时使用的速度模型.图 4b和图 4c分别是炮域和观测系统沉降法局部余弦基小波束偏移结果.仔细对比两结果发现，观测系统沉降法偏移结果较炮域偏移成像结果幅度更为均衡.虽然两种方法在理论上是等价的，但是在实际应用中，由于各种参数选取的不同，例如计算孔径和计算过程累积误差等，结果会有一定的差异，但模型中的断裂面，以及灰岩背斜等主要特征结果相对一致.

5 源-检波器观测系统沉降与炮域偏移计效率对比分析

在不考虑数据存储，传输以及可并行化程度等前提下简要对比炮域和源-检波器坐标系下观测系统沉降偏移算法的计算代价.Sava等^[26]分析了三维情况下炮域和中心点-偏移距坐标系下观测系统沉降法偏移的计算成本.在二维偏移中，炮域偏移计算量可表示为

(22)

其中N_s是炮个数，N_x为每一炮的计算孔径，N_z为计算深度，N_f为计算频率个数，两倍是因为源的场和检波器的场要分别计算.C_IC表示成像条件的计算量.假设炮之间的采样间隔为N_int，那么源-检波器观测系统沉降法计算量可表示为

(23)

N_receiver为等效检波器数量，N_x^s为观测系统沉降法每一步沉降计算孔径.炮域偏移中，检波器接收到的波场能量有可能来自炮点与检波器覆盖范围之外，例如陡倾角结构反射的能量，考虑到极端情况下深度N_zΔz处能量传播到地表扩散范围会很大，一般计算孔径N_x很大.观测系统沉降法把采集到的所有数据同时延拓，每次观测系统向下延拓深度为Δz，波场不会有较大的扩散.

炮域成像传播有效张角可以近似正比于，而沉降法中则正比于.如果保持有效张角相同，则

(24)

计算孔径N_x^s相对炮域可以小很多.当炮数量很大而且炮间隔很小，每炮检波器数量较少时，观测系统沉降法会有较大优势.

6 结论

本文将局部余弦小波束传播算子应用于观测系统沉降法偏移成像，模型算例表明该算法继承了小波束域偏移在横向速度剧烈变化情况下的成像精度.源-检波器观测系统沉降法偏移小波束算子与炮域完全相同，可以较容易地将炮域小波束传播算子应用于观测系统沉降法.在每一层对共小波束源和共小波束检波器分别依次向下延拓.观测系统沉降法把地表记录到的所有数据同时向下传播来得到地下结构的像，可以很容易的输出各种共成像道集，而小波束分解，可以提供波场和结构的局部角度相关信息，这些对模型速度分析、AVA或AVO等的研究提供了有效途径.

致谢

非常感谢谢小碧，曹军，叶月明，任浩然，陈文超关于本项工作的讨论和建议，感谢匿名评审对本文提出的修改意见.