基于兩種空間估算模型的喬木林地上碳密度估算

王海賓 鄭冬梅 王少杰 賈筱昕 許等平

(1.國(guó)家林業(yè)和草原局林產(chǎn)工業(yè)規(guī)劃設(shè)計(jì)院， 北京 100010; 2.國(guó)家林業(yè)和草原局調(diào)查規(guī)劃設(shè)計(jì)院， 北京 100714)

0 引言

森林參數(shù)的準(zhǔn)確估算是林學(xué)研究的熱點(diǎn)之一[1-2]，對(duì)森林經(jīng)營(yíng)管理、生物多樣性保護(hù)、全球氣候變化研究具有重要意義。掌握森林參數(shù)現(xiàn)狀及分布情況，有利于森林經(jīng)營(yíng)管理者科學(xué)管理森林，對(duì)在不同層面制定森林經(jīng)營(yíng)管理政策、調(diào)整經(jīng)營(yíng)管理策略具有重要作用[3]。

近年來，遙感技術(shù)得到了快速發(fā)展，因即時(shí)、快速和低成本等特點(diǎn)而使其在國(guó)土、農(nóng)林、水利、礦業(yè)等領(lǐng)域得到了廣泛應(yīng)用，取得了較好的應(yīng)用效果[4-5]。尤其是在區(qū)域和大尺度范圍的監(jiān)測(cè)上，光學(xué)遙感數(shù)據(jù)覆蓋范圍廣、獲取周期短、數(shù)據(jù)成本低廉等，具有激光雷達(dá)、微波雷達(dá)等其他遙感數(shù)據(jù)無法比擬的優(yōu)勢(shì)[6]。因此，應(yīng)用光學(xué)遙感技術(shù)結(jié)合地面樣地?cái)?shù)據(jù)對(duì)森林參數(shù)進(jìn)行估算成為必要的手段[7]。

筆者曾參與全國(guó)森林植被碳儲(chǔ)量遙感估算工作[8-9]，其中的喬木林碳儲(chǔ)量估算主要基于ZY-3多光譜影像，結(jié)合地面調(diào)查數(shù)據(jù)，應(yīng)用傳統(tǒng)的線性模型，以省為單位來估算喬木林碳儲(chǔ)量。由于覆蓋全省的ZY-3多光譜影像數(shù)據(jù)量較大，存在時(shí)相不一致的問題，同時(shí)因拍攝角度、云霧等多因素干擾，增加了估算結(jié)果的不確定性。為改善估算效果，依據(jù)估算目標(biāo)范圍以及遙感影像空間分辨率的適用性等因素，選擇Landsat 8多光譜影像，結(jié)合實(shí)地調(diào)查數(shù)據(jù)進(jìn)行了森林植被碳儲(chǔ)量估算的研究分析，以期為區(qū)域及全國(guó)尺度范圍的森林碳儲(chǔ)量估算提供更好的數(shù)據(jù)源及方法，減小估算結(jié)果的不確定性，達(dá)到增強(qiáng)監(jiān)測(cè)成果時(shí)效性、提高監(jiān)測(cè)效率、降低工作成本、改進(jìn)和優(yōu)化全國(guó)森林植被碳儲(chǔ)量估算數(shù)據(jù)源和估算方法的目的。

為挖掘Landsat 8多光譜影像在全國(guó)森林植被碳儲(chǔ)量遙感估算中的應(yīng)用潛力，同時(shí)探索更好的估算方法，本研究以位于浙江省內(nèi)的一景Landsat 8影像覆蓋范圍為研究區(qū)，以喬木林地上碳密度為研究對(duì)象，基于研究區(qū)內(nèi)的地面樣地調(diào)查數(shù)據(jù)，通過計(jì)算獲得喬木林地上碳密度數(shù)據(jù)并將其作為因變量，采用協(xié)同克里格插值法和地理加權(quán)回歸方法，構(gòu)建喬木林地上碳密度模型，并對(duì)建模結(jié)果進(jìn)行對(duì)比分析，以期為區(qū)域及全國(guó)范圍的喬木林地上碳密度及森林參數(shù)遙感快速估算提供更好的數(shù)據(jù)源和估算方法，并提高遙感估算結(jié)果的精度和可靠性。

1 研究區(qū)及數(shù)據(jù)獲取

1.1 研究區(qū)概況

