地震偏移成像基本原理

1、地震成像原理与方法地震成像原理与方法第一章 偏移成像第二章 三维叠前深度偏移成像理论与方法第三章 共聚集点偏移第四章 共反射面叠加第五章 偏移速度建模第一章第一章 偏移成像偏移成像 1.1 偏移成像的基本原理偏移成像的基本原理1.2 波动方程偏移波动方程偏移 1.3 叠前偏移叠前偏移 1.4 偏移速度分析偏移速度分析 1.5 深度偏移深度偏移 1.6 三维偏移三维偏移 1.7 二维和三维叠前深度偏移二维和三维叠前深度偏移 地震技术的发展趋势地震技术的发展趋势:1. 三维叠前深度偏移三维叠前深度偏移(3DPSDM)-地震成像地震成像 (波动方程法波动方程法3DPSDM, CRS叠加叠加, CFP

2、偏移偏移)2. 四维地震四维地震-开发地震开发地震 (VSP技术技术, P-S技术技术, 井间地震井间地震, 3D_AVO技术技术, 4D地震地震, 弹性波阻抗反演弹性波阻抗反演, 裂缝分析裂缝分析, 岩石物理岩石物理, 地震相与地震属地震相与地震属性分析性分析, 油藏描述等油藏描述等) 三大处理技术三大处理技术:反褶积、叠加、和偏移成像反褶积、叠加、和偏移成像反褶积和叠加引自其它相关学科反褶积和叠加引自其它相关学科偏移成像基于古典技术偏移成像基于古典技术偏移成像偏移成像: 1.具有地震勘探本身的特征。具有地震勘探本身的特征。 2.计算机使其研究由地震波运动学特征计算机使其研究由地震波运动学特

3、征 过渡到地震波动力学特征过渡到地震波动力学特征 3.提高地震空间分辨率和保真度提高地震空间分辨率和保真度一偏移成像的概念一偏移成像的概念 偏移偏移反偏移反偏移反射地震方法反射地震方法: 1.激发弹性波，2.记录反射波, 3.研究地质岩层结构和物性特征。是一种反散射问题是一种反散射问题。 反射地震成像分做两步反射地震成像分做两步: 1. 记录反射波，2. 处理反射波。地震偏移技术是使地震偏移技术是使反射界面最佳成像的一种技术反射界面最佳成像的一种技术。 1.1 偏移成像的基本原理偏移成像的基本原理地震偏移地震偏移: 叠前或叠前或/和叠后偏移和叠后偏移叠前偏移叠前偏移: 使使CSP道集记录或道集

4、记录或COF道集记录中道集记录中的反射波归位的反射波归位, 绕射波收敛绕射波收敛叠后偏移叠后偏移: 基于水平叠加剖面，采用爆炸反射基于水平叠加剖面，采用爆炸反射面概念实现倾斜反射层归位和绕射波收敛面概念实现倾斜反射层归位和绕射波收敛偏移原理和偏移效果偏移原理和偏移效果见下图见下图1. 偏移成像的基本概念偏移成像的基本概念偏移原理图偏移原理图4/ )tan(2ttvdx4/ )tan(1 12/122tvtdt2/1224/ )tan(1/tantantttv偏移过程的定量分析图偏移过程的定量分析图2. 发展史发展史1).古典的偏移技术(60年代前) -反射点的空间位置成像;2).早期的计算机偏

5、移技术(6070年代)-定性和概念性地对反射波运动学特征成像;3).波动方程偏移技术(70年代后)-定性或/和定量地对反射波运动学或/和动力学特征成像.1).有限差分法有限差分法波动方程偏移：70年代初期，J.Claerbout教授首先提出了用有限差分法解单程波动方程的近似式，用地面观测的地震数据重建地震波在地下传播过程中的波场，从这些传播过程的波场中提取使地震界面成像的那些数据，组成地震偏移剖面。由于这种偏移方法在计算过程中要解波动方程或其近似式，所以被称为波动方程法偏移技术波动方程法偏移技术。 2). Kirchhoff积分法积分法波动方程法偏移：70年代中期，French和Schneid

6、er等在绕射偏移法的基础上使用了波动方程解的Kirchhoff积分公式，发展为地震偏移的波动方程积分法地震偏移的波动方程积分法。使绕射偏移建立在可靠的波的基本原理上。因而改善了偏移剖面，取得了良好的效果。3).富里叶变换法富里叶变换法波动方程法偏移：70年代后期，Stolt和Gazdag等又先后提出了在频率-波数域解波动方程，外推地震波场的方法。这种方法被称为F-K域偏移方法域偏移方法。由于该方法计算简单，效率高，因而很快得到了推广。 波动方程偏移技术的发展波动方程偏移技术的发展3. 偏移方法分类偏移方法分类二基于射线理论的叠后偏移二基于射线理论的叠后偏移与叠前偏移与叠前偏移 经典的偏移方法和

