CN110322547B - 一种储层自适应四面体剖分方法 - Google Patents

Abstract

本发明提供了一种储层自适应四面体剖分方法，所述方法包括体域内自适应尺寸计算、体域内布点、四面体剖分、局部细分，本发明通过研究模型空间复杂性及其对网格尺寸的影响，详细分析了约束条件内部的距离计算方法及其与自适应网格尺寸的关系，提出复杂条件下约束体域期望网格尺寸的计算方法；在边界已剖分的基础上采用前沿推进法进行空间内四面体顶点布置；使用约束Delaunay法对已布顶点进行四面体剖分。针对网格质量优化及分辨率提高的需求，研究了储层网格细分方法及特点。

Description

一种储层自适应四面体剖分方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种储层自适应四面体剖分方法。

背景技术

储层网格化过程及其结果对于表征储层的非均质性十分重要。网格是地质属性的载体，是三维地质建模和油藏数值模拟之间连接的纽带和桥梁。储层网格的划分不仅应考虑已知数据的数学特征，同时应考虑储层的地质规律。因此，网格化过程必须与储层地质特征紧密联系在一起，使其保持所蕴含的地质意义。由于储层的非均质性，储层不同区域的属性不一致，体现在网格上则是网格大小不一致，即要求自适应大小的四面体剖分。

储层的四面体剖分是在储层构造模型及人工对象(如水力压裂缝等)的约束下实现的四面体剖分。储层构造模型由地层面、断层面、工区建模边界组成。储层自适应四面体网格剖分结果应达到以下要求：(1)准确保持约束单元形态，与约束面共形；(2)覆盖整个域空间；(3)满足尺寸自适应的要求和质量控制。由于储层的复杂性，包括形态复杂性和属性非均质性，出现了许多特殊情况，比如复杂断裂、裂缝等产生的复杂限定条件，比如由于地层不整合产生的小角度，比如薄层，比如关注区域导致的分辨率变化要求，比如网格总数的控制等。

目前储层四面体网格剖分的基本策略有两个：一是全局剖分策略，即所有单元一起剖分；一是逐一剖分策略，即每个约束单元独立剖分。其中逐一剖分策略明显优于全局剖分策略。在剖分过程中要保证相邻单元的连续性，即相邻单元共享边界离散/剖分的一致性。储层四面体剖分的常规解决方法和过程已经成熟，主要包括三个步骤：(1)根据需要重剖分约束面；(2)根据期望标准对区域内部插入点，作为剖分的顶点；(3)将区域内的所有点进行四面体剖分，约束边界上的顶点、线段、三角面等等均作为剖分的限定条件出现。常用四面体剖分方法主要有前沿推进法和Delaunay法。自适应网格剖分技术是指网格剖分算法在指定的目标条件(如尺寸)的前提下对研究区域进行网格剖分，剖分过程中网格尺寸受到目标条件的限制，剖分后网格的大小分布与预期一致。自适应四面体网格生成可以分为以下几种：(1)合并几个网格(Montenegro et al.，2009)；(2)细分一个粗网格(González-Yuste et al.，2004)；(3)约束于一个网格大小属性或密度函数(Romain et al.，2011)生成。网格大小属性可以使用到界面、曲线或中轴的符号距离计算。

发明内容

针对现有技术存在的问题，本发明提供了一种储层自适应四面体剖分方法。

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

一种储层自适应四面体剖分方法，它包括以下步骤：

(1)体域内自适应尺寸计算：首先假设一个期望网格尺寸，建立背景网格，生成区域内的网格尺寸控制点；根据期望网格尺寸，从水平和纵向检测区域的边界距离，根据边界距离判定网格尺寸，取两个网格尺寸的最小值作为区域的期望网格尺寸；

(2)体域内布点：使用前沿推进法在约束体域内进行正四面体网格布点；

(3)四面体剖分：使用约束Delaunay四面体剖分方法进行四面体剖分；

(4)局部细分：使用点插入法和边对分法进行局部细分。

进一步的，所述步骤(1)中所述期望网格尺寸的确定还包括：

表面距离的计算：分别计算水平方向、竖直方向和正交方向的距离，取最小值作为边界距离；

体内距离的计算：取垂向距离和法向距离的最小值作为三角平面片对应的边界距离。

进一步的，所述步骤(1)中采用叉树法建立规则四边形/六面体背景网格，具体为：在储层网格剖分中，给出期望网格尺寸，即剖分后的网格尺寸不能大于期望网格尺寸，以期望网格尺寸为最小网格划分空间区域，二维区域建立网格数为n*n的网格矩阵，三维区域建立n*n*n的网格矩阵，在此基础上对每个网格持续细分，直到满足要求为止。

进一步的，所述步骤(1)中对于二维平面区域采用四叉树细分方法，对于三维空间区域采用八叉树细分方法。

进一步的，所述步骤(1)还包括在四叉树背景网格中，取每个叶子节点的正方形网格，在中心处插入一点，其属性值为网格边长，代表该正方形网格范围内的网格尺寸；在八叉树背景网格中，取每个叶子节点的立方体网格，在中心处插入一点，其属性值为网格边长，代表该六面体网格范围内的网格尺寸，网格中心的插入点即为区域内的网格尺寸控制点。

进一步的，所述步骤(2)中所述前沿推进法布点具体为：

步骤1：建立表示推进前沿的三角半面片队列Q，建立顶点集合P；

步骤2：约束面的所有三角半面片加入队列Q；

步骤3：从Q中取一个三角半面片t，表示为顶点(v1，v2，v3)，令v表示三角半面片t的质心，mv表示约束体域中距v最近的网格尺寸，表示t的法向量，则待求顶点为/>和/>其中/>|.|是绝对值函数，||.||是L2-范数，由前沿方向，舍弃p2；

步骤4：判断顶点p1是否满足条件：一是在约束体域内，二是不太靠近已布点，如果p1满足条件，将p1加入集合P，将由p1和以p1为中心、以1.5×mp1为半径的球体内的所有点对构成的所有三角形半面片插入队列Q，mp1是距p1处最近的网格尺寸；

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

