CN105044771A - 基于有限差分法的三维tti双相介质地震波场数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其中该方法包括取得固体和流体应力张量、固体和流体应变张量，并转换为本构方程式；根据应力与位移的对应关系，取得几何方程式；根据本构、几何方程式、流体相对于固体的运动及应力与位移的对应关系，取得运动微分方程式；对运动微分方程两边取散度，取得地震波的第一和第二纵波方程式；对第一和第二纵波方程式，令对y的偏导数等于零，对空间偏导数和时间偏导数采用2N阶精度展开式和二阶精度中心差分格式进行差分离散，取得第一和第二差分方程式；对第一与第二差分方程式进行吸收边界条件处理，以取得对应的地震波场数值。通过本发明，以实现物理地震波场的实时传播模拟。

Description

基于有限差分法的三维TTI双相介质地震波场数值模拟方法

技术领域

本发明涉及地震勘探的技术领域，尤其涉及一种用平均入射角道集进行PP波与PS波联合AVO反演方法。

背景技术

石油和天然气是经济建设和人类生活不可缺少的重要能源，油气的供应直接影响着世界各国的经济发展步伐和人们的生活水平，因此，寻找具有工业意义的油气储集层一直是地球物理学家和地质学家的首要目标。我国从1993年开始成为石油进口国，并且进口额逐年增大，如何摆脱国民经济对油气进口的依赖，提高国内油气产量，是我国油气工业所面临的主要任务。

油气产量的提高是以地质勘探查明的油气新储量与剩余油储量的增加为前提的，这就要求地质工作者必须尽快找到具有工业生产价值的新的油气资源或剩余油气资源。石油地震勘探的主要任务就是解决以上两个问题，为石油工业的持续发展服务，几十年来，我国的地球物理工作者进行了大量的卓有成效的工作，在许多油气田的勘探与开发中发挥了重要作用。但以往利用地震勘探资料找油气都属间接找油气范畴，即由物探工作者将野外地震勘探资料的处理结果提交给地质工作者，再由地质工作者结合地质资料与其它资料确定勘探区的油气分布与储量。这种方法的优点是多学科综合利用，信息量较大，缺点是钻探成功率低，勘探周期长、成本高，经济效益差，并且物探资料中许多有效信息未能得到合理应用，造成资源浪费。利用地震资料直接寻找油气可以大大缩短勘探周期，降低勘探成本，提高钻探成功率。自从70年代初墨西哥湾利用地震亮点技术直接找气成功后，许多地球物理学家投身于该领域相关技术研究并在理论上、方法上和实践中取得了一定突破，如AVO技术、速度反演技术、地震属性技术等。上述技术在某些特定地区取得了好的效果，但在大规模推广过程中遇到了困难，主要原因是这些技术都以传统的单相介质理论作为出发点，没有充分考虑油气藏的双相特征，导致了预测或反演结果的较大误差甚至错误结果。采用基于双相介质理论的油气检测技术可以克服上述缺陷，提高油气检测精度。

双相介质理论是上世纪五十年代开始发展起来的一种新的地震波场理论。该理论假设地下介质是由固相和流相(或气相)组成的。固相是指组成地下岩石的骨架颗粒，流相是指充填于岩石孔隙中的流体。传统的地球物理方法常将地层介质加以简化，看成是纯固体(单相介质)介质，这样做在岩石孔隙度很小或孔隙中只含有束缚流体时是成立的，而当岩石孔隙度较大，而且孔隙中含有连续可动流体时，弹性理论简化就会有较大偏差，甚至得不到正确结论。实际上，大多充满流体或气体的介质(双相介质或多相介质)，由于固体介质与流体介质的相互作用，使之与单相介质的物理－力学性质有巨大差异。例如，如果充填流体是理想液体，则在物相分界面上切向应力为零；如果充填的是气体，则骨架与气体分界面上的切向应力和法向应力均应消失。含油气地层实际上是具有固体状态与流体状态(气体状态)的双相(三相)介质。实践发现，致密介质的经典模型不能很好的描述波在油层中的传播过程，需要对它进行完善。双相介质理论充分考虑了介质的结构、流体与气体的特殊性质、局部特性与整体效应的关系，双相介质模型更接近于实际，更能准确地描述实际地层结构和地层性质，更能适应越来越复杂的油气储藏勘探的实际需要。因此，研究基于双相介质理论的油气检测方法已显得非常必要。

总之，双相介质的基本理论研究已趋于成熟，而其实际应用还处于探索阶段，目前，这一理论在工业应用方面的成功实例还不多见，还没有达到成熟阶段，但理论的创新必将带来技术上的革命，随着研究工作的进一步深入，双相介质中的弹性波理论必将在工业生产中得到广泛的应用，并产生巨大的经济效益。

发明内容

本发明的主要目的在于提供一种基于有限差分法的三维TTI双相介质地震波场数值模拟方法，以实现物理地震波场的实时传播模拟。

为解决上述问题，本发明实施例提供一种基于有限差分法的三维TTI双相介质地震波场数值模拟方法，包括：取得地震波的固体应力张量、流体应力张量、固体应变张量和流体应变张量；根据应力与应变的对应关系，将所述固体应力张量、流体应力张量、固体应变张量和流体应变张量转换为所述地震波的本构方程式；根据应力与位移的对应关系，取得所述地震波的几何方程式；根据所述本构方程式、所述几何方程式、流体相对于固体的运动及应力与位移的对应关系，取得所述地震波的运动微分方程式；对运动微分方程两边取散度，取得所述地震波的第一纵波方程式，并令第一纵波方程式中的耗散系数等于零，以取得第二纵波方程式；对第一纵波方程式，令对y的偏导数等于零，对空间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，取得第一差分方程式，其中N为大于1的正整数；对第二纵波方程式，令对y的偏导数等于零，对空间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，取得第二差分方程式；对所述第一差分方程式与所述第二差分方程式进行吸收边界条件处理，以取得对应的地震波场数值。

根据本发明的技术方案，通过矩形体剖分、离散化，在时间上进行高阶近似，在边界条件上使用交错网格的吸收边界条件，实现了固体相与流体相耦合作用下的双相介质数值方程的迭代求解，即实现物理地震波场的实时传播模拟。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1a和图1b分別是采用以上差分格式计算的均匀双相各向同性介质中固相波场和流相波场的快照的示意图；

图2a和图2b分別是双相各向同性介质中固相波场和流相波场的单炮记录的示意图；

图3为边界入射波的示意图；

图4是吸收边界条件中角点处理方法的示意图；

图5a和图5b分别为采用以上边界条件计算的均匀双相各向同性介质中固相波场和流相波场的快照示意图；

图6a和图6b分別为没有采用任何边界条件处理计算得到的均匀双相各向同性介质中的固相波场和流相波场的快照示意图；

图7a和图7b分别是利用上述方法计算的双相各向同性介质中x分量波场的固相波场和流相波场的快照的示意图；

图7c和图7d分别是利用上述方法计算的双相各向同性介质中z分量波场的固相波场和流相波场的快照的示意图；

图8是交错网格示意图；

图9a和图9b分别是利用上述方法计算的双相各向同性介质中x分量波场的固相波场和流相波场的快照的示意图；

图9c和图9d分别是利用上述方法计算的双相各向同性介质中z分量波场的固相波场和流相波场的快照的示意图；

图10a和图10b分别是模型一的固相波场和流相波场的快照的示意图；

图11a和图11b分别是模型二的固相波场和流相波场的快照的示意图；

图12是两层介质的模型的示意图；

图13a和图13b分别是模型三的固相波场和流相波场的快照的示意图；

图14是交错网格的示意图；

图15a和图15b分别是固相波场和流相波场的波场快照的示意图；

图16a和图16b分别是固相波场和流相波场的波场快照的示意图；

图17a和图17b分别是固相波场和流相波场的波场快照的示意图；

图18a和图18b分别是固相波场和流相波场的波场快照的示意图；

图19a和图19b分别是固相波场和流相波场的波场快照的示意图；

图20a和图20b分别是固相波场和流相波场的波场快照的示意图；

图21a和图21b分别是固相波场和流相波场的波场快照的示意图；

图22和图23分别是x分量和z分量地震记录的示意图；

图24和图25分别是采用二阶时间差分精度、十阶空间差分精度的波场快照的示意图；

图26是根据本发明实施例的基于有限差分法的三维TTI双相介质地震波场数值模拟方法的流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，以下结合附图及具体实施例，对本发明作进一步地详细说明。

双相介质理论认为地下介质由固体和流体两部分组成，固体即为岩石的骨架，流体为充填于骨架孔隙和裂缝中的液体或气体。当孔隙和裂缝中含有两种以上的流体时也称为多相介质。地震波在固相与液相或固相与气相组成的岩石中传播时与其在单相介质中的传播规律不同。由于固体介质与流体介质的相互作用，地震波在双相介质中的传播规律变得更复杂。为方便，在研究双相介质时，一般忽略固体与流体之间的物理、化学作用以及它们在高温高压环境下的热弹性效应。

双相介质中的基本方程包括三组：本构方程、几何方程和运动微分方程。

1.1本构方程

双相介质中的应力张量分为两部分：固体应力张量和流体应力张量。固体应力张量指作用在单位容积立方体每一面的固体部分上的应力，如公式(1.1)所示：

$\left[\begin{array}{ccc}{\mathrm{\σ}}_{xx}& {\mathrm{\τ}}_{xy}& {\mathrm{\τ}}_{xz}\\ {\mathrm{\τ}}_{yx}& {\mathrm{\σ}}_{yy}& {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{zx}& {\mathrm{\τ}}_{zy}& {\mathrm{\σ}}_{zz}\end{array}\right]---\left(1.1\right)$

其中，σ_xx、σ_yy、σ_zz为固体相的正应力，τ_xy、τ_xz、τ_yx、τ_yz、τ_zx、τ_zy为固体相的切应力，且τ_xy＝τ_yx，τ_xz＝τ_zx，τ_yz＝τ_zy。

流体应力张量指作用在该单位容积立方体每一面的流体部分上的应力，如公式(1.2)所示：

$\left[\begin{array}{ccc}s& 0& 0\\ 0& s& 0\\ 0& 0& s\end{array}\right]---\left(1.2\right)$

s满足公式(1.3)，如下所示：

s＝-φp(1.3)

其中，φ为每单位截面中流体面积所占百分数，相当于有效孔隙率，p为流体压力，负号表示应力s与流体压力p方向相反。由于流体内无切应力，所以应力张量只含有主对对角分量。

双相介质中的应变张量也分为两部分：固体应变张量和流体应变张量。

作用于单元体每一固体截面的应变张量，如公式(1.4)所示：

$\left[\begin{array}{ccc}{e}_{xx}& e_{xy}& e_{xz}\\ e_{yx}& e_{yy}& e_{yz}\\ e_{zx}& e_{zy}& e_{zz}\end{array}\right]---\left(1.4\right)$

其中，e_xx、e_yy、e_zz表示固相正应变，e_xy、e_xz、e_yx、e_yz、e_zx、e_zy表示固相切应变，且e_xy＝e_yx，e_xz＝e_zx，e_yz＝e_zy。

作用于单元体每一流体截面的应变张量用下列矩阵表示，如公式(1.5)所示：

$\left[\begin{array}{ccc}\mathrm{\ϵ}& 0& 0\\ 0& \mathrm{\ϵ}& 0\\ 0& 0& \mathrm{\ϵ}\end{array}\right]---\left(1.5\right)$

其中，ε表示流相的体应变。

双相线性弹性介质的应力与应变具有线性关系，用广义Hooke定律表示，如公式(1.6)所示：

$\left[\begin{array}{c}{\mathrm{\σ}}_{xx}\\ {\mathrm{\σ}}_{yy}\\ {\mathrm{\σ}}_{zz}\\ {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{zx}\\ {\mathrm{\τ}}_{xy}\\ s\end{array}\right]=\left[\begin{array}{ccccccc}{C}_{1111}& C_{1122}& C_{1133}& C_{1123}& C_{1131}& C_{1112}& Q_1\\ C_{2211}& C_{2222}& C_{2233}& C_{2223}& C_{2231}& C_{2212}& Q_2\\ C_{3311}& C_{3322}& C_{3333}& C_{3323}& C_{3331}& C_{3312}& Q_3\\ C_{2311}& C_{2322}& C_{2333}& C_{2323}& C_{2331}& C_{2312}& Q_4\\ C_{3111}& C_{3122}& C_{3133}& C_{3123}& C_{3131}& C_{3112}& Q_5\\ C_{1211}& C_{1222}& C_{1233}& C_{1223}& C_{1231}& C_{1212}& Q_6\\ Q_1& Q_2& Q_3& Q_4& Q_5& Q_6& R\end{array}\right]\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ 2e_{yz}\\ 2e_{zx}\\ 2e_{xy}\\ \mathrm{\ϵ}\end{array}\right]---\left(1.6\right)$

其中，C_ijkl(i、j、k、l＝1，2，3，4)是固体相的弹性参数，R是流体相的弹性参数，Q＝(Q₁，Q2，Q3，Q4，Q5，Q6)^T是固相和流相耦合关系的弹性参数。

当介质为各向同性介质时，广义Hooke定律变为公式(1.7)，如下所示：

$\left[\begin{array}{c}{\mathrm{\σ}}_{xx}\\ {\mathrm{\σ}}_{yy}\\ {\mathrm{\σ}}_{zz}\\ {\mathrm{\τ}}_{yz}\\ {\mathrm{\τ}}_{zx}\\ {\mathrm{\τ}}_{xy}\\ s\end{array}\right]=\left[\begin{array}{ccccccc}{C}_{1111}& C_{1122}& C_{1133}& 0& 0& 0& Q\\ C_{2211}& C_{2222}& C_{2233}& 0& 0& 0& Q\\ C_{3311}& C_{3322}& C_{3333}& 0& 0& 0& Q\\ 0& 0& 0& C_{2323}& 0& 0& 0\\ 0& 0& 0& 0& C_{3131}& 0& 0\\ 0& 0& 0& 0& 0& C_{1212}& 0\\ Q& Q& Q& 0& 0& 0& R\end{array}\right]\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ 2e_{yz}\\ 2e_{zx}\\ 2e_{xy}\\ \mathrm{\ϵ}\end{array}\right]---\left(1.7\right)$

其中，C₁₁₂₂＝C₂₂₁₁＝C₁₁₃₃＝C₃₃₁₁＝C₂₂₃₃＝C₃₃₂₂＝A，C₁₂₁₂＝C₂₃₂₃＝C₃₁₃₁＝N，C₁₁₁₁＝C₂₂₂₂＝C₃₃₃₃＝A+2N。所以，在双相各向同性介质中只有四个独立的弹性参数。此时应力与应变的关系又可写为公式(1.8)，如下所示：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{xx}=2N{e}_{xx}+A\mathrm{\θ}+Q\mathrm{\ϵ}\\ {\mathrm{\σ}}_{yy}=2N{e}_{yy}+A\mathrm{\θ}+Q\mathrm{\ϵ}\\ {\mathrm{\σ}}_{zz}=2N{e}_{zz}+A\mathrm{\θ}+Q\mathrm{\ϵ}\\ {\mathrm{\τ}}_{xy}=2N{e}_{xy}\\ {\mathrm{\τ}}_{xz}=2N{e}_{xz}\\ {\mathrm{\τ}}_{yz}=2{\mathrm{Ne}}_{yz}\\ s=Q\mathrm{\θ}+R\mathrm{\ϵ}\end{array}\right.---\left(1.8\right)$

式中，θ表示固相体应变， $\mathrm{\θ}=divu=\frac{\∂u_x}{\∂x}+\frac{\∂u_y}{\∂y}+\frac{\∂u_z}{\∂z}=e_{xx}+e_{yy}+e_{zz};$ u为固相位移向量，u_x、u_y、u_z分别为固相位移向量u在x、y和z方向的分量；ε表示流相体应变，U为流相位移向量，U_x、U_y、U_z分别表示流相位移向量U在x、y和z方向的分量；A和N相当于单相各向同性弹性波理论中的拉梅系数，其中N＝μ；R表示使一定体积的流体流入某集合体而又使该集合体保持总体积不变所需施加在流体上的压力的一种量度；Q反映固体与流体体积变化之间的耦合性质。

1.2几何方程

双相介质中的几何方程也即双相介质中的应变与位移的关系，如公式(1.9)所示：

$\left[\begin{array}{c}\begin{array}{ccc}{e}_{xx}=\frac{\∂u_x}{\∂x}& e_{yy}=\frac{\∂u_y}{\∂y}& e_{zz}=\frac{\∂u_z}{\∂z}\\ e_{xy}=\frac{1}{2}\left(\frac{\∂u_x}{\∂y}+\frac{\∂u_y}{\∂x}\right)& e_{xz}=\frac{1}{2}\left(\frac{\∂u_x}{\∂z}+\frac{\∂u_z}{\∂x}\right)& e_{yz}=\frac{1}{2}\left(\frac{\∂u_y}{\∂z}+\frac{\∂u_z}{\∂y}\right)\\ e_{xy}=e_{yx}& e_{xz}=e_{zx}& e_{yz}=e_{zy}\end{array}\\ \begin{array}{cc}e=\frac{\∂u_x}{\∂x}+\frac{\∂u_y}{\∂y}+\frac{\∂u_z}{\∂z}& \mathrm{\ϵ}=\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\end{array}\end{array}\right]---\left(1.9\right)$

1.3运动微分方程

假设组成双相介质的固体骨架是统计各向同性的，孔隙是连通的，主孔隙具有不可渗透的孔隙壁，孔隙中充满各向同性的粘滞流体，流体可压缩，流体与固体之间存在相对位移，并可产生摩擦，流体相对固体之间的流动是Poiseuille型流动。根据达西定律可得到双相介质中的运动微分方程。

在双相介质中，流体相对于固体的运动满足广义达西定律，即如公式(1.10)所示：

$\left[\begin{array}{c}\frac{\∂s}{\∂x}+F_x^f\\ \frac{\∂s}{\∂y}+F_y^f\\ \frac{\∂s}{\∂z}+F_z^f\end{array}\right]=\left[\begin{array}{c}\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_x+{\mathrm{\ρ}}_{22}{U}_x\right)}{\∂t^2}\\ \frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_y+{\mathrm{\ρ}}_{22}{U}_y\right)}{\∂t^2}\\ \frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_z+{\mathrm{\ρ}}_{22}{U}_z\right)}{\∂t^2}\end{array}\right]+\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]\left[\begin{array}{c}\frac{\∂\left(U_x-u_x\right)}{\∂t}\\ \frac{\∂\left(U_y-u_y\right)}{\∂t}\\ \frac{\∂\left(U_z-u_z\right)}{\∂t}\end{array}\right]---\left(1.10\right)$

相应地，固体介质所满足的应力-位移关系如公式(1.10)所示：

$\left[\begin{array}{c}\frac{\∂{\mathrm{\σ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z}+F_x^s\\ \frac{\∂{\mathrm{\τ}}_{yx}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{yy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂z}+F_y^s\\ \frac{\∂{\mathrm{\τ}}_{zx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{zy}}{\∂y}+\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}+F_z^s\end{array}\right]=\left[\begin{array}{c}\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_x+{\mathrm{\ρ}}_{12}{U}_x\right)}{\∂t^2}\\ \frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_y+{\mathrm{\ρ}}_{12}{U}_y\right)}{\∂t^2}\\ \frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_z+{\mathrm{\ρ}}_{12}{U}_z\right)}{\∂t^2}\end{array}\right]-\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]\left[\begin{array}{c}\frac{\∂\left(U_x-u_x\right)}{\∂t}\\ \frac{\∂\left(U_y-u_y\right)}{\∂t}\\ \frac{\∂\left(U_z-u_z\right)}{\∂t}\end{array}\right]---\left(1.11\right)$

其中，和分别表示固相外力的三个分量，和分别表示流相外力的三个分量；[b_ij]_3×3为耗散系数矩阵，b_ij＝b_ji，其值由达西渗透系数k_ij、流体粘滞系数η和孔隙率φ决定，且满足公式(1.12)，如下所示：

$b_{ij}=\frac{{\mathrm{\η\φ}}^2}{k_{ij}}---\left(1.12\right)$

其中，ρ₁₁、ρ₂₂和ρ₁₂是质量密度参数，ρ₁₁表示单元体中固体相对流体运动时固体部分总的等效质量，ρ₂₂表示单元体中流体相对固体运动时流体部分总的等效质量，ρ₁₂表示流体和固体之间的质量耦合系数。ρ₁₁、ρ₂₂和ρ₁₂满足如公式(1.10)所示：

$\left\{\begin{array}{c}{\mathrm{\ρ}}_{11}+2{\mathrm{\ρ}}_{12}+{\mathrm{\ρ}}_{22}=\mathrm{\ρ}\\ \mathrm{\ρ}={\mathrm{\ρ}}_1+{\mathrm{\ρ}}_2\\ {\mathrm{\ρ}}_1={\mathrm{\ρ}}_{11}+{\mathrm{\ρ}}_{12}=(1-\mathrm{\φ}){\mathrm{\ρ}}_s\\ {\mathrm{\ρ}}_2={\mathrm{\ρ}}_{12}+{\mathrm{\ρ}}_{22}={\mathrm{\φ\ρ}}_f\end{array}\right.---\left(1.13\right)$

ρ表示流体和固体组成的集合体的总质量密度，ρ_s表示固体的质量密度，ρ_f表示流体的质量密度，ρ₁表示集合体每单位体积固体的质量，ρ₂表示集合体每单位体积流体的质量。

当忽略流体与固体之间相对运动所产生的耗散力时，便得到保守系中(无耗散情况下)用应力表示的运动微分方程式，如公式(1.14)所示：

$\{\begin{array}{c}\frac{\∂{\mathrm{\σ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z}+F_x^s=\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_x+{\mathrm{\ρ}}_{12}{U}_x\right)}{\∂t^2}\\ \frac{\∂s}{\∂x}+F_x^f=\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_x+{\mathrm{\ρ}}_{22}{U}_x\right)}{\∂t^2}\\ \frac{\∂{\mathrm{\τ}}_{yx}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{yy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂z}+F_y^s=\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_y+{\mathrm{\ρ}}_{12}{U}_y\right)}{\∂t^2}\\ \frac{\∂s}{\∂y}+F_y^f=\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_y+{\mathrm{\ρ}}_{22}{U}_y\right)}{\∂t^2}\\ \frac{\∂{\mathrm{\τ}}_{zx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{zy}}{\∂y}+\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}+F_z^s=\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_z+{\mathrm{\ρ}}_{12}{U}_z\right)}{\∂t^2}\\ \frac{\∂s}{\∂z}+F_z^f=\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_z+{\mathrm{\ρ}}_{22}{U}_z\right)}{\∂t^2}\end{array}---\left(1.14\right)$

