CN106094038A - 适用于tti介质的频率域有限元全吸收pml方法 - Google Patents

Abstract

本发明涉及一种适用于TTI介质的频率域有限元全吸收PML方法，构造了一种频域有限元全吸收PML边界条件，实现了TTI介质复杂井眼条件的高精度声波测井模拟，得到了与理论解吻合的数值模拟结果，为声波测井仪器制造及利用声波测井资料评价复杂油气藏地层提供坚实的数据支撑；本发明综合了C‑PML当中的X、Z向吸收边界的拉伸因子和M‑PML当中的X、Z向吸收边界的拉伸因子分别融合为完全吸收的X向、Z向吸收边界的倏逝波抑制因子和切向衰减因子。本发明的有益效果是：通过建立全吸收PML的边界条件，兼具C‑PML和M‑PML的优势，实现了频域有限元PML吸收边界在求解TTI介质时的高效性。能高效、精确的进行TTI页岩地层的声波测井数值模拟，为声波测井仪器制造及页岩气藏开发提供必不可少的理论依据。

Description

适用于TTI介质的频率域有限元全吸收PML方法

技术领域

本发明属于地球物理(测井)勘探方法，尤其涉及弹性波数值模拟领域。

背景技术

随着油气田勘探开发的进行，页岩气等非常规油气资源成为勘探开发的热点。利用数值模拟手段考查页岩地层中井眼声波传播规律对有效勘探开发页岩油气及研究仪器构造指导仪器开发制造，具有重要意义。在各种常用数值模拟方法中，有限元方法具有其独特的优势：精度高、适应复杂边界、有利于不同场的耦合以及易于处理各种类型声源，在弹性波波场模拟中具有重要的应用价值。页岩地层可以简化为VTI(垂向横向各向同性)介质。直接对VTI地层进行模拟比较复杂，通过Bond变换旋转参考坐标系，使得VTI介质转化为TTI(水平横向各向同性)介质。频域有限元方法是TTI介质中有效的数值模拟方法。

井孔声场数值模拟所模拟的地层是无限大的，受计算机容量等的限制，数值模拟只能在有限区域内进行。为了能在有限区域中模拟无限区域中的波动过程，需要在有限区域的边界处引入吸收边界。现有比较先进的两种有限元吸收边界为C-PML和M-PML吸收边界，C-PML在极端入射角入射的情况下精度较高，但是在某些极端各向异性介质中会出现数值不稳定。M-PML在极端各向异性介质中具有较高的稳定性，但需要设置较厚的PML吸收层，大大增大了计算量，且修正因子不够优化的情况下，准确度不够高。

发明内容：

本发明的目的是提出适用于TTI介质的频率域有限元全吸收PML方法的技术方案。在频域下进行,具有更好的稳定性、更高的计算效率、更好的吸收效果、声源易于处理等优点；采用该吸收边界技术可以在吸收边界很薄的情况下有较高的精度和更好的稳定性，从而提高计算效率，此外在某些极端各向异性介质中，常规PML方法无法使用，本方法依然有效。

为了实现上述目的，本发明的技术方案是：适用于TTI介质的频率域有限元全吸收PML方法，用于构造频域有限元全吸收PML边界条件，实现TTI介质复杂井眼条件的高精度声波测井模拟；

所述方法采用复拉伸坐标系的弹性波方程，所述复拉伸坐标系的弹性波方程在所述频域下表示形式为：

