李 俐 許連香 王鵬新 齊 璇 王 蕾

(1.中國農業(yè)大學土地科學與技術學院， 北京 100083； 2.農業(yè)農村部農業(yè)災害遙感重點實驗室， 北京 100083；3.中國農業(yè)大學信息與電氣工程學院， 北京 100083)

0 引言

近年來，隨著經濟全球化的深入，國際國內糧食市場競爭日益激烈，作物產量信息受到社會各界越來越多的關注。玉米是我國主要糧食作物之一[1]，對玉米進行大區(qū)域長勢監(jiān)測和單產預測，能夠幫助農業(yè)經營者和相關部門制定合理的生產計劃，并及時優(yōu)化調整玉米種植結構，同時對保障糧食安全、促進糧食貿易也具有重要意義。

目前，作物產量預測方法主要基于作物生長機理模型和經驗回歸模型[2]。作物生長機理模型通過輸入氣象、土壤、田間管理和作物品種遺傳特性等參數(shù)，可以在單點尺度上動態(tài)模擬作物生長發(fā)育的變化及產量的形成過程，適用于點尺度的作物產量預測[3]；基于數(shù)據(jù)同化方法耦合遙感觀測數(shù)據(jù)和作物生長機理模型，可實現(xiàn)區(qū)域尺度的作物產量預測[4]，但由于過程復雜、輸入?yún)?shù)眾多，在一定程度上限制了其在大面積作物產量預測方面的廣泛應用。經驗回歸模型通常選取與產量相關的特征參數(shù)，與作物單產建立統(tǒng)計相關模型，預測作物產量，是一種簡便且符合業(yè)務化運行要求的大范圍作物產量預測方法[5]。

經驗回歸模型常采用的特征參數(shù)包括歸一化植被指數(shù)( Normalized difference vegetation index，NDVI) 、比值植被指數(shù)(Ratio vegetation index，RVI)、差值植被指數(shù)(Difference vegetation index，DVI)和綠度植被指數(shù)(Greenness vegetation index，GVI) 等[6]，這些植被指數(shù)可間接反映作物的長勢和生物量信息。隨著定量遙感技術的發(fā)展，作物葉面積指數(shù)(Leaf area index，LAI)和生物量等越來越多的陸地生物物理量可通過遙感技術獲得，這類參數(shù)是作物產量形成的物質基礎，可決定產量可能達到的最高上限[7]。有研究表明，葉面積指數(shù)與光合作用、蒸騰作用以及凈初級生產力等密切相關，因此常作為監(jiān)測作物長勢與預測作物產量的特征參數(shù)[8]。目前，基于作物葉面積指數(shù)進行作物產量預測?；趩蝹€生育期或某個時段的葉面積指數(shù)，可實現(xiàn)作物收獲前一段時間的預測單產[9]。然而，作物單產是作物在整個生長階段長勢的綜合反映[10]，不同階段的作物長勢對作物單產的影響不同，僅依靠單個生育期或某個時段LAI反映的長勢信息進行單產預測會影響單產預測精度。王鵬新等[11]綜合考慮作物不同生育時期LAI對產量的影響程度，利用隨機森林回歸法賦予不同生育時期LAI不同權重，得到加權LAI，以反映作物的綜合長勢，進一步提高單產預測的準確性，但獲取作物整個生長階段(出苗期—成熟期)LAI所需時間周期較長。因此，可嘗試利用預測模型獲取夏玉米收獲前生育時期的LAI數(shù)據(jù)，并在此基礎上進行夏玉米的早期單產預測研究。

本文以河北中部平原為研究區(qū)，以夏玉米為研究對象，選取與籽粒產量密切相關的LAI作為遙感特征參數(shù)，利用求和自回歸移動平均模型(Autoregressive integrated moving average model，ARIMA)和徑向基函數(shù)(Radial basis function，RBF)神經網(wǎng)絡對該地區(qū)夏玉米的LAI進行預測，并基于LAI監(jiān)測數(shù)據(jù)和預測數(shù)據(jù)，結合河北中部平原加權LAI與產量的相關性研究成果，對研究區(qū)2016—2018年夏玉米單產進行預測，以期提高夏玉米單產預測的準確性和時效性。

