灾变模型在马尾松毛虫幼虫发生量预报中的应用

张书平， 余 燕， 毕守东， 周夏芝， 邹运鼎， 张国庆， 张 桢， 方国飞， 宋玉双

（1. 安徽农业大学 理学院， 安徽 合肥230036； 2. 安徽农业大学 林学与园林学院， 安徽 合肥230036；3. 安徽省潜山县林业局， 安徽 潜山246300； 4. 国家林业和草原局 森林病虫害防治总站， 辽宁 沈阳110034）

马尾松毛虫Dendrolimus punctatus分布于安徽、 河南、 四川、 贵州、 陕西、 云南、 江西、 江苏、 湖南、 浙江、 福建、 广东、 海南和广西等， 主要危害马尾松Pinus massoniana， 还危害黑松Pinus thunbergii， 火炬松Pinus taeda， 湿地松Pinus elliottii， 晚松Pinus rigidavar.serotina和海南松Pinus fenzeliana等松属Pinus植物。 20 世纪中叶在中国森林害虫中马尾松毛虫是发生最广、 危害面积最大、 经常猖獗成灾的害虫。 在广大丘陵地区虫害此起彼伏， 为害时如同火烧， 造成了巨大的经济效益和生态效益损失。该虫不但影响林业生产， 还危害人身健康［1-4］， 人们的林事活动中接触马尾松毛虫毒毛， 容易引发皮炎和关节肿痛。 进入21 世纪， 封山育林、 混交、 间作等措施优化了森林生态环境， 使马尾松毛虫的危害得到有效控制， 但该虫具有强大的繁殖潜力， 遇到有利的生态环境极易爆发成灾， 不能放松对它的监测。 马尾松毛虫1 a 发生2～4 代， 发生世代的多少随不同地方而异。 在河南省信阳地区1 a 发生2 代为主， 在长江流域诸省1 a 发生2～3 代， 而在广东、 广西、 福建南部1 a 发生3～4 代， 海南1 a 发生4～5代［1］。 安徽潜山县1 a 发生3 代， 即4-6 月上旬为越冬代， 6 月上旬至8 月中下旬为1 代， 8 月中下旬至12 月为2 代。 马尾松毛虫发生的预测预报是对其进行综合防治的基本工作。 研究者［5-9］分别采用不同的预测方法预测马尾松毛虫的发生量、 虫害等级、 发生类别、 发生空间格局， 为马尾松毛虫的综合防治工作提供了有力支持。 由于各地气象条件、 植被条件和地形地貌等不同， 马尾松毛虫的发生特点也不完全相同。 为了有效地防治马尾松毛虫， 本研究采用灰色理论中的灾变预测法研究了安徽省潜山县马尾松毛虫幼虫越冬代、 1 代和2 代严重发生的年份， 以期为马尾松毛虫的综合治理提供科学依据。

1 材料与方法

1.1 材料来源

马尾松毛虫资料来自国家林业和草原局森林病虫中心测报点——安徽省潜山县森林病虫防治站， 时间跨度为1989-2016 年， 其中1998 年缺失。 采用踏查和详查相结合的方法， 沿林班线、 林道、 公路、铁路等线路调查， 目测发生范围、 危害状况， 发现虫情或灾情立即设临时标准地， 采取平行线抽样法抽取20 株标准株详查。 幼虫期调查， 1～2 龄幼虫调查枯黄卷曲的枝数， 推算幼虫数， 3 龄以上的幼虫且3 m 以下的小树直接调查合计树冠上的幼虫数， 大树用“虫粪粒推算法” 调查， 幼虫越冬期间调查树干基部树皮缝中的幼虫数推算全部虫口。

1.2 灾变预测的建模方法

灰色系统理论认为， 一切随机量都是在一定范围内、 一定时段上变化的灰色量及其灰色过程， 主张从事物内部研究其发展变化规律。 对于灰色量的处理， 不是去寻求它的统计规律和概率分布， 而是从无规律的原始数据中找出规律， 即对数据通过一定方式处理后再建立模型。 因为客观系统无论怎样复杂，它总是有关联、 有序、 有整体功能的， 作为系统行为特征的数据， 总是蕴含着某种规律。

如果在x（0）中有新的数列n′＜n是符合灾变定义的数列， 则对于上灾变其中：i′为灾变年序号，V′为灾变年序号取值范围。

