基于多端元光谱混合分析方法的大兴安岭火后植被盖度恢复研究

陈 宝,刘志华,房 磊

1 中国科学院沈阳应用生态研究所,中国科学院森林生态与管理重点实验室，沈阳 110016 2 中国科学院大学,北京 100049

北方针叶林广泛分布于北半球中高纬度地区(50° N—70° N),约占陆地面积的17%和碳储量的35%[1]。火干扰是影响北方针叶林植被动态、结构和生态系统功能的重要的驱动因素。研究表明,气候变暖会增加北方针叶林中火干扰的频率、面积和强度[2],将会对北方针叶林生态系统中植物群落、植被覆盖度、生物量和景观格局造成更加显著的影响。因此,在当前全球变暖的背景下研究火干扰后植被恢复动态,对于理解干扰-生态系统-气候之间相互反馈机制,以及制定相应的火后森林恢复策略,具有十分重要的生态学意义[3]。

研究火烧迹地植被恢复主要有野外调查、模型和遥感等方法。基于野外调查的火烧迹地研究方法可准确得到生态系统结构和功能恢复状态,但费时、费力且具有时效性,难以在大时空尺度上连续监测火烧迹地植被恢复动态[4]。模型可以模拟大时空尺度上植被恢复,但模拟结果往往受模型的假设条件和输入参数影响,难以反映植被恢复的空间异质性[5]。遥感技术通过监测不同地物的光谱特征,可有效监测地表绿色植被,是监测火后植被恢复的有效技术手段[6]。随着遥感数据时间、空间和光谱分辨率的逐步提高,遥感技术已经成为研究火烧迹地植被恢复动态及其驱动力分析的主要手段之一[7]。遥感监测植被主要是由地表光谱反射特征决定的,在可见光波段内(380 nm—760 nm),叶绿素是植物光谱响应的主要影响因素,由于色素的吸收,绿色植被在可见光波段内具有较低的反射率(<20%),而在近红外波段(760 nm—800 nm),反射率急剧上升,形成“红边”现象。为了凸显绿色植被的光谱反射特征,研究人员构建了许多遥感植被指数。由于植被指数与地表植被的覆盖度和生物量具有很强的相关性,因此植被指数广泛用于监测、分析和绘制火后植被动态[8]。相关的植被指数有归一化植被指数(Normalized Difference Vegetation Index,NDVI)[9-10]、土壤调整植被指数(Soil Adjusted Vegetation Index：SAVI)、归一化短波红外指数(Normalized Difference Shortwave Infrared Index,NDSWIR)、增强型植被指数(Enhanced Vegetation Index,EVI)[11]等。

此外,混合像元在中、低分辨率遥感影像上普遍存在。火烧迹地通常是植被和基质的混合物组成,在中等分辨率的遥感影像(如Landsat)中是典型的混合像元[12]。混合像元分解(Spectral Mixture Analysis,SMA)是解决火烧迹地混合像元问题的一种有效途径。SMA是多种地表覆盖物光谱特征的加权混合值(以面积比例为权重)的理论假设,通过利用地表覆盖先验知识设定纯净光谱端元(即典型地物)并量化其光谱特征,结合多光谱遥感的多通道重复观测优势,对各光谱端元的覆盖度进行分解[13]。与传统的光谱指数方法相比,SMA可以探测一些低覆盖度部分,同时也可以减少土壤颜色和湿度的背景影响[14]。该方法产生的盖度图具有物理意义,相比植被光谱指数更具解释性[15]。

由Roberts等[16]于1998年提出的多端元光谱混合分析(Multiple Endmember Spectral Mixture Analysis,MESMA)方法考虑到了端元光谱的空间变异以及每个混合像元中地物数量的可变性。同时,MESMA在分解过程中将像元光谱视为端元光谱的线性组合,运用一套最优多端元模型分解得到每个端元的丰度图[16],多应用于植被种类制图、土地覆盖分类及信息提取等[17]。近期已有研究尝试运用MESMA技术监测火后植被恢复动态,如Fernandez等以地中海地区为研究区,以1999年至2011 年共13年的Landsat TM和Landsat ETM+ 序列数据为基础数据,使用MESMA分析的方法进行建模,监测火烧烈度以及火后恢复状况[18]。然而,MESMA能否适用于北方针叶林生态系统火后植被恢复监测尚无报道。

