基于時間序列與橫截面數據的吉林省水稻產量預測對比分析*

陳 威，祁偉彥，袁福香，李哲敏※

（1. 中國農業科學院農業信息研究所，北京100081；2. 吉林省氣象科學研究所，長春130062）

0 引言

水稻是世界三大糧食作物之一，是我國最主要的糧食作物，全國65%以上的人口以稻米為主食。我國水稻播種面積占全國糧食播種面積的27%左右，產量占全國糧食總產量的37%左右［1］。因此，水稻的穩產增產以及對水稻產量的準確預測對我國農業政策調整和保障我國糧食安全問題具有重要意義。

國內學者在糧食產量預測方面做了大量研究，有的利用糧食產量的時間序列數據進行預測，有的利用與糧食生產密切相關的多重參數的橫截面數據進行預測。時間序列數據預測常用的預測算法有自回歸滑動平均、人工神經網絡等。由于時間序列預測只需變量本身的歷史數據，不需要其他參數的數據，模型構建較為簡單，因此在糧食產量預測中應用十分廣泛。然而，時間序列預測假設糧食產量具有線性變化的規律，但糧食生產具有受多維度因素影響的復雜性，使其呈現非線性變化。因此，單純采用時間序列的方法進行糧食產量預測，預測精度不高。橫截面數據預測常用的方法有支持向量回歸、神經網絡等。在利用橫截面數據進行糧食產量預測以及支持向量回歸構建模型時參數一般依靠經驗來確定，而神經網絡訓練模型對網絡初始權重高度敏感，且對訓練集樣本的依賴性強，因此預測結果有時并不理想。

水稻產量受自然因素、生產技術因素、社會經濟因素和隨機因素等多因素影響，是線形關系和非線性關系并存的復雜性問題。使用單一模型進行水稻產量預測，當影響產量的重要因素發生變化時可能導致預測結果的不穩定。為了對比不同模型預測效果，分析各模型預測水稻產量的特點、不足及適用條件，為糧食產量預測問題模型選擇提供依據，該研究從時間序列預測和橫截面數據預測兩種角度，選取4種模型對水稻產量進行預測。在時間序列預測方面，依據吉林省水稻產量的歷史數據，基于自回歸滑動平均（Autoregressive Integrated Moving Average Model，ARIMA）模型和長短期記憶人工神經網絡（Long-Short Term Memory，LSTM）模型進行吉林省水稻產量預測。ARIMA模型在糧食產量預測中應用廣泛［2-3］，而LSTM模型屬于反饋神經網絡的范疇，在時間序列的學習和處理上有獨特優勢，在語音及圖像識別、機器翻譯等方面已有較多研究［4-6］，但在糧食產量預測中應用較少［7］。由于農作物生長產生過程中動物害蟲和病原體的侵害可引起顯著的產量潛在損失［8-9］，競爭性作物品種種植面積對水稻的產量也有重要影響，因此選取吉林省重要病蟲害發生情況、玉米和大豆種植面積產量等歷史數據作為橫截面數據預測的原始數據。在橫截面數據預測模型選擇上，選取支持向量回歸模型和多層感知器模型進行預測。支持向量回歸模型在解決小樣本問題中具有優勢，已有學者采用此模型進行農作物產量的預測［10］，而多層感知器模型作為最簡單的神經網絡模型，具有良好的容錯性和自適性，可有效處理非線性可分離問題，在語音識別、機器翻譯等方面已有大量研究［11-13］，用于糧食產量預測的報道不多。

1 模型簡介

1.1 自回歸積分滑動平均（ARIMA）模型

ARIMA模型全稱為自回歸積分滑動平均模型，是指將非平穩時間序列轉化為平穩時間序列，然后將因變量僅對它的滯后值以及隨機誤差項的現值和滯后值進行回歸所建立的模型。ARIMA模型的3個參數為p、d、q，其中p代表預測模型中采用的時序數據本身的滯后數，也叫做自回歸項數，d代表時間序列成為平穩時所作的差分次數，q代表預測模型中采用的預測誤差的滯后數，也叫做移動平均項數［14］。

1.2 長短期記憶人工神經網絡（LSTM）模型

長短期記憶人工神經網絡是一種改進的循環神經網絡（Recurrent Neural Network，RNN），避免了傳統RNN隱含層層數過多計算量龐大而產生的梯度消失和梯度爆炸的問題。LSTM的經典結構如圖1所示，Cell是計算節點的核心，通過遺忘閥門（forget gate）、輸入閥門（input gate）和輸出閥門（output gate）的打開或者關閉來判斷模型網絡的記憶態（之前網絡的狀態）在該層輸出的結果是否達到閾值，從而加入到當前該層的計算中［15］。

