GEE时序——利用sentinel-2（哨兵-2）数据进行地表物候学分析（时间序列平滑法估算和非平滑算法代码）

简介

哨兵-2A/B 串联卫星的空间分辨率高、重访时间长，有可能改进对陆地表面物候的检索。不过，生物群落和区域特征在很大程度上限制了陆表物候学算法的设计。在北极地区，这种生物群落特有的特征包括长期积雪、持续云层覆盖和生长季节短暂。在此，我们评估了哨兵-2 获取北极高分辨率 LSP 地图的可行性。我们通过在谷歌地球引擎（GEE）中简单实施阈值法，提取了 2019 年和 2020 年的季节开始和结束时间（分别为 SoS 和 EoS）。我们发现哨兵-2 和 PhenoCam 的指标之间具有高度的相似性；哨兵-2 增强植被指数（EVI）的结果最佳（SoS 的均方根误差（RMSE）和平均误差（ME）分别为 3.0 d 和 -0.3 d，EoS 的均方根误差（RMSE）和平均误差（ME）分别为 6.5 d 和 -3.8 d），尽管其他植被指数的表现类似。在同质地貌中，哨兵-2 EVI 的物候图与中分辨率成像分光仪（MODIS）提取的相同物候图相比效果良好（SoS 的 RMSE 和 ME 分别为 9.2 d 和 2.9 d，EoS 为 6.4 d 和 -0.9d）。当哨兵-2 时间序列出现长时间（>40 d）的数据缺失时，将过滤不可靠的 LSP 估计值，并激活质量标志指示器；春季和初夏的不连续性（9.2%）低于夏末和秋季（39.4%）。哨兵-2 高分辨率 LSP 地图和 GEE 物候提取方法将为植被监测提供支持，并有助于改进陆地表面模型对北极植被物候的表示。

 Normalized Difference Phenology Index，NDPI

归一化差异物候指数（Normalized Difference Phenology Index，NDPI）是一种用于测量植被物候变化的指数。物候是指植物生命周期的不同阶段，如萌芽、开花和结果等。

NDPI通过计算某一时刻的归一化差异植被指数（Normalized Difference Vegetation Index，NDVI）与该地点历史上相同时期的平均NDVI之间的差异来表示物候变化。NDVI是通过测量植被红外辐射与可见光辐射之间的比值来计算的。

NDPI的值范围从-1到1，值趋近于1表示该时期植被比历史上相同时期更为繁茂，值趋近于-1则表示植被比历史上相同时期较少。

NDPI可以用于监测植被的生长状态和季节性变化，对于农业、生态学和气象学等领域的研究具有重要意义。

 NDPI是一种植被指数，用于监测植被物候随时间的变化。它是利用遥感数据计算的，用于量化植被的物候变化。NDPI是从近红外波段和红光波段的反射率差异中导出的，它提供了有关植物生长、凋落和其他物候事件的时间信息。该指数经过归一化处理，以考虑反射率的变化，并且对研究环境因素对植被物候的影响非常有用。研究人员利用NDPI改进了对使用Sentinel-2时间序列数据估算北极地表物候的方法。

文章引用

Descals, A.; Verger, A.; Yin, G.; Pe?uelas, J. Improved Estimates of Arctic Land Surface Phenology Using Sentinel-2 Time Series. Remote Sens. 2020, 12, 3738. https://doi.org/10.3390/rs12223738

数据

我们利用哨兵-2 level-2A 数据生成了 2019 年和 2020 年的 LSP 指标，这些数据提供了 10 米、20 米和 60 米空间分辨率的每日冠顶反射率；从这些数据中，我们使用了 10 米分辨率波段 2、3、4 和 8，以及 20 米波段 12。哨兵-2A 和-2B 多光谱卫星分别于 2015 年和 2017 年发射，在赤道的重访时间为 5 d，随着纬度的增加而缩短。2019 年 5 月 1 日至 9 月 30 日期间，纬度在 70 至 75° 之间的最大重访时间为 1.7 d（平均：0.9 d）（图 1a）。将哨兵-2 估计的 LSP 指标与 500 米 MOD09GAv6 产品估计的相同指标进行了比较 

 

利用哨兵-2 进行地表物候学分析
该资源库包含谷歌地球引擎代码，该代码是为利用哨兵-2 在高纬度地区估算陆地表面物候（LSP）指标（季节开始和结束（SoS 和 EoS））而开发的。

LSP 提取方法是阈值法。

代码：

