基于MaxEnt模型的怒江金丝猴潜在适生区预测

肖 祺，刘 成，2，徐梦珍，周雄冬

（1.中国水利水电科学研究院，北京 100038；2.国际泥沙研究培训中心，北京 100048；3.清华大学 土木水利学院，北京 100084）

1 研究背景

怒江金丝猴（Rhinopithecus strykeri）是我国一级重点保护野生动物，2010年被首次发现，是目前世界上记录的第五种金丝猴，被世界自然保护联盟（IUCN）列为极危物种。怒江金丝猴首先被发现于缅甸，此后由中国科学家在高黎贡山地区发现其踪迹。种群数量仅几百只，生活在1700 m至3100 m的中山湿性常绿阔叶林和温凉性针叶林中。国内对于仰鼻猴属的研究主要集中于早年发现的川金丝猴和滇金丝猴，怒江金丝猴作为新发现物种，其相关成果多局限于生物学描述及种群动态研究[1-3]，而对其适生栖息地分布格局及变化规律缺乏科学研究。金丝猴生活在横断山区西缘，该地区地形气候条件复杂，是全球生物多样性热点区域之一[4]，同时此地还存在的水能开发对生态的影响问题[5]，对其进行适生栖息地分布研究对于保护我国生物多样性具有重要意义。近年来使用生态位模型进行生物多样性保护、入侵生物潜在分布区的研究成为新的热点[6-8]。常见的生态位模型有MaxEnt、Domain、Bioclim、ENFA和GARP等。王国峥等[9]使用4种生态位模型对我国特有孑遗植物金钱松的潜在适生区进行预测，发现MaxEnt模型结果更精确；曹向锋等[10]使用5种生态位模型对入侵生物黄顶菊的潜在适生区进行预测，筛选后采纳了MaxEnt的预测结果；段义忠等[11]在不同时段使用3种生态位模型预测，结果也表明MaxEnt模型的最适生预测与现实分布保持一致。本研究使用MaxEnt模型对怒江金丝猴现代以及未来气候情景下的适生区进行研究，为深入了解这一物种的时空分布格局、进行生态适生区保护提供理论依据。

2 研究方法

2.1 模型分析

2.1.1 模型含义 MaxEnt最大熵模型（maximum entropy model）是以最大熵理论（the theory of maximum entropy）为基础的物种地理尺度空间分布模型，是目前应用最广、表现最好的生态位模型之一[12]，主要用于预测物种分布区。本研究中使用开发者Steven Phillips提供的3.4.1版本。

2.1.2 模型参数 MaxEnt模型是一种基于机器学习的复杂模型，准确率高但可转移性差，大多数研究采用默认参数来构建模型，默认参数的设置来源于模型开发者对部分地区的266个物种数据的测试，且以模拟其现实分布为目的[13]，所以在模拟物种未来潜在适生区时，继续使用默认参数会使MaxEnt模型产生过拟合的问题，使预测结果不可靠[14]，需调用ENMeval数据包[15]对MaxEnt模型参数中的调控倍频（regularization multiplier，RM）和特征组合（feature combination，FC）进行调整，分析各种参数条件下模型的复杂度，选取最低复杂度的模型参数。特征类型分为：线型（linear）-L，二次型（quadratic）-Q，铰链型（hinge）-H，乘积型（product）-P和阈值型（threshold）-T。RM的调整区间为0.5～4.0，每次增加0.5；由于样本量较少，FC选取2种组合方式，分别为L、LQ（同时使用线型与二次型）。利用ENMeval所输出的各类数据选择最适生的RM与FC。

2.1.3 模型评价 运用赤池信息量准则（AIC）评价模型的复杂度和拟合优度，AIC越小，模型越好，通常选择AIC最小的模型；运用AUC（ROC曲线下方的面积大小）、AUCtest（训练AUC）、AUCdiff（测试AUC值之差）、orMTP（最小训练集遗漏率）和or10（10%训练遗漏率）[16-18]评估不同参数组合的拟合度，综合以上数据选择出最适合怒江金丝猴的参数组合。利用最优参数组合在MaxEnt模型中进行建模和预测适生区分布，AUC评价模型的性能好坏时取值范围为0.5～1，AUC值越大表示模型性能越好，当AUC值≥0.9时预测结果精度高且较为可靠。

