遥感异常与地球化学指标间的映射关系模型探索

2021-04-08陈三明马云来宋世鹏

矿产与地质 2021年1期

陈三明，马云来，韦龙，董霖，宋世鹏

(桂林理工大学地球科学学院，广西桂林 541006)

0 引言

前人在这方面的工作成果颇为丰富，代表性研究工作有陈涛[1]基于TM数据利用掩膜去除干扰因素、采用主成份分析、密度分割方法结合地质信息提取研究区矿化异常信息；苏一鸣[2]等采用去干扰异常主分量门限化技术，提取羟基、铁染异常；王宏[3]等利用高分一号、ETM+影像数据对区内岩性、构造进行解译，通过主成份分析法和比值法结合ETM+和ASTER影像数据对区内矿化蚀变信息进行提取；唐尧[4]等利用遥感数据与蚀变矿物的谱带特征提取矿化蚀变信息，对ETM的不同波段组合进行增强处理后确定围岩蚀变的范围和强度；汪子义[5]等利用混合调制匹配滤波方法开展基于OLI数据的蚀变矿物组合提取研究，并与Hyperion 高光谱数据提取结果进行对比验证，分析远景区的蚀变矿物组合分带特征；贺金鑫[6]等基于LandSat 8数据，以主成份分析法和比值法结合的方法增强蚀变信息，等等。但由于遥感信息的不确定性和反演多解性，单纯的利用遥感矿化蚀变信息找矿实际效果并不太理想。

笔者认为，如果将宏观层面上的遥感矿化蚀变信息与微观层面具有明确找矿指示意义的地球化学示矿指标在中观层面上建立模型联系起来，通过对比印证，就能够大大减少单一运用遥感信息带来的不确定性和反演的多解性，从而就能增加遥感找矿的可靠性。

1 研究方法及遥感数据预处理

1.1 研究方法

选取LandSat 8 数据、6种元素(Cu、Pb、Zn、Co、Ni和Cr)的化探数据为基本的信息源，先对遥感影像的多光谱波段(Band 1～Band 7)进行预处理，包括辐射校正(包括辐射定标、大气校正)、坐标精校正、高程校正，再和配准后的全色波段(Band 8)融合，然后对融合后的影像利用已知的工作区数据建立MASK掩膜，并选取相应的波段，使用PCA主成份分析法、比值法结合阈值分割法得到增强的蚀变异常信息，初步划定异常区域。最后把得到的增强后的弱化蚀变信息同化探数据分析出的弱化异常信息对比分析，在研究区地质资料及气象水文等资料相结合的基础上，通过建立地球化学指标与LandSat 8的特征波段净反射率之间的统计学映射关系，建立示矿综合地球化学指标与LandSat 8特征波段净反射率组合关系定量模型，以此为基础，对未知区域开展遥感地球化学找矿预测，最终圈定异常区域和预测区域。

具体过程见图1。

图1 技术路线图Fig.1 Technology roadmap

1.2 遥感数据预处理

这里的数据预处理指的是对遥感数据LandSat 8的处理，化探数据的处理会在下面的章节中另外阐述。数据预处理包括遥感原始数据的辐射定标、大气校正、坐标精校正、高程校正和图像融合。辐射定标是将传感器记录的电压或数字量化值转化为辐射亮度(Radiance)或者转化为地表(表观)反射率(Reflectance)和亮度温度(Brightness Temperature)的过程，方法有绝对定标和相对定标两种。本次工作直接利用ENVI 5.3的辐射定标工具(Radiometric Calibration)来完成对原始数据的定标操作。大气校正的目的是消除大气和光照等因素对地物反射的影响，获取地物反射率、辐射、地表温度等真实物理模型参数，用来消除大气中水蒸气、氧气、二氧化碳、甲烷和臭氧等对地物反射的影响，消除大气分子和气溶胶散射的影响。大多数情况下，大气校正同时也是反演地物真实反射率的过程，方法包括基于大气参数的模型校正方法以及基于遥感数据统计特征规律的回归分析方法。坐标精校正是赋予遥感数据像元的空间定位属性，根据研究区的分带本次工作统一使用西安80坐标系。高程校正是为了减少地面高程差对提取结果的影响，在高差地貌其伏较大的山地丘陵地区尤为明显。图像融合是指使用特定的数学算法，将两幅或多幅遥感影像数据，在设定好的地理坐标系中融合合成一幅新的影像的过程。本次工作采用多光谱影像数据(Band 1～Band 7)和全色波段(Band 8)融合，既保留了丰富的光谱信息，同时提高了的空间分辨率，为接下来的信息提取提高了精度。图2为预处理前后遥感影像的对比图。