研究區(qū)為浙江省內(nèi)一景Landsat 8影像覆蓋的范圍(圖1)，包括金華市、杭州市、紹興市、衢州市、麗水市和溫州市，地理坐標(biāo)為東經(jīng)118°28′20″～120°49′11″，北緯27°48′3″～29°55′26″。區(qū)內(nèi)以丘陵和山地為主，屬于亞熱帶季風(fēng)氣候，四季分明，雨熱同季。區(qū)內(nèi)土壤類型多樣，主要有紅壤、黃壤、粗骨土、水稻土、潮土、鹽土等類型。研究區(qū)森林資源豐富，樹種較多，森林類型以針葉林為主，闊葉林、針闊混交林較多，主要包括馬尾松(PinusmassonianaLamb.)、杉木(Cunninghamialanceolata(Lamb.)Hook.)、櫟類(CryptomeriafortuneiHooibrenk ex Otto et Dietr.)、木荷(SchimasuperbaGardn. et Champ.)以及其他軟闊類等優(yōu)勢(shì)樹種。

圖1 研究區(qū)內(nèi)調(diào)查樣地分布Fig.1 Distribution of survey samples in study area

1.2 數(shù)據(jù)獲取

1.2.1地面數(shù)據(jù)

地面數(shù)據(jù)獲取時(shí)間為2012年。為使樣地具有代表性，依據(jù)優(yōu)勢(shì)樹種和林分年齡(幼齡林、中齡林、近熟林、成熟林、過熟林)，設(shè)置和調(diào)查了96個(gè)樣地(圖1)，樣地形狀為正方形，面積為0.067 hm2，涵蓋針葉林、闊葉林、混交林等不同森林類型。樣地屬性因子包括樹種、胸徑(起測(cè)胸徑為5 cm)、平均樹高、郁閉度、株數(shù)、年齡、坡向、坡度、海拔、喬木林蓄積量等因子。

將獲取的樣地喬木林地上蓄積量乘以生物量換算因子(Biomass expansion factor，BEF)得到喬木林地上生物量，并轉(zhuǎn)換為每公頃生物量(t/hm2)。采用沈楚楚[10]提出的浙江省四大樹種的BEF，得到每公頃生物量分別為杉木林0.745 3 t/m3、馬尾松林0.883 9 t/m3、軟闊類1.657 2 t/m3、硬闊類1.070 5 t/m3。地面樣地喬木林地上生物量數(shù)據(jù)統(tǒng)計(jì)結(jié)果見表1。

表1 地面樣地喬木林地上生物量統(tǒng)計(jì)結(jié)果Tab.1 Statistic results of aboveground biomass of arbor forest for ground sample t/hm2

將每公頃生物量乘以樹種含碳率即可得到單位面積喬木林地上碳密度。樹種含碳率參考《土地利用、土地利用變化與林業(yè)碳匯計(jì)量監(jiān)測(cè)技術(shù)指南》，指南中的樹種含碳率可覆蓋研究區(qū)內(nèi)的各個(gè)優(yōu)勢(shì)樹種。

在獲得喬木林地上碳密度數(shù)據(jù)后，采用二倍標(biāo)準(zhǔn)差法剔除異常值，結(jié)果剩余95個(gè)地面樣地?cái)?shù)據(jù)，將其作為因變量參與喬木林地上碳密度模型的構(gòu)建。

1.2.2遙感影像及輔助數(shù)據(jù)

本研究選用一景Landsat 8 影像作為遙感數(shù)據(jù)源，覆蓋整個(gè)研究區(qū)，影像軌道號(hào)為119/040，空間分辨率為30 m，采集時(shí)間為2014年10月26日。由于樹木生長(zhǎng)較慢，本研究假設(shè)地面數(shù)據(jù)采集時(shí)間與影像獲取時(shí)間不存在差異。應(yīng)用ENVI 5.3軟件對(duì)Landsat 8多光譜影像進(jìn)行預(yù)處理，包括輻射定標(biāo)、大氣校正。本研究所采用波段序號(hào)為2～7，分別為藍(lán)、綠、紅、近紅外、短波紅外1、短波紅外2共6個(gè)波段。

本研究采用的輔助數(shù)據(jù)為覆蓋研究區(qū)1∶50 000數(shù)字高程模型(Digital elevation model，DEM)數(shù)據(jù)。

2 研究方法

2.1 空間自相關(guān)分析

