基于無人機遙感技術的大豆苗數估算研究*

李金陽，張偉，康燁，許秀英，亓立強，石文強

(1. 黑龍江八一農墾大學工程學院，黑龍江大慶，163319；2. 黑龍江省保護性耕作工程技術研究中心，黑龍江大慶，163319)

0 引言

大豆苗數估算研究對于大豆田間出苗情況的評估起著重要作用。無人機遙感技術可以實現大范圍內空間圖像等相關信息的快速獲取[1]，目前，已有學者應用無人機遙感技術開展了對不同作物苗數估算的研究，但專門針對苗期大豆苗數估算的研究較少。

Feng等[2]利用無人機獲取棉花幼苗圖像，基于深度學習算法用于棉花出苗評估，誤差為4.4%。劉志等[3]通過無人機獲取玉米幼苗圖像，對玉米幼苗的識別和分割結果有助于實際生產過程中對玉米植株的統計和出苗率的計算。Sankaran等[4]通過在無人機圖像上選擇適當的多邊形閾值區域，實現了對馬鈴薯的出苗估測，但與實際計數誤差較大。Jin等[5]基于無人機圖像，使用支持向量機及粒子群優化算法估算了苗期小麥株數和密度，方法的可重復性很好，但數據處理過程繁雜。Liu等[6]通過小麥苗期無人機圖像，使用區域定位算法計算缺苗區域長度，為有效測量小麥出苗均勻性提供了一種新方法。Gn?dinger等[7]通過對不同生長階段的玉米植株無人機圖像進行增強顏色對比度和創建閾值處理后估算株數，誤差≤5%，但雜草會影響其結果的準確性。

為了提高大豆苗數估算的準確性、解決因田間雜草影響較大帶來的出苗估測誤差較大等問題，以苗期獲取的無人機圖像為研究對象，通過分析不同植被指數、直方圖均衡化及Otsu閾值算法對大豆目標的提取結果，開展大豆苗數估算研究，以期快速評估大豆田間出苗情況，更好的服務生產實際問題。

1 材料和方法

1.1 試驗地概況

試驗在廣東省湛江市南亞熱帶作物研究所循環農業研究中心進行，試驗地塊地理位置為東經110°20′、北緯21°10′之間，土壤類型為紅壤土。

選取試驗地大小為54 m×11 m，采用1.1 m壟上三行的種植模式，即在1.1 m的壟臺上均勻種植三行大豆，充分利用壟臺的邊際效應，通過減小大豆植株間的行距來擴大種植密度進而實現密植栽培，提高單位面積的保苗株數來提高作物產量[8]。本試驗中最外側兩行大豆行距為450 mm，每兩行行距為225 mm。

1.2 圖像獲取及預處理

無人機采集平臺提供了靈活和即時的圖像處理功能，其可變的飛行高度和更好的圖像分辨率等處理選項適用于面向大豆目標的使用[9-10]。試驗中采用大疆精靈4多光譜版無人機作為圖像數據采集平臺，通過DJI GS Pro軟件設置無人機飛行參數，包括在指定區域內設置航線并完成相關數據的接收和處理過程。為了在較短的時間內獲得所需圖像，無人機采用“S”型的飛行路線。試驗航向重復率設置為80%，旁向重復率設置為70%，飛行高度設置為5 m。相機型號為P4 Multispectral Camera，相機朝向沿航線方向，拍攝時鏡頭垂直向下，采用航點懸停拍照模式。

利用無人機飛行平臺采集的圖像具有很大的重疊率，如果直接用來提取苗情信息，難以反映大豆實際生長情況，因此，對圖像進行預處理操作，即圖像拼接、圖像裁剪和圖像切割。無人機飛行過程共拍攝圖像3 650幅，獲取無人機圖像的POS數據，利用Agisoft PhotoScan軟件對圖像進行拼接處理，如圖1所示為圖像拼接流程圖。

