CN102749296A - 多光谱遥感斑岩铜矿床蚀变分带信息提取方法及装置 - Google Patents

Abstract

本发明提供了一种多光谱遥感斑岩铜矿床蚀变分带信息提取方法及装置，其中所述方法包括对遥感数据进行预处理以获得基础数据的步骤，还包括如下步骤：A、对所述基础数据进行主成分分析法处理；B、对所述基础数据用比值法处理；C、对所述基础数据用斜率法处理；D、对所述基础数据用多项式法处理；E、对所述基础数据用吸收深度法处理；F、分别根据阈值筛选所述比值法、斜率发、多项式法和吸收深度法的处理结果；G、将所述步骤A与步骤F得到的结果进行逻辑与运算，得到的运算结果为蚀变分带信息提取结果。本发明可广泛应用于矿产资源遥感数据分析，准确度高。

Description

多光谱遥感斑岩铜矿床蚀变分带信息提取方法及装置

本发明中央级科研院所基本科研业务费项目（编号：K1207）、国家自然科学基金项目（编号：41102205）；地质大调查项目（编号：1212011120994）；国家重大基础研究计划项目（编号：2007CB411406-06）；国家科技支撑专题（编号：2006BAB07B00-7-2）的联合资助。

技术领域

本发明涉及一种地质遥感数据处理方法和装置，特别是从遥感数据中提取蚀变分带信息的方法和装置。

背景技术

随着社会经济的发展，人类社会对矿产资源的需要日益增大。矿产资源已经成为影响国家安全和可持续发展的重要战略资源。而随着持续不断的开采，寻找矿产资源的难度不断加大，因此急需找矿勘查技术的创新。近些年来，遥感技术在矿产勘查中得到逐步应用，取得了重要的成果，特别是在荒漠戈壁区发现了许多矿床和矿（化）点及大量蚀变异常。遥感技术在矿产勘查中的应用是以物质电磁波理论为基础，借助遥感技术获取数据,研究化学元素在地表或其它行星表面的分布、含量及迁移的科学，它的特点是快速、大范围获取数据。除能获取地球静态参数外，遥感地球化学还可以测量地球动态参数。

矿床的形成过程是某种有用元素的逐步富集过程，而这种成矿物质通常由成矿热液搬运和富集的。近矿围岩蚀变是成矿物质逐步富集成矿过程中留下的印迹。最常见的蚀变为：硅化、绢云母化、绿泥石化、云英岩化、矽卡岩化、白云岩化、重晶石化及锰铁碳酸盐化。地质学家断言绝大多数内生矿床都伴随有其围岩的交代蚀变现象,而且蚀变带范围大于矿体分布的范围数倍至数十倍。把近矿围岩蚀变现象作为找矿标志可追索到数百年前。根据围岩蚀变发现的大型金属、非金属矿床很多，例如：北美、俄罗斯的大部分斑岩铜矿，美国的许多白钨矿和犹他州的大铝矿，我国的铜官山铜矿，西澳大利亚的大型金矿，墨西哥的大铂矿，哈萨克斯坦的刚玉矿以及世界大多数锡矿等，这些实例均证明围岩蚀变现象作为找矿标志的重要意义。绝大多数矿床的围岩蚀变具有分带，以斑岩铜矿最为典型。

利用近矿蚀变围岩指导找矿的问题十分复杂，能够肯定的是蚀变岩石的发现，可以指示找矿方向，增加找到矿床的机会。地质学家认为，尽管有蚀变岩存在不一定有矿，然而大型、特大型内生矿床一般均有强烈且较大范围的围岩蚀变。蚀变遥感异常就是通过输血手段从遥感数据量化提取出来的、用以表征有可能与矿化最相关的近矿蚀变岩石的指示信息(Шехтман.热液矿床详细构造预测图[M].北京，地质出版社，1982.)。

