CN101369024A - 一种地震波动方程生成方法及系统 - Google Patents

Abstract

本发明提供一种地震波动方程生成方法及系统，该地震波动方程生成方法包括：获取声波方程数据；获取地质参数信息；根据所述的声波方程数据和所述的地质参数信息进行任意差分精细积分，得出地震传播方程；求出稳定性条件、初始输入条件、边界处理条件；由所述的地震传播方程以及所述的稳定性条件、初始输入条件、边界处理条件生成地震波动方程。本发明有较高的精度和较好的数值稳定性，可以在几乎不增加计算量的情况下较大地提高了计算精度和稳定性。用任意差分精细积分法导出的完全匹配层边界条件对边界反射有很好的衰减作用。本发明可以为复杂区地震波传播规律研究提供了一个高精度、稳定性能好的数据体。

Description

一种地震波动方程生成方法及系统

技术领域

本发明有关于地质勘探技术领域，尤其是有关于地震勘探技术领域，具体地讲是一种地震波动方程生成方法及系统。

背景技术

在油气勘探中，能更好地、更准确地知道地质结构尤为重要，而知道了地质结构以后就可以判断出气田、油田的位置，以及地震传播的地震波动传播过程。

而目前，研究地震波动传播过程是通过求出地震波动方程来进行，即通过地震波场数值模拟技术来进行模拟。

上述的地震波场数值模拟技术是在地下介质结构和参数已知的情况下，利用理论计算的方法研究地震波在地下介质中的传播规律，并获得人工合成地震记录的一种技术。国内外在地震勘探生产实践中所采用的地震模型数值模拟方法有射线方法(射线追踪法)和波动方程数值解法(克希霍夫积分法、有限差分法、有限单元法、边界单元法和虚谱法)。

射线追踪的理论基础是在高频近似条件下，地震波场的主能量沿射线轨迹传播。主要问题在于：①难以处理介质中较强的速度变化；②难以求出多值走时中的全局最小走时；③复杂波场的振幅信息不准确；④阴影区内射线覆盖密度不足。

波动方程数值解法除了获得走时信息，还可得到波场的振幅信息，能描述在复杂介质中波的传播过程。

有限差分法是一种最常用的正演模拟方法，现已比较成熟，正向提高精度的方向发展。但差分法是从整体上考虑用一组差分方程来代替微分方程，对于复杂的结构形式及边界条件难以处理，对于高频信号分辨能力受到限制。

伪谱法做波动方程传播的正演已有一段时间，早在70年代初期就有人介绍该方法，它在气象预报、非线性波动和涡流模拟等领域得到广泛应用(最初由Kreiss和Oliger(1972)建议)。大约十年后，由Gazdag(1981)、Kosloff和Baysal(1982)最先应用于地震波动问题，虚谱法开始进入地震模拟方面。Kosloff等用虚谱法对声波方程进行了正演模拟，并与解析计算结果和超声物理模拟进行了比较，证明了方法的正确性。但是该方法处理吸收边界和自由表面边界条件比较困难，计算量和占用内存也较大。

有限元法在70年代将有限元法引入地震勘探领域，对其理论和应用做了大量研究和探索。有限元法可以从标量、矢量运动方程出发，不限于特定的坐标系，对整个区域进行灵活剖分，适用于均匀、非均匀地层所组成的复杂构造形态，但是当计算域较大时，用有限元法求解，因其空间坐标离散后未知数多，势必导致计算机耗时量大，计算成本高。目前对三维地震正演问题硬件条件尚不够。

目前三维正演模拟由于计算效率低，精度也不高，相对模型也比较简单，无法满足生产的需要。因此迫切需要有快速的、高精度的三维正演方法。

发明内容

本发明正是基于上述现有的问题而提出，其目的在于提供一种地震波动方程生成方法及系统，以快速、准确的求出地震波动方程。

为了实现上述目的，本发明的实施例提供一种地震波动方程生成方法，包括以下步骤：获取声波方程数据；获取地质参数信息；根据所述的声波方程数据和所述的地质参数信息进行任意差分精细积分，得出地震传播方程；求出稳定性条件、初始输入条件、边界处理条件；由所述的地震传播方程以及所述的稳定性条件、初始输入条件、边界处理条件生成地震波动方程。

为了实现上述目的，本发明的实施例还提供一种地震波动方程生成系统，该地震波动方程生成系统包括：声波方程获取单元，用于获取声波方程数据；地质参数获取单元，用于获取地质参数信息；传播方程生成单元，用于根据所述的声波方程数据和所述的地质参数信息进行任意差分精细积分，得出地震传播方程；条件生成单元，用于求出稳定性条件、初始输入条件、边界处理条件；波动方程生成单元，用于由所述的地震传播方程以及所述的稳定性条件、初始输入条件、边界处理条件生成地震波动方程。

本发明的任意差分精细积分方法有较高的精度和较好的数值稳定性，可以在几乎不增加计算量的情况下较大地提高了计算精度和稳定性。完全匹配边界条件可以用于波场边界上进行振幅衰减。用任意差分精细积分法导出的完全匹配层边界条件对边界反射有很好的衰减作用。本发明可以为复杂区地震波传播规律研究提供了一个高精度、稳定性能好的数据体。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：

图1A所示的是本发明的实施方式的地震波动方程生成系统的结构框图；

图1B所示的是本发明的实施方式的地质参数获取单元的结构框图；

图1C所示的是本发明的实施方式的传播方程生成单元的结构框图；

图2所示的是本发明的实施方式的地震波动方程生成方法的流程图；

图3所示的是本发明的实施方式的进行任意差分精细积分的结果示意图；

图4所示的是本发明的实施方式的进行任意差分精细积分的结果示意图；

图5所示的是本发明的实施方式的一维波动方程完全匹配边界条件应用效果对比图；

图6所示的是本发明的实施方式的完全匹配层吸收边界示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。

本发明实施例提供一种地震波动方程生成方法及系统。以下结合附图对本发明进行详细说明。

图1A所示的是本发明的实施方式的地震波动方程生成系统的结构框图、图1B所示的是本发明的实施方式的地质参数获取单元的结构框图、图1C所示的是本发明的实施方式的传播方程生成单元的结构框图。

其中，本发明的实施方式的地震波动方程生成系统包括：声波方程获取单元101，用于获取声波方程数据；地质参数获取单元102，用于获取地质参数信息；传播方程生成单元103，用于根据所述的声波方程数据和所述的地质参数信息进行任意差分精细积分，得出地震传播方程；条件生成单元104，用于求出稳定性条件、初始输入条件、边界处理条件；波动方程生成单元105，用于由所述的地震传播方程以及所述的稳定性条件、初始输入条件、边界处理条件生成地震波动方程。

上述地质参数获取单元包括：地震勘探模块106，采集地震数据；数据处理模块107，利用promax软件处理所述地震数据得出地震剖面；分析模块108，对所述地震剖面进行分析得出对应的地质参数信息。

而上述传播方程生成单元包括：有限差分模块109，用于进行时间有限差分；精细积分模块110，用于进行空间精细积分，所述传播方程生成单元进行时间有限差分和空间精细积分的交替来进行所述的任意差分精细积分。

另外，上述边界条件是指完全匹配层吸收边界条件，以在完全匹配边界条件下进行任意差分精细积分，并且求出的地震波动方程为一半解析解。

图2所示的是本发明的实施方式的地震波动方程生成方法的流程图，一下结合图2详细说明本发明的地震波动方程生成方法。

S201：获取声波方程数据；

S202：采集地震数据；

S203：利用promax软件处理所述地震数据得出地震剖面；

S204：对所述地震剖面进行分析得出对应的地质参数信息；

S205：进行时间有限差分；

S206：交替进行空间精细积分得出地震传播方程；

S207：求出稳定性条件、初始输入条件、边界处理条件；

S208：由所述的地震传播方程以及所述的稳定性条件、初始输入条件、边界处理条件生成地震波动方程。

上述边界条件是指完全匹配层吸收边界条件，以在完全匹配边界条件下进行任意差分精细积分，且上述地震波动方程为一半解析解。

以下详细说明本发明解决其技术问题所采用的技术手段：

一维波动方程任意差分精细积分原理：

本节以一维波动方程的初边值问题为例来阐述任意差分精细积分法的基本原理。

设有如下初边值问题：

$\left\{\begin{array}{c}\frac{{\∂}^{2}u}{\∂t^2}=v^2\frac{{\∂}^{2}u}{\∂x^2}\\ u(x,0)=f_1\left(x\right)\\ \frac{\∂u(x,0)}{\∂t}=f_2\left(x\right)\\ u(0,t)=g_1\left(t\right),u(L,t)=g_2\left(t\right)\end{array}\right.---(1-1)$

式中，u代表位移，t是时间变量，x是空间坐标，v是波的传播速度，L是一维区域的长度。

对(1-1)中的第一式进行空间坐标离散，在任一点j(点序号)处，有：

$\begin{array}{cc}\frac{{\∂}^2{u}_j}{\∂t^2}=v^2\frac{{\∂}^2{u}_j}{\∂x^2}{|}_{x=x_j}& (0\≤j\≤J)\end{array}---(1-2)$

这里，J是总离散点数。将函数u(x，t)在j点的邻域内作泰勒展开，得到：

$u_i=u_j+\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\frac{1}{k!}{\left(\mathrm{\Δ}{x}_i\frac{\∂}{\∂x}\right)}^k{u|}_{x=x_j}+O\left[{\left(\mathrm{\Δ}{x}_i\right)}^{m+1}\right]---(1-3)$

其中，Δx_i＝x_i-x_j，i是j点邻域内的任意点，对(1-3)式移项并取n个加权参数后再求和得：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i(u_i-u_j)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i(\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\frac{1}{k!}{\left(\mathrm{\Δ}{x}_i\frac{\∂}{\∂x}\right)}^k{u|}_{x=x_j}+O\left[{\left(\mathrm{\Δ}{x}_i\right)}^{m+1}\right])---(1-4)$

其中α_i是加权系数，取m＝n并省略掉(1-4)式右端项的高阶项O[(Δx_i)^m+1]，再令其等于(1-2)式的右端，得到：

$\left\{\begin{array}{cc}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i\mathrm{\Δ}{x}_i=0,& \underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\left(\mathrm{\Δ}{x}_i\right)}^2=2v^2\\ \underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\left(\mathrm{\Δ}{x}_i\right)}^3=0,& \underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\left(\mathrm{\Δ}{x}_i\right)}^4=0\\ \·\·\·\·\·\·\·\·\·\·\·& \·\·\·\·\·\·\·\·\·\·\·\\ \underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\left(\mathrm{\Δ}{x}_i\right)}^{m-1}=0,& \underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\left(\mathrm{\Δ}{x}_i\right)}^m=0,\end{array}\right.---(1-5)$

