基于CNN-RNN網絡的中國冬小麥估產

赫曉慧，羅浩田，喬夢佳，田智慧※，周廣勝

（1. 鄭州大學地球科學與技術學院，鄭州 450001； 2. 鄭州大學信息工程學院，鄭州 450001；3. 中國氣象科學研究院鄭州大學生態氣象聯合實驗室，鄭州 450001）

0 引 言

農業是人類社會發展和穩定的基礎，而中國是人口大國，糧食安全尤為重要[1]。小麥是中國的主要糧作物之一，據統計，2020年小麥播種面積為23.38萬hm2，總產量13 425萬t，占全國糧食總產量的1/5，因此，準確、及時地預測中國小麥產量對全國糧食安全乃至世界農業發展至關重要[2-3]。

近幾十年來，國內外學者圍繞著作物產量預測、糧食安全等問題展開了深入而廣泛的研究。許多學者基于經驗統計模型[4]和作物生長模擬模型[5]等傳統的估產模型預測作物產量，這些模型通常對于某一區域的特定農作物產量預測精度較高，但很難推廣到擁有大規模種植區域以及需要進行多農作物產量預估的地區和國家[6-8]。

隨著傳感器的發展，遙感衛星可以獲得大量高質量、高時空分辨率的圖像，為大尺度農作物產量預測提供了可能。遙感衛星捕獲的光譜信息能充分反映農作物的生長狀態，以往的研究者大多傾向于從影像中提取作物相關指標來建立與作物產量之間的關系[9-11]，如歸一化差異植被指數（Normalized Difference Vegetation Index，NDVI）、增強植被指數（Enhanced Vegetation Index，EVI）等[12-14]，由于農作物生長環境復雜，植株狀態多變，僅依靠植被指數可能無法完全反映影像中農作物的生長特征。

為了充分利用遙感圖像中的光譜信息，提升農作物估產的準確性，研究者逐漸將深度學習應用到了農作物估產中。如Qi等[15-16]通過構建兩個獨立分支的卷積神經網絡（Convolutional Neural Networks，CNN）結構分別處理高分辨率RGB圖像和多光譜影像，整合了影像中的空間- 光譜信息進行估產，并發現水稻成熟期有較高精度的產量預測；Nevavuori等[17]構建了基于無人機獲取的NDVI和RGB圖像的CNNs作物產量預測模型，結果表明基于RGB圖像的CNN模型優于NDVI數據。盡管基于CNN的方法取得了較好的結果，但是這些研究僅僅基于農作物某個特定生長期的影像對產量進行預測，如Son等[18-20]發現在農作物的生育后期的預測精度最高。然而，農作物的生長是一個生物量長時間累積的過程，包括多個生育期。已有研究表明，在作物生長期內的不同時相遙感圖像中，作物信息與其產量關系的顯著性具有較大差別[21-25]。因此，綜合利用農作物生長過程中的多時相遙感影像進行估產，有助于掌握不同時間階段的農作物生長狀況，提升產量預測精度[26-30]。雖然CNN被廣泛用于從遙感影像中提取光譜特征[31]。但是由于結構的限制，CNN無法對時間序列進行建模[32-34]。因此，僅利用CNN無法從多時相影像提取中作物生長的特征。因此，如何在提取遙感影像空譜特征的基礎上，為長時間序列的遙感影像進行建模，是提升農作物產量預測的關鍵問題。

為解決上述問題，本文將全國冬小麥主產區作為研究區域，并選取研究區2001－2018年縣級冬小麥生育期的時間序列遙感影像作為數據源，同時設計了一個端到端的深度學習估產網絡從多時相遙感影像中自動學習空間-光譜-時間聯合特征。具體來說，首先將每個時間節點的冬小麥影像獨立地嵌入卷積子網模塊，在卷積子網充分提取影像中的空間-光譜特征之后，將其按時間序列傳遞到遞歸子網單元，整合冬小麥整個生育期的時間信息。本文所提出的農作物產量預測網絡旨在能夠從時間序列遙感影像中多尺度、全方位地挖掘作物生長特征與冬小麥產量之間的復雜關系，以實現大尺度的冬小麥產量預測。

1 研究區及數據

1.1 研究區

