如何在spatstat的kppm函数中添加矢量与栅格协变量并完成滑坡建模

1. 空间对象转换优化 #

直接用as.ppp(st_as_sf(shp_lslide))易出现窗口不匹配问题，建议先提取点坐标并绑定流域窗口：

# 提取滑坡点坐标
lslide_coords <- st_coordinates(shp_lslide)
# 转换流域为owin对象
wshed_owin <- as.owin(shp_wshed)
# 创建带窗口的ppp对象（自动过滤窗口外的点）
lslide_ppp <- ppp(lslide_coords[,1], lslide_coords[,2], window = wshed_owin)

无需额外执行lslide[wshed]子集操作，ppp创建时已完成窗口过滤。

2. 协变量处理规范 #

栅格协变量 #

• 投影对齐：所有栅格必须与ppp对象投影一致，否则会出现空间错位：

  # 以流域坐标系为基准
  target_crs <- st_crs(shp_wshed)
  # 统一所有栅格投影
  elevation <- st_transform(elevation, target_crs)
  aspect <- st_transform(aspect, target_crs)
  # 其余栅格同理
  

• 分类变量处理：土地覆盖属于分类变量，转im后需转为因子，避免被当作连续变量建模：

  landCover_im <- as.im(landCover)
  landCover_im[] <- factor(landCover_im[])
  

矢量协变量（土地覆盖、岩性） #

矢量协变量需栅格化为im对象才能用于kppm，用terra包实现更高效：

# 读取岩性矢量
litho_vect <- vect("lithology.shp")
# 以高程栅格为模板栅格化岩性，保留类别字段
litho_rast <- rasterize(litho_vect, rast(elevation), field = "FormUnitNa")
# 转为im对象并设置为因子
litho_im <- as.im(litho_rast)
litho_im[] <- factor(litho_im[])

3. 模型公式修正 #

公式中-1仅适用于分类变量（去除截距，获取每个类别的系数），连续变量无需添加。修正后的公式示例：

# 全变量模型（分类变量加-1，连续变量正常保留）
fit1 <- kppm(lslide_ppp, ~ elevation + aspect + slope + planCurva + proCurva + landCover_im - 1, "Thomas")
# 仅土地覆盖的简化模型
fit2 <- kppm(lslide_ppp, ~ landCover_im - 1, "Thomas")

4. 模型筛选与显著性分析 #

• 模型比较：用卡方检验对比不同模型的拟合优度，筛选最优模型：

  anova(fit1, fit2, test = "Chisq")
  

• 显著性识别：查看summary(fit1)输出的系数p值，筛选p<0.05的土地覆盖类别与岩性类型。

5. 滑坡敏感性预测 #

确保所有协变量覆盖整个流域窗口，直接用predict()生成敏感性图：

# 整理协变量列表
covars <- list(elevation = as.im(elevation), aspect = as.im(aspect), 
               slope = as.im(slope), planCurva = as.im(planCurva), 
               proCurva = as.im(proCurva), landCover_im = landCover_im)
# 生成敏感性预测
susceptibility <- predict(fit1, newdata = covars)
# 转为terra栅格便于后续可视化与导出
susceptibility_rast <- rast(susceptibility)
plot(susceptibility_rast, main = "流域滑坡敏感性分布图")

常见错误排查 #

• 窗口不匹配：检查所有空间对象的投影与范围是否一致
• 因子水平错误：确保分类协变量的因子水平在点数据与协变量中完全匹配
• 模型不收敛：先拟合简单模型（如仅含土地覆盖），逐步添加协变量；或调整聚类模型初始参数

内容的提问来源于stack exchange，提问作者alipin ng sahod