基于Cubist 模型的天山北坡草地鼠群密度时空分布特征

潘 群， 施海洋， 张文强， 罗格平,4， 陈春波

（1.新疆大学资源与环境科学学院，新疆 乌鲁木齐 830046；2.中国科学院新疆生态与地理研究所荒漠与绿洲生态国家重点实验室，新疆 乌鲁木齐 830011；3.中国科学院大学，北京 100049；4.中国科学院中亚生态与环境研究中心，新疆 乌鲁木齐 830011）

啮齿动物是草地哺乳动物中种类和数量最多的一类动物，就物质循环和能量流动而言，其在草地生态系统中具有重要的功能和地位［1］。然而，其种群数量过高可能会对草地造成危害即草原鼠害，其危害主要表现在啃食和储藏大量牧草、挖土造丘，导致草地裸露，严重破坏植被空间结构，降低草原生产力和草地牧草产量［2-3］，打破原有“土-草-畜-鼠”的动态平衡，致使草地退化严重［4］。

新疆具有丰富的草地资源，然而近年来，由于过度放牧、草地资源开发不当、全球气候异常等，导致新疆草地鼠害持续大面积发生，已对畜牧业的可持续发展和草原生物多样性产生了不利影响。为了更及时地防治鼠害，需要加深对啮齿动物种群密度时空分布驱动机制的理解并更准确地预测其时空分布。目前，虽然已有较多对于啮齿动物空间分布格局的研究，但大多集中于研究小尺度的种群数量随时间的变化［5］或研究其潜在空间分布，却很少涉及大尺度的种群密度时空分布的变化特征及其机制。如乔雪丽等［6］利用最大熵模型（MaxEnt）对在不同环境条件下的大沙鼠在亚洲的潜在地理分布进行预测；何咏琪等［7］基于“3S”技术与地面实地调查数据，建立了草原鼠害监测模型并模拟出了青海省草原鼠害分布区；Mohammadi 等［8］模拟了伊朗气候变化下两种沙漠跳鼠的当前和未来潜在分布。

目前，环境因子已被广泛用于物种种群密度分布的驱动机制分析和时空分布预测［9-10］。对啮齿动物而言，环境条件在其种群密度波动和空间分布中扮演着重要角色，可通过直接或其相互作用对其产生重要影响［11］。啮齿动物种群密度的时空分布受地形［12］、土壤质地［13］、气温［5］、降水［14］、植被［15-16］、放牧［17］等的影响，且其具有年际和季节性变化的特征［18］。因此，通过建立啮齿动物种群密度与环境因子的关系来预测其种群密度的时空分布具有一定可行性。在衡量害鼠种群密度的指标中，有效洞口密度是最有效的指标之一［19］。该指标已广泛地在实地调查中被估计并记录，样本数相对丰富，适用于本研究中大尺度种群密度的归一性分析和预测。

天山山脉北坡属于中亚干旱地区典型的山地-绿洲-沙漠（从高地到低地）生态系统［20］，动植物的分布很大程度上取决于海拔的变化。在天山北坡，啮齿动物对草原的损害主要由分布在平原荒漠草地的大沙鼠（Rhombomys opimus）和山地草原的黄兔尾鼠（Eolagurus luteus）引起，且大沙鼠主要栖息在海拔900 m 以下的荒漠草地［21］。因此，本文结合文献中数据特点将天山北坡分为海拔≤900 m 和＞900 m 2 类区域，且主要的优势种分别为大沙鼠和黄兔尾鼠。

综上，本文以有效洞口密度反映害鼠的种群密度，从已发表的文献收集天山北坡草地小尺度的有效洞口密度的实地调查信息，利用遥感技术、统计分析等手段，对2类区域分别运用随机森林（RF）模型和Cubist 模型模拟有效洞口密度。最后，选取较优方法定量研究环境因子对有效洞口密度时空分布的影响，对1982—2015年的天山北坡草地有效洞口密度时空分布进行预测，并分析其变化特征。

