段志鑫，吳 侃，宋月婷，白志輝，李 亮，周大偉

(1.中國(guó)礦業(yè)大學(xué) 環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221116；2.武漢大學(xué) 遙感信息工程學(xué)院，湖北 武漢 430079； 3.冀中能源峰峰集團(tuán)有限公司，河北 邯鄲 056000)

相似材料模擬實(shí)驗(yàn)是研究巖體內(nèi)部壓力與巖層運(yùn)移規(guī)律的重要手段，由于其具有易于實(shí)現(xiàn)、觀測(cè)周期短、實(shí)驗(yàn)結(jié)果直觀等優(yōu)點(diǎn)，被廣泛地應(yīng)用于巖土工程和礦山開(kāi)采等領(lǐng)域[1]。在礦山開(kāi)采沉陷研究中，相似材料模擬實(shí)驗(yàn)的觀測(cè)內(nèi)容主要包括模擬巖層的位移和裂隙分布情況。目前，相似材料模型巖層位移監(jiān)測(cè)的方法比較成熟多樣，主要包括百分表法[2]、燈泡透鏡法[3]、全站儀法[4]、近景攝影測(cè)量法[5-6]、光纖光柵傳感器監(jiān)測(cè)法[7]等。但是對(duì)于相似材料模型在模擬推進(jìn)過(guò)程中裂隙發(fā)育和破裂形態(tài)信息，例如裂隙的發(fā)育高度、離層尺寸及位置等，由于觀測(cè)手段和數(shù)據(jù)處理方法的限制，研究較少，目前此類(lèi)信息的獲取方法仍以手工量測(cè)為主，不僅效率低，且精度較差。許家林等對(duì)基于圖像的相似材料模型采動(dòng)裂隙提取與量化方法進(jìn)行了大量研究，并開(kāi)發(fā)了對(duì)應(yīng)的模型裂隙處理程序FIMAGE，該方法可以較好地實(shí)現(xiàn)對(duì)背景單一的相似模型裂隙檢測(cè)，但是未考慮圖像畸變的影響[8-9]；張?jiān)迄i利用白光掃描儀獲取模型的點(diǎn)云并通過(guò)鄰域法提取模型破裂邊界，該方法數(shù)據(jù)處理效率較低，無(wú)法滿(mǎn)足實(shí)際應(yīng)用需求[10]；陳朋等將相似模型表面的三維點(diǎn)云數(shù)據(jù)轉(zhuǎn)換成深度圖像，并通過(guò)邊緣檢測(cè)算子提取模型破裂邊界，該方法觀測(cè)成本高，且對(duì)細(xì)小裂縫的提取狀況不佳[11]。

視覺(jué)測(cè)量技術(shù)是近些年工業(yè)檢測(cè)中最常用的手段之一，具有簡(jiǎn)單、高效、成本低廉等優(yōu)點(diǎn)[12]，在道路病害監(jiān)測(cè)[13]、工件缺損檢測(cè)[14]、結(jié)構(gòu)體測(cè)量[15]等領(lǐng)域得到了廣泛地應(yīng)用。針對(duì)目前相似材料模型裂隙檢測(cè)及定位方法自動(dòng)化程度低、精度較差等問(wèn)題，本文擬從二維數(shù)碼圖像出發(fā)，對(duì)基于視覺(jué)測(cè)量技術(shù)的相似材料模型裂縫提取和定位方法進(jìn)行研究，為研究相似模擬實(shí)驗(yàn)中模型覆巖層的破壞和“三帶”發(fā)育規(guī)律提供技術(shù)支持。

1 模型裂隙檢測(cè)

原始模型圖像仍然為三通道彩色圖像，直接基于該圖像進(jìn)行裂縫提取過(guò)程比較復(fù)雜，因此需要先將模型圖像轉(zhuǎn)換成灰度圖像，并進(jìn)行圖像預(yù)處理和濾波等一系列的處理流程，得到完整的模型裂隙區(qū)域信息。

1.1 圖像預(yù)處理