1 材料與方法

1.1 研究區(qū)概況

選擇河北省中部平原的保定市、石家莊市、廊坊市、滄州市和衡水市的部分或全部地區(qū)為研究區(qū)域，其覆蓋范圍為114°32′～117°36′E，36°57′～39°50′N，面積約為5.3×104km2，如圖1所示。該地區(qū)屬半干旱半濕潤大陸性季風氣候，雨熱同期，年均降水量范圍為350～700 mm，且時空分布不均，降水主要集中在夏季，占全年的65%～70%，降水量由南向北逐漸減少。該地區(qū)是我國重要的玉米生產基地[12]，冬小麥-夏玉米輪作是其典型的種植模式。

圖1 研究區(qū)域位置Fig.1 Location of study area

1.2 數(shù)據(jù)獲取與處理

1.2.1LAI數(shù)據(jù)預處理

選取2010—2018年夏玉米生育期(7—9月)，基于Terra和Aqua衛(wèi)星上的MODIS傳感器獲得的MODIS LAI產品MCD15A3H。與MOD15A2和MYD15A2產品相比，MCD15A3H產品具有較高的時空分辨率，其時間分辨率為4 d，空間分辨率為500 m，有利于作物長勢和物候的監(jiān)測。利用MRT對原始LAI影像進行鑲嵌、重采樣、投影轉換和裁剪等預處理，并且輸出投影統(tǒng)一為Lambert投影。

1.2.2Savitzky-Golay濾波平滑處理

為消除由于云、大氣等因素的噪聲影響引起的LAI數(shù)據(jù)驟降現(xiàn)象，本文應用上包絡線Savitzky-Golay(S-G)濾波對預處理后LAI進行平滑處理[13]。首先對LAI時序數(shù)據(jù)的長期變化趨勢進行擬合。再通過局部循環(huán)迭代S-G濾波使得擬合結果更接近于LAI時序數(shù)據(jù)的上包絡曲線。以饒陽縣某像素為例，圖2為其經過S-G濾波(平滑多項式次數(shù)為2，窗口尺寸為5)處理后的MODIS LAI時間序列曲線和原始LAI時間序列曲線。可以看出，原始LAI時間序列波谷到波峰的曲線呈鋸齒狀的不規(guī)則波動，而經S-G濾波處理后的LAI曲線在保持原有曲線基本形狀的基礎上保證了曲線的平穩(wěn)變化，更符合作物生長的物候特征，有利于進行時間序列LAI的變化趨勢分析及作物長勢信息提取。最后逐像素取每旬所包含的多時相LAI最大值作為該旬的LAI值，通過疊加研究區(qū)域夏玉米種植區(qū)圖(圖3)，得到研究區(qū)域夏玉米種植區(qū)2010—2018年每年7—9月以旬為單位的LAI時間序列數(shù)據(jù)。

圖2 S-G濾波后LAI時間序列與原始LAI時間序列對比Fig.2 Comparison of S-G filtered LAI and original LAI

圖3 河北中部平原夏玉米種植區(qū)(2017年)Fig.3 Summer maize planting area in central plain of Hebei Province (2017)

1.3 夏玉米種植區(qū)的獲取

研究表明，地表植被覆蓋類型不同，其LAI預測精度也不同[14]。本文主要研究夏玉米的LAI預測，為消除其他地物類型的影響，對研究區(qū)域地物類型進行了劃分，并在進行LAI時間序列預測研究時利用分類信息對夏玉米種植區(qū)進行掩膜處理，而棉花、林地和草地等其他地物類型排除在外。具體農作物分類及夏玉米種植區(qū)域提取通過文獻[15]提出的應用一階差分法和重構LAI的傅里葉變換的諧波特征識別方法實現(xiàn)，研究區(qū)域作物識別的總體精度為82.51%，夏玉米的識別精度為88.5%。

