李成榮，王 妍，張 超，劉云根，張詩文，鄧再春，錢 慧

（1.西南林業(yè)大學(xué) 林學(xué)院，昆明 650224；2.云南省山地農(nóng)村生態(tài)環(huán)境演變與污染治理重點實驗室，昆明 650224；3.西南林業(yè)大學(xué) 生態(tài)與環(huán)境學(xué)院，昆明 650224）

森林碳儲量是森林生態(tài)系統(tǒng)各碳庫中碳元素的儲備量，其主要通過森林植被吸收并固定中CO2于植被或土壤中實現(xiàn)，是森林生態(tài)系統(tǒng)結(jié)構(gòu)和森林質(zhì)量評價的重要指標之一，亦是衡量區(qū)域生態(tài)環(huán)境優(yōu)劣的重要指標之一。森林生態(tài)系統(tǒng)蘊含著巨大的碳儲量，與全球氣候變化、全球碳循環(huán)存在緊密聯(lián)系［1-2］。我國西南地區(qū)是世界上喀斯特地貌分布面積最大、最集中的地區(qū)［3］。石漠化的發(fā)生在導(dǎo)致生態(tài)功能下降的同時，使得區(qū)域碳匯能力下降，因此，石漠化地區(qū)的植被恢復(fù)應(yīng)受到廣泛關(guān)注。準確、快速地估算喀斯特地區(qū)的森林碳儲量，了解其分布與變化過程，在為區(qū)域石漠化治理、生態(tài)修復(fù)提供數(shù)據(jù)支撐的同時，可促進碳儲量計量方法的完善，精準挖掘和評價森林碳匯潛力，助力森林碳匯發(fā)展。

利用遙感手段進行區(qū)域森林碳儲量的估算是基于植物的反射光譜特征實現(xiàn)的。建立光譜信息與森林碳儲量之間的數(shù)學(xué)擬合模型已成為森林碳儲量遙感反演研究的熱點。目前，有關(guān)森林碳儲量的預(yù)測研究受不同尺度、不同方法、不同森林類型及不同生境條件的影響，總體存在較大不確定性，對森林碳儲量的預(yù)測仍缺少合理有效的方法［4］?；谶b感的森林碳儲量反演方法包括參數(shù)和非參數(shù)經(jīng)驗?zāi)Ｐ头椒?，其中，參?shù)模型包括線性、對數(shù)或冪函數(shù)等［5］形式，非參數(shù)模型包括神經(jīng)網(wǎng)絡(luò)［6］、k最近鄰域法［7］、隨機森林［8］和支持向量機［9］等，但都存在參數(shù)復(fù)雜、適用性差等缺陷，在估算森林碳儲量上有著不同的效果。王璟睿［10］通過構(gòu)建隨機森林、BP 神經(jīng)網(wǎng)絡(luò)、支持向量機、協(xié)同克里金法與地理加權(quán)回歸的方法，對廣東省森林碳密度進行研究，得出隨機森林效果最優(yōu)。梯度提升回歸樹（gradient boosting regression tree，GBRT）是一種非參數(shù)建模方法，它的基本原理是通過構(gòu)建多個弱分類器，經(jīng)過多次迭代來構(gòu)建模型［11］，GBRT 能靈活處理那些復(fù)雜的非線性關(guān)系，即使數(shù)據(jù)樣本類型不同，也能在大樣本數(shù)據(jù)預(yù)測中發(fā)揮很好的效果［12］，擁有運行速度更快、精度高等優(yōu)點，已在葉綠素含量反演［13］、NDVI 預(yù)測［14］、生物量估算［15］與樹高估算［16］等方面得到應(yīng)用。然而，利用GBRT模型對區(qū)域尺度森林碳儲量預(yù)測研究還鮮見報道。森林資源連續(xù)清查數(shù)據(jù)具有涉及范圍廣、森林類型全、調(diào)查對象豐富、測定因子易獲得、時間連續(xù)性強等優(yōu)點，涂宏濤等［17］基于云南省第九次森林資源清查數(shù)據(jù)，對碳儲量分布特征進行了研究。

