李沛婷，趙慶展，田文忠，馬永建，陳洪

(1石河子大學(xué)信息科學(xué)與技術(shù)學(xué)院/國家遙感中心新疆兵團(tuán)分部/兵團(tuán)空間信息工程技術(shù)研究中心,新疆 石河子 832003; 2 石河子大學(xué)機(jī)械電氣工程學(xué)院,新疆 石河子 832003)

機(jī)載激光雷達(dá)(light detection and ranging,LiDAR)因其有非接觸主動式和可以快速獲取物體高精度點云等優(yōu)點,已成為獲取點云的一種新穎方式,其在智慧城市、資源調(diào)查和基礎(chǔ)測繪等領(lǐng)域發(fā)揮著越來越重要的作用[1]。機(jī)載LiDAR集成了激光掃描儀、動態(tài)差分全球定位系統(tǒng)(global positioning system,GPS)、慣性導(dǎo)航系統(tǒng)(inertial navigation system,INS)數(shù)碼相機(jī)等先進(jìn)技術(shù)[2],其中定位定姿系統(tǒng)(positioning and orientation system,POS)控制確定飛行平臺的空間位置和姿態(tài),通過激光掃描儀獲取光學(xué)中心與地物目標(biāo)之間的距離信息,利用數(shù)碼相機(jī)捕捉圖像信息[3]。因為機(jī)載LiDAR采集的點云具有數(shù)據(jù)量大、冗余高、點密度不均勻等特點,所以研究點云數(shù)據(jù)處理是后期將其應(yīng)用在不同領(lǐng)域的基礎(chǔ)[4],其中點云處理包括預(yù)處理和后處理,預(yù)處理主要是點云濾波,目的是較準(zhǔn)確分離地面點與非地面點,且利用地面點云生成數(shù)字高程模型(digital elevation model,DEM)。DEM是地理信息系統(tǒng)發(fā)展的基礎(chǔ)數(shù)據(jù),其在很多應(yīng)用領(lǐng)域上都是重要的基本需求和地形產(chǎn)品[5],如生成荒漠區(qū)高精度DEM,是開展荒漠區(qū)植被、地形變換和生態(tài)系統(tǒng)分布等研究的重要基礎(chǔ)。王智等[6]結(jié)合亞洲數(shù)字高程數(shù)據(jù)和新疆土地利用數(shù)據(jù),實現(xiàn)對新疆山地-綠洲-荒漠中植被覆蓋變化時空特征研究;李向婷等[7]結(jié)合機(jī)載LiDAR測量得到的DEM數(shù)據(jù)和MODIS影像數(shù)據(jù),提取新疆荒漠稀疏植被覆蓋度。

目前,國內(nèi)外學(xué)者一方面直接使用點云的三維坐標(biāo)屬性,或者通過三維坐標(biāo)獲取點云法向量等間接特征實現(xiàn)濾波,另一方面也逐漸通過對點云回波次數(shù)和強度的深入研究來提高點云濾波精度。Pirotti F等[8]結(jié)合地面激光掃描儀的回波次數(shù),實現(xiàn)點云濾波得到DEM;Reymann等[9]針對不同試驗區(qū)的點云數(shù)據(jù),得到結(jié)合點云回波次數(shù)和回波強度屬性信息以實現(xiàn)點云分類,比僅使用點云的幾何信息更加精確;Ghamisi P等[10]采用改進(jìn)的支持向量機(jī),結(jié)合點云高度變化、回波強度、激光掃描回波和擬合的殘差實現(xiàn)點云分類;龔亮等[11]和曾靜靜等[12]利用點云強度、回波次數(shù)和高度三種屬性信息,提取機(jī)載LiDAR點云的道路,結(jié)合屬性信息可提高點云分類精度;劉凱斯[13]根據(jù)首末次回波高度值去除植被點,不僅提升點云的濾波效率,而且以回波次數(shù)、高度值和回波強度三個屬性信息作為地物分類的約束條件構(gòu)建決策樹實現(xiàn)點云分類;樊敬敬[14]提出了直接利用首尾兩次回波強度之差提取植被信息,結(jié)果表明該方法提取植被的效果較好,避免了使用兩次回波高度差提取過程中的植被和建筑物邊緣混合問題。