不计外力时，公式(1.14)变为公式(1.15)，如下所示：

$\{\begin{array}{c}\frac{\∂{\mathrm{\σ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z}=\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_x+{\mathrm{\ρ}}_{12}{U}_x\right)}{\∂t^2}\\ \frac{\∂s}{\∂x}=\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_x+{\mathrm{\ρ}}_{22}{U}_x\right)}{\∂t^2}\\ \frac{\∂{\mathrm{\τ}}_{yx}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{yy}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂z}=\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_y+{\mathrm{\ρ}}_{12}{U}_y\right)}{\∂t^2}\\ \frac{\∂s}{\∂y}=\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_y+{\mathrm{\ρ}}_{22}{U}_y\right)}{\∂t^2}\\ \frac{\∂{\mathrm{\τ}}_{zx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{zy}}{\∂y}+\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}=\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_z+{\mathrm{\ρ}}_{12}{U}_z\right)}{\∂t^2}\\ \frac{\∂s}{\∂z}=\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_z+{\mathrm{\ρ}}_{22}{U}_z\right)}{\∂t^2}\end{array}---\left(1.15\right)$

将广义Hooke定律(即公式(1.15))和应变-位移关系(即公式(1.9))代入用应力表示的运动微分方程(即(1.10)和(1.11))，并略去外力，可得到用位移表示的运动微分方程，如公式(1.16a)、(1.16b)、(1.16c)、(1.16d)、(1.16e)、(1.16f)所示：

$\begin{array}{c}{c}_{1111}\frac{{\∂}^2{u}_x}{\∂x^2}+c_{\mathrm{1122}}\frac{{\∂}^2{u}_y}{\∂x\∂y}+c_{1133}\frac{{\∂}^2{u}_z}{\∂x\∂z}+c_{1123}\left(\frac{{\∂}^2{u}_y}{\∂x\∂z}+\frac{{\∂}^2{u}_z}{\∂x\∂y}\right)+c_{1131}\left(\frac{{\∂}^2{u}_z}{\∂x^2}+\frac{{\∂}^2{u}_x}{\∂x\∂z}\right)+c_{1112}\left(\frac{{\∂}^2{u}_x}{\∂x\∂y}+\frac{{\∂}^2{u}_y}{\∂x^2}\right)\\ +c_{1211}\frac{{\∂}^2{u}_x}{\∂x\∂y}+c_{\mathrm{1222}}\frac{{\∂}^2{u}_y}{\∂y^2}+c_{1233}\frac{{\∂}^2{u}_z}{\∂y\∂z}+c_{1223}\left(\frac{{\∂}^2{u}_y}{\∂y\∂z}+\frac{{\∂}^2{u}_z}{\∂y^2}\right)+c_{1231}\left(\frac{{\∂}^2{u}_z}{\∂x\∂y}+\frac{{\∂}^2{u}_x}{\∂y\∂z}\right)+c_{1212}\left(\frac{{\∂}^2{u}_x}{\∂y^2}+\frac{{\∂}^2{u}_y}{\∂x\∂y}\right)\\ +c_{3111}\frac{{\∂}^2{u}_x}{\∂x\∂z}+c_{\mathrm{3122}}\frac{{\∂}^2{u}_y}{\∂y\∂z}+c_{3133}\frac{{\∂}^2{u}_z}{\∂z^2}+c_{3123}\left(\frac{{\∂}^2{u}_y}{\∂z^2}+\frac{{\∂}^2{u}_z}{\∂y\∂z}\right)+c_{3131}\left(\frac{{\∂}^2{u}_z}{\∂x\∂z}+\frac{{\∂}^2{u}_x}{\∂z^2}\right)+c_{3112}\left(\frac{{\∂}^2{u}_x}{\∂y\∂z}+\frac{{\∂}^2{u}_y}{\∂x\∂z}\right)\\ +Q_1\frac{\∂}{\∂x}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)+Q_6\frac{\∂}{\∂y}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)+Q_5\frac{\∂}{\∂z}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)\\ =\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_x+{\mathrm{\ρ}}_{12}{U}_x\right)}{\∂t^2}-b_{11}\frac{\∂\left(U_x-u_x\right)}{\∂t}-b_{12}\frac{\∂\left(U_y-u_y\right)}{\∂t}-b_{13}\frac{\∂\left(U_z-u_z\right)}{\∂t}\end{array}---\left(1.16a\right)$

$\begin{array}{c}{c}_{1211}\frac{{\∂}^2{u}_x}{\∂x^2}+c_{\mathrm{1222}}\frac{{\∂}^2{u}_y}{\∂x\∂y}+c_{1233}\frac{{\∂}^2{u}_z}{\∂x\∂z}+c_{1223}\left(\frac{{\∂}^2{u}_y}{\∂x\∂z}+\frac{{\∂}^2{u}_z}{\∂x\∂y}\right)+c_{1231}\left(\frac{{\∂}^2{u}_z}{\∂x^2}+\frac{{\∂}^2{u}_x}{\∂x\∂z}\right)+c_{1212}\left(\frac{{\∂}^2{u}_x}{\∂x\∂y}+\frac{{\∂}^2{u}_y}{\∂x^2}\right)\\ +c_{2211}\frac{{\∂}^2{u}_x}{\∂x\∂y}+c_{\mathrm{2222}}\frac{{\∂}^2{u}_y}{\∂y^2}+c_{2233}\frac{{\∂}^2{u}_z}{\∂y\∂z}+c_{2223}\left(\frac{{\∂}^2{u}_y}{\∂y\∂z}+\frac{{\∂}^2{u}_z}{\∂y^2}\right)+c_{2231}\left(\frac{{\∂}^2{u}_z}{\∂x\∂y}+\frac{{\∂}^2{u}_x}{\∂y\∂z}\right)+c_{2212}\left(\frac{{\∂}^2{u}_x}{\∂y^2}+\frac{{\∂}^2{u}_y}{\∂x\∂y}\right)\\ +c_{2311}\frac{{\∂}^2{u}_x}{\∂x\∂z}+c_{\mathrm{2322}}\frac{{\∂}^2{u}_y}{\∂y\∂z}+c_{2333}\frac{{\∂}^2{u}_z}{\∂z^2}+c_{2323}\left(\frac{{\∂}^2{u}_y}{\∂z^2}+\frac{{\∂}^2{u}_z}{\∂y\∂z}\right)+c_{2331}\left(\frac{{\∂}^2{u}_z}{\∂x\∂z}+\frac{{\∂}^2{u}_x}{\∂z^2}\right)+c_{2312}\left(\frac{{\∂}^2{u}_x}{\∂y\∂z}+\frac{{\∂}^2{u}_y}{\∂x\∂z}\right)\\ +Q_6\frac{\∂}{\∂x}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)+Q_2\frac{\∂}{\∂y}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)+Q_4\frac{\∂}{\∂z}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)\\ =\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_y+{\mathrm{\ρ}}_{12}{U}_y\right)}{\∂t^2}-b_{21}\frac{\∂\left(U_x-u_x\right)}{\∂t}-b_{22}\frac{\∂\left(U_y-u_y\right)}{\∂t}-b_{23}\frac{\∂\left(U_z-u_z\right)}{\∂t}\end{array}---\left(1.16b\right)$

$\begin{array}{c}{c}_{3111}\frac{{\∂}^2{u}_x}{\∂x^2}+c_{\mathrm{3122}}\frac{{\∂}^2{u}_y}{\∂x\∂y}+c_{3133}\frac{{\∂}^2{u}_z}{\∂x\∂z}+c_{3123}\left(\frac{{\∂}^2{u}_y}{\∂x\∂z}+\frac{{\∂}^2{u}_z}{\∂x\∂y}\right)+c_{3131}\left(\frac{{\∂}^2{u}_z}{\∂x^2}+\frac{{\∂}^2{u}_x}{\∂x\∂z}\right)+c_{3112}\left(\frac{{\∂}^2{u}_x}{\∂x\∂y}+\frac{{\∂}^2{u}_y}{\∂x^2}\right)\\ +c_{2311}\frac{{\∂}^2{u}_x}{\∂x\∂y}+c_{\mathrm{2322}}\frac{{\∂}^2{u}_y}{\∂y^2}+c_{2333}\frac{{\∂}^2{u}_z}{\∂y\∂z}+c_{2323}\left(\frac{{\∂}^2{u}_y}{\∂y\∂z}+\frac{{\∂}^2{u}_z}{\∂y^2}\right)+c_{2331}\left(\frac{{\∂}^2{u}_z}{\∂x\∂y}+\frac{{\∂}^2{u}_x}{\∂y\∂z}\right)+c_{2312}\left(\frac{{\∂}^2{u}_x}{\∂y^2}+\frac{{\∂}^2{u}_y}{\∂x\∂y}\right)\\ +c_{3311}\frac{{\∂}^2{u}_x}{\∂x\∂z}+c_{\mathrm{3322}}\frac{{\∂}^2{u}_y}{\∂y\∂z}+c_{3333}\frac{{\∂}^2{u}_z}{\∂z^2}+c_{3323}\left(\frac{{\∂}^2{u}_y}{\∂z^2}+\frac{{\∂}^2{u}_z}{\∂y\∂z}\right)+c_{3331}\left(\frac{{\∂}^2{u}_z}{\∂x\∂z}+\frac{{\∂}^2{u}_x}{\∂z^2}\right)+c_{3312}\left(\frac{{\∂}^2{u}_x}{\∂y\∂z}+\frac{{\∂}^2{u}_y}{\∂x\∂z}\right)\\ +Q_5\frac{\∂}{\∂x}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)+Q_4\frac{\∂}{\∂y}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)+Q_3\frac{\∂}{\∂z}\left(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}+\frac{\∂U_z}{\∂z}\right)\\ =\frac{{\∂}^2\left({\mathrm{\ρ}}_{11}{u}_z+{\mathrm{\ρ}}_{12}{U}_z\right)}{\∂t^2}-b_{31}\frac{\∂\left(U_x-u_x\right)}{\∂t}-b_{32}\frac{\∂\left(U_y-u_y\right)}{\∂t}-b_{33}\frac{\∂\left(U_z-u_z\right)}{\∂t}\end{array}---\left(1.16c\right)$

$\begin{array}{c}{Q}_1\frac{{\∂}^2{u}_x}{\∂x^2}+Q_2\frac{{\∂}^2{u}_y}{\∂x\∂y}+Q_3\frac{{\∂}^2{u}_z}{\∂x\∂z}+Q_4(\frac{{\∂}^2{u}_y}{\∂x\∂z}+\frac{{\∂}^2{u}_z}{\∂x\∂y})+Q_5(\frac{{\∂}^2{u}_z}{\∂x^2}+\frac{{\∂}^2{u}_x}{\∂x\∂z})+Q_6(\frac{{\∂}^2{u}_x}{\∂x\∂y}+\frac{{\∂}^2{u}_y}{\∂x^2})\\ +R\frac{\∂}{\∂x}(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}-+\frac{\∂U_z}{\∂z})\\ =\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_x+{\mathrm{\ρ}}_{22}{U}_x\right)}{\∂t^2}+b_{11}\frac{\∂\left(U_x-u_x\right)}{\∂t}+b_{12}\frac{\∂\left(U_y-u_y\right)}{\∂t}+b_{13}\frac{\∂\left(U_z-u_z\right)}{\∂t}\end{array}---\left(1.16d\right)$

$\begin{array}{c}{Q}_1\frac{{\∂}^2{u}_x}{\∂x\∂y}+Q_2\frac{{\∂}^2{u}_y}{\∂y^2}+Q_3\frac{{\∂}^2{u}_z}{\∂y\∂z}+Q_4(\frac{{\∂}^2{u}_y}{\∂y\∂z}+\frac{{\∂}^2{u}_z}{\∂y^2})+Q_5(\frac{{\∂}^2{u}_z}{\∂x\∂y}+\frac{{\∂}^2{u}_x}{\∂y\∂z})+Q_6(\frac{{\∂}^2{u}_x}{\∂y^2}+\frac{{\∂}^2{u}_y}{\∂x\∂y})\\ +R\frac{\∂}{\∂y}(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}-+\frac{\∂U_z}{\∂z})\\ =\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_y+{\mathrm{\ρ}}_{22}{U}_y\right)}{\∂t^2}+b_{21}\frac{\∂\left(U_x-u_x\right)}{\∂t}+b_{22}\frac{\∂\left(U_y-u_y\right)}{\∂t}+b_{23}\frac{\∂\left(U_z-u_z\right)}{\∂t}\end{array}---\left(1.16e\right)$

$\begin{array}{c}{Q}_1\frac{{\∂}^2{u}_x}{\∂x\∂z}+Q_2\frac{{\∂}^2{u}_y}{\∂y\∂z}+Q_3\frac{{\∂}^2{u}_z}{\∂z^2}+Q_4(\frac{{\∂}^2{u}_y}{\∂z^2}+\frac{{\∂}^2{u}_z}{\∂y\∂z})+Q_5(\frac{{\∂}^2{u}_z}{\∂x\∂z}+\frac{{\∂}^2{u}_x}{\∂z^2})+Q_6(\frac{{\∂}^2{u}_x}{\∂y\∂z}+\frac{{\∂}^2{u}_y}{\∂x\∂z})\\ +R\frac{\∂}{\∂z}(\frac{\∂U_x}{\∂x}+\frac{\∂U_y}{\∂y}-+\frac{\∂U_z}{\∂z})\\ =\frac{{\∂}^2\left({\mathrm{\ρ}}_{12}{u}_z+{\mathrm{\ρ}}_{22}{U}_z\right)}{\∂t^2}+b_{31}\frac{\∂\left(U_x-u_x\right)}{\∂t}+b_{32}\frac{\∂\left(U_y-u_y\right)}{\∂t}+b_{33}\frac{\∂\left(U_z-u_z\right)}{\∂t}\end{array}---\left(1.16f\right)$

当介质为双相各向同性介质时，耗散系数矩阵变为公式(1.17)，如下所示：

$\left[\begin{array}{ccc}b& 0& 0\\ 0& b& 0\\ 0& 0& b\end{array}\right]---\left(1.17\right)$

其中，b为耗散系数。

将公式(1.7)和(1.9)公式代入用应力表示的运动微分方程(即公式(1.10)和公式(1.11))，略去外力并整理，可得到双相各向同性介质中用位移表示的运动微分方程，如公式(1.18a)、(1.18b)、(1.18c)、(1.18d)、(1.18e)、(1.18f)所示：

$N{\▿}^2{u}_x+(A+N)\frac{\∂\mathrm{\θ}}{\∂x}+Q\frac{\∂\mathrm{\ϵ}}{\∂x}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}{u}_x+{\mathrm{\ρ}}_{12}{U}_x)+b\frac{\∂(u_x-U_x)}{\∂t}---\left(1.18a\right)$

$Q\frac{\∂\mathrm{\θ}}{\∂x}+R\frac{\∂\mathrm{\ϵ}}{\∂x}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}{u}_x+{\mathrm{\ρ}}_{22}{U}_x)-b\frac{\∂(u_x-U_x)}{\∂t}---\left(1.18b\right)$

$N{\▿}^2{u}_y+(A+N)\frac{\∂\mathrm{\θ}}{\∂y}+Q\frac{\∂\mathrm{\ϵ}}{\∂y}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}{u}_y+{\mathrm{\ρ}}_{12}{U}_y)+b\frac{\∂(u_y-U_y)}{\∂t}---\left(1.18c\right)$

$Q\frac{\∂\mathrm{\θ}}{\∂y}+R\frac{\∂\mathrm{\ϵ}}{\∂y}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}{u}_y+{\mathrm{\ρ}}_{22}{U}_y)-b\frac{\∂(u_y-U_y)}{\∂t}---\left(1.18d\right)$

$N{\▿}^2{u}_z+(A+N)\frac{\∂\mathrm{\θ}}{\∂z}+Q\frac{\∂\mathrm{\ϵ}}{\∂z}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}{u}_z+{\mathrm{\ρ}}_{12}{U}_z)+b\frac{\∂(u_z-U_z)}{\∂t}---\left(1.18e\right)$

$Q\frac{\∂\mathrm{\θ}}{\∂z}+R\frac{\∂\mathrm{\ϵ}}{\∂z}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}{u}_z+{\mathrm{\ρ}}_{22}{U}_z)-b\frac{\∂(u_z-U_z)}{\∂t}---\left(1.18f\right)$

公式(1.18a)、(1.18b)、(1.18c)、(1.18d)、(1.18e)、(1.18f)写成向量形式，如公式(1.19a)、(1.19b)所示：

$N{\▿}^{2}u+\▿\[(A+N)\mathrm{\θ}+Q\mathrm{\ϵ}\]=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}u+{\mathrm{\ρ}}_{12}U)+b\frac{\∂(u-U)}{\∂t}---\left(1.19a\right)$

$\▿(Q\mathrm{\θ}+R\mathrm{\ϵ})=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}u+{\mathrm{\ρ}}_{22}U)-b\frac{\∂(u-U)}{\∂t}---\left(1.19b\right)$

当忽略流体与固体之间相对运动所产生的耗散力时，令b＝0便得到保守系中(无耗散情况下)用位移表示的运动微分方程式，如公式(1.20a)、(1.20b)、(1.20c)、(1.20d)、(1.20e)、(1.20f)所示：

$N{\▿}^2{u}_x+(A+N)\frac{\∂\mathrm{\θ}}{\∂x}+Q\frac{\∂\mathrm{\ϵ}}{\∂x}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}{u}_x+{\mathrm{\ρ}}_{12}{U}_x)---\left(1.20a\right)$

$Q\frac{\∂\mathrm{\θ}}{\∂x}+R\frac{\∂\mathrm{\ϵ}}{\∂x}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}{u}_x+{\mathrm{\ρ}}_{22}{U}_x)---\left(1.20b\right)$

$N{\▿}^2{u}_y+(A+N)\frac{\∂\mathrm{\θ}}{\∂y}+Q\frac{\∂\mathrm{\ϵ}}{\∂y}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}{u}_y+{\mathrm{\ρ}}_{12}{U}_y)---\left(1.20c\right)$

$Q\frac{\∂\mathrm{\θ}}{\∂y}+R\frac{\∂\mathrm{\ϵ}}{\∂y}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}{u}_y+{\mathrm{\ρ}}_{22}{U}_y)---\left(1.20d\right)$

$N{\▿}^2{u}_z+(A+N)\frac{\∂\mathrm{\θ}}{\∂z}+Q\frac{\∂\mathrm{\ϵ}}{\∂z}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}{u}_z+{\mathrm{\ρ}}_{12}{U}_z)---\left(1.20e\right)$

$Q\frac{\∂\mathrm{\θ}}{\∂z}+R\frac{\∂\mathrm{\ϵ}}{\∂z}=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}{u}_z+{\mathrm{\ρ}}_{22}{U}_z)---\left(1.20f\right)$

公式(1.20a)、(1.20b)、(1.20c)、(1.20d)、(1.20e)、(1.20f)写成向量形式，如公式(1.21a)、(1.21b)所示：

$N{\▿}^{2}u+\▿\[(A+N)\mathrm{\θ}+Q\mathrm{\ϵ}\]=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}u+{\mathrm{\ρ}}_{12}U)---\left(1.21a\right)$

$\▿(Q\mathrm{\θ}+R\mathrm{\ϵ})=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}u+{\mathrm{\ρ}}_{22}U)---\left(1.21b\right)$

1.4双相各向同性介质中的纵波方程

假设介质在统计意义上具有各向同性特征，则纵波与横波之间无耦合作用，各自满足相应的波动方程。

记▽×u＝θ，▽×U＝ε，对公式(1.19)两边取散度，得公式(1.22)，如下所示：

$\left\{\begin{array}{c}{\▿}^2\left(P\mathrm{\θ}+Q\mathrm{\ϵ}\right)=\frac{\∂2}{\∂t^2}\left({\mathrm{\ρ}}_{11}\mathrm{\θ}+{\mathrm{\ρ}}_{12}\mathrm{\ϵ}\right)+b\frac{\∂}{\∂t}\left(\mathrm{\θ}-\mathrm{\ϵ}\right)\\ {\▿}^2\left(Q\mathrm{\θ}+R\mathrm{\ϵ}\right)=\frac{\∂2}{\∂t^2}\left({\mathrm{\ρ}}_{12}\mathrm{\θ}+{\mathrm{\ρ}}_{22}\mathrm{\ϵ}\right)-b\frac{\∂}{\∂t}\left(\mathrm{\θ}-\mathrm{\ϵ}\right)\end{array}\right.---\left(1.22\right)$

这即为双相各向同性介质中纵波传播所满足的方程的一般形式。

令b＝0便得到保守系中双相各向同性介质中的纵波方程，如公式(1.23)所示：

$\left\{\begin{array}{c}{\▿}^2(P\mathrm{\θ}+Q\mathrm{\ϵ})=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}\mathrm{\θ}+{\mathrm{\ρ}}_{12}\mathrm{\ϵ})\\ {\▿}^2(Q\mathrm{\θ}+R\mathrm{\ϵ})=\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}\mathrm{\θ}+{\mathrm{\ρ}}_{22}\mathrm{\ϵ})\end{array}\right.---\left(1.23\right)$

其中，P＝A+2N。

