CN106682237B - 一种山体边界自动提取方法 - Google Patents

Abstract

本发明提供了一种山体边界自动提取方法，获取山体边界面、山脊线和山谷线数据；基于山脊线数据，关联山脊线经过的面，初次形成山体边界；查找与山体边界具有邻接关系的面，并以山谷线为约束，关联符合条件的面；将规则约束条件外的面进行关联，最终形成地理尺度的山体边界。本发明基于山体边界面数据，以地形特征线作为约束条件，关联符合条件的面，实现了地理尺度山体边界的自动提取，提取后的山体较为完整，尺度相对合理，可直接满足地理研究的需要，因此可成为构建山体对象的基础，并为地貌学和人文学的特征分析、地质灾害后的山体重构以及地形多尺度综合提供了有效支持。

Description

一种山体边界自动提取方法

技术领域

本发明涉及地理信息技术领域，具体涉及一种山体边界自动提取方法。

背景技术

目前，在已存的地形模型中，通常采用规则采样点、三角面片、等高线等形式来表达山体，而避免直接表达山体对象。实际上，地形中存在山峰与河谷等现象，也能明确其位置和形态，可见地形中存在对象。在地貌学中，山体作为一种重要的地形对象，其构建主要有以下必要性：(1)为地貌学中地貌特征及水文学中水文特征的分析提供服务；(2)在地质灾害方面，能够在山体重构时做出迅速地改变和调整；(3)山体的对象模型，可以为地形多尺度综合提供特征数据。

传统的山体边界自动提取方法主要思路为：(1)数据预处理：将原始DEM利用差值运算进行地形倒置，形成反地形DEM，再对反地形DEM进行填洼，获取无洼地的DEM；(2)水流方向提取：对反地形无洼地的DEM进行流向分析；(3)流域分析：基于流向数据，对其进行流域分析以实现盆地边界的提取，最后经栅矢转换得到矢量流域边界，初步获取山体边界。

由于山地地形的复杂性及多尺度问题，导致上述方法在实施时，提取的山体边界较为破碎，与地理视角下的山体边界不一致，难以满足地理研究需要，需进一步手工处理，投入大，效率低。因此，对破碎的山体边界基于规则约束进行关联，最终形成地理尺度的山体对象极为关键。

发明内容

发明目的：本发明的目的在于针对现有技术的问题，基于地形特征线作为地形的骨架线来控制地形的起伏变化，提出了一种山体边界自动提取方法，通过相对破碎的山体边界面数据，以地形特征线作为约束条件，关联符合条件的面，基本实现了地理尺度山体边界的自动提取。

技术方案：本发明提供了一种山体边界自动提取方法，包括以下步骤：

(1)获取山体边界面、山脊线和山谷线数据；

(2)基于山脊线数据，关联山脊线经过的面，初次形成山体边界；

(3)查找与山体边界具有邻接关系的面，并以山谷线为约束，关联符合条件的面；

(4)将规则约束条件外的面进行关联，最终形成地理尺度的山体边界。

进一步，步骤(1)通过加载DEM数据，运用ArcGIS并基于地形表面水文分析原理获取矢量的图层数据，针对相对破碎的山体边界面图层，运用ArcGIS中的消除工具，以字段Area阈值，消除小于该阈值的图斑，不断调整阈值的大小，经多次消除，得到山体边界面要素图层。

进一步，步骤(2)包括：

(2-1)加载山体边界面要素图层、山脊线线要素图层和山谷线线要素图层，分别得到山体边界面要素集合Mount＝{mt_i|i＝1,2,…,m}、山脊线线要素集合RidgeL＝{rl_i|i＝1,2,…,n}和山谷线线要素集合ValleyL＝{vl_i|i＝1,2,…,r}；

其中，mt_i表示第i个面要素，m为面要素的个数；rl_i表示第i个山脊线线要素，n为山脊线线要素的个数；vl_i表示第i个山谷线线要素，r为山谷线线要素的个数；

(2-2)遵循山脊线线要素个数n决定地理尺度的山体边界面要素个数的原则，对集合Mount中的每一个面要素mt_i按以下步骤进行属性添加与赋值：