为了便于辨认， 特记灾变年序号i′为z（i′）或z（i）。 并作灾变号映射记z的1 次累加生成数（accumulated generating operation，则灾变预测的GM（1,1）模型有如下算式：其中：t为年份，a为发展灰数，u为内生控制灰数，为待估灰数， B 为累加矩阵， BT为转置矩阵，yN为常数项向量。

2 结果与分析

2.1 马尾松毛虫累计虫口数

马尾松毛虫每年越冬代、 1 代幼虫和2 代幼虫的累计虫口数见表1， 将其作为原始数据。

表1 马尾松毛虫越冬代、 1 代幼虫和2 代幼虫累计虫口数Table 1 Accumulative population of D. punctatus during wintering, first and second generation larvae

2.2 马尾松毛虫越冬代灾变预测

2.2.1 确定灾变阈值 根据安徽省潜山县森林病虫害防治总站提供的1983-2016 年马尾松毛虫发生情况的资料， 马尾松毛虫越冬代虫口数大于15 头·株-1的年份为大发生年， 故以每株虫口数15 头为灾变阈值。

2.2.2 对原始数列作有关映射 原始数列为x（0）={x（0）（1）， x（0）（2）， …， x（0）（28）}={42.5， 37.3， 18.5，12.2， 7.8， 18.6， 18.2， 38.6， 42.1， 6.2， 7.6， 6.2， 7.1， 15.7， 7.6， 6.4， 6.6， 6.4， 6.6， 8.4， 16.2，16.6， 7.2， 7.2， 6.4， 6.2， 7.1}。 对原始数列作灾变映射， 按ξ≥15 得灾变数列：

然后再作灾变序号映射p，p{xξ（0）}={z}，pxξ（0）={z（1′），z（2′），z（3′），z（4′），z（5′），z（6′），z（7′），z（8′），z（9′），z（10′）}={1′， 2′， 3′， 4′， 5′， 6′， 7′， 8′， 9′， 10′}=z。

由于1′=1， 2′=2， 3′=3， 4′=6， 5′=7， 6′=8， 7′=9， 8′=15， 9′=22， 10′=23， 故得到灾变日期集z={1， 2， 3， 6， 7， 8， 9， 15， 22， 23}即为建模时所需的原始数列。

2.2.3 建立灰色预测GM（1,1）模型 对灾变年序集z 建立GM（1,1）模型， 按原始数列z= {1， 2， 3， 6，7， 8， 9， 15， 22， 23}作OAG（一次累加生成算子）得：

OAG{z（0）}=z（1）={z（1）（1）， z（1）（2）， z（1）（3）， z（1）（4）， z（1）（5）， z（1）（6）， z（1）（7）， z（1）（8）， z（1）（9）， z（1）（10）}={1， 3， 6， 12， 19， 27， 36， 51， 73， 96}。

按z（1）建模：由此得出马尾松毛虫越冬代虫口数的GM（1,1）灾变预测模型：

将根据马尾松毛虫越冬代虫口数的GM（1,1）灾变预测模型得到的拟合值、 误差等列于表2。

表2 马尾松毛虫越冬代虫口数预测的拟合值与误差Table 2 Fitting value and error of population prediction of D. punctatus in overwintering generation

对观察值和拟合值间的差异进行t检验，t=0.101 7，df=18 时，t0.05=2.10，t＜t0.05， 差异不显著， 平均误差（2 个均数之差）只有观察值均数的3.40%。 实际观察值与通过GM（1,1）灾变预测模型得到的拟合值之间误差较小， 各灾变年的平均精度为84.40%， 总体精度（观察值的平均值减去误差的平均值与观察值的平均值之比） 为96.25%。 由此可以通过这个模型求得未来5 个时刻的预测值， 用所建灾变预测模型z＾（1）（k+1）=9.580 75e0.26933k-8.580 75对未来马尾松毛虫越冬代灾变年份进行预测z＾（0）（11）=z＾（1）（11）-z（1）（10）=33.434 8， 由原始序号得知， 原始序列最后一次灾变是2011 年， z（0）（10）=23， 现预测下一次灾变年的序号z＾（0）（11）=33.434 8， 即从2011 年马尾松毛虫越冬代灾变年算起， 再过10 a 即2021 年为马尾松毛虫越冬代大发生年， 同理可预测之后年份z（t+1）=33.434 8，z（t+2）=43.769 1，z（t+3）=57.297 6，z（t+4）=75.007 5，z（t+5）=98.191 4。

