基于NDWI和卷積神經(jīng)網(wǎng)絡(luò)的冬小麥產(chǎn)量估測方法

劉峻明 周 舟 和曉彤 王鵬新 黃健熙

(1.中國農(nóng)業(yè)大學(xué)土地科學(xué)與技術(shù)學(xué)院， 北京 100083； 2.中國農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院， 北京 100083)

0 引言

冬小麥?zhǔn)鞘澜缛蠹Z食作物之一，準(zhǔn)確預(yù)測冬小麥單產(chǎn)及其空間分布對保障國家糧食安全和挖掘區(qū)域可利用的農(nóng)業(yè)資源具有重要意義。目前利用遙感數(shù)據(jù)進(jìn)行冬小麥產(chǎn)量預(yù)測主要有基于作物敏感波段反射率的統(tǒng)計回歸方法[1]、基于遙感信息和作物模型的數(shù)據(jù)同化方法[2]等。統(tǒng)計回歸方法中應(yīng)用最廣泛的是植被指數(shù)法，植被指數(shù)是由多光譜數(shù)據(jù)經(jīng)線性或非線性組合構(gòu)成的對植被有一定指示意義的數(shù)值，被廣泛用于植被長勢監(jiān)測和作物估產(chǎn)領(lǐng)域。然而采用植被指數(shù)的經(jīng)驗?zāi)Ｐ屯ǔＪ腔谀骋蛔魑镌谔囟l件下建立，例如歸一化差值植被指數(shù)(Normalized difference vegetation index, NDVI)在植被達(dá)到一定覆蓋度后增長緩慢，在高植被覆蓋地區(qū)存在過飽和現(xiàn)象，往往導(dǎo)致模型缺乏普適性。數(shù)據(jù)同化技術(shù)的思想最早由文獻(xiàn)[3]提出并在作物估產(chǎn)領(lǐng)域得到應(yīng)用。通過數(shù)據(jù)同化技術(shù)使遙感數(shù)據(jù)參與到作物模型模擬過程中，可以提升模型在區(qū)域尺度的應(yīng)用精度，但作物模型的參數(shù)標(biāo)定以及數(shù)據(jù)同化系統(tǒng)的運行效率仍是當(dāng)前應(yīng)用的難點，在實際應(yīng)用中很難快速得到大范圍產(chǎn)量預(yù)測結(jié)果且預(yù)測精度還有待提高[4]。

隨著機器學(xué)習(xí)技術(shù)不斷發(fā)展，機器學(xué)習(xí)和深度學(xué)習(xí)方法被成功應(yīng)用到多個領(lǐng)域，如圖像識別、機器翻譯、信號處理等[5]。傳統(tǒng)的機器學(xué)習(xí)方法如支持向量機、決策樹及隨機森林等能夠較好地解決非線性問題，并且有較好的效果[6-9]。在作物估產(chǎn)領(lǐng)域，研究已表明[10-12]深度學(xué)習(xí)方法通常能夠得到比傳統(tǒng)機器學(xué)習(xí)方法更高的精度。其中，卷積神經(jīng)網(wǎng)絡(luò)是使用最廣泛的深度學(xué)習(xí)方法。相較于其他神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，需要的參數(shù)相對較少，且可直接將多維圖像數(shù)據(jù)作為網(wǎng)絡(luò)輸入特征，信息損失小。目前在作物產(chǎn)量預(yù)測方面，文獻(xiàn)[13]利用卷積神經(jīng)網(wǎng)絡(luò)和遙感數(shù)據(jù)對玉米產(chǎn)量進(jìn)行回歸預(yù)測，相比于支持向量回歸算法更有優(yōu)勢并取得了較好的預(yù)測結(jié)果。文獻(xiàn)[14]提出了一種新的卷積神經(jīng)網(wǎng)絡(luò)估產(chǎn)框架，使用遙感影像的直方圖信息作為模型輸入，既達(dá)到了數(shù)據(jù)降維的效果，同時也取得較高的產(chǎn)量預(yù)測精度，此方法具有遷移學(xué)習(xí)的能力[15]，在產(chǎn)量數(shù)據(jù)較少的地區(qū)同樣有較好的效果，為基于卷積神經(jīng)網(wǎng)絡(luò)的估產(chǎn)方法在區(qū)域擴(kuò)展性上提供了參考。

