一种水平岩层构造地貌的自动识别方法与流程

本发明涉及地理信息技术应用领域，尤其涉及一种基于矢量地形地质图，实现一种水平岩层构造地貌自动化识别的方法。

背景技术：

水平岩层是指沉积成岩后只有整体升降而未经倾伏和褶曲的原始水平产状的地层，也包括经过构造变动，但仍具有近水平产状的地层，如大型平卧褶皱两翼的岩层等，水平岩层构造地貌，主要包括：构造平原，丹霞地貌，方山、尖山等。一般认为，倾角小于5°的岩层即是水平岩层。水平岩层在地质图上表现为地质界线大致平行于地形等高线。水平岩层的准确识别与高效提取，对于三维地质建模应用和水平岩层构造地貌研究，具有较为重要的实用价值和研究意义。然而，传统地质应用中，水平岩层构造的识别主要通过人工读取地形地质图来进行，不仅需要一定的专业基础、读图效率低下，且识别质量因人而异。

技术实现要素：

发明目的：本发明针对现有技术存在的问题，提供一种水平岩层构造地貌自动化识别的方法。

技术方案：本发明所述的水平岩层构造地貌自动化识别的方法，包括以下步骤：

(1)将矢量地形地质图的岩层面要素转换为岩层界线线要素；

(2)将岩层界线线要素转换为离散的点，生成离散点集合；

(3)基于等高线插值生成数字高程模型的基础上，进行离散点高程赋值并计算每条岩层界线所对应点集合的高程属性值；

(4)基于二分法，根据高程属性值迭代筛选出符合判别阈值的水平岩层界线；

(5)根据水平岩层界线的属性值，匹配上、下岩层界线，生成水平岩层面状要素。

其中，步骤(1)具体包括以下步骤：

(1-1)将矢量地形地质图数据读入岩层面要素集合A＝{a_i1|i1＝1,2,...,n_a}，将等高线数据读入等高线集合E＝{e_i2|i2＝1,2,...,n_e}，设置判定标准差的阈值Y1和Y2，长度阈值Len和判定整个岩层是否水平的阈值F；其中，a_i1为包含rockid属性的岩层面要素，n_a为岩层面要素的个数，e_i2为等高线，n_e为等高线的条数，Y1≥Y2，F∈[0.75,1]；

(1-2)将岩层面要素集合离散化为线要素集合：读取岩层边界并预处理，从而将A转换为岩层界线线要素集合L₁＝{l_i|i＝1,2,...,n_l}，l_i为岩层界线线要素，包含界线左岩层要素rockid1和界线右岩层要素rockid2两个属性，n_l为岩层界线的条数。

其中，步骤(2)具体包括以下步骤：

依次读取L₁中的每个岩层界线线要素l_i，将l_i中的每个节点坐标记录到集合p_i＝{(x_ij，y_ij)|j＝1，2，..，len_i}，其中，(x_ij，y_ij)表示l_i的第j个节点坐标，len_i表示l_i所包含节点的个数；最终生成所有岩层界线离散后的点集合P＝{p_i|i＝1，2，..，n_l}；

其中，步骤(3)具体包括以下步骤：

(3-1)将矢量的等高线离散化为数据点，然后采用IDW插值方法内插生成数字高程模型DEM；

(3-2)点集合P中点的高程信息的获取：根据DEM，逐一读取点集合P中节点的坐标信息，将节点的坐标值换算成DEM的栅格行列号其中：

式中，为栅格数据的行列号；X，Y为点集合P中节点的横、纵坐标值；X₀，Y₀为栅格数据的起始坐标，size表示栅格大小。

(3-3)根据获取相应的栅格值，并将栅格值作为相应点的高程属性值记录，则将P转换为具有高程属性值的点集合Q＝{q_i|i＝1，2，..，n_l}，其中，q_i＝{(x_ij，y_ij，d_ij)|j＝1，2，..，len_i}，d_ij为q_i的第j个点的高程属性值。

其中，步骤(4)具体包括以下步骤：

(4-1)分别计算Q的子集q_i的高程属性列的标准差S_i和q_i的长度length_i，其中，

式中，表示q_i平均高程属性值；

(4-2)根据计算结果进行如下判断操作：

5)若S_i＞Y1，则q_i所代表的岩层界线不水平，跳至步骤(4-3)；

6)若S_i≤Y2，q_i转换为空间线要素，和其所属的岩层界线一并存储到集合L₂，跳至步骤(4-3)；

7)若Y2＜S_i≤Y1，且length_i＜Len，跳至步骤(4-3)；

8)若Y2＜S_i≤Y1，且length_i≥Len，采用以下公式将集合q_i分为子集合q_i1和子集合q_i2，并返回执行步骤(4-1)；

(4-3)将判别为水平岩层界线的相邻子集合并；

(4-4)读取Q的下一个子集，重复步骤(4-1)，直至i取尽1,2,..,n_l，处理完所有的子集，最后生成水平岩层界线集合L₂＝{lh_i3|i3＝1,2,...,n_lh}，lh_i3为水平岩层界线，n_lh为水平岩层界线的条数。

