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

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

(1.國家林業和草原局林產工業規劃設計院， 北京 100010; 2.國家林業和草原局調查規劃設計院， 北京 100714)

0 引言

森林參數的準確估算是林學研究的熱點之一[1-2]，對森林經營管理、生物多樣性保護、全球氣候變化研究具有重要意義。掌握森林參數現狀及分布情況，有利于森林經營管理者科學管理森林，對在不同層面制定森林經營管理政策、調整經營管理策略具有重要作用[3]。

近年來，遙感技術得到了快速發展，因即時、快速和低成本等特點而使其在國土、農林、水利、礦業等領域得到了廣泛應用，取得了較好的應用效果[4-5]。尤其是在區域和大尺度范圍的監測上，光學遙感數據覆蓋范圍廣、獲取周期短、數據成本低廉等，具有激光雷達、微波雷達等其他遙感數據無法比擬的優勢[6]。因此，應用光學遙感技術結合地面樣地數據對森林參數進行估算成為必要的手段[7]。

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

為挖掘Landsat 8多光譜影像在全國森林植被碳儲量遙感估算中的應用潛力，同時探索更好的估算方法，本研究以位于浙江省內的一景Landsat 8影像覆蓋范圍為研究區，以喬木林地上碳密度為研究對象，基于研究區內的地面樣地調查數據，通過計算獲得喬木林地上碳密度數據并將其作為因變量，采用協同克里格插值法和地理加權回歸方法，構建喬木林地上碳密度模型，并對建模結果進行對比分析，以期為區域及全國范圍的喬木林地上碳密度及森林參數遙感快速估算提供更好的數據源和估算方法，并提高遙感估算結果的精度和可靠性。

1 研究區及數據獲取

1.1 研究區概況

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

圖1 研究區內調查樣地分布Fig.1 Distribution of survey samples in study area

1.2 數據獲取

1.2.1地面數據

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

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

表1 地面樣地喬木林地上生物量統計結果Tab.1 Statistic results of aboveground biomass of arbor forest for ground sample t/hm2

將每公頃生物量乘以樹種含碳率即可得到單位面積喬木林地上碳密度。樹種含碳率參考《土地利用、土地利用變化與林業碳匯計量監測技術指南》，指南中的樹種含碳率可覆蓋研究區內的各個優勢樹種。

在獲得喬木林地上碳密度數據后，采用二倍標準差法剔除異常值，結果剩余95個地面樣地數據，將其作為因變量參與喬木林地上碳密度模型的構建。

1.2.2遙感影像及輔助數據

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

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

2 研究方法

2.1 空間自相關分析

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

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

(1)

式中n——樣本總數

I——莫蘭指數

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

wij——空間位置i和j的權重

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

2.2 自變量提取及優選

2.2.1自變量提取

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

表2 植被指數Tab.2 Vegetation indices

2.2.2自變量優選

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

2.3 建模方法

2.3.1協同克里格插值法

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

(2)

其中

式中Z*(x0)——待估點處的蓄積量估測值

λ1i、λ2i——主變量Z1和輔變量Z2實測值的權重

n、m——參與估測x0點的蓄積量和輔助變量Z2的實測值數目

2.3.2地理加權回歸法

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

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

(3)

式中x、y——空間樣本的坐標值

N——樣本數0——截距

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

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

采用地理加權回歸法構建模型時，模型帶寬的選擇至關重要，在很大程度上決定模型的預測準確度[3,14]。本研究采用最小化信息準則(Akaike information criterion，AIC)[9]來確定最佳帶寬：當AIC值取最小值時，即可獲得最優帶寬及地理加權回歸模型。AIC值計算公式為

AIC=-2lnH+2k

(4)

式中k——未知參數的數量

H——似然函數AIC——AIC值

最優模型帶寬在AIC值最小時獲得。

此外，在完成GWR模型構建后，需進一步對GWR模型進行參數的非平穩性檢驗，檢驗模型效果。地理加權回歸理論和GWR模型參數的非平穩性檢驗詳見文獻[3，14]。

2.4 模型精度評價

本文采用留一交叉驗證法對建模結果進行評價。采用的評價指標包括：決定系數(R2)、均方根誤差(Root mean square error，RMSE)、平均絕對誤差(Mean absolute error，MAE)和總預報偏差的相對誤差(Relative error，RE)。

3 研究結果

3.1 描述性統計及空間自相關性分析

3.1.1描述性統計分析

使用SPSS 20.0統計軟件，采用K-S檢驗對數據進行檢驗。經分析，喬木林地上碳密度數據服從正態分布，K-S檢驗值為0.374，其平均值為19.61 t/hm2，標準差13.29 t/hm2，中位數17.91 t/hm2，最小值為0.88 t/hm2，最大值為56.02 t/hm2。

3.1.2空間自相關分析

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

3.2 自變量優選結果

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

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

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

在相關性分析的基礎上，本研究采用RMSQ法對22個自變量因子進一步優選，結果見表4。經篩選后，RVI、RVI54、RVI64、Greenness、B3_Correlation、B6_Mean共計6個因子被優選出來。

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

3.3 模型構建

3.3.1協同克里格插值法

(1)模型擬合

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

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

為0，塊金值與基臺值的比值小于25%，表明系統具有很強的相關性。依據變異函數的評價標準可知，在使用RVI進行協同克里格插值時，孔穴效應模型可獲得最優的變異函數擬合效果，在使用RVI、RVI54插值時，K-貝塞耳模型可獲得最優的擬合結果，在使用RVI、RVI54、RVI64時，K-貝塞耳模型可獲得最優的擬合結果。

(2)精度評價

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

3.3.2地理加權回歸法

本研究使用地理加權回歸法構建喬木林地上碳密度模型時，采用與協同克里格插值相同的自變量組構建GWR模型，并對構建的GWR模型進行對比分析。

(1)帶寬選擇

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

表6 GWR模型參數計算結果Tab.6 Parameter estimation results of GWR model

圖2 GWR模型最優帶寬Fig.2 Selection of optimal bandwidth for GWR model

(2)模型擬合

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

(3)參數的非平穩性檢驗

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

表7 空間平穩性檢驗Tab.7 Stationary test of relationship

3.3.3兩種建模方法比較

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

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

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

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

4 討論

筆者分別采用了協同克里格插值和地理加權回歸對蓄積量和生物量進行了估算[8-9]，結果顯示地理加權回歸估算效果優于協同克里格插值法，但其研究僅僅采用了測樹因子作為自變量進行了估算，而這些自變量往往缺少空間連續性特征。本研究將Landsat 8多光譜影像及DEM數據提取的因子作為自變量參與建模，這些衍生的自變量因子具有空間連續性特征，可以有效避免之前研究中自變量因子不具有空間連續性特征的缺點，并且可以減少外業調查的工作量。由本文的研究結果可知，地理加權回歸模型適用于地級市及區域尺度的喬木林地上碳密度估算，因為喬木林地上碳密度存在空間非平穩性。本研究采用Landsat 8多光譜影像結合地理加權回歸法對喬木林地上碳密度處理可以獲得較好的估算效果，與前期應用ZY-3多光譜影像進行測試估算的效果相比， Landsat 8多光譜影像在空間分辨率以及多光譜波段方面具有更適合估算區域尺度喬木林地上碳密度的優勢，為區域森林植被碳儲量和森林參數的估算提供了較好的數據源和估算方法。

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

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

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

5 結束語

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