一种地表出露岩层产状的自适应判定方法与流程

本发明属于地理信息技术应用领域，具体涉及一种基于地形地质图，自适应判定地表出露岩层产状的方法。

背景技术：

岩层的产状是指岩层在空间的产出状态。除水平岩层成水平状态产出外，一切倾斜岩层的产状均以其走向、倾向和倾角表示，称为岩层产状三要素。产状作为地质测量与成图中常见的地质要素之一，在地质各领域有着广泛应用，其获取方法主要包括野外实测和间接计算两种方法。

野外的岩层产状测量方法，通常是使用地质罗盘、坡度仪等仪器，直接在出露的岩层层面上量测。但是，野外岩石的风化、层面和劈理、线理面区分困难等因素均会影响产状的测量。有时受限于岩层出露情况，只能借鉴附近区域的数据，再通过室内计算方法间接获得产状数据。产状间接计算中常用的方法，主要有“三点法”和“相邻等高线法”，其原理为通过寻找共面且不共线的多点，确定平面的形态。在计算过程中，一方面，需要人工筛选合适的点和计算方法，并通过绘制辅助点和线，对产状进行计算；另一方面，需要人为确定产状符号位置。相关计算工作的人为化处理，导致计算效率低下，且产状数据的质量也因人而异。

为此，本专利针对地表岩层，拟基于地形地质图，实现一种地表岩层产状的自适应计算方法。其重点在于自适应完成交点选择、计算方法选择和产状点位选择。

技术实现要素：

本发明所要解决的技术问题在于，克服现有技术存在的缺陷，提出了一种地表岩层产状的自适应判定方法，基于地形地质图，自适应完成交点选择、计算方法选择和产状点位选择。从而实现地表出露岩层产状的自适应判定。

本专利的主要思路如下：1)提取和预处理岩层边界线，计算岩层界线与等高线的交点；2)遍历同一岩层界线上的邻接交点，通过智能的交点选择、计算方法选择和产状点位选择等步骤，实现产状的自适应判定；3)对产状信息进行符号化和标注。

基于上述思路，本发明提出的一种地表出露岩层产状的自适应判定方法，其步骤如下：

(一)交点集合获取

步骤1.1：加载矢量格式的岩层面要素图层和等高线线图层，得到所有面要素集合Stra＝{s_i|i＝1,2,3,…,n}(其中：s_i为岩层面要素，n为岩层面要素的数量)和等高线的线要素集合CourLine＝{cl_i|i＝1,2,3,…,p}(其中：cl_i为等高线的线要素，p为等高线要素的数量)。

步骤1.2：读取集合Stra中每个岩层面要素的边界，将其预处理为单线要素，并剔除图幅边缘部分、第四纪岩层边界部分、断层界线部分以及不整合接触界线部分，得到岩层层面边界线集合SLine＝{sl_i|i＝1,2,3,…,c}，sl_i是岩层层面边界线要素，c为所有岩层保留的层面出露边界线的数量。

步骤1.3：针对任一边界线要素sl_i遍历等高线集合CourLine，求取线要素sl_i与所有等高线的交点集合InterP_i，对边界线集合SLine内所有元素进行交点计算，得到集合InterP＝InterP₁∪InterP₂∪InterP₃∪…∪InterP_l，l为与等高线相交的边界线数量。记录交点的高程H，源岩层编号SID和源等高线编号CID，在同一集合InterP_i内，交点的SID属性均为i。

(二)产状自适应获取与剔除

针对交点集合InterP_i＝{pi_j|j＝1,2,3,…,qi}，qi为交点集合InterP_i包含的交点数量，逐点进行条件判断，获取产状a＝{ρ，θ，δ}，并确定产状点位。其中：ρ为岩层走向，θ是岩层倾向，δ是岩层倾角。

由j＝1开始进行判断，具体步骤为：

步骤2.1：读取交点pi_j、pi_j+1和pi_j+2和pi_j+3(若下标超出qi，则使用pi_(j+1)％qi、pi_(j+2)％qi和pi_(j+3)％qi)，判别交点pi_j、pi_j+1、pi_j+2和pi_j+3是否满足相邻等高线法的下列选点规则；如果符合，执行步骤2.4，如果不符合执行步骤2.2：

a)交点pi_j和pi_j+3高程一致，交点pi_j+1和pi_j+2高程一致；

