基于BIOMOD的黄河源区高原鼠兔潜在分布及其影响因子

杜嘉星，孙 义，向 波，陈建军，秦 彧，侯秀敏，于红妍，宜树华

（1.冰冻圈科学国家重点实验室 / 中国科学院西北生态环境资源研究院，甘肃 兰州 730000；2.中国科学院大学，北京 100049；3.南通大学地理科学学院，江苏 南通 226007；4.南通大学脆弱生态环境研究所，江苏 南通 226007；5.重庆市气候中心，重庆 401147；6.桂林理工大学测绘地理信息学院，广西 桂林 541004；7.青海省草原总站，青海 西宁 810008）

高原鼠兔(Ochotona curzoniae)是分布在青藏高原及其周边地区特有的小型草食哺乳动物，是高寒草甸生态系统的关键种[1]。栖息地海拔为3 000-5 100 m，喜欢居住在坡位较低、土壤疏松干燥、地形开阔、距离水源地较近、植被低矮的高山草甸和亚高山草甸[2-4]。高原鼠兔对其所处生态系统的影响主要取决于其种群密度。当种群密度过高时，其啃食牧草、挖掘鼠洞，常形成生物灾害，能直接导致草地退化；而当种群密度在合理范围内，能作为青藏高原中肉食性动物的重要食物来源，其所挖掘的洞道为鸟类、两栖类动物提供了栖息场所，同时洞道的挖掘行为能增加入渗、减少径流，增加土壤水分和碳含量，对生态系统食物网结构及其能量流通和物质循环具有重要意义[5-11]。

认识高原鼠兔的空间分布及其影响因子对了解草地退化的原因、“黑土滩”的形成和高原鼠兔对栖息地的选择以及在生态系统中的作用有重要意义。传统的高原鼠兔调查方法有直接计数法、标志重捕法、指数估计法[12]。这些方法虽然准确度较高，但费时费力，难以展开大范围调查。无人机是利用无线电遥控设备和自备的程序控制装置操纵的，或者由车载计算机完全或间歇地自主操作的不载人飞机[13]。无人机应用前景广阔，如无人侦察机和靶机，以及航拍、农业植保、测绘等无人机[14-15]。在生物群落调查中，无人机具有成本低、效率高等特点，是进行生物群落快速调查的优良工具。例如，Watts等[16]用无人机影像来估计美国短吻鳄(Alligator mississippiensis)的种群数量，Israel[17]用无人机装载热红外相机来监测西方狍(Capreolus capreolus)的活动，郭新磊等[18]在青藏高原地区用无人机获取高原鼠兔的洞口数量。

物种分布模型 (species distribution models, SDMs)是利用物种的观测数据与环境数据基于统计或理论推导得出的经验模型[19]。物种数据可以是基于随机或分层野外采样的简单存在、存在/缺失或物种丰富度观测数据[20]。环境预测因子可以直接或间接影响物种，按照影响程度的大小排序，可以分为限制因子、扰动和资源[21-22]。目前，物种分布模型已经成为基础生态学和生物地理学研究的重要工具，全球变化背景下物种的分布和气候之间的关系、区域气候变化对植物群落和功能的影响、生态系统功能群和关键种的监测和预测、生态系统不同尺度多样性的管理和保护、外来物种入侵区域的预测、生态系统的关键物种的潜在分布预测和保护区规划等方面的研究上均有广泛应用[23]。BIOMOD(BIOdiversity MODelling)是一个基于R语言的用于物种分布集成预测的计算机平台，能够处理一系列的模型中的方法不确定性和物种-环境关系的检验，由Wilfried Thuiller博士于2003年提出[24-25]。它可以使用多种不同方式建模，并对模型进行评估，从而使模型精度最大化并预测物种分布。目前该模型的最新版本为2016年发布的BIOMOD2(Version:3.3-7)。本研究应用BIOMOD物种分布集成模型集成平台，结合无人机遥感技术(unmanned aerial vehicle remote sensing, UAVRS)，对高原鼠兔在黄河源区的潜在分布进行研究，揭示其空间分布的规律并探究影响其空间分布的限制因子。

1 材料与方法

1.1 研究区概况

黄河源区位于青藏高原东北部边缘(95°50′-103°30′ E, 32°30′-36°10′ N)，横跨青海、甘肃、四川三省，平均海拔 4 100 m，年平均气温-2.7 ℃，年平均降水517 mm。流域范围达12.2万km2，是黄河流域的重要产流区。研究区植被类型以高寒草甸和高山草原化草甸为主，优势种有矮生嵩草(Kobresia humilis)、垂穗披碱草(Elymus nutans)、青藏苔草(Carex moorcroftii)、异针茅(Stipa aliena)、高原早熟禾(Poa alpigena)等[26-27]。本研究区界定在河源至唐乃亥水文站的黄河流域范围。

