災(zāi)變模型在馬尾松毛蟲幼蟲發(fā)生量預(yù)報中的應(yīng)用

張書平， 余 燕， 畢守東， 周夏芝， 鄒運鼎， 張國慶， 張 楨， 方國飛， 宋玉雙

（1. 安徽農(nóng)業(yè)大學(xué) 理學(xué)院， 安徽 合肥230036； 2. 安徽農(nóng)業(yè)大學(xué) 林學(xué)與園林學(xué)院， 安徽 合肥230036；3. 安徽省潛山縣林業(yè)局， 安徽 潛山246300； 4. 國家林業(yè)和草原局 森林病蟲害防治總站， 遼寧 沈陽110034）

馬尾松毛蟲Dendrolimus punctatus分布于安徽、 河南、 四川、 貴州、 陜西、 云南、 江西、 江蘇、 湖南、 浙江、 福建、 廣東、 海南和廣西等， 主要危害馬尾松Pinus massoniana， 還危害黑松Pinus thunbergii， 火炬松Pinus taeda， 濕地松Pinus elliottii， 晚松Pinus rigidavar.serotina和海南松Pinus fenzeliana等松屬Pinus植物。 20 世紀中葉在中國森林害蟲中馬尾松毛蟲是發(fā)生最廣、 危害面積最大、 經(jīng)常猖獗成災(zāi)的害蟲。 在廣大丘陵地區(qū)蟲害此起彼伏， 為害時如同火燒， 造成了巨大的經(jīng)濟效益和生態(tài)效益損失。該蟲不但影響林業(yè)生產(chǎn)， 還危害人身健康［1-4］， 人們的林事活動中接觸馬尾松毛蟲毒毛， 容易引發(fā)皮炎和關(guān)節(jié)腫痛。 進入21 世紀， 封山育林、 混交、 間作等措施優(yōu)化了森林生態(tài)環(huán)境， 使馬尾松毛蟲的危害得到有效控制， 但該蟲具有強大的繁殖潛力， 遇到有利的生態(tài)環(huán)境極易爆發(fā)成災(zāi)， 不能放松對它的監(jiān)測。 馬尾松毛蟲1 a 發(fā)生2～4 代， 發(fā)生世代的多少隨不同地方而異。 在河南省信陽地區(qū)1 a 發(fā)生2 代為主， 在長江流域諸省1 a 發(fā)生2～3 代， 而在廣東、 廣西、 福建南部1 a 發(fā)生3～4 代， 海南1 a 發(fā)生4～5代［1］。 安徽潛山縣1 a 發(fā)生3 代， 即4-6 月上旬為越冬代， 6 月上旬至8 月中下旬為1 代， 8 月中下旬至12 月為2 代。 馬尾松毛蟲發(fā)生的預(yù)測預(yù)報是對其進行綜合防治的基本工作。 研究者［5-9］分別采用不同的預(yù)測方法預(yù)測馬尾松毛蟲的發(fā)生量、 蟲害等級、 發(fā)生類別、 發(fā)生空間格局， 為馬尾松毛蟲的綜合防治工作提供了有力支持。 由于各地氣象條件、 植被條件和地形地貌等不同， 馬尾松毛蟲的發(fā)生特點也不完全相同。 為了有效地防治馬尾松毛蟲， 本研究采用灰色理論中的災(zāi)變預(yù)測法研究了安徽省潛山縣馬尾松毛蟲幼蟲越冬代、 1 代和2 代嚴重發(fā)生的年份， 以期為馬尾松毛蟲的綜合治理提供科學(xué)依據(jù)。

1 材料與方法

1.1 材料來源

馬尾松毛蟲資料來自國家林業(yè)和草原局森林病蟲中心測報點——安徽省潛山縣森林病蟲防治站， 時間跨度為1989-2016 年， 其中1998 年缺失。 采用踏查和詳查相結(jié)合的方法， 沿林班線、 林道、 公路、鐵路等線路調(diào)查， 目測發(fā)生范圍、 危害狀況， 發(fā)現(xiàn)蟲情或災(zāi)情立即設(shè)臨時標(biāo)準地， 采取平行線抽樣法抽取20 株標(biāo)準株詳查。 幼蟲期調(diào)查， 1～2 齡幼蟲調(diào)查枯黃卷曲的枝數(shù)， 推算幼蟲數(shù)， 3 齡以上的幼蟲且3 m 以下的小樹直接調(diào)查合計樹冠上的幼蟲數(shù)， 大樹用“蟲糞粒推算法” 調(diào)查， 幼蟲越冬期間調(diào)查樹干基部樹皮縫中的幼蟲數(shù)推算全部蟲口。