①给面要素mt_i添加属性Id，Id表示山脊线编号，初始化为-1；

②读取集合RidgeL中的线要素rl_i包含的点要素集合RLPt＝{rp_j|j＝1,2,…,a}，rp_j表示线要素rl_i上第j个点要素，a为点要素的个数；

③判断集合RLPt中的点rp_j是否在面要素mt_i内，当存在任一点符合条件时，修改mt_i的编号属性Id为rl_i的下标i；

④循环执行步骤②-③，直至集合RidgeL中的线要素遍历完毕；

(2-3)根据编号属性Id的不同，对面要素mt_i分组：遍历集合Mount中的每一个面要素mt_i，当mt_i的Id等于-1时，将mt_i加入面要素集合NoUnionM；当Id等于1时，将mt_i加入面要素集合UnionS₁，当Id等于2时，将mt_i加入面要素集合UnionS₂，……，当Id等于n时，将mt_i加入面要素集合UnionS_p；遍历完毕，得到集合NoUnionM＝{nm_i|i＝1,2,…,q}和以集合为元素的集合SET＝{UnionS₁,UnionS₂,…,UnionS_p}；

其中，nm_i表示不符合关联条件的面要素，q为面要素的个数，UnionS_p表示具有相同Id的面要素集合，p为面要素集合的个数；

(2-4)面要素关联：对SET中的每一个面要素集合UnionS_p，进行面要素合并，形成新的面要素，得到关联后的面要素集合UnionM＝{um_i|i＝1,2,…,g}，um_i表示关联后的面要素，g为面要素的个数。

进一步，步骤(3)包括：

(3-1)给集合NoUnionM中的一面要素nm_i添加属性Fid、LineC和Dir；

其中，Fid表示山谷线编号，初始化为-1，LineC表示判断nm_i相对山谷线位置的有向线段，Dir表示nm_i在有向线段的左右侧，Dir取值为“left”或“right”；

(3-2)读取集合ValleyL中的一线要素vl_i包含的点要素集合VLPt＝{vp_j|j＝1,2,…,b}，vp_j表示线要素vl_i上第j个点要素，b为点要素的个数；

(3-3)循环取集合VLPt中的相邻两点vp_j和vp_j+1，若有向线段与面要素nm_i相交，将

加入有向线段集合LineSeg，得到集合

其中，

表示第k个有向线段，sp为有向线段的首点，ep为有向线段的尾点，c为有向线段的个数；

(3-4)对集合LineSeg的要素个数c进行以下判断：

若c等于0，则执行步骤(3-2)-(3-3)，继续计算有向线段；

若c等于1，修改nm_i的属性，Fid为vl_i的下标i，属性LineC等于有向线段

属性Dir的值按步骤(3-5)计算；

若c大于1，修改nm_i的属性，Fid为vl_i的下标i，属性LineC等于集合LineSeg中有向线段

的首点sp和有向线段

的尾点ep组成的有向线段

属性Dir的值按步骤(3-5)计算；

(3-5)计算属性Dir的值：

读取面要素nm_i包含的点要素集合NUnionPt＝{np_j|j＝1,2,…,d}，np_j表示面要素nm_i上第j个点要素，d为点要素的个数；

循环取集合NUnionPt中的点np_j，采用判断点在直线左右侧的方法，根据公式(1)计算有向线段

的首尾点与点np_j的面积量S(sp,ep,np_j)；

S(P₁,P₂,P₃)＝(x₁-x₃)×(y₂-y₃)-(y₁-y₃)×(x₂-x₃) (1)

若S大于0，将点np_j加入点要素集合LeftPt；反之，将点np_j加入点要素集合RightPt；

设集合LeftPt、RightPt的点要素个数分别是e、f，若e小于3或f小于3，则根据公式(2)计算Dir的值；若e、f均不小于3，则根据公式(3)计算Dir的值；

其中，S_e、S_f分别表示集合LeftPt、RightPt中的点所构成面的面积；

(3-6)循环执行步骤(3-1)-(3-5)，直至集合NoUnionM中的所有面要素完成属性的添加与赋值；

(3-7)基于规则约束条件，对经属性赋值后的面要素集合NoUnionM按以下步骤递归：

