基于STARMA模型的感病松林內(nèi)松墨天牛成蟲發(fā)生量預(yù)測

周書永，丁昌地，陳繪畫

(1.浙江省臨海市森林病蟲防治檢疫站，浙江 臨海 317000；2.浙江省臨海市桃渚鎮(zhèn)農(nóng)業(yè)綜合服務(wù)中心，浙江 臨海 317000；3.浙江省仙居縣林業(yè)局，浙江 仙居 317300)

1 引言

松墨天牛(MonochamusalternatusHope)成蟲取食馬尾松(Pinusmasoniana)、黑松(P.thunbergii)等樹種1～2年生的嫩枝作為補(bǔ)充營養(yǎng)，導(dǎo)致寄主樹勢衰弱進(jìn)而引發(fā)松墨天牛成蟲聚集產(chǎn)卵、過多的松墨天牛幼蟲鉆蛀導(dǎo)致松樹寄主死亡[1]，松墨天牛還是我國松材線蟲病Bursaphelenchusxylophilus(Steiner&Buhere,1934)Nickle的主要傳播媒介。松墨天牛一生的大部分時(shí)間都生活于寄主體內(nèi)，只在成蟲時(shí)期離開寄主，啃食寄主嫩枝用以補(bǔ)充營養(yǎng)。松墨天牛成蟲的林間數(shù)量動(dòng)態(tài)變化多有研究[2～5]，陳順立等還用神經(jīng)網(wǎng)絡(luò)方法預(yù)測松墨天牛的發(fā)生量[6]，但有關(guān)感病馬尾松林中松墨天牛成蟲的發(fā)生量尚未見報(bào)道。

時(shí)空序列是指在空間上有相關(guān)關(guān)系的多個(gè)時(shí)間序列的集合。時(shí)空序列建模則是指尋找一種分析時(shí)空序列數(shù)據(jù)的方法，對未觀測時(shí)空位置的屬性進(jìn)行預(yù)測。實(shí)現(xiàn)時(shí)空序列模型的方法有時(shí)空動(dòng)力學(xué)方法、時(shí)空相關(guān)統(tǒng)計(jì)學(xué)方法、時(shí)空地統(tǒng)計(jì)學(xué)方法等[7]，時(shí)空自相關(guān)移動(dòng)平均模型(Spatial Temporal Auto Regressive and Moving Average.STARMA)于20世紀(jì)80年代提出，它假設(shè)變量是當(dāng)前位置以及其空間鄰接位置過去幾個(gè)時(shí)期觀測值與隨機(jī)誤差項(xiàng)的線性組合，根據(jù)歷史的觀測數(shù)據(jù)建立時(shí)空相關(guān)-自相關(guān)函數(shù)得到時(shí)空相關(guān)性模型[8,9]，STARMA模型已在交通、環(huán)境、經(jīng)濟(jì)等領(lǐng)域都有成功的應(yīng)用[10,11]。本文基于3個(gè)感病松林試驗(yàn)點(diǎn)內(nèi)松墨天牛成蟲的連續(xù)時(shí)序監(jiān)測數(shù)據(jù)開展松墨天牛成蟲發(fā)生量的預(yù)測研究，根據(jù)松墨天牛成蟲誘捕數(shù)量的時(shí)間自相關(guān)性、空間相關(guān)性以及時(shí)間空間的互相關(guān)性，建立松墨天牛發(fā)生量的時(shí)空自相關(guān)移動(dòng)平均模型，對松墨天牛發(fā)生量進(jìn)行預(yù)測。

2 材料與方法

2.1 試驗(yàn)點(diǎn)概況

臨海市位于浙江省東南部，地處28°40′37″～29°04′15″N、120°49′41″～121°42′25″E，東西長為85 km，南北寬為44 km，全市總面積2208.5 km2，其中松林面積為5.9×104hm2，以馬尾松純林為主，3個(gè)試驗(yàn)點(diǎn)的概況及懸掛的誘捕器數(shù)量詳見表1。

表1 3個(gè)試驗(yàn)點(diǎn)的概況及懸掛的誘捕器數(shù)量

2.2 試驗(yàn)材料及松墨天牛成蟲林間動(dòng)態(tài)觀測

