GEE进行Theil-Sen Median斜率估计和Mann-Kendall检验

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']}, '显著趋势 ');