周媛媛，罗 斌，周 辉，赵 青，王 伟

(中国科学院武汉岩土力学研究所 测绘遥感信息工程国家重点实验室，湖北 武汉 430080)

我国煤矿开采深度以每年8～12 m的速度增加，预计在未来20年我国很多煤矿将进入 1 000～1 500 m的深度[1].超深的矿道势必有潜在的施工隐患，例如巷道变形剧烈、采场矿压剧烈、采场失稳、岩爆与冲击地压等由于岩体急剧不稳定而造成的安全事故[2].因此需要对钻孔的安全监测评估来确保相关施工人员的安全.

对于钻孔的安全评估可以视为一次管道检查过程.管道的检查方法根据是否发生直接接触，分为接触式检测和非接触式检测.非接触式的方法规避了人工检视的低效率高风险，具体分为CCD摄像法[3]、激光三角法[4]和激光投影法[5].CCD摄像法又称为CCTV(closed circuit television)法[6]，利用管道机器人配载的CCD摄像头对内壁进行拍摄，根据反馈的图像视频数据进行相应的处理分析.

机器人进行管道检测的研究如下，2009年 Bodenmann等[7]提出了基于稀疏特征对管道内壁墙进行可视化映射的管道视觉成像系统；2013年Kawasue等[8]提出了一种配备2台环形激光器和1台CCD相机进行管道检测的移动机器人.国内对于图像的拼接的算法也有许多研究进展，例如文献[9]提出基于8参数透视映射的图像拼接，优点是匹配精准，缺点是匹配效率低耗时长；文献[10]提出利用标定参数将全景图像重建成无畸变的内壁展开图像，该方法的弊端在于对需要已知直径的管道进行相应的标定测试.本文利用图像处理的方式，对于视频序列进行拼接，三维化内壁展开图，相较于基于激光扫描和基于环形光与传感器的重建方法相比，数据量小处理且效率高，可以满足管道内壁可视化检视的需求.

1 视觉系统处理流程

基于钻孔内壁可视化的需求，整个视频序列拼接处理方法分为前期机器人对钻孔的视频录制，中期对录制的视频处理，后期输出孔径的三维模型建立，视频序列拼接具体处理流程如图1所示.

本地读取录制视频，按照一定的时间间隔进行采样，得到关于钻孔内壁的序列图像.对于单帧图像，完成目标环形区间图像的展开，生成矩形条带；对于序列图像，完成帧间的特征匹配以及无缝拼接，得到钻孔内壁展开图.利用得到的内壁展开图，本文提出基于openGL建模的方法，生成三维内壁展开图图像模型，用于钻孔安全状况评估.

2 图像分割及圆心提取

2.1 自适应阈值图像分割

机器人前端负责照明的灯带满足了拍摄区域附近的亮度需求，而钻孔深部由于光线微弱而表现为黑色区域(如图2中原始图像所示)，利用图像分割可以提取出钻孔中心区域.由于依靠设定单一阈值分割无法满足不同条件下拍摄的内部序列图像，因此本文采用基于最大类间误差的OTSU[11-13]法对图像进行分割.

图像可以视为含有目标的前景区间和无用的背景区间两个部分，OSTU通过遍历迭代计算出分割阈值T，对图像进行分割后得到前景区间以及背景区间的灰度直方图，T值满足两者的灰度直方图的方差最大.

以图像I(x,y)为例，当前分割阈值为T.总像素MXN、灰度属于[0,T]的归为背景区间，统计背景像素总数为N1，平均灰度为μ1；灰度属于(T,255]的归为前景区间，统计前景像素总数为N0，平均灰度为M0.前景区间的像素比例为ω0，背景区间像素比例为ω1.图像的平均灰度为μ，当前的类间方差g的计算如下：

(1)

μ=ω0×μ0+ω1×μ1；

(2)

g=ω0(μ0-μ)2+ω1(μ1-μ)2.

(3)

采用遍历的方法最后确定分割阈值T，使类间方差g值达到最大.根据该阈值，将图像分割为圆心区域Ac和非圆心区域A0(如图2实验结果所示).

从实验结果中可以清楚看到采用随机阈值进行分割的时候无法彻底的将目标区域与背景剥离，需要一定反复的试错才能得到理想的阈值；OSTU对图像进行分割后基本实现了提取孔径中心区间的需求.初步阈值分割后，可以基本分离出圆心中心区.如图5所示，圆心区域 周边存在细小的干扰区域，对于圆心的计算产生了一定的干扰，本文采用基于八邻域跟踪[14-15]的方法，根据对闭合边缘的长度判断，对干扰区域进行剔除.

为了对目标区域边界进行编码，首先设定目标区域为二值图像，设前景为1，背景为0. 按照从左到右、从下到上的扫描方式对目标区域边缘像素点进行搜索检查，将最先检查出来的白色像素点作为轮廓跟踪的第一个轮廓点，若在图像中搜索不到黑色像素点，则搜索结束[16]，实验结果如图3所示.

