CN106023161A - 一种破碎山脊线的自动关联方法 - Google Patents

Abstract

本发明公开了一种破碎山脊线的自动关联方法。该方法包括如下步骤：(1)对破碎山脊线进行预处理，并设定关联约束条件；(2)对破碎山脊线进行空间剖分处理，得到每条山脊线的空间影响范围；(3)基于空间剖分结果，构建山脊线的ARG模型；(4)基于关联约束条件，提取符合关联约束条件的山脊线要素对，进行山脊线要素的关联；(5)对关联处理后的山脊线重构ARG模型，继续提取符合关联约束条件的山脊线要素对，进行关联处理，循环这一过程，直至山脊线图层中不存在符合关联约束条件的邻接山脊线要素对，则关联处理结束。

Description

一种破碎山脊线的自动关联方法

技术领域

本发明属于地理信息技术应用领域，具体涉及一种破碎山脊线的自动关联方法。

背景技术

山脊线是地形起伏变化的分界线，在地形表示方面起着重要作用。数字地面模型(DEM)数据中隐含了相关的山脊线信息，因而，基于DEM提取山脊线已得到了广泛的研究。

传统的基于DEM提取山脊线的研究中主要采用基于地形表面流水分析的方法。以山脊线的物理特性为依据，通过模拟和分析地形表面流水的运动状况得到分水线(山脊线)(参见《提取山脊线和山谷线的一种新方法》，黄培之，武汉大学学报：信息科学版，2001年，26(3))。由于流水模拟的方法在地形破碎的地区噪音非常大，难免会遗漏部分地形特征点，导致局部地形特征线在连接时较为困难，进而提取的山脊线易破碎并产生很多毛刺(参见《利用等高线数据提取山脊(谷)线算法研究》，靳海亮，康建荣，高井祥，武汉大学学报：信息科学版，2005年，30(9))。

为了克服这一不足，传统的提取方法采用对地形数据进行洼地填充的方式，来剔除噪音，进而减少破碎和毛刺。该方法在一定程度上提高了山脊线的连续性和正确性，但提取的山脊线仍需借助人工判别的手段对山脊线进行关联连接处理后，方可运用到具体分析和应用中，处理效率低，且易引入人为错误。

发明内容

本发明针对传统的基于DEM提取的山脊线中出现的毛刺和山脊线破碎问题，提出一种破碎山脊线的自动关联方法。该方法可实现毛刺的剔除和破碎山脊线的关联，满足自动化提取山脊线的需要。本发明方法同样适用于破碎山谷线的自动关联处理。

本发明的技术方案如下：

首先，对破碎山脊线进行预处理，并设定关联约束条件；然后，对破碎山脊线进行空间剖分处理，构建山脊线的邻接关系图(Attributed Relation Graph，简称ARG)模型；其次，基于关联约束条件，提取符合关联约束条件的山脊线要素对，并进行山脊线要素的关联；最后，对关联处理后的山脊线重构ARG模型，继续提取符合关联约束条件的山脊线要素对，进行关联处理，直至山脊线图层中不存在符合关联约束条件的邻接山脊线要素对，则关联处理结束。

上述方法的具体步骤如下：

(一)破碎山脊线的预处理

步骤11：加载山脊线图层，并将相应的山脊线要素读入集合R₁＝{r_i|i＝1，2，…，n₁}，n₁为山脊线要素的条数；

步骤12：依据山脊线图层信息，设定关联约束条件：设置相对走向约束值SThreshold、相对趋势约束值IThreshold和最短距离约束值DThreshold；

步骤13：剔除非线状要素的山脊线要素，并读入集合R₂＝{r_i|i＝1，2，…，n₂}，n₂为剔除了非线状要素的山脊线的个数；

步骤14：针对存在公共点的两山脊线要素进行去毛刺处理，并读入集合R₃＝{r_i|i＝1，2，…，n₃}，n₃为预处理后的山脊线要素的个数；具体处理过程如下：

a)将山脊线要素r_a和r_b分别读入点集PointA＝{p₁，p₂，…，p_u}和PointB＝{p₁，p₂，…，p_v}，u、v分别为点集PointA、PointB中点的个数；

b)计算山脊线要素r_a和r_b的公共交点，存入点集CrossPoint＝{p₁，p₂，…，p_c}，c为公共交点的个数；

c)若公共交点的个数c＝1，则按均值法或最大频数法计算山脊线要素r_a的走向Strike_a和山脊线要素r_b的走向Strike_b；否则，跳转下面的步骤f)；

d)若|Strike_a-Strike_b|＜SThreshold，则对山脊线要素r_a和r_b进行合并处理；否则，跳转e)；

e)对山脊线要素r_a和r_b进行打断操作，山脊线要素r_a和山脊线要素r_b均删除公共交点CrossPoint；

f)公共交点的个数c≥2，计算山脊线要素r_a的长度length_a和山脊线要素r_b的长度length_b；

若length_a＞length_b，则删除长度较短山脊线要素r_b；否则，删除山脊线要素r_a；

(二)山脊线要素的空间剖分处理

步骤21：将山脊线要素集合R₃读入点集P＝{p_i(x_i，y_i)|i＝1，2，…，m}，m为点要素的个数，并使用点集P的“山脊线要素”属性来记录该点所归属的山脊线要素；

