運用地理加權泊松模型估測天然次生林枯損量分布1)

張凌宇 劉兆剛 董靈波

(東北林業(yè)大學，哈爾濱，150040)

天然次生林具有樹種繁多、分布范圍廣等特點，在我國東北、西南和東南地區(qū)均有分布，天然次生林不僅作為用材林以及林副產(chǎn)品的重要生產(chǎn)基地，在涵養(yǎng)水源以及維持生態(tài)平衡方面等都起到了重要的作用。我國次生林面積約占全國森林面積的46.2%，對天然次生林進行全面規(guī)劃、改善林分質量，使其盡快成為重要的森林資源。林分枯損影響著森林資源的動態(tài)變化，在森林演替過程中也起著至關重要的作用，因此，對天然次生林枯損株數(shù)分布的研究具有重要的理論意義。

在20世紀60年代，我國學者就已經(jīng)開始了對林分枯損的研究，李克志[1]利用帶嶺地區(qū)4種天然林的實測數(shù)據(jù)，對紅松等4種天然林枯損率進行了計算，對各齡組的枯損原因進行了分析；曾偉生[2]利用2006年西藏地區(qū)云杉林的實測固定樣地數(shù)據(jù)，在非線性混合估算模型的基礎上，構建了云杉林消耗率、枯損率、收獲率三者一致兼容的估算模型，這一方法為西藏自治區(qū)天然云杉林的枯損量預估提供了參考依據(jù)。隨著統(tǒng)計模型在林業(yè)上的應用，一些學者開始利用全局模型對林分枯損株數(shù)進行預測和分析，張雄清等[3]利用落葉松(Larixolgensis)林分連續(xù)觀測數(shù)據(jù)，分別利用泊松回歸模型、負二項模型、零膨脹模型擬合林分枯損株數(shù)。由于在林業(yè)調(diào)查中所收集的數(shù)據(jù)往往處于不同的地理位置，無論是研究枯損概率的logistic模型，還是研究枯損株數(shù)的泊松模型，實際上都沒有考慮到空間異質性，因此，空間異質性在有關林業(yè)的各項研究中普遍存在，也往往是研究者最容易忽略的問題。

20世紀90年代，英國學者提出了地理加權回歸模型(GWR)，解決空間非平穩(wěn)性以及空間異質性問題[4-5]。在國內(nèi)林業(yè)領域，有關GWR模型的研究一般集中在森林碳儲量以及生物量預估等方面[6-8]，將地理加權回歸模型應用到林分枯損株數(shù)空間分布中的研究還未見報道。地理加權泊松模型(GWPR)是以地理加權回歸模型(GWR)為基礎擴展建立的局域泊松模型，可以用來構建以林木枯損株數(shù)為因變量的局域模型。本研究采用帽兒山實驗林場2期固定樣地數(shù)據(jù)，分別建立了模擬天然次生林枯損株數(shù)的全局泊松模型和局域模型，將為大范圍內(nèi)天然次生林枯損株數(shù)預測以及解決計數(shù)型數(shù)據(jù)的空間非平穩(wěn)性問題提供理論基礎。

1 研究區(qū)概況

帽兒山實驗林場位于黑龍江省尚志市西部，地理坐標127°30′～127°34′E、45°20′～45°25′N，總面積26 496 hm2，共劃分10個森林經(jīng)營施業(yè)區(qū)、151個林班。帽兒山實驗林場的植被類型屬于長白山植物區(qū)系，是地帶性頂級植被闊葉紅松林被人為干擾破壞后形成的較典型天然次生林。主要林分類型有楊樺林、珍貴硬闊混交林、色木(Acermono)、柞樹(Quercusmongolica)為主的硬雜木林,其它主要樹種還有水曲柳(Fraxinusmandshurica)、胡桃楸(Junglusmandshurica)、黃菠蘿(Phellodendronamurense)、楓樺(Betualcostata)、白樺(Betualplatyphlla)、山楊(Populusdavidiana)等[9]。

2 研究方法

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

