第3章Kirchhoff积分法叠前深度偏移

第三章 Kirchhoff 积分法叠前深度偏移大家知道，叠前偏移的概念早在70年代中期就提出来了，但由于叠前记录的信噪比较低，偏移的初始模型又很难选准，加之当时的计算机无法承受叠前偏移较大的计算量，直到90年代叠前偏移才开始尝试应用于油气勘探地震数据的精细处理中。

常见的叠前深度偏移方法可以分为两类：第一类是基于绕射扫描叠加原理的Kirchhoff 积分法，另一类是基于波动方程的偏移方法（如有限差分偏移方法、Fourier 偏移方法等）。

本章重点讨论Kirchhoff 积分法叠前深度偏移。

Kirchhoff 积分法叠前深度偏移被认为是一种高效实用的叠前深度偏移方法，目前主要完善三维采集和叠前深度偏移软件。

积分法具有高偏移角度、无频散、占用资源少和实现效率高的特点，并且积分法能够适应变化的观测系统和起伏的地表，优化的射线追踪法和改进的有限差分法能够在速度场变化的情况下快速准确地计算绕射波和反射波旅行时，从而使积分法能够适应复杂的构造成像。

地震偏移成像问题，经过最近十多年的研究与发展，已经基本解决了和正在解决三维偏移，叠前深度偏移和多分量地震偏移等诸问题。

但是偏移中有诸多问题尚未解决，例如真振幅偏移问题和各向异性介质中的地震偏移问题。

近年来，解决真振幅偏移问题就是偏移地震数据得到真正的振幅和相位信息，从而为岩性解释服务。

由于积分法具有许多优点，因此研究Kirchhoff 型保幅叠前深度偏移具有很高的理论价值和实用价值。

下面就变速射线追踪法计算走时、有限差分法计算走时以及Kirchhoff 型常规叠前深度偏移和保幅叠前深度偏移做详细讨论和分析。

§3.1 变速射线追踪法计算走时Kirchhoff 积分法叠前深度偏移已在实际生产中应用了多年，并解决了不少复杂构造的成像问题（Zhu & Lines, 1998）。

Kirchhoff 积分法的关键是绕射旅行时的计算，目前常用的计算方法是射线追踪法和有限差分法（Schneider, 1992, 1995）。

有限差分绕射旅行时计算是基于费马原理，可在直角坐标系或球坐标系实现，具体方法原理将在本章第二节介绍。

射线追踪法计算绕射旅行时可分为常速法和变速法，常速法很简单，在此不再赘述；下面主要介绍变速法。

考虑到地下介质在纵横向上通常是变速的，为通过射线追踪较精确地生成CMP 道集中各道的反射旅行时，下面我们基于Langan （1985）的思想，推导变速介质条件下的射线追踪方程。

