基于KD樹的規(guī)則格網(wǎng)DEM插值技術(shù)

黃 艷，盧學(xué)良

1.西安測(cè)繪研究所，陜西 西安，710054;

2.地理信息工程國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安，710054

1 引 言

數(shù)字高程模型(digital elevation model，DEM)作為一種表達(dá)地形的數(shù)學(xué)模型，自1958年由美國麻省理工學(xué)院Miller教授提出以來，由于其表達(dá)方式直觀多樣、便于更新等優(yōu)勢(shì)，已經(jīng)在測(cè)繪導(dǎo)航、農(nóng)林規(guī)劃、地學(xué)分析等領(lǐng)域得到了廣泛的應(yīng)用和研究，成為國家經(jīng)濟(jì)建設(shè)以及國防建設(shè)的重要基礎(chǔ)數(shù)據(jù)之一[1]。DEM的生產(chǎn)途徑多種多樣，目前應(yīng)用最為廣泛、生產(chǎn)效率最高的方法一般是通過航空航天攝影測(cè)量、合成孔徑雷達(dá)干涉測(cè)量或是機(jī)載激光掃描等方式獲取地形地物的點(diǎn)云數(shù)據(jù)，然后通過內(nèi)插的方式獲得規(guī)則格網(wǎng)的DEM[2，3]。在這些方法中，插值是DEM生產(chǎn)的重要環(huán)節(jié)，常用的DEM插值算法有反距離加權(quán)插值算法、謝別德插值算法和徑向基函數(shù)插值算法等[4，5]，這些算法無不是利用待插值點(diǎn)鄰域范圍內(nèi)的已知采樣點(diǎn)通過特定的加權(quán)方式完成插值。因此，鄰域采樣點(diǎn)的獲取是完成插值的前提，而為了高效獲取待插值點(diǎn)的鄰域采樣點(diǎn)，一般要對(duì)點(diǎn)云建立空間索引，常用的空間索引方法為柵格法。本文在深入分析柵格法索引的基礎(chǔ)上，引入了KD樹來實(shí)現(xiàn)鄰域點(diǎn)的快速獲取，并通過試驗(yàn)對(duì)兩種空間索引的效果進(jìn)行了對(duì)比分析。

2 DEM插值

DEM插值是根據(jù)已知采樣點(diǎn)的高程值去估計(jì)未知插值點(diǎn)高程值的過程[6]，其主要用于缺值估計(jì)、等值線內(nèi)插和離散點(diǎn)云數(shù)據(jù)的格網(wǎng)化[7]。DEM插值可以通過已知采樣點(diǎn)來估計(jì)未知插值點(diǎn)，主要是因?yàn)榈匦尉哂锌臻g異質(zhì)性和空間相關(guān)性等基本特征，使得利用一些空間位置合理的采樣點(diǎn)(一般為鄰域采樣點(diǎn))獲得對(duì)地形表面相對(duì)精確的描述成為可能。

2.1 鄰域采樣點(diǎn)的獲取

(1)柵格法索引及其存在的問題

在對(duì)大數(shù)據(jù)量的地面點(diǎn)云進(jìn)行插值時(shí)，插值點(diǎn)鄰域采樣點(diǎn)的獲取一般不宜采用全局搜索的方法。為了提高搜索效率，首先要建立點(diǎn)云數(shù)據(jù)的空間索引，常用的空間索引方法為柵格法。

柵格法的基本思想如圖1所示，首先輸入點(diǎn)云數(shù)據(jù)，在確定點(diǎn)云數(shù)據(jù)的平面范圍后，將該平面區(qū)域按照設(shè)定的尺度柵格化，柵格單元一般為正方形，然后將各個(gè)小柵格與其內(nèi)包含的點(diǎn)關(guān)聯(lián)起來。對(duì)于待插值點(diǎn)，在確定該點(diǎn)所在的柵格后，通常該柵格及其8鄰域柵格所包含的點(diǎn)即為待插值點(diǎn)的候選鄰域采樣點(diǎn)，最后經(jīng)過距離約束或其他方式即可篩選出鄰域采樣點(diǎn)。