1)读取集合UnionM中的一面要素um_i，将其加入面要素集合TmpPoly；

2)遍历集合NoUnionM，若存在面要素nm_i与um_i是邻接关系，按以下步骤计算几个约束参数，否则，跳转步骤1)：

a)读取nm_i的编号属性Fid，若Fid等于-1，则不计算d)、e)；

b)统计nm_i与集合UnionM中的面要素具有邻接关系的总面要素数count，若count等于1，不计算c)，若Fid不等于-1且count不等于1，不计算e)；

c)比较nm_i与集合UnionM中的面要素具有邻接关系的所有面要素，找到存在最大公共边界的面要素um'；

d)读取nm_i的属性值Dir，采用nm_i的属性值LineC，计算um_i的属性值Dir'；

e)计算nm_i与um_i的公共边界长度占nm_i周长的比例scale；

3)判断约束参数若符合以下四个约束条件之一，则执行操作：将nm_i加入集合TmpPoly，并从集合NoUnionM中移除nm_i；

条件1：Fid等于-1且count大于1且um'等于um_i；

条件2：Fid等于-1且count等于1；

条件3：Fid不等于-1且count大于1且um'等于um_i且Dir'等于Dir；

条件4：Fid不等于-1且count等于1且scale不小于0.5且Dir'等于Dir；

4)面要素关联：若集合TmpPoly的面要素个数大于1，则对其进行面要素的合并，形成新的面要素，并替换原来的面要素um_i，置空TmpPoly；

5)循环执行步骤1)-4)，直至集合UnionM中的面要素um_i遍历完毕。

进一步，步骤(4)包括：

(4-1)循环计算出集合NoUnionM中的面要素nm_i与集合UnionM中具有最大公共边界的面要素um_i，记录um_i的下标，得到下标集合Index＝{in_k|k＝1,2,…,w}，w为集合NoUnionM中面要素的个数；

(4-2)按以下步骤进行面要素关联：

Ⅰ、读取集合UnionM中的面要素um_i，将其加入面要素集合TmpFS；

Ⅱ、遍历下标集合Index，若um_i的下标i等于in_k，则将nm_k加入集合TmpFS；

Ⅲ、将集合TmpFS中的面要素进行合并，形成新的面要素，将其加入面要素集合UnionM'，TmpFS置空；

Ⅳ、循环执行步骤Ⅰ-Ⅲ，得到关联后的面要素集合UnionM'＝{um'_i|i＝1,2,…,h}，um'_i表示第i个山体边界面要素，h为面要素的个数，将集合UnionM'写入面要素图层。

有益效果：本发明基于山体边界面数据，以地形特征线作为约束条件，关联符合条件的面，实现了地理尺度山体边界的自动提取，提取后的山体较为完整，尺度相对合理，可直接满足地理研究的需要，因此可成为构建山体对象的基础，并为地貌学和人文学的特征分析、地质灾害后的山体重构以及地形多尺度综合提供了有效支持。

附图说明

图1为本发明方法的流程图

图2为DEM与相对破碎的山体边界面图层示意图；

图3为DEM与山脊线线图层示意图；

图4为DEM与山谷线线图层示意图；

图5为经消除处理后的山体边界面图层示意图；

图6为实施例步骤(2)中未关联面要素和关联后的面要素图层示意图；

图7为实施例步骤(3-7)中满足四种约束条件的情况示意图，其中(a)为条件2，(b)为条件1，(c)为条件3，(d)为条件4)；

图8为实施例步骤(3-7)递归过程中未关联面要素和关联后的面要素图层示意图，其中(a)-(f)分别为6次递归结果；

图9为实施例步骤(4)中关联后的山体边界面要素图层示意图；

图10为地理尺度的山体边界面要素图层与山体阴影叠加示意图。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

实施例：一种山体边界自动提取方法，如图1所示，具体操作如下：

步骤(1)具体包括：

(1-1)加载DEM数据，借助ArcGIS软件，基于地形表面水文分析原理，获取矢量的山体边界面图层(图2)、山脊线线图层(图3)和山谷线线图层数据(图4)；

