CN106123812A - 基于遥感影像获取起伏地表甘蔗种植面积的方法及装置 - Google Patents

Abstract

本发明公开了一种基于遥感影像获取起伏地表甘蔗种植面积的方法及装置，方法包括：获取甘蔗地块的图斑矢量图及DME栅格数据；根据DME栅格数据获取图斑矢量图的坡度均值；根据坡度均值与坡度阈值将甘蔗地块的图斑矢量图分为平地地块和坡度地块；获取平地地块的面积；将坡度地块进行微面元分解成多个相互邻接的平面小面元，获取平面小面元的面积，将多个平面小面元的面积相加获取坡度地块的面积；将两个地块的面积相加获取甘蔗地块的面积。本发明根据遥感获取的平面范围矢量，结合高精度数字高程模型进行甘蔗地块的地形判断，通过重要地形和边界特征的提取，实现甘蔗真实种植面积的快速拟合计算，为政策制订和生产管理决策提供准确的数据支持。

Description

基于遥感影像获取起伏地表甘蔗种植面积的方法及装置

技术领域

本发明涉及甘蔗种植支撑领域，特别涉及一种基于遥感影像获取起伏地表甘蔗种植面积的方法及装置。

背景技术

甘蔗是我国重要的糖料作物，主要分布于我国广西、云南、海南和广东等南方省份，其中广西的甘蔗种植面积及蔗糖产量长期以来均占全国60％以上，位居全国首位。近年来甘蔗种植出现较大幅度波动，准确掌握甘蔗种植面积和变化状况对于政府决策及糖企生产管理具有重大意义。当前，遥感技术已成为甘蔗种植面积监测的有力手段，然而以往遥感面积监测获取的大多是平面投影面积，在旱坡地广泛分布的甘蔗产区，这种忽略地形起伏变化的面积计算方式可能会导致较大的误差。

现有的遥感甘蔗种植面积信息提取是基于各种分类方法对遥感影像进行分类，尔后对分类图像进行专题信息的矢量化得到。由于卫星遥感影像是地物在二维平面上的映射，因此基于影像分类获得的甘蔗矢量图斑也是其真实种植面积在二维平面上的投影。而直接采用二维投影面积作为实际种植面积，未考虑图斑的地形变化情况，在旱坡地甘蔗种植广泛的甘蔗主产区，这种面积计算方式会导致结果存在较大误差。

在目前遥感甘蔗种植面积提取研究中还未发现考虑地表起伏的甘蔗真实面积获取方法。地形起伏拟合等相关技术多出现在三维建模、三维可视化等计算机技术领域，强调机器视觉展示效果，在作物种植面积相关计算中并不多见。地形起伏拟合可看作是曲面面积计算的一种，目前曲面面积计算一般是基于多项式的函数积分拟合方法，首先利用多项式对特定方向的某一地形断面的曲线形状进行拟合，然后在垂直断面方向对该函数进行积分计算，从而得到曲面面积。基于函数拟合积分的曲面面积方法虽然具有简单方便的特点，但其存在以下缺点：基于函数积分的曲面拟合适用于地形较为固定，变化不大的小区域，当区域较大或地形复杂时，多项式函数难以准确表达多变的地形表面，会给结果带来较大误差；由于DEM格网具有较多的数据冗余，图斑对象内地形剖面线的函数拟合、优化和积分需要耗费大量计算资源和运算时间，当图斑对象数量较大时，计算效率低下。

公开于该背景技术部分的信息仅仅旨在增加对本发明的总体背景的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域一般技术人员所公知的现有技术。

发明内容

本发明的目的在于提供一种基于遥感影像获取起伏地表甘蔗种植面积的方法及装置，从而克服传统作物种植面积计算方式效率低下，结果带来较大误差的缺点。

为实现上述目的，根据本发明一方面，本发明提供了一种基于遥感影像获取起伏地表甘蔗种植面积的方法，包括：

S100：获取甘蔗地块的图斑矢量图以及DME栅格数据；

S101：根据所述DME栅格数据获取所述图斑矢量图的坡度均值；根据所述坡度均值与坡度阈值将所述甘蔗地块的图斑矢量图分为平地地块和坡度地块；

S102：获取所述平地地块的面积；

