基于GF-2影像和随机森林算法的花岗伟晶岩提取

杜晓川, 娄德波, 徐林刚, 范莹琳, 张琳, 李婉悦

(1.中国地质大学(北京)地球科学与资源学院,北京 100083; 2.中国地质科学院矿产资源研究所,自然资源部成矿作用与资源评价重点实验室,北京 100037)

0 引言

随着近年来新能源、新兴材料的发展,锂等关键金属越来越受到重视[1],金属锂由于其特殊的物理化学性质,被广泛应用到诸如电池、玻璃、陶瓷等各个工业领域,具有极高的战略价值,其不仅是高科技产业必需材料,也逐渐成为日常生活中广泛使用的常规材料。当前我国国内锂业发展高度依赖矿物型锂矿特别是花岗伟晶岩型锂矿的供应[2],此类型锂矿床通常位于大型花岗岩侵入体边缘,受地质构造控制明显,通常以伟晶岩脉集群的形式产出,总体呈带状分布,单条伟晶岩脉形态多为脉状或板状[3],因此基于高分辨率影像的遥感技术是寻找和发现伟晶岩脉的有效方法[4]。

遥感领域相关技术不断发展成熟,使得遥感技术在地质信息提取的研究中得到了广泛应用,有许多利用遥感技术进行含锂伟晶岩探测的实例。如代晶晶等[5-6]在川西甲基卡地区进行了锂辉石等伟晶岩稀有金属矿床相关矿石矿物的波谱测试研究和伟晶岩遥感地质填图; Cardoso-Fernandes等[7]利用多源遥感影像数据,采用RGB彩色合成、波段比值分析、主成分分析等方法识别含锂伟晶岩的相关异常,完成蚀变填图; 姚佛军等[4]通过增强ASTER遥感图像中的差异来增强地质体的岩性差异,并突出含锂伟晶岩的遥感特征; 王海宇[8]以WorldView-3为数据源,并基于全卷积网络(fully convolutional networks, FCN)网络和U-net网络2种深度学习语义分割模型识别伟晶岩脉信息。前人通过遥感影像数据在含锂伟晶岩提取方面开展了大量的理论研究和实践,但大多数研究者侧重于使用低空间分辨影像数据,提取方法也较为传统,尽管已有基于高空间分辨率影像数据与深度学习网络模型提取伟晶岩岩脉的先例,但该方法耗时较长,效率较低,且对样本数量有一定要求。因此本文在众多机器学习算法中选用随机森林算法,其算法具有如下优点: ①数据兼容性强,数据集无须规范化,且可以直接处理高维数据; ②抗噪声能力强,受异常值影响较小; ③鲁棒性强,不容易发生过拟合现象; ④计算效率高、速度快、泛化能力强[9-11],且该方法结合高空间分辨率影像数据在提取复杂地物的效果较好。

据此,以青海省海西蒙古族藏族自治州天峻县扎卡东南部区域为研究区,提取预处理后的GF-2高空间分辨率影像数据中各类型地物的光谱特征、纹理特征、指数特征、地形特征及边缘特征等25个特征变量构建特征子集,对子集中的特征变量进行特征重要性分析,依据特征重要性得分进行特征优选,确定最优的提取花岗伟晶岩的特征组合,并进行随机森林分类,对分类结果进行精度评定,为后续该地区的含锂花岗伟晶岩找矿勘查工作提供依据。

1 研究区概况与数据源

1.1 研究区概况

研究区位于青海省海西蒙古族藏族自治州天峻县扎卡东南部,地理坐标为E99°27′16″～99°30′32″, N36°53′36″～36°54′43″(图1)。天峻县地处青海省东北部,青海湖西侧、祁连山南麓,柴达木盆地东部。境内高山纵横,山脉呈东南西北走向,地形以山地为主,高山、中低山、山谷和山间盆地相间分布。平均海拔在4 000 m以上,地处高原寒带气候,年平均气温-1.5 ℃,年降水量360 mm,干旱少雨,多风沙。

图1 研究区遥感影像

1.2 数据源及预处理

GF-2卫星于2014年8月19日发射,是我国自主研制的首颗空间分辨率优于1 m的民用光学遥感卫星。搭载有2台高分辨率的1 m全色与4 m多光谱相机,传感器的波段参数如表1所示。对大气校正后的多光谱影像以及定标后的全色影像做正射校正处理及融合处理,获得空间分辨率为1 m的彩色合成影像,再通过研究区域的矢量范围对影像数据进行裁剪,得到研究区遥感影像数据。

