復(fù)雜環(huán)境中蘋果樹識別與導(dǎo)航線提取方法

張彥斐 魏 鵬 宮金良 蘭玉彬

(1.山東理工大學(xué)農(nóng)業(yè)工程與食品科學(xué)學(xué)院， 淄博 255000； 2.山東理工大學(xué)機械工程學(xué)院， 淄博 255000)

0 引言

農(nóng)業(yè)精準化管理的首要任務(wù)在于對作物的識別提取，其結(jié)果可用于作物種植面積統(tǒng)計、長勢分析與病蟲害監(jiān)測等；農(nóng)業(yè)智慧化管理在于智能農(nóng)機具的自主作業(yè)。對于基于機器視覺導(dǎo)航作業(yè)的農(nóng)業(yè)機器人，作物行中心線的提取極其重要，這也是視覺導(dǎo)航亟待解決的關(guān)鍵問題[1-2]。

盡管蘋果產(chǎn)業(yè)的發(fā)展態(tài)勢較好，但大部分果園的管理方式仍簡單粗放[3]。目前，已有諸多學(xué)者對蘋果果實及果樹的相關(guān)信息提取開展了廣泛研究。實際上，多數(shù)果園環(huán)境存在雜草、遮蔽物、果樹冠層陰影等干擾因素，且果樹之間存在茂盛程度、高低等長勢差異。復(fù)雜環(huán)境下的果樹識別與導(dǎo)航是果園智慧化管理所面臨的核心技術(shù)問題。

目前，農(nóng)業(yè)機械自動導(dǎo)航作業(yè)的實現(xiàn)主要通過機器視覺和衛(wèi)星定位技術(shù)。其中,視覺導(dǎo)航技術(shù)的實現(xiàn)成本低、適應(yīng)性強，已成為農(nóng)業(yè)機械導(dǎo)航的研究熱點[4-7]。現(xiàn)階段田間機器視覺導(dǎo)航多通過斜向視角提取作物田壟中心線或行中心線[8-13]，然后生成導(dǎo)航參數(shù)輔助導(dǎo)航作業(yè)。該方式視野范圍較小、效率低，無法對作業(yè)區(qū)域進行全局路徑規(guī)劃。無人機遙感技術(shù)具有覆蓋范圍廣、快速準確、操作靈活等特點[14-15]；同時，無人機遙感能有效彌補衛(wèi)星遙感成像受云霧等天氣影響的缺點。部分學(xué)者借助無人機遙感影像研究玉米、甘蔗等作物大田的導(dǎo)航線提取[16-19]，但對于果園中蘋果果樹行的導(dǎo)航線提取鮮有研究。

針對具有裸露土壤、遮蔽物、果樹冠層陰影和雜草的復(fù)雜蘋果園環(huán)境，本文在對無人機多光譜遙感影像預(yù)處理的基礎(chǔ)上，對比使用多種機器學(xué)習(xí)算法識別本文最終融合圖像和普通正射影像中的蘋果樹像元并分析各精度；最后以感興趣區(qū)域劃分的方式提取果樹識別結(jié)果中各果樹行導(dǎo)航線，以期為智慧果園中智能農(nóng)業(yè)機器人的自主作業(yè)路徑規(guī)劃提供技術(shù)參考。

1 數(shù)據(jù)采集與處理

1.1 研究區(qū)域概況

本文研究區(qū)位于淄博市沂源縣山東理工大學(xué)生態(tài)無人農(nóng)場智慧果園(36°6′N，118°3′E)，如圖1所示。果園海拔約330 m，屬于溫帶半濕潤大陸性季風(fēng)氣候，地形多為丘陵，年平均氣溫12℃左右，平均降雨量700 mm左右，全年日照時長為2 100 h左右，適合蘋果樹的種植。

圖1 研究區(qū)概況Fig.1 Overview of study area

1.2 數(shù)據(jù)采集