蚀变遥感信息是隐含在岩石等背景噪声里的微弱信息，需要通过某种算法进行提取，一方面抑制遥感图像的背景噪声，另一方面增强微弱的蚀变遥感异常信息。目前，常用的提前蚀变遥感信息的方法采用主成分分析法。主成分分析法（principal componentanalysis(PCA)）主成分分析法是一种数学变换的方法,它把给定的一组相关变量通过线性变换转成另一组不相关的变量，这些新的变量按照方差依次递减的顺序排列。在数学变换中保持变量的总方差不变，使第一变量具有最大的方差，称为第一主成分，第二变量的方差次大，并且和第一变量不相关，称为第二主成分。依次类推，I个变量就有I个主成分。主成分分析法用较少的变量去解释原来资料中的大部分变异，将我们手中许多相关性很高的变量转化成彼此相互独立或不相关的变量。通常是选出比原始变量个数少，能解释大部分资料中的变异的几个新变量，即所谓主成分，并用以解释资料的综合性指标。

利用主成分分析方法从遥感数据中提取蚀变遥感信息存在的一个问题是：该方法仅能得知遥感数据中是否存在蚀变信息或某一特定矿物组合信息，而不能直接区分出遥感数据中的矿物或矿物组合的信息。

发明内容

为了解决现有提取蚀变遥感信息方法不能直接区分矿物或矿物组合的信息的问题，本发明提供了一种多光谱遥感斑岩铜矿床蚀变分带信息提取方法和装置。

本发明的技术方案如下：

多光谱遥感斑岩铜矿床蚀变分带信息提取方法，包括对遥感数据进行预处理以获得基础数据的步骤，还包括如下步骤：

A、对所述基础数据进行主成分分析法处理；

B、对所述基础数据用波段吸收特征法进行处理；所述波段吸收特征法包括比值法，即在基础数据光谱图中获得A/B的值；

其中，A为波段内吸收特征最显著的第一吸收谷的高值，B为所述第一吸收谷的最低值；

C、根据阈值筛选所述波段吸收特征法处理的结果；

D、将所述步骤A与步骤C得到的结果进行逻辑与运算，得到的运算结果为蚀变分带信息提取结果。

步骤B中所述波段吸收特征法还包括斜率法，即在所述基础数据光谱图中获得(A-B)/A的值。

步骤B中所述波段吸收特征法还包括多项式法，即在所述基础数据光谱图中获得(C-D)/C+(A-B)/A的值；其中，C为波段内吸收特征次显著的第二吸收谷的高值；D为所述第二吸收谷的最低值。

步骤B中所述波段吸收特征法还包括吸收深度法，即在所述基础数据光谱图中获得(A-B)/(A-C)的值；其中，波段内吸收特征次显著的C为第二吸收谷的高值。

步骤A中所述主成分分析法包括获取特征光谱步骤：对已知的蚀变矿物或蚀变矿物组合的光谱数据进行滤波，转换成对应于遥感数据波段的特征光谱；根据所述特征光谱对所述基础数据进行主成分分析处理。

所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法还包括如下步骤：根据步骤D得到的所述蚀变分带信息提取结果制成彩色图像。

所述对遥感数据进行预处理包括如下步骤：选取成像时间为春季或秋季的遥感图像作为原始数据；对属于ASTER数据的所述原始数据选择6、3和1波段组合，或对属于ETM数据的所述原始数据选择7、4和3波段组合。

所述对遥感数据进行预处理包括如下步骤：将所述遥感数据中不包含全部波段数据的边界区域去除。

所述对遥感数据进行预处理包括如下掩膜步骤：即将所述遥感数据中不需要进行处理的数据去除。

多光谱遥感斑岩铜矿床蚀变分带信息提取装置，包括预处理模块，用于对遥感数据进行预处理以获得基础数据；还包括如下模块：

主成分分析模块，与所述预处理模块连接，用于对所述基础数据进行主成分分析法处理；

波段吸收特征模块，与所述预处理模块连接，用于对所述基础数据用波段吸收特征法进行处理；所述波段吸收特征法包括比值法，即在基础数据光谱图中获得A/B的值；其中，A为波段内吸收特征最显著的第一吸收谷的高值，B为所述第一吸收谷的最低值；

筛选模块，与所述波段吸收特征模块连接，用于根据阈值筛选所述波段吸收特征模块处理的结果；

结果获取模块，分别与所述筛选模块和所述主成分分析模块连接，用于将所述主成分分析模块与筛选模块得到的结果进行逻辑与运算，得到的运算结果为蚀变分带信息提取结果。