1.2 災(zāi)變預(yù)測的建模方法

灰色系統(tǒng)理論認為， 一切隨機量都是在一定范圍內(nèi)、 一定時段上變化的灰色量及其灰色過程， 主張從事物內(nèi)部研究其發(fā)展變化規(guī)律。 對于灰色量的處理， 不是去尋求它的統(tǒng)計規(guī)律和概率分布， 而是從無規(guī)律的原始數(shù)據(jù)中找出規(guī)律， 即對數(shù)據(jù)通過一定方式處理后再建立模型。 因為客觀系統(tǒng)無論怎樣復(fù)雜，它總是有關(guān)聯(lián)、 有序、 有整體功能的， 作為系統(tǒng)行為特征的數(shù)據(jù)， 總是蘊含著某種規(guī)律。

如果在x（0）中有新的數(shù)列n′＜n是符合災(zāi)變定義的數(shù)列， 則對于上災(zāi)變其中：i′為災(zāi)變年序號，V′為災(zāi)變年序號取值范圍。

為了便于辨認， 特記災(zāi)變年序號i′為z（i′）或z（i）。 并作災(zāi)變號映射記z的1 次累加生成數(shù)（accumulated generating operation，則災(zāi)變預(yù)測的GM（1,1）模型有如下算式：其中：t為年份，a為發(fā)展灰數(shù)，u為內(nèi)生控制灰數(shù)，為待估灰數(shù)， B 為累加矩陣， BT為轉(zhuǎn)置矩陣，yN為常數(shù)項向量。

2 結(jié)果與分析

2.1 馬尾松毛蟲累計蟲口數(shù)

馬尾松毛蟲每年越冬代、 1 代幼蟲和2 代幼蟲的累計蟲口數(shù)見表1， 將其作為原始數(shù)據(jù)。

表1 馬尾松毛蟲越冬代、 1 代幼蟲和2 代幼蟲累計蟲口數(shù)Table 1 Accumulative population of D. punctatus during wintering, first and second generation larvae

2.2 馬尾松毛蟲越冬代災(zāi)變預(yù)測

2.2.1 確定災(zāi)變閾值 根據(jù)安徽省潛山縣森林病蟲害防治總站提供的1983-2016 年馬尾松毛蟲發(fā)生情況的資料， 馬尾松毛蟲越冬代蟲口數(shù)大于15 頭·株-1的年份為大發(fā)生年， 故以每株蟲口數(shù)15 頭為災(zāi)變閾值。

2.2.2 對原始數(shù)列作有關(guān)映射 原始數(shù)列為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}。 對原始數(shù)列作災(zāi)變映射， 按ξ≥15 得災(zāi)變數(shù)列：

然后再作災(zāi)變序號映射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āi)變?nèi)掌诩痾={1， 2， 3， 6， 7， 8， 9， 15， 22， 23}即為建模時所需的原始數(shù)列。

2.2.3 建立灰色預(yù)測GM（1,1）模型 對災(zāi)變年序集z 建立GM（1,1）模型， 按原始數(shù)列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）建模：由此得出馬尾松毛蟲越冬代蟲口數(shù)的GM（1,1）災(zāi)變預(yù)測模型：

將根據(jù)馬尾松毛蟲越冬代蟲口數(shù)的GM（1,1）災(zāi)變預(yù)測模型得到的擬合值、 誤差等列于表2。

表2 馬尾松毛蟲越冬代蟲口數(shù)預(yù)測的擬合值與誤差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 個均數(shù)之差）只有觀察值均數(shù)的3.40%。 實際觀察值與通過GM（1,1）災(zāi)變預(yù)測模型得到的擬合值之間誤差較小， 各災(zāi)變年的平均精度為84.40%， 總體精度（觀察值的平均值減去誤差的平均值與觀察值的平均值之比） 為96.25%。 由此可以通過這個模型求得未來5 個時刻的預(yù)測值， 用所建災(zāi)變預(yù)測模型z＾（1）（k+1）=9.580 75e0.26933k-8.580 75對未來馬尾松毛蟲越冬代災(zāi)變年份進行預(yù)測z＾（0）（11）=z＾（1）（11）-z（1）（10）=33.434 8， 由原始序號得知， 原始序列最后一次災(zāi)變是2011 年， z（0）（10）=23， 現(xiàn)預(yù)測下一次災(zāi)變年的序號z＾（0）（11）=33.434 8， 即從2011 年馬尾松毛蟲越冬代災(zāi)變年算起， 再過10 a 即2021 年為馬尾松毛蟲越冬代大發(fā)生年， 同理可預(yù)測之后年份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 代幼蟲災(zāi)變預(yù)測

