基于普通数码影像的露采矿山三维可视化分析

史雨川，李 浩，杨 彪，刘庆群

（1. 河海大学 地球科学与工程学院，江苏 南京 210098）

基于普通数码影像的露采矿山三维可视化分析

史雨川1，李 浩1，杨 彪1，刘庆群1

（1. 河海大学 地球科学与工程学院，江苏 南京 210098）

研究了基于普通数码影像的矿山三维可视化技术。首先将检校后的数码相机以物方控制的单模型摄影方式对目标拍摄，取得原始影像并逐一校正；然后自动进行像对定向、立体匹配，采集影像点云；最后根据影像匹配点云三维建模，并依据前期地形数据、终采境界设计图、地质勘察数据对露采矿山进行三维可视化分析，实现动态检测。

露采矿山；普通数码影像；三维可视化；近景摄影测量

本文采用近景摄影测量方法，基于普通数码影像研究露采矿山三维可视化技术，通过量测化检校的相机对矿山开采宕口进行拍摄，并对原始影像畸变校正，自动进行相对定向、绝对定向以及影像立体匹配，利用采集的匹配点云重建矿山三维模型，并根据前期地形数据、终采境界设计图进行可视化分析[1-3]，最终实现矿山影像数据采集、处理、三维重建和可视化的一体化实时动态检测。

1 数码相机量测化检校

由于普通数码相机不提供畸变系数和内方位元素，因此必须进行量测化检校，使之获取的数字影像具有可控的精度。为求得畸变系数k1、k2、k3、p1、p2，应用精密平面控制场，将构像畸变差改正模型纳入影像透视变换，并用未知数交替解算[4]。再根据三维控制场的物方坐标及经过畸变改正的控制点像片量测坐标，利用单张像片后方交会算法解算内方位元素f、x0、y0[5,6]。将解算后的内方位元素重新代入上述畸变系数求解过程，利用新的畸变系数继续解算内方位元素，如此交替迭代直到畸变系数、内方位元素变化较小时停止，流程如图1。

实验显示，普通数码影像畸变最大可达几十个像元，校正后中误差小于1个像元，且主距固定时畸变系数和内方位元素比较稳定。

2 矿山数据快速采集与三维建模

2.1 影像快速采集

图1 畸变系数和内方位元素检校交替计算流程

从露采矿山需求考虑，由于单幅影像已经可以覆盖目标，所以采用外业工作量较小的物方控制的单模型摄影方式。拍摄前，在矿山宕口布置若干个控制点，精确测定控制点的坐标，通过单张像片空间后方交会计算影像的外方位元素。采取立体拍摄的办法分别在左站和右站对区域目标进行外业拍摄，且交向角不宜太大，基线长度与拍摄距离比为1/10左右，既能保证精度，又便于影像处理，相邻像对的拍摄覆盖范围保证衔接即可，如图2。

图2 矿山摄影方法示意图

为满足像对绝对定向中模型变形改正要求，每个像对至少要有4个控制点，且均匀分布在像片四角，相邻像对共用2个像控点。控制点可以采用标志点布设在矿山上，也可以选取特征明显的点。使用免棱镜全站仪以自由设站方式精确测量其空间坐标，利用GPS测量其中任意2个以上控制点的大地坐标，使其在数据处理时控制点坐标能够转换为大地坐标系。矿

山单像对拍摄原始影像如图3所示。

图3 矿山原始影像

根据测定得到的相机畸变系数和内方位元素，利用畸变差改正模型，对原始数码影像像素重排，得到畸变校正后的影像。

2.2 像对自动定向参数解算

传统的相对定向元素求解是采用迭代法，优点是只需较少的相对定向点，且定向精度高；缺点是必须提供较准确的未知数初值，一旦初值设置不当，求解方程迭代就可能不收敛。因露采矿山采用徒手拍摄，条件多变且目标复杂，所以相对定向元素初值很难估算[7,8]。为此，利用控制点像片坐标和空间坐标，通过角锥体法单像空间后方交会计算得到左右像片外方位元素，分别设为和并依据外方位元素自动设定相对定向初值。在摄影测量坐标系中，左右片旋转矩阵为：

