杜靈瑀,馬秋禾,賁 進,王 蕊
信息工程大學地理空間信息學院,河南 鄭州 450001
格網系統(tǒng)是一種多分辨率柵格層次數據結構,本質是通過規(guī)則格網單元的分解與聚合將連續(xù)平面映射到多層次格網中,且每個單元具有包含層次關系的唯一編碼[1]。格網系統(tǒng)的離散性使其更適合計算機處理,有助于多源、異構、多尺度空間數據的融合及多要素的空間疊加和復合分析[2]。三角形格網系統(tǒng)廣泛應用于地形建模[3-5]、空間數據索引與可視化[6]、地圖綜合模型[7];矩形(四邊形)格網系統(tǒng)廣泛應用于資源環(huán)境綜合科學調查[8]、遙感影像數據存儲與管理[9]及地圖多級顯示[10]。
相比于上述格網,六邊形格網具有形態(tài)緊湊、平面覆蓋率和角度分辨率高及一致相鄰等特性[11],更有利于空間數據的組織、處理和分析。然而,六邊形不具備自相似性,格網層次關系描述及計算是研究難點之一[12]。文獻[13—14]提出“廣義平衡三進制”(generalized balanced ternary,GBT),將高分辨率六邊形單元聚合為低分辨率單元,但其相鄰層次單元的方向連續(xù)旋轉,導致相關應用較復雜[15]。文獻[1,16—17]分別設計了三孔六邊形格網系統(tǒng)編碼方案,并借助改進的GBT實現不同層次單元快速索引,其編碼運算無法完全用二進制優(yōu)化,效率不甚理想[18]。文獻[12]以“格元—格點—格邊—格心”概念模型[19]為基礎,建立了三孔六邊形格網系統(tǒng)編碼方案,并揭示其數學本質。文獻[18]針對四孔六邊形格網系統(tǒng)設計了“六邊形四元平衡結構”(hexagonal quaternary balanced structure,HQBS),但在編碼運算過程中頻繁出現“正則化”失敗回退重算的情況[12]。針對這一缺陷,文獻[20]設計了“六邊形格點四叉樹”(hexagon lattice quad tree,HLQT),其編碼運算效率優(yōu)于HQBS及PYXIS,但其碼元分布不對稱,不利于鄰近單元查詢。文獻[21]采用菱形塊四叉樹組織球面六邊形格網層次結構,并提出格網編碼方案。在實踐應用方面,ESRI公司在ArcGIS GeoAnalytics Server中采用六邊形格網系統(tǒng)分析、展示數據(http:∥www.esrichina.com.cn/openclass-2017.html)。UBER公司在其流處理系統(tǒng)中采用六邊形格網系統(tǒng)實時采集、分析車輛信息(http:∥www.infoq.com/cn/presentations/uber-stream-processing-system-and-practice)。歐空局的SMOS衛(wèi)星也采用六邊形格網存儲、處理影像數據[22]。
綜上所述,當前六邊形格網系統(tǒng)的研究已取得可喜進展,但現有成果仍存在編碼方案不完善、運算效率待提升的問題。本文引入復進制數理論,依據間隔層次單元隸屬關系,建立平面四孔六邊形格網系統(tǒng)數學模型,在此基礎上提出具有結構對稱性的等效單元編碼方案、定義編碼運算并歸納運算規(guī)則、設計編碼索引及編碼與笛卡兒坐標互換算法,嘗試提高編碼操作效率。
復進制定位計數系統(tǒng)(complex radix positional system)由計算機科學泰斗Donald E.Knuth于1960年提出[23],它的基數(進制)、數字位集合元素可以是復數,本質上是二進制、十進制等常見定位計數系統(tǒng)在復數域上的推廣。復數在復進制定位計數系統(tǒng)中的表示形式稱為復進制數(complex radix number)。在全球離散格網的相關研究中,通常采用笛卡兒坐標描述單元位置[24],該方式雖然直觀,但與編碼關聯性不強。研究表明,以復進制數形式表示單元位置,更有助于編碼設計及其運算規(guī)律推導[25-26]。
作為位置、對象及數據與格網單元的關聯點,通常利用格心代替格網單元[19]。四孔六邊形格網的剖分規(guī)律如圖1所示,第n層格心由第n-1層格心及單元各邊中心組成,相鄰層次單元邊長為2倍關系,且單元方向不隨層次改變。
圖1 四孔六邊形格網系統(tǒng)示意圖Fig.1 Diagram of aperture 4 hexagon grid system
式中,“∪”、“+”是集合間的并運算、加法運算符。
證明:設第n-1層格網單元外接圓半徑為rn-1,第n層格心集合為Ln,根據格網系統(tǒng)生成規(guī)律,Ln可表示為
所以
此時有
(1)
不同于HLQT、HQBS等四孔六邊形格網系統(tǒng)逐層編碼方案以及PYXIS三孔六邊形格網系統(tǒng)奇、偶層分別編碼方案[1],這里的定位計數系統(tǒng)建立在層次遞推關系的基礎上。隨層次增加,格心呈現向內凹陷加密的趨勢。如圖4所示,用不同大小和顏色的點表示L1—L5中的格點,圖中空缺部分可由首層相鄰單元生成的子格點填補。
由上述規(guī)則直接給出L1與L2的編碼,如圖5所示。以L1、L2和集合D的編碼為基礎,依據層次遞推關系建立Ln(n≥3)編碼,由Ln-1得到Ln中繼承子單元的過程編碼層次增加但格心位置不變,因此規(guī)定Ln(n≥3)中繼承子單元編碼由Ln-1編碼末尾添加“00”得到,如圖5L3編碼中黑色單元所示。剖分子單元的格心由第n-2層格網單元依據集合D剖分得到,但第n層與第n-2層并非相鄰層次,因此在Ln-2編碼末尾添加“00”表示跨層后再添加集合D對應的等效碼元,如圖5L3編碼中紅色單元所示。
根據上述編碼規(guī)則,每層碼元由2位字符構成,由于單元歸屬明確,單元編碼具有唯一性;復進制數集合D中的元素關于原點對稱,編碼具有對稱性,有利于編碼索引操作。
圖2 基于層次遞推的格網系統(tǒng)層次關系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
編碼記錄了格網單元相對原點的距離和方位,可通過一維編碼運算實現二維復平面內的幾何操作,達到降維目的。相較于浮點數,編碼更適合計算機存儲與處理,有利于提高運算效率。
與矢量加法類似,可利用平行四邊形法則定義編碼加法運算[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。將加法表結果統(tǒng)一記為8位,則加法運算的進位操作僅存在于相鄰基礎碼元。
表1 編碼加法查找略表
圖6 錯位相加法示意圖Fig.6 Diagram of staggered addition
對于第5層編碼加法0000010003⊕0000000003,錯位相加法得到的結果為0000010300,而實際結果為0000001000。分析可知,錯位相加法以單元0000010003為中心加上0000000003得到最終結果,而實際編碼0000001000由第4層繼承得到,其生成中心為原點單元,中心單元的不一致導致計算結果的不一致。出現該問題的編碼存在連續(xù)非“00”編碼的明顯特征,以編碼α=α1α2…α2n-1α2n為例,αk-3αk-2αk-1αk是出現連續(xù)非“00”編碼的位置,將α1α2…αk轉換為α1α2…αk-3αk-200⊕00…00αk-1αk的結果,并在末尾加上αk后的編碼αk+1…α2n可將編碼α轉化為正確編碼,避免“回退”操作。
編碼乘法對應旋轉和縮放操作[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部分繼承子單元和剖分子單元的編碼規(guī)律即可實現子單元的查找。父單元查找也分為鄰層和隔層兩種情況,對于第n(n≥3)層編碼α=α1α2…α2n-1α2n,若α2n-1α2n=00,則該單元繼承于第n-1層,將末尾編碼“00”去掉可得到鄰層父單元。若α2n-1α2n≠00,則該單元由第n-2層的單元剖分得到,將末尾4位編碼去掉可得到隔層父單元。
本文編碼方案可借助十六進制數實現。碼元{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在編碼結構和編碼運算效率方面具有較大優(yōu)勢[18,20]。因此,選擇上述兩編碼方案進行對比試驗。本文采用Visual Studio 2012開發(fā)工具實現本文編碼方案對應的編碼加法、鄰近單元查詢、坐標與編碼互換算法,并與HQBS及HLQT相應算法比較。全部代碼均編譯為Release版本,并在相同臺式機(硬件配置:Intel Core i5-4590M CPU @ 3.30 GHz,8 GB RAM;操作系統(tǒng):Windows 7 x64旗艦版)上測試。4—11層中逐層隨機生成2500個單元,不重復相加,得到3 123 750組加法運算實例;隨機生成1 000 000組坐標作為笛卡兒坐標到編碼轉換試驗數據;第4層隨機選取256個單元,5—11層以4倍遞增確定隨機選取單元個數作為編碼到坐標轉換及鄰近單元查詢試驗數據。以單位時間內完成操作的次數為效率指標,統(tǒng)計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格心與頂點混合編碼,運算規(guī)律復雜且易出現編碼“正則化”失敗回退重算的情況。HLQT逐層進行碼元相加,計算非首位碼元相加時,除了得到當前位計算結果碼元,還會在左側所有位產生進位碼元。雖然相比于HQBS和HLQT,本文方案的加法表規(guī)模較大。但任何編碼方案中的加法表均以二維數組形式存儲,通過碼元在二維數組中的位置,可直接定位相應結果,因此加法表的規(guī)模幾乎不會影響其運算效率。同時,本文方案以每兩層為基本計算單元,且進位只發(fā)生在相鄰基本單元之間,相比之下計算量更小,因而效率更高。
圖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)中帶來的諸多優(yōu)勢相比,該轉換算法約20%的效率損失可以接受。
(3) 本文方案編碼轉換為笛卡兒坐標的平均效率約為HQBS的6.29倍、HLQT的2.36倍。HQBS算法涉及較多的浮點數運算,本文方案與HLQT均采用整數代替浮點數運算,因此兩者效率均高于HQBS。相較于HLQT,由于本文方案的編碼以“00”和非“00”碼元交替形式存在,計算中只需考慮非“00”碼元,因而效率更高。
本文結合平面四孔六邊形格網結構特點,引入復進制數理論,通過間隔層次格網單元隸屬關系建立數學模型,以此為基礎提出等效編碼方案、定義編碼運算并歸納運算規(guī)則、設計編碼索引及編碼與坐標相互轉換算法。相較于同類成果,本文編碼方案具有結構對稱性,在編碼加法、鄰近單元查詢及本方案編碼轉換為笛卡兒坐標等方面效率優(yōu)勢明顯,有助于各類格網處理、分析算法的實現,具有較大的應用潛力。
本文僅在平面內建立了編碼方案并實現了相關操作,借助多面體及合適的映射關系可將本文編碼方案擴展到球(參考橢球)面,實現全球四孔六邊形離散格網系統(tǒng)的構建。在此過程中涉及的相關問題,將是筆者下一步研究的重點內容。