基于倾斜摄影测量的露天矿数值模型重构与计算实例

刘阳晓，朱万成，刘文胜，刘溪鸽，王江梅

( 1. 东北大学 资源与土木工程学院，辽宁 沈阳 110819；2. 鞍钢集团矿业有限公司，辽宁 鞍山 114001 )

数值模拟作为常用的露天矿稳定性分析方法，其模型是否精准直接影响计算结果的准确性。传统的数值计算模型，是由地质工作者现场实地采集有限的地形数据建立的。然而，由于大型露天矿山地形、地表构造较为复杂，与剥离开采、降水及植被生长等多因素关联，以实际采集有限离散数据构建简化地质模型的方法亟需改进。目前，随着工业无人机技术的日益普及，倾斜摄影测量及配套三维实景建模软件发展迅速，人们可以快速高效地获得视觉良好的高精度露天矿三维实景模型。然而，如何进一步通过三维实景模型快速重构露天矿数值计算模型，为高效准确的数值建模提供新方法，更具现实意义。

倾斜摄影测量技术通过同一飞行平台搭载多镜头，从垂直、倾斜多角度采集影像数据，利用三维实景建模软件经过空中三角测量、影像密集匹配、构建TIN模型和纹理映射后[1]生成精度为厘米级[2-3]的三维实景模型。利用无人机灵活、安全、精度高、周期短的特点，倾斜摄影测量在露天矿山的地质调查[4-5]、资源监管[6]、运输监测[7]、工程量计算[8-9]和边坡变形监测[10-11]方面均有研究和应用。针对通过倾斜摄影三维实景模型重构数值计算模型的问题，国内外学者也有一定的研究。高相波[12]等对具有地质灾害隐患的危岩体进行高精度三维激光扫描获得点云模型后，对边坡发育的植被部分点云进行降噪、网格清理等人工处理后，构造边坡NURBS曲面模型，实体化并划分网格后可导入FLAC3D软件；金爱兵[13]等选取植被覆盖较少的边坡进行摄影测量得到数字高程模型，通过Surfer、Rhino软件转化为网格模型，实现无人机测量到三维数值模型的快速建立；邓小龙[14]等利用逆向工程软件，由点云模型重构三维地质模型，用Hypermesh划分网格后用于数值计算。总体而言，国内外学者的研究多集中在倾斜摄影测量精度、模型建立及露天矿监测方面，针对倾斜摄影模型重构数值计算模型的研究，多选取没有植被干扰的边坡进行摄影测量，或对采集到的数据进行人工去噪后，利用软件完成倾斜摄影模型到数值计算模型的转换。

露天矿倾斜摄影测量，受植被、矿山设备和人员等的影响，建立的三维实景模型包含大量干扰 物，致使在进一步的模型重构与网格划分中频繁出错，无法直接建立有限元数值计算模型。并且，对倾斜摄影测量模型人工去噪又极为费时费力，模型重构效率低下。针对上述问题，本文采用数字图像处理的方法清除干扰物，生成规则化等高线，再构建实体模型并划分网格，可高效完成由三维实景模型重构数值计算模型。

1 模型重构原理及网格划分

本文数值计算模型重构方法分为两部分，第1部分是编程实现的等高线提取及基于数字图像处理的等高线规则化方法；第2部分是借助第三方软件实现的由规则化等高线建立实体模型并划分网格，生成网格模型用于数值计算。数值计算模型重构流程如图1所示。

图1 数值计算模型重构流程 Fig. 1 Reconstruction flow of numerical calculation model

1.1 等高线提取及规则化实现

1.1.1 模型预处理及等高线提取

由倾斜摄影采集数据建立的三维实景模型，多用于在场景中展示而无法导入其他三维建模软 件。为解决三维实景模型难以二次开发的问题，本文通过遍历模型三角网，提取三角网模型节点三维坐标构建点云模型。点云模型为已知坐标点的集合，模型可保存为文本，易于二次开发。

本文采用直通滤波的方法在点云模型中提取等高点云。直通滤波是对点云模型指定维度指定区间进行剪裁，保留区间内的点。将滤波维度设置在露天矿深度方向，将滤波区间设置在较小范围，则直通滤波得到的点近似为等高点。

