基于独立成分法的遥感蚀变信息提取研究_以新疆包古图矿区为例
基于独立成分法的遥感蚀变信息提取研究
——以新疆包古图矿区为例
李明明1 ，周可法2
（1 Esri中国（北京）有限公司西安代表处，陕西 西安 700075;
2 中国科学院新疆生态与地理研究所，新疆 乌鲁木齐 840000）
摘 要：近矿围岩蚀变是矿化的1个主要特征，是矿床蚀变的1个直接标志。可以用来发现和预测各类矿产，增加找到矿床的机会。遥感图像蚀变信息的提取已经成为遥感找矿的1个重要方法，但这些对找矿有指导意义的矿化蚀变信息常常受其它地物信息的干扰，而且受遥感图像的波谱分辨率和空间分辨率的制约，往往表现的很微弱。因此，国内外学者也在不断尝试各种技术方法提取这种矿化蚀变弱信息。
通过比较ICA（独立成分分析）和PCA(主成分分析)方法的相似性与差异性，从统计分析的角度看，简单ICA和PCA一样，同属多变量数据分析的线性方法。因此结合PCA在遥感蚀变信息的应用，以ICA方法替代传统的PCA方法，应用ICA的原理，提出了1种基于ICA的遥感蚀变信息提取方法。
实验中按照PCA提取蚀变信息的准则来处理经ICA处理的成分，采用ENVI软件中提供的独立成分分析模块实行。以新疆包古图矿区为例，利用Crosta方法提取蚀变信息，通过ETM+1、ETM+3、ETM+4、ETM+5和ETM+1、ETM+4、ETM+5、ETM+7的波段组合分别进行ICA分析提取铁染蚀变和羟基蚀变信息。对代表铁染蚀变的独立成分的判断准则是：ETM+3的系数应与ETM+1、ETM+4的系数相反。选择第四成分进行铁染蚀变信息的提取。对代表羟基和碳酸根离子独立成分的判断准则是：ETM+5系数应与ETM+7、ETM+4的系数符号相反，ETM+1一般与ETM+5系数符号相同。选择第四成分进行羟基蚀变信息的提取。铁染蚀变和羟基蚀变存在于绝大多数成矿岩体中，提取这两种蚀变信息基本可以确定研究区成矿岩石的分布情况。 最后，对ETM+影像提取羟基蚀变信息和铁染蚀变信息的结果进行野外实地验证。结果表明：应用ICA法提取的蚀变信息的有效信息量要大于PCA法，伪异常信息明显要少于PCA 法。为找矿提供了更为准确的信息，减少找矿成本。结果证明了独立成分分析方法的有效性。
关键字：独立成分分析；主成分分析；蚀变信息；遥感；ENVI软件
文献标识号：A
0引言
遥感蚀变异常信息提取的方法有多种，常用或较常用的有比值法、主成分分析法、光谱角法、归一化植被指数法、混合像元分解法以及它们的混合法[1]。每种方法都有其针对性和适用的范围，针对不同的区域选择应不同的方法和组合。其中，主成分分析法提取蚀变信息是相对最为广泛的。主成分分析（Principal Component Analysis，PCA）是基于信号二阶统计特性的分析方法，由于所获各主成分之间互不相关，主成分之间信息没有重复或冗余。PCA的这一基本性质在蚀变异常信息提取中被充分利用。多光谱遥感数据通过PCA所获每一主成分常常代表一定的地质意义，且互不重复，即各主成分的地质意义有其独特性。但是，由于蚀变矿物形成的影像特征在遥感图像上往往表现得很微弱或不明显，甚至被“淹没”在主体色调中。很多弱信息包含在高阶统计量中[2]。若仅仅采用PCA算法，对图像进行处理可能会忽略很多关键信息[3]。因此，若有1种数据处理方法能够有效地消除影像间的各种相关信息，特别是高阶相关信息，则对于矿物蚀变信息的提取具有非常重要的意义。
独立成分分析（Independent Component Analysis，ICA）算法是一种基于高阶统计量的去相关多元数据处理方法，不仅能够消除多变量数据中的二阶相关关系，且能够消除数据间的高阶相关。已成功地应用于语音信号处理、通信、人脸识别、图像特征提取、神经计算和医学信号处理等众多领域。并且已经将其用于遥感图像的数据变换与分析，如，对超光谱图像进行波段选择和压缩以及光谱图像的分类和变化检测[4]。鉴于ICA方法的特点及其在应用领域取得的成果，笔者在利用多光谱影像提取矿物蚀变信息的过程中引入ICA方法，使得随机变量之间相互独立，更简洁有效地表达数据并有利于后续的信息分割等处理，最终取得了较好的效果。
1 研究方法
1.1 ICA 原理
目前，在信号处理领域有1种ICA 的方法，它能从观测信号出发，估计出已知信息量很少的源信号，而所
获得的源信号是互相独立的[5]。用数学语言描述，即，通过独立成分分析，从一组观察信号x =｛x 1,x 2,…,x m ｝
分离出一组统计独立的信号源s =｛s 1,s 2,…,s n ｝的估计值s =｛y 1,y 2,…,y n ｝，y 也是统计独立的。
假定第i 个观测信号由n 个独立成分线性混合而成：
1122,1,2,,i i i in n x a s a s a s i m =+++=L L (1)
用矢量x 表示观测变量｛x 1,x 2,…,x m ｝T ,矢量S 表示源变量｛s 1,s 2,…,s n ｝T
,A (m ×n )表示混合矩阵a ij 式(1)可表示成： 1n
j j j x a s As ===∑ (2)
统计模型式(2)称为独立成分分析模型。该模型描述了观测数据是如何由信息源s 混合生成的。源变量s 是隐藏变量，不能直接观测到，而且混合矩阵A 也是未知。所有能观测到的数据只有随机变量x ，所以必须估计出混合矩阵A 和s 。
独立成分分析模型的求解基于1个简单的假设：源变量s 是统计独立的，同时非高斯分布。在估计出混合矩
阵A 后，可以进一步计算混合矩阵A 的逆，也即分离矩阵W =A -1，从而得到独立成分s 的估计y 。
y Wx = (3)
ICA的目的就是估计出分离矩阵W 。目前独立成分分析的实现方法有很多：信息最大化法、基于极大似然估
计的ICA算法、最小化互信息(负熵最大化)法、FastICA算法等[6]。笔者采用的是遥感处理软件ENVI中提供的独
立成分分析方法来处理影像。
遥感影像所表现的信息是不同地物对某个波段反射率的1种反映，多光谱遥感影像的各波段可以认为是相互独立的的各种地物的光谱信息按一定的规律混合的结果。将ICA 作用于多光谱影像，得到的每个独立成分可
集中体现某些地物的信息，使得不同地物间的分离度增大[7]。
1.2 研究方法与技术
笔者采用ICA 方法提取波段信息中的高阶信息和独立成分，以独立成分替代主成分提取铁染和羟基的蚀变信息。经ICA 处理后的各成分相对于PCA 的各成分来说是独立的，且包含了一些基于高阶统计的信息。从统计分析的角度看，简单ICA 和PCA 一样，同属多变量数据分析的线性方法。因此，可以按照PCA 提取蚀变信息的
准则来处理经ICA 处理的成分。以TM 数据为例，通过Crosta 方法[8]说明其准则，通过TM1、TM3、TM4、TM5和
TM1、TM4、TM5、TM7的波段组合分别进行主成分分析提取铁染蚀变和羟基蚀变信息[9]。由TM1、TM3、TM4、TM5
作为输入波段组合进行正交线性变换，变换后的某个新的组分图像可能集中了铁染蚀变信息。对代表铁染蚀变的主成分的判断准则是TM3的系数应与TM1、TM4的系数相反。由TM1、TM4、TM5、TM7作为输入波段进行主成分分析。对代表羟基和碳酸根离子主成分的判断准则是TM5系数应与TM7、TM4的系数符号相反，TM1一般与TM5系数符号相同。依有关地物的波谱特征，羟基和碳酸根离子信息包含于符合这判断准则的主成分内。铁染蚀变
和羟基蚀变存在于绝大多数成矿岩体中，提取这两种蚀变信息基本可以确定研究区成矿岩石的分布情况[10]。
所以，笔者的研究方法就是通过ETM+1、ETM+3、ETM+4、ETM+5和ETM+1、ETM+4、ETM+5、ETM+7波段组合分别进行ICA提取铁染蚀变和羟基蚀变信息。为了显示ICA提取结果的优越性，同时也进行PCA处理的对比试验。
基于独立成分分析法提取矿物蚀变信息的方法主要包括：图像预处理、ICA 分析、选择有效成分、异常切割、蚀变信息处理等步骤（见图1）。所有的操作步骤都是在ENVI 软件中实现的。原始数据经过一系列的图像预处理，包括几何纠正，大气校正，掩膜去背景等；再选择波段进行独立成分分析和主成分分析，根据上述的原则进行成分的选择；将选择的成分用相同的异常切割标准处理，突出蚀变信息；并与已知的矿点资料进行叠合分析，比较ICA 和PCA 方法提取的蚀变结果。
在做异常切割时，利用( X + kσ) 确定异常下限和划分异常强度等级。X 是某一成分的统计均值代表区域背景，σ是该成分的标准差。k 取值一般是（1～3），有了这一标准，切割异常时可以减少主观任意性，并使操
作较为规范化[11]。经ICA处理后得到的各成分是相互独立的，理论上可以看成是单一的地物。但是由于大气，
传感器等影响，各成分并不只是单一的地物，所以在提取异常时还是利用了阈值切割法。
图1 技术路线图
2应用实例
以包古图斑岩铜矿的Ⅱ号、Ⅴ号岩体为研究区域，以ETM影像作为数据源，展示ICA方法提取蚀变信息的效果。
包古图铜金多金属矿区位于新疆托里县境内，其地理坐标为东经84°15′00″～84°22′30″，北纬45°20′00″～45°35′00″。属中低山区，海拔高度600 m～1100 m，相对高差一般在150 m～500 m，地形总体较平缓，交通便利。该区矿产丰富，矿化蚀变信息较强，基岩裸露多，植被覆盖较少，很适合利用遥感技术方法来探测未知区的矿化信息。本研究采用的是2002年9月21日的ETM+影像，行号为145，列号为028。
该矿区围岩蚀变广泛，Ⅱ号斑岩体的蚀变矿物主要是绿泥石、方解石、埃洛石、伊利石。Ⅴ号斑岩体围岩蚀变有伊利石化、绿泥石化、绿帘石化、绢云母化、钾长石化、黑云母化及水白云母化等。钾化带中的蚀变矿物组合为石英、钾长石化、黑云母、绢云母。石英绢云母化带的蚀变矿物组合为石英、绢云母、水白云母、黄铁矿等。青盘岩化带的蚀变矿物组合为绿泥石、黝帘石、钠长石，局部还有黑云母及绢云母的叠加。矿化在花岗闪长斑岩和凝灰质粉砂岩中都有不同程度的发育[12]。
2.1铁染信息的提取
该区围岩蚀变较强，根据含铁离子矿物的特征波谱，对ETM+1、3、4、5波段进行独立成分分析，各特征向量见表1。为了比较ICA和PCA提取蚀变信息的区别，对ETM+1、3、4、5波段也做PCA处理，各成分特征向量见表2。
表1 ETM+1、ETM+3、ETM+4、ETM+5独立成分分析特征向量载荷表
ICs ETM+1 ETM+3 ETM+4 ETM+5 特征值
IC1 0.200 226 0.432 757 0.578 697 0.661 620 59.622 537
IC2 0.216 342 0.344 090 0.523 758 -0.748 650 4.812 713
IC3 -0.513 766 -0.597 408 0.615 702 0.007 704 2.727 232
IC4 0.805 698 -0.580 885 0.108 162 0.041 516 0.731 692
表2 ETM+1、ETM+3、ETM+4、ETM+5 主成分分析特征向量载荷表
PCs ETM+1 ETM+3 ETM+4 ETM+5 特征值
PC1 0.200 226 0.432 757 0.578 697 0.661 620 59.622 537
PC2 0.216 342 0.344 090 0.523 758 -0.748 650 4.812 713
PC3 -0.513 766 -0.597 408 0.615 702 0.007 704 2.727 232
PC4 0.805 698 -0.580 885 0.108 162 0.041 516 0.731 692
各独立成分和主成分的特征向量是一样的，这说明ICA和PCA在统计分析上是一致的。同时，也表明ICA可以和PCA一样可用于提取蚀变信息。经ICA和PCA 2种方法分别处理后各成分的特征向量载荷是一样的（见表1、表2），但并不意味着ICA和PCA的统计量是一样的。例如，成分PC4的直方图是光滑的曲线，而成分IC4的直方图是锯齿形的、有多个波峰和波谷（见图2），说明成分IC4里含有多类信息，是经过ICA的高阶统计分析获取的。这些微弱的蚀变信息更多的更完全的体现在成分IC4中。结合已知地质资料设定合理的阈值，可以将蚀变信息和其他噪声区别出来，从而将直方图中表示的不同类的信息提取出来。
首先根据表1和表2可以看出IC4/PC4成分主要反映了ETM+1和ETM+3的信息，正符合铁染信息在ETM+1波段具有强吸收，在ETM+3具有强反射的特点。所以，将IC4/PC4作为反映铁染信息的成分。此外，根据一般的原则IC4/PC4的特征值最小，基本上是噪声，由于围岩的蚀变信息在遥感影像上是极其微弱的，所以IC4/PC4是极有可能混有蚀变信息和噪声的成分。为了比较这2种方法提取蚀变信息的差异性，必须保证2种方法提取的矿物蚀变异常信息的标准一致，将得到的IC4和PC4 进行同一水平切割。根据研究区已有的矿点信息和蚀变带来确定阈值，采用
相同的切割标准（X+1.5σ）对其直方图进行切割，然后对切割后的二值图像进行窗口为3×3的中值滤波，最终得到与地质背景相符合的蚀变异常（见图3、图4）。
图2 ETM1、3、4、5经ICA、PCA处理后的IC4、PC4成分的直方图
经ICA和PCA提取的部分铁染蚀变信息能与矿点很好的重合，2种方法提取的部分蚀变信息是交织在一起的，说明了2种方法的共性。但是，很明显的看到在同一切割标准下，PCA方法提取的蚀变信息中包含很多的伪异常。经实地验证，红色标记大部分的带状信息都是植被。在Ⅱ号矿体处，经ICA提取的蚀变信息几乎均落在该矿体的已知蚀变带内。经PCA提取的蚀变信息有部分落在已知蚀变带内，但同ICA法比较也有较多的信息落在已知蚀变带的范围之外，且经实地验证，落在已知蚀变带外的信息都是伪异常。Ⅱ号矿体北面和南面的红色区域证明是沿河的植被（见图4）。在Ⅴ号矿体处，2种方法提取的蚀变信息落在已知蚀变带内的信息量基本一样的，这说明在该矿体处两种方法提取的效果相当。从整个研究区来看，虽然经ICA提取的蚀变信息量比PCA要少，但提取的大部分蚀变区域对应的岩石都有蚀变。而且在已知蚀变带范围内，有效地蚀变信息量要多于经PCA方法提取的。
图3 经ICA提取的铁染蚀变信息图
图4 经PCA提取的铁染蚀变信息图
2.2羟基蚀变信息提取
该区泥化蚀变较强，含羟基的蚀变矿物主要有伊利石、绢云母、绿泥石、黑云母等。这些蚀变含有OH-、CO32-离子或离子基团，它们在ETM+5波段均具有高反射，而在ETM+7波段具有强吸收的特征[9]。因此，采用独立成分分析法，对ETM+1、4、5、7波段做ICA处理来提取OH-和CO32-泥化信息，各特征向量见表3。由于独立成分和主成分分析的统计特征向量值是趋于一致的，所以只给出独立成分分析的特征向量的载荷表。
表3 ETM1、ETM4、ETM5、ETM7独立成分分析特征向量载荷表
ICs ETM1 ETM4 ETM5 ETM7 特征值
IC1 0.174 314 0.497 060 0.628 248 0.572 581 70.822 049
IC2 0.122 302 0.811 381 -0.182 095 -0.541 798 7.227 081
IC3 0.877 439 -0.064 532 -0.410 618 0.239 434 1.838 025
IC4 -0.429833 0.300 712 -0.635 247 0.566 813 0.731 899
在表3中，由于5和7波段的贡献系数最大，说明它们对IC4的贡献最大，也就是说IC4的信息主要来源于ETM+5、7波段，且ETM+5为负值，ETM+7为正值，具有相反的贡献标志。所以，IC4主要反映5和7波段信息，这与蚀变阴离子在5波段具有高反射，在7波段具有高吸收的特点相符合，同时符合判断准则。因此，该成分可以作为泥化蚀变信息的指示分量。根据研究区已有的矿点信息和蚀变带来确定阈值对选择的IC4和PC4 进行同一水平切割。本次试验采用（X+2σ）对IC4和PC4的直方图进行切割，然后对切割后的二值图像进行窗口为3×3 的中值滤波，最终得到与地质背景相符合的蚀变异常。
将提取的蚀变信息结果图与已知的矿点和蚀变带进行叠加分析，图5、图6中可以看到ICA和PCA提取的部分蚀变信息能很好的与II号和V号矿体重合；在图7中可以很清楚的发现，在II号矿体处，经ICA法提取的有效蚀变信息量要高于经PCA法；在V号矿体处，两种方法提取的信息落在已知蚀变带内的信息量是一样的，但是在该矿体西面，出现了明显的红色伪异常信息。这说明了在同一标准下，经ICA法提取的有效蚀变信息量是大于PCA法的，而且相对PCA法，ICA法提取的伪异常是较少的。
图5 经ICA提取的羟基蚀变信息图
图6 经PCA提取的羟基蚀变信息图
图7 II号、V号矿体局部放大图
3结论与讨论
独立成分分析利用了较主成分分析更高的高阶统计信息，获得独立的成分，因而，能更加全面揭示多光谱数据间的本质结构。笔者根据ICA和PCA在原理和方法的相似性和差异性出发，利用ICA方法提取矿物蚀变信息，经实验取得了较好的效果，在同一标准下比传统的PCA方法效果要好。这一方法在类似区域开展遥感蚀变信息提取工作具有一定的指导意义。
将此次经ICA提取蚀变信息的结果与该区矿产资料进行叠合分析，蚀变信息与已知矿化点吻合良好，发现区内已知的矿点和矿化蚀变带均在所提取的蚀变信息异常内。
在铁染蚀变的提取中，PCA方法提取信息量明显高于ICA方法，但是经野外实地调查，验证了PCA方法提取的铁染蚀变信息有很多是伪异常；ICA方法提取的信息有效地减少了伪蚀变异常，从而为找矿提供了更为准确的信息，减少找矿成本（见图3、图4）。
在羟基异常的提取中ICA方法提取的有效信息量明显高于PCA方法，这有利于对围岩蚀变弱信息的确定，快速准确的定位信息量最大的区域。
鉴于工作区属于浅覆盖负地形地貌, 所以异常没有大片出露, 但所检查的几处异常点在覆盖之下都有强弱不同的岩石蚀变。
因此，这表明应用ICA方法所提取的遥感蚀变信息比常规的PCA方法要有效、可靠。
但将ICA方法用于提取遥感蚀变信息，可能存在以下问题：
计算独立成分的方法不同及选择迭代次数不同，所获取的蚀变信息量可能不同。针对不同的区域选择不同的计算方法来获取适应该区域的最佳计算独立成分的方法。
在对选择的成分进行分割时，阈值的选择具有一定的局限性，根据已知的地质资料并不能完全满足要求，而且阈值的确定对提取的结果有很大的影响，如何实现快速确定阈值有待进一步的研究。
以上的这些问题需要进一步的研究改进。
参考文献
[1] 荆 凤,陈建平.矿化蚀变信息的遥感提取方法综述[J].遥感信息,2005(2):62-65,57.
[2] 万 宁,吴 飞.基于ICA 的全色影像和多光谱影像融合算法[J].计算机工程,2006,32((7):218-220.
[3] 钟家强,王润生.基于独立成分分析的多时相遥感图像变化检测[J].电子与信息学报,2006,28(6):994-998.
[4] CHANG C I, CHANG S S,J SMITH,et al.Linear spectral random mixture analysis for hyperspectral imagery[J].IEEE
trans.on Geosci.Remote Sensing,2002,40(2):375-392.
[5] C JUTTEN,J HERAUL.Independent component analysis verus principal component analysis.in Proc.Europ.Singal
Processing Conf[J].EUSIPC0,1998,88:643-646.
[6] A HYVRINEN.Fast and robust fixed-point algorithms for independent component analysis[J].IEEE Trans.on Neural
Networks，,1999,10(3):626-634.
[7] 边 境,周庆利.一种多光谱图像和全色图像融合算法[J].微电子学与计算机,2005,22(6):5-7,11.
[8] M H TANGSTANI,F MOORE.Iron oxide and hydrotyl enhancement wsing the crosta method：a case study from the Iagtos
Belt,Fars Province,Iran[J].Infornational Journal of Applied Earth Observation and Geoinformation,2000,2(2):140 -146.
[9] 韩 玲,谢秋昌,利用遥感技术对新疆西天山地区矿化蚀变信息的提取[J].遥感信息,2007(6):49-51.
[10] 张玉君,曾朝铭,陈 薇.ETM+(TM)蚀变遥感异常提取方法研究与应用——方法选择和技术流程[J].国土资源遥感,2003(2):
44-49.
[11] 张楠楠,张晓帆,周可法,等.干旱区TM 图像蚀变信息提取方法研究[J].干旱区地理,2006,29(3):398-402.
[12] 耿新霞,建 民,张玉君,等.新疆西准噶尔包古图斑岩铜矿蚀变岩石的光谱特征及其意义[J].现代地质,2008,22(1):116-
122.
基金项目：中国科学院知识创新工程重要方向项目（kzcx2-yw-107）,国家自然科学基金项目（40601103）
第一作者简介：李明明，1983年生，男，湖北潜江人，硕士研究生，主要从事地理信息系统和遥感信息提取技术与应用研究。