為了增強(qiáng)圖像裂縫區(qū)的顯著性，采用分段線性灰度變換的方法對(duì)圖像進(jìn)行增強(qiáng)處理。該方法先對(duì)圖像灰度進(jìn)行分段，然后在各個(gè)灰度區(qū)間內(nèi)，分別建立灰度線性映射關(guān)系，在增強(qiáng)興趣區(qū)的同時(shí)達(dá)到抑制背景和噪聲的目的。選取模型局部圖像進(jìn)行實(shí)驗(yàn)，原始灰度圖如圖1(a)所示，灰度變換后結(jié)果如圖1(b)所示。

圖1 圖像預(yù)處理

灰度圖像中冗余數(shù)據(jù)較多，不便于進(jìn)行圖像分析，因此需要將灰度圖像二值化。采用Otsu算法進(jìn)行閾值分割，該算法具有較好的適用性和分割效果，其核心是通過(guò)不斷對(duì)閾值迭代調(diào)整，使圖像前景與背景之間的灰度方差達(dá)到最大。設(shè)G(x，y)為二值圖像，g(x，y)為灰度圖，T0為Otsu法分割閾值，則分割數(shù)學(xué)過(guò)程如式(1)所示，二值圖像如圖1(c)所示。

(1)

1.2 多級(jí)連通域?yàn)V波

模型的二值圖像中，仍含有大量形態(tài)多樣，分布無(wú)規(guī)律的噪聲，傳統(tǒng)的數(shù)學(xué)濾波模型無(wú)法有效濾除此類(lèi)噪聲，因而，本文根據(jù)噪聲的分布和形狀特點(diǎn)，設(shè)計(jì)了一種多級(jí)連通域?yàn)V波方法。主要的步驟如下：

(1)搜索二值圖像上每個(gè)非0像素的4-鄰域，得到圖像的若干個(gè)4-鄰域連通域，建立離散像素之間的拓?fù)潢P(guān)系。

(2)基于連通域面積濾波。模型二值圖像上存在著大量面積較小的斑塊噪聲，可通過(guò)比較統(tǒng)計(jì)，設(shè)置合理的面積閾值，濾除面積較小噪聲區(qū)域。

(3)基于連通域外觀比濾波。模型裂隙一般大體呈水平向分布，而連通域的外觀比能反映其在水平方向上的細(xì)長(zhǎng)程度，因此可利用該指數(shù)對(duì)縱向分布的噪聲進(jìn)行濾除。設(shè)gi(x，y)為二值圖像的連通域，n為連通域個(gè)數(shù)，其目標(biāo)外接矩形的長(zhǎng)寬分別為L(zhǎng)i，Wi，該連通域外觀比為ri，設(shè)置濾波閾值為T(mén)r，則濾波過(guò)程如式(2)，(3)所示。

ri=Li/Wi

(2)

(3)

(4)基于連通域形狀因子濾波。由于網(wǎng)格墨跡及材料反光，二值圖像上會(huì)存在邊緣復(fù)雜的噪聲，而裂隙一般分布緊湊且邊緣相對(duì)規(guī)則，因此可以利用連通域的形狀因子進(jìn)行濾波。首先獲取每個(gè)連通域gi(x，y)的周長(zhǎng)Bi和面積Si，并根據(jù)公式(4)計(jì)算形狀因子Fi，設(shè)定合理的閾值Tf，則濾波的數(shù)學(xué)過(guò)程如式(5)所示。

(4)

(5)

(5)特殊噪聲濾除。相似模型表面布設(shè)有大量圓形攝影測(cè)量標(biāo)志點(diǎn)，經(jīng)過(guò)相機(jī)透鏡中心投影后變形為橢圓，當(dāng)這些標(biāo)志點(diǎn)和裂隙連在一起時(shí)，會(huì)難以直接通過(guò)上述三步濾波方法濾除。在盡量保證正直拍攝的情況下，投影橢圓的扁率接近0，因而，可以通過(guò)近似圓形的形狀因子指標(biāo)，即Fi≈1，尋找出獨(dú)立的標(biāo)志點(diǎn)連通域，然后以此為連通域?yàn)V波模板，基于相似度判斷，將混雜在裂隙區(qū)的標(biāo)志點(diǎn)識(shí)別并濾除。

