R语言Ordinary Kriging插值求助：克里金方差异常偏大

我收集了一组观测数据，目标是用gstat包执行普通克里金（Ordinary Kriging），将观测值插值到规则网格上。目前在变异函数参数选择上遇到困难，得到的克里金方差值远高于预期：我预期网格节点的误差在2-4cm范围内，但结果显示误差为4.2-4.5米。希望能识别并纠正方法中代码或基础概念层面的问题，以下是完整代码和实验/拟合变异函数图：

完整代码 #

# Load the sf library for spatial data handling
library(sf)
library(gstat)
library(ggplot2)
library(rnaturalearth)

# Set the seed for reproducibility
set.seed(123)

world2 <- ne_coastline(scale = "medium", returnclass = "sf")
countries <- ne_countries(scale = "medium", returnclass = "sf") 

# Create a dummy data frame
dummy_data <- data.frame(
  lat = runif(1000, min = 21, max = 27),    # Random latitude values between 21 and 27
  lon = runif(1000, min = 23, max = 27), # Random longitude values between 23 and 27
  height_m = runif(1000, min = 14, max = 28), # Random height values between 14 and 28 meters
  SD_heights = runif(1000, min = 0.02, max = 0.04) # Random SD of the heights at observed locations meters
)

# Define the extent of the grid
lat_min <- 21
lat_max <- 27
long_min <- 23
long_max <- 27

# Create a sequence of latitude and longitude values
lat_values <- seq(lat_min, lat_max, by = 0.5)
long_values <- seq(long_min, long_max, by = 0.5)

# Create the grid points
grid_point <- expand.grid(lat = lat_values, lon = long_values)

# plot 
ggplot() +
  geom_sf(data = world2, size = 1, colour = "black") +
  geom_sf(data= countries,  color= "#E5E5E5", fill = "#E5E5E5") +
  coord_sf(xlim = c(23, 27), ylim = c(21, 27)) +
  geom_point(data= grid_point, aes(x = lon, y = lat), color ="red")+ 
  geom_point(data= dummy_data, aes(x = lon, y = lat))


# Ordinary Kriging
# Compute the Empirical Variogram:
# Define the original CRS (Assuming WGS 84)
original_crs <- st_crs("+proj=longlat +datum=WGS84")
# Create a data frame with the original latitude and longitude data
data_df <- dummy_data

# Convert the data frame to a spatial data frame (sf object) with the original CRS
data_sf <- st_as_sf(data_df, coords = c("lon", "lat"), crs = original_crs)

# Define the target projection (e.g., UTM Zone 33N)
target_crs <- st_crs("+proj=utm +zone=33 +datum=WGS84 +units=m")

# Project the data to the target CRS
data_sf_projected <- st_transform(data_sf, target_crs)
coor_m <- st_coordinates(data_sf_projected) %>% as.data.frame()

# Detrend the data using linear regression
# Create a data frame with the heights and coordinates
data_combined <- data.frame(
  h = data_df$height_m,
  lon = coor_m$X,
  lat = coor_m$Y
)
# Perform linear regression
lm_model <- lm(h ~ lon + lat, data = data_combined)

# Extract the coefficients a, b, and c
a <- coef(lm_model)[1]
b <- coef(lm_model)[2]
c <- coef(lm_model)[3]

detrended_h <- resid(lm_model)

# Add the detrended variable to the spatial data frame
data_sf_projected$detrended_h <- detrended_h

# Grid points for interpolation in the target CRS
# Create a data frame for the grid points in the target CRS
interp_sf <- st_as_sf(grid_point, coords = c("lon", "lat"), crs = original_crs)

# Reproject interp_sf to the target CRS
interp_sf_reprojected <- st_transform(interp_sf, crs = target_crs)

# Set a smaller cutoff value for the empirical variogram (adjust as needed)
cutoff_value <- 300000
search_radius <- 20000

# Compute the empirical variogram with the new cutoff
variogram_model_detrended <- variogram(detrended_h ~ 1, data = data_sf_projected, width = search_radius,  cutoff = cutoff_value)
plot(variogram_model_detrended)

# Fit a variogram model to the empirical variogram
variogram_fit_detrended <- fit.variogram(variogram_model_detrended, model = vgm(psill = 17.3, model = "Exp", range = 250000, nugget = 16.2))


# Plot the empirical variogram
plot(variogram_model_detrended, plot.type = "cloud", main = "Empirical and Fitted Variograms")


# Create a data frame of the empirical variogram
empirical_variogram_df <- data.frame(dist = variogram_model_detrended$dist,
                                     gamma = variogram_model_detrended$gamma)


preds = variogramLine(variogram_fit_detrended, maxdist = max(variogram_model_detrended$dist))
head(preds)


