基于地統計學和多源遙感數據的森林碳密度估算*

嚴恩萍 趙運林 林 輝 莫登奎 王廣興,3

(1.中南林業科技大學 林業遙感大數據與生態安全湖南省重點實驗室 長沙 410004； 2.中南林業科技大學林學院 長沙 410004；3.南伊利諾伊大學地理系 卡本代爾 629012)

基于地統計學和多源遙感數據的森林碳密度估算*

嚴恩萍1,2趙運林1,2林 輝1,2莫登奎1,2王廣興1,2,3

(1.中南林業科技大學 林業遙感大數據與生態安全湖南省重點實驗室 長沙 410004； 2.中南林業科技大學林學院 長沙 410004；3.南伊利諾伊大學地理系 卡本代爾 629012)

【目的】 基于遙感影像空間分辨率和地面樣地大小不一致的現象，采用地統計學和多源遙感數據進行森林碳密度估算，為MODIS數據在區域森林碳密度估算領域的應用提供參考。【方法】 以湖南省攸縣為試驗區，首先利用基于塊的序列高斯協同模擬算法，將25.8 m×25.8 m的樣地數據分別上推到250 m × 250 m、500 m × 500 m和1 000 m×1 000 m； 然后將上推后的樣地數據分別與MOD13Q1、MOD09A1、MOD15A2數據結合，利用序列高斯協同模擬算法開展區域森林碳密度估算研究； 最后將最優結果用于湖南省森林碳密度估算。【結果】 Landsat5和MODIS數據與森林碳密度的敏感因子具有高度相似性，排在前3位的分別為1/TM3、1/TM2、1/TM1和1/Band1、1/Band4、1/Band3； 與植被指數產品MOD13Q1和MOD15A2相比，多光譜數據Landsat5和MOD09A1在攸縣森林碳密度估算方面顯示出巨大潛力，估算精度分別為82.02%和75.64%； 基于MOD09A1的序列高斯協同模擬算法具有很好的適用性，可用于湖南省森林碳密度的空間模擬，估算精度為74.07%。【結論】 采用基于塊的序列高斯協同模擬算法，可以實現由地面樣地到不同空間分辨率MODIS像元之間的轉換； 由于空間分辨率的限制，MOD09A1數據在刻畫空間細節方面不如Landsat5精細。該研究方法適用于地面調查樣地大小和遙感影像空間分辨率不一致的區域森林碳密度估算。

林業遙感； 森林資源清查； 多源遙感； 基于塊的序列高斯協同模擬； 森林碳密度

作為陸地生態系統的主體，森林貯存了陸地生態系統有機碳的76%～98%，在全球碳循環中發揮著重要作用(劉暢等, 2014; Yanetal., 2015)。森林碳蘊含著豐富的地表結構和功能信息，近年來，隨著全球氣候變化和“溫室效應”的加劇，森林碳密度的空間分布和動態變化研究日益受到重視(嚴恩萍等, 2015; Zhaoetal., 2016)。遙感技術的快速發展和國家森林資源清查體系的持續推進，基于遙感影像和地面樣地數據的估算方法日趨成熟，如回歸分析、神經網絡、KNN和空間模擬等(Tianetal., 2012; Wangetal., 2004a; 2009; Luetal., 2012; 戚玉嬌等, 2015)，但受調查成本的限制，目前森林資源調查樣地大小主要分布于10 m × 10 m～50 m × 50 m之間，然而區域森林碳密度制圖采用的遙感影像空間分辨率通常在100 m × 100 m～1 000 m × 1 000 m。如湖南省森林資源清查的樣地大小為25.8 m × 25.8 m，但是Ladsat5和MODIS數據的像元大小分別為30 m × 30 m和250 m × 250 m、500 m × 500 m、1 000 m × 1 000 m，存在明顯的遙感影像空間分辨率和地面樣地大小不一致的現象。因此，需要建立一種尺度轉換關系，將小尺度的地面樣地同大尺度的遙感像元結合起來。

目前，常用的尺度上推方法有最鄰近像元法、窗口平均法和塊克里格法。最鄰近像元法假設較大的塊由N個較小的像元組成，將離塊中心最近的像元值分配給塊，該方法雖然簡單，但是效果不好，無法真實反映較小單元數據間的細微變化； 窗口平均法是一種簡單且使用廣泛的尺度上推方法，即假設窗口內所有像元具有相同的權重，將窗口內像元的平均值賦值給較大的窗口，但忽略了樣本數據間的空間自相關性(Wangetal., 2004b)； 塊克里格法是通過計算塊內較小像元的克里金估計(冉有華等, 2009)，然后計算所有子像元平均值的一種地統計學方法，林業上常用的是序列高斯協同模擬算法。如張茂震等(2009； 2014)將森林資源清查樣地數據與Landsat5影像結合，采用序列高斯協同模擬算法繪制了臨安市和仙居縣森林碳的空間分布，發現模擬結果與地面樣地數據具有較好的空間一致性。沈希等(2011)以臨安市為研究區，比較了序列高斯協同模擬和一元二次非線性回歸算法刻畫森林碳分布的真實程度，發現前者模擬結果更接近地面樣地估算結果。但這些研究多單純采用Landsat5數據，缺少利用序列高斯協同模擬算法結合多源遙感數據開展區域森林碳密度估算的研究，特別是結合MODIS這種多空間分辨率的大尺度數據。

鑒于此，本研究首先以森林資源清查數據和Landsat5影像為信息源，利用基于塊的序列高斯協同模擬算法，將25.8 m × 25.8 m的樣地數據分別上推到250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m， 然后將MODIS數據和上推后的樣地數據結合開展森林碳密度估算，最后將最優結果用于湖南省森林碳密度的估算研究，以期為MODIS數據在區域森林碳密度估算領域的應用提供參考。

1 材料與方法

1.1 研究區概況

