基于3S技术的坡耕地划分方法与流程

基于3s技术的坡耕地划分方法
技术领域
1.本发明涉及水土保持行业坡耕地划分领域，具体地涉及基于3s技术的坡耕地划分方法。


背景技术：

2.耕地是最重要的土地资源，耕地的数量和质量决定了一个国家的人口承载力和实现可持续发展的能力。由于经济人口快速增长以及城镇化建设，导致原有耕地被侵占，生态退耕与农业结构调整也加剧了耕地减少的趋势，因此山丘区百姓对耕地需求强烈，转向坡耕地开发。由于人口密集，山地丘陵区面积占比较大，坡耕地的面积分布较为广泛，坡耕地水土流失成为了经济发展中一个突出的生态环境问题。坡耕地是水土流失的主要策源地，长期的水土流失使土地资源遭到破坏，土壤肥力逐年下降，土层减薄，土壤质地变粗，涵养水源和生态保护功能减弱，给生态环境和水利安全带来不利影响，同时制约着农民收入水平的提高和生活质量的改善，也损害了区域社会经济的可持续发展。坡度是反应耕地地表形态的自然要素特征，是评价耕地等级、利用条件和水土流失状况的重要指标，也是制定耕地利用与保护、生态退耕政策的主要依据，准确掌握坡耕地的数量及其空间分布，在制定农业发展战略、实施国土资源整治和开发以及生态建设等方面都具有重要意义。新时期水利改革发展重点工作任务之一就是全面强化水生态治理修复，包括推进坡耕地综合整治，因此，急需通过坡耕地调查，准确掌握坡耕地的数量、质量及其空间分布和土地利用等状况，为坡耕地规划和坡耕地综合整治工作提供准确的治理模式、投入机制及对策。
3.现有技术中，对于坡耕地划分的主流技术还是采用的表格统计方式，即人工采集、统计数据再进行坡耕地划分。虽然现在有各种遥感技术来进行地块的划分，但坡耕地由于其特殊性，使得遥感技术在这个领域的应用一直迟滞不前；其根本原因在于坡耕地不是平面，其坡度的大小、相邻地块之间的关系犬牙交错；另一方面不是所有的坡度都适合耕作；因此直到现在，对于坡耕地还是采用的人工统计划分的技术手段。现有技术的缺陷如下：
4.1.由于没有图件成果，更缺乏落实到图斑的矢量数据，从而导致无法为坡耕地综合整治和坡耕地综合防治规划提供详实的基础数据；
5.2.由于人工采集数据、分析数据并最终落实划分，需要很长时间，从而导致工作效率低下，且数据更新频率很低；
6.3.由于人工作业不可避免会有各种影响数据准确性的因素存在，从而使得最后的划分结果并非最合理的结果。


技术实现要素：

