空间U统计量法在遥感蚀变信息提取中的应用研究

胡 波，朱谷昌，张远飞，肖 燃

空间U统计量法在遥感蚀变信息提取中的应用研究

胡 波1，2，朱谷昌1，2，张远飞2，肖 燃3

(1.中南大学，长沙 410083;2.有色金属矿产地质调查中心，北京 100012;3.江西应用工程职业学院，萍乡 337000)

鉴于传统蚀变异常信息提取方法的局限性，试图在二维光谱特征空间(即散点图)中研究遥感蚀变信息的提取方法。在二维散点图中，两个波段灰度值的联合分布体现了各向异性的特点，聚类常呈类似椭圆形分布。对散点图中点的频数使用空间U统计量法逐个获得地物聚类的椭圆参数，将每个椭圆内的点映射到遥感图像上，经过目视解译，最终获得蚀变异常信息的空间分布。以青海省巴音山地区蚀变信息提取为例，对空间U统计量方法进行了应用，并将提取的异常区域与已获得的各种数据进行对比，发现铁化和泥化区域吻合度较高，取得了比主成分分析更好的蚀变信息提取效果。

蚀变异常信息提取;U统计量法;聚类各向异性;二维散点图;异常聚类

0 引言

当前，用于遥感蚀变异常信息提取的方法主要有基于主成分分析法和波段比值法两种。由于这两种方法采用的是相对固定的处理流程和分析手段，往往不能够因地制宜地提取遥感蚀变异常信息。张玉君［1］等人对几种传统蚀变信息提取方法做了总结，认为应用比值法得到的结果不可避免地会出现难以解释的假异常，故提出采用以主成分分析为主，辅以比值法和光谱角法的“去干扰异常主分量门限化技术流程”;张远飞［2］等人的研究发现，基于掩模+主成分分析+比值法的方法适用环境存在局限性。从光谱特征空间来说，传统的方法对于偏离出背景的蚀变信息捕捉效果较好。

笔者尝试在二维光谱特征空间中研究遥感蚀变异常信息的提取。二维散点图是二维光谱特征空间的一种平面表现形式，它所反映的是2个波段灰度值的联合分布。同类地物无论在空间上分布如何，其在散点图当中总是形成聚类，并且呈现出近似椭圆的形态;而不同地物组合在一起则形成椭圆的套叠结构。找到蚀变异常信息所构成的聚类椭圆，然后将其映射到遥感图像上，就可以获得蚀变异常信息的空间分布。

空间U统计量法是一种用于多重分形的局部奇异性分析方法，它能够获取局部位置的各向异性参数［3，4］，并且可以有效减少地物聚类椭圆之间交叠区域被误判的概率。在散点图中使用该方法可以获得各个聚类椭圆的参数，依次找到并剔除背景和各种干扰，最终提取到蚀变异常信息。

本文以青海省都兰县巴音山地区为例，采用空间U统计量法提取了铁化和泥化蚀变异常信息，并将其与该地区的已知数据进行对比，验证方法的应用效果。

1 光谱特征空间遥感数据的结构特点

1.1 光谱特征空间的相关定义

光谱特征空间是指在一幅多光谱图像中假设有n个波段，则每一个像元在各个波段上的灰度值将构成一个向量，用 X=(x1，x2，…，xn)T表示。包含所有X的n维空间称为光谱特征空间［5］。本文重点研究的是二维光谱特征空间，如TM1与TM3的组合，那么每个像元矢量即为X=(xTM1，yTM3)T。以该像元在TM1波段上的灰度值作为x坐标，在TM3波段上的灰度值作为y坐标，包含所有X的空间就是本文所说的二维光谱特征空间(即二维散点图，文中简称散点图)。散点图中(x，y)处的值代表了遥感图像中TM1=x且TM3=y这种联合分布的频数。把频数作为z轴，可以在三维空间上观察地物聚类的形态和分布特征，文中称立体散点图。

