通过BIOMOD2进行白芷未来气候适宜性预测——承接前文的后续投影实践

Riceneeder Lv2
  2026-04-02 17:57:30 2026-04-02 17:57:30 创建   2026-04-02 19:31:19 2026-04-02 19:31:19 更新    2.6k 字  10 分钟  

上一篇文章已经完成了白芷在当前气候条件下的物种分布建模、评估与集成预测。真正有价值的下一步，不只是看“现在哪里适合”，而是继续回答“未来哪些区域还适合、适合性会如何变化”这个问题。这个步骤通常被称为未来情景投影，也是物种分布模型在气候变化研究中的核心应用之一。

这篇文章就接着前文，介绍如何在已经保存好 BIOMOD2 建模对象的前提下，直接开展后续预测。整个流程的重点很简单：先准备未来气候环境变量，再把它们与训练阶段的变量顺序和名称对齐，最后调用 BIOMOD2 的投影与集成预测接口，输出未来情景下的白芷适宜性分布图。

一、前提条件

这一部分默认你已经完成了当前气候下的建模，并保存了 biomod_objects.rds，其中至少包含：

1
2
3
4
# 读取前文保存的 BIOMOD2 模型对象
objs <- readRDS("outputs/Angelica_dahurica/biomod_objects.rds")
biomod_model_out <- objs$biomod_model_out  # 单个模型的训练结果
biomod_em <- objs$biomod_em              # 集成模型的结果

换句话说，后续预测不需要重新训练模型，只需要复用已经筛选和集成好的模型对象。这样做的好处很明显：

  • 节省大量计算时间；
  • 保证未来投影与当前模型完全一致；
  • 便于批量切换不同 SSP 情景和未来时间段。

二、准备未来气候数据

未来预测最容易出错的地方，不是模型本身，而是环境变量的准备。未来栅格文件必须满足两个条件：

  1. 变量种类与当前建模阶段一致；
  2. 变量名称和顺序与训练时严格一致。

通常建议把未来环境变量单独放进一个目录，例如：

1
2
3
4
5
6
future_env/
└── ssp126_2041-2060/
    ├── bio1.tif
    ├── bio2.tif
    ├── bio3.tif
    └── ...

如果你后面还要预测 SSP245、SSP585，或者 2021-2040、2061-2080 等时间段，可以继续按同样结构扩展目录。这样脚本只需要改一个路径，就能切换情景。

三、读取已保存对象

未来预测脚本的第一步，是直接读取前文保存的模型对象：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# 设置工作目录为项目根目录
setwd("/data2/xty/zqf/gzk")

library(terra)
library(sp)
library(dismo)
library(ggplot2)
library(viridis)
library(biomod2)
library(R.utils)
library(randomForest)
library(gbm)
library(mgcv)
library(stats)

output_dir <- file.path(getwd(), "outputs", "Angelica_dahurica")
# 定义已保存对象文件的完整路径
objs_path <- file.path(output_dir, "biomod_objects.rds")
# 检查对象文件是否存在
if (!file.exists(objs_path)) {
  stop("未找到 biomod_objects.rds，请先运行当前气候脚本到保存对象步骤。")
}
# 从磁盘读取已保存的 BIOMOD2 对象
objs <- readRDS(objs_path)
biomod_model_out <- objs$biomod_model_out  # 取出单模型训练结果
biomod_em <- objs$biomod_em                # 取出集成模型结果

这段代码的意义在于：当前模型和集成模型已经被固定下来，后续只做投影，不再重新拟合。

四、加载未来环境变量

未来栅格建议使用 terra 统一处理。它比传统的 raster 更适合大栅格数据，内存占用也更友好。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 构建未来环境变量所在目录的路径
future_env_dir <- file.path(getwd(), "future_env", "ssp126_2041-2060")

# 检查目录是否存在
if (!dir.exists(future_env_dir)) {
  stop(paste("未来环境目录不存在:", future_env_dir))
}

# 列出目录中所有 .tif 栅格文件的完整路径
future_files <- list.files(future_env_dir, pattern = "\\.tif$", full.names = TRUE)
# 检查是否找到了文件
if (length(future_files) == 0) {
  stop("未来环境目录中未找到 .tif 文件")
}

