基于ASTER与Sentinel-2A融合数据的云南普朗铜矿化蚀变信息提取

陈 琪，赵志芳，姜琦刚，夏既胜，孙 涛，曾诗卉

(1.云南大学地球科学学院，云南昆明 650500；2.云南大学国际河流与生态安全研究院，云南昆明 650500；3.云南省高校国产高分卫星遥感地质工程研究中心，云南昆明 650500；4.国土资源部三江成矿作用与资源勘查利用重点实验室，云南昆明 650051；5.吉林大学地球探测科学与技术学院，吉林长春 130000)

0 引言

铜矿是《全国矿产资源规划(2016-2020年)》列出的战略性金属矿产，是国家可持续发展的重要资源。云南普朗铜矿区是我国铜矿资源找矿勘查的主要潜力地段之一，该区位于滇西北多云多雨山地高原区，山高谷深、地质条件复杂，采用传统地面调查方法开展找矿勘查工作具有较大的局限性，而不受地面自然条件制约，且具有幅宽范围大、数据采集快速、地物信息丰富等特点的遥感技术，可在普朗铜矿资源勘查中发挥重要作用。

普朗铜矿蚀变分带特征明显，从内向外依次分布为钾化硅化带(KSi)、绢英岩化带(SiSe)以及青磐岩化带(ChEp)(范玉华和李文昌，2006；刘学龙等，2013；王凯等，2016；Liu et al.,2016)，其中，ChEp是近矿标志，KSi、SiSe及其过渡地带是高品位斑岩型铜矿找矿的重要区域，该特征可为普朗铜矿资源量勘查工作提供重要指导。

利用遥感技术进行矿化蚀变信息识别已十分成熟(Rowan et al.，1977；Loughlin，1991；王日冬和邢立新，2000；吕凤军等，2006；Gabr et al.,2015；Zhang et al.,2016；温利刚等，2017；宿虎等，2020；Traore et al.,2020)。在大量学者的研究中，遥感数据源主要采用的是Landsat-8(汪子义等，2016；贺金鑫等，2019)、ASTER(武慧智等，2019)和WorldView-3等(Sun et al.,2017；Touba and Majid,2018)；矿化遥感蚀变信息识别方法主要为主成分分析法(PCA)(张玉君等，2007；赵志芳等，2012)、比值法(BR)和光谱角法(SAM)(邓婕等，2017；Meer et al.,2012)等。然而，在遥感数据源上，Landsat-8、ASTER多光谱数据虽可免费下载，但其空间分辨率较低，而WorldView-3多光谱数据虽具有较高的空间分辨率，但其价格昂贵。在矿化遥感蚀变信息提取方法上，传统方法在植被覆盖区识别效率较低，难以提取与背景信息混合的弱蚀变信息；在研究区域上，高植被覆盖区亦是找矿勘查的重要潜力地带。

鉴于此，本研究以植被覆盖较厚的普朗矿区为研究区，选取在短波红外波段具有多个波段设置的ASTER数据和在可见光近红外波段具有更多波段设置且有更高空间分辨率的Sentinel-2A数据，在通过光谱协同统一两种数据地表反射率的基础上进行图像融合处理，以获取能兼顾高空间分辨率和高光谱保真度的ASTER-Sentinel-2A融合数据集；进而利用“主成分分析方法(PCA)+ 能谱-面积(S-A)法”进行普朗铜矿区矿化遥感蚀变信息识别。通过深入分析及野外验证，以期提出普朗铜矿蚀变分带及找矿潜力地带新认识，为普朗铜矿找矿勘查提供支撑。

1 研究区概况和遥感数据源

1.1 研究区概况

位于滇西北的中甸地区分布有大量斑岩型的铜金矿床，形成了重要的铜多金属成矿带，是我国铜矿资源找矿勘查的重要潜力地带(杨岳清等，2002；Hou et al.,2007;刘学龙等，2013；李文昌等,2013)。普朗铜矿是该地区最大的斑岩矿床，位于义敦岛弧带的南部(图1内插图)，其通常被认为是甘孜-理塘洋壳碰撞所形成(Deng et al.,2014)。

