结合多尺度分割和随机森林的变质矿物提取

唐淑兰，孟 勇，王国强，卜 涛

1) 西安财经大学管理学院，西安 710100 2) 中国地质调查局西安地质调查中心，西安 710054

变质矿物的识别是变质岩研究的基础，造山带中的变质岩是造山带不同演化阶段动力学过程的直接记录，其中，中−低压变质带的分布是研究造山过程和地壳演化等地球科学问题的重要参考依据[1].

遥感技术作为对传统实地调查的补充，已广泛应用于矿产勘查和岩性识别中[2−3].一些矿物提取方法得到了快速发展：如比值法、主成分分析（Principal component analysis，PCA）、支持向量机（Support vector machine，SVM）、光谱角（Spectral angle mapping，SAM）、混合像元分解、随机森林(Random forest，RF)等[4−8].同时，图像处理技术也被结合到遥感矿物提取中：如多尺度分割、小波变换等[9−10].具有较高分辨率的ASTER传感器通过14个波段的波谱数据，提供了更加精细的矿物及岩石信息.

矿物信息在遥感影像上表现为弱信息，利用反射带和吸收带的比值，可增强各地质信息之间的光谱差异.可根据各矿物的特征性光谱特征，选择合适的ASTER波段进行比值增强.基于高分辨率遥感影像的光谱特征提取矿物信息，会出现“椒盐现象”[11]；而遥感影像上变质矿物呈现块状或条带状分布，可使用面向对象思想，通过多尺度分割得到相似的对象块，利用影像的光谱、形状、大小、色调、纹理等特征[12]，提高矿物提取精度.为了使遥感影像矿物提取更加精确，各种非参数监督分类法得到了发展，应用较为广泛的包括神经网络（Neural network, NN）、SVM和RF[13−14].RF 是一个集成分类器，将复杂的分类过程分解为一系列决策过程（树）.由于计算速度快、参数要求少、对训练数据的统计假设少、对噪声或过拟合的敏感性较低，RF分类器在遥感信息提取中得到了越来越多的关注[15].但是，由于部分矿物的光谱相似性、变质矿物与主要造岩矿物的混合性，仅基于光谱特征提取变质矿物精度有限.

本文以ASTER数据为研究对象，利用波段比值增强各变质矿物的光谱特征；接着，结合光谱特征和变差函数纹理进行多尺度分割；然后，利用RF提取变质矿物；最后，通过野外调查验证精度.本研究旨在提高变质矿物提取精度，为变质带的遥感解译提供技术参考.

1 地质背景

研究区（图1）位于甘肃省玉门市地区，大地构造隶属于北山造山带.出露的主要地层有：晚太古界、长城系、蓟县系、石炭系、二叠系、侏罗系、白垩系和第四系.发育有石炭纪和二叠纪中酸性侵入岩.其中，晚太古界主要岩性组合为：（云母）石英片岩、（石榴石、黑云母、白云母）斜长片麻岩、（石榴石、黑云母）斜长角闪片麻岩、二长片麻岩、斜长变粒岩、钠长绿帘绿泥片岩和大理岩等；长城系主要岩性为：大理岩、石英岩、长石石英岩和绿泥石石英千枚岩等；蓟县系主要为大理岩，局部见蛇纹石、透辉石和阳起石化；石炭系和二叠系主要为浅变质碎屑岩，局部可见绿泥石化，石炭系碎屑岩中夹少量灰岩；侏罗系和白垩系主要为未变质的陆源碎屑岩沉积.

图1 研究区地质简图Fig.1 Geological sketch of the study area

2 数据与方法

2.1 数据

2.1.1 ASTER 数据

