孫 雨，劉家軍，翟德高，柳振江，張方方，趙英俊，劉鵬飛，王子濤

(1.中國地質(zhì)大學(xué)(北京)地球科學(xué)與資源學(xué)院,北京 100083；2.核工業(yè)北京地質(zhì)研究院，北京 100029；3.遙感信息與圖像分析技術(shù)國家級重點(diǎn)實(shí)驗(yàn)室，北京 100029；4.中國地質(zhì)大學(xué)(北京)地質(zhì)過程與礦產(chǎn)資源國家重點(diǎn)實(shí)驗(yàn)室,北京 100083)

0 引言

高光譜技術(shù)能夠根據(jù)指紋性光譜特征進(jìn)行地物鑒別，在巖性識別和蝕變礦物提取等地質(zhì)領(lǐng)域取得了較好的應(yīng)用效果，是當(dāng)前國內(nèi)外遙感應(yīng)用研究的熱點(diǎn)之一。衛(wèi)星高光譜遙感由于視域?qū)拸V，一次成像即可獲取大面積、價格低廉的高光譜數(shù)據(jù)，經(jīng)過數(shù)據(jù)處理后能夠?yàn)榈刭|(zhì)礦產(chǎn)勘查和基礎(chǔ)地質(zhì)調(diào)查提供重要遙感信息。搭載有高光譜傳感器的衛(wèi)星主要有美國EO-1、我國的GF-5和資源一號02D衛(wèi)星，以往學(xué)者針對EO-1的Hyperion衛(wèi)星高光譜數(shù)據(jù)處理和應(yīng)用開展了大量研究工作(Hubbard and Crowley，2005；Zadeh et al.，2014；Ayoobi and Tangestani，2018；Hu et al.，2019；Souza et al.，2021)，針對我國GF-5衛(wèi)星高光譜數(shù)據(jù)也開展了一些地質(zhì)應(yīng)用研究(董新豐等，2020；Ye et al.，2020；李娜等，2021)。當(dāng)前，全世界民用在軌能夠提供短波紅外譜段衛(wèi)星高光譜數(shù)據(jù)的只有我國的資源一號02D衛(wèi)星，有關(guān)這一新型高光譜數(shù)據(jù)的地質(zhì)應(yīng)用研究還較少(李娜等，2020；魏紅艷等，2020；李根軍等，2021)，迫切需要在數(shù)據(jù)處理方法和應(yīng)用效果方面加強(qiáng)研究交流，推動“天上衛(wèi)星好用”和“地上數(shù)據(jù)用好”兩方面更好的協(xié)調(diào)發(fā)展。

本文應(yīng)用甘肅省瓜州縣頭吊泉-南大灘地區(qū)的資源一號02D高光譜數(shù)據(jù)，探索了該高光譜數(shù)據(jù)的地質(zhì)應(yīng)用前景。采用改進(jìn)的SAM填圖算法提取該地區(qū)的褐鐵礦、綠泥石、方解石、短波云母等7種蝕變礦物，分析蝕變礦物的地質(zhì)分布特征，并對新發(fā)現(xiàn)的金礦化蝕變帶進(jìn)行野外查證，證實(shí)存在地表金礦化異常，為該區(qū)找礦勘查提供了重要線索。

1 研究區(qū)地質(zhì)背景

圖1 甘肅頭吊泉-南大灘地區(qū)地質(zhì)簡圖

2 衛(wèi)星高光譜數(shù)據(jù)源

2019年9月12日11時26分，資源一號02D衛(wèi)星在我國太原衛(wèi)星發(fā)射中心成功發(fā)射升空進(jìn)入預(yù)定軌道。作為我國自主建造并成功運(yùn)行的首顆民用高光譜業(yè)務(wù)衛(wèi)星，該衛(wèi)星是國家民用空間基礎(chǔ)設(shè)施中新型對地觀測衛(wèi)星發(fā)展的又一重要成果。該星可有效獲取寬幅高光譜等遙感數(shù)據(jù)，將與后續(xù)系列衛(wèi)星組網(wǎng)形成全球領(lǐng)先的業(yè)務(wù)化對地光譜探測能力，為自然資源領(lǐng)域的衛(wèi)星高光譜數(shù)據(jù)業(yè)務(wù)化應(yīng)用提供了重要的數(shù)據(jù)保障(楊詩瑞，2019)。

