潘 萍,孫玉軍,歐陽勛志,寧金魁,馮瑞琦,汪 清

1 北京林業(yè)大學(xué)林學(xué)院,北京 100083 2 江西農(nóng)業(yè)大學(xué)林學(xué)院,南昌 330045 3 鄱陽湖流域森林生態(tài)系統(tǒng)保護(hù)與修復(fù)國家林業(yè)和草原局重點(diǎn)實(shí)驗(yàn)室,南昌 330045

森林生態(tài)系統(tǒng)作為最大的陸地碳庫,其植被和土壤分別約占全球植被碳庫的86%以上和土壤碳庫的73%以上,在減緩氣候變化和維護(hù)全球碳循環(huán)等方面發(fā)揮著重要的作用[1-2]。因此,森林生態(tài)系統(tǒng)的固碳能力受到廣泛關(guān)注。碳密度是衡量森林生態(tài)系統(tǒng)固碳能力的重要指標(biāo)之一[3-4],也是準(zhǔn)確估算森林碳儲(chǔ)量的基礎(chǔ)。由于森林生態(tài)系統(tǒng)是一個(gè)具有較強(qiáng)的空間異質(zhì)性的復(fù)雜系統(tǒng),其碳密度受氣候、地形地貌及周圍其他環(huán)境等因素的影響[5],而這些因素在空間上也往往存在相關(guān)性,在空間上相互影響、相互作用。因此,從區(qū)域尺度上探明影響森林碳密度的主要因素及構(gòu)建估測模型,可為森林碳儲(chǔ)量估算及碳匯林業(yè)的經(jīng)營管理等提供科學(xué)依據(jù)。

在區(qū)域尺度上,森林碳密度/碳儲(chǔ)量的估測主要是通過構(gòu)建碳密度/碳儲(chǔ)量與其影響因子之間的關(guān)系模型。由于傳統(tǒng)的回歸模型在構(gòu)建中忽略了各變量之間的空間異質(zhì)性,因此,空間誤差模型[6]、空間滯后模型[7]、空間濾波模型[8]、空間Durbin模型[9]、地理加權(quán)回歸模型[10]等空間模型受到許多學(xué)者的關(guān)注。如李月等[11]、劉暢等[12]采用空間誤差模型構(gòu)建碳儲(chǔ)量與影響因素關(guān)系模型時(shí)均得出其擬合精度明顯優(yōu)于最小二乘模型；Lin等[13]、歐光龍等[14]、郭含茹等[15]基于地理加權(quán)回歸模型進(jìn)行研究,均表明地理加權(quán)回歸模型比最小二乘模型在碳密度/碳儲(chǔ)量擬合及預(yù)估等方面具有明顯優(yōu)勢。然而,由于森林生態(tài)系統(tǒng)的自然地理要素及森林類型等的不同可能會(huì)造成不同區(qū)域相同森林生態(tài)系統(tǒng)或同一區(qū)域不同森林生態(tài)系統(tǒng)碳密度空間異質(zhì)性的差異,從而可能導(dǎo)致其碳密度估測最優(yōu)空間模型的不同,而目前多數(shù)學(xué)者主要集中在某一空間模型與傳統(tǒng)回歸模型在碳密度估測中的比較分析,較為缺乏對(duì)不同空間模型之間的比較研究,加強(qiáng)這方面的研究將有利于更加全面了解和掌握空間模型在森林生態(tài)系統(tǒng)碳密度估測方面的應(yīng)用。

馬尾松(Pinusmassoniana)是我國南方分布廣泛的鄉(xiāng)土針葉樹種,具有適應(yīng)性強(qiáng)、耐干旱與貧瘠的特點(diǎn)。江西省馬尾松林總面積為239.4萬hm2,占喬木林總面積的28.9%,是江西省主要的森林類型之一,在區(qū)域碳循環(huán)及應(yīng)對(duì)氣候變化中扮演著重要角色[16]。本研究以江西省馬尾松林生態(tài)系統(tǒng)為研究對(duì)象,分別利用最小二乘模型、空間誤差模型、空間滯后模型和地理加權(quán)回歸模型構(gòu)建其碳密度與影響因子之間的關(guān)系模型,比較分析選出最優(yōu)擬合模型,以豐富空間模型在森林碳密度估測中的應(yīng)用及為江西省馬尾松林碳匯林業(yè)的經(jīng)營規(guī)劃與管理提供參考依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