1.1.2 等高线规则化方法

受矿山情况复杂及实景模型存在离散三角面的影响，由点云模型提取的等高点存在噪声点而无法直接导入三维建模软件实现实体模型重构及网格划分。本文将三维点云映射至二维图像[15]，利用数字图像处理方法剔除噪声点，规则化等高线。采用点栅格化的方法将三维点云映射至二维图像，以点云的点平均距离为像素点边长，构造值为0的图像矩阵。遍历点云，利用式( 1 )确定三维点在二维图像矩阵中的位置，并将矩阵值设置为1，即得到三维点云对应的二维图像。

式中，round( )为四舍五入运算；(X，Y)为点云坐标； (I，J)为点在图像矩阵中的坐标；L为图像像素边 长。

将三维点云映射至二维图像后，通过填充图像孔洞并统计连通区域面积，识别筛除噪声点。剔除噪声点后通过检测孔洞边缘获得规则化后的等高点，将等高点按逆时针顺序连接即得到规则化等高线。

( 1 ) 形态学开操作与闭操作。形态学操作利用结构元在图像中遍历，提取特定图像分量。形态 学[16]的基本操作包括膨胀和腐蚀，膨胀操作使图像变粗，腐蚀使图像变细。膨胀和腐蚀的定义分别为

式中，A为图像；B为结构元；为结构元相对原点的反射；z为结构元的原点；∅为空集。

由式( 2 )可知，结构元B对图像A膨胀的结果是结构元B的反射与图像A有重叠时z的集合。结构元B对图像A进行腐蚀的结果是结构元B中包含于图像A中时z的集合。膨胀和腐蚀的组合使用构成开操作和闭操作，开操作和闭操作的定义为

开操作在图像中的作用是断开狭窄的间断和消除细的突出物，闭操作则是为了填补轮廓线中的断裂和消除小的孔洞。在本文研究中，将点云图像填充孔洞前，用闭操作连接细的断开部分，以避免填充溢出。在填充孔洞后，统计连通区域前，用开操作断开细的连接以保证区分等高点云与噪声点部分。利用形态学闭操作和开操作，使规则化等高点的程序更加健壮。

( 2 ) 填充孔洞及连通分量统计。填充孔洞在本文有2个作用，一是对原点云图像进行边界提取，会提取到2条边缘线，填充孔洞后可以只提取到1条边缘线；二是填充孔洞后可以根据连通区域的面积剔除离群点。孔洞填充以孔洞边界为约束条件，以孔洞内1点为种子点，通过对孔洞循环膨胀操作直至图像稳定，孔洞全部填充。孔洞填充的公式为

式中，Xk为填充孔洞后的图像；为图像补集；k＝0时，X0为初始种子点图像。

对点云图像填充孔洞后，由于噪声点与等高点孔洞属于不同连通分量，通过统计各个连通分量像素点个数，可以区分等高点和噪声点，剔除噪声 点。

( 3 ) Sobel边缘检测。边缘检测最常用的方法是检测亮度的不连续性，其由梯度的幅值衡量。边缘检测的准则是在图像中寻找梯度的幅值大于阈值的地方。梯度的幅值为

Sobel算子[17]采用水平和垂直2个方向的卷积模板与图像进行卷积运算，以3×3邻域的行和列之间的离散差计算梯度。图像3×3邻域示意图如图2( a )所示，Sobel算子水平和垂直梯度模板如图2( b )，( c )所示。

图2 3×3图像邻域及Sobel算子示意 Fig. 2 Neighborhood of image and Sobel operator

Sobel算子2个方向卷积模板与图像进行卷积运算后，可得水平和垂直方向一阶导数为

式中，Z为亮度。

由Sobel算子卷积得到水平和垂直梯度后，可得Z5处梯度公式为

如果∇f大于等于阈值T，则该位置像素为边缘像素。

( 4 ) 改进的区域生长算法连接无序点成线。传统区域生长算法是基于灰度相似性的图像分割方法，其基本思想是将具有相似性质的像素集合起来构成区域。具体地，确定种子点后，将在种子点像素周围满足相似规则的种子像素合并到区域中，直到没有满足条件的像素点可被包括，则1个区域生长完成，寻找下一个像素点进行生长[18]。本文在种子选取方法上进行改进，传统的种子随机生成，本文中种子的选取，是在1个区域生长的最后1个点按照8邻域扩散的方式查找，由此可以保证2个邻域首尾按顺序连接。

