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

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

0 引 言

冬小麥?zhǔn)侵袊蠹Z食作物之一，在農(nóng)業(yè)生產(chǎn)中占有較大比例，同時在中國的糧食安全、糧食生產(chǎn)以及農(nóng)業(yè)可持續(xù)發(fā)展中也具有至關(guān)重要的作用。遙感數(shù)據(jù)中不同波段的光譜信息以及基于此反演的植被指數(shù)可以直接反映作物生長狀況，是對作物在多種影響因素下綜合作用的體現(xiàn)。隨著遙感技術(shù)的發(fā)展，其時空分辨率逐漸提升，大范圍、高時效的地表觀測變?yōu)榭赡堋Ｄ壳袄眠b感數(shù)據(jù)進行冬小麥產(chǎn)量預(yù)測主要有基于作物敏感波段反射率的統(tǒng)計回歸方法[1-3]和基于遙感信息和作物生長模型的數(shù)據(jù)同化方法[4-6]。統(tǒng)計回歸方法通過建立地表反射率、植被指數(shù)等遙感信息與作物產(chǎn)量之間的關(guān)系模型來實現(xiàn)產(chǎn)量估測，操作簡單但存在容易過擬合的缺點[7]。數(shù)據(jù)同化技術(shù)的思想最早由Charney 等[8]提出并在作物估產(chǎn)領(lǐng)域得到應(yīng)用。遙感觀測和作物生長機理模型的數(shù)據(jù)同化方法能顯著提高模型參數(shù)估計和模型動態(tài)模擬能力，改善監(jiān)測精度，但作物生長模型需要大量參數(shù)作為輸入，例如物候特征、土壤特征、品種系數(shù)、天氣數(shù)據(jù)等，而同化狀態(tài)變量偏少嚴(yán)重影響了同化模型的泛化性能，在實際應(yīng)用中很難快速得到大范圍產(chǎn)量預(yù)測結(jié)果且預(yù)測精度還有待提高[9]。在樣本構(gòu)建方面，以往的研究中多使用區(qū)域平均值來構(gòu)建樣本，例如Jiang 等[10]使用多時相的區(qū)域地表反射率進行平均值空間聚合后建立縣級冬小麥單產(chǎn)的預(yù)測模型，但平均值反映的是樣本平均狀態(tài)，代表性較弱；Wolanin 等[11]采用區(qū)域加權(quán)平均值進行空間聚合，雖較平均值提升了樣本的代表性，但無法反映樣本總體的分布情況。區(qū)域樣本特征的構(gòu)建仍是當(dāng)前提升遙感數(shù)據(jù)大區(qū)域應(yīng)用普適性的至關(guān)重要的部分。頻率直方圖是數(shù)值數(shù)據(jù)分布的精確圖形表示，是一個連續(xù)變量概率分布的估計，能夠反映樣本總體分布情況，通過將不同波段的光譜信息和反演得到的植被指數(shù)映射到像素計數(shù)的頻率直方圖中信息損失較少，大大低于平均值，從而提取更多有效信息，進一步提高作物單產(chǎn)預(yù)測的效率。

隨著機器學(xué)習(xí)的不斷發(fā)展，研究人員開始將機器學(xué)習(xí)技術(shù)應(yīng)用到作物遙感估產(chǎn)研究中，并取得較好效果[12]。Ribeiro 等[13-14]利用植被條件指數(shù)、溫度條件指數(shù)和植被健康指數(shù)等遙感指標(biāo)以及葉面積指數(shù)、生物量等作物參數(shù)，建立人工神經(jīng)網(wǎng)絡(luò)模型來預(yù)測作物產(chǎn)量，但訓(xùn)練過程收斂速度慢且易陷入局部極值，模型參數(shù)復(fù)雜[15]。Wen 等[16]利用支持向量機（Support Vector Machine，SVM）算法對作物產(chǎn)量進行了準(zhǔn)確預(yù)測，但核函數(shù)的選取和參數(shù)的確定通?；诮?jīng)驗，對精度影響較大[17]。隨機森林（Random Forest，RF）是由Breiman于2001 年提出的一種基于分類與回歸樹的算法[18]，它在處理多維特征時對多重共線性并不敏感[19]，通過算法的隨機性可以有效避免過擬合問題，在作物產(chǎn)量預(yù)測方面相對于決策樹、支持向量機、人工神經(jīng)網(wǎng)絡(luò)等算法有顯著優(yōu)勢，其在同等運算率下具有較高的預(yù)測精度和泛化性能[20-22]。

