地球物理资料字处理(第二十五讲)

第二节有限差分法波动方程偏移
（Finite Difference Wave Equation Migration）
波动方程偏移是建立在射线理论基础上的绕射扫描偏移具有明显的优越性。

因此地震波在地下介质的传播遵循波动方程规律，而射线理论只是近似地描述地震波的传播规律。

因此，应用波动方程进行偏移，一方面偏移归位效果好，特别是为复杂的地质条件下反射层的偏移归位提供了可能性；另一方面在实现反射层偏移归位的过程中保持了反射波的特征。

1979年，J．F．Claerbout提出了用有限差分法直接解简化的波动方程，以实现地震波场的空间归位。

这种方法是对标量波动方程进行了差分，从波场延拓的观点获得地震偏移剖面。

因为用了差分网格，可以在垂向和横向改变速度参数，适用于一般的非均匀介质。

这种方法是首先对波动方程进行Claerbout坐标变换，将标量波动方程分为上行波和下行波两部分，这样可以避免用标量波动方程做波场延拓所带来的不适定性问题；然后忽略掉深度方向的二次偏导数项，得到变换后的简化波动方程；最后用有限差分法求解波动方程，使反射层偏移归位到反射界面空间真实位置。

在这一节首先介绍波动方程偏移的成像原理，然后介绍叠后有限差分法波动方程偏移的原理和差分方程的建立及其求解。

一．波动方程偏移的成像原理（Imaging Principle of Wave Equation Migration ）
无论那种偏移方法，包括射线理论的和波动理论的偏移方法，抽象出来都不过是把地表地震剖面转换成近似的地质剖面的一个过程，它包括两个步骤：波场延拓和成像。

波场的延拓就是已地表地震剖面为边界条件，求取地下各个点的波场值，也称为波场重建。

成像就是在得到了地下各个点的波场之后，从重建波场中检测出反映界面信息的波场值的过程。

J．F．Claerbout成像原理指出：存在于地壳里的反射界面，是处在下行波到达且刚刚出现上行波那一瞬间。

在应用这一成像原理时，针对不同的具体问题又有不同的形式。

1.爆炸反射界面成像原理
爆炸反射界面成像原理是最常用，最简单的一种成像原理。

它把地下反射界面想象成具有爆炸性的物质或爆炸源，爆炸源的形状，位置与反射界面的形状和位置一致。

它所产生的波动脉冲波，其强度，极性与界面的反射系数大小
和正负一致，因此称地下反射界面为爆炸反射面。

假定在0t =时刻，所有的爆炸反射界面同时起爆，反射上行波到达地面各个观测点，波的传播速度为介质速度的一半，即2
e v v =，如果利用波动方程式将地面测得到的波场记录作反时间方向传播，即向下延拓，则0t =时刻的波场值就正确地描述了地下反射界面的位置，自动实现了偏移成像。

此成像原理适用于对水平叠加地震剖面的偏移处理，因为水平叠加剖面相当于零炮检距剖面，自炮点发出的下行波到达反射点的路径与自该点反射返回地面的上行波的路径一致，由于做半速代换，所以叠后偏移方法就是利用延拓波场的0t =时刻的的值进行成像。

2. 测线下延成像原理
测线下延成像原理适用于有炮检距的的记录，常在叠前偏移时使用。

叠前偏移的基本思想就是以地表记录为边界条件。

在整条测线上交替的把炮点和检波点向下沿拓，直到炮点和检波器相当于埋到地下同一位置为止。

由于这是对地下每一点而言，就好像既埋置了震源又埋置了检波器一样，根据J . 的原始成像原理，显然这时要取0t =时刻的波场值输出。

测线下延成像原理也适用于零炮检距记录的偏移成像，但是对零炮检距记录使用爆炸反射界面成像原理更简单。

3. 波场传播的时间一致性成像原理
波场传播时间的一致性成像原理是弹性波偏移中的一种成像原理。

因为弹性波偏移同时要处理纵波、横波和转换波，因此它应有自己特殊的成像原理形式。

在弹性波偏移中，一般只向下延拓检波点而炮点不动，即延拓后的波场好象是炮点仍在地表位置而地下各点都埋置了检波器接收到的波场一样。

