规则网格速度模型的自适应非结构三角网格化方法

摘要

一种规则网格速度模型的自适应非结构三角网格化方法，该方法将网格生成过程类比为骨架平衡的过程，生成网格大小随速度自适应分布的高质量网格。此外该方法考虑速度模型的二阶梯度场的信息，将网格节点映射到规则网格模型定义的隐式物性分界面上，使非结构网格单元的边与速度模型中的物性分界面自适应贴合，实现了规则网格模型向非结构网格的准确转化，达到保存重要物性分界面的目的，为有限元的数值模拟方法提供可靠的非结构网格模型。

权利要求

1.一种规则网格速度模型的自适应非结构三角网格化方法，包括以下步骤：

步骤1：定义网格密度函数

步骤2：根据所述网格密度函数 在所述规则网格速度模型区域内随机分布初始的网格节点；

步骤3：对所述网格节点进行Delaunay非结构三角网格剖分，并将非结构三角网格的节点坐标定义为二维数组变量p＝[x z]；

步骤4：将所述非结构三角网格类比于骨架结构，定义所述非结构三角网格中每个节点所受到的骨架臂的力F(p)；

步骤5：对所述规则网格速度模型进行二阶求导，获得所述规则网格速度模型的二阶梯度场，根据所述二阶梯度场设置二阶梯度力场imgF，所述二阶梯度力场imgF正比于所述二阶梯度场，然后通过二维最相邻插值获得所述非结构三角网格的每个网格节点受到的二阶梯度场力imgF(p)；

步骤6：构建所述非结构三角网格节点坐标和节点合力关于时间的偏微分方程系统，其中所述节点合力由所述骨架臂的力F(p)和所述二阶梯度场力imgF(p)构成；

步骤7：离散化所述偏微分方程系统，获得所述非结构三角网格节点坐标的更新关系式，根据现在时刻的节点坐标及节点合力，计算下一个时刻的节点坐标，获得更新的网格节点；

步骤8：判断所述偏微分方程系统是否趋于平衡，如果所述偏微分方程系统还没有趋于平衡，返回步骤3，如果所述偏微分方程系统趋于平衡，继续到步骤9；

步骤9：针对最终更新的网格节点进行Delaunay非结构三角网格剖分，获得最终的非结构三角网格单元，根据所述规则网格速度模型二维最相邻插值获得所述最终的非结构三角网格单元的物性参数，输出所述网格单元的网格节点坐标、网格单元编号及物性参数。

2.根据权利要求1所述的规则网格速度模型的自适应非结构三角网格化方法，其中，在所述步骤1中，根据所述规则网格速度模型的速度分布以及地震波有限元法模拟中的空间采样需求，定义网格密度函数

3.根据权利要求1所述的规则网格速度模型的自适应非结构三角网格化方法，其中，在所述步骤1中，所述网格密度函数 其中vs为横波速度，fmax为在地震波有限元法中所用子波的最大频率，N表示一个最短地震波波长内所需的采样点数。

4.根据权利要求1所述的规则网格速度模型的自适应非结构三角网格化方法，其中，在所述步骤4中，根据所述骨架的拓扑结构定义非结构三角网格中每个节点所受到的骨架臂的力F(p)，所述力F(p)包括所述非结构三角网格中每个节点在水平和垂直方向上的受力：其中，Fint表示来自所述骨架臂的内部力函数f(l,l0)，Fext表示来自边界的外力，F(p)的第一列为x分量，第二列为z分量。

5.根据权利要求4所述的规则网格速度模型的自适应非结构三角网格化方法，其中，所述内部力函数f(l,l0)根据以下公式计算：其中l为所述骨架臂的实际长度，l0为期望的非结构三角网格长度，k为弹性系数。

6.根据权利要求1所述的规则网格速度模型的自适应非结构三角网格化方法，其中所述非结构三角网格节点在所述节点合力的作用下发生移动。

7.根据权利要求1所述的规则网格速度模型的自适应非结构三角网格化方法，其中，在所述步骤6中，所述偏微分方程系统为：所述偏微分方程系统的初始条件为p(0)＝p0。

