1 引言

地震数据偏移成像是当前油气勘探开发中获取地下构造图像的主要手段，尤其是复杂构造区油气勘探开发中的构造成像(Gray et al.，2001；Yilmaz，2001；Etgen et al.，2009；Leveille et al.，2011)，但随着油气勘探开发目标的日趋复杂化和精细化，地震数据偏移成像结果不仅要满足复杂区域构造成像的要求，还要能反映地下岩性的变化，使之能在岩性研究方面希望具有类似地震数据全波形反演的功效，以满足地震数据岩性处理解释，特别是复合型复杂油气藏地震数据岩性处理解释的需要.因此，要求地震数据偏移成像具有更坚实的理论基础.

地震勘探中的现代地震数据偏移成像方法发源于20世纪70年代初期，其最初目的是通过数学处理的方法把在地表观测到的地震反射波场送回到产生它的位置上去，得到反映地下反射体位置信息的偏移成像波场(Claerbout，1971，1976，1985；马在田，1989；Robein，2010；Stolt，1986；Bleistein et al.，2000；李振春等，2011).经过40多年的发展，偏移成像方法已从叠后走向了叠前，从2维走向了3维，从时间偏移成像走向了深度偏移成像，从各向同性介质走向了各向异性介质，从标量声波方程走向了矢量弹性波方程(Du et al.，2012；Chang and McMechan，1994)，从不需迭代的常规偏移成像走向了真振幅偏移成像和基于迭代的最小二乘偏移成像(Zhang et al.，2013a；Yao and Jakubowicz，2012，2013；Dong et al.，2012；Gray，1997；孙建国，2002；崔兴福等，2004；刘伊克等，2006；张宇，2006；李振春，2008；周华敏，2014)，并尝试通过偏移实现速度和波阻抗反演(Duan et al.，2013.，Zhang et al.，2013b).地震数据偏移成像的计算实现分为两大类，1)是基于Kirchhoff积分的实现，2)是基于微分方程求解的实现.在基于积分的偏移实现中又包含有一般射线的Kirchhoff积分法(French，1974；Schneider，1978；Sun，2000)和高斯射线束的Kirchhoff方法(Červený et al.，1982；Costa et al.，1989；Hill，1990).在基于微分方程求解的偏移实现中包含有基于波场方向分解的单程波波动方程方法(Gazdag，1978；Berkhout，1979；Huang and Wu，1996；Ristow and Rühl，1994；Chen and Liu，2006；陈生昌，2007)和全波波动方程方法(Hemon，1978；McMechan，1983；Baysal，1983；Whitmore，1983).