为解公式(1.23)，为方便起见，考虑平面波。设沿x方向传播的平面波的形式，如公式(1.24)所示：

$\left\{\begin{array}{c}\mathrm{\θ}=c_1\exp \[i(lx+\mathrm{\α}t)\]\\ \mathrm{\ϵ}=c_2\exp \[i(lx+\mathrm{\α}t)\]\end{array}\right.---\left(1.24\right)$

这些波的传播速度为公式(1.25)，如下所示：

V＝α/l(1.25)

为方便计，定义参考速度V_c，如公式(1.26)：

V_c＝H/ρ(1.26)

式中，H＝P+R+2Q。V_c代表流体和固体同步运动(即θ＝ε)时纵波的传播速度。

为方便起见，进一步引入下列无量纲参量，如公式(1.27)所示：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{11}=\frac{P}{H},{\mathrm{\σ}}_{22}=\frac{R}{H},{\mathrm{\σ}}_{12}=\frac{Q}{H}\\ {\mathrm{\γ}}_{11}=\frac{{\mathrm{\ρ}}_{11}}{\mathrm{\ρ}},{\mathrm{\γ}}_{22}=\frac{{\mathrm{\ρ}}_{22}}{\mathrm{\ρ}},{\mathrm{\γ}}_{12}=\frac{{\mathrm{\ρ}}_{12}}{\mathrm{\ρ}}\end{array}\right.---\left(1.27\right)$

公式(1.27)满足公式(1.28)，如下所示：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{11}{\mathrm{\σ}}_{22}-{\mathrm{\σ}}_{12}^2>0\\ {\mathrm{\γ}}_{11}{\mathrm{\γ}}_{22}-{\mathrm{\γ}}_{12}^2>0\\ {\mathrm{\σ}}_{11}+{\mathrm{\σ}}_{22}+2{\mathrm{\σ}}_{12}=1\\ {\mathrm{\γ}}_{11}+{\mathrm{\γ}}_{22}+2{\mathrm{\γ}}_{12}=1\end{array}\right.---\left(1.28\right)$

将公式(1.24)代入公式(1.23)，并令公式(1.28)，如下所示：

$z=\frac{V_c^2}{V^2},---\left(1.29\right)$

以得公式(1.30)，如下所示：

$\{\begin{array}{c}z({\mathrm{\σ}}_{11}{c}_1+{\mathrm{\σ}}_{12}{c}_2)={\mathrm{\γ}}_{11}{c}_1+{\mathrm{\γ}}_{12}{c}_2\\ z({\mathrm{\σ}}_{12}{c}_1+{\mathrm{\σ}}_{22}{c}_2)={\mathrm{\γ}}_{12}{c}_1+{\mathrm{\γ}}_{22}{c}_2\end{array}.---\left(1.30\right)$

公式(1.30)可写成以c₁、c₂为未知数的线性方程组，如公式(1.31)所示：

$\left[\begin{array}{cc}{\mathrm{z\σ}}_{11}-{\mathrm{\γ}}_{11}& {\mathrm{z\σ}}_{12}-{\mathrm{\γ}}_{12}\\ {\mathrm{z\σ}}_{12}-{\mathrm{\γ}}_{12}& {\mathrm{z\σ}}_{22}-{\mathrm{\γ}}_{22}\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=0,---\left(1.31\right)$

要求公式(1.31)有非零解，必须满足公式(1.32)，如下所示：

$\left|\begin{array}{cc}{\mathrm{z\σ}}_{11}-{\mathrm{\γ}}_{11}& {\mathrm{z\σ}}_{12}-{\mathrm{\γ}}_{12}\\ {\mathrm{z\σ}}_{12}-{\mathrm{\γ}}_{12}& {\mathrm{z\σ}}_{22}-{\mathrm{\γ}}_{22}\end{array}\right|=0,---\left(1.32\right)$

展开并整理，得公式(1.33)，如下所示：

$({\mathrm{\σ}}_{11}{\mathrm{\σ}}_{22}-{\mathrm{\σ}}_{12}^2)z^2-({\mathrm{\σ}}_{11}{\mathrm{\γ}}_{22}+{\mathrm{\σ}}_{22}{\mathrm{\γ}}_{11}-2{\mathrm{\σ}}_{12}{\mathrm{\γ}}_{12})z+({\mathrm{\γ}}_{11}{\mathrm{\γ}}_{22}-{\mathrm{\γ}}_{12}^2)=0.---\left(1.33\right)$

求解上面的方程，得到方程的两个根z₁和z₂，且z₁和z₂均大于0。因此，满足公式(1.23)的波有两个，其速度分别为公式(1.34)和公式(1.35)，如下所示：

$V_1=\sqrt{\frac{V_c^2}{z_1}}---\left(1.34\right)$

$V_2=\sqrt{\frac{V_c^2}{z_2}}---\left(1.35\right)$

由这两种波的振幅关系还可证明有下列正交关系存在，如公式(1.36)所示：

${\mathrm{\γ}}_{11}{c}_1^1{c}_1^2+2{\mathrm{\γ}}_{12}(c_1^1{c}_2^2+c_2^1{c}_1^2)+{\mathrm{\γ}}_{22}{c}_2^1{c}_2^2=0,---\left(1.36\right)$

从而可以得出，对一种纵波，固相位移和流相位移是同相的，而对另一种纵波，固相位移和流相位移是反相的。由公式(1.30)可进一步导出下列关系，如公式(1.37)所示：

$z_i=\frac{{\mathrm{\γ}}_{11}{\left(c_1^i\right)}^2+2{\mathrm{\γ}}_{12}\left(c_1^i\right)\left(c_2^i\right)+{\mathrm{\γ}}_{22}{\left(c_2^i\right)}^2}{{\mathrm{\σ}}_{11}{\left(c_1^i\right)}^2+2{\mathrm{\σ}}_{12}\left(c_1^i\right)\left(c_2^i\right)+{\mathrm{\σ}}_{22}{\left(c_2^i\right)}^2},---\left(1.37\right)$

因为γ₁₂是唯一的可以为负的系数，所以速度高的波固相位移和流相位移是同相的，速度低的波固相位移和流相位移是反相的。一般称速度高的纵波为第一类纵波(快纵波)，称速度较低的波为第二类纵波(慢纵波)。

对于有耗散的纵波方程，如公式(1.22)，同样考虑平面波，令公式(1.38)，如下所示：

$\{\begin{array}{c}\mathrm{\θ}=c_1\exp \[i(lx+\mathrm{\α}t)\]\\ \mathrm{\ϵ}=c_2\exp \[i(lx+\mathrm{\α}t)\]\end{array},---\left(1.38\right)$

将公式(1.37)代入公式(1.22)，并消去常量c₁、c₂，得到下面的关系，如公式(1.39)所示：

$(PR-Q^2)\frac{l^4}{{\mathrm{\α}}^4}-({\mathrm{R\ρ}}_{11}+{\mathrm{P\ρ}}_{22}-2{\mathrm{Q\ρ}}_{12})\frac{l^2}{{\mathrm{\α}}^2}+{\mathrm{\ρ}}_{11}{\mathrm{\ρ}}_{22}-{\mathrm{\ρ}}_{12}^2+\frac{ib}{\mathrm{\α}}\[(P+R+2Q)\frac{l^2}{{\mathrm{\α}}^2}-\mathrm{\ρ}\]=0---\left(1.39\right)$

公式(1.39)写成无量纲形式，如公式(1.40)所示：

$({\mathrm{\σ}}_{11}{\mathrm{\σ}}_{22}-{\mathrm{\σ}}_{12}^2)z^2-({\mathrm{\σ}}_{22}{\mathrm{\γ}}_{11}+{\mathrm{\σ}}_{11}{\mathrm{\γ}}_{22}-2{\mathrm{\σ}}_{12}{\mathrm{\γ}}_{12})z+({\mathrm{\γ}}_{11}{\mathrm{\γ}}_{22}-{\mathrm{\γ}}_{12}^2)+\frac{ib}{\mathrm{\α}\mathrm{\ρ}}(z-1)=0---\left(1.40\right)$

在公式(1.40)中令b＝0，便得到与公式(1.33)相同的式子。利用公式(1.33)的两个根z₁和z₂，公式(1.40)可写为公式(1.41)，如下所示：

(z-z₁)(z-z₂)+iM(z-1)＝0，(1.41)

其中，M定义为公式(1.42)，如下所示：

$M=\frac{b}{\mathrm{\α}\mathrm{\ρ}({\mathrm{\σ}}_{11}{\mathrm{\σ}}_{22}-{\mathrm{\σ}}_{12}^2)}.---\left(1.42\right)$

且定义公式(1.43)，如下所示：

$f_c=\frac{b}{2{\mathrm{\π\ρ}}_2}=\frac{b}{2\mathrm{\π}\mathrm{\ρ}({\mathrm{\γ}}_{12}+{\mathrm{\γ}}_{22})}---\left(1.43\right)$

其中，f_c称为特征频率。因为α＝2πf(f为波的频率)，所以M又可写成的函数，如公式(1.44)所示：

$M=\frac{f}{f_c}\frac{({\mathrm{\γ}}_{12}+{\mathrm{\γ}}_{22})}{({\mathrm{\σ}}_{11}{\mathrm{\σ}}_{22}-{\mathrm{\σ}}_{12}^2)}.---\left(1.44\right)$

设Z_I和Z_II为方程(1.41)的两个根，其中Z_I→1(当f＝0时)，则它对应于第一类纵波，而Z_II对应于第二类纵波。

以上是基于Biot理论的双相介质中的纵波方程，可见，不论在保守系中还是在有耗散系中，在双相介质中均存在两类纵波，即第一类纵波(快纵波)和第二类纵波(慢纵波)；对第一类纵波固相位移和流相位移是同相的，对第二类纵波固相位移和流相位移是反相的；第一类纵波与单相介质中的纵波具有相同的性质，第二类纵波的传播类似于扩散现象或热传导现象，具有很强的频散和衰减性质。

1.5双相各向同性介质中的横波方程

记▽×u＝w，▽×U＝Ω对公式(1.19)两边取旋度，得双相各向同性介质中横波所满足的波动方程，如公式(1.45)所示：

$\left\{\begin{array}{c}\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}w+{\mathrm{\ρ}}_{12}\mathrm{\Ω})+b\frac{\∂}{\∂t}(w-\mathrm{\Ω})=N{\▿}^{2}w\\ \frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}w+{\mathrm{\ρ}}_{22}\mathrm{\Ω})-b\frac{\∂}{\∂t}(w-\mathrm{\Ω})=0\end{array}\right.---\left(1.45\right)$

令b＝0，便得到保守系中双相各向同性介质中的横波方程，如公式(1.46)所示：

$\left\{\begin{array}{c}\frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{11}w+{\mathrm{\ρ}}_{12}\mathrm{\Ω})=N{\▿}^{2}w\\ \frac{{\∂}^2}{\∂t^2}({\mathrm{\ρ}}_{12}w+{\mathrm{\ρ}}_{22}\mathrm{\Ω})=0\end{array}\right.---\left(1.46\right)$

要求解方程(1.46)，消去方程组中的Ω，得到公式(1.47)，如下所示：

$N{\▿}^{2}w={\mathrm{\ρ}}_{11}(1-\frac{{\mathrm{\ρ}}_{12}^2}{{\mathrm{\ρ}}_{11}{\mathrm{\ρ}}_{22}})\frac{{\∂}^{2}w}{\∂t^2},---\left(1.47\right)$

所以，公式中仅存在一种类型的旋转波，其传播速度如公式(1.48)所示：

$V_s={\[\frac{N}{{\mathrm{\ρ}}_{11}(1-{\mathrm{\ρ}}_{12}^2/{\mathrm{\ρ}}_{11}{\mathrm{\ρ}}_{22})}\]}^{\frac{1}{2}}---\left(1.48\right)$

由公式(1.46)中的第二个式子可以知道，固体的旋转与流体的旋转成比例地耦合，如公式(1.49)所示：

$\mathrm{\Ω}=-\frac{{\mathrm{\ρ}}_{12}}{{\mathrm{\ρ}}_{22}}w,---\left(1.49\right)$

因为ρ₁₂≤0，ρ₂₂>0，所以固体与流体的旋转是同相的。

对于有耗散系中的横波方程(即公式(1.45))，为方便，考虑平面波。设沿x方向传播，固体和流体的偏振方向为z，传播的平面波的形式如公式(1.50)所示：

$\left\{\begin{array}{c}w=c_1\exp \[i(lx+\mathrm{\α}t)\]\\ \mathrm{\Ω}=c_2\exp \[i(lx+\mathrm{\α}t)\]\end{array}\right.,---\left(1.50\right)$

将公式(1.50)代入公式(1.45)，消去常量c₁和c²，得到公式(1.51)，如下所示：

Nl²/ρα²＝E_r-iE_i，(1.51)

其中，Er及Ei定义如公式(1.52)所示：

$\left\{\begin{array}{c}{E}_r=\frac{1+{\mathrm{\γ}}_{22}\frac{{\mathrm{\γ}}_{11}{\mathrm{\γ}}_{22}-{\mathrm{\γ}}_{12}^2}{{\left({\mathrm{\γ}}_{12}+{\mathrm{\γ}}_{22}\right)}^2}{\left(\frac{f}{f_c}\right)}^2}{1+{\left(\frac{{\mathrm{\γ}}_{22}}{{\mathrm{\γ}}_{12}+{\mathrm{\γ}}_{22}}\right)}^2{\left(\frac{f}{f_c}\right)}^2}\\ E_i=\frac{f}{f_c}\frac{\left({\mathrm{\γ}}_{12}+{\mathrm{\γ}}_{22}\right)}{1+{\left(\frac{{\mathrm{\γ}}_{22}}{{\mathrm{\γ}}_{12}+{\mathrm{\γ}}_{22}}\right)}^2{\left(\frac{f}{f_c}\right)}^2}\end{array}\right..---\left(1.52\right)$

波的频率为令l＝l_r+il_i，则波的相速度为公式(1.53)，如下所示：

v_r＝α/|l_r|，(1.53)

引入一个参考速度，如公式(1.54)所示：

$V_r={(N/\mathrm{\ρ})}^{\frac{1}{2}},---\left(1.54\right)$

此速度代表流体和固体之间无相对运动时旋转波的速度。

从公式(1.51)可得到公式(1.55)，如下所示：

$\frac{v_r}{V_r}=\sqrt{2}/{\[{(E_r^2+E_i^2)}^{\frac{1}{2}}+E_r\]}^{\frac{1}{2}},---\left(1.55\right)$

所以可得到公式(1.56)，如下所示：

$v_r=\sqrt{2}{V}_r/{\[{(E_r^2+E_i^2)}^{\frac{1}{2}}+E_r\]}^{\frac{1}{2}},---\left(1.56\right)$

这即为有耗散时双相各向同性介质中横波的相速度，它是频率比f/f_c和动力学参数γ_ij的函数。

所以，在双相各向同性介质中存在一类横波，这一横波在固体介质中和流体介质中同时存在，它们耦合在一起，这与单相介质是不同的，在单相介质中横波不能在流体中传播。这进一步说明了在双相介质中由于流体介质与固体介质的相互作用使地震波的传播规律发生了变化，用经典的单相介质理论解决地质问题特别是油气储集层问题并不适宜，而双相介质理论更符合实际情况，必能更加准确有效地解决油气储层问题，必将比建立在单相介质理论基础上的地震波场理论带来更高的经济效益。

2二维双相各向同性介质中纵波方程差分格式的构造

在公式(1.23)中，令对y的偏导数等于零，则得到保守系中双相各向同性介质的二维纵波方程，如公式(1.57)所示：

$\left\{\begin{array}{c}\frac{{\∂}^2\left(P\mathrm{\θ}+Q\mathrm{\ϵ}\right)}{\∂x^2}+\frac{{\∂}^2\left(P\mathrm{\θ}+Q\mathrm{\ϵ}\right)}{\∂z^2}=\frac{{\∂}^2}{\∂t^2}\left({\mathrm{\ρ}}_{11}\mathrm{\θ}+{\mathrm{\ρ}}_{12}\mathrm{\ϵ}\right)\\ \frac{{\∂}^2\left(Q\mathrm{\θ}+R\mathrm{\ϵ}\right)}{\∂x^2}+\frac{{\∂}^2\left(Q\mathrm{\θ}+R\mathrm{\ϵ}\right)}{\∂z^2}=\frac{{\∂}^2}{\∂t^2}\left({\mathrm{\ρ}}_{12}\mathrm{\θ}+{\mathrm{\ρ}}_{22}\mathrm{\ϵ}\right)\end{array}\right.---\left(1.57\right)$

假设介质分块均匀，则公式(1.57)又可写成公式(1.58)，如下所示：

$\left\{\begin{array}{c}P\left(\frac{{\∂}^2\mathrm{\θ}}{\∂x^2}+\frac{{\∂}^2\mathrm{\θ}}{\∂z^2}\right)+Q\left(\frac{{\∂}^2\mathrm{\ϵ}}{\∂x^2}+\frac{{\∂}^2\mathrm{\ϵ}}{\∂z^2}\right)=\frac{{\∂}^2}{\∂t^2}\left({\mathrm{\ρ}}_{11}\mathrm{\θ}+{\mathrm{\ρ}}_{12}\mathrm{\ϵ}\right)\\ Q\left(\frac{{\∂}^2\mathrm{\θ}}{\∂x^2}+\frac{{\∂}^2\mathrm{\θ}}{\∂z^2}\right)+R\left(\frac{{\∂}^2\mathrm{\ϵ}}{\∂x^2}+\frac{{\∂}^2\mathrm{\ϵ}}{\∂z^2}\right)=\frac{{\∂}^2}{\∂t^2}\left({\mathrm{\ρ}}_{12}\mathrm{\θ}+{\mathrm{\ρ}}_{22}\mathrm{\ϵ}\right)\end{array}\right.---\left(1.58\right)$

对空间偏导数采用任意偶数2N(N为正整数)阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，有如下的离散格式，如公式(1.59a)和(1.59b)所示：

$\begin{array}{c}P\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i,j+m\right)}{h^2}\right)+Q\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i,j+m\right)}{h^2}\right)\\ ={\mathrm{\ρ}}_{11}\frac{{\mathrm{\θ}}^{n+1}\left(i,j\right)+{\mathrm{\θ}}^{n-1}\left(i,j\right)-2{\mathrm{\θ}}^n\left(i,j\right)}{{\mathrm{\Δt}}^2}+{\mathrm{\ρ}}_{12}\frac{{\mathrm{\ϵ}}^{n+1}\left(i,j\right)+{\mathrm{\ϵ}}^{n-1}\left(i,j\right)-2{\mathrm{\ϵ}}^n\left(i,j\right)}{{\mathrm{\Δt}}^2}\end{array}---\left(1.59a\right)$

$\begin{array}{c}Q\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i,j+m\right)}{h^2}\right)+R\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i,j+m\right)}{h^2}\right)\\ ={\mathrm{\ρ}}_{12}\frac{{\mathrm{\θ}}^{n+1}\left(i,j\right)+{\mathrm{\θ}}^{n-1}\left(i,j\right)-2{\mathrm{\θ}}^n\left(i,j\right)}{{\mathrm{\Δt}}^2}+{\mathrm{\ρ}}_{22}\frac{{\mathrm{\ϵ}}^{n+1}\left(i,j\right)+{\mathrm{\ϵ}}^{n-1}\left(i,j\right)-2{\mathrm{\ϵ}}^n\left(i,j\right)}{{\mathrm{\Δt}}^2}\end{array}---\left(1.59b\right)$

将公式(1.59a)和公式(1.59b)进行整理，得公式(1.60a)和(1.60b)，如下所示：

$\begin{array}{c}{\mathrm{\ρ}}_{11}{\mathrm{\θ}}^{n+1}\left(i,j\right)+{\mathrm{\ρ}}_{12}{\mathrm{\ϵ}}^{n+1}\left(i,j\right)\\ =P\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ +Q\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ -{\mathrm{\ρ}}_{11}\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-2{\mathrm{\θ}}^n\left(i,j\right)\]-{\mathrm{\ρ}}_{12}\[{\mathrm{\ϵ}}^{n-1}\left(i,j\right)-2{\mathrm{\ϵ}}^n\left(i,j\right)\]\end{array}---\left(1.60a\right)$

$\begin{array}{c}{\mathrm{\ρ}}_{12}{\mathrm{\θ}}^{n+1}\left(i,j\right)+{\mathrm{\ρ}}_{22}{\mathrm{\ϵ}}^{n+1}\left(i,j\right)\\ =Q\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ +R\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ -{\mathrm{\ρ}}_{12}\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-2{\mathrm{\θ}}^n\left(i,j\right)\]-{\mathrm{\ρ}}_{22}\[{\mathrm{\ϵ}}^{n-1}\left(i,j\right)-2{\mathrm{\ϵ}}^n\left(i,j\right)\]\end{array}---\left(1.60b\right)$

其中，i是x方向的空间序号，j是z方向的空间序号；角标n表示时刻；h为空间差分步长(x方向和z方向的离散步长可相等，也可不相等，这里取为相等)；θⁿ(i,j)表示固相体应变θ在n时刻在(ih,jh)处的值；εⁿ(i,j)表示流相体应变ε在n时刻在(ih,jh)处的值；为差分系数，2N为差分精度；当N＝1时，有 $c_0^{\left(1\right)}=-2.0000,c_1^{\left(1\right)}=1.0000;$ 当N＝2时，有 $c_0^{\left(2\right)}=-2.5000,$ $c_1^{\left(2\right)}=1.3333,c_2^{\left(2\right)}=-0.0833;$ 当N＝4时，有 $c_0^{\left(4\right)}=-2.8472,c_1^{\left(4\right)}=1.6000,c_2^{\left(4\right)}=-2.000,$ $c_3^{\left(4\right)}=0.0254,c_4^{\left(4\right)}=-0.0018;$ △t是时间采样间隔。在实际计算中，为防止近边界网格点上差分阶数的突然降低而带来较大误差，在近边界处采用逐步降低差分阶数的方法。