采用協(xié)同克里格插值法對(duì)喬木林地上碳密度建模前，需要對(duì)使用的碳密度數(shù)據(jù)進(jìn)行描述性統(tǒng)計(jì)分析，采用柯爾莫哥洛夫-斯蜜諾夫檢驗(yàn)(K-S檢驗(yàn))[9]對(duì)數(shù)據(jù)進(jìn)行分析，檢驗(yàn)其是否符合正態(tài)分布。

此外，采用地理加權(quán)回歸法建模前，需要對(duì)喬木林地上碳密度進(jìn)行空間自相關(guān)分析，檢驗(yàn)研究數(shù)據(jù)是否隨機(jī)分布。采用莫蘭指數(shù)(Moran’sI)進(jìn)行檢驗(yàn)[11]，計(jì)算公式如下

(1)

式中n——樣本總數(shù)

I——莫蘭指數(shù)

xi、xj——空間位置i和j的變量取值

wij——空間位置i和j的權(quán)重

Moran’sI的取值范圍在[-1，1]之間，小于0表示負(fù)相關(guān)，等于0表示不相關(guān)，大于0則表示正相關(guān)。

2.2 自變量提取及優(yōu)選

2.2.1自變量提取

自變量一般選擇光學(xué)遙感數(shù)據(jù)、地形數(shù)據(jù)以及經(jīng)計(jì)算衍生的一些因子，常用的包括植被指數(shù)[12]、紋理特征[13]、主成分變換因子[3]、纓帽變換因子[3]、地形因子[14]等。本研究基于Landsat 8多光譜波段及DEM數(shù)據(jù)，計(jì)算了植被指數(shù)(9個(gè))、紋理特征(48個(gè))、主成分變換因子(3個(gè))、纓帽變換因子(3個(gè))和地形因子(3個(gè))共66個(gè)特征作為自變量參與模型構(gòu)建。其中，選用的植被指數(shù)見表2，其中BLUE、RED、NIR、SWIR1、SWIR2分別為藍(lán)、紅、近紅、短波紅外1、短波紅外2波長(zhǎng)； SAVI中，L為0.5；ARVI中，r為校正系數(shù)，取1.0；EVI中，增益系數(shù)G、土壤調(diào)整因子L、修正系數(shù)C1和C2分別為2.5、0.10、6.0和7.5。紋理特征采用灰度共生矩陣方法中的8個(gè)紋理特征(平均值、方差、均一性、對(duì)比度、相異性、熵、角二階矩、相關(guān)性)；主成分變換和纓帽變換因子分別選擇主成分分析獲取的前3個(gè)波段和纓帽變換的前3個(gè)波段(亮度、綠度和濕度)；地形因子采用DEM數(shù)據(jù)提取的海拔、坡度、坡向3個(gè)因子。按照地面實(shí)測(cè)樣地位置，提取各個(gè)樣地對(duì)應(yīng)的自變量因子。

表2 植被指數(shù)Tab.2 Vegetation indices

2.2.2自變量?jī)?yōu)選

在構(gòu)建模型時(shí)，自變量對(duì)所構(gòu)建模型的穩(wěn)定性有很大影響，有些因子對(duì)碳密度估算的影響很小，甚至沒有影響，則將其從自變量因子集中剔除，將剩余的因子組成最優(yōu)自變量集。在本研究中，分兩步對(duì)自變量進(jìn)行優(yōu)選：首先利用相關(guān)性分析提取自變量因子；其次選用平均殘差平方和(Residual mean square，RMSQ)準(zhǔn)則[15]，對(duì)自變量做進(jìn)一步優(yōu)選。

2.3 建模方法

2.3.1協(xié)同克里格插值法

協(xié)同克里格插值法(Co-Kriging，COK)在變量估測(cè)上具有較好的預(yù)測(cè)精度[8,16-18]，它是利用兩個(gè)或兩個(gè)以上的變量，將其中一個(gè)作為主變量，剩余變量作為輔助變量，將主變量的自相關(guān)性和主輔變量的交互相關(guān)性結(jié)合起來進(jìn)行無偏最優(yōu)估計(jì)[17-18]，協(xié)同克里格插值公式如下

(2)

其中

式中Z*(x0)——待估點(diǎn)處的蓄積量估測(cè)值

λ1i、λ2i——主變量Z1和輔變量Z2實(shí)測(cè)值的權(quán)重

n、m——參與估測(cè)x0點(diǎn)的蓄積量和輔助變量Z2的實(shí)測(cè)值數(shù)目

2.3.2地理加權(quán)回歸法

