多維時間序列因變量的快速定階及其害蟲發生量預測

徐鏡善 ，王 凱 ，袁哲明

（1.湖南大眾傳媒職業技術學院，湖南 長沙 410100；2.湖南農業大學植物保護學院，湖南 長沙410128；3.湖南農業大學作物種質創新與資源利用國家重點實驗室，湖南 長沙 410128）

害蟲發生量不僅受多種外在因素（如氣象，生理，生態等）影響，而且與其歷年發生動態極其相關，屬于典型的多維時間序列復雜非線性數據[1-2]。傳統自回歸模型（Autoregressive，AR）假定當年觀察值是過去p年（p為階次）觀察值的線性組合，但僅能用于一維時間序列分析且無法給出拓階上限，而地統計學（Geostatistics，GS）中半變異函數模型的后效時間長度可實現一維時間序列自動快速定階[3-4]。帶控制項的自回歸滑動平均模型（Controlled Autoregressive Integrating Moving Average，CARMA）及帶受控項的自回歸模型（Controlled Autoregressive，CAR）等已應用于多個多維時間序列分析，但由于其線性本質，導致應用受限[5-7]。基于結構風險最小的支持向量機（Support Vector Machine，SVM）以統計學習理論為基礎，較好地解決了局部最小、過學習、非線性等問題，有效應用于復雜非線性數據系統的建模預測[8-11]。研究結合地統計學半變異函數模型與支持向量回歸（Support Vector Regression，SVR），實現多維時間序列的因變量自動定階與自變量非線性篩選，并應用于2種害蟲發生量預測。

1 原理與方法

1.1 數據預處理

多維時間序列因變量yt若具有明顯的上升或下降趨勢，則需要對其進行去趨勢平穩化處理。研究采用對數線性去趨勢（Log-linear De-trending，LLD）平穩化法對因變量進行平穩化處理[12]：首先，對因變量yt取自然對數得lnyt；其次，以lnyt與時間t建立一元線性回歸方程lnyt=a+bt，求解得截距a與回歸系數b，以時間t為自變量回帶可得預測(Inyt)′；最后，對數線性去趨勢平穩化后的因變量為：

由于支持向量機程序LIBSVM[13]對自變量取值較為敏感，所以將多維時間序列自變量根據式（2）規格化到-1到1的區間：

以平穩化處理后的因變量yt與規格化后的自變量xi進行后續分析。

1.2 地統計學后效時間長度

地統計學主要研究在空間分布上既有結構性又有隨機性的自然現象，其以區域化變量理論為基礎，以半變異函數為主要工具[3-4]。對于多維時間序列因變量觀測值y(t)，t=1，…，，其實驗半變異函數值r(h)可用下式計算：

式（3）中，h是兩個觀測值之間的時間間隔，N(h)是相隔時間為h的數據對y(t)和y(t+h)的對數，y(t)和y(t+h)分別是時間t和時間t+h的觀測值。假定兩個觀測值間最大時間間隔為max(d)，為保證N(h)充分大，一般規定h≤max(d)。

半變異函數曲線常可用球狀模型擬合：

式中，c0為塊金值，c為基臺值，a為變程，即r(h)達到基臺值時的間隔距離h，它表示在h≥a以后，區域化變量的空間相關性消失。因此，對某個待測點的估計，應根據其距離a以內的已知點來進行。袁哲明等[14-15]將半變異函數模型應用于一維等間隔時間序列，研究了二化螟與三化螟種群的時間格局。半變異函數模型的變程a對應于時間序列分析的后效時間長度，根據變程a可實現多維時間序列因變量的快速定階。

1.3 因變量與自變量定階

設已平穩化處理的多維時間序列為（yt，xij），t=1，2，…，n；j=1，2，…，m；其中 n 為樣本數，m 為自變量個數。首先，根據CAR對多維時間序列原始自變量拓且僅拓一階，則訓練集樣本個數變為n-1，自變量個數變為2m；其次，對多維時間序列因變量計算半變異函數，為避免擬合半變異函數曲線帶來的擬合誤差及額外計算量，設定的時間間隔范圍內最大半變異函數值r(h)對應的h即為后效時間長度a，因變量拓展a階后訓練集樣本個數變為n-1-a，自變量個數變為2m+a。

