薛喜成， 楊艷平， 李識博， 楊 康

(西安科技大學(xué)地質(zhì)與環(huán)境學(xué)院，陜西 西安 710054)

0 引 言

采礦相似模擬試驗是將礦山巖層按比例縮小為相似模型，在該模型上進(jìn)行模擬開采[1]，然后通過一定的方法得到模型體的位移量和裂縫發(fā)育情況。目前最常用的方法有百分表方法、千分表方法、透鏡法、三維激光掃描儀法等[2]。百分表法常用于監(jiān)測特殊點[3]，優(yōu)點是精度高、簡單方便，缺點是一個百分表只能觀測單個點僅在一個方向移動時的位移。透鏡法使用經(jīng)緯儀測量[4]，觀測不方便，當(dāng)監(jiān)測點多且密時，工作量也會隨之增大。因此，迫切需要一套新的，可以準(zhǔn)確、高效獲取相似模擬模型表面巖層沉陷數(shù)據(jù)的方法。而數(shù)字散斑相關(guān)方法（DSCM）基本原理是在物體發(fā)生變形之前與物體發(fā)生變形之后其表面散斑點會隨著物體表面移動，通過對比散斑點的位移情況即可確定物體變形。該方法相對自動化程度高，且有非接觸、精度高等優(yōu)點。有限元方法常用于物體表面結(jié)構(gòu)的矩陣分析中，它的主要原理是將連續(xù)體看作多個離散的單元，根據(jù)分塊近似的思想，假定一個插值函數(shù)，近似地表示所有單元之間的規(guī)律。將有限元思想與數(shù)字散斑技術(shù)結(jié)合，并根據(jù)試驗特點改進(jìn)算法，通過匹配不同時期監(jiān)測點，實時計算得到模型形變量。在二十世紀(jì)80年代初國外科學(xué)家首先使用DSCM方法研究了剛體的變形測量，為進(jìn)一步研究提供了可靠的理論依據(jù)。Sutton[5]等提出粗－細(xì)相關(guān)搜索法，使相關(guān)搜索精度達(dá)到亞像素級。在亞像素估計方面，插值方法目前常使用三次樣條插值以及分形插值[6]等。國內(nèi)學(xué)者基于相關(guān)搜索方法，提出比相關(guān)系數(shù)和搜索方法更一般形式，該方法在確保測量精度的前提下，顯著地節(jié)省了時間并提高了速度。目前，DSCM主要用于研究材料的破壞變形，本文嘗試將DSCM應(yīng)用于采礦相似模擬試驗。

1 非量測相機(jī)檢校

相機(jī)檢校的目的是為確定相機(jī)的內(nèi)方位元素（焦距和像主點坐標(biāo)）、徑向畸變誤差和切向畸變誤差，為計算單應(yīng)變換矩陣提供基本參數(shù)[7]。本試驗使用基于稀疏矩陣的光束法相機(jī)檢校法獲得相機(jī)內(nèi)方位元素及鏡頭畸變誤差，該方法是傳統(tǒng)光束法進(jìn)行改良后的一種檢校方法。圖1為室內(nèi)相機(jī)檢校場，表1為相機(jī)檢校結(jié)果。

圖1 室內(nèi)相機(jī)檢校場

表1 相機(jī)檢校結(jié)果

2 圖像預(yù)處理

2.1 噪聲處理

DSCM 相關(guān)匹配是基于灰度的匹配，需將不同時期拍攝到的圖像轉(zhuǎn)換為灰度圖像。相似材料模型中含有大量具有高反射率的云母、沙石，會使圖像中噪聲點過多，反差不明顯[8]。而DSCM的分析技術(shù)[9]涉及數(shù)字圖像散斑點獲取 、圖像相關(guān)系數(shù)計算等，噪聲對測量精度具有很大影響。因此對于實時采集的圖像，要選擇合適的濾波函數(shù)進(jìn)行去噪處理[10]。本試驗采用Daubechies小波濾波器，小波函數(shù)φ(t)由尺度函數(shù)φ(t)加權(quán)組成。

