一种适用于交错网格有限差分的全吸收PML方法

摘要

本发明涉及一种适用于交错网格有限差分的全吸收PML方法，综合PML中的d⊥，C-PML中的α，k，M-PML中的m因素，得到与复杂的大斜度井声波测井吻合的数值模拟结果；将C-PML当中的X、Y、Z向吸收边界的拉伸因子和M-PML当中的X、Y、Z向吸收边界的拉伸因子分别融合为混合X向、Y向、Z向吸收边界的拉伸因子：本发明的有益效果是：通过参数优化，吸收C-PML和M-PML优势，实现了PML吸收边界在求解TTI介质时的高效性。在极端入射角入射的情况下计算精确，对极端各向异性介质也能保持良好的稳定性。

说明书

技术领域

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

背景技术

石油是一种重要的战略资源。从经济和技术的角度，无论是对外合作、通过国际投标购买国外油田，还是在国内进行勘探寻找新的油气区以及对老油田进行重新评价，数值模拟方法都是一个重要的关键技术，由于它可以为实际数据提供理论依据，所以无论是对国民经济的发展，还是石油储备量的合理确定，都有十分重要的意义。

声反射成像测井已经成为一种具有潜力的远探测方法，无论是对地质导向，还是对于大斜度井和水平井环境下的页岩气开采，都起到了重要的作用。在各种常用的声反射成像测井技术的数值模拟方法中，交错网格有限差分方法是一种相对高效的方法，并且可以对复杂的三维模型进行模拟，更接近真实情况。我们在研究大斜度井穿越页岩层(往往是VTI(垂向横向各向同性)介质)所产生的井孔声场时，所遇到的一个较大问题是，有限差分网格为正四边形，在大斜度井中的建模会产生锯齿状网格，导致无法避免的离散误差。通过将参考坐标系进行坐标变换，可以将VTI介质中的模拟问题转换为TTI(水平横向各向同性)介质中。

发明内容

本发明的目的是提出一种适用于交错网格有限差分的全吸收PML方法的技术方案。同时具备稳定性和准确性，使PML吸收边界在求解TTI介质时的实现高效性。

为了实现上述目的，本发明的技术方案是：一种适用于交错网格有限差分的全吸收PML方法，综合PML中的d⊥，C-PML中的α，k，M-PML中的m因素，得到与复杂的大斜度井声波测井吻合的数值模拟结果；

将C-PML当中的X向吸收边界的拉伸因子和M-PML当中的X向吸收边界的拉伸因子融合为一个混合X向吸收边界的拉伸因子：

$>s_x=k_x+\frac{d_x+m_{x/y}{d}_y+m_{x/z}{d}_z}{{\alpha }_x+\mathrm{i\omega }},$>

将C-PML当中的Y向吸收边界的拉伸因子和M-PML当中的Y向吸收边界的拉伸因子融合为一个混合Y向吸收边界的拉伸因子：

$>s_y=k_y+\frac{d_y+m_{y/x}{d}_x+m_{y/z}{d}_z}{{\alpha }_y+\mathrm{i\omega }},$>

将C-PML当中的Z向吸收边界的拉伸因子和M-PML当中的X向吸收边界的拉伸因子融合为一个混合X向吸收边界的拉伸因子：

$>s_z=k_z+\frac{d_z+m_{z/x}{d}_x+m_{z/y}{d}_y}{{\alpha }_z+\mathrm{i\omega }},$>

上述式中：d_x为PML在X向的衰减函数，d_y为PML在Y向的衰减函数，d_z为PML在Z向的衰减函数，

α_x为X向复频移函数，

α_y为Y向复频移函数，

α_z为Z向复频移函数，

k_x为X向抑制倏逝波的函数，

k_y为Y向抑制倏逝波的函数，

k_z为Z向抑制倏逝波的函数，

m_x/y为X向与Y向间的切向衰减控制因子，

m_x/z为X向与Z向间的切向衰减控制因子，

m_y/x为Y向与X向间的切向衰减控制因子，

m_y/z为Y向与Z向间的切向衰减控制因子，

m_z/x为Z向与X向间的切向衰减控制因子，

m_z/y为Z向与Y向间的切向衰减控制因子，

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

i为虚数单位。

更进一步，所述切向衰减控制因子m_x/y、m_x/z、m_y/x、m_y/z、m_z/x、m_z/y的取值范围为0.005～0.020。

本发明的有益效果是：通过参数优化，兼具C-PML和M-PML的优势，实现了PML吸收边界在求解TTI介质时的高效性。在极端入射角入射的情况下计算精确，对极端各向异性介质也能保持良好的稳定性；可准确地对TTI介质进行数值模拟，获得复杂的大斜度井模型。

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

附图说明

