張昭，陳川，李云鵬

1)新疆大學(xué)地質(zhì)與礦業(yè)工程學(xué)院，烏魯木齊，830017；2)新疆中亞造山帶大陸動力學(xué)與成礦預(yù)測重點實驗室，烏魯木齊，830017；3)新疆地礦局信息中心，烏魯木齊，830099

內(nèi)容提要: 遙感技術(shù)廣泛應(yīng)用于地質(zhì)基礎(chǔ)調(diào)查、礦產(chǎn)資源勘探、環(huán)境評估和地質(zhì)災(zāi)害調(diào)查等方面。它已從多光譜發(fā)展到高光譜階段，Landsat-8是目前最具有代表性和最常用的多光譜數(shù)據(jù)，ASTER具有高的分辨率和多波段特征，資源一號02D(ZY1-02D)衛(wèi)星是我國2019年發(fā)射的高光譜業(yè)務(wù)衛(wèi)星。為了更好地了解多源遙感數(shù)據(jù)在巖礦識別中的作用，在新疆東天山卡拉麥里地區(qū)進(jìn)行了相關(guān)研究。結(jié)果表明：Landsat-8 OLI的PCA變換結(jié)果清晰識別了研究區(qū)不同的巖性和地層；使用Landsat-8 OLI、ASTER和ZY1-02D高光譜數(shù)據(jù)，分別采取不同的圖像端元提取方法，在進(jìn)行光譜分析的基礎(chǔ)上，利用光譜角填圖(SAM)即可得到研究區(qū)的主要礦物分類圖件。通過野外驗證，應(yīng)用GIS技術(shù)進(jìn)行集成和分析，修正相關(guān)圖件后，便得到了精準(zhǔn)的礦物分類綜合圖。研究表明：多源遙感數(shù)據(jù)的集成在巖礦識別方面效果良好、前景巨大。

近年來，遙感技術(shù)在地質(zhì)構(gòu)造分析、巖石礦物識別、蝕變異常提取、礦產(chǎn)資源勘探等方面應(yīng)用日趨普遍和成熟(Sabins，1999；張玉君等，2002；李志忠等，2009；王潤生等，2011；Peyghambari et al.，2021)。新疆東準(zhǔn)噶爾卡拉麥里地區(qū)干旱的氣候條件、戈壁灘型自然景觀、良好的巖性露頭給該地區(qū)利用遙感技術(shù)進(jìn)行巖石礦物識別和開展相關(guān)礦產(chǎn)資源勘探工作提供了有利條件。國內(nèi)外學(xué)者使用多光譜遙感數(shù)據(jù)在不同的地區(qū)開展了巖礦識別的大量研究工作(李培軍等，2007；Khan et al.，2008；黃照強等，2010；Rajendran et al.，2011；別小娟等，2013；Pournamdari et al.，2014；Abdelaziz et al.，2017；楊偉光等，2018；??rtük et al.，2020；Ahmadi et al.，2021)。李培軍等(2007)使用ASTER數(shù)據(jù)，采用光譜角填圖法(Spectral Angle Mapper，SAM)對青海德爾尼蛇綠巖的主要巖石組成和蝕變礦物進(jìn)行了探測，較有效地識別了蛇綠巖的主要巖性和相關(guān)礦物。楊偉光等(2018)在西藏羅布莎地區(qū)利用WorldView-2高空間分辨率遙感影像和ETM+數(shù)據(jù)對羅布莎鉻鐵礦床進(jìn)行遙感控礦因素提取，使用了一系列光譜增強方法，評估遙感數(shù)據(jù)對鉻鐵礦成礦巖體地幔橄欖巖識別與圈定的能力和可行性，證明利用遙感技術(shù)探測潛在鉻鐵礦床的有效性。黃照強等(2010)、別小娟等(2013)也在西藏羅布莎地區(qū)利用遙感數(shù)據(jù)，采用不同的影像處理技術(shù)對蛇綠巖套和鉻鐵礦成礦相關(guān)的巖礦信息進(jìn)行了提取和分析，并取得了不錯的效果。Rajendran(2012)等在阿曼北部塞梅爾蛇綠巖地塊使用Landsat TM和ASTER數(shù)據(jù)，選用去相關(guān)拉伸(Decorrelated Stretching，DS)、不同的波段比值(Band Rationing，BR)和主成分分析(Principal Component Analysis，PCA)等圖像處理技術(shù)，成功圈定該地區(qū)蛇綠巖帶內(nèi)潛在的鉻鐵礦礦化帶。Abdelaziz(2017)等在阿富汗洛加爾地塊利用Landsat-8數(shù)據(jù)采用波段組合、主成分分析、波段比值和監(jiān)督分類技術(shù)相結(jié)合的方法對洛加爾省鉻鐵礦有關(guān)的蛇紋巖進(jìn)行了分離。??rtük(2020)等基于ASTER和Landsat-8數(shù)據(jù)分析了土耳其安納托利亞中部布爾巴斯蛇綠巖巖石分布，采用彩色波段比值、主成分分析和最小噪聲分離變換(Minimum Noise Fraction，MNF)對遙感數(shù)據(jù)進(jìn)行處理，有效地描繪了研究區(qū)域中輝石、變質(zhì)巖和輝長巖的空間分布以及一批以前未被識別的蛇紋石化純橄欖巖、輝石、變質(zhì)巖底巖和片麻巖。

圖 1 準(zhǔn)噶爾盆地地質(zhì)簡圖(a)(據(jù)Chen et al.，2004)和晚古生代額爾齊斯—卡拉麥里島弧增生雜巖帶構(gòu)造單元劃分略圖(b)(據(jù)陜西省地質(zhì)調(diào)查院滴水幅(L45C003004)1∶250, 000萬區(qū)域地質(zhì)調(diào)查報告修編，2015年)Fig. 1 Geological sketch of Junggar Basin (a) (according to Chen et al., 2004) and outline of tectonic unit division of the Late Paleozoic Erqis karamaili island arc accretionary complex belt (b) (Modified according to the 1∶250, 000 regional geological survey report of Dishui map sheet (L45C003004) of Shaanxi Geological Survey Institute, 2015)

