||
Sen斜率估计用于计算趋势值,通常与MK非参数检验结合使用,即先计算Sen趋势值,然后使用MK方法判断趋势显著性。
Theil-Sen Median方法又被称为Sen斜率估计,是一种稳健的非参数统计的趋势计算方法。该方法计算效率高,对于测量误差和离群数据不敏感,常被用于长时间序列数据的趋势分析中。
β=mean(xj−xij−i),∀j>i
式中:xj 和 xi 为时间序列数据。β大于0表示时间序列呈现上升趋势;β小于0表示时间序列呈现下降趋势。
Mann-Kendall属于非参数检验方法,与其他参数检验的方法相比,不需要样本遵从一定的分布,受异常值干扰小,更适合顺序变量。Mann-Kendall检验已经在水文、气象趋势变化相关研究中得到了大量的成功应用,用于判断径流、降水、气候等的变化趋势的显著性。
R语言使用Raster包进行栅格计算,利用trend包sen.slope函数进行sen+mk的计算。对于下面代码的阅读,一定要看代码帮助!代码如下:
library(sp)
library(raster)
library(rgdal)
library(trend)
library(terra)
#输入一个文件夹内的单波段TIFF数据,在这里是历年的NDVI年最大值
flnames <- list.files(path = './yearmax/', pattern = '.tif$')
fl <- paste0("./yearmax/", flnames)
firs <- raster(fl[1])
for (i in 2:34) {
r <- raster(fl[i])
firs <- stack(firs, r)
}
fun <- function(y){
if(length(na.omit(y)) <34) return(c(NA, NA, NA)) #删除数据不连续含有NA的像元
MK_estimate <- sens.slope(ts(na.omit(y), start = 1982, end = 2015, frequency = 1), conf.level = 0.95) #Sen斜率估计
slope <- MK_estimate$estimate
MK_test <- MK_estimate$p.value
Zs <- MK_estimate$statistic
return(c(Zs, slope, MK_test))
}
e <- calc(firs, fun) #栅格计算
e_Zs <- subset(e,1) #提取Z统计量
e_slope <- subset(e,2) #提取sen斜率
e_MKtest <- subset(e,3) #提取p值
plot(e_Zs)
plot(e_slope)
plot(e_MKtest)
writeRaster(e_Zs, "./Sen+MK/e_Zs.tif", format="GTiff", overwrite=T)
writeRaster(e_slope, "./Sen+MK/e_slope.tif", format="GTiff", overwrite=T)
writeRaster(e_MKtest, "./Sen+MK/e_MKtest.tif", format="GTiff", overwrite=T)
如果还没明白,下面二维码推荐一个培训班,有专门的R语言老师带着学习:
前面R语言计算完slope和Z值后,根据这两个结果就可以进行NDVI趋势制图了。
结合SNDVI和Z统计量划分NDVI变化趋势:
Slope划分过程
Z值重分类
使用栅格计算器将Slope和Z值计算结果相乘,最后得到趋势变化划分
栅格计算器相乘
对计算完成后的栅格属性表里面不同种类的像元数进行统计,如此即可获得下面的NDVI变化趋势统计图表。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-11-5 10:38
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社