图2 数据预处理前后遥感图像对比图Fig.2 Comparison of remote sensing images before and after data preprocessing

2 蚀变信息提取

2.1 主成份分析法与比值法相结合的蚀变信息提取

2.1.1 铁染和羟基异常信息提取

研究区内与铁离子有关的矿物主要是赤铁矿、针铁矿等，根据以上矿物的波谱特征，在0.45～0.52 μm和0.875 μm附近的波曲线具有2个吸收谷，其正好对应LandSat 8数据的Band 2和Band 5波段。0.630～0.680 μm和1.500～1.700 μm附近可见较高的反射峰，分别对应LandSat 8数据的Band 4和Band 6波段[6]。因此可以选取Band 2、Band 4、Band 5、Band 6波段进行主成份分析提取铁染异常。但是，经分析后发现这4个波段所提提取的特征向量值并不符合铁离子矿物的波谱特征。根据梁昊[7]、马威[8]等人的资料进一步研究发现，利用比值法即Band 6/Band 5，可以进一步增强铁化蚀变图像的亮度反差。由此可见，选用band 6/Band 5代替Band 6，即用Band 2、Band 4、Band 5、Band 6/Band 5波段组合进行主成份分析，能够进一步增强铁染蚀变信息。经过主成份分析后计算所得特征向量值见表1。

负值代表吸收，正值代表反射，绝对值的大小代表能力的强弱。

由表1可见，PC 1主要反映的是Band 4、5的加信息；PC 2主要反映的是Band 2的加信息和Band 5的减信息；PC 3主要反映的是Band 3的加信息和Band 2的减信息；PC 4则几乎全部反映Band 6/Band 5的加信息。

表1 铁染蚀变异常特征向量Table 1 Feature vector of iron-stained alteration anomaly

根据铁染蚀变异常的提取原理，波段特征信息在LandSat 8遥感影像数据中的Band 2、Band 5波段具有强吸收特征与Band 4或Band 6/Band 5具有较高的反射率，通常还具有Band 2、5的贡献值系数与Band 4的贡献值系数相反，Band 4和Band 6符号相同的特点。因此选择PC 3作为包含铁染蚀变异常的主分量进行分析。

同理，研究区内相关黏土矿物在LandSat 8数据的Band 2 处的反射率呈现较高的反射峰，在Band 5处可见较小的吸收谷，而黏土矿物在Band 6处的反射率较高，在Band 7处黏土矿物均出现强烈的吸收特征。因此，选取Band 2、Band 5、Band 6、Band 7这4个波段进行主成份分析，并提取研究区内的羟基异常。主成份分析后得到的特征向量值见表2。

表2 羟基异常特征向量Table 2 Characteristic vector of hydroxyl anomaly

由表2可见，PC 1主要反映的是Band 5、Band 6、Band 7的的加信息；PC 2主要反映的是Band 2、Band 5的加信息和Band 7的减信息；PC 3主要反映的是Band 2的加信息；PC 4主要反映的是Band 6的减信息和Band 7的加信息。根据含羟基类矿物的波谱特征，在LandSat 8遥感影像数据中的Band 2、Band 6波段强吸收特征，在Band 5、Band 7具有较高的反射率，通常还有Band 5、Band 7的贡献值系数与Band 2的贡献值系数相反，Band 2和Band 6符号相同的特征。最后选定PC 4作为主分量代表羟基异常，PC 4图像中的高亮度信息就是表征含羟基类矿物的蚀变异常信息。

