顧雯鈞，徐愛(ài)俊，尹建新

（1.浙江農(nóng)林大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院，浙江 杭州 311300；2.浙江農(nóng)林大學(xué) 浙江省林業(yè)智能監(jiān)測(cè)與信息技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311300；3.浙江農(nóng)林大學(xué) 林業(yè)感知技術(shù)與智能裝備國(guó)家林業(yè)與草原局重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311300）

立木胸徑和樹(shù)高是衡量立木生長(zhǎng)狀況的重要因子[1]，立木因子的獲取是森林資源調(diào)查中最重要的工作之一[2-3]。傳統(tǒng)的立木因子測(cè)量主要依靠人工或二維圖像完成。人工測(cè)量方法中，使用皮尺、測(cè)距儀等測(cè)樹(shù)工具[4]，工作量大且主觀誤差較大；電子經(jīng)緯儀、全站儀等精密設(shè)備[5-6]專業(yè)性強(qiáng)，但操作復(fù)雜。隨著計(jì)算機(jī)視覺(jué)技術(shù)的發(fā)展，基于二維圖像的測(cè)量方法[7-9]取得了一定的成果，但由于維度的限制，這類方法對(duì)于遮擋嚴(yán)重的場(chǎng)景測(cè)量精度不高。近幾年，三維重建方法被廣泛用于植物表形參數(shù)的測(cè)量[10-11]。其中地面激光掃描技術(shù)(terrestrial laser scanning, TLS)[12-13]測(cè)量精度高，但地面激光掃描儀價(jià)格昂貴，產(chǎn)生點(diǎn)云數(shù)據(jù)量大，對(duì)計(jì)算機(jī)的運(yùn)行速度、傳輸容量等有較高要求，使其在植株表型參數(shù)提取與分析方面的應(yīng)用受到了制約，難以滿足社會(huì)化的需求。運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)(structure from motion，SfM)[14]是一種通過(guò)特征點(diǎn)匹配建立不同視角的二維圖像的關(guān)系，從而恢復(fù)相機(jī)參數(shù)與目標(biāo)三維信息的算法。該算法為植株表型參數(shù)測(cè)量提供了一種新思路[15]。梁秀英等[16]使用戶外小車自動(dòng)拍攝田間玉米Zea mays植株，基于SfM 算法首先獲取玉米植株稀疏三維點(diǎn)云，然后再進(jìn)一步得到植株稠密點(diǎn)云，通過(guò)點(diǎn)云數(shù)據(jù)測(cè)得玉米株高、莖粗和葉面積等參數(shù)。HUI等[17]基于SfM算法，利用Visual SFM軟件對(duì)茄Solanum melongena、辣椒Capsicum annuum和黃瓜Cucumis sativus等3種蔬菜進(jìn)行了重建和表型解析，獲取了植物葉片數(shù)、葉長(zhǎng)、葉寬和葉面積數(shù)據(jù)。該方法除用于小型作物表型參數(shù)的測(cè)量外，在立木因子測(cè)量方面也有應(yīng)用。MILLER等[18]、徐慧丹等[19]圍繞單株立木拍攝重疊度約70%的圖片，并基于SfM算法獲取其三維點(diǎn)云模型，從而估計(jì)樹(shù)高、胸徑等參數(shù)。孫英偉等[20]在室外圍繞單株立木進(jìn)行拍攝，并基于SfM算法獲取立木三維點(diǎn)云，將點(diǎn)云最高和最低的差值作為樹(shù)高數(shù)據(jù)，對(duì)1.3 m處點(diǎn)云擬合圓提取胸徑值。這些方法雖然能夠有效替代傳統(tǒng)方法和TLS測(cè)量立木因子，但均是圍繞立木人為控制拍攝間隔獲取圖像，步驟繁瑣且精度難以保障，僅能測(cè)量單株立木，效率較低。