江西省位于我國東南部,地理坐標(biāo)為24°29′—30°04′N,113°34′—118°28′E,面積16.69×104km2。該區(qū)屬于亞熱帶溫暖濕潤季風(fēng)氣候,氣候溫和,雨量充沛,年均氣溫16.3—19.5℃,年均降雨量1675 mm。地貌以山地、丘陵為主,林地的土壤類型主要有紅壤、黃紅壤、黃壤等,成土母巖主要有花崗巖、頁巖、砂巖、板巖等。全省森林資源豐富,森林覆蓋率為63.1%,主要森林類型有常綠闊葉林、針葉林、針闊葉混交林、竹林等。

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

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

樣地?cái)?shù)據(jù)。根據(jù)研究區(qū)的森林分布圖及森林資源統(tǒng)計(jì)表,綜合考慮馬尾松林在全省的地域分布、林齡結(jié)構(gòu)、地形地貌特征等因素,于2011—2016年期間設(shè)置樣地調(diào)查,并收集2016年江西省森林資源連續(xù)清查馬尾松林樣地?cái)?shù)據(jù)(對(duì)部分缺失林下植被、凋落物等相關(guān)數(shù)據(jù)的樣地進(jìn)行補(bǔ)充調(diào)查)。樣地面積為900 m2和800 m2,樣地?cái)?shù)共計(jì)611個(gè),其樣地分布見圖1。喬木層每木調(diào)查起測胸徑為≥5 cm(胸徑<5 cm 的喬木視為灌木),測定其胸徑、樹高、林分密度等林分因子及地理坐標(biāo)、海拔、坡度、坡向、坡位、土層厚度、腐殖質(zhì)層厚度等立地因子,同時(shí)選取標(biāo)準(zhǔn)木取樣測定其干、枝、葉、根的碳含量。在樣地的上、中、下部分別布設(shè)2 m×2 m的灌木樣方,在所選灌木樣方中各設(shè)置1個(gè)1 m×1 m的草本樣方和凋落物樣方。采用“樣方收獲法”分別收集灌木(分枝、葉、根)、草本(分地上、地下部分)及凋落物并稱量鮮重,取樣用于測定其含水率和碳含量。在樣地內(nèi)選擇代表性地塊挖取100 cm深度(未達(dá)100 cm的挖至母巖)的土壤剖面,根據(jù)土層厚度按0—10 cm、10—20 cm、20—30 cm、30—50 cm、50—100 cm五個(gè)層次,對(duì)各層土壤采用環(huán)刀法(體積為100 cm3)取樣來測定土壤容重,并分別取約1 kg的土壤樣品帶回實(shí)驗(yàn)室用于土壤碳含量的測定。

圖1 樣地分布

氣象數(shù)據(jù)。本研究首先根據(jù)各樣地的地理坐標(biāo)從Wang等[17]編寫的提取亞太地區(qū)氣候數(shù)據(jù)的ClimateAP軟件中獲取1980—2016年期間各樣地每年的年降水量與年溫度數(shù)據(jù),然后計(jì)算得出各樣地的年均降水量與年均溫度。

1.2.2碳含量測定與碳密度計(jì)算

所有樣品研磨粉碎后過0.25 mm篩,采用重鉻酸鉀氧化-外加熱法來測定其碳含量[18]。喬木層、林下植被層及凋落物層碳密度為各組分生物量乘以對(duì)應(yīng)的碳含量,其中喬木層生物量采用異速生長方程W=a×Db,應(yīng)用我國立木生物量模型及碳計(jì)量參數(shù)系列行業(yè)標(biāo)準(zhǔn)計(jì)算[19-22]。土壤層碳密度的具體計(jì)算公式如下[16]：

(1)

式中,Dsoc為土壤層碳密度(t/hm2),Ci為土壤有機(jī)碳含量(g/kg),Di為土壤容重(g/cm3),Ei為土層厚度(cm),Gi為大于2 mm的石礫所占的體積百分比(%)。

