基于混沌理論的害蟲發(fā)生量非線性建模與預(yù)測

向昌盛，萬方浩

(1. 湖南工程學(xué)院計算機(jī)與通信系，湖南 湘潭 411104；2. 植物病蟲害生物學(xué)國家重點實驗室，北京 100193；3. 中國農(nóng)業(yè)科學(xué)院植物保護(hù)研究所，北京100193；4. 青島農(nóng)業(yè)大學(xué)農(nóng)學(xué)與植物保護(hù)學(xué)院，山東 青島266109)

害蟲災(zāi)害的發(fā)生是一種高度復(fù)雜的自然現(xiàn)象，人們對其規(guī)律進(jìn)行長期探索，積累了大量的歷史數(shù)據(jù)，利用這些歷史數(shù)據(jù)對害蟲發(fā)生趨勢進(jìn)行預(yù)測，能夠為害蟲防治提供有價值的參考意見；因此，害蟲發(fā)生量建模與預(yù)測具有重要的研究意義[1]。

傳統(tǒng)的害蟲預(yù)測模型主要有多元線性回歸、自回歸滑動平均、趨勢外推法等[2–3]。由于受到氣象、作物種類、栽培管理等因素的綜合影響，加上害蟲發(fā)生自身獨特的變化規(guī)律，使得害蟲發(fā)生具有規(guī)律性和周期性及不穩(wěn)定性，造成傳統(tǒng)的模型預(yù)測精度較低[4]。隨著人工智能技術(shù)發(fā)展，出現(xiàn)了基于神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)等機(jī)器學(xué)習(xí)算法的害蟲發(fā)生量預(yù)測模型，它們可以較好地描述害蟲的非線性發(fā)生趨勢，提高了預(yù)測精度[5–7]。在害蟲發(fā)生量建模過程中，神經(jīng)網(wǎng)絡(luò)基于結(jié)構(gòu)風(fēng)險最大原理，要求訓(xùn)練數(shù)量比較大，而害蟲發(fā)生量預(yù)測是一種小樣本預(yù)測問題，導(dǎo)致神經(jīng)網(wǎng)絡(luò)常出現(xiàn)“過擬合”、訓(xùn)練速度慢、陷入局部極小等缺陷；支持向量機(jī)雖然具有良好的泛化能力，但是訓(xùn)練時間相對較長，建模效率低，難以滿足害蟲發(fā)生預(yù)測實時性要求。最小二乘支持向量機(jī)(LSSVM)不僅克服了神經(jīng)網(wǎng)絡(luò)存在的缺陷，而且解決了支持向量機(jī)建模效率低的問題，泛化性能優(yōu)異[8–10]。大量研究結(jié)果表明，害蟲的發(fā)生具有弱混沌性，而當(dāng)前害蟲發(fā)生量預(yù)測建模過程中，害蟲發(fā)生的混沌性常常被忽略，導(dǎo)致害蟲預(yù)測精度有待提高[11]。近年來，隨著混沌理論研究的不斷深入，將混沌理論和機(jī)器學(xué)習(xí)方法相結(jié)合，為害蟲發(fā)生預(yù)測開辟了一種新的途徑[12]。

由于害蟲發(fā)生量的非線性、混沌、突變性等特點，至今沒有統(tǒng)一的預(yù)測模型。為了客觀、準(zhǔn)確描述害蟲發(fā)生量的變化趨勢，以提高害蟲發(fā)生量的預(yù)測精度為目標(biāo)，筆者提出一種基于混沌理論的害蟲發(fā)生量非線性預(yù)測模型(PSR–LSSVM)。首先采用混沌理論的相空間重構(gòu)(phase space reconstruction，PSR)[13]，將一維的害蟲發(fā)生量時間序列轉(zhuǎn)化成矩陣形式，挖掘隱藏于害蟲發(fā)生量歷史數(shù)據(jù)中的信息，再采用任意非線性映射和學(xué)習(xí)逼近能力的LSSVM對構(gòu)造的時間序列進(jìn)行訓(xùn)練，建立害蟲發(fā)生量預(yù)測模型，揭示害蟲復(fù)雜的發(fā)生過程，以期為害蟲發(fā)生量研究提供一種新的簡單、有效的方法。

1 害蟲發(fā)生量非線性預(yù)測模型

1.1 相空間重構(gòu)

