東北地區落葉松人工林生物量轉換與擴展因子空間自回歸模型*

何 瀟 雷相東

(中國林業科學研究院資源信息研究所 北京 100091)

森林是陸地生態系統中最大的碳匯和碳源，占整個陸地碳庫的56%，是全球碳循環的重要組成部分(董文福等， 2002)。在全球碳循環過程和森林生態系統固碳服務功能量化中，森林生物量特征(如數量、分布和動態等)被認為是森林生長和收獲預測的重要內容和評價指標之一(Canadelletal.， 2008； Panetal.， 2011)。IPCC(2003； 2006)指出，森林生物量計量可選用生物量異速生長方程法、生物量轉換與擴展因子(biomass conversion and expansion factors，BCEF)法、模型模擬法等，其中BCEF法(林分喬木層生物量與林分蓄積量之比)在區域森林生物量計量中得到廣泛應用，實現了生物量估算從樣地到區域的尺度轉換(Fangetal.， 2001； 李海奎等， 2012； Levyetal.， 2004； Liskietal.， 2006)，是IPCC(2006)重點推薦的大尺度森林生物量估算方法。一般來說，大尺度森林生物量估算是基于森林清查數據進行的(Shvidenkoetal.， 2002； 張茂震等， 2009)，但使用的BCEF大都基于少數代表性有限的區域研究(Joostenetal.， 2004)或IPCC給出的固定值，BCEF隨林分發育階段不同會有所變化(Lehtonenetal.， 2007； 羅云建等， 2009)，這將導致森林生物量估算的不確定性增加。Brown(2002)和Jalkanen等(2005)詳細描述了BCEF隨林分蓄積量的變化，蓄積量較小時具有較大的BCEF，而蓄積量較大時BCEF呈指數下降并趨于穩定。因此，為了減少使用BCEF進行森林生物量估算的不確定性，應針對特定研究區和樹種建立BCEF與森林調查因子之間的關系。

目前，很多研究探討了BCEF與蓄積量或其他林分變量之間的關系，如Fang等(2001)利用雙曲線模型表示BCEF與林分蓄積量之間的關系； Soares等(2012)采用分段雙曲線模型，并選擇優勢高作為BCEF預測變量； Peichl等(2007)和Lehtonen等(2004)建立了對年齡敏感的BCEF模型； Jagodziński等(2019)發現基于林分年齡和蓄積量的雙曲線BCEF模型精度較高； 方精云等(1996)和羅云建等(2007； 2010)計算了落葉松(Larix. spp)林的BCEF值。傳統研究方法假定各樣地間的BCEF是相互獨立的(Mateuetal.， 1998)，采用最小二乘法估計模型參數。然而，從林學的角度看，各樣地間的林分因子數據并非獨立，尤其在大區域尺度上具有明顯的空間差異，因此使用普通模型建立的森林生長收獲和預估模型在預測時將導致有偏估計(Mengetal.， 2009)； 從統計學的角度看，忽略空間自相關會造成模型殘差不符合獨立性假設，模型參數有效性和模型估計效果降低，從而增加生物量計量結果的不確定性(Westetal.， 1984； Gregorie， 1987)。

綜上可知，關于BCEF的研究存在以下問題： 1) BCEF的最優模型并不一致； 2) BCEF與林分平方平均直徑、斷面積、密度、蓄積量和平均年齡等林分因子的關系存在不確定性； 3) 大部分研究只是利用簡單相關性分析和普通回歸方法分析BCEF與林分變量之間的關系，并沒有考慮BCEF數據的空間自相關性和利用空間信息。鑒于此，本研究以東北地區3省(黑龍江、吉林和遼寧)的落葉松人工林為對象，基于最小二乘普通回歸方法和空間自回歸方法分別建立BCEF模型，并比較2種方法的模型參數、模型殘差在空間上的差異性以及模型擬合效果，以期為落葉松人工林生物量精準估算提供模型支撐依據。

1 研究區概況

以我國東北黑龍江、吉林和遼寧3省為研究區，屬溫帶季風性氣候，冬季漫長寒冷，夏季短促溫暖，土壤以黑土為主，土地肥沃。研究區主要針葉樹種有落葉松、云杉(Piceaasperata)、冷杉(Abiesfabri)、紅松(Pinuskoraiensis)等，主要闊葉樹種有白樺(Betulaplatyphylla)、水曲柳(Fraxinusmandshurica)、黃檗(Phellodendronamurense)等。對于落葉松純林，主要建群種有興安落葉松(L.gmelinii)、長白落葉松(L.olgensis)、日本落葉松(L.kaempferi)和少量華北落葉松(L.principis-rupprechtii)等。