图1是本发明一个各向同性均一模型分别采用C-PML和M-PML的效果图，上侧为C-PML的吸收效果，下侧为M-PML的吸收效果；

图2是本发明一个极端各向异性均一模型分别采用C-PML和M-PML的效果图，上侧为C-PML的吸收效果，下侧为M-PML的吸收效果；

图3是一个极端TTI各向异性井眼模型图；

图4是本发明当m＝0.000(即只采用C-PML时)的波场传播情况；

图5是本发明当m＝0.005(即采用全吸收融合PML时)的波场传播情况；

图6是当模型较大时的参考解的波形图；

图7是本发明当m＝0.000(即只采用C-PML时)的波形图；

图8是本发明当m＝0.005(即采用全吸收融合PML时)，吸收参数恰当时的波形图；

图9是本发明当m＝0.010(即采用全吸收融合PML时)，吸收参数不恰当时的波形图。

具体实施方式

声反射成像测井已经成为一种具有潜力的远探测方法，在各种常用的声反射成像测井技术的数值模拟方法中，交错网格有限差分方法是一种相对高效的方法，并且很容易应用到研究复杂的三维情况中。大斜度井穿越页岩层(往往是VTI介质)所产生的井孔声场时，由于有限差分网格为正四边形，因此大斜度井的建模会产生锯齿状，并产生无法避免的离散误差。

在本发明中，我们采用一种方法来避免这种问题，它是通过旋转参考坐标系，使得井轴总是与直角坐标系中的其中一个坐标轴保持一致，这样VTI介质将通过Bond变换矩阵转化为TTI介质。于是我们的这种方法将可以广泛适用于(1)考查VTI对称轴与井轴之间不同角度的影响；(2)考查井轴与井外地层界面不同角度的影响。

本发明主要致力于解决PML(完全匹配吸收层)吸收边界在求解TTI介质时的高效性。C-PML(卷积完全匹配吸收层)在极端入射角入射的情况下依旧计算精确，然而它在某些极端各向异性介质中会出现数值不稳定。另一方面，M-PML(多轴完全匹配吸收层)即使对极端各向异性介质也能保持卓越的稳定性，但在PML层厚度不够厚，以及修正因子不够优化的情况下，却牺牲了准确性。而过厚的PML层厚度又大幅降低了计算效率。为了达到稳定性和准确性的共赢，我们引入一种混合型吸收边界，通过参数优化，使两种PML的优势能够共存。

本发明的特点是同时反映了PML中的d⊥，C-PML中的α，k，M-PML中的d//等因素的影响，可以最大限度克服稳定性和准确性难以兼顾的难题。同时这种新型吸收边界在吸收层较薄的情况下(10个网格左右的厚度，其厚度远远小于波长)依然可以获得较好的吸收效果，故计算效率也优于PML、C-PML和M-PML的任何一个。在这种方法中，传统的PML是在此新型PML下的特例，即本方法涵括了PML、C-PML和M-PML，同时在上述PML之一不能适用的情况下，本方法依然适用。

三维交错网格有限差分方法中的速度-应力方程为：

$>-\mathrm{i\omega }{\tau }_{\mathrm{xx}}=\frac{1}{s_x}\frac{\partial v_x}{\partial x}+\frac{1}{s_y}\frac{\partial v_y}{\partial y}+\frac{1}{s_z}\frac{\partial v_z}{\partial z}$>

其中以s_x为例，它为PML吸收边界的拉伸因子。在最为传统的PML当中，而在C-PML和M-PML当中，它分别被修正为和$>s_x=1+\frac{d_x+m_{x/y}{d}_y+m_{x/z}{d}_z}{\mathrm{i\omega }}.$>同理可以得到s_y和s_z。

本发明的技术方案是：一种适用于交错网格有限差分的全吸收PML方法。综合PML中的d⊥，C-PML中的α，k，M-PML中的m因素，得到与复杂的大斜度井声波测井吻合的数值模拟结果；

将C-PML当中的X向吸收边界的拉伸因子和M-PML当中的X向吸收边界的拉伸因子融合为一个混合X向吸收边界的拉伸因子：

$>s_x=k_x+\frac{d_x+m_{x/y}{d}_y+m_{x/z}{d}_z}{{\alpha }_x+\mathrm{i\omega }};$>

将C-PML当中的Y向吸收边界的拉伸因子和M-PML当中的Y向吸收边界的拉伸因子融合为一个混合Y向吸收边界的拉伸因子：

$>s_y=k_y+\frac{d_y+m_{y/x}{d}_x+m_{y/z}{d}_z}{{\alpha }_y+\mathrm{i\omega }};$>

将C-PML当中的Z向吸收边界的拉伸因子和M-PML当中的X向吸收边界的拉伸因子融合为一个混合X向吸收边界的拉伸因子：