圖1 長短期記憶人工神經網絡（LSTM）模型結構Fig.1 Model structure of LSTM（Long-Short Term Memory）

1.3 支持向量回歸（SVR）

支持向量機（Support Vector Machine，SVM）是機器學習領域中基于統計學習理論的分類回歸方法，由AT&T貝爾實驗室的Vanpnik于20世紀90年代提出，通過構造核函數將原空間非線性問題轉換成高維空間的線性問題，可以成功處理分離問題，辨別分析回歸問題等［16］。在處理分類問題時，支持向量機相當于找到一個高維特征平面，使兩個分類集合的支持向量或者所有數據離分類平面最遠，而支持向量回歸（Support Vector Regression，SVR）則是找到一個高維回歸平面，使一個集合中的所有數據到該平面的距離最近。

1.4 多層感知器MLP模型

MLP（Multilayer Perceptron，MLR），即多層感知器，是一種常見的神經網絡模型，由輸入層、隱藏層和輸出層組成，其基本結構如圖2所示。MLP是一種監督學習模型，在模型訓練過程中，不斷提供完整的輸入和輸出，模型基于一定長度的歷史數據不斷訓練、學習，從而構建出輸入和期望的輸出之間的最優模型，而后利用得出的最優模型進行預測［17］。

圖2 多層感知器（MLP）基本結構Fig.2 Basic structure of MLP（Multilayer Perceptron）

2 水稻產量預測實證分析

2.1 數據來源

研究中所用到的實驗數據為1949—2015年吉林省玉米、大豆及水稻的種植面積、單產及總產量數據，以上數據來自吉林省統計局；吉林省稻瘟病發病頻率（1961—2009年）、玉米螟發生情況（1980—2007年）、地下害蟲及大豆蚜發生情況（1981—2006年）、食心蟲發病情況（1982—2006年）、黏蟲發生情況（1979—2006年），以上數據來自吉林省氣象局。

以1990年為基期，1990年之前的數據作為訓練集數據用于模型訓練，1990年之后的數據作為預測集數據，用于檢測模型預測效果。

2.2 評價指標

使用均方根誤差值（Root Mean Square Error，RMSE）作為評價指標對各模型模擬效果進行評價和比較。RMSE是觀測值與真值偏差的平方和觀測次數n比值的平方根，計算方法見公式（1）。RMSE值越小，表明預測值與擬合值之間偏差越小，預測性能越好，反之則預測性能越差。

2.3 從時間序列角度預測產量

（1）基于ARIMA的時間序列預測

①時間序列平穩性檢驗及平穩化處理

導入吉林省1949—1990年歷年水稻總產量原始數據，并檢驗時間序列的穩定性。若時間序列的統計學特性（例如均值、方差）隨時間保持不變，則認為該序列是平穩的。

圖3 1950—1990年吉林省水稻總產量Fig.3 Total rice yield in Jilin province from 1950 to 1990

如圖3所示，1949—1990年吉林省水稻總產量呈上升趨勢，年增長幅度不同。由移動平均值和移動標準方差曲線可以看出，移動標準方差隨時間變化較小，但移動平均值隨時間呈明顯的上升趨勢，因此該序列不是平穩序列，需要進行平穩化處理。

由于原數據值域范圍比較大，為了縮小值域，同時保留其他信息，對原始數據進行對數轉換，之后采用差分的方法，將每一年的數值與滯后10年的值作差，去除數據的趨勢性（圖4）。

圖4 原始數據平穩化處理Fig.4 Stationary process of raw data

經對數轉換和差分處理后，數據的移動均值和移動標準方差在時間軸上振幅明顯縮小，DFtest結果顯示Test Statistic的值小于Critical Value（1%）的值，表明在99%的置信度下，該時間序列數據是穩定的，如圖5。

圖5 平穩化處理后移動均值與移動標準方差曲線Fig.5 Rolling mean and standard deviation after stationary process

②構建ARIMA模型

通過自相關函數（ACF）和部分自相關函數（PACF）進行ARIMA（p，d，q）的p，q參數估計。由前文差分部分已知，原始數據經一階差分后數據已經穩定，可得到d=1。所以用一階差分化的Y’t=Yt-Yt-1作為輸入。畫出ACF，PACF的圖像（圖6）。

