一種改進(jìn)的DQG格網(wǎng)編碼與經(jīng)緯度高效轉(zhuǎn)換算法

丁加成，趙學(xué)勝

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

0 引言

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

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

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

1 DQG格網(wǎng)剖分與八分體編號(hào)

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

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

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

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

1.2 梯形區(qū)域編號(hào)與行列號(hào)間關(guān)系

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

(1)

(2)

(3)

(4)

(5)

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

(6)

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

k-1

(7)

即：

(8)

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

1.3 格元經(jīng)緯差計(jì)算

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

(9)

將式(8)帶入，可得

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

(10)

2 DQG編碼及轉(zhuǎn)化方法

2.1 編碼存儲(chǔ)方式

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

圖2 格網(wǎng)編碼結(jié)構(gòu)體定義(C++描述)Fig.2 Grid code structure definition (C++)

2.2 Morton編碼與解碼

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

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

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

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

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

2.3 經(jīng)緯度坐標(biāo)與DQG編碼間轉(zhuǎn)換算法

在DQG格元編碼和解碼過程中，需要計(jì)算格元緯度差、八分體編號(hào)及進(jìn)行八分體內(nèi)經(jīng)緯度和全球經(jīng)緯度間轉(zhuǎn)換，計(jì)算過程見文獻(xiàn)[15,16]。

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

(11)

(12)

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

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

(13)

Lc=(j+0.5)ΔL

(14)

3 實(shí)驗(yàn)分析與效率對(duì)比

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

圖5 編碼轉(zhuǎn)換速度對(duì)比Fig.5 Comparison of code conversion speed

表2 編碼轉(zhuǎn)換算法耗時(shí)Table 2 Time consumption of code conversion algorithm

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

4 結(jié)論

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