引入環(huán)境變量的香格里拉市高山松碳儲量估測*

殷唐燕,張加龍,廖易,王飛平,曹軍,和云潤,陳朝情,肖慶琳

(1.西南林業(yè)大學(xué) 林學(xué)院,云南 昆明 650224;2.西北農(nóng)林科技大學(xué) 機(jī)械與電子工程學(xué)院,陜西 楊陵 712100)

森林作為生物圈中最大的陸地碳庫,其儲存著全球約80%的地上碳儲量和40%的地下碳儲量,在全球碳循環(huán)中發(fā)揮著不可或缺的作用[1-3]。森林地上碳儲量能夠?qū)侠淼纳仲Y源保護(hù)與管理提供依據(jù),是固碳能力的重要參考指標(biāo)[4-5]。隨著“碳達(dá)峰、碳中和”目標(biāo)的提出,如何準(zhǔn)確、有效地估算森林碳儲量對于碳監(jiān)測和降低大氣CO2濃度具有重要意義。

相較于傳統(tǒng)的森林碳儲量測算方法[6],遙感估測森林碳儲量因具有快速、成本低和大范圍研究等優(yōu)點(diǎn)被廣泛應(yīng)用于林業(yè)領(lǐng)域[7-8]。常用于森林碳儲量估測的預(yù)測變量有紋理因子、植被指數(shù)和信息增強(qiáng)因子等[9-10],此外,部分學(xué)者在估測森林碳儲量時(shí)還考慮了地形和氣候等環(huán)境因素,并發(fā)現(xiàn)環(huán)境變量能在一定程度上提高模型的估測精度。蔣云姣等[11]研究發(fā)現(xiàn)海拔、坡度、亮度指數(shù)、濕度指數(shù)等是影響森林地上部分生物量的重要環(huán)境因素。何瀟等[12]通過對比林分生物量基礎(chǔ)模型和含氣候變量的林分生物量模型的精度,發(fā)現(xiàn)氣候變量能明顯改善模型效果。選擇合適的變量篩選方法可以提高模型精度、降低運(yùn)算時(shí)間和模型的復(fù)雜性[13-14]。香格里拉市位于青藏高原橫斷山區(qū),地形結(jié)構(gòu)復(fù)雜,有著多種氣候類型,區(qū)域內(nèi)樹木的生長狀況及樹種的空間分布容易受到環(huán)境因素的影響,從而影響到森林碳儲量,因此,引入環(huán)境變量進(jìn)行高山松碳儲量估測尤為重要[15]。

高山松(Pinusdensata)在香格里拉分布廣泛,具有重要的經(jīng)濟(jì)價(jià)值和生態(tài)價(jià)值,對該地區(qū)的碳收支平衡具有重要影響。目前,已經(jīng)開展了較多遙感估測高山松生物量/碳儲量的研究[16-17],王書賢等[18]對高山松生物量和遙感因子進(jìn)行空間自相關(guān)和異質(zhì)性分析,發(fā)現(xiàn)高山松變化主要受人為干擾、立地條件影響。曹軍等[19]和廖易等[20]通過建立不同變量組合的高山松估測模型,發(fā)現(xiàn)引入地形可以顯著提高高山松地上碳儲量的估測精度。關(guān)于引入氣候、土壤等環(huán)境變量構(gòu)建高山松遙感估測的相關(guān)研究較少?；诖?研究以Landsat TM/OLI影像為遙感數(shù)據(jù)源,結(jié)合香格里拉市的森林資源連續(xù)清查數(shù)據(jù)、氣象數(shù)據(jù)等,采用Pearson相關(guān)性法、Kendall’s τ相關(guān)性法、Spearman相關(guān)性法、距離相關(guān)性法和決策樹法對預(yù)測變量進(jìn)行篩選,構(gòu)建隨機(jī)森林模型,并在模型中引入環(huán)境變量(氣溫、降水、年日照時(shí)長、年均溫、人口密度、地表溫度和土壤濕度),以探究不同變量篩選方法及引入環(huán)境變量是否能提高地上碳儲量建模精度,為高山松碳儲量建模精度的提高提供新思路。

1 研究區(qū)與方法

1.1 研究區(qū)概況

香格里拉市位于云南省迪慶藏族自治州東部、青藏高原的東南端和橫斷山脈的東側(cè),地處26°52′～28°52′N,99°20′～100°19′E,平均海拔為3 459 m。受地形和西南季風(fēng)影響,該區(qū)干濕季節(jié)分明,年均溫為5.4℃,年均降水量為268～945 mm,無霜期為129～197 d。自然資源豐富,森林覆蓋率達(dá)75%,區(qū)域內(nèi)優(yōu)勢樹種有高山松、冷杉(Abiesfabri)、云杉(Piceaasperata)、云南松(Pinusyunnanensis)等[21],見圖1。

