CN101539998B - 蚀变遥感异常提取方法和系统 - Google Patents

Abstract

本发明提出了一种蚀变遥感异常提取方法和装置，用于从输入的遥感数据中提取表征与矿化相关的近矿蚀变岩石的信息，所述遥感数据包括多个波段的某地理区域内的像元值。所述方法包括下述步骤：对输入的遥感数据进行前处理；从经所述前处理后得到的遥感数据所表征的遥感图像中切割出包含已知矿床类型的第一图像区域；对该第一图像区域，利用预先确定的阈值从该第一图像区域中选择其中异常符合该阈值的第二图像区域；以及针对该第二图像区域，基于该已知矿床类型的参考谱，确定第二图像区域和该已知矿床类型之间的相似程度。

Description

蚀变遥感异常提取方法和系统

本发明受国家科技支撑计划重点项目(编号：2006BAB07B07-02)与国家重大基础研究计划项目(编号：2007CB411406)，地质调查工作项目(编号：1212010818018)与中央级公益性科研院所基本科研业务费专项资金资助项目(编号：K2007-2-7)的资助。

技术领域

本发明一般地涉及利用信息技术辅助矿床勘查，具体地涉及蚀变遥感异常提取的方法和装置。

背景技术

近矿围岩蚀变是成矿物质逐步富集成矿过程中留下的印迹。最常见的蚀变为：硅化、绢云母化、绿泥石化、云英岩化、矽卡岩化、白云岩化、重晶石化及锰铁碳酸盐化。地质学家断言绝大多数内生矿床都伴随有其围岩的交代蚀变现象。因而可以利用近矿蚀变围岩指导找矿。

所谓蚀变遥感异常信息是利用数学手段从遥感数据量化提取出来的、用以表征有可能与矿化最相关的近矿蚀变岩石的指示信息。蚀变遥感异常信息作为找矿标志已被广泛应用。

已经发展了很多利用遥感数据进行矿化信息提取的技术。赵元洪等在“波段比值的主成分复合在热液蚀变信息提取中的应用”.国土资源遥感，1991(3)中提出了波段比值的主成分复合法。1994年，何国金，胡德永，陈志军等在“从TM图像中直接提取金矿化信息”.遥感技术与应用，1995，10(3)：51-54中提出了“微量信息处理”方法。张远飞，吴健生在“基于遥感图像提取矿化蚀变信息”有色金属矿产与勘查，1999，8(6)：604-606中利用“多元数据分析+比值+主成份变换+掩膜+分类(分割)”的方法在新疆、内蒙古及江西、云南成功的提取了金矿化蚀变信息。刘成，王丹丽，李笑梅在“用混合像元线性模型提取中等植被覆盖区的粘土蚀变信息”[J].地质找矿论，2003，18(2)：131-137)中利用混合像元线性分解模型提取了卧龙泉地区粘土蚀变信息。

但是，目前所利用的这些技术提取出来的异常往往包含各种类型矿床的蚀变信息，而且包含有大量的与矿化关系不大的“面积型”异常。所谓“面积型”异常是指在一定切割水平上，所提取出来的异常明显成片，异常范围相对较大，大约10公里这个数量级。这类异常往往与矿化关系不大，通常由区域性岩石蚀变引起，例如在新疆东天山翠岭地区火山岩的绿帘石化。如果存在大量的遥感蚀变异常与矿化没有关系，将使得提取出来的遥感蚀变异常缺乏找矿的针对性，对地质矿产勘查来说，增加了分析与判别有效找矿信息的难度。

这里所说的分析与判别有效找矿信息是非常重要的。将遥感技术应用于找矿勘查，一般分为三个步骤：数据处理、数据分析和有效找矿信息判别。对应用蚀变遥感异常找矿来说，首先需要提取遥感数据中包含的蚀变遥感异常，然后分析这些异常的性质，判别其与矿化的关系。已有的蚀变遥感异常提取方法提取出来的异常信息是一个包含各种蚀变的信息，比如，利用主成分分析技术提取出来的异常是某个特征组合优先排序的结果，其中包含了各种矿床的混杂的信息。从如此混杂的信息中提取有效找矿信息是很困难和耗时的。

发明内容

针对现有技术的上述问题，提出了本发明。