地理加權(quán)回歸法(Gergraphic weighted regression，GWR)是用回歸原理研究具有空間(或區(qū)域)分布特征的兩個(gè)或多個(gè)變量之間數(shù)量關(guān)系的方法，它是基于普通全局回歸模型的擴(kuò)展，將空間影響考慮到模型構(gòu)建中，以距離權(quán)重的形式加入到模型之中，故地理加權(quán)回歸模型[3]的基本形式如下

=0(x,y)+1(x,y)X1+2(x,y)X2+…+N(x,y)Xn+ε

(3)

式中x、y——空間樣本的坐標(biāo)值

N——樣本數(shù)0——截距

ε——模型殘差，服從N(0，σ2)分布

Xn——因子變量——因變量

采用地理加權(quán)回歸法構(gòu)建模型時(shí)，模型帶寬的選擇至關(guān)重要，在很大程度上決定模型的預(yù)測(cè)準(zhǔn)確度[3,14]。本研究采用最小化信息準(zhǔn)則(Akaike information criterion，AIC)[9]來確定最佳帶寬：當(dāng)AIC值取最小值時(shí)，即可獲得最優(yōu)帶寬及地理加權(quán)回歸模型。AIC值計(jì)算公式為

AIC=-2lnH+2k

(4)

式中k——未知參數(shù)的數(shù)量

H——似然函數(shù)AIC——AIC值

最優(yōu)模型帶寬在AIC值最小時(shí)獲得。

此外，在完成GWR模型構(gòu)建后，需進(jìn)一步對(duì)GWR模型進(jìn)行參數(shù)的非平穩(wěn)性檢驗(yàn)，檢驗(yàn)?zāi)Ｐ托Ч?。地理加?quán)回歸理論和GWR模型參數(shù)的非平穩(wěn)性檢驗(yàn)詳見文獻(xiàn)[3，14]。

2.4 模型精度評(píng)價(jià)

本文采用留一交叉驗(yàn)證法對(duì)建模結(jié)果進(jìn)行評(píng)價(jià)。采用的評(píng)價(jià)指標(biāo)包括：決定系數(shù)(R2)、均方根誤差(Root mean square error，RMSE)、平均絕對(duì)誤差(Mean absolute error，MAE)和總預(yù)報(bào)偏差的相對(duì)誤差(Relative error，RE)。

3 研究結(jié)果

3.1 描述性統(tǒng)計(jì)及空間自相關(guān)性分析

3.1.1描述性統(tǒng)計(jì)分析

使用SPSS 20.0統(tǒng)計(jì)軟件，采用K-S檢驗(yàn)對(duì)數(shù)據(jù)進(jìn)行檢驗(yàn)。經(jīng)分析，喬木林地上碳密度數(shù)據(jù)服從正態(tài)分布，K-S檢驗(yàn)值為0.374，其平均值為19.61 t/hm2，標(biāo)準(zhǔn)差13.29 t/hm2，中位數(shù)17.91 t/hm2，最小值為0.88 t/hm2，最大值為56.02 t/hm2。

3.1.2空間自相關(guān)分析

本研究對(duì)喬木林地上碳密度進(jìn)行空間自相關(guān)分析，結(jié)果顯示研究區(qū)內(nèi)喬木林地上碳密度的莫蘭指數(shù)為0.083 7，呈空間正相關(guān)，并且在10%的顯著性水平(p=0.098)下Z為1.654(大于1.65)，高度顯著。依據(jù)莫蘭指數(shù)[11]的定義，當(dāng)其取值范圍在(0,1]之間時(shí)，喬木林地上碳密度在空間分布上呈現(xiàn)聚集狀態(tài)，喬木林地上碳密度具有整體空間分布特性，說明研究區(qū)內(nèi)的喬木林地上碳密度數(shù)據(jù)分布是非隨機(jī)的。

3.2 自變量?jī)?yōu)選結(jié)果

應(yīng)用SPSS 20.0軟件對(duì)提取的喬木林地上碳密度和66個(gè)自變量因子進(jìn)行皮爾森(Pearson)相關(guān)性分析，篩選出與碳密度極顯著相關(guān)(p<0.01)的自變量因子，提取的自變量因子共計(jì)22個(gè)，結(jié)果詳見表3。B1、B2、B3、B6分別表示Landsat 8影像的第2(藍(lán))、3(綠)、4(紅)、7(短波紅外2)波段；Mean表示平均值紋理，Entropy表示熵紋理，Second表示角二階矩紋理，Correlation表示相關(guān)性紋理，Variance表示方差紋理，Contrast表示對(duì)比度紋理，Diss表示相異性紋理。