圖1 研究區(qū)概況圖

1.2 數(shù)據(jù)來源與處理

1.2.1 森林資源連續(xù)清查數(shù)據(jù)

樣地?cái)?shù)據(jù)來源于1987、1992、1997、2002、2007、2012、2017年香格里拉市的森林資源連續(xù)清查數(shù)據(jù),共136個(gè)高山松樣地?cái)?shù)據(jù),樣地按照6 km×8 km的網(wǎng)格分布,面積為(28.28 m×28.28 m)0.08 hm2。

首先利用樣地平均樹高和胸徑計(jì)算平均生物量(t),再根據(jù)樣地平均生物量和樣地高山松株數(shù)計(jì)算樣地總生物量,高山松單木異速生長方程[22]見公式①。根據(jù)推算得到的每公頃生物量(t/hm2)乘以含碳率得到樣地碳儲量[23],具體見公式②。高山松的含碳率將根據(jù)《森林生態(tài)系統(tǒng)碳儲量計(jì)量指南》中推薦的高山松含碳率(0.501)確定。

W=0.073×D1.739×H0.880

①

C=W×CF

②

式中:W代表高山松地上生物量,D代表胸徑,H代表樹高,C代表碳儲量,CF代表含碳率。

在得到高山松碳儲量數(shù)據(jù)后,通過剔除碳儲量過小以及離群值后,剩余124個(gè)樣地用于構(gòu)建高山松地上碳儲量模型,詳情見表1。

表1 樣地碳儲量統(tǒng)計(jì)結(jié)果

表2 預(yù)測變量

1.2.2 Landsat數(shù)據(jù)

從http://glovis.usgs.gov/收集了1987、1992、1997、2002、2007、2012、2017年的Landsat TM/OLI遙感影像,條帶號分別為131/041、132/040、132/041,共7期21景,影像分別經(jīng)過輻射定標(biāo)、大氣校正、鑲嵌、裁剪等預(yù)處理。

1.2.3 環(huán)境變量數(shù)據(jù)

研究所選擇的環(huán)境變量包括氣象數(shù)據(jù)(年均降水、氣溫、干燥度、潛在蒸散發(fā))、人口密度、地表溫度和土壤濕度。

氣象數(shù)據(jù)來源于國家青藏高原科學(xué)數(shù)據(jù)中心(http://data.tpdc.ac.cn)。人口密度數(shù)據(jù)來源于官方網(wǎng)站(https://www.WorldPop.org/),數(shù)據(jù)集只有2002、2007、2012和2017年的人口密度數(shù)據(jù),1987、1992和1997年的人口密度數(shù)據(jù)將通過《云南省統(tǒng)計(jì)年鑒》和《中甸縣志》計(jì)算得到人口密度。為統(tǒng)一數(shù)據(jù)格式,將氣象數(shù)據(jù)通過重采樣至30 m分辨率。地表溫度采用單窗算法[24]反演得到,土壤濕度使用簡化后的溫度植被干旱指數(shù)(TVDI),利用地表溫度和歸一化植被指數(shù)(NDVI)得到[25],與土壤水分狀況直接相關(guān)[26]。

1.2.4 地形數(shù)據(jù)

從地理空間數(shù)據(jù)云平臺(http://www.gscloud.cn/)下載分辨率為30 m的DEM數(shù)據(jù),處理后提取海拔、坡度和坡向數(shù)據(jù)作為地形變量加入到預(yù)測變量進(jìn)行篩選。

1.3 研究方法

研究基于連清樣地?cái)?shù)據(jù)、Landsat TM/OLI影像數(shù)據(jù)和地形數(shù)據(jù)進(jìn)行相關(guān)處理后提取預(yù)測變量,從Landsat TM/OLI影像數(shù)據(jù)和環(huán)境數(shù)據(jù)提取七種不同環(huán)境變量,使用五種變量篩選方法對預(yù)測變量進(jìn)行篩選,確定最優(yōu)的預(yù)測變量組合直接引入環(huán)境變量建立隨機(jī)森林模型并評價(jià)精度,選擇引入環(huán)境變量后最優(yōu)的模型進(jìn)行高山松碳儲量反演制圖。技術(shù)路線見圖2。

圖2 技術(shù)路線圖

1.3.1 預(yù)測變量選擇

