CN113539051B - 基于地质图的地层界线逐点岩层产状获取方法及装置 - Google Patents

Abstract

本发明公开了一种基于地质图的地层界线逐点岩层产状获取方法及装置，具体处理步骤包括：首先，读取面状地层图层、已知产状点数据和DEM数据，根据包含关系匹配地层与产状点，并由DEM提取等高线；提取等高线数据和地层界线交点，用三点法计算交点的产状，并将各个交点产状赋值到最邻近边界点上；然后，判断各个相邻地层间连续关系和接触关系，并基于地层间的相互关系，将已知产状信息赋值到相关边界点上，每个地层基于产状修正的边界点，重新线性插值计算边界点的产状；最后，每个地层基于所有非插值获取产状的边界点，重新线性插值计算边界点产状，并完成地质图地层边界点产状数据生成。本发明对计算数据的依赖较少，自动化程度和计算精度较高。

Description

基于地质图的地层界线逐点岩层产状获取方法及装置

技术领域

本发明属于地理信息技术应用领域，尤其涉及一种基于地质图的地层界线逐点岩层产状获取方法及装置。

背景技术

岩层的产状是指岩层在空间的产出状态。除水平岩层成水平状态产出外，一切倾斜岩层的产状均以其走向、倾向和倾角表示，称为岩层产状三要素。产状作为地质测量与成图中常见的地质要素之一，在地层信息表达、地质剖面制图、三维地质建模等应用中具有重要作用。

岩层产状的获取方法主要包括野外实测和间接计算两种方法。野外的岩层产状测量方法，通常是使用地质罗盘、坡度仪等仪器，直接在出露的岩层层面上量测。但是，野外岩石的风化、层面劈理和线理面区分困难、岩层出露情况等因素均会影响产状的直接测量，且野外量测工作量较大。产状间接计算中常用的方法，主要有“三点法”和“相邻等高线法”，或者综合运用两种方法计算，其原理为通过寻找共面且不共线的多点，确定平面的形态。然而，该方法往往只能计算个别地层边界点上产状信息，且计算出产状的边界点可能在地层界线上分布不均。

为此，针对地表出露岩层，本专利拟基于地形地质图，研究实现一种基于地质图的地层界线逐点岩层产状获取方法。

发明内容

发明目的：本发明针对现有技术存在的问题，提供一种基于地质图的地层界线处岩层产状获取方法。

技术方案：本文所述的基于地质图的地层界线逐点岩层产状获取方法具体包括：

(1)读取岩层的地层面图层、DEM数据和已知产状点数据，得到地层集合Stra、地层界线集合SLine、产状点集合OPoint以及表示地层和产状点关系的地层产状点关系表SOR；

(2)判断地层集合Stra中各相邻地层间相互关系，并基于地层间的相互关系，将实测产状信息赋值到相关边界点上；

(3)提取地层界线集合SLine中所有地层界线与等高线的交点，自适应地采用三点法或四点法计算地层界线的产状，并将产状赋值到对应的边界点上；

(4)每个地层基于产状已知的边界点，将其地层界线分为若干分段，得到地层界线分段集合SLine′；

(5)对于SLine′中每个地层界线分段，根据两端边界点线性插值计算分段内部边界点的产状，存入最终边界点产状集合AT′；

(6)根据边界点的最终产状对每个地层界线分段进行识别，若识别到地层界线分段上存在异常产状边界点，则利用曲率法近似计算其上边界点真实倾向，更新最终边界点产状集合AT′，并基于异常产状边界点划分当前地层界线分段，再转到步骤(5)，直到所有地层界线分段上都不存在异常产状边界点；

(7)根据最终边界点产状集合AT′生成地质图地层边界点产状数据。

进一步的，步骤(1)包括：

(1-1)读取矢量格式的地层面图层，得到地层集合Stra＝{s_i|i＝1,2,3,…,sn}，其中s_i为第i个地层，sn为地层的数量；

(1-2)读取集合Stra中每个地层的边界线，并将其处理为单线，存入地层界线集合SLine＝{sl_i|i＝1,2,3,…,sn}，sl_i是s_i的地层界线；

(1-3)读取DEM数据，从中提取等高线存入等高线集合CourLine＝{cl_j|k＝1,2,3,…,cn},其中cl_j为第k条等高线素，cn为等高线数量；

(1-4)设定步长λ，对集合SLine中每个地层界线，获取其上所有边界点，并从DEM提取每一点所对应的高程值，存入到边界点集合SPoint＝{sp_ij(x_ij,y_ij,z_ij)|i＝1,2,3,…,sn，j＝1,2,3,…,pn}，其中sp_ij是地层界线sl_i的第j个边界点，x_ij,y_ij,z_ij分别为该边界点的横、纵坐标和高程值，pn为地层界线sl_i上边界点个数；之后在地层界线sl_i上曲率大于预设阈值的部位提取若干特征点，也作为边界点按顺序加入集合SPoint中；

