1 引 言
畜牧业直接影响了全球近10亿人的生计和粮食安全［1-2］，近几十年来，全球财富积累，人们对动物源性食品的需求增加，促进了畜牧业的发展，但同时给自然资源和环境带来了巨大的压力。不断增加的资源压力迫使人类和家畜向野生动物栖息地靠近，进而增加了由野生动物携带的传染性人畜共患病病原体在家畜和人体内出现和传播的机会［3-5］。要解决上述问题，需要获取较为精细的家畜空间分布数据。而传统的家畜统计数据多以行政区划为单元进行收集，空间分辨率低，难以提供家畜空间分布的细节信息，且无法与基于网格等基础地理单元的数据进行整合计算。因此，家畜统计数据网格化，即将以行政区为基本单元的家畜统计数据扩展到一定空间尺度的网格上，实现社会统计数据与自然数据的整合，具有重要的理论和现实意义［6］。
目前国外对家畜网格化研究方面做了很多有益探索，生产了一定的家畜网格数据。联合国粮食与农业组织（Food and Agriclture Organization， FAO）于2007年发布了世界家畜网格数据库（Gridded Livestock of the World， GLW），它使用多层线性回归方法，提供了第一份标准化的全球家畜密度分布图，空间分辨率为3′（赤道上约为5 km），数据时间为2005年［7］。2014年，Robinson等［6］对GLW进行改进，包括自动化处理、改进输入数据，得到了1 km分辨率的牛、猪、鸡的全球分布图和鸭的部分分布图（GLW2），数据时间为2006年。2016年，Nicolas等［8］分别使用随机森林和多层线性回归方法对非洲牛和亚洲鸡的统计数据进行了空间降尺度，分别得到了1 km分辨率和10 km分辨率的非洲牛和亚洲鸡的分布图，并对两种方法的表现进行评估，证明基于随机森林回归方法的预测结果比多层线性回归方法更好。2018年，Gilbert等［9］对前两个版本的世界家畜网格数据库（GLW1、GLW2）在方法上进行了改善，使用随机森林回归方法代替多层线性回归方法，得到了全球牛、水牛、马、绵羊、山羊、猪、鸡和鸭的网格分布图（GLW3），空间分辨率为0.083°（赤道上约10 km），数据时间为2010年。国内关于家畜网格化的研究较少，2017年，乔宇鑫等［10］基于克拉克指数模型的网格化技术分析了新源县畜牧活动密度的空间分布格局。随后，吴雪［11］使用随机森林对青海省家畜进行网格化并分析了青海省的草畜平衡情况。
近年来，家畜网格数据作为基础数据，陆续应用到了许多方面，主要包括自然资源合理利用，如基于家畜网格数据评估草畜平衡［12］、估算家畜耗氧量［13］、量化畜牧业用水［14］等；环境影响评估，如基于家畜网格数据量化甲烷的排放［15］；以及传染病风险评估，如基于家畜网格数据评估蓝舌病毒爆发的高风险地区［16］等。可见，家畜网格化研究具有重要的意义，然而就中国区域而言，现有的家畜网格数据存在两个弊端：①统计数据尺度粗糙：现有家畜网格数据在中国区域使用的多为省/市级家畜统计数据；②统计数据时间相对滞后，如GLW2、GLW3在中国区域用于建模的家畜统计数据时间为2000年，远早于数据时间；进而限制了家畜网格化产品在中国区域的精度。对此，本研究收集了2010年县级家畜年末存栏统计数据，以甘肃省羊（山羊和绵羊）为研究对象，利用随机森林回归算法，结合土地覆被、地形、植被、气候以及人类活动等11个环境因子数据，发展家畜分布的空间降尺度方法，获得公里网格上羊的空间分布信息。
2 研究区与数据
2.1 研究区概况
甘肃省位于中国西部地区，介于32°11′~42°57′ N、92°13′~108°46′ E之间，海拔多在1 000 m以上。甘肃省地处黄河中上游，是黄土高原、青藏高原和内蒙古高原三大高原的交汇地带，其地理位置如图1（a）所示。
图1
甘肃省境内地形复杂，以山地和高原为主，区域内地形地貌差异明显；气候类型多样，省内有亚热带季风区、暖温带及温带季风区、暖温带及温带干旱区、以及高寒山地垂直气候带［17］。结合不同的地形地貌特征和气候条件，将甘肃省87个县区划分为六个生态类型区［18］，如图1（b）所示。其中荒漠区气候高寒，风急沙大，山岩裸露，荒漠连片，是难以耕作之地，不适宜畜牧业发展；河西走廊戈壁区地势平坦，光热充足，水资源丰富，农业、畜牧业发展前景广阔；黄土高原丘陵区气候温和，但水资源缺乏，植被稀疏，水土流失较严重，畜牧业在农业中所占比重不大；黄土高原沟壑区气候温和，雨量较多，畜牧业在本区有一定的地位；秦岭中山峡谷区跨长江、黄河两大水系，水资源丰富，区内山大坡陡，地形复杂，农业生产随海拔不同呈现明显的垂直变化，畜牧业比重较低；甘南高原草原草甸区内牧草丰盛，以高山、亚高山草甸草原为主，水资源丰富，畜牧业比重大［19］。
2.2 统计数据与适宜性掩蔽
甘肃省羊（山羊和绵羊）的县级年末存栏统计数据来源于《2011年中国统计年鉴》，共87个县区。对每个县区进行适宜性掩蔽，剔除明确不适合羊生存的区域，将各县区羊的数量统计值除以掩蔽后对应县区的面积，获得调整后的各县区羊的密度数据。
适宜性掩蔽的区域选择较为保守，只排除明确不适合羊生存的区域，掩蔽的区域包括3部分（表1）：永久性水体（中国多时期土地利用土地覆被遥感监测数据集（CNLUCC）［20］中1 km×1 km像元内超过50%的部分被水体覆盖的区域）；中心城区（付晶莹等［21］的中国公里网格人口分布数据集中人口密度超过10 000/km2的区域）；以及自然保护区数据（具有严格的保护措施和对人类活动严格监管的区域），该区域的环境可能适合家畜生存，但由于严格的保护而不会存在家畜分布的情况，因此需要进行人为掩蔽，自然保护区数据来源于联合国环境规划署的世界保护区数据库（World Database on Protected Areas， WDPA）［22］。
表1 适宜性掩蔽区域
Table 1
2.3 所选取的环境因子数据
表2 所选取的环境因子
Table 2
图2
（1）林地、草地、耕地以及荒漠百分比数据：土地覆被数据来源于中国科学院资源环境科学数据中心（http：∥www.resdc.cn）的中国多时期土地利用土地覆被遥感监测数据集（CNLUCC）［20］，空间分辨率为100 m×100 m。林地、草地、耕地以及荒漠百分比数据是由该数据集计算得到的。计算方法为统计每个公里网格中各地物类别的像元个数，将其除以每个公里网格中包含的总像元数。
（2）海拔和坡度数据：海拔和坡度数据来源于国家青藏高原科学数据中心（http：∥data.tpdc.ac.cn）的中国数字高程图［23］，该数据集是根据中国1∶25万等高线和高程点生成的，空间分辨率为1 km，包括DEM、山影、坡度和坡向等数据。
（3）归一化植被指数：归一化植被指数数据来源于MODIS中国合成产品中的TERRA中国500 m NDVI 5 d合成产品MODND1F。使用最近邻法将该产品重采样为1 km×1 km分辨率，并计算每个像元的算数平均值。
（4）白天地表温度：白天地表温度数据来源于中国地表温度月值数据集MODLT1M［24］，空间分辨率为1 km，对全年12幅数据计算对应像元的时间算数平均值。
（6）城市可达性：城市可达性数据来源于Weiss等绘制的2015年全球城市可达性地图［28］，空间分辨率为1 km，该研究根据到达最近的城市中心（≥1 500人/km2）所需的出行时间（不包括航空）来定义城市可达性。
（7）人口密度：人口密度数据来源于付晶莹等［21］的2010年中国公里网格人口分布数据集。
因各种数据的投影方式、坐标系均不尽相同，故给所有数据设置相同的Krasovsky椭球体坐标和Albers投影，并统一空间分辨率为1 km×1 km。
3 研究方法
3.1 随机森林回归算法
回归决策树通过一组决策规则对预测变量（Xi）与响应变量（Y）之间的关系进行建模，其基本思想是：将特征空间递归地划分为依次更小的空间，划分时穷举每一个特征的每一个阈值，通过最小化损失函数来寻找最优切分特征和最优切分点，直到达到预设的终止条件时停止划分［29］。回归决策树的优势在于，它可以对数据中的复杂关系进行建模，并通过决策规则的自适应特性解决预测变量和响应变量之间的非线性关系。但是，当树生长过深时，可能出现过拟合现象。
为了克服此问题，Breiman［30］提出了随机森林算法（Random Forest， RF），包括随机森林分类和随机森林回归算法，其中随机森林回归的本质是多棵相互独立的回归决策树的集合，它通过使用随机子样本集与预测变量的随机子集来构建并生长每棵树。
随机森林回归的基本算法［31］：利用bootstrap方法从总训练样本集中有放回地随机抽取K个不同的训练样本数据子集，生成K棵对应的回归决策树；对于每棵决策树，在每个分裂节点上，通过最小化损失函数确定最佳切分特征集和对应的切分点，每棵树都完整生长而不需要剪枝；重复上述步骤，建立K棵决策树，组成随机森林；利用每棵决策树分别对预测数据进行预测，得到对应的单棵决策树的预测值；对得到的K棵决策树预测结果取平均值，作为最终预测结果。
3.2 家畜分布网格化估算模型
家畜分布网格化估算模型整体流程如图3所示。以甘肃省为研究区，对统一投影和空间分辨率的11个公里网格环境因子数据进行均值分区统计，得到对应的县级环境因子数据，作为构建模型的预测变量。对统计年鉴中甘肃省各县区羊的统计数据与甘肃省各县区矢量边界数据进行属性表连接，计算适宜性掩蔽后各县区的面积，各县区羊的统计数据除以对应各县区的面积得到适宜性掩蔽后的县区尺度羊的密度，并对该密度转换作log10（n+1）变换，以便归一化响应变量的分布［32］，转换后的县区尺度羊的密度值作为构建模型的响应变量。基于处理好的预测变量和响应变量，共87个训练样本，使用随机森林回归算法构建县区尺度上的随机森林回归模型，假设环境因子对羊的空间分布的影响在县级尺度上与公里网格上是一样的［8］，将空间分辨率为1 km的环境因子数据输入在县级尺度训练好随机森林回归模型中，逐像元估算甘肃省公里网格上羊的空间分布情况，并对其进行反对数变换，得到甘肃省羊的公里网格空间分布预测结果。对模型预测准确性的真正验证包括对不同像元上羊的数量进行实地调查，使用调查结果衡量预测精度。然而，家畜数据通常是按行政单元收集的，因此实验对预测结果按县区进行聚合，在县区尺度与统计数据作比较，进而对预测结果进行验证。最后，对各县区羊的统计数据与预测数据的差进行空间插值，将插值结果与预测结果相加，从而对预测结果进行校正。
图3
3.3 精度评估
分别对随机森林回归模型预测得到的公里网格上羊的空间分布情况和使用空间插值校正后羊的空间分布情况以县区为单位进行聚合，得到87个县区羊的预测数量和校正数量，与各县区羊的统计数量进行对比。选用3种评价指标来对预测结果和校正结果进行验证，分别是相关系数（R）、均方根误差（RMSE）以及相对均方根误差（RRMSE），计算公式如下：
其中：i为县区的编号；n为县区的总个数；Pi 是对应县区上随机森林预测的羊的数量或校正后羊的数量；Si 是对应县区羊的统计数量；
RRMSE是均方根误差RMSE与
4 结 果
4.1 参数选择
影响随机森林回归算法性能的参数主要包括随机森林算法中决策树的数量以及构建树时所使用的特征数量。图4展示了使用不同参数构建随机森林预测模型对应的均方根误差。综合考虑时间消耗以及模型的误差，采用决策树的数量为120，特征数量为4来构建随机森林回归预测模型。
图4
4.2 环境因子重要性分析
平均准确度降低程度（Mean Decrease Accuracy， MDA）是将一个变量的取值变成随机数时随机森林预测准确性的下降程度，其下降程度越大表明该变量重要性越大［34］。根据该方法对影响甘肃省羊的空间分布的环境因子进行重要性排序（图5）。结果表明：在甘肃省内，对羊的空间分布影响较大的因子为耕地百分比、海拔、地表温度和坡度。耕地百分比对羊的空间分布影响最大，耕地百分比与农业人口活动密切相关，与野生动物相比，家畜的分布多在人类活动区域内，即使采用放养的方式，家畜的分布也通常保持在居民点附近。其次是海拔，海拔变化会影响区域气候、降水、植被和水源分布，这些因素都会对羊的分布有一定影响，已有研究表明，地形与家畜的分布有较强的关系，家畜更多地聚集在海拔较低的区域［35］。除耕地百分比和海拔以外，地表温度和坡度对甘肃省羊的空间分布影响差异不大，地表温度略高于坡度。
图5
图5
影响羊空间分布的环境因子重要性
Fig.5
Importance of environmental factors influencing the spatial distribution of sheep
4.3 空间分布预测结果
根据3.2节中随机森林回归的公里网格上家畜分布估算模型，将甘肃省羊的县级统计数据（图6（a））空间降尺度为1 km分辨率的网格数据，结果见图6（b）。图6（c）为使用空间插值对随机森林降尺度结果进行校正后的甘肃省公里网格上羊的分布图。整体而言，甘肃省羊的平均密度为44只/km2，空间分布呈现中部聚集模式，即甘肃省中部羊的分布多，两端少。最高值出现在临夏回族自治州的东乡族自治县（219只/km2），最低值分布在陇南市的两当县境内。甘肃省羊的密度高值区主要分布在地势平坦，光热充足，水资源丰富的河西走廊戈壁区和牧草丰盛的甘南高原草原草甸区；羊的密度低值区位于研究区两端的荒漠区和秦岭中山峡谷区，其中荒漠区气候高寒，荒漠连片，人口分布稀疏，不适宜羊的生存。秦岭中山峡谷区水资源丰富，但山大坡陡，地形复杂，同样不适宜畜牧业的发展，故这两个生态类型区内羊的分布较少。为了展示更多细节，从研究区中任意选择两个区域进行局部放大，图6（d）~图6（g）分别为两个选中区域的预测结果和校正结果局部放大图，可以看出，随机森林预测结果可以在公里网格上刻画出羊的空间分布细节信息，使用空间插值对预测结果进行校正后的羊的公里网格分布图对随机森林预测结果图在局部进行了改善，使得公里网格上羊的空间分布更为合理。
图6
图6
甘肃省羊的统计数据与预测结果对比图
审图号：GS （2019）1823
Fig.6
Comparison of statistical data and predicted results of sheep in Gansu Province
4.4 模型精度评估
为验证实验中随机森林回归算法的适用性和解释性，计算随机森林回归得到的甘肃省各县区羊的数量和校正后羊的数量与统计数据的R、RMSE和RRMSE，计算结果如表3所示。可以看出模型在训练集上有较好的表现，各县区统计值与模型估计值之间的相对均方根误差为8.2%，小于10%。对随机森林预测结果（1 km×1 km）进行聚合，得到各县区羊的预测数量与对应县区羊的统计数据对比，其均方根误差分别为0.37和23.3%，同时相关系数R为0.77。对使用空间插值进行校正后的羊的空间分布数据（1 km×1 km）聚合到县区尺度，与对应县区统计数据进行对比，校正结果与统计数据之间的偏差进一步缩小，相关系数和均方根误差分别为0.88和0.24，相对均方根误差为15.1%，介于10%~20%之间，认为模型精度较好。
表3 模型预测精度评估
Table3
| 训练数据 | 0.96 | 0.13 | 8.2% |
| 预测结果 | 0.77 | 0.37 | 23.3% |
| 校正结果 | 0.88 | 0.24 | 15.1% |
将校正后的公里网格上羊的预测结果聚合到各个县区，与统计数据在各县区上做比较，如图7所示。可以看出，羊的估计值与统计值在多数县区上相差较小，但在部分县区，如：位于荒漠区的肃北县、阿克塞县、玉门市和肃州区以及位于河西走廊戈壁区南部祁连山区的肃南县和凉州区，预测结果与统计数据相差较大。这些区域大多面积较大，其中肃北县和肃南县均存在飞地的情况，即隶属于同一行政区管辖但在地理位置上不毗连，在这些县区内地理环境因子变化较大，对其整体取均值作为训练数据不能很好地代表该县区的情况，故而对预测结果造成一定的误差。
图7
图7
县区尺度羊的统计与预测对比图
Fig.7
Comparison of statistics and prediction of sheep at county scale
为进一步探讨不同生态区的模型表现，分别计算了甘肃省6个生态类型区羊的估计值与统计值之间的相关系数R，其与RMSE和RRMSE之间有显著的相关性，R越高，RMSE和RRMSE越小，故只选择R进行展示（图8）。从图中可以看出，在研究区中，有5个生态区的R大于0.5，最高值为黄土高原沟壑区，达到0.98。从区域角度来看，黄土高原沟壑区气候温和，雨量较多，适合畜牧业发展，估计值与统计数据之间的相关性最高（R=0.98）。相对地，荒漠区的县区估计值与统计数据之间的相关系数最低（R=0.23），该区域自然环境恶劣，植被稀少，荒漠连片，展现出不适宜羊分布的特征。此外，计算了各生态类型区内县区的平均面积，发现模型在不同生态区的表现与生态区内县区的规模呈现明显负相关关系，这与DaRe等［36］的研究结论相一致，训练数据的规模会影响降尺度结果。
图8
图8
不同生态区内县区面积与模型表现
Fig.8
The area of counties in different ecological zones and the performance of models
5 结 论
家畜网格化的重要性源于各种研究对高空间分辨率家畜分布数据的需求。家畜网格化数据可以为自然资源合理利用、环境生态保护、传染病风险评估和畜牧业可持续发展等研究提供数据支持。因此，发展家畜网格化技术具有重要的理论和现实意义。
为了获得细尺度上家畜的空间分布情况，探讨随机森林回归算法在家畜统计数据网格化研究中的适用性，实验选取甘肃省作为研究区，结合遥感数据的高空间分辨率和统计数据的准确可靠性优势，使用机器学习技术，开展了基于随机森林回归的羊的统计数据降尺度研究，得到了1 km尺度上羊的空间分布信息，解决了家畜分布数据以行政区为单元，难以与多源数据进行融合分析决策的问题。通过对比降尺度结果与统计数据，认为该方法在整个研究区内表现良好，能够获得精度较高的高空间分辨率羊的空间分布信息，整体能够描绘出统计单元内部羊的空间分布情况，但在不同生态区中表现有一定差异。
虽然文中的网格化实验初步展示出了不错的结果，但该方法仍存在以下缺陷：存在潜在假设，即不同尺度上环境因子和羊的分布之间的关系不变；很难获取实际数据对结果进行验证，比如以家庭为单位的家畜养殖统计数据，实地调察可以解决这一问题，但成本高昂，特别是大面积区域；另外，影响家畜羊的空间分布的因素很复杂，除去本实验中所用的环境因子，当地的农业结构、经济发展水平、政府政策等都会对预测结果产生影响。后续应进一步修正预测变量集，考虑添加更多社会经济和人类活动因子，对模型加以改进。
参考文献
The State of Food and Agriculture-Livestock in the Balance
［EB/OL］. ，
World Livestock 2011-Livestock in Food Security
［EB/OL］. ，
Smart investments in sustainable food production： revisiting mixed crop-livestock systems
［J］.
Molecular basis for the generation in pigs of influenza a viruses with pandemic potential
［J］.
The pig as a mixing vessel for influenza viruses： human and veterinary implications
［J］.
Mapping the global distribution of livestock
［J］.
Gridded Livestock of the World
［EB/OL］. ，
Using random forest to improve the downscaling of global livestock census data
［J］.
Global distribution data for cattle， buffaloes， horses， sheep， goats， pigs， chickens and ducks in 2010
［J］.
Research on gridding of livestock spatial density based on multi-source information
［J］.
基于多源数据的家畜活动密度分布格网化研究
［J］.
Analysis on the spatiotemporal changes of livestock in Qinghai grassland in recent 30 years
［D］.
Global assessment of grassland carrying capacities and relative stocking densities of livestock
［J］.
Estimation of gridded atmospheric oxygen consumption from 1975 to 2018
［J］.
Where is the planetary boundary for freshwater being exceeded because of livestock farming？
［J］.
Methane emissions from livestock in East Asia during 1961-2019
［J］.
Predicting temperature-dependent transmission suitability of bluetongue virus in livestock
［J］.
Analysis of causes and regional distribution of ecological environment change in Gansu Province
［J］.
甘肃省生态环境变化成因及区域分布分析
［J］.
Investigation on the relationship between environmental iodine level and iodine nutrition of women of childbearing age in Six Ecological Regions of Gansu Province
［D］.
甘肃省六大生态区外环境碘水平与育龄期妇女碘营养关系调查研究
［D］.
The present situation and preliminary plan of animal husbandry division in Gansu Province
［J］.
甘肃省畜牧业现状及区划的初步设想
［J］.
Multi-Period Land Use and Land Cover Remote Sensing Monitoring Data Set in China
［DB/OL］.
China's Kilometer Grid Population Distribution Dataset
［DB/OL］.
The World Database on Protected Areas
［DB/OL］.
Digital Elevation Model of China （1KM）
［DB/OL］.
Data Set of Monthly Surface Temperature in China from 2000 to 2014
［DB/OL］.
The first high-resolution meteorological forcing dataset for land process studies over China
［J］.
On downward shortwave and long-wave radiations over high altitude regions： observation and，odeling in the Tibetan Plateau
［J］.
China Meteorological Forcing Dataset （1979-2018）
［DB/OL］.
A global map of travel time to cities to assess inequalities in accessibility in 2015
［J］.
Newer classification and regression tree techniques： bagging and random forests for ecological prediction
［J］.
Random forests
［J］.
Evaluation of empirical models for predicting monthly mean horizontal diffuse solar radiation
［J］.
Using "random forest" for classification and regression
［J］.
随机森林模型在分类与回归分析中的应用
［J］.
Cattle grazing distribution patterns related to topography across diverse rangeland ecosystems of North America
［J］.
Downscaling livestock census data using multivariate predictive models： sensitivity to modifiable areal unit problem
［J］.