8.根据权利要求1所述的规则网格速度模型的自适应非结构三角网格化方法，其中，在所述步骤7中，所述非结构三角网格节点坐标的更新关系式为：pn+1＝pn+Δt·(F+imgd·n+1 nimgF)，其中p 为n+1时刻的节点坐标，p为n时刻的节点坐标，Δt为更新的时间步长，imgd为调节参数。

9.根据权利要求1所述的规则网格速度模型的自适应非结构三角网格化方法，其中，在所述步骤7中，通过向前欧拉时间离散格式离散化所述偏微分方程系统。

说明书

规则网格速度模型的自适应非结构三角网格化方法

技术领域

[0001] 本公开涉及有限元地震波数值模拟方法，特别涉及一种规则网格速度模型的自适应非结构三角网格化方法。

背景技术

[0002] 目前，在我国南方和西部地区，油气地震勘探的重点正转向丘陵、山前构造带等复杂地区。这些地区地表条件异常复杂，地形起伏剧烈，高差变化非常大，岩性速度变化大导致近地表结构不均匀性严重，同时地下构造复杂，如褶皱强烈、断层发育、构造陡峻、地层变化大等。这些特点导致这些地区地震资料信噪比低及静校正难等各方面的问题。从根本上解决这些勘探问题，需要对起伏地表条件下地震波的传播规律和波场特征进行深入的认识，而有限元法是进行复杂地表和复杂构造地区地震波数值模拟最有效的技术手段。

[0003] 在进行地震波数值模拟时，首先需要对地球物理-地质模型(包括纵横波速度、密度等)进行空间上的网格化离散。有限差分法一般使用规则网格，而有限元法则多使用非结构网格。目前在地震勘探领域中，速度模型大都是定义在规则网格上，以二维矩阵数据的形式给出，其主要原因是有限差分法仍然是目前使用最广泛的地震波数值模拟技术，且规则网格的数据结构简单，便于广泛地传播使用。但规则网格定义的速度模型存在着模型的过采样问题及缺乏显式的边界等问题。非结构网格是指没有规则拓扑关系的网格。它能够对模型进行贴体的网格剖分。特别地，三角形网格能够模拟任意复杂的地形起伏及复杂的地下构造，因此能够更好地刻画模型中的界面信息。此外，采用非结构网格能够适应地下介质的该物性参数变化，自适应的调整网格大小和密度，而不会造成过采样的问题。因此，对于具有复杂结构的模型，非结构网格拥有更强的适应能力。

[0004] 由于目前速度模型基本上是定义在规则网格上，而有限元法使用非结构网格，因此需要考虑如何将规则网格上定义的模型三角网格化。此外，采用间断Galerkin等高阶有限元数值模拟方法时，为了保证模拟精度和计算效率，对使用的非结构网格剖分质量也有较高的需求。

[0005] 对于速度等地球物理-地质模型，最感兴趣的是地下介质物性参数突变的地方。在进行地震波数值模拟时，需要充分尊重这些界面信息的存在。当采用常规思路进行网格剖分时，需要将规则网格定义速度模型转换为数字图像，然后进行图像预处理，勾勒出模型内部边界线(segment)，即物理模型向几何模型的转化过程，最后对几何模型进行网格剖分。对于简单的块状模型，因为其界面清晰，区域划分容易，采用图像处理的思路不失为一种好的方法。但是，当地下模型复杂，界面信息凌乱时，使用图像处理进行区域划分比较麻烦。并且，发明人发现，常规思路的几个步骤是相对独立的，无法形成一个独立的、自适应的网格剖分工具。

[0006] 网格的生成通常被认为是一项复杂的任务，因此网格生成器也总是以“黑匣子”的形式独立于有限元代码的编写环节。Distmesh是由Per-Olof Persson(2005)提出的一种非常简单的三角形网格生成技术。区别于传统的基于几何概念的网格生成方法，Distmesh算法基于三角形网格与钢骨架结构的类比，以及钢骨架结构的整体平衡性，能够提供高质量的网格剖分结果，同时能够物性分布自适应地改变网格密度，但Distmesh算法难以保证剖分的非结构网格与物性分界面的自适应贴体化。基于上述原因，期待一种适用于山前带等复杂(地表/地下)地区的规则网格速度模型的自适应非结构三角网格化方法。

