如何对时间仅存于TIFF文件名的栅格堆栈做逐像素线性回归分析

第一步：从文件名提取年份信息 #

年份仅存储在TIFF文件名中，需先提取并生成对应时间向量，确保和栅格堆栈图层顺序一致：

library(raster)
# 假设栅格文件名格式如"xxx_1986.tif"，请根据实际文件名格式调整正则表达式
file_names <- names(stack)
years <- as.integer(sub(".*_(\\d{4})\\.tif", "\\1", file_names))
# 验证年份数量与栅格层数匹配，避免顺序错误
stopifnot(length(years) == nlayers(stack))

第二步：修正像素级回归函数 #

你之前的函数逻辑有误，需将每个像素的时间序列值与年份做回归，同时处理NA值并应用Cochrane-Orcutt修正：

library(orcutt)

# 定义每个像素的趋势计算函数
trend_fun <- function(x) {
  # 全NA像素直接返回NA，避免报错
  if (all(is.na(x))) {
    return(NA)
  }
  # 构建以年份为自变量的线性回归模型
  lm_model <- lm(x ~ years)
  # 应用Cochrane-Orcutt修正自相关问题
  co_model <- cochrane.orcutt(lm_model)
  # 返回修正后的斜率系数（即趋势值）
  return(co_model$coefficients[2])
}

# 计算所有像素的趋势斜率
slope_raster <- calc(stack, trend_fun)

可选扩展 #

• 普通线性回归（无需修正）：如果不需要处理自相关问题，可简化函数：

simple_trend_fun <- function(x) {
  if (all(is.na(x))) return(NA)
  coef(lm(x ~ years))[2]
}
simple_slope <- calc(stack, simple_trend_fun)

• 多统计量输出：如需同时获取截距、p值等指标，返回向量即可生成多波段栅格：

multi_stat_fun <- function(x) {
  if (all(is.na(x))) return(c(NA, NA, NA))
  lm_model <- lm(x ~ years)
  co_model <- cochrane.orcutt(lm_model)
  # 返回斜率、截距、斜率对应的p值
  c(slope = co_model$coefficients[2],
    intercept = co_model$coefficients[1],
    p_value = summary(co_model)$coefficients[2, 4])
}
multi_stat_raster <- calc(stack, multi_stat_fun)

内容的提问来源于stack exchange，提问作者GeogDay