本文選取全國冬小麥主產區作為研究區域，包括14個省和自治區：河北、陜西、江蘇、浙江、安徽、山東、河南、湖北、四川、云南、陜西、甘肅、寧夏回族自治區和新疆維吾爾族自治區，共1 713個區縣。

1.2 數據源

1.2.1 MODIS影像

本文采用中等分辨率成像光譜儀（MODIS）衛星圖像，該衛星影像有較高的時間分辨率，能獲取足夠多檢測冬小麥產量的多時相影像。數據下載自Google Earth Engine（https://earthengine.google.com/）。本文主要使用3種MODIS產品數據，其時間覆蓋為2001—2018年。MOD09A1 V6產品提供了500 m和8 d的地表反射率，包括7個波段，如紅波段、藍波段和近紅外波段等，其采集的地表反射率影像可以反映作物的生長狀態和環境[35-36]；MOD11A2 V6產品可在1 km×1 km的網格中提供每8 d的晝夜地表溫度影像，晝夜地表溫度與作物冠層溫度聯系緊密，有利于進行產量預測[37]。MCD12Q1 V6產品以年為單位提供全球土地覆蓋類型，其中包括草原、農田、城市、水體等17類土地類型，其空間分辨率為500 m；數值比例表示把原始光譜值從整型轉化為浮點型的比例系數；數值范圍表示原始數據中光譜值的最大值和最小值范圍，超出數值范圍的則是無效數據。

1.2.2 產量數據

冬小麥產量數據是以區縣為單位的統計數據，其中2001－2008年的冬小麥各區縣年平均產量及播種面積數據來源于中國科學院資源學科創新平臺[38]（http://www.data.ac.cn），2009—2018年數據收集自各省、市、區縣統計年鑒。

1.3 數據預處理

為統一影像分辨率，將MOD11A2的空間分辨率重采樣為500 m，再利用MCD12Q1中農作物掩膜數據提取各區縣的農作物種植面積，由于農作物掩膜不能區分具體的農作物類型，因此需要結合作物產量數據，剔除產量較低和影像中農作物像素點較少的縣，即將縣域作物像素小于2×2的影像剔除，當影像中作物像素點太少，估產網絡難以從中提取有用信息，反而會被當成噪聲數據，進而影響網絡的性能。對MOD09A1和MOD11A2各個波段進行提取、融合，即利用GDAL庫分別讀取MOD09A1和MOD11A2影像，將其波段按順序排列，共計9個波段，包括MOD09A1的7個波段和MOD11A2的2個波段，最后輸出融合后的影像。

2 研究方法

2.1 估產樣本構建

本文選取2001—2018年的冬小麥影像數據作為樣本。一般來說，冬小麥在10月中旬播種，次年7月初收獲，由此確定影像數據的時間范圍從每年的第289天到次年的第180天，共30個時間步。冬小麥生育期包括：播種期、出苗期、分蘗期、越冬期、返青期、起身期(生物學拔節)、拔節期、孕穗期、抽穗期、開花期、灌漿期、成熟期，每個生育期約21 d，把每兩個生育期劃分一個時間節點，即每兩個生育期對應5個時間步：播種-出苗期（0～5）、分蘗-越冬期（5～10）、返青-起身期（10～15）、拔節-孕穗期（15～20）、抽穗-開花期（20～25）、灌漿-成熟期（25～30）。結合縣邊界數據，裁剪每個縣

所有時間步長t∈(0,30)的影像其中，b代表參與農作物估產的波段數，id代表區縣的國家行政區劃代碼，y代表預測年份。由于各縣級單位冬小麥種植面積不同，為統一網絡的輸入，統計影像內農作物像素點個數，取32×32作為閾值。當縣域作物像素小于閾值時，將除作物像素的其余點填充為0；當縣域作物像素大于閾值時，剔除零星分布的農作物像素，這部分像素點可能會被當成噪聲，取農作物比較集中的像素作為網絡輸入；以確保影像整體信息完整。為提升模型精度和收斂速度，對冬小麥數據集進行標準化。

式中μ表示均值，σ表示方差。計算出整個數據集的均值和方差，并應用上述公式對影像中每個像素點進行標準化。

2.2 CNN-GRU融合網絡