高光譜數(shù)據(jù)在巖礦識別和礦產(chǎn)勘查中的應(yīng)用也十分廣泛(Kruse，1996；Roy et al.，2009；Bishop et al.，2011；Govil et al.，2018；魏紅艷等，2020；Govil et al.，2021)。李根軍等(2020)應(yīng)用ZY1-02D高光譜數(shù)據(jù)提取了大理巖、二長花崗巖等巖性信息，方解石、白云石等造巖礦物信息以及綠泥石、褐鐵礦等蝕變礦物信息，與實測成果較為吻合，表明該數(shù)據(jù)對巖礦信息的識別效果較好。李娜等(2020)應(yīng)用ZY1-02D高光譜數(shù)據(jù)在哈密遙感地質(zhì)試驗場地區(qū)主要對其在巖性-構(gòu)造信息識別、礦物信息提取等方面開展了應(yīng)用評價，清晰地顯示和區(qū)分了不同地層、巖體、構(gòu)造等地質(zhì)體。連琛芹等(2020)應(yīng)用GF-5高光譜數(shù)據(jù)在植被覆蓋區(qū)進(jìn)行了蝕變信息提取研究，提取的蝕變異常與研究區(qū)地質(zhì)信息吻合較好。董新豐等(2020)，在甘肅柳園地區(qū)利用GF-5高光譜衛(wèi)星數(shù)據(jù)開展礦物填圖應(yīng)用，對其在礦產(chǎn)資源方面的應(yīng)用前景進(jìn)行了初步評價。

由以上文獻(xiàn)綜述可知，多光譜數(shù)據(jù)在巖礦識別，特別是對鉻鐵礦成礦有關(guān)的蛇綠巖的提取中應(yīng)用廣泛、效果顯著，而高光譜遙感技術(shù)在巖礦識別中也獲得了比較成功的應(yīng)用。但使用多光譜數(shù)據(jù)和高光譜數(shù)據(jù)相結(jié)合，以更好發(fā)揮其各自優(yōu)勢而實現(xiàn)巖礦識別的應(yīng)用卻不多。

筆者等以新疆東準(zhǔn)噶爾卡拉麥里地區(qū)為研究對象，利用Landsat-8和ASTER(Advanced Spaceborne Thermal Emission and Reflection Radiometer，高級星載熱輻射和反射探測器)多光譜衛(wèi)星遙感數(shù)據(jù)和資源一號02D(ZY1E衛(wèi)星)AHSI高光譜數(shù)據(jù)，采用連續(xù)最大角凸錐(SMACC)、最小噪聲分離(MNF)、純凈像元指數(shù)(PPI)、n維可視化(n-D Visualizer)、光譜角制圖法(SAM)等遙感圖像處理算法，對研究區(qū)的巖石和礦物進(jìn)行識別，并結(jié)合野外實地調(diào)查和分析，應(yīng)用GIS軟件進(jìn)一步對遙感圖像進(jìn)行巖礦信息提取和空間展布分析。從而對研究區(qū)的巖礦進(jìn)行分類，以為該地區(qū)下一步更大比例尺的地質(zhì)調(diào)查工作和鉻鐵礦的勘探工作提供參考依據(jù)。

Ⅰ—阿爾泰微陸塊；Ⅱ2—晚古生代額爾齊斯—卡拉麥里島弧增生雜巖帶: —二臺—托讓庫都克島(額爾齊斯)弧滑脫—推斷褶皺帶, —托爾特庫勒—庫蘭喀茲干周緣盆地滑脫—推覆帶, —卡姆斯特—奧塔烏克塔什(卡拉麥里)島弧褶斷帶, —卡拉麥里蛇綠混雜巖帶, —卡拉麥里被動陸緣沖斷—滑脫帶, —芒能奧依—塔爾蘇溝板內(nèi)伸展盆地褶皺帶；Ⅲ—準(zhǔn)噶爾微陸塊: Ⅲ1—東準(zhǔn)山前(間)斷(坳)陷帶Ⅰ—Altay microblock; Ⅱ2—Late Paleozoic Ertix—Karamaili island arc accretionary complex belt: intraplate extensional basin fold belt;Ⅲ—Junggar microblock: Ⅲ1—Eastern Junggar piedmont (inter) fault (depression) belt

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

新疆卡拉麥里蛇綠混雜巖帶位于準(zhǔn)噶爾盆地東北緣，是中國西北部東準(zhǔn)噶爾南緣一條重要的蛇綠混雜巖帶(圖1a)，它處于西伯利亞與哈薩克斯坦—準(zhǔn)噶爾兩大板塊的交匯部位，是中亞造山帶的一個重要構(gòu)造單元(李錦軼等，1990、2009；李現(xiàn)冰，2013；陳新蔚等，2013；樊立飛等，2014；黃崗等，2017)。卡拉麥里蛇綠混雜巖帶在構(gòu)造單元上屬于晚古生代島弧增生巖帶(Ⅱ2)，其受控于近NWW向的卡拉麥里區(qū)域性深大斷裂(圖1b)。

