2019622039陈鑫地震成像理论共28页

地震成像定理
2019622039 陈鑫
一、综述将双程波动算子分解成单程波动算子的方法和途径； 1.褶积算子分解法
用褶积算子分解，具体如下：
在非均匀介质情况下全弹性双程延拓矢量方程写为
s
s s Q A Q z
∂=∂，式中s A 是一个6⨯6的矩阵，考虑二维情况，用空间褶积代替空间求导，并假设介质λ，μ，
ρ空间导数可以忽略，则s A 变为4⨯4的矩阵
1314232431324142
000000
0s a a a a A a a a a ⎛⎫ ⎪ ⎪
=
⎪
⎪⎝⎭
(1-1-1) 式中：
1301412312402310321411420()*,()*
224()*,()()()*
22()*,()*()*,()*
iw a d x a d x a d x a iw d x d x iw a iw d x a d x iw
a d x a d x λ
λμλμ
λμλμ
ρλμλμ
ρμ
==+++=-=-++==-=-=
（1-1-2）
注，*代表褶积。

上式不能应用到流体中，因为流体介质的0μ=，而（1-1-2）式中μ为42a 的分母。

算子s A 可分解为
1
s s S A L L -=Λ （1-1-3）
式中s L ，s Λ和1s L -的表达式将在下面给出，为了使物理意义更清楚一些，先要说明，在二维情况下，如果波前的法线位于垂直平面内，则在该平面（,x z
平面）极化的波是纵波P 波和SV 横波。

此时0y V =，0zy σ=，与定义式相符，
波矢量包含P 波和SV 横波。

对水平极化的SH 横波，有0x z V V ==，0zx xz σσ==，此时波矢量将更为简化，形式上与声学方程相似，在此暂不讨论。

对质点速度而言，还可用一个标量位的梯度和一个矢量位的旋度来表示，即V θϕ=∇+∇⨯,式中θ为标量位或胀缩位，ϕ为矢量位或剪切位。

对于均匀介质，θ与ϕ彼此独立，可分别导出纵、横波方程。

在非均匀介质中，θ和ϕ相互耦合。

设s W ϕ=∇⨯，对于二维,x z 平面，选s W 的,x z 分量为零，s s y W W =，只有y 分量存在。

1
12
2s L
L L L L ⎛⎫
= ⎪-⎝⎭
（1-1-4）
00⎛⎫ΛΛ= ⎪ ⎪-Λ⎝⎭
（1-1-5）
111
12111212s
L L L L L -----⎛⎫
= ⎪-⎝⎭
（1-1-6） 其中1L ，2L ，∧，0，1L -，12L -都是2⨯2的矩阵，它们的表达式分别为：
11111
02*
()*
2()*[()2()]*P P iH d x L d x H iw d x i d x w w μμρ⎛⎫
- ⎪= ⎪--+ ⎪
⎝⎭
（1-1-7） 0211211[()2()]*2()**()**s
s iw d x i d x d x H L w w d x iH μμρ⎛⎫+-∨ ⎪= ⎪ ⎪∨⎝⎭
（1-1-8） 110
P SV iH iH ⎛⎫
-*Λ=
⎪-*⎝⎭
（1-1-9） 1
1121T
L L i ωρ
--=
Λ （1-1-10）
式中2T L 表示2L 的转置矩阵。

1
211T
L L i ωρ
-=
Λ （1-1-11）
1
10
0P SV iH iH --⎛⎫
*Λ= ⎪*⎝⎭
（1-1-12） 00000⎛⎫
= ⎪⎝⎭
（1-1-13）
()()()()()()11112
22
0211
21
1
22
02,,2,P P P P P P SV SV SV SV SV SV H H x H H H H d x d x H H x H
H
H H
d x d x δωρδλμ
ωρμ
--*=*==+*=+*==+ （1-1-14） 以上各式中上角标SV 和P 分别对应SV 横波和纵波P 。

如果定义θ-和θ+分别为上行和下行P 波位，-ψ和+ψ为上行和下行SV 波位，则
1*()
p
iH z
θθθθθθ-++-
=+∂=--∂ （1-1-15）
_1,
*()sv iH z
ϕ
ψψψϕϕ+-+∂=+=-+∂ （1-1-16）
再定义s P 为：,,,T
s P θθ++--
⎡⎤=ψ-ψ⎣⎦ （1-1-17）
可得s s s Q L P = （1-1-18） 1s s s P L Q -= （1-1-19） 容易看出s L 为合成算子，而1s L -为分解算子，将（1-1-18）带入（1-1-19）计算得
s
s s P B P z
∂=∂，式中 1
s
s
s s L B L z
-∂=∧-∂ （1-1-20） 至此就建立了非均匀介质情况下全弹性双程延拓方程与全弹性单程耦合延拓方程，利用分解算子双程弹性矢量s Q 可分解为单程耦合波矢量s P ，反之利用合成算子s P 也可合成s Q 。

2．Stolt 隐式展开法求单程波方程
地震波在地下介质中的传播，近似服从标量波动方程
22222221P P P
x z c t ∂∂∂+=∂∂∂ （1-2-1） 为了介绍动坐标变换，先考虑常速度的一维波动方程
222221P P
z c t
∂∂=∂∂ （1-2-2）
引入随时间而运动的坐标变换
z z ct
t t
'=+⎧⎨
'=⎩ （1-2-3） 由于坐标轴z '以速度c 上行（取向下为正），以速度c 上行的波在这个坐标系中便是静止的，正如两辆以速度c 并行的汽车，车内的乘客看到对方的车是静止的一样。

