高光谱遥感数据的处理及其在铀资源勘查的应用

—以广西苗儿山地区为例

刘德长 谢红接 赵英俊 张进业

0引言

高光谱遥感数据是80年代以来发展起来的新一代(第三代)遥感数据。它使遥感技术步入可以同时获取地球表面物质成分信息和空间分布特征信息的新阶段。本文以广西苗儿山地区机载成像细分红外高光谱数据为例，探讨了高光谱遥感图像的预处理和深化处理，野外实测光谱曲线库建立与典型地物的波谱特征，以及高光谱数据在铀资源勘查中的应用等。

1 高光谱遥感数据的预处理和深化处理

1.1 苗儿山地区高光谱遥感数据特点

上海技术物理所航空红外细分光谱扫描仪是以地质遥感为主要应用目的开发的，在2～2.5um光谱区间细分为12个工作波段。广西苗儿山地区高光谱遥感数据，由中南地质局230所和上海技术物理所于1988年11月2日、4日获取。平台为伊尔14型飞机，航速为250km/h，航高（绝对航高）达4500m,沿NE30。方向飞行。所选8个通道的中心波长分别为1.600μm、2.143μm、2.200μm、2.250μm、2.300μm、2.330μm、2.370μm和2.450μm。通过分析，苗儿山地区高光谱遥感数据有如下特点：

（1）由于仪器本身的原因，所有数据沿航向具周期性条带噪声；

（2）与TM和MSS等图像不同，高光谱图像在同一航带的各波段之间，沿航向和扫描方向有几个至十几个象元之位移；

（3）由于航向分辩率（8.68m）和扫描方向分辩率（27m）的不一致，象元呈长方形，其长（航向）宽（扫描向）比为27:8.68=3.11:1;

（4）由于飞行方向、扫描角度、太阳方位角和太阳高度的几何关系的变化和大气的影响，扫描图像存在边缘辐射差异；

（5）由于飞行速度不等和扫描角度的变化引起图像几何差异；

（6）第7航带数据由于记录方式的错乱，引起图像与正确方式所成的图像呈镜向对称关系；

（7）高光谱图像的一个显著特点是能够从图像中提取任意象元的光谱率曲线，但由于成像光谱仪获取的图像象元亮度值并不能代表地物的反射率值，需要进行反射率转换。

1.2 高光谱数据的预处理

在研究上述所有特点的基础上，开发了一系列数据预处理软件，并进行了条带去除、波段间配准、航向压缩、镜向变换、辐射校正、相对反射率转换等预处理，最终得到清晰的图像。

1.2.1去条带方法研究

（1）原始数据分析 原始数据存在严重的条带干扰，垂直于航线方向分布，且从第一通道到第八通道条带越来越增强。经分析，这种条带是由于在成像过程中机械或光学部件不稳定造成的，具有明暗交替的特征，并具有一定的周期性，以9个象素为周期，呈正弦波状，部分以8个象素为周期。

（2）干扰条带的消除方法 可以将干扰条带作为加性噪声处理以消除条带，要求尽量减少信息损失，一般有两种处理方法：

①平滑法，可以采用9点平滑或9×9点，但9×9窗口平滑效果不好，在航线方向上以9点平滑效果好，但也会损失有效信息。

②滤波法，从前面的分析可知干扰条带的周期为9，通过不断自动调整波形振幅，可以有效地消除干扰，研究过程中使用了滤波法。

（3）效果分析图 图1左图为原始图像，右图为去条带后的图像，对比左右两侧的图像可以明显地看出干扰条带被有效地消除了，而原始数据信息的损失比较小。条带消除后图像的列数保持不变，但行数减少，如第5航带行数由5712减为5632。

图1去条带影像效果对比图

（左侧为原图像，右侧为处理后图像）

1.2.2 波段之间配准

一般情况下，TM图像的波段与波段之间都是完全重合的，不存在象元位移，而苗儿山高光谱遥感图像存在着明显的象元位移，一般为几个→十几个象元。从5、6、7、9条航带来看，都是第2波段（2.143μm）和第3波段（2.2μm）之间象元位移量最小，如第5航带的第2，3波段则完全重合，这可能说明该仪器的第2，3通道工作状况最好。下面以第9航带为例来说明波段之间的位移大小及配准方法。以第2波段为基准，其他波段则向左右或上下位移。

以第2波段为基准

第1波段 右移1个 下移 7个

第3波段 左移2个 上移4个

第4波段 左移4个 下移14个

第5波段 左移4个 下移7个

第6波段 左移5个 下移8个

第7波段 右移2个 下移3个

