一种基于高分辨率卫星影像的露天采矿场信息提取方法与流程

本发明涉及遥感影像处理技术领域，具体地指一种基于高分辨率卫星影像的露天采矿场信息提取方法。

背景技术：

随着近年来大力发展遥感卫星事业，我国已拥有自主产权的高分辨率遥感卫星高分一号(全色影像分辨率为2m，多光谱分辨率10m)和高分二号(全色影像分辨率为0.8m，多光谱分辨率2m)，然而对于区域尺度的露天采矿场，目前运用高分辨率遥感影像进行目标识别的方法主要仍然以人工目视解译为主。遥感卫星影像空间分辨率的极大提高，能够更为清晰地表现更细致丰富的地物要素信息，遥感与面向对象图像处理算法的结合可以提高区域尺度下露天采矿场的信息提取效率和准确度。

基于国产高分辨率卫星遥感影像的露天采矿场信息提取方法的核心在于：①国产高分辨率卫星遥感影像预处理；②影像分割层次结构的建立；③遥感与机器学习算法在露天采矿场信息提取中的结合。

经过我国学者长期不懈的辛勤研究，在中低分辨率遥感影像(空间分辨率大于15m)的研究成果基础上，露天采矿场信息提取相关的研究得到进一步深化、细化，并取得了不俗的研究成果。2004年，王晓红、聂洪峰等研究人员以江西省崇义县为研究区，基于quickbird影像和spot-5两种高分辨率遥感卫星影像，使用目视室内判读的方法，开展了南方高植被覆盖山区内，复杂野外调查条件下大规模采矿区域矿山开采信息提取、非法采矿活动监测和矿山开发状况评估等相关研究，相较于过去使用中低分辨率遥感影像，判读效果更加准确，进一步节约监测调查成本(高分辨率卫星数据在矿山开发状况及环境监测中的应用效果比较，国土资源遥感，2004年第16期)。目前，我国自主研发的高分辨率遥感卫星包括两个系列：资源系列卫星和高分系列卫星。2016年，熊前进、柴小婷等研究人员使用国产资源3号卫星高分辨率遥感影像，通过最佳波段组合指数值计算确定了矿山开采状态、开采范围(资源3号卫星影像在矿山监测应用中的数据处理方法，武钢技术，2016年第54期)；而魏江龙、周颖智等研究人员以会理多金属矿区为研究区，基于国产高分一号遥感影像为数据源，通过数据预处理、建立遥感解译标志、统计分析等研究工作，提取了研究区内矿山分布和开采信息(基于高分一号数据的矿山遥感监测——以会理多金属矿区为例，有色金属:矿山部分，2016年第68期)；在此基础上，马秀强、彭令等研究人员也尝试将高分二号数据应用于湖北省大冶矿山地质环境调查中，他们通过分析地表地物光谱特征和空间特征，采用波段运算和密度分割，定性分析矿区水体污染级别，综合分析矿区土地利用状况，并对重要矿集区的矿山开发状况和环境变化进行监测(高分二号数据在湖北大冶矿山地质环境调查中的应用，国土资源遥感，2017年第29期)。

除传统的人工目视解译方法之外，进入高分辨率遥感影像时代之后，仍然有不少研究人员在继续监督分类和非监督分类在矿山遥感调查和监测中的应用方面的研究，例如陈兴杰在2017年使用国产gf-1号数据，对比了最大似然法、最小距离法、支持向量机等几种监督分类方法在土地利用和覆盖方面的应用效果，并讨论了监督分类方法在高分影像中的应用各个方面的差别(基于gf-1号卫星影像的监督分类方法比较，矿山测量，2017年第45期)；此外，代晶晶、王登红等研究人员在2013年使用ikonos影像，运用isodata非监督算法对稀土矿开采范围及周边河流污染程度进行评估(ikonos遥感数据在离子吸附型稀土矿区环境污染调查中的应用研究——以赣南寻乌地区为例，地球学报，2013年第3期)。

传统的遥感影像解译方法本质上还是从基于像素的角度来理解遥感影像的，这种方法只能反映单个像素的光谱特征，无法从整体上理解影像的特点，没有利用图像的单个对象特征以及对象之间的联系，对于高分辨率遥感图像，传统的基于像素的方法会导致影像信息的大量冗余和资源的浪费。

技术实现要素：

本发明的目的就是要解决上述背景技术存在传统遥感影像解译方法存在会导致影像信息的大量冗余以及资源的浪费的问题，提供一种基于高分辨率卫星影像的露天采矿场信息提取方法。

