CN113902580B - 一种基于随机森林模型的历史耕地分布重建方法 - Google Patents

Abstract

本发明公开了一种基于随机森林模型的历史耕地分布重建方法，包括以下步骤：S1:历史耕地数据资料查询，空间数据搜集。S2:对各时段耕地面积数据进行校准、插值和验证，以获得统一标准下的分县耕地面积数据集。S3:选取影响耕地分布的自然和社会区位因子，结合当代遥感影像，对随机森林分类模型进行训练；基于各影响因素对模型的贡献程度，筛选模型特征参数。S4:基于训练得到的随机森林分类模型以及预处理后的耕地面积，对耕地面积进行空间分配，并进行迭代。本发明优点是：实现由耕地重建垦殖率数据向布尔型数据的转变，大幅度提高了耕地空间布局重建的分辨率与精度，解决了与现有及未来高分辨率遥感影像的衔接。

Description

一种基于随机森林模型的历史耕地分布重建方法

技术领域

本发明涉及耕地重建技术领域，特别涉及一种能解决高维输入参数、高效可靠的基于随机森林模型的历史耕地空间分布重建方法。

背景技术

近百年来，由于人口激增对食物和原材料的需求不断增长，土地利用和覆被变化(LUCC)过程发生了巨大的变化。耕地是土地资源利用中最为活跃的方式之一，对农业粮食生产、全球变化、区域生态系统服务具有重大影响。长时间耕地时空数据集可以为理解区域农业发展轨迹提供一个全面的视角，对于提出未来环境或经济战略具有重要意义。同时，重建的数据集也可应用于国家和全球碳排放核算和陆地生态系统生态模拟中。然而，目前所广泛使用的是基于传统方法开发的中低分辨率农田数据集。因此，迫切需要一个长时序、高分辨率的农田数据集，以提高生态模拟精度并为未来决策服务。

近年来，在LUCC、BIOME300等研究计划的推动下，定量重建历史土地覆被的工作取得重大突破，尤其是一些区域和全球土地覆被数据集的建立。耕地空间布局的重建建立在数量重建的基础上，即将数量重建的结果在研究区范围内进行分配。其中最常用的方法是“自上而下”的静态配置模型，该模型根据土地适宜性、人口密度或当代土地利用模式将定量重建的历史LUCC数据分配到相应的空间网格中，具有代表性的全球尺度产品之一是全球环境历史数据库(HYDE)，但它从人口和消费等社会经济参数中间接得出，因此在局部地区的分辨率较低。另一种方法是通过“自下而上”的空间演化模型，在传统配置模型的基础上考虑了元胞状态或人类土地利用活动方式。在局部或区域尺度上，通过元胞自动机模型或融入了耕地利用主体行为规则的多智能体模型，通过循环迭代完成耕地空间布局的模拟。当前，国内耕地格局重建主要用到的模型或方法包括网络化重建模型、分区网络化模型、垦殖倾向指数模型、基于土地适宜性的空间分配模型、比例修正法等。

然而现有广泛使用的传统分配方法无法处理高维的数据、存在主观效应；它不仅需要对影响土地利用的因素进行选择、量化和归一化，还需要采用数学语言建立抽象函数或一系列模型来描述特征，计算成本较高，并且在时间或空间尺度上的扩展性有限。近来，机器学习算法正在成为模拟土地利用变化最强大的工具。随机森林(RF)模型是利用多棵树对样本进行训练并预测的一种分类器，该分类器可以输入大量变量，快速学习后输出高准确度的分类或回归结果，同时评估变量的重要性，不会产生过拟合问题。RF模型在处理高维输入数据集方面被证明是有效的，具有不考虑多重共线性、避免过拟合问题、运算速度快和可靠性高等优点，也已成功应用于预测短期土地利用变化。随着辅助数据类型和数量的增加，更需要RF等机器学习方法来处理这些高维、高分辨率的数据。RF算法通过客观训练，使得模拟结果从百分比形式向布尔型数据形式的转化成为可能；同时随机森林算法对变量共线性问题不敏感，在处理高维数据时无需进行降维，也无需对自变量进行标准化处理；随机森林算法的随机性大大减少了模型的过拟合现象；此外，随机森林算法可以对变量进行重要性排序，利于对影响耕地空间分布的因子贡献率进行分析解释。将随机森林为代表的机器学习方法引入到耕地空间布局建模方法中，是耕地重建领域与新兴信息技术结合的一次尝试，充实了耕地空间布局重建的方法体系。

