CN105956066B - 一种褶皱地貌类型的自动化识别方法 - Google Patents

Abstract

本发明公开了一种褶皱地貌类型的自动化识别方法，具体包括如下步骤：1)根据核部地层与横切剖面线提取原则，提取可能存在褶皱的场景条带；2)基于邻接属性关系图建模方法，对场景条带进行场景空间结构的建模与化简；3)通过判别化简后的邻接属性关系模型能否由对应的褶皱规则推导而来，实现针对不同褶皱构造模式的匹配识别；4)基于识别出的褶皱场景进行山和谷状态的进一步判别，最终确定褶皱构造地貌类型，粗略划定褶皱地貌单元。该方法针对褶皱“两翼对称”的出露规律，以结构模式匹配的方法，基本实现了在基岩山地地区，对褶皱构造地貌的自动化识别。

Description

一种褶皱地貌类型的自动化识别方法

技术领域

本发明属于地理信息技术应用领域，具体涉及一种基于地质图与DEM，通过空间结构模式匹配，实现自动化识别褶皱构造地貌单元的方法。

背景技术

褶皱几乎控制了大中型地貌的基本形态，在一定程度上反应了其力学背景，对于矿产的形成、分布等也有一定的控制作用，在地质测量、地貌形成和地貌分类等研究中具有重要意义。褶皱地貌的识别是有效支持地质图认知、构造解析、三维建模等工作的基础。

传统的地貌类型划分方法主要为人工操作，其工作量大，历时长，且受限于制图人员的技术和认识。因此，众多研究学者转而探讨自动化的地貌类型识别与地貌单元提取方法。当前常用的自动识别方法主要基于DEM或遥感影像，通过单元统计、特征提取或标识解译等方法，实现针对地貌形态或典型地貌的自动化识别。而褶皱构造地貌因其特有的“对称重复”地层组合空间结构，目前常用的基于地形信息进行单元统计或是特征提取的自动化地貌类型识别方法并不适用。

褶皱地貌类型的自动化识别与单元提取，应首先识别出褶皱构造，再对构造内的地貌进行综合判别。虽然，褶皱构造形态复杂、类型多样，且没有具体的边界；但其本质为一组受力弯曲的地层，在剥蚀后呈现出两翼对称的、模式化的特定地层组合结构。因此，本发明采用地形地质图、DEM等数据源，基于空间结构场景建模和文法模式匹配等方法，实现褶皱地貌的自动识别。

发明内容

本发明的主要目的在于：针对褶皱地貌“对称重复”的地层组合空间结构特征，基于地形与地质信息，提出一种构造地貌自动化识别方法。

本发明采用的技术方案如下：

一种褶皱地貌类型的自动化识别方法，包括如下步骤：

1)根据核部地层与横切剖面线提取原则，提取可能存在褶皱的场景条带；

2)基于邻接属性关系图建模方法，对场景条带进行场景空间结构的建模与化简；

3)通过判别化简后的邻接属性关系模型能否由对应的褶皱规则推导而来，实现针对不同褶皱构造模式的匹配识别；

4)基于识别出的褶皱场景进行山和谷形态的进一步判别，最终确定褶皱构造地貌类型，粗略划定褶皱地貌单元。

本发明的方法针对褶皱“两翼对称”的出露规律，以结构模式匹配的方法，基本实现了在基岩山地地区，对褶皱构造地貌的自动化识别。相比传统的人工识别方法，提高了效率，也为非专业人员的地质图读图识图提供了帮助。此外，相比目前主要基于地表形态进行的地貌自动识别提取方法，结合考虑了地层结构分布的信息，解决了空间结构特征明显的地貌自动识别的问题。同时，也可作为地表地貌自动识别方法的补充，为构造地貌的自动化识别提供了思路。

附图说明

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

图2示例地层年代序列表截图；

图3实施例实验数据，(a)地质体面图层数据、(b)DEM数据；

图4核部地层选取示例；

图5实施例识别场景条带结果图；

图6邻接对象判别示例；

图7实施例场景模型，图中的“·”代表地层对象点，“—”代表对象邻接关系；

图8实施例化简后场景模型，图中的“·”代表化简后地层对象点，“—”代表化简后关系；

图9实施例褶皱地貌单元。

具体实施方式

本发明流程图如图1所示，褶皱地貌类型自动识别的具体方法步骤如下：

(一)数据预处理

步骤11：本文方法是在邻接属性关系图的基础上进行的，地层对象与地层间关系详细属性列表如下：

步骤12：加载shp格式的地质面图层数据，得到所有地层面要素集合OriS＝{os_i|i＝0，1，2，...，n-1}，n为地层面要素的个数。要求os_i属性表中至少包括一个[Age]字段；

步骤13：对os_i属性赋值。os_i[Id]与os_i[SourceId]赋值为i；若os_i[Trend]不为空，则根据公式(1)将其转换为极坐标系角度；若os_i[Trend]为空，则通过计算os_i多边形的走向(参考文献：《一个求解多边形最小面积外接矩形的算法》，程鹏飞、闫浩文、韩振辉，工程图学学报，2008，(1))来代替；

os_i[Length]赋值os_i多边形最小外包矩形的长边长度值；os_i[Gravity]赋值其重心点坐标gp_i(gx_i，gy_i)，通过公式(2)，基于地层面要素os_i的顶点集OriV_i＝{ov_ij(x_ij，y_ij)|j＝0，1，2，...，m_i-1}(m_i为面要素os_i的顶点个数)来计算得到；