S103：将所述坡度地块进行微面元分解成多个相互邻接的平面小面元，获取所述平面小面元的面积，将多个所述平面小面元的面积相加以获取所述坡度地块的面积；以及

S104：将所述平地地块和所述坡度地块的面积相加获取所述甘蔗地块的面积。

优选地，上述技术方案中，步骤102具体包括：对所述平地地块形成的多边形进行任意一点与该多边形的每条边构成的多个三角形依次面积求和。

优选地，上述技术方案中，步骤S103具体包括：将所述坡度地块的多个水平面投影按地形特征点分为相互邻接的多个平面小面元；根据所述平面小面元的顶点高程建立三维空间坐标，获取所述平面小面元的面积。

优选地，上述技术方案中，对所述平面小面元内的地形特征点进行三角剖分，然后对所述平面小面元形成的多边形进行三角剖分；获取所述平面小面元进行三角剖分后每个三角形的面积，从而获取所述平面小面元的面积。

为实现上述目的，根据本发明一方面，提供了一种基于遥感影像获取起伏地表甘蔗种植面积的装置，包括：

数据获取模块：获取甘蔗地块的图斑矢量图以及DME栅格数据；

地块分解模块：根据所述DME栅格数据获取所述图斑矢量图的坡度均值；根据所述坡度均值与坡度阈值将所述甘蔗地块的图斑矢量图分为平地地块和坡度地块；

第一计算模块：获取所述平地地块的面积；

第二计算模块：将所述坡度地块进行微面元分解成多个相互邻接的平面小面元，获取所述平面小面元的面积，将多个所述平面小面元的面积相加以获取所述坡度地块的面积；以及

第三计算模块：将所述平地地块和所述坡度地块的面积相加获取所述甘蔗地块的面积。

优选地，上述技术方案中，所述第一计算模块具体包括：对所述平地地块形成的多边形进行任意一点与该多边形的每条边构成的多个三角形依次面积求和。

优选地，上述技术方案中，所述第二计算模块具体包括：将所述坡度地块的多个水平面投影按地形特征点分为相互邻接的多个平面小面元；根据所述平面小面元的顶点高程建立三维空间坐标，获取所述平面小面元的面积。

优选地，上述技术方案中，对所述平面小面元内的地形特征点进行三角剖分，然后对所述平面小面元形成的多边形进行三角剖分；获取所述平面小面元进行三角剖分后每个三角形的面积，从而获取所述平面小面元的面积。

与现有技术相比，本发明具有如下有益效果：

本发明根据遥感获取的平面范围矢量，结合高精度数字高程模型(DEM)进行甘蔗地块的地形判断，并通过重要地形和边界特征的提取，实现甘蔗真实种植面积的快速拟合计算，为政策制订和生产管理决策提供准确的数据支持。

附图说明

图1是根据本发明基于遥感影像获取起伏地表甘蔗种植面积的方法的流程图。

图2是根据本发明多边形面积计算的示意图。

图3是根据本发明坡面地块分解的示意图。

图4-1是根据本发明地形特征提取点的示意图。

图4-2是根据本发明地形特征点的Delaunay三角剖分的示意图。

图4-3是根据本发明剖分后进行边缘处理的示意图。

图4-4是根据本发明地块微面元化的示意图。

图5是根据本发明Delaunay三角网的空圆特性示意图。

图6是根据本发明最大化最小角特性示意图。

图7是根据本发明LOP优化过程示意图。

图8是根据本发明微面元拟合面积计算示意图。

图9是根据本发明三角形从平面投影至三维空间的边长换算示意图。

图10是根据本发明基于遥感影像获取起伏地表甘蔗种植面积的装置的结构图。

具体实施方式

下面结合附图，对本发明的具体实施方式进行详细描述，但应当理解本发明的保护范围并不受具体实施方式的限制。

除非另有其它明确表示，否则在整个说明书和权利要求书中，术语“包括”或其变换如“包含”或“包括有”等等将被理解为包括所陈述的组成部分，而并未排除其它组成部分。

本实施例提出一种顾及地形特征的甘蔗真实种植面积计算方法，根据遥感获取的平面范围矢量，结合高精度数字高程模型(DEM)进行甘蔗地块的地形判断，并通过重要地形和边界特征的提取，实现甘蔗真实种植面积的快速拟合计算。

