基于HJ-CCD数据的中高分辨率植被物候提取1)——以密云水库上游为例

王金英 李晓松 范文义

(中国科学院遥感与数字地球研究所(数字地球重点实验室)，北京，100101) (东北林业大学)

责任编辑:王广建。

植物物候反映了环境条件季节和年际变化的特点。传统的物候观测，记录植物在一年的生长过程中，随着气候的季节性变化而发生的(包括:发芽、展叶、开花、叶变色、落叶等)非常具体、直观的物候现象［1］，这些基于点的记录虽然具有准确详细的特点，但是覆盖面很小，不适于进行大区域气候变化研究。伴随着遥感技术而兴起的遥感物候学，将植被物候观测由植株个体提升到生态系统层面。利用遥感数据时间序列长、覆盖范围广的优势，有效地揭示了植被在时空上的改变与全球气候变化的联系［2-4］。各种观测表明，在全球气候变暖的背景下，植物生长季明显延长，引起植被生态系统的生产力、结构组成，以及土壤-植被-大气系统水、热、碳交换的改变。因此，植物物候的研究有利于提高气候-植被之间物质与能量交换的模拟精度，对于准确评估植被生产力与全球碳收支具有重要意义［5］。

受传感器时间分辨率的限制，大多数遥感物候关注的是大尺度的物候变化，监测数据源多为中低分辨率的卫星数据。有较长时间序列的AVHRR 数据有效记录了区域乃至全球物候多年的变化趋势，为全球气候变化的研究提供了有效途径［6-8］。MODIS 数据与AVHRR 数据相比，在空间分辨率上有较大提升，使遥感物候监测在相对复杂地表状态下能够实现。Ahl 等［9］利用地面观测数据，对不同的MODIS 植被指数产品用来检测植被物候的效果做了对比分析;Hudson 等［10］利用经过BRDF 校正后的MODIS 反射率数据，计算得到的植被指数，并结合MODIS 土地覆盖数据对植被春季返青期在海拔和纬度上的分异规律作了研究。然而，由于空间分辨率限制，这些物候监测结果与传统物候数据的可比性很小，所得的物候数据往往得不到有力的验证。而具有较高空间分辨率的Landsat 卫星影像由于重访周期长，即便是小区域物候监测的时间序列影像也难以得到［11-12］。综上所述，限于数据原因，遥感物候的监测尺度很难细化，这严重影响了其精度的验证，进而也阻碍了遥感物候在诸多领域应用价值的发挥。

中国于2008年9月6日发射的环境减灾小卫星，其A/B 星搭载的4 部CCD 传感器能够提供2 d重访的高时间分辨率的遥感监测数据，检验证明，CCD 传感器4 个波段中，近红外与红波段表现最好，与TM 传感器相应波段相比质量相当［13］。基于遥感的植被物候监测需要能够反映植被变化的连续时间序列数据，环境卫星的这些特点很好的满足了植被物候监测的需要。为此，本文选取密云上游地区，利用覆盖该区域的环境卫星数据，开展了30 m尺度植被物候提取的尝试，并与低分辨率提取物候及传统气象物候进行了对比，以期检验HJ-CCD 数据监测物候的价值与潜力。

1 研究区概况

密云水库上游地区包括密云、怀柔、延庆、沽源、赤城、崇礼、宣化、怀来、丰宁、兴隆、承德等12 个县的大部分或部分地区。该区位于欧亚大陆东部中纬度地区，属大陆性季风气候，四季分明。根据生态十年环境遥感监测土地覆盖分类系统［14-15］，建立该区域土地覆盖类型数据(见图1)，该区域自然植被覆盖面积比例为84.33%，耕地面积比例为13.93%。植被类型以落叶灌木林分布最广，面积比例为42.52%，其次为落叶阔叶林，面积比例为16.81%。在全球气候变化的大背景下，我国北方地区植被物候变化明显，主要表现为生长季的提前与延长［15］，这种变化会进一步引起区域生态系统生产力和物质循环的改变。密云水库作为京津地区重要的保障性水源，对其上游地区植被物候的变化监测具有重要意义。

图1 研究区域土地覆被类型图

2 研究方法

2.1 数据来源