2.3 马尾松毛虫1 代幼虫灾变预测

与上述马尾松毛虫越冬代虫口数GM（1,1）灾变预测模型相似， 由表1 根据阈值大于15 头·株-1可以得到马尾松毛虫1 代幼虫虫口数GM（1,1）灾变预测模型。

原始数列为x（0）={x（0）（1）， x（0）（2）， ...， x（0）（28）}={18.3， 42.5， 7.6， 7.8， 7.4， 38.6， 19.6， 40.1， 7.1，6.6， 7.2， 6.1， 8.6， 18.6， 7.7， 16.4， 12.0， 12.6， 12.1， 12.8， 16.2， 16.1， 7.4， 6.7， 6.7， 6.3， 6.6}。

对原始数列作灾变映射， 按ξ≥15 得灾变数列， 最终得到灾变日期集={1， 2， 6， 7， 8， 15， 17，22， 23}即为建模时所需的原始数列。

作OAG得：OAG{z（0）}=z（1）={1， 3， 9， 16， 24， 39， 56， 78， 101}。

按z（1）建模：即a＾=-0.241 87，u=4.155 67，u/a=-17.181 8，z（0）（1）=1，由此得出马尾松毛虫1 代幼虫虫口数的GM（1,1）灾变预测模型：

将根据马尾松毛虫1 代幼虫虫口数的GM（1,1）灾变预测模型得到的拟合值、 误差等列于表3。

表3 马尾松毛虫1 代幼虫虫口数预测的拟合值与误差Table 3 Fitting value and error of prediction of population number of first generation larvae of D. punctatus

对观察值和拟合值的差异进行t检验，t=0.218 6，df=16 时，t0.05=2.12，t＜t0.05， 两者差异不显著， 平均误差（2 个均数之差）只有观察值均数的6.49%。 实际观察值与通过GM（1,1）灾变预测模型得到的拟合值之间误差较小， 各灾变年的精度平均为84.85%， 总体精度为92.34%。 由此可以通过这个模型求得未来5 个时刻的预测值， 用所建灾变预测模型z＾（1）（k+1）=18.181 8e0.24187k-17.181 8 对未来马尾松毛虫一代幼虫灾变年份进行预测z＾（0）（10）=z＾（1）（10）-z（1）（9）=34.443 8， 由原始序号得知， 原始序列最后一次灾变是2011年，z（0）（9）=23， 现预测下一次灾变年的序号z＾（0）（10）=34.443 8， 即从2011 年马尾松毛虫1 代幼虫灾变年算起， 再过11 a 即2022 年为马尾松毛虫一代幼虫大发生年， 同理可以预测之后年份z（t+1）=34.443 8，z（t+2）=43.868 4，z（t+3）=55.871 7，z（t+4）=71.159 5，z（t+5）=90.630 2。

2.4 马尾松毛虫2 代幼虫灾变预测

与马尾松毛虫越冬代和1 代幼虫虫口数GM（1,1）灾变预测模型相似， 由表1 根据阈值大于15 头·株-1可以得到马尾松毛虫2 代幼虫虫口数GM（1,1）灾变预测模型。

原始数列为x（0）={x（0）（1）， x（0）（2）， ...， x（0）（28）}={34.1， 18.4， 13.2， 11.8， 26.5， 40.4， 19.6， 52.6，12.4， 11.6， 12.3， 12.6， 16.4， 18.8， 22.6， 12.6， 11.8， 12.1， 13.2， 13.8， 26.8， 16.7， 8.9， 7.2， 7.4，7.4， 8.1}。

对原始数列作灾变映射， 按ξ≥15 得灾变数列， 最终得到灾变日期集z={1， 2， 5， 6， 7， 8， 14，15， 16， 22， 23}即为建模时所需的原始数列。

作OAG得：={1， 3， 8， 14， 21， 29， 43， 58， 74， 96， 119}。

按z（1）建模：即a＾=-0.197 58，u=3.778 40，u/a=-19.123 7，z（0）（1）=1， 由此得出马尾松毛虫2 代幼虫虫口数的GM（1,1）灾变预测模型：

将根据马尾松毛虫2 代幼虫虫口数的GM（1,1）灾变预测模型得到的拟合值、 误差等列于表4。

