王 蕊，賁 進，杜靈瑀，周建彬，李祝鑫

1. 信息工程大學(xué)地理空間信息學(xué)院,河南 鄭州 450001； 2. 61287部隊，四川 成都610000

全球離散格網(wǎng)系統(tǒng)(discrete global grid system，DGGS)是一種以整個地球為研究對象、以地理空間位置數(shù)字化表達為主要特征的新型數(shù)據(jù)模型。它將地表剖分為均勻統(tǒng)一的層次格網(wǎng)單元結(jié)構(gòu)，并給每個單元賦予全局唯一的地址編碼代替?zhèn)鹘y(tǒng)地理坐標(biāo)參與運算和操作[1]。與傳統(tǒng)空間數(shù)據(jù)模型相比，其優(yōu)勢體現(xiàn)在：面向全球，能夠提供全球一致的空間基準(zhǔn)；天然多尺度，有助于處理多源、異構(gòu)地理空間數(shù)據(jù)；純編碼運算，有助于提高數(shù)據(jù)處理效率。

與三角形、四邊形相比，六邊形具有拓?fù)潢P(guān)系一致、采樣效率高等特點，有助于提高空間分析算法精度[2-3]和可視化表達效果[4-5]，其主要缺陷是不具備自相似性，無法直接建立多分辨率格網(wǎng)單元層次關(guān)系。文獻[6]提出一種多維數(shù)據(jù)結(jié)構(gòu)，二維特例可通過“廣義平衡三進制”(generalized balanced ternary，GBT)運算實現(xiàn)不同層次六邊形單元的編碼索引[7]。文獻[8]改進GBT的不足，提出HIP(hexagonal image processing)結(jié)構(gòu)，并將其應(yīng)用于影像處理。文獻[9—11]區(qū)分奇(偶)層建立了平面三孔六邊形格網(wǎng)系統(tǒng)的單元編碼索引并將其擴展到球面，在此基礎(chǔ)上，文獻[12]實現(xiàn)了傅里葉變換算法。但這類方案對奇(偶)層單元分別處理，增加了編碼運算復(fù)雜度，索引效率不高。文獻[3,13—15]采用“單元隸屬圖形”建立六邊形格網(wǎng)的編碼索引，主要不足是破壞了六邊形單元的完整性。文獻[16—17]借助格網(wǎng)單元中心與頂點建立了平面四孔六邊形格網(wǎng)系統(tǒng)的“四元平衡結(jié)構(gòu)”(hexagonal quad balanced structure，HQBS)并將其擴展到球面，主要缺陷是概念模型復(fù)雜，且會在編碼運算中會出現(xiàn)“正則化”失敗需要回退重算的情況，嚴(yán)重影響索引效率[18]。

綜上所述，單元編碼索引是六邊形全球離散格網(wǎng)系統(tǒng)研究的核心問題之一，平面格網(wǎng)系統(tǒng)層次關(guān)系描述及編碼方案設(shè)計則是該問題的突破口。現(xiàn)有成果仍然存在不足，亟須改進完善。本文根據(jù)平面四孔六邊形格網(wǎng)結(jié)構(gòu)特點，設(shè)計“六邊形格點四叉樹”(hexagon lattice quad tree，HLQT)層次編碼結(jié)構(gòu)，定義編碼運算等效代替格網(wǎng)單元的幾何操作，歸納編碼運算規(guī)律并實現(xiàn)了對應(yīng)的算法。借助編碼運算，設(shè)計并實現(xiàn)二維直角坐標(biāo)與單元編碼的相互轉(zhuǎn)換算法，通過對比試驗驗證本文方案的正確性與優(yōu)越性。

1 編碼方案設(shè)計

為了便于討論，本文將平面離散格網(wǎng)系統(tǒng)的單元中心定義為“格點”，作為單元的等效替代對象，四孔六邊形剖分產(chǎn)生的格點以及剖分規(guī)律如圖1所示：第n(n≥1)層格點是由第n-1層格點與n-1層單元各邊中點組合而成，第n-1層單元邊長是第n層單元邊長的2倍，單元方向不隨層次變化。

圖1 四孔六邊形格網(wǎng)系統(tǒng)格點與剖分規(guī)律示意圖Fig.1 Diagram of lattice points and subdivision rules for the aperture 4 hexagonal grid

格點位置通常采用笛卡爾坐標(biāo)描述[10,14,17,19-20]，這種方法雖然直觀形象、便于理解，但因采用了二維獨立變量(如x、y)，表達形式不簡潔且不易與編碼建立聯(lián)系。故本文采用復(fù)數(shù)表示格點位置[21-22]，建立平面四孔六邊形格點系統(tǒng)的唯一性描述方案[23]。