1.2 背景、干扰与蚀变异常的定义

文献［6］基于光谱特征空间给出了背景、干扰与蚀变异常的定义。背景是指在多波段遥感图像数据点集空间中的超椭球体(实际上为近似超平面体)，而游离在超椭球体之外的“小类”点集则为干扰或蚀变异常目标。从统计角度来看，遥感图像是由散点图中高频数的背景和相对低频数的异常组成。背景是图像上频数最高的区域，在立体散点图当中呈“主峰”形态;异常是相对于背景定义的，在“主峰”周边有时会出现一些相对独立的“次峰”形态，这就是异常。本次研究中蚀变异常信息提取的目标就是首先在光谱特征空间中定位蚀变异常，然后使用一定的技术方法获得其空间分布。

1.3 光谱特征空间遥感数据的几何特征

在二维散点图中，地物灰度的联合分布会在不同方向上产生变形，形成的聚类体现出各向异性的特点，往往近似一个椭圆。依照多元高斯分布等值线描述［7］:多元高斯密度的等值线是中心在均值向量m处的超椭球平面。该超椭球面的主轴平行于特征向量e1，而特征向量a是相应的方差，如图1所示。

图1 多元高斯密度等值线Fig.1 Isoline of multi－ Gaussian density

假设地物灰度的概率密度服从一定的正态分布，在2个波段组合中，这种地物的灰度概率密度服从二元高斯分布。散点图上，每类地物会形成椭圆形的聚类。

众多地区的散点图上都出现了椭圆的套叠结构，频数最高区域对应的是背景，而周边的频数相对较低，变化较缓的椭圆对应的是异常。假设每一类地物的灰度服从正态分布，那么每个波段的灰度直方图是叠加之后的多种地物灰度的分布。可以作出以下推断:

(1)在散点图当中，各类地物分别服从二元高斯分布，但由于受到地物混合影响，其联合分布形态发生变化，但仍近似于椭圆形;

(2)遥感图像是多种地物的整体映射，故在散点图中应表现为椭圆的套叠。

基于二维光谱特征空间中遥感信息呈现的近似椭圆形的空间分布特征，使用空间U统计量法获得各向异性参数，并且定位蚀变异常的聚类椭圆，最终便可提取到蚀变异常信息。

2 空间U统计量法

U统计量的构造为

式中，z0表示在分类范围内异常总体和背景总体的分类阈值，通常采用局部的均值。一般情况下，当 U［B(x，y)(r，β，θ)］＞ 0 时，认为样本属于异常总体;反之可以判定样本属于背景总体。|U|越大，则表明分类效果越好;而|U|取值接近0，则意味着样本内存有混合现象。

最佳的窗口参数应当是U值取得最大时所对应的 U*［B(x，y)(r，β，θ)］，即

此时的椭圆参数(r，β，θ)就是所求。随着窗口半径的加大，U*［B(x，y)(r，β，θ)］的值应该不断增大才有意义。如果取值减小后又继续增大，就意味着窗口已经跨越了边界，应该在值出现减小时就终止运算，所以，有必要用生成的U*－r曲线图作为辅助(图2)。

图2 U*－r曲线示例Fig.2 Instance of U*－ r curve

图2上曲线1表示异常总体，在|U*|值减小前也就是rA处取得极值;曲线2表示背景总体，在|U*|值减小前也就是rB处取得极值，当|U*|接近0时的r'，表明背景和异常有混杂现象。

根据空间U统计量的构造方法可以看出，对于一个样本，它的期望μ没有改变，而方差σ2［U］=σ2/n。故，当判断一个样本是属于背景总体还是异常总体时，显著缩小了出现判定错误的区间(图3)。

图3 用样本值和U作为统计量的分布函数Fig.3 Distribution function by statistics of sample value and U value

图3 中α表示A类样本被判定为B类的概率;β表示B类样本被判定为A类的概率。不难看出，采用U统计量之后，出现两类错误的概率显著减小了。

