基于偏移補償模型的極值點聚類苗帶識別算法研究

肖菲菲, 蔣蘋, 胡文武, 廖榮華, 張丹慧, 金生

(湖南農業大學機電工程學院, 長沙 410128)

目前，精準農業是農業智能化發展的研究熱點之一，其中田間自動導航技術是完成精準農業的關鍵環節[1]。農業機械的自動導航技術現已被廣泛應用于噴藥、施肥和除草等農業作業過程[2-3]，該項技術不僅可以降低勞動強度和農業生產成本，還可提高農機作業效率。而在農業機械的自動導航過程中，農作物的苗帶識別和定位是核心技術之一，在很大程度上影響農業機械的導航效果和作業的精準度[4-7]。苗帶是指根據農作物種植方式，在同一方向上農作物連續性的集合，在作物生長初期，苗帶為條狀；在作物生長后期，作物封行后，苗帶為塊狀。苗帶識別是指通過某種方法區分農作物不同生長狀態下的苗帶區域與相鄰苗帶間的行間間隙所在位置。通過苗帶識別，在可行走區域規劃行走路徑，使機械行走過程中避免碾壓苗帶，實現田間管理機械自主低傷苗作業，對促進精準農業發展有重要意義[8]。

當農業機械行駛在環境復雜的農田中，為了避免障礙物和壓傷農作物，需要通過傳感器進行感知識別，當前常用于苗帶檢測識別的傳感器有視覺傳感器和激光雷達等[9-14]。王姍姍等[15]通過視覺傳感器檢測水稻秧苗，采用基于特征點鄰域 Hough變換算法完成水稻苗中心線提取，獲得平均識別準確率為92%，其識別率較好，但響應速度慢，只適用于速度低于1.5 m·s-1。張博立等[16]利用機器視覺提取玉米苗帶中心線，采用動態網格與分區域聚類算法相結合進行識別，獲得平均玉米苗帶識別準確率為93.4%，有效地處理了雜草對識別效果的影響，但由于圖像相關處理問題，導致中心線檢測失敗,且目前僅適用于識別生長初期玉米作物行。姜國權等[17]通過粒子群優化算法獲取麥田作物行的特征點，利用最小二乘法進行線性擬合，運行速度相對傳統Hough算法有很大提高，但并不適用于處理雜草較多的圖像。Hu等[18]利用二維伺服平臺，通過激光雷達掃描獲得水稻苗帶點云信息，設計基于斜率虛化均值迭代聚類算法，采用斜率虛化處理點云數據，通過區域迭代聚類獲得苗帶特征，利用最小二乘法進行苗帶軌跡線性化擬合，經試驗獲得橫向最大誤差為45 mm，中位點偏差為7 mm，但未考慮激光掃描偏移問題，且存在重復運行算法時出現聚類失敗的情況。

由以上研究可知，視覺傳感器成本較低，可實時獲取作物行相對位置，因此廣泛用于作物苗帶識別，但傳感器易受外界環境因素的影響，且數據處理相對繁雜，部分系統并不能滿足真實的農田環境。而激光雷達具有易操作、分辨率高、抗干擾性強等特點，更宜用于農田環境，但目前大多研究集中于室內環境或已知結構的室外環境，對于水田的作物識別研究相對較少，且并未考慮掃描偏移對識別結果的影響。

綜上，本文提出了基于偏移補償模型的極值點聚類苗帶識別算法，重點研究激光掃描偏移對識別結果的影響。首次提出了苗帶識別的相關定義，根據獲得的苗帶點云線特征，采用極值檢測算法獲取聚類要素點，根據苗帶特征，分析激光測距儀的掃描特點，建立偏移補償模型，對獲取的聚類要素點進行偏移矯正，最后通過最小二乘法進行線性擬合，獲得苗帶軌跡信息。提高了軌跡重合精度，增強了算法運行效率，獲得的苗帶軌跡參數可為田間管理機械低傷苗直線行駛提供理論依據。

1 材料與方法

1.1 系統構成及工作原理

導航定位控制系統由植保機行駛控制單元(electronic control unit)和苗帶識別系統構成，其中苗帶識別系統包含激光雷達、舵機調平單元與慣性測量模塊，如圖1所示。

注：1—施肥斗; 2—發動機; 3—苗帶識別系統； 4—植保機行車控制單元；5—藥箱；6—噴桿。

