Question

用ETM数据提取“蚀变遥感异常”引人关注的几个技术问题

Answer 1

张玉君

（中国国上资源航空物探遥感中心，北京；返聘于中国地质科学院矿产资源研究所，北京）

摘要：本文通过模型试验讨论提取遥感异常中的两个重要问题：主分量分析的灵敏度极限（或检出限）和植被干扰容限度；还介绍为改进邻景蚀变异常的可比性所做的准归一化技术，此技术在一百多景ETM蚀变遥感异常提取中得到成功应用。

关键词：主分量分析（PCA）灵敏度极限（或检出限）；植被干扰容限度；相邻景蚀变异常可比性；准归一化

1 主分量分析（PCA）的灵敏度极限（或检出限）

当用PCA提取遥感异常时，无论是羟基异常（用TM1、TM4、TM5、TM7）或是铁染异常（用TM1、TM3、TM4、TM5）都常常出现为第4主分量（张玉君等，1998,2002,2003；李昌国等，1997），其本征值都是最小的，也就是说蚀变遥感信息仅占全景信息量的很小份额。那么应考虑一个问题：蚀变信息量少到何种程度将被PCA丢失，即PCA的检测灵敏度下限（或检出限）是多少？

为了探讨此问题，设计了如下模型实验：建四维图像（2550×2000×4），分割为（1275×1000×4）四个子图，利用在柳沟峡获得的TM图像采样数据，依左上、右上、左下、右下的次序输入金矿蚀变岩、植被、辉长岩、冰的TM1、4、5、7波段数据，并令其按列方向从0到相应波段DN值递变，递变周期为255列。若以（255×1×4）为一个单元，则全图有20000个单元，左上角的蚀变岩有5000个单元，经PCA提取，门限用3σ，得出蚀变岩占5000/20000即1/4时的PCA结果。然后对蚀变岩区依次用辉长岩取代，蚀变岩区分别保留500个单元（1275×100×4）、10个单元（1275×2×4）、5个单元（1275×1×4）、3个单元（765×1×4）和1个单元（255×1×4）；对每一次改变都进行了 PCA，门限取3σ。门限值的确定以不出现干扰（来自冰或植被）为原则，当门限降低出现干扰时意味着来自蚀变岩的信息已减弱，以致淹没于干扰中。六次PCA检出限模型实验结果列入表1：

表1　PCA检出限模型实验结果简表

续表

根据检出限模型实验结果，可以认为用PCA提取蚀变信息的检出限优于两万分之一或十万分之五。这意味着在面积为两万平方公里的图幅上，只要蚀变岩（有强有弱）的总面积不小于一平方公里，即可被PCA提取出来。此检出限模型实验结果具参考意义。

2　提取蚀变遥感异常受高植被干扰的容限度

植被是矿致遥感异常提取过程中的常见干扰因素。为了探讨植被干扰的容限度问题，设计和进行了两个模型实验。实验结果都表明，当混合像素中植被成分达到或超过50%的时候，蚀变异常就很难提取了，这一结果与澳大利亚飞机上所做测量是一致的。

为了初探高植被区的异常提取效果，并应云南地调院要求，对云南中甸高山植被茂密地区的13240景做了异常提取，还分别以普朗斑岩铜矿和红山铁矿为参考样板区，用光谱角法（SAM）对异常做了优选；经云南地调院综合地、物、化、遥选点和地面槽探查证，发现了地苏嘎Cu矿点，地表出露两个矿化斑岩体（宽度分别为6米和11米），当年即投入了钻探，初步含Cu品位为0.1～0.5%。

证实戈壁荒漠区所取得的方法技术在高植被区的应用虽受到相当限制，但仍能发挥一定作用。

3　相邻景蚀变遥感异常可比性的改进

大面积提取蚀变遥感异常工作常常需要多景图像的拼接，其中包括分层异常的拼接，日益迫切地感到不同时相的同景ETM数据或毗邻不同景ETM数据所获结果的可比性晃亟待改善。笔者进行了五项校正研究：大气径辐射校正、太阳高度角校正、日地距离校正、大气层上照度较正及增益校正，通过这五项校正明显改进了相邻景及不同时相同景信息提取结果的可比性。

