CN105447225B - 一种应用于声波有限差分数值模拟的组合吸收边界条件 - Google Patents

Abstract

一种应用于声波有限差分数值模拟的组合吸收边界条件，属于地震勘探数值模拟领域，具体包括以下步骤：基于2N(N>0)阶精度交错网格有限差分格式进行声波方程数值模拟时，首先在人工截断边界处设置L(L>N)层完全匹配层(PML)，应用PML边界条件吸收来自中心波场的边界反射波；然后对于PML外部的N层边界应用Higdon三阶吸收边界条件，以吸收PML的外边界反射。本发明方法充分利用PML边界条件和Higdon三阶吸收边界条件两者的优势，能够有效吸收人工边界内层和外层的边界反射，从而实现了高精度的有限差分数值模拟。

Description

一种应用于声波有限差分数值模拟的组合吸收边界条件

技术领域

本发明属于地震勘探数值模拟领域，具体地涉及一种应用于声波有限差分数值模拟的组合吸收边界条件。

背景技术

地震波正演模拟是一种通过将真实的地下地层介质简化为数学模型，然后用数值计算方法模拟地震波在模型中传播的方法。数值模拟是了解地震波在介质中传播规律、帮助识别实测数据中有效信息的重要手段。在基于计算机实现的地震勘探数值模拟中，受限于当前计算机的存储能力和计算能力，对于无限空间中的地质模型需要引入人工截断边界来界定计算区域，但是简单地人为截断将导致边界处产生强烈的反射干扰，因此在地震波正演模拟时通常都需要对人工边界进行特殊处理以消除人工边界的虚假反射。

目前常用的消除人工边界反射的处理方法主要有两类：第一类是最佳匹配层(PML)边界条件方法，第二类是基于单程波近似的吸收边界条件方法。基于最佳匹配层的边界条件其主要思想是在人工截断边界上设置一定厚度的“匹配层”，层内引入衰减因子，使得地震波在此区域内传播时能量迅速衰减以消除人工边界反射。基于单程波近似的吸收边界条件则是从常规双程波动方程中分解出外行波方程(只能描述外行波的传播规律)，并将外行波方程置于人工边界区域以此作为吸收边界条件，从而达到消除边界反射的目的。

PML边界条件理论上能够有效吸收各个角度的边界入射波，具有较高的边界反射吸收效率，但PML最外层是一个强反射界面，波场模拟时其最外层边界仍然会产生部分残余边界反射，从而在一定程度上影响波场模拟的精度。基于单程波近似的吸收边界条件方法主要包括Clayton-Enquist吸收边界方法和Higdon吸收边界方法，其吸收效率与边界条件的阶数有关。高阶的Clayton-Enquist吸收边界方法难以实现差分计算，因此在实际模拟时，Clayton-Enquist吸收边界方法通常只能应用到二阶吸收边界条件，其对大角度入射波的吸收效果较差。

发明内容

本发明要解决的技术问题在于提供一种应用于声波有限差分数值模拟的组合吸收边界条件，本发明充分利用PML边界条件和Higdon三阶吸收边界条件的优势(综合考虑边界反射吸收效果和计算效率，本发明采用Higdon三阶吸收边界条件，将PML边界条件与Higdon吸收边界条件有机结合，在不增加模拟波场范围大小的情况下，显著提高了人工边界反射的吸收效率。

本发明采取以下技术方案：

一种应用于声波有限差分数值模拟的组合吸收边界条件，具体包括以下步骤：

(1)在数值模拟的中心波场区域利用声波方程交错网格有限差分方法进行波场计算，确定中心波场的计算区域大小、观测系统信息；有限差分格式的空间精度为2N阶，其中N>0，N取值3或4；在中心波场的人工截断边界处构造L层完全匹配层，其中L>N，应用PML边界条件吸收来自中心波场的边界反射波；PML内的阻尼因子d(s)采用四阶指数型吸收衰减因子，其表达式为：

$d\left(s\right)=log\left(\frac{1}{R}\right)\·\frac{5V_p}{2L}{\left(\frac{L-s}{L}\right)}^4$

其中，s为PML内的计算点到PML最外层边界的距离，R为理论反射系数，V_p为地震波传播速度，L为PML厚度；

(2)在PML外部的N层边界应用Higdon三阶吸收边界条件，以吸收PML的外边界反射，其左边界的表达式为：

$B_{3}P=({\mathrm{cos\α}}_1\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})({\mathrm{cos\α}}_2\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})({\mathrm{cos\α}}_3\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})P=0$