7、早期的计算机偏移方法 都是基于射线理论射线理论 经典的偏移方法只研究到达时间。只研究到达时间。叠后偏移有圆弧切线法和线段移动法；叠前偏移包括椭圆切线法和交会法等 早期的计算机偏移方法利用了波前、绕射等地震波传播的惠更斯原理，尽管只是定性的、概念性的，但与手工操作法相比偏移剖面除了归位精度除了归位精度提高外，还考虑了波形特征。提高外，还考虑了波形特征。叠后偏移有波前模糊法、绕射曲线叠加法; 叠前偏移有Rockwell偏移叠加法和Paturet-Tariel偏移叠加法等。 1. 叠后偏移叠后偏移 叠后偏移叠后偏移: 即叠加偏移，是对叠加后的地震记录做偏移。下面介绍圆弧切线法、波前模糊法和绕射曲线(

8、面)叠加法。 1). 圆弧切线法圆弧切线法 一次反射波NMO后, 得到时间叠加剖面 由 （1.1.1） 得到视深度剖面 如果界面的倾角 =0或者很小，例如只有 或更小，则视深度界面就是真深度界面。如果界面倾角不可忽略，则应当进行倾角校正，以求出反射界面的真实位置。 0*21vtz 1 校正的做法是以地面各点为圆心，以各点下至视界面的垂直距离为半径做圆弧，其圆弧族的切线即为校正后的反射界面（v=cont）。 当速度是深度的函数时，例如 为常数时，则圆弧的圆心不位于地面上，而位于地面点的正下方某深度上。这时，圆心的深度和圆弧的半径由下式求出： ),1 (0zvv)(11)(100000tvshrt

9、vchz（1.1.2） 2). 波前模糊法波前模糊法 波前模糊法也可以称为波前切线法，它是对叠加后的地震剖面进行偏移的方法。这个方法是反推反射界面上的波场。 以地面接收点为中心，把相当于反射到达时间上的值送到以 的深度为半径的圆弧上去。如果我们把深度z仍以双程时间表示，就把反射数值送到以t为半径的圆弧上去（图1-4）。把各道上的所有反射波值都按这个原则去做，并把送到同一点的值叠加起来，就可以组成偏移剖面。把某道 上某时间t上的振幅值送到相邻各道上的时间 由下式算出： zvt20 xit21224 vxttii（1.1.3） 0 xxxii其中 用波前振幅叠加来求反射界面发出的波前实际上就是用这

