韩 岭,童 星,李俊涛 ,杨 彪
(1.湖北省地质局第八地质大队,湖北 襄阳 441000;2.河海大学 地学院,南京 211100)
露天矿产资源是重要的自然资源,对保障国家基础设施建设、国防建设以及城市化进程起着至关重要的作用[1]。为促进矿产储量的合理开发和有效保护,兼顾环境保护与可持续发展,国家对露天矿实行动态监管,加强了对矿山储量的动态监测,要求详细掌握矿山的年开采量、损失量、资源储量变化量。我国露天矿数量多、分布广,储量动态监测工作量巨大,长期盛行的全站仪或RTK地面接触式监测方法存在测点难以到达、安全隐患多、成果直观性差等不足,且这种点的数据采集方式已不能满足矿山动态监测的精度和效率要求[2]。
近年来,无人机(UAV)航测技术以其高效的面域数据采集方式,为解除矿山监测人员安全风险、提高监测效率和精度带来了良好契机。为此,国内外对露天矿无人机动态监测技术开展了大量实践研究。许志华等[3]提出了一种基于无人机影像序列重建和DTM三角网差值法的露天矿采剥量、堆放量计算方法;ESPOSITO等[4]研究了基于无人机三维点云和正射影像的意大利撒丁岛露天矿多时相的开采范围和体积变化评估方法;王果[5]、陈凯[6]、白洋[7]、谢君洋等[8]研究了基于无人机影像的矿山动态监测方法,通过两期无人机DSM叠加,分别计算矿山占地面积、填挖方、剥离量和爆堆土方量;马国超等[9]将无人机倾斜摄影与LiDAR相结合开展了露天采场安全监测融合应用研究,实现了多尺度精准监测采场安全隐患;黄立鑫等[10]利用无人机DSM研究了尾矿库洪水淹没分析和库容计算方法。然而,现有研究仅限于利用无人机航测技术进行矿山土石剥离方量估算,当矿山存在多个矿种或存在浮土层、非矿夹层时,土石方量并不能作为储量监测最终结果,还须紧密结合地质模型、终采设计、浮土层分布、矿界和采深限值等多约束条件,智能分析每类矿体的界内合法开采量与越界越层等非法开采量,过程非常复杂。因此,目前行业内储量估算还普遍依靠人工平行断面法。该方法仅利用无人机数据生成的高精度断面数据参与计算,断面之间的高精度地形数据不能有效利用。当矿界形状不规则或地形起伏较大时,将造成严重的精度损失,产生较大误差。综上,无人机地形数据、矿山地质数据等多源异构数据融合,以及储量智能计算方法亟待进一步研究解决。
本文研究矿山无人机数据与矿山地质数据的无缝集成建模与储量精细化计算方法,以充分发挥无人机技术优势,实现智能精准的露天矿储量动态监测,并为矿山智能监测技术研究与开发提供理论参考与技术框架。
本方法总体技术路线如图1所示。
图1 总体技术路线Fig.1 Basic technical flow
基于无人机航测,获取矿山各监测期的高精度地形数据,精细重建矿山多时相三维地形模型,是储量智能化动态监测的关键技术之一。
1.1.1 初始地形与现采地形三维建模
矿山多时相地形包括初始地形、各期现采地形和终采设计地形。其中,矿山初始地形是储量计算的基准。鉴于无人机摄影测量与无人机LiDAR的互补优势(如表1),矿山初始地形建模时,低植被覆盖度的矿山宜采用无人机摄影测量技术,植被覆盖度较高的矿山则采用无人机LiDAR技术。无人机数据处理时,须滤除植被才能精确还原真实地表。现有的植被滤除处理方法有两种。第一种是先生成矿区DSM,再采用DSM地形滤波算法得到矿山DEM;第二种是先对无人机三维点云数据进行地形滤波,再基于滤波后点云建立DEM。第一种方法虽然快速高效,但会造成局部真实地形数据丢失,当植被覆盖度高、地形复杂时,将产生较大的误差[11]。第二种方法直接对全部点云进行滤波,可靠性较高。因此本文采用第二种方法。
表1 无人机摄影测量与LiDAR的互补优势
无人机三维点云数据海量,滤波后密度不均匀,兼顾点云插值DEM的精度、效率和计算机内存消耗,采用基于快速排序的逐点插入法构建Delaunay三角网,再采用三角网拓扑关系引导的插值点快速定位算法和三角形插值算法,实现海量点云的DEM快速生成,如图2所示。
图2 矿山三维点云滤波与地形建模Fig.2 3D point cloud filtering and modeling
初始模型建立后,后续动态监测仅需在已开采范围内进行,未开采区域的地形可从初始DEM中获取,无须考虑植被遮盖影响。因此,无需采用无人机LiDAR,仅采用无人机摄影测量生成实景三维模型和DSM、DOM,以降低数据获取代价,监测成果也更加丰富、直观。当矿区存在稳定特征点时,还可从前期实景三维模型中提取像控点,免除现场像控点测设工作,进一步提高数据获取效率。
1.1.2 终采设计地形三维建模
终采设计图规定了矿山的合法开采范围和最终开采形态。建立矿山终采模型DEM,为越界、越层开采定量分析提供判据。建模方法如下:
1)提取终采设计图中各边坡、平台、平面的轮廓线;
2)根据轮廓线已知顶点高程或设计坡比参数,插值计算全部节点的高程;
3)对轮廓线进行三角化,得到三角网;
4)通过三角网插值,得到终采设计DEM(图3)。
图3 终采设计地形三维建模Fig.3 3D terrain modeling of final mining design
不同于常规的土石方计算,矿山常存在多矿种、非矿夹层、浮土层,需紧密结合地质资料才能智能分析各矿种的动态储量。针对常规地质体建模需要依赖大量人工操作来完成,研究两种常见地质资料的矿体智能建模方法,以实现各矿种储量分类计算。
1.2.1 DEM与地质平面图融合的矿体智能建模
地质平面图的主要内容包括二维地质界线和地质产状标注,矿体智能建模流程如下:
1)DEM叠加:将地质图与DEM叠加,如图4(a)所示,并根据DEM插值地质界线各节点的高程;
2)产状分配:按邻接关系将产状信息赋予邻近地质界线的邻近节点,进而根据已知产状的节点插值未知节点的产状;
3)地质界线偏移:根据产状和地质体底面高程,计算地质界线节点的偏移点坐标,再将偏移点连线,经自相交检测和退化处理,得到底板地质界线;
4)曲面生成:将顶面和底面地质界线分别进行三角化,构建地质体顶面和底面三角网;将顶面与底面地质界线对应节点按次序连接成三角网,构成地质体侧面;基于DEM,对顶面三角形进行高程插值和递归剖分,直到每个三角形内高程起伏小于设定阈值。三角形剖分方法如图5所示,每次沿最长边上的中线将三角形一分为二。
5)拓扑重构:对各地质体进行孤岛与孔洞检测,并重构地质体拓扑关系。最终,形成一个由三角网曲面构成的地质体,如图4(b)所示。
图4 地形地质平面图建模Fig.4 Topographic and geological plan map modeling
图5 三角形递归剖分方法Fig.5 Triangle recursive subdivision method
1.2.2 DEM与地质剖面图融合的矿体智能建模
平面图建模方法适合呈水平分布的单层地质体建模,对于包含复杂曲面夹层的垂直分布多层地质体,难以准确建模。此时地质剖面图建模方法具有互补优势。方法基本流程如下:
1)空间还原:根据起点和终点地面坐标,将二维剖面线还原到真实空间位置,如图6(a)所示;
图6 地质剖面图建模Fig.6 Geological profile modeling
2)剖面分割:逐剖面线提取特征点,并结合至高点和最低点对剖面线进行分割;
3)剖面匹配:根据相邻剖面线重心之间的距离、形状参数建立相似性测度,对相邻剖面中的剖面线进行正反匹配,确定剖面线之间的空间连接关系;
4)剖面连接:按照分割结果,将匹配成功的剖面线连接成空间三角网,得到地质体的表面模型;
5)表面剖分:对表面模型各三角形进行判断,若为裸露地表,则基于DEM对三角形进行高程插值和递归剖分,直到每个三角形内高程起伏小于设定阈值。最终生成的三维地质体模型如图6(b)所示。
考虑矿山浮土层厚度变化通常不大,采用全局或局部多项式曲面拟合算法,建立整个矿山的浮土层厚度估算模型,修正地质体地表高程。
动态储量智能化计算,即自动计算分析不同矿种的界内、越界、越层、边坡开采量,如图7所示,图中标注的动用储量类型说明见表2。
图7 精细储量计算剖面示意图Fig.7 Fine reserve calculation profile
表2 精细储量类型
设初始模型、上年度模型、现采模型分别记为DEM0,DSM1,DSM2。对指定计算范围内DSM2逐格网进行体积计算,方法如下:
1)棱柱构建:取DSM2单元格网四个角点作为四棱柱底面;根据角点平面坐标分别从DEM0、DSM1中插值出初始地表高程h0和上期地表高程h1;为消除DSM2中植被影响,取h0,h1中最小值建立四棱柱顶面;为消除已采矿石堆体影响,当棱柱上表面高程低于下表面时,则将其视为无效棱柱,不参与后续计算;
2)分别用最低开采标高平面、各矿体表面、终采DEM,对四棱柱进行求交,若交点存在,则对四棱柱进行截取分段,并判断每段对应的矿种和储量类型;
3)将四棱柱沿对角线分割成两个三棱柱进行分段体积计算;
4)按上述步骤,对计算范围内的四棱柱体积分类累加,并引入矿界平面范围,区分界内与越界开采,最终得到表1所示的精细储量;
5)引入矿石比重参数、总储量和历史储量记录,统计各矿种的总开采与保有储量。
由于无人机系统重建的DEM、DSM格网大小为cm级,因此上述算法计算结果接近于三维动态储量的严密积分值。
基于上述框架,通过VC++面向对象编程,底层开发了一套《OMRIS无人机露天矿储量动态监测系统》(如图8),实现了无人机LiDAR与摄影测量等多源异构数据融合、多时相矿山三维数字孪生、动态储量精细计算、监测成果输出以及监测信息管理一体化,已成功应用于湖北省数十座露天矿动态储量监测。以湖北省襄阳市某露天石料矿储量动态监测为例(图8),分析应用效果。
图8 无人机实景三维模型与地质体模型叠加Fig.8 Superposition of UAV real scene 3D model and geological model
该矿山平面面积约5.3万m2,采深范围260~340 m,资源总量273.8万m3。系统基于1∶1 000地形地质图、终采设计图分别建立了矿山初始DEM、终采DEM和矿体模型,并通过无人机摄影测量建立了2020、2021年两期DSM,将系统计算得出的2021年各矿种动用储量与专业地质人员采用传统平行断面法计算结果进行对比,如表3所示。计算结果表明,随着所选控制断面数量的增加,平行断面法计算结果逐步向本系统计算结果趋近。由平行断面法原理可知,随着断面数增加,其计算结果将逐步趋近真值。可见,本系统计算结果准确可靠。
表3 本文方法与平行断面法计算结果对比
效率方面,当使用6条以上剖面时,平行剖面法各种图件制作须耗时24 h以上,储量计算耗时1 h以上,主要依靠人工完成;而本系统在普通PC机上运行,数据处理总耗时小于15 min,其中储量计算采用多线程技术优化后耗时不超过30 s。
构建了一套较完整的无人机露天矿动态储量监测方法,并给出了该方法的底层实现框架。应用结果表明,本文方法切实可行、精度可靠,与现有方法相比,本文方法具有显著优势:
1)有效融合了矿山地质与无人机数据,能适应地质构造复杂、地形起伏大、矿界极不规则、植被密集覆盖等各种情形,实现了各矿体的界内、越界、越层、边坡动态储量的精细化、智能化计算分析;
2)将数字孪生技术与矿山动态监测业务深度融合,通过三维可视化呈现,克服二维图件抽象复杂的局限性,形成对矿山开采动态的直观、准确、全面感知;
3)数据处理智能化,人工干预环节少,不受计算人员主观因素影响。