本发明的技术方案为：一种基于高分辨率卫星影像的露天采矿场信息提取方法，其特征在于：包括以下步骤：1)、对拍摄有目标露天采矿场的国产高分辨率卫星影像进行预处理；

2)、计算露天采矿场和出现在卫星影像内的其他区域背景环境要素的最优分割尺度，构建卫星影像的分割层次结构；

3)、基于经过预处理后的卫星影像，计算每个分割对象的特征因子集，分层掩膜去除与特征因子集有关的区域背景环境要素；

4)、通过随机森林算法构建决策树，根据特征因子集对去除了与特征因子集有关的区域背景环境要素后剩下的混合要素进行分类，获得露天采矿场信息。

进一步的所述的步骤1)中，对卫星影像进行预处理的方法为：对拍摄有目标露天采矿场的国产高分辨率卫星影像进行辐射定标、大气校正、正射校正、图像融合处理的处理操作。

进一步的所述的步骤2)中，计算最优分割尺度的方法为：采用欧几里得二指数法来定量判断各种地质环境要素的最优分割尺度，对于无法通过欧几里得二指数法来确定区分的地质环境要素，通过影像分割的方法来划分对象单元。

进一步的所述的步骤2)中，构建卫星影像的分割层次结构的方法为：参考分割尺度值、提取地质环境要素由易到难的原则构建信息提取层，通过对不同地质环境要素采用不同的分割尺度，确定适合每种地质环境要素的分割尺度，在此基础上建立可靠的分类规则集，无法细分的所有地质环境要素设置为最底层。

进一步的所述的步骤3)中，所述的特征因子集包括光谱特征、几何特征、纹理特征和空间特征。

进一步的所述的光谱特征包括各波段灰度平均值、各波段灰度标准差、亮度值、归一化植被指数和最大差分。

进一步的所述的几何特征包括长宽比、密度、紧致度、形状指数、矩形拟合度、不对称性、圆度。

进一步的所述的纹理特征包括基于灰度共生矩阵的同质性、对比度、相异度、熵、均值、标准差、相似度以及基于灰度差分矢量的熵、均值、对比度。

进一步的所述的空间特征包括距工矿仓储用地距离和邻域接触关系。

进一步的所述的步骤4)中，通过随机森林算法构建决策树的方法为：采用随机选取分裂属性集的方法，从整个训练区数据集中选择2/3大小的训练样本集，余留的1/3数据作为袋外数据来估算随机森林算法的学习能力，重复这个过程ntree次，300≤ntree≤500，最后利用生成的ntree棵决策树的预测结果进行投票选择最终的分类结果。

本发明的优点有：1、相较传统的室内人工目视解译，这项技术方案提供了一种在大区域尺度上识别露天采矿场的方法，同时一定程度上提高了信息提取效率；

2、通过影像分割层次结构，满足了对各类地质环境要素的分层提取，同时得到每一层的特征因子规则集，为定量化分析不同地物要素在国产高分辨率卫星遥感影像上的特征规律提供了依据；

3、由于使用了随机森林分类算法，这项技术方案能够满足大数据量国产高分辨率卫星遥感影像上露天采矿场的精确提取，这是其他方法还未实现的。

本发明通过对卫星影像进行层次分割处理，利用随机森林分类算法，能够快速且精确获得露天采矿场的信息，处理方式简单、高效，解决了现有基础存在信息冗余和资源浪费的问题，具有极大的推广价值。

附图说明

图1：露天采矿场信息提取方法流程图；

图2：遥感影像预处理流程；

图3：各类地质环境要素参照样本图；

图4：本实施例建立的6个地物信息提取层次结构图；

具体实施方式

下面结合附图和具体实施例对本发明作进一步的详细说明。

如图1所示，为本实施例的基于国产高分辨率卫星遥感影像的露天采矿场信息提取方法流程图。本实施例提供了一种露天采矿场信息提取方法，包括以下步骤：

第一步：获取目标区域的原始国产高分辨率遥感影像，并进行预处理。

获取目标区域的原始国产高分辨率遥感影像。具体实施中，原始高分辨率遥感影像可选用高分一号卫星或高分二号卫星高分辨率遥感影像，他们的多光谱波段都含有4个波段的信息(分别为蓝、绿、红、近红外波段)，具体有效载荷技术指标如表1-2所示：

表1高分一号卫星有效载荷技术指标

表2高分二号卫星有效载荷技术指标