10、种方法做切线。 要求：较密的地震道和较高的信噪比，以得到满意的偏移剖面。 3). 绕射曲线（面）叠加法绕射曲线（面）叠加法 绕射曲线或绕射曲面叠加法是把地震剖面上的波场振幅值按绕射波时距曲线进行相加。因为绕射波时距曲线与所有反射波的时距曲线形状相比较，其凸率最大，故亦可称它为最大凸率法。 具体做法是，当要得到地震剖面上某个 点的偏移后的数据时，我们要计算一条以这点为顶点的绕射双曲线。它在各道上的时间t由下式算出： ),(00tx212204 vxttii（1.1.4） 0 xxxii式中 在进行偏移时我们把各道上等于上式时间t的波场值取出来叠加在 点的波场值上，这就算完成了 点的偏移处理，如图

11、1-5所示。 ),(00tx),(00tx 无论是波前模糊法还是绕射叠加法，其基本原理都是根据惠更斯原理提出来的。 波前模糊法是把一个道上的波场值送到各个道上去叠加输出道法输出道法；而绕射叠加法是把各道上的相应值取来在一道上叠加输入输入道法道法。两者都符合反射波归位和绕射波收敛的要求，而且它们的叠加值也相等。 波前弧或绕射曲线在x方向上的范围L称为偏移孔径偏移孔径。L的范围是由最大实际倾角来决定的。倾角越大，L越大；有效波越深，L也越大。L的大小可用下式来估算： tan20vtL（1.1.5） 孔径的中心，原则上应当位于 处，但也可以是不对称的。 0 x 图1-7是用绕射叠加偏移法处理前后的地

12、震剖面。从对比中可以看出，偏移后剖面上的地层层位关系得到了正确的反映。有利于地质解释。 0 x2叠前偏移叠前偏移 叠前偏移叠前偏移: 即偏移叠加，是对叠加前的多次覆盖的地震记录先偏移，再叠加。下面介绍椭圆切线法、Rockwell偏移叠加法和Paturet-Tariel偏移叠加法。 1）椭圆切线法）椭圆切线法 当给定CSP记录时，可用椭圆切线法(图1-8)。反射点(2D)位于以炮点和接收点为焦点的椭圆上，这个椭圆的方程可表示为： 1442222222ltvztvx（1.1.6） 对每个炮检距的记录上的反射波画好椭圆弧。做椭圆弧对每个炮检距的记录上的反射波画好椭圆弧。做椭圆弧族的切线即为偏移后的剖

13、面。族的切线即为偏移后的剖面。 2）Rockwell偏移叠加法偏移叠加法 Rockwell偏移叠加法实际上是叠后偏移所使用的波前模糊法的一个扩展。 具体做法：把每个记录道上任一t时刻的采样值，在以炮检距中点的地面点为原点的直角坐标系中送到以 为长轴， 为短轴的椭圆与各个地震记录道垂直线相交的各个点上去，并且与其它地震道送至该交点上的采样振幅值相加，即得偏移叠加剖面。 偏移叠加实质上是用振幅叠加来做切线的。偏移叠加实质上是用振幅叠加来做切线的。 2vt2)(222htv3）Paturet-Tariel偏移叠加法偏移叠加法 1971年Paturet-Tariel用相同炮检距的剖面进行叠前偏移，把所

14、有相同炮检距的偏移后的剖面叠加得到偏移叠加剖面。叠前偏移的原理如图1-9所示。 绕射点M所产生的绕射波到达时曲线为： 21212202202424vhxtvhxttx（1.1.7） 当炮检距 时，上式表现为： 0h21022204vxttx（1.1.8） 式中 为从M点到A点的双倍旅行时间。 和 的曲线表示在图1-9的右图中。 为了进行偏移，我们应当把 的曲线上的地震能量（即采样点振幅）送到零炮检距绕射双曲线的顶点M上去叠加。这样, 把各个相同炮检距的剖面偏移后叠加在一起即得偏移叠加剖面。 xt0t0 xtxt图图1-10. 偏移叠加剖面与叠后偏移剖面对比图(a). 水平叠加剖面(b). 叠后

15、偏移剖面; (c). 偏移叠加剖面三基于波动方程的波场外推与三基于波动方程的波场外推与地震成像原理地震成像原理 使用波动方程进行偏移，首先就是要重建反射波的原来重建反射波的原来波场。反射界面上刚刚产生的反射波，就认为是该反射面的波场。反射界面上刚刚产生的反射波，就认为是该反射面的像。像。 为进行波场外推，把波动方程分解为上行波方程和下行上行波方程和下行波方程。波方程。1上行波和下行波上行波和下行波 波动方程有两个解, 一般表示为 。在地震勘探中一般取深度方向向下为正取深度方向向下为正z的方向。的方向。向正z方向传播的地震波称为下行波下行波，即用 所表示的波。向负z方向传播的波为上行波上行波，即

16、用 代表的波。下行波即入射波，上行波为反射波下行波即入射波，上行波为反射波。 rvrti/)/(exprvrti/)/(exprvrti/)/(exp 只有在均匀各向同性完全弹性介质的情况下上行波和下行波才是分离的。分离过程如下:二维波动方程为： 22222221tuvzuxu（1.1.9） 对(1.1.9)式相对x和t做二维付里叶正变换，并进行算子分解得到： 0)()(22222222uikdzdikdzdukdzudukvdzudzzzx（1.1.10） 其中利用了波散关系： 2222vkkzx（1.1.11） 由（1.1.10）式得出： ukviuikdzudxz222（1.1.12）

17、其中，正号代表上行波方程，负号代表下行波方程。 2波场外推波场外推 正向外推正向外推就是根据波在当前位置上的振动情况向波的自然传播方向用计算手段预测出波场。反向外推反向外推是向波的自然传播方向的反方向上重建原来的波场。对一个波场应是进行正向外推还是反向外推均有物理问题决定。 1）上行波的外推）上行波的外推 dzkviuudx222（1.1.13） zkvixezuzzu222)()(积分结果为： （1.1.14） 由此得出上行波的正、反向外推式。 （1）上行波正向外推公式）上行波正向外推公式 上行波的正向外推式就是向负z方向的外推公式。从（1.1.14）式可求出为： zkvixezzuzu22

18、2)()(（1.1.15） 根据这个公式可以计算模拟反射波的地震记录计算模拟反射波的地震记录（地震图地震图）。（2）上行波反向外推公式）上行波反向外推公式 上行波的反向外推式就是向正z方向的外推公式。从（1.1.14）式可得出为： zkvixezuzzu222)()(（1.1.16） 根据这个公式可以进行地震记录的向下半空间延拓，求出地下任何一点的波场，实现地震波偏移的目的实现地震波偏移的目的。 2）下行波的外推）下行波的外推 dzkvidddx222（1.1.17） zkvixezdzzd222)()(积分结果为：（1.1.18） 据此可以得出下行波的正、反向外推公式。 （1）下行波正向外推

19、公式）下行波正向外推公式 下行波的正向外推式是指沿正z方向的外推。其外推式为： zkvixezdzzd222)()(（1.1.19） 这个方程可用来模拟下行波的地震记录模拟下行波的地震记录。（2）下行波反向外推公式）下行波反向外推公式 下行波的反向外推是指沿负z方向的外推。其外推式为： zkvixezzdzd222)()(（1.1.20） 上式可用来从下行波场进行反向求源的计算工作从下行波场进行反向求源的计算工作。 下面分析波场本身的条件对外推结果的影响 22xzkkk（1.1.21） 当 时， 为正或负的实数，这时所有外推公式中存在虚指数。说明在外推过程中波场发生相位变化。一般都能得出正确的

20、结果。 xkk zk当 时， 值为虚数： kkxzk22kkikxz（1.1.22） zkkxezuzzu22)()(（1.1.23） 波场外推时只有振幅变化，而无相位变化。当指数项取负号时，外推的波场迅速衰减，称这种波为倏逝波倏逝波。当指数项取正号时，外推波场迅速增大，这是一种实际不存在的波，只是进行波场计算时发生，我们称它为耗损波耗损波。在计算中要避免发生这种情况。 当 时，上行波的外推式可写为： kkxzkkxezuzzu22)()(（1.1.24） 此时反向外推遇到倏逝波，正向外推发生耗损波。分别表示为： zkkxezuzzu22)()(（1.1.25） zkkxezzuzu22)()

