采用地表反射率模型的Landsat时序数据集重构

李凯， 陈芸芝， 汪小钦， 陈雪娇

(1. 福州大学 空间数据挖掘与信息共享教育部重点实验室， 福建 福州 350116；2. 福州大学 卫星空间信息技术综合应用国家地方联合工程研究中心， 福建 福州 350116)

高时空分辨率遥感数据在监测土地快速变化、作物生长,以及物候参数反演等方面具有重要的作用[1].Landsat系列卫星以较长的数据、较高的空间分辨率，以及较好的传感器性能在植被覆盖变化监测、植被生长评估等方面得到广泛应用[2].但在多云地区，Landsat系列卫星影像云覆盖率较高，且晴空条件下获取的影像数据时间间隔较长、季相不一致，这些都限制了影像的有效利用，使得这类传感器难以实现高频次、周期性的监测.通过重构Landsat时序数据集，插补缺失观测值，可加强对地表动态监测的能力.噪声的影响普遍存在于遥感图像中，通过一系列方法降低时序数据中的噪声水平，可为各种研究提供可靠的时序数据集[3].朱慧等[4]利用不同方法对重庆2010-2014年中分辨率成像光谱仪(MODIS)时序数据进行重建，认为WS(whittaker smoother)方法能较好地对不同土地覆盖类型进行重建.宋春桥等[5]对藏北地区各植被类型拟合效果较好的非对称高斯函数进行拟合，重构2007-2009年MODIS NDVI(归一化植被指数)时序数据集.刘亚南等[6]在传统单Logistic模型的基础上，采用参数构建模型对秦岭样区2001-2013年MODIS NDVI时序数据进行重建.李明等[7]利用同一类地物高质量像元均值代替噪声的S-G滤波算法，对江西省2001-2003年MODIS NDVI时间序列进行重建.殷悦等[8]通过6种方法对鄱阳湖平原地区2001-2013年的地球观测系统(SPOT)和MODIS时序数据进行重构，认为对于高植被覆盖区域S-G滤波算法的去噪效果最优.为满足大范围、高精度、快速变化的遥感监测，邬明权等[9]提出多源遥感数据时空融合技术，此类方法可以克服中等空间分辨率影像时间分辨率过低的局限[10].郭文静等[11]通过光谱差异改进ESTARFM(enhanced spatial and temporal adaptive reflectance fusion model)模型，对TM(thematic mapper)和AVHRR的NDVI数据进行融合，有效地实现若尔盖高原NDVI数据集的重构.孙锐等[12]采用分类拟合的方法，计算光谱距离权重改进STARFM(spatial and temporal adaptive reflectance fusion model)模型，融合中国环境一号A/B星CCD(HJ-1CCD)，重构山东省西北部NDVI时序数据集.张猛等[13]利用STARFM算法，对Landsat OLI(operational land imager)和MODIS数据融合构建空间分辨率为30 m、时间间隔为16 d的时序数据集，用于区域植被净初级生产力(NPP)的遥感估算.Meng等[14]基于STARFM算法，提出一种新的植被指数融合模型用于融合HJ-1CCD和MODIS数据，生成NDVI时序数据集.

目前，一系列基于降噪的时序数据重构方法已经取得了很好的成果，重构的遥感数据能准确地反映地表特征的变化[15]，但更适用于时间分辨率较高的时序数据集.时空融合算法利用时间信息和空间信息进行时序数据集的重构，但存在一定的局限性.首先，时空融合技术非常依赖较低空间分辨率影像的观测值；其次，该方法始终需要清晰无云的影像.因此，对于快速发生变化的区域，使用时空融合的方法进行变化监测效率较低.针对目前时间序列数据重构存在的局限性，Zhu等[16]提出地表反射率模型，该模型不依赖高时间分辨率影像，生成任意时间的清晰无云的Landsat影像.本文以福州市区某东部面积近300 km2的区域作为研究区，采用地表反射率模型进行Landsat时序数据集重构.

1 研究区概况与数据源

1.1 研究区概况

研究区位于福建省福州市区东部，北纬为26°0′55.81″～26°8′47.41″，东经为119°21′45.16″～119°33′39.22″，面积约为300 km2.研究区位置，如图1所示.

