《测绘学报》

构建与学术的桥梁拉近与权威的距离

平面四孔六边形格网系统编码运算

王蕊1,贲进1
,杜灵瑀1,周建彬1,李祝鑫2

1.信息工程大学地理空间信息学院,河南郑州450001;2.61287部队,四川成都610000

收稿日期:2017-06-30;修回日期:2018-03-23

基金项目:国家重点研发计划(2018YFB0505301);国家自然科学基金(41671410)

第一作者简介:王蕊(1995-),女,硕士生,研究方向为空间数据模型。E-mail:wr_paper@126.com

通信作者:贲进,E-mail:benj@

EncodingandOperationforthePlanarAperture4HexagonGridSystem

WANGRui1,BENJin1,DULingyu1,ZHOUJianbin1,LIZhuxin2

Abstract:Discreteglobalgridsystemisanewdatamodelwhictem,thispaperdesignsanencodingschemenamedHexagonLatticeQuadTree(HLQT).Codeoperationsaredefined,rulesofthemaregeneralizedandbasedonthese,transformatio,HLQTovercomesthedisadvantagescaughtbyencodingschemeswhic,ut6timesthatofPYXISandabout5timesthatofHQBS,theefficiencyofthetransformalgorithmfrom2-dimensionalcoordinatestocodesisabout5timesthatofHQBS,andtheefficiencyofthetransformalgorithmfromcodesto2-dimensionalcoordinatesisabout3timesthatofHQBS.

Keywords:discreteglobalgridsystemhexagonaddressingcodecoordinatetransform

全球离散格网系统(discreteglobalgridsystem,DGGS)是一种以整个地球为研究对象、以地理空间位置数字化表达为主要特征的新型数据模型。它将地表剖分为均匀统一的层次格网单元结构,并给每个单元赋予全局唯一的地址编码代替传统地理坐标参与运算和操作[1]。与传统空间数据模型相比,其优势体现在:面向全球,能够提供全球一致的空间基准;天然多尺度,有助于处理多源、异构地理空间数据;纯编码运算,有助于提高数据处理效率。

与三角形、四边形相比,六边形具有拓扑关系一致、采样效率高等特点,有助于提高空间分析算法精度[2-3]和可视化表达效果[4-5],其主要缺陷是不具备自相似性,无法直接建立多分辨率格网单元层次关系。文献[6]提出一种多维数据结构,二维特例可通过“广义平衡三进制”(generalizedbalancedternary,GBT)运算实现不同层次六边形单元的编码索引[7]。文献[8]改进GBT的不足,提出HIP(hexagonalimageprocessing)结构,并将其应用于影像处理。文献[9-11]区分奇(偶)层建立了平面三孔六边形格网系统的单元编码索引并将其扩展到球面,在此基础上,文献[12]实现了傅里叶变换算法。但这类方案对奇(偶)层单元分别处理,增加了编码运算复杂度,索引效率不高。文献[3,13-15]采用“单元隶属图形”建立六边形格网的编码索引,主要不足是破坏了六边形单元的完整性。文献[16-17]借助格网单元中心与顶点建立了平面四孔六边形格网系统的“四元平衡结构”(hexagonalquadbalancedstructure,HQBS)并将其扩展到球面,主要缺陷是概念模型复杂,且会在编码运算中会出现“正则化”失败需要回退重算的情况,严重影响索引效率[18]。

综上所述,单元编码索引是六边形全球离散格网系统研究的核心问题之一,平面格网系统层次关系描述及编码方案设计则是该问题的突破口。现有成果仍然存在不足,亟须改进完善。本文根据平面四孔六边形格网结构特点,设计“六边形格点四叉树”(hexagonlatticequadtree,HLQT)层次编码结构,定义编码运算等效代替格网单元的几何操作,归纳编码运算规律并实现了对应的算法。借助编码运算,设计并实现二维直角坐标与单元编码的相互转换算法,通过对比试验验证本文方案的正确性与优越性。

1编码方案设计

为了便于讨论,本文将平面离散格网系统的单元中心定义为“格点”,作为单元的等效替代对象,四孔六边形剖分产生的格点以及剖分规律如图1所示:第n(n≥1)层格点是由第n-1层格点与n-1层单元各边中点组合而成,第n-1层单元边长是第n层单元边长的2倍,单元方向不随层次变化。


图1四孔六边形格网系统格点与剖分规律示意图nalgrid

