馬尾松毛蟲災情指數的方差分析周期外推預報

余燕 張書平 畢守東

摘要 為了提高馬尾松毛蟲Dendrolimus punctatus Walker災情預報結果的準確性，為提高防治效果提供科學依據， 本文采用方差分析周期外推預報法，以安徽省潛山縣的馬尾松毛蟲各代幼蟲累計發生量和當代發生面積計算求得的災情指數為依據，進行方差分析周期外推預報，并對預報結果進行驗證。預報結果的歷史符合率，1989-2016年全年災情指數預報結果歷史符合率為89.29%;預報2017年的結果與實況一致。1989-2016年的越冬代歷史符合率為89.29%，2017年預報值與實際完全一致。1997-2016年的一代預報結果歷史符合率為100%。1989-2016年的二代預報結果的歷史符合率為85.19%。2017年由于越冬代飛防使用滅幼脲，使當年一、二代蟲口大幅度降低，預報值高于實況。方差分析周期外推預報法對馬尾松毛蟲的災情指數預報是一種較理想的預報方法。

關鍵詞 馬尾松毛蟲幼蟲; 災情指數; 方差分析周期外推預報

中圖分類號： S 431.2

文獻標識碼： ADOI： 10.16688/j.zwbh.2018328

Abstract In order to improve the accuracy of disaster forecast of Dendrolimus punctatus， and provide scientific basis for improving the control efficacy， the periodic extrapolation prediction method of variance analysis was performed on the basis of the damage index calculated by the cumulative amount of larvae and the area of occurrence in each generation of D.punctatus in Qianshan county， Anhui province， and the prediction results were verified. The historical conformity rate of the annual disaster index forecast results was 89.29% for 1989-2016， and the forecasts for 2017 were in line with the actual data. The historical conformity rates of the overwintering and 2nd generations were 89.29% and 85.19% for 1989-2016， respectively， and the forecasts were in line with the actual data for 2017， while the historical conformity rate of the 1st generation was 100% in 1997-2016. In 2017， due to the application of chlorbenzuron by aerial spraying， the 1st and 2nd generations were greatly reduced and the forecast value was higher than the actual value. The method of variance analysis and periodic extrapolation is an ideal method for predicting the occurrence of D.punctatus.

Key words Dendrolimus punctatus larvae; damage index; analysis of variance periodic extrapolation prediction

馬尾松毛蟲Dendrolimus punctatus Walker分布于皖、豫、川、黔、陜、滇、贛、蘇、湘、浙、閩、粵、瓊、桂和臺灣等省，主要為害馬尾松Pinus massoniana Lamb.，還為害黑松P.thunbergii Parl.、火炬松P.taeda Linn.、濕地松 P.elliottii Engelm.、晚松P.rigida var. serotina （Michxa） Loud ex Hoopes、海南松P.fenzeliana Hand.Mazz.等松屬植物。20世紀中葉馬尾松毛蟲是我國發生最廣、危害面積最大，經常猖獗成災的森林害蟲。該蟲不但影響林業生產，還危害人體健康[14]。進入21世紀，森林管理上采用了封山育林、混交、間作等措施，優化了森林生態系統，增加了物種多樣性以及科學地運用綜合治理措施，使馬尾松毛蟲的危害得到有效控制，但該蟲具有巨大的繁殖潛能，遇到有利條件極易暴發成災，對其監測絲毫不能放松。

馬尾松毛蟲一年發生2～4代，發生世代的多少隨不同地區而異。馬尾松毛蟲發生的預測預報是對其進行綜合防治的基本工作;森林保護工作者分別采用不同的預測方法預測馬尾松毛蟲的發生量、蟲害等級以及發生類別、發生空間格局，為馬尾松毛蟲的綜合防治工作提供了有力支持[515]。由于各地氣象條件、植被條件和地形地貌的差異，馬尾松毛蟲的發生特點也不完全相同。對安徽省潛山縣馬尾松毛蟲1983-2016年的各世代幼蟲為害的災情指數鮮見有具體的系統研究報道，本研究采用方差分析周期外推預報法預報馬尾松毛蟲幼蟲為害的災情指數，以期為馬尾松毛蟲的綜合治理提供科學依據。

