丹顶鹤(Grus japonensis)为鹤科(Gruidae)、鹤属(Grus)中的一类大型涉禽, 随季节不同周期性进行迁徙, 主要分布在中国东北、蒙古东部、俄罗斯乌苏里江东岸、朝鲜、韩国以及日本北海道^[1—3]。丹顶鹤对生存环境要求极为苛刻, 不仅是湿地生物多样性保护的旗舰物种, 更是湿地生态环境是否洁净、安全的最敏感和最明显的指示者^[4—6]。近年来, 随着人类活动加剧, 全球气候变暖, 湿地生态环境不断遭到破坏, 造成了丹顶鹤迁徙生境的丧失和破碎化, 严重影响到丹顶鹤栖息地的选择与分布^[7—8]。气候变化如何影响野生动物迁徙生境的分布范围已经成为当前保护生物学界关注的热点问题^[9]。

物种分布模型(Species Distribution Models, SDMs)是研究丹顶鹤潜在生境分布的有效方法, 主要通过利用物种的样本信息和对应的环境变量信息进行关联得出物种的分布与环境变量之间的关系, 并将其应用于所研究区内对目标物种的分布进行估计^[10]。随着生态位理论的完善和统计学模型的不断发展, SDMs经历了由浅入深逐渐完善的发展过程^[11]。目前, 主要使用的物种分布模型有广义线性模型(GLM)、推进式回归树模型(GBM)、分类与回归树模型(CTA)、人工神经网络(ANN)、表面分布区分室模型(SRE)、柔性判别分析(FDA)、随机森林(RF)、最大熵值模型(MaxEnt)等^[12]。然而仅使用单一模型对物种潜在分布区进行预测时往往稳定性不强、偏差较大, 而基于多模型构建的组合模型在准确度上的表现则相对更佳^[13]。W. Thuiller提出BIOMOD计算框架^[14], 能够集成10种物种分布模型(SRE、GLM、GAM、CTA、ANN、GBM、RF、FDA、MARS、MaxEnt), 在物种分布的模拟和预测领域发挥了重要作用^[15—17]。物种分布受到多种因素的影响, 如环境、地理和生态因素等, 仅使用单一模型对珍稀水禽迁徙路径上物种潜在分布区进行模拟和预测时往往稳定性不强、偏差较大。本文使用的组合模型可以降低因为对环境变量与生境之间相互关系判断错误的误差, 综合多种生境因子的组合情况以提高预测结果的鲁棒性, 使整体预测结果更加可靠。

目前丹顶鹤的生境模拟研究多采用数据驱动的MaxEnt模型, 研究范围往往围绕丹顶鹤的繁殖地和越冬地等重要保护区开展生境的适宜性评价^[18—21], 很少通过多模型的集成来模拟迁徙路径上潜在生境的分布范围及其生境质量^[20], 气候变化对于亚洲东部丹顶鹤迁徙路线生境质量的影响有待深入的研究。使用组合模型框架来对不同的模型进行评估并且基于评价指标和专家知识选择最优的模型要比使用单一模型进行物种分布更加可靠^[14]。本文基于气候变化背景下的138个丹顶鹤样本数据, 进行模型筛选后, 利用BIOMOD2建模平台中的9种模型算法构建组合模型模拟丹顶鹤在当前气候条件下的分布并进一步预测其在未来气候中的潜在分布情况, 分析丹顶鹤在中国适宜性分布与环境变量之间的联系及潜在适宜区的整体变化趋势, 为开展丹顶鹤栖息地合理保护管理提供了科学依据, 对保护生态资源和促进丹顶鹤野生种群健康, 维护区域生物多样性具有重要意义。

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

丹顶鹤分布范围包括黑龙江、内蒙古、辽宁、河北、山东、江苏分布区。为尽可能全面分析丹顶鹤南迁时节的潜在分布区, 并精确预测丹顶鹤停歇地, 本文研究区涵盖丹顶鹤样本分布区的相邻省份, 包括吉林、山西、陕西、北京、天津、河南、安徽省份(27°787′—54°787′N, 97°179′—136°679′E)。该区域海拔自东向西逐渐升高(-57—3654m), 东部平原地带主要为亚热带大陆性季风气候区和温带大陆性季风气候区, 季风性显著, 夏季高温多雨, 随纬度升高最冷月平均温度变化显著。西部地区平均海拔较东部高, 为温带大陆性气候区, 干燥少雨。全区域海拔由低向高土地覆被变化明显, 中部广阔平原遍布耕作和管理的农业用地, 由南向北从常绿阔叶林向落叶针叶林过渡, 由东向西植被覆盖度逐渐减少。

