如何在R语言中针对三个栅格栈计算逐像素局部回归？

没问题，我来帮你搞定这个逐像素的二元线性回归需求！你需要用s1和s2作为自变量，逐像素对s的时间序列做回归，最终得到每个像素的截距、s1对应的斜率、s2对应的斜率，用terra包可以这样实现：

步骤1：加载包与示例数据 #

首先我们先加载terra包，并调用示例栅格栈来模拟你的数据（三个栅格栈的层数一致，代表相同长度的时间序列）：

library(terra)
# 加载示例栅格栈，模拟时间序列数据
s <- rast(system.file("ex/logo.tif", package="terra"))
s1 <- rast(system.file("ex/logo.tif", package="terra"))
s2 <- rast(system.file("ex/logo.tif", package="terra"))

步骤2：定义逐像素回归函数 #

我们需要定义一个函数，让它接收每个像素的所有时间点数据（s、s1、s2的时间序列值），然后执行二元线性回归计算系数。这里的回归模型是：
y = β₀ + β₁*x₁ + β₂*x₂
其中y是s的时间序列，x₁是s1的时间序列，x₂是s2的时间序列，β₀是截距，β₁、β₂分别是s1、s2对应的斜率。

# 定义逐像素回归函数
pixel_reg <- function(pixel_vals) {
  # 计算时间序列长度（三个栅格栈的层数一致）
  ts_length <- length(pixel_vals) / 3
  # 拆分出因变量和自变量的时间序列
  y <- pixel_vals[1:ts_length]          # s的时间序列
  x1 <- pixel_vals[(ts_length+1):(2*ts_length)]  # s1的时间序列
  x2 <- pixel_vals[(2*ts_length+1):(3*ts_length)]# s2的时间序列
  
  # 构造包含截距项的设计矩阵
  design_matrix <- cbind(1, x1, x2)
  # 用最小二乘法计算回归系数
  coeffs <- solve(t(design_matrix) %*% design_matrix) %*% t(design_matrix) %*% y
  
  # 将系数转为向量返回，顺序是：截距、s1斜率、s2斜率
  as.vector(coeffs)
}

步骤3：执行逐像素计算 #

把三个栅格栈组合成一个列表，用terra的app()函数逐像素应用我们定义的回归函数，得到结果栅格栈：

# 组合三个栅格栈
combined_stacks <- list(s, s1, s2)
# 逐像素计算回归系数
result_rasts <- app(combined_stacks, pixel_reg)

# 给结果图层命名，方便识别
names(result_rasts) <- c("intercept", "slope_s1", "slope_s2")

步骤4：查看结果 #

最后可以用plot()函数查看每个系数的空间分布：

plot(result_rasts)

补充说明 #

和你提供的单自变量示例不同，这里因为每个像素的自变量（s1、s2的时间序列）是变化的，所以没办法提前预计算(XtX)^-1，必须每个像素单独计算设计矩阵的逆。如果你的栅格数据量很大，可以考虑用terra::focal()或者并行计算来提升效率，但对于常规规模的数据，上面的代码已经足够好用啦。

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