地震数据的反射波传播包含入射波传播、反射、反射波传播三个过程，在形式上可用Berkhout(Berkhout，1982)提出的WRW概念模型表示，当前的基于微分方程求解的偏移方法都是根据上述的反射波场传播过程发展起来的，但由于没有特定的描述反射波传播的数学物理方程，因此当前的反射地震数据偏移成像方法还不是基于反射波传播方程建立的.利用微分方程求解的偏移方法的成像公式所基于的Claerbout“时间一致性”成像原理是一个根据入射波和反射波的走时之间的关系对有关地下反射点存在位置的概念性描述，Claerbout提出的成像公式所表述的反射波场与入射波场间的关系(即在反射面上反射波场等于反射系数乘以入射波场)是一种不满足反射波场传播方程的定量关系，这也是Claerbout偏移成像方法的主要错误所在，即利用Claerbout成像公式得到的偏移成像结果会存在相位误差，致使成像位置不准确.我们知道对于反射波，反射系数定义为反射波振幅除以入射波振幅(Sheriff，2002)，而不是反射波场除以入射波场.因此，我们认为应该根据描述地震波传播的数学物理方程所建立的反射波场与入射波场间的关系推导偏移成像方法中的成像公式，并依此修正Claerbout偏移成像方法中的错误.Kirchhoff积分类偏移成像方法是以绕射波场传播的惠更斯原理和绕射波场传播表达的Kirchhoff积分公式为基础，这类方法在地震数据偏移成像中仅考虑了地震波的传播，没有考虑地震波传播中入射波与速度不均匀的相互作用，而地震数据中的散射波和反射波正是由这种相互作用激发出来，因此我们认为基于绕射波传播理论的Kirchhoff积分类偏移成像方法与地震数据中的散射波场和反射波场的产生和传播过程是不匹配的.Bleistein等(2000)利用地震波传播的散射理论、拟微分算子理论和Fourier积分算子理论，并把Kirchhoff积分公式应用于反射波场，提出了改进的Kirchhoff积分偏移成像方法——真振幅Kirchhoff积分偏移成像方法，但是在他的改进方法中对于反射面上的反射波场与入射波场间的关系表达还是采用Claerbout成像公式中的表达.地震数据的偏移成像在数学上属于对线性反演问题的求解，所以应该根据地震波的正向传播方程来构建偏移成像方法.在当前的地震数据偏移成像中，对于散射地震数据和反射地震数据认为可以使用相同的偏移成像方法，而没有考虑散射波与反射波在传播上的不同，甚至直接用散射波传播公式来推导反射地震数据的偏移成像公式，这是不恰当的.在地震波传播过程中，入射波与地下不均匀体相互作用产生散射波，如果不均匀体的尺寸相对于地震波的波长很大，入射波与地下不均匀体相互作用就产生反射波.因此，散射地震数据的偏移成像方法要基于散射波的传播方程构建，反射地震数据的偏移成像方法要基于反射波的传播方程构建.

基于上述认识和当前地震数据偏移成像方法存在的问题，本文以地震波传播的散射理论为出发点，首先根据描述一次散射波正向传播的线性方程研究建立散射地震数据的偏移成像方法理论；再在高频近似条件下，对产生散射波场的地下速度扰动函数的空间变化进行近似得到地下反射率分布，把基于速度扰动函数的散射波传播方程退化为基于反射率函数的反射波传播方程，再根据描述一次反射波正向传播的线性方程研究建立反射地震数据的偏移成像方法理论.最后用模拟地震数据验证本文所提出的散射地震数据和反射地震数据偏移成像方法理论的正确性和有效性.

2 方法理论

为了不失一般性，本文仅考虑基于标量波动方程的地震数据偏移成像，即利用下述的波动方程描述地震波的传播，

式中，p(x_s，x，t)代表压力波场； v(x)为地下速度模型； f(t)表述震源函数； x=(x，y，z)；x_s 表示震源位置； ∇² 表示Laplace算子，有 . 对于给定速度模型 v(x)，方程(1)描述地震波场的非线性传播，即波场 p(x_s，x，t)与 v(x)是非线性关系.方程(1)所对应的Green函数 g(x_s，x，t)由下述方程描述，即

把方程(1)、(2)变换到频率域，分别得到

为了得到地震数据偏移成像所对应的线性正演问题，本文利用扰动方法建立地震波的传播方程.2.1 一次散射波的正向传播方程

为了描述速度不均匀性对波场传播的作用，利用扰动法对速度模型进行表达，即

式中，c(x)为假定的背景(或参考)速度模型； α(x)表示相对于背景速度模型的速度不均匀分布，也称为速度扰动函数.把表达式代入方程(3)可得到描述散射波场的传播方程，有式中，P_s(x，ω)为速度不均匀产生的波场扰动，也称为散射波场； P(x，ω)表示总波场，有 P(x，ω)=P_i(x，ω)+P_s(x，ω)，其中P_i(x，ω)表示背景速度模型下的波场，也称为背景波场或入射波场，有把方程(7)写为时间域形式，有把方程(6)右边的总波场用入射波场与散射波场的和表示，得到方程(8)的右端项为激发散射波场的源项，它是由入射波场与速度不均匀相互作用以及散射波场再与速度不均匀相互作用共同产生的.因此，方程(8)描述的散射波场传播同样也是非线性传播，即散射波场 P_s(x，ω)与α(x)是非线性关系.如果在散射波场的源项中仅考虑入射波场与速度不均匀的相互作用而不考虑散射波场再与速度不均匀的相互作用，即如果在速度不均匀处有 P_s(x，ω)≤P_i(x，ω)，在Born近似下方程(8)可写为把方程(9)写为时间域形式，有方程(9)所表示的波场传播是一种线性传播，散射波场仅是由入射波场与速度不均匀相互作用产生的一次散射波场，我们把方程(9)称为一次散射波正向传播的线性方程.在传播方程(9)中散射波场 P_s(x，ω)不仅与入射波场 P_i(x，ω)呈线性关系，与 α(x)也是线性关系.