3.1.2 交互式阈值拉伸及低通滤波处理

通过分析特征向量后，利用交互式拉伸工具增强弱化的铁染、羟基蚀变信息，能更直观的观测到异常。通常取统计数据的均值加二至三倍的均方差数值拉伸异常，然而这样得到的遥感异常出现所谓的“焦盐现象”，这不利于异常的空间性规律识别。此时再将异常的数据在ENVI中的数学形态学模块里，进行低通滤波处理，得到模糊影像，但遥感异常的空间性得到了显著的增强，一些碎点异常得到了有效的抑制。

3.2 蚀变信息分级处理

为了进一步直观地观测异常信息的图像，利用Surfer 软件，在对相关主分量图像进行3×3的高斯低通滤波处理和数字拉伸的基础上，结合研究区实际情况，对于不同类型的蚀变信息建立相应的定量化的标准，用来进行阈值分割提取。本文采用主分量分析门限法划分蚀变异常等级，即利用x+Kσ来划分异常等级。x为均值，σ为标准差，K值一般取1～3。根据研究区内特点，异常均分为三级，K值分别为1、1.5、2，即以σ、1.5σ、2σ为阈值划分为一级、二级和三级异常。图3、图4分别为加载底图后的铁染蚀变异常和羟基异常。

图3 加载底图后的铁染蚀变异常分布图Fig.3 Distribution of iron-stained alteration anomaly after loading of the base map

3 化探数据处理与分析

3.1 化探数据的处理和统计分析

3.1.1 成矿元素的聚类分析

在实验室内先对野外样品进行分析测试，然后利用SPSS对地球化学数据进行统计分析，根据分析结果，最终选定Co、Ni、Cr、Cu、Pb、Zn六种元素作为数据样本。通常为了描述数据衡量各组数据之间的近似源关系，通过不同的计算方法将需要研究的数据进行分类，将多个类别数据归结到少数几个大类的不同簇中，一般使用聚类分析的方法。地球化学元素本身在地表漫长的演化过程中，经过各种地质作用以及地球化学作用下形成如今的地球化学环境，其组合具有一定的解释成矿规律的能力。处理时，先把六种元素的地球化学数值做甚高值去除后标准化，然后聚类分析，得到树状图。根据六元素聚类分析树状联接图可以看出，这些元素之间存在的关联性可知大体可分为两大类。其中Pb、Co聚为一类，Cu、Ni、Zn、Cr聚为另一类。聚类分析的结果不是特别的理想，表明我们需要进行进一步的地球化学指标的提取和组合。

3.1.2 成矿元素的因子分析

因子分析是较常用的一种统计学分析方法，一般伴随使用聚类分析作为参考。元素地球化学数据因子分析结果见表3。

从相关性矩阵表3可以求得其特征值和累积百分比，见表4。从主因子的累积百分比看到，前两个特征值所代表的方差已经占总方差的83%。

根据元素统计情况及因子分析情况结合碎石图，可知一共有3个特征值超过1的主成份，这3个成份的贡献率按顺序分别是72.838%，10.470%，7.833%，可以看出最后因子的贡献率较低不足10%，因此舍弃只保留前2个主成份因子。由表4得知，前2个主成份因子的累积贡献率达到了83.308%，即代表这2个主成份因子可以解释表达了原数据的83.308%的主要内容。

表4 总方差解释Table 4 Explanation of total variance

结合因子分析旋转成份矩阵表(表5)分析，两个因子F1、F2已经提取了Cu、Pb、Zn、Co、Ni、Cr，元素的大部分信息。根据旋转成分矩阵，计算得出2个因子公式如下：