b)交点pi_j+1和pi_j+2的源等高线编号属性CID一致；

c)向量和的最小夹角不大于阈值ε，即判定两向量方位近似平行；阈值ε的设定一般不大于6°，两线夹角在此范围内可视为近似平行。阈值的设定越小，则平行要求越严格，可根据需求调整阈值以满足不同要求。

步骤2.2：读取交点pi_j、pi_j+1和pi_j+2(若下标超出qi，则使用pi_(j+1)％qi、pi_(j+2)％qi)，判别是否满足三点法的下列选点规则；如果符合，执行步骤2.6；如果不符合，执行步骤2.3:

a)交点pi_j、pi_j+1和pi_j+2中至少两点的高程不等；

b)交点pi_j、pi_j+1和pi_j+2不共线。

步骤2.3：重新选择判断起点。若原线要素sl_i为环线要素，则取d＝qi为界，若原线要素sl_i非环线要素，则取d＝qi-3为界。下标j增加1之后，若j不超过d，则返回步骤2.1；若j超过d，则交点集合InterP_i上的产状信息获取完成，执行步骤2.8。

步骤2.4：相邻等高线法计算产状。基于交点pi_j(x_j,y_j,h₁)、pi_j+1(x_j+1,y_j+1,h₂)、pi_j+2(x_j+2,y_j+2,h₂)和pi_j+3(x_j+3,y_j+3,h₁)使用相邻等高线法计算产状：

a)对于交点pi_j、pi_j+1、pi_j+2和pi_j+3，计算平面向量和分

别取和方向上的单位向量和取和向量和的方位角度ρ为岩层走向；

b)根据公式(1)计算岩层倾向θ；

c)计算交点pi_j+1和交点pi_j+2的中点，记为辅助点ap₁(x₁,y₁,h₂)，过点ap₁，沿θ方向的

直线与向量v₂交于点ap₂(x₂,y₂,h₁)，根据公式(2)计算倾角δ；

步骤2.5：相邻等高线法产状点位自动确定。记直线l为过辅助点ap₁且方向为θ的直线，l与交点pi_j+1与pi_j+2之间的sl_i边界线部分相交于点ap₃，则以ap₃坐标(x₃,y₃)作为产状标注点位坐标，存储产状信息a。返回步骤2.2。

步骤2.6：三点法计算产状。基于交点pi_j(x₁,y₁,h₁)，pi_j+1(x₂,y₂,h₂)，pi_j+2(x₃,y₃,h₃)使用三点法计算产状：

a)将交点pi_j、pi_j+1和pi_j+2的坐标代入公式(3)，得到岩层层面方程(公式(4)所示，

其中D为任意常数)的系数A、B和C；

Ax+By+Cz+D＝0 (4)

b)根据系数A、B，确定岩层倾向线方程(公式(5)所示)，其中M为任意常数；

Bx+Ay+M＝0 (5)

c)由公式(6)、(7)分别计算、修正倾向θ；

d)根据倾向θ，由公式(8)计算岩层走向；

ρ＝θ±90 (8)

e)根据系数A、B和C，由公式(9)计算岩层倾角δ；

步骤2.7：三点法产状标注点位自动确定。以交点pi_j+1的坐标作为产状点位坐标，存储产状信息a。返回步骤2.3。

步骤2.8：遍历岩层边界线sl_i对应的交点集合InterP_i上的交点，得到Att_i＝{ai_j|j＝1,2,3,…,ci}，其中ci为边界线sl_i上的产状点数量。对边界线集合SLine内所有岩层边界线sl_i进行产状计算，得到产状集合Att＝Att₁∪Att₂∪Att₃∪…Att_s，s为有效获取了产状的交点集合数量。

步骤2.9：对具明显错误的产状值进行剔除。以产状集合Att_i为单位，基于构造特征结合地质知识对错误的产状值进行剔除，得到筛选后的产状集合Att_i＝{ai_j|j＝1,2,3,…,di}，其中di为边界线sl_i上的纠错后保留的产状点数量。错误产状值剔除步骤如下：

a)褶皱具有“两翼产状较稳定，中间部分产状过渡”的构造特征。结合“V字型法则”，整体判断两翼产状倾向。具体表现为：同斜褶皱两翼产状同向，直立向斜两翼产状指向外部，直立背斜两翼产状指向内部。剔除褶皱上不满足上述褶皱构造特征的产状值；

