基于Sentinel-2数据的典型草原光合/非光合植被覆盖度估算

柴国奇，王静璞，邹学勇，王光镇，韩 柳，王周龙

(1.鲁东大学资源与环境工程学院，山东 烟台 264025； 2.北京师范大学地表过程与资源生态国家重点实验室，北京 100875)

草地作为全球分布最广的植被类型之一，在陆地生态系统碳循环中起着重要的作用[1]。草地植被从功能角度可以分为光合植被(photosynthetic vegetation，PV)和非光合植被(non-photosynthetic vegetation,NPV)，其中PV主要指进行光合作用的绿叶，NPV主要指枯草、枯叶、枯枝、枯干、凋落物等[2]。PV无疑是草地植被的重要组成部分，但它并不是唯一的组成部分[3-4]，NPV在草原生态系统中发挥着重要的作用，是评价地表植被覆盖状况的重要指标[5-6]。PV和NPV不仅能够有效减少径流和蒸发，增加土壤水分渗透和有机质含量，提高土壤湿度和温度，改善土壤质量[7-10]，而且还影响着生态系统中C储存、CO2交换量以及营养物质的流动与循环[10]。另外，NPV在火灾风险和频率、风和水侵蚀等方面也起着关键作用[11]。然而，近年来人们对草地的过度使用以及雪灾、旱灾的频繁发生造成草原大面积退化、沙化和盐碱化，水土流失和风沙危害日益严重[12]。因此，及时、准确地掌握草原光合植被覆盖度(fractional cover of photosynthetic vegetation,fPV)和非光合植被覆盖度(fractional cover of non-photosynthetic vegetation,fNPV)信息对草原的合理利用以及草原生态系统的监测和保护具有重要意义。

遥感技术的发展为快速、准确地获取大范围fPV和fNPV信息提供了新的技术手段。近几十年来，国内外利用遥感技术对fPV的估算研究已取得较大进展，但对于fNPV的估算研究较少，同步获取fPV和fNPV信息的研究更少[13]。研究表明光谱混合分析(spectral mixture analysis,SMA)是利用多光谱影像估算fNPV的一种有效方法[11,14]。然而，每个像元BS背景的光谱特征是影响fNPV估算效果的关键因素，由于一些地区BS光谱变化较大，导致了SMA估算fNPV变得较为困难[3]。多端元光谱混合分析(multiple endmember spectral mixture analysis，MESMA)是对SMA的一种改进[15]，可较好地解决包括BS在内的地物类内光谱异质性问题。目前，MESMA在fPV和fNPV的估算中已得到广泛应用，并取得了较好的估算效果[3,16-17]。如Meyer和Okin[18]基于MODIS数据建立PV、NPV和BS光谱库，利用MESMA估算博茨瓦纳西部地区的fPV和fNPV，R2分别为0.62和0.55，RMSE分别为0.09和0.16。Numata等[19]利用Landsat TM数据，较好地估算了亚马逊地区6个牧场的fPV和fNPV，提高了监测牧场退化的能力。以往研究多基于较低分辨率的数据(如MODIS)，空间异质性明显，估算精度具有不确定性；分辨率相对高的数据(如Landsat)重访周期较长，很难获取较密集的时间序列的影像数据。因此，利用时间和空间分辨率较高的多光谱影像数据将成为估算fPV和fNPV的研究重点。

Sentinel-2环境监测卫星目前已发射Sentinel-2A和Sentinel-2B两颗，其搭载的多光谱成像仪(MSI)能够提供最高10 m的空间分辨率及5 d重返周期的多光谱数据。但由于该卫星是新发射卫星，目前将Sentinel-2影像数据应用于草原地区开展fPV、fNPV监测的研究较少。为此，本研究以锡林郭勒典型草原为研究区，以Sentinel-2A多光谱影像为数据源，利用MESMA估算锡林郭勒典型草原的fPV和fNPV，并对其估算的结果进行精度评价，在此基础上，分析典型草原fPV和fNPV的时空变化。

1 数据收集与处理

1.1 研究区概况