如圖1所示，所提出的CNN-GRU估產模型由卷積子網、遞歸子網和全連接層組成。為提取多時相影像的聯合光譜-空間-時間特征，卷積子網從每個時間步的影像中提取空間-光譜特征映射，最終得到30個時間步的特征向量。在卷積子網的基礎上，遞歸子網以卷積層的特征作為輸入，利用門控循環單元（Gated Recurrent Unit，GRU）挖掘不同層次的特征間的依賴關系，整合CNN框架部分輸出的特征向量。最后利用全連接層進行產量預測。

2.2.1 基于卷積子網的空間光譜特征提取

農作物的生長主要伴隨光譜變化，因此從遙感影像中準確地挖掘光譜信息是農作物產量預測地關鍵。同時，由于農作物分布通常呈現聚集性，農作物像素之間的空間相關性也十分重要，而CNN作為深度學習的一個重要分支，其獨特的卷積結構能提取遙感影像中豐富的光譜-空間特征[39-40]。

卷積子網包括4個卷積層、4個池化層和一個全連接層，如圖2所示。輸入數據為每個時間步的冬小麥影像。在每一張影像上利用3×3的卷積核的點積之和來提取固定大小的光譜特征，并利用窗口滑動以覆蓋整個空間維度，以提取冬小麥的空間-光譜信息。本文選擇線性整流函數（Rectified Linear Unit，ReLu）作為激活函數，不僅能加速網絡收斂，還能避免梯度爆炸或梯度消失。最后輸出特征圖表示為

式中f(·)表示ReLu激活函數，xi表示上一個卷積層的輸出，wi表示第i-1層和第i層之間的權重，bi表示貝葉斯偏置項。

通過卷積層提取特征后，為減少計算量，本文選擇最大池化函數來處理從卷積運算中獲得的特征映射結果，其中池化窗口為2×2，池化步長為2。

式中xpool表示池化層的輸出。卷積子網最終將每個時間步的冬小麥影像映射為高維的空間-光譜特征{f1,f2,…,f30}，并將其作為遞歸子網的輸入。

2.2.2 基于遞歸子網的時間相關性建模

在農作物生長過程中，其植株形體特征會隨生育期變化而改變，卷積神經網絡（CNN）僅能提取單一生育期的農作物生長特征，而通過農作物不同生育期植株狀態的變化來挖掘深層次的時序特征尤為重要。循環神經網絡（Recurrent Neural Network，RNN）可以很好地關聯上下文信息，尤其擅長針對時間序列數據建模[41-42]。但傳統RNN面臨著梯度消失的問題，難以對長時間序列數據建模。因此，本文使用GRU代替傳統RNN網絡，GRU是循環神經網絡的變體，適用于挖掘長時間序列數據的內在聯系，并且有效解決了梯度消失的問題[43]。

因此，在CNN的基礎上，本文進一步使用GRU網絡構建農作物生長過程中的時間關系。GRU由重置門rt、更新門zt和隱藏狀態ht組成，將從不同時間步獲取的高維空譜特征xt作為輸入，其中隱藏狀態ht用來處理連續的冬小麥時間特征，在考慮先前冬小麥生長狀態的情況下學習當前的特征；重置門用于控制前一時刻的隱藏單元ht1-對當前輸入xt的影響；更新門確定在當前隱藏狀態ht之前，保留多少先前時刻的狀態信息。第t∈(0,30)時間步各個門的迭代參數可表示為

式中xt為第t時間步的輸入分別表示更新門、重置門、候選狀態、隱藏狀態，σ是sigmoid激活函數，表示權重矩陣和偏執向量，☉表示矩陣運算。這種結構使GRU能夠根據冬小麥先前的生長信息，結合新的輸入進行自我更新。因此，GRU可以充分挖掘各個時間步之間的內在聯系，從而有效反映農作物生長過程中的相互關系。

2.3 loss函數和評價指標

在網絡最后，將所有的時間輸出展平成一個向量，利用一個神經元的全連接層輸出冬小麥預測產量，并在每個池化層和全連接層添加比例為0.5的dropout結構，以避免模型出現過擬合現象。本文使用自適應矩估計（adaptive moment estimation，Adam）優化器訓練網絡，并通過L2損失函數[44]來計算預測產量（preid）與真實產量（realid）之間的訓練損失值，如公式所示：

