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

刘沁茹,孙 睿,*

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；相较于线性模型,非线性模型能更好地呈现森林地上生物量和地表参数间的关系；总体上,降尺度生物量呈现高值区低估,低值区高估的现象。