其中k=2–(2N–1)(N=2，3，4,···)，k值不同，權(quán)重gk的值也不同。

關(guān)于數(shù)字散斑位移場測量，金觀昌[11]等已驗證小波減噪方法的效果，該方法使得噪聲信號大幅減少，減噪效果非常明顯，能夠真實反映變形位移場的測試結(jié)果。

2.2 圖像裁剪及各時期圖像上開采進(jìn)度獲取

在基準(zhǔn)圖像上裁取布有散斑點的區(qū)域，該區(qū)域上像點的像素坐標(biāo)加上該區(qū)域左上角在整張影像上的坐標(biāo)即可得到像點在整張影像上的坐標(biāo)。在目標(biāo)圖像上獲取開采進(jìn)度點的坐標(biāo)，將模型表面位移場進(jìn)行分區(qū)處理，圖2為分區(qū)處理示意圖。

圖2 模型分區(qū)處理

3 相關(guān)函數(shù)

圖像灰度分布是物體位移和變形信息的載體，數(shù)字散斑相關(guān)就是通過比對搜索點與基準(zhǔn)點周圍一定區(qū)域灰度矩陣（相關(guān)窗口）的相關(guān)系數(shù)，根據(jù)極值條件尋找最佳匹配，即在目標(biāo)影像中識別出基準(zhǔn)影像中散斑點的匹配點[12]，原理如圖3所示。因此首先要確定表征兩幅圖像相關(guān)窗口相似程度的相關(guān)函數(shù)

圖3 數(shù)字散斑相關(guān)方法計算原理圖

歸一化協(xié)方差相關(guān)函數(shù)通過兩圖像相關(guān)窗口的灰度均方差來實現(xiàn)協(xié)方差相關(guān)函數(shù)的歸一化。

式中：f=f(x，y)——以基準(zhǔn)點為中心的子區(qū)的灰度值；

fm——以基準(zhǔn)點為中心的子區(qū)灰度平均值；

g=g(x+u，y+v)——以目標(biāo)點為中心的子區(qū)的灰度值；

gm——以目標(biāo)點為中心的子區(qū)灰度平均值；

u，v——其水平方向和垂直方向的位移值；

M——子區(qū)域的寬或高[13]。

對歸一化協(xié)方差相關(guān)函數(shù)進(jìn)一步處理，得：

該相關(guān)函數(shù)不僅繼承了原歸一化協(xié)方差相關(guān)函數(shù)的優(yōu)點，還將原歸一化協(xié)方差相關(guān)函數(shù)中的開方運(yùn)算轉(zhuǎn)換為乘法運(yùn)算[14]，很大程度上提高了運(yùn)算速度。

4 基于有限元的監(jiān)測點位移初值獲取

有限元方法將連續(xù)體看作有限個離散單元的集合體，以代替它原來的結(jié)構(gòu)，并在單元上指定結(jié)點[15]（見圖4） 。

圖4 有限元節(jié)點

將有限元思想應(yīng)用于采礦相似模擬試驗；概括起來可分為以下幾個步驟：

1） 結(jié)構(gòu)的離散化。將圖2中位移沉陷區(qū)和幾乎無位移區(qū)分割成71×31個單元，假定該單元集合體滿足一定的位移規(guī)律。取每個單元體的中心點為結(jié)點，在結(jié)點周圍一定區(qū)域內(nèi)逐點搜索，計算并比較每個點的相關(guān)系數(shù)值，得到該結(jié)點的匹配點（即相關(guān)系數(shù)最大值點）。

2）由于相似模擬試驗?zāi)Ｐ褪窍蛳鲁料?，且向左或向右偏移；另外該區(qū)域位移矢量大小自下向上遞減，自中間向兩邊遞減。因此取每個結(jié)點的搜索區(qū)域為以該結(jié)點為頂點，頂角為30°的等腰三角形，取底部最中處結(jié)點的三角形搜索區(qū)的高150個像素，設(shè)置一個遞減閾值，便可得各結(jié)點三角形搜索區(qū)的高。

