CN106353797A - 一种高精度地震正演模拟方法 - Google Patents

Abstract

本发明提出了一种高精度的地震正演模拟的方法。包括：（1）建立地质模型；（2）对地质模型离散划分为波动方程交错网格差分格式；（3）建立弹性波方程交错网格差分格式；（4）构造新的粘弹模型，建立非均匀介质粘弹性波方程数值模拟；（5）边界条件处理；（6）有限差分解法稳定性分析；（7）进行观测系统设计，开展正演模拟。本发明的方法与常规网格高阶差分相比，可以进一步提高数值模拟的精度并压制数值频散。交错网格高阶有限差分法具有很好的模拟效果，计算效率较高，用交错网格高阶有限差分实现了弹性波以及粘弹性波的数值模拟，从得出的波场快照和炮记录中分析波在各种复杂介质内部的反射、透射、绕射、散射以及能量的衰减等运动学和动力学的各种细节特征。

Description

一种高精度地震正演模拟方法

技术领域

本发明属于油气勘探地震正演模拟技术领域，是一种比较高效的地震正演模拟方法。

现有技术

地震数值模拟是地震勘探和地震学的重要基础。地震数值模拟就是在假定地下介质结构模型和相应物理参数已知的情况下，模拟研究地震波在地下各种介质中的传播规律，并计算在地面或地下各观测点所应观测到的数值地震记录的一种地震模拟方法。这种地震数值模拟方法，在地震采集、处理与解释过程中，都是一个不可缺少的内容。

有限差分法是地震数值模拟比较常用的一种方法。有限差分法对计算区域网格化，通过数值求解描述地震波传播的微分方程来模拟波的传播。该方法在地震波模拟中使用最为广泛，主要问题是计算量比较大，对计算机内存要求较高。Alterman和Karal(1968)首先将有限差分法应用于层状介质弹性波传播的数值模拟中。此后，Boore(1972)又将有限差分法用于非均匀介质地震波传播的模拟。Dalain(1986)和Mufti et al(1990)讨论了用高阶差分方法解决声波正演模拟问题。此后，Bayliss(1986)、Levander(1988)采用了空间四阶有限差分法，来模拟弹性波传播的地震记录。为了进一步模拟地震波在实际非完全弹性地层内的传播，Carcione等(1988)提出了粘滞声波在地层中传播的模拟方法。Robertsson等(1994)给出了粘弹性波有限差分模拟方法。Carcione和Helle(1999)提出了孔隙粘弹性介质中地震波传播的交错网格有限差分模拟方法。Pitarka(1999)给出了三维各向同性介质中弹性波的矩形非规则交错网格有限差分模拟方法。

发明内容。

本发明是在总结现有方法具有各自的优劣基础上，提出了一种具有很好的模拟效果，计算效率较高、高精度的地震正演模拟的方法。

本发明的技术方案：

高精度的地震正演模拟的方法，包括：

(1)建立地质模型；

(2)对地质模型离散划分为波动方程交错网格差分格式；

(3)建立弹性波方程交错网格差分格式；

(4)构造新的粘弹模型，建立非均匀介质粘弹性波方程数值模拟；

(5)边界条件处理；

(6)有限差分解法稳定性分析；

(7)进行观测系统设计，开展正演模拟。

上述方案进一步：

(1)建立地质模型；

(2)将波场和物性参数定义在不同的网格点上，差分时其它点以所求点为中心做向前或向后差分，即在正演模拟时将波场分量与弹性参数分别定义在不同的网格点上，并将向前高阶差分和向后高阶差分应用于交错网格中即可得到高阶交错网格差分格式；

(3)采用交错网格，利用二维弹性波一阶速度-应力方程，计算速度与应力的离散值，通过各项精度的离散形式，建立弹性波方程交错网格差分格式；

(4)采用线性粘弹性模型中的开尔文粘弹性介质对粘弹性波动方程数值进行模拟，非均匀粘弹性介质中，将位移方程表示成一阶速度－应力方程组，采用交错网格，得2N阶空间差分精度、二阶时间差分精度交错网格高阶有限差分格式；

(5)边界条件处理采用了完全匹配层作为吸收边界，在区域的周围加上完全匹配层；

(6)有限差分解法稳定性分析；

(7)进行观测系统设计，开展正演模拟。

上述方案更进一步：

(1)建立地质模型；

(2)对地质模型离散划分为波动方程交错网格差分格式

波动方程交错网格差分格式按照A/B/C/D交叉分布并间隔排列的方式，将波场和物性参数定义在A、B、C、D所对应的网格点上，差分时其它点以所求点为中心做向前或向后差分，即在正演模拟时将波场分量与弹性参数分别定义在不同的网格点上，

设 $A=\frac{\∂C}{\∂x}+\frac{\∂D}{\∂z},$ 则可以得

$A(m,n)=\frac{1}{\mathrm{\Δ}x}\[C(m+\frac{1}{2},n)-C(m-\frac{1}{2},n)\]+\frac{1}{\mathrm{\Δ}z}\[D(m,n+\frac{1}{2})-D(m,n-\frac{1}{2})\]---\left(1\right)$

其它点的差分格式同理，向前高阶差分和向后高阶差分可表示为下式

$D_x^{+}f\left(n\right)=\frac{1}{\mathrm{\Δ}x}\[\underset{i=1}{\overset{i=N}{\Σ}}{a}_i\[f(n+i)-f(n-i+1)\]\]---\left(2\right)$

$D_x^{-}f\left(n\right)=\frac{1}{\mathrm{\Δ}x}\[\underset{i=1}{\overset{i=N}{\Σ}}{a}_i\[f(n+i-1)-f(n-i)\]\]---\left(3\right)$

其中，表示向前差分，表示向后差分，a_i为差分权系数，将向前高阶差分和向后高阶差分应用于交错网格中即可得到高阶交错网格差分格式；

(3)建立弹性波方程交错网格差分格式

二维弹性波一阶速度-应力方程为：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂V_x}{\∂t}=\frac{{\∂\mathrm{\σ}}_{\mathrm{xx}}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{\mathrm{xz}}}{\∂z}\\ \mathrm{\ρ}\frac{V_z}{\∂t}=\frac{{\∂\mathrm{\σ}}_{\mathrm{zx}}}{\∂x}+\frac{{\∂\mathrm{\σ}}_{\mathrm{zz}}}{\∂z}\\ \frac{{\∂\mathrm{\σ}}_{\mathrm{xx}}}{\∂t}=C_{11}\frac{\∂V_x}{\∂x}+C_{13}\frac{\∂V_z}{\∂z}\\ \frac{{\∂\mathrm{\σ}}_{\mathrm{zz}}}{\∂t}=C_{13}\frac{\∂V_x}{\∂x}+C_{33}\frac{\∂V_z}{\∂z}\\ \frac{{\∂\mathrm{\σ}}_{\mathrm{xz}}}{\∂t}=C_{44}\frac{\∂V_x}{\∂z}+C_{44}\frac{\∂V_y}{\∂x}\end{array}\right.---\left(4\right)$