資源一號02D衛(wèi)星軌道高度為778 km，搭載的高光譜傳感器載荷可獲取幅寬60 km、空間分辨率優(yōu)于30 m的高光譜(AHSI)數(shù)據(jù)。資源一號02D高光譜數(shù)據(jù)譜段共166個，其中可見近紅外譜段光譜分辨率優(yōu)于10 nm，短波紅外譜段光譜分辨率優(yōu)于20 nm。

本次使用的是1景2020年10月12日獲取的資源一號02D AHSI高光譜遙感數(shù)據(jù)，產(chǎn)品級別為L1A級，產(chǎn)品號為L1A0000179333，包含GeoTiff數(shù)據(jù)文件、軌道參數(shù)RPB文件及定標(biāo)及光譜響應(yīng)RAW文件等。頭吊泉-南大灘地區(qū)位于整景高光譜數(shù)據(jù)的西北側(cè)(圖2)，面積約450 km2。

圖2 頭吊泉-南大灘地區(qū)高光譜影像位置圖

3 高光譜數(shù)據(jù)處理方法

3.1 高光譜數(shù)據(jù)預(yù)處理

資源一號02D衛(wèi)星的L1A級高光譜原始數(shù)據(jù)為帶RPC校正參數(shù)的DN值數(shù)據(jù)，需要進(jìn)行數(shù)據(jù)預(yù)處理才能用于蝕變礦物填圖，數(shù)據(jù)預(yù)處理由基礎(chǔ)預(yù)處理和大氣校正兩部分構(gòu)成。在基礎(chǔ)預(yù)處理方面，應(yīng)用當(dāng)前熱門的python編程語言，開發(fā)了具有自主知識產(chǎn)權(quán)的資源一號02D高光譜數(shù)據(jù)批量基礎(chǔ)預(yù)處理軟件模塊，具備解壓縮、正射校正、輻射定標(biāo)、波段組合、BSQ轉(zhuǎn)BIL、波段去除功能，實(shí)現(xiàn)一鍵化對指定文件夾下所有高光譜原始數(shù)據(jù)進(jìn)行批量的基礎(chǔ)預(yù)處理。軟件模塊的主要功能函數(shù)如圖3所示。

圖3 高光譜數(shù)據(jù)基礎(chǔ)預(yù)處理模塊主要功能函數(shù)

批量基礎(chǔ)預(yù)處理軟件模塊的處理思路是首先應(yīng)用ASTER GDEM V3高程數(shù)據(jù)進(jìn)行正射校正，得到具有高斯克呂格-CGCS2000投影坐標(biāo)系統(tǒng)的高光譜數(shù)據(jù)。然后，應(yīng)用原始數(shù)據(jù)自帶的定標(biāo)參數(shù)進(jìn)行輻射定標(biāo)，得到具有物理意義的輻射亮度數(shù)據(jù)，并將VNIR和SWIR數(shù)據(jù)文件合并成具有166波段的一個數(shù)據(jù)文件，再將BSQ格式轉(zhuǎn)換為BIL格式。隨后，去除高光譜數(shù)據(jù)的水汽吸收擾動明顯的波段和目視效果差的傳感器起止波段(共24個)，最終得到142個波段的高光譜數(shù)據(jù)。

資源一號02D高光譜數(shù)據(jù)的單景tar壓縮包大小為1.31 G，應(yīng)用本次開發(fā)高光譜預(yù)處理模塊的處理時長為121 s，使用ENVI 5.3商業(yè)軟件的處理時長為579 s，本次開發(fā)的模塊在耗時上約為商業(yè)軟件的1/5，可有效助力于區(qū)域大面積的資源一號02D高光譜數(shù)據(jù)快速預(yù)處理。