2.1.4 适生区划分 根据MaxEnt模型的结果，按照0～1的取值范围输出适生程度，0为最不适生，1为最适生。参考联合国政府间气候变化专门委员会第5次评估报告中的描述[19]，将适生区划分为4个等级，分别为不适生区（0～0.3）、低适生区（0.3～0.5）、中适生区（0.5～0.7）和高适生区（0.7～1.0）。

2.2 物种分布数据物种分布数据来源于两方面，一是全球生物多样性信息网络（https：//www.gbif.org/）、中国数字植物标本馆平台（http：//www.cvh.ac.cn/）网站；二是相关文献，根据怒江金丝猴相应的出没地点信息，借助 GeoNames（http：//www.geonames.org/）查询其经纬度数据。剔除冗余数据点、坐标缺失的点及坐标无意义的点，获得44个怒江金丝猴分布点。

2.3 环境变量获取及分析现代及未来气候变量数据来源于WorldClim（https：//www.worldclim.org/），包含19个气候变量，最冷季度和最干季度等气候变量来源于每月的温度和降雨量值，目的是生成更具生物意义的变量，由于不同地区出现极端环境的月份不同，所以在生成变量时不固定月份而是以最冷季度和最干季度等来表示。空间分辨率为1 km×1 km，时间分辨率分别为现代（1970—2000年），未来（a：2021—2040年；b：2041—2060年；c：2061—2080年；d：2081—2100年）。未来气候模式基于第六次耦合模式比较计划（CMIP6），使用新全球气候模式BCC-CSM 2-MR和中等未来共享社会经济路径（SSP245）数据。将19个气候变量进行聚类和皮尔逊相关性分析，在此基础上筛选生物气候变量。聚类分析将生物气候变量分为七组（图1）：（1）最冷月份最低温度、最干季度平均温度、最冷季度平均温度；（2）年降水量、最暖季度降水量、最湿月份降水量、最湿季度降水量；（3）最热月份最高温度、最湿季度平均温度、最暖季度平均温度；（4）温度变化系数、年温差；（5）等温性、昼夜温差月均值、降水量变化方差；（6）年均温；（7）最冷季度降水量、最干月份降水量、最干季度降水量。在相关性分析中，各组内变量间的相关性均大于0.9，与聚类分析的结果相符。组内挑选与其他变量最不相关的一个变量（表1），保留最重要的生态变量，单独计算各变量与其他变量的相关性，取其绝对值的和作为单变量总相关性，组内单变量总相关性最低的变量即为与其他变量最不相关的变量。筛选出7个变量（年均温、等温性、温度变化系数、最湿季度平均温度、最冷季度平均温度、最暖季度降水量、最冷季度降水量），作为预测怒江金丝猴适生区的气候变量。

图1 不同生物气候变量的聚类分析结果

表1 相关性分析结果

在气候变量的基础上加入影响怒江金丝猴生长发育的其他环境变量，包括：高程（DEM）、土地利用（LUCC）、归一化植被指数（NDVI）。高程数据来源于SRTM，土地利用数据与归一化植被指数数据来源于中科院地理科学与资源研究所（https：//www.resdc.cn/Default.aspx）。通过ArcGIS剪裁与转换得到与环境数据相同的空间分辨率。最终选择10个变量（年均温、等温性、温度变化系数、最湿季度平均温度、最冷季度平均温度、最暖季度降水量、最冷季度降水量、高程、土地利用、归一化植被指数），作为预测怒江金丝猴适生区的环境变量（表2）。

表2 环境变量及其参数

