Earth Engine 时间序列平滑及编程体会

Posted 2021-05-02 GeoTalks

撰文 | 嵩岘山人，熊氏阿回      授权发布

编辑 | GeoTalks

1 目标

基于 Google Earth Engine 平台，获取一年期的 Sentinel-2（哨兵2号） 数据，并计算相应的 NDVI 时间序列以及基于 Satizky-Golay 滤波方法的 NDVI 时间序列的平滑结果，使对于地图上每个被点击的点，都能得到对应的 NDVI 时间序列及其平滑结果的图表显示。

2 数据源

Sentinel-2 卫星由欧洲空间局(ESA)研制开发，其中 Sentinel-2A 在2015年6月23日发射升空，Sentinel-2A 搭载的高分辨率多光谱成像仪(MSI)具有13个光谱波段，时间及空间分辨率均较高（各波段空间分辨率范围是10m~60m），且多光谱成像能力强，为全球土地和植被研究提供了一个全新的视角，在监测地球植被健康状况方面有广阔应用。Sentinel-2 较高的空间和时间分辨率也适于计算 NDVI 时间序列。Sentinel-2 数据的 NDVI 计算公式如下：

Sentinel2: NDVI = (band8 - band4) / (band8 + band4)

3 Satizky-Golay 时间序列平滑方法

理论上，由于植被冠层随时间变化幅度较小，则绘制得到的 NDVI 曲线应是一条连续平滑的曲线。但在现实情况下，由于云层干扰、数据传输误差、二向性反射或地面冰雪的影响，在 NDVI 曲线中总是会有明显的突升或突降。因此在得到 NDVI 曲线之后一般还要对其进行平滑处理。

S-G 平滑方法主要是通过选取某个点附近固定个数的点拟合一个多项式，通过该多项式拟合出这个点的平滑值。该方法对于 NDVI 的尺度以及传感器的类型没有严格的要求。通过设置拟合效果指数，利用 S-G 滤波不断的进行迭代，当该指数达到局部最小时，则标志着这一序列点的拟合效果最佳。

4 GEE 实现

本编程基于 GEE 中的 javascript 实现。

/*The following codes aim at computing the NDVI time series from the  S2 images of a year and drawing the charts of the data as well as its smoothing result by the Savitzky-Golay method,for each click on the map*/

4.1 NDVI 的计算及 Satizky-Golay 平滑算法 

4.1.1 获取 Sentinel-2 数据

通过 Earth Engine 自带的 ImageCollection 函数在线获取 Sentinel-2 自2016年6月1日到2017月6月1日的一年期数据集。

//Getting the S2 image collection of a year

 var l8 = ee.ImageCollection('COPERNICUS/S2')

 .filterDate('2016-06-01', '2017-06-01');

4.1.2 计算 NDVI

计算所获数据集的 NDVI 时间序列数据， Sentinel-2 的 NDVI 计算公式见于“3. 数据源”，在实际操作中可直接借助 image.normalizedDifference 函数实现。

//Computing the NDVI by band8 and band4

var ndvi = l8.map(function(image) {

  return image.select().addBands(image.normalizedDifference(['B8', 'B4']));

});

var vis = {min: 0, max: 1, palette: [

  'FFFFFF', 'CE7E45', 'FCD163', '66A000', '207401',

  '056201', '004C00', '023B01', '012E01', '011301'

]};

Map.addLayer(ndvi, vis, 'NDVI');

Map.setCenter(-94.84497, 39.01918, 8);

4.1.3 获取对应点的 NDVI 数据

因为对于每一个在地图上被点击的点，除了计算该位置的 NDVI 之外，还要用相关算法得到其平滑后的结果，为方便后续平滑函数的处理，需要从之前的 NDVI 计算结果中提取出对应某点的 NDVI 时间序列。

//Get the NDVI time series from the clicked point region

   var y = ndvi.getRegion(point, 250).slice(1) // remove header

      .map(function(item) { return ee.Number(ee.List(item).get(-1)).divide(1e4) })

