Alpha-shape算法構(gòu)建棗樹點云三維模型

付昱興，李承明，朱 江，王寶龍，張 斌，付 威※

（1.海南大學機電工程學院，?？?570228；2.石河子大學機械電氣工程學院，石河子 832003）

0 引 言

紅棗是中國特色林果，含有豐富維生素C、多種微量鈣以及氨基酸等豐富的營養(yǎng)成分，具有極高的營養(yǎng)價值和醫(yī)藥價值[1-2]。目前，中國紅棗種植面積已占全世界紅棗種植面積的99%[3]，棗園亟待實現(xiàn)信息化管理。近年來，隨著計算機和信息技術(shù)的發(fā)展，果樹三維點云模型的研究正逐步成為國內(nèi)外農(nóng)業(yè)研究領(lǐng)域的重點和熱點[4]。利用三維重建技術(shù)構(gòu)建精確真實的棗樹三維形態(tài)，為棗樹智能化修剪提供可視化模型，能夠有效地提高棗園智能化作業(yè)水平[5]。

目前，針對農(nóng)林作物的三維重構(gòu)方法大致分為以下幾類：基于規(guī)則的建模、基于圖像的建模和基于點云的建模[6]。基于規(guī)則的建模通過研究植物實際生長過程的形態(tài)變化和拓撲結(jié)構(gòu)演變，建立植物數(shù)字模型庫。Prusinkiewicz等[7]基于龜圖解譯算法改進模擬植物生長的L系統(tǒng)（L-system），使其更好地演示現(xiàn)實中植物的生長過程；Linermann等[8]提出的植物生長系統(tǒng)X-frog通過控制生長參數(shù)，生成可視化三維植株。上述研究可以快速生成真實的三維樹木形態(tài)特征模型，但模型的植物學特征表現(xiàn)粗糙，不夠自然?；趫D像的建模方法根據(jù)圖像序列和視頻恢復(fù)三維信息，Guo等[9]通過多視角立體視覺，基于圖像序列實現(xiàn)樹模型重建；孫國祥等[10]基于相位相關(guān)原理完成了多視角圖像重建；楊長輝等[11]基于雙目相機在復(fù)雜背景下恢復(fù)重建柑橘枝干?；趫D像的建模方法較好實現(xiàn)了三維模型重建，但是需要從不同視覺采集圖像，且重建精度易受光照和復(fù)雜背景的影響?；邳c云的建模方法借助激光設(shè)備以及輔助設(shè)計進行點云獲取，在點云基礎(chǔ)上重建樹模型。李詩銳等[12]利用深度相機Kinect V2實時、準確地重建了三維場景；麥春艷等[13]和 Chattopadhyay等[14]基于多角度點云數(shù)據(jù)實現(xiàn)了對蘋果樹的建模和測量，提高了模型的真實感；Livny等[15]在點云基礎(chǔ)上使用加權(quán)生成樹構(gòu)建枝干結(jié)構(gòu)圖（Branch-Structure Graph，BSG），并通過施加軟約束實現(xiàn)重建，但該方法在枝條重疊時表現(xiàn)出不確定性；唐麗玉等[16]基于分離枝干與葉片的分治思想，根據(jù)枝干骨架點和冠層特征點拓撲結(jié)構(gòu)重建樹模型，取得了較好的重建效果；師翊等[17]在 Runions等[18]研究的基礎(chǔ)上引入角度約束改進空間殖民算法，基于實測點云重構(gòu)出高精度的樹模型。

本研究在粗精配準結(jié)合的思想基礎(chǔ)上，提出一種適合休眠期棗樹三維點云處理和重建的方法。采用微軟公司第二代深度相機 Kinect V2獲取點云數(shù)據(jù)，基于kd-tree-ICP算法與人工標記結(jié)合的方法提高配準精度，實現(xiàn)多視角的點云配準，獲得完整點云信息；通過Alpha-shape算法進行表面重構(gòu)，為棗樹智能修剪提供可視化模型和基礎(chǔ)數(shù)據(jù)。

1 材料與方法

1.1 試驗設(shè)備和信息采集

