基于国产高分卫星数据的矿山环境变化检测

王立娟， 靳 晓， 贾虎军， 唐 尧， 马国超

(1.成都理工大学环境与土木工程学院，成都 610045； 2.四川省安全科学技术研究院，成都 610045； 3.重大危险源测控四川省重点实验室，成都 610045)

0 引言

遥感技术是一门综合性对地观测技术，具有大面积同步观测、信息获取速度快、周期短、综合性强等特点，它弥补了传统矿山监测技术手段的缺陷，能够全面、客观、有效、准确和动态地反映出矿山及周边区域的情况，且不受环境、人工等因素的影响。特别是随着我国航空航天技术的不断发展，遥感影像的空间分辨率得到了大幅度提高，其中优于1 m的高空间分辨率遥感影像(以下简称“高分影像”)可为遥感技术在矿山监测方面的应用提供有力的数据支撑。

近几十年里，SPOT5，IKONOS，TM和QuickBird等遥感影像在矿山监测中得到了广泛应用[1-2]，随着我国卫星技术的迅速发展，许多学者都尝试将国产高分遥感技术融入到矿山监测技术体系中。路云阁等[3]结合西藏自治区矿山遥感监测工作的特点，提出并实现了从国产卫星遥感数据管理、增强与校正、信息提取、统计分析等一体化解决方案，为大规模矿山遥感动态监测提供了技术支持; 安志宏等[4]利用资源一号02C星遥感影像，以河北承德多金属矿区和江西寻乌稀土矿区为实验区，开展了1∶5万比例尺的矿山遥感监测应用研究，并针对所遇到的问题提出了解决方案，为02C星数据在矿山遥感监测的应用起到了示范作用; 魏江龙等[5]则以会理多金属矿区为研究区，以高分一号卫星影像为数据源，对研究区内的矿山分布和开采情况进行提取，掌握了研究区内矿产资源开发现状和矿山环境状态。然而，目前国产高分卫星影像在矿山监测方面的应用仍然以专业技术人员目视解译为主。随着矿山监测技术的智能化、定量化要求不断提高，传统的以目视解译为主的遥感技术已经不能满足矿山监测的需求，更加高效和定量化的遥感技术亟待融入到矿山监测的技术体系中。

为此，本文以我国高分二号(GF-2)遥感影像为数据源，根据矿山监测的目标，构建多源特征集，在变化向量分析法(change vector analysis,CVA)的基础上自动选择训练样本，然后利用极限学习机(extreme learning machine，ELM)提取变化信息，提出一种自动化程度比较高的面向对象的变化检测方法(CVA-ELM)，既提高了检测精度，又避免了提供先验知识。以四川省攀枝花市的典型矿山为例，开展矿山及周边环境的动态监测实验，为矿山遥感监测提供了应用范例。

1 矿山环境变化检测方法

1.1 检测目标确定

确定检测目标是遥感影像变化检测技术的首要步骤，是遥感数据源、变化检测算法和模型选择的基础。矿山遥感监测一方面是为了掌握矿山自身的生产状态，另一方面是为了了解矿山周边环境的变化情况。采场是矿山生产活动的主要场所，也是矿山监测的重要区域。除采场外，尾矿库和排土场作为矿山生产过程中的重大危险源，同属于矿山的重要要素，是矿山监测的重点关注对象。在尾矿库及周边区域的动态监测中，库区的扩张、干滩和坝体的变化、库区周边山体植被的变化、下游建筑物和重要设施的变化为主要监测目标； 而对于采场(排土场)及周边区域的动态监测，采场(排土场)的扩张、周边区域植被的变化、建筑的增加等则为重点监测目标。

1.2 影像预处理与特征提取

遥感影像变化检测的数据预处理的主要目的是使同一位置上相同地物在不同时相的影像中地理坐标和特征属性相同，而不同地物则表现出不同的特征属性。因此，需对遥感影像进行预处理，除大气校正、正射校正和融合等基本预处理过程外，几何配准和相对辐射校正是多时相遥感影像预处理的重要内容。

根据矿山环境的检测目标，本文在检测中首先提取了对尾水、干滩和蓝色屋顶的建筑比较敏感的归一化水体指数(normalized differential water index, NDWI)、对植被比较敏感的归一化植被指数(normalized differential vegetation index, NDVI)以及对高亮反射率的建筑比较敏感的亮度特征，并将它们同原始光谱数据一起作为变化检测的特征集。