根据本发明的一个方面，提出了一种蚀变遥感异常提取方法，用于从输入的遥感数据提取表征与矿化相关的近矿蚀变岩石的信息，所述遥感数据包括多个波段的地理区域内的像元值，包括下述步骤：对输入的遥感数据进行前处理；从经所述前处理后得到的遥感数据所表征的遥感图像中切割出包含已知矿床类型的第一图像区域；对该第一图像区域，利用预先确定的阈值从该第一图像区域中选择其中异常符合该阈值的第二图像区域；以及针对该第二图像区域，基于该已知矿床类型的参考谱，确定第二图像区域和该已知矿床类型之间的相似程度，其中，所述前处理步骤包括：进行边界信息剔除，由此通过去除边界信息，使得各个波段的边界处信息重合，所谓边界信息是指，在获取遥感数据时，各个波段在边界处的信息不重合区域的信息，判断有无云、水体、阴影区、白泥地、冰雪、植被、湿地、干河道、冲积扇九类常见干扰异常，如有则进行去除，去除方法为比值法、高端或低端切割法、Q值法、光谱角法中的一个或其组合；统计遥感区域内的像元值，并以直方图表示，判断该直方图是否符合正态分布，如果不满足正态分布，则把该遥感区域划分为多个次级区域，使得每个次级区域内的每个波段直方图符合正态分布；将符合正态分布的区域内的像元值统一转换到[0，255]的范围；切割掉正态分布区域中的边缘部分，对输入的数据进行主成分分析，并将输入的数据转换到通过主成分分析所确定的特征向量所表征的空间坐标系中，所述选择第二图像区域的步骤包括：如果切割第一图像区域的步骤中切割出包含已知矿床类型的第一图像区域，则利用下述公式确定第二图像区域：