于是
反变换到原坐标系，得到上行波方程：
0P P
C
z t
∂∂-+=∂∂ （1-2-4）
变换（1-2-3）不能推广到二维，为此把它改写成
///z c z c t t z c t
t c tc
z z ''=+=+⎧⎧⎨
⎨''==⎩⎩即 （1-2-5） （1-2-5）和（1-2-3）类比容易发现上行波方程在该坐标系中应为 由（1-2-5）得
22222222121,P P P P P P P z z C t z z C t z C t ∂∂∂∂∂∂∂=+=++''''''∂∂∂∂∂∂∂∂ （1-2-6）
代入（1-2-2） 即得：
0P
z ∂='∂ （1-2-7） 反变换到(,)z t 同样得到（1-2-4）
对于二维方程（1-2-1）Claerbout 把（1-2-5）推广成浮动坐标系
/t z c t x x
z z '=+⎧⎪
'=⎨⎪'=⎩
（1-2-8）
将（1-2-6）和2222
P P
x x ∂∂='∂∂代入（1-2-1）得
2222
220P P P
x z C z t ∂∂∂++=''''
∂∂∂∂ （1-2-9） 这时，（1-2-7）式不成立，但22P z ∂∂相对于22P
x
∂'∂和22P C z t ∂''∂∂很小，在一定条件下可
以忽略掉。

下面我们分析22P
z
∂∂可以从（1-2-9）中略去的条件。

设x k 和z k 为沿x 轴和z 轴的波数，则有： 即，sin x k C
ω
θ=
同理，cos z k C
ω
θ-=
于是波阵面与地面夹角为θ的平面上行波 考虑到（1-2-8）式得 则
用22P z ∂'∂同22P x ∂'∂和交叉导数2P z t ∂''
∂∂相比，分别得 可见，只要地层的反射界面的倾角θ足够小时，22P
z ∂'∂可以从（1-2-9）中略去，
从而得到Claerbout 15°上行波方程。

222
02P C P
z t x ∂∂+='''
∂∂∂ （1-2-10） Stolt 在方程（1-2-9）的基础上再对z '求导，并用（9）消去22P z ∂'∂，然后略去33
P
z ∂'∂得到三阶的45°方程
3233222
042P C P C P
z t x z t x ∂∂∂-+=''
''''∂∂∂∂∂∂ （1-2-11） 3．Claerbout 显示展开法求单程波方程
假定(,,)P x z t 的傅氏变换
存在，则
()11exp()
(exp())
22x x x
P i t ik x dxdt ik i t ik x Pdxdt x ωωπ
π
∂-−−−−→-∂⎰⎰⎰⎰分部积分即P x ∂∂的傅氏变换等于P 的傅氏变换乘以x ik ，同理P t ∂∂的傅氏变换等于P 的傅氏变换乘以i ω-，这种对应关系简记为
,x ik i x t
ω∂∂
↔↔-∂∂ （1-3-1）
设C 为常数，对（1-2-1）的x t 和作傅氏变换，由（1-3-1）得 其通解为：12exp()exp()z z P C ik z C ik z =+- 其中， 21(
)x
z Ck k C
ω
ω
=-
- （1-3-2）
第一、二项两个特解分别表示上行波和下行波，上行波显然满足方程
0z dP
ik P dz
-= （1-3-3） 为了求出它的微分形式，Claerbout 将（1-3-2）用泰勒式展开取一阶近似得到15°方程，后来F.Muir 改用精度较高的连分式展开得到45°方程和60°方程等等。

下面介绍连分式展开法。

容易验证根式21R x =-可以展成连分式
2011,11n n x R R R +==-+ （1-3-4）
事实上，由2
11x R R ∞∞=-+得221R x ∞=-，即R R ∞→在（1-3-4）中分别取0n =和
1n =，并且记z
Ck R ω
=-
和x
Ck x ω
=
便得到（1-3-2）的两个展开式
2
112x z
Ck Ck ωω⎛⎫
-=- ⎪⎝⎭
（1-3-5）
2
2()12()/2
x z
x Ck Ck Ck ωωω-=-
- （1-3-6） 将（1-3-5）代入（1-3-3）得
反变换到(,,)x z t 域，由（1-3-1）得
22222
102P P C P
t z C t x ∂∂∂-+=∂∂∂∂ （1-3-7） 在浮动坐标系（1-3-7）中应用公式（1-2-6）易知上述方程就是15°方程（1-2-10）。

用同样的方法将（1-3-6）代入（1-3-3）可得到45°方程（1-2-11）。

二、对出现的文献和教科书中的成像条件进行综述，给出每个成像条件的所依赖的物理事实。

常用的成像条件有零时间成像条件、激发时间成像条件和波阻抗成像条件等。

零时间成像条件只适用于迭后逆时偏移; 波阻抗成像条件对于声波方程逆时偏移较实用, 并且当干扰波较强时, 其误差较大; 激发时间成像条件也在声波方程的叠前逆时偏移中广泛采用。

为了不失一般性并简化表达,忽略权系数时的叠前数据0()P z 可以表示为
0E
R
∂=∂ ()()()()()00,,o m m m o p z W z z R z W z z s z -++= （2-1） 式中，(),o m W z z -和()0,m W z z +分别代表声波方程源点（接收点）到空间点的单程传播；()o s z +代表源子波，()m R z 代表反射系数分布。

在深度m z 处，一次反射的正演模型为：
()()()0m m p z R z s z -+= （2-2） 故 ()()()1
m m m R z p z s z -+= （2-3）
对于多层的情况，其正演模型是：
()()()0000,m p z x z z s z -+= （2-4） 即 ()()()1,[]m m m m x z z p z s z -+-= （2-5）
忽略掉所有的传播损失，则(),m m x z z 定义为：
()()()0,,m m m m m x z z R z x z z =+ （2-6）
通过加权求和，可以得到全反射系数矩阵()m R z 。

这样，(),m m x z z 中的每一
个频率成分都变换到了线性Radon 域，并沿着恒定射线参数进行累加。

因此，为了得到更好的反射系数，要对所有的频率进行累加，即： ()()1,m m
m R z x z z N ω
=
∑ （2-7）
式中，N 是频率元素个数。

通常情况下，所谓的“零偏移距反射系数”可以适应构造成像，而且在这种情况下只需估计全反射系数矩阵()m R z 的对角线上元素的值。

()m R z 中给个元素的值都可以通过()m p z →-和()m s z →+
的商得到，即：
()()
()
,,,1
,,,,,,j k m j k m j
k
m
p x y z R x y z N
s x y z
ωωωω-∧
+
=
∑ （2-8）
式中，ω为稳定常数。

为了保证上式计算的稳定性，成像原理的简易实现，即所谓的相关形式为： ()()()1,,,,,,,,,j k m j k m j k m R x y z p x y z s x y z N
ωωω∧
-+
⎡⎤=
⎣⎦ （2-9） 称（2-9）式为成像条件1.
对于非均匀的宏观模型，源波场被严重扭曲，所以必须考虑源场的振幅变化，使用（2-8）式来成像。

由于是非均匀介质，入射波场振幅值很小的点必然存在，因此为了保证求解的稳定性，则：
()()()()2
2,,,,,,1,,,,,,j k m j k m j k m j
k
m
p x y z s x y z R x y z N s x y z ωωωωε
-+
∧
+⎡⎤⎣⎦=+ （2-10） 称（2-10）式为成像条件2。

另外一种得到对反射系数最好宽带估计的方法如下： 对于所有的频率 ()()(),,,,,,,,,j
k
m
j k m j k m p x y z
R x y z s x y z ωωω∧
-
+= （2-11）
由
0E
R
∂=∂可以得到： ()()()()
2
,,,,,,,,,,,,j
k m
j k m j k m j
k
m
p x y z s x y z R x y z s x y z ω
ω
ωωωω-
+
∧
+
⎡⎤⎣⎦
=
∑∑ （2-12）
式中，分子代表零时刻入射场和反射场的互相关，分母代表零时刻入射场的自相关。

在实际运算过程中，必须给出方程再稳定的求解方程，即： ()()()()2
2
,,,,,,,,,,,,j
k m j k m j k m j
k
m
p x y z s x y z R x y z s x y z ω
ω
ωωωωε
-
+
∧
+
⎡⎤⎣⎦
=
+∑∑ （2-13）
称（2-13）式为成像条件3.
三、从双程波动方程出发推导以地下半空间中的Kirchhoff 绕射积分和Rayleigh 积分
(一)Kirchhoff 绕射积分方程
设有一简谐波，其圆频率为ω
将上式代入均匀介质情况下的声学波动方程
2222
1t
u
V u ∂∂=∇ （3-1）
式中波场()t z y x u u ,,,=，V 为波的长波速度，介质均匀时V 为常数。

代入后得到以下形式的亥姆霍兹方程
022=+∇u K u （3-2） 式中V K /ω=为圆波数。

假定在任意观测点P 的周围，任取一个闭合曲面S ，用W 表示闭合面所包含的体积，又设在S 面内和S 面上，u 及其各一阶导数，二阶导数都是连续的，G 是具有同样连续性的另一任意函数，按格林定理得
()
ds n u G n G u dw u G G u s W ⎰⎰⎰⎰⎰⎪⎭⎫ ⎝
⎛∂∂-∂∂-=∇-∇22 （3-3）
式中n ∂∂表示沿着闭合曲面S 向里的法线方向求导。

若G 也满足方程（3-2） 022=+∇G K G （3-4） 由（3-2）和（3-4）得 两式相减消除u GK 2得
022=∇-∇u G G u （3-5） 则有
0=⎪⎭⎫ ⎝⎛∂∂-∂∂⎰⎰dS n u G n G
u s (3-6)
令 ()r e z y x G ikr -=,, （3-7） 式中r 表示P 点到任意点()z y x ,,的距离。

根据（3-7）式，在0=r 处，函数G 有一个奇点，又因G 假定是连续可微的，P 点必须从积分域中挖掉，用半径为ε的小球面包围P 点并且开拓积分到表面S 和小球面S '之间的全部体积得
0='⎥⎦⎤
⎢⎣⎡∂∂-⎪⎪⎭⎫ ⎝⎛∂∂+⎥⎦⎤⎢⎣⎡∂∂-⎪⎪⎭⎫ ⎝⎛∂∂⎰⎰⎰⎰'----S d n u r e r e n u dS n u r e r e n u S ikr ikr S ikr ikr （3-8） 在小球面S '上，法线n
的方向与r 的方向一致，n ∂∂可换成r ∂∂，于是有
S d r u r e r iK r e u dS n u r e r e n u S ikr ikr S ikr ikr '⎥⎦⎤
⎢⎣
⎡∂∂-⎪⎭⎫ ⎝⎛---=⎥⎦⎤⎢⎣⎡∂∂-⎪⎪⎭⎫ ⎝⎛∂∂⎰⎰⎰⎰'----1 （3-9） 在小球面S '上，ε=r ，并且有Ω='d S d 2ε式中Ωd 为立体角，得
Ω⎥⎦⎤⎢⎣⎡∂∂-⎪⎭⎫ ⎝
⎛
---=⎥⎦⎤⎢⎣⎡∂∂-⎪⎪⎭
⎫ ⎝⎛∂∂⎰⎰⎰⎰----d r u e iK e u dS n u r e r e n u ikr ikr ikr ikr 21εεεε （3-11）
当包围P 点的小球无线缩小，即球半径0→ε时，（3-11）式右边一、三项为0，第三项的极限为
()dS n u r e r e n u P u S iKr iKr ⎰⎰⎥⎦
⎤
⎢⎣
⎡∂∂-⎪⎪⎭⎫ ⎝⎛∂∂=
--π
41
（3-12） 此式称为简谐波的克希霍夫积分式。

假设时间空间函数(,,,)u x y z t 是波动方程的一个解，对其作傅氏变换，
(,,,)u x y z ω，显然，u 满足方程（3-2），则
()()()dS n z y x u r e r e n z y x u P u S iKr iKr ⎰⎰⎥⎦
⎤
⎢⎣
⎡∂∂-⎪⎪⎭⎫ ⎝⎛∂∂=
--ωωπ
ω,,,,,,41
, （3-13） (,)u P ω相对ω做傅氏反变换得
()()⎰⎰⎰⎥⎦
⎤⎢⎣⎡∂∂-⎪⎪⎭⎫ ⎝⎛∂∂=
--∞
∞-s t
i ikr ikr d dSe n z y x u r e r e n z y x u ωωωππ
ω,,,,,,4121
(3-14) 交换积分次序，并用V ω代替K 得：
令[]()()ωωπ
ωd e z y x u V r t z y x u u V r t i -∞
∞
-⎰
=⎪⎭⎫ ⎝⎛
-
=,,,21,,,，并考虑到ωi 与
t
∂∂
的傅氏变换对应，则
()[]dS n u r t u n r Vr r n u t P u S ⎰⎰⎭⎬⎫
⎩
⎨⎧⎥⎦⎤⎢⎣⎡∂∂-⎥⎦⎤⎢⎣⎡∂∂∂∂-⎪⎭⎫ ⎝⎛∂∂=
11141,π
（3-15） 上式为一般情况下的克希霍夫积分公式。

（二）Rayleigh 积分
设u 为压力，n q 为法向质点速度，根据牛顿第二定律，单位体积m ∆的运动方程为
式中ρ为介质密度。

将其代入(3-12)
对()t P u ,和n q 相对t 做傅氏变换并考虑到t ∂∂和ωi 的对应关系得：
()dS r e q i r e n u P u S iKr n iKr
⎰⎰⎥⎦⎤
⎢⎣
⎡+⎪⎪⎭⎫ ⎝⎛∂∂=
-ωρπω41
, （3-16）
设闭合面由上半空间（0Z <）的半球面2S 和平面1S （0Z =）组成，这时
121112(,)[()()]4iKr iKr
n n S S i e e u P q H dS q H dS r r ωρωπ--=+++⎰⎰⎰⎰ （3-17）
设上半空间的波场值是由下半空间（0Z >）中的次生源产生的，则在计算上半空间P 点在有限时间间隔0MAX t T ≤≤内的波场值时，总可以找到R 使来自2S 面上的贡献在时间小于MAX T 时尚未达到P 点，即对给定的MAX T ，总可以找到球半径0MAX R VT =，使得0R R >的任意点有：
则111(,)()4iKr
n S i e u P q H dS r ωρωπ-=+⎰⎰ （3-18）
1H 在1S 上合1S 内部满足方程（2）
，同时在S 上满足 1()iKr
e H n r
-∂+∂=0 1S 上任一点(,,)x y z 至P 点的距离为222()()()P P P r x x y y z z =-+-+-，
则2(
)(1)/iKr
iKr P z z e e iKr r n r r
--+∂=---∂ (3-18) 对1S 面取P 点的镜像P '，222
()()()P P P r x x y y z z '=-+-+-，
取1i K r e H r '
-
='
则 2211()(1)/cos (1)/iKr iKr P H z z r H e iKr r e iKr r n r n r ϕ''--'
∂-∂∂''''==---=-+''
∂∂∂(3-19) 则 1
1(,)2iKr
n S i e u P q dS r ωρ
ωπ
-=
⎰⎰ （3-20） 此为半空间的克希霍夫积分，又称为Rayleigh I 积分。