1 数据与方法

1.1 研究区概况

天山北坡（图1）是指由横贯新疆中部天山山脉发育的河流向北流经的区域，地理位置为79°53′～96°06′E，42°50′～46°12′N，东起伊吾，西至温泉县，北连古尔班通古特沙漠。区域内自然条件的垂直分带性明显，地势南高北低，海拔由山地的5000 m左右下降到沙漠边缘的200 m左右。由于位于亚欧大陆腹地，远离海洋，属于典型的温带大陆性气候，其特点是冬季寒冷漫长，夏季干旱炎热，年均温在-6.5～11.0 ℃，降水季节变化显著，年内降水量主要集中在5、6 月，2 月最少，年降水量由南部山区中低山带的500 mm降至北部沙漠区的100 mm左右。天山北坡独特的自然条件使其成为全球温带干旱区大型山地-绿洲-沙漠生态系统的典型代表［22］。

1.2 数据来源

本研究采用的数据主要包括有效洞口密度和环境因子数据（表1）。有效洞口密度数据通过已公开发表的文献［23-45］获取，所搜集记录中仅有地名但无经纬度信息的通过GPSSPG 网站（http://www.gpsspg.com/）进行查询［46］。对得到的数据进行整理，最终得到了275条具有时空信息的有效洞口密度数据，年份跨度为1982—2015年。本文初步选择了21个可能影响有效洞口密度的环境因子，主要包括年份、月份、地形、土壤质地、植被、气温、降水量、短波辐射和放牧强度。

表1 环境因子数据来源Tab.1 Data sources of environmental factors

1.3 研究方法

首先，通过ArcGIS 10.2 软件对获得的1982—2015年的有效洞口密度时空信息进行整理，生成含有准确空间位置信息的月尺度矢量文件。随后，根据每个样本的时空信息提取对应的21 个环境因子并根据海拔将天山北坡分为≤900 m 和＞900 m 2类。最后，分别基于RF模型和Cubist模型拟合有效洞口密度，并选取较优模型定量研究有效洞口密度时空分布的变化特征及其机制。

1.3.1RF模型RF 模型是一种以决策树为基本单元的集成学习方法［56］。其基于Bagging算法，从原始样本集N中有放回地随机抽取n个样本生成新的训练样本集，通过算法训练得到n棵决策树，最后通过平均的方式整合各决策树的预测结果，便可得到RF模型的最终结果。RF 模型在每个回归树的建立过程中，对样本空间和特征空间进行随机采样，引入随机属性，进而降低了回归树模型之间的相关性，又通过组合大量的回归树，改善了模型的泛化性，使得此算法具备高效、准确等特点。本研究利用网格搜索算法［57］探索各区数据的最优参数组合，采用Python中scikit-learn包进行模型构建［58］。

1.3.2Cubist模型Cubist 模型是Quinlan 的M5′模型树的进一步发展［59］。Cubist模型在多元线性模型中得到一系列if-then 形式的规则，当一组变量符合规则的条件时，使用相应的模型计算变量的预测值［60］。Cubist算法用一系列组合起来的分段线性模型，很好地解决了非线性问题。此外，该方法还具有变量的重要性分析功能［61］。本研究使用R软件中Cubist 程序包，其主要参数为committees 和neighbors。

1.3.3 模型精度验证本文从数据集随机选择80%作为模型的训练样本集，剩余的20%作为模型验证样本集。选用均方根误差（RMSE）和决定系数（R2）作为模型精度评价指标，综合评价RF模型和Cubist模型的预测能力。计算公式如下：

式中：Yi、Yi′和分别为有效洞口密度实际值、预测值和样本平均值；n为样本数。

2 结果与分析

2.1 有效洞口密度样本数据分布特征分析