发明内容

本发明针对现有技术无法处理高维数据、主观色彩较重、变量解释性差、百分比形式的预测结果难以精确展现耕地空间分布的缺点，提供了一种基于随机森林模型的历史耕地分布重建方法，解决了现有技术存在的缺陷。

为了实现以上发明目的，本发明采取的技术方案如下：

一种基于随机森林模型的历史耕地分布重建方法，包括以下步骤：

S1:历史耕地数据资料查询，空间数据搜集。分区县搜集研究区近代时期地方志、近代和现代时期农业普查、调查资料、统计年鉴中的耕地面积记载数据；同时，搜集对耕地分布具有影响的自然和社会因素的矢量和栅格数据，包括地势地貌、气候、河流水文、土壤、自然灾害、道路、农村居民点、历史人口数据，并对数据进行预处理。

S2:耕地面积重建。对各个时段耕地面积数据进行校准、插值和验证，以获得统一标准下的分县耕地面积数据集。

S3:RF模型训练。在对历史耕地进行重建前，就研究区实际情况，选取影响耕地分布的自然和社会区位因子，结合当代遥感影像，对随机森林分类模型进行训练，以保障模型拟合优度；同时，基于各影响因素对模型的贡献程度，筛选模型特征参数。

S4:耕地历史分布重建。基于训练得到的各区县随机森林分类模型以及预处理后的耕地面积，在ArcGIS Pro下的Python3环境中，对各区县耕地面积进行空间分配，并进行迭代，最终得到各年份历史耕地空间分布。

进一步地，S3的具体子步骤如下：

S31：将所有空间数据转换成统一投影坐标系以及参考椭球体。

S32：对历史耕地分布具有影响的自然和社会区位因子进行空间量化处理。得到10个影响耕地分布的当代和历史自然和社会区位因素，包括：高程、坡度、地形起伏度、气候生产潜力、邻域耕地数量、与农村聚落的距离、与河流距离、洪涝风险、土壤侵蚀模数和土壤养分指数。

S33：训练随机森林分类模型。在随机森林分类模型中，以当代遥感解译的耕地分布作为因变量，影响耕地分布的自然和社会区位因素为自变量，对模型进行训练拟合。在训练模型前，按照8：2的比例将数据划分为训练集和测试集，便于对模型性能进行科学评估。在随机森林分类训练时，需要对n_estimators和max_features两个重要参数进行设定。n_estimators为随机森林建模过程中决策树的生长规模，随着n_estimators的增大，模型的收敛程度增加到一定水平时会趋于平稳，但模型运行时长也会相应增加。max_features为决策树的叶子节点在探寻最佳分枝时设定的变量个数。

S34：利用随机森林分类模型，通过对各影响因子对模型拟合的重要性进行排序。

进一步地，S4的具体子步骤如下：

S41：由于统计来源的耕地面积和遥感耕地面积存在差异，使用调和算法将耕地面积重建结果调整至遥感数据水平，调和算法如下：

以当代遥感影像年份作为起始年进行回溯调和，其中t1是当前年份；t2是当前年份的上一年(t2＝t1-1)；k是县份指数；

是k县在t1年的耕地面积重建结果；

是k县在t2年的耕地面积重建结果；

是k县在t1年调和后的耕地面积；

是k县在t2年调和后的耕地面积；α(k)是相对异常的权重，而1-α(k)是绝对异常的权重。

S42：将耕地起始年网格化数据、影响耕地空间分布的特征参数集、训练择优后的随机森林分类模型以及待重建年份的区县耕地面积导入到分配模型中。

