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

徐鏡善 ，王 凱 ，袁哲明

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

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

1 原理與方法

1.1 數(shù)據(jù)預處理

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

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

以平穩(wěn)化處理后的因變量yt與規(guī)格化后的自變量xi進行后續(xù)分析。

1.2 地統(tǒng)計學后效時間長度

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

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

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

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

1.3 因變量與自變量定階

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

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

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

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

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

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

1.6 參比模型

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

2 結(jié)果與分析

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

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

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

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

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

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

2.2 晚稻第五代褐飛虱發(fā)生量預測

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

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

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

3 結(jié)論

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

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

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

[3]李哈濱，王政權(quán)，王慶成.空間異質(zhì)性定量研究理論與方法[J].應用生態(tài)學報，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]鄧自立，郭一新.動態(tài)系統(tǒng)分析及其應用[M].沈陽：遼寧科學技術出版社，1985.31-130.

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

[9]鄧乃揚，田英杰.數(shù)據(jù)挖掘中的新方法－支持向量機[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īng)典分析與地統(tǒng)計分析比較研究[J].應用生態(tài)學報，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]陳 淵，袁哲明，周 瑋，等.基于地統(tǒng)計學與支持向量回歸的QSAR 建模[J].物理化學學報，2009，25（8）：1587-1592.

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

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