松樹蛀干類害蟲引誘劑(規(guī)格：300 mL/瓶)由浙江省林業(yè)有害生物防治檢疫局和中國林業(yè)科學(xué)研究院亞熱帶林業(yè)研究所研制、寧波中化化學(xué)品有限公司生產(chǎn)(2010年后由浙江黃巖鼎正化工有限公司生產(chǎn))，以及與松墨天牛引誘劑配套的小型折疊式誘捕器(2007～2011年4～9月份使用)。松墨天牛性引誘劑F1由浙江農(nóng)林大學(xué)生物農(nóng)藥高效制備國家地方聯(lián)合工程實(shí)驗(yàn)室研制，以及杭州科森農(nóng)化有限公司與浙江農(nóng)林大學(xué)生物農(nóng)藥高效制備國家地方聯(lián)合工程實(shí)驗(yàn)室研制BF-1型誘捕器(2012～2015年4～9月份使用)。2007～2015年4～9月份，在臨海市江南街道、沿江鎮(zhèn)及大田街道的馬尾松人工林和天然林內(nèi)，將誘捕器架設(shè)在距地面1.5 m、兩株相近的馬尾松樹間鐵絲上。誘捕器的間距為80～130 m，單只離林緣40～70 m，3個(gè)試驗(yàn)點(diǎn)各架設(shè)20只誘捕器。在引誘劑瓶口的包裝材料(直徑3 cm)上打0.5 cm直徑小洞3個(gè)。每隔3～4周更換1次引誘劑，收蟲的間隔時(shí)間為2～7 d，以月為單位進(jìn)行統(tǒng)計(jì)。由于4月的誘捕量太少，故將4月的誘捕量合并到5月。

2.3 STARMA模型

STARMA模型為自回歸移動(dòng)平均(Autoregressive Moving Average，ARMA)模型在空間的擴(kuò)展，該模型在考慮某觀測值所在位置時(shí)間序列的同時(shí)，也考慮其空間相鄰位置的時(shí)間序列，從而提高時(shí)空序列預(yù)測的準(zhǔn)確度。

2.3.1 STARMA模型的描述

STARMA模型是針對同一個(gè)變量的空間相關(guān)性與時(shí)間相關(guān)性建立的，Zα(t)是由空間點(diǎn)的過去時(shí)刻值和過去隨機(jī)誤差值，再加上其空間鄰近點(diǎn)的過去時(shí)刻值和過去隨機(jī)誤差值影響，空間位置的相關(guān)性以空間權(quán)值矩陣的形式表示。STARMA模型為：

(1)

式(1)中：p為時(shí)間自回歸階數(shù)；q為時(shí)間移動(dòng)平均階數(shù)；λk為第k個(gè)自回歸項(xiàng)空間階數(shù)；mk為第k個(gè)移動(dòng)平均項(xiàng)空間階數(shù)；φkl為需要估計(jì)的時(shí)間滯后k、空間滯后的自回歸系數(shù)；θkl為需要估計(jì)的時(shí)間滯后k、空間滯后的移動(dòng)平均系數(shù)；Wl為l階空間鄰接的N×N權(quán)重矩陣，εα(t)是隨機(jī)變量。

2.3.2 STARMA模型的建立

STARMA模型的建模過程包括數(shù)據(jù)預(yù)處理、模型辨識(shí)、參數(shù)估計(jì)、模型診斷以及預(yù)測分析等5步。通過數(shù)據(jù)預(yù)處理，將非平穩(wěn)時(shí)空序列轉(zhuǎn)換為平穩(wěn)時(shí)空序列；模型辨識(shí)則是初步確認(rèn)候選模型；參數(shù)估計(jì)為求解候選模型的參數(shù)；模型診斷則是根據(jù)相關(guān)指標(biāo)判斷是否接受所建立的模型；模型預(yù)測是將建立的模型進(jìn)行預(yù)測并分析模型的預(yù)測性能。

設(shè)N是二維空間的監(jiān)測點(diǎn)數(shù)，T是觀測時(shí)間的最大值，則經(jīng)T次觀測，共有總數(shù)為NT的觀測值，全部觀測值的平均值為[12]：

mz=∑∑∑Zαt/NT

(2)

空間權(quán)重矩陣是二維空間位置上變量值受其空間鄰近位置變量值影響的定量化測度。若研究空間有n個(gè)面積單元，其中任何兩個(gè)面積單元間都存在一個(gè)空間關(guān)系，則有n×n對關(guān)系，因而可以用n×n的矩陣存儲(chǔ)這n個(gè)面積單元間的空間關(guān)系?？臻g權(quán)重矩陣W的具體形式為：

(3)

空間權(quán)重矩陣的建立方法有邊界法[13,14]、拓?fù)鋱D法[9]等，用邊界法建立的權(quán)重矩陣元素為：

(4)

(5)

二維空間某單元的變量值受其鄰近空間位置變量值的影響度稱空間延遲算子，用L表示，則空間延遲算子L與空間權(quán)重矩陣W的關(guān)系為：

