做生态调查和物种保护的人十个里有八个会被问同一个问题这个物种到底分布在哪儿回答这个问题最常用的定量工具就是物种分布模型简称SDM。它用物种出现点的记录叠加温度、降水、地形这类环境变量在环境空间中拟合出物种的“偏好范围”再映射回地理空间形成一张潜在分布概率图。我之前在做一个鸟类栖息地评估项目时就是把三年野外调查整理出来的出现点全部交给R语言来跑SDM最后输出的适宜栖息地图直接用来参与了保护区调整的讨论。这条路走到今天R语言几乎是绕不开的主战场因为从数据清洗、变量筛选、模型训练到结果可视化生态学的主流工具链基本都长在R里。即便你之前只用Excel整理过物种记录只要按这篇文章的步骤把流程走一遍也能跑出一套像模像样的分布模型。文章里我会用一个完整实例串起SDM的核心流程包括数据准备、环境变量处理、建模实现、精度评估以及我亲身踩过的那些坑。1. 物种分布模型到底在解决什么问题1.1 从“在哪见过”到“可能在哪”的底层逻辑物种分布模型的核心思想并不复杂本质上就是做一次“环境偏好拟合”。你有一堆物种出现点的经纬度坐标这些点会带出一组对应的环境条件比如年均温20度、年降水量1500毫米、海拔800米。物种不会随机出现在所有地方它出现的位置一定对应着它能耐受的环境组合。模型要做的就是根据这些“出现点的环境特征”去估算环境空间里的哪些区域是它喜欢的、哪些是不喜欢的。打个比方你去相亲网站填了一堆心仪对象的条件系统根据你的“出现记录”给你推荐潜在匹配对象。物种分布模型干的就是这个活它把你物种的“择偶条件”学到手然后在整个地图上扫描哪里还有符合条件的对象。需要注意的是这里有个重要前提模型只描述“环境适宜性”不等于物种真的会出现在所有适宜区。历史上有没有过去、能不能扩散到那里、有没有天敌竞争都会影响实际分布。所以SDM的预测结果严格来讲叫“潜在分布区”或“生境适宜性”不是“实际分布区”。1.2 为什么我把SDM放进了R语言里做生态学的老底子就是统计学R语言从出生开始就长在统计和数据分析的土壤里所以生态学三大件——数据整理、空间分析、统计建模R语言全都覆盖了。具体到SDMR的优势就更明显了。一个原因是生态学家长期的积累形成了完整的包生态。dismo包里有物种分布建模的通用框架biomod2可以一口气跑十几二十个算法再集成maxnet直接在R里实现最大熵模型raster和terra负责栅格数据处理。这些包不是孤立的它们之间可以无缝衔接你用terra读一张栅格喂给maxnet训练模型再用terra把预测结果画成图全程不用切换软件。另一个原因是可复现性。野外调查的项目一做就是好几年中间人员流动频繁如果你靠Excel点鼠标操作数据一变就要全部重来。用R脚本写下来的流程换了人也能照样跑参数怎么设的、数据怎么清洗的都有迹可循。这一点对长期生态监测项目来说价值远远大于任何花哨的界面工具。1.3 一次完整SDM分析的标准流程如果只用一个流程图概括SDM项目大概是这样第一步整理物种出现点数据并做空间清洗第二步准备研究区域的环境变量栅格第三步把出现点和环境变量关联起来构建模型的训练数据第四步选择合适的算法训练模型第五步把模型预测到整个研究区域生成适宜性分布图第六步评估模型精度并确定适宜/不适宜的阈值。听起来不复杂但每一步都有讲究。数据清洗不彻底模型就会学到坐标错误或重复记录的噪音环境变量不筛选高度相关的变量会干扰参数估计背景点的选择方式会直接影响模型结果甚至比选择什么算法还重要。这篇文章后面几个部分我就是按这个流程来展开的。2. 数据准备没有好数据模型全是空中楼阁2.1 出现点数据怎么来、怎么洗出现点数据是SDM的燃料。最常见的来源一是自己的野外调查记录二是公共数据库比如全球生物多样性信息机构GBIF或者中国植物志、动物志相关数据库。GBIF的数据量很大但质量参差不齐直接下载下来用往往会踩坑。不管数据从哪来第一件事都是清洗。我自己的处理顺序大致是删除坐标缺失和坐标格式异常比如经纬度超出合理范围的记录删除明显错误坐标比如落在海里、落在研究区边界之外的去除重复记录尤其是同一格网单元内的重复点采集年份过老的记录建议单独评估因为土地利用和气候已经变了如果有条件优先保留有凭证标本或专家鉴定的记录。重复记录的问题很多人不在意但它会实实在在地扭曲模型。如果一个地方被反复调查出现了50条记录模型会误以为这是物种特别偏好的区域无形中放大了采样偏差。我一般会用空间过滤的办法比如按1公里网格去重保证模型输入的空间独立性。用R实现去重并不难先把出现点转成terra的SpatVector格式再基于栅格单元格去重。一个大原则是出现点必须是“可靠的、有代表性的”但你很难做到完美所以实操中要诚实地记录数据限制并在文章里说明。2.2 环境变量怎么选、怎么避免“变量打架”SDM的环境变量通常分成几类气候变量温度、降水、极端气候、地形变量海拔、坡度、坡向、土壤变量、植被变量以及人类活动影响变量。最常用的气候数据是WorldClim的19个生物气候变量代号bio1到bio19它们是把月均温、月降水等原始数据加工成对生物有意义的气候指标比如bio1是年均温、bio5最热月最高温、bio12年降水量、bio15降水季节性等。WorldClim数据可以在线下载也可以从官网手动下载后本地读取。R语言里有geodata包可以直接拉取但如果你网络环境一般手动下载到本地再导入反而更稳。建议用10分分辨率约18公里先跑通流程不要一上来就下载30秒分辨率的全球数据文件大不说建模时内存很容易爆掉。环境变量选择的第一个原则是“生态相关性”你选的变量要能解释这个物种为什么分布在那里。比如高山植物温度季节性、极端低温和地形就很重要对湿地水鸟而言降水和距水体距离可能排在前面。第二个原则是“避免共线性”。我见过不少案例一口气把19个bio变量全塞进模型结果是变量间高度相关模型参数很不稳定解释起来也说不清楚。常用的做法是两阶段筛选先算变量两两之间的Pearson相关系数如果绝对值大于0.7或0.8就保留一个然后对保留下来的变量再用方差膨胀因子VIF检查VIF大于10就要处理。R里有现成的函数可以做也可以自己写个循环都不复杂。2.3 研究区域边界怎么定研究区域的划定是个容易被低估的问题。你要预测的地方如果太小模型的背景点只反映了局部环境物种的完整环境偏好学不到如果预测范围太大比如拿中国全境的背景点预测一个只分布在西南山区的物种又会引入大量“模型假设上不可能存在”的过度缺失。实操上我建议研究区域最好是物种实际可达的范围可以按生态区、生物地理区来定或者结合物种的扩散能力和历史分布记录画一个合理的缓冲区。这样做既符合物种分布建模的“可达性”理论假设也能让模型的背景点更有意义。3. 用R实现一个完整的SDM流程3.1 准备环境和加载工具包为了照顾刚入门的朋友我先说明一下开发环境。我本地的配置是R 4.3以上RStudio桌面版。下面这几个包是我做SDM几乎必装的terra栅格数据处理新版R里慢慢取代raster接口更简洁速度也快geodata下载WorldClim和行政边界数据maxnet纯R实现的最大熵模型不用调Javadismo老牌SDM工具包提供很多辅助函数pROC画ROC曲线用corrplot可视化变量相关性。安装和加载代码install.packages(c(terra, geodata, maxnet, dismo, pROC, corrplot)) library(terra) library(geodata) library(maxnet) library(dismo) library(pROC) library(corrplot)如果是第一次安装这些包可能会编译较长时间建议保持网络畅通。macOS的话如果编译报错大概率要装一下Xcode Command Line ToolsWindows则要确保Rtools装好。这些环境问题属于“人人都可能遇到、但解决后一劳永逸”的典型坑。3.2 整理出现点数据这里我以自己跑过的一个项目为例结构完全脱敏但流程相同。野外调查记录的CSV长这样species,lon,lat,year Target_species,102.5132,24.8634,2021 Target_species,102.6871,24.9352,2021 Target_species,102.8423,24.6287,2022 ...读入数据后我会做一个空间点检查和去重obs - read.csv(occurrence.csv) head(obs) obs - obs[obs$species Target_species, ] obs - obs[!is.na(obs$lon) !is.na(obs$lat), ] pts - vect(obs, geom c(lon, lat), crs EPSG:4326) summary(pts)这里要特别提醒经纬度坐标统一用WGS84EPSG:4326是对的但后续建模时不一定需要投影坐标系栅格数据本身也是经纬度网格时直接跑即可。如果你要算距离、面积再考虑投影到适合研究区域的投影坐标系。去重可以用一个简单办法以环境栅格为参考按栅格单元格去重。等加载了环境变量再执行这个去重操作会更合适。3.3 获取并准备环境变量我用geodata包下载WorldClim的19个生物气候变量分辨率选10分速度快、适合概念验证bioclim_global - worldclim_global(var bio, res 10, path data/)如果你手动下载WorldClim数据通常拿到的是tif文件直接terra::rast读取即可。这一步结束后你会得到一个SpatRaster对象里面是bio1到bio19共19个图层。然后裁剪到研究范围并统一缺失值处理和变量名study_area - ext(c(98, 108, 22, 30)) # 假定研究区范围 clim_crop - crop(bioclim_global, study_area) # 给图层命名 names(clim_crop) - paste0(bio, 1:19)实际项目中我的研究区域是这样界定的先做物种分布点的缓冲区再叠加省级生态行政边界取交集作为研究区。如果你是做国家公园、保护区等特定范围直接用对应边界裁剪就行。不要贪大也不要太小合理即可。3.4 变量筛选把相关性高的变量请出去环境变量准备好以后我先做相关性检验这是建模前最有价值的一步。直接看代码set.seed(123) bg_sample - spatSample(clim_crop, size 10000, method random, na.rm TRUE) bg_sample - bg_sample[complete.cases(bg_sample), ] cor_matrix - cor(bg_sample, method pearson) corrplot(cor_matrix, method color, type upper, tl.cex 0.7)这段代码随机抽取了10000个栅格点来代表研究区的环境背景。如果corrplot图里有大块深色就说明某些bio变量高度相关。我通常会给模型选4到8个变量优先选择生态意义明晰的那些。比如我那个项目需要关注温度季节性就保留了bio4同时把和它相关度很高的bio1去掉需要反映极端低温就保留bio6和bio11。如果你愿意再严谨一些可以用VIF辅助判断。下面对候选变量做一次VIF检验的简化思路# 对候选变量计算VIF library(usdm) candidate_layers - clim_crop[[c(bio2, bio4, bio6, bio12, bio15, bio18)]] vif_result - vif(candidate_layers) print(vif_result)VIF的经验阈值是小于10严格一点会要求小于5。如果某个变量VIF超标就删掉它重新计算直到剩下的变量都能满足要求。3.5 提取出现点的环境值构造建模数据集接下来把出现点对应的环境值提取出来再生成“背景点”——也就是研究区里随机的“非出现”或者说“可用环境”样本。背景点在最大熵模型里非常重要它充当的是物种没被记录到、但环境可用的“参照系”。presence_vals - extract(clim_crop, pts, ID FALSE) presence_vals - presence_vals[complete.cases(presence_vals), ] set.seed(123) bg - spatSample(clim_crop, size 10000, method random, na.rm TRUE) bg - bg[complete.cases(bg), ] # 构造建模数据 p - c(rep(1, nrow(presence_vals)), rep(0, nrow(bg))) env_data - rbind(presence_vals, bg)背景点数量一般取出现点数量的10倍以上甚至100倍都常见。但一味增加背景点不会无限提升模型还会拖慢计算。10000个背景点对大多数项目已经够用我一般在出现点少的时候用10000点多的时候用20000左右。这里有个容易犯错的地方提取出现点环境值时如果某个出现点所在位置有环境变量缺失值比如水体被掩膜掉了这个点会被剔除。建议剔除前打印一下删了多少做到心里有数。如果删掉太多要检查是不是数据配准出了问题。3.6 用maxnet实现最大熵模型MaxEnt是SDM里最常用的算法之一它的特点是小样本下表现稳定结果相对可靠。历史上它是以Java程序运行的dismo包调用起来经常要配置路径很麻烦。maxnet包则把核心算法重写成了纯R版本直接丢掉Java依赖调用方式也更清晰。library(maxnet) p - c(rep(1, nrow(presence_vals)), rep(0, nrow(bg))) env_data - rbind(presence_vals, bg) model_maxnet - maxnet( p p, variables env_data, f maxnet.formula(p, env_data, classes lq) )这里的f是特征函数可以根据数据量选择不同的特征组合。默认的“lq”表示线性linear和二次quadratic特征不容易过拟合如果样本量大可以试试“lqph”加上乘积和铰链特征。MaxEnt的老用户对feature classes和正则化倍乘系数应该不陌生。这些参数本质上控制着模型复杂度正则化系数越大模型越平滑越不容易过拟合但太大了又会欠拟合丢失真实信号。实践经验是出现点少的时候用简单特征组合。模型训练完成后用summary看一下变量贡献。maxnet对象里有一个varmax属性直观了解哪些变量在模型中起了主要作用。3.7 把模型预测到整个研究区预测就是把训练好的模型“扫描”每个栅格单元的环境值给出一个适宜性概率。pred_rast - predict(clim_crop, model_maxnet, type cloglog) plot(pred_rast, main MaxEnt 物种分布预测)注意这里的type cloglog输出的是一个0到1之间的连续值表示相对栖息地适宜性不是严格意义上的出现概率。为什么要用cloglog而不用原始输出因为cloglog变换能更好地线性化预测结果并且与出现点密度存在数学上的关联在解释上更接近“区域内有该物种存在的概率估计”。在论文里我会把这种图标注为“生境适宜性指数”避免被误解成确切的分布概率。预测结果出来后我习惯把适宜性图按阈值二分类或者按低、中、高适宜性分三档在保护区规划时更好用。阈值怎么选下一节详细说。4. 模型评估与阈值选择别只看AUC4.1 训练集评估和交叉验证要分开看模型训练完后第一件要做的事是评估精度。最经典的指标是AUCROC曲线下面积。一个容易误用的地方很多人直接拿训练数据去算AUC结果高达0.95高兴得不得了。其实这只能证明模型“记住了”训练数据不能证明它预测能力好。更科学的做法是把数据划分成训练集和测试集或者做交叉验证。我常用的是“k折空间交叉验证”的思路而不是普通的随机划分因为物种出现点存在空间自相关临近的点环境相似随机划分会让测试集和训练集过度相似导致AUC虚高。一个简单的实现方式以5折为例set.seed(42) folds - kfold(pts, k 5) auc_vals - numeric(5) for (i in 1:5) { train_idx - which(folds ! i) test_idx - which(folds i) train_p - rep(1, length(train_idx)) bg_idx - sample(1:nrow(bg), 5000) train_data - rbind(presence_vals[train_idx, ], bg[bg_idx, ]) train_p_all - c(rep(1, length(train_idx)), rep(0, length(bg_idx))) fit - maxnet(train_p_all, train_data, f maxnet.formula(train_p_all, train_data, classes lq)) test_pres - extract(clim_crop, pts[test_idx, ], ID FALSE) test_bg - bg[sample(1:nrow(bg), 5000), ] test_data - rbind(test_pres, test_bg) test_labels - c(rep(1, nrow(test_pres)), rep(0, nrow(test_bg))) pred_vals - predict(fit, test_data, type cloglog) roc_obj - roc(test_labels, pred_vals, quiet TRUE) auc_vals[i] - as.numeric(auc(roc_obj)) } mean(auc_vals)AUC的平均值大于0.8算是可以接受能在论文里说“模型判别力较好”。这时我对模型的信心才会建立起来。4.2 TSS与阈值选择AUC是对所有可能阈值下模型表现的综合评价但在实际管理中我们需要一个具体的切割点预测值超过多少算“适宜区”。这时候就要算TSS真技能统计量TSS 敏感度 特异性 - 1取值范围-1到1越高越好。选择阈值的办法很多常用的是最大化TSS的阈值maxTSS也有人用“出现概率的10百分位”percentile 10作为阈值表示把90%的出现点都保留在“适宜区”内适合偏保守的预测。我自己的习惯是如果目标是为保护规划划定最小核心区我取TSS最大阈值如果是要找潜在扩散区域需要看得更宽泛一些我取出现点预测值10百分位阈值。计算TSS和阈值可以用一个简洁的函数find_threshold - function(obs_pres, obs_bg, thresholds seq(0.01, 0.99, by 0.01)) { tss_max - -999 t_opt - 0.5 for (t in thresholds) { tp - mean(obs_pres t) fp - mean(obs_bg t) tss - tp (1 - fp) - 1 if (tss tss_max) { tss_max - tss t_opt - t } } return(c(threshold t_opt, tss_max tss_max)) }在建模流程中先用训练集拟合再用独立的验证集计算预测值最后应用这个函数得到阈值。把阈值应用到预测栅格上就能得到一张二分类的“适宜/不适宜”图。4.3 结果可视化别再“红配绿”结果图是给合作方和评审专家看的这部分的专业度直接影响项目沟通效率。terra包画图功能不错但默认配色经常让我想吐槽。我一般用自定义配色比如从浅黄到棕褐的渐变或者使用hcl.colorsplot(pred_rast, col hcl.colors(100, palette Terrain), main 目标物种生境适宜性)如果要出符合期刊要求的图可以直接导出更高分辨率的tiff再用AI或Inkscape做最后排版。R的图形设备能输出300dpi以上的图完全够用。还有一个细节同一项目里的图范围、配色、图例要保持一致这样几期数据放一起才有可比性。5. 常见问题与排查技巧实录5.1 共线性问题静悄悄毁掉模型稳定性的元凶我早期做SDM时试过不筛选变量直接把19个bio全丢进模型。训练集AUC非常漂亮但一到空间预测就出现一些诡异的高值斑块换一批背景点结果又变了。后来做变量重要性排序时发现bio1、bio5、bio6、bio11几乎都在争抢同一段方差重要性排名每次运行都剧烈变化。排查之后才明白变量高度共线性会让模型参数极不稳定特征响应曲线也难以解释。所以现在我把“变量预筛选”放在模型训练之前的强制环节相关矩阵、VIF一定都要过一遍。这不是学术洁癖是从结果可靠性出发的刚性需求。5.2 背景点的选择比你想象的更影响结果背景点代表的是“可被使用的环境资源”不是严格意义上的“非分布区”。很多物种的实际分布受调查强度影响调查少的地方不代表物种不出现。如果背景点随机撒在全研究区可能会把“可适宜但没被发现”的区域误当作背景削弱模型区分度反过来如果只在出现点附近生成背景点模型又学不到全区域的环境特征。我实践下来的建议背景点尽量覆盖全部环境梯度与研究区范围一致数量要足够多。研究区范围原则上不小于出现点的环境空间范围否则模型预测会失真。如果做的是特定区域的精细预测可以把背景点限制在可扩散范围内但这需要很强的生态学假设支撑行文时要明确说明。5.3 小样本怎么办别硬跑复杂特征物种出现点少于30个的时候SDM的可靠性会大打折扣但很多珍稀物种偏偏就是记录极少。这种情况下我的建议是优先用简单算法和简单特征例如maxnet只保留线性特征的classes l或者用包含线性、二次的classes lq避免使用铰链、乘积这类容易过拟合的高阶特征。还可以引入“偏差校正”用出现点密度或调查强度作为偏移项减少采样偏差的干扰。模型结果解释也要更加谨慎要明确说明样本量让小尺度预测结果只能作为“方向性参考”不能直接作为精确管理单元边界。很多期刊对这类研究的接受度也正在提高前提是你把不确定性讨论充分。5.4 可复现性你的脚本就是未来的你项目结束半年后合作方突然问“能不能把之前图里阈值改一下再出一版”。如果你当时有完整脚本10分钟就能搞定如果全是手动操作基本要从头来一遍。所以我在项目一开始就建一个标准目录结构project/ ├── data/ │ ├── raw/ # 原始数据只读不写 │ ├── processed/ # 处理后的中间数据 │ └── output/ # 最终预测栅格和图表 ├── script/ # 所有R脚本 └── doc/ # 分析说明和日志脚本里用set.seed固定随机种子确保每次跑出来的结果一致所有数据导入路径用相对路径方便不同机器协作。这么做一开始会“浪费”一点时间但后面省下来的时间非常可观。这也是很多导师或团队负责人特别看重的职业习惯。5.5 预测到不同时空的陷阱外推要谨慎模型是基于当前气候条件拟合的如果你把模型直接预测到未来气候情景或完全不同的地理区域相当于拿着今天的“择偶标准”去预测几十年后的婚恋市场大概率会失真。未来气候预测不是不能做但至少要做两步一是确保未来环境变量范围不超过训练时的环境范围二是用MOP相似性分析等工具标出“外推区”在这些区域的结果可信度要打折扣。terra里有类似函数可以检测新环境的非类似区域跑一遍输出一张环境相似性图和预测图一起展示就严谨多了。写在最后的一点经验之谈物种分布模型不是终点它更像是生态学问题的一个“数据透镜”。跑通流程不难难的是每一步都带着生态学问题去思考我的数据代表什么我的背景点假设是否合理我的预测范围是否会误导决策我一直提醒自己的是R语言代码只是工具好模型依赖的是对物种、对环境的理解。如果这篇文章能让你动手把手上的出现点数据跑出第一张分布预测图那我的目的就达到了。等你有了一版完整结果再回头看那些报错信息、那些奇怪的预测斑块你会发现自己已经能跟它们“对话”了。这才是我认为做数据最有趣的地方。