數(shù)據(jù)采集時間為2021年9月30日10:00—12:30，天氣晴朗，風(fēng)力小于3級，采集對象為果園內(nèi)隨機選定的部分種植區(qū)域。研究區(qū)域面積約為7 163.1 m2，蘋果樹高度為4.5～5.5 m。該試驗采用大疆M300型無人機，搭載傳感器為長光禹辰MS600 Pro型多光譜相機。航拍前進行相機灰板定標和研究區(qū)飛行航線規(guī)劃，設(shè)定飛行參數(shù)為高度50 m，速度2.7 m/s，航向重疊度80%，旁向重疊度70%。相機可獲取6個波段數(shù)據(jù)，包括可見光波段、雙紅邊波段和近紅外波段，具體參數(shù)見表1。

表1 MS600 Pro型多光譜相機參數(shù)Tab.1 Parameters of MS600 Pro multispectral camera

1.3 圖像處理方法

1.3.1蘋果樹識別提取流程

使用無人機獲取多光譜遙感影像進行蘋果樹識別提取。首先，將獲取的遙感影像進行拼接、輻射定標等預(yù)處理得到正射影像圖(Digital orthophoto map，DOM)及數(shù)字表面模型(Digital surface model，DSM)圖像。因果樹與雜草等地物存在顯著的高程差異，故DSM圖像可作為區(qū)分兩者的敏感特征。其次，選取計算能夠顯著區(qū)分果樹與雜草的歸一化差異綠度指數(shù)(Normalized difference greenness index，NDGI)[20]、比值植被指數(shù)(Ratio vegetation index，RVI)[21]分布圖，并構(gòu)建NDGI、RVI、DSM融合圖像。此外，由于研究區(qū)域內(nèi)存在明顯裸露的土壤、灌溉設(shè)備等遮蔽物以及果樹冠層陰影，非研究對象像元與蘋果樹像元混雜，會降低提取的精度。故選取過綠植被(Excess green，EXG)指數(shù)[22]和歸一化差異冠層陰影指數(shù)(Normalized difference canopy shadow index，NDCSI)[23]以閾值分割法對融合圖像中的果樹冠層陰影、土壤及遮蔽物進行掩膜剔除處理。最后，通過使用多種機器學(xué)習(xí)算法分別對掩膜后的融合圖像與普通正射影像中的蘋果樹像元進行提取并對比分析。

1.3.2遙感影像預(yù)處理及特征選取

本文使用Yusense Map V2.1.3軟件分析無人機獲取的6波段遙感影像，依次經(jīng)過相機參數(shù)讀取、空中三角測量、波段配準、影像拼接后，通過輻射定標、影像裁剪等操作完成預(yù)處理，最終得到地面分辨率為每像素3.6 cm的DOM影像，見圖2。

圖2 普通正射影像Fig.2 Ordinary orthophoto image

在特征選取方面，考慮到該區(qū)域雜草茂盛且與蘋果樹葉片顏色相近，故選取可顯示綠度差異的NDGI；同時選取RVI解決果樹與雜草的葉片中葉綠素含量不同以及土壤對臨近植被像元的輻射影響問題，既能表征不同植被覆蓋下生物量信息與葉綠素含量的相關(guān)性，又能增強植被與土壤之間的輻射差異。

按文獻[20-23]中對應(yīng)公式獲取各植被指數(shù)分布圖，構(gòu)建蘋果樹識別提取的差異特征圖像集。

1.3.3圖像融合及后處理

將計算得到的NDGI、RVI分布圖與DSM圖像構(gòu)建多特征融合圖像進行蘋果樹的識別提取。具體融合方式為使用ArcGIS校正各圖像后進行像素匹配，采用最鄰近算法進行圖像采樣，然后使用鑲嵌工具分別將體現(xiàn)顏色、葉綠素含量、高程特征差異的NDGI、RVI、DSM圖像進行融合。

其中，為解決各特征數(shù)據(jù)的量綱差異難題以便進行加權(quán)融合，對待融合的各圖像像元值進行歸一化預(yù)處理。將歸一化的NDGI、RVI、DSM圖像中同一位置的像元求均值得到像元均值δ，計算式為

