結合地統(tǒng)計學與Landsat 8影像的喬木林地上碳儲量估算

邱新彩，鄭冬梅，王海賓，安天宇，許等平，郝月蘭，彭道黎

（1.北京林業(yè)大學 林學院，北京 100083；2.國家林業(yè)和草原局 調查規(guī)劃設計院，北京 100714；3.國家林業(yè)和草原局 林產(chǎn)工業(yè)規(guī)劃設計院，北京 100010）

森林是全球陸地生態(tài)系統(tǒng)中最大的有機碳庫，在維護生態(tài)環(huán)境和碳平衡中發(fā)揮巨大作用[1]。森林碳儲量是研究大氣碳收支的重要參數(shù)，也是評價國家固碳潛力是否有效的數(shù)據(jù)支持[2]。準確評估森林碳儲量，不僅可以推動碳循環(huán)和碳匯研究工作的進展，而且可為森林資源監(jiān)測和經(jīng)營管理、可持續(xù)發(fā)展和生態(tài)環(huán)境建設提供決策服務[3]。傳統(tǒng)的森林碳儲量估算方法，如采用森林資源一類、二類調查數(shù)據(jù)或生態(tài)樣地的測樹資料來估算碳儲量，存在費時費力、調查間隔期長、適用范圍小等缺點，遙感數(shù)據(jù)的全覆蓋、多時相、多空間分辨率等諸多優(yōu)勢特征，為碳儲量估算提供可行路徑，與傳統(tǒng)估算方法相比更具優(yōu)勢[4]。它可以與地面實測樣地數(shù)據(jù)相結合，建立碳儲量遙感估算方程，為快速、準確、實時地估算多尺度森林碳儲量提供了有效手段[5-7]。

目前國內外關于森林碳儲量研究主要集中在國家或區(qū)域尺度上，普遍通過直接或間接測定森林植被或優(yōu)勢樹種的生物量，再乘以對應的碳元素含量推算而得[8-9]。Landsat 8 作為新型的Landsat 陸地資源衛(wèi)星數(shù)據(jù)，提供了光譜信息豐富、時間分辨率較高、獲取速度快的遙感信息源，在林業(yè)監(jiān)測中獲得了較好的效果[10-12]，但應用Landsat 8 影像估算森林生物量和碳儲量的研究還不多。森林生物量和碳儲量的估算方法多以參數(shù)和非參數(shù)方法為主，其中參數(shù)方法以回歸分析為主[13-14]，非參數(shù)方法包括隨機森林、人工神經(jīng)網(wǎng)絡、K 最近鄰等的研究和應用較多[15-20]。但是這些方法大多都沒有考慮空間異質性對森林生物量和碳儲量的影響，缺乏明確的地理空間信息，被認為是非空間方法。地統(tǒng)計學因其兼顧了空間結構特征并具備完好的空間制圖功能，被廣泛地應用在地質學、生態(tài)學等領域。近年來，國內外學者采用不同的地統(tǒng)計學方法對森林變量進行空間插值及驗證分析，在森林蓄積量和生物量研究中取得了較好的估算效果[21-23]。賀鵬[24]、王海賓[25]等采用了協(xié)同克里格插值法對森林蓄積量、森林生物量進行了空間插值研究，結果均顯示了較好的適用性，獲得了較高的預測精度。在文章中，作者均未采用遙感數(shù)據(jù)結合地統(tǒng)計學方法對蓄積量和生物量進行適用性研究，但作者對后續(xù)采用遙感數(shù)據(jù)結合地統(tǒng)計學方法估算變量進行了展望。因此，本文借鑒前人的研究展望，將遙感數(shù)據(jù)與協(xié)同克里格插值法相結合來估算喬木林地上碳儲量，試驗并討論這種方法的可行性和適用性。此外，為與非空間估算方法對比，本文將常用的穩(wěn)健估計方法作為建模方法參與模型構建，并對兩種估算方法的精度進行對比分析。

1 研究區(qū)概況

