球体网格真三维地质模型漏洞的纬度带推扫填补算法*

王金鑫，曹泽宁，陈艺航，秦子龙，石 焱

(1.郑州大学 地球科学与技术学院，河南 郑州 450001；2.郑州大学 水利科学与工程学院，河南 郑州 450001)

地球剖分网格(Earth Tessellation Grid, ETG)是一种基于(椭)球体的对整个地球重力场进行递归剖分而形成的多分辨率的离散网格地球系统空间参考模型[1]，包括传统面向地球系统大型科学计算的地球物理和地球系统网格[2]以及面向时空大数据整合、集成、建模与应用的数字地球平台网格[3]两种。数字地球平台网格又包括全球离散网格(Discrete Global Grid System, DGGS)[4]和地球系统空间网格(Earth System Spatial Grid，ESSG)[5]两种。球体测地线八叉树网格(Sphere Geodesic Octree Grid，SGOG)[6-7]是ESSG模型的一种，首先进行球面大圆弧QTM四叉树剖分，然后进行径向二叉树剖分，再完成球体的三维网格分割。具有剖分规则简明、网格形体简单、排列规整、几何特征明晰的特点[6-7]。

与二维图件相比较，三维地质成果图件具有诸多优点，如基础资料利用的充分性与灵活性、成果的综合性、应用的直观性等[8]，并可动态地显示和查询地层模型的几何特征[9]，对地学工作者和地质工程师进行地质分析极为重要[10-14]。如何根据已有的离散地层数据，构建出完整的三维栅格地质模型，是采用钻孔数据建立三维地质模型的关键和难点。本文针对球体网格真三维地质建模存在的漏洞(横向与径向漏洞)修补问题，基于多边形填充原理和SGOG编码邻近搜索方法，设计了修补算法，取得了良好的效果。

1 填补算法技术路线

传统利用钻孔数据和体素(元)构建三维地质模型的方法，一般是找到钻孔数据的空间外包围盒(常为规则的几何形体)，然后进行空间剖分，逻辑上属于由外向内的方法。这种方法虽然不存在漏洞问题，但不利于复杂地质体的模型建立。球体网格概念的提出，基于网格的规则剖分和编码机制，可直接利用(或经适当加密的)离散采样点数据进行三维构模，逻辑上属于由内向外的方法。该方法十分有利于复杂地质体的三维建模，但由于采样点数量的稀疏性，常常会出现由于数据密度低而产生模型“漏洞”问题，包括横向表面漏洞(如图1(a)所示)和径向内部漏洞(如图1(b)所示)。

图1 地质真三维模型构建时产生的“漏洞”

1.1 算法总体思路

本文以SGOG网格为例，基于网格编码的拓扑推理提出纬度带推扫填补算法。在SGOG网格体系中，球面四叉树编码确定了网格的球面位置，径向二叉树编码确定了网格的径向位置。对位于同一球面位置的径向内部漏洞，可根据某位置上的最上、最下层网格直接进行漏洞编码填补。

横向上下表面漏洞的编码填补是本算法的关键。其主要思路为：首先，根据实际的应用需求，确定地质体三维建模的球面剖分层次(包括径向剖分层次)；其次，依据SGOG网格球面四叉树编码，进行网格的纬度带划分；再次，根据空间网格的球面拓扑关系以及填补规则，以地层钻孔数据圈定的空间范围作边界约束条件，遍历进行横向漏洞的编码填补；然后，对新增的横向球面网格进行径向上下面的高度位置内插并对其进行径向填充，进而完成整个地层模型的完整构建。对于地质体中多地层数据，重复进行上述流程，完成每层模型的构建，基于实际地理坐标进行集成即可。

1.2 算法关键技术

1.2.1 网格编码纬度带的划分

