基于Landsat 8数据的弓长岭矿区遥感蚀变异常信息提取

宋 坤 王恩德 付建飞 姚玉增 孙 健

(1.东北大学资源与土木工程学院,辽宁 沈阳 110819;2.中国海洋大学海底科学与探测技术教育部重点实验室,山东 青岛 266100)

辽宁鞍本地区弓长岭铁矿是我国唯一由鞍山式贫铁矿经后期热液改造形成的大型磁铁富矿,也是唯一具有工业开采价值的富铁矿[1]。国内外诸多学者在弓长岭铁矿成因方面进行了深入研究,其找矿方法以地质、地球物理及地球化学方法为主,对于遥感方法的应用较少。自20世纪70年代以来,遥感技术以其快速、无损的优势被广泛应用在地质领域[2]。相对于较传统的地质找矿方法在矿产勘探工作中的高成本、低效率等不足来说,遥感技术可在较大范围内提取研究区蚀变矿化信息,为矿产资源潜力评估提供帮助,为找矿勘探工作提供方向[3]。本研究以弓长岭铁矿为例,选取Landsat8OLI遥感数据,以羟基矿物和铁染矿物在特定波段的吸收或反射特征为理论依据,应用主成分分析法和波段比值法对遥感异常蚀变信息进行提取。通过总结并分析蚀变信息与铁矿体之间的关系,为鞍本地区及其他同类型矿床的找矿工作提供依据。

1 研究区概况

1.1 研究区地质背景

研究区位于弓长岭及其周边区域,大致以弓长岭为中心区域,向西至辽阳市东侧,向东至南芬铁矿区,呈70 km×50 km的矩形,面积约3 500 km2(图 1)。弓长岭铁矿资源十分丰富,具有储量大、分布广等特点。研究区位于华北克拉通北缘东段的辽东地区,具有下部太古宙基底与上部古生界沉积盖层的双层结构。基底岩性主要由鞍山群茨沟组、大峪沟组和樱桃园组变质岩组成[4]。我国赋存于太古宙变质岩系的铁矿最早发现于这一带,称之为鞍山式铁矿,以贫矿多、少有富矿、储量大为主要特点。据统计,辽宁省铁矿资源储量占全国总量的1/4左右,且鞍山式铁矿占辽宁省铁矿总量的96%。

图1 研究区区域地质特征Fig.1 Regional geological characteristics of the study area

1.2 矿床地质特征

鞍本地区以弓长岭矿区规模最大,也最为典型。弓长岭矿区主要为沉积变质型铁矿,磁铁富矿体一部分分布在各层条带状铁矿床中,还有部分含铁矿体产在蚀变围岩或混合岩中[5]。这表明富铁矿体可能是条带状磁铁矿体经热液改造形成[6]。矿集区呈NW—SE向展布,弓长岭一、二、三矿区以及独木矿区为铁矿床的主要分布区域(图2)。鞍山群茨沟组是该地区的主要赋矿地层,二矿区茨沟组地层出露最为完整。弓长岭背斜为控矿构造[7]。矿区内矿石类型主要为磁铁石英岩,其中二矿区内含有富铁矿石(图2)。出露的岩浆岩主要为晚太古代混合花岗岩,弓长岭矿区呈孤岛状分布在混合花岗岩中[8]。

图2 弓长岭矿区地质特征Fig.2 Geological characteristics of Gongchangling mining area

2 遥感数据处理

2.1 数据源

Landsat 8卫星于2013年成功发射,卫星上搭载了OLI陆地成像仪和TIRS热红外传感器,其空间分辨率和光谱分辨率等方面与之前的Landsat系列卫星基本保持一致。Landsat 8与Landsat 7卫星波段对应关系参见表1。Landsat 8遥感影像数据分辨率较高,且影像涉及范围广,在遥感异常信息提取中被广泛使用;同时,OLI和ETM+数据波段分配大体一致,且OLI波段范围比ETM+更窄,可更好地表现出蚀变矿物波的谱特征[9]。

表1 Landsat 8与Landsat 7卫星波段参数Table 1 Satellite bands parameters of Landsat 8 and Landsat 7

通过对TM和ETM+数据波段特征及其应用的分析研究,总结出OLI传感器中不同波段的光谱特征及其适用领域,即Band 5、Band 6和Band 7 3个波段与蚀变异常提取相关,其中,对碳酸盐以及黏土矿物等有明显吸收特征的光谱波段为Band 7。

研究区Landsat 8 OLI图像的成像时间为2015-04-20,卫星轨道号为119-31,经反复比对分析,该时段图像地物轮廓清晰、层次感强,图上无云层,地面无积雪,植被覆盖率低,地表岩石、土壤裸露较好,成像质量最佳,有效满足了本次研究的要求(图3)。研究区图像范围为东经 123°04′40″～123°58′01″,北纬40°53′14″～41°25′54″。

