气候变化下三明市主要针叶树种潜在地理分布的预测

张今朝, 刘 健, 余坤勇, 曾 琪, 郑文英, 艾婧文, 上官莎逸

(1.福建农林大学林学院；2.福建农林大学园林学院；3.福建省资源环境监测与可持续经营利用重点实验室，福建 福州 350002)

气候是决定地球上物种分布的最主要因素，而物种分布格局的变化最能直接和间接地反映气候变化[1-3].IPCC研究表明：近百年来(1880—2012年)全球平均增温0.85 ℃[4]，预计到21世纪末地表平均温度将上升2.6～4.8 ℃[5]，全球变暖已是无需争辩的事实.随着全球气候变化，生物的适宜分布区也将发生改变，预测气候变化对物种潜在分布格局的影响已成为生物地理学和全球变化研究的热点和重点领域[6-7].目前，用于物种分布预测的生态位模型主要有Garp模型、Bioclim模型、Domain模型和MaxEnt模型等.其中，MaxEnt模型以运行速度快、运行时间短、易于操作、运行结果稳定[8-9]等优势，得到广泛应用.在国内，MaxEnt模型被广泛地应用于气候变化对物种分布影响，包括粮食作物适生区的预测[10]、入侵物种的分布预测[11]、濒危物种保护[12]、药用植物潜在地理分布模拟[13]以及病虫害分布[14]等方面，并且均取得较好的模拟效果.

杉木(Cunninghamialanceolata(Lamb.) Hook.)和马尾松(PinusmassonianaLamb.)是南方主要的针叶树种，也是中国南部的主要用材树种，经济价值高，用途广.三明市是我国南方重要集体林区之一，也是我国林改唯一试点和海峡两岸现代林业合作实验区，境内森林覆盖率高，森林资源丰富，杉木、马尾松种植历史悠久.但一直以来，杉木和马尾松种植区域的选择主要依靠林农的传统种植经验，缺少合理规划，在一定程度上迟缓了树木的生长速度，降低了林地资源的利用率.因此，如何在当前及未来气候变化、经济发展背景下，减少因盲目种植带来的树种和林地资源的浪费、确保三明主要针叶树种的科学合理种植成为目前待解决的问题.为此，本研究以三明市三元区、梅列区的杉木和马尾松为研究对象，基于地理分布信息和气候环境数据，利用MaxEnt模型，模拟预测杉木、马尾松在三明市当前(1950—2000年)、未来气候情景下(2050年和2070年)潜在适宜生境分布范围和空间格局的变化规律，并依据野外调查结果验证模型预测的准确性.本研究的结果将对三明市主要针叶树种分布区的优化调整、林地资源的合理利用提供辅助资料，对应对气候变化、促进林业可持续经营具有重要的理论和现实意义.

1 研究区域概况

研究区位于福建省中部偏西隶属于三明市三元区、梅列区(26°01′N—26°25′N， 117°16′E—117°48′E)，下辖岩前镇、陈大镇、莘口镇、城东乡、徐碧乡、洋溪乡以及中村乡7个乡镇，西北与明溪县衔接，东北与沙县相邻，东南与大田县接壤，西南与永安市毗连.研究区地理位置如图1所示，地处沿海内陆山区，海拔多在200～800 m，地形以丘陵和低山为主，境内山岭起伏，连绵不绝，低山、丘陵、沟谷、盆地犬齿交错，溪流纵横，形成“风小、湿度大”的地理环境.气候温暖湿润，四季分明，雨量充沛.年平均气温14～19 ℃，年平均雨量1 500～1 900 mm，其中1—3月是低温多雨的冷湿期，5—6月为梅雨季节，7—9月为高温少雨的伏旱期.土壤质地多为沙土壤或者轻土壤，腐殖质多为灰黑色，厚度4.5～15.5 cm，土层深厚肥沃.

图1 研究区地理位置Fig.1 Geographic location of the study area

2 研究方法

2.1 数据收集

