基于CNN-RNN網(wǎng)絡(luò)的中國冬小麥估產(chǎn)

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

（1. 鄭州大學(xué)地球科學(xué)與技術(shù)學(xué)院，鄭州 450001； 2. 鄭州大學(xué)信息工程學(xué)院，鄭州 450001；3. 中國氣象科學(xué)研究院鄭州大學(xué)生態(tài)氣象聯(lián)合實(shí)驗(yàn)室，鄭州 450001）

0 引 言

農(nóng)業(yè)是人類社會(huì)發(fā)展和穩(wěn)定的基礎(chǔ)，而中國是人口大國，糧食安全尤為重要[1]。小麥?zhǔn)侵袊闹饕Z作物之一，據(jù)統(tǒng)計(jì)，2020年小麥播種面積為23.38萬hm2，總產(chǎn)量13 425萬t，占全國糧食總產(chǎn)量的1/5，因此，準(zhǔn)確、及時(shí)地預(yù)測中國小麥產(chǎn)量對(duì)全國糧食安全乃至世界農(nóng)業(yè)發(fā)展至關(guān)重要[2-3]。

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

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

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

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

1 研究區(qū)及數(shù)據(jù)

1.1 研究區(qū)

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

1.2 數(shù)據(jù)源

1.2.1 MODIS影像

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

1.2.2 產(chǎn)量數(shù)據(jù)

冬小麥產(chǎn)量數(shù)據(jù)是以區(qū)縣為單位的統(tǒng)計(jì)數(shù)據(jù)，其中2001－2008年的冬小麥各區(qū)縣年平均產(chǎn)量及播種面積數(shù)據(jù)來源于中國科學(xué)院資源學(xué)科創(chuàng)新平臺(tái)[38]（http://www.data.ac.cn），2009—2018年數(shù)據(jù)收集自各省、市、區(qū)縣統(tǒng)計(jì)年鑒。

1.3 數(shù)據(jù)預(yù)處理

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

2 研究方法

2.1 估產(chǎn)樣本構(gòu)建

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

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

式中μ表示均值，σ表示方差。計(jì)算出整個(gè)數(shù)據(jù)集的均值和方差，并應(yīng)用上述公式對(duì)影像中每個(gè)像素點(diǎn)進(jìn)行標(biāo)準(zhǔn)化。

2.2 CNN-GRU融合網(wǎng)絡(luò)

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

2.2.1 基于卷積子網(wǎng)的空間光譜特征提取

農(nóng)作物的生長主要伴隨光譜變化，因此從遙感影像中準(zhǔn)確地挖掘光譜信息是農(nóng)作物產(chǎn)量預(yù)測地關(guān)鍵。同時(shí)，由于農(nóng)作物分布通常呈現(xiàn)聚集性，農(nóng)作物像素之間的空間相關(guān)性也十分重要，而CNN作為深度學(xué)習(xí)的一個(gè)重要分支，其獨(dú)特的卷積結(jié)構(gòu)能提取遙感影像中豐富的光譜-空間特征[39-40]。

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

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

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

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

2.2.2 基于遞歸子網(wǎng)的時(shí)間相關(guān)性建模

在農(nóng)作物生長過程中，其植株形體特征會(huì)隨生育期變化而改變，卷積神經(jīng)網(wǎng)絡(luò)（CNN）僅能提取單一生育期的農(nóng)作物生長特征，而通過農(nóng)作物不同生育期植株?duì)顟B(tài)的變化來挖掘深層次的時(shí)序特征尤為重要。循環(huán)神經(jīng)網(wǎng)絡(luò)（Recurrent Neural Network，RNN）可以很好地關(guān)聯(lián)上下文信息，尤其擅長針對(duì)時(shí)間序列數(shù)據(jù)建模[41-42]。但傳統(tǒng)RNN面臨著梯度消失的問題，難以對(duì)長時(shí)間序列數(shù)據(jù)建模。因此，本文使用GRU代替?zhèn)鹘y(tǒng)RNN網(wǎng)絡(luò)，GRU是循環(huán)神經(jīng)網(wǎng)絡(luò)的變體，適用于挖掘長時(shí)間序列數(shù)據(jù)的內(nèi)在聯(lián)系，并且有效解決了梯度消失的問題[43]。

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

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

