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

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个顶点的坐标为(x_i,y_i,z_i)，其中i＝1,2,…,8，i为顶点序号，网格质心坐标为(x_c,y_c,z_c)，质心坐标计算公式为：

进一步的，所述步骤四中，将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值可分为多个体块，根据各体块包含的邻域六面体网格质心坐标计算各体块质心坐标方法为：

设当前顶点的邻域六面体网格可分为多个体块{B_k}，k表示体块序号，每个体块B_k的质心坐标为(X_k,Y_k,Z_k)，体块B_k包含的邻域六面体网格个数为n_k，体块B_k内各邻域六面体网格的质心坐标为其中i表示体块B_k包含的邻域六面体网格序号，i＝1,2,…,n_k，则体块B_k质心坐标计算公式为：

进一步的，所述步骤五中，根据当前顶点和各体块的坐标位置，以及各体块包含的邻域六面体网格个数，采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离，具体步骤为：

1、将当前顶点的每个体块看作磁铁的“同极”，各同极对顶点具有相斥的作用力，作用力采用基于重力模型的空间相互作用力SIM模型(Spatial Interaction Model)进行度量。

设当前顶点为v，顶点v坐标位置为(x_v,y_v,z_v)，体块B_k质心坐标为(X_k,Y_k,Z_k)，k表示体块序号，顶点v质量为M_v，体块B_k质量为M_k，体块B_k对顶点v的相斥作用力向量为：

式中，G、α分别为引力常数和质量指数，这里取1；为体块B_k中心到顶点v连线向量，体块中心取质心；β为距离指数，亦取1。

方向由确定，方向向量计算公式为：

2、对于每个取G＝α＝β＝1，由运动定律F＝ma可知顶点v在体块B_k作用力下的加速度向量为：

方向同样由确定。

3、在不考虑初始速度和时间的前提下，加速度代表了顶点v在作用力作用下的移动偏移量，即移动距离将每个六面体网格的质量视为1，则体块B_k的质量M_k等于体块B_k包含的网格个数n_k。对于三维六面体模型，任一顶点v的邻域六面体网格个数最多为8个，当所有邻域六面体网格均属于同一体块时，则意味着邻域六面体网格属性值均相同、属性分布是均质的，不需要进行平滑处理，因此体块B_k的邻域六面体网格个数n_k≤7，即体块B_k的质量M_k最大为7。设顶点移动方向的邻域六面体网格对角线长度diag为顶点移动距离的上限，则顶点v在最大质量和移动上限约束下的移动距离可表达为：

式中，γ为平滑系数，用于对移动距离进行约束，约束的目的是使顶点v移动后的新位置能达到最佳平滑效果；n_k为体块B_k包含的邻域六面体网格个数；为体块B_k质心到顶点v的向量。

平滑系数γ取值范围为区间[0,1]，具体取值根据平滑效果评价指标确定。

平滑效果评价指标为当前顶点与其不同体块分界面上相邻顶点连线形成的夹角余弦值的平均值。具体地说，在二维视角下，评价指标为当前顶点与其处于不同体块分界线上的相邻顶点连线形成的夹角余弦平均值；在三维视角下，评价指标为当前顶点与XOY、XOZ、YOZ三个剖面上处于不同体块分界线上的相邻顶点连线形成的夹角余弦平均值。余弦值或余弦平均值越接近于-1(即夹角越接近于180度)，平滑效果越好。

γ取值方法：平滑系数γ在区间[0,1]依据固定步长分别取值(例如γ分别取值0.01、0.02、0.03、…)，将不同取值分别代入上式计算移动距离获取相应的顶点v移动后的新位置，并根据移动后的位置分别计算相应的平滑效果评价指标，从而选择评价指标最接近于-1的γ值为最终的平滑系数值。

4、当前顶点v在体块B_k作用下的移动距离向量最终表示为：

式中，(x_v,y_v,z_v)为顶点v坐标；(X_k,Y_k,Z_k)为体块B_k质心坐标；n_k为体块B_k包含的邻域六面体网格个数；diag为顶点v移动方向的邻域六面体网格对角线长度；γ为上一步骤确定的平滑系数值；为体块B_k质心到顶点v的向量，为向量的模，表示体块B_k与顶点v的距离。

进一步的，所述步骤六中，在获取各体块单独作用力下移动方向和距离的基础上，当前顶点v在所有体块{B_k}合力作用下沿x轴、y轴、z轴三个方向的移动分量为：

