基于多光谱差分植被抑制法的蚀变信息提取

赵恒谦，高尉，杨天艺，刘帅，赵学胜

1.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083；2.煤炭资源与安全开采国家重点实验室，北京 100083

蚀变信息提取是遥感在地质领域最成功也是最重要的应用之一[1]。根据蚀变岩和正常岩石在遥感影像中特征波段的光谱响应差异，可以提取对找矿有重要指示意义的蚀变矿物信息，如铁染蚀变、羟基蚀变、碳酸岩蚀变等[2-6]。但是，已有研究主要在裸露或半裸露岩层区域开展[7]。在植被覆盖区域遥感影像中的蚀变矿物信息非常微弱，植被信息的干扰在很大程度上影响遥感蚀变信息的提取精度，制约了遥感找矿技术的进一步应用推广[8]。

为了解决植被信息对遥感蚀变信息提取的影响，许多学者以ETM+/TM和ASTER为数据源，在植被覆盖地区开展了多种遥感蚀变信息提取方法，如最佳波段+主成分分析、无损线性拉伸+抑制干扰因素+波段比值法、主成分分析+密度分割、MNF变换+矿物标识等方法[9-11]。然而，这些研究方法都是针对中低植被覆盖区，对于高植被覆盖区的蚀变信息提取效果不明显。在高植被覆盖区，遥感图像分析方法提取的岩矿信息多是小目标、弱信息[12]。因此，又有学者在分析矿物、植被、土壤地面光谱的基础上，采用Hyperion高光谱卫星数据对高植被覆盖区进行了植被抑制和遥感蚀变信息提取研究，证实基于差分的植被抑制方法可以有效地抑制植被信息并提取岩矿弱信息[13]。然而，Hyperion高光谱卫星数据幅宽较窄，无法获取大面积区域的蚀变矿物分布情况。

目前主流的遥感卫星传感器是多光谱传感器，多光谱数据空间分辨率和时间分辨率都较高，幅宽面积大，数据质量好且数据源众多[14]，其中Sentinel-2数据是近年来比较新且有代表性的卫星多光谱数据。Sentinel-2由两颗相同的卫星哨兵2号A(Sentinel-2A)与B(Sentinel-2B)组成，具有多光谱、重访周期短(双星组网5 d)、空间分辨率高、数据幅宽范围大及免费获取等优点，已经在裸露岩层区蚀变矿物信息提取研究中得到应用[15]。然而，在高植被覆盖区，如何去除植被干扰获取蚀变矿物信息仍是一个难题。因此，本文采用Sentinel-2A作为卫星多光谱数据源，以黑龙江呼玛地区作为实验区，通过植被与矿物反射率做差分计算的方法来抑制植被信息，利用主成分分析方法提取大面积高植被覆盖区的蚀变矿物信息。

1 实验区和数据处理

本研究选取黑龙江呼玛地区作为研究区。该区域位于E126°—127°/N50°—52°的大兴安岭东麓，地势呈西北高(丘陵区)、东南低(平原区)，平均海拔350 m。研究区植被以林木为主，植被覆盖率达95%以上；土地肥沃，有棕壤、黑土、草甸沼泽土和草甸土等丰富的土壤资源；同时，矿产资源也十分丰富，盛产金、铁、煤、石墨、石英、云母以及高岭土等。

