CN108053405B - 一种耕地制图方法 - Google Patents

Abstract

一种耕地制图方法，其通过对至少三个不同的耕地产品i进行处理，制作出特定分辨率M的耕地分布图，然后进行融合处理获得最终耕地图表。本发明所提供的一种耕地图表制作方法，通过对多套现有技术平台提供的耕地产品进行融合，所获得的耕地分布图总体精度和可靠性大大提升，且耕地数量较为准确，与统计数据的一致性较好。

Description

一种耕地制图方法

技术领域

本发明涉及一种遥感图像处理方法，特别是一种利用已有遥感制图数据融合制作耕地分布图的方法。

背景技术

耕地是人类社会赖以生存和发展的基础。准确及时的耕地空间分布和动态变化信息，是解决农作物生长监测、产量估算、粮食安全评估等农业土地***关注的热点问题的关键。随着卫星传感器的不断发展和计算机制图技术的进步，遥感技术已经成为包括耕地在内的地表覆盖类型空间分布信息获取的重要手段(Ban等，2015；Chen等，2013)。从全球到区域尺度，目前已有多套提供地表覆盖信息数据的技术平台，如1km分辨率的IGBP-DISCover(International Geosphere Biosphere Programme,Data and Information SystemsCover)数据集(Loveland等，2000)，300米分辨率的GlobCover和CCI-LC(Climate ChangeInitiative Land Cover)数据集(Bontemps等，2009；Defourny,2017)，30米分辨率的GlobeLand30数据集等(Chen等2015)。这些数据都可以提供耕地分布信息，即耕地数据的影像资料，也就是都可以提供相应的耕地产品，然而，上述技术平台所提供的耕地产品多是针对气候变化模型等应用需求而研发并获得的，所以不能满足粮食安全、作物监测和产量估算等农业应用的要求(Fritz等，2013)。

首先，已有的不同技术平台所提供的耕地产品中耕地类型定义不一致，如GlobCover 2009和CCI-LC采用FAO(Food and Agriculture Organization of the UnitedNations，***粮农组织)的耕地类型定义方法，MODIS C5(MODIS Collection 5)采用IGBP分类体系的耕地定义方法，而GlobeLand30和FROM-GLC根据30米遥感影像特征对耕地进行了自定义，因此不同耕地产品所提供的耕地数据之间的可比性较差(Fritz等，2014)。其次，从影像上估算的耕地面积和统计数据相差较大，遥感影像上混合像素比较普遍，现有的耕地产品仅能提供像素是耕地或者不是耕地的信息，不能提供像素内耕地比例，因此估算的耕地面积普遍高于统计数据(Schepaschenko等，2015)。

发明内容

本发明要解决的技术问题是提供一种耕地制图方法，以减少或避免前面所提到的问题。

本发明提供了一种耕地制图方法，其通过对至少三个不同的耕地产品i进行处理，制作出特定分辨率M的耕地分布图，其包括如下步骤：

步骤A：以国家为单元，估算多个所述耕地产品i在国家j的耕地面积A_i,j，对于每个耕地产品i，采用经纬度地理坐标进行数据预处理，将每个耕地产品i处理为所述特定分辨率M的耕地分布图，之后基于等面积投影原理计算每个网格对应的实际面积，然后利用以下公式计算所述耕地产品i中的耕地面积A_i,j：

A_i,j＝∑p_i,j×a_i,j

上式中，p_i,j是预处理后的耕地比例，a_i,j是对应的含有耕地像素的网格的实际面积。

步骤B：对于步骤A获得的多个所述特定分辨率M的耕地分布图中的每个相同坐标的像素，构建不同产品的权重打分表，利用以下公式计算每个耕地产品i的面积差异自适应指数Diff_i,j：

上式中，A_FAO,j是国家j的耕地统计数据。指数Diff_i,j的值越小，表明耕地产品i和统计数据越接近，那么将该耕地产品赋予较高的权重；否则，和统计数据差别越大，给该耕地产品赋予较低的权重。

步骤C:按照以下表格结构，构建耕地置信打分规则表，每个耕地产品所对应取值规则为，该耕地产品认为该像素为耕地时取值为1，否则取值为0，耕地一致量表示对同一像素同时有多少个耕地产品认为该像素为耕地。置信分是从1开始顺序增加的整数，值越大表示置信度越低。耕地一致量相同时，耕地产品权重越高，则置信度越高，相应的置信分越低。