给定背景速度模型 c(x)和速度扰动 α(x)以及震源函数 f(t)，利用方程(7-1)和方程(9-1)就可以在时间域进行一次散射波场的模拟，即

式中，x_g 为检波点位置； d_s(x_s，x_g，t)代表模拟得到的检波点位置上的一次散射波场数据.根据式(10)，一次散射波场的模拟包含有入射波场传播、入射波场的时间二阶导数与速度不均匀的相互作用产生感应源并激发出散射波场、散射波场传播三个过程.根据波动方程的线性性质，式(10)所表示的模拟计算可写为2.2 散射地震数据偏移成像方法理论

Born近似条件下得到的波场传播方程(9)、(9-1)不仅可用来进行一次散射波场模拟，同样也是散射地震数据反演和偏移成像的基础.

利用背景速度模型下的Green函数，入射波场和一次散射波场在频率域可分别表示为

令，

得到

把方程(15)写为矩阵形式，有式中，u 表示散射地震数据 D_s(x_s，x_g，ω)的矩阵； W 代表以Green函数 G(x_g，x，ω)为积分核的传播矩阵； m 为待求参数矩阵.矩阵方程(16)的最小二乘解为式中，W^-g 代表 W 的最小二乘广义逆矩阵； W^*为W 的伴随矩阵.把式(17)求得的 m和方程(12)得到的入射波场 P_i(x_s，x，ω)代入式(14)，就可得到对速度扰动 α(x)的估计，即式中Re表示取实部运算.

考虑到式(17)计算的稳定性和逆矩阵计算的复杂度以及式(18)中波场相除运算的稳定性，我们用伴随矩阵 W^* 代替式(17)中的最小二乘广义逆矩阵 W^-g，即

式(19)的运算在实质上是对波场数据 u 的伴随传播.把式(18)改写为式中 P_i^c(x_s，x，ω)为 P_i(x_s，x，ω)的复共轭.

上述式(12)、(19)和(20)所表示的运算构成了对速度扰动 α(x)的近似求解，实质就是实现了对 α(x)的偏移成像.把对 α(x)的偏移成像运算转换到时间域，有

1)地下入射波场的计算

2)地下散射波场的伴随(逆时)重建

3)对速度扰动的成像

根据式(11)所表示的一次散射波场模拟计算，则上述第一步的入射波场的计算可以改为用下述方程计算地下时间二阶导数入射波场，即

相应的上述第三步的速度扰动成像应改为

上述的公式(21)或(24)、(22)、(23)或(25)构成了本文提出的一次散射地震数据的偏移成像方法理论.

2.3 一次反射波的正向传播方程

在地震波场方法理论研究中，高频近似方法得到了广泛的应用.所谓高频近似是指在地震波的波长相对于速度不均匀性的变化尺度很小的情况下，可对地震波传播的波现象和速度扰动函数 α(x)的空间变化进行近似.在高频近似条件下，可用几何光学理论近似波动理论，前面讨论的速度不均匀体产生的散射就退化为速度分界面产生的反射和折射.在高频近似条件下，速度扰动函数 α(x)的空间变化相对于地震波的波长尺度可视为缓慢变化甚至常数.