发明内容

[0007] 本公开的目的是建立一种适用于复杂地区的能够有效匹配规则网格模型的非结构三角网格化方法。

[0008] 本公开所采用的解决方案如下：

[0009] 一种规则网格速度模型的自适应非结构三角网格化方法，包括以下步骤：

[0010] 步骤1：定义网格密度函数

[0011] 步骤2：根据所述网格密度函数 在所述规则网格速度模型区域内随机分布初始的网格节点；

[0012] 步骤3：对所述网格节点进行Delaunay非结构三角网格剖分，并将非结构三角网格的节点坐标定义为二维数组变量p＝[x z]；

[0013] 步骤4：将所述非结构三角网格类比于骨架结构，定义所述非结构三角网格中每个节点所受到的骨架臂的力F(p)；

[0014] 步骤5：对所述规则网格速度模型进行二阶求导，获得所述规则网格速度模型的二阶梯度场，根据所述二阶梯度场设置二阶梯度力场imgF，所述二阶梯度力场imgF正比于所述二阶梯度场，然后通过二维最相邻插值获得所述非结构三角网格的每个网格节点受到的二阶梯度场力imgF(p)；

[0015] 步骤6：构建所述非结构三角网格节点坐标和节点合力关于时间的偏微分方程系统，其中所述节点合力由所述骨架臂的力F(p)和所述二阶梯度场力imgF(p)构成；

[0016] 步骤7：离散化所述偏微分方程系统，获得所述非结构三角网格节点坐标的更新关系式，根据现在时刻的节点坐标及节点合力，计算下一个时刻的节点坐标，获得更新的网格节点；

[0017] 步骤8：判断所述偏微分方程系统是否趋于平衡，如果所述偏微分方程系统还没有趋于平衡，返回步骤3，如果所述偏微分方程系统趋于平衡，继续到步骤9；

[0018] 步骤9：针对最终更新的网格节点进行Delaunay非结构三角网格剖分，获得最终的非结构三角网格单元，根据所述规则网格速度模型二维最相邻插值获得所述最终的非结构三角网格单元的物性参数，输出所述网格单元的网格节点坐标、网格单元编号及物性参数。

[0019] 优选地，在所述步骤1中，根据所述规则网格速度模型的速度分布以及地震波有限元法模拟中的空间采样需求，定义网格密度函数

[0020] 优选地，在所述步骤1中，所述网格密度函数 其中vs为横波速度，fmax为在所述地震波有限元法中所用子波的最大频率，N表示一个最短地震波波长内所需的采样点数。

[0021] 优选地，在所述步骤4中，根据所述骨架的拓扑结构定义非结构三角网格中每个节点所受到的骨架臂的力F(p)，所述力F(p)包括所述非结构三角网格中每个节点在水平和垂直方向上的受力：

[0022]

[0023] 其中，Fint表示来自所述骨架臂的内部力函数f(l,l0)，Fext表示来自边界的外力，F(p)的第一列为x分量，第二列为z分量。

[0024] 优选地，所述内部力函数f(l,l0)根据以下公式计算：

[0025]

[0026] 其中l为所述骨架臂的实际长度，l0为期望的非结构三角网格长度，k为弹性系数。

[0027] 优选地，所述非结构三角网格节点在所述节点合力的作用下发生移动。

[0028] 优选地，在所述步骤6中，所述偏微分方程系统为：

[0029]

[0030] 所述偏微分方程系统的初始条件为p(0)＝p0。

[0031] 优选地，在所述步骤7中，所述非结构三角网格节点坐标的更新关系式为：pn+1＝pn+Δt·(F+imgd·imgF)，其中pn+1为n+1时刻的节点坐标，pn为n时刻的节点坐标，Δt为更新的时间步长，imgd为调节参数。

[0032] 优选地，在所述步骤7中，通过向前欧拉时间离散格式离散化所述偏微分方程系统。