浙江?。?80°01′～123°10′E，27°06′～31°11′N）地處中國東南沿海，全省陸地面積1.02×105km2，地勢從西南向東北呈階梯狀下降趨勢，西南以山地為主，平均海拔800 m；中部多為丘陵，海拔在100～500 m；東北部是低平的沖積平原，海拔在10 m 以下。浙江省屬亞熱帶季風性氣候，四季分明，雨熱同季，年平均氣溫15～18 ℃，無霜期230～270 d，年平均日照時數(shù)1 710～2 100 h，年 均 降 水 量980～2 000 mm。土壤類型多樣，主要有紅壤、黃壤、粗骨土、水稻土、潮土、鹽土等類型。常綠闊葉林為該省的地帶性植被。浙江省森林資源豐富，森林面積為6.04×106hm2，占林地面積的91.70%，森林覆蓋率達59.43%，森林蓄積達2.81×108m3，占活立木總蓄積的89.58%。

2 數(shù)據(jù)準備及處理

2.1 專項調查數(shù)據(jù)

2012年，通過外業(yè)調查獲取專項調查數(shù)據(jù)，依據(jù)林分年齡（分為幼齡林、中齡林、近成過熟林）和優(yōu)勢樹種作為控制，設置和調查了96 個樣地，盡量使其具有代表性，樣地為0.067 hm2的正方形樣地，主要包括針葉林、闊葉林、混交林等不同森林類型，專項調查樣地分布見圖1。樣地屬性因子包括樹種、胸徑（起測胸徑為5 cm）、平均高、郁閉度、坡向、坡度、海拔、株數(shù)、年齡、喬木林地上生物量和喬木林地下生物量等因子。本研究用到的因子為喬木林地上生物量因子，樣地調查數(shù)據(jù)的生物量統(tǒng)計結果見表1。

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

表1 調查樣地生物量統(tǒng)計結果Table 1 Statistic results of biomass for the survey samples in the study area

依據(jù)獲取的樣地蓄積量，乘以生物量換算因子（Biomass expansion factor，BEF）得到喬木林地上生物量，為便于計算，將喬木林地上生物量轉換為每公頃生物量（t·hm-2），每公頃生物量乘以相對應的含碳率即可得到喬木林地上碳儲量。蓄積量直接采用已計算好的樣地蓄積量數(shù)據(jù)。BEF采用沈楚楚[24]提供的浙江省4 大樹種（組）的BEF，分別為杉木林（0.745 3 t·m-3）、馬尾松林（0.883 9 t·m-3）、軟闊類（1.657 2 t·m-3）、硬闊類（1.070 5 t·m-3）。將其他樹種分別歸并到軟闊類和硬闊類樹種（組）并采用相應的BEF 進行生物量計算。本文依據(jù)《土地利用、土地利用變化與林業(yè)碳匯計量監(jiān)測技術指南》確定樹種含碳率，指南中的樹種含碳率可覆蓋研究區(qū)內的各個優(yōu)勢樹種（組）。

在獲得喬木林地上碳儲量數(shù)據(jù)后，采用二倍標準差法剔除異常值，結果剩余95 個專項調查數(shù)據(jù)，將其作為因變量用于本文Landsat 8 影像構建喬木林地上碳儲量模型。

2.2 輔助數(shù)據(jù)

本研究采用的輔助數(shù)據(jù)為覆蓋研究區(qū)的數(shù)字高程模型（Digital elevation model，DEM）數(shù)據(jù)，比例尺為1:50 000，用于提供地形等因子。

2.3 遙感數(shù)據(jù)及預處理

本研究選用一景Landsat 8 OLI（Operational Land Image，陸地成像儀傳感器）影像作為遙感數(shù)據(jù)源，覆蓋整個研究區(qū)，空間分辨率為30 米，影像軌道號為119/040，采集時間為2014年10月26日。所采用波段包括波段序號為2～7 的6 個波段，應用ENVI5.3 軟件對Landsat 8 影像進行預處理，包括輻射定標、大氣校正。

3 研究方法

3.1 建模變量設置

本研究選用植被指數(shù)（9 個）、紋理波段（48個）、主成分變換因子（3 個）、纓帽變換因子（3個）和地形因子（3 個）共66 個特征作為自變量參與模型構建。按照地面實測樣地位置，提取各個樣地對應的自變量因子。表2為本研究采用的9個植被指數(shù)。