上式是一个关于α_i的m维线性方程组，解(1-5)可求出加权系数α_i。

(对于均匀网格，m＝2的情形，有 ${\mathrm{\α}}_i={\left(\frac{v}{\mathrm{\Δx}}\right)}^2$ )。

(1-2)式可近似为：

$\frac{{\∂}^2{u}_j}{\∂t^2}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i(u_i-u_j)---(1-6)$

在t_n时刻附近的小领域内，(1-6)可表示为下式：

$\frac{{\∂}^2{u}_j}{\∂t^2}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i(u_{i,n}-u_j)---(1-7)$

(上式u的第一个下标代表空间变量，第二个下标代表时间变量，以下同)

(1-7)式右端的u_j移到方程的左端得t_n时刻小领域内常微分方程：

$\frac{{\∂}^2{u}_j}{\∂t^2}+u_j\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_{i,n}---(1-8)$

在t_n时刻小领域内，将上式右端看成是不随时间变化的常数项，其解析解为：

$u_j=c_1\exp \left(\mathrm{ia}(t-t_n)\right)+c_2\exp (-\mathrm{ia}(t-t_n))+\frac{b_n}{a^2}---(1-9)$

式中， $a^2=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i,b_n=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_{i,n}$

可以看出(1-9)式有两个未知数c₁和c₂，选择在子域[t_n-1，t_n+1]积分，假设t_n-1，t_n时刻的值已知，求出t_n+1时刻的递推公式：

$\left\{\begin{array}{c}{u}_{j,n+1}=2u_{j,n}\cos \left(\mathrm{a\Δt}\right)-u_{j,n-1}+2\frac{b_n}{a^2}[1-\cos \left(\mathrm{a\Δt}\right)]\\ u_{j,1}=u_{j,0}\cos \left(\mathrm{a\Δt}\right)+\frac{b_0}{a^2}[1-\cos \left(\mathrm{a\Δt}\right)]+\mathrm{\Δt}{f}_2\left(x_j\right)\\ u_{j,0}=f_1\left(x_j\right)\end{array}\right.---(1-10)$

为了进一步提高差分格式的精度，对(1-6)移项：

$\frac{{\∂}^2{u}_j}{\∂t^2}+u_j\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_i---(1-11)$

将(1-11)式的右端项在t_n的邻域内做泰勒展开

$\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_i=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_{i,n}+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\stackrel{\·}{u}}_{i,n}(t-t_n)+\frac{1}{2}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\stackrel{\·\·}{u}}_{i,n}{(t-t_n)}^2+O\left({(t-t_n)}^3\right)$

根据(1-11)式的右端项不同的截断，可以得到不同的计算格式，展开到t-t_n的二阶导数，得：

$\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_i=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_{i,n}+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\stackrel{\·}{u}}_{i,n}(t-t_n)+\frac{1}{2}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\stackrel{\·\·}{u}}_{i,n}{(t-t_n)}^2---(1-12)$

解常微分方程：

$\frac{{\∂}^2{u}_j}{{\∂t}^2}+u_j\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_{i,n}+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\stackrel{\·}{u}}_{i,n}(t-t_n)+\frac{1}{2}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\stackrel{\·\·}{u}}_{i,n}{(t-t_n)}^2---(1-13)$

得解析解(此即精细积分的意义)为：

$u_j=c_1\exp \left(\mathrm{ia}(t-t_n)\right)+c_2\exp (-\mathrm{ia}(t-t_n))+\frac{b_n}{a^2}$

$+\frac{e_n(t-t_n)}{a^2}+\frac{d_n{(t-t_n)}^2}{a^2}-\frac{2d_n}{a^4}---(1-14)$

上式中， $a^2=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i,b_n=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{u}_{i,n},$ $e_n=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\stackrel{\·}{u}}_{i,n},d_n=\frac{1}{2}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\stackrel{\·\·}{u}}_{i,n},,$ 对(1-14)式选择积分子域为[t_n-1，t_n+1]，利用t_n-1，t_n时刻的值，确定两个未知数c₁和c₂，t_n+1时刻的递推计算公式为：

$\left\{\begin{array}{c}{u}_{j,n+1}=2u_{j,n}\cos \left(\mathrm{a\Δt}\right)-u_{j,n-1}+2(\frac{b_n}{a^2}-\frac{2d_n}{a^4})[1-\cos \left(\mathrm{a\Δt}\right)]+\frac{2d_n}{a^2}{\left(\mathrm{\Δt}\right)}^2\\ u_{j,1}=u_{j,0}\cos \left(\mathrm{a\Δt}\right)+(\frac{b_0}{a^2}-\frac{2d_0}{a^4})[1-\cos \left(\mathrm{a\Δt}\right)]+\frac{d_0}{a^2}{\left(\mathrm{\Δt}\right)}^2+\mathrm{\Δt}{f}_2\left(x_j\right)\\ u_{j,0}=f_1\left(x_j\right)\end{array}\right.---(1-15)$

上式中a，b_n同(1-9)式，而d_n待求，下面求d_n

由于

的计算可由(1-6)式求得： $\frac{{\∂}^2{u}_i}{\∂t^2}=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k(u_k-u_i),$ 即： ${\stackrel{\·\·}{u}}_{i,n}=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k(u_{k,n}-u_{i,n})=b_{i,n}-a_i^2{u}_{i,n}$ 求得 $d_n=\frac{1}{2}\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i(b_{i,n}-a_i^2{u}_{i,n})\right]$ (1-16)

同理可证，当(1-11)式的右端项进行t-t_n的零阶展开和t-t_n的一阶展开时，得到的递推公式相同，都为：

$\left\{\begin{array}{c}{u}_{j,n+1}=2u_{j,n}\cos \left(\mathrm{a\Δt}\right)-u_{j,n-1}+2\frac{b_n}{a^2}[1-\cos \left(\mathrm{a\Δt}\right)]\\ u_{j,1}=u_{j,0}\cos \left(\mathrm{a\Δt}\right)+\frac{b_0}{a^2}[1-\cos \left(\mathrm{a\Δt}\right)]+\mathrm{\Δt}{f}_2\left(x_j\right)\\ u_{j,0}=f_1\left(x_j\right)\end{array}\right.---(1-17)$

由于(1-15)是时间二阶展开，其截断误差较小，因此理论上，(1-15)式较(1-10)具有更高的时间精度。

二维及三维波动方程的精细积分方法的导出：

精细积分法是采用差分的方法计算空间偏导数，采用子域精细积分的方法计算时间偏导数。对于二维波动方程：

$\left\{\begin{array}{c}\frac{{\∂}^{2}u}{\∂t^2}=v^2(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂z}^2})+\mathrm{\δ}(x-x_0,z-z_0)f\left(t\right)\\ u(x,z,0)=0\\ \frac{\∂u(x,z,0)}{\∂t}=0\end{array}\right.---(2-1)$

首先对上式中的第一式进行二维空间坐标离散(按均匀网格划分)，在任一离散点j(点序号)处，有：

$\frac{{\∂}^2{u}_j}{{\∂t}^2}=v^2{(\frac{{\∂}^2{u}_j}{\∂x^2}+\frac{{\∂}^2{u}_j}{{\∂z}^2})|}_{(x,z)=(x_j,z_j)}+\mathrm{\δ}(x_j-x_0,z_j-z_0)f\left(t\right)---(2-2)$

上式中，(0≤j≤J)，J是总离散点数。

对(2-2)式u在j点的x邻域作泰勒展开：

$u_i=u_j+\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\frac{1}{k!}{\left(\mathrm{\Δ}{x}_i\frac{\∂}{\∂x}\right)}^k{u|}_{x=x_j}+O\left[{\left(\mathrm{\Δ}{x}_i\right)}^{m+1}\right]$

上式中，Δx_i＝x_i-x_j，i是j的x邻域内的任意点，对上式移项并加权得：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i(u_i-u_j)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i(\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\frac{1}{k!}{\left(\mathrm{\Δ}{x}_i\frac{\∂}{\∂x}\right)}^k{u|}_{x=x_j}+O\left[{\left(\mathrm{\Δ}{x}_i\right)}^{m+1}\right])---(2-3)$

其中α_i是加权系数，m原则上可任意选取，但m>n使下式的方程组成为矛盾方程组，m<n使下式的方程组为欠定的，在此一般取m＝n(以m＝n为例来展开讨论)，

省略掉(2-3)式右端项的高阶项O[(Δx_i)^m+1]，再令其等于(2-2)式的右端项

则有：

$\left\{\begin{array}{cc}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\mathrm{\&Delta;}{x}_i=0,& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^2=2v^2\\ \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^3=0,& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^4=0\\ \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\\ \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^{m-1}=0,& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^m=0,\end{array}\right.$

上式是一个关于α_i的m维线性方程组，解之可求出加权系数α_i。

令： $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(u_i-u_j)=v^2\frac{{\&PartialD;}^2{u}_j}{\&PartialD;x^2}.$

同理，

也可以这样处理，将u在j点的z邻域作泰勒展开：

$u_p=u_j+\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{z}_p\frac{\&PartialD;}{\&PartialD;z}\right)}^k{u|}_{z=z_j}+O\left[{\left(\mathrm{\&Delta;}{z}_p\right)}^{m+1}\right]$

上式中，Δz_p＝z_p-z_j，p是j的z邻域内的任意点，对上式移项并加权得：

$\underset{p=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p(u_p-u_j)=\underset{p=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p(\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{z}_p\frac{\&PartialD;}{\&PartialD;z}\right)}^k{u|}_{z=z_j}+O\left[{\left(\mathrm{\&Delta;}{z}_p\right)}^{m+1}\right])---(2-4)$

其中α_p是加权系数，m原则上可任意选取，但m>n使下式的方程组成为矛盾方程组，m<n使下式的方程组为欠定的，在此一般取m＝n(以m＝n为例来展开讨论)，省略掉(2-4)式右端项的高阶项O[(Δz_p)^m+1])，再令其等于(2-2)式的右端项

则有：