歸一化差值水指數(shù)(Normalized difference water index, NDWI)是基于MODIS地表反射率2個近紅外波段計算得到，它能快速有效地提取植被冠層的水分含量，并及時地響應(yīng)植被冠層受水分脅迫的影響[16]。本文基于MODIS數(shù)據(jù)，以遙感植被特征直方圖分布信息作為輸入變量，應(yīng)用卷積神經(jīng)網(wǎng)絡(luò)對冬小麥產(chǎn)量進(jìn)行回歸預(yù)測，并通過分析不同遙感植被特征在冬小麥產(chǎn)量估測上的表現(xiàn)，探討NDWI在區(qū)域冬小麥產(chǎn)量估測上的應(yīng)用，同時結(jié)合田間實測霜凍害資料進(jìn)一步探究模型在霜凍害發(fā)生年份的表現(xiàn)，以期為縣域尺度農(nóng)作物產(chǎn)量預(yù)測提供有效方法。

1 研究區(qū)與實驗數(shù)據(jù)

1.1 研究區(qū)域

河南省地處31°23′0″～36°22′0″N，110°22′0″～116°38′0″E(圖1)，屬北亞熱帶濕潤氣候與暖溫帶半濕潤季風(fēng)氣候間的過渡氣候，年平均降水量為500～1 000 mm，研究區(qū)冬小麥多種植冬性、弱冬性品種，一般在9月中旬至10月上旬播種，翌年5月底至6月初成熟[17]。獨特的地理位置和氣候條件使其成為我國冬小麥霜凍害高發(fā)區(qū)之一，霜凍害年際發(fā)生率高達(dá)60%[18]，嚴(yán)重地區(qū)可減產(chǎn)60%～70%[19]，對冬小麥的高產(chǎn)穩(wěn)產(chǎn)構(gòu)成威脅。

1.2 數(shù)據(jù)來源與處理

采用MODIS 8d合成地表反射率產(chǎn)品(MOD09A1)作為遙感數(shù)據(jù)源，分別提取其紅光波段、近紅外波段、藍(lán)光波段、綠光波段的地表反射率數(shù)據(jù)作為遙感輸入特征，空間分辨率為500 m，時間分辨率為8 d(2004-10-08—2017-06-18)，應(yīng)用MODIS重投影工具軟件(Modis reprojection tool, MRT)對影像的研究區(qū)域進(jìn)行裁剪，投影轉(zhuǎn)換為WGS84坐標(biāo)系統(tǒng)，并利用MODIS土地覆蓋類型產(chǎn)品(MCD12Q1)植物功能型分類方案(Land cover type 5)識別冬小麥種植區(qū)。產(chǎn)量數(shù)據(jù)來自《河南省農(nóng)村統(tǒng)計年鑒》[20]中縣級冬小麥單產(chǎn)數(shù)據(jù)，在2005—2017年間連續(xù)種植冬小麥的縣市共有103個，共獲得1 339個有效單產(chǎn)數(shù)據(jù)。冬小麥?zhǔn)転?zāi)、成災(zāi)、絕收面積等霜凍害資料來自于農(nóng)業(yè)農(nóng)村部種植業(yè)管理司歷年自然災(zāi)害數(shù)據(jù)庫(http:∥sjcx.fldj.agri.cn/moazzys/zaiqing.aspx)。

2 研究方法

2.1 遙感植被指數(shù)計算

遙感植被指數(shù)的選取主要圍繞植物光合作用和冠層水分條件，因此對于這2個參數(shù)敏感的MODIS可見光、近紅外波段都被納入到候選波段中，其中包括計算NDVI所需的紅光波段、短波近紅外波段，以及計算NDWI所需的短波近紅外波段、長波近紅外波段。通過將提取的可見光、近紅外波段的反射率數(shù)據(jù)作為輸入特征，分別計算得到NDVI[21]、NDWI[16]、綠紅植被指數(shù)(Green red vegetation index, GRVI)[22]、綠色歸一化植被指數(shù)(Green normalized difference vegetation index, GNDVI)[23]、調(diào)整土壤亮度植被指數(shù)(Optimal soil adjusted vegetation index, OSAVI)[24]、土壤調(diào)節(jié)植被指數(shù)(Soil adjusted vegetation index, SAVI)[25]、改進(jìn)型土壤調(diào)節(jié)植被指數(shù)(Modified soil-adjusted vegetation index, MSAVI)[26]8 d間隔的時間序列，時間序列長度為32，覆蓋了冬小麥的整個生育周期(10月至次年6月中旬)。