格点位置通常采用笛卡尔坐标描述[10,14,17,19-20],这种方法虽然直观形象、便于理解,但因采用了二维独立变量(如x、y),表达形式不简洁且不易与编码建立联系。故本文采用复数表示格点位置[21-22],建立平面四孔六边形格点系统的唯一性描述方案[23]。

在复数平面上,令


,D′={0,ω′,ω′2,ω′3,ω′4,ω′5,ω′6},


,D={0,ω,ω2,ω3},经过严格证明,四孔六边形格网系统的格点集合可唯一表示为


(1)

Ln(n≥1)表示第n层格点集合。式(1)中“∑”并不是集合之间的并(∪),而是集合之间相加,即从参与计算的n个集合中任选一个元素,一共n个元素求和。图2所示为L1与L2层次关系示意图。L1为
D′表示的7个格点集合,以L1中每个格点为中心,按集合D中元素的方向各生成4个子格点构成L2,其中L2的中心子格点与L1重合。按上述思路,可以依次生成Ln(n>2)。


图2L1与L2格点层次关系示意

以上描述方案具有唯一性,生成的子单元互不重叠,每个子单元具有唯一的父单元。Ln是典型的四叉树层次结构,即随着剖分层次的深入,L1中的每个格点都分别生长为一颗独立的四叉树,这为格网编码提供了理论依据。

以上方案还表明:四孔六边形格点系统与实数域上的二进制、八进制、十进制等定位计数系统的本质完全相同,是复数域上一种特定形式的“数”[24]。该定位计数系统的进制为2,数字位集合为D′和D。对于x∈Ln,其具体形式为


(2)

或与十进制数类似,省略幂指数的底,式(2)可记为


(3)

式(3)中的小数点表示格点随剖分层次向内部生长,越来越密集。用不同大小和颜色的点表示L1~L5中的元素,它们在复平面上的分布如图3所示。


图3L1~L5在复平面上的分布~L5inthecomplexplane

为了便于计算机处理,需要对格网系统中的每个单元指定一个地址码,这一过程称为“编码”。编码方案必须满足完整性、唯一性与层次性[19]。

根据式(3),取di幂指数,称为码元,即d1=ω′e→e(1≤e≤6),di=ωe→e(1≤e≤3,2≤i≤n),则数字位集合D′替换为码元集合E′={0,1,2,3,4,5,6},数字位集合D替换为码元集合E={0,1,2,3}。得到x∈Ln的唯一编码标识e1e2…en(e1∈E′,ei∈E,2≤i≤n)。显然,该编码方案还满足完整性和层次性的要求。

L3的完整编码如图4所示,不同格点层次隶属和关联用不同大小的圆和三角形表示。将这一格点层次结构命名为“六边形格点四叉树”(hexagonlatticequadtree,HLQT)。


图4L3编码及运算示意图

2编码运算规律

HLQT编码记录了某一格点与原点(编码全部为0的格点)的距离和方位,是复数平面上“定位计数系统”的等效表示。因此格点对应复数之间的各类运算可从二维复数平面空间降维到一维编码空间中实现,且编码相较于浮点数更适合计算机处理,可有效提高计算效率。

2.1定义

在复数平面上,任意两复数的加、减法满足平行四边形法则。在HLQT编码集合
n(n≥1)中,设编码α,β∈n,定义算子C(…)为计算给定HLQT编码对应复数,算子H(…)为计算给定复数的HLQT编码,则α、β的编码加法定义为复数C(α)+C(β)对应的HLQT编码H[C(α)+C(β)],表示为γ=α⊕β=H[C(α)+C(β)]。如图4所示,α⊕β=123⊕010=233。根据归纳、验证,编码加法满足交换群的性质[25]。HLQT的减法运算满足向量相减的平行四边形法则,例如:α⊝β=123⊝010=033。⊕与⊝是逆运算,编码减法也满足交换群的性质。

2.2运算规律

令α=α1α2…αn,β=b1b2…bn,编码加法计算方式与十进制加法类似,从右侧末位码元(即an或bn)开始,向左侧首位码元(即a1或b1)逐位计算,可表示为


(4)

其中


(5)

i-1与n-i表示码元0的个数。当i=1时,e∈{0,1,2,3,4,5,6};当i>1时,e∈{0,1,2,3}。εiai⊕εibi计算结果的后n-i位都为码元0,这与十进制加法相同,计算至第i位时,其右侧n-i位均计算完毕,不参与运算。因此,式(5)在计算时可记作