圖6中，上下兩條灰線之間為置信區間，p的值為ACF第一次穿過上置信區間時的橫軸值。q的值為PACF第一次穿過上置信區間的橫軸值。所以從圖6可以得到p=2，q=2。由此，可生成模型ARIMA（2，1，2）。

取AR模 型（ARIMA（2，1，0））、MA模 型（ARIMA（0，1，2））、ARIMA模 型（ARIMA（2，1，2））這3種參數取值的3個模型與穩定化處理后的輸入數據進行擬合，比較RSS值得出，ARIMA模型（ARIMA（2，1，2））擬合度最好（圖7）。

③ARIMA模型預測

對ARIMA模型（ARIMA（2，1，2））擬合值進行差分和對數處理的逆運算（圖8a），使擬合值回到原始區間，并利用該模型對1991—2015年水稻產量進行預測（圖8b）。

圖6 自相關函數（ACF）和部分自相關函數（PACF）曲線Fig.6 Curve graph of ACF and PACF

圖7 3種參數ARIMA模型與平穩數據擬合度分析Fig.7 Fitting analysis of three parameters ARIMA model with stationary data

圖8 差分與對數的逆運算（a）及水稻產量預測（b）Fig.8 Inverse operation of difference and logarithm（a）and yield prediction of rice（b）

1949—2015年吉林省水稻產量整體呈上升趨勢，期間由于氣候因素和政策導向等導致一些年份水稻產量出現較大波動，特別是2005—2015年。ARIMA模型預測結果較好地反映了水稻產量整體上升的趨勢，但對原始數序列的波動沒有做出反應。

（2）基于LSTM的時間序列預測

①建立模型

使用Keras深度學習庫在Python中構建LSTM網絡模型。參數設置為默認參數，即輸入層有1個input，隱藏層4個神經元，輸出層為預測一個值，激活函數用sigmoid，迭代100次，batch size 為1。利用上面得到的平穩序列數據對模型進行訓練。

②模型預測

利用得到的LSTM模型對1991—2015年水稻總產量進行預測，結果如圖9。

圖9 長短期記憶人工神經網絡（LSTM）模型預測水稻產量Fig.9 Prediction of rice yield using LSTM model

可以看出，LSTM模型較預測曲線整體趨勢較為平緩，預測產量普遍低于實際產量，未能反應水稻產量整體上升的總體趨勢，對于原始數序列的波動反應不明顯。

2.4 從橫截面數據角度預測產量

實驗中用于預測水稻產量的病蟲害數據為吉林省歷年稻瘟病發病頻率（1961—2009年）、玉米螟發生情況（1980—2007年）、地下害蟲及大豆蚜發生情況（1981—2006年）、食心蟲發病情況（1982—2006年）、黏蟲發生情況（1979—2006年）。用于預測水稻產量的其他特征數據為水稻的競爭性作物種植面積與產量的歷年數據，包括1949—2015年吉林省玉米、大豆種植面積、單產及總產量數據。

實驗使用sklearn機器學習庫進行數據的處理和SVR、MLP模型訓練及預測。主要包括數據預處理、定義分類器、訓練模型、模型預測等步驟。由于原始數據由多種特征數據組成，為高維數據，為了更有效處理數據，對原始數據進行主成分分析（PCA）降維處理后，再次進行上述分析流程，并與未經降維處理的預測結果進行對比。

采用每種方法建立預測模型時，考慮兩種特征的組合：①考慮水稻生長受當年病蟲害影響，采用病蟲害數據預測當年水稻產量；②考慮競爭性作物種植面積及產量對農戶種植意向的影響，采用競爭性作物種植面積及產量數據預測下一年水稻產量。兩種特征組合可分別應用于當年及下一年水稻產量的預測任務。

（1）基于SVR的機器學習預測

對原始數據進行歸一化處理后，將1950—1990年的數據作為訓練集用于預測模型的建立，1991—2015年的數據作為測試集用于檢驗模型預測效果。

使用sklearn機器學習庫SVR模型對水稻產量進行預測。主要步驟為：對原始數據中的缺失值進行函數填充后，對數據進行歸一化處理；選擇分類器為SVR，用訓練集數據對模型進行訓練，并建立預測模型；將測試集數據于模型預測效果的檢驗。對原始數據進行主成分分析PCA降維處理后，再次進行上述分析流程。