对从文献中获取的有效洞口密度数据进行统计分析，如表2所示，共275条数据，其中最小值为1个·hm-2，最大值为1370 个·hm-2，平均值为254.23个·hm-2。海拔≤900 m地区的样点有效洞口密度范围为4～1220 个·hm-2，海拔＞900 m 地区的样点有效洞口密度在1～1370个·hm-2之间。

表2 有效洞口密度统计描述Tab.2 Statistical description of active burrow entrance densities

2.2 草原有效洞口密度拟合模型精度评价与优选

将全部环境因子作为自变量，有效洞口密度为因变量，分别构建RF模型和Cubist模型并计算R2和RMSE（图2）。结果表明，在海拔≤900 m 地区其RF模型和Cubist 模型的R2分别为0.773 和0.832，在海拔＞900 m 地区其RF 模型和Cubist 模型的R2均在0.600以上，Cubist模型均高于RF模型。同时，RF模型的RMSE均大于Cubist模型。综上，Cubist模型具有更高的拟合精度和更小的误差，因此本研究采用Cubist模型来预测天山北坡草地有效洞口密度。

图2 有效洞口密度拟合值在海拔≤900 m和>900 m地区的验证散点图Fig.2 Scatter polt of fitting values of active burrow entrance densities in areas with an elevation less than or equal to 900 m and higher than 900 m

2.3 影响因子的重要性分析

本文基于Cubist 模型各因子的重要性排序，分析环境因子对天山北坡草地有效洞口密度时空分布的影响（图3）。Cubist 模型可根据样本信息将相关性不强的输入因子剔除，结合预测因子的相对重要性，筛选出预测结果较好的影响因子。重要性显示海拔≤900 m 和海拔＞900 m 地区之间存在着显著差异。在天山北坡海拔≤900 m地区，对草地有效洞口密度分布贡献最大的是月最低气温，其次为月平均气温、叶面积指数（LAI）。而对于＞900 m地区，贡献最大的为年累积降水量，其次为放牧强度、归一化植被指数（NDVI）。气候因子、植被状况与人类活动因子对天山北坡有效洞口密度时空分布影响较大。

图3 有效洞口密度的影响因子重要性Fig.3 Importance of influencing factors for active burrow entrance densities

2.4 天山北坡草地有效洞口密度时空分布

使用此Cubist模型对研究区草地有效洞口密度时空分布进行反演（图4）。并根据Cubist 模型分别计算1982—2015年天山北坡草地有效洞口密度在4月和10 月的平均值（图5）。比较1982、1993、2004年和2015 年有效洞口密度的年际空间分布差异以及其春季（4月）和秋季（10月）的季节性变化。结果表明，有效洞口密度具有明显的区域差异，靠近人类活动频繁的农田生态系统周围草地有效洞口密度普遍较高，而北部受人类干扰较少的荒漠地区其密度较低。在海拔≤900 m 地区4 月有效洞口密度平均为354个·hm-2，10月平均为424个·hm-2；海拔＞900 m地区4月有效洞口密度平均为421个·hm-2，10月平均为407 个·hm-2。1982—2015 年海拔≤900 m地区有效洞口密度整体为波动递增的趋势，34 a间4月的有效洞口密度增加了37.25%，10 月增加了41.58%；而海拔＞900 m 地区整体为波动递减的趋势，34 a 间4 月减少了4.83%，10 月减少了26.99%，且草地有效洞口密度波动剧烈（图5）。在天山北坡海拔≤900 m地区每年10月的有效洞口密度相对于4月显著增加，体现出较强的季节差异性。

图4 1982—2015年天山北坡草地每11 a的有效洞口密度空间分布Fig.4 Spatial distributions of active burrow entrance densities in grassland on the north slope of Tianshan Mountains every 11 years from 1982 to 2015

图5 1982—2015年天山北坡草地有效洞口密度在海拔≤900 m和>900 m地区的变化趋势Fig.5 Change trend of active burrow entrance densities in grassland on the north slope of Tianshan Mountains in areas with an elevation less than or equal to 900 m and higher than 900 m from 1982 to 2015

