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

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

0 引言

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

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

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

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

1.1 研究區(qū)域

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

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

采用MODIS 8d合成地表反射率產品(MOD09A1)作為遙感數(shù)據(jù)源，分別提取其紅光波段、近紅外波段、藍光波段、綠光波段的地表反射率數(shù)據(jù)作為遙感輸入特征，空間分辨率為500 m，時間分辨率為8 d(2004-10-08—2017-06-18)，應用MODIS重投影工具軟件(Modis reprojection tool, MRT)對影像的研究區(qū)域進行裁剪，投影轉換為WGS84坐標系統(tǒng)，并利用MODIS土地覆蓋類型產品(MCD12Q1)植物功能型分類方案(Land cover type 5)識別冬小麥種植區(qū)。產量數(shù)據(jù)來自《河南省農村統(tǒng)計年鑒》[20]中縣級冬小麥單產數(shù)據(jù)，在2005—2017年間連續(xù)種植冬小麥的縣市共有103個，共獲得1 339個有效單產數(shù)據(jù)。冬小麥受災、成災、絕收面積等霜凍害資料來自于農業(yè)農村部種植業(yè)管理司歷年自然災害數(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]、調整土壤亮度植被指數(shù)(Optimal soil adjusted vegetation index, OSAVI)[24]、土壤調節(jié)植被指數(shù)(Soil adjusted vegetation index, SAVI)[25]、改進型土壤調節(jié)植被指數(shù)(Modified soil-adjusted vegetation index, MSAVI)[26]8 d間隔的時間序列，時間序列長度為32，覆蓋了冬小麥的整個生育周期(10月至次年6月中旬)。

2.2 直方圖信息提取

因各縣域遙感影像形狀、像元數(shù)量差異較大，為實現(xià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ū)間內，統(tǒng)計各區(qū)間像素百分比，得到頻率直方圖。

2.3 單產去趨勢處理

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

2.4 卷積神經網(wǎng)絡

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

2.5 模型精度評價與驗證

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

3 結果與分析

3.1 不同植被特征估產精度比較

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

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

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

3.2 不同生育階段NDWI對產量的影響

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

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

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

3.3 模型預測誤差空間分布

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

4 討論

4.1 霜凍害影響下NDWI模型預測精度

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

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

4.2 模型局限性

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

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

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

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

5 結論

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

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

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