2.4 模型参数选择在小样本情况下AIC转变成AICc，随着样本量增加，AICc收敛成AIC。所以AICc可以应用于任何样本大小的情况下。得到各个模型的AICc值后，每个AICc值与最小的AICc值相减，得到ΔAICc。通过ENMeval计算 16种不同参数设置组合的 ΔAICc值（表3）。在 16种参数组合中，ΔAICc＜2的参数组合共有2组，分别为：RM=3，RM=3.5，FC均选取 LQ。它们的 ΔAICc分别为1.8837和0.0520，2组参数组合的ΔAICc都在接受范围内，其复杂度和拟合度均较低，具有较高的可靠性。进一步对 ENMeval生成的其他评估指标进行比较，当 RM=3，FC选取 LQ时，其 AUCdiff、orMTP和or10值均为最低，表明RM=3，FC选取LQ时的拟合度更好，并且其AUCtest的值更高，表明其能更好地区分测试分布点和背景点，因此，以RM=3，FC选取LQ作为最终的设置参数。在该设置参数条件下，MaxEnt模型模拟怒江金丝猴适生区的平均AUC值为0.9995，表明模型的预测效果良好。

表3 不同参数设置下模型评估结果

2.5 主导环境变量分析Maxent模型在训练时会记录变量为适应模型做出的贡献，最终将变量贡献量转换为百分数输出为贡献率，置换重要性是由该变量在训练过程中，训练点的随机置换和由此导致的AUC值的下降幅度来决定，AUC下降幅度越大，模型对该环境变量依赖越大。不同气候变量对MaxEnt模型的贡献率和置换重要性（表2）结果显示，共有4个气候变量对模型的贡献率大于10%，分别为等温性（11.95%）、温度变化系数（28.60%）、最冷季度平均温度（16.44%）和 NDVI（17.39%），其累计贡献率达74.38%，其中最冷季度平均温度和NDVI的重要性较高，最冷季度平均温度重要性达到50.17%。为了减少气候变量之间的相互影响，以进一步解释哪些变量在模型中是最重要的，对10个变量的重要性进行刀切（Jackknife）法检测（图2）。当使用刀切法时，会创建多种模型，包括单独使用每个变量创建一个模型、依次排除每个变量并使用剩余变量创建模型和使用所有变量创建模型。结果显示，单独使用年均温、最湿季度平均温度和土地利用时，几乎没有任何增益，表明这些参数对预测怒江金丝猴适生区的帮助很小；当单独使用时，增益最高的环境变量是温度变化系数，因此它拥有最重要的信息，这与贡献率结果相符合。当忽略单个变量时，降低增益最多的环境变量是最冷季度降水量，因此它拥有其他变量中不存在的最多信息。综上所述，温度变化系数、最冷季度平均温度、最冷季度降水量和归一化植被指数是影响怒江金丝猴分布的主导变量。

图2 变量重要性的刀切法检测结果

3 结果与讨论

3.1 怒江金丝猴现代适生区模拟现代背景下，44个怒江金丝猴分布点中属于高适生区、中适生区、低适生区和不适生区的点数分别为 17个，19个，7个，1个，相应比例分别为38.64%、43.18%、15.91%和2.27%，模拟结果（图3）中我国境内现存点均在高适生范围内，表明了预测结果的准确性和MaxEnt模型用于评价怒江金丝猴适生区的可行性。怒江金丝猴在中国的现代适生区总面积约为30 790 km2，其中高适生区面积为6 207 km2，占适生区总面积的20.16%；主要适生区分布于云南西部的三江并流区域。高适生区集中于怒江两岸的高黎贡山与怒山，垂直高差较大，有着典型的山地垂直植被带和多样的动植物资源，为动植物的生存提供了自低向高的海拔条件，具有典型的亚热带气候特征。

图3 怒江金丝猴的分布点与模拟现代适生区分布

3.2 怒江金丝猴适生区预测现代至未来4个阶段，怒江金丝猴总适生区均有增长（图4）。根据模型对未来适生区的预测（图5），云南西部的适生区在2021—2040年面积最大并且是高适生区最多的时间段，与现代相比总适生区面积增加98.59%，高适生区面积增加64.91%，主要表现为澜沧江东部适生区的扩大，并且延伸至金沙江以东；2021—2080年，各类型适生区面积都呈下降趋势，向高适生区收缩，更集中于怒江与澜沧江之间的山区；2041—2060年的适生区面积与轮廓和2081—2100年相似，分别比现代增加了79.36%和79.01%。

