基于響應面方法的馬尾松毛蟲發生量混沌特性檢測及其預測

陳繪畫 王堅婭 徐志宏

(浙江省仙居縣林業局，仙居，317300) (浙江農林大學)

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

1 材料與方法

1.1 材料來源

馬尾松毛蟲資料來自浙江省仙居縣森林病蟲防治站。該縣于1982年發現有馬尾松毛蟲危害，自1983年起，每年都在相對固定的時間內調查各代蟲情的危害情況。馬尾松毛蟲在該縣一年發生2代，以幼蟲于當年9—10月越冬，翌年3—4月再出蟄危害，故將第2代舍去，按一年2次調查結果進行處理，共計56代調查數據。將1983—2008年的調查數據用于建立響應面模型，2009—2010年的調查數據用作所建立模型預報結果的檢驗。

在馬尾松毛蟲常發區和偶發區，于各鄉(鎮、街道)行政區域內，以林班為單位，在每代馬尾松毛蟲幼蟲期對發生情況按林班線、林間小道等巡視路線進行詳細線路踏查，發現有蟲情或災情，設立臨時標準地進行調查。每塊標準地面積大于0.2 hm2。在標準地內采用對角線或平行線抽樣方法，隨機選取20株標準樹，詳細調查蟲口密度、有蟲株率，目測有蟲面積;否則，按不同類型的林相、立地條件，在調查線路上隨機選取20株標準樹，詳細調查蟲口密度、有蟲株率，目測有蟲面積。

1.2 時間序列的自相關與偏自相關分析

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

1.3 響應面模型

在昆蟲生態領域，復雜的動力學可能來自簡單的種群模型［6-7］，離散世代差分方程是描述種群動態最簡單的模型，它的1級滯后的一般形式為：

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

由函數關系F的不同選擇產生了Logistic模型、Hassell模型等幾種重要的模型［17］。(1)式再增加1級時滯，則變為：

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

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

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

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

2 結果與分析

2.1 有蟲面積序列的內在動態定性診斷

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

圖1 馬尾松毛蟲有蟲面積時間序列自相關圖

圖2 馬尾松毛蟲有蟲面積時間序列偏自相關圖

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

2.2 有蟲面積響應面模型的建立及序列混沌特性的檢測

馬尾松毛蟲有蟲面積序列的ACF接近可衰減為0的阻尼正弦型振蕩，說明有蟲面積的周期模式為平穩周期，可用響應面方法重構有蟲面積的內在動態，以確定有蟲面積動態所屬確定性動態的種類［17］。

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

方程經 F 檢驗，(F=43.386) ＞ (F0.01(2，49)=5.08)，復相關系數R=0.799 4，說明模型(6)的擬合程度較高，可以由其預測馬尾松毛蟲的有蟲面積。

建立模型(6)后，就可通過迭代來確定有蟲面積動態吸引子的類型。初值N1、N2設為有蟲面積時間序列的平均值，用 Nt、Nt-1值作圖來顯示模擬的軌跡［13，17］，見圖3。圖3 的Nt、Nt-1軌跡象一個拉開又折疊的橢圓，根據文獻［17］中提供的混沌判定法則，可以確定有蟲面積序列是混沌序列。

圖3 馬尾松毛蟲有蟲面積時間序列的Nt、Nt-1軌跡

2.3 有蟲面積混沌序列的預測及檢驗

利用建立的馬尾松毛蟲有蟲面積離散世代非線性模型，預測沒有參與建模的2009—2010年4代有蟲面積(表1)。由表1得2009—2010年4代有蟲面積的平均相對誤差為10.11%，以松毛蟲災害長期(相隔2個世代或1 a以上)預報準確率60%以上為標準［19］，4代的預測成功率達100%。

表1 2009—2010年4代有蟲面積的預測值與實際值

3 結論與討論

基于非線性動力學原理，筆者在構建馬尾松毛蟲有蟲面積2級時滯離散世代非線性差分響應面模型的基礎上，判定馬尾松毛蟲有蟲面積序列為混沌序列;利用混沌序列內在的確定性和非線性，對2009—2010年4代有蟲面積實測數據資料進行預報，得到了4次預報的平均相對誤差僅為10.11%，預報成功率達100%。這為利用時間序列的可預測性構建新型預測模型提供可靠的依據。

通常認為馬尾松毛蟲具有周期性暴發的特點。通過自相關和偏自相關分析，認為馬尾松毛蟲有蟲面積動態是平穩的，周期性不顯著，具有一定的復雜性。馬尾松毛蟲各代間有蟲面積存在時間延遲效應，由偏自相關函數分析的結果可見，當代有蟲面積大小與前1代、前2代的有蟲面積大小有關，屬于2級密度相關，最佳響應面模型參數的估計結果嵌入維為2。

馬尾松毛蟲多發生在低海拔、人類活動頻繁的松林中，受到人為的干預較大。就系統的角度而言，由于林業生產受當地經濟社會等的調節作用而導致馬尾松林面積的不斷變化，這樣對于馬尾松毛蟲種群長期行為的預測是很難做到的。從資料的搜集和系統的監測角度看，本研究所采用的統計資料長度，能滿足分析的要求，但松林生態系統的不穩定性決定了馬尾松毛蟲種群動態的預測只能在短時間范圍內進行;長期的資料則由于時代的發展而采用的標準不同，難以準確反映一個地區的發生情況和長期規律性。基于這樣的資料，無論采用哪種分析方法都會有偏差。

［1］蔡邦華.關于防治松毛蟲的研究工作［J］.科學通報，1955(4)：43-45.

［2］陳曉峰，李典謨.馬尾松毛蟲種群動態及模擬研究［J］.昆蟲學報，1993，36(1)：56-62.

［3］馬小明.馬尾松毛蟲種群動態模擬研究［J］.林業科學，1994，30(1)：88-92.

［4］張真，李典謨.馬尾松毛蟲暴發機制分析［J］.林業科學，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］.生態學報，2002，22(7)：1061-1067.

［17］趙中華，沈佐銳.昆蟲種群動態非線性建模理論與應用［J］.生物數學學報，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］龐正轟.馬尾松毛蟲災害預測預報綜合技術研究［D］.北京：北京林業大學，2004.