基于Sentinel-2A的矿化蚀变异常信息提取应用

王　磊　杨　斌　李　丹　陈　财

(西南科技大学环境与资源学院　四川绵阳　621010)

遥感技术经过半个世纪的高速发展与广泛应用，利用多光谱遥感数据提取无人区矿化蚀变信息异常已经成为一种快速、精确、低成本的找矿手段[1-2]。成矿区的蚀变岩石所呈现的遥感光谱信息与金属矿床类型密不可分，提取的蚀变信息异常是一种重要的找矿线索，而遥感找矿正是应用此原理产生的技术[3]，尤其在高海拔、高裂度、高风险的“三高”地区，遥感找矿已经成为地质勘探中一项不可或缺的技术。目前，遥感找矿主要运用Landsat 8数据和Aster数据，基于这两种数据的提取方法已趋于成熟，例如利用ETM数据基于SVM和SAM方法提取矿石蚀变信息，利用ASTER数据进行蚀变岩填图和热液蚀变信息提取[4-6]。运用的提取方法主要有：基于光谱特征的比值变换法(Band Ratio)、主成分分析法(PCA)、光谱角填图法(SAM)、对应分析法(R-Q型因子分析法)、混合像元分解法(MPD)以及目前运用较少的基于核方法的高斯径向基核主成分分析法(KPCA)，通过对比分析以上方法，决定采用提取精度较高的主成分分析法。

Sentinel-2A卫星是欧空局2015年底发射的一颗多光谱遥感卫星，该卫星为“哥白尼”计划下多光谱成像任务中的首颗卫星，与另一颗2017年初发射的Sentinel-2B卫星组成双子星。Sentinel-2A卫星空间分辨率最高可达10 m，拥有13个波段，单星重返周期10 d，在未来数年内将发射Sentinel-3A等持续卫星[7]。由于Sentinel-2A数据具有空间分辨率高、多光谱、重返周期短和免费等优点，在未来遥感研究领域必将占有一席之地，因此选用Sentinel-2A遥感数据，以康定雅拉地区为例进行矿化蚀变异常信息提取，为Sentinel-2A新型光学遥感数据在遥感找矿中的应用提供参考。

1　Sentinel-2A概况

1.1　Sentinel-2A卫星简介

Sentinel-2A卫星于2015年6月23日成功发射升空，在距离地面高度786 km、倾角98.5°的太阳同步轨道上运行，卫星设计寿命为7 a(燃料可维持12 a)，单星赤道重返周期为10 d，“哨兵”双星每隔5 d对全球84°N-56°S范围内的所有地物覆盖成像一次。Sentinel-2A 卫星自发射升空后经过近半年的运行和调试后正式向全球用户提供免费数据下载，Sentinel-2A数据以大宽幅、多光谱、高分辨率等优点成为多光谱遥感应用领域关注的热点[8]，该卫星将作为欧洲SPOT-5和Landsat 8卫星的延续者，在其服役期间将提供全球陆地和海岸带的全部数据，数据主要用于全球高分辨率和高重返能力的陆地观测、生物物理变化制图、环境监测、海岸监测和内陆水域监测以及风险和灾害制图等。

1.2　光谱特性

Sentinel-2A卫星上搭载的高分辨率多光谱成像仪(MSI)是一种基于滤光片的推扫式成像仪，能有效区分可见光(VIS)、近红外(NIR)和短波红外(SWIR)区间的13个波段，成像幅宽为290 km，每10 d更新一次全球陆地表面成像数据[9-10]。Sentinel-2A光谱分辨率为15～180 nm，空间分辨率为10，20，60 m 3种，13个采集谱段中有6个可见光通道，4个近红外短波通道，3个近红外长波通道。为了直观突出Sentinel-2A遥感数据特点，将服役年限、重返周期、成像宽幅、波段数量和空间分辨率基本参数与Landsat 8进行对比分析，结果如表1，表2。

表1　Sentinel-2A和Landsat 8参数对比Table 1 Comparison of two satellite parameters

表2　Sentinel-2A和Landsat 8波段对比Table 2 Comparison of two satellite bands(to be continued)

续表2　Sentinel-2A和Landsat 8波段对比Table 2 Comparison of two satellite bands(continued)

从表1可以看出Sentinel-2A遥感数据在各个方面都有很大的优势，尤其是在影像分辨率和波谱数量上优势明显，这将有利于定量遥感的发展。

