森林生物量遙感降尺度研究

劉沁茹,孫 睿,*

1 遙感科學國家重點實驗室,北京師范大學地理科學學部,北京 100875 2 北京市陸表遙感數據產品工程技術研究中心,北京師范大學地理科學學部,北京 100875

森林生物量作為森林生態系統長期生產代謝過程中積累的產物,在評價森林固碳能力方面發揮著重要的作用,其大小和分布格局與生物多樣性、碳氧平衡、氣候變化等密切相關。目前,結合遙感數據建立定量反演模型已成為獲取森林生物量產品的主要方式。GLAS星載激光雷達可以獲取全球森林冠高,并用于森林生物量估算,目前已有利用GLAS數據生成的全球或區域1 km分辨率森林生物量產品,如Saatchi等[1]于2011年生成的1 km分辨率全球熱帶地區森林地上生物量產品、Baccini等[2]生成的1 km分辨率非洲森林地上生物量產品。但低空間分辨率產品不能細致地呈現森林生物量的空間分布格局,不能滿足小區域森林調查及動態變化監測的需要,而地面調查與機載激光雷達數據獲取成本高,工作量大。因此,有必要結合較低分辨率森林生物量產品與易獲取的較高分辨率多光譜或SAR數據進行森林生物量產品降尺度方法研究,提高森林生物量產品的空間分辨率。

降尺度方法被廣泛用于土壤學、氣候學等領域[3],它的實質是將大尺度數據分解為區域尺度數據,實現亞像元信息提取[4]。目前常用的遙感數據降尺度方法有3種：數理統計降尺度、基于調制分配降尺度和基于光譜混合模型降尺度。數理統計回歸分析通過建立某些地表參數(如植被指數、地表反照率等)和不同分辨率影像值的回歸關系實現降尺度[4-14],核心即為“關系尺度不變”。Kustas等[7]提出DisTrad方法,首次將歸一化差值植被指數NDVI用于地表溫度LST的降尺度。隨后植被覆蓋度[8]、EVI[9]和地表反照率α[11,13]等因子也逐漸被引入降尺度模型。基于調制分配降尺度目前主要用于地表溫度,一些學者通過 PBIM(像元強度調制)算法實現了降尺度[15-16]。該算法利用高分辨率反射波段影像值,對應識別低分辨率熱像元中的地形變化,根據反向陰陽坡的觀測照度差異校正低分辨率熱像元溫度,但顯然這種方法只適用于白天的LST產品降尺度,且地表類型單一。Nichol[17]在此基礎上提出了適用于夜間且地形平坦區域的EM(Emissivity Modulation)算法。基于光譜混合模型降尺度考慮了各端元組分的差異,通過建立關聯高低分辨率影像的回歸關系實現降尺度[15]。

近年來,地表溫度、土壤濕度等領域的尺度轉換方法研究成果豐富,但針對森林生物量產品降尺度的研究尚不多見。已有的研究也多是基于森林植被分布圖,給出各植被類型的生物量代表值,將其作為各植被類型在大尺度上的平均生物量降尺度到格網上,得到高分辨率生物量數據[18-19],這并非是傳統意義上的生物量測算值[20]。且這些研究多是將遙感數據和地面調查數據結合,而利用高分辨率多光譜數據對低分辨率森林生物量進行統計學降尺度的研究較少。因此,本研究擬建立多光譜地表參數和低分辨率森林地上生物量產品之間的統計關系,探究高低分辨率森林地上生物量之間的轉換辦法,生成高分辨率森林地上生物量產品。

1 研究區概況