2 研究方法

2.1 LAI預測模型

2.1.1ARIMA預測模型

ARIMA(p,d,q)模型[16-17]可看作ARMA(p,q)模型的推廣，參數(shù)p為自回歸階數(shù)，q為移動平均階數(shù)，d為差分階數(shù)，可應用于非平穩(wěn)時間序列預測。利用ARIMA模型進行預測的基本步驟[18]為：

(1)數(shù)據(jù)平穩(wěn)化處理：首先利用單位根檢驗(Augmented dickey-fuller，ADF)判斷LAI時間序列的平穩(wěn)性。若為非平穩(wěn)序列，則采用d階差分處理將其轉換為平穩(wěn)序列，對平穩(wěn)后的LAI時間序列擬合ARMA(p,q)模型

(1)

其中

(2)

式中Xt——LAI樣本時間序列數(shù)據(jù)

θi、φi——模型參數(shù)

B——后移算子

εt——白噪聲序列

(2)模型定階：ARMA(p,q)模型的參數(shù)p、q可通過平穩(wěn)時間序列的自相關函數(shù)(Autocorrelation function，ACF)和偏自相關函數(shù)(Partial autocorrelation function，PACF)確定。若序列的自相關函數(shù)拖尾，而偏自相關函數(shù)在p步截尾，則可建立AR(p)模型；若序列的偏自相關函數(shù)拖尾，而自相關函數(shù)在q步截尾，則可建立MA(q)模型；若序列的自相關函數(shù)和偏自相關函數(shù)均拖尾，則可建立ARMA(p,q)模型。然后利用最小信息準則(Akaike information criterion，AIC)進行模型優(yōu)選，逐像素確定自回歸階數(shù)p和移動平均階數(shù)q。AIC準則利用似然函數(shù)估計值最大的原則來確定適用的模型，AIC值最小的模型即為最優(yōu)模型。

(3)參數(shù)估計：模型定階后，對選定模型中的參數(shù)θi、φi進行估計。常用的估計方法有矩估計、最小二乘估計、極大似然估計等。極大似然估計可充分利用序列值的信息，精度較高，故本文利用極大似然估計法對參數(shù)θi、φi進行估計。

(4)模型檢驗：對已建立模型進行顯著性檢驗，若模型殘差序列是白噪聲，認為模型擬合顯著有效；否則說明殘差序列中仍有信息未被提取，需要重新擬合模型。

2.1.2RBF神經網(wǎng)絡模型

徑向基神經網(wǎng)絡[19-20]是一種3層前向神經網(wǎng)絡，網(wǎng)絡結構較為簡單、訓練簡潔并且收斂速度快，能夠逼近任意非線性函數(shù)。其基本思想是用徑向基函數(shù)作為隱含層節(jié)點的“基” 構成隱含層空間，利用隱含層對輸入矢量進行變換，將低維的模式輸入數(shù)據(jù)映射到高維空間，然后通過對隱含層節(jié)點的輸出加權求和得到輸出，其結構如圖4所示。

圖4 徑向基(RBF)神經網(wǎng)絡結構示意圖Fig.4 RBF neural network structure diagram

徑向基神經網(wǎng)絡中常用的徑向基函數(shù)是高斯基函數(shù)，設神經網(wǎng)絡的徑向基向量為H=[h1,h2,…,hk]T，hk為基函數(shù)，其表達式為

(3)

式中xs——隱含層第s個神經元的輸入向量

ci——隱含層第i個神經元高斯基函數(shù)的中心

σi——隱含層第i個神經元高斯基函數(shù)的方差

‖·‖——歐式范數(shù)，表示輸入向量與中心向量的距離

k——隱含層神經元個數(shù)

S——輸入向量總數(shù)