由上述研究可見,結(jié)合點云回波次數(shù)、回波強度和點云高度等屬性信息,可以快速高效地實現(xiàn)點云濾波和點云地物分類等處理。在實際應(yīng)用中,如果直接使用原始點云的回波屬性進(jìn)行濾波處理,將會大大降低運算速度和效率,因此，本文在進(jìn)行點云濾波之前,首先采用箱形圖對點云屬性的離群情況進(jìn)行分析,在去除大量離群點云的基礎(chǔ)上,再通過反復(fù)建立三角網(wǎng)來實現(xiàn)點云濾波。

1 材料與方法

1.1 無人機(jī)平臺

以瑞士AeroScout公司提供的Scout B1-100單旋翼油動無人直升機(jī)作為飛行平臺,其空機(jī)質(zhì)量50 kg,有效載荷18 kg,對應(yīng)的長寬高為3.3 m×1.0 m×1.3 m,飛行時間可達(dá)90 min。

POS系統(tǒng)是由Oxford Technical Solutions(OXTS)公司提供的Survey+2,其采用GNSS / INS緊密耦合的方法實現(xiàn)飛行器空間位置和姿態(tài)的確定,刷新頻率為100 Hz,定位精度為2 cm,橫滾和俯仰精度為0.03°。

以Riegl公司提供的VUX-1全波段激光雷達(dá)為飛行器,激光掃描儀的掃描方式為線性,波長是近紅外,掃描速度范圍10～200 scan/s,最大掃描角度為330°,激光脈沖頻率550 kHz,測量精度15 mm,不僅可以記錄16 bit的回波強度,而且可以記錄無限次的回波次數(shù)。

在數(shù)據(jù)獲取過程中,因為激光掃描儀VUX-1和Survey+2是固定連接的,所以點云精度和POS系統(tǒng)關(guān)系緊密。

1.2 數(shù)據(jù)獲取

1.2.1 研究區(qū)概況

研究區(qū)位于新疆第八師150團(tuán)5連周邊的荒漠植被區(qū)(85°58′35′E～85°59′4′E,44°57′29′N～44°58′0′N),主要以旱生和沙生灌木為主,如梭梭、駱駝蓬、麻黃、堿蓬草和駝絨藜等。于2017年7月31日,用無人機(jī)搭載傳感器在荒漠植被區(qū)上空飛行采集數(shù)據(jù),飛行速度5 m/s,飛行高度60 m,共3條航線,每條航線掃描面積86 400 m2。

機(jī)載LiDAR系統(tǒng)可獲取目標(biāo)對象的豐富屬性信息,其數(shù)據(jù)形式稱為點云[15]。本文研究以WGS84為地理坐標(biāo)系統(tǒng),以UTM_Zone_44 N為投影坐標(biāo)系獲取點云數(shù)據(jù),以LAS格式存儲數(shù)據(jù)[16]。在獲取LiDAR數(shù)據(jù)過程中,由于外界環(huán)境的影響導(dǎo)致系統(tǒng)容易產(chǎn)生噪聲點云[17],這些噪聲點多為孤立點,表現(xiàn)為明顯高于地物或者低于地表的點云[18],所以本文研究在采用激光掃描儀配套軟件調(diào)整點云數(shù)據(jù)的掃描姿態(tài)、方位和粗差等處理后,再手動去除大量噪聲點云。

1.2.2 點云可視化

截取研究區(qū)部分點云作為本文原始數(shù)據(jù),其包含研究區(qū)主要地物類型,每個點的屬性主要包括點云三維坐標(biāo)、強度、回波次數(shù)、掃描角度、GPS時間等信息[19]。

