徐鏡善 ，王 凱 ，袁哲明

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

1 原理與方法

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

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

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

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

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

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

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

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

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

1.3 因變量與自變量定階

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

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

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

1.5 一步預(yù)測(cè)及模型評(píng)價(jià)指標(biāo)

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

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

1.6 參比模型

為了體現(xiàn)研究方法的優(yōu)勢(shì)，需要與已有的模型進(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ù)測(cè)

數(shù)據(jù)來源于江蘇省通州市1973~1997年共25 a的西太平洋副高及稻縱卷葉螟二代高峰日每66.7m2蛾量Y[17]，對(duì)應(yīng)的預(yù)報(bào)自變量為副高的面積指數(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ú)立測(cè)試集。各年份的觀測(cè)值及各模型預(yù)測(cè)值見表2，均方根誤差RMSE預(yù)測(cè)指標(biāo)見表3。由表3可知，該研究預(yù)測(cè)所得RMSE值在所有參比模型中最小，說明所建模型穩(wěn)定性及外部預(yù)測(cè)能力最好。由表2可知，除了1988年、1993年、1994年這3 a的蛾量預(yù)測(cè)值與觀測(cè)值偏差較大，其他年份的預(yù)測(cè)偏差均較??；對(duì)1988年的蛾量預(yù)測(cè)，MLR與SLR線性模型的預(yù)測(cè)值比基于SVR建立的非線性模型預(yù)測(cè)值更加接近觀測(cè)值，說明了當(dāng)年數(shù)據(jù)各因素間的線性本質(zhì)；對(duì)1993年的蛾量預(yù)測(cè)，MLR與SLR預(yù)測(cè)出了負(fù)值，表明當(dāng)年數(shù)據(jù)仍用線性模型預(yù)測(cè)是不合理的。對(duì)最后3年即1995～1997年的蛾量預(yù)測(cè)，該研究方法比其他參比模型的預(yù)測(cè)偏差大幅減小。因此，該研究方法所建模型在整體上優(yōu)于其他參比模型，顯示了因變量快速定階以及一步預(yù)測(cè)法的優(yōu)勢(shì)。

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

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

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

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

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

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

3 結(jié)論

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

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

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

[3]李哈濱，王政權(quán)，王慶成.空間異質(zhì)性定量研究理論與方法[J].應(yīng)用生態(tài)學(xué)報(bào)，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]鄧自立，郭一新.動(dòng)態(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ù)報(bào)大氣污染物濃度[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]王海燕，盧 山.非線性時(shí)間序列分析及其應(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ì)分析比較研究[J].應(yīng)用生態(tài)學(xué)報(bào)，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)計(jì)學(xué)與支持向量回歸的QSAR 建模[J].物理化學(xué)學(xué)報(bào)，2009，25（8）：1587-1592.

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

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