(1)

以上描述方案具有唯一性，生成的子單元互不重疊，每個子單元具有唯一的父單元。Ln是典型的四叉樹層次結(jié)構(gòu)，即隨著剖分層次的深入，L1中的每個格點都分別生長為一顆獨立的四叉樹，這為格網(wǎng)編碼提供了理論依據(jù)。

圖2 L1與L2格點層次關(guān)系示意Fig.2 Diagram of hierarchical relationship between L1 and L2

以上方案還表明：四孔六邊形格點系統(tǒng)與實數(shù)域上的二進制、八進制、十進制等定位計數(shù)系統(tǒng)的本質(zhì)完全相同，是復(fù)數(shù)域上一種特定形式的“數(shù)”[24]。該定位計數(shù)系統(tǒng)的進制為2，數(shù)字位集合為D′和D。對于x∈Ln，其具體形式為

d2,…,dn∈D)

(2)

或與十進制數(shù)類似，省略冪指數(shù)的底，式(2)可記為

x=.d1d2…dn(d1∈D′,d2,…,dn∈D)

(3)

式(3)中的小數(shù)點表示格點隨剖分層次向內(nèi)部生長，越來越密集。用不同大小和顏色的點表示L1～L5中的元素，它們在復(fù)平面上的分布如圖3 所示。

為了便于計算機處理，需要對格網(wǎng)系統(tǒng)中的每個單元指定一個地址碼，這一過程稱為“編碼”。編碼方案必須滿足完整性、唯一性與層次性[19]。

根據(jù)式(3)，取di冪指數(shù)，稱為碼元，即d1=ω′e→e(1≤e≤6)，di=ωe→e(1≤e≤3,2≤i≤n)，則數(shù)字位集合D′替換為碼元集合E′={0,1,2,3,4,5,6}，數(shù)字位集合D替換為碼元集合E={0,1,2,3}。得到x∈Ln的唯一編碼標(biāo)識e1e2…en(e1∈E′,ei∈E,2≤i≤n)。顯然，該編碼方案還滿足完整性和層次性的要求。

圖3 L1～L5在復(fù)平面上的分布Fig.3 Distribution of L1～L5 in the complex plane

L3的完整編碼如圖4所示，不同格點層次隸屬和關(guān)聯(lián)用不同大小的圓和三角形表示。將這一格點層次結(jié)構(gòu)命名為“六邊形格點四叉樹”(hexagon lattice quad tree，HLQT)。

圖4 L3編碼及運算示意圖Fig.4 Representation of codes in the 3rd level and examples of code arithmetic

2 編碼運算規(guī)律

HLQT編碼記錄了某一格點與原點(編碼全部為0的格點)的距離和方位，是復(fù)數(shù)平面上“定位計數(shù)系統(tǒng)”的等效表示。因此格點對應(yīng)復(fù)數(shù)之間的各類運算可從二維復(fù)數(shù)平面空間降維到一維編碼空間中實現(xiàn)，且編碼相較于浮點數(shù)更適合計算機處理，可有效提高計算效率。

2.1 定義

在復(fù)數(shù)平面上，任意兩復(fù)數(shù)的加、減法滿足平行四邊形法則。在HLQT編碼集合n(n≥1)中，設(shè)編碼α,β∈n，定義算子C(…)為計算給定HLQT編碼對應(yīng)復(fù)數(shù)，算子H(…)為計算給定復(fù)數(shù)的HLQT編碼，則α、β的編碼加法定義為復(fù)數(shù)C(α)+C(β)對應(yīng)的HLQT編碼H[C(α)+C(β)]，表示為γ=α⊕β=H[C(α)+C(β)]。如圖4所示，α⊕β=123⊕010=233。根據(jù)歸納、驗證，編碼加法滿足交換群的性質(zhì)[25]。HLQT的減法運算滿足向量相減的平行四邊形法則，例如：α?β=123?010=033。⊕與?是逆運算，編碼減法也滿足交換群的性質(zhì)。

2.2 運算規(guī)律

令α=α1α2…αn，β=b1b2…bn，編碼加法計算方式與十進制加法類似，從右側(cè)末位碼元(即an或bn)開始，向左側(cè)首位碼元(即a1或b1)逐位計算，可表示為

α⊕β=a1a2…an⊕b1b2…bn=

(4)

其中

(5)

對于首位編碼碼元(i=1)，其運算規(guī)律以表1 加法表的形式給出，表中左列與首行為對應(yīng)相加碼元。由于相加結(jié)果會超出原始7個碼元集合，因此增加{10,20,30,40,50,60,100,200,300,400,500,600}12個首層格點編碼，其復(fù)數(shù)平面位置可其相加碼元的向量加法運算得出。