1 材料與方法

1.1 材料來源

馬尾松毛蟲資料來自國家林業局的森林病蟲中心測報點——安徽省潛山縣森林病蟲防治站，資料的時間跨度為1983-2016年。根據國家林業局森林病蟲防治總站制定的調查方法[16]進行調查，即采用踏查和詳查相結合的方法，對不同蟲態采取不同的觀測方法。沿林斑線、林道、公路、鐵路等線路調查，目測發生范圍、為害狀況，發現蟲情或災情立即設臨時標準地，采取平行線抽樣法抽取20株標準株詳查。卵期調查是在雌蛾羽化高峰后1～2 d調查平均卵塊數，卵塊平均粒數;幼蟲期調查，1～2齡幼蟲調查枯黃卷曲的枝數，推算幼蟲數，3齡以上的幼蟲3 m以下的小樹直接調查合計樹冠上的幼蟲數，大樹用“蟲糞粒推算法”調查，幼蟲越冬期間調查樹干基部樹皮縫中的幼蟲數推算全部蟲口。蛹期調查，在結繭盛期后2～5 d剖繭，調查雌雄性比、平均雌蛹重、死亡率、寄生率;成蟲期調查，在成蟲羽化前2～3 d及羽化期用黑光燈誘集，統計其誘集量和雌雄性比。

1.2 災情指數計算公式

為了便于分析馬尾松毛蟲幼蟲為害嚴重程度，考慮世代累計幼蟲蟲口和當代發生面積兩個因素，提出了利用災情指數為依據進行分析。災情指數：Z=xi×yi×nxk×yH×n×100%。xi為第i年某世代累計蟲口，yi為第i年某世代馬尾松毛蟲發生面積，n為平均每公頃松樹株數，k為某代或全年最大的蟲口，H為最大的發生面積。潛山縣松林面積為6.986 7萬hm2，1983年以來發生最大面積的是1996年的二代，其發生面積為33 381.4 hm2，占全縣松林面積的47.78%，最大蟲口為1996年，二代為52.6頭/株。災情指數分為以下4級：Z≤15%，輕度發生;15.1%≤Z≤30%，中等發生;30.1%≤Z≤45%，重度發生;Z>45%，特大發生。

1.3 方差分析周期外推預報法

方差分析周期外推預報是時間序列分析方法的一種[1]，吳勁松[18]、張梅等[19]、楊長登[20]、武漢中心氣象臺等[23]、康曉慧等[21]和Kwon等[22]將該法應用于氣象、農業、物理等方面的預測。如果對某種害蟲的災情指數進行等時距觀測，就得到一列觀測值y1，y2，y3，…，yn，即時間序列。如果該序列存在著周期長度為T的周期振蕩，則當把序列排成T個數值一行時，得到T列數值因為每一列的數值都是在周期振蕩的同一位相觀測到的，每列數值內即組內（如y1、yT+1、…、yKT+1）差異不大，而各列平均值之間即組間的差異則較大，即各列間的差異較各列內的大。即組間均方差大于組內均方差。反之，如果序列的周期長度為T，在排列觀測值時不是按T個數值一行，而是按t（t≠T）個數值一行，得到t列數值，由于這t列數值中每一列都是由不同位相觀測值組成的，所以，各列平均值間就不會有顯著差異，利用以上原理，如果當我們不知道序列中是否存在周期振蕩時，就可以取不同的t值（試測周期，t的最大值應取小于或等于N/2，N為樣本數），將序列排列起來，用方差分析檢驗各列平均值之間是否存在顯著差異，如果檢測結果表明t取某些數值時，各列平均值間差異顯著，則該序列為周期振蕩序列，其周期長度為t。

在實際時序中常包括若干個周期，通常取最顯著的（即F值最大的）試測周期作為第一周期，各列平均值作為第一周期分量。然后從原序列中減去第一周期分量，則得到第一周期余序列，對第一周期余序列繼續進行方差分析，確定其中各列平均值差異最顯著的周期為第二周期，相應各列平均值為第二周期分量。第一周期余序列減去第二周期分量，得到第二周期余序列，如此反復，重復試測周期，進行分析直到余序列沒有顯著周期時為止。預報時將各周期分量疊加，然后外延即可求出未來時段的災情指數預報值。