圖1 柵格索引建立及檢索流程

利用柵格法建立空間索引主要有兩個(gè)要點(diǎn)：一是柵格大小的確定;二是待插值點(diǎn)搜索范圍的確定，即待插值點(diǎn)所在柵格的鄰域柵格的范圍。其中柵格的大小與DEM的分辨率有關(guān)，一般是DEM分辨率的若干倍。而鄰域柵格的選擇范圍一般是其8鄰域，當(dāng)然也可以向外擴(kuò)充，原則上是柵格越大，鄰域的選擇范圍越小，柵格越小，鄰域的選擇范圍越大。

然而利用柵格法作為空間索引，可能存在某些點(diǎn)比部分鄰域采樣點(diǎn)距離待插值點(diǎn)更近的情況。如圖2所示，顯而易見的是doutside＜dinside，這樣很可能會(huì)導(dǎo)致計(jì)算出的待插值點(diǎn)的若干個(gè)最近鄰點(diǎn)不夠準(zhǔn)確(并非最鄰近點(diǎn))，進(jìn)而會(huì)影響最終的插值精度。

圖2 鄰域柵格的內(nèi)點(diǎn)與外點(diǎn)距離示意

目前有兩種辦法來解決上述情況。第一種方法是在DEM插值點(diǎn)設(shè)定某個(gè)距離閾值d，如果某些最鄰近采樣點(diǎn)與待插值點(diǎn)之間的距離大于d，那么它將會(huì)被排除。因?yàn)榫嚯x較遠(yuǎn)的采樣點(diǎn)可能并不能提高待插值點(diǎn)的插值精度，甚至?xí)兴鶕p害，基于這點(diǎn)考慮，將柵格的邊長設(shè)置為d，便可以規(guī)避上述問題，但這樣會(huì)明顯降低空間索引的效率。另一種方法是如果理想中的鄰域柵格的選擇范圍是其8鄰域，在實(shí)際處理中把選擇范圍擴(kuò)大到24鄰域，可以在很大程度上避免上述問題(并不能完全規(guī)避)，但同樣會(huì)降低空間索引的效率。

(2)KD樹索引

鑒于上述問題，本文在空間索引上使用了KD樹。KD樹(K-Dimension tree)最早由斯坦福大學(xué)的Jon Bentley教授在論文中提出[8]。它是一種應(yīng)用于高維空間的數(shù)據(jù)索引結(jié)構(gòu)，采用分而治之的思想將整個(gè)空間數(shù)據(jù)集在K維空間按照一定的規(guī)則順次劃分為若干個(gè)小空間。二維KD樹的構(gòu)建如圖3所示。首先紅色的線條將點(diǎn)集分成兩個(gè)部分，然后每個(gè)子塊被藍(lán)色線條分別分成兩個(gè)部分，接著按照綠色線條、紫色線條、淺藍(lán)色線條的順序依次將子塊劃分，直至子塊中的點(diǎn)集規(guī)模小于局部規(guī)模閾值。KD樹實(shí)際上就是決策樹，每次找區(qū)分度最大的一條軸線進(jìn)行二分，二分至葉子節(jié)點(diǎn)只包含一個(gè)點(diǎn)。在查找最鄰近點(diǎn)時(shí)，先找到葉子節(jié)點(diǎn)，然后回溯其祖先節(jié)點(diǎn)，以這個(gè)點(diǎn)的距離剪枝搜索，最終找到符合要求的若干個(gè)最鄰近點(diǎn)。

圖3 二維KD樹[9]

2.2 常用的插值算法

目前常用的插值算法有反距離加權(quán)插值算法、謝別德插值算法和徑向基函數(shù)插值算法等。

(1)反距離加權(quán)插值算法

