周媛媛，羅 斌，周 輝，趙 青，王 偉

(中國(guó)科學(xué)院武漢巖土力學(xué)研究所 測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430080)

我國(guó)煤礦開采深度以每年8～12 m的速度增加，預(yù)計(jì)在未來20年我國(guó)很多煤礦將進(jìn)入 1 000～1 500 m的深度[1].超深的礦道勢(shì)必有潛在的施工隱患，例如巷道變形劇烈、采場(chǎng)礦壓劇烈、采場(chǎng)失穩(wěn)、巖爆與沖擊地壓等由于巖體急劇不穩(wěn)定而造成的安全事故[2].因此需要對(duì)鉆孔的安全監(jiān)測(cè)評(píng)估來確保相關(guān)施工人員的安全.

對(duì)于鉆孔的安全評(píng)估可以視為一次管道檢查過程.管道的檢查方法根據(jù)是否發(fā)生直接接觸，分為接觸式檢測(cè)和非接觸式檢測(cè).非接觸式的方法規(guī)避了人工檢視的低效率高風(fēng)險(xiǎn)，具體分為CCD攝像法[3]、激光三角法[4]和激光投影法[5].CCD攝像法又稱為CCTV(closed circuit television)法[6]，利用管道機(jī)器人配載的CCD攝像頭對(duì)內(nèi)壁進(jìn)行拍攝，根據(jù)反饋的圖像視頻數(shù)據(jù)進(jìn)行相應(yīng)的處理分析.

機(jī)器人進(jìn)行管道檢測(cè)的研究如下，2009年 Bodenmann等[7]提出了基于稀疏特征對(duì)管道內(nèi)壁墻進(jìn)行可視化映射的管道視覺成像系統(tǒng)；2013年Kawasue等[8]提出了一種配備2臺(tái)環(huán)形激光器和1臺(tái)CCD相機(jī)進(jìn)行管道檢測(cè)的移動(dòng)機(jī)器人.國(guó)內(nèi)對(duì)于圖像的拼接的算法也有許多研究進(jìn)展，例如文獻(xiàn)[9]提出基于8參數(shù)透視映射的圖像拼接，優(yōu)點(diǎn)是匹配精準(zhǔn)，缺點(diǎn)是匹配效率低耗時(shí)長(zhǎng)；文獻(xiàn)[10]提出利用標(biāo)定參數(shù)將全景圖像重建成無畸變的內(nèi)壁展開圖像，該方法的弊端在于對(duì)需要已知直徑的管道進(jìn)行相應(yīng)的標(biāo)定測(cè)試.本文利用圖像處理的方式，對(duì)于視頻序列進(jìn)行拼接，三維化內(nèi)壁展開圖，相較于基于激光掃描和基于環(huán)形光與傳感器的重建方法相比，數(shù)據(jù)量小處理且效率高，可以滿足管道內(nèi)壁可視化檢視的需求.

1 視覺系統(tǒng)處理流程

基于鉆孔內(nèi)壁可視化的需求，整個(gè)視頻序列拼接處理方法分為前期機(jī)器人對(duì)鉆孔的視頻錄制，中期對(duì)錄制的視頻處理，后期輸出孔徑的三維模型建立，視頻序列拼接具體處理流程如圖1所示.

本地讀取錄制視頻，按照一定的時(shí)間間隔進(jìn)行采樣，得到關(guān)于鉆孔內(nèi)壁的序列圖像.對(duì)于單幀圖像，完成目標(biāo)環(huán)形區(qū)間圖像的展開，生成矩形條帶；對(duì)于序列圖像，完成幀間的特征匹配以及無縫拼接，得到鉆孔內(nèi)壁展開圖.利用得到的內(nèi)壁展開圖，本文提出基于openGL建模的方法，生成三維內(nèi)壁展開圖圖像模型，用于鉆孔安全狀況評(píng)估.

2 圖像分割及圓心提取

2.1 自適應(yīng)閾值圖像分割

機(jī)器人前端負(fù)責(zé)照明的燈帶滿足了拍攝區(qū)域附近的亮度需求，而鉆孔深部由于光線微弱而表現(xiàn)為黑色區(qū)域(如圖2中原始圖像所示)，利用圖像分割可以提取出鉆孔中心區(qū)域.由于依靠設(shè)定單一閾值分割無法滿足不同條件下拍攝的內(nèi)部序列圖像，因此本文采用基于最大類間誤差的OTSU[11-13]法對(duì)圖像進(jìn)行分割.

