CN102426625A - 基于aster卫星数据的斑岩铜典型蚀变带矿物信息提取方法 - Google Patents

Abstract

本发明涉及一种基于ASTER卫星数据的斑岩铜典型蚀变带矿物信息提取方法，即泥化-似千枚岩化蚀变带矿物组合信息提取方法，其特征在于它包含以下步骤：第一步：ASTER数据的读入；第二步：数据质量分析与控制；第三步：数据预处理；第四步：构建4波段(B1，B4，B6，B7)主成分分析模型，获取特征向量矩阵；第五步：特征向量分析，确定蚀变矿物信息所在主分量。该方法简洁、实用、可操作性强，可以半定量的提取出斑岩铜矿典型蚀变带泥化-似千枚岩化矿物组合信息，而且，所用ASTER卫星数据已覆盖全球90％以上陆地，具有在全球斑岩铜矿成矿带(区)推广应用的价值，可为寻找该类矿产提供重要的直接指示信息。

Description

基于ASTER卫星数据的斑岩铜典型蚀变带矿物信息提取方法

技术领域

本发明涉及一种基于ASTER(Advanced Spaceborne Thermal Emissionand Reflection Radiometer，即高级空间热放射与反射测量仪)卫星数据的斑岩铜典型蚀变带矿物信息提取方法，尤其是涉及一种基于ASTER卫星数据的斑岩铜矿泥化-似千枚岩化蚀变带矿物组合信息提取方法，属于遥感地质技术，适用于利用ASTER卫星数据半定量识别斑岩铜矿床典型蚀变带中泥化-千枚岩化蚀变带矿物组合。

背景技术

ASTER卫星数据是介于传统TM(Thematic Mapper，即专题制图仪)/ETM+(Enhanced Thematic Mapper Plus，即增强的专题制图仪+)等多光谱数据和高光谱数据之间的多光谱数据，国外也有专家称之为高光谱数据。高光谱数据识别矿物可以实现单矿物或矿物组合定量识别，但是，目前没有在轨运行的适用于地质矿产相关的矿物蚀变的数据源；现阶段开展的蚀变矿物信息提取，均以TM/ETM+和ASTER数据为数据源，前者提取铁染和羟基异常信息技术方法相对成熟，ASTER数据由于其包含3个可见光(0.52-0.86um)、6个短波红外(1.60-2.43um)和5个热红外(8.125-11.65um)共计14个波段，可以提取更加细微的蚀变矿物和组合信息，国内外相关专家开展了大量的针对性研究与应用，目前主要研建了包括Fe³⁺离子、Al-OH基团、CO₃ ²离子团和Mg-OH基团类矿物信息提取方法。但是，还没有形成就某一特定成矿类型典型蚀变带信息提取的技术方法。

斑岩铜矿床占全球铜矿产量的70％以上，且多为大型-特大型矿床，其具有典型的围岩蚀变分带性，从中心向外位依次为：(1)钾质蚀变带：蚀变矿物主要为黑云母和钾长石；(2)似千枚岩化带(石英-绢云母化带)：蚀变矿物主要为石英和绢云母；(3)泥质蚀变带：蚀变矿物主要为高岭石、蒙脱石、明矾石和石英；(4)青盘岩化带：蚀变矿物主要为绿泥石、绿帘石和方解石。通过前人的研究发现上述4个蚀变带并不是每个矿床都发育齐全的，但以石英、绢云母构成的似千枚岩化带，几乎在所有斑岩型铜矿中均有广泛发育，其强度、范围和矿化的规模有直接的关系。

安第斯巨型成矿带是全球最大的斑岩铜矿成矿带，该带内广泛分布大型-特大型斑岩铜矿床，铜产量占全球铜产量的40％以上，且均具有典型斑岩铜矿围岩蚀变带，矿体与似千枚岩化带密切相关。据此，只要能研发提取出似千枚岩化带类蚀变矿物或矿物组合的方法，将为斑岩型铜矿的找矿工作提供直接的指示信息。

泥质岩化带和似千枚岩化带主要蚀变矿物包括石英、绢云母、高岭石、蒙脱石、伊利石和明矾石等，后5种矿物均具有在2.2吸收谷(对应ASTER数据的B6)，B4高反射峰以及B7相对较高反射峰特征。据此，构造4个波段(B1，B 4，B 6，B 7)主成分分析模型，通过对比分析其他典型蚀变矿物波谱曲线特征，发现其他矿物波谱曲线不具备该特征，这样可以利用该4个波段主成分分析模型提取出泥化-似千枚岩化带蚀变矿物组合信息。

