吳宇豪，曹雪峰，余岸竹，孫萬忠

信息工程大學地理空間信息學院，河南 鄭州 450001

全球離散網格系統(tǒng)(discrete global grid system,DGGS)是一種面向全球的新型多分辨率數(shù)據(jù)融合與地學模擬解決方案[1-3]，它借助特定方法對地球空間進行均勻遞歸剖分，構建多分辨率的無縫不重疊網格結構，在數(shù)據(jù)操作中使用格元編碼代替?zhèn)鹘y(tǒng)的地理空間坐標[4]。空間填充曲線以分形幾何、代數(shù)等數(shù)學理論為支撐，刻畫被遍歷空間中各點在曲線上對應位置之間的關系及其隨階數(shù)的變化，具備構建抽象、嚴格、一般性的數(shù)學模型的理論潛力[5]。空間填充曲線編碼較好地顧及了格元鄰近關系表達的需求，在執(zhí)行鄰域操作時效率較高[6]。

研究人員經過廣泛研究對比后發(fā)現(xiàn)，在眾多空間填充曲線中，Hilbert曲線具有最好的空間聚簇性[7-8]，即多維空間中位置相鄰或相近的空間目標映射到Hilbert曲線上后能夠保持最佳的鄰近關系，可有效提高多維數(shù)據(jù)在磁盤等一維物理存儲中的訪問效率[9]，已被廣泛應用于網格單元編碼研究領域[10]。文獻[11]提出剖分影像金字塔模型及建塔策略，利用網格單元編碼的唯一性和Hilbert曲線空間連續(xù)性對影像塊進行索引與組織。文獻[12]采用Hilbert曲線構建菱形網格單元編碼模型，并配合提出Hilbert碼與地理坐標的相互轉換算法。文獻[13]結合其八叉樹結構，設計立體網格單元Hilbert編碼方法，并針對多分辨率數(shù)據(jù)組織要求，提出緊致Hilbert索引[14]。文獻[15]將地球圈層空間網格的研究范圍擴展至地月空間，按照類似方法構建了地月圈層網格[16]，分析給出了圈層網格Hilbert碼與“行-列-高”碼的轉換方法，鄰近圈層單元查找模型以及圈層單元拓撲模型[17]。文獻[18]通過層次嵌套細分構建了規(guī)則的多分辨率時空網格，并基于Hilbert曲線提出一種時空數(shù)據(jù)索引方法。文獻[19]面向時空數(shù)據(jù)設計了四維Hilbert碼索引，并提出一種查詢窗口與Hilbert碼段的轉換方式。文獻[20]針對時空軌跡數(shù)據(jù)設計了一種時空網格系統(tǒng)，提出基于Hilbert曲線的時空數(shù)據(jù)索引進行軌跡數(shù)據(jù)管理。

格元Hilbert碼的計算效率是影響全球離散網格系統(tǒng)中數(shù)據(jù)組織效率的關鍵之一，包括笛卡兒坐標(或地理坐標)至Hilbert碼的轉換計算以及鄰近格元Hilbert碼的計算等。笛卡兒坐標至Hilbert碼轉換計算的經典方法是由文獻[21]針對二維空間構造提出的。該算法是基于Morton碼的二進制位迭代法，算法復雜度為O(n2)。文獻[22]對n維m階Hilbert曲線的空間變換進行研究，提出基于空間坐標變換、自底向上的迭代生成算法，算法復雜度為O(nm)。文獻[23]提出狀態(tài)轉移矩陣描述二維Hilbert的層級演進關系，將坐標變換轉換為C++中的數(shù)組運算，減少了嵌套循環(huán)中迭代次數(shù)。文獻[24]研究了三維Hilbert曲線的層級演進關系，給出了三維Hilbert曲線的初始地址表與層級演進表，但是其給出的地址表形式為若干個二進制數(shù)的簡單集合，無法直接應用于編碼運算。

鄰近格元索引是空間聚類、查詢等基本空間操作的基礎[25]。計算得到鄰近格元Hilbert碼是實現(xiàn)鄰近格元索引的前提。典型的鄰近格元Hilbert碼計算基于Morton碼與Hilbert碼之間的空間變換關系實現(xiàn)[26-27]，分析其步驟可發(fā)現(xiàn)其存在兩個缺點：①需要逐層級拆解，將Hilbert碼轉換為Morton碼，耗費較長時間，單次拆解復雜度為O(n)，n為Hilbert曲線維度。②當中心格元層級為m時，需要進行2次步長為m的循環(huán)，總時間復雜度為O(2nm)，隨著格元層級的增高，效率隨之降低。文獻[28]總結了二維Hilbert曲線的層級演進關系，通過記錄中心格元在各層級中的相對位置，實現(xiàn)鄰近格元Hilbert碼計算，對比Morton碼轉換算法在效率上有較大提升，但是空間維度變化會帶來空間填充順序與演進關系的差異，導致二維Hilbert碼算法無法直接應用至更高維度曲線中。