步骤14：配置XML格式的“地层年代序列表”，示例如图2。根据os_i的[Age]属性，在“地层年代序列表”中查找到节点node_j，使node_j[name]＝os_i[Age]：

1)若node_j[RockType]为“I”，则os_i[IsIntrusion]赋值“True”，os_i[StateTag]赋值“Delete”；

2)若node_j[RockType]为“S”，则os_i[IsIntrusion]赋值“False”。再判断node_i是否是“Q”的子节点，若是，则os_i[StateTag]赋值“Delete”；反之，os_i[StateTag]赋值“Keep”；

(二)识别条带的提取

步骤21：获取集合OriS＝{os_i|i＝0，1，2，...，n-1}中每个os_i面要素的邻接地层面要素集合NeiS_i，并从NeiS_i中移除属性[StateTag]为“Delete”的元素；

步骤22：对于属性[StateTag]为“Keep”的地层面要素os_i，判断其邻接要素集合NeiS_i＝{ns_ij|j＝0，1，2，...，l_i-1}(l_i为os_i的邻接地层中属性[StateTag]为“Keep”的地层个数)是否满足以下条件：

若满足以上任一条件，则判定面要素os_i为核部地层要素，将其加入核部地层要素集合CorS＝{cs_j|j＝0，1，2，...，p-1}，p为查找到的核部地层的数量；

步骤23：针对每个核部地层cs_j，过其重心点核部地层cs_j[Gravity]坐标垂直于走向cs_j[Trend]做垂线sl_j，得到剖面线集合SecL＝{sl_j|j＝0，1，2，...，p-1}，对sl_j逐个做缓冲区，得到裁剪面要素集合CutP＝{cp_j|j＝0，1，2，...，p-1}。给裁剪面cp_j[SourceCS]属性赋值cs_j[Id]，裁剪面cp_j[CSTrend]属性赋值cs_j[Trend]；

步骤24：用裁剪面要素集合CutP对原始地层面要素集合OriS进行裁剪，得到场景地层面要素集合SceS＝{ss_k|k＝0，1，2，...，s-1}，s为识别条带上地层对象的数量；

步骤25：对面要素ss_k的属性进行赋值。对于裁剪面cp_j与原始地层面os_i裁剪而成的面要素ss_k而言，ss_k[SourceCS]和ss_k[CSTrend]继承cp_j的同名属性；ss_k[SourceID]、ss_k[Trend]、ss_k[Length]、ss_k[IsIntrusion]和ss_k[StateTag]分别继承os_i的同名属性；ss_k[Id]赋值k；ss_k[Area]和ss_k[Gravity]基于ss_k的多边形要素重新计算；

(三)识别结构建模与化简

步骤31：由上一步获得地层面要素集合SceS，根据ss_k的[Gravity]属性绘制点图层SceneObjectLayer，得到地层对象点集SceG＝{sg_i|i＝0，1，2，...，s-1}，对象点sg_i继承面要素ss_k的所有属性；

步骤32：对于任意面要素ss_i与ss_j(j>i)，判断是否满足以下条件：

若满足，则操作步骤33：

步骤33：连接地层对象点sg_i与sg_i绘制线rl_k，对线rl_k属性赋值如下：

查找地层对象点sg_i与sg_j的属性[Age]在“地层年代序列表”中对应节点的编号node_a[Id]与node_b[Id]：

1)若node_a[Id]≤node_b[Id]，rl_k[StartId]赋值为sg_i[Id]；rl_k[EndId]赋值为sg_j[Id]；rl_k[TimeDis]赋值为node_b[Id]-node_a[Id]；rl_k[AzimuthAng]赋值由对象点sg_i的坐标指向对象点sg_j的坐标的方位角ang_ab；

2)若node_a[Id]＞node_b[Id]，rl_k[StartId]赋值sg_j[Id]；rl_k[EndId]赋值sg_i[Id]；rl_k[TimeDis]赋值node_a[Id]-node_b[Id]；rl_k[AzimuthAng]赋值由对象点sg_j的坐标指向对象点sg_i的坐标的方位角ang_ba；

步骤34：得到关系线集合RelL＝{rl_k|k＝0，1，2，...，t-1}(t为条带上的邻接关系数量)后，根据关系线集合RelL建立无向邻接集合SimAList＝{sal_i|i＝0，1，2，...，s-1}。其中，集合sal_i代表属性[Id]等于“i”的点所邻接对象的属性[Id]集合。每个关系rl_k，根据属性rl_k[StartId]与rl_k[EndId]，往i＝rl_k[StartId]的sal_i中添加rl_k[EndId]，i＝rl_k[EndId]的sal_i中添加rl_k[StartId]；

步骤35：除了属性[Source]与[SourceC]相等的核部地层对象点不进行化简操作外，对每个地层对象点sg_i及其邻接面要素化简情况和处理方案如下：

1)若sg_i[StateTag]为“Keep”，且集合sal_i长度等于1。对于唯一的对象点sg_j[Id]∈sal_i，若满足sg_j[StateTag]为“Delete”，则sg_i[StateTag]赋值“Delete”，清空集合sal_i使其长度为0；

2)若sg_i[StateTag]为“Keep”，且集合sal_i长度大于0。若存在对象点sg_j，满足①sg_i[StateTag]为“Keep”；②sg_i[Id]∈sal_i；③sg_j[Age]与sg_i[Age]在“地层年代表序列表”中对应的编号node_a[Id]与node_b[Id]相等。则对对象点sg_i与sg_j进行如下的合并操作：