(1-5)读取已知产状点数据，存入产状点集合OPoint＝{op_l|l＝1，2,...,on},其中op_l为第l个产状点，on代表产状点数量；

(1-6)根据包含关系，将每个产状点匹配到相应地层，存储到地层产状点关系表SOR＝{so_l(s_i,op_l)}中，so_l为s_i和op_l的关系。

进一步的，步骤(2)包括：

(2-1)从地层集合Stra中获取两个相邻地层s_i和s_i′，若两者皆是沉积岩，转到步骤(2-2)，否则判定相邻地层s_i和s_i′的关系为无显著关系，转到步骤(2-3)；

(2-2)判断s_i和s_i′相对年代属性之差，若为1，判定相邻地层s_i和s_i′的关系为两地层连续且为整合接触关系，若大于1，则判定相邻地层s_i和s_i′的关系为不整合接触关系；

(2-3)将地层间关系存入地层间关系表SSR＝{ssr_i,i′}中，ssr_i,i′表示地层s_i和s_i′的关系；

(2-4)重复步骤(2-1)到(2-3)，直到完成所有相邻地层间关系的判断；

(2-5)从产状点集合OPoint中获取任一产状点op_l，经过op_l，沿op_l的倾向做直线el，从op_l出发，el两端遇到不整合接触的两个地层时被截断；

(2-6)el在op_l倾向方向与op_l所在地层s′_i有一交点，取地层s′_i的地层界线上距离该交点最近的一个边界点，将op_l的产状赋予该边界点，记录在参考产状集合AT中；

(2-7)获取el上一个不在s′_i上的交点，若交点在op_l倾向方向上，且同地层的另一个交点不在op_l倾向方向上，则转到步骤(2-8)；若交点和同地层的另一个交点都不在op_l倾向方向上，且该交点离op_l更近，则同样转到步骤(2-8)，否则转到(2-10)；

(2-8)若交点两侧地层为整合接触关系，则转到(2-9)；若为不整合接触关系，且交点两侧地层中，离op_l近的地层年代更大，则转到(2-9)，否则转到(2-10)；

(2-9)取地层s′_i的地层界线上距离该交点最近的一个边界点，将op_l的产状赋予该边界点，记录在参考产状集合AT中；

(2-10)重复步骤(2-5)-(2-9)，直到所有的产状点对边界点赋值完毕。

进一步的，步骤(3)包括：

(3-1)从边界点集合SPoint中读取任一地层界线sl_i的边界点，获取所有等高线与sl_i边界点的交点存入集合CSP＝{csp_m(x_m,y_m,h_m)|m＝1,2,...,csn}，其中m代表交点序号，x_m,y_m,h_m分别为该交点的横、纵坐标和高程值，csn代表交点个数；

(3-2)从集合CSP中获取任意连续四个交点csp_m-1、csp_m、csp_m+1和csp_m+2，判别这四个交点是否满足相邻等高线法的下列选点规则；

a)交点csp_m-1和csp_m+2高程一致，交点csp_m和csp_m+1高程一致；

b)csp_m和csp_m+1位于同一条等高线；

c)向量

和

的最小夹角不大于阈值ε，即判定两向量方位近似平行；

如果符合，执行步骤(3-4)，如果不符合执行步骤(3-3)：

(3-3)判断交点csp_m-1、csp_m和csp_m+1是否满足三点法的下列选点规则；

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

b)交点csp_m-1、csp_m和csp_m+1不共线；

如果符合，执行步骤(3-5)；如果不符合，执行步骤(3-6)：

(3-4)基于交点csp_m-1、csp_m、csp_m+1和csp_m+2使用四点法计算边界点产状，执行步骤(3-6)；

(3-5)基于交点csp_m-1、csp_m、csp_m+1使用三点法计算边界点产状，执行步骤(3-6)；

(3-6)循环执行步骤(3-2)-(3-5)，直到集合CSP中所有连续交点被遍历完，得到sl_i所有边界点的产状；

(3-7)循环执行步骤(3-1)-(3-6)，直到地层界线集合SLine中所有地层界线被遍历完，得到所有地层界线的所有边界点的产状。

进一步的，步骤(3-4)包括：

(3-4-1)对于交点csp_m-1、csp_m、csp_m+1和csp_m+2，计算平面向量

和

取

和

向量和的方位角度ρ为走向；

(3-4-2)根据下式计算倾向θ；

(3-4-3)计算交点csp_m-1和交点csp_m+2的中点，记为辅助点ap_m1(x′_m,y′_m，h′_m)，过点ap_m1沿θ方向作直线L，L与向量

交于点ap_m2(x″_m,y″_m,h″_m)，与地层界线sl_i交于点ap_m3(x″′_m,y″′_m,h″′_m)，根据下式计算倾角δ；