研究區位于湖南省東部的攸縣(圖1)，113°09′09″—113°51′30″E，26°46′34″—27°26′30″N。境內四季分明，雨水充足，土壤肥沃，屬中亞熱帶季風濕潤氣候常綠闊葉林帶，年均氣溫17.8 ℃，無霜期292天，年降水量1 410 mm左右。全縣現有林地166 866.7 hm2，其中森林蓄積量和森林覆蓋率分別達到311.87 m3和57.24%。

圖1 研究區位置Fig.1 Geographic position of study area

1.2 數據來源與預處理

1.2.1 固定樣地數據 現有森林碳密度估算一般通過直接或間接測定獲取森林生物量，再乘以生物量中的含碳率推算而得。本文以攸縣2009年森林資源清查數據為基礎，利用生物量回歸方程分樹種(組)計算生物量，再乘以相應樹種的含碳率系數，獲取研究區固定樣地碳密度數據。固定樣地是以4 km × 8 km抽樣間隔設置的正方形樣地，樣地面積0.067 hm2。樣地調查因子包括樣地類別、地類、齡組、優勢樹種、樹高、胸徑等。

采用表1列出的回歸模型(李海奎, 2010)，分樹種(組)計算樣地生物量，混交林生物量分別按比例[6杉(Cunninghamialanceolata)4馬(Pinusmassoniana)、5軟闊5硬闊、3.6杉2.4馬2軟闊2硬闊]計算，經濟林、灌木林分別按平均生物量23.7 t·hm-2和19.76 t·hm-2計算(吳丹等， 2011)，上述3類樹種含碳率均取0.5。對于無明確生物量回歸模型的樹種，采用近似樹種參數替代。

將研究區固定樣地和遙感影像疊加，剔除個別云覆蓋樣地，最后保留78塊樣地作為研究樣本。通過隨機抽樣將樣本數據分為建模樣本和驗證樣本，抽樣比例分別為2/3和1/3； 結合湖南省森林資源分布情況，采用分層抽樣選取1/2的固定樣地作為測算樣本，剔除少量云覆蓋樣地，最后保留2 892塊樣地作為研究樣本，其中建模樣本占3/4，驗證樣本占1/4。1.2.2 遙感數據 采用Landsat5和MODIS 2種遙感影像。Landsat5數據空間分辨率為30 m × 30 m，含云量均低于1%，影像質量較好，用于攸縣森林碳密度遙感反演的接收時間為2009年8月21日(軌道號123/41)，用于湖南省森林碳密度反演的Landsat5數據接收時間依次為2009年8月21日(軌道號123/40、123/41、123/42)、2009年8月28日(軌道號124/40、124/41、124/42)和2009年7月18日(軌道號125/40、125/41、125/42)。數據預處理包括輻射定標、大氣校正和幾何校正3種操作； 數據變換包括倒數運算、比值運算(兩波段、三波段和四波段組合間的比值運算)、植被指數運算(NDVI、EVI、ARVI和SAVI 4種指數)、主成分變換和紋理變換(均值、角二階矩、對比度、相關、相異、熵、逆差距和方差8種紋理因子)，參與運算的波段包括TM1、TM2、TM3、TM4、TM5和TM7共6個原始波段，共計生成88種基于Landsat5的遙感變量。

①WS表示樹干生物量，WP表示樹皮生物量，WB表示樹枝生物量，WL表示樹葉生物量，WT表示地上部分總生物量，WR表示地下部分生物量，D表示樣地平均胸徑，H表示樣地平均樹高。WS: stem biomass;WP: bark biomass;WB: branch biomass;WL: leaf biomass;WT: lotal aboveground biomass;D: tree diameter at breast height;H: tree height.

MODIS數據包括MOD13Q1、MOD09A1和MOD15A2共3種產品，空間分辨率分別為250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m，接收時間均為2009年8月21日(軌道號h28v06，用于湖南省森林碳密度反演的MODIS數據軌道號包括h27v05、h27v06和h28v06)，3種產品的預處理包括投影轉換、影像鑲嵌、影像裁剪以及圖像修復4個步驟，均通過遙感軟件ENVI 5.0實現。其中圖像修復采用IDL語言編程實現(崔麗華等， 2009； Xiangetal., 2013)，包括條帶去除和去云處理2部分。同Landsat5，MOD09A1數據最后選取Band3、Band4、Band1、Band2、Band6、Band7共6個波段分別進行倒數運算、比值運算、植被指數運算、主成分變換和紋理變換，共計生成88種基于MOD09A1的遙感變量； MOD13Q1和MOD15A2均包括2種植被指數產品，分別為NDVI、EVI和FPAR、LAI，遙感變量均包括原始植被指數和16種紋理因子。

為篩選合適的遙感因子，研究運用ArcGIS 10.2軟件提取地面樣地所在位置的遙感因子，利用SPSS 20.0軟件分析遙感變量和樣地森林碳密度值間的Pearson相關性，保留相關性最高的遙感因子參與后續模擬。

1.2.3 土地利用數據 根據研究區實際情況，參照《土地利用現狀分類標準》(GB/T 21010—2015)，將研究區土地利用/覆蓋分為森林、耕地、水體、草地、建設用地和其他用地6大類。在ENVI 5.0遙感軟件支持下，利用先驗知識建立訓練樣本，對處理后的Landsat5遙感影像進行監督分類，得到初始分類結果，然后結合森林資源連續清查數據和森林資源分布圖，對監督分類結果進行修正，得到覆蓋湖南省和攸縣的2009年土地利用/覆蓋數據。經檢驗，總體精度分別為85.1%和89.5%，Kappa系數均超過0.79，精度較高，滿足研究的需要。最后提取分類結果中的森林信息，將其作為掩膜數據，用于后續序列高斯協同模擬和森林碳密度制圖。

1.3 研究方法

