賈夢琦，蔡振江，胡 建，陳浩松，王 楠，劉紅星

（1. 河北農(nóng)業(yè)大學(xué) 機(jī)電工程學(xué)院，河北 保定 071000； 2. 河北農(nóng)業(yè)大學(xué) 經(jīng)濟(jì)管理學(xué)院，河北 保定 071000；3. 中央農(nóng)業(yè)廣播學(xué)校保定分校，河北 保定 071000）

河北省是我國糧食主產(chǎn)區(qū)之一［1］，保定市糧食年產(chǎn)量基本占到河北省的1/5。近年來，保障糧食安全在促進(jìn)國家和社會經(jīng)濟(jì)發(fā)展方面發(fā)揮著重要作用，糧食產(chǎn)量的預(yù)測成為了1 個重要的研究課題。由于糧食產(chǎn)量通常受到氣象、土地、人類利用活動、制度等因素影響［2-6］，對于糧食產(chǎn)量的準(zhǔn)確預(yù)測十分困難。因此，利用有針對性的模型和方法預(yù)測糧食產(chǎn)量的發(fā)展趨勢具有很強(qiáng)的現(xiàn)實意義。

近年來，眾多學(xué)者為了準(zhǔn)確掌握和預(yù)測全國各地的糧食產(chǎn)量進(jìn)行了大量的研究。王曉玲，孔思琪［7］利用灰色關(guān)聯(lián)度分析黑龍江省糧食產(chǎn)量的影響因素。趙金霞［8］針對糧食產(chǎn)量影響因素利用層次分析法，得出了糧食增量受社會經(jīng)濟(jì)因素影響較大。戎陸慶等［9］采用灰色關(guān)聯(lián)(GRA)方法和BP 神經(jīng)網(wǎng)絡(luò)模型對廣西糧食產(chǎn)量進(jìn)行了預(yù)測，將預(yù)測的相對誤差降低到1.5%以下。王春雨［10］利用遙感技術(shù)分析了孟印緬3 個地區(qū)的氣候因素和人類活動因素對糧食產(chǎn)量的影響，提出了相應(yīng)的管理政策。武單［11］依據(jù)計量地理學(xué)理論，針對湖南省糧食產(chǎn)量的多個影響因素進(jìn)行研究，分析得出湖南省糧食產(chǎn)量的主要影響因素是農(nóng)業(yè)機(jī)械總動力和化肥施用量的結(jié)論。闕斐艷［12］分析湖南省糧食產(chǎn)量變化與農(nóng)村勞動力投入的相關(guān)性發(fā)現(xiàn)，湖南省農(nóng)村勞動力轉(zhuǎn)移對糧食產(chǎn)量沒有顯著影響。李禮連等［13］對1978—2009 年湖南省糧食產(chǎn)量波動進(jìn)行分析認(rèn)為，糧食產(chǎn)量的增加與糧食作物播種面積、化肥使用量和抗災(zāi)防災(zāi)財政投入相關(guān)性較大。

筆者認(rèn)為，現(xiàn)有的預(yù)測方法對于數(shù)據(jù)量的依賴較大，在實際使用中不太方便。因此，需要找到能夠準(zhǔn)確預(yù)測小數(shù)據(jù)集的方法對糧食產(chǎn)量進(jìn)行合理的預(yù)測。本研究采用ARIMA 模型、LSTM 模型以及ARIMAGRNN 組合模型對糧食產(chǎn)量進(jìn)行預(yù)測，探究其預(yù)測效果，可為糧食產(chǎn)量的監(jiān)測和控制提供依據(jù)。

1 研究數(shù)據(jù)與方法

1.1 數(shù)據(jù)來源

本文采用的樣本數(shù)據(jù)為1996—2017 年保定市糧食產(chǎn)量相關(guān)數(shù)據(jù)，樣本數(shù)據(jù)來源于《保定市經(jīng)濟(jì)統(tǒng)計年鑒》。收集1996—2017 年保定市糧食產(chǎn)量及影響因子用于模型構(gòu)建及檢驗。

1.2 統(tǒng)計學(xué)分析

采用Excel 工作表建立保定市糧食產(chǎn)量及糧食產(chǎn)量影響因子數(shù)據(jù)庫，運用Python 軟件進(jìn)行ARIMA 模型、LSTM 模型以及ARIMA-GRNN 組合模型的構(gòu)建和預(yù)測。

1.3 ARIMA 模型的構(gòu)建

ARIMA 模型是時間序列分析模型中較為常用的1 種，主要用于發(fā)現(xiàn)和探索數(shù)據(jù)隨時間變化的規(guī)律［14］。1970 年，Box 和Jenkins 基于時間序列分析發(fā)表了差分自回歸滑動平均模型(Autoregressive Integrated Moving Average model)，簡稱ARIMA 模型［15-16］。