求解公式(1.60a)和公式(1.60b)，得公式(1.61)，如下所示：

$\{\begin{array}{c}{\mathrm{\θ}}^{n+1}(i,j)=\frac{{\mathrm{\ρ}}_{22}A-{\mathrm{\ρ}}_{12}B}{{\mathrm{\ρ}}_{11}{\mathrm{\ρ}}_{22}-{\mathrm{\ρ}}_{12}^2}\\ {\mathrm{\ϵ}}^{n+1}(i,j)=\frac{{\mathrm{\ρ}}_{11}B-{\mathrm{\ρ}}_{12}A}{{\mathrm{\ρ}}_{11}{\mathrm{\ρ}}_{22}-{\mathrm{\ρ}}_{12}^2}\end{array},---\left(1.61\right)$

其中，A和B分别定义如公式(1.62)和公式(1.63)所示：

$\begin{array}{c}A=P\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ +Q\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ -{\mathrm{\ρ}}_{11}\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-2{\mathrm{\θ}}^n\left(i,j\right)\]-{\mathrm{\ρ}}_{12}\[{\mathrm{\ϵ}}^{n-1}\left(i,j\right)-2{\mathrm{\ϵ}}^n\left(i,j\right)\]\end{array},---\left(1.62\right)$

$\begin{array}{c}B=Q\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ +R\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ -{\mathrm{\ρ}}_{12}\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-2{\mathrm{\θ}}^n\left(i,j\right)\]-{\mathrm{\ρ}}_{22}\[{\mathrm{\ϵ}}^{n-1}\left(i,j\right)-2{\mathrm{\ϵ}}^n\left(i,j\right)\]\end{array}.---\left(1.63\right)$

同样，在公式(1.22)中，令对y的偏导数等于零，则得到有耗散时双相各向同性介质中的二维纵波方程，如公式(1.64)所示：

$\left\{\begin{array}{c}\frac{{\∂}^2\left(P\mathrm{\θ}+Q\mathrm{\ϵ}\right)}{\∂x^2}+\frac{{\∂}^2\left(P\mathrm{\θ}+Q\mathrm{\ϵ}\right)}{\∂z^2}=\frac{{\∂}^2}{\∂t^2}\left({\mathrm{\ρ}}_{11}\mathrm{\θ}+{\mathrm{\ρ}}_{12}\mathrm{\ϵ}\right)+b\frac{\∂}{\∂t}\left(\mathrm{\θ}-\mathrm{\ϵ}\right)\\ \frac{{\∂}^2\left(Q\mathrm{\θ}+R\mathrm{\ϵ}\right)}{\∂x^2}+\frac{{\∂}^2\left(Q\mathrm{\θ}+R\mathrm{\ϵ}\right)}{\∂z^2}=\frac{{\∂}^2}{\∂t^2}\left({\mathrm{\ρ}}_{12}\mathrm{\θ}+{\mathrm{\ρ}}_{22}\mathrm{\ϵ}\right)-b\frac{\∂}{\∂t}\left(\mathrm{\θ}-\mathrm{\ϵ}\right)\end{array}\right..---\left(1.64\right)$

介质分块均匀时，公式(1.64)又可写为公式(1.65)，如下所示：

$\left\{\begin{array}{c}P\left(\frac{{\∂}^2\mathrm{\θ}}{\∂x^2}+\frac{{\∂}^2\mathrm{\θ}}{\∂z^2}\right)+Q\left(\frac{{\∂}^2\mathrm{\ϵ}}{\∂x^2}+\frac{{\∂}^2\mathrm{\ϵ}}{\∂z^2}\right)=\frac{{\∂}^2}{\∂t^2}\left({\mathrm{\ρ}}_{11}\mathrm{\θ}+{\mathrm{\ρ}}_{12}\mathrm{\ϵ}\right)+b\frac{\∂}{\∂t}\left(\mathrm{\θ}-\mathrm{\ϵ}\right)\\ Q\left(\frac{{\∂}^2\mathrm{\θ}}{\∂x^2}+\frac{{\∂}^2\mathrm{\θ}}{\∂z^2}\right)+R\left(\frac{{\∂}^2\mathrm{\ϵ}}{\∂x^2}+\frac{{\∂}^2\mathrm{\ϵ}}{\∂z^2}\right)=\frac{{\∂}^2}{\∂t^2}\left({\mathrm{\ρ}}_{12}\mathrm{\θ}+{\mathrm{\ρ}}_{22}\mathrm{\ϵ}\right)-b\frac{\∂}{\∂t}\left(\mathrm{\θ}-\mathrm{\ϵ}\right)\end{array}\right.---\left(1.65\right)$

采用与公式(1.58)同样的方法将公式(1.65)进行差分离散，略去推导过程，得到如下的离散格式，如公式(1.66)所示：

$\left\{\begin{array}{c}{\mathrm{\θ}}^{n+1}(i,j)=\frac{(2{\mathrm{\ρ}}_{22}+b\mathrm{\Δ}t)A^{\′}-(2{\mathrm{\ρ}}_{12}-b\mathrm{\Δ}t)B^{\′}}{(2{\mathrm{\ρ}}_{11}+b\mathrm{\Δ}t)(2{\mathrm{\ρ}}_{22}+b\mathrm{\Δ}t)-{(2{\mathrm{\ρ}}_{12}-b\mathrm{\Δ}t)}^2}\\ {\mathrm{\ϵ}}^{n+1}(i,j)=\frac{(2{\mathrm{\ρ}}_{11}+b\mathrm{\Δ}t)B^{\′}-(2{\mathrm{\ρ}}_{12}-b\mathrm{\Δ}t)A^{\′}}{(2{\mathrm{\ρ}}_{11}+b\mathrm{\Δ}t)(2{\mathrm{\ρ}}_{22}+b\mathrm{\Δ}t)-{(2{\mathrm{\ρ}}_{12}-b\mathrm{\Δ}t)}^2}\end{array}\right.,---\left(1.66\right)$

其中，A’和B’分别定义如公式(1.67)和公式(1.68)所示：

$\begin{array}{c}{A}^{\′}=2\{P\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ +Q\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ -{\mathrm{\ρ}}_{11}\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-2{\mathrm{\θ}}^n\left(i,j\right)\]-{\mathrm{\ρ}}_{12}\[{\mathrm{\ϵ}}^{n-1}\left(i,j\right)-2{\mathrm{\ϵ}}^n\left(i,j\right)\]\}\\ +b\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-{\mathrm{\ϵ}}^{n-1}\left(i,j\right)\]\mathrm{\Δ}t\\ =2A+b\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-{\mathrm{\ϵ}}^{n-1}\left(i,j\right)\]\mathrm{\Δ}t\end{array},---\left(1.67\right)$

$\begin{array}{c}{B}^{\′}=2\{Q\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\θ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ +R\left(\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\Σ}}}\frac{c_m^{\left(N\right)}{\mathrm{\ϵ}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\Δt}}^2\\ -{\mathrm{\ρ}}_{11}\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-2{\mathrm{\θ}}^n\left(i,j\right)\]-{\mathrm{\ρ}}_{22}\[{\mathrm{\ϵ}}^{n-1}\left(i,j\right)-2{\mathrm{\ϵ}}^n\left(i,j\right)\]\}\\ -b\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-{\mathrm{\ϵ}}^{n-1}\left(i,j\right)\]\mathrm{\Δ}t\\ =2B-b\[{\mathrm{\θ}}^{n-1}\left(i,j\right)-{\mathrm{\ϵ}}^{n-1}\left(i,j\right)\]\mathrm{\Δ}t\end{array}.---\left(1.68\right)$

图1a和图1b分別是采用以上差分格式计算的均匀双相各向同性介质中固相波场和流相波场的快照的示意图，差分精度为8阶。图1a和图1b中，第一纵波P₁和第二纵波P₂清晰可见，且第一纵波P₁和第二纵波P₂的速度与理论分析结果完全一致，表明本发明所述的算法是正确的。图2a和图2b分別是双相各向同性介质中固相波场和流相波场的单炮记录的示意图，采用的观测系统为：中间放炮，101道接收；采样间隔：0.5ms；道间距：30m；主频60Hz。介质弹性参数见表1。

表1均匀双相各向同性介质的弹性参数

P	Q	R	ρ11	ρ22	ρ12	b
							16.025	2.2575	3.81	1850	405	-135	0.5

P,Q,R：10⁹kg.m^-1.s^-2；ρ_ij：kg.m^-3；b：10^-6kg.m^-3.s^-1

2.1稳定性条件

上面所导出的差分格式(1.61)和(1.66)的执行是通过按时间步推进来计算地震波场在计算空间内的变化规律的，因此，此差分格式必须满足地震波场传播的因果关系。也即时间变量步长△t和空间变量步长h之间必须满足一定关系，否则，将会出现数值的不稳定，随着计算步数的增加，被计算的场量数值无限制地增大。在进行双相介质中纵波场正演计算时，时间变量步长△t和空间差分步长h应满足如下稳定性条件，如公式(1.69)所示：

$\frac{v\mathrm{\&Delta;}t}{h}<1/\sqrt{2\underset{m=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{c}_{2m-1}^{\left(N\right)}}---\left(1.69\right)$

其中N₁为不超过N的最大奇数，v为所计算空间内介质的最大速度。

2.2吸收边界条件

在地球物理学所研究的领域中，常常所研究的区域是无限大的，而在无限大的区域无法用计算机求解，必须将求解区域截断，用有界区域代替无界区域。引入的人工边界会使波到达时产生反射，必须采取特殊的计算方法来消除边界的伪反射。使向区域外传播的波不在边界上产生反射的算法也称为无反射边界条件或吸收边界条件。

关于吸收边界条件的问题，前人已作了许多研究工作。目前，主要有三种吸收边界条件方法：外推法，由Taylor等人(1969)提出；外行波的模拟法，Taflove和Brodwin(1975)使用此法处理了边界条件的反射问题；衰减法，在计算区域边界内侧的一定范围内，对地震波能量进行一定程度的衰减，从而使地震波传播到边界时的能量变得较弱，与反射能量相比微乎其微。早在1980年，Taflove(1980)直接在计算空间周围设置吸收边界层，使外行波被吸收，实际就是衰减法的思想。邵治龙等(1998)也使用了衰减法对边界条件作了处理，取得了一定的效果。

目前，应用最广泛的吸收边界条件是Clayton和Engquist(1977)，Engquist和Majda(1979)提出的吸收边界条件。它是一种外行波的模拟法，其主要思路是通过波动方程的因子分解来获得单行波方程，并由此建立吸收边界条件。这种边界条件能完全吸收垂直入射到边界的反射波，二阶以上的Clayton-Engquist条件对非垂直入射到人为边界上波的吸收只能部分吸收，而且随着入射角偏离垂直方向越大，这种方法对人为边界条件反射的吸收效果越差，当入射角接近π/2时，吸收效果已很不好。针对Clayton-Engquist吸收边界条件的局限性，罗大清等(1999)给出了加权方向校正Clayton吸收边界条件。虽然罗大清等给出的加权方向校正吸收边界条件是适用于单相介质中的声波方程的，只要将其稍加改进即可适用于双相介质中的纵波方程。在截断边界处，分别将纵波与横波波场函数固相位移θ和流相位移ε当作外行平面波，求得边界处入射波的入射角后将其代入改进后的Clayton吸收边界条件中，便可得到双相各向同性中纵波方程所满足的吸收边界条件。本章以左边界为例讨论。

要使波在到达人为左边界处(x＝0)不发生反射，也即使右行波的能量为零，地震波在该处是向左传播的单向波。为获得二维问题的单向波方程，对二维波动方程算子进行因子分解。设φ(x,z,t)是二维问题中的任意一场分量，则对无源区域有波动方程，如公式(1.70)所示：

$(\frac{{\&part;}^2}{\&part;x^2}+\frac{{\&part;}^2}{\&part;z^2}-\frac{1}{v^2}\frac{{\&part;}^2}{\&part;t^2})\mathrm{\&phi;}=0.---\left(1.70\right)$

定义算子，如公式(1.71a)和公式(1.71b)所示：

$L_2^{+}=\frac{\&part;}{\&part;x}+\frac{1}{v}\frac{\&part;}{\&part;t}\sqrt{1-S^2}---\left(1.71a\right)$

$L_2^{-}\frac{\&part;}{\&part;x}-\frac{1}{v}\frac{\&part;}{\&part;t}\sqrt{1-S^2},---\left(1.71b\right)$

其中，S定義如公式(1.72)所示：

$S=v\frac{\&part;}{\&part;z}/\frac{\&part;}{\&part;t},---\left(1.72\right)$

对公式(1.70)进行因子分解，得到公式(1.73)，如下所示：

$L_2^{+}{L}_2^{-}\mathrm{\&phi;}=0.---\left(1.73\right)$

可以证明，当把L^- ₂应用于左边界x＝0(图1-3)的φ时，φ以任意角度从Ω(Ω为计算空间)内部入射到x＝0边界的平面波都会被边界所吸收。换句话说，如公式(1.74)所示：

$\begin{array}{cc}{L}_2^{-}\mathrm{\&phi;}=0& (x=0)\end{array}---\left(1.74\right)$

就是保证从Ω内部以任意角度入射到x＝0边界的平面波φ的精确解析吸收边界条件。事实上，公式(1.71)中由于根号的存在并不适合直接进行数值计算。在实际计算中，吸收边界条件是通过对精确吸收边界条件中的根号部分取近似而得到的。对公式(1.71)中的根号部分进行二阶Taylor展开，得公式(1.75)，如下所示：

$\sqrt{1-S^2}=1-\frac{S^2}{2}---\left(1.75\right)$

将公式(1.75)代入公式(1.74)，得到左边界x＝0时二阶精度近似吸收边界条件，如公式(1.76)所示：

$(\frac{{\&part;}^2}{\&part;x\&part;t}-\frac{1}{v}\frac{{\&part;}^2}{\&part;t^2}+\frac{v}{2}\frac{{\&part;}^2}{\&part;z^2})\mathrm{\&phi;}=0---\left(1.76\right)$

公式(1.76)由于破坏了原算子的严格数学关系，将不再对入射到边界上的波完全吸收，它只能完全吸收垂直入射到边界上的波，对非垂直入射的波只能部分吸收，且随着入射角偏离垂直方向越大，其对外行波的吸收性能越差。为改善边界条件对外行波的吸收性能，本研究将罗大清等(1999)的加权方向校正吸收边界条件用到双相各向同性介质中纵波场有限差分数值模拟的边界条件中。

在地震波场中，直达波的能量最强，截断边界对直达波的伪反射也强于对其它波的伪反射，大量数值计算实例也表明了这一点。在前述部分中我们已经知道，在双相各向同性介质的纵波场中存在两类纵波，第一类纵波与单相介质中的纵波具有相同的性质，第二类纵波具有很强的频散和衰减效应，在实际中很难观测到，因此，在进行数值计算时，我们主要吸收截断边界对第一纵波的伪反射。

图3为边界入射波的示意图。如图3所示，θ为第一纵波到达截断边界上一点P处的入射角，利用公式(1.32)可以得出双相各向同性介质中第一纵波的传播速度根据加权方向校正Clayton吸收边界得到下面的左边界x＝0吸收边界条件，如公式(1.77)所示：

$(\frac{{\&part;}^2}{\&part;x\&part;t}-\frac{1}{V_{p_1}}\frac{{\&part;}^2}{\&part;t^2}+\frac{V_{p_1}}{(1+cos\mathrm{\&theta;})}\frac{{\&part;}^2}{\&part;z^2})\mathrm{\&phi;}=0---\left(1.77\right)$

与此类似，可以得到其它三个边界条件，如公式(1.78)、(1.79)和(1.80)所示：

右边界

$(\frac{{\&part;}^2}{\&part;x\&part;t}+\frac{1}{V_{p_1}}\frac{{\&part;}^2}{\&part;t^2}-\frac{V_{p_1}}{(1+cos\mathrm{\&theta;})}\frac{{\&part;}^2}{\&part;z^2})\mathrm{\&phi;}=0---\left(1.78\right)$

顶边界

$(\frac{{\&part;}^2}{\&part;z\&part;t}-\frac{1}{V_{p_1}}\frac{{\&part;}^2}{\&part;t^2}+\frac{V_{p_1}}{(1+cos\mathrm{\&theta;})}\frac{{\&part;}^2}{\&part;x^2})\mathrm{\&phi;}=0---\left(1.79\right)$

底边界

$(\frac{{\&part;}^2}{\&part;z\&part;t}+\frac{1}{V_{p_1}}\frac{{\&part;}^2}{\&part;t^2}-\frac{V_{p_1}}{(1+cos\mathrm{\&theta;})}\frac{{\&part;}^2}{\&part;x^2})\mathrm{\&phi;}=0---\left(1.80\right)$

左边界条件的二阶差分格式如公式(1.81)所示：

$\begin{array}{c}{\mathrm{\&phi;}}^{n+1}\left(0,j\right)=-{\mathrm{\&phi;}}^{n-1}\left(1,j\right)+\frac{V_{P_1}\left(0,j\right)\&times;\mathrm{\&Delta;}t-h}{V_{P_1}\left(0,j\right)\&times;\mathrm{\&Delta;}t+h}\&lsqb;{\mathrm{\&phi;}}^{n+1}\left(1,j\right)+{\mathrm{\&phi;}}^{n-1}\left(0,j\right)\&rsqb;\\ +\frac{2h}{V_{P_1}\left(0,j\right)\&times;\mathrm{\&Delta;}t+h}\&lsqb;{\mathrm{\&phi;}}^n\left(1,j\right)+{\mathrm{\&phi;}}^n\left(0,j\right)\&rsqb;+\frac{{\&lsqb;V_{P_1}\left(0,j\right)\&times;\mathrm{\&Delta;}t\&rsqb;}^2}{\left(1+\cos \mathrm{\&theta;}\right)\&lsqb;V_{P_1}\left(0,j\right)\&times;\mathrm{\&Delta;}t+h\&rsqb;\&times;h}\\ \&times;\&lsqb;{\mathrm{\&phi;}}^n\left(1,j+1\right)+{\mathrm{\&phi;}}^n\left(0,j+1\right)-2{\mathrm{\&phi;}}^n\left(1,j\right)-2{\mathrm{\&phi;}}^n\left(0,j\right)+{\mathrm{\&phi;}}^n\left(1,j-1\right)+{\mathrm{\&phi;}}^n\left(0,j-1\right)\&rsqb;\end{array}---\left(1.81\right)$

其它三个边界条件的格式可以类推，不赘述。

由公式(1.81)可以发现，用于计算左边界条件的差分格式中含有j+1和j-1项。这就是说，如果j的取值从0到N(正整数)，则在计算φⁿ⁺¹(0,0)时要用到φⁿ⁺¹(0,-1)和φⁿ⁺¹(0,N+1)的值，而这两个值属于计算空间外的值。另三个角点也有类似的问题。所以，对于矩形边界计算空间而言，上述四个边界条件不能用于角点的计算。因此需要特殊的差分格式计算角点场值。本发明采用Taflove和Umashankar(1982)给出的方法处理角点，以左上角(0,0)点为例讨论。

图4是吸收边界条件中角点处理方法的示意图。设波传播方向上距(0,0)点一个网格步长h处的外行波场值为则角点(0,0)处的值可视为沿射线传播h的结果。角点上场值φ与之间满足下列关系：

${\mathrm{\&phi;}}^{n+1}(0,0)=f_{rad}{\stackrel{\&OverBar;}{\mathrm{\&phi;}}}^{n-1}$

其中f_rad为外行散射波的衰减因子，且d_c为波从激发源点传播到取值处以空间步长h为单位的距离；是(0,0)点及在网格空间内部相邻点(0,1)、(1,0)、(1,1)处在n-1时间步场值的线性插值。若用α表示径向射线与x轴之间的夹角，则有公式(1.82)，如下所示：

$\begin{array}{c}{\stackrel{\&OverBar;}{\mathrm{\&phi;}}}^{n-1}=\left(1-\sin \mathrm{\&alpha;}\right)\left(1-\cos \mathrm{\&alpha;}\right){\mathrm{\&phi;}}^{n-1}\left(0,0\right)+\left(1-\sin \mathrm{\&alpha;}\right){\mathrm{cos\&alpha;\&phi;}}^{n-1}\left(1,0\right)\\ +\sin \mathrm{\&alpha;}\left(1-\cos \mathrm{\&alpha;}\right){\mathrm{\&phi;}}^{n-1}\left(0,1\right)+{\mathrm{sin\&alpha;cos\&alpha;\&phi;}}^{n-1}\left(1,1\right)\end{array}---\left(1.82\right)$

其它三个角点的场值计算格式可类似得到。

图5a和图5b分别为采用以上边界条件计算的均匀双相各向同性介质中固相波场和流相波场的快照示意图，t＝700ms，激发源设在网格中心，计算空间大小为3000m×3000m，空间步长为10m×10m，时间步长为0.5ms，主频60Hz，差分精度为8阶，介质弹性参数见表2。图6a和图6b分別为没有采用任何边界条件处理计算得到的均匀双相各向同性介质中的固相波场和流相波场的快照示意图，其它参数同图5a和图5b。

表2均匀双相各向同性介质的弹性参数

P	Q	R	ρ11	ρ22	ρ12	b
							32.025	2.2575	16.10	1050	1050	0	0

P,Q,R：10⁹kg.m^-1.s^-2；ρ_ij：kg.m^-3；b：10^-6kg.m^-3.s^-1

比较图5a与图5b和图6a与图6b可以看出图，5a与图5b中第一纵波P₁在人为边界的伪反射大部分被吸收，而图6a与图6b中第一纵波P₁在人为边界的伪反射很强，说明上述边界条件能很好地压制边界的伪反射。

3二维双相各向同性介质中的弹性波方程差分格式的构造

在公式(1.21)中，令对y的偏导数等于零，并进行整理，得到保守系中双相各向同性介质的二维弹性波方程，如公式(1.83a)、(1.83b)、(1.83c)和(1.83d)所示：

$\left(A+2N\right)\frac{{\&part;}^2{u}_x}{\&part;x^2}+N\frac{{\&part;}^2{u}_x}{\&part;z^2}+\left(A+N\right)\frac{{\&part;}^2{u}_z}{\&part;x\&part;z}+Q\left(\frac{{\&part;}^2{u}_x}{\&part;x^2}+\frac{{\&part;}^2{U}_x}{\&part;x\&part;z}\right)=\frac{{\&part;}^2}{\&part;t^2}\left({\mathrm{\&rho;}}_{11}{u}_x+{\mathrm{\&rho;}}_{12}{U}_x\right)---\left(1.83a\right)$

$Q(\frac{{\&part;}^2{u}_x}{\&part;x^2}+\frac{{\&part;}^2{u}_z}{\&part;x\&part;z})+R(\frac{{\&part;}^2{U}_x}{\&part;x^2}+\frac{{\&part;}^2{U}_z}{\&part;x\&part;z})=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{12}{u}_x+{\mathrm{\&rho;}}_{22}{U}_x)---\left(1.83b\right)$

