机器视觉技术估算边坡连续形变体积

叶春阳 许传华 孙国权 聂 闻，4

（1.中国科学院海西研究院泉州装备制造研究所，福建泉州362000；2.中北大学电气与控制工程学院，山西太原030000；3.中钢集团马鞍山矿山研究总院股份有限公司，安徽马鞍山243000；4.江西理工大学资源与环境工程学院，江西赣州341000）

采矿区域的不断扩张，伴随而来的是裸露边坡数量的增加，矿山边坡一旦发生形变失稳将会给矿山生产以及人员生命安全带来巨大的威胁［1-4］。因此通过采取一定的技术措施降低事故发生率对矿山生产安全具有重要意义［5-6］。边坡形变识别方法按照监测形式可分为接触式监测和非接触式监测两类。以接触式监测为代表的边坡形变识别方法主要是在边坡中埋设各种传感器设备，然后通过传感器数据来识别边坡形变区域。这种方式主要存在两方面的不足：一是需要耗费大量的人力物力，且无法对设备安装人员生命安全提供保障；二是获取的信息有限，通常一个传感器只能反映一个点的变化情况，因此这种方法获取的数据难以反映边坡整体变化情况。机器视觉是一种近几年来发展较快且比较成熟的技术，它作为一种非接触式的技术能够有效减少人为参与并提升生产效率。通常机器视觉所用的监测仪器主要有深度相机、雷达等，许多学者利用这些设备提出了许多经典的边坡数据处理方法。例如YANG等［7］利用超像素算法处理裁剪的边坡视频数据，获取了边坡形变区域的轮廓；MWANIKI等［8］利用图像增强方法处理图像数据提升边坡的识别精度；OHNISHI等［9］研究了数字相机在边坡形变位移上的应用，实现了用数字相机测量边坡位移量；刘昌军等［10］利用激光数据实现了高边坡危岩体的识别；李洪梁等［11］利用三维激光扫描技术实现了堰塞湖的应急测绘。以上方法剔除了作业人员与边坡的近距离接触环节，在一定程度上避免了人员受到生命安全威胁；此外，相机和雷达获取的数据范围更广、维度更高，从而更能反映边坡的整体趋势。

本研究通过室内模型试验，首先利用深度相机实时获取边坡形变图像数据和点云数据；然后分析相机坐标系、像素坐标系、世界坐标系、图像坐标系之间的转换关系，通过相机标定最终建立图像与点云之间的联系；最后利用背景差分法实现边坡形变区域的识别，通过前后时刻的边坡表面位移量估算边坡形变的体积。本研究旨在探索边坡形变二维数据与三维数据相结合实现边坡连续形变区域的体积估算，为边坡形变特征信息的自动提取提供参考。

1 相机测量原理

本研究分别建立世界坐标系、相机坐标系（点云坐标系）、图像坐标系以及像素坐标系，并对相机进行标定，确定点云坐标与像素坐标之间的参数关系。

1.1 世界坐标系

世界坐标系用于描述相机的空间位置，真实反映相机在三维空间的坐标关系。如图1所示，Ow、Xw、Yw、Zw分别代表世界坐标系的原点、X轴、Y轴和Z轴。

1.2 相机坐标系

相机坐标系是联系世界坐标系与图像坐标系的过渡坐标系，是空间数据过渡到平面坐标之间的桥梁。从世界坐标系变换到相机坐标系属于刚体变换，即物体不会发生形变，只需要进行平移和旋转。通常以 Oc、Xc、Yc、Zc表示相机坐标系，本研究使用的点云数据所在坐标系即为相机坐标。

1.3 图像坐标系

图像坐标系是相机坐标系沿Z轴投影而来的，常用O、X、Y表示图像坐标系的原点、X轴和Y轴，如图2所示。

如图2所示，坐标系原点O与相机坐标系原点Oc之间的距离为f，其中ΔABOc与ΔONOc相似，ΔPcBOc与ΔPNOc相似，则可得以下关系式：

进而求得相机坐标系与图像坐标系的关系为

1.4 像素坐标系

像素坐标系由图像坐标系沿X轴负方向和Y轴负方向分别平移u0和v0个单位得到，常用O0、U、V来分别表示像素坐标系的原点、X轴和Y轴。如图3所示，若矩形O0QWE为一幅图像，则点(u0,v0)为该图像的中心。

由图3可得图像坐标系与像素坐标系转换关系：

式中dx，dy分别表示图像坐标系下沿X轴和Y轴的长度系数。

1.5 相机内参计算

由相机坐标系、图像坐标系以及像素坐标系三者之间的转换关系，可以推导出相机坐标系与像素坐标系有以下关系：

式中，a，b表示相机内参。

如图4所示，选取15张不同拍摄角度下的标定板照片，根据张正友［12］标定法对相机标定。相机内参计算得：a=624.19，b=635.99。

2 边坡形变区域识别

2.1 缺失值补充