喬木層、林下植被層、凋落物層與土壤層的碳密度之和為生態(tài)系統(tǒng)碳密度。

1.2.3模型選取與評(píng)價(jià)

模型選取。(1)最小二乘模型(OLS)：是建立在整體區(qū)域上的一種因變量偏差的平方和最小的模型估計(jì)；(2)空間誤差模型(SEM)：是指模型的誤差項(xiàng)導(dǎo)致了空間變量之間的相關(guān)性,變量之間的空間相互作用存在于誤差項(xiàng),研究的主要是相鄰樣地點(diǎn)的誤差沖擊對(duì)其他樣地點(diǎn)的影響；(3)空間滯后模型(SLM)：是指被解釋變量之間的空間依賴性對(duì)模型非常重要從而導(dǎo)致了空間相關(guān),它研究的是因變量在鄰接地區(qū)的行為對(duì)整個(gè)區(qū)域其他地區(qū)行為的影響；(4)地理加權(quán)回歸模型(GWR)：是一種局部參數(shù)估計(jì)的模型,其本質(zhì)是為數(shù)據(jù)集中的每一個(gè)要素都建立獨(dú)立的方程[23],通過將空間結(jié)構(gòu)嵌入線性回歸模型中,以探測空間關(guān)系的非平穩(wěn)性,該方法的估計(jì)結(jié)果有明確的解析表示,而且還能對(duì)其結(jié)果的參數(shù)估計(jì)進(jìn)行統(tǒng)計(jì)檢驗(yàn),其計(jì)算公式為：

(2)

其中,(yi;xi1,xi2, …,xip)為在地理位置(ui,vi)處的因變量yi和自變量xi1,xi2, …,xip的觀測值(i= 1, 2, …,n)。βk(ui,vi)(k=1, 2, …,p)為第i個(gè)觀測點(diǎn)(ui,vi)處的未知參數(shù),是(ui,vi)的任意函數(shù),?i(i= 1, 2, …,n)為獨(dú)立同分布的誤差項(xiàng),通常假定其服從N(0, σ2)。

模型評(píng)價(jià)。選用決定系數(shù)(R2)、均方誤差(MSE)和赤池信息量準(zhǔn)則(AIC)3個(gè)統(tǒng)計(jì)量對(duì)模型擬合效果進(jìn)行評(píng)價(jià)。R2主要是用來評(píng)價(jià)模型的擬合優(yōu)度的指標(biāo),其值越高則表明模型的擬合效果越好,即擬合出來的模型穩(wěn)定性越好；MSE主要是用來衡量因變量實(shí)測值與模型預(yù)估值之間偏差的指標(biāo),其值越低表明模型的預(yù)估結(jié)果越接近于實(shí)際測量值,即模型的擬合能力越高；AIC可以用來衡量模型的實(shí)用性和復(fù)雜性,其值越小表明模型的擬合效果越好。

模型殘差檢驗(yàn)。選擇Moran′s I指數(shù)檢驗(yàn)?zāi)Ｐ蜌埐畹目臻g自相關(guān)性,其計(jì)算公式如下：

(3)

在R軟件中,采用多元線性逐步回歸方法對(duì)生態(tài)系統(tǒng)碳密度的影響因子進(jìn)行篩選,并構(gòu)建最小二乘模型；利用Geoda軟件建立空間誤差模型和空間滯后模型；在GWR 4.0軟件中建立地理加權(quán)回歸模型。

2 結(jié)果與分析

2.1 影響因子篩選

根據(jù)森林碳密度研究的相關(guān)文獻(xiàn),采用頻度統(tǒng)計(jì)法及樣地調(diào)查內(nèi)容選取影響生態(tài)系統(tǒng)碳密度的因子,主要包括立地因子(海拔、坡度、坡向、坡位、土層厚度、腐殖質(zhì)厚度)、植被因子(林齡、胸徑、株數(shù)密度、樹高、郁閉度、林下植被蓋度、凋落物厚度)以及氣象因子(年均溫度、年均降水量)等15個(gè)因子。對(duì)坡向、坡位定性因子進(jìn)行賦值,將坡向劃分為陽坡、半陽坡、半陰坡、陰坡、無坡向,并依次賦值為1—5；將坡位劃分為山脊、上部、中部、下部、山谷、平地、全坡,依次賦值為1—7。利用多元線性逐步回歸方法對(duì)各因子與生態(tài)系統(tǒng)碳密度進(jìn)行分析,從中篩選出與碳密度相關(guān)性顯著且各自變量之間的方差膨脹因子均小于5(即不存在多重共線性)的因子,用于數(shù)據(jù)建模。篩選出的因子分別為：海拔、坡度、土層厚度、胸徑、年均溫度、年均降水量,其統(tǒng)計(jì)量見表1。