式中，Δd_x、Δd_y、Δd_z分别为x轴、y轴、z轴方向的移动分量；γ为平滑系数；diag为顶点v移动方向的邻域六面体网格对角线长度；(x_v,y_v,z_v)为顶点v坐标；(X_k,Y_k,Z_k)为体块B_k质心坐标；n_k为体块B_k包含的邻域六面体网格个数；为体块B_k质心到顶点v的向量模。

进一步的，所述步骤七中，根据x轴、y轴、z轴三个方向的移动分量，当前顶点v(x_v,y_v,z_v)最终平滑移动的新位置(x′_v,y′_v,z′_v)为：

式中，Δd_x、Δd_y、Δd_z分别为x轴、y轴、z轴方向的移动分量；(x_v,y_v,z_v)为顶点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算法中，网格顶点v_i平滑移动到与其共面邻接顶点v_j的平均位置，其移动偏移量其中n是v_i邻接顶点的个数，w_j用于度量邻接顶点v_j在平滑移动中的权重大小，最简单的形式是采用相同权重，即w_j＝1/n。

在本发明方法中，对上述平滑思路进行扩展，以顶点周围共用该顶点的八个邻域六面体网格组成平滑窗口，将窗口内邻域六面体网格划分为不同的体块，不同体块对顶点相斥作用力的大小度量了顶点平滑移动的不同权重，所有体块作用力的合力决定了顶点最终移动的距离和方向。

利用发明方法对三维区域边界的锯齿进行平滑处理，目的是降低边界表达的噪声，使其符合实际的地质现象。对于有些顶点，例如断层、尖灭等，是不需要进行平滑处理的。当顶点需进行平滑处理时，其邻域的各体块总是将顶点外推，目的是使体块尖锐、锯齿状的边界面和边界点趋向光滑，这与地层分界面的连续性是一致的，只不过质量较大的体块将顶点外推的作用力较大，质量较小的体块将顶点外推的作用力较小，不同方向、大小的作用力形成合力最终将顶点推到相对平衡的位置，既降低了边界的锯齿效应，又保持了网格的稳定性。

发明方法的主要思路可概括为：将共用同一顶点的八个邻域六面体网格体划分为不同的体块，每个体块包含具有相同属性值的网格。将每个体块看作磁铁的“同极”，同极之间具有相斥的作用力，顶点移动的方向为体块作用力的方向，移动的距离以作用力的大小作为权重，作用力大则移动距离大，作用力小则移动距离小。所有体块合力的方向和大小决定了顶点最终移动的方向和距离。基于邻域作用力的地质六面体网格平滑方法就是根据地质属性分布的一般规律，在保证六面体网格之间三维拓扑关系不变的同时对网格顶点进行移动平滑。

如图2所示，一种基于邻域作用力的地质六面体网格平滑方法，包括以下步骤：

步骤一：根据三维六面体模型网格顶点坐标集合依次遍历三维模型每个顶点，及共用该顶点的邻域六面体网格；对于每个顶点，进行下面步骤二至步骤七处理。

步骤二：对当前顶点的不同情况进行判断以决定是否进行平滑处理，对需平滑的顶点进行下面步骤三至步骤七处理。

对当前顶点的不同情况进行判断以决定是否进行平滑处理，包括以下五种情况(如图3所示)：

1、若顶点为三维模型边界面上的点，统计共用该顶点的邻域六面体网格在边界面上的属性值分布情况，若邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；若属性值不相同，在二维边界面上采用下面步骤三至步骤七进行顶点平滑处理；

2、若顶点为三维模型的角点或边点，则不进行平滑处理；

3、若顶点为被人工标注过的断层、尖灭等不需要进行平滑的特征点，则不进行平滑处理；

断层、地层尖灭等区域覆盖的不需要进行平滑的特征顶点，应在平滑处理前，由地质工程师进行交互操作，圈定或输入一些特征点或特征区域，避免平滑时被破坏掉。

4、若顶点在三维模型内部，且邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；

5、若顶点在三维模型内部，且邻域六面体网格属性值不相同，则采用下面步骤三至步骤七进行顶点平滑处理。

步骤三：计算当前顶点各邻域六面体网格质心坐标。

对于每个邻域六面体网格，设该六面体网格8个顶点的坐标为(x_i,y_i,z_i)，其中i＝1,2,…,8，i为顶点序号，网格质心坐标为(x_c,y_c,z_c)，质心坐标计算公式为：

步骤四：将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值可分为多个体块，并计算各体块的质心坐标。