图1 普朗铜矿区地质图(据刘学龙等,2013 修改)Fig.1 Geological map of the Pulang copper mine (modified from Liu et al.,2013)Ⅰ-扬子板块；Ⅱ-甘孜-理塘寺板块结合带；Ⅲ-义敦岛弧带；Ⅳ-中咱微陆块；Ⅴ-金沙江结合带；Ⅵ-江达-维西火山弧；Ⅶ-昌都-兰坪陆块；Ⅷ-三达山-景洪火山弧；Ⅸ-澜沧江结合带；Ⅹ-保山地块；1-第四系；2-图姆沟组；3-石英二长斑岩；4-石英闪长玢岩；5-花 岗闪长斑岩；6-角岩；7-钾化带；8-绢英岩化带；9-青磐岩化带；10-蚀变界线；11-矿体Ⅰ-Yangtze Plate;Ⅱ-Ganzi-Litang plate junction belt;Ⅲ-Yidun island arc belt;Ⅳ-Zhongzan microblock;Ⅴ-Jinsha River junction zone;Ⅵ-Jiangda-Weixi volcanic arc;Ⅶ-Changdu-Lanping Block;Ⅷ-Sandashan-Jinghong volcanic arc;Ⅸ-Lancang River junction belt;Ⅹ-Baoshan Block；1-Quaternary;2-Tumugou Formation;3-quartz monzonite porphyry;4-quartz diorite porphyrite;5-granodiorite porphyry;6-hornblende; 7-potassium zone;8-sericitization zone;9-porpylitic zone;10-alteration boundary line;11-orebody

普朗铜矿的主要成矿岩体为普朗复式斑岩体，其周边分布地层为图姆沟组(图1)。区内岩浆岩主要包括三类：第一类为石英二长斑岩，根据U-Pb测年结果，其年龄约为212.8 ± 1.9 Ma (刘学龙等,2013)；第二类为石英闪长玢岩，其年龄约为219.6 ± 3.5 Ma (庞振山等,2008)，第三类为花岗闪长斑岩，其年龄约为206.3±0.7 Ma(刘学龙等,2013)。其中，石英二长斑岩中铜含量要明显高于其它斑岩，通常认为石英二长斑岩与铜矿成矿关系最为紧密(曾普胜等,2004)。

普朗铜矿从内到外蚀变分带特征较为明显，各蚀变带均发育有典型的指示性特征蚀变矿物。本研究选取了石英、绢云母、绿泥石、绿帘石四种特征矿物来研究普朗铜矿区的矿化遥感蚀变分带特征。

1.2 遥感数据源

1.2.1 ASTER数据

作为美国宇航局实施的地球观测系统计划，ASTER在1999年12月18日成功发射升空，成像幅宽为60 km*60 km。ASTER数据拥有3个谱段共计14个波段，分别为:(1) 3个空间分辨率为15 m的可见光和近红外波段(Band 1～Band 3)，其波长区间设置为0.52～0.86 μm;(2) 6个空间分辨率为30 m的短波红外波段(Band 4～Band 9)，其波长区间设置为1.6～2.43 μm;(3) 5个空间分辨率为90 m的热红外波段(Band 10～Band 14)，其波长区间设置为8.125～11.65 μm。研究区涉及1景ASTER数据，其获取时间为2006年12月22日。

1.2.2 Sentinel-2A数据

Sentinel-2A作为“哥白尼计划”中的首颗多光谱成像卫星，在2015年6月23日成功发射升空，成像幅宽为290 km。Sentinel-2A数据拥有3个谱段共计13个波段，分别为:(1) 可见光谱段内(波长范围为0.433～0.747 μm)设置有6个波段，空间分辨率为10 m(Band 2、Band 3、Band 4)、20 m(Band 5、Band 6)和60 m(Band 1);(2)近红外谱段内(波长范围为0.773～1.39 μm)设置有5个波段，空间分辨率为10 m(Band 8)、20 m(Band 7、Band 8A)、60 m(Band 9、Band 10);(3)短波红外谱段内(波长范围为1.565～2.28 μm)设置有2个波段，空间分辨率为20 m(Band 11、Band 12)。研究区涉及1景Sentinel-2A数据，其获取时间为2016年12月4日。

