祁连山北坡青海云杉的潜在分布

涂振宇，勾晓华，邹松兵

(兰州大学资源环境学院，西部环境教育部重点实验室，甘肃 兰州 730000)

祁连山位于我国青藏高原东北缘，由多条西北—东南走向的平行山脉组成。祁连山南北坡地形差异明显，南坡地势相对平缓，北坡海拔变化较大，主要分布在2 000～5 500 m的范围内[1-2]。祁连山森林生态系统具有重要的生态功能，是我国西北地区的重要生态屏障，但在21世纪初期，森林砍伐、过度放牧和开矿开采等行为不断发生[3]，森林面积不断减少，水源涵养能力降低[4]，生态调节功能逐步减弱，对西北生态环境造成了严重威胁[5]。青海云杉(Piceacrassifolia)作为祁连山森林的主要建群种[6],在维系西北地区生态平衡、水土保持、水源涵养和生物多样性等方面发挥着重要作用[7-8]。目前关于青海云杉的研究，大多集中在其生理生态特性及相关领域，虽然也有部分研究关注其分布特征，但未能进行较为准确的潜在分布预测。

潜在植被是指在没有人为干预的条件下，立地植被所能演替到与当地气候达到平衡态的顶级植被类型[9-10]。物种分布模型是利用物种分布数据以及环境数据来推断物种分布与环境之间的定量关系，构建模型，以物种实际分布及其分布区的环境特征来预测物种潜在分布情况。本研究以青海云杉为目标物种，利用最大熵模型(MaxEnt)、生态位因子分析模型(ENFA)、生物气候模型(Bioclim)及域模型(Domain)对目标物种进行潜在分布模拟。MaxEnt模型通过已知的物种分布特征及环境特征，建立约束条件，探索在该约束条件下熵最大的可能分布，以此确定物种潜在分布特征[11-13];ENFA模型通过比较物种分布区域环境与背景环境之间的差异来预测物种的潜在分布范围[14-16];Bioclim模型根据物种分布区域的环境特征构建矩形包络，以此来评价潜在适宜性[17-18];Domain模型通过比较环境空间的相似性，计算目标物种分布的适宜程度，确定其潜在分布[19-20]。随着时代进步，计算机技术及GIS的快速发展，环境数据的多样化，物种分布模型的种类不断增加且应用范围越发广泛。至今，物种分布模型已经被用于研究全球气候变化背景下的物种分布与气候环境之间的关系[21]、外来物种入侵预测[22]、物种的潜在分布预测和生态系统恢复及规划保护等诸多方面[23]。

对祁连山区域的青海云杉进行潜在分布研究，有助于了解青海云杉与环境之间的作用规律及影响机制，确定广域上青海云杉的分布可能，为区域生态恢复及规划提供重要依据，具有重要的科学与实践价值。

1 材料与方法

1.1 研究区概况

根据青海云杉主要分布区域，选择甘肃省祁连山国家自然保护区(97°23′～103°29′E，36°29′～39°42′N)作为研究区域，研究区地处甘肃省与青海省交界，主要包括甘肃省张掖市肃南裕固族自治县、民乐县、山丹县，以及武威市的天祝藏族自治县。总面积约为2.65×104km2，其中海拔2 500～4 500 m的区域面积约为2.31×104km2。祁连山国家自然保护区位于祁连山北坡，地势陡峭，海拔落差较大，最高的素珠链主峰海拔为5 564 m。研究区内降水主要集中在5—10月，年平均降水量429.2 mm，年均气温在0 ℃左右，最高气温出现在每年7月，月平均最高气温23.5 ℃，月平均最低气温-21.9 ℃。研究区气候具有明显的垂直变化特征，自下而上可划分4个垂直气候带，即浅山干草原气候带、中山森林草原气候带、亚高山灌丛草甸气候带和高山冰雪植被气候带[24]。主要土壤类型自上而下有高山寒漠土、高山草甸土、亚高山草甸土、高山草原土、灰褐土、黑钙土、栗钙土、灰钙土[25]。

1.2 数据材料

