无人机遥感技术在北京首云铁矿储量动态监测中的应用

向 杰， 陈建平， 李 诗， 赖自力， 黄浩中， 刘 静， 谢 帅

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

0 引言

无人机(unmanned aerial vehicle，UAV)遥感技术是无人机驾驶飞行器技术、遥感传感器技术、遥测遥控技术、通信技术、POS定位定姿技术和GPS差分定位技术等一系列技术的组合[1]。近年来，UAV遥感技术广泛应用于国土资源环境领域，成为国土资源遥感领域的研究热点之一。例如，利用UAV遥感开展土地资源调查与分类； 进行林业资源分布调查与资源总量估算； 开展环境监测与污染追踪调查等[2-3]。由于UAV遥感和矿山开采条件的限制，在矿山储量动态监测方面应用较少[4-5]。如何快速、准确、廉价地进行矿山监测一直是国土资源遥感领域的研究难点。

本文以北京市密云区首云铁矿为研究区，利用无人机摄影测量技术(UAV structure from motion，SfM)获取不同时相的高精度数字地表模型(digital surface model，DSM)，快速准确计算矿山动用储量，从而实现对储量的动态监控。与传统方法相比，该方法具有准确、快速、性价比高的特点，对于中小范围地区的矿山监测具有重要的推广价值。

1 研究区概况

首云铁矿属于密云铁矿的一部分，位于北京市密云区巨各庄镇境内，距北京80 km，属于较典型的丘陵地带。首云铁矿始建于1959年，1970年建成投产，是集采矿、选矿、竖炉焙烧连续生产工艺为一体的国有中型冶金矿山企业，曾隶属于北京市冶金局和首钢总公司，是原首钢铁矿冶金部重点矿山企业之一。

首云铁矿为露天开采，该矿区开采面积大约为1.5 km2，平均海拔为200 m，矿区植被稀少(图1)，这些条件都为UAV遥感调查提供了良好的环境。在2014年8月—2016年10月期间，该矿山一直处于露天开采状态，并且该矿山矿石体重和矿石品位比较均匀，通过计算地表体积变化即可计算出矿山动用储量，这为UAV进行矿山储量监测提供了前提条件。

图1 北京市密云区首云铁矿位置Fig.1 Location of Shouyun iron mine in Miyun district, Beijing

2 UAV数据获取与处理

2.1 UAV数据获取

本次工作分别于2014年8月和2016年10月对首云铁矿进行了UAV野外调查。2次UAV野外调查期间的天气晴朗无风，非常适合UAV飞行。另外，尽量选择正午期间进行飞行，减少了山体影像对结果的影响。2次UAV调查都采用台湾观天科技生产的X5型小型电动固定翼UAV，该UAV翼展为1.2 m，起飞重量为1.5 kg，有效载荷为300 g，续航时间可达40 min。UAV搭载的自驾仪和GPS可以实现按照规划航线自动飞行和定时、定距拍照； 机载POS系统可同步记录拍照时UAV姿态以用于数据后续处理； 搭载的相机为Sony QX100镜头，其基本参数如表1。

表1 UAV搭载数码相机基本参数Tab.1 Parameters of camera equipped on UAV

2次UAV调查飞行高度都设定为海拔450 m，设定航向重叠率80%，旁向重叠率60%。2014年8月的调查中获取了499张带POS信息的高清影像，实际航向重叠率81%，旁向重叠率62%； 2016年10月的调查过程中获取了674张带POS信息的高清影像，实际航向重叠率80%，旁向重叠率45%，并利用DGPS在特征明显且无变化的区域测定了7个控制点(ground control point，GCP)。UAV飞行航线如图2所示。

图2 UAV野外调查与数据获取航线Fig.2 Field surveys and data acquisition lines of UAV

2.2 UAV数据处理

UAV数据处理主要是为了得到高精度的DSM和高空间分辨率的数字正射影像图(digital orthophoto map，DOM)。近年来市场上涌现了一大批UAV遥感影像处理软件，如国内的DPGrid，PixelGrid和DPMatrix等，国外的PhotoScan，Pix4Dmapper和Inpho等[6]。

