基于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蝕變帶，具有較大找礦潛力,為普朗銅礦的找礦勘查提供了參考。