反距離加權(quán)插值算法是根據(jù)相近相似的原理，對(duì)每個(gè)采樣點(diǎn)賦予一定的權(quán)重，權(quán)重隨著采樣點(diǎn)與插值點(diǎn)之間距離的增加而減小。任一插值點(diǎn)處的值是各采樣點(diǎn)加權(quán)之和，如式(1)所示。

其中，zp為插值點(diǎn)的高程值;λi為第i個(gè)點(diǎn)的權(quán)重;di為第i個(gè)采樣點(diǎn)到插值點(diǎn)的距離;d-u為距離衰減函數(shù)，冪指數(shù)u具有隨著距離的增加而減小其他位置影響的作用，通常取值為1或2。

(2)謝別德插值算法

謝別德插值算法本質(zhì)上是一種標(biāo)準(zhǔn)的導(dǎo)數(shù)距離加權(quán)過程，權(quán)函數(shù)如式(2)所示。

其中，wi為權(quán)重;di為第i點(diǎn)距待插值點(diǎn)的距離;r為調(diào)整距離。

(3)徑向基函數(shù)插值算法

徑向基函數(shù)插值算法是一系列用于精確插值算子的統(tǒng)稱，其插值原理是任何一個(gè)表面都可以用多個(gè)曲面的線性組合去逼近。通常情況下，徑向基函數(shù)插值算法可以表述為兩部分之和，如式(3)所示。

其中，zp為插值點(diǎn)的高程值;λi為i第個(gè)點(diǎn)的權(quán)重;di為第i個(gè)采樣點(diǎn)到插值點(diǎn)的距離;φ(di)為徑向基函數(shù);fi為“趨勢(shì)”函數(shù)，是次數(shù)小于m的基本多項(xiàng)式函數(shù)。

3 試驗(yàn)與分析

3.1 試驗(yàn)方案

試驗(yàn)一：通過對(duì)前述兩種空間索引的對(duì)比分析，測(cè)試空間索引效果。試驗(yàn)環(huán)境為Intel(R)Xeon(R)CPU E5-2630 v3(2.4GHz)，32G內(nèi)存，操作系統(tǒng)是Window7 64位。試驗(yàn)數(shù)據(jù)采用某區(qū)域5m分辨率的衛(wèi)星影像通過密集匹配并濾波后得到的點(diǎn)云數(shù)據(jù)。點(diǎn)云數(shù)據(jù)有A、B兩塊，A區(qū)域的數(shù)據(jù)量為：356MB(文本文件)，含有10989051個(gè)點(diǎn);B區(qū)域的數(shù)據(jù)量為：1.14GB(文本文件)，含有36551055個(gè)點(diǎn)。DEM的采樣分辨率設(shè)置為15m，最鄰近點(diǎn)的檢索數(shù)量設(shè)置為4個(gè)，距離閾值d設(shè)置為60m(即4倍的DEM采樣分辨率)。主要統(tǒng)計(jì)的指標(biāo)有最鄰近點(diǎn)查詢的正確率、空間索引的構(gòu)建時(shí)間以及DEM插值計(jì)算時(shí)間(鄰域點(diǎn)的搜索與插值)。對(duì)于柵格法由于不同參數(shù)的設(shè)定會(huì)產(chǎn)生不同的試驗(yàn)結(jié)果，這里設(shè)定k為柵格大小與DEM分辨率的比率，n為待插值點(diǎn)所在柵格的鄰域柵格的范圍(n鄰域)。

試驗(yàn)二：利用試驗(yàn)一中表現(xiàn)效果較好的空間索引，通過反距離加權(quán)插值算法完成A區(qū)域的DEM的采樣，并通過目視評(píng)估采樣效果。

3.2 試驗(yàn)分析與結(jié)論

試驗(yàn)一的結(jié)果見表1。

表1 不同空間索引的統(tǒng)計(jì)指標(biāo)