為方便后期描述,本文將點云三維坐標(biāo)、回波強度和回波次數(shù)分別命名為[X,Y,Z,I,Re]。截取區(qū)域的覆蓋面積約為80 m×77 m,點云密度為128.224 pts/m2,點云三維坐標(biāo)范圍為{X=[20.042 5,100.697 5],Y=[80.507 5,157.875 0],Z=[275.290 0,284.039 2]},其中X和Y分別表示東和北方向的位置,Z表示點云高度,單位為m。回波強度范圍為{I=[5 767,65 535]},掃描視場角度范圍為{Angle=[-22,33]}。

圖1是以點云高度和回波強度屬性信息可視化原始激光雷達(dá)點云數(shù)據(jù)。

圖1 試驗區(qū)點云高度和回波強度圖Fig.1 Elevation and echo intensity maps of original point clouds in test area

不同激光掃描儀系統(tǒng),可以記錄不同的回波次數(shù),包括單次回波和多次回波。在多次回波中,接受第1個和最后一個回波信號的分別被稱為首次和末次回波,除了首次和末次回波外的回波信號為中間次數(shù)的回波[20]。

圖2a是以全部回波次數(shù)可視化試驗區(qū)點云,其中深棕色表示第1次回波點云,深藍(lán)表示第2次回波點云,綠色表示第3次回波點云。分析圖2a可知:

第1次回波對應(yīng)的點云個數(shù)為777 915,占點云總數(shù)98.49%,點云高度范圍為{Z=[275.29,284.039 2]}。

第2次回波獲取的點云個數(shù)為11 581,占點云總數(shù)1.49%,高度范圍為{Z=[275.35,280.17]}。

第3次回波獲取點云個數(shù)為339,占點云總數(shù)0.04%,點云高度范圍為{Z=[275.49,278.44]}。

第4次和第5次對應(yīng)的點云個數(shù)分別為20和4,高度范圍分別為{Z=[275.81,278.21]}和{Z=[275.64,278.11]}。

上述分析表明,第1次與第2次回波對應(yīng)點云的高度范圍相差很小,但第1次與第3、4、5次回波的點云高度范圍相差較大,且目視判別第3、4、5次回波對應(yīng)的點云,發(fā)現(xiàn)其均為噪聲點云。

對第2次回波點云的強度進(jìn)行可視化,結(jié)果如圖2b所示,圖中不同顏色代表不同強度值。從圖2b可知,第2次回波主要包括地面點云、電線的陰影點云和微量低矮植被點云。因為本文最后通過對比點云濾波結(jié)果分析檢測離群點云結(jié)果,即只保留地面點云,所以本文后期僅使用箱形圖方法分析第1次回波對應(yīng)的點云,再合并檢測結(jié)果和第2次回波點云,作為點云濾波處理的輸入數(shù)據(jù)集。

圖2 試驗區(qū)點云回波次數(shù)可視化Fig.2 Showed by point clouds echo times of test area

1.3 檢測離群點云

1.3.1 技術(shù)路線

利用Matlab R2017軟件編程提取點云三維坐標(biāo)、回波強度和回波次數(shù),在分析每次回波對應(yīng)點云的分布情況后,保留第1次和第2次回波對應(yīng)的點云。首先，采用箱形圖方法對第1次回波點云的高度屬性進(jìn)行離群分析,得到高度離群點云數(shù)據(jù)集D1；再在分離點云數(shù)據(jù)集D1的基礎(chǔ)上,對剩余點云進(jìn)行強度屬性離群分析,得到強度離群點云數(shù)據(jù)集D2、分離高度及強度離群點云后對應(yīng)的點云數(shù)據(jù)集D3；然后,合并數(shù)據(jù)集D3和第2次回波對應(yīng)的點云,得到點云數(shù)據(jù)集D4,并結(jié)合軟件ENVI 5.3對數(shù)據(jù)集D4進(jìn)行濾波；最后,與直接對原始點云濾波結(jié)果對比,其技術(shù)流程如圖3所示。

圖3 結(jié)合點云屬性和箱形圖檢測離群點云技術(shù)流程圖Fig.3 Flowchart of detecting of point clouds′ outliers based on combining point clouds′ attribute information and box-plot method