Sentinel-2A遥感数据有13个波段，据目前国内外研究表明Band 1 是气溶胶波段，通常用于气溶胶校正，而Band 2，Band 3，Band 4，Band 5，Band 6，Band 7，Band 8和Band 8a是提取植被的良好波段，主要用于陆地覆盖分类、植被监测以及矿物探测，Band 9 是水汽波段，用于水蒸气校正，Band 10是卷云波段，特殊情况下用于洋流探测，Band 11和Band 12适合提取雪、云、冰信息和物理参数反演。其中，4个空间分辨率为10 m的波段可使Sentinel-2A卫星数据同SPOT-5或Landsat 8任务之间保持良好的连续性，以满足用户特别是在土地覆盖基本分类方面的需求；6个空间分辨率为20 m的波段可以满足用户在土地覆盖高级分类和地球物理参数反演等方面的需求；3个空间分辨率为60 m的波段主要用于大气校正和卷云筛选。在波段分布上，在Sentinel-2A卫星数据比Landsat 8分布在可见光区间多1个波段、在近红外多3个波段、短波红外波段数相同，但在中远红外区间，Sentinel-2A无波段分布，较Landsat 8稍显不足。可以看出Sentinel-2A的优势在于可见光和近红外短波的应用，比如遥感找矿、植被监测以及陆地覆盖高级分类和定量分析等应用，并且Sentinel-2A卫星数据在波谱分辨率和空间分辨率上较Landsat 8分辨率更高，可以预见在定量遥感某些领域将取代Landsat 8而得到更广泛的应用。

1.3　数据级别介绍

根据欧空局提供的Sentinel-2A遥感数据相关说明资料并参考国内外学者的论文[11]，将Sentinel-2A遥感数据分为5个层次级别数据，如表3。

表3　Sentinel-2A数据层次级别Table 3 Sentinel-2A data hierarchy level

Sentinel-2A数据处理过程和Landsat 8数据相似，由前端传感器收集数据、压缩传回地面工作站、工作站进行数据解压、数据重排列，然后进行遥感数据特有的辐射校正、重采样、正射影像校正和大气表观反射率计算等一系列操作组成。

2　数据及预处理分析

2.1　研究区概况及数据选取

选取四川省康定雅拉局部地区为研究区域，地理坐标在北纬30°0′-30°20′、东经101°25′-102°0′之间。长期以来，由于受青藏高原急剧隆起抬升和大渡河等河流的切割作用，导致研究区地形异常复杂，其间分布有高山、坡陡和深谷等复杂地形地貌。研究区平均海拔2 500 m左右，气候属于青藏高原型大陆季风气候区，常年气候干燥、降雨量小、昼夜温差大，植被发育较差，积雪覆盖区较少，并且地处康定偏远山区，自然环境恶劣，人口稀少，交通不便，不易于现场勘察找矿。在这样的高海拔、高风险、高裂度的地区进行人工找矿成本高、工期长，人工勘探不具有现实意义，但是遥感找矿正好适用于这种植被覆盖少、海拔较高的山区。根据甘孜藏族自治州矿产资源分布图显示，目前已经探明的矿产资源有金、银、铅等贵重金属矿产，其中岩金矿和砂金矿分布居多，铜矿其次，研究区矿产资源丰富，是用于基于Sentinel-2A数据的遥感找矿的良好地区，并且研究区内Sentinel-2A遥感数据和Landsat 8遥感数据影像清晰，辨识度高，因此选取该地区作为矿化蚀变异常信息提取的试验区。

研究区实验数据采用Sentinel-2A遥感数据和Landsat 8遥感数据，其数据级别分别为L1C和L1T，影像拍摄时间为秋冬季节，植被积雪覆盖少，对提取蚀变信息影响小，两景影像含云量少，不同地物之间的层次感强，色调对比度好，纹理细节清晰。图1为研究区影像，表4列出了影像基本信息。

图1　研究区影像Fig.1　Study area image

数据类型名称波段数数据级别采集时间含云量Sentinel⁃2AS2A＿OPER＿PRD＿MSIL1C＿PDMC＿20161122T152355＿R104＿V20161121T035042＿20161121T03504213L1C2016-11-215．1%Landsat8LC81310392016277LGN0011L1T2016-10-038．3%

2.2　数据预处理