采用交错网格，设分别是速度V_x，V_z与应力σ_xx、σ_zz、σ_xz的离散值，则方程组(4)中各项精度为O(Δt²+Δx⁸)的离散形式的差分格式为

$\begin{array}{c}{U}_{i,j}^{k+1/2}=U_{i,j}^{k-1/2}+\frac{\mathrm{\Δ}t}{{\mathrm{\ρ}}_{i,j}}\{\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[R_{i+(2n-1)/2,j}^k-R_{i-(2n-1)/2,j}^k\]\\ +\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[H_{i,j+(2n-1)/2}^k-H_{i,j-(2n-1)/2}^k\]\}\end{array}---\left(5\right)$

$\begin{array}{c}{V}_{i+1/2,j+1/2}^{k+1/2}=V_{i+1/2,j+1/2}^{k-1/2}+\frac{\mathrm{\Δ}t}{{\mathrm{\ρ}}_{i+1/2,j+1/2}}\{\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[H_{i+n,j+1/2}^k-H_{i-(n-1),j+1/2}^k\]\\ +\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[T_{i+1/2,j+n}^k-T_{i+1/2,j-(n-1)}^k\]\}\end{array}---\left(6\right)$

$\begin{array}{c}{R}_{i+1/2,j}^k=R_{i+1/2,j}^{k-1}+\mathrm{\Δ}t*\{C_{11}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[U_{i+n,j+1/2}^{k-1/2}-U_{i-(n-1),j}^{k-1/2}\]\\ +C_{13}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[V_{i+1/2,j+(2n-1)/2}^{k-1/2}-V_{i+1/2,j-(2n-1)/2}^{k-1/2}\]\}\end{array}---\left(7\right)$

$\begin{array}{c}{T}_{i+1/2,j}^k=T_{i+1/2,j}^{k-1}+\mathrm{\Δ}t*\{C_{33}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}(V_{i+1/2,j+(2n-1)/2}^{k-1/2}-V_{i+1/2,j-(2n-1)/2}^{k-1/2})\\ +C_{13}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}(U_{i+n,j}^{k-1/2}-U_{i-(n-1),j}^{k-1/2})\}\end{array}---\left(8\right)$

$\begin{array}{c}{H}_{i,j+1/2}^k=H_{i,j+1/2}^{k-1}+\mathrm{\Δ}t*\[C_{44}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}(U_{i,j+n}^{k-1/2}-U_{i,j-(n-1)}^{k-1/2})\\ +C_{44}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}(V_{i+(2n-1)/2,j+1/2}^{k+1/2}-V_{i-(2n-1)/2,j+1/2}^{k+1/2})\]\end{array}---\left(9\right)$

(4)构造新的粘弹模型，建立非均匀介质粘弹性波方程数值模拟

采用线性粘弹性模型中的开尔文粘弹性介质对粘弹性波动方程数值进行模拟，

对于开尔文体模型，在直角坐标系中，二维粘弹性波动方程为：

式中，为位移向量，ρ为密度函数，λ,μ为拉梅系数，λ′,μ′为粘滞系数，θ为体应变， $\▿=(\frac{\∂}{\∂x},\frac{\∂}{\∂y},\frac{\∂}{\∂z}),\mathrm{\Δ}=(\frac{{\∂}^2}{\∂x^2},\frac{{\∂}^2}{\∂y^2},\frac{{\∂}^2}{\∂z^2}),$

非均匀粘弹性介质中，将位移方程表示成一阶速度－应力方程组

$\left\{\begin{array}{c}\frac{\∂v_x}{\∂t}=\frac{1}{\mathrm{\ρ}(x,z)}(\frac{\∂{\mathrm{\τ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z})\\ \frac{\∂v_z}{\∂t}=\frac{1}{\mathrm{\ρ}(x,z)}(\frac{\∂{\mathrm{\τ}}_{xz}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{zz}}{\∂z})\\ \frac{\∂{\mathrm{\τ}}_{xx}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_x}{\∂x}+\mathrm{\λ}\frac{\∂v_z}{\∂z}+({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})\frac{{\∂}^2{v}_x}{\∂x\∂t}+{\mathrm{\λ}}^{\′}\frac{{\∂}^2{v}_z}{\∂z\∂t}\\ \frac{\∂{\mathrm{\τ}}_{zz}}{\∂t}=\mathrm{\λ}\frac{\∂v_x}{\∂x}+(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_z}{\∂z}+{\mathrm{\λ}}^{\′}\frac{{\∂}^2{v}_x}{\∂x\∂t}+({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})\frac{{\∂}^2{v}_z}{\∂z\∂t}\\ \frac{\∂{\mathrm{\τ}}_{xz}}{\∂t}=\mathrm{\μ}\frac{\∂v_x}{\∂z}+\mathrm{\μ}\frac{\∂v_z}{\∂x}{\mathrm{\μ}}^{\′}\frac{{\∂}^2{v}_x}{\∂z\∂t}+{\mathrm{\μ}}^{\′}\frac{{\∂}^2{v}_z}{\∂x\∂t}\end{array}\right.---\left(11\right)$

其中，λ+2μ＝ρv_p ²， ${\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′}=\frac{{{\mathrm{\ρv}}_p}^2}{Q_p\mathrm{\ω}},$ μ＝ρv_s ²， ${\mathrm{\μ}}^{\′}=\frac{{{\mathrm{\ρv}}_s}^2}{Q_s\mathrm{\ω}},$

λ＝ρv_p ²-2ρv_s ²， ${\mathrm{\λ}}^{\′}=\frac{{{\mathrm{\ρv}}_p}^2}{Q_p\mathrm{\ω}}-2\frac{{{\mathrm{\ρv}}_s}^2}{Q_s\mathrm{\ω}},$

其中，v_p表示纵波速度，Q_p表示纵波品质因子，v_s表示横波速度，Q_s表示横波品质因子，ω为圆频率；

采用交错网格，得2N阶空间差分精度、二阶时间差分精度交错网格高阶有限差分格式如式：

$v_x^t(i^{+},j)=v_x^{t-1}(i^{+},j)+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x\mathrm{\ρ}}\{D_x^{+}\[{\mathrm{\τ}}_{xx}^{t-}(i,j)\]\}-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z\mathrm{\ρ}}\{D_z^{-}\[{\mathrm{\τ}}_{zz}^{t-}(i^{+},j^{+})\]\}---\left(12\right)$

$v_z^t(i,j^{+})=v_z^{t-1}(i,j^{+})+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x\mathrm{\ρ}}\{D_x^{-}\[{\mathrm{\τ}}_{xz}^{t-}(i^{+},j^{+})\]\}-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z\mathrm{\ρ}}\{D_z^{+}\[{\mathrm{\τ}}_{zz}^{t-}(i,j)\]\}---\left(13\right)$

