CN105512651B - 一种基于能谱角的岩性分类方法 - Google Patents

Abstract

本发明提供了一种基于能谱角的岩性分类方法。该方法通过将每一采样点的能谱测量值的钾、铀、钍三个分量分别进行正规化处理，正规化处理后的采样点的能谱数据转换为三维坐标系中的一个能谱向量；计算每一采样点的能谱向量与工作区内所有预设的岩性聚类标准能谱向量之间的夹角（即能谱角）余弦，选取其中（）值最小且小于等于能谱角判断阈值t的初始聚类能谱向量所对应的岩性作为该采样点的推断岩性；若（）最小值大于能谱角判断阈值t，则将该采样点的岩性推断为新的岩性类别。通过本发明能够快速开展地层岩性分类，有效减弱地表浅覆盖层对下伏基岩能谱信息的屏蔽作用，进而提高岩性分类的可靠性和准确性。

Description

一种基于能谱角的岩性分类方法

技术领域

本发明涉及利用能谱信息辅助区域地质调查领域，具体地说是一种基于能谱角的岩性分类方法。

背景技术

在森林、草原等浅覆盖区，由于植被和土壤覆盖，基岩露头少，区域地质调查野外工作难度较大。隐伏的岩体、地层、构造、矿产(化)以及地表地质现象与深部地质构造变化的关系等信息无法直接观察研究；利用转石或局部零散露头获取信息量少，可信度低，难以提高填图的质量。若通过工程揭露基岩将大大增加工作费用、延长工作时间，大规模的施工也不利于森林、草原资源和自然景观的保护。

自然界中的各种岩石中均有不同含量的天然放射性核素。岩性及岩石的生成时代与放射性核素的含量具有明显的相关关系。据此，可以应用放射性核素测量数据进行岩性分类，辅助区域地质调查工作。

目前，国内外广泛应用的航空伽玛能谱测量可以提供广域的近地表放射性核素含量信息。利用能谱数据提取岩性信息的方法主要包括通过钾(K)、铀(U)、钍(Th)核素含量等值线平面图、核素比值(如Th/K、U/K)等值线平面图等划分岩性编图单元。这种方法对于放射性特征非常明显的岩性的识别效果较好，而对于放射性含量中等，差别不大的岩性分类效果欠佳，通常需要同时对比多种能谱信息，不易操作。

之后借鉴遥感图像处理方法，发展起来了能谱图像处理方法，即将放射性能谱信息中的钾(K)、铀(U)、钍(Th)核素含量转换为红、绿、蓝三色通道，进行假彩色合成，生成能谱图像。通过能谱图像的色彩信息进行岩性分类。此方法的岩性分类效果较之前的分类方法有了一定提高。该方法对于裸岩区岩性的可识别率较高，但是在地表有薄土层的浅覆盖区，由于土层厚薄不均，对下伏岩性放射性能谱信息的屏蔽影响亦有不同，因此此种情况下应用该方法进行岩性分类，其效果不佳，会产生较多的虚假岩性信息。

发明内容

本发明的目的就是提供一种基于能谱角的岩性分类方法，该方法适用于裸岩区和浅覆盖区，采用该方法能够从能谱测量数据中提取出不同类别岩性的空间分布信息。

本发明是这样实现的：一种基于能谱角的岩性分类方法，包括如下步骤：

a、在工作区开展面积性能谱测量或收集合适精度的能谱测量(航空能谱测量或地面能谱测量)数据，每一测量点的能谱数据均包括铀(U)、钍(Th)、钾(K)三个核素含量(或分量)。

b、对所有能谱测量数据进行统计和检查，将异常值进行标记和屏蔽，使其不参与区域性数据统计和以岩性分类为目的的后续操作。所述异常值是指其内至少有一个核素含量超出对应的预设阈值范围的测量点的能谱数据。

c、通过与地质图和遥感影像图空间配准后的对比，在工作区内圈定出水体、冲洪积物、人工建筑和大型工程分布区域，对该些区域的能谱测量数据进行标记和屏蔽，使其不参与区域性数据统计和以岩性分类为目的的后续操作。

