GEE非参数趋势分析（Mk-Sen）

朱雀 2024-04-17 06:17 66阅读 0赞

趋势分析是寻找感兴趣的东西正在增加的地方，或者 减少多少。更具体地说，本教程演示了 使用非参数 Mann-Kendall 检测影像中的单调趋势 测试是否存在增加或减少的趋势以及 Sen 的斜率 量化趋势的幅度（如果存在）。本教程还显示 估计 Mann-Kendall 检验统计量的方差，Mann Kendall 检验统计量是 检验是否存在任何趋势，以及统计量的 P 值（假设 正态分布）。  
重要提示：此处介绍的方法 适用于评估单调趋势（即没有季节性的数据） 在离散数据中（即非浮点）。此外，如果应用 本教程中的方法对新数据（即地区、时间范围、来源）可以 需要调整和可视化参数以适应 特别的结果。

**时间序列数据**  
我们将使用来自MOD13A1数据集的MODIS增强植被指数（EVI）的时间序列。此图像集合的每个像素都包含一个时间序列，我们 将计算每个像素的统计信息。假设将集合筛选为 一个季节足以获得具有单调趋势的时间序列。自 检查该假设是否适用于您感兴趣的领域，添加 收集到地图并使用检查器，单击一些点并查看 控制台中显示的系列图表。根据需要调整过滤器。

var mod13 = ee.ImageCollection('MODIS/006/MOD13A1');
    var coll = mod13.select('EVI')
        .filter(ee.Filter.calendarRange(8, 9, 'month'));
    Map.setCenter(-121.6, 37.3, 10);
    Map.addLayer(coll, {
        min: -500, max: 5000, palette: ['white', 'green']}, 'coll');

此集合中的 EVI 时间序列（一个像素）示例如下所示。是 这个像素有趋势吗？更重要的是，从统计学上讲，这是一个 这个像素的显着趋势？请继续阅读，找出答案！

EVI time series in a single pixel  
![在这里插入图片描述][fbf32512bc8a4d36a704bfe0c9ba4774.png]  
**将时间序列连接到自身**  
感兴趣的非参数统计是通过检查所有可能的来计算的 时间序列中唯一值的排序对。如果有 n 个时间点 在序列中，我们需要检查 N（N-1）/2 对 （i， j），i<j，其中 i 和 j 是任意时间索引。（这里我们用（i，j）来表示 由 i 和 j 索引的一对 EVI 图像）。为此，请将集合加入到 本身带有时间过滤器。时态过滤器将通过所有 按时间顺序排列的后期图像。在联接的集合中，每个图像都将 将它之后的所有图像存储在属性中。after

var afterFilter = ee.Filter.lessThan({
        
      leftField: 'system:time_start',
      rightField: 'system:time_start'
    });
    
    var joined = ee.ImageCollection(ee.Join.saveAll('after').apply({
        
      primary: coll,
      secondary: coll,
      condition: afterFilter
    }));

**Mann-Kendall 趋势检验**  
Mann-Kendall 趋势定义为所有货币对的符号之和。 如果时间 j 的 EVI 大于时间 i 的 EVI，则为 1，如果 相反是真的，否则为零（如果它们相等）。计算方式 循环访问集合并计算集合中的每个图像和图像中的每个图像。signsign(i, j)ijafter

缩放到您感兴趣的区域并定义大致对称的颜色拉伸 使用地图图层样式设置对话框（即 和 的均值 近似为零）。红色像素呈递减趋势，绿色像素呈递减趋势 像素呈增加趋势。如下图所示 美国加利福尼亚州一个地区的 Mann-Kendall 统计量。地图图钉是 在 Googleplex 的大致位置。该点是 从中提取上述时间序列的点。我们想要 以确定此地图中的哪些像素具有显著趋势。minmax  
![在这里插入图片描述][2f5a03dbb5c4490583c38b2eec0baeed.png]  
您可能还想知道趋势的幅度或斜率 在当前背景下随时间推移的趋势。一种非参数化方式 用森的斜率来评估。

**森氏斜坡**  
Sen 斜率的计算方式与 Mann-Kendall 统计量类似。 但是，不是将所有符号对相加，而是计算斜率 对于所有对。Sen 斜率是所有这些对的中位斜率。在 在下文中，斜率是以天为单位计算的，以避免数值上的微小斜率 （这可能是由于改用纪元时间而产生的）。

要获取 Sen 的截距（如果需要），请计算所有截距并取 中位数。

var epochDate = ee.Date('1970-01-01');
    var sensIntercept = coll.map(function(image) {
        
      var epochDays = image.date().difference(epochDate, 'days').float();
      return image.subtract(sensSlope.multiply(epochDays)).float();
    }).reduce(ee.Reducer.median(), 2);
    Map.addLayer(sensIntercept, {
        min: -6600, max: 20600}, 'sensIntercept');

森的斜坡图如下所示。请注意，该模式类似于 Mann-Kendall 统计量，但不完全相同。此外，还有 哪些像素具有显着趋势的问题，这个问题将 很快得到答复。  
![在这里插入图片描述][79d73c9b6440491ead6a08144b005105.png]  
前面的示例仅用于演示计算。如果你需要 Sen 的斜率和/或截距，您还可以使用 Sen 的斜率减速器，它可能更容易（更少的代码）和更高的效率，但计算 intercept 作为通过中位数的线的 y 值。

