CN104574449B

CN104574449B - 基于dem的投影面积计算方法

Info

Publication number: CN104574449B
Application number: CN201510043819.0A
Authority: CN
Inventors: 郭春喜; 聂建亮; 王文利; 蒋光伟
Original assignee: Center For Geodetic Data Processing National Administration Of Surveying Mapping And Geoinformation
Current assignee: Center For Geodetic Data Processing National Administration Of Surveying Mapping And Geoinformation
Priority date: 2015-01-27
Filing date: 2015-01-27
Publication date: 2017-11-21
Anticipated expiration: 2035-01-27
Also published as: CN104574449A

Abstract

本发明提供了一种基于DEM的投影面积计算方法，步骤为：提出DEM数据及封闭区域边界，计算交点坐标；交点数据预处理：在得到交点坐标基础上，需要对交点数据预处理，判断交点的排列顺序，计算交点所在的行列数；边界区域的有效面积计算，按区间依次计算，中间整格网面积计算按照交点在每一行面积，与传统的技术相比，计算简单，且通过计算任意多边形的投影面反推该图斑的实际地表面积，该方法大大缩短计算时间，提高计算效率，且计算精度与构TIN方法相当。

Description

基于DEM的投影面积计算方法

技术领域

本发明涉及土地测量技术领域，特别涉及一种基于DEM的投影面积计算方法。

背景技术

随着计算机的普及与GIS的迅速发展，图斑面积成为GIS中的一项重要统计内容。2013年国务院颁布了第一次地理国情普查的通知，其中一项重要内容就是土地面积的统计计算。在面积计算过程中，必须考虑计算机的内存问题。一幅DEM数据文件包括大量格网数据，若将数据全部读入将造成计算机空间的浪费；常用的面积计算方法构造三角网，该计算方法模型严密，计算精度高，但计算复杂、计算工作量大。而基于DEM的投影面积计算方法根据投影面上图斑面积，利用实际面积与投影图斑面积的比例因子，反算图斑时间面积，从而大大降低计算量，提高计算效率。

发明内容

针对现有技术中的上述缺陷和问题，本发明实施例的目的是提供一种更好的基于DEM的投影面积计算方法，大大缩短计算时间，且计算精度与构TIN方法相当，可以大大提高计算效率。

为了达到上述目的，本发明实施例提供如下技术方案：

一种基于DEM的投影面积计算方法，其特征在于，步骤如下：

Q1：提取计算区域的DEM数据及封闭区域边界，并计算交点坐标；

Q2：交点数据预处理：在得到交点数据基础上，需要对交点数据预处理，判断交点的排列顺序，计算交点所在的行列数；

Q21：交点的排列顺序，计算区域中最北、最东、最南、最西点坐标，统计这四个点在交点文件中的排列位置，然后根据排列位置判断交点的排列顺序，顺序为顺时针或逆时针。

Q22：交点的行列数计算：采用逐行计算区域面积，需要计算交点在整幅图中的行列号，统计每一个格网内的交点个数；每个交点统计的信息包括：交点在交点文件所属于的区域；交点在DEM中的行号；交点在DEM所在的列号；交点在交点文件中的行号；交点在本行交点序列中的位置；交点的属性，0为与格网的交点，1为格网内部点；

Q3：边界区域的有效面积计算，按区间依次计算，

Q31：当多条折线有公共相交部分时，需要判断格网的几个顶点属于多边形，并将格网顶点与交点按一定重新顺序排列，交点的排列顺序为顺时针时，符合左手法则，四指从线段起点指向终点，大拇指垂直于其它四指，大拇指指向为区域内方向)；反之，符合右手法则；计算公式为：

(n为封闭区域交点总个数，i为n个交点中第i个交点)

Q32：当多条折线没有公共相交部分时，将当前格网分割成面积更小的格网，仍采用步骤Q31中的计算公式来计算格网中的不规则多边形的面积；

Q4：中间整格网面积计算：

按照Q22步骤中得到的交点的计算在每一行面积，

Q41：先判断属于区域有效整格网的位置；

Q42：对于每一个整格网，用中心点为顶点的四个三角形的面积代替四边形面积，三角形面积计算采用海伦公式

其中，a、b、c为三角形的边长；S_Δ为三角形面积；S_四边形为四边形面积。

作为上述技术方案的优选，所述Q21步骤中对排列位置判断交点的排列顺序，顺序为顺时针或逆时针，当若最北、最东、最南、最西点分别为1，20，50，88，交点排列顺序为顺时针；当最北、最东、最南、最西点分别为20，1，88，50，交点排列顺序为逆时针。

作为上述技术方案的优选，所述Q22步骤中还需要根据交点所在格网的连续性将交点分成若干区间，当第i行包含N个连续区间，则将该行交点分为N个区间，N为大于1的自然数。