(1-2)针对相对破碎的山体边界面图层，运用ArcGIS中的消除工具，以面积字段Area阈值，消除小于100m²的图斑，调整阈值的大小，经7次消除，最终消除面积小于3万m²的图斑，得到山体边界面要素图层(图5)。

步骤(2)具体包括：

(2-1)加载山体边界面要素图层、山脊线线要素图层和山谷线线要素图层，得到山体边界面要素集合Mount＝{mt_i|i＝1,2,…,437}、山脊线线要素集合RidgeL＝{rl_i|i＝1,2,...,8}和山谷线线要素集合ValleyL＝{vl_i|i＝1,2,…,7}；

(2-2)遵循山脊线线要素个数决定地理尺度的山体边界面要素个数(即为8)的原则，对集合Mount中的每一个面要素mt_i按以下步骤进行属性添加与赋值：

①给面要素mt₁添加属性Id，初始化为-1；

②读取集合RidgeL中的线要素rl₁包含的点要素集合RLPt＝{rp_j|j＝1,2,…,924}；

③遍历集合RLPt，mt₁不包含线要素rl₁上的点，编号属性Id不变；

④循环执行步骤②-③，直至集合RidgeL中的线要素遍历完毕；

(2-3)根据编号属性Id的不同，对面要素mt_i分组：遍历集合Mount中的每一个面要素mt_i，mt₁的Id等于-1，将mt₁加入集合NoUnionM，mt₃的Id等于1，将mt₃加入集合UnionS₁，……，mt₄₃₅的Id等于8，将mt₄₃₅加入集合UnionS₈……遍历完毕，得到集合NoUnionM＝{nm_i|i＝1,2,…,269}(图6灰色面要素)和以集合为元素的集合SET＝{UnionS₁,UnionS₂,…,UnionS₈}；

(2-4)面要素关联：对SET中的每一个面要素集合(UnionS₁-UnionS₈)，进行面要素合并，形成新的面要素，得到关联后的面要素集合UnionM＝{um_i|i＝1,2,…,8}(图6未填充面要素)。

步骤(3)具体包括：

(3-1)给集合NoUnionM中的一面要素nm₈₃添加属性Fid、LineC和Dir；

(3-2)读取集合ValleyL中的一线要素vl₄包含的点要素集合VLPt＝{vp_j|j＝1,2,…,859}；

(3-3)循环取集合VLPt中的相邻两点，找到vp₇₆₀和vp₇₆₁使有向线段

与面要素nm₈₃相交，将

加入有向线段集合LineSeg中，遍历完毕，得到集合

共35个有向线段；

(3-4)对集合LineSeg的要素个数35进行判断：35大于1，修改nm₈₃的属性Fid为vl₄的下标4，属性LineC等于集合LineSeg中第一个有向线段

的首点vp₇₆₀和最后一个有向线段

的尾点vp₈₅₉组成的有向线段

属性Dir的值按步骤(3-5)计算；

(3-5)计算属性Dir的值：

读取面要素nm₈₃包含的点要素集合NUnionPt＝{np_j|j＝1,2,…,217}；

循环取集合NUnionPt中的点np₁，根据公式(1)计算有向线段

的首点vp₇₆₀、尾点vp₈₅₉与点np₁这三点的面积量S(vp₇₆₀,vp₈₅₉,np₁)＝595212.70；S大于0，将点np₁加入点要素集合LeftPt，遍历完毕，得到集合LeftPt＝{np₁,np₂,…,np₂₁₇}和集合RightPt＝{np₇₈,np₇₉,…,np₁₄₃}，，点要素个数分别是162、54；

集合LeftPt、RightPt的点要素个数均大于3，则根据公式(3)计算S_e＝535959.49、S_f＝21218.17，S_e>S_f，故Dir的值为“left”；

(3-6)循环执行步骤(3-1)-(3-5)，直至集合NoUnionM中的所有面要素完成属性的添加与赋值；

(3-7)基于规则约束条件，对经属性赋值后的面要素集合NoUnionM按以下步骤递归(图8(a)-(f))：

1)读取集合UnionM中的一面要素um₁，将其加入面要素集合TmpPoly；