本研究數(shù)據(jù)來源于帽兒山實驗林場森林資源二類調(diào)查的2期固定樣地數(shù)據(jù)(2004年和2016年)，在所有的固定樣地中剔除宜林沼澤、非林地、農(nóng)田等樣地，同時檢查數(shù)據(jù)中的異常值，對于異常值予以剔除，最終選擇101塊天然次生林樣地作為研究數(shù)據(jù)，樣地面積均為0.06 hm2，數(shù)據(jù)涵蓋研究區(qū)內(nèi)的所有10個施業(yè)區(qū)。

2.2 模型變量的選擇

采用混合逐步選擇法對模型變量進行選擇，混合逐步選擇法是在常規(guī)模型的基礎上，根據(jù)一定的顯著性標準，逐步將所選變量剔除或添加進模型的方法[10]，本研究將顯著性標準設為α=0.05，最終得到海拔、坡度、林分平均胸徑、林分密度和蓄積量共5個影響次生林枯損株數(shù)分布的變量。以往國內(nèi)外研究結果表明，林分枯損株數(shù)主要與林分密度、林分平均胸徑以及單位面積蓄積量有關[11-12]，同時，林分單位面積枯損株數(shù)與土壤因子、立地條件和林分密度有著密切的關系[13]，而坡度以及海拔作為重要的地形因子影響著土壤水分以及養(yǎng)分的含量[14]，間接影響著天然次生林枯損株數(shù)分布。在所有的模型中對獨立變量均采取了標準化處理，模型各變量的統(tǒng)計量見表1。

表1 模型各變量的基本統(tǒng)計量

2.3 全局泊松模型

全局泊松模型(GP)作為分析計數(shù)型數(shù)據(jù)的一種方法，來模擬天然次生林枯損株數(shù)的分布情況，其概率密度函數(shù)如下：

(1)

式中：P(Y=y)表示在一段時間內(nèi)時間Y發(fā)生次數(shù)的概率；λ為隨機變量Y的期望和方差，且期望和方差相等，即E(Y)=λ，Var(Y)=λ；當y=0時，表示天然次生林枯損株數(shù)為0的概率；當y=1時，表示天然次生林枯損株數(shù)為1的概率；當y=k時，即P(y=k)=e-λλk/k!，表示枯損株數(shù)為k的概率。全局泊松回歸模型的形式如下：

log(E(Y))=logμ=β0+β1x1+β2x2+β3x3+β4x4+

β5x5。

(2)

式中：log(.)為一個鏈接函數(shù)；μ表示枯損株數(shù)；β1～β5為模型的回歸系數(shù)；x1～x5分別為海拔、坡度、平均胸徑、林分密度以及單位面積蓄積量。

2.4 地理加權泊松模型

地理加權泊松模型屬于全局泊松回歸模型在局域形式上的一種表達，在計算過程中可以將每塊樣地的地理坐標都納入到模型中，具體形式如下：

log(e(Y))=logμ=β0(ui,vi)+β1(ui,vi)x1+β2(ui,vi)x2+β3(ui,vi)x3+β4(ui,vi)x4+

β5(ui,vi)x5。

(3)

式中：log(.)為一個鏈接函數(shù)；β1(ui,vi)～β5(ui,vi)為樣地點(ui,vi)上的回歸系數(shù)；x1～x5分別為海拔、坡度、平均胸徑、林分密度以及單位面積蓄積量。

2.5 模型殘差的空間自相關性分析

真實值與預測值之間的差為模型殘差，通過引入莫蘭指數(shù)(MoranI)這一概念來檢驗模型殘差中的空間相關性，近年來，莫蘭指數(shù)被普遍運用到在林業(yè)研究中[6,15]。全局莫蘭指數(shù)計算公式如下：

(4)

(5)

