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

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

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

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

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

1.1 研究區(qū)概況

研究區(qū)位于呂梁市孝義市驛馬鄉(xiāng)下荊封村的某個廠區(qū)。地形表面復(fù)雜，溝谷縱橫，屬大陸性半干旱氣候。研究區(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ù)準(zhǔn)備

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

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

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

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

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

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

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

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

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

(a)植被提取前

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

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

(1)

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

ER=1.4r-g.

(2)

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

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

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

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

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

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

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

圖3 S和θ的幾何意義

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

e.重復(fù)步驟b～d直到所有未分類的點(diǎn)全部判斷完畢。

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

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

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

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

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

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

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

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

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

圖5 箱型圖觀察交集DEM排除大于兩倍標(biāo)準(zhǔn)差的點(diǎn)云變化

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

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

表4 統(tǒng)計(jì)量觀察交集DEM排除大于兩倍標(biāo)準(zhǔn)差的點(diǎn)云變化

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

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

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

可以看到,在處理重疊點(diǎn)后，兩次輸出的點(diǎn)云數(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í)值的程度，其值越小則兩份數(shù)據(jù)之間差異越小。把之前數(shù)據(jù)處理得到的植被的2 040個點(diǎn)、建筑物的205個點(diǎn)、DEM的153 238個點(diǎn)輸入下列公式中，計(jì)算結(jié)果如表5所示。

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

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

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

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

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

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

3.3 Lidar點(diǎn)云和影像點(diǎn)云的DEM驗(yàn)證

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

(a)硬化地表DEM點(diǎn)云

4 結(jié)論

1)當(dāng)不做點(diǎn)云分類時得到的均方根誤差(ERMS)達(dá)到724 mm，而做過點(diǎn)云分類之后，植被層的ERMS可以達(dá)到1 138 mm，相比不做點(diǎn)云分類的724 mm漲幅57%；但是建筑物層的ERMS計(jì)算得到108 mm，降幅達(dá)到85%；DEM通過計(jì)算為627 mm，降幅13%?？梢娋葘Ρ确治龅谝徊剑仨毾确诸?，再對比分析，使用文中點(diǎn)云數(shù)據(jù)處理方法，得到相關(guān)結(jié)論。

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

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