表3 自變量因子與喬木林地上碳密度的Pearson相關(guān)系數(shù)Tab.3 Pearson correlation coefficients between independent variables and above-ground carbon density of arbor forest

注：** 表示在 0.01 水平(雙側(cè))上顯著相關(guān)，下同。

在相關(guān)性分析的基礎(chǔ)上，本研究采用RMSQ法對(duì)22個(gè)自變量因子進(jìn)一步優(yōu)選，結(jié)果見表4。經(jīng)篩選后，RVI、RVI54、RVI64、Greenness、B3_Correlation、B6_Mean共計(jì)6個(gè)因子被優(yōu)選出來。

表4 基于RMSQ篩選的自變量因子Tab.4 Independent factors filter based on RMSQ

3.3 模型構(gòu)建

3.3.1協(xié)同克里格插值法

(1)模型擬合

基于優(yōu)選的自變量，提取出與喬木林地上碳密度相關(guān)性最高的3個(gè)自變量(RVI、RVI54、RVI64)，分別建立3個(gè)自變量組，見表5。

采用協(xié)同克里格插值法對(duì)喬木林地上碳密度進(jìn)行建模。本文對(duì)ArcGIS 10.1地統(tǒng)計(jì)分析模塊中的11種變異模型進(jìn)行了計(jì)算和分析，根據(jù)變異函數(shù)的理論和評(píng)價(jià)方法，得到最優(yōu)變異函數(shù)的擬合結(jié)果(表5)。由表5可知，在協(xié)同克里格插值中，采用基于RMSQ提取的3組自變量參與插值后，塊金值C0

為0，塊金值與基臺(tái)值的比值小于25%，表明系統(tǒng)具有很強(qiáng)的相關(guān)性。依據(jù)變異函數(shù)的評(píng)價(jià)標(biāo)準(zhǔn)可知，在使用RVI進(jìn)行協(xié)同克里格插值時(shí)，孔穴效應(yīng)模型可獲得最優(yōu)的變異函數(shù)擬合效果，在使用RVI、RVI54插值時(shí)，K-貝塞耳模型可獲得最優(yōu)的擬合結(jié)果，在使用RVI、RVI54、RVI64時(shí)，K-貝塞耳模型可獲得最優(yōu)的擬合結(jié)果。

(2)精度評(píng)價(jià)

基于3組自變量，應(yīng)用協(xié)同克里格插值法對(duì)喬木林地上碳密度進(jìn)行建模，模型精度評(píng)價(jià)指標(biāo)中，RVI、RVI54、RVI64組成的自變量組所構(gòu)建模型具有最高預(yù)測(cè)精度(R2為0.47，RMSE為9.72 t/hm2，MAE為7.41 t/hm2，RE為0.12%)，優(yōu)于RVI構(gòu)建的模型精度(R2為0.43，RMSE為10.36 t/hm2，MAE為8.22 t/hm2，RE為0.19%)和RVI、RVI54所構(gòu)建的模型精度(R2為0.47，RMSE為9.73 t/hm2，MAE為7.44 t/hm2，RE為0.09%)。由以上分析可知，應(yīng)用協(xié)同克里格插值估算喬木林地上碳密度時(shí)，選用相關(guān)性較高的自變量因子參與插值，可獲得較好的預(yù)測(cè)效果。

3.3.2地理加權(quán)回歸法

本研究使用地理加權(quán)回歸法構(gòu)建喬木林地上碳密度模型時(shí)，采用與協(xié)同克里格插值相同的自變量組構(gòu)建GWR模型，并對(duì)構(gòu)建的GWR模型進(jìn)行對(duì)比分析。

(1)帶寬選擇