1.3.1 序列高斯協同模擬 序列高斯協同模擬(sequential Gaussian co-simulation, SGCS)是以地面樣地和遙感影像為基礎，通過結合變異函數，利用隨機模擬算法估計未知參數的方法。該算法假設研究區由n個等大小的像元構成，每個像元的估計值從已有估計值和周圍樣地數據確定的條件累積分布中通過隨機抽樣獲得。該條件分布由一個統計平均數和方差確定，統計平均數和方差可以基于已有的地面樣地數據和遙感影像通過點位協同簡單克里格估計實現。具體計算公式如下：

(1)

(2)

式中：zsck(u)代表協同簡單克里格算法像元位置u處森林碳密度的預測值；δ2(sck)(u)代表協同簡單克里格算法的預測方差；z(uα)代表樣本數據，α=1,2,…,n(u)；n(u)代表給定搜索范圍內獲得的樣本數量；y(u)代表像元u處的光譜變量；λsckα、λscky分別代表協同簡單克里格算法中樣本數據和影像數據的權重；mz、my分別代表地面樣本數據和遙感影像數據的均值；Czz(0)代表地面樣本數據的方差；Czy(0)代表森林碳密度與遙感影像光譜變量的協方差； 且當h=uα-u時，Czy(h)代表森林碳密度預測值與光譜變量的交叉協方差函數。

當輸出的森林碳密度分布圖像元與固定樣地大小一致時，SGCS模擬過程為： 1) 用隨機抽樣方法設置一種遍歷每個像元的順序； 2) 從分布中隨機抽取一個像元位置值u，采用協同簡單克里格算法分別估計像元位置u處的預測值和方差，由此預測值和方差確定一個條件累積分布函數； 3) 從分布中隨機抽取一個值，將其作為隨機變量在像元位置u的實現。重復步驟1～3，直至所有像元都有估計值，進而得到一張覆蓋整個研究區的森林碳密度分布圖。將該過程執行M次，就可得到M張分布圖，通過求和取平均值獲取研究區森林碳密度的均值分布圖和方差分布圖，本研究的估計次數M全部取250。具體程序執行過程，通過修改Wang等(2004a)提出的空間估計算法實現。變異函數是序列高斯協同模擬中最重要的參數之一，研究采用VARIOWIN 3.0(Pannatier, 1996)分析樣地森林碳密度，采用Spherical標準化模型進行擬合，從中選擇最優變異函數模型。

1.3.2 基于塊的序列高斯協同模擬 為解決固定樣地大小(25.8 m × 25.8 m)與MODIS空間分辨率(250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m)不一致的問題，本研究采用Wang等(2004a)提出的尺度轉換方法——基于塊的序列高斯協同模擬(sequential Gaussian block co-simulation, SGBCS)，首先將25.8 m × 25.8 m的地面樣地分別上推到250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3種尺度，實現地面樣地數據和MODIS像元空間分辨率的匹配； 然后將上推后的樣地數據與MODIS數據結合，利用序列高斯協同模擬算法依次開展250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3種尺度的森林碳密度估算； 最后將最優算法用于湖南省森林碳密度估算，探討其在大區域森林碳密度估算領域的適用性。地面樣本的劃分同Landsat5數據，2/3數據用于建模，1/3數據用于驗證。為深入探討基于塊的序列高斯協同模擬算法的適用性，本研究將其用于湖南省森林碳密度的估算研究，其中3/4(2 169塊樣地)數據作為建模樣本，余下的1/4數據(723塊樣地)作為驗證樣本。

研究包括4種形式的森林碳密度模擬： 1) 基于固定樣地和攸縣Landsat5影像的SGCS(30 m × 30 m)； 2) 基于固定樣地和攸縣Landsat5影像的SGBCS(250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m)； 3) 基于上推樣地和攸縣MODIS影像的SGCS(250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m)； 4) 基于上推樣地和湖南省MODIS影像的SGCS(500 m × 500 m)。

(3)

式中：cov(ui,uj)表示像元ui和uj的協方差，從該累積分布中，隨機抽取一個值作為該塊的模擬值；σ2(sck)表示像元ui的協同克里格方差。

基于塊的協同克里格模擬中，每個塊的條件累積分布取決于塊內像元的預測值和克里格方差以及二者之間的協方差。基于塊的協同克里格模擬算法克服了已有尺度上推算法的缺陷，可以較好地運用連續分布的遙感影像開展森林碳密度的空間估計研究，算法的效果取決于用戶估計的精度、方差和運算時間。

1.3.3 精度評價 基于固定樣地的森林碳密度估算在由單木水平、樣地水平推算到區域水平的過程中，存在很大的不確定性，忽略這些不確定性將導致區域森林碳密度的高估或低估。本研究采用判定系數(R2)、均方根誤差(root mean square error, RMSE)和相對誤差(relative error, RE)3個指標，分別對Landsat5和MODIS像元水平的森林碳密度模擬結果進行不確定性分析。

(4)

(5)

(6)

R2可反映估測值與對應實測值之間的趨勢線擬合程度； RMSE雖能很好地反映估測模型的可靠性，但其是一個絕對數，大小與均值相關，無法直觀反映結果的準確程度；RE具有相對性，不僅考慮了樣本估計值與實測值之間誤差的大小，同時也兼顧了樣本本身的大小，其值越小模型的估測精度越高。

2 結果與分析

2.1 空間變異函數

研究針對攸縣25.8 m × 25.8 m、250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 4種尺度的樣地數據進行半方差函數擬合(圖2)。變程表示空間變異和自相關函數的最大距離，觀察圖2可知模型的變程分別為15.62、14.91、17.04和17.76 km，在此距離之內，觀測值是空間相關的； 在此范圍之外，觀測值本質上相互獨立。

圖2 森林碳密度數據的標準化變異函數γ(h)Fig.2 Spatial autocorrelation γ(h) of forest carbon using spherical model for standardized datah表示距離h is distance(km).

其標準化Spherical模型為：

(7)

(8)

(9)