图1 研究区位置Fig.1 Position of study area

由图1可知：地形以山地丘陵为主，植被类型多样，森林覆盖率达到57.8%.研究区属典型的亚热带季风气候，全年云雨天气较多，质量好的遥感影像相对较少.

1.2 数据源及其预处理

基于研究目标，选取2013-2016年云覆盖90%以下的Landsat OLI影像，收集到的影像共44景.通过对原始影像进行大气校正，消除大气吸收和散射造成的辐射量误差.大部分影像受云、云阴影等影响严重，对于少量云覆盖的影像，采用Landsat OLI自带的第9波段云检测，设定阈值对云进行掩膜；对于含云量较高的影像，采用基于对象的Fmask算法进行云和云阴影的检测.由于去云、去云阴影后影像依然存在异常值，拟合去云、去云阴影后的时序影像.将拟合影像与Fmask算法检测结果进行对比，结果表明：对于没有检测出的异常值，其在蓝光波段的值明显高于模型的拟合值[17].

2 研究方法

时序地表反射率预测模型由谐波模型和一个趋势分量组成，用于预测每个波段的地表反射率，模型较为稳定、强健，不容易受噪声(云，云阴影等)影响.Zhu等[16,18]利用地表反射率模型,重构美国新英格兰地区1982-2011年Landsat影像，并利用拟合影像进行土地覆盖的连续变化监测；利用1982-2013年Landsat TM和ETM+影像，重构美国6个地区的Landsat时序数据集，并对重构影像进行评价.

地表反射率模型由4个系数组成:a0,i用于估算i波段整体的像元值或平均值；a1,i,b1,i用于模拟i波段由物候和太阳角度差异引起的年内变化；c1,i表示i波段的年际变化，用于估计i波段的长期趋势.影像去云、去云阴影和去异常值等预处理后，剩余每个像元的有效观测值个数达到12个，保证模型的精度和稳定.对于有效观测值个数小于12大于6的像元，采用地表反射率模型；对于有效观测值小于5个的像元，使用所有有效观测值的中位数来表示地表反射率.

采用最小二乘法(OLS)拟合，时序模型为

上式中:x为日期；T为一年的天数(T=365.25).

3 结果与分析

3.1 拟合结果分析

为研究区内每个Landsat像元生成时间序列模型，重构2013-2016年时序数据集，对2013-2016年所有可用的Landsat OLI观测值拟合NDVI时间序列.选取荷叶、阔叶树、柳树、马尾松4种不同类型的植被像元，拟合结果与原始NDVI对比，如图2所示.图2中：时间序列编号1～16代表2013年，17～39代表2014年，40～61代表2015年，62～84代表2016年.

由图2可知：原始序列中大部分点的NDVI值都与拟合序列较为接近；对于季节特征变化大的荷叶，NDVI值变化范围明显大于其他3种植被，对于四季常绿的阔叶树和马尾松NDVI值全年保持在较高的水平；随着有效像元数的增加，模型拟合效果越好，拟合结果与原始NDVI值越为接近.因此，拟合结果与原始有效像元保持较好的一致性，模型能够拟合不同植被稳定的生长变化过程.

(a) 荷叶 (b) 阔叶树

(c) 柳树 (d) 马尾松图2 拟合结果与原始NDVI对比Fig.2 Comparison of fitting results with original NDVI

3.2 实测光谱验证

2016年7月22日、2016年11月16日和2016年12月23日分别对研究区内荷叶、阔叶树、柳树、马尾松4种类型植被的光谱(图1)进行测量，利用实测数据，验证模型拟合结果.

拟合结果与实测NDVI对比，如图3所示.由图3可知：阔叶树、柳树的拟合结果变化趋势与实测NDVI基本保持一致，7-11月份NDVI值逐渐增大，11-12月份缓慢降低；马尾松、荷叶的NDVI值变化幅度小，NDVI保持在较高的水平，7-11月份的拟合结果变化幅度小于实测NDVI变化幅度，11-12月份两者变化趋势基本保持一致；样点的实测值高于拟合结果的NDVI值，主要是因为在像元尺度下，存在混合信息现象，导致拟合结果的NDVI值偏低，而实测值为一个纯净点且植被本身存在覆盖不均匀.因此，不同类型植被的拟合结果与实测值保持了较好的一致性.

