CN110335357B

CN110335357B - 一种约束曲面多分辨率控制预处理方法

Info

Publication number: CN110335357B
Application number: CN201910433643.8A
Authority: CN
Inventors: 李兆亮; 熊盛青; 刘培刚; 刘钰洋
Original assignee: China University of Petroleum East China; China Aero Geophysical Survey and Remote Sensing Center for Natural Resources
Current assignee: China University of Petroleum East China; China Aero Geophysical Survey and Remote Sensing Center for Natural Resources
Priority date: 2019-05-23
Filing date: 2019-05-23
Publication date: 2023-01-24
Anticipated expiration: 2039-05-23
Also published as: CN110335357A

Abstract

本发明提供了一种约束曲面多分辨率控制预处理方法，包括网格尺寸计算，边界线重构，边界面重构，边界面网格优化，本发明研究重构边界对储层表面的逼近程度，采用能量泛函对重构过程进行控制，提出简洁的控制函数；以角点、边界线、约束面的层级递进顺序计算各级别元素的网格尺寸；根据网格尺寸控制实现边界线和约束面的重构。本发明采用模型结构改进前沿推进法实现面内布点和三角剖分；采用Voronoi图保护区实现复杂限定条件的处理；并且研究表面三角网格优化方法，并对Laplacian算法进行改进，实现空间曲面的光顺。本发明提高了约束面的质量，进而提高了网格的质量。

Description

一种约束曲面多分辨率控制预处理方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种约束曲面多分辨率控制预处理方法。

背景技术

网格的质量直接影响到有限元分析、数值模拟等实际应用的精确度、数值的稳定性和正确性，而在建模过程中网格之间的交切等操作造成网格质量逐渐变差，将影响数值的稳定性和建模的整个流程。同时，受计算机软硬件的限制，有限元分析及油藏数值模拟等数值计算对计算网格的数量有较大的限制。因此需要优化四面体网格的质量，同时控制四面体网格的数量。

在建模工作中，约束面控制着储层网格剖分质量和效果。首先，储层网格化结果需要尽可能拟合约束面，且需要尽量保留各个尺度的地质特征；其次，网格化结构受到约束面本身质量的影响；最后，约束面的质量会影响剖分算法。因此，在储层网格化之前，首先要保证约束面的质量。要实现对四面体网格的质量和数量控制，必然需要预先对约束曲面进行优化和多分辨率的控制。

如果使用的约束数据来自于构造模型的表面以及各种界面。一般需对这些数据进行拓扑处理以保证其几何一致性和拓扑一致性，但这些界面的质量可能较差，主要存在以下几种情况：

(1)约束面本身的剖分质量差，其构成三角形中狭长或扁平面片多。这与控制点的数量和分布有关，也与三角化的策略和技术有关。

(2)约束面本分的分辨率不符合要求。比如在自适应三角剖分中，一般要求变化平缓区域中三角形面积大、分辨率低，而变化剧烈区域中三角形面积小、分辨率高。但是储层地质模型中在平缓区域的三角形不能无限大，应该符合一定的预期尺寸。

(3)不同约束面之间的分辨率不一致。比如构造模型中，当不同地层面来源数据密度不一致时，经剖分后表面的分辨率可能不一致。同时，地层面与断层面之间分辨率一般不一致；地层面、断层面与裂缝面之间分辨率一般也不一致。

此外，对于要求变分辨率的情况，如多种岩体、多个沉积相、以及需关注的区域如高渗透率砂体等，要求不同分辨率，也需要特殊处理。

对于约束表面的重构属于网格重剖分技术。目前关于网格重剖分技术的研究较多，八叉树法、前沿推进法、Delaunay法等均可用于网格重剖分。这些方法的基本思想基本一致：(1)修改控制边界，比如修改边界线上的点以改变边界分段；(2)修改面内的顶点分布，使其符合分辨率要求，点之间的距离一般采用目标函数控制；(3)对所有点进行约束剖分。需要注意的是，为了保持模型的特征，控制点是不可以变动的。

发明内容

针对现有技术存在的问题，本发明提供了一种约束曲面多分辨率控制预处理方法。

为实现上述发明目的，本发明采用以下技术方案予以实现：

一种约束曲面多分辨率控制预处理方法，所述方法包括以下步骤：

S1：网格尺寸计算：计算角点网格尺寸、边界线网格尺寸、约束面网格尺寸；

S2：边界线重构：采用从角点开始的传播方法，根据网格尺寸添加点，对边界线进行重构；

S3：边界面重构：在原约束面上布置顶点，将这些顶点使用改进的前沿推进算法进行三角网格剖分，生成新的约束面；

S4：边界面网格优化：使用改进的Laplacian光顺法进行平面网格优化，再使用改进的Laplacian光顺法进行约束面网格优化。

进一步的，所述S1具体为：以角点、边界线、约束面的层级递进顺序计算角点网格尺寸、边界线网格尺寸、约束面网格尺寸；所述边界线网格尺寸的计算采用从角点开始的尺寸传播算法；所述约束面网格尺寸的计算采用从边界线开始的尺寸传播算法。

进一步的，所述边界线重构还包括边界线逼近控制的处理方法：在预处理过程中，对于重构边界线中新生成的每一条线段，通过边界线距离控制函数检测原线段的端点到线段的最大距离，如果最大距离大于一定阈值，则调整线段端点重新构造线段并检测，直到新线段满足要求为止；

所述边界线距离控制函数定义为：

f(E_ab)＝Max(||P_iE_ab||)＜ε，P_i∈L，P_i∈E_ab

其中，L表示原边界线，P_a和P_b表示新生成的线段端点，P_i表示原边界线的线段端点，E_ab表示新生成的线段，P_i’表示P_i在线段E_ab所在直线上的投影，ε表示距离阈值。

进一步的，所述S3中采用曲面前沿推进法在原约束面上布置顶点。

进一步的，所述三角网格剖分的方法为拓扑模型的改进的前沿推进算法：

步骤1：建立点集PS、边集ES、三角形集TS、临时顶点集LS、临时顶点队列LQ；

