基于不同空間模型的馬尾松林生態系統碳密度研究

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

1 北京林業大學林學院,北京 100083 2 江西農業大學林學院,南昌 330045 3 鄱陽湖流域森林生態系統保護與修復國家林業和草原局重點實驗室,南昌 330045

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

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

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

1 材料與方法

1.1 研究區概況

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

1.2 數據與方法

1.2.1數據來源

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

圖1 樣地分布

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

1.2.2碳含量測定與碳密度計算

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

(1)

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

喬木層、林下植被層、凋落物層與土壤層的碳密度之和為生態系統碳密度。

1.2.3模型選取與評價

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

(2)

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

模型評價。選用決定系數(R2)、均方誤差(MSE)和赤池信息量準則(AIC)3個統計量對模型擬合效果進行評價。R2主要是用來評價模型的擬合優度的指標,其值越高則表明模型的擬合效果越好,即擬合出來的模型穩定性越好；MSE主要是用來衡量因變量實測值與模型預估值之間偏差的指標,其值越低表明模型的預估結果越接近于實際測量值,即模型的擬合能力越高；AIC可以用來衡量模型的實用性和復雜性,其值越小表明模型的擬合效果越好。

模型殘差檢驗。選擇Moran′s I指數檢驗模型殘差的空間自相關性,其計算公式如下：

(3)

在R軟件中,采用多元線性逐步回歸方法對生態系統碳密度的影響因子進行篩選,并構建最小二乘模型；利用Geoda軟件建立空間誤差模型和空間滯后模型；在GWR 4.0軟件中建立地理加權回歸模型。

2 結果與分析

2.1 影響因子篩選

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

表1 碳密度及各因子基本統計量

2.2 模型構建

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

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

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

2.2.2地理加權回歸模型

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

表3 地理加權回歸模型的系數估計值

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

2.3 模型評價與檢驗

2.3.1模型評價

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

表4 模型評價結果

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

圖2 碳密度實測值與模型估計值散點圖

2.3.2模型殘差的空間自相關性

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

表5 模型殘差的Moran′s I及檢驗

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

圖3 模型殘差Moran′s I

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

3 討論與結論

3.1 討論

由于森林生態系統具有較強的空間異質性,所以,不少學者利用空間模型研究森林碳密度/碳儲量/生物量時均得出空間模型比最小二乘模型(OLS)的擬合效果更好的結論。如李月等[11]采用空間誤差模型(SEM)構建碳儲量與影響因素關系時得出其擬合精度明顯優于OLS模型,其主要是由于SEM減小了模型殘差的空間自相關性；劉暢等[24]利用空間滯后模型(SLM)預估林木含碳量得出SLM優于OLS模型,認為SLM模型能更好地解決模型構建中的空間效應,其穩定性更高,對各因子的解釋能力更強。

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

另外,一些學者以森林碳儲量/生物量為因變量,以遙感影像數據為自變量構建空間模型時同樣也得出了GWR模型優于其他空間模型的結論。如戚玉嬌[30]通過構建黑龍江省大興安嶺地區的森林植被地上碳儲量與TM遙感影像數據之間的關系表明GWR模型的擬合效果最優；郭含茹等[15]利用GWR模型、傳統全局回歸模型、協同克里格插值構建浙江省仙臺縣森林地上部分碳儲量和Landsat TM遙感影像數據之間的關系也表明GWR模型是有效的,閭妍宇等[31]、Wang等[32]的研究也得出同樣的結論。可見,GWR模型能有效解決森林生態系統空間異質性問題。

在區域尺度上,對森林碳密度/碳儲量/生物量空間模型構建時,多數學者僅采用某一空間模型與傳統回歸模型進行比較,而本研究通過不同空間模型比較分析,發現GWR模型擬合效果優于其他空間模型。然而,基于GWR模型的擴展模型的相應出現,如地理加權回歸Kriging模型、地理加權回歸Logistic模型、地理加權回歸泊松模型等,而這些擴展模型若用于本研究中其擬合精度是否會提高,還有待今后進一步深入探討。

3.2 結論

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