阅读说明：本技术 一种二维起伏地表下的sh波曲线网格模拟方法 (SH wave curve grid simulation method under two-dimensional undulating surface ) 是由 刘炯 刘喜武 刘宇巍 霍志周 钱恪然 张金强 于 2019-03-01 设计创作，主要内容包括：本发明涉及一种二维起伏地表下的SH波曲线网格模拟方法,包括以下步骤：步骤1：建立不规则的物理空间和规则的计算空间之间的映射关系；步骤2：将所述规则的计算空间中的SH波进行波场模拟；步骤3：将所述波场模拟的结果通过所述映射关系映射到所述物理空间中,得到所述物理空间中SH波沿起伏地表的波场模拟结果。本发明方法稳定且不会产生传统直角网格“阶梯”离散产生的数值噪音,起伏地表的新模拟精度较高。在解决自由面模拟的同时,保证了模拟的高精度。本发明较好解决了起伏地表用传统直角网格模拟时认为数值噪音误差,不会出现模拟不稳定的现象,有效保证了波场数值模拟的精度。(The invention relates to a method for simulating SH wave curve grids under a two-dimensional undulating surface, which comprises the following steps: step 1: establishing a mapping relation between an irregular physical space and a regular calculation space; step 2: carrying out wave field simulation on the SH wave in the regular calculation space; and step 3: and mapping the result of the wave field simulation to the physical space through the mapping relation to obtain a wave field simulation result of the SH waves in the physical space along the undulating surface. The method is stable, numerical noise generated by the 'ladder' dispersion of the traditional right-angle grid can not be generated, and the new simulation precision of the undulating surface is higher. The method solves the problem of free surface simulation and ensures the high precision of simulation. The method well solves the problem that numerical noise errors are considered when the undulating surface is simulated by using the traditional right-angle grid, the phenomenon of unstable simulation cannot occur, and the precision of numerical simulation of the wave field is effectively ensured.)

一种二维起伏地表下的SH波曲线网格模拟方法

技术领域

本发明属于油汽地球物理技术领域，具体涉及一种二维起伏地表下的SH波曲线网格模拟方法。

背景技术

根据波质点振动方向与波的传播方向的不同，波可以分为横波(即S波)与纵波(即P波)。质点振动方向与传播方向相同的波为纵波；质点振动方向与传播方向垂直的波称为横波。质点振动发生在与波的传播面相垂直的面内的波为SV波,质点振动发生在与波的传播面相平行的面内的波为SH波。

而起伏地表是山区地震勘探的最重要的干扰因素之一，采用数值模拟的方法来识别和消除这种干扰可以提高地震勘探的效果，因此国内外许多学者利用波动方程数值模拟的开展了起伏地表下的地震波传播规律的研究。

有限元法可以很好描述不规则地表自由面的形状，且自由面条件在有限元法模拟中可以自动满足，因此原则上它很适合自由表面复杂的问题。但是低阶的有限元法会有较强的频散，而一般的高阶有限元会产生伪波，且有限元法对内存要求高，计算量大，因此在计算规模较大的油气勘探领域中运用得相对比较少。

有限差分法和伪谱法是目前地球物理领域地震模拟中最经常用到的两种方法。传统的有限差分法和伪谱法在直角网格上进行模拟。许多研究者首将这两种模拟方法直接应用至起伏地表的地震波模拟中。它们先近似然后将自由表面上一部分空间补足，使整体成为一个规则的矩形空间，然后用直角网格来离散计算模型空间，这导致起伏地表被离散成数学模型中的“阶梯”形状；最后在规则的模型空间采用传统有限差分法或伪谱法来模拟地震波的传播。然而，采用传统直角网格对起伏地表进行阶梯离散会对模型的外形产生较大的误差，这会对物理量的模拟产生较大的误差，而且离散的“阶梯”会产生较强的认为数值散射。此外这种在起伏的自由界面的处理，一般隐含地表和上部空间的边界条件，又由于上部空间的物理性质和实际地下性质差异明显，这会导致模拟过程中容易出现不稳定的现象，使得模拟过程失败。而且界面的精度最多只能达到二阶。

针对上述传统直角坐标离散起伏地表时出现的问题，国外学者Tessemer(1992)引入Fornberg(1988)地下不规则表面的处理思想，利用坐标映射技术将不规则的自由表面映射至规则平表面上，然后进行计算。随后国内外一些学者利用相同的方式进行了起伏地表下转换波(P-SV波)的数值模拟。但这些方法只在纵向方向进行简单的坐标映射，水平方向则保持不变，当地表起伏比较剧烈时，这种处理可能会导致物理空间的曲线网格歪度过大，从而导致最后较大的误差。

