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

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

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

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

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

1 材料與方法

1.1 無(wú)人機(jī)平臺(tái)

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

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

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

在數(shù)據(jù)獲取過(guò)程中,因?yàn)榧す鈷呙鑳xVUX-1和Survey+2是固定連接的,所以點(diǎn)云精度和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日,用無(wú)人機(jī)搭載傳感器在荒漠植被區(qū)上空飛行采集數(shù)據(jù),飛行速度5 m/s,飛行高度60 m,共3條航線,每條航線掃描面積86 400 m2。

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

1.2.2 點(diǎn)云可視化

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

為方便后期描述,本文將點(diǎn)云三維坐標(biāo)、回波強(qiáng)度和回波次數(shù)分別命名為[X,Y,Z,I,Re]。截取區(qū)域的覆蓋面積約為80 m×77 m,點(diǎn)云密度為128.224 pts/m2,點(diǎn)云三維坐標(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表示點(diǎn)云高度,單位為m?；夭◤?qiáng)度范圍為{I=[5 767,65 535]},掃描視場(chǎng)角度范圍為{Angle=[-22,33]}。

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

圖1 試驗(yàn)區(qū)點(diǎn)云高度和回波強(qiáng)度圖Fig.1 Elevation and echo intensity maps of original point clouds in test area

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

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

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

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

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

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

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

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

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

1.3 檢測(cè)離群點(diǎn)云

1.3.1 技術(shù)路線

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

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

1.3.2 箱形圖原理

檢測(cè)離群點(diǎn)云的主要目的是識(shí)別不同地物類別,常用方法主要是基于統(tǒng)計(jì)、鄰近度、密度和聚類來(lái)實(shí)現(xiàn)檢測(cè)[21]，但是,在對(duì)單變量數(shù)據(jù)集進(jìn)行離群探測(cè)分析時(shí),箱形圖因其簡(jiǎn)單、直觀的優(yōu)點(diǎn),成為最常選擇的方法之一[22]。箱形圖是描述最小值、第一四分位數(shù)Q1、中位數(shù)、第三四分位數(shù)Q3、最大值和異常值的一種統(tǒng)計(jì)方法,其中心位置為中位數(shù),箱子長(zhǎng)度表示四分位數(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為異常值截距線對(duì)應(yīng)的值,BE1和BE2為極端值截距線對(duì)應(yīng)的值。

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

2 結(jié)果與分析

2.1 箱形圖分析高度屬性

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

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

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

2.2 箱形圖分析強(qiáng)度屬性

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

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

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

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

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

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

2.3 點(diǎn)云濾波

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

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

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

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

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

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

3 結(jié)論

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

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

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

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