范帅博，童晓冲，雷 毅

（信息工程大学 地理空间信息学院，河南 郑州 450000）

0 引 言

全球离散格网系统采用特定方法将地球均匀离散化，形成无缝无叠的多分辨率格网层次结构，是一种新型的空间数据组织、管理与应用模型[1-2]，其中基于正多面体剖分的离散格网是广受关注的格网系统之一[3-4]，本文的研究也正是围绕该类格网系统展开。近年来，针对不同类型的正多面体格网生产方法，国内外诸多研究都给出了不同的方法，如针对三角形QTM格网的生成方法[5-6]；针对菱形剖分的层次格网生产算法[7-8]；六边形格网的正多面体生成算法[8]。上述研究成果中，多面体格网的生成方式中关于不同层级格网单元的构造，归纳起来主要有两种思路：一种是采用逐层递归的方式进行剖分[9]；另一种是采用分层逐单元排列的方式进行生成。虽然这些方法可以生成各类不同的三角形、四边形和六边形格网，但是一般而言，在这些格网的生成过程中，存在下面两类问题：

1）需要根据目标网格的形态（三角形、四边形、六边形）、孔径（3孔或4孔）[9]、剖分频率[1]的质因数、对偶性、对心偏心性等因素[8]定制相应的格网生成算法；

2）在采用递归方法进行格网生成的过程中，高层次格网的生成需要设计有效分布式算法，算法复杂度较高。

针对上述问题，通过对现有不同类型格网单元，以及不同层级格网的分析，本文提出了一种基于单元复制方法的通用化全球离散格网生成方法，该方法采用“简单单元复制+有效区域限制”的方式，能够解决不同形态格网、不同层级格网的统一化生成问题。为了说明问题，论文以正二十面体上六边形格网剖分为例展开讨论，其他类型的格网生成可以采用通用的方式。全球离散格网的生成方法采用文献[10]中图1中的方式，本文设计算法主要集中在多尺度格网系统的构建这一过程中。

1 基本思想

由于正二十面体是中心对称图形，因此仅考虑一个三角面上格网剖分情况，我们定义为单一三角面格网系统，首先，固定网格的基本剖分形态，分析现有各种剖分类型的单一三角面格网系统，其本质区别在于单一三角面格网系统中所包含的剖分网格的数目和其相对于三角面的位置姿态有所不同。因此，从另一角度而言，我们固定网格的基本剖分形态并将其扩充至无限平面中使之无缝无叠，我们称之为“基扩充平面格网系统”。下面就算法的基本思想进行介绍：

在基扩充平面格网系统中，建立适当的离散格网斜坐标系，首先固定某一初始正三角形的顶点坐标，称该三角形为初始单一三角形；通过对初始单一三角形的边界顶点坐标进行适当的旋转、平移、缩放变换，使之与基扩充平面格网系统中的部分剖分网格单元适当地组合起来，从而生成各种剖分类型任意层次的单一三角面格网系统，这样不同组合变换参数与特定剖分类型和层次的单一三角面格网系统建立一一对应的数学关系。通过改变网格的剖分形态和起算变换参数从而生成三角形、四边形、六边形任意剖分类型任意剖分层次的单一三角面格网系统。

结合上述的思路，算法分为下面几个步骤：

1）建立基扩充格网斜坐标系；

2）确定待计算剖分格网系统的类型；

3）确定不同层级格网的有效控制边界；

4）计算不同层级格网单元节点坐标；

5）生成单一三角面格网系统坐标；

6）建立不同层级单元节点的关联关系。

具体的算法技术路线图，如图1所示。

图1 算法技术路线图Fig.1 Algorithm technology roadmap

2 算法流程

2.1 建立基扩充格网斜坐标系

基扩充格网坐标系是建立在无限平面离散格网中的斜坐标系，主要用于确定单一三角面格网系统中各种剖分形态、描述剖分类型规则排列的网格单元的中心位置，其坐标值都是整数，下文统一称格网坐标。为了便于研究，本文以六边形格网为例（后文中提到的格网皆指六边形格网），定义坐标原点位于无限平面离散格网中任一格网单元中心，水平相邻格网中心连线作为坐标系I轴，坐标系J轴按照I轴逆时针方向旋转120°方向得到。

另外，还需要建立辅助的空间直角坐标系，该坐标系的原点与X轴分别与基扩充格网斜坐标系的原点与I轴重合，Y轴按照X轴逆时针方向旋转90°得到，Z轴遵循右手螺旋法则，下文统一称直角坐标。

