航空高光谱CASI-SASI数据蚀变矿物信息提取与应用——以甘肃省敦煌市小金窝子地区为例

孙 雨，刘家军，赵英俊，翟德高，柳振江，张方方，秦 凯，田 丰，康向阳，裴承凯，刘鹏飞

(1.中国地质大学(北京)地球科学与资源学院,北京 100083；2.核工业北京地质研究院，北京 100029；3.遥感信息与图像分析技术国家级重点实验室，北京 100029；4.中国地质大学(北京)地质过程与矿产资源国家重点实验室,北京 100083；5.中核大地勘察设计有限公司,北京 100013)

0 前言

在成矿热液影响下，围绕矿体往往发育远大于矿体面积的围岩蚀变。围岩蚀变是一种重要的地质找矿判据(朱永峰和安芳，2010；Yousefi et. al，2018；刘小雨等，2020)，应用遥感技术进行围岩蚀变研究已得到重视(Sabins，1999；Goetz，2009；王润生等，2011)。航空高光谱遥感数据因拥有上百个狭窄且连续波段，可捕捉围岩蚀变矿物的诊断性吸收谱带(张紫程等，2011)，从而实现了从地物分类到地物鉴别的飞跃。相对于传统多光谱和卫星高光谱遥感数据，航空高光谱遥感数据具有光谱分辨率优于20 nm 且空间分辨率普遍优于5 m的“双高”特点，可以提供米级尺度无缝分布的蚀变矿物类型和空间定位信息，进而有效服务于多种矿产资源的找矿勘查工作(任广利等，2013；汪重午等，2014；叶发旺等，2021)。

航空高光谱数据已在地质领域金属矿床的蚀变矿物信息提取和应用方面取得了较好的实际效果，常用的航空高光谱数据有AVIRIS(甘甫平等，2003；Jain et al.，2019；Rani et al.，2020；Samani et al.，2021)、HyMap(Graham etal.，2017；Kopackova and Koucka.，2017；Salleset al.，2017；孙永彬等，2018；董新丰等，2020)、CASI-SASI(Zhao et al.，2015；刘德长等，2017；叶发旺等，2021)等。前人常用的蚀变矿物信息提取方法有光谱角(SAM)、混合调制匹配滤波(MTMF)等，更多侧重于光谱形态相似度度量，对光谱特征参量在矿物提取中作用重视程度不够，在蚀变矿物信息应用方面尚需加强对标志性蚀变矿物种类的精细识别。

因此，本文在对航空高光谱CASI-SASI数据开展预处理后，开展了高光谱蚀变矿物信息提取方法研究，主要包括数据降维、光谱增值、整体谱形-特征参量协同的矿物填图等步骤，在甘肃省敦煌市小金窝子地区开展蚀变矿物信息提取，分析总结该区蚀变矿物的分布特征和金矿化指示意义，并对蚀变矿物提取结果进行查证，为深化航空高光谱数据地质应用提供了思路借鉴。

1 地质背景

研究区位于甘肃省敦煌市正北约110 km，大地构造分区位置属于塔里木板块北山地块的北山晚古生代多旋回陆内裂谷带(杨建国等，2012)。研究区经历了长期而复杂的地质构造演化历史，成矿地质条件优越，是我国西北地区重要的Au-Cu-Fe-Pb-Zn-W-Sn-V-U-P成矿带，区内金矿(床)点密集，是西北五省重要的金矿集区之一(张发荣等，2014)。

研究区出露地层为中-下泥盆统三个井组(D1-2sg)、第四系(Q)，三个井组岩性主要为灰绿色变质砂岩、角岩化变质砂岩，第四系岩性为冲洪积砾石、砂土。研究区侵入岩体分布广泛，中酸性岩体有泥盆纪明舒井石英闪长岩(δοDM)和三叠纪金滩子二长花岗岩(ηγTJ)，基性岩体为南部的志留纪辉铜山基性杂岩(vSHT)。断裂构造为脆性断层，走向以北东向为主，其次为北东东向及北北西向。研究区有高庙、梭子山、小金窝子、石庙等6个金矿点，均分布在中酸性花岗岩体内部及外围，受近东西向深大断裂及其次级断裂、裂隙控制(图1)。