在模型帶寬的選擇上，分別采用固定距離法和自適應(yīng)方法，確定最優(yōu)帶寬。本研究基于植被指數(shù)RVI因子的模型帶寬的計(jì)算結(jié)果見表6和圖2。當(dāng)采用固定距離法時(shí)，模型帶寬為2 514 047 m，模型決定系數(shù)R2為0.46，調(diào)整后R2為0.45；采用自適應(yīng)方法選擇帶寬時(shí)，根據(jù)AIC確定的最佳帶寬可知，隨著帶寬(臨近點(diǎn)個(gè)數(shù))的增加，AIC值在最開始呈現(xiàn)大幅度減小趨勢(shì)，隨后減小的趨勢(shì)逐漸趨于平緩，在臨近點(diǎn)個(gè)數(shù)為30個(gè)時(shí)，AIC值的變化斜率明顯降低，變化幅度減小(圖2)。因此，可確定臨近點(diǎn)個(gè)數(shù)30為最優(yōu)帶寬。對(duì)比兩種方法(表6)可知，基于自適應(yīng)法選擇的帶寬具有較高的調(diào)整決定系數(shù)(0.48)和較低殘差平方和RSS(28 276.96 t2/hm4)，優(yōu)于固定距離法(調(diào)整決定系數(shù)為0.45，RS為36 238.38 t2/hm4)。因此，本研究選擇自適應(yīng)法校準(zhǔn)權(quán)重函數(shù)。由以上分析可知，基于植被指數(shù)RVI構(gòu)建的GWR模型的最優(yōu)臨近點(diǎn)個(gè)數(shù)為30個(gè)，即表示選擇回歸點(diǎn)周邊的30個(gè)點(diǎn)作為核局部帶寬中最為臨近要素的點(diǎn)，參與建模。

表6 GWR模型參數(shù)計(jì)算結(jié)果Tab.6 Parameter estimation results of GWR model

圖2 GWR模型最優(yōu)帶寬Fig.2 Selection of optimal bandwidth for GWR model

(2)模型擬合

在確定最優(yōu)帶寬后，本研究借助ArcGIS 10.1軟件中的地理加權(quán)回歸模塊，分別基于協(xié)同克里格法中采用的2種自變量組構(gòu)建GWR模型。對(duì)95個(gè)樣地碳密度實(shí)測(cè)值進(jìn)行估算和精度評(píng)價(jià)。自變量組為RVI、RVI54、RVI64構(gòu)建的模型具有最高的預(yù)測(cè)精度(R2為0.74，RMSE為6.84 t/hm2，MAE為5.13 t/hm2，RE為0.74%)，優(yōu)于RVI、RVI54構(gòu)建的模型預(yù)測(cè)精度(R2為0.71，RMSE為7.18 t/hm2，MAE為5.37 t/hm2，RE為0.65%)和RVI構(gòu)建的模型預(yù)測(cè)精度(R2為0.58，RMSE為8.63 t/hm2，MAE為6.58 t/hm2，RE為0.21%)，可以獲得較為滿意的估算效果。由此可知，在應(yīng)用GWR構(gòu)建模型時(shí)，隨著自變量個(gè)數(shù)的增加，GWR模型的擬合精度逐漸提高，可以獲得較為滿意的估算效果。因此，在構(gòu)建GWR模型時(shí)，選擇合適的解釋變量對(duì)提高GWR模型預(yù)測(cè)精度具有重要作用。

(3)參數(shù)的非平穩(wěn)性檢驗(yàn)

本研究應(yīng)用ArcGIS 10.1軟件計(jì)算了普通最小二乘(Ordinary least square，OLS)模型和GWR模型的回歸系數(shù)，結(jié)果見表7。根據(jù)文獻(xiàn)[3]的非平穩(wěn)性檢驗(yàn)標(biāo)準(zhǔn)，GWR模型各回歸系數(shù)的第1四分位(25%)和第3四分位(75%)值的變化范圍都大于OLS模型的二倍標(biāo)準(zhǔn)誤值，因此，認(rèn)為喬木林地上碳密度的空間關(guān)系是非平穩(wěn)性的。通過計(jì)算OLS和GWR模型的AIC值，結(jié)果顯示兩種模型的AIC值之差大于3，進(jìn)一步證明了喬木林地上碳密度存在明顯的空間非平穩(wěn)性。

表7 空間平穩(wěn)性檢驗(yàn)Tab.7 Stationary test of relationship

3.3.3兩種建模方法比較