假設RBF神經網(wǎng)絡隱含層和輸出層的連接權值為W=[w1,w2,…,wk]T，則輸出層的輸出為

O=WTH=w1h1+w2h2+…+wkhk

(4)

將式(3)代入式(4)，可得輸出表達式為

(5)

基于自組織選取中心的RBF神經網(wǎng)絡學習算法進行參數(shù)的求解，該方法由兩個階段組成，一是自組織學習階段：求解隱含層基函數(shù)的中心與方差；二是監(jiān)督學習階段：求解隱含層到輸出層權值的階段。具體步驟為：

(1)基于K-均值聚類方法求基函數(shù)中心[21]

初始化聚類中心。隨機選取k個訓練樣本作為聚類中心ci(i=1,2,…,k)。將輸入的訓練樣本集合按最近鄰規(guī)則分組。即按照xs與中心ci之間的歐氏距離將xs分配給中心為ci的輸入樣本集合?s。重新調整聚類中心，計算各個聚類集合中訓練樣本平均值，即新的聚類中心ci，直到新的聚類中心不再發(fā)生變化時，所得到的ci為RBF神經網(wǎng)絡最終的基函數(shù)中心。

(2)求解方差σi

本文RBF神經網(wǎng)絡的基函數(shù)為高斯基函數(shù)，其方差σi計算式為

(6)

式中cmax——所選取中心間的最大距離

(3)計算隱含層和輸出層間的權值

隱含層至輸出層之間神經元的連接權值W，可利用最小二乘法直接計算，即

(7)

RBF神經網(wǎng)絡訓練過程中，徑向基函數(shù)神經網(wǎng)絡可自動增加徑向基神經元，直到達到所設立的誤差為止。另外，徑向基函數(shù)的擴展速度對網(wǎng)絡的收斂性影響很大，一般不宜過大。經過反復試驗和調試，設定均方誤差為0.01，徑向基函數(shù)擴展速度為0.5。

將2010年7月上旬—2018年8月下旬的LAI數(shù)據(jù)作為分析建模數(shù)據(jù)，2016—2018年每年9月上旬—下旬的LAI數(shù)據(jù)作為檢驗數(shù)據(jù)。逐像素提取多旬LAI建模數(shù)據(jù)組成一維時間序列，分別作為模型的輸入數(shù)據(jù)，運用ARIMA模型對河北中部平原所有像素未來的LAI變化狀況進行多步預測[22]。RBF神經網(wǎng)絡利用相同的輸入數(shù)據(jù)，設置輸入層節(jié)點個數(shù)為3，輸出層節(jié)點個數(shù)為1[23]，基于滾動預測[24]的方法對LAI時間序列進行多步預測。在預測過程中，將得到的預測值作為下一步預測的輸入來計算出下一步的預測值，如此迭代，即完成LAI向前1、2、3步的預測。

2.2 夏玉米單產預測方法

根據(jù)該地區(qū)夏玉米主要生育期劃分，一般7月上旬至7月中旬為出苗—拔節(jié)期、7月下旬至8月上旬為拔節(jié)—抽雄期、8月中旬至9月上旬為抽雄—乳熟期、9月中旬至9月下旬為乳熟—成熟期，分別取各生育時期內所包含的多旬LAI的平均值作為該生育時期的LAI值?；谖墨I[25]各生育期LAI的權重系數(shù)：出苗—拔節(jié)期為0.148 3，拔節(jié)—抽雄期為0.360 8，抽雄—乳熟期為0.274 5，乳熟—成熟期為0.216 4。另外，由于文獻[25]中為對比分析不同特征參數(shù)與夏玉米單產的線性關系，對各生育期權重系數(shù)對應的多旬LAI值均按最大值為7，最小值為0進行了歸一化處理，值域范圍為[0,1]。因此，在使用上述權重系數(shù)計算各縣(區(qū))加權LAI時需對旬LAI進行同樣的歸一化處理。