3 讨论

本研究主要基于Cubist 模型、结合环境因子预测有效洞口密度，基本实现了草原小型啮齿动物种群密度时空量化与表达，有效地弥补了目前在大时

空尺度草原害鼠时空变化规律研究方面的不足。已有的研究因缺乏长时间序列以及较大研究区的观测，导致对小型啮齿动物种群密度区域分布差异的认识存在局限性。

在研究方法方面，已有的研究中鲜有将RF、Cubist模型等机器学习方法应用到小型啮齿动物种群密度估算的研究方面。因此，本研究尝试将不同机器学习方法在该领域进行应用，这也是本文的主要创新点之一。所构建的Cubist模型在量化环境因子对有效洞口密度的影响方面，优于RF 模型。草原小型啮齿动物生态特点是种群数量波动剧烈，有些年份较低，有些年份猛增，且其分布具有明显的区域性差异［62］。Cubist模型可以针对不同区域构建合适的应用模型，其能够同时处理连续的和离散的数据［63］，因此，相对于RF模型，Cubist模型可能能够更好地模拟草原有效洞口密度。

首先，本研究发现植被状况在鼠类种群密度的时空分布中占重要地位，但一个有趣的现象是，在海拔≤900 m地区，影响有效洞口密度时空分布的主要是LAI，而对于海拔＞900 m 地区为NDVI（图3）。这可能是由于遥感NDVI数据在监测梭梭的生长和地下生物量方面不是很有效。因为在海拔≤900 m的平原荒漠地区，大沙鼠将梭梭的嫩根和树枝作为主要食物，而黄兔尾鼠主要食用青蒿属植物。有效洞口密度与LAI 或NDVI 之间关系的差异可能受到海拔≤900 m 和＞900 m 地区中大沙鼠和黄兔尾鼠所消耗不同类型植被的影响。植被生长状况直接决定鼠类食物数量和品质［64］。植物群落的各项指标（植被高度、盖度、密度和地下生物量以及群落多样性等）都综合地影响着啮齿动物种群密度［65-67］。因本文中有效洞口密度数据是通过已公开发表的文献获取，而文献中缺乏调查有效洞口密度时所涉及的植被基本性状指标，且我们目前无法获取天山北坡植被群落的空间数据。因此，本文环境变量中的植被性状指标只涉及植被指数（NDVI、LAI），虽植被指数已广泛用来定性和定量评价植被覆盖及其生长活力［68］，但对于优势植物种类等性状指标的缺乏，可能会使环境因子的预测能力受到很大的限制。此外，使用NDVI和LAI数据较粗糙的空间分辨率可能会给此类分析带来不确定性，特别是在植被覆盖率较低的沙漠地区，这可能导致低估空间植被变异性对有效洞口密度的影响。除食物因素外，气候因子与放牧强度对啮齿动物种群密度也产生了极大的影响，并且影响会因地区而异。对比2 类地区，发现有效洞口密度对温度响应显著的区域为海拔≤900 m 地区，而降水影响显著的区域为＞900 m地区。这可能是因为在＞900 m 海拔带降水的时空异质性较高，因此对植被生长的影响更大［69］，这与天山北坡和全国范围的相关性空间研究结论是相吻合的［70-71］。此外，梭梭是大沙鼠的主要食物来源，而黄兔尾鼠主要食用青蒿属植物，且黄兔尾鼠其在阴雨刮风等气温较低天气条件下，活动性明显降低。这可能也是黄兔尾鼠种群密度对降水更敏感的原因。在自然情况下，鼠类通过与草地环境的相互作用、相互协调使之经常保持动态的平衡状态。然而，长期以来，人们为经济利益所驱动的高强度利用方式，深刻地影响鼠类及天敌动物赖以生存地草地基质条件，而群居的植食性鼠类生殖力强，种群数量增长快，往往在此类次生环境中占据优势而发生鼠害［72］。一些学者基于大面积的草地调查资料也发现了草地鼠类种群密度的空间分布与放牧强度有密切关系［73-74］。