研究使用的預(yù)測變量包括原始單波段、波段組合、植被指數(shù)、紋理因子、信息增強(qiáng)因子和地形變量,其中紋理因子提取了影像第1、2、3、4、5、7波段的3×3、5×5、7×7、9×9、11×11、13×13、15×15、17×17、19×19窗口的9種紋理。詳見圖2。采用Pearson、Spearman、Kendall’s τ、距離相關(guān)性法以及決策樹法分別選取與碳儲量顯著相關(guān)的前10的預(yù)測變量。使用ArcGIS軟件以樣地中心點(diǎn)為基礎(chǔ),提取該點(diǎn)與Landsat影像對應(yīng)位置的像元值。

最佳變量組合是提高建模精度的關(guān)鍵。Pearson相關(guān)性法能夠描述變量X與變量Y之間的線性相關(guān)性,是特征變量篩選最常用的方法,適用于正態(tài)分布的連續(xù)變量。Spearman相關(guān)性法[27]和Kendall’s τ相關(guān)性法[28]都是利用非參數(shù)的思想描述兩個(gè)有序變量之間的相關(guān)性強(qiáng)弱,對非線性關(guān)系更為敏感。距離相關(guān)性法[29]是描述任意變量間獨(dú)立性的度量關(guān)系數(shù)越大,越能說明預(yù)測變量越重要[30]。決策樹法能利用平方誤差對分裂變量進(jìn)行選擇[31],并生成二叉樹,將目標(biāo)變量區(qū)分出來,具有運(yùn)算速度快、結(jié)果易理解和解釋等優(yōu)點(diǎn)。

1.3.2 模型構(gòu)建

研究采用隨機(jī)森林(random forest,RF)方法進(jìn)行高山松碳儲量模型的構(gòu)建,將124個(gè)樣本中的80%(99個(gè))用于模型的擬合,20%(25個(gè)) 用于模型的檢驗(yàn)。

隨機(jī)森林是改進(jìn)后的決策樹算法,適用于多數(shù)分類與回歸的問題,是由一系列的決策樹組成[32]。通過Python語言的Sklearm工具包中的“Random forest”算法進(jìn)行建模分析,在構(gòu)建過程中需要對運(yùn)行參數(shù)進(jìn)行調(diào)整,最終確定參數(shù)范圍:最大迭代次數(shù)為50～300,最大深度為15～20,葉節(jié)點(diǎn)最少樣本數(shù)以及內(nèi)部節(jié)點(diǎn)再劃分所需最小樣本數(shù)采用默認(rèn)值,分別為1和2。

1.3.3 模型精度評價(jià)

使用決定系數(shù)(R2)、均方根誤差(RMSE)、相對均方根誤差(rRMSE)、預(yù)測精度(P)作為模型精度評價(jià)指標(biāo)。

③

④

⑤

⑥

2 結(jié)果與分析

2.1 預(yù)測變量篩選結(jié)果

通過對預(yù)測變量等進(jìn)行篩選,結(jié)果見表3。各方法選擇的變量組合各不相同,Pearson和距離相關(guān)性法所選擇的變量都為紋理因子,并且有7個(gè)相同的變量,相關(guān)性最高的為15×15窗口第7波段的角二階矩(0.285**);Spearman和Kendall’s τ相關(guān)性法所篩選的變量相同,其變量包含9個(gè)紋理因子和1個(gè)植被指數(shù),相關(guān)性最高的變量為ND54,其相關(guān)性分別為-0.303**和-0.203**;決策樹法篩選的變量則包含了多種類型的變量,包含7個(gè)紋理因子、1個(gè)波段組合因子、1個(gè)植被指數(shù)和1個(gè)信息增強(qiáng)因子,重要性最高的為19×19窗口第5波段的偏度(13.04%)。

表3 不同變量選擇方法及其相關(guān)性(重要性)

表4 不同變量組合模型建模效果

通過分析發(fā)現(xiàn)紋理因子相較于其他因子是較為重要的因子,通過相關(guān)性發(fā)現(xiàn),第1、5和7波段的信息隨高山松碳儲量具有重要影響,紋理信息中的偏度(SK)、角二階矩(SM)、均一性(HO)和相異性(DI)的占比較多。

2.2 高山松碳儲量估測模型構(gòu)建

2.2.1 基于不同變量組合的建模結(jié)果