其中，cosα_j(j＝1,2,3)为入射角度，v为波速，P为声波方程中的压强值；

上式有限差分格式为：

$\begin{array}{c}{P}_{m,n}^{k+1}=(r_{1,2,3}^2{P}_{m+1,n}^{k+1}+r_{1,2,3}^3{P}_{m,n}^k+r_{1,2,3}^4{P}_{m+1,n}^k+r_{1,2,3}^5{P}_{m+2,n}^{k+1}+r_{1,2,3}^6{P}_{m,n}^{k-1}+r_{1,2,3}^7{P}_{m+2,n}^k+r_{1,2,3}^8{P}_{m+1,n}^{k-1}+r_{1,2,3}^9{P}_{m+2,n}^{k-1}\\ +r_{1,2,3}^{10}{P}_{m+3,n}^{k+1}+r_{1,2,3}^{11}{P}_{m,n}^{k-2}+r_{1,2,3}^{12}{P}_{m+3,n}^k+r_{1,2,3}^{13}{P}_{m+3,n}^{k-1}+r_{1,2,3}^{14}{P}_{m+1,n}^{k-2}+r_{1,2,3}^{15}{P}_{m+2,n}^{k-2}+r_{1,2,3}^{16}{P}_{m+3,n}^{k-2})/(-r_{1,2,3}^1)\end{array}$

其中，m,n为空间离散网格点坐标，m＝0,1,…,N-1,N，n＝0,1,…,N_z-1,N_z，N_z为中心波场的纵向网格点数，左边界表达式中各系数表达式如下：

$\begin{array}{ccc}\begin{array}{c}{r}_{1,2,3}^1=l_{1,1}{s}_{2,3}^1\\ r_{1,2,3}^2=l_{1,1}{s}_{2,3}^2+l_{1,2}{s}_{2,3}^1\\ r_{1,2,3}^3=l_{1,1}{s}_{2,3}^3+l_{1,3}{s}_{2,3}^1\\ r_{1,2,3}^4=l_{1,1}{s}_{2,3}^4+l_{1,2}{s}_{2,3}^3+l_{1,3}{s}_{2,3}^2+l_{1,4}{s}_{2,3}^1\\ r_{1,2,3}^5=l_{1,1}{s}_{2,3}^5+l_{1,2}{s}_{2,3}^2\\ r_{1,2,3}^6=l_{1,1}{s}_{2,3}^6+l_{1,3}{s}_{2,3}^3\\ r_{1,2,3}^7=l_{1,1}{s}_{2,3}^7+l_{1,2}{s}_{2,3}^4+l_{1,3}{s}_{2,3}^5+l_{1,4}{s}_{2,3}^2\\ r_{1,2,3}^8=l_{1,1}{s}_{2,3}^8+l_{1,2}{s}_{2,3}^6+l_{1,3}{s}_{2,3}^4+l_{1,4}{s}_{2,3}^3\end{array}& \begin{array}{c}{r}_{1,2,3}^9=l_{1,1}{s}_{2,3}^9+l_{1,2}{s}_{2,3}^8+l_{1,3}{s}_{2,3}^7+l_{1,4}{s}_{2,3}^4\\ r_{1,2,3}^{10}=l_{1,2}{s}_{2,3}^5\\ r_{1,2,3}^{11}=l_{1,3}{s}_{2,3}^6\\ r_{1,2,3}^{12}=l_{1,2}{s}_{2,3}^7+l_{1,4}{s}_{2,3}^5\\ r_{1,2,3}^{13}=l_{1,2}{s}_{2,3}^9+l_{1,4}{s}_{2,3}^7\\ r_{1,2,3}^{14}=l_{1,3}{s}_{2,3}^8+l_{1,4}{s}_{2,3}^6\\ r_{1,2,3}^{15}=l_{1,3}{s}_{2,3}^9+l_{1,4}{s}_{2,3}^8\\ r_{1,2,3}^{16}=l_{1,4}{s}_{2,3}^9\end{array}& \begin{array}{c}{s}_{2,3}^1=l_{2,1}{l}_{3,1}\\ s_{2,3}^2=l_{2,1}{l}_{3,2}+l_{2,2}{l}_{3,1}\\ s_{2,3}^2=l_{2,1}{l}_{3,3}+l_{2,3}{l}_{3,1}\\ s_{2,3}^4=l_{2,1}{l}_{3,4}+l_{2,2}{l}_{3,3}+l_{2,3}{l}_{3,2}+l_{2,4}{l}_{3,1}\\ s_{2,3}^5=l_{2,2}{l}_{3,2}\\ s_{2,3}^6=l_{2,3}{l}_{3,3}\\ s_{2,3}^7=l_{2,2}{l}_{3,4}+l_{2,4}{l}_{3,2}\\ s_{2,3}^8=l_{2,3}{l}_{3,4}+l_{2,4}{l}_{3,3}\\ s_{2,3}^9=l_{2,4}{l}_{3,4}\end{array}\end{array}$