表5 旋转后的成份得分系数矩阵Table 5 Component score coefficient matrix after rotation

F1=0.393×Cu-0.316×Pb+0.023×Zn+0.024×Co+0.485×Ni+0.038×Cr

F2=-0.323×Cu+0.573×Pb+0.330×Zn-0.198×Co-0.262×Ni+0.684×Cr

3.2 地球化学因子关系提取

指示元素的组合关系一般包括定性组合和定量(半定量)组合[9]。定量(半定量)的组合关系主要用来表征各个元素间数学特征值及其空间分布规律状况。按组合性质分为如下几种表现形式：

1)简单含量关系：主要根据相关性质，如划分限定的阈值、利用某些指示元素之间的含量特征形式来指示含矿异常和非矿异常。

2)比值：指单元素之间的比值关系，因其在矿床地球化学异常三度空间中随空间部位动态变化形成关系式效应，故比值的关系映射动态变化具有一定的规律性。

3)多元素组合比值：分为垒乘异常及比值和垒加异常及比值，例如：(a×b)/(c×d)和(a+b)/(c+d)。这种组合关系一般是相对的，没有比较严格的组合标准。可以参考元素的活动性顺序来划分组合关系。

4)指示元素之间含量的关系：主要表现在元素对或多元素的组合含量之间的共消长关系和互消长关系。利用指示元素含量间相关性来确定地球化学异常时，将其中某些元素的相关系数用作指标。

分析提取出的两种因子的性质，观察其旋转成份矩阵可以知道：F1主要代表了Cu、Co、Ni、Cr元素，反映研究去的岩石元素共生组合，F1因子中的Cu、Co、Cr、Zn、Pb在岩浆型铜镍硫化物矿床中对成矿有明显指示作用[7]，且Ni本身为良好指示元素。F2因子主要提取了Pb、Zn、Cr元素的信息，代表着成矿需要的重要元素组合。通过多次的组合计算分析，最终选定地球化学指标因子F1(Cu-Co-Ni)与因子F2(Pb-Zn-Cr)，把F1/F2的值作为找矿地球化学敏感指标，这样就建立了遥感地球化学指标间的模型的因变量(图5)。

图5 F1、F2因子异常分布Fig.5 Anomaly distribution of F1 and F2 factors(a) 金川铜镍矿区元素弱地球化学F1(Cu-Co-NI)因子异常分布图 (b) 金川铜镍矿区元素弱地球化学F2(Pb-Zn-Cr)因子异常分布图

3.3 遥感地球化学统计学映射模型的建立

统计分析地球化学元素时，在样本微量元素与光谱反射率统计分析的元素选择中，要先对微量元素的结果进行多元统计，分析其与主成矿元素的相关性；然后筛选出研究区中和主成矿元素成矿关系密切的元素组合，对选出的元素组合进行标准化，提取出最能代表主成矿元素成矿的元素综合异常数据(用x表示)，建立地物光谱反射率与微量元素地球化学指标之间的关系函数。对上述筛选出的样本光谱反射率、元素综合异常数据进行统计分析，得出其相关函数，见公式(1)：

ρ=F(x)

(1)

其中：ρ为样本光谱反射率；F为样本光谱反射率与微量元素地球化学指标之间的统计关系(函数)。可以求出遥感影像数据灰度值与地物光谱反射率之间的方程为

L=ρ(E0·cosθ)/π·d2=DN·gain+bias，

可以求出：

ρ=(DN·gain+bias)·π·d2/E0·cosθ

(2)

其中：L为地物在大气顶部的辐射能量值；ρ为样品的光谱反射率；DN为样本的灰度均值；d为日地天文单位距离；E0为大气顶部的太阳辐照度；θ为成像时的太阳天顶角；gain和bias分别为图像的增益量与偏移量，可从遥感数据的头文件中读取。