1.3 CVA-ELM变化检测方法

变化检测方法按照是否需要提供先验知识，可以分为非监督变化检测方法与监督变化检测方法[6]2大类。非监督变化检测方法虽然不需要提供先验知识，自动化程度相对较高，但是检测结果受影像的成像环境及变化阈值的影响较大，其常用的方法有差值法、比值法、CVA法[7]和主成分分析法(principal component analysis, PCA)[8]等; 而监督变化检测方法则能够在一定程度上降低大气、传感器及周围环境对变化检测结果的影响，且不需要设置变化阈值[9]，其中决策树(decision tree, DT)[9]、支持向量机(support vector machine, SVM)[10]和ELM[11]等方法均具有良好的性能，是监督变化检测方法中常用的分类器。2种类型的变化检测方法都有自身的优势与缺点，将二者融合，取长补短，不仅能够提高变化检测方法的自动化程度，而且可以增加算法的鲁棒性。

CVA利用多特征空间中的变化矢量表示地物的变化，算法简单、易于实现，与差值法和比值法相比，能够综合考虑多个特征波段，得到的变化检测结果更客观、更接近地物的真实变化。ELM是一种用于训练单隐藏层反馈神经网络(single hidden layer feedback neural network，SLFNs)的学习算法，它对于能够无限可微的激活函数，可以随机设置其输入权值和隐藏层偏置，不需要进行迭代调整[12-13]，这使其不仅训练速度比传统的基于梯度算法快，而且避免了诸如局部最小以及过拟合等问题。ELM的显著优势是训练速度快、检测精度高，特别适用于数据量较大的影像的处理。因此，本文将CVA和ELM结合起来，用于开展大范围矿山遥感动态监测。

由于高分影像中地物往往以多个像元聚集成区域的方式呈现，以像元为单位的变化检测方法在高分影像的处理中存在一定的局限性[14-15]，而面向对象的方法以分割对象为单位进行处理，更符合地物在影像中真实形态，所以本文构建了一种面向对象的自动变化检测方法。基本流程如图1所示。

图1 面向对象的变化检测方法流程Fig.1 Flow chart of object-based automatic change detection method

主要包括以下步骤：

1)对预处理后的多时相遥感影像进行均值漂移(mean shift，MS)分割，并根据矿山的特点进行特征提取，构建对象级的多源特征集。

2)在CVA的基础上自动选择训练样本[16]。样本的自动选择方法通常建立在以下假设上： 如果2个时相的差值影像中对象的值越大，则对象属于变化类别的可能性就越大； 相反如果值越小，对象属于未变化类别的可能性就越大。因此，基于CVA的样本自动选择方法的步骤包括： 首先获取CVA差值影像，然后对所有的对象进行由大到小的排序，最后选择前α·N的对象为变化类别的样本，后α·N的对象为未变化类别的样本，其中N为差值影像中对象的个数，α为相应的阈值，α∈(0,0.5)。

3)以所构建的多源特征集作为输入影像，利用自动选择的训练样本对ELM进行训练，然后获取监测区域的变化信息。

2 实验与分析

为了验证CVA-ELM变化检测方法的优越性，并分析比较各种变化检测方法在矿山监测方面的适用性，选择矿山周边的小部分区域作为实验区域，用于评价各种方法的检测精度。选择CVA、基于像元的差值法(diff_pixel)、面向对象的差值法(diff_OB)、面向对象的变化向量分析法(CVA_OB)、面向对象的支持向量机(SVM_OB)共5种方法为CVA-ELM的对比算法，其中SVM_OB也是在CVA自动选择训练样本的基础上利用支持向量机提取变化信息的方法。

本次实验选择2015年2月24日和2016年2月14日2景GF-2遥感影像，制作了全色和多光谱数据的融合影像，影像大小为400像元×400像元，覆盖某矿山周边160 000 m2的区域。制作的B3(R),B2(G),B1(B)波段融合影像如图2所示。为了对各种变化检测方法进行定量化的精度评价，根据Google Earth与实地调查相结合，获得了此区域的真实参考变化图(图3)，其中图3(a)中白色为变化区域，黑色为未变化区域，图3(b)为参考变化图与2015年遥感影像的叠加图，红色为变化区域。

(a) 2015年2月24日 (b) 2016年2月14日

图2实验区域2期影像

