李曉紅

(山西省基礎(chǔ)地理信息院，山西太原030001)

一、引 言

機(jī)載激光雷達(dá)(light detection and ranging，LiDAR)是一種快速獲取高精度地面和地物三維信息的新技術(shù)。隨著機(jī)載激光雷達(dá)技術(shù)的高速發(fā)展，獲取的點(diǎn)云數(shù)據(jù)量越來越大。如何處理海量的點(diǎn)云數(shù)據(jù)，并從中提取有用的地形信息和地物信息，是當(dāng)前急需解決的問題。

DEM有多種表現(xiàn)形式，其中TIN和規(guī)則格網(wǎng)DEM是兩種重要的表達(dá)方式。從表達(dá)地形信息的角度而言，TIN模型的優(yōu)點(diǎn)是它能以不同層次的分辨率來描述地形表面。規(guī)則格網(wǎng)則具有存儲(chǔ)量小，數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)單，便于存儲(chǔ)、管理及分析計(jì)算等優(yōu)點(diǎn)，因此該數(shù)字地形模型得到了廣泛的應(yīng)用。規(guī)則格網(wǎng)文件結(jié)構(gòu)相對(duì)簡(jiǎn)單，存儲(chǔ)了格網(wǎng)點(diǎn)坐標(biāo)的最大值、最小值、格網(wǎng)間距、格網(wǎng)點(diǎn)數(shù)及高程值。規(guī)則格網(wǎng)可以直接從離散采樣點(diǎn)據(jù)插值得到，也可以通過不規(guī)則三角網(wǎng)(TIN)數(shù)據(jù)內(nèi)插得到。本文提出一種將LiDAR點(diǎn)云數(shù)據(jù)構(gòu)成TIN后快速內(nèi)插生成規(guī)則格網(wǎng)DEM的算法。

二、DEM的生成方法

DEM可以通過攝影測(cè)量的方式生成[1]，但是通過先獲取離散的高程點(diǎn)，再對(duì)均勻分布的格網(wǎng)點(diǎn)高程內(nèi)插生成DEM的方式更為普遍。

插值可以分為整體插值、分塊插值和逐點(diǎn)插值。整體內(nèi)插是將整個(gè)地區(qū)采用高次多項(xiàng)式來表示。由于實(shí)際地形是復(fù)雜變化多樣的，單用一個(gè)多項(xiàng)式難以表述實(shí)際地形，因此DEM內(nèi)插中一般不采用整體函數(shù)內(nèi)插。局部分塊內(nèi)插是將地形按照一定的規(guī)則進(jìn)行分塊，分塊的大小可以根據(jù)地形的復(fù)雜程度來確定。對(duì)每一塊區(qū)域采用單獨(dú)的函數(shù)模型進(jìn)行擬合內(nèi)插，不同的分塊單元可用不同的內(nèi)插函數(shù)，常用的內(nèi)插數(shù)函數(shù)有線性內(nèi)插、雙線性內(nèi)插、多項(xiàng)式內(nèi)插、樣條函數(shù)、多層曲面疊加法等。逐點(diǎn)內(nèi)插的本質(zhì)也屬于局部?jī)?nèi)插，但它與局部分塊內(nèi)插也有所不同，它是以內(nèi)插點(diǎn)為中心，確定一個(gè)指定大小的范圍，并通過這個(gè)范圍內(nèi)的點(diǎn)來選定合適的數(shù)學(xué)模型去擬合計(jì)算內(nèi)插點(diǎn)的高程。因此，逐點(diǎn)內(nèi)插的每個(gè)內(nèi)插點(diǎn)都隨著點(diǎn)的位置而變化，每個(gè)內(nèi)插的數(shù)學(xué)模型也不一樣，通常逐點(diǎn)內(nèi)插比局部分塊內(nèi)插具有更高的精度。

綜合考慮機(jī)載LiDAR點(diǎn)云濾波后數(shù)據(jù)特點(diǎn)、生成效率和精度，本文采用逐點(diǎn)內(nèi)插方法生成DEM。差值算法采用距離倒數(shù)加權(quán)插值法，它的基本原理是：對(duì)于每一個(gè)插值點(diǎn)，可以設(shè)置一個(gè)距該點(diǎn)最近的指定半徑的范圍，落在該范圍內(nèi)的點(diǎn)都被定權(quán)參與該點(diǎn)插值的計(jì)算，具體的定權(quán)方式是與距離成反比，距離插值點(diǎn)越遠(yuǎn)的點(diǎn)賦予的權(quán)值越小，反之距離越近的點(diǎn)賦予的權(quán)值越大，如圖1所示。利用TIN中Delaunay可以快速確定一個(gè)點(diǎn)一定鄰域范圍內(nèi)的點(diǎn)，以提高生成DEM效率。

三、算法流程

LiDAR點(diǎn)云生成DEM算法主要分為3個(gè)步驟：