所述波段吸收特征模块还用于在所述基础数据光谱图中获得(A-B)/A的值。

所述波段吸收特征模块还用于在所述基础数据光谱图中获得(C-D)/C+(A-B)/A的值；其中，C为波段内吸收特征次显著的第二吸收谷的高值；D为所述第二吸收谷的最低值。

所述波段吸收特征模块还用于在所述基础数据光谱图中获得(A-B)/(A-C)的值；其中，C为波段内吸收特征次显著的第二吸收谷的高值。

所述主成分分析模块还用于对已知的蚀变矿物或蚀变矿物组合的光谱数据进行滤波，转换成对应于遥感数据波段的特征光谱；根据所述特征光谱对所述基础数据进行主成分分析处理。

所述多光谱遥感斑岩铜矿床蚀变分带信息提取装置还包括图像生成模块，与所述结果获取模块连接，用于根据结果获取模块得到的所述蚀变分带信息提取结果制成彩色图像。

所述预处理模块还用于：选取成像时间为春季或秋季的遥感图像作为原始数据；对属于ASTER数据的所述原始数据选择6、3和1波段组合，或对属于ETM数据的所述原始数据选择7、4和3波段组合。

所述预处理模块还用于：将所述遥感数据中不包含全部波段数据的边界区域去除。

本发明的技术效果：

本发明采用的技术方案，除了对遥感数据进行主成分分析外，还引入了波段吸收特征法，并通过逻辑与的方式合并了两种方法的结果，兼顾了光谱的特征和深度，因此能够区分遥感数据中矿物或矿物组合的信息，实现了本发明的目的。

附图说明

图1为本发明多光谱遥感斑岩铜矿床蚀变分带信息提取方法的流程图。

图2为方解石光谱图的一个实例。

图3为本发明多光谱遥感斑岩铜矿床蚀变分带信息提取装置的原理图。

图4为采用现有主成分分析法获得的信息提取结果。

图5为本发明技术方案获得的信息提取结果。

具体实施方式

在对本发明的技术方案进行详细说明前，对涉及的名词进行解释。

ETM为Enhanced Thematic Mapper Plus的缩写，中文名称为增强型专题制图仪，是安放在美国陆地卫星上的传感器，以获取地表图像信息。

ASTER为Advanced Spaceborne Thermal Emission and Reflection Radiometer的缩写，中文名称为先进星载热发射仪和反辐射仪，是安放在美国Terra卫星上的高分辨解析地表图像传感器，其主要任务是通过14个频道获取整个地表的高分辨解析图像数据。

以下结合附图对本发明的技术方案进行详细说明。

图1显示了本发明多光谱遥感斑岩铜矿床蚀变分带信息提取方法的步骤，以下详细说明。

第一，预处理步骤。这一步骤对遥感数据（图像）进行处理，以形成可供后续步骤使用的基础数据。给出预处理步骤具体实施的一个实例如下。

（1）遥感图像数据选择

为了提高后续数据处理的准确性，必须对获得的遥感数据进行严格的筛选。要求遥感图像的时相尽可能地选择在植被发育弱、冰雪覆盖较少的季节，同时要求数据获取期间，该区域的天空云量较少。选择适合的遥感图像后，进一步地，要进行影像合成，一般对于ASTER数据一般选择6、3、1波段组合合成影像数据（即第一优化数据），ETM选择7、4、3波段组合合成影像数据（即第一优化数据）。这样可以明确看出植被、冰雪、云和水等地物。植被一般呈现绿色，冰雪是亮白色，或亮浅蓝色，云为亮白色，水为暗蓝色或暗黑色等。

（2）边界区域去除

判断遥感图像是否有边界区域。边界区域是指在缺失某个波段的信息的区域，或者说各波段信息不重合区域，这一区域往往处于遥感图像的边界。

如果判断存在上述边界区域，则去除边界区域。具体的方法就是对一个区域是否含有每一个波段信息进行判断，如果是，则附值为1；如果为否，则附值为0，生成一个二值图像，然后把每一个波段二值图像相乘形成一个新的二值图像，最后把每个波段与二值图像相乘，这样就去除了边界信息。具体公式如下：