作为上述技术方案的优选，所述Q3步骤应用于地表面积计算时，当对于高分辨的格网，且地形变化幅度不大时，用多边形的投影面积推算其地表面积，计算公式与步骤Q31中的计算公式相同，其投影面积与地表面积之间的函数关系为：其中，S_地表为地表面积，S为使用步骤Q31中的计算公式得出的多边形的投影面积，S_ABCD，S_A1B1C1D1为在高分辨率的格网内，投影面积与地表面积分别所在的四边形的面积。

作为上述技术方案的优选，所述Q41步骤判断属于区域有效整格网的位置的具体方法为，利用点与多边形的所属关系即点在多边形内或外判断方法，然后再利用计算射线与多边形的交点个数，判断点与多边形的所属关系。

本发明实施例提供的一种基于DEM的投影面积计算方法，与传统的技术相比，通过计算任意多边形的投影面反推该图斑的实际地表面积，该方法大大缩短计算时间，提高计算效率，且计算精度与构TIN方法相当。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例1的基于DEM的投影面积计算方法的流程步骤示意图。

图2为本发明实施例1的基于DEM的投影面积计算方法的被拆分成网格的移动窗口的结构示意图。

图3为本发明实施例1的基于DEM的投影面积计算方法的Q21步骤的交点逆时针排列的结构示意图。

图4为本发明实施例1的基于DEM的投影面积计算方法的Q21步骤的交点顺时针排列的结构示意图。

图5为本发明实施例1的基于DEM的投影面积计算方法的Q31步骤的多边形顶点逆时针排列的结构示意图。

图6为本发明实施例1的基于DEM的投影面积计算方法的Q31步骤的多边形顶点顺时针排列的结构示意图。

图7为本发明实施例1的基于DEM的投影面积计算方法的Q32步骤的当多条折线没有公共相交部分时，多条折线组成的多边形的结构示意图。

图8为本发明实施例1的基于DEM的投影面积计算方法的Q3步骤应用于地表面积计算时，多边形投影面积与地表面积的结构示意图。

具体实施方式

下面将结合本发明的附图，对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

如图1-8所示，本发明实施例提供的一种基于DEM的投影面积计算方法，如下：

交点数据预处理

在得到根据DEM文件内插得到的交点数据基础上，需要对交点数据预处理，判断交点的排列顺序，计算交点所在的行列数。

一)交点的排列顺序

计算区域中最北、最东、最南、最西点坐标，统计这四个点在交点文件中的排列位置，然后根据排列位置判断交点的排列顺序，顺时针或逆时针，如图3～4所示，，如若最北、最东、最南、最西点分别为1，20，50，88，交点排列顺序为顺时针；若最北、最东、最南、最西点分别为20，1，88，50，交点排列顺序为逆时针。

二)交点的行列数计算

若采用逐行计算区域面积，需要计算交点在整幅图中的行列号，统计每一个格网内的交点个数。交点数据预处理文件格式为

其中，第1列为某交点在交点文件所属于的区域；第2列为该交点在DEM中的行号；第3列为该交点在DEM所在的列号；第4列为该交点在交点文件中的行号；第5列为该交点在本行交点序列中的位置；第6列为该交点的属性，0为与格网的交点，1为格网内部点。

另外，需要根据交点所在格网的连续性将交点分成若干区间。如图2所示，第i行包含两个连续区间，因此，将该行交点分为两个区间。

边界区域的有效面积计算

边界区域的有效面积计算是面积计算的核心，同时也是面积计算的难点。在计算边界区域有效面积时，按区间依次计算。以第i行第j交点区间第k格网为例，说明边界区域有效面积计算。

一)多条折线有公共相交部分

若第k格网中有n个交点，每一个交点首尾相连，如图5～6所示。在计算不规则多边形面积时，采用

不规则多边形计算公式要求顶点坐标按顺序排列，因此，需要判断格网的几个顶点属于多边形，并将格网顶点与交点按一定重新顺序排列。

为了提高排序的可靠性，在第k格网中的交点为已排好序列的基础上，选取一个与格网相交的交点i和属于多边形的格网顶点进行排序，这样大大降低了多边形顶点重新排序造成的误判率。如图5、6中，将折线由4个交点构成的折线看作为端点为i与i+3线段，按一定法则判断格网4个顶点与线段的关系，最终判定一个顶点属于多边形。线段与格网顶点的关系满足一定的规律，但这种判断需先判断交点的排列顺序，即交点的排列顺序为顺时针时，符合左手法则(四指从线段起点指向终点，大拇指垂直于其它四指，大拇指指向为区域内方向)；反之，符合右手法则。

二)多条折线没有公共相交部分

对于由3条以上折线没有公共相交部分时，采用以上方法较确定有效区域。如图7所示，三条折线没有公共相交部分，格网内有效面积由两个多变形构成，这种图形结构采用以上方法容易造成格网顶点的误判，导致面积计算错误。因此，提出了基于更高分辨率的计算方法，将格网分成更小分辨率的格网，仍采用公式(1)计算格网中的不规则多边形的面积，这样就避免了一个格网有两个以上多边形的情况。