为了说明这种成像原理，以入射波（P 波）反射纵波（PP 波）和转换波（PS ）为例。

如图7-5-a 所示，入射纵波传播时间为t d 反射纵波传播时间为t RP 转换波传播时间为t RS 。

只有在反射点Q 处才有t d =t RP =t RS 。

即这时反射纵波PP 和转换波PS 在反射点Q 上结合在一起。

如图7-5（b ）
是弹性波的成像原理 其中
图中
11P
SQ Q R RP t v +=
21P S SQ Q R RS t v v =
+ P
SR d t v =
图75-波场传播时间一致性成像原理示意图
由于延拓到R 点处时，如果R 不在反射界面上，一般来说，RP RS d t t t ≠≠，只有当R 在反射界面上时，才有RP RS d t t t ==。

因此，在弹性波的延拓波场中，必须取延拓的纵波和横波（包括转换波）波场的传播时间与入射波传播时间一致的那一时间的值来构成偏移输出剖面。

二 150有限差分法波动方程的导出 (150 Finite Difference Wave Equation Derivation)
Claerbout 差分偏移是对水平叠加时间剖面，即自激自收地震记录进行的，根据爆炸反射界面成像原理，认为水平叠加剖面是位于地下反射界面上的炮点零时刻产生的脉冲波，沿界面的法向以半速度上行到地面被检波器记录道的上行波(,0,)p x z t =。

因此，地下任一点、任一时刻的波场(,,)p x z t 也都可以认为是上行波，都服从半速度的上行波方程。

偏移就是利用这种上行波方程，将地面记录的上行波(,0,)p x z t =作为边界条件，求零时刻的波场值(,,0)p x z t =也就是反射界面处的波场值，
偏移技术的关键，就是寻求地震波所遵循的上行波方程。

地震波在地下介质中的传播如果介质是均匀的，各项同性的和完全弹性的，则波场函数(,,)p x z t 满足二维纵波波动方程
22222221p
p p z
v t x ∂∂∂+=∂∂∂ (721)-- 其中v 为波的传播速度，在讨论中假设它是常数，在实际应用差分偏移方法时它可以随x 和z 而变。

对于叠加剖面而言，波动方程（7-2-1）描述的是双程波的传播，既有下行波也有上行波。

但是，根据爆炸反射界面成像原理，我们对双程时间的叠加剖面做偏移时，要对方程（7-2-1）做一些改造才能适用。

首先，要把双程时间改造为单程时间，这可通过做半速代换2v ⎛⎫ ⎪⎝⎭用代替v 来解决
，这时方程（7-2-1）就变为：
22222224p p p x z v t
∂∂∂+=∂∂∂ （7-2-2） 其次，，要做上、下行波分离，从方程（7-2-2）中去掉下行波，把它简化成仅适合上行波的形式。

为此，可对（7-2-2）做Claerbout 上行波坐标变换：
2x x z z z t t v ⎧⎪'=⎪'=⎨⎪'⎪=+⎩
(723)--
则相应的变换关系为：
22
22
x x ∂∂→'∂∂ (724)-- 22
22t t ∂∂→'
∂∂ (725)-- 2z t z z z t z z v t ''∂∂∂∂∂∂∂→⋅+=+''''
∂∂∂∂∂∂∂
222()z z z v t ∂∂∂∂→+=''∂∂∂∂22()()z t z z v t z t z v t z
''∂∂∂∂∂∂∂∂+⋅++''''''∂∂∂∂∂∂∂∂ 222
22
44z v z t v t ∂∂∂=++''''∂∂∂∂ (726)-- 因为坐标变换并不改变实际波场状态，因此存在
(,,)(,,)p x z t p x z t '''= (727)--
由以上变换公式可以得到方程(722)--在新坐标系(,,)x z t '''中的形式：
22222()04p v p p z t x z ∂∂∂++=''''
∂∂∂∂ 以后为了书写的方便把,,x z t '''仍写成,,x z t 的形式，则上式应为：
22222()04p v p p z t x z
∂∂∂++=∂∂∂∂ (728)-- 因此，对水平迭加剖面做偏移，可以归结为求解如下解问题：
222220min max 0()400(,0,)(,)
t t T x x x x z p v p p z t x z p p p p p p x t g x t ≤≥<>=⎧∂∂∂++⎪∂∂∂∂⎪⎪==⎨⎪==⎪==⎪⎩ 其中min max (,)x x 就是剖面的范围，T 为地震记录最大时间。

