为了有效管理全球多尺度海量空间数据，多位学者研究了基于正多面体剖分地球的方法，将地球表面剖分成面积、形状近似相等，具有多分辨率层次结构的格网单元^[1]，称为正多面体全球离散格网系统 (discrete global grid system，DGGS)。常见的正多面体全球离散格网系统包括基于正六面体的四边形格网^[2-5]、基于正八面体的三角形格网^[6-10]、基于正八面体的菱形格网^[11-12]与基于正八面体的六边形格网^[13]和基于正二十面体的三角形格网^[14-17]、基于正二十面体的菱形格网^[18-19]与基于正二十面体的六边形格网^[20-22]。正多面体全球离散格网基于正多面体与球面的映射关系进行球面递归剖分^[23]，在全球范围内是无缝的、稳定的和近似均匀的，在空间数据表达^[24-26]、地球系统模式计算^[27-28]方面得到了深入研究，已成为当前全球离散格网研究热点之一。

总的来说，基于几种正多面体剖分球面各有千秋。正六面体的每个面均是正方形，可以生成类似于经纬格网的交错格网，因而被应用于地球系统模式计算中^[4]。正八面体的6个顶点可以定位在两个极点和4个等分的赤道点上，任意的经纬度坐标能较容易地定位到其中一个面上^[29]，鉴于该特性可以较容易地实现空间数据的表达，如经典的QTM模型^[6]。一般而言，多面体的面越小，在和球面转换过程中的变形就越小^[30]。相比较而言，二十面体是几种正多面体中面最小，最接近球面的正多面体，在此基础上剖分的全球离散格网变形相对较小，具有更均匀的几何性质，有利于全球空间数据集成与地球系统模式计算，因此逐步引起相关学者的关注^{[19, 23, 28]}。

相关研究者基于正二十面体三角形格网、六边形格网的剖分方法分别设计了其相应的格网编码方案。Fekete利用正二十面体三角形格网建立了SQT编码，并在此基础上研究了球面空间上的邻近分析和可视化^[14]。袁文采用L型空间填充曲线设计了正二十面体三角形格网的面片编码模型，并给出了面片节点的生成、访问及寻址算法^[17]。Sahr受到GBT结构的启发，在球面上借助与广义平衡三进制类似的编码运算，提出了正二十面体六边形格网编码和快速索引算法^[22]。童晓冲基于全球六边形格网提出了具备层次性的金字塔结构编码方案，采用递归方法设计并实现了编码与地理坐标的转换算法^[31]。正二十面体菱形格网虽然有成熟的剖分方案^[18-19]，但鲜有对其编码模型的研究。全球离散格网通常使用格网编码代替地理坐标进行各种空间操作，因此编码模型是格网系统至关重要的组成部分。格网编码既隐含了对应单元的位置信息，又表达了比例尺和精度^[30]。鉴于此，本文拟研究基于正二十面体菱形离散格网的编码模型。

基于正二十面体的球面菱形离散格网，有别于结构简单的正八面体，其初始菱形块的边线并不贴合经纬线，且几何性质、拓扑性质上都更为复杂，这对构建正二十面体球面菱形离散格网的层次编码模型和建立其与地理坐标间的映射关系转换带来了新的挑战。但值得注意的是，球面菱形离散格网采用四叉树剖分，在基础菱形内，其几何结构与平面栅格类似，适于现有的空间填充曲线的应用。目前应用最广泛的空间填充曲线包括Z曲线、Gray曲线及Hilbert曲线，不同的空间填充曲线具有不同的空间聚集能力，因而使得索引效率存在差异^[32]。其中，Hilbert曲线被证明能够最好保持空间点的局部邻接性^[33]，即保证在多维空间上邻近的对象映射到一维线性空间上也是邻近的，已被广泛应用于空间数据库的几何索引中^[34]。

因此，本文基于Hilbert曲线构建了正二十面体球面菱形离散格网的编码模型，在此基础上研究了格网单元编码与地理坐标的相互转化算法，并设计了相关试验，对转换精度与效率进行了分析。

1 正二十面体球面菱形离散格网剖分方法

