基于單軌道OT DInSAR的礦山時序三維多量級形變重建

張樂樂 劉 波 賈世利 常英豪 王佳冰 魏增明

1 內蒙古工業大學機械工程學院，呼和浩特市愛民街49號，010051 2 內蒙古煤炭地質勘查(集團)測繪院有限公司，呼和浩特市展東路25號，010052

礦山地下煤炭資源開采會造成土地塌陷及采空區積水，進而導致耕地面積減小和農作物產量受損，地表房屋、建筑物及交通運輸道路扭曲變形，嚴重時還會導致山體開裂、崩塌、滑坡、泥石流等重大自然災害。合成孔徑雷達差分干涉測量(differential interferometric synthetic aperture radar，DInSAR)技術是基于面的大范圍測量，具有空間分辨率高、覆蓋范圍廣、成本低、快速準確及大尺度連續覆蓋等優點[1-2]，能獲取地表雷達視線向(light of sight，LOS)的高精度小量級形變。但礦山形變往往具有沉降速度快、形變梯度大、多量級和多方向等特點，而DInSAR技術僅能監測地表發生的緩慢、小量級形變，因此監測結果不能反映礦山地表真實形變。基于SAR影像的偏移量跟蹤(offset-tracking，OT)技術[3]具有不受云霧遮擋、抗噪聲能力強、不受形變梯度制約等特點，能有效監測大范圍、大梯度的地表形變，但其監測精度通常在1/10～1/30像元，并不能有效監測下沉盆地邊緣附近的小量級形變，在精確識別下沉盆地邊緣時會受到影響。因此，有學者聯合DInSAR和OT(簡稱OT DInSAR)技術分別獲取了礦山地表LOS向不同量級的形變[4-8]。但OT DInSAR技術僅能提取礦山的LOS向非連續多量級形變，不能提取礦山連續的三維多量級形變，而傳統的三維形變監測方法對SAR影像數量和獲取時間有較高要求[9-12]。基于單軌道InSAR影像提取礦山三維形變能很好地解決這一問題，但該方法目前僅實現了礦山三維非連續大量級形變提取且水平方向位移存在方向偏差的問題[13-16]。為此，本文聯合OT DInSAR技術、多量級形變先驗融合模型及概率積分法，提出一種基于單軌道OT DInSAR的礦山時序三維多量級形變重建方法，以獲取高精度、完整和真實的礦山地表形變信息。

1 礦山三維多量級形變重建模型

1.1 LOS向多量級形變先驗融合模型

基于礦山地理信息和InSAR影像信息的形變先驗知識，可構建礦山地表多量級形變融合模型，對DInSAR技術和OT技術獲取的大量級和小量級的LOS向形變進行礦山多量級形變融合。基于形變先驗知識的融合模型可表示為[7]：

1.2 三維多量級形變重建模型

根據隨機介質理論推導出的水平移動與傾斜成正比的關系，對于相同的采深H，比例系數是一個常數，則可得水平位移與沉降梯度的關系為：

(2)

式中，b為水平移動系數，H為采深,β為主要影響角，這些參數均可通過采礦資料獲取。P(x,H)為沿x方向的單元梯度，P(y,H)為沿y方向的單元梯度。

水平方向真實的疊加開采工作面的礦山地表形變網格如圖1(a)所示，圖中給出了開采工作面的位置和觀測線AA′方向的5個形變點，橫坐標對應東西向，縱坐標對應南北向，i和j代表形變圖的行和列，用來表示形變點的地理位置，i

圖1 礦山三維地表形變曲線Fig.1 Curves of three-dimensional multi-level mining deformation

利用礦山地表監測的水平移動數據、下沉數據、采深H及采空區分布，基于概率積分法模型可求取水平移動系數b和主要影響角正切tanβ。Yang等[15-16]基于單軌道InSAR數據利用李志偉等[13]的方法提取了同一礦山的三維形變，但測量結果中水平位移的正向和負向移動出現顛倒的情況，造成嚴重的方向偏差，大大降低了水平向的測量精度。本文以圖1的幾何圖為參考，給出修正后的礦山地表形變點(i,j)在水平向的位移，具體公式為：