因此,本研究以我国大兴安岭寒温带针叶林为研究对象,以传统的NDVI为参考,评价MESMA在监测火后植被恢复中的有效性,拟回答以下三个问题：(1)与传统的NDVI相比,采用MESMA所获得的植被盖度的精度如何？(2)火烧烈度对MESMA和NDVI获得的植被盖度精度是否有影响？(3)最佳端元可否用于不同年份的遥感数据,实现时间序列分析？本研究结果将有助于评价MESMA在不同生态系统中的适用性,并加深对火干扰-生态系统-气候之间的相互反馈的理解。

1 研究区与数据源

1.1 研究区

图1 研究区位置图及2014年7月1日WorldView-2(2 m分辨率)合成真彩色火烧斑块图和非植被与植被示意图Fig.1 Study area,overlaid with true color composite from a July 1,2014 WorldView-2 image,typical non-vegetation and vegetation 黄色线表示2000年火烧边界

大兴安岭呼中自然保护区(122°42′—123°18′ E,51°17′—51°56′ N)属于中纬度大陆性寒温带,夏季短暂且温暖湿润,冬季寒冷干燥而漫长。年降水量约为460 mm(>80%的降雨发生在5—10月),研究区主要森林类型为寒温性针叶林,主要树种是兴安落叶松(Larixgmelinii)(>90%)和白桦(Betulaplatyphylla),其他针叶树(如樟子松：Pinussylvesrisvar.mongolica、云杉：Piceaasperata等)和阔叶树(如山杨：Populusdavidiana和甜杨Populussuaveolens)分布较少,通常在较干燥和排水良好的地点。林下和地面植被通常是由常绿灌木(如偃松：Pinuspumila等)、落叶灌木(柴桦：Betulafruticosa等)和一些草本植物(柳兰：Chamaenerionangustifolium、灰脉苔草：Carexappendiculata等)组成[19]。

本研究选取呼中自然保护区2000年6月17日发生的8700 hm2的雷击火烧区为研究区,区域内主要植被类型为成熟落叶松林。火烧区内分为高烈度(dNBR：differenced normalized burn ratio >1.099)火烧、中烈度火烧(1.099 >dNBR >0.484)和低烈度火烧(dNBR<0.484)。火烧后无人为干扰,植被以自然恢复为主,前期研究表明更新苗数量主要受火烧烈度和地形影响[20-22]。

1.2 数据源