不难看出，此定解问题
是不适定的，因为方程中有关于z 的二阶偏导数，但在野外只能得到波场0z p =，而并不能得到0
z p
z =∂∂，所以这是一个定解条件不足的欠定问题。

为了解决这种不适定问题，Claerbout
采取了忽略方程中关于z 的二阶偏导数项的方法，从而得到了著名的15o 近似方程。

222
04p V p z t x ∂∂+=∂∂∂
（7-2-10） 因为我们进行的是时间偏移，因此要把浓度坐标z 换成用时间表示的浓度坐标
τ，这要做变换：2z V
τ= （7-2-11）
这时方程（7-2-10）变为：
222208p V p t x
τ∂∂+=∂∂∂
（7-2-12） 这就是实际中广泛使用的15o 近似方程形式。

在经过（7-2-11）变换之后，方程（7-2-8）变为
2222221082p V p p t x ττ∂∂∂++=∂∂∂∂
（7-2-13）
三 差分方程的建立及求解 (Difference Wave Equation Established & Solving )
1．150差分方程的建立(150Difference Wave Equation Established)
使用150的近似方程（7-2-12），偏移定解问题（7-2-9）可以写成
⎪⎪⎪⎩⎪⎪⎪⎨⎧=====∂∂+∂∂∂=><≥≤),(0080
02222max min
t x g p p P p p x p V t
p x x x x T t t ττ （7-2-14） 为了求解此定解问题，用差分方法建立适当的差分方程来近似微分方程，仅仅就差分结解法而言，可以用多种不同的差分格式建立差分方程。

用差分方程近似地替代偏微分方程，可分为显式方程和隐性方程两大类。

通常，显式方程求解比较容易，快速，但精度不高，对大倾角适应性差，而隐性差分方程虽然求解比较困难，比较慢，但精度高，对倾角大的反射偏移效果好。

在这里我们介绍如图7-6所示的12点差分格式建立的差分方程求解，图中点（i ，j ，k ）表示坐标,,)x t τ（中的(,,)i j k x t τ点。

(1,1,1)i j k +++表示(,,)i j k x x t t ττ+∆+∆+∆其中，k j i ,,分别表示坐标点的序号，t x ∆∆∆,,τ分别表示t x ,,τ的坐标增量。