表1 7個碼元⊕運算查找表

對于非首位碼元相加(i>1)，除了得到當(dāng)前位的計算結(jié)果碼元，還會在左側(cè)首位至第i-1位產(chǎn)生進位碼元，具體規(guī)律如下

(6)

(1) 0⊕e=e，0作為進位碼元無累加意義，1,2,3作為進位碼元有累加意義。

因此在計算過程中，不計算0并優(yōu)先計算相同碼元后，每層剩余需相加的碼元不會超過3個，剩余最多情況時，所需相加的碼元集合為{1,2,3}。

(3) 不必開辟大量內(nèi)存存放碼元，只需統(tǒng)計各碼元出現(xiàn)次數(shù)，提高計算效率。

此外，在計算過程中，還需要對首位編碼溢出的情況加以控制：

(1) 先將相加為0的編碼消除，這樣剩余編碼個數(shù)不超過3個。

(2) 因為沒有建立碼元{10,20,30,40,50,60}的加法表，為了避免其進行加法運算，對剩余編碼，先進行間隔相加，如1⊕3=2；然后，將相同編碼相加，如1⊕1=100。

圖5(a)為編碼加法運算示例2322⊕1033=100，201，其中黑色、灰色不帶下劃線和帶下劃線碼元分別表示原始相加碼元、計算結(jié)果碼元和進位碼元。圖5(b)是根據(jù)上述運算規(guī)律所設(shè)計的算法流程。

圖5 編碼加法運算示例與算法流程Fig.5 Example and flow chart for add operation

3 坐標(biāo)與編碼轉(zhuǎn)換

采用一維格網(wǎng)編碼代替?zhèn)鹘y(tǒng)二、三維笛卡爾坐標(biāo)參與運算是全球離散格網(wǎng)系統(tǒng)具有較高數(shù)據(jù)處理效率的重要原因。盡管在格網(wǎng)系統(tǒng)中，單元操作可完全通過編碼實現(xiàn)，但目前大部分空間數(shù)據(jù)仍采用地理坐標(biāo)或投影坐標(biāo)給出。為了保證現(xiàn)有數(shù)據(jù)進、出格網(wǎng)系統(tǒng)，必須建立笛卡爾實數(shù)坐標(biāo)與單元編碼碼元序列之間的對應(yīng)關(guān)系。

坐標(biāo)到編碼的轉(zhuǎn)換本質(zhì)是根據(jù)編碼規(guī)則，記錄逼近坐標(biāo)點所經(jīng)過各層次單元的編碼，逼近路徑的不同導(dǎo)致轉(zhuǎn)換效率差異。本文借助編碼運算實現(xiàn)的轉(zhuǎn)換，以矢量相加的“平行四邊形法則”為依據(jù)，避開了層次路徑的迭代，具有直觀、高效的特點。

3.1 二維直角坐標(biāo)轉(zhuǎn)換為編碼

給定二維直角坐標(biāo)已知的點p(x,y)，利用平行四邊形法則，將點p投影到兩條特定的坐標(biāo)軸上，得到兩條軸上對應(yīng)的分量編碼α1和α2，再通過編碼加法運算，計算出點p所屬單元編碼α。

圖6為二維直角坐標(biāo)到編碼的轉(zhuǎn)換示意圖，格網(wǎng)層次為2，編碼前3位為首層編碼，最后一位為第二層編碼。引入三軸坐標(biāo)系O-IJK作為投影坐標(biāo)軸，原點與O-XY一致，J軸與Y軸一致，I、J、K軸之間正向夾角120°，按順時針方向排列，將360°平面分成6個象限。

圖6 平面坐標(biāo)與編碼轉(zhuǎn)換示意Fig.6 Transformation between geographic coordinates and codes

投影坐標(biāo)軸的選擇可由點p和原點O連線與X軸的夾角θ的正切值來判斷。I、J、K軸是過一列連續(xù)排列格點的直線，位于其上格點編碼α的后n-1位編碼N{n-1}，與所在坐標(biāo)軸和其編號k有如下關(guān)系

(7)

式中，DBin(k)代表編號k的二進制數(shù)序列，編號k為單元在軸上的位置，可由投影長度計算得到。αi的首位編碼N0雖不符合二進制規(guī)律，但其與k之間的關(guān)系也很容易歸納。最后，根據(jù)編碼加法運算規(guī)律，可得

α=α1⊕α2

(8)

3.2 編碼轉(zhuǎn)換為二維直角坐標(biāo)