對河北中部平原53個縣(區(qū))2010—2015年夏玉米主要生育期的加權LAI與玉米單產進行回歸分析[25]，得到夏玉米單產的回歸模型為：Y=5 651X+4 493，P<0.001。其中，X表示加權LAI，Y表示夏玉米單產(kg/hm2)。基于2010—2015年的單產監(jiān)測結果與夏玉米實際單產結果，應用線性回歸分析的方法分析兩者之間的相關性(圖5)，可以發(fā)現(xiàn)，在樣本數(shù)n=318情況下，監(jiān)測單產與實際單產的相關系數(shù)(r)為0.487，達到顯著性水平(P<0.001)，決定系數(shù)(R2)為0.237，說明基于粒子群優(yōu)化投影尋蹤法的估產模型的精度較高。

圖5 基于粒子群優(yōu)化投影尋蹤法的監(jiān)測單產與實際單產散點圖Fig.5 Scattered plot of actual yields and estimated ones based on projection pursuit with particle swarm optimization

利用加權LAI與夏玉米單產間的線性關系，結合LAI預測數(shù)據(jù)，可在玉米收獲前1—3旬進行夏玉米單產預測。其中，向前1旬的單產預測過程為：利用7月上旬至9月中旬的LAI監(jiān)測數(shù)據(jù)和9月下旬的LAI預測數(shù)據(jù)，得到加權LAI，再根據(jù)上述產量回歸模型得到夏玉米單產1旬預測結果；向前2旬的單產預測過程為：利用7月上旬至9月上旬的LAI監(jiān)測數(shù)據(jù)和9月中旬、9月下旬的LAI預測數(shù)據(jù)，得到加權LAI，再根據(jù)產量回歸模型得到夏玉米單產2旬預測結果；向前3旬的單產預測過程為：利用7月上旬至8月下旬的LAI監(jiān)測數(shù)據(jù)和9月上旬至9月下旬的預測數(shù)據(jù)，得到加權LAI，再根據(jù)產量回歸模型得到夏玉米單產3旬預測結果。

應用相對誤差(Relative error，RE)、絕對誤差(Absolute error，AE)與均方根誤差(Root mean square error，RMSE)作為河北中部平原夏玉米單產預測精度的評價指標。首先，分別統(tǒng)計各縣(區(qū))預測單產及監(jiān)測單產的平均產量，并計算各縣(區(qū))平均產量的相對誤差，以對縣域尺度上單產預測結果進行精度評價。

逐像素計算預測單產與監(jiān)測單產間的相對誤差和絕對誤差，得到預測單產相對誤差的空間分布圖并統(tǒng)計每旬預測單產結果所有像素絕對誤差的最大值、最小值、平均值和均方根誤差，對像素尺度上單產預測結果進行精度評價。

3 結果與分析

3.1 葉面積指數(shù)預測結果及精度評價

基于ARIMA模型和RBF神經網(wǎng)絡分別逐像素對河北中部平原LAI進行預測，圖6b、6c從上至下依次為基于兩模型2017年8月下旬分別向前預測1、2、3步得到的2017年9月上旬至下旬LAI預測結果。ARIMA模型和RBF神經網(wǎng)絡1步、2步和3步預測結果都基本反映了葉面積指數(shù)監(jiān)測結果的時空分布特征。從時間上看，9月上旬至下旬屬于夏玉米的花粒期階段(包括抽雄—乳熟期和乳熟—成熟期)中后期，此時夏玉米從以營養(yǎng)生長為主轉向以生殖生長為主，生殖器官生長速度加快，而營養(yǎng)器官(葉片等)生長速率減緩，植株下部葉齡較大葉片開始枯黃，LAI開始緩慢下降。兩模型LAI預測結果均準確反映了夏玉米LAI在9月上旬至9月下旬間變化特點，葉面積指數(shù)隨時間的增長而降低，并且受降水等環(huán)境因子的影響，各像素葉面積指數(shù)的下降速率不同，降水較豐沛的地區(qū)LAI下降較慢，而降水偏少的地區(qū)LAI下降較快，主要是因為降水充足可延緩葉片衰老，這與文獻[26]的研究結果一致。以石家莊市藁城區(qū)為例，其2017年9月降水較常年偏多，兩模型LAI預測結果均反映了藁城區(qū)各像素LAI較其他同期降水偏少的縣區(qū)(新樂市、高邑縣等)下降速率慢，表明兩模型對LAI變化趨勢預測較為準確。從空間上看，預測結果顯示河北中部平原葉面積指數(shù)具有較明顯的區(qū)域特征，滄州市西部及衡水市南部地區(qū)LAI偏高，廊坊市及滄州市東部地區(qū)LAI偏低，整體特征與監(jiān)測結果較吻合。