21、(（1.1.26） 由此可见，用上行波方程进行向下波场外推永远是计算稳定的。而用上行波方程进行正向外推就可能遇到耗损波，因此有可能是不稳定的。除非在计算中不断地把 的波场滤除掉。 同理可求出 时下行波的外推式为： kkxzkkxezdzzd22)()(（1.1.27） 此时也是反向外推遇到倏逝波，正向外推遇到耗损波。 kkx3）波场外推的）波场外推的Kirchhoff积分法积分法 Kirchhoff积分法并不直接解波动方程，而是用数学方法来描述关于波的传播的惠更斯原理，从而求出空间上任一点波场值的。 Kirchhoff早在1883年就证明了，从扰动区向外某点 传播的波的t时刻的波场 可以从扰动

22、区封闭表面上的波场 以及该波场对时间和表面法线方向的导数通过积分式求出来。因此要假定 在封闭面上和封闭面内有直至二阶导数的连续性。 Kirchhoff利用了格林定理： ),(111zyxM),(111tzyxu),(crtzyxu),(tzyxuVSdSvuvnvudVuvvu)(22（1.1.28） u取为波场函数， Rv12121212)()()(zzyyxxR, 当把观测点用包含有波前面在内的封闭曲面包围起来，如图1-11（a）,（b）那样的封闭时，这样的封闭面S和它所包围的体积V作为（1.1.28）式的积分限，经过一定的推导后得出 点的正向外推波场为： ),(zyxM SdSnuRRn

23、utunRvRtzyxu11141),(111（1.1.29） 这里 的方向取封闭表面的外法线方向。 如果把观测点M移至封闭面外，则有： n0),(111tzyxu（1.1.30） （ 1.1.29）式中 ),(vRtzyxuu （1.1.29）式就是著名的Kirchhoff积分积分。它描述了物理波场描述了物理波场传播的过程传播的过程，也满足奇次波动方程，是它的积分形式解。对我们来说，也可以称它为正向外推公式正向外推公式。 注意注意: Kirchhoff积分只满足均匀介质的情况。 是推迟场 下面讨论用Kirchhoff积分进行波场反向外推问题(地震偏移)。这时，所取的封闭体积V应在波前传播方向

24、的反方向，计算点 就在这个封闭体内。根据格林定理同样可求出形式上相同的反向外推的Kirchhoff积分式： ),(321zyxM 11111141),(SdSnuRRnutunRvRtzyxu（1.1.31） 式中的u不再是推迟场，而是超前场 。 ),(vRtzyxu （1.1.31）式为用于波场反向外推的Kirchhoff积分式。它可用于上行波的反向外推，也可用于下行波的反向外推。当然，这种外推与正向外推不同，它不代表一个物理过程不代表一个物理过程，而只是一只是一种重建波场的计算过程种重建波场的计算过程。 3地震反射波场成像地震反射波场成像 从波动场的观点叙述反射波成像的一般原理。 地震成像

25、-地震偏移反射系数值-反映该反射点反射系数相对值的反射波振幅反射成像实际上就是把地面上观测到的反射波归位到产生它的反射点上去。地震偏移与地震成像在现阶段可以视为同一概念。 地震偏移成像地震偏移成像: 一是上行波场的反向外推; 二是在外推波场中提取成像值。Claerbout提出下述反射波成像原则提出下述反射波成像原则： 反射面位于这些点上，其入射波的初至与反射波的产生时间相同。 如图1-12所示 反射波成像的基本公式可写为： ),(),(),(ddtzxdtzxuzxMap（1.1.32） （1.1.32）式没有考虑反射系数随着入射角变化的情况，它实质上是相位信息的公式。或者说，它对接近法线入射

26、的情况时基本是正确的，能够反映反射系数在各点上的变化情况。 应用（1.1.32）式涉及到要选择下行波的初始时间。这是一个困难问题。我们通过假设下行波是最小相位假设下行波是最小相位而避开这个问题。我们把 作为初始时间，可推出如下的反射图象公式： dtdttzxdtzxuzxMap),(),(),(（1.1.38） 当下行波是脉冲波时, (1.1.38)式很精确。但是，如果 是 一个短延续长度的子波时, 它只是一个很好的近似成像公式。 ),(tzxd爆炸反射面的概念爆炸反射面的概念: 水平叠加剖面后的自激自收剖面等价于在反射界面上同时爆炸产生地震波，并以半速度向外传播，在地面上观测到的上行波剖面。

27、这就是Loewenthal等人首先提出来的爆炸反射面爆炸反射面的概念（图1-14）。这个概念对于理解水平叠加剖面的偏移成像是很重要的。因为它比较直观地说明了这种剖面的成像原理。比前面所述的反射波成像的一般原理要容易理解得多。在这里我们引入了两种等价。一个是水平叠加剖面等价于自激自收剖面。另一个是自激自收剖面等价于界面上同时激发在地面记录的上行波剖面。这种等价只是概念上的，实际上只有一种水平叠加剖面，并没有三种剖面。 1.2 波动方程偏移波动方程偏移 地震偏移成像技术发展到今天已经产生了各种形式的在各种域实现的方法。历史上曾经起过作用的根据几何光学原理的成像方法已经被淘汰。现在正在流行的是建立在

28、波动方程基础上的三种方法，即Kirchhoff积分法，有限差分法和积分法，有限差分法和F-K法及其各种变形法及其各种变形。这三种方法由于有相同的数理基础，因此它们的原理相同。同时，因计算方法不同，它们之间又有许多不同之处。下面讨论三种方讨论三种方法对水平叠加地震剖面的偏移法对水平叠加地震剖面的偏移。 一频率一频率-波数域波动方程偏移波数域波动方程偏移 采用爆炸反射面的理论。为了成像，要求向地面以下反向外推地震波场。假定z轴垂直向下为正，测线沿x轴，则u(x,z,0)表示偏移后的真实剖面，而u(x,0,t)是未偏移的叠加剖面。 在均匀各向同性完全弹性介质中，用半速度代替地震波传播速度，则标量波动

29、方程变为： 0)(42222222zuxuvtu（1.2.1） ukzuukxuutukkutzxuzxzx),(),(222222222（1.2.2） 对（1.2.1）式进行傅里叶变换并利用（1.2.2）式有 0)(42222zxkkv（1.2.3） 2212zxzkkkv其中正号代表上行波，负号是下行波。 1Stolt偏移法偏移法 设 为 的二维傅里叶变换，对（1.2.1）式进行上述变换得到： 0)(422222ukkvtuzx),(tkkuzx),(tzxu将（1.2.3）式代入上式有 0222utu按上行波求解，即取正值得 tizxzxekkAtkku),(),(其中A与t无关。令t=

30、0，上式变为： ),()0 ,(zxzxkkAkku从而， 是待求的偏移剖面 的傅里叶变换。 ),(zxkkA)0 ,(zxu下面讨论用水平叠加剖面 如何求出 。对 做傅里叶逆变换得： ), 0 ,(txu),(zxkkA),(tkkuzxzzkxkitizxxdkeekkAdktzxuzx)(2),(41),(令z=0，上式变为： zxktizxxdkekkAdktxux)(2),(41), 0 ,(（1.2.4） 设水平叠加剖面 的二维傅里叶变换为 ，则 ), 0 ,(txu),(xkBdtetxudxkBxktixx)(), 0 ,(),(（1.2.5） 其逆变换为： dekBdktxu

31、xktixxx)(2),(41), 0 ,(（1.2.6） 比较（1.2.4）与（1.2.6）有dkBdkkkAxzzx),(),(这样 zxzxdkdkBkkA),(),(按上行波 取正号并对 微分得 zk2222/12)/12,(),(zxzxzxzxkkvkkkvkBkkA（1.2.7） 对 做二维傅里叶逆变换得到： ),(zxkkAzxzkxkizxdkdkekkAzxuzx)(2),(41)0 ,( （1.2.8） 就是要求取的偏移剖面。 )0 ,(zxu输入零偏移距剖面 二维付氏变换, 公式(1.2.5) 得到用公式（1.2.7）把 映射成 并 标定振幅，得到二维付氏逆变换, 公式

32、(1.2.8) 得到偏移剖面图1-16 均速Stolt偏移流程图 ), 0,(tzxu）,(xkBzk）zxkkA,(）0,(tzxu 上述偏移原理见图1-15。由图1-15和（1.2.7）式可看出，在每个频率 移向新的频率 时，要乘上一个振幅比例 。通过这个频率移动，把视倾角 转换为真倾角 。 其流程见图1-16。 上述频率波数域的偏移方法称为Stolt偏移方法偏移方法。 Stolt法的偏移效果见图1-17，1-18和1-19。 2/vk2/zvkkkz/2Gazdag相移法相移法 对标量波动方程（1.2.1）相对x和t做二维傅里叶变换得到： 0)(2222ukkzux（1.2.9） 式中

33、。求解（1.2.9）式得出F-k域的向下外推公式 vk/2), 0 ,()exp(),(22xxxkukkizzku（1.2.10） 偏移成像公式是把上式变换回到空间-时间域，并取t=0时刻的波场值为成像值。即 dezkudktzxuxikxxx),(41)0,(2 Gazdag相移法的流程见图1-20。 偏移效果见图1-21。 输入零偏移距剖面), 0,(tzxu 二维付氏变换, ), 0,(zkux在新的z处计算 （kx , z , ）uu用相移算子exp（-ikzz）将地面数据向地下延拓相加所有频率（成像原理，t=0）， （kx , z , t=0）u在x方向作付氏逆变换,）0,(tzx

34、u偏移剖面图1-20 Gazdag的相移偏移法流程图 二克希霍夫积分法波动方程偏移二克希霍夫积分法波动方程偏移 前面导出了波动方程边值问题的Kirchhoff积分解。下面研究把它用于地震成像问题。 现在转写反向外推的Kirchhoff积分如下： SdSnuRRnutunRvRtzyxu11141),(（1.2.11） 图图1-22 求地震问题Kirchhoff积分解图示 因为在（1.2.11）式中需知道 ，即波场在地面法向的导数值。但是，这个导数值目前是无法观测和计算的，因此，需想法去掉含 的项。为此不再用（1.2.11）式，而从格林定理一般式： 0zununuVSdSnuwnwudVuwwu

35、)(22（1.2.12） 出发，设格林函数w为： 11rrw（1.2.13） 来代替1/R，则可以达到目的。上式中： 22020) ()()(zzyyxxr22020) ()()(zzyyxxr（1.2.14a） （1.2.14b） 把（1.2.13）式代入（1.2.12）式，得 011111122udVrudVrdSnurrrrnuVVS（1.2.15） 由此求出向下外推的Kirchhoff积分为： dAtuzrrzrrvzurrrrzutzyxuA0000111111141),(（1.2.16） 式中 A为地面的面积。 求下面的微分：33011rzzrzzrrz220011rzzrzzzr

36、rzrr（1.2.17a） （1.2.17b） 把以上二式代入（1.2.16）式，得 dAturzzrzzvzurrrzzrzzutzyxuA2203311141),(（1.2.18） 当我们把z取成地面上的点时，即z=0时，则有： 011,rrRrr（1.2.18）变为下列形式： dAtuuRvRvtzyxuAcos21),(（1.2.19） 式中 22020)()(coszyyxxzRzA为地面的面积。 （1.2.19）式又可写为： dxdytvRtyxuvRtyxuRvRvtzyxuA), 0 ,(), 0 ,(cos21),(0000（1.2.20） （1.2.20）式与下式等价： R

37、vRttztyxdAudttzyxuA)(), 0 ,(21),(000000（1.2.21） 式中 vRtt0 现在我们来证明（1.2.21）与（1.2.20）式等价。我们先对 求导，再对 求积分。该过程如下： 0z0tAAAAdxdytvRtyxuvRtyxuRvRvvRtyxuRvtvRtyxuRvdAvRttRvvRttRvtyxdAudtRvRttztyxdAudttzyxu), 0 ,(), 0 ,(cos21), 0 ,(), 0 ,(cos21)()( cos), 0 ,(21)(), 0 ,(21),(00000000000000000000（1.2.22） （1.2.22）

38、式与（1.2.20）式完全相同，因此（1.2.21）式也与（1.2.20）完全相等。由（1.2.14a）式可知 RvRttzRvRttz000（1.2.23） 将（1.2.23）式代入（1.2.21）式，得到 RvRtyxudAzRvRtttyxudtdAzRvRtttyxudAdtztzyxuAAA), 0 ,(21)(), 0 ,(21)(), 0 ,(21),(000000000000（1.2.24） 根据褶积的定义，我们把（1.2.21）式写成三维褶积符号形式，则有 000)(21*), 0 ,(),(RvRtztyxutzyxu（1.2.25） 其中 2220zyxR把（1.2.25

39、）式对x,y和t进行傅里叶变换，则可写为： ),(), 0 ,(),(zkkHkkuzkkuyxyxyx（1.2.26） 式中褶积算子H为： )(00)(21),(ykxktiyxyxeRvRtdxdydtzzkkH（1.2.27） 对t积分得到： yikxikzyxviyxeezyxedxdyzH22222221（1.2.28） 括号内积分是第一类Hankel函数，上式可写为 yikxyezykvHidyzH2222)2(0（1.2.29） 利用圆柱函数间的关系： 00iNJH(1)0式中 和 是Bessel函数和Neuman函数，把它们代入对y的积分式（1.2.29）中，则最终得到H的表达

40、式为： 0J0N222222222222222222222cossinyxyxkkvizyxkkvizyxyxyxyxekkvezikkvkkvzikkvkkvzzH （1.2.30） 把（1.2.30）式代入（1.2.26）式中，得 222exp), 0 ,(),(yxyxyxkkvizkkuzkku（1.2.31） 上式与频率-波数域的向下外推公式一致。因此，在常速介质中Kirchhoff积分法与频率-波数域的波场向下外推公式完全等价。 下面讨论利用Kirchhoff积分法对水平叠加剖面进行波动方程偏移的步骤。 将水平叠加剖面看做是炮检距为零的自激自收地震剖面u(x,y,0,t)； 利用爆

41、炸反射面的思想将自激自收剖面等效为在反射界面上同时激发产生地震波,以半速度向外传播，在地面上观测到的上行波剖面u(x,y,0,t)； 利用（1.2.24）式将单程的上行波剖面u(x,y,0,t)向下延拓，得到深度为z的面上的波场值 RvRtyxudAztzyxuA), 0 ,(21),(00（1.2.32） 根据成像原理，对所有地下点（z0）取t=0时的波场值，即可实现三维偏移成像。此时，成像值为 RvRyxudydxzzyxu), 0 ,(21)0 ,(0000（1.2.33） Kirchhoff积分法的偏移效果见图1-23。 三有限差分法波动方程偏移三有限差分法波动方程偏移 下面讨论使用有

42、限差分法对水平叠加地震剖面的偏移问题。为了把上行波方程表示为空间-时间域的表达式，需要把上行波方程表示为某种近似式。然后在空间-时间域研究其差分方程及求解问题。最后讨论一些计算方法和效果。 1上行波的空间上行波的空间-时时间域方程间域方程 为了适应介质速度的空间变化，我们要在空间空间-时间时间域域中进行偏移成像或地震图的模拟工作。首先就要把上行波方程表示在空间-时间域中，这需要用到某种根式展开。 1）二项式展开）二项式展开 下面我们将用到 这样的二项式展开，在这里我们介绍几种展开式。 21)1 (X（1）Taylor展开展开 这是一个众所周知的显式展开式，它一般表达为： 54322567128

43、516181211)1 (21XXXXXX（1.2.34） 展开条件 。 1X 如果把这种展开式用于微分算子，在不进行辅助处理时将找不到稳定的有限差分方程来解相应的微分方程。因此我们在使用二级近似以上的展开式时不能用这种展开式。 （2）连续分式展开，或称为）连续分式展开，或称为Pade展开展开 这个展开式表示为如下形式（ ）： 1X41414121)1 (21XXXXX（1.2.35） 这是一种隐式展开式。其各级展开式如下。 一级展开式： 21)1 (21XX（1.2.36a） 二级展开式： 4121)1 (21XXX（1.2.36b） 三级展开式： 414121)1 (21XXXX（1.2.

44、36c） 高级展开式可依此类推。 （3）迭代展开）迭代展开 这种隐式展开法，是把前一级的展开结果代入下一级的展开式中。设 21)1 (XR则逐次迭代展开式可表示为（ ）： 1X1,1101RRXRnn（1.2.37） 用这种展开方法得到的各级展开式如下。 一级展开式： 211XR（1.2.38a） 二级展开式： 41212212XXXXR（1.2.38b） 三级展开式： 41412122213XXXXXXR（1.2.38c） 高级近似式可依此类推出来。 从（1.2.36）和（1.2.38）公式组可以看出，后两种展开是等价的。 2）上行波的空间）上行波的空间-时间域方程时间域方程 在第一节已经求

45、出了频率-波数域的上行波方程(1.1.12)式： ukvizux222用迭代展开法展开上行波方程： uRvkviuvkvizunxx11122222221（1.2.39） 式中 1,1/101222RRvkRnxn由（1.2.39）式求出各级近似式如下。 一级近似式： uvkvizux21222（1.2.40a） 二级近似式： uvkvkvizuxx421222222（1.2.40b） 三级近似式： uvkvkvkvizuxxx21181211222444222（1.2.40c） 高级近似式可依次类推。现在，我们把（1.2.40a）式转换到空间-时间域，求出一级近似方程。 上行波方程（1.2.

46、40a）式可改写为： 0222uvkuvzuix（1.2.41） 对（1.2.41）式进行傅里叶反变换： 0),(4),(21),(21)(2)(2)( xxktixxxxktixxxktixdkdezkukvdkdezkuvdkdezkuizxxx（1.2.42） 根据傅里叶变换的微分性质： xxktixdkdezkuittzxux)(),(21),( xxktixdkdezkuttzxux)(222),(21),( xxktixxdkdezkukxtzxux)(222),(21),(（1.2.43a） （1.2.43b）（1.2.43c）把(1.2.43a)、(1.2.43b)和(1.2.

47、43c)式代入(1.2.42)式，得到： 01222222tuvxuvtzu（1.2.44） 上式就是空间-时间域的一级近似的上行波方程，常常被称为 方程方程。 15 同理可求出空间-时间域的二级及二级以上近似的上行波方程。经推导，空间-时间域的二级近似的上行波方程为： 04314233323223txuvtuvzxuvztu（1.2.45） 上式常常被称为 波动方程波动方程。 453）浮动坐标系中的单程波方程）浮动坐标系中的单程波方程上行波方程在一定的浮动坐标系中可以简化。 我们对二维波动方程： 22222221tuvzuxu（1.2.46） 进行如下的坐标转换： vzttzzxx坐标变换前

48、后波场本身是不变的，因此存在： ) , , (),(tzxutzxu（1.2.47） （1.2.48） 从（1.2.47）和（1.2.48）导出下列导数等式： 2222 xuxu（1.2.49a） 2222222212tuvtzuvzuzu（1.2.49b） 2222 tutu（1.2.49c） 将（1.2.49）各式代入（1.2.46）式中得到新坐标系中的波动方程为： 0222222zutzuvxu（1.2.50） 上式变换到频率-波数域为： 02222zuzuviukx（1.2.51） （1.2.51）式可改写为： ukvuvizx2222（1.2.52） 从上式得到下列关系式： zviz

49、它表示坐标变换前后的算子关系坐标变换前后的算子关系。因此，上行波方程可表示为： ukviuvizx21222（1.2.53） （1.2.53）式的右端项可展开为各级近似式，便得到上行波各级近似方程。 一级近似式 ukvizux22（1.2.54） 二级近似式 ukvvkizuxx4122222（1.2.55） ukvkvvkizuxxx21422423432三级近似式 （1.2.56） 高级近似式可依次类推。由以上各式用前述方法可求出空间-时间域的各级近似方程。下面给出一级和二级近似方程。 一级近似方程一级近似方程 用推导（1.2.44）式那样的方法可以求出一级近似方程为： 02222xuvt

50、zu（1.2.57） 与（1.2.44）式相比，减少了一项。从而也减少了计算时间和差减少了计算时间和差分时的时间层（少了一层）分时的时间层（少了一层）。另外，保持了计算的稳定性保持了计算的稳定性。 二级近似方程二级近似方程用推导（1.2.45）式那样的方法可以求出二级近似方程为： 0242323223txuvzxuvztu（1.2.58） 这个方程与（1.2.45）式相比也是少了一项。这也会减少计算工减少计算工作量作量。 2有限差分法地震偏移技术有限差分法地震偏移技术 如前所述，水平叠加地震剖面水平叠加地震剖面可以看做是自激自收地震剖面自激自收地震剖面；又可以看做是所有反射面同时爆炸产生波源向

51、地面传播，被地面的接收器记录的上行波剖面上行波剖面。对于这种观测结果，为了成像，要求向地面以下反向外推地震波场向地面以下反向外推地震波场。在外推过程中假设地震剖面上无任何多次波，也不存在任何规则干扰波，如折射波等。如果在剖面上存在这些波，在外推过程中也都按反射一次波处理，但它们是不能正确归位的，只能造成偏移成像剖面的干扰。因此，如果存在这些波，应当在偏移处理前把它们滤掉。 1）浮动坐标下的有限差分法地震偏移）浮动坐标下的有限差分法地震偏移 采用浮动坐标系, 只讨论一级近似的上行波二阶偏微分方程（1.2.57）的有限差分偏移问题。考虑到爆炸反射面的概念，用v/2速度代替v。这样，（1.2.57）

52、式可重新写成： 04222xuvztu（1.2.59） 这里的速度，假设它是常数，在实用中它可以随 和 而变。 (1.2.59)式存在于空间 内，给定的混合条件为： )0 ,0 ,(ZzTtx0),(Tzxu0),(tTzxu),(), 0,(txtzxu0),(tZxu（1.2.60a） （1.2.60b）（1.2.60c）（1.2.60d）式中 即为在地面所观测的地震波场。 ),(), 0,(txtzxuxz 目的目的: 通过解上述微分方程求出地面以下任何点（ ）上的曾经在该点出现过的上行波的波场值(位移或压力场振幅) 。这样的问题是适定的，可解的。 0,zx),(tzxu 下面用有限差分

53、法来解这个方程。采用对称隐式(Crank-Nicolson)的差分格式（图1-24）。 实用中，常采用下列符号： 图1-24 差分格式 )2(1)(1),(),(21, 121,21, 1221,1,21,njinjinjinjixxnjinjinjitnjiuuuxuuututzxutjznxiuu 对图1-24格式的中心点进行差分，（1.2.59）式可化为如下的差分方程： )()(21212121,1,1,njinjixxnjinjituuauu（1.2.61） 式中 2222,8xxxxxttxxtxztva由（1.2.61）式整理得： njixxtnjixxtuaua2121,1,)(

54、)(（1.2.62） 为了求出以整采样点为依据的表达式，（1.2.62）式可写成： )(2)(2,1,1,1,11,1,11,njinjixxnjinjinjinjixxnjinjiuuauuuuauu（1.2.63） 从（1.2.63）式可以解出求 的递推式。把求 做为递推结果值，不取 作为递推结果值是从物理条件（1.2.60）式考虑的，这也是计算稳定性所要求的。（1.2.63）式经过简单整理后得： ),(1,tzzxuunji1,njiu),(11,ttzzxuunjinjixxnjixxnjixxnjixxuIuIuIuI1,11,1,)()()()(（1.2.64） 式中 )0 , 1

55、 , 0(,1622Ixztva算子 （1.2.64）式写成矩阵式为： njnjnjnjAUUUBAU1111)(（1.2.65） 式中 、 为三对角矩阵。 为在 深度层和 时间层上沿 轴的波场值的列向量。 ， 和 为相应层上沿x轴的波场值列向量。 ABnjUznztjtxnjU11njU11njU这些矩阵和列向量表示如下： 2102102121A2102102121B（1.2.66） （1.2.67） TnjInjinjnjnjnjuuuuuU1,1,1, 31, 21, 11.,TnjInjinjnjnjnjuuuuuU11,11,11, 311, 211, 111.,TnjInjinjn

56、jnjnjuuuuuU,1, 3, 2, 1.,TnjInjinjnjnjnjuuuuuU1,1,1, 31, 21, 11.,（1.2.68） （1.2.69）（1.2.70）（1.2.71）如果把差分算子 取为 ，这相当于差分加权，可提高差分精度。这时差分方程（1.2.64）可改写为： xx)/(xxxxInjixxnjinjixxnjixxuIuuIuI1,11,1,)()()()(（1.2.72） 式中 。 61上式可以写成矩阵式(1.2.65) ,不过，矩阵A和B由下面的矩阵表示。 )(21)()(00)()(21)()()(21A)(2100)(21)(21B（1.2.73） （1

57、.2.74） 只要矩阵A可逆，（1.2.65）可解。得到： njnjnjnjUUUBAU11111)(（1.2.75） 矩阵A中的主对角元素是优势元素。因为 永远大于零。因此总有： 221或者 612232因而矩阵A的行列式不为零，即 0detA所以矩阵A是可逆的。 （1.2.76） 矩阵方程（1.2.65）式可用矩阵方法矩阵方法求解，也可以改写为代数方程组后用追赶法追赶法求解，要求：差分方程必须稳定要求：差分方程必须稳定。 在浮动坐标系中，成像时间是 。 vzt/22）一般坐标下的有限差分法地震偏移）一般坐标下的有限差分法地震偏移 若按一般坐标系，叠加地震剖面的有限差分偏移成像过程可用图1-

58、25说明。图1-25为（x,z,t）坐标系，地震剖面在（x,t）平面上，偏移剖面则在（x,z）平面上。有限差分偏移是按一定步长的z来外推地震剖面（x,t）,每外推一个步长，就将t=0的波场作为输出。这些输出结果就组成了偏移剖面（x,z,t=0）。 3）有限差分法逆时偏移）有限差分法逆时偏移 Baysal等人在1983年提出了有限差分逆时偏移的方法，它是从一个波场为零的（x,z）起始平面，按时间反推，并以地震剖面资料u（x,z=0,t）作为每一步进时间的边界条件（z=0）,得出时间t=0的（x,z）平面就成为偏移结果 u（x,z,t=0）（见图1-26）。 图图1-25. 地表z=0的地震剖面（

59、x,t），它向下延拓获得各个离散深度上的时间剖面。这里用粗黑箭头指出外推方向，偏移获得的剖面用t=0（根据成像原理）的（x,z）平面表示。 图图1-26. 逆时偏移从数据体底部的全零（x,z）平面开始，按时间向t=0反推，计算出不同时间的（x,z）平面切片；这些地下切片在图中用一系列水平面来表示，反推方向按粗黑箭头所示，每个时间平面（x,z）都包含有出自地震剖面的边界值（虚线表示的z=0平面上的x线 ），t=0的（x,z）平面即为偏移剖面（顶部的水平面） 3吸收边界条件吸收边界条件仍然讨论二维情况。二维波动方程： 012222222tuvzuxu（1.2.77） 一般求解域：0,0 ,),(t

60、zxtzx实际求解域：Ttbzaxatzx0 ,0 ,),(0),(tzxua0),(tzxxua0),(tbxu0),(tbxzu在求解的过程中，一般给定如下的边界条件： （1.2.78a） （1.2.78b） （1.2.78c） （1.2.78d） 这样，就在两边和底界人为地造成了边界，称之为计算边界计算边界。这样的边界不可避免地会产生边界反射边界反射。这个边界反射是很强的。在现代偏移技术中为了避免在有用的地震剖面范围内出现边界反射波，常常要在地震剖面的两边进行扩边，扩边道有时要达到96道或更多。这就会使计算量和设备资源使用量增大。为了克服这个不足，提出了使用吸收边界条件吸收边界条件的算法