本研究以河南省為研究區(qū)，以中分辨率成像光譜儀（Moderate Resolution Imaging Spectroradiometer，MODIS）地表反射率數(shù)據(jù)作為輸入變量，建立基于隨機森林算法的冬小麥遙感估產(chǎn)模型。為了降低特征提取過程的信息損失，引入頻率直方圖法構(gòu)建樣本特征集，提高不同波段地表反射率的預(yù)測精度，以期為縣域尺度冬小麥產(chǎn)量的預(yù)測提供有效方法。

1 材料與方法

1.1 研究區(qū)概況

河南省位于中國中東部（31°23′～36°22′N，110°22′～116°38′E）（圖1），是中國重要的糧食生產(chǎn)地區(qū)之一，屬北亞熱帶濕潤氣候與暖溫帶半濕潤季風(fēng)氣候間的過渡氣候，年均降水量為500～900 mm，降水季節(jié)分布不均，主要集中在夏季。該地區(qū)冬小麥多分布在中東部地區(qū)，西部山區(qū)種植面積較少，9 月中下旬至10 月上旬播種，12 月中旬后越冬，翌年2 月下旬至3 月上旬開始返青，3月下旬至4 月中旬拔節(jié)，4 月下旬至5 月中旬抽穗，5 月底至6 月初成熟。河南省冬小麥種植區(qū)分布圖和河南省高程、積溫、降水?dāng)?shù)據(jù)如圖1 所示。

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

本研究采用的遙感數(shù)據(jù)為MODIS 8 d 合成地表反射率產(chǎn)品（MOD09A1）和土地覆蓋類型產(chǎn)品（MCD12Q1）。分別提取河南省2003—2015 年間每年3—5 月MOD09A1產(chǎn)品的紅光波段（620～670 nm）、近紅外1 波段（841～876 nm）、藍光波段（459～479 nm）、綠光波段（545～565 nm）、近紅外2 波段（1 230～1 250 nm）、短波紅外1 波段（1 628～1 652 nm）、短波紅外2 波段（2 105～2 155 nm）共7 個波段信息，應(yīng)用MODIS 重投影工具（Modis Reprojection Tool，MRT）對影像的研究區(qū)域進行裁剪，投影轉(zhuǎn)換為WGS84 坐標(biāo)系統(tǒng)，并利用MCD12Q1植物功能型分類方案（Land Cover Type5）識別冬小麥種植區(qū)。產(chǎn)量數(shù)據(jù)來自《河南省農(nóng)村統(tǒng)計年鑒》[23]中縣級冬小麥單產(chǎn)數(shù)據(jù)，2003—2015 年間共有106 個連續(xù)種植冬小麥的縣市（圖1），剔除單產(chǎn)缺失年份，共獲得1 369 個有效單產(chǎn)數(shù)據(jù)作為本研究的單產(chǎn)樣本。氣象數(shù)據(jù)來自國家氣象信息中心數(shù)據(jù)共享服務(wù)平臺《中國地面氣候資料日值數(shù) 據(jù) 集 V3.0 》 （ http://data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0.html）。

1.3 研究方法

1.3.1 特征集構(gòu)造

冬小麥在返青期至灌漿期營養(yǎng)器官和生殖器官迅速生長發(fā)育，新生分蘗和根、莖、葉都將明顯生長，是決定結(jié)實率和千粒質(zhì)量最為關(guān)鍵的時期[23]。歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）是反映作物長勢的重要植被指數(shù)，歸一化水分指數(shù)（Normalized Difference Water Index，NDWI）則在識別植被冠層水分方面有良好表現(xiàn)[24]。將該時段內(nèi)的遙感特征作為產(chǎn)量影響因子，劃分為10 個時間段，分別為3 月6日—3 月13 日、3 月14 日—3 月21 日、3 月22 日—3 月29 日、3 月30 日—4 月6 日、4 月7 日—4 月14 日、4月15 日—4 月22 日、4 月23 日—4 月30 日、5 月1 日—5 月8 日、5 月9 日—5 月16 日、5 月17 日—5 月24日，分別提取各縣域10 個時段內(nèi)的8 d 合成產(chǎn)品MOD09A1 7 個波段信息，并根據(jù)式（1）和式（2）計算得到NDVI 和NDWI。