Hilbert曲線具有自相似性[5]，其體現(xiàn)為曲線在不同層級間的形狀特征是相似的，高階曲線由多條低階曲線做空間變換后相連構成。研究高低階Hilbert曲線之間的層級演進關系將會為Hilbert碼的計算效率帶來顯著提升，并為之后的空間分析帶來新的模型與方法[15,24,29]。本文首先設計八叉樹網格的Hilbert編碼結構，然后構建Hilbert曲線層級演進模型，最后以層級演進模型為基礎，分別解決笛卡兒坐標至Hilbert碼轉換計算、鄰近格元Hilbert碼計算問題。

1 格元Hilbert碼設計

Hilbert空間填充曲線與網格單元(簡稱格元)之間的連續(xù)映射關系在數(shù)學角度上可描述為一種線段T至n維數(shù)據(jù)空間Rn上的連續(xù)映射函數(shù)，即線段T中的每一等份t與n維數(shù)據(jù)空間Rn中的每一格元Q之間的存在一一映射關系，格元Q的大小由線段的等分次數(shù)m確定。在Hilbert空間填充曲線的構造過程中，m-1階曲線(低階)至m階曲線(高階)的一次演進時，將各個維度等分為2m等份，那么n維數(shù)據(jù)空間Rn則被剖分為2mn個m層級格元Q，這些m層級格元Q與線段T的2mn等分的t之間一一對應，由同一個m-1格元剖分得到的2n個m層級格元的集合被稱作Hilbert細胞[30]。

圖1 Hilbert碼與Hilbert細胞Fig.1 Hilbert code and Hilbert cell

通過記錄1至m階曲線演進過程中，格元在每一個細胞中的填充順序，可為格元設置唯一對應的Hilbert碼，其編碼形式如圖1(a)所示。將m層級格元的Hilbert碼記作hm，格元在各層級細胞中的填充順序記作i1,i2,…,im，Hilbert碼hm滿足式(1)

hm=i1⊕i2⊕…⊕im

(1)

式中，符號⊕表示左移一位后相加。Hilbert碼hm存在以下特性：①m層級格元的Hilbert碼hm共有m位；②子格元Hilbert碼與父格元Hilbert碼只相差最后一位，子格元Hilbert碼等于父格元Hilbert碼左移一位后相加上子格元的填充順序。

2 Hilbert曲線層級演進模型

2.1 Hilbert曲線狀態(tài)向量

圖2 24種填充順序Fig.2 24 filling sequences

(2)

2.2 狀態(tài)向量層級映射關系

表1 Hilbert曲線基因列表Tab.1 Hilbert curve gene list

以基因列表為基礎，可對狀態(tài)向量在層級間的映射關系進行如下推導。

由曲線構造過程可得，一個狀態(tài)向量為sk的細胞(對應一個m層級格元)，其各個子格元(m+1層級格元)的狀態(tài)向量滿足以下層級映射公式

(3)

式中，表達式sk?G[i]表示對狀態(tài)向量sk的細胞做坐標變換G[i]后重新記錄狀態(tài)向量。

圖3 狀態(tài)向量s2的層級映射Fig.3 state vector s2’s level mapping diagram

(4)

以式(4)為基礎，通過以下步驟遞推可得其余狀態(tài)向量的層級映射：①計算狀態(tài)向量s2與狀態(tài)向量sk(k≠2)之間的坐標變換G(k,2)；②將式(4)左右兩邊各做坐標變換G(k,2)，得狀態(tài)向量sk滿足以下映射關系

(5)

計算所有狀態(tài)向量的層級映射關系后，引入演進矩陣E記錄24種低階狀態(tài)向量與高階狀態(tài)向量間的映射關系，見式(6)

(6)

狀態(tài)演進矩陣E的第u行對應映射關系為

(7)

狀態(tài)矩陣S與演進矩陣E是本文進行三維Hilbert曲線編碼層級演進模型描述的重要工具：狀態(tài)矩陣S記錄了三維Hilbert曲線的填充順序，演進矩陣E記錄了三維Hilbert曲線填充順序在高低階曲線之間的變化趨勢，兩個矩陣相互配合實現(xiàn)三維Hilbert曲線層級演進模型的形式化描述。

