NDVI时间序列分析之Sen+MK分析全过程梳理
Sen斜率估计用于计算趋势值，通常与MK非参数检验结合使用，即先计算Sen趋势值，然后使用MK方法判断趋势显著性。
原理
Theil-Sen Median方法又被称为Sen斜率估计，是一种稳健的非参数统计的趋势计算方法。该方法计算效率高，对于测量误差和离群数据不敏感，常被用于长时间序列数据的趋势分析中。
β = m e d i a n ( x j − x i j − i ) , ∀ j > i \beta=median(\frac{x_j-x_i}{j-i}),∀j>i β=median(j−ixj−xi),∀j>i
式中： x j x_j xj 和 x i x_i xi 为时间序列数据。β大于0表示时间序列呈现上升趋势；β小于0表示时间序列呈现下降趋势。
Mann-Kendall属于非参数检验方法，与其他参数检验的方法相比，不需要样本遵从一定的分布，受异常值干扰小，更适合顺序变量。Mann-Kendall检验已经在水文、气象趋势变化相关研究中得到了大量的成功应用，用于判断径流、降水、气候等的变化趋势的显著性。
R语言Sen+MK并行计算
R语言使用Terra包进行栅格并行计算，利用trend包sen.slope函数进行sen+mk的计算。对于下面代码的阅读，一定要看代码帮助！代码如下：
library(terra)library(trend)#输入一个文件夹内的单波段TIFF数据，在这里是历年的NDVI年最大值flnames <- list.files(path = './ChinaYearMean/', pattern = '.tif$')fl <- paste0("./ChinaYearMean/", flnames)firs <- rast(fl)#Sen+MK计算fun_sen <- function(x){if(length(na.omit(x)) <34) return(c(NA, NA, NA)) #删除数据不连续含有NA的像元MK_estimate <- trend::sens.slope(ts(na.omit(x), start = 1982, end = , frequency = 1), conf.level = 0.95) #Sen斜率估计slope <- MK_estimate$estimateMK_test <- MK_estimate$p.valueZs <- MK_estimate$statisticreturn(c(Zs, slope, MK_test))}firs_sen = app(firs, fun_sen, cores=4 )names(firs_sen) = c("Z", "slope", "p-value")#显示输出plot(firs_sen)writeRaster(firs_sen, filename = "./firs_sen.tif", names=firs_sen@ptr[["names"]])
如果还没明白，下面二维码推荐一个培训班，有专门的R语言老师带着学习：
NDVI趋势制图
前面R语言计算完slope和Z值后，根据这两个结果就可以进行NDVI趋势制图了。
变化趋势划分
结合 S N D V I S_{NDVI} SNDVI和Z统计量划分NDVI变化趋势：
slope -0.0005~0.0005稳定区域大于或等于0.0005植被改善区域小于-0.0005为植被退化区域 Z统计量 置信水平0.05Z绝对值大于1.96显著Z绝对值小于等于1.96不显著
Slope划分
Slope被划分为三级： S N D V I ≤ − 0.0005 S_{NDVI} \le -0.0005 SNDVI≤−0.0005 植被退化 − 0.0005 ≥ S N D V I ≥ 0.0005 -0.0005 \ge S_{NDVI} \ge 0.0005 −0.0005≥SNDVI≥0.0005 植被生长稳定 S N D V I ≥ 0.0005 S_{NDVI} \ge 0.0005 SNDVI≥0.0005 植被改善 使用重分类（Reclassify）对slope进行划分 由于slope.tif文件研究区范围外的值非空，所以在这里先裁剪了一下裁剪所用矢量和栅格数据坐标系需要一致，否则范围容易出错统一使用了WGS84地理坐标系作为空间参考使用Model builder构建地理处理流 重分类结果： -1退化0稳定1改善Z值划分
对Z值进行重分类，确定显著性 ∣ Z s ∣ ≤ 1.96 \lvert Z_s\rvert \le 1.96 ∣Zs∣≤1.96 未通过95%置信度检验，不显著 ∣ Z s ∣ ≥ 1.96 \lvert Z_s\rvert \ge 1.96 ∣Zs∣≥1.96 通过95%置信度检验，显著 重分类结果： 1不显著2显著变化趋势计算
使用栅格计算器将Slope和Z值计算结果相乘，最后得到趋势变化划分
-2严重退化-1轻微退化0稳定不变1轻微改善2明显改善
对计算完成后的栅格属性表里面不同种类的像元数进行统计，如此即可获得下面的NDVI变化趋势统计图表。
