賁 進，童曉沖，元朝鵬

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

1 引 言

全球離散格網(wǎng)數(shù)據(jù)模型的基本思想是將地球遞歸剖分為形狀、面積近似相等且具有規(guī)則層次結(jié)構(gòu)的單元，每個單元在記錄位置信息的同時也表達(dá)比例尺和精度，因而具有處理多尺度數(shù)據(jù)的潛力。同一剖分產(chǎn)生的一系列不同分辨率的格網(wǎng)稱為“格網(wǎng)系統(tǒng)”，按照數(shù)學(xué)基礎(chǔ)的不同，可將其劃分為基于地理坐標(biāo)系的全球離散格網(wǎng)系統(tǒng)和基于多面體剖分的全球離散格網(wǎng)系統(tǒng)。

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

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

2 六邊形格網(wǎng)系統(tǒng)

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

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

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

圖1 孔徑為3的六邊形格網(wǎng)系統(tǒng)Fig.1 The aperture 3hexagonal grid system

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

圖2 孔徑為4的六邊形格網(wǎng)系統(tǒng)Fig.2 The aperture 4hexagonal grid system

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

3 OA4HDGGS的數(shù)學(xué)描述

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

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

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

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

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

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

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

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

圖3 格網(wǎng)上的操作示意圖Fig.3 Operations on grids

（1）對偶。對偶操作D（T）的頂點集合V［D（T）］＝C（T），當(dāng)且僅當(dāng)T中對應(yīng)單元共用一條邊時，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原點的正八面體或其在球面上的映射；格網(wǎng)集合序列遞歸定義為

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

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

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

4 OA4HDGGS的單元索引算法

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

圖5 坐標(biāo)系的定義Fig.5 Definitions of coordinate systems

定理1：h∈Hn中心（t∈Tn頂點）的格網(wǎng)坐標(biāo)滿足是h在R3中的坐標(biāo)。

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

該定理不僅建立了格網(wǎng)坐標(biāo)系和笛卡兒坐標(biāo)系之間的聯(lián)系，而且也說明（a，b，c）與Hn的單元一一對應(yīng)。

定理2：當(dāng)且僅當(dāng)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），其格網(wǎng)坐標(biāo)滿足定理3。

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

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

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

該定理描述了中心子單元的層次關(guān)系。

定理5：設(shè)h（a，b，c）∈Hn是鄰近子單元，則h有2個鄰近單元h′（d，e，f）的格網(wǎng)坐標(biāo)全是偶數(shù)，h在Hn－1上父單元是

例如，h（1，3，2）∈H2是鄰近子單元，根據(jù)定理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兩個單元的格網(wǎng)坐標(biāo)全是偶數(shù)，h在H1上有和個父單元。滿足定理5。

該定理描述了鄰近子單元的層次關(guān)系。

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

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

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

5 對比試驗與分析

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

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

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

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

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

造成以上現(xiàn)象的根本原因在于本文索引算法可以完全采用加減、移位等執(zhí)行效率極高的基本操作實現(xiàn)，而在計算機中沒有與BT對應(yīng)的數(shù)據(jù)類型，只能通過數(shù)組、結(jié)構(gòu)等模擬其運算，導(dǎo)致BT索引算法效率較低。

表2 試驗結(jié)果（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 結(jié) 論

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

［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.（陳軍，侯妙樂，趙學(xué)勝.球面四元三角網(wǎng)的基本拓?fù)潢P(guān)系描述和計算［J］.測繪學(xué)報，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.計算機程序設(shè)計藝術(shù)：半數(shù)值算法［M］.3版.蘇運霖，譯.北京：國防工業(yè)出版社，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.（劉百惠，杜荔，于濤.十進制數(shù)與對稱三進制數(shù)之間轉(zhuǎn)換的一種算法［J］.遼寧大學(xué)學(xué)報：自然科學(xué)版，1985（4）：29－35.）