否则转步骤6；

步骤6：算法结束。

进一步的，所述步骤(2)中判断p1是否位于约束体域内部的方法为：

将三维空间区域进行约束区域四面体剖分，该剖分不考虑质量问题，只要剖分成四面体即可；

判断点是否被至少一个四面体所包含，若是，则在约束体域内，否则，位于约束体域外。

进一步的，所述步骤(2)中新顶点构造所述前沿三角半面片的方法：新前沿三角半面片的三条边为原推进前沿有向边的对边和新顶点与原边端点连接构成的两条新半边，三条边方向一致，由原前沿三角半面共产生9条半边、3个三角半面片。

进一步的，所述步骤(4)中局部细分具体包括：以四面体中底面三角形形态和第四点到底面距离的大小进行分类，分为较近、适中、较远三类，根据分类选择细分方法，远近的具体标准根据实际模型确定。

进一步的，所述步骤(4)中局部细分时，常规四面体细分使用8单元细分法；不规则四面体细分，如果网格中的小角度不可改变，采用保角对分法；如果只是已有顶点不可变，采用点插入法细分；当细分四面体一边划分3段及以上，不利于邻接四面体的剖分判断，采用保角对分法细分。

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

(1)本发明提出一种新的储层自适应四面体剖分方法，通过研究模型空间复杂性及其对网格尺寸的影响，详细分析了约束条件内部的距离计算方法及其与自适应网格尺寸的关系，提出复杂条件下约束体域期望网格尺寸的计算方法；

(2)在边界已剖分的基础上采用前沿推进法进行空间内四面体顶点布置；

(3)使用约束Delaunay法对已布顶点进行四面体剖分；并采用拓扑结构设计算法；

(4)针对网格质量优化及分辨率提高的需求，研究了储层网格细分方法及特点。在网格中存在尖角而又不能消除时，提出“保角对分法”以保护小角度不再恶化，在网格质量较好时适合采用边对分法提高分辨比率，网格质量差时适合采用Delaunay点插入法。

附图说明

图1为自适应四面体剖分流程图；

图2为距离影响网格分辨率的情况示意图；

图3为曲率影响分辨率的示意图，Caumon et al.，2009；

图4为剖分面区域的四叉树背景网格示意图；

图5为空间任意点的网格尺寸计算示意图；

图6为二维表面内部距离计算示意图；

图7为体域内前沿推进法布点示意图；

图8为约束Delaunay四面体剖分示意图，其中：a、四面体abcp是约束Delaunay的；b、三角形abc是局部Delaunay的；

图9为二维Delaunay点插入法示意图；

图10为边对分法示意图；

图11为四面体边对分基本类型，Escobar，2005；

图12为四面体边对分法分解模板，liu&Joe，1996；Frank，2006；

图13为三棱柱的四面体剖分示意图；

图14为splinte四面体细化剖分示意图；

图15为四-三棱柱四面体剖分示意图；

图16为四棱柱四面体剖分示意图。

具体实施方式

下面结合附图和实施例对本发明的技术方案做进一步详细的说明。

实施例1、一种储层自适应四面体剖分方法

本发明提供了一种新的自适应四面体剖分方法，图1为自适应四面体剖分流程图。

1、模型(体域)空间复杂性及网格尺寸计算

1.1模型(体域)空间复杂性分析

在使用三维储层模型进行空间分析和数值模拟时，模型网格具有两个特点：(1)必须达到一定精度才有意义或者符合要求，即网格尺寸不是完全随意和自动的，应该有一个大概预期；(2)根据空间形态和特征，针对空间区域自动确定网格单元尺寸，控制整体网格数量，满足应用需要而又使网格数量尽量少。因此，在对模型进行四面体剖分之前，首先需要确定四面体网格的期望尺寸。

针对一个具有若干约束体域的三维储层模型，不同约束体域的网格剖分精度控制分为两种情况：一是人为指定每个约束体域的网格精度(即网格单元的期望尺寸)；另一种是根据约束体域自动确定内部的网格单元尺寸。而且人为指定的方式下，可以将多个约束体域指定相同的网格单元尺寸，也可以分别指定不同的网格单元尺寸。人为指定网格单元尺寸的方式不多解释，本发明讨论网格单元尺寸的自动计算方法。网格剖分时一般会给出全局网格单元尺寸，根据约束区域的特征计算网格单元尺寸时受到全局网格单元尺寸的约束。

衡量几何要素复杂性的四个度量方面：数量、距离、角度和曲率，其中数量是针对模型中约束体域个数的度量，由于本发明采用约束体域逐一网格化的策略，可以忽略数量这一因素。因此本发明重点从距离、角度、曲率三个方面讨论模型(体域)空间复杂度。

(1)距离

从空间整体上看，决定空间分辨率大小的第一决定因素是距离，即区域各个方向的长度范围。储层空间中的距离度量主要分为水平和纵向两个维度，一般情况下，水平范围远远大于纵向范围，因此网格精度通常由纵向距离决定。但是当水平方向因特殊原因出现小距离时，仍需要特殊处理。距离影响网格分辨率的情况如图2所示。

水平距离由界定平面区域(上顶/下底边界面)的边界线决定，如图2上排所示，大概可以划分为以下几类：a.常规型，指平面区域各方向距离相差不大，各方向对网格分辨率影响一致，或者各方向距离足够大，在各方向上对网格分辨率没有影响；b.狭长型，指平面区域一个方向距离较大，另一方向距离相对较小，网格分辨率主要由较小的距离控制；c&d.收紧型，指平面区域某个方向上整体距离较大，而在另一个方向上，部分距离较小，另一些距离较大；e.不规则型，指平面区域没有明显的方向性，而且区域形态不规则。需要说明的是，水平距离的方向与坐标系没有严格的对应关系，在划分类型时所指的方向也并不一定是坐标轴方向。比如对于河流相储层，一个方向可能是河道方向，另一个方向则是垂直河道的方向。因为河道通常是弯曲的，所以河道方向是曲线变化的，其垂直方向也在不断变化。