鑒于此，本研究以滇東地區(qū)為研究區(qū)，基于1997 年、2002 年、2007 年森林資源連續(xù)清查數(shù)據(jù)與同期遙感影像數(shù)據(jù)，建立遙感因子與森林碳儲量之間的偏最小二乘模型、KNN 模型、隨機森林模型與GBRT 模型，通過模型檢驗（R2、RMSE、rRMSE、P）選擇最優(yōu)模型進行森林碳儲量估算，結(jié)合空間自相關(guān)分析法，揭示了滇東地區(qū)1997—2022 年25 年間森林碳儲量變化特征，為森林可持續(xù)經(jīng)營與區(qū)域碳匯提升提供技術(shù)參考。

1 研究區(qū)概況

選取我國石漠化分布的典型區(qū)域——滇東地區(qū)（包括昆明市、曲靖市、玉溪市、紅河州和文山州）作為研究區(qū)，介于22o26′～27o03′N、101o16′～106o12′E 之間，土地總面積1.293×105km2，平均海拔1 651 m，屬亞熱帶高原季風(fēng)性氣候，年均降水量為1 094 mm。該區(qū)發(fā)育著各種類型的喀斯特地貌，屬于云南乃至整個西南地區(qū)石漠化最為嚴重的區(qū)域。地形以起伏的低山和丘陵為主。研究區(qū)內(nèi)天然林占絕對優(yōu)勢，植被類型主要為暖性針葉林、暖性針闊混交林、落葉闊葉林、常綠闊葉林、常綠落葉混交林和灌草叢等。

圖1 研究區(qū)地理位置及樣地分布Fig.1 Location of the study area and distribution of the plots

2 研究方法

2.1 數(shù)據(jù)收集與處理

2.1.1 樣地數(shù)據(jù)收集及森林碳儲量計算

收集位于滇東地區(qū)的第5 次（1997 年）、第6次（2002 年）和第7 次（2007 年）森林資源連續(xù)清查樣地數(shù)據(jù)。云南省森林資源連續(xù)清查采用系統(tǒng)抽樣方法，抽樣間距為6 km×8 km，方形樣地面積為0.08 hm2。以滇東地區(qū)主要的優(yōu)勢樹種（樹種組）為對象，共篩選出430個樣地（表1），進行森林生物量與森林碳儲量的計算。生物量計算采用蓄積量-生物量模型，碳儲量計算采用生物量-碳儲量模型，公式如下：

表1 滇東地區(qū)主要優(yōu)勢樹種模型轉(zhuǎn)換參數(shù)Table 1 Model transformation parameters of the main dominant tree species in eastern Yunnan

式中：B為森林生物量，t/hm2；V為森林蓄積量，m3/hm2；a、b為蓄積量與生物量之間的模型轉(zhuǎn)換參數(shù)；C為碳儲量，t/hm2；Cc為含碳系數(shù)。各優(yōu)勢樹種及其轉(zhuǎn)換參數(shù)如表1所示［18-19］。

2.1.2 遙感影像處理

研究采用的遙感影像來源于谷歌地球引擎（google earth engine，GEE）的1997—2022年Landsat TM/OLI 數(shù)據(jù)，每5年為一個時段，成像時間為5—10 月，空間分辨率為30 m×30 m。該數(shù)據(jù)集已經(jīng)過大氣校正、幾何校正和地形校正等預(yù)處理，采用去云算法對影像進行去云處理，對去云后的像元進行時間插值處理，采用中值算法合成影像［20］。DEM 數(shù)據(jù)來源于地理空間數(shù)據(jù)云，空間分辨率為30 m×30 m。

2.2 建模因子提取

Landsat TM 影像去除Band 6 TIR 波段，Landsat OLI 影像去除Band 1 Coastal 波段，以Landsat TM 影像為例，提取的遙感因子結(jié)果如下，各因子計算公式見參考文獻［21］：

（1）原始單波段因子，包括TM1（遙感影像第1 波段）、TM2（遙感影像第2 波段）、TM3（遙感影像第3 波段）、TM4（遙感影像第4 波段）、TM5（遙感影像第5波段）、TM7（遙感影像第7波段）。

（2）地形因子，包括坡度、坡向、高程。

