徐鏡善 ，王 凱 ，袁哲明

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

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

1 原理與方法

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

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

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

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

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

地統(tǒng)計學(xué)主要研究在空間分布上既有結(jié)構(gòu)性又有隨機(jī)性的自然現(xiàn)象，其以區(qū)域化變量理論為基礎(chǔ)，以半變異函數(shù)為主要工具[3-4]。對于多維時間序列因變量觀測值y(t)，t=1，…，，其實(shí)驗(yàn)半變異函數(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ù)曲線?？捎们驙钅Ｐ蛿M合：

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

1.3 因變量與自變量定階

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

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

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

1.5 一步預(yù)測及模型評價指標(biāo)

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

模型的獨(dú)立預(yù)測精度采用均方根誤差（Root Mean Square Error，RMSE）作為評價指標(biāo)：

1.6 參比模型

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

2 結(jié)果與分析

2.1 稻縱卷葉螟二代高峰日蛾量預(yù)測

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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

[9]鄧乃揚(yáng)，田英杰.數(shù)據(jù)挖掘中的新方法－支持向量機(jī)[M].北京：科學(xué)出版社，2004.

[10]馬曉光，胡 非.利用支持向量機(jī)預(yù)報大氣污染物濃度[J].自然科學(xué)進(jìn)展，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]王海燕，盧 山.非線性時間序列分析及其應(yīng)用[M].北京：科學(xué)出版社，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].應(yīng)用生態(tài)學(xué)報，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)計學(xué)與支持向量回歸的QSAR 建模[J].物理化學(xué)學(xué)報，2009，25（8）：1587-1592.

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

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