CN103616721B - 基于二阶偏微分波动方程的pml吸收边界条件的方法 - Google Patents

Abstract

本发明提供了一种基于二阶偏微分波动方程的PML吸收边界条件的方法，包括以下步骤：加载某一采样时刻的震源的纵波波场；计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，并根据所述方向的二阶偏导数对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；重复上一步骤，分别对所述采样时刻的纵波波场在三维空间另外两个方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；重复以上步骤，直至完成所述纵波波场所有采样时刻的处理。采用本发明的完全匹配层吸收边界条件的实现复杂度低、计算存储量小。

Description

基于二阶偏微分波动方程的PML吸收边界条件的方法

技术领域

本发明涉及地球物理勘探领域中一种波动方程数值模拟技术，尤其涉及一种基于二阶偏微分波动方程的PML（perfectmatchedlayer，完全匹配层）吸收边界条件的方法。

背景技术

地震波数值模拟技术在地球物理勘探领域得到广泛的应用，常用来设计地震采集观测系统、研究复杂储层的地震波响应、实现地震波逆时偏移成像和全波形反演速度建模技术等。

在波动方程数值模拟技术中，一项重要的研究内容就是吸收边界条件。在现实世界中，速度介质模型的边界可能为无穷大；采用计算机进行数值模拟实现时，只能对无限大的模型进行截断，这样就会导致地震波传到截断边界时产生不渴望的虚假反射。因此，必须对这种虚假的边界反射进行吸收或消除，一般使用人为地引入吸收边界条件的方法来解决。

目前应用于波动方程的吸收边界条件主要有三类。第一类是海绵吸收边界条件，在一定层厚的介质中对波场逐渐进行衰减，如阻尼吸收边界条件；第二类是基于波动方程的旁轴近似吸收边界条件，本质上是一种单程波，不足之处是不能很好处理入射角度较大的反射问题；第三类就是基于波动方程的PML吸收边界条件，由于不同完全匹配层之间的复波阻抗相同，波在完全匹配层之间传播时不产生任何反射只产生衰减。到目前为止，完全匹配层吸收边界条件是吸收效果最理想的吸收边界条件，在数值模拟精度较高的场合，一般主要采用完全匹配层吸收边界条件。其中，PML是通过在FDTD（Finite-DifferenceTime-Domain，时域有限差分）区域截断边界设置一种特殊的介质层，该层介质的波阻抗与相邻介质的波阻抗完全匹配，入射波将无反射地穿过分界面而进入PML。

但是，现有的完全匹配层吸收边界条件只适用于一阶偏微分波动方程数值模拟，主要采用的为分裂式完全匹配层吸收边界条件。在计算过程中，为了节省计算存储，现有技术常常在完全匹配层区域对波场进行分裂，对于三维空间问题，其需要考虑的区域包括6个面部区域、12个棱边区域和8个角点区域，导致实现起来非常复杂。

发明内容

本发明的目的在于提供一种基于二阶偏微分波动方程的PML吸收边界条件的方法，以减小在波动方程数值模拟技术中采用完全匹配层吸收边界条件的实现复杂度。

为达到上述目的，本发明提供了一种基于二阶偏微分波动方程的PML吸收边界条件的方法，包括以下步骤：

加载某一采样时刻的震源的纵波波场；

计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，并根据所述方向的二阶偏导数对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；

重复上一步骤，分别对所述采样时刻的纵波波场在三维空间另外两个方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；

重复以上步骤，直至完成所述纵波波场所有采样时刻的处理。

本发明的方法，在所述在整个三维数值模拟计算区域对所述纵波波场进行更新之后，重复以上步骤，直至完成所述纵波波场所有采样时刻的处理之前，还包括：

在整个三维数值模拟计算区域对所述纵波波场进行更新。

本发明的方法，在所述加载某一采样时刻的震源的纵波波场之前，还包括：

预先将纵波波场的三维速度模型在三维空间上进行网格离散化，并设置波动方程数值模拟所需参数；

预先将经过网格离散化处理后的三维速度模型的整个三维空间数值模拟区域划分为主区域和完全匹配层区域，并在所述完全匹配层区域的每个方向各引入两个完全匹配层辅助波场，所述主区域为纵波波场传播时能量不被吸收衰减的区域；

预先对所述完全匹配层区域进行参数设置。