$\left\{\begin{array}{cc}\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p\mathrm{\&Delta;}{z}_p=0,& \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^2=2v^2\\ \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^3=0,& \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^4=0\\ \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\\ \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^{m-1}=0,& \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^m=0,\end{array}\right.$

上式是一个关于α_p的m维线性方程组，解之可求出加权系数α_p。 $\underset{p=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p(u_p-u_j)=v^2\frac{{\&PartialD;}^2{u}_j}{\&PartialD;z^2},$ 对α_p，α_i重新排列，并统一用α_i表示。(2-2)右端的二阶导数可用相邻点的差分来表示。暂不考虑震源项δ(x-x₀，z-z₀)f(t)，那么(2-2)式可表示为：

$\frac{{\&PartialD;}^2{u}_i}{\&PartialD;t^2}=\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(u_i-u_j)$

该式的右端可近似为t_n时刻，关于u_j微分方程，因此上式可写为：

$\frac{{\&PartialD;}^2{u}_j}{{\&PartialD;t}^2}+u_j\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i=\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_i---(2-5)$

在t_n时刻的小邻域内(2-5)的右端项可展开到t-t_n的二阶项来近似，即：

$\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_i=\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_{i,n}+\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;}{u}}_{i,n}(t-t_n)+\frac{1}{2}\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;\&CenterDot;}{u}}_{i,n}{(t-t_n)}^2$

将上式代入方程(2-5)，求得的解析解为：

$u_j=c_1\exp \left(\mathrm{ia}(t-t_n)\right)+c_2\exp (-\mathrm{ia}(t-t_n))+\frac{b_n}{a^2}+\frac{e_n(t-t_n)}{a^2}+$

$\frac{d_n{(t-t_n)}^2}{a^2}-\frac{2d_n}{a^4}$

式中， $a^2=\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i,b_n=\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_{i,n},e_n=\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;}{u}}_{i,n},d_n=\frac{1}{2}\underset{i=1}{\overset{2*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;\&CenterDot;}{u}}_{i,n}$

选择积分子域[t_n-1，t_n+1]，假定t_n-1，t_n的波场已知，可确定c₁和c₂，进而可以求出t_n+1时刻的值：

$\left\{\begin{array}{c}{u}_{j,n+1}=2u_{j,n}\cos \left(\mathrm{a\&Delta;t}\right)-u_{j,n-1}+2(\frac{b_n}{a^2}-\frac{2d_n}{a^4})[1-\cos \left(\mathrm{a\&Delta;t}\right)]+\frac{2d_n}{a^2}{\left(\mathrm{\&Delta;t}\right)}^2\\ u_{j,1}=u_{j,0}\cos \left(\mathrm{\&alpha;\&Delta;t}\right)+(\frac{b_0}{a^2}-\frac{2d_0}{a^4})[1-\cos \left(\mathrm{a\&Delta;t}\right)]+\frac{d_0}{a^2}{\left(\mathrm{\&Delta;t}\right)}^2\\ u_{j,0}=0\end{array}\right.---(2-6)$

对(2-5)式的右端项进行t-t_n的零阶展开和t-t_n的一阶展开时，得到的递推公式相同，都为：

$\left\{\begin{array}{c}{u}_{j,n+1}=2u_{j,n}\cos \left(\mathrm{a\&Delta;t}\right)-u_{j,n-1}+2\frac{b_n}{a^2}[1-\cos \left(\mathrm{a\&Delta;t}\right)]\\ u_{j,1}=u_{j,0}\cos \left(\mathrm{a\&Delta;t}\right)+\frac{b_0}{a^2}[1-\cos \left(\mathrm{a\&Delta;t}\right)]\\ u_{j,0}=0\end{array}\right.---(2-7)$

从(2-15)和(2-6)式可以看出，一维与二维波动方程的ADPI(任意差分精细积分法)方法计算公式形式上完全相同，但实际上公式中的b_n和d_n是与空间维数密切相关的。

对于三维波动方程：

$\left\{\begin{array}{c}\frac{{\&PartialD;}^{2}u}{\&PartialD;t^2}=v^2(\frac{{\&PartialD;}^{2}u}{{\&PartialD;x}^2}+\frac{{\&PartialD;}^{2}u}{\&PartialD;y^2}+\frac{{\&PartialD;}^{2}u}{{\&PartialD;z}^2})+\mathrm{\&delta;}(x-x_0,y-y_0,z-z_0)f\left(t\right)\\ u(x,y,z,0)=0\\ \frac{\&PartialD;u(x,y,z,0)}{\&PartialD;t}=0\end{array}\right.---(2-8)$

具体做法为：

首先对上式中的第一式进行三维空间坐标离散(即按均匀网格划分)，在任一离散点j(点序号)处，有：

$\frac{{\&PartialD;}^2{u}_j}{{\&PartialD;t}^2}=v^2{(\frac{{\&PartialD;}^2{u}_j}{\&PartialD;x^2}+\frac{{\&PartialD;}^2{u}_j}{\&PartialD;y^2}+\frac{{\&PartialD;}^2{u}_j}{{\&PartialD;z}^2})|}_{(x,y,z)=(x_j,y_j,z_j)}+\mathrm{\&delta;}(x_j-x_0,y_j-y_0,z_j-z_0)f\left(t\right)---(2-9)$

上式中，(0≤j≤J)，J是总离散点数

(2-9)式的右端项

可用任意差分的方法求得，将u在j点的x邻域作泰勒展开：

$u_i=u_j+\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{x}_i\frac{\&PartialD;}{\&PartialD;x}\right)}^k{u|}_{x=x_j}+O\left[{\left(\mathrm{\&Delta;}{x}_i\right)}^{m+1}\right]$

上式中，Δx_i＝x_i-x_j，i是j的x邻域内的任意点，对上式移项并加权得：

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(u_i-u_j)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{x}_i\frac{\&PartialD;}{\&PartialD;x}\right)}^k{u|}_{x=x_j}+O\left[{\left(\mathrm{\&Delta;}{x}_i\right)}^{m+1}\right])---(2-10)$

其中α_i是加权系数，一般取m＝n，省略掉(2-3)式右端项的高阶项O[(Δx_i)^m+1]，再令其等于(2-2)式的右端项

则有：

$\left\{\begin{array}{cc}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\mathrm{\&Delta;}{x}_i=0,& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^2=2v^2\\ \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^3=0,& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^4=0\\ \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\\ \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^{m-1}=0,& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^m=0,\end{array}\right.$

上式是一个关于α_i的m维线性方程组，解之可求出加权系数α_i。

令： $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(u_i-u_j)=\frac{{\&PartialD;}^2{u}_j}{\&PartialD;x^2}{v}^2$

同理，

也可以这样处理，将u在j点的z邻域作Taylor展开：

$u_p=u_j+\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{z}_p\frac{\&PartialD;}{\&PartialD;z}\right)}^k{u|}_{z=z_j}+O\left[{\left(\mathrm{\&Delta;}{z}_p\right)}^{m+1}\right]$

上式中，Δz_p＝z_p-z_j，p是j的z邻域内的任意点，对上式移项并加权得：

$\underset{p=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p(u_p-u_j)=\underset{p=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p(\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{z}_p\frac{\&PartialD;}{\&PartialD;z}\right)}^k{u|}_{z=z_j}+O\left[{\left(\mathrm{\&Delta;}{z}_p\right)}^{m+1}\right])---(2-11)$

取m＝n，省略掉(2-11)式右端项的高阶项O[(Δz_p)^m+1])，再令其等于(2-9)式的右端项

则有：

$\left\{\begin{array}{cc}\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p\mathrm{\&Delta;}{z}_p=0,& \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^2=2v^2\\ \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^3=0,& \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^4=0\\ \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\\ \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^{m-1}=0,& \underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p{\left(\mathrm{\&Delta;}{z}_p\right)}^m=0,\end{array}\right.$

上式是一个关于α_p的m维线性方程组，解之可求出加权系数α_p。

令： $\underset{p=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_p(u_p-u_j)=\frac{{\&PartialD;}^2{u}_j}{\&PartialD;z^2}{v}^2.$

对做与

相同的处理，将u在j点的y邻域作泰勒展开：

$u_q=u_j+\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{y}_q\frac{\&PartialD;}{\&PartialD;y}\right)}^k{u|}_{y=y_j}+O\left[{\left(\mathrm{\&Delta;}{y}_p\right)}^{m+1}\right]$

上式中，Δy_q＝y_q-z_j，q是j的y邻域内的任意点，对上式移项并加权得：

$\underset{q=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q(u_q-u_j)=\underset{q=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q(\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{y}_q\frac{\&PartialD;}{\&PartialD;y}\right)}^k{u|}_{y=y_j}+O\left[{\left(\mathrm{\&Delta;}{y}_q\right)}^{m+1}\right])---(2-12)$

取m＝n，省略掉(2-12)式右端项的高阶项O[(Δy_q)^m+1]，再令其等于(2-9)式的右端项则有：

$\left\{\begin{array}{cc}\underset{q=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q\mathrm{\&Delta;}{y}_q=0,& \underset{q=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q{\left(\mathrm{\&Delta;}{y}_q\right)}^2=2v^2\\ \underset{q=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q{\left(\mathrm{\&Delta;}{y}_q\right)}^3=0,& \underset{q=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q{\left(\mathrm{\&Delta;}{y}_q\right)}^4=0\\ \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\\ \underset{q=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q{\left(\mathrm{\&Delta;}{y}_q\right)}^{m-1}=0,& \underset{q=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q{\left(\mathrm{\&Delta;}{y}_q\right)}^m=0\end{array}\right.$

上式是一个关于α_q的m维线性方程组，解之可求出加权系数α_q。

令： $\underset{q=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_q(u_q-u_j)=\frac{{\&PartialD;}^2{u}_j}{\&PartialD;y^2}{v}^2$

对α_p，α_q，α_i重新排列，并统一用α_i表示。则α_i是大小为3*m的系数。

(2-2)右端的二阶导数可用相邻点的任意差分来表示，暂不考虑震源，那么可表示为：

$\frac{{\&PartialD;}^2{u}_j}{\&PartialD;t^2}=\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(u_i-u_j)$

该式的右端可近似为t_n时刻，关于u_j微分方程，因此上式可写为：

$\frac{{\&PartialD;}^2{u}_j}{{\&PartialD;t}^2}+u_j\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i=\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_i---(2-13)$

将(2-13)的右端项展开到t-t_n的二阶，即：

$\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_i=\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_{i,n}+\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;}{u}}_{i,n}(t-t_n)+\frac{1}{2}\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;\&CenterDot;}{u}}_{i,n}{(t-t_n)}^2$

这时的解析解为：

$u_j=c_1\exp \left(\mathrm{ia}(t-t_n)\right)+c_2\exp (-\mathrm{ia}(t-t_n))+\frac{b_n}{a^2}+\frac{e_n(t-t_n)}{a^2}+$

$\frac{d_n{(t-t_n)}^2}{a^2}-\frac{2d_n}{a^4}$

上式中， $a^2=\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i,b_n=\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_{i,n},e_n=\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;}{u}}_{i,n},d_n=\frac{1}{2}\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;\&CenterDot;}{u}}_{i,n}$

选择积分子域[t_n-1，t_n+1]，假定t_n-1，t_n的波场已知，可确定c₁和c₂，进而可以求出t_n+1时刻的值：