步骤22：基于点集P，构建Voronoi图Vor；以图层的最小外包矩形为界，构建点集P的Voronoi图Vor＝{T_i|i＝1，２，…，m′}，T_i表示第i个Thicssen多边形，m'为构成Voronoi图Vor的Thicssen多边形T的个数；并将点集P中各点的“山脊线要素”属性，赋给包含该点的相应Thicssen多边形T_i的“山脊线要素”属性；

步骤23：基于Thicssen多边形T_i的“山脊线要素”属性合并Thicssen多边形T_i，得到每条山脊线要素r_i的影响范围；依次将Voronoi图Vor中具有相同“山脊线要素”属性值的Thicssen多边形T_i合并为一个多边形T'_i，生成图Vor'＝{T'_i|i＝1,2,…,n'}，n'为图Vor'的Thicssen多边形的个数；图Vor'即山脊线要素集合R₃基于最小外包矩形的空间影响范围集合，T'_i代表第i条山脊线要素r_i的影响范围；

(三)基于空间剖分结果，构建山脊线要素的ARG模型G＝(V，E)

步骤31：计算山脊线要素的ARG模型G＝(V，E)中节点集合V＝{v₁，v₂，…，v_n3}及其相关属性；其中，V中每个节点代表一个山脊线要素对象，E为ARG模型的边集合；具体步骤如下：

a)计算山脊线要素r_i的剖分范围T'_i的质心点c_i得到节点v_i，并存入节点集合V中；i为循环变量，i∈[1，n₃]，初值为1；

b)读取山脊线要素r_i的点集C_i＝{c_i1，c_i2，…，c_it}，t为山脊线要素r_i中点的个数；

c)基于点集C_i分别计算直线c_i1c_i2、c_i2c_i3、…、c_it-1c_it的倾斜角α₁₂、α₂₃、…、α_t-1t，并按最大频数法或均值法计算山脊线要素的走向Strike_i，Strike_i∈[0°，180°)，并存入节点v_i的属性中；

d)基于点集C_i，计算山脊线要素r_i的最小凸包的点集H_i＝{h_i1，h_i2，…，h_is}，s为点的个数，并存入节点v_i的属性中；

e)求出山脊线要素r_i的剖分范围T'_i的顶点集Vt_i＝{p_i1，p_i2，…，p_il}，并存入节点v_i的属性，l为顶点集Vt_i中点的个数；

f)若循环变量i＜n₃，则i加1，继续执行步骤a)-e)，直至所有节点及其属性均提取完成；

步骤32：基于节点集合V，计算山脊线的ARG模型G＝(V，E)的边集合w₁为边的个数，具体步骤如下：

a)读取节点集合V中的节点v_i和v_j所对应的影响范围T'_i的顶点集Vt_i和影响范围T'_j的顶点集Vt_j；其中，i为循环变量，i∈[1，n₃-1]，初值为1；j为循环变量，j∈[i+1，n₃]，初值为i+1；

b)判断顶点集Vt_i和Vt_j是否存在连续的两公共点，若存在，则山脊线要素r_i的影响范围T'_i和r_j的影响范围T'_j为邻接拓扑关系，连接节点v_i和v_j，得到ARG模型一条边e＝(v_i，v_j)，称节点v_i与v_j互为邻接节点，山脊线要素r_i和山脊线要素r_j为边e所对应的邻接山脊线要素对，并将边e存入边集合E；否则，跳转步骤c)；

c)若j≤n₃，则循环变量j加1，继续执行步骤a)和b)；否则，循环变量i加1，执行下一步骤d)；

d)若循环变量i≤n₃-1，则循环变量j＝i+1，继续执行步骤a)-c)；否则，ARG模型的边集合已提取结束；

步骤33：计算ARG模型G＝(V，E)的边集合的属性{S，I，D}，S表示该边所对应的邻接山脊线要素对的相对走向关系，I表示邻接山脊线要素对的趋势关系，D表示邻接山脊线要素对的最短距离关系；具体步骤如下：

a)读取边集合E的第k条边e_k，k为循环变量，k∈[1，w₁]，获取节点v_a和节点v_b的走向属性Strike_a和Strike_b；

b)计算由节点v_a与v_b所对应的山脊线要素r_a与r_b相对走向关系S_k＝Strike_a-Strike_b，并存入边e_k的属性中；

c)计算节点v_a和v_b构成的直线v_av_b的倾斜角α_k，α_k∈[0°，180°)，作为边e_k的相对趋势关系I_k，并存入边e_k的属性中；

d)获取边e_k的节点v_a的最小凸包的点集H_a＝{h_a1，h_a2，…，h_ap}和节点v_b的最小凸包的点集H_b＝{h_b1，h_b2，…，h_bq}，p和q均为点的个数；

e)按式(1)计算h_ai与h_bj的距离d_ij，并存入最短距离集合Distance中，其中i和j为循

环变量，初值均为1，i∈[1，p]，j∈[1，q]；

$d_{ij}=\sqrt{{(x_{ai}-x_{bj})}^2+{(y_{ai}-y_{bj})}^2}---\left(1\right)$