a)令a_i＝sg_i[Area]，a_j＝sg_j[Area]，根据公式(3)重新计算点坐标(x_i，y_i)，并赋值给sg_i[GravityPoint]。保留对象点sg_i的其他属性；

b)对象点sg_j[StateTag]赋值为“Delete”；

c)将集合sal_j整合到集合sal_i中，清空集合sal_j使其长度为0；

3)若sg_i[StateTag]为“Delete”，且集合sal_i长度等于1。对于唯一的sg_u[Id]∈sal_i。若满足sg_u[StateTag]为“Delete”，执行步骤35.2).a)～c)对对象点sg_i与sg_u进行合并操作；

4)若sg_i[StateTag]为“Delete”，且集合sal_i长度等于2。存在对象点sg_u和sg_v(v>u)，sg_u[Id]∈sal_i，sg_v[Id]∈sal_i。若满足则将sg_v[Id]加入集合sal_u，将sg_u[Id]加入集合sal_v。清空集合sal_i使其长度为0；

5)若sg_i[StateTag]为“Delete”，且集合sal_i长度大于2。至少存在对象sg_u点和sg_v(v>u)，满足sg_u[Id]∈sal_i，sg_v[Id]∈sal_i。对对象点sg_u和sg_v进行如下判断：

c)若存在对象点sg_u与sg_v满足：①sg_u[Age]与sg_v[Age]在“地层年代表序列表”中对应的编号node_a[Id]与node_b[Id]相等，②sg_u[StateTag]与sg_v[StateTag]都为“Keep”，③方位角ang_uv或ang_vu与sg_i[CSTrend]差值小于容限。则执行步骤35.2).a)～c)对对象点sg_u与sg_v进行合并操作；

d)若存在对象点sg_u与sg_v满足：①方位角ang_ui与方位角ang_iv差值不超出容限，②对于任意sg_t[Id]∈sal_u，方位角ang_ut与ang_tv的差值均超出容限，③对于任意sg_t[Id]∈sal_v，方位角ang_vt与ang_tu的差值均超出容限，④则将sg_v[Id]从集合sal_i中移除并加入集合sal_u，将sg_u[Id]从集合sal_i中移除并加入集合sal_v；

e)不再出现以上情况，则清空集合sal_i使其长度为0；

步骤36：遍历邻接关系集合SimAList，若对象点sg_i[StateTag]为“Delete”，清空集合sal_i使其长度为0；反之，若存在对象点sg_j满足sg_j[StateTag]为“Delete”，且sg_j[Id]∈sal_i，则将sg_j[Id]从集合sal_i中移除；

步骤37：完成上一步后，遍历sg_i[StateTag]为“Keep”的对象，若其集合sal_i的长度大于2，则分别计算sg_jsg_i与其邻接对象的方位角，保留与sg_i[CSTrend]最接近垂直的两个角度及其对象，将其它对象的[Id]从集合sal_i中移除；

步骤38：取点集SceG中属性[StateTag]为“Keep”的点绘制化简后的地层对象点图层SimSOLayer。得到化简后地层对象点集SimSG＝{ssg_j|j＝0，1，2，...，c-1}，c为化简后保留的地层对象的数量，即SceG中属性[StateTag]为“Keep”的对象的数量；

步骤39：根据删改后的邻接关系集合SimAList＝{sal_i|i＝0，1，2，...，s-1}绘制化简后的有向线图层SimRLLayer。对于每个对象点ssg_j，若存在ssg_u[Id]∈sal_i，则按照步骤33连接地层对象点ssg_u与ssg_j，绘制线srl_k。得到化简后的关系线集合SimRL＝{srl_k|k＝0，1，2，...，r-1}，r为化简后的邻接关系数量。

(四)褶皱构造模式匹配

步骤41：创建集合OriFold＝{of_i|i＝0，1，2，...，p-1}存放褶皱对象，初始化的褶皱对象数量与核部地层对象点的数量相等。褶皱对象的数据属性如表：

步骤42：基于化简后的地层对象点集SimSG＝{ssg_j|j＝0，1，2，...，c-1}与关系线集SimRL＝{srl_k|k＝0，1，2，...，r-1}，进行褶皱构造识别。获取ssg_j[SourceCS]＝ssg_j[SourceId]的对象点加入点集CorG，得到化简后的核部地层对象点集合CorG＝{cg_i|i＝0，1，2，...，p-1}。对每个褶皱对象of_i属性进行赋值，of_i[SourceCS]赋值cg_i[SourceCS]；of_i[IdList]中添加cg_i[Id]；

步骤43：取cg_i[Id]的值给索引编号“LeftS”与“RightS”，对化简后关系线集合SimRL中的元素进行如下规则判别，其中预先设定tolerance为容限：

1)若存在对象满足背斜规则集，则属性of_i[FoldType]赋值“Anticline”。将srl_a[EndId]插入of_i[IdList]的第一个位置，将srl_b[EndId]添加到of_i[IdList]的最后，将srl_a[EndId]与srl_b[EndId]赋值给索引编号“LeftS”与“RightS”，重复本规则集的判别直至不再满足背斜规则集，则结束判别。背斜规则集如下：

2)若存在对象满足向斜规则集，则属性of_i[FoldType]赋值“Syncline”。将srl_a[StartId]插入of_i[IdList]的第一个位置，将srl_b[StartId]添加到of_i[IdList]的最后，将srl_a[StartId]与srl_b[StartId]赋值给索引编号“LeftS”与“RightS”，重复本规则集的判别直至不再满足向斜规则集，则结束判别。向斜规则集如下：

