李 濤，常 江,胡東升，廉旭剛，呂俊沛

(1.華陽新材料科技集團有限公司，山西 陽泉 045000；2.太原理工大學 礦業(yè)工程學院，太原 030024)

目前，無人機搭載激光雷達模塊采集點云數(shù)據(jù)，獲取地表高精度三維坐標，相比數(shù)碼航空攝影測量有獨特優(yōu)勢[1]。首先，激光雷達直接獲取地物表面高精度三維坐標，不受山體、建筑物陰影的影響，單純數(shù)碼航空攝影測量在陰影區(qū)域判讀困難。其次，激光雷達的多次回波對植被具有較好的穿透性，而植被覆蓋一直是數(shù)碼航空攝影不能解決的問題。最后，攝影測量提高精度需經(jīng)過后差分技術(shù)[2-3]和SFM技術(shù)[4-5]，即通過設計控制點的增多與布設位置的調(diào)整，通過拍攝多視角照片結(jié)合特征匹配算法，DEM精度可以達到±9 cm，而原始激光點云即可達到±10 cm的高程精度。同時,有多篇文章提出Lidar獲取高精度DEM的關(guān)鍵技術(shù)，大致可以分為兩類，即通過設計改進濾波算法和通過調(diào)整采集方案。例如,提出過綠減過紅指數(shù)提取植被[6-7]、距離限制濾波提取地面點[8]等改進算法；提出從設備選擇、點云密度設計、植被覆蓋密集山區(qū)數(shù)據(jù)獲取方法、點云數(shù)據(jù)分類組合算法、空白區(qū)處理等方面[9-10]進行探討,并提出改進方法。

目前，對一種點云數(shù)據(jù)(Lidar點云或者攝影測量點云)精度研究比較多，但是缺乏Lidar點云和攝影測量點云對同一片地區(qū)DEM的對比誤差研究，定量化描述Lidar點云和攝影測量點云的誤差缺乏方法也缺乏實際應用案例。本文提出一種點云對比方法，并結(jié)合某地區(qū)的實際數(shù)據(jù)進行案例演示，為后續(xù)此類型研究提供參考方案。

1 研究區(qū)概況和數(shù)據(jù)準備

1.1 研究區(qū)概況

研究區(qū)位于呂梁市孝義市驛馬鄉(xiāng)下荊封村的某個廠區(qū)。地形表面復雜，溝谷縱橫，屬大陸性半干旱氣候。研究區(qū)大部分地區(qū)植被茂密，除了一條山路和廠區(qū)外，只有少量地表裸露，可以滿足不同地貌條件之間高程數(shù)據(jù)的對比。

地理位置東經(jīng)111.604 2°—111.608 1°，北緯36.957 2°—36.960 0°之間，高程最低947 m，最高992 m，研究區(qū)面積大約0.055 km2。如圖1所示。

圖1 研究區(qū)概況

1.2 數(shù)據(jù)準備

數(shù)據(jù)采集使用飛馬D2000無人機航測平臺，一次飛行搭載D-LIDAR2000激光雷達模塊，重疊度設置為25%。另一次飛行搭載索尼 A6000相機，有效像素2 430萬，航向和旁向重疊分別為80%和60%，采用網(wǎng)絡RTK/PPK高精度融合POS+免像控方案。數(shù)據(jù)采集時間是2021年9月20日，搭載索尼 A6000相機采集的數(shù)據(jù)，可滿足高精度1∶500地形圖精度需求。搭載D-LIDAR2000激光雷達模塊采集的數(shù)據(jù)可以直接用于本次的研究。雷達模塊參數(shù)如表1所示。航測模塊參數(shù)如表2所示。

表1 D-LIDAR2000激光雷達模塊的參數(shù)

表2 SONY a6000航測模塊的參數(shù)

2 點云數(shù)據(jù)處理

2.1 數(shù)據(jù)修正