其中，步骤(5)具体包括以下步骤：

(5-1)读取岩层面要素集合A的一个元素a_i1，用rn_i1表示a_i1的rockid属性：

5)遍历L₂中所有水平岩层界线，将rockid1属性等于rn_i1的水平岩层界线记录到集合N1_i1中；

6)遍历L₂中所有水平岩层界线，将rockid2属性等于rn_i1的水平岩层界线记录到集合N2_i1中；

7)根据公式计算得到f_i，其中，M表示N1_i1和N2_i1中所有水平岩层界线所包含的节点数，N表示岩层面a_i1所包含的节点数；

8)若N1_i1和N2_i1都不为空且f_i≥F，说明a_i1是水平岩层要素；否则，a_i1是非水平岩层要素；

(5-2)重复步骤(5-1)，直到处理完A中所有元素；

(5-3)将识别出的水平岩层面要素另存为新图层，即为水平岩层构造地貌。

有益效果：本发明与现有技术相比，其显著优点是：本专利基于矢量地形地质图，通过将岩层面转换为岩层界线、将岩层界线转换为离散的点、离散点高程赋值、水平岩层界线提取、水平岩层界线匹配、水平岩层输出等步骤，有效实现了水平岩层构造地貌的自动化识别。

附图说明

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

图2为本实施例采用的地形地质图；

图3为离散化后的点数据示意图；

图4为数字高程模型示意图；

图5为含高程值的点集合数据示意图；

图6为水平岩层界线图；

图7为水平岩层构造识别结果示意图。

具体实施方式

如图1所示，本发明的流程如图1所示，具体包括以下步骤：

(一)岩层面离散化为点集合

步骤11：读取shp格式的地形地质图，将岩层面要素读入集合A，本实施例共有9个岩层面要素；将等高线数据读入集合E，本实施例中共有18条等高线要素；设置阈值：Y1＝11.5、Y2＝11.0、len＝100、F＝0.75；

步骤12：将岩层面要素集合转换成空间线要素集合。根据专利(陈楹等.一种直立岩层的自动化识别方法.中国专利，CN106023197A.2016-10-12)中对岩层边界读取及预处理的步骤，读取岩层边界并预处理，从而将A转换为空间线要素集合L₁。本实施例中共有14条岩层界线要素，如图2所示；

步骤13：将L₁中任一岩层界线数据l_i离散成点集数据p_i，最后，生成点集合P。本实施例中，共有14条岩层界线要素，P包含14个点集，共计1307个点，数据如图3所示。

(二)点的高程赋值

步骤21：将矢量等高线离散化为数据点，然后采用IDW插值方法内插生成DEM。本实施例中，采用IDW插值方法生成的DEM如图4所示；

步骤22：点集合P中点的高程信息的获取，生成Q。本实施例中14条岩层界线要素，共计1307个含高程值的点，其部分数据内容如图5所示。

(三)水平岩层界线的提取

步骤31：读取Q中的一个子集q_i。本实施例中，点集合q₁共包含76个点，部分数据内容亦如图5所示。利用公式，计算q₁的高程属性列的标准差S₁＝1.97506257和q₁的长度length₁＝680.3654，判断可知：S₁＝1.97506257≤Y2，因此，q₁所代表的线ID为0的岩层界线是水平岩层界线，将q₁转换为空间线状要素(含相应的属性数据)存储到集合L₂中。

步骤32：将判别为水平岩层界线的相邻子集合，合并处理。本实施例中，由于不存在相邻子集合，故此处理过程跳过；

步骤33：读取下一个点集，重复步骤31，直到读取完所有的点集。本实施例中，按以上步骤依次对余下的点集进行处理得到线ID为：1,2,3,5,10,12，共六条水平岩层界线，加上之前的线ID为0的水平岩层界线合计七条，因此，L₂＝{lh_i3|i3＝1,2,...,7}，数据如图6所示。

(四)水平岩层的识别与提取

步骤41：读取岩层面要素集合A的一个元素a₁，本实施例中，rn₁＝K1：

1)遍历L₂中所有线要素，将rockid1属性等于rn₁的线的ID记录到集合N1₁中。本实施例中，N1₁包含线ID为：1、2、5共3条界线；

2)遍历L₂中所有线要素，将rockid2属性等于rn₁的线的ID记录到集合N2₁中。本实施例中，N2₁包含线ID为：3，共1条界线；

3)根据公式计算得到f_i。本实施例中，根据公式计算得到f₁＝1.0。

a₁所对应的N1₁，N2₁都不为空且f₁＞0.75，因此，a₁是水平岩层；

步骤42：重复步骤41，直到处理完A中所有数据；

步骤43：将识别出的水平岩层面成图输出。本实施例中，按照上述步骤对A进行处理。最后，识别得到两个符合要求的岩层面要素，分别是a₁和a₃，它们对应的rockid属性都是K1。识别出的岩层面数据内容如图7所示。