2.2 直方圖信息提取

因各縣域遙感影像形狀、像元數(shù)量差異較大，為實現(xiàn)樣本結(jié)構(gòu)的標(biāo)準(zhǔn)化，統(tǒng)計各縣域不同遙感植被指數(shù)在冬小麥生育期時間序列上的像素分布直方圖作為樣本特征。NDVI、GRVI、GNDVI、OSAVI、MSAVI的像素值范圍為(0，1)，NDWI的范圍為(-0.25，0.5)，SAVI的范圍為(0，1.5)，使用最大最小值法將NDWI和SAVI歸一化至0到1之間，計算式為

式中Inorm——歸一化后NDWI或SAVI值

I——NDWI或SAVI值

Imin——NDWI或SAVI的最小值

Imax——NDWI或SAVI的最大值

將所有植被指數(shù)值等間隔劃分至32個區(qū)間內(nèi)，統(tǒng)計各區(qū)間像素百分比，得到頻率直方圖。

2.3 單產(chǎn)去趨勢處理

2003—2017年河南省小麥產(chǎn)量總體呈增長趨勢(圖2)，主要得益于小麥品種的改進(jìn)、管理技術(shù)的提高和農(nóng)業(yè)政策的改革[27-29]，為去除由技術(shù)進(jìn)步、經(jīng)濟(jì)社會發(fā)展等因素引起的冬小麥單產(chǎn)變化，本文采用5 a滑動平均法計算得到其趨勢單產(chǎn)，用實際單產(chǎn)減去趨勢單產(chǎn)對研究區(qū)單產(chǎn)數(shù)據(jù)進(jìn)行去趨勢處理。

2.4 卷積神經(jīng)網(wǎng)絡(luò)

卷積神經(jīng)網(wǎng)絡(luò)是一種前饋型的神經(jīng)網(wǎng)絡(luò)，其主要組成結(jié)構(gòu)包括輸入層、卷積層、池化層、全連接層、激活函數(shù)等。以研究區(qū)逐年各縣的不同遙感植被指數(shù)的像素直方圖及其對應(yīng)的縣域冬小麥單產(chǎn)作為模型的樣本數(shù)據(jù)集，采用卷積神經(jīng)網(wǎng)絡(luò)模型進(jìn)行訓(xùn)練和驗證。本文的卷積神經(jīng)網(wǎng)絡(luò)輸入層為32×32×7的矩陣，各卷積層的卷積核個數(shù)依次是128、128、256、256、512、512、512，卷積核尺寸都是3×3，滑動步長分別為1、2、1、2、1、1、2。同時，在每一個卷積層上進(jìn)行批歸一化和線性整流函數(shù)(ReLU)激活操作，并在全連接層加入隨機失活(Dropout)操作。具體參數(shù)設(shè)置為：使用方差縮放方法初始化網(wǎng)絡(luò)權(quán)重，偏差初始化為0，初始學(xué)習(xí)率為0.001，Drouput設(shè)置為0.5，使用自適應(yīng)矩估計(Adam)優(yōu)化器，運行一次輸入的樣本數(shù)為32。

2.5 模型精度評價與驗證

采用決定系數(shù)(Coefficient of determination，R2)、均方根誤差(Root mean square error，RMSE)和平均絕對誤差(Mean absolute error，MAE)3個指標(biāo)對模型擬合程度優(yōu)劣進(jìn)行評價[30]。

3 結(jié)果與分析

3.1 不同植被特征估產(chǎn)精度比較