本研究提出基于運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)的多株立木因子測(cè)量方法，圍繞立木場(chǎng)景拍攝視頻，使用關(guān)鍵幀提取算法自動(dòng)提取三維重建所需的立木圖像，減少拍攝的復(fù)雜程度；基于SfM算法獲取立木三維點(diǎn)云后，進(jìn)行多株立木點(diǎn)云分割，實(shí)現(xiàn)立木單株點(diǎn)云的提取和樹(shù)高、胸徑的自動(dòng)測(cè)量，為多株立木因子的自動(dòng)測(cè)量提供了新思路。

1 材料與方法

1.1 立木視頻采集

以浙江農(nóng)林大學(xué)校園內(nèi)的立木為研究對(duì)象，在自然條件下使用智能手機(jī)環(huán)繞包含多株立木的場(chǎng)景拍攝一段視頻。拍攝時(shí)將手機(jī)移動(dòng)速度盡量保持均勻并減少晃動(dòng)，盡可能拍到立木的全部，并將相機(jī)鏡頭始終指向場(chǎng)景中心。由于每個(gè)場(chǎng)景中的立木數(shù)量、高度、間距均不相同，因此拍攝路線也不同。圖1為某場(chǎng)景的拍攝路線示意圖，場(chǎng)景中放置的藍(lán)色長(zhǎng)方體尺寸已知，有助于將立木三維點(diǎn)云的尺度恢復(fù)到與實(shí)際相同。

圖1 拍攝路線示意圖Figure 1 Schematic diagram of shooting route

1.2 立木三維點(diǎn)云獲取

1.2.1 視頻關(guān)鍵幀提取 視頻相比圖像能夠傳達(dá)更為豐富的信息，但是一個(gè)完整的視頻序列包含較多冗余數(shù)據(jù)，關(guān)鍵幀提取算法能夠從視頻中選擇滿足特定要求的幀，使得提取的關(guān)鍵幀序列能夠在滿足后續(xù)立木三維點(diǎn)云獲取所需基本信息的基礎(chǔ)上盡量減少冗余。將固定采樣法和相似度度量法相結(jié)合，對(duì)立木視頻進(jìn)行關(guān)鍵幀提取。先用固定采樣法提取視頻幀作為候選幀，即以視頻的第1幀作為候選幀的開(kāi)始，使用固定間隔為t幀的采樣算法At對(duì)包含L幀的視頻VL進(jìn)行候選幀提取，得到候選幀序列Fc，其中f1,f1+t,···,f1+nt為候選幀序列中的每幀，(1±nt)≤L。

式(1)中，將所有的候選幀提取完成后，以候選幀的第1幀作為關(guān)鍵幀的開(kāi)始，采用差異值哈希(dHash)算法[21]計(jì)算當(dāng)前關(guān)鍵幀與候選幀之間的相似度。該算法首先將2幅圖像分別縮小并灰度化，然后比較每行相鄰元素的大小，根據(jù)比較結(jié)果賦值1或0，得到2幅圖像的哈希值，最后利用漢明距離計(jì)算2幅圖像之間的相似度S。

式(2)中：Ii表示當(dāng)前關(guān)鍵幀，fj表示候選幀，dH(Ii)和dH(fj)分別表示2幅圖像的哈希值，D表示漢明距離。

計(jì)算出相似度S之后，設(shè)定相似度閾值S1、S2(分別取值為0.6和0.8)。當(dāng)S＞0.8時(shí)，將篩選出的關(guān)鍵幀序列進(jìn)行三維重建時(shí)立木三維模型的效果并沒(méi)有明顯提升，但是關(guān)鍵幀圖像大量增加導(dǎo)致效率降低；當(dāng)S＞0.6可以較好地排除視頻中突然出現(xiàn)的干擾(如車輛)。因此，若0.6≤S≤0.8，則將該候選幀加入到關(guān)鍵幀序列，并將該幀設(shè)置為當(dāng)前關(guān)鍵幀。然后繼續(xù)計(jì)算當(dāng)前關(guān)鍵幀與下一候選幀的相似度，直至篩選完所有候選幀，得到視頻關(guān)鍵幀序列Fe。