2 結果與分析

將潛山縣1983-2016年越冬代幼蟲累計蟲口動態列于表1，其中1998年由于前兩年馬尾松毛蟲暴發，用滅幼脲采取飛防措施導致該年沒有發生，因此統計值為0，根據蟲口變動特點，劃分為1983-1998年和1999-2016年兩個時段。越冬代幼蟲1983-2016年蟲口數量變動很大，1983-1998年除了1992年12.2頭/株，1993年為7.8頭/株和1998年0頭/株外，其他年份都在12.2～42.1之間，1983-1998年平均達24.76頭/株。1999-2016年蟲口較低，18年平均8.43頭/株。一代幼蟲蟲口數量變化與越冬代相似，1983-1998年除了1991、1992、1993、1997和1998年數量較少外，16年平均達23.16頭/株，1999-2016年平均蟲口為10.37頭/株。二代幼蟲蟲口數量1983-1998年平均為24.98頭/株，1999-2016年平均蟲口為13.35頭/株。總之1983-1997年各代幼蟲累計蟲口數量顯著高于1999-2016年的蟲口數量。

將1983-2016年馬尾松毛蟲的發生面積列于表2，1983-1998年和1999-2016年兩個時段，各年各代發生面積差別很大。越冬代發生面積：1984年高達13 600 hm2，1985年為11 541.47 hm2。1983-1998年平均發生面積為4 477.42 hm2。1999-2016年19年平均發生面積為1 929.57 hm2。一代發生面積：1983-1998年發生面積較大的年份主要集中在1984-1987年，4年平均發生面積高達10 697.5 hm2。1983-1998年平均發生面積3 899.27 hm2，1999-2016年平均發生面積1 506.41 hm2。二代發生面積：1983-1998年之間發生面積大的年份是1983、1984、1985、1987和1996年，分別為9 162.33、8 163.67、10 393.33、17 020和33 381.4 hm2。16年平均發生面積為6 605.94 hm2。1999-2016年平均發生面積為4 755.65 hm2。這些波動都體現在周期振蕩上。總之不管是發生面積還是各代幼蟲累計蟲口都是1983-1998年大，波動也大，1999-2016年的波動較小。由于波動幅度大小不一樣，因此要對各時間段逐一進行周期試測。

2.1 全年的災情指數動態

為了便于分析全年的馬尾松毛蟲災情情況，將全年越冬代、一代、二代發生面積匯總，蟲口數量采用越冬代、一代和二代之和的平均值。根據災情指數公式，1983-2016年最大的全年累計發生面積為1996年的42 224.73 hm2，年度平均蟲口最多的是1996年43.77頭/株，兩者之積作為全年災情指數的分母，計算年度災情指數列于表3，并繪于圖1。

災情指數大于45%的有1985年、1987年和1996年，屬于特大發生年。災情指數大于30%小于45%的有1984年、1986年和2010年，是重度發生年;災情指數大于15%小于30%的有1983年、1988年、1989年、1994年和2003年，屬于中度發生年;災情指數均小于15%，屬于輕發年，一共23年，占34年的67.65%，由于歷史的原因，亂砍濫伐，森林生態系統遭到嚴重破壞，所以1983-1990年災情指數一直很高。由于1992年、1993年全年平均蟲口低，發生面積小，疏于防治，致使1994年和1995年蟲口數量和發生面積都在上升，1996年大暴發，災情指數達到100%。1996年和1997年兩年連續飛防大大減小了1998年的災情指數。進入21世紀以來由于采取了封山育林等一系列科學管理措施，優化了森林生態系統，使馬尾松毛蟲各代的蟲口數量都較低。