$y_i=\underset{i=0}{\overset{n}{\∩}}(x_i>0)\·x_i$

其中n指所使用的遥感图像波段总数，i=1,…,n，xi和yi指分别指i波段去除波段前后的值。

（3）去干扰

为了消除云、水体、阴影区、白泥地、冰雪、植被、湿地、干河道、冲积扇等遥感数据中常见干扰，需要对这些干扰进行去除。去除方法选用比值法、高端或低端切割法、Q值法、光谱角法等。

1）高端或低端切割法

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

$y_i=\left(\begin{array}{c}{x}_b\begin{array}{c}>\\ <\end{array}{C}_b\end{array}\right)\&times;x_i$

其中，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(\begin{array}{c}\frac{x_a}{x_b}\begin{array}{c}>\\ <\end{array}{C}_a\end{array}\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)=\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，查表可得其夹角。

（4）判断直方图是否是正态分布

利用频域直方图，综合考虑每个波段，判断每个波段的直方图是否近似为正态分布。判断这个窗口内的各个波段直方图是否近似为正态分布。以下给出一个判断直方图是否是正态分布的计算判断方法。

设定一个M×N的窗口，这个窗口内遥感图像某一个波段像元值为x_j，k（j=1,m;k=1,n），像元均值为标准差为σ，利用偏度系数和峰度系数进行判断。

偏度系数满足公式：

$g=\frac{m\&times;n}{(m\&times;n-1)(m\&times;n-2)}\underset{j=1,k=1}{\overset{m,n}{\mathrm{\&Sigma;}}}{\left(\frac{x_{j,k}-\stackrel{\&OverBar;}{x}}{\mathrm{\&sigma;}}\right)}^3\&Element;[-{\mathrm{\&epsiv;}}_1,{\mathrm{\&epsiv;}}_1]$

其中ε₁为给定一个很小的正数。

峰度系数满足公式：

$f=\frac{m\&times;n(m\&times;n+1)}{(m\&times;n-1)(m\&times;n-2)(m\&times;n-3)}\underset{j=1,k=1}{\overset{m,n}{\mathrm{\&Sigma;}}}{\left(\frac{x_{j,k}-\stackrel{\&OverBar;}{x}}{\mathrm{\&sigma;}}\right)}^4-\frac{3{(m\&times;n-1)}^2}{(m\&times;n-2)(m\&times;n-3)}\&Element;[-{\mathrm{\&epsiv;}}_2,{\mathrm{\&epsiv;}}_2]$ 其中ε₂为给定一个很小的正数。

（5）分区处理

如果判断直方图不满足正态分布，则把图像区域划分为多个次级区域，并判断该次级区域是否满足正态分布，如果不满足，则进一步划分再次级区域，直到满足正态分布。

一般采用按照某类标准进行分类，比如地形复杂地区与地形简单地区分成不同的类，植被发育地区与植被不发育地区分出不同的类等。采用二值位图方式进行区分，公式如下：

$C_j^i=W^i\&times;x_j$

其中，Cⁱ _j为第j波段分区的第i区结果波段，Wi为第i类分类二值位图（在这个类里为1，不在这个位里为0），x_j为第j波段的原始数据值，j=1,2,…,n，n为波段总数，i=2,3，…，表示分类结果。

（6）线性拉伸

在遥感图像中开一个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}=x_i)\&times;1,$ 其中i∈[0,n]，x_j,k=x_i为逻辑操作

取p_i的最大与最小值，即max(p_i)和min(p_i)。然后把最小值作为0，最大值作为255，中间其他值按内插重新采样。公式如下：

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

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

（7）掩膜生成

掩膜作用是确定哪些数据需要参加计算，哪些不需要。要生成一个掩膜，其实就是生成一个二值图像，0值表示不需要参与计算的数据，1值表示需要参与计算的数据。生成掩膜一般采用波段的逻辑计算方法。

$\mathrm{\&tau;}=y_i\&CircleTimes;c_0$ 或 $\mathrm{\&tau;}=\left(\begin{array}{c}{x}_i\&CirclePlus;\begin{array}{c}{x}_j\\ c_1\end{array}\end{array}\right)\&CircleTimes;c_2$