圖6 2017年9月預測結果及實際監(jiān)測結果Fig.6 Forecasting and monitoring results of September 2017

圖7 兩種模型預測結果絕對誤差頻數(shù)分布圖Fig.7 Frequency distributions of absolute errors of forecasting results of two models in September 2017

對比分析兩個模型的LAI預測精度，計算得到兩模型2016—2018年1步、2步和3步預測與監(jiān)測結果的絕對誤差(預測值與監(jiān)測值的差)和絕對誤差頻數(shù)分布圖(共3×22 985個像素，圖7)。結果表明，1步預測結果的絕對誤差分布較為集中，兩模型峰值十分接近，ARIMA模型為0.04 m2/m2，RBF神經網(wǎng)絡模型為0.07 m2/m2。不同的是LAI預測結果的誤差范圍，對比頻數(shù)大于100時絕對誤差的分布范圍，ARIMA模型主要分布在[-1.86 m2/m2，1.43 m2/m2]，而RBF神經網(wǎng)絡主要分布在[-2.56 m2/m2，1.77 m2/m2]，較ARIMA模型誤差分布更為分散。隨預測步長增加，兩模型誤差范圍均呈增大趨勢，3步預測結果誤差分布較2步預測結果分散。另外，逐像素計算得到兩模型2016—2018年1步、2步和3步預測結果的平均絕對誤差和均方根誤差(表1)，結果表明，基于ARIMA模型1、2步預測結果的MAE和RMSE均低于基于RBF神經網(wǎng)絡的誤差，MAE分別降低了0.12、0.05 m2/m2，RMSE分別降低了0.18、0.14 m2/m2，3步預測結果的MAE雖較RBF神經網(wǎng)絡的MAE高0.10 m2/m2，但兩者RMSE相等。整體來看，ARIMA模型預測結果準確性和穩(wěn)定性更好，預測結果反映的LAI變化與實際情況更為吻合，更適合河北中部平原的夏玉米LAI預測。

3.2 玉米單產預測結果及精度評價

利用LAI遙感監(jiān)測結果結合基于ARIMA模型的LAI預測結果，采用粒子群優(yōu)化投影尋蹤法確定的夏玉米單產回歸模型得到2016—2018年的河北中部平原夏玉米單產監(jiān)測結果和1—3旬單產預測結果(圖8)。由圖8單產監(jiān)測結果可得，河北中部平原西部玉米單產最高，南部和北部次之，東部單產最低。2016年，河北中部平原53縣(區(qū))平均單產為6 912 kg/hm2，西部大部分縣(區(qū))高于7 000 kg/hm2，少部分縣(區(qū))接近8 000 kg/hm2；2017年，大部分地區(qū)單產約為6 800 kg/hm2，東部地區(qū)單產較低，約為6 300 kg/hm2；2018年單產略高于2017年，西部大部分縣(區(qū))單產在7 000 kg/hm2左右，東北部地區(qū)單產相對偏低，在6 300 kg/hm2左右。從預測結果來看，向前1旬、2旬和3旬預測單產均與監(jiān)測單產較吻合，預測單產時空變化規(guī)律與監(jiān)測產量一致，整體預測精度較高。