$\left\{\begin{array}{c}I=A>{\mathrm{\&epsiv;}}_{\min }\\ I=A<{\mathrm{\&epsiv;}}_{\max }\end{array}\right.$

其中，I为选取的第二图像区域，A为连通的异常面积，ε_min与ε_max分别是预定的异常下限阈值和异常上限阈值。

根据本发明的另一方面，提取出一种蚀变遥感异常提取装置，用于从输入的遥感数据提取表征与矿化相关的近矿蚀变岩石的信息，所述遥感数据包括多个波段的地理区域内的像元值，包括：对输入的遥感数据进行前处理的部件；从经所述前处理后得到的遥感数据所表征的遥感图像中切割出包含已知矿床类型的第一图像区域的部件；对该第一图像区域，利用预先确定的阈值从该第一图像区域中选择其中异常符合该阈值的第二图像区域的部件；以及针对该第二图像区域，基于该已知矿床类型的参考谱，确定第二图像区域和该已知矿床类型之间的相似程度的部件，其中，所述前处理部件进行下述操作：边界信息剔除，由此通过去除边界信息，使得各个波段的边界处信息重合，所谓边界信息是指，在获取遥感数据时，各个波段在边界处的信息不重合区域的信息，判断有无云、水体、阴影区、白泥地、冰雪、植被、湿地、干河道、冲积扇九类常见干扰异常，如有则进行去除，去除方法为比值法、高端或低端切割法、Q值法、光谱角法中的一个或其组合；统计遥感区域内的像元值，并以直方图表示，判断该直方图是否符合正态分布，如果不满足正态分布，则把该遥感区域划分为多个次级区域，使得每个次级区域内的每个波段直方图符合正态分布；将符合正态分布的区域内的像元值统一转换到[0，255]的范围；以及切割掉正态分布区域中的边缘部分，对输入的数据进行主成分分析，并将输入的数据转换到通过主成分分析所确定的特征向量所表征的空间坐标系中，其中，如果感兴趣区域初步切割部件切割出包含已知矿床类型的第一图像区域，则感兴趣区域细化选取部件利用下述公式确定第二图像区域：

$\left\{\begin{array}{c}I=A>{\mathrm{\&epsiv;}}_{\min }\\ I=A<{\mathrm{\&epsiv;}}_{\max }\end{array}\right.$

其中，I为选取的第二图像区域，A为连通的异常面积，ε_min与ε_max分别是预定的异常下限阈值和异常上限阈值。

本发明基于某一特定类型矿床的蚀变岩石具有独特的波谱特征，由此提取该类型矿床的蚀变遥感异常，并将其结果有针对性地应用于找寻该类型矿床。

利用本发明的技术方案，能够解决已有蚀变遥感异常提取方法所获取的异常分布范围广，大量的与找矿关系不大的“面积型”异常出现，找矿没有针对性等问题。本发明采用矿床分类型蚀变遥感异常提取技术获取的异常分布范围小，减少了与找矿关系不大的“面积型”异常，对该类矿床的勘查针对性强，有利于更快更准的优选找矿靶区，对于矿产勘查工作可以起到节省时间、节约人力物力，事半功倍的作用，是科学技术进步促进生产发展的新技术。

该技术的有效性在实践中得到了验证。申请人利用这种方法发现了多个矿(化)点，为当地的找矿做出了显著贡献。

附图说明

通过结合附图对本发明各个实施例的详细描述，本领域的技术人员可以更好地理解本发明的上述和其它方面、特点和优点，其中在所有附图中使用相同或相似的附图标记来表示相同或者相似的部件，在附图中：

图1示出了根据本发明一个实施例的蚀变遥感异常提取方法的流程图；

图2示出了根据本发明一个实施例的遥感数据前处理的流程图；

图3示意性地示出了根据本发明一个实施例的边界信息；

图4示意性地示出了边界信息去除前后的对比示图；

图5示出了根据本发明一个实施例的蚀变遥感异常提取装置的框图。

具体实施方式

下面将结合附图详细描述本发明的具体实施例。如果考虑到对某些相关现有技术的详细描述可能会混淆本发明的要点，则不会在这里提供其详细描述。

本发明的矿床分类型异常提取是建立在分析测试某种类型矿床岩石波谱特征的基础上的一种有针对性的蚀变遥感异常提取技术。这种技术的基础是典型矿床中大量蚀变岩石的波谱测试分析，波谱曲线特征统计，规律总结。例如，通过对新疆东天山地区黄山东镁铁质铜镍矿、沙泉子夕卡岩型铅锌矿和天目剪切带型金矿三种不同类型矿床的波谱统计结果，可以看出，三种不同类型矿床的波谱特征差别显著。正是基于上述发现，本发明所采用的分类技术不同于遥感处理常见的监督分类技术和非监督分类技术，是一种与典型矿床的空间分布及蚀变岩石的波谱特征相关的特征异常区门限化监督分类技术。

下面将参照附图，对于本发明的示例性实施例进行描述。

图1示出了根据本发明一个实施例的蚀变遥感异常提取方法100的流程图。

如图1所示，根据本发明一个实施例的蚀变遥感异常提取方法100包括：用于从输入的遥感数据提取表征与矿化相关的近矿蚀变岩石的信息，所述遥感数据包括多个波段的某地理区域内的像元值，包括下述步骤：在步骤S110，对输入的遥感数据进行前处理；在步骤S120，从经所述前处理后得到的遥感数据所表征的遥感图像中切割出包含已知矿床类型的第一图像区域；在步骤S130，对该第一图像区域，利用预先确定的阈值从该第一图像区域中选择其中异常符合该阈值的第二图像区域；以及在步骤S140，针对该第二图像区域，基于该已知矿床类型的参考谱，确定第二图像区域和该已知矿床类型之间的相似程度。

首先，需要说明的是，这里所说的遥感数据，是一种基于地面光谱特征的中低分辨率的多光谱卫星遥感数据，所获得的数据格式一般是数字的，即计算机储存的位图格式，一般从卫星地面站或其他代理机构购买，数据反映了每个波段辐射特征的灰度值。

在本发明的一个实施例中，数据选择自某个日期etm或aster数据中某个地区的数据，例如选择如下期间的数据：植被发育弱、冰雪覆盖较少的季节，同时要求数据获取期间，该区域的天空云量较少，比如选择东天山成矿带的某景etm或aster数据。etm数据与aster数据是两种多光谱遥感数据，地质学家或矿务学家常常将其用于地质找矿，这两种数据是公开的，任何人都可以得到。Etm指Enhanced Thematic Mapper Plus，增强性专题制图仪，tm指Thematic Mapper Plus专题制图仪，是安放在美国陆地卫星上的传感器，用于接受地面的数据。Aster指Advanced Spaceborne Thermal Emission and Reflection Radiometer，先进星载热发射和反射辐射仪。是安放在美国Terra卫星上的日本设计的传感器。选择数据以后我们要对数据进行影像合成，一般对于etm数据选择743波段组合合成影像数据，对于aster数据一般选择631波段组合合成影像数据。作为示例，所谓影像数据可以采用RGB格式，例如R放置6波段(或7波段)，G放置3波段(或4波段)，B放置1波段(或3波段)。影像数据用于显示遥感图像的影像。

下面，参照附图对于图1所示流程图中的各步骤进行详细说明。

一、遥感数据前处理

图2示出了根据本发明一个实施例的遥感数据前处理S110的流程图。

如图2所示，根据本发明一个实施例，在步骤S111，进行边界信息剔除，所谓边界信息是指，在获取遥感数据时，各个波段在边界处的信息不重合区域的信息，由此通过去除边界信息，使得各个波段的边界处信息重合。

在步骤S111后，前进到步骤S112。

在步骤S112，去除常见干扰异常。具体地，对于经过边界信息去除处理的遥感数据，判断有无云、水体、阴影区、白泥地、冰雪、植被、湿地、干河道、冲积扇九类常见干扰异常，如有则进行去除，去除方法为比值法、高端或低端切割法、Q值法、光谱角法中的一个或其组合。然后，前进到步骤S113。

在步骤S113，对于去除了统计某个遥感区域内的像元值，并以直方图表示，判断该直方图是否符合正态分布，如果不满足正态分布，则把该遥感区域划分为多个次级区域，使得每个次级区域内的每个波段直方图符合正态分布，这是因为正态分布能够使得分析后蚀变信息量更大，并且提取结果更准确。然后，前进到步骤S114。

在步骤S114，将符合各个正态分布的区域内的像元值统一转换到[0，255]的范围。然后，前进到步骤S115。

在步骤S115，对经步骤S115处理得到的数据进行主成分变换。

在步骤S116，切割掉正态分布区域的边缘部分。

下面对剔除边界信息步骤S111、去除常见干绕S112、进一步进行详细说明。

1剔除边界信息

图3示意性地示出了根据本发明一个实施例的边界信息。

如图3所示，三个波段红色R、绿色G、蓝色B的信息在平面(代表地球表面)的坐标(X，Y)是不重合的。比如，蓝色B波段在左侧所示部分没有信息，而红色R波段在右侧所示部分不存在信息。边界信息剔除指的是，提取出三个波段均存在信息的区域，而去除掉存在信息不重叠的边界区域。例如，在该实施例中，剔除边界信息后的区域是图中中间的灰色区域。边界信息的提取使得每一个波段对应部分都含有信息。

2、去除常见干扰

通常，遥感数据中，要避免即过滤掉云、水体、阴影区、白泥地、冰雪、植被、湿地、干河道、冲积扇等九类常见干扰异常。对于干扰异常的检测可以采用目估法，一般干扰地物能够在TM/ETM的743彩色合成图像或ASTER的631彩色合成图像有明显的特征，比如云的白色等，去除方法选用比值法、高端或低端切割法、Q值法、光谱角法等。

下面分别加以说明。

(1)高端或低端切割法

主要是利用干扰地物在遥感图像上某个波段有特征的高反射或强吸收，即某波段干扰地物有高值或低值。比如水体在TM/ETM的第7波段有低值，采用低端切割方法处理，而云在TM/ETM的第1波段有高值，采用高端切割方法处理，白泥地在TM/ETM的第3波段有高值，采用高端切割方法处理等。公式如下：

其中，i＝0，...，n，n指所使用的遥感图像波段总数，x_i和y_i指分别指i波段去除干扰信息前和后的波段值，b∈[1，...，n]，C_b是阈值，x_b是原始b波段对应的值。采用这个公式的目的在于：给定一个约束的条件，使得满足这个条件的图像数据保留下来，其他的图像数据被赋值为零。

(2)比值法

通常采用比值法可以去除阴影、水体、冰雪、白泥地等多种干扰。首先判断干扰地物的各个波段的波谱特征，比如TM/ETM图像的阴影区第1波段明显大于第7波段，因此采用对于第7波段对第1波段的比值设定一个阈值进行去除的方法，植被采用第5波段比第4波段或者第3波段比第4波段的方法等。公式如下：

$y_i=\left(\frac{x_a}{x_b}\begin{array}{c}>\\ <\end{array}{C}_a\right)\&times;x_i$

其中，i＝0，...，n，n指所使用的遥感图像波段总数，x_i和y_i指分别指i波段去除“尖锐”信息前和后的的波段值，a∈[1，...，n]，C_a是常数，x_a，x_b是原始a、b波段对应的值。这个公式的目的就是：给定一个约束的条件，使得满足此约束条件的图像数据保留下来，其他的图像数据被赋值为零。

(3)Q值法

通常用于解决雪边或湖边湿地、干河道、冲积区、薄云等干扰。定义Q值如下：

Q＝(x_a×k^a-x_b×k^b)/x_c×k^c

其中，x_a，x_b，x_c为参与主成分分析的波段，k^a，k^b，k^c为参与主成分变化的x_a，x_b，x_c对应本征向量的值。

(4)光谱角法

对于薄云等来说，采用其他手段去除干扰难度较大，可以采用光谱角法。光谱角法把每一个多维空间点以其空间向量来表征，对比空间向量角的相似性。它是一种监督分类，如果判断某图像属于薄云类，则进行删除。要求对每一类别有一个已知参考谱。此参考谱可以是地面测得存入参考谱库的，也可以从具已知条件的图面单元做感兴趣区统计，存入参考谱库。公式如下：

$\left(\widehat{\mathrm{\&alpha;},\mathrm{\&beta;}}\right)\text{=arccos}\frac{(\mathrm{\&alpha;},\mathrm{\&beta;})}{\left|\mathrm{\&alpha;}\right|\left|\mathrm{\&beta;}\right|}$

式中(α，β)为n维向量α，β的内积，按内积定义

        (α，β)＝α₁β₁+α₂β₂+...+α_nβ_n

当α，β为列向量时，(α，β)＝α’β＝β’α

$\mathrm{\&alpha;}=\left[\begin{array}{c}{\mathrm{\&alpha;}}_1\\ {\mathrm{\&alpha;}}_2\\ .\\ .\\ .\\ {\mathrm{\&alpha;}}_n\end{array}\right];$ $\mathrm{\&beta;}=\left[\begin{array}{c}{\mathrm{\&beta;}}_1\\ {\mathrm{\&beta;}}_2\\ .\\ .\\ .\\ {\mathrm{\&beta;}}_n\end{array}\right]$

${\mathrm{\&alpha;}}^{\&prime;}\mathrm{\&beta;}=[{\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&alpha;}}_n]\left[\begin{array}{c}{\mathrm{\&beta;}}_1\\ {\mathrm{\&beta;}}_2\\ .\\ .\\ .\\ {\mathrm{\&beta;}}_n\end{array}\right]$

|α|、|β|为向量α、β的长度，

$\left|\mathrm{\&alpha;}\right|=\sqrt{(\mathrm{\&alpha;},\mathrm{\&alpha;})}=\sqrt{({\mathrm{\&alpha;}}_1^2,{\mathrm{\&alpha;}}_2^2\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&alpha;}}_n^2)}={\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^2\right)}^{1/2}$

