一个世纪以来, 气候变化已对生物多样性及生态系统稳定造成了巨大影响^[1—2]。生境是生物赖以生存的基础条件和生物多样性存在的场所^[3], 而气候变化通过影响生境的数量、质量和空间配置以制约物种的分布、生存和繁衍^[4]。因此, 识别气候变化下物种生境的变化特征对生物多样性的保护有重要意义^[5]。

目前已广泛应用多种生态位模型结合气候情景数据进行物种生境的模拟和预测。最大熵模型(maximum entropy model, MaxEnt)能将物种分布和环境变量相关联, 并得出两者之间的关系, 是目前预测物种分布的最佳算法之一^[6]。气候情景是对未来气候发展变化的可能性描述, 国际耦合模式比较计划(CMIP)提供的气候模式数据可支持未来气候评估^[7], 目前已发布到第六版(CMIP6)。共享社会经济路径(SSPs)是CMIP6建立的新发展情景框架, 其设定基于现在国家、区域的实际情况和未来发展规划^[8], 相较CMIP5, CMIP6模式对中国区域的气温和降水模拟能力更佳^[9]。

大熊猫是全球珍稀野生动物和生物多样性保护的旗舰物种^[10], 也是物种生境模拟预测的热门对象^[11], 具有扩散和繁殖能力低、食物单一、分布范围窄、对生境变化敏感等特点, 极易受气候变化带来的不利影响^[12]。已有大熊猫生境研究的内容多为过往年份的生境时空变化及质量评价^{[10, 13—14]}, 影响因子探究^[15—17], 关于未来气候情景下大熊猫生境变化的研究较少, 且多基于CMIP5的数据^[18—21], 亟需基于新的气候模式评估未来大熊猫生境质量演变特征。

大熊猫国家公园以大熊猫种群及其生境保护为主要目标, 是我国生态安全屏障的关键区域, 在应对气候变化危机中起着先行者作用^[22]。大熊猫国家公园区域内地形地貌复杂、气候垂直分异大, 生境破碎化严重^[23], 而已有的大熊猫国家公园研究集中于生态旅游^[24]、自然教育^[25]、社区治理^[26]、监测体系构建^[27]等方面, 针对国家公园气候适应性管理的研究还较少^[28]。因此, 识别气候变化下的大熊猫生境变化以进行适应性保护管理, 已成为大熊猫国家公园面临的重要任务之一。

本研究以大熊猫国家公园岷山片区(以下简称为“岷山片区”)为研究区, 筛选气候、地形地貌、生物环境、人为干扰4类16项环境变量, 应用MaxEnt模型探讨影响大熊猫分布的主要环境变量, 并模拟现在(2004—2020年), 预测未来(2050s)低强迫(SSP1-2.6)和高强迫(SSP5-8.5)情景下的大熊猫生境, 分析气候变化下大熊猫生境适宜性的时空演变特征及趋势, 为区域内大熊猫的气候适应性保护和国家公园管理提出建设性建议。

1 材料与方法 1.1 研究区概况

《大熊猫国家公园体制试点实施方案(2017—2020年)》(以下简称“《试点实施方案》”)^[29]中将大熊猫国家公园分为岷山片区、邛崃山-大相岭片区、秦岭片区和白水江片区四个片区。岷山片区位于岷山山系中部, 行政区划涉及成都、德阳、绵阳、广元、阿坝藏族羌族自治州等5个市(州), 汶川、茂县、平武、九寨沟、松潘、青川等12个县(市、区), 总面积10013km², 占四川省内大熊猫国家公园面积(20177km²)的49.63%, 地理坐标为103°32′22″—105°4′38″E, 31°6′10″—33°7′56″N(图 1)。区域内地形复杂, 地势西高东低, 海拔664—5513m；气候垂直分异大, 年均温12—16℃；分布着5大类, 34个自然保护地类型, 濒危珍稀物种种类丰富。根据全国第四次大熊猫调查结果显示, 岷山片区共有大熊猫656只, 是大熊猫分布最多、适宜生境面积最大、遗传多样性最高的区域^[30]。

1.2 数据来源与预处理 1.2.1 大熊猫痕迹点