图1显示了根据本发明优选实施方式的一种基于遥感影像获取起伏地表甘蔗种植面积的方法，包括：

S100：获取甘蔗地块的图斑矢量图以及DME栅格数据；DEM栅格数据为高精度DME栅格数据；

S101：根据DME栅格数据获取图斑矢量图的坡度均值；根据坡度均值与预设的坡度阈值(如坡度阈值Slope>2°)将甘蔗地块的图斑矢量图分为平地地块和坡度地块；

S102：获取平地地块的面积；

S103：将坡度地块进行微面元分解成多个相互邻接的平面小面元，获取平面小面元的面积，将多个平面小面元的面积相加以获取坡度地块的面积；

S104：将平地地块和坡度地块的面积相加获取甘蔗地块的面积。

步骤102具体包括：对平地地块形成的多边形进行任意一点与该多边形的每条边构成的多个三角形依次面积求和。具体为，平地地块面积的计算等价于求平面多边形的面积，本方案选用了O’Rourke提出的任意多边形面积计算方法，即多边形面积等于对任意一点与多边形每条边构成的多个三角形的矢量面积依次求和得出，为便于计算，任意点取多边形外接的某一顶点；而三角形的矢量面积可由两个平面向量的外积求得，如图2所示。

继续参考图2，多边形Ω的顶点P_k(k＝1,2,…,m)沿逆时针排列，P_k+1＝P₁，坐标依次是(x₁,y₁),(x₂,y₂),…,(x_m,y_m)逆时针为正方向，顺时针为负方向。设由向量外积计算得因此任意多边形的面积公式为：

在地块面积计算时，O点取多边形外接矩形的左下角顶点，各顶点坐标为对应的地理投影坐标，单位为米。

步骤S103具体包括：

遥感获取的旱坡地甘蔗图斑地形多样，不仅包括单一坡面，也可能包含山脊、山谷，山峰等复杂地形，因而坡地地块的面积计算需考虑地块内部的地形变化特点。常见的函数曲面拟合需进行大量复杂的积分运算，不仅计算效率低，而且也难以对复杂地形进行准确描述。针对这种情况，该实施例提出了基于地形特征点提取的微面元分解与微面元三维拟合相结合的坡地面积计算方法，其中，复杂曲面的微面元分解是曲面面积计算的关键，主要包括地形特征点提取和包含离散点任意多边形的三角剖分两部分，其技术流程如图3所示，具体包括：

步骤S1031：首先将坡度地块(复杂曲面或多面体表面)的水平面投影按地形特征点拆分为相互邻接的平面小面元；

该步骤中，涉及地形特征点的提取，其具体为：

地形特征点是地形发生凹凸起伏变化的“拐点”，描述了地块图斑内部的地形起伏变化特点，地形特征点之间的区域可近似看作一平面，可用地形特征点作为顶点构建小面元，实现曲面的微面元分解。地形特征点的类型、定义及判别规则如表1所示。根据地形特征点的数学判别规则，以3×3窗口遍历地块范围内的DEM栅格，从8个方向搜寻地形特征点。值得注意的是，当山脊、山谷等线特征表现为一直线时，线上所有点均处于同一平面，此时只需保留两个端点，剔除线段中间的冗余特征点以减少拆分步骤，提高计算效率。因此，在所有初始地形特征提取完成后，重新采用3×3窗口进行迭代以实现直线型地形特征的检测与冗余点剔除。地形特征点简化过程如下：

1)将初始地形特征点在DEM栅格中的位置以链表形式存储，读取链表中第i个特征点位置作为3×3窗口的中心，将窗口定位至当前点(第一次读取时i＝1)，转入步骤2；

2)判断窗口8个方向有无同类型的地形特征点(如同为山脊点或山谷点)，若没有则i＝i+1，返回步骤1读取下一点位置；如某方向上存在同类型特征点，移动窗口并将该点作为新的中心点，执行步骤3；

3)在同一方向检查是否还存在同类型点，如不存在则i＝i+1，返回步骤1；如存在则将当前中心点标识为冗余点，并将窗口移至下一点；

4)重复执行步骤3，直至链表中所有点均得到遍历，最后将标识为冗余的点剔除，完成地形特征点的简化。