遥感数据:本文选取2010年覆盖研究区域的所有环境星数据，共50 期，时间分布(见图2)。对获取的数据进行了辐射校正和几何精纠正等预处理。辐射校正利用中国资源卫星中心公布的环境减灾星(A/B 星)各载荷在轨绝对辐射定标系数，将图像DN 值转换为辐亮度;使用ERDAS 遥感图像处理软件，以该区域的TM 影像为参考图像，对其进行几何精纠正，纠正后的几何误差控制在1 个像元内。在对图像完成预处理之后，从中确定并截取研究区域子影像，应用6S 辐射传输模型建立查找表［16］，对影像进行大气校正。由于并不是每一景影像都能全覆盖研究区域，为了最大限度的保留变化信息，我们不采取时间域的融合，而是保留所有数据的原有属性。对于不能全覆盖研究区域的影像，将未覆盖区域标记为无效值，最终保证所有影像范围上的一致性，对于云覆盖的地方也标记为无效值。

气象数据:本研究气象数据来自中国区域地面气象要素数据集中的近地面气温要素，该数据集由中国科学院青藏高原研究所开发，主要用于近地面气象与环境要素的再分析［17］。数据集以国际上现有的Princeton 再分析资料、GLDAS 资料、GEWEX-SRB 辐射资料，以及TRMM 降水资料为背景场，融合中国气象局常规气象观测数据制作而成。其时间分辨率为3 h，水平空间分辨率0.1°。利用该数据集中的温度数据，对每天的数据求平均，合成2010年研究区域月平均温度数据，作为气象物候模型的基础数据。

图2 研究区域2010年环境星影像时间分布

2.2 时间序列归一化植被指数(NDVI)重构

利用经过预处理后的环境星影像，根据环境卫星CCD 数据的参数设置，利用公式INDV= (B4+B3)/(B4-B3)计算时间序列的归一化植被指数，其中:INDV表示时间序列归一化植被指数，B3表示近红外波段反射率，B4表示红波段反射率。NDVI 反映植被生长过程的时间序列曲线理论上应该是连续平滑的，但由于传感器本身性能、数据传输过程失误、太阳光照角度、观测视角、地物双向性反射以及云、大气气溶胶等观测条件因时间而异。此外，还可能有其它随机干扰因素，得到的观测值包含很多不可预测的噪音，NDVI 时间序列呈锯齿状的不规则波动，反映季节变化趋势不明显。为了获得高质量的连续植被指数数据，需要对数据集进行重建［18］。时间序列的一个本质特征就是观测数据在时间序列上的依存性，某一特定时刻的观测值与相邻观测值之间具有一定程度的依赖性。因此，在对参数时间变化规律进行归纳总结的基础上，可以建立描述遥感时间序列的数学模型，根据模型分析时序数据变化规律，以达到重建缺失数据，提高整体数据质量的目的。

本研究采用基于滤波函数平滑思想的S -G 滤波，应用于NDVI 时间序列数据重建，主要基于以下两个前提条件［19］:①序列的变化趋势应服从植被年际动态变化的渐进特征;②由于云和大气对NDVI的影响一般为负偏移，所以，多数噪声应低于NDVI序列数据的平均值。在利用S - G 滤波方法重建NDVI 序列数据的应用中，需要确定滑动窗口的宽度以及平滑多项式的阶数，以保证NDVI 数据拟合的准确性。由于不同区域数据量的不同，为获得理想的全年逐日时间序列NDVI 数据集，滑动窗口的大小需要依据原始数据覆盖的天数来确定，理论上原始数据量越大，窗口宽度也会越大。本研究多项式的阶数取2。通过S - G 滤波模拟出全年逐日NDVI 数据后，再通过局部循环，使重建后的时序数据更接近于原有NDVI 时序数据上的包络线。

2.3 物候参数提取

基于遥感的植被物候参数描述的往往是景观或者整个生态系统的季节性变化，与传统的特定物种的物候事件(发芽，开花等)无关［20］。遥感物候研究多关注于植物生长季的开始日期和成熟日期，由此得到植物生长季长度的变化。这些参数之所以受到广泛关注，原因在于:(1)受到遥感数据时空分辨率的限制，无法进行植被周期变化的详细描述;(2)鉴于遥感物候与传统物候的区别，生长季内的宏观尺度具体描述意义不大。因此，选择对生长季开始、结束日期以及长度的关注，对景观和生态系统层面的生物量和碳循环研究具有重要意义。根据本研究区域数据时空分辨率较高的特点，考虑到植物生长季长度由返青期和休眠期决定，本文选取植被的返青期、成熟期、衰老期、休眠期四个参数作为研究对象。

