基于响应面方法的马尾松毛虫发生量混沌特性检测及其预测

陈绘画 王坚娅 徐志宏

(浙江省仙居县林业局，仙居，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.