基于稻穗幾何形態(tài)模式識別的在穗籽粒數(shù)估測

馬志宏， 貢 亮， 林 可， 毛雨晗， 吳 偉， 劉成良

(上海交通大學 機械與動力工程學院, 上海 200240)

隨著人口的增加，水稻產(chǎn)量的提高已成為農(nóng)業(yè)研究的重點問題.稻穗是水稻的結(jié)實器官，與水稻產(chǎn)量直接相關(guān)，其表型特征是水稻育種[1-2]研究中重要的量化指標，主要包括每穗籽粒數(shù)、一次枝梗數(shù)量及長度、一次枝梗籽粒數(shù)、結(jié)實率、稻穗骨架結(jié)構(gòu)[3]等.目前這些特征數(shù)據(jù)主要是通過傳統(tǒng)手工測量得到，這種方法不僅勞動強度大，而且效率低下，受人為主觀性限制，重復性較差.隨著技術(shù)的發(fā)展，自動化、高通量等技術(shù)運用于作物表型特征的提取，逐漸取代了傳統(tǒng)手工測量，可以節(jié)省大量時間和人力成本.圖像處理技術(shù)可以從高分辨率圖像自動測得各種表型參數(shù)[4-5].針對稻穗籽粒數(shù)的測量，國內(nèi)外科研人員已經(jīng)開發(fā)了一些基于圖像處理的自動化設(shè)備[6-7]，這些設(shè)備均需要先對稻穗脫粒，再結(jié)合圖像處理方法對稻穗籽粒數(shù)和其他性狀參數(shù)進行測量，會破壞稻穗形態(tài)結(jié)構(gòu)，不利于樣本保存.

基于圖像分割的方法是密集顆粒狀物體計數(shù)的常用方法，在細胞計數(shù)[8-9]上頗具效果.但是對于稻穗平面圖像，由于谷粒間互相遮擋和重疊較嚴重，谷粒之間邊緣不夠明顯，令計數(shù)難以進行，國內(nèi)外相關(guān)研究也多采用稻穗面積等特征來預(yù)測谷粒個數(shù).P-TRAP[10]是一款用于水稻表型測量的開源軟件，它采用總面積除以單個谷粒平均面積，并乘以一個修正系數(shù)加以修正的結(jié)果，作為籽粒數(shù)的測量結(jié)果，這種方法需要提前測量非重疊谷粒的平均面積，而且修正系數(shù)由人為設(shè)定，主觀因素干擾大.文獻[11]中采用X射線成像技術(shù)對水稻穗部籽粒數(shù)無損檢測， X射線具有良好的透射特性，可以較好地解決遮擋問題，但這種方法具有放射性污染.文獻[12]中采用稻穗投影面積和一次枝?？傞L度進行回歸預(yù)測，籽粒數(shù)預(yù)測誤差為 7.9%，但只能預(yù)測整株稻穗籽粒數(shù)，不能預(yù)測每個一次枝梗所含籽粒數(shù).

除了圖像處理技術(shù)，機器學習技術(shù)也越來越多地應(yīng)用到植物表型提取和測量研究領(lǐng)域[13-14].但是，利用機器學習結(jié)合圖像形態(tài)學特征來預(yù)測稻穗谷粒數(shù)方面的相關(guān)研究尚處于空白階段.因此，本文提出采用圖像特征提取和機器學習方法，訓練模型進行稻穗谷粒數(shù)預(yù)測.具體做法是，利用圖像處理方法提取稻穗一次枝梗，基于支持向量機(Support Vector Machine，SVM)訓練一次枝梗谷粒數(shù)預(yù)測模型，預(yù)測結(jié)果求和作為整株稻穗的谷粒數(shù).本文認為二次枝梗為一次枝梗的一部分，且二次枝梗通常較小，使用圖像處理算法提取二次枝梗難度較大，在預(yù)測籽粒數(shù)時不再區(qū)分一次枝梗和二次枝梗，這樣對總粒數(shù)測量并沒有影響.本文以3個品種的稻穗為研究對象，分別為粳稻日本晴、9522、秈稻9311，并通過實驗驗證了所提方法的有效性.