发明内容

本发明的目的是，针对目前具有巨大铜矿找矿潜力的斑岩型铜矿还没有提取其典型蚀变带矿物组合信息的问题，提出基于ASTER卫星数据的斑岩铜典型蚀变带矿物信息提取方法，即基于ASTER数据源4个波段(B1，B4，B6，B7)主成分分析模型。在斑岩铜矿成矿区/带均可以采用该方法，提取出典型蚀变带主要矿物组合信息，为斑岩铜矿的找矿工作提供重要的直接指示信息。

现有的TM/ETM+数据和ASTER数据开展的蚀变矿物提取技术均为基于蚀变矿物波谱曲线特征开展大类蚀变异常信息提取，还没有形成基于地质成矿系统特别是具有重大找矿意义的典型蚀变带主要蚀变信息提取方法，有必要研发基于ASTER卫星数据的斑岩铜矿泥化-似千枚岩化蚀变带矿物组合信息提取方法。首先，基于USGS标准波谱数据库重建蚀变带主要蚀变矿物对应ASTER数据B1-B9的波谱曲线，综合对比分析斑岩铜矿围岩蚀变带蚀变矿物波谱曲线特征，并构建4个波段(B1，B4，B6，B7)主成分分析模型，研发基于ASTER数据的斑岩铜矿泥化-似千枚岩化蚀变带主要矿物组合信息提取方法，为在斑岩铜矿的找矿工作提供直接指示信息。

本发明针对上述需要解决的技术问题，提出了一种基于ASTER卫星数据的斑岩铜典型蚀变带矿物信息提取方法，如图1所示，具体包括如下步骤：

步骤一：ASTER数据的读入；

步骤二：数据质量分析与控制；

步骤三：数据预处理；

步骤四：构建4波段(B1，B4，B6，B7)主成分分析模型，获取特征向量矩阵；

步骤五：特征向量矩阵分析，确定蚀变矿物信息所在主分量；

步骤六：统计出蚀变矿物信息所在主分量象元灰度值平均值和标准离差(σ)；

步骤七：根据步骤六的统计结果，确定阈值，提取蚀变矿物组合信息。

其中步骤三所述的“数据预处理”，其具体方法为，

首先，利用ASTER数据管理与分发官网提供的Crosstalk 3.0软件，消除传感器本身设置引起的串扰因素；然后，采用FLAASH(Fast Line-of-SightAtmospheric Analysis of Spectral Hypercubes)大气校正模型，完成ASTER数据大气校正；最后，消除黑边和边框、植被、冰雪、水体(海洋、湖和湿地)边缘、云等引起的干扰信息，并制作掩膜层。

其中，步骤四中所述的4波段(B1，B4，B6，B7)主成分分析模型构建的主要依据是：基于USGS(United States Geological Survey，即美国地质勘探局)标准波谱数据库重建斑岩铜矿蚀变带主要蚀变矿物对应ASTER数据B1至B9的波谱曲线(如图2)，综合对比分析斑岩铜矿围岩蚀变带蚀变矿物波谱曲线特征，确定白(绢)云母(似千枚岩化带主要蚀变矿物)、高岭石、蒙脱石、伊利石和明矾石(泥质岩化带主要蚀变矿物)等蚀变矿物均具有在2.2吸收谷(对应ASTER数据的B6，图2b)、B4高反射峰以及B7反射峰特征，而其他蚀变矿物不具备该特征，进而利用这一特征构建4波段(B1，B4，B6，B7)主成分分析方法，并通过阈值的确定，提取出斑岩铜矿床泥化-似千枚岩化蚀变带矿物组合信息。

步骤四中所述的构建4波段(B1，B4，B6，B7)主成分分析模型，获取特征向量矩阵，具体流程是：

(1)计算每一波段的均值、方差

①各波段的均值为：

第一主成分波段(B1)均值：

其中，n为第一波段影像的行列号，n²即为影像所包含的元素总数；

表示所有元素的总和。第二主成分波段(B4)均值

第三主成分波段(B6)均值

第四主成分波段(B7)均值

依次类推；

②按照样本方差计算公式：

${\mathrm{\σ}}_i^2=\frac{1}{n^2-1}\underset{j=1}{\overset{n^2}{\mathrm{\Σ}}}{(x_j-{\stackrel{\‾}{x}}_i)}^2,i=\mathrm{1,2,3,4}$

