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

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

原始数据${x^{(0)}} = \left[ {\begin{array}{*{20}{l}} {x_{(1)}^{(0)}, }&{x_{(2)}^{(0)}, }&{ \cdots , x_{(n)}^{(0)}} \end{array}} \right]$, $x_{(i)}^{(0)}$指数列x⁽⁰⁾中第i年的发生数量。令ξ为异常值, 如果$x_{(i)}^{(0)} \ge \xi $为灾变, 则称这种情况为上灾变。如果$x_{(i)}^{(0)} \le \xi $为灾变, 则称为下灾变。一般上下灾变统称灾变, ξ为灾变阈值, 或灾变值^[10]。在害虫发生量和为害损失的预测预报中, 最关心的是上灾变, 也即预测数量或损失超过某个阈值出现的时间。x⁽⁰⁾中每一个数据$x_{(\mathit{k})}^{(0)}$应该说是数值$x_{(\mathit{k})}^{(0)}$与年序号k的二元组合, 记为[k, $x_{(\mathit{k})}^{(0)}$]。

如果在x⁽⁰⁾中有新的数列$x_\xi ^{(0)} = \left[ {x_{\left( {{1^\prime }} \right)}^{(0)}, \quad x_{\left( {{2^\prime }} \right)}^{(0)}, \cdots , x_{\left( {{n^\prime }} \right)}^{(0)}} \right]$, n′＜n是符合灾变定义的数列, 则对于上灾变$x_{\left( {{i^\prime }} \right)}^{(0)}$≥ξ, V′∈[1′, 2′, …, n′], 其中:i′为灾变年序号, V′为灾变年序号取值范围。

从x⁽⁰⁾到x_ξ⁽⁰⁾是一种映射, 记为ξ:[x⁽⁰⁾]→[x_ξ⁽⁰⁾]。为了建立灾变预测数据集, 将二元组合$x_{\left( {{i^\prime }} \right)}^{(0)}$记为[i′, $x_{\left( {{i^\prime }} \right)}^{(0)}$]。

为了便于辨认, 特记灾变年序号i′为z(i′)或z(i)。并作灾变号映射p, p[x_ξ⁽⁰⁾]=(ω); p[x_ξ⁽⁰⁾]=[z(1′), z(2′), …, z(n′)]=(z), 记z的1次累加生成数(accumulated generating operation, O_AG)为O_AG:(z)→(z⁽¹⁾), 则灾变预测的GM(1, 1)模型有如下算式:d(z)⁽¹⁾/dt+az⁽¹⁾=u, $\hat{a}=\left[\begin{array}{l} {a} \\ {u} \end{array}\right], \quad \hat{a}=\left(\boldsymbol{B}^{\mathrm{T}} \boldsymbol{B}\right)^{-1} \boldsymbol{B}^{\mathrm{T}} y_{N}$, 其中:t为年份, a为发展灰数, u为内生控制灰数, $\hat a$为待估灰数, B为累加矩阵, B^T为转置矩阵, y_N为常数项向量。

灾变模型为${\hat z^{(1)}}(k + 1) = \left[ {{z^{(0)}}(1) - u/a} \right]{{\rm{e}}^{ - ak}} + u/a$, 其中:z⁽⁰⁾为原始灾变年序号, 如z⁽⁰⁾(1)为第1次发生灾变的年份, 即1989年, ${\hat z^{(1)}}(k + 1)$为灾变年序号。

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

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