為評估NDWI在冬小麥產(chǎn)量估測上的表現(xiàn)，利用經(jīng)遙感植被指數(shù)計算直方圖信息提取得到的NDVI、NDWI、GRVI、GNDVI、OSAVI、SAVI、MSAVI 7個遙感特征，以2005—2014年1 030組數(shù)據(jù)作為訓(xùn)練樣本，分別將7個遙感特征的直方圖信息作為特征集，實際單產(chǎn)作為目標(biāo)變量構(gòu)建基于CNN的回歸預(yù)測模型，以2015—2017年309組數(shù)據(jù)作為驗證樣本輸入模型，對比分析不同遙感植被指數(shù)對預(yù)測精度的影響并對去趨勢前后模型預(yù)測精度進(jìn)行比較，結(jié)果見表1。從表1中可以看出，相對于植被指數(shù)NDVI、SAVI、OSAVI、GNDVI、MSAVI、GRVI，NDWI表現(xiàn)出更好的預(yù)測效果，單產(chǎn)去趨勢前后的NDWI對產(chǎn)量的預(yù)測精度均高于NDVI、SAVI等植被指數(shù)。

表1 不同植被指數(shù)去趨勢前后驗證結(jié)果Tab.1 Verification results of different vegetation indexes

在冬小麥主要生育期，水分是限制冬小麥葉片生長的重要因素，通過影響冬小麥的光合作用來限制其子房發(fā)育灌漿后期，莖葉中的營養(yǎng)物質(zhì)通過植株體內(nèi)的水運輸?shù)阶蚜Ｖ校謺苯佑绊憼I養(yǎng)物質(zhì)的輸送，進(jìn)而影響冬小麥產(chǎn)量[31]。相比于去趨勢前模型預(yù)測結(jié)果，NDWI在單產(chǎn)去趨勢后R2提高了0.05，實測單產(chǎn)和預(yù)測單產(chǎn)的散點圖如圖3所示，趨勢線與1∶1線交于0點附近，大部分樣本聚集在1∶1線周圍，R2最高達(dá)到0.79，MAE和RMSE分別為482、637 kg/hm2，主要是由于去除了產(chǎn)量年際間社會經(jīng)濟(jì)因素的影響。

3.2 不同生育階段NDWI對產(chǎn)量的影響

為進(jìn)一步分析不同生育階段NDWI對產(chǎn)量的影響，將冬小麥全生育期劃分為6個時間段，分別為10月8日—11月25日、12月3日—2月26日、3月6日—3月30日、4月7日—4月30日、5月1日—5月17日、5月25日—6月18日，大致對應(yīng)冬小麥出苗—越冬、越冬—返青、返青—拔節(jié)、抽穗—灌漿、乳熟—成熟階段，分別以各時間段NDWI作為樣本特征輸入，去趨勢單產(chǎn)為目標(biāo)變量，驗證結(jié)果見圖4。從圖4中可以看出，抽穗—灌漿階段模型預(yù)測精度最優(yōu)，MAE和RMSE分別為552 kg/hm2和759 kg/hm2，R2最高達(dá)到0.74，說明該階段影像反映的植被狀況對產(chǎn)量的影響最大，乳熟—成熟以及返青—拔節(jié)階段次之，出苗—越冬階段預(yù)測效果相對較差。冬小麥在返青—拔節(jié)階段主要進(jìn)行營養(yǎng)生長，該階段是決定穗數(shù)和粒數(shù)的關(guān)鍵時期，但其生長特征并不能完全反映產(chǎn)量形成器官的干物質(zhì)積累過程[32]，因此該生育階段的模型精度較低。抽穗—灌漿階段有機物從營養(yǎng)器官轉(zhuǎn)移到籽粒，該階段NDWI與冬小麥千粒質(zhì)量密切相關(guān)[33]，故此階段估產(chǎn)精度最高。在冬小麥乳熟—成熟階段，冠層和莖稈的營養(yǎng)物質(zhì)向籽粒轉(zhuǎn)移，葉片中的葉綠素含量下降，與產(chǎn)量的相關(guān)性變?nèi)鮗34]，故在成熟后期估產(chǎn)模型精度下降。