1.2 数据获取

1.2.1 高原鼠兔分布数据

本研究所采用的高原鼠兔分布数据来自于2018年7-8月进行的野外调查，在研究区范围内共布设了208个工作点，使用大疆创新科技有限公司 (DJ-Innovations, DJI)生产的 Phantom 3 Professional和Mavic Pro无人机进行航拍。由无人机生态环境研究网络团队自主研发的FragMAP软件设置工作点及航线，控制无人机按照设定航线飞行并定点拍照[28]，航线的中心点设置与MODIS像元中心点相同，飞行高度为20 m，覆盖地面范围约为910 m2(35 m × 26 m)，单张照片分辨率为 4 000 像素 × 3 000像素。后期使用自主研发的Proposal Classifier软件对照片进行处理，选取合适的阈值，通过自动识别和人工手动矫正相结合的方法对高原鼠兔洞口进行提取，并将处于同一航线上的洞口数量进行平均，结合当地的有效洞口系数，获取高原鼠兔的存在/不存在点。

1.2.2 环境数据

本研究共获取了48种环境数据，包括以下几种类型。

1)气候数据：气候数据来自中国科学院青藏高原研究所阳坤团队[29-30]，包括年平均气温、气温周较差、等温性、温度季节性、最高气温、最低气温、气温年较差、最湿季平均气温、最干季平均气温、最暖季平均气温、最冷季平均气温、年降水量、最湿期降水量、最干期降水量、降水季节性、最湿季降水量、最干季降水量、最暖季降水量、最冷季降水量、年平均辐射、最大辐射、最小辐射、辐射季节性、最湿季辐射、最干季辐射、最暖季辐射、最冷季辐射、年平均湿度、最大湿度、湿度季节性、最湿季湿度、最干季湿度和最暖季湿度的0.1°格网数据，采用克里金插值的方法扩展到整个研究区，空间分辨率1 km，共33个栅格图层。

2)植被数据：植被数据来自LAADS DAAC(ladsweb.modaps.eosdis.nasa.gov)的 MODIS/Terra MOD1-3A3植被指数产品，空间分辨率1 km，时间分辨率30 d。时间范围为2008年11月-2018年10月(共11 年)，通过使用 MRT(Modis projection tool)软件进行读取并进行投影，在QGIS Desktop中进行镶嵌和裁剪，最后通过像元统计得到NDVI年平均值、NDVI最大值、NDVI最小值、NDVI年较差共4个栅格图层。

3)地形数据：数字高程数据(digital elevation model, DEM)为 SRTM 90 m DEM 数据，来自美国地质勘探局(www.usgs.gov)。使用QGIS Desktop对DEM数据提取坡度和坡向，得到海拔、坡度、坡向共3个栅格图层。

4)土壤数据：土壤数据来自SoilGrids(www.soilgrids.org)[31]1 km分辨率的土层厚度、0.3-0.6 m土层有机碳储量、0.3 m处土壤容重、0.3 m处黏土含量、0.3 m处粗屑体容积、0.3 m处淤泥含量、0.3 m处含沙量、0.3 m处土壤pH数据，共8个栅格图层。

在使用上述环境数据构建模型过程中，为了减少因子自相关对模型精度的影响，对其进行Pearson 相关性分析，保留|r| ＞ 0.7 的环境变量对中的其中一个。经过筛选后剩下年平均气温、气温周较差、等温性、温度季节性、年降水量、最湿期降水量、最干期降水量、年平均辐射、辐射季节性、海拔、坡度、坡向、NDVI年平均、NDVI年较差、有机碳储量、土壤容重、土壤pH共17种环境数据用于建模。

1.3 模型介绍

在BIOMOD中共提供了10种物种分布模型可供选择，并且可以通过不同指标计算模型精度来评估不同模型在研究区的适用性，从而筛选出最优的模型，包括广义线性模型(generalized linear model,GLM)[32]、广义增强回归模型 (generalized boosted regression models, GBM)[33-34]、广义加性模型 (generalized additive model, GAM)[35]、 分 类 树 分 析 (classification tree analysis, CTA)[36]、人工神经网络 (artificial neural network, ANN)[37]、表面分布区分室模型 (surface range envelope, SRE)[38]、弹性判别分析 (flexible discriminant analysis, FDA)[39]、多元自适应回归样条 (multiple adaptive regression splines, MARS)[40-41]、随机森林 (random forest, RF)[42]和最大熵模型 (maximum entropy model,MaxEnt)[43]。