为了提高提取蚀变信息精度，消除大气对地物波谱的影响和便于后期进行矿化蚀变信息处理，在进行矿化蚀变信息提取前，需要对Sentinel-2A影像和Landsat 8影像做图像预处理。Sentinel-2A遥感数据预处理分为大气下垫面校正、图像重采样及L2A级别的数据导出、波段融合及研究区裁剪等三大步骤。传统大气校正忽略了大气下垫面对地物波谱的影响，为了得到更精确的地物波谱信息，使用欧空局提供的SNAP(Sentinel Application Platform，版本5.0)软件中Sen2Cor模块进行大气下垫面校正[12]，其具体操作过程为：首先在SNAP中打开sentinel-2A的原始数据(L1C数据)，其次调用Sen2Cor处理模块(该模块需要用户手动安装)，设置参数Resolution为10 m(可以选择10,20,60 m)，然后进行重采样(Resampling)操作，最后导出L2A级别数据至ENVI中进行矿化蚀变信息提取。Landsat 8遥感数据预处理同样需要对影像进行大气校正、图像重采样和研究区裁剪等步骤，殊途同归，最后也需要使用ENVI进行矿化蚀变信息提取，其具体流程如图2所示。

图2　数据处理流程图Fig.2　Data processing flow chart

Sentinel-2A遥感数据预处理特殊在于需要调用SNAP软件中Sen2Cor模块进行大气下垫面(BOA)校正，并将L1A级别数据转为L2A级别数据[13-14]；经SNAP软件预处理后的Sentinel-2A遥感数据，由13个波段减少至10个波段，缺少Band 1，Band 9和Band 10共3个波段，处理后数据是10个波段单独保存的img文件，在影像空间分辨率上提高到10 m。由于进行大气下垫面校正，一定程度上减低大气散射对地物波谱的影响，间接提高了矿化蚀变异常信息提取精度。为直观了解大气下垫面校正对波谱产生的影响，在SNAP软件中，分别选取Sentinel-2A数据(L1C级别和L2A级别)中植被、裸地、岩石和积雪4种不同类型地物，用这4种典型地物光谱特征图来反映经Sen2Cor模块处理前后的地物光谱变化情况，如图3所示。

图3　地物波谱图(L1C和L2A)Fig.3　Geomorphic spectra (L1C & L2A)

从图3可知，在经过大气下垫面校正后，积雪波谱曲线变化最大，岩石波谱曲线和植被波谱曲线除在0.75 μm到0.90 μm有波动，其他地方大致相同，裸地波谱曲线从可见光到2.5 μm附近都存在着细微改变，从前后变化可以看出L2A数据的光谱特性与真实地物光谱特性更为接近。

3　矿化蚀变异常信息提取

矿化蚀变异常信息提取是以矿床围岩发生蚀变导致其光谱异常为依据而进行的提取方法，围岩矿化蚀变实际上是由于热液经过一定的区域(如断层等)后形成一个矿物组成、化学组成等不同于围岩的一个带，而在形成蚀变的过程中，羟基蚀变非常常见，因此通过提取羟基蚀变异常信息来快速圈定蚀变带，从而确定找矿靶区。根据USGS标准波普库中各矿物波谱特征曲线发现，羟基蚀变围岩中的主要代表矿物有蒙脱石(montmorillonite)、埃洛石(halloysite)、伊利石(illite)和绢云母(muscovite)等4种矿石，这4种矿物波谱特征曲线如图4所示。据矿产资源分布图反映，研究区内分布有数量不等硫矿、铜矿、水晶矿、石灰石矿、铅锌矿和金矿，在这些类型不同的矿床中，其围岩蚀变形成蚀变带时常常伴随着不同程度的羟基蚀变，故可以通过羟基蚀变信息来反演可能会发生羟基蚀变的矿床位置，从而快速精确地圈定所找矿床位置。

图4　羟基矿物波谱图Fig.4　Hydroxyl mineral spectroscopy