在本次研究中，首先需要清除數(shù)據(jù)在采集過程中或者使用相關(guān)軟件處理過程中產(chǎn)生的異常點、孤立點[11]。使用相關(guān)濾波算法時，要注意清除低點(不是正常的低矮地物)。使用點云分類工具，設置相關(guān)參數(shù)分離低點和孤立點，為了保證清除完整，可以多次執(zhí)行該命令。特殊噪點需要手動清除，使用分類工具classify above line，從側(cè)視圖中觀察發(fā)現(xiàn)噪點。

2.2 數(shù)據(jù)分層

對已經(jīng)清除噪點的兩份點云數(shù)據(jù)進行點云分類，兩份數(shù)據(jù)分別劃分為植被層、建筑物層、地面層(地面層生成DEM)，主要使用點云分類工具，設置不同的參數(shù)，保存對應層的數(shù)據(jù)。兩份數(shù)據(jù)提取相同層的時候盡量保證參數(shù)一致，使得兩份數(shù)據(jù)相同層的點云在數(shù)量和質(zhì)量上保持一致。每層數(shù)據(jù)的保存需要規(guī)定格式，方便后續(xù)的分析處理。

1)攝影測量點云提取植被層。使用點云分類工具的by vegetation index功能，因為攝影測量點云是帶有顏色(RGB值)的，其原理[6]主要是通過過綠減過紅指數(shù)(EG-ER)把RGB值當參數(shù)，通過計算結(jié)果和閾值來判斷是否是綠色，再將密度值和高程差值作為輸入特征變量[7]，使用支持向量機(SVM)算法做點云分類，來識別是否是植被。設置的min value和max value含義是規(guī)定綠色的范圍，目標找到這范圍內(nèi)的綠色點。如圖2所示。

(a)植被提取前

EG-ER作為一種改進的顏色指數(shù)通過將EG指數(shù)圖像與ER指數(shù)圖像相減，發(fā)現(xiàn)基于閾值的EG-ER指數(shù)可以較好地將植被與背景分離，計算公式如下：

EG-ER=3g-2.4r-b.

(1)

其中ER為過紅指數(shù)，計算方法如下：

ER=1.4r-g.

(2)

上式的R,G,B值是歸一化后的值，即0～255歸一到0～1之間。

2)點云提取地面點層制作DEM。使用點云分類工具的ground功能，該功能的原理是發(fā)現(xiàn)點云之間的距離和角度的關(guān)系[12]。因為不同地物構(gòu)建的三角網(wǎng)在角度、邊長、高度上有差異，通過設置相關(guān)參數(shù)，排除地表的植被和建筑，只留下地面點。分離地面點的原理有兩個，距離限制濾波[8]分類地面點和非地面點，角度濾波去除地物的側(cè)面信息。然后依據(jù)局部地形設置動態(tài)閾值，以表面擬合區(qū)域生長算法擴充地面種子點，循環(huán)迭代逐漸逼近真實地面[13]。

濾除起伏突然比較大的點，比如建筑物點，植被點等，使用角度濾除，一般計算θ值的大小，計算公式如下：

獲得地面種子點之后，便可以進行表面擬合區(qū)域生長算法擴充地面種子點，其流程可以描述為：

a.利用地面種子點構(gòu)建TIN；

b.選擇一個未分類點，尋找該點所在的三角形；

c.計算該點到三角形平面的距離S和該點與三角形3個角點的夾角，并選出最大角度，其幾何意義如圖3所示；

圖3 S和θ的幾何意義

d.根據(jù)人工設置的最大距離閾值Smax和最大角度閾值θmax，再規(guī)定區(qū)間，將該點加入地面點；

e.重復步驟b～d直到所有未分類的點全部判斷完畢。