圖像可以視為含有目標(biāo)的前景區(qū)間和無用的背景區(qū)間兩個(gè)部分，OSTU通過遍歷迭代計(jì)算出分割閾值T，對(duì)圖像進(jìn)行分割后得到前景區(qū)間以及背景區(qū)間的灰度直方圖，T值滿足兩者的灰度直方圖的方差最大.

以圖像I(x,y)為例，當(dāng)前分割閾值為T.總像素MXN、灰度屬于[0,T]的歸為背景區(qū)間，統(tǒng)計(jì)背景像素總數(shù)為N1，平均灰度為μ1；灰度屬于(T,255]的歸為前景區(qū)間，統(tǒng)計(jì)前景像素總數(shù)為N0，平均灰度為M0.前景區(qū)間的像素比例為ω0，背景區(qū)間像素比例為ω1.圖像的平均灰度為μ，當(dāng)前的類間方差g的計(jì)算如下：

(1)

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

(2)

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

(3)

采用遍歷的方法最后確定分割閾值T，使類間方差g值達(dá)到最大.根據(jù)該閾值，將圖像分割為圓心區(qū)域Ac和非圓心區(qū)域A0(如圖2實(shí)驗(yàn)結(jié)果所示).

從實(shí)驗(yàn)結(jié)果中可以清楚看到采用隨機(jī)閾值進(jìn)行分割的時(shí)候無法徹底的將目標(biāo)區(qū)域與背景剝離，需要一定反復(fù)的試錯(cuò)才能得到理想的閾值；OSTU對(duì)圖像進(jìn)行分割后基本實(shí)現(xiàn)了提取孔徑中心區(qū)間的需求.初步閾值分割后，可以基本分離出圓心中心區(qū).如圖5所示，圓心區(qū)域 周邊存在細(xì)小的干擾區(qū)域，對(duì)于圓心的計(jì)算產(chǎn)生了一定的干擾，本文采用基于八鄰域跟蹤[14-15]的方法，根據(jù)對(duì)閉合邊緣的長(zhǎng)度判斷，對(duì)干擾區(qū)域進(jìn)行剔除.

為了對(duì)目標(biāo)區(qū)域邊界進(jìn)行編碼，首先設(shè)定目標(biāo)區(qū)域?yàn)槎祱D像，設(shè)前景為1，背景為0. 按照從左到右、從下到上的掃描方式對(duì)目標(biāo)區(qū)域邊緣像素點(diǎn)進(jìn)行搜索檢查，將最先檢查出來的白色像素點(diǎn)作為輪廓跟蹤的第一個(gè)輪廓點(diǎn)，若在圖像中搜索不到黑色像素點(diǎn)，則搜索結(jié)束[16]，實(shí)驗(yàn)結(jié)果如圖3所示.

跟蹤的邊緣集合記作{boundary|bi∈boundary,i∈(0,N]}，遍歷所有邊緣，邊緣長(zhǎng)度小于設(shè)置閾值的邊緣均標(biāo)記為干擾區(qū)域進(jìn)行剔除.最后可得不含干擾區(qū)域的圓心區(qū)域邊緣Bc,圓心區(qū)域Ac內(nèi)部點(diǎn)為Pc,{Ac|Pc(xi,yi),Pc∈Ac}，將Pc(xi,yi)像素值賦值為1，對(duì)輪廓內(nèi)部區(qū)域進(jìn)行填充.

2.2 圓心提取

得到圓心區(qū)域Ac后，對(duì)于圓心提取，本文進(jìn)行了2種方法的對(duì)比嘗試，第1種方法是基于最小距離場(chǎng)法[17-18]的內(nèi)切圓圓心提取，第2種方法是基于Ac區(qū)域內(nèi)像素坐標(biāo)的平均求取圓心.

內(nèi)切圓的圓心可以視作對(duì)原始圖像的一次中軸線的骨架提取的過程，骨架對(duì)噪聲、位置和旋轉(zhuǎn)等情況具有魯棒性，可以重建原物體，對(duì)于復(fù)原物體具有代表性.一個(gè)點(diǎn)的距離變換(DT, distance transform),定義為該點(diǎn)到邊界點(diǎn)的最短距離.距離變換的大小可以衡量該點(diǎn)與中軸的靠近程度，距離變換值的越大的點(diǎn)越靠近中軸[19].本文將歐拉距離作為距離的指標(biāo),對(duì)于2點(diǎn)P1(x1,y1)和P2(x2,y2)的歐拉距離d為：

(4)

本文對(duì)進(jìn)行干擾區(qū)域剔除后的孔徑圖像的輪廓跟蹤圖像作距離變換，尋找在圓心區(qū)域內(nèi)部具有最大DT值的點(diǎn)，將其作為內(nèi)切的圓心點(diǎn)(如圖5所示).