研究區(qū)主要地層為石炭系和泥盆系，僅在西南部出露少量的三疊系和志留系地層(圖2)。石炭系地層以下石炭統(tǒng)姜巴斯套組為代表，其次是塔木崗組和上石炭統(tǒng)巴塔瑪依內(nèi)山組。其中姜巴斯套組可以分為兩段(圖2)，第一段主要為一套陸源碎屑巖、火山碎屑巖夾少量火山熔巖組合，主要巖性為凝灰質(zhì)礫巖、凝灰質(zhì)含礫長石巖屑粗砂巖、凝灰質(zhì)粗砂巖、凝灰質(zhì)粉砂巖、巖屑凝灰?guī)r夾玄武質(zhì)火山角礫巖、玄武玢巖等；第二段為一套碎屑巖組合，巖性為淺灰色砂巖、細(xì)砂巖、長石巖屑砂巖、長石砂巖、凝灰質(zhì)礫巖、凝灰質(zhì)長石粗砂巖、巖屑粗砂巖。其在北部下與下石炭統(tǒng)黑山頭組，上與上石炭統(tǒng)巴塔瑪依內(nèi)山組為不整合接觸，南部與下石炭統(tǒng)塔木崗組為斷層接觸或與卡拉麥里蛇綠構(gòu)造混雜巖帶為不整合接觸，上部與那林卡拉組為角度不整合接觸。塔木崗組主要分布于卡拉麥里山南一帶，其呈北西—南東向、近東西向展布，主要為一套陸源碎屑巖、火山碎屑巖夾火山熔巖組合，其下與泥盆系卡拉麥里組為不整合或整合接觸，上部被上石炭統(tǒng)巴塔瑪依內(nèi)山組角度不整合覆蓋。研究區(qū)出露的塔木崗組第一段巖性主要為暗灰色蝕變玄武玢巖、綠灰色流紋質(zhì)英安斑巖、凝灰質(zhì)細(xì)砂巖、安山玢巖、流紋巖、英安斑巖、灰紫色流紋斑巖、酸性角礫熔巖、紫灰色中酸性熔結(jié)凝灰?guī)r、灰綠色細(xì)粒安山質(zhì)凝灰角礫巖、灰黃色中酸性凝灰?guī)r等。巴塔瑪依內(nèi)山組主要為火山熔巖夾火山碎屑巖組合，其角度不整合于下石炭統(tǒng)塔木崗組、姜巴斯套組等之上，或與下石炭統(tǒng)那林卡拉組為斷層接觸。泥盆系的地層以北塔山組分布最為廣泛，其總體為一套火山碎屑巖、火山熔巖以及陸源碎屑巖組合，在研究區(qū)西部主要為火山碎屑巖、陸源碎屑巖夾少量熔巖組合，主要分布于卡拉麥里蛇綠構(gòu)造混雜巖帶及南部。

圖 2 卡拉麥里地區(qū)地質(zhì)礦產(chǎn)圖及采樣位置(據(jù)陜西省地質(zhì)礦產(chǎn)局區(qū)域地質(zhì)礦產(chǎn)研究院2009～2012年北塔山牧場幅(L46C003001)、滴水泉幅(L45C003004)1∶25萬地質(zhì)圖修編)Fig. 2 Geological and mineral map and sampling location of Karamaili area (revised according to the 1∶250000 geological map of Beitashan ranch (L46C03001) and Dishuiquan (L45C03004) of Regional Institute of Geology and Mineral Resources of Shaanxi Bureau of Geology and Mineral Resources from 2009 to 2012)T1js1—尖山溝組第一段；C2bt1—巴塔瑪依內(nèi)山組第一段；C2h—弧形梁組；C1n2—那林卡拉組第二段；C1t1—塔木崗組第一段；C1j2—姜巴斯套組第二段；C1j1—姜巴斯套組第一段；D2bt2—北塔山組第二段；D2bt1—北塔山組第一段；Dk—卡拉麥里組；S3-D1h—紅柳溝組；S2b—白山包組；C2ηγβb—中細(xì)粒似斑狀黑云二長花崗巖；C2ηγβa—中粗粒似斑狀黑云二長花崗巖；C2χργψc—中粒角閃石堿長花崗巖；C2χργb—中細(xì)粒鈉鐵閃石堿長花崗巖；C2χργa—中粒鈉鐵閃石堿長花崗巖；D3γδο—灰白色中粗粒英云閃長巖；C1δο—灰色細(xì)粒石英閃長巖；D3γδο—中細(xì)粒英云閃長巖；DCsi—硅質(zhì)巖、硅質(zhì)粉砂巖；DCβ—塊狀與枕狀玄武巖；DCν—堆晶輝長巖；DCσ—變質(zhì)橄欖巖、輝石巖、蛇紋巖等；δμ—閃長玢巖脈；ν—輝長巖脈T1js1—The 1st Mem. of Jianshangou Fm.; C2bt1—The 1st Mem. of Batamayineishan Fm.; C2h—Huxingliang Fm.; C1n2—The 2st Mem. of Nalinkala Fm.; C1t1—The 1st Mem. of Tamugang Fm.; C1j2—The 2st Mem. of Jiangbasitau Fm.; C1j1—The 1st Mem. of Jiangbasitau Fm.; D2bt2—The 2st Mem. of Beitashan Fm.; D2bt1—The 1st Mem. of Beitashan Fm.; Dk—Karamaili Fm.; S3—D1h—Hongliugou Fm.; S2b—Baishanbao Fm.Karamaili Fm.; C2ηγβb—Medium fine grained porphyritic biotite monzogranite；C2ηγβa—Medium coarse grained porphyritic biotite monzogranite；C2χργψc—Medium grained hornblende alkali feldspar granite；C2χργb—Medium fine grained amphibole alkali feldspar granite；C2χργa—Medium grained amphibole alkali feldspar granite；D3γδο—Gray white medium coarse grained tonalite；C1δο—Gray fine-grained quartz diorite；D3γδο—Medium fine grained tonalite；DCsi—Siliceous rock, siliceous siltstone；DCβ—Massive and Pillow Basalt；DCν—Cumulate gabbro；DCσ—Metamorphic peridotite, pyroxenite, serpentinite, etc；δμ—Diorite porphyrite dyke；ν—Gabbro dyke

研究區(qū)的侵入巖主要是晚石炭世黑云二長花崗巖和鈉鐵閃石堿長花崗巖，早石炭世石英閃長巖和英云閃長巖。

