陳繪畫(huà) 王堅(jiān)婭 徐志宏

(浙江省仙居縣林業(yè)局，仙居，317300) (浙江農(nóng)林大學(xué))

馬尾松毛蟲(chóng)(Dendrolimuspunctatus Walker)是我國(guó)南方馬尾松林(Pinus massoniana)的歷史性大害蟲(chóng)，從20世紀(jì)50年代開(kāi)始我國(guó)就有學(xué)者對(duì)馬尾松毛蟲(chóng)進(jìn)行了研究［1］，陳曉峰、馬小明、張真等則先后對(duì)馬尾松毛蟲(chóng)的種群動(dòng)態(tài)及其暴發(fā)機(jī)制進(jìn)行了研究［2-4］。美國(guó)數(shù)學(xué)生態(tài)學(xué)家May R等發(fā)表了一系列論文［5-7］，揭示了生命系統(tǒng)中混沌存在的可能性。此后，許多生態(tài)學(xué)家在研究如何檢測(cè)自然種群的混沌行為［8-15］。張真等用Lyapunov指數(shù)檢測(cè)了安徽冬至縣金寺山林場(chǎng)和福建連江縣馬尾松毛蟲(chóng)發(fā)生面積序列的混沌現(xiàn)象［16］。一般情況下，計(jì)算Lyapunov指數(shù)需要較多的數(shù)據(jù)量，而昆蟲(chóng)限于自身的生物學(xué)特性，其種群序列通常只有幾十個(gè)數(shù)據(jù)，Ellner等［14］提出3種適于這種少數(shù)據(jù)量、具噪音的低維動(dòng)態(tài)系統(tǒng)模型，即：響應(yīng)面方法(RSM)、薄板樣條法(TPS)和前饋神經(jīng)網(wǎng)絡(luò)法(FNN)，其中RSM適合于20～50數(shù)據(jù)點(diǎn)的時(shí)間序列;FNN適合于50～500數(shù)據(jù)點(diǎn)的序列。本研究利用仙居縣1983—2010年的馬尾松毛蟲(chóng)調(diào)查資料，用時(shí)間序列和響應(yīng)面法分析馬尾松毛蟲(chóng)有蟲(chóng)面積動(dòng)態(tài)變化過(guò)程中的混沌現(xiàn)象，揭示馬尾松毛蟲(chóng)種群動(dòng)態(tài)的復(fù)雜性。

1 材料與方法

1.1 材料來(lái)源

馬尾松毛蟲(chóng)資料來(lái)自浙江省仙居縣森林病蟲(chóng)防治站。該縣于1982年發(fā)現(xiàn)有馬尾松毛蟲(chóng)危害，自1983年起，每年都在相對(duì)固定的時(shí)間內(nèi)調(diào)查各代蟲(chóng)情的危害情況。馬尾松毛蟲(chóng)在該縣一年發(fā)生2代，以幼蟲(chóng)于當(dāng)年9—10月越冬，翌年3—4月再出蟄危害，故將第2代舍去，按一年2次調(diào)查結(jié)果進(jìn)行處理，共計(jì)56代調(diào)查數(shù)據(jù)。將1983—2008年的調(diào)查數(shù)據(jù)用于建立響應(yīng)面模型，2009—2010年的調(diào)查數(shù)據(jù)用作所建立模型預(yù)報(bào)結(jié)果的檢驗(yàn)。

在馬尾松毛蟲(chóng)常發(fā)區(qū)和偶發(fā)區(qū)，于各鄉(xiāng)(鎮(zhèn)、街道)行政區(qū)域內(nèi)，以林班為單位，在每代馬尾松毛蟲(chóng)幼蟲(chóng)期對(duì)發(fā)生情況按林班線、林間小道等巡視路線進(jìn)行詳細(xì)線路踏查，發(fā)現(xiàn)有蟲(chóng)情或?yàn)?zāi)情，設(shè)立臨時(shí)標(biāo)準(zhǔn)地進(jìn)行調(diào)查。每塊標(biāo)準(zhǔn)地面積大于0.2 hm2。在標(biāo)準(zhǔn)地內(nèi)采用對(duì)角線或平行線抽樣方法，隨機(jī)選取20株標(biāo)準(zhǔn)樹(shù)，詳細(xì)調(diào)查蟲(chóng)口密度、有蟲(chóng)株率，目測(cè)有蟲(chóng)面積;否則，按不同類型的林相、立地條件，在調(diào)查線路上隨機(jī)選取20株標(biāo)準(zhǔn)樹(shù)，詳細(xì)調(diào)查蟲(chóng)口密度、有蟲(chóng)株率，目測(cè)有蟲(chóng)面積。