式中ρRed為紅光波段反射率；ρNIR1為近紅外1 波段反射率；ρNIR2為近紅外2 波段反射率。

為了較為準(zhǔn)確地獲得與縣域冬小麥產(chǎn)量的對應(yīng)關(guān)系，使用多時段的區(qū)域地表反射率，分別采用均值法和頻率直方圖法進行空間聚合構(gòu)造樣本特征集，然后建立縣級冬小麥單產(chǎn)的預(yù)測模型。均值法使用各縣域10 個時間段7 個波段地表反射率的均值作為特征值，共計70 個特征變量；頻率直方圖法使用各縣域10 個時間段7 個波段地表反射率以及NDVI 和NDWI 的頻率直方圖（圖2）來構(gòu)建樣本特征集，由圖2 可知，各波段的地表反射率范圍在[0, 0.6）之間，且峰值出現(xiàn)的位置差異明顯，NDVI的范圍在[0, 1）之間，NDWI 的范圍在（–0.25，0.25）之間。分別將7 個波段的地表反射率以0.01 為間隔等間隔劃分60 等份，NDVI 和NDWI 則以0.1 為間隔等間隔劃分為10 等份，以直方圖各區(qū)間的頻率作為特征值，地表反射率共計4 200 個特征變量，NDVI 和NDWI 共計200個特征變量。

1.3.2 隨機森林算法

隨機森林是由多棵決策樹構(gòu)成的組合分類模型，利用自助法有放回地重復(fù)隨機地從原始樣本集中抽取k個樣本生成訓(xùn)練樣本集合，如此重復(fù)得到原始樣本集的L個子集用來構(gòu)建L棵決策樹，形成隨機森林[18]，模型的最終輸出由森林中的每一棵決策樹共同決定。本研究采用隨機森林回歸模型，以研究區(qū)逐年各縣的M個特征變量及其對應(yīng)的冬小麥單產(chǎn)共同構(gòu)成隨機森林的樣本集，在建模時，隨機從M個特征變量中選取m個變量（m

1.4 模型精度評價與驗證

采用決定系數(shù)（Coefficient of determination，R2）、均方根誤差（Root Mean Square Error，RMSE）和平均絕對誤差（Mean Absolute Error，MAE）3 個指標(biāo)來評價模型精度。其中，R2反映了冬小麥實測產(chǎn)量與估測產(chǎn)量的擬合優(yōu)度，取值范圍為（0,1），數(shù)值越大擬合效果越好。RMSE 和MAE 則反映冬小麥產(chǎn)量估測的精密程度。

2 結(jié)果與分析

2.1 均值法和頻率直方圖法結(jié)果對比

以2003－2012 年的7 個波段地表反射率數(shù)據(jù)和冬小麥單產(chǎn)作為訓(xùn)練樣本，分別采用均值法和頻率直方圖法構(gòu)建隨機森林冬小麥單產(chǎn)預(yù)測模型，以2013—2015 年數(shù)據(jù)作為驗證樣本，冬小麥產(chǎn)量估測結(jié)果如圖3 所示。由圖3 可知，均值法和頻率直方圖法的MAE 和RMSE 分別為660、860 kg/hm2，R2分別為0.76、0.83，均值法與頻率直方圖法均通過顯著性檢驗且達極顯著水平（P<0.001）。頻率直方圖法的擬合線更貼近于1∶1 線，在冬小麥單產(chǎn)實測值偏高的地方預(yù)測準(zhǔn)確度更高。相較于均值法，頻率直方圖法能達到更好的預(yù)測效果，從原理上看，均值法的樣本特征體現(xiàn)的是變量在不同時間段的平均狀態(tài)，而頻率直方圖法的樣本特征體現(xiàn)了不同時間段的變量概率分布，相當(dāng)于將均值展開，從而擁有更多的細(xì)節(jié)特征，在隨機森林學(xué)習(xí)過程中可以獲得變量在不同區(qū)間的影響。

2.2 不同波段地表反射率組合結(jié)果

分別以2003－2012 年的7 個不同波段的地表反射率采用直方圖頻率法構(gòu)建特征集，利用2013－2015 年產(chǎn)量數(shù)據(jù)進行驗證（表1）。其中，近紅外1 波段精度最高，R2達到0.76，RMSE 和MAE 分別為784 和636 kg/hm2，優(yōu)于其他波段；其次是近紅外2 波段和紅波段，R2分別為0.67 和0.63，藍波段最低，R2為0.58。植被的反射光譜曲線在近紅外波段呈迅速上升態(tài)勢，呈陡而近于直線的形態(tài)，在冬小麥產(chǎn)量估測和長勢分析方面有顯著優(yōu)勢[26]，紅波段和近紅外波段則分別是葉綠素以及水和二氧化碳的強吸收帶，植被的反射光譜曲線在600～700、1 000～1 300 nm 處具有波谷的形態(tài)且光譜反射率數(shù)值較低，可看成是反映冬小麥光合作用和冠層水分條件的重要波段。

表1 不同波段對冬小麥產(chǎn)量估測結(jié)果對比Table 1 Comparison of yield estimation results of winter wheat by surface reflectance of different bands

由不同波段地表反射率對冬小麥產(chǎn)量估測結(jié)果對比可知，紅光波段、近紅外1 波段、近紅外2 波段與冬小麥單產(chǎn)之間的關(guān)系較為密切。將其進行多波段組合，分別有紅光波段+近紅外1 波段組合、近紅外1 波段+近紅外2 波段組合、紅光波段+近紅外2 波段組合、紅光波段+近紅外1 波段+近紅外2 波段組合4 種組合形式，另外根據(jù)式（1）和式（2）基于紅光波段、近紅外1 波段、近紅外2 波段可計算得到NDVI 和NDWI，將兩者結(jié)合得到NDVI+NDWI 組合。以2003—2012 年的各組合數(shù)據(jù)基于頻率直方圖法構(gòu)建特征集，利用2013—2015 年產(chǎn)量數(shù)據(jù)進行驗證（表2）。由表2 可知，多個波段地表反射率組合的預(yù)測效果要顯著優(yōu)于單一波段地表反射率的效果，R2均高于0.78（P<0.001）。相比單一波段地表反射率數(shù)據(jù)，多個波段組合可獲得更多波段信息，雖然更多波段的加入可能會引入更多的噪聲，但當(dāng)信息強于噪聲時組合的貢獻會得到體現(xiàn)。NDVI、NDWI 指數(shù)的R2分別達到0.73 和0.82，而NDVI+NDWI 組合的效果最好，R2達到0.89，其MAE 和RMSE 也最小，分別為444、527 kg/hm2，通過顯著性檢驗且達極顯著水平（P<0.001）。

表2 不同波段組合對冬小麥產(chǎn)量估測結(jié)果對比Table 2 Comparison of winter wheat yield estimation results of different band surface reflectance combinations

NDVI 主要反映冬小麥長勢和冬小麥營養(yǎng)生長狀況，冬小麥返青至抽穗階段，NDVI 逐漸增加，抽穗期達到峰值，而后抽穗至成熟階段，因為葉片逐漸變黃干枯，呈下降趨勢，故NDVI 與冬小麥產(chǎn)量的相關(guān)性主要體現(xiàn)在抽穗之前。但NDVI 存在兩點主要劣勢，會影響其對產(chǎn)量的預(yù)測效果，1）在冬小麥葉片封壟后容易出現(xiàn)飽和現(xiàn)象，不能充分反映作物的長勢特征，2）在反映災(zāi)害影響方面存在明顯的滯后現(xiàn)象，尤其不能反映短時突發(fā)的作物脅迫狀況，如突然的霜凍低溫等。NDWI 主要反映的是冬小麥冠層水分狀況。在冬小麥拔節(jié)期至抽穗期，水分是限制冬小麥葉片生長的重要因素，影響冬小麥的光合作用，冬小麥子房在發(fā)育過程中需要光合作用的產(chǎn)物使其生長，若沒有足夠的光合產(chǎn)物，子房發(fā)育將中止進而造成籽粒產(chǎn)量下降，從而對冬小麥產(chǎn)量影響較大[27]。冬小麥進入乳熟期后，水分作為溶媒把莖葉中的營養(yǎng)物質(zhì)運輸?shù)阶蚜Ｖ?，水分會直接影響營養(yǎng)物質(zhì)的輸送，進而影響冬小麥產(chǎn)量。在河南省冬麥區(qū)，干旱和晚霜凍害是影響冬小麥產(chǎn)量的兩種主要農(nóng)業(yè)氣象災(zāi)害，均與水分有密切關(guān)系，干旱是冬小麥持續(xù)受水分脅迫導(dǎo)致的災(zāi)害，晚霜凍害則是由于拔節(jié)期至抽穗期氣溫低于0 引起植株體內(nèi)結(jié)冰導(dǎo)致的低溫脅迫災(zāi)害，但充足的水分條件可以有效減弱晚霜凍害造成的影響。NDWI 能快速有效地提取冬小麥冠層的水分含量，因此在水分和低溫脅迫狀態(tài)下的產(chǎn)量預(yù)測方面具有明顯優(yōu)勢。故NDVI 和NDWI 的結(jié)合可以實現(xiàn)優(yōu)勢互補，獲得更優(yōu)預(yù)測效果?；贜DVI 和NDWI 組合的冬小麥實測產(chǎn)量和估測產(chǎn)量間的散點圖如圖4 所示。

2.3 不同時段歸一化植被指數(shù)和歸一化水分指數(shù)組合預(yù)測結(jié)果

基于上述結(jié)果和分析，選取NDVI+NDWI 組合采用頻率直方圖法構(gòu)建特征集，探討在不同時段數(shù)據(jù)對產(chǎn)量預(yù)測效果的影響。

根據(jù)冬小麥關(guān)鍵生育期，將歸一化植被指數(shù)和歸一化水分指數(shù)組合預(yù)測精度時序變化曲線大致分為3 個時段（圖5），分別為3 月6 日—4 月6 日、4 月7 日—4月30 日、5 月1 日—5 月23 日。

由圖5 可知，4 月7 日—4 月30 日的預(yù)測精度高于其他2 個時段，MAE 和RMSE 分別為602 和720 kg/hm2，R2達到0.82。在3 月6 日—4 月6 日時段內(nèi)冬小麥主要進行營養(yǎng)生長，產(chǎn)量形成器官的干物質(zhì)積累過程在這一時期影像中并不能得到充分反映，模型預(yù)測精度較低，MAE和RMSE 分別為664 和808 kg/hm2，R2為0.78。而4 月7 日—4 月30 日大致處于冬小麥抽穗期，此階段冬小麥將光合作用產(chǎn)生的有機物從營養(yǎng)器官轉(zhuǎn)移到籽粒中，植被指數(shù)與最終冬小麥千粒質(zhì)量密切相關(guān)[28]。5 月1 日—5月23 日大致對應(yīng)冬小麥灌漿和成熟階段，冠層和莖稈的營養(yǎng)物質(zhì)向籽粒轉(zhuǎn)移，葉片中的葉綠素含量下降，導(dǎo)致NDVI 下降，與產(chǎn)量相關(guān)性變?nèi)酰穗A段水分對籽粒飽滿結(jié)實有重要影響，所以此階段的相關(guān)性可能主要是來自NDWI 的貢獻，但有待進一步研究。

