李运华，牛天昊

（北京航空航天大学自动化科学与电气工程学院，北京100191）

特大型矿用挖掘机在露天矿铲装作业时，由于环境恶劣、作业视场大、能见度差等原因，操作人员很难估计挖掘阻力信息，只能凭经验选择操作功率模式。因此，在挖掘作业开始前通过状态感知和仿真的手段预测挖掘阻力信息，对实现挖掘机最佳作业具有十分重要的意义［1］。

离散元素法是仿真预测挖掘阻力的重要方法，陈烨等［2］利用离散元素法建立挖掘对象的仿真模型，分析挖掘过程中的受力和耗能。在离散元仿真中需要知道矿石粒径参数信息，目前主要是通过人工取样的方式获取，费时费力。黄习敏等［3］对矿石粒径获取进行了研究，利用图像处理技术对矿石图像进行了分割和识别。本文提出利用分水岭算法（Watershed Algorithm，WA）对航拍获得的矿石照片进行分割，从而获得矿石粒径大小，并与离散元仿真相结合。该方法能解决离散元仿真中关键参数获取困难的问题，可以更加准确有效地预测挖掘阻力，为矿山安全、经济作业提供了有效的数据支持。

1 矿石粒径提取

1.1 方案设计

对无人机航拍获取到的矿石图片进行图像处理，获得真实的矿石粒径大小及分布，其实现步骤如下：

步骤1对得到的原始图像进行预处理，得到降噪后的图像；

步骤2对预处理后的图像进行形态学重构，使矿石边缘更加清晰；

步骤3对重构后的图像进行阈值分割，得到二值图像；

步骤4利用WA，可以将图像中每个矿石颗粒分割出来，并计算每一个矿石颗粒包含的像素点个数；

步骤5结合标定过的图像，推算其真实的粒径。

基于图像处理的矿石粒径提取流程图如图1所示。

图1 基于图像处理的矿石粒径提取流程图Fig.1 Flow chart of coal-ore particle size extraction based on image processing

1.2 图像预处理

采用RGB 彩色图像表征矿石信息。为了减少处理工作量，增加背景与矿石颗粒的对比度，需要将彩色图像变换为灰度图像。由于空气中的灰尘以及拍摄时的抖动等原因，图像不可避免地会出现噪声，因此，需要对转换后的灰度图进行滤波降噪。基于实验对比，这里采用中值滤波，如图2所示。

图2 中值滤波后的图像Fig.2 Image after median filtering

1.3 形态学重构

图像进行滤波后，虽然减小了噪声的影响，但是由于矿石本身纹理的存在，使得图像中仍有许多的像素极值区域，对之后的分割造成极大的影响。因此，使用形态学中的重构运算来解决这类问题［4］。首先对图像进行开操作、腐蚀操作后进行第1次重构；其次对操作后的图像进行膨胀，将膨胀后的图像取补，将上一步操作后的图像取补做掩膜，再次重构，并对结果取补。重构操作后不仅能进一步去噪，更能简化矿石纹理，保留矿石原本形状，消除多余的极值区域。

如图3 所示，经过重构后的图像，矿石自身的纹理被消除了，更加突出了矿石原本的轮廓，矿石与背景的明暗对比更加突出，对于后续分割步骤具有重要意义。

图3 形态学重构后的图像Fig.3 Image after morphological reconstruction

1.4 二值化处理

彩色图像二值化是图像分割中最常见的方法，并且在本研究中发现WA 在二值图上的表现更加良好，因此，需要对重构后的图像进行二值化处理，以便后续步骤的进行。本研究采用的是自适应阈值的二值化方法［5］，能尽量减小光照对图像的影响。该算法的原理是将图像视为所有行拼起来的单行，之后遍历图像，计算某一像素点一定范围内的所有像素点的平均值，并将所求得的平均值作为阈值，判断该像素点为黑色还是白色。

假设fs(n)是n处像素点前s个像素的和：