b)若岩层厚度均匀，则两侧产状基本为同向。对于岩层形状简单，两侧界线基本平行的岩层，首先可根据周围构造或结合“V字型法则”大致确定一侧界线上的产状；继而推测另一侧界线上的产状。剔除简单岩层上不满足其产状特征的产状值；

c)形态较平直的界线部分，若界线上产状倾向不稳定，统计界线上产状倾向在不同方向区间上的频率，保留高频并连续的产状区间数据。区间大小根据实际情况而定，在产状个数较少的情况下，可取较大的区间。取一平直界线上的的16个产状数据为例(如下表所示)，以10°为间隔对产状倾向做雷达图(如图2所示)。其产状数据倾向连续且较多分布于90°到130°之间，则保留90°到130°之间的产状数据，剔除其他低频、离散的数据值

步骤3、产状数据符号化

步骤3.1：对于产状点ai_j，过其坐标(axi_j,ayi_j)点，沿其倾向θ绘制单位长度的线，垂直于倾向θ绘制单位长度的线，完成产状符号的绘制。

步骤3.2：将倾向θ由极坐标系转为大地坐标系，对产状要素倾角δ进行标注。

本发明方法针对地表出露岩层，基于地形地质图，自适应完成交点选择、计算方法选择和产状点位选择的产状判定，实现了地表出露岩层产状的自适应判定，实用性强，适用性广。

附图说明

图1本发明方法流程图。

图2示例产状倾向雷达图

图3实施例实验数据。

图4岩层边界线和等高线交点图。

图5示例交点图。

图6示例三点法自适应选点结果图。

图7示例相邻等高线法自适应选点结果图。

图8实施例全图三点法自适应选点结果图。

图9实施例全图相邻等高线法自适应选点结果图。

图10实施例全图产状符号化和标注结果图。

图11示例实测产状数据与计算产状数据对比图。

具体实施方式

下面结合实施例和附图，对本发明方法作进一步详细说明。

实施例：选用某地1:50000地形地质图(如图2所示)数据，下面结合附图，基于shp格式的地质体面图层和等高线图层，描述本方法方法实施过程并进一步说明本发明的效果。

(一)交点集合获取

步骤1.1：加载岩层面要素图层和等高线图层，得到所有面要素集合Stra＝{s_i|i＝1,2,3,…,29}，等高线要素集合CourLine＝{cl_i|i＝1,2,3,…,21}。

步骤1.2：读取集合Stra中岩层面要素的边界线分别进行预处理，去除第四纪地层边界线及地层与第四纪接触的边界线部分，并去除图幅边界的部分，得到边界线集合SLine＝{sl_i|i＝1,2,3,…,47}。

步骤1.3：对边界线集合SLine和等高线集合CourLine求取交点，边界线要素sli上的交点存入一个集合InterP_i，得到交点集合InterP₁＝{(400875.27,3271721.62),(400957.22,3271801.49),(401056.83,3271892.77)}，InterP₂＝{(401235.55,3272047.50),(401334.05,3272128.06),(401434.51,3272198.33),(401551.99,3272129.90),(401497.99,3272039.65),(401398.55,3272039.69),(401376.85,3271795.00),(401606.74,3271790.62)}，…，InterP₄₄＝{(400144.21,3274334.45),(400102.31,3274263.53),(400051.63,3274186.33),(400050.79,3273975.12),(400131.35,3274013.82),(400196.16,3274044.96),(400278.58,3274088.92),(400354.36,3274135.09)}，交点记录等高线高程H，源岩层编号SID和源等高线编号CID，在同一集合InterP_i内，交点的SID属性一致为i。所有交点集合InterP＝InterP₁∪InterP₂∪InterP₃∪…∪InterP₄₄，包括交点398个(如图3所示)。

(二)产状自适应计算与剔除

针对交点集合自适应选点、计算及点位确定，在此设定向量平行判断的角度阈值ε为6°。以岩层面要素s₇的边界线为例，其边界包括岩层边界线sl₁₁、sl₁₂、sl₁₃和sl₁₄，与等高线分别存在16个、3个、0个和9个交点，则边界线sl₁₁、sl₁₂和sl₁₄分别对应交点集合InterP₉＝{p9_j|j＝1,2,3,…,16}、InterP₁₀＝{p10_j|j＝1,2,3}和InterP₁₁＝{p11_j|j＝1,2,3,…,9}(如图4所示)。以交点集合InterP₉为例，交点坐标值如下表所示。由交点p9₁开始进行判断和计算；