其中，Δt、Δh分别为时间、空间采样间隔；右边界、上边界和下边界类比以上推导过程得出。

本发明与现有技术相比的有益效果：

本发明采取组合吸收边界条件，将四阶指数型吸收衰减因子的PML边界条件与Higdon三阶吸收边界条件进行组合，作为声波方程有限差分数值模拟时的边界条件，即在PML外部的N层(有限差分格式的空间精度为2N)网格点上使用Higdon三阶吸收边界条件。与单纯使用PML边界条件相比，一方面PML+Higdon组合吸收边界条件能够集合二者的优势，利用Higdon吸收边界条件的特点对PML外边界的强反射进行再吸收，总体上加强了对边界反射的吸收效果，从而提高了数值模拟的精度；另一方面由于应用PML+Higdon组合吸收边界条件所设置的Higdon吸收边界为PML外部因有限差分计算而引入的N层网格点，使用该边界条件时整个波场的大小与仅使用PML边界条件时相等，因而不增加整个波场的内存消耗。

附图说明

图1为PML和Higdon吸收区域的组合示意图；

图2为应用PML边界条件时小角度入射情况下的波场；

图3为应用PML+Higdon组合吸收边界条件时小角度入射情况下的波场；

图4为分别应用PML边界条件与PML+Higdon组合吸收边界条件时，在4000m深度处波场值曲线，实线对应PML边界条件时的波场值曲线，虚线对应PML+Higdon组合吸收边界条件时的波场值曲线。

图5为应用PML边界条件时大角度入射情况下的波场；

图6为应用PML+Higdon组合吸收边界条件时大角度入射情况下的波场；

图7为分别应用PML边界条件与PML+Higdon组合吸收边界条件时，在7500m深度处波场值曲线，实线对应PML边界条件时的波场值曲线，虚线对应PML+Higdon组合吸收边界条件时的波场值曲线。

具体实施方式

下面通过实施例结合附图来对本发明的技术方案作进一步解释，但本发明的保护范围不受实施例任何形式上的限制。

实施例1

本发明采取组合吸收边界条件，在声波方程有限差分计算的过程中，将四阶指数型吸收衰减因子的PML边界条件与Higdon三阶吸收边界条件进行组合，在人工截断边界处构造组合边界区域，即在PML外部的N层(有限差分格式的空间精度为2N)网格点上使用Higdon三阶吸收边界条件，显著提高了正演模拟的边界吸收效果。

本发明的主要实施过程分为两个步骤：利用声波方程进行交错网格有限差分数值模拟，在人工截断边界处构造L层PML，应用PML边界条件吸收来自中心波场的边界反射波；在PML外部的N层区域应用Higdon三阶吸收边界条件以吸收PML的外边界反射。其具体步骤如下：

(1)确定中心波场的计算区域大小、观测系统信息等，利用基于声波方程的交错网格有限差分方法实现中心波场的数值模拟，其交错网格有限差分格式为：