纵向距离指界定约束体域的顶底边界面的距离。与水平距离类似，因地质要素形态不同，纵向距离也分为不同类型，如图2下排所示：a.规则型，指储层在纵向上各处距离比较均匀，且距离较大对网格分辨率没有影响；b.扁平型，指储层在纵向上各处距离比较均匀，但距离较小需要提高网格分辨率；c&d.收缩型，指储层在纵向上各处距离不均匀，距离较大的区域满足网格分辨率，而距离较小的区域需要提高网格分辨率；e.尖角型，这种类型主要是因为断层的原因导致在纵向上造成夹角而形成，此种类型的夹角点一定是角点。

约束体域在计算网格精度时，需要同时考虑水平距离和纵向距离变化，从而获得自适应的网格尺寸。

(2)角度与曲率

曲率衡量的是几何体(曲线、曲面)的不平坦程度，如果几何体上某点P的曲率为0时，几何体在该点处是平坦的(直线或者平面)，曲率半径无穷大；如果点P处曲率大于0(曲率不能负)，则表示点P处是不平坦的，曲率越大不平坦程度越高。曲率变化对分辨率的影响如图3所示。曲率的计算比较复杂，可以用平坦线/面对几何体的逼近程度来代替表征曲率。以边界线为例，边界线的弯曲在顶点处。在顶点P两侧距离为r处分别取点A和B，取AB中点O，则线段AB与点P的距离(PO)表示了线段AB对边界线在顶点P处的逼近程度，此距离与AP和PB的夹角直接相关。

1.2、约束体域网格尺寸计算方法

为了控制网格的分辨率，通常需要控制网格的尺寸。而目标网格尺寸需要在重构之前根据约束条件和整体分辨率来构建。目标网格尺寸可以使用一个统一属性值，这样比较简单；也可以采取背景网格的方式，建立一个代表网格尺寸分布状态的背景网格。使用背景网格法能够较快地计算插入点的网格尺寸，应用比较广泛。建立背景网格可以采用四边形/六面体，也可以采用三角形/四面体，本发明采用叉树法建立规则四边形/六面体背景网格。

叉树法的原理比较简单，分为自顶向下和自底向上两种思路。以二维四叉树为例，自顶向下的基本思想是将平面区域等分为四个子区域，判断等分后的子区域内属性是否一致，如果一致则此平面区域不再划分，否则将子区域继续等分为四个子区域，以此类推，直至子区域内属性一致或者子区域足够小为止。二维四叉树自底向上的基本思路是，先按照最小尺寸将平面区域进行n*n划分，然后判断相邻四个平级网格的属性，如果属性一致则合并为一个网格，否则此四个网格作为叶子节点。三维空间使用八叉树，每次将空间区域等分为八个子空间。叉树法建立的四叉树或八叉树通常作为空间索引结构使用，需要注意树的深度和平衡性。本发明采用叉树法进行空间划分不需考虑此问题。

(1)叉树法建立尺寸背景网格

背景网格的建立本身是一个细分的过程，本发明采用自顶向下的叉树建立思路。在储层网格剖分中，人为给出期望网格尺寸，即剖分后的网格尺寸不能大于期望网格尺寸。以期望网格尺寸为最小网格划分空间区域，二维区域建立网格数为n*n的网格矩阵，三维区域建立n*n*n的网格矩阵。在此基础上对每个网格持续细分，直到满足要求为止。空间区域内部的网格一般比较均匀，网格细分主要出现在边界及其邻近区域。

四叉树细分，见图4。对于二维平面区域采用四叉树细分方法建立背景网格。假设期望网格尺寸为1，将平面区域划分为n*n个边长为1的正方形网格。针对每一个正方形网格进行检测，如果网格内不包含边界线段，则网格作为叶子节点，不再细分；如果网格内包含边界线段且满足一定条件(如，将所包含的线段的长度的平均值取对数，结果小于网格边长的一半)，则将网格划分为四个相等的正方形子网格。再对生成的正方形子网格进行细分，以此类推，直到所有的网格都不再细分为止。

八叉树细分。对于三维空间区域，采用八叉树细分方法建立背景网格。假设期望网格尺寸为1，将空间区域划分为n*n*n个边长为l的立方体网格。针对每一个立方体网格进行检测，如果网格内不包含边界三角面片，则网格作为叶子节点，不再细分；如果网格内包含边界三角面片且满足一定条件(如，将所包含的三角面片的面积的平均值开立方，结果小于网格边长的一半)，则将网格划分为八个相等的立方体子网格。再对生成的立方体子网格进行细分，以此类推，直到所有的网格都不再细分为止。

(2)生成网格尺寸控制点

在四叉树背景网格中，取每个叶子节点的正方形网格，在中心处插入一点，其属性值为网格边长，代表该正方形网格范围内的网格尺寸。

在八叉树背景网格中，取每个叶子节点的立方体网格，在中心处插入一点，其属性值为网格边长，代表该六面体网格范围内的网格尺寸。

网格中心的插入点即为区域内的网格尺寸控制点。在网格剖分过程中，每一点处的网格尺寸受到网格尺寸控制点的约束。约束方式有两种，一种是取距离最近的约束点的尺寸，另一种是取周围所有约束点的(加权)平均。图5为空间任意点的网格尺寸计算，图中扁圆点为网格尺寸控制点，点p为尺寸待定点。点p处的网格尺寸或者取距离最近的c1点的尺寸，或者取能够影响点p的临近点c₁、c₂、c₃、c₄的(加权)平均尺寸。

1.3复杂性因素对自适应网格尺寸的影响