為確定研究區(qū)冬小麥產(chǎn)量估測的最佳時間，根據(jù)NDWI在不同生育階段的模型驗證結(jié)果，對抽穗—灌漿階段進(jìn)一步劃分為3個時間段，分別為4月7—14日、4月15—23日、4月23—30日，分別對各時間段的NDWI進(jìn)行訓(xùn)練并預(yù)測2015—2017年對應(yīng)時段的冬小麥產(chǎn)量，結(jié)果見表2。從表2中可以看出，4月23—30日的NDWI對產(chǎn)量的決定系數(shù)可達(dá)到0.72，MAE和RMSE分別為566、763 kg/hm2。這主要是由于籽粒最終產(chǎn)量主要來源于抽穗—成熟階段葉片的光合產(chǎn)物，而灌漿后期是籽粒干物質(zhì)積累最旺盛的時期，地上干生物量中籽粒比重較大[35]。

表2 抽穗—灌漿生育階段驗證結(jié)果Tab.2 Verification results of NDWI in heading—filling growth stage

3.3 模型預(yù)測誤差空間分布

對NDWI模型精度進(jìn)行逐年驗證，2015—2017年的精度驗證結(jié)果如圖5所示，每一年份分為4幅子圖，包括單產(chǎn)估測值與實測值散點圖、縣域?qū)崪y單產(chǎn)分布、模型估測單產(chǎn)分布以及誤差分布。從空間分布上來看，模型估測單產(chǎn)與實測單產(chǎn)圖中高產(chǎn)區(qū)和低產(chǎn)區(qū)分布基本一致，東部單產(chǎn)最高，中部次之，西部最低，整體呈東高西低。從誤差分布圖中可以看出，大部分區(qū)域誤差在±300 kg/hm2內(nèi)，估測誤差大于900 kg/hm2的縣主要分布在西部和北部山區(qū)與東部黃淮海平原交界處，這些縣的單產(chǎn)較低，對應(yīng)于散點圖中的低產(chǎn)部分，低產(chǎn)縣的模型估測單產(chǎn)普遍高于實測單產(chǎn)。造成此誤差的主要原因有：①低產(chǎn)區(qū)的單產(chǎn)樣本數(shù)量較少，從散點圖中可以看出，中產(chǎn)和高產(chǎn)數(shù)據(jù)的密度較大，而低產(chǎn)數(shù)據(jù)則相對較少，而樣本分布不平衡是造成機器學(xué)習(xí)和深度學(xué)習(xí)預(yù)測偏差的主要原因之一。②低產(chǎn)區(qū)域多為山區(qū)-平原過渡地帶，地勢西高東低，地形較為復(fù)雜，農(nóng)田小氣候與地形因素會對冬小麥產(chǎn)量產(chǎn)生一定的影響，從而影響到模型估測效果。

4 討論

4.1 霜凍害影響下NDWI模型預(yù)測精度

考慮到河南省霜凍害發(fā)生次數(shù)較頻繁，造成小麥減產(chǎn)，有效地預(yù)測霜凍害影響下的冬小麥產(chǎn)量，對于冬小麥災(zāi)害預(yù)警、穩(wěn)產(chǎn)高產(chǎn)具有重要意義。為探究NDWI模型在霜凍害影響下的預(yù)測效果，使用留一年法對模型精度進(jìn)行逐年驗證，結(jié)合2005—2017年農(nóng)業(yè)農(nóng)村部種植業(yè)管理司歷年自然災(zāi)害數(shù)據(jù)庫河南省霜凍害資料，驗證結(jié)果如圖6所示。由圖6可以看出，模型預(yù)測精度在2005—2017年整體上呈現(xiàn)波動變化的趨勢，有無霜凍發(fā)生年份均維持在一較高水平，在預(yù)測精度達(dá)到最高的2006年、2008年和2013年均有霜凍發(fā)生。霜年平均R2約為0.78，平均RMSE、MAE分別為682、527 kg/hm2，這表明模型在霜凍害影響下仍然能保持較好的預(yù)測效果。