本文使用均方根誤差（RMSE）、平均絕對誤差（Mean Absolute Error，MAE）和決定系數（R2）來評估模型預測精度[45-47]。

2.4 模型訓練

本文使用python對影像數據進行預處理，在基于TensorFlow的深度學習框架下進行模型的搭建、訓練和測試，并安裝CUDA等運算平臺以搭建GPU加速環境。實驗的硬件仿真環境為Intel(R) Core(TM) i7-6900K CPU，使用NVIDIA GeForce GTX 1080 GPU對模型進行加速，顯存大小為8G。

為實現全國冬小麥主產區的產量預測，本文對CNN-GRU估產模型進行訓練。在訓練階段，以全國冬小麥主產區整個生育期的遙感影像作為訓練數據，年平均產量為標簽，年份為2001－2015年，時間步t∈(0,30)。為確保無偏性，每個批次從訓練集中隨機抽取32個樣本進行訓練。初始學習率為0.000 1，分階段逐步遞減為初始值的10倍，達到迭代次數時，訓練停止。在測試階段以2016－2018年全國冬小麥主產區的年平均產量和影像作為測試數據。

3 結果與分析

3.1 數據貢獻度分析

為評估MODIS影像中地表反射率（7個波段）和地表晝夜溫度數據對農作物估產的影響程度。本文對溫度數據和反射率數據分別進行訓練，并將反射率7個波段依次去除一個，使用其余波段進行訓練。如表1所示（以RMSE為例，單位kg/hm2），其中刪除波段1表示刪除第一個波段，僅使用第2～7個波段，以此類推，最后一列波段1～7表示使用地表反射率的所有波段。結果表明：1）只利用反射率訓練的CNN_GRU估產模型RMSE年平均值分別為841 kg/hm2，優于僅使用溫度數據的預測結果，說明反射率數據相較于溫度數據包含了更多的農作物信息，如紅光波段對綠色植被具有強吸收性，對農作物的覆蓋度、生長狀態比較敏感。2）使用反射率所有波段訓練的模型效果最好，說明反射率7個波段均對均對農作物產量預測起決定性作用。在分別刪除第1、2、7個波段后，RMSE下降了8%左右，而刪除第3個波段時，RMSE僅下降了3%，說明第1（紅光波段）、2（近紅外波段）、7（短波紅外波段）相較于其他波段對農作物估產更為重要，而第3個波段（藍光波段）主要反映植被的下墊面土壤與植被之間的差異，可能相較于其他波段不能更直接的反映作物的生長狀態。

表1 地表反射率（7個波段）和地表晝夜溫度數據對農作物估產的影響Table 1 The effect of surface reflectance (7 bands) and surface day and night temperature data on crop yield estimation kg·hm-2

3.2 估產模型性能檢測

為驗證CNN-GRU估產模型的優越性，將其與使用最廣泛的多個農作物產量預測模型進行比較，如深度學習框架下的CNN、GRU和傳統機器學習算法中的SVR、DT和RF。其中所有估產模型均使用同一個數據集。但CNN難以處理長時間序列數據，則將時間序列和波段相疊加作為CNN的輸入，即（32,32,30×9）。而GRU網絡需要將輸入數據展平為時間序列向量（30,32×32×9）。此外，由于傳統機器學習模型本身的結構限制，在輸入之前需要把輸入數據展平為一維向量，但30×32×32×9=276480維度太高，容易導致“維度災難”，因此，使用影像的NDVI作為傳統機器學習模型的輸入。各個模型的關鍵參數設置如下：

1）卷積神經網絡：包含4個卷積層、4個池化層、1個全連接層，批處理大小為32，初始學習率為0.000 1。

2）門控循環單元模型：由128個神經元的GRU和全連接層組成，批處理大小為32，初始學習率為0.001。

3）支持向量回歸：使用徑向基（RBF）作為SVR的核函數，懲罰因子C=2,5,10,100，松弛變量γ=10-1,10-2,10-3。

4）決策樹：包括樹的最大深度max_depth=3,4,5,6,7，葉子節點最少樣本數min_samples_leaf=10,15,20,30，內部節點再劃分所需最小樣本數min_samples_split=5,8,10。

5）隨機森林：包括子樹的數量num_estimators=50,80,100，每棵樹的最大深度max_depth=3,6,9。