4.1.4 S-G 平滑

使用 Satizky-Golay 滤波法对 NDVI 时间序列进行平滑，输入的数据是 4.1.3 中得到的对应某点的 NDVI 时间序列 y，最终得到平滑后的时间序列 yValues 。

//The following is the realization of Savitzky-Golay soomthing method

   var window_size = 11

   var half_window = (window_size - 1)/2

   var deriv = 0

   var order = 2

   var order_range = ee.List.sequence(0,order)

   var k_range = ee.List.sequence(-half_window, half_window)

   var b = ee.Array(k_range.map(function (k) { return order_range.map(function(o) { return ee.Number(k).pow(o)})}))

   print(b)

   

   var mPI = ee.Array(b.matrixPseudoInverse())

   print(mPI)

   var impulse_response = (mPI.slice({axis: 0, start: deriv, end: deriv+1})).project([1])

   print(impulse_response)

   

   var y0 = y.get(0)

   var firstvals = y.slice(1, half_window+1).reverse().map(

     function(e) { return ee.Number(e).subtract(y0).abs().multiply(-1).add(y0) }

   )

   print(firstvals)

   

   var yend = y.get(-1)

   var lastvals = y.slice(-half_window-1,-1).reverse().map(

     function(e) { return ee.Number(e).subtract(yend).abs().add(yend) }

   )

   print(lastvals)

   

   var y_ext = firstvals.cat(y).cat(lastvals)

   print(y_ext.length())

   

   var runLength = ee.List.sequence(0, y_ext.length().subtract(window_size))

   var smooth = runLength.map(function(i) {

     return ee.Array(y_ext.slice(ee.Number(i), ee.Number(i).add(window_size))).multiply(impulse_response).reduce("sum", [0]).get([0])

   })

   print(smooth)

   

   var yValues = ee.Array.cat([y, smooth], 1);

4.2 相关图件的绘制 

4.2.1 创建基础图件

创建一些基础图件，包括图表的相关注记、用户与图形界面交互时的提示信息、经纬度信息等。

// Create a panel to hold our widgets.

var panel = ui.Panel();

panel.style().set('width', '300px');

// Create an intro panel with labels.

var intro = ui.Panel([

  ui.Label({

    value: 'Chart Inspector',

    style: {fontSize: '20px', fontWeight: 'bold'}

  }),

  ui.Label('Click a point on the map to inspect.')

]);

panel.add(intro);

// Create panels to hold lon/lat values.

var lon = ui.Label();

var lat = ui.Label();

panel.add(ui.Panel([lon, lat], ui.Panel.Layout.flow('horizontal')));

4.2.2 创建 Click 响应函数

Click 函数的目的是对用户在地图上点击的每一个点进行响应，获取对应点的经纬度坐标值，并根据该点位置进行 NDVI 时间序列的计算、平滑、图表绘制等内容。

// Register a callback on the default map to be invoked when the map is clicked.

Map.onClick(function(coords) {

  // Update the lon/lat panel with values from the click event.

  lon.setValue('lon: ' + coords.lon.toFixed(2)),

  lat.setValue('lat: ' + coords.lat.toFixed(2));

  // Add a red dot for the point clicked on.

  var point = ee.Geometry.Point(coords.lon, coords.lat);

  var dot = ui.Map.Layer(point, {color: 'FF0000'});

  Map.layers().set(1, dot);

  // The rest parts of this function are designed to be filled with the codes aiming at realizing the calulation and smoothing of the NDVI time series, as well as drawing charts for the computing results.

});

Map.style().set('cursor', 'crosshair');

// Add the panel to the ui.root.

ui.root.insert(0, panel);

4.2.3 绘制 NDVI 图表

根据计算出来的 NDVI 数据绘制相应图表，注记、字体等参数可在函数中设置。