d、为了得到空间分布均匀的采样数据，对工作区内剩余的符合参与分类要求的能谱测量值进行网格化处理。

e、对能谱网格数据中的U、Th和K含量进行正规化处理得到U′、Th′和K′。

首先，对参与分类的能谱网格数据进行统计，分别得到铀、钍、钾各分量的最大值(U_max、Th_max、K_max)和最小值(U_min、Th_min、K_min)。

对能谱网格数据中的U、Th、K含量进行正规化处理，三个核素含量的变换公式为：

上述三个公式中，U′_i、Th′_i和K′_i分别为正规化处理后第i个测量点的能谱数据中铀、钍、钾三个核素含量，U_i、Th_i和K_i分别为正规化处理前第i个测量点的能谱网格数据中铀、钍、钾三个核素含量，i＝1、2、……、n，n为工作区内参与分类的测量点数。变换后U′_i、Th′_i和K′_i的值均在0到1之间，以使三个核素含量值分布在同一量级。

f、生成能谱向量。

将U′、Th′、K′看作互相垂直的三个坐标轴，正规化处理后每个采样点的能谱值即成为三维空间中的一个能谱向量Q_i(U′_i,Th′_i,K′_i)。

将三维空间中的任意两个能谱向量之间的夹角，定义为能谱角θ。采样点间的能谱角θ与采样点能谱向量的模的大小无关，仅与不同岩性之间物质成分的差异有关。

从而利用能谱角的大小可以快速甄别出二者岩性的相似性，并在一定程度上减弱地表浅层覆土对下伏岩性的伽马射线能量的屏蔽影响。

g、参考已有地质资料和/或进行实地能谱测量、采样及岩矿鉴定，为工作区内预设初始地层岩性类别标志样本。

h、将初始地层岩性类别标志样本内的所有正规化后的U′、Th′、K′进行直方图分析，找出峰值，转化为三维直角坐标系中的三维坐标值，作为地层岩性的初始聚类标志向量。

i、计算工作区内的所有参与分类的测量点的能谱向量与各预设初始聚类标志向量之间能谱角的余弦值cosθ。

设工作区内的某采样点的能谱向量为Q_i(U′_i,Th′_i,K′_i)，某预设初始聚类标志向量为Q_j(U′_j,Th′_j,K′_j)，则Q_i(U′_i,Th′_i,K′_i)和Q_j(U′_j,Th′_j,K′_j)的能谱角余弦cosθ_ij的计算公式为：

其中，i＝1、2、……、n，n为工作区内所有参与分类的测量点的数量；j＝1、2、……、m，m为预设初始聚类标志向量数。

设能谱角函数(或称能谱角判别函数)为S，且S＝1-cosθ，则工作区内任一采样点(i)对应某一类岩性(j)的能谱角函数为S_ij。

S_ij是一个介于0到1之间的值，S_ij越小表示采样点i与岩性j的相似性越大，S_ij越大表示采样点i与岩性j的相似性越小。

设定适宜于工作区内岩性差异的能谱角判断阈值t，0＜t＜1。

j、选取每一采样点与所有预设岩性标准聚类向量的能谱角判别函数S中的最小值S_k。

若S_ik≤t，则采样点i的岩性与第k类岩性最为相似，将其岩性推断为第k类岩性。

若S_ik＞t，则采样点i的岩性与所有预设的岩性的相似度均较低，将其岩性推断为新的岩性。

k、在第一轮岩性分类结束后，对统计推断为新类别的正规化处理后的各采样点进行核素含量直方图分析，如为近似正态分布的单峰，则设立一个新的岩性聚类标志样本，取其铀、钍、钾三个核素含量的峰值，构成新的岩性聚类标志向量；如为双峰或多峰分布，则添加两个或多个新的岩性聚类标志样本，并分别取其峰值，构成两个或多个新的岩性聚类标志向量。

将形成的新的岩性聚类标志向量加入初始的岩性聚类标志向量集中，在补充后的岩性聚类标志向量集的基础上进行第二轮岩性分类。如此反复，直至再无新的岩性产生，分类结果基本合理。