三)地表面积计算

对于高分辨的格网，若地形变化幅度不大，可以用多边形的投影面积推算其地表面积，但如何确定二者可靠合理的函数关系是比较复杂的。在多次比较分析的基础上，提出了一种基于比例因子的新面积计算方法。

如图8所示，在高分辨率的格网内，四边形A₁B₁C₁D₁与四边形ABCD的面积比与不规则多边形S_地表与S的面积比一致，即

若采用公式(1)得到不规则多边形S，则S_地表为

中间整格网面积计算

在交点数据预处理中，已经交点进行区间划分。在每一行面积计算时，要判断属于区域有效整格网的位置。此处采用根据一种组合方法判断整格网的在该行位置，该方法利用点与线段的关系判断方法(线段与格网顶点判断原理一致)，然后再利用陈瑞卿在西安交通大学学报的41卷第1期提出的一种快速算法进一步判断。另外，对于每一个整格网，用中心点为顶点的四个三角形的面积代替四边形面积。三角形面积计算采用海伦公式

实际应用

以陕西省吴堡县数据为例，分析该面积量算方法的计算效果。该数据包括26704个数据文件，约1300万交点个数，最终构成了17710个图斑，共50多km2。

采用新的面积计算方法与构TIN计算方法两种方案，26704个数据文件(17710个图斑)所用计算时间相差较大，其中新的面积计算方法用了3077.329秒，构TIN算法共用了12600秒；总面积互差如下表所示。

两方案总面积互差(m2)	万分比
		31392.5045	0.62

新的面积计算方法根据图斑在投影面上的面积，利用自适应因子推算图斑的地表面积，其精度与严密构TIN方法得到的面积精度相当，计算效率大大提高。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims

1.一种基于DEM的投影面积计算方法，其特征在于，步骤如下：

Q21：交点的排列顺序,计算区域中最北、最东、最南、最西点坐标，统计这四个点在交点文件中的排列位置，然后根据排列位置判断交点的排列顺序，顺序为顺时针或逆时针；

Q22：交点的行列数计算:采用逐行计算区域面积，需要计算交点在整幅图中的行列号，统计每一个格网内的交点个数；每个交点统计的信息包括：交点在交点文件所属于的区域；交点在DEM中的行号；交点在DEM所在的列号；交点在交点文件中的行号；交点在本行交点序列中的位置；交点的属性，0为与格网的交点，1为格网内部点；

Q3：边界区域的有效面积计算,按区间依次计算，

Q31：当多条折线有公共相交部分时，需要判断格网的几个顶点属于多边形，并将格网顶点与交点按一定重新顺序排列,交点的排列顺序为顺时针时，符合左手法则,四指从线段起点指向终点，大拇指垂直于其它四指，大拇指指向为区域内方向-；反之，符合右手法则；计算公式为：

n为封闭区域交点总个数，i为n个交点中第i个交点，S为投影面积；

Q4：中间整格网面积计算：

按照Q22步骤中得到的交点的计算在每一行面积，

Q41：先判断属于区域有效整格网的位置；

Q42：对于每一个整格网，用中心点为顶点的四个三角形的面积代替四边形面积,三角形面积计算采用海伦公式

2.根据权利要求1所述的一种基于DEM的投影面积计算方法，其特征在于，所述Q21步骤中对排列位置判断交点的排列顺序，顺序为顺时针或逆时针，当若最北、最东、最南、最西点分别为1，20，50，88，交点排列顺序为顺时针；当最北、最东、最南、最西点分别为20，1，88，50，交点排列顺序为逆时针。

3.根据权利要求1所述的一种基于DEM的投影面积计算方法，其特征在于，所述Q22步骤中还需要根据交点所在格网的连续性将交点分成若干区间，当第i行包含N个连续区间，则将该行交点分为N个区间，N为大于1的自然数。

4.根据权利要求1所述的一种基于DEM的投影面积计算方法，其特征在于，所述Q3步骤应用于地表面积计算时，当对于高分辨的格网，且地形变化幅度不大时，用多边形的投影面积推算其地表面积，计算公式与步骤Q31中的计算公式相同，其投影面积与地表面积之间的函数关系为：其中,S_地表为地表面积，S为使用步骤Q31中的计算公式得出的多边形的投影面积，S_ABCD,S_A1B1C1D1为在高分辨率的格网内，投影面积与地表面积分别所在的四边形的面积。

5.根据权利要求1所述的一种基于DEM的投影面积计算方法，其特征在于，所述Q41步骤判断属于区域有效整格网的位置的具体方法为，利用点与多边形的所属关系即点在多边形内或外判断方法，然后再利用计算射线与多边形的交点个数，判断点与多边形的所属关系。