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

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

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

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

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

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

基于礦山地理信息和InSAR影像信息的形變先驗知識，可構(gòu)建礦山地表多量級形變?nèi)诤夏Ｐ?，對DInSAR技術(shù)和OT技術(shù)獲取的大量級和小量級的LOS向形變進(jìn)行礦山多量級形變?nèi)诤??；谛巫兿闰炛R的融合模型可表示為[7]：

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

根據(jù)隨機介質(zhì)理論推導(dǎo)出的水平移動與傾斜成正比的關(guān)系，對于相同的采深H，比例系數(shù)是一個常數(shù)，則可得水平位移與沉降梯度的關(guān)系為：

(2)

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

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

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

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

(3)

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

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

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

(4)

式中，α為SAR傳感器的軌道方位角，θ為SAR傳感器發(fā)射電磁波的入射角。將式(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)輸入覆蓋研究區(qū)的高分辨率InSAR影像和數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)。

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

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

2 研究區(qū)分析

2.1 區(qū)域概況

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

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

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

2.2 InSAR實驗數(shù)據(jù)和DEM數(shù)據(jù)

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

3 實驗結(jié)果及分析

3.1 OT DInSAR處理

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

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

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

圖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組觀測周期內(nèi)，最大下沉值為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大柳塔礦區(qū)LOS向多量級形變?nèi)诤辖Y(jié)果Fig.5 Multi-level deformation fusion results in Daliuta mining area from 2012-12-13 to 2013-06-07

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

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

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

(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礦山時序三維多量級形變監(jiān)測結(jié)果Fig.6 Three-dimensional multi-level mining deformation from 2012-12-13 to 2013-06-07

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

圖8 本文方法監(jiān)測結(jié)果與水準(zhǔn)測量結(jié)果對比Fig.8 Comparison of the deformation obtained by proposed method and leveling

4 結(jié) 語

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