第2種方法將屬于圓心區(qū)域內(nèi)部點(diǎn)Pc,{Ac|Pc(xi,yi),Pc∈Ac}，橫縱坐標(biāo)分別求和后求得均值，視作圓心坐標(biāo):

(5)

兩次求得圓心坐標(biāo)結(jié)果對(duì)比如下(藍(lán)色加號(hào)為內(nèi)切圓法的圓心，紅色加號(hào)為坐標(biāo)均值法的圓心，右圖為放大示意圖):

可以發(fā)現(xiàn)內(nèi)切圓圓心Pc(xc,yc)與坐標(biāo)均值圓心Pa(xa,ya)有明顯偏移(見表1):

表1 2種方法的圓心坐標(biāo)

由于坐標(biāo)均值法容易受到不規(guī)則的形狀干擾，偏向突起端，因此本文采用內(nèi)切圓法確定.

3 環(huán)形圖像的展開及校正

對(duì)機(jī)器人返回的原始拍攝視頻以Δt時(shí)間間隔進(jìn)行采樣，得到關(guān)于鉆孔內(nèi)壁的圖像序列.為直觀描述和評(píng)價(jià)鉆孔內(nèi)的巖體結(jié)構(gòu)，需要對(duì)采樣所得的目標(biāo)幀需要進(jìn)行環(huán)形圖像的展開，具體步驟如下：

步驟1 確定孔徑的中心O(x,y)以及對(duì)應(yīng)的內(nèi)徑r1和外徑r2，r1和r2決定了展開矩形條帶的寬度，確定的O(x,y)的位置如果發(fā)生偏移，展開圖像會(huì)存在畸變；

步驟2 建立環(huán)形圖像到矩形圖像的投影關(guān)系，進(jìn)行圖像展開.

展開后的矩形圖像上一點(diǎn)P(x,y),對(duì)應(yīng)的原始環(huán)形圖像上坐標(biāo)為P′(x′,y′),在環(huán)形圖像中用極坐標(biāo)表示點(diǎn)P′,圓心坐標(biāo)為(x0,y0),P′的極半徑為ρ,方位角為θ,P′坐標(biāo)可以表示為:

(6)

以O(shè)點(diǎn)為圓心進(jìn)行展開，展開圖像高度為(r2-r1),點(diǎn)P(x,y)對(duì)應(yīng)的極半徑和方位角可以表示為:

(7)

將式(7)代入式(6)，建立環(huán)形圖像與矩形圖像之間的坐標(biāo)轉(zhuǎn)換關(guān)系:

(8)

根據(jù)式(8)可以對(duì)環(huán)形圖像進(jìn)行展開，得到矩形圖像.圖像展開后不難發(fā)現(xiàn)，圖像在徑向方向上有明顯的畸變(如圖7所示)，靠近圓心端被壓縮，遠(yuǎn)離圓心端被拉伸，因此我們引入棋盤格標(biāo)定的方法，加入畸變校正參數(shù)，對(duì)展開的圖像進(jìn)行畸變校正.

原始展開的矩形圖像上的點(diǎn)坐標(biāo)記作Pi(xi,yi)，校正后矩形圖像上點(diǎn)的坐標(biāo)記作Pi′(xi′,yi′)，本文利用多項(xiàng)式擬合的方式在徑向方向上進(jìn)行畸變參數(shù)的求解:

(9)

其中a,b,c3個(gè)參數(shù)為畸變參數(shù)，本文實(shí)驗(yàn)中采樣點(diǎn)為24個(gè)，利用采樣點(diǎn)進(jìn)行二項(xiàng)式擬合后求解畸變參數(shù)，對(duì)展開矩形進(jìn)行徑向畸變校正，9徑向方向上棋盤格寬度均勻，無拉伸或壓縮現(xiàn)象，改善了徑向的畸變問題(如圖8徑向校正圖像所示)，最后獲取橫縱系數(shù)比，得到水平校正圖像,校正為標(biāo)準(zhǔn)的棋盤格圖像.

計(jì)算所得畸變校正參數(shù):[a,b,c]=[-0.007，1.582，-3.464]

通過對(duì)畸變校正后的圖像提取角點(diǎn)，提取16×4個(gè)方格對(duì)其橫縱間隔進(jìn)行誤差檢驗(yàn),結(jié)果見圖9.

