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

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

【技术实现步骤摘要】

本公开涉及有限元地震波数值模拟方法，特别涉及一种规则网格速度模型的自适应非结构三角网格化方法。
技术介绍
目前，在我国南方和西部地区，油气地震勘探的重点正转向丘陵、山前构造带等复杂地区。这些地区地表条件异常复杂，地形起伏剧烈，高差变化非常大，岩性速度变化大导致近地表结构不均匀性严重，同时地下构造复杂，如褶皱强烈、断层发育、构造陡峻、地层变化大等。这些特点导致这些地区地震资料信噪比低及静校正难等各方面的问题。从根本上解决这些勘探问题，需要对起伏地表条件下地震波的传播规律和波场特征进行深入的认识，而有限元法是进行复杂地表和复杂构造地区地震波数值模拟最有效的技术手段。在进行地震波数值模拟时，首先需要对地球物理-地质模型(包括纵横波速度、密度等)进行空间上的网格化离散。有限差分法一般使用规则网格，而有限元法则多使用非结构网格。目前在地震勘探领域中，速度模型大都是定义在规则网格上，以二维矩阵数据的形式给出，其主要原因是有限差分法仍然是目前使用最广泛的地震波数值模拟技术，且规则网格的数据结构简单，便于广泛地传播使用。但规则网格定义的速度模型存在着模型的过采样问题及缺乏显式的边界等问题。非结构网格是指没有规则拓扑关系的网格。它能够对模型进行贴体的网格剖分。特别地，三角形网格能够模拟任意复杂的地形起伏及复杂的地下构造，因此能够更好地刻画模型中的界面信息。此外，采用非结构网格能够适应地下介质的该物性参数变化，自适应的调整网格大小和密度，而不会造成过采样的问题。因此，对于具有复杂结构的模型，非结构网格拥有更强的适应能力。由于目前速度模型基本上是定义在规则网格上，而有限元法使用非结构网格，因此需要考虑如何将规则网格上定义的模型三角网格化。此外，采用间断Galerkin等高阶有限元数值模拟方法时，为了保证模拟精度和计算效率，对使用的非结构网格剖分质量也有较高的需求。对于速度等地球物理-地质模型，最感兴趣的是地下介质物性参数突变的地方。在进行地震波数值模拟时，需要充分尊重这些界面信息的存在。当采用常规思路进行网格剖分时，需要将规则网格定义速度模型转换为数字图像，然后进行图像预处理，勾勒出模型内部边界线(segment)，即物理模型向几何模型的转化过程，最后对几何模型进行网格剖分。对于简单的块状模型，因为其界面清晰，区域划分容易，采用图像处理的思路不失为一种好的方法。但是，当地下模型复杂，界面信息凌乱时，使用图像处理进行区域划分比较麻烦。并且，专利技术人发现，常规思路的几个步骤是相对独立的，无法形成一个独立的、自适应的网格剖分工具。网格的生成通常被认为是一项复杂的任务，因此网格生成器也总是以“黑匣子”的形式独立于有限元代码的编写环节。Distmesh是由Per-OlofPersson(2005)提出的一种非常简单的三角形网格生成技术。区别于传统的基于几何概念的网格生成方法，Distmesh算法基于三角形网格与钢骨架结构的类比，以及钢骨架结构的整体平衡性，能够提供高质量的网格剖分结果，同时能够物性分布自适应地改变网格密度，但Distmesh算法难以保证剖分的非结构网格与物性分界面的自适应贴体化。基于上述原因，期待一种适用于山前带等复杂(地表/地下)地区的规则网格速度模型的自适应非结构三角网格化方法。
技术实现思路
本公开的目的是建立一种适用于复杂地区的能够有效匹配规则网格模型的非结构三角网格化方法。本公开所采用的解决方案如下：一种规则网格速度模型的自适应非结构三角网格化方法，包括以下步骤：步骤1：定义网格密度函数步骤2：根据所述网格密度函数在所述规则网格速度模型区域内随机分布初始的网格节点；步骤3：对所述网格节点进行Delaunay非结构三角网格剖分，并将非结构三角网格的节点坐标定义为二维数组变量p＝[xz]；步骤4：将所述非结构三角网格类比于骨架结构，定义所述非结构三角网格中每个节点所受到的骨架臂的力F(p)；步骤5：对所述规则网格速度模型进行二阶求导，获得所述规则网格速度模型的二阶梯度场，根据所述二阶梯度场设置二阶梯度力场imgF，所述二阶梯度力场imgF正比于所述二阶梯度场，然后通过二维最相邻插值获得所述非结构三角网格的每个网格节点受到的二阶梯度场力imgF(p)；步骤6：构建所述非结构三角网格节点坐标和节点合力关于时间的偏微分方程系统，其中所述节点合力由所述骨架臂的力F(p)和所述二阶梯度场力imgF(p)构成；步骤7：离散化所述偏微分方程系统，获得所述非结构三角网格节点坐标的更新关系式，根据现在时刻的节点坐标及节点合力，计算下一个时刻的节点坐标，获得更新的网格节点；步骤8：判断所述偏微分方程系统是否趋于平衡，如果所述偏微分方程系统还没有趋于平衡，返回步骤3，如果所述偏微分方程系统趋于平衡，继续到步骤9；步骤9：针对最终更新的网格节点进行Delaunay非结构三角网格剖分，获得最终的非结构三角网格单元，根据所述规则网格速度模型二维最相邻插值获得所述最终的非结构三角网格单元的物性参数，输出所述网格单元的网格节点坐标、网格单元编号及物性参数。优选地，在所述步骤1中，根据所述规则网格速度模型的速度分布以及地震波有限元法模拟中的空间采样需求，定义网格密度函数优选地，在所述步骤1中，所述网格密度函数其中vs为横波速度，fmax为在所述地震波有限元法中所用子波的最大频率，N表示一个最短地震波波长内所需的采样点数。优选地，在所述步骤4中，根据所述骨架的拓扑结构定义非结构三角网格中每个节点所受到的骨架臂的力F(p)，所述力F(p)包括所述非结构三角网格中每个节点在水平和垂直方向上的受力：其中，Fint表示来自所述骨架臂的内部力函数f(l,l0)，Fext表示来自边界的外力，F(p)的第一列为x分量，第二列为z分量。优选地，所述内部力函数f(l,l0)根据以下公式计算：其中l为所述骨架臂的实际长度，l0为期望的非结构三角网格长度，k为弹性系数。优选地，所述非结构三角网格节点在所述节点合力的作用下发生移动。优选地，在所述步骤6中，所述偏微分方程系统为：所述偏微分方程系统的初始条件为p(0)＝p0。优选地，在所述步骤7中，所述非结构三角网格节点坐标的更新关系式为：pn+1＝pn+Δt·(F+imgd·imgF)，其中pn+1为n+1时刻的节点坐标，pn为n时刻的节点坐标，Δt为更新的时间步长，imgd为调节参数。优选地，在所述步骤7中，通过向前欧拉时间离散格式离散化所述偏微分方程系统。本公开的优点是针对复杂地表及复杂构造模型的有限元地震波数值模拟，基于Distmesh三角网格剖分算法，提出了一种规则网格速度模型自适应非结构三角网格化方法，从而为复杂地区地震波高精度有限元数值模拟提供了高质量的非结构网格模型。与传统上需要将规则网格速度模型转换为数字图像，然后进行图像预处理，勾勒出模型内部边界线，最后对几何模型进行网格剖分的方法相比，本公开的方法基于Distmesh网格算法实现速度模型的自适应网格剖分，使网格密度分布随模型中速度大小的分布变化，网格单元质量接近理想。此外，本公开的方法引入速度模型的二阶梯度场信息，在二阶梯度场从两侧指向物性分界面，可以将物性分界面附近的网格节点映射到物性分界面上，进而达到网格本文档来自技高网...