卡拉麥里蛇綠混雜巖帶位于研究區(qū)南部的卡拉麥里山一帶，且主要分布于卡拉麥里大斷裂的北側(cè)，其主要由晚古生代蛇綠巖殘片、蛇綠巖上覆巖系、外來巖片等構(gòu)造塊體和變形基質(zhì)兩部分組成(韓秉峻，2019)。蛇綠巖分布區(qū)內(nèi)出露的地層主要是泥盆系和石炭系，巖石組合以陸源碎屑巖、火山碎屑巖和熔巖為特征(趙磊等，2019)，其代表性的地層單元是中泥盆統(tǒng)北塔山組和下石炭統(tǒng)姜巴斯套組，并且這這些地層大部分逆沖到蛇綠巖之上，而局部地段可見下石炭統(tǒng)姜巴斯陶組與下伏變質(zhì)超基性巖呈不整合接觸關(guān)系。蛇綠混雜巖帶內(nèi)各肢解蛇綠巖巖塊多以斷裂為界混雜堆疊，受剪切作用明顯，巖塊原生構(gòu)造被后期構(gòu)造改造或置換，但蛇綠巖的各組成單元出露仍較為齊全，基本上包括了典型蛇綠巖所有巖石類型。在卡拉麥里山東段紅柳溝地區(qū)，組分發(fā)育不全，主要為不同程度蝕變的方輝橄欖巖、二輝橄欖巖、輝長巖、玄武巖等。中段清水一帶組分發(fā)育相對較全，包括不同蝕變程度的方輝橄欖巖、含鉻鐵礦的純橄巖、二輝橄欖巖、單斜輝石巖、堆晶輝長巖、角閃巖、輝綠巖、斜長巖、斜長花崗巖、枕狀玄武巖、塊狀玄武巖、放射蟲硅質(zhì)巖等。西段由于受露頭的影響，可見組分有蛇紋巖、蛇紋石化橄欖巖、玄武巖、輝長巖、斜長花崗巖和硅質(zhì)巖等?；|(zhì)部分主要為泥盆—石炭紀(jì)的變凝灰泥質(zhì)長石粉砂巖、變晶屑巖屑凝灰?guī)r、凝灰?guī)r、千枚巖、糜棱巖化凝灰?guī)r、糜棱巖等。研究區(qū)的礦產(chǎn)以鉻鐵礦為主，此外還有少量的金礦、赤鐵礦—褐鐵礦分布(張峰等，2014)(圖2)。鉻鐵礦礦床和礦點產(chǎn)于蛇綠巖上部的斜輝橄欖巖、純橄巖中，但多集中在純橄巖，這些巖石均已蝕變?yōu)樯呒y巖或蛇紋石化橄欖巖。

2 遙感數(shù)據(jù)及巖礦識別方法

2.1 遙感數(shù)據(jù)

2.1.1 數(shù)據(jù)簡介

本文根據(jù)研究目的和制圖需求，選用Landsat-8、ASTER多光譜遙感數(shù)據(jù)以及資源一號02D(ZY1E衛(wèi)星)AHSI高光譜數(shù)據(jù)，現(xiàn)對不同的衛(wèi)星傳感器性能特征介紹如下。

Landsat-8 OLI數(shù)據(jù)包含9個波段，其中Band 1～7和Band 9的空間分辨率為30 m，Band 8全色波段空間分辨率為15 m，成像幅寬為185 km × 185 km。本文使用了一景Landsat-8數(shù)據(jù)，其成像時間為2021年4月14日，其條帶號和行編號分別為141/29，云量覆蓋為0.28。ASTER數(shù)據(jù)含14個波段信息，分別是空間分辨率為15 m的3個可見光近紅外(VNIR)波段，其光譜范圍在0.52～0.86 μm之間，空間分辨率為30 m的6個短波紅外波段(SWIR)，其光譜范圍在1.6～2.43 μm之間，空間分辨率為90 m的8個熱紅外波段，其光譜范圍在在8.125～11.65μm之間。其成像幅寬為60 km × 60 km。本文使用了4景ASTER數(shù)據(jù)(表1)。資源一號02D AHSI高光譜數(shù)據(jù)，空間分辨率為30 m，波段數(shù)量為166個，光譜范圍為0.4～2.5 μm，成像幅寬為60 km × 60 km。光譜分辨率在可見光—近紅外為10 nm、在短波紅外為20 nm。本次使用的數(shù)據(jù)景號是ZY1E_AHSI_E89.90_N45.01_20201027_005899_L1A0000192836，成像日期為2020年10月27日，云量覆蓋為3.6%。

表1 卡拉麥里地區(qū)覆蓋的4景ASTER數(shù)據(jù)基本信息Table 1 Basic information of ASTER data of 4 scenes covered in the Karamaili area

2.1.2數(shù)據(jù)預(yù)處理

本次獲取的Landsat-8數(shù)據(jù)為L1TP級，對Landsat-8 OLI數(shù)據(jù)進(jìn)行輻射定標(biāo)、FLAASH大氣校正(主要參數(shù)見表2)。

表2 不同遙感數(shù)據(jù)FLAASH大氣校正主要參數(shù)Table 2 Main parameters of FLAASH atmospheric correction for different remote sensing

由于ASTER的VNIR和SWIR分別處于不同的傳感器子系統(tǒng)，故通過“Layer Stacking(波段疊加)”將每景的Band 1～9放在同一多光譜數(shù)據(jù)集中，并將SWIR的6個空間分辨率為30 m的波段重采樣至15 m，然后依次進(jìn)行輻射定標(biāo)、FLAASH大氣校正(主要參數(shù)見表2)。對高光譜數(shù)據(jù)ZY1-02D AHSI首先進(jìn)行輻射定標(biāo)，然后將多個波段拆分成單波段TIFF格式以便后續(xù)處理。ZY1-02D高光譜短波紅外范圍內(nèi)數(shù)據(jù)條紋現(xiàn)象明顯(圖3)，故對圖像的噪聲特征進(jìn)行分析后，使用傅立葉變換頻譜分析和反傅立葉變換的方法，采用Python編程的方法去除影像中的條紋(代碼見附錄)。將去除條紋的數(shù)據(jù)通過波段疊加(Band Stack)重新組合成166個波譜數(shù)據(jù)，然后選用FLAASH大氣校正模型對其進(jìn)行校正(主要參數(shù)見表2)。以Landsat-8的全色波段為參考影像對ZY1-02D高光譜進(jìn)行幾何校正。

圖 3 ZY1-02D AHSI去除條紋前后影像對比(選自短波紅外第133波段)Fig. 3 Image comparison before and after fringe removal of ZY1-02D AHSI (selected from the 133rd band of short wave infrared)

2.2 巖礦識別技術(shù)方法