在一定的球面剖分精度下，模型上的球面网格所在的纬度位置(并非真正意义上的纬度)是该球面三角形编码的函数，而纬度位置可以用纬度带作为行数据表示。需要说明的是，SGOG网格的边为测地线大圆弧，与纬线小圆弧不重合。这里的纬度带仅作为填补的顺序导向，不是严格意义上的网格范围，填补通过编码实现，所以SGOG网格的边不影响算法的实现。

对于给定层级SGOG剖分，4个球面三角形(属于某个上层父三角形)被划分到2条纬度带中，分别记为第‘0’纬度带和第‘1’纬度带，如图2所示。那么当该级子剖分单元的编码(这里采用修正方向编码方法[7])为‘1’时，返回当前层级的纬度带标记值index=1；当该级子剖分单元的编码为‘0’、‘2’、‘3’时，返回当前层级的纬度带标记值index=0，即

图2 上、下三角形的纬度带划分

(1)

式中，index为计算当前层级的纬度带标记值；code[n]为SGOG球面编码序列；i为当前计算编码的剖分子级(从0开始计)。

这里的纬度带是根据球面剖分层次从低纬度到高纬度按照升序的方式划分的，记低纬度带的编码为‘0’，高维度带编码为‘1’(北半球的情况)。当上一层级对应的球面三角为上三角时(即上一层级编码为‘1’)，编码为‘1’的子三角处于高纬度位置(如图2左图所示)；当上一层级对应的球面三角为下三角时(即上一层级编码为‘0’)，编码为‘1’的子三角处于低纬度位置(如图2右图所示)；故此时存在纬度带划分的正反序问题，必须进行相应的调整。当计算的上一层级球面三角为下三角时，记编码为‘1’时返回的标记值为‘0’，编码为‘0’、‘2’、‘3’时返回的标记值为‘1’，故对式(1)修改得：

(2)

通过计算SGOG网格每一层级编码的纬度带标记值，可得到一条长度(位数)与网格编码相同的纬度带标记值序列index[n]，对该序列每一位数按照如下公式计算，即可得到该SGOG球面编码所在的纬度带值，即

(3)

式中，Nlat为计算得到的纬度带编号；index[n]为纬度带标记值序列；n为返回的标记值长度；i为当前标记值位置(从0开始计)。

当计算位于不同半球的SGOG球面四叉树编码的纬度带时，SGOG球面四叉树编码所处的位置只会因为南北半球的不同而导致东西位置的改变，并不会改变编码所在的南北位置，故不会影响纬度带编号的计算，如图3所示，所以此方法适用于全球SGOG球面四叉树编码。

图3 南北半球在2级剖分层次下的各编码纬度带位置示意图

1.2.2 基于纬度带的SGOG编码填补

根据上述方法，取出每个位于同一纬度带上的所有内部和边界的SGOG球面四叉树编码，从边界球面四叉树编码出发自东向西的顺序对同一纬度带上的编码进行邻近计算(如图4所示，方向m至n即从东向西)，若发现邻近网格不存在，则对该网格进行填补。修正方向编码具有子三角形相对于其二级父三角形位置固定的特点，据此设计其邻近搜索算法(另文著述)。

图4 基于编码的填补情况示意图

图中虚线部分为模型区域拟定的边界，蓝色部分为存在的网格，则共有5种编码填补情况：

情况1：编码a左、右邻近编码均存在，则不需要进行编码填补；

情况2：编码b仅存在左邻近或仅存在右邻近或不存在左右邻近，对缺少的编码进行编码填补，并对填补后的编码继续做左(右)邻近判断，直到触碰到边界为止；

情况3：编码c位于边界上且存在另一边的邻近编码，则不需要进行编码填补；

情况4：编码d位于边界上但不存在另一边的邻近编码，对缺少的编码进行编码填补，并对填补后的编码继续做左(右)邻近判断，直到触碰到边界为止；

情况5：编码e已超出边界范围，删除超出边界范围的编码。