在文獻[22，30]中，基因列表只表示了Hilbert曲線某一種填充順序在層級間的層級演進關系，因而在運算使用中需要臨時通過多維度迭代計算，才能得出下一層級格元的填充順序。本文中，首先，求取了三維Hilbert曲線所有情況的填充順序，然后，推導所有填充順序的層級演進關系，運算使用中可直接從層級演進關系模型中提取得到當前填充順序與下一層級填充順序，因此便于設計更加簡潔直觀的算法。

3 基于層級演進模型的格元Hilbert碼計算

Hilbert曲線層級演進模型對曲線的填充順序以及其在層級間的變化關系進行描述，在此基礎上可以進行Hilbert曲線構造過程的復現(xiàn)，并設計實現(xiàn)格元Hilbert碼計算方法。

3.1 笛卡兒坐標至Hilbert碼計算

全球離散網格采用網格單元編碼代替?zhèn)鹘y(tǒng)的二維、三維笛卡兒坐標參與運算。然而現(xiàn)有的多數(shù)地理空間數(shù)據(jù)仍然采用笛卡兒坐標給出，為保證現(xiàn)有數(shù)據(jù)能夠融入全球離散網格系統(tǒng)中，需要設計相應的格元笛卡兒坐標至Hilbert碼計算的方法。

笛卡兒坐標至Hilbert碼的計算實質是對格元剖分和曲線演進的復現(xiàn)，即記錄格元在各次曲線演進過程中所處的細胞填充順序，故曲線演進過程復現(xiàn)的效率不同將導致計算效率存在差異。本文借助Hilbert曲線層級演進模型復現(xiàn)曲線演進過程，避免了各維度迭代步驟，具有直觀簡潔的特點。

以狀態(tài)矩陣S與演進矩陣E為工具，本文的格元剖分和曲線演進復現(xiàn)思路如下：從1層級格元和1階Hilbert曲線為起點，逐層級、階級向高層級、階級剖分演進，其中在各次剖分和演進過程中，先使用狀態(tài)矩陣S判斷當前格元的填充順序，后使用演進矩陣E更新下一層級細胞的狀態(tài)向量。圖4給出了計算示意圖，圖中深色格元坐標P(3,0,3)，采用狀態(tài)矩陣S與演進矩陣E進行三次復現(xiàn)后，得到Hilbert碼h3=042。計算過程的具體步驟如圖5所示。

圖4 笛卡兒坐標至Hilbert碼計算Fig.4 Cartesian coordinate to Hilbert code calculation

圖5 笛卡兒坐標至Hilbert碼計算流程Fig.5 calculation process from Cartesian coordinates to Hilbert code

文獻[22]實現(xiàn)笛卡兒坐標至Hilbert碼的轉換計算過程中，單次運算共需要對n個維度上的坐標進行迭代變換，復雜度為O(n)，算法共需要進行m次循環(huán)運算，其總時間復雜度為O(nm)。本文使用Hilbert曲線層級演進模型復現(xiàn)曲線構造過程，不存在迭代變換步驟，單次運算只需使用狀態(tài)矩陣S查詢填充順序、演進矩陣E查詢下一層級狀態(tài)向量，查詢復雜度為O(1)，全程需進行m次循環(huán)，總時間復雜度為O(m)，相比于文獻[22]迭代算法復雜度更小。

3.2 鄰近格元Hilbert碼計算

鄰近格元索引是空間聚類、鄰近分析等空間操作的基礎。本文以Hilbert曲線層級演進模型為基礎，設計鄰近格元Hilbert碼計算方法，實現(xiàn)鄰近格元索引。

在Hilbert碼包含父格元與子格元的層級關系，子格元Hilbert碼與父格元Hilbert碼只相差最后一位，子格元Hilbert碼等于父格元Hilbert碼左移一位后相加上子格元的填充順序。因此，通過Hilbert碼位運算可實現(xiàn)父子單元的查詢，位運算分為左移與右移兩種形式：①右移運算。hm?1，即m級格元Hilbert碼右移1位得到其m-1級父格元的Hilbert碼。②左移運算。hm?1，即m級格元Hilbert碼左移1位得到其m+1級子格元中第一個填充的子格元Hilbert碼，若需計算第i個填充子格元的Hilbert碼，則左移1位后將最后一位替換為i。