计算各波段的方差值，其中n为该波段影像的行列号，n²即为影像所包含的元素总数；由此得到4个波段分别对应的方差值：σ₁ ²、σ₂ ²、σ₃ ²、σ₄ ²。

(2)计算协方差矩阵：

首先，逐一计算协方差矩阵的元素，随机试验的自协方差函数为：

$C(i,j)=E\left\{(x_i(k,l)-{\stackrel{\‾}{x}}_i)(x_j(k,l)-{\stackrel{\‾}{x}}_j)\right\}$

其中x_i(k，l)为像素点(k，l)在i波段的值，

为i波段的平均值，x_j(k，l)为像素点(k，l)在j波段的值，为j波段的平均值。上述公式也可以表示为：

$C(i,j)=\frac{1}{4^2-1}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}\underset{l=1}{\overset{n}{\mathrm{\Σ}}}(x_i(k,l)-{\stackrel{\‾}{x}}_i)(x_j(k,l)-{\stackrel{\‾}{x}}_j)$

由此可以求出4个波段对应的协方差值，分别是对应第一主成分波段B1的C(1，1)，C(1，2)，C(1，3)，C(1，4)；对应第二主成分波段B4的C(2，1)，C(2，2)，C(2，3)，C(2，4)；对应第三主成分波段B6的C(3，1)，C(3，2)，C(3，3)，C(3，4)；对应第四主成分波段B7的C(4，1)，C(4，2)，C(4，3)，C(4，4)；从而得到协方差矩阵为：

$C=\left[\begin{array}{c}C\left(\mathrm{1,1}\right),C\left(\mathrm{1,2}\right),C\left(\mathrm{1,3}\right),C\left(\mathrm{1,4}\right)\\ C\left(\mathrm{2,1}\right),C\left(\mathrm{2,2}\right),C\left(\mathrm{2,3}\right),C\left(\mathrm{2,4}\right)\\ C\left(\mathrm{3,1}\right),C\left(\mathrm{3,2}\right),C\left(\mathrm{3,3}\right),C\left(\mathrm{3,4}\right)\\ C\left(\mathrm{4,1}\right),C\left(\mathrm{4,2}\right),C\left(\mathrm{4,3}\right),C\left(\mathrm{4,4}\right)\end{array}\right];$

(3)计算矩阵C的特征值与特征向量

矩阵C的特征方程是|λE-C|＝0，它的根就是C的全部特征根。对每个特征值λ₀，其次方程组(λ₀E-C)x＝0的每一个非零解都是C的对应于λ₀的特征向量。求出C的特征值为λ₁，λ₂，λ₃，λ₄，要求λ₁＞λ₂＞λ₃＞λ₄。相应的特征向量为μ₁，μ₂，μ₃，μ₄；

(4)计算图像的主分量

用相应的特征向量作为行，就得到变换矩阵A，根据主成分变换的数学原理，A是X空间的协方差矩阵C的转置矩阵(T)，即

A＝C^T

用该矩阵对每个像素进行变换，可求得图像的主分量：变换后，用4维向量表示像素的第一个元素组成了原图像的第一主分量，第二个元素组成了原图像的第二主分量，同样，第三、四个元素组成了原图像的第三、四主分量，即得到4个主成分分量PC1，PC2，PC3，PC4。

其中步骤五所述的“特征向量矩阵分析，确定蚀变矿物信息所在主分量”，其方法是：分析特征向量矩阵中的各主分量，看哪个主分量具有如下特征：即B7和B4的贡献系数与B6的贡献系数符号相反，而且B6的贡献系数绝对值相对较高。经过大量的试验表明，在第四主分量(PC4)具有该特征，进而判定所要提取的蚀变矿物组合信息所在主分量，即第四主分量。

其中步骤六所述的“统计出蚀变矿物信息所在主分量象元灰度值平均值和标准离差(σ)”时，应注意不要统计掩膜层对应的区域。其中平均值计算公式为：

$\stackrel{\‾}{X}=(X_1+X_2+...+X_n)/n$

是主分量图像中所有像素值的算术平均值

标准离差(σ)计算公式为：

$\mathrm{\σ}=\sqrt{[{(\stackrel{\‾}{X}-X_1)}^2+{(\stackrel{\‾}{X}-X_2)}^2+...+{(\stackrel{\‾}{X}-X_n)}^2]/n};$

其中步骤七所述的“确定阈值”，前人在TM/ETM+数据异常信息提取方法研究中做过相应的研究与大量的试验统计分析得出异常分级的经验值，在本方法的研制过程中，也对数十景TM/ETM+数据和百余景ASTER进行了相应的研究，发现采用1.5σ(1.5倍标准离差，以下类似)、2.0σ和2.5σ确定阈值，切割出的三级蚀变矿物组合信息，与传统经验值所确定的阈值或分级标准基本一致。