l、在分类完成后，经实地查证、采样和岩矿鉴定，确定各类岩性名称。

本发明通过将每一个测量点(或称采样点)的能谱测量值的钾(K)、铀(U)、钍(Th)三个分量分别进行正规化处理后，作为三维坐标系中的一组坐标值，将该采样点转换为三维坐标系中的一个能谱向量。计算每一个采样点的能谱向量与工作区内所有预设的岩性聚类标准能谱向量之间的夹角(即能谱角)，选取其中最小的能谱角，二者的能谱核素含量相似度最高。根据能谱角求能谱角函数S，最小的能谱角对应的能谱角函数S最小。判断最小的能谱角函数S是否大于预设能谱角判断阈值t，如果是，则将该采样点的岩性推断为新的岩性，如果否，则将此聚类标准能谱向量所对应的岩性推断为该采样点的岩性。重复上述步骤，直至工作区内所有采样点的岩性均推断完成，完成第一轮岩性的分类。在第一轮岩性分类结束后，将所有推断出新岩性的聚类标志能谱向量加入预设的岩性聚类能谱向量集中，在补充后的岩性聚类向量集的基础上进行第二轮岩性分类；如此反复，直至再无新类别岩性的产生。

本发明通过采样点的能谱向量与工作区内所有预设的岩性聚类标准能谱向量之间的夹角(即能谱角)来判断采样点的岩性与聚类标准能谱向量所对应的岩性之间的相似度，进而判定采样点的岩性为新类别的岩性或聚类标准能谱向量所对应的岩性。由于能谱角与能谱向量的模无关，仅与不同岩性之间物质成分的差异有关，从而可以减弱地表浅层覆土对下伏岩性的伽马射线能量的屏蔽影响。对于地表覆土厚度不均造成的虚假岩性信息具有较好的甄别效果，岩性分类准确率较现有的能谱含量、比值参数等信息的解译效果有了较大的提升，并可对侵入岩、火山岩生成期次提供新的佐证，辅助地质人员更深入全面地认识了解区域地质过程。

本发明具有以下优点：(1)简单易操作，分类指标物理意义明确；(2)能够降低地表覆土对下伏岩石(地层)能谱值的影响；(3)可以对照地质图、遥感影像图和野外验证情况，通过修改预设岩性类别、增减类别标志样本(即聚类标准能谱向量)、调节能谱角判断阈值等操作，多次修正分类结果，直至达到满意。

附图说明

图1是本发明的方法流程图。

具体实施方式

下面结合具体实施例，对本发明中的技术方案进行清楚、完整地描述。

如图1所示，本发明所提供的基于能谱角的岩性分类方法，包括如下步骤：

(1)、在黑龙江省小兴安岭东南段工作区，首先对区内的航空能谱测量数据进行必要的整理、统计。除航空能谱测量数据外，还收集了该工作区的地质图以及遥感影像图。

对能谱测量数据进行统计、检查，找出异常值并进行标记和屏蔽，使其不参与后续的分类处理。此处的异常值是指其内至少有一个核素含量超出对应的预设阈值范围的测量点的能谱数据。本实施例中，在测量点的能谱数据中，核素含量对应的预设阈值范围为：核素含量平均值±3倍均方差。K、U、Th三个核素，只要有一个核素含量超出其平均值加减3倍均方差的范围，则该采样点即被视为异常值而排除。异常值(区)所处区域在实地查证中通过能谱测量、采样及岩矿鉴定确定岩性。

去除异常值后，使工作区与空间配准后的地质图和遥感影像图相对照，将工作区内的水体、冲洪积物、人工建筑(如居民地等)和大型工程区域(如公路、地铁、机场等)进行圈定和屏蔽，使这些区域内的采样点不参与后续的分类处理。

(2)、对工作区内剩余的所有参与分类的采样点的能谱数据进行网格化处理。

设定合适的网格大小，通常为测线间距d的1/4，对工作区内剩余的能谱测量数据进行网格化处理，使其转换为在空间上均匀分布的采样点。网格化处理后的能谱数据可以称为能谱网格数据。