图1 甘肃敦煌小金窝子地区地质简图Fig.1 Regional geological sketch map of the Xiaojinwozi area in Dunhuang City of Gansu Province

2 高光谱数据源及预处理

本文使用的高光谱数据采集日期为2011年7月6日，采集系统为搭载在运-12飞机上的CASI-SASI航空高光谱测量系统(表1)。飞行前应用室内定标系统获取了传感器的实验室光谱定标系数，飞行时同时获取了CASI-SASI高光谱数据、激光雷达高程数据和地面差分基站数据，飞机过顶时同步获取了地面黑白定标布和定标地物的ASD实测光谱数据。

表1 CASI-SASI航空高光谱测量系统主要参数

CASI-SASI航空高光谱数据预处理包括辐射定标、几何校正和大气校正。辐射定标过程应用实验室光谱定标系数对原始数据进行处理，得到辐射亮度数据。几何校正过程应用差分基站数据和激光雷达高程数据进行像元空间位置处理，去除高程对光谱的影响，得到具有投影坐标信息的高光谱数据。大气校正过程先进行水汽影响和低信噪比波段的去除，再基于地面同步实测光谱应用经验线性法进行处理，得到了高光谱反射率数据。最后，按照研究区范围进行裁剪后得到了用于矿物填图的高光谱反射率数据。

应用高光谱CASI反射率数据，选择中心波长633 nm、557 nm、461 nm的波段分别作为红、绿、蓝波段，制作了小金窝子地区的高光谱真彩色影像图(图2)。高光谱真彩色影像图的色调对比度强且纹理清晰，断裂构造影像特征显著。研究区基岩裸露，有利于开展遥感地质工作。西北部为低山区，发育第四系冲洪积大滩；东部和南部为中低山区，中-下泥盆统三个井组地层呈暗灰色，金滩子岩体呈灰红色纹理细腻，明舒井石英闪长岩体呈深灰色纹理粗糙，不同地物的影像颜色和纹理差异明显。金矿点均位于线状影像特征清晰的断裂构造附近，金矿化明显受断裂构造控制。

图2 小金窝子地区CASI高光谱真彩色影像图Fig.2 True-color image of the Xiaojinwozi area using hyperspectral CASI data

3 高光谱蚀变矿物信息提取

高光谱蚀变矿物信息提取技术流程由数据降维、光谱增值、光谱沙漏处理、端元光谱优选、包络线去除、整体谱形-特征参量协同的矿物填图、阈值分割等步骤构成(图3)。

图3 CASI-SASI数据蚀变矿物信息提取技术流程图Fig.3 Technical workflow showing the extraction of alteration minerals based on CASI-SASI data

3.1 数据降维增值和光谱沙漏

蚀变矿物反射光谱形成的物理机理是阳离子的电子跃迁和阴离子基团的分子振动(Goetz et al.，1985；燕守勋等，2003)，数据降维的目的是结合矿物光谱地质学知识对高光谱数据进行波段筛选重组(郑向涛，2017)，挑选出蚀变矿物光谱特征明显的波段，即CASI数据中心波长为500～996 nm的波段和SASI数据中心波长为2105～2435 nm的波段，分别组成新的CASI和SASI数据。数据降维后，CASI数据的波段数由36个减少到27个，SASI数据波段数则由83个减少到21个。

考虑到SASI数据光谱采样间隔为15 nm，为了获取更为细致的光谱吸收特征，对SASI数据进行光谱增值处理，具体做法是使用三次B样条函数对高光谱数据进行波段间光谱插值增值(刘翔舸等，2010；殷瑞娟等，2013；徐宏根等，2019，孙雨等，2022a)，将光谱采样间隔由15 nm提升到5 nm，增值后的SASI高光谱数据拥有61个波段。

端元光谱的确定有选用实测光谱和选用图像光谱两类方式，本文采用光谱沙漏处理并优选后的图像光谱作为蚀变矿物端元光谱，原因是其为高光谱数据获取时的实际光谱反映，避免了实测光谱矿物不典型及光谱混合干扰等不利因素。采用ENVI V5.3软件内置的光谱沙漏工具对降维后的CASI数据和降维增值后的SASI数据进行光谱沙漏处理(李光辉等，2013)，具体包括最小噪声分离变换(Minimum Noise Fraction,MNF)、纯净像元指数计算(Pixel Purity Index,PPI)、n维可视化分析(n-Dimensional Visualizer)，得到待选图像端元。