(10)

研究區樣地間距在東西向和南北向分別為8 km和4 km，森林碳密度空間變異的動態范圍(變程)大于樣地間距，樣本具有足夠的空間信息量。

2.2 樣地數據分析

分析研究區4種尺度的樣地碳密度統計結果(表2)可知，攸縣25.8 m × 25.8 m樣地碳密度在0 ～ 39.27 t·hm-2之間變化，標準差為8.39，且變異系數大于1.00，屬于強變異，說明攸縣森林碳密度數值分布合理，相互之間存在較大差異。值得注意的是，隨著樣地空間尺度增大，樣地碳密度平均值、標準差、變異系數呈逐漸減小的趨勢，說明研究區的森林樣地空間分布不均勻。

表 2 攸縣固定樣地碳密度統計結果①

① SD: Standard deviation; CV:Coefficient of variation.

2.3 相關性分析

采用SPSS 20.0軟件分別計算樣地碳密度(包括25.8 m × 25.8 m固定樣地和3種上推樣地)與遙感變量間的Pearson相關性，結果發現基于Landsat5和MOD09A1的遙感變量與樣地碳密度的相關性較高，分別在-0.455～0.497和-0.626～0.763之間變化，表3列出了相關系數排在前6位的遙感變量。對于Landsat5數據，當顯著水平為0.01時，與森林碳密度達到顯著相關的變量有1/TM3、TM1var、Elevation、1/TM2等15個因子，其中1/TM1、1/TM2、1/TM3和TM2mean4個變量相關系數在0.450以上； 對于MOD09A1數據，當顯著水平為0.01時，與森林碳密度相關系數達到顯著的變量有1/Band1、1/Band3、1/Band7、NDVI等28個因子，其中1/Band1、1/Band3 和1/Band4 3個變量相關系數在0.700以上。

表3 樣地森林碳密度與遙感變量的Pearson相關性①

① **： 在0.01水平(雙側)上顯著相關，*： 在0.05水平(雙側)上顯著相關； TM2mean表示TM2灰度共生矩陣的平均值。Correlation is significant at the 0.01 level ** and the 0.05 level*; TM2meanis the average value of gray-level co-occurrence matrix； Band67=Band6/Band7.

MOD13Q1和MOD15A2數據衍生的遙感變量與森林碳密度的相關系數均較低，分別介于-0.230～0.312和-0.390～0.395之間，當顯著水平為0.01時，二者與森林碳密度相關系數達到顯著的變量分別為NDVImean和LAImean、LAIhom、LAIdis、FPARmean、LAIcon，其中相關系數最高的因子為NDVImean和LAImean。因此，研究依次選擇相關系數最高的4個因子： 1/TM3(Landsat5數據)、NDVImean(MOD13Q1數據)、1/Band1(MOD09A1數據)和LAImean(MOD15A2數據)，參與攸縣森林碳密度的空間模擬，相關系數分別為0.497、0.312、0.763和0.395，顯著水平均達0.000。

2.4 基于Landsat5的森林碳密度模擬

基于固定樣地和Landsat5影像，利用序列高斯協同模擬算法導出攸縣30 m × 30 m像元水平的森林碳密度分布圖，同時采用驗證樣本對其進行精度評價。結果表明，模型預測值和實測值之間具有良好的線性擬合關系(圖3a)，決定系數R2達0.82，估測精度為82.02%，且殘差圖(圖3b)散點分布均勻。張茂震等(2014)基于浙江省仙居縣固定樣地和Landsat5影像，采用序列高斯協同模擬算法開展了森林碳估算研究，結果表明估算效果較好，進一步佐證了該算法的適用性。

圖3 基于Landsat5的攸縣森林碳密度模擬值精度驗證Fig.3 Accuracy assessments for simulated values of forest carbon density in You county based on Landsat5a.預測值與實測值比較Comparison of predicted and observed values; b.殘差分析Residual analysis.

為直觀反映模擬結果與地面樣地值的吻合程度，研究將森林碳密度模擬結果與地面樣地值疊加(圖4)，結果發現攸縣森林碳密度整體上與地面樣地實測值的分布趨勢保持一致： 東部由于經營管理完善，分布大面積質量較好的針葉成熟林，森林碳密度高； 中部和西南部主要用于農業種植和城鎮建設，森林碳密度相對較低。采用ArcGIS 10.0軟件進行統計得出，攸縣2009年平均森林碳密度達15.06 t·hm-2，具有較高的實用性和可信度。

研究采用的固定樣地大小為25.8 m × 25.8 m，但用于大尺度反演的遙感數據MODIS像元大小分別為250 m × 250 m、500 m × 500 m、1 000 m ×1 000 m，本文在前期數據處理的基礎上，通過引入Wang等(2004a)設計的基于塊的序列高斯協同模擬，結合Landsat5數據和地面樣地數據，將地面調查樣地大小從25.8 m × 25.8 m分別上推到250 m × 250 m、500 m × 500 m 和1 000 m × 1 000 m(圖5)。基于尺度上推的序列高斯協同模擬程序，通過修改Dr Journel的協同克里格程序和Dr Almeida的協同模擬程序實現(Pannatier, 1996; Deutschetal., 1992)。

為進一步驗證Landsat5數據上推結果的有效性，研究采用ArcGIS10.0軟件空間分析模塊的“波段集統計”工具，分別計算3種尺度上推結果與30 m × 30 m模擬結果重采樣后森林碳密度的相關性，結果顯示相關系數分別為0.996、0.987和0.964，說明上推效果較好，可用于后續基于MODIS的碳密度估算。

圖4 攸縣基于Landsat5數據和序列高斯協同模擬算法的森林碳密度空間預測指標與樣地實測值的比較Fig.4 Spatial distribution of predicted indexes for forest carbon compared with the plot values at spatial resolution of 30 m×30 m by sequential Gaussian co-simulation using Landsat5 image in You countya.預測均值Predicted values； b.預測均方差Predicted variance； c.預測值大于均值概率Probability for predicted values larger than sample means.