$\begin{array}{c}{\mathrm{\τ}}_{xx}^{t+}(i,j)={\mathrm{\τ}}_{xx}^{t-}(i,j)+\frac{\mathrm{\Δ}t(\mathrm{\λ}+2\mathrm{\μ})}{\mathrm{\Δ}x}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t\mathrm{\λ}}{\mathrm{\Δ}x}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\\ +\frac{\mathrm{\Δ}t({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{{\mathrm{\Δt\λ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\end{array}---\left(14\right)$

$\begin{array}{c}{\mathrm{\τ}}_{zz}^{t+}(i,j)={\mathrm{\τ}}_{zz}^{t-}(i,j)+\frac{\mathrm{\Δ}t\mathrm{\λ}}{\mathrm{\Δ}x}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t(\mathrm{\λ}+2\mathrm{\μ})}{\mathrm{\Δ}x}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\\ +\frac{{\mathrm{\Δt\λ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\end{array}---\left(15\right)$

$\begin{array}{c}{\mathrm{\τ}}_{xz}^{t+}(i,j)={\mathrm{\τ}}_{xx}^{t-}(i,j)+\frac{\mathrm{\Δ}t\mathrm{\μ}}{\mathrm{\Δ}x}\{D_z^{+}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t\mathrm{\μ}}{\mathrm{\Δ}x}\{D_x^{+}\[v_z^t(i,j^{+})\]\}\\ +\frac{{\mathrm{\Δt\μ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{+}\[v_x^t(i^{+},j)\]\}+\frac{{\mathrm{\Δt\μ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{+}\[v_z^t(i,j^{+})\]\}\end{array}---\left(16\right)$

其中上标“+”表示向前差分，“-“表示向后差分；

(5)对于边界条件的处理，采用了完全匹配层边界条件

利用完全匹配层作为吸收边界是在区域ABCD四周引入完全匹配层，在区域四个边角区域1中，令d_x≠0；d_z≠0，速度V都等于角点的速度，在z方向区域2中，令d_x＝0；d_z≠0,速度V在z方向为常数，在x方向和边界的速度相等，在x方向区域3中，令d_x≠0；d_z＝0，速度V在x方向为常数，在z方向和边界的速度相等；

对于二维弹性波方程而言，这里V_x、V_z、σ_xx、σ_zz和σ_xz都进行分解，分解结果如下：

$V_x=V_x^x+V_x^z,V_z=V_z^x+V_z^z,$

${\mathrm{\σ}}_{xx}={\mathrm{\σ}}_{xx}^x+{\mathrm{\σ}}_{xx}^z,{\mathrm{\σ}}_{zz}={\mathrm{\σ}}_{zz}^x+{\mathrm{\σ}}_{zz}^z,{\mathrm{\σ}}_{xz}={\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z.$

二维弹性波PML公式：

$\left\{\begin{array}{c}\frac{\∂V_x^x}{\∂t}+d_x{V}_x^x=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xx}^x+{\mathrm{\σ}}_{xx}^z)}{\∂x}\\ \frac{\∂V_x^z}{\∂t}+d_z{V}_x^z=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z)}{\∂z}\\ \frac{\∂V_z^x}{\∂t}+d_x{V}_z^x=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z)}{\∂x}\\ \frac{\∂V_z^z}{\∂t}+d_z{V}_z^z=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{zz}^x+{\mathrm{\σ}}_{zz}^z)}{\∂z}\end{array}\right.---\left(17\right)$

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\σ}}_{xx}^x}{\∂t}+d_x{\mathrm{\σ}}_{xx}^x=C_{11}\frac{\∂(V_x^x+V_x^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xx}^z}{\∂t}+d_z{\mathrm{\σ}}_{xx}^z=C_{13}\frac{\∂(V_z^x+V_z^z)}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{zz}^x}{\∂t}+d_x{\mathrm{\σ}}_{zz}^x=C_{13}\frac{\∂(V_x^x+V_x^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{zz}^z}{\∂t}+d_z{\mathrm{\σ}}_{zz}^z=C_{33}\frac{\∂(V_z^x+V_z^z)}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{xz}^x}{\∂t}+d_x{\mathrm{\σ}}_{xz}^x=C_{44}\frac{\∂(V_z^x+L_z^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xz}^z}{\∂t}+d_z{\mathrm{\σ}}_{xz}^z=C_{44}\frac{\∂(V_x^x+V_x^z)}{\∂z}\end{array}\right.---\left(18\right)$

(6)有限差分解法稳定性分析

2L阶精度的稳定性条件如式(19)，式中，Δt为时间间隔，Δx和Δz为网格大小，v_p为纵波速度，a_l为L个差分系数，

${\mathrm{\Δtv}}_p\sqrt{\frac{1}{{\mathrm{\Δx}}^2}+\frac{1}{{\mathrm{\Δz}}^2}}\≤\frac{1}{\underset{l=1}{\overset{L}{\Σ}}\left|a_l\right|}---\left(19\right)$

(7)进行观测系统设计，开展正演模拟。

本发明的方法在波动方程数值模拟中，有限差分法是最常用的一种方法，而其中的交错网格高阶差分方法与常规网格高阶差分相比，可以进一步提高数值模拟的精度并压制数值频散。交错网格高阶有限差分法具有很好的模拟效果，计算效率较高，用交错网格高阶有限差分实现了弹性波以及粘弹性波的数值模拟，从得出的波场快照和炮记录中分析波在各种复杂介质内部的反射、透射、绕射、散射以及能量的衰减等运动学和动力学的各种细节特征。

附图说明

图1交错网格差分示意图；

图2完全匹配层吸收边界示意图；

图3本发明的一种实施例的流程图；

图4XX地质模型图；

图5弹性波正演模拟水平分量图；

图6弹性波正演模拟垂直分量图。

具体实施方式

参照附图3，本发明的具体实施实例为：

(1)建立地质模型

(2)对地质模型离散划分为波动方程交错网格差分格式

波动方程交错网格差分格式如图1所示，图1中的每个点都可以通过周围点的向前差分或向后差分来得到。当将波场和物性参数定义在不同的网格点时，那么，就可以得到交错网格的差分格式。如图中A、B、C、D所对应的点都定义在不同的网格点上，差分时其它点以所求点为中心做向前或向后差分，即在正演模拟时将波场分量与弹性参数分别定义在不同的网格点上，例如

$A=\frac{\∂C}{\∂x}+\frac{\∂D}{\∂z},$ 则可以得

$A(m,n)=\frac{1}{\mathrm{\Δ}x}\[C(m+\frac{1}{2},n)-C(m-\frac{1}{2},n)\]+\frac{1}{\mathrm{\Δ}z}\[D(m,n+\frac{1}{2})-D(m,n-\frac{1}{2})\]---\left(1\right)$

其它点的差分格式同理，向前高阶差分和向后高阶差分可表示为下式

$D_x^{+}f\left(n\right)=\frac{1}{\mathrm{\Δx}}\[\underset{i=1}{\overset{i=N}{\Σ}}{a}_i\[f(n+i)-f(n-i+1)\]\]---\left(2\right)$

$D_x^{-}f\left(n\right)=\frac{1}{\mathrm{\Δ}x}\[\underset{i=1}{\overset{i=N}{\Σ}}{a}_i\[f(n+i-1)-f(n-i)\]\]---\left(3\right)$

其中，表示向前差分，表示向后差分，a_i为差分权系数。将向前高阶差分和向后高阶差分应用于交错网格中即可得到高阶交错网格差分格式。

(3)建立弹性波方程交错网格差分格式

二维弹性波一阶速度-应力方程为：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂V_x}{\∂t}=\frac{\∂{\mathrm{\σ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{xz}}{\∂z}\\ \mathrm{\ρ}\frac{\∂V_z}{\∂t}=\frac{\∂{\mathrm{\σ}}_{zx}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{xx}}{\∂t}=C_{11}\frac{\∂V_x}{\∂x}+C_{13}\frac{\∂V_z}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{zz}}{\∂t}=C_{13}\frac{\∂V_x}{\∂x}+C_{33}\frac{\∂V_z}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{xz}}{\∂t}=C_{44}\frac{\∂V_x}{\∂z}+C_{44}\frac{\∂V_y}{\∂x}\end{array}\right.---\left(4\right)$