S43：使用随机森林分类模型估算各网格单元在重建年份上的耕地分布概率，结合该年相对于上一重建年份的耕地变化数量，通过逐年回溯迭代确定耕地变化的位置。针对耕地数量增多的情况，按耕地分布概率由高到低确定新增耕地网格，同时对潜在耕地分布范围以外的区域进行限制分配；针对耕地数量减少的情况，按耕地分布概率由低到高将网格转为非耕地。

S44：将各区县耕地空间分布重建结果进行镶嵌，得到研究区域各重建时间断面的耕地网格化空间数据集。

与现有技术相比，本发明的优点在于：

1、模拟得到高精度耕地历史空间分布数据集，实现耕地重建结果由垦殖率数据向布尔型数据的转变。

2、大幅度提高了耕地空间布局重建的分辨率与精度；

3、将耕地面积重建结果与高分辨率遥感数据进行调和统一，实现了耕地分布时空连续性，丰富了历史时期土地覆被重建的研究案例，解决了与现有及未来高分辨率遥感影像的衔接。

附图说明

图1为本发明实施例耕地空间分布网格化重建方法的流程图；

图2为本发明实施例中各区县测试集上的AUC-ROC精度分析图；

图3为本发明实施例中影响因素重要性排序；

图4为本发明实施例中1911、1933、1945、1957、1960、1980、2000和2010年沱江流域30m×30m耕地空间分布；

图5为本发明实施例中1980年重建与遥感影像耕地垦殖率对比图：a.1980年重建耕地垦殖率，b.1980年遥感影像耕地垦殖率，c.1980年重建和耕地遥感影像耕地垦殖率差异比率；

图6为本发明实施例中1980年耕地分布1km重建垦殖率与遥感数据绝对差异直方图。

图7为本发明实施例中与现有最为常用的土地适宜性分配模型重建的三种耕地数据集在1980年重建的结果对比图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下根据附图并列举实施例，对本发明做进一步详细说明。

如图1所示，一种基于随机森林模型的历史耕地分布重建方法，包括以下步骤：

S1：根据历史时期地方志、统计年鉴、农业普查等数据源，获取沱江流域27个区县1911-2010年的耕地数量、幅员面积、人口数量等数据。同时，收集提取历史和当代沱江流域土地利用数据、行政区划矢量数据、对历史耕地分布具有影响力的沱江流域自然环境及社会区位因子空间数据，包括DEM数据、气候数据、土壤数据、历史居民点地理数据、土壤侵蚀数据、气象站点数据。

S2：分析历史不同时期耕地变化特征，选取典型断面，对各个时段耕地面积数据进行校准、插值和验证，以获得统一标准下的分县耕地面积数据集。

S3：在对历史耕地进行重建前，就沱江流域实际情况，选取影响耕地分布的自然和社会区位因子，结合当代遥感影像，对随机森林分类模型进行训练，以保障模型拟合优度；同时，基于各影响因素对模型的贡献程度，筛选模型特征参数，探讨沱江流域耕地空间布局的影响机制。将耕地面积统一至遥感数据标准，以保障重建成果与后续遥感影像的衔接。

S3的具体子步骤如下：

S31：选取耕地空间布局的影响因素并进行预处理。结合沱江流域实际情况，选取高程、坡度、地形起伏度、气候生产潜力、邻域耕地数量、与农村聚落的距离、与河流距离、洪涝风险、土壤侵蚀模数和土壤养分指数作为自变量。

S32：训练随机森林分类模型。在随机森林分类模型中，以高程、坡度、地形起伏度、气候生产潜力、邻近耕地、与最近的居民点距离、与河流距离、洪涝风险、土壤侵蚀和土壤肥力作为自变量，以2017年遥感解译的耕地分布作为因变量，对模型进行训练拟合。在训练模型前，按照8：2的比例将数据划分为训练集和测试集，便于对模型性能进行科学评估。在随机森林分类训练时，采用网格搜索的方式确定nestimators的取值为200，max_features的取值为3，其他参数为默认，在该参数组合下，模型的性能较稳定且计算量适中。

S33：训练结果精度检验。通过对测试集的训练结果进行精度检验，可知该模型在各区县的拟合优度均高于0.8，为了直观地展示随机森林模型在各区县的拟合效果，分别绘制了各区县测试集上的ROC曲线图，由图2可知，各区县AUC-ROC值均高于0.85，随机森林分类效果较优。

