贡嘎山竹类资源分布模型模拟研究

李志明，王 乾,余国宝，陈绪文*

(1.四川贡嘎山国家级自然保护区石棉管理处，四川 石棉 625400;2.中国科学院成都生物研究所，四川 成都 610041)

生物多样性是区域性生态—进化过程的外在表现，是生态系统功能的基础，是维系人类社会生存和发展的根本[1]。目前全球生物多样性在以前所未有的速度消失，即使作为旗舰性保护物种的大熊猫，也受到栖息地生境破碎化的不良影响。四川盆地西缘山地是大熊猫主要分布区域，也是箭竹类竹种的分化中心和现代分布中心[4]，大熊猫保护很大程度上得益于这里适宜的生境条件和竹类资源[2]。箭竹类作为一种林下和灌丛中的木本植物资源，不仅是当地生物多样性的重要组成部分，还是重要的经济植物，其种类和生存状态对生态系统平衡有着重要的作用[3]。竹类主要靠无性克隆繁殖，营养期长达几十年，扩散能力较差[5]。而区域内竹子一旦开花多会同步进行，并随之造成大片竹林死亡，导致大熊猫出现食物短缺甚至死亡。因此竹类分布对环境反应具有滞后性和周期性，环境波动过大容易导致局部灭绝。

气候是决定物种分布的顶层驱动力，其次为生境类型和生态过程，尺度从大到小通过层层筛选，物种的实际分布范围要小于气候适宜区[9]。熊猫保护区天然竹类主要分布在林下，从上层林冠类型到地面基质状况，从种间竞争到演替阶段以及干扰历史无不对其分布起着约束作用。气候变暖速率超过物种扩散速率也会造成竹类生境范围和资源的大幅减少[6～8]。因此详细研究竹类分布与气候和环境的关系特征对于了解竹类分布状况，掌握承载力，及早发现敏感区域，开展科学合理的保护措施有重要意义。目前关于竹类的空间分布研究仅有少数描述性的资料发表，没有具体分布位点，不能确定相应的环境因素，或者涉及保密只给出整个熊猫保护区的分析结果[7]，这种结果涉及范围大，气候、植被等环境因素的空间异质性较强，难以反映物种的环境需求。

本文应用四川省第1高山——贡嘎山来研究竹类资源分布与生态模型尺度的关系。采用贡嘎山竹类调查分布位点和WorldClim生物气候数据，利用集合模型构建了贡嘎山区3种竹类的分布预测，以期为保护区分区规划和生物多样性保护工作的依据。

1 研究区域

贡嘎山是四川省第1高山，也是青藏高原东部的最高山，整个山区东部为高山峡谷地貌，向西则与川西高原相接，东西气候差异显著[10]，动植物资源极为丰富，植被类型多样，空间分布格局显著[11]。河谷地带为亚热带干旱河谷气候区，山顶则为冰雪永冻带，东坡山势陡峭，降水丰富，植被茂密，植被垂直带谱丰富，类型多样，物种丰富。山体整体处于中亚热带向北亚热带过渡区域，基带性植被基本稳定，气候梯度由山体地形地貌格局造成，南北环境差异不显著，由于整体范围不大，可以认为环境因素的响应在尺度上是不变的，因此适合分析物种的环境需求。

贡嘎山西部为山地向高原过渡的山原地貌区，气候干冷、辐射强烈，植被以亚高山针叶林和高山灌草丛为主。贡嘎山自然保护区环绕贡嘎山主峰及其周边山峰，包括石棉、泸定、康定和九龙4个县的部分地区，总面积约 4 091 km2。贡嘎山位于横断山脉东部，在地理环境上为水热转换地带，是东西向和南北向上的生物区系过渡地带，具有重要的科研价值。贡嘎山同时也是大熊猫3个种群(小相岭、大相岭和邛崃山)的连接区域，对保持熊猫遗传多样性的完整性起着关键作用。而竹类作为熊猫的食物在廊道的质量上的作用尤为重要。我们2017年的调查发现整个分布区有竹子种类8种，其中天然竹类6种，包括箭竹属3种，巴山木竹属两种，玉山竹属1种。我们选用其中分布最广的3个竹种做空间分布研究，即丰实箭竹(Fargesiaferax)、冷箭竹(Bashaniafangiana)和峨热竹(Bashaniaspanostachya)。

