CN104749628A - 一种基于弥散黏滞性波动方程的吸收边界反射方法 - Google Patents

Abstract

本发明公开了一种基于弥散黏滞性波动方程的吸收边界反射的方法，为利用该方程进行复杂介质中地震波数值模拟提供了一种处理人工边界反射波的有效途径。首先从理论上阐述了一维和二维情况下该方法用于处理人工边界反射的原理；其次，在均匀介质和复杂储层介质模型中利用有限差分法对弥散黏滞性波动方程进行数值模拟，并采用所提的吸收边界反射方法进行处理，数值模拟结果验证了该方法的有效性。该方法也适用于声波方程和Stokes方程。此外，该技术方案易于实现，可操作性强。

Description

一种基于弥散黏滞性波动方程的吸收边界反射方法

技术领域

本发明属于地球物理勘探领域，涉及一种地震波数值模拟中的吸收边界处理方法，特别涉及一种基于弥散黏滞性波动方程的吸收边界反射方法。

背景技术

地震波数值模拟技术在油气勘探中有着重要的作用，它对于认识地下介质中波的传播规律、验证所研究的各种方法技术起关键性的作用。由于实际的地下介质是无限大的，而数值模拟的计算区域是有限的，这样会造成由于计算区域被截断而产生人工边界，当地震波传到该边界时会产生人工边界反射波。因此，需要在计算区域的边界处使用合理的吸收边界条件进行处理，才能模拟实际地下介质中波的传播。关于人工边界反射的问题，学者们已经发展了很多种方法，这些方法总体分为三类：第一类是预测法—利用波动方程的不同形式如单程波方程、波场外推等近似计算边界处的波场值；第二类是阻尼衰减法—利用合适的衰减函数将临近边界处的波场值进行衰减，在非均匀介质中由于波以不同的传播速度沿着各个方向传播，因此，这类方法中衰减函数比较难确定；第三类是完全匹配层法(PML—Perfect Matched Layer)—计算区域的内部利用波动方程计算，而在计算区域的外部需要加入适当的吸收层，并利用相应的带衰减函数的吸收边界控制方程进行计算，只要吸收层内衰减函数选择的合理，这种方法就能很好地吸收各种频率来自各个方向的波。每一种吸收边界反射的方法都有各自的优点和缺点，在第一类方法中，如Clayton-Engquist方法和Higdon方法的优点是需要相对较少的存储量，但是这些方法只在有限的入射角范围能够很好的吸收边界反射波；Lindman方法在较大的角度范围能达到很好的吸收效果，但它只在边界是均匀的介质中有效；Liao方法能够在较宽的角度范围和非均匀介质边界处有很好的效果，但它需要双精度实现来保证其稳定性，也需要更多的计算存储。第二类方法中，如Kosloff方法在法向入射时比较有效，但对于非法向入射会产生明显的边界反射，同时它也需要更大的存储空间。第三类PML方法是由Berenger在电磁波模拟中提出，它是目前被认为效果最好的吸收边界条件，它的实现比较简单、稳定，能够适用于具有大范围入射角和宽频带的入射波，广泛应用于很多领域，如线性Euler方程、声学介质波场模拟、弹性波的传播、黏弹性介质中波的传播和孔隙弹性介质中波的传播等。

弥散黏滞性波动方程是为了解释实际和实验中观察到的依赖于频率变化的反射现象(如：低频伴影)而提出的，这种波动方程还可用于研究在含有流体介质中地震波的衰减和频散问题。因而，弥散黏滞性波动方程的求解在地震勘探领域有着重要的意义。利用弥散黏滞性波动方程模拟地震波在含有流体介质中的传播过程也会遇到人工边界处理问题。

发明内容

本发明的目的在于针对弥散黏滞性波动方程数值模拟中的人工边界反射问题，提出一种基于弥散黏滞性波动方程的吸收边界反射方法，即在计算区域的外部加入适当的吸收层，并利用相应的带衰减函数的吸收边界控制方程进行计算，计算区域的内部利用波动方程计算。由于声波方程和Stokes方程是弥散黏滞性波动方程的特例，所以该方法也适用于声波方程和Stokes方程。

为了实现上述目的，本发明所采用的技术方案包括以下步骤：

1)根据求解一维弥散性黏滞性波动方程的吸收边界控制方程，推广确定二维弥散性黏滞性波动方程的吸收边界反射控制方程；其中，一维弥散性黏滞性波动方程为：

$\frac{{\∂}^{2}u}{{\∂t}^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}\frac{{\∂}^{3}u}{{\∂x}^2\∂t}-{\mathrm{\υ}}^2\frac{{\∂}^{2}u}{{\∂x}^2}=0---\left(1\right)$

二维弥散性黏滞性波动方程为：

$\frac{{\∂}^{2}u}{{\∂t}^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}(\frac{{\∂}^{3}u}{{\∂x}^2\∂t}+\frac{{\∂}^{3}u}{{\∂z}^2\∂t})-{\mathrm{\υ}}^2(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂z}^2})=0---\left(2\right)$

式中：u为波场函数；γ，η分别为弥散和黏滞性衰减系数，它们是岩石的孔隙度、渗透率以及流体的密度、粘度等参数的函数；υ为非频散介质中波的传播速度；x，t分别为空间和时间变量；式(1)和式(2)中左端第一项表示惯性项，第二项为弥散耗损力即扩散项，第三项表示黏滞性阻尼，第四项为波动方程的弹性部分；

二维弥散性黏滞性波动方程的吸收边界控制方程为：

$\left\{\begin{array}{c}\frac{{\∂S}_x^{*}}{\∂t}+(d_1\left(x\right)+\mathrm{\γ})S_x^{*}+\mathrm{\γ}{d}_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_x}^{*}-\mathrm{\η}\frac{{\∂R}_x^{*}}{\∂x}-{\mathrm{\υ}}^2\frac{{\∂T}_x^{*}}{\∂x}=0\\ \frac{{\∂S}_z^{*}}{\∂t}+(d_2\left(z\right)+\mathrm{\γ})S_z^{*}+\mathrm{\γ}{d}_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_2}^{*}-\mathrm{\η}\frac{{\∂R}_z^{*}}{\∂z}-{\mathrm{\υ}}^2\frac{{\∂T}_z^{*}}{\∂z}=0\\ \frac{{\∂u}_x^{*}}{\∂t}=S_x^{*}\frac{{\∂u}_z^{*}}{\∂t}=S_z^{*}\\ \frac{{\∂u}_x^{*}}{\∂x}+\frac{{\∂u}_z^{*}}{\∂x}=T_x^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_x^{*}\\ \frac{{\∂u}_x^{*}}{\∂z}+\frac{{\∂u}_z^{*}}{\∂z}=T_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_z^{*},\\ \frac{{\∂S}_x^{*}}{\∂x}+\frac{{\∂S}_z^{*}}{\∂x}=R_x^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R_x^{*}\\ \frac{{\∂S}_x^{*}}{\∂z}+\frac{{\∂S}_z^{*}}{\∂z}=R_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R_z^{*}\end{array}\right.---\left(3\right)$

式中：d₁(x)，d₂(z)分别为沿x，z方向的衰减函数；

2)建立介质模型，在其周围加入吸收边界层，吸收边界层的厚度取值至少为一个波长；

3)选择吸收边界中的吸收衰减函数，其中吸收衰减函数包括指数函数、对数函数、正/余弦函数或其中的一种或几种的复合函数；

4)对所计算的介质模型进行网格离散；