圖1 稻穗谷粒數(shù)預(yù)測方法流程圖Fig.1 Flow chart of rice grain number prediction method

1 研究方法

由于稻穗谷粒之間遮擋和粘連情況較嚴重，而且2個粘連谷粒之間沒有明顯的邊界特征，所以基于圖像分割方法的谷粒計數(shù)難以得到高精度的結(jié)果.如果采用機器學習的方法，直接使用圖像處理方法從稻穗圖像中提取特征，訓練粒數(shù)預(yù)測模型來計算稻穗谷粒個數(shù)，則可以避免圖像分割，且隨著訓練樣本的不斷擴充，理論上模型預(yù)測精度和泛化能力可以提高.文獻[13]中采用稻穗投影面積和一次枝?？傞L度進行回歸預(yù)測，因此可以使用稻穗投影面積和一次枝梗長度作為訓練特征的一部分.選擇訓練對象時，如果將整株稻穗作為訓練對象，則存在以下問題：如果將個數(shù)預(yù)測看作分類問題，其類別會超過200個，當?shù)舅氲臉颖緮?shù)量較少時，模型容易出現(xiàn)過擬合.實際中，采集并標注超過500個稻穗樣本需要耗費很大的勞動力.考慮到每株稻穗平均包含若干個一次枝梗，每個一次枝梗所包含的谷粒數(shù)大部分在40粒以內(nèi)，如果將一次枝梗作為訓練對象，則可以增加訓練樣本數(shù)量，減少類別數(shù)目，有利于對抗過擬合.因此，本文提出如圖1所示的稻穗谷粒數(shù)預(yù)測方法，首先進行一次枝梗的提取，對每個一次枝梗分別提取特征并用SVM模型預(yù)測谷粒個數(shù)，隨后將所有預(yù)測值相加即可得到稻穗谷?？倲?shù).

2 水稻稻穗形態(tài)圖像特征提取

特征提取是機器學習的前提，也是決定模型預(yù)測精度的關(guān)鍵，如果特征選擇不妥當，比如存在過多相關(guān)性差的冗余特征，訓練結(jié)果可能出現(xiàn)過擬合等問題.觀察稻穗樣本可以看出，顏色特征跟谷粒個數(shù)幾乎沒有相關(guān)性，而形態(tài)特征則有顯著相關(guān)性.另外，顏色特征容易受到光照環(huán)境影響，而形態(tài)特征與光照環(huán)境無關(guān)，不容易受到外界環(huán)境的干擾.因此，本文只采用形態(tài)學特征訓練模型.

2.1 圖像預(yù)處理

提取形態(tài)學特征前，需要將圖像轉(zhuǎn)換為二值圖，使一次枝梗的前景像素值為1，背景的像素值為0.由于光照情況下一次枝梗周圍會產(chǎn)生陰影，為了克服陰影帶來的影響，分割出完整的一次枝梗前景，這里采用在HSV(色度、飽和度和純度)顏色空間進行二值化操作.H、S分量可以消除亮度信息的影響，考慮本研究中由于光照不均勻，一次枝梗周圍會產(chǎn)生陰影，影響一次枝梗前景二值化，故可采用S通道抑制陰影.利用OSTU自適應(yīng)閾值的方法對S通道圖像二值化，再通過開運算進行形態(tài)學濾波，去除二值圖中的噪點，同時消除一次支梗頭尾部多余的細枝，得到二值圖(見圖2(b)).隨后計算距離變換圖，距離變換是計算并標識空間點(對目標點)距離的過程，它最終把二值圖像變換為灰度圖像(其中每個像素的灰度值等于它到最近目標點的距離)，使用歐式距離變換方式，對二值圖進行處理得到圖2(c)，圖中顏色越接近黑色，表示離最近點距離越大，越接近白色，表示離最近點距離越小.它能反映圖像的某些形態(tài)學特征，如果配合骨架信息，能得到一次枝梗沿軸線上不同部位谷粒的分布密度.為此，先提取二值圖的骨架，使用細化算法進行骨架提取，得到骨架圖，骨架圖與距離變化圖相乘得到圖2(d)，記為骨架距離圖，由圖2(d)可以看出一次枝梗沿軸線上不同部位谷粒的分布密度.

