基于LiDAR點云數(shù)據(jù)的等值線繪制方法①

何 睦,李夢珍

1(江西師范大學(xué) 地理與環(huán)境學(xué)院,南昌 330000)

2(江西師范大學(xué) 計算機信息工程學(xué)院,南昌 330000)

1 引言

1.1 研究背景與意義

激光雷達即激光探測與測距技術(shù)(Light Detection And Ranging,LiDAR)是一種逐漸被廣泛應(yīng)用的測繪技術(shù)[1].激光測距系統(tǒng)可以實現(xiàn)傳感器向測量目標發(fā)射激光束,并記錄發(fā)射時刻,激光束到達目標后按照同一光路反射,再由傳感器接收,并記錄接收時刻,那么激光器至測量目標的距離即可通過公式(距離=光速×時間差/2)計算得出.將此激光器安裝在機電設(shè)備上,控制其按一定角度水平或豎直擺動,此時可稱其為激光掃描儀.再將其與高精度動態(tài)載體姿態(tài)測量系統(tǒng)(INS)、高精度動態(tài)GPS差分定位系統(tǒng)(DGPS)相結(jié)合,即可構(gòu)成一個三維激光掃描系統(tǒng),經(jīng)此系統(tǒng)可獲取到ASCII格式的三維點云數(shù)據(jù),記錄每個點三維坐標以此表達目標地物三維空間信息.同時,按照載荷平臺的不同可將其分為地面激光雷達、機載激光雷達、車載激光雷達和星載激光雷達.激光雷達較于傳統(tǒng)遙感測量手段具有幾項明顯的優(yōu)勢：(1) 直接獲取三維坐標信息;(2) 采用主動式測量,不受光照影響;(3) 復(fù)雜地形對雷達測量影響小;(4) 根據(jù)激光的同一脈沖的多次反射對數(shù)據(jù)進行分層可去除植被影響,只獲得地表數(shù)據(jù),提高精度.因此,激光雷達測量技術(shù)被廣泛用于各地學(xué)研究領(lǐng)域,機載雷達主要應(yīng)用于林業(yè)、城市規(guī)劃水利等,車載雷達主要是對地物的側(cè)面進行激光掃描,星載雷達則是用于植被、極地冰川等研究,地面雷達可用于工程測量、建筑物建模等領(lǐng)域.

將激光雷達測量技術(shù)用于測繪行業(yè)在近幾年剛開始起步,而將其用于地形測繪領(lǐng)域也是一個新的探究方向[2].傳統(tǒng)的地形測量必須攜帶GPS儀器、全站儀等眾多負重設(shè)備,還要在許多環(huán)境惡劣的野外進行測量,再經(jīng)內(nèi)業(yè)處理后才能生成地形圖,制作周期長,工作強度和危險性都較大.若能將激光雷達這種遙感測量方式應(yīng)用于現(xiàn)在測量行業(yè)不僅可以簡化外業(yè)流程,縮短工作時間,也可提高內(nèi)業(yè)自動化程度.基于此,本文探討將某處地形的LiDAR點云數(shù)據(jù)作為原始數(shù)據(jù),獲取此處地形的等高線.

1.2 國內(nèi)外研究進展

利用點云數(shù)據(jù)繪制等值線在國內(nèi)外已有一定的研究,一些學(xué)者采用軟件平臺進行了一些實驗;馮梅等[3]分別用GlobaMapper,GeoTIN,Tindem,JX-4C,VirtuoZo和TerraScan六種相關(guān)軟件提取某塊山區(qū)等高線;李鴻軼等[4]使用TerraSolid提取等高線;白志遠[5]則將ArcGIS與LP360相結(jié)合提取等高線.這些實驗均基于一些山區(qū)地形進行,且并沒有選擇大區(qū)域數(shù)據(jù)進行實驗,使用軟件進行數(shù)據(jù)的處理雖有一定的便利性,但對于等值線生成過程中的變量控制性較差,且缺乏普適性.而也有一些學(xué)者研究算法,通過各種方法生成等值線,如吳杭彬等[6]采用迭代凸包的方法提取激光掃描數(shù)據(jù)的分層多細節(jié)等值線,但要經(jīng)過多次迭代凸包,算法效率低.王宗躍等[7]建立三角網(wǎng)跟蹤等高點,采用樣條函數(shù)內(nèi)插等高線;Boissonnat等[8]提出了基于TIN的等高線跟蹤方法.這些算法均有其各自的特色和優(yōu)點,但都只是對固有數(shù)據(jù)加工生成等值線,而未生成不同比例尺下的等高線.但在實際應(yīng)用中,隨著比例尺的減小,地形圖勢必會綜合一些信息,但為了體現(xiàn)地形的形態(tài),應(yīng)保留等值線特征處,才具有清晰性和精確性.因此,探究如何控制生成等值線的變量以期繪制表達更實際的等值線,以及可以根據(jù)不同比例尺的設(shè)定,繪制相應(yīng)的等值線.