2 標(biāo)志點(diǎn)識(shí)別與定位

2.1 標(biāo)志點(diǎn)分類(lèi)

本次試驗(yàn)中采用的標(biāo)志點(diǎn)包括12位編碼點(diǎn)和非編碼點(diǎn)。其中12位編碼點(diǎn)主要由白色圓形標(biāo)志點(diǎn)和12等分的同心圓環(huán)編碼塊組成，白色圓點(diǎn)域和編碼圓環(huán)內(nèi)外徑之比為3∶7∶12，每一個(gè)30°編碼塊取黑色或白色，分別對(duì)應(yīng)0和1，如圖2(a)所示，因此12個(gè)編碼塊按照一定順序排列便可得到一個(gè)12位二進(jìn)制數(shù)，將該二進(jìn)制數(shù)進(jìn)行移位運(yùn)算，以其對(duì)應(yīng)的最小十進(jìn)制數(shù)值作為編碼值，12位編碼標(biāo)志點(diǎn)通過(guò)編碼塊組合可以包含351個(gè)編碼方案。非編碼點(diǎn)由中心白色圓域和外圍黑色同心圓環(huán)構(gòu)成，其內(nèi)外徑之比為1∶2，如圖2(b)所示。

圖2 三維視覺(jué)測(cè)量標(biāo)志點(diǎn)類(lèi)型

2.2 標(biāo)志點(diǎn)中心定位

為了確定裂縫的空間位置，需要對(duì)標(biāo)志點(diǎn)中心進(jìn)行精確定位。本文根據(jù)標(biāo)志點(diǎn)的幾何特點(diǎn)，設(shè)計(jì)了一種半徑約束下最小二乘橢圓擬合的中心定位方法，該方法的主要思路為：

(1)首先基于連通域面積濾波的方法濾除二值圖像中面積較大或者較小的噪聲。

(2)根據(jù)實(shí)驗(yàn)中圓形標(biāo)志點(diǎn)變形的特點(diǎn)，推導(dǎo)確定標(biāo)志點(diǎn)連通域外觀形狀因子的范圍為0.955

(3)利用形心法計(jì)算每個(gè)標(biāo)志點(diǎn)連通域的中心坐標(biāo)(u0，v0)，并計(jì)算連通域邊界點(diǎn)到中心的平均距離r。

(4)利用canny算子對(duì)模型灰度圖像進(jìn)行邊緣檢測(cè)，并根據(jù)圓形標(biāo)志點(diǎn)內(nèi)外徑之比為1∶2，確定約束半徑為3/2r，濾除連通域中心約束半徑以外的點(diǎn)，以此得到每個(gè)標(biāo)志點(diǎn)的內(nèi)圓邊緣點(diǎn)(xi，yi)。

(5)利用公式(6)對(duì)步驟(4)中得到的邊緣點(diǎn)進(jìn)行最小二乘法橢圓擬合，得到橢圓的圓心坐標(biāo)(xc，yc)，即為標(biāo)志點(diǎn)的中心坐標(biāo)，如圖3(a)所示。

圖3 標(biāo)志點(diǎn)識(shí)別與定位

(6)

式中，f為目標(biāo)函數(shù)；A，B，C，D，E為橢圓方程的5個(gè)系數(shù)；N為邊緣點(diǎn)個(gè)數(shù)。

2.3 編碼標(biāo)志點(diǎn)識(shí)別

編碼標(biāo)志點(diǎn)是進(jìn)行坐標(biāo)空間轉(zhuǎn)換的關(guān)鍵，因此對(duì)其進(jìn)行準(zhǔn)確識(shí)別具有重要的意義。根據(jù)編碼標(biāo)志點(diǎn)的幾何和變形特點(diǎn)，設(shè)計(jì)了一種有效的編碼點(diǎn)識(shí)別方法，該方法的主要步驟包括：

(1)在1.2節(jié)連通域面積濾波的基礎(chǔ)上，根據(jù)編碼標(biāo)志點(diǎn)編碼塊的幾何特點(diǎn)及變形特征，根據(jù)式(4)計(jì)算了編碼塊的外觀形狀因子系數(shù)變化范圍為1.186