方格的水平方向邊長(zhǎng)記作橫向間隔，方格的縱向邊長(zhǎng)記作縱向間隔.統(tǒng)計(jì)64個(gè)方格的橫向、縱向間隔，求其均值、標(biāo)準(zhǔn)差，計(jì)算相對(duì)于標(biāo)準(zhǔn)邊長(zhǎng)的相對(duì)誤差(見表2).結(jié)果顯示相對(duì)誤差在3%以內(nèi)，具備良好的校正結(jié)果.

表2 35pixel×35pixel棋盤格圖像檢測(cè)結(jié)果

4 序列圖像的拼接

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

圖像拼接就是把針對(duì)同一場(chǎng)景的相互有部分重疊的一系列圖片合成一張大的寬視角的圖像.拼接后的圖像要求最大程度地與原始圖像接近,失真盡可能小,沒有明顯的縫合線[20].圖像拼接最根本的問題在于如何解決圖像重疊區(qū)域的配準(zhǔn).本文采取的匹配策略是基于模板匹配.模板匹配根據(jù)配準(zhǔn)方法的倚重不同，主要分為兩類:第1大類是基于特征的模板配準(zhǔn)，第2大類是基于灰度值的模板配準(zhǔn)[21].

基于灰度相關(guān)的目標(biāo)匹配算法通過計(jì)算模板圖像與目標(biāo)圖像的灰度相關(guān)系數(shù)(grayscale correlation coefficient)來判斷目標(biāo)圖像中是否存在與目標(biāo)相同或相近似的圖像.

基于幾何特征的模板匹配算法對(duì)特征進(jìn)行訓(xùn)練，通過訓(xùn)練后得到的特征對(duì)目標(biāo)圖像進(jìn)行特定搜索匹配，對(duì)于存在變形、旋轉(zhuǎn)以及光照條件變化下的圖像匹配具備良好的效果.

基于灰度相關(guān)的目標(biāo)匹配算法需要遍歷圖像，實(shí)際計(jì)算量偏大.在實(shí)際應(yīng)用過程中滿足效率的需求，本文提出利用harris算法[22]提取特征點(diǎn)，建立圍繞特征點(diǎn)構(gòu)成的搜索模板，基于灰度進(jìn)行匹配，提高單純基于灰度模板匹配的計(jì)算效率，滿足實(shí)際生產(chǎn)要求.

(10)

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

式(10)中G(s)為高斯模板，gx,gy為在x, y方向上的灰度梯度，M為求得的灰度自相關(guān)函數(shù)，利用對(duì)自相關(guān)函數(shù)M特征值的計(jì)算求解點(diǎn)的特征值.利用非極大值抑制算法(non-maximum suppression, NMS)搜索局部極大值，在鄰域中選擇最優(yōu)點(diǎn).最后根據(jù)距離限制，對(duì)所有的局部極值點(diǎn)進(jìn)行排序,根據(jù)需求數(shù)目確定最終的harris角點(diǎn).提取特征點(diǎn)后，以特征點(diǎn)為中心，建立大小為21×21像素的目標(biāo)模板，用于搜索.

在模板匹配的過程中，相似度的計(jì)算分為以基于差值平方和(sum of squared differences,簡(jiǎn)稱SSD)[23]的匹配和基于相關(guān)系數(shù)(normalized cross-correlation,簡(jiǎn)稱NCC)[24]的匹配.

SSD是最基本的一種評(píng)價(jià)機(jī)制，假設(shè)原模板大小為M×N，目標(biāo)圖像為T，搜索圖像為S:

(11)

T(m,n)為目標(biāo)圖像在(m,n)處的灰度值,通過計(jì)算與搜索圖像每一點(diǎn)S(i,j)處的構(gòu)成的大小為M×N模板灰度值差的平方來評(píng)價(jià)相似度，SSD值越小，相似度越高.

NCC是通過對(duì)SSD進(jìn)行歸一化處理，得到的歸一化系數(shù):

(12)

NCC的值域?yàn)閇0,1],NCC的值越接近于1，相似度越高，NCC=1證明與原圖像一致.

4.2 誤匹配剔除

匹配的特征點(diǎn)分為目標(biāo)點(diǎn)點(diǎn)集和匹配點(diǎn)點(diǎn)集，中間存在一定的低精度的匹配點(diǎn)，因?yàn)楸疚睦肦ANSAC[25-26]算法對(duì)誤匹配點(diǎn)進(jìn)行剔除，提高序列圖像拼接精度.

RANSAC算法通過尋找一個(gè)最佳單應(yīng)性矩陣H，使得滿足該矩陣的數(shù)據(jù)點(diǎn)個(gè)數(shù)最多.H矩陣的大小為3×3，滿足歸一化要求使h33=1，H矩陣內(nèi)含8個(gè)待解參數(shù)，因此至少選取4對(duì)不共線的初始匹配對(duì)進(jìn)行參數(shù)求解：

