格网系统是一种多分辨率栅格层次数据结构，本质是通过规则格网单元的分解与聚合将连续平面映射到多层次格网中，且每个单元具有包含层次关系的唯一编码^[1]。格网系统的离散性使其更适合计算机处理，有助于多源、异构、多尺度空间数据的融合及多要素的空间叠加和复合分析^[2]。三角形格网系统广泛应用于地形建模^[3-5]、空间数据索引与可视化^[6]、地图综合模型^[7]；矩形(四边形)格网系统广泛应用于资源环境综合科学调查^[8]、遥感影像数据存储与管理^[9]及地图多级显示^[10]。

相比于上述格网，六边形格网具有形态紧凑、平面覆盖率和角度分辨率高及一致相邻等特性^[11]，更有利于空间数据的组织、处理和分析。然而，六边形不具备自相似性，格网层次关系描述及计算是研究难点之一^[12]。文献[13-14]提出“广义平衡三进制”(generalized balanced ternary，GBT)，将高分辨率六边形单元聚合为低分辨率单元，但其相邻层次单元的方向连续旋转，导致相关应用较复杂^[15]。文献[1, 16-17]分别设计了三孔六边形格网系统编码方案，并借助改进的GBT实现不同层次单元快速索引，其编码运算无法完全用二进制优化，效率不甚理想^[18]。文献[12]以“格元-格点-格边-格心”概念模型^[19]为基础，建立了三孔六边形格网系统编码方案，并揭示其数学本质。文献[18]针对四孔六边形格网系统设计了“六边形四元平衡结构”(hexagonal quaternary balanced structure，HQBS)，但在编码运算过程中频繁出现“正则化”失败回退重算的情况^[12]。针对这一缺陷，文献[20]设计了“六边形格点四叉树”(hexagon lattice quad tree, HLQT)，其编码运算效率优于HQBS及PYXIS，但其码元分布不对称，不利于邻近单元查询。文献[21]采用菱形块四叉树组织球面六边形格网层次结构，并提出格网编码方案。在实践应用方面，ESRI公司在ArcGIS GeoAnalytics Server中采用六边形格网系统分析、展示数据(http://www.esrichina.com.cn/openclass-2017.html)。UBER公司在其流处理系统中采用六边形格网系统实时采集、分析车辆信息(http://www.infoq.com/cn/presentations/uber-stream-processing-system-and-practice)。欧空局的SMOS卫星也采用六边形格网存储、处理影像数据^[22]。

综上所述，当前六边形格网系统的研究已取得可喜进展，但现有成果仍存在编码方案不完善、运算效率待提升的问题。本文引入复进制数理论, 依据间隔层次单元隶属关系，建立平面四孔六边形格网系统数学模型，在此基础上提出具有结构对称性的等效单元编码方案、定义编码运算并归纳运算规则、设计编码索引及编码与笛卡儿坐标互换算法，尝试提高编码操作效率。

1 复进制数建模

复进制定位计数系统(complex radix positional system)由计算机科学泰斗Donald E.Knuth于1960年提出^[23]，它的基数(进制)、数字位集合元素可以是复数，本质上是二进制、十进制等常见定位计数系统在复数域上的推广。复数在复进制定位计数系统中的表示形式称为复进制数(complex radix number)。在全球离散格网的相关研究中，通常采用笛卡儿坐标描述单元位置^[24]，该方式虽然直观，但与编码关联性不强。研究表明，以复进制数形式表示单元位置，更有助于编码设计及其运算规律推导^[25-26]。

作为位置、对象及数据与格网单元的关联点，通常利用格心代替格网单元^[19]。四孔六边形格网的剖分规律如图 1所示，第n层格心由第n-1层格心及单元各边中心组成，相邻层次单元边长为2倍关系，且单元方向不随层次改变。

四孔六边形格网层次关系，还可从间隔层次递推的角度分析。如图 2所示，第3层单元可由第1层单元剖分形成的12个子单元及第2层单元面积缩小得到的继承子单元构成。按此递推关系，以1、2层为基础，第n(n≥3)层单元可由第n-1层单元继承及第n-2层单元剖分得到。该层次关系避免了逐层剖分过程中同一子单元隶属多个父单元的情况，无论是继承子单元还是剖分子单元都具有唯一父单元，有利于建立复进制数模型。

令首层格网单元的格心位于复平面原点处，复数对应于与实轴夹角为3的单位向量。继承子单元面积缩小但格心位置保持不变，因此不需单独表示；剖分子单元是通过剖分产生的新单元，其格心位置可借助复进制数集合D={ω, ω², ω³, ω⁴, ω⁵, ω⁶, ω+ω², ω²+ω³, ω³+ω⁴, ω⁴+ω⁵, ω⁵+ω⁶, ω⁶+ω}描述，如图 3所示，则四孔六边形格网系统的格心集合为

式中，“∪”、“+”是集合间的并运算、加法运算符。

证明：设第n-1层格网单元外接圆半径为r_n-1，第n层格心集合为L_n，根据格网系统生成规律，L_n可表示为

因为，则

所以

式中，L₁=0，令L₂=L₁∪(L₁+D)，则此时有

式(1)表明，格心集合本质上是特定形式的定位计数系统。令l_k∈L_k，d_k∈D(k=1, 2, 3, …, n-2)，对于α∈L_n(n≥3)，当α属于继承子单元时，α=l_n-1，当α属于剖分子单元时，α=l_n-2+d_n-2。继承子单元的格心继承于邻层单元的格心，二者位置一一对应；剖分子单元由隔层单元在内部剖分得到，任意剖分子单元均唯一归属于相应隔层单元。因此，该层次关系中对于继承子单元和剖分子单元格心位置的描述均具有唯一性，便于后续编码操作。

不同于HLQT、HQBS等四孔六边形格网系统逐层编码方案以及PYXIS三孔六边形格网系统奇、偶层分别编码方案^[1]，这里的定位计数系统建立在层次递推关系的基础上。随层次增加，格心呈现向内凹陷加密的趋势。如图 4所示，用不同大小和颜色的点表示L₁-L₅中的格点，图中空缺部分可由首层相邻单元生成的子格点填补。

2 等效编码方案

第n(n≥3)层中包含继承子单元和剖分子单元，对它们分别编码。复进制数集合D对应的等效码元决定了剖分子单元编码，考虑到编码的对称性及运算的简易性，制定如下编码规则：将集合D划分为两个子集，子集{ω, ω², ω³, ω⁴, ω⁵, ω⁶}对应6个单位向量，以码元0k替换ω^k(1≤k≤6)，对应码元集合为{01, 02, 03, 04, 05, 06}。子集{ω+ω², ω²+ω³, ω³+ω⁴, ω⁴+ω⁵, ω⁵+ω⁶, ω⁶+ω}对应单位向量相加得到的向量，以码元k0替换ω^k+ω^k+1(1≤k≤6), 对应码元集合为{10, 20, 30, 40, 50, 60}，第n层原点单元编码记为个。

由上述规则直接给出L₁与L₂的编码，如图 5所示。以L₁、L₂和集合D的编码为基础，依据层次递推关系建立L_n(n≥3)编码，由L_n-1得到L_n中继承子单元的过程编码层次增加但格心位置不变，因此规定L_n(n≥3)中继承子单元编码由L_n-1编码末尾添加“00”得到，如图 5 L₃编码中黑色单元所示。剖分子单元的格心由第n-2层格网单元依据集合D剖分得到，但第n层与第n-2层并非相邻层次，因此在L_n-2编码末尾添加“00”表示跨层后再添加集合D对应的等效码元，如图 5 L₃编码中红色单元所示。

根据上述编码规则，每层码元由2位字符构成，由于单元归属明确，单元编码具有唯一性；复进制数集合D中的元素关于原点对称，编码具有对称性，有利于编码索引操作。

3 编码运算规则和索引算法

编码记录了格网单元相对原点的距离和方位，可通过一维编码运算实现二维复平面内的几何操作，达到降维目的。相较于浮点数，编码更适合计算机存储与处理，有利于提高运算效率。

3.1 编码加法运算规则

与矢量加法类似，可利用平行四边形法则定义编码加法运算^[27]，以符号⊕表示。按照层次关系，L_n(n≥3)与L_n-1和L_n-2均相关，在编码加法中，以4位(即每两层)编码为基础码元

依照平行四边形法则形成25×25加法查找表，见表 1。将加法表结果统一记为8位，则加法运算的进位操作仅存在于相邻基础码元。

以加法表为基础讨论任意层次编码加法，第n(n≥3)层编码Code₁和Code₂由2n位字符构成。由右至左按4位字符分组，不足4位的部分以“00”补足，编码可分为(⌈⌉表示向上取整)组基础码元。每组基础码元依加法表得到的8位结果仍按4位字符由右至左分为保留位编码B及进位编码J。按图 6所示错位相加法进行计算，依次将每一步中的B₁^x记入结果，直到全部计算完毕。对于J₁^m，当其不为“0000”时计入结果。需要注意的是，计算每一步中最后的进位编码(如J_m-1²、J_m-2³等)时应进行累加。以J_m-1²为例，需由B_m¹⊕J_m-1¹得到的进位编码再与J_m¹相加得到。

以第5层编码加法0000020040⊕0001004000为例，由“00”补足后化为{0040, 4000}，{0002, 0100}，{0000, 0000}3组基础码元。0040⊕4000=00400010，对应保留位编码B₁¹=0010，进位编码J₁¹=0040；同理0002⊕0100=00010030，B₂¹=0030，J₂¹=0001，0000⊕0000=00000000，B₃¹=0000，J₃¹=0000。B₁¹=0010计入结果，根据错位相加法，J₁¹⊕B₂¹=0040⊕0030=00040001，则B₁²=0001，J₁²=0004，J₂¹⊕B₃¹=0001⊕0000=00000001，则B₂²=0001，J₂²=0000⊕J₃¹=0000，B₁²=0001计入结果。最后，计算J₁²⊕B₂²=00004⊕0001=00000000，则B₁³=0000计入结果，J₁³=0000⊕J₂²=0000，不计入结果，可得编码000000010010，除去补足的“00”，最终结果编码为0000010010，与图 7所示平行四边形法则的计算结果一致。

对于第5层编码加法0000010003⊕0000000003，错位相加法得到的结果为0000010300，而实际结果为0000001000。分析可知，错位相加法以单元0000010003为中心加上0000000003得到最终结果，而实际编码0000001000由第4层继承得到，其生成中心为原点单元，中心单元的不一致导致计算结果的不一致。出现该问题的编码存在连续非“00”编码的明显特征，以编码α=α₁α₂…α_2n-1α_2n为例，α_k-3α_k-2α_k-1α_k是出现连续非“00”编码的位置，将α₁α₂…α_k转换为α₁α₂…α_k-3α_k-200⊕00…00α_k-1α_k的结果，并在末尾加上α_k后的编码α_k+1…α_2n可将编码α转化为正确编码，避免“回退”操作。

3.2 编码乘法运算规则

编码乘法对应旋转和缩放操作^[26]，以符号⊗表示。集合{01, 02, 03, 04, 05, 06}中的码元分别对应旋转0°、逆时针旋转60°、逆时针旋转120°、旋转180°、顺时针旋转120°、顺时针旋转60°。据此可得编码乘法运算查找表，见表 2。

乘法运算的码元及结果均为2位编码，将编码按2位字符分组，逐组运算并记录结果即可。以0000060002⊗02为例，依乘法表，00⊗02＝00，02⊗02＝03，02⊗06＝01，则结果编码为0000010003，与逆时针旋转60°结果一致。

3.3 编码索引算法

编码索引包含的层次关系及邻近关系查询是空间聚类、邻近分析等操作的基础。将第n层编码α=α₁α₂…α_2n-1α_2n与相应单位码元e_i=个(i=1, 2, 3, 4, 5, 6)相加可得邻近单元编码α⊕e_i(i=1, 2, 3, 4, 5, 6)。

层次单元查找分为子单元和父单元查找，按照第3部分继承子单元和剖分子单元的编码规律即可实现子单元的查找。父单元查找也分为邻层和隔层两种情况，对于第n(n≥3)层编码α=α₁α₂…α_2n-1α_2n，若α_2n-1α_2n=00，则该单元继承于第n-1层，将末尾编码“00”去掉可得到邻层父单元。若α_2n-1α_2n≠00，则该单元由第n-2层的单元剖分得到，将末尾4位编码去掉可得到隔层父单元。

4 平面笛卡儿坐标与编码转换 4.1 格网编码到平面坐标的转换

编码本质上就是复进制数的等效表达方式，因此每层中的2位编码与集合D中的复进制数一一对应，将第n(n≥1)层编码α按2位字符由左至右划分为1到n层，逐层将编码对应的复进制数及相应缩放系数相乘并累加即可得到平面坐标。设首层格网单元边长为r₁，则第n层缩放系数为。

以第5层编码0000600002为例，由左至右第3层码元60对应复进制数为ω⁶+ω，缩放系数为，第5层码元02对应复进制数为ω²，缩放系数为，则编码对应的平面坐标为(ω⁶+ω)+ω²。

4.2 平面坐标到格网编码的转换

对于复平面内任意坐标(ℝ_P, Ⅱ_P)，将其转化为斜轴坐标(i, j)，根据格网层次n及单元距离d_n可得到i，j分量对应的格网单元，如图 8坐标轴i，j上红色单元所示，同时夹角为120°的斜轴坐标系本身满足平行四边形法则，利用编码加法将i，j上的格网单元相加便可得到平面坐标对应的编码。

5 试验与分析

本文编码方案可借助十六进制数实现。码元{00, 01, 02, 03, 04, 05, 06, 10, 20, 30, 40, 50, 60}与十六进制数{0, 1, 2, 3, 4, 5, 6, A, B, C, D, E, F}对应, 可将字符编码转化为更适合计算机存储的十六进制整数，并借助二进制位运算提高计算效率。相比于其他现有六边形格网编码方案，HQBS和HLQT在编码结构和编码运算效率方面具有较大优势^{[18, 20]}。因此，选择上述两编码方案进行对比试验。本文采用Visual Studio 2012开发工具实现本文编码方案对应的编码加法、邻近单元查询、坐标与编码互换算法，并与HQBS及HLQT相应算法比较。全部代码均编译为Release版本，并在相同台式机(硬件配置：Intel Core i5-4590M CPU @ 3.30 GHz，8 GB RAM；操作系统：Windows 7 x64旗舰版)上测试。4-11层中逐层随机生成2500个单元，不重复相加，得到3 123 750组加法运算实例；随机生成1 000 000组坐标作为笛卡儿坐标到编码转换试验数据；第4层随机选取256个单元，5-11层以4倍递增确定随机选取单元个数作为编码到坐标转换及邻近单元查询试验数据。以单位时间内完成操作的次数为效率指标，统计10次测试的平均值，部分试验结果如图 9-图 12所示。

分析试验结果，得出以下结论：

(1) 本文方案编码加法的平均效率约为HQBS的10.11倍，HLQT的2.00倍。HQBS格心与顶点混合编码，运算规律复杂且易出现编码“正则化”失败回退重算的情况。HLQT逐层进行码元相加，计算非首位码元相加时，除了得到当前位计算结果码元，还会在左侧所有位产生进位码元。虽然相比于HQBS和HLQT，本文方案的加法表规模较大。但任何编码方案中的加法表均以二维数组形式存储，通过码元在二维数组中的位置，可直接定位相应结果，因此加法表的规模几乎不会影响其运算效率。同时，本文方案以每两层为基本计算单元，且进位只发生在相邻基本单元之间，相比之下计算量更小，因而效率更高。

(2) 本文方案笛卡儿坐标转换为编码的平均效率约为HQBS的4.33倍、HLQT的0.80倍。本文方案与HLQT在i、j轴上的编码分布均具有二进制数特征，可采用位运算加速，因此两者效率均高于HQBS。相较于HLQT，由于本文方案的编码具有对称性，在i、j轴上的部分编码与二进制数的转换略复杂，导致转换效率略低。笔者认为，与对称编码在测试(1)、(3)、(4)中带来的诸多优势相比，该转换算法约20%的效率损失可以接受。

(3) 本文方案编码转换为笛卡儿坐标的平均效率约为HQBS的6.29倍、HLQT的2.36倍。HQBS算法涉及较多的浮点数运算，本文方案与HLQT均采用整数代替浮点数运算，因此两者效率均高于HQBS。相较于HLQT，由于本文方案的编码以“00”和非“00”码元交替形式存在，计算中只需考虑非“00”码元，因而效率更高。

(4) 本文方案编码邻近单元查询效率约为HQBS的4.23倍、HLQT的1.83倍。HQBS和HLQT都是基于非对称四叉树设计的编码方案，邻近单元查找涉及的编码运算较复杂。同样的操作，本文编码方案只需根据中心对称的单位码元e_i=(i=1, 2, 3, 4, 5, 6)即可得到，效率更高。

6 总结与展望

本文结合平面四孔六边形格网结构特点，引入复进制数理论，通过间隔层次格网单元隶属关系建立数学模型，以此为基础提出等效编码方案、定义编码运算并归纳运算规则、设计编码索引及编码与坐标相互转换算法。相较于同类成果，本文编码方案具有结构对称性，在编码加法、邻近单元查询及本方案编码转换为笛卡儿坐标等方面效率优势明显，有助于各类格网处理、分析算法的实现，具有较大的应用潜力。

本文仅在平面内建立了编码方案并实现了相关操作，借助多面体及合适的映射关系可将本文编码方案扩展到球(参考椭球)面，实现全球四孔六边形离散格网系统的构建。在此过程中涉及的相关问题，将是笔者下一步研究的重点内容。