$\left\{\begin{array}{c}{v}_x{|}_{i,j+\frac{1}{2}}^{k+\frac{1}{2}}=v_x{|}_{i,j+\frac{1}{2}}^{k-\frac{1}{2}}+\frac{1}{{\mathrm{\ρ}}_{i,j+\frac{1}{2}}}\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z}\underset{n=1}{\overset{N}{\Σ}}{a}_n(P_{i,j+n}^k+P_{i,j-n+1}^k)\\ v_z{|}_{i+\frac{1}{2},j}^{k+\frac{1}{2}}=v_z{|}_{i+\frac{1}{2},j}^{k-\frac{1}{2}}+\frac{1}{{\mathrm{\ρ}}_{i+\frac{1}{2},j}}\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z}\underset{n=1}{\overset{N}{\Σ}}{a}_n(P_{i+n,j}^k-P_{i-n+1,j}^k)\\ P{|}_{i,j}^{k+1}=P{|}_{i,j}^k+K{|}_{i,j}^{k+\frac{1}{2}}\left[\begin{array}{c}\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z}\underset{n=1}{\overset{N}{\Σ}}{a}_n(v_z{|}_{i+n-\frac{1}{2},j}^{k+\frac{1}{2}}-v_z{|}_{i-n+\frac{1}{2},j}^{k+\frac{1}{2}})\\ +\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\underset{n=1}{\overset{N}{\Σ}}{a}_n(v_x{|}_{i,j+n-\frac{1}{2}}^{k+\frac{1}{2}}-v_z{|}_{i,j-n+\frac{1}{2}}^{k+\frac{1}{2}})\end{array}\right]\end{array}\right.$

其中，表示在空间规则离散节点i,j、时间规则离散节点k上的压强值，表示在空间x方向上位于规则的j节点、在z方向位于的半节点、时间的半节点上的速度分量，表示在空间x方向上位于的半节点、在z方向位于规则的i节点、时间的半节点上的速度分量，K为介质的物性参数，K＝ρc²，ρ为介质的密度，c为声波在介质中传播的速度，一阶导数差分系数a_n的表达式为：

$a_n=\frac{{(-1)}^{n+1}\underset{i=1,i\≠n}{\overset{N}{\Π}}{(2i-1)}^2}{(2n-1)\underset{i=1,i\≠n}{\overset{N}{\Π}}|{(2n-1)}^2-{(2i-1)}^2|},n=1,2,\mathrm{...},N$

在中心波场的人工截断边界构建组合吸收边界区域，如图1所示，图中L层PML利用PML边界条件进行计算，衰减因子采用四阶指数型吸收衰减因子，其表达式为：

$d\left(s\right)=log\left(\frac{1}{R}\right)\·\frac{5V_p}{2L}{\left(\frac{L-s}{L}\right)}^4$

其中，s为PML区域内的计算点到其最外层边界的距离，R为理论反射系数(R＝0.001)，V_p为地震波传播速度，L为PML厚度。在中心波场区域d_x＝0，d_z＝0；在PML的上下部分d_x＝0，d_z＝d(s)；在PML的左右部分d_x＝d(s)，d_z＝0；在PML的角点区域d_x＝d(s)，d_z＝d(s)；

(2)在PML外部的N层网格点上采用Higdon三阶吸收边界条件对入射到PML外边界上的波进行吸收衰减。以Higdon三阶吸收条件的左边界为例，其表达式为：

$B_{3}P=({\mathrm{cos\α}}_1\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})({\mathrm{cos\α}}_2\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})({\mathrm{cos\α}}_3\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})P=0$

其中，cosα_j(j＝1,2,3)为入射角度，v为波速，P为声波方程中的压强值；上式有限差分格式为：