ASTER具有14个波段的波谱数据，近红外及可见光3个波段，可提取铁及稀土元素；短波红外6个波段，可提取含羟基及碳酸盐化蚀变信息；热红外5个波段，可识别石榴石、黑云母及长石等.本文选择甘肃中盐池地区ASTER数据进行变质矿物提取，图像为拍摄于2003-8-15日的L1T级数据，图像清晰，无积雪、植被覆盖，有云覆盖.数据预处理包括串扰校正、辐射定标、大气定标及去云处理，大气定标采用FLAASH模块进行，去云处理是鉴于近红外波段的异常高值即为云覆盖，进行掩膜运算.

2.1.2 变质矿物光谱特征

研究区需提取的标志性矿物为：黑云母（Bi）、白云母（Mus）、角闪石（Am）、绿泥石（Chl）、石榴石（Gt）、阳起石（Act）.根据野外岩石样品光谱测试数据（图2），获得各变质矿物的吸收谱带与ASTER波段的对应关系（表1）.

表1 矿物的吸收谱带与 ASTER 波段的对应关系Table 1 Correspondence relation between the absorption bands of the minerals and ASTER bands

图2 标志性矿物反射率曲线Fig.2 Reflectance curve of marker minerals

2.2 方法

因矿物具有丛集特征，本文基于面向对象思想，采用多尺度分割和RF分类法提取变质矿物.对预处理之后的影像通过比值运算，得到各矿物增强影像；选取矿物特征，构造分类特征向量；利用RF筛选特征并提取矿物；野外采样、薄片鉴定，通过鉴定结果对矿物提取结果进行精度评价.技术流程如图3所示.

图3 技术流程Fig.3 Technical process

2.2.1 比值运算

黑云母在 8.6 μm 处（波段 11）强吸收，8.0 μm（波段 10）和9.0 μm 处（波段 12）强反射.白云母在波段1和波段6强吸收，波段5和波段7强反射.角闪石在波段7、波段9处强反射，而在波段8强吸收.绿泥石在波段1和波段8处有较强吸收谷，在波段5和波段9有强反射.石榴石在热红外9.25 μm 处（波段 12）有强吸收，10.21 μm 处（波段13）有强反射.阳起石在 2.3 μm 处（波段 8）有吸收 谷，2.2 μm (波 段 6)、2.35 μm (波段9)处有反射峰.增强各变质矿物信息的波段比值公式见表2，其中bi表示ASTER影像第i波段的光谱反射率.

表2 各矿物比值公式Table 2 Ratio formula of minerals

2.2.2 多尺度分割

比值法增强的影像有很多伪信息，需进一步处理.因矿物信息多呈块状和条带状，面向对象的分析法可利用高分辨率影像相邻像素之间的关系[16].面向对象特征提取的基础是影像分割，其目的是将影像分割成若干均匀的区域，分割性能至关重要.分割算法需指定尺度参数来控制对象的大小，多尺度策略较单尺度参数更能满足分割的需要.较为有效的方法是从多尺度分割结果中选择全局最优尺度，然后改进欠分割和过分割区域，最后结合最优分割结果与改进后的分割区域.本文先给定一个尺度参数进行过分割，然后使用区域合并细化分割结果，合并的约束条件为变差函数.比值影像为单波段影像.进行初分割时区域合并的依据为光谱相似性，D(X，Y)为两个相邻区域X和Y的平均灰度的欧式距离式（1），n、m分别为X、Y区域的像素总数.

计算每一个初分割对象块的变差函数，得到变差函数矩阵.变差函数可描述影像像素的空间相关性和变异性，提取更详细纹理信息，有3个参数：方向θ、步长d、窗口大小α[17]，可表示为：

式中，γ(θ,d)为像元在方向θ，步长为d的变差函数；N(θ,d)为影像在θ方向的步长为d的像素对的数量，θ取 0°、45°、90°、135°4个方向；F(yi)和F(yi+d)分别为yi和yi+d的灰度值.像素灰度值具有空间自相关性，自相关性随距离增大而变小，所以变差函数的步长d不应超过窗口大小α的一半.本文各对象块为不规则形状，不使用传统的a×a移动窗口，而是以各对象块覆盖范围作为各自的计算窗口.步长d的上限设为不规则对象块各方向像素数（α）的一半，即去掉小数部分保留整数部分（用|α|表示）.因不规则对象块的各方向变差函数向量维数不同，取各方向不同步长变差函数的累计平均值：