S34：影响因素的重要性排序。图3显示了在随机森林分类模型中耕地空间分布各影响因素的重要程度排序情况，结果显示邻域耕地数量的特征重要性显著高于其它9个因子的重要程度，高程、气候生产潜力、到农村聚落距离、洪涝危险性、坡度和到河流距离这6个因子的相对重要性都介于0.05-0.1之间，地形起伏度、土壤侵蚀和土壤养分指标的重要性值均在0.05以下，为保障模型的分类和拟合精度，本研究将10个特征参数予以保留。

S4：过去百年沱江流域耕地分布重建。基于训练得到的各区县随机森林分类模型以及预处理后的耕地面积，在ArcGIS Pro下的Python3环境中，通过构建网格化分配模型对各区县耕地面积进行空间分配，然后将网格化重建的结果与遥感解译图像以及典型的全球重建结果进行对比。具体包括：

S41：将耕地面积统一至遥感数据标准，以保障重建成果与后续遥感影像的衔接。

S42：将2017年耕地网格化数据、影响耕地空间分布的特征参数集、训练择优后的随机森林分类模型以及待重建年份的区县耕地面积导入到分配模型中。其中，假设自然环境因子在百年尺度上保持不变，而邻域耕地数量、距农村聚落距离和人口密度数据则根据迭代年份实时更新。

S43：使用随机森林分类模型估算各网格单元在重建年份上的耕地分布概率，结合该年相对于上一重建年份的耕地变化数量，通过逐年回溯迭代确定耕地变化的位置。针对耕地数量增多的情况，按耕地分布概率由高到低确定新增耕地网格，同时对潜在耕地分布范围以外的区域进行限制分配；针对耕地数量减少的情况，按耕地分布概率由低到高将网格转为非耕地。

S44：将各区县耕地空间分布重建结果进行镶嵌，得到沱江流域百年尺度上共8个时间断面(1911年、1933年、1945年、1957年、1960年、1980年、2000年、2010年)的耕地网格化空间数据集，空间分辨率为30m×30m，如图4所示。

步骤S5，重建结果检验。为了定量评估重建方法的科学性与重建结果的准确性，将基于耕地网格化模型分配的1980年流域耕地垦殖率数据和同期遥感耕地垦殖率空间布局进行对比(图5)，通过分析二者差异来检验模型的合理性与可行性。为了更直观地展示耕地网格化模型与遥感影像的差异分布区间，绘制了垦殖率绝对差异的直方图(图6)。

在本发明提供的这种基于随机森林模型的耕地空间分布网格化重建方法中，采用的随机森林模型可以使得模拟结果从百分比形式向布尔型数据形式的转化；在处理高维数据时无需进行降维，也无需对自变量进行标准化处理；随机森林算法的随机性大大减少了模型的过拟合现象；随机森林算法可以对变量进行重要性排序，利于对影响耕地空间分布的因子贡献率进行分析解释。在一定程度上解决了传统技术指标权重的主观色彩较重、百分比形式的预测结果难以精确展现耕地空间分布的缺点，大幅度提高了耕地网络化重建的分辨率与精度，也解决了与现有及未来高分辨率遥感影像的衔接，构建了一套新的耕地网格化重建方法体系。

与现有最为常用的土地适宜性分配模型重建的三种耕地数据集在1980年重建的结果对比(图7)。HYDE3.2数据集中相同垦殖率的网格连片分布，而本本发明由于在网格化分配时以县级耕地面积作为数量限制，局地垦殖率差异明显。与ChinaCropland数据集相比，ChinaCropland数据集的垦殖率空间分布差异较HYDE3.2数据集明显，但受5km分辨率限制，其对较小区域耕地分布情况的捕捉仍不如本研究精准。与中国传统农区数据集相比，虽然二者的垦殖率网格单元都为1km，但传统农区数据集对垦殖率的展现较极端，几乎只包含20％以下和80％以上垦殖率区间，而本研究重建出的耕地分布更精细，基于此得出的网格垦殖率更贴合区域实际。