在采用開發(fā)的算法模塊進(jìn)行基礎(chǔ)預(yù)處理后，再設(shè)置好各項(xiàng)參數(shù)后應(yīng)用ENVI 5.3中FLAASH算法模塊進(jìn)行大氣校正，得到可用于礦物填圖的高光譜反射率數(shù)據(jù)。

首先，是對學(xué)生進(jìn)行分層。根據(jù)學(xué)生之前的學(xué)習(xí)成績、班主任建議等因素，綜合評估，把學(xué)生分為A班，B班，原則上比例為1∶1。A班側(cè)重提高，B班側(cè)重基礎(chǔ)。學(xué)生分層工作中最重要的是心理輔導(dǎo)，要向?qū)W生宣傳正確的教育理念，疏導(dǎo)學(xué)生心理問題，避免學(xué)生出現(xiàn)自卑等負(fù)面情緒。

3.2 高光譜數(shù)據(jù)降維增值

高光譜數(shù)據(jù)擁有上百個波段，存在兩個處理關(guān)鍵問題。一是，如果將142個波段全部進(jìn)行計(jì)算將極大降低處理效率，而且蝕變礦物并不是在所有波段都具有明顯吸收特征；二是，在蝕變礦物反映敏感的短波紅外譜段(燕守勛等，2003)，其光譜采樣間隔在16 nm左右，與GF-5高光譜數(shù)據(jù)的8 nm相比光譜分辨率有所降低，需要考慮在光譜維度進(jìn)行升采樣以提高光譜分辨率從而可以更加精細(xì)的識別礦物。

針對第一個問題，本文采取波段篩選組合的高光譜數(shù)據(jù)降維處理方法，有針對性地對蝕變礦物有明顯吸收反射特征的波段進(jìn)行挑選后組合，使礦物特征更加集中(鄭向濤，2017)。對于VNIR高光譜數(shù)據(jù)，選擇499～997 nm譜段；對于SWIR高光譜數(shù)據(jù)，選擇2014～2451 nm譜段。

三次樣條函數(shù)作為數(shù)學(xué)分析中常用的插值方法，廣泛應(yīng)用于離散點(diǎn)插值計(jì)算(張成良和劉小泉，2006；劉翔舸等，2010；殷瑞娟等，2013)。楊萍等(2007)應(yīng)用三次樣條插值法根據(jù)8波段數(shù)值求解出106個波段(間隔5nm)的光譜響應(yīng)值，認(rèn)為該方法可以非常近似地模擬出光滑的地物光譜曲線。梁晨(2016)對同型號不同分辨率儀器實(shí)測數(shù)據(jù)進(jìn)行了三次樣條插值，解決了光譜測量點(diǎn)的匹配問題。因此，本文針對第二個問題，引入了三次樣條函數(shù)對短波紅外譜段在相鄰波段間進(jìn)行插值運(yùn)算，實(shí)現(xiàn)不破壞光譜固有信息的光譜維增值，將光譜采樣間隔由16 nm升采樣到了2 nm。

3.3 蝕變填圖端元選取

本文應(yīng)用降維增值后的高光譜數(shù)據(jù)，使用ENVI 5.3軟件的光譜沙漏(Spectral Hourglass)流程算法(李光輝等，2013)，得到蝕變礦物在VNIR譜段和SWIR譜段的圖像候選端元光譜。然后，基于專家光譜知識對候選端元光譜進(jìn)行優(yōu)選，得到VNIR譜段褐鐵礦的圖像端元光譜(圖4a)和SWIR譜段綠泥石、方解石、白云石、短波云母、中短波云母、中長波云母的圖像端元光譜(圖4b)。為了進(jìn)一步凸顯不同蝕變礦物的吸收和反射特征，對SWIR譜段圖像候選端元光譜進(jìn)行包絡(luò)線去除，得到該譜段用于蝕變礦物填圖的圖像端元光譜(圖5)。其中，短波云母端元第一特征吸收峰在2199 nm附近，中短波云母端元光譜第一特征吸收峰在2206 nm附近，中長波云母端元光譜第一特征吸收峰在2216 nm附近。