本设计实施例选取了豫西洛阳市栾川露天钼矿开采矿区2016年01至2017年10月的高分一号卫星遥感影像，用于作为对目标区域露天采矿场的高分辨率遥感影像数据。具体实施中，该步骤还包括国产高分辨率遥感影像预处理，遥感影像获取后，通常需经过预处理后才能作为信息提取的数据源使用，如图2所示，国产高分辨率遥感影像预处理流程包括辐射定标、大气校正、正射校正和图像融合等操作。

第二步：确定最优分割尺度。本实例研究中采用的是欧几里得二指数法来定量判断各种地物的最优分割尺度，具体实现过程如下。

1、在arcgis中分别对大范围地质环境背景要素(居民地、林地、耕地)选取样本作为参照多边形，并计算面积大小。各类地质环境要素参照样本图如图3所示：

2、本次实验设置融合后gf-1影像四个波段权重均为1，当林地分割尺度范围由50到450，居民地和耕地分割尺度范围由50到250；紧致度因子compactness参数为0.5；形状参数shape小于等于0.5，取值0.1、0.3、0.5进行多尺度分割实验。将试验的分割结果矢量文件导入arcgis中，并与其相对应的参照多边形进行叠置分析及计算，得到3种地质环境要素的最优分割参数表3：

表3各地环境要素最优分割参数

第三步：信息提取层层结构和信息提取规则建立。

综合地质环境类别最优分割参数及影像特征，参考分割尺度值、提取地物由易到难的原则构建信息提取层，通过对不同地物要素采用不同分割尺度，确定每种地物适合自己的分割尺度，在此基础上建立可靠的分类规则集，本实施例中建立的6个地物信息提取层次结构如图4所示，信息提取将按以下顺序进行。

本实施例中从可选用的特征因子集中选择了亮度(brightness)、密度(density)、灰度共生矩阵对比度(glcmcontrast)、近红外波段灰度均值(nirmean)、矩形拟合度(rectangularfit)、圆度(roundness)、形状指数(shapeindex)、归一化植被指数(ndvi)等特征因子。利用训练区数据集，根据分类层次逐层选取特征并建立相应的特征规则，其特征规则集见表4。

表4特征规则集

分层提取的结果将通过类层间的相关性将其他层地物继承到level1层，比如level3层要继承level4层中的林地，此时level3即为level4的父层，level3中的林地可定义为：class-relationtosubobjects>existenceof>林地；同理，此时level4层作为level3层的子层，level4层中的林地也可按子层关系定义为：class-relationtosuperobjects>existenceof>林地。通过这样的层间关系，最后获取分层掩膜的遥感影像图，其中其它4这一类包括了4大类(露天采场、尾矿库、工业广场、裸地)难以靠分层分割区分开的地质环境要素。

使用特征规则集，同样对测试区数据集进行分层分割，得到包含了道路、河流湖泊、居民地、林地、耕地和其它4等6类结果的分层分割结果图。

第四步：在分层分割结果的基础上，基于随机森林分类算法对其它4种地质环境要素进行分类，提取露天采矿场。

在豫西研究区训练区对其它4种地质环境要素进行样本选择，本节研究随机挑选了1245个样本对象，其中403个露天采场样本、292个尾矿库样本、270个工业广场样本、280个裸地样本，从特征因子集中选取23个特征参数参与模型建立，具体特征集参数如表5。

表5特征集汇总

1、光谱特征

用来描述影像对象的光谱信息，是由真实的地物和成像状态所决定的光学物理属性，与对象的灰度值相关，其包含内容及具体数学定义如下：

a、均值(mean)，由构成影像对象的所有n个像素的层值cli计算得到层平均值，即：

其中：cli——单个像素的光谱灰度值；

n——像素数量；

b、亮度(brightness)，包含光谱信息的影像层平均值的总和除以层数。对影像对象来说，就是光谱均值的均值，即：

其中：ci——单个影像层的光谱灰度均值；

nl——影像层数量；

c、最大差分(maxdiff)，计算最大差分，只要把影像对象的最大值减去其最小值即可，需要对对象所有通道的最大最小值都做比较，最后的结果再除以亮度。

c＝max[max(cli)-min(cli)]/b,i＝1,2,3,…,n

其中：cli——单个像素的光谱灰度值；

b——影像亮度值；

n——像素数量。

d、标准差(standarddeviation)，由构成一个影像对象所有通道的n个像素的像素值计算得到标准差，即：

其中：cli——单个像素的光谱灰度值；

ci——单个影像层的光谱灰度均值；

n——像素数量。