$>s_z=k_z+\frac{d_z+m_{z/x}{d}_x+m_{z/y}{d}_y}{{\alpha }_z+\mathrm{i\omega }}.$>

上述式中：d_x为PML在X向的衰减函数，d_y为PML在Y向的衰减函数，d_z为PML在Z向的衰减函数；

α_x为X向复频移函数，对应于C-PML中的X向复频移函数；

α_y为Y向复频移函数，对应于C-PML中的Y向复频移函数；

α_z为Z向复频移函数，对应于C-PML中的Z向复频移函数；

k_x为X向抑制倏逝波的函数，对应于C-PML中的X向抑制倏逝波的函数；

k_y为Y向抑制倏逝波的函数，对应于C-PML中的Y向抑制倏逝波的函数；

k_z为Z向抑制倏逝波的函数，对应于C-PML中的Z向抑制倏逝波的函数；

m_x/y为X向与Y向间的切向衰减控制因子；对应于M-PML中X向与Y向间的切向衰减控制因子；

m_x/z为X向与Z向间的切向衰减控制因子；对应于M-PML中X向与Z向间的切向衰减控制因子；

m_y/x为Y向与X向间的切向衰减控制因子；对应于M-PML中Y向与X向间的切向衰减控制因子；

m_y/z为Y向与Z向间的切向衰减控制因子；对应于M-PML中Y向与Z向间的切向衰减控制因子；

m_z/x为Z向与X向间的切向衰减控制因子；对应于M-PML中Z向与X向间的切向衰减控制因子；

m_z/y为Z向与Y向间的切向衰减控制因子；对应于M-PML中Z向与Y向间的切向衰减控制因子；

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

i为虚数单位。

所述切向衰减控制因子m_x/y、m_x/z、m_y/x、m_y/z、m_z/x、m_z/y的取值范围为0.005～0.020。

本发明在实现C-PML的基础上，将其中的d_x修正为d_x+m_x/yd_y+m_x/zd_z，并同理修改d_y和d_z为相应的形式，即实现了两种方法的有机融合。

大量实验表明，相比于传统PML和C-PML，尽管优化的α和k能够在无精度损失的前提下大幅提高吸收效果，但在TTI介质的模拟中仍不足够精确。于是需要M-PML中的m(包括m_x/y、m_x/z、m_y/x、m_y/z、m_z/x、m_z/y)取一个足够小的值，经过研究与实验，m值取0.005到0.020，就可以在最小的精度损失中达到最大的稳定性。

实施例一：

我们考虑两个较为典型的均一模型，如图1和图2所示。其中图1为各向同性介质(其广义Hooke矩阵为$>C=\left(\begin{array}{ccc}4& & \\ & 4& \\ & & 2\end{array}\right)\mathrm{GPa}$>)的波场快照，此时C-PML和M-PML均是稳定的，C-PML吸收效果较好，而M-PML产生了部分虚假反射。图2为极端各向异性介质(其广义Hooke矩阵为$>C=\left(\begin{array}{ccc}4& 7.5& \\ 7.5& 20& \\ & & 2\end{array}\right)\mathrm{GPa}$>)的波场快照，此时C-PML是不稳定的，M-PML是稳定的，C-PML吸收效果较好，而M-PML产生了部分虚假反射。此例子印证了两种PML吸收边界融合的必要性。

我们进一步考虑一个典型的声波测井模型，如图3所示。模型大小为0.6m×6.0m，井内流体内有一实心刚性仪器存在，声源采用的是中心频率为8kHz的偶极子声源。井外地层为将一个VTI地层旋转45°所得，其广义Hooke矩阵为$>C=\left(\begin{array}{ccc}17.47& 11.93& 2.14\\ 11.93& 17.47& 2.14\\ 2.14& 2.14& 4.91\end{array}\right)\mathrm{GPa}.$>另外，我们在这个问题中，吸收边界层只占10个网格的厚度，使得其厚度远远小于波长，这对我们的数值模拟结果是一个非常大的挑战。

图4和图5为每1ms间隔所得到的波场快照，可见当m＝0.000时(即C-PML情形)，模型中出现了数值不稳定，而当m＝0.005时，模拟效果较好。

我们进一步考查模拟波形的结果。由于当TI介质的对称轴与井轴不一致，难以求取解析解，因此本专利与参考解进行相互验证。参考解是通过将模型设置的较大，从而使得PML吸收边界不影响模拟结果而求取的。如图6-9所示。其中当m＝0.000时，数值不稳定导致T＝3.0ms时出现了较大的误差。当m＝0.005时，两者吻合较好。而m＝0.010时，两者吻合度略差于m＝0.005的情况，说明该混合吸收边界的参数优化对于数值模拟的准确性是至关重要的。

本发明的方法已被应用于一个实际的三维大斜度井问题中，并取得了良好的效果。