為對(duì)比兩種建模方法的估算效果，本研究計(jì)算了兩種方法預(yù)測(cè)的喬木林地上碳密度的變異系數(shù)，并對(duì)兩種模型的預(yù)測(cè)結(jié)果進(jìn)行對(duì)比分析，采用R2、殘差平方和(Residual sum of squares，RSS)以及均方根誤差(RMSE)、變異系數(shù)(Coefficient of variation，CV)等指標(biāo)對(duì)其進(jìn)行驗(yàn)證(表8)。由表8可知，GWR模型可獲得最好的預(yù)測(cè)效果(R2=0.74，RSS=2 500.11 t2/hm4，RMSE為5.13 t/hm2)，其次為協(xié)同克里格插值(R2=0.47，RSS=8 975.45 t2/hm4，RMSE為9.72 t/hm2)，說明GWR模型的預(yù)測(cè)精度較好。經(jīng)統(tǒng)計(jì)，實(shí)測(cè)碳密度的變異系數(shù)為0.677 7。由表8可知，GWR模型的變異系數(shù)最大，為0.537 2，優(yōu)于協(xié)同克里格插值法的0.496 8，GWR模型預(yù)測(cè)的喬木林地上碳密度保留了0.537 2的空間異質(zhì)性特征。由以上分析可知GWR模型采用了局部回歸方式，在一定程度上提高了模型的預(yù)測(cè)精度，較好地保留了喬木林地上碳密度的空間異質(zhì)性，與協(xié)同克里格插值法相比，具有較高的空間模擬預(yù)測(cè)精度。

表8 兩種模型估算精度比較Tab.8 Comparison of estimation accuracy of two models

兩種模型的喬木林地上碳密度預(yù)測(cè)偏差分布如圖3所示，協(xié)同克里格插值法的預(yù)測(cè)偏差為-19.54～27.79 t/hm2，GWR模型的預(yù)測(cè)偏差為-14.52～18.30 t/hm2，由此可知GWR模型的殘差分布區(qū)間小于協(xié)同克里格插值法，說明GWR模型具有更好的預(yù)測(cè)效果。

圖3 兩種模型的殘差分布Fig.3 Distribution of residuals of two models

4 討論

筆者分別采用了協(xié)同克里格插值和地理加權(quán)回歸對(duì)蓄積量和生物量進(jìn)行了估算[8-9]，結(jié)果顯示地理加權(quán)回歸估算效果優(yōu)于協(xié)同克里格插值法，但其研究?jī)H僅采用了測(cè)樹因子作為自變量進(jìn)行了估算，而這些自變量往往缺少空間連續(xù)性特征。本研究將Landsat 8多光譜影像及DEM數(shù)據(jù)提取的因子作為自變量參與建模，這些衍生的自變量因子具有空間連續(xù)性特征，可以有效避免之前研究中自變量因子不具有空間連續(xù)性特征的缺點(diǎn)，并且可以減少外業(yè)調(diào)查的工作量。由本文的研究結(jié)果可知，地理加權(quán)回歸模型適用于地級(jí)市及區(qū)域尺度的喬木林地上碳密度估算，因?yàn)閱棠玖值厣咸济芏却嬖诳臻g非平穩(wěn)性。本研究采用Landsat 8多光譜影像結(jié)合地理加權(quán)回歸法對(duì)喬木林地上碳密度處理可以獲得較好的估算效果，與前期應(yīng)用ZY-3多光譜影像進(jìn)行測(cè)試估算的效果相比， Landsat 8多光譜影像在空間分辨率以及多光譜波段方面具有更適合估算區(qū)域尺度喬木林地上碳密度的優(yōu)勢(shì)，為區(qū)域森林植被碳儲(chǔ)量和森林參數(shù)的估算提供了較好的數(shù)據(jù)源和估算方法。

準(zhǔn)確估算森林參數(shù)的關(guān)鍵環(huán)節(jié)包括兩部分：高質(zhì)量的遙感數(shù)據(jù)和地面數(shù)據(jù)以及估算方法。在地面數(shù)據(jù)確定的基礎(chǔ)上，遙感數(shù)據(jù)源和估算方法成為決定估算效果的關(guān)鍵因素。首先，在遙感數(shù)據(jù)上，光學(xué)遙感數(shù)據(jù)在森林參數(shù)的估算中得到了廣泛的研究和應(yīng)用，但只反映森林的水平信息，而目前激光雷達(dá)數(shù)據(jù)能夠提供森林垂直信息(如樹高)，并且在很多研究中顯示了很好的適用性[19-20]。因此，在以后的研究和推廣應(yīng)用中可以引入激光雷達(dá)數(shù)據(jù)，將其與光學(xué)遙感數(shù)據(jù)進(jìn)行協(xié)同，估算森林參數(shù)，獲取更高的估算精度。其次，在估算方法上，應(yīng)用較為廣泛的方法包括參數(shù)法、非參數(shù)法和機(jī)理模型3種主要方法[21]，每種方法都有自己的優(yōu)勢(shì)和劣勢(shì)，參數(shù)法大多是通過回歸模型來進(jìn)行估算，非參數(shù)法雖然可提高估算精度，但其無法反映估算變量和自變量間的內(nèi)在關(guān)系，機(jī)理模型雖然考慮了估算變量和自變量間的內(nèi)在機(jī)理關(guān)系，但其模型機(jī)理往往較為復(fù)雜，在中大尺度范圍內(nèi)進(jìn)行應(yīng)用往往需要大量的數(shù)據(jù)支撐。由此可知，應(yīng)用較為簡(jiǎn)單的估算方法獲得精度更高的估算效果成為中大尺度森林參數(shù)估算的關(guān)鍵。