由公式(1)、(2)组成方程组，可以推导出DN与X之间的对应关系：

(3)

由此通过F(x)函数，建立了遥感影像灰度值与地球化学异常之间的判别函数关系，也即提取出的遥感地球化学异常模型。最后将遥感数据各波段的灰度值代入式(3)，就可以得到代表综合元素异常值的一组数据，其值的高低代表所求遥感地球化学异常的强弱，x的高值区就是通过遥感数据获取的遥感地球化学异常区。

4 应用研究

4.1 金川矿区测区地质地球化学概况

金川矿床位于甘肃西部龙首山脉东段北坡，是世界第三大铜镍硫化物矿床。大地构造位置上位于华北地台阿拉善地块西南缘的龙首山隆起带内，北依准噶尔晚古生代褶皱系，南邻北祁连褶皱带，东部为华北地台，西隔祁连早古生代褶皱系与柴达木地块相望。

金川矿区南延测区出露从太古代到古生代的一系列地层，老地层的原始沉积建造已被彻底改造。受区域性构造作用影响，各岩性带呈近NW-SE向带状平行展布。测区出现多期强烈的基性—超基性岩浆活动为特征，孕育着超大型铜镍矿床的白家嘴子超基性杂岩体位于普查区北面，辉绿岩体主要呈脉状或岩枝产于各时代的地层中。

测区地球化学统计表明，Cu、Pb、Zn等金属元素在各时代蚀变火山岩中含量最高，在变质碎屑岩及含碳变质岩中较高，在泥质及碳酸盐岩中含量普遍较低。区内Cu、Pb、Zn等金属元素的富集成矿与火山活动及区域变质作用密切相关。印支期侵入岩Pb元素含量普遍较高，为地壳克拉克值的2倍以上；海西晚期钾长花岗岩Pb含量均较高，为地壳克拉克值的2倍；二长花岗岩Pb、Zn、As、Ni、Co、Mo含量均较高，超过了地壳克拉克值。海西中期各类岩体的Pb元素含量普遍偏高，其他元素含量则相对偏低。超基性岩Ni、Co元素含量是克拉克值的200倍以上，其他元素含量则较低。加里东期斜长花岗岩中的Pb、As含量普遍大于地壳克拉克值，Cu、Zn只在局部地段含量较高：石英闪长岩As、Ni、Co含量都高于地壳克拉克值数倍以上。不同地区不同期次的中性、酸性侵入岩中，Cu、Pb、Zn元素含量普遍较高，并呈正消长关系[10]。

4.2 遥感地球化学异常指标的提取结果分析及异常评价

根据LandSat 8遥感地球化学异常综合模型进行信息提取的空间展布特征，依据异常大小进行密度分割，本区可划分为2个成矿远景预测区，分别编号为1、2号(图6)。第1个预测大区(一级预测区)之中分2个小预测区，分别标记为1-1、1-2、，第2个预测大区(二级预测区)只有1个别标记为2-1。其中，第1预测区1-1为已开采的矿床，从预测的结果来看说明本文建立的基于灰色关联度的高光谱遥感地球化学综合预测模型达到了较为理想的结果。1-2预测区则与NW向的基性超基性岩脉的空间分布较为一致，且形成NW向的条状轴向分布格局。地层岩性主要为上元古界震旦系烧火筒群青灰色微晶白云岩、含砾绢云千枚岩等。构造以断裂为主，F2、F3和F4断裂总体走向NW，异常位于3组断裂围限的震旦系中。岩浆岩主要为基性的辉绿岩脉，少量超基性辉橄岩脉，特点是岩脉出露面积较大，形态与走向与全区一致，辉绿岩脉一般长1500～1800 m，宽30～50 m，其余小岩脉呈零星分布。第2预测区2-1为团状分布格局，似近SN向展布，与矿区的整体构造走向近乎垂直，需要进一步综合其他勘查手段作进一步的分析研究。