通過表中的指標(biāo)統(tǒng)計(jì)，可以看出，點(diǎn)云數(shù)量和插值點(diǎn)數(shù)的比例并非理論上的9∶1(DEM分辨率和影像分辨率的比例為3∶1)，這是因?yàn)閷?shí)際插值的范圍并非恰好是點(diǎn)云數(shù)據(jù)的平面分布范圍，而是矩形包圍盒的范圍，因此，插值點(diǎn)的數(shù)量會(huì)比理論上要多。對(duì)于柵格法，空間索引的構(gòu)建時(shí)間與點(diǎn)云數(shù)量成正比，與柵格的大小沒有關(guān)系，而插值的計(jì)算時(shí)間與檢索的鄰域范圍成正相關(guān)關(guān)系，與插值點(diǎn)數(shù)成正比。當(dāng)k=2、n=8時(shí)，正如前述的分析的那樣，并不能保證最鄰近點(diǎn)檢索的正確性，而對(duì)于k=2、n=24和k=4、n=8的情況，試驗(yàn)結(jié)果與理論相符(與距離閾值d的設(shè)定有關(guān))，能保證最鄰近點(diǎn)檢索的完全正確。對(duì)于KD樹，空間索引構(gòu)建時(shí)間與點(diǎn)云的數(shù)量成正比，插值計(jì)算時(shí)間與插值點(diǎn)數(shù)成正比，并且KD樹能保證檢索的最鄰近點(diǎn)完全正確。對(duì)比兩種空間索引，在保證正確率的前提下，可以得出的結(jié)論是KD的效率要優(yōu)于柵格法。

根據(jù)上述結(jié)論這里選取KD樹作為空間索引進(jìn)行試驗(yàn)二。試驗(yàn)二的結(jié)果如圖4所示。

圖4 DEM采樣效果

圖4中左上角為原始點(diǎn)云，左下角為采樣后的DEM，右側(cè)分別為原始點(diǎn)云和采樣后DEM的局部放大圖，從圖中可以看出基于KD樹的DEM插值的效果較為理想。

4 結(jié) 語

本文對(duì)DEM插值中常用的空間索引進(jìn)行了詳細(xì)的介紹和分析，指出了其中存在的問題，并在此基礎(chǔ)上引入了KD樹。試驗(yàn)證明，KD樹在效率方面具有優(yōu)勢(shì)，并且基于KD樹的DEM插值的效果很好。當(dāng)然柵格法的特點(diǎn)是實(shí)現(xiàn)簡(jiǎn)單、靈活，對(duì)于不同的數(shù)據(jù)特點(diǎn)、參數(shù)設(shè)置和成果要求，柵格法同樣可以加以運(yùn)用。

[1]張亞南.DEM分辨率確定與尺度轉(zhuǎn)換方法研究[D].南京：南京師范大學(xué)，2014.

[2]張海霞，陳向陽，趙文普等.高分辨率遙感影像自動(dòng)生成DEM的研究[J].測(cè)繪與空間地理信息，2015(7)：103-105.

[3]湯國安.我國數(shù)字高程模型與數(shù)字地形分析研究進(jìn)展[J].地理學(xué)報(bào)，2014，69(9)：1305-1325.

[4]張錦明.DEM插值算法適應(yīng)性研究[D].鄭州：信息工程大學(xué)，2012.

[5]王耀革.DEM建模與不確定性分析[D].鄭州：信息工程大學(xué)，2009.

[6]盧華興，劉學(xué)軍，湯國安.DEM誤差模型研究[C].鄭州：第三屆地理信息系統(tǒng)全國博士生學(xué)術(shù)論壇，2008.

[7]李新，程國棟，盧玲.空間內(nèi)插方法比較[J].地球科學(xué)進(jìn)展，2000，15(3)：260-265.

[8]Bentley JL.Multidimensional Binary Search Trees Used for Associative Searching[J].Communications of the ACM，1975，18(9)：509-517.

[9]Muja M，Lowe D G.Scalable Nearest Neighbor Algorithms for High Dimensional Data[J].IEEE Transactions on Pattern Analysis&Machine Intelligence，2014，36(11)：2227-2240.