3 应用实例

以青海省都兰县巴音山铜、铅、锌、银多金属矿集区为实验区，该区位于柴北缘火山弧裂陷构造带、东昆仑北构造带及鄂拉山构造带三者交汇处，那里岩浆活动强烈且岩类复杂，成矿物质来源广泛，成矿作用类型多种，矿产种类较多。

该区出露的地层主要有古元古界沙柳河群(Pt1dk)片岩、麻岩，古生界上奥陶统滩涧山群(O3tn)以基性、中基性、酸性火山岩及其碎屑岩类为特点的细碧角斑岩类，间夹长英质碎屑岩、硅质岩、大理岩透镜体的岩性组合，变质后一般为绿片岩和绢云母石英片岩、千枚岩及大理岩等。其上不整合覆盖有泥盆系上统砂岩、砾岩和千枚岩。

本区构造以柴北缘火山弧裂陷—裂谷构造带及其后期继承性近东西向构造线为主导格局。区内构造较为复杂，为矿床定位造就了良好条件。除下古生界广泛发育双峰式火山岩系外，加里东晚期和海西晚期均有花岗岩体沿东西向断裂或层间断裂带侵入，为成矿物质流体的活动提供了热动力条件。

综合分析，阿尔茨托山东段的巴音山一带是重要的多金属成矿带。该地区有着详实的野外考察资料和丰富的工作总结，拥有充分的地质、物探、化探及遥感等数据储备，因而选择该地区作为已知地区进行新蚀变信息提取方法的实验和评价。

传统方法使用 TM1、TM3、TM4、TM5等波段组合进行主成分提取，并通过判定确定PC4是代表铁化蚀变的主分量，利用最优彩色密度分割的方法，最终取得铁化蚀变的区域。通过与各种数据比较，发现提取到的铁化蚀变区域与巴音山顶、东山及乌龙滩南部等几个大型矿区吻合良好，但是在沙柳河东岸并处于山脊线西侧的几个小型矿点，与之对应的蚀变却全部出现在山脊线东侧，偏离较大。同样，使用TM1、TM4、TM5、TM7等波段提取泥化蚀变信息，其聚类椭圆特点与铁化信息的类似。

［8］，笔者对于巴音山地区采用TM4与TM3波段组合提取铁化异常信息，用TM7与TM5波段组合提取泥化异常信息。

从TM4和TM3的立体散点图(图4)不难看出，此波段组合下的联合灰度分布有明显的聚类特征。图像中的“主峰”对应背景信息，其周边相对低矮的“次峰”对应异常信息。

图4 TM3与TM4波段组合立体散点图Fig.4 3D scatter plot by TM3 and TM4

图5 是在TM4和TM3散点图中提取各类地物的流程，其中蓝—黄—绿—红—品红的过渡表示频数由低至高。

通过使用空间U统计量法，首先可以得到聚类椭圆I。这是图像上频数最高的部分(即背景)，为避免其对之后提取聚类的影响，故先将其剔除。通过剔除局部最高的峰值，之前被高峰所掩盖的较低峰的区域变成最高峰，变为红色显示出来。当剔除背景椭圆I之后，图像中另外两个椭圆逐渐显露出来。如此通过不断地“捕捉椭圆—去除椭圆—捕捉椭圆”，层层剥离背景和干扰，最终使得遥感蚀变异常这类弱信息“水落石出”。

图5 巴音山TM3与TM4组合散点图中的聚类提取流程Fig.5 Extraction process of clusters in scatter plot by TM3 and TM4 in Bayinshan