根據以上流程，利用病蟲害數據預測當年水稻產量，得到預測曲線SVR1和SVR1（主成分分析降維）。利用競爭性作物種植面積及產量數據預測下一年水稻產量，得到預測曲線SVR2和SVR2（主成分分析降維）（圖10）。

圖10 利用支持向量機（SVR）模型預測水稻產量Fig.10 Prediction of rice yield using SVR model

SVR模型預測曲線呈波動上升，預測產量普遍低于實際產量，在一定程度上反應了原始數序列的波動情況。由分析結果可知，采用PCA對影響因素進行降維處理后，SVR模型預測準確性均有所提高。

（2）基于MLP神經網絡的機器學習預測

使用sklearn機器學習庫MLP模型對水稻產量進行預測。主要步驟為：對原始數據中的缺失值進行函數填充后，對數據進行歸一化處理；選擇分類器為MLR，用訓練集數據對模型進行訓練，并建立預測模型；將測試集數據于模型預測效果的檢驗。對原始數據進行主成分分析PCA降維處理后，再次進行上述分析流程。

根據以上流程，利用病蟲害數據預測當年水稻產量，得到預測曲線MLP1和MLP1（主成分分析降維）。利用競爭性作物種植面積及產量數據預測下一年水稻產量，得到預測曲線MLP2和MLP2（主成分分析降維）（圖11）。

圖11 利用多層感知器（MLP）模型預測水稻產量Fig.11 Prediction of rice yield using MLP model

MLP模型預測產量普遍低于實際產量，預測結果隨原始數序列的波動出現劇烈波動。采用PCA對影響因素進行降維處理后，預測曲線波動性減緩，MLP1模型預測準確性提高，而MLP2模型預測準確性降低。

3 結果與討論

利用吉林省病蟲害和玉米大豆產量等歷史數據，基于4種模型對吉林省水稻產量進行預測，預測值與實際值的RMSE見表1?；贏RIMA模型和LSTM模型的時間序列預測RMSE相對橫截面預測RMSE較小，但預測值與實際值之間仍存在差距。在橫截面數據預測中，原始數據經主成分分析PCA降維處理后，SVR1、SVR2和MLP1模型預測的RMSE值相對未經降維處理均有所下降，表明在數據維度過高、數據量較小情況下，降維處理在多數情況下可提升模型預測性能。

近年來我國對糧食生產實施了良種補貼、農機購置補貼、農民直補等一系列農業支持政策，加之生物技術、農業工程和管理技術的發展和農田水利基礎設施的不斷完善，使糧食持續增產潛力大大提升。利用ARIMA和LSTM進行時間序列預測時，只考慮產量隨時間的變化趨勢，沒有考慮其他因素的變化對趨勢延續性的影響，導致對近年產量的預測與實際情況相比偏低。

表1 不同方法預測結果RMSE值Table 1 RMSE of predicting results with different methods

農業生產是復雜的經濟與自然再生產過程，糧食生產過程受到氣象條件、作物品種、農資投入、田間管理、市場供需、國家政策等多方面的影響［18］。而利用SVR和MLP模型進行橫截面數據預測時，模型訓練使用的原始數據為水稻重要病蟲害發生情況、大豆種植面積與產量、玉米種植面積與產量等數據，致使模型主要反映水稻產量與病蟲害發生情況及競爭性作物種植情況的關系，沒有考慮其他影響產量的氣象因素、投入管理因素、社會經濟因素，導致模型預測產量與實際產量相比偏低。

由時間序列模型與橫截面數據預測模型對比的實證分析結果可知，在進行水稻產量預測時，如缺乏影響產量的重要因素的數據，則應用時間序列模型分析結果較好。反之，如掌握影響產量的重要因素數據，則宜采用橫截面數據預測模型。在進行橫截面數據預測時，如分析影響因素較多，可采用PCA對影響因素進行降維處理，通常能夠提高預測結果的準確性。此外，在進行橫截面數據預測中，如訓練數據較少時，采用回歸模型比神經網絡模型預測結果準確性高，且神經網絡模型中，在訓練數據量較小時，算法表現不夠穩定。在進行橫截面數據進行水稻產量預測時，利用競爭性作物種植面積及產量數據預測下一年水稻產量的結果優于利用病蟲害數據預測當年水稻產量。

在后續研究中，擬通過補充氣象數據、農資投入等數據，提高基于SVR和MLP模型的橫截面數據預測準確性，并通過構建時間序列預測和橫截面預測組合模型，進一步提高水稻產量預測的精度。