其中，(x_ai，y_ai)和(x_bj，y_bj)分别为h_ai和h_bj的坐标；

f)若循环变量j≤q，循环变量j加1，继续执行上一步骤e)；否则，循环变量i加1，执行下一步骤g)；

g)若循环变量i≤p，则循环变量j＝1，继续执行步骤e)；否则，距离d_ij计算结束，执行下一步骤h)；

h)读取最短距离集合Distance，取集合Distance中的最小值D_k＝min{d₁₁，d₁₂，…，d_pq}作为边e_k的最短距离关系，并存入边e_k的属性中；

i)读取边集合E的下一条边e_k+1，循环以上步骤，直至边集合E中每条边的相对走向关系S、相对趋势关系I和最短距离关系D均计算结束；

(四)山脊线的关联处理

步骤41：遍历ARG模型的边集合E＝{e₁，e₂，…，e_w1}中的每条边e及属性{S，I，D}，并获取边e的节点v_a和节点v_b的Strike_a和Strike_b；

步骤42：对每条边e均进行关联约束条件的判断，若均满足以下条件：

$\left\{\begin{array}{c}\left|S\right|<SThreahold\\ \left|{\mathrm{Strike}}_a-I\right|<Ithreshold\\ \left|{\mathrm{Strike}}_b-I\right|<Ithreshold\\ D<DThreshold\end{array}\right.$

则边e符合关联约束条件，节点v_a和v_b所对应的山脊线要素r_a和r_b满足关联约束条件，并将边e记为e'存入边集合E'；若不满足，则跳转下一步骤43；

步骤43：循环执行步骤42，直至所有边e均进行关联约束条件的判断后，得到边集合w₂为边集合E′内的边数；

步骤44：基于边集合进行山脊线要素的关联处理，具体步骤如下：

a)读取边集合E′的第i条边e'_i，并提取边e'_i的节点v_a和v_b，其中i为循环变量，初值为w₂；

b)读取并提取边集合E'的第j条边e'_j的节点v_c和v_d，其中j为循环变量，初值为w₂-1；

c)若节点v_a、v_b、v_c和v_d，均满足v_a≠v_c，v_a≠v_d，v_b≠v_c，v_b≠v_d，则边e'_i与边e'_j间不存在共同关联，循环变量j减1，执行步骤e)；否则执行步骤d)；

d)由于边e'_i和边e'_j存在共同关联，将e'_j存入集合L_i(L_i存储与e'_i存在共同合并的边e')循环变量j减1，并跳转步骤e)；

e)若j≥1，则跳转步骤b)，否则执行下一步骤f)；

f)将边e'_i存入集合L_i，依次读取L_i中边e'，并按式(2)分别计算每条边e'所对应的节点v_e和v_f对应的山脊线要素r_e和r_f的差异度Difference，并存入集合Sum中；

Difference＝p₁×|S|+p₂×(|Strike_e-I|+|Strike_e-I|)+p₃×D(2)

其中，p₁，p₂，p₃分别代表相对走向关系S、相对趋势关系I、最短距离关系D的权值，p₁、p₂、p₃均∈(0，1)之间，且满足p₁+p₂+p₃＝1；

g)遍历集合Sum，并获取值最小的差异度Difference所对应的边e'，对边e'的节点v_e和v_f所对应的山脊线要素r_e和r_f进行关联处理；

h)循环变量i减1，若i≥2，则循环变量j＝i-1，继续执行步骤a)-步骤g)；否则，输出关联后的山脊线要素集合LM₁＝{r_i|i＝1，2，…，l₁}，l₁为第一次关联处理后的山脊线的条数；

步骤45：循环执行步骤(一)至(四)，直至循环K次后该山脊线要素集合LM_K不存在符合关联约束条件的邻接山脊线要素对，则连接处理结束，输出山脊线要素集合LM_K＝{r_i|i＝1，2，…，l_K}。

本发明提出一种破碎山脊线的自动关联方法。该方法以关联约束为前提条件，提取符合关联约束条件的邻接山脊线要素对进行关联处理。若某一山脊线要素与多个山脊线要素间均满足关联约束条件，即存在共同关联对象时，则计算该山脊线要素与其所有邻接山脊线要素的差异度，选取与其差异度最小的邻接山脊线要素作为关联对象，进行关联连接。本方法同样适用于破碎山谷线的自动关联。与传统的方法相比，本方法采用关联处理的方式，在一定程度上克服了由于DEM数据抗躁性较差，而引起的地形特征点连接困难的问题，并对毛刺进行了剔除处理，进而满足了自动化提取山脊线的需要。

附图说明

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

图2庐山山脊线数据；

图3预处理后的庐山山脊线数据；

图4庐山山脊线图层的空间剖分结果图V'；

图5庐山山脊线的ARG模型；

图6 ARG模型边E的属性表；

图7为实施例中第1次关联处理后的庐山山脊线图层的空间剖分结果图；

图8为实施例中第6次关联处理后的庐山山脊线图层的空间剖分结果图；

图9为实施例中第6次关联处理后的符号化效果图。

注：图3-图5均为第一次关联处理过程图。

具体实施方式

