Google Earth Engine:对NDVI进行惠特克平滑算法进行长时序分析

简介

函数

ee.Array.identity(size)

Arguments:

Returns: Array

transpose(axis1, axis2)

Arguments:

Returns: Array

matrixMultiply(image2)

Arguments:

Returns: Image

matrixSolve(image2)

Arguments:

Returns: Image

arrayFlatten(coordinateLabels, separator)

Arguments:

Returns: Image

arrayReduce(reducer, axes, fieldAxis)

Arguments:

Returns: Image

代码

结果

简介

惠特克(GEE)平滑算法是一种用于时间序列预测的统计方法，特别适用于非线性、非平稳和非高斯的数据。该算法基于广义估计方程，通过最小化残差的平方和来拟合数据并找到最佳的平滑曲线。

GEE平滑算法的主要思想是在时间序列数据中引入一个平滑函数来描述数据的趋势和周期性变化。该平滑函数由一系列基函数的线性组合组成，其中每个基函数具有不同的频率和振幅。通过调整基函数的权重，可以得到最佳的平滑曲线，以最大程度地拟合数据。

在实际应用中，GEE平滑算法通常与其他统计方法结合使用，例如自回归移动平均模型(ARIMA)或指数平滑法。通过将GEE平滑算法与其他方法相结合，可以进一步提高时间序列的预测准确度和稳定性。

总的来说，GEE平滑算法是一种针对非线性、非平稳和非高斯数据的时间序列预测方法，通过引入一个平滑函数来描述数据的趋势和周期性变化，以最大程度地拟合数据。它在实际应用中通常与其他统计方法结合使用，以进一步提高预测的准确度和稳定性。

函数

ee.Array.identity(size)

Creates a 2D identity matrix of the given size.

创建一个给定大小的二维标识矩阵。

Arguments:

size (Integer):

The length of each axis.

Returns: Array

transpose(axis1, axis2)

Transposes two dimensions of an array.

平移数组的两个维度。

Arguments:

this:array (Array):

Array to transpose.

axis1 (Integer, default: 0):

First axis to swap.

axis2 (Integer, default: 1):

Second axis to swap.

Returns: Array

matrixMultiply(image2)

Returns the matrix multiplication A * B for each matched pair of bands in image1 and image2. If either image1 or image2 has only 1 band, then it is used against all the bands in the other image. If the images have the same number of bands, but not the same names, they're used pairwise in the natural order. The output bands are named for the longer of the two inputs, or if they're equal in length, in image1's order. The type of the output pixels is the union of the input types.

返回图像 1 和图像 2 中每对匹配波段的矩阵乘法 A * B。如果图像 1 或图像 2 中只有一个波段，则该波段将与另一幅图像中的所有波段相对应。如果图像中的条带数量相同，但名称不同，则按自然顺序成对使用。输出波段以两个输入波段中较长的一个命名，如果两个输入波段长度相等，则按图像 1 的顺序命名。输出像素的类型是输入类型的组合。

Arguments:

this:image1 (Image):

The image from which the left operand bands are taken.

image2 (Image):

The image from which the right operand bands are taken.

Returns: Image

matrixSolve(image2)

Solves for x in the matrix equation A * x = B, finding a least-squares solution if A is overdetermined for each matched pair of bands in image1 and image2. If either image1 or image2 has only 1 band, then it is used against all the bands in the other image. If the images have the same number of bands, but not the same names, they're used pairwise in the natural order. The output bands are named for the longer of the two inputs, or if they're equal in length, in image1's order. The type of the output pixels is the union of the input types.

求解矩阵方程 A * x = B 中的 x，如果 A 对图像 1 和图像 2 中每对匹配的波段都是过确定的，则找到最小二乘法解。如果图像 1 或图像 2 中只有一个波段，则使用该波段与另一幅图像中的所有波段进行比对。如果图像中的波段数相同，但名称不相同，则按自然顺序成对使用。输出波段以两个输入波段中较长的一个命名，如果两个输入波段长度相等，则按图像 1 的顺序命名。输出像素的类型是输入类型的组合。