2)遍历集合NoUnionM，存在面要素nm₂与um₁是邻接关系，按以下步骤计算2个约束参数：

a)读取nm₂的编号属性Fid等于-1；

b)统计nm₂与集合UnionM中的面要素具有邻接关系的总面要素数count等于1；

3)约束参数符合条件2(图7(a))，则将nm₂加入集合TmpPoly，并从集合NoUnionM中移除nm₂；

继续取集合NoUnionM中的面要素，分别列举符合条件1、3和4的情况：

①存在面要素nm₁₀₈与um₁是邻接关系，计算3个约束参数：

a)读取nm₁₀₈的编号属性Fid等于-1；

b)统计nm₁₀₈与集合UnionM中的面要素具有邻接关系的总面要素数count等于2；

c)比较nm₁₀₈与集合UnionM中的面要素具有邻接关系的2个面要素um₁、um₄，找到存在最大公共边界的面要素um₁；

判断约束参数符合条件1的情况(图7(b))；

②存在面要素nm₁₀₅与um₁是邻接关系，计算5个约束参数：

a)读取nm₁₀₅的编号属性Fid等于4；

b)统计nm₁₀₅与集合UnionM中的面要素具有邻接关系的总面要素数count等于2；

c)比较nm₁₀₅与集合UnionM中的面要素具有邻接关系的2个面要素um₁、um₅，找到存在最大公共边界的面要素um₁；

d)读取nm₁₀₅的属性值Dir为“left”，采用nm₁₀₅的属性值LineC，计算um₁的属性值Dir'为“left”；

判断约束参数符合条件3的情况(图7(c))；

③存在面要素nm₁₀₆与um₁是邻接关系，计算5个约束参数：

a)读取nm₁₀₆的编号属性Fid等于1；

b)统计nm₁₀₆与集合UnionM中的面要素具有邻接关系的总面要素count等于1；

c)读取nm₁₀₆的属性值Dir为“right”，采用nm₁₀₆的属性值LineC，计算nm₁₀₆的属性值Dir'为“right”；

d)计算nm₁₀₆与um₁的公共边界长度占nm₁₀₆周长的比例scale等于0.59；

判断约束参数符合条件4的情况(图7(d))；

4)面要素关联：集合TmpPoly的面要素个数为36，对其进行面要素的合并，形成新的面要素，并替换原来的面要素um₁，置空TmpPoly；

5)循环执行步骤1)-4)，直至集合UnionM中的面要素um_i遍历完毕。

步骤(4)具体包括：

(4-1)循环计算出集合NoUnionM中的面要素nm₁与集合UnionM中具有最大公共边界的面要素um₅，记录um₅的下标5，遍历完毕，得到下标集合Index＝{5,1,…,7}，共7个下标；

(4-2)按以下步骤进行面要素关联：

Ⅰ、读取集合UnionM中的面要素um₁，将其加入面要素集合TmpFS；

Ⅱ、遍历下标集合Index，um₁的下标1等于in₂，则将nm₂加入集合TmpFS；

Ⅲ、将集合TmpFS中的面要素进行合并，形成新的面要素，将其加入面要素集合UnionM'，TmpFS置空；

Ⅳ、循环执行步骤Ⅰ-Ⅲ，得到关联后的面要素集合UnionM'＝{um'_i|i＝1,2,…,8}，将集合UnionM'写入面要素图层(图9)。

将以上实施例最终得到的地理尺度的山体边界面要素图层，与相应的山体阴影图进行叠加(图10)发现：山体边界基本符合地形情况，效果良好。

Claims (5)

1.一种山体边界自动提取方法，其特征在于：包括以下步骤：

(1)获取山体边界面、山脊线和山谷线数据；

(2)基于山脊线数据，关联山脊线经过的山体边界面，初次形成山体边界；

