1. 物种分布模型到底在解决什么问题1.1 从一张分布点到一张预测图的思考方式转变做生态学或者保护生物学相关工作的人一定都有过这种经历手头攒了一批某个物种的GPS记录点可能是自己跑野外收的也可能是从标本数据库里扒下来的然后呢大部分人的下一步操作是打开地图软件把这些点撒上去看一眼“哦这个物种主要分布在西南山区”然后写进报告里。但如果有人问你为什么这些点出现在这里而不是出现在隔壁那个环境看起来差不多的山头或者问你气候变化以后这些点还会存在吗你就卡住了。物种分布模型Species Distribution Model以下简称SDM本质上就是把“分布点在哪里”这个问题翻译成“分布点背后的环境条件是什么”。它的核心逻辑并不复杂收集物种的出现记录提取这些记录点位置上的环境数据温度、降水、海拔、植被等等然后建立一个统计或机器学习模型去拟合“什么样的环境条件组合下这个物种会出现”最后把这个拟合结果投射回整个研究区域得到一张连续的概率分布图。我第一次用R做SDM的时候最震撼的其实不是模型的数学细节而是思考方式的转变从“记录点在地图上的位置”转向“环境变量在空间上的取值”。分布点只是采样环境条件才是真正的解释变量。这个转变一旦完成后面很多操作变量筛选、模型评价、空间预测都会顺畅很多。1.2 为什么我推荐用R而不是其他工具市面上能做SDM的工具不少MaxEnt软件是很多人入门的选择界面友好、点几下鼠标就能出图。但用过一段时间你就会发现MaxEnt的GUI版本在数据处理、批量建模、结果提取这些环节上非常不灵活。ARCGIS也能做但收费且脚本化程度低。相比之下R语言在这个领域几乎是事实标准生态学顶刊上发表的SDM相关文章绝大多数都在方法部分提到了R主要原因有三点第一R的包生态太完整了。从数据清洗dplyr、sf、环境数据下载raster、terra、geodata、模型算法maxnet、randomForest、gbm、biomod2、到结果可视化ggplot2、tmap全部贯通。一个项目从原始数据到发表级图表可以不走任何其他软件。第二可重复性。GUI操作没法记录完整的分析日志但R脚本可以。你自己的项目三个月以后回看还能知道每一步做了什么同行评审要求提供分析代码也能直接给出去。这是现代科研合作和出版的基本要求。第三灵活性和可扩展性。批量跑30个物种的模型写个循环就行。比较5种算法的表现biomod2一套搞定。把模型投影到2070年的气候情景下换个环境数据集重新predict就行。这些操作在GUI软件里要么做不到要么费半天劲。当然R的学习曲线比点鼠标要陡一些但SDM这个领域值得你付出这个学习成本。接下来我按自己实际建模的流程一步步拆解每个关键环节怎么做、为什么这么做以及哪些坑我替你踩过了。2. 建模前的数据准备最容易翻车的环节2.1 物种出现数据你收集到的只是“半个真相”先说一个很多新手会忽略的事实SDM永远存在“假缺失”的问题。你拿到的出现点代表的是“在这里发现了这个物种”但不代表“没出现在别的地方”。所以模型本质上是在回答一个不太对称的问题——“当前记录点所对应的环境条件能不能把这里和整个背景区域区分开”。理解了这一点你就明白为什么出现数据的质量直接影响模型成败。出现数据的主要来源是全球生物多样性信息网络GBIF这个数据库收录了全球几十亿条物种观测和标本记录下载很方便。但直接用GBIF原始数据进行建模几乎是灾难级的操作因为数据里的坑实在太多了坐标漂移和错误比如经纬度恰好落在海洋里、落在国界外、或者是个整数坐标比如36.0000103.0000这些大概率是地理参考不准的记录。重复记录同一个标本可能被多次数字化录入或者一条采集记录被重复上传。采样偏差这是最隐蔽也最致命的——人们倾向于在道路边、保护区边界、研究站附近采集标本所以出现点在空间上高度聚集模型会把“容易被人类采样”这个特征学进去。我常用的清洗流程先下载数据后用sftime或sf包转成空间对象检查坐标系然后过滤掉坐标缺失、坐标为0、落在海洋里的点可以用gshhg海岸线数据做反向裁剪再做空间稀疏化去重。空间稀疏化的逻辑很简单假设你的环境数据分辨率是1公里那么同一公里网格内保留一条记录就够了。这个操作可以用spThin包或者用sf包自己写一个按网格抽样的函数。注意稀疏化之后一定要看一眼保留了多少个点。一般来说用于建模的独立出现点最好不少于30个低于这个数模型容易过拟合输出结果解释起来很勉强。我见过有人拿7个点跑MaxEnt结果响应曲线完全是锯齿状根本没法用。除了坐标清洗还有一个常被忽略的操作检查物种鉴定是否可靠。GBIF里同一个物种名下可能混了亚种、近似种甚至错误鉴定记录。对于生态需求差别大的类群比如一种耐旱、一种喜湿但形态相似不分清楚就直接建模出来的模型解释性会很差。2.2 环境变量选对变量比调参更管用建模的第二根支柱是环境变量。常用的环境数据源是WorldClim提供了历史上1970-2000和未来气候情景下的19个生物气候变量BIO1-BIO19包括年均温、昼夜温差、季节性降水这些经过生物学转化的指标。为什么不用原始气象数据而是用这些衍生变量因为物种对“年均温”这种平均值不敏感但对“最暖季度降水量”“最干月降水量”这种极限和季节性指标更敏感这些变量更贴近生态位的内涵。拿到数据后筛选变量这步一定要做仔细。很多人贪全把19个变量全部丢进模型结果变量共线性严重模型参数解释混乱重要的变量被噪声掩盖。操作上我建议两条腿走路相关性分析先用raster或terra包对变量做相关矩阵分析把Pearson相关系数绝对值大于0.7的变量视为高度相关同一组内只保留你认为生态意义更重要或更能代表机制的一个。比如BIO1年均温和BIO5最暖月最高温经常相关度很高那你需要根据自己的研究问题决定保留哪个。VIF方差膨胀因子检验用usdm包的vifstep函数可以对保留后的变量组合做多重共线性检验自动迭代剔除VIF超过阈值通常取10的变量。我个人的习惯是先用相关性矩阵粗筛到6-8个变量再用VIF检验确认“最终名单”。变量数量不宜过多特别是当出现点只有几十个时每多加一个变量都在消耗模型自由度。地层上还有一个很少人提但很重要的点变量的空间分辨率要统一。如果你的海拔数据是90米的SRTM气候数据是1公里的WorldClim植被数据是500米的MODIS那么建模前必须通过重采样统一到同一个网格系统。最省事的做法是先定义你研究区域的栅格模板分辨率、范围、坐标系然后把所有变量重采样到这张模板上。R里用terra包的project和resample函数就能完成。2.3 伪缺失点没有“没找到”也要造前面提过SDM的任务区分为“出现”和“背景”这就引出一个必须处理的问题物种分布模型需要知道“不出现”的区域在哪里。但真实数据里几乎没有可靠的“确认不出现”记录除非做了系统调查并有阴性记录绝大多数时候你需要自己生成伪缺失点全称是pseudo-absence或background points。这里存在一个方法学上的派系之分presence-only方法如MaxEnt认为你只需要提供背景点背景点的含义是“环境状态在某一小格被随机采样到的机会”物种出现与否未指认presence-absence方法如广义线性模型、随机森林则需要明确的“假缺失”点它们会被当作0参与建模。我的经验是生成伪缺失点要注意三点数量要足一般是出现点数量的10到50倍。原因是环境网格中有少量不适合的格子可能未被采样若背景点太少模型对“不适宜区域”的辨识力不足。采样策略要合理最简单的是在研究区域范围内完全随机撒点。但更好的做法是需要与出现点在空间上保持一些距离比如在出现点周围划一个1公里或5公里的缓冲区缓冲区内的点不作为伪缺失点这样避免模型只学到了“离记录点近就是宜栖”的伪规律。重复抽样最好重复生成3-5套背景点分别建模后取平均或委员会投票这样降低单次随机抽样带来的偶然性。biomod2在这方面可以直接设置NbRunEval和NbRepetitions。3. 建模方案与核心代码实现3.1 建模策略单个算法还是集成模型SDM的算法库丰富到了让新手选择困难的地步广义线性模型、广义加性模型、随机森林、梯度提升机、最大熵模型MaxEnt、人工神经网络、支持向量机等等。我实际建模的体感是算法之间在历史气候投影下的结果差异远小于变量筛选和出现数据质量带来的差异。换句话说多数情况下你不需要过度纠结选哪个算法。但如果你需要发表文章或任务比较重要我的建议是用集成建模Ensemble Modeling的思路跑多个算法比如随机森林、梯度提升机、MaxEnt然后计算AUC和TSS评估每个算法的表现给表现好的算法更高权重做加权平均预测。biomod2包把这个流程封装得很好也是目前生态学文章里最常见的SDM建模框架之一。如果你只是想快速得到一张像样的预测图maxnet包MaxEnt的R实现是最省事的选择它只需要出现点和背景点的环境数据就能拟合出一个不错的模型。3.2 完整的R建模流程示例这一节我给出一个从环境数据到预测图的完整R代码流程用的数据是我造的一个模拟场景假设我们要对“某类山地植物”在西南地区做分布建模环境数据就用WorldClim的六个变量。你想跑真实数据时把读入文件路径换一下就行。第一步准备环境变量栅格library(terra) # 设置工作目录把下载好的环境变量tif放在这里 setwd(D:/sdm_demo) # 假设已经有以下6个变量tif文件 vars - c(bio1.tif, bio4.tif, bio7.tif, bio12.tif, bio15.tif, bio18.tif) env_stack - rast(vars) # 统一检查分辨率、范围和坐标系 env_stack这里有个细节如果你从WorldClim下载数据后解压得到的变量分辨率完全一致那直接用rast()读取即可。但如果你混合了不同来源的数据需要先做重采样。重采样的方法是以其中一个栅格为模板# 将env_stack中的所有变量统一到第一个变量的分辨率 env_unified - resample(env_stack, env_stack[[1]])第二步准备出现点和背景点library(sf) # 读入物种出现记录CSV里至少要有x经度和y纬度两列 occ - read.csv(species_occurrence.csv) occ_sf - st_as_sf(occ, coords c(x, y), crs 4326) # 提取每个出现点对应的环境变量值 occ_env - terra::extract(env_unified, occ_sf) occ_data - cbind(occ, occ_env) # 剔除环境变量有缺失的记录 occ_data - na.omit(occ_data)生成背景点有两种常用方式。一种是用terra包在整张栅格上随机采样# 在非NA的栅格区域随机生成10000个背景点 set.seed(123) bg_points - spatSample(env_unified, size 10000, method random, na.rm TRUE, as.points TRUE, xy TRUE) bg_env - terra::extract(env_unified, bg_points) bg_data - cbind(as.data.frame(bg_points), bg_env) bg_data - na.omit(bg_data)另一种是带空间缓冲的伪缺失点离所有出现点至少1公里# 给出现点做一个1km缓冲区 buffer_sf - st_buffer(occ_sf, dist 1000) # 在环境图层上随机采样10000点 candidate_pts - spatSample(env_unified, size 100000, method random, na.rm TRUE, as.points TRUE, xy TRUE) candidate_sf - st_as_sf(candidate_pts, coords c(x, y), crs 4326) # 去掉落在缓冲区内的点 inside - st_intersects(candidate_sf, buffer_sf, sparse FALSE) bg_final - candidate_sf[!apply(inside, 1, any), ] # 如果过滤后点太少可以增大初始采样量第三步建模这里我先演示用maxnet包跑一个快速的MaxEnt模型后面再演示biomod2集成建模。library(maxnet) # 准备模型输入出现/背景标记 环境变量矩阵 set.seed(123) bg_sample - bg_data[sample(1:nrow(bg_data), 1000), ] # 合并出现与环境数据注意列名统一 presence - rep(1, nrow(occ_data)) absence - rep(0, nrow(bg_sample)) all_data - rbind(occ_data[ , c(bio1, bio4, bio7, bio12, bio15, bio18)], bg_sample[ , c(bio1, bio4, bio7, bio12, bio15, bio18)]) labels - c(presence, absence) # 拟合模型 fit - maxnet(p labels, variables all_data, regmult 1.0) # 正则化参数默认1.0调大更平滑 # 查看变量贡献 plot(fit)maxnet的regmult参数是控制模型平滑度的默认1.0。如果出现点很少建议调大到2或3减少过拟合如果出现点很多可以适当降低到0.5。这个参数的经验性很强需要多看响应曲线来判断。第四步预测和出图# 对整个研究区域做预测 pred - predict(env_unified, fit, type cloglog) # 保存为GeoTIFF writeRaster(pred, species_distribution_prediction.tif, overwrite TRUE) # 用ggplot2画图 library(ggplot2) library(terra) library(tidyterra) ggplot() geom_spatraster(data pred) scale_fill_viridis_c() labs(fill 适生概率) theme_minimal()到这里一个最小可行的SDM流程就跑通了。3.3 升级版用biomod2做集成建模如果你目标是要写文章建议直接用biomod2。它的核心优势是自动化处理模型评价、多次重复、多算法集成并统一输出格式。下面代码是基于当前版本biomod24.2.x的用法与老版本有些差异但注释够清晰library(biomod2) # 读取环境变量和物种数据 myExpl - stack(env_unified) # biomod2的4.2版本还依赖raster对象 # 构建建模数据对象 myRespName - SpeciesA myRespXY - occ_data[ , c(x, y)] myResp - rep(1, nrow(occ_data)) # 出现点全部标记1 # 生成伪缺失点 set.seed(123) bg_coords - as.data.frame(bg_final, xy TRUE)[, c(x, y)] abs_points - sample(1:nrow(bg_coords), 10 * nrow(occ_data)) myRespXY_bg - bg_coords[abs_points, ] myResp_bg - rep(0, length(abs_points)) # 合并出现背景 myRespXY_total - rbind(myRespXY, myRespXY_bg) myResp_total - c(myResp, myResp_bg) # 构建BIOMOD_FormatingData myBiomodData - BIOMOD_FormatingData( resp.var myResp_total, expl.var myExpl, resp.xy myRespXY_total, resp.name myRespName, PA.nb.rep 3, # 重复3次伪缺失采样 PA.nb.absences 1000, # 每次1000个伪缺失点 PA.strategy random ) # 设置模型选项 myBiomodOption - BIOMOD_ModelingOptions( MAXENT.Phillips list(regmult 1), RF list(ntree 500), GBM list(n.trees 2000, interaction.depth 3) ) # 建模 myBiomodModelOut - BIOMOD_Modeling( myBiomodData, models c(RF, GBM, MAXENT.Phillips), models.options myBiomodOption, NbRunEval 5, # 5次重复 DataSplit 80, # 80%训练20%评价 Prevalence 0.5, VarImport 3, # 计算变量重要度重复3次 do.full.models TRUE, seed.val 123 ) # 模型评价 myBiomodModelEval - get_evaluations(myBiomodModelOut) AUC_values - myBiomodModelEval[AUC, Testing.data, , ,] TSS_values - myBiomodModelEval[TSS, Testing.data, , ,] # 展示评价结果 print(AUC_values) print(TSS_values)biomod2的评价逻辑非常清晰每次运行会用80%数据建模型、20%数据做测试重复5次所以你需要关注的是AUC和TSS在多次重复中的均值和波动。AUC大于0.8、TSS大于0.6这个模型基本算合格。跑完评价后做集成投影# 构建集成模型取AUC权重平均 myBiomodEM - BIOMOD_EnsembleModeling( modeling.output myBiomodModelOut, chosen.models all, em.by all, eval.metric c(TSS), eval.metric.quality.threshold c(0.6), models.eval.meth c(TSS), prob.mean TRUE, # 输出平均概率 prob.mean.weight TRUE # 按TSS值加权平均 ) # 做空间投影 myBiomodProj - BIOMOD_Projection( modeling.output myBiomodModelOut, new.env myExpl, proj.name current, selected.models all, binary.meth TSS ) # 绘制预测图结果会输出到工作目录 plot(myBiomodProj)这段代码跑完后你会在工作目录下得到当前气候条件下的集成适生区图。每一个小步骤都有对应的输出文件非常适合后续整理成论文方法附件。3.4 把模型投射到未来气候情景小心“环境外推”很多SDM项目最后都要做未来预测比如“2050年该物种的适生区会怎么变化”。这个需求听起来简单但有个基础性的问题未来气候条件可能超出建模时所用历史环境变量的取值范围。举个例子你建模时所用的年均温范围是5-20℃但未来某个高排放情景下部分地区年均温达到23℃模型在处理这个超过训练范围的新数据时就等于在“盲猜”。这种预测外推到训练数据未覆盖的区域时极不可靠学术界称之为“no-analog climate”无类比气候问题。R里可以用extrapolation相关的工具检查外推风险比如ecospat包的ecospat.mess函数library(ecospat) # 计算MESS指标识别哪些区域的环境变量值超出训练范围 mess_result - ecospat.mess( env_stack, # 当前环境变量应转换为矩阵 unique(occ_data[ , vars]) # 训练时用到的环境取值 )结果中负值越大的区域代表外推风险越高预测结果的可信度越需要打折。发表在期刊上的文章审稿人经常要求作者报告MESS或类似的“漂移”检验结果。4. 常见问题与模型调试实录4.1 采样偏差陷阱模型学的是你的采样路线我处理过一个案例某研究团队用某种两栖动物模型预测保护区适生区结果模型图上一团一团的高适生区恰好都围绕几条公路沿线。排查数据后发现GBIF上这个物种的数据大多来自沿路调查而远离公路的森林区域几乎采样空白。模型学到的是“公路附近环境适宜”自然结果变形。解决办法有三个层次最省事用空间稀疏化后再对背景点做“采样偏差矫正”——比如在生成背景点时也偏向公路附近形成一个与出现点相似的采样概率分布。这种思路在maxnet中可以直接传入“采样密度栅格”配合实现。更彻底结合系统调查数据主动补采一些距离道路远的区域但这需要野外人力现实约束大。折中只把结果解释为“当前调查偏好的条件下该物种可能的分布”在文章中明确说出采样偏差的局限性而不是斩钉截铁地画“适生区边界”。4.2 模型AUC总是0.95先别高兴新手最容易犯的毛病看到评价集AUC接近0.95甚至0.99高兴得不行觉得模型非常完美。实质上AUC过高往往意味着数据存在泄漏或空间自相关同一地点的采样记录被同时分到了训练集和测试集空间距离极近模型其实在背这些点的位置特征而不是环境特征。我在实操中必做的一个操作数据分块spatial block。不要用随机抽样切分训练和测试改成按空间位置分块比如把研究区域切成4x4的网格在每个网格内取样一部分网格整体作为训练另一部分作为测试。R里的blockCV包专门干这个library(blockCV) # 基于出现点的空间坐标划分训练/测试块 cv_scheme - spatialBlock( speciesData occ_sf, species Presence, theRange 50000, # 分块大小一般设为与研究尺度匹配的距离单位米 k 5, selection systematic ) # 看每个块的分配情况 plot(cv_scheme)用空间分块的方式得到的AUC才是更可靠的模型泛化能力指标。如果分块后AUC明显下降说明模型存在空间自相关过拟合需要平衡变量、增加正则化或降低模型复杂度。AUC阈值参考多说一句0.7-0.8还算可接受0.8-0.9是合理超过0.9平均而言已经属于“有点好到不真实”的范围尤其当出现点较少时尤其要警惕。4.3 投影时报错的排查思路我接触过很多R语言用户的SDM报错下面这个速查表覆盖了最常见的几类报错信息或核心现象可能原因解决方法Error: cannot open file ...读取tif的路径中含中文或特殊字符路径全改为英文确认文件存在用list.files()查看出现点和背景点的环境变量全是NA点坐标在栅格范围外或坐标系不一致crs()检查先用st_transform统一坐标系eval(ext)报错或包函数找不到biomod2与raster版本冲突更新或固定包版本检查library()加载顺序预测栅格大部分为极低值或全0环境变量范围不匹配或变量名冲突检查names()确认用于建模的变量和预测时的变量名完全一致出现点过少导致maxnet拟合失败数据量不足或正则化过小增大regmult到2以上减少模型变量数量spatSample产生大量NA背景点栅格中有大片NA如海洋或缺失区用mask裁剪后再采样其中最常见的隐藏bug是变量名不一致。你建模时用的数据框列名是bio1但投影时栅格图层名可能是wc2.1_30s_bio_1。predict函数遇到的变量名和训练时完全不同模型输出的结果没有意义但代码不会报错。所以每次投影前用names(env_stack)对照一遍训练数据的列名是负责任的做法。4.4 关于模型迁移和未来预测的一些更现实的建议即使你把模型做得很好迁移到未来气候情景时也仍然有两个容易被忽视的问题第一未来气候下物种有演化适应和扩散能力但模型默认“生态位保守”。也就是说模型假设物种对环境的需求在未来不变物种也不会迁到模型认为不适生的区域之外。这是一种简化但也是一种必要的简化。在写讨论时必须正面说明这个假设的存在和影响。第二选择合适的气候情景和全球气候模式GCM。不同GCM对同一情景的模拟结果差异很大一篇严谨的文章至少要跑多个GCM然后展示一致性和不确定性范围。WorldClim未来数据里提供了不同GCM、不同SSP共享社会经济路径下的选择共享路径比如SSP126、SSP245、SSP585分别对应低、中、高排放。建议至少用3个GCM与2个SSP做组合最后集成取交集或者分级结果。4.5 一个小但实用的建议多做痕迹“软检验”最后一个建议来自我踩过多次坑之后的体会模型做好之后先不要急着下结论。拿预测图和物种已知的生物学知识做一次“软检验”。你研究的物种有没有已知的分布上限预测图里适生区会不会跑到了海拔6000米的极高山地有没有已知的分布下限比如某些物种不会出现在极端干旱的荒漠区如果预测图中出现大面积超出物种已知忍耐范围的区域很可能是变量筛选或扩样偏差导致的问题需要回到数据环节排查。我个人的习惯是把预测图对照全球生物多样性信息网络上的全部有效记录包括建模用的和未用的用肉眼检查一遍已知分布点是否落在高适生区有没有大量历史记录点被模型预测成完全不适宜如果有检查是不是和出现点坐标精度、变量分辨率有关。这个过程虽然“不严谨”但非常有效新手和老手都在用只是老手不会告诉你。5. 写在最后的一点体会物种分布模型做多了以后我的最大体会是这个模型漂亮与否80%取决于数据清洗和变量筛选20%才是算法调参和代码实现。很多人一上来就急着跑模型花在研究环境变量相关性上的时间几乎没有最后得到一张很漂亮但解释不清的预测图。回到文章开头的话题——SDM的精髓就在于把“分布点”转化为“环境条件”如果你环境数据这步没有想清楚后面所有输出都是空中楼阁。我自己每次接一个新的SDM项目都会先花两到三天时间只做一件事整理出现点数据、检查数据的空间分布模式、筛选环境变量、做出变量相关性矩阵然后再开始建模。这个习惯帮我省掉了大量返工时间。如果你准备开始自己的SDM项目我也建议你从数据准备和变量理解做起R代码只要多练习几次就能流畅但判断“变量选得对不对”“数据是否足够可靠”的能力才是真正需要长期积累的本事。