τ为生成的掩膜，y_i为去除“尖锐”信息的数据，x_i，x_j为原始数据，为关系运算符（包括<、≤、≯、>、≥、≮、≠、＝等），为数学运算符（包括±、×、÷等），c₁，c₂，c₃为常数。

（8）形成基础数据

利用掩膜对每一个波段的数据进行处理，一般采用遥感数据的波段加上掩膜形成一种去干扰数据。一般是掩膜值为1值，遥感数据值不变，掩膜值为0值，遥感数据值变为0值。具体计算方法为：

y(i,j)=x(i,j)&τ

式中，y(i,j)为处理后的像素值，x(i,j)为处理前的像素值；i=1，2，...,n,j=1,2,....,n，为图像中的第i行第j行的数值。&为逻辑与运算符。经过处理后形成的基础数据为去除了干扰的数据。

第二，主成分分析法处理步骤。

在这一步骤中，先利用滤波函数把遥感数据的连续光谱变换成对应于遥感数据波段的特征光谱，具体方式如下：

蚀变矿物或蚀变矿物组合的波长向量其为X=[x₁,x₂，.....,x_n],对应光谱特征向量Y=[y₁,y₂,....,y_n]。

遥感图像每一点的图像采样光谱特征向量为

A（i,j)=[a₁(i,j),a₂(i,j),....,a_n(i,j)]。

那么蚀变矿物波谱特征导向算法如下：

B（i,j)=A(i,j)Y'=[a₁(i,j)*y₁,a₂(i,j)*y₂,.....,a_n(i,j)*y_n]

式中，B(i,j)为处理后的遥感图像每一点的图像采样光谱，A(i,j)为处理前的遥感图像每一点采样光谱，Y为蚀变矿物或蚀变矿物组合对应波长向量其为X=[x₁,x₂，.....,x_n]的光谱特征向量。蚀变矿物或蚀变矿物组合对应波长向量与遥感图像对应的波长向量一致。

然后，根据上述特征光谱，利用现有的主成分分析进一步处理（即在遥感数据中以特征光谱作为提取的主成分信息）：

主分量分析或主分量变换的原理为：第一步是移动坐标原点，使平均值为零。在这一步骤后，可将坐标旋转，使一个坐标轴与数据具有最大分布的方向相符合，这个旋转后的新轴即第一主分量，它占有总变异的第一大分额。垂直于它的另一个坐标轴则代表其余变异的方向，这就是第二主分量。在两维以上的多维空间里，这样的处理将继续进行，以确定一组直角坐标轴，这些轴逐渐将全部变异分配（消耗）掉，它并不能全部包含在一个次一级主分量中，而是有多少个原始参数就会有几个主分量。各主分量变异值的总和与变换前的变异值总和相等，这就是信息量守恒。

原始有几个波段的数据就影射到几个新的主分量上。各主分量是本征向量线性相加组合而成的。在数学中，就是找出一些新变量ξ₁，ξ₂，……，ξ_p，使它们是X的线性函数，而且彼此不相关，即

${\mathrm{\&xi;}}_i=\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{L}_{\mathrm{ik}}{X}_k$

实际上，就是求出p2个常数L_ik(i，k=1,…，p)按矩阵表示：

CL=λL

式中：L为本征向量，各L_ik是此本征向量的分量；λ为C矩阵的本征值。λ和L有以下特点：

称为迹，或总变异，对应于不同λ的L(即各主分量)线性不相关；且正交。

从线性代数中知道协方差矩阵C的本征多项式为det(λI-C)，此本征多项式的根λ都是协方差矩阵C的本征值。

计算过程为：

求协方差矩C

$C=\frac{1}{N}{\mathrm{XX}}^{\&prime;}$

求本征值λ

｜λI｜-C=0

求本征向量L

(λI-C)L=0

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

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

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

对于求出的本征向量，按照与参与主成分分析的各个波段进行对应，考虑符合蚀变异常特征的那个本征向量，一般为第4个向量。对应关系如下表：

如果某一异常的特征为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经过符合转化后的结果。