2 研究基礎(chǔ)與方法

2.1 研究基礎(chǔ)

LiDAR系統(tǒng)在提供數(shù)據(jù)更新便利的同時,帶來的是點云數(shù)據(jù)的海量性.數(shù)據(jù)量越大對內(nèi)存要求越高,因此便會影響實際生成等值線的效率.同時,若想針對不同比例尺地形情況下,均有其對應(yīng)的等值線形態(tài),也需要對源數(shù)據(jù)進行處理,判斷哪些點對地形特征描述更為重要,若前期對此重要性進行排序,則后期即可根據(jù)比例尺確定某一個節(jié)點,以此為界,則即可在特定比例尺下,繪制特征表述完全的等值線.若要在保證其精度的同時能夠刪去冗余點,費立凡、何津[9]等曾在2006年提出了三維Douglas-Peucker算法,使用此法綜合隨機的三維離散點,無論進行多大尺度的綜合,保留的點均是地貌特征點,相對次要的點將被刪除.2013年,何津、費立凡等[10]提出了基于三維Douglas-Peucker算法綜合三維離散點,再回放等高線,但其根據(jù)某個確定的閾值選取點集合,閾值的確定有一定的主觀性.在本文中,對三維Douglas-Peucker算法進行改進,第一步不是選取所需要的點子,而是對點子進行重要性排隊,對于地形的描述而言,越重要的點越排在前列.在排序之后根據(jù)不同尺度下等高線所需要的綜合程度確定所選點子的數(shù)量.在這種方法里確定取舍點的分界線可以稱其為“剁尾刀”的定位,在多尺度輸出的情況下可以提高幾十倍的效率,并且需要輸出的尺度越多,“剁”得越多,效率越高.

使用三維Douglas-Peucker算法可以實現(xiàn)三維點數(shù)據(jù)的綜合,但要繪制成等值線要經(jīng)過點數(shù)據(jù)的柵格化,進而繪制.于彩霞等[11]提出將點云數(shù)據(jù)先柵格化后提取海岸線防止海岸線的一些曲折及破碎問題.基于此,本文選擇將點云數(shù)據(jù)先插值生成DEM,進而繪制等值線.

2.2 研究方法

(1) 在原有三維Douglas-Peucker算法的基礎(chǔ)上進行改進,加速算法運行的速度以提高效率.將原始算法中對點子的選取改為對點子的排序,在排序之后對其進行選取提高了控制性.傳統(tǒng)的分裂點選取點子的方法是根據(jù)地形圖制圖要求確定閾值大小,進而確定點到面的距離,大于閾值的則保留,小于閾值的刪除.而現(xiàn)今,如何確定每個點子的特征重要性,如何客觀地對其排序,到底保留哪些點可以回放出特征完全的等值線,即為本文的研究內(nèi)容.

(2) 在能夠?qū)崿F(xiàn)點子重要性排隊后,保留原始數(shù)據(jù)的多少比率下的數(shù)據(jù)才能生成符合要求的等值線,也就是探究已知量比例尺與未知量數(shù)據(jù)保留比率之間的關(guān)系.而在點子數(shù)據(jù)生成等值線過程時,柵格的分辨率也會影響最終等值線繪制的綜合程度.所以探究生成柵格分辨率與等值線效果的關(guān)系也成為另一個內(nèi)容.

圖1 技術(shù)流程圖

2.3 前期數(shù)據(jù)準備