L0zα=INzα,L1zα=W1zβ,…,Lhzα=Whzβ

(6)

式(6)中：I為單位矩陣，h為矩陣階數(shù)。

3 結(jié)果與分析

3.1 數(shù)據(jù)預(yù)處理

將3個(gè)感病松林試驗(yàn)點(diǎn)內(nèi)2007～2013年每年5～9月的松墨天牛成蟲連續(xù)誘捕數(shù)據(jù)按下式進(jìn)行數(shù)據(jù)預(yù)處理，去掉原始觀測數(shù)據(jù)中的線性趨勢和周期性趨勢：

Zα(t)=Yα(t)+Tα(t)+Cα(t)

(7)

式(7)中：Yα(t)為二階平穩(wěn)隨機(jī)函數(shù)，Tα(t)為線性趨勢，Cα(t)為周期性趨勢。線性趨勢和周期性趨勢的去除方法詳見文獻(xiàn)[15,16]。

3.2 STARMA模型的辨識(shí)及參數(shù)估計(jì)

在建立STARMA模型時(shí)，通常根據(jù)時(shí)空自相關(guān)和偏相關(guān)函數(shù)來確定模型的階數(shù)，3個(gè)試驗(yàn)點(diǎn)觀測值預(yù)處理后的時(shí)空自相關(guān)和偏相關(guān)函數(shù)見表2、3。

表2 3個(gè)試驗(yàn)點(diǎn)的時(shí)空自相關(guān)函數(shù)

表3 3個(gè)試驗(yàn)點(diǎn)的時(shí)空偏自相關(guān)函數(shù)

從表2、表3可以看出，無論是時(shí)空自相關(guān)函數(shù)還是時(shí)空偏自相關(guān)函數(shù)，江南街道和沿江鎮(zhèn)在時(shí)間上都延遲1期后呈截尾趨勢，大田街道則在時(shí)間上都延遲3期后呈截尾趨勢。根據(jù)引誘劑誘殺松墨天牛的實(shí)際，結(jié)合文獻(xiàn)13的研究，空間權(quán)值矩陣的階取2。因此可以判斷江南街道和沿江鎮(zhèn)的數(shù)據(jù)均為一個(gè)1期時(shí)間延遲和2階空間延遲的自相關(guān)過程，大田街道則是一個(gè)3期時(shí)間延遲和2階空間延遲的自相關(guān)過程。對數(shù)據(jù)預(yù)處理后的3個(gè)試驗(yàn)點(diǎn)內(nèi)2007～2013年每年5～9月的松墨天牛成蟲連續(xù)誘捕數(shù)據(jù)按式(1)進(jìn)行最小二乘估計(jì)，得到大田街道的方程表達(dá)式為：

z(t)=0.2120+0.0847z(t-1)-0.1185z(t-2)-0.0727z(t-3)+0.0428W1·z(t-1)-0.0784W1·z(t-2)-0.0576W1·z(t-3)+0.2339W2·z(t-1)-0.1456W2·z(t-2)+0.0842W2·z(t-3)

(8)

江南街道的表達(dá)式為：

z(t)=-0.0985+0.1184z(t-1)-0.0063W1·z(t-1)+0.2641W2·z(t-1)

(9)

沿江鎮(zhèn)的表達(dá)式為：

z(t)=-0.2221+0.1155z(t-1)+0.0325W1·z(t-1)+0.0407W2·z(t-1)

(10)

3.3 STARMA模型的檢驗(yàn)

對模型(8)～(10)的參數(shù)進(jìn)行顯著性檢驗(yàn)，刪除模型中參數(shù)檢驗(yàn)不顯著的自變量，再次進(jìn)行最小二乘估計(jì)和模型參數(shù)顯著性檢驗(yàn)，得到所有自變量都顯著的模型為[17]：

大田街道：

z(t)=0.2202-0.0676z(t-1)-0.1274z(t-2)+0.0910z(t-3)-0.1798W2·z(t-2)+0.2445W2·z(t-3)

(11)

江南街道：

z(t)=-0.0985+0.1169z(t-1)+0.2609W2·z(t-1)

(12)

沿江鎮(zhèn)：

z(t)=-0.2201+0.1462z(t-1)

(13)

對模型(11)～(13)進(jìn)行F檢驗(yàn)，檢驗(yàn)結(jié)果為：大田街道，F(xiàn)=14.5927>F0.01(5,500)=3.05；江南街道，F(xiàn)=25.1445>F0.01(2,500)=4.65；沿江鎮(zhèn)，F(xiàn)=16.4453>F0.01(1,500)=6.69，說明所設(shè)的模型(11)～(13)正確。模型(11)～(13)的殘差時(shí)空自相關(guān)函數(shù)和時(shí)空偏自相關(guān)函數(shù)計(jì)算結(jié)果列于表4、表5。