经过上述提取流程，最终确定了6个聚类椭圆。图6是散点图上各个椭圆的位置。表1列出了各个椭圆的参数。将各椭圆映射到图像上进行目视解译，对聚类椭圆 I、III、IV、VI分别定性判断为背景、浅层植被、裸露河床、深厚植被。椭圆II所对应的主要是沿沙柳河、乌龙滩东侧地区以及巴音山东侧河流沿岸地区，该地区径流密布(基本干涸);椭圆V对应沿巴音山的山脊线两侧分布的区域、东山顶部以及哈莉哈德山顶部区域。

图6 TM3与TM4提取的6个椭圆Fig.6 6 ovals extracted from TM3 and TM4

表1 TM3与TM4散点图中的聚类椭圆参数Tab.1 Cluster ovals’parameters of scatter plot by TM3 and TM4

图7中黄色曲线框中的异常区域对应椭圆V，其余对应椭圆II。这两个椭圆所对应地区共同特点是矿产丰富且有多处已探明的矿床。从已知矿点分布来看，椭圆II和椭圆V同属蚀变异常的信息聚类，但却属于不同的聚类椭圆。据考察资料显示，这些地区蚀变均以硅化、绢云母化、绿泥石化及褐铁矿化等为主，蚀变矿物一致，所以蚀变异常出现两个聚类不是由于蚀变矿物不同造成的，而是受河流泥沙类物质影响。研究区域河流基本干涸，河床裸露，河道中的泥沙类物质反射率较高，使得裸露河床出现在散点图右上方(椭圆IV)。野外记录显示，河道周边确实存在一定量的蚀变物质并且与泥沙类物质混合，这就造成沿沙柳河及乌龙滩局部铁化蚀变异常区域反射率较高。而没有受到泥沙类物质影响的巴音山顶部、东山顶部和哈莉哈德山顶部反射率相对较低。从整体形态来看，椭圆II和椭圆V同在散点图对角线上侧，符合Fe3+在TM3波段反射率高，在TM4波段反射率低的物理特征。

另外，从图6中还可以看出，有两个植被的聚类，浅层植被对应植被覆盖稀疏地区，深层植被代表植被茂盛地区。浅层植被地区在散点图中的位置连接着深层植被与背景，充分反映了其植物、岩石及土壤混杂的特点。而对于深层植被地区，因为植被十分茂盛几乎遮蔽了其余地物的反射光，故其光谱特征较为独立，距离背景椭圆I较远。

将椭圆II、椭圆V映射到遥感图像上便获得了铁化蚀变异常区域。对比该地区物探和化探数据，发现提取结果与之高度对应，并且与野外记录相符。图7是提取到的铁化蚀变异常与已知矿点分布叠加图。图上蓝色—黄色—绿色—红色—品红过渡色表示蚀变异常由弱至强。经统计，在21个已知矿点当中，有19个毗邻或者存在于铁化异常区，命中率达90.4%。距离铁化蚀变区域较远的有2个，分别是吉给申沟铅矿床(野外勘察显示该处矿体规模较小，为共生矿产)和藏碑沟铅、锌、银多金属矿点(未见详细记录)。

图7 铁化蚀变与已知矿点分布叠加图Fig.7 Stacking chart of ironed alteration and known deposits

利用TM7与TM5影像组合的散点图提取泥化蚀变异常过程与上述提取铁化信息的类似。泥化信息也出现了两个聚类椭圆，与铁化蚀变信息的聚类椭圆呈同样的特点，此处不再赘述。

将蚀变异常信息映射在遥感图像上，经统计，在21个矿点中，有20个在泥化异常区内，命中率达95.2%，仅吉给申沟铅矿床未出现在蚀变区域内。

叠加两类异常，发现在矿点周围基本都有异常分布，从而肯定了提取结果的正确性。在沙柳河东岸、巴音山北部、哈莉哈德山及东山南部等地区(图7中红色框范围)，泥化和铁化均有广泛分布，但是矿体分布有限，需要进一步工作查明原因。而乌龙滩北部等其余蓝框范围内虽然有大片蚀变分布，但这是由于河流淤积造成的蚀变矿物聚集，考察意义不大。