(3-4-4)将倾向θ、倾角δ作为ap_m3同一地层界线上最邻近的边界点的产状，存储到参考产状集合AT＝{α_ij(θ,δ)}中，α_ij表示地层界线sl_i的第j个边界点的产状。

进一步的，步骤(3-5)包括：

(3-5-1)根据下式得到局部岩层层面方程Ax+By+Cz+D＝0的系数A、B和C，其中{A，B,C}即为三个边界点所确定三角面的法线，D为任意常数；

(3-5-2)根据系数A、B，得到岩层倾向线方程，其中M为任意实数；

Bx+Ay+M＝0

(3-5-3)采用下式计算倾向θ；

(3-5-4)采用下式计算倾角δ；

(3-5-5)将倾向θ、倾角δ作为csp_m同一地层界线上最邻近的边界点的产状存储到参考产状集合AT中。

进一步的，步骤(4)包括：

(4-1)以产状已知的边界点为分割点，将地层界线集合SLine中所有地层界线划分为若干地层界线分段，存入地层界线分段集合SLine′＝{sl_ir|i＝1，2,3,…,sn，r＝1,2,3,…,ln_i′}，sl_ir表示sl_i第i个分段，sn为地层的数量，sl_r为地层分段要素，ln_i′为sl_i上分段数；

(4-2)按照地层界线分段，将地层边界点集合划分为若干个连续的边界点子集合SPoint′，每个子集合中，第一个边界点前和最后一个边界点后的两个边界点为产状已知的边界点。

进一步的，步骤(5)包括：

(5-1)获取SLine′中任一地层界线分段对应的边界点子集合SPoint′，基于第一个边界点前和最后一个边界点后的两个边界点的产状，对每个边界点子集合SPoint′中的边界点产状线性插值；

(5-2)按照边界点和产状对应关系，将插值后的边界点产状存入最终边界点产状集合AT′＝{α′_ij(θ，δ)}中，α′_ij表示地层界线sl_i的第j个边界点的最终产状。

进一步的，步骤(6)包括：

(6-1)读取SLine′中每个地层界线分段对应的边界点子集合SPoint′，依次获取其所有边界点，存入二维点集S2Point＝{s2p_u(x_u，y_u)|u＝1，2，...，n}，其中s2p_u代表第u个边界点，n代表边界点个数；对于每个二维点集S2Point，执行如下步骤(6-2)到(6-8)；

(6-2)获取二维点集S2Point中任一边界点s2p_u，及其相邻边界点s2p_u-1和s2p_u+1，由集合AT′获取s2p_u的倾向θ，若θ指向地层内部，则判断s2p_u属于异常波动部位点，转到步骤(6-3)，否则转到步骤(6-7)；

(6-3)若三点所成夹角开口朝内，将s2p_u的倾向翻转180°，并转到步骤(6-7)，否则转到步骤(6-4)；

(6-4)若s2p_u的倾向与s2p_u到其任一相邻边界点的向量的夹角小于三点所成夹角的补角，则转到步骤(6-5)，用曲率法计算s2p_u的倾向；

(6-5)计算边界点s2p_u与其相邻两个边界点的外接圆圆心坐标otp_u(x′_u,y′_u)；

(6-6)根据下式计算s2p_u处倾向；

(6-7)循环步骤(6-2)-(6-6)，直到二维点集S2Point被遍历完，识别并计算出所有异常波动部位边界点处的地层倾向，并更新到参考产状集合AT中；

(6-8)基于步骤(6-7)中新存入产状的边界点，重新将地层界线分段，更新集合SLine′，并转到步骤(5)，直到所有地层界线分段上都不存在异常产状。

步骤(7)包括：

(7-1)对于边界点sp_ij，过其坐标(x_ij,y_ij,z_ij)点，沿其倾向θ绘制单位长度的线，垂直于倾向θ绘制单位长度的线，完成产状符号的绘制；

(7-2)将倾向θ由极坐标系转为大地坐标系，对产状倾角δ进行标注。

本发明所述的基于地质图的地层界线逐点岩层产状获取装置，包括处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述方法。

有益效果：本发明与现有技术相比，其显著优点是：本发明充分运用了地质图产状信息，采用间接法计算地层整体产状趋势。该方法既避免了对钻孔、剖面等地质调查数据的过度依赖，又实现了地层界线上逐点的岩层产状获取。

附图说明

图1是本实施例中采用的地表DEM数据；

图2是实施例区域地质图数据；

图3是实施例产状点数据表；

图4是本发明方法总流程图；

图5是本发明方法基于四点发或三点法产状计算流程图；

图6是地层间实测产状赋值边界点示意图；

图7是曲率法计算产状示意图；

图8是实施例逐点产状示意图；

图9是实施例实测产状赋值结果图；

图10是实施例基于四点法产状计算结果图；