與上述馬尾松毛蟲越冬代蟲口數(shù)GM（1,1）災(zāi)變預(yù)測模型相似， 由表1 根據(jù)閾值大于15 頭·株-1可以得到馬尾松毛蟲1 代幼蟲蟲口數(shù)GM（1,1）災(zāi)變預(yù)測模型。

原始數(shù)列為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}。

對原始數(shù)列作災(zāi)變映射， 按ξ≥15 得災(zāi)變數(shù)列， 最終得到災(zāi)變?nèi)掌诩?{1， 2， 6， 7， 8， 15， 17，22， 23}即為建模時所需的原始數(shù)列。

作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 代幼蟲蟲口數(shù)的GM（1,1）災(zāi)變預(yù)測模型：

將根據(jù)馬尾松毛蟲1 代幼蟲蟲口數(shù)的GM（1,1）災(zāi)變預(yù)測模型得到的擬合值、 誤差等列于表3。

表3 馬尾松毛蟲1 代幼蟲蟲口數(shù)預(yù)測的擬合值與誤差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 個均數(shù)之差）只有觀察值均數(shù)的6.49%。 實際觀察值與通過GM（1,1）災(zāi)變預(yù)測模型得到的擬合值之間誤差較小， 各災(zāi)變年的精度平均為84.85%， 總體精度為92.34%。 由此可以通過這個模型求得未來5 個時刻的預(yù)測值， 用所建災(zāi)變預(yù)測模型z＾（1）（k+1）=18.181 8e0.24187k-17.181 8 對未來馬尾松毛蟲一代幼蟲災(zāi)變年份進行預(yù)測z＾（0）（10）=z＾（1）（10）-z（1）（9）=34.443 8， 由原始序號得知， 原始序列最后一次災(zāi)變是2011年，z（0）（9）=23， 現(xiàn)預(yù)測下一次災(zāi)變年的序號z＾（0）（10）=34.443 8， 即從2011 年馬尾松毛蟲1 代幼蟲災(zāi)變年算起， 再過11 a 即2022 年為馬尾松毛蟲一代幼蟲大發(fā)生年， 同理可以預(yù)測之后年份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 代幼蟲災(zāi)變預(yù)測

與馬尾松毛蟲越冬代和1 代幼蟲蟲口數(shù)GM（1,1）災(zāi)變預(yù)測模型相似， 由表1 根據(jù)閾值大于15 頭·株-1可以得到馬尾松毛蟲2 代幼蟲蟲口數(shù)GM（1,1）災(zāi)變預(yù)測模型。

原始數(shù)列為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}。

對原始數(shù)列作災(zāi)變映射， 按ξ≥15 得災(zāi)變數(shù)列， 最終得到災(zāi)變?nèi)掌诩痾={1， 2， 5， 6， 7， 8， 14，15， 16， 22， 23}即為建模時所需的原始數(shù)列。

作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 代幼蟲蟲口數(shù)的GM（1,1）災(zāi)變預(yù)測模型：

將根據(jù)馬尾松毛蟲2 代幼蟲蟲口數(shù)的GM（1,1）災(zāi)變預(yù)測模型得到的擬合值、 誤差等列于表4。

對觀察值和擬合值的差異進行t檢驗，t=0.195 2，df=20 時，t0.05=2.09，t＜t0.05， 兩者差異不顯著， 平均誤差只占觀察值均數(shù)的5.4%。 實際觀察值與通過GM（1,1）災(zāi)變預(yù)測模型得到的擬合值之間誤差較小，各災(zāi)變年的精度平均為84.08%， 總體平均精度為94.09%。 由此可以通過這個模型求得未來5 個時刻的預(yù)測值， 用所建災(zāi)變預(yù)測模型對未來馬尾松毛蟲2 代幼蟲災(zāi)變年份進行預(yù)測z＾（0）（12）=z＾（1）（12）-z（1）（11）=31.704 2， 由原始序號得知， 最后一次災(zāi)變是2011 年，z（0）（11）=23， 現(xiàn)預(yù)測下一次災(zāi)變年的序號z＾（0）（12）=31.704 2， 即從2011 年馬尾松毛蟲2 代幼蟲災(zāi)變年算起， 再過9 a 即2020 年為馬尾松毛蟲2 代幼蟲大發(fā)生年， 同理可預(yù)測之后年份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 代幼蟲蟲口數(shù)預(yù)測的擬合值與誤差Table 4 Fitting value and error of prediction of population number of second generation larvae of D. punctatus