$\begin{array}{c}{P}_{m,n}^{k+1}=(r_{1,2,3}^2{P}_{m+1,n}^{k+1}+r_{1,2,3}^3{P}_{m,n}^k+r_{1,2,3}^4{P}_{m+1,n}^k+r_{1,2,3}^5{P}_{m+2,n}^{k+1}+r_{1,2,3}^6{P}_{m,n}^{k-1}+r_{1,2,3}^7{P}_{m+2,n}^k+r_{1,2,3}^8{P}_{m+1,n}^{k-1}+r_{1,2,3}^9{P}_{m+2,n}^{k-1}\\ +r_{1,2,3}^{10}{P}_{m+3,n}^{k+1}+r_{1,2,3}^{11}{P}_{m,n}^{k-2}+r_{1,2,3}^{12}{P}_{m+3,n}^k+r_{1,2,3}^{13}{P}_{m+3,n}^{k-1}+r_{1,2,3}^{14}{P}_{m+1,n}^{k-2}+r_{1,2,3}^{15}{P}_{m+2,n}^{k-2}+r_{1,2,3}^{16}{P}_{m+3,n}^{k-2})/(-r_{1,2,3}^1)\end{array}$

其中，m,n为空间离散网格点坐标，m＝0,1,…,N-1,N，n＝0,1,…,N_z-1,N_z，N_z为中心波场的纵向网格点数，左边界表达式中各系数表达式如下：

$\begin{array}{ccc}\begin{array}{c}{r}_{1,2,3}^1=l_{1,1}{s}_{2,3}^1\\ r_{1,2,3}^2=l_{1,1}{s}_{2,3}^2+l_{1,2}{s}_{2,3}^1\\ r_{1,2,3}^3=l_{1,1}{s}_{2,3}^3+l_{1,3}{s}_{2,3}^1\\ r_{1,2,3}^4=l_{1,1}{s}_{2,3}^4+l_{1,2}{s}_{2,3}^3+l_{1,3}{s}_{2,3}^2+l_{1,4}{s}_{2,3}^1\\ r_{1,2,3}^5=l_{1,1}{s}_{2,3}^5+l_{1,2}{s}_{2,3}^2\\ r_{1,2,3}^6=l_{1,1}{s}_{2,3}^6+l_{1,3}{s}_{2,3}^3\\ r_{1,2,3}^7=l_{1,1}{s}_{2,3}^7+l_{1,2}{s}_{2,3}^4+l_{1,3}{s}_{2,3}^5+l_{1,4}{s}_{2,3}^2\\ r_{1,2,3}^8=l_{1,1}{s}_{2,3}^8+l_{1,2}{s}_{2,3}^6+l_{1,3}{s}_{2,3}^4+l_{1,4}{s}_{2,3}^3\end{array}& \begin{array}{c}{r}_{1,2,3}^9=l_{1,1}{s}_{2,3}^9+l_{1,2}{s}_{2,3}^8+l_{1,3}{s}_{2,3}^7+l_{1,4}{s}_{2,3}^4\\ r_{1,2,3}^{10}=l_{1,2}{s}_{2,3}^5\\ r_{1,2,3}^{11}=l_{1,3}{s}_{2,3}^6\\ r_{1,2,3}^{12}=l_{1,2}{s}_{2,3}^7+l_{1,4}{s}_{2,3}^5\\ r_{1,2,3}^{13}=l_{1,2}{s}_{2,3}^9+l_{1,4}{s}_{2,3}^7\\ r_{1,2,3}^{14}=l_{1,3}{s}_{2,3}^8+l_{1,4}{s}_{2,3}^6\\ r_{1,2,3}^{15}=l_{1,3}{s}_{2,3}^9+l_{1,4}{s}_{2,3}^8\\ r_{1,2,3}^{16}=l_{1,4}{s}_{2,3}^9\end{array}& \begin{array}{c}{s}_{2,3}^1=l_{2,1}{l}_{3,1}\\ s_{2,3}^2=l_{2,1}{l}_{3,2}+l_{2,2}{l}_{3,1}\\ s_{2,3}^2=l_{2,1}{l}_{3,3}+l_{2,3}{l}_{3,1}\\ s_{2,3}^4=l_{2,1}{l}_{3,4}+l_{2,2}{l}_{3,3}+l_{2,3}{l}_{3,2}+l_{2,4}{l}_{3,1}\\ s_{2,3}^5=l_{2,2}{l}_{3,2}\\ s_{2,3}^6=l_{2,3}{l}_{3,3}\\ s_{2,3}^7=l_{2,2}{l}_{3,4}+l_{2,4}{l}_{3,2}\\ s_{2,3}^8=l_{2,3}{l}_{3,4}+l_{2,4}{l}_{3,3}\\ s_{2,3}^9=l_{2,4}{l}_{3,4}\end{array}\end{array}$