图3 研究区原始遥感影像Fig.3 Original remote sensing image of the study area

2.2 遥感影像预处理

在遥感成像过程中,遥感影像因受其搭载平台、传感器误差以及大气环境、太阳高度角、地形等不同因素影响,得到的测量值与地物本身的光谱反射率存在误差。为了消除或修正辐射畸变,得到真实的地物反射率,须对原始影像数据进行辐射校正,通常包括传感器校正和大气校正两类[10]。

对传感器误差的校正是通过辐射定标进行,将传感器中的电压或数字量化值(Digital Number,DN)换算为大气外层的表面反射率。大气校正是通过将大气顶层的辐射亮度值转换为地表反射率消除或修正误差,从而获取地物的真实反射率(图4)。上述处理均使用ENVI4.7软件的辐射定标和FLAASH大气校正模块,参照遥感数据的参数设置完成校正。

图4 大气校正前后遥感图像及其波谱曲线Fig.4 Remote sensing images and their pop curves before and after atmospheric correction

3 研究方法

不同的矿物、岩石及矿化蚀变信息都具有其特有的光谱特征[11]。依据蚀变矿物的不同波谱特征,可以有效提取遥感蚀变异常信息[12]。通常情况下,蚀变异常信息提取的主要相关矿物离子包括铁离子、羟基或碳酸根等离子(表2)。

表2 对岩石反射光谱特征起主导作用的离子和基团的光谱特征[11]Table 2 Spectral characteristics of ions and groups that dominate the spectral characteristics of rock reflections

分析表2可知:

(1)含铁(Fe2+、Fe3+)离子。以含Fe3+的褐铁矿、赤铁矿等蚀变矿物为主,其中含Fe2+的蚀变矿物相对较少。如图5所示,Fe3+矿物在0.40～0.55μm和0.85～0.95μm波段表现为明显的吸收特征,对应OLI 2(0.45～0.51μm)波段、OLI 3(0.53～0.59μm)波段和OLI 5(0.85～0.88μm)波段的强吸收特征,以及在OLI 4(0.64～0.67μm)波段的强反射特征。

图5 含铁类蚀变矿物反射光谱曲线Fig.5 Reflectance spectral curves of iron dyeing altered minerals

(2)含羟基(OH-)或碳酸根(CO32-)离子。该类离子蚀变矿物主要为黏土矿物,如图6所示,在1.4 μm及2.2～2.4μm波段存在吸收特征,对应OLI 7(2.11～2.29μm)波段产生低反射值,在 OLI 6(1.57～1.65μm)波段产生高反射值。含碳酸根矿物主要存在5个特征吸收谱带,分别在2.35μm、2.55μm波段处表现出较强吸收,在1.9、2.0、2.16 μm波段处表现出相对较弱吸收,对应OLI 7(2.11～2.29μm)波段形成吸收特征。

图6 含羟基类蚀变矿物反射光谱曲线Fig.6 Reflectance spectral curves of hydroxyl altered minerals

3.1 波段运算法

波段运算法包括波段的加减运算与比值运算。波段加法通过将多幅图像亮度值相加求平均,达到图像视觉增强效果、边缘增强效果。图像减法运算则可以去背景,突出研究对象。蚀变围岩与周围岩石之间的亮度值存在差异,可通过波段加减法运算增强矿化蚀变信息[13]。

为了增强蚀变信息,还可以运用波段比值法。根据每个象元在两个不同波段上的亮度比值形成新的比值图像[14],可消除地物反射或避免光照引起的阴影干扰,突出重要信息。

比值法通常包括简单比值法、差和比值法以及截取比值法等。结合弓长岭矿区实际情况,突出异常信息,本研究运用简单比值法对图像进行处理[15]。简单比值法运算公式为

式中,DNm(x,y),DNn(x,y)分别为像元 (x,y)在m和n波段上的亮度值;Rmn(x,y)为输出的亮度比值。

OLI遥感影像通常可对铁氧化物、氢氧化物类,羟基矿物类,碳酸盐、泥化类3类蚀变矿物进行识别[16]。在遥感异常信息提取过程中增强蚀变信息的主要OLI波段比值为:①Band 4/Band 2,主要增强铁氧化物、氢氧化物类蚀变;②Band 6/Band 4,主要增强硅化类蚀变;③Band 6/Band 7,主要增强碳酸盐化、泥化类蚀变。

本研究利用ENVI 4.7软件中的“band math”工具进行上述相关波段进行运算,结果如图7所示。

图7 波段比值法图像处理结果Fig.7 Image processing results by band ratio method

3.2 主成分分析法