5)利用有限差分法对二维弥散黏滞性波动方程进行数值模拟，在计算区域内部利用式(4)弥散黏滞性波动方程的有限差分格式计算，在吸收边界层内利用式(3)计算，以消除人工边界反射波；其中，式(4)如下：

$\begin{array}{c}{u}_{j,m}^{n+1}=[2-4a-\mathrm{\γ}\left(\mathrm{\Δt}\right)-4b]u_{j,m}^n-[1-\mathrm{\γ}\left(\mathrm{\Δt}\right)-4a]u_{j,m}^{n-1}\\ -a(u_{j+1,m}^{n-1}+u_{j-1,m}^{n-1}+u_{j,m+1}^{n-1}+u_{j,m-1}^{n-1})+(a+b)\left[(u_{j+1,m}^n+u_{j-1,m}^n+u_{j,m+1}^n+u_{j,m-1}^n)\right]\end{array}---\left(4\right)$

式中： 表示第n时间步在网格点(x_j,z_m)处的波场值；n为时间采样点，j为空间x方向的采样点，m为空间z方向的采样点；

震源函数采用爆炸源，在空间上采用高斯函数，时间上采用Ricker子波，震源函数的形式为：

s(x,z,t)＝g(x,z)·f(t)

式中：

f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)

$g(x,z)=\frac{\exp \{-[{(x-x_0)}^2+{(z-z_0)}^2\left]\right\}}{{\mathrm{\β}}^2}$

式中：f₀表示Ricker子波的中心频率，模型计算中β为常数；(x₀,z₀)表示震源的空间位置。

所述步骤1)中，根据求解一维弥散性黏滞性波动方程的吸收边界算法，推广确定二维弥散性黏滞性波动方程的吸收边界反射控制方程的具体方法如下：

1-1)一维弥散黏滞性波动方程的求解：

对式(1)一维弥散黏滞性波动方程关于时间t做Fourier变换，整理后得到：

$\frac{{\∂}^2\hat{u}(x,\mathrm{\ω})}{{\∂x}^2}-\frac{\left(\mathrm{i\ω}\right)[\left(\mathrm{\ω}\right)+\mathrm{\γ}]}{[\left(\mathrm{i\ω}\right)\mathrm{\η}+{\mathrm{\υ}}^2]}\hat{u}(x,\mathrm{\ω})=0---\left(5\right)$

式中：是u(x,t)关于时间t的Fourier变换；ω为角频率；

式(5)是关于x的常微分方程，其特征值λ为：

${\mathrm{\λ}}^2=\frac{\left(\mathrm{i\ω}\right)[\left(\mathrm{i\ω}\right)+\mathrm{\γ}]}{[\left(\mathrm{i\ω}\right)\mathrm{\η}+{\mathrm{\υ}}^2]}---\left(6\right)$

由式(6)得到

${\mathrm{\λ}}_k=r^{\frac{1}{2}}[\cos \left(\frac{\mathrm{\θ}+2\mathrm{k\π}}{2}\right)+i\sin \left(\frac{\mathrm{\θ}+2\mathrm{k\π}}{2}\right)],k=\mathrm{1,2}---\left(7\right)$

式中：为λ_k的模；θ为相位角，由arctanθ＝I/R确定；λ_k的实部R和虚部I分别为：

$R=\frac{{-\mathrm{\ω}}^2[{\mathrm{\υ}}^2-\mathrm{\γ\η}]}{[{\mathrm{\υ}}^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2]}---\left(8\right)$

$I=\frac{\mathrm{\ω}[{\mathrm{\γ\υ}}^2+{\mathrm{\ω}}^2]}{[{\mathrm{\υ}}^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2]}---\left(9\right)$

所以式(1)解的形式为：

$u(x,t)=\frac{1}{2\mathrm{\π}}\∫\hat{u}(x,\mathrm{\ω})e^{\mathrm{i\ωt}}\mathrm{d\ω}=\frac{1}{2\mathrm{\π}}\∫[c_1{e}^{{\mathrm{\λ}}_{1}x}+c_2{e}^{{\mathrm{\λ}}_{2}x}]e^{\mathrm{i\ωt}}\mathrm{d\ω}---\left(10\right)$

式中：c₁，c₂为常数，由初始条件确定；

当衰减参数γ为零时，式(1)是即为Stokes方程，而这种情况下式(10)是即为Stokes方程的解。当衰减参数γ，η均为零时，式(10)式即为一维声波方程的解；

将式(1)写成等价的偏微分方程组，然后构造与其相应的辅助方程组并引入坐标变换，得到在新坐标系下的方程组，进而得到式(1)的解是指数衰减的，从而利用构建带衰减函数的方程组用于吸收人工边界反射波。

与式(1)等价的一阶偏微分方程组为：