设R1、R2分别为左右片的旋转矩阵，则右片相对于左片的旋转矩阵为：

经解算得出连续法相对定向初值：

初值确定之后，还需要若干对同名像点。利用Forstner算子自动提取左右片中的特征点，利用相对定向初值建立的概略核线作为约束条件，然后进行特征点匹配，自动剔除小区域内视差异常的错误匹配点，得出正确的同名点。将同名像点坐标与初值代入迭代公式，自动进行相对定向元素解算。

为确定相对定向所得几何模型在物方坐标系中的绝对位置，需进行空间相似变换，即像对绝对定向。空间相似变换是7个绝对定向元素的非线性函数，仍采用迭代法进行解算。由于连续法相对定向中采用左片的像空间坐标系为基准，所以模型空间方位由左片外方位元素确定，故绝对定向元素初值为取2个控制点的实际边长和对应模型长度比值作为模型比例系数的初值，绝对定向点自动选取之前在计算影像外方位元素时已经人工选取的实测控制点，最后自动进行绝对定向元素的迭代解算。

2.3 矿山影像立体匹配

由于露采矿山摄影距离较近，且开采表面崎岖变化较大，所以拍摄位置的变化容易引起较大的目标构像变形。其中包括摄影比例尺的差异、影像断裂、同名子影像形状变化和位移，而且不同视点会导致目标的背景不同，因此常规的单基线立体匹配的可靠性和精度较低[9]。为提高匹配准确性，设计了一种基于种子点约束的金字塔多级概率松弛匹配方法。

获取种子点的目的是用其约束搜索范围以提高匹配精度及降低匹配代价。方法是对于待匹配的格网点根据其周围的种子点估计其右片的对应位置，将左右影像的种子点对看作匹配的线段，种子点对之间的格网点根据右片相匹配的线段及格网点个数均分估计待匹配点的初始位置，并由前一点的相关系数修改估计的初始值位置，以适应视差的变化。依次遍历每个格网点，得到每个格网点的候选匹配点[10,11]。流程如图4，匹配结果如图5。

图4 影像匹配流程图

图5 匹配效果图

图6 纹理映射后的单像对模型

2.4 坐标转换与三维建模

匹配完成后，即根据相对定向、绝对定向参数和匹配结果生成影像点云，对目标进行三维重建。考虑到控制点坐标是通过全站仪自由设站获得的，因此需要将采集点坐标从自由坐标系转换到大地坐标系，转换时只需考虑水平方向的旋转和原点平移。根据实测的2个以上控制点的大地坐标、自由坐标，计算坐标系旋转参数和平移参数，然后根据式（5）将采集点坐标（X，Y，Z）转换到大地坐标（Xt，Yt，Zt），通过插值生成DEM。纹理映射后的单像对模型如图6所示，矿山开采宕口三维模型如图7所示。

图7 分级渲染后的矿山开采宕口三维模型

3 矿山三维可视化分析

3.1 精度评价

为评价精度，精确测量布置的8个检查点坐标，将摄影测量计算的空间坐标与外业实测坐标进行对比统计。使用Canon G5数码相机对南京浦口某矿山进行拍摄，拍摄距离约100 m，对比结果如表1所示。

表1 检查点精度/m

统计检查点坐标中误差为：mx=0．080 m，my=0．059 m，mz=0．075 m。实验表明，使用本文方法采集点云坐标相对精度高于1‰，满足一般矿山精度要求。

3.2 储量分析

矿山三维模型底面积计算根据终采边界拐点的XY平面坐标计算得到，即计算终采边界拐点所构成平面区域的面积，亦为终采范围内模型在水平面上的投影面积。矿山三维模型表面积计算是基于DEM进行的，将DEM每个格网分成两个对角三角形，遍历每个三角形，其面积之和即为格网的表面积，亦为模型的表面积。体积计算时，将矿山三维模型分解成若干个三棱柱体，计算每个三棱柱体在最低开采高程以上的体积，然后对这若干个体积累加求和，即得到模型体积。