1.3.2 箱形圖原理

檢測離群點云的主要目的是識別不同地物類別,常用方法主要是基于統(tǒng)計、鄰近度、密度和聚類來實現(xiàn)檢測[21]，但是,在對單變量數(shù)據(jù)集進(jìn)行離群探測分析時,箱形圖因其簡單、直觀的優(yōu)點,成為最常選擇的方法之一[22]。箱形圖是描述最小值、第一四分位數(shù)Q1、中位數(shù)、第三四分位數(shù)Q3、最大值和異常值的一種統(tǒng)計方法,其中心位置為中位數(shù),箱子長度表示四分位數(shù)的間距IQR,離群值是根據(jù)異常值截距線和極端值截距線得到,其中異常值截距線和極端值截距線的定義如式(1)所示。

IQR=Q3-Q1，

B1=Q3+1.5×IQR,B2=Q1-1.5×IQR，

BE1=Q3+3×IQR,BE2=Q1-3×IQR。

(1)

式(1)中，B1和B2為異常值截距線對應(yīng)的值,BE1和BE2為極端值截距線對應(yīng)的值。

本文離群點云是由屬性大于B1和小于B2,或者大于BE1和小于BE2的點云組成。

2 結(jié)果與分析

2.1 箱形圖分析高度屬性

利用回波次數(shù)去除噪聲點云后,采用箱形圖對第1次回波對應(yīng)的點云高度進(jìn)行離群統(tǒng)計分析,得到如圖4所示的點云高度離群統(tǒng)計結(jié)果,其中縱坐標(biāo)表示點云高度值。

圖4 第1次回波點云對應(yīng)的高度箱形圖及高度離群點云可視化Fig.4 Box-plot and outlier point clouds of elevation values corresponding to the first echo point cloud data

圖4a為第1次回波次數(shù)對應(yīng)的高度箱形圖,經(jīng)過分析得出:點云高度Z為278.976 7是整個研究區(qū)高度臨界點,且將高度范圍在{Z=[278.976 7,284.039 2]}之間的點云歸為高度離群點,對高度離群點云進(jìn)行三維可視化得到灰度圖4b,令該數(shù)據(jù)集為D1,其對應(yīng)的回波強度范圍為{I=[5 767,65 535]}。將高度范圍在{Z=[275.290 0,278.976 7]}之間的點云歸為分離高度離群后的點云,回波強度范圍為{I=[7 143,58 369]},此時第一四分位數(shù)對應(yīng)的點云高度為276.628 2,即第1次回波中包含25%的點云高度的范圍為275.290 0～276.628 2;第三四分位數(shù)對應(yīng)的點云高度為277.567 9,第1次回波包含75%的點云高度的范圍為277.567 9～278.976 7。

2.2 箱形圖分析強度屬性

通過上述對點云高度屬性數(shù)據(jù)進(jìn)行箱形圖分析,可以去除大量電線、車輛和微量高植被地物點云數(shù)據(jù)集D1，之后在分離高度離群點云的基礎(chǔ)上,對剩余點云進(jìn)行強度數(shù)據(jù)離群分析,得到分離高度離群點云后對應(yīng)的強度離群點云數(shù)據(jù)集D2,再去除點云數(shù)據(jù)集D1和D2得到剩余點云數(shù)據(jù)為D3。

不同數(shù)據(jù)集可視化結(jié)果如圖5所示。

圖5 分離第1次回波高度離群后對應(yīng)的點云強度分析結(jié)果可視化Fig.5 Visualization of analyzing point clouds echo intensity after separating elevation values outlier point clouds from first echo

