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

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

背景技术：

目前，在已存的地形模型中，通常采用规则采样点、三角面片、等高线等形式来表达山体，而避免直接表达山体对象。实际上，地形中存在山峰与河谷等现象，也能明确其位置和形态，可见地形中存在对象。在地貌学中，山体作为一种重要的地形对象，其构建主要有以下必要性：(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)发现：山体边界基本符合地形情况，效果良好。