3.2 端元光谱优选和包络线去除

绘制待选图像端元的光谱曲线，在综合考虑整体形态、特征吸收峰波长位置、吸收峰数量、吸收峰深度关系等因素基础上，基于矿物光谱专家知识逐一分析各条光谱曲线，优选出地质意义明确的矿物图像端元光谱。从CASI数据得到了褐铁矿图像端元光谱(图4a)，SASI数据得到了绿泥石、绿帘石、短波云母、中短波云母、中长波云母、长波云母图像端元光谱(图5a、5c)。

包络线去除将地物光谱转换到相同的基准背景下，可以有效突出矿物光谱曲线的吸收和反射特征，减弱不同光照条件下反射率强度差异并有效抑制噪声，为矿物精细识别奠定了有利基础(Clark and Roush，1984；Hecker et al.，2019)。本文对由SASI数据得到的绿泥石等6种矿物图像端元光谱进行了包络线去除，得到了包络线去除后的矿物图像端元光谱(图5b、5d)，突出了矿物的光谱吸收特征信息。

褐铁矿图像端元光谱与美国USGS光谱库的褐铁矿光谱较为相似(图4)。整体谱形表现为在500～750 nm先急剧再缓慢上升，750～1000 nm先缓慢下降再持续上升。特征参量表现为存在1个860 nm附近的宽缓吸收峰和1个750 nm附近的反射峰。

图4 CASI数据褐铁矿图像端元光谱(a)与USGS光谱库中褐铁矿光谱(b)Fig.4 Comparison of limonite spectra of the image endmember from CASI data(a)and from the USGS spectral library(b)

绿泥石和绿帘石图像端元光谱的整体谱形和光谱特征参量均较相似(图5)，主要差异在于包络线去除前主吸收峰两侧反射峰的连线形态(表2)。

图5 SASI数据包络线去除前后的蚀变矿物图像端元光谱Fig.5 Image endmember spectra of alteration minerals from SASI data before and after continuum removal

表2 绿泥石和绿帘石的图像端元光谱特征一览表

绢云母图像端元光谱包括4个亚类，分别为主吸收峰位于2200 nm处的短波云母、2205 nm处的中短波云母、2210 nm处的中长波云母、2215 nm处的长波云母(图5c)，在包络线去除后特征更为明显(图5d)。不同绢云母亚类图像端元光谱的整体形态均表现为在2100～2210 nm附近呈总体下降趋势，2210～2345 nm先上升再下降，2345 nm以后多为持续上升。光谱特征参量均表现为存在形态均基本对称的双吸收峰，主吸收峰位于2200～2215 nm附近，次吸收峰位于2345 nm附近。

3.3 整体谱形-特征参量协同的矿物填图

整体谱形是光谱相似性度量中最常用的参数，广泛应用于光谱角填图方法(SAM，Kruse et al.，1993)等矿物填图方法中，可以较好地区分Al-OH矿物和Mg-OH矿物等矿物大类，但难以准确精细区分矿物。光谱特征参量有吸收峰数量、吸收峰位置、吸收峰深度关系等，可用于精细区分矿物种类，但未充分发挥高光谱波段众多的优势。针对二者的优点和缺陷，本文提出了整体谱形-特征参量协同的矿物填图方法，实现了光谱整体谱形的全局特征与光谱特征参量的局部特征的优势互补。

整体谱形-特征参量协同的矿物填图方法将具有n个波段的光谱视为n维空间向量，采用计算向量夹角余弦弧度值的方式进行整体谱形度量，弧度值在(0,1)区间范围内的像元才参与得分计算，弧度值越小则得分越高。对于满足整体谱形度量条件的像元，根据设定的矿物种类判别条件进行特征参量度量(孙雨等，2022b)，满足全部判别条件的像元才能确定矿物种类，得到全部像元的矿物种类得分灰度图。最后，设定合适的阈值进行低端分割(一般阈值θ<0.05)，得到蚀变矿物分布图。

4 高光谱蚀变矿物信息地质涵义分析