根據圖1反映的信息對1989年以來全年災情指數進行周期試測，t=3，F=0.558;t=4，F=0.429;t=5，F=0.56;t=6，F=0.709;t=7，F=4.765;t=8，F=1.356;t=9，F=0.755;t=10，F=0.625;t=11，F=1.812;t=12，F=1.072;t=13，F=0.959;t=14，F=1.158。查F表，F0.01（df=6， 21）=3.81，只有t=7時，F>F0.01（6，21），組間差異極顯著。將其第一周期t=7的周期分量列于表4，求得第一周期余序列也列于表4，對第一周期余序列繼續周期試測，t=3，F=2.137;t=4，F=1.094;t=5，F=1.083;t=6，F=1.479;t=8，F=0.994;t=9，F=1.026;t=10，F=0.945;t=11，F=1.218;t=12，F=1.026;t=13，F=0.786;t=14，F=1.35，F值均小于F0.05，差異均不顯著，經評判，1989年、1994年和2003年與實況不符，28年的歷史符合率為89.29%。根據方差分析周期外推預測結果2017年為45.78%，2017年全年發生面積21 149.533 3 hm2，年平均蟲口為5.033 3頭/株，災情指數實況為57.57%，預測結果與實況相符，均是特大發生年。

2.2 1989-2016年潛山縣馬尾松毛蟲越冬代災情指數

通過周期試測1989-2016年有周期振蕩現象存在，越冬代最大蟲口為1989年42.5頭/株，越冬代最大發生面積為1984年13 600 hm2;將兩者乘積作為災情指數計算時的分母，將求得各年災情指數列于表5，并繪于圖2，進行周期試測，t=3，F=0.755;t=4，F=0.207;t=5，F=0.702;t=6，F=0.826;t=7，F=5.323;t=8，F=0.918;t=9，F=0.441;t=10，F=0.948;t=11，F=0.879;t=12，F=0.725;t=13，F=1.106，t=14，F=5.213。查F表，F0.01（6，21）=3.81，只有t=7時，F>F0.01（6，21），組間差異極顯著。將其第一周期t=7的周期分量列于表5，求得第一周期余序列也列于表5，對第一周期余序列繼續周期試測，t=3，F=1.577;t=4，F=0.662;t=5，F=1.647;t=6，F=1.381;t=8，F=0.498;t=9，F=0.93;t=10，F=1.196;t=11，F=1.102;t=12，F=0.813;t=13，F=0.478;t=14，F=1.429。t由3至12為周期，其F值均小于F0.05，表明第一周期余序列沒有周期振蕩現象存在，第一周期分量就是預報值，歷史符合率為89.29%。2017年越冬代發生面積為5 521.47 hm2，發生蟲口數量為12.8頭/株，災情指數為12.23%。2018年越冬代發生面積為5 790 hm2，發生蟲口數量為1.8頭/株，災情指數為1.8%，預測災情指數2017年為12.90%，2018年為24.75%，2017年預測值與實況相符，2018年預測值高一級，其原因是2017年越冬代進行了飛防，壓低了2018年的蟲口，使預測值大于實況值。

2.3 馬尾松毛蟲一代幼蟲災情指數

通過周期試測1997-2016年有周期振蕩現象存在，一代最大蟲口數量為1986年的46.1頭/株，一代最大發生面積為1987年的14 733.33 hm2。將兩者乘積作為災情指數計算時的分母，將求得的各年災情指數列于表6，并繪于圖3，進行周期試測，t=3時，F=0.529;t=4時，F=1.72;t=5時，F=0.488;t=6時，F=0.636;t=7時，F=1.393;t=8時，F=11.474;t=9時，F=1.242;t=10時，F=0.999。只有t=8為周期的F>F0.01（7，20）=3.7，t=8的周期分量和余序列也列于表6，繼續對t=8為周期的余序列進行周期試測，t=3時，F=0.151;t=4，F=0;t=5，F=0.797;t=6，F=0.479;t=7，F=0.882;t=9，t=0.848;t=10，t=2.327。表明第一周期余序列的t為3、4、5、6、7、9、10時F值都小于F0.05，即沒有周期振蕩現象存在，第一周期分量就是預報值。對其進行評判，28年的歷史符合率為100%。2017年一代發生面積為9 855.2 hm2，2017年一代蟲口數量為1.2頭/株，實際災情指數2017年為1.74，預測外推結果依次為0.31。2018年發生面積為2 333.33 hm2，2018年蟲口數量為1.6頭/株，實際災情指數2018年為0.5，預測外推結果為2.325。2017年和2018年災情實況和預測值均為輕度發生，預報結果準確。

