第1章地震偏移成像基础

第一章 地震偏移成像基础地震偏移技术是现代地震勘探数据处理的三大基本技术之一。它是在过去的古典技术上发展起来的，其它两大技术都是从其它相关学科引进到地震中来的。所以，偏移技术具有地震勘探本身的特征。但是，地震偏移方法本身由于使用计算机而引起了许多革命性的变化。这就是把它从研究简单的探测目标的几何图形进而发展成研究反射界面空间的波场特征、振幅变化和反射率等。本章主要介绍地震偏移成像技术的基础知识。首先给出偏移成像的概念；第二节介绍有限差分法的基础知识；第三节叙述基于波动方程的波场外推与地震成像原理；第四节讨论波场外推的Kirchhoff积分法；第五节简单分析Born近似和Rytov近似；最后阐述基于DeWolf近似、薄板近似、屏近似和相屏传播算子计算反向散射波场的方法。§1.1 偏移成像的概念反射地震方法是根据在地面上以一定方式进行弹性波激发，并在地面的一定范围（孔径）内记录来自地下弹性分界面的反射波来研究地下地质岩层结构及其物性特征的一种方法。因此，也可以把它看做是一种反散射问题。就反射地震观测方式的特点，它的成像问题要分做两步，第一步是按照一定的方式记录到达地面的反射波，第二步用计算机按一定的计算方法对观测数据进行处理，使之成为反映地下地质分层面位置及反射系数值的反射界面的像。而地震偏移技术就是在第二步过程使反射界面最佳地成像的一种技术。地震偏移可在叠前做也可在叠后做。叠前偏移是把共炮点道集记录或共偏移距道集记录中的反射波归位到产生它们的反射界面上并使绕射波收敛到产生它的绕射点上。在把反射波回投到反射界面上和绕射波收敛到绕射点上时要去掉传播过程的效应， 如扩散与衰减等。最后得到能够反映界面反射系数特点的并正确归位了的地震波形剖面，即偏移剖面。叠后偏移是在水平叠加剖面的基础上进行的，针对水平叠加剖面上存在的倾斜反射层不能正确地归位和绕射波不能完全收敛的问题， 采用了爆炸反射面的概念来实现倾斜反射层的正确归位和绕射波的完全收敛。地震偏移的效果见图 1-1和图1-2。地震偏移的类型见表 1-1。地震偏移技术在二十世纪六十年代以前是用手工操作的一种制图技术， 只是用来求得反射点的空间位置，而不考虑反射波的特点。它是一种古典的偏移方法。早期的计算机偏移方法是在古典的偏移方法的基础上提出来的。其中有的成功了，有的失败了。成功的是那些符合波的传播特征的方法。尽管这些方法使用了波前、绕射等地震波传播的惠更斯原理，但只是定性的、概念性的。偏移剖面的质量虽然能够满足最基本的要求，但归位的精度和成像时的波形特征都不是很准确的。因此，研究更有效的地震偏移方法是很迫切的。二十世纪七十年代初J.Claerbout教授首先提出了用有限差分法解单程波动方程的近似式，用地面观测的地震数据重建地震波在地下传播过程中的波场，从这些传播过程的波场中提取使地震界面成像的那些数据，组成地震偏移剖面。由于这种偏移方法在计算过程中要解波动方程或其近似式，所以被称为波动方程法偏移技术。以后，French和Schneider等在绕射偏移法的基础上使用了波动方程解的Kirchhoff积分公式，发展为地震偏移的波动方程积分法。使绕射偏移建立在可靠的波的基本原理上。因而改善了偏移剖面，取得了良好的效果。图1-1(a)共中心点叠加剖面， (b)偏移剖面，(c)明显的绕射波 D、偏移前的倾斜同相轴 B和偏移后的倾斜同相轴A。偏移使倾斜同相轴 B归位到它的真实地下界面 A，并使绕射波 D收敛到其顶点 P。点画线指出了盐丘的边界。图1-2偏移前(a)及偏移后(b)弯曲反射界面 (向斜和背斜 )的形状。详细情况见正文 (模型据 Union OilCompany)在二十世纪七十年代后期， Stolt和Gazdag等又先后提出了在频率 -波数域解波动方程，外推地震波场的方法。这种方法被称为F-K域偏移方法。由于该方法计算简单，效率高，因而很快得到了推广。上述三种波动方程偏移方法是同时并存的，因为它们各有自己的特点，因而不能用一个方法来取代其它方法。使用时视具体条件和要求决定采用何种方法。波动方程偏移方法在最近20年间迅速发展并不断完善，许多人对此做出了有益的贡献。其中，Loewenthal等人的爆炸反射面的概念对于理解叠加剖面的偏移成像具有很大价值，Hubral，Larner等人提出的深度偏移的概念具有很大意义，Berkhout提出的偏移过程是一个空间褶积的概念对于偏移的横向分辨力的理解很有益处，马在田院士提出的高阶方程的分裂算法对提高有限差分法偏移的精度有很大贡献，Yilmaz等提出的双平方根法为解决叠前偏移奠定了基础。现在仍有许多学者还在探索波动方程偏移技术，以期更加完善该方法。表1-1偏移方法分类类型论述叠加解释人员总是想要的剖面。法向射线深度转换严格适用于没有构造倾角且速度只随深度变化的情况。时间偏移适用于叠加剖面上有绕射波或构造倾角以及速度有垂向变化的情况；速度的横向变化不大时也能用。深度偏移用于叠加剖面上有构造倾角和强横向变速的情况。叠前部分偏移叠后偏移适用于叠加剖面与零炮检距剖面等价的情况，但不适合具（PSPM）有不同叠加速度的地层倾角不一致或强横向变速的地区，叠前部分偏移(倾角时差(DMO)校正)能够为叠后偏移提供更好的叠加剖面，但叠前部分偏移只解决具有不同叠加速度的地层倾角不一致的问题。叠前全时间偏移输出偏移剖面，不产生未经偏移的中间叠加剖面，所以不太受欢迎，因为解释人员普遍喜欢既有叠加剖面又有偏移剖面的解释方式。但不论如何这是解决倾角不一致地层问题的最精确方法。叠前部分偏移是这种处理方法的一种简化。叠前深度偏移用于严重横向变速的情况，这时已无法作合适的叠加处理。三维叠后时间偏移用于叠加剖面上出现来自射线平面以外的倾斜同相轴(即垂直测线方向)的情况，这是叠后最常用的一种三维偏移方法。三维叠后深度偏移用来解决与三维地下复杂构造有关的强横向变速问题。三维叠前时间偏移用于叠前部分偏移不适用且叠加剖面上有横向倾斜层反射的情况。三维叠前深度偏移只要计算机机时允许，并且又能精确知道三维速度模型，这是人人乐于接受的处理方法。§1.2 有限差分法的基础知识在计算机上进行数值运算，使用的是离散的和有限的数值，而不是连续的和无限的函数，为此要为离散数值的计算建立基本方法。最基本和广泛使用的方法就是有限差分法，借助这一方法可以研究连续物理系统的性质，近似地、但相当精确地解出各种数理方程问题。本节将概要地介绍有限差分法的基础知识，供以后各章解偏移方程使用。一．差分方程的建立与求解1．有限差分的概念我们感兴趣的是用有限差分法解各类微分方程，因此，要把导数用有限差分来近似，所以我们这里只研究用有限差分近似导数的方法。（1）一元函数差分法当一个函数u和它的各阶导数是变量x的单值的、有限的和连续的函数时，可以用泰勒原理展开为：u(xx)u(x)xdu1x2d2u1x3d3u（1-1）dx2dx26dx3和u(xx)u(x)xdu1x2d2u1x3d3u（1-2）dx2dx26dx3以上二式相加，得：u(xx)u(xx)2u(x)x2d2u(x4)（1-3）dx2式中O(x4)表示包含有x的四阶和高于四阶以上的项。从（1-1）~（1-3）式我们可以导出（图1-3）下列表示式：一阶向前差分：xu1[u(xx)u(x)]xdu(x)（1-4）dx一阶向后差分：xu1[u(x)u(xx)]xdu(x)（1-5）dx一阶中心差分：xu1[u(xx)u(xx)]2xdu(x2)dx二阶差分：21xxuxux2[u(xx)d2u2（1-7）dx2(x)（2）多元函数差分法图1-3一元函数差分法

（1-6）2u(x) u(x x)]图1-4二元函数差分网格设有一个二元函数u(x,t)，我们用网格把它们离散（图1-4），令x=ix，t=jti和j为整数。用ui,j=u(ix,jt)表示各网格点上的函数值。现在我们用泰勒级数展开下面各点之值为：ui1,ju[(i1)x,jt]u(xix,tj)ui,jx(u)i,j1x2(2u2)i,j1x3(3u3)i,j（1-8）x2x6xui1,ju[(i1)x,jt]u(xix,tj)ui,ju122u13(3u（1-9）x()i,jx(2)i,jxx3)i,jx2x6由此可以求出对 x的一阶和二阶差分为：xu1(ui1,jui,j)xu（1-10）(x)x1xu(ui,jui1,j)xu（1-11）(x)xxu1ui1,j)(ui1,j2 x(x2)（1-12）xxxux2u1(ui1,j2ui,jui1.j)2x2u(x2)（1-13）x2同理可求出对t的一阶和二阶差分：tu1(ui,j1ui,j)tu(t)（1-14）ttu1(ui,jui,j1)tu(t)（1-15）ttu1(ui,j1ui,j1)2tu(t2)（1-16）tttut2u12(ui,j12ui,jui,j1)t2u(t2)（1-17）t22．建立差分方程建立差分方程的方法很多，有积分法，物理量守恒法，变分法和最小平方法等。对于常系数的微分方程来说， 积分法是最简便的。用积分法构造差分方程的过程见参考文献 [1]。3．差分方程的格式差分方程的格式基本可分为两大类：即显式格式和隐式格式。在实际工作中又可以衍生出许多格式，甚至可以用显式与隐式联合形式的差分格式。在这里我们仅以抛物型偏微分方程为例说明常用的显式格式和 Crank-Nicolson格式。1）显式格式设有抛物型偏微分方程：u2u0（1-18）tx2经推导可具体写出差分方程为：ui,j1ui,jx2t(ui1,j2ui,jui1,j)（1-19）t所用的差分网格如图1-5所示。令，则（1-19）式可写为：x2ui,j1ui,j(ui1,j2ui,jui1,j)（1-20）其中ui,j，ui1,j和ui1,j为已知值。从第j时间层上的已知值， 可用（1-20）式直接计算出第 j 1时间层的值。以此类推，可解出全时间上的物理量值。图1-5显式差分格式 图1-6隐式差分格式2）Crank-Nicolson隐式格式虽然显式法计算上很简单，但它有一个严重的缺点，即时间步长 t一定要很小，必须满足0t1，才能保持计算上稳定并达到必要的精度。Crank-Nicolson在1947年x22提出了一个使所有的有限r值都满足计算要求（收敛性和稳定性）的方法。他们把2u/x2项用第j和第j+1时间行上的平均差分来逼近。求出下列的差分方程（图1-6）：ui,j1ui,jui1,j12ui,j1ui1,j1ui1,j2ui,jui1,j（1-21）t2x2x2由此得出下列等式：ui1,j1 2(1 )ui,j1 ui1,j1 ui1,j 2(1 )ui,j ui1,j （1-22）其中，t2。（1-22）式的左边包含有3个未知数值，右边3个值是已知的。x如果在每个时间层上有I个格点，则方程（1-22）根据零时间层(j=0)计算第一时间层(j=1)的网格点上的函数值ui1时可列出I个联立的方程组。由于j=1时间层上的初值和边值是给出的，因此方程组表现为三对角线阵。对于这样的方程组，可用代数的追赶法求解，也可以用矩阵方程来分析。4．差分方程的解法1）追赶法1-22）方程组可用（1-23）式表示：iui 1 iui iui 1 i其中，0 i I。补充的边界条件为：0u1 0u0 0IuI 1 IuI I

1-23）1-24）我们要求出未知值ui。用递推法求解，即已知i点上的ui值可推出i+1点上的ui1值。这就要先求出中间助算未知数xi和yi，其关系式为：ui1xiuiyi（1-25）为此要求出xi,yi。（1-25）式对任何i点都是成立的。我们把它代入（1-23）式，得：i(xiui yi) iui iui 1 i或者写为：uiiui1iiyi（1-26）ixiixiii这个公式与（1-25）式相同，它把i点和i1点用（1-27）式：uixi1ui1yi1（1-27）连结起来。由此求出下列等式：xi1i，yi1iiyi（1-28）xixiiiii用这个关系式可以求出所有的xi和yi的值。这样，（1-25）和（1-28）就给出解三对角线方程组的两个递推过程。利用（1-28）式，我们在网格上从iI到i0求出所有的xi和yi值。在求出所有的xi和yi之后，从i0向iI用（1-25）式求出所有的ui1值。在iI点确定初始值xI1和yI1。在i0点从边界条件确定u01值。把第I点的边界条件（1-24）与递推关系式（1-25）比较，则得：xI1IIyI1I（1-29）I在用递推关系求出所有的xi和yi，一直到x0和y0之后，则用（1-24）和（1-25）式可求出未知函数值u01为：u01x000y0（1-30）0x00从u01出发用（1-25）式可求出所有的ui1值。根据上述递推原理可求出所有时间层上的ui,j值。（2）矩阵法现在从矩阵观点来分析三对角线方程组的求解问题。令w为已知量，需要求出的向量为u,它们满足矩阵方程：Auw（1-31）其中A是三对角矩阵。我们引入平移的矩阵算子L和L，它们对向量u的分量有下列关系：L{ui}{ui1}（1-32）L{ui}{ui1}（1-33）由于A是三对角矩阵，故（1-31）可写为：(PLQRL)uw（1-34）其中P，Q和R是对角矩阵。（1-34）式给出双向递推过程，可分为两步。首先有下列关系式：LuXuy（1-35）求对角矩阵X和向量y。把（1-35）式与对角阵P相乘，写为下列形式：(PLPXO)uPy（1-36）式中，O为零阵。从（1-34）式减去（1-36）式，得(QPX)uRLuwPy（1-37）或者写为：u(QPX)1RLu(QPX)1(wPy)（1-38）1-38）式与（1-35）式在形式上是一致的。为了使它与所讨论的问题相联系，必须使它们等价。因此得出下列等式：LX(QPX)1RLy(QPX)1(wPy)（1-39）方程（1-35）和（1-39）组成两步递推过程。这些结果与标量方法所得结果一致。因为P、Q、R和X是对角矩阵，所以它们每列的对角线上的元素相应地是 i、i、i和xi。显然，矩阵结果与代数方程是一致的。二．收敛性、相容性和稳定性为了使差分方程的解能够相当精确地逼近于相应的微分方程的解，我们必须讨论它所需要的条件。这是不同的，但又互相关联的两个条件。第一个是近似的差分方程的精确解向微分方程的解收敛的条件， 第二个是在解差分方程过程中任何误差不可无界地增长的条件。1．收敛性令u是变量x和t的偏微分方程的精确解，u?是用来近似偏微分方程的差分方程的精确解。当t和xuu，就说有限差分方程是收敛的。趋近于零时，?趋近于某个固定点上的e u u?之差称为离散误差。讨论收敛性问题的简单例子请看参考文献[1]。一般来说收敛性问题常常是难以研究的。因为离散误差的最终表达式经常含有未知导数项，无法估计其边界值。但是，对于线性偏微分方程中的抛物型与双曲型方程来说可以用稳定性和相容性来研究收敛性。对此存在Lax等价原理。Lax等价原理该定理说，当给出了一个适定的初值问题和它的有限差分方程时，如果该差分方程满足相容性条件，则稳定性是收敛性的充分而且必要的条件。关于这个定理的证明请看参考文献 [8]。下面我们给出相容性和稳定性的定义， 研究相容性和稳定性问题的实例和方法请看参考文献 [1]。2．相容性设F(?)0代表在(i,j)点上的差分方程。如果把?用u来代替，则()就称为i,jFi,juu在网格点(i,j)上的局部截断误差。如果在网格步长趋于零时这个误差亦趋于零，则称差分方程与偏微分方程相容。3．稳定性在解差分方程过程中，除了离散误差之外，还有舍入误差。这就是差分方程的理论精确解ui,j与计算机输出的实际值Ni,j之差。全部误差有离散误差和总舍入误差组成。??Ni,j)（1-40）i,jui,jNi,j(ui,jui,j)(ui,j离散误差在一定的差分格式和差分网格步长下是可以控制的。而总的舍入误差不仅与差分方程有关，而且与计算机舍入误差有关，也与算法有关。差分解法是以步进方式工作的。在逐步推进过程中这种误差也会逐渐积累。 这种误差积累是保持有界还是恶性增长以致把解淹没在积累的误差中，这就是数值解的稳定性问题。数值稳定性是差分格式的必须条件。也是差分方程收敛性的必要条件。因此，在确定一个差分格式作为计算手段之前，一定要研究该差分格式的稳定性。常用的有两种研究稳定性的方法。一个是 Fourier分析法，这是一种比较简单的，但不考虑边界条件的方法。另一种方法是矩阵法。它把方程组表示为矩阵，并考查与矩阵有关的本征值。此外还有能量法。§1.3 基于波动方程的波场外推与地震成像原理使用波动方程进行偏移，首先就是要重建反射波的原来波场。反射界面上刚刚产生的反射波，就认为是该反射面的像。为了使用波动方程公式进行波场外推，一般是把波动方程分解为上行波方程和下行波方程。这是由于波动方程一般式既满足上行波也满足下行波。我们在成像中使用的是上行波方程。在这一节里我们将讨论上行波与下行波方程，波场外推方法和地震成像的原理。一．上行波和下行波波动方程有两个解，一般表示为 {exp[ i (t r/v)]}/r，其中 r是从源到观测点的距离，是正实数。在地震勘探中一般取深度方向向下为正 z的方向。向正 z方向传播的地震波称为下行波，即用{exp[i(tr/v)]}/r所表示的波。向负z方向传播的波为上行波，即用{exp[i(tr/v)]}/r代表的波。下行波即入射波，上行波为反射波。为了利用波动方程外推波场要求把它分解为上行波方程和下行波方程，然后才能利用它们进行波场外推。否则将难以进行外推的计算。从Claerbout导出的结果可以看出，既使只存在z方向介质不均匀的情况下，上行波和下行波也是相互耦合的。只有在均匀各向同性完全弹性介质的情况下上行波和下行波才是分离的。由于现代地震研究的主要是纵波，因此，我们可以用声波方程来代替弹性波方程。设u和w表示x和z方向上的质点位移分量，v表示纵波的传播速度。考虑到在地震勘探中我们研究的是无源问题，因此在均匀各向同性完全弹性介质中有如下的二维波动方程：2u2u12u（1-41）x2z2v2t2对（1-41）式相对x和t做二维Fourier正变换，并进行算子分解得到：2~22~2~2~dd~du(duikz)(dz22kx)udz2kzu(dzikz)u0（1-42）vdz其中利用了波散关系：kx2kz22v2（1-43）式中，表示圆频率，kx和kz分别表示水平方向和深度方向上的波数。由（1-42）式得出：~~22~du（1-44dzikzui2kxu）v其中，正号代表上行波方程，负号代表下行波方程，~(kx,z,)表示波场ux,z,t的二维uFourier正变换。二．波场外推无论上行波还是下行波都可以用数学手段进行外推，外推方向可以是正向的，也可以是反向的。所谓正向外推就是根据波在当前位置上的振动情况向波的自然传播方向用计算手段预测出波场。反向外推是向波的自然传播方向的反方向上重建原来的波场。这样，对每一种波来说，无论是上行波，还是下行波，都可以进行正向外推和反向外推。对一个波场应是进行正向外推还是反向外推均有物理问题决定。1．上行波的外推把（1-44）式中的上行波方程改写为：~22dui（1-45）~v2kxdzu~1-45）式两边取积分，其积分限为z和zz，有其中u表示上行波波场。对（zz~22zzduz~iv2kxzdzu积分结果为：~2kx2z)izu(zev2（1-46）~u(z)由此得出上行波的正、反向外推式。1）上行波正向外推公式上行波的正向外推式就是向负z方向的外推公式。从（1-46）式可求出为：2kx2~~iv2zz)e（1-47）u(z)u(z根据这个公式可以计算模拟反射波的地震记录（地震图）。2）上行波反向外推公式上行波的反向外推式就是向正z方向的外推公式。从（1-46）式可得出为：2kx2z~~iv2（1-48）u(zz)u(z)e根据这个公式可以进行地震记录的向下半空间延拓，求出地下任何一点的波场，实现地震波偏移的目的。2．下行波的外推把（1-44）式中的下行波方程改写为：~2ddi2dz（1-49）~2kxdv~1-49）式进行积分，仍取积分限为z和zz，得其中d表示下行波波场。对（~2z)ikx2zd(zev2（1-50）~d(z)据此可以得出下行波的正、反向外推公式。1）下行波正向外推公式下行波的正向外推式是指沿正z方向的外推。其外推式为：2~~ikx2zv2（1-51）d(zz)d(z)e这个方程可用来模拟下行波的地震记录。2）下行波反向外推公式下行波的反向外推是指沿负z方向的外推。其外推式为：2~~z)eiv2kx2z（1-52）d(z)d(z1-52）式可用来从下行波场进行反向求源的计算工作。从上面的外推公式中可以看出，上行波的正向外推公式与下行波的正向外推公式在形式上是一致的。上行波的反向外推公式与下行波的反向外推公式在形式上也是一致的。尽管它们在形式上一致，但外推方向不同，是正相反的（从坐标系来看），同时各自的初始条件也不相同。这是应当注意的。现在来分析波场本身的条件对外推结果的影响。这里所说的波场本身的条件就是外推方程中指数项的根号项：kzk2kx2（1-53）的情况，式中k222。当kkx时，kz为正或负的实数，这时所有外推公式中存在虚v指数。说明在外推过程中波场发生相位变化。一般都能得出正确的结果。当kxk时，kz值为虚数：kzikx2k2（1-54）当把它代入各个外推公式中时，将得到如下形式的外推方程：~~kx2k2z（1-55）u(zz)u(z)e式中 z为正值，即 z zn1 zn。从（1-55）式可以看出，波场外推时只有振幅变化，而无相位变化。当指数项取负号时，外推的波场迅速衰减，称这种波为倏逝波。当指数项取正号时，外推波场迅速增大，这是一种实际不存在的波，只是进行波场计算时发生，我们称它为耗损波。在计算中要避免发生这种情况。现在我们来看一看上行波和下行波外推时产生波场损耗增大的情况。当kx k时，上行波的外推式可写为：~z)2k2zu(zkx~e（1-56）u(z)从（1-56）式可以看出，此时反向外推遇到倏逝波，正向外推发生耗损波。分别表示为：~z)~u(zu(z)e~~z)eu(z)u(z

kx2 k2 zk2x k2 z

1-57）1-58）由此可以得出结论，用上行波方程进行向下波场外推永远是计算稳定的。而用上行波方程进行正向外推就可能遇到耗损波，因此有可能是不稳定的。除非在计算中不断地把kxk的波场滤除掉。同理可求出kxk时下行波的外推式为：~d(zz)ekx2k2z（1-59）~d(z)此时也是反向外推遇到倏逝波，正向外推遇到耗损波。上面是在波数-频率域进行的波场外推的讨论。很容易把讨论的结论转用于时间-空间域，这只要把（1.3.4）式进行某种展开，代入外推式中，并经过Fourier反变换回到时间-空间域中即可得到相应的外推方程。除了用波动方程分解后的上行波和下行波方程进行波场的正、反向外推之外，还可以用Kirchhoff积分公式进行波场的正反向外推。三．地震反射波场成像在这里我们将从波动场的观点来叙述反射波成像的一般原理。关于反射波成像的概念和几何做图的各种方法请见参考文献 [1]。在介质的速度参数为已知的条件下，确定反射图象的任务就是求反射点的空间位置及其反射系数。由于我们现在还无法求出确切的反射系数值，成像的反射系数这一要求实际上是用能反映该反射点反射系数相对值的反射波振幅来表示的。因此，在目前阶段，反射成像实际上就是把地面上观测到的反射波归位到产生它的反射点上去。能做到这一点就算实现了成像。这实际就是地震偏移问题。因此地震偏移与地震成像在现阶段可以视为同一概念。为了实现地震偏移成像， 首先要进行上行波场的反向外推。 外推后求出的各点波场值，有的是来自本点的反射波，有的是该点下方许多点上的反射波。因此，要在外推波场中提取成像值。Claerbout提出下述反射波成像原则：反射面位于这些点上，其入射波的初至与反射波的产生时间相同。如图1-7所示，入射波从O点发射，反射波可设想从虚源O*点出射。在p1点上入射波先到达，反射波后到达。在p3点上设想的反射波先到达，入射波后到达。在这两点上反射波与入射波不同时，因此它们都不是反射图象的位置。只有在p2点上入射波和反射波才同时到达，而该点正好位于反射界面上。在这点图1-7入射波与反射波上的反射系数为反射波振幅u(x,z,td)除以入射波的振幅d(x,z,td)。因此，反射波成像的基本公式可写为：Map(x,z) u(x,z,td)d(x,z,td) （1-60）1-60）式没有考虑反射系数随着入射角变化的情况，它实质上是相位信息的公式。或者说，它对接近法线入射的情况时基本是正确的，能够反映反射系数在各点上的变化情况。应用（1-60）式涉及到要选择下行波的初始时间。这是一个困难问题。我们通过假设下行波是最小相位而避开这个问题。我们把td作为初始时间，则上行波和下行波对这个时间写成z变换形式为：U(Z)u0u1Zu2Z2（1-61a）D(Z)d0d1Zd2Z2（1-61b）如果D（Z）是最小相位序列，则图象函数可以表示为：Map(x,z)U(z)D(z)d（1-62）把1/D(z)展开，我们得到：Map(x,z)u01f1(ui,di)zf2(ui,di)z2dd0式中f1,f2为系数。由于z的各幂次项在单位圆上摆动，取0~2的积分限时它们等于零。在不考虑一个常系数的情况下上式等于：Map(x,z)u0d0（1-63）由于（1-60）式本质上也是等于u0d0。因此，（1-62）式与（1-60）式等价。（1-62）式比（1-60）式有优点。前者不需要考虑拾取发生的时间，而是使用u和d的全频谱。但是要求下行波是最小相位的。我们把（1-62）式的分子和分母都乘上下行波的复共轭函数D*，则有：Map(x,z)U(x,z,)D*(x,z,)d（1-64）D(x,z,)D*(x,z,)考虑到分母中的谱密度DD*不含有相位信息，因此反射图象公式可写为：Map(x,z)U(x,z,)D*(x,z,)d（1-65）根据Parseval公式，（1-65）右端项有：1 U(x,z, )D*(x,z, )d u(x,z,t)d*(x,z,t)dt2考虑到下行波是实数序列，并且不考虑积分号前的常系数，则反射图象可写为：Map(x,z) u(x,z,t)d(x,z,t)dt （1-66）当下行波是脉冲波时， （1-66）式是很精确的。但是，如果 d(x,z,t)是一个短延续长度的子波时，它只是一个很好的近似成像公式。为了验证上面反射成像公式的正确性，我们着重从相位角度来讨论它。为此取任意反射界面 f(x,z)，或可写为：Zf(x)（1-67）取一单震源在O点激发（图1-8），设激发波为脉冲，则从O点向界面发射的入射波显然可写为（设介质是均匀的，速度为v）：d(x,z,t)A0(tx2z2v)（1-68）这个任意界面的反射波方程是很复杂的。为了求出反射波的方程，我们用射线法的虚震源来导出。在弯曲界面的情况下虚震源是变动的。在二维情况下它是一条曲线。取Z f(x)界面上任意点 R的导数为：dz（1-69）f'(xR)dxR求其切线方程为：zzRf'(x)(xxR)（1-70）通过O点做垂直切线的直线方程为：xz （1-71）f'(x)联立解（1-70）和（1-71）式，则得二直线的交点为：

图1-8任意反射界面的入射波和反射波xcf'2(xR)xRf(xR)f'(xR)1f'2(xR)（1-72a）f(xR)f'(xR)xR（1-72b）zc2(xR)1f'由此求出虚震源 O*(x0,z0)的坐标方程为：x02xcu(xR)（1-73a）z02zcw(xR)（1-73b）现在我们来研究Zf(x)界面的反射波的传播方程。为此，首先要求出反射波的波前面的表达式，即要求出下列方程的包络：xu(xR)2zw(xR)2v2t2（1-74）把（1-74）式对xR微分，得xu(xR)u'(xR)zw(xR)w'(xR)0（1-75）联立解（1-74）和（1-75）方程可求出包络，其参数式为：xu(xR)vtw'(xR)（1-76a）u'2(xR)w'2(xR)zw(xR)vtu'(xR)（1-76b）u'2(xR)w'2(xR)由此求出反射波波前面的方程为（式中一撇表示导数）：tpxu(xR)r(xR)zw(xR)r(xR)vw'(xR)（1-77）vu'(xR)式中：r(xR)u'2'2(xR)w(xR)反射成像式为：Map(x,z)u(x,z,ttp)(tx2z2v)dt（1-78）从而求得成像位置的公式为：tpx2z2v（1-79）把（1-77）式代入（1-79）式，经过整理后，得v2(x)x2z22r2(x)xz2u(x)1v2(x)w(x)r2(x)xv2(x)u'(x)w'(x)u'(x)w'(x)2w(x)11u(x)r2(x)zu2(x)1v2(x)（1-80）v2(x)u'(x)w'(x)w2(x)1v102(x)式中：v(x)u'(x)，r(x)u'2'2w'(x)(x)w(x)为了直观验证（ 1-80）式，我们设任意界面为：z f(x) h则：f'(x) 0将它代入（1-80）式，求出界面成像位置方程为：4hz 4h2 0z h与原设界面方程一致。这说明反射成像的普遍公式（ 1-66）是正确的。以上是反射波成像一般原理的讨论。在地震勘探中常常要对水平叠加的地震记录进行偏移成像处理。水平叠加剖面可以看做是沿测线单点自激自收的观测剖面。在水平层均匀介质覆盖的情况下水平叠加地震剖面和沿测线进行的单点自激自收的剖面完全等价。 在非均匀介质和非水平层界面的情况下水平叠加剖面与自激自收剖面只是近似相等。 由于经济和效率的原因，地震勘探从来没有使用过自激自收观测方式。所以常常就把水平叠加地震剖面看做是自激自收的地震剖面。自激自收剖面的特点就是各个点上的入射波和反射波沿着同一条法线射线传播的。因此，在地面上观测到的反射波可以认为是在反射界面上同时激发地震波在地面上观测的结果。也就是自激自收地震剖面又与在反射界面上同时爆炸产生地震波，并以半速度向外传播，在地面上观测到的上行波是等价的。这就是 Loewenthal等人首先提出来的爆炸反射面的概念（图 1-9）。这个概念对于理解水平叠加剖面的偏移成像是很重要的。因为它比较直观地说明了这种剖面的成像原理。比前面所述的反射波成像的一般原理要容易理解得多。在这里我们引入了两种等价。一个是水平叠加剖面等价于自激自收剖面。另一个是自激自 图1-9爆炸反射面的概念收剖面等价于界面上同时激发在地面记录的上行波剖面。这种等价只是概念上的，实际上只有一种水平叠加剖面，并没有三种剖面。§1.4波场外推的Kirchhoff 积分法Kirchhoff 积分法并不直接解波动方程， 而是用数学方法来描述关于波的传播的惠更斯原理，从而求出空间上任一点波场值的。而这个波场值正好满足波动方程。Kirchhoff积分最初是为了描述波场从一个波前向传播方向上任何一点传播结果而导出的，它描述的正是一个实际物理过程，这也正是我们所说的正向外推，这样的正向外推的Kirchhoff公式对任何方向传播的波都是适用的。Kirchhoff早在1883年就证明了，从扰动区向外某点M(x1,y1,z1)传播的波的

t时刻的波场

u(x1,y1

,z1,t)

，可以从扰动区封闭表面上的波场ru(x,y,z,t )以及该波场对时间和表面法线方向的导数通过积分式求出来。因此要假定cu(x,y,z,t)在封闭面上和封闭面内有直至二阶导数的连续性。Kirchhoff 利用了格林定理：(u2vv2u)dVuvvudS（1-81）VSnv式中u，v为具有直至二阶导数的连续性的函数。在我们所讨论的情况下，u取为波场函数，v1。其中R2(xx1)2(yy1)2(zz1)2。R1-10（a）,（b）那样的封当把观测点用包含有波前面在内的封闭曲面包围起来，如图闭时，这样的封闭面S和它所包围的体积V作为（1-81）式的积分限，经过一定的推导后得出M(x,y,z)点的正向外推波场为：u(x1,y1,z1,t)11Ruu11u4SvRntnRRdS（1-82）n这里n的方向取封闭表面的外法线方向。如果把观测点M移至封闭面外，则有：u(x1,y1,z1,t)0（1-83）（1-82）式中u(,,,R)，v为传v播速度。在实用中所取的计算封闭面一定要把需计算波场的那一点或所有计算点包围在图1-10计算正向外推波场的封闭面内。同时，计算点M一定要满足处在波前推进方向的前方这一物理要求。否则，既使能够计算出某个数值，也不代表那一点的波场。1-82）式就是著名的Kirchhoff积分。它描述了物理波场传播的过程，也满足齐次波动方程，是它的积分形式解。对我们来说，也可以称它为正向外推公式，它既可以用于上行波的正向外推，也可以用于下行波的正向外推。只是外推方向由物理条件而定。这一点与前一节把波动方程分解为上行波方程和下行波方程后它们的外推公式在形式上是相同的情况完全符合。还要指出，（1-82）式的v是介质的传播速度，在积分式的推导中令其为常数。因此Kirchhoff积分也只满足均匀介质的情况。下面讨论用 Kirchhoff 积分进行波场反向外推问题。对偏移来说，需要利用 Kirchhoff积分进行反向外推。 这时，我们所取的封闭体积 V应当在波前传播方向的反方向， 计算点(x1,y2,z3)就在这个封闭体内。根据格林定理同样可求出形式上相同的反向外推的Kirchhoff积分式：11Ruu11uu(x1,y1,z1,t)vRntnRRdS'（1-84）4S'n式中的[[u]]不再是推迟场，而是超前场(,,,R)。uxyztv（1-84）式为用于波场反向外推的Kirchhoff积分式。它可用于上行波的反向外推，也可用于下行波的反向外推。当然，这种外推与正向外推不同，它不代表一个物理过程，而只是一种重建波场的计算过程。1.5Born近似和Rytov近似本节首先讨论 Born近似，然后讨论 Rytov基本特点。为以下讨论的方便，给出如下的目标函数

近似，最后讨论

Born

近似和

Rytov

近似的O(r)v01（1-85）v(r)则有uk02u2k02O(r)u(r)k02O2(r)u(r)（1-86）在波场弱散射的假设条件下讨论Born近似和Rytov近似表达式。一．Born近似设波函数具有形式解uum（1-87）m0其中，um关于O(r)是m阶的，把式（1-87）代入式（1-86）比较O(r)的阶数，得um满足方程组u0k02u00（1-88）uk2u12k02O(r)u0（1-89）10222（1-90）um1k0um12k0O(r)umk0O2(r)um1m=1,2,3,。从上面的方程组可得各阶Born近似u1(r)2k02drG(rr)O(r)u0(r)（1-91）式中，u0为入射波场，um1(r)2k02drG(rr)O(r)um(r)k02drG(rr)O2(r)um1(r)（1-92）特别是u2(r)2k02drG(rr)O(r)u1(r)k02drG(rr)O(r)u0(r)（1-93）二．Rytov近似设波函数ue，则方程（1-86）可写为(222（1-94）)2k02k0O(r)k0O2(r)设相位有级数解，即m（1-95）m0其中，m关于O(r)是m阶的。将式（1-95）代入式（1-94）。比较O(r)的阶数，得到012

(0)2k0200（1-96）22k020（1-97）01O(r)202(1)2k02O2(r)0（1-98）这里，方程（ 1-96）对应方程（ 1-94）在O（r）=0时的情形，即无扰动情形。令u0 e0，则入射波u0满足u0 k02u0 0 （1-99）设入射波u0为u0e0，0ikr（1-100）利用方程（1-97）、（1-99）和自由空间的 Green函数G，有2k02r)O(r)u0(r)（1-101）1drG(ru0(r)同理，利用式（1-98）和（1-100），得到21drG(rr)[k02O2(r)(1)2]u0(r)（1-102）u0(r)一般地，对m3，有1m1mdrG(rr)[pmp]u0(r)（1-103）u0(r)p1其中，m称为第m阶的Rytov近似。在上面表达式的推导过程中，利用式（1-100），则关于的一阶Rytov近似的方程式（1-97）就可以成为(1u0)22（1-104）k01u02k0O(r)u0比较式（1-104）和式（1-89），除相差一个因子u0外，方程形式完全相同，但这并不意味着Born近似和Rytov近似是等同的。应该指出，Born近似利用的是波场的振幅，而Rytov近似则利用的是波场的相位。利用波场相位的Rytov近似在某些假定之下可与旅行时反演相比较。三．Born近似和Rytov近似的基本特点在研究Born近似时，假设总波场是由入射波场u0与散射波场us叠加而成，即uu0us（1-105）这时方程（1-85）变为usk02us2k02O(r)u(r)k02O2(r)u(r)（1-106）利用自由空间的Green函数G，得到散射波场us为us(r)2k02drG(rr)O(r)(u0us)（1-107）Born近似是在假定usu0之下，在目标函数O(r)内测不到散射场，用u0代替u u0 us，即us(r)2drG(rr)O(r)u0(r)（1-108）2k0但是当目标O(r)的尺度较大时，目标O(r)内的波场uu0us是不能简单地用入射波u0来代替的。如假定目标O(r)是半径为a的圆柱体，在n(v0v1)/v1条件下，沿k方向传播的平面波u0eikr在圆柱体内部并不等于入射波u0，特别在通过圆柱体对称轴的、柱体内部的路径上的波场为u0ei(1n)kr，尤其是通过圆柱体的波场的相位变化是2an4na（1-109）v0其中圆频率2/，为波长。Born近似有效的一个必要条件是入射波通过目标时的波长变化小于，即4na/或na1（1-110）4式（1-110）表明，Born近似的有效区域依赖于目标尺度的大小。在目标尺度较小时，Born近似有效。关于Rytov近似，令相位0s（1-111）则式（1-97）变成s20s(s)22k02O(r)k02O2(r)0（1-112）求解方程（1-112），得s1drG(rr)[2k02O(r)k02O2(r)()2]u0(r)（1-113）u0Rytov近似是在假定(s)222(r)22（1-114）koO2koO(r)2koO(r)下进行的，这就要求( s)2 ko2n即n1（1-115）(s)242式（1-115）表明，Rytov近似对目标尺度的依赖性不像 Born近似那样紧密； Rytov近似的有效性依赖于波长下复相位的变化情况。换言之， s越小，Rytov近似效果越好，这个特点是Born近似所没有的。§1.6DeWolf近似、薄板近似和屏近似基于DeWolf近似和相屏单程传播，我们可以在地面反射地震观测中计算一次反射波或反向散射波的合成地震记录。具体实现时一般有两种不同的算法：一种是直接利用多次正向散射-单次反向散射（ MFSB）近似，其中计算反向散射波的步长等于网格间隔，但计算正向传播的步长使用远远大于网格间隔的屏间隔；另一种是对正向和反向散射波的屏近似，这样两种步长都能采用屏间隔，从而大大提高了计算效率。在复杂不均匀介质（尤其是3-D介质）中，快速模拟方法和算法是发展复杂构造地震成像、反演和解释方法的关键。有限差分和有限元法从原理上能模拟在任意不均匀介质中波的传播，但计算工作量大、效率低。在大型3-D弹性波问题的情况下更是不可想象。下面给出了一种基于多次正向散射-单次反向散射（Multiple-forescatteringsingle-backscatteringMFSB）近似，即DeWolf近似，在地面反射地震观测中计算反向散射场的方法。为快速实现该方法，我们给出了双域算法。当不均匀体尺度大于地震波主波长时，这一理论可进一步通过薄板近似和屏散射近似得到合理近似和简化。薄板近似和屏近似能减少计算工作量、提高计算效率。当介质体内的不连续性不是很强或不均匀体的参数扰动不太大时，在不均匀体和谐振散射之间的混响通常可以忽略不计。但正向散射的累加效应一般不能忽略掉。实际上，对于大型不均匀体介质或长传播距离，多次散射对正演模拟和反问题来讲是很重要的。在这种情况下，Born近似不满足，但可应用 DeWolf近似。一．Lipmann-Schwinger 方程从标量波动方程出发2( 2 )u(x) 0 （1-116）v2(x)定义F(x)v021s2