孔徑為4的全球六邊形格網系統索引方法

賁 進，童曉沖，元朝鵬

1.信息工程大學測繪學院，河南鄭州450052；2.西安測繪信息技術總站，陜西西安710054

1 引 言

全球離散格網數據模型的基本思想是將地球遞歸剖分為形狀、面積近似相等且具有規則層次結構的單元，每個單元在記錄位置信息的同時也表達比例尺和精度，因而具有處理多尺度數據的潛力。同一剖分產生的一系列不同分辨率的格網稱為“格網系統”，按照數學基礎的不同，可將其劃分為基于地理坐標系的全球離散格網系統和基于多面體剖分的全球離散格網系統。

地理坐標系是最常使用的球面坐標系，在此基礎上建立的經緯格網系統得到廣泛應用。這類格網系統大都存在高緯地區單元變形嚴重的缺陷，而且矩形格網單元臨近關系的定義并不嚴格［1］?；诙嗝骟w剖分的全球離散格網系統是經緯格網的推廣，在空間分析、元胞自動機、氣候模擬等方面均有應用［2－4］。

與經緯格網相比，基于多面體剖分的全球離散格網系統具有以下優勢：① 以整個地球為研究對象且兩極不存在畸變，更適合處理全球尺度的問題；② 對地球表面同構離散化，有助于數據統一建模和重組；③ 格網具有層次性，在結構上支持多尺度數據表達；④ 空間運算可以采用編碼運算完成，更適合計算機處理。

2 六邊形格網系統

在構建格網常用的3種幾何圖形或剖分方法（三角形、四邊形和六邊形）中，六邊形格網的空間覆蓋效率和角度分辨率最高，單元一致相鄰且無共頂點的角鄰近單元，這些特性都有利于空間數據的處理和分析［1］。

然而，六邊形全球離散格網系統的相關問題一直都是學術界的研究難點，這是因為：① 六邊形不具備自相似性，一個大六邊形不能完全分解為若干小六邊形，導致不同分辨率格網之間的層次關系描述困難；② 采用六邊形不能完全覆蓋全球，在多面體頂點處必然出現其他形狀的單元，需要分別處理；③ 六邊形剖分的方法較多，幾乎不可能建立適合各種情況的“統一剖分模型”，只能分別討論。

孔徑（aperture），即相鄰層次格網單元的面積之比，是描述格網系統或剖分方法的重要指標?？讖皆叫?，相鄰層次格網的分辨率變化越均勻，對多分辨率應用越有利。六邊形剖分的最小孔徑是3（圖1），許多學者對這類格網系統進行了較深入的研究。文獻［5—6］將平面上的廣義平衡三進制（general balanced ternary，GBT）擴展到球面，提出基于正二十面體的六邊形格網編碼、索引算法；文獻［7—8］從理論上證明采用平衡三進制（balanced ternary，BT）可有效描述基于正八面體的六邊形格網層次關系；文獻［9］采用復平面的向量組合建立基于正二十面體的六邊形格網編碼方案；文獻［10］研究基于正二十面體的六邊形格網三角化算法?？讖綖?的六邊形格網系統存在的主要不足是單元的方向隨剖分層次交替變化，這導致相關算法較復雜。

圖1 孔徑為3的六邊形格網系統Fig.1 The aperture 3hexagonal grid system

孔徑為4的六邊形剖分方法不唯一，對應的格網系統均具有單元方向固定、層次結構簡單的優點。圖2是典型的孔徑為4的六邊形剖分，每個單元僅有1或2個父單元，這一特性有助于建立高效的格網索引。但是，目前學術界對這類格網系統的研究尚不深入。

圖2 孔徑為4的六邊形格網系統Fig.2 The aperture 4hexagonal grid system

以基于正八面體的、孔徑為4的六邊形格網系統（octahedron－based aperture 4hexagonal discrete global grid system，OA4HDGGS）為研究對象，從集合論的角度建立不同層次六邊形格網和三角形格網之間的遞推關系，借助三角形單元頂點建立六邊形單元層次關系并據此設計單元索引算法，通過對比試驗驗證該方法的效率。

3 OA4HDGGS的數學描述

從集合論的角度思考，格網系統是一系列不同分辨率格網的集合，也是一些單元的集合。盡管直接描述OA4HDGGS有難度，但是六邊形可以分解為三角形，能夠建立三角形格網集合之間的關聯就相當于建立了六邊形格網集合之間的關聯，這是描述OA4HDGGS的基本思想。為簡化表達，下文出現的符號及其約定見表1。