1.2.2 立木三維點(diǎn)云生成 在立木視頻關(guān)鍵幀序列(即有序圖像集)提取的基礎(chǔ)上，基于運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)算法生成立木三維點(diǎn)云，基本流程如圖2所示。①特征點(diǎn)提取與匹配。通過(guò)尺度不變特征變換算法(scaleinvariant feature transformation，SIFT)[22]對(duì)圖像進(jìn)行特征點(diǎn)提取，并計(jì)算每2幅圖像特征點(diǎn)描述子之間的歐幾里得距離來(lái)完成特征點(diǎn)匹配，得到所有圖像的特征匹配關(guān)系，其中n為關(guān)鍵幀圖像的數(shù)量，m中儲(chǔ)存了2幅圖像間的特征點(diǎn)匹配對(duì)，其中i為特征點(diǎn)匹配對(duì)的數(shù)量。②使用2幅圖像進(jìn)行重建初始化。首先將第1幅圖像與第2幅圖像中的匹配點(diǎn)對(duì)x和x′代入對(duì)極幾何約束條件[式(3)]求得基礎(chǔ)矩陣F，該矩陣代表了2幅圖像匹配點(diǎn)之間的空間幾何關(guān)系。

圖2 生成立木三維點(diǎn)云流程圖Figure 2 Flow chart of creating 3D point cloud of trees

再通過(guò)已知的相機(jī)內(nèi)參K和式(4)求得本質(zhì)矩陣E和，對(duì)E進(jìn)行奇異值分解后得到第2幅圖像相對(duì)第1幅圖像的相機(jī)旋轉(zhuǎn)矩陣R和平移向量T。

在此基礎(chǔ)上，可通過(guò)式(5)求得2個(gè)相機(jī)的投影矩陣P1和P2：

通過(guò)P1和P2建立二維匹配點(diǎn)x和x′及其對(duì)應(yīng)的三維空間點(diǎn)X的關(guān)系：

并將匹配點(diǎn)x和x′的齊次坐標(biāo) [u,v, 1]T及 [u′,v′, 1]T代入上式，得：

式(7)中：Pi1、Pi2、Pi3是Pi的行向量，對(duì)A進(jìn)行奇異值分解后得到二維圖像特征點(diǎn)的三維坐標(biāo)。

③加入新的圖像。根據(jù)新加入的一幅圖像與已重建的上一幅圖像上特征點(diǎn)的匹配關(guān)系，獲得已知的三維點(diǎn)X對(duì)應(yīng)新圖像上的二維點(diǎn)x，將對(duì)應(yīng)的二維點(diǎn)和三維點(diǎn)代入式(8)可求得投影矩陣Pi，然后重復(fù)使用步驟②中的方法獲得新圖像上其余未重建的匹配點(diǎn)的三維坐標(biāo)。

④利用光束平差(bundle adjustment，BA)[23]算法進(jìn)行集束調(diào)整，減小相機(jī)姿態(tài)R、T和重建三維點(diǎn)坐標(biāo)的誤差。將所有三維點(diǎn)通過(guò)重新映射的方式到得到其在圖像上的二維點(diǎn)坐標(biāo)，然后計(jì)算與其真實(shí)二維點(diǎn)之間的距離，并將該距離最小化以完成優(yōu)化。即：

式(9)中：xij是第i個(gè)特征點(diǎn)在第j幅圖像上的坐標(biāo)，Pj為投影矩陣，Xi為第i個(gè)特征點(diǎn)的三維坐標(biāo)。每重建一幅圖像進(jìn)行一次集束調(diào)整，當(dāng)所有圖像重建完成后，再次進(jìn)行集束調(diào)整，獲得優(yōu)化后的稀疏三維點(diǎn)云。

⑤三維點(diǎn)云稠密化。采用多視角立體集群算法(cluster multi-view stereo，CMVS)和多視角拼接算法(patch-based multi-view stereo，PMVS)對(duì)重建出來(lái)的稀疏點(diǎn)云進(jìn)行稠密化[24]。

1.3 立木點(diǎn)云數(shù)據(jù)處理