Fig.2Imagesofexperimentregionintwoperiods

(a) 参考变化 (b) 参考变化叠加

图3实验区域参考变化

Fig.3Referencechangemapofexperimentregion

对原始遥感影像进行校正，然后在实验区域的影像上均匀地选择30个同名点，配准误差均控制在0.5个像元之内。相对辐射校正以2015年2月24日遥感影像作为基准，采用线性回归分析法完成。

面向对象的变化检测方法以分割对象作为处理单位，为了得到相同的对象，采用MS分割方法对2个时相的影像进行复合分割，得到影像处理的分割对象。根据矿山的特点及监测目标，利用遥感影像提取特征构建多源特征集。特征集构建的完备程度直接影响变化检测的结果。在本实验中提取NDVI，NDWI以及亮度特征与原始的光谱特征共同构建多源特征集，各个特征影像如图4所示。

(a) 2015年2月24日 NDVI (b) 2015年2月24日 NDWI (c) 2015年2月24日亮度特征

(d) 2016年2月14日 NDVI (e) 2016年2月14日 NDWI (f) 2016年2月14日亮度特征

图4实验区域特征影像

Fig.4Featureimagesofexperimentregion

依照1.3中的方法进行训练样本的自动选取。为了控制训练样本的数量在一个合理的范围，参数α的范围被设定为[0.05,0.15]。由于SVM的检测性能好，鲁棒性强，因此被用于选择最佳参数。图5显示了不同参数选取的训练样本利用SVM进行变化检测被错误检测的像元数量，当α被设定为0.12时，被错误检测的像元数量最少。

图5 不同参数选择的训练样本得到错误检测像元的数量Fig.5 Number of error detected pixels which were obtained by different training samples with different parameters

表1为各种变化检测方法的精度评价，包括总体精度、Kappa系数、虚检率和漏检率。由于ELM的输入权重矩阵和隐含层偏差是随机获取的，所以CVA-ELM的各项精度为10次运算结果的平均值。

表1 各种方法的变化检测精度Tab.1 Accuracy of different change detection methods

从表1中可以看出，2种面向对象的自动变化检测方法(SVM-OB和CVA-ELM)的检测精度均高于其他方法，其中CVA-ELM的精度更高，且其运算时间远远低于SVM-OB。与基于像元的方法相比较，CVA-OB的总体精度和Kappa系数均低于CVA，而diff-OB的总体精度和Kappa系数均高于diff-pixel。这主要是因为非监督变化检测方法需要设置变化阈值，CVA-OB与CVA均使用最大期望算法(expectation maximization algorithm，EM)自动获取变化阈值，EM算法更适用于基于像元的图像处理方法； diff-OB和diff-pixel为人工获取变化阈值，不同的方法变化阈值也不同，而目前缺少针对面向对象的变化阈值自动获取算法，这也限制了面向对象的非监督变化检测方法的发展和应用。

图6为实验区域不同变化检测方法的检测结果，与图3的参考变化图相比，从视觉上看，CVA-ELM的检测结果(图6(f))与参考变化图最接近，这验证了CVA-ELM的良好检测性能和在矿山安全监测的适用性。CVA的检测结果(图6(a))与CVA-OB的检测结果(图6(b))相比，后者的椒盐噪声减少了很多，但是漏检的像元也比较多。同样，diff-pixel(图6(c))和diff-OB(图6(d))的检测结果相比，后者检测结果中虚检的像元数目较少，体现了面向对象的图像处理方法的优势。

(a) CVA (b) CVA-OB (c) diff-pixel

(d) diff-OB (e) SVM-OB(f) CVA-ELM

图6各种方法的变化检测结果

Fig.6Changemapobtainedbydifferentmethods

3 应用案例

四川省攀枝花市米易县地处攀西钒钛磁铁矿区的核心部位，拥有丰富的矿产资源，县域内采场、排土场和尾矿库数量多且分布分散。本文分别以攀枝花市米易县域内的尾矿库、排土场及采场为对象，开展矿山环境的变化检测。本次检测选择2015年2月24日和2016年2月14日的GF-2遥感影像，经过正射校正、辐射校正、几何纠正、影像融合及配准等预处理，以融合影像为数据源，利用CVA-ELM变化检测方法实施矿山及周边环境的动态监测。

3.1 尾矿库及周边环境变化检测

