秦凯，陈建平，赵英俊，朱玲，杨国防

(1. 核工业北京地质研究院 遥感信息与图像分析技术国家级重点实验室，北京 100029；2.中国地质大学(北京) 国土资源与高新技术研究中心，北京 100083)

0 引言

赤铁矿、针铁矿和褐铁矿等含氧化铁矿物分布广泛，在可见-近红外(VNIR)波段产生诊断性的吸收光谱特征[1-2]。这些吸收特性来源于电子过程，主要包括两种类型：第一种是两个不同离子间的价间电荷转移(intervalence charge transfer，IVCT)，例如三价铁和氧离子的电荷转移。在光谱上，主要表现为480～550 nm的强吸收。第二种是离子的晶体场或能级的电荷转移。对于蚀变矿物岩石光谱，三价铁离子晶体场转移(crystal field transitions，CFT)会在630～715 nm和850～1 000 nm之间产生吸收特征。亨特和阿什莉观察到，铁氧化物的含量与在850～1 000 nm之间的晶体场吸收光谱的深度特征有关[1]，这也被Cudahy和Ramanaidou的实验工作进行了验证，得出结论：Fe2O3的含量与850～1 000 nm的吸收深度相关[2-3]。

可是，采用航空高光谱技术获取的岩矿光谱还与传感器采样间隔、噪声和地形等因素有关。这些因素都可能影响吸收深度，降低氧化铁定量评估的精度。在之前的研究中，Rodger等设计了一种简单二次多项式拟合方法，可以对不同光谱分辨率的高光谱数据进行光谱吸收峰位置和深度信息的提取[4]。Murphy等指出地形会影响氧化铁矿物的吸收深度[5-6]。连续统去除(continuum removal)方法可以有效去除地形对光谱吸收特征的影响[7-9]。除此之外，不同岩性的理化特性可能会对氧化铁在光谱响应上具有不同的效果，所以对于氧化铁的光谱反演，需要不同岩性、不同矿物的混合光谱的先验知识。张琴等以美国地质调查局(USGS)、喷气推进实验室(JPL)和约翰霍普金斯大学(JHU)提供的标准数字光谱库为基础数据，基于去连续统后小波包分解高频和低频部分运用多元逐步回归法、偏最小二乘法建模建立了Fe2O3含量的13个高光谱反演模型[10]。闫柏琨等以亚利桑那州国家大学(Arizona State University，ASU)光谱库为数据源，研究了硅酸盐矿物发射率光谱克里斯琴森特征与SiO2含量的定量关系[11]。

前人的研究大多基于实验室光谱或者地面光谱，针对航空高光谱数据获取、处理和分析的特点，基于赵英俊等在祁连山黑沟铁矿等典型矿区的航空高光谱矿物填图和矿产预测研究基础[12]，结合连续统去除、二次拟合、多元回归分析等方法，构建了基于航空高光谱数据的地表氧化铁含量评估方法，通过铁离子晶体场转移光谱吸收特征和氧化铁含量的回归分析，获得了研究区氧化铁矿物含量分布图。

1 方法

通过借鉴前人方法和结合航空高光谱特点，提出了一套地表氧化铁含量的反演方法，该方法的流程：首先基于大气特征谱线的定标参数进行数据预处理；其次利用连续统去除和二次拟合的方法获取铁离子CFT光谱特征；最后根据地面实测氧化铁含量进行回归分析。流程图如图1所示。

图1 航空高光谱数据的地表氧化铁含量评估流程图

1.1 数据预处理

航空高光谱预处理就是对成像光谱仪获取的原始辐射数据进行定标和大气辐射校正。预处理的目的主要是消除大气和光照等因素对地物反射的影响，获得地表反射率数据。航空高光谱数据的反射率计算主要基于大气辐射传输模型来计算，其中大气成分的参数主要通过地面实测和吸收线评估获取，水汽成分的评估可以采用大气预处理微分吸收法逐像元的评估水汽含量，二氧化碳可以采用地面实测数据，气溶胶采用乡村模型。各参数确定之后，就可以计算出航空高光谱图像数据的各像元的反射率光谱。

由于航空数据采集的特殊性，每次成像光谱仪的安装和飞行过程的差异，会造成实验室定标参数的误差，成像光谱仪本身的smile效应等因素也可能会造成反射率数据的错误。可以根据气体吸收特征等方法迭代进行定标校正，提高反射率的精度。

1.2 连续统去除