圖5 攸縣森林碳密度基于塊的序列高斯協同模擬預測結果Fig.5 Spatial distributions of simulated forest carbon density using sequential Gaussian block co-simulation algorithm and Landsat5 imagea. 250 m × 250 m; b.500 m × 500 m; c.1 000 m × 1 000 m

2.5 基于MODIS的森林碳密度模擬

利用固定樣地點分別提取250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3種水平的森林碳密度，同時結合MOD13Q1、MOD09A1和MOD15A2影像，開展相應空間分辨率的森林碳密度估算。表4列出了3種MODIS數據森林碳密度模擬值的精度檢驗結果。

分析表4可知，MOD09A1數據的模擬效果最好，R2達0.78，其次是MOD13Q1數據(R2=0.57)，MOD15A2數據效果最差(R2僅0.53)； 均方根誤差(RMSE)和估測精度與決定系數(R2)表現一致，其中模擬誤差分別為4.65、2.56和4.26 t·hm-2，估測精度分別為58.58%、75.64%和54.84%。3種檢驗指標共同表明多光譜數據MOD09A1的估算結果明顯優于植被指數產品MOD13Q1和MOD15A2，在區域森林碳密度估算方面顯示出明顯的優勢。

同Landsat5數據，研究基于MOD09A1影像的森林碳密度模擬結果包括3部分： 碳密度模擬平均值分布(圖6a)、碳密度模擬平均方差分布(圖6b)、碳密度模擬值大于均值的概率分布(圖6c)，空間分辨率為500 m × 500 m。

表4 3種MODIS數據森林碳密度模擬值的精度檢驗結果

圖6 攸縣基于MOD09A1影像和序列高斯協同模擬算法的森林碳密度空間預測指標分布Fig.6 Spatial distribution of predicted indexes for forest carbon density using the sequential Gaussian co-simulation algorithm based on MOD09A1 image in You Countya.預測均值Predicted values; b.預測均方差Predicted variances; c.預測值大于均值概率Probability for predicted values larger than sample means.

分析可知，基于MOD09A1數據模擬的森林碳密度空間分布(圖6a)與研究區內基于Landsat5模擬的森林碳密度具有較好的一致性，即樣地模擬值較高的局部區域，相應MODIS數據的模擬值也較高。這說明利用MODIS數據結合基于塊的序列高斯協同模擬算法，可以真實再現研究區森林碳密度的空間分布格局，模擬結果的方差分布(圖6b)和模擬結果大于平均數的概率分布(圖6c)，與森林碳密度模擬結果的平均值具有較好的空間一致性。

2.6 湖南省森林碳密度模擬

為進一步探討基于塊的序列高斯協同模擬算法的適用性，研究將其用于湖南省森林碳密度的估算研究。首先借助ArcGIS10.0軟件的sample工具，抽取研究區3/4數據(2 169塊樣地)開展基于Landsat5影像和固定樣地的森林碳密度估算研究； 然后采用預留的1/4數據(723塊樣地)進行精度驗證(圖7)。結果表明，模型預測值和實際觀測值之間具有良好的線性擬合關系，決定系數(R2)達0.84，估測精度為87.02%，且殘差散點分布均勻，說明模型的擬合效果較好，可用于后續基于MOD09A1的森林碳密度估算。

在前期數據處理的基礎上，采用基于塊的序列高斯協同模擬算法，將湖南省固定樣地從25.8 m × 25.8 m上推到500 m × 500 m，利用固定樣地點提取相應水平的森林碳密度值(圖8a)。同時結合MOD09A1影像，采用序列高斯協同模擬算法開展湖南省森林碳密度的估算研究； 樣本數據的劃分同Landsat5。

采用決定系數(R2)、均方根誤差(RMSE)和估測精度(EA)3個指標檢驗湖南省森林碳密度的估算精度，具體檢驗結果見表5。分析可知，基于MOD09A1數據的湖南省森林碳密度模擬效果較好，決定系數(R2)達0.75，估測精度為74.06%，進一步佐證了基于塊的序列高斯協同模擬算法的有效性，可用于區域森林碳密度的估算。

圖7 湖南省森林碳密度樣地模擬值精度驗證Fig.7 Accuracy assessments for simulated values of forest carbon density in Hunan Provincea.預測值與實測值比較Comparison of predicted and observed values; b.殘差分析Residual analysis.

圖8 湖南省森林碳密度空間分布Fig.8 Spatial distribution map of forest carbon density in Hunan Provincea.上推樣點空間分布Spatial distribution for up-scaling plots using spatial block co-simulation algorithm;b.模擬值空間分布 Spatial distribution for simulated values using sequential Gaussian co-simulation algorithm.

表5 基于MOD09A1數據的湖南省森林碳密度精度檢驗結果

圖8b顯示了基于MOD09A1影像的湖南省森林碳密度分布，分析可知，湖南省森林碳密度與地面樣地上推值保持一致的分布趨勢： 高值區主要分布在湘東、湘南的郴州、湘西的懷化、吉首、張家界等地，湘中地區的森林碳密度較小。因為研究區東、南、西三面環山，中部分割為丘陵性河谷盆地，北部是沉積平原和吞噬湖泊，特殊的地形條件使得湘東、西、南三面山區分布著大面積的森林植被，由于經營管理完善，森林質量較好，碳密度大； 湘北和湘中地區主要用于農田種植和水體覆蓋，森林碳密度相對較小。沅江以及資水附近的大面積連續森林碳密度高值分布，充分反映了森林植被隨地表水分布的規律； 而中南部和東北部洞庭湖附近的大面積連續低值分布，則反映出人類活動對森林植被覆蓋的影響，森林碳密度整體上與研究區植被分布規律一致。該結論進一步說明利用基于塊的序列高斯協同模擬算法可以再現湖南省森林碳密度的空間分布格局。