ARIMA(p,d,q)模型中對于p、d、q的確定最為關(guān)鍵，AR 為自回歸模型，p為自回歸平均階數(shù)；p階自回歸過程的定義為：

yt為當(dāng)前值，μ為常數(shù)項，p為自回歸模型的階數(shù)，γi為自相關(guān)系數(shù)，εt為誤差。

MA 為滑動平均模型，q為滑動平均階數(shù)［17］；q階自回歸過程的定義為：

yt為當(dāng)前值，μ為常數(shù)項，q為滑動平均模型的階數(shù)，θi為自相關(guān)系數(shù)，εt-i為誤差。

d為將非平穩(wěn)時間序列變?yōu)槠椒€(wěn)序列的差分運算次數(shù)。

ARIMA 建模過程具體分為5 個步驟：特征分析、模型定階、模型參數(shù)優(yōu)化、模型驗證和模型預(yù)測，如圖1 所示。

（1）特征分析。ARIMA 模型最關(guān)鍵的是需要將非平穩(wěn)時間序列變?yōu)槠椒€(wěn)時間序列。平穩(wěn)時間序列需要均值和方差不發(fā)生明顯的變化，方差越大，數(shù)據(jù)波動越明顯。方差計算公式如下：

σ2為總體方差，X為變量，μ為總體均值，N為總體個數(shù)。方差次數(shù)即ARIMA（p,d,q）中d的階數(shù)。

（2）模型識別。ARIMA(p,d,q)模型中p、q的選擇，根據(jù)自相關(guān)函數(shù)ACF (Auto correlation function)和偏自相關(guān)函數(shù)PACF(Partial autocorrelation function)圖來確定，其中p和q的取值不唯一。ACF 自相關(guān)函數(shù)反映了數(shù)據(jù)在不同時間序列之間的相關(guān)性，自相關(guān)函數(shù)ACF 公式如下：

其中，ρk為自相關(guān)函數(shù)，t為滯后時間，k為滯后期數(shù)；若k=3 時，不僅描述了yt和yt-3 之間的相關(guān)性，還描述了yt-1 和yt-2 之間的相關(guān)性。自相關(guān)函數(shù)ACF(k)的取值范圍為［-1,1］，1 為正相關(guān)，-1為負(fù)相關(guān)，0 為不相關(guān)。

偏自相關(guān)函數(shù)PACF 描述的是預(yù)測值與實際值之間的線性相關(guān)性。當(dāng)k=3 時，描述的是yt和yt-3之間的相關(guān)性，剔除了yt-1 和yt-2 的影響。

（3）模型參數(shù)估計。ARIMA 模型選擇參數(shù)時，赤池信息準(zhǔn)則(AIC)可以彌補(bǔ)ACF 和PACF 定階的主觀性。AIC的定義為：

其中L為最大似然函數(shù)，N為數(shù)據(jù)個數(shù)，k為變量個數(shù)。參數(shù)的選擇需運用赤池信息準(zhǔn)則(AIC)進(jìn)行比較擇優(yōu)，在同樣本的模型估計中，AIC值越小模型越好，說明該模型雖然使用的參數(shù)較少但擬合度卻足夠。

（4）模型檢驗。若要對時間序列進(jìn)行白噪聲檢驗，需把非平穩(wěn)序列經(jīng)過n次差分運算變?yōu)槠椒€(wěn)序列。若白噪聲檢驗P＞0.05，證明該時間序列擬合效果良好。

（5）模型應(yīng)用。運用上述得到的最優(yōu)模型進(jìn)行預(yù)測，將預(yù)測值與實際值進(jìn)行對比，計算平均相對誤差。

圖1 ARIMA 建模過程示意圖Fig. 1 Schematic diagram of ARIMA modeling process

1.4 ARIMA—GRNN 組合模型的構(gòu)建

廣義回歸神經(jīng)網(wǎng)絡(luò)(GRNN)是在1991 年被Specht正式提出的基于非線性回歸分析的前饋型神經(jīng)網(wǎng)絡(luò)模型。GRNN 針對較少的時間序列數(shù)據(jù)，具有良好的預(yù)測效果，且在處理非線性數(shù)據(jù)方面優(yōu)于ARIMA 模型。通過建立ARIMA-GRNN 組合模型既能識別線性數(shù)據(jù)也能識別非線性數(shù)據(jù)，理論上可以提高模型的預(yù)測精度。ARIMA-GRNN 神經(jīng)網(wǎng)絡(luò)在糧食產(chǎn)量預(yù)測的運算中，主要有3 個步驟［18］：