对比主成分分析方法的提取效果，本文方法提取结果更加准确，尤其是在沙柳河东岸的几个矿点，异常与矿点吻合度较高。

4 结论

(1)从光谱特征空间分析入手，地物分布特征表现的更加直观。在散点图中，原本空间上分散的地物能够形成聚类，并且在立体散点图中形成局部的“高峰”。通过观察二维散点图和立体散点图的图形特征，蚀变异常信息提取可以做到有的放矢。

(2)蚀变异常信息在散点图上可能出现不止一个聚类，这取决于蚀变矿物的光谱以及诸如矿物粒度、地物混合情况、含水量及阴影等综合因素的影响。有必要加强对研究地区地理现象的整理和分析，以避免蚀变异常提取过程中错提、漏提现象。

(3)应用空间U统计量法提取蚀变异常信息是一种新尝试，还要进一步加强光谱混合机制和遥感数据空间结构的研究，加深对地物聚类规律的把握，以提升蚀变异常信息提取的可信度。

参考文献:

［1］ 张玉君，曾朝铭，陈 薇.ETM+(TM)蚀变遥感异常提取方法研究与应用——方法选择和技术流程［J］.国土资源遥感，2003(2):44－50.

［2］ 张远飞，吴健生.基于遥感图像提取矿化蚀变信息［J］.有色金属矿产与勘查，1999，8(6):604 －606.

［3］ Cheng Q M，Agterberg F P，Bonham － Carter G F.A Spatial Analysis Method for Geochemical Anomaly Separation［J］.Journal of Geochemical Exploration，1996，56(3):183 －195.

［4］ 陈志军.多重分形局部奇异性分析方法及其在矿产资源信息提取中的应用［D］.武汉:中国地质大学，2007.

［5］ 中国知网.CNKI概念知识元库［DB/OL］.http//define.cnki.net/，2006.

［6］ 张远飞，吴德文，朱谷昌，等.遥感蚀变信息检测中背景与干扰问题的研究［J］.国土资源遥感，2008(2):22－26.

［7］ Therrien C W.Discrete Random Signals and Statistical Signal Processing［M］//周宗潭，董国华，徐馨，等.独立成分分析.北京:电子工业出版社，2007:24.

［8］ 张玉君，杨建民，陈 薇.ETM+(TM)蚀变遥感异常提取方法研究与应用——地质依据和波谱前提［J］.国土资源遥感，2002(4):30－36.

The Application of Spatial U－static Method to the Extraction of Alteration Anomalies

HU Bo1，2，ZHU Gu － chang1，2，ZHANG Yuan － fei2，XIAO Ran3
(1.Central South University，Changsha 410083，China;2.China Non － ferrous Metals Resource Geological Survey，Beijing 100012，China;3.Jiangxi Application Engineering Vocational College，Pingxiang 337000，China)

The authors tried to extract alteration anomalies in spectral characteristic space(scatter plot)in view of the limitation of the traditional methods.The scatter plot takes on an anisotropic feature in associated distribution of RS data’s grey scale.The distribution is usually combined by oval clusters.Parameters of oval clusters are acquired sequentially by applying U －Statistic method in the frequency of the scatter plot.Through mapping the points inside the oval into RS image and interpreting visually，the spatial distribution of alteration anomalies is obtained eventually.In this paper，this new method was described with the instance of Bayinshan area in Qinghai province.Other data acquired were also comparatively studied，and it is found that the anomalies of ferruginization and argillation are consistent well with each other.This new method has a better performance than PCA in the study area.

Alteration information extraction;U－static method;Anisotropy of cluster;2D scatter plot;Anomaly cluster

TP 753

A

1001－070X(2011)03－0071－06

2010－11－30;

2011－02－05

胡 波(1986－)，男，硕士研究生，地图学与地理信息系统专业，主要从事GIS和RS的应用研究。

(责任编辑:刁淑娟)