空间区域网格剖分中所谓的自适应分辨率，具有两层含义。一种是指在同一约束体域中，由于区域几何要素的限制，导致区域内网格尺寸不能全部达到期望网格尺寸而造成的分辨率不一致，呈现逐渐过渡的样式；另一种是指在相邻的不同区域中，由于不同区域几何要素的各自限制，导致不同区域中的期望网格尺寸不一致的情况，而且相邻区域的分辨率在边界处要达到一致。

在具体讨论之前，先假设一个整体期望网格尺寸，设为L(边长)。那么计算每个约束体域的期望网格尺寸就是检验L对该约束体域的适应性问题。

在前文中，本发明将模型(体域)的空间复杂性归为距离和角度两个因素考量，复杂度越高的区域越是需要高分辨率网格剖分。此处，角度因素同样以距离方式处理。

(1)以距离确定网格尺寸

首先从水平和纵向检测区域的边界距离，即区域的局部直径D1、D2，D1为水平方向距离、D2为垂直方向距离。所谓局部直径，在二维区域指边界线构成线段之间的距离，在三维区域指边界面构成平面片之间的距离。以L为标准尺度对边界距离进行衡量，依据空间距离的分类，可能会出现以下情形：①水平距离和纵向距离均远大于标准尺度L，则该区域内的期望网格尺寸即为L；②水平距离远大于L而纵向距离小于L，则着重考察纵向距离来判定期望网格尺寸，反之亦然；③水平距离或者纵向距离部分远大于L，而部分小于L，则统计距离大于和小于L的区域所占的比重，来判定期望网格尺寸。

对于边界距离整体小于L的区域，显然期望网格尺寸L不适用，需要减小尺寸，网格尺寸计算分水平和纵向分别进行，以水平距离为例，方法如下：取边界距离的平均值与L进行比较，如果平均值的平方根小于L/2，则取L的一半作为新L，再次将边界距离平均值的平方根与L比较，以此类推，直到不满足条件为止，此时L即为该区域水平距离决定的网格尺寸。同理求得纵向距离决定的网格尺寸后，取两个网格尺寸的最小值作为区域的期望网格尺寸。

对于边界距离部分大于L、部分小于L的情况，除了判断边界距离与L的差异程度外，还需要判断距离小于L区域所占的比重。如果边界距离小于L的区域比重较大，则区域的期望网格尺寸以小距离区域决定的网格尺寸为准；如果比重较小，则期望网格尺寸采用L。

(2)表面距离计算方法

由于区域边界形态可能不规则，因此计算区域内某一部分的边界距离时，需要确定对应的边界元素，对应元素之间的距离作为这一部分的边界距离，即局部直径。由于边界形态不规则，计算局部距离的对应元素不能事先确定。通过分析发现，如图6所示二维边界示例，边界线上某一线段的对应元素仅可能在三个方向：水平方向、竖直方向和正交方向。分别计算三个方向的距离，取最小值作为区域的局部边界距离(局部直径)。

线段之间的距离是两条线段上最近点之间的距离，三角面片之间的距离是两个三角面片上最近点之间的距离。为简化计算，可以取线段对应端点之间的距离的平均值作为线段之间的距离，取三角面片之间对应顶点的距离的平均值作为三角面片之间的距离。

由于某些线段端点在边界线上的对应点不一定是端点，因此距离计算方法修正如下：

水平方向距离计算。过线段两个端点A和B分别向AB正方向(右手定则)作水平射线，交边界于另两点C和D，线段AC和BD的长度均值为d1，即对于线段AB，其水平方向局部直径为d1。

竖直方向距离计算。过线段两个端点A和B分别向AB正方向作竖直射线，交边界于另两点E和F，线段AE和BF的长度均值为d2，即对于线段AB来说，其竖直方向局部直径为d2。

正交方向距离计算。过线段两个端点A和B分别向AB正方向垂直于线段AB作射线，交边界于另两点G和H，线段AG和BH的长度均值为d3，即对于线段AB来说，其正交方向局部直径为d3。

如果交点CD之间存在顶点P，通过顶点P作AC的平行线交AB于点O，则d1为线段AC、BD、OP的长度平均值。竖直方向和正交方向边界距离同样处理。

(3)体内距离计算方法

三维约束体域中的边界距离采用三角面片之间的距离计算，与二维中的线段类似，以三角面片的三个顶点为基础计算平均距离。边界面上的距离可以通过二维方法计算出来，因此此处只需要计算体内距离。沉积储层一般是水平的，即体内距离是其垂向距离。考虑褶皱、反转等变形情况，同时计算三角面片的法向距离(右手定则)。取垂向距离和法向距离的最小值作为该三角平面片对应的局部直径(局部距离)。方法如下：

垂向距离计算。过三角形三个顶点A、B、C分别作纵向作射线，交边界面于三点E、F、G，则线段AE、BF、CG的平均长度即为三角形ABC对应的垂向距离。

法向距离计算。过三角形三个顶点A、B、C分别垂直三角形ABC平面作射线，交边界面于三点E、F、G，则线段AE、BF、CG的平均长度即为三角形ABC对应的垂向距离。

计算垂向距离时，如果在EFG区域内有原边界面的顶点P，则从P平行AE向三角形ABC作射线，交三角形ABC于点O，则垂向距离为线段AE、BF、CG、OP的平均长度。对法向距离同样方法处理。

2、体域内布点

储层三维结构模型四面体剖分后生成的四面体网格，应该是与构造模型共形的，即构造模型的表面三角网格正是四面体网格的外表面。用于储层模型网格重构的方法不同，最终得到的四面体网格的质量也是不一样的。

在四面体剖分算法中，重构的构造模型边界面采用带方向的半边结构表示，并且根据体域形态建造了表示网格尺寸的背景立方体网格。采用逐一剖分的策略，利用体域的边界三角网格作为约束，将构造模型中的各个约束体域一个一个地进行四面体剖分。

四面体剖分过程分为两个步骤：第一步，在约束体域内布点，生成最终的四面体顶点；第二步，在约束体域边界面三角网格的约束下对生成的顶点集合进行四面体剖分，生成四面体网格。