表1 碳密度及各因子基本統(tǒng)計(jì)量

2.2 模型構(gòu)建

2.2.1最小二乘模型、空間誤差模型與空間滯后模型

最小二乘模型、空間誤差模型與空間滯后模型的參數(shù)估計(jì)結(jié)果見表2,坡度與碳密度的系數(shù)為負(fù)值,這表明碳密度隨著坡度的增大呈減少的趨勢；海拔、土層厚度、胸徑、年均溫度、年均降水量與碳密度的系數(shù)均為正值,表明碳密度與這些因子之間存在正相關(guān),其中系數(shù)最大的為胸徑,最小的為年均降水量。海拔、坡度、土層厚度和胸徑對(duì)碳密度的影響均達(dá)顯著性水平(P<0.05)。

表2 最小二乘、空間誤差及空間滯后模型的擬合結(jié)果

2.2.2地理加權(quán)回歸模型

地理加權(quán)回歸模型的參數(shù)估計(jì)結(jié)果見表3,海拔、土層厚度、胸徑與碳密度的系數(shù)平均值分別為0.0817、0.4616、5.0531,其第1四分位數(shù)至最大值的范圍分別為0.0074—0.2920、0.3081—1.1967、3.5514—10.3628,可看出這些因子75%的系數(shù)均為正值,說明碳密度與海拔、土層厚度、胸徑3個(gè)因子之間存在正相關(guān)關(guān)系；坡度與碳密度的系數(shù)平均值為-2.4527,其最小值至第3四分位數(shù)的范圍為-4.4625至-1.8330,表明坡度與碳密度之間為負(fù)相關(guān)關(guān)系,即碳密度隨著坡度的增大呈逐漸減少的趨勢。年均溫度、年均降水量與碳密度的系數(shù)平均值為1.2671和0.0237,但50%的系數(shù)為正值,另外50%部分為負(fù)值,說明氣象因子對(duì)碳密度的影響比較多變。

表3 地理加權(quán)回歸模型的系數(shù)估計(jì)值

從模型參數(shù)估計(jì)結(jié)果可知,4種模型均表明碳密度與坡度呈負(fù)相關(guān),與海拔、土層厚度、胸徑呈正相關(guān)。最小二乘模型、空間誤差模型和空間滯后模型的結(jié)果均表明海拔、坡度、土層厚度、胸徑對(duì)碳密度有顯著影響(P<0.05),年均溫度和年均降水量對(duì)碳密度的影響并不顯著(P>0.05),而地理加權(quán)回歸模型的結(jié)果則反映出氣象因子對(duì)碳密度的影響比較多變。

2.3 模型評(píng)價(jià)與檢驗(yàn)

2.3.1模型評(píng)價(jià)

據(jù)表4可得,4種模型的決定系數(shù)(R2)最大的為地理加權(quán)回歸模型,其次為空間誤差模型與空間滯后模型,最小的為最小二乘模型；地理加權(quán)回歸模型的均方誤差(MSE)與赤池信息準(zhǔn)則(AIC)值均為4種模型中最小,而最小二乘模型的AIC與MSE值為最大。由此可得,模型擬合效果由高到低依次為地理加權(quán)回歸模型>空間誤差模型>空間滯后模型>最小二乘模型。

表4 模型評(píng)價(jià)結(jié)果