综上可知，本发明重建结果与各典型重建数据集在空间分布格局上总体趋势一致，同时克服了典型数据集在细节刻画方面的分辨率限制，更适合作为区域尺度的土地利用基础数据集，以支撑区域的气候与生态环境效应模拟等研究。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的实施方法，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (1)

1.一种基于随机森林模型的历史耕地分布重建方法，其特征在于，包括以下步骤：

S1:历史耕地数据资料查询，空间数据搜集；分区县搜集研究区近代时期地方志、近代和现代时期农业普查、调查资料、统计年鉴中的耕地面积记载数据；同时，搜集对耕地分布具有影响的自然和社会因素的矢量和栅格数据，包括地势地貌、气候、河流水文、土壤、自然灾害、道路、农村居民点、历史人口数据，并对数据进行预处理；

S2:耕地面积重建；对各个时段耕地面积数据进行校准、插值和验证，以获得统一标准下的分县耕地面积数据集；

S3:RF模型训练；在对历史耕地进行重建前，就研究区实际情况，选取影响耕地分布的自然和社会区位因子，结合当代遥感影像，对随机森林分类模型进行训练，以保障模型拟合优度；同时，基于各影响因素对模型的贡献程度，筛选模型特征参数；

S3的具体子步骤如下：

S31：将所有空间数据转换成统一投影坐标系以及参考椭球体；

S32：对历史耕地分布具有影响的自然和社会区位因子进行空间量化处理；得到10个影响耕地分布的当代和历史自然和社会区位因素，包括：高程、坡度、地形起伏度、气候生产潜力、邻近耕地、与最近的居民点距离、与河流距离、洪涝风险、土壤侵蚀和土壤肥力；

S33：训练随机森林分类模型；在随机森林分类模型中，以当代遥感解译的耕地分布作为因变量，影响耕地分布的自然和社会区位因素为自变量，对模型进行训练拟合；在训练模型前，按照8：2的比例将数据划分为训练集和测试集，便于对模型性能进行科学评估；在随机森林分类训练时，需要对n_estimators和max_features两个重要参数进行设定；n_estimators为随机森林建模过程中决策树的生长规模，随着n_estimators的增大，模型的收敛程度增加到一定水平时会趋于平稳，但模型运行时长也会相应增加；max_features为决策树的叶子节点在探寻最佳分枝时设定的变量个数；

S34：利用随机森林分类模型，通过对各影响因子对模型拟合的重要性进行排序；

S4:耕地历史分布重建；基于训练得到的各区县随机森林分类模型以及预处理后的耕地面积，在ArcGIS Pro下的Python3环境中，对各区县耕地面积进行空间分配，并进行迭代，最终得到各年份历史耕地空间分布；

S4的具体子步骤如下：

S41：由于统计来源的耕地面积和遥感耕地面积存在差异，使用调和算法将耕地面积重建结果调整至遥感数据水平，调和算法如下：

以当代遥感影像年份作为起始年进行回溯调和，其中t1是当前年份；t2是当前年份的上一年，t2＝t1-1；k是县份指数；

是k县在t1年的耕地面积重建结果；

是k县在t2年的耕地面积重建结果；

是k县在t1年调和后的耕地面积；

是k县在t2年调和后的耕地面积；α(k)是相对异常的权重，而1-α(k)是绝对异常的权重；

S42：将耕地起始年网格化数据、影响耕地空间分布的特征参数集、训练择优后的随机森林分类模型以及待重建年份的区县耕地面积导入到分配模型中；

S43：使用随机森林分类模型估算各网格单元在重建年份上的耕地分布概率，结合该年相对于上一重建年份的耕地变化数量，通过逐年回溯迭代确定耕地变化的位置；针对耕地数量增多的情况，按耕地分布概率由高到低确定新增耕地网格，同时对潜在耕地分布范围以外的区域进行限制分配；针对耕地数量减少的情况，按耕地分布概率由低到高将网格转为非耕地；

S44：将各区县耕地空间分布重建结果进行镶嵌，得到研究区域各重建时间断面的耕地网格化空间数据集。

Images