在细致进行这五项校正时，避开了难以做到针对与气溶胶的密度及水蒸气的浓度相关的大气精校正，故称之为原始数据的准归一化处理（未达到严格意义上的归一化）。原始数据经准归一后，获同一比例尺的视反射率值（每一DN代表0.2%或0.25%视反射率）。

3.1径辐射校正

笔者试用了两种方法即波段相关分析法和直方图最小值法进行径辐射校正，结果对比列入表2。波段相关分析法只有当感兴趣区内仅存在一种岩性时才比较准确，故选用直方图最小值法。

表2　径辐射校正对比表

3.2　日地距离校正

地球以每秒29.79公里的平均速度沿椭圆形轨道绕太阳公转，此椭圆轨道的扁率（也称偏心率）为0.0167，长半径为149597870km，与地球赤道面相交成23°26′的角度。日地距离d定义为日心到地心的距离，取日地平均距离为天文单位距离（AU），但此距离也是变化的，1968～1983年地球与太阳的平均距离为149600000km,1984年以后为149597870km。每年约在1月4日为近日点（perihelion），约在7月4日为远日点（aphelion），近日点距太阳147100000km，远日点距太阳152100000km。

为了研究日地距离影响，笔者从天文志获得1998～2003年六年内日地距离数值表，每半小时一个点。太阳辐照度与d²成反比；过去认为日地距离变化带来的太阳辐照度变化约为5%，一般无须校正；为了更准确地作评价，特计算了1998～2003年六年日地距离变化带来的太阳辐照度变化平均为6.68%，见表3，由于此值达到近7%，故认为进行日地距离校正是必要的。

表3　日地距离最大校正系数计算表

Landsat 7的轨道为太阳同步轨道（朱述龙等，2000），轨道倾角为98.2°，周期为98.9′，过赤道时刻为地方平均时上午10时，1/4周期为24.725′，我们西部任务区分布在北纬27°～北纬44°之间，过境时间约为地方时上午9时50分左右，故可取地方时上午10时。日地距离表中所给时间为世界时，或称格林尼治时（UT），即0时区时间（经度0°子午线处的民用时），表示为GTM。为了查找某一景图的日地距离，首先确定其经度属哪个时区，以时区标准经度两侧各7.5°为界，然后计算当地上午10时对应的格林尼治时，用此格林尼治时即可查到所要求的日地距离。例如，某景ETM图像位于东经90°＋7.5°范围内，它距格林尼治有6个时区，故应查表中该口上午（10-6）=4时的日地距离d。

3.3　增益校正

在美国的Landsat 7计划中，为保证对地面数据的正常获取及对数据的充分利用，Landsat 7对ETM^＋传感器的信号处理部分进行了重新设计，使其可以在两种状态下工作，即高增益状态和低增益状态。在图像上表现为像元亮度的变化。对于1、2、3、4、5、7波段增益的改变分为三个组（林友明，2003）：波段1、2、3一同变化，波段4单独设置，波段5、7一同调整。传感器增益设置的规律是，按季节在某一纬度范围内对传感器的增益进行设置。当一景图某一波段增益设置发生变化时，在图像上该波段的亮度值便发牛阶跃性变化。据悉此变化多发生在一景图像的下四分之一处，北京地面站对有增益改变的数据进行了处理，将增益改变后的像元亮度调整到增益改变前的水平上，因此用户不会看到有阶跃性变化的图像。

但是当研究毗邻图景的可比性问题时，就必须考虑它们的增益。

3.4　太阳高度角引起的辐射误差校正及太阳高度角的计算

太阳高度角引起的畸变校正是将太阳光线倾斜照射时获取的图像校正为太阳光线垂直照射时获取的图像，太阳以高度角φ斜射时得到的图像g（x,y）与直射时得到的图像f（x,y）有下式中的关系，且各波段图像可采用相同的太阳高度角φ进行校正。

张玉君地质勘查新方法研究论文集

太阳方位角随成像季节、地理纬度的变化而变化，通常太阳方位角引起的图像辐射值变化只对图像细部特征产生影响，故可不进行太阳方位角的校正。

由于太阳高度角的影响，在图像上产生地形地物阴影；太阳高度角引起的畸变校正并不能消除地形地物的阴影，我们的对策是将全阴影区及大部半阴影区在预处理时作为干扰区加以剔除，所以在太阳高度角校正时不涉及这些地区。

一般情况下太阳高度角可从数据头文件中得到，在缺失太阳高度角数据时，可用ENVI图像处理软件获得。