本发明的优点及功效在于：该方法简洁、实用、可操作性强，可以半定量的提取出斑岩铜矿典型蚀变带泥化-似千枚岩化矿物组合信息，而且，所用ASTER卫星数据已覆盖全球90％以上陆地，具有在全球斑岩铜矿成矿带(区)推广应用的价值，可为寻找该类矿产提供重要的直接指示信息。

附图说明

图1为本发明涉及的一种基于ASTER卫星数据的斑岩铜矿泥化-似千枚岩化蚀变带矿物组合信息提取方法的实现流程；

图2(a)、(b)、(c)为本发明涉及的基于USGS标准波谱数据库重建蚀变带主要蚀变矿物对应ASTER数据B1-B9的波谱曲线。

具体实施方式

为了更好的说明本发明涉及的基于ASTER卫星数据的斑岩铜矿泥化-似千枚岩化蚀变带矿物组合信息提取方法，在秘鲁南部斑岩铜矿矿集区阿里帕基地区进行试验，具体步骤如下：

第一步，读入ASTER数据(ID：ASTL1A 0607271504480607300378)，制作B1至B9的数据集；

第二步，数据质量分析与控制，采用B2、B3、B1(RGB)假彩色合作，同时检查B4至B9，发现数据质量良好；

第三步，数据预处理，

(1)利用Crosstalk 3.0软件消除串扰效应；

(2)采用FLAASH大气校正模型，对ASTER数据进行大气校正；

(3)干扰地物信息去除，采用B3与B6与B7逻辑“与”运算，获取边框和黑边数据，B6/B1获取水体数据，B3/B2获取植被数据，然后将这3类数据再做逻辑“与”运算获取需要去除的干扰地物信息，并制作掩膜层。

第四步，构建4波段(B1，B4，B6，B7)主成分分析模型，进行主成份分析，获取特征向量矩阵(如下表所示)；

ASTER 1、4、6、7主成分分析特征矩阵

	Band 1	Band 4	Band 6	Band 7
					PC 1	0.524873	0.490004	0.467668	0.515453
PC 2	0.783112	-0.593528	-0.166453	-0.082177
					PC 3	-0.329839	-0.636640	0.459523	0.524151
PC 4	-0.049497	0.048070	-0.736490	0.672920

第五步，分析特征向量矩阵，依据上述判别方法，可见PC 4为蚀变矿物信息所在主分量；

第六步，统计出蚀变矿物信息所在主分量象元灰度值平均值和标准离差(σ)，应注意不要统计掩膜层对应的区域。其中平均值计算公式为：

$\stackrel{\‾}{X}=(X_1+X_2+...+X_n)/n$

是主分量图像中所有像素值的算术平均值

标准离差(σ)计算公式为：

$\mathrm{\σ}=\sqrt{[{(\stackrel{\‾}{X}-X_1)}^2+{(\stackrel{\‾}{X}-X_2)}^2+...+{(\stackrel{\‾}{X}-X_n)}^2]/n};$

根据上述的计算公式，计算出PC4主分量象元灰度值的平均值为127.60和标准离差(σ)为41.14；

第七步，根据步骤六的计算结果，确定需提取的一、二和三级蚀变矿物组合信息的阈值分别为190、210和231，据此分别提取出一、二和三级蚀变矿物组合信息。

试验结果表明，在已知的塞罗维德(Cerra Verde)大型斑岩铜矿床周边具有很有的分级蚀变矿物组合信息，与已知地质背景资料吻合，另外在其南部和东南部环形构造边缘也出现了类似的分级蚀变矿物组合信息，从高分辨卫星(QuickBrid卫星)影像图上，可以清楚的显示出个别地区已开始露天采矿，大部分区域正在开展系统的地面勘探工程，证明所研发的基于ASTER卫星数据的斑岩铜矿泥化-似千枚岩化蚀变带矿物组合信息提取方法可靠，具有很好的应用价值。

Claims (2)

1.一种基于ASTER卫星数据的斑岩铜典型蚀变带矿物信息提取方法，即泥化-似千枚岩化蚀变带矿物组合信息提取方法，其特征在于它包含以下步骤：

第一步：ASTER数据的读入；

第二步：数据质量分析与控制；

第三步：数据预处理；

第四步：构建4波段(B1，B4，B6，B7)主成分分析模型，获取特征向量矩阵；