将各方向变差函数累计平均值顺序连接得到最终纹理特征向量：

设两个相邻对象块纹理向量分别为γa和γb，通过卡方距离式（5）来评价两块的相似度，当距离小于设定值时合并这两块.其中，k为变差函数向量的维度，这里k=4.

通过分割指数（Global scare，GS）选择最优多尺度分割结果，GS 指数是衡量对象块内部同质性和块间异质性的参数，该值越小说明分割效果越好.

2.2.3 随机森林 (RF)

RF通过集成多个弱分类器（树），采用平均或投票法得到最终分类结果，精度和泛化能力较高，擅于处理高维数据.具体过程：每颗树采用有放回的抽样（Bagging）随机选择原始样本集的子样本集，利用分 类和回 归 树 （Classification and regression trees，CART)算法训练二叉决策树，构建弱分类器，对各分类器的分类结果采用多数投票法输出结果.原始样本集的2/3用于分类，1/3用于验证.RF有2个重要参数：分类的特征数量和决策树的数量.没有被Bagging 采用的数据称为袋外 (Out of bag，OOB)数据，利用OOB预测结果平均错误率来表征不同特征的重要性.RF特征数量一般选输入变量总数的平方根[18].决策树数量的上限一般设置为1000.

根据已有地质资料，以多尺度分割之后的各对象块为单位进行特征选择，共选择特征3类，9个维度：

（1）光谱特征：各矿物在ASTER影像的光谱范围具有特征性光谱特征，选择各对象块的平均灰度（Sp）及标准差（De）作为光谱特征.

（2）纹理特征：选择各对象块的变差函数纹理4个（Va1、Va2、Va3、Va4分别代表 0°、45°、90°、135°变差函数纹理）.

（3）几何特征：选择面积特征（Ar）、形状指数（Sh）、走向特征（Tr）.Ar是指对象块的像素总数(T).对象边缘的平滑性指标即 Sh：[19]，其中，s为形状指数，l为对象块的边缘长度.矿物分布具有丛集性，走向特征能标志矿物的分布情况.走向特征选多个特征向量中，特征值较大的向量的方向.

利用RF筛选特征，首先使用OOB检测某一棵决策树（L）的正确分类个数L（a），然后对OOB中特征a进行扰动，再次检测L的正确分类个数L（a^），|L（a^）−L（a）|为a对L的重要性，最后计算a对所有决策树的重要性的平均值，即为a对整个森林的重要性[20].根据重要性排序，降低特征空间维度.

2.2.4 精度评价

进行两次精度评价，一次通过OOB计算精度，评价算法稳定性、并筛选树数；第二次通过野外实地调查、采样、薄片鉴定得到矿物鉴定结果，用混淆矩阵评价矿物提取精度，分别统计制图精度（Production accuracy，PA）、用户精度（User accuracy，UA）、总精度（Overall accuracy，OA）和Kappa系数.

3 结果与分析

3.1 图像增强

经过比值运算，目标矿物增强的同时其他光谱相似矿物的信息也增强了（图4）.(b12+b10)/b11同时增强了黑云母信息和中基性斜长石矿物信息；(b5+b7)/b6增强了白云母及高岭土等黏土矿物信息，对伊利石、蒙脱石等亦有增强；(b6+b9)/(b8+b7)增强了黑云母和角闪石信息，含绿泥石的变质岩也呈现高值特征；(b1+b9)/b8增强了绿泥石信息，对碳酸盐岩矿物也有较好的区分作用；b13/b12增强石榴石信息的同时，突出了碱性长石族矿物；(b6+b9)/b8可增强阳起石信息，也能突出黑云母、角闪石信息.另外，由于不同岩性所含矿物成分百分比不同，导致比值运算结果有较多不确定信息.但是，增强影像依然可为后续工作提供数据基础.