[0033] 本公开的优点是针对复杂地表及复杂构造模型的有限元地震波数值模拟，基于Distmesh三角网格剖分算法，提出了一种规则网格速度模型自适应非结构三角网格化方法，从而为复杂地区地震波高精度有限元数值模拟提供了高质量的非结构网格模型。

[0034] 与传统上需要将规则网格速度模型转换为数字图像，然后进行图像预处理，勾勒出模型内部边界线，最后对几何模型进行网格剖分的方法相比，本公开的方法基于Distmesh网格算法实现速度模型的自适应网格剖分，使网格密度分布随模型中速度大小的分布变化，网格单元质量接近理想。此外，本公开的方法引入速度模型的二阶梯度场信息，在二阶梯度场从两侧指向物性分界面，可以将物性分界面附近的网格节点映射到物性分界面上，进而达到网格单元边界与物性分界面贴合的目的。本公开的方法实现了规则网格速度模型的自适应非结构三角网格化，克服了传统方法需要人工提取物性分界面、耗时耗力的缺点，能够为有限元地震波数值模拟算法，特别是间断Galerkin等一类高阶算法提供良好的非结构三角网格模型。

附图说明

[0035] 通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显。

[0036] 图1显示根据示例性实施例的规则网格速度模型的自适应非结构三角网格化方法的流程图。

[0037] 图2显示根据示例性实施例的二维规则网格速度模型及其网格密度要求。

[0038] 图3(a)和图3(b)分别显示根据示例性实施例的规则网格速度模型的二阶梯度场分布及其局部放大图。

[0039] 图4显示根据示例性实施例的在物性分界面两侧的二阶梯度场分布函数。

[0040] 图5显示未使用根据示例性实施例的二阶梯度场映射时算法生成的非结构三角网格。

[0041] 图6显示根据示例性实施例的二阶梯度场映射生成的非结构三角网格单元。

[0042] 图7显示示例性实施例中的Marmousi模型。

[0043] 图8和图9分别显示通过根据示例性实施例的规则网格速度模型的自适应非结构三角网格化方法，对图7的Marmousi模型进行网格剖分的结果及其局部放大图。

具体实施方式

[0044] 下面将参照附图更详细地描述本公开的优选实施方式。虽然附图中显示了本公开的优选实施方式，然而应该理解，可以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

[0045] 本公开提供一种规则网格速度模型的自适应非结构三角网格化方法，包括以下步骤：

[0046] 步骤1：根据所述规则网格速度模型的速度分布以及地震波有限元法模拟中的空间采样需求，定义网格密度函数

[0047] 步骤2：根据所述网格密度函数 在所述规则网格速度模型区域内随机分布初始的网格节点；

[0048] 步骤3：对所述网格节点进行Delaunay非结构三角网格剖分，并将所述非结构三角网格的节点坐标定义为二维数组变量p＝[x z]；

[0049] 步骤4：将剖分的非结构三角网格类比于骨架结构，骨架结构即为非结构三角网格，根据所述骨架的拓扑结构定义所述非结构三角网格中每个节点所受到的骨架臂的力F(p)，F代表臂力场；

[0050] 步骤5：对所述规则网格速度模型进行二阶求导，获得所述规则网格速度模型的二阶梯度场，根据二阶梯度场设置二阶梯度力场imgF，二阶梯度力场imgF正比于二阶梯度场，其中二阶梯度力场imgF定义在所述规则网格上，通过二维最相邻插值获得非结构三角网格的每个网格节点受到的二阶梯度场力imgF(p)；

[0051] 步骤6：引入时间变量，构建所述非结构三角网格节点坐标和节点合力关于时间的偏微分方程系统，其中所述臂力场F和所述二阶梯度力场imgF共同构成了所有非结构三角网格节点的合力场，所述骨架臂的力F(p)和所述二阶梯度场力imgF(p)构成所述非结构三角网格节点受到的节点合力，非结构三角网格节点在所述合力的作用下发生移动；