$(A+2N)\frac{{\&part;}^2{u}_z}{\&part;z^2}+N\frac{{\&part;}^2{u}_z}{\&part;x^2}+(A+N)\frac{{\&part;}^2{u}_x}{\&part;x\&part;z}+Q(\frac{{\&part;}^2{U}_z}{\&part;z^2}+\frac{{\&part;}^2{U}_x}{\&part;x\&part;z})=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{11}{u}_z+{\mathrm{\&rho;}}_{12}{U}_z)---\left(1.83c\right)$

$Q(\frac{{\&part;}^2{u}_z}{\&part;z^2}+\frac{{\&part;}^2{u}_x}{\&part;x\&part;z})+R(\frac{{\&part;}^2{U}_z}{\&part;z^2}+\frac{{\&part;}^2{U}_x}{\&part;x\&part;z})=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{12}{u}_z+{\mathrm{\&rho;}}_{22}{U}_z)---\left(1.83d\right)$

同样，在公式(1.19)中，令对y的偏导数等于零，并进行整理，得到有耗散时双相各向同性介质的二维弹性波方程，如公式(1.84a)、(1.84b)、(1.84c)和(1.84d)所示：

$\begin{array}{c}\left(A+2N\right)\frac{{\&part;}^2{u}_x}{\&part;x^2}+N\frac{{\&part;}^2{u}_x}{\&part;z^2}+\left(A+N\right)\frac{{\&part;}^2{u}_x}{\&part;x\&part;z}+Q(\frac{{\&part;}^2{U}_x}{\&part;x^2}+\frac{{\&part;}^2{U}_z}{\&part;x\&part;z})\\ =\frac{{\&part;}^2}{\&part;t^2}\left({\mathrm{\&rho;}}_{11}{u}_x+{\mathrm{\&rho;}}_{12}{U}_x\right)+b\frac{\&part;}{\&part;t}\left(u_x-U_x\right)\end{array}---\left(1.84a\right)$

$Q(\frac{{\&part;}^2{u}_x}{\&part;x^2}+\frac{{\&part;}^2{u}_z}{\&part;x\&part;z})+R(\frac{{\&part;}^2{U}_x}{\&part;x^2}+\frac{{\&part;}^2{U}_z}{\&part;x\&part;z})=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{12}{u}_x+{\mathrm{\&rho;}}_{22}{U}_x)-b\frac{\&part;}{\&part;t}\left(u_x-U_x\right)---\left(1.84b\right)$

$\begin{array}{c}(A+2N)\frac{{\&part;}^2{u}_z}{\&part;z^2}+N\frac{{\&part;}^2{u}_z}{\&part;x^2}+(A+N)\frac{{\&part;}^2{u}_x}{\&part;x\&part;z}+Q(\frac{{\&part;}^2{U}_z}{\&part;z^2}+\frac{{\&part;}^2{U}_x}{\&part;x\&part;z})\\ =\frac{{\&part;}^2}{\&part;t^2}\left({\mathrm{\&rho;}}_{11}{u}_z+{\mathrm{\&rho;}}_{12}{U}_z\right)+b\frac{\&part;}{\&part;t}\left(u_z-U_z\right)\end{array}---\left(1.84c\right)$

$Q(\frac{{\&part;}^2{u}_z}{\&part;z^2}+\frac{{\&part;}^2{u}_x}{\&part;x\&part;z})+R(\frac{{\&part;}^2{U}_z}{\&part;z^2}+\frac{{\&part;}^2{U}_x}{\&part;x\&part;z})=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{12}{u}_z+{\mathrm{\&rho;}}_{22}{U}_z)-b\frac{\&part;}{\&part;t}(u_z-U_z)---\left(1.84d\right)$

3.1二维双相各向同性介质中的弹性波方程的高阶差分格式

对空间偏导数采用任意偶数2N₁(N₁为正整数)阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，公式(1.83a)、(1.83b)、(1.83c)和(1.83d)有如下的离散格式，如公式(1.85a)、(1.85b)、(1.85c)和(1.85d)所示：

$\begin{array}{c}\left(A+2N\right)\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i+m,j\right)}{{\mathrm{\&Delta;x}}^2}+N\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i,j+m\right)}{{\mathrm{\&Delta;z}}^2}\\ +\left(A+N\right)\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_z^n\left(i+m_1,j+m_2\right)}{\mathrm{\&Delta;}x\mathrm{\&Delta;}z}\\ +Q\left(\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_x^n\left(i+m,j\right)}{{\mathrm{\&Delta;x}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_z\left(i+m_1,j+m_2\right)}{\mathrm{\&Delta;}x\mathrm{\&Delta;}z}\right)\\ ={\mathrm{\&rho;}}_{11}\frac{u_x^{n+1}\left(i,j\right)+u_x^{n-1}\left(i,j\right)-2u_x^n\left(i,j\right)}{{\mathrm{\&Delta;t}}^2}+{\mathrm{\&rho;}}_{12}\frac{U_x^{n+1}\left(i,j\right)+u_x^{n-1}\left(i,j\right)-2U_x^n\left(i,j\right)}{{\mathrm{\&Delta;t}}^2}\end{array}---\left(1.85a\right)$

$\begin{array}{c}Q\left(\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i+m,j\right)}{{\mathrm{\&Delta;x}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_z^n\left(i+m_1,j+m_2\right)}{\mathrm{\&Delta;}x\mathrm{\&Delta;}z}\right)\\ +R\left(\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_x^n\left(i+m,j\right)}{{\mathrm{\&Delta;x}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_z\left(i+m_1,j+m_2\right)}{\mathrm{\&Delta;}x\mathrm{\&Delta;}z}\right)\\ ={\mathrm{\&rho;}}_{12}\frac{u_x^{n+1}\left(i,j\right)+u_x^{n-1}\left(i,j\right)-2u_x^n\left(i,j\right)}{{\mathrm{\&Delta;t}}^2}+{\mathrm{\&rho;}}_{22}\frac{U_x^{n+1}\left(i,j\right)+U_x^{n-1}\left(i,j\right)-2U_x^n\left(i,j\right)}{{\mathrm{\&Delta;t}}^2}\end{array}---\left(1.85b\right)$

$\begin{array}{c}\left(A+2N\right)\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i,j+m\right)}{{\mathrm{\&Delta;z}}^2}+N\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i+m,j\right)}{{\mathrm{\&Delta;x}}^2}\\ +\left(A+N\right)\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_x^n\left(i+m_1,j+m_2\right)}{\mathrm{\&Delta;}x\mathrm{\&Delta;}z}\\ +Q\left(\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_z^n\left(i,j+m\right)}{{\mathrm{\&Delta;z}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_x\left(i+m_1,j+m_2\right)}{\mathrm{\&Delta;}x\mathrm{\&Delta;}z}\right)\\ ={\mathrm{\&rho;}}_{11}\frac{u_z^{n+1}\left(i,j\right)+u_z^{n-1}\left(i,j\right)-2u_z^n\left(i,j\right)}{{\mathrm{\&Delta;t}}^2}+{\mathrm{\&rho;}}_{12}\frac{U_z^{n+1}\left(i,j\right)+u_z^{n-1}\left(i,j\right)-2U_z^n\left(i,j\right)}{{\mathrm{\&Delta;t}}^2}\end{array}---\left(1.85c\right)$

$\begin{array}{c}Q\left(\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i,j+m\right)}{{\mathrm{\&Delta;z}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_x^n\left(i+m_1,j+m_2\right)}{\mathrm{\&Delta;}x\mathrm{\&Delta;}z}\right)\\ +R\left(\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_z^n\left(i,j+m\right)}{{\mathrm{\&Delta;z}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_x\left(i+m_1,j+m_2\right)}{\mathrm{\&Delta;}x\mathrm{\&Delta;}z}\right)\\ ={\mathrm{\&rho;}}_{12}\frac{u_z^{n+1}\left(i,j\right)+u_z^{n-1}\left(i,j\right)-2u_z^n\left(i,j\right)}{{\mathrm{\&Delta;t}}^2}+{\mathrm{\&rho;}}_{22}\frac{U_z^{n+1}\left(i,j\right)+U_z^{n-1}\left(i,j\right)-2U_z^n\left(i,j\right)}{{\mathrm{\&Delta;t}}^2}\end{array}---\left(1.85d\right)$

其中，△x为x方向空间差分步长，△z为z方向空间差分步长；为空间二阶导数的偶数阶精度差分系数、和为空间一阶导数的偶数阶精度差分系数，且取△x＝△z＝△s，进一步化简和整理上式，得到公式(1.86a)、(1.86b)、(1.86c)和(1.86d)，如下所示：

$\begin{array}{c}{\mathrm{\&rho;}}_{11}{u}_x^{n+1}\left(i,j\right)+{\mathrm{\&rho;}}_{12}{U}_x^{n+1}\left(i,j\right)\\ =\left(A+2N\right)\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2+N\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i,j+m\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\\ +\left(A+N\right)\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_z^n\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\\ +Q\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_x^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_z\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}\&rsqb;{\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{11}\&lsqb;u_x^{n-1}\left(i,j\right)-2u_x^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{12}\&lsqb;U_x^{n-1}\left(i,j\right)-2U_x^n\left(i,j\right)\&rsqb;\end{array}---\left(1.86a\right)$

$\begin{array}{c}{\mathrm{\&rho;}}_{12}{u}_x^{n+1}\left(i,j\right)+{\mathrm{\&rho;}}_{22}{U}_x^{n+1}\left(i,j\right)\\ =Q\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_z^n\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\&rsqb;\\ +R\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_x^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_z\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}\&rsqb;{\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{12}\&lsqb;u_x^{n-1}\left(i,j\right)-2u_x^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{22}\&lsqb;U_x^{n-1}\left(i,j\right)-2U_x^n\left(i,j\right)\&rsqb;\end{array}---\left(1.86b\right)$

$\begin{array}{c}{\mathrm{\&rho;}}_{11}{u}_z^{n+1}\left(i,j\right)+{\mathrm{\&rho;}}_{12}{U}_z^{n+1}\left(i,j\right)\\ =\left(A+2N\right)\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2+N\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i,j+m\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\\ +\left(A+N\right)\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_x^n\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\\ +Q\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_z^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_x\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}\&rsqb;{\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{11}\&lsqb;u_z^{n-1}\left(i,j\right)-2u_z^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{12}\&lsqb;U_z^{n-1}\left(i,j\right)-2U_z^n\left(i,j\right)\&rsqb;\end{array}---\left(1.86c\right)$

$\begin{array}{c}{\mathrm{\&rho;}}_{12}{u}_z^{n+1}\left(i,j\right)+{\mathrm{\&rho;}}_{22}{U}_z^{n+1}\left(i,j\right)\\ =Q\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_x^n\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\&rsqb;\\ +R\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_z^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_x\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}\&rsqb;{\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{12}\&lsqb;u_z^{n-1}\left(i,j\right)-2u_z^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{22}\&lsqb;U_z^{n-1}\left(i,j\right)-2U_z^n\left(i,j\right)\&rsqb;\end{array}---\left(1.86d\right)$

解公式(1.86a)、(1.86b)、(1.86c)和(1.86d)，得公式(1.87)，如下所示：

$\left\{\begin{array}{c}{u}_x^{n+1}\left(i,j\right)=\left({\mathrm{\&rho;}}_{22}A-{\mathrm{\&rho;}}_{12}B\right)/\left({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2\right)\\ U_x^{n+1}\left(i,j\right)=\left({\mathrm{\&rho;}}_{11}B-{\mathrm{\&rho;}}_{12}A\right)/\left({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2\right)\\ u_z^{n+1}\left(i,j\right)=\left({\mathrm{\&rho;}}_{22}C-{\mathrm{\&rho;}}_{12}D\right)/\left({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2\right)\\ U_z^{n+1}\left(i,j\right)=\left({\mathrm{\&rho;}}_{11}D-{\mathrm{\&rho;}}_{12}C\right)/\left({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2\right)\end{array}\right.---\left(1.87\right)$

其中，A、B、C、D分别定义如公式(1.88a)、(1.88b)、(1.88c)和(1.88d)所示：

$\begin{array}{c}A=\left(A+2N\right)\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2+N\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i,j+m\right)}{{\mathrm{\&Delta;z}}^2}{\mathrm{\&Delta;t}}^2\\ +\left(A+N\right)\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_z^n\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\\ +Q\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_x^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_z\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}\&rsqb;{\mathrm{\&Delta;t}}^2\\ ={\mathrm{\&rho;}}_{11}\&lsqb;u_x^{n-1}\left(i,j\right)-2u_x^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{12}\&lsqb;U_x^{n-1}\left(i,j\right)-2U_x^n\left(i,j\right)\&rsqb;\end{array}---\left(1.88a\right)$

$\begin{array}{c}B=Q\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_x}^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_z^n\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\&rsqb;\\ +R\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_x^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_z\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}\&rsqb;{\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{12}\&lsqb;u_x^{n-1}\left(i,j\right)-2u_x^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{22}\&lsqb;U_x^{n-1}\left(i,j\right)-2U_x^n\left(i,j\right)\&rsqb;\end{array}---\left(1.88b\right)$

$\begin{array}{c}C=\left(A+2N\right)\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2+N\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i,j+m\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\\ +\left(A+N\right)\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_x^n\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\\ +Q\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_z^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_x\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}\&rsqb;{\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{11}\&lsqb;u_z^{n-1}\left(i,j\right)-2u_z^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{12}\&lsqb;U_z^{n-1}\left(i,j\right)-2U_z^n\left(i,j\right)\&rsqb;\end{array}---\left(1.88c\right)$

$\begin{array}{c}D=Q\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{u_z}^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{u}_x^n\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}{\mathrm{\&Delta;t}}^2\&rsqb;\\ +R\&lsqb;\underset{m=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N_1\right)}{U}_z^n\left(i+m,j\right)}{{\mathrm{\&Delta;s}}^2}+\underset{m_1=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\underset{m_2=-N_1}{\overset{N_1}{\mathrm{\&Sigma;}}}\frac{c_{m_1}^{\left(N_1\right)}{c}_{m_2}^{\left(N_1\right)}{U}_x\left(i+m_1,j+m_2\right)}{{\mathrm{\&Delta;s}}^2}\&rsqb;{\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{12}\&lsqb;u_z^{n-1}\left(i,j\right)-2u_z^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{22}\&lsqb;U_z^{n-1}\left(i,j\right)-2U_z^n\left(i,j\right)\&rsqb;\end{array}---\left(1.88d\right)$

这就是保守系中双相各向同性介质中弹性波方程的有限差分法差分格式。

对于有耗散时双相各向同性介质中的弹性波方程有限差分格式，其推导过程与保守系中的弹性波方程类似，本发明略去其推导过程，其差分格式如公式(1.89)所示：

$\left\{\begin{array}{c}{u}_x^{n+1}\left(i,j\right)=\frac{\left(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t\right)A^{\&prime;}-\left(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t\right)B^{\&prime;}}{\left(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t\right)\left(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t\right)-{\left(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t\right)}^2}\\ U_x^{n+1}\left(i,j\right)=\frac{\left(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t\right)B^{\&prime;}-\left(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t\right)A^{\&prime;}}{\left(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t\right)\left(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t\right)-{\left(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t\right)}^2}\\ u_z^{n+1}\left(i,j\right)=\frac{\left(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t\right)C^{\&prime;}-\left(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t\right)D^{\&prime;}}{\left(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t\right)\left(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t\right)-{\left(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t\right)}^2}\\ U_z^{n+1}\left(i,j\right)=\frac{\left(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t\right)D^{\&prime;}-\left(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t\right)C^{\&prime;}}{\left(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t\right)\left(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t\right)-{\left(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t\right)}^2}\end{array}\right.---\left(1.89\right)$

其中，A’、B’、C’、D’分别定义如公式(1.90a)、(1.90b)、(1.90c)和(1.90d)所示：

$A^{\&prime;}=2A+b\&lsqb;u_x^{n-1}(i,j)-U_x^{n-1}(i,j)\&rsqb;\mathrm{\&Delta;}t---\left(1.90a\right)$

$B^{\&prime;}=2B+b\&lsqb;u_x^{n-1}(i,j)-U_x^{n-1}(i,j)\&rsqb;\mathrm{\&Delta;}t---\left(1.90b\right)$

$C^{\&prime;}=2C+b\&lsqb;u_z^{n-1}(i,j)-U_z^{n-1}(i,j)\&rsqb;\mathrm{\&Delta;}t---\left(1.90c\right)$

$D^{\&prime;}=2D+b\&lsqb;u_z^{n-1}(i,j)-U_z^{n-1}(i,j)\&rsqb;\mathrm{\&Delta;}t---\left(1.90d\right)$

3.2吸收边界条件

采用与双相各向同性介质中纵波方程相同的加权方向校正吸收边界条件计算截断边界处的地震波场值，四个角点的计算也采用与双相各向同性介质中纵波方程角点计算相同的方法，数值计算表明此方法能够吸收大部分边界反射，达到较好的效果(图1-14)。

图7a和图7b分别是利用上述方法计算的双相各向同性介质中x分量波场的固相波场和流相波场的快照的示意图，图7c和图7d分别是利用上述方法计算的双相各向同性介质中z分量波场的固相波场和流相波场的快照的示意图，t＝400ms，介质弹性参数见表3。在双相各向同性介质中存在着两类纵波和一类横波，在图7a和图7b中，三种波均清晰可见，且其传播速度与理论分析结果一致，表明本发明的算法是正确的。

表3模型四弹性参数

A,N,Q,R：10⁹kg.m^-1.s^-2；ρ_ij：kg.m^-3；b：10^-6kg.m^-3.s^-1

3.3二维双相各向同性介质中的应力-速度弹性波方程

计算精度和运算效率是波场正演模拟的一个很重要方面，弹性波的正演比纵波的正演计算量大，因此其运算效率相对更低。为提高运算效率，将由应力表示的弹性波方程，如公式(1.15)，进行变形，化为一阶偏微分方程，即化为一阶应力-速度弹性波方程，其表达式如公式(1.91)所示：

$\left\{\begin{array}{c}\frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;x}+\frac{\&part;{\mathrm{\&tau;}}_{xy}}{\&part;y}+\frac{\&part;{\mathrm{\&tau;}}_{xz}}{\&part;z}=\frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{11}{v}_x+{\mathrm{\&rho;}}_{12}{V}_x\right)\\ \frac{\&part;s}{\&part;x}=\frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{12}{v}_x+{\mathrm{\&rho;}}_{22}{V}_x\right)\\ \frac{\&part;{\mathrm{\&tau;}}_{yx}}{\&part;x}+\frac{\&part;{\mathrm{\&sigma;}}_{yy}}{\&part;y}+\frac{\&part;{\mathrm{\&tau;}}_{yz}}{\&part;z}=\frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{11}{v}_y+{\mathrm{\&rho;}}_{12}{V}_y\right)\\ \frac{\&part;s}{\&part;y}=\frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{12}{v}_y+{\mathrm{\&rho;}}_{22}{V}_y\right)\\ \frac{\&part;{\mathrm{\&tau;}}_{zx}}{\&part;x}+\frac{\&part;{\mathrm{\&tau;}}_{zy}}{\&part;y}+\frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;z}=\frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{11}{v}_z+{\mathrm{\&rho;}}_{12}{V}_z\right)\\ \frac{\&part;s}{\&part;z}=\frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{12}{v}_z+{\mathrm{\&rho;}}_{22}{V}_z\right)\end{array}\right.---\left(1.91\right)$

其中，v_x、v_y、v_z为固相质点振动速度在x、y、z三个方向的分量(m/s)；V_x、V_y、V_z为流相质点振动速度在x、y、z三个方向的分量(m/s)。令对y的偏导数等于零，得到二维应力-速度方程，如公式(1.92)所示：

$\left\{\begin{array}{c}\frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;x}+\frac{\&part;{\mathrm{\&tau;}}_{xz}}{\&part;z}=\frac{\&part;}{\&part;t}({\mathrm{\&rho;}}_{11}{v}_x+{\mathrm{\&rho;}}_{12}{V}_x)\\ \frac{\&part;s}{\&part;x}=\frac{\&part;}{\&part;t}({\mathrm{\&rho;}}_{12}{v}_x+{\mathrm{\&rho;}}_{22}{V}_x)\\ \frac{\&part;{\mathrm{\&tau;}}_{zx}}{\&part;x}+\frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;z}=\frac{\&part;}{\&part;t}({\mathrm{\&rho;}}_{11}{v}_z+{\mathrm{\&rho;}}_{12}{V}_z)\\ \frac{\&part;s}{\&part;z}=\frac{\&part;}{\&part;t}({\mathrm{\&rho;}}_{12}{v}_z+{\mathrm{\&rho;}}_{22}{V}_z)\end{array}\right.---\left(1.92\right)$

在公式(1.92)中，由前述已知τ_xz＝τ_zx，所以共有8个未知量，由应力-应变关系，如公式(1.8)，可以再补充四个应力-速度关系，表达式如公式(1.93)所示：