若通過右移運算可尋找到公共父格元，則可采用左移運算得到鄰近格元。進行位運算的同時，需要對Hilbert曲線采取相應的退化與演進操作，以確定格元的填充順序i以及父子格元的狀態(tài)向量。

s(C(v,i))=sE[k-1][i]

(8)

s(v)=sr(k)+1

(9)

式中，r(k)表示演進矩陣E第i列中值為k的元素所在行數(shù)。獲得父格元狀態(tài)向量后，先依據(jù)填充順序判斷當前格元所處的子格元編號，后判斷指定鄰近方向DirN上鄰近格元與中心格元是否同屬于一個父格元(判斷依據(jù)見表2)。若同屬于一個父格元，則查詢得到公共父格元，下一步開始演進過程。

表2 鄰近格元編號Tab.2 Neighbor grid elements number

圖6 鄰近格元Hilbert碼計算 Fig.6 Hilbert code calculation of neighbor grid elements

圖7 鄰近格元Hilbert碼計算流程Fig.7 Calculation flow of Hilbert codes with neighbor grid elements

在文獻[26]的轉換算法中，為了生成鄰近格元Hilbert碼，需要通過2m次復雜度為O(n)的拆解轉換步驟，其時間復雜度為O(2nm)，因此隨著格元層級m的增高，其計算效率將會逐漸降低。在本文計算方法中，采用Hilbert曲線層級演進模型實現(xiàn)了Hilbert曲線的退化與演進，無需將Hilbert碼轉換為其他類型的碼或坐標，算法的時間復雜度為O(2k)，k為中心格元與公共父格元之間的層級差。本文方法單次計算復雜度為O(1)，相較文獻[26]轉換算法更低，且循環(huán)次數(shù)與格元層級無關，故計算效率更高。

4 試驗與分析

為測試上述以Hilbert曲線層級演進模型為基礎的算法效率，本文采用Visual C++ 2015開發(fā)工具分別實現(xiàn)了笛卡兒坐標至Hilbert碼的迭代計算方法[22]與本文計算方法，以及Morton碼轉換鄰近格元Hilbert碼計算方法[26]與本文鄰近格元Hilbert碼計算方法。全部程序編譯為Release版本，并在計算機(CPU Intel Core i7-7700K 雙核4.2 GHz，內存64 GB，硬盤7200 r/min)上進行測試。

4.1 笛卡兒坐標至Hilbert碼計算效率測試

為測試本文基于Hilbert曲線層級演進模型的笛卡兒坐標至Hilbert碼計算效率，設計試驗進行測試：首先，在不同計算次數(shù)、相同階數(shù)情況下，測試本文算法完成計算所需時間，并與迭代算法進行對比。然后，在不同階數(shù)、相同計算次數(shù)情況下，測試本文算法完成計算所需時間，并與迭代算法進行對比。測試具體數(shù)值設置如下。

(1) 隨機生成5×106個格元P(x,y,z)，使用本文計算方法和迭代計算方法生成格元對應的5-15階Hilbert碼，記錄兩種算法計算耗時，結果如圖8(a)所示。

(2) 隨機生成{1×105，2×105，…,20×105}個格元P(x,y,z)，使用兩種計算方法生成格元對應的10階Hilbert碼，記錄兩種算法計算耗時，結果如圖8(b)所示。

圖8 測試結果Fig.8 Test results

由測試結果可知：

(1) 圖8(a)中，在計算不同階數(shù)Hilbert碼時，本文方法所需的轉換時間隨轉換階數(shù)的增加而呈上升趨勢，在階數(shù)15時，進行1次轉換所需的平均時間約為0.8 μs；對比迭代方法，在不同階數(shù)時，本文方法的計算耗時更少，效率更高，如表3所示在各組試驗中，本文算法效率提升7%至23%。

表3 不同階數(shù)5×106個格元計算測試結果對比Tab.3 Comparison of calculation results of 5×106 grid elements with different orders

(2) 圖8(b)中，在計算不同次數(shù)的Hilbert碼時，本文方法所需的轉換時間隨轉換次數(shù)的增加而呈上升趨勢，進行1次轉換所需的平均時間約為0.7 μs；對比迭代方法，本文方法的計算耗時更少，效率更高，見表4。在各組試驗中，本文算法效率提升8%至15%。

表4 不同次數(shù)測試結果對比Tab.4 Comparison of test results of different times

分析差異產生原因，本文方法以Hilbert碼的層級演進模型為基礎，采用狀態(tài)矩陣S與演進矩陣E復現(xiàn)Hilbert曲線演進過程，避免了多維度迭代的步驟，算法過程更加直接簡潔，因而獲得更高的效率。本文方法能夠以微妙級時間為代價完成1次笛卡兒坐標至Hilbert碼的轉換，對比于迭代算法，計算耗時更少，效率更高。