Arguments:

this:image1 (Image):

The image from which the left operand bands are taken.

image2 (Image):

The image from which the right operand bands are taken.

Returns: Image

arrayFlatten(coordinateLabels, separator)

Converts a single-band image of equal-shape multidimensional pixels to an image of scalar pixels, with one band for each element of the array.

将等形多维像素的单波段图像转换为标量像素图像，阵列中的每个元素对应一个波段。

Arguments:

this:image (Image):

Image of multidimensional pixels to flatten.

coordinateLabels (List):

Name of each position along each axis. For example, 2x2 arrays with axes meaning 'day' and 'color' could have labels like [['monday', 'tuesday'], ['red', 'green']], resulting in band names'monday_red', 'monday_green', 'tuesday_red', and 'tuesday_green'.

separator (String, default: "_"):

Separator between array labels in each band name.

Returns: Image

arrayReduce(reducer, axes, fieldAxis)

Reduces elements of each array pixel.

减少每个阵列像素的元素。

Arguments:

this:input (Image):

Input image.

reducer (Reducer):

The reducer to apply.

axes (List):

The list of array axes to reduce in each pixel. The output will have a length of 1 in all these axes.

fieldAxis (Integer, default: null):

The axis for the reducer's input and output fields. Only required if the reducer has multiple inputs or outputs.

Returns: Image

代码

//加载研究区
var geometry = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[113.44773227683973, 38.6708907304602],
          [113.44773227683973, 38.64783484482313],
          [113.47588474266004, 38.64783484482313],
          [113.47588474266004, 38.6708907304602]]], null, false);

// 将 qa 位图像转换为标志的辅助函数
function extractBits(image, start, end, newName) {
    // 计算我们需要提取的比特。
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    // 返回提取的质量保证位的单波段图像，并为该波段命名。
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
}

// 在输入矩阵上获取指定阶次的差分矩阵的函数。将矩阵和阶次作为参数
function getDifferenceMatrix(inputMatrix, order){
    var rowCount = ee.Number(inputMatrix.length().get([0]));
    var left = inputMatrix.slice(0,0,rowCount.subtract(1));
    var right = inputMatrix.slice(0,1,rowCount);
    if (order > 1 ){
        return getDifferenceMatrix(left.subtract(right), order-1)}
    return left.subtract(right);
};

// 将数组图像解包为图像和波段
// 以数组图像、图像 ID 列表和乐队名称列表为参数
function unpack(arrayImage, imageIds, bands){
    
    function iter(item, icoll){
        
        function innerIter(innerItem, innerList){
            return ee.List(innerList).add(ee.String(item).cat("_").cat(ee.String(innerItem)))}
        
        var temp = bands.iterate(innerIter, ee.List([]));
        return ee.ImageCollection(icoll)
            .merge(ee.ImageCollection(ee.Image(arrayImage).select(temp,bands).set("id",item)))}

    var imgcoll  = ee.ImageCollection(imageIds.iterate(iter, ee.ImageCollection([])));
    return imgcoll}


// 用于计算回归结果的反对数比率并转换回百分比单位的函数
function inverseLogRatio(image) {
  var bands = image.bandNames();
  var t = image.get("system:time_start");
  var ilrImage = ee.Image(100).divide(ee.Image(1).add(image.exp())).rename(bands);
  return ilrImage.set("system:time_start",t);
}