(2)根據(jù)2.2節(jié)中所求取各個(gè)標(biāo)志點(diǎn)白色圓域的圓心坐標(biāo)及橢圓方程，以一定的擴(kuò)充系數(shù)σ對(duì)每個(gè)擬合橢圓進(jìn)行擴(kuò)充，并獲取擴(kuò)充橢圓的邊緣點(diǎn)坐標(biāo)及其灰度值，根據(jù)編碼點(diǎn)的設(shè)計(jì)比例，確定橢圓擴(kuò)充系數(shù)為σ=3。

(3)為了消除投影變形的影響，根據(jù)式(7)對(duì)擴(kuò)充橢圓進(jìn)行仿射變換，將其變換成圓，如圖3(c)所示。

(7)

式中，(x，y)和(x’，y’)分別為變換前后邊緣點(diǎn)坐標(biāo)；a，b分別為擴(kuò)充橢圓的長(zhǎng)短半軸；α為橢圓傾角。

(4)按照順時(shí)針順序?qū)A邊緣進(jìn)行采樣，尋找灰度分別為0和1的相鄰點(diǎn)，取其中點(diǎn)作為分界點(diǎn)，并計(jì)算前后相鄰兩個(gè)分界點(diǎn)之間的圓心角及該角內(nèi)包含的12位編碼塊數(shù)量。

(5)根據(jù)圓心角兩側(cè)灰度變化規(guī)律，確定夾角內(nèi)編碼塊的灰度為0或1，據(jù)此得到一個(gè)12位的二進(jìn)制編碼序列，由于采樣起點(diǎn)不同，可能得到不同的二進(jìn)制序列，因而需要對(duì)非0二進(jìn)制序列進(jìn)行移位運(yùn)算，取其對(duì)應(yīng)的最小十進(jìn)制數(shù)值作為該編碼標(biāo)志點(diǎn)的編碼序號(hào)。由于本實(shí)驗(yàn)采用的編碼規(guī)則與XJTUDP近景攝影測(cè)量系統(tǒng)不同，因此編碼值也不相同，但仍是一一對(duì)應(yīng)關(guān)系。點(diǎn)位識(shí)別效果如圖3(d)所示。

3 模型裂隙定位

3.1 畸變校正

普通數(shù)碼相機(jī)在成像的過(guò)程中都會(huì)存在一定程度的幾何失真，即相機(jī)畸變，畸變的存在會(huì)使圖像上裂隙輪廓點(diǎn)及標(biāo)志點(diǎn)的真實(shí)位置產(chǎn)生誤差，因此在進(jìn)行坐標(biāo)空間轉(zhuǎn)換之前需要進(jìn)行畸變校正，建立像點(diǎn)的坐標(biāo)改正模型，如式(8)所示[16]。

(8)

式中，(Δx，Δy)為像點(diǎn)坐標(biāo)改正值；(x，y)為像點(diǎn)實(shí)際量測(cè)坐標(biāo)；(x0，y0)為像主點(diǎn)坐標(biāo)；k1，k2，k3為徑向畸變系數(shù)；p1，p2為偏心畸變系數(shù)；r為像點(diǎn)的畸變半徑。相機(jī)畸變系數(shù)可通過(guò)相機(jī)標(biāo)定得到。

3.2 直接線性變換

經(jīng)過(guò)本文方法檢測(cè)獲取的裂隙輪廓坐標(biāo)和標(biāo)志點(diǎn)坐標(biāo)均為像方坐標(biāo)，無(wú)法直接用于分析研究，需要將像方坐標(biāo)轉(zhuǎn)換為物方坐標(biāo)。為了確定裂隙在模型表面的實(shí)際位置，首先對(duì)模型裂隙圖像上的編碼標(biāo)志點(diǎn)進(jìn)行識(shí)別和中心定位，然后利用近景攝影測(cè)量系統(tǒng)測(cè)出標(biāo)志點(diǎn)的物方坐標(biāo)，進(jìn)而建立裂隙坐標(biāo)空間的直接線性變換模型。