紋理是一個非常復雜的空間屬性，易受到傳感器觀測角度以及觀測對象所處環(huán)境的不同而產(chǎn)生較大的變化，此外紋理特征的度量易受到窗口大小、方向以及步長的影響。為探究Landsat 8 多光譜波段紋理特征的應用潛力，本研究采用灰度共生矩陣方法分別提取Landsat8 多光譜波段的8種紋理特征，8 種紋理特征的計算見表3。

表2 研究中用到的植被指數(shù)?Table 2 Vegetation indices used in the study

表3 研究中用到的紋理特征?Table 3 Texture features used in the study

影像的多光譜波段間常存在較高的相關性，圖像變換可以將圖像信息壓縮到少數(shù)幾個有效波段，去除冗余信息。目前較為常用的有主成分分析法和纓帽變換法，在以往的研究中，主成分分析法的前3 個波段和纓帽變換的前3 個波段應用較為廣泛[10,27]。本研究采用主成分分析獲取的前3 個波段和纓帽變換的前3 個波段（亮度Brightness、綠度Greenness 和濕度Wetness）共計6 個波段變換因子作為自變量，通過分析確定相關自變量參與碳儲量模型的構建。此外，應用DEM 數(shù)據(jù)提取對應的海拔、坡度、坡向3 個因子，參與模型構建。

3.2 建模變量優(yōu)選

采用皮爾森相關性分析法對自變量進行篩選，提取出極顯著相關的因子，進一步采用方差膨脹因子（Variance inflation factor，VIF）法篩選自變量，將篩選的自變量用于模型構建。VIF 是指解釋變量之間存在多重共線性時的方差與不存在多重共線性時的方差之比，為容忍度的倒數(shù)，VIF 越大，顯示共線性越嚴重[28]。經(jīng)驗判斷方法表明：當0＜VIF＜10，不存在多重共線性；當10≤VIF＜100，存在較強的多重共線性；當VIF≥100，存在嚴重多重共線性。依據(jù)相關的文獻資料[29-30]，本研究設置0＜VIF＜10，容忍度＞0.1 來進行控制，對自變量做進一步提取。

3.3 穩(wěn)健估計

穩(wěn)健估計是實際情況偏離假設模型的不敏感性估計，主要目的是在建模中尋求樣本數(shù)據(jù)分布的主體，不遷就異常樣地，采用一定的抗差函數(shù)削弱其使用比例，以此來孤立異常樣地和抵制其不利影響，從而使估測模型具有一定的“穩(wěn)定性”和“抗干擾性”[31-32]。穩(wěn)健估計實質上是求非線性函數(shù)的極值解，采用迭代函數(shù)法求解，具體步驟及原理詳見文獻[32]。

3.4 協(xié)同克里格插值

協(xié)同克里格法（Co-kriging，COK）在變量估測上具有較好的預測精度[24-25,33]，它是利用兩個或兩個以上的變量，將其中一個作為主變量，剩余變量作為輔助變量，將主變量的自相關性和主輔變量的交互相關性結合起來進行無偏最優(yōu)估計[30]，詳細的理論介紹可參考文獻[24]和[25]。協(xié)同克里格插值的公式如下：

式中，Z*(x0)為待估點處的蓄積量估測值；λ1i和λ2i分別是主變量Z1和輔變量Z2實測值的權重，且n和m分別是參與估測x0點的蓄積量和其它輔助變量Z2的實測值數(shù)目。

3.5 模型驗證

考慮到覆蓋單景Landsat8 影像的專項調查樣地數(shù)量較少的情況，本文采用留一交叉驗證法進行評價。通過對比模型預測的喬木林地上碳儲量與實測碳儲量值，比較不同建模方法的預測精度，以評價不同模型的估算性能。本研究采用的評價指標有：決定系數(shù)（R2）、均方根誤差（Root mean square error，RMSE）、平均絕對偏差（Mean absolute error，MAE）和總預報偏差的相對誤差（RE），R2越大，模型的擬合精度越好，RMSE、MAE 和RE 值越小，表示模型預測精度越高。4 個評價指標計算公式如下：

式中：yi為碳儲量實測值，為碳儲量預測值，為碳儲量實測值平均值，n為樣本個數(shù)。

4 結果與分析

4.1 描述性統(tǒng)計分析