2.4 模型誤差空間分布

對NDVI 和NDWI 組合模型精度進行驗證，2013—2015 年的精度驗證結(jié)果如圖6 所示，包含實測冬小麥產(chǎn)量圖、預(yù)測冬小麥產(chǎn)量圖以及冬小麥產(chǎn)量相對誤差分布圖。由圖6 可知，實測冬小麥產(chǎn)量和預(yù)測冬小麥產(chǎn)量呈現(xiàn)基本一致的空間分布現(xiàn)象，相對誤差在±20%范圍內(nèi)，西低東高，西部地區(qū)冬小麥實測產(chǎn)量各年平均保持在3 000～4 000 kg/hm2，低于河南省平均產(chǎn)量30%～40%左右，模型預(yù)測單產(chǎn)偏高，相對誤差平均分布在–10%～20%之間；東部平原地區(qū)的冬小麥實測產(chǎn)量各年均保持在7 000～8 000 kg/hm2范圍內(nèi)，預(yù)測產(chǎn)量與其較為一致，相對誤差平均在–5%左右。總體上，相對誤差由西北向東南逐漸遞減，這主要是地形、溫度、降水綜合作用的影響。河南省地勢整體上西高東低，東部是大片的平原區(qū)，西部是山區(qū)，受地形影響，冬小麥產(chǎn)量在地勢較高的地方基本分布在3 000～4 000 kg/hm2，東部平原區(qū)大致處于6 000～8 000 kg/hm2范圍內(nèi)，由西向東呈過渡增加的趨勢。溫度也是影響冬小麥生長發(fā)育的重要條件，冬小麥在關(guān)鍵生育期需要大于0 的積溫2 200～3 800 ℃，由河南省年均積溫分布可知，西北部山區(qū)關(guān)鍵生育期大于0的積溫小于3 000 ℃，不能充分滿足冬小麥正常生長所需，對冬小麥產(chǎn)量影響較大。河南省累計降水量呈顯著的緯向分布，由南向北逐漸減少，但因具有良好的灌溉條件，南北降水的差異對冬小麥產(chǎn)量從空間分布上看未形成顯著制約。南部信陽地區(qū)，熱量、水分條件均較為充足，但因該地區(qū)實行稻麥輪作制，前茬水稻使農(nóng)田長期浸水而影響土壤耕性，導(dǎo)致后茬冬小麥根系吸水吸肥性能降低，故造成低產(chǎn)。綜上，在低產(chǎn)地區(qū)冬小麥預(yù)測產(chǎn)量普遍高于實測產(chǎn)量，而在中高產(chǎn)地區(qū)則較實測產(chǎn)量偏低，造成此誤差的主要原因是隨機森林算法預(yù)測冬小麥產(chǎn)量由多棵決策樹結(jié)果取平均值得到，估產(chǎn)模型本身傾向于樣本的平均狀態(tài)。