由程函方程可推出如下的射线方程])(1[])(1[r v ds r d r v ds d r ∇= （3-1） 其中，)(r v 是波速，r是空间位置，s 是与路径长度有关的仿射参数。

路径长度l 由（3-2）式给出⎰''=s s d r n l 0)( （3-2）其中， dsr d r n =)( （3-3）)(r n 是位置r 处指向射线传播方向的矢量，射线旅行时为s d r v r n s t s'''=⎰0)()()( （3-4） 若)(n v 为正实数且时不变，则有1)(=r n（3-5） 由射线方程，经推导得出 s d s d ])r (v 1[)r (v s d )r (v )r (v n ˆr )s (r s 0s 0r s 0000'''''∇'+''+=⎰⎰⎰''' （3-6） 其中，0r 是炮点(s=0)相对于原点的位置，0nˆ 为s=0处射线方向的单位矢量。

当射线通过恒定速度梯度或恒定慢度梯度介质传播时，可建立方程（3-6）的简单近似解。

如果射线通过恒定速度梯度介质传播，速度场可表示为r v r v ⋅+=λ*)( （3-7）其中，*v 是原点处的速度，λ 为速度梯度)]([r v r ∇。

借助于Taylor 展开()1(*<⋅v r λ)，由（3-6）、（3-3）和（3-4）分别求得())(O ]n ˆ[s v 6n ˆv 2s )]n ˆ(v 2s 1[s n ˆr )s (r 32023200020000λ+⋅λ-λ-λ-⋅λ++= （3-8）())(O ]n ˆ[s v 2n ˆv s )]n ˆ(v s 1[n ˆ)s (n 320222000000λ+⋅λ-λ-λ-⋅λ+= （3-9） ())(O }v 2s )n ˆ(]n ˆ[v 6s 1{v s )s (t 3002022020λ+⋅λ-⋅λ+λ+= （3-10） 其中，)1(11*0*0v r v v ⋅-=λ, 0v 是射线入射点处的速度。

方程（3-8）、（3-9）、和（3-10）就是恒定速度梯度介质条件下的射线追踪方程。

§3.2 有限差分法计算走时一．二维有限差分绕射旅行时计算方法在Kirchhoff 叠前深度偏移中，绕射走时的计算精度直接影响了成像精度。

为此，我们利用Eikonal 方程、费马原理和有限差分近似，基于二维介质在规则网格上进行绕射走时的有限差分计算，它无需走时的内插，并且把它扩展到三维是很简单的。

1．方形网格情况下的绕射走时有限差分计算在高频近似下，二维几何路径及二维波前的传播符合射线追踪的程函方程，()2022,z x s z t x t =⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂ （3-11） 方程（3-11）建立了走时梯度与速度模型的基本定量关系。

坐标轴是x 和z ，s 0是慢度（速度的倒数）。

方程（3-11）中的两个微分项能够应用有限差分来近似表示：()312021t t t t xx t --+∆=∂∂ （3-12a ） ()321021t t t t zz t --+∆=∂∂ （3-12b ） 其中t 0、t 1、t 2和t 3的含义如图3-1所示，x ∆和z ∆分别是横向和纵向采样间距，方形网格情况下，h z x =∆=∆。

在平面波近似下，已知t 0、t 1和t 2可以求t 3的值。

即已知 四个点中的任意三个都可以把另一点计算出来。

2122003)()(2t t hs t t --+= （3-13）为获取待求行或列上第一个点处的旅行时，利用有限差分得出如下计算公式 图3-1 方形网格走时计算示意图2122003)(25.0)(t t hs t t --+= （3-14）其中，3210,,,t t t t 分别是内行或内列上旅行时的相对极小值、两侧的旅行时、和待求行或列上的旅行时。

2．矩形网格情况下的绕射走时有限差分计算为提高绕射走时有限差分计算方法的适应性，考虑具有不同纵、横向采样间距的矩形网格。

矩形网格走时计算见图3-2a 。

基于方程（3-11）和（3-12），经推导可以得出：经过（3-15）式，“十”字形的四个顶点值都可以求取出来，这样十字形的各个顶点构成一个环，下面要做的是如何把该环进行外推，求取整个模型上的走时。

我们把距离前一个环上的极小和极大走时点最近的点定义为该环的相对极小和极大走时点。

对于沿z 方向边界上的 图3-2b ‘十’字型中, O 点是震源点的位置, ‘十’字型各点的走时已知, 各顶点的走时可由相邻三个点的走时求出。

O A B C D图3-2a 矩形网格中水平距离为x ∆, 垂直距离为z ∆,慢度s 0为常数, t 0、t 1和t 2的走时已知, 可以求t 3的走时值。

S 0t 3 t 2 x ∆ z ∆t 0、 t 1 （3-15） 22221202212022102)1(3)()(2)()(z x t t s z x z x t t t z t t t x t a ∆+∆--∆+∆∆∆+-+∆+-+∆=t 0、 t 1 z ∆S 0x ∆ t 2 t 3第一个点可由方程（3-11）的非中心有限差分求出，即方程（3-16a ）， )4)((222120)1(3z t t s x t t b ∆--∆+= （3-16a ） t 3是待求C 和D 的走时，t 0是图3-2b 中点O 的走时（相对极小值走时），t 1和t 2、分别为A 和B 点的走时。

同理可求得沿x 方向边界上的第一个点的走时为：)4)((2221200)1(3x t t s z t t b ∆--∆+= （3-16b ） 式（3-15）和（3-16）考虑了平面透射波的情况。

当计算下一环走时时，在四条边上依次进行（当震源点在地下时，震源下面的点的走时不必要计算）。

逐个从相对极小走时点出发，沿边界向前、向后逐点扫描计算走时，直到遇到极大走时点或顶点为止。

当所有四个边上的网格点的走时都计算出来后，再计算四个顶点的走时。

在扫描过程中，相对极大值被计算了两次，选其中较小的一个作为该点的走时，这等价于几何射线从两边到达同一个点，我们只考虑了初至射线。

另外，当（3-15）式右边根号下出现负值时有， ⎩⎨⎧∆+∆+=x s s t z s s t t a ),min(),min(202101)2(3 （3-17） 当（3-16）式方程的右边根号下出现负值时，我们选取： l s s t t b ∆+=),min(100)2(3 （3-18）当向z 方向外推时，z l =，当向x 方向外推时x l =。

式（3-17）和（3-18）考虑了首波出现的情况。

若考虑顶点的散射波，则有：0220)3(3s z x t t ∆+∆+= （3-19）式（3-19）考虑了散射波传播的情况。

依据费马原理，在考虑了平面透射波、首波和散射波情况下的最终初至走时为：),,min()3(3)2(3)1(33t t t t = （3-20）具体的实现步骤如下：如图3-3所示，从源点到各边的长度分别为R 1、R 2、R 3和R 4，以R 1<R 2<R 3<R 4为例，按照上面的方法一环一环地逐渐向外扩展。

首先由源点O 进行矩阵扩展直到R 1环，如图3-3虚线所示；然后基于R 1环向左、向上和向右外推直到R 2环，如图2细实线所示；再由R 2环向上和向右外推至R 3环，如图3-3中较粗实线所示；最后由R 3环向上外推至R 4环，如图3-3中最粗实线所示。

扩展矩阵与方阵所不同的是不需要整环外推，而是只对未到达边界的各边进行外推，直至各边都到达边界为止，这样整个模型上的图3-3 扩展矩阵图 R 1、R 2、R 3和R 4分别为源点O 到下边、左边、上边和右边的距离。

先由源点O进行矩阵扩展直到R 1环，如虚线所示；然后基于R 1环向左、向上和向右外推直到R 2环；再由R 2环向上和向右外推至R 3环；最后由R 3环向上外推至R 4环，这样整个模型上的走时就可计算出来。

R 1 R 2 R 3R 4 O最终初至走时（公式（3-20））就都计算出来了。