何 瀟 雷相東

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

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

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

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

1 研究區(qū)概況

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

2 數(shù)據(jù)與方法

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

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

2.2 研究方法

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

(1)

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

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

2.2.2 普通回歸模型 由于當(dāng)前開發(fā)的同步自回歸模型(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的空間自相關(guān)性(Luetal.， 2011； Kisslingetal.， 2008)。SEM考慮自變量的空間自相關(guān)性，即某一空間的因變量與相同空間內(nèi)的自變量、相鄰區(qū)域的自變量有關(guān)； SLM考慮因變量的空間自相關(guān)性，即某一空間的因變量與相同空間內(nèi)的自變量、相鄰區(qū)域的因變量有關(guān)。SAR 模型表達(dá)式(Anselin，1988)如下:

SEMY=Xβ+λWu+ε；

(2)

SLMY=ρWY+Xβ+ ε。

(3)

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

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

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

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

(4)

(5)

(6)

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

(7)

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

3 結(jié)果與分析

3.1 BCEF以及各林分因子的空間自相關(guān)性

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

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

3.2 備選模型優(yōu)選

由表3可知，不同林分變量和模型對BCEF的擬合度和解釋能力不同，不同林分變量對應(yīng)的最優(yōu)模型也不一樣。從模型形式上看，指數(shù)模型的參數(shù)估計值均有異常，說明該模型不適合描述BCEF； S形曲線模型的評價指標(biāo)均不理想，R2最大不超過0.50，RMSE均超過0.01 mg·m-3，rRMSE超過1.30%； 異速生長模型、對數(shù)模型和雙曲線模型擬合結(jié)果較好，模型的R2大多在0.397～0.958之間。從自變量上看，Dg是擬合BCEF的最優(yōu)變量，模型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模型，其擬合效果較好，除指數(shù)模型、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中可看出普通回歸模型的殘差存在明顯的空間自相關(guān)性，其中模型14和20的MI差別不大，模型2的MI最大。可見，需要考慮空間效應(yīng)對BCEF的影響。

表3 基于最小二乘法的模型參數(shù)估計值與模型評價指標(biāo)①

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

3.3 SAR模型擬合結(jié)果與評價

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

表4 以Dg為自變量的SAR模型參數(shù)估計值與模型評價指標(biāo)

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

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

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

4 討論

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

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

5 結(jié)論

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