$\left\{\begin{array}{c}\frac{\∂S}{\∂t}(x,t)+\mathrm{fS}(x,t)-\mathrm{\η}\frac{\∂R}{\∂x}(x,t)-{\mathrm{\υ}}^2\frac{\∂T}{\∂x}(x,t)=0\\ \frac{\∂u}{\∂t}(x,t)=S(x,t)\\ \frac{\∂u}{\∂x}(x,t)=T(x,t)\\ \frac{\∂S}{\∂x}=R(x,t)\end{array}\right.---\left(11\right)$

其中， $S=\frac{\∂u}{\∂t},T=\frac{\∂u}{\∂x},R=\frac{\∂S}{\∂x}$ 为中间变量；

构建与式(1)和式(11)相对应的含有衰减函数的吸收边界控制方程组为：

$\left\{\begin{array}{c}\frac{\∂S^{*}}{\∂t}(x,t)+(d\left(x\right)+\mathrm{\γ})S^{*}(x,t)+\mathrm{\γd}\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*S^{*}(x,t)\\ -\mathrm{\η}\frac{{\∂R}^{*}}{\∂x}(x,t)-{\mathrm{\υ}}^2\frac{{\∂T}^{*}}{\∂x}(x,t)=0\\ \frac{{\∂u}^{*}}{\∂t}(x,t)=S^{*}(x,t)\\ \frac{{\∂u}^{*}}{\∂x}(x,t)=T^{*}(x,t)+d\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T^{*}(x,t)\\ \frac{{\∂S}^{*}}{\∂x}=R^{*}(x,t)+d\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R^{*}(x,t)\end{array}\right.---\left(12\right)$

式中：为中间变量；sgn(t)为符号函数；“*”表示时间域卷积；

关于空间变量x引入坐标变换

$\tilde{x}=x-\frac{i}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}---\left(13\right)$

将式(13)代入式(12)，再关于t做Fourier变换，得到：

$\left\{\begin{array}{c}\left(\mathrm{i\ω}\right){\hat{S}}^{*}(\tilde{x},\mathrm{\ω})+\mathrm{\γ}{\hat{S}}^{*}(\tilde{x},\mathrm{\ω})-\mathrm{\η}\frac{\∂{\hat{R}}^{*}(\tilde{x},\mathrm{\ω})}{\∂\tilde{x}}-{\mathrm{\υ}}^2\frac{\∂{\hat{T}}^{*}(\tilde{x},\mathrm{\ω})}{\∂\tilde{x}}=0\\ \left(\mathrm{i\ω}\right){\hat{u}}^{*}(\tilde{x},\mathrm{\ω})={\hat{S}}^{*}(\tilde{x},\mathrm{\ω})\\ \frac{\∂{\hat{u}}^{*}}{\∂\tilde{x}}(\tilde{x},\mathrm{\ω})={\hat{T}}^{*}(\tilde{x},\mathrm{\ω})\\ \frac{\∂{\hat{S}}^{*}}{\∂\tilde{x}}(\tilde{x},\mathrm{\ω})={\hat{R}}^{*}(\tilde{x},\mathrm{\ω})\end{array}\right.---\left(14\right)$

再对(14)式关于角频率ω做逆Fourier变换，得到：

$\left\{\begin{array}{c}\frac{{\∂S}^{*}}{\∂t}(\tilde{x},t)+\mathrm{\γ}{S}^{*}(\tilde{x},t)-\mathrm{\η}\frac{\∂R^{*}}{\∂\tilde{x}}(\tilde{x},t)-{\mathrm{\υ}}^2\frac{\∂T^{*}}{\∂\tilde{x}}(\tilde{x},t)=0\\ \frac{{\∂u}^{*}}{\∂t}(\tilde{x},t)=S^{*}(\tilde{x},t)\\ \frac{{\∂u}^{*}}{\∂\tilde{x}}(\tilde{x},t)=T^{*}(\tilde{x},t)\\ \frac{{\∂S}^{*}}{\∂\tilde{x}}(\tilde{x},t)=R^{*}(\tilde{x},t)\\ \end{array}\right.---\left(15\right)$

式(15)和式(11)的形式完全相同，不同在于空间坐标与x，前者为复数，后者为实数；由式(1)的解可知，式(15)的解为：

$\left\{\begin{array}{c}{u}^{*}(\tilde{x},t)=\frac{1}{2\mathrm{\π}}\∫\hat{u}(\tilde{x},\mathrm{\ω})e^{\mathrm{i\ωt}}\mathrm{d\ω}\\ =\frac{1}{2\mathrm{\π}}{\∫e}^{\mathrm{i\ωt}}[c_3{e}^{{\mathrm{\λ}}_{1}x}{e}^{-\frac{{\mathrm{i\λ}}_1}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}}+c_4{e}^{{\mathrm{\λ}}_{2}x}{e}^{-\frac{{\mathrm{i\λ}}_2}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}}]\mathrm{d\ω}\end{array}\right.---\left(16\right)$

式中：c₃，c₄为常数，它们由初始条件确定；

从式(16)可以明显看出，方程组(15)的解沿x方向指数衰减，衰减形态和衰减大小由d(x)决定；因此，只要恰当的选择该函数使得式(13)在计算区域内部为式(1)式的解，而在边界处波函数呈现快速衰减，从而达到吸收边界反射波的目的；

1-2)二维弥散黏滞性波动方程吸收边界算法原理

与一维情况相同，引入中间变量S_x，S_z，T_x，T_z，R_x，R_z，u_x，u_z，并令 $S_x=\frac{{\∂u}_x}{\∂t},S_z=\frac{{\∂u}_z}{\∂t},R_x=\frac{{\∂S}_x}{\∂x},R_z=\frac{{\∂S}_z}{\∂z},T_x=\frac{{\∂u}_x}{\∂x},T_z=\frac{{\∂u}_z}{\∂z},$ 则与式(2)二维弥散黏滞性波动方程等价的偏微分方程组的形式为：

$\left\{\begin{array}{c}\frac{{\∂S}_x}{\∂t}+\mathrm{\γ}{S}_x-\mathrm{\η}\frac{{\∂R}_x}{\∂x}-{\mathrm{\υ}}^2\frac{{\∂T}_x}{\∂x}=0\\ \frac{{\∂S}_z}{\∂t}+\mathrm{\γ}{S}_z-\mathrm{\η}\frac{{\∂R}_z}{\∂z}-{\mathrm{\υ}}^2\frac{{\∂T}_z}{\∂z}=0\\ \frac{{\∂u}_x}{\∂t}=S_x,\frac{{\∂u}_z}{\∂t}=S_z\\ \frac{{\∂S}_x}{\∂x}+\frac{{\∂S}_z}{\∂x}=R_x,\frac{{\∂S}_x}{\∂z}+\frac{{\∂S}_z}{\∂z}=R_z\\ \frac{{\∂u}_x}{\∂x}+\frac{{\∂u}_z}{\∂x}=T_x,\frac{{\∂u}_x}{\∂z}+\frac{{\∂u}_z}{\∂z}=T_z\end{array}\right.---\left(17\right)$

同样构建与式(17)对应的带有衰减函数的辅助方程组为：

$\left\{\begin{array}{c}\frac{{\∂S}_x^{*}}{\∂t}+(d_1\left(x\right)+\mathrm{\γ})S_x^{*}+\mathrm{\γ}{d}_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_x}^{*}-\mathrm{\η}\frac{{\∂R}_x^{*}}{\∂x}-{\mathrm{\υ}}^2\frac{{\∂T}_x^{*}}{\∂x}=0\\ \frac{{\∂S}_z^{*}}{\∂t}+(d_2\left(z\right)+\mathrm{\γ})S_z^{*}+\mathrm{\γ}{d}_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_2}^{*}-\mathrm{\η}\frac{{\∂R}_z^{*}}{\∂z}-{\mathrm{\υ}}^2\frac{{\∂T}_z^{*}}{\∂z}=0\\ \frac{{\∂u}_x^{*}}{\∂t}=S_x^{*}\frac{{\∂u}_z^{*}}{\∂t}=S_z^{*}\\ \frac{{\∂u}_x^{*}}{\∂x}+\frac{{\∂u}_z^{*}}{\∂x}=T_x^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_x^{*}\\ \frac{{\∂u}_x^{*}}{\∂z}+\frac{{\∂u}_z^{*}}{\∂z}=T_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_z^{*},\\ \frac{{\∂S}_x^{*}}{\∂x}+\frac{{\∂S}_z^{*}}{\∂x}=R_x^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R_x^{*}\\ \frac{{\∂S}_x^{*}}{\∂z}+\frac{{\∂S}_z^{*}}{\∂z}=R_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R_z^{*}\end{array}\right.---\left(18\right).$

所述步骤3)中，吸收衰减函数采用正余弦函数，其形式为：

d₁(x)＝M[1-cos(2πx/l)]

d₂(z)＝M[1-cos(2πz/l)]

式中：M为幅值常量，l为吸收边界层的厚度。

与现有技术相比，本发明具有以下有益效果：

本发明给出了弥散黏滞性波动方程的吸收边界控制方程，为利用弥散黏滞性波动方程进行地震波有限差分数值模拟提供了一种有效地处理人工边界反射波的方法；弥散黏滞性波动方程可用于描述复杂含流体介质中地震波的传播特性，在复杂介质中，只要在吸收边界层内将本发明所得的吸收边界控制方程中相应的介质参数(地震波的传播速度、弥散衰减系数和黏滞性衰减系数)利用计算区域临近边界处的值进行替换，且吸收边界层内衰减函数选择的合理，所得方法也能够处理复杂介质中由于人工边界而产生的人工边界反射波。本发明应用较为广泛，对于与弥散黏滞性波动方程的形式相同的方程都能适用。特别地，由于声波方程和Stokes方程是弥散黏滞性波动方程的特例，因而，本发明所得方法也可用于声波方程和Stokes方程的有限差分数值模拟。

图及表说明

图1本发明中均匀介质模型PML计算区域示意图；

图2均匀介质模型下，当γ＝0.0，η＝0.02时的边界处理结果对比图；其中，(a)为利用本文吸收边界条件处理后的波场快照(320ms)，(b)为未做边界处理的波场快照(320ms)，(c)为均匀介质中在(350m,700m)处接受到的地震记录，实线为经过边界处理后的结果，虚线为未做边界处理的结果；

