多源遥感数据在新疆卡拉麦里地区岩矿识别中的应用

张昭，陈川，李云鹏

1)新疆大学地质与矿业工程学院，乌鲁木齐，830017；
2)新疆中亚造山带大陆动力学与成矿预测重点实验室，乌鲁木齐，830017；
3)新疆地矿局信息中心，乌鲁木齐，830099

内容提要: 遥感技术广泛应用于地质基础调查、矿产资源勘探、环境评估和地质灾害调查等方面。它已从多光谱发展到高光谱阶段，Landsat-8是目前最具有代表性和最常用的多光谱数据，ASTER具有高的分辨率和多波段特征，资源一号02D(ZY1-02D)卫星是我国2019年发射的高光谱业务卫星。为了更好地了解多源遥感数据在岩矿识别中的作用，在新疆东天山卡拉麦里地区进行了相关研究。结果表明：Landsat-8 OLI的PCA变换结果清晰识别了研究区不同的岩性和地层；
使用Landsat-8 OLI、ASTER和ZY1-02D高光谱数据，分别采取不同的图像端元提取方法，在进行光谱分析的基础上，利用光谱角填图(SAM)即可得到研究区的主要矿物分类图件。通过野外验证，应用GIS技术进行集成和分析，修正相关图件后，便得到了精准的矿物分类综合图。研究表明：多源遥感数据的集成在岩矿识别方面效果良好、前景巨大。