2 數據與方法

2.1 數據來源

數據來源于黑龍江、吉林和遼寧3省第6～8次全國森林資源連續清查中落葉松人工純林固定樣地。根據一類調查規程，將落葉松蓄積量占比65%以上的林分稱為落葉松純林。總計1 005塊樣地，刪除株數小于20株的樣地，以BCEF為目標變量，采用箱線圖法剔除異常數據，共得到788塊有效樣地數據，緯度在39.68°—53.19°N之間，經度在122.30°—133.02°E之間。基于每木檢尺數據得到每塊樣地的林分斷面積(basal area,BA)、平方平均直徑(quadratic mean diameter,Dg)和密度(density,N)等基本因子,林分平均高(mean height,H)為3～5株平均木的樹高平均值。樣地基本因子統計量見表1。觀測數據林分平均年齡(mean age,Age)在8～52年之間，數據覆蓋廣，有代表性。

2.2 研究方法

2.2.1 樣地BCEF 計算 BCEF(mg·m-3)是林分生物量(mg·hm-2)與林分蓄積量(m3·hm-2)之間的轉換系數，計算公式如下：

(1)

式中： BCEFi表示第i塊樣地的生物量轉換與擴展因子；bij表示第i塊樣地第j株樹的總生物量(干、枝、葉、根生物量之和)；vij表示第i塊樣地第j株樹的材積；ni表示第i塊樣地內的林木株數。

樣地BCEF計算需將樣地內每株林木的生物量和材積累加，得到樣地尺度的生物量和蓄積量。對于落葉松，本研究采用LY/T 2654—2016《立木生物量模型及碳計量參數——落葉松》(國家林業局， 2016)中的一元生物量模型和一元材積模型，基于單木胸徑計算相應的生物量和材積； 由于標準中生物量和材積模型是來自同一套建模樣本所建立的兼容模型，因此樣地BCEF計算能夠保持穩定。對于其他伴生樹種，單木生物量采用文獻(陳傳國等， 1989； Wang， 2006)中的模型計算，并利用國家森林資源連續清查樣木表中提供的材積結果(其根據各省頒布的材積表估算，鑒于伴生樹種數量較少，能夠保持BCEF的穩定性)。

2.2.2 普通回歸模型 由于當前開發的同步自回歸模型(SAR)為線性模型(Mengetal.， 2009)，因此本研究首先將5種常見BCEF非線性模型通過變量變換進行線性化(模型具體形式見表2)，作為備選模型，然后建立BCEF普通回歸模型，自變量包括V、Dg、BA、H、N和Age等。

表2 線性化的備選模型①

2.2.3 空間自回歸(SAR)模型 運用空間誤差模型(spatial error model,SEM)和空間滯后模型(spatial lag model,SLM)2種空間自回歸方法處理不同采樣點間BCEF的空間自相關性(Luetal.， 2011； Kisslingetal.， 2008)。SEM考慮自變量的空間自相關性，即某一空間的因變量與相同空間內的自變量、相鄰區域的自變量有關； SLM考慮因變量的空間自相關性，即某一空間的因變量與相同空間內的自變量、相鄰區域的因變量有關。SAR 模型表達式(Anselin，1988)如下:

SEMY=Xβ+λWu+ε；

(2)

SLMY=ρWY+Xβ+ ε。

(3)

式中：Y表示因變量組成的n×1維向量；X表示自變量矩陣，本研究中是由林分因子變量或其變量變換后組成的n×(p+1)維矩陣，p為自變量個數；β表示模型待估參數(含截距)；ε表示模型殘差；W表示空間加權矩陣，本研究對不同空間加權矩陣進行測試，只有在相接鄰近矩陣條件下模型收斂，最終采用的空間加權矩陣為相接鄰近矩陣；u是誤差項，Wu為空間誤差效應，用來表示某一樣地的BCEF不僅與同一樣地的自變量有關，還與相接鄰樣地的自變量有關，即空間誤差效應；WY為空間滯后效應，用來表示某一樣地的BCEF不僅與同一樣地的自變量有關，還與相接鄰樣地的BCEF有關，即空間滯后效應；WY和Wu均遵循空間自回歸過程；λ和ρ分別表示空間自回歸模型中的空間自相關誤差項的待估參數和因變量空間滯后項的待估參數。

