李金陽(yáng)，張偉，康燁，許秀英，亓立強(qiáng)，石文強(qiáng)

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

0 引言

大豆苗數(shù)估算研究對(duì)于大豆田間出苗情況的評(píng)估起著重要作用。無(wú)人機(jī)遙感技術(shù)可以實(shí)現(xiàn)大范圍內(nèi)空間圖像等相關(guān)信息的快速獲取[1]，目前，已有學(xué)者應(yīng)用無(wú)人機(jī)遙感技術(shù)開(kāi)展了對(duì)不同作物苗數(shù)估算的研究，但專(zhuān)門(mén)針對(duì)苗期大豆苗數(shù)估算的研究較少。

Feng等[2]利用無(wú)人機(jī)獲取棉花幼苗圖像，基于深度學(xué)習(xí)算法用于棉花出苗評(píng)估，誤差為4.4%。劉志等[3]通過(guò)無(wú)人機(jī)獲取玉米幼苗圖像，對(duì)玉米幼苗的識(shí)別和分割結(jié)果有助于實(shí)際生產(chǎn)過(guò)程中對(duì)玉米植株的統(tǒng)計(jì)和出苗率的計(jì)算。Sankaran等[4]通過(guò)在無(wú)人機(jī)圖像上選擇適當(dāng)?shù)亩噙呅伍撝祬^(qū)域，實(shí)現(xiàn)了對(duì)馬鈴薯的出苗估測(cè)，但與實(shí)際計(jì)數(shù)誤差較大。Jin等[5]基于無(wú)人機(jī)圖像，使用支持向量機(jī)及粒子群優(yōu)化算法估算了苗期小麥株數(shù)和密度，方法的可重復(fù)性很好，但數(shù)據(jù)處理過(guò)程繁雜。Liu等[6]通過(guò)小麥苗期無(wú)人機(jī)圖像，使用區(qū)域定位算法計(jì)算缺苗區(qū)域長(zhǎng)度，為有效測(cè)量小麥出苗均勻性提供了一種新方法。Gn?dinger等[7]通過(guò)對(duì)不同生長(zhǎng)階段的玉米植株無(wú)人機(jī)圖像進(jìn)行增強(qiáng)顏色對(duì)比度和創(chuàng)建閾值處理后估算株數(shù)，誤差≤5%，但雜草會(huì)影響其結(jié)果的準(zhǔn)確性。

為了提高大豆苗數(shù)估算的準(zhǔn)確性、解決因田間雜草影響較大帶來(lái)的出苗估測(cè)誤差較大等問(wèn)題，以苗期獲取的無(wú)人機(jī)圖像為研究對(duì)象，通過(guò)分析不同植被指數(shù)、直方圖均衡化及Otsu閾值算法對(duì)大豆目標(biāo)的提取結(jié)果，開(kāi)展大豆苗數(shù)估算研究，以期快速評(píng)估大豆田間出苗情況，更好的服務(wù)生產(chǎn)實(shí)際問(wèn)題。

1 材料和方法

1.1 試驗(yàn)地概況

試驗(yàn)在廣東省湛江市南亞熱帶作物研究所循環(huán)農(nóng)業(yè)研究中心進(jìn)行，試驗(yàn)地塊地理位置為東經(jīng)110°20′、北緯21°10′之間，土壤類(lèi)型為紅壤土。

選取試驗(yàn)地大小為54 m×11 m，采用1.1 m壟上三行的種植模式，即在1.1 m的壟臺(tái)上均勻種植三行大豆，充分利用壟臺(tái)的邊際效應(yīng)，通過(guò)減小大豆植株間的行距來(lái)擴(kuò)大種植密度進(jìn)而實(shí)現(xiàn)密植栽培，提高單位面積的保苗株數(shù)來(lái)提高作物產(chǎn)量[8]。本試驗(yàn)中最外側(cè)兩行大豆行距為450 mm，每?jī)尚行芯酁?25 mm。

1.2 圖像獲取及預(yù)處理

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