基于上述的认识，在高频近似下，定义背景速度模型下入射波传播方向 n 的波数 ，并在波数域把速度扰动 α(x)沿入射波传播方向的导数定义为反射率函数，即

因此在高频近似下，利用式(26)的定义，散射波传播方程(6)就退化为反射波传播方程，有相应的一次散射波传播方程(9)退化为一次反射波传播方程，即方程(28)所描述的反射波场传播是一种线性传播，反射波场 P_r(x，ω)不仅与入射波场 P_i(x，ω)呈线性关系，与 r(x)也是线性关系.把方程(28)写成时间域形式，有

方程(28)、(28-1)分别是本文推导出的高频近似条件下描述一次反射波正向传播的线性方程的频率域形式和时间域形式.对于式(28)、(28-1)中的反射率函数 r(x)也可定义为速度扰动 α(x)函数沿入射波传播方向 n 的空间导数，即

由式(26)或(29)定义的反射率函数 r(x)是一个沿着速度分界面的分布函数. r(x)不同于无量纲的反射系数 r(x)，它具有长度倒数的量纲.

给定背景速度模型 c(x)和反射率分布函数 r(x)以及震源函数 f(t)，利用方程(7)和方程(28)就可以进行一次反射波场的模拟，即

式中，x_g 为检波点位置； d_r(x_s，x_g，t)代表模拟得到的检波点位置上的一次反射波场数据.根据式(30)，一次反射波场的模拟包含有入射波场传播、入射波场的时间一阶导数与反射率的相互作用产生感应源并激发出反射波场、反射波场传播三个过程.类似于式(11)，式(30)可写为2.4 反射地震数据偏移成像方法理论

2.3节在高频近似条件下得到的波场传播方程(28)、(28-1)不仅可用来进行一次反射波场模拟，同样也是反射地震数据反演和偏移成像的基础.

应用与散射地震数据偏移成像方法理论中相类似的公式推导，我们可得到对 r(x)进行成像的时间域计算公式，即

1)地下入射波场的计算

2)地下反射波场的伴随(逆时)重建

3)对反射率的成像

根据式(31)所表示的一次反射波场模拟计算，则上述第一步的入射波场的计算可以改为用下述方程计算地下时间一阶导数入射波场，即

相应的上述第三步的反射率成像应改为上述的公式(32)或(35)、(33)、(34)或(36)构成了本文提出的一次反射地震数据的偏移成像方法理论.为了对比，下述为常规的(基于Claerbout偏移成像方法理论的)反射地震数据逆时偏移成像的计算公式.

1)地下入射波场的计算

2)地下反射波场的伴随(逆时)重建

3)对反射率的成像

由对比可看出，本文提出的反射地震数据偏移成像方法计算公式与常规的反射地震数据偏移成像方法计算公式的主要差异在成像公式，即式(34)与(39)的不同.由于成像公式(34)式中含有震源波场的一阶时间导数，因此常规偏移成像与本文提出的偏移成像有90°的相位差；再者由于导数运算可提高分辨率，所以本文提出的偏移成像较常规偏移成像在分辨率方面有提高.对于反射地震数据的偏移成像，建议使用本文提出的反射地震数据偏移成像方法及其计算公式，因为Claerbout的偏移成像公式不仅存在90°的相位误差，致使成像位置不准确，而且所得到的偏移成像结果分辨率低.

对比本文提出的散射地震数据偏移成像方法计算公式(21)、(22)和(23)与常规的反射地震数据偏移成像方法计算公式(37)、(38)和(39)可看出，两者的差异也是成像公式的不同，成像公式(23)中含有震源波场的二阶时间导数，因此常规偏移成像与本文提出的散射地震数据偏移成像有180°的相位差(即反极性)；再者由于导数运算可提高分辨率，所以本文提出的散射地震数据偏移成像较应用常规偏移成像方法进行散射地震数据的偏移成像在分辨率方面有提高.对于散射地震数据的偏移成像，建议使用本文提出的散射地震数据偏移成像方法及其计算公式，因为应用常规偏移成像方法进行散射地震数据的偏移成像不仅存在180°的相位误差，出现反极性，而且分辨率低.如果用本文提出的散射偏移成像方法对反射数据进行偏移成像，则得到的反射面偏移成像结果会出现-90°的相位差；反之，如果用本文提出的反射偏移成像方法对散射数据进行偏移成像，则得到的散射体偏移成像结果会出现90°的相位差.