格點本質(zhì)上就是復(fù)平面上一種特殊形式的數(shù)，編碼是其另一種有效的表達方式，因此編碼與格點的復(fù)數(shù)平面坐標(biāo)可以相互轉(zhuǎn)化。將編碼按復(fù)進制展開，計算其在復(fù)數(shù)平面上的實部坐標(biāo)x和虛部坐標(biāo)y，得到格點的二維直角坐標(biāo)(x,y)。

假設(shè)格網(wǎng)層次為n，某單元編碼α=e1e2…en，根據(jù)碼元與數(shù)字位間的關(guān)系有

C(α)=C(e1e2…en)=

2nω′e1+2n-1ωe2+…+2ωen=

(9)

4 試驗與分析

為了驗證以上算法的效率，筆者采用Visual C++ 2012開發(fā)工具分別實現(xiàn)了HLQT、PYXIS[11]和HQBS[16-17]的編碼加法，以及HLQT和HQBS的二維直角坐標(biāo)與編碼相互轉(zhuǎn)換算法。由于HLQT與HQBS方案編碼總數(shù)會隨著層次升高急劇增加，因此在加法算法中這兩種方案隨機選取固定個數(shù)編碼，兩兩相加。全部程序編譯為Release版本，并在一臺兼容機(軟件配置：Windows 7 x64旗艦版+SP1；硬件配置：Intel Core i5-6500 CPU@3.20 GHZ，8 G RAM，KingSton 120 GB SSD)上進行測試。統(tǒng)計各算法6—11層的運算效率，加法運算與轉(zhuǎn)換算法效率分別為單位時間內(nèi)作加法運算的次數(shù)和完成轉(zhuǎn)換的坐標(biāo)或編碼個數(shù)，取10次測試的平均值，結(jié)果見表2、表3、圖7、圖8所示。

表2 3種編碼加法的運算效率統(tǒng)計結(jié)果

圖7 編碼加法運算效率對比Fig.7 Efficiency comparison of code addition

文獻[17]驗證了HQBS編碼運算效率要優(yōu)于PYXIS，這與圖7的結(jié)果基本一致。因此本文在轉(zhuǎn)換算法中，只對比HQBS與HLQT的效率。

圖8 二維直角坐標(biāo)與編碼轉(zhuǎn)換算法效率對比Fig.8 Efficiency comparison of transform algorithms

表3 HQBS與HLQT兩種轉(zhuǎn)換效率測試結(jié)果

試驗結(jié)果表明：

(1) HLQT編碼加法的平均運算效率約為PYXIS算法的6倍、HQBS算法的5倍。這是因為PYXIS奇(偶)分層編碼，導(dǎo)致加法查找表比較復(fù)雜，而HLQT無須區(qū)分奇(偶)層，查找表規(guī)模較小(首位為7×7，其余位為4×4)；HQBS單元中心與頂點混合編碼，導(dǎo)致運算規(guī)則復(fù)雜且易出現(xiàn)編碼“正則化”失敗需要退回重算的情況，而HLQT只對格點編碼，運算規(guī)律更直接，因而效率較高。

(2) HLQT的二維直角坐標(biāo)向編碼轉(zhuǎn)換的平均運算效率大約是HQBS的5倍。這主要得益于HLQT的高率編碼運算，且轉(zhuǎn)換算法中位于I、J、K3軸上參與加法運算的HLQT編碼在碼元結(jié)構(gòu)上具有二進制數(shù)特征，采用效率極高的二進制運算直接得到計算結(jié)果。

(3) HLQT的編碼向二維直角坐標(biāo)轉(zhuǎn)換的平均運算效率是HQBS的3倍。因為HQBS算法涉及較多浮點數(shù)運算，而HLQT將浮點數(shù)運算處理為等效的整數(shù)運算，提高了運算效率。因碼元個數(shù)隨層次升高線性增加，二者效率均呈線性化下降趨勢。

5 總結(jié)與展望

本文根據(jù)平面四孔六邊形格網(wǎng)系統(tǒng)結(jié)構(gòu)特點，采用復(fù)進制數(shù)描述單元層次關(guān)系并設(shè)計了HLQT層次編碼結(jié)構(gòu)。該結(jié)構(gòu)只對格點編碼，碼元具有明確幾何意義，編碼方案和編碼運算具有嚴(yán)密的理論基礎(chǔ)，克服了PYXIS奇(偶)分層編碼、HQBS單元中心與頂點混合編碼的不足，概念模型簡單，編碼運算和編碼轉(zhuǎn)換效率較高。為了建立四孔六邊形全球離散格網(wǎng)的編碼運算體系，下一步還需研究HLQT在球面上的擴展方案。