# 用 terra::rast() 读取所有 .tif 文件成为栅格堆（不读入内存，提高效率）
future_stack <- terra::rast(future_files)

如果未来文件本身没有正确命名，可以再补一层处理：

1
2
3
4
5
6
7
8
# 从已训练模型中提取训练时使用的变量名称列表
train_var_names <- biomod_model_out@expl.var.names

# 检查未来栅格是否缺少名称，若缺少则从文件名提取
if (is.null(names(future_stack)) || any(names(future_stack) == "")) {
  # 从文件路径中提取文件名（不含扩展名），作为栅格层的名称
  names(future_stack) <- tools::file_path_sans_ext(basename(future_files))
}

这一段看起来简单，但很关键。BIOMOD2 在投影时并不是只看“文件里有多少层”，而是看这些层是否能与训练时的变量一一对应。名称对不上，后面就会报错；顺序不对，结果也可能失真。

五、检查变量一致性

投影前最好先做一次显式校验：

1
2
3
4
5
6
7
8
9
# 找出在训练变量中出现、但在未来栅格中缺少的变量
missing_vars <- setdiff(train_var_names, names(future_stack))
# 如果发现缺失的变量，立即停止并报错
if (length(missing_vars) > 0) {
  stop(paste("未来环境缺少变量:", paste(missing_vars, collapse = ", ")))
}

# 按训练阶段的变量顺序重新排列未来栅格堆中的层
future_stack <- future_stack[[train_var_names]]

这里的逻辑非常直接：

  • setdiff() 用来找出训练阶段有、未来数据里没有的变量；
  • future_stack[[train_var_names]] 用来按训练时的变量顺序重新排列。

只要这一步做对，后续投影就会稳很多。

六、开展未来投影

接下来就进入核心步骤：单模型投影和集成预测。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# 定义投影的名称（用于区分不同的气候情景和时间段）
proj_name <- "future_SSP126_2041-2060"

# 第一步：用所有单个模型进行未来环境投影
biomod_proj_future <- BIOMOD_Projection(
  bm.mod = biomod_model_out,      # 已训练的单模型集合
  new.env = future_stack,         # 未来环境变量栅格堆
  proj.name = proj_name,          # 投影名称
  models.chosen = "all",          # 使用所有模型
  metric.binary = "TSS",          # 基于 TSS 指标的二元阈值
  nb.cpu = 1                      # 使用单个 CPU（根据需要可调整）
)

# 清理内存
gc()

# 第二步：基于集成模型进行集成投影预测
biomod_em_proj_future <- BIOMOD_EnsembleForecasting(
  bm.em = biomod_em,              # 已训练的集成模型
  bm.proj = biomod_proj_future,   # 单模型投影结果
  proj.name = proj_name,          # 投影名称
  models.chosen = "all",          # 使用所有集成方法
  metric.binary = "TSS",          # 基于 TSS 的阈值
  nb.cpu = 1                      # 单 CPU 计算
)

这一步的含义可以这样理解：

  • BIOMOD_Projection() 先把训练好的单个模型投到未来环境中；
  • BIOMOD_EnsembleForecasting() 再基于集成模型生成最终的综合结果。

如果你前面在集成建模阶段使用了 EMmean、EMcv 和 EMwmean，这里输出的也会是对应的一套集成投影结果。

这一步对于计算资源的要求可能比较高，尤其是当未来环境变量较多、模型数量较多时。建议在服务器上运行，并且根据实际情况调整 nb.cpu 参数以加速计算。

七、导出结果文件

完成预测后，建议把结果统一写到单独的输出目录里，避免和当前气候结果混在一起：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# 生成输出目录路径（按投影名称组织）
output_path <- file.path(output_dir, paste0("proj_", proj_name))
# 创建输出目录（若不存在）
dir.create(output_path, recursive = TRUE, showWarnings = FALSE)

# 将 terra 的 SpatRaster 对象解包成能直接写入栅格文件的格式
unpacked_val <- terra::unwrap(biomod_em_proj_future@proj.out@val)

# 导出完整的集成预测结果（包含所有层）
writeRaster(
  unpacked_val,
  filename = file.path(output_path, "Angelica_SSP126_2041-2060_ensemble.tif"),
  overwrite = TRUE
)