（3）紋理特征因子，包括均值（mean）、方差（variance）、均一性（homogeneity）、對比度（contrast）、相關(guān)性（correlation）、相異性（dissimilarity）、角二矩（angular second moment）、熵（entropy），釆用基于灰度共生矩陣的方法，窗口大小為7×7，每個原始單波段因子對應(yīng)8個紋理特征值。

（4）植被指數(shù)因子，包括歸一化植被指數(shù)（normalized difference vegetation index，NDVI）、比值植被指數(shù)（ratio vegetation index，RVI）、增強型植被指數(shù)（enhanced vegetation index，EVI）、土壤調(diào)節(jié)植被指數(shù)（soil adjusted vegetation index，SAVI）、大氣阻抗植被指數(shù)（atmospherically resistant vegetation index，ARVI）、差值植被指數(shù)（difference vegetation index，DVI）、非線性指數(shù)（nonlinear vegetation index，NLI）、垂直植被指數(shù)（perpendicular vegetation index，PVI）、土壤調(diào)整比值植被指數(shù)（soil adjusted ratio vegetation，SARV）、修正簡單比值指數(shù)（modified simple ratio，MSR）、轉(zhuǎn)換型植被指數(shù)（transformed vegetation index，TVI）、綠光歸一化差值植被指數(shù)（green normalized difference vegetation index，GNDVI）、改良土壤調(diào)整植被指數(shù)（soil adjusted vegetation index，SAVI）、修正型土壤調(diào)節(jié)植被指數(shù)（modified soil adjusted vegetation index，MSAVI）、波段組合因子（VIS123、Albedo、TM3/Albedo）、歸一化差值含水指數(shù)（normalized difference moisture index，NDMI）、可見光大氣阻被指數(shù)（visible atmospherically resistant index，VARI）、改進的歸一化差值水分指數(shù)（modified normalized difference water index，MNDWI）。

（5）主成分變換因子，包括第一主成分（PCA1）、第二主成分（PCA2）、第三主成分（PCA3）。

（6）纓帽變換因子，包括亮度（BVI）、綠度（GVI）、濕度（WVI）。

2.3 估算模型構(gòu)建

（1）偏最小二乘回歸（partial least-squares regression，簡稱PLS回歸），該模型應(yīng)用了主成分分析的基本思想，并考慮了主成分與因變量之間的相關(guān)性，本研究利用Minitab構(gòu)建偏最小二乘模型。

（2）KNN（k-nearest neighbor）法，又稱基準樣本法［7］，是一種根據(jù)最鄰近的一個或多個樣本屬性共同決定未知樣本屬性的決策分類方法，運用Python 編程軟件構(gòu)建KNN 模型，模型訓(xùn)練前需要進行調(diào)整的參數(shù)包括鄰居數(shù)n_neighbors、權(quán)重weights、度量metric以及算法algorithm。

（3）隨機森林（random forest，RF），是一種基于多個分類樹的集成學(xué)習(xí)算法［8］，運用Python編程軟件構(gòu)建隨機森林模型，模型訓(xùn)練前需要進行調(diào)整的參數(shù)包括學(xué)習(xí)器數(shù)量n_estimators、決策樹最大深度max_depth、內(nèi)部節(jié)點分裂所需的最小樣本數(shù)min_samples_split 以及葉節(jié)點最小樣本數(shù)min_samples_leaf。

（4）梯度提升回歸樹（gradient boosted regression trees，GBRT）［11-12］，與隨機森林相似，是一種集成方法，通過合并多個決策樹，經(jīng)過多次迭代來構(gòu)建模型，本研究基于Python 編程軟件實現(xiàn)。具體，對于一個回歸問題y=F(x)，以x∈Rd輸入與y∈R輸出組成，給定訓(xùn)練集數(shù)據(jù)Dtrain:={(x1,y1),…,(xN,yN)}和損失函數(shù)L[y,F(x) ]，完成算法迭代訓(xùn)練后可表示為：

式中：F(x)為初始模型預(yù)測值；M為迭代次數(shù)；為第次迭代的弱學(xué)習(xí)器數(shù)量，其通過優(yōu)化目標函數(shù)所得；α為學(xué)習(xí)率；J為第m次迭代的回歸樹中葉子節(jié)點個數(shù)；fm為GBRT中單棵決策樹的模型；Cmj為回歸樹中葉子節(jié)點的最優(yōu)輸出值；Rmj為葉子節(jié)點所屬特征空間，當(dāng)x與Rmj屬于同一特征子空間，。