在计算光谱特征并反演物理量之前，需进行光谱的连续统去除。利用直线连接光谱起始与终止波长处反射率，搜寻整个谱段上反射率大于连接直线且差值最大的波段，用该波段的反射率代替连接直线在该波段的值，生成连接起始波长、终止波长与该波段的折线，上述过程反复迭代循环，直至在所有波段上，反射率均小于或等于生成的折线，最终形成的折线即为连续统。在可见至短波红外波谱范围内，矿物的吸收特征可以看作为矿物吸收和背景吸收2个部分。其中背景吸收(连续统)是叠加在矿物吸收特征之上，会造成吸收峰位置的偏移。因此，采用去连续的方法来减少背景信息和混合矿物对光谱的影响。连续统去除的方法可以减少地形因素的影响，岩石的反射光谱和连续统光谱都包含相同的地形影响因子，在比值算法计算下，地形影响被消除。

1.3 吸收峰位置和深度计算

本文采用二次拟合的方法计算铁离子吸收峰位置和吸收深度，该方法不受高光谱设备的光谱采样间隔、信噪比等因素影响，波长和连续统去除后的反射率可以用公式(1)来表示：

ρcr(λ)=aλ2+bλ+c

(1)

公式中拟合系数[a，b，c]可以将连续统去除后的3个波段数据通过直接带入法计算。

(2)

(3)

(4)

对公式(1)计算一阶导数，计算导数为零值的波长，将直接带入法计算出的拟合系数带入导数公式。

(5)

(6)

λab为吸收峰位置的波长，进一步求解吸收深度：

(7)

1.4 回归分析

根据地面采样分析测试的Fe2O3含量与对应光谱的吸收特征，基于线性回归的方法建立特定吸收波段吸收深度与氧化铁含量的回归方程。应用回归方程对航空成像光谱数据进行氧化铁含量的反演。

2 实验与分析

2.1 研究区概况

黑沟铁矿位于甘肃省嘉峪关市镜铁山东南14 km处，有简易公路可以到达，地理位置：97°57′31″E，39°18′24″N。

黑沟铁矿区出露地层为长城系桦树沟组和蓟县系镜铁山群，由浅变质的杂色千枚岩、石英岩及铁矿层组成，总体上为一套复式含铁建造，呈北西—南东向展布，为西部桦树沟矿区含铁岩系东延部分。矿区内岩浆岩仅见石英闪长玢岩及辉绿岩。石英闪长玢岩见于矿区南部，侵入灰绿色千枚岩及铁矿层中，与岩层、矿层走向基本平行，局部斜切，呈脉状产出。辉绿岩为岩墙产出，见于矿区西部，侵入灰绿色千枚岩中(图2)。

注：此图据甘肃省地质局镜铁山地质资料修编。图2 黑沟矿床地质图

黑沟矿区本身为一向斜构造，位于桦树沟复向斜东南端，轴向北西西，总体上向南东倾没。矿区范围内，黑沟向斜两端则向中间倾伏，向斜两翼次一级小型褶皱发育，次一级褶皱轴向与主向斜轴向基本一致，属于同期构造。

矿区内断裂构造有北东东和北北东两组，断层面多倾向北北西或北西西。破碎带不甚发育且延长不远，断层整体上南东盘相对升起，在垂直陷落的同时，断层北西盘普遍向南西方向位移，对于矿体破坏不大。矿区内岩层及矿层中，节理发育，其中以铁矿层及灰色千枚岩中的节理尤为发育。

2.2 航空高光谱数据获取和预处理

采用空中国王飞机搭载CASI/SASI航空高光谱传感器，在2014年7月完成了祁连山黑沟铁矿研究区的航空高光谱数据。飞行前完成了成像光谱仪的实验室定标，所获取的航空高光谱数据的参数见表1。完成了工作区CASI/SASI数据的辐射校正、正射校正、几何校正、色差处理、影像拼接和大气校正等数据预处理工作。主要包括6个步骤：辐射校正；传感器姿态数据处理；GPS定位数据处理；姿态数据与定位数据时间同步与集成；精细的几何校正；地形校正；大气校正和光谱重建。其中，采用的是5 m分辨率卫星DEM数据完成了正射校正。大气成分的参数主要通过地面实测和吸收线评估获取，其中水汽成分的评估采用大气预处理微分吸收法逐像元的评估水汽含量，二氧化碳采用地面实测数据357.5 ppm(CO2混合比)，气溶胶采用乡村模型。大气参数确定之后计算出研究区CASI/SASI的各像元反射率光谱。

表1 航空高光谱基本技术参数

2.3 地面光谱数据获取与分析测试

采用可见光-短波红外ASD FieldSpec Pro FR光谱仪采集了研究区地面光谱数据。典型岩石的光谱曲线见图3。

在黑沟铁矿开展了多条地质剖面的调查和光谱测量，获取了地面典型地质单元的光谱数据，用于矿物高光谱矿物信息提取和分析。以黑沟矿区的“3790开采平台”高光谱地面测量剖面为例，主要测量了千枚岩绢云母粉砂质板岩(QLYZ-7A QLYZ-12A)、赤铁矿化砂岩(QLYZ-7B QLYZ-8B)、绿泥粉砂质板岩(QLYZ-8A)、铁白云岩(QLYZ-9B、QLYZ-10B、QLYZ-11A)和铁矿铁矿石(QLYZ-11B)等岩石的光谱。有以下特征：