其次，天山北坡海拔＞900 m靠近人类活动频繁的农田生态系统周围的草地区域，其有效洞口密度普遍较高，这可能是因为黄兔尾鼠主要栖息在丘陵及荒漠草原，农垦后聚集在农田及水渠附近［75］。黄兔尾鼠数量年际变动较为剧烈，往往在中低密度水平连续维持若干年之后，种群密度突然升高，而后又处于低潮期。低潮后，在自然条件适宜的年份密度又开始增大（图5）。由本研究可以发现，黄兔尾鼠数量的年际变动可能主要归因于降水量、放牧强度和植被状况，但鉴于黄兔尾鼠种群密度与植被-降水-人类活动的关系复杂、时空不确定性大，这一现象的产生机制，还需要进一步探讨。

最后，环境因子和有效洞口密度是一个动态变化的过程，遥感技术能够及时进行大范围监测了解环境各个时期的变化。有效洞口密度的时空动态在响应环境因子和人类活动的复杂作用时展现出叠加效应，需要从多视角研究小型啮齿动物种群密度时空变化。本文选择了1982—2015 年的数据进行建模，能够反映啮齿动物种群密度的年际和季节的变化情况，但对于偶尔极高的繁殖率，环境因子的预测能力可能受到很大的限制。此外，鼠类种群年龄结构、天敌种群分布以及人为因素［76-77］在鼠洞密度分布的预测中具有重要影响，这在未来的研究中需要更多关注。

4 结论

本研究收集了1982—2015 年的天山北坡有效洞口密度的小尺度时空分布数据和环境因子数据，并基于Cubist 模型估算天山北坡草地有效洞口密度，量化了环境要素对小型啮齿动物种群密度的影响，取得了较好的研究精度和效果，为区域草原鼠害预防预测工作提供了新思路、新方法。得到以下结论：

（1）利用Cubist 模型和RF 模型对天山北坡草地有效洞口密度进行预测，海拔≤900 m地区拟合效果R2分别为0.832和0.773，海拔＞900 m地区R2分别为0.667 和0.612，Cubist 模型拟合精度明显优于RF模型。

（2）基于Cubist 模型的影响因子重要性分析结果，植被状况、气象因子和放牧强度是天山北坡有效洞口密度时空分布的主要环境驱动因素。在天山北坡内的海拔≤900 m 和＞900 m 地区中，有效洞口密度的驱动机制存在着显著差异。

（3）在海拔≤900 m地区，气温对有效洞口密度的贡献度比在海拔＞900 m地区更高，而降水量对其影响更低。在海拔≤900 m地区，影响有效洞口密度时空分布的主要是LAI，而对于海拔＞900 m 地区为NDVI。这可能是受到大沙鼠和黄兔尾鼠所消耗不同类型植被的影响。

（4）有效洞口密度的时空分布具有明显的年份和季节波动，1982—2015 年海拔≤900 m 地区有效洞口密度整体为波动递增的趋势，34 a 间4 月有效洞口密度增加了37.25%，10月增加了41.58%；而海拔＞900 m 地区整体为波动递减的趋势，34 a 间4 月减少了4.83%，10月减少了26.99%。

（5）在海拔＞900 m 地区，草地有效洞口密度波动剧烈，有效洞口高密度区主要分布在靠近人类活动的牧区和绿洲区。这可能主要归因于降水量、放牧强度和植被状况的叠加效应，但鉴于黄兔尾鼠种群密度与植被-降水-人类活动的关系复杂、时空不确定性大，这一现象的产生机制，还需要进一步探讨。未来可通过大数据分析手段、提高数据的精度或多视角考虑种群密度的影响因子，尤其是植被基本性状指标，并与空间异质性探测方法结合，多层次深化本研究方向和主题。