$\left\{\begin{array}{c}{u}_{j,n+1}=2u_{j,n}\cos \left(\mathrm{a\&Delta;t}\right)-u_{j,n-1}+2(\frac{b_n}{a^2}-\frac{2d_n}{a^4})[1-\cos \left(\mathrm{a\&Delta;t}\right)]+\frac{2d_n}{a^2}{\left(\mathrm{\&Delta;t}\right)}^2\\ u_{j,1}=u_{j,0}\cos \left(\mathrm{a\&Delta;t}\right)+(\frac{b_0}{a^2}-\frac{2d_0}{a^4})[1-\cos \left(\mathrm{a\&Delta;t}\right)]+\frac{d_0}{a^2}{\left(\mathrm{\&Delta;t}\right)}^2\\ u_{j,0}=0\end{array}\right.---(2-14)$

上式中， $a^2=\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i,b_n=\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{u}_{i,n},d_n=\frac{1}{2}\underset{i=1}{\overset{3*m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\stackrel{\&CenterDot;\&CenterDot;}{u}}_{i,n}.$

当对(2-6)式的右端项进行t-t_n的零阶展开和t-t_n的一阶展开时，得到的递推公式相同，都为

$\left\{\begin{array}{c}{u}_{j,n+1}=2u_{j,n}\cos \left(\mathrm{a\&Delta;t}\right)-u_{j,n-1}+2\frac{b_n}{a^2}[1-\cos \left(\mathrm{a\&Delta;t}\right)]\\ u_{j,1}=u_{j,0}\cos \left(\mathrm{a\&Delta;t}\right)+\frac{b_0}{a^2}[1-\cos \left(\mathrm{a\&Delta;t}\right)]\\ u_{j,0}=0\end{array}\right.---(2-15)$

运用公式(2-14)进行数值求解时以下几个问题必须注意：

①对

项，公式(2-14)在邻域[t_n-1，t_n+1]内二阶展开，随着展开的阶数增高，数值求解的精度也会相应提高。当然计算公式的复杂性和计算花费的时间也会增加。

②d_n的计算涉及到求

用预报--校正方法可提高计算精度，首先用(2-15)式预测u_j，n+1，用下列中心差分计算公式预测

和

$u_{j,n+1}^{\&prime;}=(u_{j,n+1}-u_{j,n-1})/\left(2\mathrm{\&Delta;t}\right)$

${\stackrel{\&CenterDot;\&CenterDot;}{u}}_{j,n}=(u_{j,n+1}^{\&prime;}-2\&times;u_{j,n}+u_{j,n-1})/\left({\mathrm{\&Delta;}}^{2}t\right)---(2-16)$

$d_n=1/2\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}(u_{j,n+1}^{\&prime;}-2*u_{j,n}+u_{j,n-1})/\left({\mathrm{\&Delta;}}^{2}t\right)---(2-17)$

(2-7)的第一个公式简化为校正公式：

$u_{j,n+1}=u_{j,n+1}^{\&prime;}+\frac{2d_n}{a^2}{\left(\mathrm{\&Delta;t}\right)}^2-\frac{4d_n}{a^4}[1-\cos \left(\mathrm{a\&Delta;t}\right)]---(2-18)$

预报一校正方法比较实用，精度较高，且节约了计算量，本文使用了这种方法。

③要使得ADPI方法成为实用的地震正演技术，一方面，需要用理论解析解作为参照，证明方法的正确性；另一方面，必须考虑如何确定将实际半无穷空间波场转化为有限区域波场的边界条件。同时要研究方法的数值计算稳定性，并与实际科研生产中正在使用的其它地震正演技术对比，说明ADPI方法的精度。

围绕问题③展开ADPI方法的进一步研究，下面的内容将逐一回答上述的计算公式的正确性、边界条件、稳定性、与数值计算的精确度问题。

ADPI方法稳定性分析：

采用Von Neumann方法来分析(1-15)式的数值稳定性。首先介绍下面的引理。

引理：实系数二次方程μ²-bμ-c＝0的根，按其模小于或等于1的充分且必要条件是|b|≤1-c，|c|≤1。

为了便于与有限差分法比较，这里考虑均匀网格、当m＝2时，可以得到：

$\left|\mathrm{\&Delta;}{x}_i\right|=\mathrm{\&Delta;x},{\mathrm{\&alpha;}}_i={\left(\frac{v}{\mathrm{\&Delta;x}}\right)}^2=\mathrm{\&alpha;},b_n=\mathrm{\&alpha;}(u_{j-1,n}+u_{j+1,n}),a^2=2\mathrm{\&alpha;}$

令u_j，n＝Vⁿe^ikjΔx，代入(1-15)式，得到下式：

$V^{n+1}{e}^{\mathrm{ikj\&Delta;x}}=2\cos \left(\mathrm{a\&Delta;t}\right)V^n{e}^{\mathrm{ikj\&Delta;x}}+\frac{2\mathrm{\&alpha;}(V^n{e}^{\mathrm{ik}(j-1)\mathrm{\&Delta;x}}+V^n{e}^{\mathrm{ik}(j+1)\mathrm{\&Delta;x}})(1-\cos \left(\mathrm{a\&Delta;t}\right))}{a^2}-V^{n-1}{e}^{\mathrm{ikj\&Delta;x}}$

根据欧拉公式e^ikΔx＝cos(kΔx)+i sin(kΔx)，消去公因子e^ikjΔx得

其中λ＝cos(aΔt)，β＝cos(kΔx)。由于(4-1)式是三层格式，Vⁿ⁺¹与Vⁿ、V^n-1有关，故令：

Vⁿ＝1×Vⁿ+0×V^n-1    (4-2)

将(4-1)式和(4-2)式合并写成

$\left[\begin{array}{c}{V}^{n+1}\\ V^n\end{array}\right]=\left[\begin{array}{cc}2\mathrm{\&lambda;}+2(2\mathrm{\&beta;}-{\mathrm{\&beta;}}^2)(1-\mathrm{\&lambda;})+2\mathrm{\&alpha;\&beta;}(\mathrm{\&beta;}-1){\left(\mathrm{\&Delta;t}\right)}^2& -1\\ 1& 0\end{array}\right]\left[\begin{array}{c}{V}^n\\ V^{n-1}\end{array}\right]---(4-3)$

式(4-3)中过渡矩阵为：

$G=\left[\begin{array}{cc}2\mathrm{\&lambda;}+2(2\mathrm{\&beta;}-{\mathrm{\&beta;}}^2)(1-\mathrm{\&lambda;})+2\mathrm{\&alpha;\&beta;}(\mathrm{\&beta;}-1){\left(\mathrm{\&Delta;t}\right)}^2& -1\\ 1& 0\end{array}\right]$

G的特征值μ满足：

μ²-bμ-c＝0

其中b＝2λ+2(2β-β²)(1-λ)+2αβ(β-1)(Δt)²，c＝-1

当 $0\&le;\mathrm{a\&Delta;t}\&le;\frac{\mathrm{\&pi;}}{2}$ 时，(aΔt)²≈2(1-cos(aΔt))

b≈2λ+2β(1-λ)

${\left(\frac{b}{2}\right)}^2\&ap;{\mathrm{\&lambda;}}^2+{\mathrm{\&beta;}}^2{(1-\mathrm{\&lambda;})}^2+2\mathrm{\&lambda;\&beta;}(1-\mathrm{\&lambda;})\&le;{\mathrm{\&lambda;}}^2+{(1-\mathrm{\&lambda;})}^2+2\mathrm{\&lambda;}(1-\mathrm{\&lambda;})$

$={(\mathrm{\&lambda;}+1-\mathrm{\&lambda;})}^2=1$

所以|b|≤2＝1-c，|c|＝1

由引理知，当 $0\&le;\mathrm{a\&Delta;t}\&le;\frac{\mathrm{\&pi;}}{2}$ 时，G的特征值小于等于1。

定理1：由Von Neumann条件知(1-15)式稳定。即当 $0\&le;\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;x}}\&le;\frac{\mathrm{\&pi;}}{2\sqrt{2}}=1.11089$ 时，四点二阶精细积分波动方程解法稳定。

对零阶展开和一阶展开得到的计算格式(1-17)，稳定性条件仍为 $0\&le;\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;x}}\&le;\frac{\mathrm{\&pi;}}{2\sqrt{2}}=1.11089.$

对于(1-1)式的有限差分格式，由参考文献可知其稳定性条件为 $0\&le;\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;x}}\&le;1,$ 对比可知精细积分法比有限差分法更为稳定。

下面用数值求解的例子直观地说明了ADPI方法的稳定性比有限差分格式更好。图3表明，当dt＝0.00025s， $\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;x}}=1.228$ 时，有限差分法不稳定，结果发散。而ADPI方法稳定，结果正确；图4表明，当dt＝0.00075s， $\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;x}}=1.165$ 有限差分法不稳定，结果发散。而ADPI方法稳定，结果正确。

例1：dt＝0.00025s， $\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;x}}=1.228.$

其中，图3的左图有限差分法的结果，不稳定，解发散；

图3的右图二阶ADPI法的结果，稳定，解正确。

例2：dt＝0.00075s， $\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;x}}=1.165.$

其中，图4的左图有限差分法的结果，不稳定，解发散

图4的右图二阶ADPI法的结果，稳定，解正确

边界条件分析：

利用波动方程模拟地震记录时，吸收边界条件是非常重要的。在一个有限的区域中，利用数值方法进行计算时，如果不对边界进行处理，就会产生边界反射。现有许多方法用以消除边界反射。Cerjan等人通过在计算边界附近引入损耗介质来衰减向外传播的波。由于在两个具有不同吸收系数的损耗介质的界面上会产生反射，这就要求有较多的吸收层，才能有较好的吸收效果，而较多的吸收层意味着较大的计算量。另一种比较著名的吸收边界条件是Clayton吸收边界条件，该吸收边界条件在旁轴近似理论的基础上导出，在特定的入射角和频率范围内，具有较好的吸收效果。Berenger针对电磁波传播情况，给出了一种高效的完全匹配层吸收边界条件，并在理论上证明该方法可以完全吸收来自各个方向、各种频率的电磁波，而不产生任何反射。这里将完全匹配层的思想引入到三维波动方程ADPI解法中，并推导了三维波动方程ADPI解法的完全匹配层吸收边界条件。

①一维波动方程完全匹配层原理

在一维声波介质中，波的传播可以用如下方程来描述：

$\frac{1}{V^2\left(x\right)}\frac{{\&PartialD;}^{2}p(x,t)}{{\&PartialD;t}^2}=\frac{\&PartialD;p(x,t)}{\&PartialD;x^2}---(5-1)$