在主成分分析前，经过处理后使得每个波段的直方图为正态分布，转化后的异常主分量（即某本征向量）直方图也为正态分布，我们利用正态分布的有关理论进行异常的切割。正态分布公式如下：

$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σ*SF;或L=127.5+k*127.5/4；H=L+1

式中H、L分别为切割高、低门限值；k为倍数；σ为标准离差（在PCA（主成分分析）报告中会给出的）；SF为比例因子；k和SK由主分量分析报告给出。

第三，波段吸收特征法处理。

利用如下波段吸收特征法之一或多个进行处理：

比值法，即在基础数据光谱图中获得A/B的值。

斜率法，即在所述基础数据光谱图中获得(A-B)/A的值。

多项式法，即在所述基础数据光谱图中获得(C-D)/C+(A-B)/A的值。

吸收深度法，即在所述基础数据光谱图中获得(A-B)/(A-C)的值。

上述波段吸收特征法中，A为波段内吸收特征最显著的第一吸收谷的高值，B为该第一吸收谷的最低值，C为波段内吸收特征最显著的第二吸收谷的高值；D为第二吸收谷的最低值。图2是方解石的光谱，B处为光谱吸收特征最显著吸收谷的最低值，A为该吸收谷的高值（高值为吸收谷两侧的值），A-B为吸收深度。D为吸收特征第二显著的吸收谷的最低值，C为该吸收谷的一个高值。根据三个特征可以定义多项式、斜率、比值和吸收深度。

根据波段吸收特征法来进行信息提取，优势在于兼顾了光谱的特征和深度两个方面，因而使得本发明可以区分遥感信息中的矿物信息或矿物组合的信息。

第四，筛选波段吸收特征法处理结果。

对于波段吸收特征法提取出来的异常信息，根据光谱特征的关系，采用异常阈值法进行选取的。该方法对提取的异常采用一个上限和一个下限进行约束。公式如下：

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

其中，I为门限化结果，A为经过波段吸收处理后的值，ε_min与ε_max分别是选取的异常下限和异常上限。阈值的选取通常根据经验进行确定，例如根据选取数据量的5%-60%选择阈值的上、下限。

斑岩铜矿是指与斑状结构的中酸性浅成或超浅成的小侵入体-如花岗闪长斑岩、石英二长斑岩、石英斑岩等—有关，具有钾氢蚀变矿物晕和铜、钼、银、铅锌和硫地球化学晕的岩浆期后中-高温热液形成的细脉侵染型硫化物铜（钼或金）矿床（高合明，1995；袁见齐等，1985）。根据多龙矿集区遥感蚀变分带模型、土屋斑岩铜矿遥感蚀变分带模型、德兴斑岩铜矿遥感蚀变分带模型，由于岩浆岩期和热液期其都具有钾氢蚀变矿物晕，在岩体和围岩形成羟基蚀变矿物、钾化，中酸性岩浆大量的石英易形成硅化，从岩浆岩体中心向围岩外围，岩浆热液温度逐步降低，形成的蚀变矿物温度也从中高温向中低温过渡，地表矿物再次次生蚀变，一般岩体及其靠近岩体的围岩形成的是钾化、硅化、黑云母化、黑云母次生蚀变矿物绿泥石，在这个蚀变分带中由于硫含量高热液与钙质围岩反应也常形成硬石膏矿物；在其外围温度稍低，主要形成绢云母化、硅化、高岭石化；在绢英岩化外围是泥化分带，主要的蚀变矿物为高岭石化、高岭石化次生蚀变为埃洛石化、绿泥石化、伊利石化和蒙脱石化，最外围是青磐岩化带，形成碳酸盐化、绿泥石化、绿帘石化、蒙脱石化等蚀变。根据前文的矿物光谱特征分析和USGS提供的谱线，对斑岩铜矿的蚀变分带各个分带情况进行统计和总结，分析各种蚀变矿物的谱，根据其特征来总结各个分带的吸收特征，见下表。

斑岩铜矿蚀变分带矿物光谱特征

第五，结果合并步骤。