在本節中，使用不同的估產模型分別預測2016—2018年冬小麥產量，為確保比較的公平性，所有估產模型均在相同的軟硬件環境下訓練和測試。分析表2可知：1）所提出的CNN-GRU估產模型年平均RMSE為818.3 kg/hm2，相較于CNN、GRU、SVR、RF和DT分別降低了20.13%、18.81%、29.51%、34.84%和36.57%，年平均MAE為560 kg/hm2，相較于其他模型分別降低了25.46%、22.39%、37.09%、39.95%和42.26%。表明CNN-GRU模型能準確地預測冬小麥產量。但對于部分低產量的區縣（如新疆維吾爾自治區的伊寧縣和和碩縣）和部分高產量的區縣（陜西省西安市的灞橋區和河南省南陽市的新野縣）CNN-GRU模型預測誤差偏大（如圖3所示）。產生這種誤差的原因可能是冬小麥數據集中低產量和高產量的樣本數占偏低，導致模型對于部分低產量和高產量的地區泛化能力較差，預測的誤差較高；對于低產量的新疆維吾爾自治區海拔高、環境復雜等不確定性因素導致預測精度降低；對于高產量的灞橋區和新野縣來說，地處平原，均為冬小麥氣候適宜區，種植面積較大，而模型輸入中沒有表達區域異質性的參數，導致低估了部分高產量地區的冬小麥產量。2）CNN和GRU的性能相當，年平均RMSE分別為920、996 kg/hm2，均高于所提出的CNN-GRU估產模型，表明單一的神經網絡在作物估產領域效果不佳。CNN能有效地提取影像中空間-光譜特征，但作物生長隨時間變化，農作物在不同的生育期有不同的形態特征，時間信息的缺失嚴重影響了CNN的整體性能；而GRU能整合長時間序列的作物生長信息，但不能捕捉復雜的空間關系。因此所提出的CNN-GRU估產模型可以有效地提取冬小麥生長期的空間-光譜-時間信息，并將其應用于冬小麥產量預測中。3）深度學習網絡在性能上均優于傳統的機器學習模型，在冬小麥產量預測這一復雜問題的研究中，深度學習網絡更能挖掘冬小麥影像中的時空特征。

表2 不同估產模型性能對比Table 2 Performance comparison of different yield estimation models

3.3 不同生育期冬小麥預測

本節分別對作物累計生育期的時間步（0～5，0～10，0～15，0～20，0～25，0～30）和單一生育期的時間步（0～5，＞5～10，＞10～15，＞15～20，＞20～25，＞25～30）進行單產預測，以評估CNN-GRU網絡模型預測冬小麥產量的最佳時間窗。

如圖4所示，CNN-GRU估產模型在生育期初期性能較差，2016－2018年平均RMSE和MAE分別為1 109、840.6 kg/hm2，這是由于估產模型在生育期初期缺乏足夠的作物生長信息和環境信息，但隨著時間的推移，越來越多的影像信息被整合到估產網絡中，模型的性能逐步提高。1）單一生育期的預測精度很低，相較于全生育期影像（0～30）低了25%至28%，這是由于單一生育期僅使用了部分時間步的影像數據，缺乏前后生育期冬小麥的生長狀態，不能充分發揮CNN-GRU估產網絡整合長時間序列影像時空特征的優勢。2）返青-起身期和拔節-孕穗期兩個時間步的RMSE相差最大，約10%左右，拔節-孕穗期相較于返青-起身期預測精度明顯提升。在返青-起身期時，麥苗開始生長春生第一葉，莖部第一節開始伸長但為完全露出地面。而在拔節期，冬小麥植株莖部露出地面1.5～2 cm，生長特征明顯。表明CNN-GRU模型能有效地捕捉冬小麥的各個生育期的特點及其生長狀態的變化。3）在灌漿-成熟期時，CNN-GRU估產模型精度最高，2016－2018年的平均RMSE和MAE分別為817、556 kg/hm2，R2為0.76，均優于其他各個生育期。在冬小麥生育期后期，所有的作物生長信息均被整合到CNN-GRU估產模型中，該估產模型能精確地預測全國冬小麥主產區產量。4）2016－2018年冬小麥拔節-孕穗期和抽穗-開花期估產模型預測精度相當，年平均RMSE分別為849、823 kg/hm2，僅比成熟期低3.9%和0.7%，說明CNN-GRU估產模型在冬小麥生育期的中后期已經能準確地預測其產量，即所提出的估產模型有能力提前2個月預測全國冬小麥主產區產量。