在小金窝子地区，基于航空高光谱CASI-SASI数据提取出了褐铁矿、绿泥石、绿帘石、短波云母、中短波云母、中长波云母、长波云母共7种蚀变矿物信息。应用基于空间坐标的多源地学信息复合方法，将高光谱蚀变矿物、地质岩性、断裂构造叠加到CASI真彩色影像上，编制了小金窝子地区高光谱遥感地质综合图(图6)。

图6 小金窝子地区高光谱遥感地质综合图Fig.6 Comprehensive map showing hyperspectral geological information of the Xiaojinwozi area

高光谱CASI-SASI数据提取的褐铁矿矿物在高庙、梭子山、小金窝子、石庙地段均有分布。褐铁矿在高庙地段整体呈密集星点状，局部为团块状产出于三叠纪金滩子二长花岗岩(ηγTJ)，在梭子山北侧呈密集斑块状产出于地势低矮、剥蚀强烈的浅肉红色金滩子二长花岗岩内，在梭子山南侧和东侧呈斑块状、星点状产出于金滩子二长花岗岩内，在小金窝子地段呈云朵状、小团块状产出于泥盆纪明舒井石英闪长岩(δοDM)。在石庙西-石庙南一线，星点状分布的褐铁矿产出于金滩子二长花岗岩内，束带状分布的褐铁矿发育在花岗岩脉及石英脉上。

绿泥石矿物主要分布在高庙南、梭子山南、石庙南地段。高庙南地段的绿泥石呈北东向定向展布，产出于中-下泥盆统三个井组(D1-2sg)，绿泥石的定向分布与地层走向一致，反映了层间断裂的热液活动。梭子山南地段的绿泥石呈斑块状分布，产出于志留纪辉铜山基性杂岩(vSHT)，由辉石等暗色矿物蚀变而成。绿泥石在石庙正南地段呈斑点状产出于三个井组地层内，在石庙西南地段呈点状产出于辉铜山基性杂岩边部。

绿帘石矿物主要分布在梭子山、小金窝子北、石庙北地段。梭子山地段的绿帘石与北东向段断裂关系密切，呈束状、细带状定向分布于三个井组地层内。小金窝子北地段的绿帘石呈小斑块状产出于三个井组地层及明舒井石英闪长岩内，北侧的金滩子二长花岗岩则基本不发育绿帘石矿物。绿帘石在石庙正北地段呈不规则小块状展布，在石庙北东侧因受北北西向断裂控制而呈北北西向束状分布。

短波云母主要分布在高庙-梭子山一线南侧、石庙东地段。短波云母在高庙-梭子山一线南侧的呈细小块状产出于金滩子二长花岗岩的局部地段，在石庙东地段呈小块状、斑点状产出于金滩子二长花岗岩内。中短波云母广泛分布于研究区内，存在两类分布样式。第一类为定向分布特征不明显的斑块状区域性分布，以梭子山南、石庙北东和石庙东地段为代表。中短波云母在梭子山南地段呈总体东西向的飘带状、椭圆状产出于金滩子二长花岗岩内，在石庙北东地段呈斑块状产出于明舒井石英闪长岩内，在石庙东地段呈块状产出于金滩子二长花岗岩内，反映了岩体成岩热液蚀变活动。第二类为具定向性的团块状、条带状局域性分布，以高庙-梭子山一线、小金窝子-石庙南一线、石庙-石庙南一线为代表。中短波云母在高庙-梭子山一线呈北东东向细条带状产出于金滩子二长花岗岩内，与北东东向断裂及节理关系密切；在小金窝子-石庙南一线呈近东西向的细脉状产出于石英脉、花岗岩脉上；在石庙-石庙南一线呈北北西向的细条带状产出于明舒井石英闪长岩及金滩子二长花岗岩内，与北北西向的断裂关系密切。中长波云母主要分布在高庙-梭子山北侧、小金窝子西地段。中长波云母在高庙-梭子山北侧地段呈带状产出于金滩子二长花岗岩内，在小金窝子西地段呈块状产出于金滩子二长花岗岩及三个井组地层内。长波云母主要分布在高庙北东侧，呈北西向、北东向细条状产出于金滩子二长花岗岩内。

