1880 —2012年以来，全球地表平均温度升高了0.85 (± 0.2)℃，预计到21世纪末仍将上升0.3~4.8 ℃，(IPCC，2014)。气候变暖作为影响植被地理分布的主要因子之一(Davis et al., 2001；Pearson et al., 2003)，对植被分布格局产生极大的影响(Esteve-Selma et al., 2012；Guiterrez et al., 2014)。探究树种分布与气候因子之间的关系并预测当前及未来不同气候情景的树种分布格局变化已成为生物地理学和生态学研究的热点(Wang et al., 2016；Chakraborty et al., 2016；Dyderski et al., 2018)。

大尺度树种分布格局是多个气候因子共同作用的结果。韩国北汉山地区蒙古栎(Quercus mongolica)分布与年均气温、温暖指数和潜在蒸散量高度相关(Yang，2013)，而中国东北地区蒙古栎分布则主要受季节性降水量和年均气温的影响(贾翔等，2017)；在欧洲中部喀尔巴阡山脉西北部地区，影响栎属(Quercus)植被分布的主要限制因子是5—7月份降水，而东北部地区则是冬季和早春的气温(Nechita et al., 2017)；Li等(2016a)研究表明最湿月份降水量、大于0 ℃年积温和最冷月最低气温是影响油松(Pinus tabulaeformis)分布的主要气候因子，而在更大的地理尺度(Wang et al., 2016)年均气温和年均降水量则是限制油松分布的主要因子；云南地区的针叶树种分布主要受冬季低温影响(Li et al., 2016b)，而东北地区则与年均降水量和生物寒冷指数相关(Leng et al., 2008a)。由此可见，不同地区影响树种分布的主导气候因子不同，阐明特定地区、不同优势树种适宜分布格局及其与气候因子的关系对制定有效的经营策略十分必要。

物种分布模型是评估气候变化对物种分布影响的主要方法(王娟等，2006；Guisan et al., 2000)，其主要利用环境数据和物种的分布数据, 依据特定的算法估计物种的生态位，以概率的形式反映物种对生境的偏好程度(李国庆等，2013)。随机森林模型(Wang et al., 2016; Leng et al., 2008a)、生物气候模型(Ge et al., 2017)、广义线性模型(Guisan et al., 2000)和最大熵(MaxEnt)模型(Zhang et al., 2013；Yang et al., 2013)等是目前应用广泛的物种分布模型。与其他模型相比，MaxEnt模型具有样本量要求低(陈新美等，2012)、变量处理灵活(只需要发生数据)、降噪效果好、预算精度高等优点(许仲林等，2015)。因此，本研究采用MaxEnt模型进行研究区主要优势树种适宜分布区预测。

华北落叶松(Larix principis-rupprechtii)作为速生针叶树种，是华北亚高山地区的主要造林树种(赵海燕等，2015)，在河北省主要分布在承德、张家口和保定等地的山区；油松为中国特有树种，耐低温、干旱和瘠薄，是北方地区最主要的造林树种之一；蒙古栎喜光、耐寒冷，是我国温带地区阔叶林或针阔混交林的主要树种之一(任宪威等，1997)。以上3个树种在河北地区分布最广泛，也是当地最主要的造林和用材树种。本研究以华北落叶松、油松和蒙古栎为研究对象，探讨当前3个树种的适宜分布格局及其主导气候因子，预测未来气候变化情景下这3个树种的适宜分布区变化，以期为评估气候变化对森林的影响及制定适宜森林经营策略提供理论依据。

河北省(113°27′—119°50′E，36°05′—42°40′N)横跨华北、东北两大地区，环抱北京。地势西北高，东南低，由西北向东南倾斜。地貌复杂多样，有坝上高原、燕山山地、太行山山地和河北平原等地貌单元。气候属于温带大陆性季风气候，年日照时数2 303 h，全年无霜期81~204天，年均降水量300~800 mm，1月平均气温在3 ℃以下，7月平均气温18~27 ℃。森林类型主要属于暖温带落叶阔叶林，亦有少量的寒温带针叶林。主要树种有蒙古栎、油松、侧柏(Platycladus orientalis)、华北落叶松和云杉(Picea asperata)等(任宪威等, 1997)。