2 材料与方法

竹类分布位点数据来自于带GPS定位功能的相机拍摄的竹类微距照片，照片拍摄于2017年4月和9月的贡嘎山竹类调查工作中。共拍摄两万余张照片，从中选出1500个竹类位点。去除重叠的位点后，共筛选出320个竹种分布点，各竹种的分布位点在预测图中标出，以资参考。调查区域基本涵盖贡嘎山区可通行的区域。

气候数据来自于WorldClim数据库[12]的分辨率约1 km的空间插值生成的生物气象数据，共包括19个变量，由于数据之间可能存在共线性问题，造成模型预测效果较差，因此需要事先将共线性较强的数据剔除。依据成对Pearson相关系数矩阵对这19个生物气候变量进行选择，去除平均相关系数绝对值大于0.7的变量[7]，最后保留6个变量，即bio2、bio3、bio4、bio12、bio15、bio17变量，其中前3个变量与温度有关，后3个与降水有关。所有分析使用biomod2包[13]在CRAN R软件中进行。Biomod2包中包括了多种模型分析，其中集合模型(Ensemble Modeling)将模型中的10个模型(GLM,GBM,GAM,CTA,ANN,SRE,FDA,MARS,RF,MAXENT)作为子模型进行物种分布建模，然后利用投票方式，以票数过半来决定物种的在每个空间子区上的分布情况。这10个模型包括广义线性模型、广义增强模型、广义加性模型、分类树分析、人工神经网络、生物气候模型、柔性判别分析、多重适应性回归样条、随机森林以及最大熵模型等。这些模型各自在某一方面表现较强，综合的结果通常较单一模型表现更优。

调查的竹类分布位点数据被分成两部分，其中80%用于模型训练，其余部分用于模型质量评估，采用5折交叉验证，计算模型预测精度、敏感性和专一性，TSS(true skill statitic)[14]进行模型比较。气候数据为栅格数据，空间分辨率30″，在贡嘎山区约890m，坐标系为WGS84经纬度格式。通过竹类分布位点计算对应栅格点上的各气象要素的值分布情况，然后应用各种概率模型计算最优的反应曲线。最后反演到整个区域，找到相应的分布范围。

3 结果

单个模型中以广义线性模型表现最差，随机森林模型表现最好，TSS值大于0.85，说明模型的敏感性和专一性都较高。集合模型通过对10种模型的预测投票，以少数服从多数的方法，博采众家之长，得到预测结果，TSS值均大于0.90，表现效果超出单个模型的最好效果。

预测的竹类基本沿着沟谷分布，包含着实际分布点位，各竹种之间呈现一定的空间分布格局(图1)。丰实箭竹分布海拔较低，集中分布在泸定县磨西镇和石棉县草科乡地段，南部九龙县洪坝乡也有一定的分布。康定县城附近有零星的分布。沿贡嘎山主峰经线以西区域均无丰实箭竹分布。冷箭竹分布区域较为狭窄，主要分布在磨西镇和康定县城附近海拔较高的有林地段，从康定折多塘至贡嘎山主峰一线以西地段均无冷箭竹分布。峨热竹分布较为广泛，除西南部和东部海拔较低处外，基本都有分布，其中草科乡以西的沟谷中、康定县城附近、洪坝乡附近分布较多，在保护区之外康定县朋布西和六巴的东部地区分布较广。

图1 使用集体模型构建的贡嘎山区3种竹类的空间分布预测图(A为丰实箭竹，B为冷箭竹，C为峨热竹)，白色方块构成的区域为预测分布区，红色叉为实际调查分布位点，灰白色线包围区域为贡嘎山自然保护区范围，图片背景为海拔高程图。Fig.1 The predicted distribution map of three bamboo species using the ensemble model(A:Fargesia ferax; B:Bashania fangiana; C:Bashania spanostachya).The area composed of white squares are predicted distribution area of the specific species.The crosses are actual present records of the species.The natural reserve of the Gongga Mountains are surrounded by the pale white lines.The background is a digital elevation map.