使用时间序列平滑法估算 SoS 和 EoS（20 天组成和三次插值法）

此代码分为9各部分

SECTION 1 - Define parameters and setup

SECTION 2 - Define function for cubic interpolation

SECTION 3 - Prepare Sentinel-2 data

SECTION 4 - Generate composites

SECTION 5 - Increase window size for empty composites

SECTION 6 - Interpolate time series

SECTION 7 - Estimate phenology metrics

SECTION 8 - Plot results

SECTION 9 - Display legends  

第 1 节 - 定义参数和设置

第 2 节 - 定义三次插值函数

第3节 - 准备哨兵-2数据

第4节 - 生成合成数据

第 5 节 - 增加空合成的窗口大小

第6节 - 插值时间序列

第 7 节 - 估算物候指标

第 8 节 - 绘制结果

第 9 节 - 显示图例  

 平滑后的阈值法。在提取 LSP 指标之前，对时间序列数据进行了平滑处理，这是 LSP 估计中的常见做法，以减少时间序列数据的噪声和不连续性[21]。选择处理步骤的标准是在不影响物候曲线再现的情况下，在 GEE 中实施这些步骤的可行性。过度平滑时间序列可能会导致不切实际的生长季节再现。首先，我们每 20 天应用一个平均半径为 10 d 的移动平均窗口；如果 20 d 合成窗口中的某一像素因缺乏有效观测而为空，则窗口大小增加到 40 d。阈值是根据插值时间序列的振幅而不是每日观测数据估算的，然后分别根据插值时间序列中超过动态阈值的第一天和最后一天来估算SoS和EoS。

///
// Land Surface Estimation (LSP) with Sentinel-2 in the Arctic
//
// This is a demo code for the estimation of the start and end of season (SoS and EoS) with the threshold method with
// time series smoothing (20-day composition and cubic interpolation).
//
// The LSP method is the threhsold method
//
// Adrià Descals - a.descals@creaf.uab.cat
// CREAF - Centre de Recerca Ecològica i Aplicacions Forestals

var textSections = '' +
  "\n SECTION 1 - Define parameters and setup " +
  "\n SECTION 2 - Define function for cubic interpolation " +
  "\n SECTION 3 - Prepare Sentinel-2 data " +
  "\n SECTION 4 - Generate composites  " +
  "\n SECTION 5 - Increase window size for empty composites  " +
  "\n SECTION 6 - Interpolate time series " +
  "\n SECTION 7 - Estimate phenology metrics " +
  "\n SECTION 8 - Plot results " +
  "\n SECTION 9 - Display legends " 
print(textSections)
 
 

//_______________________________________________________________________________________________________________________
// SECTION 1 - Define parameters and setup
Map.setOptions('satellite')

var geom = /* color: #ff0000 */ee.Geometry.Point([17.07398, 67.7578]);
var point = geom // inspect results for this location
var roi = point.buffer(100000).bounds() // display results around the point of interest

var vegIndex = 'evi' // Specify the vegetation index (VI): 'ndvi', 'evi', 'gcc', or 'ndpi'.
var th = 0.5 // Define the percentage of amplitude for the estimation of the threshold 
var threshMin = 0.2 // minimum NDVI value for the reclassification of snow values
var scale = 50 // scale of the analysis
var year1 = 2019 // year of processing

var lag = 20 // temporal resolution of the composites

var startDate = year1+'-04-15'
var endDate = year1+'-12-20'

Map.centerObject(roi,7)


//_______________________________________________________________________________________________________________________
// SECTION 2 - Define function for cubic interpolation

var cubicInterpolation = function(collection,step){ 

var listDekads = ee.List.sequence(1, collection.size().subtract(3), 1);

var colInterp = listDekads.map(function(ii){

  var ii = ee.Number(ii);
  
    var p0 = ee.Image(collection.toList(10000).get(ee.Number(ii).subtract(1)));
    var p1 = ee.Image(collection.toList(10000).get(ii));
    var p2 = ee.Image(collection.toList(10000).get(ee.Number(ii).add(1)));
    var p3 = ee.Image(collection.toList(10000).get(ee.Number(ii).add(2)));
  
    var diff01 = ee.Date(p1.get('system:time_start')).difference(ee.Date(p0.get('system:time_start')), 'day');
    var diff12 = ee.Date(p2.get('system:time_start')).difference(ee.Date(p1.get('system:time_start')), 'day');
    var diff23 = ee.Date(p3.get('system:time_start')).difference(ee.Date(p2.get('system:time_start')), 'day'