丁加成，趙學勝

(中國礦業(yè)大學(北京)地球科學與測繪工程學院，北京 100083)

0 引言

球面退化四叉樹格網(Degenerate Quadtree Grid,DQG)[1]是一種建立在經緯度格網基礎上的變間隔經緯度格網模型，具有徑向對稱性、方向一致性及平移相和性等優(yōu)點[2]，其四邊形網格的投影為矩形[3]，適用于數(shù)字地形的構建和遙感影像的管理等[4,5]，是全球地表建模及可視化分析[6]的基礎框架之一。格網利用編碼表達空間對象，而現(xiàn)有空間數(shù)據(jù)多以經緯度的形式存在[7]。因此，編碼與經緯度間的轉換效率是制約格網應用的重要因素[8]。隨著空間信息獲取技術的進步和更新頻率的加快，時空數(shù)據(jù)量急劇增加[8,9]，數(shù)據(jù)類型也日益豐富[10]，對格網編碼與經緯度間的轉換效率提出了更高要求。

國內外學者提出了諸多格網編碼與經緯度轉換方法。例如：ZOT (Zenithal Ortho Triangular)投影法[11]按曼哈頓距離進行轉換，效率較高，但編碼失去方向性，不便于鄰近搜索；等三角投影法(Equal Triangles Projection,ETP)[12]計算精度較高，但涉及投影操作及冪指數(shù)運算，效率較低；行列逼近法[7]在保證算法效率的同時，生成的編碼具有方向性，但存在半個格網誤差；三向互化算法[13]有效降低了轉換誤碼率，其轉換速度介于行列逼近法和ZOT算法之間；王謙等[8]對行列逼近法和三向互化算法進行改進，提高了編碼轉換精度；張玉梅等[14]提出的球面菱形格網與經緯度轉換算法不涉及投影操作，計算速度較快。上述算法在編碼與經緯度間轉換時存在遞歸過程，轉換效率有待提升。針對DQG格網編碼，崔馬軍等[15]最早利用DQG和經緯度格網的對應關系提出相應的轉換算法，實現(xiàn)簡單，但其采用四進制編碼，轉換效率需進一步提升；余接情等[16]對DQG進行三維擴展，提出一種球體退化八叉樹立體格網SDOG，采用二進制方式存儲編碼，并設計對應的單層次和多層次格網編碼方式，但需進行多次浮點運算，使用逐位處理進行編碼/解碼，特別是處理海量數(shù)據(jù)時，其效率仍需進一步提升。

為適應海量格網數(shù)據(jù)的實時計算與分析需求，本文在現(xiàn)有算法基礎上進行改進：用二進制代替四進制編碼，并推導出格元行號和經度差間關系，優(yōu)化格元列號計算方式，以期提高編(解)碼效率，為海量時空大數(shù)據(jù)的高效轉換和實時分析提供可能。

1 DQG格網剖分與八分體編號

隨著緯度升高，DQG格元的經度差自適應變大，在極點周圍退化為三角形。為分析DQG格元的分布特征，構建DQG格元的行號、列號、格元經緯差以及剖分層次間關系，本文引出梯形區(qū)域的概念。

1.1 梯形區(qū)域的定義

余接情等[16]定義緯域為同一層域中經度差相同的格元集合。由于DQG中沒有層域的概念，緯度的計算過程不涉及層域，本文定義八分體內經度差和緯度差均相同的格元集合為梯形區(qū)域。對同一剖分層次，格元緯度差相同，因此，梯形區(qū)域的定義可簡化為八分體內經度差相同的格元集合，八分體由多個梯形區(qū)域組成。如圖1所示，將梯形區(qū)域按緯度由高到低進行編號，最大編號kmax=剖分層次level，包含極點的退化三角形編號為0。圖中右側數(shù)字表示各梯形區(qū)域內格元的行數(shù)，除0號區(qū)域外，梯形區(qū)域內格元的列數(shù)是行數(shù)的兩倍。

圖1 梯形區(qū)域的分布和編號Fig.1 Distribution and serial number of trapezoid areas

1.2 梯形區(qū)域編號與行列號間關系

表1 梯形區(qū)域的行數(shù)和列數(shù)Table 1 Number of rows and columns of trapezoid areas

(1)

(2)

(3)

(4)

(5)

將式(3)和式(4)代入式(5)，當k=0時，0≤i≤0，即i=0。當k>0時有以下關系成立：

(6)

對不等式取對數(shù)可得：

k-1

(7)

即：

(8)

由式(8)可知，格元所在的梯形區(qū)域編號由格元行號決定，與輸入點緯度和剖分層次無關，反映了DQG格元的分布規(guī)律。將i=0代入式(8)，得k=0，即式(8)在i=0時仍然適用。

1.3 格元經緯差計算

八分體的經度范圍為[0,π/2]，由行內格元平分，結合式(2)可計算i行格元經度差：

(9)

將式(8)帶入，可得

ΔLi=π/2log2(i+1)+1

(10)

2 DQG編碼及轉化方法

2.1 編碼存儲方式

本文以自定義結構體存儲格網編碼，用C++語言描述(圖2)。其中，整型成員m_morton存放格元行列號對應的二進制Morton碼，無符號字符成員m_oct_level的高三位存放八分體編號，低五位存放剖分層次，每個編碼占據(jù)9 B存儲空間。編碼能表達的最大剖分層次為32，對應格網單元最大邊長小于4 mm[2]，能滿足多數(shù)空間分析和表達的精度要求。

圖2 格網編碼結構體定義(C++描述)Fig.2 Grid code structure definition (C++)

2.2 Morton編碼與解碼