激光雷達安裝于舵機平臺上，慣性測量模塊采集姿態數據，舵機調平單元通過采集慣性測量模塊的數據，控制舵機進行俯仰、翻滾、旋轉動作，進而實現激光雷達底座的俯仰角、偏航角、橫滾角可控調節，并通過CAN接口與植保機ECU單元進行數據交換。ECU單元采集激光雷達掃描的點云數據與舵機平臺角度數據后，通過處理獲得水稻苗帶信息，并識別可行駛區間，控制田間管理機械行駛于苗帶行間間隙中，減少作業過程中對水稻的損害。

1.2 苗帶識別系統硬件與數據采集軟件

水稻苗帶識別系統的基本結構(圖2)由三部分構成，分別為激光雷達、二維舵機調平單元和支撐平臺。本研究選取日本HOKUYO的UTM-30LX二維激光雷達，測量范圍為100～30 000 mm，100～10 000 mm的測距誤差為±30 mm，掃描角度范圍-135°～135°，掃描周期40 Hz，角度分辨率為0.25°。舵機調平單元采用達盛舵機RDS3115，實現X、Y軸的水平調節控制，在履帶底盤底部安裝旋轉調節舵機，實現X、Y軸的旋轉調節，完成激光雷達俯仰角、橫滾角、航向角的調節控制。為實現舵機角度初始化校準和激光雷達中位調整，采用核心處理器STM32F103配合MPU6050六軸姿態傳感器，實現舵機角度調節。

注：1—控制單元; 2—支撐平臺; 3—旋轉調節舵機； 4—激光雷達；5—俯仰調節舵機；6—翻滾調節舵機。

激光雷達通過軟件UrgBenri Standard V1.8.1進行數據采集，數據存儲格式為ubh，使用Lenovo CPU1.60GHz在Win10環境下運用Matlab R2016b進行數據分析與算法研究。

1.3 點云數據獲取與苗帶軌跡提取分析

1.3.1點云數據獲取為確定激光束發出位置，確保激光雷達采集數據的準確性，本研究在數據采集前通過最小二乘平差法求得真實距離與激光雷達所測距離之間的數學關系[19-20]對測量儀進行標定，以消除激光雷達自身掃描誤差。苗帶數據采集于2020年7月湖南農科院試驗田中，此時水稻處于移栽后15 d，當天天氣預報風速為2 km·h-1，環境溫度32 ℃，目測水稻葉片處于靜止狀態。采用卷尺測量水稻高度、行株距和激光雷達距離水田平面的高度等參數，采用5組連續15株水稻苗信息求均值的方法，獲得水稻高度和行距株距。通過測試計算，得到水稻苗與水面高度平均值為360 mm，平均行距為220 mm，激光雷達掃描中心點到水田水平面的垂直距離為960 mm。在實際稻田環境數據采集過程中，二維舵機調平單元的俯仰掃描初始角度為9°，步進角度0.5°，終止角度17°，激光雷達掃描角度范圍為-16°～16°，掃描寬度對應5行苗帶，以此試驗條件獲取苗帶點云數據。

1.3.2坐標轉化激光雷達測距顯示參數主要為當前激光雷達掃描的起止角、步長、掃描線數，導出數據為距離與角度的極坐標數據，系統參數還包括當前激光雷達的俯仰角度數據。圖3為苗帶識別系統工作時的相關參數，以水田水平面為參考，以激光雷達激光旋轉中心為原點，垂直方向為Z軸，激光雷達初始光軸與機具直行重合方向為Y軸，垂直前進方向為X軸，激光雷達俯仰角度θ，激光測距長度為L，對應激光線掃描線角度為γ，Y軸對應激光雷達中心位置偏移初始角度為γ0。為方便數據分析，對激光極坐標數據進行三維坐標轉換，基于上述坐標轉換模型獲得被測點相對激光雷達當前位置的相對直角坐標表達式(1)。

圖3 掃描系統相關參數

(1)

1.3.3苗帶軌跡獲取方法為獲得最終苗帶軌跡，本研究系統將激光雷達和慣性測量模塊采集的數據融合，獲得被測目標點的位置坐標，運用Matlab R2016b進行數據分析與算法研究。根據農業機械作業環境和激光雷達點云數據特征，確定本文數據處理流程為數據預處理、點云極值點提取、苗帶特征點聚類和軌跡獲取與評價四個部分。數據預處理主要包括點云數據格式轉換、數據去噪和數據坐標轉換等步驟；點云極值點提取主要實現苗帶特征點的提取；偏移補償均值聚類主要進行苗帶特征點聚類；軌跡獲取與評價主要完成聚類點的線性擬合，并對獲取的軌跡線進行有效性評價。