万年沟尾矿库是目前四川省的第一大库，是比较有代表性的大型尾矿库。本文以万年沟尾矿库为实验区开展尾矿库及周边环境的变化检测。检测区域影像大小为3 579像元×2 884像元，平面面积约10 840 956 m2，覆盖了万年沟尾矿库及周边的大范围区域。2个时期的遥感影像如图7所示。

(a) 2015年2月24日 (b) 2016年2月14日

图7万年沟尾矿库及周边区域2个时期遥感影像

Fig.7ImagesofWanniangoutailingpondandthesurroundingareasintwoperiods

图8(a)展示了万年沟尾矿库及周边区域2015—2016年间的变化图，其中白色区域为变化区域，黑色区域为未发生变化区域，红色线条为变化区域的边界。为了更加直观地展示变化区域的分布，判断变化区域的类别，将变化信息叠加在2015年2月24日的遥感影像上得到变化检测叠加图，如图8(b)所示。可以看出万年沟尾矿库及周边区域的主要变化为尾矿库扩张、干滩变化、坝体变化、植被减少及建筑的变化等，变化面积约为372 047 m2，约占总面积的3.43%。

(a) 变化检测 (b) 变化信息叠加

图8万年沟尾矿库及周边区域2015—2016年变化

Fig.8ChangemapofWanniangoutailingpondandthesurroundingareasfrom2015to2016

3.2 排土场及周边环境变化检测

以威龙州排土场为实验区开展排土场及周边环境的变化检测。检测区域影像大小为1 584像元×1 181像元，平面面积大约为1 870 704 m2，2个时期的遥感影像如图9所示。

(a) 2015年2月24日 (b) 2016年2月14日

图9威龙州排土场及周边区域2期遥感影像

Fig.9ImagesofWeilongzhoudumpandthesurroundingareasintwoperiods

图10中分别展示了威龙州排土场及周边区域的2015—2016年间的变化图及变化信息叠加图。

(a) 变化检测 (b) 变化信息叠加

图10威龙州排土场2015—2016变化

Fig.10ChangemapofWeilongzhoudumpandthesurroundingareasfrom2015to2016

由图9和图10可知，威龙州排土场及周边区域在2015—2016年间的变化主要为排土场的扩张、植被的减少以及排土场下游房屋的减少，变化面积约为86 901 m2，约占总面积的4.64%。威龙州排土场下游区域为政府划定的拆迁区域，图10(b)直观显示了此区域已拆迁和未拆迁的建筑，为政府监管和矿山企业安全生产提供了有效的数据支撑。

3.3 采场及周边环境变化检测

以冰花兰采场作为实验检测区域开展采场及周边环境变化检测，检测区域高分影像的大小为1 957像元×2 458像元，平面面积约4 810 306 m2，2个时期的遥感影像如图11所示。

(a) 2015年2月24日 (b) 2016年2月14日

图11冰花兰采场及周边区域2期影像

Fig.11ImagesofBinghualanstopeandthesurroundingareasintwoperiods

图12为冰花兰采场在2015—2016年间的变化图。从图中可以看出，冰花兰采场在其西南方向开辟了新的开采区域，而原有的开采区域未发生明显变化，基本停止开采工作。从2015年至2016年间冰花兰采场及周边区域的主要变化为采场的扩张、周边尾矿库的扩张以及植被的变化，变化面积约为189 596 m2，约占总面积的3.94%。

(a) 变化检测 (b) 变化信息叠加

图12冰花兰采场2015—2016年变化

Fig.12ChangemapofBinglanhuastopeandthesurroundingareasfrom2015to2016

4 结论

1)本文根据矿山监测面向定量化和自动化的发展要求，将遥感影像变化检测技术应用于矿山及周边环境动态监测中，并比较分析了几种常用的变化检测方法在矿山监测中的适用性。

2)根据矿山环境的特点，以国产高分二号卫星影像为数据源，构建一种自动化程度比较高的面向对象的变化检测方法。通过实验证明，本文所构建的面向对象的变化检测方法与常用的几种变化检测方法相比，具有自动化程度高，检测性能好的优势，更适用于矿山监测应用。

3)立足于矿山生产安全和矿山生态环境监测需要，明确矿山监测的重点关注对象和监测目标，并分别以攀枝花市米易县的典型尾矿库、排土场和采场为例，开展矿山及周边环境的变化检测实验，为矿山实施大规模遥感动态监测提供了范例。