圖2為碳密度實(shí)測值與4種模型估計(jì)值之間的散點(diǎn)圖。4種模型的估計(jì)值與碳密度實(shí)測值的配對(duì)點(diǎn)均離散分布在1:1的參考線附近,當(dāng)碳密度實(shí)測值較小時(shí),最小二乘模型、空間誤差模型及空間滯后模型擬合得出的估計(jì)值普遍偏高,而地理加權(quán)回歸模型的估計(jì)值比其他3種模型更為接近實(shí)測值；當(dāng)碳密度實(shí)測值較大時(shí),最小二乘模型、空間誤差模型及空間滯后模型擬合得出的估計(jì)值均普遍偏低,但地理加權(quán)回歸模型的估計(jì)值與實(shí)測值偏低的程度明顯小于其他3種模型。由此可得,地理加權(quán)回歸模型的擬合結(jié)果優(yōu)于其他3種模型。

圖2 碳密度實(shí)測值與模型估計(jì)值散點(diǎn)圖

2.3.2模型殘差的空間自相關(guān)性

由表5可得,最小二乘模型、空間誤差模型和空間滯后模型的Moran′s I均為正值,且P<0.05；地理加權(quán)回歸模型的Moran′s I為負(fù)值,且P>0.05。由此得出,最小二乘模型、空間誤差模型和空間滯后模型的模型殘差均存在顯著的空間自相關(guān)性,而地理加權(quán)回歸模型的模型殘差的空間自相關(guān)性不顯著,這表明地理加權(quán)回歸模型可有效降低模型殘差的空間自相關(guān)性。

表5 模型殘差的Moran′s I及檢驗(yàn)

4種模型殘差的Moran′s I隨距離的變化見圖3。最小二乘模型、空間誤差模型和空間滯后模型殘差的Moran′s I均為正值,即3種模型殘差的空間相關(guān)性為正相關(guān),且模型殘差的Moran′s I在帶寬較小時(shí),其降低幅度均較大,隨著距離的增加其降低幅度逐漸減弱,逐漸趨近于0；地理加權(quán)回歸模型殘差的Moran′s I為負(fù)值,即其模型殘差的空間相關(guān)性為負(fù)相關(guān),其模型殘差的Moran′s I在帶寬較小時(shí),其增加幅度較大,隨著距離的增大其增加幅度逐漸減弱,也逐漸趨近于0。模型殘差的Moran′s I由大到小表現(xiàn)為最小二乘模型>空間滯后模型>空間誤差模型>地理加權(quán)回歸模型。由此也可得出,與其他模型相比,地理加權(quán)回歸模型有效降低了模型殘差的空間自相關(guān)性。

圖3 模型殘差Moran′s I

綜上分析,4種模型中,地理加權(quán)回歸模型決定系數(shù)(R2)最大,均方誤差(MSE)與赤池信息準(zhǔn)則(AIC)均最小,且能有效降低模型殘差的空間自相關(guān)性,表明地理加權(quán)回歸模型的擬合效果優(yōu)于其他3種模型。

3 討論與結(jié)論

3.1 討論

由于森林生態(tài)系統(tǒng)具有較強(qiáng)的空間異質(zhì)性,所以,不少學(xué)者利用空間模型研究森林碳密度/碳儲(chǔ)量/生物量時(shí)均得出空間模型比最小二乘模型(OLS)的擬合效果更好的結(jié)論。如李月等[11]采用空間誤差模型(SEM)構(gòu)建碳儲(chǔ)量與影響因素關(guān)系時(shí)得出其擬合精度明顯優(yōu)于OLS模型,其主要是由于SEM減小了模型殘差的空間自相關(guān)性；劉暢等[24]利用空間滯后模型(SLM)預(yù)估林木含碳量得出SLM優(yōu)于OLS模型,認(rèn)為SLM模型能更好地解決模型構(gòu)建中的空間效應(yīng),其穩(wěn)定性更高,對(duì)各因子的解釋能力更強(qiáng)。