3 数值试验

上文的方法理论与公式推导都是三维的形式，为了数值计算的简便，我们用两个二维模型的合成数据进行试验以验证上述所提出的散射地震数据和反射地震数据偏移成像方法的正确性与有效性.为了对比，我们在试验中还给出了用式(37)、(38)、(39)所表示的常规逆时偏移成像的结果，即在模型的偏移成像试验中分别给出了常规逆时偏移成像的结果、散射偏移成像的结果和反射偏移成像的结果.由于所有的数据模拟和各种偏移成像计算都是二维形式，因此方法试验结果对比的结论是可靠的.

3.1 简单模型试验

简单模型的速度分布如图 1所示，在该模型中既包含有大小不一的9个高速、低速散射体，用以验证散射地震数据偏移成像方法的效果，又包含有1个水平反射面，用以验证反射地震数据偏移成像方法的效果.模型的水平方向尺寸5 km，深度方向尺寸2 km，网格间距都是10 m.模型中的5个高速散射体(为自左向右的1、3、5、7、9号散射体，即图中的黑点)大小自左向右分别为3个网格点、4个网格点、5个网格点、6个网格点和8个网格点，高速散射体的速度都为2500 m·s^-1，模型中的4个低速散射体(为自左向右的2、4、6、8号散射体，即图中的白点)大小自左向右分别为3个网格点、4个网格点、5个网格点和6个网格点，低速散射体的速度都为1600 m·s^-1，模型的上层速度为2000 m·s^-1，下层速度为3000 m·s^-1，速度分界面位于nz=150网格点.采用中间放炮两边接收方式进行观测，炮间距20 m，每炮401道，共250炮，道间距10 m，时间采样点数2500，采样间隔1 ms，震源采用主频为20 Hz的雷克子波.

分别应用常规逆时偏移成像方法和本文提出的散射地震数据偏移成像方法、反射地震数据偏移成像方法对图 1模型的模拟地震数据进行偏移成像.图 2、3和4分别为常规逆时偏移成像结果、散射偏移成像结果和反射偏移成像结果.

把图 2、3、4与图 1的模型比较可看出，对于散射体的成像，图 3的散射偏移成像效果最好，不仅分辨率高，而且相位(极性)、位置正确；图 2所示的散射体成像结果分辨率低，虽然成像位置准确但有180°的相位差，出现了反极性，即把1、3、5、7、9号5个高速散射体成像为低速散射体，把2、4、6、8号4个低速散射体成像为高速散射体；图 4所示的散射体成像结果存在90°的相位误差，致使散射体成像位置对应不准确.对于反射面的成像，图 4的反射偏移成像结果最好，不仅分辨率高，而且相位(极性)、位置正确；图 2所示的常规偏移成像的反射面成像结果存在90°的相位误差，致使反射面成像位置对应不准确；图 3所示的散射偏移成像的反射面成像结果有-90°的相位误差，致使反射面成像位置对应也不准确.为仔细对比各种偏移成像方法对散射体和反射面成像的分辨率、相位(极性)和成像位置，图 5为穿过第5个散射体和反射面的第232道的3种偏移成像结果的波形对比图，左边为常规逆时偏移成像的，中间为散射偏移成像的，右边为反射偏移成像的.由图 5的波形对比可看出，散射偏移成像方法对散射体的成像不仅位置准确，极性也正确，波峰反映高速散射体；反射偏移成像方法对反射面的成像不仅位置准确、分辨率高，相位(极性)也正确，波峰反映高速反射面；常规偏移成像方法对散射体和反射面的成像都存在相位上的误差，对散射体成像有180°的相位差，把高速散射体成像为低速散射体，对反射面成像有90°的相位差，致使反射面成像位置不准，而且分辨率低.这个模型的数值试验结果与本文的方法理论结果是一致的.