本发明以前沿推进法在约束体域内进行正四面体网格布点，需要注意的是，此处“正四面体”表示四面体的四条边尽量一致而不是完全相等，下文中出现时意义相同，不再说明。从约束体域的约束面开始，针对约束面上的每一个三角面片，向约束体域内部计算潜在的、能够构成正四面体的顶点。在每一个计算过程中，如果潜在点满足了一定条件，则添加该点，否则舍弃。

体域内四面体剖分的布点算法采用队列数据结构(算法1)，见图7。

算法1体域内前沿推进法布点

步骤1建立表示推进前沿的三角半面片队列Q，建立顶点集合P。

步骤2约束面的所有三角半面片加入队列Q；

步骤3从Q中取一个三角半面片t，表示为顶点(v1，v2，v3)，令v表示三角半面片t的质心，mv表示约束体域中距v最近的网格尺寸，表示t的法向量。则待求顶点为/>和其中/>|.|是绝对值函数，||.||是L2-范数，见图7。由前沿方向，舍弃p2。

步骤4判断顶点p1是否满足条件：一是在约束体域内，二是不太靠近已布点。如果p1满足条件，将p1加入集合P，将由p1和以p1为中心、以1.5×mp1为半径的球体内的所有点对构成的所有三角形半面片插入队列Q，mp1是距p1处最近的网格尺寸。

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

否则转步骤6；

步骤6算法结束。

点在平面多边形内的判定方法已经非常成熟，比如射线判交法。而点在三维空间内的判定比较麻烦，原因在于空间形态的复杂性导致可能出现的特殊情况非常多。有人提出将三维区域投影到平面进行判断的方法，但是投影平面不容易选择，目前没有成熟的直接判定方法。

为了测试p1是否位于约束体域内部，本发明提出一种简化三维空间，进而进行判断的方法。方法分为两步：第一步将三维空间区域进行约束区域四面体剖分，该剖分不考虑质量问题，只要剖分成四面体即可；第二步判断点是否被至少一个四面体所包含，若是则在约束体域内，否则位于约束体域外。

为了测试p1是否距离已布顶点集合P中另一个顶点太近，建立一个中心为p1、半径为γmp1的球体，其中γ∈[0，1]由用户给出，作为p1的排除区域。如果排除区域内没有任何其他点，则表示顶点p1距离其他顶点不太近。

新顶点构造前沿三角半面片的方法如下：新前沿三角半面片的三条边为原推进前沿有向边的对边和新顶点与原边端点连接构成的两条新半边，三条边方向一致(同为顺时针或者逆时针)。由原前沿三角半面共产生9条半边、3个三角半面片。

3、约束Delaunay四面体剖分

假设空间区域内存在一个分段线性复合形(PLC)，记为P。将P中顶点进行Delaunay四面体剖分，其结果网格可能不会覆盖P中的边和面。Lee&Lin(1986)提出平面内的约束Delaunay四面体剖分(CDT，constrained Delaunay tetrahedralization)方法，Shewchuk(1998；2008)把CDT推广到三维及更高维。约束Delaunay四面体剖分是Delaunay四面体剖分的一种变形，必须涵盖P中的边和面。

设面片f∈X，点p和q位于这个平面片f的两侧，线段pq与平面片f相交，见图8，由于平面片f的分隔，点p和q相互不可见。如果一个四面体t的顶点在X内，其外接球不包含X的其他点，则该四面体是符合约束Delaunay(constrained Delaunay)的，这些点相对四面体t内部的任何点是可见的，见图8a。

假设T是P的一个四面体网格剖分结果，t、t₁、t₂是T中的四面体单元，s是T中的一个三角形。如果s只属于t，或者是t₁和t₂的公共面，而且t₁和t₂的其他顶点不在s的外接球中，则称s是局部Delaunay(locally Delaunay)的。换言之，t₁的外接球不包含t₂的其他顶点，t₂的外接球不包含t₁的其他顶点，见图8b。如果T中每个三角形都是局部Delaunay的，并且P中的任意限定面都不包含T内的三角形，则T是约束Delaunay四面体剖分(CDT)的。

约束Delaunay四面体剖分(CDT)的基本准则与Delaunay四面体剖分是一样的，不过略有不同。对CDT而言，边界多边形内的三角形可以不是局部Delaunay的。总体而言，CDT保持了Delaunay四面体剖分的大部分有利特性。值得注意的是，CDT产生的四面体网格内的单元，包括四面体、三角形和边，不一定符合Delaunay准则，这是对于特殊现象的妥协。

对于一个任意的分段线性复合形P，不一定能够进行约束Delaunay四面体剖分。为确保CDT的存在，必须加入Steiner点，加入的Steiner点集记为S；加入Steiner点集后，约束条件从P变为P∪S，记为P’；最终的约束四面体Delaunay剖分(CDT)是P’的CDT，记为SCDT；与限定Delaunay四面体剖分相比，SCDT过程中添加的Steiner点通常更少；特别是当P中存在一些小角度现象时，限定Delaunay四面体剖分非常困难，而且需要大量的Steiner点；CDT可以很容易地解决这种问题(Si，2015)。

4、储层四面体网格局部细分方法

四面体网格细化的目标是提高四面体网格分辨率，提高四面体网格质量。以一定规则在四面体网格中插入新点，新点与原有节点组成新四面体构成空间的细分。通常有两种方式，一种是将点插入四面体内部，称为点插入法；一种是插入四面体边上，称为分边法；如果插入点位于边的中心，则称为边对分法。插入面上的方法非常少见。Frank(2006)对于四面体网格局部细分方法有较为详细的总结。

4.1、常用细分方法

1)点插入法

点插入算法基于Delaunay三角剖分。四面体网格的Delaunay准则是，在一个四面体的外接球内没有其他四面体的顶点。Chew(1989)提出基于Delaunay准则的二维网格细化方法，Shewchuk(1998；2008)将其扩展到三维。点插入法采用Bowyer-Watson算法实现，使用以下步骤来执行Delaunay细化三角剖分：