4.2 鄰近格元Hilbert碼計算效率測試

為測試本文基于Hilbert曲線層級演進模型的鄰近格元Hilbert，碼計算效率，設計試驗進行測試：首先，在相同層級、不同格元個數(shù)情況下，測試本文鄰近格元Hilbert碼計算方法完成計算所需時間，并與迭代算法進行對比。隨后，在相同格元個數(shù)、不同格元層級情況下，測試本文鄰近格元Hilbert碼計算方法完成計算所需時間，并與迭代算法進行對比。測試試驗具體數(shù)值設置如下。

(1) 在層級15情況下，隨機選取當前層級的{1×106,2×106,…,8×106}個格元，使用Morton碼轉換鄰近格元Hilbert碼計算方法與本文鄰近格元Hilbert碼計算方法進行某一隨機方向上鄰近格元Hilbert碼計算，記錄兩種算法計算耗時，結果如圖9(a)所示。

(2) 在層級5—15下，隨機選取當前層級中的1×106個格元，并分別采用兩種計算方法進行某一隨機方向上鄰近格元Hilbert碼計算，記錄兩種算法計算耗時，結果如圖9(b)所示。

由測試結果可知：

(1) 圖9(a)中，相同格元層級情況下，兩種方法對應曲線隨著格元個數(shù)的增加均勻上升，表明兩種方法的計算所需時間均與格元個數(shù)呈正相關；對比同一格元數(shù)下，兩種方法所需的計算時間，本文方法的計算時間均要少于轉換方法，各組測試中本文算法效率提高4.0至4.5倍，見表5。

圖9 測試結果Fig.9 Test results

表5 不同格元數(shù)測試結果對比Tab.5 Comparison of test results of different grid elements

(2) 圖9(b)中，相同格元個數(shù)情況下，轉換方法對應曲線隨層級增加呈現(xiàn)較為均勻的增長趨勢，說明其計算效率隨層級的增加而線性降低；本文方法對應曲線較為平緩，隨著層級的增加存在較為微小的波動，但計算耗時整體處于2 s以內，表明本文方法效率受層級影響較小。對比圖9(b)中同一層級下，兩種方法所需的計算時間，本文方法所需時間均小于轉換方法。本文算法效率提升倍數(shù)隨層級增加有擴大趨勢，在層級15情況下，效率提升了4.4倍，見表6。

表6 不同層級測試結果對比Tab.6 Comparison of test results at different levels

分析差異產生原因，轉換方法需要將Hilbert碼通過層級拆解轉換為Morton碼后，才能實現(xiàn)鄰近格元Hilbert碼計算，其時間復雜度與層級呈線性相關；本文計算方法則是通過層級演進模型復現(xiàn)了填充順序在層級間的變化關系，只需尋找公共父格元即可實現(xiàn)鄰近格元Hilbert碼計算，因而效率更高，且算法時間復雜度只與中心格元與公共父格元的層級差相關，計算效率不隨中心格元的層級發(fā)生顯著改變。綜上所述，在鄰近格元Hilbert碼計算效率上，本文方法的效率優(yōu)于轉換方法。

5 總結與展望

高效、簡便的網格單元編碼運算及索引是全球離散網格系統(tǒng)的核心問題。作為網格編碼研究的重要工具，Hilbert曲線網格編碼實現(xiàn)了坐標等效降維表達，但是如何直接計算Hilbert曲線不同層級之間的變換、如何將格元多維空間結構與關系轉換到一維Hilbert碼的相關運算上仍然需要深入研究。本文對八叉樹網格Hilbert曲線層級演進關系進行建模，采用狀態(tài)矩陣S與演進矩陣E復現(xiàn)三維Hilbert曲線的演進過程。以層級演進模型為理論基礎，設計笛卡兒坐標至Hilbert碼計算方法以及鄰近格元Hilbert碼計算方法，避免現(xiàn)有算法中煩瑣的迭代步驟，提高了Hilbert碼層級屬性的利用率，計算效率較高。下一步還需要研究如何在全球離散網格系統(tǒng)上運用Hilbert曲線層級演進模型這樣的工具，包括：從局部歐氏空間到全球流形的變化帶來方向、距離等度量不同，不同的全球離散網格系統(tǒng)在剖分方法上的不一致帶來空間填充曲線具體實現(xiàn)等問題。