3.2 Sigbee2A模型试验

为了进一步验证本文提出的散射偏移成像方法和反射偏移成像方法的效果，我们用国际标准的测试偏移成像方法的Sigbee2A模型进行试验.该模型不仅含有高速盐体、精细的反射面，还包含许多高速散射体，十分适合对本文方法的验证.图 6为用于偏移成像试验的背景速度模型.模型大小为3201×1201个网格点，水平和垂直方向的网格间距均为7.62 m，该模拟数据的观测系统为海上拖缆观测系统.该地震数据含有500炮，炮间距为45.72 m，拖缆长度为7932.42 m，检波点间距为22.86，每一炮的道数不均，最少为57道，最多为348道，时间采样间隔为8 ms，时间采样点数为1500.

分别应用常规逆时偏移成像方法和本文提出的散射地震数据偏移成像方法、反射地震数据偏移成像方法对Sigbee2A模型数据进行偏移成像.图 7、8和9分别为常规逆时偏移成像结果、散射偏移成像结果和反射偏移成像结果.

对比图 7、8和9可看出，对于高速散射体的成像，散射偏移成像方法得到的散射体的偏移成像结果要好于常规偏移成像方法和反射偏移成像方法的偏移成像结果；对于反射面的成像，反射偏移成像方法得到的反射面的偏移成像结果要好于常规偏移成像方法和散射偏移成像方法的偏移成像结果.为了仔细对比3种偏移成像方法的效果，我们截取模型和3种偏移成像结果的左边部分进行放大，如图 10、11、12和13所示.图 10为真实模型的局部放大图，图 11、12和13分别为常规逆时偏移成像结果、散射偏移成像结果和反射偏移成像结果的局部放大图.对比放大图 10、11、12和13可看出，本文提出的散射偏移成像方法对高速散射体成像、反射偏移成像方法对反射面成像，均取得了理想的效果，成像结果不仅分辨率高，而且位置准确、相位(极性)正确.散射偏移成像方法有利于散射体的成像，反射偏移成像方法有利于反射面的成像.

4 结论

地震数据偏移成像方法的研究应以描述地震波传播的数学物理方程为基础，不同波的地震数据偏移成像方法需要利用不同波的传播方程.对散射体产生的散射地震数据的偏移成像应采用基于散射波传播方程建立的散射地震数据偏移成像方法，对反射面产生的反射地震数据的偏移成像应采用基于反射波传播方程建立的反射地震数据偏移成像方法.在散射波传播方程的基础上，本文利用高频近似推导出了基于反射率函数的反射波传播方程，这是对反射波传播方法理论的发展.当前的逆时偏移成像方法中的波场正向外推和反向外推都没问题，主要问题是成像公式，因为其成像公式是基于Claerbout的时间一致性成像概念导出的，而不是基于散射波场传播或反射波场传播方程导出的.当前的逆时偏移成像方法对反射面的成像存在90°的相位误差，致使反射面成像位置不准确，对散射体的成像存在180°的相位误差，致使散射体成像结果出现反极性.本文得到的成像公式是基于波场传播方程导出，对于散射和反射有不同的成像公式，散射成像得到的偏移成像结果是无量纲的速度相对扰动，反射成像得到的偏移成像结果是长度倒数量纲的反射率.对于反射面的偏移成像，本文提出的反射地震数据偏移成像方法可得到比常规偏移成像方法更好的偏移成像结果，即准确的相位(极性)、成像位置和更高的分辨率.对于散射体的偏移成像，本文提出的散射地震数据偏移成像方法可得到比应用常规偏移成像方法更好的偏移成像结果，即准确的相位(极性)和更高的分辨率.