采用深度相機Kinect V2（微軟公司，美國）作為信息采集設(shè)備，基于飛行時間原理（Time of Flight，TOF）獲取點云數(shù)據(jù)。相比于V1版本，Kinect V2作為升級版[19]，抗干擾能力更強，精度更高，可達0.02 mm，其參數(shù)如表1所示。

表1 深度相機Kinect V2參數(shù)Table 1 Parameters of depth camera Kinect V2

試驗地點為新疆阿拉爾13團棗園（81°28'E，40°34'N），試驗時間為2019年11月份，選擇晴朗無風的天氣，在16：00—17：00室外光照較弱的時段，使用 Kinect V2在自然環(huán)境下對駿棗樹進行拍攝，傳感器與樹干中心的水平距離距離約為2 m，試驗場景如圖1a所示。試驗基于雙站信息采集策略，每180°采集1次數(shù)據(jù)，對樹進行正反視角拍攝（圖1b）。

圖1 數(shù)據(jù)采集Fig.1 Data acquisition

使用微軟公司提供的應(yīng)用程序開發(fā)工具 Kinect for Windows SDK以及應(yīng)用程序編程接口（Application Programming Interface，API）獲得彩色圖像坐標系、深度圖像坐標系和相機坐標系之間的坐標映射關(guān)系，融合果樹的彩色和深度信息，計算得到目標場景在相機坐標系中具有顏色的三維坐標信息。為減少傳感器系統(tǒng)噪聲與波動，將單個視角下的多幀點云數(shù)據(jù)疊加融合獲得彩色點云圖像（圖2）。本試驗采用Visual Studio 2015 集成化開發(fā)環(huán)境、OpenGL圖形標準庫、點云庫（Point Cloud Library，PCL）和MATLAB 2018b完成。

圖2 彩色原始點云圖像Fig.2 Original color point cloud image

1.2 點云數(shù)據(jù)預(yù)處理

彩色圖像與深度圖像經(jīng)過信息融合生成了大量點云，點云中夾雜了噪聲點，會對后續(xù)的點云配準工作造成較大影響，需要進行去噪處理。

噪聲主要分為背景點、地面噪聲以及樹周圍的離群散點。針對噪聲數(shù)據(jù)的不同特點，采用不同的去噪方法。背景點基于目標的三維坐標信息，根據(jù)感興趣區(qū)域（Region Of Interest，ROI）的深度信息d(x，y)去除果樹背景信息[20]。通過設(shè)定果樹的深度范圍，將取值范圍外的點云視為無效點去除。根據(jù)試驗條件，卷尺測得的棗樹最近端及最遠端與傳感器的有效距離位于 1 300～3 000 mm之間，根據(jù)式（1）得到去除背景后的點云（圖3a）。

式中D(x，y)為處理后的深度信息，mm；d(x，y)為原始深度信息，mm；x、y為圖像像素橫縱坐標。

針對圖3的地面噪聲和離群散點，分別采用直通濾波和統(tǒng)計濾波器（statistical outlier removal）去除[21]。直通濾波的具體做法為以30 mm為步長，在棗樹豎直高度方向建立直方圖統(tǒng)計點云數(shù)量。目標個體點云數(shù)量顯著少于地面區(qū)域點云數(shù)量，從地面大尺度噪聲向ROI過渡時出現(xiàn)點云數(shù)量驟減，基于此特點確定坐標范圍，分割地面與ROI。

圖3 背景去噪后點云Fig.3 Point cloud after background denoising

離群散點密度稀疏，圍繞枝干離散分布，統(tǒng)計濾波器基于k近鄰原理實現(xiàn)，首先查找每個點pi(xi,yi,zi)(i=1,2,…,n)的k個近鄰點pj(xj,yj,zj) (j=1,2,…,k)，并計算該點到k個近鄰點的平均距離（di，mm），然后計算全局平均距離平均值（μ，mm）和標準差（σ，mm），計算如式（2）所示：

式中k為近鄰點個數(shù)；n為點云中點的數(shù)量。