表1 GF-2全色多光谱相机波段参数

2 研究方法

2.1 特征提取及样本选择

花岗伟晶岩(图2(a))在分布位置及产出形态上受构造控制作用明显,研究区内断裂沿NW-SE向分布,因此花岗伟晶岩大多沿NW-SE向展布,呈条状、脉状产出; 在影像色调上,其多呈灰白色或略带土黄色; 由于研究区内海拔较高、风沙较多、风化能力较强,且花岗伟晶岩晶体较大,具有较强的抗风化性,因而在高度上要突出于周围地物,并且其表面纹理粗糙,边界较明显。研究区内其他地物类型主要分为道路(图2(b))、干涸河流(图2(c))、冰雪(图2(d))、阳坡裸地(图2(e))、阴坡裸地(图2(f))5类,其中阳坡裸地与阴坡裸地以研究区内最高海拔线为界线(图1),海拔沿界线南北两侧逐渐降低,加上卫星拍摄角度影响,使得二者地物的色调差别较大。阳坡裸地坡向朝南,色调以土黄色、棕色为主; 阴坡裸地坡向朝北,色调以棕黑色为主; 干涸河流的形状与边界较为明显,宽度较窄; 冰雪主要分布在雪线较低的阴坡裸地上,色调以亮白色为主,形状为长条状,展布大多为南北向; 道路色调为墨青色、边界明显。基于不同地物的特点,本文选择提取地物的光谱特征、纹理特征、指数特征、地形特征、边缘特征作为随机森林分类所需的特征变量。光谱特征包括GF-2高空间分辨率影像数据4个波段的灰度特征(B,G,R,NIR)、影像灰度平均值特征(Img Mean)、影像灰度标准差特征(Img Std),以及对影像原始4个波段进行限制对比度自适应直方图均衡化的CLAHE特征(CLAHE B,CLAHE G,CLAHE R,CLAHE NIR)。

(a) 花岗伟晶岩样本 (b) 道路样本 (c) 干涸河流样本 (d) 冰雪样本 (e) 阳坡裸地样本 (f) 阴坡裸地样本

纹理特征是一种全局特征,对噪声具有较强抵御能力,能刻画出图像区域所对应地物表面的性质,具有旋转不变的特性[12],且已有研究表明纹理特征有助于改善地表信息提取的精度[13-16]。本文在提取纹理特征前先对影像数据进行主成分分析变换,该方法通过应用正交变换减少相关数据的冗余,减少各组分间的相关性,尽可能将有用信息集中到新的组分中,且各新组分间互不相关,达到信息增强、降维及减少计算量的目的[17]。其中影像数据的第一主成分(PC1)包含所有波段中98.5%的特征信息,因此将PC1作为影像的光谱特征的同时,仅对PC1采用灰度共生矩阵方法提取纹理特征,选取均值(Mean)、方差(Variance)、均质性(Homogeneity)、对比度(Contrast)、相异性(Dissimilarity)、熵(Entropy)、相关性(Correlation)、角二阶矩(angular second moment,ASM)等8种彼此相关性弱的纹理信息特征统计量作为纹理特征。

指数特征是通过一定公式的波段运算所获取的灰度影像。指数特征可以消除地形差异影响,突出地物特征[18]。基于阴坡裸地存在阴影的特点,本文提取用于突出阴影区域与明亮区域对比的阴影指数(shadow index,SI)特征,其计算公式为:

SI=(B+G+R+NIR)/4

,

(1)

式中B,G,R,NIR分别为蓝光、绿光、红光及近红外波段的反射率。基于研究区内裸地范围较大,本文提取突出裸地区域亮度的裸地区域指数(bare area index,BAI)特征,其计算公式为:

BAI=(B-NIR)/(B+NIR)

。

(2)

地形特征是描述地势变化的特征,由于研究区内地势起伏较大,且地物与地形因素有一定关联,因此有必要提取坡度(Slope)、坡向(Aspect)、高程等地形特征。本文从美国地质调查局获取研究区30 m空间分辨率的ASTER GDEM数字高程模型(digital elevation model,DEM)数据,并分别提取坡度、坡向,DEM地形特征。

边缘特征是影像中特性(如像素灰度、纹理等)分布不连续处,影像周围特性有阶跃变化的像素集合,其中研究区内花岗伟晶岩、干涸河流、冰雪、道路的边缘特征较明显,边缘灰度梯度变化较大。Canny边缘检测能较精确估算出图像边缘的强度、梯度方向,具有定位准确、单边响应和信噪比高等优势[19]。将从GF-2高空间分辨率影像中提取出的光谱特征、纹理特征、指数特征、地形特征与边缘特征等25个特征变量构成初始特征空间。