(a) 荷叶 (b) 阔叶树

(c) 柳树 (d) 马尾松 图3 拟合结果与实测NDVI对比Fig.3 Comparison of fitting results with measured NDVI

3.3 原始影像验证

选择不同季节且影像质量较好的原始影像与拟合影像进行均方根误差分析，拟合结果与原始影像均方根误差，如表1所示.表1中：RMSE为均方根误差；Blue，Green，Red，NIR为影像波段的名称.

由表1可知：可见光波段的RMSE较小(以地表反射率为单位，约为0.01);近红外波段的RMSE值较大，这是因为研究区的大部分被植被覆盖，植被在近红外波段能产生较高的反射率;拟合结果与原始影像非常相近，差异很小，拟合精度较高.

分析2015年第270天的拟合结果与原始影像各波段相关性，拟合结果与原始影像散点图，如图4所示.图4中：x为拟合结果；y为原始影像;R为相关系数.由图4可知：拟合结果与原始影像的散点分布集中，4个波段都集中在直线y=x附近；各波段的相关系数均高于0.9，表明两者具有较强的相关性,整体拟合精度较高.

(a) Blue (b) Green

(c) Red (d) NIR图4 拟合结果与原始影像散点图Fig.4 Fitting results and original image scatter image

为了进一步验证模型的能力，拟合4个不同季节(2015年第270天；2016年第49天；2016年第97天；2016年第209天)影像，拟合结果与原始影像对比图，如图5所示.由图5可知：Landsat OLI波段5，4，3合成标准假彩色图像，植被在影像中大致呈红色；2015年第270天与2016年第209天的拟合影像红色区域面积明显高于2016年第49天与2016年第97天，夏秋两季植被覆盖明显高于春冬两季;拟合影像空间信息与地表覆盖的变化过程均与同一时间的原始影像保持较高的一致性，并且能够清晰地表达植被的物候差异；对于真实Landsat影像被云和云阴影覆盖的地方，拟合影像仍然可以提供清晰的观测.

(a) 2015年第270天(原始)(b) 2016年第49天(原始)(c) 2016年第97天(原始)(d) 2016年第209天(原始)

(e) 2015年第270天(拟合)(f) 2016年第49天(拟合)(g) 2016年第97天(拟合)(h) 2016年第209天(拟合)图5 拟合结果与原始影像对比图Fig.5 Comparison of fitting results with original image

拟合结果与原始影像NDVI差值图，如图6所示.由图6可知：2015年第270天、2016年第49天、2016年第209天这3个时间点的大部分区域NDVI差异较小，极小部分差异较大的区域主要集中在地表覆盖变化幅度大的区域；图5(c)影像存在云的影响，原始影像红光波段的反射率整体偏高，导致拟合结果与原始影像NDVI的差值偏大；拟合结果与原始影像的差异较小，说明模型的拟合效果较好.

(a) 2015年第270天 (b) 2016年第49天

(c) 2016年第97天 (d) 2016年第209天图6 拟合结果与原始影像差值图Fig.6 Difference of fitting results with original image

4 结论

采用地表反射率模型，基于2013-2016年云覆盖90%以下的所有Landsat OLI影像，构建Landsat遥感时序数据集.通过地表反射率模型获得的拟合结果与真实影像之间整体具有较高的相关性，各波段的相关系数均高于0.9，均方根误差较低，可见光波段在0.01左右，近红外波段在0.02左右.对不同类型植被均能取得较好的拟合效果，不同类型植被拟合NDVI与实测NDVI变化趋势基本保持一致.拟合结果能清晰表现出植被的物候差异，与原始影像之间的差异较小，其空间信息及地表覆盖的变化过程也与原始影像保持较高的一致性.通过地表反射率模型提供的密集清晰无云影像可用于植被动态变化监测研究，进一步提高植被变化监测的准确性和时效性.

该模型需要大量清晰的Landsat OLI影像进行精确的时间序列模型预测，对于长期多云和雪覆盖的区域，可能没有足够的清晰影像满足模型的要求，无法进行精确的时间序列模型预测，下一步将针对这方面进行改进，提出更加全面的时序重建方法.