步骤44：对所有核部地层对象点cg_i完成规则匹配后，可知初始褶皱对象集合OriFold＝{of_i|i＝0，1，2，...，p-1}中包括了不符合褶皱规则的褶皱对象，遍历OriFold，若of_i[IdList]的长度小于3，则将of_i从中OriFold移除，最终得到褶皱单元集合Fold＝{f_i|i＝0，1，2，...，q-1}，q为符合褶皱规则的褶皱单元数量；

步骤45：对于褶皱单元f_i，依次找到属性[Id]值等于f_i[IdList][0]至f_i[IdList][ListCount-1]的化简前地层对象点sg_j，并将sg_j[SourceId]添加到f_i[SIdList]中；

步骤46：对集合Fold中两两褶皱的重复部分进行舍取。对于褶皱单元f_i，遍历褶皱单元f_j(j＝0，1，2，...，q-1,j≠i)：

1)若存在某个褶皱单元f_j(j>i)，满足f_i[SIdList]与f_j[SIdList]完全相同，或是f_i[SIdList]与f_j[SIdList].Reverse()完全相同，则从褶皱单元集合Fold中移除褶皱单元f_j；

2)比较f_i[SIdList]数组与每个f_j[SIdList]数组中的元素，若存在相等的元素，记录其在f_i[SIdList]数组中的下标。对于所有记录的下标，取最接近ListCount/2且小于ListCount/2的下标a和最接近ListCount/2且大于ListCount/2的下标b，保留f_i[IdList]数组f_i[IdList][a]至f_i[IdList][b]的部分，保留f_i[SIdList]数组f_i[SIdList][a]至f_i[SIdList][b]的部分；

步骤46：对褶皱单元集合Fold＝{f_i|i＝1，2，3，...，q}中的褶皱单元f_i其余属性进行赋值。获取属性[SourceId]与[SourceCS]相等，且等于f_i[SourceCS]属性的地层对象点sg_j，f_i[Trend]赋值sg_j[Trend]，f_i[Gravity]赋值sg_j[Gravity]，f_i[a]赋值sg_j[Length]；获取属性[Id]等于f_i[IdList][0]和f_i[IdList][ListCount-1]的地层对象点sg_left与sg_right，计算sg_left[Gravity]与sg_left[Gravity]间的距离d，f_i[b]赋值d。

(五)褶皱地貌类型确定

步骤51：加载同一投影坐标系下的DEM数据；

步骤52：根据规则识别得到的褶皱单元集合Fold＝{f_i|i＝0，1，２，...，q-1}，对于褶皱单元f_i，找到属性[Id]为f_i[IdList][0]、f_i[IdList][ListCount/2]和f_i[IdList][ListCount-1]的地层对象点sg_left、sg_center与sg_right，根据地层对象点的[Gravity]属性，从栅格数据上获取高程分别为L、C和R；

步骤53：若L-C＞0&R-C＞0，则f_i[Landform]赋值“Valley”；反之，f_i[Landform]赋值“Mountain”；

步骤54：以f_i[Gravity]为中心点，以f_i[a]为长轴，以f_i[b]为短轴，以f_i[Trend]为旋转角度，绘制椭圆表示褶皱地貌单元，得到褶皱地貌类型图层FoldLandformLayer。

实施例

下面结合附图并通过描述一个基于地形地质信息自动识别褶皱地貌类型并提取地貌单元的实例，来进一步说明本发明的效果。本实例选择1：50000的庐山地区的地质体面图层和DEM图层，如图3所示。

具体实施过程如下：

(一)数据预处理

步骤11：加载shp格式的地质体面图层数据“地质分区”，读取地层面要素读入集合OriS＝{os_i|i＝1，2，...，137}，本实施例中共包含138个地层面要素；

步骤12：加载的“地质分区”图层基本参数为：宽13034m，高9930m，图幅区域为点集{(399008.2,3271576.7),(412171.3,3271576.7),(412171.3,3281630.9),(399008.2,3281630.9)}围成的矩形。识别条带宽带可以默认取值为图幅长度的1/100，也可以人工设置，本实施例中设置取值为200m。设置化简和褶皱匹配时的角度容限tolerance为15°；

步骤12：“地质分区”图层的属性表中包括[Age]和[Trend]两个字段，以面要素os₂₇和面要素os₁₀₄为例，os₂₇自带属性包括os₂₇[Age]为“Z_1l^3”，os₂₇[Trend]为“51”，os₁₀₄自带属性包括os₁₀₄[Age]为“Z_1l^1”，os₁₀₄[Trend]为“47”。通过计算，对其部分属性进行赋值如下：

步骤13：在加载的“地层年代序列表”中os₂₇[Age]对应条目node₃，os₁₀₄[Age]对应条目node₁，对os_i部分属性赋值如下：

(二)识别条带的提取

步骤21：获取OriS中每个面要素os_i的邻接地层要素集合NeiS_i，以要素os₂₇和要素os₁₀₅为例，要素os₂₇的邻接地层要素集合NeiS₂₇有4个要素，面要素os₁₀₄的邻接地层要素集合NeiS₁₀₄有4个要素，如图4所示。各个邻接地层要素的属性以及对应在“地层年代序列表”的条目的分别如下：

则依据[StateTag]属性，将面要素os₁₁₆和os₁₃₄从NeiS₂₇中移除，将面要素os₁₁₆从NeiS₁₀₄中移除；

步骤22：根据规则选取核部地层，以面要素os₂₇和os₁₀₄为例，规则判别如下：