$\left\{\begin{array}{c}-\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{u}=\frac{\∂{\widetilde{\mathrm{\τ}}}_{xx}}{\∂x}+\frac{\∂{\widetilde{\mathrm{\τ}}}_{xz}}{\∂z}\\ -\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{w}=\frac{\∂{\widetilde{\mathrm{\τ}}}_{zx}}{\∂x}+\frac{\∂{\widetilde{\mathrm{\τ}}}_{zz}}{\∂z}\end{array}\right.$

在PML区域内，加入拉伸函数，将上述方程表示为：

$\left\{\begin{array}{c}-\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{u}=\frac{1}{s_x}\frac{\∂{\widetilde{\mathrm{\τ}}}_{xx}}{\∂x}+\frac{1}{s_z}\frac{\∂{\widetilde{\mathrm{\τ}}}_{xz}}{\∂z}\\ -\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{w}=\frac{1}{s_x}\frac{\∂{\widetilde{\mathrm{\τ}}}_{zx}}{\∂x}+\frac{1}{s_z}\frac{\∂{\widetilde{\mathrm{\τ}}}_{zz}}{\∂z}\end{array}\right.$

所述PML是完全匹配层，s_x和s_z为PML吸收边界内的拉伸函数，s_x和s_z分别构造x和z方向的两个拉伸函数分量，s_x和s_z的表达式如下：

$s_x=k_x+\frac{d_x+{\mathrm{md}}_z}{{\mathrm{\α}}_x+i\mathrm{\ω}}$

$s_z=k_z+\frac{d_z+{\mathrm{md}}_x}{{\mathrm{\α}}_z+i\mathrm{\ω}}$

上述式中：

为x方向的实数衰减因子，

为x方向复频移因子，

为x方向倏逝波抑制因子，

L_x为吸收层x方向的厚度；

m为切向衰减因子；

为z方向的实数衰减因子，

为z方向复频移因子，

为z方向倏逝波抑制因子，

L_z为吸收层z方向的厚度，

ω为角频率，为频率的2π倍，

i为虚数单位。

更进一步，所述实数衰减因子复频移因子α₀＝2ω，倏逝波抑制因子中k₀取值范围为1～20，切向衰减因子m取值范围为0.005～0.02。

本发明的主要优点是：频率域有限元的全吸收PML，兼具有普通PML，C-PML和M-PML的优点，实现了TTI介质复杂井眼条件的高精度声波测井模拟，获得了在TTI介质中进行声波测井有限元模拟时的高效性、稳定性，可以设置非常薄的吸收边界，节省计算时间和计算内存需求，在极端各向异性介质中也可保持良好的稳定性；可以准确的实现TTI介质声波测井数值模拟，获得复杂介质大斜度井中的高精度模拟结果；得到了与理论解吻合的数值模拟结果，为声波测井仪器制造及利用声波测井资料评价复杂油气藏地层提供坚实的数据支撑。

下面结合附图和实施实例对本发明作详细描述。

附图说明

图1是本发明一个各向同性均匀模型分别采用从左到右M-PML、C-PML和全吸收PML的波场传播效果图；

图2是本发明一个极端各向异性TTI介质中m＝0.000和m＝0.02时的波场快照图，左侧为m＝0.000，右侧为m＝0.02时的波场快照图；

图3是均匀各向同性模型的参考解的波形图；

图4是本发明当m＝0.000时的波形图；

图5是本发明当m＝0.02(即采用全吸收PML时)，吸收效果较好时的波形图；

具体实施方式：

下面对本发明做进一步的详细说明。

频域有限元法具有更稳定、计算效率更高、PML效果更好以及更容易处理各种类型的声源等优点，在地球物理正演模拟尤其是在页岩气地层声波测井模拟中具有重要的应用价值。频域有限元法的基本思路是通过有限元方法求解频域下的弹性波方程，获得频谱并通过傅里叶反变换获得时域下的波形为了同时满足稳定性和准确性的要求，本发明构造了一种全吸收PML吸收边界，采用最优化的参数，使得新构造的PML兼具C-PML和M-PML两种吸收边界的优势。

本发明主要在二维直角坐标系下完成。在该坐标系下弹性波在弹性介质中的传播满足弹性波方程，在时间域的表达式为：

其中ρ为弹性介质的密度，为位移向量，C为广义虎克矩阵，在二维直角坐标系下的表达式为：

$C=\left[\begin{array}{ccc}{C}_{11}& C_{13}& C_{15}\\ C_{13}& C_{33}& C_{35}\\ C_{15}& C_{35}& C_{55}\end{array}\right]$

弹性波在频域下的形式满足弹性波方程

在二维直角坐标系下，上式化为

$\left\{\begin{array}{c}-\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{u}=\frac{\∂{\widetilde{\mathrm{\τ}}}_{xx}}{\∂\tilde{x}}+\frac{\∂{\widetilde{\mathrm{\τ}}}_{xz}}{\∂\tilde{z}}\\ -\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{w}=\frac{\∂{\widetilde{\mathrm{\τ}}}_{zx}}{\∂\tilde{x}}+\frac{\∂{\widetilde{\mathrm{\τ}}}_{zz}}{\∂\tilde{z}}\end{array}\right.$

其中各应力变量的表达式为：

${\widetilde{\mathrm{\τ}}}_{xx}=C_{11}\frac{\∂\tilde{u}}{\∂\tilde{x}}+C_{13}\frac{\∂\tilde{w}}{\∂\tilde{z}}+C_{15}(\frac{\∂\tilde{w}}{\∂\tilde{x}}+\frac{\∂\tilde{u}}{\∂\tilde{z}})$

${\widetilde{\mathrm{\τ}}}_{zz}=C_{13}\frac{\∂\tilde{u}}{\∂\tilde{x}}+C_{33}\frac{\∂\tilde{w}}{\∂\tilde{z}}+C_{35}(\frac{\∂\tilde{w}}{\∂\tilde{x}}+\frac{\∂\tilde{u}}{\∂\tilde{z}})$

${\widetilde{\mathrm{\τ}}}_{xz}={\widetilde{\mathrm{\τ}}}_{xz}=C_{15}\frac{\∂\tilde{u}}{\∂\tilde{x}}+C_{35}\frac{\∂\tilde{w}}{\∂\tilde{z}}+C_{55}(\frac{\∂\tilde{w}}{\∂\tilde{x}}+\frac{\∂\tilde{u}}{\∂\tilde{z}})$

引入复拉伸坐标系

$\frac{\∂}{\∂\tilde{n}}=\frac{1}{s_n}\frac{\∂}{\∂n}(n=x,z)$

可将弹性波方程在频域下表示为如下形式：

$\left\{\begin{array}{c}-\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{u}=\frac{1}{s_x}\frac{\∂{\widetilde{\mathrm{\τ}}}_{xx}}{\∂x}+\frac{1}{s_z}\frac{\∂{\widetilde{\mathrm{\τ}}}_{xz}}{\∂z}\\ -\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{w}=\frac{1}{s_x}\frac{\∂{\widetilde{\mathrm{\τ}}}_{zx}}{\∂x}+\frac{1}{s_z}\frac{\∂{\widetilde{\mathrm{\τ}}}_{zz}}{\∂z}\end{array}\right.$

进行坐标系复拉伸后的应力变量为

${\widetilde{\mathrm{\τ}}}_{xx}=C_{11}\frac{1}{s_x}\frac{\∂\tilde{u}}{\∂x}+C_{13}\frac{1}{s_z}\frac{\∂\tilde{w}}{\∂z}+C_{15}(\frac{1}{s_x}\frac{\∂\tilde{w}}{\∂x}+\frac{1}{s_z}\frac{\∂\tilde{u}}{\∂z})$

${\widetilde{\mathrm{\τ}}}_{zz}=C_{13}\frac{1}{s_x}\frac{\∂\tilde{u}}{\∂x}+C_{33}\frac{1}{s_{xz}}\frac{\∂\tilde{w}}{\∂z}+C_{35}(\frac{1}{s_x}\frac{\∂\tilde{w}}{\∂x}+\frac{1}{s_z}\frac{\∂\tilde{u}}{\∂z})$

${\widetilde{\mathrm{\τ}}}_{xz}={\widetilde{\mathrm{\τ}}}_{zx}=C_{15}\frac{1}{s_x}\frac{\∂\tilde{u}}{\∂x}+C_{35}\frac{1}{s_z}\frac{\∂\tilde{w}}{\∂z}+C_{55}(\frac{1}{s_x}\frac{\∂\tilde{w}}{\∂x}+\frac{1}{s_z}\frac{\∂\tilde{u}}{\∂z})$

上式即为带有PML的弹性波方程在频域下的表达式。在非PML区域，s_x和s_z的取值为1；在PML区域，s_x和s_z即为本发明构造的PML完全匹配层的x方向和z方向的拉伸函数。s_x和s_z的表达式如下：

$s_x=k_x+\frac{d_x+{\mathrm{md}}_z}{{\mathrm{\α}}_x+i\mathrm{\ω}}$

$s_z=k_z+\frac{d_z+{\mathrm{md}}_x}{{\mathrm{\α}}_z+i\mathrm{\ω}}$

上述式中：

为x方向的实数衰减因子，

为x方向复频移因子，

为x方向倏逝波抑制因子，

L_x为吸收层x方向的厚度；

m为切向衰减因子；

为z方向的实数衰减因子，

为z方向复频移因子，

为z方向倏逝波抑制因子，

L_z为吸收层z方向的厚度，

ω为角频率，为频率的2π倍，

i为虚数单位。

其中，实数衰减因子复频移因子α₀＝2ω，倏逝波抑制因子中k₀取值范围为1～20，切向衰减因子m取值范围为0.005～0.02。

为了便于推导弹性波方程在频域下的等效积分弱形式，将各等式的左右两边均乘以变量s_xs_z，可得

$\left\{\begin{array}{c}-\mathrm{\ρ}{\mathrm{\ω}}^2{s}_x{s}_z\tilde{u}=\frac{\∂\left(s_z{\widetilde{\mathrm{\τ}}}_{xx}\right)}{\∂x}+\frac{\∂\left(s_x{\widetilde{\mathrm{\τ}}}_{xz}\right)}{\∂z}\\ -\mathrm{\ρ}{\mathrm{\ω}}^2{s}_x{s}_z\tilde{w}=\frac{\∂\left(s_z{\widetilde{\mathrm{\τ}}}_{xz}\right)}{\∂x}+\frac{\∂\left(s_x{\widetilde{\mathrm{\τ}}}_{zz}\right)}{\∂z}\end{array}\right.$

将上式与势函数在弹性介质区域Ω_E作内积，并对公式两边代入格林公式，整理可得

其中L_E代表弹性介质区域的边界，代表弹性介质边界的外法向方向。

下面求解流固耦合问题。声波方程在频域下的等效积分弱形式

声波场和位移场在流固边界处在时间域满足如下公式：

${\mathrm{\ρ}}_f\frac{{\∂}^{2}u}{\∂t^2}=-\frac{\∂p}{\∂x},{\mathrm{\ρ}}_f\frac{{\∂}^{2}w}{\∂t^2}=-\frac{\∂p}{\∂z}$

其中ρ_f为流体介质的密度，在频率域的表达式为：

${\mathrm{\ρ}}_f{\mathrm{\ω}}^2\tilde{u}=\frac{\∂\tilde{p}}{\∂x},{\mathrm{\ρ}}_f{\mathrm{\ω}}^2\tilde{w}=\frac{\∂\tilde{p}}{\∂z}$

故上式最终化为：

其中L_A→E代表流体区域和弹性介质区域的交界面。

由于流体区域和弹性介质区域的交界面需要满足法向应力连续，切向应力为零的边界条件，故

${\widetilde{\mathrm{\τ}}}_{xx}={\widetilde{\mathrm{\τ}}}_{zz}=-\tilde{p}$

${\widetilde{\mathrm{\τ}}}_{xz}=0$

最终化为：

其中L_E→A代表弹性介质区域和流体区域的交界面。

大量数据模拟实例表明，相比于传统PML和C-PML，优化的α和k能够在无精度损失的前提下极大的提高吸收效果。为了实现在TTI介质的模拟中足够精确，令m取一个足够小的值(通常为0.005到0.02之间)，就可以在不损失精度的前提下达到最大的稳定性。

实施例一：

考虑一个TTI声波测井模型，其声源采用的是中心频率为3kHz的单极子声源。井外地层为将一个VTI地层旋转45°所得，其广义Hooke矩阵为在此问题中，吸收边界层只占8个网格的厚度，远远小于波长，这对数值模拟方法是一个非常大的挑战。

图1从左到右分别为采用M-PML，C-PML和本发明的完全吸收PML吸收边界的波场图。图中可以看到，针对该模型，M-PML和C-PML吸收边界下，随着波场的传播有非常明显的边界反射和不稳定，图中虚线所示区域，而采用本发明提出的完全吸收PML边界条件，有更好的吸收效果和更好的稳定性。

接下来考察m取值不同时的吸收效果对比。图2为前文所述TTI介质模型中m取值不同时的频域响应。左图为m＝0.000时，3000Hz频率的响应特征，在这种情况下，在左上角可观察到明显的虚假反射；而在右图m＝0.02的情况下，可看到虚假反射得到了很好的消除，具有更好的吸收效果。

另外考查接收器接收到的波形特征。其实现方法为以100Hz为间隔，共取100个频率点，得到从100Hz到10000Hz的频率响应特征，并通过傅里叶反变换计算得到时域波形。图3为时域理论参考解的波形。图4为当m＝0.000时的时域波形特征，此时可观察到全场存在明显的虚假反射，在6ms前后，八个接收器上都有不同的来自界面的反射信号；而在图5当m＝0.02时的时域波形特征中，全场无虚假反射，且与参考解吻合较好。说明该全吸收PML吸收边界的参数优化对于数值模拟的准确性是非常重要的。

本发明的方法将已被用于一个实际的页岩TTI地层大斜度井问题中，并取得了良好的效果。

Claims (2)

1.适用于TTI介质的频率域有限元全吸收PML方法，用于构造频域有限元全吸收PML边界条件，实现TTI介质复杂井眼条件的高精度声波测井模拟；其特征在于：

所述方法采用复拉伸坐标系的弹性波方程，所述复拉伸坐标系的弹性波方程在所述频域下表示形式为：

$\left\{\begin{array}{c}-\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{u}=\frac{\∂{\widetilde{\mathrm{\τ}}}_{xx}}{\∂x}+\frac{\∂{\widetilde{\mathrm{\τ}}}_{xz}}{\∂z}\\ -\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{w}=\frac{\∂{\widetilde{\mathrm{\τ}}}_{zx}}{\∂x}+\frac{\∂{\widetilde{\mathrm{\τ}}}_{zz}}{\∂z}\end{array}\right.$

在PML区域内，加入拉伸函数，将上述方程表示为：

$\left\{\begin{array}{c}-\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{u}=\frac{1}{s_x}\frac{\∂{\widetilde{\mathrm{\τ}}}_{xx}}{\∂x}+\frac{1}{s_z}\frac{\∂{\widetilde{\mathrm{\τ}}}_{xz}}{\∂z}\\ -\mathrm{\ρ}{\mathrm{\ω}}^2\tilde{w}=\frac{1}{s_x}\frac{\∂{\widetilde{\mathrm{\τ}}}_{zx}}{\∂x}+\frac{1}{s_z}\frac{\∂{\widetilde{\mathrm{\τ}}}_{zz}}{\∂z}\end{array}\right.$

所述PML是完全匹配层，s_x和s_z为PML吸收边界内的拉伸函数，s_x和s_z分别构造x和z方向的两个拉伸函数分量，s_x和s_z的表达式如下：

$s_x=k_x+\frac{d_x+{\mathrm{md}}_z}{{\mathrm{\α}}_x+i\mathrm{\ω}}$

$s_z=k_z+\frac{d_z+{\mathrm{md}}_x}{{\mathrm{\α}}_z+i\mathrm{\ω}}$

上述式中：

为x方向的实数衰减因子，

为x方向复频移因子，

为x方向倏逝波抑制因子，

L_x为吸收层x方向的厚度；

m为切向衰减因子；

为z方向的实数衰减因子，

为z方向复频移因子，

为z方向倏逝波抑制因子，

L_z为吸收层z方向的厚度，

ω为角频率，为频率的2π倍，

i为虚数单位。

2.根据权利要求1所述的适用于TTI介质的频率域有限元全吸收PML方法，其特征在于，所述实数衰减因子复频移因子α₀＝2ω，倏逝波抑制因子中k₀取值范围为1～20，切向衰减因子m取值范围为0.005～0.02。