δ=(σNDGI+σRVI+σDSM)/3

(1)

式中σNDGI、σRVI、σDSM——NDGI、RVI、DSM圖像中歸一化處理得到的各點像元值

在3個特征中，雜草與蘋果樹在DSM上的差異最顯著，故在像元均值δ的基礎(chǔ)上添加2倍的DSM圖像歸一化像元值σDSM以突出顯示蘋果樹像元所在的區(qū)域，有

a=δ+2σDSM

(2)

式中a——融合圖像中各點像元值

將各圖像同一位置的像元值按式(1)、(2)計算后作為融合圖像對應(yīng)位置的像元值。獲取的初始融合圖像如圖3所示。

圖3 初始融合圖像Fig.3 Initial fusion image

由于裸露的土壤、遮蔽物等非植被對象的EXG指數(shù)較綠色植被明顯小，為提高蘋果樹的識別提取精度，采用閾值分割法去除非植被對象像元。將EXG指數(shù)分布圖像中各像素值進行數(shù)據(jù)直方圖統(tǒng)計，以等步長累加的方式遍歷直方圖橫坐標，選取明確區(qū)分非植被對象的臨界點數(shù)值作為分割閾值。對NDCSI圖像采用相同方法確定樹冠陰影的最佳分割閾值。

使用ArcGIS重分類工具分別對NDCSI、EXG指數(shù)圖像以最佳分割閾值為斷點進行二值化，然后建立掩膜所需的矢量文件。將矢量文件與融合圖像進行柵格配準后裁除所有非植被對象。

如圖4所示，將掩膜處理后的融合圖像進行濾波去噪、邊界清理等處理后作為果樹識別的最終圖像，用于探究不同分類方法對普通正射影像和最終融合圖像中的蘋果樹識別提取效果。

圖4 最終融合圖像Fig.4 Final fusion image

如圖5所示，研究區(qū)內(nèi)原為樹冠陰影以及裸露土壤、灌溉設(shè)備等遮蔽物的像元已被較好剔除，大大降低了非植被混雜像元對果樹提取的干擾，有效提高識別精度。

圖5 非植被像元剔除效果對比Fig.5 Comparison of removal effects of non-vegetation pixels

2 蘋果樹提取

2.1 蘋果樹樣本集合構(gòu)建與識別提取方法

針對研究區(qū)多光譜遙感正射影像和最終融合圖像，分別隨機選取60個樣本作為蘋果樹識別訓(xùn)練集，40個樣本作為驗證集。其中，樣本采樣區(qū)域尺寸均為40像素×30像素，且訓(xùn)練樣本與驗證樣本不重疊。分別采用支持向量機(Support vector machine，SVM)法、最大似然(Maximum like-lihood classification，MLC)法和隨機森林(Random forest，RF)法3種方法[24-28]提取2種圖像中的蘋果樹像元，對比分析其實際效果。

2.2 蘋果樹識別提取精度評價與效果分析

對各方法的識別提取結(jié)果，采用計算混淆矩陣的方式給出總體精度和Kappa系數(shù)對識別精度進行評價，以制圖精度和用戶精度對識別準確性進行評判[29-30]。

如圖6、7所示，MLC法對最終融合圖像中果樹像元的提取效果最優(yōu)，其用戶精度、制圖精度、總體分類精度、Kappa系數(shù)分別為88.57%、93.93%、93.00%、0.882 4。

圖6 3種方法對2種圖像中蘋果樹的分類精度比較Fig.6 Statistical comparison of classification accuracy of apple trees in two images by three methods

圖7 3種方法對2種圖像中蘋果樹的Kappa系數(shù)比較Fig.7 Statistical comparison of Kappa coefficients of apple trees in two images by three methods

