基于PLS模型的高光谱遥感地质岩性成分反演分析

2021-01-07种绍龙

矿山测量 2020年6期

种绍龙

(甘肃省地质矿产勘查开发局测绘勘查院，甘肃兰州 730060)

遥感是通过高空传感器收集地面的电磁波特性，并根据相关特性进行识别、提取和分析地物的技术，已广泛应用于金属矿产勘查与预测、矿化蚀变信息提取、地质构造解释等工程领域中[1-4]。

高光谱遥感诞生于上世纪80年代，具有分辨率高、图谱合一、信息量大等一系列优势，在农业生产、灾害预报、环境监测以及地质探查领域获得了较为成熟的应用，在地质勘查方面，可进行岩性识别、蚀变矿物填图及信息提取、区域性矿产资源调查等多项工作[5-6]。由于高光谱遥感数据具有波段多、信息量大、相关性强、冗余多、信噪比低等特点，在进行计算分析时，需要对数据进行预处理和特征提取，同时尽可能在保证精度的前提下，降低计算维度，实现图像的简易精准化处理[7-8]。在岩性定量识别方面，专家学者提出了主成分分析法、多元线性回归、典型相关分析等单一或者多种判别法组合应用的先例[9-10]，取得了一定的成果，但是上述方法均有各自的局限之处，如多元线性回归法无法体现复杂变量之间的相关性、典型相关分析仅能描述两组变量间的线性相关性、主成分分析在数据量不足时的解释性较差等缺点，因此，有必要对高光谱遥感的降维计算方法进行改进，以解决上述方法的局限性。

本文拟采用偏最小二乘法(PLS)对高光谱遥感数据进行建模和降维计算，通过对研究区玄武岩和花岗岩的成分反演分析，证明该方法的可行性和准确性。

1 研究区概况

研究区位于甘肃酒泉北山-花牛山-柳园镇一带，面积约为500 km2，平均海拔1 800 m，属温带大陆性干旱高原气候，多年平均气温0℃，多年平均降雨量50 mm，受构造作用和岩浆活动影响，区内出露大面积为玄武岩和花岗岩，研究区内成矿地质条件优越，包含金、铜、铅、锌、铁、银、长石等20余种矿产，因此，开展该区域岩性成分的反演分析，有利于成矿带的预测和确定，研究区地理及地质状况如图1所示。

图1 研究区地理及地质概况

2 研究技术路线

研究主要分为两条线路，一是通过机载CASI、SASI、TASI等获得遥感数据，通过影像预处理和数据降维后，再基于PLS建立相应的分析模型，其中影像预处理主要是通过数据重采样、大气校正和几何校正等几个步骤，数据重采样采用双线性内插法，几何校正采用二次多项式法，大气校正采用FLAASH辐射传输校正模型；数据降维处理主要包括主成分分析和最小噪声分离变换法，后者相比前者降低了对噪声的敏感性，提高图像质量，因此文中选取最小噪声分离变换法对数据进行降维处理。二是通过野外岩石样本采集，获取ASD光谱测试数据和主量成分测试数据，经光谱重采样后，与PLS建模分析结果叠加综合分析后得到研究区的岩性成分反演。如图2所示。

图2 研究思路

3 偏最小二乘法回归模型反演分析法

3.1 方法介绍

偏最小二乘法(Partial Least Square，简称PLS)是众多多元统计方法中的一种，可以对多因变量对多自变量的问题进行回归分析，是集多元线性回归、典型相关分析和主成分分析三种方法优缺点为一体的方法，通过建立变量数据组之间的回归模型，实现遥感数据的降维处理，其基本原理为：

假设经ASD岩石样本数量为n个，岩石主量元素成分数为m，岩石样品在不同波段处的数据为p，分别构建标准化变量矩阵X和Y：

X={n1……nm}×m

(1)

Y={n1……nm}×p

(2)

求解矩阵X和Y的第一主成分值t1和u1，要求使两个矩阵的方差和相关度均达到最大值，即有：

Max{Cov(Xt1,Yu1)}=max

(3)

(4)