1)因os₂₇[StateTag]为“Keep”，对集合NeiS₂₇进行判断：面要素os₂₇对应的node[Id]为3，NeiS₂₇中的面要素os₁₀₅其node[Id]为2小于os₂₇，而面要素os₂₃其node[Id]为6大于os₂₇，不满足规则；此外，NeiS₂₇长度为2，不满足规则。因此，要素os₂₇不是核部地层；

2)因os₁₀₄[StateTag]为“Keep”，对集合NeiS₁₀₄进行判断：面要素os₁₀₄对应的node[Id]为1，NeiS₁₀₄中的面要素os₁₀₅、os₁₂₉和os₁₃₀对应的node[Id]分别为2、1和2，都大于或完全等于os₁₀₄，满足规则。因此，要素os₁₀₄满足核部地层选取原则，加入核部地层面要素集合CorS；

对所有属性[StateTag]为“Keep”的地层os_i进行规则判别，得到核部地层要素集合CorS＝{cs_l|l＝0，1，2，...，40}，其中cs₃₅即为os₁₀₄；

步骤23：按照属性对核部地层要素集合CorS中的所有核部地层面要素绘制剖面线，得到其裁剪面要素集合CutP＝{cp_l|l＝0，1，2，...，40}。以核部地层cs₃₅为例，在cs₃₅的重心点(402597.7,3277015.0)处，沿着133°方向绘制直线，与图幅边界相交，得到直线sl₃₅<(399008.2,3280248.602),(408616.3,3271576.7)>。沿直线以200m为宽，做缓冲区，得到裁剪面cp₃₅，根据核部地层要素cs₃₅的属性[Id]和[Trend]给裁剪面cp₃₅赋值，cp₃₅[SourceCS]为“104”，cp₃₅[CSTrend]为“43”；

步骤24：用裁剪面要素集合CutP中的41个面要素与原始地层面要素OriS中的138个面要素分别进行求交，得到场景地层面要素集合SceS＝{ss_k|k＝0，1，2，...，839}，绘制图层如图5所示。以裁剪面cp₃₅为例，其与12个原始地层面要素os₀、os₈₉、os₉₀、os₉₆、os₁₀₃、os₁₀₄、os₁₀₅、os₁₂₃、os₁₂₄、os₁₂₅、os₁₂₉和os₁₃₀相交，求交并进行拆分多部件要素后得到16个场景地层面要素{ss₇₂₂,ss₇₂₃，ss₇₂₄，ss₇₂₅，ss₇₂₆，ss₇₂₇，ss₇₂₈，ss₇₂₉，ss₇₃₀，ss₇₃₁，ss₇₃₂，ss₇₃₃，ss₇₃₄，ss₇₃₅，ss₇₃₆，ss₇₃₇}；

步骤25：对SceS内的所有要素属性进行赋值，如面要素ss₇₂₂为cp₃₅与os₀相交得到的，则其属性赋值如下：

(三)识别结构建模与化简

步骤31：根据面要素集合SceS＝{ss_i|i＝0，1，2，...，839}中要素的ss_k[Gravity]属性绘制点图层，得到地层对象点集SceG＝{sg_i|i＝0，1，2，...，839}。对象点sg_k继承场景面要素ss_k的所有属性；

步骤32：对于同一识别条带内邻接的地层面要素ss_i与ss_j(j>i)，需要连接对应的对象点sg_i与sg_j。以要素ss₇₂₁、ss₇₂₂、ss₇₂₃、和ss₇₂₄为例(如图6所示)，依据规则进行如下判断：

1)如面要素ss₇₂₁与ss₇₂₂，ss₇₂₁[SourceCS]与ss₇₂₂[SourceCS]不相等，则不符合条件；

2)如面要素ss₇₂₂与ss₇₂₃满足ss₇₂₂[SourceCS]＝ss₇₂₃[SourceCS]，且面要素ss₇₂₂与ss₇₂₃有重合边界点，则连接对象点sg₇₂₂与sg₇₂₃，得到线rl₇₃₀<sg₇₂₂[Gravity],sg₇₂₃[Gravity]>；

3)如面要素ss₇₂₂与ss₇₂₄满足ss₇₂₂[SourceCS]＝ss₇₂₄[SourceCS]，但面要素ss₇₂₂与ss₇₂₄没有重合的边界点，则不符合条件；

步骤33：完成对SceS的遍历，绘制线图层，得到关系线集合RelL＝{rl_k|k＝0，1，2，...，847}，共有847对在识别场景条带中邻接的地层对象，线要素代表对象间的关系。对关系线要素rl_k的属性进行赋值，以属性[SourceCS]为104的地层对象为例，其属性如下表：

连接上表间的邻接对象，对关系线要素属性赋值如下表：

结合对象与关系的图层与属性信息，完成“地质分区”图层中识别场景条带上邻接关系属性模型的构建，点线图层如图7所示；

步骤34：基于地层对象点集合SceG和关系线集合RelL构建无向邻接集合SimAList＝{sali|i＝0，1，2，...，839}，集合sali中存储对象sgi邻接的对象的属性[Id]序列。以上述表所示的对象与关系为例，可得到：

sal₇₂₂＝{723，725，732}； sal₇₂₃＝{722，731}；

sal₇₂₄＝{732，733，734}； sal₇₂₅＝{722，732}；

sal₇₂₆＝{735，737}； sal₇₂₇＝{728，737}；

sal₇₂₈＝{727，731}； sal₇₂₉＝{732}；