使用地面激光雷達掃描技術(shù)獲取地形數(shù)據(jù),首先進行野外數(shù)據(jù)的采集,考察測區(qū)環(huán)境并確定掃描儀、測站數(shù)和標靶的位置.選定合適的測站數(shù)既要保證最終的數(shù)據(jù)能表達測區(qū)完整地形,也要保證數(shù)據(jù)冗余程度低.之后進行數(shù)據(jù)的初始處理,要逐步進行點云匹配、消除噪聲、坐標變換、圖像分割.當點云數(shù)據(jù)來自不同測站點時,需要將其匹配到同一個坐標系,通用的方法是通過GPS獲取掃描站點的地理坐標,并將點云數(shù)據(jù)轉(zhuǎn)換到同一個地理坐標系中.配準完成后,手工對圖像進行重采樣,消除噪聲如電線桿、植被等,最后可以對圖像進行分割,選取出感興趣的區(qū)域.

本文所選用的ASCII格式的點云數(shù)據(jù),獲取到的是使用單測站所測得的礦山數(shù)據(jù),因此其表現(xiàn)的是礦山的單側(cè)面.可將此文件在ArcMap中轉(zhuǎn)化為Shapefile文件,手動添加X、Y、Z字段,并使用自動計算器計算每個點的三維坐標并賦值.可將此類數(shù)據(jù)在ArcScene進行三維顯示.

3 算法設(shè)計與實現(xiàn)

3.1 算法的設(shè)計

(1) 創(chuàng)建基面

遍歷Shapefile文件中所有的點,獲取每個點的三維坐標,進行依次比較后貯存X,Y,Z三方向的最大值與最小值.記Xmin,Xmax分別為X方向上的最小值和最大值,Y方向與Z方向以此類推.則O(Xmin,Ymin,Zmin)、A(Xmax,Ymin,Zmin)、B(Xmin,Ymax,Zmin)分別構(gòu)成OA和OB兩個矢量線,以此構(gòu)成首基面OAB;C(Xmax,Ymin,Zmax)、D(Xmin,Ymin,Zmax)、A(Xmax,Ymin,Zmin)分別構(gòu)成CD和CA兩個矢量線,以此構(gòu)成首基面CAD;E(Xmax,Ymax,Zmin)、F(Xmax,Ymax,Zmax)、A(Xmax,Ymin,Zmin)分別構(gòu)成EA和EF兩個矢量線,以此構(gòu)成首基面EFA;OAB,CDA,EFA為三個正交的首基面,如圖2所示.

圖2 3個正交首基面

(2) 單首基面(一合一)3D點子重要性排序

單首基面對3D點子重要性排序的方法依據(jù)點到面的距離,而三維點到面的距離與二維點到線的距離相對應(yīng),因此利用二維曲線結(jié)構(gòu)圖更為直觀的優(yōu)點說明三維數(shù)據(jù),如圖3所示.

圖3 利用二維曲線對三維數(shù)據(jù)重要性排序說明的結(jié)構(gòu)圖

二維下對點子進行重要性排隊的步驟如下(流程見圖 4)：

Setp 1.無條件地將曲線的首末點(即點1和點2)選取(已被選取點標記為紅色,下同);

將連接試驗曲線首末兩點的直線作為首基線,考察原曲線介于首末點之間的所有點,可找到點線距為最大值的點(即點3),故點3被選取;如圖4(a).

圖4 點子重要性排序示意圖

Setp 2.介于已選點1和3之間、3與2之間,分別考察原曲線介于局部首末點之間的所有點,可找出最大值的點(即點4和點19標記為黑色),將點4插入1和3之間,將點19插入3和2之間,如圖4(b).

Setp 3.在被插入的中間點點集中(此處指的是點4和點19),由于點4的點線距大于點19的,所以點4被選取(已被選取的點4被標成紅色),如圖4(c).

Setp 4.介于1與4之間、4與3之間,分別考察原曲線介于局部首末點之間的所有點,可找出點線距為最大值的點(即點5、15),并插入各區(qū)間,如圖4(d).

Setp 5.在被插入的中間點點集中(此處指的是點5、15、19等黑色標號所指),由于點5的點線距為最大,所以點5被選取(點5被標成紅色),如圖4(e).