利用主成分分析法提取的信息包含了多种蚀变矿物的混合信息，而波段吸收特征法提取的信息定义为蚀变矿物的蚀变异常信息，那么根据蚀变矿物吸收特征提取的信息与主成分分析提取的信息进行叠加，限定蚀变矿物的信息，同时把蚀变分带中不同矿物的蚀变信息进行叠加处理。公式如下：

info¹=info_PCA&info¹ _RA&info¹ _SL&info¹ _PO&info¹ _AB

式中中，info¹为遥感蚀变矿物信息，info_PCA为主成分分析提取的信息，info¹ _RA为比值法提取的信息，info¹ _SL为斜率法提取的信息，info¹ _PO为多项式法提取的信息，info¹ _AB为相关吸收法提取的信息，&为逻辑与运算符。

info=info¹&info²&...&infoⁿ

info为遥感蚀变分带信息，info¹，info²，...，infoⁿ为蚀变分带所包含蚀变矿物的信息，&为逻辑与运算符。这样就实现了分类型蚀变遥感异常信息的提取。

第六，生成图像。

经过以上基本处理后，提取了各个蚀变分带遥感信息，底图采用RGB合成一幅假彩色图像，各遥感蚀变分带信息利用不同的色彩叠加在RGB底图上。这样做出来的图像会反映遥感蚀变分带的精细信息，更加适合人的视觉习惯，简单、易用。

图4是采用现有主成分分析法获得的结果。图5是采用本发明方法处理同一遥感数据后获得的结果。与图4相比，图5中显示了从遥感数据中能够提取出不同矿物或矿物组合的信息，显示了本发明处理结果的优势。

图3显示了本发明多光谱遥感斑岩铜矿床蚀变分带信息提取装置的原理。

图3所示多光谱遥感斑岩铜矿床蚀变分带信息提取装置包括：预处理模块；与预处理模块连接的主成分分析模块、波段吸收特征模块；与波段吸收特征模块连接的筛选模块；分别与筛选模块和主成分分析模块连接的结果获取模块；及与结果获取模块连接的图像生成模块。

其中，预处理模块执行的步骤与图1中“预处理”步骤相同。主成分分析模块执行的步骤与图1中“主成分分析法处理”步骤相同。波段吸收特征模块执行的步骤与图1中“波段吸收特征法处理”步骤相同。筛选模块执行的步骤与图1中“筛选波段吸收特征法处理结果”步骤相同。结果获取模块执行的步骤与图1中“结果合并”步骤相同。图像生成模块执行的步骤与图1中“生成图像”步骤相同。

值得注意的是，以上所述仅为本发明的较佳实施例，并非因此限定本发明的专利保护范围，本发明还可以采用等同技术进行替换。故凡运用本发明的说明书及图示内容所作的等效变化，或直接或间接运用于其他相关技术领域均同理皆包含于本发明所涵盖的范围内。

Claims (17)

1.多光谱遥感斑岩铜矿床蚀变分带信息提取方法，包括对遥感数据进行预处理以获得基础数据的步骤，其特征在于还包括如下步骤：

A、对所述基础数据进行主成分分析法处理；

B、对所述基础数据用波段吸收特征法进行处理；所述波段吸收特征法包括比值法，即在基础数据光谱图中获得A/B的值；

其中，A为波段内吸收特征最显著的第一吸收谷的高值，B为所述第一吸收谷的最低值；

C、根据阈值筛选所述波段吸收特征法处理的结果；

D、将所述步骤A与步骤C得到的结果进行逻辑与运算，得到的运算结果为蚀变分带信息提取结果。

2.根据权利要求1所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法，其特征在于：步骤B中所述波段吸收特征法还包括斜率法，即在所述基础数据光谱图中获得(A-B)/A的值。

3.根据权利要求1所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法，其特征在于：步骤B中所述波段吸收特征法还包括多项式法，即在所述基础数据光谱图中获得(C-D)/C+(A-B)/A的值；其中，C为波段内吸收特征次显著的第二吸收谷的高值；D为所述第二吸收谷的最低值。

4.根据权利要求1所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法，其特征在于：步骤B中所述波段吸收特征法还包括吸收深度法，即在所述基础数据光谱图中获得(A-B)/(A-C)的值；其中，波段内吸收特征次显著的C为第二吸收谷的高值。