主成分分析法(Principal Components Analysis,PCA)是较常用的蚀变信息提取方法。通过对反应蚀变信息的多个波段进行特征统计的多维正交线性变换(K-L变换),生成1组新的组分图像[17]。主要目的是集中多波段中有用的信息到尽可能少的新组分图像中,确保地物信息之间互不干扰再去相关性[18]。同时,选用计算结果中的新波段作为多波段原始图像,可实现增强重要信息的目的[19]。

主成分分析法在遥感蚀变异常信息提取中应用广泛,被称为CROSTA法。根据前人研究成果,对于OL I数据,常用OLI 2、OLI 5、OLI 6、OLI 7这4个波段提取羟基蚀变异常,用OLI 2、OLI 4、OLI 5、OLI 6这4个波段提取铁染遥感蚀变异常。

铁染异常主分量的判断准则为含铁染异常矿物OLI 4与OLI 2及OLI 5的系数符号相反、OLI 4与OLI 6的系数符号相同[20]。同时,还可以根据不同波段在不同主分量上的载荷情况来辅助确定异常主成分(该主分量在OLI 4波段载荷较大,在OLI 2波段载荷较小)。本研究将与上述准则对应的主分量称为铁染异常主分量。

羟基异常主分量的判断准则为含羟基异常矿物的OLI 6系数应与OLI 7及OLI 5的系数符号相反,OLI 2一般与OLI 6系数符号相同[21-22],且该主分量在OLI 6、OLI 7波段载荷较大。本研究将与上述准则对应的主分量称为羟基异常主分量。

3.2.1 铁染蚀变异常信息提取

根据上述PCA方法提取原理,对研究区图像的OLI 2、OLI4、OLI 5、OLI6波段进行主成分分析,得到4个主分量,并对其进行统计分析,结果见表3。

表3 OLI 2、4、5、6波段提取铁染异常特征向量矩阵及各主分量Table 3 Eigenvector matrix and eigenvalues of each principal component of iron dyeing anomaly with OLI 2,4,5 and 6

表3主要体现了4个波段的亮度信息,其特征向量值在4个波段中均为负值,不符合判断准则。同时,根据上述铁染异常主分量判断准则,发现并没有完全符合条件的主成分。但是第4主成分(PC4)最接近判断准则,其在 OLI 4上的载荷系数(0.609 854)与在OLI 2上的载荷系数(-0.782 576)符号相反,而在OLI 5上的载荷系数(0.002 123)虽然是正,但数值很小接近于0。考虑到获取图像和预处理时可能会带入一些误差,因此可以将PC4视为满足判断准则的主分量。从含铁矿物光谱特征来看,含铁蚀变岩石在OLI 2强吸收、OLI 4强反射,因此含铁染信息的主分量在OLI 2和OLI4显示为强吸收和强反射且在这两个波段上具有较大的载荷,可以看出PC4主分量满足判断准则,代表了铁染蚀变信息。

通过分析PC4灰度图像发现,与实际矿点对比发现异常区域为低值(显示为暗),因此对PC4数据取相反数,使得铁染异常区显示为高值(亮),取反后的-PC4灰度图像如图8所示。

图8 铁染异常主分量DN值取反图像Fig.8 Principal component DN value inverse image of iron dyeing anomaly

3.2.2 羟基蚀变异常信息提取

根据上述PCA方法提取原理,对研究区图像的OLI 2、OLI5、OLI6、OLI7波段进行主成分分析,得到4个主分量,并对其进行了统计分析,结果见表4。

表4 OLI 2、5、6、7 提取羟基异常特征向量矩阵及各主分量Table 4 Eigenvector matrix and eigenvalues of each principal component of hydroxyl anomaly with OLI 2,5,6 and 7

由表4可知:根据上述羟基异常主分量判断准则,该主分量在OLI 6、OLI 7波段载荷较大。同时,可以发现 PC4主分量符合判定准则,其 OLI 6(0.617 392)系数与OLI 7(-0.543 045)及OLI 5(-0.331 771)系数符号相反,OLI2(0.462447)与OLI 6(0.617 392)系数符号相同。PC4特征向量载荷系数绝对值较大的是OLI 6和OLI 7,分别为0.617 392和0.543 045,二者符号相反,而其他波段载荷较小,表明PC4主成分的信息主要来自OLI 6和OLI 7,与黏土矿物和碳酸盐矿物存在有关,且与OLI 6高反射和OLI7强吸收相符,突出了泥化蚀变信息。因此,本研究选择PC4作为羟基信息。结合研究区实际情况发现灰度图像中亮的区域为羟基异常区(图 9)。

图9 羟基异常主分量图像Fig.9 Principal component image of hydroxyl anomaly

4 弓长岭矿区蚀变异常信息提取

4.1 铁染异常