1.4 基于偏移補償模型的極值點聚類

1.4.1苗帶模型建立由于田間管理機械直線行駛僅需要獲取以水平面為參考的平面坐標，因此，對三維數據平面化，獲得基于激光雷達為原點的XY坐標定位數據。

(2)

為了方便計算苗帶位置，獲得苗帶函數，根據激光雷達掃描測距方法建立理想苗帶模型，如圖4所示，長條框為苗帶，激光測距儀距離植株距離b，苗帶寬度d1,相鄰苗帶間距為d2。為驗證模型的正確性，使用5個紙箱進行模擬試驗，中心間距500 mm，UrgBenri Standard V1.8.1采集圖像。由圖像對比可知 ，建立的模型與模擬實驗獲得的點云線相對應，圖像在X軸方向都具有周期性，且為具有規律性的波浪峰，模型結果成立。

圖4 數據采集模型

從數據采集模型(圖4)可知，紅色線為激光雷達采集的點云數據線，圖像在X方向具有周期性，周期為d1+d2，且圖像關于Y軸對稱，根據模型圖和已知參數計算模型函數，函數共分為三段，計算得出函數關系式(3)。

(3)

1.4.2建立偏移補償模型根據建立的模型可知，由于苗帶的物理遮擋，導致部分位置的點云數據無法正確采集。在掃描中間苗帶時(圖5)，掃描線經過Q點后，右側面的點云數據由于遮擋無法采集，使得在QN線上的點云數據產生了偏移，所以需要建立偏移補償模型進行矯正。在激光掃描范圍內，以原點為中心左右兩側最高峰為大波峰，最低峰稱為小波峰，Y小于小波峰以下的數據相對較為準確，偏移量較小，可通過聚類算法直接獲得苗帶聚類點；Y在小波峰和大波峰之間的數據偏移較為嚴重，需對獲得的苗帶聚類點進行偏移矯正?；谝陨蠁栴}，建立偏移補償模型，如圖5所示，在縱向小波峰到大波峰范圍(GA≤Y≤NA)內，獲得偏移量ΔX的矯正范圍。

圖5 偏移補償模型

在EM線上，通過苗帶模型函數，在點E、M處分別得到如下公式。

(4)

(5)

(6)

在QN線上，通過苗帶模型函數，在點H、N處分別得到如下公式。

t×tanγ1=GA=MB

(7)

(8)

(9)

則苗帶間隙最大偏移量計算如下。

(10)

因此，偏移量范圍為0≤ΔX,基于此建立的偏移模型在同一波中，取低于小波峰的點云偏移量趨向0，在小波峰與大波峰之間的點云偏移量為ΔX，獲得矯正模型(圖6)，表達如式(11)所示。

圖6 矯正模型

(11)

1.4.3聚類要素點的獲取根據建立的模型(圖4)可知，每條點云數據線在X方向都存在規律性的波浪峰，而波峰波谷所在位置與苗帶和苗帶行間間隙相關，然而波峰波谷位置又對應函數中的極值所在位置，因此，可通過計算函數的極值所在位置得到與苗帶位置相關信息。本研究采用極大值檢測算法進行極值檢測，尋找苗帶位置，由式(12)確定。

s=find{diff[sign(diff(w))]<0}+1

(12)

式中，diff指差分運算，sign為符號函數，find表示查找點云數據線中滿足條件的數據位置[21]。

當采樣點滿足式(12)時，其點云數據在鄰域范圍內呈先增加后減小的趨勢，即局部極大值點，獲得波峰位置；將原數據的縱坐標值取反，再利用函數查找波谷位置。因農田環境復雜，水稻葉片具有下垂交疊等特征，且此特征也可作為判別苗帶位置相關因素，為避免信息缺失，所以閾值范圍取最小值到最大值，使用Matlab 繪制多角度點云數據線的波峰波谷圖，進而獲得苗帶聚類特征信息。

1.4.4改進K-means均值聚類均值聚類算法的基本思想是以空間中k個點為中心聚類，對最接近他們的對象進行歸類，通過迭代的方法不斷調整各分類中心位置，直到收斂[22]。聚類算法中聚類數量與初始聚類中心的選擇會直接影響聚類結果，通過改進K-means聚類算法提高苗帶點云數據的聚類效果。通過矯正函數獲得最佳聚類個數K,具體流程如圖7。初始中心的選擇通過密度公式(13)確定，在r∈(0,d1+d2)范圍內，定義對應運算區域的聚類點位置X，計算xi內的點密度ρ(xi),當密度達到最大時，在對應范圍內采用歐式距離度量方式選取距離xi最遠的點作為第(i+1)類初始聚類中心xi+1，并且除去xi、xi+1，以此類推，重復上述步驟再獲得下一個聚類中心，不斷迭代，直至獲得第K個聚類中心結束。