研究区位于内蒙古自治区锡林郭勒盟中北部(114°46′-118°25′ E、43°15′-46°02′ N),海拔800～1 400 m(图1)。地处中温带半干旱大陆性季风气候区，四季分明，气温偏低，年均温0～3 ℃，无霜期90～120 d；雨季短促，年均降水量为295 mm，由东南向西北递减，主要集中在夏秋两季。研究区地貌以波状高平原为主，由东南向西北方向倾斜，东、南部多低山丘陵；西、北部地形平坦，为高原草地；土壤以风沙土为主，部分地区有栗钙土、粽钙土，植被群落中主要有羊草(Leymuschinensis)、大针茅(Stipagrandis)、克氏针茅(S.krylovii)等。

1.2 野外数据获取

2017年9月28日至10月3日在锡林浩特市、东乌珠穆沁旗和西乌珠穆沁旗开展了野外地面观测试验，采集的数据为实测样地PV、NPV和BS盖度数据。本研究共获得51个样地数据(图1)，具体采集步骤为：首先，在研究区相对平坦均质的地区选择一个20 m×20 m(对应Sentinel-2A影像像元)大小的样地，在样地中心通过手持GPS记录相应的地理位置以及样地植被状况等信息；其后，在每个样地内布设3～5个0.5 m×0.5 m的样方，取样地内所有样方盖度的平均值作为该样地最终的盖度。样方盖度测算采用照相法，利用数码相机于样方中心点正上方1.5 m处垂直拍摄2张照片，并记录每个样方对应的照片编号。

图1 研究区及采样点位置Fig. 1 Study areas and sampling location

在测算样方盖度前，首先通过ENVI 5.3软件切除照片中非样方部分，剩余样方部分用于进行植被覆盖度的测算，根据照片的RGB值分出PV-NPV-BS三端元并建立训练样本。然后利用监督分类中的最大似然分类(Maximum likelihood classification,MLC)算法对照片进行分类处理，并结合目视解译验证分类结果的精度。最后统计出各样方的fPV、fNPV、fBS，取2次结果的平均值作为样方的最终盖度[20]。

1.3 遥感影像获取与预处理

Sentinel-2A卫星属于欧洲“哥白尼计划”中Sentinel-2系列卫星的第一颗光学遥感卫星，已于2015年6月23日发射升空。Sentinel-2系列卫星能够用于时间、空间分辨率较高的全球陆地观测，可提供植被、土壤和水覆盖等图像，还可用于重大自然灾害的监测。Sentinel-2A卫星携载MSI多光谱成像仪，可覆盖13个光谱波段，观测幅宽达到290 km，单颗卫星重访周期为10 d，具有10、20、60 m 3种不同的空间分辨率，部分参数信息如表1所列。

表1 Sentinel-2A卫星数据部分参数信息Table 1 The characteristic of Sentinel-2A