1.4 模型检验

为了评价模型的预测精度，将样本划分为两个子集，一个子集包含全部样本的70%，另一个包含30%。样本的70%用于建模，30%用于检验，并引入真实技巧统计值 (true skill statistics, TSS)和 AUC(area under curve)两种指标来对模型进行评估。

1.4.1 TSS

TSS = Sensitivity + Specificity - 1 = TPR - FPR。其中，Sensitivity为敏感性，Specificity为特异性，TPR(true positive rate)为真阳性率，FPR(false positive rate)为假阳性率[44]。TSS 的取值范围在 (0, 1)区间，TSS的值越接近于1，表示真阳性率与假阳性率的差值越大，模型的效果就越好。

1.4.2 AUC

AUC来源于受试者工作特征曲线(receiver operating characteristic curve, ROC)，被定义为 ROC 曲线下与坐标轴围成的面积[45]。ROC曲线是反映敏感性和特异性连续变量的综合指标，ROC曲线上每个点反映着对同一信号刺激的感受性。AUC的取值范围一般在(0.5, 1)区间内。越接近于1的模型表示预测效果越好，而越接近于0.5的模型表示越接近于随机猜测，模型没有预测价值。因此，AUC值越大的模型效果越好。

2 结果与分析

2.1 高原鼠兔分布

在实地调查中发现，如选用是否存在高原鼠兔洞口作为判别高原鼠兔是否存在的依据，则完全没有洞口的样本量非常少，故在本研究中选用是否存在1个或更多有效洞口作为判别依据。根据Yi等[46]在青藏高原地区不同草地类型的调查，该地区的有效洞口系数在0.22～0.42。根据泊松分布，10个洞口中存在1个或更多有效洞口的概率P＞ 0.96，可以认为10个洞口中至少存在1个活动的洞口。故对于洞口个数大于10的样地认为有高原鼠兔存在，小于10则认为不存在。通过对黄河源区野外获取的航拍照片进行洞口提取，共获得了82个存在点数据和126个不存在点数据，并根据航拍工作点的定位进行绘图(图1)。

图1 黄河源区采样点Figure 1 Sampling points in the source region of the Yellow River Basin

2.2 模型精度评估

将10种模型每个重复运行多次，每次随机选取70%的样本进行建模，30%的样本对建模的结果进行评估。10种模型中除GAM由于获取估计参数困难运行失败外，其余9种模型均成功运行，将9种模型运行多次得到的TSS与AUC两个参数进行描述性统计(表1)，并绘制箱形图(图2)。在9种模型中，RF的整体表现最佳，TSS平均值为0.57，AUC平均值为0.83，变异系数仅为11%和5%，说明在多次运行的过程中模型能较准确地预测并能保持稳定性。次优的模型为GBM、MARS和FDA，这3种模型均能一定程度上预测高原鼠兔的分布。其余模型均表现一般，其中SRE最差，TSS为0.09，AUC 为 0.54，接近于随机分类 (AUC = 0.5)，可以认为该模型基本无法在该地区准确预测。

2.3 环境因子重要性

使用模型内置的因子重要性计算函数，将包含全部环境因子的集合定义为ref数据集，而随机剔除单个环境因子后的集合定义为shuffled数据集，使用两个集合进行预测，并计算预测结果的简单相关性(默认为Pearson相关性)。计算公式：Importance = 1 - cor(pred_ref, pred_shuffled)。剔除掉表现较差的SRE以及运行失败的GAM，统计其余8种模型的各个环境因子的重要性求平均并按大小排序(表2)。

表1 9种模型TSS与AUC描述性统计Table 1 Descriptive statistics of the TSS and AUC of nine models

图2 9种模型TSS与AUC比较Figure 2 Comparison of the TSS and AUC of nine models

表2 8种模型环境因子重要性Table 2 Importance of the environment factors in eight models

2.4 预测结果

根据以上对模型精度的评估结果，采用准确度和稳定性均较好的RF。在多次重复运行结果中选取预测效果最佳的一次结果 (TSS = 0.74, AUC = 0.92)，对结果进行绘图(图3)。