本发明的方法是适用于破碎山脊线和山谷线的自动关联处理，下面结合附图，给出破碎山脊线的实施例，来进一步说明本发明的效果。具体流程如图1所示。

本实例的实验区域位于江西北部，东经115°57'-116°11'，北纬29°33'-29°35'。实验数据为基于地形表面流水分析原理，借助ArcGIS软件，提取得到实验区域的山脊线矢量数据。

以庐山破碎山脊线的关联为例，具体的实施过程如下：

(一)破碎山脊线的预处理

步骤11：加载山脊线图层，如图2所示。并将相应的山脊线要素读入集合R₁中。本实例中R₁＝{r_i|i＝1，２，…，974}，共974个山脊线要素；

步骤12：依据山脊线图层信息设定关联约束条件：相对走向约束值SThreshold、相对趋势约束值IThreshold和最短距离约束值DThreshold。本实例中，相对走向约束值SThreshold设为25°、相对趋势约束值IThreshold设为45°和最短距离约束值DThreshold设为150m；

步骤13：剔除非线状要素的山脊线要素，本实例中共剔除113条山脊线要素。并读入集合R₂，共861条山脊线要素；

步骤14：针对存在公共点的两山脊线要素进行去毛刺处理，并读入集合R₃＝{r_i|i＝1，2，…，n₃}。本实例中，以有公共点的山脊线要素r₇₉₅和山脊线要素r₇₉₆为例，具体步骤如下：

a)将山脊线要素r₇₉₅和r₇₉₆分别读入点集PointA＝{p₁，p₂，…，p₇}和PointB＝{p₁，p₂，…，p₁₄}。PointA中点的个数u＝7，PointB中点的个数v＝14；

b)计算山脊线要素r₇₉₅和r₇₉₆的公共交点，存入集合CrossPoint，公共交点个数c＝1；

c)公共交点的个数c＝1，按均值法计算山脊线要素r₇₉₅的走向Strike₇₉₅＝41.03°，按最大频数法计算山脊线要素r₇₉₆的走向Strike₇₉₆＝28.08°；

d)|Strike₇₉₅-Strike₇₉₆|＝12.95°＜SThreshold＝25°，故对山脊线要素r₇₉₅和r₇₉₆进行合并处理；

对有公共点的山脊线要素均进行步骤14操作后，得到R₃＝{r_i|i＝1，２，…，598}，共n₃＝598条预处理后的山脊线要素，如图3所示。

(二)山脊线要素的空间剖分处理

步骤21：将山脊线要素集合R₃读入点集P＝{p_i(x_i，y_i)|i＝1，2，…，4352}，共4352个点要素，并使用点集P的“山脊线要素”属性来记录该点所归属的山脊线要素。如点集P中的点要素1，点要素2，点要素3，点要素4，点要素5均归属于山脊线r₁；

步骤22：基于点集P，构建Voronoi图Vor。以图层的最小外包矩形为界，本实例中山脊线要素图层的最小外包矩形的顶点坐标为{(399064.06，3274424.25)，(404729.06，3274424.25)，(404729.06，3270731.75)，(399064.06，3270731.75)，(399064.06，3274424.25)}。以最小外包矩形为界为界，得到Voronoi图Vor，Vor＝{T_i|i＝1，２，…，4352}共有4352个Thicssen多边形构成。并将点集P中各点的“山脊线要素”属性，赋给包含该点的相应Thicssen多边形的“山脊线要素”属性，如如Thicssen多边形T₁、T₂、T₃、T₄、T₅分别包含点要素1、点要素2、点要素3、点要素4、点要素5，则T₁、T₂、T₃、T₄、T₅“山脊线要素”的属性值均为1；

步骤23：基于Thicssen多边形的“山脊线要素”属性合并Thicssen多边形，得到每条山脊线要素r_i的剖分范围。以“山脊线要素”属性值为1为例，则“山脊线要素”属性值均为1的Thicssen多边形T₁、T₂、T₃、T₄、T₅，合并为一个多边形T'₁，多边形T'₁代表山脊线要素r₁的影响范围。依次合并后，本实例最终得到合并后的图Vor′＝{T′_i|i＝1，２，…，598}，如图4所示。图Vor'共由598个多边形组成，每个多边形对应一个山脊线要素的影响范围。

(三)基于空间剖分结果，构建山脊线要素的ARG模型G＝(V，E)

步骤31：计算山脊线要素的ARG模型G＝(V，E)中节点集合V＝{v₁，v₂，…，v₅₉₈}及其相关属性。其中，V中每个节点代表一个山脊线要素对象，E为ARG模型的边集合。本实例中，以山脊线r₁为例，具体步骤如下：

a)计算山脊线要素r₁的剖分范围T'₁的质心点c₁得到节点v₁＝{402037.969，3274311.5}，并依次存入节点集合V中；

b)读取山脊线要素r₁的点集C₁＝{c₁₁，c₁₂，…，c₁₅}，点的个数t＝5；

c)基于点集c₁＝{c₁₁，c₁₂，…，c₁₅}，按均值法计算直线c₁₁c₁₂、c₁₂c₁₃、c₁₃c₁₄、c₁₄c₁₅的倾斜角α₁₂＝123.77°、α₂₃＝145.35°、α₃₄＝161.92°、α₄₅＝172.48°，并存入节点v₁的属性中；