Setp 6.如此循環(huán)往復(fù),其依照的口訣便是：① 點子選上變紅點;② 紅點之間插黑點;③ 黑點之間要PK;④ PK得勝變紅點;⑤ 全是紅點就結(jié)束;否則轉(zhuǎn)Step 2不要停.黑點逐個變成紅點,變成紅點的先后順序也就是其重要性順序.最終點子順序如圖4(f).

將上述排序方法轉(zhuǎn)化為三維,以O(shè)AB首基面為例,則O點為原點分別與上述方法中的局部首末點相連構(gòu)成一個矢量面,此時計算點線距轉(zhuǎn)化為計算點面距.若以O(shè)AB為一個首基面,設(shè)此面的一般方程式為其中a、b、c、d均為已知量.點到此基面的距離由公式(1)計算：

其中,a、b、c、d為OAB基面的一般方程式的系數(shù),X、Y、Z為文件中點的三維坐標.

(3) 三合一3D點子排序

按照上一步的方法,獲得3組分別以O(shè)AB、CDA、EFA為基面得到的順序點隊列I、II、III,在這些隊列中,從左到右,點的重要性依次降低,對每個隊列中的點依次編號.則點子的序號越大,說明此點的重要性越小.以I為依據(jù),依次讀取其中的點對應(yīng)的序號,并在II和III中找到相同的點,獲取其序號.將同一點在此3個隊列中的序號相加得到數(shù)值sum,表示其重要性,數(shù)值越小則表示其重要性越大.按照值從小往大的順序依次排列其對應(yīng)的點,將這些點按此順序存入到隊列中,由此將體現(xiàn)3個方向特征重要性的點序列融合成一個點隊列.

(4) 剁尾刀

以三合一過程獲取的點隊列為依據(jù),確定所要保留的點率,由公式(2)得到最終保留點的數(shù)量.

其中,count為源數(shù)據(jù)點的數(shù)量,r為保留比例,取值范圍為0～1,n為最終保留點的數(shù)量.此時確定的數(shù)量即為要截取點子列表中的數(shù)據(jù)量,從隊列的左端開始往右獲取剁尾刀的保留點界限,直接得到左側(cè)應(yīng)保留的點集合.

3.2 算法實現(xiàn)

依托 ESRI的 ArcGIS平臺,使用 C#語言,在VisualStudio 2010環(huán)境下,進行Add-In定制開發(fā)ArcMap擴展桌面應(yīng)用程序功能.將原始格式為ASCII點云數(shù)據(jù)經(jīng)預(yù)處理后的得到shapefile格式的點云數(shù)據(jù)后創(chuàng)建X、Y、Z坐標字段.將shapefile格式的點云數(shù)據(jù)作為源數(shù)據(jù),設(shè)定固定比例的特征保留率,經(jīng)上述算法處理后得到結(jié)果點云數(shù)據(jù).

利用COM組件技術(shù),實現(xiàn)ArcObjects開放的特定接口,實現(xiàn)3D Analyst Tools工具欄下的Topo To Raster和Contour功能,繪制等值線,再實現(xiàn)可視化接口使其直接在ArcMap中顯示.在實現(xiàn)點數(shù)據(jù)綜合后運用二次開發(fā)技術(shù)實現(xiàn)等值線繪制的自動一體化過程.

4 算法設(shè)計與實現(xiàn)

測試數(shù)據(jù)區(qū)域為LiDAR系統(tǒng)單測站所采集的礦山點云數(shù)據(jù),其中共有29 677個點.在點云數(shù)據(jù)經(jīng)過本算法生成等值線的過程中,點集合保留的比率和DEM生成的分辨率是算法的兩個可控變量,成為影響最終等值線的生成結(jié)果的主要因素.點云數(shù)據(jù)在低精度DEM時生成等值線時,海量數(shù)據(jù)的優(yōu)勢得不到體現(xiàn),繪制等值線時與傳統(tǒng)方法繪制差異不顯著,在1 m至5 m精度范圍內(nèi)可以較好的反映地形信息[12,13],故選擇在保留原始數(shù)據(jù)100%的基礎(chǔ)上,分別生成分辨率為1.0 m、1.5 m、2.5 m的DEM對比分析繪制等值線差異.按梯度選擇出對比分析保留10%、50%、100%的點數(shù),生成DEM分辨率為1.5 m,繪制的等值線結(jié)果差異.最終對比保留點率和分辨率這兩個參數(shù)作為變量對回放等值線的影響.