# 导出集成预测的第一层（通常是 EMmean，集成模型的平均值）
if (nlyr(unpacked_val) >= 1) {
  writeRaster(unpacked_val[[1]], file.path(output_path, "Angelica_SSP126_2041-2060_EMmean.tif"), overwrite = TRUE)
}
# 导出第二层（通常是 EMcv，集成模型的变异系数，反映模型间的差异）
if (nlyr(unpacked_val) >= 2) {
  writeRaster(unpacked_val[[2]], file.path(output_path, "Angelica_SSP126_2041-2060_EMcv.tif"), overwrite = TRUE)
}
# 导出第三层（通常是 EMwmean，集成模型的加权平均值）
if (nlyr(unpacked_val) >= 3) {
  writeRaster(unpacked_val[[3]], file.path(output_path, "Angelica_SSP126_2041-2060_EMwmean.tif"), overwrite = TRUE)
}

# 打印完成提示信息
message("未来气候预测完成，结果已保存至: ", output_path)

导出后，你会得到一个总的集成结果文件，以及若干拆分后的子结果。这样后续无论是做制图、阈值分类，还是和当前气候结果做差值分析，都会方便很多。

八、可选的自动提醒

如果你经常在本地或服务器上跑长时间任务，可以顺手加一个 Bark 推送，脚本跑完后直接收到通知：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# Bark 推送服务的基础 URL（需替换为自己设备的 Key）
bark_base_url <- "https://api.day.app/你自己设备的BarkKey/"
# 推送消息的标题
bark_title <- "BIOMOD2"
# 推送消息的内容体
bark_body <- paste0("未来气候预测完成，结果已保存至: ", output_path)

# 组合完整的推送 URL（将标题和内容进行 URL 编码）
bark_url <- paste0(
  bark_base_url,
  utils::URLencode(bark_title, reserved = TRUE),
  "/",
  utils::URLencode(bark_body, reserved = TRUE)
)

# 使用 curl 命令发送 HTTP 请求到 Bark 服务（不显示输出）
try({
  system(paste("curl -fsS", shQuote(bark_url)), ignore.stdout = TRUE, ignore.stderr = TRUE)
}, silent = TRUE)  # 若发送失败也不中断脚本

这部分不是必须，但对于长时间批量投影很实用。

九、后续可以怎么扩展

如果你想把这套流程继续做得更完整，下一步通常有三种扩展方向：

  1. 批量跑多个 SSP 情景，比如 SSP126、SSP245、SSP585；
  2. 比较不同时间段的变化，比如 2031-2050、2041-2060、2061-2080；
  3. 把未来结果和当前结果做差值分析，输出适宜区扩张、收缩和迁移方向。

其中第三步往往最有论文价值，因为它不仅能告诉你“未来哪里适合”，还可以进一步解释“适宜区是增加了还是减少了”。

十、总结

白芷的未来适宜性预测，本质上是把前文训练好的 BIOMOD2 模型，迁移到气候变化情景下重新评估其空间适宜度。只要确保模型对象已保存、未来环境变量与训练变量完全一致，再用 BIOMOD_Projection() 和 BIOMOD_EnsembleForecasting() 依次完成单模型和集成投影，就能得到相对可靠的未来分布结果。

AIGC声明

本文内容使用AI根据源代码生成，可能存在不准确或不完整的地方，请读者自行核实相关信息。

  • 标题: 通过BIOMOD2进行白芷未来气候适宜性预测——承接前文的后续投影实践
  • 作者: Riceneeder
  • 创建于 : 2026-04-02 17:57:30
  • 更新于 : 2026-04-02 19:31:19
  • 链接: https://gankun.cn.lu/posts/2026-04-02/
  • 版权声明: 版权所有 © Riceneeder，禁止转载。
评论
目录
通过BIOMOD2进行白芷未来气候适宜性预测——承接前文的后续投影实践
  1. 一、前提条件
  2. 二、准备未来气候数据
  3. 三、读取已保存对象
  4. 四、加载未来环境变量
  5. 五、检查变量一致性
  6. 六、开展未来投影
  7. 七、导出结果文件
  8. 八、可选的自动提醒
  9. 九、后续可以怎么扩展
  10. 十、总结