地球物理资料数字处理(第二十五讲 ).doc

第二节 有限差分法波动方程偏移（Finite Difference Wave Equation Migration） 波动方程偏移是建立在射线理论基础上的绕射扫描偏移具有明显的优越性。因此地震波在地下介质的传播遵循波动方程规律，而射线理论只是近似地描述地震波的传播规律。因此，应用波动方程进行偏移，一方面偏移归位效果好，特别是为复杂的地质条件下反射层的偏移归位提供了可能性；另一方面在实现反射层偏移归位的过程中保持了反射波的特征。1979年，JFClaerbout提出了用有限差分法直接解简化的波动方程，以实现地震波场的空间归位。这种方法是对标量波动方程进行了差分，从波场延拓的观点获得地震偏移剖面。因为用了差分网格，可以在垂向和横向改变速度参数，适用于一般的非均匀介质。这种方法是首先对波动方程进行Claerbout坐标变换，将标量波动方程分为上行波和下行波两部分，这样可以避免用标量波动方程做波场延拓所带来的不适定性问题；然后忽略掉深度方向的二次偏导数项，得到变换后的简化波动方程；最后用有限差分法求解波动方程，使反射层偏移归位到反射界面空间真实位置。在这一节首先介绍波动方程偏移的成像原理，然后介绍叠后有限差分法波动方程偏移的原理和差分方程的建立及其求解。一 波动方程偏移的成像原理（Imaging Principle of Wave Equation Migration ）无论那种偏移方法，包括射线理论的和波动理论的偏移方法，抽象出来都不过是把地表地震剖面转换成近似的地质剖面的一个过程，它包括两个步骤：波场延拓和成像。波场的延拓就是已地表地震剖面为边界条件，求取地下各个点的波场值，也称为波场重建。成像就是在得到了地下各个点的波场之后， 从重建波场中检测出反映界面信息的波场值的过程。JFClaerbout成像原理指出：存在于地壳里的反射界面，是处在下行波到达且刚刚出现上行波那一瞬间。在应用这一成像原理时，针对不同的具体问题又有不同的形式。1. 爆炸反射界面成像原理爆炸反射界面成像原理是最常用，最简单的一种成像原理。它把地下反射界面想象成具有爆炸性的物质或爆炸源，爆炸源的形状，位置与反射界面的形状和位置一致。它所产生的波动脉冲波，其强度，极性与界面的反射系数大小和正负一致，因此称地下反射界面为爆炸反射面。假定在时刻，所有的爆炸反射界面同时起爆，反射上行波到达地面各个观测点，波的传播速度为介质速度的一半，即，如果利用波动方程式将地面测得到的波场记录作反时间方向传播，即向下延拓，则时刻的波场值就正确地描述了地下反射界面的位置，自动实现了偏移成像。此成像原理适用于对水平叠加地震剖面的偏移处理，因为水平叠加剖面相当于零炮检距剖面，自炮点发出的下行波到达反射点的路径与自该点反射返回地面的上行波的路径一致，由于做半速代换，所以叠后偏移方法就是利用延拓波场的时刻的的值进行成像。2. 测线下延成像原理测线下延成像原理适用于有炮检距的的记录，常在叠前偏移时使用。叠前偏移的基本思想就是以地表记录为边界条件。在整条测线上交替的把炮点和检波点向下沿拓，直到炮点和检波器相当于埋到地下同一位置为止。由于这是对地下每一点而言，就好像既埋置了震源又埋置了检波器一样，根据J .F.Claerbout 的原始成像原理，显然这时要取时刻的波场值输出。测线下延成像原理也适用于零炮检距记录的偏移成像，但是对零炮检距记录使用爆炸反射界面成像原理更简单。3. 波场传播的时间一致性成像原理 波场传播时间的一致性成像原理是弹性波偏移中的一种成像原理。因为弹性波偏移同时要处理纵波、横波和转换波，因此它应有自己特殊的成像原理形式。在弹性波偏移中，一般只向下延拓检波点而炮点不动，即延拓后的波场好象是炮点仍在地表位置而地下各点都埋置了检波器接收到的波场一样。为了说明这种成像原理，以入射波（P波）反射纵波（PP波）和转换波（PS）为例。如图7-5-a所示，入射纵波传播时间为t反射纵波传播时间为t转换波传播时间为t。只有在反射点Q处才有t=t=t。即这时反射纵波PP和转换波PS在反射点Q上结合在一起。如图7-5（b）是弹性波的成像原理 其中图中图波场传播时间一致性成像原理示意图由于延拓到点处时，如果不在反射界面上，一般来说，只有当在反射界面上时，才有。因此，在弹性波的延拓波场中，必须取延拓的纵波和横波（包括转换波）波场的传播时间与入射波传播时间一致的那一时间的值来构成偏移输出剖面。二 15有限差分法波动方程的导出 (15 Finite Difference Wave Equation Derivation)Claerbout差分偏移是对水平叠加时间剖面，即自激自收地震记录进行的，根据爆炸反射界面成像原理，认为水平叠加剖面是位于地下反射界面上的炮点零时刻产生的脉冲波，沿界面的法向以半速度上行到地面被检波器记录道的上行波。因此，地下任一点、任一时刻的波场也都可以认为是上行波，都服从半速度的上行波方程。偏移就是利用这种上行波方程，将地面记录的上行波作为边界条件，求零时刻的波场值也就是反射界面处的波场值，偏移技术的关键，就是寻求地震波所遵循的上行波方程。地震波在地下介质中的传播如果介质是均匀的，各项同性的和完全弹性的，则波场函数满足二维纵波波动方程 其中为波的传播速度，在讨论中假设它是常数，在实际应用差分偏移方法时它可以随和而变。对于叠加剖面而言，波动方程（7-2-1）描述的是双程波的传播，既有下行波也有上行波。但是，根据爆炸反射界面成像原理，我们对双程时间的叠加剖面做偏移时，要对方程（7-2-1）做一些改造才能适用。首先，要把双程时间改造为单程时间，这可通过做半速代换来解决，这时方程（7-2-1）就变为： （7-2-2）其次，要做上、下行波分离，从方程（7-2-2）中去掉下行波，把它简化成仅适合上行波的形式。为此，可对（7-2-2）做Claerbout上行波坐标变换： 则相应的变换关系为： 因为坐标变换并不改变实际波场状态，因此存在 由以上变换公式可以得到方程在新坐标系中的形式：以后为了书写的方便把仍写成的形式，则上式应为： 因此，对水平迭加剖面做偏移，可以归结为求解如下解问题：其中就是剖面的范围，T为地震记录最大时间。不难看出，此定解问题是不适定的，因为方程中有关于的二阶偏导数，但在野外只能得到波场，而并不能得到，所以这是一个定解条件不足的欠定问题。为了解决这种不适定问题，Claerbout采取了忽略方程中关于z的二阶偏导数项的方法，从而得到了著名的15o近似方程。 （7-2-10）因为我们进行的是时间偏移，因此要把浓度坐标z换成用时间表示的浓度坐标，这要做变换： （7-2-11）这时方程（7-2-10）变为： （7-2-12）这就是实际中广泛使用的15o近似方程形式。在经过（7-2-11）变换之后，方程（7-2-8）变为 （7-2-13）三 差分方程的建立及求解 (Difference Wave Equation Established & Solving )115差分方程的建立(15Difference Wave Equation Established)使用15的近似方程（7-2-12），偏移定解问题（7-2-9）可以写成 （7-2-14）为了求解此定解问题，用差分方法建立适当的差分方程来近似微分方程，仅仅就差分结解法而言，可以用多种不同的差分格式建立差分方程。用差分方程近似地替代偏微分方程，可分为显式方程和隐性方程两大类。通常，显式方程求解比较容易，快速，但精度不高，对大倾角适应性差，而隐性差分方程虽然求解比较困难，比较慢，但精度高，对倾角大的反射偏移效果好。 在这里我们介绍如图7-6所示的12点差分格式建立的差分方程求解，图中点（i，j，k）表示坐标中的点。表示其中，分别表示坐标点的序号，分别表示的坐标增量。 图7-6 12点的差分格式示意图12点差分网格的重点为表示点。其余12点的波场函数均在点用aylor级数展开。在这里不介绍其详细的展开过程，感兴趣的读者可以参考有关的参考文献。把aylor级数展开式中的高于四阶的小量忽略掉，并用12个点的适当的线形组合可以近似地表示方程（7-2-12）的偏导数；= - + （7-2-15）= + + + （7-2-16） 将（7-2-15）和（7-2-16）代入定解问题中的偏微分方程，并令： ，则得到： -=从而得到15有限差分法偏移的差分表达式：=- （7-1-17）根据上述差分方程，并考虑到定解问题中给定的初使条件和边界条件，可以在整条测线上沿时间深度的方向从开始，由到逐层向下延拓，直到所需要的最大时间深度为止。而在由时间深度向下一个时间深度延拓时，则是沿时间轴的相反方向从最大时间出发，由到直至为止，在时间深度上取时的波场函数值为输出值，最后得到偏移后结果如图7-7所示从中可以看到，偏移时对于时间变量是从最大观测时间来做逆时递推的，这也体现出延拓是使波反向传播的实质。另外，尽管我们所使用15的近似波动方程，只适用常速介质，但在有限差分计算中却可以处理速度的从纵向和横向变化，只要在一个空间差分间隔和一个延拓步长为一个常数就可以。215差分方程的求解(15Difference Wave Equation Solving)求解差分方程（7-2-18）有不同解法，可以利用褶积算法进行求解。方程（7-2-18）右边的系数实际上是对的褶积因子，其分子是褶积因子，分母是反褶积因子。为了用递推方法求出反褶积因子对这个系数用 变换，设=其分子为=分母为=分子，分母的 变换分别为：= =是Z的二次多项式，它有两个跟并且。设其中模小于1的根为则 = （7-2-19）上式中的分子是褶积因子的变化，分母是反褶积因子的变换。令 （7-2-20） （7-2-21）则 （7-2-22）由（7-2-20）式中，表示沿轴的算子，在差分方程中，右边第一项的的褶积因子作用于的褶积结果为： = + （7-2-23） + 由（7-2-21）式可知：或 （7-2-24）设 则由（7-2-24）式得到 ：=令上式两边的幂级数系数相等，得到因此，差分方程（7-2-18）中，右边第一项的的褶积因子作用于的褶积结果为：上式中是求取的递推形式，根据和计算的值。 同样，根据（7-2-22）可以知道：或 因此差分方程（7-2-18）中，右边的第一项的的褶积因子作用的结果为： （7-2-27）上式是求取的递推形式，根据和的值计算出S值，其中 然后，从上面计算出的（7-2-18）式第一项值中减去第二值就可以求得差分方程（7-2-18）式的解，得到的值。 上面介绍的是采用12点中心差分格式建立的差分方程及其求解过程，用差分方程近似微分方程求解，一个很重要的问题是要求差分格式是稳定。只有如此，才能通过差分方程求得微分方程的数值解。在偏移中 ，所谓稳定性问题就是误差是否随差分方程的递推求解过程增大的问题。如果误差不增大，或规定增的不超过某种比较慢的速率，那这种差分格式就是稳定的，经过理论证明，本节所采用的差分格式是稳定的。 15有限差分偏移方法被理论和实践证明，存在倾角的局限性，即对比较大的角度（）的倾斜层偏移效果不好，主要表现是：较大倾角的地震反射波不能正确归位，波形特征失真和出现频散。为了解决这个问题，在15偏移的基础上，又发展了多种适应大倾角的波动方程偏移方法。A Berkhour和R。Stolt先后推导出适用更大倾角的反射界面的二阶近似45方程和三阶近似60方程，它们分别为三阶和四阶偏微分方程，用同样的方法可以推导出更高阶的近似方程。但是利用这些近似方程的偏移要受到倾角的限制，并且由于偏微分方程阶数随着倾角的增大而提高，应用有限差分法来直接求解高阶方程的困难的。八十年代初，马在田提出了高阶方程的偏移的分裂法，把高阶方程分裂成一组二阶微分方程，以避免高阶微商差分带来的困难。原则上可以实现任意倾角的差分偏移，具有重要的理论意义。但是在大倾角的偏移方法的计算量要比15有限差分偏移的计算量大的多。另外，还有许多种大倾角偏移方法，如线