步骤2：将边界离散后小线段的端点及区域内所有布设的顶点加入PS，以初始化点集PS；

步骤3：将边界离散后生成的所有有向线段加入ES，以初始化边集ES；此时，边集ES就是前沿Ω；

步骤4：取边集ES中的边Ei，设Ei的端点为A和B，则Ei可表示为

以Ei 为基础，即选定Ei作为选定前沿，搜索临近Ei且位于Ei左侧的顶点Pi加入LS，使用右手准则；

步骤5：取LS中所有的点Pi，计算Pi与Ei组成的顶角∠APiB的大小，将Pi按照顶角从大到小排序，使用最大最小角准则存入队列LQ；

步骤6：取LQ中一点Pi，与边Ei组成三角形ABPi；

步骤7：判断新三角形是否满足以下要求：

③边

与已生成的边不能相交；

④新三角形面积的必须大于零；

如果不满足要求转步骤6；如果满足要求转步骤8；

步骤8：将新顶点与选定前沿端点连成的有向边

加入到边集ES中，生成前沿Ω的新元素，并将新三角形ABPi加入到三角形集TS中，删除原选定的前沿Ei；

步骤9：如果边集ES为空，转步骤10；如果边集ES不为空；转步骤4；

步骤10：算法结束。

进一步的，所述S3还包括：使用改进的前沿推进算法进行三角网格剖分结束后，根据拓扑模型，即对应边关系，合并半边将其转化为普通三角网格。

进一步的，所述S3中对于复杂限定条件的区域剖分，具体步骤包括：

(1)利用Voronoi图建立限定条件的保护区，小角度问题在保护区内解决；

(2)保护区以外的区域按照普通过程进行剖分；

(3)两部分区域合并作为最终的网格剖分结果。

进一步的，所述步骤S3还包括约束面逼近控制的处理方法：对于重构曲面中新生成的每一个面片，根据约束面的距离控制函数检测原约束面的面片顶点到面片的最大距离，如果最大距离大于一定阈值，则调整面片顶点位置，构造新面片并检测，直到新面片满足要求为止；

所述约束面的距离控制函数定义为：

f(T_abc)＝Max(||P_iT_abc||)＜ε，P_i∈S，P_i∈T_abc

其中，S表示原约束面，P_a、P_b、P_c表示新生成的面片顶点，P_i表示原约束面的顶点，T_abc表示新生成的面片，P_i’表示P_i在三角形T_abc平面内的投影，ε表示距离阈值。

进一步的，所述S4中所述改进的Laplacian算法：

步骤1：建立待优化顶点数组Q1、Q2，建立顶点优化数据E，优化计数器Count，以约束面所有顶点初始化数组Q1；

步骤2：取Q1中一个顶点p(xi,yi，zi)，

步骤2.1：如果点p为控制点，则点p存入数组Q2对应位置，转步骤2；否则转步骤2.2；

步骤2.2：自拓扑模型TM中查找所有以p为起点的有向边，这些有向边的终点即为顶点p的邻接点；

步骤3：计算所有邻接点的平均坐标点q(xj，yj，zj)，将q存入数组Q2对应位置，作为优化后的p点；

步骤4：如果Q1为空，则转步骤5；否则转步骤2；

步骤5：计算每个顶点的优化能量，即优化前与优化后的距离，e＝f(p，q)＝||p-q||，存入数组E对应位置；

步骤6：交换数组Q1和Q2的元素；

步骤7：取Q1中一个顶点p(xi，yi，zi)，

步骤7.1：如果顶点的优化能量e＜ε，其中ε表示距离阈值，则点p存入数组Q2 对应位置，Count加1，转步骤7；否则转步骤7.2；

步骤7.2：自拓扑模型TM中查找所有以p为起点的有向边，这些有向边的终点即为顶点p的邻接点；

步骤8：计算所有邻接点的平均坐标点q(xj，yj，zj，)，将q存入数组Q2对应位置，作为优化后的p点；

步骤9：如果Q1为空，则转步骤10；否则步骤7；

步骤10：如果Count＝顶点数，转步骤11；否则转步骤5；

步骤11：输出结果网格，算法结束。

进一步的，所述S4还包括针对曲面的改进方法：在曲面条件下，使用顶点p的邻接点计算的平均坐标点q不在原曲面上，q与邻接点构成的三角面片与原曲面之间产生一定距离，且此距离的最大值是点p与点q分别所处的平行平面之间的距离，称为“最小距离”，当点p与点q的“最小距离”超过一定限度时，将其平行移动到与点q垂直的位置。

与现有技术相比，本发明的优点和有益效果是：

(1)本发明研究重构边界对储层表面的逼近程度，采用能量泛函(距离)对重构过程进行控制，提出简洁的控制函数；

(2)以角点、边界线、约束面的层级递进顺序计算各级别元素的网格尺寸；根据网格尺寸控制实现边界线和约束面的重构。

(3)本发明采用模型结构改进前沿推进法实现面内布点和三角剖分，本算法中采用半边进行三角网格剖分，每条半边只用一次，而且具有方向限定，简化了推进过程中的检测方法及判定数量；采用Voronoi图保护区实现复杂限定条件的处理；

(4)并且研究表面三角网格优化方法，并对Laplacian算法进行改进，实现空间曲面的光顺。

附图说明

图1为本发明的约束曲面多分辨率控制预处理的流程示意图；

图2为角点的网格尺寸示意图；

图3为角点处边界线夹角比较小时边界线上的网格尺寸计算；

图4为前沿推进法布点示意图；

图5为限定条件的保护区示例，Branets et al.，2009；

图6为Voronoi保护区及三角剖分a；

图7为Voronoi保护区及三角剖分b；

图8为复杂限定条件保护区a；

图9为复杂限定条件保护区b；

图10为不同半径时保护圆相交情况；

图11为空间中曲线逼近。

具体实施方式

下面结合实施例对本发明作进一步详细的描述，但发明的实施方式不限于此。

实施例1、约束曲面的逼近程度控制