圖2 一次枝梗形態(tài)學特征圖Fig.2 Morphological features of primary branches

具體實現(xiàn)過程使用基于Python 3.5語言編寫的程序.使用OpenCV(版本2.4.9)工具進行HSV顏色空間轉(zhuǎn)換、二值化、開運算、輪廓提取、細化、距離變換等操作.其中，二值化閾值采用函數(shù)Threshold計算，方式參數(shù)為OSTU，形態(tài)學操作的結(jié)構(gòu)元素設(shè)置半徑為3的圓盤結(jié)構(gòu)，距離變換方式為歐氏距離，細化操作使用Hilditch算法[15]實現(xiàn).

2.2 特征提取

對于一次枝梗形態(tài)學特征，最直觀的是投影面積和一次枝梗長度特征[12]，在此基礎(chǔ)上，又提取了二值圖前景的輪廓和最小外接矩形，計算出周長，最小外接矩形長、寬、長寬比，矩形度(區(qū)域面積與最小外接矩形面積之比)，骨架曲線長度這些特征.另外，根據(jù)一次枝梗的形狀特點，設(shè)計了骨架距離平均值和骨架距離方差，前者計算沿骨架上每一點的距離值平均值，容易發(fā)現(xiàn)，這個平均值一定程度上反映了一次枝梗的平均半徑.設(shè)二值圖為B(x,y)，距離變換圖為D(x,y)，骨架圖為T(x,y)，骨架距離圖為H(x,y)，則H(i,j)=D(i,j)T(i,j)，i=1,2,…,w，j=1,2,…,h，其中w是圖像寬度，h是圖像高度.特征計算公式如下：

(1) 投影面積.投影面積可以直接從一次枝梗二值(圖2(b))反映，這里以一個像素點代表單位面積，投影面積則是對二值圖求和

(1)

(2) 周長.周長即輪廓長度，設(shè)一個閉合輪廓坐標點序列為{(x1,y1),(x2,y2),…,(xn,yn)},記xn+1=x1,yn+1=y1，其中n為序列長度，則周長為

(2)

(3) 最小外接矩形長寬比，主軸長度.提取最小外接矩形，得到其最小外接矩形2個邊長，記為l1，l2，則長寬比為

f3=max{l1/l2,l2/l1}

(3)

使用最小外接矩形長邊作為主軸長度

f4=max{l1,l2}

(4)

(4) 矩形度.矩形度表示區(qū)域面積與其最小外接矩形面積的比值

f5=f1/(l1l2)

(5)

(6)

(6) 骨架距離均值.遍歷骨架上點，求骨架點所對應(yīng)的距離值的均值，公式為

(7)

以上所有特征都具有旋轉(zhuǎn)不變性，平移不變性，但是某些特征不具有縮放不變性.因此需要保證訓練樣本的像素比和測試樣本的像素比相同，否則同一個一次枝梗在不同像素比的圖像中，提取的特征可能會有較大差異.本文實驗均為同一平臺，所得到圖像的像素比均為 0.118 mm/像素.而實際應(yīng)用過程中，圖像可能來自其他成像平臺，對于這些圖像，提取特征前需要先將圖像等比例縮放到固定的像素比，使得像素比的值等于訓練集圖像的像素比 0.118 mm/像素，才能用本文算法得到正確的結(jié)果.

3 籽粒數(shù)預(yù)測SVM模型

籽粒數(shù)預(yù)測的目標是一個非負整數(shù)，在有限樣本下，這個非負整數(shù)有上界U，此時粒數(shù)預(yù)測可以看成一個分類問題，其類別數(shù)為U個.本文算法的模式識別部分選擇作為分類器， SVM是建立在統(tǒng)計學習理論VC維理論和結(jié)構(gòu)風險最小化原理基礎(chǔ)上的機器學習方法.它在解決小樣本、高維和非線性模式識別問題中表現(xiàn)出許多獨特的優(yōu)勢，而且具有堅實的理論基礎(chǔ)[16].