那么n处像素点的值T(n)是置1（黑色）还是0（白色），则取决于该点像素值与其前s个像素的平均值的百分之t的大小关系，计算公式为

对重构后的图像进行自适应阈值的二值化处理，如图4所示。

图4 二值化处理后的图像Fig.4 Image after binarization processing

1.5 分水岭算法

WA 在处理存在黏连、重叠情况的图像时，有着良好的效果，因此，本文选择WA 分割矿石图像。WA 把图像看成测地学上的地貌［6］，图像中每一点像素的灰度值表示该点的海拔高度，通过模拟泛洪过程求取分水岭的分割线。具体过程为：假设在每个区域最小值的位置上打一个洞并且让水以均匀的上升速率从洞中涌出，从低到高淹没整个地形。当处在不同的盆地中的水将要聚合在一起时，修建的大坝将阻止聚合。这些大坝的边界就是分水岭的分割线。

模拟泛洪过程是个递归过程，其初始条件为

式中：Xhmin为图像中灰度值为最小值的像素点。

递归过程为

式中：h、hmin和hmax依次为灰度值、最小灰度值和最大灰度值；Xh+1为灰度值，即海拔高度为h+1 的所有像素点；minh+1为此点属于新产生盆地最小值点，即在h+1 海拔高度又产生了新盆地；Xh∩Xh+1为Xh与Xh+1相交；CXh为Xh所在的盆地。

通过此递归过程，将图像中的所有像素点划分到每个盆地中，最后所有同时属于2 个及以上盆地的相邻像素点所连成的线便是每块矿石的边缘轮廓线。图5为利用WA分割矿石图像的分割线。

图5 分水岭脊线图Fig.5 Segmentation line map obtained by WA

为了更好地展示分割效果，在原图中将分割出的不同区域用不同的颜色覆盖，绘制出伪彩色图，如图6 所示。从伪彩色图中可以看出，大部分的矿石颗粒都被一种颜色完整地覆盖，说明WA 可以较好地完成矿石图像的分割。

图6 伪彩色图展示效果Fig.6 Pseudo-color map

1.6 获得矿石粒径

对分割出来的矿石区域进行连通域标定，读取出每一区域的面积。图像处理中图像的面积单位为像素点的个数，因此，用像素点个数统计矿石面积分布情况，如图7所示。

图7 矿石面积分布Fig.7 Ore area distribution

像素点数与实际尺寸的换算可以采用以下方法：利用无人机飞到相同高度，对标准正方形进行拍摄，标定用的标准方形颜色与背景颜色要有明显差异。获取标准方形在图像中所占像素点个数［7］，用像素点个数与标准方形面积相除即可获得比例系数Z。标准方形如图8所示。

图8 标准方形及其标定图Fig.8 Standard square and its calibration image

若第i个矿石颗粒的面积为Si，则通过简化计算第i个矿石的面积径di为

式中：Di为通过图像得到的第i个矿石颗粒的实际近似粒径。

通过此方法，可以获取矿石图像中每个矿石所占像素点，并完成像素点与实际粒径的转换，获取矿石的粒径大小及分布信息。

2 挖掘阻力仿真

2.1 方案设计

在挖掘机工作装置结构与物料参数已知的基础上，在Edem 仿真软件中，利用离散元素法计算出对于该物料1 次挖掘作业过程中所受挖掘阻力值，并根据阻力值的大小建议合适的操作功率模式，达到既满足作业需求又节省能量的效果，具体流程如图9所示。

图9 挖掘阻力仿真流程图Fig.9 Flow chart of excavation resistance simulation

2.2 离散元素法

离散元素法的基本思路是将仿真对象看成刚性的或可变形的物体，当受到外力时，该物体会获得1 个加速度，在1 个仿真步长内可通过计算得到新的位置；在下一个仿真步长时，物体在新位置又会产生碰撞，根据牛顿第二定律又可以计算出新的外力，进而确定新的位置，如此迭代计算直至仿真结束或平衡［8］。本文所研究的是不可变形、滑动摩擦情况较少的矿石，根据这些特点，在模拟矿山挖掘作业时，选择Hertz-Mindlin（no slip）接触模型［9］，该模型是离散元素法中应用最为广泛的模型，在力的计算方面精确且高效［10］。

