基于STWR模型的森林病蟲影響因素研究

蘇少強，闕 翔，嚴(yán)宣輝，何中聲，劉金福，李夢航，黃朝法，劉 海

(1 福建農(nóng)林大學(xué) a 計算機與信息學(xué)院，b福建省統(tǒng)計信息研究中心，c 生態(tài)與資源統(tǒng)計福建省高校重點實驗室，福建 福州 350002； 2 福建師范大學(xué) 數(shù)字福建環(huán)境監(jiān)測物聯(lián)網(wǎng)實驗室，福建 福州 350117； 3 福建省林業(yè)調(diào)查規(guī)劃院，福建 福州 350003)

森林病蟲害已成為林木生長和固碳的最大威脅之一[1]， 2006-2010年，我國年均森林病蟲害發(fā)生面積1 151.8萬hm2，造成林業(yè)年均直接經(jīng)濟(jì)損失約145億元[2]。加強森林病蟲害防治工作，減少我國森林病蟲害損失的關(guān)鍵在于厘清影響森林病蟲害發(fā)生關(guān)鍵因素及其時空變化特征。

氣候變化、人類活動能夠影響森林病蟲害發(fā)生率及其空間分布與演變。其中氣候變化是森林病蟲害規(guī)?；把杆贁U散的重要因素[3]，如長期干旱或澇災(zāi)會降低樹木的抗逆性，進(jìn)而導(dǎo)致森林更易遭受病蟲害影響[4]；暖冬可能加重了楊干象、美國白蛾、松縱坑切梢小蠹等蟲害的發(fā)生[5]；降水量是導(dǎo)致馬尾松毛蟲發(fā)生的主因之一[6]；氣溫升高和降水增多易使楊樹潰瘍病更嚴(yán)重[7]。另外，人類活動也會影響森林環(huán)境，如改變?nèi)斯ち峙c天然林的構(gòu)成[8]，進(jìn)而影響森林受病蟲害攻擊的可能性及存活率[9]，加劇松材線蟲、舞毒蛾等的擴散[10-11]。

國內(nèi)外學(xué)者對森林病蟲害的驅(qū)動因素開展了廣泛研究，取得了豐碩的成果。如Marini等[12]采用線性效應(yīng)模型分析得出小蠹蟲的局部暴發(fā)主要受夏季降雨不足和氣溫升高主導(dǎo)；張鵬霞等[13]以溫濕系數(shù)滑動均值、降水量變率等氣候因子分析得出，氣候變暖與環(huán)境干旱化驅(qū)動江西省森林病蟲害發(fā)生。Bao等[14]通過松毛蟲發(fā)生與不同尺度標(biāo)準(zhǔn)降水指數(shù)的相關(guān)性研究表明，長期干旱比短期干旱對松毛蟲暴發(fā)的影響更大。此外，氣候變化與病蟲害暴發(fā)周期也具有一定的同步性[15]。然而，森林病蟲害的發(fā)生及規(guī)模具有一定的復(fù)雜性和不確定性，多數(shù)研究仍以全局的線性模型(如OLS)量化氣候環(huán)境等因素的影響程度[16-17]，忽略了對森林病蟲害防治研究具有重要價值的空間異質(zhì)性和時變信息。地理加權(quán)回歸模型(GWR)引入空間位置信息，能很好地解決空間異質(zhì)性的問題，但對森林病蟲害發(fā)生具有的時間非平穩(wěn)性尚缺少考慮。

為更好地處理森林病蟲害分析過程中空間異質(zhì)性和時間異質(zhì)性問題，本研究采用基于數(shù)值差異變化率衰減賦權(quán)的時空加權(quán)回歸(STWR)模型，通過夏季均溫、冬季均溫、月均降水量和夜間燈光均值等變量，構(gòu)建各變量對福建省森林病蟲害發(fā)生影響的非平穩(wěn)性關(guān)系，揭示氣候變化、人類活動等因素對森林病蟲害發(fā)生的影響，為森林病蟲害的綜合防治提供理論支持。

1 材料與方法

1.1 研究區(qū)概況