图3均匀介质模型下，当γ＝10.0，η＝0.0时的边界处理结果对比图；其中，(a)为利用本文吸收边界条件处理后的波场快照(320ms)，(b)为未做边界处理的波场快照(320ms)，(c)为均匀介质中在(350m,700m)处接受到的地震记录，实线为经过边界处理后的结果，虚线为未做边界处理的结果；

图4均匀介质模型下，当γ＝20.0，η＝0.02时的边界处理结果对比图；其中，(a)为利用本文吸收边界条件处理后的波场快照(320ms)，(b)为未做边界处理的波场快照(320ms)，(c)为均匀介质中在(350m,700m)处接受到的地震记录，实线为经过边界处理后的结果，虚线为未做边界处理的结果；

图5复杂储层介质模型示意图；

图6弥散黏滞性波在复杂储层介质模型下，采用所提边界反射方法处理后不同时刻的波场快照；

图7复杂储层介质模型中采用本发明所提吸收边界条件处理后的VSP地震记录。

具体实施方式

下面结合附图对本发明做进一步详细的说明：

本发明吸收边界反射方法的原理如下：

一维和二维弥散黏滞性波动方程的数学描述为：

$\frac{{\∂}^{2}u}{{\∂t}^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}\frac{{\∂}^{3}u}{{\∂x}^2\∂t}-{\mathrm{\υ}}^2\frac{{\∂}^{2}u}{{\∂x}^2}=0---\left(1\right)$

$\frac{{\∂}^{2}u}{{\∂t}^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}(\frac{{\∂}^{3}u}{{\∂x}^2\∂t}+\frac{{\∂}^{3}u}{{\∂z}^2\∂t})-{\mathrm{\υ}}^2(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂z}^2})=0---\left(2\right)$

式中：u为波场函数；γ，η分别为弥散和黏滞性衰减系数，它们是岩石的孔隙度、渗透率以及流体的密度、粘度等参数的函数；υ为非频散介质中波的传播速度；x，t分别为空间和时间变量。(1)和(2)式中左端第一项表示惯性项，第二项为弥散耗损力即扩散项，第三项表示黏滞性阻尼，第四项为波动方程的弹性部分。我们注意到，(1)和(2)既包含了惯性项也包含了扩散项，也就是耦合了惯性流与扩散流，或者说耦合了低频扩散与高频波现象。

首先以一维弥散黏滞性波动方程为例说明该吸收边界反射方法的原理，其次，将其推广为二维情况。

1.一维弥散黏滞性波动方程吸收边界算法原理

一维弥散黏滞性波动方程的求解：

对一维弥散黏滞性波动式(1)关于时间t做Fourier变换，整理后得到：

$\frac{{\∂}^2\hat{u}(x,\mathrm{\ω})}{{\∂x}^2}-\frac{\left(\mathrm{i\ω}\right)[\left(\mathrm{\ω}\right)+\mathrm{\γ}]}{[\left(\mathrm{i\ω}\right)\mathrm{\η}+{\mathrm{\υ}}^2]}\hat{u}(x,\mathrm{\ω})=0---\left(3\right)$

式中：是u(x,t)关于时间t的Fourier变换；ω为角频率。

(3)式是关于x的常微分方程，其特征值λ为：

${\mathrm{\λ}}^2=\frac{\left(\mathrm{i\ω}\right)[\left(\mathrm{i\ω}\right)+\mathrm{\γ}]}{[\left(\mathrm{i\ω}\right)\mathrm{\η}+{\mathrm{\υ}}^2]}---\left(4\right)$

由(4)得到

${\mathrm{\λ}}_k=r^{\frac{1}{2}}[\cos \left(\frac{\mathrm{\θ}+2\mathrm{k\π}}{2}\right)+i\sin \left(\frac{\mathrm{\θ}+2\mathrm{k\π}}{2}\right)],k=\mathrm{1,2}---\left(5\right)$

式中：为λ_k的模；θ为相位角，由arctanθ＝I/R确定；λ_k的实部R和虚部I分别为：

$R=\frac{{-\mathrm{\ω}}^2[{\mathrm{\υ}}^2-\mathrm{\γ\η}]}{[{\mathrm{\υ}}^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2]}---\left(6\right)$

$I=\frac{\mathrm{\ω}[{\mathrm{\γ\υ}}^2+{\mathrm{\ω}}^2]}{[{\mathrm{\υ}}^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2]}---\left(7\right)$

所以式(1)解的形式为：

$u(x,t)=\frac{1}{2\mathrm{\π}}\∫\hat{u}(x,\mathrm{\ω})e^{\mathrm{i\ωt}}\mathrm{d\ω}=\frac{1}{2\mathrm{\π}}\∫[c_1{e}^{{\mathrm{\λ}}_{1}x}+c_2{e}^{{\mathrm{\λ}}_{2}x}]e^{\mathrm{i\ωt}}\mathrm{d\ω}---\left(8\right)$

式中：c₁，c₂为常数，由初始条件确定。

当衰减参数γ为零时，(1)是即为Stokes方程，而这种情况下(8)是即为Stokes方程的解。当衰减参数γ，η均为零时，(8)式即为一维声波方程的解。

将式(1)写成等价的偏微分方程组，然后构造与其相应的辅助方程组并引入坐标变换，得到在新坐标系下的方程组，进而得到式(1)的解是指数衰减的，从而利用构建带衰减函数的方程组用于吸收人工边界反射波。

与式(1)等价的一阶偏微分方程组为：