光照不均会导致部分点云数据缺失。数据的缺失一方面会影响数据的完整性，另一方面可能导致计算结果与真实结果相差较大。插值技术能够有效解决因数据缺失造成的误差较大问题。因此本研究采用最近邻插值算法［13-14］估计缺失数据。最近邻插值算法的原理是根据距离缺失值最近的k个点的值估算缺失值，其定义如下：

式中，P为缺失值的估计值；pi为缺失值附近第i点的值。

本研究首先由式（4）确定点云数据中缺失值的位置，然后由式（5）对缺失值进行插值估算，其中k设置为5。

2.2 坡体提取

将坡体数据与背景数据分离有利于坡体中形变区域的识别。因此在进行边坡形变区域识别前需要对数据进行预处理。首先在点云3个坐标轴方向设置阈值，获取粗提取的边坡数据；然后利用坡体与背景的颜色差异对坡体进行精细提取；最后保留的数据即为最终的坡体数据。依次处理每个时刻的边坡点云数据以获得完整时间序列下的坡体数据。

2.3 边坡形变区域识别

边坡发生形变时会有两个特征发生显著性变化：一是边坡形变区域颜色会发生改变；二是边坡表面深度会发生改变。合理地利用边坡形变时的两个特征能够有效地识别边坡形变区域。本研究首先用颜色特征处理边坡数据获取带有噪声的边坡形变区域；然后利用边坡表面深度特征细化上述结果，获取更加精细的边坡形变区域。

在外界光源稳定的情况下，边坡形变区域颜色特征不会发生变化或者发生微小变化，而形变区域颜色变化明显。基于以上信息，本研究采用背景差分算法［15-17］剥离边坡非形变区域数据，得到带有噪声的边坡形变区域数据。

利用边坡深度特征细化带有噪声的边坡形变区域。随着时间推移，边坡表面某一区域并未发生形变，但是该区域在后续时刻的深度会发生细微变化。造成这种细微变化原因是边坡表面土壤颗粒发生轻微移动。为防止边坡表面土壤颗粒发生轻微移动的区域被误识别为形变区域，设置一个阈值来避免该种情况发生。阈值选取规则为：选取初始时刻处理的点云数据，记该时刻点云平均深度为h0；选取视觉可见的坡体发生形变的第一张点云数据，记该时刻点云平均深度为h；变化量d=h-h0则为阈值。当边坡点云变化量大于d，则视为发生形变；反之，则未发生形变。经过上述两个步骤处理后的区域即为最终识别的边坡形变区域。

3 试验验证

3.1 物理模型试验

为获取边坡形变数据，在室内进行模拟降雨导致边坡形变破坏的试验，并用深度相机对边坡形变过程进行拍摄记录。图5为试验现场图，试验模型箱长1.6 m、宽0.8 m、高0.8 m。试验共分为两组，A组降雨为65 mm/h，B组降雨为30 mm/h，两组坡角均为45°。

3.2 结果分析

全部程序由C++编程语言编写。如图6所示，第1列为原始点云，第2列为粗提取的边坡形变区域结果，第3列为粗提取后细化的边坡形变区域结果。从第1列数据可以看出边坡发生形变后，坡体向内坍塌，此时该形变区域的Z坐标将发生变化。第2列数据显示基于阈值分割后坡体上部未发生显著形变的部分区域仍然被标记为形变区域，造成该种情况的原因是模拟的降雨影响了相机对该区域坐标的计算。噪声剔除后，如第3列所示，可以看见优化后的区域只剩下坡底被标记的形变区域。

为了定量分析本研究方法的可靠性，以目标分割系数Dice［18］作为评价指标定量分析所提方法的性能。Dice是描述识别结果与原始破坏区域的重合一致性程度，其公式为

其中，G表示原始破坏区域，由人工分割获取；P为本研究方法识别结果。Dice取值范围为0～1。当Dice值越接近1时，表明识别结果越好；当Dice值越接近0时，识别结果越差。

由式（6）计算的图6中A组、B组数据Dice值表明，A组数据Dice得分为91.24%，B组数据Dice分值为93.53%，两者平均Dice为92.39%。计算结果充分证明本研究方法能够精确提取边坡形变区域。

通过上述方法对每组数据进行处理，获取每个时刻标记的边坡形变区域，然后利用识别的边坡形变区域的坐标与初始时刻坐标进行作差后积分即可获得随时间变化的累计边坡体积。如图7所示，图中两条曲线为本研究方法求取的两组试验中边坡累计体积随时间的变化曲线。在相同的坡脚情况下，降雨大小为65 mm/h的A组试验边坡比降雨大小为30 mm/h的B组试验边坡更早产生形变。

4 结 语

针对边坡连续形变信息难以精确提取的问题，提出了一种基于深度信息的边坡形变区域识别方法。该方法基于识别的形变区域前后时刻深度信息的变化能够准确估算边坡形变的体积。由于该方法仍以试验数据为主，因此采用真实数据对该方法进行进一步优化，是下一步的工作方向。