圖7 獲取聚類點數流程

ρ(xi)={x∈X|dist(xi,x)≤r}

(13)

1.4.5苗帶軌跡處理為便于獲取苗帶軌跡，獲得更多用于線性擬合的聚類點，同時降低數據采樣維度，采用縱向等區間劃分增加數據獲取量，獲得更為準確的軌跡參數，對Y軸向范圍進行n等分區間聚類。通過分層聚類后，獲得分段范圍內的聚類點。

(14)

式中，k為單次聚類數量，(x,y)為聚類點對應的坐標位置,n為Y軸等分處理的次數，從而獲得nk個聚類點。則對應苗帶行間間隙點為(xv,yv),矩陣為O(n)。

(15)

(16)

(17)

(18)

基于偏移模型進行矯正獲得聚類矯正后的矩陣T。

(19)

從而獲得對應Y軸向的s組聚類點矩陣M(s)。

(20)

對上述矩陣坐標點進行線性擬合獲取對應苗帶行間間隙的軌跡方程。對聚類點坐標矩陣通過最小二乘法進行線性擬合，獲得Ak和Bk的值，從而獲取苗帶對應的理論軌跡。

(21)

(22)

對應苗帶間隙軌跡表達式如下。

y=Akx+Bk

(23)

1.5 苗帶軌跡評價

參考文獻[23]的評價指標，在獲得理論苗帶軌跡后，與對應水稻苗帶真實軌跡進行對比，采用測量軌跡與真實軌跡平行度和中心偏差距離兩個指標評價算法的平行度(式24)和中心偏差(式25)。

Err_max=|xk-x′k|max,yk∈(a～b)

(24)

Err_mid=|xk-x′k|max,yk=(a+b)/2

(25)

式中，xk為對應第k行x的真實值，x′k為對應第k行x的測量值。

2 結果與分析

2.1 測量儀標定誤差分析

通過5組不同距離的數據，對激光測距儀進行標定，標定后實際距離與激光雷達測量結果見表1?？梢钥闯?，隨著激光雷達所測距離由0.645增加到14.389 m,數據標定前與標定后的誤差也逐漸增加，但標定后的誤差明顯低于標定前的。當被測物體距激光雷達14.389 m時，標定后的誤差為0.011 m，相對標定前的誤差0.037 m明顯減少。說明測距儀使用前有必要進行標定。同樣，完成二維舵機平臺與激光雷達的聯合標定。

表1 雷達標定前后測量值誤差

2.2 試驗數據模型構建

農田環境中，作物通常采用彼此平行的種植方式[24]，采集試驗區域對應有效測量范圍的水平投影數據。以激光雷達為坐標原點，為消除側偏導致的方位差，激光雷達前視方向與中間列位置基本平行，為滿足前視行駛需要，設定前視水平投影距離為3 030 mm，前視水平投影探測有效范圍為2 000～3 030 mm，有效長度為1 030 mm，左右水平投影探測有效寬度為-500～+500 mm。測試環境的真實苗帶軌跡如式(26)。

(26)

則測試環境的可行駛苗帶行間間隙軌跡如式(27)。

(27)

2.3 苗帶點云數據分析

2.3.1苗帶三維點云顯示及分析通過Matlab R2016b對激光測距儀采集的點云數據進行預處理，繪制苗帶軌跡三維立體圖，如圖8所示。

圖8 MATLAB轉換后的點云數據

根據田間采集的實況與三維點云圖外觀對比分析可知，對應的苗帶信息清晰有規律，且相鄰苗帶行間間隙的初始軌跡鮮明。箭頭標明了苗帶行間間隙軌跡所在位置，但受葉片傾斜折彎和交叉下垂以及激光雷達掃描方式的影響，行間間隙寬度較窄，且隨著前視距離的增大間隙軌跡不易識別。多角度掃描軌跡線在靠近X軸有五個突出的波谷，為激光雷達掃描的第一行水稻苗所呈現出的形狀，且每根點云數據線都存在這種規律性的波浪峰狀特征。