$\left\{\begin{array}{c}\frac{\∂S}{\∂t}(x,t)+\mathrm{fS}(x,t)-\mathrm{\η}\frac{\∂R}{\∂x}(x,t)-{\mathrm{\υ}}^2\frac{\∂T}{\∂x}(x,t)=0\\ \frac{\∂u}{\∂t}(x,t)=S(x,t)\\ \frac{\∂u}{\∂x}(x,t)=T(x,t)\\ \frac{\∂S}{\∂x}=R(x,t)\end{array}\right.---\left(9\right)$

其中， $S=\frac{\∂u}{\∂t},T=\frac{\∂u}{\∂x},R=\frac{\∂S}{\∂x}$ 为中间变量。

构建与(1)式、(9)式相对应的含有衰减函数的辅助方程组为：

$\left\{\begin{array}{c}\frac{\∂S^{*}}{\∂t}(x,t)+(d\left(x\right)+\mathrm{\γ})S^{*}(x,t)+\mathrm{\γd}\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*S^{*}(x,t)\\ -\mathrm{\η}\frac{{\∂R}^{*}}{\∂x}(x,t)-{\mathrm{\υ}}^2\frac{{\∂T}^{*}}{\∂x}(x,t)=0\\ \frac{{\∂u}^{*}}{\∂t}(x,t)=S^{*}(x,t)\\ \frac{{\∂u}^{*}}{\∂x}(x,t)=T^{*}(x,t)+d\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T^{*}(x,t)\\ \frac{{\∂S}^{*}}{\∂x}=R^{*}(x,t)+d\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R^{*}(x,t)\end{array}\right.---\left(10\right)$

式中：为中间变量；sgn(t)为符号函数；“*”表示时间域卷积。

关于空间变量x引入坐标变换

$\tilde{x}=x-\frac{i}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}---\left(11\right)$

将(11)式代入(10)式，再关于t做Fourier变换，得到：

$\left\{\begin{array}{c}\left(\mathrm{i\ω}\right){\hat{S}}^{*}(\tilde{x},\mathrm{\ω})+\mathrm{\γ}{\hat{S}}^{*}(\tilde{x},\mathrm{\ω})-\mathrm{\η}\frac{\∂{\hat{R}}^{*}(\tilde{x},\mathrm{\ω})}{\∂\tilde{x}}-{\mathrm{\υ}}^2\frac{\∂{\hat{T}}^{*}(\tilde{x},\mathrm{\ω})}{\∂\tilde{x}}=0\\ \left(\mathrm{i\ω}\right){\hat{u}}^{*}(\tilde{x},\mathrm{\ω})={\hat{S}}^{*}(\tilde{x},\mathrm{\ω})\\ \frac{\∂{\hat{u}}^{*}}{\∂\tilde{x}}(\tilde{x},\mathrm{\ω})={\hat{T}}^{*}(\tilde{x},\mathrm{\ω})\\ \frac{\∂{\hat{S}}^{*}}{\∂\tilde{x}}(\tilde{x},\mathrm{\ω})={\hat{R}}^{*}(\tilde{x},\mathrm{\ω})\end{array}\right.---\left(12\right)$

再对(12)式关于角频率ω做逆Fourier变换，得到：

$\left\{\begin{array}{c}\frac{{\∂S}^{*}}{\∂t}(\tilde{x},t)+\mathrm{\γ}{S}^{*}(\tilde{x},t)-\mathrm{\η}\frac{\∂R^{*}}{\∂\tilde{x}}(\tilde{x},t)-{\mathrm{\υ}}^2\frac{\∂T^{*}}{\∂\tilde{x}}(\tilde{x},t)=0\\ \frac{{\∂u}^{*}}{\∂t}(\tilde{x},t)=S^{*}(\tilde{x},t)\\ \frac{{\∂u}^{*}}{\∂\tilde{x}}(\tilde{x},t)=T^{*}(\tilde{x},t)\\ \frac{{\∂S}^{*}}{\∂\tilde{x}}(\tilde{x},t)=R^{*}(\tilde{x},t)\\ \end{array}\right.---\left(13\right)$

(13)式和(9)式的形式完全相同，不同在于空间坐标与x，前者为复数，后者为实数。由式(1)式的解(8)式可知，(13)式的解为：

$\left\{\begin{array}{c}{u}^{*}(\tilde{x},t)=\frac{1}{2\mathrm{\π}}\∫\hat{u}(\tilde{x},\mathrm{\ω})e^{\mathrm{i\ωt}}\mathrm{d\ω}\\ =\frac{1}{2\mathrm{\π}}{\∫e}^{\mathrm{i\ωt}}[c_3{e}^{{\mathrm{\λ}}_{1}x}{e}^{-\frac{{\mathrm{i\λ}}_1}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}}+c_4{e}^{{\mathrm{\λ}}_{2}x}{e}^{-\frac{{\mathrm{i\λ}}_2}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}}]\mathrm{d\ω}\end{array}\right.---\left(14\right)$

式中：c₃，c₄为常数，它们由初始条件确定。

从(14)式可以明显看出，方程组(13)式的解沿x方向指数衰减，衰减形态和衰减大小由d(x)决定。因此，只要恰当的选择该函数使得(13)式在计算区域内部为式(1)式的解，而在边界处波函数呈现快速衰减，这样既可达到吸收边界反射波的目的。

2.二维弥散黏滞性波动方程吸收边界算法原理

与一维情况类似，引入中间变量S_x，S_z，T_x，T_z，R_x，R_z，u_x，u_z，并令 $S_x=\frac{{\∂u}_x}{\∂t},S_z=\frac{{\∂u}_z}{\∂t},R_x=\frac{{\∂S}_x}{\∂x},R_z=\frac{{\∂S}_z}{\∂z},T_x=\frac{{\∂u}_x}{\∂x},T_z=\frac{{\∂u}_z}{\∂z},$ 则与二维弥散黏滞性波动方程(2)式等价的偏微分方程组的形式为：

$\left\{\begin{array}{c}\frac{{\∂S}_x}{\∂t}+\mathrm{\γ}{S}_x-\mathrm{\η}\frac{{\∂R}_x}{\∂x}-{\mathrm{\υ}}^2\frac{{\∂T}_x}{\∂x}=0\\ \frac{{\∂S}_z}{\∂t}+\mathrm{\γ}{S}_z-\mathrm{\η}\frac{{\∂R}_z}{\∂z}-{\mathrm{\υ}}^2\frac{{\∂T}_z}{\∂z}=0\\ \frac{{\∂u}_x}{\∂t}=S_x,\frac{{\∂u}_z}{\∂t}=S_z\\ \frac{{\∂S}_x}{\∂x}+\frac{{\∂S}_z}{\∂x}=R_x,\frac{{\∂S}_x}{\∂z}+\frac{{\∂S}_z}{\∂z}=R_z\\ \frac{{\∂u}_x}{\∂x}+\frac{{\∂u}_z}{\∂x}=T_x,\frac{{\∂u}_x}{\∂z}+\frac{{\∂u}_z}{\∂z}=T_z\end{array}\right.---\left(15\right)$

同样构建与(15)式对应的带有衰减函数的辅助方程组为：

$\left\{\begin{array}{c}\frac{{\∂S}_x^{*}}{\∂t}+(d_1\left(x\right)+\mathrm{\γ})S_x^{*}+\mathrm{\γ}{d}_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_x}^{*}-\mathrm{\η}\frac{{\∂R}_x^{*}}{\∂x}-{\mathrm{\υ}}^2\frac{{\∂T}_x^{*}}{\∂x}=0\\ \frac{{\∂S}_z^{*}}{\∂t}+(d_2\left(z\right)+\mathrm{\γ})S_z^{*}+\mathrm{\γ}{d}_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_2}^{*}-\mathrm{\η}\frac{{\∂R}_z^{*}}{\∂z}-{\mathrm{\υ}}^2\frac{{\∂T}_z^{*}}{\∂z}=0\\ \frac{{\∂u}_x^{*}}{\∂t}=S_x^{*}\frac{{\∂u}_z^{*}}{\∂t}=S_z^{*}\\ \frac{{\∂u}_x^{*}}{\∂x}+\frac{{\∂u}_z^{*}}{\∂x}=T_x^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_x^{*}\\ \frac{{\∂u}_x^{*}}{\∂z}+\frac{{\∂u}_z^{*}}{\∂z}=T_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_z^{*},\\ \frac{{\∂S}_x^{*}}{\∂x}+\frac{{\∂S}_z^{*}}{\∂x}=R_x^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R_x^{*}\\ \frac{{\∂S}_x^{*}}{\∂z}+\frac{{\∂S}_z^{*}}{\∂z}=R_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R_z^{*}\end{array}\right.---\left(16\right)$

式中：d₁(x)，d₂(z)分别为沿x，z方向的衰减函数。

同一维情形，可知(16)式的解随传播距离指数衰减，衰减形态由函数d₁(x)，d₂(z)决定。在边界处平行于x方向选择合理的衰减函数d₂(z)，平行于z方向选择合理的衰减函数d₁(x)，而在与x和z平行方向的交接区域，沿x、z方向的衰减函数都要使用。利用这种方法计算波场使得在边界处能达到好的吸收效果，其计算区域示意图如图1所示。

本发明的基础是地震波的传播及其数值模拟方法，本发明基于弥散黏滞性波动方程提出了一种吸收边界反射的方法，具体实施步骤分别为：

1)介质模型的吸收边界

对于所建立的介质模型，在其周围加入合理的吸收边界层，吸收边界层的厚度与地震波的波长有关，一般情况下，吸收边界的厚度取值至少为一个波长。

本发明中均匀介质模型PML计算区域示意图如图1所示，“0”表示内部计算区域，“1”、“2”和“3”表示PML吸收层，在这些吸收层内所加的吸收衰减函数不同，如：“1”区域的衰减函数需沿着z方向呈衰减趋势，“2”区域的衰减函数需沿着x方向呈衰减趋势，“3”区域的衰减函数需在x，z方向都具有衰减特性。

2)衰减函数的选择