3 結(jié) 論

1）采用頻率直方圖法構(gòu)建樣本特征結(jié)合植被指數(shù)建立冬小麥遙感估產(chǎn)模型的預(yù)測效果顯著優(yōu)于均值法（P<0.001）。

2）在單個地表反射率波段中近紅外1 波段精度最高，多波段或多指數(shù)特征組合預(yù)測可以進一步提高預(yù)測精度。歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）和歸一化水分指數(shù)（Normalized Difference Water Index，NDWI）的組合效果最優(yōu)，決定系數(shù)達到0.89，平均絕對誤差和均方根誤差分別為 444、527 kg/hm2，該組合在4 月7 日—4 月30 日相對于其他時間段對產(chǎn)量預(yù)測的影響更大。

3）冬小麥預(yù)測產(chǎn)量和實測產(chǎn)量的相對誤差均分布在±20%范圍內(nèi)，呈現(xiàn)西高東低的分布趨勢，模型預(yù)測結(jié)果在低產(chǎn)地區(qū)普遍高于實測單產(chǎn)，而在中高產(chǎn)地區(qū)較實測產(chǎn)量偏低。

本研究構(gòu)建的基于隨機森林算法的冬小麥遙感估產(chǎn)模型采用頻率直方圖和植被指數(shù)綜合結(jié)合構(gòu)建特征變量集，但頻率直方圖的計算需要以足夠多的像元為基礎(chǔ)，故以縣域作為基本空間單元來計算輸入變量的頻率直方圖，最終結(jié)果也只能反映到縣域，后續(xù)研究會進一步考慮將冬小麥產(chǎn)量反映到像元尺度上。