本次DOM的获取采用PhotoScan V1.2.5进行流程化处理，生成了高清的DOM。为方便后续计算，DSM获取采用如下步骤： ①在PhotoScan V1.2.5中预处理，其过程包括添加影像—加载POS点信息—对齐影像—生成网格—添加GCP—优化相机对齐方式—建立密集点云—生成点云数据； ②将导出的点云数据导入到cloud compare中，去除掉植被或一些噪声点； ③将处理好的点云数据导入到ArcGIS中，采用自然邻域法进行插值，生成空间分辨率为1 m的DSM[7]。

要进行地表变化的计算，2期DSM之间的配准极为重要，因为它会直接影响最后的研究精度。为此，本文首先利用DGPS测量的7个GCP数据，在ArcGIS中对2016年数据进行配准； 其次，对比2014年与2016年的DOM，在已配准完成的2016年数据上，选取特征明显且没有发生改变的点作为2014年的GCP； 最后，利用选取出的10个GCP配准2014年数据。采用这种方法不仅很好地解决了2014年控制点缺失的问题，也使得2期数据互相匹配，提高了后期地表变化计算的精度。

2014年和2016年2期UAV摄影测量获取的点云及DSM如图3所示。

(a) 2014年点云 (b) 2016年点云

(c) 2014年DSM (d) 2016年DSM

图32期UAV摄影测量获取的点云与DSM

Fig.3CloudpointandDSMoftwoperiodsusingUAVdata

通过以上处理过程，获得了2期高清DOM与配准过后的1 m空间分辨率DSM。2016年DSM误差采用已经测定的7个GCP来计算均方根误差(root mean square error，RMSE)得到。由于2014年没有测量GCP，因此基于配准完成的2016年数据，在没有发现变化的稳定区域挑选出10个GCP来计算DSM误差。2014年8月DOM的空间分辨率可达到0.081 m, 获得的DSM误差为0.44 m； 2016年10月DOM的空间分辨率达到0.072 m，获得的DSM误差为0.24 m。2期DSM的误差都属于亚m级，能满足进行矿山储量估算的精度要求。

3 动用储量计算

地形变化算法(DSM of Difference，DoD)最早运用于地形地貌领域。20世纪90年代早期，瑞士地貌学家Lane等便开始利用多期DSM相减来定量描述地形的变化[8]。近年来，获取DSM (或DEM)的技术手段不断丰富，如： 卫星立体像对、合成孔径雷达、SfM和激光雷达等，DSM精度在不断提高，DoD算法也不断改进完善[9]。本研究通过SfM，分别获取了2014年与2016年首云露天矿山的DSM，采用2期DSM相减的方式即可表示矿山的开采量。由于该矿山矿石密度与矿石品位都比较均匀，与矿山开采量相乘即可得到储量变化，进而动态监测露天矿山的储量变化。

但是，当用DSM来表示实际的地表模型时存在不确定性δZ，即

Z真实=ZDSM±δZ,

(1)

式中：Z真实表示实际的高程；ZDSM表示DSM所表示的高程；δZ表示对于单幅DSM的误差，而要计算地形变化必须要考虑到2期DSM的传播误差，通常δZ受到很多因素的影响，如测量误差、采样密度、类型和插值方法等。

2003年英国地貌学家Brasington等对于传播误差δDoD提出的公式[10]被广为接受，即

(2)

式中：δznew表示较新DSM的误差；δzold表示较老DSM的误差。

为了更加准确地计算矿山动用储量，消除传播误差对于DoD结果的影响显得尤为重要，本文采用了最小检出限法和转换概率法2种方法进行计算。

3.1 最小检出限法

最小检出限法是将传播误差δDoD作为最小检出限，即2期DSM相减，高程值变化小于δDoD视为是由于误差造成，不是由于矿山开采造成的，只有当高程值变化大于δDoD时，才视其为真正的开采所致[8]。2014年与2016年这2期DSM直接相减得到原始DoD如图4(a)所示。如前所述，2014年DSM误差为0.44 m，2016年DSM误差为0.24 m，根据公式(2)可求得传播误差为0.5 m。将传播误差0.5 m作为最小检出限法的DoD结果如图4(b)所示。

(a) 原始DoD (b) 最小检出限法DoD

图4最小检出限法计算地形变化结果

Fig.4ResultsofDSMofdifferencebyusingthelimitofdetectionmethod