第五步：特征向量分析，确定蚀变矿物信息所在主分量；

第六步：统计出蚀变矿物信息所在主分量象元灰度值的平均值和标准离差(σ)；

第七步：根据第六步的统计结果，确定阈值，提取蚀变矿物组合信息。

2.根据权利要求1所述的基于ASTER卫星数据的斑岩铜矿典型蚀变带矿物信息提取方法，其特征在于：步骤四中所述的构建4波段(B1，B4，B6，B7)主成分分析模型，获取特征向量矩阵，具体流程是：

(1)计算每一波段的均值、方差

①各波段的均值为：

第一主成分波段(B1)均值：

其中，n为第一波段影像的行列号，n²即为影像所包含的元素总数；

表示所有元素的总和；第二主成分波段(B4)均值

第三主成分波段(B6)均值

第四主成分波段(B7)均值

依次类推；

②按照样本方差计算公式：

${\mathrm{\σ}}_i^2=\frac{1}{n^2-1}\underset{j=1}{\overset{n^2}{\mathrm{\Σ}}}{(x_j-{\stackrel{\‾}{x}}_i)}^2,i=\mathrm{1,2,3,4}$

计算各波段的方差值，其中n为该波段影像的行列号，n²即为影像所包含的元素总数；由此得到4个波段分别对应的方差值：σ₁ ²、σ₂ ²、σ₃ ²、σ₄ ²；

(2)计算协方差矩阵：

首先，逐一计算协方差矩阵的元素，随机试验的自协方差函数为：

$C(i,j)=E\left\{(x_i(k,l)-{\stackrel{\‾}{x}}_i)(x_j(k,l)-{\stackrel{\‾}{x}}_j)\right\}$

其中x_i(k，l)为像素点(k，l)在i波段的值，

为i波段的平均值，x_j(k，l)为像素点(k，l)在j波段的值，

为j波段的平均值；上述公式也可以表示为：

$C(i,j)=\frac{1}{4^2-1}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}\underset{l=1}{\overset{n}{\mathrm{\Σ}}}(x_i(k,l)-{\stackrel{\‾}{x}}_i)(x_j(k,l)-{\stackrel{\‾}{x}}_j)$

由此可以求出4个波段对应的协方差值，分别是对应第一主成分波段B1的C(1，1)，C(1，2)，C(1，3)，C(1，4)；对应第二主成分波段B4的C(2，1)，C(2，2)，C(2，3)，C(2，4)；对应第三主成分波段B6的C(3，1)，C(3，2)，C(3，3)，C(3，4)；对应第四主成分波段B7的C(4，1)，C(4，2)，C(4，3)，C(4，4)；从而得到协方差矩阵为：

$C=\left[\begin{array}{c}C\left(\mathrm{1,1}\right),C\left(\mathrm{1,2}\right),C\left(\mathrm{1,3}\right),C\left(\mathrm{1,4}\right)\\ C\left(\mathrm{2,1}\right),C\left(\mathrm{2,2}\right),C\left(\mathrm{2,3}\right),C\left(\mathrm{2,4}\right)\\ C\left(\mathrm{3,1}\right),C\left(\mathrm{3,2}\right),C\left(\mathrm{3,3}\right),C\left(\mathrm{3,4}\right)\\ C\left(\mathrm{4,1}\right),C\left(\mathrm{4,2}\right),C\left(\mathrm{4,3}\right),C\left(\mathrm{4,4}\right)\end{array}\right];$

(3)计算协方差矩阵C的特征值与特征向量

协方差矩阵C的特征方程是|λE-C|＝0，它的根就是C的全部特征根，对每个特征值λ₀，其次方程组(λ₀E-C)x＝0的每一个非零解都是C的对应于λ₀的特征向量；求出C的特征值为λ₁，λ₂，λ₃，λ₄，要求λ₁＞λ₂＞λ₃＞λ₄；相应的特征向量为μ₁，μ₂，μ₃，μ₄；

(4)计算图像的主分量

用相应的特征向量作为行，就得到变换矩阵A，根据主成分变换的数学原理，A是X空间的协方差矩阵C的转置矩阵(T)，即

A＝C^T

用该矩阵对每个像素进行变换，可求得图像的主分量，变换后，用4维向量表示像素的第一个元素组成了原图像的第一主分量，第二个元素组成了原图像的第二主分量，同样，第三、四个元素组成了原图像的第三、四主分量，即得到4个主成分分量PC1，PC2，PC3，PC4。

Images