在研究中，首先設定一個零假設，即當莫蘭指數(shù)的值為零時，假設模型殘差服從隨機分布，對于本研究而言，即天然次生林枯損株數(shù)的空間隨機分布。通過計算莫蘭指數(shù)，給出一個可以判斷出是否拒絕零假設的閾值，稱之為Z值，如果Z值在-1.96～1.96的區(qū)間內(nèi)，則P>0.05(α=0.05)，表示不能拒絕模型殘差服從隨機分布這一零假設，認為模型殘差表現(xiàn)出的模式在很大程度上是隨機空間過程產(chǎn)生的結果；反之，如果不在這一區(qū)間內(nèi)，則認為模型殘差的空間模式不是隨機過程產(chǎn)生的，在這種情況下可以拒絕零假設，所表現(xiàn)出來的空間分布模式為具有統(tǒng)計性顯著的聚類或分散模式，這種方法可以判斷n個樣點之間的觀測值是否存在空間自相關性[6]。

2.6 模型評價

本研究使用均方誤差(MSE)以及赤池信息準則(AIC)值來比較2種模型的擬合效果，當MSE值越小時，說明模型的精度越高，模型的擬合效果越好。具體公式如下：

(6)

赤池信息準則是在熵的概念基礎上建立的一種衡量統(tǒng)計模型擬合效果優(yōu)良的指標[16]，在模型比較時，AIC值越小的模型，說明該模型擬合效果越好。具體公式如下：

AIC=2p-2ln(L)。

(7)

式中：L為最大似然函數(shù)值；p為變量個數(shù)。

本文采用SAS9.4通過混合逐步選擇法對影響天然次生林枯損株數(shù)的變量進行篩選，并建立全局泊松模型；采用GWR4.0建立局域GWPR模型，采用ROOKCASE軟件來計算全局莫蘭指數(shù)以及局域莫蘭指數(shù)，利用反距離權重法(IDW)對局域模型各變量對應參數(shù)的空間分布圖、局域莫蘭指數(shù)的空間分布圖、局域Z值的空間分布圖以及全局和局域模型枯損株數(shù)的空間分布圖分別進行繪制。

3 結果與分析

3.1 模型擬合

3.1.1 全局泊松模型和局域模型的擬合結果

表2分別給出了全局泊松模型參數(shù)估計值、標準誤差、顯著性檢驗以及局域模型系數(shù)的擬合結果。對于全局泊松模型來說，所有的模型系數(shù)均通過了α=0.05水平下的顯著性檢驗，海拔、坡度、林分密度以及單位面積蓄積量的參數(shù)估計值均為正值，說明這4個變量與天然次生林枯損株數(shù)呈正相關關系，相反地，只有林分平均胸徑的參數(shù)估計值為負值，說明林分平均胸徑與天然次生林枯損株數(shù)呈負相關關系。在所有影響天然次生林枯損株數(shù)的因子中，林分因子對枯損株數(shù)的影響最大，其中林分平均胸徑對次生林枯損株數(shù)的影響最大，其次是林分密度，最后是單位面積蓄積量，海拔和坡度等地形因子對于次生林枯損株數(shù)的影響要小于林分因子。

表2 全局泊松模型和局域模型的參數(shù)估計值

全局泊松模型只有一組固定的模型參數(shù)，局域模型的參數(shù)則是連續(xù)變化的，表2給出了5組模型參數(shù)統(tǒng)計量，分別為最小值、Q1(四分位數(shù)下限值)、中值、Q3(四分位數(shù)上限值)以及最大值。

由表2可知，局域模型參數(shù)隨著地理位置的不同產(chǎn)生了一定的變化范圍，在全局泊松模型中β0=1.668 3，β0的參數(shù)范圍為-0.127 2～3.043 2，同理，在全局泊松模型中β1=0.452 7，β2=0.404 7，β3=-1.427 2，β4=1.128 3，β5=0.966 4，對應局域模型中β1的范圍為-1.435 9～1.013 4，β2的范圍為-0.238 9～2.044 1，β3的范圍為-3.362 9～2.004 8，β4的范圍為-1.575 7～5.192 5，β5的范圍為-2.967 6～2.829 3。

3.1.2 局域模型各參數(shù)的空間分布