表1 符號約定Tab.1 The convention of signs

設正八面體中心和單位球（半徑為1）球心均在R3的原點，在正八面體表面或球面取頂點坐標為（x1，x2，x3）的三角形單元t的中心，即

取三角形格網T中所有單元的中心構成一個集合，記為C（T）＝｛θ（t）｜t∈T｝。

在正八面體表面或球面，取頂點坐標為（x1，x2，…，x6）的六邊形單元h的中心，即

h邊的集合記為ε（h），取任意一邊e（xi，xi＋1）（i＝1，2，…，5）的中點，記為

相鄰各邊中點連線構成的集合記為τ（h），即τ（h）＝｛l［μ（ei），μ（ei＋1）］｜ei∈h，ei＋1∈h｝，l［μ（ei），μ（ei＋1）］表示端點為μ（ei）、μ（ei＋1）的線段，下同。

h中心和各邊中點的連線構成的集合記為ρ（h），即ρ（h）＝｛l［θ（h），μ（ei）］｜ei∈h｝。取六邊形格網H中所有單元邊的中點構成一個集合，記為M（H）＝｛μ（h）｜h∈H｝。連接H中所有單元邊的中點構成一個集合，記為J（H）＝｛τ（h）｜h∈H｝。連接H中所有單元的中心和邊的中點構成一個集合，記為R（H）＝｛ρ（h）｜h∈H｝。為了建立三角形格網T和六邊形格網H之間的關聯，定義如下操作，如圖3。

圖3 格網上的操作示意圖Fig.3 Operations on grids

（1）對偶。對偶操作D（T）的頂點集合V［D（T）］＝C（T），當且僅當T中對應單元共用一條邊時，D（T）中兩個頂點之間有邊連接。

（2）中心剖分。中心剖分S（H）的頂點集合記為V［S（H）］＝C（H）∪M（H），S（H）邊的集合E［S（H）］＝R（H）∪J（H）。

在球面或八面體表面定義集合序列（Tn，Hn）（前4層如圖4）：T0是位于R3原點的正八面體或其在球面上的映射；格網集合序列遞歸定義為

式中，Hn即OA4HDGGS，n≥1是剖分層次。

采用數學歸納法不難證明如下結論：① 三角形格網Tn共有2×9×4n個單元；②格網Hn中共9×4n－4有個六邊形單元，頂點處共有6個四邊形單元；③Tn的頂點集合Vn∶＝V（Tn）嚴格滿足V0?V1?V2?…?Vn；④Hn的中心集合Cn∶＝C（Hn）嚴格滿足H0?H1?H2?…?Hn；⑤h∈Hn與Hn＋1中7個單元相交，其中有1個中心子單元，6個鄰近子單元；⑥h∈Hn與Hn－1中1個或2個單元相交，當h是中心子單元時是1個單元，當h是鄰近子單元時是2個單元，這些n－1層上的單元稱為h的父單元。

圖4 （Tn，Hn）在單個三角面上的示意圖Fig.4 （Tn，Hn）on a single triangular face

4 OA4HDGGS的單元索引算法

設t∈T0的頂點位于R3中（±1，0，0）、（0，±1，0）、（0，0，±1），h∈Hn的笛卡兒坐標不是整數，這給單元快速索引造成不便。為了解決這一問題，建立三軸格網坐標系，其單元坐標h（a，b，c）均為整數（圖5）。在三軸格網坐標系下，Hn滿足如下定理（具體證明與文獻［8］類似，本文不再贅述）。

圖5 坐標系的定義Fig.5 Definitions of coordinate systems

定理1：h∈Hn中心（t∈Tn頂點）的格網坐標滿足是h在R3中的坐標。

例如，h（2，2，2）∈H2滿足2＋2＋2＝3× 22－1＝6，其笛卡兒坐標

該定理不僅建立了格網坐標系和笛卡兒坐標系之間的聯系，而且也說明（a，b，c）與Hn的單元一一對應。

定理2：當且僅當Vn∶＝V（Tn）中兩頂點A（a1，a2，a3）和B（b1，b2，b3）滿足｜ai－bi｜≤1（i＝1，2，3）時，A、B兩點之間有邊相連。

