一种地质图专题信息提取方法

一种地质图专题信息提取方法
【专利摘要】本发明根据图例信息，针对地质图上各种要素颜色、形状、纹理、空间关系等特征的不同，采用颜色分割、边界追踪、模式识别等技术，提取出地质图中所蕴含的与矿产有关的主要信息，如地层、岩浆岩、构造、矿点等信息，形成专题信息图层。其中地质图专题信息提取方法包括：地质图数据的预处理方法，地质图专题信息提取方法，信息输出方法。解决现有地质图信息提取的难点问题，即以半结构化的标准地质图为研究对象，以地质图信息提取为目的，结合图例信息，提取出地质图中所蕴含的与矿产有关的主要信息，如地层、岩浆岩、断层、矿点等信息，构建专题信息图层，形成一套半结构化数据(图像)转换为结构数据的技术。
【专利说明】
_种地质图专题信息提取方法
技术领域
[0001] 本发明涉及一种地质图专题信息提取方法，尤其涉及一种自动提取出栅格格式 (JPG、TIFF等)地质图中所蕴含的与矿产有关的主要信息，并形成矢量的专题信息图层的方 法。
【背景技术】
[0002] 专题信息是指地质图中的某一类信息，如地层、岩浆岩体、蚀变带、断层、地质界 线、矿点等。
[0003] 图像几何校正是指通过一系列的数学模型来改正和消除由于各种因素引起的原 始图像上各地物的几何位置、形状、尺寸等特征与在参照系统中的表达要求不一致时产生 的变形(如缩放、旋转、平移、仿射等几何畸变），得到正射影像或者近似正射影像的过程。
[0004] 图像去噪是现实中的数字图像在数字化和传输过程中常受到成像设备与外部环 境噪声干扰等影响，称为含噪图像或噪声图像。减少数字图像中噪声的过程称为图像去噪。
[0005] 颜色分割是由于扫描后的彩色地质图同一颜色的RGB不完全相同，通过颜色分割 使相近颜色融合，产生大色块的效果。模糊C-均值算法(FCM)就是一种颜色分割方法。
[0006] 拓扑关系是指满足拓扑几何学原理的各空间数据间的相互关系。即用结点、弧段、 多边形和岛所表示的实体之间的邻接、关联、包含和连通关系。
[0007]模式识别是指对表征事物或现象的各种形式的(数值的、文字的和逻辑关系的)信 息进行处理和分析，以对事物或现象进行描述、辨认、分类和解释的过程。
[0008] 图像匹配就是把不同的传感器或者同一传感器在不同时间、不同成像条件下对同 一景物获取的两幅或多幅图像进行对比，从而找到该组图像中的共有景物，或者是根据已 知模板图像，在另一幅图中寻找相应图像的过程。序贯相似性检测算法就是图像匹配的一 种方法。
[0009] 数字地球就是利用信息技术，将地球表面上每一点的特征信息数字化。地质图是 表示地质现象及构造特征的专题地图，反映了地质现象的空间展布关系。它是一个区域地 质研究的重要成果，融合了丰富的地质专家知识，能够反映一个地区整体的地貌地质特征。 大数据时代的到来使思维发生改变，可通过大数据分析挖掘方法，快速将半结构化地质图 件，转变为结构化数据，用于分析(如成矿预测分析）。
[0010] 由于受多种因素的影响，扫描地形图图像的质量千差万别。对于质量清晰，没有退 化变形的图像来说，准确的进行地图跟踪矢量化及符号识别，其难度将大大降低。但是对于 质量不甚理想，退化变形严重的图像来说，理想状况下的追踪识别方法将失去功效。目前， 大部分地图要素识别算法研究都是基于二值分色图像展开了，而如何得到清晰完整的二值 分色图像，正是彩色地图识别的一个难点。其次，大部分算法在设计时都没有充分考虑地图 的复杂性，忽视了地图上各种图形文字混杂、不同要素相互覆盖，甚至粘连、断裂的情况，导 致算法不够稳健，不具有鲁棒性。
[0011] 总之，现在的地图信息提取技术或者矢量化，都是针对地形图展开研究的，还没有 人做过关于地质信息提取方面的研究。针对栅格格式的地质图，如果应用传统技术，即人工 将这些图件矢量化，变成结构化数据，这无疑是一项巨大的工程，需要花费大量的时间、精 力及金钱。而现在的地图信息提取技术都是针对地形图展开研究的，与地质图颜色的多样 性相比，地形图最多只有四色，即为红、黄、青、黑四种颜色，所以地形图的分色(颜色分割） 更容易，信息提取也就更为方便。市面上现在的矢量化软件，也只能处理比较简单的图像， 对于地质图中各种要素相互叠压覆盖，矢量化效果很差，根本无法达到实用的效果。

【发明内容】

[0012] 本发明为了解决现有地质图信息提取的难点问题，即以半结构化的标准地质图为 研究对象，以地质图信息提取为目的，结合图例信息，提取出地质图中所蕴含的与矿产有关 的主要信息，如地层、岩浆岩、断层、矿点等信息，构建专题信息图层，形成一套半结构化数 据(图像)转换为结构数据的技术。
[0013] 本发明所述的标准地质图主要是包括主图、图例等部分的地质图，图例与地质图 具有完备性和一致性的原则，即地质图上各种符号和颜色所代表内容与指标在图例中都有 所说明。
[0014] 本发明的主要内容为:根据图例信息，针对地质图上各种要素颜色、形状、纹理、空 间关系等特征的不同，采用颜色分割、边界追踪、模式识别等技术，提取出地质图中所蕴含 的与矿产有关的主要信息，如地层、岩浆岩、构造、矿点等信息，形成专题信息图层。
[0015] 本发明的地质图专题信息提取方法包括:第一，地质图数据的预处理方法，包括图 像几何校正、彩色图像去噪、主图区域和图例区域定位;第二，地质图专题信息提取方法:包 括提取面要素信息、提取线要素信息、提取点要素信息，构建各专题信息图层;第三，信息输 出方法:将提取的信息与GIS软件相连，形成矢量化的专题信息图层。
[0016] (1)地质图像的预处理方法
[0017] 对彩色地质图像进行预处理，通过几何校正对失真图像进行的复原性处理;通过 自适应中值滤波方法对扫描地质图进行彩色图像去噪，在减少图像中的噪声的同时尽可能 的保留图像细节，不对图像引起模糊现象;根据地质图的二值化后区域连通性特点确定地 质图主图区域和图例区域的定位。通过对地质图数据进行预处理，为下一步专题信息的提 取作好准备。
[0018] (2)地质图专题信息提取方法
[0019] 将经过预处理的地质图，分别同时地进行提取面要素信息、提取线要素信息、提取 点要素信息，即面要素信息、线要素信息、点要素信息的提取是并行完成的。
[0020] 第一，提取面要素信息，即提取地层、蚀变带、岩体等各种面要素信息，包括:边界 白色化处理、颜色分割、连通体构建、颜色分割后的小色斑消除的处理、边界处白色区域修 复。
[0021] 第二，提取线要素信息，根据线要素信息的特点，将其分为两类，一类是地质边界、 岩体边界、蚀变带边界等面要素区域边界线信息提取;另一类是非区域边界线信息提取，如 断层线等。地质边界、岩体边界、蚀变带边界等面要素区域边界线通过区域边界线的追踪及 拓扑关系的建立、区域轮廓线的矢量化及光滑、断层线的提取，而断层线通过颜色特征来提 取。
[0022]第三，点要素信息的提取方法，对于矿点，在标准地质图中，图中所出现的矿点都 会在地质图图例中用相同（包括大小、颜色等）的符号加以说明，所以首先通过自动提取出 图例中的矿点符号，然后以提取的矿点符号为模板，图像匹配，在地质图中找到相应矿点的 位置并定位，最后在相应位置生成矿点，形成矿点信息图层。
[0023] (3)信息输出方法
[0024]将提取的面要素信息、线要素信息、点要素信息图层与现有GIS软件相结合，形成 矢量化的专题信息图层。
【附图说明】
[0025] 图1整体技术路线图
[0026]图2彩色图像的自适应滤波流程图 [0027]图3面要素信息提取的流程图 [0028]图4颜色分割步骤流程图 [0029]图5面要素区域边界追踪流程图 [0030]图6断层线提取的流程图
【具体实施方式】
[0031] 一地质图图像预处理方法
[0032]地质图图像预处理方法是指地质图图像的几何校正、彩色图像去噪以及主图区域 和图例区域定位等预处理方法。
[0033] 具体为：
[0034] 1几何校正
[0035] STEP 1:对于畸变图像中的任意一点（Xi，yi)，根据原图坐标系和畸变图像坐标系 [xu} =k{x.\\) 之间的关系式j f 、.，_找出对应畸变图像中的对应点（11^，71^) = [111(11，71)，112(叉1， yi)]；
[0036] STEP 2:由于通常情况下，（XUl，ym)并不是整数点，所以采用双线性插值方法找出 和(xm，yUl)最接近的整数点(^乂山根据畸变图像⑴^丨处^邻域内的斗个像素点的灰 度值，采用双线性插值方法计算出该点的灰度值g (Xi，yi)的计算公式 容(wX1-([工+办]，~).'([工+知]+1，[少+办])其中&- + + Aj(1 - Ax) f{[x + Sx],[y + Sy] +1) + AxAyf ([x + &] +1,[y + S)i] +1) '、 x_[x+5x], A y = y+y-[y+5y]。
[0037] 2彩色图像去噪
[0038] STEP 1:对彩色数字图像进行分解，即将彩色图像分解为R、G、B三个分量；
[0039] STEP 2:分别对R、G、B三个分量的图像进行灰度图像的自适应中值滤波处理；
[0040] STEP 3:将三个分量的去噪图像进行组合，得到的图像就是去噪后的彩色图像。 [0041] 3主图区域和图例区域定位是将二值化后地质图中最大的连通区域确定为地质图 的主图区域范围。
[0042]二地质图专题信息提取方法
[0043]将经过预处理的地质图，并行地进行面要素信息、线要素信息、点要素信息的提 取。
[0044] 1提取面要素信息
[0045]第一，对经过预处理的地质图进行边界白色化处理，包括:①采对地质图进行边缘 查找;②对边缘图像进行一次粗化处理;③将原图像与边界图像相减，从而在原图中，将边 界处全部进行白色化处理。其目的是为了减少边缘的过渡色、混淆色、色散像素，同时也减 少聚类数目，以提高颜色分割的效果。具体为：
[0046] A边缘查找，具体步骤为： _ f.(-, , \
[0047] STEP 1:根据梯度向量公式丽g(明=+ (%):和梯度方向公式+，,0 = arctan ? , 计算每一点(x，y)处的局部梯度mag(Vf)及边缘方向a(x，y)，从而确定边缘点；
[0048] STEP 2:确定高阈值TdP低阈值TsmST〗），用这两个阈值对STEP 1中出现的边缘 点做阈值化处理，其中，大于1^的边缘象素称为强边缘象素，TjPT22间的边缘象素则称为 弱边缘象素；
[0049] STEP 3:最后，将邻域8-连接的弱边缘象素加到强边缘象素中，并进行边缘连接， 得到边缘图像。
[0050] B粗化处理，对边缘查找的图像中的背景进行细化处理，然后通过对细化结果求补 来达到对边缘图像的粗化处理。
[0051] 第二，对经过边界白色化处理的地质图进行颜色分割，具体步骤为：
[0052] STEP 1:相关参数的初始化。根据图例的中面要素区域的颜色来初始化聚类的类 别数c及聚类中心￥(()) = {^^2丄^。};设定迭代停止的阈值6，初始化迭代计数器为1 = 0和 最大迭代次数lmax = 100;设加权指数m = 2。
[0053] STEP 2:当迭代次数1 < 100时，根据聚类中心V(1) = {V1，v2，L，Vc}计算隶属度矩阵U (1)。如果禮*0(即^^1|))，则利用公式
计算像素点pk到所有聚类中 心的if(i = l，2，L，C)的隶属度，得到隶属度矩阵Uw。其中，dik(j = l，2，L，n;k=l，2，L，n) 为第k个像素点pk到第i个聚类中心#(丨=1，2丄，(3)的距离。如果，以=0(即则 u][] =1 pk=v、:、 n ^ a 。在Lab颜色空间中，定义任意两个像素点颜色的欧式距离为： wJ/LO 其13且jVk
[0054] STEP 3:对图像的颜色重新进行聚类划分。在图像的隶属度矩阵U(1)中，从每列中 选择隶属度最大的点作为对应点的归属域，根据公式
重新计算每个类 的聚类中心va+1)。
[0055] STEP 4:如果| |V(1+1)_V(1)| |<e，则迭代停止并输出最终的划分矩阵U和聚类中心 V，否则令1 = 1+1，转向执行STEP 2。
[0056] STEP 5:去模糊化。根据模糊划分矩阵U提供的像素对各聚类中心的隶属度来确定 该像素点最后的归属。采用最大隶属转换方法，用C k表示第k个像素点的归属类别，SP Q. =arg;{max{//,/r}} V/,VA: 〇
[0057]第三，对颜色分割后的地质图进行连通体构建
[0058]对颜色分割后的地质图图像，根据颜色可形成单色版的图层，然后进行扫描串的 构建，接着根据扫描串的邻接关系进而形成连通体，并计算每个连通体的面积和黑白比（即 连通体的最小外接矩形内前景像素个数与背景像素个数之比）。
[0059]第四，分析颜色分割后噪声的特点，以连通体为基本处理单元，针对噪声区域与非 噪声区域两类连通体特征的不同来消除小色斑(即噪声）：具体步骤为：
[0060] STEP 1:通过面积特征去除全部的孤立点状噪声和小面积片状噪声以及部分条带 状噪声；
[0061] STEP 2:利用黑白比去除余下的条带状噪声；
[0062] STEP 3:得到无噪声的面区域图像，完成颜色分割后面区域噪声的消除。
[0063]第五，对小色斑消除后的地质图进行边界处白色区域修复，具体步骤为：
[0064] STEP 1:扫描一次整张地质图的白色像素，标记出某个白色像素的8-邻域内有两 个或两个以上为前景像素(非白色）的像素点；
[0065] STEP 2:将该像素点的颜色置换为4-邻域内颜色最多的那个颜色；
[0066] STEP 3:重复STEP 1和STEP 2,层层修复，直到所有白色区域完全被修复为止。 [0067] 2提取线要素信息
[0068]地质图上的线要素基本可以分为两类，一类是地质边界、蚀变带边界等面要素区 域边界线;另一类是非区域边界线，如断层线等。
[0069] (1)提取区域边界线信息
[0070]第一:提取边界点和结点
[0071]根据地质图图像中相邻四个栅格单元颜色值的异同来判断交汇点处是边界点还 是结点。相邻四个单元栅格颜色值的不同，分为三类15种情况:不存在边界点或结点、只可 能存在边界点、只可能存在结点。逐行分析4个栅格单元交汇点的类型，快速地提取出图像 中所有的坐标点，并记录他们的行列坐标、类型以及连接信息。
[0072] 第二:跟踪弧段
[0073]在栅格图像中，存在且只存在这两种类型的弧段:开放弧段和封闭弧段。封闭弧段 是只包含有边界点的弧段，余下的为全为开放弧段(至少包含有两个结点的弧段）。按照先 跟踪开放弧段然后再跟踪封闭弧段的顺序来处理。
[0074] A跟踪开放弧段，开放弧段的跟踪从一个结点开始，具体步骤如下：
[0075]在整幅图像中找到一个未标记过(没有被跟踪过)的结点，从该结点的0方向开始， 根据边界点的连接信息进而可以跟踪得到另一个新的边界点，如此下去，当跟踪到某个结 点时，说明该弧段已经跟踪完毕。同理，依次对该结点的1、2、3三个方向上的弧段进行跟踪， 从而完成该结点处四个方向上全部弧段的跟踪生成。在跟踪完一条弧段后，为了避免后面 在跟踪其它弧段时，对已跟踪弧段的重复跟踪，应该清除此弧段上所有边界点的连接信息， 使这些边界变成无效点，即表示不会再有弧段与它们相连接。如果该弧段只由两个的结点 组成，那么就清除该方向上的连接信息即可，使此方向上变成无邻接点。当一个结点上下左 右四个方向上的弧段都跟踪完成后，该结点也变成为无效点。
[0076] B跟踪封闭弧段，跟踪完成全部开放弧段之后，剩下的有效点都为边界点（且只有 边界点）。由于封闭弧段是一个闭合的弧段，从有效的且ID最小的边界点开始跟踪，弧段跟 踪完成之后必定回到该点，完成该条封闭弧段的生成。对跟踪完的封闭弧段的边界点，要及 时清除连接信息。
[0077] 第三:组建多边形
[0078]在建立弧段和结点的连接关系的基础上，沿弧段方向，从而组建多边形。
[0079] A建立结点和弧段之间的连接关系，即显式地记录每一个结点与哪几条弧段有连 接关系。一个结点最多能与4(即右、上、左、下4个方向）条弧段相连接。实际上，在跟踪弧段 的同时就建立起了结点与弧段的连接信息。
[0080] B由弧段组建多边形，利用弧段组建多边形，分成两类:开放弧段和封闭弧段组建 多边形。
[0081] 首先，开放弧段组建多边形。开放弧段为两个多边形之间的公用边界，沿弧段方 向，两边面域的属性不同，所以，沿弧段方向，采用"属性匹配"搜索多边形的连接弧段，可以 得到两个不同的多边形。在此过程中，多边形的邻接关系同时也会被记录。当一个多边形组 建完成后，应该给这个多边形赋上相应的属性值。在开放弧段组建多边形的过程中应该注 意:①弧段方向的左边=逆弧段方向的右边，弧段方向的右边=逆弧段方向的左边;②每条 开放弧段都具有两次搜索机会，即分别采用沿弧段方向左、右两边的属性为匹配模板进行 搜索弧段。因此，每条弧段具有两个搜索标志（即弧段方向的左边搜索标志和右边搜索标 志）。在组建多边形的过程中，每生成一个新的多边形后，对组成该多边形的全部弧段，更新 它们的左边搜索标志或右边搜索标志，具体做法为:在一个多形组建过程中，若使用左边属 性连接弧段，则该弧段的左边搜索标志置为false;若使用右边属性连接弧段，则该弧段的 右边搜索标志置为false。当一条弧段的左右两边搜索标志都为false时，依次轮到搜索该 弧段时，则跳过此弧段，不再对它进行搜索。这种设置双向标志的作用在于保证每条弧段只 使用两次，且不会产生重复多边形，避免判断重复搜索，加快了处理速度。
[0082] 采用"属性匹配"法搜索完开放弧段后，再找出图中的封闭弧段来组建多边形。由 于封闭弧段全部由边界点组成，且自身封闭，所以只需要找到封闭弧段中ID最小的边界点， 从此点开始，以弧段右边的属性为匹配模板搜索一次即可。
[0083]第四：建立拓扑关系
[0084] 根据地质图的特点，拓扑关系的建立要存储多边形之间的相邻和包含关系。在建 立多边形的相邻关系之前首先要解决多边形自交问题;在存储多边形的包含关系前要解决 "岛"问题，建立复杂的"岛"多边形。具体方法包括；
[0085] 第一步骤:多边形自交处理方法，在追踪弧段组建多边形时，若遇到对角相同的情 况，当采用弧段方向左边(或右边)的属性为匹配模板搜索到对角点〇时，就向左边(或右边） 转90度后继续搜索下一条弧段。解决多边形的自交问题。
[0086] 第二步骤："岛"现象的处理方法，需要处理共边岛问题和多重嵌套岛问题。
[0087] A在处理共边岛问题时，根据图论的思想，将坐标点（边界点和结点）当作顶点，而 弧段当作边。如果能从一个多边形的任一个顶点出发，不经过组成多边形的弧段而到达另 一个顶点，这就说明此多边形是多余的，应该被删除。通过这种方法就可以完全地解决共边 岛问题。
[0088] B在解决共边岛问题后，要建立多边形之间的包含关系。为了提取处理效率，标记 出图像中的所有连通体，并检查每个连通体中是否存在岛屿多边形。采用此方法建立多边 形的包含关系不用对多重嵌套岛进行特殊考虑。具体步骤为：
[0089] STEP 1:任意找到一个多边形(假定为多边形A)，如果它没有被处理过，则选择此 多边形左上角的第一个像素点为种子，采用区域生长法构建它的连通体；否则，执行STEP 2,直到所有多边形都被处理完为止。
[0090] STEP 2:提取出组成此多边的所有点（包括边界点和结点），找出与这些点有关的 多所有多边形，并将这些多边形放入一个动态容器中。
[0091] STEP 3:检验这些多边形是否包含中多边形A中。因为这些多边形与多边形A只存 在包含或者邻接关系，所以本文首先检查这些多边形与多边形A是否具有邻接关系（多边形 的邻接关系在组成多边形时已经记录）。如果一个多边形与多边形A具有邻接关系，则它不 可能包含于多边形A中，所以将其排除。当排除完所有与多边形A邻接的多边形之后，剩余的 多边形全部包含于多边形A中。记录他们的包含关系，并清空动态容器，释放内存空间。
[0092] 第五，轮廓线的矢量化及光滑，首先，在轮廓线上选取初始的分段点;然后，根据斜 率的大小在轮廓线上依次选取其它分段点；最后，用一系列的点来表示物体的轮廓线。具体 为：
[0093] STEP 1:初始分段点的选取及折线拟合，对于每一条轮廓线，首先计算轮廓线上每 个点的左、右两侧的方向矢量，然后计算该轮廓线上每一个点的变化率，最后根据轮廓线上 点的变化率对其进行初步分段，并对每一段都采用最小二乘法进行折线段的拟合，得到一 组折线段。这样，一系列离散点就拟合成为一组折线段。
[0094] STEP 2:分段光滑
[0095] A自动判定折线段是否需要用曲线进一步拟合
[0096] 为了使曲线达到光滑连续的效果，在一条轮廓线的所有折线段中，一部分折线段 需要进一步拟合成曲线。将这些折线段分成两类，即一类是需要进一步进行曲线拟合的，另 一类则不需要采用曲线进行拟合。
[0097]对于由任意两个相邻的初始分段点Pa和Pb组成的一个折线段AB，首先假定采用直 线对其进行拟合，然后计算该折线段的拟合误差，根据拟合误差来判定它是否需要用曲线 进一步进行拟合。设直线PaPb的方程为：ax+by+c = 0,则折线段AB的拟合误差S定义如下：
>其中，cU为轮廓线上的点Pjlj折线段AB的垂直距离。令点坐标为( Xl，yi)，则
3若S越大，则说明拟合出直线AB的那些点越应该被拟合成曲线。
[0098] B直线或曲线拟合
[0099]当轮廓线某段(例如初始分段点Pa和Pb之间）为直线时，按上述方法进行分段，最后 所得分段点就为Pa和Pb两个端点，只需要用直线连接Pa和Pb即可。当轮廓线需要被拟合成曲 线时，采用Bezier曲线对其进行曲线拟合。
[0100] (2)提取非区域边界线信息
[0101] A断层识别
[0102] 通过图例信息得到断层线的RGB值，从而实现地质图断层的提取从标准彩色地质 图中识别断层线的方法，其具体步骤：
[0103] STEP 1:首先自动识别彩色地质图中的断层图例，并根据其图例得到断层线的RGB 值。
[0104] STEP 2:根据STEP 1中得到的断层线RGB值，然后提取出地质图中RGB值满足基其 一定范围内的的像素点。
[0105] STEP 3:对根据RGB颜色值提取出的断层线进行处理。分析断层线的特点，采用连 通体构建方法并根据面积特征和黑白比特征来去除伪断层线。
[0106] B 细化
[0107] 对识别出的断层线进行细化处理。细化算法的具体步骤：
[0108] STEP 1:判断是否全部像素点都被处理过，若全部像素点都被处理过，退出循环， 算法结束;否则执行STEP 2;
[0109] STEP 2:水平方向上的细化。
[0110]如果中心像素点P为黑像素时，进入下一步，否则就结束当前的循环；
[0111]判断中心像素点P点的左右两个邻点P3、P4,如果P3、P4中至少有一个等于1，则取得 中心像素点P点的8-邻域内的各个像素，按公式 25+P6*26+P7*27计算Index值，并查索引表Template。如果Template[Index] = 1，则P = 0，即P 点可以被删除；
[0112] 继续取下一个像素点。
[0113] STEP 3:垂直方向上的细化。
[0114] 如果中心像素点P为黑像素时，进入下一步，否则就结束当前的循环；
[0115] 判断中心像素点P点的上下两个邻点P^Ps，如果中至少有一个等于1，则取得 中心像素点P点的8-邻域内的各个像素，按公式 25+P6*26+P7*27计算Index值，并查索引表Template，如果Template[ Index] = 1，贝丨JP = 0，即P 点可以被删除；
[0116] 继续取下一个像素点。
[0117] C 剪枝
[0118]经过细化处理后图像中的线条已变成了单像素宽的骨架图像，保留下来的有些地 方存在多余的分支短线(又称为毛刺），为原始图像中不存的分支，对其进行剪枝处理。如考 查某几条线的交叉点像素（即8-邻域内黑象素的个数大于2)，依次跟踪每个方向上邻接的 黑像素点，当跟踪到达线端点或者另一个交叉点为止，统计其跟踪过的黑像素个数。如果在 某一个方向上黑像素的总个数小于某一阈值，则说明它是一个短"枝"，不构成一条线。因 此，必须将这个方向上搜索过的像素点的值由1变为〇(交叉点的像素值不变）。如果在某一 个方向上黑像素的总个数大于等于该阈值，则此方向上的像素值不变。
[0119] D断点的连接及矢量化
[0120] 两条断层线是否应该相连，主要取决于这两条曲线的走势及断点之间的距离。用d 来表示断层线两断点之间的距离，用断层线的断点处切线的夹角9来表示这两条曲线的走 势，用P来表示两条断层线连接的概率，则P与0、d之间的函数关系为:pzhlcose+h/d。其 中，是比例因子。取Ai = 0.6, = 0.4。断点连接完成后，轮廓矢量化及光滑方法对断 层线进行矢量化，得到其矢量数据。
[0121] 3提取点要素信息
[0122] 由于在标准彩色地质图中，图例与地质图具有完备性和一致性的原则，即图中所 有出现的矿点都会在地质图图例中用相同（包括大小、颜色等）的符号加以说明。因此，首先 提取出图例中的矿点符号，然后以提取的矿点符号为模板，通过序贯相似性检测算法进行 图像匹配，在地质图中找到相应矿点的位置，最后在相应位置生成矿点，形成矿点信息图 层。具体为：
[0123] 第一，提取矿点模板，分析矿点图例的特征，将其二值化后，通过面积特征和黑白 比特征将矿点图例全部自动识别出来，用于后续矿点图像匹配的模板图像。
[0124] 第二，识别矿点，采用改进的序贯相似性检测算法来识别地质图中的矿点，即在待 匹配图像的每个像素点处以一定的方式(如随机不重复的顺序)选择像元，并对模板和待匹 配图像在此像素点的差值进行累加，当累加值大于指定阈值时，则说明该像素点不是匹配 点，则停止对此像素点的计算，转而对下一个像素点的搜索，直到找到最佳匹配位置为止。 具体步骤为：
[0125] STEP 1:定义绝对误差函数：=以(M)-免;(U)-+ 印、/)|，其中，
S为待匹配图像，它的大小为MX N;T为模板图像，其大小为mXn;在待匹配图像中搜索模板图像时的移动窗口为SU，它的大 小也为m X n，i和j为模板覆盖区域(子图）的左上角像素点在待匹配图像S中的坐标，称其为 参考点。
[0126] STEP 2:确定自适应阈值Tk。自适应阈值Tk的构造方法为:在待搜索图像S(x，y)中， 计算它的每一个像素点（i，j)处的D( i，j)值，令阈值Tk的初始值为D( i，j);若在计算像素点 (1，」）处的值0(1，」）的过程中，0(1，」）的累加值超过了阈值1\时，则停止像素点（1，」）处的 搜索，转而计算下一个像素点；若计算完像素点（i，j)处的值D(i，j)且其值大于阈值T k时， 更新阈值Tk，即将D(i，j)的值赋给阈值Tk，并且将此像素点的位置坐标记录为匹配点。
[0127] STEP 3:在子图S1，j (k，1)中随机选取像素点，计算这个像素点与模板T中对应点的 误差值，然后将此差值与其它所有像素点的差值进行累加。当累加到r次且误差总和超过阈 值T k时，则累加停止并且记下累加次数r。序贯相似性检测算法的检测曲面定义为：
[0128] STEP 4:由于在I (i，j)值大的像素点（i，j)处，需要多次累加才能使总差误超过阈 值Tk，所以将此像素点作为匹配点。
[0129] 在实际的应用中，应该注意：
[0130] (1)由于矿点一般为圆形点状物，所以应使用圆形模板，采用序贯相似性检测算法 来搜索地质图中的矿点。因为一般情况下地质图中同一种矿不止一个，所以通常最佳匹配 位置也不止一个，因此匹配的相似程度大于一定阈值的位置都认为是矿点位置。
[0131] (2)采用粗一细结合的搜索方式寻找匹配点。在图像的(M-m+l)X(N-n+l)个像素 点中寻找所有匹配点时，由于地质图中矿点较少，大多数为非匹配点，所以先用粗网格的方 式搜索匹配点，并记录其粗匹配点位置;然后再以粗匹配点为矩心，大小为mXn( m，n为模板 图像的宽和高）的区域内进行逐点精匹配，寻找所有可能匹配点。这种策略不用计算所有匹 配点，可提高搜索效率。
[0132] 4信息输出方法
[0133] 将提取的面要素信息、线要素信息、点要素信息图层与现有GIS软件相结合，形成 矢量化的专题信息图层。
【主权项】
1. 一种地质图专题信息提取方法，其特征在于： 提取面要素信息，包括边界白色化处理、颜色分割、连通体构建、小色斑消除、边界处白 色区域修复； 提取线要素信息，包括提取区域边界线信息和提取非区域边界线信息； 提取点要素信息，包括生成矿点模板、匹配图像、定位矿点； 三种方法是独立运行互不影响的方法。2. 按照权利要求1所述的方法，其特征在于:其所述的边界白色化处理，具体为： STEP 1:对地质图进行边缘查找； STEP 2:对边界图像进行一次粗化处理； STEP 3:将原图像与边界图像相减，从而在原图中，将边界处全部进行白色化处理。3. 按照权利要求1所述的方法，其特征在于:其所述的颜色分割，具体为： STEP 1:相关参数的初始化。根据图例的中面要素区域的颜色来初始化聚类的类别数c 及聚类中心￥(<3) = {^^2丄，1};设定迭代停止的阈值^初始化迭代计数器为1 = 0和最大迭 代次数Imax = IOO;设加权指数m = 2。 STEP 2:当迭代次数1 < 100时，根据聚类中心V(1) = {V1，V2，L，v。}计算隶属度矩阵U(1)。 如果尹OWPiVfvp)像素点Pk到所有聚类中心的 <1(1 = 1，2，1^)的隶属度，得到隶属度矩阵1](1)。其中，(^(」=1，2丄，11;1^=1，2丄，11)为第1^ 个像素点Pk到第i个聚类中心vf ( i = 1，2，L，c)的距离。如果，=0 (即Λ = vf )，则>在1^13颜色空间中，定义任意两个像素点颜色的欧式距离为：STEP 3:对图像的颜色重新进行聚类划分。在图像的隶属度矩阵U(1)中，从每列中选择 隶属度最大的点作为对应点的归属域：重新计算每个类的聚 类中心V(1+1)。 STEP 4:如果I |V(1+1)-V(1)| |<e，则迭代停止并输出最终的划分矩阵U和聚类中心V，否 则令1 = 1+1，转向执行STEP 2。 STEP 5:去模糊化。根据模糊划分矩阵U提供的像素对各聚类中心的隶属度来确定该像 素点最后的归属。采用最大隶属转换方法，用Ck表示第k个像素点的归属类别，即Ck = argi {max{yik} I Vi5V^r t,4. 按照权利要求1所述的方法，其特征在于:其所述的连通体构建，具体为:对颜色分割 后的地质图图像，根据颜色可形成单色版的图层，然后进行扫描串的构建，接着根据扫描串 的邻接关系进而形成连通体，并计算每个连通体的面积和黑白比。5. 按照权利要求1所述的方法，其特征在于:其所述的小色斑消除，具体为： STEP 1:通过面积特征去除全部的孤立点状噪声和小面积片状噪声以及部分条带状噪 声； STEP 2:利用黑白比去除余下的条带状噪声； STEP 3:得到无噪声的面区域图像，完成颜色分割后面区域噪声的消除。6. 按照权利要求1所述的方法，其特征在于:其所述的边界处白色区域修复，具体为： STEP 1:扫描一次整张地质图的白色像素，标记出某个白色像素的8-邻域内有两个或 两个以上为前景像素(非白色）的像素点； STEP 2:将该像素点的颜色置换为4-邻域内颜色最多的那个颜色； STEP 3:重复STEP 1和STEP 2,层层修复，直到所有白色区域完全被修复为止。7. 按照权利要求1所述的方法，其特征在于:其所述的提取区域边界线信息，具体为： STEP 1:提取边界点和结点，逐行分析4个栅格单元交汇点的类型，快速地提取出图像 中所有的坐标点，并记录他们的行列坐标、类型以及连接信息； STEP 2:跟踪弧段，先跟踪开放弧段然后再跟踪封闭弧段； STEP 3:组建多边形，在建立弧段和结点的连接关系的基础上，采用"属性匹配"搜索多 边形的连接弧段组建多边形，在组成多边形的同时，解决多边形自交问题； STEP 4:建立拓扑关系； STEP 5:轮廓线的矢量化及光滑。8. 按照权利要求7所述的方法，其特征在于：其所述的区域轮廓线的矢量化及光滑方 法，具体为： STEP 1:初始分段点的选取及折线拟合，对于每一条轮廓线，首先计算轮廓线上每个点 的左、右两侧的方向矢量，然后计算该轮廓线上每一个点的变化率，最后根据轮廓线上点的 变化率对其进行初步分段，并对每一段都采用最小二乘法进行折线段的拟合，得到一组折 线段。 STEP 2:分段光滑，自动判定折线段是否需要用曲线进一步拟合。当判定轮廓线的某段 为直线时，直接用直线连接即可。当轮廓线需要被拟合成曲线时，采用Bezier曲线对其进行 曲线拟合。9. 按照权利要求1所述的方法，其特征在于:其所述的提取非区域边界线信息，具体为： STEP 1:断层识别，通过图例信息得到断层线的RGB值，实现断层线的粗提取，从地质图 中识别断层线； STEP 2:细化，对识别出的断层线进行细化处理； STEP 3:剪枝，经过细化处理，图像中的线条已变成了单像素宽的骨架图像，但有些地 方存在多余的分支短线，需要对其进行剪枝处理； STEP 4:断点的连接及矢量化。10. 按照权利要求1所述的方法，其特征在于:其所述的提取矿点模板，具体为：分析矿 点图例的特征，将其二值化后，通过面积特征和黑白比特征将矿点图例全部自动识别出来， 用于后续矿点图像匹配的模板图像。11. 按照权利要求1所述的方法，其特征在于:还包括地质图的预处理方法，对地质图图 像依次进行几何校正、彩色图像去噪、主图区域和图例区域定位。12.按照权利要求1所述的方法，其特征在于:还包括信息输出方法，将提取的面要素信 息、线要素信息、点要素信息与GIS软件相结合，形成矢量化的专题信息图层。
【文档编号】G06F17/30GK105931295SQ201610551178
【公开日】2016年9月7日
【申请日】2016年7月13日
【发明人】陈建平, 徐彬, 徐炀, 郑啸
【申请人】中国地质大学（北京）