1.4 支持向量回歸非線性變量篩選

SVR核函數的選擇一般是經驗性的，徑向基核函數在多數數據集上比其他核函數（如線性核、二項式核等）表現優異[16]，因此研究采用徑向基核作為SVR建模核函數。對于訓練集（yt，xij），t=1，2，…，n-1-a；j=1，2，…，2m+a，首先對 2m+a 個自變量經10次交叉驗證搜尋最優SVR參數，對應均方誤差（Mean Squared Error，MSE）記為 MSE0；其次，依次剔除第 j個自變量得 MSEj，若 min[MSEj]＜MSE0，則剔除最小MSE對應的自變量，進入下一輪篩選（此時自變量個數變為2m+a-1），反之，非線性變量篩選過程結束；最后，保留自變量用于后續SVR建模。為方便后文參比，將該自變量篩選方法命名為Support Vector Regression-Nonlinear Variable Screening（SVR-NVS）。

1.5 一步預測及模型評價指標

此時訓練集為（yt，xij），t=1，2，…，n-1-a；j=1，2，…，s；s為保留自變量個數。以徑向基核函數為基礎，對訓練集經10次交叉驗證搜尋最優參數，根據確定的最優參數建立SVR模型。對第1個待測樣本，依訓練集確定的階次與保留自變量構建新的測試樣本，以SVR模型預測該樣本，根據式（1）反推得最終預測值。預測第2個待測樣本時，第1個測試樣本需加入到訓練集中，此為一步預測。由于訓練集有改變，所以需要重新拓階及變量篩選。

模型的獨立預測精度采用均方根誤差（Root Mean Square Error，RMSE）作為評價指標：

1.6 參比模型

為了體現研究方法的優勢，需要與已有的模型進行參比。研究參比了4種模型，其中包括基于SVR的非線性模型SVR-NVS與SVR，參比線性模型包括多元線性回歸（Multiple Linear Regression，MLR）與逐步線性回歸（Stepwise Linear Regression，SLR）。所有模型均在MATLAB（2012a）環境下實現。

2 結果與分析

2.1 稻縱卷葉螟二代高峰日蛾量預測

數據來源于江蘇省通州市1973~1997年共25 a的西太平洋副高及稻縱卷葉螟二代高峰日每66.7m2蛾量Y[17]，對應的預報自變量為副高的面積指數（X1），強度指數（X2），西伸脊點位置（X3，經度），脊線位置（X4，經度），北界位置（X5，緯度）（表 1）。

表1 稻縱卷葉螟二代高峰日蛾量

研究以1973～1987年15 a的數據為訓練集，1988～1997年間10 a數據作為獨立測試集。各年份的觀測值及各模型預測值見表2，均方根誤差RMSE預測指標見表3。由表3可知，該研究預測所得RMSE值在所有參比模型中最小，說明所建模型穩定性及外部預測能力最好。由表2可知，除了1988年、1993年、1994年這3 a的蛾量預測值與觀測值偏差較大，其他年份的預測偏差均較小；對1988年的蛾量預測，MLR與SLR線性模型的預測值比基于SVR建立的非線性模型預測值更加接近觀測值，說明了當年數據各因素間的線性本質；對1993年的蛾量預測，MLR與SLR預測出了負值，表明當年數據仍用線性模型預測是不合理的。對最后3年即1995～1997年的蛾量預測，該研究方法比其他參比模型的預測偏差大幅減小。因此，該研究方法所建模型在整體上優于其他參比模型，顯示了因變量快速定階以及一步預測法的優勢。

表2 稻縱卷葉螟二代高峰日蛾量觀測值及預測值

表3 參比模型均方根誤差RMSE值

2.2 晚稻第五代褐飛虱發生量預測