利用遙感數(shù)據(jù)識別地層巖性和礦物信息時，由于不同傳感器光譜及空間分辨率的差異，致使對目標(biāo)地物的識別存在明顯的專屬性，因此筆者等選用Landsat-8 OLI、ASTER和ZY1-02D AHSI 3種遙感數(shù)據(jù)，配合開展地層巖性和礦物識別，以期達(dá)到對研究區(qū)的巖石地層和相關(guān)礦物進(jìn)行識別分類的目的(圖4)。

圖 4 多源遙感數(shù)據(jù)地層巖性和巖石礦物識別流程Fig. 4 Identification process of formation lithology, rock and mineral from multi-source remote sensing data

2.2.1數(shù)據(jù)降維

為有效識別出遙感影像海量信息中的有效數(shù)據(jù)，需首先對遙感影像進(jìn)行降維(Dimension reduction)處理，即將高維數(shù)據(jù)轉(zhuǎn)換為低維數(shù)據(jù)而不丟失信息的過程稱為降維(Dimension reduction)。目前，最常用的數(shù)據(jù)降維方法包括主成分分析法和最小噪聲分離法。

主成分分析(Principal Component Analysis，PCA)是一種圖像增強技術(shù)，它將高度相關(guān)的原始圖像數(shù)據(jù)轉(zhuǎn)換為一組主成分或PC波段的不相關(guān)變量，通過抑制在所有波段占主導(dǎo)地位的輻照度效應(yīng)來分離噪聲成分，從而達(dá)到降低數(shù)據(jù)維度的目的。

本次研究對Landsat-8 OLI數(shù)據(jù)進(jìn)行PCA變換之后，顯示前3個主成分包含了所有波段中99.54 %的信息(表3)，表示了數(shù)據(jù)集中的域變化，并突出顯示了所有波段的最常見特征。因此對前3個主成分PC3、PC2、PC1進(jìn)行假彩色合成，并結(jié)合地質(zhì)圖和野外勘查便得到了研究區(qū)的地層巖性信息(表4)。

表3 Landsat-8 OLI主成分分析結(jié)果統(tǒng)計表Table 3 Statistics of Landsat-8 OLI principal component analysis results

表4 Landsat-8 OLI主成分分析PC3、PC2、PC1假彩色合成解譯卡拉麥里地區(qū)主要地層巖性特征表Table 4 Main stratigraphic lithologic characteristics of landsat-8 oli principal component analysis PC3, PC2 and PC1 false color synthetic interpretation for Karamaili area

最小噪聲分離(Minimum Noise Fraction, MNF)是將一幅多波段圖像的主要信息集中在前面幾個波段中，主要作用是判斷圖像數(shù)據(jù)維數(shù)、分離數(shù)據(jù)中的噪聲,減少后處理中的計算量。它由Green等(1988)開發(fā)，該方法是由兩個獨立的主成分分析旋轉(zhuǎn)組成的線性變換。第一次旋轉(zhuǎn)利用噪聲協(xié)方差矩陣的主分量對數(shù)據(jù)中的噪聲進(jìn)行去相關(guān)和重縮放，使得有噪聲的數(shù)據(jù)具有單位方差，并且沒有波段間相關(guān)性。在第二階段，通過分析特征值和相關(guān)圖像來計算固有維數(shù)。因此，與以噪聲為主的圖像相關(guān)的單位特征值可以分離出來，從而增強光譜處理結(jié)果。此次研究對ASTER的9個波段和ZY1-02D AHSI經(jīng)過篩選的152個波段進(jìn)行MNF變換，以達(dá)到降維、抑制噪聲和信息增強的目的。

2.2.2圖像端元提取

混合像元廣泛存在于圖像中，影響遙感圖像的分類精度以及目標(biāo)探測效果(藍(lán)金輝等，2018)。若一個像素對應(yīng)的地面區(qū)域內(nèi)只包含一種特征地物，則稱此像素為純像元，而端元即為數(shù)據(jù)中代表類別特征的理想化純像元，提取遙感數(shù)據(jù)中的純像元技術(shù)稱為光譜端元提取(齊濱，2012；陳晉等，2016)。常用方法有基于連續(xù)最大角凸錐技術(shù)、純凈像元指數(shù)技術(shù)和n維可視化技術(shù)(王茂芝等，2015)。

基于連續(xù)最大角凸錐(sequential maximum angle convex cone，SMACC)使用凸錐模型表示矢量數(shù)據(jù)，直接從多光譜和高光譜數(shù)據(jù)中提取圖像端元。它通過識別多維空間中的極值點，并形成一個包含第一個端元的凸錐。然后，將約束斜投影應(yīng)用于現(xiàn)有圓錐體，圓錐體將增加以包括下一個端元。簡單來說，SMACC首先查找最亮的像素，然后查找與最亮像素差異最大的像素，然后查找與前兩個像素差異最大的像素。該過程將繼續(xù)進(jìn)行，直到投影衍生出一個已存在于凸錐內(nèi)的端元或達(dá)到指定數(shù)量的端元。在這項研究中，對經(jīng)過大氣校正后的Landsat-8 OLI數(shù)據(jù)使用SMACC方法進(jìn)行圖像端元自動提取，端元波譜提取數(shù)量設(shè)置為10，誤差容限值保持默認(rèn)值0，分離端元波譜的約束條件選擇“Posivity only”，即把每個波長的正值端元波譜作為約束條件。在合并相似端元波譜后，共提取出6個圖像端元的波譜。

純凈像元指數(shù)( Pixel Purity Index，PPI)和n維可視化(n- Dimensionsl Visualizer )用于收集ASTER和ZY1-02D AHSI圖像中的端元波譜。通過計算多波譜和高光譜圖像的純凈像元指數(shù)(PPI)，進(jìn)而在圖像中尋找波譜最“純”的像元。這可以通過連續(xù)收集n維散點圖中的極值來確定(Kruse，1998；Govil et al.，2021)。n-D Visualizer用于定位、識別、聚集數(shù)據(jù)集中最純的像元，從而獲取純凈的端元波譜，以生成光譜庫文件，進(jìn)而對地物進(jìn)行分類。因為經(jīng)過MNF變換后前幾個波段將包含90%以上的信息，所以選取ASTER的前三個波段，ZY1-02D AHSI的前十個波段進(jìn)行PPI的提取。在執(zhí)行PPI時，ASTER和ZY1-02D AHSI的迭代次數(shù)均設(shè)置為10000，前者的閾值系數(shù)設(shè)置為2.5，后者為3.0。針對ASTER和ZY1-02D AHSI數(shù)據(jù)分別提取了7個和9個圖像端元的波譜。