近年来，遥感技术在地质构造分析、岩石矿物识别、蚀变异常提取、矿产资源勘探等方面应用日趋普遍和成熟(Sabins，1999；
张玉君等，2002；
李志忠等，2009；
王润生等，2011；
Peyghambari et al.，2021)。新疆东准噶尔卡拉麦里地区干旱的气候条件、戈壁滩型自然景观、良好的岩性露头给该地区利用遥感技术进行岩石矿物识别和开展相关矿产资源勘探工作提供了有利条件。国内外学者使用多光谱遥感数据在不同的地区开展了岩矿识别的大量研究工作(李培军等，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数据，采用光谱角填图法(Spectral Angle Mapper，SAM)对青海德尔尼蛇绿岩的主要岩石组成和蚀变矿物进行了探测，较有效地识别了蛇绿岩的主要岩性和相关矿物。杨伟光等(2018)在西藏罗布莎地区利用WorldView-2高空间分辨率遥感影像和ETM+数据对罗布莎铬铁矿床进行遥感控矿因素提取，使用了一系列光谱增强方法，评估遥感数据对铬铁矿成矿岩体地幔橄榄岩识别与圈定的能力和可行性，证明利用遥感技术探测潜在铬铁矿床的有效性。黄照强等(2010)、别小娟等(2013)也在西藏罗布莎地区利用遥感数据，采用不同的影像处理技术对蛇绿岩套和铬铁矿成矿相关的岩矿信息进行了提取和分析，并取得了不错的效果。Rajendran(2012)等在阿曼北部塞梅尔蛇绿岩地块使用Landsat TM和ASTER数据，选用去相关拉伸(Decorrelated Stretching，DS)、不同的波段比值(Band Rationing，BR)和主成分分析(Principal Component Analysis，PCA)等图像处理技术，成功圈定该地区蛇绿岩带内潜在的铬铁矿矿化带。Abdelaziz(2017)等在阿富汗洛加尔地块利用Landsat-8数据采用波段组合、主成分分析、波段比值和监督分类技术相结合的方法对洛加尔省铬铁矿有关的蛇纹岩进行了分离。Çörtük(2020)等基于ASTER和Landsat-8数据分析了土耳其安纳托利亚中部布尔巴斯蛇绿岩岩石分布，采用彩色波段比值、主成分分析和最小噪声分离变换(Minimum Noise Fraction，MNF)对遥感数据进行处理，有效地描绘了研究区域中辉石、变质岩和辉长岩的空间分布以及一批以前未被识别的蛇纹石化纯橄榄岩、辉石、变质岩底岩和片麻岩。

图 1 准噶尔盆地地质简图(a)(据Chen et al.，2004)和晚古生代额尔齐斯—卡拉麦里岛弧增生杂岩带构造单元划分略图(b)(据陕西省地质调查院滴水幅(L45C003004)1∶250, 000万区域地质调查报告修编，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)

高光谱数据在岩矿识别和矿产勘查中的应用也十分广泛(Kruse，1996；
Roy et al.，2009；
Bishop et al.，2011；
Govil et al.，2018；
魏红艳等，2020；
Govil et al.，2021)。李根军等(2020)应用ZY1-02D高光谱数据提取了大理岩、二长花岗岩等岩性信息，方解石、白云石等造岩矿物信息以及绿泥石、褐铁矿等蚀变矿物信息，与实测成果较为吻合，表明该数据对岩矿信息的识别效果较好。李娜等(2020)应用ZY1-02D高光谱数据在哈密遥感地质试验场地区主要对其在岩性-构造信息识别、矿物信息提取等方面开展了应用评价，清晰地显示和区分了不同地层、岩体、构造等地质体。连琛芹等(2020)应用GF-5高光谱数据在植被覆盖区进行了蚀变信息提取研究，提取的蚀变异常与研究区地质信息吻合较好。董新丰等(2020)，在甘肃柳园地区利用GF-5高光谱卫星数据开展矿物填图应用，对其在矿产资源方面的应用前景进行了初步评价。

由以上文献综述可知，多光谱数据在岩矿识别，特别是对铬铁矿成矿有关的蛇绿岩的提取中应用广泛、效果显著，而高光谱遥感技术在岩矿识别中也获得了比较成功的应用。但使用多光谱数据和高光谱数据相结合，以更好发挥其各自优势而实现岩矿识别的应用却不多。

笔者等以新疆东准噶尔卡拉麦里地区为研究对象，利用Landsat-8和ASTER(Advanced Spaceborne Thermal Emission and Reflection Radiometer，高级星载热辐射和反射探测器)多光谱卫星遥感数据和资源一号02D(ZY1E卫星)AHSI高光谱数据，采用连续最大角凸锥(SMACC)、最小噪声分离(MNF)、纯净像元指数(PPI)、n维可视化(n-D Visualizer)、光谱角制图法(SAM)等遥感图像处理算法，对研究区的岩石和矿物进行识别，并结合野外实地调查和分析，应用GIS软件进一步对遥感图像进行岩矿信息提取和空间展布分析。从而对研究区的岩矿进行分类，以为该地区下一步更大比例尺的地质调查工作和铬铁矿的勘探工作提供参考依据。

Ⅰ—阿尔泰微陆块；
Ⅱ2—晚古生代额尔齐斯—卡拉麦里岛弧增生杂岩带: —二台—托让库都克岛(额尔齐斯)弧滑脱—推断褶皱带, —托尔特库勒—库兰喀兹干周缘盆地滑脱—推覆带, —卡姆斯特—奥塔乌克塔什(卡拉麦里)岛弧褶断带, —卡拉麦里蛇绿混杂岩带, —卡拉麦里被动陆缘冲断—滑脱带, —芒能奥依—塔尔苏沟板内伸展盆地褶皱带；
Ⅲ—准噶尔微陆块: Ⅲ1—东准山前(间)断(坳)陷带Ⅰ—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

新疆卡拉麦里蛇绿混杂岩带位于准噶尔盆地东北缘，是中国西北部东准噶尔南缘一条重要的蛇绿混杂岩带(图1a)，它处于西伯利亚与哈萨克斯坦—准噶尔两大板块的交汇部位，是中亚造山带的一个重要构造单元(李锦轶等，1990、2009；
李现冰，2013；
陈新蔚等，2013；
樊立飞等，2014；
黄岗等，2017)。卡拉麦里蛇绿混杂岩带在构造单元上属于晚古生代岛弧增生岩带(Ⅱ2)，其受控于近NWW向的卡拉麦里区域性深大断裂(图1b)。

研究区主要地层为石炭系和泥盆系，仅在西南部出露少量的三叠系和志留系地层(图2)。石炭系地层以下石炭统姜巴斯套组为代表，其次是塔木岗组和上石炭统巴塔玛依内山组。其中姜巴斯套组可以分为两段(图2)，第一段主要为一套陆源碎屑岩、火山碎屑岩夹少量火山熔岩组合，主要岩性为凝灰质砾岩、凝灰质含砾长石岩屑粗砂岩、凝灰质粗砂岩、凝灰质粉砂岩、岩屑凝灰岩夹玄武质火山角砾岩、玄武玢岩等；
第二段为一套碎屑岩组合，岩性为浅灰色砂岩、细砂岩、长石岩屑砂岩、长石砂岩、凝灰质砾岩、凝灰质长石粗砂岩、岩屑粗砂岩。其在北部下与下石炭统黑山头组，上与上石炭统巴塔玛依内山组为不整合接触，南部与下石炭统塔木岗组为断层接触或与卡拉麦里蛇绿构造混杂岩带为不整合接触，上部与那林卡拉组为角度不整合接触。塔木岗组主要分布于卡拉麦里山南一带，其呈北西—南东向、近东西向展布，主要为一套陆源碎屑岩、火山碎屑岩夹火山熔岩组合，其下与泥盆系卡拉麦里组为不整合或整合接触，上部被上石炭统巴塔玛依内山组角度不整合覆盖。研究区出露的塔木岗组第一段岩性主要为暗灰色蚀变玄武玢岩、绿灰色流纹质英安斑岩、凝灰质细砂岩、安山玢岩、流纹岩、英安斑岩、灰紫色流纹斑岩、酸性角砾熔岩、紫灰色中酸性熔结凝灰岩、灰绿色细粒安山质凝灰角砾岩、灰黄色中酸性凝灰岩等。巴塔玛依内山组主要为火山熔岩夹火山碎屑岩组合，其角度不整合于下石炭统塔木岗组、姜巴斯套组等之上，或与下石炭统那林卡拉组为断层接触。泥盆系的地层以北塔山组分布最为广泛，其总体为一套火山碎屑岩、火山熔岩以及陆源碎屑岩组合，在研究区西部主要为火山碎屑岩、陆源碎屑岩夹少量熔岩组合，主要分布于卡拉麦里蛇绿构造混杂岩带及南部。

图 2 卡拉麦里地区地质矿产图及采样位置(据陕西省地质矿产局区域地质矿产研究院2009～2012年北塔山牧场幅(L46C003001)、滴水泉幅(L45C003004)1∶25万地质图修编)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—巴塔玛依内山组第一段；
C2h—弧形梁组；
C1n2—那林卡拉组第二段；
C1t1—塔木岗组第一段；
C1j2—姜巴斯套组第二段；
C1j1—姜巴斯套组第一段；
D2bt2—北塔山组第二段；
D2bt1—北塔山组第一段；
Dk—卡拉麦里组；
S3-D1h—红柳沟组；
S2b—白山包组；
C2ηγβb—中细粒似斑状黑云二长花岗岩；
C2ηγβa—中粗粒似斑状黑云二长花岗岩；
C2χργψc—中粒角闪石碱长花岗岩；
C2χργb—中细粒钠铁闪石碱长花岗岩；
C2χργa—中粒钠铁闪石碱长花岗岩；
D3γδο—灰白色中粗粒英云闪长岩；
C1δο—灰色细粒石英闪长岩；
D3γδο—中细粒英云闪长岩；
DCsi—硅质岩、硅质粉砂岩；
DCβ—块状与枕状玄武岩；
DCν—堆晶辉长岩；
DCσ—变质橄榄岩、辉石岩、蛇纹岩等；
δμ—闪长玢岩脉；
ν—辉长岩脉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

研究区的侵入岩主要是晚石炭世黑云二长花岗岩和钠铁闪石碱长花岗岩，早石炭世石英闪长岩和英云闪长岩。

卡拉麦里蛇绿混杂岩带位于研究区南部的卡拉麦里山一带，且主要分布于卡拉麦里大断裂的北侧，其主要由晚古生代蛇绿岩残片、蛇绿岩上覆岩系、外来岩片等构造块体和变形基质两部分组成(韩秉峻，2019)。蛇绿岩分布区内出露的地层主要是泥盆系和石炭系，岩石组合以陆源碎屑岩、火山碎屑岩和熔岩为特征(赵磊等，2019)，其代表性的地层单元是中泥盆统北塔山组和下石炭统姜巴斯套组，并且这这些地层大部分逆冲到蛇绿岩之上，而局部地段可见下石炭统姜巴斯陶组与下伏变质超基性岩呈不整合接触关系。蛇绿混杂岩带内各肢解蛇绿岩岩块多以断裂为界混杂堆叠，受剪切作用明显，岩块原生构造被后期构造改造或置换，但蛇绿岩的各组成单元出露仍较为齐全，基本上包括了典型蛇绿岩所有岩石类型。在卡拉麦里山东段红柳沟地区，组分发育不全，主要为不同程度蚀变的方辉橄榄岩、二辉橄榄岩、辉长岩、玄武岩等。中段清水一带组分发育相对较全，包括不同蚀变程度的方辉橄榄岩、含铬铁矿的纯橄岩、二辉橄榄岩、单斜辉石岩、堆晶辉长岩、角闪岩、辉绿岩、斜长岩、斜长花岗岩、枕状玄武岩、块状玄武岩、放射虫硅质岩等。西段由于受露头的影响，可见组分有蛇纹岩、蛇纹石化橄榄岩、玄武岩、辉长岩、斜长花岗岩和硅质岩等。基质部分主要为泥盆—石炭纪的变凝灰泥质长石粉砂岩、变晶屑岩屑凝灰岩、凝灰岩、千枚岩、糜棱岩化凝灰岩、糜棱岩等。研究区的矿产以铬铁矿为主，此外还有少量的金矿、赤铁矿—褐铁矿分布(张峰等，2014)(图2)。铬铁矿矿床和矿点产于蛇绿岩上部的斜辉橄榄岩、纯橄岩中，但多集中在纯橄岩，这些岩石均已蚀变为蛇纹岩或蛇纹石化橄榄岩。

2.1 遥感数据

2.1.1 数据简介

本文根据研究目的和制图需求，选用Landsat-8、ASTER多光谱遥感数据以及资源一号02D(ZY1E卫星)AHSI高光谱数据，现对不同的卫星传感器性能特征介绍如下。

Landsat-8 OLI数据包含9个波段，其中Band 1～7和Band 9的空间分辨率为30 m，Band 8全色波段空间分辨率为15 m，成像幅宽为185 km × 185 km。本文使用了一景Landsat-8数据，其成像时间为2021年4月14日，其条带号和行编号分别为141/29，云量覆盖为0.28。ASTER数据含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数据(表1)。资源一号02D AHSI高光谱数据，空间分辨率为30 m，波段数量为166个，光谱范围为0.4～2.5 μm，成像幅宽为60 km × 60 km。光谱分辨率在可见光—近红外为10 nm、在短波红外为20 nm。本次使用的数据景号是ZY1E_AHSI_E89.90_N45.01_20201027_005899_L1A0000192836，成像日期为2020年10月27日，云量覆盖为3.6%。

表1 卡拉麦里地区覆盖的4景ASTER数据基本信息Table 1 Basic information of ASTER data of 4 scenes covered in the Karamaili area

2.1.2数据预处理

本次获取的Landsat-8数据为L1TP级，对Landsat-8 OLI数据进行辐射定标、FLAASH大气校正(主要参数见表2)。

表2 不同遥感数据FLAASH大气校正主要参数Table 2 Main parameters of FLAASH atmospheric correction for different remote sensing

由于ASTER的VNIR和SWIR分别处于不同的传感器子系统，故通过“Layer Stacking(波段叠加)”将每景的Band 1～9放在同一多光谱数据集中，并将SWIR的6个空间分辨率为30 m的波段重采样至15 m，然后依次进行辐射定标、FLAASH大气校正(主要参数见表2)。对高光谱数据ZY1-02D AHSI首先进行辐射定标，然后将多个波段拆分成单波段TIFF格式以便后续处理。ZY1-02D高光谱短波红外范围内数据条纹现象明显(图3)，故对图像的噪声特征进行分析后，使用傅立叶变换频谱分析和反傅立叶变换的方法，采用Python编程的方法去除影像中的条纹(代码见附录)。将去除条纹的数据通过波段叠加(Band Stack)重新组合成166个波谱数据，然后选用FLAASH大气校正模型对其进行校正(主要参数见表2)。以Landsat-8的全色波段为参考影像对ZY1-02D高光谱进行几何校正。

图 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 岩矿识别技术方法

利用遥感数据识别地层岩性和矿物信息时，由于不同传感器光谱及空间分辨率的差异，致使对目标地物的识别存在明显的专属性，因此笔者等选用Landsat-8 OLI、ASTER和ZY1-02D AHSI 3种遥感数据，配合开展地层岩性和矿物识别，以期达到对研究区的岩石地层和相关矿物进行识别分类的目的(图4)。

图 4 多源遥感数据地层岩性和岩石矿物识别流程Fig. 4 Identification process of formation lithology, rock and mineral from multi-source remote sensing data

2.2.1数据降维

为有效识别出遥感影像海量信息中的有效数据，需首先对遥感影像进行降维(Dimension reduction)处理，即将高维数据转换为低维数据而不丢失信息的过程称为降维(Dimension reduction)。目前，最常用的数据降维方法包括主成分分析法和最小噪声分离法。

主成分分析(Principal Component Analysis，PCA)是一种图像增强技术，它将高度相关的原始图像数据转换为一组主成分或PC波段的不相关变量，通过抑制在所有波段占主导地位的辐照度效应来分离噪声成分，从而达到降低数据维度的目的。

本次研究对Landsat-8 OLI数据进行PCA变换之后，显示前3个主成分包含了所有波段中99.54 %的信息(表3)，表示了数据集中的域变化，并突出显示了所有波段的最常见特征。因此对前3个主成分PC3、PC2、PC1进行假彩色合成，并结合地质图和野外勘查便得到了研究区的地层岩性信息(表4)。

表3 Landsat-8 OLI主成分分析结果统计表Table 3 Statistics of Landsat-8 OLI principal component analysis results

表4 Landsat-8 OLI主成分分析PC3、PC2、PC1假彩色合成解译卡拉麦里地区主要地层岩性特征表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)是将一幅多波段图像的主要信息集中在前面几个波段中，主要作用是判断图像数据维数、分离数据中的噪声,减少后处理中的计算量。它由Green等(1988)开发，该方法是由两个独立的主成分分析旋转组成的线性变换。第一次旋转利用噪声协方差矩阵的主分量对数据中的噪声进行去相关和重缩放，使得有噪声的数据具有单位方差，并且没有波段间相关性。在第二阶段，通过分析特征值和相关图像来计算固有维数。因此，与以噪声为主的图像相关的单位特征值可以分离出来，从而增强光谱处理结果。此次研究对ASTER的9个波段和ZY1-02D AHSI经过筛选的152个波段进行MNF变换，以达到降维、抑制噪声和信息增强的目的。