通常，一個SVM基本分類器只能完成二分類任務(wù)，對于對分類任務(wù)，一般采用多個SVM集成預(yù)測，通常采用的訓練策略有一對余型SVM、一對一型SVM等，如果類別數(shù)為n，一對余型SVM需要訓練n-1個SVM基本分類器，一對一型SVM需要訓練n(n-1)/2個SVM基本分類器.考慮到粒數(shù)預(yù)測值為連續(xù)整數(shù)的特點，本研究設(shè)計了一種集成SVM預(yù)測模型訓練策略.設(shè)訓練集中一次枝梗粒數(shù)的最大值為k粒，最小值為1粒，則設(shè)計這樣一個分類器，它由(k-1)個二分類SVM組成，第i個SVM的輸出表示粒數(shù)是否超過i，其中i=1,2,…,k-1.訓練第i個SVM分類器Fi時，正樣本對應(yīng)大于i粒的一次枝梗，負樣本對應(yīng)少于或等于i粒的一次枝梗，則一次枝梗谷粒數(shù)預(yù)測值為

(8)

相比于一對余型SVM的訓練策略，上述方法需要的SVM個數(shù)同樣為n-1.不同的是，只有對位于兩頭的SVM訓練時，訓練集樣本不平衡問題才會比較顯著，所以這種策略相對一對余型SVM可以更好地解決樣本集不平衡問題，相對于一對一型SVM有預(yù)測速度的優(yōu)勢.

訓練對象是一次枝梗圖像(見圖2(a))，需要將每個一次枝梗摘下，分別拍照并手工計數(shù)標記.以9個一次枝梗為一組，將其按3×3的九宮格位置均勻擺放在背景板上，背景板尺寸為300 mm×240 mm，每個小方格可容納的一次枝梗長度為100 mm，可容納大多數(shù)長度的一次枝梗.使用500萬像素的工業(yè)相機，調(diào)整相機高度和焦距，使相機的視野剛好完整包含背景板，隨后進行一次拍照，采集圖像.然后，翻轉(zhuǎn)每個一次枝梗，再進行一次拍照，這樣，總共得到2張照片，故一組一次枝?？梢缘玫?8個訓練樣本，在擺放枝梗的同時，手工數(shù)出籽粒個數(shù)，按圖像名記錄在表格，隨后經(jīng)過圖像裁剪，得到所有一次枝梗圖像的訓練樣本，根據(jù)表格標記對應(yīng)的籽粒數(shù).實驗所用稻穗為3個不同的品種混合而成，分別為粳稻日本晴、9522、秈稻9311.實驗中總共采集了 2 688 幅稻穗一次枝梗圖像，裁剪后分辨率為800 pixel×600 pixel，每幅一次枝梗圖像僅包含1個一次枝梗.采集所用的相機為中國大恒公司生產(chǎn)的MER-500-7UM/UC-L型工業(yè)數(shù)字相機，實驗軟件基于Python 3.5，計算機的配置為處理器 Intel(R) Core(TM)i7-6700@3.40 GHz，內(nèi)存8 GB.

訓練時，采用徑向基核函數(shù)，將 2 688 個樣本分為訓練集和測試集，大小分別為 2 190 幅，498幅，最大一次枝梗谷粒數(shù)為26粒，最小為1粒，故k=26.在訓練時，使用交叉驗證策略，將訓練集分為5份，其中 1 752 幅用來訓練，438幅用來驗證，以交叉驗證的平均預(yù)測誤差為優(yōu)化目標，尋找最優(yōu)懲罰參數(shù)C和核函數(shù)參數(shù)σ.

圖3 稻穗形態(tài)特征圖和一次枝梗提取效果Fig.3 Extraction effect of morphological features of panicle and primary branches

為保留有效特征，減少不重要特征數(shù)量，避免過擬合，本文使用ReliefF算法[17]先對特征進行重要性分析和排序，并根據(jù)特征權(quán)重，選擇若干組特征組合分別進行訓練，特征權(quán)重和不同特征組合的訓練結(jié)果如表1所示.由表1可知，誤差最小的2組特征分別為組合4(f1,f2,f4,f5,f6)和組合5(f1,f2,f3,f4,f5,f6)，兩者誤差相等，但前者比后者少一個特征f3，可以減少一點計算量，故最終選擇前者作為訓練特征.