各类型地物选取25个影像样本构成训练样本集,图2为32像素×32像素的影像样本块展示。为提高分类检测精确度,选取原则依据样本要尽可能包含各个类别地物的各种不同类型特征,且单个样本尽量仅包含单个类别地物。根据随机森林实验需求,将样本集按照2∶8的比例随机划分为测试集和训练集。

2.2 随机森林算法及参数组合优化

随机森林算法最早是由Breiman[20]提出的一种机器学习算法,本文依据Gini系数最小原则构建分类和决策树(classification and regression tree,CART),并基于多棵决策树对样本进行训练,根据训练得到的模型,对未知待测样本类别进行预测的一种监督学习分类算法[21]。

本文使用Python语言编程调用sklearn库进行随机森林模型训练,采用随机搜索交叉验证(RandomizedSearchCV)与网格搜索交叉验证(GridSearchCV)[22]的双重方式进行参数优化。单一的GridSearchCV是对每一组排列组合的参数进行遍历,最后选取出最好的一组参数,但由于参数排列组合数量过多,逐个遍历面临相当耗时与维度灾难等问题,因此在进行GridSearchCV前先通过RandomizedSearchCV(设置随机搭配超参数组合次数(n_iter)为200次,交叉验证折数(cv)为5次)对取值范围内的参数进行有限次的随机排列组合,并对其进行遍历,从中选取出一组较优参数组合,再通过进行10折的GridSearchCV遍历选取出最优参数组合。本文最优参数组合设置如表2所示。

表2 最优参数组合

2.3 特征重要性评估与特征优选

特征集中包含着不同属性、不同尺度的多样化的众多特征变量,过多的特征变量会造成数据的冗余、分类时间过长、以及分类精度下降等问题[23],因此需要通过特征选择筛选出对模型结果影响较大的特征变量构成新的特征子集。特征选择主要是通过一定的准则对特征重要性排序,或者在分类精度保证的情况下选择一组数量最小,且最具有代表性的特征子集[24]。本文先通过调用sklearn库中基于Gini 系数原理进行特征重要性计算的feature_importances参数对每个特征进行重要性评估,对得到的特征重要性得分按降序排列(图3),然后采用顺序前向特征选择法[25],在每一次迭代中选择当前特征集中特征重要性得分最高的特征变量进行模型训练,并获取相应模型的测试集精度,依次增加特征变量数量直至全部纳入,从中选取出测试集精度最高、特征变量数量相对最少的特征组合。

图3 特征重要性排序

2.4 形态学运算

基于最佳参数组合与最优特征变量组合的随机森林模型对研究区各地物类型进行分类,对分类结果进行先膨胀后腐蚀的形态学闭运算操作,闭运算在填充同类型地物间的细微裂隙,使地物间更连接且更平滑的同时,还可以消除孤立噪声点的影响。

3 实验结果与精度分析

3.1 特征优选评价

分别对特征空间中的25个特征变量通过顺序向前特征选择法逐一训练后的结果进行评估,生成特征变量个数与测试集精度关系图(图4),由图4可知: ①由首个特征变量逐个增加至16个特征变量参与模型训练期间,测试集精度整体呈现迅速上升至平稳上升趋势,由最初的84.01%上升至96.24%; ②特征变量数量在17～25个期间,由于特征变量数量过多,产生数据冗余、训练耗时增加等问题,导致分类器出现一定的过拟合现象,造成测试集精度有些许下降,由96.24%下降至95.94%,由此可知并不是特征变量数量越多测试集精度就越高。根据测试集精度结果,选取前16个特征变量组合作为最优特征集,并基于该特征集进行随机森林分类。

图4 特征变量个数与测试集精度关系

3.2 特征变量分析

遥感影像中不同波段对于分类不同地物具有重要意义[24],波段间反射率值的差异有助于区分不同类型地物,通过计算不同特征变量在不同类型地物的反射率值,并对反射率值进行归一化处理,得到各类型地物波谱曲线(图5)。光谱特征(图5(a))中的NIR波段的地物波谱曲线反射率值差异较大,可以很好地区分各类型地物,并且NIR波段的特征重要性得分最高; 纹理特征(图5(b))中Mean波段的各类型地物反射率值有所差异,具有区分地物的能力,花岗伟晶岩的纹理特征在大多数波段反射率值都较高,原因为风化剥蚀等作用使得花岗伟晶岩岩石表面产生较多纹理; 指数特征(图5(c))中阴坡裸地的阴影指数反射率值较低,与其他类型地物有明显区分; 地形特征(图5(d))中,由于研究区内高程起伏较大,高程范围由3 671～4 029 m,高差相差358 m,使得DEM波段中不同类型地物反射率值差异较大。