**Mann-Kendall 统计量的方差**  
计算 Mann-Kendall 统计量的方差时，由于 数据中可能存在联系（即 等于零）。 计算这些关系可能会有点棘手，需要基于数组 前向差分。请注意，您可以在计算中评论您是否不关心领带和 想要忽略该更正。sign(i, j).subtract(groupFactorSum)kendallVariance

// Values that are in a group (ties).  Set all else to zero.
    var groups = coll.map(function(i) {
        
      var matches = coll.map(function(j) {
        
        return i.eq(j); // i and j are images.
      }).sum();
      return i.multiply(matches.gt(1));
    });
    
    // Compute tie group sizes in a sequence.  The first group is discarded.
    var group = function(array) {
        
      var length = array.arrayLength(0);
      // Array of indices.  These are 1-indexed.
      var indices = ee.Image([1])
          .arrayRepeat(0, length)
          .arrayAccum(0, ee.Reducer.sum())
          .toArray(1);
      var sorted = array.arraySort();
      var left = sorted.arraySlice(0, 1);
      var right = sorted.arraySlice(0, 0, -1);
      // Indices of the end of runs.
      var mask = left.neq(right)
      // Always keep the last index, the end of the sequence.
          .arrayCat(ee.Image(ee.Array([[1]])), 0);
      var runIndices = indices.arrayMask(mask);
      // Subtract the indices to get run lengths.
      var groupSizes = runIndices.arraySlice(0, 1)
          .subtract(runIndices.arraySlice(0, 0, -1));
      return groupSizes;
    };
    
    // See equation 2.6 in Sen (1968).
    var factors = function(image) {
        
      return image.expression('b() * (b() - 1) * (b() * 2 + 5)');
    };
    
    var groupSizes = group(groups.toArray());
    var groupFactors = factors(groupSizes);
    var groupFactorSum = groupFactors.arrayReduce('sum', [0])
          .arrayGet([0, 0]);
    
    var count = joined.count();
    
    var kendallVariance = factors(count)
        .subtract(groupFactorSum)
        .divide(18)
        .float();
    Map.addLayer(kendallVariance, {
        min: 1700, max: 85000}, 'kendallVariance');

**显著性检验**  
Mann-Kendall 统计量对于适当大 样品。假设我们的样本足够大且不相关。在这些之下 假设，Mann-Kendall 统计量的真实均值为零，并且 方差如上计算。要计算标准正态统计量 （）， 将统计量除以其标准差。z 统计量的 P 值 （观察到这种极值的概率）是 1 - P（|z| < Z）。为 对是否存在任何趋势（正或负）的双侧测试 95% 置信水平，将 P 值与 0.975 进行比较。（或者，比较 Z\* 的 z 统计量，其中 Z\* 是逆分布函数 0.975）。z

标准正态分布函数可以在 Earth Engine 中计算得出 error 函数 .分布函数及其逆函数 如下所示以供参考。在这里，我们使用分布函数得到 1 - P（|z| < Z） 并将其与 0.975 进行比较。erf()

// Compute Z-statistics.
    var zero = kendall.multiply(kendall.eq(0));
    var pos = kendall.multiply(kendall.gt(0)).subtract(1);
    var neg = kendall.multiply(kendall.lt(0)).add(1);
    
    var z = zero
        .add(pos.divide(kendallVariance.sqrt()))
        .add(neg.divide(kendallVariance.sqrt()));
    Map.addLayer(z, {
        min: -2, max: 2}, 'z');
    
    // https://en.wikipedia.org/wiki/Error_function#Cumulative_distribution_function
    function eeCdf(z) {
        
      return ee.Image(0.5)
          .multiply(ee.Image(1).add(ee.Image(z).divide(ee.Image(2).sqrt()).erf()));
    }
    
    function invCdf(p) {
        
      return ee.Image(2).sqrt()
          .multiply(ee.Image(p).multiply(2).subtract(1).erfInv());
    }
    
    // Compute P-values.
    var p = ee.Image(1).subtract(eeCdf(z.abs()));
    Map.addLayer(p, {
        min: 0, max: 1}, 'p');
    
    // Pixels that can have the null hypothesis (there is no trend) rejected.
    // Specifically, if the true trend is zero, there would be less than 5%
    // chance of randomly obtaining the observed result (that there is a trend).
    Map.addLayer(p.lte(0.025), {
        min: 0, max: 1}, 'significant trends');

下面显示了“有效”像素（白色）的映射，或 通过测试的像素。如果我们的假设是正确的 并且您对计算感到满意，您可能希望处理这些 像素不同，用作蒙版。请注意，像素 上面显示的时间序列（下图中的红色）没有 显著趋势。p.lte(0.025)p.lte(0.025)  
![在这里插入图片描述][c500882243a545b3a0278689f89410a4.png]  
有人称这个过程为显著性测试。“重要”的像素 满足条件，并假定具有实数 趋势。其他像素被假定为没有趋势，并从 进一步分析。有些人认为这个过程是邪教 （Ziliak 和 McCloskey 2009）。 由您决定！p.lte(0.025)

[fbf32512bc8a4d36a704bfe0c9ba4774.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/7a290e04503b45ff9fa6bf0dbab84fa8.png
[2f5a03dbb5c4490583c38b2eec0baeed.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/097ce767c9004969803a6147af0f3223.png
[79d73c9b6440491ead6a08144b005105.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/de6ae6420c204d0d9c019039df12b5e5.png
[c500882243a545b3a0278689f89410a4.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/581e9017570f45669712b75a513fbd03.png