羟基蚀变异常信息提取所采用的方法是主成分分析法(PCA)，该方法主要采用“Crosta(克罗斯塔)”原理，是一种根据地物波谱特征和特征向量矩阵来提取目标地物信息的方法。结合矿化蚀变异常信息提取过程中主成分分析法相关理论[15]，发现通常采用具有代表蚀变基团波谱“高低起伏”的4个波段组合进行主成分分析，比如利用Landsat 7数据时，使用TM1，TM4，TM5，TM7这4个波段组合进行主成分分析提取羟基蚀变异常信息。经主成分分析后产生的4个分量，第一个分量囊括绝大部分信息，后面分量信息依次递减，噪声依次增多，通常认为PC1主要是植被等地物信息，PC2或者PC3是矿化蚀变异常信息或云量、阴影信息等，PC4主要是噪声，含有少量特殊蚀变信息。

3.1　Sentinel-2A蚀变信息提取

在基于Sentinel-2A影像提取羟基蚀变信息过程中，常见羟基蚀变矿物在Sentinel-2A的11波段即波长1.6 μm附近有强反射峰，在12波段即波长2.2 μm附近有强吸收谷，因此选择Sentinel-2A数据Band 2,Band 4,Band 11和Band 12这4个波段组合进行主成分分析提取羟基蚀变信息。其具体实现技术途径分为3个大步骤：首先，将下载的L1C“哨兵2A”载入SNAP软件，调用Sen2Cor模块将其处理为L2A数据(也可以直接解压下载数据压缩包，然后在DOS窗口输入代码调用Sen2Cor模块，代码为：L2A_Process xxx --resolution=10；xxx为解压文件路径)，接着在SNAP软件中将L2A数据在OUTPUT菜单下将其输出为 .ENVI格式的数据；然后将得到的数据在ENVI软件中打开，并将需要的若干波段重组得到 .dat数据；最后点击ENVI中的PCA Rotation菜单下的Forward PCA Rotation New Statistics and Rotate子菜单，并设置输入数据、统计文件和数据输出路径、参数默认，对 .dat数据进行处理分析，找到诊断波段并设置合理阈值，接着密度分割、信息提取等操作。羟基蚀变信息主分量选择标准是Band 11和Band 12符号相反，且贡献值应该最大；主成分特征向量矩阵结果如表5所示，根据主分量选择标准可知，选用PC3作为表征羟基蚀变信息增强图像，并对其进行密度分割，提取结果如图5所示。

表5　主成分分析特征向量(Sentin1-2A数据)Table 5 Feature vector of principal component analysis (sentinel-2A data)

图5　Sentinel-2A蚀变信息Fig.5　Sentinel-2A alteration information

图5中红色区域代表Sentinel-2A影像提取的羟基蚀变异常信息，经初步目视解译，所提取的蚀变信息大部分位于裸露岩石区，整体成条带状分布，提取的羟基矿化蚀变信息影像清晰，效果良好，经过ArcGIS面积统计，Sentinel-2A提取的矿化蚀变信息面积约为151.18 km2，约占总面积的6.25%，符合蚀变信息分布规律。

3.2　Landsat 8蚀变信息提取

在基于Landsat 8影像进行PCA法提取羟基蚀变信息异常过程中，根据张玉君等知名学者的研究[16]，采用PCA[2,4,6,7]这4个波段做主成分分析，羟基蚀变信息的主分量选择标准是Band 6和Band 7符号相反，且贡献值应该最大；由表6主成分特征向量矩阵结果及羟基蚀变异常表征分量具有的特性可知，选用PC3作为表征羟基蚀变信息增强图像，并对其进行密度分割，提取结果如图6。

表6　主成分分析特征向量(Landsat 8数据)Table 6 Feature vector of principal component analysis(Landsat 8 data)

图6　Landsat 8蚀变信息Fig.6　Landsat 8 alteration information

图6中红色区域代表Landsat 8影像提取的羟基蚀变异常信息，经初步目视解译，所提取的蚀变信息也大部分位于裸露岩石区，整体成条带状分布，并且和图5中Sentinel-2A影像提取结果分布相似，提取羟基矿化蚀变信息影像清晰，效果良好，经过ArcGIS面积统计，Landsat 8提取的矿化蚀变信息面积为109.95 km2，约占总面积的4.54%，符合分布规律。

3.3　蚀变信息精度评价

在研究区范围内，对比Sentinel-2A影像和Landsat 8影像提取羟基蚀变信息，并结合该地区地质图(图略)进行验证。