前人应用光谱测量技术在金矿化与绢云母化蚀变岩光谱吸收峰位置关系方面开展了一些研究，发现焦家金矿田绢英岩化蚀变阶段的绢云母矿物结构中的Al-OH吸收峰位置存在2185～2225 nm之间漂移(范潇等，2015)。邵雪维等学者对胶东新城金矿蚀变岩的光谱测量发现金矿化岩心段的绢云母Al-OH吸收峰位置在2205～2208 nm(邵雪维等，2021)。本文总结了高光谱CASI-SASI数据提取的4类绢云母亚类在各个金矿区的分布情况(表3)，认为高光谱数据在金矿点提取出的蚀变矿物主要为褐铁矿、中短波云母及绿帘石、中长波云母等，多数金矿点均出现了褐铁矿+中短波云母的蚀变矿物组合，可以作为甘肃省敦煌市小金窝子地区金矿的找矿指示标志。

表3 小金窝子地区金矿点高光谱蚀变-岩性-构造要素一览表

5 高光谱蚀变矿物信息查证

综合采用野外地质调查和ASD光谱实测的方式，选择小金窝子地区蚀变矿物发育的代表性地点，对提取的蚀变矿物种类进行了查证，共设计了16个查证点，对7类蚀变矿物均部署了查证点(图6)。

野外地质调查发现，褐铁矿化蚀变岩呈褐黄色，发育在二长花岗岩内，全岩普遍发生蚀变(图7a)。绿泥石化蚀变岩呈灰绿色，产于深灰色-灰褐色辉长岩小断裂及裂隙附近，未蚀变原岩表面覆盖有暗色沙漠漆(图7b)。绿帘石化蚀变岩呈黄绿色，发育在三个井组变质地层内(图7c)。绢云母化蚀变岩呈灰黄色(图7d)，由于绢云母矿物颗粒细小，肉眼难以直接分辨，主要是根据ASD光谱仪获得的实测光谱吸收特征判别所属绢云母种类。

图7 小金窝子地区蚀变岩石野外照片Fig.7 Field photos showing altered rocks in the Xiaojinwozi area

应用美国ASD光谱仪，开展了蚀变岩样品的光谱实测，各种蚀变矿物反射光谱如图8所示。褐铁矿光谱表现为752 nm附近反射峰和867 nm附近吸收峰(图8a)；绿泥石和绿帘石光谱均存在清晰的双吸收，前者主吸收峰两侧反射肩连线呈下倾状(图8b)，而后者呈上扬状(图8c)；各种类型的绢云母光谱均存在2198～2216 nm附近的主吸收峰和2350 nm附近的次吸收峰(图8d)。蚀变岩的实测光谱与矿物图像端元光谱及USGS光谱库中相应矿物光谱十分相似，从光谱角度证实了本次提取蚀变矿物的存在。

图8 小金窝子地区蚀变矿物实测光谱Fig.8 Field measured spectra of alteration minerals in the Xiaojinwozi area

6 结论

(1)提出了光谱增值、整体谱形-特征参量协同等数据处理新思路，建立了航空高光谱CASI-SASI数据蚀变矿物信息提取的技术流程，在甘肃省敦煌市小金窝子地区提取出了褐铁矿、绿泥石、绿帘石、短波云母、中短波云母、中长波云母和长波云母共7种蚀变矿物，编制了高光谱遥感地质综合图。

(2)航空高光谱CASI-SASI数据提取的大面积区域性分布的蚀变矿物多为区域热液作用产物，小范围定向分布的蚀变矿物与断裂、岩脉等成矿热液活动关系密切，如定向分布的绿泥石和绿帘石矿物反映了与断裂构造关系密切的热液活动，具定向性的团块状、条带状局域性分布的中短波云母矿物与北东东向、北北西向断裂构造及近东西向酸性脉体关系密切，在本区多个金矿点均表现为褐铁矿、中短波云母等蚀变矿物。

(3)基于航空高光谱CASI-SASI数据提取的褐铁矿+中短波云母的蚀变组合可以作为小金窝子地区中酸性岩体内金矿的找矿指示标志，这一规律在北山南带敦煌地区具有普遍性。全国其他成矿带的金矿或其他矿种是否也存在类似规律，尚需开展更多高光谱遥感地质研究进行分析和对比，以期发现普适规律更好指导基础地质科研和地质矿产勘查工作。