d)基于点集C₁，计算山脊线要素r₁的最小凸包的点集H₁＝{h₁₁，h₁₂，h₁₃，h₁₄}，并存入节点v₁的属性中；

e)求出山脊线要素r₁的剖分范围T'₁的顶点集Vt₁＝{p₁₁，p₁₂，…，p₁₂₀}，共20个点，并存入节点v₁的属性；

循环以上步骤，直至598个节点及其属性均提取完成；

步骤32：计算ARG模型G＝(V，E)的边集合本实例中，以节点v₁＝{402037.969，3274311.5}和v₂＝{401976.2，3274374.25}为例，具体步骤如下：

a)读取节点集合V中的节点v₁和v₂所对应的空间影响范围T'₁的顶点集Vt₁＝{p₁₁，p₁₂，…，p₁₂₀}和影响范围T'₂的顶点集Vt₂＝{p₂₁，p₂₂，…，p₂₂₁}。此时，循环变量i＝1，j＝2；

b)顶点集Vt₁和Vt₂中存在连续的两公共点，则山脊线要素r₁的影响范围T'₁和r₂的影响范围T'₂为邻接拓扑关系，连接节点v₁和v₂，得到ARG模型的第一条边e＝(v₁,v₂)。称节点v₁与v₂互为邻接节点，山脊线要素r₁和山脊线要素r₂为边e所对应的邻接山脊线要素对，并将边e存入边集合E；

c)循环变量j＝2≤n₃＝598，则循环变量j加1，继续执行步骤a)和b)，直至ARG模型的边集合E均提取结束。本实例中，边集合E＝{e₁，e₂，…，e₁₆₉₅}，共w＝1695条边，如图5所示。

步骤33：计算ARG模型G＝(V，E)的边集合E＝{e₁，e₂，…，e₁₆₉₅}的属性{S，I，D}，S表示该边所对应的邻接山脊线要素对的相对走向关系，I表示邻接山脊线要素对的趋势关系，D表示邻接山脊线要素对的最短距离关系。以边e₁为例，具体步骤如下：

a)读取边集合E的第1条边e₁，获取节点v₁和节点v₂的走向属性和Strike₁＝150.88°和Strike₂＝16.58°；

b)计算节点v₁与v₂所对应的山脊线要素r₁与r₂相对走向关系S＝Strike_a-Strike_b＝134.30°，并存入边e₁的属性中；

c)计算由节点v₁和v₂构成的直线v₁v₂的倾斜角α＝134.55°，作为边e₁的山脊线相对趋势关系I，并存入边e₁的属性中；

d)获取边e₁的两节点v₁和v₂的最小凸包的点集H₁＝{h₁₁，h₁₂，h₁₃，h₁₄}和H₂＝{h₂₁，h₂₂，…，h₂₉}，p＝4，q＝9；

e)按式(1)依次计算d₁₁，存入最短距离Distance集合，循环变量i＝1，j＝1：

$d_{11}=\sqrt{{(402009.063-102014.063)}^2+{(3274339.25-3274351.75)}^2}=13.46;$

f)循环变量j＝1≤q＝9，循环变量j加1，继续执行上一步骤e)，本实例中，最短距离集合Distance＝{13.46，23.49，…，115.30}；

g)读取最短距离集合Distance，取集合Distance中的最小值D＝min{13.46，23.49，…，115.30}＝13.46作为边e₁的最短距离关系，D代表节点v₁与v₂所对应的山脊线要素r₁与r₂的最短距离；

h)读取边集合E的下一条边e₂，循环以上步骤，直至边集合E的每条边的相对走向关系S、相对趋势关系I和最短距离关系D均计算结束。本实例的边集合E＝{e₁，e₂，…，e₁₆₉₅}的属性{S，I，D}计算结果如图6所示。

(四)山脊线的关联处理

步骤41：遍历ARG模型的边集合E＝{e₁，e₂，…，e₁₆₉₅}，本实例以边e₁、e₁₅₁₇、e₁₆₆₆为例，进行关联处理。获取边e₁、e₁₅₁₇和e₁₆₆₆及其属性{S₁，I₁，D₁}、{S₁₅₁₇，I₁₅₁₇，D₁₅₁₇}和{S₁₆₆₆，I₁₆₆₆，D₁₆₆₆}。并获取边e₁的节点v₁的Strike₁和节点v₂的Strike₂、边e₁₅₁₇的节点v₅₀₀的Strike₅₀₀和节点v₅₆₀的Strike₅₆₀和边e₁₆₆₆的节点v₅₅₉的Strike₅₅₉和节点v₅₆₀的Strike₅₆₀；

步骤42：对边e₁、e₁₅₁₇、e₁₆₆₆进行关联约束的判断，边e₁计算结果如下：

边e₁₅₁₇的计算结果如下：

边e₁₆₆₆的计算结果如下：

相对走向关系S₁＞相对走向阈值SThreahold，则边e₁不符合关联约束条件，则节点v₁和v₂所对应的山脊线要素对r₁和r₂不能进行关联处理，跳转下一步骤43；边e₁₅₁₇和边e₁₆₆₆满足关联约束条件，可进行关联处理，将边e₁₅₁₇和边e₁₆₆₆均存入边集合E'；