使用SPSS 20.0 統(tǒng)計軟件，對單位面積喬木林地上的碳儲量進行了描述性統(tǒng)計，并對其進行K-S檢驗，結果顯示喬木林地上碳儲量數(shù)據(jù)服從正態(tài)分布。表4是單位面積喬木林地上碳儲量的描述性統(tǒng)計和正態(tài)分布性檢驗結果。

表4 喬木林地上碳儲量描述性統(tǒng)計結果Table 4 Descriptive statistics of carbon storage in the arbor forest

4.2 自變量篩選結果

應用SPSS 20.0 軟件，對提取的專項調查數(shù)據(jù)的喬木林地上碳儲量和66 個自變量因子進行皮爾森相關性分析，篩選出與碳儲量極顯著相關（P＜0.01）的自變量因子共計22 個，見表5，包括NDVI、RVI、RVI54、RVI64、SAVI、NLI、ARVI、PCA2、Greenness、海拔、B1_Mean、B2_Mean、B3_Mean、B3_Entropy、B3_Second、B3_Correlation、B6_Mean、B6_Variance、B6_Contrast、B6_Diss、B6_Entropy、B6_Second。

對提取出的極顯著相關的22 個自變量因子進一步優(yōu)選，生成基于VIF 方法篩選的自變量集，篩選后的7 個自變量因子見表6，包括RVI、NLI、海拔、B1_Mean、B3_Correlation、B6_Mean、B6_Variance。

4.3 穩(wěn)健估計

將篩選的自變量集構建喬木林地上碳儲量模型，應用MATLAB2014a 計算各個模型的參數(shù)，最后得到穩(wěn)健估計法構建的喬木林地上碳儲量估算模型，其數(shù)學表達式為C=0.000 1×RVI+2.100 6×NLI+0.000 1×海拔-0.000 2×B1_Mean+0.000 5×B3_Mean-0.000 4×B6_Mean-0.000 5×B6_Variance，C 為喬木林地上碳儲量，其他因子為對應的自變量因子。

表5 自變量因子與喬木林地上碳儲量的Pearson 相關系數(shù)?Table 5 Pearson correlation coefficients between independent variables and arbor forest above ground carbon storage

表6 基于VIF 法篩選的自變量因子?Table 6 Independent factors filter based on the VIF method

4.4 協(xié)同克里格插值

基于VIF 法篩選的自變量，提取出與喬木林地上碳儲量相關性最高的1 個自變量（RVI）和相關性最低的2 個自變量（B6_Mean、B3_Correlation），采用協(xié)同克里格插值對喬木林地上碳儲量進行估算。為取得更好地估算效果，本研究對ArcGIS 10.1 地統(tǒng)計分析模塊中的11 種變異模型進行了計算和分析，根據(jù)變異函數(shù)的理論和評價方法[25]，得到最優(yōu)變異函數(shù)的擬合結果，變異函數(shù)擬合結果見表7。由表7可知，在協(xié)同克里格插值中，采用提取的3 個變量參與插值后，塊金值為0，塊金值/基臺值小于25%，表明系統(tǒng)具有很強的相關性，依據(jù)變異函數(shù)的評價標準可知，在使用RVI、B6_Mean、B3_Correlation 時，K-貝塞耳模型可獲得最優(yōu)的擬合結果。為獲得較高的預測精度，本研究選擇VIF 中的RVI、B6_Mean、B3_Correlation 3 個變量擬合的最優(yōu)變異函數(shù)（K-貝塞耳模型）作為最優(yōu)變異函數(shù)，最后采用協(xié)同克里格法對喬木林地上碳儲量進行插值，并對插值效果進行評價。

表7 基于VIF 提取的3 個自變量擬合的最優(yōu)方差函數(shù)模型及參數(shù)Table 7 Optimal variance function model and parameters fitted by three different variables extracted by VIF

4.5 模型評價及對比

應用2 種方法對喬木林地上的碳儲量進行估算，并采用4 種評價指標對預測結果進行比較，結果見表8。由表8可知，協(xié)同克里格插值法的決定系數(shù)（R2=0.45）高于穩(wěn)健估計（R2=0.42），表示協(xié)同克里格插值的擬合精度較好。協(xié)同克里格插值法的均方根誤差、平均絕對偏差和總預報偏差的相對誤差（RMSE=9.88 t·hm-2，MAE=7.75 t·hm-2，RE=0.24%）均明顯低于穩(wěn)健估計法（RMSE=10.15 t·hm-2，MAE=8.03 t·hm-2，RE=0.27%），表示協(xié)同克里格模型預測精度更高。由以上分析可知，協(xié)同克里格插值法對喬木林地上碳儲量的估算性能更強，具有較高的空間模擬預測精度。