1.2 時(shí)間序列的自相關(guān)與偏自相關(guān)分析

時(shí)間序列的自相關(guān)函數(shù)(ACF)主要用于對(duì)昆蟲(chóng)種群動(dòng)態(tài)的定性分析。令馬尾松毛蟲(chóng)各代調(diào)查的有蟲(chóng)面積為Nt，Lt≡log Nt，然后計(jì)算Lt與Lt-τ的相關(guān)系數(shù)得到自相關(guān)函數(shù)，其中時(shí)滯τ(τ=1，2，…)為2次調(diào)查間的時(shí)間間隔。再以τ為函數(shù)對(duì)自相關(guān)系數(shù)作圖，該圖形形狀能顯示種群動(dòng)態(tài)是平穩(wěn)的還是周期的信息。偏自相關(guān)函數(shù)(PACF)則是在排除其他數(shù)據(jù)點(diǎn)影響下各時(shí)滯數(shù)據(jù)點(diǎn)的相關(guān)系數(shù)，PACF能夠確定模型的滯后代數(shù)［16-18］。

1.3 響應(yīng)面模型

在昆蟲(chóng)生態(tài)領(lǐng)域，復(fù)雜的動(dòng)力學(xué)可能來(lái)自簡(jiǎn)單的種群模型［6-7］，離散世代差分方程是描述種群動(dòng)態(tài)最簡(jiǎn)單的模型，它的1級(jí)滯后的一般形式為：

式中：Nt為t世代時(shí)的種群密度;t為種群世代數(shù);εt為引起種群變化的外部因子;F為描述上代種群密度影響下一代種群密度的函數(shù)。

由函數(shù)關(guān)系F的不同選擇產(chǎn)生了Logistic模型、Hassell模型等幾種重要的模型［17］。(1)式再增加1級(jí)時(shí)滯，則變?yōu)椋?/p>

令種群世代轉(zhuǎn)換率rt≡log(Nt/Nt-1)，種群密度的對(duì)數(shù)轉(zhuǎn)換為：式中：θ1，θ2為冪轉(zhuǎn)換系數(shù)。由此得到8參數(shù)的非線性響應(yīng)面模型為：

式中：a0，a1，a2，a11，a22，a12為確定模型(4)的參數(shù)，其中 a0表示常數(shù)項(xiàng)待定系數(shù)，ai表示一次項(xiàng)待定系數(shù)，aii表示二次項(xiàng)待定系數(shù)，aij表示交叉項(xiàng)待定系數(shù)，i，j=1，2，3，…，n;εt為隨機(jī)誤差。

式(4)的另一種表現(xiàn)形式為：

式中：k1、k2和 k3是 a0，a1，a2，a11，a22，a12的函數(shù)［15］;d 為模型的階。

2 結(jié)果與分析

2.1 有蟲(chóng)面積序列的內(nèi)在動(dòng)態(tài)定性診斷

對(duì)1983—2008年的馬尾松毛蟲(chóng)有蟲(chóng)面積序列進(jìn)行自相關(guān)和偏自相關(guān)分析(圖1、圖2)，其ACF接近可衰減為0的阻尼正弦型振蕩，在滯后6代時(shí)其自相關(guān)系數(shù)還達(dá)到95%的顯著水平，但這并不表明3～4 a前的有蟲(chóng)面積還能直接影響到當(dāng)前的有蟲(chóng)面積，更大的可能則是通過(guò)中間連年的影響才使得3～4 a前的有蟲(chóng)面積與現(xiàn)在的有蟲(chóng)面積的自相關(guān)系數(shù)達(dá)到95%的顯著水平，因而用ACF無(wú)法得到模型的滯后代數(shù)。

圖1 馬尾松毛蟲(chóng)有蟲(chóng)面積時(shí)間序列自相關(guān)圖

圖2 馬尾松毛蟲(chóng)有蟲(chóng)面積時(shí)間序列偏自相關(guān)圖

使用PACF可以解決模型的滯后代數(shù)難以確定的問(wèn)題。從圖2可以看出，PACF在滯后l代及2代時(shí)均達(dá)到95%的顯著水平。由Bartlett的2法則，得：(|RPACF［2］|=0.328 2)＞(2/=0.280 1)(RPACF［2］為滯后2代的偏自相關(guān)系數(shù))。因此，可以確定離散世代差分方程模型的滯后代數(shù)為2［18］。

2.2 有蟲(chóng)面積響應(yīng)面模型的建立及序列混沌特性的檢測(cè)