储量计算是基于高精度矿山三维模型自动进行的，因此依次建立以下矿山DEM：①将矿山开采前的原始栅格地形图矢量化，插值得到矿区原始DEM0；②将设计好的矿区终采境界图数字化，生成终采境界DEM1；③根据数码影像采集的点云数据生成矿山开采宕口DEM2，获取目前已开采区域R；④根据平面位置将DEM0与DEM2叠加，用区域R内的DEM0数据替换为DEM对应区域的数据，得到矿区当前DEM[12]。

23

储量计算具体步骤如下：①分别计算矿界内DEM0、DEM1、DEM3的体积V0、V1、V3；②计算总浮土层表面积S0，即最低开采高程以上DEM0的表面积；③计算已剥离浮土层表面积S1，即最低开采高程以上、开采区域R内的DEM0的表面积；④总储量T0=（V0-V1-S0×h）×a×t，已采量T1=（V0-V3-S1×h）×a×t，剩余储量T2=T0-T1。其中，a为含矿率，t为矿石比重，h为浮土层厚度。利用每期采集的矿山开采宕口DEM不断更新，可以清晰地反映出矿山储量动态，精确分析超开采情况。

以南京浦口某矿山为例，对开挖宕口共拍摄5个像对，使用免棱镜全站仪测量23个控制点，并通过GPS测量其中2个点的大地坐标作为转换点，外业影像获取及数据采集约45 min，内业数据处理耗时1 h左右。矿区面积102 900 m2，矿石比重2．8 t/m3，含矿率78%，浮土层厚度2 m。储量计算结果如表2所示，对比外业实测计算已采储量2 123 806 t，相对精度达到1．5‰，完全满足矿山要求。

表2 储量计算结果/t

[1] 朱煜峰． 概述我国矿山测量技术的新进展[J]． 中国矿业, 2004, 13(4): 7-8

[2] 段奇三． 徕卡HDS 8800三维激光扫描仪在露天矿中的应用[J]．测绘通报, 2011(12): 79-80

[3] 杨彪, 严丽英, 李浩． 基于普通数码相机的露采矿储量动态检测系统[J]． 工程勘察, 2013(3): 54-59

[4] 陈舒, 李浩, 黄河． 普通数码相机三种检校方法的量测精度评价[J]． 遥感信息, 2012, 27(5): 73-76

[5] Lucchese L, Mitra S K． Using Saddle Points for Subpixel Feature Dectecyion in Camerapore, a Calibration Targets[C]． 2002 Asia Pacific Conference on Circutis and Systems, Singapore, 2002

[6] Lucchese L． Geometric Calibration of Digital Cameras Through Multi-view Rectification[J]． Image and Vision Computing, 2005, 23(5):517-539

[7] 杨立君, 马明栋, 苗立志． 非量测相机近景数字影像相对定向方法研究[J]． 大地测量与地球动力学, 2012, 32(4): 135-138

[8] Tsai R Y, A Versatile Cemera Calibration Technique for Highaccurancy 3D Machine Vision Metrology Using Off-theshelf Cameras and Lenses[J]． IEEE Journal of Robotics and Automation, 1987, 3(4): 323-344

[9] 陈新玺． 多基线普通数码影像近景摄影测量技术研究[D]． 南京: 河海大学, 2006

[10] Zhang L． Automatic Digital Surface Model(DSM) Generation from Linear Array Images[D]． Zurich: Swiss Federal Institute of Technology, 2005

[11] Lowe D G． Distinctive Image Features from Scale-Invariant Interest Points[J]． International Journal of Computer Vision, 2004, 60(2): 91-110

[12] 蒋锐, 宋焕斌． 基于三维栅格数据的露采矿山储量动态监测研究与应用[J]． 矿产与地质, 2009, 23(5): 469-472

P234.1

B

1672-4623（2014）04-0017-04

10.11709/j.issn.1672-4623.2014.04.006

史雨川，硕士，研究方向为数字摄影测量。

2013-08-27。

项目来源：国家自然科学基金资助项目（51079053）。