离散构造模型是一个理想的分段逼近连续对象，离散模型比所有的接近连续对象更准确。离散表面的精确度是由点的精确度和点的密度来决定的。实际上，构造模型应具有最小分辨率从而能够反映出期望几何结构的复杂性。但是，对几何复杂程度的理解与数据特征相关。在维持逻辑边界的定义下，表面的分辨率可以修改。

因此，模型的分辨率应至少达到如下程度：模型和数据之间的不匹配是在数据的不确定性范围之内；考虑到解释输入、模拟推理以及确保网格质量的需要，模型的分辨率也可以更高。在许多情况下，地质表面实行可变分辨率或自适应分辨率有用而且合适。因此，在对原始的构造模型进行质量控制和不同分辨率控制时仍然需要保持对原始曲面的合理逼近。本发明提出了在离散的构造模型质量控制和不同分辨率控制时仍然保持对原始曲面的合理逼近的控制函数。

Zhao et al.(2000)为了表示空间曲面Γ与离散点集PS的逼近程度，定义了曲面能量泛函数：

其中，d(X)＝distance(X，S)为曲面上任意一点X到离散点集PS的距离，ds为曲面上点X邻域所在的微小面元。

本发明中，作为输入的约束曲面本身是离散的，使用三角网表示。因此，对于约束曲面的逼近即重构三角网格对原三角网格的拟合问题。约束曲面的重构需要解决网格尺寸相近与形态一致的问题。

重构的约束曲面与原曲面形态不一致的情形是由于两曲面相应元素之间距离过大导致的。本发明根据离散型构造曲面的特点，分别解决边界线和边界面的逼近问题。

(1)边界线逼近控制

在预处理过程中，对于重构边界线中新生成的每一条线段，检测线段端点间原边界线上各个点到线段的最大距离，如果最大距离大于一定阈值，则调整线段端点重新构造线段并检测，直到新线段满足要求为止。与新线段距离最大的点一定是原边界线的线段端点，因此，只要检测原线段的端点即可。

设L表示原边界线，P_a和P_b表示新生成的线段端点，P_i表示原边界线的线段端点，E_ab表示新生成的线段，P_i’表示P_i在线段E_ab所在直线上的投影，ε表示距离阈值，则边界线距离控制函数定义为：

f(E_ab)＝Max(||P_iE_ab||)＜ε，P_i∈L，P_i∈E_ab

(2)约束曲面逼近控制

在预处理过程中，对于重构曲面中新生成的每一个面片，检测面片范围内原曲面上的各点到面片的最大距离，如果最大距离大于一定阈值，则调整面片顶点位置，构造新面片并检测，直到新面片满足要求为止。距离面片最大的点一定是原约束面的面片顶点，因此，只要检测原约束面的面片顶点即可。

设S表示原约束面，P_a、P_b、P_c表示新生成的面片顶点，P_i表示原约束面的顶点， T_abc表示新生成的面片，P_i’表示P_i在三角形T_abc平面内的投影，ε表示距离阈值，则约束面的距离控制函数定义为：

f(T_abc)＝Max(||P_iT_abc||)＜ε，P_i∈S，P_i∈T_abc

实施例2、约束曲面多分辨率控制预处理方法

对于由多组边界面和边界线构成的构造模型，进行网格重剖分的经典策略是逐一重网格化。

本发明采取逐一重网格化的策略，先重构线性几何要素，即每条边界线及限定条件(线)，然后在曲面空间里重构面几何要素，即每个约束面，不再参数化到平面上。在重构几何要素的过程中，必须计算其元素尺寸，在此基础上进行网格划分。

2.1、网格尺寸计算方法

自适应分辨率网格剖分中，单一网格元素的大小受到多个因素影响，最主要的因素是网格所在区域的形态特征(如曲面的曲率、构造模型的断距、断层线等)以及预期网格大小。对于具体区域而言，通常是以预期网格为前提、根据区域形态计算区域内每个位置的网格尺寸。为网格剖分需要，本发明就复杂约束体域中控制要素，包括角点、边界线(接触线)、约束面，讨论其尺寸控制问题。

根据构造模型的局部特征，以控制要素的尺寸作为基础，使用元素到控制要素的距离函数来计算其他元素尺寸。控制要素可以是角点、边界线、约束面。如果存在多个特征要素，每个特征要素在元素位置都会有一个控制尺寸δ_i，则元素的最终尺寸由所有尺寸的最小值决定，即

δ(x)＝Min{δ₁，δ₂，…，δ_n}

(1)角点网格尺寸

角点是一种控制点，本身没有大小，也不可细分，角点的网格尺寸指的是角点附近区域的网格大小。角点网格大小受多种因素制约，图2是角点的网格尺寸示意图，如图2所示：①如果是独立区域的角点，其网格大小与区域期望网格大小相同；②如果多个区域共享角点，则角点网格大小为多个区域预期网格尺寸的最小值；③如果控制点处几何形态特殊，如夹角很小，则加入其他控制因素，如密度。另外控制点的尺寸可以人工给定。角点尺寸一般小于约束面区域内的期望网格大小，尤其是角度较小时，必须压缩尺寸，以保证网格质量。

设角点处边界线夹角为α，若α＞π/3，则不考虑角度对角点尺寸的影响；若α较小，则需要特殊计算，其方法与边界线的尺寸计算相关。

(2)边界线网格尺寸计算方法

两个面共享的接触线可能由一条或者多条边界线组成，每条边界线由两个端点和连接端点的一个或者多个小线段构成，每条边界线的端点是不可改动的控制点(角点)。边界线上的小线段将是约束面中组成面片的边。边界线网格的大小实际是指以边界线中线段为边的网格的大小，最终转化为网格边的长度，具体到边界线来讲则是边界线上线段的长度。边界线网格大小(即将来的边长)受到区域网格大小、角点附近网格尺寸的限制。边界线的主体部分，通常为边界线中段、占边界线的大部分，网格尺寸与约束面的预期网格尺寸一致。如果边界线为多个约束面共享，则取各约束面期望网格尺寸的最小值。边界线的两端部分受到角点附近网格尺寸的制约。