上述原理在提取過程中體現(xiàn)為設置不同的參數(shù)，Max building size設置需要根據(jù)實地情況填寫，必須大于當?shù)氐淖畲蠼ㄖ娣e。Iteration angle和Iteration distance兩個參數(shù)是三角網(wǎng)的參數(shù)，代表角度和距離。城區(qū)設置一般偏小，比如4°和0.6 m；山區(qū)設置一般偏大，比如6°和1.4 m。如果得到過濾點云密度太低，需要選擇加密點云ground層，因為ground層的點云密度越大，使用Export lattice Model工具插值得到規(guī)則格網(wǎng)DEM越精確。

2.3 DEM點云數(shù)據(jù)處理

按照2.2描述進行數(shù)據(jù)分層，一共得到6個文件，分別是Lidar點云的building.txt，tree.txt，DEM.txt和攝影測量點云的building.txt，tree.txt，DEM.txt。然后主要使用python語言，寫代碼實現(xiàn)數(shù)據(jù)的分析處理。

主要思想是，由于兩份數(shù)據(jù)的坐標系相同，同一地區(qū)范圍得到的同一種地物類型的點云，必定存在相同坐標點，兩幅點云的相同坐標點的高程差應該是差異較小。但是由于數(shù)據(jù)獲取方式不同、數(shù)據(jù)處理方式不同和數(shù)據(jù)類型不同，兩幅點云的高程存在不容忽視的差異，所以需要通過以下方法量化分析二者差異。如圖4所示。

圖4 算法實現(xiàn)方法圖

使用Lidar點云和影像點云的DEM.txt來舉例：通過上述處理，得到Lidar點云大約20萬個點，攝影測量點云17萬個點，使用python的pandas庫可以實現(xiàn)找相同坐標點的操作，并且可以最大程度保留原始數(shù)據(jù)[14-15]。經(jīng)過上述操作，保留大約16萬對相同坐標點，是原始數(shù)據(jù)Lidar點云數(shù)量的80%，是原始數(shù)據(jù)攝影測量點云的95%。然后對新表的列名重新命名，進而計算Lidar點云高程減去攝影測量點云高程的差值。如表3所示。

表3 兩份DEM的交集數(shù)據(jù)以及高程差

從Lidar-photo這一列中，可以看到，兩份DEM的高程差是普遍存在的。由于軟件處理植被點云可能不完整，或者由于無人機飛行初始和結(jié)尾狀態(tài)的不穩(wěn)定性[9]，導致出現(xiàn)的高程差比較大，稱之為粗差點，這些點是需要排除的。

選擇剔除的方法一般是認為兩倍標準差是數(shù)據(jù)的合理范圍，大于兩倍標準差的可以認為是粗差，選擇剔除。如圖5所示。

圖5 箱型圖觀察交集DEM排除大于兩倍標準差的點云變化

箱型圖中黑色圓圈點表示大于1.5倍標準差的數(shù)據(jù)，矩形方框表示數(shù)據(jù)的75%集中的區(qū)間，綠色的線表示平均值。下表是處理前后的關(guān)鍵統(tǒng)計量的變化，平均值變化較小，進一步說明排除大于兩倍標準差的方法是合理的。

經(jīng)過上述處理，剔除5%的粗差點云，得到的點云數(shù)量占原來的95%。如表4所示。

表4 統(tǒng)計量觀察交集DEM排除大于兩倍標準差的點云變化

從數(shù)據(jù)中可以看出，已經(jīng)把處理前最大高差19 m的數(shù)據(jù)排除，但存在最大高差5 m或者4 m，而選擇保留是因為在實際中可能存在5 m或4 m的高差?，F(xiàn)在兩份DEM高差數(shù)據(jù)大部分的高差有幾十厘米甚至更小，從第一四分位(25%)為-0.583 75 m和第三四分位(75%)為0.147 m可以看出。因此，待研究的數(shù)據(jù)整體是比較合理的，研究結(jié)果是具有參考意義的。

2.4 植被和建筑點云數(shù)據(jù)處理