其中，Δt、Δh分别为时间、空间采样间隔。下面给出右边界、上边界和下边界的有限差分格式。

右边界表达式：

$\begin{array}{c}{P}_{Nx+N+2L+m,n}^{k+1}=(r_{1,2,3}^2{P}_{Nx+N+2L+m-1,n}^{k+1}+r_{1,2,3}^3{P}_{Nx+N+2L+m,n}^k+r_{1,2,3}^4{P}_{Nx+N+2L+m-1,n}^k+r_{1,2,3}^5{P}_{Nx+N+2L+m-2,n}^{k+1}\\ +r_{1,2,3}^6{P}_{Nx+N+2L+m,n}^{k-1}+r_{1,2,3}^7{P}_{Nx+N+2L+m-2,n}^k+r_{1,2,3}^8{P}_{Nx+N+2L+m-1,n}^{k-1}+r_{1,2,3}^9{P}_{Nx+N+2L+m-2,n}^{k-1}\\ +r_{1,2,3}^{10}{P}_{Nx+N+2L+m-3,n}^{k+1}+r_{1,2,3}^{11}{P}_{Nx+N+2L+m,n}^{k-2}+r_{1,2,3}^{12}{P}_{Nx+N+2L+m-3,n}^k+r_{1,2,3}^{13}{P}_{Nx+N+2L+m-3,n}^{k-1}\\ +r_{1,2,3}^{14}{P}_{Nx+N+2L+m-1,n}^{k-2}+r_{1,2,3}^{15}{P}_{Nx+N+2L+m-2,n}^{k-2}+r_{1,2,3}^{16}{P}_{Nx+N+2L+m-3,n}^{k-2})/(-r_{1,2,3}^1)\end{array}$

其中，m,n为空间离散网格点坐标，m＝0,1,…,N-1,N，n＝0,1,…,N_z-1,N_z，N_z为中心波场的纵向网格点数，L为PML层数，各系数表达式与左边界相同；

上边界表达式：

$\begin{array}{c}{P}_{m,n}^{k+1}=(r_{1,2,3}^2{P}_{m,n+1}^{k+1}+r_{1,2,3}^3{P}_{m,n}^k+r_{1,2,3}^4{P}_{m,n+1}^k+r_{1,2,3}^5{P}_{m,n+2}^{k+1}+r_{1,2,3}^6{P}_{m,n}^{k-1}+r_{1,2,3}^7{P}_{m,n+2}^k+r_{1,2,3}^8{P}_{m,n+1}^{k-1}+r_{1,2,3}^9{P}_{m,n+2}^{k-1}\\ +r_{1,2,3}^{10}{P}_{m,n+3}^{k+1}+r_{1,2,3}^{11}{P}_{m,n}^{k-2}+r_{1,2,3}^{12}{P}_{m,n+3}^k+r_{1,2,3}^{13}{P}_{m,n+3}^{k-1}+r_{1,2,3}^{14}{P}_{m,n+1}^{k-2}+r_{1,2,3}^{15}{P}_{m,n+2}^{k-2}+r_{1,2,3}^{16}{P}_{m,n+3}^{k-2})/(-r_{1,2,3}^1)\end{array}$

其中，m,n为空间离散网格点坐标，m＝0,1,…,N_x-1,N_x，n＝0,1,…,N-1,N，N_x为中心波场的横向网格点数，各系数表达式与左边界相同；

下边界表达式为：