由于地块边界矢量大多由遥感分类结果进行矢量化获得，矢量化过程可能引入过多冗余点，因此也需要进行“拐点”的检测和简化。矢量边界简化过程简介如下：

第1步，选定闭环矢量中的头两个节点P₁，P₂为拐点，依次计算余下节点P_k(k＝3，4，...，n)和第2点连线L_PkP2与前两点连线L_P1P2的夹角∠P_kP₂P₁，取序列中∠P_kP₂P₁大于预设角度阈值(该实施例取175°；当三点共线时为180°)的第一个点作为新拐点P_f1，剔除序列中排在拐点之前的其它节点；

第2步，拐点位置后移一位，即P1＝P2，P2＝P_f1，重新计算P_f1之后余下序列中各点与L_P1P2的夹角，获取新拐点及去除多余节点；

第3步，重复执行第2步过程直至序列中节点数为0，完成矢量边界的简化。

最后将矢量边界点集与地形特征点集作比对，剔除地形特征点集中与边界位置重合的点，使得地形特征点始终处于边界内部，同时将该点作为边界拐点添加至边界拐点链表中，避免后继处理中边界为凹多边形时的三角剖分错误和提高三角剖分效率。

表13×3窗口搜索地形特征点

13	12	9
			10	15	9
9	12	13

步骤S1032：根据各平面小面元的顶点高程建立三维空间坐标，对平面小面元内的地形特征点进行三角剖分，然后对平面小面元形成的多边形进行三角剖分；获取平面小面元进行三角剖分后每个三角形的面积，从而获取平面小面元的面积。在最大程度保留原有地形特征的同时，有效提高计算效率和精度。

经简化后的地形特征点和边界矢量如图4-1所示，其平面投影为包含离散点(地形特征点)的多边形。由于三角形在构造和面积计算方面的便捷性，故以三角形作为微面元剖分单元。考虑到凹多边形边界在剖分中的复杂性(边界顶点与地形特征点一起剖分时，可能出现凹多边形的剖分结果为凸多边形的情况)，将三角剖分过程分为针对地形特征点和针对边界顶点两部分进行处理。

首先对多边形内部的离散地形特征点进行Delaunay三角剖分，生成互不相交的三角网(图4-2)；然后对三角网凸轮廓与边界构成的环状区域(图4-3)进行任意多边形的Delaunay三角剖分，最终生成图4-4所示地块图斑的完整三角网结构。以下对离散点和任意多边形的Delaunay剖分作简要介绍。

(1)离散点Delaunay三角剖分

点集的三角剖分须满足三个条件：1.除了端点，平面图中的边不包含点集中的任何点；2.没有相交边；3.所有的面都是三角形，且所有三角面的合集是散点集的凸包；Delaunay三角剖分是三角剖分的特殊情形，除了具备以上3种性质，还符合两项重要准则，一是空圆特性，三角网中任意四点不能共圆，即Delaunay三角网中任一三角形的外接圆范围内不会有其它点存在，如图5所示；二是最大化最小角特性，在两个相邻的三角形构成凸四边形的对角线，在相互交换后，六个内角的最小角不再增大，使得Delaunay三角网是最接近于规则化的三角网，最大程度避免狭长型三角形的出现，如图6所示。

Delaunay三角网具有一系列优异特性：1.最接近，以最近邻的三点构成三角形，且各三角形的边不相交；2.唯一性，不论从区域何处开始构建，最终都将得到一致的结果；3.最优性，任意两个相邻三角形形成的凸四边形的对角线如果可以互换的话，那么两个三角形六个内角中最小的角度不会变大；4.最规则，将三角网中每个三角形按最小角作升序排列，则Delaunay三角网的排列得到的数值最大；5.区域性，新增、删除、移动某一个顶点时只会影响临近的三角形；6.具有凸多边形的外壳；Delaunay三角网最外层边界形成一个凸多边形的外壳。

平面Delaunay三角网的构造有多种成熟算法，如逐点插入法、三角网生长法、分治算法等，由于前两种算法时间复杂度较高，这里采用基于递归的分治算法实现地形特征点的三角剖分，其基本步骤是：

1)把点集V按横、纵坐标升序排列，然后递归地执行以下2-6步骤；

2)将点集V分为近似相等的子集V_L和V_R；