痕迹点数据来源于2015年发布的《四川省第四次大熊猫调查报告》^[31](以下简称“《四调报告》”)。通过ArcGIS进行数字化处理后得到461个大熊猫痕迹点。为避免空间自相关, 根据每只大熊猫一般需要3.89—6.40km²的活动领域^[32], 以1.20km作为半径对大熊猫痕迹点做进一步筛选, 最终得到326个大熊猫痕迹点用于MaxEnt模型计算。

1.2.2 环境变量

影响大熊猫生境的环境变量较多^[33], 而环境变量的数量、各变量的空间相关性和维度量会影响MaxEnt模型预测结果的准确性, 输入的环境变量过多、过复杂易导致模型的预测结果不可靠^[34]。基于大熊猫生态习性、相关文献总结、研究区实际情况和环境变量复杂度^{[16, 21]}, 筛选得到4类16项影响大熊猫生境的环境变量(表 1)。

气候数据来源于Worldclim2.1数据库(https://www.worldclim.org/)的19个生物气候变量(Bio1—Bio19), 其因具有较强的生物学意义已被广泛应用于物种生境的预测^[35]。现在气候时段选择数据集可提供的最新历史数据(1970—2000年), 未来气候情景选择2050年(2041—2060的平均值)的SSP1-2.6(2100年辐射强迫达到2.6W/m²)和SSP5-8.5(2100年辐射强迫达到8.5W/m²), 现在与未来气候数据的分辨率均为1km。为避免气候变量间的空间共线性影响, 采用相关性/贡献度法(CVC)和主成分分析法(PCA)排除相关系数大于0.80和贡献率小于1%的变量, 最终筛选得到7个变量。

地形地貌数据来源于地理空间数据云(https://www.gscloud.cn/)的DEM数据, 分辨率为30m, 由ArcGIS分析计算后得到海拔、坡度、坡向。土地利用数据来源于全球地表覆盖数据GlobeLand30(http://www.globallandcover.com/), 分辨率为30m。NDVI来源于中国资源环境科学与数据中心(https://www.resdc.cn/)的年度NDVI分布数据集, 分辨率为1km。水系、道路数据来源于全国地理信息资源目录服务系统(https://www.webmap.cn/)的1 ∶ 25万全国基础地理数据库。矿山和水电站数据来源于《四调报告》^[29]。除气候变量外, 其余环境变量数据年份均为2015年。

在ArcGIS 10.2中, 统一所有数据的坐标系为WGS_1984_UTM_Zone_48N, 分辨率为30m, 基于研究区边界进行裁剪, 并对水系、道路、矿山和水电站数据进行欧式距离分析, 最终输出MaxEnt模型所需要的ASCII格式文件。

1.3 研究方法 1.3.1 MaxEnt模型构建及评估

应用MaxEnt3.4.4进行建模, 输入大熊猫痕迹点及环境变量数据, 随机选择75%的痕迹点构建模型, 剩余25%的痕迹点用于验证模型精度。以受试者工作特征曲线下面积(area under curve, AUC)来评估模型结果^[36], AUC值的取值范围在0.5—1之间, 越接近1模型结果精度越高, 一般认为AUC值大于0.70即模型达到可接受的性能^[37]。构建MaxEnt模型需调整参数以达到更精确的模拟结果, 其中特征函数组合和正则化系数是优化模型的两个重要参数^[38]。模型提供5种特征函数：线性特征(linear feature, L)、二次特征(quadratic feature, Q)、乘积特征(product feature, P)、阈值特征(threshold feature, T)、合页特征(hinge feature, H), 所选函数越多模型越复杂^[39]。研究以AUC值为指示, 多次试验后选择LQPH特征函数模式组合, 正则化参数设置为2, 其余参数设为默认, 并选择刀切法(Jackknife)结合模型自行输出的环境贡献率以探讨各环境变量的重要性。

1.3.2 生境适宜性划分及评估

MaxEnt模型输出的结果是基于物种的存在概率逻辑值(P)生成的连续栅格图层, P值范围为0—1, 越接近1, 表示该物种在该区域的适生概率越大, 通过设定阈值可将连续结果转化为离散的适宜和不适宜的分类结果^[40]。MaxEnt模型的11种阈值判定标准中, TPT平衡阈值(Balance training omission, predicted area and threshold value)为平衡训练遗漏、预测面积和阈值基础之上一个较为保守的二分类阈值, 能确保最小的遗漏率^[40], 因此选择TPT平衡阈值进行适宜生境和不适宜生境的划分, 再利用ArcGIS的自然断点法将适宜生境划分为高、中、低三个等级。现在气候、SSP1-2.6情景、SSP5-8.5情景的TPT平衡阈值分别为0.15、0.12、0.10(表 2)。因大熊猫生境空间分布不规则, 研究以质心定义适宜生境的中心点来表征大熊猫适宜生境的迁移趋势^[41], 并通过统计大熊猫痕迹点在不同气候情景及适宜性等级生境中的出现频率, 以分析气候变化对大熊猫活动范围的影响。

2 结果与分析 2.1 模型结果评估

现在气候、SSP1-2.6情景和SSP5-8.5情景的测试集AUC最高值分别为0.81、0.80、0.83, 与之对应的训练集AUC值分别为0.80、0.83、0.81, 均大于0.80, 说明模型预测结果良好。

2.2 影响大熊猫分布的环境变量

环境变量贡献率和刀切法检验结果显示, 贡献率最高的前3项依次为最暖月最高温(Bio5)、海拔、最暖季度降水量(Bio18), 前10项的累计贡献率分别为96.90%。刀切法结果前3项为最暖月最高温(Bio5)、海拔、最暖季度降水量(Bio18), 其次为降水量季节性变化方差(Bio15)和距道路距离, 总体与贡献率结果基本相同。各类环境变量中, 气候类变量的累计贡献率为71.20%, 其次为地形地貌类(16.10%)、人为干扰类(8.60%)、生物环境类(4.10%), 且贡献率前10项中属于气候类变量的有5项(图 2), 表明气候类变量在MaxEnt模型模拟和预测研究区大熊猫生境方面起着重要作用, 通过未来气候数据预测大熊猫生境变化较可靠。

2.3 气候变化下大熊猫生境适宜性的时空变化 2.3.1 大熊猫生境现状

现在气候下, 大熊猫总适宜生境面积为6797.21km², 占研究区总面积的67.88%。其中高适宜占比25.84%, 中适宜19.60%, 低适宜25.84%。高适宜生境主要位于研究区中部及东北部, 包括平武县、北川羌族自治县及青川县西部；中、低适宜生境主要位于研究区南部及北部, 包括茂县东部, 安州区、绵竹市、什邡市、彭州市、都江堰市北部, 九寨沟县南部；不适宜生境占比32.12%, 主要位于西北部、东部和南部的研究区边界附近(图 3)。综上表明, 岷山片区目前整体生境质量较好, 但不适宜生境的比例也不容忽视, 保护工作仍有待进一步加强。

2.3.2 气候变化对大熊猫生境适宜性的影响

大熊猫生境面积变化结果显示, 与现在气候相比, SSP1-2.6情景下, 岷山片区大熊猫总适宜生境面积为6659.22 km², 占比66.51%, 缩减率为2.03%, 其中高、中适宜生境面积缩减率依次为32.66%、0.56%, 低适宜生境面积有所增加, 增加率为31.94%。SSP5-8.5情景下, 研究区大熊猫总适宜生境面积为4504.66km², 占比44.99%, 缩减率为33.73%, 其中高、中、低适宜生境面积缩减率依次为58.77%、33.85%、4.55%(图 4)。总体上, 未来气候情景下的研究区大熊猫的适宜生境面积呈下降趋势, SSP5-8.5情景下气候对大熊猫生境造成的影响更大。

大熊猫生境空间分布变化结果显示, SSP1-2.6情景下, 适宜生境丧失区域主要位于研究区东部, 包括平武县中部、北川羌族自治县西北部, 茂县东部；SSP5-8.5情景下, 除东部外, 研究区北部和南部也有大量生境丧失, 包括北川羌族自治县西部、茂县东部、绵竹市、什邡市北部及都江堰市北部。未来两种气候情景下岷山片区西北部均有新增生境, 主要包括松潘县东部、平武县西北部及彭州市北部(图 5)。

质心变化结果显示, 相较现在气候, 未来大熊猫适宜生境质心(32° 16′ 58″N, 104° 7′ 55″E)均向西北方向移动, SSP5-8.5情景(32° 19′ 14″N, 104°4′ 23″E)的变动范围相较SSP1-2.6情景(32° 18′ 9″N, 104° 5′ 28″E)更大。同一气候情景下, 不同适宜等级生境质心的迁移方向也存在差异。SSP1-2.6情景下, 各等级适宜生境的质心呈向研究区中部靠近的趋势, 其中高适宜生境质心向西南方向迁移, 中、低适宜生境质心向东北方向迁移；SSP5-8.5情景下, 各等级适宜生境的质心都有一定程度的北移, 其中高、低适宜生境质心向西北方向迁移, 中适宜生境质心向东北方向迁移(图 6)。总体来看, 大熊猫适宜生境呈向西、西北部扩展的趋势。

痕迹点出现频率结果显示, 现在气候下, 大熊猫痕迹点与适宜生境基本重合, 痕迹点在适宜生境出现的频率为97.85%, 其中高、中适宜生境出现的频率分别为67.79%、24.23%, 表明大熊猫更倾向于在高、中适宜生境中活动。SSP1-2.6情景下, 痕迹点仍有大部分位于适宜生境内, 出现频率为91.41%, 其中出现在高适宜生境的频率减少了22.70%, 中、低适宜生境的频率分别增加了3.38%、12.88%。SSP5-8.5情景下, 大部分痕迹点超出了高、中适宜生境范围, 出现在适宜生境的频率下降至82.81%, 其中痕迹点出现在高、中适宜生境的频率分别减少了43.56%和2.76%, 出现在低适宜生境和不适宜生境的频率增加了28.22%和15.64%(图 7)。综上表明, 未来气候情景下显著的适宜生境空间分布变化与质量下降对大熊猫的活动范围有较大限制。

3 讨论 3.1 结果准确性评估

模型结果精度验证表明, 通过MaxEnt模型对大熊猫生境进行模拟和预测具有较高的准确性。在构建未来大熊猫生境预测模型时, 除气候数据外, 其他环境数据均与现在气候保持一致, 理由如下：首先, MaxEnt模型的优势为即使物种数据分布信息或分布区的环境变量不完整, 也能较为精准地预测分析未来物种的潜在分布区^[42]；其次, 研究结果显示气候类变量对MaxEnt模型预测的贡献率远大于其他变量类型；同时, 大熊猫国家公园实行非常严格的保护, 其他环境变量的变化均受到一定程度的控制, 因此采用该种方法具有一定的合理性, 研究结果也较可靠。

3.2 环境变量对岷山地区大熊猫分布的影响

气候类变量对岷山片区大熊猫分布的综合影响程度最大, 与已有结果相符^{[16, 43]}, 但相关研究认为距河流的距离^[14]、降雨量变化方差^[16]、最高温^[44]等为最主要的影响项。研究认为主要有三方面原因：一是已有研究多采用发布较早的CMIP5数据作为气候变量源数据, 且忽略了人为干扰的影响, 在变量自相关性的排除标准上也存在差异^{[16, 21]}。二是MaxEnt模型在默认设置时的模拟效果不佳, 需对参数进行调整, 目前较多研究都未调整参数, 复杂的函数关系易导致过拟合从而影响结果^[43]。三是已有研究表明气候变化对不同山系、地域的影响程度不同^[44], 如晏婷婷等^[19]表明气候变化对邛崃山系的影响较小。本研究的气候类变量采用CMIP6数据, 并将人为干扰纳入环境变量综合分析, 调整MaxEnt模型参数至精度最佳, 得到影响岷山片区大熊猫分布的主要环境变量为最热月最高温(Bio5)、海拔、最暖季度降水量(Bio18)和距道路距离。

3.3 气候变化下岷山片区大熊猫生境的时空变化特征

气候变化对岷山片区大熊猫生境的影响主要体现在面积、质量和空间分布上。生境面积上, 未来气候SSP1-2.6和SSP5-8.5情景下岷山片区的大熊猫适宜生境面积均有所减少, 但由于气候情景、模型选择和研究范围的影响, 已有研究关于大熊猫生境的面积减少比例存在差异^[20]。本研究通过选择最新气候模式中的“最佳情景”和“最坏情景”进行生境预测, 可提供气候变化下2050s岷山片区大熊猫适宜生境面积的变化范围。生境质量上, 未来大量高、中适宜生境向低适宜和不适宜生境转化, 整体质量呈下降趋势, 其中新增生境的现有植被以高山灌丛、草甸为主, 较难演替成大熊猫生存所需的森林^[45], 南部新增生境区域有较多耕地、草地分割林地斑块, 因此未来新增生境以低、中适宜性为主。空间分布上, 适宜生境呈向西、西北方向高海拔地区移动的趋势, 与刘艳萍^[44]和Liu J等^[46]研究结果相符, 一方面是未来岷山片区整体气候变得干暖, 高纬度、高海拔地区更适宜大熊猫种群生存, 另一方面研究认为东部适宜生境的大量减少也是迫使大熊猫质心迁移的主要原因之一。

3.4 未来大熊猫生境保护建议

综上, 当前大熊猫国家公园建设的保护工作已初见成效, 但仍需加强对气候变化下大熊猫生境丧失和新增区的关注, 研究基于结果并结合《试点实施方案》、《大熊猫国家公园总体规划(2023—2030年)》(征求意见稿)(以下简称“《总体规划(意见稿)》”)和《四川省大熊猫国家公园管理办法》提出气候变化下大熊猫生境的保护建议：(1)识别优先保护区。优先保护区的确认有利于提高保护效率^[47], 岷山片区未来适宜生境丧失区主要涉及王朗、九寨沟、唐家河、宝顶沟和龙溪虹口等自然保护区, 需进行优先保护；小寨子沟和雪宝顶自然保护区生境相对稳定, 可作为大熊猫应对气候危机的避难所。(2)优化生态网络。构建景观生态网络是保护野生动物和提高区域气候适应能力的有效途径^[48], 目前岷山片区已规划3条大熊猫迁徙廊道以沟通主要的局域种群, 分别为土地岭、施家堡和黄土梁廊道, 气候变化下3条廊道附近区域均有大量生境丧失, 需加强对廊道连通性的评估与保护；开展对西北部、南部新增生境的资源勘探工作, 筛选适宜的生境进行抚育以构建网络；同时大熊猫具有季节性迁徙的习性^[49], 未来可进一步对种群内部迁徙廊道进行系统规划。(3)调整管控分区。《总体规划(意见稿)》中将大熊猫国家公园分为核心保护区和一般保护区, 为应对未来气候变化, 可对分区进行优化调整。核心保护区除基本覆盖现有的高、中适宜生境区和局域种群的聚集区外, 可将未来适宜生境的新增区并入管理, 区域内进一步限制采矿、放牧等人类活动, 控制新建工程的数量及其对生态环境的干扰, 并需建立动态、全方位的气候及生境质量监测体系；一般控制区主要覆盖靠近人类聚落、道路的低适宜及不适宜生境区, 可结合区域保护现状进行一定程度的生态开发, 促进区域“自然-经济-社会”综合协调发展。

4 结论

研究基于MaxEnt模型得到影响岷山片区大熊猫分布的主要环境变量为最暖月最高温(Bio5)、海拔、最暖季度降水量(Bio18)和距道路距离, 并通过模拟和预测得到气候变化下大熊猫生境适宜性的时空演变趋势：未来大熊猫适宜生境面积、质量整体呈下降趋势, 适宜生境质心向西、西北方向的高海拔区域移动, 大熊猫活动范围逐渐缩小。因数据源及模型模拟精度的限制, 本文仍存在不足, 如现在气候数据的时间与定义时间有偏差；同时可食竹是影响大熊猫生境选择的重要环境变量, 但因其受气候变化影响的不确定性及对模型精度的影响, 未考虑将其代入模型, 后续研究可着眼于数据精度提升, 环境变量优化及进一步探讨气候变化对大熊猫生理的影响机制, 构建更完善的物种预测模型。