1.2 数据来源及预处理

准确的物种分布数据是构建物种分布模型的基础, 丹顶鹤的存在点位数据获取自中国观鸟记录中心(http://www.birdreport.cn/)、全球生物多样性信息服务网络平台(Global Biodiversity Information Facility, GBIF)(http://www.gbif.org/)及公开发表的文献^[22—26], 为保证数据的精确性和模型构建的科学性, 本研究筛选掉重复、地点有误的点位数据, 仅保留南迁时节(10—11月)有精确位置信息的数据, 最终保留138个丹顶鹤分布样点(图 1)。

完备且准确的环境因子是物种分布模型构建的基本要求^[27], 本文选择应用最为广泛的五种环境变量即气候变量(Bioclimatic variables系列数据包括19个具有生态学意义的气候变量^[28])、地形变量、植被变量、水文变量、人为干扰变量构建物种分布模型(表 1)。地形变量(海拔、坡度、坡向)由地理空间数据云(https://www.gscloud.cn/)提供的数字高程模型(DEM)获得。归一化植被指数(Normalized Difference Vegetation Index, NDVI)通过GEE(Google Earth Engine, GEE)云平台提取所得。水文、人为干扰变量通过对河流、道路矢量数据(来源于中国国家基础地理信息数据平台)进行距离制图获得。

同时本文设定的地形、植被、水文、人为等因子在未来SSP1.2-6气候情境下不发生改变, 仅研究气候变化对丹顶鹤未来分布格局的影响。本文当前(1970—2000年)和未来(2021—2100年)气候变量均来源于WorldClim数据库(http://www.worldclim.org/)。未来环境变量采用2021—2040年(2030s)、2041—2060(2050s)、2061—2080(2070s)、2081—2100(2090s)四个阶段的SSP1.2-6气候情景数据。为避免环境变量间的相关性过高导致模型过拟合^[29], 采用R软件包进行相关系数(Person)分析, 对相关系数>0.8的多组环境变量进行筛选, 每组仅选取一个与丹顶鹤生境特征关系关联紧密的变量参与建模, 并检查筛选所得的环境变量的相关系数(r＜0.8)^[30]。最终筛选得到Bio 1年平均气温、Bio 2平均气温日较差、Bio 3等温性、Bio 4气温季节性变动系数、Bio 13最湿月降水量、Bio 14最干月降水量、ASL海拔、FVC植被覆盖度、Disriver离河流距离、Disroad离道路距离共10个环境变量参与丹顶鹤物种分布模型建模(表 1)。

1.3 基于BIOMOD 2的物种分布单模型构建

本研究依据物种的存在点位数据, 采用BIOMOD 2算法中可成功运行的9个单模型进行建模。由于缺乏缺失样本, 因此在建模过程中随机生成5组伪缺失样本, 每组各包含500个伪缺失点位数据。并将参与建模的样本随机分为80%的训练样本和20%的验证样本, 前者作为训练集, 后者用于精度评估, 将每个单模型运算算法重复10次。最后使用接受者操作特征曲线的曲线下面积(area under ROC curve, AUC)、真实技巧统计值(true skill statistics, TSS)及Kappa系数对9个模型进行拟合精度的评价。

1.4 基于BIOMOD 2的物种分布组合模型的构建

本研究采用组合物种分布模型来预测丹顶鹤适宜生境的空间分布格局, 以此来减少物种分布模型运算过程及伪样本点生成过程中的不确定性^[31]。综合考虑BIOMOD 2所包括的9种物种分布模型的所有450个单模型运算结果(仅保留AUC>0.85的单模型运算结果), 采用各模型的AUC值作为权重进行加权平均, 构建组合物种分布模型。具体构建过程：①将9种单模型运算结果进行归一化处理, 即将每个模型的运算结果除以其最大值。②利用每种单模型的归一化AUC均值与所有9种模型的归一化AUC均值之和的比值作为各模型的权重(式1)。③将9种单模型的归一化AUC均值乘以相应的权重, 并求和, 以此构建组合物种分布模型(式2)。

(1)

式中, w_j是第j个模型的权重; r_j是第j个单模型的归一化AUC均值; h表示AUC>0.7的单模型运算结果的个数。

(2)

式中, y_i是第i个栅格的丹顶鹤适宜生境指数; w_j是第j个模型的权重; x_ij是第j个模型中的第i个栅格的数值。y_i的取值范围[0, 1], y_i越接近1, 表明该栅格单元的环境条件越适宜丹顶鹤停歇^[27]。

最终利用ArcGIS 10.8对输出的组合模型结果进行适宜性分级, 划分为5个适宜性等级：0—0.2为不适宜分布区、0.2—0.4为低适宜分布区、0.4—0.6为中适宜分布区、0.6—0.8为高适宜分布区、0.8—1为最适宜分布区^[32]。

1.5 模型准确性评价

本文采用AUC、TSS及Kappa系数等3种评估指标来评估模型的准确性^[33]。AUC为接受者操作特征曲线(receiver operating characteristic, ROC)的曲线下面积, 不受分布点发生率(prevalence)和判断阈值的影响^[34]。TSS表示预测结果区分“是”和“否”的能力, 不依赖于分布点发生率, 但是受到阈值的影响^[29], 其计算公式如下：

(3)

式中, TPR为真阳性率(true positive rate), FPR为假阳性率(false positive rate), TSS值∈(-1, 1), 值越接近1表明准确性越高^[35]。

Kappa系数指预测相对于随机发生的准确率, 受分布点发生率和阈值的影响^[29], 其计算公式如下：

(4)

式中, (expectedcorrect)_ram=1/N[(TP+FN)(TP+FP)+(TN+FN)(TN+FP)], TP、TN、FP、FN分别为真阳性、真阴性、假阳性、假阴性计数值^[29]。

2 结果与分析 2.1 丹顶鹤物种分布模型精度评估

基于138个正样本数据模拟丹顶鹤的潜在适宜生境, 并依据AUC值、TSS值和Kappa系数评估单模型和组合模型的模拟精度(图 2)。研究结果表明：9种单模型的模拟精度均较好(AUC>0.7、TSS>0.5、Kappa>0.5), 其中RF随机森林算法在当下条件模拟精度最高(AUC>0.95、TSS>0.9、Kappa>0.85), 且GBM、MARS、MAXENT、FDA及GLM均具有较高模拟精度(AUC>0.9、TSS>0.8、Kappa>0.8), 说明单模型组合构建的模型能较好地预测丹顶鹤的潜在适宜生境。集成9种单模型所形成的组合模型AUC值、TSS值及Kappa系数分别为0.99、0.95、0.96, 较单模型预测精度更高, 说明组合模型能更精确地模拟丹顶鹤的潜在适宜生境。

2.2 物种分布模型的环境变量重要性评估

不同模型中环境变量的贡献率如表 3所示, 依据各个环境变量的贡献率获取影响丹顶鹤分布的主要变量, 相对贡献率由高到低依次是海拔、平均气温日较差、年平均气温、最干月降水量、等温性、气温季节性变动系数、最湿月降水量、植被覆盖度、离道路距离和离河流距离。其中平均气温日较差和海拔是影响丹顶鹤分布的主导因子。

2.3 丹顶鹤当前潜在生境分布格局模拟

结合BIOMOD 2组合模型按相等间隔分类的预测结果(图 3)和丹顶鹤潜在适宜生境的面积统计图(图 4)分析可知, 当前丹顶鹤潜在适宜生境即可作为丹顶鹤迁徙停歇地的面积为3.02×10⁵km²。主要分布在黑龙江省扎龙湿地保护区、黑龙江下游沿岸地区、环渤海地区、山东省黄河三角洲自然保护区、江苏省盐城湿地保护区及河南、安徽等地。其中, 最适宜生境面积为0.29×10⁵km², 占总适宜生境面积的9.41%, 主要分布在山东省西北部、环渤海沿海地区及江苏省盐城湿地保护区。高适宜生境面积为0.25×10⁵km², 占总适宜生境面积的8.11%, 主要分布在黑龙江省扎龙湿地保护区、吉林省向海湿地保护区、环渤海内陆地区及江苏、安徽等地。中适宜生境面积为1.30×10⁵km², 占总适宜生境面积的42.2%, 主要分布在黑龙江下游沿岸地区及辽宁、河北、山东、河南、安徽、江苏等内陆地区。丹顶鹤的潜在适宜生境在吉林和辽宁之间出现了明显的不连续现象, 但在其余地区连通性较好。

2.4 气候变化对未来丹顶鹤潜在生境适宜性的影响

本研究采用相同气候背景(SPP1-2.6)下的环境指标(本文设定的地形、植被、水文、人为等因子在未来SSP1.2-6气候情境下不发生改变), 通过构建组合模型预测未来4个不同时期的丹顶鹤潜在适宜生境(图 5), 分别将未来预测结果和当前预测结果根据适宜和不适宜进行二值化, 由此得到不同时期丹顶鹤潜在适宜生境的面积变化(图 6、图 7、表 4)分析丹顶鹤未来分布格局。组合模型预测结果表明, 丹顶鹤未来的潜在适宜生境分布格局和当前分布格局相近, 但潜在适宜生境面积缩减幅度较大, 相对当前分布格局的范围变化为-8.94%—-14.39%。在2021—2040阶段, 潜在适宜生境总面积为2.60×10⁵km², 较当前潜在适宜生境总面积缩减了0.65×10⁵km², 缩减率达21.67%, 缩减面积显著。主要原因可能是该时期全球变暖趋势加剧, 导致丹顶鹤适宜生境面积急剧减少。在2041—2060阶段, 潜在适宜生境总面积为2.58×10⁵km², 较当前生境总面积缩减了0.69×10⁵km², 扩增了0.25×10⁵km², 面积缩减趋势仍较为明显。在2061—2080阶段, 潜在适宜生境总面积为2.75×10⁵km², 较当前生境总面积缩减了0.51×10⁵km², 扩增了0.24×10⁵km², 扩增率为7.98%, 缩减率为16.93%。该阶段丹顶鹤适宜生境面积缩减率最小, 可能与人类为保护环境所作的相关举措息息相关, 为丹顶鹤营造了适宜的生存环境。在2081—2100阶段, 潜在适宜生境总面积为2.60×10⁵km², 较当前生境总面积缩减了0.59×10⁵km², 扩增了0.17×10⁵km²。缩减和扩增面积相较之变化幅度较为稳定。

总体上, 丹顶鹤未来分布格局与当前分布格局相似, 未来潜在适宜生境面积明显缩小的区域主要分布在胶东半岛和环渤海区域, 最、高、中适宜性栖息地面积均呈减少的趋势, 缩减幅度大于扩增幅度。

3 讨论

气候变化正在加速影响大型涉禽类候鸟的迁徙与分布格局, 其中最明显的表现形式即繁殖地与迁徙停歇地分布范围的缩减^[36]。有研究表明气候变化会导致东亚地区丹顶鹤到2050年丧失约30%的繁殖地适生区^[19]。本研究发现气温日较差和海拔是影响丹顶鹤迁徙停歇地生境适宜性的主导因子。受全球变暖的影响, 研究区的平均气温日较差逐年增加。日平均最高气温的增加导致地表潜在蒸发量增大, 湿地水资源蓄积量减少, 适栖生境面积进一步萎缩。日平均最低气温的下降和极端低温的出现, 会使得部分丹顶鹤的食物种类在寒冷的天气期间变得不那么容易获得, 造成丹顶鹤取食困难, 给迁徙期物种生存带来严峻的挑战^[37]。在高海拔地形的阻挡作用下, 丹顶鹤的迁徙廊道也受到了严重的影响。受山东丘陵和长白山脉的地形影响, 丹顶鹤迁徙地的适宜生境在山东省中部和南部、吉林省南部与辽宁省北部之间出现了明显的不连续现象。总体来看, 丹顶鹤的分布格局未来不会发生较大变化, 环渤海地区和黄河三角洲地区仍然是丹顶鹤主要的南迁的停歇地, 山东省中部和南部、吉林省与辽宁省之间栖息地的间断未来不会消失。但潜在适宜生境面积缩减幅度较大, 相对当前分布格局的范围变化为-14.39%—-8.94%, 高度适宜栖息地和中度适宜栖息地的面积会持续减少, 维持现有适宜栖息地和恢复原来的栖息地比较困难, 栖息地对气候的变化抵抗力较弱。在SSP1.2-6情境下全球的CO₂按照严格的标准削减, 在2050年达到了零排放, 但在如此理想的条件下丹顶鹤的栖息地适宜面积仍然有缩减, 因此未来丹顶鹤适宜栖息地的变化不容乐观, 未来的气候变化对丹顶鹤的生存构成一定的威胁。

丹顶鹤迁徙生境对温度极为敏感, 受气候变化的影响较大, 极端天气频发将会影响正常的生物迁徙节律。本研究发现在SSP1.2-6的气候情境下, 丹顶鹤秋季迁徙路线上适宜栖息地的面积整体呈收缩趋势, 但空间上的变化趋势表现出明显差异。受气候变暖影响, 齐齐哈尔的扎龙湿地北部、辽宁省胶东丘陵东缘、长白山脉南缘等部分地区的优质栖息地的面积将会少量增加。而部分湿地分布区的丹顶鹤适栖生境面积将会缩减, 如盐城市生物圈保护区南缘, 丹顶鹤的种群可能会迁出目前为其设计的保护区, 从而影响保护区的生境适宜性^[38]。由于丹顶鹤大陆种群的繁殖范围分布在中国、俄罗斯和蒙古的边界, 气候变化也可能导致一些种群北上俄罗斯和蒙古, 从而改变每个国家保护这一濒危物种的作用^[7]。因此, 如何协调好边境丹顶鹤跨境迁徙, 成为中国丹顶鹤保护的主要问题之一^[39]。

预测气候变化对野生动物分布格局潜在影响的关键技术之一即构建物种分布模型。与单一模型相比本文使用的组合模型可以降低因为对环境变量与生境之间相互关系判断错误的误差。如表 3所示, GBW和GTA模型并未检测出气温日较差对于丹顶鹤生境适宜性的贡献率, RF和GTA模型过低的预测了海拔对于丹顶鹤迁徙生境的影响。而BIOMOD组合模型同其它单一模型均指示气温日较差和海拔是影响丹顶鹤迁徙生境适宜性的最重要的因子。这说明BIOMOD组合模型为多模型结果的整合的综合, 可有效避免单一模型不足, 提高预测结果的准确率, 减少由单一建模方法引起的偏差。当使用的物种分布数据和模拟重复次数较多时, 组合模型精度更高、更稳定^[33]。候鸟迁徙时期其适宜栖息地变化随温度, 初级生产力等变化较快。未来将针对秋季迁徙期的不同阶段, 进一步获取较高时间分辨率的环境因子, 以提高生境适宜性评价及模拟结果的可靠性。除本研究重点考虑的气候、地形因素外, 还包括隐蔽条件、食物来源和人为干扰程度等。因此, 在利用物种分布模型预测野生丹顶鹤潜在分布时, 应全面考虑上述因素。丹顶鹤主要以鱼、虾等水生动物以及水生植物为食, 对湿地空间的依赖导致其空间分布受气候和食物资源变化的双重影响^[40—41]。湿地多生长对环境变化敏感的湿生植物, 温度降水的剧烈变化会导致大量敏感型植被退化消失, 湿生植被逐步被演替, 进而改变湿地的生态功能^[42]。因此, 预测该物种在未来气候中的潜在分布情况, 对其生境的保护管理政策的制定与未来优先保护区的划定具有重要的参考意义。

4 结论

基于本文通过对多变量的遴选, 系统地模拟了在SSP1.2-6气候情境下丹顶鹤当下及未来的迁徙和栖息地分布情况, 并对丹顶鹤分布的区域及变化成因进行了分析, 得出主要结论如下：

(1) 不同物种分布模型中的气候变量的贡献率不同, 相对贡献率由高到低依次是海拔、平均气温日较差、年平均气温、最干月降水量、等温性、气温季节性变动系数、最湿月降水量、植被覆盖度、离道路距离和离河流距离, 其中平均气温日较差和海拔是影响丹顶鹤分布的主导因子。

(2) 当前丹顶鹤潜在适宜生境即可作为丹顶鹤迁徙停歇地的面积为3.02×10⁵km²。主要分布在黑龙江省扎龙湿地保护区、黑龙江下游沿岸地区、环渤海地区、山东省黄河三角洲自然保护区、江苏省盐城湿地保护区及河南、安徽等地。丹顶鹤的潜在适宜生境在吉林和辽宁之间出现了明显的不连续现象, 但在其余地区连通性较好。

(3) SSP1.2-6气候背景下, 相对当前分布格局的范围变化为-8.94%—-14.39%。丹顶鹤未来分布格局与当前分布格局相似, 未来潜在适宜生境总面积较当前明显缩减且缩减幅度较大, 最、高、中适宜生境面积变化均呈缩减趋势, 缩减幅度远大于扩增幅度。