本研究以30 m分辨率的Landsat ETM+ 为主要数据源,该数据来源于美国地质调查局(USGS)(http://glovis.usgs.gov/),选取该火烧迹地2010年6月22日和2014年6月1日无云的遥感影像数据,同时选取2014年7月1日2 m分辨率的WorldView-2高分影像,以及提取火烧烈度图作为最后精度验证的辅助数据。

2 研究方法

研究方法主要包括四个步骤：① Landsat ETM+影像预处理；② 多端元光谱混合分析(MESMA)程序：包括创建合适端元光谱库,选择最佳端元,对Landsat影像预处理之后的图像进行光谱混合分解,得到不同端元的盖度图像[23]；③ 基于传统的NDVI方法计算植被覆盖度；④ 精度验证：运用WorldView-2影像(2 m)计算Landsat像元尺度上(30 m)植被盖度作为验证数据,比较NDVI获得的植被盖图和MESMA获得的植被盖图。在此基础上,将选择的最佳端元运用到2010年Landsat影像上,对影像进行分解,比较2014和2010年Landsat获得的分数影像的分解精度,评估同一套端元光谱运用于不同年份的有效性。

2.1 图像预处理

USGS上获取的Landsat数据,已经经过辐射校正和几何校正,基本满足研究需要,可以直接使用。首先对选取的影像进行波段合成,再利用WorldView-2(2 m分辨率)高分影像和火烧斑块边界进行研究区裁剪,采用WGS84椭球体、UTM(51N)作为投影坐标系。

2.2 多端元光谱混合分析(MESMA)程序

多端元光谱混合分解是线性光谱混合分解的变形,端元选择决定了影像的分解精度[21]。光谱端元的选取通常有两种方法：从遥感图像中提取的“图像端元”(Image Endmember)和光谱库中得到的“参考端元”(Reference Endmember)[24-26]。图像端元可从遥感图像中运用一定的指标提取得到,具有两个优点：① 他们与图像数据的尺度相同；② 比较容易获得。参考端元是直接从已经存在的端元库中选取得到,具有以下优点：可以获得纯净的端元光谱；基于参考端元分解得到的丰度图更接近真实的土地覆盖情况。本研究以图像端元为主,选取最终端元光谱库。

根据研究区的地物特征,选取植被和非植被的感兴趣区域,共计1027个像元,生成1027条光谱曲线,从而构成初始端元光谱库。最终端元光谱库的选择主要有三种方法：基于计数的终端选择(Count-based Endmember Selection,CoB)、端元均方根误差(Endmember Average RMSE,EAR)、最小平均光谱角(Minimum Average Spectral Angle,MASA)[27-30]。使用CoB选择最佳端元时,选择拥有最高的类内(in_CoB)和最低的类外(out_CoB)的端元,即这个端元可以模拟同一类型内的所有光谱,但不模拟这一类之外的任何光谱。使用EAR时,选择一个类内最低RMSE(Root Mean Square Error)的最终端元。同时使用MASA时,选择具有最低平均光谱角的最终端元。

森林生态系统当中,植被、非植被是地表一定存在的两种地物类型,但由于地形、光照等因素的影响会造成一定程度的阴影,因此在遥感影像上表现为植被、非植被和阴影三种端元。本研究选取绿色植被和非植被端元光谱库,在每个Landsat像元上构建绿色植被(Green Vegetation,GV)-非植被(None Green Vegeattaion,NGV)-阴影(Shade)三端元模型,利用MESMA获得每个像元GV、NGV、Shade的丰度图像,并通过像元水平上RMSE值来评估模型的分解性能。本研究选取丰度图像的范围在0—1之间,最大允许阴影丰度值为0.8,RMSE值小于0.025。如果几个模型同时满足这些条件,选择最低RMSE的模型来获得MESMA分解的植被覆盖度(GVMESMA)。由于阴影的存在会给地物的分解结果带来一定的误差,所以我们在MESMA的结果上再进行阴影归一化。

NDVI计算获得的植被盖度：

FVC=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(1)

式中：FVC为植被盖度；NDVIsoil为裸地的归一化植被指数；NDVIveg为100%植被覆盖条件下的归一化植被指数。理论上100%植被盖度NDVIveg为100,裸地NDVIsoil为0。由于受到地形、气候、植被类型等多种因素的制约,在同一幅影像中NDVIveg和NDVIsoil都不能取固定值。以往研究中[31]将NDVI数值直方图与实地数据相结合,从NDVI统计直方图中读取累积频率为95%和5%时的NDVI值作为NDVIveg(0.847)和NDVIsoil(0.224)。运用(1)式计算Landsat影像水平上的植被盖度(GVNDVI)。

2.3 精度评价

本研究利用2014年WorldView-2数据计算的植被盖度(GVworldview)来验证MESMA获得的植被盖度(GVMESMA)和NDVI获得的植被盖度(GVNDVI)。GVworldview的计算过程如下：(1)计算2 m分辨率水平上WorldView-2数据的NDVI,绘制NDVI数值分布曲线,得到NDVIWorldView-2(NDVI >0.6)；(2)利用移动窗口法,计算30 m分辨率上平均植被盖度(GVWorldView-2)；(3)计算GVWorldView-2与GVMESMA和GVNDVI之间的相关关系,评价MESMA和NDVI获得的植被盖度。

2.4 2010年图像分解

为评估MESMA能否有效监测长时间序列的火后植被恢复,本研究将2014年影像中提取的最终端元光谱运用到2010年Landsat影像分解中,分解得到植被盖度影像GVMESMA2010。比较用同一套端元光谱分解得到的2014年和2010年植被盖度影像及两幅影像的RMSE。

3 结果与分析

图2显示了从每个光谱库中选取的一些最终端元的实例。利用CoB、EAR和MASA三个指标,每个指标选取最佳的2条植被光谱曲线和2条非植被光谱曲线,分别构成了6条植被光谱曲线和6条非植被光谱曲线的终端光谱库。运用这些终端光谱对影像进行光谱混合分解,得到不同端元的盖度图像。

图2 植被端元和非植被端元的光谱曲线Fig.2 Spectral signature for green vegetation and none green vegetation endmemberGV(green vegetation)表示植被；NGV(none green vegetation)表示非植被

图3显示了利用2014年Landsat ETM+影像进行多端元光谱混合分析后分解得到的阴影归一化植被和非植被盖度图。植被丰度图像显示了高植被盖度恢复区中存在明显的空间异质性(白色区域),可能与火烧烈度、生物物理条件和微气候有关。

图3 MESMA获得的2014年Landsat ETM+的植被(a)和非植被盖度图(b)Fig.3 Fraction image for green vegetation (a) and none green vegetation (b) on 2014 Landsat ETM+ image by MESMA method

图4为利用2014年高分影像为验证数据,分别对MESMA植被覆盖度和NDVI计算所得的植被覆盖度(2014年)进行精度验证的结果。结果表明MESMA获得的植被盖度和NDVI获得的植被盖度精度相近(MESMA:R2=0.691;NDVI:R2=0.700)。MESMA和NDVI获得的植被盖度在低盖度状态(<0.5)均有高估植被盖度的趋势(如数据点在1∶1线之上)。火烧烈度对两种方法的精度均有影响,中等烈度下的精度(MESMA:R2=0.719;NDVI:R2=0.695)都高于低(MESMA:R2=0.580;NDVI:R2=0.629)、高烈度下的精度(MESMA:R2=0.548;NDVI:R2=0.596)。中等烈度下,MESMA的精度比NDVI的精度略有提高；而在低、高等火烧烈度下,MESMA植被提取精度低于基于NDVI的植被盖度提取。总的来说,在大兴安岭地区火后植被恢复过程中,MESMA和NDVI方法获得的植被覆盖度差异不大(图4)。

图4 利用Worldview获得的植被盖度与MESMA和NDVI获得植被盖度在不同火烧烈度下精度评价Fig.4 Comparison of fractional image by WorldView-2 with those by MESMA and NDVI under different fire severity classes 黑线表示1∶1线,红线和蓝线分别表示NDVI和MESMA的验证结果

通过将从2014年Landsat影像中提取的最佳端元运用到2010年Landsat影像分解中,获得阴影归一化后的植被和非植被盖度图。从植被盖度图中可以清楚的看到植被和非植被区域。对比2014年(图3)和2010年植被盖度图像中的植被区域(白色区域),可清楚的观察到图3中的植被区域明显多于图5,而2014年的植被盖度图像中非植被区域(黑色区域)少于2010年。表明从2010年到2014年,绿色植被覆盖面积有所增加。比较用同一套光谱端元计算所产生的2014年与2010年的RMSE(2014std=0.003,2014mean=0.002,2010std=0.007,2010mean= 0.006),均值和标准差都在0—0.25之间,表明2014年的端元光谱运用到2010年影像上是可行的。因此,最佳端元可用于不同时相的遥感影像,可为后续的植被恢复时间序列研究提供依据(图3，图5)。

图5 利用2014年Landsat ETM+图像中获得的最佳端元运用到2010年Landsat图像得到的阴影归一化后(a)植被盖度图和(b)非植被盖度图像Fig.5 Fraction image for green vegetation (a) and none green vegetation (b) on 2010 Landsat ETM+ image by MESMA method using endmember from 2014 Landsat ETM+ image

4 讨论

遥感光谱指数和光谱混合分解是目前监测大时空尺度上火后植被恢复的最主要方法之一[32]。本研究考虑到不同端元光谱特征在空间上的异质性,利用6个绿色植被和6个非植被端元进行影像分解,获得植被盖度图像。利用MESMA和传统的NDVI方法计算植被盖度,并利用高分影像对两种方法进行精度对比,评价MESMA在火后植被盖度提取以及植被恢复研究中的适用性。理论上MESMA获得的植被覆盖度应该比NDVI获得的植被覆盖度精度更高[32],但本研究结果表明,MESMA方法所获得植被覆盖度图精度与基于NDVI的植被盖度精度相近。造成这种结果的主要原因可能有：(1)本研究植被盖度恢复在Landsat像元尺度上还没有达到饱和点,因此还在NDVI的有效监测范围内；(2)分解端元选择存在不合理性,造成分解精度没有预期效果。鉴于 NDVI已经广泛应用于北方针叶林生态系统火后植被恢复状况,并得到较好的研究结果。因此,就本研究结果来说,两种方法所获得的精度相近且获得较好的效果,表明在火烧后的15年内,MESMA和NDVI都是研究植被盖度恢复的有效手段。本研究选取分辨率为2 m的WorldView-2高分影像为验证数据,尽管与实地测量数据相比有一定的局限性,但相对于30 m分辨率的Landsat数据,高分影像植被盖度提取更为精确,在一定意义上可以作为验证数据。尽管如此,在今后的研究中,将进行野外考察,获取实地测量数据来验证我们的结果。

通过对MESMA结果的分析,我们认为最终端元的选取在MESMA方法中起着重要的作用,是准确估算丰度图像的关键步骤。但是,由于不同研究区土地覆盖类型多样,终端光谱库应当有足够数量的端元来代表每类地物的光谱变化,可能造成端元间的互相干扰增大,提高端元选择的难度[33]；此外,提取端元数量的增加会使端元分解模型数量增大,加重分解的计算负荷,从而降低解混的效率和精度。因此,一套好的最终端元一般包括线性独立性、光谱代表性和空间通用性等标准,这也是多端元光谱混合分析方法操作过程中的难点。为尽可能的避免这些问题的出现,本研究利用CoB、EAR和MASA三个指标来选择最佳端元,构建终端光谱库。

森林火烧烈度[34]是指森林火灾对森林生态系统的影响以及破坏程度。本研究发现,火烧烈度对两种方法的分解精度均有影响,中等烈度下的两种方法的精度均高于低、高烈度下的精度。以往的研究发现火烧烈度的大小直接影响森林生态系统中植被演替、森林火后动态、碳库的动态变化、枯枝落叶分解过程以及C/N循环等各种生态过程,间接影响地表植被的光谱特征[35],这就导致了MESMA和NDVI在不同火烧烈度状态下对于地表植被的光谱响应能力存在差异,在分解过程中才会呈现火烧烈度影响分解精度的情况。

本研究表明,MESMA与NDVI方法均适用于研究大兴安岭地区火后植被盖度提取。为了研究时间序列下的植被恢复,同时也为简化光谱的构建,本研究将同一套终端光谱应用于不同年份的影像分解,并得到这一套终端光谱适用于不同年份的一组影像,这可为今后时间序列植被恢复的研究提供依据。

尽管MESMA在研究森林火后植被恢复过程中存在一定的难度和问题,但从理论上说,MESMA是研究火后植被恢复的有效手段。加大对MESMA在森林生态系统中火后植被恢复的应用研究,可为以后的火后植被恢复研究提供强有力的技术手段。本研究主要用于植被与非植被的分解,并且也得到较为理想的结果。在今后的研究中,要进一步探究MESMA在森林生态系统中不同植被种类分解中应用的适用性,并扩展到其他的生态系统,为利用多光谱影像监测火后系统恢复提供新的视角。

5 结论

(1)MESMA方法获得的植被盖度(R2=0.691)与传统的NDVI获得的植被盖度(R2=0.700)精度无统计差异,中烈度下获得的植被覆盖精度高于低、高火烧烈度。

(2)为验证同一端元能否运用到不同时相的Landsat影像中,本研究将从2014年影像中获取的最佳端元运用到2010年影像中获得植被盖度图,结果表明2014年与2010年得到的RMSE(均方根误差)均值分别为0.0015和0.0065,说明最佳端元可用于不同时相的影像分解。

(3)本研究表明MESMA方法可有效监测北方针叶林中火后植被盖度恢复,并可运用于时间序列遥感影像监测植被恢复动态。