(3)查找与山体边界具有邻接关系的山体边界面，并以山谷线为约束，关联符合条件的面；其中关联符合约束条件为根据编号属性Id的不同，对面要素mt_i分组：遍历集合Mount中的每一个面要素mt_i，Mount为山体边界面要素集合；当mt_i的Id等于-1时，将mt_i加入面要素集合NoUnionM；当Id等于1时，将mt_i加入面要素集合UnionS₁，当Id等于2时，将mt_i加入面要素集合UnionS₂，……，当Id等于n时，将mt_i加入面要素集合UnionS_n；遍历完毕，得到集合NoUnionM＝{nm_i|i＝1,2,…,q}和以集合为元素的集合SET＝{UnionS₁,UnionS₂,…,UnionS_n}；其中，nm_i表示不符合关联条件的面要素，q为面要素的个数，UnionS_n表示具有相同Id的面要素集合，n为面要素集合的个数；

(4)将规则约束条件外的山体边界面进行关联，最终形成地理尺度的山体边界。

2.根据权利要求1所述的山体边界自动提取方法，其特征在于：步骤(1)通过加载DEM数据，运用ArcGIS并基于地形表面水文分析原理获取矢量的相对破碎的山体边界面图层、山脊线线图层和山谷线线图层数据，针对相对破碎的山体边界面图层，运用ArcGIS中的消除工具，参照山体边界面要素的面积字段Area的值设置阈值Value，消除面积小于阈值Value的山体边界面要素，不断调整阈值的大小，经多次消除，得到山体边界面要素图层。

3.根据权利要求2所述的山体边界自动提取方法，其特征在于：步骤(2)包括：

(2-1)加载山体边界面要素图层、山脊线线要素图层和山谷线线要素图层，分别得到山体边界面要素集合Mount＝{mt_i|i＝1,2,…,m}、山脊线线要素集合RidgeL＝{rl_i|i＝1,2,…,n}和山谷线线要素集合ValleyL＝{vl_i|i＝1,2,…,r}；

其中，mt_i表示第i个面要素，m为面要素的个数；rl_i表示第i个山脊线线要素，n为山脊线线要素的个数；vl_i表示第i个山谷线线要素，r为山谷线线要素的个数；

(2-2)遵循山脊线线要素个数n决定地理尺度的山体边界面要素个数的原则，对集合Mount中的每一个面要素mt_i按以下步骤进行属性添加与赋值：

①给面要素mt_i添加属性Id，Id表示山脊线编号，初始化为-1；

②读取集合RidgeL中的线要素rl_i包含的点要素集合RLPt＝{rp_j|j＝1,2,…,a}，rp_j表示线要素rl_i上第j个点要素，a为点要素的个数；

③判断集合RLPt中的点rp_j是否在面要素mt_i内，当存在任一点符合条件时，修改mt_i的编号属性Id为rl_i的下标i；

④循环执行步骤②-③，直至集合RidgeL中的线要素遍历完毕；

(2-3)面要素关联：对SET中的每一个面要素集合UnionS_n，进行面要素合并，形成新的面要素，得到关联后的面要素集合UnionM＝{um_i|i＝1,2,…,g}，um_i表示关联后的面要素，g为面要素的个数。

4.根据权利要求3所述的山体边界自动提取方法，其特征在于：步骤(3)包括：

(3-1)给集合NoUnionM中的一面要素nm_i添加属性Fid、LineC和Dir；

其中，Fid表示山谷线编号，初始化为-1，LineC表示判断nm_i相对山谷线位置的有向线段，Dir表示nm_i在有向线段的左右侧，Dir取值为“left”或“right”；

(3-2)读取集合ValleyL中的一线要素vl_i包含的点要素集合VLPt＝{vp_j|j＝1,2,…,b}，vp_j表示线要素vl_i上第j个点要素，b为点要素的个数；

(3-3)循环取集合VLPt中的相邻两点vp_j和vp_j+1，若有向线

与面要素nm_i相交，将

加入有向线段集合LineSeg，得到集合

其中，

表示第k个有向线段，sp为有向线段的首点，ep为有向线段的尾点，c为有向线段的个数；

(3-4)对集合LineSeg的要素个数c进行以下判断：

若c等于0，则执行步骤(3-2)-(3-3)，继续计算有向线段；

若c等于1，修改nm_i的属性，Fid为vl_i的下标i，属性LineC等于有向线段

属性Dir的值按步骤(3-5)计算；

若c大于1，修改nm_i的属性，Fid为vl_i的下标i，属性LineC等于集合LineSeg中有向线段

