一种层状岩质边坡类型的自动划分方法与流程

本发明属于地理信息技术应用技术，具体涉及一种层状岩质边坡类型的自动划分方法。

背景技术：

边坡是自然或人工形成的斜坡，不仅是人类工程活动中最基本的地质环境之一，还是工程建设中最常见的工程形式。依据不同的分类标准，可将边坡划分为不同的类型：按地层岩性，可分为土质边坡和岩质边坡；按岩层结构，可分为整体状、块状、层状、碎裂状和散体状结构5大类型等等。而在各种类型的岩质边坡中，层状岩质边坡被认为是分布最广，结构形式最简单的一类边坡，因此对其开展的相关研究也较多。

对层状岩质边坡稳定性影响较大的边坡结构参数，主要包括：边坡与岩层走向之间的夹角、倾向之间的夹角以及岩层倾角等。据边坡与岩层走向的关系，首先，可将层状岩质边坡划分成垂向坡(走向垂直)、切向坡(走向斜交)、近水平坡(走向平行)；进一步，据边坡与岩层倾向的关系，可把近水平边坡划分为水平坡(岩层倾角约0°)、顺向坡(倾向相同)、反倾坡(倾向相反)。综上所述，层状岩质边坡大致可划分为以下5种类型：①水平坡；②顺向坡；③切向坡；④垂向坡；⑤反倾坡。

自动化划分层状岩质边坡类型的前提，是获取边坡边界数据及岩层产状数据。边坡边界数据可基于dem数据源，通过提取山体边界生成(李梦圆,李安波,解宪丽,等.一种山体边界自动提取方法[p].中国:201710049100.7)；而岩层产状数据可基于地形地质图，自适应完成岩层产状的计算(李安波,陈楹,姚蒙蒙,等.一种地表出露岩层产状的自适应判定方法[p].中国:201610912428.2)。

层状岩质边坡的类型划分，是正确分析其变形破坏机理、进行稳定性评价及采用合理的支护设计方法的前提条件。

然而，现有边坡类型划分的相关研究，通常采用人工或半自动化的划分方法，费时费力，且自动化程度低，难以满足大范围边坡类型划分的需要。

技术实现要素：

发明目的：本发明的目的在于解决现有技术中存在的不足，提供一种层状岩质边坡类型的自动划分方法，本发明基于山体边界面数据、山脊线数据及产状数据，进而实现层状岩质边坡类型的自动划分。

技术方案：本发明的一种层状岩质边坡类型的自动划分方法，具体依次包括以下步骤：

(1)预处理山脊线，再基于山脊线分割山体边界，以形成边坡单元；

(2)计算边坡产状与岩层产状信息；

(3)根据层状岩质边坡划分标准，对边坡进行分类。

进一步的，所述步骤(1)具体包括以下过程：

(1-1)加载山体边界面图层、山脊线图层和产状点图层，分别得到山体边界面要素集合mount＝{mi|i＝1,2,…,n}、山脊线线要素集合ridgel＝{rli|i＝1,2,…,a}和产状点要素集合occur＝{oi|i＝1,2,…,b}；

其中，mi表示第i个面要素，n为面要素的个数，rli表示第i个山脊线线要素，a为山脊线线要素的个数，oi表示第i个产状点要素，其属性包括倾向[tendency]和倾角[dipangle]，b为产状点要素的个数；

(1-2)计算山脊线走向：针对每一个山脊线线要素rli，将其读入点集pi＝{pj|j＝1,2,…,t}，t为点要素的个数，依次取相邻两点pj和pj+1，计算pj对象点指向pj+1对象点的方位角，得到方位角集合azimuth＝{azi|i＝1,2,…,t-1}，azi∈[0°,180°)；再按最大频数法筛选方位角：

(1-3)两端延伸山脊线；

(1-4)分割山体成边坡：针对每一个面要素mi，查找其与线要素rli交点前一点的下标e1和e2：此处面要素和线要素的交点即是指rli首尾两点p1和pt；将面要素mi读入点集ci＝{cj|j＝1,2,…,s}，依次取相邻两点cj、cj+1，若首点p1在线段上，则e1＝j；若尾点pt在线段上，则e2＝j；再比较e1与e2的大小，s为点要素的个数。

进一步的，所述步骤(1-2)中筛选方位角的方法为：