2.3 loss函數(shù)和評(píng)價(jià)指標(biāo)

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

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

2.4 模型訓(xùn)練

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

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

3 結(jié)果與分析

3.1 數(shù)據(jù)貢獻(xiàn)度分析

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

表1 地表反射率（7個(gè)波段）和地表晝夜溫度數(shù)據(jù)對(duì)農(nóng)作物估產(chǎn)的影響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 估產(chǎn)模型性能檢測

為驗(yàn)證CNN-GRU估產(chǎn)模型的優(yōu)越性，將其與使用最廣泛的多個(gè)農(nóng)作物產(chǎn)量預(yù)測模型進(jìn)行比較，如深度學(xué)習(xí)框架下的CNN、GRU和傳統(tǒng)機(jī)器學(xué)習(xí)算法中的SVR、DT和RF。其中所有估產(chǎn)模型均使用同一個(gè)數(shù)據(jù)集。但CNN難以處理長時(shí)間序列數(shù)據(jù)，則將時(shí)間序列和波段相疊加作為CNN的輸入，即（32,32,30×9）。而GRU網(wǎng)絡(luò)需要將輸入數(shù)據(jù)展平為時(shí)間序列向量（30,32×32×9）。此外，由于傳統(tǒng)機(jī)器學(xué)習(xí)模型本身的結(jié)構(gòu)限制，在輸入之前需要把輸入數(shù)據(jù)展平為一維向量，但30×32×32×9=276480維度太高，容易導(dǎo)致“維度災(zāi)難”，因此，使用影像的NDVI作為傳統(tǒng)機(jī)器學(xué)習(xí)模型的輸入。各個(gè)模型的關(guān)鍵參數(shù)設(shè)置如下：

1）卷積神經(jīng)網(wǎng)絡(luò)：包含4個(gè)卷積層、4個(gè)池化層、1個(gè)全連接層，批處理大小為32，初始學(xué)習(xí)率為0.000 1。

2）門控循環(huán)單元模型：由128個(gè)神經(jīng)元的GRU和全連接層組成，批處理大小為32，初始學(xué)習(xí)率為0.001。

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

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

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

在本節(jié)中，使用不同的估產(chǎn)模型分別預(yù)測2016—2018年冬小麥產(chǎn)量，為確保比較的公平性，所有估產(chǎn)模型均在相同的軟硬件環(huán)境下訓(xùn)練和測試。分析表2可知：1）所提出的CNN-GRU估產(chǎn)模型年平均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模型能準(zhǔn)確地預(yù)測冬小麥產(chǎn)量。但對(duì)于部分低產(chǎn)量的區(qū)縣（如新疆維吾爾自治區(qū)的伊寧縣和和碩縣）和部分高產(chǎn)量的區(qū)縣（陜西省西安市的灞橋區(qū)和河南省南陽市的新野縣）CNN-GRU模型預(yù)測誤差偏大（如圖3所示）。產(chǎn)生這種誤差的原因可能是冬小麥數(shù)據(jù)集中低產(chǎn)量和高產(chǎn)量的樣本數(shù)占偏低，導(dǎo)致模型對(duì)于部分低產(chǎn)量和高產(chǎn)量的地區(qū)泛化能力較差，預(yù)測的誤差較高；對(duì)于低產(chǎn)量的新疆維吾爾自治區(qū)海拔高、環(huán)境復(fù)雜等不確定性因素導(dǎo)致預(yù)測精度降低；對(duì)于高產(chǎn)量的灞橋區(qū)和新野縣來說，地處平原，均為冬小麥氣候適宜區(qū)，種植面積較大，而模型輸入中沒有表達(dá)區(qū)域異質(zhì)性的參數(shù)，導(dǎo)致低估了部分高產(chǎn)量地區(qū)的冬小麥產(chǎn)量。2）CNN和GRU的性能相當(dāng)，年平均RMSE分別為920、996 kg/hm2，均高于所提出的CNN-GRU估產(chǎn)模型，表明單一的神經(jīng)網(wǎng)絡(luò)在作物估產(chǎn)領(lǐng)域效果不佳。CNN能有效地提取影像中空間-光譜特征，但作物生長隨時(shí)間變化，農(nóng)作物在不同的生育期有不同的形態(tài)特征，時(shí)間信息的缺失嚴(yán)重影響了CNN的整體性能；而GRU能整合長時(shí)間序列的作物生長信息，但不能捕捉復(fù)雜的空間關(guān)系。因此所提出的CNN-GRU估產(chǎn)模型可以有效地提取冬小麥生長期的空間-光譜-時(shí)間信息，并將其應(yīng)用于冬小麥產(chǎn)量預(yù)測中。3）深度學(xué)習(xí)網(wǎng)絡(luò)在性能上均優(yōu)于傳統(tǒng)的機(jī)器學(xué)習(xí)模型，在冬小麥產(chǎn)量預(yù)測這一復(fù)雜問題的研究中，深度學(xué)習(xí)網(wǎng)絡(luò)更能挖掘冬小麥影像中的時(shí)空特征。