图4 比值增强结果.（a）Bi；（b）Mus；（c）Am；（d）Chl；（e）Gt；（f）ActFig.4 Ratio enhancement results: (a) Bi; (b) Mus; (c) Am; (d) Chl;(e) Gt; (f) Act

3.2 多尺度分割

本文多尺度分割流程（图5）如下：

图5 多尺度分割流程Fig.5 Multiscale segmentation process

（1）对目标矿物增强之后的影像进行单尺度分割.从左上角到右下角遍历影像，相邻像素灰度距离小于0.002时合并两区域，循环遍历，直到达到尺度限制.分割尺度参数设为0.2，紧致度参数设为0.5，形状参数设为0.1.

（2）遍历初分割之后的影像，计算各对象块的变差函数纹理，构造整体影像的变差函数矩阵.

（3）遍历各对象块，计算当前对象块和相邻对象块的纹理相似性，当两个纹理向量的卡方距离小于0.0003时，合并两对象块，修改矩阵标识.循环执行，直到选定最优分割结果.

（4）计算每次合并之后的整体影像的GS值，记录GS值的变化趋势，输出GS值最小时的分割结果.

随着多尺度分割过程中合并次数增多，各分割结果的GS值均呈现出最小值（图6）.黑云母、白云母、角闪石、绿泥石、石榴石、阳起石分别经过1次、5次、6次、4次、1次、2次合并，GS值最小.GS值最小时的分割结果即为最优分割结果（图7）.增强黑云母和石榴石信息时采用了ASTER热红外波段影像，热红外波段分辨率较低，影像分割收敛更快[21].

图6 GS 值变化趋势Fig.6 Change trend of GS values

图7 多尺度分割结果.（a）Bi；（b）Mus；（c）Am；（d）Chl；（e）Gt；（f）ActFig.7 Multiscale segmentation results: (a) Bi; (b) Mus; (c) Am; (d) Chl;(e) Gt; (f) Act

3.3 RF 提取矿物

3.3.1 样本选择

根据野外调查和目视解译结果，选择目标矿物对象作为样本集.6种矿物样本集如表3，随机选择每一样本集的2/3样本数进行分类，1/3进行精度验证.因角闪石分布较为分散，覆盖面积大，选择的样本最多；而石榴石和阳起石只有少量出现，选择的样本较少.

表3 矿物样本数Table 3 Mineral samples

3.3.2 特征选择

通过OOB检测误差率，进行特征的重要程度归一化排名（图8）.由于比值运算增强了各矿物的光谱特征，光谱特征的重要程度最高；变差函数统计了地质体空间相关性，研究区在45°方向的区域化变量的相关性最强；岩石由多个矿物组分，岩性分布多为条带状或块状分布，面积大小较大程度影响分类结果.随着特征维数的增加，计算效率降低.树数不同时，按照重要程度排名选择不同特征个数进行分类，各结果表明随着特征个数增加提取精度先增大后变小，综合看3个特征量时矿物提取精度最高，如表4所示.表4显示了树数为100时不同特征个数的提取精度，相比于采用9个特征，采用3个特征的各矿物提取精度分别提高了4.3543%、2.4350%、4.6721%、3.8256%、4.9135%及1.6645%.经综合判断，最终选择Sp、Va2、及Ar 3个特征量提取矿物.

图8 特征重要程度排名Fig.8 Ranking of feature importance

表4 不同特征数量矿物提取精度Table 4 Extraction precision of minerals with different characteristic numbers

3.3.3 决策树个数选择