(1-2-1)通过参数t及360°来计算方位角间隔n，n＝360°/t；其中，t为360的约数(t常用取值为36或72)；

(1-2-2)将集合azimuth中的方位角分布在各均匀区间[0°,n)、[n,2n)、[2n,3n)、…、[(t-1)×n,360°)；

(1-2-3)计算频数最大的区间内的方位角的均值，并将该均值作为rli的大致走向φi，φi∈[0°,180°)。

进一步的，所述步骤(1-3)中包括以下2个环节：

(1-3-1)延伸山脊线尾端：依次取线要素rli尾端三个点点pt-2、pt-1和pt，计算相邻两点间距离和的均值d；再倒序遍历点集pi，依次取两点点pj和pj-1，其中j为循环变量，初值为t，做出如下判断：

a)若这两点均不在山体边界面要素mi内，则移除点pj；

b)若仅点pj-1在面要素mi内，则移除点pj，且在j+1位置处将线段与面要素mi的交点插入点集pi，并终止循环；

c)若这两点均在面要素mi内，则按公式(1)计算新增点的坐标p'(x,y)，若点p'在面要素mi内，将点p'加入点集pi，j＝j+1；反之，将线段与面要素mi的交点加入点集pi，并终止循环；

其中，pj.x和pj.y分别为点pj的横坐标和纵坐标，φi为线要素rli的大致走向；

(1-3-2)延伸山脊线首端：步骤同延伸山脊线尾端的方法即步骤(1-3-1)。

进一步的，所述步骤(1-4)中e1与e2的大小比较方法为：

(1-4-1)若e1<e2，遍历点集ci，取点cj，进行以下判断及操作：

i、若j<e1或j>e2，将点cj加入点集rslope；

ii、若j＝e1，将点cj及点集pi中的所有点加入点集rslope；

iii、若e1<j<e2，将点cj加入点集lslope；

iv、若j＝e2，将点cj及点集pi中的逆序点加入点集lslope；

(1-4-2)若e1>e2，处理步骤同(1-4-1)；

点集ci遍历完毕，得到左侧坡点集lslope＝{lsi|i＝1,2,…,q}和右侧坡点集rslope＝{rsi|i＝1,2,…,r}，其中，q为点要素的个数，r为点要素的个数。

进一步的，所述步骤(2)的具体过程为：

(2-1)计算边坡产状：山脊线的走向作为左、右侧坡的大致走向，即ρl'＝ρr'＝φi；而左、右侧坡的倾向可按公式(2)计算：

其中，ρl'和θl'为左侧坡的走向和倾向，ρr'和θr'是指右侧坡的走向和倾向；

(2-2)计算岩层产状：

(2-2-1)左侧坡的岩层产状计算步骤如下：

a、遍历产状点要素集合occur，筛选在左侧坡点集lslope内的所有产状点；

b、读取左侧坡内产状点的倾向，按最大频数法查找频数最大的区间内的倾向，再进一步查找对应的产状点，得到左侧坡的优势产状点集loccp＝{lpi|i＝1,2,…,h}，h为优势产状点的个数；

c、遍历点集loccp，计算所有点的倾向均值、倾角均值分别作为左侧坡岩层的倾向θl、倾角δl；按公式(3)计算左侧坡岩层的走向ρl，ρl∈[0°,180°)；

(2-2-2)计算右侧坡岩层的产状(ρr,θr,δr)，步骤同(2-2-1)，ρr、θr和δr分别是右侧坡岩层的走向、倾向和倾角。

进一步的，所述步骤(3)具体包括：

(3-1)划分层状岩质边坡：

(3-1-1)设置左侧坡的类型属性：分别计算左侧坡与左侧坡岩层走向差的绝对值α、倾向差的绝对值β，若β>180°，则β＝360°-β；设走向一致性容差为u，倾向相同或相反容差为w，按如下判断条件设置左侧坡类型：

i、若60°<α<120°，则记类型属性type为“垂向坡”；

ii、若u≤α≤60°或120°≤α≤180°-u，则记类型属性type为“切向坡”；

iii、若0°≤α<u或180°-u<α<180°，可分为3种情况：当左侧坡倾角满足0°≤δl≤5°时，则记类型属性type为“水平坡”；当倾向差绝对值β满足0°≤β<w时，则记类型属性type为“顺向坡”；当β满足180°-w<β≤180°时，则记类型属性type为“反倾坡”；