圖4 優(yōu)選后的蝕變礦物圖像端元光譜Fig.4 Selected image endmember spectra of alteration minerals

圖5 SWIR譜段數(shù)據(jù)包絡(luò)線去除后的蝕變礦物圖像端元光譜Fig.5 Image endmember spectra of alteration minerals in SWIR data after continuum removal

3.4 蝕變礦物填圖與制圖

在VNIR譜段，褐鐵礦等蝕變礦物的光譜吸收峰較為寬緩。針對這一特點(diǎn)，將優(yōu)選后的圖像端元光譜作為參考光譜，以各個像元的圖像光譜為待匹配光譜。在SWIR譜段，云母等蝕變礦物的光譜吸收峰較為尖銳，而且往往存在不止一個光譜吸收峰，因此將包絡(luò)線去除后的圖像端元光譜作為參考光譜，以各個像元包絡(luò)線去除后的圖像光譜為待匹配光譜。在選擇好參考光譜和待匹配光譜后，采用本文建立的改進(jìn)的光譜角填圖方法進(jìn)行蝕變礦物填圖。

改進(jìn)的光譜角填圖方法是針對傳統(tǒng)光譜角填圖(Spectral Angle Mapping)方法僅度量全局譜形相似度的缺陷，在對SAM方法增加一系列特征吸收-反射峰位判別條件后，逐像元計(jì)算待匹配光譜與參考光譜的相似度。本文應(yīng)用當(dāng)前流行的python程序開發(fā)語言，開發(fā)了改進(jìn)的光譜角填圖算法，可以批量進(jìn)行高光譜蝕變礦物填圖。

高光譜蝕變礦物填圖計(jì)算完成后，得到礦物匹配得分灰度圖。像元灰度值為弧度值，灰度值越小則與參考光譜匹配程度越高，進(jìn)行低值端閾值切割并填充對應(yīng)的顏色，得到各種單礦物分布圖，并以空間坐標(biāo)信息為鏈接，將所有蝕變礦物全部疊加到巖性-構(gòu)造地質(zhì)圖上，在添加比例尺等圖飾信息后，編制頭吊泉-南大灘高光譜遙感地質(zhì)綜合圖(圖6)。然后，應(yīng)用成像時間為2021年1月29日的亞米級(0.8 m)空間分辨率的GF-2國產(chǎn)遙感數(shù)據(jù)制作了真彩色高分影像圖，在疊加了褐鐵礦和中短波云母蝕變礦物后，編制了南大灘-西紅灘地段蝕變高清分布圖(圖7)。

圖6 頭吊泉-南大灘地區(qū)高光譜遙感地質(zhì)綜合圖Fig.6 Comprehensive map showing hyperspectral geological information of the Toudiaoquan-Nandatan area

圖7 南大灘-西紅灘地段金礦化蝕變帶分布圖Fig.7 Distribution of gold mineralization alteration belt in the Nandatan-Xihongtan section

4 結(jié)果與討論

4.1 蝕變礦物分布特征

4.2 金礦化蝕變帶展布特征

許多類型的礦產(chǎn)都與熱液蝕變關(guān)系密切，是開展礦產(chǎn)資源勘查的有利信號(朱永峰和安芳，2010，杜斌等，2021)。蝕變礦物是成礦熱液蝕變活動的直觀證據(jù)，以往多光譜遙感技術(shù)提取的羥基等礦物族級別的蝕變礦物為金屬礦產(chǎn)勘查提供了重要的遙感判據(jù)(王欽軍等，2017；Yousefi et. al，2018；劉小雨等，2020)，高光譜遙感更是提供了傳統(tǒng)地質(zhì)方法難以獲取的無縫像元尺度的蝕變單礦物信息(孫雨等，2015；連琛芹等，2020；吳德文等，2021)，為金礦、銅礦、鈾礦、鎢鉬礦等金屬礦種的找礦勘查圈定了范圍更小的找礦遠(yuǎn)景區(qū)(汪重午等，2014；賀金鑫等，2017；劉德長等，2017；葉發(fā)旺等，2021)。