從圖1A可以看出，在研究區(qū)的東南部地區(qū)以及東北部小范圍內(nèi)，海拔參數(shù)的估計值為負值，說明在這些區(qū)域內(nèi)，隨著海拔的逐漸增高，天然次生林枯損株數(shù)逐漸減少，而在研究區(qū)的中部地區(qū)以及西部地區(qū)，海拔參數(shù)的估計值為正值，說明在這些區(qū)域內(nèi)枯損株數(shù)隨著海拔高度的增加而增加。

圖1 GWPR模型各變量對應參數(shù)的空間分布

圖1B中，在研究區(qū)的北部以及南部地區(qū)很小范圍內(nèi)的坡度參數(shù)值為負值，所占區(qū)域的比例很小，在大部分范圍內(nèi)，隨著坡度的逐漸增加，林分枯損株數(shù)逐漸增加。

圖1C中，西部地區(qū)及東南部分地區(qū)林分平均胸徑的參數(shù)估計值為正，說明在這兩個區(qū)域內(nèi)，林分平均胸徑與枯損株數(shù)呈正相關關系，在其他區(qū)域內(nèi)林分平均胸徑與枯損株數(shù)呈負相關關系，從總體上看，在整個研究區(qū)范圍內(nèi)，自西南向東北方向林分平均胸徑參數(shù)估計值的變化范圍從正到負。

在圖1D中，林分密度的參數(shù)估計值除了在東北部地區(qū)極小范圍內(nèi)為負值外，在其他大部分地區(qū)林分密度的參數(shù)估計值均為正值，說明在這些區(qū)域內(nèi)，隨著林分密度的增加，天然次生林枯損株數(shù)也隨之增加。

在圖1E中，林分單位面積蓄積量的參數(shù)估計值自東南至西北方向均為負值，參數(shù)估計值為負值的區(qū)域約占整個研究區(qū)域面積的50%，說明在此研究區(qū)域內(nèi)，隨著蓄積量的增加，林分枯損株數(shù)也隨之增加，在其他區(qū)域內(nèi)林分單位面積蓄積量與枯損株數(shù)呈正相關關系?？偟膩碚f，在局域模型下各參數(shù)分級的細化程度較為清晰，同時形成了較好范圍內(nèi)的局域化空間分布效果。

3.2 模型擬合優(yōu)度比較

根據(jù)兩種模型的均方誤差和赤池信息準則值的計算結果可知，全局泊松模型赤池信息準則值為702.14，局域泊松模型的赤池信息準則值為192.28，說明局域模型在擬合數(shù)據(jù)的優(yōu)良性上產(chǎn)生了較好的效果；同時局域模型的均方誤差為14.64，全局泊松模型的均方誤差值為37.28，說明局域模型的擬合效果要好于全局模型。

3.3 空間自相關性分析

3.3.1 模型殘差的全局空間自相關性

圖2顯示了處于不同空間滯后距離(Lag)狀態(tài)下，全局泊松模型和局域模型殘差的莫蘭指數(shù)(MoranI)相關圖。從總體上看，在空間滯后距離大致相同的情況下，全局泊松模型的莫蘭指數(shù)值要明顯大于局域模型，表明在大部分情況下，全局泊松模型的全局空間相關性要明顯高于局域模型；隨著空間滯后距離的增加，全局泊松模型的莫蘭指數(shù)值逐漸減小，而局域模型的莫蘭指數(shù)值基本在0點上下浮動，起伏沒有全局泊松模型明顯，說明在降低模型殘差全局空間自相關性的效果上，局域模型要明顯好于全局泊松模型。此外，全局泊松模型(GP)在空間滯后距離為14 km時，產(chǎn)生最小的莫蘭指數(shù)值，局域模型(GWPR)在空間滯后距離為5 km時產(chǎn)生最小的莫蘭指數(shù)值。

圖2 全局泊松模型和局域模型殘差的莫蘭指數(shù)相關圖

3.3.2 模型殘差的局域空間自相關性

圖3給出了全局泊松模型和局域模型殘差的莫蘭指數(shù)值的空間分布圖。