采用交错网格，设分别是速度V_x，V_z与应力σ_xx、σ_zz、σ_xz的离散值，则方程组(4)中各项精度为O(Δt²+Δx⁸)的离散形式的差分格式为

$\begin{array}{c}{U}_{i,j}^{k+1/2}=U_{i,j}^{k-1/2}+\frac{\mathrm{\Δ}t}{{\mathrm{\ρ}}_{i,j}}\{\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[R_{i+(2n-1)/2,j}^k-R_{i-(2n-1)/2,j}^k\]\\ +\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[H_{i,j+(2n-1)/2}^k-H_{i,j-(2n-1)/2}^k\]\}\end{array}---\left(5\right)$

$\begin{array}{c}{V}_{i+1/2,j+1/2}^{k+1/2}=V_{i+1/2,j+1/2}^{k-1/2}+\frac{\mathrm{\Δ}t}{{\mathrm{\ρ}}_{i+1/2,j+1/2}}\{\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[H_{i+n,j+1/2}^k-H_{i-(n-1),j+1/2}^k\]\\ +\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[T_{i+1/2,j+n}^k-T_{i+1/2,j-(n-1)}^k\]\}\end{array}---\left(6\right)$

$\begin{array}{c}{R}_{i+1/2,j}^k=R_{i+1/2,j}^{k-1}+\mathrm{\Δ}t*\{C_{11}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[U_{i+n,j+1/2}^{k-1/2}-U_{i-(n-1),j}^{k-1/2}\]\\ +C_{13}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[V_{i+1/2,j+(2n-1)/2}^{k-1/2}-V_{i+1/2,j-(2n-1)/2}^{k-1/2}\]\}\end{array}---\left(7\right)$

$\begin{array}{c}{T}_{i+1/2,j}^k=T_{i+1/2,j}^{k-1}+\mathrm{\Δ}t*\{C_{33}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}(V_{i+1/2,j+(2n-1)/2}^{k-1/2}-V_{i+1/2,j-(2n-1)/2}^{k-1/2})\\ +C_{13}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}(U_{i+n,j}^{k-1/2}-U_{i-(n-1),j}^{k-1/2})\}\end{array}---\left(8\right)$

$\begin{array}{c}{H}_{i,j+1/2}^k=H_{i,j+1/2}^{k-1}+\mathrm{\Δ}t*\[C_{44}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}(U_{i,j+n}^{k-1/2}-U_{i,j-(n-1)}^{k-1/2})\\ +C_{44}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}(V_{i+(2n-1)/2,j+1/2}^{k+1/2}-V_{i-(2n-1)/2,j+1/2}^{k+1/2})\]\end{array}---\left(9\right)$

(4)构造新的粘弹模型，建立非均匀介质粘弹性波方程数值模拟

由于实际介质并非是完全弹性的，所以粘弹性比弹性更接近实际介质的应力-应变关系。采用线性粘弹性模型中的开尔文(kelvin)粘弹性介质对粘弹性波动方程数值模拟进行了研究。

对于开尔文体模型，在直角坐标系中，二维粘弹性波动方程为：

式中，为位移向量，ρ为密度函数，λ,μ为拉梅系数，λ′,μ′为粘滞系数，θ为体应变， $\▿=(\frac{\∂}{\∂x},\frac{\∂}{\∂y},\frac{\∂}{\∂z}),\mathrm{\Δ}=(\frac{{\∂}^2}{\∂x^2},\frac{{\∂}^2}{\∂y^2},\frac{{\∂}^2}{\∂z^2}).$

非均匀粘弹性介质中。可将位移方程表示成一阶速度－应力方程组