2.4 馬尾松毛蟲二代幼蟲災情指數預報

通過周期試測1989-2016年有周期振蕩現象存在，二代最多蟲口數量為1996年的52.6頭/株，二代最大發生面積為1996年的33381.4 hm2。將兩者乘積作為災情指數的分母，求得的各年災情指數列于表7，并繪于圖4，進行周期試測，t=3，F=0.395;t=4，F=0.639;t=5，F=0.482;t=6，F=0.702;t=7，F=4.057;t=8，F=1.484;t=9，F=1.23;t=10，F=0.692;t=11，F=1.23;t=12，F=1.062;t=13，F=0.935;t=14，F=8.295。t為3、4、5、6、7、8、9、10、11、12、13、14時的F值只有t=14時，F>F0.01（13，14）=3.72，即差異極顯著，對t=14周期的余序列繼續進行周期試測，t=3，F=1.922;t=4，F=0.345;t=5，F=0.733;t=6，F=1.118;t=7，F=0;t=8，F=0.831;t=9，F=1.401;t=10，F=0.838;t=11，F=1.487;t=12，F=0.806;t=13，F=0.988。t為3、4、5、6、7、8、9、10、11、12、13時其F值均小于F0.05，表明第一周期余序列沒有周期振蕩存在，第一周期分量就是預報值。按照災情指數的分級，歷史符合率為100%，其中2003年和2010年為特大發生，其余均為輕度發生。2017年二代發生面積5 772.87 hm2，二代幼蟲蟲口數量1.1頭/株，災情指數為0.36，外推預測結果災情指數為12.935，實況和預報值均為輕度發生，外推災情指數高于實際發生的災情指數。由于2017年越冬代飛防使用了滅幼脲，所以使當年二代蟲口數量大幅下降，使實況災情指數降低。1989年以來預測結果歷史符合率為85.19%，只有1994、1996、2003和2010年的預測結果略低于實況。

3 小結與討論

本文采用方差分析周期外推預報法以安徽省潛山縣馬尾松毛蟲各代幼蟲累計發生量和當代發生面積為依據計算災情指數，對災情指數進行周期試測，并對越冬代、一代、二代和全年的災情指數預報值和實況值進行評判和比較，結果是越冬代的歷史符合率為89.29%，2017年的預報值與實況一致。一代歷史符合率為100%，二代的歷史符合率為85.19%。2017年的預報值比實況高一級，全年災情預報結果歷史符合率為89.29%。方差分析周期外推預報法較之于BP神經網絡法，馬爾科夫鏈法，災變預測法，列聯表分析法優，計算方便，預報結果較準確，且計算時采用分析的數據不是其他生態因子數據，對一些與害蟲發生有關的不易獲取的生態因子數據，作為預報要素，這些數據與預報量有關時則不易進行預報，而方差分析周期外推預報法可以克服這些困難，但方差分析周期外推預報法法只適用于預報因子有周期振蕩現象時的預報，這是其局限性。

對于一個縣域范圍總體上的蟲害災情評判，到目前為止還沒有一個較理想的標準，特別是森林害蟲。本文采用災情指數，基于危害蟲態的數量以及發生面積等因子計算災情指數。本文是基于安徽省潛山縣1983-2016年馬尾松毛蟲各代幼蟲發生數量和發生面積的歷史資料進行研究的，其他蟲害也可參考此法，災情指數通過外推可預報未來的發生嚴重程度，以便提前做好預防工作，變被動為主動，盡量減少蟲害損失。本研究發現越冬代的振蕩周期是7年，一代是8年，二代是14年，為什么三代的震蕩周期不一致？分析其原因，可能主要包括：害蟲的數量變化一是受氣象等物理環境因素影響，二是受管理措施影響，三是受天敵的種類和數量多少影響，這三種因子對越冬代、一代、二代的作用可能是不一樣的，因此出現越冬代、一代和二代的災情指數振蕩周期不一致的現象;其不同代幼蟲受外界因素的具體影響情況以及馬尾松毛蟲響應的內在因素尚需開展系統研究。