表1 ARIMA模型和RBF神經網(wǎng)絡模型預測誤差的統(tǒng)計分析Tab.1 Statistical results of forecasting errors of ARIMA model and RBF neural network model m2/m2

圖8 河北中部平原夏玉米產量監(jiān)測和預測結果Fig.8 Monitoring and forecasting yields of summer maize in central plain of Hebei

3.2.1縣域尺度單產預測結果精度評價

計算得到河北中部平原各縣(區(qū))夏玉米2016—2018年的監(jiān)測單產與向前1—3旬預測單產之間的相對誤差及其統(tǒng)計直方圖見圖9，對縣域尺度夏玉米單產預測結果精度進行評價。結果表明，向前1旬的單產預測精度較高，2016—2018年各縣(區(qū))向前1旬的相對誤差十分接近，均分布在1%以內，最大值為2016年定興縣的0.79%。向前2旬較向前1旬的相對誤差有所增加，誤差分布相對分散，最大值為2018年獻縣的3.53%。向前3旬的相對誤差分布與向前2旬相似，并且誤差分布不均勻，誤差最小為2017年海興縣的1.21%，而最大為2018年獻縣的3.73%。從年際變化來看，各縣(區(qū))向前1—3旬的單產預測結果在不同年份的相對誤差均十分接近，誤差波動最大的為向前3旬中獻縣，2018年較2017年增大了3.64個百分點，其余相對誤差的波動值均小于該值。總體來說，向前1旬、2旬和3旬的玉米單產預測結果均與監(jiān)測結果十分接近，預測精度雖然會隨預測步長增加而降低，但整體預測精度較高，各縣(區(qū))單產預測結果相對誤差均在4%以內，說明基于ARIMA模型的LAI預測數(shù)據(jù)可以較好地反映玉米生育后期的長勢變化及干物質向籽粒轉移的能力。

為進一步驗證縣域尺度夏玉米單產預測的精度，基于2016—2018年的單產預測結果與夏玉米監(jiān)測結果，應用線性回歸分析的方法分析兩者之間的相關性，見圖10，可以發(fā)現(xiàn)，在樣本數(shù)n=159的情況下，向前1旬、2旬和3旬預測單產與監(jiān)測單產均呈顯著的正相關(P<0.001)，決定系數(shù)(R2)分別為0.998、0.960和0.947，說明基于該方法的縣域尺度預測單產精度較高。

圖9 河北中部平原各縣(區(qū))玉米單產預測的相對誤差Fig.9 Statistical histograms of relative error of forecasting maize yield in each county (district) in central plain of Hebei Province

圖10 河北中部平原預測單產與監(jiān)測單產散點圖Fig.10 Scattered plots of forecasting yields and monitoring ones in central plain of Hebei Province

3.2.2像素尺度單產預測結果精度評價

逐像素計算河北中部平原預測單產與監(jiān)測單產的相對誤差，得到2016—2018年3種預測單產相對誤差的空間分布圖(圖11)，對像素尺度夏玉米單產預測結果精度進行評價。結果表明，2016—2018年向前1旬預測單產相對誤差分布相似，分別有89.8%、92.0%和93.2%像素相對誤差小于1%，表明向前1旬預測單產精度較好；隨預測時間的增加，2016—2018年向前2旬和向前3旬單產預測結果相對誤差均呈增加趨勢，且2016年和2017年向前2旬和3旬單產預測結果相對誤差空間分布相似，2018年向前2旬和3旬單產預測結果相對誤差較2016年和2017年偏大，其向前3旬單產預測結果具有較大不確定性，但仍有90.3%像素相對誤差均小于5%，表明不同年份間夏玉米向前1旬、2旬和3旬單產預測結果精度存在一定差異，但整體來看不同年份間單產預測精度均較高。