此外，相對于普通正射影像，本文構(gòu)建的最終融合圖像使3種方法的識別精度均得到有效提升。其中，總體分類精度提升均超10個百分點，最低為13個百分點；對RF法的用戶精度提升幅度最大，為27.12個百分點；對SVM法的制圖精度提升幅度最大，為9.03個百分點；對SVM法的Kappa系數(shù)提升幅度最大，為22.55%，且對其余2種方法的提升均在20%以上。

圖8為對普通正射影像和最終融合圖像中蘋果樹像元提取中Kappa系數(shù)最高的MLC法提取結(jié)果。對比圖8a、8b可見，針對具有較高提取精度的相同算法，本文構(gòu)建最終融合圖像的方式可更好地從包含裸露土壤、遮蔽物、雜草、陰影等非果樹對象的背景中有效提取蘋果樹像元，實現(xiàn)復(fù)雜環(huán)境下高精度提取蘋果樹像元。

圖8 MLC法對圖像中蘋果樹的提取結(jié)果Fig.8 Extraction results of two different images by MLC method

3 蘋果樹行中心線提取

智能機器人在果園自主作業(yè)時，可根據(jù)果樹行中心線進行輔助定位和導(dǎo)航作業(yè)。目前作物行中心線的提取方法主要采用標準或改進Hough變換的方法。但前者易受農(nóng)田環(huán)境影響，存在計算成本高、峰值檢測難等弊端；后者雖降低了計算功耗，但在復(fù)雜環(huán)境下的準確性和實時性有待提升[31-32]。基于最小二乘法擬合作物行特征點獲取中心線方法具有實時性好、準確性高等優(yōu)點[33-34]。

3.1 提取流程與行特征點識別

蘋果樹行中心線提取流程為：在效果最優(yōu)的MLC法對最終融合圖像中蘋果樹像元提取結(jié)果的基礎(chǔ)上，降噪處理后進行二值化、形態(tài)學(xué)轉(zhuǎn)換，再以行為單位劃分感興趣區(qū)域提取果樹輪廓，計算各輪廓質(zhì)心作為行特征點，最終以最小二乘法擬合各行特征點得到果樹行中心線。

由于MLC法對融合圖像中蘋果樹像元提取的精度限制，提取結(jié)果在果樹像元附近往往存在雜草及離散噪聲像元。為保證導(dǎo)航線提取精度，讀取并灰度化結(jié)果圖像，采用高斯核參數(shù)為(5,5)對圖像進行高斯降噪。設(shè)定初始閾值為0并使用OTSU算法重新對圖像進行二值化分割。為避免臨近果樹的雜草像元被錯誤提取的影響，設(shè)定閉運算的卷積核參數(shù)為(9,9)，對圖像先腐蝕再膨脹，然后設(shè)定開運算的卷積核參數(shù)為(7,7)對圖像進行先膨脹再腐蝕的形態(tài)學(xué)轉(zhuǎn)換處理，以便去除輪廓前景中的孔洞和黑點，保證提取輪廓的完整性。

由于研究區(qū)域內(nèi)各行果樹的分布及長勢存在差異，本文選擇平行四邊形掩膜窗口作為感興趣區(qū)域，對各行果樹逐一進行輪廓提取。求取各果樹輪廓的特征矩mpq后，再計算各輪廓質(zhì)心橫縱坐標(Cx，Cy)，并將其作為果樹行中心線提取的特征點。計算式為

(3)

(4)

式中x、y——輪廓中各點橫、縱坐標

M、N——輪廓中各點最大橫、縱坐標

f(x,y)——輪廓中各點像素值

p、q——階數(shù)

如圖9所示，各果樹行長勢的差異經(jīng)形態(tài)學(xué)變換處理后轉(zhuǎn)換為輪廓數(shù)量的差異，長勢稀疏的果樹行得到的輪廓數(shù)量更多，即獲取的特征點數(shù)量更多，故本文方法更利于對長勢分布差異明顯且長勢稀疏的果樹行中心線的提取。

圖9 果樹行特征點提取結(jié)果Fig.9 Fruit tree row feature point extraction results