(3)

式中，μE=bH/(tanβ·ΔE)，μN=bH/(tanβ·ΔN)；dW(i,j)為垂直向形變，dE(i,j)為東西向形變，dN(i,j)為南北向形變。

根據SAR傳感器的成像幾何原理，LOS向形變可分解為三維方向的投影，具體公式為：

dLOS(i,j)=dW(i,j)cosθ-dN(i,j)·

(4)

式中，α為SAR傳感器的軌道方位角，θ為SAR傳感器發射電磁波的入射角。將式(3)代入式(4)得到：

dLOS(i,j)=[L1(i,j)L2(i,j)L3(i,j)]·

[dW(i,j)dW(i,j+1)dW(i+1,j)]T

(5)

其中，

(6)

由此可得每個像元的垂直向形變，然后代入式(3)可得東西向和南北向的形變。礦山三維多量級形變重建的步驟(圖2)主要為：

圖2 礦山三維多量級形變解算流程Fig.2 Calculating flowchart of three-dimensional multi-level mining deformation

1)輸入覆蓋研究區的高分辨率InSAR影像和數字高程模型(digital elevation model，DEM)數據。

4)將礦山地表形變點(i,j)對應的LOS向形變值DLOS(i,j)、α、θ、b、H及tanβ代入式(3)和式(4)，利用式(5)和式(6)進行回代處理可得dW(i,j)。

5)將dW(i,j)代入式(3)可得dN(i,j)和dE(i,j)，最終得到基于單軌OT DInSAR重建的礦山三維多量級形變監測結果。

2 研究區分析

2.1 區域概況

大柳塔煤礦地處內蒙古自治區與陜西省的交界處，地理位置為39°13′53″～39°21′32″N、110°12′23″～110°22′54″E。井田位于黃土高原北側和毛烏素沙漠東南側，地勢中部高、東西低，南低北高，平均海拔為1 200 m，礦山內大部分地貌屬風沙堆積地貌，植被稀少。

研究區為52303工作面，該區域煤層儲量大、埋藏深度淺且煤炭質量優良，工作面及觀測線地理位置如圖3所示，圖中包含52303工作面和52304工作面的地理位置及沿工作面走向(2′2″)和傾向(1′1″)布設的地表移動觀測站位置，研究區煤層傾角小于2°，觀測站下方煤層傾角近似于水平煤層。52303工作面開采方向為自東向西沿煤層走向回采，工作面為刀把式，采煤方法為走向長壁式采煤法，全部垮落法頂板管理。設計采高為6.5 m，52303工作面面寬為301 m，走向推進長度為4 443.3 m，開采始于2012-12-06，止于2014-03，水平移動系數b=0.32，平均采深H=230 m，主要影響角正切tanβ=1.24。

圖3 工作面位置Fig.3 Location of the workface

2.2 InSAR實驗數據和DEM數據

選用TerraSAR-X(TSX)雷達衛星的單視斜距復影像數據產品，雷達波入射角為42.43°，軌道方位角為189.53°，像元尺寸為0.91×0.86。表1為TSX影像參數，為去除地形相位的影響，DEM采用SRTM的30 m分辨率數據產品。

3 實驗結果及分析

3.1 OT DInSAR處理

利用OT DInSAR技術獲取52303工作面開采造成的LOS向不同量級形變。在DInSAR數據處理中，距離向和方位向采用3×5的多視因子來消除斑點噪聲，網格尺寸約為5 m。為提高SAR影像對的相干性，濾波方法選取Goldstein濾波，由于X波段的SAR影像對去相關效應高度敏感，故濾波窗口設為32×32。解纏方法采用最小費用流法，解纏閾值設為0.3，對相干性大于0.3的點進行解纏處理，反之則舍棄。在OT數據處理中，為獲取與DInSAR處理結果相同的像素尺寸，多視因子同樣采用3×5；綜合考慮計算的性能和效率，互相關窗口大小取128×128，過采樣因子取4。