5 結(jié)論與討論

5.1 算法實現(xiàn)

由點子數(shù)據(jù)繪制等值線的過程中,影響最終結(jié)果的主要因素為點子的數(shù)量和DEM的分辨率,點云數(shù)據(jù)中的點子數(shù)量和DEM的分辨率與繪制等值線結(jié)果的精度成正比,為對比這兩個因素對等值線繪制的影響,先固定一個變量再變化另一變量的不同情況進行對比分析得到其影響程度.

(1) 不同分辨率

導(dǎo)入源數(shù)據(jù),為保證點子數(shù)量不變,故使用3DDP算法保留100%的點集合分別轉(zhuǎn)換成網(wǎng)格邊長1 m、1.5 m、2.5 m的DEM數(shù)據(jù),進而繪制出等值線,如圖5、6、7所示.

圖5 由網(wǎng)格邊長為1 m的DEM轉(zhuǎn)成的等高線

圖6 由網(wǎng)格邊長為1.5 m的DEM轉(zhuǎn)成的等高線

圖7 由網(wǎng)格邊長為2.5 m的DEM轉(zhuǎn)成的等高線

通過實驗可以看出,利用重要性排序后的固定點子,轉(zhuǎn)換成邊長變量為1.5 m的柵格DEM,回放出的等高線比較符合工程測量與繪圖的要求(見圖6).由前三幅圖中可見,根據(jù)綜合后的點數(shù)據(jù)所繪制的等高線能夠很好地表現(xiàn)所測范圍的地形,且大體趨勢一致.隨著生成DEM的格網(wǎng)增大,分辨率逐漸降低,綜合程度增大,等值線趨于平滑,且碎屑多邊形逐漸減少,提升主觀視覺效果.將這3種情況疊合在一起的結(jié)果如圖8所示,圖幅的中間部位三者吻合度較高,在圖幅的上端和下端稍有差異,但并不影響整體的地形表達效果.

圖8 三者疊加

(2) 不同比率保留點集合

輸入源數(shù)據(jù),使用3D-DP算法,根據(jù)設(shè)定的10%、50%、100%比率分別保留2968、14 839、29 677個點數(shù)據(jù),插值為分辨率為1.5 m的DEM,進而繪制成等高距為5 m的等值線.三類比率對應(yīng)的等值線如圖9、10、11所示,將此三種數(shù)據(jù)疊加如圖12所示.

圖9 保留100%點的回放等高線

通過實驗可以看出,利用保留率變量為100%的最重要的點子,轉(zhuǎn)換成固定邊長為1.5 m的柵格DEM,回放出的等高線比較符合工程測量與繪圖的要求(見圖10).由圖11等高線可以看出,隨著保留點子數(shù)量的減少,即使是只保留至10%的點子,依舊能夠表現(xiàn)所測地區(qū)的基本地貌形態(tài).且由等高線對比分析時可以看出保留不同程度點子所繪制的等值線差異不大,可說明通過3D-DP算法綜合的點子能夠保留下最能表達地形特征的點,在一定限度的點數(shù)下,可對表達同一地形的冗余點進行刪除.

圖10 保留50%點的回放等高線

圖11 保留10%點的回放等高線

圖12 三者疊加

實驗一和實驗二分別降低分辨率和減少點子數(shù)量,均加大了綜合程度,由實驗可知減少點子數(shù)量相對于降低分辨率所繪制的等值線吻合度更高,差異更小.由此可以得出,通過3D-DP算法保留特征點可以滿足不同綜合程度下的等值線繪制,也因此驗證3D-DP中點子的順序排隊和“剁尾刀”的方法易于通過給定點子保留百分比控制等高線的綜合程度,與之前國內(nèi)外其他研究相比[2-8],本方法充分利用了點云數(shù)據(jù)的優(yōu)勢,且體量更小,易于調(diào)用,具有更好的普適性,在比例尺為1：500-1：20 000的范圍內(nèi)可以實現(xiàn)小范圍內(nèi)多尺度地圖制圖綜合的自動化,可輸出不同比例尺下的等值線,與已有平臺兼容良好,可嵌套在ArcGIS平臺中,結(jié)果可直接在ArcMap中顯示,對于多尺度的等高線輸出更為方便、高效.