利用無(wú)人機(jī)飛行平臺(tái)采集的圖像具有很大的重疊率，如果直接用來(lái)提取苗情信息，難以反映大豆實(shí)際生長(zhǎng)情況，因此，對(duì)圖像進(jìn)行預(yù)處理操作，即圖像拼接、圖像裁剪和圖像切割。無(wú)人機(jī)飛行過(guò)程共拍攝圖像3 650幅，獲取無(wú)人機(jī)圖像的POS數(shù)據(jù)，利用Agisoft PhotoScan軟件對(duì)圖像進(jìn)行拼接處理，如圖1所示為圖像拼接流程圖。

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

由于拼接時(shí)圖像邊緣數(shù)據(jù)畸變較大，存在較多異常值，這會(huì)影響算法設(shè)計(jì)和模型的建立。因此，需要在拼接的圖像上手動(dòng)裁剪，剔除邊緣區(qū)域，選擇合適的研究區(qū)。為便于后續(xù)試驗(yàn)開(kāi)展，將研究區(qū)切割成24份大小一致的區(qū)域，如圖2所示。

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

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

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

1.3 大豆目標(biāo)的識(shí)別與分割

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

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

圖像分割作為圖像處理中的關(guān)鍵一環(huán)，對(duì)圖像分析結(jié)果有很大的影響，圖像分割通過(guò)把圖像劃分為若干個(gè)區(qū)域，從中提取出感興趣區(qū)域,是精準(zhǔn)農(nóng)業(yè)、植物表型研究的關(guān)鍵技術(shù)[11]，可以高效、無(wú)損的獲取農(nóng)作物信息，實(shí)時(shí)了解作物生長(zhǎng)動(dòng)態(tài)，更好的管理作物。通過(guò)采用多種方法對(duì)大豆目標(biāo)進(jìn)行提取，選擇出可用于大豆目標(biāo)分割的最佳方法。

1.3.1 基于植被指數(shù)的大豆目標(biāo)識(shí)別與分割

目標(biāo)識(shí)別和分割的關(guān)鍵是將識(shí)別對(duì)象與背景分離，植被指數(shù)在農(nóng)業(yè)作物鑒定中應(yīng)用廣泛[12]。大豆田間可見(jiàn)光圖像由豆苗、雜草、滴灌管和土壤四部分構(gòu)成，豆苗和雜草呈綠色、土壤呈紅棕色、滴灌管呈黑色，根據(jù)顏色特征可采用植被指數(shù)對(duì)圖像進(jìn)行識(shí)別與分割，提取圖像中的大豆目標(biāo)。表1列舉了各候選植被指數(shù)公式。

表1 植被指數(shù)列表Tab. 1 Vegetation indices list

1.3.2 基于直方圖均衡化的大豆目標(biāo)識(shí)別與分割

為避免分別對(duì)彩色圖像的R、G、B三個(gè)分量做直方圖均衡化導(dǎo)致結(jié)果圖像色彩失真的現(xiàn)象發(fā)生，將RGB原色空間轉(zhuǎn)換為HSV顏色空間后分別對(duì)H、S、V三個(gè)分量進(jìn)行直方圖均衡化操作，如圖6所示為RGB原色空間與HSV顏色空間轉(zhuǎn)換示意圖。

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

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

(1)

(2)

v=max

(3)

1.3.3 基于Otsu閾值算法的大豆目標(biāo)識(shí)別與分割

選取Otsu閾值算法來(lái)自動(dòng)確定閾值，進(jìn)行大豆目標(biāo)的分割。Otsu閾值算法即最大類(lèi)間方差法，是一種自適應(yīng)的閾值提取方法。這種方法的核心思想是通過(guò)一個(gè)最佳閾值，使得目標(biāo)和背景的類(lèi)間方差最大，不需要人為設(shè)定其他參數(shù)，是一種自動(dòng)選擇閾值的方法，計(jì)算過(guò)程簡(jiǎn)單、穩(wěn)定[14]。其快速運(yùn)算、錯(cuò)分概率小的特點(diǎn)，適合于大豆目標(biāo)的識(shí)別和提取，如式(4)、式(5)所示。

gm=w0·g0+w1·g1

(4)

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

(5)

式中：w0——目標(biāo)的像元比例；

w1——背景的像元比例；