以图4示例，设顶点v周围八个邻域六面体网格的离散属性值代表不同的沉积微相类型，如用数值1、2、3分别代表坝主体、坝缘、前三角洲类型，相同的属性值指示相同的沉积微相类型，据此可将八个邻域六面体网格划分为三个体块B₁、B₂、B₃，分别代表上述三种微相类型。图4为了便于展示，将顶点v与邻域六面体网格拉开了距离。

设当前顶点的邻域六面体网格可分为多个体块{B_k}，每个体块B_k的质心坐标为(X_k,Y_k,Z_k)，体块B_k包含的邻域六面体网格个数为n_k，体块B_k内各邻域六面体网格的质心坐标为其中i表示体块B_k包含的邻域六面体网格序号，i＝1,2,…,n_k，则体块B_k质心坐标计算公式为：

步骤五：采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离。

1、将当前顶点的每个体块看作磁铁的“同极”，各同极对顶点具有相斥的作用力，作用力采用基于重力模型的空间相互作用力SIM模型(Spatial Interaction Model)进行度量。

设当前顶点为v，顶点v坐标位置为(x_v,y_v,z_v)，体块B_k质心坐标为(X_k,Y_k,Z_k)，k表示体块序号，顶点v质量为M_v，体块B_k质量为M_k，体块B_k对顶点v的相斥作用力向量为：

式中，G、α分别为引力常数和质量指数，这里取1；为体块B_k中心到顶点v连线向量，体块中心取质心；β为距离指数，亦取1。

方向由确定。图4中三个体块对顶点v的作用力大小分别由上式(3)确定，方向分别由向量确定，则合力实际上是三个作用力向量之和。

方向向量计算公式为：

2、对于每个取G＝α＝β＝1，由运动定律F＝ma可知顶点v在体块B_k作用力下的加速度向量为：

方向同样由确定。在不考虑初始速度和时间的前提下，加速度代表了顶点v在作用力作用下的移动偏移量(即移动距离)则在合力作用下总的移动偏移量为各向量之和。

公式(5)说明大小与体块B_k的质量M_k成正比，与距离成反比。若将每个六面体网格的质量简单地视为1，则大小实际上与体块B_k包含的六面体网格个数n_k成正比，表明当体块包含的网格个数较多时，对顶点的作用力较大，顶点移动距离较长，如图4中的体块B₃；反之，则对顶点的作用力较小，顶点移动距离较短，如图4中的体块B₂。

下面对网格体块分布的几种具体情况进行图示分析：