第五代褐飛虱是晚稻生長期間的主要害蟲，提高褐飛虱發生量預測準確度，及時采取有效防治措施，對于晚稻豐收至關重要。數據來源于浙江省杭州市蕭山區1974～2005年晚稻第五代褐飛虱的發生量[18]。

以1974～2000年多維時間序列數據為訓練集，2001～2005年為獨立測試樣本，各年份的觀測值及各個模型預測值見表4，均方根誤差RMSE值見表3。由表4可知，該研究方法對各年份的預測精度明顯優于其他模型。對2001年褐飛虱的發生量預測值比其他SVR非線性模型的預測值精確了將近一個數量級，也優于SLR模型，與MLR模型預測準確度相當；對2002年、2003年的褐飛虱發生量進行預測，該研究方法預測值較其他4個參比模型明顯接近真實觀測值。

表4 晚稻第五代褐飛虱發生量觀測值及預測值

3 結論

以地統計學半變異函數模型的后效時間長度為指導，實現了多維時間序列因變量的自動、快速定階，并結合SVR非線性篩選自變量，以一步預測法成功應用于2種害蟲發生量的預測。該方法具有地統計學半變異函數模型對因變量定階自動化、快速的優點，適合時間序列分析；同時具有SVR結構風險最小、非線性、有效避免過擬合與局部最小、外部預測能力優異等優點，適合多因素非線性回歸分析。經試驗表明，該研究方法可有效應用于害蟲預測預報。

[1]吳承禎，洪 偉.林木生長的多維時間序列分析[J].應用生態學報，1999，10（4）：395-398.

[2]周立陽，費惠新，張孝羲.多維時間序列分析在稻縱卷葉螟長期預測預報上的試用[J].植物保護學報，1995，22（1）：1-6.

[3]李哈濱，王政權，王慶成.空間異質性定量研究理論與方法[J].應用生態學報，1998，9（6）：651-657.

[4]Liebhold A M，Rossi R E，Kemp W P.Geostatistics and geographic information systems in applied insect ecology[J].Annu.Rev.Entomo，1993，（38）：303-327.

[5]Box Q E P，Jenkins G M.Time series analysis：forecasting and control[M].San Francisco：Holden-day Press，1970.

[6]Hannan E J.The estimation of the order of an ARMA process[J].Ann.Statist，1980，8（5）：1071-1081.

[7]鄧自立，郭一新.動態系統分析及其應用[M].沈陽：遼寧科學技術出版社，1985.31-130.

[8]V Vapnik.The nature of statistical learning theory[M].New York:Springer Verlag Press，2001.

[9]鄧乃揚，田英杰.數據挖掘中的新方法－支持向量機[M].北京：科學出版社，2004.

[10]馬曉光，胡 非.利用支持向量機預報大氣污染物濃度[J].自然科學進展，2004，14（3）：349-353.

[11]Ping F P，Wei C H.Support vector machines with simulated annealing algorithms in electricity load forecasting[J].Energy Conversion and Management，2005，46：2669-2688.

[12]王海燕，盧 山.非線性時間序列分析及其應用[M].北京：科學出版社，2006.139-141.

[13]Chang C C，Lin C J.LIBSVM:a library for support vector machines[J].ACM Transactions on Intelligent Systems and Technology(TIST)，2011，2（3）：27.

[14]袁哲明，付 威，李方一.二化螟種群空間格局的經典分析與地統計分析比較研究[J].應用生態學報，2004，15（4）：610-614.

[15]Yuan Z M，Wang Z，Hu X Y.Geostatistical analysis on the temporal patterns of the Yellow Rice Borer.Tryporyza incertulas[J].Rice Science，2005，12（3）：207-212.

[16]陳 淵，袁哲明，周 瑋，等.基于地統計學與支持向量回歸的QSAR 建模[J].物理化學學報，2009，25（8）：1587-1592.

[17]汪四水，張孝羲.害蟲預測建模中的因子綜合[J].南京農業大學學報，2000，23（2）：35-38.

[18]陳水校.簡易多級法預測晚稻第五代褐飛虱發生量[J].江蘇農業科學，2006，（4）：48-50.