图4 怒江金丝猴现代以及未来各适生区面积

图5 未来气候条件下怒江金丝猴的适生区分布预测

3.3 讨论本研究中，综合变量贡献率、重要性以及刀切法检测，发现温度变化系数、最冷季度平均温度、最冷季度降水量和归一化植被指数是影响怒江金丝猴分布的主导变量。怒江金丝猴为典型的森林树栖动物，生活在1700 m至3100 m的原始森林中，其分布受到温度、降水及栖息地植被情况的多重影响[1]，目前已知存在点高黎贡山区域海拔高差较大，且存在明显的山地植被垂直分布带[20]，而生境中的植物对怒江金丝猴的觅食行为产生影响，冬季的低温及降雪令食物变得匮乏，使其向低海拔地区移动生活，上述分析表明与低温期相关的变量是决定怒江金丝猴地理空间分布的重要环境因素。

相较于已知存在区怒江以西的高黎贡山，模拟现代高适生区多出了怒江和金沙江之间的怒山山脉。实际情况下在这两个区域并未观测到怒江金丝猴的存在，在生境相似的情况下模拟结果会偏向于存在可能性高[21]，怒江以东的区域由于怒江、澜沧江的隔离作用，即使西南山区生境相似，也并没有怒江金丝猴生活在此区域，同为仰鼻猴属的滇金丝猴分布于澜沧江与金沙江之间云岭山脉两侧，主要活动在2500 m至4500 m高度的云杉、冷杉林中，两个物种之间没有发现交流迹象。在仅考虑下垫面情况下，模拟结果显示在台湾地区也存在怒江金丝猴的高适生区，但由于地理上的隔离怒江金丝猴不可能出现在台湾，陈之瑞等[22]总结了中国西南地区与台湾植物的间断分布现象，表示两者植物区系上有极大的相似性，台湾的中央山脉最高主峰近4000 m，从植被和海拔气候都与高黎贡山有十分相近的条件，这是模拟在台湾也有高适生区分布的客观原因，在进行珍稀物种保护时可考虑让台湾地区成为异地迁移保护的备选区。

SSP245为综合考虑了SSP-RCP的未来气候情景，其中共享社会经济路径为 “中间路线”的情景（SSP2），即在整个21世纪都将延续历史发展的情景，典型浓度排放路径为温室气体中等排放（RCP4.5）。在此情景下，未来中国西南部气温有显著的增长趋势，降水在不同时段内保持平稳变化[23]。受到全球变暖的影响，对低温敏感的怒江金丝猴适生区可能会向外扩张，由澜沧江向东扩大，在近现代经历一个高速扩张的时期，适生区达到最大，接着下降到较为稳定的2041—2060年和2081—2100年时段，该时期相较于现代南北方向上缩减，东西方向上扩大。综合未来气候情景下的适生区预测，怒江金丝猴可以长久的生活在高黎贡山区域，并且适生区有所扩大，可以在怒江以东以寻找其踪迹，为该种金丝猴的种群保育工作和保护区划分提供了依据。

4 结论

本研究利用MaxEnt优化模型模拟了怒江金丝猴现代以及未来气候情景下4个时段的适生区。模型的AUC值为0.9995，模型表现良好。结果表明，影响怒江金丝猴生存的主导变量为温度变化系数、最冷季度平均温度、最冷季度降水量和归一化植被指数。怒江金丝猴在当前气候条件下适生区总面积为30 790 km2，高适生区面积为6 207 km2，主要集中在云南的高黎贡山区域，怒江与澜沧江之间的怒山也是最适生怒江金丝猴生存的地方。未来条件下其适生区增加，2021—2040年为增加最多的时间段，相较于现代增加了总适生区98.59%，西南山区部分向东扩张延伸至金沙江。本研究为怒江金丝猴这一极危物种的种群研究和保护区建立提供了科学依据。