3)在V_L和V_R中生成三角网；

4)用Lawson提出的局部优化算法LOP优化所生成的三角网，是指成为D-三角网；

5)找出连接V_L和V_R中两个凸壳的底线和顶线；

6)由底线至顶线合并V_L和V_R中两个三角网。

通过分治的思路将问题简化，把点集划分到足够小，使其易于生成三角网，然后把子集中的三角网合并生成最终的三角网，用LOP算法保证其成为D-三角网。

Lawson提出的局部优化过程LOP(Local Optimization Procedure)过程步骤如下所示，一般三角网经过LOP处理，即可确保成为Delaunay三角网：

1)将两个具有公共边的三角形合并成一个多边形；

2)以最大空圆准则作检查，看其第四个顶点是否在三角形外接圆内；

3)如果在，修正对角线即将对角线对调，完成局部优化过程的处理。LOP处理过程如下图7所示：

(2)任意多边形Delaunay三角剖分

地形特征点的Delaunay三角剖分完成后，其凸包外壳与甘蔗图斑边界之间形成含岛的多边形，此时内外边界区域的三角剖分实质上等同于对任意多边形(包括含岛的凸/凹多边形)的三角剖分，内环凸壳的边可由是否具有两个邻接三角形的判别准则从边结构链表中进行逐一遍历得出。任意多边形的Delaunay三角剖分有多种成熟的构建方法，这里采用基于有向边的任意多边形Delaunay三角剖分改进算法，其基本思想是：首先判断当前多边形任一边界点p是否在当前边的左侧，若是则计算p点与当前边的夹角，并添加到候选节点列表中，然后用快速排序法对候选节点链表按照夹角从大到小排列，当候选节点链表中与当前边夹角最大，且与当前边可是的候选节点即为能和当前边构成有效的Delaunay三角形的顶点。所以取排序后的第一个候选节点p，判断p点是否与当前边L₁₁L₁₂可视。若可视则p点即为能和当前边构成合法Delaunay三角形的节点，若不可使则取下一个候选节点，找到与当前边可是的候选节点即为所求的顶点。

(3)微面元面积拟合计算

微面元的真实面积通过对平面投影下的三角网进行三维空间坐标换算和面积累加获得，其流程如图8所示。从三角形结构链表中每次取出一个三角形，读取三角形三个顶点的x,y,z坐标值，首先根据x,y坐标计算三角形每条边在投影水平面的边长，然后结合各点的高程z值计算每条边在三维空间的真实长度(几何换算关系如图9所示)，接着利用海伦公式根据已知的三条边边长计算三角形面积：p＝(a+b+c)/2；其中a,b,c分别为三角形边长，p为半周长；最后，通过对链表中所有三角形的三维空间面积进行累计求和，得出地块图斑的真实面积。

本发明根据遥感获取的平面范围矢量，结合高精度数字高程模型(DEM)进行甘蔗地块的地形判断，并通过重要地形和边界特征的提取，实现甘蔗真实种植面积的快速拟合计算，为政策制订和生产管理决策提供准确的数据支持。

图10显示了根据本发明优选实施方式的一种基于遥感影像获取起伏地表甘蔗种植面积的装置，

数据获取模块10：获取甘蔗地块的图斑矢量图以及DME栅格数据；

地块分解模块20：根据DME栅格数据获取图斑矢量图的坡度均值；根据坡度均值与坡度阈值将甘蔗地块的图斑矢量图分为平地地块和坡度地块；

第一计算模块30：获取平地地块的面积；

第二计算模块40：将坡度地块进行微面元分解成多个相互邻接的平面小面元，获取平面小面元的面积，将多个平面小面元的面积相加以获取坡度地块的面积；以及

第三计算模块50：将平地地块和坡度地块的面积相加获取所述甘蔗地块的面积。

第一计算模块30具体包括：对平地地块形成的多边形进行任意一点与该多边形的每条边构成的多个三角形依次面积求和。

第二计算模块40具体包括：将坡度地块的多个水平面投影按地形特征点分为相互邻接的多个平面小面元；根据平面小面元的顶点高程建立三维空间坐标，获取平面小面元的面积。更为具体地，对平面小面元内的地形特征点进行三角剖分，然后对平面小面元形成的多边形进行三角剖分；获取平面小面元进行三角剖分后每个三角形位于三维空间坐标中的面积，从而获取平面小面元的面积。