式中：u(x，t)表示位移；V(x)表示速度。通过引入中间变量A(x，t)此方程也可写为两个一阶偏微分方程组的形式：

$\left\{\begin{array}{c}\frac{\&PartialD;p(x,t)}{\&PartialD;t}=V^2\left(x\right)\frac{\&PartialD;A(x,t)}{\&PartialD;x}\\ \frac{\&PartialD;A(x,t)}{\&PartialD;t}=\frac{\&PartialD;p(x,t)}{\&PartialD;x}\end{array}\right.---(5-2)$

容易证明，式(5-1)和式(5-2)是等价的。由式(5-2)构建新的方程组得：

$\left\{\begin{array}{c}\frac{\&PartialD;p^{*}(x,t)}{\&PartialD;t}+d\left(x\right)p^{*}(x,t)=V^2\left(x\right)\frac{\&PartialD;A^{*}(x,t)}{\&PartialD;x}\\ \frac{\&PartialD;A^{*}(x,t)}{\&PartialD;t}+d\left(x\right)A^{*}(x,t)=\frac{\&PartialD;p^{*}(x,t)}{\&PartialD;x}\end{array}\right.---(5-3)$

式中：p^*(x，t)表示新方程的解；d(x)为一个不随时间变化的边界衰减函数。

下面研究式(5-2)的解p(x，t)和式(5-3)的解p^*(x，t)的关系。首先对式(5-3)做关于时间的Fourier变换得：

$\left\{\begin{array}{c}\mathrm{i\&omega;}{\hat{p}}^{*}(x,\mathrm{\&omega;})+d\left(x\right){\hat{p}}^{*}(x,\mathrm{\&omega;})=V^2\left(x\right)\frac{\&PartialD;{\hat{A}}^{*}(x,\mathrm{\&omega;})}{\&PartialD;x}\\ \mathrm{i\&omega;}{\hat{A}}^{*}(x,\mathrm{\&omega;})+d\left(x\right){\hat{A}}^{*}(x,\mathrm{\&omega;})=\frac{\&PartialD;{\hat{u}}^{*}(x,\mathrm{\&omega;})}{\&PartialD;x}\end{array}\right.---(5-4)$

其中，

分别代表p^*(x，ω)、A^*(x，ω)的Fourier变换。

对x做坐标变换：

$\tilde{x}\left(x\right)=x-\frac{i}{\mathrm{\&omega;}}{\&Integral;}_0^{x}d\left(s\right)\mathrm{ds}---(5-5)$

则有：

$\frac{\&PartialD;\tilde{x}\left(x\right)}{\&PartialD;x}=1-\frac{\mathrm{id}\left(x\right)}{\mathrm{\&omega;}}---(5-6)$

$\frac{\&PartialD;{\hat{A}}^{*}\left(x\right)}{\&PartialD;x}=\frac{\&PartialD;{\hat{A}}^{*}\left(\tilde{x}\right)}{\&PartialD;\tilde{x}}\frac{\&PartialD;\tilde{x}\left(x\right)}{\&PartialD;x}$

$=\frac{\mathrm{i\&omega;}+d\left(x\right)}{\mathrm{i\&omega;}}\frac{\&PartialD;{\hat{A}}^{*}\left(\tilde{x}\right)}{\&PartialD;\tilde{x}}---(5-7)$

$\frac{\&PartialD;{\hat{p}}^{*}\left(x\right)}{\&PartialD;x}=\frac{\&PartialD;{\hat{p}}^{*}\left(\tilde{x}\right)}{\&PartialD;\tilde{x}}\frac{\&PartialD;\tilde{x}\left(x\right)}{\&PartialD;\tilde{x}}$

$=\frac{\mathrm{i\&omega;}+d\left(x\right)}{\mathrm{i\&omega;}}\frac{\&PartialD;{\hat{p}}^{*}\left(\tilde{x}\right)}{\&PartialD;\tilde{x}}---(5-8)$

将式(5-7)、式(5-8)分别代入式(5-4)并整理得：

$\left\{\begin{array}{c}\mathrm{i\&omega;}{\hat{p}}^{*}(\tilde{x},\mathrm{\&omega;})=V^2\left(\tilde{x}\right)\frac{\&PartialD;{\hat{A}}^{*}(\tilde{x},\mathrm{\&omega;})}{\&PartialD;\tilde{x}}\\ \mathrm{i\&omega;}{\hat{A}}^{*}(\tilde{x},\mathrm{\&omega;})=\frac{\&PartialD;{\hat{p}}^{*}(\tilde{x},t)}{\&PartialD;\tilde{x}}\end{array}\right.---(5-9)$

对式(5-9)做关于时间的Fourier逆变换，可得：

$\left\{\begin{array}{c}\frac{\&PartialD;p^{*}(\tilde{x},t)}{\&PartialD;t}=V^2\left(\tilde{x}\right)\frac{\&PartialD;A^{*}(\tilde{x},t)}{\&PartialD;\tilde{x}}\\ \frac{\&PartialD;A^{*}(x,t)}{\&PartialD;t}=\frac{\&PartialD;p^{*}(\tilde{x},t)}{\&PartialD;\tilde{x}}\end{array}\right.---(5-10)$

现在来比较式(5-2)和式(5-10)，两对式子(方程组)在形式上完全相同，所以他们有相同形式的解，但是在不同的空间坐标中，假设式(5-2)的解是p(x，t)，则式(5-10)的解就是

即式(5-3)的解为：

$p(x-\frac{i}{\mathrm{\&omega;}}{\&Integral;}_0^{x}d\left(s\right)\mathrm{ds},t)$

也就是说，可以先求解式(5-2)，然后再对它们的解做坐标变换，就可得到式(5-3)的解。现在考察式(5-3)解的性质。式(5-2)在均匀介质中有如下形式的特解：

p(x，t)＝p₀exp[-i(k_xx-ωt)]     (5-11)

由上述讨论可知，对应的式(5-3)的解为：

$p^{*}(x,t)=p_0\exp [-i(k_x\tilde{x}-\mathrm{\&omega;t})]---(5-12)$

将式(5-5)代入上式，并简化可得：

$p^{*}(x,t)=p_0\exp [-i(k_{x}x-\mathrm{\&omega;t})]\&times;\exp [-\frac{k_x}{\mathrm{\&omega;}}{\&Integral;}_0^{x}d\left(s\right)\mathrm{ds}]---(5-13)$

上述解的振幅比为：

$\frac{\left|\right|p^{*}(x,t)\left|\right|}{\left|\right|p(x,t)\left|\right|}=\exp [-\frac{k_x}{\mathrm{\&omega;}}{\&Integral;}_0^{x}d\left(s\right)\mathrm{ds}]---(5-14)$

此式说明，式(5-3)的解相对式(5-2)的解是衰减的，而d(x)起到了一个衰减系数的作用，波场按传播距离的呈指数衰减，衰减速度很快。因而不会产生边界反射。

下面用一个简单的示例说明一维波动方程完全匹配的效果。

图5是一维波动方程完全匹配边界条件应用效果对比图。

其中，(a)未加衰减量的一维波动方程正演；(b)加入衰减量后的一维波动方程正演；(c)进一步加大衰落减量的一维波动方程正演。

从图5可以看出加入吸收系数后振幅开始衰减，衰减系数越大，振幅衰减越快。示例表明完全匹配层边界条件是可以用于波场边界上进行振幅衰减的，而衰减函数d(x)选择可以灵活调整完全匹配层边界条件。

②二维波动方程完全匹配层原理

对下列的二维波动方程：

$\frac{1}{V^2(x,z)}\frac{{\&PartialD;}^{2}p(x,z,t)}{{\&PartialD;t}^2}=\frac{{\&PartialD;}^{2}p(x,z,t)}{{\&PartialD;x}^2}+\frac{{\&PartialD;}^{2}p(x,z,t)}{{\&PartialD;z}^2}---(5-15)$

其中：p(x，z，t)是位移函数；V(x，z)是介质的速度，此方程可写为等价的一阶偏微分方程组的形式：

$\left\{\begin{array}{c}p=u_1+u_2\\ \frac{{\&PartialD;u}_1}{\&PartialD;t}=V^2(x,z)\frac{\&PartialD;A_1}{\&PartialD;x}\\ \frac{{\&PartialD;u}_2}{\&PartialD;t}=V^2(x,z)\frac{{\&PartialD;A}_2}{\&PartialD;z}\\ \frac{{\&PartialD;A}_1}{\&PartialD;t}=\frac{{\&PartialD;u}_1}{\&PartialD;x}+\frac{{\&PartialD;u}_2}{\&PartialD;x}\\ \frac{{\&PartialD;A}_2}{\&PartialD;t}=\frac{{\&PartialD;u}_1}{\&PartialD;z}+\frac{{\&PartialD;u}_2}{\&PartialD;z}\end{array}\right.---(5-16)$

其中：u₁、u₂、A₁、A₂是引入的中间变量。

相应的完全匹配层控制方程为：

$\left\{\begin{array}{c}{p}^{*}={u^{*}}_1+{u^{*}}_2\\ \frac{\&PartialD;{u^{*}}_1}{\&PartialD;t}+d_1\left(x\right){u^{*}}_1=V^2(x,z)\frac{\&PartialD;{A^{*}}_1}{\&PartialD;x}\\ \frac{\&PartialD;{u^{*}}_2}{\&PartialD;t}+d_2\left(z\right){u^{*}}_2=V^2(x,z)\frac{\&PartialD;{A^{*}}_2}{\&PartialD;z}\\ \frac{\&PartialD;{A^{*}}_1}{\&PartialD;t}+d_1\left(x\right){A^{*}}_1=\frac{\&PartialD;{u^{*}}_1}{\&PartialD;x}+\frac{\&PartialD;{u^{*}}_2}{\&PartialD;x}\\ \frac{\&PartialD;{A^{*}}_2}{\&PartialD;t}+d_2\left(z\right){A^{*}}_2=\frac{\&PartialD;{u^{*}}_1}{\&PartialD;z}+\frac{\&PartialD;{u^{*}}_2}{\&PartialD;z}\end{array}\right.---(5-17)$

和一维情况相同，式(5-17)的解也是衰减的。d₁(x)和d₂(z)分别为x方向和z方向的衰减系数，也就是说d₁(x)起到衰减x方向传播的波、d₂(z)起到衰减z方向传播的波的作用。对于任意方向传播的波，可以通过矢量分解，分解成x方向和z方向传播的波，分别进行衰减。波场是按传播距离的指数规律衰减，衰减速度很快。当衰减系数d₁(x)、d₂(z)随空间位置变化时，不会在介质中产生任何反射。这些性质使得这种介质特别适合用作波动方程的边界吸收介质。

图6是完全匹配层吸收边界示意图。利用完全匹配层作为吸收边界的基本做法是在所研究区域的四周引入完全匹配层。如图6所示，区域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方向和边界的速度相等。这样在计算边界的周围都有完全匹配层吸收介质，波由区域内通过边界传播到完全匹配层时，不会产生任何反射。波在完全匹配层中传播时，不会产生反射，并且按传播距离的指数规律衰减。当波传播到完全匹配层的边界时，波场近似为零，也不会产生反射。

同理，可以推出到三维波动方程完全匹配层的边界，其形式与二维类同。

③二维波动方程任意差分精细积分法求解完全匹配层计算公式

用任意差分精细积分进行完全匹配边界计算，需要对(5-17)式方程组在计算区域内和边界用同样方法求解，若仅在边界上的用完全匹配层解(5-17)式，就将出现波场的不相容性，从而出现数值计算引起的波场反射。为了消除这种反射，必须推导方程组(5-17)式的任意差分精细积分完全匹配边界计算公式。

由(5-17)式，令v＝u^* ₁，w＝u^* ₂，p_x＝A^* ₁，p_z＝A^* ₂整理后得到：

$\frac{\&PartialD;}{\&PartialD;t}\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]+d\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& V& V\\ 0& V& 0& 0\\ 0& 0& 0& 0\end{array}\right]\frac{\&PartialD;}{\&PartialD;x}\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]+\left[\begin{array}{cccc}0& 0& V& V\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ V& 0& 0& 0\end{array}\right]\frac{\&PartialD;}{\&PartialD;z}\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]---(5-18)$

上式的右端项 $\frac{\&PartialD;}{\&PartialD;x}\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]$ 可在点的x邻域作泰勒展开：

$\left[\begin{array}{c}{v}_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{zi}}\end{array}\right]=\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right]+\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{x}_i\frac{\&PartialD;}{\&PartialD;x}\right)}^k{\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]}_{x=x_j}+O\left[{\left(\mathrm{\&Delta;}{x}_i\right)}^{m+1}\right]---(5-19)$