(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，则体块B_k的质量M_k等于体块B_k包含的网格个数n_k。对于三维六面体模型，任一顶点v的邻域六面体网格个数最多为8个，当所有邻域六面体网格均属于同一体块时，则意味着邻域六面体网格属性值均相同、属性分布是均质的，不需要进行平滑处理，即前述步骤二中第4种情况，因此，体块B_k的邻域六面体网格个数n_k≤7，即体块B_k的质量M_k最大为7。设顶点移动方向的邻域六面体网格对角线长度diag为顶点移动距离的上限，则顶点v在体块B_k作用下基于最大质量和移动上限约束的移动距离可表达为：

式中，γ为平滑系数，用于对移动距离进行约束，约束的目的是使顶点v移动后的新位置能达到最佳平滑效果；n_k为体块B_k包含的邻域六面体网格个数；为体块B_k质心到顶点v的向量，k表示体块序号。

当邻域六面体网格划分为体块后，公式(6)中每个体块的质心位置(X_k,Y_k,Z_k)、方向向量包含的网格个数n_k及网格对角线长度diag都是确定的，因此平滑系数γ的选择对平滑效果至关重要。

平滑系数γ取值范围为区间[0,1]，0表示移动距离为0，即完全不进行平滑；1表示移动距离受限于移动上限diag约束。当γ取0、1之间的小数时，则表明受限于γ·diag约束，γ具体取值应根据顶点移动后的平滑效果评价指标确定。

平滑效果评价指标为当前顶点与其不同体块分界面上相邻顶点连线形成的夹角余弦值的平均值。

具体地说，在二维视角下，评价指标为当前顶点与其处于不同体块分界线上的相邻顶点连线形成的向量夹角余弦值，如图12所示，向量夹角越接近于180度，即夹角连线越接近于直线，余弦值越接近于-1，平滑效果越好。

在三维视角下，对于分界面上的平滑顶点，需分别计算顶点与XOY、XOZ、YOZ三个剖面上处于不同体块分界线上的相邻顶点连线形成的向量夹角余弦值，取其平均值作为评价指标，来评价整体的平滑效果，其值越接近于-1，平滑效果越好。

γ取值方法：平滑系数γ在区间[0,1]依据固定步长分别取值(例如γ分别取值0.01、0.02、0.03、…)，将不同取值分别代入上式计算移动距离获取相应的顶点v移动后的新位置，并根据移动后的位置分别计算相应的平滑效果评价指标，从而选择评价指标最接近于-1的γ值为最终的平滑系数值。

4、当前顶点v在体块B_k作用下的移动距离向量最终表示为：

式中，(x_v,y_v,z_v)为顶点v坐标；(X_k,Y_k,Z_k)为体块B_k质心坐标；n_k为体块B_k包含的六面体网格个数；diag为顶点v移动方向的邻域六面体网格对角线长度；γ为上一步骤确定的平滑系数值；为体块B_k质心到顶点v的向量，为向量的模，表示体块B_k与顶点v的距离。

步骤六：在获取各体块单独作用力下移动方向和距离的基础上，计算当前顶点在所有体块合力作用下沿x轴、y轴、z轴三个方向的移动分量。

当前顶点v在所有体块{B_k}合力作用下沿x轴、y轴、z轴三个方向的移动分量为：

式中，Δd_x、Δd_y、Δd_z分别为x轴、y轴、z轴方向的移动分量；γ为平滑系数；diag为顶点v移动方向的邻域六面体网格对角线长度；(x_v,y_v,z_v)为顶点v坐标；(X_k,Y_k,Z_k)为体块B_k质心坐标，k表示体块序号；n_k为体块B_k包含的邻域六面体网格个数；为体块B_k质心到顶点v的向量模。

步骤七：根据x轴、y轴、z轴三个方向的移动分量计算当前顶点最终的平滑位置。

当前顶点v(x_v,y_v,z_v)最终平滑移动的新位置(x′_v,y′_v,z′_v)为：

式中，Δd_x、Δd_y、Δd_z分别为x轴、y轴、z轴方向的移动分量；(x_v,y_v,z_v)为顶点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 (8)

1.一种基于邻域作用力的地质六面体网格平滑方法，其特征在于,包括以下步骤：

步骤一：根据三维六面体模型网格顶点坐标集合依次遍历三维模型每个顶点，及共用所述顶点的邻域六面体网格；对于每个顶点，进行下面步骤二至步骤七处理；

步骤二：对当前顶点的不同情况进行判断，以决定是否进行平滑处理，对需平滑的顶点进行下面步骤三至步骤七处理，并忽略不需要平滑处理的顶点；

步骤三：计算当前顶点各邻域六面体网格质心坐标；

步骤四：将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值分为多个体块，并根据各体块包含的邻域六面体网格质心坐标计算各个体块的质心坐标；

步骤五：根据当前顶点和各体块的坐标位置，以及各体块包含的邻域六面体网格个数，采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离；

步骤六：在获取各体块单独作用力下移动方向和距离的基础上，计算当前顶点在所有体块合力作用下沿x轴、y轴、z轴三个方向的移动分量；

步骤七：根据x轴、y轴、z轴三个方向的移动分量计算当前顶点最终的平滑位置。

2.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤二中，对当前顶点的不同情况进行判断以决定是否进行平滑处理，包括以下五种情况：

(1)若顶点为三维模型边界面上的点，统计共用该顶点的邻域六面体网格在边界面上的属性值分布情况，若邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；若属性值不相同，在二维边界面上进行顶点平滑处理；

(2)若顶点为三维模型的角点或边点，则不进行平滑处理；

(3)若顶点为被人工标注过的断层、尖灭不需要进行平滑的特征点，则不进行平滑处理；

(4)若顶点在三维模型内部，且邻域六面体网格属性值均相同，则认为属性分布是均质的，不进行平滑处理；

(5)若顶点在三维模型内部，且邻域六面体网格属性值不相同，则进行顶点平滑处理。

3.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤三中，计算当前顶点各邻域六面体网格质心坐标方法为：

对于每个邻域六面体网格，设该邻域六面体网格8个顶点的坐标为(x_i,y_i,z_i)，其中i＝1,2,…,8，i为顶点序号，网格质心坐标为(x_c,y_c,z_c)，质心坐标计算公式为：

4.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤四中，将属性值相同的邻域六面体网格看作同一体块，则当前顶点的邻域六面体网格根据其属性值可分为多个体块，根据各体块包含的邻域六面体网格质心坐标计算各体块质心坐标方法为：

设当前顶点的邻域六面体网格分为多个体块{B_k}，k表示体块序号，每个体块B_k的质心坐标为(X_k,Y_k,Z_k)，体块B_k包含的邻域六面体网格个数为n_k，体块B_k内各邻域六面体网格的质心坐标为其中i表示体块B_k包含的邻域六面体网格序号，i＝1,2,…,n_k，则体块B_k质心坐标计算公式为：

5.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤五中，根据当前顶点和各体块的坐标位置，以及各体块包含的邻域六面体网格个数，采用基于重力模型的空间相互作用力模型计算当前顶点在各体块单独作用力下的移动方向及距离，具体步骤为：

(1)将当前顶点的每个体块看作磁铁的同极，各同极对顶点具有相斥的作用力，作用力采用基于重力模型的空间相互作用力SIM模型(Spatial Interaction Model)进行度量，

设当前顶点为v，顶点v坐标位置为(x_v,y_v,z_v)，体块B_k质心坐标为(X_k,Y_k,Z_k)，k表示体块序号，顶点v质量为M_v，体块B_k质量为M_k，体块B_k对顶点v的相斥作用力向量为：

式中，G、α分别为引力常数和质量指数，这里取1；为体块B_k中心到顶点v连线向量，体块中心取质心；β为距离指数，亦取1；

方向由确定，方向向量计算公式为：

(2)对于每个取G＝α＝β＝1，由运动定律F＝ma可知顶点v在体块B_k作用力下的加速度向量为：

方向同样由确定；

(3)在不考虑初始速度和时间的前提下，加速度代表了顶点v在作用力作用下的移动偏移量，即移动距离将每个六面体网格的质量视为1，则体块B_k的质量M_k等于体块B_k包含的网格个数n_k，对于三维六面体模型，任一顶点v的邻域六面体网格个数最多为8个，当所有邻域六面体网格均属于同一体块时，则意味着邻域六面体网格属性值均相同、属性分布是均质的，不需要进行平滑处理，因此体块B_k的邻域六面体网格个数n_k≤7，即体块B_k的质量M_k最大为7，设顶点移动方向的邻域六面体网格对角线长度diag为顶点移动距离的上限，则顶点v在最大质量和移动上限约束下的移动距离表达为：

式中，γ为平滑系数，用于对移动距离进行约束，约束的目的是使顶点v移动后的新位置能达到最佳平滑效果；n_k为体块B_k包含的邻域六面体网格个数；为体块B_k质心到顶点v的向量，k表示体块序号；

(4)当前顶点v在体块B_k作用下的移动距离向量最终表示为：

式中，(x_v,y_v,z_v)为顶点v坐标；(X_k,Y_k,Z_k)为体块B_k质心坐标；n_k为体块B_k包含的邻域六面体网格个数；diag为顶点v移动方向的邻域六面体网格对角线长度；γ为上一步骤确定的平滑系数值；为体块B_k质心到顶点v的向量，为向量的模，表示体块B_k与顶点v的距离。

6.根据权利要求5所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述平滑系数γ取值范围为区间[0,1]，具体取值根据平滑效果评价指标确定，平滑效果评价指标为当前顶点与其不同体块分界面上相邻顶点连线形成的夹角余弦值的平均值，γ取值方法：平滑系数γ在区间[0,1]依据固定步长分别取值，将不同取值分别代入上式计算移动距离获取相应的顶点v移动后的新位置，并根据移动后的位置分别计算相应的平滑效果评价指标，从而选择评价指标最接近于-1的γ值为最终的平滑系数值。

7.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤六中，在获取各体块单独作用力下移动方向和距离的基础上，当前顶点v在所有体块{B_k}合力作用下沿x轴、y轴、z轴三个方向的移动分量为：

式中，Δd_x、Δd_y、Δd_z分别为x轴、y轴、z轴方向的移动分量；γ为平滑系数；diag为顶点v移动方向的邻域六面体网格对角线长度；(x_v,y_v,z_v)为顶点v坐标；(X_k,Y_k,Z_k)为体块B_k质心坐标，k表示体块序号；n_k为体块B_k包含的邻域六面体网格个数；为体块B_k质心到顶点v的向量模。

8.根据权利要求1所述的基于邻域作用力的地质六面体网格平滑方法，其特征在于：所述步骤七中，根据x轴、y轴、z轴三个方向的移动分量，当前顶点v(x_v,y_v,z_v)最终平滑移动的新位置(x′_v,y′_v,z′_v)为：

式中，Δd_x、Δd_y、Δd_z分别为x轴、y轴、z轴方向的移动分量；(x_v,y_v,z_v)为顶点v的当前坐标位置；(x′_v,y′_v,z′_v)为顶点v的新坐标位置。