1. 从两张图看穿整片林子为什么要把空间格局和群落稳定性放在一起分析做森林生态研究的人都有一个共同的痛点样方数据好不容易收回来了除了算个Shannon-Wiener指数画张柱状图好像就不知道还能干什么了。更尴尬的是如果有人问一句你这几个样方的多样性差异是真实存在的还是因为样方离得近才像的你大概率会愣住。这个问题不是刁难它是森林生态学的核心命题之一——生物多样性的空间格局到底受什么控制以及这种格局能不能维持群落的稳定。R语言恰好是解决这一连串问题的高效工具链。它能把零散分布在样方记录表里的物种多度数据、地理坐标、环境因子数据串成一条完整的分析流水线从α多样性样方内的物种丰富度到β多样性样方间的物种组成差异再到群落时间稳定性多度在时间轴上的波动程度每一步都有成熟的R包和方法支持。这不是什么高深的新范式噱头而是把生态学里一直存在但常常被分开处理的问题真正合并起来你既要描述格局也要解释格局的形成机制还要评估格局变化对群落稳定性的后果。这篇内容适合谁看正在做森林固定样地复查数据的硕士生和博士生准备用已有数据重新挖掘故事的研究人员以及想系统梳理R语言生态分析流程的入门者。我会从数据整理讲起一路走到格局分析和稳定性评估涉及的代码逻辑、参数选择、可视化呈现方式都会结合实际样地数据展开。看完你至少能回答两件事我的样地数据里多样性在空间上是怎么分布的这种分布方式对群落的稳定性意味着什么2. 数据的底子打不好后面全是空中楼阁2.1 样地数据结构长什么样才算能用先别急着跑代码。森林样地的原始数据通常长这样每一行是一条样方-物种-多度记录配套一个样方坐标经纬度或UTM坐标、一个样方大小比如20m×20m、若干环境因子海拔、坡度、坡向、土壤pH、林冠开度等。如果你想做时间稳定性分析还需要同一批样方在不同年份的重复调查数据至少两期三期以上更理想。我见过太多人拿到数据第一件事就是算多样性指数结果算到一半发现样方ID有重复、物种名有同物异名、多度单位有的记株数有的记盖度。R语言里做数据处理的第一步永远不是分析而是把这张表清洗成每一行唯一、每一列明确的标准长数据格式。我的建议是先用tidyr和dplyr做一次全量体检library(tidyverse) # 假设你的原始数据叫 raw_data.csv df - read.csv(raw_data.csv, stringsAsFactors FALSE) # 查看缺失值、重复样方、异常多度 summary(df) df %% count(plot_id) %% filter(n 1) # 检查样方是否重复 df %% filter(abundance 0) # 多度是否出现0或负值这里有个很容易忽略的细节多度为0的记录在长数据里应不应该保留我的习惯是分析之前剔除但保留一份原始备份。因为有些多样性指数比如Chao1需要知道样方里未被观测到的物种的估计这个信息是从仅有1株的物种数和仅有2株的物种数推算的0记录不会参与计算但会影响你对数据覆盖率的判断。先把原始数据存档再在分析副本里操作这是R项目的第一条铁律。2.2 坐标和其他环境因子的预处理不标准化就等着结果骗你样方坐标到手之后先确认坐标系是否统一。森林样地数据经常混着GPS的经纬度WGS84和RTK测的UTM坐标混在一起算距离会得到荒谬的结果。R语言里处理坐标系常用sf包但如果你只是算样方间的地理距离用geosphere包的distm函数就够了library(geosphere) # coords 是 data.frame包含经纬度两列 d_geo - distm(coords, fun distGeo) # 返回距离矩阵单位米环境因子方面别把所有变量一股脑全塞进模型。先做相关性筛查海拔和温度、坡度和土壤含水量经常高度相关保留生态学意义最强的一个即可。连续变量在进入多元模型之前要标准化z-score否则量纲大的变量会主导距离计算和回归系数。这里用vegan::decostand(x, method standardize)是最省事的做法。2.3 空间自相关的第一印象先画一张样方分布图在动手算任何指数之前先把样方位置画出来用颜色或气泡大小表示物种丰富度。这一步治好了我无数次拿着数据硬跑模型的毛病。样方如果集中在山谷、沿沟系分布那你的空间格局很可能是取样设计造成的伪格局如果样方覆盖了整个坡面那分析结果才有生态解释价值。用ggplot2加ggspatial做底图是最快的library(ggplot2) ggplot(coords, aes(x easting, y northing, size richness)) geom_point(alpha 0.7) coord_equal() theme_minimal()如果样方数量比较多超过30个还可以叠加一个简易的插值表面用gstat包的idw做反距离加权插值快速看看多样性高值区在空间上是连片还是离散。这一步的输出不一定会进论文但它决定了你后面选哪种空间分析方法。3. 多样性空间格局的分析主线从α到β再到γ的拆解逻辑3.1 α多样性的多面性丰富度、Shannon、Simpson到底该怎么选α多样性是每个样方内部的物种多样性最常用的三个指数各有脾气。物种丰富度Species richness就是数数直观但容易被稀有种带偏Shannon指数对常见种敏感数值越大说明不确定性越高也就是多样性越高Simpson指数对优势种敏感反映的是随机抽两个个体它们是不同物种的概率。没有哪个指数绝对正确关键看你问什么问题。如果你关心的是群落对干扰的响应Simpson指数往往更稳因为它主要受优势种影响而优势种通常在干扰下最先变化如果你关心的是生境异质性对物种共存的影响物种丰富度更直接。我的建议是三者都算放进一个数据框里后续分析各取所需。vegan的diversity函数一条命令搞定library(vegan) # comm 是样方-物种多度矩阵行是样方列是物种 rich - specnumber(comm) shannon - diversity(comm, index shannon) simpson - diversity(comm, index simpson) alpha_df - data.frame(plot_id rownames(comm), rich, shannon, simpson)还有个容易被忽略的步骤稀释曲线rarefaction。如果不同样方的总多度差异很大直接比较物种丰富度会不公平——你采到100个个体的样方当然比只采到20个个体的样方物种更多。用vegan::rrarefy把所有样方稀释到相同多度再比较或者至少画一条稀疏曲线看看曲线是否到达平台期。这叫抽样努力校正论文里不写这个审稿人一眼就能看穿你的丰富度比较有问题。3.2 β多样性拆解物种更替和嵌套谁在主导样方间的差异β多样性说的是样方之间的物种组成差异。同样一组样方β多样性高可能有两种完全不同的生态过程第一种是环境筛选导致物种更替turnover比如阳坡和阴坡的树种完全不同第二种是排序过程物种贫乏的样方只是物种丰富样方的一个子集nestedness比如小样方的物种全都能在大样方里找到。两种过程的管理意义完全不同前者说明保护需要涵盖多种生境后者说明保护一个大斑块就够了。所以在算β多样性时我建议不要只报一个总的Sørensen相异度而是用betapart包把它分解成更替和嵌套两个部分library(betapart) # comm_pa 是0/1的物种存在-不存在矩阵 beta_multi - beta.multi(comm_pa) beta_pair - beta.pair(comm_pa, index.family sorensen) # beta.SOR 是总相异beta.SIM 是更替部分beta.NES 是嵌套部分跑完之后你会得到一个很有信息量的比例关系如果更替部分占总β多样性的80%以上说明样方梯度上的环境过滤很强接下来应该重点分析环境梯度和物种组成的关系如果嵌套部分占比高则提示取样强度不足或者存在明显的面积效应。这一步直接决定后面怎么做空间格局分析。3.3 空间格局的正式检验变差函数、Morans I和Mantel检验三件套这可能是森林生态研究里最常被用错的一组工具。先说它们各自解决什么问题Morans I检验某个变量在空间上是否呈现显著的聚集或分散。取值接近1是正自相关高值和高值相邻接近-1是负自相关高值和低值相邻接近0是随机分布。spdep包里有现成实现但关键在于怎么定义相邻——用距离阈值还是邻接矩阵。变差函数Variogram描述两点之间的半方差随距离的变化。如果半方差随着距离增加而升高然后趋于平稳说明存在空间结构如果一直是纯块金效应nugget说明样方间没有空间自相关。gstat包可以拟合理论变差函数模型。Mantel检验检验两个距离矩阵之间的相关性。比如物种组成相异度矩阵和地理距离矩阵是否相关或者物种相异度矩阵和环境差异矩阵是否相关。vegan::mantel用起来很简单但要意识到它只检验线性相关且对距离矩阵的非独立性问题比较敏感。以一个200m×200m的森林样地网格为例20个10m×10m样方。先用spdep构建距离阈值邻接权重library(spdep) # coords_mat 是样方坐标矩阵 nb - dnearneigh(coords_mat, d1 0, d2 50) # 50米内算邻居 lw - nb2listw(nb, style W) moran_test - moran.test(alpha_df$shannon, lw) print(moran_test)如果Morans I显著为正说明多样性高值样方确实在空间上扎堆。下一步你就可以理直气壮地做空间插值或空间回归而不是假装样方之间互相独立。很多经典生态数据在空间上高度自相关如果无视这一点直接做普通线性回归自由度被高估P值会系统性偏小——也就是假阳性率飙升。这个坑我在自己第一篇森林数据分析文章里踩过后来补了空间自相关检验结论方向没变但显著性和效应量都更可信了。3.4 驱动因素的空间分层变差分解Variation Partitioning的R实现知道多样性存在空间格局之后下一个问题是这个格局是环境因子驱动的还是纯粹由空间过程扩散限制、中性过程驱动的变差分解把多样性数据的总变异拆成四块环境因子单独解释的部分、空间变量单独解释的部分、环境和空间共同解释的部分、残差。这个分解逻辑特别像线性回归里的方差分析但在生态学语境下有明确的理论意义。做法分三步。第一步用vegan::rda做冗余分析得到环境模型的校正R²RsquareAdj第二步生成空间变量PCNM/dbMEM特征向量这相当于把样方坐标转化成一组正交的空间距离变量adespatial包里的dbmem函数可以直接做第三步用varpart函数做分解library(vegan) library(adespatial) # env 是标准化后的环境因子数据框 # 生成空间特征向量 spa - dbmem(coords_mat, thresh 50) # 阈值与邻接距离一致 # 变差分解 vp - varpart(comm_hellinger, env, as.data.frame(spa)) plot(vp)Hellinger转化是这一步容易被忽略的前提。原始多度数据里有大量零值直接做RDA会被少数优势种主导。decostand(comm, method hellinger)先把多度做平方根标准化让稀有物种也能在分析里发声。分解结果出来后如果环境和空间共同解释的部分占了大头通常说明环境因子本身就有空间结构这时需要进一步用偏RDA或者空间滤波方法区分机制。4. 群落稳定性不是一句空话如何用多期数据量化稳不稳4.1 稳定性的三重定义抵抗力、恢复力、时间不变异性群落稳定性在生态学里是个筐什么都能往里装。森林生态研究中最实用的定义是这三类时间不变异性Temporal invariance物种多度或群落属性总多度、多样性在多年间的波动幅度。波动越小越稳定。抵抗力Resistance受到扰动后群落属性偏离原状态的程度。偏离越小抵抗力越强。恢复力Resilience扰动后回到原状态的速度。回得越快恢复力越强。对固定样地数据来说最容易计算的是时间不变异性。一个经典的指标是群落总多度在年份间的变异系数CVCV 标准差/均值越小越稳定。另一个是从种群角度出发的群落时间稳定性定义为群落总多度的均值除以标准差μ/σ这其实就是CV的倒数。别小看这个简单的比值它有一个理论底线如果所有物种的多度完全同步波动那么群落稳定性等于所有物种种群稳定性的加权平均如果物种间异步波动群落稳定性可以超过任何一个单一物种的稳定性。这就是著名的多样性-稳定性关系的机制核心——保险效应和补偿动态。4.2 从种群到群落同步性指数怎么算要想判断一个群落的稳定性是靠少数优势种撑着还是靠物种间的异步补偿撑着需要计算物种间的同步性。ecodiversity包或自写函数都能计算这里给一个基于Loreau de Mazancourt (2008)的简化实现思路# comm_year 是多年份的样方-物种多度矩阵每行是一个样方在某年的数据 # 计算物种多度时间序列的同步性 synchrony - function(mat) { # mat 行是时间年份列是物种 var_total - var(rowSums(mat)) sum_var - sum(apply(mat, 2, var)) var_total / sum_var # 同步性1表示完全同步趋近0表示完全异步 }这个比值超过0.5基本上可以判断群落动态由同步波动主导稳定性来源单一低于0.3说明物种波动存在较强的异步补偿群落有内在稳定器。做这个分析要注意两个数据前提第一物种鉴定在多年间必须稳定同一个物种在不同年份被鉴定成不同名字同步性会离谱第二多度记录方法要一致第一年记胸径≥1cm的个体第二年却记胸径≥5cm的个体这种尺度不匹配会直接毁掉时间序列分析。4.3 用混合效应模型检验多样性是否提升了稳定性这是整条分析链里最容易被审稿人追问的部分。简单做法是算每个样方的多样性指数或者起始年份的多样性和时间稳定性μ/σ然后做线性回归。但样方之间不是独立的同一座山的不同坡向、不同海拔会引入空间结构。我的建议是用线性混合效应模型把空间区块或者坡度级别作为随机效应library(lme4) # stab_df 每行是一个样方包含 diversity、stability、block、elevation等 m1 - lmer(stability ~ diversity elevation (1|block), data stab_df) summary(m1)这里有个容易掉进去的陷阱不要在同一段数据上既算多样性又算稳定性。比如你只有2015和2020两期数据物种多样性用的是2020年的稳定性也用这两年的多度变化来算那多样性预测稳定性就成了自己预测自己。正确的做法是用较早年份比如2015的多样性预测之后年份2015-2020的稳定性或者用起始状态预测后续动态。时间错位是观察性研究中给因果推断留的一线生机错过就真的是相关分析了。5. 把格局和稳定性串成故事可视化与结果整合的实战经验5.1 地图上讲故事多样性热力插值图的R语言做法分析做完了图表不能拉胯。森林生态研究里最出效果的可视化有两类空间分布图和排序图。空间分布图方面我推荐先用ggplot2画样方散点图再用gstat的普通克里金插值生成连续表面library(gstat) # 用样方数据拟合变差函数 v - variogram(shannon ~ 1, data alpha_df, locations ~eastingnorthing) vfit - fit.variogram(v, model vgm(Sph)) # 生成预测网格grid 是覆盖整个样地范围的规则网点 krige_result - krige(shannon ~ 1, alpha_df, grid, model vfit)插值结果再用geom_raster画成热力图叠加样方点和等高线一张论文级别的空间格局图就出来了。这里的一个实操心得是插值图一定要标注样方实际位置否则读者会误以为每个像素都有实测数据支撑。透明度调低一点让底图地形信息透出来信息量会更丰富。5.2 排序图揭示环境驱动机制RDA双序图的正确打开方式RDA双序图biplot是展示样方-物种-环境三者关系的经典方式。箭头表示环境因子点代表样方或物种箭头方向与点在同侧说明正相关。用vegan::ordiplot或ggplot2手动绘制都行。手动绘制的好处是可以把样方按坡向或海拔分组着色rda_result - rda(comm_hel ~ elevation slope soil_pH canopy_open, data env) ord - ordiplot(rda_result, type none) scores_sites - as.data.frame(scores(ord, display sites)) scores_env - as.data.frame(scores(ord, display bp)) # 然后把 scores_sites 合并样方分组信息用 ggplot2 画图画这种图时切记环境箭头的长度代表该因子对排序轴的贡献大小不是生态效应的绝对大小箭头之间的夹角代表因子间的相关性锐角正相关、钝角负相关、直角无相关。这些细节写图注的时候要交代清楚。5.3 把三块结果缝合成一个逻辑闭环最后这一步属于升华环节也是最容易写出论文故事的部分。把前面算出来的多样性分解、空间格局检验和稳定性分析结果放在一起通常能形成三个层次的结论第一层森林群落的α多样性在空间上存在显著正自相关高多样性样方集中在特定生境斑块第二层β多样性以物种更替为主说明环境异质性驱动了物种组成的空间分异空间变量能解释相当一部分多样性变异第三层在控制空间结构后多样性依然与群落时间稳定性正相关且这种关系部分由物种异步性介导。这三个层次不是并列的而是一条因果链环境异质性塑造了多样性的空间格局而多样性的空间分布又影响了群落应对环境波动的能力。分析到这一步文章的叙事逻辑就完整了格局描述、机制推断、功能后果。6. 这套流程里的七个常见坑每一个我都替你踩过6.1 样方大小不一致还硬比较不同样方面积直接比较物种丰富度是经典错误。如果样方大小确实不一致必须做稀释或者用单位面积的物种数物种密度。但要注意物种密度对面积尺度的依赖很强最好在方法部分明确说明。6.2 多样性指数选了一堆但不知道回答什么问题Shannon、Simpson、Pielou均匀度、Fisher α……全放进去算一遍不等于分析深入。每个指数对应一个生态学问题选2-3个跟你的假说直接相关的就好。做多元分析时如果多个指数高度相关先做主成分提取或直接删掉冗余指标。6.3 忽视零膨胀对β多样性的影响森林数据里零太多了。一个20m×20m样方里常见树种就那么十几种稀有种只出现在少数样方。betapart处理0/1数据时对零很敏感如果数据里全是偶见种导致的你有我没有β多样性会被显著高估。解决办法是适当剔除极其稀有的物种比如只在1个样方里出现过的或者在方法里讨论稀有物种的处理策略。6.4 空间权重矩阵的选择不当dnearneigh的距离阈值设定会直接影响Morans I的结果。阈值太大会让所有样方都互为邻居检验变成全局平均比较阈值太小会出现孤立样方。实操上我会画一个样方间距离分布直方图把阈值定在距离分布的某个自然断裂点或四分之一分位数附近并做敏感性分析——换几个阈值看看结论是否稳健。6.5 用全部环境变量做变差分解不分主次环境数据经常收集了几十个变量一股脑全放进去RDA校正R²是虚高的。先做VIF检查剔除共线性强的变量VIF10的删掉或者用前向选择vegan::ordistep筛选出真正有独立解释能力的子集。6.6 多期数据的时间尺度不匹配有的样方2015年做了有的到2019年才做时间序列长度不一致直接算μ/σ做比较会偏向短序列的样方因为短序列标准差本身偏小。统一到相同观测次数或者用线性混合模型处理不平衡数据。6.7 只看P值不看效应量生态学数据样本量有限P值很容易在0.05边缘反复横跳。报告效应量比如R²、标准化回归系数、部分RDA的解释方差百分比通常比P值更有信息量。我用rstatix包或者直接手算效应量跑完模型先看量级再看显著性。7. 一个可复用的脚本流程从原始样地数据到格局-稳定性结论最后给一套可以直接改路径就跑的简化流程基于我自己的样地数据重构去掉了敏感信息只保留核心结构和注释# 森林样地多样性空间格局与稳定性分析 # 依赖包 library(tidyverse) library(vegan) library(betapart) library(spdep) library(gstat) library(adespatial) library(lme4) # 1. 读入数据长格式 raw - read.csv(forest_plot_data.csv, stringsAsFactors FALSE) # 列名约定plot_id, year, species, abundance, easting, northing, # elevation, slope, soil_pH, canopy_open # 2. 转成样方-物种矩阵取最近年份比如2018 comm - raw %% filter(year 2018) %% select(plot_id, species, abundance) %% pivot_wider(names_from species, values_from abundance, values_fill 0) # 3. α多样性 alpha - data.frame( richness specnumber(comm), shannon diversity(comm, shannon), simpson diversity(comm, simpson) ) # 4. β多样性及更替/嵌套分解 comm_pa - decostand(comm, method pa) beta_pari - beta.pair(comm_pa, index.family sorensen) # 5. 空间自相关检验 coords - raw %% filter(year 2018) %% select(plot_id, easting, northing) %% distinct() nb - dnearneigh(as.matrix(coords[, c(easting, northing)]), 0, 100) lw - nb2listw(nb, style W) moran.test(alpha$shannon, lw) # 6. 变差分解 env - raw %% filter(year 2018) %% select(plot_id, elevation, slope, soil_pH, canopy_open) %% distinct() comm_hel - decostand(comm, hellinger) spa - dbmem(as.matrix(coords[, c(easting, northing)]), thresh 100) vp - varpart(comm_hel, env[, -1], as.data.frame(spa)) plot(vp) # 7. 时间稳定性需要多年数据 stab - raw %% group_by(plot_id, year) %% summarise(total_abund sum(abundance), .groups drop) %% group_by(plot_id) %% summarise(mean_abund mean(total_abund), sd_abund sd(total_abund), stability mean_abund / sd_abund) # 8. 多样性-稳定性模型随机效应区块 mod_dat - alpha %% mutate(plot_id rownames(comm)) %% left_join(stab, by plot_id) %% left_join(raw %% select(plot_id, block) %% distinct(), by plot_id) model - lmer(stability ~ shannon elevation (1|block), data mod_dat) summary(model)脚本不是万能的但流程骨架是通用的。你在自己的数据上跑的时候最需要替换的是样方间距离阈值、环境变量集和年份定义。这些参数没有标准答案完全取决于你的样地尺度和采样设计。回头说说我自己用这套流程的经历。最早我只想算一个β多样性觉得加个空间自相关检验就够严谨了。后来发现审稿人问你的多样性格局对群落稳定性有什么影响我才被迫补上多期数据分析。补完之后才发现真正有意思的结论来自格局和稳定性的交叉部分——空间上聚集的多样性高值区往往也是时间上波动最小的区域这背后是环境缓冲和物种补偿的共同作用。如果你手里的数据恰好也有多年重复调查千万别只停留在多样性描述层面把格局和稳定性串起来分析你会看到一张完全不同的森林图景。
R语言森林生态数据分析：从空间格局到群落稳定性的完整流程
1. 从两张图看穿整片林子为什么要把空间格局和群落稳定性放在一起分析做森林生态研究的人都有一个共同的痛点样方数据好不容易收回来了除了算个Shannon-Wiener指数画张柱状图好像就不知道还能干什么了。更尴尬的是如果有人问一句你这几个样方的多样性差异是真实存在的还是因为样方离得近才像的你大概率会愣住。这个问题不是刁难它是森林生态学的核心命题之一——生物多样性的空间格局到底受什么控制以及这种格局能不能维持群落的稳定。R语言恰好是解决这一连串问题的高效工具链。它能把零散分布在样方记录表里的物种多度数据、地理坐标、环境因子数据串成一条完整的分析流水线从α多样性样方内的物种丰富度到β多样性样方间的物种组成差异再到群落时间稳定性多度在时间轴上的波动程度每一步都有成熟的R包和方法支持。这不是什么高深的新范式噱头而是把生态学里一直存在但常常被分开处理的问题真正合并起来你既要描述格局也要解释格局的形成机制还要评估格局变化对群落稳定性的后果。这篇内容适合谁看正在做森林固定样地复查数据的硕士生和博士生准备用已有数据重新挖掘故事的研究人员以及想系统梳理R语言生态分析流程的入门者。我会从数据整理讲起一路走到格局分析和稳定性评估涉及的代码逻辑、参数选择、可视化呈现方式都会结合实际样地数据展开。看完你至少能回答两件事我的样地数据里多样性在空间上是怎么分布的这种分布方式对群落的稳定性意味着什么2. 数据的底子打不好后面全是空中楼阁2.1 样地数据结构长什么样才算能用先别急着跑代码。森林样地的原始数据通常长这样每一行是一条样方-物种-多度记录配套一个样方坐标经纬度或UTM坐标、一个样方大小比如20m×20m、若干环境因子海拔、坡度、坡向、土壤pH、林冠开度等。如果你想做时间稳定性分析还需要同一批样方在不同年份的重复调查数据至少两期三期以上更理想。我见过太多人拿到数据第一件事就是算多样性指数结果算到一半发现样方ID有重复、物种名有同物异名、多度单位有的记株数有的记盖度。R语言里做数据处理的第一步永远不是分析而是把这张表清洗成每一行唯一、每一列明确的标准长数据格式。我的建议是先用tidyr和dplyr做一次全量体检library(tidyverse) # 假设你的原始数据叫 raw_data.csv df - read.csv(raw_data.csv, stringsAsFactors FALSE) # 查看缺失值、重复样方、异常多度 summary(df) df %% count(plot_id) %% filter(n 1) # 检查样方是否重复 df %% filter(abundance 0) # 多度是否出现0或负值这里有个很容易忽略的细节多度为0的记录在长数据里应不应该保留我的习惯是分析之前剔除但保留一份原始备份。因为有些多样性指数比如Chao1需要知道样方里未被观测到的物种的估计这个信息是从仅有1株的物种数和仅有2株的物种数推算的0记录不会参与计算但会影响你对数据覆盖率的判断。先把原始数据存档再在分析副本里操作这是R项目的第一条铁律。2.2 坐标和其他环境因子的预处理不标准化就等着结果骗你样方坐标到手之后先确认坐标系是否统一。森林样地数据经常混着GPS的经纬度WGS84和RTK测的UTM坐标混在一起算距离会得到荒谬的结果。R语言里处理坐标系常用sf包但如果你只是算样方间的地理距离用geosphere包的distm函数就够了library(geosphere) # coords 是 data.frame包含经纬度两列 d_geo - distm(coords, fun distGeo) # 返回距离矩阵单位米环境因子方面别把所有变量一股脑全塞进模型。先做相关性筛查海拔和温度、坡度和土壤含水量经常高度相关保留生态学意义最强的一个即可。连续变量在进入多元模型之前要标准化z-score否则量纲大的变量会主导距离计算和回归系数。这里用vegan::decostand(x, method standardize)是最省事的做法。2.3 空间自相关的第一印象先画一张样方分布图在动手算任何指数之前先把样方位置画出来用颜色或气泡大小表示物种丰富度。这一步治好了我无数次拿着数据硬跑模型的毛病。样方如果集中在山谷、沿沟系分布那你的空间格局很可能是取样设计造成的伪格局如果样方覆盖了整个坡面那分析结果才有生态解释价值。用ggplot2加ggspatial做底图是最快的library(ggplot2) ggplot(coords, aes(x easting, y northing, size richness)) geom_point(alpha 0.7) coord_equal() theme_minimal()如果样方数量比较多超过30个还可以叠加一个简易的插值表面用gstat包的idw做反距离加权插值快速看看多样性高值区在空间上是连片还是离散。这一步的输出不一定会进论文但它决定了你后面选哪种空间分析方法。3. 多样性空间格局的分析主线从α到β再到γ的拆解逻辑3.1 α多样性的多面性丰富度、Shannon、Simpson到底该怎么选α多样性是每个样方内部的物种多样性最常用的三个指数各有脾气。物种丰富度Species richness就是数数直观但容易被稀有种带偏Shannon指数对常见种敏感数值越大说明不确定性越高也就是多样性越高Simpson指数对优势种敏感反映的是随机抽两个个体它们是不同物种的概率。没有哪个指数绝对正确关键看你问什么问题。如果你关心的是群落对干扰的响应Simpson指数往往更稳因为它主要受优势种影响而优势种通常在干扰下最先变化如果你关心的是生境异质性对物种共存的影响物种丰富度更直接。我的建议是三者都算放进一个数据框里后续分析各取所需。vegan的diversity函数一条命令搞定library(vegan) # comm 是样方-物种多度矩阵行是样方列是物种 rich - specnumber(comm) shannon - diversity(comm, index shannon) simpson - diversity(comm, index simpson) alpha_df - data.frame(plot_id rownames(comm), rich, shannon, simpson)还有个容易被忽略的步骤稀释曲线rarefaction。如果不同样方的总多度差异很大直接比较物种丰富度会不公平——你采到100个个体的样方当然比只采到20个个体的样方物种更多。用vegan::rrarefy把所有样方稀释到相同多度再比较或者至少画一条稀疏曲线看看曲线是否到达平台期。这叫抽样努力校正论文里不写这个审稿人一眼就能看穿你的丰富度比较有问题。3.2 β多样性拆解物种更替和嵌套谁在主导样方间的差异β多样性说的是样方之间的物种组成差异。同样一组样方β多样性高可能有两种完全不同的生态过程第一种是环境筛选导致物种更替turnover比如阳坡和阴坡的树种完全不同第二种是排序过程物种贫乏的样方只是物种丰富样方的一个子集nestedness比如小样方的物种全都能在大样方里找到。两种过程的管理意义完全不同前者说明保护需要涵盖多种生境后者说明保护一个大斑块就够了。所以在算β多样性时我建议不要只报一个总的Sørensen相异度而是用betapart包把它分解成更替和嵌套两个部分library(betapart) # comm_pa 是0/1的物种存在-不存在矩阵 beta_multi - beta.multi(comm_pa) beta_pair - beta.pair(comm_pa, index.family sorensen) # beta.SOR 是总相异beta.SIM 是更替部分beta.NES 是嵌套部分跑完之后你会得到一个很有信息量的比例关系如果更替部分占总β多样性的80%以上说明样方梯度上的环境过滤很强接下来应该重点分析环境梯度和物种组成的关系如果嵌套部分占比高则提示取样强度不足或者存在明显的面积效应。这一步直接决定后面怎么做空间格局分析。3.3 空间格局的正式检验变差函数、Morans I和Mantel检验三件套这可能是森林生态研究里最常被用错的一组工具。先说它们各自解决什么问题Morans I检验某个变量在空间上是否呈现显著的聚集或分散。取值接近1是正自相关高值和高值相邻接近-1是负自相关高值和低值相邻接近0是随机分布。spdep包里有现成实现但关键在于怎么定义相邻——用距离阈值还是邻接矩阵。变差函数Variogram描述两点之间的半方差随距离的变化。如果半方差随着距离增加而升高然后趋于平稳说明存在空间结构如果一直是纯块金效应nugget说明样方间没有空间自相关。gstat包可以拟合理论变差函数模型。Mantel检验检验两个距离矩阵之间的相关性。比如物种组成相异度矩阵和地理距离矩阵是否相关或者物种相异度矩阵和环境差异矩阵是否相关。vegan::mantel用起来很简单但要意识到它只检验线性相关且对距离矩阵的非独立性问题比较敏感。以一个200m×200m的森林样地网格为例20个10m×10m样方。先用spdep构建距离阈值邻接权重library(spdep) # coords_mat 是样方坐标矩阵 nb - dnearneigh(coords_mat, d1 0, d2 50) # 50米内算邻居 lw - nb2listw(nb, style W) moran_test - moran.test(alpha_df$shannon, lw) print(moran_test)如果Morans I显著为正说明多样性高值样方确实在空间上扎堆。下一步你就可以理直气壮地做空间插值或空间回归而不是假装样方之间互相独立。很多经典生态数据在空间上高度自相关如果无视这一点直接做普通线性回归自由度被高估P值会系统性偏小——也就是假阳性率飙升。这个坑我在自己第一篇森林数据分析文章里踩过后来补了空间自相关检验结论方向没变但显著性和效应量都更可信了。3.4 驱动因素的空间分层变差分解Variation Partitioning的R实现知道多样性存在空间格局之后下一个问题是这个格局是环境因子驱动的还是纯粹由空间过程扩散限制、中性过程驱动的变差分解把多样性数据的总变异拆成四块环境因子单独解释的部分、空间变量单独解释的部分、环境和空间共同解释的部分、残差。这个分解逻辑特别像线性回归里的方差分析但在生态学语境下有明确的理论意义。做法分三步。第一步用vegan::rda做冗余分析得到环境模型的校正R²RsquareAdj第二步生成空间变量PCNM/dbMEM特征向量这相当于把样方坐标转化成一组正交的空间距离变量adespatial包里的dbmem函数可以直接做第三步用varpart函数做分解library(vegan) library(adespatial) # env 是标准化后的环境因子数据框 # 生成空间特征向量 spa - dbmem(coords_mat, thresh 50) # 阈值与邻接距离一致 # 变差分解 vp - varpart(comm_hellinger, env, as.data.frame(spa)) plot(vp)Hellinger转化是这一步容易被忽略的前提。原始多度数据里有大量零值直接做RDA会被少数优势种主导。decostand(comm, method hellinger)先把多度做平方根标准化让稀有物种也能在分析里发声。分解结果出来后如果环境和空间共同解释的部分占了大头通常说明环境因子本身就有空间结构这时需要进一步用偏RDA或者空间滤波方法区分机制。4. 群落稳定性不是一句空话如何用多期数据量化稳不稳4.1 稳定性的三重定义抵抗力、恢复力、时间不变异性群落稳定性在生态学里是个筐什么都能往里装。森林生态研究中最实用的定义是这三类时间不变异性Temporal invariance物种多度或群落属性总多度、多样性在多年间的波动幅度。波动越小越稳定。抵抗力Resistance受到扰动后群落属性偏离原状态的程度。偏离越小抵抗力越强。恢复力Resilience扰动后回到原状态的速度。回得越快恢复力越强。对固定样地数据来说最容易计算的是时间不变异性。一个经典的指标是群落总多度在年份间的变异系数CVCV 标准差/均值越小越稳定。另一个是从种群角度出发的群落时间稳定性定义为群落总多度的均值除以标准差μ/σ这其实就是CV的倒数。别小看这个简单的比值它有一个理论底线如果所有物种的多度完全同步波动那么群落稳定性等于所有物种种群稳定性的加权平均如果物种间异步波动群落稳定性可以超过任何一个单一物种的稳定性。这就是著名的多样性-稳定性关系的机制核心——保险效应和补偿动态。4.2 从种群到群落同步性指数怎么算要想判断一个群落的稳定性是靠少数优势种撑着还是靠物种间的异步补偿撑着需要计算物种间的同步性。ecodiversity包或自写函数都能计算这里给一个基于Loreau de Mazancourt (2008)的简化实现思路# comm_year 是多年份的样方-物种多度矩阵每行是一个样方在某年的数据 # 计算物种多度时间序列的同步性 synchrony - function(mat) { # mat 行是时间年份列是物种 var_total - var(rowSums(mat)) sum_var - sum(apply(mat, 2, var)) var_total / sum_var # 同步性1表示完全同步趋近0表示完全异步 }这个比值超过0.5基本上可以判断群落动态由同步波动主导稳定性来源单一低于0.3说明物种波动存在较强的异步补偿群落有内在稳定器。做这个分析要注意两个数据前提第一物种鉴定在多年间必须稳定同一个物种在不同年份被鉴定成不同名字同步性会离谱第二多度记录方法要一致第一年记胸径≥1cm的个体第二年却记胸径≥5cm的个体这种尺度不匹配会直接毁掉时间序列分析。4.3 用混合效应模型检验多样性是否提升了稳定性这是整条分析链里最容易被审稿人追问的部分。简单做法是算每个样方的多样性指数或者起始年份的多样性和时间稳定性μ/σ然后做线性回归。但样方之间不是独立的同一座山的不同坡向、不同海拔会引入空间结构。我的建议是用线性混合效应模型把空间区块或者坡度级别作为随机效应library(lme4) # stab_df 每行是一个样方包含 diversity、stability、block、elevation等 m1 - lmer(stability ~ diversity elevation (1|block), data stab_df) summary(m1)这里有个容易掉进去的陷阱不要在同一段数据上既算多样性又算稳定性。比如你只有2015和2020两期数据物种多样性用的是2020年的稳定性也用这两年的多度变化来算那多样性预测稳定性就成了自己预测自己。正确的做法是用较早年份比如2015的多样性预测之后年份2015-2020的稳定性或者用起始状态预测后续动态。时间错位是观察性研究中给因果推断留的一线生机错过就真的是相关分析了。5. 把格局和稳定性串成故事可视化与结果整合的实战经验5.1 地图上讲故事多样性热力插值图的R语言做法分析做完了图表不能拉胯。森林生态研究里最出效果的可视化有两类空间分布图和排序图。空间分布图方面我推荐先用ggplot2画样方散点图再用gstat的普通克里金插值生成连续表面library(gstat) # 用样方数据拟合变差函数 v - variogram(shannon ~ 1, data alpha_df, locations ~eastingnorthing) vfit - fit.variogram(v, model vgm(Sph)) # 生成预测网格grid 是覆盖整个样地范围的规则网点 krige_result - krige(shannon ~ 1, alpha_df, grid, model vfit)插值结果再用geom_raster画成热力图叠加样方点和等高线一张论文级别的空间格局图就出来了。这里的一个实操心得是插值图一定要标注样方实际位置否则读者会误以为每个像素都有实测数据支撑。透明度调低一点让底图地形信息透出来信息量会更丰富。5.2 排序图揭示环境驱动机制RDA双序图的正确打开方式RDA双序图biplot是展示样方-物种-环境三者关系的经典方式。箭头表示环境因子点代表样方或物种箭头方向与点在同侧说明正相关。用vegan::ordiplot或ggplot2手动绘制都行。手动绘制的好处是可以把样方按坡向或海拔分组着色rda_result - rda(comm_hel ~ elevation slope soil_pH canopy_open, data env) ord - ordiplot(rda_result, type none) scores_sites - as.data.frame(scores(ord, display sites)) scores_env - as.data.frame(scores(ord, display bp)) # 然后把 scores_sites 合并样方分组信息用 ggplot2 画图画这种图时切记环境箭头的长度代表该因子对排序轴的贡献大小不是生态效应的绝对大小箭头之间的夹角代表因子间的相关性锐角正相关、钝角负相关、直角无相关。这些细节写图注的时候要交代清楚。5.3 把三块结果缝合成一个逻辑闭环最后这一步属于升华环节也是最容易写出论文故事的部分。把前面算出来的多样性分解、空间格局检验和稳定性分析结果放在一起通常能形成三个层次的结论第一层森林群落的α多样性在空间上存在显著正自相关高多样性样方集中在特定生境斑块第二层β多样性以物种更替为主说明环境异质性驱动了物种组成的空间分异空间变量能解释相当一部分多样性变异第三层在控制空间结构后多样性依然与群落时间稳定性正相关且这种关系部分由物种异步性介导。这三个层次不是并列的而是一条因果链环境异质性塑造了多样性的空间格局而多样性的空间分布又影响了群落应对环境波动的能力。分析到这一步文章的叙事逻辑就完整了格局描述、机制推断、功能后果。6. 这套流程里的七个常见坑每一个我都替你踩过6.1 样方大小不一致还硬比较不同样方面积直接比较物种丰富度是经典错误。如果样方大小确实不一致必须做稀释或者用单位面积的物种数物种密度。但要注意物种密度对面积尺度的依赖很强最好在方法部分明确说明。6.2 多样性指数选了一堆但不知道回答什么问题Shannon、Simpson、Pielou均匀度、Fisher α……全放进去算一遍不等于分析深入。每个指数对应一个生态学问题选2-3个跟你的假说直接相关的就好。做多元分析时如果多个指数高度相关先做主成分提取或直接删掉冗余指标。6.3 忽视零膨胀对β多样性的影响森林数据里零太多了。一个20m×20m样方里常见树种就那么十几种稀有种只出现在少数样方。betapart处理0/1数据时对零很敏感如果数据里全是偶见种导致的你有我没有β多样性会被显著高估。解决办法是适当剔除极其稀有的物种比如只在1个样方里出现过的或者在方法里讨论稀有物种的处理策略。6.4 空间权重矩阵的选择不当dnearneigh的距离阈值设定会直接影响Morans I的结果。阈值太大会让所有样方都互为邻居检验变成全局平均比较阈值太小会出现孤立样方。实操上我会画一个样方间距离分布直方图把阈值定在距离分布的某个自然断裂点或四分之一分位数附近并做敏感性分析——换几个阈值看看结论是否稳健。6.5 用全部环境变量做变差分解不分主次环境数据经常收集了几十个变量一股脑全放进去RDA校正R²是虚高的。先做VIF检查剔除共线性强的变量VIF10的删掉或者用前向选择vegan::ordistep筛选出真正有独立解释能力的子集。6.6 多期数据的时间尺度不匹配有的样方2015年做了有的到2019年才做时间序列长度不一致直接算μ/σ做比较会偏向短序列的样方因为短序列标准差本身偏小。统一到相同观测次数或者用线性混合模型处理不平衡数据。6.7 只看P值不看效应量生态学数据样本量有限P值很容易在0.05边缘反复横跳。报告效应量比如R²、标准化回归系数、部分RDA的解释方差百分比通常比P值更有信息量。我用rstatix包或者直接手算效应量跑完模型先看量级再看显著性。7. 一个可复用的脚本流程从原始样地数据到格局-稳定性结论最后给一套可以直接改路径就跑的简化流程基于我自己的样地数据重构去掉了敏感信息只保留核心结构和注释# 森林样地多样性空间格局与稳定性分析 # 依赖包 library(tidyverse) library(vegan) library(betapart) library(spdep) library(gstat) library(adespatial) library(lme4) # 1. 读入数据长格式 raw - read.csv(forest_plot_data.csv, stringsAsFactors FALSE) # 列名约定plot_id, year, species, abundance, easting, northing, # elevation, slope, soil_pH, canopy_open # 2. 转成样方-物种矩阵取最近年份比如2018 comm - raw %% filter(year 2018) %% select(plot_id, species, abundance) %% pivot_wider(names_from species, values_from abundance, values_fill 0) # 3. α多样性 alpha - data.frame( richness specnumber(comm), shannon diversity(comm, shannon), simpson diversity(comm, simpson) ) # 4. β多样性及更替/嵌套分解 comm_pa - decostand(comm, method pa) beta_pari - beta.pair(comm_pa, index.family sorensen) # 5. 空间自相关检验 coords - raw %% filter(year 2018) %% select(plot_id, easting, northing) %% distinct() nb - dnearneigh(as.matrix(coords[, c(easting, northing)]), 0, 100) lw - nb2listw(nb, style W) moran.test(alpha$shannon, lw) # 6. 变差分解 env - raw %% filter(year 2018) %% select(plot_id, elevation, slope, soil_pH, canopy_open) %% distinct() comm_hel - decostand(comm, hellinger) spa - dbmem(as.matrix(coords[, c(easting, northing)]), thresh 100) vp - varpart(comm_hel, env[, -1], as.data.frame(spa)) plot(vp) # 7. 时间稳定性需要多年数据 stab - raw %% group_by(plot_id, year) %% summarise(total_abund sum(abundance), .groups drop) %% group_by(plot_id) %% summarise(mean_abund mean(total_abund), sd_abund sd(total_abund), stability mean_abund / sd_abund) # 8. 多样性-稳定性模型随机效应区块 mod_dat - alpha %% mutate(plot_id rownames(comm)) %% left_join(stab, by plot_id) %% left_join(raw %% select(plot_id, block) %% distinct(), by plot_id) model - lmer(stability ~ shannon elevation (1|block), data mod_dat) summary(model)脚本不是万能的但流程骨架是通用的。你在自己的数据上跑的时候最需要替换的是样方间距离阈值、环境变量集和年份定义。这些参数没有标准答案完全取决于你的样地尺度和采样设计。回头说说我自己用这套流程的经历。最早我只想算一个β多样性觉得加个空间自相关检验就够严谨了。后来发现审稿人问你的多样性格局对群落稳定性有什么影响我才被迫补上多期数据分析。补完之后才发现真正有意思的结论来自格局和稳定性的交叉部分——空间上聚集的多样性高值区往往也是时间上波动最小的区域这背后是环境缓冲和物种补偿的共同作用。如果你手里的数据恰好也有多年重复调查千万别只停留在多样性描述层面把格局和稳定性串起来分析你会看到一张完全不同的森林图景。