(3-1-2)右侧坡的类型属性设置步骤同(3-1-1)，直至所有左、右侧坡的类型划分完成；

(3-2)将划分后的边坡面要素另存为新图层。

有益效果：本发明基于山体边界面数据、山脊线数据及岩层产状数据，结合边坡划分标准，获取边坡产状与岩层产状信息，实现了层状岩质边坡类型的自动划分；相比现有的人工-半自动划分方法，自动化程度高，客观性强，能够满足大范围层状岩质边坡的类型划分需要。

附图说明

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

图2为实施例中山体边界面图层与山脊线图层示意图；

图3为实施例中地层线图层与产状点图层示意图；

图4为实施例中延伸处理后的山脊线图层示意图；

图5为实施例中分割山体边界后形成的边坡面图层示意图；

图6为实施例中划分类型的边坡面图层示意图。

具体实施方式

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

实施例：

本实例的实验数据源于庐山北部山区dem数据，格网大小为2162×2260。基于dem数据源，提取山脊线(如图2)、山谷线及等高线数据；再以反地形dem提取盆地边界面数据，然后基于山脊线和山谷线自动提取山体边界面数据(如图2)；最终结合等高线与地质体面数据，计算岩层产状信息(如图3)。下面结合附图，并通过描述一个自动划分层状岩质边坡类型的实例，来进一步说明本发明的效果。

如图1所示，本实施例中，层状岩质边坡类型的自动划分方法包括以下步骤：

步骤(1)，具体包括：

(1-1)如图2和图3所示，加载山体边界面图层、山脊线图层和产状点图层，分别得到山体边界面要素集合mount＝{mi|i＝1,2,…,8}、山脊线线要素集合ridgel＝{rli|i＝1,2,…,8}和产状点要素集合occur＝{oi|i＝1,2,…,507}；

(1-2)计算山脊线走向：以线要素rl1为例，将其读入点集p1＝{pj|j＝1,2,…,924}，依次取相邻两点p1、p2，计算p1对象点指向p2对象点的方位角为68°；再取相邻两点p2、p3，计算p2对象点指向p3对象点的方位角为73°；……；待点集p1最后两点p923、p924的方位角计算完毕，得到方位角集合azimuth＝{68°,73°,…,44°}(方位角个数为923)；再按最大频数法筛选方位角：

1)通过参数t＝36和360°来计算方位角间隔为：360°/36＝10°；

2)将集合azimuth中的方位角分布在各均匀区间[0°,10°)、[10°,20°)、[20°,30°)、…、[350°,360°)；

3)统计频数最大的区间为[50°,60°)(频数为113)，取该区间内方位角的均值54°作为rl1的大致走向；

(1-3)两端延伸山脊线：该步骤包括以下2个环节：

①延伸山脊线尾端：取线要素rl1尾端3点p922、p923和p924，计算两点间距离的均值为7.71205216370225；再倒序遍历点集p1，取两点p924、p923，两点均在面要素m1内，则按公式(1)计算新增点的坐标p'：

p'(405511.09571363218,3278760.43547032245)；

其中，p'在面要素m1内，将p'加入点集p1，j＝925；再取点集p1尾端两点p925、p924，两点均在面要素m1，则计算新增点的坐标，该点不在面要素m1内，将该点与点p925构成的线段与面要素m1的交点加入点集p1，终止循环；

②延伸山脊线首端：取线要素rl1首端3点p1、p2和p3，计算两点间距离的均值为5.9827602150087351；遍历点集p1，取两点p1、p2，两点均在面要素m1内，则按公式(1)计算新增点的坐标p'：

p'(401633.04610488308,3273824.3561243751)；

其中，p'在面要素m1内，将p'插入点集p1首端，j不变；再取点集p1首端两点p1、p2，两点均在面要素m1内，则计算新增点的坐标，该点在面要素m1内，将该点插入点集p1首端，j不变；继续上述操作，直至第7次取点集p1首端两点时，两点均在面要素m1内，但新增点不在面要素m1内，将新增点与点p1构成的线段与面要素m1的交点插入p1首端，终止循环；