2 遥感图像融合和热液蚀变信息提取方法

2.1 遥感图像融合方法

2.1.1 光谱协同

ASTER数据与Sentinel-2A数据由不同卫星不同传感器拍摄获取，其波谱响应函数各不相同，若将两种数据直接进行图像融合，势必会影响融合数据的光谱保真度。因此，本研究在图像融合之前开展了ASTER数据与Sentinel-2A数据的光谱协同，使两种数据的地表反射率保持在同一尺度上。

对比分析两种数据的波段设置发现，ASTER数据的Band 1、Band 2、Band 3、Band4等4个波段与Sentinel-2A数据的Band 3、Band 4、Band 8、Band 11等4个波段的波长范围设置重叠，其余波段的波长范围设置不重叠。因此，本研究将以上两种数据的光谱协同分为两步：(1)第一步是波长重叠波段的光谱协同，因同名地物应具有相同像元值，光谱协同处理即转化为两种数据像元值的回归拟合处理，通过构建回归方程将自变量的遥感数据像元值统一到因变量遥感数据像元值的同一尺度上；(2)第二步是波长不重叠波段的光谱协同，可通过计算其与波长不重叠波段的相关性，选取与其相关性最大的波长不重叠波段所构建的回归方程进行光谱协同。

2.1.2 遥感图像融合

本研究采用的融合方法为高通滤波法(HPF)(Ghassemian,2016)。其原理为：分别对高空间分辨率的全色数据和低空间分辨率的多光谱数据进行高通滤波和低通滤波，从而增强全色数据的高频信息(空间细节信息)和多光谱数据的低频信息(光谱信息)；在此基础上，进一步再将两种数据加权求和，以获取最终的融合数据，使其不仅具有高空间分辨率，亦具有高光谱保真度 (图2)。

图2 基于高通滤波的融合方法示意图Fig.2 Schematic diagram of image fusion based on high pass filtering (HPF)

2.2 矿化蚀变提取方法

2.2.1 比值+主成分分析法

比值法是代数运算的原理，不同蚀变矿物的波谱曲线存在不同的吸收谷、反射峰特征，通过反射峰位置与吸收谷位置反射率值得比值运算，可有效增强各地质信息之间的差异，从而从而达到提取蚀变信息的目的。比值法是矿化蚀变信息增强最常用、也是非常有效的一种方法。

主成分分析法在多光谱遥感数据矿化遥感蚀变填图中已被广泛运用，其本质是对数据进行降维处理。多光谱遥感影像设置有多个波段，通过对其开展主成分变换，可使得重复冗余原各波段转换为新的互不相干的波段集，从而各地物信息独立地分布在新的各主成分分量中，这就使主成分分析具有了分离信息的作用。

本研究将比值法与主成分分析法结合运用，充分发挥其各自的优势，用于提取石英蚀变矿物，其方法流程为：(1)通过比值运算计算石英指数(QI)、碳酸盐指数(CI)、铁镁质指数(MI)(Boardman,1992)；(2)将QI、CI、MI作为主成分分析的三个分量进行变换，根据特征向量的贡献值，选取主分量提取石英蚀变矿物(Pour et al.,2019)。

2.2.2 主成分分析+能谱-面积法

能谱-面积法是多重分形理论为应用。矿化蚀变与成矿作用相伴发生，矿化蚀变信息的记录数据是非线性结构的。此类数据在空间域上往往表现为局部不均一性，但其在频率域空间中却具有多重分形分布的特点。因此，本研究拟基于多重分形分布的特点进行矿化遥感蚀变信息识别。其关系式表达如下：