跟踪的边缘集合记作{boundary|bi∈boundary,i∈(0,N]}，遍历所有边缘，边缘长度小于设置阈值的边缘均标记为干扰区域进行剔除.最后可得不含干扰区域的圆心区域边缘Bc,圆心区域Ac内部点为Pc,{Ac|Pc(xi,yi),Pc∈Ac}，将Pc(xi,yi)像素值赋值为1，对轮廓内部区域进行填充.

2.2 圆心提取

得到圆心区域Ac后，对于圆心提取，本文进行了2种方法的对比尝试，第1种方法是基于最小距离场法[17-18]的内切圆圆心提取，第2种方法是基于Ac区域内像素坐标的平均求取圆心.

内切圆的圆心可以视作对原始图像的一次中轴线的骨架提取的过程，骨架对噪声、位置和旋转等情况具有鲁棒性，可以重建原物体，对于复原物体具有代表性.一个点的距离变换(DT, distance transform),定义为该点到边界点的最短距离.距离变换的大小可以衡量该点与中轴的靠近程度，距离变换值的越大的点越靠近中轴[19].本文将欧拉距离作为距离的指标,对于2点P1(x1,y1)和P2(x2,y2)的欧拉距离d为：

(4)

本文对进行干扰区域剔除后的孔径图像的轮廓跟踪图像作距离变换，寻找在圆心区域内部具有最大DT值的点，将其作为内切的圆心点(如图5所示).

第2种方法将属于圆心区域内部点Pc,{Ac|Pc(xi,yi),Pc∈Ac}，横纵坐标分别求和后求得均值，视作圆心坐标:

(5)

两次求得圆心坐标结果对比如下(蓝色加号为内切圆法的圆心，红色加号为坐标均值法的圆心，右图为放大示意图):

可以发现内切圆圆心Pc(xc,yc)与坐标均值圆心Pa(xa,ya)有明显偏移(见表1):

表1 2种方法的圆心坐标

由于坐标均值法容易受到不规则的形状干扰，偏向突起端，因此本文采用内切圆法确定.

3 环形图像的展开及校正

对机器人返回的原始拍摄视频以Δt时间间隔进行采样，得到关于钻孔内壁的图像序列.为直观描述和评价钻孔内的岩体结构，需要对采样所得的目标帧需要进行环形图像的展开，具体步骤如下：

步骤1 确定孔径的中心O(x,y)以及对应的内径r1和外径r2，r1和r2决定了展开矩形条带的宽度，确定的O(x,y)的位置如果发生偏移，展开图像会存在畸变；

步骤2 建立环形图像到矩形图像的投影关系，进行图像展开.

展开后的矩形图像上一点P(x,y),对应的原始环形图像上坐标为P′(x′,y′),在环形图像中用极坐标表示点P′,圆心坐标为(x0,y0),P′的极半径为ρ,方位角为θ,P′坐标可以表示为:

(6)

以O点为圆心进行展开，展开图像高度为(r2-r1),点P(x,y)对应的极半径和方位角可以表示为:

(7)

将式(7)代入式(6)，建立环形图像与矩形图像之间的坐标转换关系:

(8)

根据式(8)可以对环形图像进行展开，得到矩形图像.图像展开后不难发现，图像在径向方向上有明显的畸变(如图7所示)，靠近圆心端被压缩，远离圆心端被拉伸，因此我们引入棋盘格标定的方法，加入畸变校正参数，对展开的图像进行畸变校正.

原始展开的矩形图像上的点坐标记作Pi(xi,yi)，校正后矩形图像上点的坐标记作Pi′(xi′,yi′)，本文利用多项式拟合的方式在径向方向上进行畸变参数的求解:

(9)

其中a,b,c3个参数为畸变参数，本文实验中采样点为24个，利用采样点进行二项式拟合后求解畸变参数，对展开矩形进行径向畸变校正，9径向方向上棋盘格宽度均匀，无拉伸或压缩现象，改善了径向的畸变问题(如图8径向校正图像所示)，最后获取横纵系数比，得到水平校正图像,校正为标准的棋盘格图像.

计算所得畸变校正参数:[a,b,c]=[-0.007，1.582，-3.464]

通过对畸变校正后的图像提取角点，提取16×4个方格对其横纵间隔进行误差检验,结果见图9.

方格的水平方向边长记作横向间隔，方格的纵向边长记作纵向间隔.统计64个方格的横向、纵向间隔，求其均值、标准差，计算相对于标准边长的相对误差(见表2).结果显示相对误差在3%以内，具备良好的校正结果.

表2 35pixel×35pixel棋盘格图像检测结果

4 序列图像的拼接

4.1基于模板匹配的特征提取

图像拼接就是把针对同一场景的相互有部分重叠的一系列图片合成一张大的宽视角的图像.拼接后的图像要求最大程度地与原始图像接近,失真尽可能小,没有明显的缝合线[20].图像拼接最根本的问题在于如何解决图像重叠区域的配准.本文采取的匹配策略是基于模板匹配.模板匹配根据配准方法的倚重不同，主要分为两类:第1大类是基于特征的模板配准，第2大类是基于灰度值的模板配准[21].

基于灰度相关的目标匹配算法通过计算模板图像与目标图像的灰度相关系数(grayscale correlation coefficient)来判断目标图像中是否存在与目标相同或相近似的图像.