sal₇₃₀＝{733}； sal₇₃₁＝{723，728}；

sal₇₃₂＝{722，724，725，729}； sal₇₃₃＝{724，730}；

sal₇₃₄＝{724}； sal₇₃₅＝{726，738}；

sal₇₃₆＝{735}； sal₇₃₇＝{726}；

步骤35：借助无向邻接集合，根据条件对对象点集进行化简，得到的化简后地层对象点集SimSG＝{ssg_g|g＝0，1，2，...，369}。以[SourceCS]为104的条带为例，对其邻接关系属性模型进行化简后，得到化简后对象点属性如下：

基于表数据绘制化简后的地层对象点图层；

步骤36：化简后的示例无向邻接集合展示如下：

sal₇₂₂＝{725，728}； sal₇₂₃＝{}； sal₇₂₄＝{}；

sal₇₂₅＝{722}； sal₇₂₆＝{735，737}； sal₇₂₇＝{728，737}；

sal₇₂₈＝{722，727}； sal₇₂₉＝{}； sal₇₃₀＝{}；

sal₇₃₁＝{}； sal₇₃₂＝{}； sal₇₃₃＝{}；

sal₇₃₄＝{}； sal₇₃₅＝{726，736}； sal₇₃₆＝{735}；

sal₇₃₇＝{726}；

步骤37：根据删改后的无向邻接集合绘制化简后的线图层，得到关系线要素集合SimRL＝{srl_k|k＝0，1，2，...，307}。如邻接集合sal₇₂₂，包含元素“725”和“728”。则连接对象点sg₇₂₂与sg₇₂₈得到关系线srl₂₆₀，连接对象点sg₇₂₂与sg₇₂₇得到关系线srl₂₆₁，再如步骤33重新赋值属性。示例条带化简后的对象间关系如下表所示，其中关系线srl₂₅₅表示的为化简后新增的关系：

化简后的邻接关系属性模型点线图层如图8所示；

(四)褶皱构造模式匹配

步骤41：基于化简后的地层对象点集SimSG＝{ssg_j|j＝0，1，2，...，369}与关系线集SimRL＝{srl_k|k＝0，1，2，...，307}，进行褶皱构造识别。通过查找ssg_j[SourceCS]＝ssg_j[SourceId]的对象点获取核部地层对象点加入点集CorG，得到核部地层对象点集CorG＝{cg_i|i＝0，1，2，...，40}，并创建对应的褶皱对象集合OriFold＝{of_i|i＝0，1，2，...，40}。以上述表所示的对象为例，其中，对象点ssg₃₁₆满足[SourceCS]属性与[SourceId]属性相等，cg₃₅即为ssg₃₁₆，建立对应的褶皱对象of₃₅；

步骤42：将每个核部地层对象点cg_i的信息加入对应的褶皱对象of_i中，对属性of_i[SourceCS]赋值cg_i[SourceCS]；of_i[IdList]中添加cg_i[Id]。以上述表所示的对象与关系为例，核部地层对象点为cg₃₅，根据cg₃₅的属性[SourceCS]和[Id]给of₃₅赋值，of₃₅[SourceCS]赋值104，of₃₅[IdList]数组中添加727；

步骤43：遍历化简后地层对象点集SimSG与关系线集SimRL，将规则匹配的褶皱对象记录存放进OriFold集合中，以上述表所示的对象与关系为例，将cg₃₅的[Id]属性赋值给索引编号“LeftS”与“RightS”，LeftS＝727，RightS＝727。基于此进行递归判断如下：

1)存在关系srl₂₆₄，满足srl₂₆₄[StartId]＝LeftS，存在关系srl₂₆₅，满足srl₂₆₅≠srl₂₆₄，且srl₂₆₅[StartId]＝RightS；

2)srl₂₆₄[TimeDis]与srl₂₆₅[TimeDis]相等，都为1；

3)|srl₂₆₄[AzimuthAng]-srl₂₆₅[AzimuthAng]|＝179.68，|179.68-180|＝0.32<15；

4)则关系srl₂₆₄与srl₂₆₅满足背斜规则集，of₃₅[FoldType]赋值“Anticline”；将srl₂₆₄[EndId]与srl₂₆₅[EndId]分别插入of₃₅[IdList]的第一个位置和最后一个位置，得到of₃₅[IdList]＝{728,727,737}；将728和737分别赋值给索引编号“LeftS”与“RightS”，继续判别；

5)存在关系srl₂₆₁，满足srl₂₆₁[StartId]＝LeftS，存在关系srl₂₆₃，满足srl₂₆₃≠srl₂₆₁，且srl₂₆₃[StartId]＝RightS；

6)srl₂₆₁[TimeDis]＝2，srl₂₆₃[TimeDis]＝1，两者不相等，不满足判别规则，退出判别；

基于每个cg_i重复规则判别，修改褶皱对象of_i的属性；

步骤44：对于初始的再找对象集合OriFold＝{of_i|i＝0，1，2，...，40}中的每个褶皱对象of_i要素，存在of_i[IdList]的长度小于3，将of_i从OriFold中移除，最终得到褶皱单元Fold＝{f_i|i＝0，1，2，...，8}，分别对应属性[SourceCS]为34、64、65、66、94、101、102、103和104的识别场景条带。根据数组f_i[IdList]中的值查找属性与之相等的化简前地层对象点sg_j，将sg_j[SourceId]的值添加到f_i[SIList]中，对褶皱单元f_i的属性赋值如下：

步骤45：对褶皱单元f_i的重复部分进行处理，示例褶皱单元f₅与f₇的对象范围修改如下：