福建省為全國六大林區(qū)省份之一，氣候?qū)賮啛釒ШＱ笮约撅L(fēng)氣候，境內(nèi)丘陵連綿，河谷、盆地穿插其間，山地、丘陵占全省總面積的80%以上，森林覆蓋率為66.8%，常見森林病蟲害種類為馬尾松毛蟲、剛竹毒蛾、竹蝗、毛竹枯梢病等[19]。

1.2 數(shù)據(jù)來源

根據(jù)相關(guān)研究，人類活動與森林病蟲害的擴散程度密切相關(guān)[11]，而夜間燈光是反映人類活動強度的重要指標(biāo)[18]，可考慮將夜間燈光均值作為森林病蟲害發(fā)生的影響因素。氣候變化也是劃分森林病蟲害分布的關(guān)鍵因素[24]，夏季均溫、冬季均溫與森林病蟲害發(fā)生的線性關(guān)系極為顯著[13],月均降水量越低的年份蟲害發(fā)生面積越大[25]。同時基于數(shù)據(jù)的可獲得性原則，本研究以森林病蟲害發(fā)生率為因變量，選取夜間燈光均值、夏季均溫、冬季均溫、月均降水量作為自變量，采集2008-2018年的數(shù)據(jù)。福建省林地面積數(shù)據(jù)來源于中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心(http：//www.resdc.cn/)，蟲害統(tǒng)計數(shù)據(jù)來源于福建省林業(yè)廳，氣象數(shù)據(jù)來源于WorldClim(https：//www.worldclim.org/)，兩種夜間燈光數(shù)據(jù)DMSP/OLS和NPP/VIIRS來源于美國海洋大氣局(NOAA)(https：//www.ngdc.noaa.gov)。在構(gòu)建STWR模型前，針對不同年份的夜間燈光數(shù)據(jù)進(jìn)行Z-Score標(biāo)準(zhǔn)化處理，剔除數(shù)據(jù)缺失縣域。所有數(shù)據(jù)均采用最大最小規(guī)范化處理，相關(guān)變量描述如表1。

表1 森林病蟲害發(fā)生率及影響因素描述Table 1 Description of forest diseases and insect pests and their driving factors

1.3 森林病蟲害時空變化特征分析

1.3.1 Mann-Kendall趨勢檢驗 Mann-Kendall(M-K)趨勢檢驗法假設(shè)時間序列數(shù)據(jù)Yn是n個彼此獨立、隨機分布的序列，原假設(shè)H0為序列無趨勢；對于所有的l,k≤n且yl和yk的分布是不相同的，其特點是數(shù)據(jù)不受少數(shù)異常值的干擾，不需服從某些特定分布，檢驗統(tǒng)計量為[20]：

(1)

(2)

式中：yl為某縣域第l年的森林病蟲害發(fā)生率，yk為該縣域第k年的森林病蟲害發(fā)生率，n為年份數(shù)。

當(dāng)n≥8時，統(tǒng)計量S大致服從正態(tài)分布，其均值為0，方差為：

(3)

標(biāo)準(zhǔn)化正態(tài)統(tǒng)計變量z通過下式計算：

(4)

在雙側(cè)趨勢檢驗中，在給定的置信水平0.05上，若|z|≥1.64，則時間序列數(shù)據(jù)變化趨勢具有顯著上升或下降的趨勢。其中若z>0，則序列呈上升趨勢；z<0，則序列呈下降趨勢。

1.3.2 空間聯(lián)系局部指標(biāo) 空間聯(lián)系局部指標(biāo)LISA揭示局部區(qū)域與相鄰區(qū)域的自相關(guān)性，其本質(zhì)是將全局Moran’sI分解到各區(qū)域單元[21]。LISA主要用來分析縣域i和j之間的空間自相關(guān)性，計算公式如下：

(5)

式中：yi、yj分別為i和j縣域的森林病蟲害發(fā)生率；n是縣域數(shù)；m是縣域i的鄰近縣域個數(shù)；u是所有觀測值的平均值；Pij為i和j之間的空間權(quán)重矩陣，若i和j相鄰，則值為1，不相鄰則值為0。進(jìn)一步采用Z檢驗對Ii進(jìn)行統(tǒng)計檢驗，判斷其顯著性。通過Ii來區(qū)分局部集聚模式，包括：“高-高”集聚是指Ii顯著為正且yi-u>0,表示森林病蟲害發(fā)生率高的縣域被森林病蟲害發(fā)生率高的縣域包圍；“低-低”集聚是指Ii顯著為正且yi-u<0,表示森林病蟲害發(fā)生率低的縣域被森林病蟲害發(fā)生率低的縣域包圍；“低-高”集聚是指Ii顯著為負(fù)且yi-u<0,表示森林病蟲害發(fā)生率低的縣域被森林病蟲害發(fā)生率高的縣域包圍；“高-低”集聚是指Ii顯著為負(fù)且yi-u>0,表示森林病蟲害發(fā)生率高的縣域被森林病蟲害發(fā)生率低的縣域包圍。

1.4 森林病蟲害影響因素的時空異質(zhì)性分析

回歸分析是探究自變量與因變量之間相互依賴關(guān)系的分析方法，普通線性回歸分析是其中的一種。因森林病蟲害發(fā)生與氣候變化、人類活動強度密切相關(guān)，因此采用普通線性回歸(OLS)模型可定量分析病蟲害發(fā)生率與其影響因素的數(shù)量關(guān)系。OLS模型的計算公式如下：

(6)

式中：xik和y分別是第i個觀察點第k個自變量和第i個觀察點的因變量，β0和βi分別代表著截距和系數(shù)，k表示自變量的個數(shù)，ε是回歸方程的隨機誤差項。

地理加權(quán)回歸(GWR)模型是對普通線性回歸型(OLS)模型的擴展，其將地理位置引入回歸參數(shù)中，允許回歸參數(shù)隨空間位置變化而改變，GWR的一般計算公式為：

(7)

式中：β0(ui,vi)為(ui,vi)的常數(shù)項，yi、xik分別是在(ui,vi)上第i個觀察點的因變量和在(ui,vi)上第i個觀察點第k個自變量，βk(ui,vi)為(ui,vi)上的回歸系數(shù)，εi為誤差項。GWR允許系數(shù)βk(ui,vi)在空間上變化來度量變量之間的非平穩(wěn)性關(guān)系，其系數(shù)的估計可表示為：

(8)

式中：W(ui,vi)是n×n的對角距陣，其對角線上的Wij表示第j個觀察點對第i個回歸點的影響權(quán)重，可通過空間核函數(shù)(如Gaussian和Bi-Square核)根據(jù)距離計算得到。

時空加權(quán)回歸(STWR)[22]模型是在地理加權(quán)回歸(GWR)框架[23]的基礎(chǔ)上引入時變關(guān)鍵信息，其基本形式為：

(9)

STWR模型采用一種基于時間變化的數(shù)值差異率衰減賦權(quán)策略，通過構(gòu)建加權(quán)平均的時空核函數(shù)來計算時空權(quán)重矩陣，其計算公式為：

(10)

2 結(jié)果與分析

2.1 福建省森林病蟲害時空變化特征分析

由圖1可知，福建省縣域森林病蟲害發(fā)病率較高的地區(qū)主要分布于東部沿海，如福州市、莆田市、泉州市、廈門市、三明市等，龍巖市先增后減的趨勢較為明顯。

采用空間聯(lián)系局部指標(biāo)和M-K趨勢檢驗分析森林病蟲害發(fā)生率的時空變化特征。依據(jù)福建省各縣域空間聯(lián)系局部指標(biāo)的計算結(jié)果，福建省縣域森林病蟲害發(fā)病率存在明顯的集聚特征，主要為“高-高”集聚和“低-低”集聚，且空間集聚特征隨時間變化相對穩(wěn)定(圖2)。2008-2014年，“低-低”集聚主要分布在龍巖市和漳州市北部的南靖縣、華安縣，集聚中心逐步轉(zhuǎn)移至福州市，但覆蓋范圍明顯減少；“高-高”集聚主要集中在泉州市，這與縣域森林病蟲害發(fā)病率的分布格局基本一致，且呈一定空間溢出效應(yīng)，在莆田市城廂區(qū)、荔城區(qū)有所增加。

根據(jù)縣域森林病蟲害發(fā)病率的M-K趨勢檢驗結(jié)果(圖3-a)，2008-2018年福建省有16個縣域森林病蟲害發(fā)生率顯著下降，12個縣域顯著上升。其中莆田市、泉州市、漳州市存在顯著上升趨勢，南平市、三明市有顯著下降趨勢，說明福建省縣域森林病蟲害發(fā)生率存在時間異質(zhì)性。結(jié)合福建省森林覆蓋情況(圖3-b)，森林面積處于中等水平(286～1 861 km2)的縣域，如泉州市安溪縣、永春縣和南平市延平區(qū)、政和縣等，其森林病蟲害發(fā)生率存在顯著上升趨勢。綜上，福建省縣域森林病蟲害發(fā)生率變化過程具有時空異質(zhì)性。

2.2 森林病蟲害影響因素的時空異質(zhì)性分析

考慮到縣域森林病蟲害發(fā)病率具有時空異質(zhì)性，森林病蟲害影響因素的時空異質(zhì)性分析采用OLS、GWR(空間核函數(shù)為Adaptive Bi-square)和STWR(模型的評價準(zhǔn)則為CV交叉驗證，時空核函數(shù)為Bi-square，帶寬搜索方式為Interval)3種模型進(jìn)行建模，并對其性能進(jìn)行對比。模型的因變量為森林病蟲害發(fā)生率，自變量為夜間燈光均值、夏季均溫、冬季均溫、月均降水量，時間范圍為2008-2018年。鑒于年份較多，本研究選取2008，2010，2012，2014，2016，2018年數(shù)據(jù)進(jìn)行建模分析，結(jié)果如表2所示。

表2 OLS、GWR、STWR 3種模型的擬合性能比較Table 2 Lomparison of fitting performance of OLS,GWR and STWR models

由表2可知，2008年GWR、STWR的R2分別為0.480和0.497，Sigma分別為0.100和0.099，SSE分別為0.637和0.617，表明GWR、STWR模型的擬合效果較為接近，可能是由于STWR模型的擬合期數(shù)較少。從2010年開始，STWR模型的擬合效果顯著提升，2010-2018年，STWR的R2分別比GWR分別提高了0.244，0.434，0.466，0.392和0.549，AICc值分別比GWR降低了34.576，65.933，104.499，65.005和121.572，表明STWR的擬合效果優(yōu)于OLS和GWR模型。其中，2012和2016年STWR的擬合效果略有下降，可能由于該年份數(shù)據(jù)量較少。

采用留一法交叉驗證(LOOCV)比較2018年的STWR和GWR預(yù)測的誤差，結(jié)果(圖4)表明，STWR預(yù)測的誤差比GWR小的點數(shù)為48個(圖4-a大部分紅色點都在藍(lán)色點的下方)，占比64%，說明STWR預(yù)測接近真實值的概率比GWR更高。GWR和STWR的預(yù)測平均誤差、預(yù)測均方誤差分別為0.052 66，0.007 36和0.039 66，0.004 70，STWR的平均誤差、均方誤差均更小，預(yù)測更加精確。對比箱線圖(圖4-b)，STWR預(yù)測誤差的集中程度更高，優(yōu)于GWR。

綜上，STWR擬合及預(yù)測性能優(yōu)于GWR和OLS。為了進(jìn)一步分析各自變量對森林病蟲害發(fā)生率的影響程度，將夜間燈光均值、夏季均溫、冬季均溫、月均降水量的回歸系數(shù)展示在地圖上，根據(jù)地圖分析各自變量對森林病蟲害發(fā)生率(因變量)影響的時空變化過程。

2.2.1 夜間燈光均值對森林病蟲害發(fā)生的影響 由圖5可知，夜間燈光均值對森林病蟲害發(fā)生起正向促進(jìn)作用，正向影響的范圍隨時間推移有一定程度縮小。2008-2018年夜間燈光均值呈正向作用的縣域總數(shù)大于呈負(fù)向作用的縣域總數(shù)，2008年正向縣域數(shù)占89.33%，2014年、2018年呈正向作用的縣域數(shù)(42)及負(fù)向作用的縣域數(shù)(33)較為接近，表明人類活動可能在早期城市建設(shè)中易造成森林病蟲害暴發(fā)，隨政府加大對保護(hù)森林資源的重視程度，其對森林病蟲害的發(fā)生表現(xiàn)出抑制作用，且該作用有增長的趨勢?？臻g上，2014-2018年夜間燈光均值回歸系數(shù)的正向影響主要集中在南平市和寧德市，其原因可能是該區(qū)域近年來受制于區(qū)位條件的限制，在建設(shè)新型城鎮(zhèn)化的同時，帶來了環(huán)境污染和生態(tài)承載飽和等問題[26]。

2.2.2 夏季均溫對森林病蟲害發(fā)生的影響 夏季均溫對森林病蟲害存在正向促進(jìn)作用(圖6)。2008-2018年夏季均溫正向影響縣域數(shù)遠(yuǎn)多于負(fù)向影響縣域數(shù)。從系數(shù)值上看，夏季均溫回歸系數(shù)最高值(3.420)大于冬季均溫回歸系數(shù)最大值(1.422)，反映夏季均溫的正向影響強度更大；2008-2018年夏季均溫正向影響的縣域數(shù)分別為48，45，51，39，42，34，而冬季均溫正向影響縣域數(shù)為28，27，34，49，50，35，表明2008-2014年夏季均溫正向影響比冬季均溫分布廣泛；而2016-2018年呈相反現(xiàn)象。值得關(guān)注的是2010年，負(fù)向影響在新羅區(qū)、漳浦縣較為顯著，可能源于新羅區(qū)和漳浦縣當(dāng)年夏季均溫分別為24.736，24.958 ℃，低于2009年(25.057，25.304 ℃)和2011年(25.156，25.410 ℃)。溫度升高縮短昆蟲世代時間，7月高溫引起馬尾松毛蟲發(fā)育期提前和世代數(shù)增加，夏季高溫易促使森林蟲害有更高的概率暴發(fā)[27]。

2.2.3 冬季均溫對森林病蟲害發(fā)生的影響 總體上，冬季均溫對福建省森林病蟲害發(fā)生率的正向作用逐漸增強，而負(fù)向影響作用先減少后增加(圖7)。2008年負(fù)向影響縣域數(shù)超過了正向影響縣域數(shù)，隨時間增長，2010和2012年正向與負(fù)向影響區(qū)域基本持平，2014和2016年正向影響縣域數(shù)分別大約占據(jù)福建省的2/3及3/4，2018年負(fù)向影響縣域數(shù)(40)多于正向影響(35)。2014年負(fù)向影響主要集中在寧德市蕉城區(qū)、福安市、霞浦縣，而3個縣域當(dāng)年的冬季均溫分別為24.863，25.159，25.636 ℃，低于2013年(25.974，26.804，26.710 ℃)和2015年(24.970，25.815，25.797 ℃)，說明較低的冬季溫度可能使害蟲發(fā)育緩慢、新陳代謝下降、產(chǎn)卵量降低，甚至死亡[28]。

2.2.4 月均降水量對森林病蟲害發(fā)生的影響 月均降水量對森林病蟲害主要為負(fù)向抑制作用(圖8)。STWR擬合結(jié)果顯示,2008，2010，2016和2018年月均降水量的回歸系數(shù)均值為-0.136 8，且負(fù)向影響縣域數(shù)量均多于正向影響縣域數(shù)，2012和2014年兩者數(shù)量基本一致，表明干旱的環(huán)境可能加劇森林蟲害發(fā)生。2018年正向影響在龍巖市上杭縣、武平縣、連城縣極為顯著，3個縣域當(dāng)年的月均降水量分別為134.206,139.605,97.398 mm，多于2017年的119.155,123.149,97.398 mm，表明長時間的強度降水會降低害蟲的種群密度，從而減少蟲害發(fā)生面積[29]。

2.2.5 森林病蟲害高發(fā)生率縣域的影響因素分析 為探究福建省森林病蟲害高發(fā)生率縣域影響因素的變化特征，選取2008-2018年森林病蟲害發(fā)生率最高的5個縣域，繪制各變量影響系數(shù)隨時間變化的箱線圖(圖9)。

夜間燈光均值為最大的負(fù)向影響因子，月均降水量的負(fù)向作用較為顯著。2010、2018年夜間燈光均值的回歸系數(shù)總體分布為-0.1～0，最小值為-4.66，夜間燈光均值負(fù)向影響較為突出的原因可能是政府不斷加強對科學(xué)防治森林病蟲害的重視程度，人類活動的抑制作用日益突出。2008-2018年中，4年月均降水量的回歸系數(shù)為負(fù)值，主要分布在-0.4～0，這與降水通常抑制森林病蟲害發(fā)生的結(jié)論相符。夏季均溫、冬季均溫的回歸系數(shù)在零值附近分布較為均勻，但夏季均溫的回歸系數(shù)大于冬季均溫，表明夏季均溫的影響強于冬季均溫。

3 結(jié) 論

本研究通過M-K趨勢檢驗及空間聯(lián)系局部指標(biāo)(LISA)探究福建省縣域森林病蟲害發(fā)病率的時空分布格局，發(fā)現(xiàn)其在時空上存在強烈的局部非平穩(wěn)關(guān)系：在石獅市、晉江市等有“高-高”聚集的空間特征，空間范圍上呈先縮小后擴張的趨勢；從時間上看，16個縣域森林病蟲害發(fā)病率有顯著下降趨勢，12縣域有顯著上升趨勢。采用STWR、GWR、OLS模型對比分析各變量對縣域森林病蟲害發(fā)病率的影響，發(fā)現(xiàn)OLS及GWR模型擬合效果一般，考慮時空異質(zhì)性的STWR模型擬合結(jié)果優(yōu)勢較明顯，且STWR的預(yù)測能力優(yōu)于GWR。STWR可以更精細(xì)地測度各變量與縣域森林病蟲害發(fā)病率之間的時空非平穩(wěn)關(guān)系，具體包括：

(1)影響因素從強到弱依次為夜間燈光均值、月均降水量、夏季均溫、冬季均溫。其中夜間燈光均值、夏季均溫、冬季均溫對森林病蟲害發(fā)生起正向促進(jìn)作用，月均降水量多為負(fù)向抑制作用。

(2)4個因素對森林病蟲害的影響機制具有時空分異性：2014-2018年夜間燈光均值的正向影響位于南平市和寧德市；2010年夏季均溫在新羅區(qū)、漳浦縣的負(fù)向影響最顯著；2014年冬季均溫負(fù)向影響顯著集中在寧德市蕉城區(qū)、福安市、霞浦縣；2018年月均降水量在龍巖市上杭縣、武平縣、連城縣的正向影響極為顯著。

(3)森林病蟲害高發(fā)生率的縣域，夜間燈光均值的負(fù)向影響最大，月均降水量的負(fù)向作用較為顯著，夏季均溫、冬季均溫的回歸系數(shù)在零點附近波動，分布較為均勻。

STWR模型基于新的時空核函數(shù)解釋并分析局部時空非平穩(wěn)性，描繪了各變量對縣域森林病蟲害發(fā)生率的時空變化過程，對確定森林病蟲害的主導(dǎo)因素、因時因地制定防治策略、提高防治水平有一定的參考價值。森林病蟲害發(fā)生是一個復(fù)雜的時空變化過程，受人類活動、氣候、森林發(fā)育情況等多種要素共同影響。然而，由于數(shù)據(jù)獲取的有限性和專業(yè)知識限制，本研究僅選取了夜間燈光均值、夏季均溫、冬季均溫、年均降水量4個指標(biāo)進(jìn)行建模分析，后續(xù)研究中將考慮林分、光照、森林類型等更多因素。此外，STWR模型在病蟲害的未來發(fā)病率預(yù)測能力上仍需更深入的研究。