发明内容

为了解决上述问题，本发明提供了一种二维起伏地表下的SH波曲线网格模拟方法，该方法稳定且不会产生传统直角网格“阶梯”离散产生的数值噪音，起伏地表的新模拟精度较高。在解决自由面模拟的同时，保证了模拟的高精度。本发明较好解决了起伏地表用传统直角网格模拟时认为数值噪音误差，不会出现模拟不稳定的现象，有效保证了波场数值模拟的精度。

根据本发明提供了一种二维起伏地表下的SH波曲线网格模拟方法，包括以下步骤：

步骤1：建立不规则的物理空间和规则的计算空间之间的映射关系；

步骤2：将所述规则的计算空间中的SH波进行波场模拟；

步骤3：将所述波场模拟的结果通过所述映射关系映射到所述物理空间中，得到所述物理空间中SH波沿起伏地表的波场模拟结果。

在一个实施例中，所述步骤1中，通过生成曲线网格建立所述映射关系。

在一个实施例中，所述曲线网格的生成方法包括代数法、保角变换法、微分方程法。

在一个实施例中，所述曲线网格的生成的步骤包括：

(1)建立泊松方程组；

(2)根据泊松方程组，得到生成曲线网格的控制方程组；

(3)通过求解生成曲线网格的控制方程组的数值，建立所述物理空间和所述计算空间的映射关系。

在一个实施例中，所述步骤(1)中的泊松方程组为：

式中，x、z为不规则的物理空间中物理坐标(xoz)的自变量，ξ、η为规则的计算空间中曲线网格坐标(ξoη)的自变量；P、Q为坐标控制函数，用于实现物理坐标(xoz)中网格分布的控制；

利用物理坐标和曲线网格坐标ξ＝ξ(x,z)、η＝η(x,z)之间的映射关系改变所述泊松方程组中的自变量(x,z)为(ξ,η)，其中，ξ＝ξ(x,z)、η＝η(x,z)为将方程中的自变量x、z表达的函数用ξ、η来表示。

在一个实施例中，所述步骤(2)具体为：利用复合函数求导法则，对泊松方程组进行求导，得到生成曲线网格的控制方程组：

α·x_ξξ-2β·x_ξη+γ·x_ηη+J²(P·x_ξ+Q·x_η)＝0

α·z_ξξ-2β·z_ξη+γ·z_ηη+J²(P·z_ξ+Q·z_η)＝0

式中，

β＝x_ξx_η+z_ξz_η，x_ξ为x关于ξ的一阶导数；x_η为x关于η的一阶导数；x_ξξ为x关于ξ的二阶导数；x_ξη为x关于ξ和η的二阶偏导数；x_ηη为x关于η的二阶导数；z_ξ为z关于ξ的一阶导数；z_η为z关于η的一阶导数；z_ξξ为z关于ξ的二阶导数；z_ξη为z关于ξ和η的二阶偏导数；z_ηη为z关于η的二阶导数；自变量ξ，η定义在规则矩形区域上；J为从曲线网格坐标到物理坐标的雅克比行列式：

所述步骤(3)中，可以通过有限差分法、有限元法或伪谱法对曲线网格的控制方程组进行数值求解，得到ξ＝ξ(x,z)，η＝η(x,z)。

在一个实施例中，所述步骤2中，又具体包括以下步骤：

步骤2.1：建立SH波在物理坐标系下传播的控制方程组；

步骤2.2：根据物理坐标系下传播的控制方程组，得到SH波在曲线网格坐标系下传播的控制方程组；

步骤2.3：在规则的所述计算空间中利用数值方法，对曲线网格坐标系下传播的控制方程组进行波场模拟。

在一个实施例中，在步骤2.1中，所述SH波在物理坐标系下传播的控制方程组的建立步骤为：根据弹性力学原理，建立二维x-z空间中固体介质的SH波传播方程组：

其中，V表示固体介质沿垂直于x-z平面的y方向的速度；τ_xy、τ_yz是应力分量；f_y是沿y方向的外力；t为时间；ρ代表固体介质的密度；μ是固体介质的剪切模量。

在一个实施例中，所述步骤2.2中，依据复合函数求导法则，对物理坐标系下传播的控制方程组进行求导，得到曲线网格坐标系下的SH波方程组：