1.3.1 點(diǎn)云數(shù)據(jù)預(yù)處理 由于立木場(chǎng)景三維點(diǎn)云數(shù)據(jù)量大，夾雜了許多噪聲點(diǎn)，且點(diǎn)云模型與真實(shí)世界三維空間的方向、尺度上也存在不一致，因此需要對(duì)點(diǎn)云模型進(jìn)行下采樣、去噪、比例校正、坐標(biāo)校正、方向校正等預(yù)處理，以減少后續(xù)算法運(yùn)行時(shí)間。①采用三維體素柵格法對(duì)點(diǎn)云進(jìn)行下采樣。首先為輸入的點(diǎn)云數(shù)據(jù)創(chuàng)建由大量指定邊長(zhǎng)的微小立方體包圍的集合，即三維體素柵格；然后在每個(gè)體素內(nèi)計(jì)算所有點(diǎn)的重心，利用重心代替該體素內(nèi)的其余點(diǎn)，這樣所有的三維點(diǎn)云最終都由各個(gè)體素內(nèi)的重心表示。②去噪。使用直通濾波器對(duì)X、Y、Z閾值范圍進(jìn)行設(shè)定，劃分出感興趣的立木范圍區(qū)域，從而減少背景噪聲點(diǎn)。接著采用統(tǒng)計(jì)濾波方法去除離群噪聲點(diǎn)，計(jì)算直通濾波后的每個(gè)點(diǎn)pi，計(jì)算其與k(取50)個(gè)近鄰點(diǎn)的平均距離(pi)，假設(shè)得到的結(jié)果服從高斯分布，得平均距離均值(μ)和標(biāo)準(zhǔn)差(σ)，如果這些點(diǎn)的平均距離(pi)大于閾值Cσ(C為常量，取1)，則將被視為離群值去除。③坐標(biāo)校正。首先進(jìn)行點(diǎn)云模型與實(shí)際的比例校正，由放置物體實(shí)際尺寸hr和放置物體的模型尺寸hm得到坐標(biāo)縮放比例r=hr/hm，在進(jìn)行立木因子測(cè)量時(shí)可按比例將其換算到與實(shí)際相同的尺度。④方向校正。計(jì)算三維點(diǎn)云模型質(zhì)心點(diǎn)的三維坐標(biāo)，并基于主成分分析法(principal component analysis，PCA)[25]將點(diǎn)云模型進(jìn)行質(zhì)心化，求解點(diǎn)云模型協(xié)方差矩陣的特征值和特征向量，特征向量為立木三維點(diǎn)云模型的主方向，根據(jù)該主方向創(chuàng)建旋轉(zhuǎn)矩陣M，將原始點(diǎn)云Po以質(zhì)心點(diǎn)為中心轉(zhuǎn)換到主方向的位置，得到新的三維點(diǎn)云模型Pn=MPo。

1.3.2 多株立木點(diǎn)云分割 預(yù)處理后的三維點(diǎn)云模型中包含多株目標(biāo)立木及地面，為此需要進(jìn)行點(diǎn)云分割，分別提取單株立木三維點(diǎn)云。首先使用RANSAC算法將地面點(diǎn)云擬合成一個(gè)平面，然后提取屬于該平面內(nèi)的點(diǎn)云，從而將地面點(diǎn)云與立木點(diǎn)云進(jìn)行初步分割。

在初步分割后的點(diǎn)云中，同一立木的點(diǎn)云總是緊鄰的，而不同立木的點(diǎn)云之間存在間距，因此，根據(jù)空間點(diǎn)的歐氏距離采用條件歐幾里得聚類算法進(jìn)行多株立木三維點(diǎn)云的分割。對(duì)于三維點(diǎn)云中的某一點(diǎn)P，通過(guò)KD樹(shù)找到三維點(diǎn)云中離點(diǎn)P最近的k個(gè)點(diǎn)，這些點(diǎn)與P點(diǎn)的距離若小于閾值d則聚類到集合Q中。然后從集合Q中選取點(diǎn)P以外的點(diǎn)，重復(fù)上述過(guò)程，直到集合Q中點(diǎn)的數(shù)量不再增加，則聚類過(guò)程結(jié)束。最后將多株立木三維點(diǎn)云分割到不同的點(diǎn)云集合中，從而提取出單株立木三維點(diǎn)云。

1.4 立木因子測(cè)量