Sentinel-2A影像已经经过辐射定标和几何校正处理，投影坐标系为UTM/WGS-84，是Level-1C大气顶表观反射率数据产品，包括陆地、水、云掩模等数据。通过欧空局(ESA)的数据共享网站(https://scihub.copernicus.eu/dhus/#/home)获取到2017年4月8日、5月18日、7月17日、8月26日和10月5日共5期Sentinel-2A影像。

影像预处理主要通过ESA提供的SNAP 6.0(sentinel application platform)和ENVI 5.3软件平台完成。首先，利用SNAP 6.0的Sen2Cor(Sentinel to correction)插件对所有L1C数据进行大气校正处理，输出为20 m分辨率的9个波段(Band 2、Band 3、Band 4、Band 5、Band 6、Band 7、Band 8a、Band 11和Band 12) L2A地表反射率数据；然后，利用ENVI 5.3将9个波段数据合成多光谱影像，并进行拼接裁剪处理；最后，用阈值法将影像上的水体掩膜。

2 研究方法

2.1 MESMA模型

MESMA是在线性光谱混合分析(liner spectral mixture analysis, LSMA)的基础上发展而来的一种混合像元分解方法，它针对每一个像元建立一系列LSMA模型，然后根据最小的RMSEs约束为每个像元挑选出最佳模型进行分解，可有效解决同物异谱问题[15-16]。每个LSMA模型是基于这样一种假设：混合像元的光谱特征是由像元内不同地物的光谱特征线性组合而成[21]。可用以下公式表达：

(1)

(2)

评价模型用均方根误差(root mean square error,RMSEs)表示：

(3)

式中：M为遥感影像波段数。

2.2 PV、NPV及BS端元光谱库建立

2.2.1端元类型确定 对于线性光谱混合分析，端元的选取是混合像元分解的关键[22]。由于波段之间相关性的影响，选取的端元过多会增加模型对端元选择的敏感性[23]；而如果端元数量太少，则无法准确解释更多的光谱变化[22]。此外，不同地区的端元选取方式也存在一定差异，如城市地区常用植被-不透水层-土壤三端元模型[24]；而非城市地区，多采用植被-土壤-阴影(或干枯植被)模型[25]。由于本研究区主要分布在相对平坦均质的草原地区，考虑到模型的解混精度以及地表覆盖特点，最终选取PV、NPV及BS 3个端元。

2.2.2最优端元选取 MESMA需要PV、NPV及BS每类地物的端元光谱(从影像本身、地面实测或实验室获得)建立端元光谱库[26]。在大气校正后的Sentinel-2A影像中，影像光谱与地面实测光谱存在一定差异，为消除这一差异影响，本研究使用纯净像元指数法(pixel purity index,PPI)从影像本身直接获取某地物纯净像元光谱作为端元光谱，以保证与影像的尺度一致[25,27]。利用端元均方根误差(endmember average root mean square error，EAR)[23]、最小平均光谱角(minimum average spectral angle,MASA)[28]选取光谱库中的最优端元光谱。具体步骤如下，首先对Sentinel-2A影像进行MNF(minimum noise fraction)变换，取变换后的前6个波段计算PPI指数；然后人机交互式选择PV、NPV和BS每类地物中较为纯净的端元光谱创建光谱库，利用“VIPER Tools”计算光谱库中每条光谱的EAR、MASA；最后将选取最优端元光谱分类并建立PV、NPV和BS光谱库。

2.3 MESMA的混合像元分解

MESMA模型结构较为简单，物理意义明确、计算方便，本研究以建立的PV、NPV和BS光谱库为基础，对Sentinel-2A影像进行分解。由于地物阴影在遥感成像过程中的存在[25]，MESMA将阴影作为一个端元考虑，不仅提高了fPV和fNPV估算精度，也更符合实际情况[3]。因此，本研究主要运行PV、NPV、BS和阴影四端元模型，采用以下约束条件：端元的比例约束为-0.05～1.05，最大阴影比例为0.8，最大RMSEs为0.025。选择产生RMSEs最低的模型计算的端元丰度，然后将阴影归一化处理，最后得到每个像元PV、NPV、BS端元的丰度，即为fPV、fNPV和fBS。

2.4 精度检验

采用均方根误差(RMSE)、平均误差(mean error, ME)和决定系数(R2)来评价模型估算精度。

(4)

(5)

(6)

3 结果与分析

3.1 最优端元选取结果

在“VIPER Tools”计算结果中，选择EAR、MASA 值较低的端元为最优端元，最终选取了19条最优端元光谱，分别是6个PV端元、8个NPV端元和5个BS端元(图2)。

PV、NPV、BS端元内部光谱差异相对较大(图2)，如果不考虑地物类内光谱异质性的存在，不可避免的会对混合光谱分解结果的合理性产生影响。从平均端元光谱图中看到PV在VIS-NIR波段(400～1 100 nm)反射率差异明显，因此，PV相较与NPV和BS相对容易区分。然而，NPV与BS在VIS-NIR波段不仅没有特别的特征，而且光谱曲线非常相似，其区分则相对较为困难。尽管如此，NPV和BS光谱还是存在一些区别，BS在1 610 nm(Sentinel-2A数据Band11)和2 190 nm(Sentinel-2A数据Band12)处具有相似的反射率，而NPV在2 190 nm处比在1 610 nm处的反射率低得多。因此，利用SWIR波段的光谱特征来区分NPV成为一种可能。

图2 PV、NPV、BS端元的光谱曲线Fig. 2 Reflectance spectral of PV, NPV and BS endmember

PV、NPV、BS分别为光合植被、非光合植被和裸土。

PV, NPV and BS were photosynthetic vegetation, non-photosynthetic vegetation, and bare soil, respectively.

3.2 fPV和fNPV的估算结果

根据选取的PV、NPV、BS最优端元光谱，利用MESMA对10月5日的Sentinel-2A影像进行分解，估算出fPV、fNPV。研究区内超过97.7%的像元(掩膜水体后)能成功分解，模型分解精度RMSEs的平均值为0.004 6，误差相对较小。从绘制的fPV、fNPV和fBS的空间分布图(图3)中可以清晰地看出，总体的分布特征与本次野外调查基本一致，从西到东fPV、fNPV在逐渐增大，整体呈现东高西低的趋势。由于数据采集时间正处于牧草的收割期(物候条件为黄枯期)，研究区主要以NPV和BS为主，fNPV高的地区主要分布在降水较多的东部，fBS高的地区主要分布在降水相对较少的西部，而PV则多分布在森林和河流湖泊周边地区。未分解的像元约占2.3%，主要分布在东南部的丘陵地区，这是因为山区阴面植被受山体遮挡，阴影现象严重。还有少量分布在城市地区，城市景观复杂多样，像元混合严重，因此对城市区分解时可能会产生比较大的误差。

图3 fPV、fNPV和fBS的空间分布图Fig. 3 Map of fPV, fNPV and fBS distribution

3.3 估算精度检验

为了验证MESMA方法对fPV、fNPV和fBS的估算效果，本研究采用研究区51个样地的实测数据进行精度评价。其估算的fPV、fNPV与实测fPV、fNPV之间具有较好的相关性，R2分别为0.72和0.61(图4)。对于fPV来说，估算的RMSE为5.52%，ME为1.49%，散点分布在参考线的Y=X的左上方较多，出现了一定程度的高估；对于fNPV来说，估算的RMSE为9.37%，散点基本分布在参考线的Y=X的周围，但ME值为正，说明也出现了一定程度的高估。估算结果的高相关性，表明利用MESMA可以较好地估算fPV、fNPV；而估算的fNPV精度低于fPV表明，PV相较于NPV较容易区分，由于NPV与BS端元光谱特征相似，导致NPV与BS误分的可能性相对较大。

3.4 fPV和fNPV的季节性变化

本研究选择研究区内受人为以及牲畜活动影响较小的草地(图5)作为感兴趣区，利用MESMA估算5期Sentinel-2A影像的fPV、fNPV，并绘制感兴趣区fPV、fNPV和fBS的空间分布图(图5)。感兴趣区的fPV和fNPV以均值代表，通过分析fPV和fNPV(图6)发现，MESMA估算的fPV和fNPV的季节性变化符合草原植被物候发育特征。4月初分布大量上一年遗留的NPV，几乎没有PV，总体呈现荒芜状态，fNPV约为56%，fPV仅为2%；4月中旬以后，牧草返青期到来，20 d后牧草进入快速生长期，fPV初步增加，fNPV相应减少；5月18日时，fPV已经增加到8.5%，fNPV减少到27.75%；在7月中旬左右，大多数牧草生长达到最大值，7月17日时，牧草已经比较旺盛，fPV(约62.18%)明显多于fNPV；8月份牧草相继进入成熟期，在8月26日时，NPV已基本消失，fNPV仅为2.5%，fPV增加到78.89%；9月份左右牧草陆续开始枯黄，到10月5日时，大部分牧草已经枯黄，草地分布大量NPV，基本已恢复到枯草状态，fNPV高达60.67%，fPV约为2.5%。

图4 fPV、fNPV和 fBS估测值与实测值散点分布图Fig. 4 Scatterplot of estimated values and measured values of fPV, fNPV and fBS

4 讨论

本研究结果表明，以Sentinel-2A多光谱数据为数据源，在PV、NPV和BS三组分共存的情况下，利用MESMA方法能够同时估算出fPV和fNPV，并能取得较好的估算效果，这与Meyer和Okin[18]、Numata等[19]使用其他多光谱数据的研究结果相一致。利用MESMA估算fPV和fNPV有两大优点。一是，MESMA计算所有波段的反射率信息，减少了BS背景的影响[17]；二是，MESMA较好地解决了混合像元同物异谱的问题，能够较准确地估算fPV和fNPV。然而，Okin等[3]发现，在一些特殊区域MESMA并不能取得较好的效果，由于NPV与BS端元光谱特征相似，存在MESMA为了得到光谱拟合更小的RMSEs而造成实际的PV、NPV和BS端元组合被模型中的另一种组合所取代的可能，使得估算结果出现偏差[4,26]。因此，如何确定最合适的PV、NPV和BS端元或者如何优化MESMA来提高fPV和fNPV估算精度需要进一步的分析。

图5 不同季节fPV、fNPV和 fBS的空间分布图Fig. 5 Map of fPV, fNPV and fBS distribution in different seasons

图6 感兴趣区fPV、fNPV的季节性变化图Fig. 6 Seasonal variation of fPV and fNPV in the studied region

纯净端元的提取对于混合像元的分解精度起着至关重要的作用[4,12,29]。本研究利用MNF变换、PPI计算能够提取影像中比较纯净的像元，并且可以确定纯净像元的空间位置，但在实际的遥感影像上很难找到较纯净的像元，在纯净像元较少时，容易产生较大的误差。如果能将纯净的实测光谱端元与影像端元进行结合，充分利用实测光谱端元调整和校正影像端元，兼顾两种端元的优势来构建端元光谱库，可进一步提高MESMA的精度[30]。

本研究探讨了MESMA方法应用到Sentinel-2A多光谱遥感数据估算典型草原地区fPV和fNPV的潜力，目前在国内外采用该数据估算草原fPV和fNPV还鲜有研究。Sentinel-2A数据相较于MODIS、Landsat数据拥有更高的空间分辨率，在一定程度上会降低混合像元的空间异质性以及尺度效应对估算结果的影响，从而进一步提高估算精度。本研究通过分析锡林郭勒草原fPV和fNPV季节性变化发现，其季节性变化同草原植被物候发育特征相吻合，这与李涛等[12]、王光镇等[29]分别采用NDVI-CAI三元线性混合模型和NDVI-DFI像元三分模型估算fPV和fNPV的研究结果一致。因此，借助于Sentinel-2数据较高空间分辨率和时间分辨率(5 d)的优势可以进一步应用于典型草原fPV和fNPV的时空动态变化分析。

值得注意的是，由于野外实测数据与影像数据存在时间上不同步、样方盖度计算具有主观性的问题，这些数据误差也在一定程度上影响fPV和fNPV的估算精度。另外，本研究仅对单期遥感影像的估算结果进行了精度检验，样地数量相对较少且缺乏系统的地面真实性检验工作，模型的普适性并未得到充分的验证。因此，在后续的研究中将选取固定观测站点进行植被覆盖度连续观测，获取更多不同植被生长状况下的地面采样数据，以更好地评价MESMA估算典型草原fPV和fNPV的精度。

5 结论

本研究基于Sentinel-2A多光谱数据，利用MESMA方法估算锡林郭勒典型草原的fPV和fNPV，并通过实地观测数据进行评价，主要得到以下结论：

1)对于锡林郭勒典型草原来说，PV、NPV和BS的光谱库端元存在明显的光谱异质性特征，确定最优的端元组合对于混合像元的分解至关重要，是成功估算fPV和fNPV的关键。

2)在PV、NPV、BS 3种组分都存在的情况下，MESMA可有效估算典型草原地区的fPV和fNPV，fPV估算的RMSE为5.52%(R2=0.72)，fNPV估算的RMSE为9.37%(R2=0.61)。

3)以Sentinel-2A MSI多光谱影像为数据源，利用MESMA估算的fPV和fNPV的季节性变化符合草原植被物候发育特征，可应用于典型草原fPV和fNPV的时空动态变化快速监测。