周媛媛,羅 斌,周 輝,趙 青,王 偉
(中國科學(xué)院武漢巖土力學(xué)研究所 測繪遙感信息工程國家重點實驗室,湖北 武漢 430080)
我國煤礦開采深度以每年8~12 m的速度增加,預(yù)計在未來20年我國很多煤礦將進入 1 000~1 500 m的深度[1].超深的礦道勢必有潛在的施工隱患,例如巷道變形劇烈、采場礦壓劇烈、采場失穩(wěn)、巖爆與沖擊地壓等由于巖體急劇不穩(wěn)定而造成的安全事故[2].因此需要對鉆孔的安全監(jiān)測評估來確保相關(guān)施工人員的安全.
對于鉆孔的安全評估可以視為一次管道檢查過程.管道的檢查方法根據(jù)是否發(fā)生直接接觸,分為接觸式檢測和非接觸式檢測.非接觸式的方法規(guī)避了人工檢視的低效率高風(fēng)險,具體分為CCD攝像法[3]、激光三角法[4]和激光投影法[5].CCD攝像法又稱為CCTV(closed circuit television)法[6],利用管道機器人配載的CCD攝像頭對內(nèi)壁進行拍攝,根據(jù)反饋的圖像視頻數(shù)據(jù)進行相應(yīng)的處理分析.
機器人進行管道檢測的研究如下,2009年 Bodenmann等[7]提出了基于稀疏特征對管道內(nèi)壁墻進行可視化映射的管道視覺成像系統(tǒng);2013年Kawasue等[8]提出了一種配備2臺環(huán)形激光器和1臺CCD相機進行管道檢測的移動機器人.國內(nèi)對于圖像的拼接的算法也有許多研究進展,例如文獻[9]提出基于8參數(shù)透視映射的圖像拼接,優(yōu)點是匹配精準(zhǔn),缺點是匹配效率低耗時長;文獻[10]提出利用標(biāo)定參數(shù)將全景圖像重建成無畸變的內(nèi)壁展開圖像,該方法的弊端在于對需要已知直徑的管道進行相應(yīng)的標(biāo)定測試.本文利用圖像處理的方式,對于視頻序列進行拼接,三維化內(nèi)壁展開圖,相較于基于激光掃描和基于環(huán)形光與傳感器的重建方法相比,數(shù)據(jù)量小處理且效率高,可以滿足管道內(nèi)壁可視化檢視的需求.
基于鉆孔內(nèi)壁可視化的需求,整個視頻序列拼接處理方法分為前期機器人對鉆孔的視頻錄制,中期對錄制的視頻處理,后期輸出孔徑的三維模型建立,視頻序列拼接具體處理流程如圖1所示.
本地讀取錄制視頻,按照一定的時間間隔進行采樣,得到關(guān)于鉆孔內(nèi)壁的序列圖像.對于單幀圖像,完成目標(biāo)環(huán)形區(qū)間圖像的展開,生成矩形條帶;對于序列圖像,完成幀間的特征匹配以及無縫拼接,得到鉆孔內(nèi)壁展開圖.利用得到的內(nèi)壁展開圖,本文提出基于openGL建模的方法,生成三維內(nèi)壁展開圖圖像模型,用于鉆孔安全狀況評估.
機器人前端負責(zé)照明的燈帶滿足了拍攝區(qū)域附近的亮度需求,而鉆孔深部由于光線微弱而表現(xiàn)為黑色區(qū)域(如圖2中原始圖像所示),利用圖像分割可以提取出鉆孔中心區(qū)域.由于依靠設(shè)定單一閾值分割無法滿足不同條件下拍攝的內(nèi)部序列圖像,因此本文采用基于最大類間誤差的OTSU[11-13]法對圖像進行分割.
圖像可以視為含有目標(biāo)的前景區(qū)間和無用的背景區(qū)間兩個部分,OSTU通過遍歷迭代計算出分割閾值T,對圖像進行分割后得到前景區(qū)間以及背景區(qū)間的灰度直方圖,T值滿足兩者的灰度直方圖的方差最大.
以圖像I(x,y)為例,當(dāng)前分割閾值為T.總像素MXN、灰度屬于[0,T]的歸為背景區(qū)間,統(tǒng)計背景像素總數(shù)為N1,平均灰度為μ1;灰度屬于(T,255]的歸為前景區(qū)間,統(tǒng)計前景像素總數(shù)為N0,平均灰度為M0.前景區(qū)間的像素比例為ω0,背景區(qū)間像素比例為ω1.圖像的平均灰度為μ,當(dāng)前的類間方差g的計算如下:
(1)
μ=ω0×μ0+ω1×μ1;
(2)
g=ω0(μ0-μ)2+ω1(μ1-μ)2.
(3)
采用遍歷的方法最后確定分割閾值T,使類間方差g值達到最大.根據(jù)該閾值,將圖像分割為圓心區(qū)域Ac和非圓心區(qū)域A0(如圖2實驗結(jié)果所示).
從實驗結(jié)果中可以清楚看到采用隨機閾值進行分割的時候無法徹底的將目標(biāo)區(qū)域與背景剝離,需要一定反復(fù)的試錯才能得到理想的閾值;OSTU對圖像進行分割后基本實現(xiàn)了提取孔徑中心區(qū)間的需求.初步閾值分割后,可以基本分離出圓心中心區(qū).如圖5所示,圓心區(qū)域 周邊存在細小的干擾區(qū)域,對于圓心的計算產(chǎn)生了一定的干擾,本文采用基于八鄰域跟蹤[14-15]的方法,根據(jù)對閉合邊緣的長度判斷,對干擾區(qū)域進行剔除.
為了對目標(biāo)區(qū)域邊界進行編碼,首先設(shè)定目標(biāo)區(qū)域為二值圖像,設(shè)前景為1,背景為0. 按照從左到右、從下到上的掃描方式對目標(biāo)區(qū)域邊緣像素點進行搜索檢查,將最先檢查出來的白色像素點作為輪廓跟蹤的第一個輪廓點,若在圖像中搜索不到黑色像素點,則搜索結(jié)束[16],實驗結(jié)果如圖3所示.
跟蹤的邊緣集合記作{boundary|bi∈boundary,i∈(0,N]},遍歷所有邊緣,邊緣長度小于設(shè)置閾值的邊緣均標(biāo)記為干擾區(qū)域進行剔除.最后可得不含干擾區(qū)域的圓心區(qū)域邊緣Bc,圓心區(qū)域Ac內(nèi)部點為Pc,{Ac|Pc(xi,yi),Pc∈Ac},將Pc(xi,yi)像素值賦值為1,對輪廓內(nèi)部區(qū)域進行填充.
得到圓心區(qū)域Ac后,對于圓心提取,本文進行了2種方法的對比嘗試,第1種方法是基于最小距離場法[17-18]的內(nèi)切圓圓心提取,第2種方法是基于Ac區(qū)域內(nèi)像素坐標(biāo)的平均求取圓心.
內(nèi)切圓的圓心可以視作對原始圖像的一次中軸線的骨架提取的過程,骨架對噪聲、位置和旋轉(zhuǎn)等情況具有魯棒性,可以重建原物體,對于復(fù)原物體具有代表性.一個點的距離變換(DT, distance transform),定義為該點到邊界點的最短距離.距離變換的大小可以衡量該點與中軸的靠近程度,距離變換值的越大的點越靠近中軸[19].本文將歐拉距離作為距離的指標(biāo),對于2點P1(x1,y1)和P2(x2,y2)的歐拉距離d為:
(4)
本文對進行干擾區(qū)域剔除后的孔徑圖像的輪廓跟蹤圖像作距離變換,尋找在圓心區(qū)域內(nèi)部具有最大DT值的點,將其作為內(nèi)切的圓心點(如圖5所示).
第2種方法將屬于圓心區(qū)域內(nèi)部點Pc,{Ac|Pc(xi,yi),Pc∈Ac},橫縱坐標(biāo)分別求和后求得均值,視作圓心坐標(biāo):
(5)
兩次求得圓心坐標(biāo)結(jié)果對比如下(藍色加號為內(nèi)切圓法的圓心,紅色加號為坐標(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)切圓法確定.
對機器人返回的原始拍攝視頻以Δt時間間隔進行采樣,得到關(guān)于鉆孔內(nèi)壁的圖像序列.為直觀描述和評價鉆孔內(nèi)的巖體結(jié)構(gòu),需要對采樣所得的目標(biāo)幀需要進行環(huán)形圖像的展開,具體步驟如下:
步驟1 確定孔徑的中心O(x,y)以及對應(yīng)的內(nèi)徑r1和外徑r2,r1和r2決定了展開矩形條帶的寬度,確定的O(x,y)的位置如果發(fā)生偏移,展開圖像會存在畸變;
步驟2 建立環(huán)形圖像到矩形圖像的投影關(guān)系,進行圖像展開.
展開后的矩形圖像上一點P(x,y),對應(yīng)的原始環(huán)形圖像上坐標(biāo)為P′(x′,y′),在環(huán)形圖像中用極坐標(biāo)表示點P′,圓心坐標(biāo)為(x0,y0),P′的極半徑為ρ,方位角為θ,P′坐標(biāo)可以表示為:
(6)
以O(shè)點為圓心進行展開,展開圖像高度為(r2-r1),點P(x,y)對應(yīng)的極半徑和方位角可以表示為:
(7)
將式(7)代入式(6),建立環(huán)形圖像與矩形圖像之間的坐標(biāo)轉(zhuǎn)換關(guān)系:
(8)
根據(jù)式(8)可以對環(huán)形圖像進行展開,得到矩形圖像.圖像展開后不難發(fā)現(xiàn),圖像在徑向方向上有明顯的畸變(如圖7所示),靠近圓心端被壓縮,遠離圓心端被拉伸,因此我們引入棋盤格標(biāo)定的方法,加入畸變校正參數(shù),對展開的圖像進行畸變校正.
原始展開的矩形圖像上的點坐標(biāo)記作Pi(xi,yi),校正后矩形圖像上點的坐標(biāo)記作Pi′(xi′,yi′),本文利用多項式擬合的方式在徑向方向上進行畸變參數(shù)的求解:
(9)
其中a,b,c3個參數(shù)為畸變參數(shù),本文實驗中采樣點為24個,利用采樣點進行二項式擬合后求解畸變參數(shù),對展開矩形進行徑向畸變校正,9徑向方向上棋盤格寬度均勻,無拉伸或壓縮現(xiàn)象,改善了徑向的畸變問題(如圖8徑向校正圖像所示),最后獲取橫縱系數(shù)比,得到水平校正圖像,校正為標(biāo)準(zhǔn)的棋盤格圖像.
計算所得畸變校正參數(shù):[a,b,c]=[-0.007,1.582,-3.464]
通過對畸變校正后的圖像提取角點,提取16×4個方格對其橫縱間隔進行誤差檢驗,結(jié)果見圖9.
方格的水平方向邊長記作橫向間隔,方格的縱向邊長記作縱向間隔.統(tǒng)計64個方格的橫向、縱向間隔,求其均值、標(biāo)準(zhǔn)差,計算相對于標(biāo)準(zhǔn)邊長的相對誤差(見表2).結(jié)果顯示相對誤差在3%以內(nèi),具備良好的校正結(jié)果.
表2 35pixel×35pixel棋盤格圖像檢測結(jié)果
4.1基于模板匹配的特征提取
圖像拼接就是把針對同一場景的相互有部分重疊的一系列圖片合成一張大的寬視角的圖像.拼接后的圖像要求最大程度地與原始圖像接近,失真盡可能小,沒有明顯的縫合線[20].圖像拼接最根本的問題在于如何解決圖像重疊區(qū)域的配準(zhǔn).本文采取的匹配策略是基于模板匹配.模板匹配根據(jù)配準(zhǔn)方法的倚重不同,主要分為兩類:第1大類是基于特征的模板配準(zhǔn),第2大類是基于灰度值的模板配準(zhǔn)[21].
基于灰度相關(guān)的目標(biāo)匹配算法通過計算模板圖像與目標(biāo)圖像的灰度相關(guān)系數(shù)(grayscale correlation coefficient)來判斷目標(biāo)圖像中是否存在與目標(biāo)相同或相近似的圖像.
基于幾何特征的模板匹配算法對特征進行訓(xùn)練,通過訓(xùn)練后得到的特征對目標(biāo)圖像進行特定搜索匹配,對于存在變形、旋轉(zhuǎn)以及光照條件變化下的圖像匹配具備良好的效果.
基于灰度相關(guān)的目標(biāo)匹配算法需要遍歷圖像,實際計算量偏大.在實際應(yīng)用過程中滿足效率的需求,本文提出利用harris算法[22]提取特征點,建立圍繞特征點構(gòu)成的搜索模板,基于灰度進行匹配,提高單純基于灰度模板匹配的計算效率,滿足實際生產(chǎn)要求.
(10)
I=det(M)-ktr2(M),k=0.04.
式(10)中G(s)為高斯模板,gx,gy為在x, y方向上的灰度梯度,M為求得的灰度自相關(guān)函數(shù),利用對自相關(guān)函數(shù)M特征值的計算求解點的特征值.利用非極大值抑制算法(non-maximum suppression, NMS)搜索局部極大值,在鄰域中選擇最優(yōu)點.最后根據(jù)距離限制,對所有的局部極值點進行排序,根據(jù)需求數(shù)目確定最終的harris角點.提取特征點后,以特征點為中心,建立大小為21×21像素的目標(biāo)模板,用于搜索.
在模板匹配的過程中,相似度的計算分為以基于差值平方和(sum of squared differences,簡稱SSD)[23]的匹配和基于相關(guān)系數(shù)(normalized cross-correlation,簡稱NCC)[24]的匹配.
SSD是最基本的一種評價機制,假設(shè)原模板大小為M×N,目標(biāo)圖像為T,搜索圖像為S:
(11)
T(m,n)為目標(biāo)圖像在(m,n)處的灰度值,通過計算與搜索圖像每一點S(i,j)處的構(gòu)成的大小為M×N模板灰度值差的平方來評價相似度,SSD值越小,相似度越高.
NCC是通過對SSD進行歸一化處理,得到的歸一化系數(shù):
(12)
NCC的值域為[0,1],NCC的值越接近于1,相似度越高,NCC=1證明與原圖像一致.
匹配的特征點分為目標(biāo)點點集和匹配點點集,中間存在一定的低精度的匹配點,因為本文利用RANSAC[25-26]算法對誤匹配點進行剔除,提高序列圖像拼接精度.
RANSAC算法通過尋找一個最佳單應(yīng)性矩陣H,使得滿足該矩陣的數(shù)據(jù)點個數(shù)最多.H矩陣的大小為3×3,滿足歸一化要求使h33=1,H矩陣內(nèi)含8個待解參數(shù),因此至少選取4對不共線的初始匹配對進行參數(shù)求解:
(13)
其中,s代表尺度參數(shù),(x,y)代表目標(biāo)點坐標(biāo),(x′,y′)代表匹配點坐標(biāo).計算出單應(yīng)性矩陣后利用這個模型測試所有數(shù)據(jù),計算滿足這個模型數(shù)據(jù)點的個數(shù)與投影誤差(即代價cost函數(shù)),若此模型為最優(yōu)模型,則對應(yīng)的代價函數(shù)最小.
(14)
機器人采用USBFHD01M攝像頭進行拍攝,視頻處理基于VS2015開發(fā)環(huán)境進行了方法的實現(xiàn),利用QT實現(xiàn)了界面的編寫,界面如圖11所示,左側(cè)顯示當(dāng)前讀取的視頻,藍色圓環(huán)為當(dāng)前選取的目標(biāo)圓環(huán),用于隨后的展開以及圖像匹配,右側(cè)顯示拼接的結(jié)果.下端顯示處理的幀數(shù)以及對應(yīng)的匹配點數(shù),計算幀間的位移誤差和圓心的偏移,進行序列圖像拼接.
表3 測試視頻相關(guān)信息
此次測試將打印的紋理圖案貼于PC管內(nèi)壁,拍攝了一段含有紋理信息的鉆孔內(nèi)壁視頻. 從內(nèi)壁展開圖中截取紋理部分拼接圖片,與用手機拍攝的實際紋理圖樣進行特征匹配以及灰度直方圖匹配 (如圖12、13所示)
表4 拼接圖片精度評價
通過對原始拼接圖像進行特征匹配,得到38對有效匹配對,利用匹配對生成仿射變換模型,并利用灰度直方圖進行灰度配準(zhǔn),計算2幅圖像的MSE以及PSNR值作為精度評價指標(biāo)從變換矩陣可以看出拼接圖片與原始圖片之間只存在小尺度的變化和旋轉(zhuǎn),目視結(jié)果和PSNR結(jié)果可以說明拼接圖像基本還原了原始紋理分布,不足之處在于測試所用攝像頭分辨率偏低,對于精度造成了一定的影響.
本文采用基于openGL立體建模貼圖的方法,將得到的鉆孔內(nèi)壁展開圖像還原為三維圓柱體,得到三維內(nèi)壁展開圖像,利用“↑”、“↓”實現(xiàn)觀察視角的切換,利用“←”、“→”控制當(dāng)前模型的自動旋轉(zhuǎn)速度,利用平移面板的方向鍵控制模型所處的平面位置,實現(xiàn)全面的360度視角對內(nèi)壁展開圖進行觀察.
鉆孔爬行機器人視覺系統(tǒng)完美地規(guī)避了人工進行施工安全檢查的局限性和危險性,依靠機器人進入狹長細小的鉆孔進行巖體安全的檢查和評價,保障施工安全,對于實際生產(chǎn)有著重要意義.本文提出的基于圖像特征的序列圖像拼接具備快速精確的特征,利用機器人返回的拍攝視頻即可實現(xiàn)對鉆孔內(nèi)壁的還原,針對現(xiàn)有可視化的需求,提供了三維化的內(nèi)壁圖像處理,可以用于多視角觀察鉆孔內(nèi)壁狀況,比二維平面展開圖像更為直觀精準(zhǔn).
機器人現(xiàn)階段針對的孔徑為平滑狹長的孔徑,希望在未來的研究階段,能進一步針對復(fù)雜孔徑進行更深入的研究,滿足多樣化的施工需求.