若平均距離di在標準范圍 [μ-γσ，μ+γσ]之外，則該點定義為離散點，將其從ROI中去除。其中，γ為縮放因子，根據(jù)點云密度確定，經(jīng)驗取值為[0,1]，k越大且γ值越小時，去噪效果越明顯，分別選取k為20、40和80，γ為0.2、0.5和1.0，分9組進行去噪試驗，部分試驗對比效果如圖4所示。由圖4a和圖4c可知，出現(xiàn)過度去噪，部分枝條點云被當作噪點去除，造成局部點云缺失；圖4b在保留枝干點云的前提下去除了較多的離散點，去噪效果最佳。

圖4 離散點去除結(jié)果對比Fig.4 Comparison of the removal results of outliers

1.3 點云配準

冠層分枝之間相互遮擋，單視角下的點云不能呈現(xiàn)分枝的完整形態(tài)。為保證數(shù)據(jù)的完整性，獲得棗樹完整的三維信息，采用雙站信息采集策略，背面點云與正面點云采用相同的去噪處理，并對兩站點云進行配準處理[22]。最近點迭代（Iterative Closest Point，ICP）算法廣泛應(yīng)用于圖像與點云的配準[23]，當兩站點云空間位置接近或者重疊區(qū)域較大時可以很好完成配準過程，但點云位置相距較遠時配準的穩(wěn)定性較差。因此，基于標靶球初配準，使兩站點云在空間位置接近，再利用 ICP算法搜索鄰近點間的對應(yīng)關(guān)系實現(xiàn)精確配準。

1.3.1 點云初配準

由于兩站的點云數(shù)據(jù)不在同一坐標系下，初配準過程是將原坐標系下的點（xs，ys，zs）通過坐標平移和旋轉(zhuǎn)變?yōu)槟繕俗鴺讼迪碌狞c（xg，yg，zg），如式（3）所示：

式中 [xg yg zg]T為目標坐標系；[xs ys zs]T為原坐標系；R為旋轉(zhuǎn)矩陣；λ為尺度因子；[txtytz]T為平移矩陣。

標靶球作為多站點云掃描數(shù)據(jù)拼接的公共同名點，其球心是點云拼接中的連接點和控制點，當兩站點云有3個公共點時可進行三維坐標轉(zhuǎn)換[24]。球心的提取需要滿足一定的精度，在使用最小二乘法擬合之前對球點云預(yù)處理[25]。首先，使用三維點云處理軟件MeshLab對標靶球點云分割提取，獲得的球點云（圖5a）；其次，針對點云中的噪點，使用統(tǒng)計濾波對其進行去噪處理，減少噪點帶來的擬合誤差，紅色為噪點（圖5b）；最后，采用下采樣方法精簡點云數(shù)量并保持點云形狀特征，提高球心的計算速度和估算精度（圖5c）。

圖5 球點云預(yù)處理Fig.5 Ball point cloud preprocessing

兩站點云的初始位姿，點云位置相距較遠（圖6a）。為了提高配準過程的精度，使用人工標記進行初配準，方法如下：

1）將3個標靶球放置在棗樹周圍，球直徑為12 cm，使用MeshLab軟件對3個標靶球點云分割提??；

2）對球點云預(yù)處理，使用最小二乘法擬合球心，將球心坐標（X，Y，Z）作為3個公共點（圖6b）；

3）根據(jù)姚吉利[24]的研究方法，將3個球心坐標代入式（4）～式（6），計算旋轉(zhuǎn)矩陣R、平移矩陣[txtytz]T以及尺度因子λ。

4）根據(jù)旋轉(zhuǎn)矩陣、平移矩陣和尺度因子實現(xiàn)原坐標系下點云的轉(zhuǎn)換，初配準結(jié)果如圖6c所示。

式中λ為尺度因子，其值為對應(yīng)球心連線邊長的平均值；a、b、c為反對稱矩陣中的獨立元素；XS21、YS21、ZS21分別為原坐標系中球心2與球心1的X、Y、Z坐標差值；XS31、ZS31分別為原坐標系中球心3與球心1的X、Z坐標差值；XG21、YG21、ZG21分別為目標坐標系中球心2與球心1的X、Y、Z坐標差值；XG31、YG31、ZG31分別為目標坐標系中球心3與球心1的X、Y、Z坐標差值；式中[txtytz]T為平移矩陣，[XG1YG1ZG1]T為目標坐標系球心1坐標；[XS1YS1ZS1]T為原坐標系球心1坐標。