$\left\{\begin{array}{c}\frac{\∂v_x}{\∂t}=\frac{1}{\mathrm{\ρ}(x,z)}(\frac{\∂{\mathrm{\τ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z})\\ \frac{\∂v_z}{\∂t}=\frac{1}{\mathrm{\ρ}(x,z)}(\frac{\∂{\mathrm{\τ}}_{xz}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{zz}}{\∂z})\\ \frac{\∂{\mathrm{\τ}}_{xx}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_x}{\∂x}+\mathrm{\λ}\frac{\∂v_z}{\∂z}+({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})\frac{{\∂}^2{v}_x}{\∂x\∂t}+{\mathrm{\λ}}^{\′}\frac{{\∂}^2{v}_z}{\∂z\∂t}\\ \frac{\∂{\mathrm{\τ}}_{zz}}{\∂t}=\mathrm{\λ}\frac{\∂v_x}{\∂x}+(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_z}{\∂z}+{\mathrm{\λ}}^{\′}\frac{{\∂}^2{v}_x}{\∂x\∂t}+({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})\frac{{\∂}^2{v}_z}{\∂z\∂t}\\ \frac{\∂{\mathrm{\τ}}_{xz}}{\∂t}=\mathrm{\μ}\frac{\∂v_x}{\∂z}+\mathrm{\μ}\frac{\∂v_z}{\∂x}{\mathrm{\μ}}^{\′}\frac{{\∂}^2{v}_x}{\∂z\∂t}+{\mathrm{\μ}}^{\′}\frac{{\∂}^2{v}_z}{\∂x\∂t}\end{array}\right.---\left(11\right)$

其中，λ+2μ＝ρv_p ²， ${\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′}=\frac{{{\mathrm{\ρv}}_p}^2}{Q_p\mathrm{\ω}},$ μ＝ρv_s ²， ${\mathrm{\μ}}^{\′}=\frac{{{\mathrm{\ρv}}_s}^2}{Q_s\mathrm{\ω}},$

λ＝ρv_p ²-2ρv_s ²， ${\mathrm{\λ}}^{\′}=\frac{{{\mathrm{\ρv}}_p}^2}{Q_p\mathrm{\ω}}-2\frac{{{\mathrm{\ρv}}_s}^2}{Q_s\mathrm{\ω}}.$

其中，v_p表示纵波速度，Q_p表示纵波品质因子，v_s表示横波速度，Q_s表示横波品质因子，ω为圆频率。

采用交错网格，得2N阶空间差分精度、二阶时间差分精度交错网格高阶有限差分格式如式：

$v_x^t(i^{+},j)=v_x^{t-1}(i^{+},j)+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x\mathrm{\ρ}}\{D_x^{+}\[{\mathrm{\τ}}_{xx}^{t-}(i,j)\]\}-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z\mathrm{\ρ}}\{D_z^{-}-\[{\mathrm{\τ}}_{zz}^{t-}(i^{+},j^{+})\]\}---\left(12\right)$

$v_z^t(i,j^{+})=v_z^{t-1}(i,j^{+})+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x\mathrm{\ρ}}\{D_x^{-}\[{\mathrm{\τ}}_{xz}^{t-}(i^{+},j^{+})\]\}-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z\mathrm{\ρ}}\{D_z^{+}\[{\mathrm{\τ}}_{zz}^{t-}(i,j)\]\}---\left(13\right)$