本研究物候参数的提取采用分段Logistic 拟合法［21］，对重建后的NDVI 曲线进行拟合，建立NDVI 的‘S’型变化曲线，公式为:y(t)=c/(1+ea+bt)+d。式中，t 为时间，即一年中的第几天;y(t)为时间t 的NDVI 值;c是拟合期间的NDVI 最大值;d 是NDVI 背景初值;a、b 为拟合参数。对研究区域采取逐像元计算各参数初值，自动灵活的确定拟合系数。拟合曲线的曲率计算公式为:k= -b2cz(1-z)2(1+z)3/((1+z)4+(bcz)2)3/2，式中:k为曲率，z=ea+bt。根据曲率极值确定的物候参数。在曲线上升阶段，曲率的最大和最小值分别为返青期和成熟期开始时间，在曲线下降阶段，曲率的最小和最大值分别为凋落期和休眠期的开始日期。

2.4 可靠性分析

遥感数据由于其独特的观测角度，其反映近地面状况的真实性往往需要通过一定的方式进行评估。由于观测角度的不同，遥感物候和地面观测物候可比性非常小，原因在于二者所获取的记录不仅有量上的差异，还有遥感像元的低分辨率造成其内部具有的异质性。基于气象数据的物候模型，在一定程度上解决了这个问题，因而被众多学者直接用来检验遥感物候的可靠性［22-23］，或者通过其建立遥感物候与传统物候的联系［24-25］。本文选取基于温度阈值法的气象物候模型，设定5 ℃为阈值［26］，求得该区域植被的气象物候。由于植被返青期的对区域气候变化的响应更敏感，在此我们仅提取位于每个气象物候像元中心点的遥感物候返青期，对二者进行相关分析，来评价遥感物候的可靠性。

3 结果与分析

3.1 NDVI 重建曲线

采用S-G 滤波对时间序列NDVI 进行重建，根据每个像元点的原始数据覆盖天数n，最终滤波窗口大小设定为n/4，采用二次多项式做窗口内拟合，得到研究区域逐日NDVI 数据。从不同植被类型的NDVI 曲线(见图3)可以看出，重建后的NDVI 曲线既消除了噪声，又反映了植被的季节变化特点;不同类型的曲线特点有明显区别，这与不同植被类型物候的差异性相契合。

图3 典型植被NDVI 曲线

3.2 环境星物候参数

基于环境星数据对研究区域内的物候参数反演结果(见图4)，利用计算出的返青期和休眠期之间的差，得到生长季长度。对不同植被类型各个物候参数进行统计，得到各类型的平均物候参数(见表2)。由表2可以看出，研究区域内各植被类型的物候参数不尽相同，其根本原因在于植被基因的差异［27］。落叶针叶林和草地两种植被类型春季返青期开始较早，平均日期为4月23 -24日，草地的生长期最长，平均为222 d，各植被类型开始凋落的时间为9月中下旬。草地的各物候期开始时的NDVI值和其他类型有较大差别，这些差别对于利用基于NDVI 值阈值法的遥感物候模型有一定的参考价值。

表2 不同植被类型的物候参数

3.3 可靠性验证

本文通过提取位于温度像元中心点的所有返青期的遥感物候值，建立遥感物候与气象物候之间的散点图，利用趋势线描述二者之间的相关性(见图5)。图左侧为该区域各尺度遥感物候与气象物候整体像元点之间的关系，右侧为该研究区域分布最广泛的落叶灌木林的遥感物候与气象物候的关系图。对于落叶灌木林像元点，这里参考土地覆盖数据，仅提取30 m 尺度下的纯净像元点作为验证数据。由图5可知，各个尺度的遥感物候与气象物候均呈现正相关关系，分辨率越高的遥感数据与气象数据的相关性越大;30 m 尺度下的遥感物候较之其它大尺度与气象物候呈现较强的相关性;各尺度下落叶灌木林的遥感物候与气象物候的相关性均没有整体的好，这说明了复杂地表覆盖下的像元异质性对遥感数据与地面数据的对接有较明显的影响。

图4 遥感物候结果