求出α，β的内积和长度，便可求出cos，查表可得其夹角。

3、分区以获得正态分布的区域

首先，要对图像的波段直方图是否符合正态分布进行判断。

这里图像的波段直方图即为按波段值频率统计图。

具体地，遥感图像某一个波段像元xj，k(j＝1，m；k＝1，n)区间为[x0，xn]，统计这个窗口内的直方图，其中m，n为大于1的整数)，公式如下：

其中i∈[0，n]，xj，k＝xi为逻辑操作

如本领域公知的，判断波段直方图的技术例如为利用偏度系数和峰度系数进行判断。

需要对图像进行分区，使得分区后每个区内图像的每个波段直方图呈正态分布，这是因为每一个波段都近似服从正态，能够使得异常提取结果较优化。

4、区域内像元值规范化处理

首先利用直方图图示进行目估，如果遥感图像为一个M×N的窗口，这个窗口内遥感图像某一个波段像元x_j，k(j＝1，m；k＝1，n)区间为[x₀，x_n]，统计这个窗口内的直方图，公式如下：

$p_i=\underset{j=1,k=1}{\overset{m,n}{\mathrm{\&Sigma;}}}(x_{j,k}=A_i)\&times;1,$

其中i∈[0，u]，x_j，k＝A_i为逻辑操作，判断在坐标j，k处的像元值X是否等于A_i.