模型訓(xùn)練前需要進行調(diào)整的參數(shù)包括學(xué)習(xí)器數(shù)量n_estimators、學(xué)習(xí)率learning_rate、子樣本subsample、最大深度max_depth、葉節(jié)點最小樣本數(shù)min_samples_leaf。輸出要素包括交叉驗證結(jié)果、模型決定系數(shù)R2、均方根誤差RMSE、相對均方根誤差rRMSE、總體精度P及特征重要性。

2.4 檢驗與制圖方法

模型檢驗指標：模型決定系數(shù)（R2），均方根誤差（RMSE），相對均方根誤差（rRMSE），總體精度（P）［22］。

制圖分析方法：經(jīng)驗貝葉斯克里金法（empirical bayesian kring，EBK）［23］、Moran’s I 指數(shù)法［24］、Getis-Ord Gi*統(tǒng)計量［25］。

3 結(jié)果與分析

3.1 森林碳儲量估算模型構(gòu)建及檢驗

基于計算所得的森林碳儲量與提取的建模因子，利用SPSS進行Pearson相關(guān)性分析，最終篩選出與森林碳儲量極顯著相關(guān)的遙感因子，共66個，參與四種碳儲量估算模型的構(gòu)建，相關(guān)性分析結(jié)果如表2所示。

表2 Pearson相關(guān)性分析Table 2 Pearson correlation analysis

（續(xù)表2）

利用表2 篩選的遙感因子進行GBRT 模型的構(gòu)建，采用隨機抽樣的方法，以其中80%的樣地數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，以其余20%的樣地數(shù)據(jù)作為檢驗數(shù)據(jù)。在模型學(xué)習(xí)之前，設(shè)置學(xué)習(xí)器數(shù)量（n_estimators）、學(xué)習(xí)率（learning_rate）、子樣本（subsample）、最大深度（max_depth）、最小樣本數(shù)（min_samples_leaf）。模型參數(shù)的確定需遵循經(jīng)驗。學(xué)習(xí)器數(shù)量的增加會提升擬合能力，但同時亦會提高過擬合概率，通常在預(yù)測值收斂的情況下越小越好；學(xué)習(xí)率一般小于0.1；最大深度一般不超過15；子樣本與最小樣本數(shù)需根據(jù)樣本量的大小調(diào)整，樣本量較多時可適當(dāng)增加，避免出現(xiàn)過擬合。模型訓(xùn)練完成，根據(jù)輸出的特征重要性（相關(guān)性）結(jié)果，進行二次訓(xùn)練，最終確定的建模因子如表3所示。GBRT模型經(jīng)過20次交叉驗證篩選的最優(yōu)參數(shù)組合為：n_estimators：110；learning_rate：0.05；subsample：0.5；max_depth：8；min_samples_leaf：1。

表3 特征重要性（相關(guān)性）Table 3 Feature importance（correlation）

4 種方法精度驗證結(jié)果如表4所示，其中，GBRT模型的R2最高，為0.94；RMSE最低，為3.98 t/hm2；rRMSE最低，為13.44 t/hm2；P最高，為82.57%；其次為隨機森林模型，偏最小二乘回歸模型表現(xiàn)最差。綜上，選擇效果最優(yōu)的GBRT模型進行森林碳儲量遙感估算與制圖分析。以1997年為例，采用1 km×1 km漁網(wǎng)處理，共設(shè)有156 806個樣點，進行研究區(qū)森林碳儲量的估算，最后采用EBK法進行研究區(qū)森林碳儲量遙感制圖。

表4 估算模型精度評價Table 4 Accuracy evaluation of the estimation models

3.2 滇東地區(qū)森林碳儲量變化特征分析