由圖3A可知，在本研究區(qū)域的中部地區(qū)，全局泊松模型產(chǎn)生了較大范圍且取值為正值的局域莫蘭指數(shù)值，說明模型殘差在該區(qū)域內(nèi)傾向于出現(xiàn)相同觀測值的聚類模式；在圖3B中，局域模型則在相同區(qū)域內(nèi)產(chǎn)生了取值為負且范圍較小的局域莫蘭指數(shù)值，說明在局域模型下模型殘差出現(xiàn)了不同觀測值的聚類模式，不同觀測值的聚類模式可以很好的揭示空間異質性，產(chǎn)生了較為理想的效果。

圖3 不同模型殘差的局域莫蘭指數(shù)的空間分布

圖4給出了在α=0.05顯著性水平下，全局泊松模型和局域模型殘差局域莫蘭指數(shù)的Z值空間分布。GWPR模型(圖4B)的局域Z值在大部分情況下均在-1.96～1.96，在P>0.05的情況下，各樣地點間無顯著性差異；全局泊松模型(圖4A)的Z值呈現(xiàn)顯著性的區(qū)域，說明局域模型在降低局域空間自相關性的效果上要明顯好于全局泊松模型。

圖4 不同模型殘差的局域莫蘭指數(shù)的Z值空間分布

3.4 枯損株數(shù)的空間分布

由圖5可知，從總體上看，兩種模型在對天然次生林枯損株數(shù)空間分布的模擬上與樣地內(nèi)實際枯損株數(shù)的空間分布大致相同，即研究區(qū)的東北部地區(qū)林木枯損株數(shù)較多，在該地區(qū)內(nèi)兩種模型的林木枯損株數(shù)均在15株以上；在研究區(qū)的東南部地區(qū)，全局泊松模型的林木枯損株數(shù)在10～15株，而對應局域模型的枯損株數(shù)為8～10株；此外，兩種模型的擬合偏差存在差異，全局泊松模型對樣地內(nèi)枯損株數(shù)模擬的擬合偏差在6株左右，局域模型的擬合偏差在3～4株，要明顯低于全局泊松模型，同時在分級的細化程度上，局域模型也要好于全局泊松模型。

圖5 不同模型林木枯損量的空間分布

4 結論與討論

采用全局泊松模型和局域模型(GWPR)對帽兒山實驗林場101塊天然次生林樣地內(nèi)的林木枯損株數(shù)進行了模擬，并對2種模型的擬合效果以及空間自相關性進行了分析。對于全局泊松模型，林分平均胸徑與次生林枯損株數(shù)呈負相關關系，且影響最大，對于局域模型來說，隨著地理位置的不同，模型參數(shù)產(chǎn)生一定的變化范圍，各個參數(shù)對于枯損株數(shù)的影響大小也隨之變化；在模型擬合方面，GWPR模型擁有很高的擬合精度，AIC值明顯減小，得到了較好的模型參數(shù)的局域化空間分布效果；在解決空間自相關性上，GWPR模型殘差無論是在全局空間自相關性還是局域空間自相關性上都要小于全局模型，同時產(chǎn)生了不同觀測值少量聚類這一理想的空間分布模式；在對枯損株數(shù)空間分布的模擬上，局域模型的擬合偏差也要小于全局模型。

在林業(yè)數(shù)據(jù)研究中，一般地勢地形不同都會導致研究的數(shù)據(jù)之間存在空間差異，因此空間自相關性在有關林業(yè)的各項研究中普遍存在，同時，泊松模型屬于計數(shù)模型的范疇，在以往的研究中將GWR模型應用到林分枯損株數(shù)空間分布中的研究還未見報道。通過2個模型的空間自相關性分析發(fā)現(xiàn)，全局泊松模型的空間自相關性隨著空間滯后距離的增加逐漸下降，局域模型的莫蘭指數(shù)值在0點上下浮動，這是由于局域模型在每塊樣地內(nèi)都可以得到一組模型參數(shù)，每組模型參數(shù)都會形成一個GWPR模型，最后綜合成一個模型，局域模型在運行過程中充分考慮了地形因子變化對林木枯損株數(shù)的影響，從根本上解決了空間異質性的問題，而全局泊松模型的所有樣地只有一組模型參數(shù)，在可視化分布圖的平滑性以及精度上要小于局域模型(GWPR)。

