平面四孔六边形格网系统复进制数建模及编码运算

杜灵瑀，马秋禾，贲 进，王 蕊

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

格网系统是一种多分辨率栅格层次数据结构，本质是通过规则格网单元的分解与聚合将连续平面映射到多层次格网中，且每个单元具有包含层次关系的唯一编码[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倍关系，且单元方向不随层次改变。

图1 四孔六边形格网系统示意图Fig.1 Diagram of aperture 4 hexagon grid system

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

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

所以

此时有

(1)

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

2 等效编码方案

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

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

图2 基于层次递推的格网系统层次关系Fig.2 Hierarchical relation of grid system based on hierarchical recursion

图4 L1—L5在复平面上的分布Fig.4 Distribution of L1—L5 in the complex plane

图3 集合D在复平面上的表示Fig.3 Representation of set D in the complex plane

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

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

3.1 编码加法运算规则

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

图5 编码示意图Fig.5 Diagram of code

{0000,0001,0002,0003,0004,0005,0006}

{0010,0020,0030,0040,0050,0060}

{0100,0200,0300,0400,0500,0600}

{1000,2000,3000,4000,5000,6000}

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

表1 编码加法查找略表

图6 错位相加法示意图Fig.6 Diagram of staggered addition

对于第5层编码加法0000010003⊕0000000003，错位相加法得到的结果为0000010300，而实际结果为0000001000。分析可知，错位相加法以单元0000010003为中心加上0000000003得到最终结果，而实际编码0000001000由第4层继承得到，其生成中心为原点单元，中心单元的不一致导致计算结果的不一致。出现该问题的编码存在连续非“00”编码的明显特征，以编码α=α1α2…α2n-1α2n为例，αk-3αk-2αk-1αk是出现连续非“00”编码的位置，将α1α2…αk转换为α1α2…α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位编码，将编码按2位字符分组，逐组运算并记录结果即可。以0000060002⊗02为例，依乘法表，00⊗02=00，02⊗02=03，02⊗06=01，则结果编码为0000010003，与逆时针旋转60°结果一致。

3.3 编码索引算法

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

4 平面笛卡儿坐标与编码转换

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

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

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所示。

图7 第5层编码加法示意图Fig.7 Diagram of code addition in the 5th level

图8 笛卡儿坐标到编码转换示意图Fig.8 Diagram of Cartesian coordinates to code

图9 编码加法效率对比Fig.9 Efficiency comparison of code addition

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

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

图10 笛卡儿坐标到编码转换效率对比Fig.10 Efficiency comparison of Cartesian coordinates to code

图11 编码到笛卡儿坐标转换效率对比Fig.11 Efficiency comparison of code to Cartesian coordinates

图12 邻近单元查询效率对比Fig.12 Efficiency comparison of search for neighboring cells

(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”码元，因而效率更高。

6 总结与展望

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

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