图11是实施例基于三点法产状计算结果图；

图12是实施例异常产状修正结果图；

图13是实施例边界点产状计算结果图；

具体实施方式

下面对本发明技术方案作进一步详细的说明，本实施例选取了湖北省秭归县发育在古生代地层中的黄陵穹窿。数据源为该区域30m分辨率的DEM数据(图1)和1:20万比例尺的地质图(图2)，地质图上包括两个地层，分别是陡山沱组和灯影组，该实验数据采用的投影坐标系为WGS_1984_UTM_Zone_49N。下面结合附图，并通过描述一个具体的实施例，来进一步说明。

本实施例提供了一种基于地质图的地层界线逐点岩层产状获取方法，如图4所示，具体包括如下步骤：

(1)读取岩层的地层面图层、DEM数据和已知产状点数据(图3)，得到地层集合Stra、地层界线集合SLine、产状点集合OPoint以及表示地层和产状点关系的地层产状点关系表SOR；该步骤具体包括：

(1-1)读取矢量格式的地层面图层，得到地层集合Stra＝{s_i|i＝1,2,3,…,sn}，其中s_i为第i个地层，sn为地层的数量；

(1-2)读取集合Stra中每个地层的边界线，并将其处理为单线，存入地层界线集合SLine＝{sl_i|i＝1，2，3,…,sn}，sl_i是s_i的地层界线；

(1-3)读取DEM数据，从中提取等高线存入等高线集合CourLine＝{cl_k|k＝1,2,3,…,cn},其中cl_k为第k条等高线素，cn为等高线数量；

(1-4)设定步长λ＝100m，对集合SLine中每个地层界线，按照设定步长获取其上所有边界点，并从DEM提取每一点所对应的高程值，存入到边界点集合SPoint＝{sp_ij(x_ij,y_ij,z_ij)|i＝1,2,3,…,sn，j＝1,2,3,…,pn}，其中sp_ij是地层界线sl_i的第j个边界点，x_ij,y_ij,z_ij分别为该边界点的横、纵坐标和高程值，pn为地层界线sl_i上边界点个数；之后在地层界线sl_i上曲率大于预设阈值的部位提取若干特征点，也作为边界点按顺序加入集合SPoint中；

(1-5)读取已知产状点数据，存入产状点集合OPoint＝{op_l|l＝1,2,...,on},其中op_l为第l个产状点，on代表产状点数量；

(1-6)根据包含关系，将每个产状点匹配到相应地层，存储到地层产状点关系表SOR＝{so_l(s_i,op_l)}中，so_l为s_i和op_l的关系。

(2)判断地层集合Stra中各相邻地层间相互关系，并基于地层间的相互关系，将实测产状信息赋值到相关边界点上；该步骤具体包括：

(2-1)从地层集合Stra中获取两个相邻地层s_i和s_i′，若两者皆是沉积岩，转到步骤(2-2)，否则判定相邻地层s_i和s_i′的关系为无显著关系，转到步骤(2-3)；

(2-2)判断s_i和s_i′相对年代属性之差，若为1，判定相邻地层s_i和s_i′的关系为两地层连续且为整合接触关系，若大于1，则判定相邻地层s_i和s_i′的关系为不整合接触关系；例如，实施例中陡山沱组和灯影组互为整合接触关系；

(2-3)将地层间关系存入地层间关系表SSR＝{ssr_i,i′}中，ssr_i,i′表示地层s_i和s_i′的关系；

(2-4)重复步骤(2-1)到(2-3)，直到完成所有相邻地层间关系的判断；

(2-5)从产状点集合OPoint中获取任一产状点op_l，经过op_l，沿op_l的倾向做直线el，从op_l出发，el两端遇到不整合接触的两个地层时被截断；

(2-6)el在op_l倾向方向与op_l所在地层s′_i有一交点，取地层s′_i的地层界线上距离该交点最近的一个边界点，将op_l的产状赋予该边界点，记录在参考产状集合AT中；

(2-7)获取el上一个不在s′_i上的交点，若交点在op_l倾向方向上，且同地层的另一个交点不在op_l倾向方向上，则转到步骤(2-8)；若交点和同地层的另一个交点都不在op_l倾向方向上，且该交点离op_l更近，则同样转到步骤(2-8)，否则转到(2-10)；

(2-8)若交点两侧地层为整合接触关系，则转到(2-9)；若为不整合接触关系，且交点两侧地层中，离op_l近的地层年代更大，则转到(2-9)，否则转到(2-10)；

(2-9)取地层s′_i的地层界线上距离该交点最近的一个边界点，将op_l的产状赋予该边界点，记录在参考产状集合AT中；

(2-10)重复步骤(2-5)-(2-9)，直到所有的产状点对边界点赋值完毕。赋值结果如图9所示。