7.本发明针对上述问题，提供基于3s技术的坡耕地划分方法，其目的在于以非常准确的掌握坡耕地的数量、质量及其空间分布和土地利用等状况，为坡耕地规划和坡耕地综合整治工作提供准确的基础数据和科学的技术支撑；脱离人工作业，极大缩短工作時長，极大增高效率；使得最后的划分结果在理论上可以做到最优解。
8.为解决上述问题，本发明提供的技术方案为：
9.一种基于3s技术的坡耕地划分方法，包含以下步骤：
10.s100.采集需要进行坡耕地划分操作的调查区域的基本数据；所述基础数据包含dem、国土二调数据和dom影像；所述国土二调数据包含二调土地利用图；
11.s200.将所述dem和所述dom影像分别按所述调查区域的dem栅格主比例尺进行分幅操作，得到dem拼接图；所述dem拼接图包含栅格；所述栅格的比例尺与调查区域的dem栅格主比例尺相同；
12.s300.根据所述dem拼接图生成坡度分级栅格图；所述坡度分级栅格图包含每个所述栅格的坡度值；
13.s400.将所述坡度分级栅格图转换为矢量坡度分级图；
14.s500.从所述二调土地利用图中提取地类图斑；所述地类图斑包含耕地图斑和园地图斑；
15.s600.根据所述地类图斑和所述矢量坡度分级图制作得到耕地坡度分级分布图；所述耕地坡度分级分布图即为所述坡耕地划分方法的最终结果。
16.优选地，s100与s200之间还包含以下步骤：
17.s150.根据人工预设的数据质量标准，对所述基础数据进行数据质量检查操作，筛选出不符合所述数据质量标准的基础数据；然后根据人工预设的数据完整性标准，对所述基础数据进行数据完整性检查操作，筛选出不符合所述数据完整性标准的基础数据；
18.s151.对筛选出来的不符合数据质量标准的基础数据和不符合数据完整性标准的基础数据进行数据修正操作，使不符合数据质量标准的基础数据达到数据质量标准，且使不符合数据完整性标准的基础数据达到数据完整性标准；
19.s160.根据人工预设的调查要求，对dem的数据格式、dem的坐标系统、 dem的高程基准、dem的投影带、dom影像的数据格式、dom影像的坐标系统、dom影像的高程基准、dom影像的投影带、国土二调数据的数据格式、国土二调数据的坐标系统、国土二调数据的高程基准和国土二调数据的投影带进行调查要求检查操作，筛选出不符合所述调查要求的dem的数据格式、dem 的坐标系统、dem的高程基准、dem的投影带、dom影像的数据格式、dom 影像的坐标系统、dom影像的高程基准、dom影像的投影带、国土二调数据的数据格式、国土二调数据的坐标系统、国土二调数据的高程基准和国土二调数据的投影带；
20.s161.对筛选出来的不符合调查要求的dem的数据格式、dem的坐标系统、 dem的高程基准、dem的投影带、dom影像的数据格式、dom影像的坐标系统、dom影像的高程基准、dom影像的投影带、国土二调数据的数据格式、国土二调数据的坐标系统、国土二调数据的高程基准和国土二调数据的投影带进行数据转换操作，使不符合调查要求的dem的数据格式、dem的坐标系统、 dem的高程基准、dem的投影带、dom影像的数据格式、dom影像的坐标系统、dom影像的高程基准、dom影像的投影带、国土二调数据的数据格式、国土二调数据的坐标系统、国土二调数据的高程基准和国土二调数据的投影带达到调查要求。
21.优选地，s300中所述根据所述dem拼接图生成坡度分级栅格图具体包含以下步骤：
22.s310.根据所述dem拼接图生成坡度图；
23.s320.根据所述坡度图、国土二调的分级标准和水土保持最新要求，生成坡度分级图；所述坡度分级图中的每个栅格都具有1个坡度级；
24.s330.对所述坡度分级图进行栅格优化操作，生成优化后的坡度分级栅格图；所述栅格优化操作中使用的比例尺与调查区域的dem栅格主比例尺相同。
25.优选地，s400中所述将所述优化后的坡度分级栅格图转换为矢量坡度分级图具体包含以下步骤：
26.s410.对所述优化后的坡度分级栅格图进行矢量化操作，生成坡度分级矢量化数据；
27.s420.对所述坡度分级矢量化数据进行优化操作，得到优化后的坡度分级矢量化数据。
28.优选地，s600中所述根据所述地类图斑和所述矢量坡度分级图制作得到耕地坡度分级分布图，具体包含以下步骤：
29.s601.将所述地类图斑与所述矢量坡度分级图叠加，得到每个耕地图斑的多个坡度级分布；
30.s610.统计得到每个耕地图斑不同坡度级下的面积；
31.s620.提取面积最大的耕地图斑的坡度级；
32.s630.将提取到的面积最大的耕地图斑的坡度级赋予所有耕地图斑，得到所述耕地坡度分级分布图。
33.优选地，s330中所述对所述坡度分级图进行栅格化操作，生成所述坡度分级栅格图，具体包含以下步骤：
34.s331.移除所述坡度分级图中孤立的像元；
35.s332.对每个栅格的边缘进行平滑处理操作，使每个栅格的每个边缘都规整平滑；
36.s333.对每个进行过所述平滑处理操作的栅格进行区域分组操作，得到由多个区域组成的坡度分级图；每个所述区域具有唯一的分区值；每个所述区域包含多个像元分组；每个所述像元分组包含像元；
37.s334.将所有像元的数量小于人工预设的像元数阈值的像元分组提取出来；
38.s335.统计每个像元的数量小于人工预设的像元数阈值的像元分组的坡度级；然后以每个像元的数量小于人工预设的像元数阈值的像元分组作为掩膜，与每个像元的数量小于人工预设的像元数阈值的像元分组所对应的坡度级叠加，得到每个像元的数量小于人工预设的像元数阈值的像元分组的合并修改值；
39.s336.重复s334～s335，直至所述坡度分级图不再发生变化，即将不再发生变化的坡度分级图作为所述坡度分级栅格图输出。
40.优选地，s420中所述对所述坡度分级矢量化数据进行优化操作，得到优化后的坡度分级矢量化数据，具体包含以下步骤：
41.s421.消除面积小于人工预设的图斑面积阈值的地类图斑；
42.s422.对所述坡度分级矢量化数据所对应的矢量坡度分级图中的面轮廓中的尖角进行平滑处理操作。
43.优选地，所述坡度值按下式表达：
44.45.其中，slopedegrees为中心像元的坡度，单位为度；为中心像元开始在水平方向的变化率，为中心像元开始在垂直方向上的变化率；atan为可返回数字的反正切值。
46.优选地，所述坡度级包含≤2
°
、2
°
～6
°
、6
°
～15
°
、15
°
～25
°
和≥25
°
。
47.本发明与现有技术对比，具有以下优点：
48.1.由于本发明采用了3s技术进行坡耕地划分，落实了耕地图斑的坡耕地矢量图，从而可以非常准确的掌握坡耕地的数量、质量及其空间分布和土地利用等状况，为坡耕地规划和坡耕地综合整治工作提供准确的基础数据和科学的技术支撑；
49.2.由于本发明采用了遥感、高清成像技术，辅以计算机分析，脱离了人工采集数据、分析数据、划分耕地，从而极大缩短了工作時長，效率极大增高，且数据更新频率也得到了极大加快，甚至可以做到实时更新；
50.3.由于本发明无需人工作业，避免了有各种人为影响数据准确性的因素，从而使得最后的划分结果在理论上可以做到最优解。
附图说明
51.图1本发明具体实施例的基于3s技术的坡耕地划分的流程示意图；
52.图2本发明具体实施例的dem 3
×
3局部移动窗口的示意图；
53.图3本发明具体实施例的利用arcgis进行坡度分析的示意图；
54.图4本发明具体实施例的利用arcgis生成坡度图的示意图；
55.图5本发明具体实施例的众数滤波处理的示意图；
56.图6本发明具体实施例的应用了“众数滤波”处理后的栅格前后对比示意图；
57.图7本发明具体实施例的进行boundary clean的示意图；
58.图8本发明具体实施例的应用了“边界清理”处理后的栅格前后对比的示意图；
59.图9本发明具体实施例的进行区域合并识别聚类的示意图；
60.图10本发明具体实施例的应用了“区域分组”处理后的栅格前后对比示意图；
61.图11本发明具体实施例的在区域分组后为每个区域(断开连接的区域)分配一个唯一的分区值的示意图；
62.图12本发明具体实施例的用con函数选择小于阈值分组并赋值“0”示意图；
63.图13本发明具体实施例的小于阈值部分赋值“0”后的栅格的示意图；
64.图14本发明具体实施例的提取被con函数赋值为0的值的示意图；
65.图15本发明具体实施例的提取小于阈值的区域的示意图；
66.图16本发明具体实施例的进行焦点统计功能计算小于阈值的区域的坡度级数值的示意图；
67.图17本发明具体实施例的按坡度级就低合并至图斑的像元值的示意图；
68.图18本发明具体实施例的按坡度级就低合并后的栅格图像的示意图；
69.图19本发明具体实施例的模型构建器进行栅格综合的批量处理的示意图；
70.图20本发明具体实施例的矢量化后的坡度分级栅格图的示意图；
71.图21本发明具体实施例的对发明实施例中选择面积小于3000m2的地类图斑的示意图；
72.图22本发明具体实施例的对矢量化后的数据进行细小图斑消除的示意图；
73.图23本发明具体实施例的对矢量化后的数据进行界面平滑处理操作的示意图；
74.图24本发明具体实施例的平滑后的矢量坡度分级图的示意图；
75.图25本发明具体实施例的从二调土地利用图中提取耕地和园地图斑的示意图；
76.图26本发明具体实施例的使用tabulateintersection进行坡度提取的示意图；
77.图27本发明具体实施例的进行tabulateintersection的提取结果的示意图；
78.图28本发明具体实施例的使用summarystatistics进行地类图斑内面积最大的坡度提取的示意图；
79.图29为汇总本发明具体实施例中每个地类图斑内的面积最大值的示意图；
80.图30为选择本发明具体实施例地类图斑内分布面积最大的数据的示意图；
81.图31为获取本发明具体实施例地类图斑内分布面积最大的坡度级的示意图；
82.图32为本发明具体实施例耕地坡度分级分布图的示意图。
具体实施方式
83.下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本技术所附权利要求所限定的范围。
84.所謂3s技术是遥感技术(remotesensing，rs)、地理信息系统(geographyinformationsystems，gis)和全球定位系统(globalpositioningsystems，gps)的统称，是空间技术、传感器技术、卫星定位与导航技术和计算机技术、通讯技术相结合，多学科高度集成的对空间信息进行采集、处理、管理、分析、表达、传播和应用的现代信息技术。
85.如图1所示，一种基于3s技术的坡耕地划分方法，包含以下步骤：
86.s100.采集需要进行坡耕地划分操作的调查区域的基本数据；基础数据包含dem、国土二调数据和dom影像；国土二调数据包含二调土地利用图。
87.需要说明的是，调查区域的基本数据还包含基础资料和相关资料；其中基础资料包含社会经济数据、气象水文数据、地质地貌数据、土壤植被数据、土地利用数据、水土流失数据、水土保持数据；相关资料包含环境保护数据、国土数据、林业数据、建设数据、农业数据。
88.本具体实施例中，dom影像的精度为2.0m。
89.s150.根据人工预设的数据质量标准，对基础数据进行数据质量检查操作，筛选出不符合数据质量标准的基础数据；然后根据人工预设的数据完整性标准，对基础数据进行数据完整性检查操作，筛选出不符合数据完整性标准的基础数据。
90.s151.对筛选出来的不符合数据质量标准的基础数据和不符合数据完整性标准的基础数据进行数据修正操作，使不符合数据质量标准的基础数据达到数据质量标准，且使不符合数据完整性标准的基础数据达到数据完整性标准。
91.s160.根据人工预设的调查要求，对dem的数据格式、dem的坐标系统、dem的高程基准、dem的投影带、dom影像的数据格式、dom影像的坐标系统、dom影像的高程基准、dom影像的投影带、国土二调数据的数据格式、国土二调数据的坐标系统、国土二调数据的高程基准和国土二调数据的投影带进行调查要求检查操作，筛选出不符合调查要求的dem的数据
格式、dem的坐标系统、dem的高程基准、dem的投影带、dom影像的数据格式、dom 影像的坐标系统、dom影像的高程基准、dom影像的投影带、国土二调数据的数据格式、国土二调数据的坐标系统、国土二调数据的高程基准和国土二调数据的投影带。
92.s161.对筛选出来的不符合调查要求的dem的数据格式、dem的坐标系统、 dem的高程基准、dem的投影带、dom影像的数据格式、dom影像的坐标系统、dom影像的高程基准、dom影像的投影带、国土二调数据的数据格式、国土二调数据的坐标系统、国土二调数据的高程基准和国土二调数据的投影带进行数据转换操作，使不符合调查要求的dem的数据格式、dem的坐标系统、 dem的高程基准、dem的投影带、dom影像的数据格式、dom影像的坐标系统、dom影像的高程基准、dom影像的投影带、国土二调数据的数据格式、国土二调数据的坐标系统、国土二调数据的高程基准和国土二调数据的投影带达到调查要求。
93.需要说明的是，s150～s161的作用在于对基础数据，尤其是其中的dem、 dom影像和国土二调数据的质量、完整性等进行检查，将不符合调查要求的数据格式、坐标系统、高程基准、投影带等进行转换；这样做的意义在于前期格式化操作，使所有要使用的数据在比例尺、数据格式、坐标系统、高程基准、投影带等各处处于同样的格式和标准，这样才能进入后面的操作。
94.需要进一步说明的是，实际工作中，为提高运算速度，可将栅格面积较大的区域按照下一级行政区划范围进行影像的分割处理，进行坡耕地的提取。
95.s200.将dem和dom影像分别按调查区域的dem栅格主比例尺进行分幅操作，得到栅格化dem和栅格化dom影像；dem栅格主比例尺由人工预设；然后将栅格化dem和栅格化dom影像拼接起来，得到dem拼接图；dem拼接图包含栅格；栅格的比例尺与调查区域的dem栅格主比例尺相同。
96.s300.根据dem拼接图生成坡度分级栅格图；坡度分级栅格图包含每个栅格的坡度值。
97.需要说明的是，坡度值取决于表面从中心像元开始在水平方向和垂直方向上的变化率(增量)。
98.本具体实施例中，坡度值按式(1)表达：
[0099][0100]
其中，slopedegrees为中心像元的坡度，单位为度；为中心像元开始在水平方向的变化率，为中心像元开始在垂直方向上的变化率；atan为可返回数字的反正切值。
[0101]
需要说明的是，利用式(1)计算出每个栅格的坡度值，生成坡度分级栅格图。
[0102]
如图2所示，计算时采用dem 3
×
3窗口利用坡度计算模型计算坡度值。通过典型地貌样区(山区、丘陵区、平原区)不同算法提取的坡度统计特征对比分析，采取三阶反距离平方权差分算法提取的坡度较为准确，符合实际，本算法也是arcgis坡度分析中默认采取的算法。三阶反距离平方权差分算法按式(2) 表达：
[0103][0104]
其中，e
i
(i＝1，2，
…
，8)分别表示中心点e周围栅格点的高程；g为相邻网格点的水平和垂直距离。
[0105]
本具体实施例中，s300具体包含以下步骤：
[0106]
s310.根据dem拼接图生成坡度图。
[0107]
如图3所示，本具体实施例中，生成坡度图包含以下步骤：
[0108]
s311.用arcmap加载dem数据，使用空间分析中的slope功能：arctoolboxs —spatial analyst tools—surface—slope。
[0109]
s312.输入栅格地形数据，选择输出路径，输出单位默认degree，得到坡度图。
[0110]
s320.根据坡度图、国土二调的分级标准和水土保持最新要求，生成坡度分级图；坡度分级图中的每个栅格都具有1个坡度级。
[0111]
本具体实施例中，坡度级将土地分为了5个等级，包含≤2
°
、2
°
～6
°
、6
°
～ 15
°
、15
°
～25
°
和≥25
°
。
[0112]
如图4所示，本具体实施例中，生成坡度分级图包含以下步骤：
[0113]
s321.使用空间分析中的slope功能，输入栅格地形数据，将分类级别改为5。
[0114]
s322.依次输入2、6、15、25、90。
[0115]
s323.确认后选择输出路径，得到坡度分级图；
[0116]
s330.对坡度分级图进行栅格优化操作，生成优化后的坡度分级栅格图；栅格优化操作中使用的比例尺与调查区域的dem栅格主比例尺相同；这一步的作用在于。在坡度分级图的生成过程中经常会生成许多孤立的小型数据区域，另外在不同坡度类别的边界上通常会生成无法真实表示边界的锯齿状边界，可利用栅格综合工具对其进行平滑处理，生成最终的分类的栅格。栅格综合工具常见的工具有边界清理、众数滤波、蚕食、收缩、细化等工具。
[0117]
本具体实施例中，具体包含以下步骤：
[0118]
s331.移除坡度分级图中孤立的像元。
[0119]
本具体实施例中，使用“众数滤波”工具移除孤立的像元；具体来说：
[0120]
要移除已分类影像中单个的孤立的像元，可应用众数滤波工具。
[0121]
众数滤波工具可根据像元邻域内的众数值来替换像元。此工具需要满足两个条件才能发生替换。首先，相同值的邻近像元的数量必须多到可以成为众数值，或者至少一半的像元必须具有相同值(视指定的参数而定)。并且这些像元在滤波器内核周围必须是连续的。即，如果指定的是众数参数(majority)，则四分之三或八分之五的已连接像元必须具有相同的值；如果指定的是半数参数 (half)，则需要四分之二或八分之四的已连接像元具有相同的值。其次，那些像元必须与指定的滤波器的中心相邻。如果将替换阈值设置为half，并且两个值作为相等部分出现，则当处理的像元值与半数中某一像元值相同时将不会发生替换。half选项比majority选项允许的过滤范围广泛。第二个条件与像元的空间连通性有关，目的是将像元的空间模式的破坏程度降到最低。使用四个相邻像元数(four)会保留矩
形区域的拐角。使用八个相邻像元(eight)将使矩形区域的拐角变得平滑。将相邻像元数设置为八时，相邻的定义是共享一条边，角像元在所有相邻像元均具有相同值时才能发生更改，而边像元需要三个相邻像元(包括边上的像元)具有相同值才发生更改。将相邻像元数设置为四时，相邻的定义是共享一个角，边或角像元始终要求存在两个匹配的相邻像元才能发生替换。如果不满足这些条件，将不会进行替换，像元的值也将保持不变。运行几次众数滤波后，输出的栅格将会稳定下来(不再变化)。
[0122]
如图5和图6所示，在本次划分中，使用空间分析中的major filter功能，输入栅格坡度数据，经过反复几次众数滤波应用之后，许多较小的像元组就已经消失了。
[0123]
s332.对每个栅格的边缘进行平滑处理操作，使每个栅格的每个边缘都规整平滑。
[0124]
本具体实施例中，使用“边界清理”工具对区域进行平滑处理操作；具体来说：
[0125]
要平滑区域之间的边界，可使用边界清理工具。边界清理工具主要用于清理区域间不规整的边缘。该工具使用扩展和收缩的方法在相对较大的范围上清理边界。最初，优先级较高的区域在各个方向上覆盖其邻近的优先级较低的区域，覆盖大小为一个像元。然后，它们收缩回至那些没有完全被相同值的像元包围的像元。除内部像元以外的任何像元均可替换。此外，还可以替换区域中的小岛屿(可被视为与区域共用边界)。可保留的最小区域为3
×
3的像元块。因此，可能替换狭窄的区域。
[0126]
对于不按大小进行排序(no_sort)的默认方法，较大的值具有较高的优先级。descend—以大小的降序顺序对区域进行排序。总面积较大的区域具有较高的优先级，可以扩展到总面积较小的若干区域。ascend—以大小的升序顺序对区域进行排序。总面积较小的区域具有较高的优先级，可以扩展到总面积较大的若干区域。
[0127]
如图7和图8所示，在本次划分中，使用空间分析中的boundary clean功能，输入栅格数据，通过放大和缩小边界，较大区域将会侵占较小区域，更多较小、较窄的像元组就已经消失了。
[0128]
s333.对每个进行过平滑处理操作的栅格进行区域分组操作，得到由多个区域组成的坡度分级图；每个区域具有唯一的分区值；每个区域包含多个像元分组；每个像元分组包含像元。
[0129]
需要说明的是，这个步骤s333的作用在于识别聚类，将无分类的像元进行重新识别、分配。
[0130]
本具体实施例中，使用“区域合并”工具进行区域分组操作；具体来说：
[0131]
主滤波和边界清理工具仅可通过为误分类的像元分配直接邻域中出现频率最高的值的方式来清理单个像元或由少量误分类的像元组成的极小型像元聚类。但假设存在这样一个特定大小的阈值，该阈值以下的各相似像元分组将被视为过小，而不会对后续分析具有任何意义。同时，这些聚类应融入到周围组中。在本项目中，对属于同种坡度级别的毗连聚类而言，如果其图上像元数量小于3个，则将被视为不具备分析意义。然而，无法对这些孤立区域单独进行处理，这是由于这些区域具有与整个区域相同的坡度级值。
[0132]
要解决此问题，可使用区域合并工具。此工具会为输入栅格(已分类的影像)中的每个区域分配唯一标识符。区域是指具有相同值的一组连续像元。假设某单一区域(zone)由两个不相连的区域(region)组成。区域合并工具将把该区域分为两个新区域，各区域都具有唯一标识(区域)值。原始区域值将作为 link字段保存在输出属性表中。
[0133]
分区(zone)由栅格中所有具有相同值的像元组成。区域(region)是相同分区 (zone)类型的一组连续像元。分区(zone)可以由若干个不相连的区域(region)组成。需要分别处理区域(region)时，必须将每个区域识别为独立的实体。区域分组工具用于为栅格中的每个区域指定新值。通过扫描过程来指定值，扫描时从栅格的左上角开始，然后从左向右移动，再从上向下移动。遇到新区域时，为其指定一个唯一的值。此过程将继续执行到所有区域都分配到了一个值为止。系统将会为每个区域分配唯一编号。扫描的第一个区域接收值为1，第二个区域接收值为2，依此类推，直到所有区域均已赋值。扫描将按从左至右、从上至下的顺序进行。分配给输出分区的值取决于系统对其进行扫描的时间。
[0134]
默认情况下，为输出数据增加链接字段(python中的add_link)选项已启用。这将在输出栅格的属性表中创建名为link的项，其保留输入栅格的每个像元的原始值。link字段用于追踪每个新创建的查询或分析区域的来源。
[0135]
输出时，包含排除值的像元位置会接收零，这些区域便不会与现有nodata 像元位置相混淆。由于区域分组的编号从值1开始，因此被排除在重新分组操作之外的像元将被视为背景。这些背景像元可被重新分类或处理为任何其他值。使用条件函数工具，可将包含排除值的位置轻松地转换为nodata。
[0136]
number_neighbors在评估像元间的连接时使用的相邻像元数。four仅当具有相同值的像元与上下左右四个最邻近像元中的每一个像元直接连接时，才会定义这些像元之间的连通性。如果两个具有相同值的像元彼此只是对角线连接，则其不会被视为相连接。这是默认设置。eight仅当具有相同值的像元位于彼此最近的8像元邻域内(8个最邻近像元)时，才会定义这些像元间的连通性。其中包括彼此之间的上下左右或对角线连接。
[0137]
如图9和图10所示，使用空间分析中的region group功能，输入栅格数据。
[0138]
如图11所示，在区域分组之后，为每个区域(断开连接的区域)分配一个唯一的分区值。其输出栅格的属性表如下，属性表中link项为原地表坡度级， count表示该分组内的像元数量。
[0139]
s334.将所有像元的数量小于人工预设的像元数阈值的像元分组提取出来。
[0140]
如图13所示，本具体实施例中，使用条件函数(con函数)，选择分组内像元数小于等于3的分组。
[0141]
如图12所示，将其value值赋值为0，其余值不变。
[0142]
条件函数根据像元值在指定的条件语句中被评估为“真”还是“假”，条件函数工具允许控制每个像元的输出值。如果像元被判定为“真”，它将获得一类值；如果像元被判定为“假”，它将获得另一类值。当像元被判定为“真”时，它所获得的输出值由输入条件为真时所取的栅格数据或常数值指定。当像元被判定为“假”时，它所获得的输出值由输入条件为假时所取的栅格数据或常数值指定。
[0143]
如图14和图15所示，使用extract by attributes提取刚刚被con函数赋值为 0的值，也就是像元数小于3个的分组。
[0144]
s335.统计每个像元的数量小于人工预设的像元数阈值的像元分组的坡度级；然后以每个像元的数量小于人工预设的像元数阈值的像元分组作为掩膜，与每个像元的数量小于人工预设的像元数阈值的像元分组所对应的坡度级叠加，得到每个像元的数量小于人工预设的像元数阈值的像元分组的合并修改值。
[0145]
如图16所示，本具体实施例中，采用焦点统计功能计算小于像元数阈值的区域的坡度级的数值。
[0146]
本具体实施例采取的栅格综合规则之一为，将选中的像元数量小于3个的分组的图斑按坡度级就低不就高原则并入邻近图斑。因此需计算小于阈值的像元位置的合并值，即该位置邻域中的最小值。焦点统计工具可执行用于计算输出栅格数据的邻域运算，各输出像元的值是其指定邻域范围内所有输入像元值的函数。运算该函数可得出统计数据，例如最大值、平均值或者邻域内所有值的总和。从概念上讲，此算法在执行过程中将访问栅格中的每个像元，并且根据识别出的邻域范围计算出指定的统计数据。要计算统计数据的像元称为待处理像元。待处理像元的值以及所识别出的邻域中的所有像元值都将包含在邻域统计数据的计算中。各邻域可以重叠，因此某一邻域中的像元也可以包含在其他待处理像元的邻域中。邻域可以是环形(圆环)、圆形、矩形或楔形，通过使用核文件，也可自定义邻域形状，以及在计算统计数据之前将不同的权重分配给邻域中的各个特定像元。邻域内可以计算的统计量有均值、众数、最大值、中值、最小值、少数、范围、标准差、总和以及变异度。在计算中忽略nodata 选项可控制邻域窗口内nodata像元的处理方式。选中此选项时(data选项)，输出像元值的计算将会忽略邻域中的所有nodata像元。
[0147]
如图16所示，在进行焦点统计之前，需对要采取计算的像元点赋空值，以免对邻域分析结果产生影响，采用赋空值计算；对赋空值之后的栅格图采取focalstatistics，输入栅格数据，neighborhood选择irregular(不规则)，kernel file选取txt文本，文本中输入{3 3；1 1 1；1 0 1；1 1 1}行列式，表示周围3
×
3邻域的统计值权重。statistics type选择minimum(最小值)。
[0148]
如图17所示，选择前面步骤s334提取的需要修改合并到周围的像元分组作为掩膜，与焦点统计之后的值叠加，最后得到局部需要修改合并后的值。
[0149]
s336.重复s334～s335，直至坡度分级图不再发生变化，即将不再发生变化的坡度分级图作为坡度分级栅格图输出。步骤s336的作用在于反复迭代，得到稳定的输出结果。
[0150]
本具体实施例中，步骤s336中再次采用con函数进行运算，将如图16所示的修改后了像元值的栅格与如图13所示的提取小于阈值产生空值的栅格进行合并。
[0151]
如图18所示，通过操作，小于阈值的像元组已经消失。使用由区域合并工具生成的结果中的链接项，可将已分类影像的原始区域值重新分配给通过区域合并工具创建的各个区域。经过几次迭代运算后，栅格值趋向于稳定，最终形成的土地的坡度分级栅格图。
[0152]
需要进一步说明的是，如图19所示，在对栅格进行批量操作时，可选择模型建构器(model builder)进行处理。模型构建器是一个用来创建、编辑和管理模型的应用程序。模型是将一系列地理处理工具串联在一起的工作流，它将其中一个工具的输出作为另一个工具的输入。也可以将模型构建器看成是用于构建工作流的可视化编程语言。模型构建器除了有助于构造和执行简单工作流外，还能通过创建模型并将其共享为工具来提供扩展arcgis功能的高级方法。采用模型构建器可以大大减轻工作的负担，提高arcgis的工作效率。
[0153]
s400.将优化后的坡度分级栅格图转换为矢量坡度分级图；这一步的作用在于将坡度分级栅格图进行矢量化处理，形成坡度分级矢量数据，即矢量坡度分级图。
[0154]
步骤s400具体包含以下步骤：
[0155]
s410.对优化后的坡度分级栅格图进行矢量化操作，生成坡度分级矢量化数据。
[0156]
本具体实施例中，采用raster to polygon进行栅格转矢量计算；如图20所示为转换后的图。
[0157]
s420.对坡度分级矢量化数据进行优化操作，得到优化后的坡度分级矢量化数据；这一步骤的作用在于进行图的优化，涉及但不限于对坡度分级矢量化数据进行图斑综合处理、界线平滑处理、拓扑重建处理、数据裁切处理。
[0158]
步骤s420具体包含以下步骤：
[0159]
s421.消除面积小于人工预设的图斑面积阈值的地类图斑；这一步的作用在于消除细小的地类图斑。
[0160]
步骤s421具体包含以下步骤：
[0161]
s421.1将面积小于阈值的图斑合并至周边图斑中。经过栅格制图综合后，仍有少量小图斑残留，需要进一步的处理。使用消除工具通过将面与具有最大面积或最长公用边界的邻近面合并来消除面。要消除的要素由应用于面图层的选择内容决定。必须在之前的步骤中使用按属性选择图层或按位置选择图层或者通过查询地图图层来确定选择内容。
[0162]
s421.2如图21所示，使用select by attributes选择面积小于3000m2的图斑。
[0163]
s421.3如图22所示，然后将面积小于阈值的图斑合并。
[0164]
s422.对坡度分级矢量化数据所对应的矢量坡度分级图中的面轮廓中的尖角进行平滑处理操作。步骤s422的作用在于对坡度矢量图中的面轮廓中的尖角进行平滑处理以使制图更加美观。
[0165]
需要进一步说明的是，arcgis提供了两种平滑方法：指数核的多项式近似 (paek)方法(python中的paek)可根据平滑容差对面进行平滑处理。每个面经过平滑处理后，其折点都可能比之前多。平滑容差参数可控制计算新折点时用到的“移动”路径的长度。长度越短，保留的细节越多，处理时间也越长。贝塞尔插值方法(python中的bezier_interpolation)对面进行平滑处理时无需使用容差，而应通过创建近似的贝塞尔曲线来匹配输入面。
[0166]
本具体实施例中，采取多项式近似(paek)方法进行面的平滑，smoothingtolerance填写参数经多次调整后；在应用过程中，优选200，因为这时平滑效果较好。
[0167]
如图23所示，生成的局部矢量化图。
[0168]
如图24所示，为平滑后的本具体实施例矢量坡度分级图。
[0169]
s500.从二调土地利用图中提取地类图斑；地类图斑包含耕地图斑和园地图斑；步骤s500的作用在于从二调土地利用图(dltb)中提取耕地和园地图斑，作为坡度赋值对象；具体来说，如图25所示，使用属性表里的按属性选择功能，选择地类编码为011(水田)、012(水浇地)、013(旱地)的字段导出。
[0170]
s600.根据地类图斑和矢量坡度分级图制作得到耕地坡度分级分布图；耕地坡度分级分布图即为坡耕地划分方法的最终结果；具体包含以下步骤：
[0171]
s601.将地类图斑与矢量坡度分级图叠加，得到每个耕地图斑的多个坡度级分布；步骤s601还包含确定地类图斑的坡度级，并在地类图斑的层中对耕地图斑赋耕地坡度级属性，为后续步骤做准备。
[0172]
s610.统计得到每个耕地图斑不同坡度级下的面积；步骤s610使用tabulate
intersection(交集制表)进行坡度提取。
[0173]
需要说明的是，如图26所示，为提高运算速度，采取交集制表功能对坡度图和耕地图进行叠加，交集制表的目的是统计某块耕地图斑内各坡度级别的面积，采用交集制表可不用直接分割耕地图斑，仅以表格的形式记录统计结果。
[0174]
如图27所示，对每个耕地图斑和园地图斑(bsm)内各个坡度级别(gridcode) 的面积(shape_area)进行汇总。
[0175]
s620.提取面积最大的耕地图斑的坡度级；步骤s620使用summary statistics (汇总统计数据)进行地类图斑内面积最大的耕地图斑的坡度级提取；具体来说：
[0176]
在交集制表过程中，获得了每个图斑内的不同坡度级的面积，根据规则，当调查的耕地图斑涉及两个以上坡度级时，面积最大的坡度级为该耕地图斑的坡度级；因此还需获取单个图斑内的面积最大的坡度级，采用summary statistics (汇总统计数据)进行图斑内面积最大的坡度提取。summary statistics输出表将由包含统计运算结果的字段组成。
[0177]
如图28和图29所示，该工具可用于以下统计运算：总和、平均值、最大值、最小值、范围、标准差、计数、第一个和最后一个。
[0178]
s630.将提取到的面积最大的耕地图斑的坡度级赋予所有耕地图斑，得到耕地坡度分级分布图；步骤s630具体包含以下步骤：
[0179]
s631.如图30所示，建立耕地坡度汇总与交集制表的关联，使用join功能，使用标识码(bsm)作为关联要素，然后在属性表中按属性选择。
[0180]
s632.如图31所示，选择地类面积中与坡度汇总面积最大面积相同的数据，导出表格，再将导出的表格与地类图斑表格建立关联(join)。
[0181]
如图32所示，经过s100到s632，最终可以得到各地类图斑的坡度等级，得到的耕地坡度分级分布图即为本发明的坡耕地划分方法的最终结果。
[0182]
需要说明的是，到s632为止，得到的耕地坡度分级分布图已经是非常准确的，可以直接使用了；但耕地相关的问题在我国属于红线，因此必须经过实地检查、复核工作的检验，才能最终作为决策参考，报送决策部门及领导使用。
[0183]
因此，本具体实施例还包含以下步骤：
[0184]
s700对s600得到的耕地坡度分级分布图进行现场复核；具体包含以下步骤：
[0185]
s710现场复核工作底图和调查表制作；具体来说：
[0186]
选择坡耕地面积占国土面积比例较大以及陡坡耕地占耕地总面积比例较大的区域作为典型区域进行调查，兼顾区域地形地貌、气候、土地利用不同分布情况，选取代表性强的区域作为调查区域，每个调查区域根据不同坡度级(≤2
°
、2
°
～6
°
、6
°
～15
°
、15
°
～25
°
和≥25
°
)；在每个坡度级选择3～5个面积较大的坡耕地图斑，作为调查点。调查点在各个调查县范围内均匀分布，且处于交通较为便利的位置，以免现场复核时无法到达。按照有关制图规范，将遥感影像、划分后的坡度分级矢量图及其他辅助性的地理信息图层，如乡镇界、道路、居民点等图层叠加，制作调查区域和调查点现场复核工作底图。调查县工作底图打印a1大小，便于全局观看，了解和掌握整个县市坡耕地分布详情；调查点工作底图打印a4大小，为调查图斑的局部放大，方便现场对照核查。
[0187]
制作坡耕地调查表，主要填写内容包括经纬度、土地利用类型、坡度、耕作层厚度、种植作物类型、土壤侵蚀状况、水保措施等信息。用a4纸打印多份；需要注意的是，不能采用
网络手段填写调查表，以倒逼调查人员进入实地复核。如表1所示，为坡耕地调查表的例表，实际工作中可以有多种表现形式，但内容不能少于表1：
[0188]
表1.坡耕地调查表
[0189][0190]
s720现场调查复核；具体来说：
[0191]
根据室内确定的调查区域和调查点，进行坡耕地图斑实地调查，核实图斑实际土地利用情况，运用坡度仪进行坡度量测，填写坡耕地调查表，采用无人机和摄像机拍摄现场照片和视频作为佐证材料等，照片包括全景照片和局部坡度量测照片。
[0192]
s730复核结果整理；具体来说：
[0193]
将现场调查的结果进行归档整理，包括表格和照片，与调查点底图图斑属性进行对应，判别土地利用和坡度解译的准确性。
[0194]
本具体实施例中，现场复核共调查了大悟县、南漳县、罗田县、通山县、郧西县、秭归县、咸丰县等共10余个县市，共计135个调查点，经过复核验证，耕地原始坡度与本发明提取的坡度较为一致，准确率达到95.7％。
[0195]
上述数据表明本发明的方法在坡耕地划分方面具有很强的适用性和准确性。
[0196]
在上述的详细描述中，各种特征一起组合在单个的实施方案中，以简化本公开。不应该将这种公开方法解释为反映了这样的意图，即，所要求保护的主题的实施方案需要比清楚地在每个权利要求中所陈述的特征更多的特征。相反，如所附的权利要求书所反映的那样，本发明处于比所公开的单个实施方案的全部特征少的状态。因此，所附的权利要求书特此清楚地被并入详细描述中，其中每项权利要求独自作为本发明单独的优选实施方案。
[0197]
为使本领域内的任何技术人员能够实现或者使用本发明，上面对所公开实施例进行了描述。对于本领域技术人员来说；这些实施例的各种修改方式都是显而易见的，并且本文定义的一般原理也可以在不脱离本公开的精神和保护范围的基础上适用于其它实施例。因此，本公开并不限于本文给出的实施例，而是与本技术公开的原理和新颖性特征的最广范围相一致。
[0198]
上文的描述包括一个或多个实施例的举例。当然，为了描述上述实施例而描述部件或方法的所有可能的结合是不可能的，但是本领域普通技术人员应该认识到，各个实施例可以做进一步的组合和排列。因此，本文中描述的实施例旨在涵盖落入所附权利要求书的保护范围内的所有这样的改变、修改和变型。此外，就说明书或权利要求书中使用的术语“包含”，该词的涵盖方式类似于术语“包括”，就如同“包括，”在权利要求中用作衔接词所解释的那样。此外，使用在权利要求书的说明书中的任何一个术语“或”是要表示“非排它性的或者”。
[0199]
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含
在本发明的保护范围之内。
[0200]
最后，应当指出，以上实施例仅是本发明较有代表性的例子。显然，本发明不限于上述实施例，还可以有许多变形。凡是依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均应认为属于本发明的保护范围。
[0201]
在此，需要说明的是，上述技术方案的描述是示例性的，本说明书可以以不同形式来体现，并且不应被解释为限于本文阐述的技术方案。相反，提供这些说明将使得本发明公开将是彻底和完整的，并且将向本领域技术人员充分传达本说明书所公开的范围。此外，本发明的技术方案仅由权利要求的范围限定。
[0202]
用于描述本说明书和权利要求的各方面公开的形状、尺寸、比率、角度和数字仅仅是示例，因此，本说明书和权利要求的不限于所示出的细节。在以下描述中，当相关的已知功能或配置的详细描述被确定为不必要地模糊本说明书和权利要求的重点时，将省略详细描述。
[0203]
在使用本说明书中描述的“包括”、“具有”和“包含”的情况下，除非使用否则还可以具有另一部分或其他部分，所用的术语通常可以是单数但也可以表示复数形式。
[0204]
应该指出，尽管在本说明书可能出现并使用术语“第一”、“第二”、“顶部”、“底部”、“一侧”、“另一侧”、“一端”、“另一端”等来描述各种不同的组件，但是这些成分和部分不应受这些术语的限制。这些术语仅用于区分一个成分和部分和另一个成分和部分。例如，在不脱离本说明书的范围的情况下，第一部件可以被称为第二部件，并且类似地，第二部件可以被称为第一部件，顶部和底部的部件在一定情况下，也可以彼此对调或转换；一端和另一端的部件可以彼此性能相同或者不同。
[0205]
此外，在构成部件时，尽管没有其明确的描述，但可以理解必然包括一定的误差区域。
[0206]
在描述位置关系时，例如，当位置顺序被描述为“在...上”、“在...上方”、“在... 下方”和“下一个”时，除非使用“恰好”或“直接”这样的词汇或术语，此外则可以包括它们之间不接触或者接触的情形。如果提到第一元件位于第二元件“上”，则并不意味着在图中第一元件必须位于第二元件的上方。所述部件的上部和下部会根据观察的角度和定向的改变而改变。因此，在附图中或在实际构造中，如果涉及了第一元件位于第二元件“上”的情况可以包括第一元件位于第二元件“下方”的情况以及第一元件位于第二元件“上方”的情况。在描述时间关系时，除非使用“恰好”或“直接”，否则在描述“之后”、“后续”、“随后”和“之前”时，可以包括步骤之间并不连续的情况。本发明的各种实施方案的特征可以部分地或全部地彼此组合或者拼接，并且可以如本领域技术人员可以充分理解的以各种不同地构造来执行。本发明的实施方案可以彼此独立地执行，或者可以以相互依赖的关系一起执行。