需要指出，當λ=0時，式(2)轉化為普通線性回歸模型； 當ρ=0時，式(3)轉化為普通線性回歸模型。

本研究采用R語言spdep包的lagsarlm函數和errorsarlm函數估計SAR模型參數。

2.2.4 模型評價指標 采用決定系數(determination coefficient,R2)、均方根誤差(root mean square error, RMSE)和相對均方根誤差(relative RMSE, rRMSE)對模型進行評價，具體計算公式如下：

(4)

(5)

(6)

2.2.5 空間自相關檢驗指標 莫蘭指數(Moran’s I, MI)是衡量某一隨機變量空間自相關程度的指標，主要根據變量位置和變量值計算(Paradis， 2009)。其值分布在[-1，1]之間，若MI>0，表示空間正相關性，即隨著空間分布位置(距離)的聚集，相關性越顯著，其值越大，正相關性越明顯，呈現出聚類模式； 若MI<0，表示空間負相關性，即隨著空間分布位置(距離)的離散，相關性變得顯著，其值越小，空間差異越大，呈現出離散模式； 若MI=0，表示在空間呈隨機模式。計算公式如下：

(7)

若BCEF的MI接近0，表明BCEF不存在空間自相關性，若BCEF的MI大于或小于0，表明BCEF存在明顯的空間自相關性，此時須考慮空間效應對BCEF的影響。同理，若BCEF模型殘差的MI接近0，表明BCEF模型殘差不存在空間自相關性，使用普通線性回歸模型即可； 若BCEF模型殘差的MI大于或小于0，表明BCEF模型殘差存在明顯的空間自相關性，此時須建立BCEF的SAR模型。

3 結果與分析

3.1 BCEF以及各林分因子的空間自相關性

莫蘭指數(MI)可在不同距離水平上計算，從圖1可以看出，BCEF存在明顯的空間自相關性，當空間距離較小(<350 km)時，呈現出聚集模式，即BCEF相似的樣地聚集在一起，BCEF高值之間或較小值之間相互連接； 處于中等距離(350～700 km)時，呈現出分離模式，即BCEF相異的樣地聚集在一起，BCEF高值與較小值之間相互連接，隨著距離增加，空間自相關性逐漸減弱。這說明當空間距離較小時，同一省內的落葉松林BCEF屬性相似，隨著空間距離增加，各省之間的BCEF屬性差異逐漸體現出來，最終趨向隨機分布。V和Dg等也存在空間自相關性，MI在±0.15之間，Dg的空間分布模式與BCEF非常相似，BCEF的空間自相關性可能是由Dg等林分特征因子驅動的； 而N則不存在空間自相關性，MI僅在±0.06之間。整理建模數據發現，黑龍江省中北部多為興安落葉松，東部多為長白落葉松，吉林省多為長白落葉松，遼寧省多為興安落葉松、長白落葉松、日本落葉松和少量華北落葉松，不同省落葉松樹種分布有所不同，可見，BCEF以及各林分因子的空間自相關性與樹種有很大關系，建立BCEF模型時應該考慮空間效應對其產生的影響。

圖1 BCEF以及各林分因子的莫蘭指數隨距離的變化趨勢

3.2 備選模型優選