步骤2.1：由交点p9₁开始，先对交点p9₁、p9₂、p9₃和p9₄进行相邻等高线法的选点规则判断，并不符合相邻等高线法的选点规则；则对交点p9₁、p9₂和p9₃进行三点法的选点规则判断，并不符合三点法的选点规则，则继续判断交点p9₂。

步骤2.2：由交点p9₂开始，对交点p9₂、p9₃、p9₄和p9₅进行相邻等高线法的选点规则判断，并不符合相邻等高线法的选点规则；则对交点p9₂、p9₃和p9₄进行三点法的选点规则判断，三点满足如下条件：

1)三点的高程分别为250、250和300，高程不全相等；

2)三点不共线；

因此，三点满足三点法选点规则，进行三点法计算。

步骤2.3：将交点p9₂、p9₆和p9₇代入公式计算，得到所在平面方程的系数：

计算得到产状结果：

步骤2.4：产状符号位置取交点p9₃的坐标，继续判断交点p9₃。

步骤2.5：由交点p9₃开始，对交点p9₃、p9₄、p9₅和p9₆进行相邻等高线法的选点规则判断，并不符合相邻等高线法的选点规则；则对交点p9₃、p9₄和p9₅进行三点法的选点规则判断，三点满足三点法选点规则，进行三点法计算，得到产状结果：

产状符号位置取交点p9₄的坐标，继续判断交点p9₄。

由交点p9₄开始，对交点p9₄、p9₅、p9₆和p9₇进行相邻等高线法的选点规则判断，并不符合；则对交点p9₄、p9₅和p9₆进行三点法的选点规则判断，三点满足三点法选点规则，进行三点法计算，得到产状结果：

产状符号位置取交点p9₅的坐标，继续判断交点p9₅。

步骤2.6：由交点p9₅开始，对交点p9₅、p9₆、p9₇和p9₈进行相邻等高线法的选点规则判断，四点满足：

1)交点p9₅和交点p9₈高程一致为350m；交点p9₆和交点p9₇高程一致为400m；

2)交点p9₆和交点p9₇的源等高线编号属性CID一致；

3)交点p9₅指向交点p9₈的方位角为214.75°，交点p9₆指向交点p9₇的方位角为217.06°，两方位角相差2.31°，小于阈值ε。

因此，四点满足相邻等高线法，利用相邻等高线法进行产状计算。计算得到产状结果：

步骤2.7：计算得到产状符号位置坐标(400450.92,3272412.92)，继续判断交点p9₅、p9₆和p9₇是否满足三点法选点规则，直至遍历完成。

步骤2.8：交点集InterP₉产状计算完成，得到边界线sl₁₁的产状集合Att₆＝{a6_j|j＝1,2,3,…,15}。其中，使用三点法计算得到的有13个，其三点组合形式如图5所示；使用相邻等高线法计算得到的有2个，其四点组合形式如图6所示。

步骤2.9：对所有InterP_i计算产状，除去因不符合条件而没有计算产状值的交点集，最后得到Att＝Att₁∪Att₂∪…∪Att₂₆。考虑位置重复且产状值一致的产状点为同一个，则共包含产状点151个。其中，使用三点法计算得到的有132个，其三点组合形式如图7所示；使用相邻等高线法计算得到的有19个，其四点组合形式如图8所示。通过对错误的产状计算值进行剔除，得到114个产状计算值(如图9所示)。

(三)产状数据符号化

加载实测产状数据与计算产状数据进行对比(如图10所示)，实测数据包含32个实测产状值。取实测值300m以内的计算值，或与实测值在同一岩层界线上的计算值与实测值作比较，对比结果如下表：

基于上表计算产状值与实测产状值的误差值，将倾向、倾角划分为“正确”、“尚可”和“错误”三个等级，统计其数量。

通过32个产状实测值与其附近的计算值的对比，存在倾向或倾角计算结果误差较大的产状点有5个，倾向倾角完全计算错误的点位0个，其余数据误差在可接受范围内，总体正确率为84.4％。