步骤D：基于步骤C构建的像素耕地置信打分规则表，根据行政区域(例如：国家、省、地市三级)分多级对数据进行处理，在每一级耕地统计数据所对应的行政区划相应的地图范围，从置信度较高的分值1开始，挑选相应分值的像素，计算分值1对应的耕地面积，然后和统计数据做比较，如果小于所述行政区域范围的耕地统计数据，则继续累加置信度较高的2对应的耕地面积，直到和统计数据的耕地面积比较吻合为止，以此确定出该区域的耕地融合初步结果。提供多级耕地统计数据，根据该步骤得到多级耕地融合初步结果。

步骤E：按照行政区划级别的从低至高的顺序，逐步整合步骤D的初步融合结果，得到最终的耕地制图结果。整合规则为：

(1)在某个下级行政区域内，根据步骤D的初步融合结果，分别计算上一级行政区域的地图中对应本级行政区域范围的融合结果的耕地估算面积，以及本级行政区域融合结果的耕地估算面积，如果上一级行政区域的耕地估算面积＜本级耕地估算面积，则选取本级行政区域的耕地融合结果作为上一级行政区域中对应本级行政区划的地图范围的耕地融合结果，否则，原上一级行政区域中对应本级行政区划的地图范围的融合结果不变。

(2)对于某个下级行政区域，如果没有统计数据，在上一级行政区划的地图中按照以下公式计算该下级行政区域的统计数据：

剩余区域统计数值(即没有统计数据的所有下级行政区划)＝上一级耕地统计数据-【步骤(1)所获得的各对应下一级耕地融合结果所获得的估算面积之和】，

(3)根据步骤(2)所获得的剩余区域统计数值，利用步骤D的方法对该区域的像素重新进行校正，获得该区域的耕地融合结果，并获得最终的耕地制图数据。

本发明所提供的一种耕地图表制作方法，通过对多套现有技术平台提供的耕地产品进行融合，所获得的耕地分布图总体精度和可靠性大大提升，且耕地数量较为准确，与统计数据的一致性较好。

附图说明

以下附图仅旨在于对本发明做示意性说明和解释，并不限定本发明的范围。其中，

图1为根据本发明的一个具体实施例的一种耕地制图方法的耕地数据初步融合过程的原理示意图；

图2为图1的方法中所得到的多级融合结果的整合规则原理示意图；

图3为根据图1的方法和图2的规则获得的中国耕地融合结果的耕地示意图；

图4为图3的耕地融合结果和五套耕地产品数据的精度对比分析示意图；

图5显示了对比分析图3的耕地融合结果和五套耕地产品数据的耕地估算面积和统计数据的均方根误差和相关系数的示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现说明本发明的具体实施方式。

本发明提供了一种耕地制图方法，其通过对至少三个不同的耕地产品i进行处理，制作出特定分辨率M的耕地分布图，其包括如下步骤：

步骤A：以国家为单元，估算多个所述耕地产品i在国家j的耕地面积A_i,j，对于每个耕地产品i，采用经纬度地理坐标进行数据预处理，将每个耕地产品i处理为所述特定分辨率M的耕地分布图，之后基于等面积投影原理计算每个网格对应的实际面积，然后利用以下公式计算所述耕地产品i中的耕地面积A_i,j：

A_i,j＝∑p_i,j×a_i,j

上式中，p_i,j是预处理后的耕地比例，a_i,j是对应的含有耕地像素的网格的实际面积。

步骤B：对于步骤A获得的多个所述特定分辨率M的耕地分布图中的每个相同坐标的像素，构建不同产品的权重打分表，利用以下公式计算每个耕地产品i的面积差异自适应指数Diff_i,j：

上式中，A_FAO,j是国家j的耕地统计数据。指数Diff_i,j的值越小，表明耕地产品i和统计数据越接近，那么将该耕地产品赋予较高的权重；否则，和统计数据差别越大，给该耕地产品赋予较低的权重。

步骤C:按照以下表格结构，构建耕地置信打分规则表，每个耕地产品所对应取值规则为，该耕地产品认为该像素为耕地时取值为1，否则取值为0，耕地一致量表示对同一像素同时有多少个耕地产品认为该像素为耕地。置信分是从1开始顺序增加的整数，值越大表示置信度越低。耕地一致量相同时，耕地产品权重越高，则置信度越高，相应的置信分越低。