2.2.3光譜分析

光譜分析(Spectral Analyst)根據(jù)礦物的光譜特征來進(jìn)行礦物識別。此次研究使用ENVI軟件中的光譜分析算法，如二進(jìn)制編碼(Binary Encoding)、光譜角填圖(Spectral Angle Mapper)和光譜特征擬合(Spectral Feature Fitting)，對提取的圖像端元光譜與USGS波譜庫典型礦物波譜進(jìn)行匹配，從而識別出端元光譜的目標(biāo)礦物種類。將每種匹配算法的權(quán)重設(shè)置為1，則總得分為3，總得分越靠近3，說明圖像端元光譜與典型礦物波譜匹配度越高(表5)。研究區(qū)幾種典型礦物圖像端元波譜和USGS標(biāo)準(zhǔn)礦物波譜庫波譜對比見圖5。

表5 卡拉麥里地區(qū)不同遙感數(shù)據(jù)礦物識別時光譜分析匹配得分表Table 5 Matching scores of spectral analysis for mineral identification of different remote sensing data in Karamaili area

圖 5 卡拉麥里地區(qū)典型礦物圖像端元波譜(紅色)和標(biāo)準(zhǔn)波譜庫波譜(藍(lán)色)對比圖。 OLI數(shù)據(jù)：(a)橄欖石、(b)黑云母、(c)綠泥石；ASTER數(shù)據(jù)(d)蛇紋石、(e)輝石、(f)綠泥石；AHSI數(shù)據(jù)：(g)橄欖石、(h)高嶺石、(i)綠泥石Fig. 5 Comparison diagram of typical mineral image endmemberspectral (red) and standard spectral library (blue) in Karamaili area. OLI data: (a) olivine, (b) biotite, (c) chlorite; ASTER data: (d) serpentine, (E) pyroxene, (f) chlorite; AHSI data: (g) olivine, (H) kaolinite, (I) chlorite

2.2.4波譜角填圖

波譜角填圖(Spectral Angle Mapper，SAM)本質(zhì)上是一種監(jiān)督分類技術(shù)，它將圖像光譜與參考光譜相匹配，在多維特征空間中確定圖像光譜和參考光譜之間的角度。角度越小，說明圖像光譜和參考光譜之間的相似性越高。光譜角的值介于0到π/2之間，根據(jù)Kruse等人(1993)給出的公式計算：

(1)

其中θ是圖像光譜和參考光譜之間的角度，n是光譜波段數(shù)；t是實際光譜(圖像光譜)的反射率，r是參考光譜的反射率。

本次研究對Landsat-8 OLI經(jīng)過SMACC提取端元波譜和ASTER、ZY1-02D AHSI經(jīng)過MNF、PPI和n維可視化提取端元波譜后，分別對三種數(shù)據(jù)進(jìn)行光譜分析匹配，然后采用SAM進(jìn)行研究區(qū)礦物的識別制圖(圖6)。

圖 6 (a)Landsat-8 OLI礦物識別分類圖；(b)ASTER礦物識別分類圖；(c)ZY1-02D AHSI礦物識別分類圖Fig. 6 (a) Landsat-8 OLI mineral identification and classification diagram; (b) ASTER mineral identification classification map; (c) ZY1-02D AHSI mineral identification classification diagram

3 結(jié)果制圖、野外查驗及分析討論

3.1 OLI巖性地層制圖和分析討論

經(jīng)過PCA(主成分分析，principal component analysis)變換后，圖像色彩豐富、可獲取信息量大幅度提升。根據(jù)不同地層巖性在主成分分析后的假彩色合成影像(R：PC3，G：PC2，B：PC1)所表現(xiàn)出來的顏色、紋理、形狀、水系等特征，參考地質(zhì)資料便可識別出研究區(qū)的全部地層和絕大多數(shù)巖性(表5)。但是由于蛇綠巖混雜巖構(gòu)成的復(fù)雜性、OLI數(shù)據(jù)光譜分辨率和空間分辨率的限制，并不能識別出堆晶輝長巖、塊狀玄武巖等小型的巖漿巖和硅質(zhì)巖等沉積巖。對巖性組成非常相似的地層，如北塔山組第二段和北塔山組第一段中的砂巖和粉砂巖，需要經(jīng)過野外勘查才能正確劃分其組分。

3.2 Landsat-8 OLI、ASTER和ZY1-02D AHSI礦物識別制圖和分析討論

本研究針對Landsat-8 OLI大氣校正之后的圖像，使用SMACC方法提取圖像端元波譜，針對ASTER和ZY1-02D AHSI大氣校正、條紋去除、MNF處理之后的圖像，采用PPI和n維可視化工具提取端元波譜，并以USGS標(biāo)準(zhǔn)礦物波譜庫的波譜為參考，進(jìn)行光譜分析以確定提取的端元波譜所代表的礦物種類(表5)，最后把端元波譜作為已知波譜，采用SAM方法和圖像波譜進(jìn)行比對分析，以達(dá)到礦物分類識別(圖6)的目的。