根据以上规则，对同一纬度带的编码填补工作流程如下：(1)记在同一纬度带下所有存在的网格编码集合为D，根据SGOG网格邻近关系将集合D中各元素Di进行位置顺序的排列；(2)获取在该纬度带上边界的位置信息，规定填补方向(由东向西或由西向东)，以由东向西为例，在填补时每个纬度带上可能存在多个边界(总数为偶数)，则根据从东向西的顺序分别对边界赋予“进”与“出”的标签；(3)对“进-出”标签内的网格编码Di根据图4中的前4种情况进行填补，对“出-进”标签内的网格编码Di根据图4中的第五种情况，认为其已超出给定的边界范围外则删去。

递归上述方法完成对所有纬度带内编码的填补即完成模型横向模型填补，对横向网格还需进行高度位置的内插，进而完成整个模型的空间重构。

另外，对于存在空洞(真实存在的洞)的复杂地质体，必须首先给出空洞的数据边界(空洞部分的地层数据)，然后将其与球体网格匹配(栅格化)，得出内边界后再进行填补。

综上，填补算法流程图如图5所示。

图5 基于SGOG方向修正编码填补算法流程图

2 实验与分析

2.1 填补实例

为了验证本文填补算法的可行性，选取郑州市航空港地区某地质层三维构模作为研究对象，构建横向17层剖分层次、径向25层SGOG剖分层次的地质模型(建模前钻孔数据经过适当加密)，模型填补效果如图6所示。共填补球面横向漏洞5 078个，其约占研究区域总面积的2.5%(研究区域总面积约643.5 km2，漏洞面积约为16.4 km2)。

图6 航空港区模型填补前后比较图

由于上述模型的边界较为简单，故截取地形边界较为复杂的长沙市长沙县(山地丘陵地形)的DEM数据(来自于共享网站：http：//srtm.csi.cgiar.org，将DEM数据随机删去20%来模拟漏洞)作为地层数据进行地表真三维建模来验证算法有效性(模型剖分参数为：SGOG横向16层，径向25层)。研究区域表层横向网格总数量129 830个，其中人为删除(模拟漏洞)数量21 638个，使用本文算法共填补漏洞网格21 312个，填补产生误差小于2%，满足三维模型的可视化及空间分析等需求，填补效果如图7所示。

图7 长沙市长沙县模型填补前后比较图

根据图6、图7对比看出，无论地形复杂与否，本算法均达到了良好的修补效果。

对航空港地区的地质模型，给定2个空洞的边界，其修补效果如图8所示，可以看出本文算法在避开固有空洞的前提下对模型中的漏洞填补效果良好，并未出现对固有空洞造成错误填补的情况。

图8 存在固有空洞地质模型的漏洞填补示意图

2.2 算法的时间复杂度测算

实验环境：Intel(R) Xeon(R) Silver 4110 CPU @ 2.10GHz，RAM 64.0GB，GPU NVIDIA Quadro P2000，Win10专业版，使用MySQL数据库对编码数据的储存和调用。实验的基本原理为：选取12个存在漏洞的地质模型，对这些模型进行漏洞填补，记录每个模型填补的漏洞个数与填补耗费的时间，绘制填补数量与时间关系曲线，如图9所示。可以看出，网格填补数量(个)与工作时间(s)呈线性正相关，本文填补算法的时间复杂度为T(n) =O(f(n))。

图9 填补算法效率趋势图

3 结 语

本文以SGOG网格为例，就利用钻孔数据生成真三维地层模型的漏洞填补问题，提出了纬度带推扫填补算法。该算法充分利用球体网格的剖分规则及其空间拓扑关系，以网格编码为核心，实现地质实体的遍历填补。研究表明：(1)利用规则球体网格建立真三维地质实体模型是一种可行的方法，尤其是对大区域和全球尺度地质空间，具有较大的精度优势(考虑地球曲率，无投影变形)；(2)SGOG网格兼具TIN与Grid的特点，在真三维地表(地壳)可视化表达与空间分析中具有较大应用潜力；(3)本文提出的方法对所有球体网格具有普适性，填补效果良好。本研究对新一代数字地球平台建设及时空大数据管理与应用具有参考价值。