表1 不同特征組合的預(yù)測結(jié)果Tab.1 Prediction results of different feature combinations

4 一次枝梗提取和測量

本文的訓練對象是一次枝梗圖像，而實際中測量對象往往是整株稻穗.稻穗在自然狀態(tài)下具有收縮內(nèi)斂的趨勢，直接成像往往會導致大量谷粒被遮擋，因此實際中，通常需要將稻穗平展鋪開再進行成像，如圖3(a)所示.為了提取稻穗上所有的一次枝梗，需要剔除莖稈和芒，文獻[10,12]中直接進行形態(tài)學開操作的方式去除直徑較細的莖稈和芒，開操作首先對二值圖進行腐蝕運算，再進行膨脹運算，而腐蝕運算會將所有半徑較小的區(qū)域剔除，半徑較大

的區(qū)域則保留，目的性較差.開操作結(jié)構(gòu)元素半徑過小，不能剔除較粗的莖稈，半徑過大容易過度腐蝕，導致一次枝梗斷裂，如圖3(c)所示.所以采用這種方法通常需要二次判別篩選，將非目標區(qū)域剔除.本文根據(jù)稻穗形狀特點，提出基于局部距離方差的一次枝梗區(qū)域提取方法.這里，局部距離方差指以稻穗骨架上一點為中心，一定長寬的圓形領(lǐng)域內(nèi)，所有骨架點在距離變換圖上的距離均方差，如圖4所示.

ω(i,j,r)=sgn(r2-i2-j2)

(9)

(10)

E(x,y,r)=

(11)

E(x,y,r)T(x,y)]2

(12)

式中：x,y為骨架點坐標；r為鄰域半徑；H(x,y)為骨架距離圖(見圖3(e))；ω(i,j,r)表示圓形模板，當i,j在以(0,0)為中心，半徑為r的圓內(nèi)時，輸出為1，其余為0；E(x,y,r)為骨架點半徑為r的領(lǐng)域內(nèi)距離均值；S(x,y,r)為骨架點半徑為r的領(lǐng)域內(nèi)距離方差.

圖4 局部距離方差示意圖Fig.4 Schematic of local distance variance

經(jīng)過實驗，取閾值α，以r半徑的鄰域內(nèi)距離方差小于α，并且當前點距離值小于β時，將這點周圍以1.5倍的距離值為鄰域的區(qū)域置為背景像素值，否則，不操作.如圖4所示，(x,y)表示當前需要處理的像素點坐標，點(x,y)的局部距離方差是指，在以此點為圓點，半徑為r的圓內(nèi)，計算所有骨架上的像素點對應(yīng)的距離值的方差，圖4中灰色區(qū)域為稻穗前景，細線條狀像素點表示骨架，其顏色代表像素點對應(yīng)在距離變換圖中的距離值，越接近黑色，表示距離值越小，越接近白色，表示距離值越大.當遍歷完整個骨架像素點時，所有非目標物都被剔除，如圖3(f)所示，后面采用連通域查找算法即可提取所有一次枝梗的二值圖.經(jīng)過實驗，α取6.5，β取10，r取40時效果最佳.

算法基于局部距離方差的一次枝梗分離算法

Input: image,α,β,r#輸入稻穗圖像和參數(shù)

1:B=binarization(image) #圖像預(yù)處理和二值化

2:D=distanceTransform(B) #距離變換

3:T=thining(B) #細化

4:H=DT#矩陣D和矩陣T對應(yīng)元素相乘，得到矩陣H

5: Forxfrom 1 to width

6: Foryfrom 1 to height

7: IfH[x][y] > 0 andH[x][y]<βthen

8:s=S(x,y,r) #計算以(x,y)為中心的局部距離方差，見式(12)

9: Ifs<αthen

10: drawCircle(B,x,y,1.5×H[x][y])