首先对OLI 2、4、5、6 主成分分析得到的铁染异常主分量-PC4数据进行数学统计(表5),利用门限法分级,多次试验确定最优分割(表6),得到如图10所示的铁异常分级结果。

表5 铁染异常主分量-PC4数据数学统计信息Table 5 Mathematical statistics information of the principal component-PC4 data of iron dyeing anomaly

表6 铁染异常主分量异常分级Table 6 Classification of principal component of iron dyeing anomaly

图10 弓长岭地区铁异常图Fig.10 Iron anomaly map of Gongchangling Area

4.2 羟基异常

对OLI 2、5、6、7主成分分析得到的羟基异常主分量PC4数据进行数学统计(表7),利用门限法分级,多次试验确定最优分割(表8),得到如图11所示的羟基异常分级结果。

表7 羟基异常主分量PC4数据数学统计信息Table 7 Mathematical statistics information of the principal component PC4 data of hydroxyl anomaly

表8 羟基异常主分量异常分级Table 8 Classification of principal component of hydroxyl anomaly

图11 弓长岭地区羟基异常图Fig.11 Hydroxyl anomaly map of Gongchangling Area

针对弓长岭矿区范围进行详细的蚀变异常分级,观察其在小范围内(矿区)是否能表达出更多、更有规律的信息。将主成分分析得到的铁染和羟基异常主分量进行裁剪,只保留弓长岭矿区图像。对裁剪得到的图像数据重新进行数学统计,然后利用门限法重新分级[23]。蚀变信息特征中大部分矿区以外的背景值、假异常等干扰因素可通过处理去除,结果如图12和图13所示。

图12 弓长岭矿区铁异常图Fig.12 Iron anomaly map of Gongchangling Mning Area

图13 弓长岭矿区羟基异常图Fig.13 Hydroxyl anomaly map of Gongchangling Mining Area

4.3 分析验证

主成分分析法对蚀变信息的提取效果比较理想,研究区内主要矿区都有异常显示,能明显的与非异常区分。铁和羟基的蚀变异常区域重叠较好,也吻合了研究区内矿化区域的铁和羟基含量均较高的特征。为了进一步验证遥感蚀变异常信息与区域矿化信息的吻合程度,本研究将铁和羟基蚀变异常进行叠加处理,并将研究区内已知矿点投影到异常叠加图上,分析蚀变异常与矿点的吻合程度,结果如图14所示。

图14 蚀变异常叠加图和部分已知矿点投影Fig.14 Alteration anomaly superposition map and partial known ore occurrences projection map

由图14可知:投影的7个超大型矿床和6个大型矿床中有5个超大型矿床和4个大型矿床位于强蚀变异常区范围内,吻合程度较好,2个大型矿床处于轻微蚀变异常区范围内,可能是矿床围岩蚀变作用较弱或为隐伏矿床所致。2个超大型矿床不在明显的蚀变异常区内,吻合程度较差,分别为大台沟矿床和思山岭矿床。大台沟矿床2011年详查显示埋深达到1 100～1 400m,含铁岩系分布在太古宙鞍山群,上覆辽河群及近水平的青白口系、震旦系等沉积盖层[24]。思山岭铁矿埋深达400 m以上,赋矿岩系为鞍山群茨沟组,上覆沉积盖层为青白口系和震旦系岩层[25]。两个矿床的变质蚀变信息均被沉积盖层覆盖,影响蚀变信息提取结果。除了以上两个超大埋深的隐伏矿床外,其余矿床与蚀变异常区域基本吻合,证明主成分分析法能有效地提取矿化蚀变异常信息,在同类型矿床中具有较好的应用前景。

5 结 论

(1)以Landsat 8遥感数据为基础,通过提取弓长岭矿区蚀变矿物信息,在同类型BIF铁矿床中建立了一种快速准确的遥感找矿勘查方法。

(2)采用主成分分析法和波段比值法分别对研究区以及弓长岭矿区进行蚀变异常信息提取,在研究区内主成分分析法可以有效地抑制干扰信息,提取铁和羟基蚀变信息效果较好。弓长岭矿区铁异常主要分布在矿区中心且向南延展。羟基异常范围较铁染异常面积大,且矿区南部有较多弱异常分布,其中有部分信息为地表覆盖下的假异常。

(3)提取的蚀变异常信息与研究区内已知矿点对比发现,异常区与矿点吻合较好。蚀变信息集中区域为对找矿有指示意义的区域,反映出本研究方法对于鞍本地区及其他同类型矿床研究有一定的参考价值。但该方法在提取地表覆盖层较厚的隐伏矿床蚀变信息过程中,对假异常的识别和蚀变弱信息的增强仍存在不足。下一步应结合鞍本地区的其他矿床信息进行补充和完善,进一步提升该方法的普适性。