$\begin{array}{c}{P}_{m,Nz+N+2L+n}^{k+1}=(r_{1,2,3}^2{P}_{m,Nz+N+2L+n-1}^{k+1}+r_{1,2,3}^3{P}_{m,Nz+N+2L+n}^k+r_{1,2,3}^4{P}_{m,Nz+N+2L+n-1}^k+r_{1,2,3}^5{P}_{m,Nz+N+2L+n-2}^{k+1}\\ +r_{1,2,3}^6{P}_{m,Nz+N+2L+n}^{k+1}+r_{1,2,3}^7{P}_{m,Nz+N+2L+n-2}^k+r_{1,2,3}^8{P}_{m,Nz+N+2L+n-1}^{k-1}+r_{1,2,3}^9{P}_{m,Nz+N+2L+n-2}^{k-1}\\ +r_{1,2,3}^{10}{P}_{m,Nz+N+2L+n-3}^{k+1}+r_{1,2,3}^{11}{P}_{m,Nz+N+2L+n}^{k-2}+r_{1,2,3}^{12}{P}_{m,Nz+N+2L+n-3}^k+r_{1,2,3}^{13}{P}_{m,Nz+N+2L+n-3}^{k-1}\\ +r_{1,2,3}^{14}{P}_{m,Nz+N+2L+n-1}^{k-2}+r_{1,2,3}^{15}{P}_{m,Nz+N+2L+n-2}^{k-2}+r_{1,2,3}^{16}{P}_{m,Nz+N+2L+n-3}^{k-2})/(-r_{1,2,3}^1)\end{array}$

其中，m,n为空间离散网格点坐标，m＝0,1,…,N_x-1,N_x，n＝0,1,…,N-1,N，N_x为中心波场的横向网格点数，L为PML层数，各系数表达式与左边界相同。

通过对比图2、3以及图5、6可以看出，使用组合吸收边界条件后中心波场内的边界反射比仅使用PML边界条件时明显要弱。结合图4和图7所提供的单道数据对比可以看出，无论小角度边界入射还是大角度边界入射，当使用组合吸收边界条件时边界反射的幅值要小于仅使用PML边界条件的情况。因此，在构造相同大小的吸收层的前提下，本发明的一种应用于声波方程有限差分数值模拟的组合吸收边界条件的边界吸收效果相比仅使用PML边界条件的效果得到了显著提高，应用该方法能够有效减少声波方程有限差分数值模拟中的边界反射。

Claims (1)

1.一种应用于声波有限差分数值模拟的组合吸收边界条件的方法，其特征在于所述的方法具体包括以下步骤：

(1)在数值模拟的中心波场区域利用声波方程交错网格有限差分方法进行波场计算，确定中心波场的计算区域大小、观测系统信息；有限差分格式的空间精度为2N阶，其中N>0，N取值3或4；在中心波场的人工截断边界处构造L层完全匹配层，其中L>N，应用PML边界条件吸收来自中心波场的边界反射波；PML内的阻尼因子d(s)采用四阶指数型吸收衰减因子，其表达式为：

$d\left(s\right)=log\left(\frac{1}{R}\right)\·\frac{5V_p}{2L}{\left(\frac{L-s}{L}\right)}^4$

其中，s为PML内的计算点到PML最外层边界的距离，R为理论反射系数，V_p为地震波传播速度，L为PML厚度；

(2)在PML外部的N层边界应用Higdon三阶吸收边界条件，以吸收PML的外边界反射，其左边界的表达式为：

$B_{3}P=({\mathrm{cos\α}}_1\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})({\mathrm{cos\α}}_2\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})({\mathrm{cos\α}}_3\frac{\∂}{\∂t}-v\frac{\∂}{\∂x})P=0$

其中，cosα_j(j＝1,2,3)为入射角度，v为波速，P为声波方程中的压强值，B₃P表示变量P应用Higdon三阶吸收边界条件时左边界的表达式；

上式有限差分格式为：

$\begin{array}{c}{P}_{m,n}^{k+1}=(r_{1,2,3}^2{P}_{m+1,n}^{k+1}+r_{1,2,3}^3{P}_{m,n}^k+r_{1,2,3}^4{P}_{m+1,n}^k+r_{1,2,3}^5{P}_{m+2,n}^{k+1}+r_{1,2,3}^6{P}_{m,n}^{k-1}+r_{1,2,3}^7{P}_{m+2,n}^k+r_{1,2,3}^8{P}_{m+1,n}^{k-1}+r_{1,2,3}^9{P}_{m+2,n}^{k-1}\\ +r_{1,2,3}^{10}{P}_{m+3,n}^{k+1}+r_{1,2,3}^{11}{P}_{m,n}^{k-2}+r_{1,2,3}^{12}{P}_{m+3,n}^k+r_{1,2,3}^{13}{P}_{m+3,n}^{k-1}+r_{1,2,3}^{14}{P}_{m+1,n}^{k-2}+r_{1,2,3}^{15}{P}_{m+2,n}^{k-2}+r_{1,2,3}^{16}{P}_{m+3,n}^{k-2})/(-r_{1,2,3}^1)\end{array}$