从Sentinel-2A数据蚀变异常信息提取结果(图5)和Landsat 8数据蚀变异常信息提取结果(图6)叠加分析可以看出(图7)：两种羟基蚀变异常信息集中分布主要有6个不同区域，从研究区区域构造分析提取出的矿化蚀变信息的走向大体上与主要的线性构造走向一致，呈东南向西北沿断裂带分布，主要分布地点在折多山单元和龙布沟单元，主要岩性为粗粒似斑状黑云母二长花岗岩、中粒黑云母二长花岗岩，岩脉主要有二长花岗岩脉、基性岩脉和辉绿岩脉等，且基本上分布在石炭系、二叠系、泥盆系地层中。通过统计，Sentinel-2A提取结果比Landsat 8提取面积多41.23 km2，较Landsat 8提取精度提高了37.49%。通过地质图分析和查阅该地区的矿产资源分布图证明，Sentinel-2A数据像其他传统数据一样能够精确进行遥感找矿应用。

图7　两种异常信息分析Fig.7　Analysis of Two Kinds of Abnormal Information

4　结束语

通过对四川省康定雅拉地区Sentinel-2A遥感数据进行图像预处理、波谱特征分析和主成分分析的矿化蚀变异常信息提取分析发现，Sentinel-2A遥感数据比传统的Landsat 8遥感数据光谱分辨率和空间分辨率更高，利用Sentinel-2A遥感数据提取出的矿化蚀变信息和Landsat 8提取结果吻合，提取精度高，结合该地区地质图验证，地层分布、区域构造走向基本吻合。

[1]代晶晶. 埃塞俄比亚西部岩浆熔离型铁矿遥感找矿模型[J]. 遥感技术与应用,2012,27(3):380-386.

[2]徐云峰,唐菊兴,郑文宝,等. 扎西康锌多金属矿床外围遥感找矿预测[J]. 金属矿山,2014,(2):91-95.

[3]丁暄,童庆禧,郑兰芬,等. 光谱遥感找矿理论、技术与应用[J]. 地球化学,1992,(1):1-8.

[4]王守志,邢立新,仲波,等. 基于Landsat 8 OLI和GF-1 PMS数据融合的铁染蚀变信息提取[J]. 遥感技术与应用,2016,31(5):950-957.

[5]傅文杰,洪金益,朱谷昌. 基于SVM遥感矿化蚀变信息提取研究[J]. 国土资源遥感,2006,(2):16-19，82.

[6]杨斌,李茂娇,王世举,等. ASTER数据在塔什库尔干地区矿化蚀变信息的提取[J]. 遥感信息,2015,30(4):109-114.

[7]龚燃. 哨兵-2A光学成像卫星发射升空[J]. 国际太空,2015,(8):36-40.

[8]DRUSCH M, DEL BELLO U, CARLIER S, et al. Sentinel-2: ESA’s optical high-resolution mission for GMES operational services. Remote Sens Environ,2012, 120:25-36.

[9]岳桢干. 欧洲Sentinel-2A卫星即将大显身手——“哥白尼”对地观测计划简介(上)[J]. 红外,2015,36(8):34-48.

[10] 时丽娜,藉雪峰,王星. 利用Sentinel-2A数据的西藏阿里冰崩范围提取[J]. 内蒙古煤炭经济,2017,(2):157-160.

[11] 郑阳,吴炳方,张淼. Sentinel-2数据的冬小麦地上干生物量估算及评价[J]. 遥感学报,2017,21(2):318-328.

[12] MUELLER-WILM U. Sentinel-2 MSI—Level-2A prototype processor installation and user manual[Z].2016.

[13] HAGOLLE O, HUC M, VILLA P D, et al. A multi-temporal and multi-spectral method to estimate aerosol optical thickness over land, for the atmospheric correction of FormoSat-2, Landsat, VENuS and sentinel-2 images[J]. Remote Sens, 2015, 7(3):2668.

[14] ZHU Z, WANG S, WOODCOCK C E. Improvement and expansion of the Fmask algorithm: cloud, cloud shadow, and snow detection for Landsats 4-7, 8, and Sentinel-2 images[J]. Remote Sens Environ, 2015,159:269-77.

[15] 佟家兴,杜海鹰,朱福琴,等. 融合主成分分析与光谱角匹配的金铜矿遥感探测方法[J]. 金属矿山,2016,(11):119-123.

[16] 杨金中,方洪宾,张玉君,等. 中国西部重要成矿带遥感找矿异常提取的方法研究[J]. 国土资源遥感,2003,(3):50-53.