顧雯鈞，徐愛俊，尹建新

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

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

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

1 材料與方法

1.1 立木視頻采集

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

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

1.2 立木三維點云獲取

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

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

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

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

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

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

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

在此基礎上，可通過式(5)求得2個相機的投影矩陣P1和P2：

通過P1和P2建立二維匹配點x和x′及其對應的三維空間點X的關系：

并將匹配點x和x′的齊次坐標 [u,v, 1]T及 [u′,v′, 1]T代入上式，得：

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

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

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

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

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

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

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

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

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

1.4 立木因子測量

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

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

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

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

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

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

基于軟件平臺Visual Studio 2017和計算機視覺Opencv庫提取立木視頻關鍵幀圖像并基于運動恢復結構算法獲取立木三維點云；基于C++開源點云庫(point cloud library，PCL)對點云數(shù)據(jù)進行預處理及分割操作。

2 結果與分析

2.1 立木三維點云處理結果

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

圖5 立木三維點云處理結果Figure 5 3D point cloud processing results of trees

2.2 立木因子測量精度分析

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

表1 樹高和胸徑的測量值與實際值Table 1 Measured and actual values of tree height and DBH

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

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

3 結論與討論

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

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