Sentinel-2A卫星于2015年6月23日发射升空，设计寿命为7年，轨道高度786 km[16]。Sentinel-2A卫星上携带的高分辨率多光谱成像仪(MSI)，覆盖可见光、近红外和短波红外13个光谱段(表1)，卫星重访周期为10 d，成像幅宽达290 km。本研究使用的影像获取时间为2017年10月26日，产品等级为已几何精校正的L1C级，下载自欧洲太空局官网(https：//scihub.copernicus.eu/dhus/#/home)。采用整幅Sentinel-2A影像作为高植被覆盖区总体研究范围(约11 100 km2，图1)，其中Hyperion数据所覆盖的重合区域为重点研究区[13](约337 km2，图1黄框区)。

图1 本研究高植被覆盖区范围(Sentinel-2A数据)Fig.1 High vegetation coverage area in this study (Sentinel-2A data)

为了保证蚀变矿物信息提取的精度，使用ENVI 5.5软件对原始数据进行预处理，包括辐射定标、大气校正、波段融合以及实验范围裁剪。

表1 Sentinel-2多光谱传感器波段信息

2 蚀变信息提取方法

对于已知矿物类别的地区，可采用光谱角匹配和支持向量机方法直接提取蚀变信息。但对基础地质工作少、地面作业困难且有植被因素影响的地区，更多是使用比值与主成分分析方法来进行蚀变信息的识别与提取。

将植被抑制后的影像结合主成分分析作散点图，提取蚀变矿物信息。主成分变换是目前应用最多的特征压缩和特征选择的经典方法之一[17]。对于高植被覆盖区蚀变矿物信息提取，前提工作是去除植被干扰。刘彦丽[13]通过对比分析植被与土壤/矿物光谱曲线特征差异，提出了对波段做差分运算的差分植被抑制方法。另外，ENVI软件中也提供了一个简单快速的植被抑制工具Vegetation Suppression。

2.1 植被抑制方法

2.1.1 多光谱差分植被抑制法

差分植被抑制法是指结合特定传感器波段设置，选择植被反射率相同而土壤/矿物反射率差异较大的波段组合，将波段组合的反射率做差分计算，从而对影像中的植被信息进行抑制，为通过主成分分析散点图提取蚀变信息奠定基础。在高光谱植被差分抑制中，由于传感器具有较高光谱分辨率，很多波段组合都符合差分植被抑制法的要求，如1 336～742 nm、2 072～548 nm、2 314～701 nm、1 699～721 nm、2 203～681 nm、2 183～671 nm等。由于波段组合较多，只选取部分最优差分波段组合进行展示，如图2(a)所示。

图2 植被抑制差分波段选择示意图Fig.2 The schematic diagram of the selection of differential bands for vegetation suppression

值得注意的是，经差分植被抑制法处理之后，还需要对差分后的波段进行主成分分析并绘制散点图，才能提取蚀变矿物信息，因此至少需要3对差分波段组合才能满足多光谱差分植被抑制法的波段设置要求。对于最优化差分波段组合来讲，植被光谱在一对组合内的两个波段中应具有相近的反射率。将实验区Sentinel-2多光谱传感器的植被和土壤/矿物光谱曲线进行对比，如图2(b)所示。可以看到，Sentinel-2传感器具备3对波段组合满足最优化差分波段设置条件，从而可以通过MDVS处理得到3对差分波段(664～496 nm，2 202～560 nm，945～783 nm)。相比之下，Landsat 7/8、ASTER等其他多光谱卫星传感器不满足最优化差分波段设置要求(表2)。在实际应用中，也可以结合不同传感器波段设置，选取反射率接近的波段组合进行差分处理。本研究中，采用波段设置条件最优的Sentinel-2传感器数据进行分析和处理。

表2 常用多光谱传感器植被抑制差分波段选取结果

2.1.2 ENVI植被抑制工具

ENVI软件中提供了植被抑制工具Vegetation Suppression。该工具利用遥感影像的红光波段和近红外波段消除或削弱多光谱和高光谱影像中的植被光谱信息，对图像进行植被抑制[18]。利用植被抑制工具Vegetation Suppression对影像进行植被抑制后，波段数不会改变。研究使用的植被抑制工具是ENVI 5.5 版本中的Vegetation Suppression。

2.2 蚀变信息提取方法

对MDVS法植被抑制后的结果再进行主成分变换，分离植被与蚀变信息。选择贡献值符号相反且绝对值最大的两个主成分分量，做2D散点图并圈出异常范围。散点图中所圈定的异常范围(蚀变信息)分布在遥感影像中，经过波谱分析确定蚀变矿物类型。

对于ENVI植被抑制工具处理得到的结果，羟基蚀变矿物在Sentinel-2A影像的波长1 600 nm附近(Band 11)有强反射峰，而在波长2 200 nm附近(Band 12)有强吸收谷[19]。因此，选择Sentinel-2A数据的Band 2、4、11、12这四个波段组合进行主成分分析，提取羟基蚀变信息。Sentinel-2A数据中铁离子的光谱特征表现为：在Band 3的反射率最低，在Band 4的反射率相对Band 3高，在Band 8a的反射率下降且呈吸收谷特征，在Band 11的反射率快速上升形成反射峰，随后反射率又逐渐下降[20]。因此，提取铁染蚀变异常信息选取Sentinel-2A 数据的Band 3、4、8a、11这四个波段进行主成分分析。主成分变换后，同样使用贡献值符号相反且绝对值最大的两个主分量做散点图，最后圈定出蚀变矿物信息。

2.3 蚀变信息种类匹配及精度验证

蚀变信息提取后，对异常光谱进行匹配。光谱匹配的常用方法有光谱角匹配法(SAM)、光谱特征拟合法(SFF)以及二维编码法等。本研究使用SFF法，该方法的特点是对选定的匹配光谱特征区间内的光谱整体波形、吸收位置及吸收深度作对比，从而判断光谱的匹配程度[21]。经SFF光谱匹配后，识别出各异常代表的蚀变类型。

重点研究区地面验证点的矿物类型，是通过在研究区采集岩石样品、实验室切片、镜下分析鉴定而得到的(表3)。将影像上提取的蚀变信息光谱与USGS光谱库中光谱进行匹配，结合表3已知验证点位的矿物种类进行提取精度的验证。重点研究区的地面验证点位分布如图3所示。

表3 重点研究区地面验证点矿物种类[13]

图3 重点研究区地面验证点位分布图Fig.3 Distribution map of field verification points in key experimental areas

3 数据处理结果与分析

3.1 植被抑制结果和分析

多光谱差分植被抑制共选择3对波段，将包含植被光谱特征的绿峰和红谷附近对应的2对波段(2 202～560 nm)和(664～496 nm)做散点图；选择Band 12(2 202 nm)和Band4(664 nm)做ENVI植被抑制工具的植被抑制后散点图；将两幅散点图与原始影像散点图做直观对比，分析MDVS法与ENVI植被抑制工具法的植被抑制效果(图4)。

红圈为蚀变异常信息图4 植被抑制前后散点图对比Fig.4 Comparison of scatter plots before and after vegetation suppression

结果表明，原始影像的散点图呈扁椭圆形状，植被信息和矿物蚀变信息混杂一起，总体表现为植被背景信息。而植被反射率相同且土壤/矿物反射率差异大的两个波段做差分计算后，所求出的新波段间散点图相比于原始影像散点图，其异常信息在数据空间分布上突出明显，范围更大(图中红色虚线所圈范围)。ENVI植被抑制工具进行植被抑制的散点图相比原始数据，在数据空间分布上变化较小，但也有一小部分的异常范围，抑制植被有一定的效果。与ENVI植被抑制工具相比，MDVS方法得到的散点图中的异常信息更突出且更多，植被抑制效果更好。

3.2 蚀变信息提取结果

3.2.1 MDVS蚀变信息提取结果

利用MDVS方法将Sentinel-2数据选出的3对波段做差分，并计算各差分波段的主成分特征值，结果见表4。

表4 MDVS处理后差分波段的主成分特征值Tab.4 Principal component eigenvalues of the differential bands after MDVS processing

从表4可以看出，主成分1(PC1)与主成分2(PC2)特征值绝对值较大，且差分波段1和差分波段3特征值正负号相反，能够较好地反映不同差分波段的信息，因此适合选择这两个主成分进行散点图分析。

在理想状况下，散点图形状呈三角形，拐点分布在其3个顶点。在实际情况中，散点图形状并非呈三角形，其拐点为散点图周围的逸散部分，而蚀变异常信息往往分布在这些逸散部分。因此，在PC1和PC2的二维散点图中，先圈定外层的信息设为一级异常，再圈定内层的信息设为二级异常，分别求一级异常和二级异常的平均光谱，然后对这两个光谱进行光谱匹配，根据匹配分数确定是否为同一级异常，最后逐渐往更内层圈定。这样可将不同位置的逸散部位圈定出来，并使用不同颜色标注，即可得到不同类型的蚀变异常信息，如图5所示。

红、绿、蓝3种颜色代表光谱不同的异常信息图5 MDVS结果主成分散点图Fig.5 The scatter plot of principal components of MDVS results

将从散点图圈定的不同蚀变异常信息标记到遥感影像中，可以得到蚀变信息分布图(图6)。蚀变信息提取结果图中的红、绿、蓝颜色，与相应散点图中的异常信息颜色分别对应。为了进一步确定不同颜色区域代表的蚀变异常种类，需要对各个区域样本光谱与已有光谱库进行匹配，对其矿物类型进行识别。

红、绿、蓝3种颜色代表光谱不同的异常信息图6 重点研究区MDVS蚀变信息提取结果Fig.6 Results of alteration information extraction using MDVS method in the key experimental area

通过与地面验证信息进行对比，发现地面验证点位1、2、3恰好分别位于蚀变信息分布图中红、绿、蓝色3种不同异常信息范围内。选取验证点位的Sentinel-2A光谱与USGS矿物光谱库中的所有光谱重采样到相同分辨率的光谱，再进行SFF法光谱匹配，确定了3种颜色的异常信息所对应的蚀变矿物类型，结果如图7所示。从图7可以看到，验证点卫星数据光谱匹配结果和USGS光谱库重采样光谱特征一致，匹配效果非常好。将光谱匹配结果和表2中地面考察结果进行对比发现，点1、点2提取的矿物种类识别准确，点3提取结果相对较差，这可能与卫星空间分辨率和地面考察在尺度上有很大差异有关。

图7 验证点光谱与USGS光谱库重采样 光谱匹配结果对比Fig.7 Spectra at field verification points and the matched spectra in the resampled USGS spectral library

黏土类矿物蚀变含有羟基离子，铁染蚀变含有铁离子，这些离子在遥感图像上表现出独特的光谱特征，从而使蚀变矿物区别于非蚀变矿物。在本文的重点研究区围岩蚀变类型中，绢云母化和黏土化(高岭土化)都含有羟基团，可归类为羟基蚀变；而角闪石是铁的还原化物，归类为铁染蚀变。图6中提取结果的绿色和蓝色异常信息归为羟基蚀变，红色异常信息为铁染蚀变。因此，在多光谱影像上最终的蚀变信息可归为羟基蚀变和铁染蚀变两大类，如图8所示。

图8 重点研究区MDVS法蚀变信息提取结果Fig.8 Results of alteration information extraction by MDVS method in key experimental areas

3.2.2 ENVI植被抑制蚀变信息提取结果

根据羟基蚀变和铁染蚀变矿物在Sentinel-2A影像上不同波段具有不同的光谱吸收特征，选择具有对应波谱特征的波段做主成分变换，其结果见表5和表6。

表5 羟基蚀变提取波段的主成分特征贡献值Tab.5 Principal component eigenvectors of hydroxyl alteration extraction bands

表6 铁染蚀变提取波段的主成分特征贡献值

从表5可以看出，PC2满足羟基蚀变信息主成分分量选择标准(Band 11波段和Band 12波段特征值正负号相反，且贡献率最大)。对PC2进行密度分割提取羟基蚀变信息，如图9(a)图像中的蓝色部分。

同理，根据表6可以确定铁染蚀变信息主成分分量，只有PC2满足铁染蚀变信息主成分分量选择标准(Band 8a波段与Band 11波段特征值正负号相反，且贡献率最大)。对PC2进行密度分割提取铁染蚀变信息，如图9(b)图像中的红色区域。

图9 重点研究区ENVI植被抑制工具蚀变信息提取结果Fig.9 Results of alteration information extraction by ENVI vegetation suppression tool method in key experimental areas

为了更直观对比两种不同方法提取的重点研究区蚀变矿物信息结果，现将提取的羟基蚀变和铁染蚀变信息叠加进行考察，如图10所示。

图10 重点研究区蚀变信息叠加结果Fig.10 Overlapping results of alteration information in the key experimental area

利用ENVI植被抑制工具Vegetation Suppression，根据影像的红波段和近红外波段信息对Sentinel-2A影像进行植被抑制，然后选取特征波段做主成分变换，提取蚀变信息，如图10(b)所示。可以看到，蚀变信息的分布位置与差分植被抑制后提取的蚀变信息分布较一致，点1、2和3也包含差分植被抑制法得到的蚀变信息，但信息范围较小。用Sentinel-2A提取的蚀变信息与已知的实验区Hyperion影像提取结果(图11)作对比，可以看出，蚀变信息大部分被提取出来，并且MDVS法进行植被抑制后提取的蚀变范围和位置更准确，信息量也丰富。

图11 重点研究区Hyperion影像矿物 信息提取结果[13]Fig.11 Extraction of mineral information from Hyperion image in key experimental areas[13]

为了进一步验证Sentinel-2A数据和MDVS植被抑制方法的适用性，对包括实验区在内的更大范围研究区影像进行了蚀变矿物信息提取，结果如图12所示。从图12可以看到，利用反射率波段做差分的MDVS植被抑制方法，能够在比Hyperion数据幅宽更大的多光谱影像中进行更大空间范围的蚀变信息提取，且结果与高光谱数据相一致，这也证明了该方法在多光谱影像上的适用性和Sentinel-2数据提取蚀变信息在幅宽上的优势。

黄色虚线框内范围为重点研究区图12 高植被覆盖区MDVS蚀变信息提取叠加结果Fig.12 Overlapping result of the alteration information extraction in the high vegetation coverage area using MDVS

4 结 论

针对大范围高植被覆盖区遥感找矿面临的技术难题，在差分法植被抑制原理的基础上，结合多光谱传感器的波段设置进行差分波段选择，取得了以下成果：

(1) 基于Sentinel-2传感器波段设置特点，提出了多光谱差分植被抑制法(MDVS)，为实现高植被覆盖区蚀变矿物信息提取由高光谱数据到多光谱数据处理推广，提供了重要的技术前提。

(2) 典型实验区Sentinel-2A影像数据分析结果表明，MDVS法可以取得比已有方法更好的植被抑制效果，且基于该方法处理后的数据能够有效提取蚀变信息。

(3) 获取了大范围高植被覆盖区的蚀变信息分布图，实现了高植被覆盖区蚀变矿物信息提取推广到多光谱传感器这一重要突破，对于遥感找矿技术的发展具有重要意义。

此次实验中的高植被覆盖区影像的蚀变信息混合了大量植被信息，光谱波形和纯净矿物有很大差异。如何在提取蚀变异常大类的基础上进一步匹配识别具体矿物种类，仍有待进一步研究。另外，村庄、道路等信息在散点图中也容易被识别为异常信息，如何将其与蚀变信息进行区分需要展开更多实验分析。总之，在MDVS方法的基础上，如何进一步改善高植被覆盖区蚀变信息提取的精度，将是下一步研究工作的重点。