表8 2 種建模方法的預測精度比較Table 8 Comparison of prediction accuracy of two models

5 結論與討論

本研究基于Landsat 8 影像和DEM 數(shù)據(jù)提取植被指數(shù)、紋理特征、主成分波段、纓帽變換以及地形因子，結合研究區(qū)專項調查樣地數(shù)據(jù)和輔助數(shù)據(jù)，采用皮爾森相關系數(shù)和VIF法提取自變量，基于篩選的自變量，采用穩(wěn)健估計和協(xié)同克里格插值法構建喬木林地上碳儲量模型，并對模型精度進行對比分析。結果表明協(xié)同克里格插值法構建模型的精度要優(yōu)于穩(wěn)健估計。

在自變量篩選上，采用了兩步法對自變量進行了篩選，以期為建模提供更好的自變量因子。皮爾森相關系數(shù)法是用來反映因變量和自變量間的相關程度，通過剔除不相關的自變量因子，提取出與因變量相關性較高的自變量因子用于建模。但在實際的應用中，采用遙感影像衍生出的自變量因子間往往存在多重共線性問題，如何減小這些問題對模型穩(wěn)定性的影響，也是變量篩選的關鍵因素之一。方差膨脹因子法在這方面顯示了其優(yōu)勢，通過考察給定的解釋變量被方程中其他所有解釋變量所解釋的程度，以此來判斷其是否存在多重共線性，通過一定的標準剔除影響較大的自變量因子，減小多重共線性影響。因此，本文將兩種方法結合使用，既兼顧了因變量與自變量間的相關性，又兼顧了自變量間的共線性問題，在一定程度上可提高所構建模型的穩(wěn)定性和適用性。

研究結果顯示協(xié)同克里格插值法的估算效果要優(yōu)于穩(wěn)健估計法，這是因為協(xié)同克里格插值法將主變量的自相關性和主輔變量的交互相關性結合起來，對喬木林地上碳儲量進行無偏最優(yōu)估計，可獲得較好的估算精度。此外，本文將協(xié)同克里格插值和遙感影像相結合的方式，考慮了遙感自變量的空間連續(xù)特征，與前人只考慮不具有空間連續(xù)性特征的自變量進行插值相比，更具有優(yōu)勢，可以說是對文獻[22]和[23]的探索式研究。

采用協(xié)同克里格插值法構建喬木林地上碳儲量模型時，本研究選擇了相關性最高的RVI 自變量和相關性最低的B6_Mean、B3_Correlation 兩個自變量，這是因為在采用與喬木林地上碳儲量相關性最好的3 個自變量插值時，由于變量間存在的共線性問題，導致插值精度較低，在試驗不同3 個自變量插值時，也顯示了同樣的效果。當選擇RVI、B6_Mean、B3_Correlation 進行插值時，可獲得最高的精度。因此，在采用協(xié)同克里格插值法時，選取相關性較高且共線性較弱的自變量因子對因變量進行插值具有較好的估測效果和適用性。

在自變量因子上，地形因子與喬木林地上碳儲量顯示了較好的相關性，在經(jīng)過兩步篩選后，海拔因子仍被優(yōu)選出來作為建模變量，說明海拔因子能夠反映喬木林地上碳儲量的空間垂直分布特征。但本研究在采用協(xié)同克里格插值法構建模型時，并未將海拔因子參與建模。本研究采用的DEM 數(shù)據(jù)比例尺為1:50 000，屬于較小比例尺數(shù)據(jù)，在一定程度上均化了地理空間特征信息，未能充分挖掘出地理空間特征信息對建模精度的貢獻。因此，如何進一步挖掘地形因子信息、采用適宜比例尺DEM 數(shù)據(jù)作為建模自變量，充分利用其特有的空間特征信息，提高所構建模型的穩(wěn)定性和精度等問題有待于進一步研究。