# Create the ggplot plot
ggplot() +
  geom_point(data = empirical_variogram_df, aes(x = dist/1e3, y = gamma, color = "Empirical model"), alpha = 0.6) +
  geom_line(data = preds, aes(x = dist/1e3, y = gamma, color ="Analytical model"))+
  labs(title = "Empirical and Fitted Variograms",
       x = "Distance [km]",
       y = "Semivariance") +
  theme_bw()+
  scale_color_manual(name= "", breaks=c('Empirical model', 'Analytical model'),
                     values=c('Empirical model'='red', 'Analytical model'='blue'))+
  guides(color = guide_legend(override.aes = list(shape = c(16, NA), 
                                                  linetype = c(NA, 1)))) +
  theme(#plot.title = element_text(color = "black", face="bold", size = 24),
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    axis.title = element_text(color = "black", size = 11,
                              margin = margin(t = 0, r = 20, b = 0, l = 0)),
    axis.text = element_text(color="black", size = 10),
    # panel.background = element_blank(),
    # panel.grid = element_blank(),
    legend.title = element_text(color = "black",face="bold",  size = 11),
    legend.title.align = 0.5 ,
    legend.text = element_text(color = "black", size = 10), 
    legend.key.height = unit(0.7, 'cm'), 
    #legend.key.width = unit(2.5, 'cm'),
    legend.position = "bottom", 
    aspect.ratio = 1/2)

# Perform kriging

kriged_detrended <- krige(detrended_h ~ 1, locations = data_sf_projected, newdata = interp_sf_reprojected, vgm(psill = 17.3, model = "Exp", range = 250000, nugget = 16.2))
coor_kriged_detrended <- as.data.frame(st_coordinates(kriged_detrended))
kriged_detrended_df <- data.frame(
  h = kriged_detrended$var1.pred,
  h_var = kriged_detrended$var1.var,
  lon = coor_kriged_detrended$X,
  lat = coor_kriged_detrended$Y)

kriged_predicted_h <- kriged_detrended_df$h + (coef(lm_model)[1] + coef(lm_model)[2]* kriged_detrended_df$lon+ coef(lm_model)[3] * kriged_detrended_df$lat)
head(kriged_predicted_h)

summary(sqrt(kriged_detrended_df$h_var))

问题排查与修正建议 #

1. 模拟数据无空间相关性（核心问题） #

你的dummy_data用runif生成完全随机的height_m值，数据不存在任何空间自相关性。这种情况下变异函数的块金值接近基台值，克里金无法通过空间邻近性降低预测误差，自然会产生大的方差。

• 修正：如果是测试代码，生成带空间相关性的模拟数据：

  # 生成具有空间相关性的模拟数据
  set.seed(123)
  sim_grid <- expand.grid(x = seq(23,27,0.1), y = seq(21,27,0.1))
  sim_grid_sf <- st_as_sf(sim_grid, coords = c("x","y"), crs = original_crs)
  sim_grid_proj <- st_transform(sim_grid_sf, target_crs)
  # 生成符合Exp模型的空间数据
  sim_data <- gstat::krige(z~1, locations = sim_grid_proj, newdata = sim_grid_proj,
                           model = vgm(psill=4, model="Exp", range=50000, nugget=0.0016))
  # 抽取1000个样本作为观测数据
  dummy_data <- sim_data[sample(nrow(sim_data),1000),]
  dummy_data$lat <- st_coordinates(st_transform(dummy_data, original_crs))[,2]
  dummy_data$lon <- st_coordinates(st_transform(dummy_data, original_crs))[,1]
  dummy_data$height_m <- dummy_data$var1.pred
  dummy_data$SD_heights <- runif(1000, 0.02, 0.04)
  

2. 变异函数拟合不合理 #

你手动指定的变异函数参数与实验变异函数匹配度极差，块金值几乎等于基台值，完全没捕捉到空间结构。

• 修正：
  • 让fit.variogram自动尝试多种模型：

    variogram_fit_detrended <- fit.variogram(variogram_model_detrended, 
                                             model = vgm(c("Exp", "Sph", "Gau")),
                                             fit.method = 2)
    

  • 调整变异函数计算参数：研究区域UTM下宽度约400km，将width设为5000-10000（5-10km），cutoff设为200000（200km，区域对角线一半），让分组更合理。

3. 去趋势步骤冗余 #

你的随机数据不存在空间趋势，去趋势反而会引入额外误差。

• 修正：先通过Moran's I检验或空间可视化判断数据是否存在趋势，无趋势则直接对原始height_m做普通克里金。

4. 克里金方差概念混淆 #

你混淆了观测误差（SD_heights，2-4cm）和克里金预测误差：

• 观测误差是测量时的随机误差，属于块金值的一部分；
• 克里金方差是插值的空间不确定性，单位与原始数据一致（米），开根号后为米级误差。
• 修正：若要将观测误差纳入模型，可将块金值设为mean(SD_heights)^2，或使用协同克里金结合SD_heights作为辅助变量。

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