基于中国自然标本馆(http://www.cfh.ac.cn)、生物多样性数据库(https://www.gbif.org)、中国数字植物标本馆(http://www.cvh.ac.cn)、中国科学院植物研究所植物标本馆(http://pe.ibcas.ac.cn)和国家标本平台(http://www.nsii.org.cn)记录的标本数据，以及河北省森林资源第四次二类调查数据和第八次一类清查数据，提取出华北落叶松、油松和蒙古栎天然分布点数据。同时，为避免分布点过度聚集造成权重过高，为每个发生点构建5 km的缓冲区，去除缓冲区中叠加的样地(戚鹏程, 2009)，最终确定华北落叶松、油松和蒙古栎发生点分别为103，287和425个(图 1)。

图 1 河北省3个树种分布点 Fig. 1 Distribution points of three tree species in Hebei Province

采用高分辨率气候模型ClimateAP生成研究区气候数据。ClimateAP模型是加拿大不列颠哥伦比亚大学Wang等(2012b)在其研发的ClimateWNA模型基础上建立的适用于亚太地区的气候模型。该模型采取双线性插值和动态局部回归的方法，将气候数据降尺度为无尺度的点数据，通过经纬度坐标和海拔生成亚太地区任意点的无尺度月、季、年和未来某时间节点的高精度气候数据集。植物的耐寒性、完成生活史所需的热量供应和可利用的水资源是限制植物地理分布的主要因子(Woodward, 1987)。基于此，本研究利用物种发生点的经纬度和海拔信息以及ClimateAP模型生成了10个气候因子数据：年均气温、最热月平均气温、最冷月平均气温、大于5 ℃年积温、年均降水量、湿季降水量、干季降水量、哈格里夫斯参考蒸散量和小于0 ℃年积温，分辨率1 km × 1 km，投影坐标系Xian_1980_GK_Zone_20。

IPCC第五次报告中提出4种典型浓度路径情景：RCP2.6、RCP4.5、RCP6.0和RCP8.5，用以对未来气候变化进行模拟(IPCC，2014)。其中RCP2.6代表严格控制温室气体排放气候变化情景，RCP4.5和RCP6.0代表中度控制温室气体排放气候变化情景，RCP8.5代表温室气体高度排放气候变化情景，且RCP4.5的优先性大于RCP6.0(Moss et al., 2010；应凌霄等，2016)。本研究最终选取RCP2.6、RCP4.5和RCP8.5这3种情景，分别代表低、中、高3个温室气体排放情景模拟未来气候变化。

采用MaxEnt模型对华北落叶松、油松和蒙古栎当前及未来气候变化情景下的适宜分布区进行预测。为评估模型精度，随机选取样本点的75%作为训练集，25%作为测试集。采用受试者工作特征曲线(ROC)分析法进行模型精度验证，其曲线下面积(AUC)值的大小可以很好地反映模型预测效果，取值0~1。AUC值越大表明预测精度越高(王运生等，2007)：0.5~0.6为较差，0.6~0.7为一般，0.7~0.8为较准确，0.8~0.9为很准确，0.9~1为极准确(Vanagas，2004；Swets，1988)。同时为避免数据过度拟合，将正则化乘数值设为0.1(Li et al., 2016b)。在建立模型过程中，采取刀切法(Jackknife)和响应曲线分析各气候因子的贡献率及其对各树种分布的影响。模型结果以Logistic格式输出为ASCII文件，利用ArcGIS10.2转换为栅格数据，取值0~1。采用“10%训练存在逻辑阈值”(Zhang et al., 2013)作为适宜分布阈值界点，低于此阈值划分为不适宜区，高于此阈值划分为适宜区。

未来气候变化的不确定性是评估气候变化对森林生态系统和树种分布影响的最大挑战(Wang et al., 2012a)。为尽可能减小这种不确定性对评估结果的影响，本研究基于3个气候变化情景对未来2个时期(2040—2069年和2070—2099年)3个树种的适宜区进行预测，将3个气候变化情景预测树种分布适宜区重叠部分作为该树种的分布适宜区，以提高预测精度。

3个树种分布模型的ROC分析表明(图 2)，模型预测精度优于随机预测。华北落叶松、油松训练数据和测试数据的AUC值都大于0.9，达到极准确水平；蒙古栎训练数据和测试数据的AUC值都大于0.85，达到很准确水平。表明MaxEnt模型对3个树种均具有较好的预测能力，可用于华北落叶松、油松和蒙古栎的适宜分布区预测。

图 2 MaxEnt模型精度检验的ROC曲线 Fig. 2 Receiver operating characteristics curves of test for the MaxEnt models

3个树种当前主要分布在河北省燕山和太行山沿线地区(图 3)。华北落叶松适宜分布阈值为0.2，分布区位于海拔较高地带(＞1 000 m)，主要分布在燕山地区西北部和坝上高原东南部，太行山北段的高海拔地区亦有少量分布，总面积相对较小，仅占整个研究区的13.46%；蒙古栎的适宜分布阈值为0.32，主要分布于燕山和太行山中、高海拔地带。河流冲刷形成的河漫滩平原和低海拔地区不适宜蒙古栎生长，分布区较破碎，但分布面积相对较大，占研究区总面积的20.54%；油松适宜分布阈值为0.35，适宜分布区和蒙古栎的大体一致，但油松的范围更广一些，分布的海拔下限明显低于蒙古栎，整个燕山和太行山地区都适宜油松生长，适宜分布面积最大，占研究区总面积的23.84%。

图 3 华北落叶松、蒙古栎、油松当前的适宜分布区 Fig. 3 Current suitable distribution area of L. principis-rupprechtii, Q. mongolica, and P. tabulaeformis

刀切法分析表明(表 1)，影响华北落叶松分布的主要气候因子包括小于0 ℃积温和湿季降水量，这2个气候因子累计贡献率为80.5%；影响蒙古栎分布的主要气候因子包括最热月均气温、哈格里夫斯水分亏缺和湿季降水量，这3个气候因子累计贡献率为89.4%；影响油松分布的主要气候因子包括最热月平均气温、湿季降水量、大于5 ℃年积温和年均气温，这4个主要气候因子累计贡献率为83.5%。以上各树种主要气候因子的累计贡献率均超过80%，能够较好的解释3个树种的分布格局。

表 1 气候因子贡献率 Tab.1 Contribution rate of climatic factors

表 1 气候因子贡献率 Tab.1 Contribution rate of climatic factors % 气候因子 Climatic factor 华北落叶松 L. principis-rupprechtii 蒙古栎 Q. mongolica 油松 P. tabulaeformis 年均气温Mean annual temperature 4.5 0.6 6.3 最热月平均气温Mean warmest month temperature 2.8 53.5 43.6 最冷月平均气温Mean coldest month temperature 1.2 2.0 0.5 大于5 ℃年积温Annual temperature above 5 ℃ 2.6 0.4 13.7 年均降水量Mean annual precipitation 1.1 0.3 0.9 湿季降水量Mean precipitation of wet quarter 16.1 16.5 19.9 干季降水量Mean precipitation of dry quarter 2.8 2.8 4.5 哈格里夫斯参考蒸发量Hargreaves reference evaporation 2.7 1.5 4.7 哈格里夫斯水分亏缺Hargreaves climatic moisture deficit 1.8 19.4 4.8 小于0 ℃积温Annual temperature below 0 ℃ 64.4 3.0 1.1

响应曲线分析表明，3个树种均有适宜其生长的气候条件(图 4)。华北落叶松发生概率随大于0 ℃年积温(1 200~1 950 ℃)的增加先增加到极值后骤降，而随着湿季降水量(＞350 mm)的增加，则增加到一定程度后趋于平稳。蒙古栎发生概率随哈格里夫斯水分亏缺(150~360 mm)和最热月平均气温(19.3~24.8 ℃)的增加先增加后下降，而随湿季降水量(＞385 mm)的增加持续增加后趋于平稳。与之相似，油松发生概率随湿季降水量(455~730 mm)、最热月平均气温(22~25.5 ℃)和年均气温(4~11.5 ℃)的增加，先增加到一定程度后迅速下降，而随着大于5 ℃年积温(＞2 300 ℃)的增加，持续增加后趋于平稳。

图 4 华北落叶松、蒙古栎、油松发生概率与主要气候因子的关系 Fig. 4 Relationships between major climatic variables and the presence probability of the L. principis-rupprechtii, Q. mongolica, and P. tabulaeformis

本研究基于3个气候变化情景进行适宜分布区一致性预测。2040—2069年，3个树种适宜分布区预测具有较高的一致性(图 5)，表明树种地理分布变化的不确定性在21世纪中期相对较低；然而，预测的一致性在2070—2099年降低(图 6)，表明预测时间序列越长，预测结果的不确定性越高。

图 5 2040—2069年华北落叶松、蒙古栎、油松适宜分布区 Fig. 5 Suitable distribution area of L. principis-rupprechtii, Q. mongolica, and P. tabulaeformis in 2040—2069 period

由表 3可知，在气候变化情景下，未来时期3个树种适宜分布范围在经纬度上均无明显的迁移(经纬度变化＜1°)，但都有向高海拔地区扩散的趋势，以蒙古栎最为显著，在2040—2069年海拔上升至1 233 m，在2070—2099年上升至1 366 m。2040— 2069年，华北落叶松适宜分布范围向西大幅扩散至整个坝上高原地区，面积显著增大，由占总面积的13.46%上升至15.03%；蒙古栎分布面积变化很小(面积缩小幅度＜1%)；油松分布面积显著缩小，由23.84%降为18.82%。在2070—2099年，3个树种的适宜分布面积都显著缩小，且缩小幅度都大于3%(图 5)。

表 2 华北落叶松、蒙古栎、油松当前和未来适宜分布区变化^① Tab.2 Current and future appropriate distribution changes of L. principis-rupprechtii, Q. mongolica, and P. tabulaeformis

表 2 华北落叶松、蒙古栎、油松当前和未来适宜分布区变化① Tab.2 Current and future appropriate distribution changes of L. principis-rupprechtii, Q. mongolica, and P. tabulaeformis 变化因子 Variables factor 华北落叶松 L. principis-rupprechtii 蒙古栎 Q. mongolica 油松 P. tabulaeformis 当前Current 2040—2069 2070—2099 当前Current 2040—2069 2070—2099 当前Current 2040—2069 2070—2099 经度Longitude(°) 117.0 116.4 116.6 117.1 116.5 116.3 117.4 116.8 116.8 纬度Latitude(°) 41.5 41.6 41.7 40.8 41.2 41.3 40.5 40.8 40.9 海拔Altitude/m 1248 1367 1394 954 1233 1366 707 973 1025 面积Area(%) 13.46 15.03 11.61 20.54 19.98 14.12 23.84 18.82 14.30 ①经纬度和海拔为区域内平均值Latitude longitude, and altitude are the average values in the study area.

研究表明，华北落叶松、蒙古栎和油松均有适宜生长的气候条件和分布区，且适宜分布区随气候变化而改变。气候变化的不确定性是预测树种适宜分布区变化的最大挑战(Wang et al., 2012a)，而温室气体排放水平的不确定性是造成气候变化不确定性的重要原因(IPCC，2014)。据IPCC报告(2014)可知，RCP2.6、RCP4.5和RCP8.5情景下地表温度预计到21世纪中期可能比当前升高0.3~0.7 ℃，差异较小；而到21世纪末期，3种情景下可能分别比当前升高0.3~1.7，1.1~2.6和2.6~4.8 ℃，差异变大。未来树种适宜分布区变化的一致性程度随气候变化特征差异的变大而减小。不同气候变化情景对树种预测适宜分布区的影响在2040—2069年较小(图 5)，但在21世纪末期较大(图 6)，故本研究只对21世纪中期的树种适宜分布区变化和相应机制进行探讨。

图 6 2070— 2099年华北落叶松、蒙古栎和油松适宜分布区 Fig. 6 Suitable distribution area of L. principis-rupprechtii, Q. mongolica, and P. tabulaeformis in 2070-2099 period

在2040—2069年，华北落叶松适宜分布范围向西大幅扩散至整个坝上高原地区，面积显著增大(图 5)，可能由于华北落叶松有一定的耐湿、耐旱和耐瘠薄能力，且对不良气候影响抵抗力强。Leng等(2008b)研究表明，在东北地区气候变暖条件下华北落叶松将会向北扩散，适宜分布面积显著增大，与本研究结果的扩散方向有差异。本研究中华北落叶松向西扩散可能受河北省小于0 ℃年积温变化影响，随着气候变暖，河北省小于0 ℃年积温整体变小，华北落叶松小于0 ℃年积温的适宜区向西移动，造成华北落叶松向西迁移。未来时期，河北地区发生干旱的频率加大(Li et al., 2017)，加上华北落叶松适宜分布区变化和其较耐干旱及对土壤适应性广的生物学特性，在坝上高原地区造林中适当扩大华北落叶松的营造面积，并进行多功能高效培育，是在气候变化背景下优化树种及精准提高森林质量和生产力的关键。

贾翔(2017)等研究表明，在东北地区以蒙古栎等树种为主的阔叶红松林在气候变暖条件下将会向大兴安岭等海拔较高地区扩散且适宜分布区面积变化很小，本研究结果与其相似，这可能与蒙古栎对气候、土壤适应性广，但不耐水湿有关。气候变暖使河北省温度升高，气温带沿海拔梯度移动，降水量增加，气温和降水的变化影响互相抵消，从而造成蒙古栎适宜分布区在海拔梯度变化显著，面积变化很小。扩散后蒙古栎的适宜分布区仍较破碎也很可能与此有关。燕山地区河流多与山脉直交，切穿山地形成南北向河谷，河谷地区比较湿润，造成这些地区不适宜蒙古栎生长，蒙古栎适宜分布区破碎。当前及未来时期，蒙古栎适宜分布范围和油松重叠区域显著缩小，而与华北落叶松重叠区域则显著增大(图 5)。在气候变化情景下，通过人工调控措施营造华北落叶松与蒙古栎混交林，开展不同混交比例和混交模式配置试验，有助于培育气候适生区健康、稳定森林生态系统。

当前，油松分布格局和蒙古栎相似，但在2040—2069年，油松适宜分布范围大幅度缩小(收缩幅度＞5%)，远远超过蒙古栎的收缩幅度(＜1%)，很可能由于油松喜干冷的气候，而未来气候变化以暖湿为主。油松适宜分布区沿东南边缘收缩，但并没有像蒙古栎一样向西北扩散(图 5)。影响油松和蒙古栎分布的最主要气候因子均是最热月平均气温，且两者对最热月平均气温的响应曲线相似(图 4)。与蒙古栎相似，随气候变暖，油松适宜分布区海拔明显升高，受河北省地势影响，沿东南边缘收缩；但油松的湿季降水量的适生范围明显小于蒙古栎且其对较小的大于5 ℃年积温耐受性很低(图 4)，导致油松在较高海拔地区不适宜生存。Wang等(2016)研究表明，油松适宜分布区将会大幅度扩张，本研究结果与其相反，这可能是研究区尺度差异造成的。Wang等(2016)对油松适宜分布区的模拟建立在整个中国区域，本研究仅对河北省境内的油松适宜分布区进行模拟。在气候变化背景下，河北省油松适宜分布范围大幅度缩小应当重视，因为无论从经济价值还是生态功能来看，油松都是河北省最重要的树种之一。本研究结果应为决策者提出预警，要适度调整研究区的油松经营策略，以应对气候变化的影响。

本研究采取3个不同的气候变化情景进行河北省3个树种的适宜分布区一致性预测，这在一定程度上降低了未来气候变化的不确定性对预测结果的影响。预测结果可为气候变化下优化河北省主要森林资源经营策略提供必要的科学依据。本研究仅基于气候因子进行了3个主要树种的适宜分布预测，没有考虑土壤(Li et al., 2012)、病虫害(Poland et al., 2006)、地形和人为干扰(Magnani et al., 2007)等因素的影响。在气候变化背景下，综合考虑多因素影响是今后进一步努力研究的方向。

2040—2069年，随气候变化，华北落叶松、油松和蒙古栎的适宜分布区可能发生改变，都有向海拔较高地区扩散的趋势，其中华北落叶松适宜分布区面积显著增加，油松显著缩小，蒙古栎变化较小。华北落叶松和蒙古栎有向西北方向扩散的趋势，两者适宜分布区重叠面积显著增加，而蒙古栎和油松适宜分布区重叠面积显著缩小。未来张家口和承德地区的造林和生态修复工作应考虑到这种变化趋势，有利于适应气候变化，提高森林生产力，构建健康稳定的森林生态系统。