為了研究時間序列數(shù)據(jù)的混沌性，Packard[14]提出相空間重構(gòu)理論，可以從實測時間序列中的某一分量了解非線性動力系統(tǒng)相空間的幾何特性，并在高維相空間中恢復(fù)混沌吸引子。Takens 等[15]證明了重構(gòu)系統(tǒng)與原始系統(tǒng)在系統(tǒng)特征上具有等價關(guān)系。對于一維混沌時間序列x(i)，i=1,2,…,n 進(jìn)行相空間重構(gòu)，重構(gòu)后的m 維狀態(tài)向量可表示為：

式中：N=n–(m–1)τ 為相點的個數(shù)；m 為嵌入維數(shù)；τ 為延遲時間。

1.2 最小二乘支持向量機(jī)

對于訓(xùn)練集{(xi,yi)}，i=1,2,…,n，LSSVM 通過非線性映射函數(shù)Φ(·)將其映射到高維特征空間進(jìn)行線性回歸，則有：

式中：ω 為權(quán)值向量；b 為偏置量。

綜合考慮預(yù)測性能和建模效率，式(2)可以變?yōu)榧s束條件的優(yōu)化問題：

式中：γ 為正則化參數(shù)；ei為預(yù)測誤差。

引入拉格朗日乘子[16]，將式(3)變?yōu)閷ε紗栴}，即：

式中：αi為拉格朗日乘子。

選擇RBF 核函數(shù)作為LSSVM 核函數(shù)，定義為：

最后，LSSVM 回歸模型為：

從LSSVM 建模過程可知，LSSVM 預(yù)測性能與其參數(shù)相關(guān)，因此采用粒子群優(yōu)化算法對LSSSVM參數(shù)進(jìn)行選擇。

1.3 害蟲發(fā)生量的非線性預(yù)測模型

設(shè)害蟲發(fā)生量的歷史數(shù)據(jù)時間序列為{x1,x2,…,xn}，xn是預(yù)測目標(biāo)值，建立害蟲發(fā)生量輸入x={xn–1,xn–2,…,xn–m}與輸出y={xn}之間的非線性映射關(guān)系：Rm→R，m 為映射維數(shù)。通過m 對害蟲發(fā)生量時間序列{x1,x2,…,xn}進(jìn)行樣本重構(gòu)，得到LSSVM 的學(xué)習(xí)樣本為：

則害蟲發(fā)生量時間序列的LSSVM 回歸函數(shù)為：

式中：t=m+1,m+2,…,n。

基于LSSVM 的害蟲發(fā)生量時間序列預(yù)測模型為：

2 應(yīng)用實例

2.1 數(shù)據(jù)源

為避免單個數(shù)據(jù)集的偶然性，選取云南省普洱市思茅區(qū)1985—2006年松毛蟲發(fā)生面積(表1)和浙江省仙居縣1983—1989年松毛蟲發(fā)生面積(表2)作為研究對象。浙江省仙居縣松毛蟲發(fā)生面積每年采集3次數(shù)據(jù)，收集時間分別為4月(出蟄)、7月(第一代)和 9月(第二代)。分別取 1983—2003、1983—1988年數(shù)據(jù)作為訓(xùn)練集，各數(shù)據(jù)集最后3個數(shù)據(jù)作為測試集。

表1 1985—2006年普洱市思茅區(qū)松毛蟲發(fā)生面積 Table 1 Occurrence area of Dendrolimus punctatus in Simao, Puer in 1985–2006 hm2

表2 1983—1989年浙江省仙居縣松毛蟲發(fā)生面積 Table 2 Occurrence area of Dendrolimus punctatus in Xianju county, Zhejiang in 1983–1989 hm2

2.2 歸一化處理

由于LSSVM 的核函數(shù)值通常依賴特征向量的內(nèi)積，而較大的屬性值影響訓(xùn)練效率，需在建模之前，對數(shù)據(jù)進(jìn)行歸一化處理，將其值縮放到[0,1]之間，具體為：

式中：ix 是原始數(shù)據(jù)；ix'為歸一化后的數(shù)據(jù)；xmax和 xmin分別代表數(shù)據(jù)的最大值和最小值。

2.3 延遲時間和嵌入維的確定

由于害蟲發(fā)生是一種弱混沌時間序列，因此松毛蟲發(fā)生面積延遲時間τ=1，采用G–P 法[17]選擇嵌入維數(shù)m。

1) 設(shè)延遲時間的初值為τ=1，嵌入維數(shù)的初值為m=1。

2) 選擇合適的臨界距離r，根據(jù)式(11)計算Cn(r)，采用2個向量最大分量差作為向量距離。

式中：M 為相點的個數(shù)；r 為臨界距離；θ 為Heaviside 單位函數(shù)。

3) 用最小二乘法擬合logC(r)n–logr 曲線中的直線段。

4) 增加嵌入維數(shù)，返回步驟2)。