3）搜索位于位移沉陷區(qū)的結(jié)點。在這些結(jié)點的三角形搜索區(qū)內(nèi)搜索得到匹配點，相減即可得到該結(jié)點的位移。

4）插值函數(shù)。為了能用位移沉陷區(qū)內(nèi)所有結(jié)點位移表征位移沉陷區(qū)內(nèi)所有監(jiān)測點的位移，必須用數(shù)學(xué)公式近似表示結(jié)點位移的規(guī)律，也就是假定結(jié)點位移是結(jié)點坐標(biāo)的某種函數(shù)，即插值函數(shù)。

本試驗插值函數(shù)采用多項式函數(shù)，迭代求解得到插值函數(shù)中的待定系數(shù)。通過這些待定系數(shù)完全可以描述散斑圖子區(qū)的變形規(guī)律，從而求得監(jiān)測點的位移初值。一般選擇多項式：

對于子區(qū)中任意一點Q(x，y)，變形后為Q*(x*，y*)，位移u*，v*，可得：

由此可得Q點變形后Q*的位置：

5）從圖2中的裁剪區(qū)左上角點向右、向下每隔m個像素取點作為監(jiān)測點，m可設(shè)置為閾值，將監(jiān)測點分為兩部分，一部分為位于位移沉陷區(qū)的監(jiān)測點，位移初值由第4步得到的插值函數(shù)插值求得；一部分為位于幾乎無位移區(qū)的監(jiān)測點，位移初值為0，監(jiān)測點在裁剪區(qū)的坐標(biāo)加上裁剪區(qū)左上角在整張影像上的坐標(biāo)即為監(jiān)測點在整張影像上的像點坐標(biāo)。至此，所有監(jiān)測點在基準(zhǔn)影像上的像點坐標(biāo)以及它相對于目標(biāo)影像位移的初值已求得。

5 數(shù)字散斑相關(guān)搜索方法

在匹配點坐標(biāo)初值給定的基礎(chǔ)上，還需要使用搜索算法獲得匹配點坐標(biāo)的整像素級精確值，本試驗使用以下搜索算法。

5.1 新三步搜索法（NTSS）

三步搜索法[16]（three-step search，TSS）常用于視頻處理領(lǐng)域，既保持了較高的精度，同時大大提高了運(yùn)算速率。但是TSS三階段搜索方法存在影響準(zhǔn)確性的缺點，不易進(jìn)入局部最優(yōu)區(qū)域。因此基于三步法提出了新三步搜索法（new three-step search，NTSS）（圖5）。在三步搜索法第一步的基礎(chǔ)上增加步長為1的小范圍搜索，此時共計 17 個搜索點。若最大值點在小范圍內(nèi)，以最大值點為初值再做一次小范圍搜索（重復(fù)點可不再計算），得到最終匹配結(jié)果，此時算法結(jié)束，否則搜索步驟同三步法。

圖5 新三步法（NTSS）搜索示意圖

新三步法在保持了三步法運(yùn)算速率的同時，新增小范圍搜索可以適應(yīng)位移較小的情況，大大降低了搜索結(jié)果陷入局部最優(yōu)的概率。

5.2 菱形搜索法

菱形搜索法[17]（DS）以大小兩個菱形模板協(xié)同搜索，大菱形模板有9個像素點，小菱形模板有5個像素點。給定一個初值作為大菱形搜索中心，計算大菱形模板上9個點的相關(guān)系數(shù)。如果最大值點位于大菱形中心，此時計算包括中心點在內(nèi)的小菱形上5個點的相關(guān)系數(shù)，比較得到的最大值點即為最終匹配結(jié)果，否則繼續(xù)按照大菱形模板搜索，直至最大值點位于大菱形中心。其中，大菱形模板搜索又分兩種情況，若最大值點位于菱形四角，只需再計算5個點的相關(guān)系數(shù)，其余5點為重復(fù)點；若最大值點位于菱形四邊，則需再計算3個點的相關(guān)系數(shù)。搜索過程如圖6所示。

圖6 菱形（DS）搜索示意圖