(3)、对能谱网格数据中的U、Th、K含量进行正规化处理以得到U′、Th′、K′；

首先统计已屏蔽了异常值、水体、冲洪积物、人工建筑和大型工程区域后的所有参与分类的能谱网格数据中铀(U)、钍(Th)、钾(K)三个核素含量各自对应的最大值U_max、Th_max、K_max，以及三个核素含量各自对应的最小值U_min、Th_min、K_min，根据这六个最值并依据下面三个公式对参与分类的采样点的能谱网格数据中的U、Th、K含量进行正规化处理。

上面三个公式中，i＝1、2、……、n，n(为正整数)为所有参与分类的采样点的个数。U_i、Th_i和K_i分别为网格化处理后、正规化处理前第i个测量点的能谱网格数据中铀、钍、钾三个核素含量，U′_i、Th′_i和K′_i分别为正规化处理后第i个测量点的能谱数据中铀、钍、钾三个核素含量。将U′、Th′、K′看作互相垂直的三个坐标轴，正规化处理后的每一测量点的能谱数据构成三维空间的一个能谱向量Q_i(U′_i,Th′_i,K′_i)。

(4)、确定初始岩性编图单元。

对照已有区域地质图的地层和侵入岩体进行初步分区，将工作区划分为15类编图单元，生成带有岩性(地层)名称、代码等属性的编图单元区文件。

将岩性编图单元分区与正规化处理后的能谱数据作空间相交运算，使得每一采样点的能谱数据都带有所属的编图单元属性信息。

(5)、生成初始岩性聚类向量集。

对15类编图单元各分区内正规化处理后的采样点的能谱数据进行统计分析，具体是：对各类编图单元内所有正规化处理后的采样点的能谱数据进行铀、钍、钾三个核素含量的直方图分析，若三个直方图均为近似正态的单峰分布，则在三个直方图中分别选取铀、钍、钾三个核素含量的峰值，构成该类地层或侵入岩体单元的初始聚类标志样本核素含量值；若三个直方图中至少有一个直方图为双峰或多峰分布，则通过组合在三个直方图中分别选取铀、钍、钾三个核素含量的峰值，构成两个或多个初始聚类标志样本核素含量值。

经统计分析，将工作区的初始岩性聚类标志样本扩充为38个，这38个初始岩性聚类标志样本所对应的能谱向量即构成了初始岩性聚类向量集(或称预设岩性聚类向量集)。设定38个初始聚类标志向量组T[T1，T2，……，T38]，通过岩性能谱资料分析和必要的野外查证，建立初始岩性类别对照表，如T1—早奥陶世黑云母二长花岗岩，T2—晚奥陶世中细粒花岗闪长岩等。

(6)、计算每一采样点能谱向量与初始岩性聚类向量集中每一向量的夹角(能谱角)余弦。

计算工作区内所有参与岩性分类的采样点的能谱向量与38个初始聚类标志向量之间的能谱角θ的余弦值cosθ，得到每个采样点Q_i对应于38个初始类别的能谱角余弦组cosθ_ij[cosθ_i1，cosθ_i2，……，cosθ_i38]。

(7)、选取能谱角判别函数S＝1-cosθ的最小值S_k。

设能谱角判别函数为S，任一采样点Q_i(U′_i,Th′_i,K′_i)与任一岩性类别Q_j(U′_j,Th′_j,K′_j)之间的判别函数S_ij由下式计算而得：

上式中，i＝1、2、……、n，n为工作区内参与岩性分类的采样点数；j＝1、2、……、m，m为初始岩性聚类标志向量数，本实施例中m＝38。

本实施例中，对于任一采样点Q_i，其与38个初始岩性聚类标志向量所形成的38个能谱角判别函数为S_i[S_i1，S_i2，……，S_i38]，选取这38个能谱角判别函数中的最小值S_ik＝min[S_i]。

(8)、判断能谱角判别函数的最小值S_k是否大于预设能谱角判断阈值t。

本实施例中设定能谱角判断阈值t＝0.2。

如果S_ik≤0.2，则说明采样点Q_i的岩性与第k个岩性类别相似度最高，将所对应的岩性作为该采样点的推断岩性。