$\begin{array}{c}{\mathrm{\τ}}_{xx}^{t+}(i,j)={\mathrm{\τ}}_{xx}^{t-}(i,j)+\frac{\mathrm{\Δ}t(\mathrm{\λ}+2\mathrm{\μ})}{\mathrm{\Δ}x}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t\mathrm{\λ}}{\mathrm{\Δ}x}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\\ +\frac{\mathrm{\Δ}t({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{{\mathrm{\Δt\λ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\end{array}---\left(14\right)$

$\begin{array}{c}{\mathrm{\τ}}_{zz}^{t+}(i,j)={\mathrm{\τ}}_{zz}^{t-}(i,j)+\frac{\mathrm{\Δ}t\mathrm{\λ}}{\mathrm{\Δ}x}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t(\mathrm{\λ}+2\mathrm{\μ})}{\mathrm{\Δ}x}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\\ +\frac{{\mathrm{\Δt\λ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\end{array}---\left(15\right)$

$\begin{array}{c}{\mathrm{\τ}}_{xz}^{t+}(i,j)={\mathrm{\τ}}_{xx}^{t-}(i,j)+\frac{\mathrm{\Δ}t\mathrm{\μ}}{\mathrm{\Δ}x}\{D_z^{+}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t\mathrm{\μ}}{\mathrm{\Δ}x}\{D_z^{+}\[v_z^t(i,j^{+})\]\}\\ +\frac{{\mathrm{\Δt\μ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{+}\[v_x^t(i^{+},j)\]\}+\frac{{\mathrm{\Δt\μ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{+}\[v_z^t(i,j^{+})\]\}\end{array}---\left(16\right)$

其中上标“+”表示向前差分，“-“表示向后差分。

(5)边界条件处理

对于边界条件的处理，采用了完全匹配层边界条件。

利用完全匹配层作为吸收边界的基本做法是在所研究区域的四周引入完全匹配层。如图1所示，区域ABCD为所要研究的区域，即我们要在此区域中研究波的传播问题。在区域的周围加上完全匹配层，在区域1中，令d_x≠0；d_z≠0，速度V都等于角点的速度。在区域2中，令d_x＝0；d_z≠0,速度V在z方向为常数，在x方向和边界的速度相等。在区域3中，令d_x≠0；d_z＝0，速度V在x方向为常数，在z方向和边界的速度相等。这样在计算边界的周围都有完全匹配层吸收介质，波由区域内通过边界传播到完全匹配层时，不会产生任何反射。波在完全匹配层中传播时，不会产生反射，并且按传播距离的指数规律衰减。当波传播到完全匹配层的边界时，波场近似为零，也不会产生反射。

对于二维弹性波方程而言，这里V_x、V_z、σ_xx、σ_zz和σ_xz都进行分解，这是由于弹性波方程中，每个等式右边同时包含了沿x方向的偏微分和沿z方向的偏微分。分解结果如下：

$V_x=V_x^x+V_x^z,V_z=V_z^x+V_z^z,$

${\mathrm{\σ}}_{xx}={\mathrm{\σ}}_{xx}^x+{\mathrm{\σ}}_{xx}^z,{\mathrm{\σ}}_{zz}={\mathrm{\σ}}_{zz}^x+{\mathrm{\σ}}_{zz}^z,{\mathrm{\σ}}_{xz}={\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z.$

二维弹性波PML公式：

$\left\{\begin{array}{c}\frac{\∂V_x^x}{\∂t}+d_x{V}_x^x=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xx}^x+{\mathrm{\σ}}_{xx}^z)}{\∂x}\\ \frac{\∂V_x^z}{\∂t}+d_z{V}_x^z=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z)}{\∂z}\\ \frac{\∂V_z^x}{\∂t}+d_x{V}_z^x=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z)}{\∂x}\\ \frac{\∂V_z^z}{\∂t}+d_z{V}_z^z=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{zz}^x+{\mathrm{\σ}}_{zz}^z)}{\∂z}\end{array}\right.---\left(17\right)$

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\σ}}_{xx}^x}{\∂t}+d_x{\mathrm{\σ}}_{xx}^x=C_{11}\frac{\∂(V_x^x+V_x^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xx}^z}{\∂t}+d_x{\mathrm{\σ}}_{xx}^z=C_{13}\frac{\∂(V_z^x+V_z^z)}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{zz}^x}{\∂t}+d_x{\mathrm{\σ}}_{zz}^x=C_{13}\frac{\∂(V_x^x+V_x^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{zz}^z}{\∂t}+d_z{\mathrm{\σ}}_{zz}^z=C_{33}\frac{\∂(V_z^x+V_z^z)}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{xz}^x}{\∂t}+d_x{\mathrm{\σ}}_{xz}^x=C_{44}\frac{\∂(V_z^x+V_z^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xz}^z}{\∂t}+d_z{\mathrm{\σ}}_{xz}^z=C_{44}\frac{\∂(V_x^x+V_x^z)}{\∂z}\end{array}\right.---\left(18\right)$

(6)有限差分解法稳定性分析

稳定性问题是数值求解波动方程的基本问题，由于差分计算过程中数值参数选择不合理，可能产生无物理意义的按指数增大的数值计算结果，造成模拟结果网格频散严重，影响对问题的分析，严重时会造成溢出而使计算中断。因此，对一种数值解法，要知道使计算稳定的离散参数的取值范围，也就是要分析解法的稳定性。

2L阶精度的稳定性条件如式(19)，式中，Δt为时间间隔，Δx和Δz为网格大小，v_p为纵波速度。a_l为L个差分系数。由公式可见，在各向同性弹性介质中，纵波速度和横波速度不同时包含在稳定性条件中，稳定性条件与泊松比无关。

${\mathrm{\Δtv}}_p\sqrt{\frac{1}{{\mathrm{\Δx}}^2}+\frac{1}{{\mathrm{\Δz}}^2}}\≤\frac{1}{\underset{l=1}{\overset{L}{\Σ}}\left|a_l\right|}---\left(19\right)$

(7)进行观测系统设计，开展正演模拟。

本次研究对XX地质模型进行测试如图4。图5和图6分别为弹性波方程交错网格高阶有限差分正演模拟的地面地震正演模拟水平分量和垂直分量。由图可以看出弹性波方程正演结果不但存在反射纵波，还有反射横波、直达横波以及多种转换波存在，其波场更加丰富。

Claims (3)

1.高精度的地震正演模拟的方法，其特征是包括：

(1)建立地质模型；

(2)对地质模型离散划分为波动方程交错网格差分格式；

(3)建立弹性波方程交错网格差分格式；

(4)构造新的粘弹模型，建立非均匀介质粘弹性波方程数值模拟；

(5)边界条件处理；

(6)有限差分解法稳定性分析；

(7)进行观测系统设计，开展正演模拟。

2.根据权利要求1所述的高精度的地震正演模拟的方法，其特征是包括：

(1)建立地质模型；

(2)将波场和物性参数定义在不同的网格点上，差分时其它点以所求点为中心做向前或向后差分，即在正演模拟时将波场分量与弹性参数分别定义在不同的网格点上，并将向前高阶差分和向后高阶差分应用于交错网格中即可得到高阶交错网格差分格式；

(3)采用交错网格，利用二维弹性波一阶速度-应力方程，计算速度与应力的离散值，通过各项精度的离散形式，建立弹性波方程交错网格差分格式；

(4)采用线性粘弹性模型中的开尔文粘弹性介质对粘弹性波动方程数值进行模拟，非均匀粘弹性介质中，将位移方程表示成一阶速度－应力方程组，采用交错网格，得2N阶空间差分精度、二阶时间差分精度交错网格高阶有限差分格式；

(5)边界条件处理采用了完全匹配层作为吸收边界，在区域的周围加上完全匹配层；

(6)有限差分解法稳定性分析；

(7)进行观测系统设计，开展正演模拟。

3.根据权利要求1或2所述的高精度的地震正演模拟的方法，其特征是：

(1)建立地质模型；

(2)对地质模型离散划分为波动方程交错网格差分格式

波动方程交错网格差分格式按照A/B/C/D交叉分布并间隔排列的方式，将波场和物性参数定义在A、B、C、D所对应的网格点上，差分时其它点以所求点为中心做向前或向后差分，即在正演模拟时将波场分量与弹性参数分别定义在不同的网格点上，

设 $A=\frac{\∂C}{\∂x}+\frac{\∂D}{\∂z},$ 则可以得

$A(m,n)=\frac{1}{\mathrm{\Δ}x}\[C(m+\frac{1}{2},n)-C(m-\frac{1}{2},n)\]+\frac{1}{\mathrm{\Δ}z}\[D(m,n+\frac{1}{2})-D(m,n-\frac{1}{2})\]---\left(1\right)$

其它点的差分格式同理，向前高阶差分和向后高阶差分可表示为下式

$D_x^{+}f\left(n\right)=\frac{1}{\mathrm{\Δ}x}\[\underset{i=1}{\overset{i=N}{\Σ}}{a}_i\[f(n+i)-f(n-i+1)\]\]---\left(2\right)$

$D_x^{-}f\left(n\right)=\frac{1}{\mathrm{\Δ}x}\[\underset{i=1}{\overset{i=N}{\Σ}}{a}_i\[f(n+i-1)-f(n-i)\]\]---\left(3\right)$

其中，表示向前差分，表示向后差分，a_i为差分权系数，将向前高阶差

分和向后高阶差分应用于交错网格中即可得到高阶交错网格差分格式；

(3)建立弹性波方程交错网格差分格式

二维弹性波一阶速度-应力方程为：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂V_x}{\∂t}=\frac{\∂{\mathrm{\σ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{xz}}{\∂z}\\ \mathrm{\ρ}\frac{\∂V_z}{\∂t}=\frac{\∂{\mathrm{\σ}}_{zx}}{\∂x}+\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{xx}}{\∂t}=C_{11}\frac{\∂V_x}{\∂x}+C_{13}\frac{\∂V_z}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{zz}}{\∂t}=C_{13}\frac{\∂V_x}{\∂x}+C_{33}\frac{\∂V_z}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{xz}}{\∂t}=C_{44}\frac{\∂V_x}{\∂z}+C_{44}\frac{\∂V_y}{\∂x}\end{array}\right.---\left(4\right)$