1．將原始點(diǎn)云文件轉(zhuǎn)換為點(diǎn)云流文件

原始點(diǎn)云文件不能滿足流算法的要求，需要將原始點(diǎn)云文件轉(zhuǎn)換為點(diǎn)云流文件[2]?？筛鶕?jù)點(diǎn)云數(shù)據(jù)內(nèi)在的空間拓?fù)湫畔?，將點(diǎn)云數(shù)據(jù)平均分塊、排序插入標(biāo)識(shí)后重新寫入點(diǎn)云流文件，以便流處理算法能夠應(yīng)用。

圖1 逐點(diǎn)內(nèi)插法生成DEM

點(diǎn)云流文件的建立主要有3個(gè)步驟：

1)遍歷原始點(diǎn)云文件，計(jì)算外包矩形，將外包矩形劃分為2k×2k個(gè)分塊矩形，構(gòu)建四叉樹。

2)遍歷原始點(diǎn)云文件，計(jì)算每個(gè)分塊中包含點(diǎn)云的數(shù)目。

3)再次遍歷原始點(diǎn)云文件，對(duì)點(diǎn)云重新排序，并在每個(gè)區(qū)塊的最后一點(diǎn)插入結(jié)束標(biāo)識(shí)。

2．構(gòu)建不規(guī)則三角網(wǎng)

采用改進(jìn)后的逐點(diǎn)插入法[2]對(duì)點(diǎn)云流文件進(jìn)行Delaunay三角剖分。從點(diǎn)云流文件中讀取一個(gè)點(diǎn)時(shí)，采用逐點(diǎn)插入法進(jìn)行Delaunay三角剖分，并更新當(dāng)前的三角網(wǎng)鏈表(用于存儲(chǔ)當(dāng)前已經(jīng)生成的三角網(wǎng))。當(dāng)從點(diǎn)云文件中讀取到結(jié)束標(biāo)識(shí)時(shí)，表明當(dāng)前分塊的點(diǎn)云已經(jīng)遍歷完畢，將四叉樹(將所有分塊作為葉節(jié)點(diǎn)構(gòu)建)中的該分塊標(biāo)記為遍歷完畢。判斷當(dāng)前三角網(wǎng)鏈表中的三角形是否處于固定狀態(tài)(當(dāng)遍歷一個(gè)分塊后，如果三角形的外包圓沒有與未被遍歷的分塊相交時(shí)，認(rèn)為該三角形處于固定狀態(tài))，如果處于固定狀態(tài)，則輸出到三角網(wǎng)文件，釋放內(nèi)存;如果處于活動(dòng)狀態(tài)(如果三角形的外包圓與未被遍歷的分塊相交，則認(rèn)為三角形處于活動(dòng)狀態(tài))，則繼續(xù)保留在內(nèi)存中，等待下一分塊遍歷完畢后繼續(xù)判斷，直到所有分塊讀取完畢，完成三角剖分。

3．插值生成規(guī)則格網(wǎng)DEM

與傳統(tǒng)的DEM插值方法不同，本文所采用的插值方法不是依次尋找規(guī)則格網(wǎng)點(diǎn)的領(lǐng)域點(diǎn)進(jìn)行內(nèi)插，而是通過遍歷三角網(wǎng)文件中的三角形，判斷三角形是否包含規(guī)則格網(wǎng)點(diǎn)。如果包含規(guī)則格網(wǎng)點(diǎn)，則利用該三角形的3個(gè)頂點(diǎn)采用距離倒數(shù)加權(quán)插值法計(jì)算格網(wǎng)點(diǎn)高程，并將該格網(wǎng)點(diǎn)的行、列號(hào)及高程寫入臨時(shí)文件，釋放內(nèi)存。當(dāng)三角網(wǎng)文件遍歷完畢后，對(duì)生成的所有格網(wǎng)點(diǎn)臨時(shí)文件重新組織排序，寫入到統(tǒng)一的DEM格式文件中。這樣做的優(yōu)勢(shì)是內(nèi)存中既不保留三角網(wǎng)數(shù)據(jù)，也不保留生成的格網(wǎng)點(diǎn)數(shù)據(jù)，大大降低了程序內(nèi)存的使用量。算法流程如圖2所示。

圖2 算法流程圖

四、試驗(yàn)結(jié)果與分析

本文采用了大小為100 MB、包含300萬個(gè)點(diǎn)的LiDAR數(shù)據(jù)進(jìn)行試驗(yàn)，原始點(diǎn)云數(shù)據(jù)如圖3所示。

圖3 原始點(diǎn)云數(shù)據(jù)

首先將原始點(diǎn)云文件轉(zhuǎn)換為點(diǎn)云流文件。從原始點(diǎn)云文件轉(zhuǎn)換為點(diǎn)云流文件共耗時(shí)34 s，內(nèi)存最大使用量2．5MB。