本发明的方法，所述波动方程数值模拟所需参数包括采样网格数、三维空间方向的采样步长、三维空间方向的采样点数、时间采样点数、时间采样步长、地震子波类型和主频；所述主区域的波动方程为：

$\frac{\∂Q}{\∂t}={\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{\∂z^2})+f\left(t\right),$

其中，P为纵波波场，Q为引入的所述纵波波场P对时间变量t的一阶导数波场，υ为声波波速，x、y、z为三维空间方向，f(t)为地震子波。

本发明的方法，所述x方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}{=\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{{\∂z}^2})-{\mathrm{\υ}}^2(\frac{{\∂E}_x}{\∂x}+H_x);$

所述x方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{{\∂E}_x}{\∂t}+{\mathrm{\σ}}_x{E}_x={\mathrm{\σ}}_x\frac{\∂p}{\∂x}\\ \frac{{\∂H}_x}{\∂t}+{\mathrm{\σ}}_x{H}_x={\mathrm{\σ}}_x(\frac{{\∂}^{2}p}{{\∂x}^2}-\frac{{\∂E}_x}{\∂x})\end{array}\right.;$

所述y方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}={\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{{\∂z}^2})-{\mathrm{\υ}}^2(\frac{{\∂E}_y}{\∂y}+H_y);$

所述y方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{{\∂E}_y}{\∂t}+{\mathrm{\σ}}_y{E}_y={\mathrm{\σ}}_y\frac{\∂P}{\∂y}\\ \frac{{\∂H}_y}{\∂t}+{\mathrm{\σ}}_y{H}_y={\mathrm{\σ}}_y(\frac{{\∂}^{2}P}{{\∂y}^2}-\frac{{\∂E}_y}{\∂y})\end{array}\right.;$

所述z方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}={\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{{\∂z}^2}){-\mathrm{\υ}}^2(\frac{{\∂E}_z}{\∂z}+H_z);$

所述z方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{{\∂E}_z}{\∂t}+{\mathrm{\σ}}_z{E}_z={\mathrm{\σ}}_z\frac{\∂P}{\∂z}\\ \frac{{\∂H}_z}{\∂t}+{\mathrm{\σ}}_z{H}_z={\mathrm{\σ}}_z(\frac{{\∂}^{2}P}{\∂z^2}-\frac{\∂E_z}{\∂z})\end{array}\right.;$

其中，E_x和H_x为所述x方向引入的完全匹配层辅助波场，E_y和H_y为所述y方向引入的完全匹配层辅助波场，E_z和H_z为所述z方向引入的完全匹配层辅助波场。σ_x、σ_y和σ_z分别对应所述x方向、y方向和z方向的衰减因子。

本发明的方法，所述预先对完全匹配层区域进行参数设置包括：

设置完全匹配层吸收边界条件的层数、最大吸收衰减值和衰减函数。

本发明的方法，所述计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，具体为：

采用有限差分法或伪谱法计算x方向、y方向或z方向的二阶偏导数。

本发明的方法，所述根据所述方向的二阶偏导数，对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减，具体为：

根据所述方向的二阶偏导数，在所述整个三维空间数值模拟区域更新所述方向的纵波波场的时间一阶导数波场；

根据所述方向的二阶偏导数，在所述方向的完全匹配层区域更新所述方向的完全匹配层辅助波场；

结合更新后的所述方向的完全匹配层辅助波场，对更新后的所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。

本发明的方法，所述根据所述方向的二阶偏导数，在所述整个三维空间数值模拟区域更新所述方向的纵波波场的时间一阶导数波场，具体为：

根据x方向的二阶偏导数以及公式在所述整个三维空间数值模拟区域更新所述x方向的纵波波场的时间一阶导数波场；或，

根据y方向的二阶偏导数以及公式在所述整个三维空间数值模拟区域更新所述y方向的纵波波场的时间一阶导数波场；或，

根据z方向的二阶偏导数以及公式在所述整个三维空间数值模拟区域更新所述z方向的纵波波场的时间一阶导数波场；

其中，Δt为时间变化量。

本发明的方法，所述根据所述方向的二阶偏导数，在所述方向的完全匹配层区域更新所述方向的完全匹配层辅助波场，具体为：

根据x方向的二阶偏导数对所述E_x和所述H_x进行更新计算，具体通过公式：