由礦物識別分類圖可知，不管是多光譜數(shù)據(jù)Landsat-8 OLI和ASTER，還是高光譜數(shù)據(jù)ZY1-02D AHSI，對于一些造巖礦物，如長石、石英、黑云母的識別特別有效，但是由于這些礦物往往具有相似的光譜特征，因此本研究只能局限于將這些礦物的組合識別出來，如黑云母、石英、長石組合，鉀長石、石英組合。此外由于蛇綠巖混雜巖構(gòu)造的復(fù)雜性、“異物同譜”與“同物異譜”現(xiàn)象的存在以及遙感本身獲取的就是混合像元，因此并不能單獨識別某一種礦物，僅能識別出橄欖石和輝石、蛇紋石和斜方輝石的礦物組合及其它一些礦物組合。Landsat-8 OLI和ZY1-02D AHSI對綠泥石的識別效果好，ASTER對橄欖石和輝石組合、蛇紋石和高嶺石識別效果佳。而ZY1-02D AHSI由于高的光譜分辨率，不管是礦物的識別種類還是效果都優(yōu)于前兩種多光譜數(shù)據(jù)。但是它卻未能將方解石、白云石等Landsat-8 OLI和ASTER可識別的礦物鑒定出來。因此將多源遙感數(shù)據(jù)進(jìn)行集成更有利于研究區(qū)礦物的識別，不同的遙感數(shù)據(jù)源可起到互相印證和互補改正的作用。

3種不同數(shù)據(jù)識別出的礦物分布范圍差異可能

是由于遙感數(shù)據(jù)本身的特性引起的，如Landsat-8 OLI和ZY1-02D高光譜數(shù)據(jù)雖然空間分辨率都是30m，但是后者的光譜分辨率遠(yuǎn)遠(yuǎn)高于前者。這就意味后者對礦物之間微小差距的探測更有效。而前者由于光譜分辨率低，對一些波譜特征在局部變化不大的礦物識別更有效。當(dāng)然還有一些光照差異、地形起伏等因素的影響也會造成不同數(shù)據(jù)識別礦物性能的差異。

3.3 野外采樣及驗證

在對研究區(qū)的地層巖性分布狀況進(jìn)行深入了解的基礎(chǔ)上，結(jié)合Landsat-8 OLI數(shù)據(jù)PCA方法識別出的巖性地層特征，以及Landsat-8 OLI、ASTER和ZY1-02D AHSI數(shù)據(jù)SAM礦物識別分類圖(圖6)，進(jìn)行了野外勘查和采樣工作。根據(jù)野外實地調(diào)查，Landsat-8 OLI巖性地層分類結(jié)果和3種數(shù)據(jù)的礦物分類結(jié)果都得到了很好的驗證，但是由于“異物同譜”和“同譜異物”現(xiàn)象的存在以及遙感器性能的差異，分類時仍然將一些巖性礦物信息錯歸或遺漏。本次研究共采集了80個巖石樣品(圖2)，并將采樣點作為實際驗證點，結(jié)果發(fā)現(xiàn)69個巖石樣品的野外鑒定命名和不同遙感數(shù)據(jù)礦物識別分類結(jié)果相吻合，其準(zhǔn)確度達(dá)到了86.25%。研究區(qū)典型的巖石標(biāo)本照片見圖7。

圖7 野外勘查采集的巖石標(biāo)本：(a)橄欖巖；(b)含黏土礦物的礫巖；(c)石英巖;(d)蛇紋巖；Fig. 7 Rock samples collected from field exploration: (a) peridotite; (b) conglomerate containing clay minerals; (c) quartzite; (d) serpentinite;

3.4 礦物識別分類結(jié)果集成與分析

經(jīng)過3.2不同遙感數(shù)據(jù)礦物識別結(jié)果的討論和分析可知，Landsat-8 OLI、ASTER和ZY1-02D AHSI 3種數(shù)據(jù)對研究區(qū)不同的礦物識別都起著至關(guān)重要的作用，因此考慮將多源遙感數(shù)據(jù)分類結(jié)果進(jìn)行集成與分析，以達(dá)到最優(yōu)的礦物識別分類效果。本研究中，應(yīng)用GIS軟件對3種數(shù)據(jù)的礦物分類結(jié)果以及野外實地勘查驗證資料進(jìn)行綜合分析和數(shù)據(jù)集成。首先將遙感柵格分類圖轉(zhuǎn)換為shape矢量格式以便于進(jìn)行空間分析運算。隨后以地層巖性資料為基礎(chǔ)、以野外實地勘查采樣結(jié)果為標(biāo)準(zhǔn)，對多源遙感數(shù)據(jù)礦物分類圖進(jìn)行空間求交和合并運算以篩選出正確的礦物分布信息，對礦物分類圖件進(jìn)行調(diào)整(圖8)。

把多源遙感數(shù)據(jù)巖礦識別分類結(jié)果和野外地質(zhì)資料進(jìn)行集成分析，可彌補一些漏掉的礦物信息，并糾正一些錯誤的礦物分類結(jié)果。主要是對Landsat-8 OLI識別的橄欖石的分布范圍進(jìn)行了調(diào)整，對將一些黏土礦物而被錯誤識別為白云石的范圍進(jìn)行了刪減。對ASTER識別的蛇紋石和斜長石范圍進(jìn)行調(diào)整，經(jīng)野外驗證，ASTER提取的蛇紋巖信息多為蛇紋石化的橄欖巖，而一少部分為橄欖巖。對ZY1-02D AHSI識別的高嶺石、伊利石、蒙脫石等黏土礦物范圍進(jìn)行了調(diào)整。經(jīng)過修改調(diào)整后的礦物分類綜合圖清晰地顯示了研究區(qū)主要礦物的空間展布范圍，特別是對超基性巖區(qū)域礦物的構(gòu)成有了清晰的認(rèn)知，這對該區(qū)域與超基性巖有關(guān)的鉻鐵礦的勘探工作提供了重要的礦物分類圖件。

4 結(jié)論

筆者等將多光譜遙感數(shù)據(jù)Landsat-8和ASTER，以及ZY1-02D高光譜數(shù)據(jù)進(jìn)行集成，在新疆東天山卡拉麥里地區(qū)，針對不同數(shù)據(jù)進(jìn)行相應(yīng)預(yù)處理后，對Landsat-8 OLI采用PCA方法進(jìn)行地層巖性識別；對Landsat-8 OLI、ASTER和ZY1-02D AHSI數(shù)據(jù)使用MNF方法進(jìn)行數(shù)據(jù)降維和分離噪聲后，相應(yīng)采用SMACC、PPI結(jié)合n維可視化兩種提取圖像端元波譜方法，在光譜分析后，使用SAM方法進(jìn)行研究區(qū)礦物識別。并在野外地質(zhì)勘查的基礎(chǔ)上，對多源遙感數(shù)據(jù)礦物分類結(jié)果進(jìn)行集成分析而得到研究區(qū)精準(zhǔn)的礦物綜合分類圖。結(jié)論如下：