// Create an NDVI chart.

  var ndviChart = ui.Chart.image.series(ndvi, point, ee.Reducer.mean(), 500);

  ndviChart.setOptions({

    title: 'NDVI Over A Year',

    vAxis: {title: 'NDVI'},

    hAxis: {title: 'date', format: 'MM-yy', gridlines: {count: 7}},

  });

  panel.widgets().set(2, ndviChart);

4.2.4 绘制平滑图表

根据 S-G 方法平滑后的 NDVI 数据绘制相应图表，注记、字体等参数可在函数中设置。

//Create the smoothing chart by using the yValues gotten from the S-G method

   var smoothing_chart = ui.Chart.array.values(yValues, 0).setSeriesNames(['Raw', 'Smoothed']).

   setOptions(

     {

       title: 'Smoothing NDVI Over Time',

       vAxis: {title: 'NDVI value'},

       hAxis: {title: 'date', format: 'MM-yy', gridlines: {count: 7}}

       }

   );

   panel.widgets().set(3, smoothing_chart);

4.3 运行结果示例 

运行上述代码，在地图上点击某一点，GEE 就会自动计算 Sentinel-2 在该点处一年期的 NDVI 时间序列以及对应的 S-G 方法的平滑结果，并在左侧栏中绘制出相应图表。

图1 GEE 运行结果示例

5 编程体会

5.1 初识地球引擎 

这是我第一次使用 Google Earth Engine 完成一项任务（同时也是第一次接触 JavaScript 语句），但并不是第一次听说 GEE 了。我清楚得记得第一次听说 GEE 是在上学期的一次讲座上，那是美国堪萨斯大学的李新功教授在南大做的一次讲座，主要就是介绍他利用 GEE 所开展的研究工作，具体的内容已经不记得了，但那场讲座给我留下了两点印象：一是 GEE 对数据的计算能力非普通计算机可比；二是它的应用在国内尚未普及，离本科学生更是遥远。

孰料在将近一年之后，因缘际会，我也开始使用 GEE 处理数据，我想能接触这款利器既是我的幸运，也从侧面反映了 GEE 应用的扩散趋势。

GEE 既是对我来说既是新的工具，也代表着一门陌生的语言——JavaScript。从大一初学 C 开始，到后来的 C++、SQL、C#，到这学期正在学习的 Python和Matlab，我逐渐积累了两条经验：

一是“工具”思想，也就是把编程语言当做工具来用（相应也要找到好的指南类的书），根据实战中的目的来查找相关“工具”，这样能够很快入门，如果仅是将其作为教科书式的知识去学习的话，耗时长，而且不一定能够实战。

二是“黑箱”思想：对于许多数据处理的中间过程来说，如果有封装好的包或函数，可以直接拿来用，未必要穷究内部原理，这样会提高效率。依着这两条个人的经验，我先对 GEE 做了概览，对其界面、语法等有一大致了解，GEE 中带有文档、教程、论坛等资源，另外在知乎等网站上也有关于GEE 的专题，这些为了解 GEE 提供了很好的入门。

5.2 结语 

入门后紧接着就是针对 NDVI 的计算和平滑进行编程了，由于并没有学过 JavaScript，因此“工具”和“黑箱”是我进行解难的初始信条，我开始寻找和读图像、计算 NDVI、S-G 滤波法、图表绘制有关的函数或代码，但其中仍然遇到不少困难，印象最深的大概就是在 NDVI 计算和 S-G 平滑两大部分接口处所遇到的问题了，虽然最后只凭借两行代码解决了问题，但“一沙一世界”，在此处的数次反复也折射出不少问题，加深了我对 GEE 的认识。

因为编程前一部分计算得到的 NDVI 实际上是整个图像范围内的序列，然后凭借 EE 自带的绘图函数将所要描绘的数据读出并绘制图表。

但是 S-G 平滑算法要处理的对象是一个一维数组，所以要先从之前得到的 NDVI 影像序列中得到数组。