本研究所使用的数据主要分为4部分，即海拔数据、气温数据、降水数据及物种分布数据。研究区的数字高程模型(digital elevation model，DEM)来源于国家青藏高原数据中心(https://data.tpdc.ac.cn/)的祁连山地区30 m DEM数据集(2018)，数据垂直精度20 m、水平精度30 m。降水和气温数据来自国家青藏高原数据中心(https://data.tpdc.ac.cn/)的中国1 km分辨率逐月降水量数据集(1901—2017)，以及中国1 km分辨率逐月平均气温数据集(1901—2017)。物种分布数据根据中国植被分布图、遥感影像及野外调查数据，在Google Earth Enging平台划出准确分布区域，最终确定642处分布点。

1.3 数据处理

对影响青海云杉生理生态特征的环境要素进行分析[2,26-27]：考虑青海云杉分布的地形特征，选择海拔、坡度、坡向3个地形要素；考虑其分布区域的全年水热配置情况，选择1月平均气温(最冷月)、7月平均气温(最热月)、年均气温、1月平均降水(降水最少)、7月平均降水(降水最多)以及年降水6个气候要素，共遴选9个要素作为影响树种分布的环境因子。

对气候数据与地形要素数据进行处理，重采样至相同栅格大小(100 m×100 m)的环境要素图层。并随机选取70%的物种分布数据生成训练样本，将剩余30%的分布数据与分布点之外随机生成的“拟不存在点”构成检验样本。提取青海云杉分布区的环境特征，并通过标准差椭圆界定其环境阈值。

以受试者工作特性曲线(ROC)的下面积值(AUC)作为评价指标，对模型预测结果进行评价，AUC值为0～1.0。当AUC值大于0.9时，模型预测性能优异；当AUC值为0.7～0.9时，预测效果一般。

2 结果与分析

2.1 青海云杉分布区域环境特征分析

将青海云杉分布数据与环境数据叠置分析，提取青海云杉分布区域的环境要素，并对各要素数据处理分析，进行高斯正态分布拟合，确定其分布区域各环境要素的分布特征，且各要素拟合的决定系数(R2)均大于0.9，符合正态分布特征。对坡向数据进行适当调整，处理后的数据分布范围为 -180°～+180°，0°为正北方向，90°为正东方向。以均值左右2倍标准差为界(约95%的青海云杉分布其中)划分适宜区，各环境要素分布特征见表1。

表1 青海云杉分布区域环境要素统计特征Table 1 Statistical characteristics of environmental elements in the distribution area of Picea crassifolia

为界定多要素影响下的青海云杉分布环境阈值，将地形、气温、降水3方面的要素两两组合，地形方面以海拔为例，温度选择7月平均气温，降水因素选择7月平均降水，以正态分布曲率最大处为界(均值外)，构建标准差椭圆如图1所示。

图1 海拔-气温-降水要素影响下青海云杉分布的标准差椭圆Fig.1 The standard deviation ellipses of Picea crassifolia distribution under the influence of altitude-temperature-precipitation

由图1可以看出，标准差椭圆能有效表示青海云杉分布区域环境分布特征，在7月降水充沛时，随着海拔升高、气温降低的同时，青海云杉分布区的降水在逐渐增加。

海拔与7月平均气温的边界函数为：

式中：xdem为青海云杉分布区域的海拔，m；ytem_07表示其分布区域的7月平均气温，℃；t为参数，0≤t≤2π。

海拔与7月平均降水之间的边界函数为：

式中：xdem为青海云杉分布区域的海拔，m;ymonth_7表示其分布区域7月的平均降水量，mm;t为参数，0≤t≤2π。

7月平均降水量与7月平均气温的边界函数为：

式中：xmonth_7为青海云杉分布区域的7月平均降水量，ytem_07表示其分布区域的7月平均气温，t为参数，0≤t≤2π。

2.2 青海云杉潜在分布预测分析

利用Maxent、ENFA、Bioclim、Domain 4种物种分布模型对青海云杉潜在分布区进行预测，物种分布数据与预测结果叠置分析，以5%为阈值来划分是否适宜(即将95%的青海云杉样点划分到适宜区)，在适宜区内将适宜性等间距划分为低适宜区、中度适宜和高度适宜(如在MaxEnt模型中适宜区的适宜性范围为0.1～1.0，则将适宜性0.1～0.4 范围的划分为低适宜区，≥0.4～0.7划分为中适宜区，≥0.7～1.0划分为高适宜区)，确定青海云杉潜在适宜区分布，结果如图2所示。

图2 4种模型预测青海云杉适宜区面积统计图Fig.2 Statistical map of suitable areas of Picea crassifolia predicted by four models

4种预测模型中，ENFA模型预测的适宜区面积最大，其低适宜区及中度适宜区面积明显高于其他模型，结合预测图像可以发现ENFA模型的低适宜区呈带状分布，该模型是通过主成分分析方法，用少数变量代替全部变量进行模拟，此过程可能会忽略某种要素关键作用，如在现实中青海云杉大致分布在阴坡及半阴坡，而模型结果呈带状分布，这就与实际分布存在一定差距。Domain模型预测的各类适宜区面积均偏小，其预测的高适宜区仅40.7 km2，与其他3种模型的预测结果差距过大，这是由于Domain模型是通过计算环境空间相似性构成非连续的空间包络进行模拟，会存在一定的过拟合情况,这也是该模型与其他模型相比预测面积偏小的原因。Bioclim模型预测结果适宜区的面积偏大，这是由于其算法是根据影响要素构建矩形包络，而根据图1可知其要素两两组合后散点更接近椭圆而非矩形，通过矩形进行预测会过大估计其适宜区面积。MaxEnt、ENFA、Bioclim、Domain模型的AUC值分别为0.977、0.904、0.936、0.959，对比4种模型的初步预测结果并考虑其AUC值，MaxEnt模型的预测效果最佳，其AUC值在4种预测模型中最大，且预测结果与实际分布情况较为符合，故选取Maxent模型作为最优预测模型。

青海云杉潜在分布总面积约2 849 km2，研究区内的现有森林分布面积约1.702×103km2，其中与模型预测的潜在分布重合面积约848.8 km2，占总预测面积的30%，低、中、高3种适宜区内的现有森林分布占各适宜区的百分比分别为26%、42%、53%(图3)，这表明仍有许多地方是青海云杉潜力分布区。

图3 3种潜在适宜区内实际森林分布面积占比Fig.3 The proportion of actual forest distribution areas in the three potential suitable areas

随着人类活动加剧(如盗砍、采伐和放牧等)以及自然灾害频发(如干旱、雪压)等种种因素的影响，以青海云杉为优势树种的森林生态系统受到严重威胁[28]。1989—1998年，祁连山保护区内的裸地增加了20.7%，森林植被面积减少了33.5%[29]，这与祁连山内部牧民过度放牧有着密切联系。别强等[30]通过遥感技术等手段计算了祁连山国家自然保护区内的青海云杉的面积变化，得出1960—1990年祁连山国家自然保护区内青海云杉面积减少了11.1%。这些研究均表明祁连山国家自然保护区内仍有较大的森林恢复潜力。

MaxEnt模型预测结果显示，青海云杉潜在适宜区大多分布在祁连山国家自然保护区的中东部，越偏向西北，适宜区越少，适宜程度越低，祁连山中部黑河附近水分充沛，适宜区域明显增多。中度及高度适宜区大多分布在肃南裕固族自治县东部、民乐县南部、天祝藏族自治县北部，且河流附近适宜区分布较多。

2.3 环境变量贡献度分析

利用刀切法检测9种环境要素对模型的增益效果，结果见图4。考虑单要素和其他要素的训练增益以及各要素之间的相关性后，在9种环境要素中，对青海云杉分布贡献最大的3种要素分别为海拔、坡向和年降水量，其中海拔百分比贡献达到40.8%。综合考虑认为，地形因素对青海云杉分布影响最大，其次是降水因素，气温要素对其分布影响较小。从图4可见，坡向的增益贡献度比7月气温小，坡度似乎最小。要素对模型的增益贡献度除了单独考察该要素的训练增益，还需要考察除了该变量外其他所有变量对模型的增益效果，另外还需要考虑要素之间的相关性，其中海拔与气温、降水高度相关，综合考虑，计算结果为坡向的贡献度大于7月气温，坡度的贡献度大于其他气温、降水要素。

图4 青海云杉潜在分布环境变量重要性Jackknife检验Fig.4 Jackknife test of the importance of environmental variables in potential distributions of Picea crassifolia

3 结 论

1)通过对青海云杉分布区的环境要素进行提取，青海云杉主要分布在海拔2 414～3 753 m的阴坡及半阴坡，分布地区年降水量集中在298.26～540.16 mm，年均气温保持在-1.67～4.85 ℃。

2)通过4种生态位物种分布模型对青海云杉潜在分布区进行预测，MaxEnt模型预测结果最佳。结果显示，青海云杉潜在分布区总面积约2 849 km2，对比研究区内实际森林分布，与模型预测的潜在分布重合面积约848.8 km2，占总预测面积的30%，说明仍有较大的森林恢复潜力，且地形因素对青海云杉分布影响最大。

3)通过模型预测，适宜青海云杉生长的中度及高度适宜区主要分布在肃南裕固族自治县东部、民乐县南部、天祝藏族自治县北部，河流两岸适宜区较多，且适宜区从东南往西北逐渐减少。