式中，f_y(ξ,η)为用曲线坐标(ξoη)中沿空间y方向的外力；

为物理坐标系变换到曲线网格坐标系下的度量系数，其中，

在一个实施例中，所述步骤2.3中的数值方法包括有限差分或伪谱法。

与现有技术相比，本发明的优点在于：本发明引入不规则模型面的曲线网格剖分技术，提出了一般曲线坐标下的SH波模拟方法。本发明方法首先通过一种新的曲线网格剖分法，将实际起伏地表物理空间光滑映射至规则的矩形计算空间，建立起实际物理空间和规则计算空间的对应关系，然后在规则计算空间中对SH波进行模拟，最后将计算空间的模拟结果映射回实际物理空间中，得到SH波沿不规则起伏表面的波场传播过程。本发明方法稳定且不会产生传统直角网格“阶梯”离散产生的数值噪音，起伏地表的新模拟精度较高。在解决自由面模拟的同时，保证了模拟的高精度。本发明较好解决了起伏地表用传统直角网格模拟时认为数值噪音误差，不会出现模拟不稳定的现象，有效保证了波场数值模拟的精度。

附图说明

在下文中将基于实施例并参考附图来对本发明进行更详细的描述。其中：

图1为本发明不规则物理空间曲线网格离散原理图；

图2为本发明二维倾斜地表模型图(实施例一)；

图3为本发明倾斜地表的曲线网格离散图(实施例一)；

图4为本发明倾斜地表0.55秒时刻的波场快照图(实施例一)；

图5至图7分别为本发明3个虚拟检波器数值模拟结果和理论结果对比(实施例一)；

图8为本发明不规则地堑模型图(实施例二)；

图9为本发明不规则地堑模型的曲线网格离散图(实施例二)；

图10为本发明地堑模型0.2秒时刻的波场快照图(实施例二)。

在附图中相同的部件使用相同的附图标记。附图并未按照实际的比例绘制。

具体实施方式

下面将结合附图对本发明作进一步说明。借此对本发明如何应用技术手段解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不存在冲突，各个实施例中所提到的各项技术特征均可以任意方式组合起来。本发明并不局限于文中公开的特定实施例，而是包括落入权利要求的范围内的所有技术方案。

本发明提供了一种二维起伏地表下的SH波曲线网格模拟方法，包括以下步骤：

步骤1：建立不规则的物理空间和规则的计算空间之间的映射关系；

步骤2：将所述规则的计算空间中的SH波进行波场模拟；

步骤3：将所述波场模拟的结果通过所述映射关系映射到所述物理空间中，得到所述物理空间中SH波沿起伏地表的波场模拟结果。

其中，所述步骤1中，所述映射关系的建立就是曲线网格的生成过程，通过生成曲线网格建立所述映射关系。将曲线坐标网格技术应用至波场模拟中。即，首先将不规则的起伏地表(物理空间)的计算区域用曲线网格进行离散(如图1所示)，对不规则的空间用曲线网格进行离散(如图1左图所示)，同时建立起不规则的物理空间中物理坐标(xoz)和规则的计算空间中曲线网格计算坐标(ξoη)(如图1右图所示)之间的映射关系。波场模拟在规则的计算空间中进行，最终模拟结果会根据映射关系返回至物理空间中，从而得到的不规则的物理空间中SH波沿起伏地表的波场模拟结果。

所述曲线网格的生成方法包括代数法、保角变换法、微分方程法。

在一个实施例中，利用解泊松方程组(属微分方程法)生成曲线网格的方法，具体地，所述曲线网格的生成的步骤包括：

(1)建立泊松方程组；

(2)根据泊松方程组，得到生成曲线网格的控制方程组；

(3)通过求解生成曲线网格的控制方程组的数值，建立所述物理空间和所述计算空间的映射关系。

在一个实施例中，所述步骤(1)中的泊松方程组为：

式中，x、z为不规则的物理空间中物理坐标(xoz)的自变量，ξ、η为规则的计算空间中曲线网格坐标(ξoη)的自变量；P、Q为坐标控制函数，用于实现物理坐标(xoz)中网格分布的控制；通过对P、Q坐标控制函数的选择，实现不同曲线坐标分布，P、Q可以按需求选择。

在一个实施例中，所述步骤(2)具体为：利用复合函数求导法则，对泊松方程组进行求导，得到生成曲线网格的控制方程组：