实际调查竹类的海拔分布以丰实箭竹最低，中位数海拔为 2 200 m；其次为峨热竹，中位数海拔为 3 100 m；最高的是冷箭竹，中位数海拔为 3 400 m(图2)。3个竹种中均存在海拔分布特别低的异常值，但竹种分布海拔之间t统计检验存在显著的区别 (p<0.001)。

图2 调查竹种的实际海拔分布箱形图。Fig.2 The boxplot of elevational distribution of the surveyed bamboo species.

但是预测区域的海拔分布与实际调查分布有一定的区别，例如丰实箭竹的海拔中位数为 2 340 m，比实际提高140 m。预测的海拔最低值 1 050 m，最高值达到 3 770 m，远远超过实际调查的海拔高度范围。峨热竹预测区的海拔中位数为 3 910 m，比调查的中位数提高810 m，海拔幅度在 1 740 m～4 630 m之间，也远超调查的范围。冷箭竹预测区海拔中位数为 3 320 m，与实际调查的中位数最为接近，但是预测的海拔范围(2 170 m～4 630 m)与实际调查情况差异也很大。

在保护区内，各竹种的预测分布面积以峨热竹最大，达到231.7 km2，其次为丰实箭竹，达到170.6 km2，冷箭竹面积最小，为134.9 km2。保护区中这3种竹类的分布面积只占13%。

4 讨论

4.1 气候异质性

贡嘎山东西坡降水呈现不同的模式[10,11]，东坡受东南季风影响强烈，随着海拔升高，降水量逐步增多。由于贡嘎山地形抬升，形成地形雨，在迎风坡的一定海拔范围内降水随海拔增加而增加。其增加幅度从泸定县城到磨西镇为每百米131 mm，从磨西镇到海螺沟索道(海拔 3 000 m)为67 mm，总降水量达到约 2 000 mm[15]。贡嘎山西部地区受西南季风和东南季风的双重影响，降水模式为双峰形式，最大降水量约出现在6月和9月。贡嘎山西南部的九龙县城降水量虽大于贡嘎山北部的康定县城，但也只有同海拔的东部区域的一半。除了降水之外，同海拔地区的温度也有显著的区别，东部显著低于西部地区。海螺沟年均温4.2℃，九龙县为10.7℃。贡嘎山东坡林线位置多在海拔 3 600 m左右，而九龙县汤古瓦灰山林线海拔在 4 300 m左右，两者相差700 m。而林线的位置被认为具有相同的最热月温度或者生长季温度[16]，因此东部气温的垂直递减速率显著大于西部。

我们用相关系数法去除了大部分气候变量，只保留了共线性较弱的6个变量，这样可以增强模型的稳健性和移植性，对附近区域的竹类分布预测有指示意义。例如我们发现康定朋布西和六巴以东区域预测有大片的峨热竹分布，将来的实地探查可以验证模型的可靠性。

这种水热分布上的显著区别在山势陡峭、地形复杂的东部区域往往造成复杂的地形气候模式，用分辨率粗糙的气候模型模拟会掩盖许多特殊气候模式。目前用于预测未来气候变化的数据空间分辨率较为粗糙，基本采用WorldClim的约1km精度数据，其精度受限于站点密度和地形特征[10]，这就造成了显著的取样偏差。在平原和丘陵地区使用不会造成太大的误差，但是在川西高山峡谷区由于气象站点稀少、地形起伏剧烈以及地形气候效应较为突出的地方，使用粗分辨率数据往往造成较大的误差。例如海拔每升高1 km，温度大约下降6.5℃[16]，按照45°的山坡坡度计算，1km水平分辨率的温度数据在峡谷区的误差就已经超过了IPCC估计的未来几十年的温度升高幅度。在这种尺度下的模拟研究可能会夸大变化范围，误导生物多样性的研究和保护。而且这种预测区域所覆盖地形海拔不可避免会超过实际的分布海拔，其预测的可信度降低。因此粗分辨率数据只能给出平均的环境特征，不能反映镶嵌分布的面积较小的特殊生境特征。