g0——目標(biāo)的平均灰度；

g1——背景的平均灰度；

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

S——類(lèi)間方差。

1.3.4 雜草去除及苗數(shù)估算

田地中存在少量的雜草，基于無(wú)人機(jī)可見(jiàn)光圖像進(jìn)行苗期大豆目標(biāo)提取時(shí)，豆苗和雜草因其相似的顏色特征，很難直接分開(kāi)，這使得分割后的圖像中會(huì)存在一些雜草噪聲，因此，需要進(jìn)一步對(duì)分割后的圖像進(jìn)行處理，才可用于后續(xù)的大豆苗情指標(biāo)提取。

根據(jù)圖像的空間分辨率，同時(shí)結(jié)合圖像獲取時(shí)大豆的生長(zhǎng)情況及對(duì)真彩色圖像的解譯結(jié)果，選擇0.16 cm2為閾值，剔除面積小于此值的大豆目標(biāo)。該操作主要是對(duì)圖像中的異常值進(jìn)行剔除，即對(duì)面積過(guò)小的大豆目標(biāo)、雜草進(jìn)行剔除。

結(jié)合圖像中豆苗目標(biāo)大、雜草目標(biāo)小的特點(diǎn)，嘗試采用腐蝕膨脹及開(kāi)、閉運(yùn)算等形態(tài)學(xué)方法繼續(xù)去除大豆目標(biāo)中的雜草噪聲。腐蝕膨脹即通過(guò)一定的結(jié)構(gòu)元與圖像進(jìn)行相應(yīng)的邏輯運(yùn)算，如果圖像中完全包含腐蝕算子則保留該區(qū)域，否則丟棄該區(qū)域，這使得圖像中目標(biāo)會(huì)向內(nèi)收縮，丟棄掉小的細(xì)節(jié)和噪聲；而膨脹與此相反，會(huì)使目標(biāo)邊界向外擴(kuò)張，起到放大細(xì)節(jié)的作用[15]。通過(guò)選擇合適的結(jié)構(gòu)元進(jìn)行腐蝕操作去除雜草噪聲，再通過(guò)膨脹操作放大大豆目標(biāo)的細(xì)節(jié)特征，此類(lèi)操作并不會(huì)明顯改變物體的原先輪廓，其具體公式如式(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)——結(jié)構(gòu)元素；

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

DB——結(jié)構(gòu)元素的定義域。

通過(guò)連通區(qū)域分析對(duì)大豆苗數(shù)進(jìn)行估算，連通區(qū)域分析在圖像處理的過(guò)程中較為常用，通過(guò)對(duì)二值圖像中的白色像素進(jìn)行標(biāo)記，讓每個(gè)單獨(dú)的連通區(qū)域形成一個(gè)被標(biāo)識(shí)的塊，進(jìn)一步的就可以獲取這些塊兒的輪廓、外接矩形、質(zhì)心、不變矩等幾何參數(shù)，自動(dòng)計(jì)算大豆苗數(shù)。

2 結(jié)果與分析

2.1 大豆目標(biāo)識(shí)別與分割結(jié)果分析

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

(a) H分量提取效果

(b) S分量提取效果

(c) V分量提取效果

(d) H均衡化后直方圖

(e) S均衡化后直方圖

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

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

對(duì)24個(gè)區(qū)域圖像分別使用不同的植被指數(shù)、直方圖均衡化、Otsu閾值算法及本文所提出的方法進(jìn)行分析，選取其中1個(gè)區(qū)域?qū)Υ蠖鼓繕?biāo)識(shí)別與分割效果如圖8所示。

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

結(jié)果表明，采用Otsu閾值算法結(jié)合過(guò)綠指數(shù)用于大豆目標(biāo)的識(shí)別與分割方法是可行的，且根據(jù)此種方法分離綠色作物的可行性已有研究[17]。

(a) GNDVI分割

(b) LCI分割

(c) NDRE分割

(d) NDVI分割

(e) OSAVI分割

(f) 直方圖均衡化

2.2 雜草去除及苗數(shù)估算結(jié)果分析

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

(a) 去雜草前

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