为每个训练集构建验证树（从50到500，步长为50）.通过OOB进行精度验证，各矿物提取精度最高时，对应的树数都不同（图9）.黑云母、白云母、角闪石、绿泥石、石榴石、阳起石精度最高为0.8404、0.7914、0.7634、0.6836、0.6791、0.6165，精度最高时对应的数的棵树分别为350、300、50、200、100、50.树数对矿物提取精度的影响在可控范围之内，从50棵到500棵引起的精度变化最大的是石榴石，精度变化幅度为0.0471，最小的是阳起石，变化幅度为0.0177.

图9 决策树个数选择Fig.9 Number selection of decision trees

3.3.4 矿物提取结果

变质矿物提取结果（图10）主要分布在晚太古代变质地层中，与实际地质情况吻合.受多期次区域构造和热液活动的改造，该区域是变质矿物集中分布区，也是矿化蚀变集中分布区域，发育绿帘石化、绿泥石化、阳起石、碳酸盐化、矽卡岩化及绢云母化.结果显示，黑云母主要分布在西北部和西南部的晚太古代变质地层中，相关岩性有：黑云斜长片麻岩、二云斜长片麻岩、二云二长片麻岩、斜长角闪岩、二云石英片岩和黑云母大理岩等.白云母主要分布在晚太古代变质岩区，少量出现在二叠纪花岗岩区，相关岩性有：石榴石白云母斜长片麻岩、二云斜长片麻岩、二云石英片岩、二云石英片岩和含白云母石英岩等.角闪石主要分布在晚太古代斜长角闪岩和斜长角闪片麻岩区域，少量在蓟县纪大理岩区.绿泥石矿物信息主要在晚太古代钠长绿帘绿泥片岩和黑云斜长片麻岩中，少量出现在长城纪绿泥石千枚岩和二叠纪的钠长绿泥石英千枚岩中.石榴石矿物只出现在西北部的晚太古代地层中，主要赋存于石榴石白云母斜长片麻岩和石榴黑云斜长片麻岩中.阳起石仅在蓟县纪蛇纹透辉白云质大理岩中有零星分布.

图10 矿物提取结果.（a）Bi；（b）Mus；（c）Am；（d）Chl；（e）Gt；（f）ActFig.10 Mineral extraction results: (a) Bi; (b) Mus; (c) Am; (d) Chl;(e) Gt; and (f) Act

3.4 野外调查精度评价

通过42个点的野外调查、采样和实验室薄片鉴定等步骤，对本文的矿物信息提取结果进行精度评价.由于野外勘查条件限制，实际采样点未能均匀分布在矿物显现区，根据岩性分布情况进行了补充采样.部分样品薄片鉴定结果如表5和图11.

图11 部分样品显微照片.（a）云母石英片岩；（b）石榴石白云母斜长片麻岩；（c）蛇纹石透辉石白云质大理岩；（d）钠长石绿泥石石英千枚岩Fig.11 Micrograph of some samples: (a) mica quartz schist; (b) garnet muscovite plagioclase gneiss; (c) serpentine diopside dolomitic marble; (d) albite chlorite quartz phyllite

表5 部分样品薄片鉴定结果Table 5 Identification results of some sample slices

通过野外实测生成测试集，评价模型的泛化能力.精度评价结果（表6）显示，提取精度较高的为黑云母、白云母、角闪石，较低的为绿泥石、石榴石、阳起石.黑云母在黑云母大理岩及黑云斜长片麻岩中质量分数最高可达20%，且分布较广，实际采样点多覆盖在这些岩性上，总体精度为85.4088%、Kappa为0.7779.白云母在石榴石白云母斜长片麻岩及二云石英片岩中质量分数最高可达30%，二云斜长片麻岩中质量分数也达到了10%，提取精度为 84.7640、Kappa为 0.7833.角闪石在斜长角闪岩中质量分数最高可达到75%，斜长角闪片麻岩中质量分数最高可达到55%，矿物特征较为鲜明，提取精度为 85.7308、Kappa为0.7748.绿泥石在绿泥石英千枚岩、含绢云母石英岩、斜长角闪岩和斜长角闪片麻岩等岩性中均有体现，但质量分数一般不超过15%（仅在极少量的钠长绿泥石英千枚岩中质量分数达到了48%），因此，容易被其他高含量矿物干扰，提取精度为70.6933、Kappa为0.5938.石榴石在石榴石黑云斜长片麻岩中质量分数为5%、在石榴石白云母斜长片麻岩中质量分数为3%，在其他岩性中均未出现，所以提取时样本训练程度不足导致它的总体精度仅为65.5992、Kappa为0.5462.阳起石仅在蛇纹透辉白云质大理岩中出现了5%，提取精度为66.7509、Kappa为0.5560.