參考文獻

[1] 侯陶謙.中國松毛蟲[M].北京：科學出版社，1987：188191.

[2] 鄒運鼎，程扶玖，查光濟.松針內含物與馬尾松毛蟲生存發育關系的研究[J].林業科學，1990，26（2）：142148.

[3] 蕭剛柔.中國森林害蟲（第二版）[M].北京：林業出版社，1992：948953.

[4] 張真，李典謨.馬尾松毛蟲暴發機制分析[J].林業科學，2008，44（1）：140150.

[5] 張愛兵，陳建，王正軍，等.BP網絡模型和LOGIT模型在森林害蟲測報上的應用初報——以安徽省潛山縣馬尾松毛蟲為例[J].生態學報，2001，21（12）：21592165.

[6] PARK Y S，CEREGHINO R，COMPIN A. Applications of artificial neural networks for patterning and predicting aquatic insect species richness in running waters [J]. Ecological Modelling， 2003，160（3）： 265280.

[7] 賈春生.利用馬爾可夫鏈方法測報馬尾松毛蟲發生級別[J].東北林業大學學報， 2006， 34（5）：2122.

[8] ZHANG W J，ZHONG X Q，LIU G H. Recognizing spatial distribution patterns of grassland insects： neural network approaches [J].Stochastic Environmental Research and Risk Assessment，2008， 22（2）： 207216.

[9] 陳繪畫，王堅婭，徐志宏.基于響應面方法的馬尾松毛蟲發生量混沌特性檢測及其預測[J].東北林業大學學報，2011，39（9）：9496.

[10]田萬銀，徐華潮.浙江沿海防護林馬尾松毛蟲的預測預報模型[J].環境昆蟲學報，2012，34（4）：401406.

[11]費海澤，王鴻斌，孔祥波等.馬尾松毛蟲發生相關氣象因子篩選及預測[J].東北林業大學學報，2014，41（1）：136140.

[12]許章華，李聰慧，劉健.馬尾松毛蟲害等級的Fisher判別分析[J].農業機械學報，2014，45（6）：275283.

[13]王慶，畢猛，杜婷.基于氣象因子的馬尾松毛蟲發生率空間格局研究[J].林業科學研究，2016，29（2）：256260.

[14]余燕，李尚，王振興，等.馬尾松毛蟲幼蟲發生嚴重程度的預測研究[J].安徽農業大學學報，2017，44（5）：882893.

[15]周夏芝，王振興，余燕，等.馬尾松毛蟲幼蟲高峰期發生量的預測模型研究[J].應用昆蟲學報，2017，54（6）：10311043.

[16]國家林業局森林病蟲害防治總站.林業有害生物監測預報技術[J].北京：中國林業出版社，2013：117118.

[17]張孝羲，翟保平，牟極元，等.昆蟲生態及預測預報[M].第3版.北京：中國農業出版社，2002：262265.

[18]吳勁松.用方差分析周期及隨機時間序列法作赫章6-8月總降雨量預報[J].貴州氣象，1999，23（1）：1314.

[19]張梅，陳玉光，楊冰.方差分析周期疊加法預測農作物生長季積溫[J].現代農業科技，2017（21）：237240.

[20]楊長登.用方差分析周期疊加外推法預報年降水量[J].貴州氣象，1998，22（1）：2325.

[21]康曉慧，陳浩，張梅.3種時間序列分析模型在水稻稻瘟病預測中的應用[J].西北農林科技大學學報（自然科學版），2011，39（6）：173184.

[22]KWON J H. Crystal graphs and the combinatorics of young tableaux [M]. Handbook of Algebra. Vol 6， NorthHolland， New York， 2009： 473504.

[23]武漢中心氣象臺，武漢大學數學系.方差分析周期外推法在長期預報中的應用[J].數學學報，1974，17（3）：156163.

（責任編輯：田 喆）