前 言 一
波动是物质运动的重要表现形式.由于物质特征、结构的不同，波在其内部传播具有不同的特征，包含的信息也不尽相同.波在连续介质中的传播问题一直受到极大的重视，波传问题的研究也具有十分重要的工程应用意义，其中一个重要的应用领域是在地球物理方面.地震波勘探（简称为地震勘探）方法是石油工业界使用最为广泛的勘探方法，这种方法利用人工震源激发地震波，利用地震仪把地震波传播的情况记录下来，然后进行处理和分析.地震勘探方法依据地球岩石的弹性假设以及波的传播特性进行勘探和相应的数据处理.弹性波在介质中传播时，其路径、振幅、相位以及波形都会随着介质的物理性质和几何性质的变化而发生变化，地震勘探正是利用波的这些变化规律，根据接收到的资料来推断地下介质的性质，从而达到油气勘探的目的.
由于历史的原因，地下介质的分布极为复杂，在几何结构上，有断层、间断点，在物性组成上有岩石、流体、裂隙.正是由于地质条件的复杂性，要求建立更为复杂、尽可能真实的模型，在此模型基础上建立波动问题的控制方程(波动方程).随着地震勘探进展的不断深入，对勘探精度的要求也越来越高，勘探的目标也逐步从构造性油气藏转向裂缝性油气藏.裂缝性油气藏是指油气的储集空间和渗滤通道主要为裂缝及其连通的孔隙、溶洞的油气储集层，含油气的地层实际上是具有固体状态和流体状态的双相介质.这种介质中的弹性波速度、振幅衰减均受到孔隙介质参数的影响.目前对裂缝性储层中波传播的动力学特征尚处于积极研究阶段，而利用地震波勘探手段寻找裂缝性油气藏是当今勘探地球物理界世界性的热点问题.从提取地层岩性信息的角度看，建立和发展更为符合实际的双相介质波动理论和求解方法，充分阐明并解释各种物理现象以及波传特征具有重要的理论和实际意义.
弹性波动方程有多种求解方法，20世纪50年代初Thomson (1950)和Haskell (1953)提出了传播矩阵法；1968年Alterman和Karal (1968) 提出了数值求解弹性波动方程的有限差分方法；1970年Aki和Larner提出了谱分析的Aki-Larner方法；1975年Smith提出了有限元法；1977年Cerveny对传统的射线方法进行了进一步的研究；1978年Hong和Helmberger提出了横向非均匀介质中的glorified optics方法；1983年Lee和Largstort 提出了处理三维问题的principal curvatures方法；1989年Michel Bouchon等人提出了“边界积分方程-离散波数法”，用来研究波在具有非均匀界面的多层介质中的传播；1990年P.Berg等人提出了谱方法，用来模拟波在弹性介质中的传播.传播矩阵方法和谱方法具有计算速度快、占用内存少、计算精度高等优点，但只适用于简单地层模型（如水平层状）；基于变分和积分原理的有限单元法和边界元法，优点在于其稳定性、收敛性以及边界适应性，但是有限单元法在进行大规模的线性方程组求解，特别是用于波动响应计算的情况下，计算量和内存占用量都非常大，使用不是很广泛；基于惠更斯原理的射线追踪方法计算速度很快，但是这种旅行时（travel time）正演方法缺少振幅信息.有限差分方法是求解双曲型偏微分方程的最常用数值方法之一，通过网格差分近似波动方程微分算子，能够得到完整的弹性波场信息，而且具有编程简单、运算速度快等优点，可以用来分析处理各种复杂地质构造中的波动问题.但是，这种方法也面临计算精度低、复杂界面处理难度大等问题，局部物理和几何参数的变化要求加密整个模型网格，导致计算量的大大增加.算法耗时大是目前几乎所有波动方程模拟方法的特点，这一缺点对一次模拟是可以忍受的，但是大规模反演问题则无法实现.这也是目前采用波动方程作为模型的反演算法仍未走向实际的根本原因.发展一种网格划分灵活、对复杂模型适应性强而且计算精度高的快速不规则网格有限差分方法，已经成为力学、地球物理学等领域国内外学者研究的热点.
二
人类利用的主要矿物与能源资源大多数都埋藏在具有复杂几何构造和复杂物理属性的地球岩石圈中.在工程应用中，人们往往把含油气地层中复杂地质构造和裂缝系统简化为几种简单的模型（如水平层状、垂直裂缝系统等），把多孔介质当作无孔隙固体处理或进行经验修正，这就使得实际资料得不到充分、合理的解释.因此，建立和发展更符合地质实际的非均匀双相各向异性介质波动理论和快速求解方法已是当今地震勘探领域面临的课题.
1951年Gassmann提出了一种流体对弹性波在孔隙介质中传播影响的等效介质理论，把流体饱和双相介质作为某种等效单相介质处理，等效介质的弹性特征由骨架弹性特征和流体弹性特征求取.1956年Biot发表了《流体饱和多孔隙固体的弹性波传播理论》等两篇著名文章，从而奠定了双相介质中弹性波传播规律的基础.Biot理论假定：由弹性各向同性固体构成骨架，连接骨架孔隙体积的空间内充满可压缩性流体，流体相对于固体流动遵循Darcy定律，属于Poiseuill型流动，流体与固体相对运动时存在质量耦合效应，弹性波波场必须明显大于构成孔隙结构的颗粒的最大尺度.Biot理论指出，在流体饱和孔隙介质中传播的弹性波将产生快、慢两种纵波，快纵波类似于无孔隙介质中的普通纵波，对应于固体和同相运动的流体；慢纵波对应于固体和异相运动的流体.1980年Plona用实验证明了慢纵波的存在，为Biot理论提供了有力证据.Amos Nur总结了孔隙岩石中波动传播的特点，较详细地阐述了流体饱和度、裂缝密度、孔隙度、孔隙流体压力和围压、裂隙与孔隙空间的几何形态等因素的影响.这些影响因素对地震勘探信息处理和解释具有重要的意义.
在针对裂缝性储层勘探中，早期油气勘探中起到了非常重要作用的各向同性地层假设已远远不能满足需要.1981年Stuart Crampin详细论述了各向异性介质中弹性波传播问题，指出在各向异性介质中存在横波分裂现象，据此可以推断裂缝发育程度和展布方向.近年来的理论研究与实际观测均表明：具有水平对称轴的横向各向同性（即横观各向同性HTI介质）是描述垂直分布、币状、定向裂缝性地层的最为简单有效的模型，因而在裂缝性油气藏勘探过程中广泛使用.HTI介质可以由5个独立弹性常数来描述，1986年Thomson将5个参数用纵、横波速度和3个无量纲参数 (ε,γ和δ) 描述.在任意垂直面内，这5个弹性常数还可以近似用拉梅常数和两个无量纲参数描述，利用多相介质本构关系，这两个无量纲参数可以指示勘探最关心的参数：裂缝密度和裂缝内含物类型.因此，各向异性介质中弹性波传播特征研究变得更为复杂，而且更具有实际应用价值.
有限差分方法是常用的一种数值解法，它是在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值.该方法广泛应用于数学、力学、地球物理、电磁学等领域，是求解波动方程问题的一种重要方法.根据求解波动方程的不同形式，有限差分方法可以分为二阶差分法和一阶差分法，根据网格特征，可以分为规则网格差分法和不规则网格有限差分法.
地震波正演模拟中有限差分方法的较早研究见于1968年Alterman和Karal的研究，他们给出了均匀各向同性介质二维二阶弹性波动方程有限差分方法，之后Boore(1972)和Alford (1974)等人先后在有限差分方面做了许多工作.Kelly 于1976年提出了非均匀各向同性介质二维二阶弹性波动方程有限差分方法.这些学者采用的都是以位移作为未知变量的波动方程，包含位移对时间和空间的二阶偏导数.直接对位移表示的二阶波动方程进行时间和空间的离散，要涉及三个时间步和空间步的物理量，公式烦琐且计算量大，而且对泊松比变化比较大的模型稳定性差.
1980年Aki和Rechard提出了一维速度应力有限差分法，并对不同的差分格式进行了稳定性分析.1984年Virieux 提出了模拟SH波在非均匀介质中传播的速度-应力有限差分法，1986年Virieux又进一步提出了P-SV波在二维非均匀介质中传播的速度-应力有限差分法.该方法将弹性波动方程用速度和应力来表示，对波动方程作降阶处理，波动方程中位移对时间和空间的二阶导数降为速度和应力对于时间和空间的一阶导数，并采用交错网格进行计算，取得了良好效果.这种方法的重要优点是，对任意泊松比变化的模型正演都很稳定，适合于模拟固体、流体交界面附近的弹性波传播问题；数值频散和数值各向异性都很小，而且对泊松比的变化不敏感；采用具有高精度的差分算子，大大降低了单位波长内的网格点数，提高了计算速度.1995年，Dai给出了非均匀各向同性孔隙介质二维一阶弹性波动方程的有限差分方法，对复杂介质弹性波传播进行了进一步的研究.
以上工作全部基于笛卡儿坐标系中的规则网格，在模拟实际地层中的断层、低速层、孔洞、井眼等复杂构造时，用普通网格模拟曲线界面时出现“阶梯状”边界，在地形构造模型中产生虚假绕射波.另外，局部几何、物理参数的变化也会要求加密整个模型网格，导致计算量的大大增加，限制了有限差分方法在大规模数值模拟中的应用.需要解决的关键问题是，在提高模拟复杂介质模型计算精度的同时，保持算法较高的计算速度.因此，出现了不规则网格有限差分方法.
G.H. Shortley首先于1938年研究了Laplace 方程的不规则网格有限差分法，Mozoc在1989年提出了均匀各向同性介质中具有垂直变步长不规则网格的二维二阶弹性波动方程有限差分法，Jastram和Tessmer, Falk等给出了交错网格上的不规则网格差分方法，Tessmer, Hestholm和Ruud用变形的矩形网格模拟曲线边界，1999年Ivo Oprsal研究了非均匀各向同性介质中矩形不规则网格的二维二阶弹性波动方程有限差分法.Pitarka提出了各向同性介质中矩形不规则网格的有限差分方法，Nordstrm导出了曲线坐标下变形网格的高阶差分方法.不规则网格有限差分方法网格划分灵活，能够较好地模拟各种复杂几何界面情况，与相同规模的规则网格有限差分方法相比，计算精度大大提高；在相同精度要求下，这种方法计算速度快，可节省大量内存，成为有限差分方法极具潜力的发展方向.目前不规则网格有限差分方法还面临一些问题，仍限制在矩形不规则网格范围内，无法突破结构化网格的限制；以二阶弹性波动方程不规则网格有限差分方法为主，公式烦琐，对于复杂模型正演的稳定性和精度还未达到理想的水平.建立一种具有高精度差分算子、高适应网格剖分能力和高计算速度的不规则网格有限差分方法是未来的发展方向.
三
本书的目的就是针对具有复杂几何构造和复杂物质属性模型的弹性波传播问题，研究具有高精度、低运算成本、适应性强的有限差分方法.本书首先从均匀各向同性介质中弹性波动方程的基本理论出发，给出波动方程的一般形式及其求解方法.然后，分别讨论了波动方程的交错网格有限差分方法、不规则网格有限差分方法，通过严格公式推导建立不同格式的有限差分方程，给出了震源和边界条件的处理方法；针对各向异性介质、双相孔隙介质等复杂情况逐步展开探讨，对各种差分格式作了稳定性和数值频散分析，导出了稳定性条件.在理论分析基础上，本书还给出各种不同复杂介质模型的数值算例，并在书中提供相关源程序代码，便于读者迅速理解并掌握波动方程有限差分数值方法.在各章结尾还给出了建议进一步阅读的文献，并在全书结尾给出相关术语的索引，方便读者在全书中迅速查找关心的内容.由于作者水平和经验有限，书中难免出现不妥或者谬误之处，敬请读者和专家批评指正，提出宝贵意见.
本书的主要工作得到了国家自然科学基金（项目批准号： 10402015）和中国石油天然气集团公司中青年创新基金（项目号： 05E7010）的支持，同时也得到了清华大学校内基础研究基金的大力支持，在此表示感谢！