延伸操作完成后，线要素rl1包含点集p1的点要素个数为933。待集合ridgel中的其余线要素按步骤(1-3)执行完毕后，生成延伸处理后的山脊线线要素图层，如图4所示；

(1-4)分割山体成边坡：针对面要素m1，查找其与线要素rl1交点(即rl1首尾点p1、p933)前一点的下标e1、e2：将面要素m1读入点集c1＝{cj|j＝1,2,…,2078}，依次取相邻两点，当相邻两点为c29、c30时，尾点p933在线段上，则e2＝29；当相邻两点为c880、c881时，首点p1在线段上，则e1＝880；比较e1与e2的大小可知：

e1>e2，遍历点集c1，将c1-c28这些点加入点集lslope；然后将点c29加入点集lslope，同时将点集p1中的点逆序加入点集lslope；再将c30-c879这些点加入点集rslope；再次将c880加入点集rslope，同时将点集p1中的点依次加入点集rslope；最后将c881-c2078这些点加入点集lslope；

点集c1遍历完毕，得到左侧坡点集lslope(点要素个数为2193)和右侧坡点集rslope(点要素个数为1751)；

面要素集合mount中的其余面要素均按步骤(1-4)进行山体的分割，待全部划分左、右侧边坡后，生成边坡面要素图层，如图5所示。

步骤(2)，具体包括：

(2-1)计算边坡产状：山脊线的走向可作为左、右侧坡的大致走向。针对

山体边界面要素m1及山脊线线要素rl1，其对应左、右侧坡的走向为ρl'＝ρr'＝φ1＝

54°，而左、右侧坡的倾向可按公式(2)计算可得θl'＝144°、θr'＝324°；其余山体

边界面要素对应的左、右侧坡的产状均按上述方法计算；

(2-2)计算岩层产状：该步骤包括以下2个环节：

①左侧坡的岩层产状计算步骤如下：

a、遍历产状点要素集合occur，筛选在左侧坡点集lslope内的所有产状点(共167个产状点)；

b、读取左侧坡内产状点的倾向，按最大频数法查找频数最大的区间内的倾向，再进一步查找对应的产状点，得到左侧坡的优势产状点集loccp＝{lp1,lp2,…,lp32}；

c、遍历点集loccp，计算所有点的倾向均值、倾角均值分别作为左侧坡岩层的倾向θl＝313°、倾角δl＝75°；按公式(3)计算左侧坡岩层的走向ρl＝43°；

②右侧坡的岩层产状计算步骤如下：

a、遍历产状点要素集合occur，筛选在右侧坡点集rslope内的所有产状点(共69个产状点)；

b、读取右侧坡内产状点的倾向，按最大频数法查找频数最大的区间内的倾向，再进一步查找对应的产状点，得到右侧坡的优势产状点集roccp＝{rp1,rp2,…,rp9}；

c、遍历点集roccp，计算所有点的倾向均值、倾角均值分别作为右侧坡岩层的倾向θr＝114°、倾角δr＝80°；按公式(3)计算左侧坡岩层的走向ρr＝24°；

按步骤(2-2)继续计算其余左、右侧坡的岩层产状，直至所有边坡的岩层产状计算完毕。

步骤(3)，具体包括：

(3-1)划分层状岩质边坡：该步骤包括以下2个环节：

①设置左侧坡的类型属性：分别计算左侧坡与左侧坡岩层走向差的绝对值α＝11°、倾向差的绝对值β＝169°；设走向一致性容差、倾向相同或相反容差均为10°，由判断条件知：10°<α<60°，满足条件ii，故左侧坡的类型属性type＝“切向坡”；

②设置右侧坡的类型属性：分别计算右侧坡与右侧坡岩层走向差的绝对值α＝30°、倾向差的绝对值β＝210°，β>180°，则β＝360°-210°＝150°，由判断条件知：10°<α<60°，满足条件ii，故右侧坡的类型属性type＝“切向坡”；

继续执行上述操作来设置其余左、右侧坡的类型，直至完成所有边坡的类型属性赋值；

(3-2)如图6所示，将划分后的边坡面要素另存为新图层。

通过上述实施例可以看出，本专利结合边坡划分标准，并以边坡结构参数为切入点，从地学应用角度出发，形成一种层状岩质边坡类型的自动化划分方法，高效精准且自动划分程度高。