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

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

2 材料與方法

2.1 試驗點概況

臨海市位于浙江省東南部，地處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個試驗點的概況及懸掛的誘捕器數(shù)量詳見表1。

表1 3個試驗點的概況及懸掛的誘捕器數(shù)量

2.2 試驗材料及松墨天牛成蟲林間動態(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)合工程實驗室研制，以及杭州科森農(nóng)化有限公司與浙江農(nóng)林大學(xué)生物農(nóng)藥高效制備國家地方聯(liá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個試驗點各架設(shè)20只誘捕器。在引誘劑瓶口的包裝材料(直徑3 cm)上打0.5 cm直徑小洞3個。每隔3～4周更換1次引誘劑，收蟲的間隔時間為2～7 d，以月為單位進行統(tǒng)計。由于4月的誘捕量太少，故將4月的誘捕量合并到5月。

2.3 STARMA模型

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

2.3.1 STARMA模型的描述

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

(1)

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

2.3.2 STARMA模型的建立

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

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

mz=∑∑∑Zαt/NT

(2)

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

(3)

空間權(quán)重矩陣的建立方法有邊界法[13,14]、拓撲圖法[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個感病松林試驗點內(nèi)2007～2013年每年5～9月的松墨天牛成蟲連續(xù)誘捕數(shù)據(jù)按下式進行數(shù)據(jù)預(yù)處理，去掉原始觀測數(shù)據(jù)中的線性趨勢和周期性趨勢：

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

(7)

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

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

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

表2 3個試驗點的時空自相關(guān)函數(shù)

表3 3個試驗點的時空偏自相關(guān)函數(shù)

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

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)

江南街道的表達式為：

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

(9)

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

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

(10)

3.3 STARMA模型的檢驗

對模型(8)～(10)的參數(shù)進行顯著性檢驗，刪除模型中參數(shù)檢驗不顯著的自變量，再次進行最小二乘估計和模型參數(shù)顯著性檢驗，得到所有自變量都顯著的模型為[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)進行F檢驗，檢驗結(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)的殘差時空自相關(guān)函數(shù)和時空偏自相關(guān)函數(shù)計算結(jié)果列于表4、表5。

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

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

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

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

(14)

(15)

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

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

表4 3個試驗點殘差的時空自相關(guān)函數(shù)

表5 3個試驗點殘差的時空偏自相關(guān)函數(shù)

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

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

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