表6 精度评价Table 6 Accuracy evaluation

测试样本的拟合决定系数R2被用来评价模型的预测准确性[22].各采样点实测值和预测值对比，Bi、Mus、Am 及 Chl的R2分别为 0.971、0.962、0.966及0.965，预测值与实测值基本吻合；Gt及Act的R2分别为0.935及0.932，拟合效果不佳，但是由于在研究区石榴石和阳起石分布非常少，所以提取结果依然具备实用价值.

3.5 与其他方法提取精度对比

将本文方法与比值+阈值分割法、比值+SVM法提取的总体精度进行对比（表7）.各方法采用的比值运算公式一致，Bi、Mus、Am、Chl、Gt及Act的分割阈值分别为：1.4、3.1773、1.1406、3.01、1.405及2.25.采用SVM分类时的样本及特征与本文方法一致，核函数选择径向基函数(Radial basis function，RBF)，参数 σ分别取测试错误率最小的 14、11、13、19、8及 12，C统一取值为 9.黑云母、角闪石及绿泥石出露较多且较为集中，经多尺度分割之后精度提高幅度较大，本文方法比比值+阈值分割法分别提高了8.0542%、8.7624%及6.0339%；比比值+SVM法分别提高了4.2747%、4.6306%及4.2401%.而白云母、石榴石、阳起石出露较少，精度提高幅度不大，本文方法比比值+阈值分割法分别提高了3.481%、3.757%及3.9662%；比比值+SVM法分别提高了1.777%、2.4106%及2.5158%.

表7 本文方法与其他方法提取精度对比Table 7 Comparison of extraction accuracy between the present method and other methods

4 结论

根据变质矿物的特征性光谱特征进行比值运算增强ASTER影像，并基于光谱特征和变差函数纹理进行多尺度分割，然后，通过RF提取变质矿物信息，最后经过野外验证进行精度评价.结果表明，黑云母、白云母、角闪石等在ASTER影像上具有鉴定性特征，提取精度可分别达到85.4088%、84.7640%、85.7308%；而绿泥石、石榴石、阳起石作为次要矿物，鉴定时被主要造岩矿物干扰，提取精度达到60%以上.

基于野外调查和ASTER影像处理进行变质矿物的提取可有效提高变质岩野外调查效率和精度.本文方法可较为准确的鉴定ASTER的变质矿物；可为其他遥感影像提取矿物提供借鉴；也可作为辅助地质调查的有效手段.与其他类似研究相比，基于变差函数的多尺度分割能增强形态特征对矿物信息的区分能力；RF对训练数据的统计假设少、对矿物混合导致的噪声不敏感、分类变异性低，对多种矿物组成的岩性具有较强的分析能力.

下一步研究方向：①比值运算增强目标矿物的同时也增强了其他具有相似光谱特征的矿物，导致误分类，可通过原始岩性进行增强约束；②多尺度分割对对象内部的同质性要求较高，然而实际像元为多种矿物的混合，过分割现象较多，可通过解混像元突出主要信息；③RF的分类精度取决于树的棵树，可通过寻优算法代替人工判定树的棵树.