(a) 光谱特征波谱曲线 (b) 纹理特征波谱曲线

3.3 分类结果与精度分析

3.3.1 随机森林分类结果与精度分析

基于随机森林算法对研究区的GF-2高空间分辨率影像数据进行花岗伟晶岩提取,其各类型地物分类结果如图6所示。此次实验结果共提取出84处花岗伟晶岩,其形状大小不一,以脉状产出为主,成群出现,且轮廓较为明显,与其他类型地物区分较好。在研究区内随机均匀生成1 000个地面验证点,其中花岗伟晶岩验证点104个,道路验证点19个,冰雪验证点51个,阳坡裸地验证点414个,干涸河流验证点87个,阴坡裸地验证点325个,通过计算地面验证点与真实地物类型之间的混淆矩阵获得分类结果的精度统计(表3)评价,其中总体精度达到93.1%; Kappa系数达到0.902; 花岗伟晶岩、冰雪、阳坡裸地及阴坡裸地共4种地物的用户精度和生产者精度都高于93%、道路高于81%、干涸河流高于73%。花岗伟晶岩由于其自身的色调、纹理,以及分布地形海拔都较高等优势,在分类结果上效果较好,用户精度高达94.24%、生产者精度高达98%,以上结果说明基于多特征的随机森林算法在提取花岗伟晶岩信息中有较好的适用性。

表3 分类结果的精度统计

图6 研究区地物分类结果

3.3.2 CLAHE特征层讨论

本文新引入的4个CLAHE特征层中,CLAHE B与CLAHE R作为特征变量参与了随机森林分类,为验证CLAHE B和CLAHE R对分类结果产生的影响,故对优选的16个特征变量中去除CLAHE B和CLAHE R这2个特征变量,只保留余下的14个特征变量参与随机森林分类,并通过计算验证相同的1 000个地面验证点与真实地物类型之间的混淆矩阵,并获得分类结果的精度(表4)评价,其中总体精度为90.4%,较最佳分类结果的总体精度下降2.7百分点; Kappa系数为0.864,较最佳分类结果的Kappa系数下降0.038; 各地物类型的用户精度与生产者精度都有不同程度的下降。

表4 14个特征的分类结果精度统计

CLAHE方法的优越性在于可以通过改变不同地物间的对比度,以突出色调上的差异,相比普通的直方图均衡化方法只增强全局影像对比度,CLAHE方法在此基础上引入了自适应性与限制性,自适应性改进了增强全局对比度中只考虑整体影像区域,而忽略局部影像区域的情况,进而注重影像的局部像素差异,在突出局部影像对比度的同时获得更多细节特征; 限制性是在自适应性基础上的进一步改进,目的是为了避免均衡化过程中产生的色调不连续与对比度过度增强等问题。以精度下降最为明显的干涸河流为例,且为了效果更加直观,故采用真彩色的方式展示(图7),CLAHE真彩色样本影像中(图7(b)),干涸河流的亮度要明显亮于阳坡裸地,且对比度更加突出,边缘轮廓更加清晰,2类地物在色调上有明显差异,有利于分类精度的提升,证实了CLAHE方法在区分色调较为相近地物时的有效性,同时CLAHE特征增强了花岗伟晶岩与其他类型地物间亮度与对比度的差异,有助于提升花岗伟晶岩的提取效果,使其用户精度和生产者精度分别上升0.97百分点与1百分点。

(a) 真彩色影像 (b) CLAHE真彩色影像

4 结论

1)本文基于GF-2高空间分辨率影像数据结合特征优选的随机森林算法提出花岗伟晶岩的提取方法。地物分类结果总体精度为93.10%,Kappa系数为0.902; 花岗伟晶岩的用户精度及生产者精度分别为94.24%和98.00%,证实了该方法的有效性。

2)特征变量个数持续增多会产生数据冗余及计算量过大等问题,导致过拟合现象发生。此次研究在特征优化中确定特征变量个数为16个时分类精度最高。

3)本文新引入的限制对比度自适应直方图均衡化(CLAHE)特征层中的CLAHE B和CLAHE R的特征重要性得分较高,并参与随机森林分类。CLAHE特征变量具有提升影像对比度、突出不同类型地物间色调差异等优势,有助于区分色调相近的地物,有利于地物分类精度的提升。