气候变化背景下中国樟子松潜在分布预测

张世林，高润红，高明龙，韩淑敏，张文英，赵 静

（1.内蒙古农业大学 林学院，内蒙古 呼和浩特 010019；
2.包头市林业和草原工作站，内蒙古 包头 014030）

全球气候持续变暖会导致森林产生不同程度的响应。基于2050年气候变化情景，THOMAS 等[1]研究发现：中等气候变化背景下，覆盖地球表面20%样本区域将有15%～37%的物种面临灭绝。物种分布模型(SDM)常用来模拟物种潜在分布及适应生境变化[2]。目前，主要的物种分布模型有生物气候模型(Bioclim)、药用植物全球产地生态适宜性区划信息系统(GMPGIS)、气候动态分析模型(Climex)、规则集预测的遗传算法模型(GARP)和最大熵模型(MaxEnt)[3-5]，其中最大熵模型较其他同类模型具有预测结果准确性高、灵活性强等优点[6]。因此，最大熵模型在一些植物的潜在分布区研究中有广泛应用，如银杉Cathayaargyrophylla[7]、祁连圆柏Juniperusprzewalskii[8]、沙冬青Ammopiptanthusmongolicus[9]等。

樟子松Pinussylvestrisvar.mongolica作为欧洲赤松P.sylvestris的一个地理变种，属松科Pinaceae 松属Pinus，因其较强的耐寒、耐旱、耐贫瘠土壤和较速生等优良特性被广泛用于三北工程等人工林体系的建设[10-11]。早在1955年，樟子松作为防风固沙林首次成功引种到章古台。迄今为止，樟子松已先后在中国三北工程区的13 个省(市、自治区)590 余县引种成功。然而，年平均气温达约7 ℃时，樟子松人工林无法自然更新；
因章古台、毛乌素沙地和塞罕坝等地受降水限制[12-17]，樟子松出现生产力水平降低和枯死等现象。根据《三北工程六期规划(2021—2035年)》，樟子松仍是三北工程中最重要的常绿针叶造林树种。为防止樟子松引种后出现严重衰退或更新困难等情况，本研究利用优化后MaxEnt 模型预测不同气候环境下樟子松潜在分布，以期开展未来樟子松引种和优良种质资源培育工作。

1.1 樟子松数据来源

2018—2021年内蒙古自治区林业和草原监测规划院通过实地调查，在内蒙古自治区范围内樟子松分布区域共获取130 个分布数据；
结合中国数字植物标本馆(CVH)、中国国家标本植物平台(NSII)、全球生物多样性信息网(GBIF) 及其他相关文献共获取樟子松分布点数据115 个，经检索汇总后共保留245 条记录。将文献记载及实地调研数据整合，配合ArcGIS 软件对样点经纬度进行筛选剔除。为避免樟子松群体分布得不均匀，以经纬度2.5′×2.5′为1 个网格单元，每个单元筛选1 条信息记录，截至2021年最终保留200 个样点进行后续拟合分析。

1.2 环境因子变量及筛选