移项加权得

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(\left[\begin{array}{c}{v}_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{zi}}\end{array}\right]-\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right])=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{x}_k\frac{\&PartialD;}{\&PartialD;x}\right)}^k{\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]}_{x=x_j}+O\left[{\left(\mathrm{\&Delta;}{x}_i\right)}^{m+1}\right])---(5-20)$

令(5-20)右端展开并等于(5-18)关于x的偏导数得：

$\left\{\begin{array}{cc}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\mathrm{\&Delta;}{x}_i=1& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^2=0\\ \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^3=0& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^4=0\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^{m-1}=0& \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\left(\mathrm{\&Delta;}{x}_i\right)}^m=0\end{array}\right.$

解该方程组即可求出加权系数α_i即：

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(\left[\begin{array}{c}{v}_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{zi}}\end{array}\right]-\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right])=\frac{\&PartialD;}{\&PartialD;x}\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]---(5-21)$

同理 $\frac{\&PartialD;}{\&PartialD;z}\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]$ 在j点的z点邻域作泰勒展开：

$\left[\begin{array}{c}{v}_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{zl}}\end{array}\right]=\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right]+\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{z}_l\frac{\&PartialD;}{\&PartialD;x}\right)}^k{\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]}_{z=z_j}+O\left[{\left(\mathrm{\&Delta;}{z}_l\right)}^{m+1}\right]---(5-22)$

移项加权得：

$\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l(\left[\begin{array}{c}{v}_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{zl}}\end{array}\right]-\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right])=\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l(\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{k!}{\left(\mathrm{\&Delta;}{z}_l\frac{\&PartialD;}{\&PartialD;z}\right)}^k{\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]}_{x=x_j}+O\left[{\left({\mathrm{\&Delta;z}}_l\right)}^{m+1}\right])---(5-23)$

令上式右端展开并等于关于z的偏导数得：

$\left\{\begin{array}{cc}\underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\mathrm{\&Delta;}{z}_l=1& \underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l{\left(\mathrm{\&Delta;}{z}_l\right)}^2=0\\ \underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l{\left(\mathrm{\&Delta;}{z}_l\right)}^3=0& \underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l{\left(\mathrm{\&Delta;}{z}_l\right)}^4=0\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l{\left(\mathrm{\&Delta;}{z}_l\right)}^{m-1}=0& \underset{l=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l{\left(\mathrm{\&Delta;}{z}_l\right)}^m=0\end{array}\right.$

解该方程组即可求出加权系数α_l即：

$\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l(\left[\begin{array}{c}{v}_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{zl}}\end{array}\right]-\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right])=\frac{\&PartialD;}{\&PartialD;z}\left[\begin{array}{c}v\\ w\\ p_x\\ p_z\end{array}\right]---(5-24)$

整理得：

$\frac{\&PartialD;}{\&PartialD;t}\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right]+d\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right]=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& V& V\\ 0& V& 0& 0\\ 0& 0& 0& 0\end{array}\right]\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(\left[\begin{array}{c}{v}_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{zi}}\end{array}\right]-\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right])$

$+\left[\begin{array}{cccc}0& 0& V& V\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ V& 0& 0& 0\end{array}\right]\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l(\left[\begin{array}{c}{v}_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{zl}}\end{array}\right]-\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right])---(5-25)$

因为： $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i=0,$ 所以有：

$\frac{\&PartialD;}{\&PartialD;t}\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right]+d\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right]=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& V& V\\ 0& V& 0& 0\\ 0& 0& 0& 0\end{array}\right]\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{zi}}\end{array}\right]$

$+\left[\begin{array}{cccc}0& 0& V& V\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ V& 0& 0& 0\end{array}\right]\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{v}_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{zl}}\end{array}\right]---(5-26)$

设： $\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& V& V\\ 0& V& 0& 0\\ 0& 0& 0& 0\end{array}\right]=A;$ $\left[\begin{array}{cccc}0& 0& V& V\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ V& 0& 0& 0\end{array}\right]=B;$

$P=\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right];$ $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{i,n}\\ w_{i,n}\\ p_{\mathrm{xi},n}\\ p_{\mathrm{zi},n}\end{array}\right]=b_{i,n};$ $\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{l,n}\left[\begin{array}{c}{v}_{l,n}\\ w_{l,n}\\ p_{\mathrm{xl},n}\\ p_{\mathrm{zl},n}\end{array}\right]=b_{l,n}$

d(s)＝0.5×pml_max×[1-cos(2×π×s)/pml_thick)]     (5-27)

参数s的选择根据波的传播方向而定，若是向x方向传播，则s为x。若为z方向则s选为z。pml_max为衰减系数，其大小的选定可根据边界大小而定，但不可过大，会造成计算的不稳定，pml_thick为完全匹配层的宽度。

即 $\frac{\&PartialD;P}{\&PartialD;t}+d\left(s\right)P=A{b}_{i,n}+B{b}_{l,n}---(5-28)$

解此方程得：

$P=C{e}^{-\&Integral;d\left(s\right)\mathrm{dt}}+e^{-\&Integral;d\left(s\right)\mathrm{dt}}\&Integral;(A{b}_{i,n}+B{b}_{l,n})e^{\&Integral;d\left(s\right)\mathrm{dt}}\mathrm{dt}---(5-29)$

整理得

$P_j=C{e}^{-d\left(s\right)(t-t_n)}+\frac{1}{d\left(s\right)}{e}^{-d\left(s\right)(t-t_n)}(A{b}_{i,n}+B{b}_{l,n})e^{d\left(s\right)(t-t_n)}---(5-30)$

选择积分子域[t_n，t_n+1]求出t_n+1时刻的值

进一步提高精度可将 $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{zi}}\end{array}\right]$ 展开到t-t_n的二阶导数项：

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{zi}}\end{array}\right]=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{zi}}\end{array}\right]+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v\&prime;}_i\\ {w\&prime;}_i\\ {p\&prime;}_{\mathrm{xi}}\\ {p\&prime;}_{\mathrm{zi}}\end{array}\right](t-t_n)+\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_i^{\&prime;\&prime;}\\ w_i^{\&prime;\&prime;}\\ p_{\mathrm{xi}}^{\&prime;\&prime;}\\ p_{\mathrm{zi}}^{\&prime;\&prime;}\end{array}\right]{(t-t_n)}^2---(5-32)$ 同理将 $\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{v}_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{zl}}\end{array}\right]$ 展开为t-t_n的二阶导数项得：

$\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{v}_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{zl}}\end{array}\right]=\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{l,n}\left[\begin{array}{c}{v}_{l,n}\\ w_{l,n}\\ p_{\mathrm{xl},n}\\ p_{\mathrm{zl},n}\end{array}\right]+\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{v}_{l,n}^{\&prime;}\\ w_{l,n}^{\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;}\end{array}\right](t-t_n)+\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{l,n}^{\&prime;\&prime;}\\ w_{l,n}^{\&prime;\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;\&prime;}\end{array}\right]{(t-t_n)}^2---(5-33)$

代入并整理(5-18)得：

$\frac{\&PartialD;}{\&PartialD;t}\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right]+d\left[\begin{array}{c}{v}_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{zj}}\end{array}\right]=A\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{i,n}\\ w_{i,n}\\ p_{\mathrm{xi},n}\\ p_{\mathrm{zi},n}\end{array}\right]+A\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{i,n}^{\&prime;}\\ w_{i,n}^{\&prime;}\\ p_{\mathrm{xi},n}^{\&prime;}\\ p_{\mathrm{zi},n}^{\&prime;}\end{array}\right](t-t_n)+A\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{i,n}^{\&prime;\&prime;}\\ w_{i,n}^{\&prime;\&prime;}\\ p_{\mathrm{xi},n}^{\&prime;\&prime;}\\ p_{\mathrm{zi},n}^{\&prime;\&prime;}\end{array}\right]{(t-t_n)}^2$

$+B\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{l,n}\left[\begin{array}{c}{v}_{l,n}\\ w_{l,n}\\ p_{\mathrm{xl},n}\\ p_{\mathrm{zl},n}\end{array}\right]+B\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{v}_{l,n}^{\&prime;}\\ w_{l,n}^{\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;}\end{array}\right](t-t_n)+B\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{l,n}^{\&prime;\&prime;}\\ w_{l,n}^{\&prime;\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;\&prime;}\end{array}\right]{(t-t_n)}^2---(5-34)$