和进行；或，

根据y方向的二阶偏导数对所述E_y和所述H_y进行更新计算，具体通过公式：

和进行；或，

根据z方向的二阶偏导数对所述E_z和所述H_z进行更新计算，具体通过公式：

和进行。

本发明的方法，所述结合更新后的所述方向的完全匹配层辅助波场，对更新后的所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减，具体为：

根据更新后的E_x和H_x，以及公式对更新后的x方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；或，

根据更新后的E_y和H_y，以及公式对更新后的y方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；或，

根据更新后的E_z和H_z，以及公式对更新后的z方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。

本发明的方法，在整个三维数值模拟计算区域对所述纵波波场进行更新，具体为：

在整个三维数值模拟计算区域根据公式P←P+ΔtQ对所述纵波波场进行更新。

本发明的方法，所述和均通过有限差分法计算得到。

本发明的方法，所述的完全匹配层的层数的取值范围为4～10层。

本发明的方法，所述的完全匹配层最大吸收衰减值的取值范围为300～2000。

本发明的方法，所述的完全匹配层的衰减函数为其中h代表x方向、y方向或z方向的完全匹配层距离主区域外边界的距离，d为完全匹配层的厚度，m的取值为2或3。

本发明提出的完全匹配层吸收边界条件不需要在完全匹配层区域对纵波波场进行分裂，对虚假波场的衰减在空间x、y和z方向分别进行的，由于完全匹配层中的波动方程在三个坐标轴方向是解耦的，当计算某个面部区域的时候，已经包含了对应的四个棱边区域、四个棱角区域和一个中央区域，或者说四个棱边区域、四个角点区域和一个中央区域在对应的坐标轴方向上具有相同的衰减因子参数，因此可以一起参与计算，因此不需要单独处理交叉重叠部分，从而大大降低了实现复杂度。

附图说明

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

图1为本发明实施例的基于二阶偏微分波动方程的PML吸收边界条件的方法流程图；

图2为现有技术中三维波场传播x-z方向在没有应用吸收边界条件时的波场快照；

图3为现有技术中三维波场传播x-z方向在应用传统吸收边界条件时的波场快照；

图4为三维波场传播x-z方向在应用本发明实施例的完全匹配层吸收边界条件时的波场快照；

图5为现有技术中三维波场传播y-z方向在没有应用吸收边界条件时的波场快照；

图6为现有技术中三维波场传播y-z方向在应用传统吸收边界条件时的波场快照；

图7为三维波场传播y-z方向在应用在本发明实施例的完全匹配层吸收边界条件时的波场快照；

图8为现有技术中没有应用吸收边界条件的三维合成数据单炮记录图；

图9为现有技术中使用传统吸收边界条件的三维合成数据单炮记录图；

图10使用本发明实施例的完全匹配层吸收边界条件的三维合成数据单炮记录图。

具体实施方式

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

下面结合附图，对本发明的具体实施方式作进一步的详细说明。

请参阅图1所示，本发明实施例的基于二阶偏微分波动方程的PML吸收边界条件的方法包括以下步骤：

步骤S1，预先将纵波波场的三维速度模型在三维空间上进行网格离散化，并设置波动方程数值模拟所需参数。所需参数包括采样网格数、三维空间方向的采样步长、三维空间方向的采样点数、时间采样点数、时间采样步长、地震子波类型（例如雷克子波）和主频（例如峰值频率）；所采用的波动方程为：

$\frac{{\∂}^{2}P}{{\∂t}^2}=\frac{1}{{\mathrm{\υ}}^2}(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{{\∂z}^2})+f\left(t\right)$

其中υ为各向同性介质的纵波速度，P为纵波波场，f(t)为雷克子波，x、y、z为三维空间方向。为了便于波场更新和减小计算存储，引入纵波波场P对时间变量t的一阶导数波场Q，满足如下方程：

$\frac{\∂P}{\∂t}=Q$

从而得到等价形式的波动方程：

$\frac{\∂Q}{\∂t}=\frac{1}{{\mathrm{\υ}}^2}(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{{\∂z}^2})+f\left(t\right)$

该方程是主区域的波动方程，地震波场在该区域传播时将不产生衰减。