采用交错网格，设分别是速度V_x，V_z与应力σ_xx、σ_zz、σ_xz的离散值，则方程组(4)中各项精度为O(Δt²+Δx⁸)的离散形式的差分格式为

$\begin{array}{c}{U}_{i,j}^{k+1/2}=U_{i,j}^{k-1/2}+\frac{\mathrm{\Δ}t}{{\mathrm{\ρ}}_{i,j}}\{\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[R_{i+(2n-1)/2,j}^k-R_{i-(2n-1)/2,j}^k\]\\ +\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[H_{i,j+(2n-1)/2}^k-H_{i,j-(2n-1)/2}^k\]\}\end{array}---\left(5\right)$

$\begin{array}{c}{V}_{i+1/2,j+1/2}^{k+1/2}=V_{i+1/2,j+1/2}^{k-1/2}+\frac{\mathrm{\Δ}t}{{\mathrm{\ρ}}_{i+1/2,j+1/2}}\{\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[H_{i+n,j+1/2}^k-H_{i-(n-1),j+1/2}^k\]\\ +\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[T_{i+1/2,j+n}^k-T_{i+1/2,j-(n-1)}^k\]\end{array}---\left(6\right)$

$\begin{array}{c}{R}_{i+1/2,j}^k=R_{i+1/2,j}^{k-1}+\mathrm{\Δ}t*\{C_{11}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}\[U_{i+n,j+1/2}^{k-1/2}-U_{i-(n-1),j}^{k-1/2}\]\\ +C_{13}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}\[V_{i+1/2,j+(2n-1)/2}^{k-1/2}-V_{i+1/2,j-(2n-1)/2}^{k-1/2}\]\}\end{array}---\left(7\right)$

$\begin{array}{c}{T}_{i+1/2,j}^k=T_{i+1/2,j}^{k-1}+\mathrm{\Δ}t*\{C_{33}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}(V_{i+1/2,j+(2n-1)/2}^{k-1/2}-V_{i+1/2,j-(2n-1)/2}^{k-1/2})\\ +C_{13}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}(U_{i+n,j}^{k-1/2}-U_{i-(n-1),j}^{k-1/2})\end{array}---\left(8\right)$

$\begin{array}{c}{H}_{i,j+1/2}^k=H_{i,j+1/2}^{k-1}+\mathrm{\Δ}t*\[C_{44}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}z}{C}_n^{\left(N\right)}(U_{i,j+n}^{k-1/2}-U_{i,j-(n-1)}^{k-1/2})\\ +C_{44}\underset{n=1}{\overset{N}{\Σ}}\frac{1}{\mathrm{\Δ}x}{C}_n^{\left(N\right)}(V_{i+(2n-1)/2,j+1/2}^{k+1/2}-V_{i-(2n-1)/2,j+1/2}^{k+1/2})\]\end{array}---\left(9\right)$

(4)构造新的粘弹模型，建立非均匀介质粘弹性波方程数值模拟

采用线性粘弹性模型中的开尔文粘弹性介质对粘弹性波动方程数值进行模拟，

对于开尔文体模型，在直角坐标系中，二维粘弹性波动方程为：

式中，为位移向量，ρ为密度函数，λ,μ为拉梅系数，λ′,μ′为粘滞系数，θ为体应变， $\▿=(\frac{\∂}{\∂x},\frac{\∂}{\∂y},\frac{\∂}{\∂z}),\mathrm{\Δ}=(\frac{{\∂}^2}{\∂x^2},\frac{{\∂}^2}{\∂y^2},\frac{{\∂}^2}{\∂z^2}),$

非均匀粘弹性介质中，将位移方程表示成一阶速度－应力方程组

$\{\begin{array}{c}\frac{\∂v_x}{\∂t}=\frac{1}{\mathrm{\ρ}(x,z)}(\frac{\∂{\mathrm{\τ}}_{xx}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂z})\\ \frac{\∂v_z}{\∂t}=\frac{1}{\mathrm{\ρ}(x,z)}(\frac{\∂{\mathrm{\τ}}_{xz}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{zz}}{\∂z})\\ \frac{\∂{\mathrm{\τ}}_{xx}}{\∂t}=(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_x}{\∂x}+\mathrm{\λ}\frac{\∂v_z}{\∂z}+({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})\frac{{\∂}^2{v}_x}{\∂x\∂t}+{\mathrm{\λ}}^{\′}\frac{{\∂}^2{v}_z}{\∂z\∂t}\\ \frac{\∂{\mathrm{\τ}}_{zz}}{\∂t}=\mathrm{\λ}\frac{\∂v_x}{\∂x}+(\mathrm{\λ}+2\mathrm{\μ})\frac{\∂v_z}{\∂z}+{\mathrm{\λ}}^{\′}\frac{{\∂}^2{v}_x}{\∂x\∂t}+({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})\frac{{\∂}^2{v}_z}{\∂z\∂t}\\ \frac{\∂{\mathrm{\τ}}_{xz}}{\∂t}=\mathrm{\μ}\frac{\∂v_x}{\∂z}+\mathrm{\μ}\frac{\∂v_z}{\∂x}+{\mathrm{\μ}}^{\′}\frac{{\∂}^2{v}_x}{\∂z\∂t}+{\mathrm{\μ}}^{\′}\frac{{\∂}^2{v}_z}{\∂x\∂t}\end{array}---\left(11\right)$

其中， $\mathrm{\λ}+2\mathrm{\μ}={{\mathrm{\ρv}}_p}^2,{\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′}=\frac{{{\mathrm{\ρv}}_p}^2}{Q_p\mathrm{\ω}},\mathrm{\μ}={{\mathrm{\ρv}}_s}^2,{\mathrm{\μ}}^{\′}=\frac{{{\mathrm{\ρv}}_s}^2}{Q_s\mathrm{\ω}},$

$\mathrm{\λ}={{\mathrm{\ρv}}_p}^2-2{{\mathrm{\ρv}}_s}^2,{\mathrm{\λ}}^{\′}=\frac{{{\mathrm{\ρv}}_p}^2}{Q_p\mathrm{\ω}}-2\frac{{{\mathrm{\ρv}}_s}^2}{Q_s\mathrm{\ω}},$

其中，v_p表示纵波速度，Q_p表示纵波品质因子，v_s表示横波速度，Q_s表示横波品质因子，ω为圆频率；

采用交错网格，得2N阶空间差分精度、二阶时间差分精度交错网格高阶有

限差分格式如式：

$v_x^t(i^{+},j)=v_x^{t-1}(i^{+},j)+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x\mathrm{\ρ}}\{D_x^{+}\[{\mathrm{\τ}}_{xx}^{t-}(i,j)\]\}-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z\mathrm{\ρ}}\{D_z^{-}\[{\mathrm{\τ}}_{zz}^{t-}(i^{+},j^{+})\]\}---\left(12\right)$