A(>S′)∞S-β

(1)

式(1)中，S可视为频率域空间中拟提取的蚀变信息所在主分量的能谱密度，A是大于某一能谱密度值(S0)所对应的面积，β值表示为分形特征。对S于A进行双对数变换，在双对数(log-log)图上，会形成不同的两条线段，其表征为蚀变信息与背景信息不同的分形关系，两条线段的交点即为蚀变信息与背景信息的分割阈值。通过该阈值即可实现蚀变信息与背景信息的分离(Chen et al.,2018;Chen et al.,2019)。

本研究通过主成分分析将原波段重新组合为新的组分，使蚀变矿物信息独立地分布在新的某一主成分分量中，进而对该分量运用S-A法，以增强提取蚀变矿物信息。

3 结果与分析

3.1 遥感图像融合结果

3.1.1 光谱协同

(1)波长重叠波段光谱协同

本研究采用Sentinel-2A数据Band 3与ASTER数据Band 1进行光谱协同示范，将Sentinel-2A数据Band 3像元数值为自变量n，ASTER数据Band 1对应像元数值为因变量m，开展一元线性回归拟合，其拟合结果如下：

m=1.369n+1481.5

(2)

基于该回归方程对Sentinel-2A数据Band 3进行光谱协同，统计结果发现(表1)，原始Sentienl-2A数据Band 3的最小值(Min)、最大值(Max)、均值(Mean)、标准差(SD)等各项指标，与ASTER数据Band 1存在明显的尺度差异，而通过光谱协同后，两种数据的各项指标在数值尺度上基本趋于统一。同理，对Sentinel-2A数据Band 4、Band 8、Band 11进行了一元线性回归拟合处理，实现了两种影像数据在波长重叠波段的地表反射率统一。

表1 原始ASTER数据Band 1与光谱协同前后Sentinel-2A数据Band 3参数统计表

(2)波长不重叠波段光谱协同

Sentinel-2A数据剩下的不重叠波段，其光谱协同则需采用波长重叠波段所构建的一元线性回归方程。回归方程的选取原则为波段间相关系数最高，以减少系统误差。通过相关性计算，Band 1和Band 2光谱协同选择Band 3的回归方程；Band 5光谱协同选择Band 4的回归方程；Band 6、Band 7、Band 8A、Band 9光谱协同选择Band 8的回归方程；Band 10、Band 12光谱协同选择Band 11的回归方程。

3.1.2 遥感图像融合

(1)Sentinel-2A数据不同分辨率波段的图像融合

在光谱协同统一地表反射率后，进而开展图像融合处理。Sentinel-2A数据在Band 1，Band 5，Band 6，Band 7，Band 8A，Band 9，Band 10，Band 11，Band 12具有较低的空间分辨率，本研究通过将其与高空间分辨率Band 2，Band 3，Band 4，Band 8进行图像融合处理以提高其空间分辨率。

通过相关性计算，确定具有最大相关性作为对应融合波段的选取原则，其波段选取详细情况见表2。图3为研究区Sentinel-2A数据Band 5，Band 6，Band 7波段经HPF融合前后效果对比图。由图可知，图像融合后在地物边界、空间纹理细节等方面均达到了信息增强的效果，且较好地保持了光谱保真度。同理，采用该融合方法对Sentinel-2A数据Band 1，Band 8A，Band 9，Band 10，Band 11，Band 12进行了图像融合处理。

表2 Sentinel-2A图像融合处理波段对应表

图3 原始Sentinel-2A图像与基于HPF融合后的对比图Fig.3 Comparison of original Sentinel-2A image and image after fusion based on HPF methoda-原始数据，R(5)、G(6)、B(7)波段组合(20m)；b-HPF融合数据，R(5)、G(6)、B(7)波段组合(10m)a-original data，R(5),G(6)and B(7)band combination(20m)；b-HPF fusion data，R(5)，G(6)and B(7)band combination(10m)

(2)ASTER数据与Sentinel-2A数据的图像融合