步骤S2，预先将经过网格离散化处理后的三维速度模型的整个三维空间数值模拟区域划分为主区域和完全匹配层区域，并在完全匹配层区域的每个方向各引入两个完全匹配层辅助波场。主区域即为纵波波场传播时能量不被吸收衰减的区域。完全匹配层区域又可划分为六个部分：模型的顶部、模型的底部、模型的左侧、模型的右侧、模型的前部和模型的后部。其中：

x方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}{=\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{{\∂z}^2})-{\mathrm{\υ}}^2(\frac{{\∂E}_x}{\∂x}+H_x);$

x方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{{\∂E}_x}{\∂t}+{\mathrm{\σ}}_x{E}_x={\mathrm{\σ}}_x\frac{\∂p}{\∂x}\\ \frac{{\∂H}_x}{\∂t}+{\mathrm{\σ}}_x{H}_x={\mathrm{\σ}}_x(\frac{{\∂}^{2}p}{{\∂x}^2}-\frac{{\∂E}_x}{\∂x})\end{array}\right.;$

y方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}{=\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{{\∂z}^2})-{\mathrm{\υ}}^2(\frac{{\∂E}_y}{\∂y}+H_y);$

y方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{{\∂E}_y}{\∂t}+{\mathrm{\σ}}_y{E}_y={\mathrm{\σ}}_y\frac{\∂P}{\∂y}\\ \frac{{\∂H}_y}{\∂t}+{\mathrm{\σ}}_y{H}_y={\mathrm{\σ}}_y(\frac{{\∂}^{2}P}{{\∂y}^2}-\frac{{\∂E}_y}{\∂y})\end{array}\right.;$

z方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}={\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{{\∂x}^2}+\frac{{\∂}^{2}P}{{\∂y}^2}+\frac{{\∂}^{2}P}{{\∂z}^2}){-\mathrm{\υ}}^2(\frac{{\∂E}_z}{\∂z}+H_z);$

z方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{{\∂E}_z}{\∂t}+{\mathrm{\σ}}_z{E}_z={\mathrm{\σ}}_z\frac{\∂P}{\∂z}\\ \frac{{\∂H}_z}{\∂t}+{\mathrm{\σ}}_z{H}_z={\mathrm{\σ}}_z(\frac{{\∂}^{2}P}{\∂z^2}-\frac{\∂E_z}{\∂z})\end{array}\right.;$

其中，E_x和H_x为x方向引入的完全匹配层辅助波场，E_y和H_y为y方向引入的完全匹配层辅助波场，E_z和H_z为z方向引入的完全匹配层辅助波场。σ_x、σ_y和σ_z分别对应x方向、y方向和z方向的衰减因子。

步骤S3，预先对完全匹配层区域进行参数设置。具体包括设置完全匹配层吸收边界条件的层数、最大吸收衰减值和衰减函数。其中，完全匹配层的层数的取值范围一般设为4～10层；完全匹配层最大吸收衰减值的取值范围一般设为300～2000；完全

其中h代表x方向、y方向或z方向的完全匹配层距离主区域外边界的距离。d为完全匹配层的厚度。m的取值一般为2或3。

步骤S4，加载某一采样时刻的震源的纵波波场。所使用的公式为：P(s_x,s_y,s_z;t)=P(s_x,s_y,s_z;t)+f(t)，其中(s_x,s_y,s_z)为震源（或炮点）的空间位置。

步骤S5，计算采样时刻的纵波波场在三维空间x方向的二阶偏导数，并根据x方向的二阶偏导数对x方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。具体包括如下分步骤:

步骤S501，采用有限差分法或伪谱法计算采样时刻的纵波波场在三维空间x方向的二阶偏导数

步骤S502，根据x方向的二阶偏导数以及公式在整个三维空间数值模拟区域更新x方向的纵波波场的时间一阶导数波场，其中，Δt为时间变化量。

步骤S503，根据x方向的二阶偏导数对E_x和H_x进行更新计算，具体通过公式：和进行。

步骤S504，根据更新后E_x和H_x，以及公式对更新后的x方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。

步骤S6，计算采样时刻的纵波波场在三维空间y方向的二阶偏导数，并根据y方向的二阶偏导数对y方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。具体包括如下分步骤:

步骤S601，采用有限差分法或伪谱法计算采样时刻的纵波波场在三维空间y方向的二阶偏导数