有關研究結果表明，選取不同的空間尺度可能對空間自相關性的大小產(chǎn)生一定的影響[17-18]，因此，在以后的研究中，應適當結合研究區(qū)的面積，選取不同空間尺度下的局域模型，對天然次生林的枯損株數(shù)進行探索。此外，本文是以天然次生林為一個整體進行研究的，所選取的天然次生林樣地大多為軟闊林、硬闊混交林以及軟闊混交林。由于天然次生林樹種繁多，因此，不同的森林類型可能會導致林分枯損株數(shù)的不同，在今后的研究中，可以按照不同林型對天然次生林枯損株數(shù)的空間分布進行探討。

[1] 李克志.小興安嶺帶嶺林區(qū)天然林材積枯損率的初步研究[J].林業(yè)科學,1958(1):97-111.

[2] 曾偉生.西藏天然云杉林枯損率與采伐率模型系統(tǒng)研究[J].林業(yè)科學研究,2008,21(3):353-356.

[3] 張雄清,雷淵才,雷相東,等.基于計數(shù)模型方法的林分枯損研究[J].林業(yè)科學,2012，48(8):54-61.

[4] FOTHERINGHAM A S, CHARLTON M, BRUNSDON C. The geography of parameter space: an investigation of spatial non-stationarity[J]. Geographical Information Systems,1996,10(5):605-627.

[5] BRUNSDON C, FOTHERINGHAM A S, CHARLTON M E. Geographically Weighted Regression: A method for exploring spatial nonstationarity[J]. Geographical Analysis,1996,28(4):281-298.

[6] 劉暢.黑龍江省森林碳儲量空間分布研究[D].哈爾濱:東北林業(yè)大學,2014.

[7] 劉暢,李鳳日,甄貞.空間誤差模型在黑龍江省森林碳儲量空間分布的應用[J].應用生態(tài)學報,2014,25(10):2779-2786.

[8] 戚玉嬌.大興安嶺森林地上碳儲量遙感估算與分析[D].哈爾濱:東北林業(yè)大學,2014.

[9] 盧軍.帽兒山天然次生林樹冠結構和空間優(yōu)化經(jīng)營[D].哈爾濱:東北林業(yè)大學,2008.

[10] 王濟川.Logistic回歸模型：方法與應用[M].北京：高等教育出版社,2001.

[12] STEPHENSON N L, FRANKLIN J F. Causes and implications of the correlation between forest productivity and tree mortality rates[J]. Ecological Monographs,2011,81(4):527-555.

[13] 王岳,王海燕,楊小娟,等.基于土壤因子的近天然落葉松云冷杉林枯損模型研究[J].福建農(nóng)林大學學報(自然科學版),2015,44(4):379-383.

[14] 劉鑫,畢華興,李笑吟,等.晉西黃土區(qū)基于地形因子的土壤水分分異規(guī)律研究[J].土壤學報,2007,44(3):411-417.

[15] 劉暢,李鳳日,賈煒瑋,等.基于局域統(tǒng)計量的黑龍江省多尺度森林碳儲量空間分布變化[J].應用生態(tài)學報,2014,25(9):2493-2500.

[16] 宋喜芳,李建平,胡希遠.模型選擇信息量準則AIC及其在方差分析中的應用[J].西北農(nóng)林科技大學學報(自然科學版),2009，37(2):88-92.

[17] MA Z H, BENJAMI Z, WILLIAM F, et al. Use of localized descriptive statistics for exploring the spatial pattern changes of bird species richness at multiple scales[J]. Applied Geography,2012,32(2):185-194.

[18] GUO L G, MA Z M, ZHANG L Z. Comparison of bandwidth selection in application of geographically weighted regression: A case study[J]. Canadian Journal of Forest Research,2008,38(9):2526-2534.