假设边界线期望网格尺寸和两端角点的尺寸已确定，m表示边界线上的目标网格尺寸(即约束面内的期望网格大小)，令m_a，m_b分别表示两个角点A和B的给定尺寸。在边界线上，期望网格大小在角点A附近时受到m_a控制，在角点B附近时受到m_b控制，而远离两个角点时受到m控制。为达到这一目标需要分别定义两个角点的控制范围。

考虑网格形态，相邻网格的尺寸变化不能太大，控制在α∈[0，1]内。如果m_a≥m，第一段的理论长度为αm_a，第二段的理论长度为α²m_a，依次类推，直到αⁿ m_a≈m；如果m_a≤m，第一段的理论长度为(2-α)m_a，第二段的理论长度为(2-α)²m_a，依次类推，直到(2-α)ⁿ m_a≈m。所以，以m_a≥m为例，定义d_a为角点A的影响距离：

其中n＝|log(m/m_a)/logα|，|·|是绝对值函数。

第i段的长度变化为：

设p为边界线上到角点A距离为d的点，有：

则点p处的网格尺寸为：

每条交线有两个角点，所以选用这两个结果的最小值：min(f(d，m_a，m)，f(d，m_b，m))。

为了控制从角点到角点控制距离的网格尺寸的变化趋势，Botella et al.(2014)引入因子β，修改网格尺寸的变化控制函数。对于每个输入的线段端点，计算相应的网格大小(线段长度)：

其中d是输入线段端点和角点之间的距离。

当角点处边界线夹角比较大时，角点尺寸可以直接采用网格期望尺寸，边界线的尺寸计算方法比较简单，可以采用平均法，也可以从一个角点开始采用期望尺寸计算到另一个角点，或者中部采用期望尺寸而两角点附近略微调整。

接下来讨论特殊情况，角点处边界线夹角比较小时边界线上的网格尺寸计算方法。

设期望网格尺寸为Size，角点处边界线夹角为α，角点对边长度为1，则边界线之间距离变为Size的位置是l＝Size的位置。这时从角点到该位置的尺寸为

且夹角越小r越长。因锐角α不可能改善，因此可以引入其他度量标准，如密度。将网格区域细分，如采用边对分方法，将锐角三角形分解为若干个三角形。图3为角点处边界线夹角比较小时边界线上的网格尺寸计算，如图3所示为一种分解方式，将r一分为二，原三角形分为三个三角形。连续细分直到锐角附近网格密度达到要求为止。由于小角度的存在，此处三角形的质量不会太好。

当角点为多个相邻区域共用时，边界线在角点处的网格尺寸受相邻区域共同控制，取所有以该角点为顶点的边界线夹角的最小值控制的尺寸。

(3)约束面网格尺寸计算方法

在边界线网格尺寸计算过程中，本发明采用了从角点开始的尺寸传播算法。在约束面网格尺寸计算过程中，同样可以采用从边界线开始的尺寸传播算法。得到的网格大小应该与已经计算出来的边界线网格尺寸和给定的约束面期望网格尺寸是一致的。

假设S表示约束面，Cs表示约束面S的所有边界线的集合，p表示约束面S上一个顶点，c表示边界线集合Cs的中一条边界线，p_c表示边界线c上距离p点最近的一个顶点，dc表示顶点p到顶点pc之间的距离，m_c表示pc点的网格尺寸，ms表示约束面的期望网格尺寸。则顶点p的网格尺寸是所有边界线上到顶点p最近的顶点控制的顶点p的网格尺寸的最小值：

m_p＝min{f(d_c，m_c，m_s)}，c∈C_s (3)

其中式(3)中函数f由式(1)或者式(2)给出。

根据边界线网格尺寸和约束面期望网格尺寸计算约束面顶点网格尺寸，(算法1)。

算法1计算约束面网格尺寸

步骤1建立边界线队列LQ并初始化，建立最近顶点队列PQ，建立网格尺寸数组 DQ；

步骤2计算每条边界线中距离顶点p最近的顶点p_c：

步骤2.1从LQ中取一条边界线c；

步骤2.2依次取LQ中所有顶点，计算其到顶点p的距离d_c；

步骤2.3判断距离顶点p最近的顶点p_c存入队列PQ；

步骤2.4如果队列LQ为空，则转步骤3；否则转步骤2.1。

步骤3计算顶点p的网格尺寸：

步骤3.1从PQ中取一个顶点p_c；

步骤3.2根据顶点p_c的网格尺寸及约束面期望网格尺寸计算顶点p的网格尺寸，存入DQ；

步骤3.3若PQ为空则转步骤3.4；否则转步骤3.1；

步骤3.4取DQ中最小值作为顶点p的网格尺寸。

2.2、边界线重构方法

边界线网格尺寸计算完成之后，可以根据计算出来的边界线网格尺寸对边界线进行重构。与边界线网格尺寸计算方法类似，可以采用从角点开始的传播方法，根据网格尺寸添加点，对边界线进行重构。其中新加点处的网格大小等于输入网格中最近的网格尺寸，这也将适用于约束面和约束体域的网格尺寸。

设开始角点A的网格尺寸是m_A，p表示边界线上的一个点，与A的曲线距离为m_A，网格尺寸为m_p。p点的网格大小等于已计算网格尺寸中距离该点最近的网格尺寸，这也将适用于面和体的网格尺寸。则从角点A开始添加的第一个点q位于与角点A曲线距离

的位置。p点仅用于计算要使用的实际网格尺寸。第二个添加点从q点开始计算，以此类推，直到没有更多的点可以添加。根据开始的角点不同，边界线重构结果可能是不同的。

假设边界线的总长度为L，全边界线期望网格尺寸为m，为加速边界线重构，重构边界线可以分解为L/m条线段。如果不能整除则取整为round(L/m)条线段，此时网格边长的曲线距离为

利用1来对边界线进行规则采样，不需要平滑，而且重构之后所有段长度在m/2和3m/2之间。