起初我在计算全图像域的 NDVI 影像序列之外，专门写了计算某一点的 NDVI 影像序列的代码，然后查阅了 EE 中关于 Image 和 Array 的帮助文档，认为只要利用 .toArray() 函数即可将前边的影像序列转化为可利用的数组，但调试过程中并非如此。

之后改用 .toRegion函数对 NDVI 影像序列做处理（对于参数的查找当时也是破费功夫的事），该函数的处理结果 print 出来表现为二维数组的形式，也就是除了所需的一列 NDVI 数据之外，还包含各个波段的信息，现在虽未达到目的地，但是看来也相当接近了——只要用循环遍历的方法将二维数组中那一列 NDVI 信息读进一个一维数组就成功创建接口了。但是在写完遍历函数之后发现并非如此，对于前边的“二维数组”在使用数组下标引用， print 出来的元素却是显示 undefined，几次检查代码并无语法错误，这让我疑惑前边得到的真是一个二维数组吗？

//最初的接口部分的代码

var y1=ndvi_click.getRegion(

     {

       geometry: point,

       scale: 90

     }

   );

   print(y1.length())

   print(y1)//print(y1)输出的是类似二维数组的形式

   print(typeof y1)//显示y1的type 是 object，对象

   print(y1[3])

   

   //Get the NDVI array

   var y=[];

   for(var i=1;i<y1.length;i++){

     for(var j=0;j<y1[i].length;j++){

       if(j==(y1[i].length-1)){

         y.push(y1[i][j])

       }

     }

   }

   print(y.length())//显示 undefined

   print(y)//显示 undefined

事实上，将 .getRegion() 函数得到的变量 y1 用 print(typeof y1) 查询得到的属性是 object， 而非 array 或 list，查阅 GEE 的帮助文档中的 Client & Server部分，对 object 的解释是It's an ee.computedObject, which means it's a proxy object for something on the server. 也就是说 EE 在服务器端封装的一个数据类型，和客户端自己定义的二维数组应该还是有差别的，并且 EE 的建议是客户端和服务器端定义的变量最好不要混用，否则会造成混乱或者效率问题，在后续处理中，可以用 .getInfo()方法可以得到二维数组提出其中的 NDVI 序列，但是仍然不能适应 S-G 平滑代码中 EE 自带的一些函数。最终，还是找到 EE 自带的函数求出该点对应的 NDVI 序列从而解决了这个接口问题。

这个过程一方面让我认识到 EE 的语言并不等价于 JavaScript，其有自己的数据结构，有的时候如果搞不清楚这个问题会对整个编程产生影响；另一方面，我感到优秀的工具书或资源序列是多么地重要，在编程过程中引用的一些代码，如果能有工具书查询到其中 EE 中各种函数的参数和返回值的类型，就会有助于理解编程过程中出现的各种问题，当然，这和对 EE 的熟悉程度也有关系，希望日后在这方面能有进步。

参考资料

[1] Gorelick N, et al. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment. 2017.

[2] Ma M G and Veroustraete F.2006.Reconstructing pathfinder AVHRR land NDVI time-series data for the Northwest of China. Advances in Space Research, 37: 835—840.

[3] 赵英时. 2003. 遥感应用分析原理与方法. 北京：科学出版社.

[4] 边金虎, 等. MODIS植被指数时间序列Savitzky-Golay滤波算法重构. 遥感学报. 1007-4619 (2010) 04-725-17.

[5] 岳桢干. 欧洲Sentinel-2A卫星即将大显身手——“哥白尼”对地观测计划简介（上）. 红外. Vol. 36, No. 8, Aug 2015.

[6] 知乎专题：GEE开发. https://zhuanlan.zhihu.com/p/29620366. 

[7] 题图来源：http://www.esa.int/spaceinimages/Images/2017/12/Reindeer_Island.

---

嵩岘山人，南京大学地理信息科学系在读

熊氏阿回，美国NASA研究中心，现居美国硅谷

链接：https://zhuanlan.zhihu.com/p/32618854

来源：知乎

感谢作者与GeoTalks读者分享