2.3 挖掘机工作装置

本文以太原重工WK 系列中的正铲电铲挖掘机WK-35 为研究对象，其工作装置主要由鞍座、推压轴、起重臂、斗杆、顶部滑轮、铲斗组成。挖掘作业时，牵引绳拉动铲斗，使铲斗绕推压轴进行圆周运动。依照WK-35 挖掘机产品手册，在3dsmax 中建立电铲挖掘机工作装置三维模型，将其导入到Edem软件中，并设定其运动姿态变化。

2.4 设置材料参数

颗粒材料相关参数分为两类：一类为颗粒材料自身的参数，即泊松比、密度、剪切模量；另一类为颗粒材料间接触参数，即静摩擦系数、滚动摩擦系数、碰撞恢复系数。在本研究中，作业对象的材料是矿石，挖掘机工作装置的材料为钢铁，岩（矿）石物性资料可获得相关参数值如表1和表2所示。

表1 材料自身参数Tab.1 Material parameters

表2 材料间接触参数Tab.2 Contact parameters between materials

2.5 设置矿石粒径分布

利用第1 章中的技术方案，可以获得矿石粒径大小及分布情况这一关键参数，该参数对于破碎后露天矿矿岩模型的构建至关重要。本研究以抚顺东某露天矿场的数据为参考，得到的数据如表3所示。

表3 矿石粒径与分布表Tab.3 Ore particle size and distribution table

确定上述参数后，在Edem 中生成破碎后的露天矿矿岩模型，如图10所示。

图10 破碎后露天矿矿岩模型Fig.10 Open-pit mine rock model after blasting

2.6 仿真结果及分析

仿真过程以斗杆垂直于地面，铲斗略微插入矿堆中为仿真起始点；挖掘作业时斗杆与铲斗绕推压轴做匀速圆周运动；以铲斗完全离开矿堆，斗杆与地面接近平行为仿真终止点，仿真模拟挖掘动作。仿真模拟挖掘动作如图11所示。

图11 仿真模拟挖掘动作Fig.11 Simulated excavation

仿真结束后将铲斗分为斗与斗齿，进行受力分析。如图12 所示，斗齿在作业时受到多个方向上的力，其受到的合力值较斗受到的合力大很多，且不能体现出牵引绳给予的拉力，因此，本文对斗齿所受的力不展开讨论。斗受到的合力值一直处于合理的区间范围内，并且当铲斗被提升至空中时，合力值稳定。在被挖掘物料的重力附近，斗所受到的合力可以代表牵引绳所提供的拉力。因此，选取斗所受到的合力值作为操作功率模式选择的依据。

图12 斗与斗齿受力情况对比Fig.12 Comparison of force between bucket and bucket teeth

图13为1次挖掘过程中铲斗受到的合力情况，其峰值约为800 kN。根据《WK-35 产品介绍》可知，WK-35 电铲挖掘机最大提升力为2 150 kN。在考虑到安全作业的情况下，操作员可以将输出力定为900 kN，这时建议操作员可以将电机调成中等功率输出。

图13 斗的受力情况Fig.13 The force acted on the bucket

综上所述，在掌握了作业对象的关键参数后，在Edem 中进行离散元仿真，可以准确地获得针对该对象作业时挖掘阻力的大小，并根据挖掘阻力值指导操作员进行操作功率模式选择。

3 结论

（1）本文所提出的通过航拍矿山获取矿石图像后，对图像进行处理获得矿石真实的粒径大小的方案，可以有效解决在挖掘阻力仿真中，矿石粒径大小及分布信息获取困难的问题。

（2）本文为操作员提供了感知作业阻力的有效方法。利用已知的作业对象参数，在Edem 软件在构建模型，利用离散元素法进行仿真，获得挖掘阻力信息，为操作员选择操作功率模式提供了理论依据。