步骤D：根据步骤C构建的像素耕地置信打分规则表，根据行政区域(例如：国家、省、地市三级)分多级对数据进行处理，在每一级耕地统计数据所对应的行政区划相应的地图范围，从置信度较高的分值1开始，挑选相应分值的像素，计算分值1对应的耕地面积，然后和统计数据做比较，如果小于所述行政区域范围的耕地统计数据，则继续累加分值较高的2对应的耕地面积，直到和统计数据的耕地面积比较吻合为止，以此确定出该区域的耕地融合初步结果。提供多级耕地统计数据，根据该步骤得到多级耕地融合初步结果。

步骤E：按照行政区划级别的从低至高的顺序，逐步整合步骤D的初步融合结果，得到最终的耕地制图结果。整合规则为：

(1)在某个下级行政区域内，根据步骤D的初步融合结果，分别计算上一级行政区域的地图中对应本级行政区域范围的融合结果的耕地估算面积，以及本级行政区域融合结果的耕地估算面积，如果上一级行政区域的耕地估算面积＜本级耕地估算面积，则选取本级行政区域的耕地融合结果作为上一级行政区域中对应本级行政区划的地图范围的耕地融合结果，否则，原上一级行政区域中对应本级行政区划的地图范围的融合结果不变。

(2)对于某个下级行政区域，如果没有统计数据，在上一级行政区划的地图中按照以下公式计算该下级行政区域的统计数据：

剩余区域统计数值(即没有统计数据的所有下级行政区划)＝上一级耕地统计数据-【步骤(1)所获得的各对应下一级耕地融合结果所获得的估算面积之和】，

(3)根据步骤(2)所获得的剩余区域统计数值，利用步骤D的方法对该区域的像素重新进行校正，获得该区域的耕地融合结果，并获得最终的耕地制图数据。

下面以2010年的中国区域的耕地制图为例，具体说明本发明的方法。

在本具体实施例中，选用五套2010基准年的全球地表覆盖数据提供的耕地产品进行融合，包括GlobeLand30、CCI-LC、GlobCover 2009、MODS C5和Unified cropland，各耕地产品的详细信息如下表所示。