例如，對A（2，2，2）∈T2滿足定理2的頂點有：B1（1，3，2）、B2（1，2，3）、B3（2，1，3）、B4（3，1，2）、B5（3，2，1）、B6（2，3，1），它們在T2中和A均有邊相連。

再如，對A（0，0，6）∈T2滿足定理2的頂點有：B1（1，0，5）、B2（0，1，5）、B3（－1，0，5）、B4（0，－1，5），它們在T2中和A均有邊相連。

該定理表明Hn中一個六邊形單元有6個鄰近單元，一個四邊形單元（位于八面體的6個頂點上）有4個鄰近單元。

定理3：hA（a1，a2，a3）∈Hn的鄰近單元hB（b1，b2，b3）可表示為｜ai－bi｜≤1（i＝1，2，3）。

例如，hA（2，2，2）∈H2的鄰近單元是：hB1（1，3，2）、hB2（1，2，3）、hB3（2，1，3）、hB4（3，1，2）、hB5（3，2，1）、hB6（2，3，1），其格網坐標滿足定理3。

因為Tn的頂點即Hn的單元中心，所以定理3與定理2等價。

定理4：h（a，b，c）∈Hn在Hn＋1上的中心子單元是h′（2a，2b，2c），當且僅當a、b、c同為偶數時，h（a，b，c）是中心子單元，其在Hn－1上的父單元是

例如，h（2，2，2）∈H2在H3上的中心子單元是h′（4，4，4），h是的中心子單元，滿足定理4。

該定理描述了中心子單元的層次關系。

定理5：設h（a，b，c）∈Hn是鄰近子單元，則h有2個鄰近單元h′（d，e，f）的格網坐標全是偶數，h在Hn－1上父單元是

例如，h（1，3，2）∈H2是鄰近子單元，根據定理3，h的鄰近單元是：h1（0，4，2）、h2（0，3，3）、h3（1，2，3）、h4（2，2，2）、h5（2，3，1）、h6（1，4，1），其中h1和h4兩個單元的格網坐標全是偶數，h在H1上有和個父單元。滿足定理5。

該定理描述了鄰近子單元的層次關系。

一般來說，離散全球格網上的基本索引操作包括：① 查找某一單元的鄰近單元；② 查找某一單元的子單元；③ 查找某一單元的父單元。在定理1～5的支撐下，很容易設計出OA4HDGGS上的單元索引算法。根據定理3，可計算出h（a，b，c）∈Hn的鄰近單元；根據定理3和定理4，可計算出h（a，b，c）∈Hn在Hn＋1上的中心子單元和鄰近子單元；根據定理4和定理5，可計算出h（a，b，c）∈Hn在Hn－1上的父單元。

在算法實現過程中，可采用“編碼壓縮”技術以節省數據存儲空間。在區間［－2m－1，2m－1）內的任意整數都可以用m個二進制位表示。對于h（a，b，c）∈Hn，因為a、b、c≤3×2n－1，所以每個坐標需要用n＋1個二進制位記錄。又根據定理1，c＝±（3×2n－1－｜a｜－｜b｜），因此在準確記錄a、b的前提下，只需再用一個二進制位記錄c的正負即可。這樣共需2n＋3個二進制位對Hn中的單元進行編碼。在32位Windows XP操作系統上，如果直接采用內置的int64整數類型，能夠支持30層格網的編碼，此時單元分辨率已達到毫米級，對于地球空間信息的表達和處理已足夠。

在上述索引算法中，鄰近單元查找只涉及整數（二進制數）的自加、自減運算，在計算機中的執行效率較高。父（子）單元查找也只涉及整數（二進制數）的2倍數運算，可通過位運算實現，執行效率更高。

5 對比試驗與分析

為了驗證上述索引算法的可行性和優越性，將其與文獻［8］提出的基于BT的索引算法進行對比。

第一組試驗首先在正八面體第一卦限的三角面上進行指定層次（5～18）、孔徑為3的六邊形剖分，包括三角面邊界和頂點上的單元，第n層格網上的單元數目是（n是奇數）是偶數）。然后采用BT算法查找該層次上每個單元的鄰近單元、子單元和父單元，統計平均執行效率。BT的轉換及運算采用文獻［11］中的算法實現。

第二組試驗在同一三角面上進行指定層次（5～18）、孔徑為4的六邊形剖分，包括三角面邊界和頂點上的單元，第n層格網上的單元數目是（3·2n＋2）（3·2n＋4）。然后采用本文的索引算法進行查找并統計平均執行效率。