式中,w1为XTYYTX对应的最大特征值的特征向量；c1为YTXXTY对应最大特征值的特征向量。

求解得到第一主成分t1和u1后，利用偏最小二乘法分别对X～t1和Y～u1进行回归分析，假设这一回归方程的相关性达到理想精度，则可以停止计算，若回归方程的相关性达不到相关精度要求，则进行第二轮主成分的提取和计算，即重复公式(3)～(4)，直至回归方程的精度达到要求为止。若最终提取的主成分个数为m个，偏最小二乘回归模型就会通过实施Yk对所有主成分进行回归分析，然后再反演表达成Yk关于原变量X1,X2,…，XP的相关回归方程，其中,k=1,2，…，p。

3.2 反演分析流程

在利用偏最小二乘法进行岩性成分反演分析时，主要包括以下三个主要步骤：(1)采用原始反射率法或者去连续统法等方法对光谱数据进行预处理，消除各种非目标因素对于光谱特征的影响；(2)选择适宜的主成分个数，既要保证足够的表达精度和预测方程具有较小的残差平方和，又要尽量降低计算难度，缩短计算过程；(3)建立预测方程并进行对应评价分析；(4)根据模型计算结果对影像图进行处理，得到成分含量分布图。

4 反演结果分析

4.1 预处理方法抉择

由于光谱数据受地形、光照等的影响较为明显，因此在进行回归分析之前需要对其进行预处理，分别使用原始反射率法和去连续统法对研究区玄武岩和花岗岩光谱数据进行建模，并选择不同的主成分个数进行分析，得到的误差结果如表1所示。从表中可以看到：误差值随主成分个数的增加而降低，在玄武岩光谱数据下，采用原始反射率法预处理后的误差均大于1，而采用去连续统法处理后，当主成分个数为9个时，其误差已小于1，故玄武岩的光谱数据建议采用去连续统法进行预处理，同时主成分个数为9个即可；在花岗岩光谱数据下，采用原始反射率法和去连续统法的误差均小于1，但去连续统法的误差相对更小，故建议也采用去连续统法对数据进行预处理，且主成分个数仅需2个。综上所述：去连续统法不仅可以减少计算过程，同时还能提升岩性反演精度。

表1 不同处理方法误差分析结果

4.2 反演预测值与实测值情况

选定主成分个数之后，分别对玄武岩和花岗岩进行PLS建模分析，其中，玄武岩的ASD样本个数为26个，花岗岩的ASD样本个数为30个。通过建模分析分别获得研究区玄武岩和花岗岩的氧化物预测值和实测值的相关关系，如图3所示。从图中可以看到：玄武岩和花岗岩岩性成分的实测值和预测值的拟合度较高(拟合度R>0.8)，其中，玄武岩采用PLS模型反演分析的岩性各氧化物成分含量平均拟合度达到0.91，花岗岩采用PLS模型反演分析的岩性各氧化物成分含量平均拟合度达到0.93，表明本文提出的PLS回归模型反演分析法具有较高的计算精度，可在实际地质勘探工程中予以合理应用；玄武岩的Si成分含量最高，约为50%左右，Ti成分含量最低，约为1%左右，花岗岩的Si成分含量最高，达到70%，Mg和Ti成分含量最低，仅为0.5%左右。

图3 反演预测值与实测值关系

4.3 岩性成分分布情况

根据玄武岩和花岗岩各氧化物的模型计算结果，对研究区的影像进行叠加计算，然后去除交叉区域，可获得不同主量矿物元素的分布情况，如图4所示。从图中可以看到：研究区AL元素含量约为15%左右，呈均匀分布、几乎没有富集区域；Ca元素主要集中在研究区西南地区，其富集程度达到20%以上，东部地区的Ca元素含量相对较少；Fe元素的分布情况与Ca元素相似，但含量相对较低，约为10%左右，中南局部地区出现Fe元素的富集情况；Mg元素也主要集中于西南地区，其含量达到10%以上；Si元素与其它元素分布不同，其主要分布于研究区中东部地区，含量达到70%以上，西南地区的Si元素含量约为40%左右；Ti元素主要集中于研究区西南地区，但其含量仅为3%左右。