3 討論

3.1 相關性分析

作為反映變量之間相關性密切程度的統計指標，森林碳密度和遙感因子之間的相關系數與遙感影像的空間分辨率和傳感器有關。本研究中4種遙感數據(Landsat5、MOD13Q1、MOD09A1和MOD15A2)衍生的變量均與森林碳密度具有較高的相關性，且相關系數均在0.312以上。Landsat5衍生變量與森林碳密度的相關系數在-0.455～0.497之間波動，且1/TM3與森林碳密度的相關系數最高(達0.497)，這與Wang等(2014)的研究結果一致(相關系數0.312)； 比較而言，MOD09A1衍生變量與森林碳密度的相關系數在-0.626～0.763之間波動，其中1/Band1與森林碳密度的相關系數最高。從某種程度上說，這一結論與Landsat5數據的研究結果相似，主要原因是Landsa5的 TM3波段與MOD09A1的Band1波段具有相似的波譜區間，分別在0.63～0.69 μm和0.62～0.67 μm之間波動。此外，MOD13Q1數據的衍生變量NDVImean和MOD15A2數據的LAImean分別與森林碳密度的相關性最高，因此二者分別用于250 m × 250 m和1 000 m × 1 000 m 2種空間分辨率的森林碳密度制圖。

3.2 不同尺度森林碳密度模擬

為克服地面樣地大小和遙感影像空間分辨率不一致的缺陷，研究采用Wang等(2004a)設計的基于塊的序列高斯協同模擬算法，將樣地碳密度從25.8 m × 25.8 m分別上推到250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3種尺度，結果發現該算法可以很好地刻畫森林碳密度的空間分布格局。同250 m × 250 m的MOD13Q1和1 000 m × 1 000 m的MOD15A2產品相比，多光譜數據Landsat5和MOD09A1在森林碳密度模擬方面顯示出明顯的優勢，因為來自可見光、近紅外和中紅外波段的光譜變量可以提取指示森林碳密度信息的光譜變量(Chenetal., 2002)。Landsat5數據的模擬結果與樣地實測值之間具有很好的一致性，且由其和固定樣地數據上推的250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3種尺度的森林碳密度精度可靠，可作為相應尺度森林碳密度模擬的真值數據。雖然由MOD09A1數據模擬的森林碳密度在區域和國家尺度是可以接受的，但不能用于局部地區森林碳密度的高精度模擬，主要是因為MOD09A1數據在刻畫空間細節方面不如Landsat5數據精細； 另外，受Landsat5影像質量和序列高斯協同模擬算法數據量的限制，本文只是通過隨機抽樣獲取部分樣地上推到500 m × 500 m，然后將上推后的樣地與MOD09A1數據結合，開展湖南省森林碳密度的空間模擬，利用Landsat5和固定樣地數據開展湖南省森林碳密度的高精度模擬將是下一步的研究重點。

3.3 不同信息源森林碳密度模擬

研究采用基于塊的序列高斯協同模擬算法，將地面樣地(包括250 m × 250 m的地面樣地和3種上推后的地面樣地)分別與不同空間分辨率的遙感影像結合，開展多種空間尺度的森林碳密度模擬，該算法克服了地面樣地大小與遙感影像空間分辨率不一致的缺陷。雖然MODIS數據提供的模擬結果可用于大尺度的森林碳密度制圖，但會損失部分森林碳密度的估計精度和空間變異性。比較而言，Landsat5數據在模擬森林碳密度的空間分布和變異性方面更為精確，比MOD09A1數據呈現的森林冠層結構更為復雜，這與Yang等(2007)的研究結論一致。另外，Landsat5的TM3波段和MOD09A1的Band1波段在提供森林碳密度的預測變量方面顯示出較大的潛力。

3.4 森林碳密度模擬精度的影響因素分析

正確估算森林碳密度是研究全球氣候變化的關鍵途徑。本研究應用序列高斯協同模擬算法對攸縣和湖南省森林碳密度的空間分布進行了估算，同利用光譜值和森林參數建立的回歸模型相比，序列高斯協同模擬算法能夠更多地考慮森林參數同光譜值之間的非線性依賴關系(Trotteretal., 1997)。SGCS是基于像元的估計方法，且模擬結果的最小值接近0，更能真實反映非林地的實際情況(沈希等, 2011)。森林碳密度的模擬精度受多種因素的影響，如模擬次數L、影像特征波段的選取和采樣密度等：1) 相關研究表明，模擬次數L對模擬結果精度有重要影響，模擬次數越多，其平均值就越接近該位置的數學期望(郭含茹等, 2016)。本文當L達到250后，各像元位置森林碳密度的估計值趨于穩定，因此將其作為最終的模擬次數。2) 由于影響SGCS森林碳密度最優波段的選擇往往隨研究區和數據源的不同而變化，因此最優波段需要通過具體試驗才能獲取。本文采用的Landsat5和3種MODIS數據，雖然研究區相同，但是最優波段因數據源的不同而變化。3) 眾所周知，隨著采樣密度的增加，抽樣誤差隨之減小，但到達一定程度后，逐漸趨于穩定。本文主要研究森林資源清查樣地與不同傳感器數據相結合的模擬，然而關于不同采樣密度下的森林碳密度估計考慮不多。郭含茹等(2016)通過分析4種采樣密度對浙江省仙居縣森林碳密度模擬精度的影響，指出利用SGCS、SGBCS和森林資源清查樣地開展森林碳密度估計時，可適當降低對采樣密度的要求。

4 結論