通过对比图4(a)和(b)可以看出，最小检出限法所摒弃的信息主要分布在矿区中部和北部。从DOM上可以清晰地看到，该区域主要是道路与厂房，确实属于稳定非开采区域。矿山的主采区分布在矿坑南缘，尤其是中南部区域，开采造成的高程变化达到了80余m。矿山的堆积区主要分布在矿坑中部以及西北部，高程变化在20 m左右。

3.2 转换概率法

同样也可以将DSM的不确定性转化成概率值，通过用户自主定义置信区间的方式来去除不确定性对最终结果的影响[10-11]。假设传播误差的临界误差为U，则式(2)可以转换为式(3)， 基于t检验来确定一个置信区间。

(3)

(4)

式中|Znew-Zold|可以由原始DoD的绝对值得到。

因此，可以将原始DoD转化成如图5(a)所示的T值分布图。然后，可以通过将t统计量与其累积分布函数相关联来计算纯粹由于误差而引起的高程变化，利用t的累计分布函数将其转化为高程变化概率分布，通过设定合适的概率作为阈值来更加准确地估计动用储量。本研究按照惯例选取95%作为置信区间，来求取转换概率法DoD，结果如图5(b)所示。

(a) T值分布 (b) 转换概率法DoD

图5转换概率法计算地形变化结果

Fig.5ResultsofDSMofdifferencebyusingtheconvertedprobabilitymethod

图5结果可以看出，其主采区和堆积区分布与最小检出限法结果相同，摒弃的信息也主要分布在道路与厂房区域。但是，转换概率法摒弃了更多的信息量，其空值区域范围较大。

为了定量对比2种改进方法与原始DoD的计算结果，除了对比体积变化，还通过查询首云铁矿2014年8月—2016年10月间的开采报告，该地区铁矿石平均密度大约为4.8 t/m3，平均含矿率约为10%，矿石平均品位约为30%。将2种方法所计算得到的总体积变化分别乘以平均密度、平均含矿率和矿石平均品位即可估算出动用储量，计算结果如表2所示。

表2不同方法计算储量变化结果对比

Tab.2Comparisonofresultsbydifferentmethods

方法堆积/m3开挖/m3总体积变化/m3储量变化/t原始DoD1 325 479-14 633 968-13 308 489-1 916 422最小检出限法(0.5 m)1 264 146-14 559 008-13 294 862-1 914 460最小检出限法摒弃量61 333-74 960-13 627-1 962转换概率法(95%)1 262 525-14 556 857-13 294 332-1 914 384转换概率法摒弃量62 954-77 111-14 157-2 039

通过表2可以看到，储量体积变化比较接近(13 294 332～13 294 862 m3)，并且都在一定程度上排除了误差对于体积变化的干扰(13 627～14 157 m3)； 转换概率法采用95%作为置信区间所摒弃的信息，要大于最小检出限法以0.5 m为检出限所摒弃的信息。查询实际开采资料发现，2014年8月—2016年10月期间，首云铁矿的实际开采量约为180万t，本次研究2种方法的计算结果与矿山实际开采量都非常接近，其误差约为6.4%。

4 结论

本文将无人机摄影测量的方法运用于露天矿的储量动态监测中，以获取多时相、高空间分辨率的DSM为基础，利用最新的地形变化算法快速计算出矿山的体积变化，进而估算出矿山动用储量，为实现矿山科学管理提供了新的思路。主要结论如下：

1)通过2种地形变化算法计算出的动用储量体积，通过乘以该矿山平均含矿率、矿石平均体重以及矿石平均品位，估算的该期间动用储量与矿山实际开采量较为接近，误差仅为6.4%，证实了该方法的实用性。

2)由于金属矿产的品位、含矿率和密度具有一定的差异性，该方法的精度还有待进一步提高，目前更加适用于露天开采的煤矿或各类非金属矿产。另外，由于无人机续航能力的限制，该方法比较适合中小范围矿山的持续监测。

3)无人机不仅可以获得高空间分辨的DSM，还能方便快速地获取高清的DOM，基于DOM与DSM，还能进行矿山环境监测、灾害巡查、矿山施工管理和边坡设计与监测等。无人机摄影测量技术的运用为矿山监测和科学管理提供了全新的思路。

志谢： 感谢北京市国土局和密云区首云铁矿的大力协助，以及中国冶金地质总局遥感院王乾工程师对DGPS测量的指导。