在自變量的提取上，本文采用了Landsat 8多光譜影像的植被指數(shù)及紋理特征，文獻(xiàn)[22]表明，Landsat系列影像的全色波段在遙感定量估算中具有較好的適用性，在后續(xù)研究中，建議進(jìn)一步測(cè)試分析全色波段紋理特征在模型構(gòu)建中的適用性。此外，地形因子與喬木林地上碳密度顯示了較好的相關(guān)性，由此推測(cè)與其相關(guān)的立地條件、氣象因子及生態(tài)因子等因素也可能成為碳密度估算的影響因素，這部分因子對(duì)估算結(jié)果的影響有待進(jìn)一步挖掘分析和研究。在自變量的優(yōu)選上，本研究分兩步完成變量?jī)?yōu)選：一是通過相關(guān)性分析提取自變量，目的是減少共線性因子的影響；二是通過平均殘差平方和準(zhǔn)則進(jìn)一步優(yōu)選自變量，目的是去除對(duì)模型構(gòu)建影響較小的因子，使模型的殘差平方和更小，獲得擬合效果更優(yōu)的模型。本研究結(jié)果顯示，應(yīng)用相關(guān)性分析與平均殘差平方和準(zhǔn)則相結(jié)合的方法來優(yōu)選自變量可獲得較好的估算結(jié)果。

在本研究中還存在以下不足：①受限于數(shù)據(jù)使用權(quán)限，本研究未對(duì)多光譜影像特征對(duì)估算結(jié)果的影響進(jìn)行分析。②本研究假設(shè)影像數(shù)據(jù)獲取時(shí)間和地面數(shù)據(jù)獲取時(shí)間不存在時(shí)間差，同時(shí)受限于地面調(diào)查數(shù)據(jù)量及調(diào)查時(shí)間的限制，只采用了留一交叉驗(yàn)證法進(jìn)行驗(yàn)證。③相關(guān)性分析結(jié)合平均殘差平方和準(zhǔn)則優(yōu)選變量進(jìn)行建模的效果是否優(yōu)于單一變量篩選方法建模的結(jié)果，有待于進(jìn)一步研究分析。在后續(xù)研究中應(yīng)對(duì)存在的不足進(jìn)行深入分析，減小估算結(jié)果的不確定性，并在全國(guó)或更大尺度范圍內(nèi)，結(jié)合地形、環(huán)境、氣候等諸多因素，做進(jìn)一步的研究和推廣應(yīng)用，為區(qū)域及全國(guó)森林植被碳儲(chǔ)量及森林參數(shù)遙感估算提供更優(yōu)的方法，提高估算結(jié)果的精度和可靠性。

5 結(jié)束語

基于Landsat 8多光譜影像和DEM數(shù)據(jù)，提取了遙感及地形因子共計(jì)66個(gè)，結(jié)合地面樣地調(diào)查數(shù)據(jù)，采用皮爾森相關(guān)系數(shù)法和RMSQ法對(duì)自變量進(jìn)行優(yōu)選，分別采用協(xié)同克里格插值和地理加權(quán)回歸法構(gòu)建喬木林地上碳密度模型，并對(duì)結(jié)果進(jìn)行了分析對(duì)比。結(jié)果顯示，地理加權(quán)回歸法構(gòu)建的估算模型精度(R2為0.74，RMSE為6.84 t/hm2，MAE為5.13 t/hm2，RE為0.74%)優(yōu)于協(xié)同克里格插值法(R2為0.47，RMSE為9.72 t/hm2，MAE為7.41 t/hm2，RE為0.12%)，并且較好地保留了估算變量的空間異質(zhì)性，變異系數(shù)分別為0.537 2、0.496 8，說明應(yīng)用Landsat 8多光譜影像結(jié)合地理加權(quán)回歸模型可獲得較好的估算效果。