本发明根据遥感获取的平面范围矢量，结合高精度数字高程模型(DEM)进行甘蔗地块的地形判断，并通过重要地形和边界特征的提取，实现甘蔗真实种植面积的快速拟合计算，为政策制订和生产管理决策提供准确的数据支持。

前述对本发明的具体示例性实施方案的描述是为了说明和例证的目的。这些描述并非想将本发明限定为所公开的精确形式，并且很显然，根据上述教导，可以进行很多改变和变化。对示例性实施例进行选择和描述的目的在于解释本发明的特定原理及其实际应用，从而使得本领域的技术人员能够实现并利用本发明的各种不同的示例性实施方案以及各种不同的选择和改变。本发明的范围意在由权利要求书及其等同形式所限定。

Claims (8)

1.一种基于遥感影像获取起伏地表甘蔗种植面积的方法，其特征在于，包括：

S100：获取甘蔗地块的图斑矢量图以及DME栅格数据；

S101：根据所述DME栅格数据获取所述图斑矢量图的坡度均值；根据所述坡度均值与坡度阈值将所述甘蔗地块的图斑矢量图分为平地地块和坡度地块；

S102：获取所述平地地块的面积；

S103：将所述坡度地块进行微面元分解成多个相互邻接的平面小面元，获取所述平面小面元的面积，将多个所述平面小面元的面积相加以获取所述坡度地块的面积；以及

S104：将所述平地地块和所述坡度地块的面积相加获取所述甘蔗地块的面积。

2.根据权利要求1所述的基于遥感影像获取起伏地表甘蔗种植面积的方法，其特征在于，步骤102具体包括：

对所述平地地块形成的多边形进行任意一点与该多边形的每条边构成的多个三角形依次面积求和。

3.根据权利要求1所述的基于遥感影像获取起伏地表甘蔗种植面积的方法，其特征在于，步骤S103具体包括：

将所述坡度地块的多个水平面投影按地形特征点分为相互邻接的多个平面小面元；

根据所述平面小面元的顶点高程建立三维空间坐标，获取所述平面小面元的面积。

4.根据权利要求3所述的基于遥感影像获取起伏地表甘蔗种植面积的方法，其特征在于，

对所述平面小面元内的地形特征点进行三角剖分，然后对所述平面小面元形成的多边形进行三角剖分；

获取所述平面小面元进行三角剖分后每个三角形的面积，从而获取所述平面小面元的面积。

5.一种基于遥感影像获取起伏地表甘蔗种植面积的装置，其特征在于，包括：

数据获取模块：获取甘蔗地块的图斑矢量图以及DME栅格数据；

地块分解模块：根据所述DME栅格数据获取所述图斑矢量图的坡度均值；根据所述坡度均值与坡度阈值将所述甘蔗地块的图斑矢量图分为平地地块和坡度地块；

第一计算模块：获取所述平地地块的面积；

第二计算模块：将所述坡度地块进行微面元分解成多个相互邻接的平面小面元，获取所述平面小面元的面积，将多个所述平面小面元的面积相加以获取所述坡度地块的面积；以及

第三计算模块：将所述平地地块和所述坡度地块的面积相加获取所述甘蔗地块的面积。

6.根据权利要求5所述的基于遥感影像获取起伏地表甘蔗种植面积的装置，其特征在于，所述第一计算模块具体包括：

对所述平地地块形成的多边形进行任意一点与该多边形的每条边构成的多个三角形依次面积求和。

7.根据权利要求5所述的基于遥感影像获取起伏地表甘蔗种植面积的装置，其特征在于，所述第二计算模块具体包括：

将所述坡度地块的多个水平面投影按地形特征点分为相互邻接的多个平面小面元；

根据所述平面小面元的顶点高程建立三维空间坐标，获取所述平面小面元的面积。

8.根据权利要求7所述的基于遥感影像获取起伏地表甘蔗种植面积的装置，其特征在于，

对所述平面小面元内的地形特征点进行三角剖分，然后对所述平面小面元形成的多边形进行三角剖分；

获取所述平面小面元进行三角剖分后每个三角形的面积，从而获取所述平面小面元的面积。