2013年4月15—23日，河南省商丘地區(qū)發(fā)生春季霜凍事件，氣溫驟降15℃以上，此時小麥正處于籽粒形成的關(guān)鍵時期，受凍后會導(dǎo)致明顯的缺?，F(xiàn)象，對當(dāng)?shù)囟←湲a(chǎn)量造成極大影響。根據(jù)田間霜凍害調(diào)查結(jié)果可知，該年商丘地區(qū)平均穗粒數(shù)減少率達(dá)到40%左右，平均減產(chǎn)率達(dá)到39.6%左右[36]。本次春霜事件為探索NDWI在冬小麥低溫脅迫下的變化特征并進(jìn)一步驗證霜凍害影響提供了一個理想的案例。冬小麥發(fā)生凍害時，冠層含水量上升，植株體內(nèi)發(fā)生了失水情形。當(dāng)極端低溫超過了小麥的耐寒能力時，植株細(xì)胞原生質(zhì)體內(nèi)以及細(xì)胞間隙間的水分發(fā)生放熱凝固現(xiàn)象，植株體內(nèi)的水由于固結(jié)失去流動性，無法將營養(yǎng)物質(zhì)運送到各個器官，且細(xì)胞結(jié)構(gòu)發(fā)生不可逆的損壞，對水分的控制能力下降，當(dāng)白天氣溫上升時，植株體內(nèi)的結(jié)冰開始融化，水分開始外滲故造成失水[37]。圖7為2013年3月6日至5月17日商丘地區(qū)冬小麥種植區(qū)域的NDWI和NDVI的時間序列變化曲線，NDWI與NDVI整體變化趨勢相同，均呈先上升后下降的變化態(tài)勢。在返青初期，NDWI值較低，4月15日至4月23日和5月1日至5月9日間，NDWI有2次明顯的上升過程，結(jié)合氣象站提供的日最低溫度和日降水量數(shù)據(jù)，5月6—9日間有明顯降水現(xiàn)象，可以解釋在此期間NDWI數(shù)值的突然升高，而4月15—23日間商丘地區(qū)出現(xiàn)了大幅度降溫現(xiàn)象，4月21日最低氣溫達(dá)到0℃，低溫持續(xù)近3 d，同時期內(nèi)降水量較低，因此降水對于植被水分條件的影響可以忽略。在沒有降水輸入的情況下，NDWI數(shù)值的異常升高可認(rèn)為是由凍害引起的。

4.2 模型局限性

從估產(chǎn)驗證結(jié)果來看，NDWI能夠很好地反映植被最終的生長狀態(tài)，其R2與NDVI相比，提高了0.06，且在霜凍害影響下仍能保持較好的預(yù)測效果，但本研究仍存在一些不足以及值得進(jìn)一步探索的地方：

(1)使用植被指數(shù)的直方圖信息作為模型輸入，模型在產(chǎn)量數(shù)據(jù)較少的地區(qū)同樣有較好的效果，為基于CNN的估產(chǎn)方法在區(qū)域上的擴(kuò)展性提供了參考。但直方圖信息的提取需要樣本區(qū)域內(nèi)有足夠多的有效像元個數(shù)，因此可能并不適用于像元尺度上的單產(chǎn)估測。

(2)采用的去趨勢方法對各年份研究區(qū)全省范圍的趨勢產(chǎn)量進(jìn)行計算，而非單獨計算各年份各縣的趨勢，而不同縣的單產(chǎn)增長趨勢存在差異，故單產(chǎn)在空間上的變異性仍然存在。

(3)不同生育階段劃分的時間長度不一致，例如越冬—返青階段時間約為90 d，而返青—拔節(jié)階段、抽穗—灌漿階段時間跨度相對較小，不足30 d，所包含的影像序列長度不一，模型存在不確定性。

5 結(jié)論

(1)NDWI能夠很好地反映植被最終的生長狀態(tài)，在冬小麥生育早期的產(chǎn)量預(yù)測上表現(xiàn)出更好的預(yù)測效果，R2最高可達(dá)0.79，MAE和RMSE分別為482、637 kg/hm2，因此適合作為冬小麥估產(chǎn)指標(biāo)。

(2)NDWI在抽穗—灌漿階段對冬小麥最終產(chǎn)量影響最大，NDWI在4月23—30日時間段內(nèi)對產(chǎn)量的決定系數(shù)可達(dá)到0.72，綜合對比可知抽穗—灌漿階段NDWI對冬小麥最終產(chǎn)量影響較大。

(3)模型預(yù)測精度在霜凍發(fā)生年份，R2最高可達(dá)0.83，這表明模型在霜凍害影響下仍能保持較好的預(yù)測效果。