这里P_i表示在该M*N的窗口内，值Ai出现的次数。

我们取p_i的最大与最小值，即max(p_i)和min(p_i)，把出现频次最小的像元值表示为X(min(p_i))，把出现频次最大的像元值表示为X(max(p_i))。然后像元值拉伸到[0，255]的范围，其中中间值通过内插重新采样来获得。公式如下：

$y_{j,k}=\frac{255\&times;[X_{j,k}-X\left(\min \left(p_i\right)\right)]}{[X\left(\max \left(p_i\right)\right)-X\left(\min \left(p_i\right)\right)]}$

其中，y_j，k为原始图像某一波段像元x_j，k拉伸后的值，j＝1，m；k＝1，n。

通过上述规范化处理，把数据从原始范围拉伸到[0，255]的范围，可以粗略去除大气程辐射影响。

5、主成分分析

主成分分析或主成分变换既是一种信息提取方法，也是一种坐标变换方法。主成分分析的一般操作步骤为：第一步，移动坐标原点，使平均值为零；第二步，将坐标旋转，使一个坐标轴与数据具有最大分布的方向相符合，这个旋转后的新轴即第一主成分，它占有总变异(anomalies)的第一大分额。垂直于它的另一个坐标轴则代表其余变异的方向，这就是第二主成分。在多维空间里，这样的处理将继续进行，以确定一组直角坐标轴，这些轴逐渐将全部变异分配(消耗)掉。各主成分变异值的总和与变换前的变异值总和相等，即所谓的信息量守恒。

在N波段数据坐标轴转置的时候，协方差矩阵也将被变换，变换之后各波段之间的协方差变为零。

各点与其重心距离的平方和即为本征值λ₁，λ₂，...之和，这个和可表示为S。在一定意义上，可以说第一个成分“构成”的变异与总变异的比例是该成分的本征值与该和S之比λ₁/S，前两个成分“构成”的变异与总变异之比是(λ₁+λ₂)/S，如此等等。有时为了方便譬如可以说“前4个成分构成变异的p％”。

实际上，某主成分的本征值λ就是若消去该主成分后，对相应的本征向量所引入的均方误差值。

6、本征向量与蚀变异常特征之间符合程度的判断与处理

对于求出的本征向量，按照与参与主成分分析的各个波段进行对应，对应关系如下：

根据经验所得，符合蚀变异常特征的本征向量一般为第4个向量，因为第4个本征向量通常满足下述关系：V_c4的绝对值＞V_b4的绝对值并且V_c4的绝对值＞V_d4的绝对值，并且，V_c4与V_b4和V_d4符号相反，V_b4与V_d4符号相同。而该关系和地面波谱特征所反映的蚀变异常相符。

如果某一异常的特征为V_a4＞V_b4＜V_c4＞V_d4，那么V_a4、V_c4一定与V_b4、V_d4的符号相反，而V_a4与V_c4、V_b4与V_d4的符号相同。用于异常切割的本征向量4要求V_c4为正号，如果为负号，需要经过转换变成正号，公式如下：

$V_{c4}^T=(-1)\&times;V_{c4}$