$\left\{\begin{array}{c}\frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;t}=2N\frac{\&part;v_x}{\&part;x}+A\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\\ \frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;t}=2N\frac{\&part;v_z}{\&part;z}+A\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\\ \frac{\&part;{\mathrm{\&tau;}}_{xz}}{\&part;t}=N\left(\frac{\&part;v_x}{\&part;z}+\frac{\&part;v_z}{\&part;x}\right)\\ \frac{\&part;s}{\&part;t}=Q\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\end{array}\right.---\left(1.93\right)$

联系公式(1.92)和公式(1.93)，便得到保守系中用应力-速度表示的弹性波方程，公式(1.94)所示：

$\left\{\begin{array}{c}\frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{11}{v}_x+{\mathrm{\&rho;}}_{12}{V}_x\right)=\frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;x}+\frac{\&part;{\mathrm{\&tau;}}_{xz}}{\&part;z}\\ \frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{12}{v}_x+{\mathrm{\&rho;}}_{22}{V}_x\right)=\frac{\&part;s}{\&part;x}\\ \frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{11}{v}_z+{\mathrm{\&rho;}}_{12}{V}_z\right)=\frac{\&part;{\mathrm{\&tau;}}_{zx}}{\&part;x}+\frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;z}\\ \frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{12}{v}_z+{\mathrm{\&rho;}}_{22}{V}_z\right)=\frac{\&part;s}{\&part;z}\\ \frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;t}=2N\frac{\&part;v_x}{\&part;x}+A\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\\ \frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;t}=2N\frac{\&part;v_z}{\&part;z}+A\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\\ \frac{\&part;{\mathrm{\&tau;}}_{xz}}{\&part;t}=N\left(\frac{\&part;v_x}{\&part;z}+\frac{\&part;v_z}{\&part;x}\right)\\ \frac{\&part;s}{\&part;t}=Q\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\end{array}\right.---\left(1.94\right)$

同理，可以得到有耗散时用应力-速度表示的弹性波方程，公式(1.95)所示：

$\left\{\begin{array}{c}\frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{11}{v}_x+{\mathrm{\&rho;}}_{12}{V}_x\right)+b\left(v_x-V_x\right)=\frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;x}+\frac{\&part;{\mathrm{\&tau;}}_{xz}}{\&part;z}\\ \frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{12}{v}_x+{\mathrm{\&rho;}}_{22}{V}_x\right)-b\left(v_x-V_x\right)=\frac{\&part;s}{\&part;x}\\ \frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{11}{v}_z+{\mathrm{\&rho;}}_{12}{V}_z\right)+b\left(v_z-V_z\right)=\frac{\&part;{\mathrm{\&tau;}}_{zx}}{\&part;x}+\frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;z}\\ \frac{\&part;}{\&part;t}\left({\mathrm{\&rho;}}_{12}{v}_z+{\mathrm{\&rho;}}_{22}{V}_z\right)-b\left(v_z-V_z\right)=\frac{\&part;s}{\&part;z}\\ \frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;t}=2N\frac{\&part;v_x}{\&part;x}+A\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\\ \frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;t}=2N\frac{\&part;v_z}{\&part;z}+A\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\\ \frac{\&part;{\mathrm{\&tau;}}_{xz}}{\&part;t}=N\left(\frac{\&part;v_x}{\&part;z}+\frac{\&part;v_z}{\&part;x}\right)\\ \frac{\&part;s}{\&part;t}=Q\left(\frac{\&part;v_x}{\&part;x}+\frac{\&part;v_z}{\&part;z}\right)+Q\left(\frac{\&part;V_x}{\&part;x}+\frac{\&part;V_z}{\&part;z}\right)\end{array}\right.---\left(1.95\right)$

3.4二维双相各向同性介质中的弹性波方程的交错网格高阶差分格式

为提高计算精度，采用交错网格中的高阶差分格式进行计算，如图8所示的交错网格示意图。

时间导数采用2阶精度差分，空间导数采用2N(N为正整数)阶精度差分，保守系中的双相各向同性介质二维弹性波方程差分格式为公式(1.96a)、(1.96b)、(1.96c)、(1.96d)、(1.96e)、(1.96f)、(1.96g)和(1.96h)，如下所示：

$v_x^{n+\frac{1}{2}}(i,j)=({\mathrm{\&rho;}}_{22}A-{\mathrm{\&rho;}}_{12}B)/({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2)---\left(1.96a\right)$

$V_x^{n+\frac{1}{2}}(i,j)=({\mathrm{\&rho;}}_{11}B-{\mathrm{\&rho;}}_{12}A)/({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2)---\left(1.96b\right)$

$v_z^{n+\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2})=({\mathrm{\&rho;}}_{22}C-{\mathrm{\&rho;}}_{12}D)/({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2)---\left(1.96c\right)$

$V_z^{n+\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2})=({\mathrm{\&rho;}}_{11}D-{\mathrm{\&rho;}}_{12}C)/({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2)---\left(1.96d\right)$

$\begin{array}{c}{\mathrm{\&sigma;}}_{xx}^n\left(i+\frac{1}{2},j\right)={\mathrm{\&sigma;}}_{xx}^{n-1}\left(i+\frac{1}{2},j\right)\\ +\left(A+2N\right)\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_x^{n-\frac{1}{2}}\left(i+n,j\right)-v_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\&CenterDot;\mathrm{\&Delta;}t\\ +A\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\\ +Q(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_x^{n-\frac{1}{2}}\left(i+n,j\right)-V_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z})\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.96e\right)$

$\begin{array}{c}{\mathrm{\&sigma;}}_{zz}^n\left(i+\frac{1}{2},j\right)={\mathrm{\&sigma;}}_{zz}^{n-1}\left(i+\frac{1}{2},j\right)\\ +\left(A+2N\right)\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\\ +A\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_x^{n-\frac{1}{2}}\left(i+n,j\right)-v_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\&CenterDot;\mathrm{\&Delta;}t\\ +Q(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_x^{n-\frac{1}{2}}\left(i+n,j\right)-V_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z})\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.96f\right)$

$\begin{array}{c}{\mathrm{\&tau;}}_{xz}^n\left(i,j+\frac{1}{2}\right)={\mathrm{\&tau;}}_{xz}^{n-1}\left(i,j+\frac{1}{2}\right)\\ +N(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_x^{n-\frac{1}{2}}\left(i,j+n\right)-v_x^{n-\frac{1}{2}}\left(i,j+1-n\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2}\left(2n-1\right),j+\frac{1}{2}\right)-v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2}\left(2n-1\right),j+\frac{1}{2}\right)\&rsqb;}{\mathrm{\&Delta;}x})\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.96g\right)$

$\begin{array}{c}{s}^n\left(i+\frac{1}{2},j\right)=s^{n-1}\left(i+\frac{1}{2},j\right)\\ +Q(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_x^{n-\frac{1}{2}}\left(i+n,j\right)-v_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z})\&CenterDot;\mathrm{\&Delta;}t\\ +R(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_x^{n-\frac{1}{2}}\left(i+n,j\right)-V_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z})\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.96h\right)$

其中，A、B、C、D分别定义如公式(1.97a)、(1.97b)、(1.97c)和(1.97d)所示：

$\begin{array}{c}A={\mathrm{\&rho;}}_{11}{v}_x^{n-\frac{1}{2}}\left(i,j\right)+{\mathrm{\&rho;}}_{12}{V}_x^{n-\frac{1}{2}}\left(i,j\right)\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;{\mathrm{\&sigma;}}_{xx}^n\left(i+\frac{1}{2}\left(2n-1\right),j\right)-{\mathrm{\&sigma;}}_{xx}^n\left(i-\frac{1}{2}\left(2n-1\right),j\right)\&rsqb;}{\mathrm{\&Delta;}x}\&CenterDot;\mathrm{\&Delta;}t\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;{\mathrm{\&tau;}}_{xz}^n\left(i,j+\frac{1}{2}\left(2n-1\right)\right)-{\mathrm{\&tau;}}_{xz}^n\left(i,j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.97a\right)$

$\begin{array}{c}B={\mathrm{\&rho;}}_{12}{v}_x^{n-\frac{1}{2}}\left(i,j\right)+{\mathrm{\&rho;}}_{22}{V}_x^{n-\frac{1}{2}}\left(i,j\right)\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;s^n\left(i+\frac{1}{2}\left(2n-1\right),j\right)-s^n\left(i-\frac{1}{2}\left(2n-1\right),j\right)\&rsqb;}{\mathrm{\&Delta;}x}\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.97b\right)$

$\begin{array}{c}C={\mathrm{\&rho;}}_{11}{v}_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\right)+{\mathrm{\&rho;}}_{12}{V}_x^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\right)\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;{\mathrm{\&tau;}}_{xz}^n\left(i+n,j+\frac{1}{2}\right)-{\mathrm{\&tau;}}_{xz}^n\left(i+1-n,j+\frac{1}{2}\right)\&rsqb;}{\mathrm{\&Delta;}x}\&CenterDot;\mathrm{\&Delta;}t\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;{\mathrm{\&sigma;}}_{zz}^n\left(i+\frac{1}{2},j+n\right)-{\mathrm{\&sigma;}}_{zz}^n\left(i+\frac{1}{2},j+1-n\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.97c\right)$

$\begin{array}{c}D={\mathrm{\&rho;}}_{12}{v}_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\right)+{\mathrm{\&rho;}}_{22}{V}_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\right)\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;s^n\left(i+\frac{1}{2},j+n\right)-s^n\left(i+\frac{1}{2},j+1-n\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.97d\right)$

类似地，有耗散时的双相各向同性介质中的二维弹性波方程差分格式，如公式(1.98a)、(1.98b)、(1.98c)、(1.98d)、(1.98e)、(1.98f)、(1.98g)和(1.98h)所示：

$v_x^{n+\frac{1}{2}}\left(i,j\right)=\frac{(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)A^{\&prime;}-(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)B^{\&prime;}}{(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t)(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)-{(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)}^2}---\left(1.98a\right)$

$V_x^{n+\frac{1}{2}}(i,j)=\frac{(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)B^{\&prime;}-(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)A^{\&prime;}}{(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t)(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)-{(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)}^2}---\left(1.98b\right)$

$v_z^{n+\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2})=\frac{(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)C^{\&prime;}-(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)D^{\&prime;}}{(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t)(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)-{(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)}^2}---\left(1.98c\right)$

$V_z^{n+\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2})=\frac{(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)D^{\&prime;}-(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)C^{\&prime;}}{(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t)(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)-{(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)}^2}---\left(1.98d\right)$

$\begin{array}{c}{\mathrm{\&sigma;}}_{xx}^n\left(i+\frac{1}{2},j\right)={\mathrm{\&sigma;}}_{xx}^{n-1}\left(i+\frac{1}{2},j\right)\\ +\left(A+2N\right)\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_x^{n-\frac{1}{2}}\left(i+n,j\right)-v_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\&CenterDot;\mathrm{\&Delta;}t\\ +A\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\\ +Q(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_x^{n-\frac{1}{2}}\left(i+n,j\right)-V_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z})\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.98e\right)$

$\begin{array}{c}{\mathrm{\&sigma;}}_{zz}^n\left(i+\frac{1}{2},j\right)={\mathrm{\&sigma;}}_{zz}^{n-1}\left(i+\frac{1}{2},j\right)\\ +\left(A+2N\right)\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\\ +A\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_x^{n-\frac{1}{2}}\left(i+n,j\right)-v_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\&CenterDot;\mathrm{\&Delta;}t\\ +Q(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_x^{n-\frac{1}{2}}\left(i+n,j\right)-V_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z})\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.98f\right)$

$\begin{array}{c}{\mathrm{\&tau;}}_{xz}^n\left(i,j+\frac{1}{2}\right)={\mathrm{\&tau;}}_{xz}^{n-1}\left(i,j+\frac{1}{2}\right)\\ +N(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_x^{n-\frac{1}{2}}\left(i,j+n\right)-v_x^{n-\frac{1}{2}}\left(i,j+1-n\right)\&rsqb;}{\mathrm{\&Delta;}z}\&CenterDot;\mathrm{\&Delta;}t\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2}\left(2n-1\right),j+\frac{1}{2}\right)-v_z^{n-\frac{1}{2}}\left(i-\frac{1}{2}\left(2n-1\right),j+\frac{1}{2}\right)\&rsqb;}{\mathrm{\&Delta;}x})\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.98g\right)$

$\begin{array}{c}{s}^n\left(i+\frac{1}{2},j\right)=s^{n-1}\left(i+\frac{1}{2},j\right)\\ +Q(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_x^{n-\frac{1}{2}}\left(i+n,j\right)-v_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-v_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z})\&CenterDot;\mathrm{\&Delta;}t\\ +R(\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_x^{n-\frac{1}{2}}\left(i+n,j\right)-V_x^{n-\frac{1}{2}}\left(i+1-n,j\right)\&rsqb;}{\mathrm{\&Delta;}x}\\ +\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_n^{\left(N\right)}\&lsqb;V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2}\left(2n-1\right)\right)-V_z^{n-\frac{1}{2}}\left(i+\frac{1}{2},j-\frac{1}{2}\left(2n-1\right)\right)\&rsqb;}{\mathrm{\&Delta;}z})\&CenterDot;\mathrm{\&Delta;}t\end{array}---\left(1.98h\right)$

其中，A’、B’、C’、D’分别定义如公式(1.99a)、(1.99b)、(1.99c)和(1.99d)所示：

$A^{\&prime;}=2A-b\&lsqb;v_x^{n-\frac{1}{2}}(i,j)-V_x^{n-\frac{1}{2}}(i,j)\&rsqb;\mathrm{\&Delta;}t---\left(1.99a\right)$

$B^{\&prime;}=2B+b\&lsqb;v_x^{n-\frac{1}{2}}(i,j)-V_x^{n-\frac{1}{2}}(i,j)\&rsqb;\mathrm{\&Delta;}t---\left(1.99b\right)$

$C^{\&prime;}=2C-b\&lsqb;v_z^{n-\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2})-V_z^{n-\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2})\&rsqb;\mathrm{\&Delta;}t---\left(1.99c\right)$

$D^{\&prime;}=2D+b\&lsqb;v_z^{n-\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2})-V_z^{n-\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2})\&rsqb;\mathrm{\&Delta;}t---\left(1.99d\right)$

在公式(1.96a)～(1.99d)中，△x、△z为空间差分步长(m)，且这里取△x＝△z；△t为时间差分步长(s)；i、j为空间离散点；n为时间离散点；为一阶空间导数的2N阶精度差分系数，可通过待定系数法求得，当2N＝4时，有 $C_1^{\left(2\right)}=1.125,C_2^{\left(2\right)}=-0.04166667;$ 当2N＝6时，有 $C_1^{\left(3\right)}=1.171875,C_2^{\left(3\right)}=-0.06510416,$ $C_3^{\left(3\right)}=0.0046875;$ 当2N＝8时，有 $C_1^{\left(4\right)}=1.196289,C_2^{\left(4\right)}=-0.0797526,$ $C_3^{\left(4\right)}=0.009570313,C_4^{\left(4\right)}=-0.0006975447;$ 当2N＝10时，有 $C_1^{\left(5\right)}=1.211243,$ $C_2^{\left(5\right)}=-0.08972168,C_3^{\left(5\right)}=0.001384277,C_4^{\left(5\right)}=-0.00176566,C_5^{\left(5\right)}=\mathrm{0.0001186795.}$

3.5吸收边界条件

吸收边界条件仍然采用加权方向校正吸收边界条件，四个角点的计算也采用与式(1.82)相同的方法，但是交错网格的特殊性使得有些项并不需要计算截断边界的吸收条件，这是因为它们不在网格的最外层，视网格的具体设计而定。数值计算表明此方法能够吸收大部分边界反射，达到较好的效果。

图9a和图9b分别是利用上述方法计算的双相各向同性介质中x分量波场的固相波场和流相波场的快照的示意图，图9c和图9d分别是利用上述方法计算的双相各向同性介质中z分量波场的固相波场和流相波场的快照的示意图，t＝400ms，所用差分参数与图9a、图9b、图9c和图9d完全相同，介质弹性参数也完全一致，见表3。对比图7a、图7b、图7c与图7d和图9a、图9b、图9c与图9d，发现在采用相同差分参数的条件下，采用交错网格计算的精度比规则网格高，从而频散现象得到了较好的压制；在具体的算法实现过程中，交错网格比规则网格节省大量内存；通过计算发现，交错网格算法比规则网格算法速度快，大大提高了运算效率。所以，在数值求解弹性波方程时，采用交错网格求解可以在保证精度的基础上，尽可能地采用较大空间步长，提高运算效率，是一种理想的数值模拟方法。

4波场模拟试验

为研究波在双相介质中的传播规律，本发明利用上述正演方法对几个双相-单相介质模型进行了正演计算。

模型一为两层介质，上层介质为单相介质，下层介质为双相介质，激发源设在单相介质中，介质弹性参数见表4。

表4模型一弹性参数

P,Q,R：10⁹kg.m^-1.s^-2；ρ_ij：kg.m^-3；b：10^-6kg.m^-3.s^-1

图10a和图10b分别是模型一的固相波场和流相波场的快照的示意图，t＝1200ms。在图10a和图10b中，P₁是第一纵波直达波，P₁₁是反射的第一纵波，P₂₁是透射的第一纵波，P₂₂是透射的第二纵波。从图10a和图10b中可见，在固相波场中存在第一纵波的反射和透射以及透射的第二纵波，在流相波场中存在透射的第一纵波和第二纵波。说明纵波由单相介质入射到单相-双相介质分界面上，在单相介质中将产生反射纵波，在双相介质中将产生透射的第一类纵波和第二类纵波。同时也说明第二类纵波只在双相介质中存在。

模型二为两层介质，上层介质为双相介质，下层介质为单相介质，激发源设在双相介质中，介质弹性参数见表5。

表5模型二弹性参数

P,Q,R：10⁹kg.m^-1.s^-2；ρ_ij：kg.m^-3；b：10^-6kg.m^-3.s^-1

图11a和图11b分别是模型二的固相波场和流相波场的快照的示意图，t＝2000ms。在图11a和图11b中，P₂是第二纵波直达波，P₁₁₁和P₁₁₂分别是第一类纵波反射的第一类纵波和第二类纵波，P₂₁₁是第一类纵波透射的第一类纵波；P₁₂₁和P₁₂₂分别是第二类纵波反射的第一类纵波和第二类纵波，P₂₂₁是第二类纵波透射的第一类纵波。从图11a和图11b中可见，在双相介质中纵波激发，将产生第一类纵波和第二类纵波(在图中第一类纵波P₁已传出计算空间外)，第一类纵波和第二类纵波在到达双相-单相介质分界面时，将分别产生波的反射和透射，第一类纵波和第二类纵波均在双相介质中产生反射的第一类纵波和第二类纵波以及在单相介质中产生透射的第一类纵波。进一步说明在双相介质中存在两类纵波——第一类纵波和第二类纵波，且第二类纵波不能在单相介质中传播。

模型三为两层介质，上层介质为双相介质，下层介质仍为双相介质，激发源设在上层介质中，如图12所示，介质弹性参数见表6。

表6模型三弹性参数

P,Q,R：10⁹kg.m^-1.s^-2；ρ_ij：kg.m^-3；b：10^-6kg.m^-3.s^-1

图13a和图13b分别是模型三的固相波场和流相波场的快照的示意图，t＝2500ms。在图13a和图13b中，P₁、P₂、P₁₁₁、P₁₁₂、P₂₁₁、P₁₂₁、P₁₂₂和P₂₂₁的含义同上文，P212和P222分别为第一类纵波和第二类纵波在双相介质2中透射的第二类纵波。图13a和图13b中，有些波由于能量太弱很难分辨出，但通过固相波场和流相波场的对比，仍可确定波的存在及位置。所以在双相介质中纵波源激发可产生两类纵波，这两类纵波在到达两种双相介质的分界面时，遵从斯奈尔定律，会发生波的反射和透射，分别产生反射的两类纵波和透射的两类纵波。

通过这三个理论模型分析可以得出，双相介质中由于流体的存在，固体和流体的相互作用，使双相介质中地震波的传播规律发生了变化，出现了第二类纵波，地震波场变得复杂。地下含油气介质是典型的双相介质，采用常规的单相介质理论有时不能很好地解决油气储层问题。因此，分析双相介质中地震波的特征，寻找基于双相介质理论的油气属性提取方法对于提高寻找油气层的精度和准确率具有重要意义。

5双相TTI介质中地震波方程数值模拟

双相介质中的固相和流相的应力-应变关系满足如公式(1.6)所示的广义Hooke定律，为方便，可以将公式(1.6)写为公式(1.100)，如下所示：

$\left[\begin{array}{c}{\mathrm{\&sigma;}}_{xx}\\ {\mathrm{\&sigma;}}_{yy}\\ {\mathrm{\&sigma;}}_{zz}\\ {\mathrm{\&sigma;}}_{yz}\\ {\mathrm{\&sigma;}}_{xz}\\ {\mathrm{\&sigma;}}_{xy}\\ s\end{array}\right]\left[\begin{array}{ccccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}& Q_1\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}& Q_2\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}& Q_3\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}& Q_4\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}& Q_5\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}& Q_6\\ Q_1& Q_2& Q_3& Q_4& Q_5& Q_6& R\end{array}\right]\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ e_{yz}\\ e_{xz}\\ e_{xy}\\ \mathrm{\&epsiv;}\end{array}\right]---\left(1.100\right)$

公式(1.100)简写为公式(1.101)，如下所示：

$\left[\begin{array}{c}\mathrm{\&sigma;}\\ s\end{array}\right]=\left[\begin{array}{cc}C& Q\\ Q^T& R\end{array}\right]\left[\begin{array}{c}e\\ \mathrm{\&epsiv;}\end{array}\right]---\left(1.101\right)$