數(shù)據(jù)處理辦法和2.3一樣，但是需要注意因為Lidar點云的穿透性，是存在同一個x、y坐標有多個高程值，所以第一步需要對保存的植被層刪除重疊點，保留重疊的多個點的最高值，把最高值認為是植被高度。對于建筑物點云的重疊情況，選擇重疊點云的平均值作為建筑物的高程。

可以看到,在處理重疊點后，兩次輸出的點云數(shù)量是不一樣的。比如，植被層，第一次輸出9 429 398個，第二次輸出9 352 880個，重疊率0.81%。

3 數(shù)據(jù)對比分析

3.1 平均絕對誤差(EMA)的比較

平均絕對誤差，全稱是Mean Absolute Error，EMA的值越小，說明影像高程數(shù)據(jù)越靠近Lidar高程數(shù)據(jù)。EMAP[17]是EMA的變形，它是一個百分比值，表示偏離真實值的程度，其值越小則兩份數(shù)據(jù)之間差異越小。把之前數(shù)據(jù)處理得到的植被的2 040個點、建筑物的205個點、DEM的153 238個點輸入下列公式中，計算結(jié)果如表5所示。

表5 植被、建筑、DEM的EMA和EMAP的對比表

數(shù)據(jù)結(jié)果說明，不論影像點云還是Lidar點云，建筑物的精度是最高的。植被和DEM差異比較小，但是由于DEM是地面點合成的，從物理的角度來說確實應該比植被的精度要高，所以可以說明這份數(shù)據(jù)處理過程是合理的，其對比分析具有一定參考價值。

3.2 均方根誤差(ERMS)的比較

均方根誤差(ERMS)反映測量數(shù)據(jù)偏離真實值的程度，ERMS越小，表示測量精度越高。把Lidar點云數(shù)據(jù)當作真值，把影像點云數(shù)據(jù)當作觀測值。同樣把之前數(shù)據(jù)處理得到的植被的2 040個點、建筑物的205個點、DEM的153 238個點輸入到上述代碼中，得到以下計算結(jié)果，如表6所示。

表6 植被、建筑、DEM的ERMS的對比表

從ERMS結(jié)果中仍然可以看出，建筑物點云的數(shù)據(jù)精度最高，和其他兩種數(shù)據(jù)的差異是比較大的，植被和DEM比較接近，DEM稍優(yōu)于植被。所以，再次證明不論Lidar點云還是影像點云，在做精度對比的時候一定要先將地物進行分類。

3.3 Lidar點云和影像點云的DEM驗證

可以明顯看到在平坦或者規(guī)則地區(qū)，Lidar的DEM和影像的DEM高度重合，但是在山地或者植被覆蓋的地方二者差異還是比較明顯的，如圖6所示。

(a)硬化地表DEM點云

4 結(jié)論

1)當不做點云分類時得到的均方根誤差(ERMS)達到724 mm，而做過點云分類之后，植被層的ERMS可以達到1 138 mm，相比不做點云分類的724 mm漲幅57%；但是建筑物層的ERMS計算得到108 mm，降幅達到85%；DEM通過計算為627 mm，降幅13%。可見精度對比分析第一步，必須先分類，再對比分析，使用文中點云數(shù)據(jù)處理方法，得到相關(guān)結(jié)論。

2)經(jīng)過傳統(tǒng)攝影測量得到的影像點云，按照文中的分類方法，分類得到的硬化地表地區(qū)的點云，其精度是可以認為比較接近Lidar點云，所以在建筑不太密集的城區(qū)，可以繼續(xù)發(fā)揮傳統(tǒng)攝影測量在成本和效率方面的優(yōu)勢。

3)在實際應用案例中，經(jīng)過該方法處理后，DEM點云數(shù)量損失為5%，而植被和建筑物點云數(shù)量損失稍大，可以說明該對比方法的合理性；而且數(shù)據(jù)處理過程中，在同一坐標系下多次找兩份點云數(shù)據(jù)的相同坐標點，確保對比方法的嚴密性，得到的結(jié)論具有實際參考意義。