(3)提取地层界线集合SLine中所有地层界线与等高线的交点，自适应地采用三点法或四点法计算地层界线的产状，并将产状赋值到对应的边界点上；例如，实施例中陡山沱组和灯影组的各点如图8所示，该步骤具体包括：

(3-1)从边界点集合SPoint中读取任一地层界线sl_i的边界点，获取所有等高线与sl_i边界点的交点存入集合CSP＝{csp_m(x_m,y_m,h_m)|m＝1,2,...,csn}，其中m代表交点序号，x_m,y_m,h_m分别为该交点的横、纵坐标和高程值，csn代表交点个数；

(3-2)从集合CSP中获取任意连续四个交点csp_m-1、csp_m、csp_m+1和csp_m+2，判别这四个交点是否满足相邻等高线法的下列选点规则；

a)交点csp_m-1和csp_m+2高程一致，交点csp_m和csp_m+1高程一致；

b)csp_m和csp_m+1位于同一条等高线；

c)向量

和

的最小夹角不大于阈值ε，即判定两向量方位近似平行；

如果符合，执行步骤(3-4)，如果不符合执行步骤(3-3)：

(3-3)判断交点csp_m-1、csp_m和csp_m+1是否满足三点法的下列选点规则；

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

b)交点csp_m-1、csp_m和csp_m+1不共线；

如果符合，执行步骤(3-5)；如果不符合，执行步骤(3-6)：

(3-4)基于交点csp_m-1、csp_m、csp_m+1和csp_m+2使用四点法计算边界点产状，如图5所示，执行步骤(3-6)；

(3-5)基于交点csp_m-1、csp_m、csp_m+1使用三点法计算边界点产状，执行步骤(3-6)；

(3-6)循环执行步骤(3-2)-(3-5)，直到集合CSP中所有连续交点被遍历完，得到sl_i所有边界点的产状；基于四点法、三点法计算产状结果如图10和图11所示；

(3-7)循环执行步骤(3-1)-(3-6)，直到地层界线集合SLine中所有地层界线被遍历完，得到所有地层界线的所有边界点的产状。

步骤(3-4)包括：

(3-4-1)对于交点csp_m-1、csp_m、csp_m+1和csp_m+2，计算平面向量

和

取

和

向量和的方位角度ρ为走向；

(3-4-2)根据下式计算倾向θ；

(3-4-3)计算交点csp_m-1和交点csp_m+2的中点，记为辅助点ap_m1(x′_m，y′_m，h′_m)，过点ap_m1沿θ方向作直线L，L与向量

交于点ap_m2(x″_m,y″_m,h″_m)，与地层界线sl_i交于点ap_m3(x″′_m,y″′_m,h″′_m)，根据下式计算倾角δ；

(3-4-4)将倾向θ、倾角δ作为ap_m3同一地层界线上最邻近的边界点的产状，存储到参考产状集合AT＝{α_ij(θ,δ)}中，α_ij表示地层界线sl_i的第j个边界点的产状。

步骤(3-5)包括：

(3-5-1)根据下式得到局部岩层层面方程Ax+By+Cz+D＝0的系数A、B和C，其中{A,B,C}即为三个边界点所确定三角面的法线，D为任意常数；

(3-5-2)根据系数A、B，得到岩层倾向线方程，其中M为任意实数；

Bx+Ay+M＝0

(3-5-3)采用下式计算倾向θ；

(3-5-4)采用下式计算倾角δ；

(3-5-5)将倾向θ、倾角δ作为csp_m同一地层界线上最邻近的边界点的产状存储到参考产状集合AT中。

(4)每个地层基于产状已知的边界点，将其地层界线分为若干分段，得到地层界线分段集合SLine′；该步骤具体包括：

(4-1)以产状已知的边界点为分割点，将地层界线集合SLine中所有地层界线划分为若干地层界线分段，存入地层界线分段集合SLine′＝{sl_ir|i＝1,2,3,…,sn，r＝1,2,3,…,ln_i′}，sl_ir表示sl_i第i个分段，sn为地层的数量，sl_r为地层分段要素，ln_i′为sl_i上分段数；

(4-2)按照地层界线分段，将地层边界点集合划分为若干个连续的边界点子集合SPoint′，每个子集合中，第一个边界点前和最后一个边界点后的两个边界点为产状已知的边界点。

(5)对于SLine′中每个地层界线分段，根据两端边界点线性插值计算分段内部边界点的产状，存入最终边界点产状集合AT′；该步骤具体包括：

(5-1)获取SLine′中任一地层界线分段对应的边界点子集合SPoint′，基于第一个边界点前和最后一个边界点后的两个边界点的产状，对每个边界点子集合SPoint′中的边界点产状线性插值；

(5-2)按照边界点和产状对应关系，将插值后的边界点产状存入最终边界点产状集合AT′＝{α′_ij(θ,δ)}中，α′_ij表示地层界线sl_i的第j个边界点的最终产状。

