CN106875487B - 一种基于邻域作用力的地质六面体网格平滑方法 - Google Patents

一种基于邻域作用力的地质六面体网格平滑方法 Download PDF

Info

Publication number
CN106875487B
CN106875487B CN201710116290.XA CN201710116290A CN106875487B CN 106875487 B CN106875487 B CN 106875487B CN 201710116290 A CN201710116290 A CN 201710116290A CN 106875487 B CN106875487 B CN 106875487B
Authority
CN
China
Prior art keywords
vertex
block
neighborhood
hexahedron
smoothing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201710116290.XA
Other languages
English (en)
Other versions
CN106875487A (zh
Inventor
王海起
刘玉
彭佳琦
陈冉
桂丽
翟文龙
费涛
闫滨
车磊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN201710116290.XA priority Critical patent/CN106875487B/zh
Publication of CN106875487A publication Critical patent/CN106875487A/zh
Application granted granted Critical
Publication of CN106875487B publication Critical patent/CN106875487B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/10—Constructive solid geometry [CSG] using solid primitives, e.g. cylinders, cubes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Image Generation (AREA)

Abstract

本发明公开了一种基于邻域作用力的地质六面体网格平滑方法，包括以下步骤：遍历三维六面体模型每个顶点；判断顶点是否需进行平滑处理；对需平滑的顶点，将共有该顶点的邻域六面体网格根据属性值划分为不同的体块，属性值相同的网格属于同一体块；将每个体块看作磁铁的“同极”，各同极对顶点具有相斥的作用力，作用力采用基于重力模型的空间相互作用力模型进行度量，计算顶点在各体块单独作用力下的移动方向和距离；进而计算顶点在所有体块合力作用下沿x轴、y轴、z轴三个方向的移动分量；从而获取顶点最终的平滑位置。本发明在平滑过程中既保持了网格原有的三维拓扑关系，又降低了三维地质区域边界表达的锯齿效应，使平滑结果符合实际的地质现象。

Description

一种基于邻域作用力的地质六面体网格平滑方法
技术领域
本发明涉及三维地质体网格平滑方法领域，尤其是一种基于邻域作用力的地质六面体网格平滑方法。
背景技术
真三维地质体建模的模型体元通常采用六面体、四面体、棱柱体等，其中，六面体网格具有组织形式规则、网格单元数量及重划分次数较少、计算效率高等优点，在三维地质体、油藏体建模中应用广泛，并被大多数商业地质建模、油藏建模和数值模拟软件所支持。然而，这种源自矩形网格的体元在对三维区域边界进行离散化网格表达时存在锯齿效应。高精度的精细地质模型可以降低锯齿效应的影响、提高表征精度，但即使对于一个中等规模的地质体或油藏体，精细的地质表征往往具有百万级甚至千万级的网格规模，受限于实际实现技术，精细地质模型和油藏数值模拟器可以支持的网格规模之间存在着较大差距，对三维网格进行粗化就是将精细地质模型粗化到油藏数值模拟器能够接受的网格规模。然而粗化后的六面体体元较大，对圈闭、油层等区域边界的表达由于锯齿效应的存在不符合实际的地质现象，影响模型的描述精度、可视化效果及后期模拟计算的准确性。
对于六面体模型边界锯齿现象，可以采用边界网格重构或加密方法解决，这种方式可以构造出精细的边缘形状，如图1所示，右图对边界面上凹凸不平的六面体网格进行加密和重构，达到了精细的平滑效果，但对于地质六面体模型组织结构来说，这种方法打乱了原有结构的规则行列网格检索方式，为查找和数值模拟计算增加了难度。另一种解决途径通过适当移动六面体的顶点进行平滑，这种方式只能进行粗略的平滑，难以刻画精细的边缘形状，但保留了原有的网格组织形式，网格数和顶点数均没有发生改变。
目前对于网格顶点平滑，经典的平滑算法包括Laplacian平滑、Taubin平滑、平均曲率法等。Laplacian平滑的核心是将网格内部节点的位置移动到与该节点共面节点组成的多面体的体心处，算法简单。Taubin算法在Laplacian算法的基础上引入了滤波器及权系数，可抑制拉普拉斯算子引起的变形收缩。平均曲率法则遵循曲面曲率变化均匀即为光滑的原则。上述方法各有优缺点，不少学者亦在此基础上进行了大量的优化改进，但是对于三维地质网格模型，这些平滑方法在应用时会产生一些问题。
三维地质网格模型不同于其他建筑或工具的三维模型。图1所示为某工程零件的三维六面体网格模型，这类模型所有网格属性相同，只需对表面因形状产生的锯齿进行平滑。而地质网格模型，模型表面并不需要平滑，模型内部因为网格属性值不同发生的锯齿现象需要平滑，若分别将相同属性值的网格提取出来对表面进行平滑，平滑后的拓扑关系可能被破坏。
综上所述，如何能借助三维六面体模型网格顶点坐标点阵的查找优势，避免破坏原有的网格结构，又能保证网格之间的三维拓扑关系，是三维地质体网格平滑的关键问题。
发明内容
本发明技术解决问题：为解决现有技术存在的不足，提供一种基于邻域作用力的地质六面体网格平滑方法，在平滑过程中既保持了网格原有的三维拓扑关系，又降低了三维地质区域边界表达的锯齿效应，使平滑结果符合实际的地质现象。
为实现上述目的，本发明采用下述技术方案：一种基于邻域作用力的地质六面体网格平滑方法，遍历三维六面体模型每个顶点，判断顶点是否需进行平滑处理，对需平滑的顶点，将共有该顶点的邻域六面体网格根据属性值划分为不同的体块，属性值相同的网格属于同一体块，将每个体块看作磁铁的“同极”，各体块对顶点具有相斥的作用力，作用力采用基于重力模型的空间相互作用力模型进行度量，计算顶点在各体块单独作用力下的移动方向和距离，进而计算顶点在所有体块合力作用下沿x轴、y轴、z轴三个方向的移动分量，从而获取顶点最终的平滑位置。
具体包括以下步骤：
步骤一：根据三维六面体模型网格顶点坐标集合依次遍历三维模型每个顶点，及共用该顶点的邻域六面体网格；对于每个顶点，进行下面步骤二至步骤七处理；
步骤二：对当前顶点的不同情况进行判断以决定是否进行平滑处理，对需平滑的顶点进行下面步骤三至步骤七处理，并忽略不需要平滑处理的顶点；
步骤三：计算当前顶点各邻域六面体网格质心坐标；
步骤四：将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值可分为多个体块，并根据各体块包含的邻域六面体网格质心坐标计算各个体块的质心坐标；
步骤五：根据当前顶点和各体块的坐标位置，以及各体块包含的邻域六面体网格个数，采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离；
步骤六：在获取各体块单独作用力下移动方向和距离的基础上，计算当前顶点在所有体块合力作用下沿x轴、y轴、z轴三个方向的移动分量；
步骤七：根据x轴、y轴、z轴三个方向的移动分量计算当前顶点最终的平滑位置。
进一步的，所述步骤二中，对当前顶点的不同情况进行判断以决定是否进行平滑处理，包括以下五种情况：
1、若顶点为三维模型边界面上的点，统计共用该顶点的邻域六面体网格在边界面上的属性值分布情况，若邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；若属性值不相同，在二维边界面上采用上述的步骤三至步骤七进行顶点平滑处理；
2、若顶点为三维模型的角点或边点，则不进行平滑处理；
3、若顶点为被人工标注过的断层、尖灭等不需要进行平滑的特征点，则不进行平滑处理；
4、若顶点在三维模型内部，且邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；
5、若顶点在三维模型内部，且邻域六面体网格属性值不相同，则采用上述的步骤三至步骤七进行顶点平滑处理。
进一步的，所述步骤三中，计算当前顶点各邻域六面体网格质心坐标方法为：
对于每个邻域六面体网格，设该六面体网格8个顶点的坐标为(xi,yi,zi)，其中i＝1,2,…,8，i为顶点序号，网格质心坐标为(xc,yc,zc)，质心坐标计算公式为：
Figure BDA0001235690040000031
进一步的，所述步骤四中，将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值可分为多个体块，根据各体块包含的邻域六面体网格质心坐标计算各体块质心坐标方法为：
设当前顶点的邻域六面体网格可分为多个体块{Bk}，k表示体块序号，每个体块Bk的质心坐标为(Xk,Yk,Zk)，体块Bk包含的邻域六面体网格个数为nk，体块Bk内各邻域六面体网格的质心坐标为
Figure BDA0001235690040000041
其中i表示体块Bk包含的邻域六面体网格序号，i＝1,2,…,nk，则体块Bk质心坐标计算公式为：
Figure BDA0001235690040000042
进一步的，所述步骤五中，根据当前顶点和各体块的坐标位置，以及各体块包含的邻域六面体网格个数，采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离，具体步骤为：
1、将当前顶点的每个体块看作磁铁的“同极”，各同极对顶点具有相斥的作用力，作用力采用基于重力模型的空间相互作用力SIM模型(Spatial Interaction Model)进行度量。
设当前顶点为v，顶点v坐标位置为(xv,yv,zv)，体块Bk质心坐标为(Xk,Yk,Zk)，k表示体块序号，顶点v质量为Mv，体块Bk质量为Mk，体块Bk对顶点v的相斥作用力向量
Figure BDA0001235690040000043
为：
Figure BDA0001235690040000044
式中，G、α分别为引力常数和质量指数，这里取1；
Figure BDA0001235690040000045
为体块Bk中心到顶点v连线向量，体块中心取质心；β为距离指数，亦取1。
Figure BDA0001235690040000046
方向由
Figure BDA0001235690040000047
确定，方向向量
Figure BDA0001235690040000048
计算公式为：
Figure BDA0001235690040000049
2、对于每个
Figure BDA00012356900400000410
取G＝α＝β＝1，由运动定律F＝ma可知顶点v在体块Bk作用力下的加速度向量
Figure BDA00012356900400000411
为：
Figure BDA00012356900400000412
Figure BDA00012356900400000413
方向同样由
Figure BDA00012356900400000414
确定。
3、在不考虑初始速度和时间的前提下，加速度
Figure BDA00012356900400000415
代表了顶点v在作用力
Figure BDA00012356900400000416
作用下的移动偏移量，即移动距离
Figure BDA00012356900400000417
将每个六面体网格的质量视为1，则体块Bk的质量Mk等于体块Bk包含的网格个数nk。对于三维六面体模型，任一顶点v的邻域六面体网格个数最多为8个，当所有邻域六面体网格均属于同一体块时，则意味着邻域六面体网格属性值均相同、属性分布是均质的，不需要进行平滑处理，因此体块Bk的邻域六面体网格个数nk≤7，即体块Bk的质量Mk最大为7。设顶点移动方向的邻域六面体网格对角线长度diag为顶点移动距离的上限，则顶点v在最大质量和移动上限约束下的移动距离
Figure BDA0001235690040000051
可表达为：
Figure BDA0001235690040000052
式中，γ为平滑系数，用于对移动距离
Figure BDA0001235690040000053
进行约束，约束的目的是使顶点v移动后的新位置能达到最佳平滑效果；nk为体块Bk包含的邻域六面体网格个数；
Figure BDA0001235690040000054
为体块Bk质心到顶点v的向量。
平滑系数γ取值范围为区间[0,1]，具体取值根据平滑效果评价指标确定。
平滑效果评价指标为当前顶点与其不同体块分界面上相邻顶点连线形成的夹角余弦值的平均值。具体地说，在二维视角下，评价指标为当前顶点与其处于不同体块分界线上的相邻顶点连线形成的夹角余弦平均值；在三维视角下，评价指标为当前顶点与XOY、XOZ、YOZ三个剖面上处于不同体块分界线上的相邻顶点连线形成的夹角余弦平均值。余弦值或余弦平均值越接近于-1(即夹角越接近于180度)，平滑效果越好。
γ取值方法：平滑系数γ在区间[0,1]依据固定步长分别取值(例如γ分别取值0.01、0.02、0.03、…)，将不同取值分别代入上式计算移动距离
Figure BDA0001235690040000055
获取相应的顶点v移动后的新位置，并根据移动后的位置分别计算相应的平滑效果评价指标，从而选择评价指标最接近于-1的γ值为最终的平滑系数值。
4、当前顶点v在体块Bk作用下的移动距离向量
Figure BDA0001235690040000056
最终表示为：
Figure BDA0001235690040000057
式中，(xv,yv,zv)为顶点v坐标；(Xk,Yk,Zk)为体块Bk质心坐标；nk为体块Bk包含的邻域六面体网格个数；diag为顶点v移动方向的邻域六面体网格对角线长度；γ为上一步骤确定的平滑系数值；
Figure BDA0001235690040000058
为体块Bk质心到顶点v的向量，
Figure BDA0001235690040000059
为向量
Figure BDA00012356900400000510
的模，
Figure BDA00012356900400000511
表示体块Bk与顶点v的距离。
进一步的，所述步骤六中，在获取各体块单独作用力下移动方向和距离的基础上，当前顶点v在所有体块{Bk}合力作用下沿x轴、y轴、z轴三个方向的移动分量为：
Figure BDA0001235690040000061
式中，Δdx、Δdy、Δdz分别为x轴、y轴、z轴方向的移动分量；γ为平滑系数；diag为顶点v移动方向的邻域六面体网格对角线长度；(xv,yv,zv)为顶点v坐标；(Xk,Yk,Zk)为体块Bk质心坐标；nk为体块Bk包含的邻域六面体网格个数；
Figure BDA0001235690040000062
为体块Bk质心到顶点v的向量
Figure BDA0001235690040000063
模。
进一步的，所述步骤七中，根据x轴、y轴、z轴三个方向的移动分量，当前顶点v(xv,yv,zv)最终平滑移动的新位置(x′v,y′v,z′v)为：
Figure BDA0001235690040000064
式中，Δdx、Δdy、Δdz分别为x轴、y轴、z轴方向的移动分量；(xv,yv,zv)为顶点v的当前坐标位置；(x′v,y′v,z′v)为顶点v的新坐标位置。
本发明的有益效果：本发明对需平滑顶点的邻域六面体网格根据属性值划分为不同的体块，根据体块对顶点的作用合力方向确定顶点平滑移动的方向，通过移动方向的邻域六面体网格对角线长度限制顶点移动的最大距离，使顶点平滑距离不超出邻域六面体网格的范围，从而在平滑过程中既保持了网格原有的拓扑关系，又降低了三维区域边界表达的锯齿效应，使其符合实际的地质现象。
附图说明
图1是六面体网格边界加密重构图，左图是原始图形，右图是对表面加密重构后的图形，椭圆虚线框用于标记对比的部分；
图2是本发明流程图；
图3是三维模型表面顶点不同位置分布示意图；
图4是邻域作用力示意图；
图5是数量比例为7:1的两个邻域体块作用力平滑示意图；
图6是数量比例为5:3的两个邻域体块作用力平滑示意图；
图7是数量比例为4:4的两个邻域体块作用力平滑示意图；
图8是数量比例为5:2:1的三个邻域体块作用力平滑示意图；
图9是5×5网格中间凸起网格的平滑效果图；
图10是剖面情况下两种属性的平滑效果图；
图11是邻域六面体网格的其他体块分布情况；
图12是不同夹角对应的余弦值；
图13是平滑系数γ＝0的平滑效果图；
图14是平滑系数γ＝0.275的平滑效果图；
图15是平滑系数γ＝0.6的平滑效果图；
图16是平滑系数γ和平滑效果评价指标的关系曲线；
图17是连续台阶状锯齿现象；
图18是直角转弯锯齿现象。
具体实施方式
下面结合附图对本发明进行详细说明：
在表面平滑的Laplacian算法中，网格顶点vi平滑移动到与其共面邻接顶点vj的平均位置，其移动偏移量
Figure BDA0001235690040000071
其中n是vi邻接顶点的个数，wj用于度量邻接顶点vj在平滑移动中的权重大小，最简单的形式是采用相同权重，即wj＝1/n。
在本发明方法中，对上述平滑思路进行扩展，以顶点周围共用该顶点的八个邻域六面体网格组成平滑窗口，将窗口内邻域六面体网格划分为不同的体块，不同体块对顶点相斥作用力的大小度量了顶点平滑移动的不同权重，所有体块作用力的合力决定了顶点最终移动的距离和方向。
利用发明方法对三维区域边界的锯齿进行平滑处理，目的是降低边界表达的噪声，使其符合实际的地质现象。对于有些顶点，例如断层、尖灭等，是不需要进行平滑处理的。当顶点需进行平滑处理时，其邻域的各体块总是将顶点外推，目的是使体块尖锐、锯齿状的边界面和边界点趋向光滑，这与地层分界面的连续性是一致的，只不过质量较大的体块将顶点外推的作用力较大，质量较小的体块将顶点外推的作用力较小，不同方向、大小的作用力形成合力最终将顶点推到相对平衡的位置，既降低了边界的锯齿效应，又保持了网格的稳定性。
发明方法的主要思路可概括为：将共用同一顶点的八个邻域六面体网格体划分为不同的体块，每个体块包含具有相同属性值的网格。将每个体块看作磁铁的“同极”，同极之间具有相斥的作用力，顶点移动的方向为体块作用力的方向，移动的距离以作用力的大小作为权重，作用力大则移动距离大，作用力小则移动距离小。所有体块合力的方向和大小决定了顶点最终移动的方向和距离。基于邻域作用力的地质六面体网格平滑方法就是根据地质属性分布的一般规律，在保证六面体网格之间三维拓扑关系不变的同时对网格顶点进行移动平滑。
如图2所示，一种基于邻域作用力的地质六面体网格平滑方法，包括以下步骤：
步骤一：根据三维六面体模型网格顶点坐标集合依次遍历三维模型每个顶点，及共用该顶点的邻域六面体网格；对于每个顶点，进行下面步骤二至步骤七处理。
步骤二：对当前顶点的不同情况进行判断以决定是否进行平滑处理，对需平滑的顶点进行下面步骤三至步骤七处理。
对当前顶点的不同情况进行判断以决定是否进行平滑处理，包括以下五种情况(如图3所示)：
1、若顶点为三维模型边界面上的点，统计共用该顶点的邻域六面体网格在边界面上的属性值分布情况，若邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；若属性值不相同，在二维边界面上采用下面步骤三至步骤七进行顶点平滑处理；
2、若顶点为三维模型的角点或边点，则不进行平滑处理；
3、若顶点为被人工标注过的断层、尖灭等不需要进行平滑的特征点，则不进行平滑处理；
断层、地层尖灭等区域覆盖的不需要进行平滑的特征顶点，应在平滑处理前，由地质工程师进行交互操作，圈定或输入一些特征点或特征区域，避免平滑时被破坏掉。
4、若顶点在三维模型内部，且邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；
5、若顶点在三维模型内部，且邻域六面体网格属性值不相同，则采用下面步骤三至步骤七进行顶点平滑处理。
步骤三：计算当前顶点各邻域六面体网格质心坐标。
对于每个邻域六面体网格，设该六面体网格8个顶点的坐标为(xi,yi,zi)，其中i＝1,2,…,8，i为顶点序号，网格质心坐标为(xc,yc,zc)，质心坐标计算公式为：
Figure BDA0001235690040000081
步骤四：将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值可分为多个体块，并计算各体块的质心坐标。
以图4示例，设顶点v周围八个邻域六面体网格的离散属性值代表不同的沉积微相类型，如用数值1、2、3分别代表坝主体、坝缘、前三角洲类型，相同的属性值指示相同的沉积微相类型，据此可将八个邻域六面体网格划分为三个体块B1、B2、B3，分别代表上述三种微相类型。图4为了便于展示，将顶点v与邻域六面体网格拉开了距离。
设当前顶点的邻域六面体网格可分为多个体块{Bk}，每个体块Bk的质心坐标为(Xk,Yk,Zk)，体块Bk包含的邻域六面体网格个数为nk，体块Bk内各邻域六面体网格的质心坐标为
Figure BDA0001235690040000091
其中i表示体块Bk包含的邻域六面体网格序号，i＝1,2,…,nk，则体块Bk质心坐标计算公式为：
Figure BDA0001235690040000092
步骤五：采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离。
1、将当前顶点的每个体块看作磁铁的“同极”，各同极对顶点具有相斥的作用力，作用力采用基于重力模型的空间相互作用力SIM模型(Spatial Interaction Model)进行度量。
设当前顶点为v，顶点v坐标位置为(xv,yv,zv)，体块Bk质心坐标为(Xk,Yk,Zk)，k表示体块序号，顶点v质量为Mv，体块Bk质量为Mk，体块Bk对顶点v的相斥作用力向量
Figure BDA0001235690040000093
为：
Figure BDA0001235690040000094
式中，G、α分别为引力常数和质量指数，这里取1；
Figure BDA0001235690040000095
为体块Bk中心到顶点v连线向量，体块中心取质心；β为距离指数，亦取1。
Figure BDA0001235690040000096
方向由
Figure BDA0001235690040000097
确定。图4中三个体块对顶点v的作用力
Figure BDA0001235690040000098
大小分别由上式(3)确定，方向分别由向量
Figure BDA0001235690040000099
确定，则合力
Figure BDA00012356900400000910
实际上是三个作用力向量之和。
方向向量
Figure BDA00012356900400000911
计算公式为：
Figure BDA00012356900400000912
2、对于每个
Figure BDA0001235690040000101
取G＝α＝β＝1，由运动定律F＝ma可知顶点v在体块Bk作用力下的加速度向量
Figure BDA0001235690040000102
为：
Figure BDA0001235690040000103
Figure BDA0001235690040000104
方向同样由
Figure BDA0001235690040000105
确定。在不考虑初始速度和时间的前提下，加速度
Figure BDA0001235690040000106
代表了顶点v在作用力
Figure BDA0001235690040000107
作用下的移动偏移量(即移动距离)
Figure BDA0001235690040000108
则在合力
Figure BDA0001235690040000109
作用下总的移动偏移量为各向量
Figure BDA00012356900400001010
之和。
公式(5)说明
Figure BDA00012356900400001011
大小与体块Bk的质量Mk成正比，与距离
Figure BDA00012356900400001012
成反比。若将每个六面体网格的质量简单地视为1，则
Figure BDA00012356900400001013
大小实际上与体块Bk包含的六面体网格个数nk成正比，表明当体块包含的网格个数较多时，对顶点的作用力较大，顶点移动距离较长，如图4中的体块B3；反之，则对顶点的作用力较小，顶点移动距离较短，如图4中的体块B2。
下面对网格体块分布的几种具体情况进行图示分析：
(1)顶点邻域八个六面体网格存在两种不同属性值，如图5所示，浅色网格个数为1，深色网格个数为7，左图为原始结构，中图深浅两色箭头为两种力的作用方向，最深色箭头(即黑色箭头)为合力的作用方向，右图为顶点沿合力方向平移后的效果。
(2)顶点邻域八个六面体网格存在3个浅色网格和5个深色网格，作用力效果和平滑效果如图6所示。
(3)顶点邻域八个六面体网格存在4个浅色网格和4个深色网格，如图7所示，上下两侧的作用力相互抵消，不需要平滑。
(4)顶点邻域八个六面体网格存在3种不同属性值，作用力效果和平滑效果如图8所示。
(5)如图9所示的5×5网格中心存在一凸起网格，左图箭头示意各顶点的移动方向，右图为顶点移动平滑后的效果。
(6)在剖面情况下，图10左图为处理前的网格，箭头示意相应顶点移动的方向，右图为顶点移动后的网格，可以看出，该方法同样也适用于二维网格的平滑。
对于共用顶点的八个六面体网格，若体块包含的网格数较少，例如只有1个六面体网格(如图5所示的浅色网格)，则意味着该网格是孤立的，是一个凸起，需要平滑凸起边界；反过来看，深色六面体网格的数量较多，可以看作成包围趋势，则意味着深色网格存在一个凹陷，需要平滑凹陷边界。在本发明提出的算法中，借助“同性相斥”的作用力思想，给网格数量较多的体块一个较大的力，给网格数量较少的体块一个较小的力，从而合力的方向和大小决定了顶点移动的方向和距离，这种方法的结果也符合前述的平滑规律。
对于具有更多属性值的网格、或者属性值交错杂乱分布的情况，如图11所示，每个六面体网格几乎都相当于一个异变(凸起或凹陷)，这种情况一般不需要进行平滑处理，应保持其原有的特征。对应到本发明方法，若属性分布杂乱，则会产生四面八方的力，这些力互相抵消、或完全抵消，则顶点仅需做细微的移动、或不移动，因此在复杂情况下，发明方法也符合一般的平滑规律。实际上，在地质体中，地质属性的分布一般是连续、有一定规律的，属性在一个顶点周围杂乱分布的情况较为少见。
3、将每个六面体网格的质量视为1，则体块Bk的质量Mk等于体块Bk包含的网格个数nk。对于三维六面体模型，任一顶点v的邻域六面体网格个数最多为8个，当所有邻域六面体网格均属于同一体块时，则意味着邻域六面体网格属性值均相同、属性分布是均质的，不需要进行平滑处理，即前述步骤二中第4种情况，因此，体块Bk的邻域六面体网格个数nk≤7，即体块Bk的质量Mk最大为7。设顶点移动方向的邻域六面体网格对角线长度diag为顶点移动距离的上限，则顶点v在体块Bk作用下基于最大质量和移动上限约束的移动距离
Figure BDA0001235690040000111
可表达为：
Figure BDA0001235690040000112
式中，γ为平滑系数，用于对移动距离
Figure BDA0001235690040000113
进行约束，约束的目的是使顶点v移动后的新位置能达到最佳平滑效果；nk为体块Bk包含的邻域六面体网格个数；
Figure BDA0001235690040000114
为体块Bk质心到顶点v的向量，k表示体块序号。
当邻域六面体网格划分为体块后，公式(6)中每个体块的质心位置(Xk,Yk,Zk)、方向向量
Figure BDA0001235690040000115
包含的网格个数nk及网格对角线长度diag都是确定的，因此平滑系数γ的选择对平滑效果至关重要。
平滑系数γ取值范围为区间[0,1]，0表示移动距离
Figure BDA0001235690040000116
为0，即完全不进行平滑；1表示移动距离
Figure BDA0001235690040000117
受限于移动上限diag约束。当γ取0、1之间的小数时，则表明
Figure BDA0001235690040000118
受限于γ·diag约束，γ具体取值应根据顶点移动后的平滑效果评价指标确定。
平滑效果评价指标为当前顶点与其不同体块分界面上相邻顶点连线形成的夹角余弦值的平均值。
具体地说，在二维视角下，评价指标为当前顶点与其处于不同体块分界线上的相邻顶点连线形成的向量夹角余弦值，如图12所示，向量夹角越接近于180度，即夹角连线越接近于直线，余弦值越接近于-1，平滑效果越好。
在三维视角下，对于分界面上的平滑顶点，需分别计算顶点与XOY、XOZ、YOZ三个剖面上处于不同体块分界线上的相邻顶点连线形成的向量夹角余弦值，取其平均值作为评价指标，来评价整体的平滑效果，其值越接近于-1，平滑效果越好。
γ取值方法：平滑系数γ在区间[0,1]依据固定步长分别取值(例如γ分别取值0.01、0.02、0.03、…)，将不同取值分别代入上式计算移动距离
Figure BDA0001235690040000121
获取相应的顶点v移动后的新位置，并根据移动后的位置分别计算相应的平滑效果评价指标，从而选择评价指标最接近于-1的γ值为最终的平滑系数值。
4、当前顶点v在体块Bk作用下的移动距离向量
Figure BDA0001235690040000122
最终表示为：
Figure BDA0001235690040000123
式中，(xv,yv,zv)为顶点v坐标；(Xk,Yk,Zk)为体块Bk质心坐标；nk为体块Bk包含的六面体网格个数；diag为顶点v移动方向的邻域六面体网格对角线长度；γ为上一步骤确定的平滑系数值；
Figure BDA0001235690040000124
为体块Bk质心到顶点v的向量，
Figure BDA0001235690040000125
为向量
Figure BDA0001235690040000126
的模，
Figure BDA0001235690040000127
表示体块Bk与顶点v的距离。
步骤六：在获取各体块单独作用力下移动方向和距离的基础上，计算当前顶点在所有体块合力作用下沿x轴、y轴、z轴三个方向的移动分量。
当前顶点v在所有体块{Bk}合力作用下沿x轴、y轴、z轴三个方向的移动分量为：
Figure BDA0001235690040000128
式中，Δdx、Δdy、Δdz分别为x轴、y轴、z轴方向的移动分量；γ为平滑系数；diag为顶点v移动方向的邻域六面体网格对角线长度；(xv,yv,zv)为顶点v坐标；(Xk,Yk,Zk)为体块Bk质心坐标，k表示体块序号；nk为体块Bk包含的邻域六面体网格个数；
Figure BDA0001235690040000129
为体块Bk质心到顶点v的向量
Figure BDA00012356900400001210
模。
步骤七：根据x轴、y轴、z轴三个方向的移动分量计算当前顶点最终的平滑位置。
当前顶点v(xv,yv,zv)最终平滑移动的新位置(x′v,y′v,z′v)为：
Figure BDA00012356900400001211
式中，Δdx、Δdy、Δdz分别为x轴、y轴、z轴方向的移动分量；(xv,yv,zv)为顶点v的当前坐标位置；(x′v,y′v,z′v)为顶点v的新坐标位置。
下面通过实施例进行说明：
采用的网格为单层规则六面体网格，网格划分为63×63×1(即XOY平面为63×63个网格，Z轴方向为1层)，每个六面体网格大小为1m×1m×1m，平滑系数γ不同取值的平滑效果，及评价指标变化情况分别如图13～图16所示。
结合图16的平滑效果评价指标曲线变化情况，平滑系数γ为0表示没有进行平滑，如图13所示；随着γ值逐渐增大，平滑发生作用，评价指标也趋向-1；当平滑系数γ增大到某个值(本算例为0.275)时，曲线到达拐点位置，评价指标达到最为接近-1的值(本算例为-0.3013)，此时平滑效果最好，如图14所示；之后，随着γ取值继续增大，平滑性能下降，从图15(γ为0.6)可以看出，边界出现过度平滑现象。
对规则六面体网格，采用其他算例进行相同的计算，最佳平滑系数(曲线拐点位置)变化不大，γ基本在0.26～0.29之间变化，对比实际图形可发现，当区域边界出现较多的连续台阶状锯齿现象时(如图17所示)，拐点位置更接近于γ＝0.26，当区域边界出现较多的直角转弯现象时(如图18所示)，拐点位置更接近于γ＝0.29。实际应用时，对于规则六面体网格，可针对具体模型进行实际计算确定最佳的平滑系数γ值，也可简单地取区间[0.26,0.29]中值0.275为γ近似值。
对于不规则六面体网格，由于实际地层分界面和地质模型的多样性，上述的区间并不适用，平滑系数γ需在区间[0,1]依据固定步长分别取值，分别计算相应的平滑效果评价指标以获取最优的γ值，即前述步骤五所述的取值方法。但拐点现象是存在的，因为评价指标为平滑顶点与相邻顶点连线形成的夹角余弦值或余弦平均值，夹角总是在[0,360°]之间变化，其余弦曲线总是具有拐点位置(如图12、图16所示)。
上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。
提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (6)

1.一种基于邻域作用力的地质六面体网格平滑方法，其特征在于,包括以下步骤：
步骤一：根据三维六面体模型网格顶点坐标集合依次遍历三维模型每个顶点，及共用所述顶点的邻域六面体网格；对于每个顶点，进行下面步骤二至步骤七处理；
步骤二：对当前顶点的不同情况进行判断，以决定是否进行平滑处理，对需平滑的顶点进行下面步骤三至步骤七处理，并忽略不需要平滑处理的顶点；
步骤三：计算当前顶点各邻域六面体网格质心坐标；
步骤四：将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值分为多个体块，并根据各体块包含的邻域六面体网格质心坐标计算各个体块的质心坐标；
步骤五：根据当前顶点和各体块的坐标位置，以及各体块包含的邻域六面体网格个数，采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离；
步骤六：在获取各体块单独作用力下移动方向和距离的基础上，计算当前顶点在所有体块合力作用下沿x轴、y轴、z轴三个方向的移动分量；
步骤七：根据x轴、y轴、z轴三个方向的移动分量计算当前顶点最终的平滑位置；
所述步骤二中，对当前顶点的不同情况进行判断以决定是否进行平滑处理，包括以下五种情况：
(11)若顶点为三维模型边界面上的点，统计共用该顶点的邻域六面体网格在边界面上的属性值分布情况，若邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；若属性值不相同，在二维边界面上进行顶点平滑处理；
(12)若顶点为三维模型的角点或边点，则不进行平滑处理；
(13)若顶点为被人工标注过的断层、尖灭不需要进行平滑的特征点，则不进行平滑处理；
(14)若顶点在三维模型内部，且邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；
(15)若顶点在三维模型内部，且邻域六面体网格属性值不相同，则进行顶点平滑处理；
所述步骤五中，根据当前顶点和各体块的坐标位置，以及各体块包含的邻域六面体网格个数，采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离，具体步骤为：
(21)将当前顶点的每个体块看作磁铁的同极，各同极对顶点具有相斥的作用力，作用力采用基于重力模型的空间相互作用力模型进行度量，
设当前顶点为v，顶点v坐标位置为(xv,yv,zv)，体块Bk质心坐标为(Xk,Yk,Zk)，k表示体块序号，顶点v质量为Mv，体块Bk质量为Mk，体块Bk对顶点v的相斥作用力向量
Figure FDA0002361132600000021
为：
Figure FDA0002361132600000022
式中，G、α分别为引力常数和质量指数，这里取1；
Figure FDA0002361132600000023
为体块Bk中心到顶点v连线向量，体块中心取质心；β为距离指数，亦取1；
Figure FDA0002361132600000024
方向由
Figure FDA0002361132600000025
确定，方向向量
Figure FDA0002361132600000026
计算公式为：
Figure FDA0002361132600000027
(22)对于每个
Figure FDA0002361132600000028
取G＝α＝β＝1，由运动定律F＝ma可知顶点v在体块Bk作用力下的加速度向量
Figure FDA0002361132600000029
为：
Figure FDA00023611326000000210
Figure FDA00023611326000000211
方向同样由
Figure FDA00023611326000000212
确定；
(23)在不考虑初始速度和时间的前提下，加速度
Figure FDA00023611326000000213
代表了顶点v在作用力
Figure FDA00023611326000000214
作用下的移动偏移量，即移动距离
Figure FDA00023611326000000215
将每个六面体网格的质量视为1，则体块Bk的质量Mk等于体块Bk包含的网格个数nk，对于三维六面体模型，任一顶点v的邻域六面体网格个数最多为8个，当所有邻域六面体网格均属于同一体块时，则意味着邻域六面体网格属性值均相同、属性分布是均质的，不需要进行平滑处理，因此体块Bk的邻域六面体网格个数nk≤7，即体块Bk的质量Mk最大为7，设顶点移动方向的邻域六面体网格对角线长度diag为顶点移动距离的上限，则顶点v在最大质量和移动上限约束下的移动距离
Figure FDA00023611326000000216
表达为：
Figure FDA00023611326000000217
式中，γ为平滑系数，用于对移动距离
Figure FDA00023611326000000218
进行约束，约束的目的是使顶点v移动后的新位置能达到最佳平滑效果；nk为体块Bk包含的邻域六面体网格个数；
Figure FDA00023611326000000219
为体块Bk质心到顶点v的向量，k表示体块序号；
(24)当前顶点v在体块Bk作用下的移动距离向量
Figure FDA0002361132600000031
最终表示为：
Figure FDA0002361132600000032
式中，(xv,yv,zv)为顶点v坐标；(Xk,Yk,Zk)为体块Bk质心坐标；nk为体块Bk包含的邻域六面体网格个数；diag为顶点v移动方向的邻域六面体网格对角线长度；γ为上一步骤确定的平滑系数值；
Figure FDA0002361132600000033
为体块Bk质心到顶点v的向量，
Figure FDA0002361132600000034
为向量
Figure FDA0002361132600000035
的模，
Figure FDA0002361132600000036
Figure FDA0002361132600000037
表示体块Bk与顶点v的距离。
2.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤三中，计算当前顶点各邻域六面体网格质心坐标方法为：
对于每个邻域六面体网格，设该邻域六面体网格8个顶点的坐标为(xi,yi,zi)，其中i＝1,2,…,8，i为顶点序号，网格质心坐标为(xc,yc,zc)，质心坐标计算公式为：
Figure FDA0002361132600000038
3.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤四中，将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值可分为多个体块，根据各体块包含的邻域六面体网格质心坐标计算各体块质心坐标方法为：
设当前顶点的邻域六面体网格分为多个体块{Bk}，k表示体块序号，每个体块Bk的质心坐标为(Xk,Yk,Zk)，体块Bk包含的邻域六面体网格个数为nk，体块Bk内各邻域六面体网格的质心坐标为
Figure FDA0002361132600000039
其中i表示体块Bk包含的邻域六面体网格序号，i＝1,2,…,nk，则体块Bk质心坐标计算公式为：
Figure FDA00023611326000000310
Figure FDA00023611326000000311
4.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述平滑系数γ取值范围为区间[0,1]，具体取值根据平滑效果评价指标确定，平滑效果评价指标为当前顶点与其不同体块分界面上相邻顶点连线形成的夹角余弦值的平均值，γ取值方法：平滑系数γ在区间[0,1]依据固定步长分别取值，将不同取值分别代入上式计算移动距离
Figure FDA0002361132600000041
获取相应的顶点v移动后的新位置，并根据移动后的位置分别计算相应的平滑效果评价指标，从而选择评价指标最接近于-1的γ值为最终的平滑系数值。
5.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤六中，在获取各体块单独作用力下移动方向和距离的基础上，当前顶点v在所有体块{Bk}合力作用下沿x轴、y轴、z轴三个方向的移动分量为：
Figure FDA0002361132600000042
式中，Δdx、Δdy、Δdz分别为x轴、y轴、z轴方向的移动分量；γ为平滑系数；diag为顶点v移动方向的邻域六面体网格对角线长度；(xv,yv,zv)为顶点v坐标；(Xk,Yk,Zk)为体块Bk质心坐标，k表示体块序号；nk为体块Bk包含的邻域六面体网格个数；
Figure FDA0002361132600000043
为体块Bk质心到顶点v的向量
Figure FDA0002361132600000044
模。
6.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤七中，根据x轴、y轴、z轴三个方向的移动分量，当前顶点v(xv,yv,zv)最终平滑移动的新位置(x′v，y′v，z′v)为：
Figure FDA0002361132600000045
式中，Δdx、Δdy、Δdz分别为x轴、y轴、z轴方向的移动分量；(xv,yv,zv)为顶点v的当前坐标位置；(x′v，y′v，z′v)为顶点v的新坐标位置。
CN201710116290.XA 2017-03-01 2017-03-01 一种基于邻域作用力的地质六面体网格平滑方法 Expired - Fee Related CN106875487B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710116290.XA CN106875487B (zh) 2017-03-01 2017-03-01 一种基于邻域作用力的地质六面体网格平滑方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710116290.XA CN106875487B (zh) 2017-03-01 2017-03-01 一种基于邻域作用力的地质六面体网格平滑方法

Publications (2)

Publication Number Publication Date
CN106875487A CN106875487A (zh) 2017-06-20
CN106875487B true CN106875487B (zh) 2020-03-31

Family

ID=59168836

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710116290.XA Expired - Fee Related CN106875487B (zh) 2017-03-01 2017-03-01 一种基于邻域作用力的地质六面体网格平滑方法

Country Status (1)

Country Link
CN (1) CN106875487B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109783965B (zh) * 2019-01-25 2022-08-02 中国空气动力研究与发展中心计算空气动力研究所 一种结构网格自动分块加密方法
CN110335357B (zh) * 2019-05-23 2023-01-24 中国自然资源航空物探遥感中心 一种约束曲面多分辨率控制预处理方法
CN111445479B (zh) * 2020-03-20 2023-08-29 东软医疗系统股份有限公司 医学图像中感兴趣区域的分割方法及装置
CN114255188B (zh) * 2021-12-23 2022-06-24 中国矿业大学(北京) 三维地质岩性网格模型体绘制的边界线性平滑方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007094628A (ja) * 2005-09-28 2007-04-12 Japan Science & Technology Agency 粒子系を用いた陰関数曲面のポリゴン化方法
CN102612682A (zh) * 2009-11-12 2012-07-25 埃克森美孚上游研究公司 用于储层建模和模拟的方法和设备
CN102654923A (zh) * 2012-04-20 2012-09-05 浙江大学 一种由内向外的六面体网格生成方法
CN103914571A (zh) * 2014-04-25 2014-07-09 南京大学 一种基于网格分割的三维模型检索方法
CN105590335A (zh) * 2014-10-23 2016-05-18 富泰华工业（深圳）有限公司 点云网格精细化系统及方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007094628A (ja) * 2005-09-28 2007-04-12 Japan Science & Technology Agency 粒子系を用いた陰関数曲面のポリゴン化方法
CN102612682A (zh) * 2009-11-12 2012-07-25 埃克森美孚上游研究公司 用于储层建模和模拟的方法和设备
CN102654923A (zh) * 2012-04-20 2012-09-05 浙江大学 一种由内向外的六面体网格生成方法
CN103914571A (zh) * 2014-04-25 2014-07-09 南京大学 一种基于网格分割的三维模型检索方法
CN105590335A (zh) * 2014-10-23 2016-05-18 富泰华工业（深圳）有限公司 点云网格精细化系统及方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A data transformation method for visualizing the statistical information based on the grid;KIM M 等;《Journal of Korea Spatial Information Society》;20151231;第31-40页 *
三维网格模型平滑技术研究;陈琰;《中国优秀硕士学位论文全文数据库（信息科技辑）》;20050515;第I138-275页 *

Also Published As

Publication number Publication date
CN106875487A (zh) 2017-06-20

Similar Documents

Publication Publication Date Title
CN106875487B (zh) 一种基于邻域作用力的地质六面体网格平滑方法
CN104680573B (zh) 一种基于三角网格简化的纹理映射方法
US8537158B2 (en) Parallel triangle tessellation
CN106407605A (zh) 一种三维服装的粒子化计算机动态仿真方法
EP1840841A1 (en) Evolutionary direct manipulation of free form deformation representations for design optimisation
CN111382530A (zh) 学习用于推测实体cad特征的神经网络
CN113077553A (zh) 一种基于表面属性的三维模型分割方法
US11763048B2 (en) Computer simulation of physical fluids on a mesh in an arbitrary coordinate system
Wang et al. A triangular grid generation and optimization framework for the design of free-form gridshells
JP5916758B2 (ja) Ｇｐｕ上でのｃａｄモデルの直接的なレンダリング
EP1675068B1 (en) Evolutionary optimisation and free form deformation
CN106096130A (zh) 一种基于拉普拉斯坐标的不同材料衣物仿真及优化方法
CN104574508A (zh) 一种面向虚拟现实技术的多分辨率模型简化方法
Sanchez et al. Efficient evaluation of continuous signed distance to a polygonal mesh
Shchurova A methodology to design a 3D graphic editor for micro-modeling of fiber-reinforced composite parts
CN118154752A (zh) 基于三维离散数据对物体着色的方法、介质及设备
JP2881389B2 (ja) 自動メッシュ生成方法及びシステム
Krishnamurthy et al. Accurate GPU-accelerated surface integrals for moment computation
Rösch et al. Interactive visualization of implicit surfaces with singularities
CN104240301B (zh) 地质曲面重构方法及设备
CN107767458B (zh) 不规则三角网曲面几何拓扑一致分析方法及系统
CN101477707A (zh) 一种插值给定若干闭合曲线的曲面造型方法
CN101655832B (zh) 一种基于标量场梯度的物理变形方法
CN114972674B (zh) 一种面向切片数据的三维体积重构方法
Jabłoński et al. Real-time rendering of continuous levels of detail for sparse voxel octrees

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200331

Termination date: 20210301

CF01 Termination of patent right due to non-payment of annual fee