2.2.2图像端元提取

混合像元广泛存在于图像中，影响遥感图像的分类精度以及目标探测效果(蓝金辉等，2018)。若一个像素对应的地面区域内只包含一种特征地物，则称此像素为纯像元，而端元即为数据中代表类别特征的理想化纯像元，提取遥感数据中的纯像元技术称为光谱端元提取(齐滨，2012；
陈晋等，2016)。常用方法有基于连续最大角凸锥技术、纯净像元指数技术和n维可视化技术(王茂芝等，2015)。

基于连续最大角凸锥(sequential maximum angle convex cone，SMACC)使用凸锥模型表示矢量数据，直接从多光谱和高光谱数据中提取图像端元。它通过识别多维空间中的极值点，并形成一个包含第一个端元的凸锥。然后，将约束斜投影应用于现有圆锥体，圆锥体将增加以包括下一个端元。简单来说，SMACC首先查找最亮的像素，然后查找与最亮像素差异最大的像素，然后查找与前两个像素差异最大的像素。该过程将继续进行，直到投影衍生出一个已存在于凸锥内的端元或达到指定数量的端元。在这项研究中，对经过大气校正后的Landsat-8 OLI数据使用SMACC方法进行图像端元自动提取，端元波谱提取数量设置为10，误差容限值保持默认值0，分离端元波谱的约束条件选择“Posivity only”，即把每个波长的正值端元波谱作为约束条件。在合并相似端元波谱后，共提取出6个图像端元的波谱。