由于菱形搜索算法具有循環(huán)持續(xù)搜索的優(yōu)勢，適合大位移搜索，而新三步法由于步長限制，適合小位移搜索。因此，本試驗對于位移沉陷區(qū)的監(jiān)測點使用菱形搜索法，幾乎無位移區(qū)的監(jiān)測點使用新三步搜索法。

5.3 亞像素估計算法

常見的亞像素估計算法有灰度插值法、曲面擬合法等[18]。梁園[19]等已對數(shù)字散斑相關(guān)亞像素位移測量的各種方法進(jìn)行試驗對比與分析。灰度插值法相對其他亞像素估計法更簡單易懂，且結(jié)果較為穩(wěn)定。

本試驗采用灰度插值法的雙三次多項式插值。該方法插值后圖像分布比較光滑，散斑質(zhì)量較好時，分辨率可以達(dá)到0.03個像素。下式為對圖像進(jìn)行卷積運(yùn)算的核函數(shù)[20]：

其中d為插值點與周圍最鄰近的16個點的距離分量。圖像中任一點(x，y)的灰度值為可以表達(dá)為：

式中：g(i，j)——16個最鄰近整數(shù)像素的灰度值；

dxi和dyi——插值點與16個整像素點的距離在x和y方向上的分量。

6 單應(yīng)矩陣獲取監(jiān)測點物方坐標(biāo)

在需要得到模型位移的物方坐標(biāo)時，可使用單應(yīng)變換原理，將像點坐標(biāo)映射到物方平面。兩個不同角度拍攝的兩幅圖像之間存在著映射關(guān)系，用3×3矩陣表示，即單應(yīng)變換矩陣H[21]。本文首次嘗試用被攝物體的二維表面來代替兩幅影像中的其中一幅。已知控制點物方坐標(biāo)X和控制點像點坐標(biāo)x，根據(jù)公式（9）即可求得像平面到物平面的單應(yīng)變換矩陣。

單應(yīng)變換求解物方平面坐標(biāo)的具體步驟如下：

1） 相機(jī)檢校。通過基于稀疏矩陣的光束法相機(jī)檢校方法獲取相機(jī)的內(nèi)方位元素及畸變參數(shù)，為求解H提供已知值。

2） 控制點量測。利用高精度全站儀逐點測量控制點的物方坐標(biāo)，通過刺點得到控制點的像點坐標(biāo)。

3） 監(jiān)測點物方坐標(biāo)解算。求出每張影像的單應(yīng)變換矩陣H，根據(jù)數(shù)字散斑相關(guān)方法得出的監(jiān)測點像方坐標(biāo)，完成監(jiān)測點物方坐標(biāo)求解。

7 試驗結(jié)果分析

7.1 模擬散斑圖系統(tǒng)驗證

為了驗證本方法的可行性，首先使用模擬散斑圖進(jìn)行驗證，模擬散斑圖是用數(shù)值方法生成兩張或多張具有一定變形量的散斑圖。本試驗使用兩組具有不同變形情況的散斑圖像，一組具有豎直向下的均勻位移，另一組具有豎直斜向下的均勻位移（與豎直方向的偏角小于30°）。具有豎直向下均勻位移的散斑圖如圖7(a)，(b)所示，具有豎直斜向下均勻位移的散斑圖如圖7(c)，(d)所示，圖8(a)，(b)為模擬散斑圖位移矢量，圖9為具有不同位移矢量的模擬散斑圖精度結(jié)果。

圖7 模擬散斑圖

圖8 位移矢量

圖9 不同位移模擬散斑圖測量精度

由模擬結(jié)果可以看出，該方法位對模擬散斑圖的位移測量結(jié)果精度較高，分辨率達(dá)到0.01像素。而實際散斑圖位移變形更為復(fù)雜，計算量更大，還需使用實際散斑圖對該方法進(jìn)行系統(tǒng)驗證。

7.2 實際散斑圖系統(tǒng)驗證

基于上述理論，利用計算機(jī)語言編程來實現(xiàn)位移量的計算，表2為從第40次開挖到第77次開挖煤礦采空區(qū)上層覆巖的少部分位移數(shù)據(jù)。