由表3可知，不同林分變量和模型對BCEF的擬合度和解釋能力不同，不同林分變量對應的最優模型也不一樣。從模型形式上看，指數模型的參數估計值均有異常，說明該模型不適合描述BCEF； S形曲線模型的評價指標均不理想，R2最大不超過0.50，RMSE均超過0.01 mg·m-3，rRMSE超過1.30%； 異速生長模型、對數模型和雙曲線模型擬合結果較好，模型的R2大多在0.397～0.958之間。從自變量上看，Dg是擬合BCEF的最優變量，模型2、14和20均有很好擬合能力，R2超過0.94，RMSE均不超過0.003 mg·m-3，rRMSE不超過0.35%。H是第2個解釋BCEF較多的變量，模型4和16的R2均在0.70以上，RMSE均不超過0.007 mg·m-3，rRMSE不超過0.91%。V是第3個解釋BCEF較多的變量，模型1和13的R2均在0.60以上，RMSE不超過0.01 mg·m-3，rRMSE也不大。最后是以Age為自變量的BCEF模型，其擬合效果較好，除指數模型、S形曲線模型和雙曲線模型外，其余模型的R2均在0.51～0.55之間，RMSE均不超過0.01 mg·m-3。BA和N對BCEF的解釋能力較差，R2均不超過0.50，RMSE在0.01 mg·m-3以上。綜合分析，各自變量對BCEF的擬合效果排序為Dg>H>V>Age。雖然模型2、14和20均有很好擬合能力，但從圖2中可看出普通回歸模型的殘差存在明顯的空間自相關性，其中模型14和20的MI差別不大，模型2的MI最大。可見，需要考慮空間效應對BCEF的影響。

表3 基于最小二乘法的模型參數估計值與模型評價指標①

圖2 基于普通回歸BCEF模型殘差的莫蘭指數隨距離的變化趨勢

3.3 SAR模型擬合結果與評價

優選的普通回歸模型是以Dg為自變量的異速生長模型、對數模型和雙曲線模型(模型2、14和20)，本研究基于上述模型建立對應的SAR模型。由表4可知，與普通回歸模型相比，SAR模型的評價指標均有所提升，但模型32和34的參數估計值不顯著，只有雙曲線模型在2種空間自回歸方法下參數估計值均顯著。與模型20相比，基于SEM的模型35的R2提高3%，RMSE和rRMSE分別降低33%和35%，基于SLM的模型36的R2僅提高0.08%，RMSE和rRMSE基本沒有變化，圖3顯示出SEM比SLM的擬合效果更好。SEM中的λ和SLM中的ρ均顯著，說明BCEF具有空間自相關性。結合模型35和36殘差的MI變化說明SEM可更好消除空間自相關性，模型殘差的MI均不超過0.02，而SLM和普通最小二乘回歸模型殘差的MI趨勢基本一致。因此，基于SEM的模型30更適合描述BCEF變化。從參數估計值來看，模型35和20的參數相差不是很大，a相差0.002，b相差0.018，說明模型形式和變量選擇是合適的。綜上，雙曲線模型是最適合的模型形式，基于SEM建立的以Dg為唯一自變量的BCEF模型具有最好穩定性。

表4 以Dg為自變量的SAR模型參數估計值與模型評價指標

圖3 BCEF空間自回歸模型殘差的莫蘭指數隨距離的變化趨勢

由于Dg和H是對BCEF解釋最多的自變量，實際工作中也容易調查，而V則需要利用材積表推算，工作量很大，另外人工林的Age是非常容易獲取的指標，因此本研究嘗試以Dg、H和Age為自變量建立多元SAR模型描述BCEF變化，且選擇最穩定的雙曲線模型形式。由表5可知，模型39、40、42、43和44均存在參數不顯著的現象，在多元模型中，與Age相關的參數均不顯著，可見Age僅適合在單方程中描述BCEF變化。從模型評價指標來看，SEM均優于SLM。模型37的評價指標最好，但是與H相關的參數符號發生了變化： 模型22中，BCEF隨H增加而降低，但是模型37中BCEF隨H增加而增加。檢查建模數據的散點圖后發現，BCEF隨H增加呈明顯降低趨勢，模型37雖然在統計上最優，但是并不符合實際規律，因此也不能被采用。綜上，多元模型不適合用來描述BCEF變化。

表5 以Dg、H和Age為自變量的多元SAR模型參數估計值與模型評價指標

4 討論