1.3.2 點云精配準

ICP算法通過多次迭代搜索近鄰點逐漸逼近最優(yōu)值，為防止陷入局部最優(yōu)解和迭代時間過長，采用k維樹分割空間數(shù)據(jù)結(jié)構(gòu)精配準。點云J和點云K為2片待配準的點云，各有NJ和NK個點，ICP算法在迭代計算中搜索點云J和點云K中的最近點，并計算歐氏距離平方和的最小值，從而得到兩點集間的剛性變換值。v=[r，t]表示剛性變換值，其中r為最優(yōu)旋轉(zhuǎn)變換的單位特征矩陣，t表示平移矩陣，將該值應(yīng)用到點云J進行坐標變換。點云J中的每一個點，在點云K中搜索與其最近的點形成點對，并計算全部點對歐氏距離平方和的平均值f（v）如式（7）所示：

式中Ji為點云J中任一點；Ki為點云K中與Ji距離最近的點；Rr為剛性變換的旋轉(zhuǎn)矩陣。

其中，用單位四元數(shù)表示最優(yōu)旋轉(zhuǎn)矩陣r=[r0r1r2r3]，r0≥0，且r02+r12+r22+r32=1，旋轉(zhuǎn)矩陣Rr如式（8）所示：

經(jīng)過ICP算法多次迭代，f(v)在給定閾值下收斂為最小值，此時的變換矩陣為最優(yōu)解，如式（9）所示：

由圖6c可知，兩站點云經(jīng)過初配準，部分區(qū)域在空間位置上接近或重疊。計算點云的曲面法向量和曲率，對于點云J中的每一個非平面點Ji，在點云K中搜索與其曲率相近的點，通過對Ji設(shè)定歐式距離閾值（Dp，mm）確定搜索范圍，將閾值內(nèi)曲率相近點構(gòu)成配對點對；再使用kd-tree建立高維索引樹形數(shù)據(jù)結(jié)構(gòu)，對配對點對數(shù)據(jù)篩選重組，加快配準速度；最后，使用ICP算法求得最優(yōu)剛性變換值，將該值應(yīng)用到點云J實現(xiàn)魯棒的配準（圖6d）。

1.4 表面重構(gòu)

三維形狀由配準后的點云表示，不能很好體現(xiàn)棗樹原有的表面凹凸信息。因此，通過Alpha-shape算法對點云進行網(wǎng)格剖分，重建三維表面。算法本質(zhì)為滾球法，通過控制不同半徑的滾球，求取包絡(luò)，確定整個點集的表面[26-28]。Alpha-shape算法步驟如下：

1）配準的點云是一個點集Q={qj}，(j=1,2,…n)，從中任意選取一個點q1，與之距離＜2α的點構(gòu)成新的點集Q1，從點集Q1中任取一組點q2和q3，求出過點q2、q3和q1且半徑為α的球的球心o和o′；

2）遍歷點集Q1，依次求出其他點到球心o和o′的距離集合l和l′。如果l和l′中有一個集合的距離均＞α，則可以判斷出點q1、q2和q3是邊緣輪廓點，連接三點構(gòu)成邊界三角形；反之，則不是邊緣輪廓點，停止遍歷，執(zhí)行步驟3）；

3）選擇點集Q1中的下一組點按步驟1和2進行判斷，直到點集Q1中的所有點判斷結(jié)束輸出面片集合Δ，集合中暴露的三角面片形成局部凸包δSα(P)，計算如式（10）所示：

4）選擇點集Q中下一個點按步驟 1）～3）進行判斷，直到點集Q中的所有點判斷結(jié)束，多個局部凸包δSα(Q)重構(gòu)表面Sα(Q)。

圖6 兩站點云配準過程Fig.6 A pair of point cloud registration process

將精配準后的點云作為輸入，使用Alpha-shape算法，通過不斷調(diào)整半徑值α取值，將點云三角面片化，橢圓框線內(nèi)為局部放大圖，大量三角面片構(gòu)成了模型表面（圖7）。

圖7 樹模型重構(gòu)Fig.7 Reconstruction of tree model

2 結(jié)果與分析

2.1 點云配準誤差分析