圖5為52303工作面開采造成的礦山地表LOS向多量級形變，其中圖5(a)～5(f)為DInSAR測量的小量級形變與OT測量的大量級形變的融合結果。可以看出，多量級形變不連續且存在較多空值點，這些點多是由DInSAR測量結果中的失相干點和OT測量結果中低精度測量點構成。圖5(g)～5(l)為利用IDW估計后的LOS向連續多量級形變測量結果。可以看出，沉降的方向與該工作面由東向西的開采方向一致，沉降范圍與開采周期直接相關。但由于LOS向的測量結果只能反映雷達視線向的形變方向和大小，不能精確獲取三維方向的多量級形變，基于LOS向的多量級形變還需重建礦山三維多量級形變。

3.2 礦山三維多量級形變重建結果分析

圖6為52303工作面在2012-12-13～2013-06-07開采造成的時序三維形變，其中2012-12-13～2013-01-04垂直向的最大沉降值為2.31 m，東西向及南北向的最大形變值分別為3.3 m和2.27 m；2013-01-04～02-06垂直向的最大沉降值為3.73 m，東西向及南北向的最大形變值分別為2.24 m和6.42 m；2013-02-06～03-22垂直向的最大沉降值為3.79 m，東西向及南北向的最大形變值分別為1.86 m和5.74 m；2013-03-22～04-02垂直向的最大沉降值為2.75 m，東西向及南北向的最大形變值分別為1.65 m和4.23 m；2013-04-02～05-05垂直向的最大沉降值為4.74 m，東西向及南北向的最大形變值分別為2.55 m和5.59 m；2013-05-05～06-07垂直向的最大沉降值為2.31 m，東西向及南北向的最大形變值分別為2.33 m和2.95 m。針對該工作面的6組觀測周期內，最大下沉值為4.74 m，最大下沉速率為0.144 m/d。

(a)和(g)：2012-12-13～2013-01-04；(b)和(h)：2013-01-04～2013-02-06；(c)和(i)：2013-02-06～2013-03-22；(d)和(j)：2013-03-22～2013-04-02；(e)和(k)：2013-04-02～2013-05-05；(f)和(l)：2013-05-05～2013-06-07圖5 2012-12-30～2013-06-07大柳塔礦區LOS向多量級形變融合結果Fig.5 Multi-level deformation fusion results in Daliuta mining area from 2012-12-13 to 2013-06-07

分析每組影像對的三維多量級形變最大值后發現，2013-01-04～05-05南北向的形變最大值不符合開采沉陷規律，出現了失真。造成這種現象的原因是由于地表塌陷引起了地裂縫，在形變重建算法中計算水平向移動形變時利用了相鄰兩點的垂直向形變梯度(式(3))，而嚴重的地裂縫會產生過大的非線性垂直向形變梯度，進而造成水平向的形變失真。這樣的失真點多發生于沉降漏斗底部和塌陷式沉降的邊緣，但失真點的數量是有限的，通過對形變圖的分析可以人為地識別出這些失真點，不影響對水平向地表形變的分析。

圖7(a)為52303工作面在2012-12-13～2013-06-07垂直向的累積沉降量，隨著時間的推移，沉降區域的面積逐漸增大，最大沉降量為5.1 m，最大沉降速度為0.066 m/d。綜合開采條件下，由于工作面快速推進，使得煤層上覆巖層各層次下沉速度均加快，相對懸空的時間減少，移動變形集中，引起嚴重的地表多量級形變，對當地環境造成嚴重破壞。圖7(b)為52303工作面垂直向形變的三維顯示結果，能夠直觀反映該工作面地表每個點的累積沉降情況，便于對研究區地表損壞程度進行評估與治理。為驗證52303工作面三維多量級形變監測結果的精度，與同期地面在觀測線2′2″和1′1″的水準監測數據進行對比。由于實際情況中很難獲取同步的InSAR與水準數據，其中高程水準數據獲取時間為2012-08-19～2013-03-19，水平位移水準數據獲取時間為2012-03-20～2013-03-28，而52303工作面開采時間為2012-12-06，本文假設2012-12-06之前52303工作面地表不發生地表位移，因此水準數據與InSAR影像獲取時間分別在垂直向和水平向有3 d和6 d的時間間隔。

