落叶松(Larix)作为我国主要造林树种, 具有成林速度快、适应性强和水土保持效果好等优点。同时, 落叶松还在维持森林生态系统稳定以及发挥生态系统服务功能上具有重要作用^[1—3]。落叶松的分布受气候、地形等环境因子的影响较大, 随着气候变化的持续, 包括落叶松在内许多物种都在适生区的分布和面积上发生了变化。另外, 相比于天然林, 人工林树种单一、结构相对简单且稳定性差, 因此对气候变化更为敏感^[4—5]。在我国广泛种植的落叶松树种包括兴安落叶松(Larix gmelinii.)、长白落叶松(Larix olgensis.)、华北落叶松(Larix principis-rupprechtii.)和日本落叶松(Larix kaempferi.), 它们分布的地理范围跨度较大且生境各异, 适生区受气候变化的影响也不同。因此, 探测影响落叶松主要树种的环境因子, 揭示已有气候变化对落叶松人工林适生区分布的影响, 对提出积极、有效的林业应对气候变化策略管理具有十分重要的理论与现实意义。

物种分布模型是模拟不同气候情景下物种潜在适生区分布的主要工具, 通过分析物种分布与环境之间的关系, 以概率的形式表达物种对环境的喜好程度^[6]。主要物种分布模型包括生物气候模型(BIOCLIM)、最大熵模型(MaxEnt)和随机森林模型(RF)等^[7]。其中, MaxEnt基于物种分布信息与对应点位的环境数据, 通过最大熵原理和机器学习描述物种在生态空间上的适宜性, 从而模拟物种存在的地理区域^[8]。MaxEnt能较好地处理变量之间的交互关系, 并且在样本量较小的情况下依然保持较高的精度, 因而得到了广泛应用^[9]。刘攀峰等使用MaxEnt模型及空间分析技术研究了杜仲(Eucommia ulmoides)的生态习性及其在我国的潜在适生区, 并取得了可信度较高的结果^[10], 除此之外利用MaxEnt模型模拟蒙古扁桃(Amygdalus mongolica)、白栎(Quercus fabri)、胡杨(Populus euphratica)、马尾松(Pinus massoniana)、红松(Pinus koraiensis)、软枣猕猴桃(Actinidia arguta)、核桃(Juglans regia L.)等树种的潜在适生区时也取得了较好的效果^[11—17]。

目前有关落叶松人工林适生区与气候变化的研究还较少, 赵学鹏等分析了不同海拔高度下长白落叶松对气候变化的响应特征^[18]。杨志香等基于8个气候变量研究了兴安落叶松的气候适应性^[19]。张喜娟等预测并分析了兴安落叶松在未来气候情景下适生区范围及面积的变化^[20]。申家朋等对日本落叶松在我国的潜在适生区及其在未来情景下的变化开展了预测和分析^[21]。但上述研究都是面向未来气候变化情景下的预测分析, 还没有开展历史气候变化与落叶松适生区分布的研究。过去的几十年正是我国大规模进行植树造林的时期, 理解造林同时期气候变化对落叶松人工林在中国分布的影响, 有助于推动落叶松人工林应对气候变化的机理研究并为落叶松制定合理的种植及可持续经营方案提供科学依据。

综上, 本研究基于过去90年的气候数据, 通过MaxEnt模型模拟并分析了落叶松人工林在中国适生区的空间分布及其受气候变化影响发生的变化趋势。考虑到落叶松不同树种间的生境差异较大, 将落叶松分为兴安落叶松、长白落叶松、华北落叶松和日本落叶松4个树种, 然后按照先整体再分树种的方式开展模拟和分析。主要工作包括：(1)探测影响落叶松人工林分布的主要环境因子, 分析不同落叶松树种对主要环境因子的适应范围。(2)分析当前(1991—2020年)气候条件下落叶松人工林适生区的空间分布特征。(3)对比1931—1960年、1961—1990年以及当前气候条件下, 不同气候历史时期落叶松人工林适生区的空间分布与变化, 探究适生区分布对气候变化的响应机制。

1 数据和方法 1.1 物种分布数据