圖1 圖像拼接流程圖Fig. 1 Image mosaic flow chart

由于拼接時圖像邊緣數據畸變較大，存在較多異常值，這會影響算法設計和模型的建立。因此，需要在拼接的圖像上手動裁剪，剔除邊緣區域，選擇合適的研究區。為便于后續試驗開展，將研究區切割成24份大小一致的區域，如圖2所示。

圖2 圖像切割示意圖Fig. 2 Image cutting schematic

在圖像切割的過程中，產生的椒鹽噪聲會對圖像的分割與識別以及特征提取等操作帶來影響。為了保留圖像的邊緣輪廓、使圖像清晰，使用中值濾波算法對采集的圖像進行處理，去除椒鹽噪聲。如圖3所示為切割后原始圖像，如圖4所示為中值濾波效果圖。

圖4 中值濾波效果圖Fig. 4 Median filtering renderings

1.3 大豆目標的識別與分割

本文采用的算法流程圖如圖5所示。

圖5 算法流程圖Fig. 5 Algorithm flowchart

圖像分割作為圖像處理中的關鍵一環，對圖像分析結果有很大的影響，圖像分割通過把圖像劃分為若干個區域，從中提取出感興趣區域,是精準農業、植物表型研究的關鍵技術[11]，可以高效、無損的獲取農作物信息，實時了解作物生長動態，更好的管理作物。通過采用多種方法對大豆目標進行提取，選擇出可用于大豆目標分割的最佳方法。

1.3.1 基于植被指數的大豆目標識別與分割

目標識別和分割的關鍵是將識別對象與背景分離，植被指數在農業作物鑒定中應用廣泛[12]。大豆田間可見光圖像由豆苗、雜草、滴灌管和土壤四部分構成，豆苗和雜草呈綠色、土壤呈紅棕色、滴灌管呈黑色，根據顏色特征可采用植被指數對圖像進行識別與分割，提取圖像中的大豆目標。表1列舉了各候選植被指數公式。

表1 植被指數列表Tab. 1 Vegetation indices list

1.3.2 基于直方圖均衡化的大豆目標識別與分割

為避免分別對彩色圖像的R、G、B三個分量做直方圖均衡化導致結果圖像色彩失真的現象發生，將RGB原色空間轉換為HSV顏色空間后分別對H、S、V三個分量進行直方圖均衡化操作，如圖6所示為RGB原色空間與HSV顏色空間轉換示意圖。

圖6 RGB原色空間與HSV顏色空間轉換示意圖Fig. 6 Schematic diagram of the conversion between RGB primary color space and HSV color space

HSV顏色空間是一種包含色調(Hue)、飽和度(Saturation)、亮度(Value)的視覺顏色模型[13]。其通常用一個圓錐體來表示，圓錐體內部每個點代表一種特定的顏色，角度位置表示顏色的色調，到中軸線的距離表示飽和度，高度表示亮度。RGB轉HSV公式如式(1)～式(3)所示。

(1)

(2)

v=max

(3)

1.3.3 基于Otsu閾值算法的大豆目標識別與分割

選取Otsu閾值算法來自動確定閾值，進行大豆目標的分割。Otsu閾值算法即最大類間方差法，是一種自適應的閾值提取方法。這種方法的核心思想是通過一個最佳閾值，使得目標和背景的類間方差最大，不需要人為設定其他參數，是一種自動選擇閾值的方法，計算過程簡單、穩定[14]。其快速運算、錯分概率小的特點，適合于大豆目標的識別和提取，如式(4)、式(5)所示。

gm=w0·g0+w1·g1

(4)

S=w0(g0-gm)2+w1(g1-gm)2

(5)

式中：w0——目標的像元比例；

w1——背景的像元比例；

g0——目標的平均灰度；

g1——背景的平均灰度；

gm——影像總的平均灰度；

S——類間方差。

1.3.4 雜草去除及苗數估算