本文所采用的剖分方法^[19]以球体的内接正二十面体作为球面格网划分的基础，在两个极点各放置一个顶点，其中一条通过北极的边线投影和0°经线重合 (图 1)。经纬度坐标下南北极两处存在极点奇异性，因而需将球面上的几何位置转为以球心为原点的三维直角坐标系 (图 2) 形式表达^[19]。合并相邻的南北向三角形形成10个球面菱形，在三维直角坐标系下对球面菱形进行细化，采用大圆弧平分法确定4条边的中点，点连线的球面投影即把球面菱形分为4个小球面菱形，依次类推，完成对整个球面近似均匀的剖分 (图 3)。基于此方法所得到的格网整体分布均匀，格网单元面积、角度形变小、排列紧致^[19]。

2 正二十面体球面菱形离散格网编码模型

本文所提出的编码模型包括区位码和Hilbert码两部分组成，可以采用如下CODE的形式表示

式中，D表示区位码，用于索引正二十面体球面菱形格网的十个初始菱形块，采用二进制码表示，在计算机上占2个字节 (16位)。H为Hilbert码，以每个初始菱形块左顶点对应的格网单元为起始点建立Hilbert填充曲线，对某个基菱形上的格网单元进行索引，同样以二进制表示，每层于末尾增加2位 (bit)，因而隐含了格网的剖分层级、父子单元的Hilbert码。本文规定以正二十面体0号基菱形为例，由Hilbert码便能定位格网单元在基菱形内的位置，如图 5中第2层的菱形块对应的格网编码为00001101。

3 格网编码与地理坐标的相互转换

陆锋根据Hilbert曲线的递归构造过程，提出了一种基于空间层次分解的Hilbert空间排列码生成算法，使用可迭代规则对空间进行逐步细化^[35]。其中，Hilbert曲线分为4种子象限形态，这决定了存在4种形式的状态转移向量，即{1，2，3，4}，{1，4，3，2}，{3，2，1，4}和{3，4，1，2}(图 6)。球面菱形离散格网在初始菱形块内可以看作是倾斜的平面栅格，通过旋转平面栅格上的状态可以得到球面菱形块上Hilbert曲线的4种状态转移向量 (图 7)。在空间层次分解时，结合象限的Hilbert码和状态转移向量可以推导出子象限的Hilbert码和状态转移向量，其中状态转移向量表明子象限中曲线的旋转与反射。

本文借鉴Hilbert曲线的空间层次分解方法，设计了地理坐标与格网编码相互转换的算法。

地理坐标到格网编码的转换实质是任意球面点在格网上确定其所属格网单元的过程。本文首先根据待求点所在的地图比例尺查找确定其所在的目标层次L，在判断待求点是否在指定剖分层次的基础上，通过平面方程判别法逐层求取待定点在下一剖分层次中所属的子菱形Hilbert码，最终得到地理坐标所属的球面菱形格网单元的编码。具体的步骤如下 (图 8)。

(1) 采用文献^[19]中的经纬度与空间之间坐标的转换方法，将待定点P的地理坐标转换至空间直角坐标系 (图 2)，初始化状态转移向量V，Hilbert码H，剖分计数l：H=0，l=0；

(2) 判断地理坐标的落点区域，即正二十面体菱形格网的初始菱形块d，得到对应的二进制区位码：D=d，并以该初始菱形块为当前基菱形；

(3) 若当前菱形的剖分层次l为目标层次L，则当前菱形块即为地理坐标在指定剖分层次上对应的格网单元，将二进制区位码D与Hilbert码H组合成编码CODE并输出，算法结束；否则，进入下一步；

(4) 以当前菱形边界V₁V₂的中点M₁、V₀V₃的中点M₃以及球心O构成平面OM₁M₃，求其三维空间中的平面方程。由于平面OM₁M₃经过球心O，则可得平面OM₁M₃的法向量。根据法线方程可以直接给出平面OM₁M₃的点法式方程x_N×x+y_N×y+z_N×z=0。此时问题便转化为判断点和平面的位置关系，据此即可判断待定点P和平面OM₁M₃的位置关系。将空间直角坐标P代入平面方程，若大于0，则在平面OM₁M₃上方，反之在平面下方；

(5) 采用的同样的方法，再对另外一个方向上菱形边界V₀V₁的中点M₀、V₂V₃的中点M₂和球心O构成平面OM₀M₂进行位置判断 (图 9)，得出P点在菱形块内所处的象限，由象限号和状态转移向量V通过表 1获得该象限Hilbert码的二进制增量k和新状态转移量v，重新设定状态转移向量V=v，将增量k添加在Hilbert码尾部，更新Hilbert码和剖分计数：H=Hk，l=l+1，返回步骤 (3)。

