Alpha-shape算法構建棗樹點云三維模型

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

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

0 引 言

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

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

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

1 材料與方法

1.1 試驗設備和信息采集

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

表1 深度相機Kinect V2參數Table 1 Parameters of depth camera Kinect V2

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

圖1 數據采集Fig.1 Data acquisition

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

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

1.2 點云數據預處理

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

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

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

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

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

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

式中k為近鄰點個數；n為點云中點的數量。

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

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

1.3 點云配準

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

1.3.1 點云初配準

由于兩站的點云數據不在同一坐標系下，初配準過程是將原坐標系下的點（xs，ys，zs）通過坐標平移和旋轉變為目標坐標系下的點（xg，yg，zg），如式（3）所示：

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

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

圖5 球點云預處理Fig.5 Ball point cloud preprocessing

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

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

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

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

4）根據旋轉矩陣、平移矩陣和尺度因子實現原坐標系下點云的轉換，初配準結果如圖6c所示。

式中λ為尺度因子，其值為對應球心連線邊長的平均值；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算法通過多次迭代搜索近鄰點逐漸逼近最優值，為防止陷入局部最優解和迭代時間過長，采用k維樹分割空間數據結構精配準。點云J和點云K為2片待配準的點云，各有NJ和NK個點，ICP算法在迭代計算中搜索點云J和點云K中的最近點，并計算歐氏距離平方和的最小值，從而得到兩點集間的剛性變換值。v=[r，t]表示剛性變換值，其中r為最優旋轉變換的單位特征矩陣，t表示平移矩陣，將該值應用到點云J進行坐標變換。點云J中的每一個點，在點云K中搜索與其最近的點形成點對，并計算全部點對歐氏距離平方和的平均值f（v）如式（7）所示：

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

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

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

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

1.4 表面重構

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

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

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

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

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

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

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

圖7 樹模型重構Fig.7 Reconstruction of tree model

2 結果與分析

2.1 點云配準誤差分析

為對點云配準算法精度定量分析，將配準點云對應點之間的平均歐氏距離作為點云配準誤差的評估指標。本研究配準算法和 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以內，平均配準誤差為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 重構模型精度分析

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

2.3 建模試驗驗證

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

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

圖9 建模驗證結果Fig.9 Modeling verification results

3 結 論

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

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

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