function whittakerSmoothing(imageCollection, isCompositional, lambda){
  // 快速配置以设置默认值
  if (isCompositional === undefined || isCompositional !==true) isCompositional = false;
  if (lambda === undefined ) lambda = 10;
  // 程序启动  

  var ic = imageCollection.map(function(image){
     var t = image.get("system:time_start");
    return image.toFloat().set("system:time_start",t);
  });

  var dimension = ic.size();
  var identity_mat = ee.Array.identity(dimension);
  var difference_mat = getDifferenceMatrix(identity_mat,3);
  var difference_mat_transpose = difference_mat.transpose();
  var lamda_difference_mat = difference_mat_transpose.multiply(lambda);
  var res_mat = lamda_difference_mat.matrixMultiply(difference_mat);
  var hat_matrix = res_mat.add(identity_mat);

  
  // 备份原始数据
  var original = ic;

  // 获取原始图像属性
  var properties = ee.List(ic.iterate(function(image, list){
    return ee.List(list).add(image.toDictionary());
  },[]));
  
  var time = ee.List(ic.iterate(function(image, list){
    return ee.List(list).add(image.get("system:time_start"));
  },[]));
  
  // 如果数据是合成的
  // 计算图像在 0 到 100 之间的对比率。首先
  // 夹在 delta 和 100-delta 之间，其中 delta 是一个很小的正值。
  if (isCompositional){
    ic = ic.map(function(image){
      var t = image.get("system:time_start");
      var delta = 0.001;
      var bands = image.bandNames();
      image = image.clamp(delta,100-delta);
      image = (ee.Image.constant(100).subtract(image)).divide(image).log().rename(bands);
      return image.set("system:time_start",t);
    });
  }

  var arrayImage = original.toArray();
  var coeffimage = ee.Image(hat_matrix);
  var smoothImage = coeffimage.matrixSolve(arrayImage);
  
  var idlist = ee.List(ic.iterate(function(image, list){
    return ee.List(list).add(image.id());
  },[]));
  var bandlist = ee.Image(ic.first()).bandNames();
  
  var flatImage = smoothImage.arrayFlatten([idlist,bandlist]);
  var smoothCollection = ee.ImageCollection(unpack(flatImage, idlist, bandlist));
  
  if (isCompositional){
    smoothCollection = smoothCollection.map(inverseLogRatio);
  }
  // 通过添加后缀fitted获得新的乐队名称
  var newBandNames = bandlist.map(function(band){return ee.String(band).cat("_fitted")});
  // 重新命名平滑图像中的波段
  smoothCollection = smoothCollection.map(function(image){return ee.Image(image).rename(newBandNames)});
  
  // 一个非常笨的方法，可以flatten谷歌地球引擎生成的 ID，这样就可以将两张图片合并为图表了
  var dumbimg = arrayImage.arrayFlatten([idlist,bandlist]);
  var dumbcoll = ee.ImageCollection(unpack(dumbimg,idlist, bandlist));
  var outCollection = dumbcoll.combine(smoothCollection);
  
  var outCollectionProp = outCollection.iterate(function(image,list){
      var t = image.get("system:time_start")
    return ee.List(list).add(image.set(properties.get(ee.List(list).size())));
  },[]);

  var outCollectionProp = outCollection.iterate(function(image,list){
    return ee.List(list).add(image.set("system:time_start",time.get(ee.List(list).size())));
  },[]);


  var residue_sq = smoothImage.subtract(arrayImage).pow(ee.Image(2)).divide(dimension);
  var rmse_array = residue_sq.arrayReduce(ee.Reducer.sum(),[0]).pow(ee.Image(1/2));
  
  var rmseImage = rmse_array.arrayFlatten([["rmse"],bandlist]);
  
  return [ee.ImageCollection.fromImages(outCollectionProp), rmseImage];
}



var ndvi =ee.ImageCollection("NOAA/VIIRS/001/VNP13A1").select('NDVI').filterDate("2019-01-01","2019-12-31");
// 去除屏蔽像素
ndvi = ndvi.map(function(img){return img.unmask(ndvi.mean())});

var ndvi =  whittakerSmoothing(ndvi)[0];

// 添加图表
print(ui.Chart.image.series(
  ndvi.select(['NDVI', 'NDVI_fitted']), geometry, ee.Reducer.mean(), 500)
    .setSeriesNames(['NDVI', 'NDVI_fitted'])
    .setOptions({
      title: 'smoothed',
      lineWidth: 1,
      pointSize: 3,
}));