步骤S602，根据y方向的二阶偏导数以及公式在整个三维空间数值模拟区域更新y方向的纵波波场的时间一阶导数波场，其中，Δt为时间变化量。

步骤S603，根据y方向的二阶偏导数对E_y和H_y进行更新计算，具体通过公式：和进行。

步骤S604，根据更新后的E_y和H_y，以及公式对更新后的y方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。

步骤S7，计算采样时刻的纵波波场在三维空间z方向的二阶偏导数，并根据z方向的二阶偏导数对z方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。具体包括如下分步骤:

步骤S701，采用有限差分法或伪谱法计算采样时刻的纵波波场在三维空间z方向的二阶偏导数

步骤S702，根据z方向的二阶偏导数以及公式在整个三维空间数值模拟区域更新z方向的纵波波场的时间一阶导数波场，其中，Δt为时间变化量。

步骤S703，根据z方向的二阶偏导数对E_z和H_z进行更新计算，具体通过公式：和进行。

步骤S704，根据更新后的E_z和H_z，以及公式对更新后的z方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。

步骤S8，在整个三维数值模拟计算区域根据公式P←P+ΔtQ对纵波波场进行更新。

步骤S9，重复以上步骤，在时间方向进行迭代，直至完成纵波波场所有采样时刻的处理。

本发明实施例中，涉及到在完全匹配层区域计算的一阶空间偏导数 和均可以通过有限差分法计算得到。

如果不采用任何吸收边界条件，从x-z方向传播的波场快照（参见图2）、y-z方向传播的波场快照（参见图5）以及在地面上接收的三个合成炮记录（参见图8）可以看出，截断边界处产生的虚假反射非常强烈，很难分清有效传播的波场和截断边界处产生的虚假波场，波场显得非常复杂，影响波场研究和分析。

如果采用传统的吸收边界条件，例如10个层厚的阻尼吸收衰减边界条件，从x-z方向传播的波场快照（参见图3）、y-z方向传播的波场快照（参见图6）以及在地面上接收的三个合成炮记录（参见图9）可以看出，截断边界处产生的虚假反射很大程度上得到了压制，但仍然残留着能量较强的截断边界处引起的虚假反射，影响辨认能量较弱的有效的散射波，地震记录的质量有待进一步提高。

而采用了本发明实施例中提出的完全匹配层吸收边界条件后，无论是从从x-z方向传播的波场快照（参见图4）、y-z方向传播的波场快照（参见图7）还是在地面上接收的三个合成炮记录（参见图10）可以看出，几乎察觉不到任何由于截断边界引起的虚假反射，地震记录的质量得到很大程度的改善，由此可见，本发明实施例提出的吸收边界条件对虚假反射的衰减效果非常好。

我们知道，传统的分裂式完全匹配层吸收边界条件的实现涉及到完全匹配层六个面的交叉区域：两个面的交叉区域即棱边区域和三个面的交叉区域即棱角区域。每个面部区域由四个棱角区域、四个棱边区域和一个中央区域组成。由于每个区域的衰减因子参数不相同，在计算某个面部区域时候需要分别讨论四个棱边区域、四个棱角区域和一个中央区域；纵波波场的三维速度模型共涉及到十二个棱边区域、八个棱角区域和六个面部的中央区域，总共需要讨论26种情形，编程实现起来比较复杂。本发明实施例所使用的完全匹配层吸收边界条件在三个坐标轴方向是解耦的、相互不受影响，在计算机编程数值实现过程中，当计算某个面部区域的时候，已经包含了对应的四个棱边区域、四个棱角区域和一个中央区域，或者说四个棱边区域、四个角点区域和一个中央区域在对应的坐标轴方向具有相同的衰减因子参数，因此可以一起参与计算。所以，本文的完全匹配层吸收边界条件不涉及到交叉区域的单独计算，总共只需要讨论6种情形，这大大降低了编程实现复杂性。总之，本发明实施例基于二阶偏微分波动方程的非分列式方案与现有技术中基于一阶偏微分波动方程的分列式方案相比，其计算效率较高，实现复杂度较低。

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

Claims (12)

1.一种基于二阶偏微分波动方程的PML吸收边界条件的方法，其特征在于，包括以下步骤：

预先将纵波波场的三维速度模型在三维空间上进行网格离散化，并设置波动方程数值模拟所需参数；