經(jīng)過5種篩選方法得到的變量組合分別進(jìn)行RF建模。不同變量組合的RF建模結(jié)果如圖4所示,基于決策樹法選擇的變量組合建立的RF模型精度(R2=0.845,RMSE=10.076 t/hm2,rRMSE=29.251%,P=0.747)明顯優(yōu)于其他變量組合,相較于Pearson相關(guān)性法,R2提高了4.20%,RMSE降低了11.31%,rRMSE降低了1.92%,P提高了0.02%,其預(yù)測得到的高山松碳儲量與連清樣地的計(jì)算得到的碳儲量更為接近(圖3)。基于決策樹法的選擇的變量組合建模精度最高,其余變量組合建模精度依次為Spearman、距離相關(guān)、Kendall’s τ、Pearson相關(guān)性法。

圖4 不同海拔范圍內(nèi)地表溫度對碳儲量的影響

由此可知不同變量篩選方法能夠影響模型精度,使用決策樹法篩選的變量組合的建模精度明顯優(yōu)于其他方法。因此,采用決策樹法選擇的變量組合進(jìn)行引入環(huán)境變量的建模。

2.2.2 引入不同環(huán)境變量的建模結(jié)果

引入環(huán)境變量后模型效果見表5。從建模結(jié)果看,引入地表溫度的模型R2提升效果最好,為0.893,RMSE為8.371 t/hm2,rRMSE為22.863%,P為0.812。其余依次為降水量、土壤濕度、潛在蒸散發(fā)和干燥度,引入人口密度的模型提升效果最低,R2為0.863,RMSE為9.516 t/hm2,rRMSE為27.396%,P為0.761。相較于決策樹法篩選的變量組合模型效果,引入地表溫度后模型R2提高了4.80%,RMSE降低了1.71 t/hm2,rRMSE降低了6.388 %,P提高了6.5%。

表5 引入不同環(huán)境變量建模效果

將海拔分為1 500～3 000 m、3 000～3 200 m、3 200～3 400 m、3 400～3 700 m、3 700～3 900 m共5個(gè)梯度探究范圍內(nèi)高山松碳儲量和地表溫度的關(guān)系(圖4)。在1 500～3 000 m、3 000～3 200 m、3 400～3 700 m和3 700～3 900 m海拔范圍內(nèi),地表溫度和高山松碳儲量呈明顯正相關(guān)關(guān)系,地表溫度增加,高山松碳儲量也隨之增加,在3 200～3 400 m內(nèi)地表溫度對高山松碳儲量的影響不明顯,可以忽略地表溫度對高山松碳儲量的影響。總體而言,地表溫度對高山松碳儲量的影響呈正相關(guān)。

引入環(huán)境變量后,各模型效果都有不同程度的提升,能夠在一定程度上減少模型的誤差。因此選擇基于決策樹法篩選的變量組合并引入地表溫度進(jìn)行香格里拉市高山松碳儲量反演。

2.3 引入地表溫度的高山松碳儲量反演結(jié)果

通過引入地表溫度的RF模型結(jié)合已有的香格里拉市高山松分布范圍對1987—2017年共7年的高山松碳儲量進(jìn)行反演制圖(圖5)。高山松在香格里拉市的分布范圍較廣,碳儲量變化呈現(xiàn)“增加—減少—增加”趨勢,1987—2002年碳儲量逐漸增加,2002—2012年碳儲量出現(xiàn)下降,2012—2017年碳儲量出現(xiàn)大幅增加。最終得到1987、1992、1997、2002、2007、2012、2017年香格里拉市高山松碳儲量估測值分別為3 375.823×104t、3 500.993×104t、3 563.118×104t、3 610.929×104t、3 555.959×104t、3 628.858×104t和4 027.089×104t,與岳彩榮[33]基于遙感估測香格里拉市高山松碳儲量估測結(jié)果(4 895.30×104t)具有可對比性。30年間森林碳儲量是處于上升趨勢,共增加了651.266×104t。將1987和2017年的碳儲量變化數(shù)值分為四個(gè)等級,分別為顯著減少(<-20 t/hm2)、少量減少(-20～0 t/hm2)、少量增加(0～20 t/hm2)和顯著增加(>20 t/hm2),從下圖可以看出碳儲量變化整體呈現(xiàn)少量增加。

圖5 各年份高山松碳儲量分布

3 討論與結(jié)論

3.1 討論

(1)決策樹法篩選變量組合的可行性

