在R中使用高斯滤波提取指定半径内像素值的均值与标准差

嘿，这个需求在栅格数据分析里挺常见的，我来给你一步步拆解怎么实现，用现在主流的terra包（比旧的raster包效率更高）就能轻松搞定：

第一步：把矩阵转换成栅格对象（或读取栅格文件） #

如果你已经有了内存里的矩阵，先把它转成terra的SpatRaster对象；如果是从栅格文件（比如TIFF）读取，直接用rast()函数就行：

# 加载terra包（如果没安装先运行install.packages("terra")）
library(terra)

# 情况1：已有矩阵my_matrix，转成栅格
my_matrix <- matrix(rnorm(100*100), nrow=100) # 示例随机矩阵
r <- rast(my_matrix)

# 情况2：从文件读取栅格
# r <- rast("your_raster_file.tif")

第二步：生成高斯权重窗口 #

我们需要根据你指定的半径x生成对应的高斯核（窗口）。高斯核的大小是2*x+1（确保中心是目标像素），然后计算每个窗口位置的权重并归一化（这样计算的是加权均值）：

# 自定义生成高斯核的函数
create_gaussian_kernel <- function(radius, sigma = radius/2) {
  size <- 2 * radius + 1
  # 生成窗口内的坐标网格
  coords <- expand.grid(x = 1:size, y = 1:size)
  center <- radius + 1
  # 计算每个点到中心的距离平方
  dist_sq <- (coords$x - center)^2 + (coords$y - center)^2
  # 计算高斯权重
  kernel <- exp(-dist_sq / (2 * sigma^2))
  # 归一化，让权重总和为1
  kernel <- kernel / sum(kernel)
  # 转成矩阵格式
  matrix(kernel, nrow = size, ncol = size)
}

# 示例：生成半径为3的高斯核（7x7窗口）
gaussian_kernel <- create_gaussian_kernel(radius = 3)

这里的sigma参数可以调整，通常设为radius/2或radius/3，控制高斯权重的衰减速度，你可以根据需求微调。

第三步：计算加权均值与标准差矩阵 #

用terra的focal()函数，对每个像素应用高斯窗口，分别计算加权均值和加权标准差：

1. 加权均值矩阵 #

# 计算加权均值，因为核已经归一化，直接用加权求和即可
mean_raster <- focal(r, w = gaussian_kernel, 
                     fun = function(x) weighted.mean(x, w = gaussian_kernel, na.rm = TRUE))

# 转成矩阵格式（如果需要）
mean_matrix <- as.matrix(mean_raster)

2. 加权标准差矩阵 #

标准差需要用加权方差推导，我们先写一个自定义的加权标准差函数：

# 自定义加权标准差函数（处理NA值）
weighted_sd <- function(x, w) {
  # 移除NA的像素和对应权重
  na_mask <- !is.na(x)
  x_clean <- x[na_mask]
  w_clean <- w[na_mask]
  
  if (length(x_clean) == 0) return(NA)
  
  # 计算加权均值
  mu <- sum(x_clean * w_clean) / sum(w_clean)
  # 计算加权方差（无偏估计）
  var_weighted <- sum(w_clean * (x_clean - mu)^2) / (sum(w_clean) - sum(w_clean^2)/sum(w_clean))
  # 返回标准差
  sqrt(var_weighted)
}

# 计算标准差栅格
sd_raster <- focal(r, w = gaussian_kernel, 
                   fun = function(x) weighted_sd(x, w = gaussian_kernel))

# 转成矩阵格式
sd_matrix <- as.matrix(sd_raster)

一些注意事项 #

• 如果你的栅格有大量NA值，确保na.rm=TRUE（或者在自定义函数里处理），避免结果全是NA
• 窗口半径x不要太大，否则计算速度会变慢，同时要确保窗口大小不超过栅格的尺寸
• 如果还在用旧的raster包，操作逻辑几乎一致，只是把rast()换成raster()，focal()用法相同

内容的提问来源于stack exchange，提问作者CaffèSospeso