从格网编码到地理坐标的转换实际上是一个格网单元的编码到格网单元中心点的球面经纬坐标的转换。本文中格网单元的中心点即为菱形单元短轴的中点。转换的基本思想是根据编码进行单元的逐层递归逼近，直至定位到编码所对应的菱形单元，再求出中点即可 (图 10)。具体过程如下：

(1) 根据正二十面体与球面的对应关系获取格网编码中区位码D对应的初始菱形顶点坐标V₀、V₁、V₂、V₃；

(2) 设定初始状态转移向量V，初始化格网剖分计数：l=0，目标层次L为Hilbert码M位数的二分之一；

(3) 若当前菱形的层次l为目标层次L，输出菱形单元中点的经纬度坐标，算法结束；否则，进入下一步；

(4) 取Hilbert码中的第2l+1、2l+2位二进制码组成Hilbert增量k，由增量k和状态转移向量V根据表 2确定单元在当前菱形块下所属的象限号，依据状态转移向量V与得到的象限号通过表 1获得该象限内的新状态转移向量v，重新设定状态转移向量：V=v；

(5) 获取当前象限对应菱形块的四个顶点，更新顶点坐标V₀、V₁、V₂、V₃和剖分计数l：l=l+1。返回步骤 (3)。

4 试验与分析

有多种评价标准可以用来衡量地理坐标与格网编码间相互转换的精度，最常见的是绝对精度^[31]。绝对精度指的是将地理坐标形式的点或球面编码转换为对应的格网编码或地理坐标后，再将其转换回输入形式得到输出结果，对比输入和输出的坐标点球面距离或者编码值，如果二者在一个格网剖分单元内，则转换算法符合格网绝对精度的要求。

为了分析本文提出的编码转换算法效率与精度，选择不同数量的 (2000、1万、10万、100万) 随机坐标点在不同层次下 (12、14、16、19、21层，对应比例尺分别为1:400万、1:100万、1:25万、1:5万、1:1万) 进行编码与地理坐标的相互转换试验，与White提出的Morton码编码方案^[18]进行比较，并以绝对精度^[31]为编码转换精度评判标准。主要硬件配置环境为：CPU为Intel (R) Core (TM) i7-4770 3.10 GHz，8 GB内存，1 TB硬盘，采用VC++ 2013作为基础开发平台。试验结果如下 (表 3)。

从图 11可以发现，地理坐标转格网编码后的绝对精度误差随着剖分层级的增加而迅速减小，综合对应剖分层级 (12、14、16、19、21层) 下菱形单元的边长表 4来看，转换的绝对精度始终控制在一个格网单元内。从表 3可知，格网编码转地理坐标后格网编码值保持不变，前后所表示的格网单元始终为同一个单元。综上分析，本文提出的坐标编码转换算法以递归逼近方法为基础，相当于进行了一次格网单元局部的递归剖分，因而可以保证转换算法满足绝对精度的要求。

将本文所提出的编码方案与White在Z曲线基础上提出的编码方案^[18]进行地理坐标与格网编码转换的对比试验。从图 12可以发现，随着剖分层级的增加，基于两种编码模型的地理坐标与网格编码的转换效率均逐步下降，转换耗时呈现线性增长，而两者之间的效率表现几近一致。由于格网的剖分中心、菱形块的棱边中点与球心构成两个扇面，子菱形归属只需要与这两个扇面进行2次位置判断 (图 9)，因而均能够在顾及转换精度的前提下实现坐标向编码的快速转换。

值得注意的是，相比White在Z曲线基础上提出的编码方案，本文所提出的基于Hilbert曲线的编码模型对于空间数据的局部邻接性保持更好，因此在后续的空间数据索引、拓扑分析等方面具有明显的优势。

5 结论

本文在正二十面体球面菱形离散格网剖分方法的基础上，考虑到在球面上难以直接构建全球连续的菱形格网编码，利用基础菱形内部与二维栅格的相似性，采用分区编码与二进位的Hilbert曲线结合的方案，构建了一种具有层次结构的正二十面体球面菱形离散格网编码模型。基于此，本文借鉴平面栅格上基于空间层次分解的Hilbert码的生成算法，采用迭代细化的规则方法构造球面Hilbert曲线的同时实现了经纬度向格网编码的转换，并进一步实现了对球面菱形离散格网编码的反向解码。研究表明，由于球面菱形离散格网在结构上与平面规则格网具有相似性，因此基于Hilbert曲线的编码模型能够在地理坐标与格网编码转换方面兼具效率与精度，可以支撑全球海量空间数据建模、集成管理以及各类空间分析。