$v_z^t(i,j^{+})=v_z^{t-1}(i,j^{+})+\frac{\mathrm{\Δt}}{\mathrm{\Δ}x\mathrm{\ρ}}\{D_x^{-}\[{\mathrm{\τ}}_{xz}^{t-}(i^{+},j^{+})\]\}-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z\mathrm{\ρ}}\{D_z^{+}\[{\mathrm{\τ}}_{zz}^{t-}(i,j)\]\}---\left(13\right)$

$\begin{array}{c}{\mathrm{\τ}}_{xx}^{t+}(i,j)={\mathrm{\τ}}_{xx}^{t-}(i,j)+\frac{\mathrm{\Δ}t(\mathrm{\λ}+2\mathrm{\μ})}{\mathrm{\Δ}x}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t\mathrm{\λ}}{\mathrm{\Δ}x}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\\ +\frac{\mathrm{\Δ}t({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{{\mathrm{\Δt\λ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\end{array}---\left(14\right)$

$\begin{array}{c}{\mathrm{\τ}}_{zz}^{t+}(i,j)={\mathrm{\τ}}_{zz}^{t-}(i,j)+\frac{\mathrm{\Δ}t\mathrm{\λ}}{\mathrm{\Δ}x}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t(\mathrm{\λ}+2\mathrm{\μ})}{\mathrm{\Δ}x}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\\ +\frac{{\mathrm{\Δt\λ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{-}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{-}\[v_z^t(i,j^{+})\]\}\end{array}---\left(15\right)$

$\begin{array}{c}{\mathrm{\τ}}_{xz}^{t+}(i,j)={\mathrm{\τ}}_{xx}^{t-}(i,j)+\frac{\mathrm{\Δ}t\mathrm{\μ}}{\mathrm{\Δ}x}\{D_z^{+}\[v_x^t(i^{+},j)\]\}+\frac{\mathrm{\Δ}t\mathrm{\μ}}{\mathrm{\Δ}x}\{D_x^{+}\[v_z^t(i,j^{+})\]\}\\ +\frac{{\mathrm{\Δt\μ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_z^{+}\[v_x^t(i^{+},j)\]\}+\frac{{\mathrm{\Δt\μ}}^{\′}}{\mathrm{\Δ}x}{D}_t^{-}\{D_x^{+}\[v_z^t(i,j^{+})\]\}\end{array}---\left(16\right)$

其中上标“+”表示向前差分，“-“表示向后差分；

(5)对于边界条件的处理，采用了完全匹配层边界条件

利用完全匹配层作为吸收边界是在区域ABCD四周引入完全匹配层，在区域四个边角区域1中，令d_x≠0；d_z≠0，速度V都等于角点的速度，在z方向区域2中，令d_x＝0；d_z≠0,速度V在z方向为常数，在x方向和边界的速度相等，在x方向区域3中，令d_x≠0；d_z＝0，速度V在x方向为常数，在z方向和边界的速度相等；

对于二维弹性波方程而言，这里V_x、V_z、σ_xx、σ_zz和σ_xz都进行分解，分解结果如下：

$V_x=V_x^x+V_x^z,V_z=V_z^x+V_z^z,$

${\mathrm{\σ}}_{xx}={\mathrm{\σ}}_{xx}^x+{\mathrm{\σ}}_{xx}^z,{\mathrm{\σ}}_{zz}={\mathrm{\σ}}_{zz}^x+{\mathrm{\σ}}_{zz}^z,{\mathrm{\σ}}_{xz}={\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z.$

二维弹性波PML公式：

$\left\{\begin{array}{c}\frac{\∂V_x^x}{\∂t}+d_x{V}_x^x=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xx}^x+{\mathrm{\σ}}_{xx}^z)}{\∂x}\\ \frac{\∂V_x^z}{\∂t}+d_z{V}_x^z=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z)}{\∂z}\\ \frac{\∂V_z^x}{\∂t}+d_x{V}_z^x=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{xz}^x+{\mathrm{\σ}}_{xz}^z)}{\∂x}\\ \frac{\∂V_z^z}{\∂t}+d_z{V}_z^z=\frac{1}{\mathrm{\ρ}}\frac{\∂({\mathrm{\σ}}_{zz}^x+{\mathrm{\σ}}_{zz}^z)}{\∂z}\end{array}\right.---\left(17\right)$

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\σ}}_{xx}^x}{\∂t}+d_x{\mathrm{\σ}}_{xx}^x=C_{11}\frac{\∂(V_x^x+V_x^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xx}^z}{\∂t}+d_z{\mathrm{\σ}}_{xx}^z=C_{13}\frac{\∂(V_z^x+V_z^z)}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{zz}^x}{\∂t}+d_x{\mathrm{\σ}}_{zz}^x=C_{13}\frac{\∂(V_x^x+V_x^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{zz}^z}{\∂t}+d_z{\mathrm{\σ}}_{zz}^z=C_{33}\frac{\∂(V_z^x+V_z^z)}{\∂z}\\ \frac{\∂{\mathrm{\σ}}_{xz}^x}{\∂t}+d_x{\mathrm{\σ}}_{xz}^x=C_{44}\frac{\∂(V_z^x+V_z^z)}{\∂x}\\ \frac{\∂{\mathrm{\σ}}_{xz}^z}{\∂t}+d_z{\mathrm{\σ}}_{xz}^z=C_{44}\frac{\∂(V_x^x+V_x^z)}{\∂z}\end{array}\right.---\left(18\right)$

(6)有限差分解法稳定性分析

2L阶精度的稳定性条件如式(19)，式中，Δt为时间间隔，Δx和Δz为网格大小，v_p为纵波速度，a_l为L个差分系数，

${\mathrm{\Δtv}}_p\sqrt{\frac{1}{{\mathrm{\Δx}}^2}+\frac{1}{{\mathrm{\Δz}}^2}}\≤\frac{1}{\underset{l=1}{\overset{L}{\Σ}}\left|a_l\right|}---\left(19\right)$

(7)进行观测系统设计，开展正演模拟。