表2 第40次開挖到第77次開挖部分位移數(shù)據(jù)像素

為了方便從計算結(jié)果中總結(jié)出上覆巖層的移動規(guī)律，需要將計算出的各點位移變化量用坐標(biāo)曲線形式表示出來，在表達(dá)過程中，位移變化量用小箭頭表示。其中，子區(qū)域運(yùn)動方向為小箭頭所指的方向，小箭頭的長短表示子區(qū)域變形程度。圖10(a)到(d)分別表示第43次開挖工作面，第50次開挖工作面，第54次開挖工作面以及第77次開挖工作面的采空區(qū)上覆巖層變形情況。

圖10 不同時期工作面位移

為了驗證本系統(tǒng)的適應(yīng)度和容錯率，以第40次開挖，開挖長度為151 cm的影像為基準(zhǔn)影像，以第41次到最后一次開挖的影像作為目標(biāo)影像，計算目標(biāo)影像的位移矢量及位移場。計算得到第43次、第50次、第54次、第77次開挖，開挖進(jìn)度分別為165 cm、196 cm、210 cm、283 cm 的目標(biāo)影像的位移場，所得到的工作面物方位移分辨率可達(dá)0.001 m以上。結(jié)果如圖11所示。

本實驗以第40次開挖工作面為基準(zhǔn)影像，即第40次開挖進(jìn)度點為臨界點，第40次開挖進(jìn)度點之前的模型開采工作面應(yīng)相對基準(zhǔn)點上移，而第40次開挖之后的工作面則相對基準(zhǔn)點下移，由圖11可以觀察到本次實驗所得的應(yīng)變結(jié)果與實際的礦山開采過程中上層巖體的形變及位移情況相符。

圖11 不同時期開挖位移場

峰值信噪比（peak signal noise ratio，PSNR）常用來反映視頻處理的精度，一般情況下PSNR≥40為良好，20≤PSNR＜30為一般，PSNR＜20為差。本文引入該指標(biāo)來反映各開挖階段試驗結(jié)果的精度。將計算得到的位移矢量結(jié)果加到基準(zhǔn)影像上得到重構(gòu)影像，計算重構(gòu)影像與目標(biāo)影像的峰值信噪比，判斷該值是否符合精度要求。圖12為各個開挖階段影像的峰值信噪比?？梢钥闯鲭S著開挖的不斷進(jìn)行，模型表面散斑點位移矢量持續(xù)增大，模型表面破壞變形加劇，峰值信噪比也呈下降趨勢，即隨著模型表面破壞變形的加劇，本試驗方法的測量精度呈下降趨勢，但依舊可以反映礦山巖層相似模擬模型表面的沉陷特點與位移趨勢，總體精度符合試驗要求。

圖12 實際散斑圖測量精度

8 結(jié)束語

1）使用白灰和墨水對采礦相似模擬試驗?zāi)Ｐ瓦M(jìn)行散斑處理，制作方法簡單，且散斑圖案清晰，反差明顯，滿足試驗的要求。

2）位移場的精度依賴于位移初值的選取。使用有限元插值法獲取監(jiān)測點位移的初值，有限元劃分越密，位移初值的精度也就越高，相應(yīng)的計算量也會越大。根據(jù)模型沉陷特點及有限元位置設(shè)計的相似三角形搜索區(qū)域可減少大部分計算量，且精度保持不變。

3）對模型表面進(jìn)行分區(qū)處理，根據(jù)不同區(qū)域位移大小的不同，選擇相應(yīng)的搜索算法。選擇精度較高的灰度插值法進(jìn)行亞像素估計，試驗結(jié)果表明，未開采區(qū)上方的位移趨近于零，新開采區(qū)上方呈現(xiàn)金字塔形斷裂式位移，隨著開挖的進(jìn)行，已開挖區(qū)位移也會持續(xù)增大。

4）根據(jù)模型沉陷特點改進(jìn)的 DSCM 技術(shù)用于相似模擬模型的變形觀測是可行的，且具有非接觸、簡單易行、信息量大、精度高等優(yōu)點。