(13)

其中，s代表尺度參數(shù)，(x,y)代表目標(biāo)點(diǎn)坐標(biāo)，(x′,y′)代表匹配點(diǎn)坐標(biāo).計(jì)算出單應(yīng)性矩陣后利用這個(gè)模型測(cè)試所有數(shù)據(jù)，計(jì)算滿足這個(gè)模型數(shù)據(jù)點(diǎn)的個(gè)數(shù)與投影誤差(即代價(jià)cost函數(shù))，若此模型為最優(yōu)模型，則對(duì)應(yīng)的代價(jià)函數(shù)最小.

(14)

5 開發(fā)測(cè)試及結(jié)果分析

機(jī)器人采用USBFHD01M攝像頭進(jìn)行拍攝，視頻處理基于VS2015開發(fā)環(huán)境進(jìn)行了方法的實(shí)現(xiàn)，利用QT實(shí)現(xiàn)了界面的編寫，界面如圖11所示,左側(cè)顯示當(dāng)前讀取的視頻，藍(lán)色圓環(huán)為當(dāng)前選取的目標(biāo)圓環(huán)，用于隨后的展開以及圖像匹配，右側(cè)顯示拼接的結(jié)果.下端顯示處理的幀數(shù)以及對(duì)應(yīng)的匹配點(diǎn)數(shù)，計(jì)算幀間的位移誤差和圓心的偏移，進(jìn)行序列圖像拼接.

表3 測(cè)試視頻相關(guān)信息

5.1 拼接測(cè)試結(jié)果分析

此次測(cè)試將打印的紋理圖案貼于PC管內(nèi)壁，拍攝了一段含有紋理信息的鉆孔內(nèi)壁視頻. 從內(nèi)壁展開圖中截取紋理部分拼接圖片，與用手機(jī)拍攝的實(shí)際紋理圖樣進(jìn)行特征匹配以及灰度直方圖匹配 (如圖12、13所示)

表4 拼接圖片精度評(píng)價(jià)

通過對(duì)原始拼接圖像進(jìn)行特征匹配，得到38對(duì)有效匹配對(duì)，利用匹配對(duì)生成仿射變換模型，并利用灰度直方圖進(jìn)行灰度配準(zhǔn)，計(jì)算2幅圖像的MSE以及PSNR值作為精度評(píng)價(jià)指標(biāo)從變換矩陣可以看出拼接圖片與原始圖片之間只存在小尺度的變化和旋轉(zhuǎn)，目視結(jié)果和PSNR結(jié)果可以說明拼接圖像基本還原了原始紋理分布，不足之處在于測(cè)試所用攝像頭分辨率偏低，對(duì)于精度造成了一定的影響.

5.2 3維內(nèi)壁展開圖圖像

本文采用基于openGL立體建模貼圖的方法，將得到的鉆孔內(nèi)壁展開圖像還原為三維圓柱體，得到三維內(nèi)壁展開圖像，利用“↑”、“↓”實(shí)現(xiàn)觀察視角的切換，利用“←”、“→”控制當(dāng)前模型的自動(dòng)旋轉(zhuǎn)速度，利用平移面板的方向鍵控制模型所處的平面位置，實(shí)現(xiàn)全面的360度視角對(duì)內(nèi)壁展開圖進(jìn)行觀察.

6 結(jié)語

鉆孔爬行機(jī)器人視覺系統(tǒng)完美地規(guī)避了人工進(jìn)行施工安全檢查的局限性和危險(xiǎn)性，依靠機(jī)器人進(jìn)入狹長(zhǎng)細(xì)小的鉆孔進(jìn)行巖體安全的檢查和評(píng)價(jià)，保障施工安全，對(duì)于實(shí)際生產(chǎn)有著重要意義.本文提出的基于圖像特征的序列圖像拼接具備快速精確的特征，利用機(jī)器人返回的拍攝視頻即可實(shí)現(xiàn)對(duì)鉆孔內(nèi)壁的還原，針對(duì)現(xiàn)有可視化的需求，提供了三維化的內(nèi)壁圖像處理，可以用于多視角觀察鉆孔內(nèi)壁狀況，比二維平面展開圖像更為直觀精準(zhǔn).

機(jī)器人現(xiàn)階段針對(duì)的孔徑為平滑狹長(zhǎng)的孔徑，希望在未來的研究階段，能進(jìn)一步針對(duì)復(fù)雜孔徑進(jìn)行更深入的研究，滿足多樣化的施工需求.