直接線性變換(DLT)不需要事先測(cè)定相機(jī)的內(nèi)外方位元素，而可以直接建立點(diǎn)位的物方坐標(biāo)與像方坐標(biāo)之間的線性變換關(guān)系，因此在非量測(cè)相機(jī)攝影測(cè)量中應(yīng)用較多[17]?？紤]本文的監(jiān)測(cè)目標(biāo)為二維相似材料模型，模型裂隙與標(biāo)志點(diǎn)均布設(shè)在同一平面，因此只需要測(cè)量不少于4個(gè)編碼標(biāo)志點(diǎn)的物方坐標(biāo)，便可以求解出二維DLT單應(yīng)變換矩陣及其變換參數(shù)，進(jìn)而通過(guò)變換模型將裂隙輪廓坐標(biāo)轉(zhuǎn)換到物方空間，具體的變換模型如式(9)所示。

(9)

式中，(xu，yu)為編碼標(biāo)志點(diǎn)去畸變像方坐標(biāo)；(X，Y)為對(duì)應(yīng)標(biāo)志點(diǎn)物方坐標(biāo)；L1～L8為二維DLT變換的8個(gè)參數(shù)。

4 實(shí)驗(yàn)結(jié)果及分析

4.1 模型裂隙檢測(cè)

本次試驗(yàn)對(duì)模擬西部某礦區(qū)開(kāi)采的相似材料模型進(jìn)行裂隙檢測(cè)與定位研究，模型的幾何尺寸為3000mm×1500mm×300mm，正面布設(shè)有大量編碼標(biāo)志點(diǎn)與非編碼標(biāo)志點(diǎn)，試驗(yàn)采用普通數(shù)碼相機(jī)對(duì)模型裂隙進(jìn)行觀測(cè)。選取2期(11.09，11.18)數(shù)據(jù)進(jìn)行試驗(yàn)驗(yàn)證，試驗(yàn)效果如圖4所示，其中圖4(a)為原始裂隙圖像，圖4(b)為裂隙二值圖像，圖4(c)為裂隙檢測(cè)與定位結(jié)果。

圖4 模型裂隙檢測(cè)結(jié)果

由上述實(shí)驗(yàn)結(jié)果可以看出，利用本文方法可以較好地檢測(cè)并提取出模型的裂隙，檢測(cè)結(jié)果與原圖中裂縫位置基本一致。但是，受圖像分辨率以及光照條件的影響，對(duì)于一些細(xì)小裂隙，會(huì)存在一定的“漏檢”現(xiàn)象。

4.2 模型裂隙定位結(jié)果分析

在視覺(jué)測(cè)量中，受攝影方式、圖像分辨率、標(biāo)志點(diǎn)像方定位精度、標(biāo)志點(diǎn)實(shí)測(cè)精度、裂隙像方坐標(biāo)精度等因素的影響，相似材料模型裂隙定位會(huì)存在一定的誤差。受限于現(xiàn)有的技術(shù)手段，實(shí)驗(yàn)中難以直接對(duì)裂隙的位置精確實(shí)測(cè)并進(jìn)行精度分析，而在模型裂隙提取準(zhǔn)確的前提下，其定位精度等于周?chē)蔷幋a點(diǎn)的定位精度，因此本文通過(guò)對(duì)裂隙周?chē)姆蔷幋a標(biāo)志點(diǎn)進(jìn)行定位誤差分析，來(lái)反映裂隙的位置精度。由于XJTUDP近景攝影測(cè)量系統(tǒng)單點(diǎn)測(cè)量精度可達(dá)0.03～0.1mm[18]，完全可以滿(mǎn)足礦山相似材料模型監(jiān)測(cè)的要求，本文將XJTUDP系統(tǒng)測(cè)量的點(diǎn)位坐標(biāo)值作為實(shí)測(cè)值，與本文方法的定位結(jié)果進(jìn)行對(duì)比分析。實(shí)驗(yàn)以119，15，459和349共4個(gè)編碼點(diǎn)作為控制點(diǎn)，并隨機(jī)選取模型圖像中50個(gè)點(diǎn)進(jìn)行定位誤差分析，得到統(tǒng)計(jì)曲線如圖5所示。