根据以上位移量截取完全重合的图像部分如表1，得到第9航带最高配准后图像大小为列240，行5738。

表1 第9航带各波段截取部分

	列数	行数
第1波段	9+2：256-6	8：6125-11
第2波段	9+3：256-5	15：6125-4
第3波段	9+5：256-3	19：6125-0
第4波段	9+4：256-4	1：6125-18
第5波段	9+7：256-1	8：6125-11
第6波段	9+8：256-17	7：6125-12
第7波段	9+1：256-17	13：6125-7

（说明列数的9是由于每幅图像都有9象元的亮边，所以都被裁掉）

1.2.3 航向压缩和扫描方向拉伸处理

机载成像光谱仪沿飞行方向的分辩率由飞机的飞行速度决定，其速度为250公以/小时（即250000m/3600s），扫描率为每秒8线，2条线之间的距离为25000/3600/8=8.68，即航向分辨率为8.68m。扫描方向分辨率可由瞬时视场（6毫弧度）和飞行高度（4500米）给出，即（6/1000）×4500=27m。由此可见，图像上每个象元不是呈正方形，而是呈长方形。长、宽分别为航向和扫描方向，长：宽=27：8.68≈3.11:1。从视觉效果来看，显示出沿航向拉长或沿扫描方向压缩（扁）的现象。通地沿航向的压缩或沿扫描方向的拉伸处理（数据重采样方法用三次样条法），可得到27×27或8.68×8.68的正方形象元。 1.2.4 边缘辐射校正处理

由于飞行方向、扫描角度、太阳方位角和太阳高度的几何关系的变化及大气的影像，扫描图像存在边缘辐射差异。这种差异可以通过几何光学模型和大气模型予以校正，但由于模型十分复杂，影响因素很多不易操作，因此，通常采用统计方法进行图像的边缘辐射校正。设有图像矩阵Anm

其中n为行数（沿飞行方向），m为列数（沿扫描方向），对图像每列求均值得到向量：

反映了整个图像边缘辐射差异的平均趋势，作者以图像中机下点的辐射量为标准对其余各列像元进行归一化处理，得差值向量：

校正后得图像为：

1.2.5 正切校正处理

由于扫描角度的变化，图像中每个象元对应的地面分辨率随着变化，这种几何差异的校正，通常使用正切校正方法。

h NP1 P1 NP2 P2 O 扫描角度变化如图2：

图2 扫描角度变化示意图

其中0为传感器，NP1，NP2为实际扫描面，P1、P2为地面。传感器面向正下方等角速地收集地面数据， 为瞬时视场角，h为传感器与地面之间的距离。

正切校正算法如下：

设图像象元为等间隔网格上（即正方形象元）的点，则在扫描方向图像网格，宽度h× 。从机下点开始计算，第i点网格与机下点的距离为i×h× 。实际上，地面网格点与机下点的距离为h×tan（i× ）。作者按照网格宽度h× 对图像进行重采样，设重采样后图像网格点i对应图像网格点K，则有：

i·h·＝h·tan（k·）

求出K：

通常K不会落在网格点上，即K不是整数，可以得出i，满足：

≤K≤ ＋1

采用距离倒数插值求出K点的象元值：

校正后的图像象元值DN（i）=DN（K）

显然通过上述方法校正后，图像的列数将增加（图3）。

图3 正切校正前后对比影像图

（左图为原图像；右图为校正后图像）

1.2.6 镜像变换处理

从己经处理过的4条航带数据来看,只有第7航带记录错乱,显示出来的图像与实际图像呈镜向对应关系。可用下式表达：

f（m，n）= f（M-m，n）

m=0,1,…,M-1； n=0,1,…,N-1

其中m为列数，n为行数，M、N为该图像总的列数和总的代数，插值方法用的是三次样条法（图4）。

图4 镜像转换前后对比影像图

（左图为原图像：右图为转换后图像）

1.2.7 相对反射率转换处理

高光谱遥感图像不同于多光谱遥感图像的一个根本特点在于从高光谱图像立方体里能够提取出精确的地物光谱反射率曲线。据此，可以直接识别地物。因此，准确地提取影像中的光谱信息是图像分析的基础，也是成像光谱技术研究的一项重要内容。

成像光谱仪在0.4-2.5μm工作波长范围内主要光谱能量来自地球表面对太阳光的反射。由于太阳光受大气吸收影响，太阳光的地面辐照度随波长变化大。因此，成像光谱仪所获得的图像DN值不能代表地物的反射率值，需要进行反射率转换。

采用了2种方法进行相对反射率转换，下面介绍其原理：