(1)千枚岩绢云母粉砂质板岩的光谱曲线显示具有3个特征吸收位置，分别位于2 210 nm、2 345 nm和2 440 nm，其主要吸收峰位于2 210 nm，与绢云母矿物的铝羟基基团振动有关。

图3 典型蚀变岩石光谱吸收特征

(2)赤铁矿化砂岩由于其赤铁矿的三价铁离子的电子作用，光谱显示特征吸收峰位于880 nm，同时在其左肩存在770 nm的强反射峰。

(3)绿泥粉砂质板岩的光谱显示了铁羟基2 255 nm和2 360 nm的吸收特征组合。

(4)黑沟矿区碳酸盐主要为铁白云石，碳酸根特征吸收位置分别位于2 330 nm和二价铁的吸收谷从1 000～1 300 nm。

(5)铁矿石主要矿物为镜铁矿及少量磁铁矿，无明显吸收特征，该样品表面存在少量褐铁矿化现象。

从以上光谱分析可以得出，矿区内的岩性光谱差异较大，赤铁矿、绢云母、绿泥石、白云石等矿物光谱特征明显，以上光谱特征位置、深度和组合关系可以作为航空高光谱信息提取的依据和参考。对测量点进行了地面采样，并开展了地面样品的薄片鉴定(图4)和分析测试，部分测试结果见表2。

图4 样品光谱薄片鉴定结果

表2 典型矿床剖面样品主要化学成分含量%

2.4 吸收特征计算与回归分析

选择850～1 000 nm波谱对航空高光谱数据进行连续统去除，对航空和地面数据采用二次拟合方法，计算航空光谱数据的吸收深度位置和相对应的吸收深度，利用线性回归方法，对采样点的光谱吸收深度和Fe2O3含量进行反演。图5显示了不同采样区的反演结果。

铁矿层的主要岩石含矿砂岩、含矿硅质岩，图5(a)为12件样品的光谱与氧化铁含量的相关关系，铁矿石光谱铁离子晶体场吸收深度特征与Fe2O3含量具有明显的线性关系，相关系数R2=0.871 4；在铁矿层周边的褐铁矿化、赤铁矿化白云岩、赤铁矿化千枚岩、赤铁矿化砂岩采样18件，图5(b)显示的光谱铁离子晶体场吸收深度特征与Fe2O3含量具有明显的线性关系，相关系数R2=0.833 4。矿床区域千枚岩氧化铁含量较低，线性相关度低，相关系数R2=0.181 8(图5(c))。

图5 不同采样区吸收深度 (吸收峰位置900 nm)与氧化铁含量相关关系

3 实验结果

对黑沟铁矿的航空高光谱数据进行了数据预处理、连续统去除、二次拟合等处理，根据地面和航空高光谱数据二次拟合计算的吸收峰位置和深度结果，获取了850～1 000 nm的吸收峰位置和深度的铁离子晶体场转移光谱特征。通过对不同采样区域的回归结果分析，引入其他波段的光谱变量，开展了多元回归迭代分析，最终，建立了研究区的航空高光谱地表氧化铁含量的评估公式：

Y=59.42+122.94X1-237.49X2

(8)

式中：Y为氧化铁含量预测值;X1为吸收峰位置900 nm的光谱吸收深度;X2为2 100～2 280 nm的平均光谱。回归系数R2=0.965 2，标准误差为4.44wt%。表3为部分采样点回归分析的预测值和残差。

根据回归分析的预测结果，计算了研究区地表氧化铁含量丰度(图6)，可以更加清晰地展示出黑沟铁矿该处矿体的分布规律，氧化铁含量的区域特征与矿体分布形态基本一致，丰度较高的矿物信息表现为长轴北西向的封闭环形。

表3 氧化铁含量回归分析预测结果及误差

图6 黑沟铁矿区赤铁矿丰度图

4 结束语

本文面向航空高光谱遥感影像，充分利用氧化铁矿物的吸收光谱特征，提出了一种地表氧化铁含量的评估方法。实验结果表明该方法能有效地反演出地表氧化铁含量。本文方法在祁连山黑沟铁矿进行了应用，对矿区内的铁含量评估较为准确，该方法仍然需要地面的高光谱数据和地质背景知识，如何实现在未知区域上快速分析和评估地表氧化铁含量将是进一步研究的方向。

致谢：特此感谢核工业北京地质研究院刘德长研究员、中国国土资源航空物探遥感中心王润生教授的对工作的指导和支持。