如果S_ik＞0.2，则说明采样点Q_i的岩性与所有预设的岩性的相似度均较低，将该采样点的岩性类别定为新的岩性类别。

按照步骤(6)、(7)、(8)，对所有采样点的岩性均判断完毕后，第一轮岩性分类结束。

(9)、在第一轮岩性分类结束后，对照地质图和遥感影像图对分类结果进行分析，对于过度分类或欠分类的地区可适当调整不同岩性地层的能谱角判断阈值t，然后再进行下一轮的岩性分类。

(10)、在第二轮岩性分类之前，首先将第一轮岩性分类后所有能谱角判别函数最小值均大于能谱角判断阈值的采样点，根据区域地质和地貌特征，将其归为新的岩性类别。

(11)、统计分析所有推断为新类别的正规化处理后的采样点的铀、钍、钾三个核素含量直方图，若三个直方图均为近似正态的单峰分布，则设立一个新的岩性聚类标志样本，在三个直方图中分别取铀、钍、钾三个核素含量的峰值，构成新的岩性聚类向量；若三个直方图中至少有一个直方图为双峰或多峰分布，则添加两个或多个新的岩性聚类标志样本，并通过组合方式在三个直方图中分别选取铀、钍、钾三个核素含量的峰值，构成两个或多个新的岩性聚类向量。

(12)、将这些所形成的新的岩性聚类向量补入初始的岩性聚类标志样本集中，在补充后的岩性聚类标志样本集的基础上进行第二轮岩性分类。

(13)、第二轮岩性分类结束后，可以再对照地质图和遥感影像图对分类结果进行分析，对于过度分类或欠分类的地区可适当调整不同岩性地层的能谱角判断阈值t，然后再进行下一轮的岩性分类。

(14)、如此反复，直至再无新的岩性类别产生，分类结果基本合理。

将基于能谱角的岩性分类解译推断结果，编制成岩性能谱解译图，对于其中与原有地质图有较大差别的区域，作为野外查证的目标进行调查。

最终，通过基于能谱角的计算机自动岩性分类、野外查证、采样和岩矿鉴定，将所选实施例的工作区划分为45类岩性填图单元，经过野外实地查证，分类准确率较采用现有的能谱单核素含量等值图、能谱核素比值等值图和能谱假彩色图像等方法的岩性分类准确率均有较大幅度提高。

通过本发明中的方法能够快速开展地层岩性分类，有效减弱地表浅覆盖层对下伏基岩能谱信息的屏蔽作用。本发明作为一种定量计算方法，减少了人为操作的不确定因素，提高了岩性分类的可靠性和系统性。

上述实施例仅为清楚地描述了本发明的具体实施过程，并非对实施过程的限定。对于所属领域的普通技术人员来说，在上述说明的基础上，还可以做出其他不同形式的变化或变动，这里无需也无法对所有实施方式予以穷举。由此所引申的显而易见的变化或变动仍处于本发明的保护范围之内。

Claims (5)

1.一种基于能谱角的岩性分类方法，其特征是，包括如下步骤：

a、获取工作区内的能谱测量数据，每一测量点的能谱数据均包括铀、钍、钾三个核素含量；

b、在所有能谱测量数据中，去除异常值；所述异常值是指其内至少有一个核素含量超出对应的预设阈值范围的测量点的能谱数据；

c、在工作区内圈定水体、冲洪积物、人工建筑和工程分布区域，并将这些区域内的测量点的能谱数据去除；

d、对工作区内剩余的测量点的能谱数据进行网格化处理，使其转换为在空间上均匀分布的能谱数据；网格化处理后的能谱数据称为能谱网格数据；

e、对网格化处理后的能谱网格数据中的U、Th、K含量进行正规化处理以得到U′、Th′、K′；

所述正规化处理是指使测量点的能谱网格数据中的三个核素含量分别按如下公式进行变换：

