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

(1.中國農(nóng)業(yè)大學(xué)土地科學(xué)與技術(shù)學(xué)院， 北京 100083； 2.農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)災(zāi)害遙感重點(diǎn)實(shí)驗(yàn)室， 北京 100083；3.中國農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院， 北京 100083)

0 引言

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

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

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

1.2.2Savitzky-Golay濾波平滑處理

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

圖2 S-G濾波后LAI時(shí)間序列與原始LAI時(shí)間序列對(duì)比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預(yù)測精度也不同[14]。本文主要研究夏玉米的LAI預(yù)測，為消除其他地物類型的影響，對(duì)研究區(qū)域地物類型進(jìn)行了劃分，并在進(jìn)行LAI時(shí)間序列預(yù)測研究時(shí)利用分類信息對(duì)夏玉米種植區(qū)進(jìn)行掩膜處理，而棉花、林地和草地等其他地物類型排除在外。具體農(nóng)作物分類及夏玉米種植區(qū)域提取通過文獻(xiàn)[15]提出的應(yīng)用一階差分法和重構(gòu)LAI的傅里葉變換的諧波特征識(shí)別方法實(shí)現(xiàn)，研究區(qū)域作物識(shí)別的總體精度為82.51%，夏玉米的識(shí)別精度為88.5%。

2 研究方法

2.1 LAI預(yù)測模型

2.1.1ARIMA預(yù)測模型

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

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

(1)

其中

(2)

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

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

B——后移算子

εt——白噪聲序列

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

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

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

2.1.2RBF神經(jīng)網(wǎng)絡(luò)模型

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

圖4 徑向基(RBF)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖Fig.4 RBF neural network structure diagram

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

(3)

式中xs——隱含層第s個(gè)神經(jīng)元的輸入向量

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

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

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

k——隱含層神經(jīng)元個(gè)數(shù)

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

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

O=WTH=w1h1+w2h2+…+wkhk

(4)

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

(5)

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

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

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

(2)求解方差σi

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

(6)

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

(3)計(jì)算隱含層和輸出層間的權(quán)值

隱含層至輸出層之間神經(jīng)元的連接權(quán)值W，可利用最小二乘法直接計(jì)算，即

(7)

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

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

2.2 夏玉米單產(chǎn)預(yù)測方法

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

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

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

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

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

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

3 結(jié)果與分析

3.1 葉面積指數(shù)預(yù)測結(jié)果及精度評(píng)價(jià)

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

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

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

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

3.2 玉米單產(chǎn)預(yù)測結(jié)果及精度評(píng)價(jià)

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

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

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

3.2.1縣域尺度單產(chǎn)預(yù)測結(jié)果精度評(píng)價(jià)

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

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

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

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

3.2.2像素尺度單產(chǎn)預(yù)測結(jié)果精度評(píng)價(jià)

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

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

4 討論

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

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

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

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

5 結(jié)論

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

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