其中，m,n为空间离散网格点坐标，m＝0,1,…,N-1,N；n＝0,1,…,N_z-1,N_z；N_z为中心波场的纵向网格点数；表示在k+1时刻边界上空间坐标为m,n的变量P的值，P为声波方程中的压强值；左边界表达式中各系数表达式如下：

$\begin{array}{ccc}\begin{array}{c}{r}_{1,2,3}^1=l_{1,1}{s}_{2,3}^1\\ r_{1,2,3}^2=l_{1,1}{s}_{2,3}^2+l_{1,2}{s}_{2,3}^1\\ r_{1,2,3}^3=l_{1,1}{s}_{2,3}^3+l_{1,3}{s}_{2,3}^1\\ r_{1,2,3}^4=l_{1,1}{s}_{2,3}^4+l_{1,2}{s}_{2,3}^3+l_{1,3}{s}_{2,3}^2+l_{1,4}{s}_{2,3}^1\\ r_{1,2,3}^5=l_{1,1}{s}_{2,3}^5+l_{1,2}{s}_{2,3}^2\\ r_{1,2,3}^6=l_{1,1}{s}_{2,3}^6+l_{1,3}{s}_{2,3}^3\\ r_{1,2,3}^7=l_{1,1}{s}_{2,3}^7+l_{1,2}{s}_{2,3}^4+l_{1,3}{s}_{2,3}^5+l_{1,4}{s}_{2,3}^2\\ r_{1,2,3}^8=l_{1,1}{s}_{2,3}^8+l_{1,2}{s}_{2,3}^6+l_{1,3}{s}_{2,3}^4+l_{1,4}{s}_{2,3}^3\end{array}& \begin{array}{c}{r}_{1,2,3}^9=l_{1,1}{s}_{2,3}^9+l_{1,2}{s}_{2,3}^8+l_{1,3}{s}_{2,3}^7+l_{1,4}{s}_{2,3}^4\\ r_{1,2,3}^{10}=l_{1,2}{s}_{2,3}^5\\ r_{1,2,3}^{11}=l_{1,3}{s}_{2,3}^6\\ r_{1,2,3}^{12}=l_{1,2}{s}_{2,3}^7+l_{1,4}{s}_{2,3}^5\\ r_{1,2,3}^{13}=l_{1,2}{s}_{2,3}^9+l_{1,4}{s}_{2,3}^7\\ r_{1,2,3}^{14}=l_{1,3}{s}_{2,3}^8+l_{1,4}{s}_{2,3}^6\\ r_{1,2,3}^{15}=l_{1,3}{s}_{2,3}^9+l_{1,4}{s}_{2,3}^8\\ r_{1,2,3}^{16}=l_{1,4}{s}_{2,3}^9\end{array}& \begin{array}{c}{s}_{2,3}^1=l_{2,1}{l}_{3,1}\\ s_{2,3}^2=l_{2,1}{l}_{3,2}+l_{2,2}{l}_{3,1}\\ s_{2,3}^3=l_{2,1}{l}_{3,3}+l_{2,3}{l}_{3,1}\\ s_{2,3}^4=l_{2,1}{l}_{3,4}+l_{2,2}{l}_{3,3}+l_{2,3}{l}_{3,2}+l_{2,4}{l}_{3,1}\\ s_{2,3}^5=l_{2,2}{l}_{3,2}\\ s_{2,3}^6=l_{2,3}{l}_{3,3}\\ s_{2,3}^7=l_{2,2}{l}_{3,4}+l_{2,4}{l}_{3,2}\\ s_{2,3}^8=l_{2,3}{l}_{3,4}+l_{2,4}{l}_{3,3}\\ s_{2,3}^9=l_{2,4}{l}_{3,4}\end{array}\end{array}$

其中，Δt、Δh分别为时间、空间采样间隔；右边界、上边界和下边界类比以上推导过程得出。