（1）歸一化處理：農(nóng)業(yè)數(shù)據(jù)數(shù)量級復(fù)雜多樣，需把數(shù)據(jù)歸至［0,1］之間，避免因數(shù)據(jù)數(shù)量級差別較大導(dǎo)致GRNN 神經(jīng)網(wǎng)絡(luò)預(yù)測誤差增大。歸一化公式如下：

x為原始數(shù)據(jù)，xmin和xmax分別為原始數(shù)據(jù)的最小值、最大值，xi為歸一化處理后的數(shù)據(jù)值。

（2）確定數(shù)據(jù)樣本，包括輸入樣本、輸出樣本、訓(xùn)練樣本和預(yù)測樣本。本研究的輸入樣本為保定市1996—2014 年糧食產(chǎn)量ARIMA 模型的擬合值，輸出樣本為該時間段的實際值。

（3）構(gòu)建GRNN 神經(jīng)網(wǎng)絡(luò)模型：在模型構(gòu)建的過程中，需確定光滑因子spread這個參數(shù)。

2 結(jié)果與分析

2.1 ARIMA 模型的建立

2.1.1 保定市糧食產(chǎn)量時間序列分布 對保定市

1996—2017 年糧食產(chǎn)量（圖2）數(shù)據(jù)進(jìn)行單位根檢驗ADF，結(jié)果為ADF=0.025，證明該序列為非平穩(wěn)序列。

圖2 保定市1996—2017 年糧食產(chǎn)量Fig.2 Food production in Baoding city from 1996 to 2017

2.1.2 差分運算 由于時間序列為非平穩(wěn)序列，需對其進(jìn)行差分運算使之轉(zhuǎn)換為平穩(wěn)序列。將原序列進(jìn)行2 階差分后（圖3），均值圍繞0 不斷波動，可以看出經(jīng)過二階差分的時間序列為平穩(wěn)序列。

圖3 2 階差分圖Fig.3 2nd order difference diagram

2.1.3 模型定階 根據(jù)圖4 可以看出，該時間序列的自相關(guān)函數(shù)(ACF)值在1 階后均沒有超過置信區(qū)間（藍(lán)色區(qū)域）；偏自相關(guān)函數(shù)(PACF)值在2 階后迅速衰減至接近0，可以考慮p=0 或1、q=0、1 或2。

圖4 差分后糧食產(chǎn)量數(shù)列的ACF 及PACF 圖Fig. 4 ACF and PACF of the grain output sequence after difference

2.1.4 模型參數(shù)優(yōu)化 建立ARIMA 模型后，參數(shù)根據(jù)0、1、2 從低階到高階逐個選擇，需運用AIC信息準(zhǔn)則進(jìn)行比較擇優(yōu)，當(dāng)AIC值最小時模型最優(yōu)。由表1 可以看出，ARIMA(1,2,1)的AIC值小于其余模型，因此ARIMA(1,2,1)為選出的最優(yōu)模型。

表1 最優(yōu)模型的選取Table 1 Selection of optimal model

2.1.5 模型驗證 將ARIMA(1,2,1)進(jìn)行白噪聲檢驗后，模型的p值=0.828＞0.05（表2）。

表2 參數(shù)檢驗Table 2 Parameter test

證明ARIMA(1,2,1)模型對該時間序列擬合效果良好。

2.1.6 模型預(yù)測 運用上述得到的ARIMA(1,2,1)模型對保定市1996—2014 年糧食產(chǎn)量進(jìn)行擬合，2015—2017 年糧食產(chǎn)量進(jìn)行預(yù)測，模型的預(yù)測結(jié)果和相對誤差如表3 所示。

表3 保定市糧食產(chǎn)量預(yù)測值Table 3 Predicted value of grain output in Baoding City

2.2 糧食產(chǎn)量預(yù)測ARIMA-GRNN 組合模型的建立

運用上述ARIMA 模型預(yù)測糧食產(chǎn)量各影響因素的變化趨勢并結(jié)合GRNN 多元回歸模型對糧食產(chǎn)量進(jìn)行預(yù)測。