所有程序均在Windows XP SP3系統上采用Visual C＋＋2008SP1編譯成Release版本后在一臺兼容機（處理器Intel Pentium Dual 3.0GHz；內存1.0GB DDR2－667SDRAM）上運行，結果如表2、圖6、圖7所示。

通過以上對比試驗不難發現，本文的索引算法與BT算法相比具有明顯優勢：首先，執行效率非常高，平均約是BT算法的600倍（表2）；其次，執行效率穩定（圖6），而BT算法的效率隨格網層次的增加急劇下降（圖7）。

造成以上現象的根本原因在于本文索引算法可以完全采用加減、移位等執行效率極高的基本操作實現，而在計算機中沒有與BT對應的數據類型，只能通過數組、結構等模擬其運算，導致BT索引算法效率較低。

表2 試驗結果（10層～18層）Tab.2 The result of experiments（level 10～level 18）

圖6 本文算法的效率曲線Fig.6 The efficiency curves of our algorithms

圖7 BT算法的效率曲線Fig.7 The efficiency curves of BT algorithms

6 結 論

采用編碼代替浮點數進行運算是全球離散格網數據模型的特色之一，因此單元編碼方案和對應索引算法的設計非常重要?；谡嗣骟w的、孔徑為4的六邊形離散全球格網的幾何性質，提出一種索引方法。與現有成果相比，這種編碼索引方案的優點體現在：① 采用集合論描述格網剖分，理論基礎嚴密，表達形式簡潔；② 建立了三角格網和六邊形格網之間的關聯，便于格網的三角化顯示；③單元編碼與其空間坐標的對應關系簡單，轉換方便；④ 單元層次、鄰近關系明確，索引算法簡單高效。

［1］ SAHR K，WHITE D，KIMERLING A J.Geodesic Discrete Global Grid Systems［J］.Cartography and Geographic Information Science，2003，30（2）：121－134.

［2］ CHEN Jun，HOU Miaole，ZHAO Xuesheng.Describing and Computing Model of the Topological Relation in Spherical Surface Quaternary Triangular Mesh［J］.Acta Geodaetica et Cartographica Sinica，2007，36（2）：176－180.（陳軍，侯妙樂，趙學勝.球面四元三角網的基本拓撲關系描述和計算［J］.測繪學報，2007，36（2）：176－180.）

［3］ KIESTER R，SAHR K.Planar and Spherical Hierarchical，Multi－resolution Cellular Automata［J］.Computers，Environment and Urban Systems，2008，32（3）：204－213.

［4］ HEIKES R，RANDALL D.Numerical Integration of the Shallow－water Equations on a Twisted Icosahedral Grid.Part I：Basic Design and Results of Tests［J］.Monthly Weather Review，1995，123（6）：1862－1880.

［5］ GIBSION L，LUCAS D.Spatial Data Processing Using Generalized Balanced Ternary［C］∥Proceedings of IEEE Computer Society Conference on Pattern Recognition and Image Processing.Las Vegas：IEEE，1982：566－571.

［6］ SAHR K.Discrete Global Grid Systems：A New Class of Geospatial Data Structure［D］.Oregon：The Graduate School of the University of Oregon，2005.

［7］ DONALD E K.The Art of Computer Programming：Seminumerical Algorithms［M］.3rd ed.Translated by SU Yunlin.Beijing：National Defence Industrial Press，2002.（DONALD E K.計算機程序設計藝術：半數值算法［M］.3版.蘇運霖，譯.北京：國防工業出版社，2002.）

［8］ VINCE A.Indexing the Aperture 3Hexagonal Discrete Global Grid［J］.Journal of Visual Communicatin and Image Representation，2006，17（6）：1227－1236.

［9］ ZHENG X Q.Efficient Fourier Transforms on Hexagonal Arrays［D］.Florida：University of Florida，2007.

［10］ MATTHEW G.Triangulation of a Hierarchical Hexagon Mesh［D］.Kingston：Queen＇s University，2009.

［11］ LIU Baihui，DU Li，YU Tao.A Method of Conversion between Decimal Number and Symmetrical Ternary Number［J］.Journal of Liaoning University，Natural Sciences Edition，1985（4）：29－35.（劉百惠，杜荔，于濤.十進制數與對稱三進制數之間轉換的一種算法［J］.遼寧大學學報：自然科學版，1985（4）：29－35.）