步骤46：对褶皱单元f_i的其余属性赋值如下：

(五)褶皱地貌类型确定

步骤51：若未加载栅格数据，则加载高程栅格数据“lushanDEM.bgd”；

步骤52：对于褶皱单元f_i，找到属性[Id]为f_i[IdList][0]、f_i[IdList][ListCount/2]和f_i[IdList][ListCount-1]的地层对象点，根据其[Gravity]属性，从栅格数据上获取高程分别为L、C和R。对于示例褶皱单元f₈，分别获取对象点sg₇₂₇、sg₇₂₈与sg₇₃₇的高程为L、C和R。则，L＝321.1m，C＝409.2m，R＝730.3m；

步骤53：L、C和R不满足L-C＞0&R-C＞0，则f₈[Landform]赋值“Mountain”；

步骤54：以f_i[Gravity]为中心点，f_i[a]为长轴，f_i[b]为短轴，f_i[Trend]为旋转角度，绘制椭圆表示褶皱地貌单元。以褶皱单元f₈为例，在点(402538.610,3277066.956)处，以10900.1m为长轴，1235m为短轴，逆时针旋转43°，得到代表褶皱单元f₈的地貌单元。基于所有褶皱单元f_i的属性绘制地貌单元，得到褶皱地貌类型图层，如图9所示。

Claims (5)

1.一种褶皱地貌类型的自动化识别方法，其特征在于，包括如下步骤：

1)根据核部地层与横切剖面线提取原则，提取可能存在褶皱的场景条带；

2)基于邻接属性关系图建模方法，对场景条带进行场景空间结构的建模与化简；

3)通过判别化简后的邻接属性关系模型能否由对应的褶皱规则推导而来，实现针对不同褶皱构造模式的匹配识别；

4)基于识别出的褶皱场景进行山和谷形态的进一步判别，最终确定褶皱构造地貌类型，粗略划定褶皱地貌单元。

2.根据权利要求1所述的一种褶皱地貌类型的自动化识别方法，其特征在于，该方法的具体实现过程为：

(一)数据预处理

步骤11：加载shp格式的地质面图层数据，得到所有地层面要素集合OriS＝{os_i|i＝0，1，2，...，n-1}，n为地层面要素的个数，面要素os_i的属性包括原地层编号[SourceId]、地层年代[Age]、走向[Tend]、重心点坐标[Gravity]、岩类[RockType]、出露面积[Area]和出露长度[Length]；

步骤12：对面要素os_i的属性赋值：要求地层年代[Age]属性不为空，走向[Tend]属性为极坐标方位角；若面要素os_i属性不完整，对其进行补充；其中，岩类[RockType]属性结合事先生成的XML格式的地层年代序列表进行赋值，并根据需要标记为化简地层与保留地层两种状态；

(二)识别条带的提取

步骤21：获取集合OriS＝{os_i|i＝0，1，2，...，n-1}中每个面要素os_i的邻接地层面要素集合NeiS_i，并从集合NeiS_i中移除状态标记为化简地层的元素，得到面要素os_i的邻接地层面要素集合NeiS_i＝{ns_ij|j＝0，1，2，...，l_i-1}，l_i为面要素集合NeiS_i中状态为保留地层的元素数量；

步骤22：对于状态为保留地层的面要素os_i，判断其邻接地层面要素集合NeiS_i是否满足以下条件：

若满足以上任一条件，则判定面要素os_i为核部地层，将其加入核部地层要素集合CorS＝{cs_j|j＝0，1，2，...，p-1}，p为查找到的核部地层的数量；

步骤23：针对每个核部地层要素cs_j，过其重心点坐标[Gravity]垂直于走向[Tend]做缓冲区，得到裁剪面要素集合CutP＝{cp_j|j＝0，1，2，...，p-1}；裁剪面cp_j记录属性包括核部地层编号[SourceCS]和核部地层走向[CSTrend]；

步骤24：用裁剪面要素集合CutP对集合OriS进行裁剪，得到场景地层面要素集合SceS＝{ss_k|k＝0，1，2，...，s-1}，s为识别条带上地层对象的数量；由裁剪面cp_j与面要素os_i裁剪得到的面要素ss_k，其出露面积[Area]、出露长度[Length]和重心点坐标[Gravity]重新计算，编号[id]属性赋值为k，其余属性继承自裁剪面cp_j与面要素os_i；

(三)识别结构建模与化简

步骤31：由上一步获得地层面要素集合，根据面要素ss_k的重心点坐标[Gravity]属性新建点，得到地层对象点集SceG＝{sg_i|i＝0，1，2，...，s-1}，对象点sg_i继承面要素ss_k的所有属性；

步骤32：针对对象点集，基于邻接属性关系图的方式建模，对满足空间拓扑关系为邻接的地层对象间建立空间方位关系与时间距离关系，并绘制线要素，得到关系集合RelL＝{rl_k|k＝0，1，2，...，t-1}，t为条带上的邻接关系数量；关系属性记录为：起点对象Id[StartId]、终点对象Id[EndId]、空间方位关系[AzimuthAng]和时间距离关系[TimeDis]，定义[StartId]为较老的地层对象的编号，[EndId]为较新的地层对象的编号，[AzimuthAng]为起点对象指向终点对象的方位角，[TimeDis]为邻接对象在地层年代序列表中对应节点的[Id]值的差值的绝对值；