3.4 CNN-GRU模型魯棒性檢驗

為對提出的CNN-GRU估產模型進行魯棒性檢驗，將2008－2018年估產樣本逐年作為驗證集，其余年份作為訓練集進行訓練，如表3所示。結果表明，2008－2018年CNN-GRU估產模型的RMSE和MAE平均值為727、534 kg/hm2，R2為0.805。在不同訓練集和驗證集下，預測產量與統計產量值離散程度基本一致，表明CNN-GRU估產模型有較高的魯棒性。

表3 CNN-GRU模型逐年精度評估Table 3 Annual accuracy evaluation of CNN-GRU model

綜合分析發現，2010年產量預測誤差最大，RMSE和MAE分別為982、752 kg/hm2，R2僅0.617；2012年模型精度最高，RMSE和MAE分別為604、446 kg/hm2，R2為0.867。2009、2011、2013和2014年模型預測精度較高，RMSE和MAE在650、500 kg/hm2左右；2016－2018年預測效果一般，RMSE在820 kg/hm2左右。不同年份模型性能不同的主要原因有：1）在2010年，中國西南地區遭遇特大干旱，北方冬小麥主產區出現持續低溫天氣，在極端天氣情況下，統計數據極值差異較大[41]。如圖5a所示，2010年全國冬小麥主產區縣級尺度產量預測出較多離群值，而多數離群值的預測產量都大于統計產量，導致模型整體精度降低。2）2012年的RMSE和MAE分別為604、446 kg/hm2，R2為0.867。如圖5b所示，預測產量與統計產量接近，僅有零星幾個離群值。3）2016－2018年的收集的產量數據較少，可用于訓練和測試的區縣數量僅約650個，而其余每年的區縣數據約1 500個。并且隨著基因技術和種植技術的發展，近年來作物單產飛速提升，2016－2018年的低產量地區數量較少，導致數據分布不均衡。而2009、2011、2013和2014年的產量分布比較均衡，導致模型在不同年份預測的精度不同。

4 結 論

作物產量預測在農業管理中具有重要意義。近年來的研究表明，遙感是一種有效的估產方法。本文提出了一種基于多時相遙感影像的CNN-GRU估產模型，利用CNN從多時相遙感影像中提取豐富的空間-光譜特征，在此基礎上，GRU自適應學習冬小麥生育期各階段之間的時間依賴。試驗證明該模型能從多時相遙感影像中提取聯合空間-光譜-時間特征，挖掘作物生長特征與產量之間的復雜關系，并能夠對收獲前期的冬小麥產量進行準確預測。結論如下：

1）深度學習技術在產量預測上有巨大潛力，對比結果表明深度學習方法優于傳統機器學習方法。同時，與單一的深度學習網絡CNN和GRU相比，本文所提出的CNN-GRU估產模型精準預測了2016－2018年中國冬小麥主產區的產量，其精度均優于其他估產模型。說明該模型能有效地提取反映作物生長信息的空間-光譜-時間特征，解決了單一神經網絡提取作物生長特征不充分問題，進一步提高了大范圍冬小麥單產估算效率。

2）基于混合神經網絡的遙感估產模型魯棒性較高。本文對2008－2018年數據集逐年樣本訓練，所提出的CNN-GRU估產模型RMSE和MAE平均值分別為727、534 kg/hm2，R2平均值為0.805。該模型在不同數據集上預測產量與統計產量值離散程度基本一致，CNN-GRU估產模型性能穩定，預測精度高，可以作為一種高效的方法為全國提供常規可行的冬小麥產量監測和預報。

3）在冬小麥成熟期CNN-GRU估產模型性能最高，但拔節-孕穗期和抽穗-開花期同樣取得較高精度，RMSE分別為849、823 kg/hm2，R2均高于0.7，與成熟期預測結果幾乎沒有差異。說明本文提出的CNN-GRU估產模型有能力提前2個月預測縣級尺度的全國冬小麥主產區產量，為冬小麥生育早期產量預測提供了新的思路和方法。