2.2.1 ARIMA-GRNN 模型輸入的確定 從《保定市經(jīng)濟(jì)統(tǒng)計年鑒》中選取16 個候選輸入變量，包括保定市受災(zāi)面積、農(nóng)業(yè)機(jī)械總動力、農(nóng)林牧漁業(yè)總產(chǎn)值、當(dāng)年實際機(jī)耕地面積、當(dāng)年機(jī)械播種面積、糧食播種面積、糧食播種單產(chǎn)、小麥播種面積、小麥播種單產(chǎn)、小麥總產(chǎn)量、玉米播種面積、玉米播種單產(chǎn)、玉米總產(chǎn)量、農(nóng)用化肥施用量、農(nóng)藥使用量、有效灌溉面積。因此當(dāng)年機(jī)械播種面積、小麥總產(chǎn)量、玉米播種面積、農(nóng)用化肥施用量、玉米總產(chǎn)量（圖5）為使用皮爾遜相關(guān)性分析（圖6）篩選出的糧食產(chǎn)量預(yù)測的5 個“最佳”輸入變量作為GRNN 神經(jīng)網(wǎng)絡(luò)的輸入。

圖5 影響因子Fig. 5 Influence factor

圖6 變量皮爾遜相關(guān)系數(shù)矩陣Fig. 6 Variable Pearson correlation coefficient matrix

2.2.2 ARIMA-GRNN 建模 本研究將1996—2014年糧食產(chǎn)量5 個影響因子作為輸入樣本，將2015—2017 年糧食產(chǎn)量預(yù)測值作為輸出樣本。選擇2015—2017 年保定市糧食產(chǎn)量預(yù)測值和實際值作為測試樣本。通過對GRNN 模型進(jìn)行訓(xùn)練，確定的最優(yōu)光滑因子spread=0.02。

2.3 預(yù)測模型分析

為了與ARIMA 模型、ARIMA-GRNN 組合模型進(jìn)行性能比較，將LSTM 模型也應(yīng)用于相同的訓(xùn)練和測試數(shù)據(jù)，結(jié)果列于表4。結(jié)果表明，GRNN 模型預(yù)測2015、2016、2017 年保定市糧食產(chǎn)量分別為503.8、505.6、503.7 萬t，ARIMA-GRNN 組合模型預(yù)測的的平均相對誤差為0.47%；ARIMA 模型預(yù)測的平均相對誤差為0.96%；LSTM 模型預(yù)測的平均相對誤差為2.20%。相比之下，ARIMA-GRNN 組合模型在預(yù)測方面表現(xiàn)優(yōu)異。

表4 模型擬合及預(yù)測誤差Table 4 Model fitting and prediction error

3 結(jié)論與討論

本研究使用1 組糧食產(chǎn)量影響因素和3 個機(jī)器學(xué)習(xí)模型分析了數(shù)據(jù)驅(qū)動的糧食作物產(chǎn)量預(yù)測方法的性能。預(yù)測結(jié)果表明，ARIMA-GRNN 組合模型可以作為糧食作物產(chǎn)量預(yù)測的有利方法。本研究通過對已有文獻(xiàn)的研究發(fā)現(xiàn)，在眾多學(xué)者對糧食產(chǎn)量的預(yù)測中，陳鼎玉等人利用ARIMA 模型預(yù)測糧食產(chǎn)量，將預(yù)測誤差控制在4.06%左右［19］。郭亞菲,樊超, 閆洪濤通過主成分分析和粒子群優(yōu)化神經(jīng)網(wǎng)絡(luò)組合模型，預(yù)測我國糧食產(chǎn)量的平均相對誤差在2%左右［20］。本研究利用Python 對保定市糧食產(chǎn)量預(yù)測進(jìn)行研究發(fā)現(xiàn)，ARIMA-GRNN 組合模型在準(zhǔn)確性和穩(wěn)定性方面都優(yōu)于其他算法，而ARIMA模型和LSTM 模型則沒有優(yōu)勢，尤其是對于小特征空間的數(shù)據(jù)集。ARIMA-GRNN 組合模型較好地擬合了糧食產(chǎn)量的動態(tài)發(fā)展過程，并預(yù)測了保定市近3 年糧食產(chǎn)量的發(fā)展趨勢。不僅為糧食生產(chǎn)研究提供了準(zhǔn)確有效的預(yù)測方法，還為其他領(lǐng)域的預(yù)測研究提供借鑒。

本研究的不足之處在于，僅選取了16 個影響糧食產(chǎn)量的因子，未考慮到氣溫、降水、日照時間、土壤濕度等因素對糧食產(chǎn)量的影響，以更準(zhǔn)確、全面地反映多種因素對糧食產(chǎn)量的影響。這種方法有助于多尺度全方面地對糧食產(chǎn)量進(jìn)行準(zhǔn)確的預(yù)測，可進(jìn)一步從數(shù)據(jù)中挖掘經(jīng)濟(jì)效益。

在未來的工作中，本課題組將擴(kuò)展各種數(shù)據(jù)集，更新氣候預(yù)測模型和生長階段模型，并進(jìn)一步檢驗所提出模型的有效性和準(zhǔn)確性。