隨著遙感技術的發展，基于固定樣地數據的森林碳密度遙感制圖得到了廣泛應用，然而地面樣地大小和遙感影像空間分辨率不一致的缺陷阻礙了制圖精度的提高。鑒于此，本研究以地面固定樣地和Landsat5影像為數據源，采用基于塊的序列高斯協同模擬算法，將25.8 m × 25.8 m的地面樣地分別上推到250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m； 然后將上推樣地和MODIS影像結合，探討MODIS數據在區域森林碳密度遙感估算領域的潛力。結果表明： 1) 采用基于塊的序列高斯協同模擬算法，可以實現由地面樣地到不同空間分辨率MODIS像元之間的轉換； 2) Landsat5和MODIS數據與森林碳密度的敏感因子具有高度的相似性，排在前3位的分別為TM3、TM2、TM1和Band1、Band4、Band3； 3) 同植被指數產品MOD13Q1和MO15A2相比，多光譜數據Landsat5和MOD09A1在模擬森林碳密度方面顯示出巨大潛力，因為來自可見光、近紅外和中紅外波段的光譜變量可以提取指示森林碳密度信息的遙感因子，模擬精度分別82.02%和75.64%，但受空間分辨率的限制，MOD09A1數據在刻畫空間細節方面不如Landsat5精細； 4) 為進一步探討基于塊的序列高斯協同模擬方法的適用性，研究將其推廣到基于MOD09A1數據的湖南省森林碳密度模擬，結果表明估算效果較好，估測精度為74.06%，可用于區域森林碳密度的模擬。

崔麗華, 劉善軍, 張艷博, 等. 2009. IDL程序去除MODIS遙感影像鑲嵌中的黑色條帶. 河北理工大學學報： 自然科學版，31(3): 79-82.

(Cui L H, Liu S J, Zhang Y B,etal. 2009. Using IDL to remove black belt from the mosaic of MODIS remote sensing images. Journal of Hebei Polytechnic University： Natural Science, 31(3): 79-82. [in Chinese])

郭含茹, 張茂震, 徐麗華, 等. 2016. 不同采樣密度下縣域森林碳儲量仿真估計. 生態學報, 36(14): 4373-4385.

(Guo H R, Zhang M Z, Xu L H,etal. 2016. Simulation of regional forest carbon storage under different sampling desities. Acta Ecologica Sinica, 36(14): 4373-4385. [in Chinese])

李海奎. 2010. 中國森林植被生物量和碳儲量評估. 北京： 中國林業出版社.

(Li H K. 2010. Estimation and evaluation of forest biomass carbon storage in China. Beijing: China Forestry Publishing House. [in Chinese])

劉 暢, 李鳳日, 甄 貞. 2014. 空間誤差模型在黑龍江省森林碳儲量空間分布的應用. 應用生態學報, 25(10): 2779-2786.

(Liu C, Li F R, Zhen Z. 2014. Prediction of spatial distribution of forest carbon storage in Heilongjiang Province using spatial error model. Chinese Journal of Applied Ecology, 25(10): 2779-2786. [in Chinese])

戚玉嬌, 李鳳日. 2015. 基于KNN 方法的大興安嶺地區森林地上碳儲量遙感估算. 林業科學, 51(5): 46-55.

(Qi Y J, Li F R. 2015.Remote sensing estimation of aboveground forest carbon storage in Daxing’an Mountains based on KNN method. Scientia Silvae Sinicae, 51(5): 46-55. [in Chinese])

冉有華, 李 新. 2009. 基于塊克里金的土壤水分點觀測向像元尺度的尺度上推研究. 冰川凍土, 31(2): 275-283.

(Ran Y H, Li X. 2009. Up Scaling of point soil moisture measurements to oixel averages based on block Kriging. Journal of Glaciology and Geocryology, 31(2): 275-283.[in Chinese])

沈 希, 張茂震, 祁祥斌. 2011. 基于回歸與隨機模擬的區域森林碳分布估計方法比較. 林業科學, 47(6): 1-8.

(Shen X, Zhang M Z, Qi X B. 2011. Comparison of regional forest carbon estimation methods based on regression and stochastic simulation. Scientia Silvae Sinicae, 47(6): 1-8.[in Chinese])

吳 丹, 邵全琴, 劉紀遠, 等. 2011.1985—2030年江西泰和縣森林植被碳儲量的時空動態. 應用生態學報, 2011, 22(1): 41-46.

(Wu D, Shao Q Q, Liu J Y,etal. 2011. Spatiotemporal dynamics of forest carbon storage in Taihe County of Jiangxi Province in 1985—2030. Chinese Journal of Applied Ecology, 22(1): 41-46. [in Chinese])

嚴恩萍, 林 輝, 王廣興, 等. 2015. 基于MODIS混合像元分解的湖南省森林碳密度反. 應用生態學報, 26(11): 3433-3442.

(Yan E P, Lin H, Wang G X,etal. 2015.Estimation of human forest carbon density based on spectral mixture analysis of MODIS data. Chinese Journal of Applied Ecology, 26(11): 3433-3442. [in Chinese])

張茂震, 王廣興, 葛宏立, 等. 2014. 基于空間仿真的仙居縣森林碳分布估算. 林業科學, 50(11): 13-22.

(Zhang M Z, Wang G X, Ge H L,etal. 2014. Estimation of forest carbon distribution for Xianju county based on spatial simulation. Scientia Silvae Sinicae, 50(11): 13-22. [in Chinese])

張茂震, 王廣興, 周國模, 等. 2009. 基于森林資源清查、衛星影像數據與隨機協同模擬尺度轉換方法的森林制圖. 生態學報, 29(6): 2919-2928.

(Zhang M Z, Wang G X, Zhou G M,etal. 2009. Mapping of forest carbon by combining forest inventory data and satellite images with co-simulation based up-scaling method. Acta Ecologica Sinica, 29(6): 2919-2928. [in Chinese])

Chen J M, Pavlic G, Brown L,etal. 2002. Derivation and validation of Canada-wide coarse-resolution leaf area index maps using high-resolution satellite imagery and ground measurements. Remote sensing of environment, 80(1): 165-184.