2 種時間序列在不同嵌入維數(shù)下的關(guān)聯(lián)維數(shù)，如圖1 所示。從圖1 可知，當(dāng)嵌入維數(shù)m=2 時，普洱市思茅區(qū)的松毛蟲發(fā)生面積關(guān)聯(lián)維數(shù)達(dá)到飽和狀態(tài)，當(dāng)嵌入維數(shù)m=5 時，浙江省仙居縣的松毛蟲發(fā)生面積關(guān)聯(lián)維數(shù)達(dá)到基本飽和狀態(tài)，表明2 種時間序列的最佳嵌入維數(shù)分別為m=2 和m=5。

圖1 不同嵌入維數(shù)下的關(guān)聯(lián)維數(shù) Fig.1 Correlation dimension under different embedding dimensions

2.4 松毛蟲發(fā)生量時間序列重構(gòu)

采用上述延遲時間τ 和嵌入維數(shù)m 對2 種松毛蟲發(fā)生面積時間序列進(jìn)行重構(gòu)，并進(jìn)行歸一化處理，結(jié)果見表3 和表4。

表3 歸一化后重構(gòu)的訓(xùn)練集(普洱市思茅區(qū)) Table 3 Training set after normalize and reconstruction (Simao, Puer)

表4 歸一化后重構(gòu)的訓(xùn)練集(浙江省仙居縣) Table 4 Training set after normalize and reconstruction (Xianju county，Zhejiang)

3 結(jié)果與分析

3.1 PSR–LSSVM 的預(yù)測結(jié)果

將普洱市思茅區(qū)松毛蟲發(fā)生面積重構(gòu)后的訓(xùn)練樣本輸入LSSVM 學(xué)習(xí)，粒子群算法找到模型參數(shù)為：γ=100，σ=0.675。建立普洱市思茅區(qū)松毛蟲發(fā)生面積預(yù)測模型，對測試集進(jìn)行預(yù)測，同樣地，將浙江省仙居縣松毛蟲發(fā)生面積重構(gòu)后的訓(xùn)練樣本輸入LSSVM 學(xué)習(xí)，粒子群算法找到模型參數(shù)為：γ=1 000，σ=1.65。建立浙江省仙居縣松毛蟲發(fā)生面積預(yù)測模型，對測試集進(jìn)行預(yù)測，結(jié)果列于表5。由表5 可知，對普洱市思茅區(qū)松毛蟲發(fā)生面積的預(yù)測平均相對誤差為0.9%，對浙江省仙居縣松毛蟲發(fā)生面積的預(yù)測平均相對誤差為2.4%。

表5 PSR–LSSVM 的松毛蟲發(fā)生面積預(yù)測結(jié)果 Table 5 Prediction results of Dendrolimus punctatus occurrence areas based on PSR–LSSVM

綜上所述，對混沌理論與非線性預(yù)測算法LSSVM 進(jìn)行融合建立的預(yù)測模型，可以對松毛蟲發(fā)生面積進(jìn)行準(zhǔn)確預(yù)測，主要是由于一維松毛蟲發(fā)生面積時間序列中包含著吸引子的結(jié)構(gòu)，而重構(gòu)相空間能夠估計出系統(tǒng)演化的信息，使預(yù)測結(jié)果更符合松毛蟲發(fā)生面積的實際情況，同時利用非線性預(yù)測能力強(qiáng)的LSSVM 對重構(gòu)后的數(shù)據(jù)進(jìn)行建模，并采用粒子群優(yōu)化算法優(yōu)化LSSVM 參數(shù)，進(jìn)一步提高了松毛蟲發(fā)生面積的精度，同時該模型要求的樣本不大，克服了一般混沌預(yù)測模型要求大樣本的缺陷，通用性較好，適合于小樣本的害蟲預(yù)測。

3.2 與其他模型預(yù)測結(jié)果對比

為了使PSR–LSSVM 的預(yù)測結(jié)果更具說服力，采用線性預(yù)測模型(ARIMA)、PSR–BPNN 神經(jīng)網(wǎng)絡(luò)模型(PSR–BPNN)進(jìn)行對比試驗。模型性能評價指標(biāo)為平均絕對百分誤差(MAPE)和均方根誤差(RMSE)。表6 的結(jié)果表明：

1) BP 神經(jīng)網(wǎng)絡(luò)和LSSVM 預(yù)測比線性模型ARIMA 預(yù)測精度有大幅度提高。由于松毛蟲發(fā)生具有非線性特征，BP 神經(jīng)網(wǎng)絡(luò)、LSSVM 能夠更好地反映松毛蟲發(fā)生量的變化趨勢，有效提高了松毛蟲發(fā)生量的預(yù)測精度。