現(xiàn)有的DQG格網編碼/解碼方式多采用逐位處理實現(xiàn)[15-17]，效率不高。使用查找表進行編碼和解碼時，一次可處理多個二進制位，計算效率較高，故本文利用查找表[18]，結合位運算實現(xiàn)Morton碼的編碼和解碼。二維情況下，可建立大小為2m的查找表(m為映射位數(shù))。當m較小時，需進行映射的次數(shù)較多，效率較低；當m較大時，需建立的查找表較大，占用存儲空間較大。綜合考慮，本文取m=8。

編碼過程是將各維編號的二進制位按m位為一組分為若干組，每組利用查找表進行編碼映射，并將結果進行拼接，然后將各維編號經過映射和拼接后的二進制位進行移位，并進行按位或運算；解碼過程是對Morton碼的二進制位按m位為一組分為若干組，從每組二進制位中分解出各維的二進制位并進行拼接，即得到各維的二進制表示。圖3和圖4展示了映射位數(shù)m=4，對行號i=26和列號j=15進行編碼和解碼的過程。

圖3 查找表編碼過程Fig.3 Morton encoding using look up table

圖4 查找表解碼過程Fig.4 Morton decoding using look up table

此外，還可以使用位操作指令集(BMI2)中的并行位存儲指令pdep和并行位提取指令pext實現(xiàn)Morton編碼和解碼。經測試，其計算效率和查找表法基本相同。

2.3 經緯度坐標與DQG編碼間轉換算法

在DQG格元編碼和解碼過程中，需要計算格元緯度差、八分體編號及進行八分體內經緯度和全球經緯度間轉換，計算過程見文獻[15,16]。

2.3.1 經緯度坐標向DQG編碼轉換算法 輸入數(shù)據(jù)為經緯度坐標(L,B)、剖分層次level，輸出數(shù)據(jù)為經緯度對應的DQG格元編碼。1)由經緯度計算其所在八分體編號oct；2)由oct計算八分體內經緯度(L′,B′)；3)利用剖分層次level計算格元緯度差ΔB；4)由ΔB按式(11)計算格元行號i；5)按式(10)計算i行格元經度差ΔL；6)按式(12)計算格元列號j；7)將i和j的二進制交錯成Morton編碼，存放在DQG格網編碼成員m_morton中，設置成員m_oct_level的高三位為八分體編號，低五位為剖分層次。

(11)

(12)

2.3.2 DQG編碼向經緯度坐標轉換算法 輸入數(shù)據(jù)為DQG格元編碼，輸出數(shù)據(jù)為編碼對應格元的中心點經緯度。1)提取編碼成員m_oct_level的高三位為八分體編號oct，低五位為剖分層次level；2)對編碼成員m_morton進行解碼，得到行號i和列號j；3)計算格元緯度差ΔB，并按式(13)計算格元中心點的八分體內緯度Bc；4)由式(10)計算i行格元經度差ΔL，并按式(14)計算格元中心點的八分體內經度Lc；5)由oct和編碼格元中心點的八分體內經緯度(Lc,Bc)，計算對應的全球經緯度。

Bc=π/2-(i+0.5)ΔB

(13)

Lc=(j+0.5)ΔL

(14)

3 實驗分析與效率對比

實驗運行環(huán)境為虛擬機，硬件配置為RAM 3 GB，CPU AMD R5 4600H，核心數(shù)4，磁盤空間40 GB，操作系統(tǒng)為Ubuntu 20.04 LTS，編譯器版本GCC 9.3.0。實驗隨機生成分布于全球的100萬個經緯度坐標，分別測試本文方法與四進制DQG[15]、單層二維SDZ[16]編碼和解碼的效率(圖5)，并計算不同剖分層次下各轉換算法的平均轉換時間(表2)。

圖5 編碼轉換速度對比Fig.5 Comparison of code conversion speed

表2 編碼轉換算法耗時Table 2 Time consumption of code conversion algorithm

由表2可以看出，四進制DQG[15]編碼和解碼耗時分別是本文方法的23.38倍和16.93倍，而單層二維SDZ[16]編碼和解碼耗時分別是本文方法的6.28倍和2.89倍。若開啟編譯器最高優(yōu)化級別編譯程序，使用本文方法對100萬個隨機生成的地理坐標進行編碼耗時17 618 μs，將編碼轉換為地理坐標耗時10 358 μs。在虛擬機的配置環(huán)境下，每秒可實現(xiàn)約5 676萬個經緯度坐標的編碼及約7 364萬個編碼的解碼，基本能滿足海量地理坐標與格網編碼間的高效轉換。

4 結論

格網編碼與經緯度間轉換效率是影響海量空間數(shù)據(jù)集成和分析的重要因素。本文通過分析DQG格網剖分的特征，推導出格元經度差和行號間的數(shù)學關系；使用查找表法進行二進制DQG編碼與經緯度的轉換，提高了計算效率。與現(xiàn)有轉換算法不同，本文算法可根據(jù)格元所在行號直接計算格元經度差。實驗表明：本文算法的編碼和解碼的平均速度是四進制DQG算法的20.15倍，是單層二維SDZ算法的4.58倍；單機情況下可進行每秒千萬級的格網編碼與解碼，為海量格網數(shù)據(jù)的實時分析奠定基礎。但本文方法也存在一定不足，如編碼存儲的格元剖分層次有一定限制，轉換過程中的浮點數(shù)運算會降低計算效率等。下一步研究考慮使用字節(jié)數(shù)組存儲編碼以提高編碼的最大剖分層次，并使用異構計算進一步提高格元編碼與經緯度的轉換效率。