【技术保护点】
一种规则网格速度模型的自适应非结构三角网格化方法，包括以下步骤：步骤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非结构三角网格剖分，获得最终的非结构三角网格单元，根据所述规则网格速度模型二维最相邻插值获得所述最终的非结构三角网格单元的物性参数，输出所述网格单元的网格节点坐标、网格单元编号及物性参数。...

【技术特征摘要】
1.一种规则网格速度模型的自适应非结构三角网格化方法，包括以下步骤：步骤1：定义网格密度函数步骤2：根据所述网格密度函数在所述规则网格速度模型区域内随机分布初始的网格节点；步骤3：对所述网格节点进行Delaunay非结构三角网格剖分，并将非结构三角网格的节点坐标定义为二维数组变量p＝[xz]；步骤4：将所述非结构三角网格类比于骨架结构，定义所述非结构三角网格中每个节点所受到的骨架臂的力F(p)；步骤5：对所述规则网格速度模型进行二阶求导，获得所述规则网格速度模型的二阶梯度场，根据所述二阶梯度场设置二阶梯度力场imgF，所述二阶梯度力场imgF正比于所述二阶梯度场，然后通过二维最相邻插值获得所述非结构三角网格的每个网格节点受到的二阶梯度场力imgF(p)；步骤6：构建所述非结构三角网格节点坐标和节点合力关于时间的偏微分方程系统，其中所述节点合力由所述骨架臂的力F(p)和所述二阶梯度场力imgF(p)构成；步骤7：离散化所述偏微分方程系统，获得所述非结构三角网格节点坐标的更新关系式，根据现在时刻的节点坐标及节点合力，计算下一个时刻的节点坐标，获得更新的网格节点；步骤8：判断所述偏微分方程系统是否趋于平衡，如果所述偏微分方程系统还没有趋于平衡，返回步骤3，如果所述偏微分方程系统趋于平衡，继续到步骤9；步骤9：针对最终更新的网格节点进行Delaunay非结构三角网格剖分，获得最终的非结构三角网格单元，根据所述规则网格速度模型二维最相邻插值获得所述最终的非结构三角网格单元的物性参数，输出所述网格单元的网格节点
\t坐标、网格单元编号及物性参数。2.根据权利要求1所述的规则网格速度模型的自适应非结构三角网格化方法，其中，在所述步骤1中，根据所述规则网格速度模型的速度分布以及地震波有限元法模拟中的空间采样需求，定义网格密度函数3.根据权利要求1所述的...

【专利技术属性】
技术研发人员：薛昭，佘德平，杨丽，
申请(专利权)人：中国石油化工股份有限公司，中国石油化工股份有限公司石油物探技术研究院，
类型：发明
国别省市：北京;11