1、 叠前深度偏移成像叠前深度偏移成像理论与方法理论与方法一一.Kirchhoff积分法叠前深度偏移积分法叠前深度偏移二二.波动方程法叠前深度偏移波动方程法叠前深度偏移 进入九十年代，地震勘探逐渐向中深层和复杂构造区域的精细勘探精细勘探发展。在这种形势下，地震偏移作为一种重要的地震数据处理手段，迫切需要新的发展以解决复杂地质体成像问题。实践证明常规的叠后时间偏移在复杂介质下，很难取得好的成像效果。于是，人们提出使用叠前深度偏移叠前深度偏移。常规的叠前深度偏移算法可以分为两类：一类是基于射线基于射线追踪的追踪的Kirchhoff积分算法积分算法，另一类是基于波场延拓的波动方基于波场延拓的波动方程解法
2、程解法（如有限差分偏移、相移偏移和广义屏偏移等）。由于Kirchhoff积分偏移计算效率非常高，它在石油工业领域得到了广泛的应用。然而，利用程函方程求解初至旅行时的Kirchhoff积分偏移在成像效果上不如一些波动方程偏移方法。这主要是因为通过焦散区焦散区(在光学中，相位函数的光滑性或单值性被破坏的点的集合，此时振幅趋于无穷大，出现奇点）的射线由于存在多路径多路径而产生的后至波可能携带重要的能量，这部分信息对深度成像而言是不应当忽略的。尽管基于完全Green函数射线理论（例如在计算中考虑多至走时多至走时及相应的振幅）的Kirchhoff偏移的成像质量有所提高。但不幸的是，这类算法非常复杂，计算
3、代价特别大，且对焦散区并非完全有效。第一节第一节 KirchhoffKirchhoff积分法叠前深度偏移积分法叠前深度偏移 Kirchhoff积分法叠前深度偏移已在实际生产中应用了多年，并解决了不少复杂构造的成像问题。Kirchhoff积分法的关键是绕射旅行时的计算绕射旅行时的计算，目前常用的计算方法是射线追踪法和有射线追踪法和有限差分法限差分法。射线追踪法计算绕射旅行时可分为常速法和变速常速法和变速法法，常速法很简单；变速法将在下面介绍。有限差分绕射旅行时计算是基于费马原理基于费马原理，可在直角坐标系或球坐标系实现。一一Kirchhoff积分法叠前深度偏移成像公式积分法叠前深度偏移成像公式
4、为适应复杂地质构造和岩性成像实现保幅处理保幅处理。可采用如下的考虑传播效应的Kirchhoff积分法叠前深度偏移成像公式rsrrssrsrrrsdxtxxxxxxuxxxAxxnxxR);,();(;();();();(（1-1-1）式中 是记录面（线）；分别表示成像点、震源点和接收点；分别表示震源点到成像点和成像点到接收点的旅行时；A是几何扩散因子（振幅加权因子）；是记录面 的外法线方向；u是记录波场；R是反射系数。rsxxx,rs,n二二维有限差分绕射旅行时计算方法二二维有限差分绕射旅行时计算方法 在Kirchhoff叠前深度偏移中，绕射走时的计算精度直接影响了成像精度。为此，我们利用Ei
5、konal方程和费马原理方程和费马原理在规则网格上进行绕射走时的有限差分计算，无需走时的内插。二维空间的Eikonal方程为 222),()()(zxsztxt（1-2-2）其中，(x,z)是空间坐标，s是慢度，（1-2-2）式的两个偏导数可由有限差分近似)(213120tttthxt)(213210tttthxt（1-2-3a）（1-2-3b）其中，分别是A，B1，B2，C1处的旅行时，Bi处的旅行时由下式确定 4210,tttt)(2ABiissht（1-2-4）把（1-2-3a）和（1-2-3b）代入（1-2-2）得 为获取待求行或列上第一个点处的旅行时，利用有限差分得出如下计算公式 2
6、12203)()(2tthstt212203)(25.0)(tthstt（1-2-5）（1-2-6）其中，分别是内行或内列上旅行时的相对极小值、两侧的旅行时、和待求行或列上的旅行时。4210,tttt（1-2-5）和（1-2-6）式是基于平面波近似得出的有限差分旅行时计算公式。另外也可基于局部圆弧波前近似基于局部圆弧波前近似求取旅行时计算公式。考虑到计算效率和精度要求，可采用基于平面波近似的有基于平面波近似的有限差分旅行时计算公式限差分旅行时计算公式。图1-1.有限差分旅行时计算示意图 AB1C1C4C3C2B2B3B4 我们可按图1-1的方阵逐步扩展计算旅行时。其具体的做法是：第一，先由（1
7、-2-4）式计算出“十”形端点的值（假设中心点已知，各点慢度已知，网格间距为h），然后由（1-2-5）式计算拐点的值。第二，第一圈计算完后，任取该圈的一边找出极小点，然后以此点为起点由（1-2-6）式计算出下一条边的起点，再利用（1-2-5）式把整条边上的值计算出来。第三，用同样的方法把四条边上的值都计算出来。第四，依此类推，可以把整个模型上所有网格点的到达时计算出来。利用有限差分法计算绕射走时的利用有限差分法计算绕射走时的Kirchhoff叠前深度偏移叠前深度偏移方法的优点是运算速度快，能够处理较复杂横向速度变化的情方法的优点是运算速度快，能够处理较复杂横向速度变化的情况。并且能解决射线追踪
8、所产生的问题况。并且能解决射线追踪所产生的问题。三迎风有限差分三维旅行时计算三迎风有限差分三维旅行时计算 Kirchhoff积分法三维叠前深度偏移的核心是复杂介质中的积分法三维叠前深度偏移的核心是复杂介质中的旅行时计算旅行时计算。三维旅行时算法的效率和精度直接决定了成像方法的应用范围和效果。考虑到偏移速度场的复杂性和三维叠前深度偏移的大计算量，必须采用稳健高效的三维旅行时计算方法。目前大部分的Kirchhoff积分法三维叠前深度偏移仅利用初至波旅行时。实质上，波前面有时会数次通过速度空间中的某一点，出现多值旅行时多值旅行时，使得初至波旅行时与携带最大能量的波前的到达时不一致，此时仅用初至波旅行
9、时成像显然是不合适的。然而，利用初至波旅行时的Kirchhoff积分法偏移实现简单、计算效率高，且成像精度能满足大多数情况下的实际要求。因此，可采用稳健高效的迎风有限差分三维旅行时计算方法采用稳健高效的迎风有限差分三维旅行时计算方法。目前常用的三维旅行时计算方法大致分为两类：射线追踪射线追踪旅行时计算和直接解程函方程的旅行时计算旅行时计算和直接解程函方程的旅行时计算。射线追踪需要从程函方程导出射线方程，然后导出射线坐标、方位、和旅行时方程。依此，可用解初值问题的试射法或解边值问题的两点射线追踪法。当地下存在复杂地质结构时，这两种射线追踪法的应用会遇到如下三个困难三个困难：一是入射角的微小变化会
10、导致最终结果的很大变化；二是由于阴影区和焦散区的存在，即使从源点出发的射线角度均匀且密集，得到的旅行时场也会存在很多大小不一的空白区和多值区，此时一般的插值方法是不能解决问题的；三是速度场变化剧烈时，无法得到正确的射线路径。直接解程函方程的三维旅行时计算法能克服上述困难，这种计直接解程函方程的三维旅行时计算法能克服上述困难，这种计算法会因不同的解法而适应不同的速度场算法会因不同的解法而适应不同的速度场。鉴于计算方法的稳定性和效率以及对复杂速度场的适应性，可采用直接解程函方采用直接解程函方程的稳健高效的迎风有限差分三维旅行时计算方法程的稳健高效的迎风有限差分三维旅行时计算方法。迎风差分格式就是考
11、虑到波的传播规律构造出的差分格式迎风差分格式就是考虑到波的传播规律构造出的差分格式。该方法较高的计算效率正是源于迎风差分格式，因为它不需要搜索旅行时的相对极小值点或多次计算同一点的旅行时来保证计算过程满足因果律、计算旅行时满足费马原理。另外，迎风有限差分法适合于向量并行巨型机；改进型（无量刚形式）的有限差分方程明显提高了计算方法的稳定性；下面给出的稳定性条件自适应地决定了用于外推的径向（三维旅行时的计算采用了球坐标系）步骤；计算网格的数目、大小、形状、和边界会随径向半径和速度场的复杂性自适应选取；考虑到原坐标系是笛卡尔坐标系，旅行时计算网格是笛卡尔网格，因此球坐标系下的计算网格的边界截断也是自
12、适应的。在笛卡尔坐标系下的程函方程为 2222),()()()(zyxsztytxt（1-2-7）球坐标系下的程函方程为 22)sin(1212),()()()(22zyxstrtrrt（1-2-8）其中t是旅行时，S是慢度。定义 tvturtw,（1-2-9）考虑到计算方法的稳定性，令 rvvruu,（1-2-10）基于有限差分近似，最后得出 2/12222sin),(),(),(),(rrvrrurrsrrw（1-2-11）经推导有如下的CFL稳定性条件 vuwrr22sinsin21（1-2-12）其中 rrr/)(上述方程在全部 和 上可实现向量化和并行化实现向量化和并行化。在满足（1
13、-1-12）式的条件下，利用上述方程可实现三维旅行时场的实现三维旅行时场的计算计算。四四Marmousi模型叠前深度偏移试算模型叠前深度偏移试算 为了检验该叠前深度偏移算法，可采用美国SEG用于检验叠前深度偏移成像方法性能的Marmousi模型模型。模型速度场如图1-2所示，维数为497x750,速度场道间距为12.5m，最大深度为3000m，深度采样间隔为4.0m。该速度场的横向变化非常剧烈，因而很适合用于叠前深度偏移方法试验。该模型的正演模拟炮记录为SEG/Marmousi模型数据。共有240炮，每炮96道接收，为右边放炮方式，最小炮检距为-200m，最大炮检距为-2575m，道间距为25
14、m，道长为750个样点，时间采样率为4ms。Kirchhoff叠前深度偏移结果见图1-3。结果表明：尽管在反射振尽管在反射振幅上有些失真，但主要构造特征还是清晰可见幅上有些失真，但主要构造特征还是清晰可见。图1-3.基于二维有限差分绕射旅行时计算的M a r m ous i 模型Kirchhoff叠前深度偏移结果 图1-2.Marmousi模型速度场 第二节第二节 波动方程法叠前深度偏移波动方程法叠前深度偏移 基于波动方程的递归偏移方法相对基于波动方程的递归偏移方法相对Kirchhoff积分偏移有着积分偏移有着得天独厚的优势得天独厚的优势。首先，它们是从全方程导出从全方程导出的，不是基于高频近
15、似的渐进解。因此，波动方程偏移方法潜在地比较精确、稳健。其次，当向下延拓方法向下延拓方法可用于波场延拓而又不增加计算维数时（例如零偏移距数据），这类方法比Kirchhoff积分法的效率还高。波场延拓算子是递归偏移方法的关键波场延拓算子是递归偏移方法的关键，它们一般是从单程波从单程波方程推导出来方程推导出来的。通常大家较为关心的单程波方程偏移的典型算法是有限差分方法和付立叶方法有限差分方法和付立叶方法。前者容易处理速度的横向前者容易处理速度的横向变化变化，但其缺点在于存在频散和成像倾角限制。后者后者不存在频散且对水平层状介质能精确成像，不过，它只适用于水平层状只适用于水平层状地层地层。若采用共轭
16、梯度法优化傍轴近似方程的系数采用共轭梯度法优化傍轴近似方程的系数，可给出一种频率空间域的有限差分波场延拓方法；在反射系数估算意义下，可推出叠前深度偏移的成像公式。并且发现，当深度偏移算子中的折射项方程采用时移处理折射项方程采用时移处理，并且和成像计算中关于坐标变关于坐标变换复原的时移处理合并换复原的时移处理合并在一起时，计算可大大简化计算可大大简化。为了利用付立叶方法的优势，许多地球物理学家提出在双域双域（即频率波数域和频率空间域）进行地震波成像处理，先后提出了分步付立叶偏移、付立叶有限差分偏移等算法分步付立叶偏移、付立叶有限差分偏移等算法。下面详细介绍这两种偏移算法的基本原理，分析波场延拓算
17、子的相对误差，而且还对付立叶有限差分偏移算子给出优化改进的思路。波的传播和成像问题是有着密切联系的。De Hoop（1988）等人为了研究波在随机介质中的传播问题，基于波的散射理基于波的散射理论提出了早期的屏方法论提出了早期的屏方法。Wu R S.&Huang L.J.等在近十年来发展了广义屏方法广义屏方法，并且用于叠前深度偏移成像。下面对比各种广义屏方法，而且讨论它们之间的相互关系和稳定性问题。首先在Kirchhoff-Helmholtz积分意义下，从方法原理上把前面提到的几种递归偏移方法统一起来几种递归偏移方法统一起来，并阐述它们之间的相互联系相互联系。随后借助于脉冲响应和复杂模型偏移的数
18、值计算结果，进一步对比频率空间域有限差分法、付立叶有限差分法、分步傅立叶（相屏）法和广义屏方法各自的特点各自的特点，这些可为我们解决不同的实际问题选择合适的偏移方法提供依据。最后对波动方程三维叠前深度偏移的思路给出探索性观点探索性观点。概述概述 地震波成像地震波成像在油气勘探中占据重要位置。它的作用是使反射波或衍射波返回到产生它们的地下位置，从而得到地下地质构造的精确成像。从二十世纪60年代偏移过程由计算机实现以来，已从常规偏移即叠后时间偏移叠后时间偏移发展到了目前的叠前深度偏移叠前深度偏移。偏移方法的研究和应用是受油气勘探的实际需求驱动的，同时它又受到人们对偏移成像的认识程度和计算机处理能力
19、的制约。常规偏移（即叠后时间偏移）在以往的油气勘探过程中起到了重要作用，但随着勘探难度的提高，在构造较为复杂的地区，基于常构造较为复杂的地区，基于常规偏移的处理方法再也难见成效规偏移的处理方法再也难见成效。究其原因，一方面是由于常规处理是先叠加后偏移一方面是由于常规处理是先叠加后偏移，水平叠加过程受层状介质假设制约受层状介质假设制约，在复杂地质构造条件下，这种叠加过程很难实现同相叠加，这样会对波场产生破坏，所以用这种失真了的叠后数据去进行偏移处理难以取得好的成像效果就很自然了。为了克服非同相叠加给后续偏移带来的麻烦，人们提出使用叠前偏移叠前偏移，即先偏移处理使波场归位，再把同一地下点的偏移波场
20、相叠加。另一方面是由于时间偏移是建立在均匀介另一方面是由于时间偏移是建立在均匀介质或水平层状介质的速度模型的基础上的质或水平层状介质的速度模型的基础上的，当速度存在横向变化，或速度分界面不是水平层状的情况下，常规偏移不能满足斯奈耳定律，因此不能进行正确的反射波的偏移成像。为了解决这个问题，出现了深度偏移深度偏移。这样，叠前深度偏移就可以弥叠前深度偏移就可以弥补常规偏移的不足补常规偏移的不足。但是，卓有成效的叠前深度偏移仍然处在探索中。大家知道，叠前偏移的概念早在70年代中期就提出来了，但由于叠前记录的信噪比较低，偏移的初始模型又很难选准，加之当时的计算机无法承受叠前偏移较大的计算量，直到90年
21、代叠前偏移才开始尝试应用于油气勘探地震数据的精细处理中。常见的叠前深度偏移算法可以分为两类：第一类是基于绕射扫基于绕射扫描叠加原理的描叠加原理的Kirchhoff积分算法积分算法，另一类是基于波动方程的偏基于波动方程的偏移方法移方法（如有限差分偏移方法、付立叶偏移方法等）。起初出于对计算机处理能力方面的考虑，人们首先想到的是使用一种快捷的算法来实现叠前偏移过程。由于在扫描叠加偏移原理基础上，基于射线追踪的Kirchhoff积分方法正好具有这方面的优势，它很快成为了叠前偏移方法研究的重点，并且很快就推出了应用软件。工业应用中Kirchhoff积分法叠前偏移在某些地区有时会取得常规偏移难以得到的成
22、像效果。但是，该类方法本身存在很明显的缺陷。例如射线追踪前需要对速度场进行平滑，在速度分布过于复杂的区域，会出现焦散或阴影焦散或阴影，计算出来的旅行时场也就不准确。后来，为了提高旅行时场的精度，又发展出了有限差分法直接求解程函方程的Kirchhoff积分偏移方法。但一般仅能计算初至旅行时，无法处理在复杂速度场中存在的多至走时现象，从而影响Kirchhoff积分偏移在复杂地质体（如盐丘、超覆等）的成像效果。从近十年Kirchhoff积分偏移实际应用可以证明这一点。如果应用完全射线理论的Green函数，在计算时求解所有的到达时和相应的振幅值，可以改善该方法的成像质量，但其计算效率又会大打折扣。由于
23、波动方程偏移方法基本不存在Kirchhoff积分偏移这些困难，因此，近年来人们对波动方程叠前偏移进行了大量的研究。波动方程叠前深度偏移成像解决的是强横向变速条件下复杂地质体的地震波成像问题。近年来一些专家学者在已有方法的基础上提出几种效果较好的叠前偏移方法。大体上讲，这些方法基本上主要分为两类，其一为有限差分偏移方法有限差分偏移方法，另一类为付立叶偏移方法付立叶偏移方法。两类偏移方法各有特点，它们可以分开使用，也可以联合使用联合使用（如所谓的混合偏移混合偏移）。波动方程有限差分算法借助于差分计算，把速度、密度等介质参数的影响体现在差分计算的矩阵方程中，因此，这类方法能自动适应速度场的任意变化。
24、这类方法中比较典型的是波动方程有限差分法逆时偏移波动方程有限差分法逆时偏移。它基于双程波方程双程波方程的，既能适应速度场任意变化，又不存在倾角限既能适应速度场任意变化，又不存在倾角限制制。但它具有一些难以克服的缺点，如计算效率非常低、对计算机内存要求太高等等，这就迫使人们尽量使用单程波方程进尽量使用单程波方程进行波场延拓和成像行波场延拓和成像。为了实现深度偏移，通常按Hatton(1981)提出的思路，通过引入参考速度引入参考速度，可以把上、下行波方程分裂为上、下行波方程分裂为绕射项和折射项绕射项和折射项。绕射项方程绕射项方程实际上就是时间偏移方程时间偏移方程，它可使绕射波收敛使绕射波收敛；折
25、射项方程折射项方程的作用是校正由于横向变速引起的校正由于横向变速引起的时差时差。上、下行波方程是由双程波方程分解近似得到的，它们存在偏移倾角限制，即在某种传播角度范围内才能较准确地描述地震波的传播过程。于是张关泉（1986）、Holberg,O.(1988)等人提出对单程波方程的系数进行优化，尽量提高低阶方程的成像精度。王华忠（1997）用共轭梯度法共轭梯度法优化常规的方程，并提出了一种在时空域进行的有限差分深度偏移算法。付立叶偏移算法付立叶偏移算法一般借助快速付立叶变换来进行波场延拓计算，最早期的付立叶偏移方法应当是相移偏移相移偏移。这种方法不存在偏移倾角限制，没有频散，而且计算效率非常高。
26、但它是基于层内常速假设基于层内常速假设条件的，不能处理横向速度变化地层成像问题。为了改进相移偏移方法，人们提出在相移的基础上增加对横向速度扰动引起的时差的校正处理。这些方法的基本思路是进行速度场分裂速度场分裂，即把复杂的介质速度场分裂为“常速背景常速背景+层层内变速扰动内变速扰动”，然后针对分裂后的速度场分别进行波场延拓处理。例如，Stoffa(1990)提出分步付立叶偏移方法分步付立叶偏移方法，常速度背景常速度背景对应的波场用相移法对应的波场用相移法进行波场延拓，然后针对变速扰动引起的变速扰动引起的时差进行时移校正时差进行时移校正。该方法适用于非剧烈变速情况下的深度偏移成像。Ristow&(
27、1994)提出的付立叶有限差分深度偏移成像方法付立叶有限差分深度偏移成像方法是在分步付立叶方法的基础上，加上一个有限差分项对二阶以上速加上一个有限差分项对二阶以上速度扰动引入的时差进行校正度扰动引入的时差进行校正。该方法在剧烈变速情况下也具有非常好的成像效果。90年代初，Wu R S 与de Hoop等在波动方程Green函数解法基础上，通过一系列近似处理手段，发展成了较实用的广义屏算法广义屏算法。这类算法既可用于研究波（声波或弹性波）的传播问题，又可用于地震波场成像。该类方法的基本思路是，将速度场分解为层内常速背景和层内变速扰动速度场分解为层内常速背景和层内变速扰动。对背背景场景场相当于解常
28、速的声波方程，可通过相移实现相移实现；对变速扰动，可认为这种非均匀性相当于散射源（二次源），入射波场作用于这些散射源上，由此产生散射波场。一般假设波场延拓层厚度较小，在薄板近似薄板近似处理下，对散射场计算式的积分核再采用不同的近似方法，如Born近似、近似、De Wolf近似或近似或Rytov近似近似等等，可以得到不同的散射波场表达式。随后在此基础上，Wu R.S和Huang L.J.等人把该类方法发展成了广义屏、相屏以及局部Born近似的屏方法和局部Rytov近似的屏方法。接着他们针对上述方法的一些不足之处做了进一步的扩展和改进，提出了一系列新的屏算子，如扩展的局部Born近似和Rytov近
29、似的屏方法、拟线性局部Born近似的屏方法，等等。我们在第一部分第一部分简要介绍基于共炮道集的波动方程叠前深度偏移的基本思路，接着从第二部分第二部分把优化系数的旁轴近似方程转入频率空间域，用有限差分法进行波场延拓，并在反射系数估算意义下推导出叠前深度偏移的成像条件，接着进行数值试验（脉冲响应计算、凹槽模型叠后深度偏移和Marmousi模型叠前深度偏移测试）。第三部分到第五部分第三部分到第五部分依次叙述分步傅立叶法、付立叶有限差分法和广义屏法偏移的基本原理、实现方法和数值试验结果。第六部分第六部分阐述前面几种叠前深度偏移方法的相互联系，并把它们从脉冲响应、偏移成像效果、稳定性和计算效率等方面做比
30、较详细的对比，然后对三维波动方程叠前深度偏移做初步的探讨。第一部分 基于共炮集的波动方程叠前深度偏移概论基于共炮集的波动方程叠前深度偏移概论 基于共炮集的波动方程叠前深度偏移思路是基于共炮集的波动方程叠前深度偏移思路是，首先对每一炮进行单炮偏移成像，然后再把各炮成像结果在对应地下位置上叠加，从而得到整个成像剖面。对于每一炮对于每一炮，标准的波动方程叠前深度偏移可以分为三步分为三步：震源波场的正向延拓、炮集记录震源波场的正向延拓、炮集记录波场的反向延拓和应用成像条件求取成像值波场的反向延拓和应用成像条件求取成像值（Clearbout,1971）。为了方便叙述基于共炮集的波动方程叠前深度偏移的基本
31、过 程，我 们 引 入 基 于 单 程 波 方 程 的 波 场 传 播 算 子（Berkhout,1987）。以频率域二维波场为例，如果对震源波场 和 炮集记录波场做如下定义：);z,x(su);z,x(sv 1、：它是炮点 处频谱为 的点源激发产生的震源波场，有 （1-1）2、：它是点 处激发，排列接收到的记录波场，该波场可以写成：（1-2）这里 含有一非零道，即在接收点 处的记录道，它满足：（1-3）);0,x(sus)(f)(f)sx();0,x(su);0,x(svsdr);0,x(r,sv);0,x(sv);0,x(r,svr);0,x(sv)rx();0,x(r,sv3、：它表示在
32、深度 处的正向延拓波场，如果引入表征波场从地面传播到深度 的传播算子 ，则有：（1-4）4、：它表示记录波场 在深度 的反向延拓波场：（1-5）其中 为记录波场的反向传播算子。因为波场传播算子 描述上行波从深度 到地面的传播过程，故 描述了（向上传播的）记录波场从地面到深度的反向延拓过程（Berkhout,1987）。);z,x(su0z z)z0(W);0,x(su)z0(W);z,x(su);z,x(r,sv);0,x(r,svz);0,x(r,sv1)0z(W);z,x(r,sv1)0z(W)0z(Wz1)0z(Wz 和 分别称为下行波和上行波的深度外推算子。实际计算过程中，逐层实现上、
33、下行波的波场延拓和求取成像值。层 内波场延拓如图1所示。)z0(W10z(W)zz,z();z,x(su);z,x(r,sv)z0(W1)0z(Wz);zz,x(su);zz,x(r,sv图1 叠前深度偏移波场延拓示意图 注意到所有波场空间上离散分布在采样间隔为 的地震道上，故空间 函数可以用长度为 、振幅为 的加窗函数表示，而且以上公式中的积分可以用离散求和替代。用传统的偏移公式，我们可以得到对点 处震源在 处的单一记录道的叠前深度偏移结果具有如下形式：（1-6）式中 表示复共轭，表示取复数的实部。这种频率域的成像公式相当于时间域震源波场正向延拓值与记录波场反向延拓值的互相关。由（1-6）式
34、，可以得到共炮集数据的叠前深度偏移成像公式：rrr/1sr);z,x(r,sv);z,x(su)z,x(r,sm （1-7）如果采用的是有限差分算法进行波场延拓计算，则上式中关于 的求和并非是显式的，它一般在对整个单炮记录波场的反向延拓中自动实现。因而，共炮集记录叠前深度偏移公式可表述为：（1-8）这里 为炮 整个记录波场的反向延拓波场值。rr);z,x(r,sv);z,x(su)z,x(sm);z,x(sv);z,x(su)z,x(sm);z,x(svs 从计算角度而言，成像过程是很简单的步骤，波场外推算波场外推算子的数学形式和计算实现才是地震波偏移成像的核心子的数学形式和计算实现才是地震波
35、偏移成像的核心。目前所有基于波动方程的波场延拓算子不外乎有：波动方程有限差分波动方程有限差分波场延拓算子和付立叶波场延拓算子波场延拓算子和付立叶波场延拓算子。前一类算子既可以在时间空间域又可以在频率空间域用有限差分方法有限差分方法实现波场延拓计算，只是在波动方程在频率空间域频率空间域的形式更简单，差分计算和成像更方便。第二部分讲述的频率空间域有限差分法叠前深度偏移成像即是这类方法。后一类算子中最熟悉也是最简单的要数相移算子，它在频率波数域计算实现。然而，当速度横向变化时，关于空间坐标的付立叶变换不再成立，这迫使我们在处理横向变速介质中的波的传播和成像问题时退回到空间域。因此，在解决这些问题时，
36、一般基于速度场分裂，对背景场和扰动场分开处理，在频率波数域和频率空间域（或双域daul-domain）交替进行波场延拓计算。如第三、四、五部分讲述的分步傅立叶方法、付立叶有限差分方法和广义屏方法分步傅立叶方法、付立叶有限差分方法和广义屏方法都是采用这类双域延拓算子双域延拓算子。第二部分 频率空间域频率空间域有限差分法波动方程叠前深度偏移有限差分法波动方程叠前深度偏移 这一部分介绍一种基于共炮道集的波动方程叠前深度偏移基于共炮道集的波动方程叠前深度偏移算子算子，它可以用来进行上下行波的波场延拓上下行波的波场延拓。在基于反射系数估算的成像条件下，可实现叠前深度偏移成像实现叠前深度偏移成像。由于算子
37、为优算子为优化系数的旁轴近似方程化系数的旁轴近似方程，它具有方程阶数低方程阶数低且能对陡倾角成像对陡倾角成像的特征，并采用有限差分法波场延拓，能适应速度的任意横向适应速度的任意横向变化变化，故该算法可使复杂地质体精确成像使复杂地质体精确成像。因为一切计算在频率-空间域进行，相对于纯粹的时间空间域有限差分算法有计计算效率高、成像方便算效率高、成像方便的优点。下面进行了常速和变速介质的脉脉冲响应冲响应测试，并对凹陷模型凹陷模型爆炸反射记录进行了叠后深度偏移，对Marmousi模型模型进行了叠前深度偏移处理。试验结果表明该算法是比较优越比较优越的，且具有可供进一步挖掘的潜力进一步挖掘的潜力。一一 概
38、概 述述 叠前深度偏移是一种对复杂地质构造成像的重要的和有效叠前深度偏移是一种对复杂地质构造成像的重要的和有效的工具的工具。已有的方法不是在时间域时间域就是在频率域频率域进行。频率域算法的波场延拓既可以在波数域进行，也可以在波数域和空间域交替进行。一种比较典型的叠前时间域算法叠前时间域算法是波动方程逆时波动方程逆时偏移方法偏移方法（E.G.Chang&McMechan 1990），它一般通过有限差分法实现波场延拓。另一种时间域算法是基于射线追踪的基于射线追踪的Kirchhoff积分方法积分方法（E.G.Hu&McMechan,1986）。由于全方程逆时偏移的有限差分法是非常耗时非常耗时的，故其
39、叠前深度方法在目前生产中是难以承受的。而基于射线追踪的Kirchhoff叠前成像是一种高频近似高频近似方法，且在利用射线追踪计算绕射旅行时时，常需对速度场进行平滑处理，故实际生产中对复杂介质条件鲜见成效。Gazdag(1978)提出的频率波数域相移法频率波数域相移法具有方便快捷的优点，但它是基于层内常速假设层内常速假设的，不适应介质速度的横向变化不适应介质速度的横向变化。Stoffa(1990)在速度场分裂思想的基础上，提出了分步傅立叶叠分步傅立叶叠后偏移方法后偏移方法。随后，该方法又推广到叠前情况，对较强的速度对较强的速度横向变化都可适应横向变化都可适应。Ristow(1994)也是基于速度
40、场分裂，在分步傅立叶方法的基础上，增加了对介质速度的二阶以上扰动的校正处理。该算法称为傅立叶有限差分方法傅立叶有限差分方法，被公认对复杂地对复杂地质体具有较好的成像效果质体具有较好的成像效果。Wu.R.S.,Huang.L.J.和Jin S.W.（1992,1996,1998）提出的基于散射理论的广义屏或相屏方法基于散射理论的广义屏或相屏方法也对强变速介质具有非常好的成像效果对强变速介质具有非常好的成像效果，但不管是基于de.Wolf近似还是Born近似或 Rytov近似，要么假设散射场相对于入射场较小，要么假设场的变化较小。显然，这些假设条件就相当于限定速度场变化不能任意复杂。加之多数屏方法
41、在复杂介质条件下，并非绝对可靠，并受其稳定条件限制受其稳定条件限制。下面给出的频率空间域有限差分算法频率空间域有限差分算法首先从计算成本计算成本上考虑，以单程波方程作为波场延拓算子，但同时又对单程波方程进行有理分式逼近（王华忠，1997），使其在垂向附近较大角度范围内能尽量准确地描述地震波的传播特征，即尽量提高方程的提高方程的偏移倾角偏移倾角。由于采用的隐式差分格式无条件稳定无条件稳定的，故该算法对任何频率成分及延拓步长都不受限制。有限差分波场延拓算法的优势在于对速度的横向变化有较强的适应能力对速度的横向变化有较强的适应能力。另外，我们在频率域进行波场延拓，可以使差分方程简单化使差分方程简单化
42、，方便计算方便计算，也便于求取成像值便于求取成像值。还可仅对有限频带范围内的地震信号进行可仅对有限频带范围内的地震信号进行波场延拓和成像波场延拓和成像。接下来首先介绍了频率空间域波场延拓算子及其差分计算公式，然后在广义反射系数估算的基础上，提出了几种相近的成像条件，然后探讨了折射项方程和关于浮动坐标变换复原的合并处理，最后凹陷模型及Marmousi模型试验算例演示了该方法的可行性。二二 频率空间域波场延拓算子频率空间域波场延拓算子 从三维声波方程出发，(2-1)222222222tuv1zuyuxu其中 为介质速度。为了处理横向变速问题，对正向传播的下行波进行如下的浮动坐标变换：(2-2)其中
43、 为参考速度。则新坐标系下的波动方程为：(2-3)z,y,x(vv cdzttzzyyxxc0zutzuc2yuxu2222222 在频率-波数域中，方程（2-3）的形式为：(2-4)其中 为圆频率，为波场的频率波数域形式。将（2-4）式改写为：(2-5)从上面得到下列关系式：(2-6)0zuzuic2u)kk(222y2xuu)kkv(u)ciz(2y2x222zciz它表示坐标变换前后的算子关系。因此，下行波正向延拓方程可表示为：(2-7)为了使该方程适应深度偏移的需要，（2-7）式特做如下变换：(2-8)上式分解为 （2-9a）（2-9b）u)kkv(iu)ciz(212x2x22uvi
44、u)kkv(iuciuvizu212y2x22uviu)kkv(izu212y2x22u)c1v1(izu称（2-9a）式为绕射项方程绕射项方程，它可使绕射波收敛使绕射波收敛，(2-9b)为折射折射项方程项方程，它的作用是校正由于横向变速引起的时差校正由于横向变速引起的时差。用有限差分法求解（2-9a）式，需要把其中的根式进行展开，例如泰勒展开、连分式展开或其他逼近方法。本文采用有理分式逼近。为便于实现，我们采用一阶有理分式形式，为提高方程适应陡倾角地层的能力，特对分式系数采用共轭梯度法进行优化（王华忠，1997）。文中单程波方程具有 方程的形式，但可使高达 的陡倾角地层精确成像。为方便书写，
45、特把 记作 ，用一阶有理分式优化展开(2-9a）式，z,y,xz,y,x4570 （2-10）其中：整理（2-10）式，并舍弃小项 ，得到：（2-11）且uxvbbxvaaviu)viz(22221022221077778.0b,00000.1b99203.0a,02252.1a1010uvb)ab(i00030)yuxu(ibv)yuxu(zavzu222222222244957.0b,42614.0a（2-11）和（2-9b）式联立，得到频率空间域任意变速情况下的下行波场深度外推方程。对上行波的反向外推，采用如下坐标变换：（2-12）同理可得到下列关系式：（2-13）cdzttzzyyxx
46、zciz它表示坐标变换前后的算子关系。因此，上行波深度外推方程可表示为：（2-14）为了使该方程适应深度偏移的需要，（2-14）式特做如下变换：（2-15）上式分解为 （2-16a）（2-16b）u)kkv(iu)ciz(212y2x22uviu)kkv(iuciuvizu212y2x22uviu)kkv(izu212y2x22u)c1v1(izu对（2-16a）按下行波一样的逼近方式，并整理得到如下频率空间域的上行波方程：（2-17）将（2-17）式和（2-16b）式联立，可作为上行波波场深度外推公式。三三 上、下行波波动方程差分公式推导上、下行波波动方程差分公式推导 对下行波方程，（2-1
47、1）式用有限差分法求解。令 0)yuxu(ibv)yuxu(zavzu2222222222nj,inj,injiU)(U),z,y,x(U),z,y,x(U （2-18a）（2-18b）（2-18c）zUUzUnj,i1nj,i2UUUnj,i1nj,i21nj,i222y222x22222222yyxx1yxyx 对应于 ，对应于 ，且算子：22x2xxT22y2yyT)1,2,1(TTyx按上面差分公式，（2-11）式可变为：（2-19）化简可得：（2-20）02UU)yTxT(ivbzUU)yTxT(vaTTInj,i1nj,i2y2xnj,i1nj,i2y2x22yyxxnj,iyy2
48、y1yxx2x1x1nj,iyy2y1yxx2x1xUT)i(T)i(IUT)i(T)i(I其中：（2-21）而 ，为差分格式中用来降低频散，提高差分精度的调整系数，一般取 或稍小一些的数。对下行波方程中的折射项方程（2-9b）,容易得到其差分表达式：（2-22）其中：（2-23）由（2-20）、（2-21）、（2-22）及（2-23）式，可进行下行波的深度外推。2y1222y12x1222x1y2zbv,yavx2zbv,xavyx61)UU(i1i1Unj,i1nj,i1nj,i)cv1(v2z 对于上行波，我们注意到其绕射项方程只与下行波绕射项方程在第三项上符号相反，取 为 即可，则仅影
49、响 前的符号。所以易得到其差分格式为：（2-24）对于折射项，其差分方程写为：（2-25）按如上差分方程，就可进行上行波的深度外推了。bb2nj,iyy2y1yxx2x1x1nj,iyy2y1yxx2x1xUT)i(T)i(IUT)i(T)i(I)UU(i1i1Unj,i1nj,i1nj,i四四 成像条件成像条件 叠前偏移成像过程总体上包括两步两步，第一步是上、上、下行波的深度延拓下行波的深度延拓，即将震源波场在时间的正方向向下延拓（震源激发时刻为零时刻），将震源激发产生的记录波场沿时间的反方向向下延拓；第二步应用成像原理，提取成像值应用成像原理，提取成像值。波场延拓方法上面已有叙述，下面着重
50、推导成像条件。设均匀介质中仅在深度 处有一反射界面。设震源波场某一频率 的波场分量为 ，相应频率的地面记录波场为 ，则 （2-26）其中 和 分别为从震源点到深度点 以及从深度点 到地面接受点的传播算子，它们定量化了全部的传播影响。为反射系数矩阵，它由 处的不均匀性引起，它刻划了该处的弹性反射特征，成像的主要目的就是要得到这些反射系数。mz),z(S0),z(P0),z(S)z,z(W)z(R),z,z(W),z(P00mmm00),z,z(Wm0),z,z(W0mmzmz)z(Rmmz对于多层情况，若忽略多次波效应，第一层的波场可描述为：（2-27）式中 为第一层内引起的振幅和相位的变化。延
