三維Hilbert曲線層級(jí)演進(jìn)模型與編碼計(jì)算

吳宇豪，曹雪峰，余岸竹，孫萬(wàn)忠

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

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

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

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

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

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

1 格元Hilbert碼設(shè)計(jì)

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

圖1 Hilbert碼與Hilbert細(xì)胞Fig.1 Hilbert code and Hilbert cell

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

hm=i1⊕i2⊕…⊕im

(1)

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

2 Hilbert曲線層級(jí)演進(jìn)模型

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

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

(2)

2.2 狀態(tài)向量層級(jí)映射關(guān)系

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

以基因列表為基礎(chǔ)，可對(duì)狀態(tài)向量在層級(jí)間的映射關(guān)系進(jìn)行如下推導(dǎo)。

由曲線構(gòu)造過(guò)程可得，一個(gè)狀態(tài)向量為sk的細(xì)胞(對(duì)應(yīng)一個(gè)m層級(jí)格元)，其各個(gè)子格元(m+1層級(jí)格元)的狀態(tài)向量滿足以下層級(jí)映射公式

(3)

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

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

(4)

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

(5)

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

(6)

狀態(tài)演進(jìn)矩陣E的第u行對(duì)應(yīng)映射關(guān)系為

(7)

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

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

3 基于層級(jí)演進(jìn)模型的格元Hilbert碼計(jì)算

Hilbert曲線層級(jí)演進(jìn)模型對(duì)曲線的填充順序以及其在層級(jí)間的變化關(guān)系進(jìn)行描述，在此基礎(chǔ)上可以進(jìn)行Hilbert曲線構(gòu)造過(guò)程的復(fù)現(xiàn)，并設(shè)計(jì)實(shí)現(xiàn)格元Hilbert碼計(jì)算方法。

3.1 笛卡兒坐標(biāo)至Hilbert碼計(jì)算

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

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

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

圖4 笛卡兒坐標(biāo)至Hilbert碼計(jì)算Fig.4 Cartesian coordinate to Hilbert code calculation

圖5 笛卡兒坐標(biāo)至Hilbert碼計(jì)算流程Fig.5 calculation process from Cartesian coordinates to Hilbert code

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

3.2 鄰近格元Hilbert碼計(jì)算

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

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

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

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

(8)

s(v)=sr(k)+1

(9)

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

表2 鄰近格元編號(hào)Tab.2 Neighbor grid elements number

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

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

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

4 試驗(yàn)與分析

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

4.1 笛卡兒坐標(biāo)至Hilbert碼計(jì)算效率測(cè)試

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

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

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

圖8 測(cè)試結(jié)果Fig.8 Test results

由測(cè)試結(jié)果可知：

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

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

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

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

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

4.2 鄰近格元Hilbert碼計(jì)算效率測(cè)試

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

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

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

由測(cè)試結(jié)果可知：

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

圖9 測(cè)試結(jié)果Fig.9 Test results

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

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

表6 不同層級(jí)測(cè)試結(jié)果對(duì)比Tab.6 Comparison of test results at different levels

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

5 總結(jié)與展望

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