基于几何特征的模板匹配算法对特征进行训练，通过训练后得到的特征对目标图像进行特定搜索匹配，对于存在变形、旋转以及光照条件变化下的图像匹配具备良好的效果.

基于灰度相关的目标匹配算法需要遍历图像，实际计算量偏大.在实际应用过程中满足效率的需求，本文提出利用harris算法[22]提取特征点，建立围绕特征点构成的搜索模板，基于灰度进行匹配，提高单纯基于灰度模板匹配的计算效率，满足实际生产要求.

(10)

I=det(M)-ktr2(M),k=0.04.

式(10)中G(s)为高斯模板，gx,gy为在x, y方向上的灰度梯度，M为求得的灰度自相关函数，利用对自相关函数M特征值的计算求解点的特征值.利用非极大值抑制算法(non-maximum suppression, NMS)搜索局部极大值，在邻域中选择最优点.最后根据距离限制，对所有的局部极值点进行排序,根据需求数目确定最终的harris角点.提取特征点后，以特征点为中心，建立大小为21×21像素的目标模板，用于搜索.

在模板匹配的过程中，相似度的计算分为以基于差值平方和(sum of squared differences,简称SSD)[23]的匹配和基于相关系数(normalized cross-correlation,简称NCC)[24]的匹配.

SSD是最基本的一种评价机制，假设原模板大小为M×N，目标图像为T，搜索图像为S:

(11)

T(m,n)为目标图像在(m,n)处的灰度值,通过计算与搜索图像每一点S(i,j)处的构成的大小为M×N模板灰度值差的平方来评价相似度，SSD值越小，相似度越高.

NCC是通过对SSD进行归一化处理，得到的归一化系数:

(12)

NCC的值域为[0,1],NCC的值越接近于1，相似度越高，NCC=1证明与原图像一致.

4.2 误匹配剔除

匹配的特征点分为目标点点集和匹配点点集，中间存在一定的低精度的匹配点，因为本文利用RANSAC[25-26]算法对误匹配点进行剔除，提高序列图像拼接精度.

RANSAC算法通过寻找一个最佳单应性矩阵H，使得满足该矩阵的数据点个数最多.H矩阵的大小为3×3，满足归一化要求使h33=1，H矩阵内含8个待解参数，因此至少选取4对不共线的初始匹配对进行参数求解：

(13)

其中，s代表尺度参数，(x,y)代表目标点坐标，(x′,y′)代表匹配点坐标.计算出单应性矩阵后利用这个模型测试所有数据，计算满足这个模型数据点的个数与投影误差(即代价cost函数)，若此模型为最优模型，则对应的代价函数最小.

(14)

5 开发测试及结果分析

机器人采用USBFHD01M摄像头进行拍摄，视频处理基于VS2015开发环境进行了方法的实现，利用QT实现了界面的编写，界面如图11所示,左侧显示当前读取的视频，蓝色圆环为当前选取的目标圆环，用于随后的展开以及图像匹配，右侧显示拼接的结果.下端显示处理的帧数以及对应的匹配点数，计算帧间的位移误差和圆心的偏移，进行序列图像拼接.

表3 测试视频相关信息

5.1 拼接测试结果分析

此次测试将打印的纹理图案贴于PC管内壁，拍摄了一段含有纹理信息的钻孔内壁视频. 从内壁展开图中截取纹理部分拼接图片，与用手机拍摄的实际纹理图样进行特征匹配以及灰度直方图匹配 (如图12、13所示)

表4 拼接图片精度评价

通过对原始拼接图像进行特征匹配，得到38对有效匹配对，利用匹配对生成仿射变换模型，并利用灰度直方图进行灰度配准，计算2幅图像的MSE以及PSNR值作为精度评价指标从变换矩阵可以看出拼接图片与原始图片之间只存在小尺度的变化和旋转，目视结果和PSNR结果可以说明拼接图像基本还原了原始纹理分布，不足之处在于测试所用摄像头分辨率偏低，对于精度造成了一定的影响.

5.2 3维内壁展开图图像

本文采用基于openGL立体建模贴图的方法，将得到的钻孔内壁展开图像还原为三维圆柱体，得到三维内壁展开图像，利用“↑”、“↓”实现观察视角的切换，利用“←”、“→”控制当前模型的自动旋转速度，利用平移面板的方向键控制模型所处的平面位置，实现全面的360度视角对内壁展开图进行观察.

6 结语

钻孔爬行机器人视觉系统完美地规避了人工进行施工安全检查的局限性和危险性，依靠机器人进入狭长细小的钻孔进行岩体安全的检查和评价，保障施工安全，对于实际生产有着重要意义.本文提出的基于图像特征的序列图像拼接具备快速精确的特征，利用机器人返回的拍摄视频即可实现对钻孔内壁的还原，针对现有可视化的需求，提供了三维化的内壁图像处理，可以用于多视角观察钻孔内壁状况，比二维平面展开图像更为直观精准.

机器人现阶段针对的孔径为平滑狭长的孔径，希望在未来的研究阶段，能进一步针对复杂孔径进行更深入的研究，满足多样化的施工需求.