预先将经过网格离散化处理后的三维速度模型的整个三维空间数值模拟区域划分为主区域和完全匹配层区域，并在所述完全匹配层区域的每个方向各引入两个完全匹配层辅助波场，所述主区域为纵波波场传播时能量不被吸收衰减的区域；

预先对所述完全匹配层区域进行参数设置；

加载某一采样时刻的震源的纵波波场；

计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，并根据所述方向的二阶偏导数对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；

重复上一步骤，分别对所述采样时刻的纵波波场在三维空间另外两个方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；

在整个三维数值模拟计算区域对所述纵波波场进行更新；

重复以上步骤，直至完成所述纵波波场所有采样时刻的处理；其中：

所述波动方程数值模拟所需参数包括采样网格数、三维空间方向的采样步长、三维空间方向的采样点数、时间采样点数、时间采样步长、地震子波类型和主频；所述主区域的波动方程为：

$\frac{\∂Q}{\∂t}={\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{\∂x^2}+\frac{{\∂}^{2}P}{\∂y^2}+\frac{{\∂}^{2}P}{\∂z^2})+f\left(t\right),$

其中，P为纵波波场，Q为引入的所述纵波波场P对时间变量t的一阶导数波场，υ为声波波速，x、y、z为三维空间方向，f(t)为地震子波；

所述x方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}={\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{\∂x^2}+\frac{{\∂}^{2}P}{\∂y^2}+\frac{{\∂}^{2}P}{\∂z^2})-{\mathrm{\υ}}^2(\frac{\∂E_x}{\∂x}+H_x);$

所述x方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{\∂E_x}{\∂t}+{\mathrm{\σ}}_x{E}_x={\mathrm{\σ}}_x\frac{\∂p}{\∂x}\\ \frac{\∂H_x}{\∂t}+{\mathrm{\σ}}_x{H}_x={\mathrm{\σ}}_x(\frac{{\∂}^{2}p}{\∂x^2}-\frac{\∂E_x}{\∂x})\end{array}\right.;$

所述y方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}={\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{\∂x^2}+\frac{{\∂}^{2}P}{\∂y^2}+\frac{{\∂}^{2}P}{\∂z^2})-{\mathrm{\υ}}^2(\frac{\∂E_y}{\∂y}+H_y);$

所述y方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{\∂E_y}{\∂t}+{\mathrm{\σ}}_y{E}_y={\mathrm{\σ}}_y\frac{\∂P}{\∂y}\\ \frac{\∂H_y}{\∂t}+{\mathrm{\σ}}_y{H}_y={\mathrm{\σ}}_y(\frac{{\∂}^{2}P}{\∂y^2}-\frac{\∂E_y}{\∂y})\end{array}\right.;$

所述z方向完全匹配层区域的波动方程为：

$\frac{\∂Q}{\∂t}={\mathrm{\υ}}^2(\frac{{\∂}^{2}P}{\∂x^2}+\frac{{\∂}^{2}P}{\∂y^2}+\frac{{\∂}^{2}P}{\∂z^2})-{\mathrm{\υ}}^2(\frac{\∂E_z}{\∂z}+H_z);$

所述z方向完全匹配层区域的辅助方程为：