(1)向一个四面体的重心或者外接球的中心插入一个新点；

(2)从网格中删除因新点而违反Delaunay准则的四面体，产生一个空腔；

(3)将空腔的表面面片与插入点连接建立新的四面体。

对每一个插入点，临近的四面体受到影响。点插入方法的一个明显的缺点是，元素的大小只略有下降，算法收敛缓慢。二维Delaunay三角网格细化如图9所示，将外接圆替换为外接球，点插入法可以直接扩展到三维，对四面体网格进行细化。

不过，点插入法破坏了原四面体网格的邻接关系。如果在空间中存在限定条件，则不能保证限定条件的保持性。

2)边对分法

二维空间中的边对分通过分割三角形的边把三角形分成两个三角形。Plaza&Carey(1996)提出了在三维空间该算法的扩展及其在四面体的应用。每当四面体网格的一条边被分割，网格中共享这条边的所有四面体必须被分裂成两个四面体，见图10。对于局部网格细化，Arnold et al.(2000)提出了一种递归细化方案以取得一致性。这项工作已有的实现是使用最长边分割以实现共形网格。边对分的一个基本优点是容易实现，但其缺点是网格质量缺乏保证。

边对分法中，四面体边上的加点方式可以分为以下几种基本类型：

类型1如果四面体六条边上各有一个插入点。此时，四面体被分为八个小四面体，见图11a，划分方法由Bank et al.(1983)提出。

类型2如果四面体中有4或5条边上有插入点。则在另1、2条边上增加插入点，四面体被细分非八个小四面体。之前，没有新节点的边必须被标记。

类型3如果四面体中有3条边上有插入点，

(1)如果3条边不位于同一个面，那么在其他3条边上增加插入点，四面体被细分为八个小四面体；

(2)如果3条边位于四面体同一个面，那么创建4个小四面体，见图11b。

类型4对于四面体中有2条边上有新插入点的情况，

(1)如果2条边不位于同一个面，则构造4个小四面体，见图11c。

(2)如果2条边位于同一个面，则创建3个小四面体，见图11d。其中最长的边固定作为参考，以利用一些情形下最长边平分的性质。

类型5如果四面体中只有1条边上有插入点，细分为2个小四面体，见图11e。

类型6如果四面体没有新节点，则不用划分，见图11f。

一旦四面体细分被实现，需要搜索周边四面体，在它们的边上可能有0到6个新节点被引入，必须处理以保证网格的一致性。

新四面体网格的一致性采用以下策略来保证：如果任何小四面体需要产生并且它的质量度量小于规定的细分阈值qθ，那么将这个小四面体的父四面体加入到类型1的四面体组中。这样，新产生的网格中的元素的最小质量可以提高。

3)基于模板的8单元细分法

Plaza&Rivara(2003)提出一种连续的边对分连同预定义分区模板以分解四面体为八个单形。模板定义了一定的分解模式，分别分裂三角形或四面体的。Liu&Joe(1996)提出分解四面体为八个类似的单形，见图12a。因此每个面被分成四个相似的三角形。这种细化模式在第一步产生四个四面体和一个八面体。为了分解八面体成四个四面体，需要选择一个分裂边。这个决定应该基于质量度量。如果细化具有局部特性，则需要额外的模板以获得共形网格，见图12b、图12c、图12d。

已有实现使用如图12所示的分解模板。所有指定分解的单纯形被图12a的分解模板处理。分解区域的边界也必须细化，以实现共形网格。必须应用于单纯形的模板由分割的边数定义。图12的模板仅定义如下分裂模式：一条边被分割(图12b)，边上的两边分裂不共享相同的面(图12c)，边上的三边分裂共享相同的面(图12d)，六边分裂的主分解(图12a)。任何其他配置必须从给定的分解模板推导。边共面的双边分裂单纯形，按照图12d的模板分解成四个单形。四或五边分裂的单形，根据图12a的模板分解为八个四面体。

4)边对分法分析

(1)8单元细分(图11a、图12a)。8单元细分法将四面体分为8个四面体，其中原四面体的4个顶点处，各产生一个原四面体相似的小四面体，其角度与原四面体对应一致，边长为原四面体对应边长的一半；中间部分也分为四个小四面体，其角度与原四面体接近，四元素共享边长度与连接方式有关，大于连接边的一般，其他边为对应边的一半。

(2)四单元细分一(图11b、图12d)。这种细分方式，将原四面体四面体的三个面的角一分为二，底面边长减半，新生成的四个四面体边长比加大、角度变小、变狭长。

(3)四单元细分二(图11c、图12c)。这种细分方式，将原四面体每个面、每条边进行了对分，新生成的四面体角度减少，变扁平。

(4)三单元细分(图11d)。这种细分方式，将原四面体一分为二之后，又将其中一半一分为二，新生成的四面体角度减小、变狭长。

(5)二单元细分(图11e、图12b)。这种细分方式，将原四面体一分为二，新生成的四面体角度减小、变狭长。

可见，在边对分四面体细分方法中，只有8单元细分，也就是所有边都对分的情况下，新生成的小四面体会保持与原四面体角度和形态接近，且分辨率提高迅速；其他情况均导致角度变小、四面体边长比大幅变化的情形。因此，对于正四面体或者边长比较小的四面体，其细分适合采用8单元细分法；而对于狭长、扁平等四面体则需要特别考虑。

4.2、“保角对分”法

细分的目的是提高分辨率和网格质量，尤其是对于特殊限定的四面体网格。提高分辨率以缩小网格边长方式实现，提高网格质量则有几个标准：(1)网格尽量均匀，趋向正四面体；(2)小角度网格比例尽量低；(3)尽量不增加小角度网格。因此，对于新加点位置，需要处理。