(6)根据边界点的最终产状对每个地层界线分段进行识别，若识别到地层界线分段上存在异常产状边界点，则利用曲率法近似计算其上边界点真实倾向，如图7所示，更新最终边界点产状集合AT′，如图12所示，并基于异常产状边界点划分当前地层界线分段，再转到步骤(5)，直到所有地层界线分段上都不存在异常产状边界点；该步骤具体包括：

(6-1)读取SLine′中每个地层界线分段对应的边界点子集合SPoint′，依次获取其所有边界点，存入二维点集S2Point＝{s2p_u(x_u，y_u)|u＝1，2，...，n}，其中s2p_u代表第u个边界点，n代表边界点个数；对于每个二维点集S2Point，执行如下步骤(6-2)到(6-8)；

(6-2)获取二维点集S2Point中任一边界点s2p_u，及其相邻边界点s2p_u-1和s2p_u+1，由集合AT′获取s2p_u的倾向θ，若θ指向地层内部，则判断s2p_u属于异常波动部位点，转到步骤(6-3)，否则转到步骤(6-7)；

(6-3)若三点所成夹角开口朝内，将s2p_u的倾向翻转180°，并转到步骤(6-7)，否则转到步骤(6-4)；

(6-4)若s2p_u的倾向与s2p_u到其任一相邻边界点的向量的夹角小于三点所成夹角的补角，则转到步骤(6-5)，用曲率法计算s2p_u的倾向；

(6-5)计算边界点s2p_u与其相邻两个边界点的外接圆圆心坐标otp_u(x′_u,y′_u)；

(6-6)根据下式计算s2p_u处倾向；

(6-7)循环步骤(6-2)-(6-6)，直到二维点集S2Point被遍历完，识别并计算出所有异常波动部位边界点处的地层倾向，并更新到参考产状集合AT中；

(6-8)基于步骤(6-7)中新存入产状的边界点，重新将地层界线分段，更新集合SLine′，并转到步骤(5)，直到所有地层界线分段上都不存在异常产状。

(7)根据最终边界点产状集合AT′生成地质图地层边界点产状数据，如图13所示；该步骤具体包括：

(7-1)对于边界点sp_ij，过其坐标(x_ij,y_ij,z_ij)点，沿其倾向θ绘制单位长度的线，垂直于倾向θ绘制单位长度的线，完成产状符号的绘制；

(7-2)将倾向θ由极坐标系转为大地坐标系，对产状倾角δ进行标注。

本实施例中仅基于GDAL开源代码提供的影像处理接口读取DEM数据和导入导出矢量数据，该方法也可以使用SuperMap、QGIS等GIS软件的接口。

本实施例还提供一种基于地质图的地层界线逐点岩层产状获取装置，包括处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述方法。

Claims (8)

1.一种基于地质图的地层界线逐点岩层产状获取方法，其特征在于该方法具体包括：

(1)读取岩层的地层面图层、DEM数据和已知产状点数据，得到地层集合Stra、地层界线集合SLine、产状点集合OPoint以及表示地层和产状点关系的地层产状点关系表SOR；

具体包括：

(1-1)读取矢量格式的地层面图层，得到地层集合Stra＝{s_i|i＝1,2,3,…,sn}，其中s_i为第i个地层，sn为地层的数量；

(1-2)读取集合Stra中每个地层的边界线，并将其处理为单线，存入地层界线集合SLine＝{sl_i|i＝1,2,3,…,sn}，sl_i是s_i的地层界线；

(1-3)读取DEM数据，从中提取等高线存入等高线集合CoutLine＝{cl_k|k＝1,2,3,…,cn},其中cl_k为第k条等高线素，cn为等高线数量；

(1-4)设定步长λ，对集合SLine中每个地层界线，获取其上所有边界点，并从DEM提取每一点所对应的高程值，存入到边界点集合SPoint＝{sp_ij(x_ij,y_ij,z_ij)|i＝1,2,3,…,sn，j＝1,2,3,…,pn}，其中sp_ij是地层界线sl_i的第j个边界点，x_ij,y_ij,z_ij分别为该边界点的横、纵坐标和高程值，pn为地层界线sl_i上边界点个数；之后在地层界线sl_i上曲率大于预设阈值的部位提取若干特征点，也作为边界点按顺序加入集合SPoint中；

(1-5)读取已知产状点数据，存入产状点集合OPoint＝{op_l|l＝1,2,...,on},其中op_l为第l个产状点，on代表产状点数量；

(1-6)根据包含关系，将每个产状点匹配到相应地层，存储到地层产状点关系表SOR＝{so_l(s_i,op_l)}中，so_l为s_i和op_l的关系；

(2)判断地层集合Stra中各相邻地层间相互关系，并基于地层间的相互关系，将实测产状信息赋值到相关边界点上；