圖5 模型非編碼點(diǎn)定位誤差

綜合統(tǒng)計(jì)數(shù)據(jù)，按照公式(10)計(jì)算得到檢查點(diǎn)在X，Y方向上坐標(biāo)值相對(duì)于XJTUDP系統(tǒng)的中誤差分別為：mx=0.840mm，my=0.714mm。XJTUDP系統(tǒng)點(diǎn)位測(cè)量誤差最大為0.1mm，考慮該系統(tǒng)自身誤差，本文方法在X，Y方向上定位誤差小于0.940mm，完全可以滿(mǎn)足一般相似材料模擬實(shí)驗(yàn)中裂隙定位的要求。利用本文方法對(duì)11.18期模型圖像中的裂隙發(fā)育最高點(diǎn)定位，結(jié)果如表1所示，其中裂隙編號(hào)C1，C2，C3等見(jiàn)圖4(c)(11.18)。

(10)

4.3 裂隙檢測(cè)結(jié)果應(yīng)用

由4.1可知，通過(guò)本文的模型裂隙檢測(cè)方法，可以得到只包含裂隙的二值圖像，而每一條裂隙都對(duì)應(yīng)二值圖像上的一個(gè)連通域，對(duì)每個(gè)連通域進(jìn)行邊界遍歷，便可得到每條裂隙所對(duì)應(yīng)的輪廓點(diǎn)像方坐標(biāo)集合；然后利用本文2.2和2.3中方法對(duì)模型圖像上攝影測(cè)量編碼點(diǎn)進(jìn)行識(shí)別與中心坐標(biāo)檢測(cè)，并根據(jù)XJTUDP系統(tǒng)測(cè)得的編碼點(diǎn)實(shí)際坐標(biāo)，解算該張圖像的直接線性變換矩陣；最后通過(guò) DLT變換將裂隙輪廓點(diǎn)的像方坐標(biāo)轉(zhuǎn)換到真實(shí)物方坐標(biāo)空間，此時(shí)裂隙的輪廓尺寸、位置等信息均為實(shí)際值。為了便于模型裂隙的觀察與量測(cè)分析，可以將裂隙輪廓點(diǎn)和開(kāi)切眼的真實(shí)坐標(biāo)導(dǎo)入到CAD中，并通過(guò)閉合曲線對(duì)每一條裂隙的輪廓點(diǎn)連接，便可得到如圖6所示的三期裂隙分布圖。后續(xù)通過(guò)CAD查詢(xún)功能可方便快捷地獲取每條裂隙的尺寸、面積、分布角度等信息，為研究巖層破壞空間、“三帶”發(fā)育規(guī)律等研究提供直觀的依據(jù)。

圖6 模型裂隙分布

5 結(jié) 論

(1)對(duì)基于視覺(jué)測(cè)量的相似材料模型裂隙檢測(cè)與定位方法進(jìn)行了研究，通過(guò)模型圖像預(yù)處理和多級(jí)連通域?yàn)V波實(shí)現(xiàn)了對(duì)裂隙的準(zhǔn)確檢測(cè)與提??；設(shè)計(jì)了有效的編碼點(diǎn)識(shí)別與中心定位方法，并基于DLT變換實(shí)現(xiàn)了對(duì)模型裂隙位置的精確定位。

(2)試驗(yàn)研究表明，在裂隙檢測(cè)準(zhǔn)確的前提下，本文方法在X，Y方向上的裂隙定位誤差小于0.940mm，可以較為準(zhǔn)確高效地確定模型的裂隙發(fā)育位置；同時(shí)，本文方法具有成本低廉、可實(shí)時(shí)連續(xù)觀測(cè)等特點(diǎn)，能夠?yàn)橄嗨撇牧夏Ｐ椭袔r層破壞和 “三帶”發(fā)育規(guī)律研究提供數(shù)據(jù)支持，對(duì)于礦山開(kāi)采沉陷等的物理模擬研究具有一定的實(shí)用價(jià)值。