從表4、表5看出，模型(11)～(13)的殘差時(shí)空自相關(guān)函數(shù)值和時(shí)空偏自相關(guān)函數(shù)值都較小，在時(shí)間上或空間上都不存在顯著的自相關(guān)，說明所選擇的模型較好地表示了時(shí)空數(shù)據(jù)序列，因而所確定的模型和估計(jì)的參數(shù)可以預(yù)測各試驗(yàn)點(diǎn)內(nèi)松墨天牛的誘集數(shù)量。

3.4 STARMA模型預(yù)測結(jié)果的檢驗(yàn)

設(shè)所建模型某月松墨天牛誘捕量的預(yù)報(bào)值為X1j，實(shí)際誘捕量為X2j，則：

Dj=X1j-X2j,j=1,2,…,n

(14)

(15)

服從自由度為n-1的大分布，其中

因此，在顯著性水平為α?xí)r，可通過比較|t|與tn-1(α/2)，來檢驗(yàn)假設(shè)H0：δ=0和H1：δ≠0。即若|t|tn-1(α/2)，則預(yù)報(bào)值與實(shí)測值有顯著差異[18]。3個(gè)試驗(yàn)點(diǎn)未參與建模的2014～2015年每年5～9月松墨天牛成蟲的預(yù)測結(jié)果及檢驗(yàn)見表6。

表4 3個(gè)試驗(yàn)點(diǎn)殘差的時(shí)空自相關(guān)函數(shù)

表5 3個(gè)試驗(yàn)點(diǎn)殘差的時(shí)空偏自相關(guān)函數(shù)

表6 3個(gè)試驗(yàn)點(diǎn)未參與建模的2014～2015年每年5～9月松墨天牛成蟲的預(yù)測結(jié)果

從表6可以看出，所建立的模型預(yù)測3個(gè)試驗(yàn)點(diǎn)未參與建模的2014～2015年每年5～9月的松墨天牛成蟲，大田街道預(yù)測準(zhǔn)確的6次、不準(zhǔn)確的4次，預(yù)測成功率為60%；江南街道預(yù)測準(zhǔn)確的7次、不準(zhǔn)確的3次，預(yù)測成功率為70%；沿江鎮(zhèn)預(yù)測準(zhǔn)確的8次、不準(zhǔn)確的2次，預(yù)測成功率為80%，平均預(yù)測成功率為70%。

4 小結(jié)與討論

本文在考慮感病馬尾松林內(nèi)各松墨天牛誘捕點(diǎn)自身過去與現(xiàn)在、未來的時(shí)間關(guān)系的同時(shí)，還考慮周圍相鄰誘捕點(diǎn)的過去與該誘捕點(diǎn)未來誘捕量變化的空間關(guān)系，從整個(gè)感病松林出發(fā)，綜合評價(jià)網(wǎng)絡(luò)時(shí)空拓?fù)潢P(guān)系，將 STARMA 方法應(yīng)用到松墨天牛誘捕數(shù)量的預(yù)測領(lǐng)域中，并借助浙江省臨海市2007～2015年3個(gè)感病松林試驗(yàn)點(diǎn)來驗(yàn)證方法的可行性，預(yù)測3個(gè)試驗(yàn)點(diǎn)未參與建模的2014～2015年每年5～9月的松墨天牛成蟲誘集數(shù)量，平均預(yù)測成功率為70%。較高的預(yù)測成功率表明了此方法用于松墨天牛誘捕數(shù)量預(yù)測的可靠性，為松墨天牛的可持續(xù)控制提供科學(xué)的依據(jù)。

STARMA方法可以較好地將時(shí)間和空間屬性融入到模型中，對網(wǎng)絡(luò)條件下松墨天牛誘捕數(shù)量預(yù)測的適應(yīng)性較強(qiáng)，借助于STARMA 方法進(jìn)行預(yù)測，其平均預(yù)測效果較好，但是STARMA 方法還有待改進(jìn)，如在確定空間權(quán)重矩陣過程中，只考慮了誘捕點(diǎn)間的相鄰與否，沒有考慮相鄰誘捕點(diǎn)間距離的遠(yuǎn)近。在今后的研究中，可以在模型中加入空間距離對誘捕效果的影響，保證模型參數(shù)盡可能地反映實(shí)際情況，提高模型預(yù)測的成功率。