NDVI时间序列分析原理与实现(CV和Sen+MK趋势分析)

订正：结合论文和R语言帮助，确定sen's slope的代码原理是median，不是均值，前面老版本的推文有谬误，特此更正！

总结两种常用的NDVI时间序列栅格分析方法，变异系数和Sen+MK趋势分析。简单介绍原理和实现代码。

实验数据均为年际变化数据，在CV中使用的是多波段ENVI数据，Sen+MK趋势分析中是一个文件夹中的多幅单波段影像。

变异系数（Coefficient of Variation, CV）

原理

变异系数是衡量观测序列值变异程度的一个统计量，可以很好地反映空间数据在时间序列上变化的差异程度，评价数据时间序列的稳定性。计算公式为：

式中： 为标准差， 为算术平均值。CV越大，数据波动越大。

R语言实现代码

需要注意的是，这里输入的是一个多波段的影像，ENVI格式。

raster('./test/testdata.enp', band = 2) -> test
#先生成一个tibble，不然后面会报错
test[] %>% 
  as_tibble()

#每个波段读取为数据框
lapply(1:24, function(x) {   #波段数需要根据情况修改
  raster('./test/testdata.enp',     #输入多波段影像
         band = x) %>% 
    `[` %>% 
    as_tibble() %>% 
    mutate(id = as.numeric(row.names(.)))
}) %>% 
  bind_rows() -> df

df
df %>% 
  arrange(id) %>% 
  nest(-id) -> nested_df

nested_df$data[1][[1]] %>% 
  pull(value) %>% 
  cv()  #变异系数计算

nested_df %>% 
  mutate(cv = map_dbl(data, function(x) cv(x$value))) %>% 
  dplyr::select(-data) %>% 
  pull(cv) -> cv
cv %>% 
  raster(vals = ., nrow = 24, ncol = 31,     #像元行列数
         crs = "+proj=longlat +datum=WGS84 +no_defs",
         ext = extent(test)) -> cv_raster

plot(cv_raster)

# 输出
cv_raster %>% 
  writeRaster("cv", overwrite = TRUE)

某区域NDVI变异系数图

Sen+MK趋势分析

Sen斜率估计用于计算趋势值，通常与MK非参数检验结合使用，即先计算Sen趋势值，然后使用MK方法判断趋势显著性。

原理

Theil-Sen Median方法又被称为Sen斜率估计，是一种稳健的非参数统计的趋势计算方法。该方法计算效率高，对于测量误差和离群数据不敏感，常被用于长时间序列数据的趋势分析中。

式中： 和 为时间序列数据。β大于0表示时间序列呈现上升趋势；β小于0表示时间序列呈现下降趋势。

R语言实现代码

注意的是，这里输入的是存储在文件夹中的多幅单波段TIFF影像。

library(sp)
library(raster)
library(rgdal)
library(trend)
#输入一个文件夹内的单波段TIFF数据，在这里是历年的NDVI年最大值
flnames <- list.files(path = './yearmax/', pattern = '.tif$')
fl <- paste0("./yearmax/", flnames)
firs <- raster(fl[1])
for (i in 2:35) {
  r <- raster(fl[i])
  firs <- stack(firs, r)
}

fun <- function(y){
  if(length(na.omit(y)) <35) return(c(NA, NA, NA))   #删除数据不连续含有NA的像元
  av <- mean(y, na.rm=T)
  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
  return(c(av, slope, MK_test))
}

e <- calc(firs, fun)   #栅格计算
e_mean <- subset(e,1)  #提取均值图层
e_slope <- subset(e,2)   #提取sen斜率
e_MKtest <- subset(e,3)   #提取p值

plot(e_mean)
plot(e_slope)
plot(e_MKtest)

writeRaster(e_mean, "e_mean.tif", format="GTiff", overwrite=T)
writeRaster(e_slope, "e_slope.tif", format="GTiff", overwrite=T)
writeRaster(e_MKtest, "e_MKtest.tif", format="GTiff", overwrite=T)

Sen斜率估计栅格

参考文献

1. https://blog.csdn.net/qq_25102303/article/details/106927712
2. https://jingyan.baidu.com/article/ab69b270e00f6d6da6189f42.html
3. https://www.rdocumentation.org/packages/trend/versions/1.1.4/topics/sens.slope