1.2 实体模型构建及网格划分

Rhino作为专业的三维建模软件，可快速实现由线生成面及由面重构实体模型。本文由实体模型提取等高线后，通过对等高线的规则化处理，很大程度上避免了在Rhino软件计算中的布尔运算错误和划分网格时存在自相交面的错误。通过建模软件Rhino可以由等高线生成露天矿表面模型，由表面模型边框曲线挤出闭合曲面可快速建立三维实体模型。对实体模型利用Griddle插件划分网格后，即生成可导入数值模拟的网格模型。

2 模型重构方法应用实例

2.1 三维实景模型建立

大孤山矿区总面积176.98 hm2，开采深度至-350 m。采用科比特M6无人机搭载5镜头从垂直方向及前后左右与地面45°方向进行一定重叠率的拍摄，采集得到图像共31 110张及对应无人机POS( Positioning and Orientation System定位定姿系统 )数据。外业获得倾斜摄影的影像数据、POS信息和像控点坐标后，导入三维建模软件ContextCapture进行空中三角测量，多视角影像密集匹配，构建TIN模型和3D纹理映射后，生成高精度的三维场景模型，如图3所示。

图3 大孤山露天矿三维实景模型 Fig. 3 Real scene terrain model of Dagushan open-pit mine

2.2 等高线提取及规则化处理

2.2.1 大孤山模型预处理及等高点提取

通过分析三维实景模型的数据结构，利用三维渲染引擎OSG库的节点访问器NodeVisitor遍历实景模型的三角面节点，提取节点三维坐标构建点云模型。大孤山点云模型点数量8 000余万，为减少运算成本，提高计算效率，采用体素化网格的方法对模型下采样，在保持点云形态不变的情况下，将模型点数量由8 000余万减少至532余万，下采样后，大孤山露天矿点云模型如图4所示。

图4 大孤山露天矿点云模型 Fig. 4 Point cloud model of Dagushan open-pit mine

大孤山露天矿深度方向即为Z轴方向，最小值为-351.92 m，最大值为50 m。点平均间距为 0.95 m。以Z轴为滤波字段，滤波区间设置1.5 m，相邻滤波区间设置为3 m，对模型直通滤波处理后获得等高点云模型，如图5所示。

图5 等高点云示意 Fig. 5 Contour line of point cloud model

2.2.2 等高点规则化

受点云模型中冗余点与离群点的影响，未规则化的等高线在Rhino®软件中构建实体模型出现布尔运算错误及网格划分存在自相交面错误，无法完成网格模型重构。本文通过点栅格化的方法，将三维点云映射至二维图像后，通过数字图像处理的方法规则化等高线。

在Beijing1954坐标系下，大孤山点云模型X范围[-608.36，1 192.35]，Y范围[-792.4，858.43]。以大孤山模型-269 m高程等高点( 图6( a ) )为例，构建大小为1 892×1 735的全0图像矩阵，通过式( 1 )，令L＝0.95，Xmin＝-608.36，Ymin＝-792.4，将三维点云映射至二维图像，如图6( b )所示。对于等高点图像，填充孔洞后( 图6( c ) )，规则化等高点与噪声点区域不连 通，属于不同连通分量。通过统计图像中连通区域分布，生成连通分量标记矩阵( 图6( d ) )来区分噪声 点。

-269 m点云图像，连通分量标记结果如图6( d )所示。点云图像共包含8个连通分量，等高点构成的连通分量( 图6( d )粉色区域 )面积最大，7个噪声点构成的连通分量( 图6( d )红、蓝、黄、绿、粉红、紫、青区域 )面积小。统计各连通分量像素点个数，区分并剔除噪声点，得到无噪声点云图像( 图6( e ) )。经过边缘检测，提取边缘点即得到规则化等高点图像( 图 6( f ) )。为检验规则化等高线的精度，将规则化等高线与等高点云对比( 图6( g ) )可以看出，规则化等高线保持原点云模型形态。将无序的等高点规则化处理后，通过改进的区域生长算法，将点按逆时针方向连接为线即为规则化等高线( 图7 )。