e、归一化植被指数(ndvi)，提取植被指数的原理是根据地物光谱反射率的差异，通过比值等运算以突出植被特征、提取植被类型或估算绿色生物量。基于绿色植物叶子的细胞结构在近红外波段具有高反射，而其叶绿素在红光波段具有强吸收性的物理特征，绝大部分植被指数都利用红波段和近红外波段的反射率信息，提取植被相关信息。归一化植被指数的值越高，植被覆盖率越高；反之，则越接近无植被发育的不透水面，即：

其中：nir——近红外波段；

r——红波段；

2、形状特征

用来描述影像对象本身的形状信息，是由真实的地物所决定的空间几何属性，包含内容及数学定义如下：

a、长宽比(length/width)，长宽比等于协方差矩阵的特征值的比值，较大的特征值是分数的分子，即：

r＝eig1(s)/eig2(s)，eig1(s)＞eig2(s)

其中：eig1(s)和eig2(s)分别是协方差矩阵的特征值；

b、密度(density)，密度d可以表示为影像对象面积除上它的半径，使用密度来描述影像对象的紧致程度。在像素栅格的图形中理想的紧致形状是一个正方形，一个影像对象的形状越接近正方形，它的密度就越高，这里n是构成影像对象的像素数量，半径采用协方差矩阵来近似计算，即：

其中：n——像素数量；

c、紧致度(compactness)，对象的长k乘以宽q，再除以对象的像素个数n，即；

p＝k×q/n

其中：p——紧致度；

k——对象的长度；

q——对象的宽度；

n——对象的像素个数；

d、形状指数(shapeindex)，数学上形状指数是影像对象的边界长度除上它的面积的平方根的4倍。使用形状指数s可以描述影像对象边界的光滑度，影像对象越破碎，则它的形状指数越大，即：

其中：s——可以描述影像对象边界的光滑度；

e——影像对象的边界长度；

a——影像对象的面积；

e、矩形拟合度(rectangularfit)，计算矩形适合性的第一步是生成一个和考虑对象相同面积的矩形，计算过程中考虑对象的长宽比，此后，比较处于矩形外的对象面积和处于矩形内的未被对象填充面积比，如果值为0，意味着这个矩形不合适；如果值为1，说明能很好地适合对象；

f、不对称性(asymmetry)，一个影像对象越长，它的不对称性越高，对于一个影像对象来说，可近似一个椭圆，不对称性可表示为椭圆的短轴s和长轴l的长度比，随着不对称性的增加而特征值增加，即：

k＝1-s/l

其中：k——不对称性；

s——影像对象短轴长度；

l——影像对象长轴长度；

g、圆度(roundness)，计算分割单元接近圆形的程度，最大半径与最小半径之差为0时，圆度为0，圆度越大表示形状越偏离圆形，

r＝wmax-wmin

其中：wmax——分割单元最大半径；

wmin——分割单元最小半径；

3、纹理特征

纹理在高分辨率遥感影像分了中占有较重要的位置，用来描述影像对象的光谱在空间分布上的规律反映了地物表面的颜色、几何形状和灰度的某种变化，最普遍使用的描述纹理的方法有灰度共生矩阵(greylevelconcurrencematrix，glcm)和灰度差分向量矩阵(greyleveldifferencevector，gldv)，在此基础上使用较多的纹理统计量包括以下几类：

a、均质性(homogenity)，表示矩阵中较大像素值集中于对角线的程度，即：

其中：hom——灰度矩阵均质性；

pi,j——灰度矩阵中某个像元的灰度值；

i——灰度矩阵行数；

j——灰度矩阵列数；

b、反差(contrast)，衡量的是影像中局部的变化数量，当领域内最大值和最小值之间的差异增大，它也随之指数级增大，表示了图像视觉效果的清晰与否，即：

其中：con——灰度矩阵反差性；

pi,j——灰度矩阵中某个像元的灰度值；

i——灰度矩阵行数；

j——灰度矩阵列数；

c、相异度(dissmilarity)，与反差呈线性相关，局部反差越大，相异性越大，即：

其中：

dis——灰度矩阵相异度；

pi,j——灰度矩阵中某个像元的灰度值；

i——灰度矩阵行数；

j——灰度矩阵列数；

d、熵(entroy),熵是图像所具有的信息量的度量，熵值大小与图像信息量的大小成正比，熵值越大表明图像的信息量越丰富，说明融合影像的质量越好，即：

其中：pi——灰度共生矩阵每一行的所有像元灰度值总和；

