气候变化情景下毛竹潜在分布及动态预测
2021-07-13陈禹锦罗喻才杨光耀张文根
陈禹锦 罗喻才 于 芬 杨光耀 张文根
(江西省竹子种质资源与利用重点实验室 江西农业大学 南昌 330045)
气候影响植物生长，决定植被类型及其地理分布，同时植被是反映气候变化最鲜明的标志[1]。基于地理信息系统(GIS)技术和最大熵原理构建的MaxEnt模型，可预测物种在未来气候环境下的生境适宜分布区及其扩张或收缩情景[2-3]。近年来，该模型在群落生态学与生物多样性研究中得到了广泛关注与应用[4]。例如：气候变化对物种分布的影响研究，如海南凤仙花(Impatienshainanensis)[5]、南酸枣(Choerospondiasaxillaris)[6]、小黄花茶(Camellialuteoflora)[7]等；物种入侵监测研究，如预测红火蚁(Solenopsisinvicta)在中国的适生区分布[8]；自然保护区的规划设计，如白冠长尾雉(Syrmaticusreevesii)保护区功能区的重新划定[9]等。
毛竹[Phyllostachysedulis(Carriere) J. Houzeau][10]是中国栽培历史久、种植面积广、利用价值高的竹种，在材料、食品、康养、绿色环保等领域均有着巨大的应用潜力[11-13]。然而，毛竹入侵降低森林群落中乔木和灌木的丰富度，破坏森林植被，改变土壤养分循环模式，降低碳储量和改变土壤微生物结构，在中国南方已成为突出的生态学问题[14-16]。目前，关于毛竹环境适宜性分布、适生区预测的研究报道较少[17-18]。本文首次利用Maxent模型并结合GIS技术对气候变化情景下毛竹的潜在适生区进行了预测，建立了毛竹生境的变化情景，可为毛竹种植、引种以及入侵防治提供参考。
1 材料与方法
1.1 毛竹地理分布点数据
基于中国现有分布区内有代表性的34个毛竹居群[19]，结合中国数字植物标本馆(CVH: http://www.cvh.ac.cn/)、国家标本资源共享平台(http://www.nsii.org.cn/2017/home.php)和记录有毛竹分布地点的文献资料，经筛选和剔除，获得毛竹分布点数据316条，利用谷歌Earth卫星图拾取经纬度坐标并导入Excel中，保存为.csv格式。
1.2 环境变量的采集与处理
1.2.1 环境变量采集
从Worldclim(www.worldclim.org)获取当前时段(1981—2000年)及未来时段(2021—2040年、2061—2080年)的19个气候生物变量，分别为：年平均温(bio1)、昼夜温差月均值(bio2)、等温性(bio3)、温度季节变化(bio4)、最暖月最高温度(bio5)、最冷月最低温度(bio6)、年均温变化范围(bio7)、最湿季度平均温度(bio8)、最干季度平均温度(bio9)、最暖季度平均温度(bio10)、最冷季度平均温度(bio11)、年均降水量(bio12)、最湿月降水量(bio13)、最干月降水量(bio14)、降水量变异系数(bio15)、最湿季度降水量(bio16)、最干季度降水量(bio17)、最暖季度降水量(bio18)和最冷季度降水量(bio19)。
未来气候变量是基于联合国政府间气候变化专门委员会(IPCC)发表的第6次评估报告对全球未来气候变化所进行的预估而得出的。大气环流模型采用中国国家气候中心开发的BCC_CSM模型[20]。未来气候情景包括了4种典型浓度路径(SSPs)，分别为：SSPs126、SSPs245、SSPs370、SSPs585。本文选用SSPs126、SSPs245、SSPs585分别代表低、中、高温室气体排放情景。
土壤因子来自于寒旱区科学大数据中心(http://data.casnw.net/portal/)的中国土壤数据集v1.1(China Soil Map Based Harmonized World Soil Database, v1.1)。土壤质地和结构、容重、水分含量、酸碱度、钙离子含量等因素[21-23]均为限制毛竹生长的因子，因此选用土壤质地分类(bio20)、土壤淤泥含量(bio21)、土壤容重(bio22)、土壤基本饱和度(bio23)、土壤酸碱度(bio24)和土壤碳酸盐含量(bio25)数据。
地形因子选用海拔(bio26)、坡度(bio27)和坡向(bio28)数据。DEM高程数据来自中国科学院地理空间数据云(http://www.gscloud.cn)。高程图像通过Arcgis 10.2中spatial analyze工具提取，坡度、坡向图层由高程图层通过Arcgis 10.2中3D Analyst工具提取。
1.2.2 环境变量处理
所有环境变量使用WGS-1984地理坐标系，按照来源于国家基础地理信息系统(http：//nfgis.nsdi.gov.cn/)的中国行政区划图剪裁，空间分辨率调整为2.5′(22.5 km2)。为避免环境变量之间的过度拟合[24]，保证模型运行的准确性，将28个环境变量及毛竹分布点数据导入MaxEnt 3.4.1中进行预模拟，初步得出各环境变量的贡献率。使用Arcgis 10.2中的重采样工具将28个环境变量与毛竹分布点耦合并导出Excel表，在Spss中采用双变量相关并勾选Pearson检验进行相关性分析。若2个环境变量的相关性≥0.9，则保留预模拟中贡献率较大的环境变量。最终，确定17个环境变量为主要影响毛竹生长的环境因子，并进行最终预测得出各环境变量的贡献率及置换重要值(表1)。环境变量统一采用ASCⅡ格式。
表1 环境变量对毛竹分布预测的贡献率及置换重要值Tab.1 The contribution and permutation importance of environmental variables to the predicted distribution of P. edulis.
1.3 毛竹分布区预测方法
将毛竹分布点数据和筛选出的17个环境变量分别导入MaxEnt 3.4.1的物种分布(Samples)模块和环境变量(Environmental layers)模块。选用75%的毛竹分布点作为训练集用以建立模型，25%作为测试集用以测试模型，分布值以逻辑斯蒂(Logistic)格式输出。重复模拟10次，最终输出的ASCⅡ格式图层为10次重复的平均值。Liu等[25]对现有的13种选择阈值的方法进行过深入研究，本文选用最大特殊敏感性加特殊性(Maximum training sensitivity plus specificity)为该模型的阈值。
1.4 模型精确度评价
使用刀切法(Jackknife)计算各环境变量对毛竹分布区域影响的贡献率。以假阳性(1-特异度)为横坐标，真阳性(灵敏度)为纵坐标绘制接受者操作特征曲线(ROC曲线)。模型的优劣由ROC曲线下方的面积AUC值决定。AUC的数值范围为0～1，数值越大表示模型准确度越高。根据Swets的标准[26]，AUC值为0.5～0.6时模型不存在预测能力，AUC值为0.6～0.7时模型预测能力一般，AUC值为0.7～0.8时模型预测能力良好，而AUC值>0.9时模型准确性非常高。
1.5 适生等级划分
选取MaxEnt重复运算10次所得到的最大特殊敏感性加特殊性平均值0.198为毛竹适生区分布概率的阈值p，p<0.198为非适生区，p>0.198为潜在适生区。将预测结果导入Arcmap 10.2中，使用重采样工具以自然断点分级法将图层划分为非适生区(p<0.198)、低适生区(0.198
1.6 未来不同气候背景下毛竹分布区域的动态变化
为对比未来气候情景下毛竹适生区面积与当前适生区面积的变化趋势，以0.198作为阈值，指示适生(p>0.198)与非适生(p<0.198)2种状态，绘制二进制图层。将当前与未来毛竹适生区的二进制图层导入Arcmap中，进行叠加、运算，获得毛竹适生区的动态变化图。
2 结果与分析
2.1 模型精确度检验
采用ROC曲线下方面积AUC值检验模型精度。经过10次运算，当前适生区模型的平均AUC值为0.912(图1)。根据Swets的标准可知，AUC值>0.9时模型准确性非常高，说明模型能够准确预测毛竹的适生区域。
图1 毛竹MaxEnt模型的ROC曲线Fig.1 ROC curve of the MaxEnt model in P. edulis
2.2 影响毛竹分布的气候因子
基于贡献率、置换重要值和刀切法检验并选择限制毛竹在中国潜在地理分布的环境变量因子。贡献率较大的前4个环境因子变量分别为最干月降水量(bio14，73.3%)、年均降水量(bio12，5.9%)、最冷月最低温(bio6，5.2%)和温度年较差(bio7，3.5%)。置换重要值的大小表示模型对该变量依赖性的强弱，置换重要值较大的前4个环境因子变量为海拔(bio26，29.9%)、最干月降水量(bio14，12.9%)、最冷月最低温(bio6，11.5%)和年平均气温(bio1，9.5%)。根据刀切法所得数据(图2)可知，训练模型得分较高的前4个环境因子变量分别为最冷月最低温(bio6)、最干月降水量(bio14)、年均降水量(bio12)和年平均气温(bio1)。
图2 各环境变量对毛竹当前潜在分布模型的Jackknife检验得分Fig.2 Jackknife test scores of environmental variables on the current potential distribution model of P. edulis.
由此可知，主导毛竹潜在地理分布的是降水变量(最干月降水量、年均降水量)、温度变量(最冷月最低温、温度年较差、年平均温)和地形变量(海拔)。土壤因子虽然贡献率、置换重要值较小，但对毛竹的生长有不可替代的作用。绘制以上6个环境变量的响应曲线(图3)以表示毛竹分布概率与环境变量之间的关系。当毛竹适生区分布概率大于阈值(0.198)时，限制毛竹分布的主导环境因子所在区间分别为最干月降水量(10.0～183.1 mm)、年均降水量(852.4～3 128.6 mm)、最冷月最低温(-4.7～15.7 ℃)、年平均温(13.0～23.8 ℃)、温度年较差(<34.5 ℃)和海拔(<1 315.0 m)。
图3 毛竹对主要环境因子的响应曲线Fig.3 Response curve of P. edulis to major environmental factors.
2.3 当前毛竹潜在分布区
在当前气候条件下，毛竹实际地理分布点均位于潜在分布区范围内(图4)，预测结果与实际分布情况相吻合。毛竹潜在适生区的总面积为191.61万km2，主要位于亚热带季风气候区。分布范围为18.31°—37.47° N，93.46°—122.47° E，南起于海南岛的南部，北至陕西南部(越过秦岭，分布于秦岭北麓)、河南南部及山东沿海一带；西起于西藏东南部、四川的东部和云南中南部，东至江苏、浙江、福建沿海一带及台湾地区。其中，高适生区面积为24.75万km2，主要位于总适生区的东部和西北部，即江西中北部、湖北西部、安徽南部、江苏沿海区域、重庆西部和四川东南部等地，呈间断分布；中适生区面积为89.43万km2，主要位于总适生区的东部，即江西南部、福建大部、浙江、安徽中南部，呈连续分布；低适生区面积为77.43万km2，主要位于总适生区北部、西部和南部的边缘区域，呈连续分布。
图4 毛竹地理分布点及当前潜在适生区预测Fig.4 Geographical distribution and prediction of current potential suitable area for P. edulis.
2.4 未来毛竹潜在适生区及其动态变化
2.4.1 未来毛竹潜在适生区面积变化
在BCC_CSM大气环流模型、3种温室气体排放情景、未来2个时间段条件下，模拟得到6个毛竹未来分布模型(图5)。与当前潜在适生区面积相比(表2)，未来潜在适生区面积均有不同程度地缩减，不同等级适生区面积的变化幅度较大。高适生区面积在2021—2040年时间段内将大幅缩减，到2061—2080年时间段内将进一步缩减(缩减幅度>90%)，并且高适生区面积随温室气体排放情景由低到高变化而依次增大。中适生区面积将在2021—2040年时间段内扩张，到2061—2080年时间段内中适生区面积明显随温室气体排放情景由低到高变化而减小。低适生区面积将在2021—2040年时间段内呈现先扩张后收缩的变化趋势，在2061—2080年时间段随温室气体排放情景由低到高变化明显呈现扩张趋势(扩张幅度>30%)。
图5 未来气候情景下的毛竹潜在适生区Fig.5 Potential suitable area of P. edulis under future climate scenarios
表2 不同气候变化情景下毛竹适生区面积变化Tab.2 Suitable area change of P. edulis under different climate change scenarios
2.4.2 毛竹未来潜在适生区范围变化
从毛竹潜在适生区分布图得知，在2021—2040年时间段内，在3种温室气体排放情景下，毛竹的中、低适生区范围无明显变化；高适生区变化范围随温室气体排放情景由低到高变化依次扩大，无明显迁移趋势。而在2061—2080年时间段内，广东、广西境内的低适生区范围向北移动，且随温室气体排放情景由低到高变化，低适生区南边界向北移动范围越广。四川东部、重庆西部、湖南南部、广东北部、广西北部的中适生区范围随温室气体排放情景由低到高变化而依次减小，中适生区呈破碎化分布。高适生区范围大幅度缩减：在SSPs126模式下，仅在浙江沿海区域有小范围分布；在SSPs245模式下，仅在江苏南部、浙江北部及沿海地区有小范围分布；在SSPs585模式下，仅在浙江北部有小范围分布。从当前到2061—2080年时间段，毛竹高适生区有向东迁移的变化趋势。
2.4.3 未来毛竹潜在适生区动态分布
毛竹潜在适生区在未来气候情景下呈现向高纬度地区转移、整体北移的变化趋势。由动态分布图(图6)可知，毛竹适生区的南部除海南岛存在面积扩张现象以外，其余地区面积均缩减，南界整体呈现北移趋势。适生区的北部面积向北扩张，北界北移。其中，西北部有极小程度的扩张，东北部(山东、河南)适生区面积明显扩张，东、西部适生区范围无明显变化趋势。
图6 未来气候情景下毛竹潜在适生区动态变化Fig.6 Dynamic change of potential suitable area for P. edulis under future climate scenarios
2021—2040年时间段内SSPs126排放路径下，以及2061—2080年时间段内SSPs245、585排放路径下，毛竹潜在适生区在新疆伊宁市西部地区出现小范围的扩张。
3 讨论
本文基于MaxEnt预测模型，利用当前气候、土壤、地形数据及毛竹分布地点预测毛竹当前潜在适生区范围，并在3种温室气体排放情景下预测毛竹在未来2个时间段内(即2021—2040年，2061—2080年)潜在适生区范围及动态变化过程。由上文所得结果可知，毛竹的地理分布主要受到水热条件的制约，限制毛竹分布最重要的环境变量为最干月降水量、年降水量、最冷月最低温、年平均温、温度年较差和海拔。黄韬等[27]研究表明，越冬期间的降水天数及温度是毛竹育笋的重要限制因子，若越冬期气温低、降水量小，则对竹笋萌长不利，从而导致来年毛竹成株数量减少。该研究与本文预测结果相一致，可验证最干月降水量、最冷月最低温是限制毛竹地理分布的重要环境因素。汪阳东等[28]研究认为，年降水量对毛竹生长有显著作用，年平均温对毛竹生长有极显著作用，与本文结果也相一致。海拔通过影响土壤理化性质、土壤酶活性、光照条件和气温进而影响毛竹的生长，从而限制毛竹的垂直地理分布。毛竹在山地地貌条件下分布的海拔上限是1 200～1 500 m[29]，与本文海拔响应曲线所得出的上限1315 m大致相同。
毛竹当前潜在适生区范围为18.31°—37.47° N，93.46°—122.47° E，面积为191.61万km2。在未来气候情景下，毛竹适生区面积部分缩减，高适生区面积缩减幅度大，适生区范围整体呈现北移的趋势。王佳鑫等[17]所得的当前毛竹适生区面积为246万km2，差异主要在于未考虑地形变量，因此我们推断地形变量尤其是海拔变量对毛竹地理分布有明显的限制作用。与张雷[18]的预测结果相比，本文中当前潜在毛竹适生区范围更广，在未来气候情景下北移幅度明显较小。造成结果差异的原因主要有：1)2001年中国植被图所反映的是20世纪80和90年代的植被分布现状，而毛竹根系发达，入侵速度快，地理分布区域在较短时间尺度内会有较大变化，从2001年中国植被分布图中提取的毛竹分布数据不能很好地反映当前毛竹实际地理分布；2)未考虑土壤变量和地形变量对毛竹生长的影响。金佳鑫等[17]研究表明，土壤因素对于毛竹潜在适生区在未来气候情景下北移的幅度具有明显的限制作用。
预测结果对开展实际工作具有指示作用。毛竹当前潜在地理分布区主要位于中国东南部，是典型的亚热带季风气候区，该区域雨量充沛，气候温暖，水热条件良好，与毛竹的生长习性相适应。当前实际地理分布点范围远小于毛竹潜在适生区范围，并且我国自20世纪在云南、陕西、山东等地开展的毛竹引种实验均已取得成功[30-32]，说明在当前潜在适生区内进行毛竹的引种栽培具有一定的可行性。在未来气候情景下毛竹适生区具有北移的可能性，建议根据本文所得出的毛竹适生区变化趋势及适生区分类等级并结合具体的社会、经济、环境因素选择适宜地区进行毛竹引种栽培。
毛竹是典型的克隆植物。通过竹鞭进行无性繁殖向其他林分扩张[33]，由阔叶林过渡为竹阔混交林最终成为毛竹纯林，威胁其他物种的生存，降低了其生态功能[34]。目前，毛竹入侵产生的生态危害已成为热点研究问题。本研究获得毛竹当前及未来的潜在分布范围及适生程度，结合毛竹分布范围的时间尺度与空间尺度可知：在全球气候变暖的大前提下，毛竹适生区范围有随时间推移而北移程度加剧的变化趋势，可为未来适生区面积扩张的地理区域提供预警。应在毛竹当前潜在适生区内及未来可能扩张的区域内密切监测毛竹扩张动向，结合该区域的森林植被类型、物种组成、群落结构等特性开展有针对性的管理、清除工作，保护该地区生态功能免受破坏。总之，应权衡毛竹引种带来的经济效益与毛竹扩张造成的生态危害，丰富毛竹入侵管控手段，控制竹林边界，局限毛竹生长范围，实现竹林的精细化管理，在实现经济价值的同时保持生态环境良好、生态功能稳定。