分析滇東地區(qū)森林碳儲量動態(tài)變化（圖2）可知，1997—2022 年滇東地區(qū)碳密度呈先減少后增加再減少的趨勢，總體由1997年的21.5 t/hm2增長為2022 年的22.33 t/hm2，增幅為3.9%。1997 年滇東地區(qū)森林碳儲量分布面積占比最大為10～20 t/hm2，其主要分布于曲靖市、昆明市的東部、玉溪市的東部、紅河州的北部和文山州的西部，到2022 年該面積占比呈大幅度減少，主要向著20～30 t/hm2范圍轉(zhuǎn)移。1997年森林碳儲量在30 t/hm2以上的區(qū)域主要分布于昆明市的西南部、玉溪市的西南部與中部、紅河州的南部、文山州的東部與北部，到2022 年該范圍占比除曲靖市和文山州呈增加趨勢外，昆明市、玉溪市、紅河州呈減少趨勢。

圖2 滇東地區(qū)1997—2022年森林碳儲量動態(tài)變化Fig.2 Dynamic changes of forest carbon storage in eastern Yunnan from 1997 to 2022

分析滇東地區(qū)森林碳儲量的時間變化過程（圖3）可知，1997—2022年滇東地區(qū)森林總碳儲量呈先降低后升高再降低的變化趨勢，最低值出現(xiàn)于2002年，最高值出現(xiàn)于2017年。1997—2022年森林總碳儲量總體呈上升趨勢，增幅為4.2%；1997—2002 年呈下降趨勢，降幅為1.8%；2002—2017 年呈上升趨勢，增幅為7.9%；2017—2022 年呈下降趨勢，降幅為1.7%。1997—2022年昆明市、曲靖市、紅河州和文山州森林總碳儲量總體均呈上升趨勢，增幅分別為4.7%、9.7%、0.26%和13.3%；玉溪市則呈下降趨勢，降幅為12.6%，主要由于植被覆蓋退化所導(dǎo)致［26］。20世紀90年代以來，隨著人口增多導(dǎo)致土地依附性增強，再加上不合理的開墾耕作方式、能源獲取以及對森林的亂砍濫伐［27］，1997—2002 年滇東地區(qū)森林總碳儲量呈現(xiàn)下降趨勢。21 世紀以來，隨著退耕還林還草、集體林權(quán)制度改革以及石漠化綜合治理等生態(tài)工程的有效實施［27］，2002 年以后，森林總碳儲量開始呈現(xiàn)增長趨勢。

圖3 滇東地區(qū)森林碳儲量時間變化Fig.3 Variation of carbon storage in eastern Yunnan with time

3.3 滇東地區(qū)森林碳儲量空間異質(zhì)性特征分析

空間自相關(guān)結(jié)果顯示，1997 年、2002 年、2007 年、2012 年、2017 年、2022 年滇東地區(qū)的全局莫蘭指數(shù)分別為0.361、0.365、0.347、0.367、0.288 和0.260，Z得分分別為54.424、55.129、52.439、55.338、43.477 和39.247，P值均為0，說明有99%的概率認為森林碳儲量不是隨機分布的，隨機分布的可能性小于1%，即森林碳儲量在相鄰空間單元上呈聚集狀態(tài)，區(qū)域森林碳儲量具有空間異質(zhì)性；從時間角度出發(fā)，1997—2012 年全局莫蘭指數(shù)趨于平穩(wěn)，由于2009—2011 年連續(xù)三年干旱，全局莫蘭指數(shù)自2012年起開始呈下降趨勢。

分析滇東地區(qū)空間聚類動態(tài)演變過程（圖4）可知，在聚集類型上，面積占比最大為99%置信度下的冷熱點區(qū)域，其主要分布在不顯著過渡地帶邊緣。在聚集特征上，冷點區(qū)域與熱點區(qū)域的空間分布一般是不相鄰的，且通過不顯著空間類別進行過渡，說明森林碳儲量外在分布上一般是漸變式的，而內(nèi)在呈現(xiàn)高值（或低值）的聚集狀態(tài)。在地理位置上，1997 年滇東地區(qū)森林碳儲量高值聚類主要分布在曲靖市的東南部、昆明市的西南與西北部、紅河州的南部、文山州的廣南縣與富寧縣，玉溪市則有大面積分布，到2022 年，森林碳儲量高值聚類在玉溪市呈現(xiàn)大幅度減少，主要在峨山縣、新平縣和元江縣，文山州則呈現(xiàn)增長趨勢，主要在南部片區(qū)，其他各地州均有小幅度減少；碳儲量低值聚類主要分布在昆明市的東部、紅河州的中部、文山州的西部，曲靖市有大面積分布，到2022 年森林碳儲量低值聚類面積在曲靖市與昆明市呈減少趨勢，其他各地州則趨于平穩(wěn)。