纯净像元指数( Pixel Purity Index，PPI)和n维可视化(n- Dimensionsl Visualizer )用于收集ASTER和ZY1-02D AHSI图像中的端元波谱。通过计算多波谱和高光谱图像的纯净像元指数(PPI)，进而在图像中寻找波谱最“纯”的像元。这可以通过连续收集n维散点图中的极值来确定(Kruse，1998；
Govil et al.，2021)。n-D Visualizer用于定位、识别、聚集数据集中最纯的像元，从而获取纯净的端元波谱，以生成光谱库文件，进而对地物进行分类。因为经过MNF变换后前几个波段将包含90%以上的信息，所以选取ASTER的前三个波段，ZY1-02D AHSI的前十个波段进行PPI的提取。在执行PPI时，ASTER和ZY1-02D AHSI的迭代次数均设置为10000，前者的阈值系数设置为2.5，后者为3.0。针对ASTER和ZY1-02D AHSI数据分别提取了7个和9个图像端元的波谱。

2.2.3光谱分析

光谱分析(Spectral Analyst)根据矿物的光谱特征来进行矿物识别。此次研究使用ENVI软件中的光谱分析算法，如二进制编码(Binary Encoding)、光谱角填图(Spectral Angle Mapper)和光谱特征拟合(Spectral Feature Fitting)，对提取的图像端元光谱与USGS波谱库典型矿物波谱进行匹配，从而识别出端元光谱的目标矿物种类。将每种匹配算法的权重设置为1，则总得分为3，总得分越靠近3，说明图像端元光谱与典型矿物波谱匹配度越高(表5)。研究区几种典型矿物图像端元波谱和USGS标准矿物波谱库波谱对比见图5。