步骤43：读取边集合E的每条边，循环执行步骤42，直至所有边均进行关联约束的判断后，本实例得到边集合E′＝{e′₁，e′₂，…，e′₁₉₆}，共196条边；

步骤44：基于边集合E′＝{e′₁，e′₂，…，e′₁₉₆}，进行山脊线要素的关联处理。本实例存在共同关联的边e'₁₈₈和边e'₁₆₆为例，具体步骤如下：

a)读取边集合E′＝{e′₁，e′₂，…，e′₁₉₆}中的第188条边e'₁₈₈，并提取边e'₁₈₈的节点v₅₅₉和v₅₆₀。此时，循环变量i＝188；

b)读取并提取边集合E'的第166条边e'₁₆₆的节点v₅₀₀和v₅₆₀。此时，循环变量j＝166。

c)由于节点v₅₀₀、v₅₆₀与v₅₅₉、v₅₆₀间存在共同节点v₅₆₀，均不满足以下条件：v₅₀₀≠v₅₅₉，v₅₀₀≠v₅₆₀，v₅₆₀≠v₅₅₉，v₅₆₀＝v₅₆₀，则边e'₁₈₈与边e'₁₆₆间存在共同关联，执行步骤d)；

d)将边e'₁₆₆存入集合L₁₈₈，循环变量j减1，跳转步骤e)。

e)j≥1，跳转步骤b)，搜寻与边e'₁₈₈存在共同合并的所有边，本实例中，仅边e'₁₆₆与边e'₁₈₈存在共同合并；

f)将边e'₁₈₈存入集合L₁₈₈，依次读取L₁₈₈中边e'₁₆₆与边e'₁₈₈，并按式(2)依次计算边e'₁₆₆所对应的节点v₅₀₀和v₅₆₀对应的山脊线要素r₅₀₀和r₅₆₀间的相似度Difference₁₆₆和边e'₁₈₈所对应的节点v₅₅₉和v₅₆₀对应的山脊线要素r₅₅₉和r₅₆₀间的差异度Difference₁₈₈。本实例中，取p₁＝0.6，p₂＝0.2，p₃＝0.2，计算得到Difference₁₆₆＝0.6×0.14+0.2×(10.79+10.65)+0.2×7.5＝5.87和并存入集合Sum中；

g)遍历集合Sum，Difference₁₆₆＝5.87＜Difference₁₈₈＝38.06，故对边e'₁₆₆的节点v₅₀₀和v₅₆₀所对应的山脊线要素r₅₀₀和r₅₆₀进行关联处理；

h)循环变量i减1，i＝187≥2，则循环变量j＝187-1＝186，继续执行步骤a)-步骤g)，直至关联结束。本实例中，共关联了154条山脊线要素，关联后的山脊线要素集合LM₁＝{r_i|i＝1，2，…，444}，共对154条边所对应的山脊线要素进行了关联处理；

步骤45：基于第一次关联后山脊线要素集合LM₁＝{r_i|i＝1，2，…，444}，空间剖分结果V₁，如图7所示。循环执行步骤(一)-(四)，本实例中，直至循环6次后该山脊线要素集合LM₆中已不存在符合关联约束条件的邻接山脊线要素对，输出山脊线要素集合LM₆＝{r_i|i＝1，2，…，351}。山脊线要素集合LM₆的空间剖分结果V₆(如图8所示)，山脊线要素长度的符号化表示如图9所示。第2、3、4、5、6次关联处理时分别关联了52条、21条、10条、8条和2条山脊线要素。

Claims (3)

1.一种破碎山脊线的自动关联方法，其特征在于，包括如下步骤：首先，对破碎山脊线进行预处理，并设定关联约束条件；然后，对破碎山脊线进行空间剖分处理，构建山脊线的邻接关系图，即ARG模型；其次，基于关联约束条件，提取符合关联约束条件的山脊线要素对，并进行山脊线要素的关联；最后，对关联处理后的山脊线重构ARG模型，继续提取符合关联约束条件的山脊线要素对，进行关联处理，直至山脊线图层中不存在符合关联约束条件的邻接山脊线要素对，则关联处理结束。

2.根据权利要求1所述的一种破碎山脊线的自动关联方法，其特征在于，方法的具体步骤如下：

(一)破碎山脊线的预处理

步骤11：加载山脊线图层，并将相应的山脊线要素读入集合R₁＝{r_i|i＝1，2，...，n₁}，n₁为山脊线要素的条数；

步骤12：依据山脊线图层信息，设定关联约束条件：设置相对走向约束值SThreshold、相对趋势约束值IThreshold和最短距离约束值DThreshold；

步骤13：剔除非线状要素的山脊线要素，并读入集合R₂＝{r_i|i＝1，2，...，n₂}，n₂为剔除了非线状要素的山脊线的个数；

步骤14：针对存在公共点的两山脊线要素进行去毛刺处理，并读入集合R₃＝{r_i|i＝1，2，...，n₃}，n₃为预处理后的山脊线要素的个数；

(二)山脊线要素的空间剖分处理