图5 遥感物候与气象物候关系图

4 结论与讨论

本研究以中国环境星CCD 数据为源数据，经过对图像严格的预处理，求得HJ -NDVI(尺度30 m)数据集，在最大限度挖掘时序数据的优势基础上，结合NDVI 重建技术，获取2010年研究区域NDVI日数据。采取分段Logistic 模型，计算出研究区域返青期、成熟期、衰落期和休眠期四个遥感物候参数。与传统气象物候模型得到的数据进行对比验证表明，基于环境星数据得到的30 m 尺度的地表物候参数与之具有较高的契合度。时间序列的环境星数据打破了遥感物候监测尺度难以细化的局限，弥补了遥感物候监测在海拔变化较大以及地表覆被复杂情况下的不足，同时，也为地面物候监测数据的获取提供了有效参考。

为了比较不同监测尺度对遥感物候结果的影响，本文利用同种方法计算了研究区域不同尺度下的植被返青期，通过局部细节和整体图来进行具体说明。图6a 来自利用环境数据做出的30 m 遥感物候;图6b 中数据是利用MODIS 植被指数产品MOD13Q1 中的250 m NDVI 数据集，基于本文所述方法得到;图6c 中数据来自于MODIS 全球土地覆盖变化产品MOD12Q2 中的500 m 物候参数集。从不同尺度物候结果图形可以看出，该地区的植被返青期在第120 天均有较集中分布，体现了落叶灌木为主要植被类型的区域特点。各数据返青期监测结果在整体上有较高的一致性;在500 m 尺度下，研究区域植被物候呈较严格正态分布;在250 m、30 m 尺度上分布较为分散，二者直方图有较明显的两个波峰，第一个波峰覆盖区域主要是东南部水库和河流周围的草地和灌木，第二个波峰为其余区域各种植被的物候组成。不同尺度下返青期标准差有较大差距，尺度越小其标准差越大，说明研究区域的地表异质性较高，小尺度下的物候数据才能反映这种特点。

图6 不同尺度返青期对比

为了进一步探索复杂地表状态下各尺度遥感物候数据反映地表真实环境的能力，我们选取了两个500 m 尺度的典型像元(见图7)，对不同尺度下它们的返青期做详细探讨。像元a 内部海拔高度变化范围为1 504～1 703 m，具有地势起伏较大的特点，植被类型为单一的落叶灌木林;像元b 内部包含三种植被类型，植被分布具有不均匀的特点，地势平缓。由图7可知，这两个像元所覆盖的地表变化在500 m 和250 m 尺度下几乎无反映，而30 m 尺度下的物候对这两种地表的复杂性均有较详细体现。然而30 m 尺度下物候与海拔高度的变化并不呈现明显的比例关系，植被类型的界限也不明显，这说明了自然状态下的植被物候与海拔植被类型之间没有严格的数学关系可描述，精细尺度的遥感物候对于自然植被覆盖的地表研究是非常必要的。

图7 不同尺度遥感物候对地面状况的记录

本研究区域环境星CCD 数据多集中在7—8月份和11—12月份，分别属于该区域夏季和冬季，春秋季节的数据略显不足，这也是导致遥感物候波动较大的原因之一。为了更好地发掘HJ -CCD 数据在植被物候监测方面的价值，今后需要寻求时间上均匀覆盖的时序数据。此外，由于基于气象数据的物候模型限制，我们在此仅比较了气象物候与遥感物候返青期之间的联系，为了更好的验证遥感物候的可靠性与实用性，未来的研究应寻求更完善的模型或者地面实测数据对遥感物候结果进行评价。

［1］ 竺可桢，宛敏渭.物候学［M］.北京:科学出版社，1980.

［2］ Goward S N，Markham B，Dye D G，et al.Normalized difference vegetation index measurements from the ad-vanced very high resolution radiometer［J］.Remote Sensing of Environment，1991，35(2/3):257 -277.

［3］ Fischer A.A model for the seasonal variations of vegetation indices in coarse resolution data and its inversion to extract crop parameters［J］.Remote Sensing of Environment，1994，48(2):220-230.

［4］ Verbesselt J，Hyndman R，Newnham G，et al.Detecting trend and seasonal changes in satellite image time se-ries［J］.Remote Sensing of Environment，2010，114(1):106 -115.