步骤33：根据关系集合RelL建立无向邻接集合SimAList＝{sal_i|i＝0，1，2，...，s-1}；集合sal_i代表属性[Id]值为i的对象点所邻接的对象点的[Id]值集合；

步骤34：遍历判断对象点sg_i的原地层编号[SourceId]与核部地层编号[SourceCS]属性是否相等，若不相等，则进行化简判断；

步骤35：遍历无向邻接集合SimAList，对于化简后状态为化简地层的对象点sg_i，清空集合sal_i使其长度为0，并将其[Id]值从其他集合中移除；对于状态为保留地层的对象点sg_i，若集合sal_i的长度大于2，则保留与对象点sg_i形成的方位角最接近垂直于核部地层走向[CSTrend]的两个对象，将其它对象的[Id]值从集合sal_i中移除；

步骤36：由修改后的无向邻接集合SimAList得到化简后点集SimSG＝{ssg_j|j＝0，1，2，...，c-1}和化简后的关系集合SimRL＝{srl_k|k＝0，1，2，...，r-1}，其中，c为化简后保留的地层要素的数量，r为化简后的邻接关系数量；

(四)褶皱构造模式匹配

步骤41：创建集合OriFold＝{of_i|i＝0，1，2，...，p-1}存放初始褶皱对象，其数量与核部地层对象数量相等；褶皱判别对象of_i的属性继承自核部地层对象cg_i，并新增属性如下：褶皱地层的编号数组[IdList]、褶皱原始地层的编号数组[SIList]、褶皱单元长轴[a]、褶皱单元短轴[b]、褶皱类型[FoldType]、地貌类型[Landform]；

步骤42：基于化简后的点集SimsG与关系集合SimRL，进行褶皱构造识别；先获取化简后的核部地层点集合CorG＝{cp_i|i＝0，1，2，...，p-1}；一个cg_i对应一个褶皱判别对象of_i；

步骤43：判别关系集合SimRL中是否递归存在关系srl_a与关系srl_b总满足背斜规则集或向斜规则集；

步骤44：对于规则匹配的褶皱判别对象of_i，将其加入褶皱单元集合Fold＝{f_i|i＝0，1，2，...，q-1}，q为符合褶皱规则的褶皱单元数量；对褶皱单元f_i属性赋值：褶皱类型[FoldType]取决于其匹配的规则集，褶皱地层的编号数组[IdList]记录组成褶皱的对象点编号，褶皱原始地层的编号数组[SIList]记录组成褶皱的原始地层编号，其余属性继承自对应的核部地层对象；

步骤45：对集合Fold中两两褶皱的重复部分进行舍取，根据去重复后的褶皱单元f_i的范围，对属性褶皱单元长轴[a]和褶皱单元短轴[b]进行赋值；

(五)褶皱地貌类型确定

步骤51：加载同一投影坐标系下的DEM数据；

步骤52：对于元素f_i，从栅格数据上获取褶皱范围最左、核部与最右对象点的高程分别为L、C和R；

步骤53：若L-C＞0&R-C＞0，则地貌类型为谷；反之，地貌类型为山；绘制椭圆表示褶皱地貌单元。

3.根据权利要求2所述的一种褶皱地貌类型的自动化识别方法，其特征在于，步骤12中，若面要素os_i属性不完整，缺失的走向[Tend]属性通过计算面要素os_i多边形的走向来代替。

4.根据权利要求2所述的一种褶皱地貌类型的自动化识别方法，其特征在于，步骤34中，对对象点sg_i进行化简判断及处理的具体方案如下：

1)对于地层对象点sg_i与其邻接对象点sg_j，其中j＞i，若两者都为保留地层且地层年代[Age]属性在地层年代序列表中编号[Id]相等，或两者都为化简地层，则对对象点sg_i与sg_j进行如下的合并操作：

a)以对象点sg_i与sg_j的出露面积[Area]属性为权重重新计算重心点坐标[Gravity]，保留对象点sg_i的其他属性；

b)对象点sg_j的状态标记为化简地层；

c)将集合sal_j整合到集合sal_i中，清空集合sal_j使其长度为0；

2)若地层对象点sg_i为化简地层，且其邻接地层对象中包含至少两个状态为保留地层的对象点sg_u和sg_v，其中v＞u：

a)若对象点sg_u与sg_v满足：①对象点sg_u与sg_v的地层年代[Age]属性在地层年代表序列表中对应的节点的编号[Id]相等，②对象点sg_u指向sg_v的方位角ang_uv或对象点sg_v指向sg_u的方位角ang_vu，与对象点sg_i的核部地层走向[CSTrend]差值小于容限，则执行步骤1)中的a)至c)对对象点sg_u与sg_v进行合并操作；

b)若对象点sg_u与sg_v满足：①在对象点sg_u指向对象点sg_v的方位角ang_uv的角度阈值与对象点sg_v指向对象点sg_u的方位角ang_vu的角度阈值范围内不存在其他状态为保留地层的对象，②对象点sg_u与sg_v不是邻接对象，则将对象点sg_v的[Id]值从集合sal_i中移除并加入集合sal_u，将对象点sg_u的[Id]值从集合sal_i中移除并加入集合sal_v；

3)不再出现以上情况，则清空集合sal_i使其长度为0。

5.根据权利要求2所述的一种褶皱地貌类型的自动化识别方法，其特征在于，步骤43中，背斜规则集如下：

向斜规则集如下：

其中，索引编号LeftS与RightS分别为当前往左方向与右方向上满足规则的关系的终点对象的编号，以及将要查找的关系的起点对象的编号[Id]；tolerance为角度容限。