[0052] 步骤7：离散化所述偏微分方程系统，获得所述非结构三角网格节点坐标的更新关系式pn+1＝pn+Δt·(F+imgF)，其中pn+1为n+1时刻的节点坐标，pn为n时刻的节点坐标，Δt为更新的时间步长，imgd为调节参数，以便根据现在时刻的节点坐标及节点合力，计算下一个时刻的节点坐标，获得更新的网格节点；

[0053] 步骤8：判断所述偏微分方程系统是否趋于平衡，如果所述偏微分方程系统还没有趋于平衡，返回步骤3，针对更新的网格节点，重新进行Delaunay非结构三角网格剖分，计算非结构三角网格节点的节点合力，计算下一个更新的节点坐标；如果所述偏微分方程系统趋于平衡，继续到步骤9；

[0054] 步骤9：针对最终更新的网格节点进行Delaunay非结构三角网格剖分，获得最终的非结构三角网格单元，根据所述规则网格速度模型二维最相邻插值获得所述最终的非结构三角网格单元的物性参数(例如速度、密度等)，输出所述网格单元的网格节点坐标、网格单元编号及物性参数。

[0055] 优选地，所述网格密度函数 其中vs为横波速度，fmax为在所述地震波有限元法中所用子波的最大频率，N表示一个最短地震波波长内所需的采样点数。

[0056] 优选地，所述力F(p)包括所述非结构三角网格的节点在水平和垂直方向上的受力：

[0057]

[0058] 其中，Fint表示来自所述骨架臂的内部力函数f(l,l0)，Fext表示来自边界的外力，F(p)的第一列为x分量，第二列为在z分量。

[0059] 优选地，所述内部力函数f(l,l0)根据以下公式计算：

[0060]

[0061] 其中l为所述骨架臂的实际长度，l0为期望的非结构三角网格长度，k为弹性系数。

[0062] 优选地，所述偏微分方程系统为：

[0063]

[0064] 所述偏微分方程系统的初始条件为p(0)＝p0。

[0065] 优选地，通过向前欧拉时间离散格式离散化所述偏微分方程系统。

[0066] 下面参考示例性实施例描述根据本公开的的规则网格速度模型的自适应非结构三角网格化方法，所述方法包括以下步骤(见图1)：

[0067] 步骤1：根据规则网格速度模型的速度分布以及地震波有限元法模拟中的空间采样需求，定义网格密度函数 其中，vs为横波速度，fmax为地震波有限元法模拟中所用子波的最大频率，N表示一个最短地震波波长内所需的采样点数。图2显示了根据示例性实施例的规则网格速度模型及其网格密度要求。

[0068] 步骤2：根据定义的网格密度函数 在规则网格速度模型区域内随机分布初始的网格节点。

[0069] 步骤3：对初始的网格节点进行Delaunay非结构三角网格剖分，并将非结构三角网格的节点坐标定义为二维数组变量p＝[x z](其中x，z分别表示水平方向和垂直方向坐标)。

[0070] 步骤4：将剖分的非结构三角网格类比于骨架结构。根据骨架结构的(即非结构三角网格)拓扑结构定义骨架中每个骨架节点(即网格节点)所受到的骨架臂的力F(p)，F代表臂力场。力F(p)包括骨架节点(即网格节点)在水平和垂直方向上的受力：

[0071]

[0072] 其中，Fint表示来自骨架臂的内部力函数f(l,l0)，Fext表示来自边界的外力，F的第一列为x分量，第二列为z分量。内部力函数f(l,l0)可以模拟线性弹簧的作用力，且只考虑排斥力，而没有吸引力。内部力函数f(l,l0)根据以下公式计算：

[0073]

[0074] 其中l为骨架臂的实际长度，l0为期望的非结构三角网格长度，k为弹性系数。

[0075] 步骤5：对规则网格速度模型进行二阶求导，获得其二阶梯度场，并根据速度模型的二阶梯度场设置相应的二阶梯度力场imgF，在实施例中，二阶梯度力场imgF正比于二阶梯度场。