1.4.1 立木高度提取 將立木在z軸方向的最高點(diǎn)與最低點(diǎn)的差值默認(rèn)為其高度，如圖3所示。由于點(diǎn)云坐標(biāo)校正后單株立木的生長(zhǎng)方向均與z軸的正方向一致，且立木底部平面位于oxy平面，所以提取出單株立木點(diǎn)云后，遍歷所有點(diǎn)，找到單株立木z軸坐標(biāo)的最大值與最小值，兩者差值的絕對(duì)值即為單株立木的高度h(h=|Zmax-Zmin|)。

圖3 立木高度提取示意圖Figure 3 Schematic diagram of tree height extraction

1.4.2 立木胸徑提取 將立木位于距地面1.3 m處斷面直徑的最大值和最小值的平均視為其胸徑，如圖4所示。在單株立木的三維點(diǎn)云中，利用校正后的比例r求得1.3 m處點(diǎn)云坐標(biāo)：z=1.3/r。對(duì)該處做3 cm厚度的切片，然后再將點(diǎn)云切片投影至oxy平面，利用非線性最小二乘法進(jìn)行橢圓擬合求解橢圓一般方程式：Ax2+Bxy+Cy2+Dx+Ey+1=0。其中：A、B、C、D、E是待求解參數(shù)。求得橢圓一般方程式后，根據(jù)式(10)和式(11)求得橢圓的短軸a和長(zhǎng)軸b：

圖4 胸徑提取示意圖Figure 4 Schematic diagram of tree DBH extraction

然后將橢圓的短軸與長(zhǎng)軸的平均值作為立木的胸徑 (Φ)：Φ=(a+b)/2。

1.5 數(shù)據(jù)庫(kù)及工具

基于軟件平臺(tái)Visual Studio 2017和計(jì)算機(jī)視覺(jué)Opencv庫(kù)提取立木視頻關(guān)鍵幀圖像并基于運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)算法獲取立木三維點(diǎn)云；基于C++開(kāi)源點(diǎn)云庫(kù)(point cloud library，PCL)對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理及分割操作。

2 結(jié)果與分析

2.1 立木三維點(diǎn)云處理結(jié)果

立木視頻關(guān)鍵幀提取后，基于運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)算法獲取的立木場(chǎng)景三維點(diǎn)云如圖5A所示?？梢钥闯鲈撛既S點(diǎn)云不僅包含了4株立木，還包含了地面、背景及離群噪聲點(diǎn)。由圖5B所示：進(jìn)行點(diǎn)云下采樣之后點(diǎn)云數(shù)據(jù)減少了約3倍，這并不影響立木三維點(diǎn)云的整體形狀，且經(jīng)過(guò)預(yù)處理后的點(diǎn)云成功濾除了背景及大量離群噪聲點(diǎn)，點(diǎn)云模型的坐標(biāo)方向也得到了校正(該場(chǎng)景中種植立木的花壇長(zhǎng)度已知，無(wú)需在場(chǎng)景中放置用于比例校正的物體)。圖5C是將地面點(diǎn)云與立木點(diǎn)云進(jìn)行初步分割后的結(jié)果，圖5D是點(diǎn)云通過(guò)條件歐幾里得聚類算法分割后得到的4株立木。

圖5 立木三維點(diǎn)云處理結(jié)果Figure 5 3D point cloud processing results of trees

2.2 立木因子測(cè)量精度分析

由表1可見(jiàn)：經(jīng)過(guò)計(jì)算得出樹(shù)高測(cè)量的絕對(duì)誤差為0.03～0.24 m，測(cè)量值與實(shí)際值之間的相對(duì)誤差為0.63%～3.93%，平均相對(duì)誤差為1.96%。樹(shù)高測(cè)量值總體上比實(shí)際值低，主要是因?yàn)闃?shù)冠頂部的枝葉稀少，導(dǎo)致生成的三維點(diǎn)云在樹(shù)冠頂部較為稀疏，在點(diǎn)云處理階段被視為噪聲點(diǎn)去除，并且在一些枝葉較細(xì)或者高度較高的立木上，由于視頻拍攝不夠清晰導(dǎo)致特征點(diǎn)匹配數(shù)量不足，使得立木三維重建時(shí)存在點(diǎn)云部分缺失的現(xiàn)象。立木的胸徑測(cè)量的絕對(duì)誤差為0.20～0.85 cm，測(cè)量值與實(shí)際值之間的相對(duì)誤差為1.45%～6.92%，平均相對(duì)誤差為3.19%。測(cè)量誤差主要是由于拍攝時(shí)立木之間的遮擋，導(dǎo)致部分立木樹(shù)干重建效果不佳，并且在使用非線性最小二乘法進(jìn)行橢圓擬合時(shí)，由于利用的點(diǎn)云數(shù)量較少且分布不均勻，使得擬合不夠準(zhǔn)確。