3.5　地物视反射率的计算

文献（朱述龙等，2000）给出 p＝πLd²/（E₀×Sinφ）

式中，p为地物反射率；d为日地天文单位距离，可从日地距离表中查到；φ为太阳高度角，可从头文件中得到；E_。为太阳辐照度，文献（冯钟葵，2002）给出了一组参考值，其单位为watts/（meter²×μm）；L为地物在大气顶部的辐射亮度（L=goin×DN+bias），DN为像元值，goin（增益）和bios（偏置）可以从头文件中得到（冯钟葵，2001）。

如果暂时不考虑d和φ，则ρ_k＝π×L/E₀，若以QLMAX和QLMIN分别表示各波段可能的最大和最小像元值（冯钟葵，2002），QL_MAX=255,QL_MIN=1；ρ_kMAX=L_MAXхπ/E₀，ρ_kMIN=L_MIN×π/E₀;L_MAX和L_MIN分别表示各波段的最大和最小光谱辐射值；QL_MAx、QL_MIN、L_MAX和L_MIN均可在地面站Londsat-7产品的Metodata（MTL）文件中查到，由于无论是高增益或是低增益条件下，偏置是不变的，所以L_MIN不随增益改变而改变；各波段LMAx的仅随增益（高或地）的调节在两个水平间变化，如表4所列为低增益状态，表5所列为高增益状态，可看到L_MIN和L_MAX的变化规律。

表4　低增益状态下视反射率换算系数

表5　高增益状态下视反射率换算系数

G_k（%）＝（ρ_kMAX－ρk_MIN）×100/（QL_MAX-QL_MIN）＝（L_MAX-L_MIN）×π×100/[E₀×（QL_MAx-QL_MN）]

若K为视反射率换算系数，则 K=G_k（%）×d²/（Sinφ×G）

B＝ρ_kMIN/G_k,B取整，得B₀。

若某景ETM^＋数据每一像素的DN值经（DN-B₀）×K处理，则该图像准归一为每一灰阶代表G（%）视反射率，即此时灰阶格值为G（%）视反射率

此外还采取化探异常分级的办法，对所提取的蚀变遥感异常做门限化（Thresholding）处理，门限值取数被标准离差，即以σ做为尺度，用数倍σ值做为阈值，限定异常水平，得到羟基异常和铁染异常各三级二值图像。此技术也有助于改进相邻景蚀变遥感异常的可比性。

无论准归一化或是以σ尺度进行异常分级，在执行中都归结为便捷的填表操作，有利于工程化完成，此技术已应用于一百多景 ETM^＋的处理。彩版附图15展示两个时相不同增益（14231H, 14231L）及毗邻景（14231）重叠区色调差异很大，经准归一化后，不仅使三景图色调更为相近，并可用接近于统一规定的2σ获取二级异常，可比性有了明显改善。

参考文献

[1]张玉君，杨建民.基岩裸露区蚀变遥感信息的提取方法[J].国土资源遥感，1998,（2）:46～53

[2]张玉君，杨建民，陈薇.ETM⁺（TM）蚀变遥感异常提取方法研究与应用—地质依据和波谱前提[J].国t资源遥感，2002（4）:30～36

[3]张玉君，曾朝铭，陈薇.ETM⁺（TM）蚀变遥感异常提取方法研究与应用一方法选择和技术流程[J].国土资源遥感，2003,（2）:44～49

[4]李昌国，张玉君.试用主分量分析方法提取澜沧江兰坪地区铜矿化蚀变遥感信息[J].国土资源遥感，1997,（1）:20～30

[5][日]遥感研究会编.刘勇卫，贺雪鸿译.遥感精解.北京：测绘出版社，1993

[6]朱述龙.张占睦.遥感图像获取与分析.北京：科学出版社，2000

[7]冯钟葵.Landsat 7 ETM^＋数据的增益设置及增益改变[J].（中国遥感卫星地面站）用户简讯，2001（3）, 8～9

[8]冯钟葵，关于地物辐射值的计算问题[J].（中国遥感卫星地面站）用户简讯，2002（2）,1

[9]林友明.Landsat-7增益和偏置的计算[J].（中国遥感卫星地面站）用户简讯，2003（3）,7

原载《长春遥感应用会议“国土资源遥感技术发展文集”》，2006。