其中，V_c4 ^T为V_c4经过符号转化后的结果。

7、异常切割

因为在主成分变换前，每个波段的直方图为正态分布，所以主成分变换后的异常主分量直方图也为正态分布。经过主成分变换后，异常集中在某一主分量上，利用切割水平把高值圈定出，作为提取出来的异常值，我们可以利用正态分布的有关理论进行异常的切割。

正态分布公式如下：

$f\left(X\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}{e}^{-\frac{x^2}{2{\mathrm{\&sigma;}}^2}}$

其中，X为随机变量，σ称为标准误差。对于主成分分析的多元分析来说，σ称为标准离差，定义如下：

$\mathrm{\&sigma;}=\sqrt{\frac{{\mathrm{\&Sigma;}}_{i=1}^n(x_i-\stackrel{\&OverBar;}{x})}{n}}$

n为样本数，

为均值，x_i为每个样本的值。在做异常切割或数据切割时便可借用σ这个表征正态分布曲线的尺度。例如，主分量分析结果可以把均值(X)理解为代表区域背景，以(X+kσ)作为异常下限和划分异常强度等级。一般取±4σ作为最小与最大限值。

切割异常时有了这一尺度可以减少主观任意性，异常分级是按下列公式计算的：

L＝127.5+k*SK；或L＝127.5+Sk*127.5/4；       H＝L+1

式中H、L分别为切割高、低门限值；k为倍数；SK为比例因子；和SK由主分量分析报告给出，或者根据经验确定。

二、感兴趣区域判断和选取

1、感兴趣区域初步判断和切割

感兴趣区域初步判断是指确定预定的已知矿床蚀变出露范围与异常范围的吻合程度，例如，可以利用GIS平台看真实面积是否吻合，利用异常切割确定异常分级，并切割出包含已知矿床异常的一块小图。所谓已知矿床类型可以是用户已知的感兴趣矿床，或是利用ASTER或高光谱数据提取的典型矿床类型。

切割方式可表示为如下公式：

$\left\{\begin{array}{c}{A}^q\&Element;A\\ A^q\stackrel{a}{\&DoubleLeftRightArrow;}M\end{array}\right.$

A^q为从图A切割下来的一块小图，A^q包括已知矿M。

2、门限化感兴趣区选取

首先，讨论一下提出利用门限值进行感兴趣区域选取的理论基础。

发明人经研究发现，在异常一定值的范围内，矿床蚀变反映在蚀变遥感异常一定水平切割上有一定的上下限值，可以认为小于某一个特定的切割水平的像素值是由主成分分析变化的噪声引起的，而认为大于某一个特定的切割水平的像素是有大范围的蚀变引起的，往往与矿化关系不大。因此，基于上述分析，在选择上限和下限为某一特征参数时，能够取一个感兴趣区波谱平均的值作为分类的一个特征波谱，可以在此基础上进行分类技术的处理。

利用通过上述感兴趣区域判断所切割下来的小图进行门限化感兴趣区选取，该方法采用异常阈值法进行选取的。对其应用双门限化感兴趣区选取。该方法对提取的异常采用一个上限和一个下限进行约束。公式如下：

$\left\{\begin{array}{c}I=A>{\mathrm{\&epsiv;}}_{\min }\\ I=A<{\mathrm{\&epsiv;}}_{\max }\end{array}\right.$

其中，I为选取的感兴趣区，A为连通的异常面积，ε_min与ε_max分别是选取的异常下限和异常上限，其可以根据经验进行选取，也可以利用自适应学习方法来进行选取。

三、相似度确定或分类

认为相似的地物有相似的光谱特征。因此对应某种矿床类型的异常采用监督分类得办法提取与其相似的异常信息。采用光谱角方法或利用最大似然法对感兴趣区内平均光谱特征进行监督分类，或者得到感兴趣区与已知矿床类型之间的相似度。下面分别加以说明。

1、光谱角法

光谱角法把每一个多维空间点以其空间向量来表征，对比空间向量角的相似性。光谱角法是一种监督分类。要求对每一类别有一个已知参考谱。此参考谱可以是地面测得存入参考谱库的，也可以从具已知条件的图面单元做感兴趣区统计，存入参考谱库。公式如下：

$\left(\widehat{\mathrm{\&alpha;},\mathrm{\&beta;}}\right)\text{=arccos}\frac{(\mathrm{\&alpha;},\mathrm{\&beta;})}{\left|\mathrm{\&alpha;}\right|\left|\mathrm{\&beta;}\right|}$

式中(α，β)为n维向量α，β的内积，按内积定义

                 (α，β)＝α₁β₁+α₂β₂+...+α_nβ_n

当α，β为列向量时(α，β)＝α’β＝β’α

$\mathrm{\&alpha;}=\left[\begin{array}{c}{\mathrm{\&alpha;}}_1\\ {\mathrm{\&alpha;}}_2\\ .\\ .\\ .\\ {\mathrm{\&alpha;}}_n\end{array}\right];$ $\mathrm{\&beta;}=\left[\begin{array}{c}{\mathrm{\&beta;}}_1\\ {\mathrm{\&beta;}}_2\\ .\\ .\\ .\\ {\mathrm{\&beta;}}_n\end{array}\right]$

${\mathrm{\&alpha;}}^{\&prime;}\mathrm{\&beta;}=[{\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&alpha;}}_n]\left[\begin{array}{c}{\mathrm{\&beta;}}_1\\ {\mathrm{\&beta;}}_2\\ .\\ .\\ .\\ {\mathrm{\&beta;}}_n\end{array}\right]$

|α|、|β|为向量α、β的长度，

$\left|\mathrm{\&alpha;}\right|=\sqrt{(\mathrm{\&alpha;},\mathrm{\&alpha;})}=\sqrt{({\mathrm{\&alpha;}}_1^2,{\mathrm{\&alpha;}}_2^2\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&alpha;}}_n^2)}={\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^2\right)}^{1/2}$