式中σ为固相应力向量，σ_xx、σ_yy、σ_zz为固相正应力，σ_xy、σ_yz、σ_xz为固相切应力；s为流相的有效压力；C为固体骨架的弹性常数张量矩阵，R表示孔隙流体的弹性参数，Q表示固体体积与流体体积变化之间的耦合关系，Q＝(Q₁,Q₂,Q₃,Q₄,Q₅,Q₆)^T。那么，

对于双相各向同性介质，Q＝(Q,Q,Q,0,0,0)^T，可得公式(1.102)，如下所示：

$C=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{12}& 0& 0& 0\\ c_{12}& c_{11}& c_{12}& 0& 0& 0\\ c_{12}& c_{12}& c_{11}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{44}& 0\\ 0& 0& 0& 0& 0& c_{44}\end{array}\right]---\left(1.102\right)$

对于横向各向同性介质(VTI)介质，Q＝(Q₁,Q₁,Q₃,0,0,0)^T，可得公式(1.103)，如下所示：

$C=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& 0\\ c_{12}& c_{11}& c_{13}& 0& 0& 0\\ c_{13}& c_{13}& c_{33}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{44}& 0\\ 0& 0& 0& 0& 0& c_{66}\end{array}\right]---\left(1.103\right)$

对于方位各向异性(HTI)介质，Q＝(Q₁,Q₂,Q₁,0,0,0)^T，可得公式(1.104)，如下所示：

$C=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& 0\\ c_{12}& c_{22}& c_{12}& 0& 0& 0\\ c_{13}& c_{12}& c_{11}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{55}& 0\\ 0& 0& 0& 0& 0& c_{44}\end{array}\right]---\left(1.104\right)$

对于正交各向异性介质，Q＝(Q₁,Q₂,Q₃,0,0,0)^T，可得公式(1.105)，如下所示：

$C=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& 0\\ c_{12}& c_{22}& c_{23}& 0& 0& 0\\ C_{13}& c_{23}& c_{33}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{55}& 0\\ 0& 0& 0& 0& 0& c_{66}\end{array}\right]---\left(1.105\right)$

对于任意倾斜横向各向同性(TTI)介质，(Q₁,Q₂,Q₃,Q₄,Q₅,Q₆)^T可得公式(1.106)，如下所示：

$C=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right]---\left(1.106\right)$

5.1双相各向异性介质中的一阶应力-速度弹性波方程

双相各向异性介质中的一阶应力-速度弹性波方程，如公式(1.107)和(1.108)所示：

$({\mathrm{\&rho;}}_{12}^2-{\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22})\frac{\&part;v}{\&part;t}=({\mathrm{\&rho;}}_{12}+{\mathrm{\&rho;}}_{22})Bv-({\mathrm{\&rho;}}_{12}+{\mathrm{\&rho;}}_{22})BV+{\mathrm{\&rho;}}_{12}\&dtri;s-{\mathrm{\&rho;}}_{22}\&dtri;\&CenterDot;\mathrm{\&sigma;}---\left(1.107\right)$

$({\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2)\frac{\&part;V}{\&part;t}=({\mathrm{\&rho;}}_{11}+{\mathrm{\&rho;}}_{12})Bv-({\mathrm{\&rho;}}_{11}+{\mathrm{\&rho;}}_{12})BV+{\mathrm{\&rho;}}_{11}\&dtri;s-{\mathrm{\&rho;}}_{12}\&dtri;\&CenterDot;\mathrm{\&sigma;}---\left(1.108\right)$

对于固相，其应力-速度关系如公式(1.109)所示：

$\left[\begin{array}{c}\frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;t}\\ \frac{\&part;{\mathrm{\&sigma;}}_{yy}}{\&part;t}\\ \frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;t}\\ \frac{\&part;{\mathrm{\&sigma;}}_{yz}}{\&part;t}\\ \frac{\&part;{\mathrm{\&sigma;}}_{xz}}{\&part;t}\\ \frac{\&part;{\mathrm{\&sigma;}}_{xy}}{\&part;t}\end{array}\right]\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& C_{14}& C_{15}& C_{16}\\ C_{12}& C_{22}& C_{23}& C_{24}& C_{25}& C_{26}\\ C_{13}& C_{23}& C_{33}& C_{34}& C_{35}& C_{36}\\ C_{14}& C_{24}& C_{34}& C_{44}& C_{45}& C_{46}\\ C_{15}& C_{25}& C_{35}& C_{45}& C_{55}& C_{56}\\ C_{16}& C_{26}& C_{36}& C_{46}& C_{56}& C_{66}\end{array}\right]\left[\begin{array}{c}\frac{\&part;v_x}{\&part;x}\\ \frac{\&part;v_y}{\&part;y}\\ \frac{\&part;v_z}{\&part;z}\\ \frac{\&part;v_y}{\&part;z}+\frac{\&part;v_z}{\&part;y}\\ \frac{\&part;v_z}{\&part;x}+\frac{\&part;v_x}{\&part;z}\\ \frac{\&part;v_x}{\&part;y}+\frac{\&part;v_y}{\&part;x}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\&epsiv;Q}}_1\\ {\mathrm{\&epsiv;Q}}_2\\ {\mathrm{\&epsiv;Q}}_3\\ {\mathrm{\&epsiv;Q}}_4\\ {\mathrm{\&epsiv;Q}}_5\\ {\mathrm{\&epsiv;Q}}_6\end{array}\right]---\left(1.109\right)$

对于流相，则如公式(1.110)所示：

$\frac{\&part;s}{\&part;t}=Q_1\frac{\&part;v_x}{\&part;x}+Q_2\frac{\&part;v_y}{\&part;y}+Q_3\frac{\&part;v_z}{\&part;z}+Q_4(\frac{\&part;v_y}{\&part;z}+\frac{\&part;v_z}{\&part;y})+Q_5(\frac{\&part;v_z}{\&part;x}+\frac{\&part;v_x}{\&part;z})+Q_6(\frac{\&part;v_x}{\&part;y}+\frac{\&part;v_y}{\&part;x})+R\mathrm{\&epsiv;}---\left(1.110\right)$

其中， $\mathrm{\&epsiv;}=\frac{\&part;v_X}{\&part;x}+\frac{\&part;v_y}{\&part;y}+\frac{\&part;v_Z}{\&part;z}.$

5.2双相TTI介质中一阶应力-速度方程差分格式

同样，采用如图14所示的交错网格，相应的弹性波场分量和弹性参数空间位置见表7。

推导的双相TTI介质的2L阶空间差分精度、二阶时间差分精度的高阶有限差分格式，如公式(1.111)、(1.112)、(1.113)、(1.114)和(1.115)所示：