原始数列为x⁽⁰⁾={x⁽⁰⁾(1), x⁽⁰⁾(2), …, x⁽⁰⁾(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得灾变数列:

ξ{x⁽⁰⁾}={42.5, 37.3, 18.5, 18.6, 18.2, 38.6, 42.1, 15.7, 16.2, 16.6}={x⁽⁰⁾(1), x⁽⁰⁾(2), x⁽⁰⁾(3), x⁽⁰⁾(6), x⁽⁰⁾(7), x⁽⁰⁾(8), x⁽⁰⁾(9), x⁽⁰⁾(15), x⁽⁰⁾(22), x⁽⁰⁾(23)}={x⁽⁰⁾(1′), x⁽⁰⁾(2′), x⁽⁰⁾(3′), x⁽⁰⁾(4′), x⁽⁰⁾(5′), x⁽⁰⁾(6′), x⁽⁰⁾(7′), x⁽⁰⁾(8′), x⁽⁰⁾(9′), x⁽⁰⁾(10′)}=x_ξ⁽⁰⁾。

然后再作灾变序号映射p, p{x_ξ⁽⁰⁾}={z}, px_ξ⁽⁰⁾={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}即为建模时所需的原始数列。

对灾变年序集z建立GM(1, 1)模型, 按原始数列z={1, 2, 3, 6, 7, 8, 9, 15, 22, 23}作O_AG(一次累加生成算子)得:

O_AG{z⁽⁰⁾}=z⁽¹⁾={z⁽¹⁾(1), z⁽¹⁾(2), z⁽¹⁾(3), z⁽¹⁾(4), z⁽¹⁾(5), z⁽¹⁾(6), z⁽¹⁾(7), z⁽¹⁾(8), z⁽¹⁾(9), z⁽¹⁾(10)}={1, 3, 6, 12, 19, 27, 36, 51, 73, 96}。

按z⁽¹⁾建模:$\hat a = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right)^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{y_N} = \left[ {\begin{array}{*{20}{c}} { - 0.269\;33}\\ {2.311\;06} \end{array}} \right]$, 即$\hat a$=-0.269 33, u=2.311 06, u/a=-8.580 75, z⁽⁰⁾(1)=1, 由此得出马尾松毛虫越冬代虫口数的GM(1, 1)灾变预测模型:${\hat z^{(1)}}(k + 1)$=9.580 75e^{0.269 33k}-8.580 75。

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

对观察值和拟合值间的差异进行t检验, t=0.101 7, d_f=18时, t_0.05=2.10, t＜t_0.05, 差异不显著, 平均误差(2个均数之差)只有观察值均数的3.40%。实际观察值与通过GM(1, 1)灾变预测模型得到的拟合值之间误差较小, 各灾变年的平均精度为84.40%, 总体精度(观察值的平均值减去误差的平均值与观察值的平均值之比)为96.25%。由此可以通过这个模型求得未来5个时刻的预测值, 用所建灾变预测模型${\hat z^{(1)}}(k + 1)$=9.580 75e^{0.269 33k}-8.580 75对未来马尾松毛虫越冬代灾变年份进行预测${\hat z^{(0)}}(11) = {\hat z^{(1)}}(11)$-z⁽¹⁾(10)=33.434 8, 由原始序号得知, 原始序列最后一次灾变是2011年, z⁽⁰⁾(10)=23, 现预测下一次灾变年的序号${\hat 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。

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

原始数列为x⁽⁰⁾={x⁽⁰⁾(1), x⁽⁰⁾(2), ..., x⁽⁰⁾(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}即为建模时所需的原始数列。

作O_AG得:O_AG{z⁽⁰⁾}=z⁽¹⁾={1, 3, 9, 16, 24, 39, 56, 78, 101}。

按z⁽¹⁾建模:$\hat a = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right)^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{y_N} = \left[ {\begin{array}{*{20}{c}} { - 0.241\;87}\\ {4.155\;67} \end{array}} \right]$, 即$\hat a$=-0.241 87, u=4.155 67, u/a=-17.181 8, z⁽⁰⁾(1)=1, 由此得出马尾松毛虫1代幼虫虫口数的GM(1, 1)灾变预测模型:${\hat z^{(1)}}(k + 1)$=18.181 8e^{0.241 87k}-17.181 8。

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

对观察值和拟合值的差异进行t检验, t=0.218 6, d_f=16时, t_0.05=2.12, t＜t_0.05, 两者差异不显著, 平均误差(2个均数之差)只有观察值均数的6.49%。实际观察值与通过GM(1, 1)灾变预测模型得到的拟合值之间误差较小, 各灾变年的精度平均为84.85%, 总体精度为92.34%。由此可以通过这个模型求得未来5个时刻的预测值, 用所建灾变预测模型${\hat z^{(1)}}(k + 1)$=18.181 8e^{0.241 87k}-17.181 8对未来马尾松毛虫一代幼虫灾变年份进行预测${\hat z^0}(10) = {\hat z^{(1)}}(10)$-z⁽¹⁾(9)=34.443 8, 由原始序号得知, 原始序列最后一次灾变是2011年, z⁽⁰⁾(9)=23, 现预测下一次灾变年的序号${\hat 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。

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

原始数列为x⁽⁰⁾={x⁽⁰⁾(1), x⁽⁰⁾(2), ..., x⁽⁰⁾(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}即为建模时所需的原始数列。

作O_AG得:O_AG{z⁽⁰⁾}=z⁽¹⁾={1, 3, 8, 14, 21, 29, 43, 58, 74, 96, 119}。

按z⁽¹⁾建模:$\hat a = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right)^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{y_N}\left[ {\begin{array}{*{20}{c}} { - 0.197\;58}\\ {3.778\;40} \end{array}} \right]$, 即$\hat a$=-0.197 58, u=3.778 40, u/a=-19.123 7, z⁽⁰⁾(1)=1, 由此得出马尾松毛虫2代幼虫虫口数的GM(1, 1)灾变预测模型:${\hat z^{(1)}}(k + 1)$=20.123 7e^{0.197 58k}-19.123 7。

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

对观察值和拟合值的差异进行t检验, t=0.195 2, d_f=20时, t_0.05=2.09, t＜t_0.05, 两者差异不显著, 平均误差只占观察值均数的5.4%。实际观察值与通过GM(1, 1)灾变预测模型得到的拟合值之间误差较小, 各灾变年的精度平均为84.08%, 总体平均精度为94.09%。由此可以通过这个模型求得未来5个时刻的预测值, 用所建灾变预测模型${\hat z^{(1)}}(k + 1)$=20.123 7e^{0.197 58k}-19.123 7对未来马尾松毛虫2代幼虫灾变年份进行预测${\hat z^{(0)}}(12) = {\hat z^{(1)}}(12)$-z⁽¹⁾(11)=31.704 2, 由原始序号得知, 最后一次灾变是2011年, z⁽⁰⁾(11)=23, 现预测下一次灾变年的序号${\hat 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。

利用灰色灾变模型预测法对安徽省潜山县1989-2016年马尾松毛虫幼虫越冬代、1代和2代累计虫口的灾变进行预测, GM(1, 1)灾变预测模型分别是${\hat z^{(1)}}(k + 1)$=9.580 75e^{0.269 33k}-8.580 75, ${\hat z^{(1)}}(k + 1)$=18.181 8 e^{0.241 87k}-17.181 8和${\hat z^{(1)}}(k + 1)$=20.123 7e^{0.197 58k}-19.123 7。根据预测模型求得的已知年份的拟合值与观察值进行t检验, t=0.101 7, 0.218 6和0.195 2, 均小于t_0.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%。灰色灾变模型预测法预测周期长, 预测结果精度高, 是较理想的预测方法。对灰色灾变模型的应用在原始序列数据中应有一定的时间长度, 这样精度才能提高, 再者灾变模型随着时间的推移, 也需要不断更新和修正。