［5］ 李荣平，周广胜，张慧玲.植物物候研究进展［J］.应用生态学报，2006，17(3):541 -544.

［6］ Stockli R，Vidale P L.European plant phenology and climate as seen in a 20-year AVHRR land-surface parame-ter dataset［J］.International Journal of Remote Sensing，2004，25(17):3303 -3330.

［7］ Bacour C，Bréon F M，Maignan F.Normalization of the directional effects in NOAA-AVHRR reflectance measurements for an improved monitoring of vegetation cycles［J］.Remote Sensing of Environ-ment，2006，102(3/4):402 -413.

［8］ Maignan F，Bréon F M，Bacour C，et al.Interannual vegetation phenology estimates from global AVHRR meas-urements:Comparison with in situ data and applications［J］.Remote Sensing of Environment，2008，112(2):496 -505.

［9］ Ahl D E，Gower S T，Burrows S N，et al.Monitoring spring canopy phenology of a deciduous broadleaf forest using MODIS［J］.Remote Sensing of Environment，2006，104(1):88 -95.

［10］ Hudson Dunn A，De Beurs K M.Land surface phenology of north american mountain environments using moderate resolution imaging spectroradiometer data［J］.Remote Sensing of Environment，2011，115(5):1220 -1233.

［11］ Walker J J，De Beurs K M，Wynne R H，et al.Evaluation of landsat and MODIS data fusion products for analy-sis of dryland forest phenology［J］.Remote Sensing of Environment，2012，117:381 -393.

［12］ Melaas E K，Friedl M A，Zhu Z.Detecting interannual variation in deciduous broadleaf forest phenology using Landsat TM/ETM+ data［J］.Remote Sensing of Environment，2013，132:176-185.

［13］ 刘睿，孙九林，王卷乐，等.环境与灾害监测预报小卫星CCD数据质量评价［J］.地球科学进展，2011，26(9):971 -979.

［14］ 全国生态环境十年变化评估数据库系统.Wps1.gscloud.cn

［15］ 陈效逑，谭仲军，徐成新.利用植物物候和遥感资料确定中国北方的生长季节［J］.地学前缘，2000，S2:196.

［16］ 杜鑫，陈雪洋，蒙继华，等.基于6S 模型的环境星CCD 数据大气校正［J］.国土资源遥感，2010，22(2):22 -25，29.

［17］ Chen Y Y，Yang K，He J，et al.Improving land surface temperature modeling for dry land of China［J］.Journal of Geophysical Research，2011，116(D20):1 -15.

［18］ 顾娟，李新，黄春林.NDVI 时间序列数据集重建方法述评［J］.遥感技术与应用，2006，21(4):391 -395.

［19］ Chen J，Jönsson P，Tamura M，et al.A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golayfilter［J］.Remote Sensing of Environment，2004，91(3):332 -344.

［20］ 武永峰，李茂松，宋吉青.植物物候遥感监测研究进展［J］.气象与环境学报，2008，24(3):51 -58.

［21］ Zhang X Y，Friedl M A，Schaaf C B，et al.Monitoring vegetation phenology using MODIS［J］.Remote Sensing of Environment，2003，84(3):471 -475.

［22］ Justice C O，Townshend J G，Holben B N，et al.Analysis of the phenology of global vegetation using meteoro-logical satellite data［J］.International Journal of Remote Sensing，1985，6(8):1271-1318.

［23］ Reed B C，Schwartz M D，Xiao X.Remote sensing phenology［M］.New York:Phenology of ecosystem proc-esses，2009.

［24］ Schwartz M D，Reed B C，White M A.assessing satellite-derived start-of-season measures in the conterminous USA［J］.International Journal of Climatology，2002，22(14):1793 -1805.

［25］ Kathuroju N，White M A，Symanzik J，et al.On the use of the advanced very high resolution radiometer for development of prognostic land surface phenology models［J］.Ecological Modelling，2007，201(2):144 -156.

［26］ Shen M G，Tang Y H，Chen J，et al.Specification of thermal growing season in temperate China from 1960 to 2009［J］.Climatic Change，2012，114(3/4):783 -798.

［27］ Zhong L H，Hawkins T，Biging G，et al.A phenology-based approach to map crop types in the San Joaquin Val-ley，California［J］.International Journal of Remote Sensing，2011，32(22):7777 -7804.