令： $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{i,n}\\ w_{i,n}\\ p_{\mathrm{xi},n}\\ p_{\mathrm{zi},n}\end{array}\right]=b_{i,n};$ $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{i,n}^{\&prime;}\\ w_{i,n}^{\&prime;}\\ p_{\mathrm{xi},n}^{\&prime;}\\ p_{\mathrm{zi},n}^{\&prime;}\end{array}\right]=e_{i,n};$ $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{i,n}^{\&prime;\&prime;}\\ w_{i,n}^{\&prime;\&prime;}\\ p_{\mathrm{xi},n}^{\&prime;\&prime;}\\ p_{\mathrm{zi},n}^{\&prime;\&prime;}\end{array}\right]=d_{i,n};$

$\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{l,n}\left[\begin{array}{c}{v}_{l,n}\\ w_{l,n}\\ p_{\mathrm{xl},n}\\ p_{\mathrm{zl},n}\end{array}\right]=b_{l,n};$ $\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{v}_{l,n}^{\&prime;}\\ w_{l,n}^{\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;}\end{array}\right]=e_{l,n};$ $\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{v}_{l,n}^{\&prime;\&prime;}\\ w_{l,n}^{\&prime;\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;\&prime;}\end{array}\right]=d_{l,n}$

即： $\frac{\&PartialD;P}{\&PartialD;t}+d\left(s\right)P=A{b}_{i,n}+A{e}_{i,n}(t-t_n)+\frac{1}{2}A{d}_{i,n}{(t-t_n)}^2$

$+B{b}_{l,n}+B{e}_{l,n}(t-t_n)+\frac{1}{2}B{d}_{l,n}{(t-t_n)}^2---(5-35)$

解此方程得：

$P=C{e}^{-d\left(s\right)(t-t_n)}+e^{-d\left(s\right)(t-t_n)}[\frac{1}{d\left(s\right)}A{b}_{i,n}{e}^{d\left(s\right)(t-t_n)}$

$+A\frac{1}{d\left(s\right)}{e}_{i,n}(t-t_n)e^{d\left(s\right)(t-t_n)}-A\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}{e}^{d\left(s\right)(t-t_n)}$

$+\frac{1}{2}A\frac{1}{d\left(s\right)}{d}_{i,n}{(t-t_n)}^2{e}^{d\left(s\right)(t-t_n)}-A\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}(t-t_n)e^{d\left(s\right)(t-t_n)}$

$+A\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}{e}^{d\left(s\right)(t-t_n)}+\frac{1}{d\left(s\right)}B{b}_{i,n}{e}^{d\left(s\right)(t-t_n)}$

$+B\frac{1}{d\left(s\right)}{e}_{i,n}(t-t_n)e^{d\left(s\right)(t-t_n)}-B\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}{e}^{d\left(s\right)(t-t_n)}$

$+\frac{1}{2}B\frac{1}{d\left(s\right)}{d}_{i,n}{(t-t_n)}^2{e}^{d\left(s\right)(t-t_n)}-B\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}(t-t_n)e^{d\left(s\right)(t-t_n)}$

$+B\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}{e}^{d\left(s\right)(t-t_n)}]---(5-36)$

代换C得任意差分精细积分进行完全匹配边界递推求解计算公式：

$P_{j,n+1}=(P_{j,n}-A\frac{1}{d\left(s\right)}{b}_{i,n}-A\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}-A\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}$

$-B\frac{1}{d\left(s\right)}{b}_{i,n}-B\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}-B\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n})e^{-d\left(s\right)\mathrm{\&Delta;t}}$

$+[\frac{1}{d\left(s\right)}A{b}_{i,n}+A\frac{1}{d\left(s\right)}{e}_{i,n}\mathrm{Vt}-A\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}$

$+\frac{1}{2}A\frac{1}{d\left(s\right)}{d}_{i,n}{\left(\mathrm{Vt}\right)}^2-A\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}\left(\mathrm{Vt}\right)+A\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}$

$+[\frac{1}{d\left(s\right)}B{b}_{i,n}+B\frac{1}{d\left(s\right)}{e}_{i,n}\mathrm{Vt}-B\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}$

$+\frac{1}{2}B\frac{1}{d\left(s\right)}{d}_{i,n}{\left(\mathrm{Vt}\right)}^2-B\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}\left(\mathrm{Vt}\right)+B\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}]---(5-37)$

④三维波动方程任意差分精细积分法求解完全匹配层计算公式

根据以上二维方程推导同理可得三维方程完全匹配边界条件的求解公式。

三维波动方程等价的一阶偏微分方程组可由下式表示：

$\frac{\&PartialD;}{\&PartialD;t}\left[\begin{array}{c}u\\ v\\ w\\ p_x\\ p_y\\ p_z\end{array}\right]+d\left[\begin{array}{c}u\\ v\\ w\\ p_x\\ p_y\\ p_z\end{array}\right]=\left[\begin{array}{cccccc}0& 0& 0& V& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\frac{\&PartialD;}{\&PartialD;x}\left[\begin{array}{c}u\\ v\\ w\\ p_x\\ p_y\\ p_z\end{array}\right]$

$+\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& V& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\frac{\&PartialD;}{\&PartialD;y}\left[\begin{array}{c}u\\ v\\ w\\ p_x\\ p_y\\ p_z\end{array}\right]+\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& V\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\end{array}\right]\frac{\&PartialD;}{\&PartialD;z}\left[\begin{array}{c}u\\ v\\ w\\ p_x\\ p_y\\ p_z\end{array}\right]---(5-38)$

令： $A=\left[\begin{array}{cccccc}0& 0& 0& V& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right];$ $B=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& V& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right];$

$D=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& V\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\end{array}\right];$ $P=\left[\begin{array}{c}u\\ v\\ w\\ p_x\\ p_y\\ p_z\end{array}\right];$

作泰勒展开并整理得：

$\frac{\&PartialD;}{\&PartialD;t}\left[\begin{array}{c}{u}_j\\ v_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{yj}}\\ p_{\mathrm{zj}}\end{array}\right]+d\left[\begin{array}{c}{u}_j\\ v_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{yj}}\\ p_{\mathrm{zj}}\end{array}\right]=\left[\begin{array}{cccccc}0& 0& 0& V& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i(\left[\begin{array}{c}{u}_i\\ v_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{yi}}\\ p_{\mathrm{zi}}\end{array}\right]-\left[\begin{array}{c}{u}_j\\ v_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{yj}}\\ p_{\mathrm{zj}}\end{array}\right])$

$+\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& V& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l(\left[\begin{array}{c}{u}_l\\ v_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{yl}}\\ p_{\mathrm{zl}}\end{array}\right]-\left[\begin{array}{c}{u}_j\\ v_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{yj}}\\ p_{\mathrm{zj}}\end{array}\right])$

$+\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& V\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ V& V& V& 0& 0& 0\end{array}\right]\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_m(\left[\begin{array}{c}{u}_m\\ v_m\\ w_m\\ p_{\mathrm{xm}}\\ p_{\mathrm{ym}}\\ p_{\mathrm{zm}}\end{array}\right]-\left[\begin{array}{c}{u}_j\\ v_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{yj}}\\ p_{\mathrm{zj}}\end{array}\right])---(5-39)$

令 $b_{i,n}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{u}_i\\ v_i\\ w_i\\ p_{\mathrm{xi}}\\ p_{\mathrm{yi}}\\ p_{\mathrm{zi}}\end{array}\right];$ $b_{l,n}=\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{u}_l\\ v_l\\ w_l\\ p_{\mathrm{xl}}\\ p_{\mathrm{yl}}\\ p_{\mathrm{zl}}\end{array}\right];$ $b_{m,n}=\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_m\left[\begin{array}{c}{u}_m\\ v_m\\ w_m\\ p_{\mathrm{xm}}\\ p_{\mathrm{ym}}\\ p_{\mathrm{zm}}\end{array}\right];$

解得

$P_j=C{e}^{-d\left(s\right)(t-t_n)}+\frac{1}{d\left(s\right)}{e}^{-d\left(s\right)(t-t_n)}(A{b}_{i,n}+B{b}_{l,n}+D{b}_{m,n})e^{d\left(s\right)(t-t_n)}---(5-40)$

选择积分子域[t_n，t_n+1]确定C得

进一步提高精度后得：

$\frac{\&PartialD;}{\&PartialD;t}\left[\begin{array}{c}{u}_j\\ v_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{yj}}\\ p_{\mathrm{zj}}\end{array}\right]+d\left[\begin{array}{c}{u}_j\\ v_j\\ w_j\\ p_{\mathrm{xj}}\\ p_{\mathrm{yj}}\\ p_{\mathrm{zj}}\end{array}\right]=A\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{u}_{i,n}\\ v_{i,n}\\ w_{i,n}\\ p_{\mathrm{xi},n}\\ p_{\mathrm{yi},n}\\ p_{\mathrm{zi},n}\end{array}\right]+A\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{u}_{i,n}^{\&prime;}\\ v_{i,n}^{\&prime;}\\ w_{i,n}^{\&prime;}\\ p_{\mathrm{xi},n}^{\&prime;}\\ p_{\mathrm{yi},n}^{\&prime;}\\ p_{\mathrm{zi},n}^{\&prime;}\end{array}\right](t-t_n)+\frac{1}{2}A\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{u}_{i,n}^{\&prime;\&prime;}\\ v_{i,n}^{\&prime;\&prime;}\\ w_{i,n}^{\&prime;\&prime;}\\ p_{\mathrm{xi},n}^{\&prime;\&prime;}\\ p_{\mathrm{yi},n}^{\&prime;\&prime;}\\ p_{\mathrm{zi},n}^{\&prime;\&prime;}\end{array}\right]{(t-t_n)}^2$

$+B\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{u}_{l,n}\\ v_{l,n}\\ w_{l,n}\\ p_{\mathrm{xl},n}\\ p_{\mathrm{yl},n}\\ p_{\mathrm{zl},n}\end{array}\right]+B\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{u}_{l,n}^{\&prime;}\\ v_{l,n}^{\&prime;}\\ w_{l,n}^{\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;}\\ p_{\mathrm{yl},n}^{\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;}\end{array}\right](t-t_n)+\frac{1}{2}B\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{u}_{l,n}^{\&prime;\&prime;}\\ v_{l,n}^{\&prime;\&prime;}\\ w_{l,n}^{\&prime;\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;\&prime;}\\ p_{\mathrm{yl},n}^{\&prime;\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;\&prime;}\end{array}\right]{(t-t_n)}^2$