吸收边界中的吸收衰减函数可以有不同的形式，如指数函数、对数函数以及正/余弦函数等。本发明中利用正余弦函数，其形式为：

d₁(x)＝M[1-cos(2πx/l)]    (17)

d₂(z)＝M[1-cos(2πz/l)]    (18)

式中：M为幅值常量，l为吸收边界层的厚度。

3)模型离散化

对所计算的介质模型进行网格离散，本发明采用矩形网格对模型进行离散。

4)有限差分正演模拟

本发明中利用有限差分法对二维弥散黏滞性波动方程进行数值模拟，在吸收边界层内利用(16)式计算，在计算区域内部利用弥散黏滞性波动方程的有限差分格式(19)式计算。

$\begin{array}{c}{u}_{j,m}^{n+1}=[2-4a-\mathrm{\γ}\left(\mathrm{\Δt}\right)-4b]u_{j,m}^n-[1-\mathrm{\γ}\left(\mathrm{\Δt}\right)-4a]u_{j,m}^{n-1}\\ -a(u_{j+1,m}^{n-1}+u_{j-1,m}^{n-1}+u_{j,m+1}^{n-1}+u_{j,m-1}^{n-1})+(a+b)\left[(u_{j+1,m}^n+u_{j-1,m}^n+u_{j,m+1}^n+u_{j,m-1}^n)\right]\end{array}---\left(19\right)$

式中： $a=\frac{\mathrm{\η\Δt}}{h^2},b=\frac{{\mathrm{\υ}}^2{\left(\mathrm{\Δt}\right)}^2}{h^2}.$

震源函数也可采用不同的形式，针对所要模拟的波场选择合理的震源函数，如爆炸源、单力源，单力偶源，双力偶源等。本发明采用爆炸源，在空间上采用高斯函数，时间上采用Ricker子波，震源函数的形式为：

s(x,z,t)＝g(x,z)·f(t)    (20)式中：

f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)

$g(x,z)=\frac{\exp \{-[{(x-x_0)}^2+{(z-z_0)}^2\left]\right\}}{{\mathrm{\β}}^2}---\left(21\right)$

式中：f₀表示Ricker子波的中心频率，模型计算中β为常数；(x₀,z₀)表示震源的空间位置。

数值模拟结果分析

本发明利用有限差分法对均匀介质模型进行模拟，并将利用本发明中吸收边界条件处理后的波场与未加边界处理的波场进行对比，验证了该吸收边界方法的有效性。在此基础上，将所发明的方法应用于复杂介质模型。

1.均匀介质模型

均匀介质模型的大小为[0，1400m]×[0，1400m]，速度为υ＝2200m/s，震源位于模型中心即空间坐标为(700m，700m)处，Ricker子波的中心频率f₀＝30Hz，时间采样间隔为dt＝0.001s，时间延迟50ms，空间采样间隔为dx＝dz＝3.5m，最大计算时间t_max＝800ms，吸收边界层的厚度l＝20。该模型中分三种情况进行数值分析：(1)γ＝0.0，η＝0.02，(2)γ＝10.0，η＝0.0，(3)γ＝20.0，η＝0.02。

图2—图4(a)、(b)分别给出了弥散和黏滞性衰减参数取值不同时利用所提边界条件处理方法处理后和未做边界处理的波场快照；图2—图4(c)分别给出了在(350m，700m)处这三种情形下的地震记录。

如图2—图4所示，均匀介质模型中在不同衰减参数下用发明所提的吸收边界方法处理后和未做边界处理的数值模拟结果对比可以明显看出，采用边界处理方法后，来自边界的反射波能被很好地吸收。

2.复杂储层介质模型

复杂储层介质模型如图5所示，模型大小为[0，900m][0，900m]，两个均匀泥岩介质层中间为50m厚的砂岩储层，层位于深度412.5m处，在砂岩储层中间为10m厚的优质储层，且优质储层均匀地分为三段，这三段储层的声波波阻抗相同，但衰减参数不同。介质参数如表1所示。震源位于(450m，9m)处，Ricker子波中心频率为f₀＝30Hz，空间采样间隔为dx＝dz＝1.5m，时间采样间隔为dt＝0.00018s，时间延迟50.4ms。

表1 复杂储层介质模型参数

图6为在该储层介质模型下，利用有限差分方法对弥散黏滞性波动方程进行数值模拟得到的不同时刻的波场快照。图7为VSP记录。

从图6和图7很明显地看到：采用本发明提出的吸收边界处理方法可以有效地吸收复杂储层介质模型中的边界反射波，说明了该方法的有效性。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (3)

1.一种基于弥散黏滞性波动方程的吸收边界反射方法，其特征在于，包括以下步骤：

1)根据求解一维弥散性黏滞性波动方程的吸收边界控制方程，推广确定二维弥散性黏滞性波动方程的吸收边界反射控制方程；其中，一维弥散性黏滞性波动方程为：

$\frac{{\∂}^{2}u}{\∂t^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}\frac{{\∂}^{3}u}{\∂x^2\∂t}-{\mathrm{\υ}}^2\frac{{\∂}^{2}u}{\∂x^2}=0---\left(1\right)$

二维弥散性黏滞性波动方程为：

$\frac{{\∂}^{2}u}{\∂t^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}(\frac{{\∂}^{3}u}{\∂x^2\∂t}+\frac{{\∂}^{3}u}{\∂z^2\∂t})-{\mathrm{\υ}}^2(\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂z^2})=0---\left(2\right)$

式中：u为波场函数；γ，η分别为弥散和黏滞性衰减系数，它们是岩石的孔隙度、渗透率以及流体的密度、粘度等参数的函数；υ为非频散介质中波的传播速度；x，t分别为空间和时间变量；式(1)和式(2)中左端第一项表示惯性项，第二项为弥散耗损力即扩散项，第三项表示黏滞性阻尼，第四项为波动方程的弹性部分；

二维弥散性黏滞性波动方程的吸收边界控制方程为：

$\left\{\begin{array}{c}\frac{\∂S_x^{*}}{\∂t}+(d_1\left(x\right)+\mathrm{\γ})S_x^{*}+\mathrm{\γ}{d}_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_x}^{*}-\mathrm{\η}\frac{\∂R_x^{*}}{\∂x}-{\mathrm{\υ}}^2\frac{\∂T_x^{*}}{\∂x}=0\\ \frac{\∂S_z^{*}}{\∂t}+(d_2\left(z\right)+\mathrm{\γ})S_z^{*}+\mathrm{\γ}{d}_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_z}^{*}-\mathrm{\η}\frac{\∂R_z^{*}}{\∂z}-{\mathrm{\υ}}^2\frac{\∂T_z^{*}}{\∂z}=0\\ \frac{\∂u_x^{*}}{\∂t}=S_x^{*},\frac{\∂u_z^{*}}{\∂t}=S_z^{*}\\ \frac{\∂u_x^{*}}{\∂x}+\frac{\∂u_z^{*}}{\∂x}=T_x^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_x^{*}\\ \frac{\∂u_x^{*}}{\∂z}+\frac{\∂u_z^{*}}{\∂z}=T_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_z^{*}\\ \frac{\∂S_x^{*}}{\∂x}+\frac{\∂S_z^{*}}{\∂x}=R_s^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R_x^{*}\\ \frac{\∂S_x^{*}}{\∂z}+\frac{\∂S_z^{*}}{\∂z}=R_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_z^{*}\end{array}\right.---\left(3\right)$