本文選擇了美國馬里蘭州的兩個研究區(圖1)。研究區1位于藍嶺森林區的南山(South Mountain),海拔在500 m以上,77°20′—77°38′W,39°24′—39°41′N,地表覆蓋類型主要包括落葉闊葉林、農用地等；氣溫表現為冬季低于0℃而夏季保持在20℃左右。研究區2位于切薩皮克灣附近的平原地區,76°37′—76°59′W,38°58′—39°19′N,地表覆蓋類型主要包括混交林、農用地植被混合和城市,此外還有小范圍的落葉闊葉林；氣候狀況表現為夏熱冬溫,年平均降水量約為1000 mm。研究區1和研究區2分別位于山區和平原,地形狀況差異大；且兩個研究區森林分布格局有差異,研究區1的森林分布相較于研究區2更加密集。因此,本文對兩個研究區進行了對比研究。

2 數據及研究方法

2.1 數據及數據處理

本文使用到的數據主要包括CMS(Carbon Monitoring System)森林地上生物量產品、GEOCARBON森林地上生物量產品、Landsat-5 TM數據和ASTER GDEMV2(ASTER Global Digital Elevation Model V2)產品。

CMS森林地上生物量產品(數據下載網址https://daac.ornl.gov/)覆蓋的空間范圍為美國馬里蘭州,空間分辨率為30 m,時間分辨率為2011年一次性估測。該森林地上生物量產品以機載激光雷達數據和實地采樣數據作為原始數據,基于隨機森林回歸模型運算生成[21],單位為Mg/hm2。經投影轉換后,根據研究區經緯度范圍對數據進行裁剪。

GEOCARBON森林地上生物量產品(數據下載網址https://www.wur.nl/en.htm)空間分辨率為0.01°(約1 km),該產品通過融合Avitabile等[22]的熱帶地區生物量分布圖和Santoro等[23]的北方森林生物量分布圖生成,產品僅覆蓋森林地區,單位為Mg/hm2。其中,北緯30°以北的森林地上生物量是通過ENVISAT ASAR數據反演森林蓄積量GSV,然后結合IPCC(Intergovernmental Panel on Climate Change)的生物量轉換與擴展因子(BCEF)以及材積與生物量轉換關系估算得到的。

Landsat-5 TM數據(行號:015,列號:033；下載于http://earthexplorer.usgs.gov/)成像時間為2011年8月22日,影像包含7個波段,研究中用到的波段1—5(可見光、近紅外和中紅外波段)和波段7(中紅外波段)空間分辨率為30 m。研究區內影像無云量。數據級別為L1T,尚未進行輻射定標和大氣校正,因此采用FLAASH大氣校正模塊做進一步處理,最終生成反射率數據。

研究中使用的GDEMV2產品(數據下載于http://www.gscloud.cn/)是基于ASTER數據計算生成的,數據時期為2009年,空間分辨率為30 m。根據研究區經緯度范圍選擇4幅DEM影像進行拼接(N38W077/N38W078/N39W077/N39W078),并進行投影轉換。利用ArcGIS軟件生成坡度分布圖并裁剪。

2.2 森林地上生物量降尺度方法

本文的降尺度思路有兩種：第一種是對1 km分辨率GEOCARBON AGB產品直接降尺度；第二種采用先升尺度后降尺度的方法,對30 m分辨率CMS AGB產品進行升尺度模擬生成的低分辨率生物量數據降尺度。在直接降尺度過程中,不同源AGB產品之間存在的系統性誤差可能會導致降尺度驗證結果不準確,不能突出降尺度方法,而本文的研究目的是對森林地上生物量降尺度進行方法性探討,因此同時采用了利用升尺度模擬數據降尺度的思路。

利用統計回歸對森林地上生物量降尺度的關鍵在于建立適用于高分辨率森林地上生物量反演的模型。基于關系尺度不變的原理,將在低分辨率下線性/非線性回歸擬合的生物量反演模型進行殘差修復,并直接用于高分辨率下的生物量反演,從而利用降尺度手段獲取高分辨率森林地上生物量產品。

森林地上生物量的大小和分布格局受多種因素影響,直接因素包括林分類型、植被地域性和非地域性[24]、植被覆蓋密度等,間接因素包括地形、氣候要素以及人類活動等。因此選取能體現上述影響因素的多個地表參數作為自變量進行回歸。

本文分別采用了多元線性逐步回歸和非線性回歸方法構建降尺度模型。多元線性逐步回歸方法將多個地表參數同時用于構建方程,在構建過程中不斷引入和剔除自變量,在避免多重共線性的基礎上實現有效地表參數選取和最優回歸方程構建；非線性回歸方法僅采用和森林地上生物量相關性最高的單個地表參數作為自變量,建立冪函數模型。

兩種降尺度思路都需要采用聚合平均方式將30 m地表參數升尺度到1 km分辨率。在第二種降尺度思路中,還要將已有30 m分辨率森林地上生物量產品轉換到1 km。考慮到研究區中森林地塊和非森林地塊混合分布,1 km尺度下的生物量像元多為混合像元,純像元比例少。而本文建立的降尺度模型是針對森林的,由于非森林植被對光譜指數也有很大貢獻,如果直接建立低分辨率生物量和光譜指數間的關系而不進行像元篩選,則可能會存在低分辨率圖像上某像元森林生物量很低但光譜指數仍很高(非森林植被的貢獻所致)的現象,因此在建立關系時應盡量選擇比較純的森林像元。計算1 km網格空間上對應的所有30 m網格中森林像元個數占總像元個數的比例,并設定森林像元比例閾值,大于該閾值的1 km生物量像元作為因變量被選中,本研究設定森林像元比例閾值為0.6。而當降尺度模型在30 m尺度上估算森林生物量時,由于研究已經假定30 m尺度上都是純像元,因此該尺度上無需考慮森林像元比例。

多元線性逐步回歸降尺度和非線性回歸降尺度的思路基本相同,僅在步驟(1)建立回歸方程時有差異。多元線性逐步回歸降尺度模型建立時是將所有自變量輸入模型,在構建模型的過程中篩選自變量；而非線性回歸降尺度模型的自變量是在對每一個自變量和因變量進行相關分析后,選出的相關系數最高的單個自變量。降尺度步驟具體如下：

(1)將篩選后的1 km生物量和對應1 km地表參數值作為因變量和自變量數據,建立B1km和自變量之間的回歸方程：

(1)

式中,B1km、VI1km、ρ1km、S1km分別為1 km分辨率的森林地上生物量、植被指數、地表反射率和坡度。

(2)

(3)

式中,B30m、VI30m、ρ30m、S30m分別為30 m分辨率的森林地上生物量、植被指數、地表反射率和坡度。

(3)利用上述模型疊合高分辨率預估生物量和高分辨率殘差,獲得30 m分辨率生物量空間分布數據,從而實現整個降尺度過程。圖2為森林地上生物量降尺度流程。

圖2 降尺度流程圖Fig.2 The process of downscaling modelCMS:碳監測系統,Carbon Monitoring System

2.3 地表特征參數計算

選取了4種類型共20個地表特征參數作為降尺度回歸方程的自變量(表1),包括6個波段的地表反射率；多波段組合/植被指數；主成分變換分量和纓帽變換分量；地形因子。

表1 研究選取的地表特征參數Table 1 Land surface characteristic parameters used in the study

NDVI:歸一化差值植被指數,Normalized Difference Vegetation Index;SAVI:土壤調節植被指數,Soil-Adjusted Vegetation Index;RVI:比值植被指數,Ratio Vegetation Index;VI3:中紅外植被指數,Mid-infrared Vegetation Index;fc:植被覆蓋度,Fraction of forest coverage;STM:吸收反射與反照度因子,The ratio of absorption-band reflectance to albedo;TM452:第2/4/5波段組合,Combination of the second,fourth and fifth band;PC1:第一主成分,The first principal component;PC2:第二主成分,The second principal component;PC3:第三主成分,The third principal component;G:綠度,Greenness;B:亮度,Brightness;W:濕度,Wetness;S:坡度,Slope

表1中,r1、r2、r3、r4、r5、r7分別是TM影像的藍光、綠光、紅光、近紅外和兩個中紅外波段的地表反射率。SAVI計算公式中L為土壤調節系數,隨植被蓋度而變化,本文取常量0.5。

3 結果與討論

3.1 降尺度模型

通過多元線性逐步回歸和非線性回歸,得到森林地上生物量的估算方程(表2)。考慮到模型的普適性,建模過程中將兩個區域的自變量和因變量分別合并,建立了同時適用于兩個研究區的降尺度模型。進入擬合方程的各自變量系數都滿足P值小于0.05的條件,通過顯著性水平檢驗。

表2 森林地上生物量估算方程Table 2 The regression equations for forest aboveground biomass estimation

3.2 降尺度結果定性評價與分析

通過對兩個研究區降尺度30 m森林地上生物量數據和CMS 30 m森林地上生物量數據對比可以發現(圖3、圖4),利用模擬低分辨率數據降尺度后的30 m分辨率森林地上生物量空間分布和CMS森林地上生物量分布狀況大致相同,圖3b、3c和圖4b、4c中的生物量高值區和低值區與圖3a、圖4a吻合較好。而GEOCARBON AGB產品降尺度后的30 m分辨率森林地上生物量空間分布和CMS分布狀況相差較大,降尺度AGB數值分布集中,空間分布差異不明顯。此外,兩個研究區降尺度生物量分布圖中均存在較為明顯的斑塊結構,相較于CMS森林地上生物量數據而言其內部特征不夠平滑,這主要是由于在降尺度過程中未對殘差進行平滑處理,而是直接將其疊加到預估森林地上生物量分布圖上。

圖3 研究區1 CMS森林地上生物量分布和降尺度森林地上生物量分布比較Fig.3 Comparison between spatial distribution of CMS and downscaled forest aboveground biomass of study area 1(a)CMS結果;(b)模擬AGB多元線性逐步回歸降尺度結果;(c)模擬AGB非線性回歸降尺度結果;(d)GEOCARBON AGB多元線性逐步回歸降尺度結果;(e)GEOCARBON AGB非線性回歸降尺度結果

圖4 研究區2 CMS森林地上生物量分布和降尺度森林地上生物量分布比較Fig.4 Comparison between spatial distribution of CMS and downscaled forest aboveground biomass of study area 2 (a)CMS結果;(b)模擬AGB多元線性逐步回歸降尺度結果;(c)模擬AGB非線性回歸降尺度結果;(d)GEOCARBON AGB多元線性逐步回歸降尺度結果;(e)GEOCARBON AGB非線性回歸降尺度結果

3.3 降尺度結果定量評價與分析

對降尺度結果定量精度評價采用3個指標：相關系數r、均方根誤差RMSE、相對誤差δ。在降尺度30 m森林地上生物量分布圖和CMS 30 m森林地上生物量分布圖中均勻抽選對應像元點。篩選出其中的森林像元,得到散點圖圖5、圖6。

圖5 研究區1 CMS森林地上生物量和降尺度森林地上生物量散點圖 Fig.5 Scatter plots of CMS forest aboveground biomass versus downscaled forest aboveground biomass of study area 1 (a)模擬AGB多元線性逐步回歸降尺度;(b)模擬AGB非線性回歸降尺度;(c)GEOCARBON AGB多元線性逐步回歸降尺度;(d)GEOCARBON AGB非線性回歸降尺度

圖6 研究區2 CMS森林地上生物量和降尺度森林地上生物量散點圖 Fig.6 Scatter plots of CMS forest aboveground biomass versus downscaled forest aboveground biomass of study area 2 (a)模擬AGB多元線性逐步回歸降尺度;(b)模擬AGB非線性回歸降尺度;(c)GEOCARBON AGB多元線性逐步回歸降尺度;(d)GEOCARBON AGB非線性回歸降尺度

從表3、表4可以看出,研究區1和研究區2降尺度生物量的均值都高于CMS生物量均值。標準差反映森林地上生物量的空間變異性,生物量分布的標準差相對大小表現為：SDCMS>SD模擬AGB>SDGEOCARBONAGB,且非線性回歸降尺度生物量的標準差比多元線性逐步回歸的結果更低,說明在呈現生物量空間分布差異和細節變化方面,多元線性逐步回歸降尺度結果優于非線性回歸降尺度結果,這一結論也可從圖3b和3c、圖4b和4c的對比中看出,相較于多元線性逐步回歸降尺度生物量分布圖,非線性回歸降尺度生物量分布圖中的斑塊結構更加明顯,從而損失了較多細節。從散點圖(圖5、圖6)可看出,在生物量較小的情況下,降尺度生物量普遍高于CMS生物量；在生物量較大的情況下,降尺度生物量普遍低于CMS生物量。

非線性回歸降尺度的精度優于線性回歸降尺度,這可能是由于森林地上生物量和植被指數等自變量之間關系較為復雜,相較于線性模型,非線性模型能更好地描述它們之間的關系。此外,對比文中選擇的兩個研究區的降尺度狀況,研究區1比研究區2效果更好,這可能是由于研究區1的森林分布較研究區2更為集中,導致研究區1的森林純像元比例比研究區2更高,而本文采用的降尺度模型建立方法是利用森林像元比例閾值篩選建立降尺度模型的較純低分辨率像元,這種降尺度方法可能更適合于森林像元連續分布的研究區。

為了進一步探究兩個研究區降尺度結果存在差異的原因,本文基于景觀生態學原理,利用Fragstats軟件在斑塊類型水平上選擇部分景觀指數,分析兩個研究區森林地類之間的空間結構差異。斑塊類型層次景觀指數一般從斑塊的面積、斑塊形狀指數、斑塊密度等方面對斑塊不同類型的特征進行分析[28]。本文僅將兩個研究區分為森林和非森林兩種地類,方便對比兩個研究區森林地類分布的聚集程度。

表3 降尺度方法精度評價指標Table 3 Accuracy evaluation indexes of the downscaling method

表4 CMS森林地上生物量和降尺度森林地上生物量統計度量Table 4 Statistical measures of CMS forest aboveground biomass and downscaled forest aboveground biomass

研究選擇了斑塊密度、分離度指數和聚集指數共3個指標。斑塊密度反映了景觀破碎程度,值越大則破碎化程度越高。分離度指數是指景觀內某一景觀類型中斑塊之間的分離程度,值越小,表明景觀在區域內分布越集中,景觀較簡單。聚集指數基于同類型斑塊像元間公共邊界長度來計算,當某類型中所有像元間不存在公共邊界時,該類型的聚合程度最低[28]。

研究區1和研究區2森林地類的斑塊數分別為1362和1620。其他評價指標計算結果見表5。

表5 兩個研究區森林地類的景觀指數對比Table 5 Comparison between landscape indexes of forest in two study areas

從表5可看出,研究區1森林地類的斑塊密度和分離度指數都小于研究區2,且聚集指數大于研究區2,這說明研究區1的森林分布相較于研究區2更加聚集,破碎化程度更低,這與降尺度結果保持一致,也從側面說明本文的降尺度方法可能更適合于森林像元連續分布的研究區。

3.4 討論

本文的主要目的是探究森林地上生物量降尺度的方法,為避免不同源產品之間存在的系統性誤差影響降尺度效果,采用升尺度模擬低分辨率生物量數據和已有低分辨率生物量產品兩種數據源用于降尺度。但由于30 m和1 km之間尺度跨度大,且研究區地表覆蓋類型復雜,低分辨率生物量分布圖上多為混合像元,森林純像元少,且低分辨率分布圖的像元生物量偏低。盡管本文根據森林像元比例閾值,篩選出較純的低分辨率森林地上生物量像元用于擬合降尺度模型,但仍然不能完全避免混合像元對降尺度模型建立及精度評價的影響。混合像元引入對降尺度模型的影響主要表現在：對于由不同地表覆蓋類型組合(如農用地和森林、城市和森林)的低分辨率混合像元而言,假定其森林像元比例相同,則可能會出現森林地上生物量數值相同但對應光譜值(如某波段地表反射率)差異很大的情況,這對模型構建結果及降尺度結果有較大影響。因此,純像元選取是降尺度模型構建的關鍵環節。優化純像元選取策略及降尺度方法是下一步研究的重點。

從3.3節中的各種特征值指標、精度評價指標可以看出,本文中利用模擬低分辨率生物量數據建立的降尺度模型較好地實現了森林地上生物量從1 km空間分辨率到30 m空間分辨率的轉換,但GEOCARBON生物量產品降尺度的結果較差,原因主要有以下幾點：(1)數據源及反演方法不同,降尺度檢驗數據CMS產品與GEOCARBON產品之間存在系統偏差,圖7中1 km分辨率的GEOCARBON森林地上生物量分布和升尺度后CMS森林地上生物量分布在空間分布格局上本身存在較大差異,GEOCARBON生物量分布值普遍低于CMS升尺度結果；(2)GEOCARBON產品本身精度有限,產品自身精度會直接影響到降尺度結果。北半球溫帶地區GEOCARBON森林生物量是利用森林蓄積量(GSV)轉換估算的,GSV估算依賴于雷達后向散射測量,在密集的成熟林中信號易飽和導致GSV被低估,森林生物量被低估[23]。圖5、圖6散點圖中GEOCARBON森林地上生物量產品降尺度結果在高值區都存在明顯的低估現象。

圖7 模擬AGB和GEOCARBON AGB分布比較/(Mg/hm2)Fig.7 Comparison between spatial distribution of simulated and GEOCARBON forest aboveground biomass (a)研究區1,模擬AGB分布;(b)研究區1,GEOCARBON AGB分布;(c)研究區2,模擬AGB分布;(d)研究區2,GEOCARBON AGB分布

研究建立的降尺度模型,是基于森林地上生物量與光譜指數數據間關系尺度不變這一假設的。但事實上,由于低分辨率下建立的線性/非線性模型不一定完全適用于高分辨率尺度,這一假設也會給降尺度結果帶來誤差,從而無法準確判斷在轉換函數的尺度效應和降尺度模型本身的有效性之中,究竟哪一個因素對降尺度結果精度的影響更大。因此,接下來的研究可進一步探究降尺度模型隨尺度變化對森林生物量估算造成的誤差。

4 結論

本文以美國馬里蘭州的兩個區域為研究區,基于CMS 30 m分辨率和GEOCARBON 1 km分辨率森林地上生物量產品以及TM等數據源,通過升尺度模擬低分辨率生物量數據和直接使用低分辨率產品兩種方式,分別嘗試建立了多光譜地表參數和低分辨率森林地上生物量數據之間的統計關系,基于關系尺度不變的前提,將森林地上生物量數據的空間分辨率從1 km降尺度到30 m,并對降尺度結果進行精度評價和誤差來源分析。

模擬數據降尺度后的30 m分辨率森林地上生物量空間分布和CMS森林地上生物量分布狀況大致相同,RMSE=59.2—65.5 Mg/hm2,相關系數約為0.7；其降尺度結果優于GEOCARBON產品直接降尺度結果RMSE=75.3—79.9 Mg/hm2；相較于線性模型,非線性模型能更好地呈現森林地上生物量和地表參數間的關系；總體上,降尺度生物量呈現高值區低估,低值區高估的現象。