图6 -269 m等高线规则化处理流程 Fig. 6 Flow diagram of regulating contour line at -269 m

图7 大孤山露天矿规则化等高线 Fig. 7 Regularized contour line of Dagushan open-pit mine

2.3 实体模型建立重构及网格划分

规则化的等高线，利用Rhino®软件，可快速由线生成矿山壳体模型，拉伸生成封闭模型后，利用Griddle插件划分网格后，生成FLAC3D格式的网格模型，如图8所示。模型中面体单元数为687 052，节点数为132 713。

图8 大孤山露天矿三维网格模型 Fig. 8 Three dimensional grid model of Dagushan open-pit mine

3 模型精度验证及数值计算实例

3.1 模型精度验证

完成数值计算模型的重构后，对模型精度 进行分析。以三维实景模型为准确地形模型参 照，对比本文重构网格模型与传统方法建立的网 格模型精度。分别对三维实景模型、重构网格模 型及传统网格模型在X与Y方向的1/4，1/2，3/4处 取6条剖线，如图9( a )所示，剖线5的对比如图9( b )所示。

图9 大孤山露天矿代表性剖面选取与对比 Fig. 9 Selection and comparison of representative sections in Dagushan open-pit mine

每个模型获得6个剖面的高程向量，通过计算网格模型与实体模型剖面高程向量的欧式距离，判断剖面相似性。欧式距离定义为

式中，z1，z2为剖面高程向量；dist(z1，z2)为2向量的欧式距离。

欧式距离越小，代表剖面相似度越高，距离大则剖面相似度低。令z1为三维实景模型剖面高程向量，令z2分别为本文重构网格模型和传统网格模型的高程向量，计算模型剖面间欧式距离( 表1 )。由表1可知，本文重构模型与实景模型剖面欧式距离小于传统网格模型与实景模型剖面的欧式距离，即本文重构模型与三维实景模型剖面更加相似。这是因为，传统数值计算模型的建立由于人工采集地形数据有限而被大大简化，与实际矿山地形不符。本文基于无人机采集的大量矿山地形数据重构网格模型与实际矿山地形更加符合，准确度高于传统数值计算模型。

表1 模型剖面欧式距离 Table 1 Euclidean distance of profiles

3.2 数值计算

为验证重构的网格模型在数值模拟中的三维应用，将模型导入FLAC3D进行力学分析计算。计算模型大小为2 000 m×2 000 m×890 m。数值模拟采用线弹性本构模型，根据对露天矿取岩样后进行室内物理力学试验结果，采用计算模型参数为：岩体体积模量(K)为4.44 GPa；剪切模量(G)为3.33 GPa；密度(ρ)为2 300 kg/m3；重力加速度为9.8 m/s2。在模型前后边界施加沿X方向固定水平位移约束，左右边界施加沿Y方向的固定水平位移约束，底部施加沿Z方向的固定垂直位移约束。通过迭代计算，当体系最大不平衡力与典型内力的比率小于给定值( 1.0×10-5)时计算收敛，得到在重力作用下生成的初始地应力场。其中，垂直地应力分量计算结果如图10所示。

图10 大孤山露天矿垂直地应力分布规律 Fig. 10 Distribution of vertical stress in Dagushan open-pit mine

沿该露天矿最深标高处的点以下、岩体内部一点P为研究对象，其中，此处P恰取为最深标高点与整体模型底部距离的中点。通过数值计算后处理，得到该点的垂直地应力分量σv＝5.0 MPa；而通过岩体容重( 或密度 )进行地应力理论计算：σv＝ρgh＝2 300×9.8×222×10-6＝5.004 MPa；二者之间的相对误差不超过0.8%。

4 结 论

( 1 ) 提出了针对三维实景模型的等高线提取及规则化方法。采用直通滤波的方法提取等高点，通过数字图像处理的方法剔除噪声点及冗余点，生成规则化等高线。

( 2 ) 借助第三方软件，由规则化等高线可快速实现实体模型重构及网格划分，生成数值计算模型，提高数值建模效率。

( 3 ) 通过与实测的三维地质模型进行剖面的对比分析，发现本文方法相比传统数值建模，模型更加符合露天矿实际地形。