馬尾松毛蟲(chóng)有蟲(chóng)面積序列的ACF接近可衰減為0的阻尼正弦型振蕩，說(shuō)明有蟲(chóng)面積的周期模式為平穩(wěn)周期，可用響應(yīng)面方法重構(gòu)有蟲(chóng)面積的內(nèi)在動(dòng)態(tài)，以確定有蟲(chóng)面積動(dòng)態(tài)所屬確定性動(dòng)態(tài)的種類［17］。

擬合如方程(4)那樣的模型，一般的做法是采取線性回歸形式，即采用最小二乘法對(duì) θ1、θ2的所有組合{-1.0，-0.5，0，…，2.5，3.0}作方程(4)的回歸擬合，選取使剩余平方和最小的各參數(shù)值作為對(duì)方程(4)的估計(jì)［13，17］。然而，上述擬合方法有可能遺漏θ值小于-1或大于3的組合［15］。為避免遺漏θ值組合，借助MATLAB軟件，筆者采用非線性回歸方法對(duì)方程(5)進(jìn)行直接擬合，模型階數(shù)取3、滯后2代，擬合的1983—2008年有蟲(chóng)面積的非線性模型為：

方程經(jīng) F 檢驗(yàn)，(F=43.386) ＞ (F0.01(2，49)=5.08)，復(fù)相關(guān)系數(shù)R=0.799 4，說(shuō)明模型(6)的擬合程度較高，可以由其預(yù)測(cè)馬尾松毛蟲(chóng)的有蟲(chóng)面積。

建立模型(6)后，就可通過(guò)迭代來(lái)確定有蟲(chóng)面積動(dòng)態(tài)吸引子的類型。初值N1、N2設(shè)為有蟲(chóng)面積時(shí)間序列的平均值，用 Nt、Nt-1值作圖來(lái)顯示模擬的軌跡［13，17］，見(jiàn)圖3。圖3 的Nt、Nt-1軌跡象一個(gè)拉開(kāi)又折疊的橢圓，根據(jù)文獻(xiàn)［17］中提供的混沌判定法則，可以確定有蟲(chóng)面積序列是混沌序列。

圖3 馬尾松毛蟲(chóng)有蟲(chóng)面積時(shí)間序列的Nt、Nt-1軌跡

2.3 有蟲(chóng)面積混沌序列的預(yù)測(cè)及檢驗(yàn)

利用建立的馬尾松毛蟲(chóng)有蟲(chóng)面積離散世代非線性模型，預(yù)測(cè)沒(méi)有參與建模的2009—2010年4代有蟲(chóng)面積(表1)。由表1得2009—2010年4代有蟲(chóng)面積的平均相對(duì)誤差為10.11%，以松毛蟲(chóng)災(zāi)害長(zhǎng)期(相隔2個(gè)世代或1 a以上)預(yù)報(bào)準(zhǔn)確率60%以上為標(biāo)準(zhǔn)［19］，4代的預(yù)測(cè)成功率達(dá)100%。

表1 2009—2010年4代有蟲(chóng)面積的預(yù)測(cè)值與實(shí)際值

3 結(jié)論與討論

基于非線性動(dòng)力學(xué)原理，筆者在構(gòu)建馬尾松毛蟲(chóng)有蟲(chóng)面積2級(jí)時(shí)滯離散世代非線性差分響應(yīng)面模型的基礎(chǔ)上，判定馬尾松毛蟲(chóng)有蟲(chóng)面積序列為混沌序列;利用混沌序列內(nèi)在的確定性和非線性，對(duì)2009—2010年4代有蟲(chóng)面積實(shí)測(cè)數(shù)據(jù)資料進(jìn)行預(yù)報(bào)，得到了4次預(yù)報(bào)的平均相對(duì)誤差僅為10.11%，預(yù)報(bào)成功率達(dá)100%。這為利用時(shí)間序列的可預(yù)測(cè)性構(gòu)建新型預(yù)測(cè)模型提供可靠的依據(jù)。

通常認(rèn)為馬尾松毛蟲(chóng)具有周期性暴發(fā)的特點(diǎn)。通過(guò)自相關(guān)和偏自相關(guān)分析，認(rèn)為馬尾松毛蟲(chóng)有蟲(chóng)面積動(dòng)態(tài)是平穩(wěn)的，周期性不顯著，具有一定的復(fù)雜性。馬尾松毛蟲(chóng)各代間有蟲(chóng)面積存在時(shí)間延遲效應(yīng)，由偏自相關(guān)函數(shù)分析的結(jié)果可見(jiàn)，當(dāng)代有蟲(chóng)面積大小與前1代、前2代的有蟲(chóng)面積大小有關(guān)，屬于2級(jí)密度相關(guān)，最佳響應(yīng)面模型參數(shù)的估計(jì)結(jié)果嵌入維為2。