在物种分布的海拔预测上，使用粗分辨率的气候数据导致海拔出现较大偏差，特别是在地形陡峭的地段，预测范围往往超出竹类可能分布的海拔上限。因此粗分辨率数据只能估计平原等与平均状态较为接近的区域，在所研究的竹类中冷箭竹和丰实箭竹的实际海拔分布的中位数与预测海拔的中位数分布最接近。但是和峨热竹的海拔中位数差异达到810 m，这种情况可能和西北部大片出现的新预测分布区的海拔较高有很大关系。

4.2 物种特征

竹子的分布和植被也有一定的关系。竹子通常分布在稀疏的林下或者林窗内，要求一定的光照，因此其分布状况和森林的演替阶段有一定的关联。密林下或者人工林下由于林内光线较弱，竹类通常较难生存。冷箭竹和峨热竹生长在高海拔的针阔混交林内，特别是桦木云冷杉林下，竹子常是灌木层的优势种，呈散生状态，高度可达2 m左右。丰实箭竹常生长在针阔混交林和落叶阔叶林内，竹丛密集，竹体高大，通常高度在3 m以上。在高山灌丛和高山草地中未见竹子分布。在高海拔段的稀疏冷云杉林中，香柏、红杉、柳树、高山栎、杜鹃等物种的盖度较大，冷箭竹和峨热竹通常呈低矮的禾草状，并且较为稀疏，当地俗称“麦秧子”，不注意很难识别。

竹子对温度和光照的需求是有范围限制的，在高海拔地段，通常一个因素就可以决定其生存。我们在九龙汤古调查发现，高海拔段阳坡虽然光照充分，温度也相对较高，在和有竹类分布的阴坡同一海拔高度的阳坡或者半阳坡上却没有竹类分布。模型中如果加入植被和坡向等数据可以进一步缩小预测偏差。

物种的生存除了要满足气候条件外(基础生态位)，还需要满足土壤、植被、生物互作等因素，因此实际分布区小于基础生态位[9]。而且生物气候关系在外推时还有关系发生变化的可能，距离越远，这种变化的危险性越大。收集高质量的物种分布数据、采用更高精度的气候数据、加入更多的环境变量对于准确预测物种实际生态位会起到积极的推动作用。

4.3 模型的适用性

本文采用集合模型来预测物种分布，包括10种模型，每种模型都将数据拆成两份，其中一份用于训练，另一份用于验证，其结果可以通过敏感性和专一性的测量来进行有效性判别。数据的划分是采用随机抽样的方法生成的，对多次抽样生成的数据集进行模拟，然后求平均值，得到该种模型相应的分布预测。相比于单用某一模型，集合模型通过模型投票决定最终的预测结果，这样可以防止模型的过拟合以及弥补模型结构差异带来的缺陷，从而获得一个较为稳定的预测结果。竹类分布海拔存在一些异常值，与未检测到相反，这些是确定存在的物种分布点，而未检测到却是实际存在，但是没有发现的情况。因此可能存在海拔分布的连续性，但是由于没有调查或者调查了却没有发现而造成海拔分布的间断性。集合模型将这些异常值多做抛弃处理，可以看到检测点处没有预测分布。这种问题可能是取样偏差造成的，均匀布设样地开展取样工作可以减少这种问题。

集合模型的基础模型都是属于相关模型，通过关联找到植物与气候之间的反应曲线，进而预测到全区域适宜生态位。我们的研究只是考虑了气象因素，而这种数据集的相关性随着区域发生改变，例如Li等[7]在整个熊猫分布区的分布研究中，发现有8个气象因素满足要求，而我们则只找到6个。说明前者环境异质性更强，大区域中变量方差变大，导致相关性降低。但是局部区域变量之间的相关性依然很强，这会导致模型预测效果变差。因此将研究区域细分成相对均匀的小块做模型分析，能够更加准确地预测物种的分布。