以線性化BCEF模型為基礎模型，以不同林分變量擬合普通模型，選擇出適合描述BCEF的模型形式和自變量，結果發現，異速生長模型、對數模型和雙曲線模型均可以用來描述BCEF，其中雙曲線模型效果最好，而指數模型和S形曲線模型對BCEF擬合效果不佳。不同模型的特征不同(如拐點、漸近線等)，目前用來描述BCEF變化的最優模型研究結果并不一致，如羅云建等(2010)和左舒翟等(2014)與本研究思路一樣，均選擇多種模型形式分別擬合，從中挑選出最優模型形式和解釋能力最高的變量，其結果也不盡相同。在眾多林分變量中，本研究發現BCEF與Dg、H和V的關系十分密切。但是對BCEF與Dg之間關系的研究結論并不一致，如左舒翟等(2014)發現Dg對BCEF有顯著影響，BCEF隨Dg增加而增加，后趨于穩定； 而李海奎等(2017)發現Dg對碳計量參數無顯著影響。本研究發現H對BCEF有顯著影響，BCEF隨H增加而降低，與左舒翟等(2014)、李海奎等(2017)和竹萬寬等(2020)的研究結果一致。BCEF與V之間的良好關系是實現生物量由林分尺度向大區域尺度擴展的基礎(Fangetal.， 2001)，Jagodziński等(2019)建立BCEF模型時也發現最佳預測指標是V，與本研究結果相似。雖然Age對BCEF的擬合效果也不錯，但在一元模型情況下，由于年齡本身與立地無關，因此該類模型具有較大不確定性，從模型參數值上發現，BCEF與Age呈J形關系，即BCEF隨林分平均年齡增加而減小，最后趨于穩定，與Jagodziński等(2019)的研究結果類似。另外，本研究發現BCEF與BA和N的關系并不密切，與Soares等(2012)的結果一致； 但竹萬寬等(2020)、李海奎等(2017)發現N對碳計量參數有顯著影響。產生上述差異的原因主要有2點： 1) 本研究和左舒翟等(2014)只考慮單一變量對BCEF的影響，李海奎等(2017)是在綜合多個變量作用的基礎上進行篩選，某些林分變量之間存在嚴格的數學函數關系，因此選擇其中一個變量時會導致另一個變量失去作用； 2) 研究尺度、樹種、林木起源以及區域不同。

從BCEF本身的空間自相關性以及普通回歸模型的殘差存在空間自相關性來看，由于普通回歸模型并沒有考慮BCEF本身以及林分變量之間空間自相關性的影響，因此有必要通過空間自回歸方法降低空間自相關。本研究從普通回歸模型中挑選出合適的模型形式和解釋變量，利用2個SAR模型構建基于空間自相關的東北地區落葉松人工林BCEF模型。在2個SAR模型中，SEM優于SLM，但二者模型評價指標均優于普通回歸模型，與Kissling 等(2008)、Meng 等(2009)和婁明華等(2016)的研究結論類似，不過SLM對模型殘差空間自相關性的消除不明顯。模型30殘差的MI變得很小，說明殘差在空間上不再聚集，可提高模型擬合效果，λ小于0說明每塊樣地的Dg均受相鄰樣地較恒定的顯著負向影響，進而影響BCEF。使用SEM時，需要提供樣地坐標，利用空間位置重新確定鄰接對象，求得空間權重，再利用模型參數對BCEF進行預測，這就要求在進行外業調查工作時要保證樣地的地理坐標測量嚴格按照相關規范執行。雖然SEM不能完全消除模型殘差的空間自相關性，但可顯著降低空間自相關性，對BCEF的估計更為準確，從而能在大區域尺度上得到更加可靠的生物量估計值。需要說明的是，本研究結果只適用于落葉松人工純林，對于其他森林類型能否得到相同結論還不能確定。此外，還有2個問題值得關注： 1) 人工林和純林的結構相對單一，在混交林中如何應用空間自回歸模型是一個需要深入研究的問題，即如何將不同樹種分開對待而又不改變其空間作用； 2) 當前的空間自回歸模型是基于線性形式開發的，要得到BCEF和林分變量之間的非線性關系，只能通過變量變換，這里涉及到轉換為原尺度上時是否需要進行矯正，更有效的手段應該是開發基于非線性形式的空間自回歸模型。

5 結論

異速生長模型、對數模型和雙曲線模型均可以用來描述我國東北地區落葉松人工林生物量轉換與擴展因子，其中雙曲線模型最穩定。林分平方平均直徑是解釋轉換與擴展BCEF最多的變量，其次是林分平均高，然后是蓄積量和年齡，而林分斷面積和密度的擬合效果均不佳。 同時多元模型容易發生參數估計值不合理或者參數不顯著等問題，因此應選擇林分平方平均直徑來模擬生物量轉換與擴展因子的變化。由于生物量轉換與擴展因子存在空間自相關性，空間自回歸模型評價指標優于普通回歸模型，因此應選擇單變量雙曲函數空間誤差模型降低空間自相關性。