图7-6 12点的差分格式示意图
12点差分网格的重点A 为11(,,)22i j k ++表示)2
1,21,(t t x k j i ∆+∆+ττ点。

其余12点的波场函数),,(t x p τ均在A 点用T aylor 级数展开。

在这里不介绍其详细的展开过程，感兴趣的读者可以参考有关的参考文献。

把T aylor 级数展开式中的高于四阶的小量忽略掉，并用12个点的适当的线形组合可以近似地表示方程（7-2-12）的偏导数；
t p ∂∂∂τ2=t ∆∆τ1[]),1,1()1,1,1(6
1k j i p k j i p +--++-⎩⎨⎧ -
6
1[](1,,1)(,1,)p i j k p i j k -+-+ +[]2(,1,1)(,1,)3
p i j k p i j k ++-+ []2(,,1)(,,)3
p i j k p i j k -+- []1(1,1,1)(1,1,)6p i j k p i j k ++++-++ []1(1,,1)(1,,)6p i j k p i j k ⎫-
++-+⎬⎭ （7-2-15）
t p ∂∂∂τ2=241x ∆π[]{(1,,)2(,,)(1,,)p i j k p i j k p i j k --++
+[]{(1,1,)2(,1,)(1,1,)p i j k p i j k p i j k -+-++++
+[](1,,1)2(,,1)(1,,1)p i j k p i j k p i j k -+-++++
+
[]}(1,1,1)2(,1,1)(1,1,1)p i j k p i j k p i j k -++-++++++
（7-2-16）
将（7-2-15）和（7-2-16）代入定解问题中的偏微分方程，并令：
),0,1,0(=I (1,2,1)T =--
2232v t x τα∆∆=∆，61=β ⎪⎩
⎪⎨⎧+-=),,1(),,(),,1(),,(k j i p k j i p k j i p k j x p
则得到：
[])1,1,()(+++-k j x p T I βα-[]()(,,1)I T P x j k αβ+-+
[]()(,,)I T p x j k αβ+-+=[]()(,1,)I T P x j k αβ+-+
从而得到150有限差分法偏移的差分表达式：
),1,(k j x p +=()()I T I T
αβαβ-++-[]),,()1,1,(k j x p k j x p +++ -[])1,,(+k j x p （7-1-17） 根据上述差分方程，并考虑到定解问题中给定的初使条件和边界条件，可以在整条测线上沿时间深度τ的方向从0=τ开始，由ττ∆=j j 到ττ∆+=+)1(1j j 逐层向下延拓，直到所需要的最大时间深度m ax τ为止。

而在由时间深度j τ向下一个时间深度1+j τ延拓时，则是沿时间轴t 的相反方向从最大时间m ax t t =出发，由t k t k ∆+=+)1(1到k t k t =∆直至1+=j k t τ为止，在时间深度1j τ+上取1+=j k t τ时的波场函数值为输出值，最后得到偏移后结果),,(ττx p 如图7-7所示从中可以看到，偏移时对于时间变量t 是从最大观测时间来做逆时递推的，这也体现出延拓是使波反向传播的实质。

另外，尽管我们所使用150的近似波动方程，只适用常速介质，但在有限差分计算中却可以处理速度的从纵向和横向变化，只要在一个空间差分间隔和一个延拓步长为一个常数就可以。

2．150差分方程的求解(150Difference Wave Equation Solving)
求解差分方程（7-2-18）有不同解法，可以利用褶积算法进行求解。

方程（7-2-18） 右边的系数T
I T I )()(βαβα-++-实际上是对x 的褶积因子，其分子是褶积因子，分母是反褶积因子。

为了用递推方法求出反褶积因子对这个系数用Z 变换，设
)(x ϕ=T
I T I )()(βαβα-++- 其分子为T I )(βα+-=[],12(),()αβαβαβ+-++
分母为()I T αβ+-=[](),12(),()αβαβαβ--+---
分子，分母的Z 变换分别为：
)(1Z Φ=2()[12()]()Z Z αβαβαβ++-+++
)(2Z Φ=2()[12()]()Z Z αβαβαβ--++--- =(1)()Z Z αβρρρ
--- )(2Z Φ是Z 的二次多项式，它有两个跟21,ρρ并且121ρρ⋅=。

设其中模小于
1的根为ρ
则
=Φ)(Z 22
)()](21[)()()](21[)(Z
Z Z Z βαβαβαβαβαβα---++--+++-++ ={}))(1()()](21[)(2
ρρβαρβαβαβα---+++-++Z Z Z Z （7-2-19）
上式中的分子是褶积因子的Z 变化，分母是反褶积因子的Z 变换。

令 {}β
αρβαβαβα-+++-++=2)()](21[)()(Z Z Z F （7-2-20） Z
Z F Z R ρ-=1)()( （7-2-21） 则ρ-=
Z Z R Z S )()( （7-2-22） 由（7-2-20）式中，)(Z F 表示沿x 轴的算子，
12()(,,)wT αβαβαβρρραβαβαβ
+-++=--- 在差分方程(7218)--中，右边第一项的()F Z 的褶积因子WT 作用于
[(,1,1)(,,)]p x j k p x j k +++的褶积结果为：
i F =),)(21,(ρβαβαρβαβαρβαβα-+-+--+(1,1,1)(1,,)(,1,1)(,,)(1,1,1)(1,,)p i j k p i j k p i j k p i j k p i j k p i j k -+++-⎧⎫⎪⎪+++⎨⎬⎪⎪+++++⎩⎭
=β
αβα-+ρ[]),,1()1,1,1(k j i p k j i p -+++- +
[]),,()1,1,()(21k j i p k j i p +++-+-ρβαβα （7-2-23） +[]),,1()1,1,1(k j i p k j i p +++++-+ρβ
αβα (0,1,2)i N =⋅⋅⋅ 由（7-2-21）式可知：
Z Z R Z R Z F )()()(ρ-=
或 Z Z R Z F Z R )()()(ρ-= （7-2-24）
设 2012()i i R Z r r Z r Z r Z =+++⋅⋅⋅ 2012()i i F Z f f Z f Z f Z =+++
+⋅⋅⋅ 则由（7-2-24）式得到 ：
2012i i r r Z r Z r Z +++⋅⋅⋅
=201022()()()i i i f f r Z f r Z f r Z ρρρ+++++
++⋅⋅⋅ 令上式两边的Z 幂级数系数相等，得到
00f r =
011r f r ρ+=
1-+=i i i r f r ρ
因此，差分方程（7-2-18）中，右边第一项的()P Z 的褶积因子作用于[]),,()1,1,(k j x p k j x p +++的褶积结果为：
1-+=i i i R F R ρ
上式中是求取i R 的递推形式，根据1i R -和i F 计算R 的值。

同样，根据（7-2-22）可以知道：
)()()(Z S Z Z S Z R ρ-=
或 ()()()S Z Z R Z S Z ρ=+
因此差分方程（7-2-18）中，右边的第一项的()S Z 的褶积因子作用的结果为：
i i i S R S ρ+=-1 (0,1,2)i N =⋅⋅⋅ （7-2-27）
上式是求取1i S -的递推形式，根据i S 和i R 的值计算出S 1-i 值，其中N N R S = 然后，从上面计算出的（7-2-18）式第一项i S 值中减去第二)1,,(+k j x p 值就可以求得差分方程（7-2-18）式的解，得到),1,(k j x p +的值。

上面介绍的是采用12点中心差分格式建立的差分方程及其求解过程，用差
分方程近似微分方程求解，一个很重要的问题是要求差分格式是稳定。

只有如此，才能通过差分方程求得微分方程的数值解。

在偏移中，所谓稳定性问题就是误差是否随差分方程的递推求解过程增大的问题。

如果误差不增大，或规定增的不超过某种比较慢的速率，那这种差分格式就是稳定的，经过理论证明，本节所采用的差分格式是稳定的。

150有限差分偏移方法被理论和实践证明，存在倾角的局限性，即对比较大的角度（25）的倾斜层偏移效果不好，主要表现是：较大倾角的地震反射波不能正确归位，波形特征失真和出现频散。

为了解决这个问题，在150偏移的基础上，又发展了多种适应大倾角的波动方程偏移方法。

A．Berkhour和R。

Stolt先后推导出适用更大倾角的反射界面的二阶近似450方程和三阶近似600方程，它们分别为三阶和四阶偏微分方程，用同样的方法可以推导出更高阶的近似方程。

但是利用这些近似方程的偏移要受到倾角的限制，并且由于偏微分方程阶数随着倾角的增大而提高，应用有限差分法来直接求解高阶方程的困难的。

八十年代初，马在田提出了高阶方程的偏移的分裂法，把高阶方程分裂成一组二阶微分方程，以避免高阶微商差分带来的困难。

原则上可以实现任意倾角的差分偏移，具有重要的理论意义。

但是在大倾角的偏移方法的计算量要比150有限差分偏移的计算量大的多。

另外，还有许多种大倾角偏移方法，如线性变换波动方程偏移方法，共轭偏移方法，函数代换法等等。

如图7-8a 所表示是一个实际的水平叠加地震剖面，剖面上有平直反射也有倾斜断层面反射，如图7-8b和7-8c分别使用 150有限差分法和线性变换中的方程偏移方法所做的偏移剖面。

可见在两个偏移剖面上，绕射尾巴收敛到了绕射点，但是，对于标有字母E的单面波，150偏移损失其上部分，这是由于上部比下部分更陡，而在线性变换波动方程偏移剖面上，这个断面波就完全成像了。

图（a）水平叠加剖面
图b150偏移剖面
图c 线性变换波动方程偏移剖面
图7-8 水平叠加地震剖面及其偏移结果对比图。