1
12
1
1(,)cos 2iKr S iKr
u P u
e dS r
ωϕπ
-+=
⎰⎰ （3-21） 此为半空间压力场的克希霍夫积分公式，称为Rayleigh II 积分。

四、从双程波动方程出发推导出程函方程和射线方程
1. 程函方程的推导
程函方程在数学上是用WKBJ 方法在某种近似条件下把二阶的波动方程近似为一阶方程；在物理上是用集合光学射线理论近似波动理论。

近似的条件是波速的空间梯度不大。

程函方程在几何光学上就是光程函数所满足的方程。

当波速为常数时，波动方程
具有以下形式的解
()cos cos cos f x y z ψαβγ=++ （4-1-1）
若(),,x y z 是传播方向()0cos ,cos ,cos l αβγ=上的一点，波从原点传播到
(),,x y z 的时间（简称为走时）是
()cos ,cos ,cos /x y z c ταβγ=，显然τ满足方程
2
2
2
21x y z c
τττ⎛⎫∂∂∂⎛⎫⎛⎫
++= ⎪ ⎪ ⎪∂∂∂⎝⎭⎝⎭⎝⎭ （4-1-2）
当速度c 可变时，假定在有限时间间隔[]0,T 内波动方程的解为(),,,x y z t ψ，把它展成Fourier 级数
式中i A 是圆频率为i ω的谐波的振幅，i φ是谐波i 的相位，i τ为谐波i 的走时。