预测结果显示，高原鼠兔分布概率较高的区域主要有玛多县扎陵湖乡、玛沁县优云乡、同德县尕巴松多镇、河南县优干宁镇以及达日县与甘德县交界处的上贡麻乡，而在玛曲县东南部、若尔盖县西部与红原县北部三县之间的山区、玛多县玛查里镇东南部以及同德县和泽库县北部则分布概率较小。

3 讨论与结论

本研究使用无人机于2018年7-8月在黄河源区进行野外调查，获取了3 000余张航拍照片，从中提取了208个高原鼠兔的存在/不存在点数据，从气候、植被、地形和土壤4类共48种环境因子筛选出17种作为解释变量，使用BIOMOD集成模型平台对该区域高原鼠兔的潜在分布进行预测。本研究方法相比较传统的地面调查效率更高且更适合大范围调查，预测的结果能为鼠害防治工作提供科学依据，同时能为今后的高原鼠兔研究提供新思路。

BIOMOD中集成的10种物种分布模型中表现最好的是RF，多次运行的TSS平均值为0.57，AUC平均值为0.83，变异系数分别为11%和5%，说明该模型在黄河源区能较为准确地预测高原鼠兔的分布，并且在多次运行中仍能保持较好的稳定性。BIOMOD在此区域的应用能降低预测的不确定性和误差，提高预测的精度。

高原鼠兔的分布受多种环境因子影响，重要性较高的几个因子有海拔、NDVI年平均、坡度、土壤pH、年降水量、坡向以及年平均辐射。高原鼠兔的分布与海拔高度相关，其适合生长在高海拔低氧地区，王淯等[4]的研究表明，高原鼠兔多分布在海拔 3 000-5 100 m、地势较为平坦且坡位较低的区域。高原鼠兔的分布同时受坡度和坡向的影响，卫万荣[47]的研究表明，84.2%高原鼠兔喜好将栖息地选择在坡度在0°～5°的区域，而在6°～10°和11°～15°的区域分布分别只占13.2%和2.6%；对坡向的选择更倾向于南或偏南向；在山地地貌上占总数的100%，在所有地貌上占41.3%。这有利于增加冬季对太阳辐射的获取，提升洞内温度。NDVI能够在一定程度上反映植被的盖度和生物量，马波等[48]的研究认为高原鼠兔洞口数与NDVI存在显著正相关的线性关系，高原鼠兔有选择地利用NDVI较高、植被较好的生境。高原鼠兔的挖掘活动能够改变当地的土壤状况，Guo等[49]的研究认为高原鼠兔的洞穴对土壤pH、含水量、含沙量等土壤理化性质均有影响，能增加土壤渗透率，并提高土壤养分。同时高原鼠兔的分布受降水量和辐射的影响，郭新磊[50]的研究表明，高原鼠兔洞口密度随降水量和辐射的增高表现为先增大后减少，降水量在600～700 mm区间和辐射在240～250 J·m-2区间的地区高原鼠兔洞口密度达到最大。以上影响高原鼠兔空间分布的因子与前人的研究相一致，说明BIOMOD不仅能很好地的模拟高原鼠兔的空间分布，在影响因子的确定上也有着良好的效果。

图3 基于RF的黄河源区高原鼠兔分布预测Figure 3 Prediction of the distribution of plateau pika, based on RF in the source region of the Yellow River Basin

高原鼠兔分布除了与环境因子有关外，还与人类活动密切相关，如放牧、害鼠防治、招鹰架的布置等。张兴禄[51]在玛曲高寒草甸植被区所做试验证明轮牧制度可有效抑制高原鼠兔数量增长，而在连续放牧制度下，低放牧率下高原鼠兔的数量则明显上升，而高放牧率下由于牛羊践踏、采食过重，不适合高原鼠兔生存，数量则有所下降。梁杰荣等[52]在海北高寒草甸生态系统定位研究站调查了高原鼠兔的种群恢复过程，发现害鼠防治后残鼠数量按Logistics曲线增长；当防治效果高达90%时，其数量约在两年内也恢复至原水平。侯秀敏[53]在泽库巴滩所做试验表明，招鹰架的布置能有效吸引高原鼠兔的天敌，降低高原鼠兔的密度。当鹰架布置的间距为 500 m × 500 m 时，控制面积为每25 hm2一根，鹰架的利用率达到了98%，样地内两年间高原鼠兔的有效洞口数量降低了26.47%。

致谢：本研究使用的气候数据由青藏高原多圈层模拟与数据同化中心 (Data Assimilation and Modeling Center for Tibetan Multi-spheres)开发，中国科学院青藏高原研究所阳坤研究员提供，特在此表示诚挚的谢意。