對(duì)點(diǎn)云流文件進(jìn)行Delaunay三角剖分生成TIN，耗時(shí)9 s，內(nèi)存最大使用量2．6MB。生成的TIN如圖4所示。由于基于流的算法讀入數(shù)據(jù)進(jìn)行處理后馬上釋放內(nèi)存，只保留很少一部分未固定的數(shù)據(jù)在內(nèi)存中，所以內(nèi)存的使用量非常小。

DEM生成參數(shù)進(jìn)行設(shè)置如下：

1)采樣間隔設(shè)置：設(shè)置采樣間隔為5．0m。

2)TIN邊長(zhǎng)閾值：設(shè)置為50m(邊長(zhǎng)過長(zhǎng)的三角形內(nèi)插獲得的高程值誤差很大)。

3)BIL文件設(shè)置：高程所放比例為1．0。

圖4 三角剖分后三角網(wǎng)數(shù)據(jù)

由TIN文件內(nèi)插生成DEM共耗時(shí)3 s，內(nèi)存最大使用量2．3 MB。DEM的內(nèi)插同樣采用基于流的算法，判斷完輸入三角形后立即釋放內(nèi)存，顯著提高了算法效率，并減少了內(nèi)存使用量。生成的規(guī)則格網(wǎng)DEM如圖5所示。

圖5 DEM生成渲染圖

試驗(yàn)結(jié)果表明，采用本文算法生成規(guī)則格網(wǎng)DEM的效率非常高，內(nèi)存使用量遠(yuǎn)遠(yuǎn)低于點(diǎn)云文件本身的大小，驗(yàn)證了流算法用于海量點(diǎn)云數(shù)據(jù)處理的有效性。同時(shí)，本文算法生成的DEM準(zhǔn)確性非常高，并且能在多種復(fù)雜的地形中保持優(yōu)良的性能。

五、結(jié)束語

本文提出了一種LiDAR點(diǎn)云快速生成DEM的算法，首先將海量點(diǎn)云數(shù)據(jù)分塊，對(duì)每個(gè)分塊分別采用改進(jìn)的增量式Delaunay三角剖分算法構(gòu)建不規(guī)則三角網(wǎng);然后遍歷三角網(wǎng)文件，找到被三角形包含的格網(wǎng)點(diǎn)，采用距離倒數(shù)加權(quán)插值法計(jì)算格網(wǎng)點(diǎn)高程，存入臨時(shí)文件;最后，將所有臨時(shí)文件排序生成規(guī)則格網(wǎng)DEM。采用基于流的處理算法，實(shí)現(xiàn)了點(diǎn)云分塊處理，解決了普通計(jì)算機(jī)處理大數(shù)據(jù)量點(diǎn)云內(nèi)存不足的問題。試驗(yàn)表明，該算法在處理大數(shù)據(jù)量點(diǎn)云時(shí)，效率高，內(nèi)存占用率低，生成規(guī)則格網(wǎng)DEM的準(zhǔn)確性也非常高，能夠應(yīng)用于實(shí)際生產(chǎn)應(yīng)用中。

[1]MAUNE D F．Digital Elevation Model Technologies and Applications：The DEM Users Manual[M]．[S．l．]：ASPRS，2007．

[2]ISENBURG M，LIU Y，SHEWCHUK J，et al．Streaming Computation of Delaunay Triangulations[C]∥Proceedings of ACM SIGGRAPH 2006．New York：ACM，2006：1049-1056．

[3]ISENBURG M，LINDSTROM P．Streaming Meshes[C]∥Proceedings of IEEE Visualization．Minneapolis：IEEE Computer Society，2005：231-238．

[4]ISENBURGM，LIU Y，SHEWCHUK J，et al．Generating Raster DEM from Mass Points via TIN Streaming[C]∥Proceedings of GIScience’06．Berlin：[s．n．]，2006：186-198．

[5]ISENBURG M，LINDSTROM P，GUMHOLD S，et al．Large Mesh Simplification Using Processing Sequences[C]∥ Proceedings of the 14th IEEE Visualization．Washington DC：[s．n．]，2003．

[6]AGARWAL P K，ARGE L，DANNER A．From LiDAR to Grid DEM：A Scalable Approach[C]∥Proceedings of International Symposium on Spatial Data Handling．[S．l．]：[s．n．]，2006．

[7]CIGNONIP，MONTANIC，ROCCHINIC，et al．External Memory Management and Simplification of Huge Meshes[J]．IEEE Transactions on Visualization and Computer Graphics，2003，9(4)：525-537．

[8]AGARWAL P K，ARGE L，YIK．I/O-Efficient Construction of Constrained Delaunay Triangulations[C]∥Proceedings of the 13th European Symposium on Algorithms．Mallorca， Spain： Springer Verlag， 2005：355-366．