因为上述添加点过程中采用的是曲线距离，而添加点之后度量使用的是直线距离，如果曲线的曲率较大时，有可能出现新线段不能逼近原曲线的问题，需要进行边界逼近程度控制。如实施例1中所述，在每次添加新顶点之前，计算原边界顶点到新线段的最大距离，如果超出允许范围，则用该顶点作为添加的新顶点，并重新进行边界逼近程度控制，直至达到要求为止，将新顶点插入新边界线。

解决边界线曲率影响问题的一种策略是，在约束条件预处理过程中，将边界线在大曲率位置添加控制点，即通过添加角点的方式将边界线拆分为若干部分，以此来排除由曲率带来的问题。

3、约束面重构方法

约束面的重构是重新划分组成约束面的三角面片的过程。可以分为两个阶段：一是在原约束面上布置顶点，二是在布点之后将这些顶点进行三角剖分，生成新的约束面。重构之后的约束面与原约束面形态尽量一致，且不能丢失小特征。本发明采用前沿推进法布点，采用Delaunay法和前沿推进法两种方法进行三角剖分，并且采用 Voronoi图保护区技术处理约束边界和限定条件中出现的小角度问题。

1)前沿推进法布点

用前沿推进法计算网格顶点以进行布点的思路是，从约束面的边界线开始，向约束面内部进行计算。在每一个计算过程中，如果满足了一定条件，则添加一个点。

布点算法(算法2)采用队列这种数据结构。队列的特点是只能从队列头部取数据、从队列尾部添加数据，即所谓的先进先出。队列保证了先生成的网格前沿先进入下一轮推断，实现了“按层推进”的原则。本算法希望产生等边三角形布点，需要说明的是，此处“等边三角形”表示三角形的三条边尽量一致而不是完全相等。

算法2约束面前沿推进法布点

步骤1建立表示推进前沿的线段队列Q，建立临时交点队列L，建立布点队列P。

步骤2边界线的所有有向线段加入队列Q；

步骤3从Q中取一个有向线段，表示为顶点对(v1，v2)，令v表示(v1，v2)的中点，mv表示约束面中距v最近的网格尺寸。以v为圆心、

为半径、正交于(v1，v2)作圆，则该圆与有向线段(v1，v2)左侧区域约束面的所有交点即可能为我们需要的点，加入队列L。|·|是绝对值函数，||·||表示距离，前沿推进法布点示意图如图4 所示。

步骤4取L中的一个元素q，判断q是否满足条件：一是不太靠近已布点，二是满足距离控制函数。如果满足条件，将q加入队列P；设mq是距q最近的网格尺寸，以q 为中心、半径为1.5×mq作球，球体内的所有点与q构成的有向边插入到队列Q；

步骤5如果L不为空，则转步骤4；

否则转步骤6；

步骤6如果Q不为空，则转步骤3；

否则转步骤7；

步骤7算法结束。

算法中新生成前沿有向边的构造规则为，以与q点构成三角形的对边方向为正(逆时针或者顺时针)，对边顶点与q点构成的反方向(顺时针或者逆时针)有向边。以边

为例，生成的前沿有向边为

和

2)三角网格剖分

在边界内布点结束之后，可以利用采用多种方法实现三角网格的剖分。本发明提出一种新的拓扑模型的改进的前沿推进算法(算法3)。

算法3改进的前沿推进法三角剖分

步骤1建立点集PS、边集ES、三角形集TS、临时顶点集LS、临时顶点队列LQ。

步骤2将边界离散后小线段的端点及区域内所有布设的顶点加入PS，以初始化点集 PS。

步骤3将边界离散后生成的所有有向线段加入ES，以初始化边集ES。此时，边集ES就是前沿Ω。

步骤4取边集ES中的边Ei，设Ei的端点为A和B，则Ei可表示为

以Ei为基础，即选定Ei作为选定前沿，搜索临近Ei且位于Ei左侧的顶点Pi加入LS(右手准则)。

步骤5取LS中所有的点Pi，计算Pi与Ei组成的顶角∠APiB的大小，将Pi按照顶角从大到小排序，存入队列LQ(使用最大最小角准则)。

步骤6取LQ中一点Pi，与边Ei组成三角形ABPi。

步骤7判断新三角形是否满足以下要求：

①边

与已生成的边不能相交；

②新三角形面积的必须大于零。

如果不满足要求转步骤6；如果满足要求转步骤8。

步骤8将新顶点与选定前沿端点连成的有向边

加入到边集ES中，生成前沿Ω的新元素，并将新三角形ABPi加入到三角形集TS中，删除原选定的前沿Ei。

步骤9如果边集ES为空，转步骤10；如果边集ES不为空；转步骤4。

步骤10算法结束。

在本算法中采用半边进行三角网格剖分，每条半边只用一次，而且具有方向限定，简化了推进过程中的检测方法及判定数量。但是生成的三角网格是半边表示的，还要根据拓扑模型，即对应边关系，合并半边将其转化为普通三角网格。

3)利用保护技术简化复杂约束条件处理

对于复杂限定条件的区域剖分，限定条件的处理是最为困难的，尤其是复杂限定条件导致的小角度问题。本文采用两步法，首先利用Voronoi图建立限定条件(包括边界线)的保护区，小角度问题在保护区内解决；其次保护区以外的区域按照普通过程进行剖分；最后两部分区域合并作为最终的网格剖分结果。

由于Voronoi网格和Delaunay网格之间的对偶性，Voronoi网格的生成问题可以转换为Delaunay三角剖分的问题，而且通过Delaunay网格生成Voronoi网格的算法速度最快；同样Delaunay网格的生成问题也可以转换为Voronoi网格的生成问题来处理。 Voronoi图在网格的剖分、重剖分以及优化方面具有非常大的优势。

(1)Voronoi保护区

Branets et al.(2009)提出用保护区的概念将约束Voronoi网格生成问题转化为三角化问题。在限定条件PSLG P1周围建立保护区，保护区的边界构成新的限定条件 PSLGP2。P2的顶点顶点作为Voronoi单元的中心，并且利用Voronoi网格的边来精确表达内部特征。如此，三角剖分区域分为PSLG P2内部和外部，如图5所示。两部分分别构建的Delaunay三角网共享PSLG P2，因此整体区域实现Delaunay三角网剖分。