在基扩充格网斜坐标系中，定义格网单元的半径是R，初始单一三角形的顶点格网坐标分别是A0(0,0)，B0(1,0)，C0(1,1)。对于格网坐标(i,j)的任一单元的7个节点中，中心节点用O(i,j)标识，6个边界节点从右上角按逆时针方向依次用A(i,j)、B(i,j)、C(i,j)、D(i,j)、E(i,j)、F(i,j)标识。具体情况如图2所示。

图2 基扩充格网斜坐标系和空间直角坐标系的定义Fig.2 The def i nition of extended grid skew coordinate system and space cartesian coordinate system

通过基扩充格网斜坐标系和空间直角坐标系之间的几何变换关系，得出基扩充格网斜坐标系下格网坐标是(i,j)的任意格网单元7个节点的空间直角坐标分别如式（1）：

2.2 确定待计算剖分格网系统的类型

本文以六边形格网系统为例，六边形格网的情况比较复杂，目前常见的六边形剖分有3孔剖分和4孔剖分，其中3孔剖分常见的有1种，4孔剖分常见有4种，具体剖分情况如图3所示。

图3 孔径为3和4的六边形剖分格网Fig.3 Hexagonal grid with 3 and 4 apertures

对当前常见六边形剖分的各类单一三角面格网系统的形态进行分析，发现其本质是：剖分网格的几何形态是一致的，唯一不同的是三角面内所包含的剖分网格的数目和这些网格相对于三角面的位置姿态有所不同，并且不同层级的网格，不同剖分类型的网格，具体情况都是类似的。因此，单一三角面格网系统的形态在基扩充格网斜坐标系中主要由起始点位置S=（i0,j0）、起始边长度LN、起始边与I轴正向夹角θ3个因素共同决定。

其中起始点位置仅有两种情况：起始于格网中心节点和格网边界节点。起始边长度LN的不同取值可以生成不同尺度的格网系统。起始边与X轴正向夹角θ有3种情况：0°、30°、60°。六边形格网在基扩充格网斜坐标系中的具体剖分情况如图4所示。

图4 六边形格网的各种剖分类型Fig.4 Various partition types of Hexagonal grid

通过对图3和图4剖分类型的对比分析，以三角形边界处的格网形态为（具体情况已在图3灰色区域中标注）依据，可以得到这样的结论：图4a中的剖分方法对应常见的4孔1类剖分，4孔2类剖分，3孔偶数层剖分；图4b中的剖分方法对应常见的4孔3类剖分，3孔奇数层剖分；图4c中的剖分方法对应常见的4孔4类剖分；图4d中的剖分方法由于三角形边界处的格网形态与格网单元间不具有明显的几何对称关系，导致相邻三角面间跨面单元的归属问题[11]变的异常复杂，因此文章对此种剖分方法不做过多研究。

2.3 确定多尺度格网有效控制边界

为便于说明问题，本文以图4a中的剖分方法为例，通过对基扩充格网斜坐标系下该类剖分的单一三角面格网系统进行分析，发现多尺度格网生成的关键因素取决于起始边长LN。换言之，取不同的起始边长LN可以生产该剖分方法下不同层级的格网系统，因此研究起始边长LN的不同取值与格网层级间的数学关系是十分必要的。

正二十面体上的20个顶点间具有对称性、平等性、自反性。因此在固定剖分方法下，任一三角形边界顶点处的格网形态必然是一致的，本文以此为依据研究发现起算边长LN的取值可以表示为：LN=M×，公式中M∈Z+，R是六边形格网单元的半径。

为了体现格网系统不同层次间的关系，通过对该类剖分层级格网的研究，起算边长LN的具体取值与格网层级的关系表示为：当LN=3×2N-1×时，取适当的N值，对应该类剖分的第N层格网边界长度。

因此，定义边长起始节点所在格网单元为S=(i0,j0)，改变起算变换参数F=(S,LN,θ)=[(i0,j0)，3×2N-1×，60°]，得到该剖分方法下第N层剖分三角形的3个顶点格网坐标分别为AN(i0,j0)，BN(i0+3×2N-1,j0)，CN(i0+3×2N-1),j0+3×2N-1)。然后通过公式(1)可以快速确定多尺度格网边界顶点在空间直角坐标系下的坐标。具体结果如公式（2）：

2.4 计算多尺度格网单元格网坐标

通过上述研究发现，在给定任意起始节点所在格网单元格网坐标S=(i0,j0)和格网剖分层次N，在基扩充格网斜坐标系中，可以快速确定该层次剖分中任一格网单元的格网坐标(i,j)，其中格网坐标具体取值如公式（3）。第N层多尺度格网边界在基扩充格网斜坐标系中具体情况如下图5所示。因此，该层次剖分中的任一格网单元的7个节点的空间直角坐标皆可通过公式（1）求得。