式中：d₁(x)，d₂(z)分别为沿x，z方向的衰减函数；

2)建立介质模型，在其周围加入吸收边界层，吸收边界层的厚度取值至少为一个波长；

3)选择吸收边界中的吸收衰减函数，其中吸收衰减函数包括指数函数、对数函数、正/余弦函数或其中的一种或几种的复合函数；

4)对所计算的介质模型进行网格离散；

5)利用有限差分法对二维弥散黏滞性波动方程进行数值模拟，在计算区域内部利用式(4)弥散黏滞性波动方程的有限差分格式计算，在吸收边界层内利用式(3)计算，以消除人工边界反射波；其中，式(4)如下：

$\begin{array}{c}{u}_{j,m}^{n+1}=[2-4a-\mathrm{\γ}\left(\mathrm{\Δt}\right)-4b]u_{j,m}^n-[1-\mathrm{\γ}\left(\mathrm{\Δt}\right)-4a]u_{j,m}^{n-1}\\ -a(u_{j+1,m}^{n-1}+u_{j-1,m}^{n-1}+u_{j,m+1}^{n-1}+u_{j,m-1}^{n-1})+(a+b)\left[(u_{j+1,m}^n+u_{j-1,m}^n+u_{j,m+1}^n+u_{j,m-1}^n)\right]\end{array}---\left(4\right)$

式中： 表示第n时间步在网格点(x_j,z_m)处的波场值；n为时间采样点，j为空间x方向的采样点，m为空间z方向的采样点；

震源函数采用爆炸源，在空间上采用高斯函数，时间上采用Ricker子波，震源函数的形式为：

s(x,z,t)＝g(x,z)·f(t)

式中：

f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)

$g(x,z)=\frac{\exp [-[{(x-x_0)}^2+{(z-z_0)}^2\left]\right\}}{{\mathrm{\β}}^2}$

式中：f₀表示Ricker子波的中心频率，模型计算中β为常数；(x₀,z₀)表示震源的空间位置。

2.根据权利要求1所述的基于弥散黏滞性波动方程的吸收边界反射方法，其特征在于：所述步骤1)中，根据求解一维弥散性黏滞性波动方程的吸收边界算法，推广确定二维弥散性黏滞性波动方程的吸收边界反射控制方程的具体方法如下：

1-1)一维弥散黏滞性波动方程的求解：

对式(1)一维弥散黏滞性波动方程关于时间t做Fourier变换，整理后得到：

$\frac{{\∂}^2\hat{u}(x,\mathrm{\ω})}{\∂x^2}-\frac{\left(\mathrm{i\ω}\right)[\left(\mathrm{i\ω}\right)+\mathrm{\γ}]}{[\left(\mathrm{i\ω}\right)\mathrm{\η}+{\mathrm{\υ}}^2]}\hat{u}(x,\mathrm{\ω})=0---\left(5\right)$

式中：是u(x,t)关于时间t的Fourier变换；ω为角频率；

式(5)是关于x的常微分方程，其特征值λ为：

${\mathrm{\λ}}^2=\frac{\left(\mathrm{i\ω}\right)[\left(\mathrm{i\ω}\right)+\mathrm{\γ}]}{[\left(\mathrm{i\ω}\right)\mathrm{\η}+{\mathrm{\υ}}^2]}---\left(6\right)$

由式(6)得到

${\mathrm{\λ}}_k=r^{\frac{1}{2}}[\cos \left(\frac{\mathrm{\θ}+2\mathrm{k\π}}{2}\right)+i\sin \left(\frac{\mathrm{\θ}+2\mathrm{k\π}}{2}\right)],k=\mathrm{1,2}---\left(7\right)$

式中：为λ_k的模；θ为相位角，由arctanθ＝I/R确定；λ_k的实部R和虚部I分别为：

$R=\frac{-{\mathrm{\ω}}^2[{\mathrm{\υ}}^2-\mathrm{\γ\η}]}{[{\mathrm{\υ}}^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2]}---\left(8\right)$

$I=\frac{\mathrm{\ω}[\mathrm{\γ}{\mathrm{\υ}}^2+{\mathrm{\ω}}^2\mathrm{\η}]}{[{\mathrm{\υ}}^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2]}---\left(9\right)$

所以式(1)解的形式为：

$u(x,t)=\frac{1}{2\mathrm{\π}}\∫\hat{u}(x,\mathrm{\ω})e^{\mathrm{i\ωt}}\mathrm{d\ω}=\frac{1}{2\mathrm{\π}}\∫[c_1{e}^{{\mathrm{\λ}}_{1}x}+c_2{e}^{{\mathrm{\λ}}_{2}x}]e^{\mathrm{i\ωt}}\mathrm{d\ω}---\left(10\right)$

式中：c₁，c₂为常数，由初始条件确定；

当衰减参数γ为零时，式(1)是即为Stokes方程，而这种情况下式(10)是即为Stokes方程的解；当衰减参数γ，η均为零时，式(10)式即为一维声波方程的解；

将式(1)写成等价的偏微分方程组，然后构造与其相应的辅助方程组并引入坐标变换，得到在新坐标系下的方程组，进而得到式(1)的解是指数衰减的，从而利用构建带衰减函数的方程组用于吸收人工边界反射波；

与式(1)等价的一阶偏微分方程组为：