為對點云配準算法精度定量分析，將配準點云對應(yīng)點之間的平均歐氏距離作為點云配準誤差的評估指標。本研究配準算法和 k維樹最近點迭代（k dimensional-tree-Iterative Closest Point，kd-tree-ICP）算法的配準性能比較（表2）。由表2可知，本算法在引入初配準后，整體點云在配準計算時間上明顯增加，但同時因為初配準的引入使得配準誤差大幅度降低，從2.06 cm降至0.75 cm，配準精度提升了63.6%。雖然本算法在計算時間上有所增加，但配準的精度得到了保證。為進一步驗證本算法的可靠性，分別使用本算法和 kd-tree-ICP算法對 6棵棗樹進行測試得到兩者的配準誤差對比（圖8），本算法的配準誤差更小，配準誤差均控制在1.0 cm以內(nèi)，平均配準誤差為0.76 cm。

表2 點云配準算法性能分析Table 2 Performance analysis of point cloud registration algorithm

圖8 本算法和k維樹最近點迭代算法的配準誤差對比Fig.8 Comparison of registration errors between proposed algorithm and kd-tree-ICP algorithm

2.2 重構(gòu)模型精度分析

為定量評價重建的各枝條模型的精度，使用鑫瑞 10分度游標卡尺SRUNV測量棗樹主干和主枝的直徑，測量位置為圖7中的標記點I～V（白色實線）。使用MeshLab軟件測量重構(gòu)模型對應(yīng)標記點的直徑，計算實測值與模型值之間的絕對誤差（表3）。誤差總體控制在7%以內(nèi)，重建模型的精度比較高，可以直接反映休眠期落葉棗樹枝干的重建精度。

2.3 建模試驗驗證

對 3棵棗樹的配準點云進行表面重構(gòu)，進一步對驗證建模試驗，建模結(jié)果如圖9所示。建模后的模型表面具有粗糙感，在外觀上更接近自然樹，觀察效果更加直觀。在圖9b和圖9d的矩形框內(nèi)，出現(xiàn)了不同程度的枝條缺失現(xiàn)象，主要由于光照影響導致部分枝條紅外線不能正常反射到傳感器，造成了局部點云缺失，枝條斷節(jié)。在圖9f矩形框的放大圖中，重建出了少量孔洞，其原因為曲面局部樣本的獲取主要基于歐式距離，通常偏向局部樣本點分布密集區(qū)域，當區(qū)域內(nèi)同時包含樣點分布稀疏區(qū)和樣點分布密集區(qū)，則在該處容易重建出孔洞。

表3 樹模型實測值與重構(gòu)值對比Table 3 Comparison of measured value and reconstructed value of tree model

圖9 建模驗證結(jié)果Fig.9 Modeling verification results

3 結(jié) 論

1）通過點云直方圖設(shè)定分割閾值和統(tǒng)計濾波器進行離群局外點判定，實現(xiàn)了感興趣區(qū)域與地面分離以及濾波去噪處理，從而獲得了準確的單株果樹點云數(shù)據(jù)。

2）提出了人工標記與改進的最近點迭代（Iterative Closest Point，ICP）算法結(jié)合的配準方法，試驗結(jié)果表明，通過人工標記初配準拉近兩站點云的空間位置，可以大幅度提高 ICP算法的點云配準精度。通過對多株棗樹進行測試，配準誤差均控制在1.0 cm以內(nèi)，平均配準誤差為0.76 cm，驗證了配準算法的穩(wěn)定性。在下一步的研究中，嘗試提取點云特征關(guān)鍵點，根據(jù)關(guān)鍵點特征值的空間映射關(guān)系實現(xiàn)兩站點云的自動初配準，解決耗時長的問題。

3）使用 Alpha-shape算法對配準點云進行面片化處理，實現(xiàn)曲面重構(gòu)。重構(gòu)模型枝干的相對誤差總體控制在7%以內(nèi)，平均誤差為2.2 mm，表明模型精度較高，可以直接反映休眠期棗樹枝干的重建模型表面具有粗糙感，真實感較強，在外觀上更加接近真實樹，可作為棗樹智能修剪的可視化基礎(chǔ)數(shù)據(jù)。