Deutsch C V, Journel A G. 1992. GSLIB: geostatistical software library and user’s guide. 2nded. Oxford University Press.

Lu D S, Chen Q, Wang G X,etal. 2012. Aboveground forest biomass estimation with Landsat and lidar data and uncertainty analysis of the estimates. International Journal of Forestry Research, Article ID 436537, 16 pages, doi:10.1155/2012/436537.

Pannatier Y. 1996. Variowin: software for spatial statistics analysis in 2D. New York, Springer.

Tian X, Su Z, Chen E,etal. 2012. Estimation of forest above-ground biomass using multi-parameter remote sensing data over a cold and arid area. International Journal of Applied Earth Observation and Geoinformation, 14(1): 160-168.

Trotter C, Dymond J, Goulding C. 1997. Estimation of timer volume in a coniferous planation forest using Landsat TM. International Journal of Remote Sensing, 18(10): 2209-2223.

Wang G X, Gertner G Z, Anderson A B. 2004a.Spatial-variability-based algorithms for scaling-up spatial data and uncertainties. IEEE Transactions on Geoscience and Remote Sensing, 42(9): 2004-2015.

Wang G X, Gertner G Z, Anderson A B. 2004b. Mapping vegetation cover change using geostatistical methods andbi-temporal Landsat TM images. IEEE Transactions on Geoscience and Remote Sensing, 42(3): 632-643.

Wang G X, Tonny O, Zhang M Z,etal. 2009. Mapping and spatial uncertainty analysis of forest vegetation carbon by combining national forest inventory data and satellite images. Forest Ecology and Management, 258(7): 1275-1283.

Wang G X, Zhang M Z. 2014. Upscaling with conditional co-simulation for mapping above-ground forest carbon. Scale Issues in Remote Sensing,doi:10.1002/9781118801628.ch06.

Xiang H B, Liu J S, Cao C X,etal. 2013. Algorithms for moderate resolution imaging spectroradiometer cloud-free image compositing. Journal of Applied Remote Sensing, 7(1): 073486.doi:10.1117/1.JRS.7.073486.

Yan E P, Lin H, Wang G X,etal. 2015. Improvement of forest carbon estimation by integration of regression models and spectral unmixing of Landsat data. IEEE Geoscience and Remote Sensing Letters, 12(9): 2003-2007.

Yang P, Ryosuke S, Wu W B,etal. 2007. Evaluation of MODIS land cover and LAI products in cropland of North China Plain using in situ measurements and Landsat images. IEEE Transactions on Geoscience and Remote Sensing, 45(10): 3087-3097.

Zhao P, Lu D, Wang G,etal. 2016. Examining spectral reflectance saturation in Landsat imagery and corresponding solutions to improve forest aboveground biomass estimation. Remote Sensing, doi:10.3390/rs8060469.

(責任編輯 石紅青)

Estimation of Forest Carbon Density Based on Geostatistics and Multi-Resource Remote Sensing Data

Yan Enping1,2Zhao Yunlin1,2Lin Hui1,2Mo Dengkui1,2Wang Guangxing1,2,3

(1.KeyLaboratoryofForestryRemoteSensingBigData&EcologicalSecurityforHunanProvinceCentralSouthUniversityofForestry&TechnologyChangsha410004; 2.CollegeofForestry,CentralSouthUniversityofForestry&TechnologyChangsha410004; 3.DepartmentofGeography,SouthernIllinoisUniversityCarbondaleIl629012USA)

【Objective】 Due to the inconsistency of spatial resolutions between sample plots and image pixels, in this study the estimation of forest carbon density was conducted by using geostatistics method and multi-resource remote sensing data, which aims to provide reference for the application of MODIS data in the regional estimation of forest carbon density.【Method】 Firstly, the spatial block co-simulation algorithm was employed to scale up the sample plots of forest carbon density in You county of Hunan Province from the spatial resolution of 25.8 m×25.8 m to the spatial resolutions of 250 m × 250 m, 500 m × 500 m and 1 000 m × 1 000 m respectively. Then, MODIS images with three spatial resolutions corresponding to those mentioned above, were applied to map forest carbon density for this county using sequential Gaussian co-simulation algorithm. Finally, the best model was applied in the estimation of forest carbon density for Hunan Province.【Result】 There were highly similarities for sensitive factors of forest carbon density between Landsat5 and MOD09A1 data, according to theresult of Pearson product moment correlations, the top three sensitive factors were 1/TM3, 1/TM2, 1/TM1 for Landsat5 and 1/Band1, 1/Band4, 1/Band3 for MOD09A1 respectively; compared to the vegetation product of MOD13Q1 and MOD15A2, multi-spectral data such as Landsat5 and MOD09A1 showed great potential in the simulation of forest carbon density with the accuracy of 82.02% and 75.64%, respectively; there is good application for the sequential Gaussian co-simulation algorithms based on the image of MOD09A1, which can be used in the spatial simulation of forest carbon density for Hunan Province with the accuracy of 74.07%.【Conclusion】 The spatial block co-simulation algorithm could be used to realize the conversion of spatial resolutions from sample plots to the pixels of MODIS images. It was also found that the MODIS derived maps were more smoothed than those from Landsat5 due to the limitation of spatial resolutions, especially in the terms of capturing the spatial variability of forest carbon density. The adopted method was well suited for regional estimation of forest carbon density based on the combination of forest inventory sample plot data and remotely sensed images, especially for the areas that the plot sizes and images pixels were inconsistent.

forestry remote sensing； forest resource inventory； multi resource remote sensing； sequential Gaussian block co-simulation； forest carbon density

10.11707/j.1001-7488.20170708

2016-03-21；

2016-09-15。

國家“十二五”高技術發展研究計劃項目(2012AA102001)； 國家自然科學基金面上項目(31470643)； 林學重點學科開放基金擬資助項目(2016YB08)。

S718.5; S757.2

A

1001-7488(2017)07-0072-13

*林輝為通訊作者。