$+D\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_m\left[\begin{array}{c}{u}_{m,n}\\ v_{m,n}\\ w_{m,n}\\ p_{\mathrm{xm},n}\\ p_{\mathrm{ym},n}\\ p_{\mathrm{zm},n}\end{array}\right]+D\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_m\left[\begin{array}{c}{u}_{m,n}^{\&prime;}\\ v_{m,n}^{\&prime;}\\ w_{m,n}^{\&prime;}\\ p_{\mathrm{xm},n}^{\&prime;}\\ p_{\mathrm{ym},n}^{\&prime;}\\ p_{\mathrm{zm},n}^{\&prime;}\end{array}\right](t-t_n)+\frac{1}{2}D\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_m\left[\begin{array}{c}{u}_{m,n}^{\&prime;\&prime;}\\ v_{m,n}^{\&prime;\&prime;}\\ w_{m,n}^{\&prime;\&prime;}\\ p_{\mathrm{xm},n}^{\&prime;\&prime;}\\ p_{\mathrm{ym},n}^{\&prime;\&prime;}\\ p_{\mathrm{zm},n}^{\&prime;\&prime;}\end{array}\right]{(t-t_n)}^2$

令： $b_{i,n}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{u}_{i,n}\\ v_{i,n}\\ w_{i,n}\\ p_{\mathrm{xi},n}\\ p_{\mathrm{yi},n}\\ p_{\mathrm{zi},n}\end{array}\right];$ $e_{i,n}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{u}_{i,n}^{\&prime;}\\ v_{i,n}^{\&prime;}\\ w_{i,n}^{\&prime;}\\ p_{\mathrm{xi},n}^{\&prime;}\\ p_{\mathrm{yi},n}^{\&prime;}\\ p_{\mathrm{zi},n}^{\&prime;}\end{array}\right];$ $d_{i,n}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\left[\begin{array}{c}{u}_{i,n}^{\&prime;\&prime;}\\ v_{i,n}^{\&prime;\&prime;}\\ w_{i,n}^{\&prime;\&prime;}\\ p_{\mathrm{xi},n}^{\&prime;\&prime;}\\ p_{\mathrm{yi},n}^{\&prime;\&prime;}\\ p_{\mathrm{zi},n}^{\&prime;\&prime;}\end{array}\right];$

$b_{l,n}=\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{u}_{l,n}\\ v_{l,n}\\ w_{l,n}\\ p_{\mathrm{xl},n}\\ p_{\mathrm{yl},n}\\ p_{\mathrm{zl},n}\end{array}\right];$ $e_{l,n}=\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{u}_{l,n}^{\&prime;}\\ v_{l,n}^{\&prime;}\\ w_{l,n}^{\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;}\\ p_{\mathrm{yl},n}^{\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;}\end{array}\right];$ $d_{l,n}=\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_l\left[\begin{array}{c}{u}_{l,n}^{\&prime;\&prime;}\\ v_{l,n}^{\&prime;\&prime;}\\ w_{l,n}^{\&prime;}\\ p_{\mathrm{xl},n}^{\&prime;\&prime;}\\ p_{\mathrm{yl},n}^{\&prime;\&prime;}\\ p_{\mathrm{zl},n}^{\&prime;\&prime;}\end{array}\right];$

$b_{m,n}=\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_m\left[\begin{array}{c}{u}_{m,n}\\ v_{m,n}\\ w_{m,n}\\ p_{\mathrm{xm},n}\\ p_{\mathrm{ym},n}\\ p_{\mathrm{zm},n}\end{array}\right];$ $e_{m,n}=\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_m\left[\begin{array}{c}{u}_{m,n}^{\&prime;}\\ v_{m,n}^{\&prime;}\\ w_{m,n}^{\&prime;}\\ p_{\mathrm{xm},n}^{\&prime;}\\ p_{\mathrm{ym},n}^{\&prime;}\\ p_{\mathrm{zm},n}^{\&prime;}\end{array}\right];$ $d_{m,n}=\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_m\left[\begin{array}{c}{u}_{m,n}^{\&prime;\&prime;}\\ v_{m,n}^{\&prime;\&prime;}\\ w_{m,n}^{\&prime;\&prime;}\\ p_{\mathrm{xm},n}^{\&prime;\&prime;}\\ p_{\mathrm{ym},n}^{\&prime;\&prime;}\\ p_{\mathrm{zm},n}^{\&prime;\&prime;}\end{array}\right]$

即： $\frac{\&PartialD;P_j}{\&PartialD;t}+d(s{)P}_j=A{b}_{i,n}+A{e}_{i,n}(t-t_n)+\frac{1}{2}A{d}_{i,n}{(t-t_n)}^2$

$+B{b}_{l,n}+B{e}_{l,n}(t-t_n)+\frac{1}{2}B{d}_{l,n}{(t-t_n)}^2$

$+D{b}_{m,n}+D{e}_{m,n}(t-t_n)+\frac{1}{2}D{d}_{m,n}{(t-t_n)}^2---(5-42)$

$P_j=C{e}^{-d\left(s\right)(t-t_n)}+[\frac{1}{d\left(s\right)}A{b}_{l,n}+A\frac{1}{d\left(s\right)}{e}_{i,n}(t-t_n)-A\frac{1}{{\left[d\left(s\right)\right]}^2}$

$+\frac{1}{2}A\frac{1}{d\left(s\right)}{d}_{i,n}{(t-t_n)}^2-A\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}(t-t_n)+A\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}$

$+\frac{1}{d\left(s\right)}B{b}_{i,n}+B\frac{1}{d\left(s\right)}{e}_{i,n}(t-t_n)-B\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}$

$+\frac{1}{2}B\frac{1}{d\left(s\right)}{d}_{i,n}{(t-t_n)}^2-B\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}(t-t_n)+B\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}$

$+\frac{1}{d\left(s\right)}D{b}_{i,n}+D\frac{1}{d\left(s\right)}{e}_{i,n}(t-t_n)-D\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}$

$+\frac{1}{2}D\frac{1}{d\left(s\right)}{d}_{i,n}{(t-t_n)}^2-D\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}(t-t_n)+D\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}]---(5-43)$

$P_{j,n+1}=(P_{j,n}-A\frac{1}{d\left(s\right)}{b}_{i,n}-A\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}-A\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}$

$-B\frac{1}{d\left(s\right)}{b}_{i,n}-B\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}-B\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}-D\frac{1}{d\left(s\right)}{b}_{i,n}$

$-D\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}-D\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n})e^{-d\left(s\right)\mathrm{\&Delta;t}}$

$+[\frac{1}{d\left(s\right)}A{b}_{i,n}+A\frac{1}{d\left(s\right)}{e}_{i,n}\mathrm{Vt}-A\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}$

$+\frac{1}{2}A\frac{1}{d\left(s\right)}{d}_{i,n}{\left(\mathrm{Vt}\right)}^2-A\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}\left(\mathrm{Vt}\right)+A\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}$

$+[\frac{1}{d\left(s\right)}B{b}_{i,n}+B\frac{1}{d\left(s\right)}{e}_{i,n}\mathrm{Vt}-B\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}$

$+\frac{1}{2}B\frac{1}{d\left(s\right)}{d}_{i,n}{\left(\mathrm{Vt}\right)}^2-B\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}\left(\mathrm{Vt}\right)+B\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}]$

$+[\frac{1}{d\left(s\right)}D{b}_{i,n}+D\frac{1}{d\left(s\right)}{e}_{i,n}\mathrm{Vt}-D\frac{1}{{\left[d\left(s\right)\right]}^2}{e}_{i,n}$

$+\frac{1}{2}D\frac{1}{d\left(s\right)}{d}_{i,n}{\left(\mathrm{Vt}\right)}^2-D\frac{1}{{\left[d\left(s\right)\right]}^2}{d}_{i,n}\left(\mathrm{Vt}\right)+D\frac{1}{{\left[d\left(s\right)\right]}^3}{d}_{i,n}]---(5-44)$

公式(5-44)是三维波动方程完全匹配边界条件的求解递推公式，在波场内部其解满足方程(2-1)。波场在边界传送过程中由于d(s)的衰减作用使得波场振幅很快衰减从而达到去除边界反射波的目的。

本发明的任意差分精细积分方法有较高的精度和较好的数值稳定性，可以在几乎不增加计算量的情况下较大地提高了计算精度和稳定性。完全匹配边界条件可以用于波场边界上进行振幅衰减。用任意差分精细积分法导出的完全匹配层边界条件对边界反射有很好的衰减作用。本发明可以为复杂区地震波传播规律研究提供了一个高精度、稳定性能好的数据体。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种地震波动方程生成方法，其特征在于，所述方法包括以下步骤：

获取声波方程数据；

获取地质参数信息；

根据所述的声波方程数据和所述的地质参数信息进行任意差分精细积分，得出地震传播方程；

求出稳定性条件、初始输入条件、边界处理条件；

由所述的地震传播方程以及所述的稳定性条件、初始输入条件、边界处理条件生成地震波动方程。

2.根据权利要求1所述的地震波动方程生成方法，其特征在于，所述获取地质参数信息的步骤包括：

采集地震数据；

利用promax软件处理所述地震数据得出地震剖面；

对所述地震剖面进行分析得出对应的地质参数信息。

3.根据权利要求1所述的地震波动方程生成方法，其特征在于，所述进行任意差分精细积分的步骤包括：

进行时间有限差分和空间精细积分的交替来进行所述的任意差分精细积分。

4.根据权利要求1所述的地震波动方程生成方法，其特征在于，所述边界条件是指完全匹配层吸收边界条件，以在完全匹配边界条件下进行任意差分精细积分。

5.根据权利要求1所述的地震波动方程生成方法，其特征在于，所述地震波动方程为一半解析解。

6.一种地震波动方程生成系统，其特征在于，该地震波动方程生成系统包括：

声波方程获取单元，用于获取声波方程数据；

地质参数获取单元，用于获取地质参数信息；

传播方程生成单元，用于根据所述的声波方程数据和所述的地质参数信息进行任意差分精细积分，得出地震传播方程；

条件生成单元，用于求出稳定性条件、初始输入条件、边界处理条件；

波动方程生成单元，用于由所述的地震传播方程以及所述的稳定性条件、初始输入条件、边界处理条件生成地震波动方程。

7.根据权利要求6所述的地震波动方程生成系统，其特征在于，所述地质参数获取单元包括：

地震勘探模块，采集地震数据；

数据处理模块，利用promax软件处理所述地震数据得出地震剖面；

分析模块，对所述地震剖面进行分析得出对应的地质参数信息。

8.根据权利要求6所述的地震波动方程生成系统，其特征在于，所述传播方程生成单元包括：

有限差分模块，用于进行时间有限差分；

精细积分模块，用于进行空间精细积分，

所述传播方程生成单元进行时间有限差分和空间精细积分的交替来进行所述的任意差分精细积分。

9.根据权利要求6所述的地震波动方程生成系统，其特征在于，所述边界条件是指完全匹配层吸收边界条件，以在完全匹配边界条件下进行任意差分精细积分。

10.根据权利要求6所述的地震波动方程生成系统，其特征在于，所述地震波动方程为一半解析解。

Images