边对分法中，除了8单元细分以外，其他细分在对分边的同时还细分了边的对角，造成越是细分角度越小的问题。假如该角的度数本身就小，则进一步加剧了小角度的趋势。大量小角度的存在必然造成网格质量的下降。因此对于特殊四面体，不能简单的采用边对分法，需要特别处理，期望在减小边长的同时，不要过多地降低角度。

本发明以四面体中底面三角形形态和第四点到底面距离的大小进行分类，分为较近、适中、较远三类，分别讨论处理方法。远近的具体标准根据实际模型确定。极近与极远的情形做退化处理，此处不讨论。

1)底面为近正三角形

底面为正三角形或者接近正三角形时，根据第四点划分为三类远近位置进行，分为近、中、远三种情况，以第四点为依据。三种情况分别处理如下：

(1)如果距离较近，则第四点构成三个面的夹角较大，取底面三边中点，以四单元细分二方式细分；

(2)如果距离适中，此时四面体边长相近，角度适中，则采用8单元细分方式细分；

(3)如果距离较远，则第四顶点构成的角度较小，如果采用边对分，则继续减小角度；如果采用8单元细分，不能减少奇异四面体。因此可以采取特殊办法。根据以第四点为端点的三边中最短边与底面中最长边的比值n＝L/l，将三边均分为n份，则原四面体分解为1个与原四面体相似的小四面体和n-1份三棱柱。将每个三棱柱分为三个四面体，方式如下：取三棱柱侧面中角度最大的顶点，与对边作三角形，将三棱柱分为一个四面体和一个四棱锥；取四棱锥底面较短对角线与对点作三角形，将四棱锥分为两个四面体，见图13。

2)底面为钝角三角形

底面为钝角三角形时，无论第四点距离远近，先取钝角对边中点，以边对分法中的2单元细分方式细分。

3)底面为狭长锐角三角形

当底面为狭长锐角三角形时，必有一边较短而另两边相对较长，否则将会构成钝角三角形，其较长边至少是短边长度的2倍以上。同样根据第四点与底面的距离分为近、中、远三种情况。

(1)第四点与底面距离较近或者适中

如果第四点与短边相近，将第四点与短边构成的三角形看做底面，则四面体转化为底面为近正三角形、而第四点距离较远的情况。

如果第四点与底面锐角顶点距离相近，处理如下：取第四点与底面两长边中点构成三角形将原四面体分解与一个四面体与一个五面体；取底面较大角的顶点与对面两边中点构成三角形将五面体分解为一个四面体与一个五面体；取原四面体四条长边的中点构造四边形将五面体分为两个四棱锥；取四边形较短对角线与四面形相对顶点构成三角形，将四棱锥分为两个四面体，见图14。

(2)第四点与底面距离比较远

这里作一个规定，以第四点为端点的边与底面长边长度相仿，再远的情况不作考虑。在这种情况下，以底面短边为基准对其他长边进行划分。将原四面体分为n个小四面体(n为层数)，n-1个四-三棱柱(底面为四边形，顶面为三角形)，n-2个四棱柱。

对于四-三棱柱ABCDA’B’C’，见图15，取底面四边形较短对角线AC，与顶点B’构造三角形；再取侧面对角线A’D、中轴对角线B’D，分别与顶点A’、C、C’构成三角形，将将四-三棱柱分为五个小四面体。

对于四棱柱，见图16，其剖分结果比四-三棱柱多一个小四面体，见图16。

4)保角对分

考虑到细分四面体与邻接四面体的关系，如果一边划分3段及以上，不利于邻接四面体的剖分判断，因此可以考虑保角对分，即在分割边上仍然取中点，但是不做对角分割，以避免小角的进一步减小。如果底面为狭长锐角三角形，第四点距离底面较远，则取较长的5边中点，将原四面体分为2个小四面体，1个四-三棱柱(底面为四边形，顶面为三角形)。用前文方法将四-三棱柱将分为四个小四面体。

在保角对分方法下，细化后的四面体内不会产生比原四面体最小角更小的角，且不会产生比原四面体边长比更大的边长比。

4.3、细分过程及算法

储层四面体网格的细分随网格本身的状态和细分目的不同差异较大。如果四面体网格区域内没有特殊限定，则完全可以使用Delaunay点插入法，虽然其细分速度较慢，却可以校正四面体形态，提高网格质量；如果四面体网格区域内有限定条件，为保证限定条件不被破坏，使用边对分法比较好。在边对分细分法中，8单元细分法对四面体所有边进行对分，不需要判断四面体的形态，最为简单。

1)常规四面体细分

常规四面体细分认为区域内四面体形态比较好，不存在或者很少存在扁平、狭长等特殊四面体。为提高其分辨率，对四面体网格进行整体细分，可以全部使用8单元细分法，其处理相对简单(算法2)。

算法2常规四面体网格整体细分

步骤1初始化四面体队列Q，期望边长L；

步骤2从Q中取一四面体T，判断T的边长LT是否满足条件LT＜L，如不满足，对T进行8单元细分，生成新单元T1...T8，存入Q；

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

否则，转步骤4；

步骤4算法结束。

如果是对部分区域进行四面体网格细分，为保持网格一致性，需要对区域周围受到影响的四面体网格进行细分，直到没有受影响的四面体为止(算法3)。

算法3常规四面体网格部分细分

步骤1初始化四面体队列Q1，初始化期望边长L，初始化邻接四面体队列Q2

步骤2从Q1中取一四面体T，判断T的边长LT是否满足条件LT＜L，如不满足，对T进行8单元细分，生成新单元T1...T8，存入Q1；

步骤3如果Q1不为空，则转步骤2；

否则，转步骤4；

步骤4从Q2中取一四面体T，判断T的对分边，并根据对分边的位置分别细分T1...Tn，n＜＝8；

步骤5判断T的对分边关联的邻接四面体，记入队列Q2；

步骤6如果Q2为空，则转步骤4；

否则，转步骤7；

步骤7算法结束。

2)不规则四面体细分