用本文方法計算水準觀測線1′1″垂直向、東西向及南北向形變，并與水準測量結果進行對比，結果如圖8(a)～8(c)所示。通過對比可知，該觀測線垂直向、東西向及南北向的多量級形變與水準測量結果的形變趨勢較為一致，RMSE分別為0.188 0 m、0.224 3 m和0.207 4 m。沉降漏斗中心有部分觀測點與水準測量結果相比整體偏大，主要原因有：1)由于水準測量數據和SAR影像獲取時間存在6 d的時間間隔，該監測期間是礦山的主要沉降階段，礦山形變速率快、形變梯度大，導致與水準測量結果偏差大。2)沉降中心的形變是由OT技術獲取的，該技術雖可獲取大量級形變，但精度比DInSAR技術低，故會導致大量級形變點存在誤差。用本文方法計算水準觀測線2′2″垂直向、東西向及南北向形變，并與水準測量結果進行對比，結果如圖8(d)～8(f)所示。可以看出，觀測線2′2″上垂直向、東西向及南北向的多量級形變與水準測量結果的形變趨勢較為一致，RMSE分別為0.070 9 m、0.134 6 m和0.081 6 m。水準觀測線2′2″針對52303工作面相鄰的52304工作面布設，在觀測周期內只能監測到沉降邊緣的小量級形變，在實際情況中，由于52304工作面在2013-03開采完畢，前期開采期間勢必會對相鄰的52303工作面上方的地表造成影響，會有小量級甚至中量級形變發生，因此本文方法重建的垂直向、東西向和南北向的多量級形變監測結果的RMSE要小于0.188 0 m、0.224 3 m和0.207 4 m，滿足礦山三維多量級形變監測的精度要求。

(a)～(c):2012-12-13～2013-01-04;(d)～(f):2013-01-04～02-06;(g)～(i):2013-02-06～03-22;(j)～(l):2013-03-22～04-02;(m)～(o):2013-04-02～05-05; (p)～(r): 2013-05-05～06-07; 第1列為垂直向形變，第2列為東西向形變，第3列為南北向形變圖6 2012-12-13～2013-06-07礦山時序三維多量級形變監測結果Fig.6 Three-dimensional multi-level mining deformation from 2012-12-13 to 2013-06-07

圖7 2012-12～2013-06累積形變監測結果Fig.7 The cumulative deformation from December 2012 to June 2013

圖8 本文方法監測結果與水準測量結果對比Fig.8 Comparison of the deformation obtained by proposed method and leveling

4 結 語

針對礦山地表三維多量級形變獲取難的問題，本文提出一種礦山三維多量級形變重建方法。利用該方法對研究區2012-12-13～2013-06-07期間7景X波段的TerraSAR-X數據進行時序處理，獲取大柳塔礦山52303工作面的時序三維多量級形變，與水準測量結果對比可知，垂直向、東西向及南北向的RMSE分別小于0.188 0 m、0.224 3 m和0.207 4 m，監測結果與水準測量數據形變趨勢一致。本文提出的方法可滿足監測礦山沉降漏斗邊緣地帶cm級到沉降中心m級的多量級形變監測需求，將為礦山三維多量級形變測量提供一種大范圍、高分辨率、高精度、三維方向、低成本的監測技術，同時為InSAR技術在礦山災害監測應用的推廣提供技術保障。