本研究采用ASTER数据的所有波段与Sentinel-2A数据Band 2，Band 3，Band 4，Band 8进行融合，以提高ASTER数据空间分辨率至10m。通常波段设置间隔越近其相关性越大，基于此选取了两种图像融合处理的对应波段，详见表3。

表3 ASTER数据与Sentinel-2A数据进行融合处理的波段对应表

图4为研究区ASTER数据Band 4，Band 5，Band 6分别作为R、G、B通道进行彩色合成的图像融合前后对比图。图像融合后不仅在地物边界、空间纹理细节等方面提高了信息强度，同时较好地保持了光谱保真度。同理，基于以上图像融合方法对ASTER数据剩余的Band 1，Band 2，Band 3，Band 7，Band 8，Band 9进行了图像融合处理。

图4 ASTER数据Band 4，Band 5，Band 6彩色合成图像HPF融合前后对比图Fig.4 Comparison of original ASTER image and image after HPF fusiona-原始数据，R(4)、G(5)、B(6)波段组合(30 m)；b-HPF融合数据，R(4)、G(5)、B(6)波段组合(10 m)a-original data，R(4),G(5)and B(6)band combination(30 m)；b-HPF fusion data，R(4),G(5)and B(6)band combination(10 m)

基于以上处理，获取了最终具有高空间分辨率及高光谱保真度的ASTER-Sentinel-2A融合数据集。

3.2 矿化遥感蚀变信息提取结果

3.2.1 比值+主成分分析法提取结果

利用多光谱遥感数据开展矿化遥感蚀变信息识别的理论基础是蚀变矿物在遥感数据上的波段响应特征。根据石英标准波谱曲线在ASTER数据上的波谱响应特征(图5)，石英在ASTER数据的热红外波段Band 10(简称A10，下同)和A12有强反射峰，在ASTER数据的A11有强吸收峰，可通过比值运算计算QI(Band 11*Band 11/Band 10*Band 12)值，增强石英特征；同时根据文献(Boardman,1992)基于比值运算计算了CI(Band 13/Band 14)、MI(Band 12/Band 13)值。

对QI、CI、MI进行主成分变换，根据特征向量矩阵(表4)，PC2分量中QI(0.725769)为正且贡献值较大，CI(-0.178896)、MI(-0.664270)均为负，因此PC2分量中的亮像素可表示为石英(图7)。

表4 ASTER热红外波段比值运算-主成分分析的特征向量矩阵

3.2.2 主成分分析+能谱-面积法提取结果

根据绢云母、绿泥石、绿帘石等特征矿物的标准波谱曲线在ASTER-Sentinel-2A融合数据集上的光谱响应特征(图5)，绢云母在ASTER-Sentinel-2A融合数据的Sentinel-2A Band1(简称S1，下同)和ASTER Band 6(简称A6，下同)均有强吸收峰，在ASTER-Sentinel-2A融合数据的S8和A7均有明显的反射峰，因此确定绢云母的特征响应光谱为S1、S8、A6、A7。同理，确定绿泥石的特征响应光谱为 S3、S4、A5、A8；确定绿帘石的特征响应光谱为A1、S8、A5、A9。

图5 蚀变矿物的标准波谱曲线(USGS波谱库)与遥感数 据波段对应图Fig.5 Spectral curves of alteration minerals (from the United States Geological Survey (USGS)) corresponding to bands of remote sensing data

基于上述分析确定的各特征蚀变矿物的诊断性波段进行主成分变换，根据主成分变换结果中各分量的贡献，确定矿化遥感蚀变信息所处分量。其中，绢云母、绿泥石、绿帘石所处的分量均为PC4。将以上确定的主分量从空间域经过傅里叶变换，转换至频率域，应用S-A法生成双对数图(图6)。在双对数图中可拟合两条线段，其分别代表矿化遥感蚀变信息和背景信息，而在两线段相交位置所对应的能谱密度值即为分割阈值。由图可知，绢云母、绿泥石、绿帘石与背景信息的分割阈值分别为0.75、1.15、1.4。基于该阈值进一步开展傅里叶逆变换，即可识别相应的矿化遥感蚀变信息(图7)。