#在圖B上以點(x,y)為中心,1.5×H[x][y]為半徑畫圓，填充灰度為0

11: End if

12: End if

13: End for

14: End for

Output:B

由于這一步已經(jīng)進行了距離變換和細化操作，后面提取一次枝梗形態(tài)學特征時不必重復做這2個步驟，一定程度上節(jié)省了時間.所以，對于整株稻穗的谷粒計數(shù)算法偽代碼如上述算法所示，通過二值化，形態(tài)學預(yù)處理，基于局部距離方差的一次枝梗區(qū)域提取，對每個一次枝梗提取特征，并預(yù)測籽粒個數(shù)，最后對所有預(yù)測結(jié)果求和，得出整株稻穗的籽粒預(yù)測數(shù).從圖3(f)可以看出，本文算法可以有效去除支稈，完整保留一次枝梗，并且不會造成一次枝梗斷裂(見圖3(c))，可以保證一次枝梗形態(tài)學特征提取的準確度，提高一次枝梗谷粒數(shù)的預(yù)測精度.

5 實驗結(jié)果與分析

實驗選取62株完整的稻穗來測試本文算法的預(yù)測效果，這些稻穗樣本品種來源與訓練樣本相同，不同之處為訓練樣本是將稻穗剪切成若干一次枝梗，而測試樣本則為完整稻穗.每株稻穗包含的一次枝梗數(shù)目為7～18個，總谷粒數(shù)在40～200粒.稻穗鋪開后在同一實驗平臺(即采集訓練樣本的平臺)下進行圖像采集，并手工數(shù)粒進行標記.

為了體現(xiàn)算法的預(yù)測穩(wěn)定性，將測試樣本隨機分成3組，每組約20個稻穗，記為樣本1， 2， 3.實驗時，分別對3批樣本單獨測試1遍，然后對3批樣本合起來測試1遍.表2顯示了測試結(jié)果，可以看出最大平均誤差為 6.72%，優(yōu)于文獻[12]的測量誤差.為了體現(xiàn)本文算法的優(yōu)越性，同時使用一對余型SVM進行訓練和預(yù)測.由表2可見，一對余型SVM所得到的誤差略高于本文算法，說明本文提出的SVM訓練策略較好地解決了數(shù)據(jù)不平衡問題.

圖5所示為一組稻穗測試樣本的谷粒數(shù)預(yù)測結(jié)果和一次枝梗提取效果.結(jié)果顯示，大部分稻穗一次枝梗均能準確提取，但是當2個一次枝梗之間沒有充分分開，有粘連情況，則會將2個一次枝梗誤判為1個一次枝梗，因此稻穗擺放時需要保證每個一次枝?；ゲ徽尺B.

表2 稻穗谷粒數(shù)預(yù)測實驗結(jié)果Tab.2 Experimental results of panicle grain number prediction

圖5 稻穗測試樣本的谷粒數(shù)預(yù)測結(jié)果和一次枝梗提取效果Fig.5 The prediction results and primary branch extraction effects of rice grain numbers of test samples

6 結(jié)語

本文提出的使用形態(tài)學特征進行模型訓練來預(yù)測谷粒數(shù)的方法，在實驗中最大平均誤差為 6.72%，說明本文方法可以成功應(yīng)用在稻穗谷粒數(shù)預(yù)測上；而形態(tài)學特征對顏色和光照環(huán)境不敏感，魯棒性較高，若稻穗發(fā)生霉變或其他原因造成稻穗顏色變化，只要形狀不發(fā)生變化，依然能得到相同測量結(jié)果；此外，基于局部距離方差的一次枝梗提取方法能夠較好地避免一次枝梗斷裂的問題，保證了一次枝梗特征提取的準確度和穩(wěn)定性.不足之處是必須人工將稻穗充分鋪開，否則會出現(xiàn)較大誤差.同時，人工鋪開稻穗耗費時間和體力，因此本文方法只適合小批量樣本的測量.后續(xù)工作中，稻穗的鋪開方式還有進一步優(yōu)化空間，采用自動化或半自動化方式鋪開稻穗，可以減少人力勞動成本，提高本文方法的實用性.