田地中存在少量的雜草，基于無人機可見光圖像進行苗期大豆目標提取時，豆苗和雜草因其相似的顏色特征，很難直接分開，這使得分割后的圖像中會存在一些雜草噪聲，因此，需要進一步對分割后的圖像進行處理，才可用于后續的大豆苗情指標提取。

根據圖像的空間分辨率，同時結合圖像獲取時大豆的生長情況及對真彩色圖像的解譯結果，選擇0.16 cm2為閾值，剔除面積小于此值的大豆目標。該操作主要是對圖像中的異常值進行剔除，即對面積過小的大豆目標、雜草進行剔除。

結合圖像中豆苗目標大、雜草目標小的特點，嘗試采用腐蝕膨脹及開、閉運算等形態學方法繼續去除大豆目標中的雜草噪聲。腐蝕膨脹即通過一定的結構元與圖像進行相應的邏輯運算，如果圖像中完全包含腐蝕算子則保留該區域，否則丟棄該區域，這使得圖像中目標會向內收縮，丟棄掉小的細節和噪聲；而膨脹與此相反，會使目標邊界向外擴張，起到放大細節的作用[15]。通過選擇合適的結構元進行腐蝕操作去除雜草噪聲，再通過膨脹操作放大大豆目標的細節特征，此類操作并不會明顯改變物體的原先輪廓，其具體公式如式(6)、式(7)[16]所示。

(f?B)(x,y)=min{f(x+i,y+j)-B(i,j)|(x+i,y+j)∈Df(i,j)∈DB}

(6)

(f⊕B)(x,y)=max{f(x-i,y-j)+B(i,j)|(x-i,y-j)∈Df(i,j)∈DB}

(7)

式中：f?B——腐蝕；

f⊕B——膨脹；

f(x,y)——輸入圖像；

B(i,j)——結構元素；

Df——輸入元素的定義域；

DB——結構元素的定義域。

通過連通區域分析對大豆苗數進行估算，連通區域分析在圖像處理的過程中較為常用，通過對二值圖像中的白色像素進行標記，讓每個單獨的連通區域形成一個被標識的塊，進一步的就可以獲取這些塊兒的輪廓、外接矩形、質心、不變矩等幾何參數，自動計算大豆苗數。

2 結果與分析

2.1 大豆目標識別與分割結果分析

顏色直方圖有助于判斷、校正和優化圖像的亮度和對比度，有助于評估圖像的質量，并為進一步處理收集更多信息。對H、S、V分量進行提取并做均衡化處理后的直方圖如圖7所示。

(a) H分量提取效果

(b) S分量提取效果

(c) V分量提取效果

(d) H均衡化后直方圖

(e) S均衡化后直方圖

(f) V均衡化后直方圖圖7 H，S和V分量提取效果及均衡化后直方圖結果Fig. 7 Extraction effect of H,S and V components and histogram results after equalization

如圖7(d)、圖7(e)、圖7(f)所示，對直方圖均衡化后的圖像做了歸一化處理，當灰度值=0時，表示黑色；當灰度值=1時，表示白色；當灰度值∈(0,1)時，表示像素灰度在黑白之間變化。

對24個區域圖像分別使用不同的植被指數、直方圖均衡化、Otsu閾值算法及本文所提出的方法進行分析，選取其中1個區域對大豆目標識別與分割效果如圖8所示。

盡管在RGB圖像中各部分具備顏色差異，但通過分析GNDVI、LCI、NDRE、NDVI、OSAVI五種指數對圖像的分割效果：土壤噪聲被完全抑制，大豆目標與雜草因顏色相近，難以區分，滴灌管分割效果顯著，但未能完全去除，不利于后續采用形態學算法提取大豆目標；進行直方圖均衡化的分割效果顯示：大豆目標與土壤、滴灌管之間存在一定的間隔區域，圖像直方圖變成近似均勻分布，圖像對比度顯著增強，顏色特征明顯，但圖像中噪聲較多；如圖8(g)所示，Otsu閾值算法對大豆目標提取效果明顯，冗余噪聲較少；基于本文所提出的Otsu閾值算法結合過綠指數的方法對大豆目標的識別與分割效果如圖8(h)所示，其提取綠色植物圖像效果較好，陰影、滴灌管和土壤噪聲等均能較明顯的被抑制，植物圖像更為突出，分割更為完整。