變量篩選的結(jié)果好壞會直接影響模型的質(zhì)量,選擇合適的建模變量是準(zhǔn)確估測碳儲量的關(guān)鍵步驟。本研究基于5種變量篩選方法選擇預(yù)測變量,從不同的變量篩選方法來看,發(fā)現(xiàn)決策樹法篩選出的變量組合模型精度優(yōu)于其他篩選方法,這與Jiang等[34]和王迪等[35]的研究結(jié)果一致,機(jī)器學(xué)習(xí)算法在變量篩選中具有可解釋強(qiáng)、使用方便和準(zhǔn)確率高的優(yōu)點(diǎn)。在實(shí)際運(yùn)用中,沒有最好的方法,只有最合適的方法。決策樹法篩選的變量結(jié)果包含7個(gè)紋理因子,1個(gè)植被指數(shù)因子,1個(gè)信息增強(qiáng)因子和1個(gè)波段組合因子,紋理因子中有6個(gè)為紋理信息中的偏度(SK),紋理因子普遍表現(xiàn)出與高山松碳儲量較高的相關(guān)性,這與相關(guān)學(xué)者的研究[36-37]一致,紋理因子相較于其他變量能夠較為準(zhǔn)確的描述高山松的形態(tài)特征,對提高高山松的建模精度具有積極影響。

(2)光學(xué)遙感在估算森林碳儲量方面的局限性

研究利用RF模型進(jìn)行高山松碳儲量估測,結(jié)果表明,RF模型預(yù)測碳儲量與樣地實(shí)測碳儲量結(jié)果較為接近,已有大量學(xué)者研究[38-39]證明,RF模型在高山松的資源估測方面具有很強(qiáng)的適用性。但與相關(guān)學(xué)者的研究[40-41]出現(xiàn)相同問題,在碳儲量預(yù)測中存在低值高估和高值低估的現(xiàn)象,主要是由于Landsat數(shù)據(jù)為光學(xué)數(shù)據(jù),空間分辨率較大,相較于高分辨率數(shù)據(jù)在區(qū)域尺度上的建模精度較低,羅洪斌等[42]以LiDAR和Landsat8 OLI對橡膠林地上生物量進(jìn)行估測,發(fā)現(xiàn)聯(lián)合LiDAR和Landsat8 OLI的模型估測精度較高,且Landsat8 OLI數(shù)據(jù)的模型精度低于單一的LiDAR數(shù)據(jù)。

(3)建模過程中引入環(huán)境變量的必要性

在本研究中,引入不同的環(huán)境變量后建模效果都得到了提升,特別是在引入地表溫度后,建模效果提升最為顯著。這主要是由于地表溫度影響了土壤中的有機(jī)質(zhì)分解速率,從而影響了樹木的生長,已有研究發(fā)現(xiàn),隨著海拔的增加,溫度逐漸降低,有機(jī)質(zhì)的分解速率也隨之降低[43-44],地表溫度、濕度等環(huán)境變量與生物量之間存在密切關(guān)系[45],香格里拉地區(qū)位于低緯度、高海拔地帶,地勢復(fù)雜,山地縱橫,因此地表溫度對高山松的碳儲量的影響更為明顯。之前的研究也表明,氣候變量和植被指數(shù)等因子相較于反射率等因子對于森林生物量的估測更為重要[15],在構(gòu)建高山松生物量模型的過程中引入氣候變量,高山松生物量模型的R2和P均有所提高[46]。在本研究中,僅選擇了部分環(huán)境變量進(jìn)行分析與高山松碳儲量之間的關(guān)系,在未來的研究可以進(jìn)一步探究更多環(huán)境變量對高山松碳儲量的影響。

3.2 結(jié)論

研究從Landsat TM/OLI遙感影像和地形數(shù)據(jù)中提取預(yù)測變量,結(jié)合樣地?cái)?shù)據(jù),建立RF模型估測香格里拉市高山松碳儲量,并引入環(huán)境變量,對比模型精度,通過最優(yōu)模型反演高山松碳儲量。結(jié)果表明:(1)在篩選出的變量中,紋理因子與高山松碳儲量顯著相關(guān),其相關(guān)性最高;(2)不同變量篩選方法所選擇的變量組合對建模效果有較大影響,其中通過決策樹法篩選的變量組合模型表現(xiàn)最好,預(yù)測得到的高山松碳儲量與連清樣地計(jì)算得到的碳儲量更為接近,相較于精度最低的Pearson相關(guān)性法,R2提高了4.20%,RMSE降低了11.31%,rRMSE降低了1.92%,P提高了0.02%;(3)引入地表溫度的RF模型效果明顯優(yōu)于引入其他環(huán)境變量的RF模型,其R2為0.893,相較于基于決策樹法選擇的變量組合的模型擬合精度提高了4.80%;(4)高山松碳儲量在30年間的變化十分明顯,共增加了651.266×104t,其變化呈現(xiàn)為“增加—減少—增加”趨勢。