本研究选取空间分辨率均为30″的19 个生物气候因子和1 个地形因子作为环境因子(表1)。气候数据来源于世界气候数据库(https://www.worldclim.org)，地图数据来源于国家基础地理信息中心(http://www.ngcc.cn/ngcc/)的中国标准地图。

表1 全球气候数据库的环境因子Table 1 Environmental factors in the Global Climate Database

未来气候数据选取第6 次国际耦合模式计划(CMIP6) 的国家气候中心中等分辨率气候系统模式(BCC-CSM2-MR)。此模式为中国未来的气候模式。为更均衡合理模拟未来时期(2050s 和2100s)樟子松潜在分布，同时考虑社会经济共享路径(SSP) 和代表性浓度路径(RCP)，选择SSP-RCP2.6 (SSP126)、SSP-RCP4.5 (SSP245)和SSP-RCP8.5 (SSP585)等3 个情景对樟子松未来分布区进行预测。SSP126 情景使用可持续发展路径，代表温室气体排放量处于较低水平；
SSP245 情景使用中度发展路径，代表温室气体排放量处于中等水平；
SSP585 情景使用以化石燃料为主的发展路径，代表温室气体排放处于高水平[18]。

为避免20 个生物气候变量间的相似性过高导致模型过度拟合，通过方差膨胀因子(VIF)和Pearson 相关性检验，结合MaxEnt 模型预实验中各个因子贡献率，选择相关系数小于0.8、VIF 小于10 且对樟子松分布具有较高贡献率和生态学意义的因子。最终保留8 个气候因子和1 个地形因子进行建模，包括年平均气温、最热月最高气温、最冷月最低气温、年降水量、最冷季度平均气温、气温季节变动系数、最干季度平均气温、降水量季节性变化和海拔。

1.3 模型和方法

本研究利用MaxEnt 模型预测当前和未来时期樟子松在中国的潜在分布。将樟子松200 个分布点和9 个环境因子导入MaxEnt 3.4.1 软件，java 环境下运行，随机选取25%的样本作为测试数据集，75%的样本作为训练数据集[19]，以Bootstrap 模式，重复迭代10 次，建模结果以Logistic 形式输出。

MaxEnt 模型准确度可通过受试者工作特征曲线(ROC)进行检验。用该曲线下的面积(AUC)表示所建模型的精度，AUC 通常为0.5～1.0，其值越大，表明模型预测越精准。其中AUC 在0.5～0.6 时表明模型预测失败，0.6～0.7 时表明模型预测效果较差，0.7～0.8 时表明模型预测效果一般，0.8～0.9 时表明模型预测效果良好，0.9～1.0 时表明预测效果极好[20]，可以较准确地反映物种潜在分布。

1.4 模型优化

本模型采用优化方法为Checkerboard 2 法，调整调控倍率(RM)和特征组合(FC) 2 个参数改变模型正规化水平。通过改变特征数量、组合方式以及RM 数值，将MaxEnt 模型中的5 种特征：线性特征(L)、二次型特征(Q)、片段化特征(H)、乘积型特征(P)和阈值型特征(T)，与16 个RM 值(0.5～8.0，以0.5 为间隔梯度)排列组合出148 种组合方式。利用R 语言中的ENMeval 程序包对其进行优化测试，结果中各组合的ΔAICc 值和10%测试遗漏率越低，表示模型预测精准度越高[21-22]。

1.5 数据的分析和处理

运用ArcGIS 将模型运行后的数据进行适宜区划分和可视化处理，根据软件中提供的多种分类方法预测樟子松适宜性阈值。根据冉巧等[7]对孑遗植物银杉进行适宜区划分的方法，本研究使用自然断点法划分樟子松适宜指数。参照张伟萍等[8]对祁连圆柏适宜区划分的等级，将樟子松生境划分为不适宜区、较不适宜区、一般适宜区和高度适宜区(适宜指数分别为0～0.1、＞0.1～0.2、＞0.2～0.6、＞0.6～1.0)。对比樟子松不同时期的适宜区分布差异，获得樟子松不同适宜区空间分布图；
利用ArcGIS 10.4.1 软件中SDMtoolbox 2.0 工具包，计算樟子松在不同时期的适生区质心位置得到樟子松分布区的迁移方向，并计算其质心迁移距离。

2.1 模型优化及准确性分析

根据200 个樟子松分布点和9 个环境因子利用MaxEnt 模型对樟子松的潜在分布区域进行模拟预测。在MaxEnt 为默认参数设置时，RM 为1.0，FC 为LQHPT，ΔAICc 为29.77。当RM 为1.5，FC 为LQHPT 时，ΔAICc 为0，此时模型最优，10% 的训练遗漏率降低4.63% (表2)。因此选取RM 为1.5，FC 为LQHPT 为模型最终参数。该参数下10 次模拟训练AUC 均值为0.964 (图1)，说明预测结果精确。

图1 MaxEnt 模型的ROC 响应曲线Figure 1 ROC response curve under MaxEnt model

表2 不同参数MaxEnt 数模型优化结果Table 2 Environmental results of MaxEnt model under different parameter settings

2.2 当前中国樟子松分布

樟子松当前气候条件下的潜在适生区集中在53.55°N 以南，135.09°E 以西，总适生区面积(一般适生区面积与高度适生区面积之和)为645 104.2 km2，主要位于大兴安岭地区、内蒙古中部、黑龙江西北部，在山西、河北、吉林、辽宁也有少量分布[23]，其中，高度适生区面积为147 638.9 km2，一般适生区面积为497 465.3 km2。樟子松高度适应区绝大部分在内蒙古东北部，黑龙江西北部呈零星分布。

2.3 未来气候情景下樟子松分布预测

未来时期不同气候情景下，樟子松适生区面积较当前时期有不同程度减少。2050s，樟子松总适生区在SSP126、SSP245 和SSP585 气候情景下降35.7%、52.5%和82.6%；
2100s 与2050s 相比，总适生区面积在SSP126、SSP245 和SSP585 情景下呈下降趋势，降幅分别为62.1%、37.4% 和100.0%。SSP126-2050s 情景下，樟子松高度适生区面积相较当前大幅下降，面积约为当前的6.6%，仅存在大兴安岭较高海拔地区；
SSP245-2100s 气候情景下，樟子松高度适生区将完全消失。

在3 种气候条件下樟子松适生区面积均下降，SSP245 和SSP585 气候情景下樟子松适生区对气候的响应最为明显 (表3和表4)。SSP126 气候情景下，2050s 樟子松适生区面积变化幅度最小，适生区面积下降230 225.7 km2，降幅为35.7%，主要变化区域位于呼伦贝尔市西南部；
2100s 适生区面积下降487 968.8 km2，降幅为75.6%，主要在兴安盟和锡林郭勒盟的东北部出现较大面积消失。SSP245 气候情景下，2050s 樟子松适生区面积下降338 449.4 km2，降幅为52.5%，其变化区域主要位于内蒙古自治区呼伦贝尔市西东两侧地区和赤峰市；
2100s 其适生区面积下降453 211.8 km2，降幅为70.3%，内蒙古自治区东部和黑龙江省的大部分地区樟子松适生区消失。SSP585 气候情景下，樟子松适生区面积变化幅度最大，2050s 樟子松适生区面积下降532 882.0 km2，比当前减少82.6%，仅在呼伦贝尔市中部有樟子松分布；
2100s 樟子松适生区将在中国完全消失。

表3 中国不同气候情境下樟子松适生区面积Table 3 Suitable growing area of P. sylvestris var.mongolica under different climates scenarios in China

表4 不同时期樟子松适生区空间格局变化Table 4 Dynamic changes in the suitable area for P. sylvestris var.mongolica under different combination of climates cenarios

2.4 质心迁移

在空间分布格局方面，未来时期不同气候情景下樟子松适生区的迁移位置存在差异，但总体上呈现出向西北和西南方向迁移(表5)。当前时期樟子松适生区的质心在内蒙古兴安盟阿尔山市(47.09°N，120.25°E)。气候情景为SSP126-2100s 时，樟子松适生区向西北方向迁移，此时樟子松适生区质心位于内蒙古自治区呼伦贝尔市鄂温克族自治旗(48.41°N，120.17°E)，迁移距离为161.47 km；
气候情景为SSP585-2050s 时，樟子松向西北方向迁移，此时其适生区质心位于内蒙古自治区呼伦贝尔市新巴尔虎左旗(47.88°N，119.06°E)，迁移距离为126.36 km；
气候情景为SSP245-2100s 时，樟子松向西南方向迁移，此时其适生区质心位于内蒙古自治区锡林郭勒盟东乌珠穆沁旗(45.66°N，117.39°E)，迁移距离为270.80 km。未来时期不同气候变化情景下，全球增温增湿使得樟子松在中国适生区的质心整体向西北和西南迁移，质心随着适生区分布范围消减可能消失。

表5 不同气候情景下樟子松适生区质心变化Table 5 Core distributional shifts under different climate scenario/periods for P. sylvestris var. mongolica

2.5 环境因子分析

根据Jackknife 法分析结果(表6)，樟子松潜在地理分布贡献率超过5%的环境因子分别为最冷季度平均气温(22.1%)、降水量季节性变化(19.8%)、最冷月最低气温(18.1%)、气温季节变动系数(17.8%)、最干季度平均气温(9.0%)、最热月最高气温(6.1%)。这6 个环境因子总贡献率达92.9%。

表6 参与建模的环境因子贡献率及置换重要值Table 6 Percentage contribution and permutation importance of environment variables for P. sylvestris var.mongolica in the MaxEnt model

2.6 未来分布区生态特征

樟子松分布主要环境因子与潜在分布区的关系如表7所示。SSP126 气候情景下，2050s 和2100s 樟子松生境适宜度比当前降低了0.38 和0.54；
SSP245 气候情景下，2050s 和2100s 樟子松适宜度比当前降低了0.46 和0.55；
SSP585 气候情景下，樟子松生境适宜度急剧下降，2050s 下降了0.55，到了2100s 下降0.67，降幅为91%。

表7 樟子松适生区环境因子结果分析Table 7 Results analysis of major environmental variables in P. sylvestris var.mongolica suitable area

200 个樟子松分布点的年平均气温与其生境适宜度变化正好相反，2050s，SSP126、SSP245 和SSP585 情景下年平均气温分别增加19.5%、56.6% 和101.3%；
2100s 与2050s 相比，SSP126、SSP245和SSP585 情景下年平均气温分别增加15.9%、48.1%和90.7%。

200 个樟子松分布点的年降水量在SSP126、SSP245 和SSP585 情景下均为递增趋势，且年降水量与年平均气温变化趋势相同，都与樟子松物种适宜度变化相反。2050s，SSP126、SSP245 和SSP585 情景下年降水量分别增加5.67、20.32 和44.58 mm；
2100s 与2050s 相比，SSP126、SSP245 和SSP585 情景下年降水量分别增加29.58、18.27 和45.62 mm。

3.1 讨论

3.1.1 优化后模型对樟子松当前适生区预测 本研究基于气候和地形因子的研究，应用ENMeval 数据包中的Checkerboard 2 法来优化模型[23]。作为北方主要造林树种，樟子松集中分布在中国北方地区，具有喜光、耐寒、耐旱、耐贫瘠和抗病能力强的生长特性。在中国的大兴安岭北部，西起莫尔道戈、金河、根河，东到新林、呼玛以北有连续的分布，在南伊图里河的银河、免渡河、阿尔山等地呈带状或块状分布着樟子松纯林以及其与兴安落叶松Larixgmelinii构成的混交林[24]。自1955年樟子松被成功引种到辽宁省阜新市彰武县章古台以来，已推广到科尔沁沙地、浑善达克沙地、毛乌素沙地等地。本研究基于优化的最大熵模型预测结果显示：当前樟子松高度适生区也位于上述地区，预测结果与实际高度一致，表明经过优化后的MaxEnt 模型用于樟子松适生区的预测结果精准可靠。

3.1.2 环境变量对樟子松地理分布限制 Jackknife 法分析表明：年平均气温、最热月最高气温、最冷月最低气温、年降水量、最冷季度平均气温、气温季节变动系数、最干季度平均气温、降水量季节性变化和海拔是影响樟子松分布的主要环境因子。当樟子松存在概率大于0.6 时，对应的环境因子最适合植物生长。主导环境因子分析表明：气温和降水相关因子对樟子松分布贡献率达到98.3%，对其分布起决定性作用。这与赵晓彬[25]研究中樟子松适生于年平均气温-3.7～-1.5 ℃、年降水量375 mm 的结论相符。本研究预测樟子松的最适宜海拔为800 m 左右，与戴继先等[26]研究发现樟子松天然分布在300～900 m的阳坡或向阳坡的结论一致。水热条件一定程度上决定了樟子松潜在地理分布，地形在一定程度上也会对水热资源进行二次分配，因此，气温、降水和地形对樟子松分布格局的相互作用都不能忽视。

3.1.3 气候变化对樟子松地理分布的影响 2050s 和2100s 不同二氧化碳排放情景下，高温高湿使樟子松适生区不断减少。2100s 樟子松适生区显著减少，其降幅明显大于2050s。未来不同的气候情景下，樟子松适生区的空间分布变化随着气候变暖加剧，整体向西北高纬度和西南季节性降水量增加的地区收缩。SSP126 情景下，2050s 适生区质心迁移距离最短，说明气候小幅变暖对适生区分布范围影响不大。SSP245 情景下樟子松适生区质心点受季节性降水量影响整体向西南方向降水量充沛地区迁移，2100s 与2050s 相比其质心点迁移距离较大。SSP585 情景下温室气体排放增加，2050s 樟子松适生区质心点明显向西北高纬度地区迁移，2100s 适生区在中国消失。张日升等[27]研究表明：随着温室气体排放等环境问题加剧，年平均气温达7 ℃左右时樟子松将出现生长不良，导致大面积退化。吴祥云等[28]也曾提出：樟子松人工林引种到年平均气温为7.1 ℃的章古台会发生“早衰”现象，导致面积减退。未来气候情景下，温室气体排放，气温升高导致樟子松适生区面积逐步北移，在中国分布减少甚至至消失。综上所述，樟子松迁移趋势与CHEN[29]研究中国寒温带森林树种地理分布在未来全球变暖下往高纬度和降水量充沛地区迁移的趋势相符。

3.2 结论

本研究通过对MaxEnt 模型的优化，发现樟子松潜在分布最优模型的特征组合为片段化特征(LQHPT)、乘积型特征(P)和阈值性特征(H)，调控倍率为1.5。此模型复杂度低，预测结果较准确，预测樟子松适生区范围与实际分布情况基本一致。中国樟子松当前高度适生区集中在47.2°～53.3°N，118.8°～121.9°E 的大兴安岭北部高纬度地区。不同气候情景下，樟子松对水热变化最为敏感，特别是在年平均气温升高5.0 ℃左右时将导致其在中国的分布区发生显著变化。

通过樟子松不同适生区的划分，根据各自区划气候特点提出以下保护对策：①在不同地区选取生长良好的樟子松单株进行种质资源监测、保护、采集和培育，培育最适合当地气候条件的优良品种。②未来在不适合樟子松生长地区引种时，应避免引种樟子松。③提前进行系统的种源区划研究，培育出不同特性的樟子松并在本地区推广。应该在年平均气温为-3.7～1.6 ℃，降水量为375～425 mm 的东北等高纬度地区引种樟子松，以发挥其最大生态效益。

猜你喜欢适生区樟子松平均气温未来气候条件下当归适生区预测及时空变化分析中国中医药信息杂志(2023年6期)2023-06-13基于组合模型预测蛇足石杉在中国的适生区西北植物学报(2022年3期)2022-04-18气候变化下濒危植物半日花在中国的潜在分布植物资源与环境学报(2020年2期)2020-04-20乌兰县近38年气温变化特征分析新农民(2019年2期)2019-08-16巴拉圭瓜多竹适生区分布研究竹子学报(2018年2期)2018-10-25从全球气候变暖大背景看莱州市30a气温变化农家科技下旬刊(2017年12期)2018-04-16塞罕坝樟子松幼林抚育与管理现代园艺(2018年2期)2018-03-15初探北方樟子松栽培关键技术现代园艺(2018年2期)2018-03-15我眼中的樟子松北极光(2018年12期)2018-03-071981—2010年拐子湖地区气温变化特征及趋势分析现代农业科技(2018年22期)2018-01-15