（1）统计法 这是针对没有进行同步地面反射率测量采用的一种方便、快捷的方法，常用的有对数残差转换法和内在平均相对反射率法（IARR）。研究过程编制了对数残差转换法软件。

设有高光谱图像DN（I,j,b）,i=1,2,…M为行数，j=1,2….N为列数，b=1,2,…k为波段数。

在i行j列有光谱量均值：

波段b的灰度均值：

现假设：

①整个图像区域中各类成份足够杂乱；

②对征一波段的整个区域取平均，所得的光谱与地面吸收特征无关，则可用对数残差转换，得相对反射率；

logRy＝logDN（i,j,b）一logSij一logDb

（2）地面定标法 在飞行的同时，用野外光谱仪实地测量典型地物的光谱，与成像光谱图像中对应地点象元的DN值相对照，求出图像DN值到地面反射率的转换公式。理论与实践都表明图像DN值与地面反射率有以下近似的关系：

其中DN为波段b给定象元的数值，p（λ）为波段b的波长λ对应的地面反射率，Ab为乘积项，反映了大气传输及仪器设备的放大比；Bb为偏移项，反映了大气辐射值及仪器的零点偏移。Ab、Bb可以通过地面定标数据及相应位置的图像DN值来确定。

设有N类地面定标点的反射率，则有一组方程：

 i=1，2,…，n

求此方程组的最小二乘方解得：

本次使用的数据是1988年11月初的飞行数据，但当时没有做同步地面测量。作者于1995年11月初进行了补充测量，一方面在于了解工作区的典型地物的光谱特征，另一方面也想做一次似同步野外定标实验。全部测量数据于1995.11.4-11日获得，这期间天气晴朗，采样时间11-16点，大体与1988.11.2、4日的时相相当。以第五航带为例，当时选取了资江水，张家矿区7号带、烟竹河及公路交叉处的路边白沙及附近正常花岗岩作为定标点。

表2 野外定标点及其图像坐标表

表3 各点在对应波段上的反射率值（%）和图像DN值（0-255）表

这样依此3点可以得到一条回归曲线,可将图像上任一点的DN值转化为相对反射率值,并编制了相应的软件。上述3点转换后的相对反射率值（表4）如下：

表4 转换后的相对反射率值表

坐标点	1600 nm	2143 nm	2200 nm	2250 nm	2300 nm	2330 nm	2450 nm
1061 1153	3.568	3.278	4.441	4.092	2.828	0.050	7.351
981 956	69	65	67	66	54	54	77
979 957	31	29	22	30	39	39	32

2. 波谱特征

野外地物波谱测量是遥感应用的基础，而高光谱分辩率数据具有图谱合一的特征，因此波谱测量对航空高光谱数据的处理和解释尤其重要。一般情况下，在飞行时应同步进行野外地物波谱测量，为野外定标及数据处理和解释提供依据。但由于条件限制，如前所述飞行时没有同步进行野外地物波谱测量。为了弥补这方面的不足，研究过程中拟同步性地补充进行了野外地物波谱测量。在EVNI软件上建立了野外实测光谱曲线库，并研究了该区典型地物的波谱特征。

2.1 地物波谱测量及影响岩石反射光谱特征的主要因素

2.1.1 野外地物谱测量

利用美国地球物理环境研究公司（GER）研制的IRIS红外智能波谱仪对苗儿山地区的典型地物进行了波谱测试，共获取了150多条波谱曲线，基本上反映了研究区各典型地物的波谱特征，并为野外定标提供了基本数据。

野外光谱测量是在1995年11月4日至11月11日进行的。如前所述，这个时间段与当年获取数据的时间相当，便于地面所测光谱信息与航测的光谱信息做对应分析。在测录中用GPS在1：5的地形图上做了较准确的定位，并且对每一个测点的地物类型做了细致的描述，使以地面波谱测量为基础的处理比较可靠。为了分析测得的光谱曲线，先分析一下影响岩石反射光谱特征的主要因素。

2.1.2 影响岩石反射光谱特征的主要因素

⑴岩石成分的影响 由于岩石的成因不同，岩石的成分也随之发生变化，岩石呈现出不同颜色。例如，酸性中花岗岩以长英质矿物为主，岩石呈浅色调；基性岩中的辉长岩以铁镁质矿物为主，岩石呈暗色调。前者的光谱反射率高于后者。从基性岩到酸性岩，它们之间的反射率可在70%范围内变化，深色岩石的反射率要低于浅色岩。可见岩石光谱的反射率的高低与其颜色有关，更确切地说应与岩石的成份有关。