在空間模型的運(yùn)用中,更多學(xué)者采用了地理加權(quán)回歸模型(GWR),且其研究結(jié)果幾乎都表明GWR優(yōu)于其他模型。如Liu等[25]利用OLS模型、線性混合模型和GWR模型擬合黑龍江省喬木層碳儲(chǔ)量與影響因子間的關(guān)系時(shí)表明了GWR模型在擬合精度以及在預(yù)估方面都明顯優(yōu)于其他兩種模型；Lin等[13]利用OLS模型和GWR模型擬合福建省將樂縣喬木層碳密度與影響因子之間的關(guān)系,結(jié)果也表明了GWR模型在擬合精度及降低模型殘差自相關(guān)性方面都明顯優(yōu)于OLS模型。此外,歐光龍等[14]、Zhang等[26]的研究均表明GWR模型為最優(yōu)模型。而本研究通過對(duì)OLS、SEM、SLM和GWR 4種模型的模型評(píng)價(jià)指標(biāo)及殘差檢驗(yàn)的綜合分析,得出OLS模型的擬合效果最差,GWR模型最好。這可能主要是因?yàn)樵趨^(qū)域尺度上,森林碳密度本身存在明顯空間效應(yīng),而OLS模型是以忽略數(shù)據(jù)的空間效應(yīng)作為前提下的無偏估計(jì)；SEM、SLM模型與OLS模型相比,雖然減小了模型殘差的空間自相關(guān)性,提高了擬合效果,但由于沒有直接將這種空間自相關(guān)性納入到擬合過程中,難于解決空間異質(zhì)性問題[27]；而GWR模型適用于空間非平穩(wěn)性的研究[28],它在對(duì)數(shù)據(jù)集所有位置點(diǎn)的參數(shù)進(jìn)行估計(jì)時(shí)考慮了其空間非平穩(wěn)性,針對(duì)每一個(gè)坐標(biāo)位置點(diǎn)都有相對(duì)應(yīng)的參數(shù)[29],將模型殘差的空間自相關(guān)性直接納入到擬合過程中,有效降低了模型殘差的空間自相關(guān)性,所以其擬合效果最優(yōu)。

另外,一些學(xué)者以森林碳儲(chǔ)量/生物量為因變量,以遙感影像數(shù)據(jù)為自變量構(gòu)建空間模型時(shí)同樣也得出了GWR模型優(yōu)于其他空間模型的結(jié)論。如戚玉嬌[30]通過構(gòu)建黑龍江省大興安嶺地區(qū)的森林植被地上碳儲(chǔ)量與TM遙感影像數(shù)據(jù)之間的關(guān)系表明GWR模型的擬合效果最優(yōu)；郭含茹等[15]利用GWR模型、傳統(tǒng)全局回歸模型、協(xié)同克里格插值構(gòu)建浙江省仙臺(tái)縣森林地上部分碳儲(chǔ)量和Landsat TM遙感影像數(shù)據(jù)之間的關(guān)系也表明GWR模型是有效的,閭妍宇等[31]、Wang等[32]的研究也得出同樣的結(jié)論?？梢?GWR模型能有效解決森林生態(tài)系統(tǒng)空間異質(zhì)性問題。

在區(qū)域尺度上,對(duì)森林碳密度/碳儲(chǔ)量/生物量空間模型構(gòu)建時(shí),多數(shù)學(xué)者僅采用某一空間模型與傳統(tǒng)回歸模型進(jìn)行比較,而本研究通過不同空間模型比較分析,發(fā)現(xiàn)GWR模型擬合效果優(yōu)于其他空間模型。然而,基于GWR模型的擴(kuò)展模型的相應(yīng)出現(xiàn),如地理加權(quán)回歸Kriging模型、地理加權(quán)回歸Logistic模型、地理加權(quán)回歸泊松模型等,而這些擴(kuò)展模型若用于本研究中其擬合精度是否會(huì)提高,還有待今后進(jìn)一步深入探討。

3.2 結(jié)論

本研究以江西省為研究區(qū),通過多元線性逐步回歸方法篩選出與馬尾松林生態(tài)系統(tǒng)碳密度相關(guān)性顯著且不存在多重共線性的影響因子為海拔、坡度、土層厚度、胸徑、年均溫度和年均降水量。采用OLS、SEM、SLM和GWR模型構(gòu)建生態(tài)系統(tǒng)碳密度與其影響因子之間的關(guān)系模型,通過對(duì)模型擬合效果的評(píng)價(jià)及殘差檢驗(yàn)的綜合分析,得出GWR模型擬合效果最優(yōu),更適用于江西省馬尾松林生態(tài)系統(tǒng)碳密度的估測。