表2 不同估產(chǎn)模型性能對(duì)比Table 2 Performance comparison of different yield estimation models

3.3 不同生育期冬小麥預(yù)測

本節(jié)分別對(duì)作物累計(jì)生育期的時(shí)間步（0～5，0～10，0～15，0～20，0～25，0～30）和單一生育期的時(shí)間步（0～5，＞5～10，＞10～15，＞15～20，＞20～25，＞25～30）進(jìn)行單產(chǎn)預(yù)測，以評(píng)估CNN-GRU網(wǎng)絡(luò)模型預(yù)測冬小麥產(chǎn)量的最佳時(shí)間窗。

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

3.4 CNN-GRU模型魯棒性檢驗(yàn)

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

表3 CNN-GRU模型逐年精度評(píng)估Table 3 Annual accuracy evaluation of CNN-GRU model

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

4 結(jié) 論

作物產(chǎn)量預(yù)測在農(nóng)業(yè)管理中具有重要意義。近年來的研究表明，遙感是一種有效的估產(chǎn)方法。本文提出了一種基于多時(shí)相遙感影像的CNN-GRU估產(chǎn)模型，利用CNN從多時(shí)相遙感影像中提取豐富的空間-光譜特征，在此基礎(chǔ)上，GRU自適應(yīng)學(xué)習(xí)冬小麥生育期各階段之間的時(shí)間依賴。試驗(yàn)證明該模型能從多時(shí)相遙感影像中提取聯(lián)合空間-光譜-時(shí)間特征，挖掘作物生長特征與產(chǎn)量之間的復(fù)雜關(guān)系，并能夠?qū)κ斋@前期的冬小麥產(chǎn)量進(jìn)行準(zhǔn)確預(yù)測。結(jié)論如下：

1）深度學(xué)習(xí)技術(shù)在產(chǎn)量預(yù)測上有巨大潛力，對(duì)比結(jié)果表明深度學(xué)習(xí)方法優(yōu)于傳統(tǒng)機(jī)器學(xué)習(xí)方法。同時(shí)，與單一的深度學(xué)習(xí)網(wǎng)絡(luò)CNN和GRU相比，本文所提出的CNN-GRU估產(chǎn)模型精準(zhǔn)預(yù)測了2016－2018年中國冬小麥主產(chǎn)區(qū)的產(chǎn)量，其精度均優(yōu)于其他估產(chǎn)模型。說明該模型能有效地提取反映作物生長信息的空間-光譜-時(shí)間特征，解決了單一神經(jīng)網(wǎng)絡(luò)提取作物生長特征不充分問題，進(jìn)一步提高了大范圍冬小麥單產(chǎn)估算效率。

2）基于混合神經(jīng)網(wǎng)絡(luò)的遙感估產(chǎn)模型魯棒性較高。本文對(duì)2008－2018年數(shù)據(jù)集逐年樣本訓(xùn)練，所提出的CNN-GRU估產(chǎn)模型RMSE和MAE平均值分別為727、534 kg/hm2，R2平均值為0.805。該模型在不同數(shù)據(jù)集上預(yù)測產(chǎn)量與統(tǒng)計(jì)產(chǎn)量值離散程度基本一致，CNN-GRU估產(chǎn)模型性能穩(wěn)定，預(yù)測精度高，可以作為一種高效的方法為全國提供常規(guī)可行的冬小麥產(chǎn)量監(jiān)測和預(yù)報(bào)。

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