基于數字散斑相關方法的礦山相似模擬位移監測

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

(西安科技大學地質與環境學院，陜西 西安 710054)

0 引 言

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

1 非量測相機檢校

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

圖1 室內相機檢校場

表1 相機檢校結果

2 圖像預處理

2.1 噪聲處理

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

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

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

2.2 圖像裁剪及各時期圖像上開采進度獲取

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

圖2 模型分區處理

3 相關函數

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

圖3 數字散斑相關方法計算原理圖

歸一化協方差相關函數通過兩圖像相關窗口的灰度均方差來實現協方差相關函數的歸一化。

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

fm——以基準點為中心的子區灰度平均值；

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

gm——以目標點為中心的子區灰度平均值；

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

M——子區域的寬或高[13]。

對歸一化協方差相關函數進一步處理，得：

該相關函數不僅繼承了原歸一化協方差相關函數的優點，還將原歸一化協方差相關函數中的開方運算轉換為乘法運算[14]，很大程度上提高了運算速度。

4 基于有限元的監測點位移初值獲取

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

圖4 有限元節點

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

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

2）由于相似模擬試驗模型是向下沉陷，且向左或向右偏移；另外該區域位移矢量大小自下向上遞減，自中間向兩邊遞減。因此取每個結點的搜索區域為以該結點為頂點，頂角為30°的等腰三角形，取底部最中處結點的三角形搜索區的高150個像素，設置一個遞減閾值，便可得各結點三角形搜索區的高。

3）搜索位于位移沉陷區的結點。在這些結點的三角形搜索區內搜索得到匹配點，相減即可得到該結點的位移。

4）插值函數。為了能用位移沉陷區內所有結點位移表征位移沉陷區內所有監測點的位移，必須用數學公式近似表示結點位移的規律，也就是假定結點位移是結點坐標的某種函數，即插值函數。

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

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

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

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

5 數字散斑相關搜索方法

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

5.1 新三步搜索法（NTSS）

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

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

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

5.2 菱形搜索法

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

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

由于菱形搜索算法具有循環持續搜索的優勢，適合大位移搜索，而新三步法由于步長限制，適合小位移搜索。因此，本試驗對于位移沉陷區的監測點使用菱形搜索法，幾乎無位移區的監測點使用新三步搜索法。

5.3 亞像素估計算法

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

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

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

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

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

6 單應矩陣獲取監測點物方坐標

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

單應變換求解物方平面坐標的具體步驟如下：

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

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

3） 監測點物方坐標解算。求出每張影像的單應變換矩陣H，根據數字散斑相關方法得出的監測點像方坐標，完成監測點物方坐標求解。

7 試驗結果分析

7.1 模擬散斑圖系統驗證

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

圖7 模擬散斑圖

圖8 位移矢量

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

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

7.2 實際散斑圖系統驗證

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

表2 第40次開挖到第77次開挖部分位移數據像素

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

圖10 不同時期工作面位移

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

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

圖11 不同時期開挖位移場

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

圖12 實際散斑圖測量精度

8 結束語

1）使用白灰和墨水對采礦相似模擬試驗模型進行散斑處理，制作方法簡單，且散斑圖案清晰，反差明顯，滿足試驗的要求。

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

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

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