图6 基于“PCA+S-A法”提取蚀变矿物的logS-logA图Fig.6 logS-logA diagram of alteration minerals extracted by PCA+S-A method

图7 基于遥感数据提取的蚀变矿物的验证与分析图Fig.7 Verification and analysis of alteration minerals extracted from remote sensing data1-石英；2-绢云母；3-绿泥石；4-绿帘石；5-原蚀变界线；6-新划分蚀变界线；7-野外验证点；8-野外验证路线；9-钻孔位置1-quartz；2-sericite；3-chlorite；4-epidote；5-original alteration boundary；6-newly redivided alteration boundary；7-field confirmation point； 8-field inspection route；9-drilling hole

针对矿化遥感蚀变信息提取结果，2018年8月前往普朗斑岩型铜矿区开展了野外实地验证工作，采集了普朗斑岩型铜矿中位于钾化硅化带、绢英岩化带、青磐岩化带的86个手标本，并进行了显微镜下岩矿鉴定。通过验证发现，蚀变分带特征矿物的遥感识别精度较高，总体精度超过80%(表5)。

表5 典型点矿化遥感蚀变信息提取及野外验证一览表

3.3 蚀变矿物分布特征及找矿潜力地段分析

图7中黑色线条是原蚀变填图工作中圈定的蚀变分带边界线，通过将其与本次识别的矿化遥感蚀变信息叠加分析可以发现，位于西部的普朗铜矿主采区，本次识别的石英、绢云母、绿泥石、绿帘石信息的分布基本符合原蚀变填图工作圈定的蚀变分带。其中，石英主要分布在矿体中心，绢云母主要分布在矿体边缘，绿泥石、绿帘石主要分布在矿体外围，少量绿泥石、绿帘石分布在矿体上，已有研究发现矿化区绿泥石、绿帘石的叠加有利于增强铜矿化(Cao et al.,2019)。然而，在原蚀变填图工作圈定为ChEp分带的东部区域，本次研究工作识别有大量表征钾化硅化和绢英岩化的石英、绢云母信息，与以往认识不一致。结合绢云母和石英表征SiSe，绿泥石和绿帘石表征ChEp，石英及KSi通常位于ChEp、SiSe中心的特征，对东部区域重新划分了蚀变分带。同时，再进一步根据KSi、SiSe及其过渡地带是高品位斑岩型铜矿找矿的重要标志，推测东部区域重新划分的KSi、SiSe具有较好的找矿潜力。

针对圈定的找矿潜力地带，收集了4个钻孔资料。通过分析发现，4个钻孔所处位置均存在有铜矿资源量，铜品位约为0.1%～0.4%，且从钻孔Ⅰ位置到钻孔Ⅳ位置(本次划分的潜力地带边缘向中心靠近)，其铜矿资源量富集情况越来越好，可进一步印证该找矿潜力地带具有较好地找矿前景。

4 结论

本研究通过光谱协同处理，对Sentinel-2A数据与ASTER数据的地表反射率进行了尺度上的统一，基于此进一步采用HPF融合方法对以上两种多光谱数据开展了图像融合处理，获取了兼具高空间分辨率和高光谱保真度的ASTER-Sentinel-2A融合数据；并利用“比值+主成分分析法”与“主成分分析+能谱-面积法”识别了普朗铜矿区的矿化遥感蚀变信息分布特征，通过深入分析及野外验证，对普朗铜矿有了新的认识：(1)在原蚀变填图工作中判定为ChEp蚀变带的东部区域，本次研究识别了大量表征钾化硅化和绢英岩化的石英、绢云母信息，并根据其分布特征对该区域重新进行了蚀变带的划分；(2)推测东部区域中新划分的KSi和SiSe蚀变带，具有较大找矿潜力,为普朗铜矿的找矿勘查提供了参考。