步骤21：将山脊线要素集合R₃读入点集P＝{p_i(x_i，y_i)|i＝1，2，...，m}，m为点要素的个数，并使用点集P的“山脊线要素”属性来记录该点所归属的山脊线要素；

步骤22：基于点集P，构建Voronoi图Vor；以图层的最小外包矩形为界，构建点集P的Voronoi图Vor＝{T_i|i＝1，2，...，m′}，T_i表示第i个Thicssen多边形，m'为构成Voronoi图Vor的Thicssen多边形T的个数；并将点集P中各点的“山脊线要素”属性，赋给包含该点的相应Thicssen多边形T_i的“山脊线要素”属性；

步骤23：基于Thicssen多边形T_i的“山脊线要素”属性合并Thicssen多边形T_i，得到每条山脊线要素r_i的影响范围；依次将Voronoi图Vor中具有相同“山脊线要素”属性值的Thicssen多边形T_i合并为一个多边形T'_i，生成图Vor'＝{T'_i|i＝1,2,...,n'}，n'为图Vor'的Thicssen多边形的个数；图Vor'即山脊线要素集合R₃基于最小外包矩形的空间影响范围集合，T'_i代表第i条山脊线要素r_i的影响范围；

(三)基于空间剖分结果，构建山脊线要素的ARG模型G＝(V，E)

步骤31：计算山脊线要素的ARG模型G＝(V，E)中节点集合V＝{v₁，v₂，...，v_n3}及其相关属性；其中，V中每个节点代表一个山脊线要素对象，E为ARG模型的边集合；具体步骤如下：

a)计算山脊线要素r_i的剖分范围T'_i的质心点c_i得到节点v_i，并存入节点集合V中；i为循环变量，i∈[1，n₃]，初值为1；

b)读取山脊线要素r_i的点集C_i＝{c_i1，c_i2，...，c_it}，t为山脊线要素r_i中点的个数；

c)基于点集C_i分别计算直线c_i1c_i2、c_i2c_i3、...、c_it-1c_it的倾斜角α₁₂、α₂₃、...、α_t-1t，并按最大频数法或均值法计算山脊线要素的走向Strike_i，Strike_i∈[0°，180°)，并存入节点v_i的属性中；

d)基于点集C_i，计算山脊线要素r_i的最小凸包的点集H_i＝{h_i1，h_i2，...，h_is}，s为点的个数，并存入节点v_i的属性中；

e)求出山脊线要素r_i的剖分范围T'_i的顶点集Vt_i＝{p_i1，p_i2，...，p_il}，并存入节点v_i的属性，l为顶点集Vt_i中点的个数；

f)若循环变量i＜n₃，则i加1，继续执行步骤a)-e)，直至所有节点及其属性均提取完成；

步骤32：基于节点集合V，计算山脊线的ARG模型G＝(V，E)的边集合w₁为边的个数，具体步骤如下：

a)读取节点集合V中的节点v_i和v_j所对应的影响范围T'_i的顶点集Vt_i和影响范围T'_j的顶点集Vt_j；其中，i为循环变量，i∈[1，n₃-1]，初值为1；j为循环变量，j∈[i+1，n₃]，初值为i+1；

b)判断顶点集Vt_i和Vt_j是否存在连续的两公共点，若存在，则山脊线要素r_i的影响范围T'_i和r_j的影响范围T'_j为邻接拓扑关系，连接节点v_i和v_j，得到ARG模型一条边e＝(v_i，v_j)，称节点v_i与v_j互为邻接节点，山脊线要素r_i和山脊线要素r_j为边e所对应的邻接山脊线要素对，并将边e存入边集合E；否则，跳转步骤c)；

c)若j≤n₃，则循环变量j加1，继续执行步骤a)和b)；否则，循环变量i加1，执行下一步骤d)；

d)若循环变量i≤n₃-1，则循环变量j＝i+1，继续执行步骤a)-c)；否则，ARG模型的边集合已提取结束；

步骤33：计算ARG模型G＝(V，E)的边集合的属性{S，I，D}，S表示该边所对应的邻接山脊线要素对的相对走向关系，I表示邻接山脊线要素对的趋势关系，D表示邻接山脊线要素对的最短距离关系；具体步骤如下：

a)读取边集合E的第k条边e_k，k为循环变量，k∈[1，w₁]，获取节点v_a和节点v_b的走向属性Strike_a和Strike_b；

b)计算由节点v_a与v_b所对应的山脊线要素r_a与r_b相对走向关系S_k＝Strike_a-Strike_b，并存入边e_k的属性中；

c)计算节点v_a和v_b构成的直线v_av_b的倾斜角α_k，α_k∈[0°，180°)，作为边e_k的相对趋势关系I_k，并存入边e_k的属性中；

d)获取边e_k的节点v_a的最小凸包的点集H_a＝{h_a1，h_a2，...，h_ap}和节点v_b的最小凸包的点集H_b＝{h_b1，h_b2，...，h_bq}，p和q均为点的个数；

e)按式(1)计算h_ai与h_bj的距离d_ij，并存入最短距离集合Distance中，其中i和j为循环变量，初值均为1，i∈[1，p]，j∈[1，q]；