将上式代入波动方程其实部为 若
2
2
i i i
A A c ω∆则：1c τ∇= （4-1-3）
称之为程函方程，它表明走时的空间变化率的绝对值等价于介质的慢度。

(),,i x y z const τ=的三维空间曲面为谐波i ω的波阵面或波前，波阵面的法线就是
射线，于是空间每一点(),,x y z 的射线方向为
,,i i i
i x y z τττ
τ∆
⎛⎫∂∂∂−−→∇ ⎪∂∂∂⎝⎭
（4-1-4） 波通过空间任意两点A 和B 的走时为 其中θ为积分路径ds 和射线方向的夹角。

由（4-1-4）看出τ∇的方向为射线的方向，因此（4-1-3）为
01l c τ∇=或01
l c
τ∇=
其中0l 为射线的单位向量。

于是沿射线方向求导数得 这样得到程函方程的另一形式
011d l dl C c ⎛⎫⎛⎫
=∇ ⎪ ⎪⎝⎭⎝⎭
（4-1-5） 此式有称为几何光学的镜像方程式，其中l 为仿射参数。

若记(),,r x y z =，l 取为弧长s 由 得：
0dr
l ds
= 故（4-1-5）可写成：
11d dr ds c ds c ⎛⎫⎛⎫
=∇ ⎪ ⎪⎝⎭
⎝⎭ （4-1-6）
2. 射线方程的推导
把地球划分为矩形网络并且编号之后，给定每个网络元中心处一个慢度值
i ω，速度梯度λ→
可以由相邻网格中心慢度的差商算出，如 其中x ∆为网格水平步长。