的首点sp和有向线段

的尾点ep组成的有向线段

属性Dir的值按步骤(3-5)计算；

(3-5)计算属性Dir的值：

读取面要素nm_i包含的点要素集合NUnionPt＝{np_j|j＝1,2,…,d}，np_j表示面要素nm_i上第j个点要素，d为点要素的个数；

循环取集合NUnionPt中的点np_j，采用判断点在直线左右侧的方法，根据公式(1)计算有向线段

的首尾点与点np_j的面积量S(sp,ep,np_j)；

S(P₁,P₂,P₃)＝(x1-x3)×(y2-y₃)-(y₁-y3)×(x₂-x₃) (1)

其中，(x₁,y₁)、(x₂,y₂)、(x₃,y₃)分别表示点P₁、P₂、P₃的坐标值；

若S大于0，将点np_j加入点要素集合LeftPt；反之，将点np_j加入点要素集合RightPt；

设集合LeftPt、RightPt的点要素个数分别是e、f，若e小于3或f小于3，则根据公式(2)计算Dir的值；若e、f均不小于3，则根据公式(3)计算Dir的值；

其中，S_e、S_f分别表示集合LeftPt、RightPt中的点所构成面的面积；

(3-6)循环执行步骤(3-1)-(3-5)，直至集合NoUnionM中的所有面要素完成属性的添加与赋值；

(3-7)基于规则约束条件，对经属性赋值后的面要素集合NoUnionM按以下步骤递归：

1)读取集合UnionM中的一面要素um_i，将其加入面要素集合TmpPoly；2)遍历集合NoUnionM，若存在面要素nm_i与um_i是邻接关系，按以下步骤计算几个约束参数，否则，跳转步骤1)：

a)读取nm_i的编号属性Fid，若Fid等于-1，则不计算d)、e)；

b)统计nm_i与集合UnionM中的面要素具有邻接关系的总面要素数count，若count等于1，不进行c)，若Fid不等于-1且count不等于1，不计算e)；

c)比较nm_i与集合UnionM中的面要素具有邻接关系的所有面要素，找到存在最大公共边界的面要素um'；

d)读取nm_i的属性值Dir，采用nm_i的属性值LineC，计算um_i的属性值Dir'；

e)计算nm_i与um_i的公共边界长度占nm_i周长的比例scale；

3)判断约束参数若符合以下四个约束条件之一，则执行操作：将nm_i加入集合TmpPoly，并从集合NoUnionM中移除nm_i；

条件1：Fid等于-1且count大于1且um'等于um_i；

条件2：Fid等于-1且count等于1；

条件3：Fid不等于-1且count大于1且um'等于um_i且Dir'等于Dir；

条件4：Fid不等于-1且count等于1且scale不小于0.5且Dir'等于Dir；

4)面要素关联：若集合TmpPoly的面要素个数大于1，则对其进行面要素的合并，形成新的面要素，并替换原来的面要素um_i，置空TmpPoly；

5)循环执行步骤1)-4)，直至集合UnionM中的面要素um_i遍历完毕。

5.根据权利要求4所述的山体边界自动提取方法，其特征在于：步骤(4)包括：

(4-1)循环集合NoUnionM中的每一面要素nm_i，从集合UnionM中找到与nm_i具有最大公共边界的面要素um_i，记录um_i的下标，得到下标集合Index＝{in_k|k＝1,2,…,w}，w为集合NoUnionM中面要素的个数；

(4-2)按以下步骤进行面要素关联：

Ⅰ、读取集合UnionM中的面要素um_i，将其加入面要素集合TmpFS；

Ⅱ、遍历下标集合Index，若um_i的下标i等于in_k，则将nm_k加入集合TmpFS；

Ⅲ、将集合TmpFS中的面要素进行合并，形成新的面要素，将其加入面要素集合UnionM'，TmpFS置空；

Ⅳ、循环执行步骤Ⅰ-Ⅲ，得到关联后的面要素集合UnionM'＝{unm_i|i＝1,2,…,h}，unm_i表示第i个山体边界面要素，h为面要素的个数，将集合UnionM'写入面要素图层。

Images