【GEE】遥感数据趋势分析Sen+mk

秒速五厘米 2024-04-17 06:17 82阅读 0赞

![在这里插入图片描述][906213cc805542b6ac0e1702f4a2d62b.png]

Map.centerObject(table);
    
    // 定义时间范围
    var stary = 2001, endy = 2023;
    //NDVI图像集合
    var NDVICL = ee.ImageCollection(ee.List.sequence(stary, endy).map(function(year) {
        
      // 定义每年的开始和结束日期
      var startd = ee.Date.fromYMD(year, 1, 1);
      var endd = ee.Date.fromYMD(year, 12, 31);
      // 加载NDVI数据集，按日期过滤，并计算最大值
      return ee.ImageCollection('MODIS/006/MOD13A1')
        .filterDate(startd, endd)
        .select('NDVI')
        .max()
        .addBands(ee.Image.constant(year).toFloat().rename('year'));
    }));
    print(NDVICL)
    // 计算Sen斜率并剪裁
    var senSlope = NDVICL.select(['NDVI', 'year'])
        .reduce(ee.Reducer.sensSlope())
        .clip(table);
    // 可视化Sen斜率
    var vis = {
        bands: ['slope'], min: -0.5, max: 0.5, palette: ['red', 'white', 'green']};
    Map.addLayer(senSlope.select('slope'), vis, "Sen斜率 ");
    
    // MK检验 
    // create ones image
    var ones = ee.Image(1);
    // create zeros image
    var zeros = ee.Image(0);
    var eps = 0.01; // Set a precision value
    // 创建一个列表，包含年度NDVI图像的list
    var listofimg = NDVICL.toList(NDVICL.size());
    // 计算MK检验趋势期望均值
    var mkm = ee.ImageCollection(ee.List.sequence(0, listofimg.size().subtract(2)).map(
        function (xi) {
        
            return ee.ImageCollection.fromImages(
                ee.List.sequence(ee.Number(xi).add(1), listofimg.size().subtract(1)).map(
                    function (xj) {
        
                        var diff_ij = ee.Image(listofimg.get(xj)).select('NDVI')
                            .subtract(ee.Image(listofimg.get(xi)).select('NDVI'));
                        var diff = diff_ij.where(diff_ij.abs().lt(eps), zeros);
                        return diff.where(diff.gt(0), ones)
                            .where(diff.lt(0), ones.multiply(-1)).rename('ES');
                         })).sum();
        })).sum();
    // 时序长度 
    var n = ee.ImageCollection(NDVICL)
        .reduce('count')
        .select('NDVI_count')
        .rename('count');
    // 计算方差 var(S) = (n * (n - 1) * (2 * n + 5)) / 18
    var varS = n.multiply(
        n.subtract(1).multiply(
            n.multiply(2).add(5)))
        .divide(18).rename('varS');
    
    // 估计z统计量
    var mkz = mkm.where(mkm.abs().lt(eps), 0);
    var zcore = mkz
        .where(mkz.gt(eps), mkz.subtract(1).divide(varS.sqrt()))
        .where(mkz.lt(-eps), mkz.add(1).divide(varS.sqrt())).rename('zcore')
        .clip(table);
    
    // 100(1 - alpha)% 双边置信区间的Sen's Slope
    //alpha = 0.05，则Z[1 - alpha / 2] = Z[0.975] = 1.960
    // 计算Sen斜率的显著性
    var signif = senSlope.select(0).lt(0).and(zcore.abs().gte(1.960)).rename('signif');
    var siginc = senSlope.select(0).gt(0).and(zcore.abs().gte(1.960)).rename('increase');
    
    signif.updateMask(signif);
    siginc.updateMask(siginc);
    
    signif = signif.multiply(-1).add(siginc).clip(table);
    
    // 可视化MK Z统计量
    Map.addLayer(zcore, {
        min: -3, max: 3, palette: ['blue', 'white', 'red']}, 'Z值 ');
    
    // 可视化显著的趋势变化
    Map.addLayer(signif, {
        min: -1, max: 1, palette: ['#FFA07A', 'white', '#3CB371']}, '显著趋势 ');

[906213cc805542b6ac0e1702f4a2d62b.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/9b4a4932f03841118995b6156ca56806.png