圖4 滇東地區(qū)1997—2022年空間聚類動態(tài)演變特征Fig.4 Dynamic evolution characteristics of spatial clustering in eastern Yunnan from 1997 to 2022

4 討論與結(jié)論

4.1 討論

本文基于樣地數(shù)據(jù)與遙感影像數(shù)據(jù)，運用四種模型對滇東森林碳儲量進行估算，結(jié)果表明GBRT模型估算結(jié)果較好，但存在以下不足。一是技術(shù)方法方面：（1）基于樣地清查數(shù)據(jù)采用生物量法估算森林植被碳儲量時，由于不同類型優(yōu)勢樹種參數(shù)存在差異導(dǎo)致低于實際值，此外，不同優(yōu)勢樹種含碳率也是一個重要指標，構(gòu)建科學(xué)準確的指標系統(tǒng)，勢在必行。（2）本研究使用的Landsat TM/OLI影像數(shù)據(jù)，不同傳感器之間存在一定誤差，影像分辨率偏低，且對去云后的影像采取時間插值處理，對預(yù)測結(jié)果有一定影響。（3）本研究所使用的森林資源連續(xù)清查樣地數(shù)據(jù)偏陳舊，固定樣地偏少，后續(xù)研究可考慮使用較新的數(shù)據(jù)；本研究結(jié)果與涂宏濤等［17］基于云南省第九次（2017 年）森林資源清查數(shù)據(jù)研究結(jié)果相比，受研究尺度、研究方法、地質(zhì)差異與森林覆蓋率等條件的影響，導(dǎo)致2017 年滇東地區(qū)碳密度22.74 t/hm2低于云南省碳密度44.96 t/hm2。（4）森林生態(tài)系統(tǒng)碳儲量包括地上部分與地下部分，本研究只考慮了地上部分的森林植被層，導(dǎo)致本研究昆明市碳儲量結(jié)果低于李俊等［28］利用PLUS 與InVEST 模型估算的昆明市土地利用碳儲量結(jié)果，后續(xù)研究需進一步完善。（5）梯度提升回歸樹模型在森林碳儲量估算方面的研究較少，其可行性與適用性有待進一步探討。二是樣地本身與外界干擾方面：（1）人口密度大，林地面積占比??；退耕還林面積少，造林地可納入森林覆蓋率計算的面積??；森林覆蓋率［29］。（2）中幼林面積比重較大；立地條件差，林木生長慢；天然次生林面積大，林分質(zhì)量不高［29］。（3）生態(tài)環(huán)境建設(shè)投入經(jīng)費有限，導(dǎo)致造林保存率過低，影響森林面積增加；自然災(zāi)害對林業(yè)生產(chǎn)造成一定破壞；森林管護責(zé)任落實不到位，森林面積減?。?0］。

4.2 結(jié)論

（1）GBRT 模型的R2、RMSE、rRMSE、P均優(yōu)于偏最小二乘回歸模型、KNN 模型、隨機森林模型，其可作為估算區(qū)域尺度森林碳儲量的一種可靠途徑。

（2）1997 年滇東地區(qū)每公頃碳儲量占比最大范圍為10～20 t，到2022 年每公頃碳儲量在20 t以下的區(qū)域均向20～30 t 范圍轉(zhuǎn)移，碳密度呈增長趨勢，增幅為3.9%。

（3）1997—2022 年滇東地區(qū)森林總碳儲量總體呈上升趨勢，增幅為4.2%；除玉溪市森林總碳儲量呈現(xiàn)下降趨勢外，其他各地州均呈現(xiàn)不同幅度的增長。

（4）滇東地區(qū)森林碳儲量具有空間異質(zhì)性。2012 年以后，空間自相關(guān)程度開始下降；1997—2022 年，文山州碳儲量高值聚類分布呈現(xiàn)增長，其他地州呈現(xiàn)減少，低值聚類分布在曲靖和昆明呈現(xiàn)減少，其他地州趨于平穩(wěn)。