結果表明，采用Otsu閾值算法結合過綠指數用于大豆目標的識別與分割方法是可行的，且根據此種方法分離綠色作物的可行性已有研究[17]。

(a) GNDVI分割

(b) LCI分割

(c) NDRE分割

(d) NDVI分割

(e) OSAVI分割

(f) 直方圖均衡化

2.2 雜草去除及苗數估算結果分析

根據圖像的空間分辨率選擇0.16 cm2為閾值，剔除面積小于此值的大豆目標，圖9所示為部分區域異常值剔除效果圖。

(a) 去雜草前

(b) 去雜草后圖9 去雜草效果圖Fig. 9 Remove weeds renderings

使用Matlab中連通區域標記函數bwlabel，通過一次遍歷膨脹圖像，并記下每一行(或列)中連續的團(run)和標記的等價對，通過等價對對原來的圖像進行重新標記，得到大豆出苗情況。將基于Otsu閾值算法結合過綠指數分割圖像轉換為二值化圖像，通過形態學腐蝕、膨脹操作去除雜草噪聲效果及連通區域分析效果如圖10所示。

(a) 腐蝕操作

(b) 膨脹操作

(c) 連通區域計算出苗數量圖10 形態學操作及連通區域計算結果Fig. 10 Morphological operation and calculation results of connected regions

大豆播種時間為2020年12月1日，播種后5 d左右開始出苗，試驗獲取的是田間自然生長狀態下的圖像數據，圖像采集時間為2020年12月31日，每2 d獲取一次。研究區大豆為一穴一粒的播種方式，每穴有一株幼苗長出即為出苗，在無人機獲取大豆田間圖像時，人工對區域內大豆進行統計，以獲得大豆苗數，根據式(8)出苗統計誤差公式計算實測苗數和估算苗數誤差，實測苗數和估算苗數試驗結果如表2所示。

(8)

式中：W——誤差；

R——實際測得的苗數；

A——估算的苗數。

通過計算得出實測苗數與估算苗數的平均誤差為0.43%，根據表2所示試驗結果，通過Origin 2019b軟件構建如圖11所示的線性回歸模型。由線性回歸模型得出R2為0.909 4，線性回歸模型擬合度相對較高，能準確地通過無人機遙感技術結合苗數估算方法進行大豆苗數估算。

表2 試驗結果Tab. 2 Results of experiments

圖11 實測苗數與估算苗數的比較Fig. 11 Comparison of measured and estimated seedling numbers

3 結論

以苗期獲取的大豆無人機圖像為研究對象，最終選擇Otsu閾值算法結合過綠指數分割大豆圖像，根據空間分辨率剔除圖像中異常值，通過形態學腐蝕、膨脹操作去除雜草噪聲后采用連通區域分析估算大豆苗數。

1) 使用Otsu閾值算法結合過綠指數分割大豆圖像，噪聲能較明顯的被抑制，豆苗圖像更為突出，分割更為完整，形態學操作降低了雜草帶來的影響，提高了識別精度。

2) 基于無人機遙感技術的估算苗數與田間實測苗數建立的線性回歸模型中，相關系數R2為0.909 4，由出苗統計誤差公式計算的平均誤差為0.43%，提出的苗數估算方法誤差較低，能夠快速準確識別大豆苗數。

3) 研究可快速評估大豆最終田間出苗，為田間大豆出苗率的測定提供了有力的支持，為后續對播種效果和質量評估工作奠定了基礎，為農業生產者進行田間管理提供參考決策意見。