馬尾松毛蟲(chóng)多發(fā)生在低海拔、人類活動(dòng)頻繁的松林中，受到人為的干預(yù)較大。就系統(tǒng)的角度而言，由于林業(yè)生產(chǎn)受當(dāng)?shù)亟?jīng)濟(jì)社會(huì)等的調(diào)節(jié)作用而導(dǎo)致馬尾松林面積的不斷變化，這樣對(duì)于馬尾松毛蟲(chóng)種群長(zhǎng)期行為的預(yù)測(cè)是很難做到的。從資料的搜集和系統(tǒng)的監(jiān)測(cè)角度看，本研究所采用的統(tǒng)計(jì)資料長(zhǎng)度，能滿足分析的要求，但松林生態(tài)系統(tǒng)的不穩(wěn)定性決定了馬尾松毛蟲(chóng)種群動(dòng)態(tài)的預(yù)測(cè)只能在短時(shí)間范圍內(nèi)進(jìn)行;長(zhǎng)期的資料則由于時(shí)代的發(fā)展而采用的標(biāo)準(zhǔn)不同，難以準(zhǔn)確反映一個(gè)地區(qū)的發(fā)生情況和長(zhǎng)期規(guī)律性?；谶@樣的資料，無(wú)論采用哪種分析方法都會(huì)有偏差。

［1］蔡邦華.關(guān)于防治松毛蟲(chóng)的研究工作［J］.科學(xué)通報(bào)，1955(4)：43-45.

［2］陳曉峰，李典謨.馬尾松毛蟲(chóng)種群動(dòng)態(tài)及模擬研究［J］.昆蟲(chóng)學(xué)報(bào)，1993，36(1)：56-62.

［3］馬小明.馬尾松毛蟲(chóng)種群動(dòng)態(tài)模擬研究［J］.林業(yè)科學(xué)，1994，30(1)：88-92.

［4］張真，李典謨.馬尾松毛蟲(chóng)暴發(fā)機(jī)制分析［J］.林業(yè)科學(xué)，2008，44(1)：140-150.

［5］May R M.Biological populations with nonoverlapping generations：stable points，stable cycles，and chaos［J］.Science，1974，186：645-647.

［6］May R M.Simple mathematical models with very complicated dynamics［J］.Nature，1976，261：459-467.

［7］May R M，Oster G F.Bifurcations and dynamic complexity in simple ecological models［J］.The American Naturalist，1976，110：573-599.

［8］Hassell M P，Lawtto JH，May R M.Patterns of dynamical behaviour in single-species populations［J］.Journal of Animal Ecology，1976，45(2)：471-486.

［9］Schaffer WM.Order and chaos in ecological systems［J］.Ecology，1985，66(1)：93-106.

［10］Schaffer W M，Kot M.Chaos in ecological systems：the coals that Newcastle forgot［J］.Trends in Ecology and Evolution，1986，1(3)：58-63.

［11］Berryman A A，Millstein J A.Are ecological systems chaotic-and if not，why not? ［J］.Trends in Ecology and Evolution，1989，4(1)：26-28.

［12］Bellows T SJr.The descriptive properties of some models for density dependence［J］.Journal of Animal Ecology，1981，50：139-156.

［13］Turchin P，Taylor A D.Complex dynamics in ecological time series［J］.Ecology，1992，73(1)：289-305.

［14］Ellner S，Turchin P.Chaos in a noisy world：new methods and evidence from time-series analysis［J］.The American Naturalist，1995，145(3)：343-375.

［15］Perry JN，Woiwod I P，Hanski I.Using response-surface methodology to detect chaos in ecological time series［J］.Oikos，1993，68(2)：329-339.

［16］張真，李典謨，查光濟(jì).馬尾松毛蟲(chóng)種群動(dòng)態(tài)的時(shí)間序列分析及復(fù)雜性動(dòng)態(tài)研究［J］.生態(tài)學(xué)報(bào)，2002，22(7)：1061-1067.

［17］趙中華，沈佐銳.昆蟲(chóng)種群動(dòng)態(tài)非線性建模理論與應(yīng)用［J］.生物數(shù)學(xué)學(xué)報(bào)，2001，16(4)：439-447.

［18］Berryman A，Turchin P.Identifying the density-dependent structure underlying ecological time series［J］.Oikos，2001，92(2)：265-270.

［19］龐正轟.馬尾松毛蟲(chóng)災(zāi)害預(yù)測(cè)預(yù)報(bào)綜合技術(shù)研究［D］.北京：北京林業(yè)大學(xué)，2004.