对观察值和拟合值的差异进行t检验，t=0.195 2，df=20 时，t0.05=2.09，t＜t0.05， 两者差异不显著， 平均误差只占观察值均数的5.4%。 实际观察值与通过GM（1,1）灾变预测模型得到的拟合值之间误差较小，各灾变年的精度平均为84.08%， 总体平均精度为94.09%。 由此可以通过这个模型求得未来5 个时刻的预测值， 用所建灾变预测模型对未来马尾松毛虫2 代幼虫灾变年份进行预测z＾（0）（12）=z＾（1）（12）-z（1）（11）=31.704 2， 由原始序号得知， 最后一次灾变是2011 年，z（0）（11）=23， 现预测下一次灾变年的序号z＾（0）（12）=31.704 2， 即从2011 年马尾松毛虫2 代幼虫灾变年算起， 再过9 a 即2020 年为马尾松毛虫2 代幼虫大发生年， 同理可预测之后年份z（t+1）=31.704 2，z（t+2）=38.629 8，z（t+3）=47.068 4，z（t+4）=57.350 3，z（t+5）=69.878 2。

表4 马尾松毛虫2 代幼虫虫口数预测的拟合值与误差Table 4 Fitting value and error of prediction of population number of second generation larvae of D. punctatus

3 小结与讨论

利用灰色灾变模型预测法对安徽省潜山县1989-2016 年马尾松毛虫幼虫越冬代、 1 代和2 代累计虫口的灾变进行预测， GM（1,1）灾变预测模型分别是e0.24187k-17.181 8 和根据预测模型求得的已知年份的拟合值与观察值进行t检验，t=0.101 7， 0.218 6 和0.195 2， 均小于t0.05， 差异均不显著， 总体平均精度达96.25%、92.34%和94.09%， 各灾变年精度平均值为84.40%、 84.85%和84.08%， 误差极小。 说明灾变预测法对马尾松毛虫幼虫发生量灾变的预测预报是一种较理想的预报方法。 从2011 年算起， 预测越冬代灾变年为2021 年， 1 代灾变年为2022 年， 2 代灾变年为2020 年。

害虫种群的消长， 其本质是害虫种群与环境之间进行物质交换和能量传递的结果， 它既含有已知的信息， 又含有未知的或非确知的信息， 实际上是一个灰色系统。 灰色系统灾变预测一般是指某种大于阈值的灾害在哪一年份出现， 灰色系统灾变预测已经运用到国民经济各个方面的预测， 在植物保护方面也被大量利用。 在农业害虫的预测上， 预测结果与实际情况吻合度很高［11-13］。 灰色系统灾变预测结果在林业害虫的预测上也有报道， 吴秋芳等［14］对河南省安阳县杨树Populus苗圃地下害虫进行灰色灾变预测，平均精度达90.40%。 宋丁权等［15］用灰色理论中的灾变模型以及灰色区间模型对马尾松毛虫灾级出现的年份进行了模拟、 预测和分析， 也取得了较好的效果。 与中国马尾松毛虫发生的其他省区相比， 安徽省马尾松毛虫发生的偶灾区， 由于安徽省马尾松林面积较大， 加之还有一定面积的湿地松和火炬松的纯林和混交林， 为马尾松毛虫的猖獗提供了发生的基础条件， 安徽省地处亚热带和温带的过渡带， 温、 湿、降雨等气候条件均适宜于马尾松毛虫的发生， 马尾松毛虫发生的森林环境， 常灾区绝大部分为海拔100～300 m 的丘陵地区， 偶灾区为海拔300～500 m 的山区， 海拔500 m 以上的山区很少有马尾松毛虫发生， 安徽省为海拔400～800 m 的偶灾区［1］。 马尾松毛虫在遇到适宜的发生条件， 极易爆发成灾， 所以对马尾松毛虫猖獗成灾的预报有着重要的意义。 本研究预测了安徽省潜山县马尾松毛虫越冬代、 1 代和2代幼虫发生的严重年份， 平均精度分别为96.25%、 92.34%和94.09%。 灰色灾变模型预测法预测周期长， 预测结果精度高， 是较理想的预测方法。 对灰色灾变模型的应用在原始序列数据中应有一定的时间长度， 这样精度才能提高， 再者灾变模型随着时间的推移， 也需要不断更新和修正。