三维叠前深度偏移成像理论与方法-课件3

1、 叠前深度偏移成像理论与方法一. Kirchhoff积分法叠前深度偏移二. 波动方程法叠前深度偏移 进入九十年代，地震勘探逐渐向中深层和复杂构造区域的精细勘探发展。在这种形势下，地震偏移作为一种重要的地震数据处理手段，迫切需要新的发展以解决复杂地质体成像问题。实践证明常规的叠后时间偏移在复杂介质下，很难取得好的成像效果。于是，人们提出使用叠前深度偏移。 常规的叠前深度偏移算法可以分为两类：一类是基于射线追踪的Kirchhoff积分算法，另一类是基于波场延拓的波动方程解法（如有限差分偏移、相移偏移和广义屏偏移等）。由于Kirchhoff积分偏移计算效率非常高，它在石油工业领域得到了广泛的应用。然

2、而，利用程函方程求解初至旅行时的Kirchhoff积分偏移在成像效果上不如一些波动方程偏移方法。这主要是因为通过焦散区(在光学中，相位函数的光滑性或单值性被破坏的点的集合，此时振幅趋于无穷大，出现奇点）的射线由于存在多路径而产生的后至波可能携带重要的能量，这部分信息对深度成像而言是不应当忽略的。尽管基于完全Green函数射线理论（例如在计算中考虑多至走时及相应的振幅）的Kirchhoff偏移的成像质量有所提高。但不幸的是，这类算法非常复杂，计算代价特别大，且对焦散区并非完全有效。 第一节 Kirchhoff积分法叠前深度偏移 Kirchhoff积分法叠前深度偏移已在实际生产中应用了多年，并解决

3、了不少复杂构造的成像问题。Kirchhoff积分法的关键是绕射旅行时的计算，目前常用的计算方法是射线追踪法和有限差分法。射线追踪法计算绕射旅行时可分为常速法和变速法，常速法很简单；变速法将在下面介绍。有限差分绕射旅行时计算是基于费马原理，可在直角坐标系或球坐标系实现。一Kirchhoff积分法叠前深度偏移成像公式 为适应复杂地质构造和岩性成像实现保幅处理。可采用如下的考虑传播效应的Kirchhoff积分法叠前深度偏移成像公式（1-1-1） 式中 是记录面（线）； 分别表示成像点、震源点和接收点； 分别表示震源点到成像点和成像点到接收点的旅行时；A是几何扩散因子（振幅加权因子）； 是记录面 的外

4、法线方向；u是记录波场；R是反射系数。 二二维有限差分绕射旅行时计算方法 在Kirchhoff叠前深度偏移中，绕射走时的计算精度直接影响了成像精度。为此，我们利用Eikonal方程和费马原理在规则网格上进行绕射走时的有限差分计算，无需走时的内插。二维空间的Eikonal方程为 （1-2-2） 其中，(x,z)是空间坐标，s是慢度，（1-2-2）式的两个偏导数可由有限差分近似 （1-2-3a） （1-2-3b） 其中， 分别是A，B1，B2，C1处的旅行时，Bi处的旅行时由下式确定 （1-2-4） 把（1-2-3a）和（1-2-3b）代入（1-2-2）得 为获取待求行或列上第一个点处的旅行时，利

5、用有限差分得出如下计算公式 （1-2-5） （1-2-6） 其中， 分别是内行或内列上旅行时的相对极小值、两侧的旅行时、和待求行或列上的旅行时。 （1-2-5）和（1-2-6）式是基于平面波近似得出的有限差分旅行时计算公式。另外也可基于局部圆弧波前近似求取旅行时计算公式。 考虑到计算效率和精度要求，可采用基于平面波近似的有限差分旅行时计算公式。 图1-1. 有限差分旅行时计算示意图 AB1C1C4C3C2B2B3B4 我们可按图1-1的方阵逐步扩展计算旅行时。其具体的做法是：第一，先由（1-2-4）式计算出“十”形端点的值（假设中心点已知，各点慢度已知，网格间距为h），然后由（1-2-5）式计

6、算拐点的值。第二，第一圈计算完后，任取该圈的一边找出极小点，然后以此点为起点由（1-2-6）式计算出下一条边的起点，再利用（1-2-5）式把整条边上的值计算出来。第三，用同样的方法把四条边上的值都计算出来。第四，依此类推，可以把整个模型上所有网格点的到达时计算出来。 利用有限差分法计算绕射走时的Kirchhoff叠前深度偏移方法的优点是运算速度快，能够处理较复杂横向速度变化的情况。并且能解决射线追踪所产生的问题。 三迎风有限差分三维旅行时计算 Kirchhoff积分法三维叠前深度偏移的核心是复杂介质中的旅行时计算。三维旅行时算法的效率和精度直接决定了成像方法的应用范围和效果。考虑到偏移速度场的

7、复杂性和三维叠前深度偏移的大计算量，必须采用稳健高效的三维旅行时计算方法。目前大部分的Kirchhoff积分法三维叠前深度偏移仅利用初至波旅行时。实质上，波前面有时会数次通过速度空间中的某一点，出现多值旅行时，使得初至波旅行时与携带最大能量的波前的到达时不一致，此时仅用初至波旅行时成像显然是不合适的。然而，利用初至波旅行时的Kirchhoff积分法偏移实现简单、计算效率高，且成像精度能满足大多数情况下的实际要求。因此，可采用稳健高效的迎风有限差分三维旅行时计算方法。 目前常用的三维旅行时计算方法大致分为两类：射线追踪旅行时计算和直接解程函方程的旅行时计算。射线追踪需要从程函方程导出射线方程，然

8、后导出射线坐标、方位、和旅行时方程。依此，可用解初值问题的试射法或解边值问题的两点射线追踪法。当地下存在复杂地质结构时，这两种射线追踪法的应用会遇到如下三个困难：一是入射角的微小变化会导致最终结果的很大变化；二是由于阴影区和焦散区的存在，即使从源点出发的射线角度均匀且密集，得到的旅行时场也会存在很多大小不一的空白区和多值区，此时一般的插值方法是不能解决问题的；三是速度场变化剧烈时，无法得到正确的射线路径。直接解程函方程的三维旅行时计算法能克服上述困难，这种计算法会因不同的解法而适应不同的速度场。鉴于计算方法的稳定性和效率以及对复杂速度场的适应性，可采用直接解程函方程的稳健高效的迎风有限差分三维

9、旅行时计算方法。 迎风差分格式就是考虑到波的传播规律构造出的差分格式。该方法较高的计算效率正是源于迎风差分格式，因为它不需要搜索旅行时的相对极小值点或多次计算同一点的旅行时来保证计算过程满足因果律、计算旅行时满足费马原理。另外，迎风有限差分法适合于向量并行巨型机；改进型（无量刚形式）的有限差分方程明显提高了计算方法的稳定性；下面给出的稳定性条件自适应地决定了用于外推的径向（三维旅行时的计算采用了球坐标系）步骤；计算网格的数目、大小、形状、和边界会随径向半径和速度场的复杂性自适应选取；考虑到原坐标系是笛卡尔坐标系，旅行时计算网格是笛卡尔网格，因此球坐标系下的计算网格的边界截断也是自适应的。 在笛

10、卡尔坐标系下的程函方程为 （1-2-7） 球坐标系下的程函方程为 （1-2-8） 其中t是旅行时，S是慢度。定义 （1-2-9） 考虑到计算方法的稳定性，令 （1-2-10） 基于有限差分近似，最后得出 （1-2-11） 经推导有如下的CFL稳定性条件 （1-2-12） 其中 上述方程在全部 和 上可实现向量化和并行化。在满足（1-1-12）式的条件下，利用上述方程可实现三维旅行时场的计算。 四Marmousi模型叠前深度偏移试算 为了检验该叠前深度偏移算法，可采用美国SEG用于检验叠前深度偏移成像方法性能的Marmousi模型。模型速度场如图1-2所示，维数为497x750, 速度场道间距为

11、12.5m，最大深度为3000m，深度采样间隔为4.0m。该速度场的横向变化非常剧烈，因而很适合用于叠前深度偏移方法试验。该模型的正演模拟炮记录为SEG/Marmousi模型数据。共有240炮，每炮96道接收，为右边放炮方式，最小炮检距为-200m，最大炮检距为-2575m，道间距为25m，道长为750个样点，时间采样率为4ms。Kirchhoff叠前深度偏移结果见图1-3。结果表明：尽管在反射振幅上有些失真，但主要构造特征还是清晰可见。 图1-3. 基于二维有限差分绕射旅行时计算的Marmousi模型Kirchhoff叠前深度偏移结果 图1-2. Marmousi模型速度场 第二节 波动方程

12、法叠前深度偏移 基于波动方程的递归偏移方法相对Kirchhoff积分偏移有着得天独厚的优势。首先，它们是从全方程导出的，不是基于高频近似的渐进解。因此，波动方程偏移方法潜在地比较精确、稳健。其次，当向下延拓方法可用于波场延拓而又不增加计算维数时（例如零偏移距数据），这类方法比Kirchhoff积分法的效率还高。 波场延拓算子是递归偏移方法的关键，它们一般是从单程波方程推导出来的。通常大家较为关心的单程波方程偏移的典型算法是有限差分方法和付立叶方法。前者容易处理速度的横向变化，但其缺点在于存在频散和成像倾角限制。后者不存在频散且对水平层状介质能精确成像，不过，它只适用于水平层状地层。 若采用共轭

13、梯度法优化傍轴近似方程的系数，可给出一种频率空间域的有限差分波场延拓方法；在反射系数估算意义下，可推出叠前深度偏移的成像公式。并且发现，当深度偏移算子中的折射项方程采用时移处理，并且和成像计算中关于坐标变换复原的时移处理合并在一起时，计算可大大简化。 为了利用付立叶方法的优势，许多地球物理学家提出在双域（即频率波数域和频率空间域）进行地震波成像处理，先后提出了分步付立叶偏移、付立叶有限差分偏移等算法。下面详细介绍这两种偏移算法的基本原理，分析波场延拓算子的相对误差，而且还对付立叶有限差分偏移算子给出优化改进的思路。 波的传播和成像问题是有着密切联系的。De Hoop（1988）等人为了研究波在

14、随机介质中的传播问题，基于波的散射理论提出了早期的屏方法。Wu R S.&Huang L.J.等在近十年来发展了广义屏方法，并且用于叠前深度偏移成像。下面对比各种广义屏方法，而且讨论它们之间的相互关系和稳定性问题。 首先在Kirchhoff-Helmholtz积分意义下，从方法原理上把前面提到的几种递归偏移方法统一起来，并阐述它们之间的相互联系。随后借助于脉冲响应和复杂模型偏移的数值计算结果，进一步对比频率空间域有限差分法、付立叶有限差分法、分步傅立叶（相屏）法和广义屏方法各自的特点，这些可为我们解决不同的实际问题选择合适的偏移方法提供依据。最后对波动方程三维叠前深度偏移的思路给出探索性观点。

15、 概述 地震波成像在油气勘探中占据重要位置。它的作用是使反射波或衍射波返回到产生它们的地下位置，从而得到地下地质构造的精确成像。 从二十世纪60年代偏移过程由计算机实现以来，已从常规偏移即叠后时间偏移发展到了目前的叠前深度偏移。偏移方法的研究和应用是受油气勘探的实际需求驱动的，同时它又受到人们对偏移成像的认识程度和计算机处理能力的制约。常规偏移（即叠后时间偏移）在以往的油气勘探过程中起到了重要作用，但随着勘探难度的提高，在构造较为复杂的地区，基于常规偏移的处理方法再也难见成效。究其原因，一方面是由于常规处理是先叠加后偏移，水平叠加过程受层状介质假设制约，在复杂地质构造条件下，这种叠加过程很难实

16、现同相叠加，这样会对波场产生破坏，所以用这种失真了的叠后数据去进行偏移处理难以取得好的成像效果就很自然了。为了克服非同相叠加给后续偏移带来的麻烦，人们提出使用叠前偏移，即先偏移处理使波场归位，再把同一地下点的偏移波场相叠加。另一方面是由于时间偏移是建立在均匀介质或水平层状介质的速度模型的基础上的，当速度存在横向变化，或速度分界面不是水平层状的情况下，常规偏移不能满足斯奈耳定律，因此不能进行正确的反射波的偏移成像。为了解决这个问题，出现了深度偏移。这样，叠前深度偏移就可以弥补常规偏移的不足。但是，卓有成效的叠前深度偏移仍然处在探索中。 大家知道，叠前偏移的概念早在70年代中期就提出来了，但由于叠

17、前记录的信噪比较低，偏移的初始模型又很难选准，加之当时的计算机无法承受叠前偏移较大的计算量，直到90年代叠前偏移才开始尝试应用于油气勘探地震数据的精细处理中。常见的叠前深度偏移算法可以分为两类：第一类是基于绕射扫描叠加原理的Kirchhoff积分算法，另一类是基于波动方程的偏移方法（如有限差分偏移方法、付立叶偏移方法等）。 起初出于对计算机处理能力方面的考虑，人们首先想到的是使用一种快捷的算法来实现叠前偏移过程。由于在扫描叠加偏移原理基础上，基于射线追踪的Kirchhoff积分方法正好具有这方面的优势，它很快成为了叠前偏移方法研究的重点，并且很快就推出了应用软件。 工业应用中Kirchhoff

18、积分法叠前偏移在某些地区有时会取得常规偏移难以得到的成像效果。但是，该类方法本身存在很明显的缺陷。例如射线追踪前需要对速度场进行平滑，在速度分布过于复杂的区域，会出现焦散或阴影，计算出来的旅行时场也就不准确。后来，为了提高旅行时场的精度，又发展出了有限差分法直接求解程函方程的Kirchhoff积分偏移方法。但一般仅能计算初至旅行时，无法处理在复杂速度场中存在的多至走时现象，从而影响Kirchhoff积分偏移在复杂地质体（如盐丘、超覆等）的成像效果。从近十年Kirchhoff积分偏移实际应用可以证明这一点。如果应用完全射线理论的Green函数，在计算时求解所有的到达时和相应的振幅值，可以改善该方

19、法的成像质量，但其计算效率又会大打折扣。由于波动方程偏移方法基本不存在Kirchhoff积分偏移这些困难，因此，近年来人们对波动方程叠前偏移进行了大量的研究。 波动方程叠前深度偏移成像解决的是强横向变速条件下复杂地质体的地震波成像问题。近年来一些专家学者在已有方法的基础上提出几种效果较好的叠前偏移方法。大体上讲，这些方法基本上主要分为两类，其一为有限差分偏移方法，另一类为付立叶偏移方法。两类偏移方法各有特点，它们可以分开使用，也可以联合使用（如所谓的混合偏移）。 波动方程有限差分算法借助于差分计算，把速度、密度等介质参数的影响体现在差分计算的矩阵方程中，因此，这类方法能自动适应速度场的任意变化

20、。 这类方法中比较典型的是波动方程有限差分法逆时偏移。它基于双程波方程的，既能适应速度场任意变化，又不存在倾角限制。但它具有一些难以克服的缺点，如计算效率非常低、对计算机内存要求太高等等，这就迫使人们尽量使用单程波方程进行波场延拓和成像。为了实现深度偏移，通常按Hatton(1981)提出的思路，通过引入参考速度，可以把上、下行波方程分裂为绕射项和折射项。绕射项方程实际上就是时间偏移方程，它可使绕射波收敛；折射项方程的作用是校正由于横向变速引起的时差。上、下行波方程是由双程波方程分解近似得到的，它们存在偏移倾角限制，即在某种传播角度范围内才能较准确地描述地震波的传播过程。于是张关泉（1986）

21、、Holberg,O.(1988)等人提出对单程波方程的系数进行优化，尽量提高低阶方程的成像精度。王华忠（1997）用共轭梯度法优化常规的方程，并提出了一种在时空域进行的有限差分深度偏移算法。 付立叶偏移算法一般借助快速付立叶变换来进行波场延拓计算，最早期的付立叶偏移方法应当是相移偏移。这种方法不存在偏移倾角限制，没有频散，而且计算效率非常高。但它是基于层内常速假设条件的，不能处理横向速度变化地层成像问题。为了改进相移偏移方法，人们提出在相移的基础上增加对横向速度扰动引起的时差的校正处理。这些方法的基本思路是进行速度场分裂，即把复杂的介质速度场分裂为“常速背景+层内变速扰动”，然后针对分裂后的

22、速度场分别进行波场延拓处理。例如，Stoffa (1990) 提出分步付立叶偏移方法，常速度背景对应的波场用相移法进行波场延拓，然后针对变速扰动引起的时差进行时移校正。该方法适用于非剧烈变速情况下的深度偏移成像。 Ristow&(1994)提出的付立叶有限差分深度偏移成像方法是在分步付立叶方法的基础上，加上一个有限差分项对二阶以上速度扰动引入的时差进行校正。该方法在剧烈变速情况下也具有非常好的成像效果。90年代初， Wu R S 与de Hoop等在波动方程Green函数解法基础上，通过一系列近似处理手段，发展成了较实用的广义屏算法。这类算法既可用于研究波（声波或弹性波）的传播问题，又可用于地

23、震波场成像。该类方法的基本思路是，将速度场分解为层内常速背景和层内变速扰动。对背景场相当于解常速的声波方程，可通过相移实现；对变速扰动，可认为这种非均匀性相当于散射源（二次源），入射波场作用于这些散射源上，由此产生散射波场。 一般假设波场延拓层厚度较小，在薄板近似处理下，对散射场计算式的积分核再采用不同的近似方法，如Born近似、De Wolf近似或Rytov近似等等，可以得到不同的散射波场表达式。随后在此基础上，Wu R.S和Huang L.J.等人把该类方法发展成了广义屏、相屏以及局部Born近似的屏方法和局部Rytov近似的屏方法。接着他们针对上述方法的一些不足之处做了进一步的扩展和改进

24、，提出了一系列新的屏算子，如扩展的局部Born近似和Rytov近似的屏方法、拟线性局部Born近似的屏方法，等等。 我们在第一部分简要介绍基于共炮道集的波动方程叠前深度偏移的基本思路，接着从第二部分把优化系数的旁轴近似方程转入频率空间域，用有限差分法进行波场延拓，并在反射系数估算意义下推导出叠前深度偏移的成像条件，接着进行数值试验（脉冲响应计算、凹槽模型叠后深度偏移和Marmousi模型叠前深度偏移测试）。第三部分到第五部分依次叙述分步傅立叶法、付立叶有限差分法和广义屏法偏移的基本原理、实现方法和数值试验结果。第六部分阐述前面几种叠前深度偏移方法的相互联系，并把它们从脉冲响应、偏移成像效果、稳

25、定性和计算效率等方面做比较详细的对比，然后对三维波动方程叠前深度偏移做初步的探讨。 第一部分 基于共炮集的波动方程叠前深度偏移概论 基于共炮集的波动方程叠前深度偏移思路是，首先对每一炮进行单炮偏移成像，然后再把各炮成像结果在对应地下位置上叠加，从而得到整个成像剖面。对于每一炮，标准的波动方程叠前深度偏移可以分为三步：震源波场的正向延拓、炮集记录波场的反向延拓和应用成像条件求取成像值（Clearbout, 1971）。 为了方便叙述基于共炮集的波动方程叠前深度偏移的基本过程，我们引入基于单程波方程的波场传播算子（Berkhout,1987）。以频率域二维波场为例，如果对震源波场 和 炮集记录波场

26、做如下定义： 1、 ：它是炮点 处频谱为 的点源激发产生的震源波场，有 （1-1） 2、 ：它是点 处激发，排列接收到的记录波场，该波场可以写成： （1-2）这里 含有一非零道，即在接收点 处的记录道，它满足： （1-3） 3、 ：它表示在深度 处的正向延拓波场，如果引入表征波场从地面传播到深度 的传播算子 ，则有： （1-4）4、 ：它表示记录波场 在深度 的反向延拓波场： （1-5）其中 为记录波场的反向传播算子。因为波场传播算子 描述上行波从深度 到地面的传播过程，故 描述了（向上传播的）记录波场从地面到深度的反向延拓过程（Berkhout,1987）。 和 分别称为下行波和上行波的深度

27、外推算子。实际计算过程中，逐层实现上、下行波的波场延拓和求取成像值。层 内波场延拓如图1所示。 图1 叠前深度偏移波场延拓示意图 注意到所有波场空间上离散分布在采样间隔为 的地震道上，故空间 函数可以用长度为 、振幅为 的加窗函数表示，而且以上公式中的积分可以用离散求和替代。 用传统的偏移公式，我们可以得到对点 处震源在 处的单一记录道的叠前深度偏移结果具有如下形式： （1-6）式中 表示复共轭， 表示取复数的实部。这种频率域的成像公式相当于时间域震源波场正向延拓值与记录波场反向延拓值的互相关。由（1-6）式，可以得到共炮集数据的叠前深度偏移成像公式： （1-7）如果采用的是有限差分算法进行波

28、场延拓计算，则上式中关于 的求和并非是显式的，它一般在对整个单炮记录波场的反向延拓中自动实现。因而，共炮集记录叠前深度偏移公式可表述为： （1-8）这里 为炮 整个记录波场的反向延拓波场值。 从计算角度而言，成像过程是很简单的步骤，波场外推算子的数学形式和计算实现才是地震波偏移成像的核心。目前所有基于波动方程的波场延拓算子不外乎有：波动方程有限差分波场延拓算子和付立叶波场延拓算子。前一类算子既可以在时间空间域又可以在频率空间域用有限差分方法实现波场延拓计算，只是在波动方程在频率空间域的形式更简单，差分计算和成像更方便。第二部分讲述的频率空间域有限差分法叠前深度偏移成像即是这类方法。后一类算子中

29、最熟悉也是最简单的要数相移算子，它在频率波数域计算实现。然而，当速度横向变化时，关于空间坐标的付立叶变换不再成立，这迫使我们在处理横向变速介质中的波的传播和成像问题时退回到空间域。因此，在解决这些问题时，一般基于速度场分裂，对背景场和扰动场分开处理，在频率波数域和频率空间域（或双域daul-domain）交替进行波场延拓计算。如第三、四、五部分讲述的分步傅立叶方法、付立叶有限差分方法和广义屏方法都是采用这类双域延拓算子。 第二部分 频率空间域有限差分法波动方程叠前深度偏移 这一部分介绍一种基于共炮道集的波动方程叠前深度偏移算子，它可以用来进行上下行波的波场延拓。在基于反射系数估算的成像条件下，

30、可实现叠前深度偏移成像。由于算子为优化系数的旁轴近似方程，它具有方程阶数低且能对陡倾角成像的特征，并采用有限差分法波场延拓，能适应速度的任意横向变化，故该算法可使复杂地质体精确成像。因为一切计算在频率-空间域进行，相对于纯粹的时间空间域有限差分算法有计算效率高、成像方便的优点。下面进行了常速和变速介质的脉冲响应测试，并对凹陷模型爆炸反射记录进行了叠后深度偏移，对Marmousi模型进行了叠前深度偏移处理。试验结果表明该算法是比较优越的，且具有可供进一步挖掘的潜力。 一 概 述 叠前深度偏移是一种对复杂地质构造成像的重要的和有效的工具。已有的方法不是在时间域就是在频率域进行。频率域算法的波场延拓

31、既可以在波数域进行，也可以在波数域和空间域交替进行。一种比较典型的叠前时间域算法是波动方程逆时偏移方法（E.G.Chang&McMechan 1990），它一般通过有限差分法实现波场延拓。另一种时间域算法是基于射线追踪的Kirchhoff积分方法（E.G. Hu&McMechan, 1986）。由于全方程逆时偏移的有限差分法是非常耗时的，故其叠前深度方法在目前生产中是难以承受的。而基于射线追踪的Kirchhoff叠前成像是一种高频近似方法，且在利用射线追踪计算绕射旅行时时，常需对速度场进行平滑处理，故实际生产中对复杂介质条件鲜见成效。Gazdag(1978)提出的频率波数域相移法具有方便快捷的

32、优点，但它是基于层内常速假设的，不适应介质速度的横向变化。Stoffa(1990)在速度场分裂思想的基础上，提出了分步傅立叶叠后偏移方法。随后，该方法又推广到叠前情况，对较强的速度横向变化都可适应。Ristow(1994)也是基于速度场分裂，在分步傅立叶方法的基础上，增加了对介质速度的二阶以上扰动的校正处理。该算法称为傅立叶有限差分方法，被公认对复杂地质体具有较好的成像效果。Wu. R.S., Huang. L.J.和Jin S.W.（1992,1996,1998）提出的基于散射理论的广义屏或相屏方法也对强变速介质具有非常好的成像效果，但不管是基于de. Wolf近似还是Born近似或 Ryt

33、ov近似，要么假设散射场相对于入射场较小，要么假设场的变化较小。显然，这些假设条件就相当于限定速度场变化不能任意复杂。加之多数屏方法在复杂介质条件下，并非绝对可靠，并受其稳定条件限制。 下面给出的频率空间域有限差分算法首先从计算成本上考虑，以单程波方程作为波场延拓算子，但同时又对单程波方程进行有理分式逼近（王华忠，1997），使其在垂向附近较大角度范围内能尽量准确地描述地震波的传播特征，即尽量提高方程的偏移倾角。由于采用的隐式差分格式无条件稳定的，故该算法对任何频率成分及延拓步长都不受限制。有限差分波场延拓算法的优势在于对速度的横向变化有较强的适应能力。另外，我们在频率域进行波场延拓，可以使差

34、分方程简单化，方便计算，也便于求取成像值。还可仅对有限频带范围内的地震信号进行波场延拓和成像。 接下来首先介绍了频率空间域波场延拓算子及其差分计算公式，然后在广义反射系数估算的基础上，提出了几种相近的成像条件，然后探讨了折射项方程和关于浮动坐标变换复原的合并处理，最后凹陷模型及Marmousi模型试验算例演示了该方法的可行性。二 频率空间域波场延拓算子从三维声波方程出发， (2-1) 其中 为介质速度。为了处理横向变速问题，对正向传播的下行波进行如下的浮动坐标变换： (2-2)其中 为参考速度。则新坐标系下的波动方程为： (2-3) 在频率-波数域中，方程（2-3）的形式为： (2-4)其中

35、为圆频率， 为波场的频率波数域形式。将（2-4）式改写为： (2-5)从上面得到下列关系式： (2-6) 它表示坐标变换前后的算子关系。因此，下行波正向延拓方程可表示为： (2-7)为了使该方程适应深度偏移的需要，（2-7）式特做如下变换： (2-8)上式分解为 （2-9a） （2-9b） 称（2-9a）式为绕射项方程，它可使绕射波收敛，(2-9b)为折射项方程，它的作用是校正由于横向变速引起的时差。用有限差分法求解（2-9a）式，需要把其中的根式进行展开，例如泰勒展开、连分式展开或其他逼近方法。本文采用有理分式逼近。为便于实现，我们采用一阶有理分式形式，为提高方程适应陡倾角地层的能力，特对分

36、式系数采用共轭梯度法进行优化（王华忠，1997）。文中单程波方程具有 方程的形式，但可使高达 的陡倾角地层精确成像。为方便书写，特把 记作 ，用一阶有理分式优化展开(2-9a）式， （2-10）其中： 整理（2-10）式，并舍弃小项 ，得到： （2-11）且（2-11）和（2-9b）式联立，得到频率空间域任意变速情况下的下行波场深度外推方程。对上行波的反向外推，采用如下坐标变换： （2-12）同理可得到下列关系式： （2-13） 它表示坐标变换前后的算子关系。因此，上行波深度外推方程可表示为： （2-14）为了使该方程适应深度偏移的需要，（2-14）式特做如下变换： （2-15）上式分解为 （

37、2-16a） （2-16b） 对（2-16a）按下行波一样的逼近方式，并整理得到如下频率空间域的上行波方程： （2-17） 将（2-17）式和（2-16b）式联立，可作为上行波波场深度外推公式。 三 上、下行波波动方程差分公式推导对下行波方程，（2-11）式用有限差分法求解。令 （2-18a） （2-18b） （2-18c） 对应于 ， 对应于 ，且算子：按上面差分公式，（2-11）式可变为： （2-19）化简可得： （2-20） 其中： （2-21）而 ，为差分格式中用来降低频散，提高差分精度的调整系数，一般取 或稍小一些的数。对下行波方程中的折射项方程（2-9b）,容易得到其差分表达式：

38、（2-22）其中： （2-23）由（2-20）、（2-21）、（2-22）及（2-23）式，可进行下行波的深度外推。 对于上行波，我们注意到其绕射项方程只与下行波绕射项方程在第三项上符号相反，取 为 即可，则仅影响 前的符号。所以易得到其差分格式为： （2-24）对于折射项，其差分方程写为： （2-25）按如上差分方程，就可进行上行波的深度外推了。 四 成像条件 叠前偏移成像过程总体上包括两步，第一步是上、下行波的深度延拓，即将震源波场在时间的正方向向下延拓（震源激发时刻为零时刻），将震源激发产生的记录波场沿时间的反方向向下延拓；第二步应用成像原理，提取成像值。波场延拓方法上面已有叙述，下面着

39、重推导成像条件。 设均匀介质中仅在深度 处有一反射界面。设震源波场某一频率 的波场分量为 ，相应频率的地面记录波场为 ，则 （2-26）其中 和 分别为从震源点到深度点 以及从深度点 到地面接受点的传播算子，它们定量化了全部的传播影响。 为反射系数矩阵，它由 处的不均匀性引起，它刻划了该处的弹性反射特征，成像的主要目的就是要得到这些反射系数。 对于多层情况，若忽略多次波效应，第一层的波场可描述为： （2-27）式中 为第一层内引起的振幅和相位的变化。延拓到深度 后，得到： （2-28）因此， （2-29）如果不计透射损失， 可定义为： （2-30）它既含有深度 处的反射特性信息 ，又含有由其他

40、深度引起的振幅和相位的畸变 。De Bruin(1992)指出，这种畸变在拉冬变换域内是随机分布的，即与频率 有关，而反射系数矩阵则是与频率 无关的，为了得到反射系数的好的估计值，我们按所有频率求和， 会相干加强，而 则会干涉抵消。于是， （2-31）其中N为频率成分个数。（2-31）式可理解为叠前波场和震源波场反褶积后作线性拉冬变换，然后取原点处的数据。由（2-29）式和（2-31）式得： （2-32）或者写成： （2-33）其中 为 的共轭复数，（2-33）式即为共炮集叠前偏移的成像条件，被称作“相关成像条件”（Clearbout, 1971）。通常也可用（2-33）式的简化形式计算成像值

41、，如下式： （2-34）如果地下介质的非均匀性非常强，速度变化非常剧烈，延拓过程中震源波场会存在严重畸变，这时，用（2-33）式成像比用（2-34）式要准确一些。但为了保证成像计算稳定，我们把（2-33）式改造成下面形式： （2-35）这里， 为稳定常数。利用（2-33）式、（2-34）式或（2-35）式，可以对频率域延拓后的上、下行波波场相关成像。 如果采用了浮动坐标变换，成像条件的形式则相应有所变化。对上行波，浮动坐标变换为 ，延拓后的波场在浮动坐标系和原坐标系下是相同的，即： （2-36）上式转入频率域，可写成： （2-37）进一步变换得到： （2-38） 由（2-38）式可得到： （2

42、-39）同理，对下行波，采用了浮动坐标变换 ，则 （2-40）其中 （2-41）上面（2-39）式、（2-40）式及（2-41）式表示相关成像之前，需要对浮动坐标系下延拓后的上、下行波 作相应的时移处理。 对于叠后深度偏移，仅对爆炸反射记录波场按上行波算子进行延拓和成像，按原坐标系下零时刻成像条件，其成像公式可写为： （2-42）当然，上行波延拓是在浮动坐标系下进行的，因此成像前也要进行坐标复原处理。 五. 折射项及浮动坐标变换复原处理的合并 解决横向变速时引入了浮动坐标变换，波场延拓是在浮动坐标系下进行的，在成像前，要做浮动坐标复原处理。而且我们注意到对折射项方程，既可用有限差分法计算，也可

43、以用时移（相移）实现。如果我们用时移（相移）方式处理折射项方程，并与成像前关于浮动坐标复原的时移（相移）处理放在一起，由于两者均是成像前要进行的计算过程，因此若把这两步时移（相移）合为一步，我们会惊奇地发现，总的时移（相移）计算式大为简化。下面对这一发现做具体说明。 对下行波折射项方程，可变形为： （2-42）则有： （2-43）其中 为浮动坐标系下有限差分延深度延拓一步后得到的波场。如果回到原坐标系下，还要针对坐标变换按（2-40）式做时移处理，得到用于相关成像的原坐标下的波场为： （2-44） 上式意思是说，浮动坐标系下每层有限差分处理后，完全可以把折射项时移处理和针对坐标变换的时移处理合

44、在一起，得到各层总的时移公式。 （2-45） 对上行波，由其折射方程（2-16b）按（2-39）式处理，同样得到折射项方程的时移和针对坐标变换的时移的总和： （2-46） 单纯从数学上看，用时移处理来计算折射项方程，并把这种时移处理与针对浮动坐标变换的时移处理合在一起，不再需要求参考速度，这可以使计算大大简化。若从物理的角度看，这又意味着什么呢？下面以下行波为例，加以解释。下行波折射项方程为： （2-47）上式关于时间作付立叶变换，可得到： （2-48） 上式说明，沿着轨迹： （2-49）的波场值是不变的，保持一个常数。或者说，对于每个积分步长 来说，方程（2-48）把数据移动一个时间，这个时

45、间值为： （2-50）然后看对于浮动坐标所做的处理，由公式（2-2）： （2-51）得到延拓步内针对坐标变换的时移： （2-52）（2-50）式及（2-52）式相加，得到两步时移的总和： （2-53）若将上式再还原到微分形式，得到： （2-54）作用于差分处理后的波场 （2-55） 则有： （2-56）这完全和（2-45）式是吻合的。 综上论证，上面叙述的数学上的简化是有物理意义的。由于差分公式及总的时移公式中不再含有参考速度，故计算中也就不用求取它。折射校正时移和坐标复原时移合在一起处理，比二者分开处理要好得多。因此，我们用时移代替差分处理折射项方程，正好可以利用这一特征。 六 单炮叠前深度

46、偏移流程 基于共炮道集的叠前深度偏移是对每一炮分别成像，然后把所有炮的成像值在相应空间位置叠加，最后得到整个地下的成像剖面。对某一炮，在每一步深度延拓过程中，先分别对震源波场和当前炮记录波场按各自的延拓公式计算，然后依据成像公式求取成像值。接着以延拓后的波场作为下一层延拓的初值，进行同样的延拓和成像计算。下面这个流程图（图2-1）直观地反映了频率空间域有限差分法单炮叠前深度偏移成像的过程。 关于t做FFT变换 该炮叠前记录关于t做FFT变换 当前炮震源模拟记录按叠前成像条件，由延拓后的频率域上下波场相关求和 ，得到深度的成像值层 内下行波有限差分处理层 内下行波总时移处理层 内上行波有限差分处

47、理层 内上行波总时移处理图2-1 频率空间域有限差分法单炮叠前深度偏移流程 七 数值试例 下面从脉冲响应、凹陷模型叠后深度偏移处理和Marmousi模型的叠前成像处理来验证频率空间域有限差分算子的性能。1. 脉冲响应测试 首先，对比该偏移算子在常速度介质和不同速度扰动程度的介质（取不同的值）中脉冲响应曲线。 脉冲放置在 ， 处。以下各图中虚线表示理想的脉冲响应曲线，为一半圆，实线为实际介质中的脉冲响应曲线。下面所有情况的介质速度都为 ，参考速度 取不同值，可以反映不同的速度扰动程度。它们的脉冲响应曲线可以反映算子适应横向速度变化的能力。如图 2-2 （a）示，当取 时，参考速度与实际速度相同，

48、即为均匀常速度介质，此时的脉冲响应曲线与理论曲线在大约 以内的传播角度重合得较好，然后随着传播角度增大，实际响应曲线开始内收，且频散也随之增大。当方程系数为优化值时，常速脉冲响应如图2-2（b）所示。容易发现该脉冲曲线与理论曲线的重合程度比优化处理前要好。这时最大偏移倾角大约为 。 如果取 即 时，速度的扰动程度非常强，达到 ，这时方程系数优化前的脉冲响应曲线如图2-2（c）所示，它和相应的常速度介质中的响应曲线（图（a）基本一样。而优化系数方程的变速脉冲响应曲线如图2-2（d）所示，它和图（b）几乎不存在差异，但与图（c）相比，其偏移倾角明显提高了。由此我们可以得出如下两条结论：一是优化旁轴

49、近似方程的系数，可以提高算子的偏移倾角；二是变速脉冲响应与常速脉冲响应曲线一样，这表明频率空间域的有限差分算子能自动适应速度场任意的横向变化。我们采用的是隐式差分格式，是无条件稳定的。但正如大家所知，该方法与常规单程波方程有限差分法偏移一样，存在角度限制和频散影响。 ( a ) XWFD常速脉冲响应 （ b ）优化系数的XWFD常速脉冲响应 图2-2 频率空间域有限差分算子常速介质中的脉冲响应（ ） ( c ) XWFD变速脉冲响应 ( d ) 优化系数的XWFD变速脉冲响应 图2-3 频率空间域有限差分算子在剧烈变速介质中的脉冲响应 （ ） 2. 凹陷模型叠后深度偏移试验 如果我们仅用上行波

50、算子进行爆炸反射记录波场的深度外推，并用叠后偏移成像条件求取成像值，则可实现叠后自激自收数据的偏移成像。我们设计了一种横向速度变化剧烈的凹陷模型，用以检验频率空间域有限差分偏移算子对速度场的适应能力。 我们采用的爆炸反射记录的速度模型如图 2-3 所示，该模型为三层，速度分别为： 、 和 。第一个速度界面中间呈凹槽型。此处速度横向变化程度非常剧烈，达到 。该速度场的维数为250 x100，道间距为15.0m, 深度采样间隔为12.0m。采用波动方程有限差分法正演模拟，得到如图 2-4 所示的合成的爆炸反射记录。共250道，每道时间样点为1000，时间采样率为1ms。从图中可见，由于第二层中速度

51、场存在剧烈的横向变化，故自激自收剖面中对应第二层界面的同相轴产生明显的畸变。图 2-5 为叠后深度偏移结果。断点绕射波收敛效果非常好，断面整齐清晰，且各层界面形态和位置与速度模型完全吻合，达到了预期的偏移效果。 图2-4 爆炸反射模型速度场示意图 图2-5 合成的爆炸反射记录 图2-6 频率空间域有限差分法叠后深度偏移结果 3. Marmousi模型叠前深度偏移试验 为了检验该叠前深度偏移算法，我们采用了美国SEG用于检验叠前深度偏移成像方法性能的Marmousi模型。模型速度场如图2-6 所示，维数为497x750, 速度场道间距为12.5m，最大深度为3000.0m，深度采样间隔为4.0m

52、。该速度场的横向变化非常剧烈，因而用于叠前深度偏移方法试验非常实用。该模型的正演模拟炮记录由SEG/Marmousi模型数据光盘提供。共有240炮，每炮96道接收，为右边放炮方式，最小炮检距为200m，最大炮检距为2575m，道间距为25m，道长为750，时间采样率为4ms。 我们选取第120炮进行单炮叠前深度偏移成像，为了考虑到偏移孔径问题，特依据观测方式特征，左右两边分别加零道数5道和75道，图 2-7 显示了加零道后的第120炮记录波场，它作为该炮延拓成像的输入。当然，加零道只是一种较简单实用的处理，更合理的办法是按炮检互易的原理，将原单边放炮记录抽成双边放炮记录，然后再延拓成像。波场延

53、拓之前，对炮记录做了频率带通滤波处理，滤波窗为3-8-60-70Hz。第120炮的叠前深度偏移成像剖面如图 2-8 所示。可见单炮偏移能够把该炮覆盖区域的地质层位的像大致成出来，只是一些同相轴能量较弱。所以，多炮偏移结果的叠加处理可以实现有效信号的同相加强，提高信噪比，从而得到整个地下的成像剖面。 多炮偏移剖面如图2-9 和图 2-10 所示，前者滤波窗为0-8-100- 120Hz，后者滤波窗为0-6-60-70Hz。总体而言，该方法对Marmousi模型的叠前深度偏移成像效果很好。浅层的三套断层在成像剖面上各自的形态得到了清晰的展示，断层面整齐明晰。中部的背斜构造也准确地反映出来了。大家最

54、感兴趣的是深层2500m附近嵌在底部中间背斜地层中的小低速体，在成像剖面中具有清晰的形态和边界。从图中每些部位可以发现，采用后者窄频带带通滤波后的波场，其成像结果信噪比要高一些，反映地质层位的同相轴连续性好且能量较高。由此说明，在成像前适当的滤波处理是必要的。而且这样频率加窗，也可以减少用于延拓和成像频率成分，计算量大大降低。 不过，从偏移结果看来，频率空间域有限差分偏移结果总的信噪比还不算很高。这主要受选用的方程和差分格式的频散影响。因此，如果我们进一步提高方程的精度，采用特殊手段来减少频散，无疑会进一步改善频率空间域有限差分法叠前深度偏移成像的效果。这些正是我们今后的改进方向。 图2-7

55、Marmousi模型示意图图2-8 Marmousi模型第120炮记录 图2-9 Marmousi模型第120炮 单炮XWFD偏移剖面 图2-10 Marmousi模型XWFD偏移剖面（0-8100-120Hz） 图2-11（0-6-60-70Hz） xwfd偏移结果 xwfd偏移剖面低通滤波后结果 八 结论与讨论 从频率空间域有限差分偏移方法的基本原理、脉冲响应和模型数据试验结果均表明：波动方程有限差分法偏移能自动适应速度场纵横向的任意变化。这一点是当前所有付立叶偏移方法，如分步傅立叶法、付立叶有限差分法和广义屏方法所不具备的。傍轴近似方程的系数做优化处理，可使较低阶数的单程波方程的最大偏移

56、倾角得到提高。从各种脉冲响应曲线看出，文中提出的偏移算子在约 以内传播角度上具有较高的成像精度。速度横向变化剧烈的凹陷模型偏移取得了非常好的成像效果，Marmousi模型偏移结果表明该算法无论在构造成像，还是在振幅保持上都非常可靠。该方法可进一步在方程形式、差分格式等方面加以改进，减少频散，提高成像精度。 第三部分 分步傅立叶法波动方程叠前深度偏移 在相移偏移方法的基础上，把速度场分解为常速背景和变速扰动两部分，对常速背景在频率波数域采用相移处理，对层内的变速扰动，在频率空间域采用时移校正（第二次相移）。该偏移算法数值上通过了脉冲响应测试、叠后深度偏移和复杂模型叠前深度偏移验证，说明它在较复杂

57、地质条件下是一种稳定快速的叠前深度偏移算法。 一 概 述 偏移方法由于波场延拓不同而相互区别。双程波波动方程有限差分法逆时偏移可以适应速度场的纵横向的任意变化，且不存在偏移倾角限制。但从经济可行性上考虑，人们一般采用单程波方程的有限差分法偏移。这种用于波场延拓的单程波方程是舍弃了高阶项的近似方程，方程的阶数、空间采样率以及差分计算是采用显格式还是隐格式，都会直接影响计算的精度和稳定性。另外有限差分计算还存在频散影响。而相移法偏移（Stolt,1978;Gazdag,1978）是一种典型的付立叶偏移方法，它在频率波数域求解微分方程，计算是精确的,绝对稳定的，由于借助于快速付立叶变换，该算法的运行

58、效率非常高。然而，频率波数域的相移处理是基于层内常速度假设的，不能正确处理横向速度有变化的地震波成像问题。Gazdag &Sguazzero (1984)提出用“相移加内插（PSPI）”来克服相移法这一困难。即在每一层选取多个常速度作为参考速度，每个参考速度按相移法求取延拓波场，然后把各个延拓波场依据实际速度与参考速度的关系函数做内插，得到实际的延拓波场值。这种偏移方法同样是绝对稳定的，但其计算量随所取常速个数呈倍数关系增长，且也仅能适应速度场较缓慢的横向变化。 为了利用付立叶偏移的优势，进一步提高偏移方法适应速度横向变化的能力，Stoffa（1990）在相移偏移的基础上，提出一种新的深度偏移

59、方法，即分步傅立叶法。该方法基于速度场分裂的思想，把整个速度场视为常速背景和变速扰动的叠加。在逐层波场延拓时，针对常速背景采用相移处理，即在频率波数域实现，针对层内的变速扰动，在频率空间域采用时移校正。该方法继承了相移法的优点，同也能适应速度场的中等程度的横向变化。且与相移法深度偏移比较，每层在计算上仅多出一次反傅立叶变换和一次时移校正，比“相移加内插”法要节省得多。 下面依次从方法原理、相对误差分析、实现流程和数值试例等方面对分步傅立叶叠前深度偏移方法加以介绍，最后得出相应的结论。 二 分步傅立叶偏移方法的基本原理恒密度介质中的压缩波的传播特征可用如下方程描述： （3-1）其中， 代表压力值

60、， 是介质速度。将（3-1）式变换到频率域，得 （3-2）其中， 为圆频率， 为波场的频率域形式： （3-3）设 为介质慢度，若将慢度场场分解为两部分： （3-4）其中 为背景慢度场分量，它在层内是一个常数。 为层内扰动慢度分量。定义 为参考慢度。将（3-4）式代入（3-2）式得： （3-5） 其中： （3-6） 通过引进一个源项 ，（3-2）式的齐次方程就转换成了（3-5）式的非齐次方程。依据地震波场的叠加原理，方程（3-5）的解可以表示成： （3-7）其中前者背景慢度引起的波场，它为整个波场的主值部分，后者为波场的扰动项。由于 为（3-5）式所对应的齐次方程的解，故满足： （3-8） 由相