使用Matlab中連通區(qū)域標(biāo)記函數(shù)bwlabel，通過(guò)一次遍歷膨脹圖像，并記下每一行(或列)中連續(xù)的團(tuán)(run)和標(biāo)記的等價(jià)對(duì)，通過(guò)等價(jià)對(duì)對(duì)原來(lái)的圖像進(jìn)行重新標(biāo)記，得到大豆出苗情況。將基于Otsu閾值算法結(jié)合過(guò)綠指數(shù)分割圖像轉(zhuǎn)換為二值化圖像，通過(guò)形態(tài)學(xué)腐蝕、膨脹操作去除雜草噪聲效果及連通區(qū)域分析效果如圖10所示。

(a) 腐蝕操作

(b) 膨脹操作

(c) 連通區(qū)域計(jì)算出苗數(shù)量圖10 形態(tài)學(xué)操作及連通區(qū)域計(jì)算結(jié)果Fig. 10 Morphological operation and calculation results of connected regions

大豆播種時(shí)間為2020年12月1日，播種后5 d左右開(kāi)始出苗，試驗(yàn)獲取的是田間自然生長(zhǎng)狀態(tài)下的圖像數(shù)據(jù)，圖像采集時(shí)間為2020年12月31日，每2 d獲取一次。研究區(qū)大豆為一穴一粒的播種方式，每穴有一株幼苗長(zhǎng)出即為出苗，在無(wú)人機(jī)獲取大豆田間圖像時(shí)，人工對(duì)區(qū)域內(nèi)大豆進(jìn)行統(tǒng)計(jì)，以獲得大豆苗數(shù)，根據(jù)式(8)出苗統(tǒng)計(jì)誤差公式計(jì)算實(shí)測(cè)苗數(shù)和估算苗數(shù)誤差，實(shí)測(cè)苗數(shù)和估算苗數(shù)試驗(yàn)結(jié)果如表2所示。

(8)

式中：W——誤差；

R——實(shí)際測(cè)得的苗數(shù)；

A——估算的苗數(shù)。

通過(guò)計(jì)算得出實(shí)測(cè)苗數(shù)與估算苗數(shù)的平均誤差為0.43%，根據(jù)表2所示試驗(yàn)結(jié)果，通過(guò)Origin 2019b軟件構(gòu)建如圖11所示的線(xiàn)性回歸模型。由線(xiàn)性回歸模型得出R2為0.909 4，線(xiàn)性回歸模型擬合度相對(duì)較高，能準(zhǔn)確地通過(guò)無(wú)人機(jī)遙感技術(shù)結(jié)合苗數(shù)估算方法進(jìn)行大豆苗數(shù)估算。

表2 試驗(yàn)結(jié)果Tab. 2 Results of experiments

圖11 實(shí)測(cè)苗數(shù)與估算苗數(shù)的比較Fig. 11 Comparison of measured and estimated seedling numbers

3 結(jié)論

以苗期獲取的大豆無(wú)人機(jī)圖像為研究對(duì)象，最終選擇Otsu閾值算法結(jié)合過(guò)綠指數(shù)分割大豆圖像，根據(jù)空間分辨率剔除圖像中異常值，通過(guò)形態(tài)學(xué)腐蝕、膨脹操作去除雜草噪聲后采用連通區(qū)域分析估算大豆苗數(shù)。

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

2) 基于無(wú)人機(jī)遙感技術(shù)的估算苗數(shù)與田間實(shí)測(cè)苗數(shù)建立的線(xiàn)性回歸模型中，相關(guān)系數(shù)R2為0.909 4，由出苗統(tǒng)計(jì)誤差公式計(jì)算的平均誤差為0.43%，提出的苗數(shù)估算方法誤差較低，能夠快速準(zhǔn)確識(shí)別大豆苗數(shù)。

3) 研究可快速評(píng)估大豆最終田間出苗，為田間大豆出苗率的測(cè)定提供了有力的支持，為后續(xù)對(duì)播種效果和質(zhì)量評(píng)估工作奠定了基礎(chǔ)，為農(nóng)業(yè)生產(chǎn)者進(jìn)行田間管理提供參考決策意見(jiàn)。