[0076] 二阶梯度场在速度均匀的地方二阶梯度为零，在速度变化的界面处两侧会形成近似垂直指向界面的矢量，如图3(a)和图3(b)所示。二阶梯度场具有良好的特性，某处界面两侧的二阶梯度场分布函数，在界面上二阶梯度接近于零，在界面附近二阶梯度达到最大，远离界面梯度变为零，如图4所示。二阶梯度力场定义在规则网格上，通过二维最相邻插值获得非结构三角网格的每个网格节点受到的二阶梯度场力imgF(p)。

[0077] 步骤6：引入时间变量t，构建非结构三角网格节点坐标和节点合力关于时间的偏微分方程系统：

[0078]

[0079] 其中方程系统的初始条件为p(0)＝p0。臂力场F和二阶梯度力场imgF共同构成了所有骨架(网格)节点的合力场，骨架臂的力F(p)和二阶梯度场力imgF(p)构成非结构三角网格节点受到的节点合力，节点在合力的作用下发生移动。

[0080] 步骤7：离散化步骤6中的偏微分方程系统，上面的方程系统可以采用简单的向前欧拉时间离散格式。在tn＝nΔt时刻，逼近解pn≈p(tn)。通过下式，根据现在时刻网格节点坐标及节点合力求取下一个时刻的网格节点坐标：

[0081] pn+1＝pn+Δt·(F+imgd·imgF)

[0082] 其中pn+1为n+1时刻的节点坐标，pn为n时刻的节点坐标，Δt为更新的时间步长，imgd为调节参数，它使算法在网格质量与物性分界面贴合程度之间作出折中。

[0083] 如图5所示，当不采用二阶梯度场对网格节点进行映射时，网格并没有对模型内界面进行贴体的剖分，界面处网格凹凸不平。加入二阶梯度场信息后基本上实现了网格沿界面的贴体剖分，特别是对于正弦函数界面以及楔形体界面处的完美贴合充分说明了根据示例性实施例的方法的有效性，见图6。

[0084] 步骤8：判断所述偏微分方程系统是否趋于平衡，如果所述偏微分方程系统还没有趋于平衡，返回步骤3，针对更新的网格节点，重新进行Delaunay非结构三角网格剖分，计算非结构三角网格节点受到的骨架臂的力与二阶梯度场力的合力，计算下一个更新的节点坐标；如果所述偏微分方程系统趋于平衡，继续到步骤9。当偏微分方程系统趋于平衡时满足F(p)＝0，网格节点不再移动。此时网格趋于最佳，同时相应网格单元边界基本上与速度模型中的物性分界面贴合。

[0085] 步骤9：针对最终更新的网格节点进行Delaunay非结构三角网格剖分，获得最终的非结构三角网格单元，根据所述规则网格速度模型二维最相邻插值获得所述最终的非结构三角网格单元的物性参数，输出所述网格单元的网格节点坐标、网格单元编号及物性参数。

[0086] 图7显示示例性实施例中使用的Marmousi模型，图8和图9分别显示通过根据示例性实施例的规则网格速度模型的自适应非结构三角网格化方法，对图7的Marmousi模型进行网格剖分的结果及其局部放大图。从图8和图9可见，整体上，网格质量良好，网格对于重要界面勾勒清新，实现了沿重要物性界面进行贴体剖分的要求，这一点从局部放大图上看得更加清晰。图9中右上角的网格密度明显要比其他部分大，这也充分显示了非结构网格能够根据地下介质的分布自适应决定网格密度的特点。

[0087] 本公开的方法应用于地震波有限元数值模拟前处理的网格建模过程。本公开的方法基于Distmesh非结构三角网格剖分算法，将网格生成过程类比为骨架平衡的过程，生成网格大小随速度自适应分布的高质量网格。此外，该方法考虑速度模型的二阶梯度场的信息，将网格节点映射到规则网格模型定义的隐式物性分界面上，使非结构网格单元的边与速度模型中的物性分界面自适应贴合，实现了规则网格模型向非结构网格的准确转化，达到保存重要物性分界面的目的，为有限元的数值模拟方法提供可靠的非结构网格模型。

[0088] 以上已经描述了本公开的实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各施例。在不偏离所说明的实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释实施例的原理和实际应用，或者使本技术领域的其它普通技术人员能理解本文披露的实施例。