3 小結(jié)與討論

利用灰色災(zāi)變模型預(yù)測法對安徽省潛山縣1989-2016 年馬尾松毛蟲幼蟲越冬代、 1 代和2 代累計蟲口的災(zāi)變進行預(yù)測， GM（1,1）災(zāi)變預(yù)測模型分別是e0.24187k-17.181 8 和根據(jù)預(yù)測模型求得的已知年份的擬合值與觀察值進行t檢驗，t=0.101 7， 0.218 6 和0.195 2， 均小于t0.05， 差異均不顯著， 總體平均精度達96.25%、92.34%和94.09%， 各災(zāi)變年精度平均值為84.40%、 84.85%和84.08%， 誤差極小。 說明災(zāi)變預(yù)測法對馬尾松毛蟲幼蟲發(fā)生量災(zāi)變的預(yù)測預(yù)報是一種較理想的預(yù)報方法。 從2011 年算起， 預(yù)測越冬代災(zāi)變年為2021 年， 1 代災(zāi)變年為2022 年， 2 代災(zāi)變年為2020 年。

害蟲種群的消長， 其本質(zhì)是害蟲種群與環(huán)境之間進行物質(zhì)交換和能量傳遞的結(jié)果， 它既含有已知的信息， 又含有未知的或非確知的信息， 實際上是一個灰色系統(tǒng)。 灰色系統(tǒng)災(zāi)變預(yù)測一般是指某種大于閾值的災(zāi)害在哪一年份出現(xiàn)， 灰色系統(tǒng)災(zāi)變預(yù)測已經(jīng)運用到國民經(jīng)濟各個方面的預(yù)測， 在植物保護方面也被大量利用。 在農(nóng)業(yè)害蟲的預(yù)測上， 預(yù)測結(jié)果與實際情況吻合度很高［11-13］。 灰色系統(tǒng)災(zāi)變預(yù)測結(jié)果在林業(yè)害蟲的預(yù)測上也有報道， 吳秋芳等［14］對河南省安陽縣楊樹Populus苗圃地下害蟲進行灰色災(zāi)變預(yù)測，平均精度達90.40%。 宋丁權(quán)等［15］用灰色理論中的災(zāi)變模型以及灰色區(qū)間模型對馬尾松毛蟲災(zāi)級出現(xiàn)的年份進行了模擬、 預(yù)測和分析， 也取得了較好的效果。 與中國馬尾松毛蟲發(fā)生的其他省區(qū)相比， 安徽省馬尾松毛蟲發(fā)生的偶災(zāi)區(qū)， 由于安徽省馬尾松林面積較大， 加之還有一定面積的濕地松和火炬松的純林和混交林， 為馬尾松毛蟲的猖獗提供了發(fā)生的基礎(chǔ)條件， 安徽省地處亞熱帶和溫帶的過渡帶， 溫、 濕、降雨等氣候條件均適宜于馬尾松毛蟲的發(fā)生， 馬尾松毛蟲發(fā)生的森林環(huán)境， 常災(zāi)區(qū)絕大部分為海拔100～300 m 的丘陵地區(qū)， 偶災(zāi)區(qū)為海拔300～500 m 的山區(qū)， 海拔500 m 以上的山區(qū)很少有馬尾松毛蟲發(fā)生， 安徽省為海拔400～800 m 的偶災(zāi)區(qū)［1］。 馬尾松毛蟲在遇到適宜的發(fā)生條件， 極易爆發(fā)成災(zāi)， 所以對馬尾松毛蟲猖獗成災(zāi)的預(yù)報有著重要的意義。 本研究預(yù)測了安徽省潛山縣馬尾松毛蟲越冬代、 1 代和2代幼蟲發(fā)生的嚴重年份， 平均精度分別為96.25%、 92.34%和94.09%。 灰色災(zāi)變模型預(yù)測法預(yù)測周期長， 預(yù)測結(jié)果精度高， 是較理想的預(yù)測方法。 對灰色災(zāi)變模型的應(yīng)用在原始序列數(shù)據(jù)中應(yīng)有一定的時間長度， 這樣精度才能提高， 再者災(zāi)變模型隨著時間的推移， 也需要不斷更新和修正。