通过查阅中国数字植物标本馆(www.cvh.ac.cn)、物种多样性数据平台(http://www.especies.cn/)等, 收集落叶松主要树种在中国的分布点, 进而通过查找落叶松人工林相关文献, 根据文献中记录的落叶松人工林生长良好的区域, 筛选长势稳定的地点且保证与气候数据所对应的时期一致, 得到落叶松人工林的准确分布位置。通过谷歌地球确定分布记录的经纬度, 按照1个公里网格内只保留1条有效分布点的原则去除重复点, 共获得308条落叶松人工林的地理分布记录, 其中包括华北落叶松(Larix principis-rupprechtii.)132条, 日本落叶松(Larix kaempferi.)81条, 兴安落叶松(Larix gmelinii.)52条, 长白落叶松(Larix olgensis.)43条。

1.2 环境数据

通过国家青藏高原科学数据中心(http://data.tpdc.ac.cn/zh-hans/)下载的中国1km分辨率逐月平均气温、月最高温度、月最低温度及逐月降水量数据^[22—25], 据此计算出19个气候变量指标(表 1)驱动MaxEnt模型运行^[26]。通过地理空间数据云(http://www.gscloud.cn/home)下载数字高程数据(DEM)并利用ArcGIS 10.8中的空间分析工具计算坡度和坡向。采用1991—2020年的气候数据模拟当前气候条件下适生区的分布, 以30年为时间间隔计算1961—1990年和1931—1960年的气候数据, 数据的空间分辨率重采样为30″。分析时所用中国地图的底图来源于国家基础地理信息中心(http://www.ngcc.cn/ngcc/)。

1.3 模型建立

使用MaxEnt Version3.4.4k(https://www.gbif.org/tool/81279/maxent)进行适生区模拟分析。随机选取75%的落叶松人工林分布点作为训练集, 另外25%的分布数据作为测试集, 设置最大迭代次数为500次, 最大背景点数量设置为10000, 采用刀切法(Jackknife)评价环境变量对模型的重要性, replicate设置为10, 使模型重复运行10次, 结果取10次运行的平均值。模型输出的结果用ArcGIS 10.8进行分析。

1.4 模型精度检验

采用受试者工作特征曲线(ROC)的线下面积(AUC)来衡量模型的模拟精度。一般认为, AUC值大于0.9说明模型模拟效果“极准确”^[27—28]。在本研究中, 兴安落叶松人工林的AUC值为0.955, 长白落叶松人工林的AUC值为0.977, 华北落叶松人工林的AUC值为0.962, 日本落叶松人工林的AUC值为0.955, 落叶松人工林总体的AUC值为0.92(图 1), 说明模拟效果都达到了“极准确”的程度。

2 结果

利用MaxEnt模型将落叶松分布记录与环境变量数据进行处理分析, 得到不同类型落叶松的生态适宜度区划结果, 将结果导入ArcGIS, 根据适宜度重分类为4类：0—0.2、0.2—0.4、0.4—0.6和0.6—1, 分别对应着非适生区, 低适生区, 中适生区和高适生区^[29—30]。

2.1 影响落叶松人工林分布的主要因子

通过MaxEnt模型的Jackknife检验、模型贡献率和置换重要值来综合评估环境因子的重要程度, 根据环境因子变量响应曲线来判断落叶松人工林的存在概率与环境因子之间的关系, 当落叶松人工林存在概率大于0.5时, 认为所对应的环境因子的值有利于落叶松人工林的生长^[31—32]。

总体来看, 影响落叶松人工林分布的主要环境因子为最湿月降水量、海拔、最湿季度平均温度和温度季节变化方差, 落叶松人工林分布对主要环境因子的适宜范围如下, 最湿月降水量在100mm至200mm之间, 海拔介于200—2800m, 最湿季度平均温度在11℃至20℃之间, 温度的季节变化方差范围为50—70和100—250, 累计贡献率为73.8%(表 2, 图 2)。其中, 影响华北落叶松人工林分布的主要环境因子为海拔(1000—3000m, 括号内为适宜范围, 下同)、最湿月降水量(90—140mm)、温度季节变化方差(80—160)和最干季度降水量(20—40mm), 它们的累计贡献率为76.2%。影响日本落叶松的主要环境因子为最干季度降水量(40—120mm)、海拔(300—2800m)、温度季节变化方差(40—75)和最湿季度平均温度(13—23℃), 累计贡献率为59.6%。影响长白落叶松人工林分布的主要环境因子为温度季节变化方差(147—208)、年降水量(600—900mm)、最冷季度降水量(20—60mm)和坡度(1—6°), 累计贡献率为85.6%。影响兴安落叶松人工林分布的主要环境因子为温度季节变化方差、最湿月降水量、坡度和最湿季度平均温度, 温度季节变化方差介于160—300, 最湿月降水量介于110—175mm, 坡度介于1—7°, 最湿季度平均温度介于11—18℃(图 2, 表 2)。

2.2 落叶松人工林在1921—2020年的适生区

落叶松人工林在1991—2020年的适生区主要分布在我国的北方地区, 集中分布在长白山脉、大兴安岭、燕山、小陇山、太行山脉以及秦岭-大巴山等地区(图 3)。其中, 长白落叶松和兴安落叶松的分布比较相似, 都是集中在我国东北地区, 主要环境因子也都是温度的季节变化。相比于长白落叶松, 兴安落叶松人工林能在温度变化更大, 降水更少的高纬度生存, 主要集中在内蒙古东北部、黑龙江北部、吉林东部以及辽宁东北部地区, 潜在分布区面积为75.29×10⁴km², 其中高适生区面积为34.32×10⁴km²(表 3, 图 3)。长白落叶松集中分布在东北地区的长白山一带, 在大兴安岭区域也有部分适生区, 面积为39.52×10⁴km², 其中高适生区集中分布在长白山地区, 面积为19.05×10⁴km²(表 3, 图 3)。华北落叶松主要分布在我国的华北地区北纬33°到45°之间, 主要包括河北北部、甘肃东南部、陕西南部、内蒙古东南部以及山西大部分地区, 适生区的总面积约为68.43×10⁴km²。其中高适生区面积为18.08×10⁴km², 占适生区总面积的26.42%, 主要集中在燕山、太行山脉以及秦岭等地区(表 3, 图 3)。日本落叶松的分布区域比较偏南, 主要分布在北纬28°到37°之间以及辽东地区, 潜在适生区的面积为75.52×10⁴km², 其中高适生区面积为26.52×10⁴km², 集中分布在辽宁东部、吉林东南部和秦岭-大巴山地区, 其面积为潜在适生区面积35.12%(表 3, 图 3)。

2.3 落叶松人工林在历史气候下的适生区

落叶松人工林的适生区面积从1931—2020年不断增加, 其中1961—1990年间较1931—1960年间的高适生区增加了4.77×10⁴km², 当前气候条件下高适生区面积又进一步增加到53.57×10⁴km²(表 3)。

兴安落叶松人工林的高适生区在1931—1960年期间集中在大、小兴安岭及长白山地区, 适生区还包括河北北部、山西北部及辽宁西部地区, 适生区的面积为75.17×10⁴km²。到1961—1990年期间高适生区整体向北移动, 内蒙古东北部和黑龙江北部的大兴安岭地区由中低适生区转变为高适生区, 高适生区面积增加了8.81×10⁴km²(表 3, 图 4, 图 5)。与1961—1990年期间相比, 当前兴安落叶松人工林适生分布区的范围继续向北移动, 有部分适生区移出了我国北界, 高适生区的面积减少了0.93×10⁴km²(表 3, 图 3)。

1931—1960年长白落叶松适生区面积为20×10⁴km², 主要分布在长白山地区, 其中高适生区面积为5.63×10⁴km², 占28.15%。1961—1990年, 适生区面积减少了3.04×10⁴km², 辽宁西部地区的适生区大约向北移动了0.8°, 并且由零散的低适生区退化成了非适生区, 高适生区分布面积也显著减少, 退化为中适生区且呈破碎化分布(表 3, 图 4, 图 5)。与1961—1990年期间相比, 当前气候下长白落叶松人工林的高适生区由破碎化分布转变为聚集分布, 且面积持续增加, 增加的地区主要在吉林东部的长白山脉和黑龙江东南部, 但低适生区和中适生区的面积都有所减少。高适生区的纬度范围集中在北纬40°—48.5°之间。长白落叶松人高适生区的面积有所增加, 分布范围有向北移动的趋势。与1931—1960年期间相比, 当前气候下高适生区的范围向北移动了5°(图 3)。

与1931—1960年相比, 华北落叶松在1961—1990年期间适生区的面积增加了3.71×10⁴km², 面积增加的地区主要是河北北部和山西东部的太行山脉一带, 而高适生区范围变化不大(表 3, 图 4, 图 5)。与1961—1990年期间相比, 当前气候下适生区分布范围进一步扩大, 陕西北部和山西大部分地区由非适生区变为低适生区, 河北北部和甘肃东部的中低适生区转变为高适生区, 高适生区面积增加了10.61×10⁴km², 分布北界大约增加了3°(图 3), 这是适生区面积增加的主要原因。

日本落叶松在1931—1960年期间适生区面积为38.15×10⁴km², 高适生区面积为11.91×10⁴km², 从1931—1960年到1961—1990年, 日本落叶松适生区的面积增加了43.07%, 其中高适生区面积为7.56×10⁴km², 主要增加的区域是秦岭-大巴山一带, 而辽宁东部的适生区范围则有所减小(表 3, 图 4, 图 5)。较1961—1990年期间, 当前气候下适生区面积进一步地增加了20.94×10⁴km², 这些区域主要是甘肃小陇山林区、山东的泰山和沂蒙山地区以及吉林东南部和辽宁东部山地。高适生区在辽东山区呈先减少后增加的趋势, 在甘肃东南部小陇山林区显著增加, 而在四川、重庆、湖北地区则是先增加后减少的趋势。位于辽宁东部的高适生区向北移动了2°, 位于陕西南部、湖北西部、四川北部和重庆地区的高适生区则没有明显的偏移, 但是高适生区的分布范围却向高海拔地区集中(图 3)。

3 讨论

我国近百年来温度持续增加, 东北地区的升温趋势十分明显, 尤其是20世纪80年代以后升温最为剧烈。降水量的变化虽未表现出明显的趋势性, 但具有一定的阶段性, 1936—1959年降水较为丰富, 1974—1983年降水则开始减少。总体来上随着气温的增加, 东北地区出现暖干化的趋势^[33]。唐国利等研究了中国近百年的年均温变化趋势, 发现年平均温度的上升幅度为(0.78±0.27)℃^[34]。生长季热量不足是限制落叶松向北扩展的关键因子, 因此温度的升高有助于落叶松的适生范围向北扩张, 适生面积也有了大幅度提升。而兴安落叶松的适生范围基本不变, 原因是随着适生区向北移动, 部分兴安落叶松的适生范围移出我国北界。有类似研究表明, 气候变化会导致北半球植被向高纬度地区迁移^[35]。李峰等和冷文芳等使用未来情景的气象数据同样发现落叶松适生区有向高纬移动的趋势^[36—37], 这说明落叶松适生区这种向北扩展的趋势是伴随着气候变化而长期持续的。

长白落叶松耐寒、喜生长在湿凉的环境中^[38—39], 本研究发现长白落叶松适生区主要集中在长白山地区, 制约长白落叶松分布的主要气象因子为温度季节变化方差、年降水量和最冷季度降水量(图 6, 图 7), 近些年来高适生区呈现先减少后增加的趋势, 长白山地区这一时期影响长白落叶松分布的年降水量也是呈现先减少后增加的趋势, 高适生区减少的主要原因是在1961—1990年期间长白山东北部地区的年降水量减少至600mm以下, 低于长白落叶松对年降水量需求的阈值, 因此在本时期长白山北部地区的高适生区面积减少。当前长白落叶松适生区增加的主要原因是长白山地区温度和降水的增加, 温度季节变化方差的减小使得适生区北部边缘地区变得适宜生存, 长白山地区年降水量在600—900mm之间的区域也在显著增加, 使得适生区北界向北扩张, 因此长白落叶松适生区的向北移动的现象是与气候变暖的趋势相适应的。兴安落叶松人工林的高适生区主要分布在我国的大、小兴安岭及部分长白山地区^[40], 影响兴安落叶松分布的主要气候因子为温度季节变化方差、最湿月降水量和最湿季度平均温度(图 6, 图 7), 过去90年温度季节变化方差持续减小, 且适宜其生存的温度季节变化方差的范围向高纬度移动, 导致在大兴安岭地区的高适生区面积显著增加, 张广才岭地区高适生区面积减小。虽然兴安落叶松适生区范围在近些年来移动幅度不大, 但是适生区域内的适宜性程度发生了显著的变化, 大兴安岭北部最湿月降水量增加和最湿季度平均温度的增加同样导致兴安落叶松在该区域高适生区增加。随着东北地区暖干化的持续, 兴安落叶松的高适生区将持续北移^[41], 导致有部分移出我国北界, 这是我国兴安落叶松适生区减少的主要原因。冷文芳等模拟预测兴安落叶松适生区与气候的关系时发现, 气温每上升1℃, 兴安落叶松在我国的适生区面积将减少12%, 温度对兴安落叶松分布的影响远大于降水^{[19, 20, 37]}。

最干季度降水量是制约日本落叶松人工林分布的关键要素, 但影响适生区分布的其他气象因子不集中, 导致它们的个体贡献率较低。原因是日本落叶松适应性强, 在我国北方、中部、西部地区都有引种种植^[42], 导致适生区的南北跨度较大, 影响因素复杂。在辽东山区高适生区呈先减少后增加的趋势, 与最湿季度平均温度和最干季度降水量变化一致, 导致1961—1990年期间辽东地区适生区减少的原因主要是最湿季度平均温度的降低和最干季度降水量的减少。1961—1990年期间日本落叶松增加的主要原因是辽东地区北部区域的最干季度降水量显著增加, 达到40mm—120mm的范围, 满足日本落叶松的生存需求, 使日本落叶松在辽东地区的适生区向东北方向迁移, 面积不断增加, 近30年来随着气候变暖的加剧, 辽东地区最湿季度平均温度表现出增加的趋势, 这同样能解释适生区面积的增加。除了温度和降水外, 地形因素也是影响日本落叶松地理分布格局的关键因素, 位于秦岭-大巴山地区的高适生区向高海拔区域收缩, 这是由于该地区多为高海拔山地, 地形复杂, 导致该地区的日本落叶松无法向高纬度迁移, 因此向高海拔区域收缩。

刘宪钊等研究发现, 只考虑地形因子和气候因子时, 华北落叶松主要影响因子为海拔、坡度和最湿润月的降水量^[43]。本研究同样发现海拔是影响华北落叶松分布的最重要因子, 表现为华北落叶松主要分布在海拔相对较高的山区。影响华北落叶松分布最重要的气候因子为最湿月降水量(图 6、图 7), 在暖干的气候趋势下, 华北地区整体表现出温度升高和降水减少的趋势, 但是在高海拔地区随着温度不断增加导致降水的季节差异较大, 最干季度降水量表现出减少的趋势, 最湿月降水量却呈现出增加的趋势, 导致适宜华北落叶松生存的区域增加, 使得华北落叶松在河北坝上和山西太行山地区等高海拔区域高适生区的面积不断增加。由于海拔能够通过影响水热条件间接影响华北落叶松的生存, 其他气象因子所占比重较低, 因此气候变化对华北落叶松的影响是相对较小的。在气候变化的背景下, 四种不同落叶松人工林的适生区变化规律有所不同, 适生区面积变化也不一致。根据研究结果推测, 华北落叶松、日本落叶松和长白落叶松可以继续扩大造林面积, 但兴安落叶松应当减少在大兴安岭东麓及长白山地区的种植。

MaxEnt模型假定物种与环境之间达到一种平衡的状态^[44], 即此时的环境指标被设定为最适合该物种生存的状态, 但这并不意味着物种与环境达到了事实上的最佳状态。而且植物本身由于迁移能力较差, 导致有些植物不能到达适合的生境进而对模拟结果造成了一定影响。MaxEnt模型使用19个气候变量做计算时, 对高海拔地区的模拟结果可能会有一些偏差^[45]。尽管还存在上述不足, MaxEnt依然是模拟物种潜在适生区最常用且可靠的方法。

4 结论

本研究使用过去90年的气象数据结合地形数据和落叶松人工林的地理分布数据, 对中国四种主要落叶松人工林适生区的空间分布进行了模拟和空间变化分析。主要结论如下：

(1) 总体来看, 影响落叶松人工林分布的主要环境因子由大到小依次为：最湿月降水量、海拔、最湿季度平均温度和温度的季节变化方差, 但不同落叶松树种间的主要影响因子有较大的差异, 海拔和最干季度降水量分别是华北落叶松和日本落叶松最主要的影响因子, 而温度季节变化方差是长白落叶松和兴安落叶松最主要的影响因子。

(2) 落叶松人工林的适生区主要分布在我国北方地区。其中, 华北落叶松主要集中在燕山、太行山脉以及秦岭等地区; 日本落叶松主要集中分布在辽东地区和秦岭-大巴山地区; 长白落叶松集中分布在长白山脉; 而兴安落叶松集中在大兴安岭及长白山部分地区。

(3) 从90年的变化趋势来看, 落叶松适生区的范围整体上呈现出向北移动的趋势, 同时面积也在增加。其中华北落叶松和日本落叶松的适生区面积一直呈增加的趋势。长白落叶松适生区呈现先减少后增加的趋势。兴安落叶松在近些年有适生区减少的趋势。