3.2 中心線提取效果與分析

采用人工擬合的果樹行中心線作為比較基準[13]，評判本文方法的準確性。依據(jù)實際作業(yè)場景需求，設(shè)定本文方法提取的中心線與人工擬合中心線之間的夾角誤差小于5°時為有效提取，且依據(jù)兩者之間夾角評價提取的準確性。同時，將本文方法與改進的Hough變換算法在相同圖像預(yù)處理基礎(chǔ)上的提取結(jié)果進行對比分析。所用計算機CPU為Intel(R) Xeon(R) E5-1650 v4，主頻為3.60 GHz，內(nèi)存64 GB，操作系統(tǒng)為Windows 11專業(yè)工作站版64位系統(tǒng)，IDE平臺為Spyder 4，編程語言為Python。

圖10 果樹中心線提取結(jié)果Fig.10 Extraction results of fruit tree row center lines by different methods

圖10a為改進Hough變換算法對圖像各感興趣區(qū)域內(nèi)經(jīng)降噪和形態(tài)學(xué)處理的果樹行中心線的提取結(jié)果。可見，即便經(jīng)過一系列圖像預(yù)處理和感興趣區(qū)域劃分，該算法仍無法滿足需求，且10次測試的平均用時為4.988 9 s。

由圖10b可知，使用最小二乘法擬合特征點得到的中心線角度偏差均小于5°，平均角度偏差為0.597 5°，10次測試平均用時為0.402 3 s。

圖11 各行導(dǎo)航線提取偏差對比Fig.11 Comparison of extraction deviation of each row of navigation line

圖11為各果樹行兩線角度偏差的可視化比較結(jié)果。其中，橫坐標為圖10b中從左至右的各果樹行編號。結(jié)合圖10b、11進行整體分析，相較于人工擬合導(dǎo)航線結(jié)果，本文方法對長勢較差的第6～9果樹行提取的導(dǎo)航線角度偏差較大，但均小于2.25°；其余各行導(dǎo)航線偏差大部分小于平均偏差，小部分位于其附近。此外，對各行提取偏差進行統(tǒng)計計算，樣本標準差僅為0.135 0°，說明本文方法具有較好的穩(wěn)定性。

4 結(jié)論

(1)提出一種能夠從包含裸露土壤、遮蔽物、果樹冠層陰影和雜草的復(fù)雜環(huán)境中精準提取蘋果樹像元的方法。

(2)根據(jù)無人機多光譜影像得到DOM、DSM圖像，計算NDGI、RVI分布圖，將歸一化的NDGI、RVI、DSM圖像中同一位置的像元求均值后，添加2倍的σDSM構(gòu)建融合圖像。綜合利用NDCSI、EXG指數(shù)剔除融合圖像中的非植被像元構(gòu)建最終融合圖像。采用SVM、MLC、RF方法分別提取最終融合圖像和普通正射影像中的蘋果樹像元，并計算混淆矩陣評價各方法的識別精度。

(3)試驗結(jié)果表明，MLC法對最終融合圖像中果樹像元的提取效果最優(yōu)，其用戶精度、制圖精度、總體分類精度、Kappa系數(shù)分別為88.57%、93.93%、93.00%、0.882 4。較普通正射影像，本文構(gòu)建的融合圖像對3種方法的識別精度均有提升。其對RF法的用戶精度提升幅度最大，為27.12個百分點；對SVM法的制圖精度提升幅度最大，為9.03個百分點；對3種方法的總體分類精度提升幅度最低為13個百分點；對SVM法的Kappa系數(shù)提升幅度最大，為22.55%，且對其余2種方法的提升也均在20%以上。

(4)以感興趣區(qū)域劃分的方式提取果樹識別結(jié)果中的各果樹行導(dǎo)航線，其平均角度偏差為0.597 5°，10次測試平均用時僅0.402 3 s且具有較好的穩(wěn)定性，可大大提高農(nóng)業(yè)智能裝備的自動導(dǎo)航精度和實時性。