⑵岩石结构、构造的影响 岩石在生成过程中由于物理化学条件不同、生成不同程度、不同结晶度的岩石。岩石的结构、构造反映岩石中不同矿物颗粒大小及它们之间的关系。通常，岩石的结晶度和粒度不同造成同类岩石光谱反射率的不均匀性，如粗粒花岗岩的光谱反射率变化要大于细粒花岗岩。

⑶蚀变和风化作用的影响 在内营力和外营力作用下，岩石发生蚀变或风化，生成蚀变和风化矿物，并引起岩石结构、构造和色调等的变化。

2.2 苗儿山地区典型地物波谱库的建立及地物特征分析

2.2.1 典型地物波谱库的建立

从野外获取的150多条光谱曲线按类别可分成11个大类，并以此建成11个库。它们分别是正常花岗岩、花岗岩风化而成的沙土、花岗岩弱蚀变破碎物、花岗岩强水云母化、花岗岩与草混合物、矿石、硅化带、稻田、水体、各类植被、混合植被等。

为了使实测曲线（0.4-2.5um）能与从高光谱图像立方体中提取的曲线（1.6-2.45um）进个对比分析和匹配研究,须将原曲线按高光谱图像的七个波段及相应的带宽（100nm或50nm）进行重采样。 2.2.2 地物波谱特征分析

根据建成的库，对其波谱特征进行了分析。取各库的平均光谱曲线，如S1、S2、S3、S4，分别代表正常花岗岩、花岗岩风化而成的砂土、花岗岩弱蚀变破碎物、花岗岩强水云母化等。

S1代表正常花岗岩 总体反射率较低，在20%左右，曲线较平坦。

S2代表风化后成的沙土 总体反射率在40-500%，在2200处有较明显的吸收峰。

S3弱蚀变花岗岩 总体反射率比S2稍低，在2200nm处有较明显的吸收峰。

S4强蚀变花岗岩 总体反射率明显升高，特别是在低波段区间反射率升高更为明显，在2200nm处有很明显的吸收峰，880nm处有较明显的吸收峰。

S5花岗岩型铀矿石 总体反射率同正常花岗岩S1，在20%左右，但不同的是在2200nm处有明显吸收峰。

S6花岗岩与各种草混合后的曲线 在670-760nm之间有明显的陡坡，反映植被的典型特征，即所谓的“叶绿素红边”。植被所占比例越多，这条斜边越长，反之，越短。

S7硅化带的特征曲线 总体反射率<40%，在670nm、2200nm有吸收峰（图5-6）。

图5 硅化带的光谱特征曲线

图6 重采样后硅化带的特征曲线

S8稻田反射率特征曲线 其中，S81与植被曲线（S10）有非常近似的特点，在670-760nm之间有明显的陡坡，在1500nm处有明显吸收峰，在1650nm和2250nm处有明显反射峰，反映了田中稻草梗的反射率特点；S82反映的主要是田中泥土的反射率特征，由于它们源于花岗岩风化破碎物，故近似于花岗岩特征曲线。

S9水体 资江水，梅溪水的特征曲线，反射率低。

S10各种植被及混合植被的特征 由于叶绿素在670nm红波段具强吸收，在760nm近红外波段具高反射的特点，便得各种植被在670-760nm之间具有明显的陡坡，即所谓的“叶绿素红边”。它是健康植被的一个特征，如果植被或农作物由于受到如重金属、碳氢化合物、干旱、病虫害、缺氮等因素引起胁迫（stress），则会使这些植被或农作物的特征曲线产生微小的变化，即所谓的“红边蓝移”现象。利用“红边蓝移”现象可以预测病虫害等的发生，并及时采取补救措施。本区的某些植物，如乍子树、竹叶、松树、杂草等都具有此种现象（图7），说明上述植物受到胁迫。

图 7 植物受到胁迫后的光谱特征曲线

（上图为乍子树、竹叶、松树；下图为杂草）

3 高光谱数据的地质应用及效果

3.1 高光谱分辨率图像上硅化带的识别及效果

高光谱分辨率数据的一个最重要用途是识别蚀变矿物。尽管由于原始数据和研究区植被覆盖方面的原因，使应用效果受到了一定的限制，但在硅化带的识别方面仍取得了成功。

图8为第5航带第2波段（中心波长2.143μm）图像（已经各种处理和精校正），从该图像上可以清楚地看到许多近东西向的暗色调线性体，这些线性体在对应区域的TM图像上没有或仅略有显示，在原有的大比例尺地质图上也只有少数几条被填绘出来。通过对比图像上直接提取的已知带（硅化断裂带）和未知带的光谱曲线，发现两者具有一致的曲线特征。因此，推测这些未知的线性体可能是硅化断裂带。