α·x_ξξ-2β·x_ξη+γ·x_ηη+J²(P·x_ξ+Q·x_η)＝0

α·z_ξξ-2β·z_ξη+γ·z_ηη+J²(P·z_ξ+Q·z_η)＝0

式中，

β＝x_ξx_η+z_ξz_η，x_ξ为x关于ξ的一阶导数；x_η为x关于η的一阶导数；x_ξξ为x关于ξ的二阶导数；x_ξη为x关于ξ和η的二阶偏导数；x_ηη为x关于η的二阶导数；z_ξ为z关于ξ的一阶导数；z_η为z关于η的一阶导数；z_ξξ为z关于ξ的二阶导数；z_ξη为z关于ξ和η的二阶偏导数；z_ηη为z关于η的二阶导数；自变量ξ，η定义在规则矩形区域上，其中，所谓规则矩形区域为：自变量ξ，η的范围是设计在矩形区域上的，是规则的，其边界不是非规则的。J为从曲线网格坐标到物理坐标的雅克比行列式：

所述步骤(3)中，可以通过各种方法(包括有限差分法、有限元法、伪谱法等)对曲线网格的控制方程组进行数值求解，得到ξ＝ξ(x,z)，η＝η(x,z)，也即是物理坐标点和曲线网格坐标点的映射关系。利用物理坐标和曲线网格坐标之间的映射关系ξ＝ξ(x,z)、η＝η(x,z)改变SH波方程组中的自变量(x,z)为(ξ,η)。然后在曲线网格坐标中求解坐标控制方程，再利用上述泊松方程组得到的物理坐标点和曲线网格坐标点的映射关系，得到物理空间中的待求波场量。

在一个实施例中，所述步骤2中，又具体包括以下步骤：

步骤2.1：建立SH波在物理坐标系下传播的控制方程组；

步骤2.2：根据物理坐标系下传播的控制方程组，得到SH波在曲线网格坐标系下传播的控制方程组；

步骤2.3：在规则的所述计算空间中利用数值方法，对曲线网格坐标系下传播的控制方程组进行波场模拟。所述步骤2.3中的数值方法包括有限差分或伪谱法。

在一个实施例中，在步骤2.1中，所述SH波在物理坐标系下传播的控制方程组的建立步骤为：根据弹性力学原理，建立二维x-z空间中固体介质的SH波传播方程组：

其中，V表示固体介质沿垂直于x-z平面的y方向的速度；τ_xy、τ_yz是应力分量；f_y是沿y方向的外力；t为时间；ρ代表固体介质的密度；μ是固体介质的剪切模量。

在一个实施例中，所述步骤2.2中，依据复合函数求导法则，对物理坐标系下传播的控制方程组进行求导，得到曲线网格坐标系下的SH波方程组：

式中，f_y(ξ,η)为用曲线坐标(ξoη)中沿空间y方向的外力；

为物理坐标系变换到曲线网格坐标系下的度量系数，这些度量系数可以通过如下公式在曲线网格建立时数值求取得到，其中，

实施例一

为对二维倾斜的地表进行SH波波场模拟，其中，图2为本发明二维倾斜地表模型图，所示图2中安置了3个虚拟检波器A、B、C来记录空间点的位移。

图3是倾斜地表的曲线网格离散图。

模拟中震源采用30Hz的Ricker子波；模型左、右、下边采用吸收边界。图4是倾斜地表在0.55秒时刻的波场快照图；图4中S表示直达横波，SS表示经地表的反射横波。

图5至图7是用曲线网格数值模拟得到的检波器位移记录和理论结果的对比。可以看出，数值解和理论解符合精度高，这证明本发明方法的模拟精度高。

实施例二

图8是不规则地堑模型图。图9是不规则地堑模型的曲线网格离散图。图10是地堑模型在0.2秒时刻的波场快照图。可以看到通过曲线网格的模拟，在地表附近未见明显的数值散射，有效保证了波场数值模拟的精度。

通过实施例一和实施例二的实验结果表明，本发明方法稳定且不会产生传统直角网格“阶梯”离散产生的数值噪音，起伏地表的新模拟精度较高。在解决自由面模拟的同时，保证了模拟的高精度。本发明较好解决了起伏地表用传统直角网格模拟时认为数值噪音误差，不会出现模拟不稳定的现象，有效保证了波场数值模拟的精度。

虽然已经参考如上优选实施例对本发明进行了描述，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。