(1)在利用遙感技術(shù)進(jìn)行巖礦信息提取時，特別是礦物識別時，應(yīng)充分發(fā)揮多源遙感數(shù)據(jù)之間互相印證和補充細(xì)化的作用，以達(dá)到更精準(zhǔn)的礦物分類識別結(jié)果。

(2)在對卡拉麥里地區(qū)巖礦信息提取的過程中，經(jīng)野外實地勘查，不同的數(shù)據(jù)源都達(dá)到了較好的效果，3種數(shù)據(jù)對一些造巖礦物的組合識別效果都很好，但是每種數(shù)據(jù)都有其自己的獨特的優(yōu)勢和作用，Landsat-8識別出的白云石、ASTER識別出的橄欖石和輝石組合，ZY1-02D識別出的蛇紋石、斜方輝石組合都是其它兩種數(shù)據(jù)所未識別出來的。

(3)圖像端元波譜的提取在礦物識別中起著至關(guān)重要的作用，針對不同的數(shù)據(jù)源分別采用SMACC、PPI結(jié)合n維可視化來提取端元波譜得到了較為精準(zhǔn)的礦物識別分類圖件。

(4)巖礦識別工作應(yīng)充分發(fā)揮“3S”技術(shù)的優(yōu)勢，以遙感(RS)為數(shù)據(jù)源，在綜合地質(zhì)資料和遙感成果的基礎(chǔ)上，發(fā)揮全球定位系統(tǒng)(GPS)的優(yōu)勢，在野外進(jìn)行實地勘查驗證，在地理信息系統(tǒng)(GIS)軟件中，對不同數(shù)據(jù)源提取的巖礦信息、地質(zhì)綜合資料、野外勘查驗證結(jié)果進(jìn)行集成和綜合，從而得到精準(zhǔn)的巖礦分類識別結(jié)果。

附錄：ZY1-02D高光譜數(shù)據(jù)條紋去除Python代碼

import cv2

from osgeo import gdal

from matplotlib import pyplot as plt

import numpy as np

import os

def tif_read(tifpath, bandnum):

"""

Use GDAL to read data and transform them into arrays.

:param tifpath:tif文件的路徑

:param bandnum:需要讀取的波段

:return:該波段的數(shù)據(jù)，narray格式。len(narray)是行數(shù)，len(narray[0])列數(shù)

"""

image = gdal.Open(tifpath) #打開該圖像

if image == None:

print(tifpath + "該tif不能打開!")

return

lie = image.RasterXSize #柵格矩陣的列數(shù)

hang = image.RasterYSize #柵格矩陣的行數(shù)

im_bands = image.RasterCount #波段數(shù)

im_proj = image.GetProjection() #獲取投影信息

im_geotrans = image.GetGeoTransform() #仿射矩陣

#print(′該tif:{}個行，{}個列，{}層波段, 取出第{}層.′.format(hang, lie, im_bands, bandnum))

band = image.GetRasterBand(bandnum) # Get the information of band num.

band_array = band.ReadAsArray(0, 0, lie, hang) # Getting data from zeroth rows and 0 columns

# band_df = pd.DataFrame(band_array)

del image #減少冗余

return band_array, im_proj, im_geotrans

#def tif_write(filename, im_data, im_proj, im_geotrans):

def tif_write(filename, im_data):

"""

gdal數(shù)據(jù)類型包括

gdal.GDT_Byte,

gdal.GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,

gdal.GDT_Float32, gdal.GDT_Float64

:param filename:存出文件名

:param im_data:輸入數(shù)據(jù)

:param im_proj:投影信息

:param im_geotrans:放射變換信息

:return: 0

"""

if ′int8′ in im_data.dtype.name: #判斷柵格數(shù)據(jù)的數(shù)據(jù)類型

datatype = gdal.GDT_Byte

elif ′int16′ in im_data.dtype.name:

datatype = gdal.GDT_UInt16

else:

datatype = gdal.GDT_Float32

#判讀數(shù)組維數(shù)

if len(im_data.shape) == 3:

im_bands, im_height, im_width = im_data.shape

else:

im_bands, (im_height, im_width) = 1, im_data.shape #多維或1.2維

#創(chuàng)建文件

driver = gdal.GetDriverByName("GTiff")

dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

# dataset.SetGeoTransform(im_geotrans)

# dataset.SetProjection(im_proj)

if im_bands == 1:

dataset.GetRasterBand(1).WriteArray(im_data)

else:

for i in range(im_bands):

dataset.GetRasterBand(i + 1).WriteArray(im_data[i])

del dataset

def mquzao(img,filename):

#img = (img ～ np.amin(img)) / (np.amax(img) ～ np.amin(img)) * 255

#img_back_int = img.astype(np.uint8)

f = np.fft.fft2(img)

fshift = np.fft.fftshift(f)

mask = np.ones(img.shape, dtype=np.uint8)

mask[1015:1034, 0: 970] = 0

mask[1015:1034, 1031:] = 0

fshift = fshift * mask #濾波

f_ishift = np.fft.ifftshift(fshift)

img_back = np.fft.ifft2(f_ishift)

img_back2 = np.abs(img_back)

# img_back3 = (img_back2 ～ np.amin(img_back2)) / (np.amax(img_back2) ～ np.amin(img_back2)) * 255

# img_back_int = img_back3.astype(np.uint8)

pathout = ′OutPutImg\′

isExists = os.path.exists(pathout)

#判斷結(jié)果

if not isExists:

os.makedirs(pathout)

filename1 = pathout+filename

tif_write(filename1, img_back2)

return img_back2

inputpath = r′D:othersgaoguangtuqutiaowenInPutImg′

tiflist = os.listdir(inputpath)

for temptif in tiflist:

tifpath = inputpath+′\′+ temptif

img1,_,_ = tif_read(tifpath, bandnum=1)

mquzao(img1,filename=temptif)

print(′done′)