作者对分析结果进行了野外检验，证实它们大多数确系硅化断裂带，其中马家南边的那条硅化带规模很大，宽10～20m，主要为白色石英脉带（铀含量为26 x10－6～45义10—6），带内有残留碎裂花岗岩（铀含量40x10－6，具红化、高岭土化和弱硅化），围岩为碎裂花岗岩（铀含量为60 x10－6一70 x10－6）。硅化带内残留碎裂花岗岩的铀有丢失，而在围岩中有富集现象。在豆诈山岩体及其周围的高光谱图中，解译出的线性构造包括了该区地质图上填出的主要硅化断裂带，其中有些硅化断裂带直接控制了铀矿化。

图8 第5航第2波段（2.143μm）图像

3.2 高光谱分辨率图像与其他多源信息的综合分析及效果

将高光谱数据、TM数据与已有地质图等多源信息进行融合，可形成一幅集多种信息于一体的新类型图像（图9）。

图9 综合解译图

在该图像上，地质图上已有的主要断裂构造和硅化带都已被解译出来，同时，还解译出很多线状构造和环状构造，而这些构造在地质图上是没有的。如fl—沙子江-茶坪断裂构造，f2一断裂构造，f3一断裂构造，f4一牛大江硅化带等；环状构造有4个，分别以C1、C2、C3、C4表示。 图9 综合解译图 其中以Cl规模最大。本区是岩浆岩分布区，环状构造的形成一般与岩浆的侵人活动有关。4个环状构造都发育在两期岩体的交界部位，其中C2和C3是豆诈山岩体与香草坪岩体之接触部位，C4是牛塘坳岩体（γ53－1）、香草坪岩体和苗儿山岩体之接触部位，可能为更晚期的岩浆侵人活动所致。C1不仅规模最大，且结构复杂，该环形构造西南部处于豆诈山岩体与香草坪岩体之接触带；东北部受香草坪岩体、苗儿山岩体和牛塘傲岩体三者的影响；东南部基本上是香草坪岩体和苗儿山岩体的分界线；西北部分还同时与小环C3相交。此环形构造还是本区四条大硅化断裂带（F1, F2, F31, F32）的夹持部位。因此，该环形构造区具有很好的构造、铀源（γ51、γ52－2）和成矿的先决条件，双滑江矿床和白毛冲矿床就在此环形构造范围内。前者产于该环形构造的环边上。本次研究预测的一处铀矿化远景区（P1）就在该环的西北部分。由此可以认为，该环形构造的展布范围应该是很有成矿前景的地区。

经过上述各种处理及综合解译，结合前人对本区成矿机理的认识，本次研究预测了三个成矿有利区，通过验证，其中的两个地区放射性测量值偏高，并有铀矿化显示。

参考文献

1．郑兰芬 王晋年，1992，成像光谱遥感技术及成像光谱信息提取的分析与研究，《环境遥感》，Vol. 7，No. 10。

2．王晋年郑兰芬，1996，成像光谱图像光谱吸收鉴别模型与矿物填图研究，《环境遥感》》，Vol．11，No. 10。

3．王向军，1995，成像光谱信息可视化表达系统，中国科学院遥感应用研究所硕士学位论文。 4．王钦敏，1995，卫星遥感技术发展和应用动态，《遥感信息》，Vol．40, No. 4。

5．R. J. P. Lyon, 1996，风化及其它类荒漠漆表面层对高光谱分辨率反射率的影响（一），《环境遥感》，Vol. 11，No. 2。

6．Joseph C.H, Chein-I Chang, 1994, Hyperspectral Image，Classification and Dimenesionality Reduction：An Orthogonal Subspace Projection Approach,《IEEE Transaction on Geoscience and Remote Sensing》，Vo1.32, No.4（779-784）。

7．Geotz, A.，and Herring, M:，1989, The High Resolusion Imaging Spectrometer （HIRIS）for EOS，《IEEE Transaction on Geoscience and Remote Sensing》，Vol.27, No.2。

8．Mazer, A.S.,et al, 1988, Image Processing Software for Imaging Spectrometry Data Analysis.《Remote Sensing of Environment》,Vol.24，No. l（201-210）。

9．William H.F.，et al，1994, Retrieval of Apparent Surface Reflectance from AVRIS data：A Comparison of Empirical line, Radiative Transfer, and Spectral Mixture Methods,《Remote Sensing of Environment》》，Vo1.47（311-321）。

10．ENVI User's Guide, 1996,Research Systems, Inc.,USA

11．A Hyperspectral Remote Sensing Primer, 1996, NASA, USA。