吳煥麗 崔可旺, 張 馨 薛緒掌 鄭文剛 王 巖
(1.河北工業(yè)大學(xué)電子信息工程學(xué)院, 天津 300401; 2.北京農(nóng)業(yè)信息技術(shù)研究中心, 北京 100097)
植被覆蓋度是反映作物個體/群體動態(tài)變化的重要參考量,是一種衡量作物生長狀況的主要指標(biāo)[1-2],并與植物的光合有效吸收(APAR)相關(guān),在一定程度上反映了植物的光合和蒸騰作用[3]。WANJURA等[4]認為遙感植被指數(shù)受覆蓋度的影響大于干/鮮生物量或葉面積指數(shù)(LAI)等其他參數(shù)。因此,估測農(nóng)作物的覆蓋度在農(nóng)業(yè)生產(chǎn)與科研上有非常重要的意義。
方雨晨等[5]研究了不同覆蓋度下小麥農(nóng)田土壤對歸一化植被指數(shù)(NDVI)的影響。王春梅等[6]研究了葉面積指數(shù)和覆蓋度的時空變異性。對于農(nóng)作物覆蓋度的估測,有以下幾種方法[7-9]:①有經(jīng)驗的農(nóng)學(xué)家依據(jù)肉眼進行估測,缺點是這種方法具有很強的主觀性,不同的人估測,結(jié)果可能會有較大的差別。②在正午前后太陽角度較高時,利用直尺等儀器進行測量,缺點是易受天氣因素影響,并且勞動強度過大。③采用相機或攝像頭垂直拍攝圖像進行逐點判斷,優(yōu)點是準(zhǔn)確度很高,可實現(xiàn)在線或遠程監(jiān)測,缺點是信息采集受環(huán)境影響,依賴圖像處理算法。近年來,作物表型研究在預(yù)估作物發(fā)育狀況方面成為研究熱點,具有動態(tài)監(jiān)測作物性狀的特點,其中無損圖像分析處理已成為獲取植物覆蓋度重要手段。
獲取植被覆蓋度最關(guān)鍵是圖像背景分割,準(zhǔn)確地從背景中分割圖像是減少表型分析誤差的關(guān)鍵一步。有一些算法可以對植物葉片進行分類和分割,LIU等[10]提出了一種使用三次分水嶺標(biāo)記的自動葉提取算法。DEEPIKA等[11]提出了一種花卉植物分割方法,分別應(yīng)用梯度特征的金字塔直方圖和金字塔局部二元模式提取花卉植物的形狀和紋理特征。國內(nèi)許多學(xué)者也對背景分割做了相關(guān)研究。萬園潔等[12]提出基于改進的全卷積網(wǎng)絡(luò)圖像語義分割方法,改善小麥圖像分割中過分割以及欠分割現(xiàn)象。孫立新等[13]提出一種改進果蠅算法優(yōu)化模糊均值聚類算法的圖像分割算法。黃芬等[14]針對經(jīng)典圖像分割算法精度低、過分割等問題提出了基于HSI空間H分量的K-means算法,對小麥冠層圖像進行分割。何建斌等[15]采用K-means算法和數(shù)學(xué)形態(tài)學(xué)相結(jié)合的方法進行分割,充分利用了小麥植株顏色和背景顏色的差異。上述研究分割對象多為大葉植物或作物種植密度較為稀疏的植物,多是從復(fù)雜背景中將作物目標(biāo)分離出來,很少對拍攝對象進行連續(xù)試驗。
為克服細葉作物分割精度低、圖像分割結(jié)果未得到應(yīng)用的問題,本文提出一種融合小波分析按比例萎縮去噪算法、改進果蠅算法以及K-means算法的圖像分割算法,以獲取小麥覆蓋度,通過多幅圖像的試驗對其性能進行檢驗,并通過連續(xù)試驗探索覆蓋度與作物系數(shù)之間的關(guān)系。
數(shù)據(jù)采集地點為北京市昌平區(qū)小湯山國家精準(zhǔn)農(nóng)業(yè)研究基地,選取24套小型蒸滲儀上種植的小麥進行圖像采集,如圖1所示。試驗時間為2017年9月30日—2018年6月5日。采集點共24個小區(qū),編號為1~24,每4個小區(qū)為一組循環(huán),每個小區(qū)面積為0.7 m×1 m,分別按照參比蒸散的30%、60%、90%以及大水漫灌的方式進行灌溉,選擇大水漫灌區(qū)域(圖1中紅色區(qū)域)的小麥進行圖像數(shù)據(jù)采集和處理,編號分別為4、8、16、20、24。
圖1 數(shù)據(jù)采集地點示意圖Fig.1 Sketch of data collection site1.農(nóng)大212品種小麥 2.中麥1062品種小麥 3.小麥保護帶 4.地下小型蒸滲儀 5.蒸滲儀地下監(jiān)測室
在2018年3—5月期間,共拍攝10次,使用Canon EOS 7D相機(5 184像素×3 456像素)和DM10型紅外熱像儀(160像素×120像素),垂直于地面放置,相機距離地面1.3 m、DM10型紅外熱像儀距離地面1.5 m處進行拍攝。為降低相機幾何畸變給圖像帶來的影響,在拍攝時盡量避免使用鏡頭焦距的最廣角端或最遠攝端,并使用較小的光圈。每隔7 d對5塊小麥地分別在09:30、11:00、13:00、14:30、16:00 5個時刻拍攝5組數(shù)據(jù),每次重復(fù)拍攝3幅圖像以消除系統(tǒng)誤差,每塊區(qū)域共采集150幅高像素圖像以及150幅低像素圖像,均以JPEG格式存儲到計算機中。因天氣條件的變化,大部分圖像是在晴天拍攝,也有部分是陰天與有風(fēng)天氣下拍攝。在Matlab R2017b試驗平臺實現(xiàn)本文圖像分割算法對數(shù)據(jù)的處理。
圖像預(yù)處理是整個圖像處理的關(guān)鍵部分,對圖像分割和圖像分析都有著至關(guān)重要的作用。由于待處理圖像的采集環(huán)境是自然環(huán)境,受光照影響較大,針對這種狀況,本文圖像預(yù)處理的重點為圖像光照分布。
選擇恰當(dāng)?shù)牟噬臻g是進行有效分割的基礎(chǔ)[16-17]。而面向設(shè)備的RGB顏色空間由于其分量間高度線性相關(guān),使其不適用于直接分割彩色圖像。由于在HSV色彩空間中的各分量值固定在[0,1]區(qū)間內(nèi),有利于目標(biāo)對象特征參數(shù)的提取[18],可以避免果蠅算法隨機搜索到范圍之外。為便于果蠅算法的計算,本文將圖像由RGB色彩空間轉(zhuǎn)換到HSV色彩空間,實現(xiàn)分割空間的選擇,通過Matlab函數(shù)庫中的rgb2hsv以及hsv2rgb函數(shù)實現(xiàn)圖像在RGB與HSV色彩空間之間的轉(zhuǎn)換。
K-means算法的基本思想是以空間中k個點為中心進行聚類,對最靠近他們的對象進行歸類。通過迭代法逐次更新各聚類中心的值,直至得到最優(yōu)的聚類結(jié)果[19]。
在K-means算法中,可以人為輸入?yún)?shù)k,然后將輸入的n個數(shù)據(jù)對象劃分為k個聚類,以便使所獲得聚類滿足同一聚類中的對象相似度較高;而不同聚類中的對象相似度較低[20-22]。假設(shè)把樣本集分為c個類別,算法描述如下:
(1)適當(dāng)選擇c個類的初始中心。
(2)在第k次迭代中,對任意一個樣本求其到c個聚類中心的距離,將該樣本歸到距離最短的中心所在的類。
(3)利用均值法更新該類的中心值。
(4)對于所有的c個聚類中心,重復(fù)步驟(2)、(3),如果聚類中心值保持不變,則迭代結(jié)束,否則繼續(xù)迭代。
此算法存在的主要問題是易陷入局部最優(yōu),因此對其進行改進。
使用K-means算法作為基礎(chǔ)算法,提出兩點改進:①去噪和改進的自適應(yīng)步長果蠅算法解決局部最優(yōu)問題,得到最優(yōu)化的初始聚類中心。②用Matlab函數(shù)庫中的kmeans函數(shù)完成圖像的K-means算法分割,并將之與本文算法進行比較。自適應(yīng)果蠅K均值聚類(IFOA-K-means)算法流程圖如圖2所示。
圖2 算法流程圖Fig.2 Flow chart of algorithm
IFOA-K-means算法步驟如下:
(1)首先對圖像進行預(yù)處理。為降低異常點對后續(xù)K-means算法的影響,選擇小波比例萎縮去噪方法進行圖像去噪處理,并將圖像色彩空間由RGB空間轉(zhuǎn)換為HSV空間。
(2)初始化聚類中心選擇。利用自適應(yīng)步長果蠅算法找到最優(yōu)化的初始聚類中心。
(3)分割圖像。利用K-means算法對小麥HSV彩色圖像進行聚類分割,當(dāng)準(zhǔn)則函數(shù)收斂時,迭代中止,調(diào)用Matlab函數(shù)庫中的kmeans函數(shù)實現(xiàn)。
(4)邊緣處理。利用形態(tài)學(xué)的開閉操作對分割圖像的灰度圖像進行校正,使圖像分割效果更好。
圖3 分割算法流程圖Fig.3 Flow chart of split algorithm
(5)將本文小麥圖像分割算法與基于遺傳算法的最大類間方差法[23-24]以及傳統(tǒng)K-means算法做對比。分割算法流程圖如圖3所示。
由噪聲小波系數(shù)[25]特性可知,如果能夠根據(jù)小波系數(shù)的局部特征及時地調(diào)整去噪處理方法,使其具有很強的自適應(yīng)性,去噪效果將會得到很大改進。而局部調(diào)整的極限就是對每個小波系數(shù)都采用不同的處理方法,小波比例萎縮法的自適應(yīng)性就很好地做到了這一點。本文所處理的圖像環(huán)境是自然環(huán)境[26],受光照影響較大,因此本文提出只針對V分量做去噪處理,具體步驟如下:
(1)對加噪后圖像y(i,j)(設(shè)圖像尺寸為m×n)做正交小波變換得到小波系數(shù)Y(i,j),1≤i≤m,1≤j≤n。
(2)選用正方形窗口,把所估計的小波系數(shù)Y(i,j)放在窗的中央。窗口尺寸可以選擇3×3、5×5、7×7等。2(i,j)的計算公式為
(1)
式中Ω(i,j)——包含Y(i,j)的窗口內(nèi)小波系數(shù),窗口尺寸為M×M
在對邊緣系數(shù)估計時,還需要根據(jù)窗口尺寸對邊緣進行延拓,常用的方法是對稱延拓。
(2)
其中
median(·)——中值函數(shù)
(3)
圖4 去噪算法效果圖Fig.4 Effect diagrams of denoising algorithm
對采集到的圖像進行V分量去噪和H、S、V3個分量同時去噪,由處理效果圖可以看出,單獨對V分量去噪比同時對H、S、V3個分量去噪平滑程度大,且沒有出現(xiàn)新的噪點,而同時對3個分量去噪導(dǎo)致圖像上出現(xiàn)了不規(guī)則的黑色噪點,主要集中在紅色記號筆圈出區(qū)域,小波按比例萎縮去噪算法模型無法同時使3個分量的去噪效果都達到最佳,因而在去噪后仍有部分噪點,影響之后的圖像背景分割效果。
果蠅算法是一種基于果蠅覓食行為推演出的尋求全局優(yōu)化的新方法[27-29]。在基本果蠅算法中,迭代尋優(yōu)采用固定的步長,不利于算法的收斂和穩(wěn)定,如果步長取值過小,會降低算法的收斂速率,也會使得收斂的精度降低,但步長過長又會導(dǎo)致果蠅算法跳過最優(yōu)解,穩(wěn)定性降低,極易出現(xiàn)振蕩等問題,因而本文對果蠅算法的步長進行了改進。
為提高果蠅算法的收斂速率、精度以及穩(wěn)定性能,本文提出一種根據(jù)聚類誤差自適應(yīng)調(diào)整步長的果蠅算法,根據(jù)聚類誤差的大小來調(diào)節(jié)下一次迭代的步長,調(diào)整原則如下:①如果當(dāng)前迭代得到的聚類中心在圖像上的聚類誤差(果蠅的味道濃度判定值)越小,則步長應(yīng)該越小,所以這里使用當(dāng)前的最小味道濃度(loss)作為下一次迭代的參考量。②如果當(dāng)前迭代的loss大于上次迭代的loss,則表明上次迭代得到的聚類中心是一個比較優(yōu)的值(可能是局部最優(yōu),也可能是全局最優(yōu)),單純從這一點來看,后面的迭代應(yīng)該減小步長。
原則①可以保證在loss較大時快速進行大范圍搜索,便于找到理想的極小值點,因此有助于跳出不理想的局部極值點。原則②可以保證在理想的局部極值點處有更精確的收斂(即使用更小的步長)。
進入迭代尋優(yōu),迭代步長根據(jù)兩個自適應(yīng)步長迭代原則變換,如果最佳味道濃度低于前一迭代最佳味道濃度[30],或者當(dāng)前迭代次數(shù)大于最大迭代次數(shù),則終止迭代。本文自適應(yīng)果蠅K均值算法與基本果蠅算法運行對比,如圖5所示。
圖5 果蠅算法對比Fig.5 Comparison chart of fruit fly algorithm
由圖5可知,本文自適應(yīng)果蠅K均值算法比基本果蠅算法的收斂速率更大,濃度判定值更小,即聚類誤差更小。
5.1.1分割效果
為驗證本文算法的性能,將本文算法與傳統(tǒng)的K-means算法以及基于遺傳算法的最大類間方差法作對比。
(1)不同時期分割效果
為方便看出分割效果的區(qū)別,使用紅色記號筆對上述圖像相同位置做了標(biāo)記,圖6~8分別為3月29日09:30、4月11日09:30、5月3日09:30拍攝的4區(qū)小麥圖像。由圖6d、7d可以看出,在小麥返青期和拔節(jié)期利用最大類間方差法雖然分出了小麥的大致形態(tài),但由于小麥初期有小部分泛黃區(qū)域,因而在選擇閾值進行分割時,將不屬于作物的部分也劃分為前景。該算法易受強烈光照的影響,如圖8d所示右下角部分小麥未被識別出。在小麥返青期和拔節(jié)期利用K-means算法可以在一定程度上分割出小麥圖像,但是在小麥較為密集的地方出現(xiàn)模糊現(xiàn)象,而且在小麥抽穗期,分割失真情況也比較嚴(yán)重;而利用本文算法分割小麥圖像,在麥底和麥尖處的分割都較為清晰,邊緣也得到了較好的處理,避免了K-means算法處理圖像時出現(xiàn)的模糊現(xiàn)象。
圖6 返青期3種分割算法效果Fig.6 Effect diagrams of three segmentation algorithms for return period of wheat
圖7 拔節(jié)期3種分割算法效果Fig.7 Effect diagrams of three segmentation algorithms for jointing stage of wheat
圖8 抽穗期3種分割算法效果Fig.8 Effect diagrams of three segmentation algorithms for heading stage of wheat
(2)不同天氣分割效果
選取晴天(712.6 W/m2)、陰天(18.0 W/m2)、有風(fēng)天氣(7.4 m/s)下小麥圖像,圖9為4月16日13:00拍攝的4區(qū)小麥圖像,圖10為4月17日16:00拍攝的8區(qū)小麥圖像,圖11為4月10日09:30拍攝的16區(qū)小麥圖像,分別使用K-means、IFOA-K-means、最大類間方差法3種算法對圖像做背景分割。使用最大類間方差法處理時,陰天效果最佳,晴天出現(xiàn)強烈光照或有較明顯陰影區(qū)域?qū)Ψ指钚Ч泻艽蟮挠绊?,如圖9d所示。使用K-means算法處理時,均可在不同天氣下將小麥背景去除,但是在很多葉子的細節(jié)部分,出現(xiàn)黑色孔洞,而IFOA-K-means算法在不同天氣下對細節(jié)部分處理的效果均較好,如圖9c、10c、11c紅色記號筆圈出區(qū)域。
圖9 晴天3種分割算法效果Fig.9 Effect diagrams of three segmentation algorithms under sunny day
圖10 陰天3種分割算法效果Fig.10 Effect diagrams of three segmentation algorithms under cloudy day
圖11 有風(fēng)天氣3種分割算法效果Fig.11 Effect diagrams of three segmentation algorithms under windy day
通過上述不同時期、不同天氣下的分割效果圖可以看出,最大類間方差法在處理圖像時易受到天氣的影響,并且在小麥生長初期階段,由于小麥、土壤、灌水帶之間的分割閾值不易確定,很容易將不是小麥的部分劃分為前景。而K-means算法相對來說處理效果比較穩(wěn)定,只是在有些細節(jié)部分的處理效果圖有孔洞或模糊現(xiàn)象;本文算法在處理不同時期、不同天氣的圖像時都有較好的穩(wěn)定性以及分割效果,小麥葉子細節(jié)部分也較為清晰。
5.1.2不同尺寸圖像峰值信噪比數(shù)據(jù)對比
由于最大類間方差法易受小麥生育期、天氣等因素的影響,因而后續(xù)分析主要針對本文算法以及K-means算法。峰值信噪比(PSNR)[31]是評價圖像的客觀指標(biāo),值越大,說明分割圖像越逼真。本文計算了拔節(jié)期兩種不同尺寸圖像用不同算法分割的運行時間、信噪比,結(jié)果如表1所示。
表1 不同算法性能評估結(jié)果Tab.1 Performance evaluation result of different algorithms
由表1可以看出,圖像尺寸會直接影響算法運行的結(jié)果。對于同一算法,低像素圖像由于總的像素點少,在運行時間上比高像素圖像縮短了70%以上,但是在峰值信噪比上,比高像素圖像降低了40%以上。本文算法在處理高像素圖像時,運行時間比K-means算法縮短了53%,峰值信噪比提高了21%,在對低像素圖像進行處理時,由于像素較低導(dǎo)致了迭代速率減小,所以本文算法運行時間略高于傳統(tǒng)算法,但是峰值信噪比提高了20%。
植被覆蓋度指包括喬木、灌木、草和農(nóng)作物在內(nèi)的所有植被的冠層、枝葉在生長區(qū)域地面的垂直投影面積占研究統(tǒng)計區(qū)域面積的百分比[32]。圖像采集時為降低圖像幾何畸變應(yīng)盡量避免使用鏡頭焦距的最廣角端或最遠攝端,并使用較小的光圈。經(jīng)圖像處理后將小麥的圖像背景去除,利用分割后的灰度圖像計算小麥像素點數(shù)占總像素點數(shù)的比例,即小麥垂直投影面積占總面積的比例,所求結(jié)果便是研究區(qū)域的小麥覆蓋度。
5.2.1不同方法提取的小麥覆蓋度
選取小麥具有代表性的3個生長時期圖像數(shù)據(jù),通過采用AutoCAD軟件對圖像中的綠葉進行實際勾繪計算覆蓋度[33-34],將其作為小麥覆蓋度真值與IFOA-K-means、K-means算法處理結(jié)果做對比,選擇農(nóng)大212品種(16區(qū))和中麥1062品種(4區(qū))各1個小區(qū),每個小區(qū)每個生長時期選取正午圖像各2幅,對比結(jié)果如表2所示。
表2 不同算法處理的小麥覆蓋度對比Tab.2 Comparison of wheat coverage with different algorithms %
注:同行數(shù)據(jù)后不同字母表示處理間差異顯著。
由表2可知,IFOA-K-means算法與覆蓋度真值無顯著性差異,K-means算法與覆蓋度真值存在顯著性差異,本文算法結(jié)果更接近于覆蓋度真值,準(zhǔn)確率在90%以上。
5.2.2基于IFOA-K-means算法的覆蓋度結(jié)果
對4區(qū)和16區(qū)小麥覆蓋度結(jié)果進一步分析,根據(jù)2個小區(qū)日覆蓋度均值做出覆蓋度變化趨勢圖,如圖12所示,詳細數(shù)據(jù)如表3所示。2個小麥品種的增長趨勢大致相同,拔節(jié)期生長迅速,并在5月中旬前后,覆蓋度達到最大值后開始下降。返青期初期,由于4區(qū)小麥較16區(qū)小麥距遮雨棚更近,生長時光照強度低于16區(qū)小麥,導(dǎo)致初期中麥1062品種小麥(4區(qū))覆蓋度低于農(nóng)大212品種小麥(16區(qū));而后期,4區(qū)小麥距水泵更近,可能水量較16區(qū)稍多,因而在灌水后4區(qū)小麥生長速率超過16區(qū)??梢姡庹蘸退侄紩π←湼采w度的變化產(chǎn)生影響。
圖12 日覆蓋度變化趨勢Fig.12 Changing trend of canopy coverage
區(qū)域3月22日3月29日4月3日4月11日4月16日平均值標(biāo)準(zhǔn)差平均值標(biāo)準(zhǔn)差平均值標(biāo)準(zhǔn)差平均值標(biāo)準(zhǔn)差平均值標(biāo)準(zhǔn)差4區(qū)25.652.0030.601.7748.890.5663.471.5669.421.4216區(qū)31.881.7934.221.7240.812.0262.860.37區(qū)域4月25日5月3日5月10日5月18日5月23日平均值標(biāo)準(zhǔn)差平均值標(biāo)準(zhǔn)差平均值標(biāo)準(zhǔn)差平均值標(biāo)準(zhǔn)差平均值標(biāo)準(zhǔn)差4區(qū)82.470.9585.741.2388.750.5481.800.8776.350.8516區(qū)83.602.1485.932.3687.721.8285.321.1277.214.64
5.2.3日覆蓋度、蒸滲儀、氣象數(shù)據(jù)聯(lián)合分析
根據(jù)圖像獲取的日覆蓋度數(shù)據(jù)、蒸滲儀獲取的實際蒸散量數(shù)據(jù)、氣象數(shù)據(jù)計算得到的參比蒸散量數(shù)據(jù)可得到日覆蓋度與作物系數(shù)之間的關(guān)系,即由圖像數(shù)據(jù)便可推測作物系數(shù)。
根據(jù)蒸滲儀數(shù)據(jù)以及氣象數(shù)據(jù)得到作物系數(shù)為
Kc=ETc/ET0
(4)
式中ETc——由蒸滲儀土體電壓數(shù)據(jù)計算得到的作物實際蒸散量
ET0——由氣象數(shù)據(jù)計算得到的參比蒸散量
Kc——作物系數(shù)
由3—5月獲取的5個小區(qū)圖像數(shù)據(jù)計算得到的日覆蓋度與對應(yīng)時間計算得到的作物系數(shù)進行擬合,得到兩者之間的關(guān)系
f(x)=12.97x3-18.17x2+8.537x-0.818 1
(5)
式中x——日覆蓋度f(x)——作物系數(shù)
決定系數(shù)R2為0.953 1,擬合曲線如圖13所示。
圖13 日覆蓋度與作物系數(shù)關(guān)系曲線Fig.13 Relationship curve of canopy coverage and crop coefficient
根據(jù)式(5)便可通過圖像數(shù)據(jù)推測作物系數(shù)。日覆蓋度在20%~50%時,曲線平緩,作物系數(shù)隨覆蓋度的增長緩慢上升,在50%~90%之間時,作物系數(shù)隨覆蓋度的增長而快速上升,即拔節(jié)期之后,小麥生長比較迅速,也是灌水的最佳時期,作物管理者也需在這個時期多加注意小麥覆蓋度的變化,及時做出相應(yīng)處理。
(1)從不同天氣、不同時期的去除效果、算法運行速率以及覆蓋度的計算等方面,對最大類間方差法、K-means算法以及本文提出的IFOA-K-means算法作對比,得到本文算法均優(yōu)于前2個算法,計算覆蓋度準(zhǔn)確率在90%以上。
(2)光照和水分都會對小麥覆蓋度產(chǎn)生不同程度的影響,在返青期,光照充足會使小麥迅速生長;在拔節(jié)后期,水分對小麥覆蓋度影響較大,如果缺水,將嚴(yán)重影響小麥的品質(zhì)和產(chǎn)量。因而,在小麥生長過程中,要兼顧這兩個因素。
(3)小麥覆蓋度和作物系數(shù)具有相似的變化規(guī)律,前期由小變大,成熟期時達到最大值后逐漸減小。利用Matlab軟件平臺建立由覆蓋度預(yù)測作物系數(shù)的模型,擬合方程的決定系數(shù)R2為0.953 1。