圖5a是分離高度離群點后點云強度箱形圖,縱坐標(biāo)表示點云回波強度值。由圖5a可知,此時強度臨界值是31 325,故將回波強度范圍{I=[7 143,31 325]}歸為分離高度離群后對應(yīng)的強度離群點云數(shù)據(jù)集D2,點云高度范圍為{Z=[275.490 5,278.972 0]};將回波強度范圍{I=[31 325,58 369]}歸為分離高度和強度離群后對應(yīng)的點云數(shù)據(jù)集D3,其點云個數(shù)為766 380,點云高度范圍為{Z=[275.290 0,278.976 8]}。另外，第一四分位數(shù)和第三四分位數(shù)分別對應(yīng)的點云強度值是41 571和48 408,即分離高度離群后對應(yīng)25%的點云強度的范圍為7 143～41 571,對應(yīng)75%的點云強度的范圍為48 408～58 369。

圖5b為強度離群點云數(shù)據(jù)集D2三維可視化結(jié)果,從圖5b可知，點云數(shù)據(jù)集D2主要包括大量植被點云、微量道路點云和微量車輛點云。

圖5c為分離高度和強度離群后的點云數(shù)據(jù)集D3的可視化結(jié)果,主要包含了研究區(qū)大量的地面點云。為了進(jìn)一步提高點云濾波精度,得到研究區(qū)全部地面點云,將點云數(shù)據(jù)集D3和第2次回波點云合并得到數(shù)據(jù)集D4。

2.3 點云濾波

在軟件ENVI 5.3 ENVI LiDAR模塊的支持下,通過反復(fù)建立三角網(wǎng)分別對點云數(shù)據(jù)集D4和原始點云進(jìn)行濾波,圖6是濾波后地面點云可視化結(jié)果。

圖6a是對數(shù)據(jù)集D4濾波得到的283 852個地面點云,點云高度范圍為{Z=[275.291 5,278.312 7]},回波強度范圍為{I=[8 322,58 369]},對應(yīng)濾波消耗的時間為14.1 s。

圖6b是對原始點云濾波獲取的281 301個地面點云,點云高度范圍為{Z=[275.294 5,278.363 0]},回波強度范圍為{I=[7 470,57 976]},對應(yīng)濾波消耗的時間為46.3 s。

對比兩次濾波的結(jié)果(圖6a、b)可知,基于點云回波次數(shù)、高度和強度屬性信息,采用箱形圖方法對原始激光雷達(dá)點云進(jìn)行濾波前處理,有效保留了地面點云2551個,濾波消耗時間降低了32.2 s，而且提高了點云濾波速度。

圖6 點云數(shù)據(jù)集D4和試驗區(qū)原始點云濾波結(jié)果可視化Fig.6 Filtering of point cloud datasets D4 and original point cloud data in test area respectively

為了分析對比兩次地面點云是否存在差異,本文將上述地面點云導(dǎo)入SPSS軟件中,對點云高度進(jìn)行統(tǒng)計分析,得到點云數(shù)據(jù)集D4與原始點云分別濾波后對應(yīng)地面點云高度中位數(shù)、眾數(shù)、標(biāo)準(zhǔn)差、方差、偏度相差和偏度標(biāo)準(zhǔn)誤差都為276.900、276.512、0.610、0.372、0.003和0.005,即對2種點云數(shù)據(jù)集進(jìn)行濾波得到的地面點云無明顯差異,因此,本文采用簡單易實現(xiàn)的箱形圖方法對點云屬性進(jìn)行濾波前處理是可行且有效的。

3 結(jié)論

本文在點云濾波前,結(jié)合點云屬性信息,利用箱形圖對原始激光雷達(dá)點云進(jìn)行濾波前預(yù)處理,通過對比處理前后的濾波結(jié)果,得到以下結(jié)論:

(1)結(jié)合點云回波次數(shù),可以有效剔除多次回波點云363個。

(2)結(jié)合箱形圖先后分析點云高度和點云強度的離群情況,且利用反復(fù)建立的三角網(wǎng)實現(xiàn)濾波,可獲取地面點云766 380個,濾波消耗時間減少了32.2 s，表明通過本文方法可提高點云濾波速度。

(3)本文采用簡單易實現(xiàn)的箱形圖方法對點云屬性進(jìn)行濾波前處理是可行且有效的。