上述三个公式中，U′_i、Th′_i和K′_i分别为正规化处理后第i个测量点的能谱数据中铀、钍、钾三个核素含量，U_i、Th_i和K_i分别为正规化处理前第i个测量点的能谱网格数据中铀、钍、钾三个核素含量，U_max、Th_max和K_max分别为所有能谱网格数据中铀、钍、钾三个核素含量的最大值，U_min、Th_min和K_min分别为所有能谱网格数据中铀、钍、钾三个核素含量的最小值；i＝1、2、……、n，n为步骤d中工作区内剩余的测量点的数量；

f、将U′、Th′、K′看作互相垂直的三个坐标轴，正规化处理后的每一测量点的能谱数据构成三维空间的一个能谱向量Q_i(U′_i,Th′_i,K′_i)；

g、使第i个测量点的能谱向量Q_i(U′_i,Th′_i,K′_i)与预设岩性聚类向量集中的每个向量Q_j(U′_j,Th′_j,K′_j)分别求能谱角余弦cosθ_ij，并根据能谱角余弦cosθ_ij求能谱角函数S_ij；能谱角余弦cosθ_ij和能谱角函数S_ij的计算公式分别如下：

其中，i＝1、2、……、n，n为步骤d中工作区内剩余的测量点的数量；j＝1、2、……、m，m为预设岩性聚类向量集中向量的个数；

h、对于第i个测量点的能谱向量Q_i(U′_i,Th′_i,K′_i)，比较其与预设岩性聚类向量集中m个向量的能谱角函数S_i1、S_i2、……、S_im之间的大小，并从中选出最小值，记录该最小值为S_ik；

i、判断步骤h中最小值S_ik是否大于预设能谱角判断阈值t，如果是，则将第i个测量点的岩性推断为新类别的岩性；如果否，则将第i个测量点的岩性推断为第k类岩性；

j、重复步骤g～i，直至n个测量点的岩性均推断完毕，此时完成第一轮岩性的分类；

k、在第一轮岩性分类结束后，将所有推断为新类别岩性的正规化处理后的能谱数据进行铀、钍、钾三个核素含量的直方图分析，若三个直方图均为单峰分布，则在三个直方图中分别选取铀、钍、钾三个核素含量的峰值，构成一个新的岩性聚类能谱向量；若三个直方图中至少有一个直方图为双峰或多峰分布，则通过组合在三个直方图中分别选取铀、钍、钾三个核素含量的峰值，构成两个或多个新的岩性聚类能谱向量；

l、将步骤k中形成的新的岩性聚类能谱向量补入预设岩性聚类向量集中，重复步骤g～j，在补充后的岩性聚类向量集的基础上进行第二轮岩性分类；如此反复，直至再无新类别岩性的产生。

2.根据权利要求1所述的基于能谱角的岩性分类方法，其特征是，步骤g中预设岩性聚类向量集是通过如下方式形成的：

①、对照工作区已有的地质图，将工作区划分为若干类编图单元；

②、对各类编图单元内所有正规化处理后的能谱数据进行铀、钍、钾三个核素含量的直方图分析，若三个直方图均为单峰分布，则在三个直方图中分别选取铀、钍、钾三个核素含量的峰值，构成一个初始聚类标志样本向量；若三个直方图中至少有一个直方图为双峰或多峰分布，则通过组合在三个直方图中分别选取铀、钍、钾三个核素含量的峰值，构成两个或多个初始聚类标志样本向量；

③、将步骤②中所形成的所有初始聚类标志样本向量集合起来即形成预设岩性聚类向量集。

3.根据权利要求1所述的基于能谱角的岩性分类方法，其特征是，在每一轮岩性分类结束后，根据工作区已有的地质图和遥感影像图对分类结果进行分析，并根据分析结果调整预设能谱角判断阈值t。

4.根据权利要求1所述的基于能谱角的岩性分类方法，其特征是，步骤b中异常值是指其内至少有一个核素含量超出对应的预设阈值范围的测量点的能谱数据；在测量点的能谱数据中，核素含量对应的预设阈值范围为：核素含量平均值±3倍均方差。

5.根据权利要求1所述的基于能谱角的岩性分类方法，其特征是，在分类完成后，经实地查证、采样和岩矿鉴定，确定各类岩性名称。