射线长度s 的函数的射线位置
()''''''''
0000''()1()0(0)()(0)()s s s c r s r s r ds c r s ds ds c r c r s l →→→→→→→⎛⎫⎛⎫
⎪ ⎪⎛⎫⎝⎭ ⎪=++∇ ⎪⎛⎫⎛⎫⎝⎭ ⎪ ⎪
⎪ ⎪⎝⎭⎝
⎭⎝⎭⎰⎰⎰ （4-2-1） 其中()0r →为震源位置，0(0)l →
为入射线在0s =时的单位向量。

为了求解（4-2-1）式，考虑在一个网格区域内，速度的梯度为常数即 于是，
00c r c r r r λ→→→→→⎛⎫⎛⎫⎛⎫
=+- ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭
（4-2-2）
为简单起见，设00r →
=，并记()()''00,,c r c c r s c c r s c →→→⎛⎫⎛⎫⎛⎫
=== ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭。

对 01c r
λ→→
+按Taylor 级数展开并取第二项得
20000111r r c c c c r λλλ→→→→→→⎡⎤
⎛⎫⎢⎥
⎪=-+⎢⎥ ⎪+⎢⎥⎝⎭⎣⎦
（4-2-3） 为了计算（4-2-1）式中的第二项和第三项，首先计算
()()'''2'''0''0002
s s s d r s s r s ds ds l ds →→
→⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦
⎰⎰⎰ （4-2-4）
1d r r d r
→→
→
∇==
（4-2-5） 整理得：
()'
''20000
0012s s l c ds r s s s c c c λλ
→→→
→⎡⎤⎢⎥=+=+⎢⎥⎣⎦
⎰⎰ （4-2-6）
则：'
2
2'
'''
0''2
001......22s
s
s s c ds ds r l s c c c λλλλ→
→
→→→→⎡⎤
⎛⎫⎛⎫∇==--- ⎪ ⎪⎢⎥⎝⎭⎝⎭
⎣
⎦⎰⎰
（4-2-7） 将（4-2-6）、（4-2-7）带入（4-2-1）得：
()224
20002
00001224l s s
l s s l c c c λλλλλ→→→→→→
→→⎡⎤
⎛⎫⎛⎫⎢⎥ ⎪=+--- ⎪⎢⎥ ⎪⎝⎭⎝⎭
⎢⎥⎣⎦
（4-2-8） 方程（4-2-8）的s 与微分得到射线穿出网络的方向 ()22
200000001l d r s
l s
l s l ds c c λλλλ→→→
→→→
→
→→⎡⎤
⎛⎫⎛⎫⎢⎥ ⎪=+-
- ⎪⎢⎥ ⎪⎝⎭⎝⎭
⎢⎥⎣⎦
（4-2-9） 旅行时间为：
2220202000126l s s s l c c c λλλ→→→
→→⎛⎫⎛⎫⎛⎫ ⎪ ⎪⎛⎫⎝⎭ ⎪ ⎪=-+- ⎪ ⎪⎝⎭ ⎪
⎪ ⎪⎝⎭⎝⎭
（4-2-10） 通过一个网格的路程为：
()()()()0
00t s
dt
s V s dt V ds V t s ds
===⎰⎰
（4-2-11） 五、综述利用有限差分法解波动方程的基本原理并详细给出你所了解的两种有限差分格式。

基本原理：
有限差分法就是一种用数值表给出定解问题的方法，其主要思想就是用差商代替微商。

当x ∆很小时，
()()()
du x u x x u x dx x
+∆-≈
∆ 对于偏微分形式，写成差商形式与常微分形式是类似的。

因而，对于一个偏微分方程，我们可以通过这种近似把它写成差分方程的形式。

当然，这种近似所带来的误差总是存在的，但只要它所产生的误差不要超出我们想要限制的范围，那么这种方法对我们所讨论的问题就是可行的。

两种波动方程的有限差分格式：
差分方程的格式基本可分两大类，即显示格式和隐式格式。

在实际工作又可以衍生出很多种格式，甚至可以用显示与隐式联合形式的差分格式。

在这里我们仅以抛物型偏微分方程为例说明常用的显示格式和Grank-Nicolson 格式。

（1）显示格式
设有抛物型偏微分方程：
220u u
t x
σ∂∂-=∂∂ （5-1）
具体的差分方程为： (),1,1,,1,2
2i j i j i j
i j i j t
u u u
u u x σ++-∆=+
-+∆
（5-2） 所用的差分网格如图（5-1）所示。

令
2
t
r x
σ∆=∆，则（5-2）式可以写为：
(),1,1,,1,2i j i j i j i j i j u u r u u u +-+=+-+ （5-3） 其中,i j u ，1,i j u -和1,i j u +为已知值。

从第j 时间层上的已知值，可用（5-3）式直接计算第1j +时间层的值。

以此类推，可解出全时间上的物理量值。

（2）Grank-Nicolson 隐式格式
虽然显示法计算上很简单，但它有一个严重的缺点，即时间步长t ∆一定要很小，必须满足2
1
02
t
x σ∆<
≤
∆，才能保持计算上稳定并达到必要的精度。

Grank- Nicolson 在1947年提出了一个使所有的有限r 值都满足计算要求（收敛性与稳定性）的方法。

他们把22u x ∂∂项用第j 和1j +时间行上的平均差分式来逼近。

求出下列差分方程：
,1,1,1,11,1
1,,11,2
2
222i j i j
i j i j i j i j i j i j u u u u u u u u t
x x σ++++-+++---+-+⎧⎫
=
+
⎨⎬∆∆∆⎩
⎭
（5-4） 由此得出下列等式：
()()1,1,11,11,,1,2121i j i j i j i j i j i j ru r u ru ru r u ru -++++-+-++-=+-+ （5-6） 其中2
t
r x σ∆=
∆。

（5-6）式的左边包含有3个未知数值，右边3个值是已知的。

如果在每个时间层上有I 个格点，则方程（5-6）根据零时间层()0j =计算第一时间层()1j =的网格点上的函数值1i u -时可列出I 个联立方程组。

由于0j =时间层上的初值和边值是给出的，因此方程组表现为三对角线阵。

对于这样的方程组，可用代数的追赶法求解，也可以用矩阵方程来分析。

六、综述屏或广义屏近似的基本原理
广义屏传播算子理论和方法的发展过程：目前被广为应用的射线-Kirchhoff 成像方法常被称为Kirchhoff偏移方法，这是一种建立在射线理论基础上的深度偏移方法。

其局限之一是Kirchhoff方法的分辨率会随着深度增加而逐渐变差，从而导致对深部结构的分辨率降低，这一现象源于利用射线解来近似格林函数时菲涅尔带的影响。

局限之二是Kirchhoff偏移方法所携带的成像信息中缺乏正确
τ-）有限差分法。

发展了双域偏移技术，的振幅信息。

后来发展了时-空域（x
这些方法中没有离散化造成的格点频散，并且计算量比时-空域有限差分少得多，这些双域算子可以被视为超广角单向波传播算子，它们忽略了在非均匀性之间往返传播的交混回响，但是正确处理了与多次前向散射有关的波动现象。

利用不同逼近方法可以得到建立在不同近似基础上的广义屏传播算子，初期的推导使用了局部Born近似和De Wolf近似，之后单向波广义屏传播算子的导出被置于更正规的Hamilton路径积分之上，这一推导为对传播算子进行精度分析和作进一步推广奠定了数学基础。

广义屏算子的三种不同形式：拟屏近似，广义屏级数展开和Pade展开逼近。

拟屏方法基于弱散射近似，因此，垂直慢度算子可以被分解成背景部分和扰动部分，扰动部分可以从局部Born近似导出。

广义屏级数法按速度扰动的扰动度参数和光滑度参数将慢度算子垂直分量中的扰动部分展成级数，扰动度参数中的首项与局部Born近似解相当。

加入更多的高阶项，可以在提高广角精度的同时也涉及更大的计算量。

广角Pade屏方法从对垂直慢度算子进行平滑近似出发，无需借助小扰动近似，因此可以更自然地处理具有强速度反差的介质，但是，它所涉及的近似与传统方式对平方根算子进行展开时所用的局部均匀近似相当，在尖锐的间断面附近可能存在较大的误差。

进行平滑近似（即高频渐进）之后，垂直慢度算子被展开成Pade级数使用有限差分方法来进行广角校正。

七、推导均匀介质中的零炮检距F-K偏移公式并给出相应的实现流程；
二维均匀介质水平叠加厚零炮距地震坡面的偏移问题，此时，纵波方程应为
012
22222=∂∂-∂∂+∂∂t P
V z P x P （7-1） 式中()t z x P P ,,=。

设()ω,,z k P x 是()t z x P ,,的二维傅氏变化 ()()()dxdt e t z x P z k P x k t i X X ⎰⎰
∞∞-∞
∞
---=ωω,,,, （7-2）
而 ()()()ωωπ
ωd dk e z k
P t z x P X x k t i X
X -∞∞-∞
∞
-⎰⎰=
,,41,,2
（7-3）
对（7-1）式相对x,t 做二维傅氏变换，并考虑到22x ∂∂和22t ∂∂与()2
x iK 和()2
ωi 的对应关系，（7-1）式变为
P k P k V z P Z X 2
2222
2-=⎥⎦⎤⎢⎣⎡--=∂∂ω (7-4) 由于假定V 不随z 变化，在z 方向上z k 为常数，上式的解为
()z ik z ik X Z Z e C e C z k P -+=21,,ω （7-5） 式中21C C 、 为待定常数。

零炮检距情况下的偏移问题常采用爆炸反射界面成像原理，此时只考虑上行波，而上行波与（7-5）式中的z ik Z e -项对应，于是方程（7-5）简化为
()z ik X Z e C z k P -=2,,ω （7-6） 当0=z 时，得()ω,0,2X k P C =，2C 就是地面观测记录（零炮检距观测情况）或叠加后的记录 ()t z x P ,0,=的傅氏变换，这是已知的。

于是得
()()z ik X X Z e k P z k P -=ωω,0,,, （7-7） 将（7-7）代入（7-3）式
()()()ωωπωd dk e k P t z x P X z k x k t i x Z X --∞∞-∞
∞
-⎰⎰
=
,0,41,,2 (7-8)
根据爆炸反射界面成像原理，反射点存在于t=0的地方，令t=0,得偏移剖面。

()()()ωωπ
d dk
e k P z x P X z k x k i x Z X --∞
∞-∞
∞
-⎰⎰
=
,0,410,,2
（7-9）
为简单起见，将()0,,z x P 和()ω,0,x k P 分别记为()z x P ,和()ω,x k P 得 ()()()ωωπd dk e k P z x P X z k x k i x Z X +-∞
∞-∞
∞
-⎰⎰
=
,41,2
（7-10）
为了使（9）式成为真正的二维傅氏变换积分，必须将d ω变换成z dk ，将()ω,X k P 变成()Z X k k P ,。

由频散关系式： 22
2
2
X Z
k V
k -=ω 得： Z Z
X
Z dk k
k Vk d 22
+=
ω （7-11）
于是有： ()()()
⎰⎰
∞∞-+-∞
∞
-+=
Z X Z
X
Z z k x k i x dk dk k
k Vk e k P z x P Z X 22
2
,41,ωπ
（7-12）
根据频散公式，通过插值可将()ω,X k P 映射为()Z X k k B ,，这时（7-12）式变为： ()()
()⎰⎰
∞∞-+-∞
∞
-+=
Z
X z k x k i Z
X
Z x dk dk e k
k Vk k B z x P Z X 22
2
,41,ωπ
（7-13）
（7-13）式为F-K 偏移原理公式。

实现流程：
（1） 将叠加地震剖面()t x P ,做二维傅氏变换得()ω,X k P ； （2） 通过插值映射将()ω,X k P 变为()Z X k k B ,； （3） 用
22
Z
X
Z k
k Vk +乘()Z X k k B ,；
（4） 对步骤（3）的结果做二维傅氏反变换得偏移剖面(,)P x z 。

八、推导均匀介质中的Kirchhhoff 偏移公式并给出实现流程；
式中1S 代表0z =的平面，波场值(),,0,u x y t r V -定义在该平面上，(),,,P P P u x y z t 为平面外某点P 处所得的由(),,0,u x y t r V -产生的波场值，r 为平面上任一点至
P 点的距离。

针对地震资料偏移处理的实际问题，此方程应进一步改造。

首先地震勘探测线往往布置成规则的测网，设x 方向的测线与y 方向的测线垂直，则积分11
S S d ⎰⎰可改为x
y
dxdy ⎰
⎰，其次，在偏移问题中，测量数据定义在0z =的平面上，
地下不均匀点（反射或绕射）的波场为待求，待求波场在时间上应当较0z =平面上的记录的波场超前r V ，而不是像上式滞后r V ，为此须将()
,,0,u x y t r V -
改为(),,0,u x y t r V +这样上式变为 ()11,,,,,0,2P P P x y r
u x y z t u x y t dxdy z r V
π∂⎛⎫
=
+
⎪∂⎝⎭
⎰⎰ （8-1） 对于共中心点叠后记录，使用爆炸反射界面成象原理，偏移速度取真实速度之半，并在0t =时成象，得偏移剖面为 ()112,,,0,,0,2P P P x y r u x y z u x y dxdy z r V π∂⎛⎫
=
⎪∂⎝⎭
⎰⎰ （8-2） 式中()()12
2
22
P P P r x x y y z ⎡⎤=-+-+⎣⎦。

设x ∆为相邻共中心点间距，y ∆为相邻两测线间的距离，z ∆为深度采样间隔，并定义
x x α=∆， y y β=∆， z z γ=∆ （8-3） 式中,,m n l 和,,αβγ为正整数，则（8-2）的离散式为 ()12,,,,0,2r u m n l u r x y V Z αβαβπ∆⎡⎤
⎛⎫=
∆∆ ⎪⎢⎥∆⎝
⎭⎣⎦∑∑ （8-4） 式中()()()1
2
2
2
2222
r m x n y l z αβγ⎡⎤=-∆+-∆+-∆⎣⎦
，z ∆∆表示差商。

如果不考虑远场条件，先将t r V -改为t r V +，再运用爆炸反射界面成像原理（0t =成像，速度减半）得
()1
11,,,0,,02P P P x y r r u x y z u x y dxdy n r Vr n t V π
⎡∂∂∂⎤∂⎛⎫⎛⎫
=
- ⎪ ⎪⎢⎥∂∂∂⎝⎭⎝
⎭⎣⎦⎰⎰ （8-6） 令2P T z V =,02t r V =,并考虑到
2P z r r
VT r n z r
∂∂=-=-=-∂∂，()31112r VT n r z z r r r ∂∂∂∂⎛⎫⎛⎫=-=-= ⎪ ⎪∂∂∂∂⎝⎭⎝⎭ （8-7）
再利用中心差分
()()()000,,0,,,0,,,0,2u x y t t u x y t t u x y t t t
+∆--∆∂
=∂∆ （8-8）
式中t ∆为时间采样间隔，将（8-8）、（8-7）代入（8-6）
()()0002,,0,,,0,t u x y t t u x y t dxdy t ⎤∆--∆+⎥⎦
（8-9） 对（8-9）式离散求和，即得偏移结果。

九、综述逆时偏移的基本原理和实现流程；
逆时深度偏移是另一种波动理论的深度偏移方法，它的思路比较新颖。

在以
前所述各类偏移方法中，都把地表水平叠加时间剖面作为z 方向的边界条件，在
深度方向上延拓，而在逆时偏移方法中，计算是从叠加剖面的最后一个采样点开
始，向着负时间方向延拓。

另一方面，此方法还把近几年才应用的虚谱算法(也
称傅氏变换法)引进来，以解决空间导数的计算问题。

为把波动方程变化成仅适合于上行波的形式，J.Gazdag 于1981年首先提出
了90度无层间反射的单程上行波波动方程。

()()122222
2222,,22V x z V x z u u u u t x z x z ⎛⎫∂∂∂∂∂=+=+ ⎪∂∂∂∂∂⎝⎭ （9-1） 对上式的时间t 引入中心差分格式近似导数得到
即
()()()()()()122222,,1,,1,,,u x z n t u x z n t tV x z u x z n t x z ⎛⎫∂∂-∆=+∆-∆+∆ ⎪∂∂⎝⎭
（9- 2） 若设波场()(),,1u x z n t +∆和(),,u x z n t ∆已知，那么现在的问题就是求上式右边第二项的值。

一般来说，此项的值是难以求得的，但可用虚谱法得到
()(){}
122si ,,XZ Z X Z XZ F gn k i k k F u x z n t -=+∆⎡⎤⎣⎦ （9-3） 其中实现流程如下：
XZ F 和1XZ
F -分别是关于x y 、的正、反二维傅氏变换算符；()si Z gn k 为符号函数。

因此，（9-2）式可以写成
()()()()()(),,1,,1,,,u x z n t u x z n t tV x z M x z n t -∆=+∆-∆∆ （9-4）
此时若再规定延拓起始波场值(),,u x z T t +∆和(),,2u x z T t +∆（T 为地表时间
剖面的最大记录时间）为零，则可通过（9-3）和（9-4）式，完成从最大记录时
间T 到0t =的逆时延拓计算，根据爆炸反射界面成像原理，取0t =时刻的深度
剖面(),,0u x z 输出就得到了逆时深度偏移的最终偏移剖面。

需要指出的是，根据
波场值(),,u x z T 和(),,u x z T t +∆就可以计算出()T t -∆时刻的波场值
(),,u x z T t -∆，由于这个波场是向反射界面反向传播的，因此，用这种方法求得
的值(),0,u x T t -∆在物理上显然是不合理的。

但0z =处的波场值，正好可由地
表水平叠加时间剖面()0,u x t ()0t T ≤≤提供。

逆时深度偏移方法的优点是无倾角限制并允许速度在纵向和横向上变化，这
种方法也不会产生数值频散现象，较好地解决了此问题。

它是一种比较理想、很
有进一步研究价值的深度偏移方法。

实现流程：
（1） 利用双程波方程对震源波场进行正向外推，并保存外推波场；
（2） 利用逆时双程波方程对接收波场进行反向外推，每反向外推一步，应用成
像条件进行求和，得到局部成像数据体；
（3） 将所有炮集的逆时偏移结果进行叠加，得到最终的叠前深度偏移成像结
果。

十、试根据有关文献建立F-K 、Kirchhoff 及有限差分偏移之间
的关系
三种波动方程偏移方法都是建立在波场的反向外推技术之上的。

因此，外推
的算式相同，这就说明方法的原理是相同的。

做进一步的推导，分别得出三种波
动方程偏移方法的波场向下外推公式。

1、F-K 域的波场向下外推公式
()()2
222,,,exp ,,0,x y x y x y u k k z iz k k u k k v ωωω⎛⎫=-- ⎪ ⎪⎝⎭
（10-1） 2、Kirchhoff 积分法的波场向下外推公式
首先给出Rayleigh 积分公式的改造形式： 1
111(,,,)(,,0,)2p p p S r u x y z t u x y t dS z r V π∂=-∂⎰⎰ （10-2）
式子中1S 代表z=0的平面，波场值(,,0,)r u x y t V -
定义在该平面上，(,,,)p p p u x y z t 为平面外某点P 处所得到的由(,,0,)r u x y t V
-产生的波场值，r 为平面上任一点至P 点的距离。

那么根据单位脉冲函数()t δ的筛选性质可以将（10-2）式改写为： 000()1(,,,)(,,,)2p p p r t t V u x y z t u x y z t dt dxdy z r
δπ--∂=∂⎰⎰⎰ （10-3） 上式右边的三重积分实际上是个三维褶积，因此，可以写成： ()1(,,,)[(,,,)]2p p p p p p r t V u x y z t u x y z t z r δπ'
-∂=*'
∂ （10-4） 其中2
2
21/2[()]p p p r x y z z '=++-
再利用两褶积函数求导数的性质，（10-4）式可以变为
()11(,,,)(,,,)[](,,,)(,,1,)22p p p p p p p p p p p r t V u x y z t u x y z t u x y z t h x y z z t z r δππ
'-∂=*=*-'∂ （10-5） 再对上式中的,,p p x y t 做三维Fourier 变换得：
假设(,,,)p p p u x y z t 和(,,,)p p u x y z t 的三维Fourier 变换分别为(,,,)
xp yp p u k k z ω和(,,,)xp yp u k k z ω对（10-5）式两边做三维Fourier 变换，得
(,,,)(,,,)z ik z p p p p u x y z z u x y z e ωω-∆+∆= （10-6） （10-6）式与F-K 波动方程偏移以及相位移波动方程偏移中使用的延拓公式完全
一致。

3、有限差分法的波场外推公式： ()
22ˆˆˆˆˆˆˆˆˆˆˆ(,,,)(,,0,)exp ˆ2x y x y x y v u k k z u k k i z k k v ωωωω⎡⎤⎛⎫=∆-- ⎪⎢⎥⎝⎭⎣⎦ （10-7） （10-7）式与F-K 域的外推式相似，只是相位部分不等。

有限差分法中外推式的
相位是F-K 域外推式中相位的近似式。

这正是方程相位误差引起的。

可见三种外
推式是完全相当的。