2) 相對于PSR–BP 神經(jīng)網(wǎng)絡(luò)(PSR–BPNN)，LSSVM 的預(yù)測精度更高，說明LSSVM 能較好地克服BP 神經(jīng)網(wǎng)絡(luò)存在的過擬合、易陷入局部最優(yōu)等缺陷，對于小樣本害蟲時間序列預(yù)測，優(yōu)勢更加明顯。

3) 將混沌理論和LSSVM 用于害蟲發(fā)生量預(yù)測，誤差較小，結(jié)果接近真實值，說明PSR–LSSVM能夠更加精確反映害蟲發(fā)生變化趨勢和變化特點，預(yù)測結(jié)果更加可靠。

表6 各模型的預(yù)測誤差對比 Table 6 Comparison for prediction error of different models

[1] 張孝羲．昆蟲生態(tài)及預(yù)測預(yù)報[M]．3 版．北京：中國農(nóng)業(yè)出版社，2002．

[2] 許曉風(fēng)，馬飛，丁宗澤，等．褐飛虱發(fā)生的相空間線性回歸預(yù)測模型[J]．昆蟲學(xué)報，2002，45(4)：548–551．

[3] 岑冠軍，黃壽山，肖莉，等．ARIMA 模型在小菜蛾幼蟲種群動態(tài)中的應(yīng)用[J]．華南農(nóng)業(yè)大學(xué)學(xué)報，2008，29(1)：109–114．

[4] 張真，李典謨，查光濟(jì)．馬尾松毛蟲種群動態(tài)的時間序列分析及復(fù)雜性動態(tài)研究[J]．生態(tài)學(xué)報，2002，22(7)：1061–1067．

[5] 陳順立，張華峰，張潮巨，等．神經(jīng)網(wǎng)絡(luò)在松墨天牛發(fā)生量預(yù)報中的應(yīng)用[J]．福建林學(xué)院學(xué)報，2006，26(2)：6–9．

[6] 石晶晶，劉占宇，張莉麗，等．基于支持向量機(jī)的稻縱卷葉螟危害水稻高光譜遙感識別[J]．中國水稻科學(xué)，2009，23(3)：331–334．

[7] 向昌盛，周子英，張林峰．支持向量機(jī)在害蟲發(fā)生量預(yù)測中的應(yīng)用[J]．生物信息學(xué)，2011，9(1)：28–31．

[8] 向昌盛，周子英，武麗娜．基于 ARIMA 和 DSVM組合模型的松毛蟲發(fā)生面積預(yù)測[J]．湖南農(nóng)業(yè)大學(xué)學(xué)報：自然科學(xué)版，2010，36(4)：430–434．

[9] 譚泗橋，林雪梅，陳淵，等．基于地統(tǒng)計學(xué)的多維時間序列模型及其在生態(tài)學(xué)中的應(yīng)用[J]．湖南農(nóng)業(yè)大學(xué)學(xué)報：自然科學(xué)版，2009，35(4)：433–436．

[10] 向昌盛，袁哲明．最小二乘支持向量機(jī)在害蟲預(yù)測中的應(yīng)用[J]．湖南科技大學(xué)學(xué)報：自然科學(xué)版，2012，27(2)：71–74．

[11] 向昌盛，周子英．ARIMA 與SVM組合模型在害蟲預(yù)測中的應(yīng)用[J]．昆蟲學(xué)報，2010，53(9)：1055–1060．

[12] 朱軍生，翟保平，劉英智．基于小波分解的害蟲發(fā)生非平穩(wěn)時間序列分析和預(yù)測[J]．南京農(nóng)業(yè)大學(xué)學(xué)報，2011，34(3)：61–66．

[13] 馬飛，許曉風(fēng)，張夕林，等．相空間重構(gòu)與神經(jīng)網(wǎng)絡(luò)融合預(yù)測模型及其在害蟲測報中的應(yīng)用[J]．2002，22(8)：1297–1301．

[14] Packard N H，Crutchifield J P，F(xiàn)armer J D，et al. Geometry from a time series[J]．Phys Rev Lett，1980，45(6)：712–716．

[15] Takens F．Determing strange attractors in turbulence [J]．Lecture Notes in Mathematics，1988，898：361–381．

[16] An S J，Liu W Q，Venkatesh S．Fast cross validation algorithms for least squares support vector machines and kernel ridge regression[J]．Pattern Recognition，2007，40(2)：2154–2162．

[17] Grassberger P，Privacies I．Characterization of strange attractors[J]．Physical Review Letters，1983，50(5)：346 – 349．