表1 樹(shù)高和胸徑的測(cè)量值與實(shí)際值Table 1 Measured and actual values of tree height and DBH

圖6為樹(shù)高和胸徑測(cè)量值與實(shí)際值的線性回歸分析圖，使用本研究方法測(cè)量的立木樹(shù)高和胸徑與實(shí)際值之間均存在著顯著的線性相關(guān)關(guān)系，樹(shù)高相關(guān)系數(shù) (R2)為 0.987 9，均方根誤差 (RMSE)為 0.133 3 m；胸徑R2為 0.962 1，RMSE為 0.533 7 cm。

圖6 樹(shù)高和胸徑的測(cè)量值與實(shí)際值的回歸分析Figure 6 Regression analysis of measured and actual values of tree height and DBH

3 結(jié)論與討論

針對(duì)當(dāng)前基于三維點(diǎn)云的立木因子測(cè)量方法效率低或成本高的問(wèn)題，本研究提出了一種基于運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)的多株立木因子測(cè)量方法，使用普通智能手機(jī)實(shí)現(xiàn)了自然環(huán)境下多株立木樹(shù)高和胸徑的靈活、有效測(cè)量。本研究表明：使用本方法測(cè)量樹(shù)高的平均相對(duì)誤差為1.96%，相關(guān)系數(shù)為0.987 9，均方根誤差為 0.133 3 m；測(cè)量胸徑的平均相對(duì)誤差為 3.19%，相關(guān)系數(shù)為 0.962 1，均方根誤差為 0.533 7 cm。本方法測(cè)量樹(shù)高的平均相對(duì)誤差均低于文獻(xiàn)[26-27]，具有更高的樹(shù)高測(cè)量精度；測(cè)量胸徑的平均相對(duì)誤差低于文獻(xiàn)[26]，略高于文獻(xiàn)[27]，但是文獻(xiàn)[26-27]分別需要使用三維激光掃描儀和ToF手機(jī)，而本方法僅需使用普通智能手機(jī)，具有更好的普適性。

本方法具有較高的精度，但是測(cè)量胸徑和樹(shù)高仍有一定的誤差，可能原因?yàn)椋孩贅?shù)冠頂部的枝葉稀少，導(dǎo)致生成的三維點(diǎn)云在樹(shù)冠頂部較為稀疏，在點(diǎn)云處理階段被視為噪聲點(diǎn)去除，導(dǎo)致樹(shù)高測(cè)量值總體上比實(shí)際值低。②在一些枝葉較細(xì)或者高度較高的立木上，由于視頻拍攝不夠清晰導(dǎo)致特征點(diǎn)匹配數(shù)量不足，使得立木三維重建時(shí)點(diǎn)云部分缺失。③采用非線性最小二乘法擬合橢圓時(shí)由于樹(shù)干點(diǎn)云的數(shù)量較少或分布不均導(dǎo)致立木胸徑計(jì)算結(jié)果不理想。因此，可嘗試采用高分辨率無(wú)人機(jī)等設(shè)備進(jìn)行樹(shù)冠頂部視頻的拍攝，以獲得立木清晰且全面的圖像，同時(shí)進(jìn)一步研究立木圖像特征點(diǎn)提取與匹配算法，增加樹(shù)干與樹(shù)冠頂部枝葉的特征點(diǎn)匹配對(duì)數(shù)量，以提高重建立木的三維點(diǎn)云質(zhì)量。