对于首位编码码元(i=1),其运算规律以表1加法表的形式给出,表中左列与首行为对应相加码元。由于相加结果会超出原始7个码元集合,因此增加{10,20,30,40,50,60,100,200,300,400,500,600}12个首层格点编码,其复数平面位置可其相加码元的向量加法运算得出。

表17个码元⊕运算查找表

⊕012345600123456111002020610222020030301332303004040440340400505556045050060661010560600

对于非首位码元相加(i>1),除了得到当前位的计算结果码元,还会在左侧首位至第i-1位产生进位码元,具体规律如下


(6)

每一位进位码元再根据式(6)或表1相加,直至得到当前位最终计算结果码元。在计算过程中,层次i(1≤i<n)对应加法位上进位码元个数为


,每一层所做加法运算的次数也呈2的指数倍增长。为了提高运算效率,必须减少加法运算的次数。分析加法运算规律发现:

(1)0⊕e=e,0作为进位码元无累加意义,1,2,3作为进位码元有累加意义。

(2)当e1,e2∈{1,2,3},计算


e2(n≥2)。若e1=e2,产生一位非零进位;若e1≠e2,产生n-1位非零进位。

因此在计算过程中,不计算0并优先计算相同码元后,每层剩余需相加的码元不会超过3个,剩余最多情况时,所需相加的码元集合为{1,2,3}。

(3)不必开辟大量内存存放码元,只需统计各码元出现次数,提高计算效率。

此外,在计算过程中,还需要对首位编码溢出的情况加以控制:

(1)先将相加为0的编码消除,这样剩余编码个数不超过3个。

(2)因为没有建立码元{10,20,30,40,50,60}的加法表,为了避免其进行加法运算,对剩余编码,先进行间隔相加,如1⊕3=2;然后,将相同编码相加,如1⊕1=100。

图5(a)为编码加法运算示例2322⊕1033=100,201,其中黑色、灰色不带下划线和带下划线码元分别表示原始相加码元、计算结果码元和进位码元。图5(b)是根据上述运算规律所设计的算法流程。


图5编码加法运算示例与算法流程

3坐标与编码转换

采用一维格网编码代替传统二、三维笛卡尔坐标参与运算是全球离散格网系统具有较高数据处理效率的重要原因。尽管在格网系统中,单元操作可完全通过编码实现,但目前大部分空间数据仍采用地理坐标或投影坐标给出。为了保证现有数据进、出格网系统,必须建立笛卡尔实数坐标与单元编码码元序列之间的对应关系。

坐标到编码的转换本质是根据编码规则,记录逼近坐标点所经过各层次单元的编码,逼近路径的不同导致转换效率差异。本文借助编码运算实现的转换,以矢量相加的“平行四边形法则”为依据,避开了层次路径的迭代,具有直观、高效的特点。

3.1二维直角坐标转换为编码

给定二维直角坐标已知的点p(x,y),利用平行四边形法则,将点p投影到两条特定的坐标轴上,得到两条轴上对应的分量编码α1和α2,再通过编码加法运算,计算出点p所属单元编码α。

图6为二维直角坐标到编码的转换示意图,格网层次为2,编码前3位为首层编码,最后一位为第二层编码。引入三轴坐标系O-IJK作为投影坐标轴,原点与O-XY一致,J轴与Y轴一致,I、J、K轴之间正向夹角120°,按顺时针方向排列,将360°平面分成6个象限。


图6平面坐标与编码转换示意

投影坐标轴的选择可由点p和原点O连线与X轴的夹角θ的正切值来判断。I、J、K轴是过一列连续排列格点的直线,位于其上格点编码α的后n-1位编码N{n-1},与所在坐标轴和其编号k有如下关系


(7)

式中,DBin(k)代表编号k的二进制数序列,编号k为单元在轴上的位置,可由投影长度计算得到。αi的首位编码N0虽不符合二进制规律,但其与k之间的关系也很容易归纳。最后,根据编码加法运算规律,可得


(8)

3.2编码转换为二维直角坐标

格点本质上就是复平面上一种特殊形式的数,编码是其另一种有效的表达方式,因此编码与格点的复数平面坐标可以相互转化。将编码按复进制展开,计算其在复数平面上的实部坐标x和虚部坐标y,得到格点的二维直角坐标(x,y)。

假设格网层次为n,某单元编码α=e1e2…en,根据码元与数字位间的关系有


(9)

4试验与分析