$\left\{\begin{array}{c}\frac{\∂S}{\∂t}(x,t)+\mathrm{fS}(x,t)-\mathrm{\η}\frac{\∂R}{\∂x}(x,y)-{\mathrm{\υ}}^2\frac{\∂T}{\∂x}(x,t)=0\\ \frac{\∂u}{\∂t}(x,t)=S(x,t)\\ \frac{\∂u}{\∂x}(x,t)=T(x,t)\\ \frac{\∂S}{\∂x}=R(x,t)\end{array}\right.---\left(11\right)$

其中， $S=\frac{\∂u}{\∂t},T=\frac{\∂u}{\∂x},R=\frac{\∂S}{\∂x}$ 为中间变量；

构建与式(1)和式(11)相对应的含有衰减函数的吸收边界控制方程组为：

$\left\{\begin{array}{c}\frac{\∂S^{*}}{\∂t}(x,t)+(d\left(x\right)+\mathrm{\γ})S^{*}(x,t)+\mathrm{\γd}\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*S^{*}(x,t)\\ -\mathrm{\η}\frac{\∂R^{*}}{\∂x}(x,t)-{\mathrm{\υ}}^2\frac{\∂T^{*}}{\∂x}(x,t)=0\\ \frac{\∂u^{*}}{\∂t}(x,t)=S^{*}(x,t)\\ \frac{\∂u^{*}}{\∂x}(x,t)=T^{*}(x,t)+d\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T^{*}(x,t)\\ \frac{\∂S^{*}}{\∂x}=R^{*}(x,t)+d\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R^{*}(x,t)\end{array}\right.---\left(12\right)$

式中：为中间变量；sgn(t)为符号函数；“*”表示时间域卷积；

关于空间变量x引入坐标变换

$\tilde{x}=x-\frac{i}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}---\left(13\right)$

将式(13)代入式(12)，再关于t做Fourier变换，得到：

$\left\{\begin{array}{c}\left(\mathrm{i\ω}\right){\hat{S}}^{*}(\tilde{x},\mathrm{\ω})+\mathrm{\γ}{\hat{S}}^{*}(\tilde{x},\mathrm{\ω})-\mathrm{\η}\frac{\∂{\hat{R}}^{*}(\tilde{x},\mathrm{\ω})}{\∂\tilde{x}}-{\mathrm{\υ}}^2\frac{\∂{\hat{T}}^{*}(\tilde{x},\mathrm{\ω})}{\∂\tilde{x}}=0\\ \left(\mathrm{i\ω}\right){\hat{u}}^{*}(\tilde{x},\mathrm{\ω})={\hat{S}}^{*}(\tilde{x},\mathrm{\ω})\\ \frac{\∂{\hat{u}}^{*}}{\∂\tilde{x}}(\tilde{x},\mathrm{\ω})={\hat{T}}^{*}(\tilde{x},\mathrm{\ω})\\ \frac{\∂{\hat{S}}^{*}}{\∂\tilde{x}}(\tilde{x},\mathrm{\ω})={\hat{R}}^{*}(\tilde{x},\mathrm{\ω})\end{array}\right.---\left(14\right)$

再对(14)式关于角频率ω做逆Fourier变换，得到：

$\left\{\begin{array}{c}\frac{\∂S^{*}}{\∂t}(\tilde{x},t)+\mathrm{\γ}{S}^{*}(\tilde{x},t)-\mathrm{\η}\frac{\∂R^{*}(\tilde{x},t)}{\∂\tilde{x}}-{\mathrm{\υ}}^2\frac{\∂T^{*}}{\∂\tilde{x}}(\tilde{x},t)=0\\ \frac{\∂u^{*}}{\∂t}(\tilde{x},t)=S^{*}(\tilde{x},t)\\ \frac{\∂{\hat{u}}^{*}}{\∂\tilde{x}}(\tilde{x},t)={\hat{T}}^{*}(\tilde{x},t)\\ \frac{\∂{\hat{S}}^{*}}{\∂\tilde{x}}(\tilde{x},t)={\hat{R}}^{*}(\tilde{x},t)\end{array}\right.---\left(15\right)$

式(15)和式(11)的形式完全相同，不同在于空间坐标与x，前者为复数，后者为实数；由式(1)的解可知，式(15)的解为：

$\begin{array}{c}{u}^{*}(\tilde{x},t)=\frac{1}{2\mathrm{\π}}\∫\hat{u}(\tilde{x},\mathrm{\ω})e^{\mathrm{i\ωt}}\mathrm{d\ω}\\ =\frac{1}{2\mathrm{\π}}\∫e^{\mathrm{i\ωt}}[c_3{e}^{{\mathrm{\λ}}_{1}x}{e}^{-\frac{i{\mathrm{\λ}}_1}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}}+c_4{e}^{{\mathrm{\λ}}_{2}x}{e}^{-\frac{i{\mathrm{\λ}}_2}{\mathrm{\ω}}\underset{0}{\overset{x}{\∫}}d\left(s\right)\mathrm{ds}}]\mathrm{d\ω}\end{array}---\left(16\right)$

式中：c₃，c₄为常数，它们由初始条件确定；

从式(16)可以明显看出，方程组(15)的解沿x方向指数衰减，衰减形态和衰减大小由d(x)决定；因此，只要恰当的选择该函数使得式(13)在计算区域内部为式(1)式的解，而在边界处波函数呈现快速衰减，从而达到吸收边界反射波的目的；

1-2)二维弥散黏滞性波动方程吸收边界算法原理

与一维情况相同，引入中间变量S_x，S_z，T_x，T_z，R_x，R_z，u_x，u_z，并令 $S_x=\frac{\∂u_x}{\∂t},S_z=\frac{\∂u_z}{\∂t},R_x=\frac{\∂S_x}{\∂x},R_z=\frac{\∂S_z}{\∂z},T_x=\frac{\∂u_x}{\∂x},T_z=\frac{\∂u_z}{\∂z},$ 则与式(2)二维弥散黏滞性波动方程等价的偏微分方程组的形式为：

$\left\{\begin{array}{c}\frac{\∂S_x}{\∂t}+\mathrm{\γ}{S}_x-\mathrm{\η}\frac{\∂R_x}{\∂x}-{\mathrm{\υ}}^2\frac{\∂T_x}{\∂x}=0\\ \frac{\∂S_z}{\∂t}+\mathrm{\γ}{S}_z-\mathrm{\η}\frac{\∂R_z}{\∂z}-{\mathrm{\υ}}^2\frac{\∂T_z}{\∂}=0\\ \frac{\∂u_x}{\∂t}=S_x,\frac{\∂u_z}{\∂t}=S_z\\ \frac{\∂S_x}{\∂x}+\frac{\∂S_z}{\∂x}=R_x,\frac{\∂S_x}{\∂z}+\frac{\∂S_z}{\∂z}=R_z\\ \frac{\∂u_x}{\∂x}+\frac{\∂u_z}{\∂x}=T_x,\frac{\∂u_x}{\∂z}+\frac{\∂u_z}{\∂z}=T_z\end{array}\right.---\left(17\right)$

同样构建与式(17)对应的带有衰减函数的辅助方程组为：

$\left\{\begin{array}{c}\frac{\∂S_x^{*}}{\∂t}+(d_1\left(x\right)+\mathrm{\γ})S_x^{*}+\mathrm{\γ}{d}_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_x}^{*}-\mathrm{\η}\frac{\∂R_x^{*}}{\∂x}-{\mathrm{\υ}}^2\frac{\∂T_x^{*}}{\∂x}=0\\ \frac{\∂S_z^{*}}{\∂t}+(d_2\left(z\right)+\mathrm{\γ})S_z^{*}+\mathrm{\γ}{d}_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*{S_z}^{*}-\mathrm{\η}\frac{\∂R_z^{*}}{\∂z}-{\mathrm{\υ}}^2\frac{\∂T_z^{*}}{\∂z}=0\\ \frac{\∂u_x^{*}}{\∂t}=S_x^{*},\frac{\∂u_z^{*}}{\∂t}=S_z^{*}\\ \frac{\∂u_x^{*}}{\∂x}+\frac{\∂u_z^{*}}{\∂x}=T_x^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_x^{*}\\ \frac{\∂u_x^{*}}{\∂z}+\frac{\∂u_z^{*}}{\∂z}=T_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_z^{*}\\ \frac{\∂S_x^{*}}{\∂x}+\frac{\∂S_z^{*}}{\∂x}=R_s^{*}+d_1\left(x\right)\frac{\mathrm{sgn}\left(t\right)}{2}*R_x^{*}\\ \frac{\∂S_x^{*}}{\∂z}+\frac{\∂S_z^{*}}{\∂z}=R_z^{*}+d_2\left(z\right)\frac{\mathrm{sgn}\left(t\right)}{2}*T_z^{*}\end{array}\right.---\left(18\right).$

3.根据权利要求1所述的基于弥散黏滞性波动方程的吸收边界反射方法，其特征在于：所述步骤3)中，吸收衰减函数采用正余弦函数，其形式为：

d₁(x)＝M[1-cos(2πx/l)]

d₂(z)＝M[1-cos(2πz/l)]

式中：M为幅值常量，l为吸收边界层的厚度。