求出α，β的内积和长度，便可求出cos，查表可得其夹角。

由此，根据所得到夹角可判断是否包含已知矿床，获得得到相似度。

2、最大似然分类法

最大似乎分类，也称为贝叶斯(Bayes)分类，是基于概率统计的监督分类。

利用贝叶斯理论，假设为

属于模式类ω_i的后验概率。

根据贝叶斯理论可以计算出对象属于每一类的概率，并根据此结果将对象划分到最可能的类别中。决策分类依据为(二维情况)：

如果对象后验概率公式成立，则归为ω₁，反之归为ω₂中。

或者，如果是确定相似度，则属于各个类别的概率即反映了相似的程度信息。

通过前述处理，实现了对于识别遥感异常数据的分类处理和相似度确定处理，即为探矿提供了重要指导性信息。

四、进一步应用

除了上述处理外，本发明还可以提供进一步的信息。

比如，利用主成分分析提取的蚀变遥感异常信息包含了一个混合了的几种类型矿床的蚀变遥感异常信息，而利用监督分类技术提取出来的是与某一种类型矿床光谱特征相似的异常信息。那么这两种异常信息的交集，即既包含主成分分析提取的信息也包含监督分类提取的信息，为分类型蚀变遥感异常的信息。也就是说，因为主成分分析提取出来的是蚀变的信息，而光谱角技术提取出来的是包含与某种类型矿相关的信息，因此，即包含某种类型矿的信息又包含蚀变的信息，是某种类型矿床蚀变的信息。这种信息可以作为一种重要的找矿参数。公式如下：

info＝info_PCA.AND.info_SAM

其中，info为分类型蚀变遥感异常信息，info_PCA为主成分分析提取的信息，info_SAM为监督分类提取得信息。这样就实现了分类型蚀变遥感异常信息的提取。

五、蚀变遥感异常提取装置

图5示出了根据本发明一个实施例的蚀变遥感异常提取装置的框图。

如图5所示，蚀变遥感异常提取装置500，用于从输入的遥感数据提取表征与矿化相关的近矿蚀变岩石的信息，所述遥感数据包括多个波段的某地理区域内的像元值，包括下述部件：用于对输入的遥感数据进行前处理的部件510；用于从经所述前处理后得到的遥感数据所表征的遥感图像中切割出包含已知矿床类型的第一图像区域的感兴趣区域初步切割部件520；用于对该第一图像区域，利用预先确定的阈值从该第一图像区域中选择其中异常符合该阈值的第二图像区域的感兴趣区域细化选取部件530；以及针对该第二图像区域，基于该已知矿床类型的参考谱，确定第二图像区域和该已知矿床类型之间的相似程度的相似度确定部件540。

通过本发明的分类型矿床蚀变遥感异常信息的提取技术，能够去掉与某种类型矿化信息没有关系的干扰信息，而且只针对某种类型矿床进行蚀变遥感异常的提取，有了明确的针对性。本发明尤其适合应用于干旱半干旱地区进行矿产普查的遥感扫面、利用这种技术进行靶区的优选以及找矿，有极其重要的经济价值。而对于植被覆盖区，能用这技术进行特征导向的找矿勘察。本发明的应用与扩展应用主要有以下几个方面：

*提取出来的异常信息对找矿勘查有针对性

*利用异常信息在GIS平台进行多元找矿要素综合，优选找矿靶区

*分析不同成矿区带上异常空间分布，进行成矿潜力预测

*利用分布在不同区带上异常分布结合矿床学理论进行区域成矿规律研究

虽然上面已经结合附图示出并描述了本发明的一些实施例，但是本领域的技术人员应当理解，在不偏离本发明的原则和精神的情况下，可以对这些实施例做出变化和改变，所做的变化和改变仍然落在本发明及其等价物的保护范围之内。

Claims (6)

1.一种蚀变遥感异常提取方法，用于从输入的遥感数据提取表征与矿化相关的近矿蚀变岩石的信息，所述遥感数据包括多个波段的地理区域内的像元值，包括下述步骤：