$\begin{array}{c}{v}_x^{n+1}\left(i^{+},j,k\right)=v_x^n\left(i^{+},j,k\right)\left\{1+\&lsqb;D_2\left(i^{+},j,k\right)+D_3\left(i^{+},j,k\right)\&rsqb;b_{11}\left(i^{+},j,k\right)\mathrm{\&Delta;}t\right\}\\ -V_x\left(i^{+},j,k\right)\&times;\&lsqb;D_2\left(i^{+},j,k\right)+D_3\left(i^{+},j,k\right)\&rsqb;b_{11}\left(i^{+},j,k\right)\mathrm{\&Delta;}t\\ +v_y^n\left(i,j^{+},k\right)\&lsqb;D_2\left(i,j^{+},k\right)+D_3\left(i,j^{+},k\right)\&rsqb;b_{12}\left(i,j^{+},k\right)\mathrm{\&Delta;}t\\ -V_y^n\left(i,j^{+},k\right)\&lsqb;D_2\left(i,j^{+},k\right)+D_3\left(i,j^{+},k\right)\&rsqb;b_{12}\left(i,j^{+},k\right)\mathrm{\&Delta;}t\\ +v_z^n\left(i,j,k^{+}\right)\&lsqb;D_2\left(i,j,k^{+}\right)+D_3\left(i,j,k^{+}\right)\&rsqb;b_{13}\left(i,j,k^{+}\right)\mathrm{\&Delta;}t\\ -V_z^n\left(i,j,k^{+}\right)\&lsqb;D_2\left(i,j,k^{+}\right)+D_3\left(i,j,k^{+}\right)\&rsqb;b_{13}\left(i,j,k^{+}\right)\mathrm{\&Delta;}t\\ +D_3\left(i^{+},j,k\right)L_x^{+}\&lsqb;s\left(i,j,k\right)\&rsqb;\mathrm{\&Delta;}t-D_2\left(i^{+},j,k\right)\\ \&times;\left\{L_x^{+}\&lsqb;{\mathrm{\&sigma;}}_{xx}^{n+}\left(i,j,k\right)\&rsqb;+L_y^{-}\&lsqb;{\mathrm{\&sigma;}}_{xy}^{n+}\left(i^{+},j^{+},k\right)\&rsqb;+L_z^{-}\&lsqb;{\mathrm{\&sigma;}}_{xz}^{n+}\left(i^{+},j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +f_x\left(i^{+},j,k\right)\end{array}---\left(1.111\right)$

$\begin{array}{c}{V}_x^{n+1}\left(i^{+},j,k\right)=V_x^n\left(i^{+},j,k\right)\left\{1+\&lsqb;D_1\left(i^{+},j,k\right)+D_3\left(i^{+},j,k\right)\&rsqb;b_{11}\left(i^{+},j,k\right)\mathrm{\&Delta;}t\right\}\\ -v_x\left(i^{+},j,k\right)\&times;\&lsqb;D_1\left(i^{+},j,k\right)+D_3\left(i^{+},j,k\right)\&rsqb;b_{11}\left(i^{+},j,k\right)\mathrm{\&Delta;}t\\ +V_y^n\left(i,j^{+},k\right)\&lsqb;D_1\left(i,j^{+},k\right)+D_3\left(i,j^{+},k\right)\&rsqb;b_{12}\left(i,j^{+},k\right)\mathrm{\&Delta;}t\\ -v_y^n\left(i,j^{+},k\right)\&lsqb;D_1\left(i,j^{+},k\right)+D_3\left(i,j^{+},k\right)\&rsqb;b_{12}\left(i,j^{+},k\right)\mathrm{\&Delta;}t\\ +V_z^n\left(i,j,k^{+}\right)\&lsqb;D_1\left(i,j,k^{+}\right)+D_3\left(i,j,k^{+}\right)\&rsqb;b_{13}\left(i,j,k^{+}\right)\mathrm{\&Delta;}t\\ -v_z^n\left(i,j,k^{+}\right)\&lsqb;D_1\left(i,j,k^{+}\right)+D_3\left(i,j,k^{+}\right)\&rsqb;b_{13}\left(i,j,k^{+}\right)\mathrm{\&Delta;}t\\ -D_1\left(i^{+},j,k\right)L_x^{+}\&lsqb;s\left(i,j,k\right)\&rsqb;\mathrm{\&Delta;}t+D_3\left(i^{+},j,k\right)\\ \&times;\left\{L_x^{+}\&lsqb;{\mathrm{\&sigma;}}_{xx}^{n+}\left(i,j,k\right)\&rsqb;+L_y^{-}\&lsqb;{\mathrm{\&sigma;}}_{xy}^{n+}\left(i^{+},j^{+},k\right)\&rsqb;+L_z^{-}\&lsqb;{\mathrm{\&sigma;}}_{xz}^{n+}\left(i^{+},j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t+f_x^f\left(i^{+},j,k\right)\end{array}---\left(1.112\right)$

$\begin{array}{c}{\mathrm{\&sigma;}}_{xx}^{n+}\left(i,j,k\right)={\mathrm{\&sigma;}}_{xx}^{n-}\left(i,j,k\right)+C_{11}\left(i,j,k\right)L_x^{-}\&lsqb;v_x^n\left(i^{+},j,k\right)\&rsqb;\mathrm{\&Delta;}t\\ +C_{12}\left(i,j,k\right)L_y^{-}\&lsqb;v_y^n\left(i,j^{+},k\right)\&rsqb;\mathrm{\&Delta;}t+C_{13}\left(i,j,k\right)L_z^{-}\&lsqb;v_z^n\left(i,j,k^{+}\right)\&rsqb;\mathrm{\&Delta;}t\\ +C_{14}\left(i,j,k\right)\left\{L_z^{-}\&lsqb;v_y^n\left(i,j^{+},k\right)\&rsqb;+L_y^{-}\&lsqb;v_z^n\left(i,j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +C_{15}\left(i,j,k\right)\left\{L_z^{-}\&lsqb;v_x^n\left(i^{+},j,k\right)\&rsqb;+L_x^{-}\&lsqb;v_z^n\left(i,j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +C_{16}\left(i,j,k\right)\left\{L_y^{-}\&lsqb;v_x^n\left(i^{+},j,k\right)\&rsqb;+L_x^{-}\&lsqb;v_y^n\left(i,j^{+},k\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +Q_1\left(i,j,k\right)\left\{L_x^{-}\&lsqb;V_x\left(i^{+},j,k\right)\&rsqb;+L_y^{-}\&lsqb;V_y\left(i,j^{+},k\right)\&rsqb;+L_z^{-}\&lsqb;V_z\left(i,j,k^{+}\right)\&rsqb;\right\}\end{array}---\left(1.113\right)$

$\begin{array}{c}{\mathrm{\&sigma;}}_{xz}^{n+}\left(i^{+},j,k^{+}\right)={\mathrm{\&sigma;}}_{xz}^{n-}\left(i^{+},j,k^{+}\right)+C_{15}\left(i^{+},j,k\right)L_x^{-}\&lsqb;v_x^n\left(i^{+},j,k\right)\&rsqb;\mathrm{\&Delta;}t\\ +C_{25}\left(i,j^{+},k\right)L_y^{-}\&lsqb;v_y^n\left(i,j^{+},k\right)\&rsqb;\mathrm{\&Delta;}t+C_{35}\left(i,j,k^{+}\right)L_z^{-}\&lsqb;v_z^n\left(i,j,k^{+}\right)\&rsqb;\mathrm{\&Delta;}t\\ +C_{45}\left(i,j^{+},k^{+}\right)\{L_z^{-}\&lsqb;v_y^n\left(i,j^{+},k\right)\&rsqb;+L_y^{-}\&lsqb;v_z^n\left(i,j,k^{+}\right)\&rsqb;\}\mathrm{\&Delta;}t\\ +C_{55}\left(i^{+},j,k^{+}\right)\left\{L_z^{+}\&lsqb;v_x^n\left(i^{+},j,k\right)\&rsqb;+L_x^{+}\&lsqb;v_z^n\left(i,j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +C_{56}\left(i^{+},j^{+},k\right)\left\{L_y^{-}\&lsqb;v_x^n\left(i^{+},j,k\right)\&rsqb;+L_x^{-}\&lsqb;v_y^n\left(i,j^{+},k\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +Q_5\left(i^{+},j,k^{+}\right)\left\{L_x^{-}\&lsqb;V_x\left(i^{+},j,k\right)\&rsqb;+L_y^{-}\&lsqb;V_y\left(i,j^{+},k\right)\&rsqb;+L_z^{-}\&lsqb;V_z\left(i,j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t\end{array}---\left(1.114\right)$

$\begin{array}{c}{S}^{n+}\left(i,j,k\right)=S^{n-}\left(i,j,k\right)+Q_1\left(i,j,k\right)L_x^{-}\&lsqb;v_x\left(i^{+},j,k\right)\&rsqb;\mathrm{\&Delta;}t\\ +Q_2\left(i,j,k\right)L_y^{-}\&lsqb;v_y\left(i,j^{+},k\right)\&rsqb;\mathrm{\&Delta;}t+Q_3\left(i,j,k\right)L_z^{-}\&lsqb;v_z\left(i,j,k^{+}\right)\&rsqb;\mathrm{\&Delta;}t\\ +Q_4\left(i,j^{+},k^{+}\right)\left\{L_z^{-}\&lsqb;v_y^n\left(i,j^{+},k\right)\&rsqb;+L_y^{-}\&lsqb;v_z^n\left(i,j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +Q_5\left(i^{+},j,k^{+}\right)\left\{L_z^{+}\&lsqb;v_x^n\left(i^{+},j,k\right)\&rsqb;+L_x^{+}\&lsqb;v_z^n\left(i,j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +Q_6\left(i^{+},j,k^{+}\right)\left\{L_y^{-}\&lsqb;v_x^n\left(i^{+},j,k\right)\&rsqb;+L_x^{-}\&lsqb;v_y^n\left(i,j^{+},k\right)\&rsqb;\right\}\mathrm{\&Delta;}t\\ +R(i,j,k)\left\{L_x^{-}\&lsqb;V_x\left(i^{+},j,k\right)\&rsqb;+L_y^{-}\&lsqb;V_y\left(i,j^{+},k\right)\&rsqb;+L_z^{-}\&lsqb;V_z\left(i,j,k^{+}\right)\&rsqb;\right\}\mathrm{\&Delta;}t\end{array}---\left(1.115\right)$

其中： $D_1=\frac{{\mathrm{\&rho;}}_{11}}{{{\mathrm{\&rho;}}_{12}}^2-{\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}},D_2=\frac{{\mathrm{\&rho;}}_{22}}{{{\mathrm{\&rho;}}_{12}}^2-{\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}},D_3=\frac{{\mathrm{\&rho;}}_{12}}{{{\mathrm{\&rho;}}_{12}}^2-{\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}}$

$L_{x_i}^{+}\&lsqb;u\left(x_i\right)\&rsqb;=\frac{1}{{\mathrm{\&Delta;x}}_i}\underset{l=1}{\overset{L}{\&Sigma;}}{a}_l\{u(x_i+{\mathrm{l\&Delta;x}}_i)-u\&lsqb;x_i-(l-1){\mathrm{\&Delta;x}}_i\&rsqb;\}$

$L_{x_i}^{-}\&lsqb;u\left(x_i\right)\&rsqb;=\frac{1}{{\mathrm{\&Delta;x}}_i}\underset{l=1}{\overset{L}{\&Sigma;}}{a}_l\{u(x_i+(l-1\left){\mathrm{\&Delta;x}}_i\right)-u\&lsqb;x_i-{\mathrm{l\&Delta;x}}_i\&rsqb;\}$

其余σ_xy、σ_zz、σ_yz、v_y、v_z、V_y和V_z的差分格式同理可得。

表7双相TTI介质弹性波场分量和弹性参数位置

5.3数值模拟算例

5.3.1均匀双相TTI介质模型

均匀双相TTI介质模型大小为500m×500m，模型参数见表8，正演所用网格大小△x＝△z＝5m，采样间隔△t＝0.5ms，完全匹配层的层数δ＝50，理论反射系数R＝0.001。震源为主频40Hz的Ricker子波，位于模型中心。

表8均匀双相TTI介质模型参数

注：C_ij,Q_i,R的单位为10⁹·kg·m^-1·s^-2,ρ_ij的单位为kg·m^-3；极化角＝45°，方位角＝60°

仍然采用二阶时间差分精度、十阶空间差分精度对上述模型做正演模拟，得到的波场快照，分别如图15a和图15b(300ms)、图16a和图16b(400ms)、图17a和图17b(600ms)所示。图15a、图15b、图16a和图16b中，出现了三个波，分别是准快P波、准慢P波和准横波，与单相TTI介质相比，多了一个准慢P波，这是因为在双相介质中由于固体和流体的相互作用而产生了慢纵波。在图16a和图16b中，准快P波已传到模型边界，边界吸收效果很好。图17a和图17b中，准快P波已传至模型外部。

5.3.2两层介质模型

上层为均匀单相各向同性介质，下层为均匀双相TTI介质，模型大小为500m×500m，两层介质分界面在模型中间，模型参数见表9，正演所用的网格大小△x＝△z＝5m，采样间隔△t＝0.5ms，完全匹配层的层数δ＝20，理论反射系数R＝0.001。震源为主频40Hz的Ricker子波，位于网格坐标的(150，250)处。

依然采用二阶时间差分精度、十阶空间差分精度对上述模型做正演模拟，得到的波场快照，分别如图18a和图18b(250ms)、图19a和图19b(300ms)、图20a和图20b(400ms)、图21a和图21b(600ms)所示，图22和图23分别是x分量和z分量地震记录。

表1-9两层介质模型参数

注：C_ij,Q_i,R的单位为10⁹·kg·m^-1·s^-2,ρ_ij的单位为kg·m^-3；极化角＝45°，方位角＝60°

综合分析图18a～图23可以得知，在单相介质中纵波源激发，当纵波传播到单相各向同性介质和双相TTI介质的分界面时，发生反射和透射，反射出纵波和横波，透射出准快纵波、准慢纵波和准横波。

5.3.3三维均匀单相TTI介质模型

三维均匀单相TTI介质模型大小为200m×200×200m，模型参数见表10，正演所用参数为：网格大小△x＝△y＝△z＝2m，采样间隔△t＝0.5ms，完全匹配层的层数δ＝50，理论反射系数R＝0.001。震源为主频35Hz的Ricker子波，位于模型中心。

表10三维均匀单相TTI介质模型参数

注：C_ij,Q_i,R的单位为10⁹·kg·m^-1·s^-2,ρ_ij的单位为kg·m^-3；极化角＝30°，方位角＝45°

依然采用二阶时间差分精度、十阶空间差分精度对上述模型做正演模拟，得到的波场快照，如图24(350ms)和图25(380ms)所示。

综合分析图24和图25，图24中，P波传到了边界面上，图25中S波未传到边界面上。从图24和图25中，可以看出在xoy面、xoz面和yoz面上不论是P波还是S波，其波前面均不是圆，在三个面上均明显地看到了横波分裂。

以上，大略说明了本发明之实施例所需要运用到的相关公式，以下将提供对应的实施例来进行说明。根据本发明的实施例，提供了一种基于有限差分法的三维TTI双相介质地震波场数值模拟方法。

图26是根据本发明实施例的基于有限差分法的三维TTI双相介质地震波场数值模拟方法的流程图。

步骤S2602，取得地震波的固体应力张量、流体应力张量、固体应变张量和流体应变张量。其中，所述固体应力张量如公式(1.1)所示，所述流体应力张量如公式(1.2)所示，所述固体应变张量如公式(1.4)所示，流体应变张量如公式(1.5)所示。

步骤S2604，根据应力与应变的对应关系，将所述固体应力张量、流体应力张量、固体应变张量和流体应变张量转换为所述地震波的本构方程式。其中，所述本构方程式根据应力与应变的对应关系而推得如公式(1.6)所示，并且当介质为各向同性介质时，进一步推得如公式(1.7)所示。

步骤S2606，根据应力与位移的对应关系，取得所述地震波的几何方程式。其中，几何方程式如公式(1.9)所示。

步骤S2608，根据所述本构方程式、所述几何方程式、流体相对于固体的运动及应力与位移的对应关系，取得所述地震波的运动微分方程式。其中，所述流体相对于固体的运动如公式(1.10)所示，所述应力与位移的对应关系如公式(1.11)所示，接着将公式(1.7)和公式(1.9)代入公式(1.10)和公式(1.11)，略去外力并整理，可得到如公式(1.18a)～(1.18f)所示的运动微分方程式，且所述运动微分方程式进一步转换为向量形式，则如公式(1.19)a和(1.19b)所示。

步骤S2610，对运动微分方程两边取散度，取得所述地震波的第一纵波方程式，并令第一纵波方程式中的耗散系数等于零，以取得第二纵波方程式。其中，第一纵波方程式如公式(1.22)所示，而第二纵波方程式如公式(1.23)所示。

步骤S2612，对第一纵波方程式，令对y的偏导数等于零，对空间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，取得第一差分方程式，其中N为大于1的正整数。其中，对第一纵波方程式，令y的偏导数等于零，可获得公式(1.64)和(1.65)，再对公式(1.65)采用对空间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，以取得如公式(1.66)所示的第一差分方程式。

步骤S2614，对第二纵波方程式，令对y的偏导数等于零，对空间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，取得第二差分方程式。其中，对第二纵波方程式，令y的偏导数等于零，可获得公式(1.57)和(1.58)，再对公式(1.58)采用对空间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，以取得如公式(1.61)所示的第二差分方程式。

步骤S2616，对所述第一差分方程式与所述第二差分方程式进行吸收边界条件处理，以取得对应的地震波场数值。其中，例如利用前述2.2节所述的吸收边界条件进行处理，以取得对应的地震波场数值，而对应的模拟结果可如图5a、图5b、图6a与图6b所示。

综上所述，根据本发明的技术方案，采用通过矩形体剖分、离散化，在时间上进行高阶近似，在边界条件上使用交错网格的吸收边界条件，实现了固体相与流体相耦合作用下的双相介质数值方程的迭代求解，即实现物理地震波场的实时传播模拟。本发明不仅可以可以模拟三维，还可以模拟二维介质模型。另外，本发明不仅可以输出模拟实际野外地震数据采集的单炮记录，还可以输出波场的时间切片，从而可以进行地震波场的三维重建。

以上所述仅为本发明的实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (9)

1.一种基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，包括：

取得地震波的固体应力张量、流体应力张量、固体应变张量和流体应变张量；

根据应力与应变的对应关系，将所述固体应力张量、流体应力张量、固体应变张量和流体应变张量转换为所述地震波的本构方程式；

根据应力与位移的对应关系，取得所述地震波的几何方程式；

根据所述本构方程式、所述几何方程式、流体相对于固体的运动及应力与位移的对应关系，取得所述地震波的运动微分方程式；

对运动微分方程两边取散度，取得所述地震波的第一纵波方程式，并令第一纵波方程式中的耗散系数等于零，以取得第二纵波方程式；

对第一纵波方程式，令对y的偏导数等于零，对空间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，取得第一差分方程式，其中N为大于1的正整数；

对第二纵波方程式，令对y的偏导数等于零，对空间偏导数采用2N阶精度展开式进行差分离散，对时间偏导数采用二阶精度中心差分格式进行差分离散，取得第二差分方程式；

对所述第一差分方程式与所述第二差分方程式进行吸收边界条件处理，以取得对应的地震波场数值。

2.根据权利要求1所述的基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，所述固体应力张量满足如下公式：

$\left[\begin{array}{ccc}{\mathrm{\&sigma;}}_{xx}& {\mathrm{\&tau;}}_{xy}& {\mathrm{\&tau;}}_{xz}\\ {\mathrm{\&tau;}}_{yx}& {\mathrm{\&sigma;}}_{yy}& {\mathrm{\&tau;}}_{yz}\\ {\mathrm{\&tau;}}_{zx}& {\mathrm{\&tau;}}_{zy}& {\mathrm{\&sigma;}}_{zz}\end{array}\right],$

其中，σ_xx、σ_yy、σ_zz为固体相的正应力，τ_xy、τ_xz、τ_yx、τ_yz、τ_zx、τ_zy为固体相的切应力，且τ_xy＝τ_yx，τ_xz＝τ_zx，τ_yz＝τ_zy；

所述流体应力张量满足如下公式：

$\left[\begin{array}{ccc}s& 0& 0\\ 0& s& 0\\ 0& 0& s\end{array}\right],s=-\mathrm{\&phi;}p,$

其中，φ为每单位截面中流体面积所占百分数，p为流体压力，负号表示应力s与流体压力p方向相反；

所述固体应变张量满足如下公式：

$\left[\begin{array}{ccc}{e}_{xx}& e_{xy}& e_{xz}\\ e_{yx}& e_{yy}& e_{yz}\\ e_{zx}& e_{zy}& e_{zz}\end{array}\right],$

其中，e_xx、e_yy、e_zz表示固相正应变，e_xy、e_xz、e_yx、e_yz、e_zx、e_zy表示固相切应变，且e_xy＝e_yx，e_xz＝e_zx，e_yz＝e_zy；

所述流体应变张量满足如下公式：

$\left[\begin{array}{ccc}\mathrm{\&epsiv;}& 0& 0\\ 0& \mathrm{\&epsiv;}& 0\\ 0& 0& \mathrm{\&epsiv;}\end{array}\right],$

其中，ε表示流相的体应变。

3.根据权利要求2所述的基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，所述本构方程式满足如下公式：

$\left[\begin{array}{c}{\mathrm{\&sigma;}}_{xx}\\ {\mathrm{\&sigma;}}_{yy}\\ {\mathrm{\&sigma;}}_{zz}\\ {\mathrm{\&tau;}}_{yz}\\ {\mathrm{\&tau;}}_{zx}\\ {\mathrm{\&tau;}}_{xy}\\ s\end{array}\right]=\left[\begin{array}{ccccccc}{C}_{1111}& C_{1122}& C_{1133}& 0& 0& 0& Q\\ C_{2211}& C_{2222}& C_{2233}& 0& 0& 0& Q\\ C_{3311}& C_{3322}& C_{3333}& 0& 0& 0& Q\\ 0& 0& 0& C_{2323}& 0& 0& 0\\ 0& 0& 0& 0& C_{3131}& 0& 0\\ 0& 0& 0& 0& 0& C_{1212}& 0\\ Q& Q& Q& 0& 0& 0& R\end{array}\right]\left[\begin{array}{c}{e}_{xx}\\ e_{yy}\\ e_{zz}\\ 2e_{yz}\\ 2e_{zx}\\ 2e_{xy}\\ \mathrm{\&epsiv;}\end{array}\right],$

其中，C_ijkl(i、j、k、l＝1，2，3，4)是固体相的弹性参数，R是流体相的弹性参数，Q＝(Q₁，Q2，Q3，Q4，Q5，Q6)^T是固相和流相耦合关系的弹性参数，C₁₁₂₂＝C₂₂₁₁＝C₁₁₃₃＝C₃₃₁₁＝C₂₂₃₃＝C₃₃₂₂＝A，C₁₂₁₂＝C₂₃₂₃＝C₃₁₃₁＝N，C₁₁₁₁＝C₂₂₂₂＝C₃₃₃₃＝A+2N，

其中， $\left\{\begin{array}{c}{\mathrm{\&sigma;}}_{xx}=2N{e}_{xx}+A\mathrm{\&theta;}+Q\mathrm{\&epsiv;}\\ {\mathrm{\&sigma;}}_{yy}=2N{e}_{yy}+A\mathrm{\&theta;}+Q\mathrm{\&epsiv;}\\ {\mathrm{\&sigma;}}_{zz}=2N{e}_{zz}+A\mathrm{\&theta;}+Q\mathrm{\&epsiv;}\\ {\mathrm{\&tau;}}_{xy}=2N{e}_{xy}\\ {\mathrm{\&tau;}}_{xz}=2N{e}_{xz}\\ {\mathrm{\&tau;}}_{yz}=2{\mathrm{Ne}}_{yz}\\ s=Q\mathrm{\&theta;}+R\mathrm{\&epsiv;}\end{array}\right.,$

θ表示固相体应变， $\mathrm{\&theta;}=divu=\frac{\&part;u_x}{\&part;x}+\frac{\&part;u_y}{\&part;y}+\frac{\&part;u_z}{\&part;z}=e_{xx}+e_{yy}+e_{zz};$ u为固相位移向量，u_x、u_y、u_z分别为固相位移向量u在x、y和z方向的分量；ε表示流相体应变，U为流相位移向量，U_x、U_y、U_z分别表示流相位移向量U在x、y和z方向的分量；A和N相当于单相各向同性弹性波理论中的拉梅系数，其中N＝μ；R表示使一定体积的流体流入某集合体而又使该集合体保持总体积不变所需施加在流体上的压力的一种量度；Q反映固体与流体体积变化之间的耦合性质。

4.根据权利要求3所述的基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，所述几何方程式满足如下公式：

$\left[\begin{array}{c}\begin{array}{ccc}{e}_{xx}=\frac{\&part;u_x}{\&part;x}& e_{yy}=\frac{\&part;u_y}{\&part;y}& e_{zz}=\frac{\&part;u_z}{\&part;z}\\ e_{xy}=\frac{1}{2}\left(\frac{\&part;u_x}{\&part;y}+\frac{\&part;u_y}{\&part;x}\right)& e_{xz}=\frac{1}{2}\left(\frac{\&part;u_x}{\&part;z}+\frac{\&part;u_z}{\&part;x}\right)& e_{yz}=\frac{1}{2}\left(\frac{\&part;u_y}{\&part;z}+\frac{\&part;u_z}{\&part;y}\right)\\ e_{xy}=e_{yx}& e_{xz}=e_{zx}& e_{yz}=e_{zy}\end{array}\\ \begin{array}{cc}e=\frac{\&part;u_x}{\&part;x}+\frac{\&part;u_y}{\&part;y}+\frac{\&part;u_z}{\&part;z}& \mathrm{\&epsiv;}=\frac{\&part;U_x}{\&part;x}+\frac{\&part;U_y}{\&part;y}+\frac{\&part;U_z}{\&part;z}\end{array}\end{array}\right].$

5.根据权利要求4所述的基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，所述流体相对于固体的运动满足如下公式：

$\left[\begin{array}{c}\frac{\&part;s}{\&part;x}+F_x^f\\ \frac{\&part;s}{\&part;y}+F_y^f\\ \frac{\&part;s}{\&part;z}+F_z^f\end{array}\right]=\left[\begin{array}{c}\frac{{\&part;}^2\left({\mathrm{\&rho;}}_{12}{u}_x+{\mathrm{\&rho;}}_{22}{U}_x\right)}{\&part;t^2}\\ \frac{{\&part;}^2\left({\mathrm{\&rho;}}_{12}{u}_y+{\mathrm{\&rho;}}_{22}{U}_y\right)}{\&part;t^2}\\ \frac{{\&part;}^2\left({\mathrm{\&rho;}}_{12}{u}_z+{\mathrm{\&rho;}}_{22}{U}_z\right)}{\&part;t^2}\end{array}\right]+\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]\left[\begin{array}{c}\frac{\&part;\left(U_x-u_x\right)}{\&part;t}\\ \frac{\&part;\left(U_y-u_y\right)}{\&part;t}\\ \frac{\&part;\left(U_z-u_z\right)}{\&part;t}\end{array}\right];$

所述应力与位移的对应关系满足如下公式：

$\left[\begin{array}{c}\frac{\&part;{\mathrm{\&sigma;}}_{xx}}{\&part;x}+\frac{\&part;{\mathrm{\&tau;}}_{xy}}{\&part;y}+\frac{\&part;{\mathrm{\&tau;}}_{xz}}{\&part;z}+F_x^s\\ \frac{\&part;{\mathrm{\&tau;}}_{yx}}{\&part;x}+\frac{\&part;{\mathrm{\&sigma;}}_{yy}}{\&part;y}+\frac{\&part;{\mathrm{\&tau;}}_{yz}}{\&part;z}+F_y^s\\ \frac{\&part;{\mathrm{\&tau;}}_{zx}}{\&part;x}+\frac{\&part;{\mathrm{\&tau;}}_{zy}}{\&part;y}+\frac{\&part;{\mathrm{\&sigma;}}_{zz}}{\&part;z}+F_z^s\end{array}\right]=\left[\begin{array}{c}\frac{{\&part;}^2\left({\mathrm{\&rho;}}_{11}{u}_x+{\mathrm{\&rho;}}_{12}{U}_x\right)}{\&part;t^2}\\ \frac{{\&part;}^2\left({\mathrm{\&rho;}}_{11}{u}_y+{\mathrm{\&rho;}}_{12}{U}_y\right)}{\&part;t^2}\\ \frac{{\&part;}^2\left({\mathrm{\&rho;}}_{11}{u}_z+{\mathrm{\&rho;}}_{12}{U}_z\right)}{\&part;t^2}\end{array}\right]-\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]\left[\begin{array}{c}\frac{\&part;\left(U_x-u_x\right)}{\&part;t}\\ \frac{\&part;\left(U_y-u_y\right)}{\&part;t}\\ \frac{\&part;\left(U_z-u_z\right)}{\&part;t}\end{array}\right]$

其中，和分别表示固相外力的三个分量，和分别表示流相外力的三个分量，[b_ij]_3×3为耗散系数矩阵，b_ij＝b_ji，其值由达西渗透系数k_ij、流体粘滞系数η和孔隙率φ决定，且满足ρ₁₁、ρ₂₂和ρ₁₂是质量密度参数，其中ρ₁₁表示单元体中固体相对流体运动时固体部分总的等效质量，ρ₂₂表示单元体中流体相对固体运动时流体部分总的等效质量，ρ₁₂表示流体和固体之间的质量耦合系数。

6.根据权利要求5所述的基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，所述运动微分方程式满足如下公式：

$N{\&dtri;}^{2}u+\&dtri;\&lsqb;(A+N)\mathrm{\&theta;}+Q\mathrm{\&epsiv;}\&rsqb;=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{11}u+{\mathrm{\&rho;}}_{12}U)+b\frac{\&part;(u-U)}{\&part;t},$

$\&dtri;(Q\mathrm{\&theta;}+R\mathrm{\&epsiv;})=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{12}u+{\mathrm{\&rho;}}_{22}U)-b\frac{\&part;(u-U)}{\&part;t},$

其中，b为耗散系数。

7.根据权利要求6所述的基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，所述第一纵波方程式满足如下公式：

$\left\{\begin{array}{c}{\&dtri;}^2\left(P\mathrm{\&theta;}+Q\mathrm{\&epsiv;}\right)=\frac{{\&part;}^2}{\&part;t^2}\left({\mathrm{\&rho;}}_{11}\mathrm{\&theta;}+{\mathrm{\&rho;}}_{12}\mathrm{\&epsiv;}\right)+b\frac{\&part;}{\&part;t}\left(\mathrm{\&theta;}-\mathrm{\&epsiv;}\right)\\ {\&dtri;}^2\left(Q\mathrm{\&theta;}+R\mathrm{\&epsiv;}\right)=\frac{{\&part;}^2}{\&part;t^2}\left({\mathrm{\&rho;}}_{12}\mathrm{\&theta;}+{\mathrm{\&rho;}}_{22}\mathrm{\&epsiv;}\right)-b\frac{\&part;}{\&part;t}\left(\mathrm{\&theta;}-\mathrm{\&epsiv;}\right)\end{array}\right.,$

所述第二纵波方程式满足如下公式：

$\left\{\begin{array}{c}{\&dtri;}^2(P\mathrm{\&theta;}+Q\mathrm{\&epsiv;})=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{11}\mathrm{\&theta;}+{\mathrm{\&rho;}}_{12}\mathrm{\&epsiv;})\\ {\&dtri;}^2(Q\mathrm{\&theta;}+R\mathrm{\&epsiv;})=\frac{{\&part;}^2}{\&part;t^2}({\mathrm{\&rho;}}_{12}\mathrm{\&theta;}+{\mathrm{\&rho;}}_{22}\mathrm{\&epsiv;})\end{array}\right.$

其中， $\&dtri;\&times;u=\mathrm{\&theta;},\&dtri;\&times;U=\mathrm{\&epsiv;}.$

8.根据权利要求7所述的基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，所述第一差分方程式满足如下公式：

$\left\{\begin{array}{c}{\mathrm{\&theta;}}^{n+1}(i,j)=\frac{(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)A^{\&prime;}-(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)B^{\&prime;}}{(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t)(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)-{(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)}^2}\\ {\mathrm{\&epsiv;}}^{n+1}(i,j)=\frac{(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t)B^{\&prime;}-(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)A^{\&prime;}}{(2{\mathrm{\&rho;}}_{11}+b\mathrm{\&Delta;}t)(2{\mathrm{\&rho;}}_{22}+b\mathrm{\&Delta;}t)-{(2{\mathrm{\&rho;}}_{12}-b\mathrm{\&Delta;}t)}^2}\end{array}\right.,$

其中，

$\begin{array}{c}{A}^{\&prime;}=2\{P\left(\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&theta;}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&theta;}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\&Delta;t}}^2\\ +Q\left(\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&epsiv;}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&epsiv;}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{11}\&lsqb;{\mathrm{\&theta;}}^{n-1}\left(i,j\right)-2{\mathrm{\&theta;}}^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{12}\&lsqb;{\mathrm{\&epsiv;}}^{n-1}\left(i,j\right)-2{\mathrm{\&epsiv;}}^n\left(i,j\right)\&rsqb;\}\\ +b\&lsqb;{\mathrm{\&theta;}}^{n-1}\left(i,j\right)-{\mathrm{\&epsiv;}}^{n-1}\left(i,j\right)\&rsqb;\mathrm{\&Delta;}t\\ =2A+b\&lsqb;{\mathrm{\&theta;}}^{n-1}\left(i,j\right)-{\mathrm{\&epsiv;}}^{n-1}\left(i,j\right)\&rsqb;\mathrm{\&Delta;}t\end{array},$

$\begin{array}{c}{B}^{\&prime;}=2\{Q\left(\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&theta;}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&theta;}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\&Delta;t}}^2\\ +R\left(\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&epsiv;}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&epsiv;}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{12}\&lsqb;{\mathrm{\&theta;}}^{n-1}\left(i,j\right)-2{\mathrm{\&theta;}}^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{22}\&lsqb;{\mathrm{\&epsiv;}}^{n-1}\left(i,j\right)-2{\mathrm{\&epsiv;}}^n\left(i,j\right)\&rsqb;\}\\ -b\&lsqb;{\mathrm{\&theta;}}^{n-1}\left(i,j\right)-{\mathrm{\&epsiv;}}^{n-1}\left(i,j\right)\&rsqb;\mathrm{\&Delta;}t\\ =2B-b\&lsqb;{\mathrm{\&theta;}}^{n-1}\left(i,j\right)-{\mathrm{\&epsiv;}}^{n-1}\left(i,j\right)\&rsqb;\mathrm{\&Delta;}t\end{array},$

其中，i是x方向的空间序号，j是z方向的空间序号，角标n表示时刻；h为空间差分步长，θⁿ(i,j)表示固相体应变θ在n时刻在(ih,jh)处的值；εⁿ(i,j)表示流相体应变ε在n时刻在(ih,jh)处的值，为差分系数，2N为差分精度，b为耗散系数。

9.根据权利要求8所述的基于有限差分法的三维TTI双相介质地震波场数值模拟方法，其特征在于，所述第二差分方程式满足如下公式：

$\left\{\begin{array}{c}{\mathrm{\&theta;}}^{n+1}\left(i,j\right)=\frac{{\mathrm{\&rho;}}_{22}A-{\mathrm{\&rho;}}_{12}B}{{\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2}\\ {\mathrm{\&epsiv;}}^{n+1}\left(i,j\right)=\frac{{\mathrm{\&rho;}}_{11}B-{\mathrm{\&rho;}}_{12}A}{{\mathrm{\&rho;}}_{11}{\mathrm{\&rho;}}_{22}-{\mathrm{\&rho;}}_{12}^2}\end{array}\right.,$

其中，

$\begin{array}{c}A=P\left(\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&theta;}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&theta;}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\&Delta;t}}^2\\ +Q\left(\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&epsiv;}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&epsiv;}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{11}\&lsqb;{\mathrm{\&theta;}}^{n-1}\left(i,j\right)-2{\mathrm{\&theta;}}^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{12}\&lsqb;{\mathrm{\&epsiv;}}^{n-1}\left(i,j\right)-2{\mathrm{\&epsiv;}}^n\left(i,j\right)\&rsqb;\end{array},$

$\begin{array}{c}B=Q\left(\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&theta;}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&theta;}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\&Delta;t}}^2\\ +R\left(\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&epsiv;}}^n\left(i+m,j\right)}{h^2}+\underset{m=-N}{\overset{N}{\mathrm{\&Sigma;}}}\frac{c_m^{\left(N\right)}{\mathrm{\&epsiv;}}^n\left(i,j+m\right)}{h^2}\right){\mathrm{\&Delta;t}}^2\\ -{\mathrm{\&rho;}}_{12}\&lsqb;{\mathrm{\&theta;}}^{n-1}\left(i,j\right)-2{\mathrm{\&theta;}}^n\left(i,j\right)\&rsqb;-{\mathrm{\&rho;}}_{22}\&lsqb;{\mathrm{\&epsiv;}}^{n-1}\left(i,j\right)-2{\mathrm{\&epsiv;}}^n\left(i,j\right)\&rsqb;\end{array}.$