本文将Voronoi保护区的概念从限定条件推广到整个约束边界，在保护区内实现Voronoi图网格，保护区外采用前沿推进法布点，最终生成Delaunay三角网格。

对边界及限定条件进行细分以后，以折线上每个顶点为圆心作圆，称为保护圆，使得保护圆两两相交，并连接交点作为交线，然后根据保护圆和交线构造多边形，即Voronoi图，如图6所示。Voronoi图的构造方法如下：连接交线端点构造新边，如果边长大于半径，则取对应圆弧的中点，形成新边，直至所有边长均不大于半径为止；端点圆内只有一条交线，则取垂直平分交线的直线与圆的交点构造新边，重复截断过程。构造结束后，所有的新生成边构成折线的保护区。将折线端点与Voronoi图的顶点连接构造三角形，形成保护区内的Delaunay三角剖分。在后续剖分过程中，保护区内不允许插入点。从图7中示例可见，折线中相邻线段长度差距较大时，以线段为边的三角形质量稍差(可以通过平分法优化)，相邻线段长度近似时，三角形整体质量较好。

(2)复杂限定条件的Voronoi保护区

对于复杂的限定条件，如多条限定线交叉于一点的情况，处理思路类似。如图8中，四条限定线(标记为线1到线4，图中仅画出交叉部分)相交于一点O。以交叉点为圆心，一定长度(如交叉点网格尺寸)为半径作圆，称为交叉点的保护圆，作为交叉点处的保护区。交叉点附近区域的三角剖分必须在保护圆内实现。在保护区外布点时，保护圆内不能插入任何点。

在交叉点保护区内布点时，首先以一定距离(如交叉点网格尺寸或者最小角两点距离)在限定线两侧保护圆上布点。为了忠于相交的限定线，并且与限定线的保护区连续，限定线两侧的点对需要关于限定线对称，如点对A-B，B-C，D-E，F-G，则限定线垂直平分点对连线，分别交于点p1、p2、p3、p4。之后在保护圆远离限定线的位置布点，如H、I，使得两点之间的距离不大于半径。连接布点及交叉点，则构成保护区的三角剖分，见图8。

对于如线1和线2之间形成的尖角，只在角平分线上放一个点B，同时减小点对 AB、BC的距离，以保证垂直平分关系，但是这种剖分方法进一步减小夹角，降低剖分质量。同时其他限定线附近也可能出现小角度问题。为此采用边对分法优化。

以线1为例，如图8所示，线1与弦AB交于点p1，线2与弦BC交于点p2，删除边Op1两侧的三角形，取Op1中点p’1、Op2中点p’2、OI中点I’，连接A p’1、B p’1、B p’2、p’1 p’2、AI’、C p’2，将线1和线2附近的三角形以边半分法优化。同时继续处理受影响的邻接三角形，如图9所示。只考虑保护圆内的三角剖分，在优化过程中应取限定线与保护圆的交点构造新三角形；由于要与限定线保护区连续，所以取限定线与两侧点对连线的交点。

(3)保护圆的半径

Voronoi图保护区构造的主要困难在于如何确定保护圆的半径和间距，涉及到局部和全局的几何约束以及目标网格密度和标准。

首先考虑间距问题。因为保护圆以边界线段的端点为圆心，所以保护圆之间的间距(圆心的距离)实际上是边界线段的长度。

其次考虑半径问题。保护圆的生成有几个要求：①同一限定条件的相邻保护圆必须相交；②保护圆不能互相包含；③相邻保护圆的交点构成的线段最好是圆心线段的垂直平分线；④不同限定条件的保护圆必须远离，距离最少能容下一个最近的保护圆，本条主要针对两个约束离得太近却未相交的情况。

对于符合网格区域预期尺寸的均匀分段，设段长为r。如果半径小于r/2，则两圆相离；如果半径为r/2，则两圆相切，只有一个交点，无法构成边，如图10所示；如果半径大于r，则包入其他圆心点，交点连线过长；因此可以取3r/4。

对于边界分段不均匀的部分，由前文可知，α∈[0，1]表示相邻段的长度变化。设相邻分段长度分别为r0＝r/(2--α)、r1＝r、r2＝(2--α)r。由于2-α∈(1，2)，因此r0/2＜3r1/4＜r0， r1/2＜3r2/4＜r1，满足①②两条要求；并且相邻保护圆的交点构成的线段与圆心线段的交点在圆心线段的中点附近，偏差为(r2-r1)/2＝(1-α)r/2，满足第③条要求。

无论分段是否均匀，对第④条要求均需要检测，如果不满足，则需要适当调整保护圆的半径。实际上在进行分段细分和交叉点保护区控制之后，这种情况极少出现。

4、约束曲面网格优化技术

初始生成的三角网格，不可避免地包含一些质量较差的单元。为提高三角网格的质量，可以通过优化算法处理。总观网格优化算法，大致分为两类：改变网格微观拓扑结构的拓扑优化法和移动顶点位置的网格光顺法。拓扑优化法包括不改变顶点的对角线交换法和增删顶点的细分法，拓扑优化只改变顶点之间的连接；网格光顺法通过移动顶点以改变三角形的边长和角度，以改善三角形的形态。拓扑优化法限于初始三角网格的顶点分布，不能从根本上提升网格质量，从效果来看不及网格光顺法。本发明主要讨论网格光顺类的约束面三角网格优化技术。

1)Laplacian光顺法

①基于网格拓扑模型的算法改进

Field(1988)提出一种简单有效的优化方法，Laplacian光顺法。Laplacian光顺法的思路是考察每个顶点的邻接点，计算邻接点的平均坐标，将待考察点移动到平均坐标位置。本文基于前文建立的网格拓扑模型提出改进的Laplacian算法(算法4.4)。算法4.4改进的Laplacian算法

步骤1建立待优化顶点数组Q1、Q2，建立顶点优化数据E，优化计数器Count，以约束面所有顶点初始化数组Q1；

步骤2取Q1中一个顶点p(xi，yi，zi)，

步骤2.1如果点p为控制点，则点p存入数组Q2对应位置，转步骤2；否则转步骤2.2；