为了验证以上算法的效率,笔者采用VisualC++2012开发工具分别实现了HLQT、PYXIS[11]和HQBS[16-17]的编码加法,以及HLQT和HQBS的二维直角坐标与编码相互转换算法。由于HLQT与HQBS方案编码总数会随着层次升高急剧增加,因此在加法算法中这两种方案随机选取固定个数编码,两两相加。全部程序编译为Release版本,并在一台兼容机(软件配置:Windows7x64旗舰版+SP1;硬件配置:IntelCorei5-6500CPU@3.20GHZ,8GRAM,KingSton120GBSSD)上进行测试。统计各算法6—11层的运算效率,加法运算与转换算法效率分别为单位时间内作加法运算的次数和完成转换的坐标或编码个数,取10次测试的平均值,结果见表2、表3、图7、图8所示。

表23种编码加法的运算效率统计结果

层次算法PYXIS
HQBS
HLQT单元个数运算效率/(次/ms)单元个数运算效率/(次/ms)单元个数运算效率/(次/ms)61261893
40961442
5000954074039136330001354500080118116058853000126550006734935839135130001260500058921010546987430001212500049461132050313403000117750004355

图7编码加法运算效率对比


图8二维直角坐标与编码转换算法效率对比

表3HQBS与HLQT两种转换效率测试结果

层次算法HQBS
HLQT(x,y)→α
α→(x,y)(x,y)→α
α→(x,y)坐标个数运算效率/(个/ms)编码个数运算效率/(个/ms)坐标个数运算效率/(个/ms)编码个数运算效率/(个/ms)62880002043
409615457
22202112755
40963978373348001958163841465122202110224163843961382142891840655361382322202184496553638448930500017292621441305122202176312621443655010565929153410485761245122202170581048576356231167500015484194304116832220216565419430434164

文献[17]验证了HQBS编码运算效率要优于PYXIS,这与图7的结果基本一致。因此本文在转换算法中,只对比HQBS与HLQT的效率。

试验结果表明:

(1)HLQT编码加法的平均运算效率约为PYXIS算法的6倍、HQBS算法的5倍。这是因为PYXIS奇(偶)分层编码,导致加法查找表比较复杂,而HLQT无须区分奇(偶)层,查找表规模较小(首位为7×7,其余位为4×4);HQBS单元中心与顶点混合编码,导致运算规则复杂且易出现编码“正则化”失败需要退回重算的情况,而HLQT只对格点编码,运算规律更直接,因而效率较高。

(2)HLQT的二维直角坐标向编码转换的平均运算效率大约是HQBS的5倍。这主要得益于HLQT的高率编码运算,且转换算法中位于I、J、K3轴上参与加法运算的HLQT编码在码元结构上具有二进制数特征,采用效率极高的二进制运算直接得到计算结果。

(3)HLQT的编码向二维直角坐标转换的平均运算效率是HQBS的3倍。因为HQBS算法涉及较多浮点数运算,而HLQT将浮点数运算处理为等效的整数运算,提高了运算效率。因码元个数随层次升高线性增加,二者效率均呈线性化下降趋势。

5总结与展望

本文根据平面四孔六边形格网系统结构特点,采用复进制数描述单元层次关系并设计了HLQT层次编码结构。该结构只对格点编码,码元具有明确几何意义,编码方案和编码运算具有严密的理论基础,克服了PYXIS奇(偶)分层编码、HQBS单元中心与顶点混合编码的不足,概念模型简单,编码运算和编码转换效率较高。为了建立四孔六边形全球离散格网的编码运算体系,下一步还需研究HLQT在球面上的扩展方案。

【引文格式】王蕊,贲进,杜灵瑀,等.平面四孔六边形格网系统编码运算[J].测绘学报,2018,47(7):1018-1025.DOI:10.11947/

院士论坛|杨元喜:弹性PNT基本框架

李德仁当选国际宇航科学院院士,集5院院士于一身!

论文推荐|邸凯昌:视觉SLAM技术的进展与应用

院士论坛|李德仁:遥感双院士的中国梦

论文推荐|王密:高分辨率光学卫星影像高精度在轨实时云检测的流式计算

《测绘学报》编委金双根当选欧洲科学院院士

学术前沿|唐新明:雷达卫星自动成图的精密干涉测量关键技术

《测绘学报》“数字摄影测量与机器视觉专辑”在CPGIS2018北京论坛发布

听李德仁、杨元喜、龚健雅三位院士讲述“我的科研故事”

机器视觉|晏磊:航空遥感平台通用物理模型及可变基高比系统精度评价

权威|专业|学术|前沿

微信投稿邮箱|song_qi_fan@163.com

进群请备注:姓名+单位+稿件编号