具体包括：

(2-1)从地层集合Stra中获取两个相邻地层s_i和s_i′，若两者皆是沉积岩，转到步骤(2-2)，否则判定相邻地层s_i和s_i′的关系为无显著关系，转到步骤(2-3)；

(2-2)判断s_i和s_i′相对年代属性之差，若为1，判定相邻地层s_i和s_i′的关系为两地层连续且为整合接触关系，若大于1，则判定相邻地层s_i和s_i′的关系为不整合接触关系；

(2-3)将地层间关系存入地层间关系表SSR＝{ssr_i,i′}中，ssr_i,i′表示地层s_i和s_i′的关系；

(2-4)重复步骤(2-1)到(2-3)，直到完成所有相邻地层间关系的判断；

(2-5)从产状点集合OPoint中获取任一产状点op_l，经过op_l，沿op_l的倾向做直线el，从op_l出发，el两端遇到不整合接触的两个地层时被截断；

(2-6)el在op_l倾向方向与op_l所在地层s′_i有一交点，取地层s′_i的地层界线上距离该交点最近的一个边界点，将op_l的产状赋予该边界点，记录在参考产状集合AT中；

(2-7)获取el上一个不在s′_i上的交点，若交点在op_l倾向方向上，且同地层的另一个交点不在op_l倾向方向上，则转到步骤(2-8)；若交点和同地层的另一个交点都不在op_l倾向方向上，且该交点离op_l更近，则同样转到步骤(2-8)，否则转到(2-10)；

(2-8)若交点两侧地层为整合接触关系，则转到(2-9)；若为不整合接触关系，且交点两侧地层中，离op_l近的地层年代更大，则转到(2-9)，否则转到(2-10)；

(2-9)取地层s′_i的地层界线上距离该交点最近的一个边界点，将op_l的产状赋予该边界点，记录在参考产状集合AT中；

(2-10)重复步骤(2-5)-(2-9)，直到所有的产状点对边界点赋值完毕；

(3)提取地层界线集合SLine中所有地层界线与等高线的交点，自适应地采用三点法或四点法计算地层界线的产状，并将产状赋值到对应的边界点上；

(4)每个地层基于产状已知的边界点，将其地层界线分为若干分段，得到地层界线分段集合SLine′；

(5)对于SLine′中每个地层界线分段，根据两端边界点线性插值计算分段内部边界点的产状，存入最终边界点产状集合AT′；

(6)根据边界点的最终产状对每个地层界线分段进行识别，若识别到地层界线分段上存在异常产状边界点，则利用曲率法近似计算其上边界点真实倾向，更新最终边界点产状集合AT′，并基于异常产状边界点划分当前地层界线分段，再转到步骤(5)，直到所有地层界线分段上都不存在异常产状边界点；

(7)根据最终边界点产状集合AT′生成地质图地层边界点产状数据。

2.根据权利要求1所述的基于地质图的地层界线逐点岩层产状获取方法，其特征在于：步骤(3)包括：

(3-1)从边界点集合SPoint中读取任一地层界线sl_i的边界点，获取所有等高线与sl_i边界点的交点存入集合CSP＝{csp_m(x_m,y_m,h_m)|m＝1,2,...,csn}，其中m代表交点序号，x_m,y_m,h_m分别为该交点csp_m的横、纵坐标和高程值，csn代表交点个数；

(3-2)从集合CSP中获取任意连续四个交点csp_m-1、csp_m、csp_m+1和csp_m+2，判别这四个交点是否满足相邻等高线法的下列选点规则；

a)交点csp_m-1和csp_m+2高程一致，交点csp_m和csp_m+1高程一致；

b)csp_m和csp_m+1位于同一条等高线；

c)向量

和

的最小夹角不大于阈值ε，即判定两向量方位近似平行；

如果符合，执行步骤(3-4)，如果不符合执行步骤(3-3)：

(3-3)判断交点csp_m-1、csp_m和csp_m+1是否满足三点法的下列选点规则；

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

b)交点csp_m-1、csp_m和csp_m+1不共线；

如果符合，执行步骤(3-5)；如果不符合，执行步骤(3-6)：

(3-4)基于交点csp_m-1、csp_m、csp_m+1和csp_m+2使用四点法计算边界点产状，执行步骤(3-6)；

(3-5)基于交点csp_m-1、csp_m、csp_m+1使用三点法计算边界点产状，执行步骤(3-6)；

(3-6)循环执行步骤(3-2)-(3-5)，直到集合CSP中所有连续交点被遍历完，得到sl_i所有边界点的产状；

(3-7)循环执行步骤(3-1)-(3-6)，直到地层界线集合SLine中所有地层界线被遍历完，得到所有地层界线的所有边界点的产状。

3.根据权利要求2所述的基于地质图的地层界线逐点岩层产状获取方法，其特征在于：步骤(3-4)包括：