e、均值(mean)，代表了影像中各类地物的平均反射强度，是遥感影像中所有像元亮度值的平均值，即：

其中：pi,j——灰度矩阵中某个像元的灰度值；

n——灰度矩阵像元个数；

f、标准差(standarddeviation)，标准差处理的是特定的参考和相邻像素的组合，它和原始影像中简单的灰度级别的标准差不一样，即：

其中：

si,j——灰度矩阵的标准差

pi,j——灰度矩阵中某个像元的灰度值；

ui,j——灰度矩阵中某个选定作为参考的特定像元的灰度值；

i——灰度矩阵行数；

j——灰度矩阵列数；

4、空间特征

空间特征作为辅助信息提取的特征要素，通常使用的是除遥感影像之外收集获取的电子化土地利用规划图、地形图和地质图等。通过分析露天采矿场与特定特征地物要素之间的空间关系(如平面距离、高差、接触关系等)辅助提高露天采矿场信息提取精度。

对剩下来的混合要素进行分类，从而得到露天采矿场信息。随机森林算法在构建决策树的时候，采用随机选取分裂属性集的方法，从整个训练区数据集中选择2/3大小的训练样本集，余留的1/3数据作为袋外数据(out-of-bag,oob)来估算随机森林算法的学习能力，这个过程被称作oob估计；重复这个过程ntree次，ntree∈[300,500]且为正整数，最后利用生成的ntree棵决策树的预测结果进行投票选择最终的分类结果。除此之外，随机森林中单棵决策树的分类精度还受到了随机抽样的属性集维度mtry的影响，即单次从全体特征集中随机抽取的特征因子数目和类别，通常mtry等于特征因子总数的开方，随机森林模型的决策树数目为300，建立本实施例中随机森林模型并参与运算。

参考该区域土地利用规划图和数字高程模型，使用空间分析的技术手段来消除分类器逻辑误差，得到测试区露天采场(l1)、工业广场(l2)、尾矿库(l3)、裸地(l4)等四类地质环境要素信息提取结果，同时并采用混淆矩阵和kappa系数进行精度评价，对信息提取结果和测试区实际值进行对比。

需要说明的是，本发明实施例提供的精度评价方法只用于本发明实施例提供的遥感影像地物分类方法的精度评价，在实际施行中，由于实际地物分类数据为未知的，不需要对精度进行检验。

目前有两种方式用于分类精度验证：①混淆矩阵，②roc曲线。本文采用混淆矩阵的方法来机型精度验证，评价精度的指标包括混淆矩阵、总体分类精度、kappa系数、每类图斑的制图精度和用户精度。

通过计算经本设计实施例提供的露天采矿场信息提取方法得到的结果，其总体分类精度、kappa系数、用户精度和生产者精度表明本设计实施例的露天采矿场信息提取方法具有很高精度，

表6遥感影像地物信息提取混淆矩阵

其中：

(1)用户(使用者)精度

用户(使用者)精度表示在该次分类中，在分类图上，落在该类别上的检验点，被正确分类为该类别的比率。

(2)制图(生产者)精度

制图(生产者)精度表示在此次分类中，该类别的地面真实参考数据被正确分类的概率。

(3)总体精度

总体精度指所有正确分类的土地覆盖类别的检验点数所占总抽取的检核点数的百分比；即在混淆矩阵中对角线的所有数值和除以全部样本的总和。

式中，pii代表混淆矩阵的元素，n代表类别个数的总和。

(4)kappa系数

它是通过把所有真实参考的像元总数(n)乘以混淆矩阵对角线(xkk)的和，再减去各类中真实参考像元数与该类中被分类像元总数之积之后，再除以像元总数的平方减去各类中真实参考像元总数与该类中被分类像元总数之积对所有类别求和的结果

式中，n表示类别，n代表类别个数的总和(此指检验点数)，xii表示混淆矩阵对角线元素，xi+表示类别列总和，x+i表示类别的行总和。

通过采用上述方法对本实施例的提取方法进行精度验算，证明本实施例提供的信息提取方法具有良好的精度，能够准确将卫星影像中的露天矿场信息提取出来。

可应用于实际应用当中。本设计实施例提供的基于国产高分辨率卫星遥感影像的露天采矿场信息提取方法，通过对国产高分辨率卫星遥感影像的处理和影响对象层次结构的构建，逐层掩膜分解去除区域地质环境背景信息，最终基于随机森林分类算法从混合地物类别中提取出露天采矿场信息。最终信息提取结果精度较高，具有良好实用性。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等同物界定。