眾多研究表明，金礦化與硅化、黃鐵礦化、褐鐵礦化、絹云母化、綠泥石化等蝕變關(guān)系密切。范瀟等(2015)對焦家金礦田進(jìn)行了構(gòu)造蝕變測量，蝕變巖帶的地球化學(xué)測試表明絹英巖化蝕變階段是礦化元素富集的主要階段。礦物填圖中使用的“云母”系指具有不同結(jié)構(gòu)和粒度的絹云母、白云母族礦物，通常統(tǒng)稱為“絹云母”?！敖佋颇浮钡亩嗣骟w晶體結(jié)構(gòu)主要受Tschemark替代過程影響，從而造成礦物結(jié)構(gòu)中的Al-OH吸收峰位置發(fā)生漂移，一般介于2185～2225 nm之間。邵雪維等(2021)對膠東新城金礦深部蝕變巖進(jìn)行了光譜測量，研究結(jié)果表明金礦化巖心段的絹云母Al-OH吸收峰位置在2205～2208 nm。其他礦種圍巖蝕變中的絹云母礦物Al-OH吸收峰偏移也有類似的規(guī)律。如郭娜等(2012)對西藏甲瑪斑巖銅多金屬礦進(jìn)行了蝕變礦物光譜測量，結(jié)果表明礦體主要與Al-OH吸收峰位置處于2205～2208 nm的絹云母共存。周巖等(2019)對福建紫金山銅金礦區(qū)進(jìn)行了巖心掃描，發(fā)現(xiàn)與銅鉬礦體關(guān)系密切的絹云母Al-OH吸收峰位置大于2205 nm，以此找礦標(biāo)志發(fā)現(xiàn)了肉眼難以識別的新礦體。徐清俊等(2016)對新疆白楊河鈾礦進(jìn)行了全孔巖心光譜測量與分析，證實(shí)鈾礦化與赤鐵礦化、高鋁絹云母(Al-OH吸收峰位置在2192～2203 nm)和中鋁絹云母(Al-OH吸收峰位置在2204～2212 nm)有關(guān)，三者共同存在的地段是鈾成礦有利地段。

金礦化普遍與黃鐵絹英巖化關(guān)系密切，即與黃鐵礦(以及進(jìn)一步蝕變形成的褐鐵礦)、絹云母、石英礦物密不可分。在譜段范圍400～2500 nm的資源一號02D高光譜數(shù)據(jù)上，代表硅化的石英和代表黃鐵礦化的黃鐵礦均不具有明顯光譜特征，無法直接探測識別出來；褐鐵礦、絹云母等蝕變礦物具有光譜特征，可以直接探測識別出來。因此，金礦化在高光譜蝕變圖上可表現(xiàn)為褐鐵礦+絹云母。在詳細(xì)分析高光譜蝕變礦物分布圖結(jié)合地質(zhì)巖性背景，在研究區(qū)南部識別出了1條可能的金礦化蝕變帶。該礦化蝕變帶分布于南大灘-西紅灘地段，發(fā)育褐鐵礦和中短波云母蝕變礦物(圖7)。蝕變帶總體近東西向展布，斷續(xù)延伸約15 km，所處地質(zhì)體為海西中期及晚期花崗巖。

4.3 金礦化蝕變帶驗(yàn)證情況

針對高光譜技術(shù)在南大灘-西紅灘地段識別出的金礦化蝕變帶，綜合采用野外地質(zhì)踏勘、ASD光譜測量和地球化學(xué)樣品采集及分析測試等技術(shù)手段，驗(yàn)證了該帶地表存在金礦化現(xiàn)象。