(3-4-1)对于交点csp_m-1、csp_m、csp_m+1和csp_m+2，计算平面向量

和

取

和

向量和的方位角度ρ为走向；

(3-4-2)根据下式计算倾向θ；

(3-4-3)计算交点csp_m-1和交点csp_m+2的中点，记为辅助点ap_m1(x′_m,y′_m,h′_m)，过点ap_m1沿θ方向作直线L，L与向量

交于点ap_m2(x″_m,y″_m,h″_m)，与地层界线sl_i交于点ap_m3(x″′_m,y″′_m,h″′_m)，根据下式计算倾角δ；

(3-4-4)将倾向θ、倾角δ作为ap_m3同一地层界线上最邻近的边界点的产状，存储到参考产状集合AT＝{α_ij(θ,δ)}中，α_ij表示地层界线sl_i的第j个边界点的产状。

4.根据权利要求2所述的基于地质图的地层界线逐点岩层产状获取方法，其特征在于：步骤(3-5)包括：

(3-5-1)根据下式得到局部岩层层面方程Ax+By+Cz+D＝0的系数A、B和C，其中{A,B,C}即为三个边界点所确定三角面的法线，D为任意常数；

(3-5-2)根据系数A、B，得到岩层倾向线方程，其中M为任意实数；

Bx+Ay+M＝0

(3-5-3)采用下式计算倾向θ；

(3-5-4)采用下式计算倾角δ；

(3-5-5)将倾向θ、倾角δ作为csp_m同一地层界线上最邻近的边界点的产状存储到参考产状集合AT中。

5.根据权利要求1所述的基于地质图的地层界线逐点岩层产状获取方法，其特征在于：步骤(4)包括：

(4-1)以产状已知的边界点为分割点，将地层界线集合SLine中所有地层界线划分为若干地层界线分段，存入地层界线分段集合SLine′＝{sl_ir|i＝1,2,3,…,sn，r＝1,2,3,…,ln_i′}，sl_ir表示sl_i第i个分段，sn为地层的数量，sl_r为地层分段要素，ln_i′为sl_i上分段数；

(4-2)按照地层界线分段，将地层边界点集合划分为若干个连续的边界点子集合SPoint′，每个子集合中，第一个边界点前和最后一个边界点后的两个边界点为产状已知的边界点。

6.根据权利要求1所述的基于地质图的地层界线逐点岩层产状获取方法，其特征在于：步骤(5)包括：

(5-1)获取SLine′中任一地层界线分段对应的边界点子集合SPoint′，基于第一个边界点前和最后一个边界点后的两个边界点的产状，对每个边界点子集合SPoint′中的边界点产状线性插值；

(5-2)按照边界点和产状对应关系，将插值后的边界点产状存入最终边界点产状集合AT′＝{α′_ij(θ,δ)}中，α′_ij表示地层界线sl_i的第j个边界点的最终产状。

7.根据权利要求1所述的基于地质图的地层界线逐点岩层产状获取方法，其特征在于：步骤(6)包括：

(6-1)读取SLine′中每个地层界线分段对应的边界点子集合SPoint′，依次获取其所有边界点，存入二维点集S2Point＝{s2p_u(x_u,y_u)|u＝1,2,...,n}，其中s2p_u代表第u个边界点,n代表边界点个数；对于每个二维点集S2Point，执行如下步骤(6-2)到(6-8)；

(6-2)获取二维点集S2Pount中任一边界点s2p_u，及其相邻边界点s2p_u-1和s2p_u+1，由集合AT′获取s2p_u的倾向θ，若θ指向地层内部，则判断s2p_u属于异常波动部位点，转到步骤(6-3)，否则转到步骤(6-7)；

(6-3)若三点所成夹角开口朝内，将s2p_u的倾向翻转180°，并转到步骤(6-7)，否则转到步骤(6-4)；

(6-4)若s2p_u的倾向与s2p_u到其任一相邻边界点的向量的夹角小于三点所成夹角的补角，则转到步骤(6-5)，用曲率法计算s2p_u的倾向；

(6-5)计算边界点s2p_u与其相邻两个边界点的外接圆圆心坐标otp_u(x′_u,y′_u)；

(6-6)根据下式计算s2p_u处倾向；

(6-7)循环步骤(6-2)-(6-6)，直到二维点集S2Point被遍历完，识别并计算出所有异常波动部位边界点处的地层倾向，并更新到参考产状集合AT中；

(6-8)基于步骤(6-7)中新存入产状的边界点，重新将地层界线分段，更新集合SLine′，并转到步骤(5)，直到所有地层界线分段上都不存在异常产状。

8.一种基于地质图的地层界线逐点岩层产状获取装置，包括处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于：所述处理器执行所述程序时实现权利要求1-7中任意一项所述的基于地质图的地层界线逐点岩层产状获取方法。

Images