5.2 研究結(jié)論與展望

(1) 對于綜合后的等高線,與實測數(shù)據(jù)對比以后,確定其回放的等高線可以滿足對應(yīng)比例尺下的地形圖要求,可以實現(xiàn)小范圍內(nèi)地圖制圖綜合程度控制的自動化.

(2) 本文使用了3D D-P算法對數(shù)據(jù)進行綜合,并可以選擇一定分辨率內(nèi)的點子數(shù)據(jù)集合以減少運算時間,提高等值線生成效率.但是,若要面對大批量數(shù)據(jù)自動生成等值線時,必須對數(shù)據(jù)進行分塊處理,除點位數(shù)據(jù)保留比率與DEM分辨率之外,本算法還應(yīng)對比不同類型的典型地形之間的等值線繪制效率與精度的差異進行比較,來進一步提高算法的科學(xué)性和完整性,以及引入并行運算技術(shù),來實現(xiàn)大范圍高分辨率、多尺度等高線的地圖制圖自動化.

(3) 在僅有三維坐標的情況下生成等值線,若要真實地描述地形,在考慮地貌結(jié)構(gòu)線的參與下,以規(guī)則約束等值線的生成會使結(jié)果更貼合實際地貌.特別是處理山區(qū)中地貌與水系之間的協(xié)調(diào)套合,使其真實體現(xiàn)地形.

1 馬鴻超.激光雷達測量技術(shù)在地學(xué)中若干應(yīng)用.地球科學(xué)—中國地質(zhì)大學(xué)學(xué)報,2011,36(2)：347-354.

2 劉青峰,劉賀,王東.3維激光掃描技術(shù)在地形測量中的應(yīng)用.測繪與空間信息,2013,36(5)：153-155.

3 馮梅,鐘斌.基于LiDAR點云自動生成等高線的方法研究.測繪與空間地理信息,2012,35(6)：87-90,93.

4 李鴻軼,任淑娟,雷蕾.LiDAR數(shù)據(jù)提取等高線的方法研究.測繪通報,2013,(5)：59-60.

5 林洪文,涂丹,李國輝.基于統(tǒng)計背景模型的運動目標檢測方法.計算機工程,2003,29(16)：97-99,108.[doi：10.3969/j.issn.1000-3428.2003.16.040]

6 吳杭彬,劉春.激光掃描數(shù)據(jù)的等值線分層提取和多細節(jié)表達.同濟大學(xué)學(xué)報(自然科學(xué)版),2009,37(2)：267-271,276.

7 王宗躍,馬宏超,彭檢貴,等.基于LiDAR數(shù)據(jù)生成光滑等高線.武漢大學(xué)學(xué)報·信息科學(xué)版,2010,35(11)：1318-1321.

8 Boissonnat JD.Shape reconstruction from planar cross sections.Computer Vision,Graphics,and Image Processing,1988,44(1)：1-29.[doi：10.1016/S0734-189X(88)80028-8]

9 費立凡,何津,馬晨燕,等.3維Douglas-Peucker算法及其在DEM自動綜合中的應(yīng)用研究.測繪學(xué)報,2006,35(3)：278-284.

10 何津,費立凡,黃麗娜,等.三維Douglas-Peucker算法的等高線間接綜合方法研究.測繪學(xué)報,2013,42(3)：467-473.

11 于彩霞,許軍,許堅,等.一種從LiDAR點云中提取海岸線的新方法.測繪通報,2015,(5)：66-68.

12 劉洪鵠,龍翼,嚴冬春,等.不同分辨率DEM提取三峽庫區(qū)地形參數(shù)的精度研究.長江科學(xué)院院報,2010,27(11)：21-24.[doi：10.3969/j.issn.1001-5485.2010.11.005]

13 國家測繪局.CH/T 9008.2-2010基礎(chǔ)地理信息數(shù)字成果1：500 1：1000 1：2000 數(shù)字高程模型.北京：測繪出版社,2010.