步骤2.2自拓扑模型TM中查找所有以p为起点的有向边，这些有向边的终点即为顶点p的邻接点。

步骤3计算所有邻接点的平均坐标点q(xj，yj，zj)，将q存入数组Q2对应位置，作为优化后的p点；

步骤4如果Q1为空，则转步骤5；否则转步骤2；

步骤5计算每个顶点的优化能量，即优化前与优化后的距离，e＝f(p，q)＝||p-q||，存入数组E对应位置；

步骤6交换数组Q1和Q2的元素；

步骤7取Q1中一个顶点p(xi，yi，zi)，

步骤7.1如果顶点的优化能量e＜ε(ε表示距离阈值)，则点p存入数组Q2对应位置，Count加1，转步骤7；否则转步骤7.2；

步骤7.2自拓扑模型TM中查找所有以p为起点的有向边，这些有向边的终点即为顶点p的邻接点。

步骤8计算所有邻接点的平均坐标点q(xj，yj，zj，)，将q存入数组Q2对应位置，作为优化后的p点；

步骤9如果Q1为空，则转步骤10；否则步骤7；

步骤10如果Count＝顶点数，转步骤11，；否则转步骤5；

步骤11输出结果网格，算法结束。

②针对曲面的改进

传统的Laplacian算法，包括改进的Laplacian算法1，没有考虑曲面的逼近问题。在曲面条件下，使用顶点p的邻接点计算的平均坐标点q可能不在原曲面上，此时q 与邻接点构成的三角面片与原曲面之间产生一定距离，且此距离的最大值是点p与点q 分别所处的平行平面之间的距离，称为“最小距离”。为保持原曲面的形态，当点p与点q的“最小距离”超过一定限度时，则点p不能直接移动到点q的位置，本文将其平行移动到与点q垂直的位置。

如图11所示，以线为例说明。点p的邻接点分别是A和C。以点A和点C计算平均坐标点O，点O与点p的距离为D、“最小距离”为d，点O在过点p且平行于AC 的直线上的投影为O’，则OO’长度为d。如果d＞ε，则点p只可以移动到点O’

如此，对顶点的优化能量计算进行改造，以解决逼近曲面的要求，公式如下：

e2＝d＝2*S/||AC||

其中，D＝||p-O||表示点p与点O之间的欧氏距离，即优化能量；d表示点p与点 O之间的“最小距离”；e1表示点p与点O’之间的距离，e2表示点O与点O’之间的距离；S为ABC构成的三角形的面积，已知三顶点坐标可用海伦公式求得。

针对曲面逼近的问题，与线的情形类似，如果点p与平均坐标点O的“最小距离”d＞ε，则点p移动到点O在过点p且“最小距离”为d的平面的投影位置。由于点p的邻接点不一定共面，因此取点O的法向垂直于过点p的切面，则过点O的平面方程为：

A(x-x₀)+B(y-y₀)+C(z-z₀)＝0

其中(x0，y0，z0)为点O坐标，(A，B，C)为法向，点p到平面的距离为：

其中，(x1，y1，z1)为点p坐标。

2)基于质心Voronoi图的优化方法

Voronoi图又称为泰森多边形，以其良好的特性应用于多种领域。Voronoi图是根据一组离散点进行空间划分的结果，通常称离散点为种子点，记为S(Aurenhammer， 1991)。每个种子点p∈S对应一个Voronoi单元，在此单元内的点距离距离这个种子最近。在三维欧氏空间，种子点p的Voronoi单元p定义为

空间中任意对象Ω与种子点S的Voronoi图的交集称为种子S对Ω的受限Voronoi图。

Voronoi单元Vp与Ω的交集称为p对Ω的受限Voronoi单元，定义为V_p∩Ω＝Vp∩Ω。由于Voronoi图是研究空间的一个分割，包含在空间中的任何对象Ω的受限Voronoi 图是对象的一个分割。

连接相邻离散点构成三角形，取三角形各边的垂直平分线相交构成的图形就是Voronoi图，交点即Voronoi图的顶点。一般情况下，种子点与Voronoi图的质心是不一致的。如果每个种子点与其Voronoi图的质心重合，则此Voronoi图称为质心Voronoi 图(Duet al.，1999)。如果每个种子点位于其受限Voronoi图的质心，则此受限Voronoi 图称为受限质心Voronoi图，Vp∩Ω(Du et al.，2003；Yan et al.，2009)。

受限Voronoi图可以通过计算优化种子的位置。Lloyd算法迭代移动每个种子到它的受限Voronoi单元的质心。Liu et al.(2009)提出更快的Newton-like算法，实际结果表明Newton-like法是适用的(Yan et al.，2009)。

利用Centroidal Voronoi Tessellation(CVT)将给定数量的点适当地布于约束表面和边界线附近。在恢复约束单元的同时，根据这些点的Voronoi图与约束单元的交集，运用拓扑控制以确定最终表面的顶点和三角形。可以采用Voronoi图法对所布点位置进行优化。方法如下：

首先，放置固定数量的种子，以保证种子是约束面的良好采样，如前文采用前沿推进法所布点。种子的数量决定了模型将被重新划分的分辨率，已知目标边长，可以用约束面的面积除以目标边长的平方根计算所需种子数量。种子对约束面的受限 Voronoi图分割约束面并将每个种子关联到约束面的一小部分。

然后，采用Lloyd法或者Newton-like法对种子进行迭代计算，将种子移动到其受限Voronoi图的质心，或者接近质心。

迭代计算结束后，将会获得质量较好的受限Voronoi单元，其对偶三角形质量也较好，同时每个单元的面积相近(Pellerin et al.，2014)。

将前文采用前沿推进法所布顶点作为种子，进行CVT优化，优化后的顶点分布更为合理，能够消除前沿推进法布点过程中在前沿冲突时产生的问题。注意保护区内的顶点不参与优化。

储层构造模型的边界重构是储层网格化的基础。储层构造模型边界元素的质量不但控制模型对地质对象的表达准确程度，而且严重影响储层网格的剖分质量。