$\left\{\begin{array}{c}\frac{\∂E_z}{\∂t}+{\mathrm{\σ}}_z{E}_z={\mathrm{\σ}}_z\frac{\∂P}{\∂z}\\ \frac{\∂H_z}{\∂t}+{\mathrm{\σ}}_z{H}_z={\mathrm{\σ}}_z(\frac{{\∂}^{2}P}{\∂z^2}-\frac{\∂E_z}{\∂z})\end{array}\right.;$

其中，E_x和H_x为所述x方向引入的完全匹配层辅助波场，E_y和H_y为所述y方向引入的完全匹配层辅助波场，E_z和H_z为所述z方向引入的完全匹配层辅助波场；σ_x、σ_y和σ_z分别对应所述x方向、y方向和z方向的衰减因子。

2.根据权利要求1所述的方法，其特征在于，所述预先对完全匹配层区域进行参数设置包括：

设置完全匹配层吸收边界条件的层数、最大吸收衰减值和衰减函数。

3.根据权利要求1所述的方法，其特征在于，所述计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，具体为：

采用有限差分法或伪谱法计算x方向、y方向或z方向的二阶偏导数。

4.根据权利要求3所述的方法，其特征在于，所述根据所述方向的二阶偏导数，对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减，具体为：

根据所述方向的二阶偏导数，在所述整个三维空间数值模拟区域更新所述方向的纵波波场的时间一阶导数波场；

根据所述方向的二阶偏导数，在所述方向的完全匹配层区域更新所述方向的完全匹配层辅助波场；

结合更新后的所述方向的完全匹配层辅助波场，对更新后的所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。

5.根据权利要求4所述的方法，其特征在于，所述根据所述方向的二阶偏导数，在所述整个三维空间数值模拟区域更新所述方向的纵波波场的时间一阶导数波场，具体为：

根据x方向的二阶偏导数以及公式在所述整个三维空间数值模拟区域更新所述x方向的纵波波场的时间一阶导数波场；或，

根据y方向的二阶偏导数以及公式在所述整个三维空间数值模拟区域更新所述y方向的纵波波场的时间一阶导数波场；或，

根据z方向的二阶偏导数以及公式在所述整个三维空间数值模拟区域更新所述z方向的纵波波场的时间一阶导数波场；

其中，Δt为时间变化量。

6.根据权利要求4所述的方法，其特征在于，所述根据所述方向的二阶偏导数，在所述方向的完全匹配层区域更新所述方向的完全匹配层辅助波场，具体为：

根据x方向的二阶偏导数对所述E_x和所述H_x进行更新计算，具体通过公式：

$E_x\←e^{-{\mathrm{\σ}}_x\mathrm{\Δ}t}{E}_x+(1-e^{-{\mathrm{\σ}}_x\mathrm{\Δ}t})\frac{\∂P}{\∂x}$ 和 $H_x\←e^{-{\mathrm{\σ}}_x\mathrm{\Δ}t}{H}_x+(1-e^{-{\mathrm{\σ}}_x\mathrm{\Δ}t})(\frac{{\∂}^{2}P}{\∂x^2}-\frac{\∂E_x}{\∂x})$ 进行；或，

根据y方向的二阶偏导数对所述E_y和所述H_y进行更新计算，具体通过公式：

$E_y\←e^{-{\mathrm{\σ}}_y\mathrm{\Δ}t}{E}_y+(1-e^{-{\mathrm{\σ}}_y\mathrm{\Δ}t})\frac{\∂P}{\∂y}$ 和 $H_y\←e^{-{\mathrm{\σ}}_y\mathrm{\Δ}t}{H}_y+(1-e^{-{\mathrm{\σ}}_y\mathrm{\Δ}t})(\frac{{\∂}^{2}P}{\∂y^2}-\frac{\∂E_y}{\∂y})$ 进行；或，

根据z方向的二阶偏导数对所述E_z和所述H_z进行更新计算，具体通过公式：

和 $H_z\←e^{-{\mathrm{\σ}}_z\mathrm{\Δ}t}{H}_z+(1-e^{-{\mathrm{\σ}}_z\mathrm{\Δ}t})(\frac{{\∂}^{2}P}{\∂z^2}-\frac{\∂E_z}{\∂z})$ 进行。

7.根据权利要求4所述的方法，其特征在于，所述结合更新后的所述方向的完全匹配层辅助波场，对更新后的所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减，具体为：

根据更新后的E_x和H_x，以及公式对更新后的x方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；或，

根据更新后的E_y和H_y，以及公式对更新后的y方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；或，

根据更新后的E_z和H_z，以及公式对更新后的z方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。

8.根据权利要求4所述的方法，其特征在于，在整个三维数值模拟计算区域对所述纵波波场进行更新，具体为：

在整个三维数值模拟计算区域根据公式P←P+ΔtQ对所述纵波波场进行更新。

9.根据权利要求6所述的方法，其特征在于，所述和均通过有限差分法计算得到。

10.根据权利要求2所述的方法，其特征在于，所述的完全匹配层的层数的取值范围为4～10层。

11.根据权利要求2所述的方法，其特征在于，所述的完全匹配层最大吸收衰减值的取值范围为300～2000。

12.根据权利要求2所述的方法，其特征在于，所述的完全匹配层的衰减函数为其中h代表x方向、y方向或z方向的完全匹配层距离主区域外边界的距离，d为完全匹配层的厚度，m的取值为2或3。