杉木和马尾松在三明市的分布记录共通过3种途径获取，包括中国数字植物标本馆(http://www.cvh.ac.cn/)、国内公开发表的相关论文和森林二类调查小班数据等.其中杉木和马尾松的分布点分别选取50、100、150、200、250和300个.

环境数据来源于世界气候数据库(WORLDCLIM，http:/ /www.worldclim.org)，该数据库采用插值法将全球气象站点数据(1950—2000年)生成世界气候栅格数据，并衍生出对物种分布有重要影响的19个生物气候变量(Bioclim variables)[7,15]，变量名称为Bio1-Bio19，包括年均温(Bio1)、昼夜温差月均值(Bio2)、等温性(Bio3)、温度季节性变化标准差(Bio4)、最暖月最高温(Bio5)、最冷月最低温(Bio6)、温度年较差(Bio7)、最湿季均温(Bio8)、最干季均温(Bio9)、最暖季均温(Bio10)、最冷季均温(Bio11)、年降水量(Bio12)、最湿月降水量(Bio13)、最干月降水量(Bio14)、降水量季节性变异系数(Bio15)、最湿季降水量(Bio16)、最干季降水量(Bio17)、最暖季降水量(Bio18)和最冷季降水量(Bio19).研究采用该数据库中的 3 期气候数据，分别为当前气候(1950—2000年)及未来气候情景(2040—2060年、2060—2080年)，空间分辨率为30″×30″(约0.85 km×0.85 km).未来气候选用了IPCC第五次评估报告发布的BCC.CSM1-1模式[16]代表温室气体排放浓度中等的rcp45场景，其对温度和降水各变量的模拟和预测都有较高的准确性[17-18].研究利用GIS软件将下载的19个生物气候变量进行裁剪，得到各气候图层数据.

2.2 模型分析

基于杉木和马尾松的分布数据和环境因子采用MaxEnt(v3.3.3)分析其分布与环境因子的关系，并预测BCC-CSM1-1rcp45场景下杉木和马尾松21世纪50年代和70年代的分布状况.首先，随机选取75%的分布点作为训练集(training data)，剩余的25%作为测试集(testing data)[19].采用接受者操作特性曲线(receiver operating characteristic curve, ROC)分析法进行模型精度检验，当前，ROC曲线分析法在物种潜在分布预测模型评价中的应用十分广泛[20-21].ROC曲线基于非阈值依赖判断(threshold-independent evaluation)模型精度，即以预测结果的每一个值作为可能的判断阈值，由此计算得到相应的灵敏度和特异度[22].以特异度(1-specificity)即假阳性率为横坐标，以灵敏度(1-omission rate)即真阳性率为纵坐标绘制成ROC曲线，AUC值是ROC曲线与横坐标围成的面积值(the area under the ROC curve)，其大小能很好地说明模型模拟值的准确性，取值范围为[0,1]，此值越大表示模型判断力越强，模型预测结果越准确[23].AUC评估标准为：0.50～0.60(预测失败)； 0.60～0.70(较差)；0.70～0.80(一般)；0.80～0.90(好)；0.90～1.0(非常好)[24].同时运用刀切法(Jackknife)确定各环境因子的权重，筛选出影响研究区杉木和马尾松分布的主导环境因子[25]，刀切法常被用来分析各环境因子对预测结果的影响程度，或用来确定影响物种分布的关键性环境因素[26-27].最后，利用GIS软件将生成的ASCII格式数据转换为栅格数据，并进行重分类，依据适宜生境评价指数[1,7,24,28]，将研究区生境分为3类：即不适宜生境(0～0.25)、低适宜生境(0.25～0.50)、适宜生境(0.5～1.0).同时，统计未来气候情境下各适宜生境的面积及占研究区面积的比例，得出气候变化条件下杉木和马尾松的适宜生境分布格局.

3 结果与分析

3.1 样本量对模型预测精度的影响

有研究表明[29-30]：物种的分布数据量会影响模型模拟的准确度.一般情况下，样本量增加，模型模拟精确度也会增加，当样本量达到一定临界值时，模型精度的增加幅度会慢慢减少直至不再增加，此时模型将达到最大准确度[31].因此，研究分别选取样本量50、100、150、200、250和300个进行模拟预测，探讨样本量的大小对模型预测精度的影响，结果如表1所示.

由表1知，随着样本量的变化，AUC值出现先增加后降低的趋势，但取值变化不大，说明本研究中样本量的多少对模型模拟的精度影响较小.为使模型模拟达到最大准确度，研究选取模拟精度最高的样本量，即杉木选取200个，马尾松选取150个来进行树种的适生区预测.而对于使用不同模型进行模拟验证，需要在今后的研究中进行探讨，以确保和提高预测的准确性.

表1 不同样本模拟的AUC预测精度Table 1 AUC accuracy of different sample simulations

3.2 杉木和马尾松潜在分布区的预测结果

表2 模型AUC预测精度Table 2 Prediction precision of AUC

研究基于当前、2050年、2070年的气候数据预测三明市主要针叶树种的潜在分布及变化趋势，利用ROC曲线作为模型精度的评价指标，结果训练集和测试集的AUC值均大于0.9，如表2所示，表明构建模型的预测精度高，可以用于分析杉木、马尾松地理分布与气候变化的关系.

图2为气候变化下杉木和马尾松的潜在分布区预测结果，由图2a、2b可知，杉木当前的适宜生境分布区主要集中在研究区的中西部，包括岩前镇和城东乡大部、陈大镇南部、徐碧乡西部、莘口镇北部以及中村乡西北部的部分地区.马尾松的适宜生境分布区主要集中在研究区的中部，包括城东乡和徐碧乡大部、陈大镇中南部、莘口镇北部、洋溪乡西部、岩前镇东北部以及中村乡西北部.杉木的不适宜生境分布区主要集中在研究区的东南部以及北部的零星地区，而马尾松不适宜生境分布范围较小，只在陈大镇北部和中村乡南部有零星分布.通过与收集到的三明森林二类调查小班对比发现，杉木和马尾松预测的潜在适宜生境分布区与实际分布状况基本吻合，杉木和马尾松的主体种植区域大体包含在适宜生境内.对各级适宜生境的面积及比例进行统计(表3)，杉木和马尾松在研究区的适宜生境面积分别占研究区总面积的39.00%、31.63%，低适宜生境分别占37.46%、62.41%，不适宜生境分别占23.54%、5.96%.从预测的适宜生境结果来看，研究区杉木的适宜生境分布区面积最大，杉木生长区域较宽；而马尾松的低适宜生境分布面积达到了62.14%，约是适宜生境分布面积的2倍，马尾松的适宜生长区域相对较窄.

a和b代表当前的气候条件.图2 MaxEnt模型预测当前的杉木和马尾松的潜在分布区Fig.2 Potential distribution of Chinese fir and Pinus massoniana predicted by MaxEnt model

树种 时期适宜生境面积/km2百分比/%低适宜生境面积/km2百分比/%不适宜生境面积/km2百分比/%杉木当前 322.0839.00309.3137.46194.4023.542050年238.9428.94329.4239.89257.4331.172070年247.4129.96329.2139.87249.1730.17马尾松当前 261.1431.63515.4162.4149.245.962050年208.7825.28477.0757.77139.9416.952070年214.0925.93517.0262.6094.6811.47

3.3 影响杉木和马尾松分布的主导气候因子

基于刀切法得出各环境因子对杉木和马尾松分布影响的贡献率(表4).影响杉木和马尾松潜在分布的主要环境因子是最湿季降水量(Bio16)、最冷季降水量(Bio19)、等温性(Bio3).统计环境变量贡献率(前7位)，发现与温度有关的环境因子累积计贡献率分别为25.7%、20.6%，与降水有关的环境因子累积贡献率分别为63.8%、77.3%，降水的累积贡献率远远大于温度的累积贡献率，说明降水是影响杉木和马尾松分布的最主要因素，这与赵羿涵的研究结果一致[32].其次，影响杉木和马尾松分布的环境变量前3位相同，且总贡献率均超过了60%，同时，通过单变量响应曲线分析得到各主导气候因子的阙值(存在概率>0.5)：Bio16分别为765～800、770～800 mm；Bio19分别为201～208、200～207 mm；Bio3分别为31.8～32.8、31.8～32.7，说明杉木和马尾松最适生态条件较为接近，分布区出现重叠，这与迟健[33]的研究结果较一致.

表4 影响杉木和马尾松潜在分布区的主要环境变量Table 4 Dominant environmental variables influencing potential distribution of Chinese fir and Pinus massoniana

3.4 气候变化对杉木和马尾松分布范围的预测

根据杉木和马尾松气候适宜性划分标准，得到2050、2070年不同时期下杉木和马尾松气候适宜性分布图，如图3所示.由图2和图3可知，到2050年，杉木的生境分布区有所变化，其中研究区北部陈大镇、岩前镇适宜生境分布区向南退减，徐碧乡向东扩展，研究区东南部的不适宜生境向北扩展；到2070年，杉木的各等级生境适宜性分区变化不大.到2050年，马尾松适宜生境分布区向研究区南部缩减，而研究区东南部和北部不适宜生境分布区增加较多.由表3可知，随着时间变化，到2050年，杉木和马尾松的总适生区面积(适宜生境面积和低适宜生境面积)呈减小的趋势.然而，到2070年，杉木和马尾松的总适宜区面积略有增加.杉木和马尾松边界的变化以及收缩和扩展，展现了物种对气候变化的响应过程.从单个树种来看，到2050年，杉木和马尾松适宜生境分布区面积分布下降了10.06%和6.35%，而低适宜生境分布区面积，杉木增加了2.43%，而马尾松降低了4.64%.

a、b表示21世纪50年代的未来气候情景；c、d表示21世纪70年代的未来气候情景.图3 气候变化下杉木和马尾松在三明的潜在分布格局Fig.3 Patterns of potential distribution for Chinese fir and Pinus massoniana in Sanming City under climate change in different periods

单因子响应曲线对比分析表明(表4)，影响杉木和马尾松分布的第一环境因素均为最湿季降水量(Bio16).三明雨季出现在4—6月，这3个月降水最多，达到770 mm，占全年降水量的45%左右，说明4—6月份的降水量多少对杉木和马尾松的生长有较大的影响[34].到2050年，虽然杉木和马尾松降水因子比重仍然较大，分别占总贡献率的60.4%和73.7%，但是比重分别下降了5.7%和4.7%.值得注意的是，温度因子所占比重增大.21世纪中叶前，长江及其以南地区依次增温可达1.2～1.6 ℃，而同期降水变化不大[35].结合表3面积的变化可以推测出温度的升高对杉木和马尾松的生长起到了抑制作用.2050年到2070年，降水因子比重分别增加了18.5%和12.4%，杉木和马尾松适宜生境面积均有所增加，21世纪中期以后，中国南方温度略有升高，降水量增幅较大，达到5%[35].综上所述，杉木和马尾松适宜生境面积呈增长趋势的主要原因是：降水的增加对杉木和马尾松的生长具有促进作用，而温度升高的抑制作用小于降水的促进作用.

4 讨论

本研究通过搜集整理杉木和马尾松在三明的地理分布资料，以及当前、2050和2070年的气候数据，利用基于最大熵理论的 MaxEnt 模型并结合GIS软件，综合分析了杉木和马尾松在三明的适生区(适宜生境和低适宜生境)分布的变化趋势，并分析了影响杉木和马尾松分布的主要环境因子.结果表明，在当前气候条件下，杉木和马尾松预测的潜在适宜生境分布区与实际分布状况基本吻合，杉木主要集中分布在研究区的中西部，马尾松主要集中分布在研究区中部，这验证了MaxEnt模型预测杉木和马尾松适宜分布区的可靠性.同时ROC曲线分析法验证结果表明MaxEnt模型的AUC值均大于0.9(非常接近1)，预测结果具有较高的精度.根据Jackknife检验，一些降雨和温度相关的环境变量被认为对杉木和马尾松的潜在适宜分布有重要影响.根据环境因子的贡献率确定了影响杉木和马尾松的主导气候因子：最湿季降水量、最冷季降水量、等温性，其贡献率分别达到64.9%和77%.在未来气候情景下，气温持续升高和降水量缓慢增加，使得杉木和马尾松的适生区之间发生边界变化，杉木和马尾松的适宜生境面积先缩小后增加，但适宜生境的主体分布区变化不大，在21世纪中叶以前，对于杉木和马尾松的种植要尽量避开研究区的北部以及东部，21世纪中叶以后随着降水的增多，适宜生境面积的扩大，种植区域可适当扩展.同时最湿季降水量作为杉木和马尾松分布的决定因子，这对于杉木和马尾松在林业生产以及有效的经营管理方面具有重要参考价值.

研究采用全球尺度的Worldclim数据预测杉木和马尾松在三明市两区的潜在分布，空间尺度的不一致会对研究结果产生一定影响.因此，研究通过预测结果与森林调查的二类小班数据进行对比，发现预测结果与实际分布基本一致，在一定程度上验证了研究预测结果的精度，说明预测结果可信.此外，影响杉木和马尾松分布的因子有很多，除气候因素外，还包括土地利用、土壤以及地形等.鉴于目前没有合适的不同温室气体排放情景下土壤以及土地利用等数据，且气候变化对土壤以及土地利用变化的影响尚不明确，因此本文未作考虑.为了精确的模拟气候变化条件下杉木和马尾松的潜在分布区，在今后研究中应当全面考虑杉木和马尾松空间分布对其影响因子的响应模式，进一步揭示适于杉木和马尾松分布的生物气候特性.

[1] 郭彦龙,卫海燕,路春燕,等.气候变化下桃儿七潜在地理分布的预测[J].植物生态学报,2014,38(3):249-261.

[2] DESCOMBERS P, WISZ M S, LEPRIEUR F, et al. Forecasted coral reef decline in marine biodiversity hotspots under climate change[J]. Global Change Biology, 2015,21(7):2 479-2 487.

[3] ALLEN J L, LENDEMER J C. Climate change impacts on endemic, high-elevation lichens in a biodiversity hotspot[J]. Biodiversity & Conservation, 2016,25(3):555-568.

[4] OFIPCC W G I. Climate change 2013: The physical science basis[J]. Contribution of Working, 2013,43(22):866-871.

[5] 沈永平,王国亚.IPCC第一工作组第五次评估报告对全球气候变化认知的最新科学要点[J].冰川冻土,2013,35(5):1 068-1 076.

[6] BELLARD C, BERTELSMEIER C, LEADLEY, et al. Impacts of climate change on the future of biodiversity[J]. Ecology Letters, 2012,15(4):365-377.

[7] 高文强,王小菲,江泽平,等.气候变化下栓皮栎潜在地理分布格局及其主导气候因子[J].生态学报,2016,36(14):4 475-4 484.

[8] ORTEGAHUERTA M A, TOWNSEND P A. Modelado de nichos ecológicos y predicción de distribuciones geográficas: comparación de seis métodos[J]. Revista Mexicana De Biodiversidad, 2008,79(1):205-216.

[9] PHILLIPS S J, ANDERSON R P, SCHAPIRE R E. Maximum entropy modeling of species geographic distributions[J]. Ecological Modelling, 2006,190(3):231-259.

[10] 孙敬松,周广胜.利用最大熵法(MaxEnt)模拟中国冬小麦分布区的年代际动态变化[J].中国农业气象,2012,33(4):481-487.

[11] 李丽鹤,刘会玉,林振山,等.基于MAXENT和ZONATION的加拿大一枝黄花入侵重点监控区确定[J].生态学报,2017,37(9):3 124-3 132.

[12] 雷军成,王莎,王军围,等.未来气候变化对我国特有濒危动物黑麂适宜生境的潜在影响[J].生物多样性,2016,24(12):1 390-1 399.

[13] 张东方,张琴,郭杰,等.基于MaxEnt模型的当归全球生态适宜区和生态特征研究[J].生态学报,2017,37(15):5 111-5 120.

[14] 唐继洪,程云霞,罗礼智,等.基于Maxent模型的不同气候变化情景下我国草地螟越冬区预测[J].生态学报,2017,37(14):4 852-4 863.

[15] HIJIMANS R J, CAMERON S E, PARRA J L, et al. Very high resolution interpolated climate surfaces for global land areas[J]. International Journal of Climatology, 2005,25(15):1 965-1 978.

[16] XIN X G, WU T W, LI J L, et al. How well does BCC_CSM1. 1 reproduce the 20th century climate change over China?[J]. Atmospheric and Oceanic Science Letters, 2013,6(1):21-26.

[17] 石彦军,任余龙,王式功,等.BCC_CSM气候模式对中国区域气候变化模拟能力的检验[J].高原气象,2012,31(5):1 257-1 267.

[18] 应凌霄,刘晔,陈绍田,等.气候变化情景下基于最大熵模型的中国西南地区清香木潜在分布格局模拟[J].生物多样性,2016,24(4):453-461.

[19] MORENO R, ZAMORA R, MOLINA J R, et al. Predictive modeling of microhabitats for endemic birds in South Chilean temperate forests using maximum entropy (maxent)[J]. Ecological Informatics, 2011,6(6):364-370.

[20] 许志东,丁国华,刘保东,等.假苍耳的地理分布及潜在适生区预测[J].草业学报,2012,21(3):75-83.

[21] 张路.基于MAXENT模型预测齿裂大戟在中国的潜在分布区[J].生物安全学报,2015,24(3):194-200.

[22] 王运生,谢丙炎,万方浩,等.ROC曲线分析在评价入侵物种分布模型中的应用[J].生物多样性,2007,15(4):365-372.

[23] HANLEY J A, MCNEIL B J. The meaning and use of the area under a receiver operating characteristic (ROC) curve[J]. Radiology, 1982,143(1):29-36.

[24] 胡理乐,张海英,秦岭,等.中国五味子分布范围及气候变化影响预测[J].应用生态学报,2012,23(9):2 445-2 450.

[25] 张海娟,陈勇,黄烈健,等.基于生态位模型的薇甘菊在中国适生区的预测[J].农业工程学报,2011,27(13):413-418.

[26] PETERSON A T, COHOON K P. Sensitivity of distributional prediction algorithms to geographic data completeness[J]. Ecological Modelling, 1999,117(1):159-164.

[27] 许国.Jackknife 估计及其应用浅探[J].山东师范大学学报(自然科学版),2000,15(4):454-455.

[28] 路春燕.基于GIS与Fuzzy的野生植物生境适宜性评价与区划研究[D].西安:陕西师范大学,2012.

[29] 孙立,徐海根,吴军,等.气候变化对我国白枕鹤繁殖地分布的影响[J].生态与农村环境学报,2014,30(5):595-600.

[30] GRAHAM C H, FERRIER S, HUETTMAN F, et al. New developments in museum-based informatics and applications in biodiversity analysis[J]. Trends in Ecology & Evolution, 2004,19(9):497-503.

[31] WISZ M S, HIJMANS R J, LI J, et al. Effects of sample size on the performance of species distribution models[J]. Diversity & Distributions, 2008,14(5):763-773.

[32] 赵羿涵,刘兰妹,高露双,等.气候变化后杉木人工林单木径向生长规律研究[J].西南林业大学学报,2013,33(5):9-16.

[33] 迟健.马尾松和杉木生态习性的比较[J].湖南林业科技,1992,19(2):31-35.

[34] 姜亚林.三明沙溪下游地区年降水量系列代表性及多年变化分析[J].吉林水利,2008(7):43-46.

[35] 江志红,张霞,王冀.IPCC-AR4模式对中国21世纪气候变化的情景预估[J].地理研究,2008,27(4):787-799.