不规则四面体指比较狭长或者扁平的四面体。在特殊区域，比如薄层或者小角度，如果开始按照大边长、低分辨率方法剖分，则容易造成大量不规则的四面体。在这些区域，一方面要控制不规则四面体的继续增加，一方面是对已产生的不规则四面体进行质量改善性细分，即在细化过程中保持四面体的质量不再下降。假设网格中的小角度不可改变，采用保角对分法，即在保证小角的前提下，尽量细分边，使得细化后四面体的边长与原四面体中较短边接近；如果只是已有顶点不可变，可以采用点插入法细分，虽然迭代次数较多，但是效果明显优于边对分法。

储层网格化的目的是将空间进行合适划分，以便于表现空间内的属性分布及变化。理想的划分应该是根据空间形态和属性变化自适应地调整网格大小，并且保持网格本身形态良好、网格之间变化连续且和缓。

针对储层特点，本发明提出一种新的储层自适应四面体剖分方法。通过研究模型空间复杂性及其对网格尺寸的影响，详细分析了约束条件内部的距离计算方法及其与自适应网格尺寸的关系，提出复杂条件下约束体域期望网格尺寸的计算方法；在边界己剖分的基础上采用前沿推进法进行空间内四面体顶点布置；使用约束Delaunay法对已布顶点进行四面体剖分；并采用拓扑结构设计算法。

此外，针对网格质量优化及分辨率提高的需求，研究了储层网格细分方法及特点，Delaunay点插入法细分质量好、边对分法细分速度快。在网格中存在尖角而又不能消除时(如尖角位于控制点处)，提出“保角对分法”以保护小角度不再恶化。在网格质量较好时适合采用边对分法提高分辨比率，网格质量差时适合采用Delaunay点插入法。

Claims (6)

1.一种储层自适应四面体剖分方法，其特征在于：它包括以下步骤：

(1)体域内自适应尺寸计算：首先假设一个期望网格尺寸，建立背景网格，生成区域内的网格尺寸控制点；根据期望网格尺寸，从水平和纵向检测区域的边界距离，根据边界距离判定网格尺寸，取两个网格尺寸的最小值作为区域的期望网格尺寸；

所述期望网格尺寸的确定还包括：

表面距离的计算：分别计算水平方向、竖直方向和正交方向的距离，取最小值作为边界距离；

体内距离的计算：取垂向距离和法向距离的最小值作为三角平面片对应的边界距离；

(2)体域内布点：使用前沿推进法在约束体域内进行正四面体网格布点；

所述前沿推进法布点具体为：

步骤a：建立表示推进前沿的三角半面片队列Q，建立顶点集合P；

步骤b：约束面的所有三角半面片加入队列Q；

步骤c：从Q中取一个三角半面片t，表示为顶点(v1,v2,v3)，令v表示三角半面片t的质心，mv表示约束体域中距v最近的网格尺寸，表示t的法向量，则待求顶点为/>和其中/>|.|是绝对值函数，||.||是L2-范数，由前沿方向，舍弃p2；

步骤d：判断顶点p1是否满足条件：一是在约束体域内，二是不太靠近已布点，如果p1满足条件，将p1加入集合P，将由p1和以p1为中心、以1.5×mp1为半径的球体内的所有点对构成的所有三角形半面片插入队列Q，mp1是距p1处最近的网格尺寸；

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

否则转步骤6；

步骤f：算法结束；

所述判断p1是否位于约束体域内部的方法为：

将三维空间区域进行约束区域四面体剖分，该剖分不考虑质量问题，只要剖分成四面体即可；

判断点是否被至少一个四面体所包含，若是，则在约束体域内，否则，位于约束体域外；

(3)四面体剖分：使用约束Delaunay四面体剖分方法进行四面体剖分；

(4)局部细分：使用点插入法和边对分法进行局部细分；

所述局部细分具体包括：以四面体中底面三角形形态和第四点到底面距离的大小进行分类，分为较近、适中、较远三类，根据分类选择细分方法，远近的具体标准根据实际模型确定。

2.根据权利要求1所述的储层自适应四面体剖分方法，其特征在于：所述步骤(1)中采用叉树法建立规则四边形/六面体背景网格，具体为：在储层网格剖分中，给出期望网格尺寸，剖分后的网格尺寸不能大于期望网格尺寸，以期望网格尺寸为最小网格划分空间区域，二维区域建立网格数为n*n的网格矩阵，三维区域建立n*n*n的网格矩阵，在此基础上对每个网格持续细分，直到满足要求为止。

3.根据权利要求2所述的储层自适应四面体剖分方法，其特征在于：所述步骤(1)中对于二维平面区域采用四叉树细分方法，对于三维空间区域采用八叉树细分方法。

4.根据权利要求3所述的储层自适应四面体剖分方法，其特征在于：所述步骤(1)还包括在四叉树背景网格中，取每个叶子节点的正方形网格，在中心处插入一点，其属性值为网格边长，代表该正方形网格范围内的网格尺寸；在八叉树背景网格中，取每个叶子节点的立方体网格，在中心处插入一点，其属性值为网格边长，代表该六面体网格范围内的网格尺寸，网格中心的插入点即为区域内的网格尺寸控制点。

5.根据权利要求1所述的储层自适应四面体剖分方法，其特征在于：所述步骤(2)中新顶点构造所述前沿三角半面片的方法：新前沿三角半面片的三条边为原推进前沿有向边的对边和新顶点与原边端点连接构成的两条新半边，三条边方向一致，由原前沿三角半面共产生9条半边、3个三角半面片。

6.根据权利要求1所述的储层自适应四面体剖分方法，其特征在于：所述步骤(4)中局部细分时，常规四面体细分使用8单元细分法；不规则四面体细分，如果网格中的小角度不可改变，采用保角对分法；如果只是已有顶点不可变，采用点插入法细分；当细分四面体一边划分3段及以上，不利于邻接四面体的剖分判断，采用保角对分法细分。