图5 基扩充格网斜坐标系中第N层多尺度格网边界Fig.5 Multi-scale grid boundary in the N-th Layer in extended grid skew coordinate system

2.5 计算单一三角面格网系统单元的节点坐标

定义在第N层多尺度格网中，满足公式(3)的格网坐标为(i,j)的任意格网单元其7个节点中，设中心节点在空间直角坐标系中X轴Y轴Z轴上的分量分别表示为。从右上角按逆时针方向定义的6个边界节点在空间直角坐标系中X轴Y轴Z轴上的分量依次表示为：

通过多尺度格网和初始单一三角形的边界顶点在空间直角坐标中的几何关系，给出如下的坐标变换关系：

在公式（4）中，K表示缩放系数，H表示旋转矩阵，[x0y0z0]T表示平移向量，[x y z]T和[X Y Z]T分别表示变换前后的坐标值。具体有以下两类情况：

1）当Z=0时，仅考虑单一三角面格网系统单元（平面），其中，H是二维旋转矩阵；

2）当Z≠0时，考虑直接从初始单一三角面（平面）到空间三角面的变换，通过该方法直接生成二十面体表面的格网。

本文针对上述剖分方法，研究表明，给定起算变换参数F=(S,LN,θ)=((i0,j0),3×2N-1×60°))，通过相似缩放变换，其中，得出在基扩充格网斜坐标系下，与第N层多尺度格网中满足公式(3)的格网坐标为(i,j)的格网单元相应的单一三角面格网系统中单元的7个节点的空间直角坐标(定义为A'(i,j)，B'(i,j)，C'(i,j)，D'(i,j)，E'(i,j)，F'(i,j)，O'(i,j))分别表示如下：

2.6 建立不同层级单元节点的关联关系

建立不同层级格网单元节点间的关联关系，有助于空间数据的快速检索和应用分析的高效计算[12]。这里的关联关系指的是不同层级格网之间节点之间的对应关系，目标是为了建立层级节点之间的联系，即计算第N层格网坐标(I,J)的格网单元中，其7个节点在第N+1层格网中的具体格网坐标与节点位置。

图2所示的基扩充格网斜坐标系中，定义单一三角面格网系统的第N层格网单元中，格网坐标是(I,J)单元的7个节点位置分别用A(N,I,J)，B(N,I,J)，C(N,I,J)，D(N,I,J)，E(N,I,J)，F(N,I,J)，O(N,I,J)表示。7个节点间的位置排序遵循前述原则。通过对相邻层次格网的子单元、父单元的几何位置关系分析，得出如下结论：

公式（6）中，A(N,I,J)～B(N+1,2I+1,2J)表示第N层中格网坐标是(I,J)单元的边界节点A与第N+1层中格网坐标是(2I+1，2J)单元的边界节点B对应，公式（6）中其他情况类同。

3 结果与分析

根据前文的算法，本节针对平面单一三角面格网系统，计算得到其中各类参数，即3.1～3.6中的关键参数，见表1。

表1 各类剖分格网起算参数、有效控制边界、单元节点空间直角坐标Tab.1 Initial parameters of various partition grid types, eあective control of the boundary, unit node space cartesian coordinates

说明：表格中公式（5）见2.5，表格中公式（7）、（8）、（9）、（10）、（11）如下。

在公式（9），（10），（11）中，其他6个边界节点的空间直角坐标和中心节点O的情况类似，只需要用字母A，B，C，D，E，F分别替换上述公式中格网坐标是(i，j)单元的节点字母O即可。

另外，本论文还给出通过本方法直接转换到二十面体表面格网的情况，只需要令式（4）中的Z≠0，然后结合正二十面体的几何属性信息，建立合适的数学变换关系，可以完成从初始单一三角面（平面）到空间三角面的变换，即而从基扩充平面格网系统直接生成正二十面体表面的格网。

本文以4孔1类剖分为例，进行了下面的实验，具体实验结果如图6所示。

图6 正二十面体上4孔1类剖分第4、5、6、7层格网系统Fig.6 Divide the 4th, 5th, 6th, 7th grid system based on class 1 with 4-apertures on the icosahedron

4 结束语

本文提出了一种基于单元复制的通用化格网系统生成新算法，并以六边形离散格网为例开展了验证，实验证明该方法具有很好的通用性，适合多种类型格网生成过程的复用。该算法的核心是“简单单元复制+有效区域控制”，可以通过改变起算变换参数，生成任意形态任意类型的多尺度格网，避免了传统定制算法的局限性，实现了统一化生成各类球面格网系统的目标。相对于递归剖分方法而言，本算法在高层次格网生成中，进一步降低了算法的复杂度，并且为不同剖分类型全球离散格网系统间互操作性[13-14]问题研究提供了一种解决思路。