逐像素計算河北中部平原2016—2018年預測單產與監(jiān)測單產的絕對誤差，并統(tǒng)計其絕對誤差的最大值、最小值、均值和均方根誤差(表2)。結果表明，2016—2018年1旬預測單產RMSE分布于27.47～32.17 kg/hm2，且各像素絕對誤差均不大于288 kg/hm2，表明向前1旬的單產預測精度較高且十分穩(wěn)定。誤差隨預測時間的增加而增大，3旬預測單產的各項誤差指標較前兩旬都有明顯的增加，個別像素絕對誤差最大達到1 285 kg/hm2,但考慮到統(tǒng)計像素較多(每旬22 985個)，且絕大多數(shù)像素絕對誤差較低，因此可認為整體預測精度較高，表明基于該方法可實現(xiàn)對研究區(qū)域內夏玉米收獲前3旬進行準確的單產預測。

4 討論

針對以往單產預測模型采用單一生育期或單一時段的特征參數(shù)進行單產預測，本研究綜合考慮夏玉米主要生育期的特征參數(shù)，并且選擇精度較高的基于ARIMA模型的LAI預測數(shù)據(jù)代替夏玉米生育后期的LAI監(jiān)測數(shù)據(jù)，實現(xiàn)了夏玉米收獲前1—3旬進行單產預測，整體預測精度較高，各縣(區(qū))單產預測結果相對誤差均在4%以內。LAI預測精度的關鍵是預測模型的選擇，在今后的研究中可嘗試其他預測模型，如遺傳算法優(yōu)化的RBF神經網(wǎng)絡、ARIMA-RBF組合預測模型等，以期進一步提高LAI模型預測精度。除此之外，夏玉米研究區(qū)識別的精度也會在一定程度上對LAI預測精度產生影響，本文采用的夏玉米種植區(qū)提取方法在一定程度上保證了算法精度，在今后的研究中可通過提高夏玉米識別精度進一步減少混合地物對LAI預測精度的影響。

圖11 河北中部平原夏玉米單產預測相對誤差空間分布Fig.11 Relative error spatial distributions of forecasting summer maize yield in central plain of Hebei Province

表2 河北中部平原夏玉米單產絕對誤差分布統(tǒng)計Tab.2 Absolute errors of monitoring and forecasting yields of summer maize in central plain of Hebei kg/hm2

夏玉米單產的預測精度除了受LAI預測精度影響外，還受各生育期LAI的權重以及加權LAI與夏玉米單產間的回歸方程影響，不同的賦權方法可能導致夏玉米單產預測精度也不同，合理的賦權方法對保證基于該方法的夏玉米單產預測精度十分重要。另外，影響夏玉米單產的因素有很多，除了受到作物自身的生理因素影響外，還受到作物生長的生態(tài)環(huán)境條件的綜合影響，如水分、養(yǎng)分、溫度、光照等，它們可能成為最終產量形成的脅迫因子，使得產量發(fā)生增減的波動。今后的研究中，需要綜合考慮多種因素對作物最終單產的影響，提高作物單產預測模型的普適性和準確性。

5 結論

(1)基于ARIMA模型的LAI 1、2步預測結果的MAE和RMSE均低于基于RBF神經網(wǎng)絡的誤差，MAE分別降低了0.12、0.05 m2/m2，RMSE分別降低了0.18、0.14 m2/m2，3步預測結果的MAE雖較RBF神經網(wǎng)絡高0.10 m2/m2，但兩者RMSE相等。整體來看，ARIMA模型預測結果的準確性和穩(wěn)定性更好，更適合于河北中部平原的夏玉米LAI預測。

(2)基于粒子群優(yōu)化投影尋蹤法確定的單產回歸模型及ARIMA模型，對2016—2018年河北中部平原進行向前1旬、2旬和3旬夏玉米的單產預測，結果表明，無論是縣域尺度還是像素尺度，各旬單產預測結果均與監(jiān)測結果十分接近，并且隨預測時間的增加預測結果的不確定性增大，但整體預測精度較高，2016—2018年縣域尺度預測單產與監(jiān)測單產間最大相對誤差僅為3.73%。