表5 卡拉麦里地区不同遥感数据矿物识别时光谱分析匹配得分表Table 5 Matching scores of spectral analysis for mineral identification of different remote sensing data in Karamaili area

图 5 卡拉麦里地区典型矿物图像端元波谱(红色)和标准波谱库波谱(蓝色)对比图。

OLI数据：(a)橄榄石、(b)黑云母、(c)绿泥石；
ASTER数据(d)蛇纹石、(e)辉石、(f)绿泥石；
AHSI数据：(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)本质上是一种监督分类技术，它将图像光谱与参考光谱相匹配，在多维特征空间中确定图像光谱和参考光谱之间的角度。角度越小，说明图像光谱和参考光谱之间的相似性越高。光谱角的值介于0到π/2之间，根据Kruse等人(1993)给出的公式计算：

(1)

其中θ是图像光谱和参考光谱之间的角度，n是光谱波段数；
t是实际光谱(图像光谱)的反射率，r是参考光谱的反射率。

本次研究对Landsat-8 OLI经过SMACC提取端元波谱和ASTER、ZY1-02D AHSI经过MNF、PPI和n维可视化提取端元波谱后，分别对三种数据进行光谱分析匹配，然后采用SAM进行研究区矿物的识别制图(图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.1 OLI岩性地层制图和分析讨论

经过PCA(主成分分析，principal component analysis)变换后，图像色彩丰富、可获取信息量大幅度提升。根据不同地层岩性在主成分分析后的假彩色合成影像(R：PC3，G：PC2，B：PC1)所表现出来的颜色、纹理、形状、水系等特征，参考地质资料便可识别出研究区的全部地层和绝大多数岩性(表5)。但是由于蛇绿岩混杂岩构成的复杂性、OLI数据光谱分辨率和空间分辨率的限制，并不能识别出堆晶辉长岩、块状玄武岩等小型的岩浆岩和硅质岩等沉积岩。对岩性组成非常相似的地层，如北塔山组第二段和北塔山组第一段中的砂岩和粉砂岩，需要经过野外勘查才能正确划分其组分。

3.2 Landsat-8 OLI、ASTER和ZY1-02D AHSI矿物识别制图和分析讨论

本研究针对Landsat-8 OLI大气校正之后的图像，使用SMACC方法提取图像端元波谱，针对ASTER和ZY1-02D AHSI大气校正、条纹去除、MNF处理之后的图像，采用PPI和n维可视化工具提取端元波谱，并以USGS标准矿物波谱库的波谱为参考，进行光谱分析以确定提取的端元波谱所代表的矿物种类(表5)，最后把端元波谱作为已知波谱，采用SAM方法和图像波谱进行比对分析，以达到矿物分类识别(图6)的目的。

由矿物识别分类图可知，不管是多光谱数据Landsat-8 OLI和ASTER，还是高光谱数据ZY1-02D AHSI，对于一些造岩矿物，如长石、石英、黑云母的识别特别有效，但是由于这些矿物往往具有相似的光谱特征，因此本研究只能局限于将这些矿物的组合识别出来，如黑云母、石英、长石组合，钾长石、石英组合。此外由于蛇绿岩混杂岩构造的复杂性、“异物同谱”与“同物异谱”现象的存在以及遥感本身获取的就是混合像元，因此并不能单独识别某一种矿物，仅能识别出橄榄石和辉石、蛇纹石和斜方辉石的矿物组合及其它一些矿物组合。Landsat-8 OLI和ZY1-02D AHSI对绿泥石的识别效果好，ASTER对橄榄石和辉石组合、蛇纹石和高岭石识别效果佳。而ZY1-02D AHSI由于高的光谱分辨率，不管是矿物的识别种类还是效果都优于前两种多光谱数据。但是它却未能将方解石、白云石等Landsat-8 OLI和ASTER可识别的矿物鉴定出来。因此将多源遥感数据进行集成更有利于研究区矿物的识别，不同的遥感数据源可起到互相印证和互补改正的作用。

3种不同数据识别出的矿物分布范围差异可能

是由于遥感数据本身的特性引起的，如Landsat-8 OLI和ZY1-02D高光谱数据虽然空间分辨率都是30m，但是后者的光谱分辨率远远高于前者。这就意味后者对矿物之间微小差距的探测更有效。而前者由于光谱分辨率低，对一些波谱特征在局部变化不大的矿物识别更有效。当然还有一些光照差异、地形起伏等因素的影响也会造成不同数据识别矿物性能的差异。

3.3 野外采样及验证

在对研究区的地层岩性分布状况进行深入了解的基础上，结合Landsat-8 OLI数据PCA方法识别出的岩性地层特征，以及Landsat-8 OLI、ASTER和ZY1-02D AHSI数据SAM矿物识别分类图(图6)，进行了野外勘查和采样工作。根据野外实地调查，Landsat-8 OLI岩性地层分类结果和3种数据的矿物分类结果都得到了很好的验证，但是由于“异物同谱”和“同谱异物”现象的存在以及遥感器性能的差异，分类时仍然将一些岩性矿物信息错归或遗漏。本次研究共采集了80个岩石样品(图2)，并将采样点作为实际验证点，结果发现69个岩石样品的野外鉴定命名和不同遥感数据矿物识别分类结果相吻合，其准确度达到了86.25%。研究区典型的岩石标本照片见图7。

图7 野外勘查采集的岩石标本：(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 矿物识别分类结果集成与分析

经过3.2不同遥感数据矿物识别结果的讨论和分析可知，Landsat-8 OLI、ASTER和ZY1-02D AHSI 3种数据对研究区不同的矿物识别都起着至关重要的作用，因此考虑将多源遥感数据分类结果进行集成与分析，以达到最优的矿物识别分类效果。本研究中，应用GIS软件对3种数据的矿物分类结果以及野外实地勘查验证资料进行综合分析和数据集成。首先将遥感栅格分类图转换为shape矢量格式以便于进行空间分析运算。随后以地层岩性资料为基础、以野外实地勘查采样结果为标准，对多源遥感数据矿物分类图进行空间求交和合并运算以筛选出正确的矿物分布信息，对矿物分类图件进行调整(图8)。

把多源遥感数据岩矿识别分类结果和野外地质资料进行集成分析，可弥补一些漏掉的矿物信息，并纠正一些错误的矿物分类结果。主要是对Landsat-8 OLI识别的橄榄石的分布范围进行了调整，对将一些黏土矿物而被错误识别为白云石的范围进行了删减。对ASTER识别的蛇纹石和斜长石范围进行调整，经野外验证，ASTER提取的蛇纹岩信息多为蛇纹石化的橄榄岩，而一少部分为橄榄岩。对ZY1-02D AHSI识别的高岭石、伊利石、蒙脱石等黏土矿物范围进行了调整。经过修改调整后的矿物分类综合图清晰地显示了研究区主要矿物的空间展布范围，特别是对超基性岩区域矿物的构成有了清晰的认知，这对该区域与超基性岩有关的铬铁矿的勘探工作提供了重要的矿物分类图件。

笔者等将多光谱遥感数据Landsat-8和ASTER，以及ZY1-02D高光谱数据进行集成，在新疆东天山卡拉麦里地区，针对不同数据进行相应预处理后，对Landsat-8 OLI采用PCA方法进行地层岩性识别；
对Landsat-8 OLI、ASTER和ZY1-02D AHSI数据使用MNF方法进行数据降维和分离噪声后，相应采用SMACC、PPI结合n维可视化两种提取图像端元波谱方法，在光谱分析后，使用SAM方法进行研究区矿物识别。并在野外地质勘查的基础上，对多源遥感数据矿物分类结果进行集成分析而得到研究区精准的矿物综合分类图。结论如下：

(1)在利用遥感技术进行岩矿信息提取时，特别是矿物识别时，应充分发挥多源遥感数据之间互相印证和补充细化的作用，以达到更精准的矿物分类识别结果。

(2)在对卡拉麦里地区岩矿信息提取的过程中，经野外实地勘查，不同的数据源都达到了较好的效果，3种数据对一些造岩矿物的组合识别效果都很好，但是每种数据都有其自己的独特的优势和作用，Landsat-8识别出的白云石、ASTER识别出的橄榄石和辉石组合，ZY1-02D识别出的蛇纹石、斜方辉石组合都是其它两种数据所未识别出来的。

(3)图像端元波谱的提取在矿物识别中起着至关重要的作用，针对不同的数据源分别采用SMACC、PPI结合n维可视化来提取端元波谱得到了较为精准的矿物识别分类图件。

(4)岩矿识别工作应充分发挥“3S”技术的优势，以遥感(RS)为数据源，在综合地质资料和遥感成果的基础上，发挥全球定位系统(GPS)的优势，在野外进行实地勘查验证，在地理信息系统(GIS)软件中，对不同数据源提取的岩矿信息、地质综合资料、野外勘查验证结果进行集成和综合，从而得到精准的岩矿分类识别结果。

附录：ZY1-02D高光谱数据条纹去除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:该波段的数据，narray格式。len(narray)是行数，len(narray[0])列数

"""

image = gdal.Open(tifpath) #打开该图像

if image == None:

print(tifpath + "该tif不能打开!")

return

lie = image.RasterXSize #栅格矩阵的列数

hang = image.RasterYSize #栅格矩阵的行数

im_bands = image.RasterCount #波段数

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数据类型包括

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:输入数据

:param im_proj:投影信息

:param im_geotrans:放射变换信息

:return: 0

"""

if ′int8′ in im_data.dtype.name: #判断栅格数据的数据类型

datatype = gdal.GDT_Byte

elif ′int16′ in im_data.dtype.name:

datatype = gdal.GDT_UInt16

else:

datatype = gdal.GDT_Float32

#判读数组维数

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维

#创建文件

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)

#判断结果

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′)