5.根据权利要求1至4之一所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法，其特征在于：步骤A中所述主成分分析法包括获取特征光谱步骤：对已知的蚀变矿物或蚀变矿物组合的光谱数据进行滤波，转换成对应于遥感数据波段的特征光谱；根据所述特征光谱对所述基础数据进行主成分分析处理。

6.根据权利要求5所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法，其特征在于：还包括如下步骤：根据步骤D得到的所述蚀变分带信息提取结果制成彩色图像。

7.根据权利要求6所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法，其特征在于：所述对遥感数据进行预处理包括如下步骤：选取成像时间为春季或秋季的遥感图像作为原始数据；对属于ASTER数据的所述原始数据选择6、3和1波段组合，或对属于ETM数据的所述原始数据选择7、4和3波段组合。

8.根据权利要求7所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法，其特征在于：所述对遥感数据进行预处理包括如下步骤：将所述遥感数据中不包含全部波段数据的边界区域去除。

9.根据权利要求8所述多光谱遥感斑岩铜矿床蚀变分带信息提取方法，其特征在于：所述对遥感数据进行预处理包括如下掩膜步骤：即将所述遥感数据中不需要进行处理的数据去除。

10.多光谱遥感斑岩铜矿床蚀变分带信息提取装置，包括预处理模块，用于对遥感数据进行预处理以获得基础数据；其特征在于还包括如下模块：

主成分分析模块，与所述预处理模块连接，用于对所述基础数据进行主成分分析法处理；

波段吸收特征模块，与所述预处理模块连接，用于对所述基础数据用波段吸收特征法进行处理；所述波段吸收特征法包括比值法，即在基础数据光谱图中获得A/B的值；其中，A为波段内吸收特征最显著的第一吸收谷的高值，B为所述第一吸收谷的最低值；

筛选模块，与所述波段吸收特征模块连接，用于根据阈值筛选所述波段吸收特征模块处理的结果；

结果获取模块，分别与所述筛选模块和所述主成分分析模块连接，用于将所述主成分分析模块与筛选模块得到的结果进行逻辑与运算，得到的运算结果为蚀变分带信息提取结果。

11.根据权利要求10所述多光谱遥感斑岩铜矿床蚀变分带信息提取装置，其特征在于：所述波段吸收特征模块还用于在所述基础数据光谱图中获得(A-B)/A的值。

12.根据权利要求10所述多光谱遥感斑岩铜矿床蚀变分带信息提取装置，其特征在于：所述波段吸收特征模块还用于在所述基础数据光谱图中获得(C-D)/C+(A-B)/A的值；其中，C为波段内吸收特征次显著的第二吸收谷的高值；D为所述第二吸收谷的最低值。

13.根据权利要求10所述多光谱遥感斑岩铜矿床蚀变分带信息提取装置，其特征在于：所述波段吸收特征模块还用于在所述基础数据光谱图中获得(A-B)/(A-C)的值；其中，C为波段内吸收特征次显著的第二吸收谷的高值。

14.根据权利要求10至13之一所述多光谱遥感斑岩铜矿床蚀变分带信息提取装置，其特征在于：所述主成分分析模块还用于对已知的蚀变矿物或蚀变矿物组合的光谱数据进行滤波，转换成对应于遥感数据波段的特征光谱；根据所述特征光谱对所述基础数据进行主成分分析处理。

15.根据权利要求14所述多光谱遥感斑岩铜矿床蚀变分带信息提取装置，其特征在于：还包括图像生成模块，与所述结果获取模块连接，用于根据结果获取模块得到的所述蚀变分带信息提取结果制成彩色图像。

16.根据权利要求15所述多光谱遥感斑岩铜矿床蚀变分带信息提取装置，其特征在于：所述预处理模块还用于：选取成像时间为春季或秋季的遥感图像作为原始数据；对属于ASTER数据的所述原始数据选择6、3和1波段组合，或对属于ETM数据的所述原始数据选择7、4和3波段组合。

17.根据权利要求16所述多光谱遥感斑岩铜矿床蚀变分带信息提取装置，其特征在于：所述预处理模块还用于：将所述遥感数据中不包含全部波段数据的边界区域去除。

Images