2.3.2點云極值分析通過對每根點云線進行極值檢測，獲得波峰波谷數據點位置，獲取聚類要素點，如圖9所示。

圖9 波峰波谷數據點

注：圖中不同顏色的圈圈代表不同的分類，*代表此類型的聚類點。

表2 最終聚類中心

2.3.4苗帶軌跡提取通過聚類后進而獲得水平矩陣R(n)。

(28)

(29)

則對應苗帶行間間隙點為(xv,yv),矩陣為O(n)。

(30)

(31)

(32)

對應獲得水平聚類點矩陣Q。

根據實際采集情況，基于偏移模型計算獲得偏移量范圍為0≤ΔX≤174 mm，獲得聚類矯正后的矩陣T。

(33)

通過平面投影處理，獲得如圖11的分段均值聚類分布圖，其中黑點為聚類中心點，紅點為矯正點。

注：圖中不同顏色的圈圈代表不同的分類，*代表此類型的聚類點。

獲取Y軸向對應的4組聚類點矩陣M1～M4。

(34)

對應上述縱向4組苗帶行間間隙聚類點，在Matlab下進行線性擬合處理，獲得軌跡矩陣如式(35)，對應軌跡如圖12所示。

圖12 苗帶行間間隙軌跡對比

(35)

2.4 軌跡結果分析

在y∈(2 000,3 000) 范圍內，計算獲得軌跡與實際軌對應橫向最大偏差Err_max和中位點偏差Err_mid。

Err_max=[16111614]

(36)

Err_mid=[6461]

(37)

表4分別統計了采用兩種不同算法對苗帶點云數據處理后的X軸最大橫向偏差、中位點最大偏差和平均耗時。

表4 算法檢測偏差及檢測時間對比

從表5中數據可得出，在前視距離2 000～3 000 mm范圍內,通過斜率虛化聚類算法對點云數據進行處理，獲得軌跡與實際苗帶軌跡X軸橫向最大偏差為45 mm, 中位點最大偏差7 mm，平均耗時1.64 s。本文采用基于偏移補償模型的極值點聚類苗帶識別算法獲得X軸橫向最大偏差控制在16 mm以內，中位點最大偏差為6 mm，運算耗時1.50 s，相比斜率虛化聚類算法X軸向最大偏減小了29 mm，中心點最大偏差減小1 mm，耗時減少0.14 s，檢測效果明顯提高，所檢測苗帶行間間隙軌跡與實際軌跡具有較好的重合度，滿足水田管理機械免壓苗行駛的需求(輪寬10 cm)。

3 討論

苗帶軌跡的識別是實現田間自動導航技術關鍵。現有對作物的檢測大多針對旱地作物,如蔬菜、玉米、小麥等，對于水田的作物識別研究相對較少。本文首次提出了苗帶識別的概念，利用二維激光雷達與舵機調平單元相結合，以獲得三維掃描功能完成苗帶點云數據獲取，進行水稻苗帶識別研究，避免了視覺傳感器易受外界環境影響，數據處理的算法復雜度較高等問題，重點研究激光掃描偏移對識別結果的影響，為解決激光雷達掃描偏移和提高軌跡識別精度等問題，提出了基于偏移補償模型的極值點聚類苗帶識別算法。根據獲得的點云線特征，利用極值檢測算法獲得聚類要素點，分析激光雷達的苗帶掃描特征，建立了偏移補償模型，結合改進的K-means聚類算法進行偏移矯正聚類，最后通過最小二乘進行線性擬合，從而實現苗帶軌跡可行駛區域的識別。

通過對前視距離為2 000～3 000 mm范圍內的苗帶軌跡進行識別，獲得最大橫向偏差低于16 mm,中位值偏差控制在6 mm以內。相對斜率虛化聚類算法，考慮了激光掃描偏移等問題，為提高聚類效果，改善原始K-means聚類算法中聚類數量與初始聚類中心需自行設定等問題，結合苗帶特性對其進行了改進，試驗結果表明，基于偏移補償模型的極值點聚類苗帶識別算法軌跡識別精度更高，平均運算耗時更少，算法魯棒性更強，可為田間管理機械進行低傷苗自主導航行駛提供理論依據。但本文算法進行了Y軸縱向等分聚類，在選取不同等分數量時，可能會出現聚類效果不佳現象，后期需進一步改進。文中采集數據的試驗田為機插秧的標準化水田，測量對象無缺失，測量現場環境良好，但實際農田環境中有諸多不確定性因素影響，需要在后期進行進一步的研究。