劉茂華,韓梓威,陳一鳴,劉正軍,韓顏順
(1. 沈陽建筑大學(xué) 交通工程學(xué)院, 遼寧 沈陽 110168; 2. 中國測繪科學(xué)研究院, 北京 100089)
樹種分類是生態(tài)環(huán)境、林業(yè)測量、遙感等眾多行業(yè)和領(lǐng)域的重點研究課題,因為樹種識別對生態(tài)系統(tǒng)評估、生物多樣性監(jiān)測和森林資源利用起著至關(guān)重要的作用[1-2]。傳統(tǒng)的樹種分類方法通常是利用高光譜遙感技術(shù),通過樹木的光譜信息區(qū)分不同的樹種。然而,高光譜數(shù)據(jù)冗余度大,且存在著“同譜異物”和“同物異譜”的現(xiàn)象,忽略了樹木的三維結(jié)構(gòu)信息[3]。激光雷達(dá)(Light Detection And Ranging, LiDAR)是一種先進(jìn)的主動遙感技術(shù),它能夠高精度地快速獲取地表目標(biāo)的高度信息和三維結(jié)構(gòu)信息,且具有抗干擾能力強、低空探測性能好等優(yōu)點[4-5]。隨著無人機技術(shù)的快速發(fā)展,無人機LiDAR為快速、準(zhǔn)確的植被精細(xì)分類提供了有力的技術(shù)支撐[6]。
近年來,深度學(xué)習(xí)技術(shù)在三維數(shù)據(jù)上取得了很好的發(fā)展,根據(jù)輸入網(wǎng)絡(luò)的表達(dá)方式可以概括為以下三種:①基于體素的方法一般是將三維模型劃分為三維網(wǎng)格后再進(jìn)行3D卷積,代表模型有VoxNet、3D ShapeNets、O-CNN等[7-9];②基于多視圖的方法是先將三維模型通過投影等方法得到二維圖像,再利用圖像領(lǐng)域的方法進(jìn)行處理,代表模型有MVCNN、FPNN、DeepPano等[10-12];③基于點的方法是網(wǎng)絡(luò)能夠適應(yīng)三維數(shù)據(jù)的特點,直接對點云進(jìn)行處理,代表模型有PointNet、PointNet++、PointCNN等[13-15]。
由于深度學(xué)習(xí)在特征提取上的優(yōu)勢,近年來有學(xué)者將該技術(shù)運用到基于LiDAR數(shù)據(jù)的樹種識別。Guan等[16]基于體素實現(xiàn)了濾波和單木分割,然后將樹的點密度用波形表示,再利用深度玻爾茲曼機(Deep Boltzmann Machine, DBM)實現(xiàn)了對城市樹種的分類;Zou等[17]利用柵格化的方法在樣本空間中劃分體素,將每個體素格網(wǎng)中的點云數(shù)累積,投影到包含樹木輪廓的二維圖像上,每30°重復(fù)投影,并用一個深度信念網(wǎng)絡(luò)(Deep Belief Network, DBN)來生成高級特征,從而進(jìn)行樹種分類;Hamraz等[18]將點云轉(zhuǎn)化為數(shù)字表面模型(Digital Surface Model, DSM)和2D視圖,再利用深度卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Networks, CNN)對針葉樹和落葉樹分類。
然而,這些研究都是先將點云進(jìn)行規(guī)則化處理,例如轉(zhuǎn)換成體素或二維圖像,而沒有直接利用三維數(shù)據(jù)進(jìn)行特征提取或分類。體素劃分不僅會呈指數(shù)級地增加計算成本,還會受到分辨率的限制,存在局部信息丟失的現(xiàn)象。多視圖投影的方法不能充分利用點云的三維結(jié)構(gòu)信息,在點分類和場景理解等3D任務(wù)上具有較大的局限性。因此,本文基于森林樣區(qū)內(nèi)獲取的無人機LiDAR點云數(shù)據(jù),提出了一種新的樹種分類方法,該方法將無序的點云數(shù)據(jù)作為神經(jīng)網(wǎng)絡(luò)的輸入,直接利用其三維信息進(jìn)行特征提取,實現(xiàn)白樺樹和落葉松兩個樹種的分類。
研究區(qū)域為中國河北省承德市圍場滿族蒙古族自治縣內(nèi)的塞罕壩國家森林公園(42°24′N,117°19′E)的一部分,海拔高度為1 500~2 067 m,植被覆蓋度約為80%。研究區(qū)內(nèi)包括白樺和落葉松兩類樹種共計約1 790棵,兩類樹主要以一條較寬的土路為分隔,其中白樺樹主要分布在路的西北部,約870棵;落葉松分布在路的東南部,約920棵。
實驗數(shù)據(jù)由大疆無人機搭載的Riegl公司miniVUX 1U激光雷達(dá)掃描儀于2019年8月24日獲取,航線間距40 m,航高50 m,飛行速度5 m/s,激光雷達(dá)掃描角90°~270°,掃描線速度100 m/s?;谠摕o人機LiDAR系統(tǒng)獲得的原始數(shù)據(jù)如圖1所示,兩類樹種點云樣例數(shù)據(jù)如圖2所示,樣區(qū)內(nèi)各樹種的結(jié)構(gòu)特征見表1。
圖1 無人機LiDAR原始數(shù)據(jù)Fig.1 Airborne LiDAR raw data
(a) 白樺(a) Birch (b) 落葉松(b) Larch圖2 基于無人機LiDAR系統(tǒng)的點云樣本示例Fig.2 Example of point cloud sample based on airborne LiDAR system
表1 樣本樹木結(jié)構(gòu)特征
對于獲取的原始無人機LiDAR點云數(shù)據(jù),首先進(jìn)行預(yù)處理,包括去噪和地面濾波。其次,對預(yù)處理后的森林點云數(shù)據(jù)進(jìn)行單木分割,獲取單木點云。再次,對單木點云添加標(biāo)簽,制作數(shù)據(jù)集。訓(xùn)練集用于訓(xùn)練神經(jīng)網(wǎng)絡(luò)模型,測試集用于對得到的模型測試樹種分類精度。最后,對可能的分類精度影響因素進(jìn)行分析,并做對比實驗?;跈C載LiDAR數(shù)據(jù)的三維深度學(xué)習(xí)樹種分類工作流程如圖3所示。
圖3 樹種分類工作流程Fig.3 Work flow of tree species classification
1.3.1 去噪與地面濾波
由于森林環(huán)境的復(fù)雜性,獲取機載LiDAR數(shù)據(jù)的過程中會產(chǎn)生粗差,因此必須先對原始點云數(shù)據(jù)去噪。對每一個點搜索相同個數(shù)的鄰域點,經(jīng)多次實驗發(fā)現(xiàn),當(dāng)鄰近點個數(shù)為10時,去噪效果最好。然后計算該點到鄰域點的距離平均值Dmean及其中值m和標(biāo)準(zhǔn)差σ。計算最大距離Dmax:
Dmax=m+Kσ
(1)
其中,K為標(biāo)準(zhǔn)差倍數(shù),設(shè)置為5。若Dmean大于Dmax,則認(rèn)為該點為噪點,將其去除。
為消除地形影響,需要先對去噪后的點云數(shù)據(jù)進(jìn)行地面濾波,以獲得更準(zhǔn)確的單木分割結(jié)果。本文利用改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波算法[19]對地面點與非地面點進(jìn)行分離。首先對LiDAR點云以1 m的尺寸劃分網(wǎng)格,取每個網(wǎng)格的最低點進(jìn)行光柵化,再利用形態(tài)學(xué)開操作選取潛在地面種子點Gpotential。然后,利用平移平面擬合法濾除Gpotential中的非地面點,進(jìn)而獲得準(zhǔn)確的地面種子點Gseeds。最后,建立TIN模型,執(zhí)行先向下再向上的迭代加密,得到地面點Gresult。地面點和非地面點的分類結(jié)果如圖4所示,其中藍(lán)色點為地面點Gresult,白色點是樹木點,刪除地面點即完成濾波。
圖4 地面點與非地面點的分類結(jié)果Fig.4 Classification results of ground points and non-ground points
1.3.2 單木分割
為了制作深度學(xué)習(xí)的數(shù)據(jù)集,向每棵樹添加代表樹種的標(biāo)簽,需要對原始的森林?jǐn)?shù)據(jù)進(jìn)行單木分割。單木提取存在一定的錯誤率,例如欠分割、多樹粘連等情況,所以使用單一方法分割得到的高質(zhì)量單木數(shù)量較少。因不同分割方法得到的單木三維形狀有所差異,可作為數(shù)據(jù)增強的一種方法,故本實驗分別采取改進(jìn)的分水嶺分割算法[20]和基于點云距離的分割方法[21]對原始的森林點云數(shù)據(jù)進(jìn)行分割,以篩選出數(shù)量足夠多的高質(zhì)量單木點云。改進(jìn)的分水嶺分割算法首先從原始的激光脈沖生成數(shù)字表面模型和數(shù)字高程模型,進(jìn)而得到冠層高度模型;然后,通過在具有可變窗口大小的冠層最大值模型中搜索局部最大值來檢測樹梢,其窗口大小由冠部大小和樹高度之間的回歸曲線預(yù)測間隔的下限確定,該步驟可有效減少對樹頂點的判別錯誤;最后,利用Marker-controlled分水嶺算法完成單木分割,如圖5所示?;邳c云距離的分割將點按順序從高到低進(jìn)行分類,間距大于指定閾值的點排除在目標(biāo)樹之外,間距小于閾值的點根據(jù)最小間距規(guī)則進(jìn)行分類,分割結(jié)果如圖6所示。通過這兩種算法,共分割出單木3 598棵,經(jīng)人工篩選,刪除過分割、欠分割以及粘連嚴(yán)重的樹,保留1 200棵高質(zhì)量樹作為實驗數(shù)據(jù)。
圖5 基于改進(jìn)的分水嶺分割的單木分割結(jié)果Fig.5 Results of individual tree segmentation based on improved watershed segmentation
圖6 基于點云距離的單木分割結(jié)果Fig.6 Results of individual tree segmentation based on point cloud distance
1.3.3 樣本集準(zhǔn)備
深度神經(jīng)網(wǎng)絡(luò)需要輸入樣本具有相同的點云數(shù),故將所有單木均勻采樣2 048個點,并將每棵樹零均值歸一化至單位球內(nèi),以解決點的平移不變性。向每棵樹添加標(biāo)簽,其中白樺樹記為0,落葉松記為1。為了避免因種植區(qū)域相同而導(dǎo)致分類器可能過于適合研究地點的情況,從不同的空間中劃分?jǐn)?shù)據(jù)。將整個數(shù)據(jù)集按照8 ∶2的比例隨機分成訓(xùn)練集和測試集,其中訓(xùn)練數(shù)據(jù)取自研究區(qū)西側(cè),測試數(shù)據(jù)選取研究區(qū)的東側(cè)樹。最后,將包含(x,y,z)坐標(biāo)的原始點云數(shù)據(jù)、標(biāo)簽值和歸一化數(shù)據(jù)轉(zhuǎn)換成HDF5格式。詳見表2。
表2 訓(xùn)練集與測試集
網(wǎng)絡(luò)模型分為特征映射模塊、對稱函數(shù)模塊和分類器模塊,如圖7所示。對每一個單木點云取固定值N個點,每個點有三個維度的特征(x,y,z),即點的三維坐標(biāo)。在特征映射模塊中,輸入的是N×3的矩陣,經(jīng)過三層權(quán)重共享的多層感知器(Multi-Layer Perceptron, MLP)將每個點映射到1 024維的高維空間,得到N×1 024的特征矩陣。在對稱函數(shù)模塊,通過對稱函數(shù)max pooling得到1×1 024維的全局特征。最后,在分類器模塊,以全局特征作為輸入,通過兩個全連接層和一個softmax分類器得到兩個樹種的分類概率。
圖7 神經(jīng)網(wǎng)絡(luò)模型Fig.7 Neural network model
1)特征映射模塊:多層感知器是一種前向結(jié)構(gòu)的人工神經(jīng)網(wǎng)絡(luò),層與層之間全連接,在輸入與輸出的向量之間能夠進(jìn)行非線性映射。將輸入層向量記為X,則隱藏層的輸出為f(W1X+b1),其中W1為權(quán)重(也叫連接系數(shù)),b1為偏置,f是激活函數(shù)。在該網(wǎng)絡(luò)模型中,感知器的每一層都使用了ReLU激活函數(shù)和批處理歸一化,輸出大小依次為64、128、1 024。
2)對稱函數(shù)模塊:點云是一個無序點的集合,而點的順序不會影響這個集合本身。如果把點云表示為一個N行D列的二維矩陣(其中N代表點的個數(shù),D代表點的維度),對該矩陣做行變換則一共有N!種變換方式,而這N!種置換代表的是同一點集,因此網(wǎng)絡(luò)需要擁有對點云的置換不變性。對稱函數(shù)能夠很好地解決這個問題,例如對某一點的集合X(x1,x2,…,xn),取該集合的最大值、最小值、平均值都與點的順序無關(guān)。但是,直接對每個點做對稱操作,例如取最大值,會得到整個點集的最遠(yuǎn)點值;取平均值,會得到點集的重心,這樣會損失點集有意義幾何信息。因此,在網(wǎng)絡(luò)結(jié)構(gòu)中,首先利用MLP對每個點做相同的高維映射,通過高維空間中信息的冗余性避免對稱操作導(dǎo)致的信息丟失。該結(jié)構(gòu)實質(zhì)上是h、g、Υ的函數(shù)組合:
f(x1,x2,…,xn)=Υg(h(x1),h(x2),…,h(xn))
(2)
其中,h代表高維映射,g代表對稱函數(shù),Υ代表MLP整合信息。
3)分類器模塊:對于得到的全局特征向量,將其輸入兩個全連接層(Fully Connected layers, FC)中,并在最后一個層(輸出維度為256)中使用0.7的dropout防止過擬合。在softmax分類損失中加入一個正則化損失(權(quán)重為0.001),使矩陣接近正交。
實驗使用NVIDIA Tesla V100-PCIE-16 GB顯卡,在Ubuntu 18.4系統(tǒng)下搭配cuda 10.0框架和cuDNN 7.4.1加速庫訓(xùn)練。編譯語言使用Python3.6,深度學(xué)習(xí)框架為Tensorflow-gpu-1.13.1。模型使用Adam優(yōu)化器,初始學(xué)習(xí)率為0.001,最小學(xué)習(xí)率設(shè)定為0.000 01,通過指數(shù)衰減的方式實現(xiàn)學(xué)習(xí)率衰減,衰減率為0.7,動量為0.9,共訓(xùn)練200個周期,每20個周期學(xué)習(xí)率除以2,批處理歸一化的衰減率從0.5開始,逐漸增加到0.99,訓(xùn)練時長18 min。為了更好地提高模型的泛化能力,在訓(xùn)練過程中,通過沿上軸隨機旋轉(zhuǎn)物體來動態(tài)地增強點云,并通過均值為0和標(biāo)準(zhǔn)偏差為0.02的高斯噪聲來抖動每個點的位置。
為了探究單木分割效果對分類結(jié)果的影響,同時使本文方法和實驗結(jié)果更具實際意義,本文又對數(shù)據(jù)集的結(jié)構(gòu)做了改變。保持原數(shù)據(jù)集單木個數(shù)不變,隨機抽取10%、20%、30%、40%、50%的樹木樣本,并加入同等數(shù)量的非完整分割單木。如圖8所示,替換的樣本有過分割和欠分割的單木,非完整分割的程度及替換次序為隨機抽取,未經(jīng)人工干預(yù)。
(a) 白樺(a) Birch
為了證明本實驗方法的優(yōu)越性,與以下三種方法進(jìn)行比較。
方法1:參考Lin等[22]提出的方法,計算點云數(shù)據(jù)的9個點分布特征參數(shù)、13個冠內(nèi)結(jié)構(gòu)參數(shù)、11個樹外結(jié)構(gòu)參數(shù),共計33個特征參數(shù)輸入一個三層的MLP分類器中進(jìn)行樹種分類。
方法2:參考Guan等[16]提出的方法,如圖9所示,取預(yù)定義的波形維數(shù)n為50,對分割得到的單木進(jìn)行n個剖面的垂直剖分,將每個剖面的點統(tǒng)計值歸一化至0~1內(nèi),以生成波形圖。將得到的波形圖輸入一個兩層的DBM模型中生成高級特征,進(jìn)行樹種分類。
(a) 白樺樹樣本及其波形圖(a) Sample of white birch tree and its waveform
方法3:參考Zou等[17]提出的方法,在樣本空間中對單木做柵格化劃分,計數(shù)每個網(wǎng)格塊中的點數(shù),沿y軸在xz平面上積累網(wǎng)格,得到類似灰度圖像的圖片。對單木點云以相同的旋轉(zhuǎn)度數(shù)重復(fù)投影,得到的二維圖像作為樹木的低級特征輸入一個三層的DBN,再使用softmax分類器對同一棵樹的所有投影圖像的分類結(jié)果進(jìn)行投票。參與投票的每幅圖像使用相同的權(quán)重,并通過贏得大部分選票的類別結(jié)果來判定樹種。圖10為每旋轉(zhuǎn)36°后的柵格化示意,藍(lán)色虛線框內(nèi)為投影結(jié)果。實驗中分類最優(yōu)結(jié)果產(chǎn)生在旋轉(zhuǎn)角度為10°時。
圖10 單木的柵格化與投影Fig.10 Rasterization and projection of individual tree
實驗分類結(jié)果如表3所示,樹種分類的總體準(zhǔn)確率(Overall Accuracy, OA)為86.7%,kappa系數(shù)為0.73;生產(chǎn)者精度(Producer′s Accuracy,PA)高于85.5%,用戶精度(User Accuracy,UA)高于85.0%。圖11顯示了數(shù)據(jù)集中非完整分割單木所占比例對分類結(jié)果的影響。由圖可以看出,當(dāng)少量樣本質(zhì)量變低時(如占比10%),分類準(zhǔn)確率并沒有受到較大影響,說明該網(wǎng)絡(luò)的魯棒性較好。當(dāng)?shù)唾|(zhì)量樣本變多時(如占比30%以上),網(wǎng)絡(luò)的分類精度會迅速下降,說明此時數(shù)據(jù)集的質(zhì)量對網(wǎng)絡(luò)的學(xué)習(xí)產(chǎn)生了較大的負(fù)面影響。當(dāng)?shù)唾|(zhì)量樣本占總樣本的50%左右時,網(wǎng)絡(luò)的分類準(zhǔn)確率也降至50%附近,此時對于一個二分類問題來說,相當(dāng)于網(wǎng)絡(luò)能在樣本中學(xué)到的東西是非常少的,近乎隨機猜測。
表3 樹種分類結(jié)果的混淆矩陣
圖11 數(shù)據(jù)集質(zhì)量對分類準(zhǔn)確率的影響Fig.11 Influence of data set quality on classification accuracy
由此可以得出,造成分類錯誤的原因主要有以下兩點:一是白樺和落葉松兩個樹種在結(jié)構(gòu)形態(tài)上具有較高的相似性;二是實驗區(qū)內(nèi)的樹木比較密集,單木分割的效果會直接影響最后的分類結(jié)果。
本文方法與1.5節(jié)所述三種方法的分類結(jié)果如表4所示。其中方法1的分類精度最低,總體精度為81.7%,kappa系數(shù)為0.63,該方法雖然統(tǒng)計了樹木的33個不同的結(jié)構(gòu)參數(shù),但有限參數(shù)加淺層學(xué)習(xí)的分類方法仍限制了其分類精度。
表4 四種方法的分類精度比較結(jié)果
方法2的總體準(zhǔn)確度為84.2%,kappa系數(shù)為0.68,優(yōu)于方法1,這是因為該方法通過波形數(shù)據(jù)更好地表示了不同樹種獨特的幾何結(jié)構(gòu),且使用DBM模型生成了高級特征,減少了人工解譯的錯誤。
方法3的總體精度為85.6%,kappa系數(shù)為0.71,優(yōu)于前兩個方法。相較于方法2對每棵樹生成唯一的波形表示,每10°重復(fù)投影對樹木信息的保留更加全面,且大大增加了樣本量,有利于深度學(xué)習(xí)模型的學(xué)習(xí),因此得到了更高的分類精度。
實驗表明,本文方法優(yōu)于以上三種方法,其原因是:①與將三維數(shù)據(jù)轉(zhuǎn)化成二維表達(dá)形式的方法相比,基于點的深層神經(jīng)網(wǎng)絡(luò)能夠有效減少數(shù)據(jù)轉(zhuǎn)化過程中的信息損失,很好地保留樹木的三維結(jié)構(gòu)信息;②本文提出的網(wǎng)絡(luò)模型能夠有效地提取LiDAR數(shù)據(jù)的高維特征,有利于提高樹種分類精度。
為了避免點云數(shù)據(jù)幾何信息的損失,本文在對稱操作之前,通過參數(shù)共享的MLP將每個點映射到一個高維空間。為了探究點的特征維度對分類結(jié)果的影響,對使用128、256、512、1 024、2 048維特征的分類準(zhǔn)確度進(jìn)行了比較(如圖12所示)。結(jié)果表明,隨著特征維度的增加,總體分類準(zhǔn)確度提高,點云的幾何信息得到了更好的保留。當(dāng)特征維度達(dá)到約1 000時,分類性能達(dá)到最好。當(dāng)特征維度超過1 000以后,分類性能無明顯提升。
圖12 特征維度對分類準(zhǔn)確率的影響Fig.12 Influence of feature dimensions on classification accuracy
本研究對分割得到的單木數(shù)據(jù)分別進(jìn)行五次均勻下采樣,如圖13所示,采樣點個數(shù)從128增加至2 048的過程中,分類準(zhǔn)確度提升了13.1%。當(dāng)每棵樹的采樣點數(shù)量超過2 000時,網(wǎng)絡(luò)的分類性能趨于飽和。點密度的增加促使分類效果大幅提升的原因主要是白樺和落葉松兩個樹種相似度較高,增加點密度能夠保留更多的幾何結(jié)構(gòu)信息,使得深度神經(jīng)網(wǎng)絡(luò)能夠更好地學(xué)習(xí)。而當(dāng)點密度增加到一定數(shù)值時,網(wǎng)絡(luò)能學(xué)到的信息趨于飽和,過大的點密度只會增加信息的冗余程度而很難提升分類性能。
圖13 密度對分類準(zhǔn)確率的影響Fig.13 Effect of dot density on classification accuracy
本研究提出了一種新的基于三維深度學(xué)習(xí)的樹種分類方法。該方法對預(yù)處理后的點云數(shù)據(jù)做單木分割,將得到的單木點云作為模型的輸入,利用由特征映射模塊、對稱函數(shù)模塊和分類器模塊組成的深度神經(jīng)網(wǎng)絡(luò)自動提取樹木的高維特征并實現(xiàn)樹種分類。通過實驗得到以下結(jié)論:
1)本研究在樹種分類任務(wù)上取得了較好的效果,在包含白樺樹和落葉松兩個樹種的數(shù)據(jù)集上實現(xiàn)了86.7%的總體精度和0.73的kappa系數(shù)。
2)實驗結(jié)果表明,相比于特征參數(shù)法、波形圖法和柵格化投影法,本文方法在特征提取上存在優(yōu)勢,有助于提高樹種分類精度。
3)隨著點的特征維度和點密度的增加,分類準(zhǔn)確率不斷提高,但存在提高的臨界值。在該數(shù)據(jù)集中,對每棵樹取2 048個點,每個點取1 024維特征時,模型的分類效果最好。