对输入的遥感数据进行前处理；

从经所述前处理后得到的遥感数据所表征的遥感图像中切割出包含已知矿床类型的第一图像区域；

对该第一图像区域，利用预先确定的阈值从该第一图像区域中选择其中异常符合该阈值的第二图像区域；以及

针对该第二图像区域，基于该已知矿床类型的参考谱，确定第二图像区域和该已知矿床类型之间的相似程度，

其中，所述前处理步骤包括：

进行边界信息剔除，由此通过去除边界信息，使得各个波段的边界处信息重合，所谓边界信息是指，在获取遥感数据时，各个波段在边界处的信息不重合区域的信息，

判断有无云、水体、阴影区、白泥地、冰雪、植被、湿地、干河道、冲积扇九类常见干扰异常，如有则进行去除，去除方法为比值法、高端或低端切割法、Q值法、光谱角法中的一个或其组合；

统计遥感区域内的像元值，并以直方图表示，判断该直方图是否符合正态分布，如果不满足正态分布，则把该遥感区域划分为多个次级区域，使得每个次级区域内的每个波段直方图符合正态分布；

将符合正态分布的区域内的像元值统一转换到[0，255]的范围；

切割掉正态分布区域中的边缘部分，

对输入的数据进行主成分分析，并将输入的数据转换到通过主成分分析所确定的特征向量所表征的空间坐标系中，

所述选择第二图像区域的步骤包括：

如果切割第一图像区域的步骤中切割出包含已知矿床类型的第一图像区域，则利用下述公式确定第二图像区域：

$\left\{\begin{array}{c}I=A>{\mathrm{\&epsiv;}}_{\min }\\ I=A<{\mathrm{\&epsiv;}}_{\max }\end{array}\right.$

其中，I为选取的第二图像区域，A为连通的异常面积，ε_min与ε_max分别是预定的异常下限阈值和异常上限阈值。

2.如权利要求1所述的蚀变遥感异常提取方法，其中，该确定相似程度的步骤采用光谱角制图法或者贝叶斯分类法进行相似程度确定。

3.如权利要求1所述的蚀变遥感异常提取方法，其中，所述已知矿床类型是用户已知的感兴趣矿床，或是利用ASTER或高光谱数据提取的典型矿床类型。

4.一种蚀变遥感异常提取装置，用于从输入的遥感数据提取表征与矿化相关的近矿蚀变岩石的信息，所述遥感数据包括多个波段的地理区域内的像元值，包括下述部件：

用于对输入的遥感数据进行前处理的前处理部件；

用于从经所述前处理后得到的遥感数据所表征的遥感图像中切割出包含已知矿床类型的第一图像区域的感兴趣区域初步切割部件；

用于对该第一图像区域，利用预先确定的阈值从该第一图像区域中选择其中异常符合该阈值的第二图像区域的感兴趣区域细化选取部件；以及

针对该第二图像区域，基于该已知矿床类型的参考谱，确定第二图像区域和该已知矿床类型之间的相似程度的相似度确定部件，

其中，所述前处理部件进行下述操作：

边界信息剔除，由此通过去除边界信息，使得各个波段的边界处信息重合，所谓边界信息是指，在获取遥感数据时，各个波段在边界处的信息不重合区域的信息，

判断有无云、水体、阴影区、白泥地、冰雪、植被、湿地、干河道、冲积扇九类常见干扰异常，如有则进行去除，去除方法为比值法、高端或低端切割法、Q值法、光谱角法中的一个或其组合；

统计遥感区域内的像元值，并以直方图表示，判断该直方图是否符合正态分布，如果不满足正态分布，则把该遥感区域划分为多个次级区域，使得每个次级区域内的每个波段直方图符合正态分布；

将符合正态分布的区域内的像元值统一转换到[0，255]的范围；以及

切割掉正态分布区域中的边缘部分，

对输入的数据进行主成分分析，并将输入的数据转换到通过主成分分析所确定的特征向量所表征的空间坐标系中，

其中，如果感兴趣区域初步切割部件切割出包含已知矿床类型的第一图像区域，则感兴趣区域细化选取部件利用下述公式确定第二图像区域：

$\left\{\begin{array}{c}I=A>{\mathrm{\&epsiv;}}_{\min }\\ I=A<{\mathrm{\&epsiv;}}_{\max }\end{array}\right.$

其中，I为选取的第二图像区域，A为连通的异常面积，ε_min与ε_max分别是预定的异常下限阈值和异常上限阈值。

5.如权利要求4所述的蚀变遥感异常提取装置，其中，所述相似度确定部件采用光谱角制图法或者贝叶斯分类法进行相似程度确定。

6.如权利要求4所述的蚀变遥感异常提取装置，其中，所述已知矿床类型是用户已知的感兴趣矿床，或是利用ASTER或高光谱数据提取的典型矿床类型。

Images