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

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

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

1 礦山三維多量級(jí)形變重建模型

1.1 LOS向多量級(jí)形變先驗(yàn)融合模型

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

1.2 三維多量級(jí)形變重建模型

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

(2)

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

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

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

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

(3)

式中，μE=bH/(tanβ·ΔE)，μN(yù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)

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

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

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

4)將礦山地表形變點(diǎn)(i,j)對(duì)應(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重建的礦山三維多量級(jí)形變監(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ū)域煤層儲(chǔ)量大、埋藏深度淺且煤炭質(zhì)量優(yōu)良，工作面及觀測線地理位置如圖3所示，圖中包含52303工作面和52304工作面的地理位置及沿工作面走向(2′2″)和傾向(1′1″)布設(shè)的地表移動(dòng)觀測站位置，研究區(qū)煤層傾角小于2°，觀測站下方煤層傾角近似于水平煤層。52303工作面開采方向?yàn)樽詵|向西沿煤層走向回采，工作面為刀把式，采煤方法為走向長壁式采煤法，全部垮落法頂板管理。設(shè)計(jì)采高為6.5 m，52303工作面面寬為301 m，走向推進(jìn)長度為4 443.3 m，開采始于2012-12-06，止于2014-03，水平移動(dòng)系數(shù)b=0.32，平均采深H=230 m，主要影響角正切tanβ=1.24。

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

2.2 InSAR實(shí)驗(yàn)數(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 實(shí)驗(yàn)結(jié)果及分析

3.1 OT DInSAR處理

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

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

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

圖6為52303工作面在2012-12-13～2013-06-07開采造成的時(shí)序三維形變，其中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。針對(duì)該工作面的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向多量級(jí)形變?nèi)诤辖Y(jié)果Fig.5 Multi-level deformation fusion results in Daliuta mining area from 2012-12-13 to 2013-06-07

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

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

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

4 結(jié) 語

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