野外地質(zhì)踏勘發(fā)現(xiàn)該蝕變帶位于第四系大灘的北側(cè)，蝕變原巖為褐紅色中粗粒塊狀花崗巖，發(fā)育石英脈，可見明顯的硅化及褐鐵礦化現(xiàn)象(圖8a、8b)，絹云母化難以肉眼辨識。通過對蝕變巖的ASD光譜測量結(jié)果分析，顯示同時存在褐鐵礦和中短波云母礦物(圖8c)，褐鐵礦體現(xiàn)在實(shí)測光譜曲線中880 nm附近的寬緩吸收峰和765 nm附近的反射峰，短波云母體現(xiàn)在實(shí)測光譜曲線中2304 nm附近的第一吸收峰、2350 nm附近的第二吸收峰和2445 nm附近的第三吸收峰，不僅與USGS標(biāo)準(zhǔn)光譜庫中褐鐵礦、絹云母礦物的標(biāo)準(zhǔn)光譜曲線特征吻合，也與用于高光譜蝕變礦物提取的圖像端元光譜特征吻合較好，證實(shí)高光譜提取的蝕變礦物確實(shí)存在，高光譜填圖結(jié)果可信。

圖8 南大灘-西紅灘地段金礦化蝕變帶綜合驗(yàn)證結(jié)果Fig.8 Comprehensive verification results for alteration belt of gold mineralization at Nandatan-Xihongtan section

在南大灘-西紅灘地段東部的西紅灘地區(qū)(圖7)開展了地表撿塊地球化學(xué)樣品采集和室內(nèi)分析測試，測試結(jié)果(圖8d)顯示6個樣品的金元素含量高于0.3×10-6，最高值可達(dá)1.839×10-6，這還僅是地表巖石的金元素含量，深部巖石的金含量還會增高。多數(shù)樣品的二氧化硅含量高于80%甚至超過90%，遠(yuǎn)超花崗巖平均值(66%)，表明蝕變巖中存在較強(qiáng)的硅化現(xiàn)象?？侳eO含量也在2%上下浮動，為含鐵礦物黃鐵礦、褐鐵礦的顯示。此外，銀元素也存在高值異常。值得關(guān)注的是，南大灘-西紅灘地段的西側(cè)還存在相似的巖性-蝕變地段(圖7)，后續(xù)將開展相關(guān)礦化查證工作。

5 結(jié)論

(1)針對資源一號02D遙感數(shù)據(jù)特點(diǎn)開展了高光譜數(shù)據(jù)預(yù)處理、降維增值和SAM光譜匹配方法改進(jìn)，開發(fā)了基于python語言的資源一號02D高光譜數(shù)據(jù)批量基礎(chǔ)預(yù)處理、改進(jìn)的SAM匹配的軟件模塊。在甘肅頭吊泉-南大灘地區(qū)進(jìn)行了高光譜蝕變礦物填圖，成功提取了褐鐵礦、綠泥石、方解石、白云石、短波云母、中短波云母、中長波云母共7種蝕變礦物。

(2)資源一號02D高光譜數(shù)據(jù)在該地區(qū)提取的蝕變礦物分布具有明顯規(guī)律性。一是蝕變礦物與地質(zhì)體巖性關(guān)系密切，多數(shù)與斷裂構(gòu)造并無明顯直接聯(lián)系；二是蝕變礦物并不是均勻分布在地質(zhì)體中，而是集中在特定地質(zhì)體的局部地段。方解石和白云石碳酸鹽類蝕變礦物的分布區(qū)可以指示碳酸鹽類沉積/變質(zhì)地層的出露范圍，蝕變的定向展布揭示了地層的走向。

(3)根據(jù)資源一號02D高光譜數(shù)據(jù)提取的蝕變礦物，結(jié)合地質(zhì)背景信息分析，新發(fā)現(xiàn)1條長約15 km的金礦化蝕變帶。東部的西紅灘地段經(jīng)野外地質(zhì)調(diào)查、光譜測量、地化采樣測試等綜合驗(yàn)證存在地表金礦化異常，西部地段有待后續(xù)開展相關(guān)查證工作。資源一號02D高光譜遙感數(shù)據(jù)為地質(zhì)找礦提供了成功案例，具有較好的地質(zhì)深入應(yīng)用潛力。

致謝：感謝中國資源衛(wèi)星應(yīng)用中心提供的資源一號02D高光譜數(shù)據(jù)，感謝論文評審專家對本文提出的寶貴意見。