$d_{ij}=\sqrt{{(x_{ai}-x_{bj})}^2+{(y_{ai}-y_{bj})}^2}---\left(1\right)$

其中，(x_ai，y_ai)和(x_bj，y_bj)分别为h_ai和h_bj的坐标；

f)若循环变量j≤q，循环变量j加1，继续执行上一步骤e)；否则，循环变量i加1，执行下一步骤g)；

g)若循环变量i≤p，则循环变量j＝1，继续执行步骤e)；否则，距离d_ij计算结束，执行下一步骤h)；

h)读取最短距离集合Distance，取集合Distance中的最小值D_k＝min{d₁₁，d₁₂，...，d_pq}作为边e_k的最短距离关系，并存入边e_k的属性中；

i)读取边集合E的下一条边e_k+1，循环以上步骤，直至边集合E中每条边的相对走向关系S、相对趋势关系I和最短距离关系D均计算结束；

(四)山脊线的关联处理

步骤41：遍历ARG模型的边集合E＝{e₁，e₂，...，e_w1}中的每条边e及属性{S，I，D}，并获取边e的节点v_a和节点v_b的Strike_a和Strike_b；

步骤42：对每条边e均进行关联约束条件的判断，若均满足以下条件：

$\left\{\begin{array}{c}\left|S\right|<SThreahold\\ |{\mathrm{Strike}}_a-I|<Ithreshold\\ |{\mathrm{Strike}}_b-I|<Ithreshold\\ D<DThreshold\end{array}\right.$

则边e符合关联约束条件，节点v_a和v_b所对应的山脊线要素r_a和r_b满足关联约束条件，并将边e记为e'存入边集合E'；若不满足，则跳转下一步骤43；

步骤43：循环执行步骤42，直至所有边e均进行关联约束条件的判断后，得到边集合w₂为边集合E′内的边数；

步骤44：基于边集合进行山脊线要素的关联处理，具体步骤如下：

a)读取边集合E′的第i条边e'_i，并提取边e'_i的节点v_a和v_b，其中i为循环变量，初值为w₂；

b)读取并提取边集合E'的第j条边e'_j的节点v_c和v_d，其中j为循环变量，初值为w₂-1；

c)若节点v_a、v_b、v_c和v_d，均满足v_a≠v_c，v_a≠v_d，v_b≠v_c，v_b≠v_d，则边e'_i与边e'_j间不存在共同关联，循环变量j减1，执行步骤e)；否则执行步骤d)；

d)由于边e'_i和边e'_j存在共同关联，将e'_j存入集合L_i(L_i存储与e'_i存在共同合并的边e')循环变量j减1，并跳转步骤e)；

e)若j≥1，则跳转步骤b)，否则执行下一步骤f)；

f)将边e'_i存入集合L_i，依次读取L_i中边e'，并按式(2)分别计算每条边e'所对应的节点v_e和v_f对应的山脊线要素r_e和r_f的差异度Difference，并存入集合Sum中；

Difference＝p₁×|S|+p₂×(|Strike_e-I|+|Strike_e-I|)+p₃×D (2)

其中，p₁，p₂，p₃分别代表相对走向关系S、相对趋势关系I、最短距离关系D的权值，p₁、p₂、p₃均∈(0，1)之间，且满足p₁+p₂+p₃＝1；

g)遍历集合Sum，并获取值最小的差异度Difference所对应的边e'，对边e'的节点v_e和v_f所对应的山脊线要素r_e和r_f进行关联处理；

h)循环变量i减1，若i≥2，则循环变量j＝i-1，继续执行步骤a)-步骤g)；否则，输出关联后的山脊线要素集合LM₁＝{r_i|i＝1，2...，l₁}，l₁为第一次关联处理后的山脊线的条数；

步骤45：循环执行步骤(一)至(四)，直至循环K次后该山脊线要素集合LM_K不存在符合关联约束条件的邻接山脊线要素对，则连接处理结束，输出山脊线要素集合LM_k＝{r_i|i＝1，2，...，l_K}。

3.根据权利要求2所述的一种破碎山脊线的自动关联方法，其特征在于，步骤14的具体处理过程如下：

a)将山脊线要素r_a和r_b分别读入点集PointA＝{p₁，p₂，...，p_u}和PointB＝{p₁，p₂，...，p_v}，u、v分别为点集PointA、PointB中点的个数；

b)计算山脊线要素r_a和r_b的公共交点，存入点集CrossPoint＝{p₁，p₂，...，p_c}，c为公共交点的个数；

c)若公共交点的个数c＝1，则按均值法或最大频数法计算山脊线要素r_a的走向Strike_a和山脊线要素r_b的走向Strike_b；否则，跳转下面的步骤f)；

d)若|Strike_a-Strike_b|＜SThreshold，则对山脊线要素r_a和r_b进行合并处理；否则，跳转e)；

e)对山脊线要素r_a和r_b进行打断操作，山脊线要素r_a和山脊线要素r_b均删除公共交点CrossPoint；

f)公共交点的个数c≥2，计算山脊线要素r_a的长度length_a和山脊线要素r_b的长度length_b；若length_a＞length_b，则删除长度较短山脊线要素r_b；否则，删除山脊线要素r_a。