地质对象的形态通过边界面表达，而边界处(尤其是交界处，包括内部特征面) 通常是地质对象形态、属性变化最剧烈的位置，相对而言内部通常比较平缓。同样，边界面上不同区域的变化程度也可能不同。储层网格化的两个重点，一是对地质对象表面形态的逼近，另一个是对变化程度的反映；其关键之一即在于对边界面的刻画。

本发明研究重构边界对储层表面的逼近程度，采用能量泛函(距离)对重构过程进行控制，提出简洁的控制函数；以角点、边界线、约束面的层级递进顺序计算各级别元素的网格尺寸；根据网格尺寸控制实现边界线和约束面的重构。采用前文设计的模型结构改进前沿推进法实现面内布点和三角剖分；采用Voronoi图保护区实现复杂限定条件的处理；并且研究表面三角网格优化方法，并对Laplacian算法进行改进，实现空间曲面的光顺。

以上实施例仅用以说明本发明的技术方案，而非对其进行限制；尽管参照前述实施例对本发明进行了详细的说明，对于本领域的普通技术人员来说，依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本发明所要求保护的技术方案的精神和范围。

Claims

1.一种约束曲面多分辨率控制预处理方法，其特征在于：所述方法包括以下步骤：

S4：边界面网格优化：使用改进的Laplacian光顺法进行平面网格优化，再使用改进的Laplacian光顺法进行约束面网格优化；

所述三角网格剖分的方法为拓扑模型的改进的前沿推进算法：

步骤4：取边集ES中的边Ei，设Ei的端点为A和B，则Ei表示为

以Ei为基础，即选定Ei作为选定前沿，搜索临近Ei且位于Ei左侧的顶点Pi加入LS，使用右手准则；

步骤6：取LQ中一点Pi，与边Ei组成三角形ABPi；

步骤7：判断新三角形是否满足以下要求：

①边

与已生成的边不能相交；

②新三角形面积的必须大于零；

如果不满足要求转步骤6；如果满足要求转步骤8；

步骤8：将新顶点与选定前沿端点连成的有向边

步骤10：算法结束；

所述S4中所述改进的Laplacian光顺法：

步骤2：取Q1中一个顶点p(xi,yi,zi)；

步骤3：计算所有邻接点的平均坐标点q(xj,yj,zj)，将q存入数组Q2对应位置，作为优化后的p点；

步骤4：如果Q1为空，则转步骤5；否则转步骤2；

步骤5：计算每个顶点的优化能量，即优化前与优化后的距离，e＝f(p,q)＝||p-q||，存入数组E对应位置；

步骤6：交换数组Q1和Q2的元素；

步骤7：取Q1中一个顶点p(xi，yi，zi)；

步骤7.1：如果顶点的优化能量e<ε，其中ε表示距离阈值，则点p存入数组Q2对应位置，Count加1，转步骤7；否则转步骤7.2；

步骤9：如果Q1为空，则转步骤10；否则步骤7；

步骤10：如果Count＝顶点数，转步骤11；否则转步骤5；

步骤11：输出结果网格，算法结束。

2.根据权利要求1所述的约束曲面多分辨率控制预处理方法，其特征在于：所述S1具体为：以角点、边界线、约束面的层级递进顺序计算角点网格尺寸、边界线网格尺寸、约束面网格尺寸；所述边界线网格尺寸的计算采用从角点开始的尺寸传播算法；所述约束面网格尺寸的计算采用从边界线开始的尺寸传播算法。

3.根据权利要求1所述的约束曲面多分辨率控制预处理方法，其特征在于：所述边界线重构还包括边界线逼近控制的处理方法：在预处理过程中，对于重构边界线中新生成的每一条线段，通过边界线距离控制函数检测原线段的端点到线段的最大距离，如果最大距离大于一定阈值，则调整线段端点重新构造线段并检测，直到新线段满足要求为止；

所述边界线距离控制函数定义为：

f(E_ab)＝Max(||P_iE_ab||)<ε，P_i∈L，P_i’∈E_ab

其中，L表示原边界线，P_i表示原边界线的线段端点，E_ab表示新生成的线段，P_i’表示P_i在线段E_ab所在直线上的投影，ε表示距离阈值。

4.根据权利要求1所述的约束曲面多分辨率控制预处理方法，其特征在于：所述S3中采用曲面前沿推进法在原约束面上布置顶点。

5.根据权利要求4所述的约束曲面多分辨率控制预处理方法，其特征在于：所述S3还包括：使用改进的前沿推进算法进行三角网格剖分结束后，根据拓扑模型，即对应边关系，合并半边将其转化为普通三角网格。

6.根据权利要求1所述的约束曲面多分辨率控制预处理方法，其特征在于：步骤S3中对于复杂限定条件的区域剖分，具体步骤包括：

(2)保护区以外的区域按照普通过程进行剖分；

(3)两部分区域合并作为最终的网格剖分结果。

7.根据权利要求1所述的约束曲面多分辨率控制预处理方法，其特征在于：所述步骤S3还包括约束面逼近控制的处理方法：对于重构曲面中新生成的每一个面片，根据约束面的距离控制函数检测原约束面的面片顶点到面片的最大距离，如果最大距离大于一定阈值，则调整面片顶点位置，构造新面片并检测，直到新面片满足要求为止；

所述约束面的距离控制函数定义为：

f(T_abc)＝Max(||P_iT_abc||)<ε，P_i∈S，P_i’∈T_abc

其中，S表示原约束面，P_i表示原约束面的顶点，T_abc表示新生成的面片，P_i’表示P_i在三角形T_abc平面内的投影，ε表示距离阈值。

8.根据权利要求1所述的约束曲面多分辨率控制预处理方法，其特征在于：所述S4还包括针对曲面的改进方法：在曲面条件下，使用顶点p的邻接点计算的平均坐标点q不在原曲面上，q与邻接点构成的三角面片与原曲面之间产生一定距离，且此距离的最大值是点p与点q分别所处的平行平面之间的距离，称为“最小距离”，当点p与点q的“最小距离”超过一定限度时，将其平行移动到与点q垂直的位置。