提供国家、省和市三级统计数据。其中，国家级耕地面积统计数据来自FAOSTAT数据库(http://www.fao.org/faostat/en/#data/GC)，该数据中收集了全球200多个国家连续多年的耕地面积统计数据。次国家级(即省级和地市级)统计数据来自国际食物政策研究所(International Food Policy Research Institute，IFPRI)，该数据由某一地区的作物收获面积除以复种指数得到。在全球尺度，国家级统计数据覆盖了全球所有国家，省级统计数据覆盖了全球大部分区域，市级统计数据覆盖了全球部分区域。由于在某些区域只收集到部分的作物收获面积以及复种指数信息，因此省级和市级的耕地面积统计数据有可能小于实际的耕地面积(You等，2014)。

根据上面提供的五个耕地产品和三级统计数据，采用本发明的方法制作中国区域的耕地分布图，具体步骤为：

步骤A：估算各耕地产品i在中国区域j的耕地面积。对于每个耕地产品，采用WGS84地理坐标，将每个耕地产品处理为特定分辨率M为500米(即0.41667°)的耕地分布图，在重采样过程中计算耕地比例。以250米分辨率的Unified Cropland为例，从250米重采样到500米就是计算500米网格内，对应的四个250米网格的耕地平均比例p_i,j，即：

按照等面积投影原理计算每个网格对应的实际面积a_i,j：

其中R是地球半径(≈6371KM)，λ₁和λ₂是该网格对应的经度范围，

和

是该网格对应的纬度范围。

利用以下公式计算所述耕地产品i中的耕地总面积A_i,j：

A_i,j＝∑p_i,j×a_i,j

步骤B，利用以下公式计算每个耕地产品i的在每个相同坐标的像素的面积差异自适应指数Diff_i,j：

上式中，A_FAO,j是中国区域j的耕地统计数据。指数Diff_i,j的值越小，表明耕地产品i和统计数据越接近，那么将该耕地产品赋予较高的权重；否则，和统计数据差别越大，给该耕地产品赋予较低的权重。

步骤C，按照以下表格结构构建像素耕地置信打分规则表，

上表中，每个耕地产品所对应的列的取值规则为，该耕地产品认为该像素为耕地时取值为1，否则取值为0。耕地一致量表示对同一像素同时有多少个耕地产品认为该像素为耕地，置信分为从1开始顺序增加的整数。具体来说，假设某个区域内，耕地产品#1的自适应指数Diff_i,j最小，设置较高的权重，其次是#2、#3和#4，耕地产品#5的自适应指数值Diff_i,j最大，设置最低的权重。根据这5套耕地产品数据建立的表说明，对于某一像素，当5套耕地产品认为都是耕地，即耕地一致量是5的时候，该像素的耕地可靠性最好，在打分表中设置分值为1；当耕地一致量是4，即4套认为某像素是耕地，根据输入产品的权重，对不同的组合设置打分表分值2-6。类似的，建立所有5个耕地产品组合的打分表，整体分值在1到32之间，分值越小，认为是耕地的可能性越大；反之，则耕地可能性越小。

步骤D，对于所提供的国家、省级和地市级这三级统计数据，根据步骤C构建的像素耕地置信打分规则表，进行耕地数据融合，图1为根据本发明的一个具体实施例的一种耕地图表制作方法的耕地数据初步融合过程的原理示意图；参见图1所示，基于数据预处理和自适应打分表建立后，可得到基于5套耕地产品的耕地分值分布图(a)和耕地平均比例图(b)，利用等面积投影得到每个网格对应的实际面积(c)。从置信度较高的分值1开始，挑选相应分值的像素，计算分值1对应的耕地面积，然后和统计数据做比较，如果小于统计数据，则继续累加置信度较高的2对应的耕地面积，直到和统计数据的耕地面积比较吻合为止，以此确定出该区域的最优分值组合。提取最优分值组合对应的耕地像素，即得到统计数据校正后的耕地分布。因本实施例选用了三级统计数据，因此能得到国家级、省级和地市级三层统计数据校正的耕地融合初步结果。

步骤E，根据步骤D所得到的三级耕地融合初步结果，从低至高逐步融合，得到最终的耕地融合结果。融合流程如下，首先将省级融合结果和地市级融合结果按照规则整合，得到次国家级融合结果，然后按照相同的规则融合次国家和国家级的结果，得到最终的耕地融合结果。

具体融合规则为：

(1)在某个下级行政区域内，根据步骤D的初步融合结果，计算上一级行政区域的地图中对应下级行政区域范围的融合结果的耕地估算面积，以及下一级行政区域的融合结果的耕地估算面积，如果上一级耕地估算面积＜下一级耕地估算面积，则选取下一级耕地融合结果作为上一级行政区域中对应本级行政区划的地图范围的耕地融合结果，，否则，原上一级行政区域中对应本级行政区划的地图范围的融合结果不变。

(2)对于某个下级行政区域，如果没有统计数据，将按照以下公式计算该下级行政区域的统计数据：

剩余区域统计数值(即没有统计数据的所有下级行政区划)＝上一级耕地统计数据-【步骤(1)所获得的各对应下一级耕地融合结果的估算面积之和】，

(3)根据步骤(2)所获得的剩余区域统计数值，利用步骤D的方法对该区域的像素重新进行校正，获得该区域的耕地融合结果，并获得最终的耕地制图数据。

图2为图1的方法中所得到的多级融合结果的整合规则原理示意图，图中整体范围表示国家行政区划，国家内包含的各子区域，即次国家行政区划(省级)用字母a～f表示。假设国家级耕地面积统计数据是2900，结合国家级融合结果图，估算各次国家行政区域内的耕地估算面积。参见次国家级融合结果图，由于该级统计数据的不完整性，使得某些子区域没有融合结果，如子区域a，在两级数据融合时，遵循如下规则：

If 国家级融合面积<次国家融合面积

Then 融合结果＝次国家融合面积

Return 融合数据校正过程

with 统计面积＝国家级耕地面积-较大的次国家级耕地面积

由于子区域e和f的次国家融合面积高于国家级融合面积，因此在最后的两级结果整合图中，这两个区域以次国家级融合为准；对于其它区域a～d，则用国家级耕地面积减去已确定子区域，以此得到的耕地面积为基准，再次进行剩余区域的统计数据校正，最后得到两级融合数据的整合结果。如果在所有子区域，国家级融合面积大于次国家融合面积，则以国家级融合为最终结果。类似的，省级和地市级的融合规则同上所述。

图3为根据图1的方法和图2的规则获得的中国耕地融合结果的耕地示意图；参见图3所示，从耕地比例分布图中，可看出在平原和盆地区域，如华北平原、三江平原和四川盆地等地，地理景观一致性较好，各种产品都能比较准确的提取出耕地，因此这些区域耕地比例较高，置信度较好。在西北和华南的山区和丘陵地带，地形比较复杂，耕地和其它地物的混合比较普遍，这些区域遥感分类具有一定的不确定性，因此耕地比例相对较低，置信度较差。

利用基于检验样本的混淆矩阵对融合之后中国区域的耕地结果进行精度评价。耕地样本1403个，非耕地样本1421个，精度评价矩阵如下表所示。其中耕地和非耕地的错分率分表是21.87％和20.83％，漏分率分别是20.81％和21.89％，融合结果的总体精度是78.65％，Kappa系数是0.57。

从分区和总体两方面，对耕地融合结果和5套耕地产品数据进行精度对比分析。图4为图3的耕地融合结果和五套耕地产品数据的精度对比分析示意图；参见图4所示，融合结果的总体精度高于提供数据的5套耕地产品。其中，GlobeLand30在5套输入数据中总体精度最高，级76.27％，其次是MODIS C5(76.22％)和CCI-LC(74.22％)，GlobCover 2009的耕地数据精度最低，即70.50％。在各个分区中，除了华中和华南区域融合数据的精度略低外，其它区域融合结果精度都高于原始输入数据。由此可见，本发明所提供的融合方法能够较好的兼顾输入数据的优势，提高耕地制图的总体精度和可靠性。

除提高制图精度外，本发明所提供的方法的另一个有益效果是提高了耕地制图和统计数据的一致性。下面从相关性分析和均方根误差(root mean square error，RMSE)两方面，评估5套耕地产品的输入数据和融合数据的耕地数量与统计数据的总体离散度和吻合度。图5对比分析了图3的耕地融合结果和五套耕地产品数据的耕地估算面积和统计数据的均方根误差(RMSE)和相关系数的示意图。RMSE的值越大，表示和统计数据的离散度越高。从图5可看出，本发明所提供的融合方法的结果的数据点大多数分布在1:1线上及其附近，融合结果的RMSE最小，和统计数据的离散度最低；融合结果的相关系数最高，估算的各省耕地面积和统计数据的吻合度最好。因此耕地数量较为准确，与统计数据的一致性较好。

本发明所提供的一种耕地图表制作方法，通过对多套现有技术平台提供的耕地产品进行融合，所获得的耕地分布图总体精度和可靠性大大提升，且耕地数量较为准确，与统计数据的一致性较好。

本领域技术人员应当理解，虽然本发明是按照多个实施例的方式进行描述的，但是并非每个实施例仅包含一个独立的技术方案。说明书中如此叙述仅仅是为了清楚起见，本领域技术人员应当将说明书作为一个整体加以理解，并将各实施例中所涉及的技术方案看作是可以相互组合成不同实施例的方式来理解本发明的保护范围。

以上所述仅为本发明示意性的具体实施方式，并非用以限定本发明的范围。任何本领域的技术人员，在不脱离本发明的构思和原则的前提下所作的等同变化、修改与结合，均应属于本发明保护的范围。

本发明中所引用的参考文献包括：

[1]Ban,Y.,Gong,P.,&Giri,C.Global land cover mapping using Earthobservation satellite data:recent progresses and challenges.ISPRS Journal ofPhotogrammetry and Remote Sensing,2015,103,1–6.

[2]Bontemps,S.；Defourny,P.；Bogaer,t E.V.；Arino,O.；Kalogirou,V.；Perez,J.R.GLOBCOVER 2009.Products Description and Validation Report.Availableonline:https://core.ac.uk/download/pdf/11773712.pdf(accessed on 30April2017).

[3]Chen,J.,Chen,J.,Liao,A.,Cao,X.,Chen,L.,Chen,X.,He,C.,Han,G.,Peng,S.,Lu,M.,Zhang,W.,Tong,X.,&Mills,J.Global land cover mapping at 30mresolution:A POK-based operational approach.ISPRS Journal of Photogrammetryand Remote Sensing,2015,103,7-27.

[4]Chen,J.,Lu,M.,Chen,X.,Chen,J.,Chen,L.A spectral gradientdifference based approach for land cover change detection.ISPRS Journal ofPhotogrammetry and Remote Sensing,2013,85,1-12.

[5]Defourny,P.；Kirches,G.；Brockmann,C.；Boettcher,M.；Peters,M.；Bontemps,S.；Lamarche,C.；Schlerf,M.；Santoro,M.Land Cover CCI:Product UserGuide Version 2.(https://www.esa-landcover-cci.org/).

[6]Fritz,S.；See,L.；Mccallum,I.；You,L.；Bun,A.；Moltchanova,E.；etal.Mapping global cropland and field size.Glob.Chang.Biol.2014,21,1980–1992.

[7]Fritz,S.；See,L.；You,L.；Justice,C.；Becker‐Reshef,I.；Bydekerke,L.；etal.The need for improved maps of global cropland.Eos Trans.AGU 2013,94,31–32.

[8]Loveland,T.R.；Reed,B.C.；Brown,J.F.；Ohlen,D.O.；Zhu,Z.；Yang,L.；etal.Development of a global land cover characteristics database and IGBPDISCover from 1km AVHRR data.Int.J.Remote Sens.2000,21,1303–1330

[9]Schepaschenko,D.；See,L.；Lesiv,M.；Mccallum,I.；Fritz,S.；Salk,C.；etal.Development of a global hybrid forest mask through the synergy of remotesensing,crowdsourcing and FAO statistics.Remote Sens.Environ.2015,162,208–220.

[10]You L,Wood S,Wood-Sichra U,et al.Generating global cropdistribution maps:From census to grid.Agricultural Systems,2014,127:53-60.

Claims (1)

1.一种耕地制图方法，其通过对至少三个不同的耕地产品i进行处理，制作出特定分辨率M的耕地分布图，其包括如下步骤：

步骤A：以国家为单元，估算多个所述耕地产品i在国家j的耕地面积A_i,j，对于每个耕地产品i，采用地理坐标进行数据预处理，将每个耕地产品i处理为所述特定分辨率M的耕地分布图，之后基于等面积投影原理计算每个网格对应的实际面积，然后利用以下公式计算所述耕地产品i中的耕地面积A_i,j：

A_i,j＝∑p_i,j×a_i,j

上式中，p_i,j是预处理后的耕地比例，a_i,j是对应的含有耕地像素的网格的实际面积，

步骤B：对于步骤A获得的多个所述特定分辨率M的耕地分布图中的每个相同坐标的像素，构建不同产品的权重打分表，利用以下公式计算每个耕地产品i的面积差异自适应指数Diff_i,j：

上式中，A_FAO,j是国家j的耕地统计数据，指数Diff_i,j的值越小，表明耕地产品i和统计数据越接近，那么将该耕地产品赋予较高的权重；否则，和统计数据差别越大，给该耕地产品赋予较低的权重，

步骤C:按照以下表格结构，构建耕地置信打分规则表，每个耕地产品所对应取值规则为，该耕地产品认为该像素为耕地时取值为1，否则取值为0，耕地一致量表示对同一像素同时有多少个耕地产品认为该像素为耕地，置信分是从1开始顺序增加的整数，值越大表示置信度越低，耕地一致量相同时，耕地产品权重越高，则置信度越高，相应的置信分越低，

步骤D：根据步骤C构建的像素耕地置信打分规则表，根据行政区域分多级对数据进行处理，在每一级耕地统计数据所对应的行政区划相应的地图范围，从置信度较高的分值1开始，挑选相应分值的像素，计算分值1对应的耕地面积，然后和统计数据做比较，如果小于所述行政区域范围的耕地统计数据，则继续累加分值较高的2对应的耕地面积，直到和统计数据的耕地面积比较吻合为止，以此确定出该区域的耕地融合初步结果，提供多级耕地统计数据，根据该步骤得到多级耕地融合初步结果，

步骤E：按照行政区划级别的从低至高的顺序，逐步整合步骤D的初步融合结果，得到最终的耕地制图结果，整合规则为：

(1)在某个下级行政区域内，根据步骤D的初步融合结果，分别计算上一级行政区域的地图中对应本级行政区域范围的融合结果的耕地估算面积，以及本级行政区域融合结果的耕地估算面积，如果上一级行政区域的耕地估算面积＜本级耕地估算面积，则选取本级行政区域的耕地融合结果作为上一级行政区域中对应本级行政区划的地图范围的耕地融合结果，否则，原上一级行政区域中对应本级行政区划的地图范围的融合结果不变，

(2)对于某个下级行政区域，如果没有统计数据，在上一级行政区划的地图中按照以下公式计算该下级行政区域的统计数据：

剩余区域统计数值＝上一级耕地统计数据-【步骤(1)所获得的各对应下一级耕地融合结果所获得的估算面积之和】，

(3)根据步骤(2)所获得的剩余区域统计数值，利用步骤D的方法对该区域的像素重新进行校正，获得该区域的耕地融合结果，并获得最终的耕地制图数据。

Images