-
马尾松毛虫Dendrolimus punctatus分布于安徽、河南、四川、贵州、陕西、云南、江西、江苏、湖南、浙江、福建、广东、海南和广西等, 主要危害马尾松Pinus massoniana, 还危害黑松Pinus thunbergii, 火炬松Pinus taeda, 湿地松Pinus elliottii, 晚松Pinus rigida var. 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代严重发生的年份, 以期为马尾松毛虫的综合治理提供科学依据。
HTML
-
马尾松毛虫每年越冬代、1代幼虫和2代幼虫的累计虫口数见表 1, 将其作为原始数据。
序号 年份 累计虫口数/(头·株-1) 越冬代 1虫 2虫 1 1989 42.5 18.3 34.1 2 1990 37.3 42.5 18.4 3 1991 18.5 7.6 13.2 4 1992 12.2 7.8 11.8 5 1993 7.8 7.4 26.5 6 1994 18.6 38.6 40.4 7 1995 18.2 19.6 19.6 8 1996 38.6 40.1 52.6 9 1997 42.1 7.1 12.4 10 1998 - - - 11 1999 6.2 6.6 11.6 12 2000 7.6 7.2 12.3 13 2001 6.2 6.1 12.6 14 2002 7.1 8.6 16.4 15 2003 15.7 18.6 18.8 16 2004 7.6 7.7 22.6 17 2005 6.4 16.4 12.6 18 2006 6.6 12 11.8 19 2007 6.4 12.6 12.1 20 2008 6.6 12.1 13.2 21 2009 8.4 12.8 13.8 22 2010 16.2 16.2 26.8 23 2011 16.6 16.1 16.7 24 2012 7.2 7.4 8.9 25 2013 7.2 6.7 7.2 26 2014 6.4 6.7 7.4 27 2015 6.2 6.3 7.4 28 2016 7.1 6.6 8.1 说明:"-"表示数据缺失 Table 1. Accumulative population of D. punctatus during wintering, first and second generation larvae
-
根据安徽省潜山县森林病虫害防治总站提供的1983-2016年马尾松毛虫发生情况的资料, 马尾松毛虫越冬代虫口数大于15头·株-1的年份为大发生年, 故以每株虫口数15头为灾变阈值。
-
原始数列为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得灾变数列:
ξ{x(0)}={42.5, 37.3, 18.5, 18.6, 18.2, 38.6, 42.1, 15.7, 16.2, 16.6}={x(0)(1), x(0)(2), x(0)(3), x(0)(6), x(0)(7), x(0)(8), x(0)(9), x(0)(15), x(0)(22), x(0)(23)}={x(0)(1′), x(0)(2′), x(0)(3′), x(0)(4′), x(0)(5′), x(0)(6′), x(0)(7′), x(0)(8′), x(0)(9′), x(0)(10′)}=xξ(0)。
然后再作灾变序号映射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}即为建模时所需的原始数列。
-
对灾变年序集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)建模:$\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(0)(1)=1, 由此得出马尾松毛虫越冬代虫口数的GM(1, 1)灾变预测模型:${\hat z^{(1)}}(k + 1)$=9.580 75e0.269 33k-8.580 75。
将根据马尾松毛虫越冬代虫口数的GM(1, 1)灾变预测模型得到的拟合值、误差等列于表 2。
序号 大于阈值的年份 虫口数/(头·株-1) 精度/% 观察值 拟合值 误差 x(1) 1989 1 1.00 0 100.00 x(2) 1990 2 2.96 -0.96 51.94 x(3) 1991 3 3.88 -0.88 70.78 x(4) 1994 6 5.07 0.93 84.58 x(5) 1995 7 6.64 0.36 94.91 x(6) 1996 8 8.70 -0.70 91.29 x(7) 1997 9 11.38 -2.38 73.50 x(8) 2003 15 14.90 0.10 99.36 x(9) 2010 22 19.51 2.49 88.68 x(10) 2011 23 25.54 -2.54 88.95 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个时刻的预测值, 用所建灾变预测模型${\hat z^{(1)}}(k + 1)$=9.580 75e0.269 33k-8.580 75对未来马尾松毛虫越冬代灾变年份进行预测${\hat z^{(0)}}(11) = {\hat z^{(1)}}(11)$-z(1)(10)=33.434 8, 由原始序号得知, 原始序列最后一次灾变是2011年, z(0)(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(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)建模:$\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(0)(1)=1, 由此得出马尾松毛虫1代幼虫虫口数的GM(1, 1)灾变预测模型:${\hat z^{(1)}}(k + 1)$=18.181 8e0.241 87k-17.181 8。
将根据马尾松毛虫1代幼虫虫口数的GM(1, 1)灾变预测模型得到的拟合值、误差等列于表 3。
序号 大于阈值的年份 虫口数/(头·株-1) 精度/% 观察值 拟合值 误差 x(1) 1989 1 1.00 0 100.00 x(2) 1990 2 4.97 -2.97 48.75 x(3) 1994 6 6.34 -0.34 94.40 x(4) 1995 7 8.07 -1.07 84.72 x(5) 1996 8 10.28 -2.28 71.53 x(6) 2003 15 13.09 1.91 87.27 x(7) 2005 17 16.67 0.33 98.07 x(8) 2010 22 21.23 0.77 96.52 x(9) 2011 23 27.04 -4.04 82.42 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个时刻的预测值, 用所建灾变预测模型${\hat z^{(1)}}(k + 1)$=18.181 8e0.241 87k-17.181 8对未来马尾松毛虫一代幼虫灾变年份进行预测${\hat z^0}(10) = {\hat z^{(1)}}(10)$-z(1)(9)=34.443 8, 由原始序号得知, 原始序列最后一次灾变是2011年, z(0)(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(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得:OAG{z(0)}=z(1)={1, 3, 8, 14, 21, 29, 43, 58, 74, 96, 119}。
按z(1)建模:$\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(0)(1)=1, 由此得出马尾松毛虫2代幼虫虫口数的GM(1, 1)灾变预测模型:${\hat z^{(1)}}(k + 1)$=20.123 7e0.197 58k-19.123 7。
将根据马尾松毛虫2代幼虫虫口数的GM(1, 1)灾变预测模型得到的拟合值、误差等列于表 4。
序号 大于阈值的年份 虫口数/(头·株-1) 精度/% 观察值 拟合值 误差 x(1) 1989 1 1.00 0 100.00 x(2) 1990 2 4.40 -2.40 19.80 x(3) 1993 5 5.36 -0.36 92.88 x(4) 1994 6 6.53 -0.53 91.23 x(5) 1995 7 7.95 -0.95 86.40 x(6) 1996 8 9.69 -1.69 78.89 x(7) 2002 14 11.81 2.19 84.33 x(8) 2003 15 14.38 0.62 95.90 x(9) 2004 16 17.53 -1.53 90.46 x(10) 2010 22 21.36 0.64 97.07 x(11) 2011 23 26.02 -3.02 86.87 Table 4. Fitting value and error of prediction of population number of second generation larvae of D. punctatus
对观察值和拟合值的差异进行t检验, t=0.195 2, df=20时, t0.05=2.09, t<t0.05, 两者差异不显著, 平均误差只占观察值均数的5.4%。实际观察值与通过GM(1, 1)灾变预测模型得到的拟合值之间误差较小, 各灾变年的精度平均为84.08%, 总体平均精度为94.09%。由此可以通过这个模型求得未来5个时刻的预测值, 用所建灾变预测模型${\hat z^{(1)}}(k + 1)$=20.123 7e0.197 58k-19.123 7对未来马尾松毛虫2代幼虫灾变年份进行预测${\hat z^{(0)}}(12) = {\hat z^{(1)}}(12)$-z(1)(11)=31.704 2, 由原始序号得知, 最后一次灾变是2011年, z(0)(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。