基于沉陷對稱特征的近水平煤層開采InSAR三維位移反演模型

湯伏全，董龍凱，王宗良，黃景偲

(西安科技大學 測繪科學與技術學院，陜西 西安 710054)

InSAR技術已成為煤礦開采沉陷變形監測的重要手段[1]。煤礦開采沉陷引起的地表三維位移包括垂直下沉、沿地表移動盆地走向主斷面及傾向主斷面的水平位移(或轉換為南北向及東西向水平位移)，而InSAR影像的干涉處理多采用一種SAR數據，如升軌數據或降軌數據，主要獲得地表點在雷達視線方向(LOS)的一維位移量，這種InSAR LOS向位移是煤礦地表點下沉和水平向位移沿雷達視線方向投影后的疊加結果[2-4]。其中，下沉方向始終垂直向下，水平位移方向指向采空區中央，移動盆地內地表點的移動方向隨位置而變化，因而在InSAR監測結果中通過單個的LOS 向一維位移，并不能直接分解出地表點在垂直與水平方向上的三維位移值。為此，本文根據InSAR LOS向位移分解的幾何模型及近水平煤層開采地表移動的基本特征[5]，提出一種基于開采沉陷對稱原理將InSAR LOS向位移分解地表三維位移的新方法。

1 InSAR LOS向位移與三維位移的幾何關系

地表點沿雷達視線方向的位移(LOS向位移)由該點在東西向、南北向和垂直向三維位移在視線方向上的投影疊加而成,它們之間的幾何關系如圖1所示[6]。

圖1中θ為衛星入射方位角;αh為衛星飛行方向與北方向的夾角(順時針方向為正)[7]。由圖1的空間幾何關系可得地表點沿雷達視線方向的(LOS)位移dLOS與三維位移的關系式

dLOS=Wcosθ-sinθ[UNcos(αh-

3π/2)+UEsin(αh-3π/2)](1)

式中，W,UN,UE分別是地面點在垂直、南北、東西3個方向上的位移。LOS位移因地表點真實的三維位移引起并由InSAR處理得到，其方向指向衛星為正[7]。

圖1 LOS向位移與地面點三維位移的幾何關系Fig.1 Geometric relationships between LOS displacement and 3D displacement of ground points

2 現有的礦區INSAR三維位移反演方法分析

由式(1)可知，當地表點的三維位移均為未知量時，由單個LOS 位移方程無法確定3個未知量。一些學者在利用InSAR位移數據提取礦區地表下沉值W時，直接將式(1)中的水平向位移分量UN，UE忽略不計，計算出下沉值W的近似公式:

(2)

由于地下開采引起的地表移動特征與常規的區域沉降不同，礦區移動盆地內的地表點均發生顯著的水平位移量，一般相當于下沉量的0.3～0.5倍，尤其在移動盆地邊緣地帶，其水平位移量甚至會大于下沉量。這種將LOS位移簡單地換算成垂直下沉的方法，不僅無法獲得對開采沉陷有重要影響的水平向位移，而且所得到的下沉量及其分布可能與實際情況相差甚遠。

一些文獻利用3幅具有顯著不同幾何形狀成像的SAR影像，分別提取3個雷達視線方向的一維位移，通過3個不同的LOS位移觀測方程聯合解算未知的三維位移量，即

dLOS1=Wcosθ1-sinθ1[UNcos(αh1-3π/2)+

UEsin(αh1-3π/2)]

dLOS2=Wcosθ2-sinθ2[UNcos(αh2-3π/2)+

UEsin(αh2-3π/2)]

dLOS3=Wcosθ3-sinθ3[UNcos(αh3-3π/2)+

UEsin(αh3-3π/2)](3)

式中,θi為第i組同步衛星入射方位角;αhi為第i組衛星飛行方向與北方向的夾角，i=1,2,3，分別是地面點在垂直、南北、東西3個方向上的位移。由于該方法需要至少3個不同幾何形狀成像的干涉像對，利用不同軌道參數處理得到的三組LOS位移值之間相互獨立，它們本身存在較大的系統性誤差。因此，利用式(3)解算出的地表點三維位移量不可避免地存在顯著誤差。

近年來，有學者依據現有的開采沉陷概率積分法預計模型，利用地表移動穩定后的下沉與水平位移之間存在的概率積分函數關系，將地表點水平位移變量表達為下沉變量的函數式，實現一維LOS向位移的反演[9]。根據概率積分法模型，開采引起的地表水平位移與垂直下沉梯度(即地表傾斜值)成正比，其數學關系如式(4)所示

UN=BΔWN

UE=BΔWE(4)

式中，ΔWN，ΔWE分別為下沉W沿東西、南北向的下沉梯度(即傾斜值);B為概率積分法中描述水平移動與下沉關系的參數，由地質采礦條件或經驗值確定[2]。將約束條件方程(4)與式(1)聯合求解可得三維位移值。上述模型曾在國際知名期刊Journal of Geodesy上發表[2]，為礦區InSAR三維形變場的求解提供了一種新方法。

上述模型假定開采沉陷中地表任意點的下沉變化和水平向位移分量之間均存在相同的(函數)比例關系，因此地表三維位移量中只有下沉量為未知變量，水平位移則視為下沉的關聯變量[10]。但大量的地表移動實測資料表明，地表移動盆地中下沉與水平位移之間的關系很復雜，且難以用概率積分法函數模型來描述。因此，這種反演方法基于下沉和水平位移之間存在確定的函數關系，因而存在較大的局限性。同時，在開采過程中的地表動態變形中，垂直位移和水平位移的動態分布特征十分復雜，現有的開采沉陷理論尚未建立動態垂直位移和水平位移的函數模型。

3 基于近水平煤層開采沉陷對稱的三維位移反演模型

3.1 近水平煤層開采沉陷的對稱分布特征

根據現有的開采沉陷理論，在常規的單一長壁工作面近水平煤層開采條件下，由于煤層傾角接近為零，上覆巖層和地表的幾何特征與物理性質均表現為各向同性，開采引起的覆巖與地表沉陷盆地則呈現對稱分布特征，其對稱軸為過采空區中央且平行(或垂直)于工作面推進方向的走向主斷面和傾向主斷面，地表移動盆地內任意點的下沉和水平位移均關于走向主斷面和傾向主斷面呈對稱分布[11-13]。因此，在煤礦布設地表移動觀測站時，常常利用上述開采沉陷的對稱性，僅在工作面開切眼或終采線一側設置半條觀測線，利用所得到的半個移動盆地主斷面上的觀測值，來描述整個地表移動盆地的特征[14]。

根據水平煤層開采沉陷的基本規律，在不同開采充分程度下地表主斷面上的下沉與水平位移分布如圖2所示。圖2表明，無論在非充分開采和超充分開采條件下，地表下沉和水平位移均以盆地中心o(采空區中心正上方最大下沉點)為對稱分布[5]。表移動盆地內任意點的水平位移矢量分布如圖3所示，所有點的水平位移方向均指向盆地中心o。其中，走向主斷面Z1Z2和傾向主斷面Q1Q2相交于o點，將地表移動盆地劃分為I，II，III，IV四個區域,對稱區I和III內任意兩個對稱點P1和P3的下沉量相等，水平位移則大小相等而方向相反，對稱區II和IV內的兩個對稱點P2和P4也是如此。根據這一特征可以建立地表任意2個對稱點上的LOS向位移與三維形變之間的函數模型，并進行疊加和差分處理，從而解算出三維形變分量。

圖2 近水平煤層開采地表移動主斷面上的下沉與水平位移特征曲線Fig.2 Surface movement curve of main section of near horizontal seam mining

3.2 地表移動穩定后的InSAR三維位移反演模型

如圖3所示，對于近水平煤層開采地表移動盆地的對稱區I和III中任意一組對稱點P1，P3，兩者的LOS向位移與地表水平位移的幾何關系如圖4所示。P1和P3水平位移方向與北方向的夾角分別為ω1，ω3，兩者相差180°。

圖3 地表移動盆地對稱區域的下沉和水平位移矢量圖Fig.3 Vector map of subsidence and horizontal displa-cement in the basin symmetry area ofsurface movement

圖4 采空區地表移動平面示意Fig.4 Sketch map of the surface movement of the goaf

圖4中，地表移動對稱點P1與P3的LOS位移dLOS1，dLOS3與垂直向、南北向、東西向位移的關系由式(1)可得:

dLOS1=cosθW1-sinθcos(αh-3π/2)UN1-

sinθsin(αh-3π/2)UE1

dLOS3=cosθW3-sinθcos(αh-3π/2)UN3-

sinθsin(αh-3π/2)UE3(5)

式(5)中，各符號的意義與式(1)相同。對稱點P1和P3的水平位移的大小相等，設為SP，在南北、東西方向的投影UN1，UN3，UE1，UE3由式(6)確定:

UN1=Spcosω1，UE1=Spsinω1

UN3=Spcosω3，UE3=Spsinω3(6)

式中，ω1，ω3分別為P1，P3與移動盆地中心o點之連線與北方向的夾角，由于ω1與ω3相差180°，所以無論在南北方向還是東西方向上，對稱點P1和P3的水平位移分量都是大小相等而方向相反，且其下沉量相等，則由式(6)可得:

UN3=-UN1

UE3=-UE1

W3=W1(7)

將式(6)，(7)代入式(5)可得:

dLOS1=cosθW1-sinθcos(αh-3π/2)Spcosω1-

sinθsin(αh-3π/2)Spsinω1(8)

dLOS3=cosθW1+sinθcos(αh-3π/2)Spcosω1+

sinθsin(αh-3π/2)Spsinω1(9)

將式(8)與式(9)相加后可得到對稱點P1，P3的下沉值為

(10)

將式(8)與式(9)相減后可得到對稱點P1，P3的水平位移值:

Sp=(dLOS1-dLOS3)/{-2sinθ[(cos(αh-

3π/2)cosω1+sin(αh-3π/2)sinω1]}(11)

將式(11)代入式(6)可以得到對稱點P1，P3在南北，東西方向的水平位移分量:

UN1=(dLOS1-dLOS3)cosω1{-2sinθ[cos(αh-

3π/2)cosω1+sin(αh-3π/2)sinω1]}=-UN3(12)

UE1=(dLOS1-dLOS3)sinω1{-2sinθ[cos(αh-

3π/2)cosω1+sin(αh-3π/2)sinω1]}=-UN3(13)

依據單個InSAR干涉對獲取得LOS位移值，由式(10)～(13)可解算出礦區地表移動盆地內任意點的三維位移量。

3.3 開采沉陷InSAR-LOS位移疊加分析

地表開采沉陷是一個持續時間達1～2 a的動態過程，地表穩定后的三維位移值是多期監測的沉陷疊加的結果。在InSAR礦區形變監測中，目前常采用差分干涉處理方法獲得不同時段的LOS位移，然后進行疊加得到最終的地表形變量。進行InSAR處理時要求兩幅SAR影像之間的時間基線不能太長，以免引起影像的失相干。目前現有的SAR衛星影像的重復周期一般為15～45 d[13]。因此，礦區地表移動盆地內的三維位移是通過多景SAR進行數據處理后并疊加得到的，地表任意點的LOS位移，相當于將每個時間段內的LOS位移疊加進行矢量求和。如果將同一SAR衛星在不同時段的衛星入射方位角和飛行方向視為常量，則LOS位移的疊加計算公式如下:

∑LOS=cosθ∑wi-sinθcos(αh-3π/2)

∑uNi-sinθsin(αh-3π/2)∑uEi(14)

式中,i=1～n(n為總時段數)Wi,UNi,UEi為第i時段內地表點垂直下沉、水平南北向、水平東西向的三維位移增量。式(14)中右邊第1項為地表點各時段下沉增量Wi之和在LOS向的投影值。由于下沉方向始終為垂直向下，則該項與式(1)中右邊第1項Wcosθ完全相同。式(14)的后2項為地表點各時段的南北向、東西向水平位移矢量之和在LOS向的投影值。

在工作面推進過程中，地表移動盆地中央任意點P的水平位移經歷著復雜變化，其位移軌跡如圖5所示。當工作面由遠處推進到位置1時，開采影響波及到地表點P，該點發生下沉伴隨著指向工作面推進邊界的水平位移，其方向與工作面推進方向相反;當工作推進過P點下方后，該點下沉速度達到最大值，而水平位移方向則發生逆轉，與工作面推進方向相同[5]。當工作面遠離P點后，回采工作面的影響逐漸減小，直至地表移動穩定。在常規地質采礦條件下，盆地中央地表點P在經歷下沉和動態水平位移后，基本穩定在其初始位置的正下方，總的水平位移值趨于零。對于移動盆地的其它點，無論處于開切眼一側還是終采線一側，移動過程中始終保持方向不變[5],其大小由零增加至穩定后的水平位移分量UN,UE，式(14)實質上與式(1)相同。因此，利用多時段InSAR處理的LOS位移相疊加后，所得到的地表穩定后的移動變形仍可以按照式(10)～(13)的反演模型進行三維位移分解。

圖5 盆地中中央地表點移動軌跡Fig.5 Central surface movement trajectory in the basin

3.4 地表移動過程中的InSAR三維位移反演模型

3.4.1 反演函數模型

在工作面推進過程中，地表動態沉陷盆地中開切眼一側和推進邊界一側對應點的下沉和水平位移不再具備對稱特征，上述關于穩定后的地表三維位移反演公式(10)～(13)不適用于動態三維位移反演。但是，由于近水平煤層的傾角接近為零，開采沉陷動態過程中的下沉和水平位移始終是關于走向主斷面對稱的。如圖6所示，工作面推進過程中地表點的移動變形始終以走向主斷面Z1Z2為對稱軸線。對稱點P1,P2的下沉值、沿走向方向的水平位移值均相等，沿傾向的水平位移值大小相等而方向相反。圖6中P1點水平位移方向指向動態的地表移動盆地中心o′，該方向與北方向的夾角為ω1。由于地表沉陷具有一定的滯后性，上述動態的移動盆地中心o′并不位于動態的采空區中央位置，一般會偏向開切眼一側。設工作面推進方向與北方向的夾角為α，則對稱點P1,P2的水平位移方向與Z1Z2的夾角均為(ω1-α)。

設走向、傾向的水平位移分別為UT,UI，沿南北、東西向的水平位移分量UN,UE，則兩位移分量之間的轉換關系為:

UT=UEsinα+UNcosα

UI=UEcosα-UNsinα

UN=UIsinα+UTcosα

UE=UIcosα-UTsinα(15)

圖6 開采過程中地表水平移動示意Fig.6 Horizontal movement of the surface in the mining

式(15)中α為工作面推進方向與北方向的夾角，即工作面推進方向的方位角。設圖6中對稱點P1與P2的LOS位移為dLOS1,dLOS2，參照式(1)可得LOS位移與垂直向、南北向、東西向位移W1,W2,UN1,UN2,UE1,UE2的關系式:

dLOS1=cosθW1-sinθcos(αh-3π/2)

UN1-sinθsin(αh-3π/2)UE1

dLOS2=cosθW2-sinθcos(αh-3π/2)

UN2-sinθsin(αh-3π/2)UE2(16)

根據近水平煤層開采地表動態移動的對稱特征，可得如下關系式:

UT1=UT2

UI1=-UI2

W1=W2(17)

將式(15),(17)代入式(16)并整理得:

dLOS1=cosθW1+(Bsinα+Ccosα)UI1+

(Bcosα-Csinα)UT1

dLOS2=cosθW1-(Bsinα+Ccosα)UI1+

(Bcosα-Csinα)UT1(18)

式中B=-sinθcos(αh-3/2π),C=-sinθsin(αh-3/2π)。將式(18)中的兩式進行求差得到傾向水平位移分量UI1,UI2為

(19)

由圖6和式(19)可得走向水平位移分量UT1,UT2為

式中,水平位移的方位角ω1為P1點與動態地表移動盆地中心o′之連線的方位角，將式(18)中的兩個等式相加得到:

dLOS1+dLOS2=2cosθW1+2(Bcosα-Csinα)UT1(21)

將式(20)代入式(21)即的下沉分量W1,W2為

UN1=(dLOS1-dLOS2)[sinαtan(α-ω1)+

cosα]/2tan(α-ω1)(Bsinα+Ccosα)

UN2=(dLOS1-dLOS2)[-sinαtan(α-ω1)+

cosα]/2tan(α-ω1)(Bsinα+Ccosα)

UE1=(dLOS1-dLOS2)[cosαtan(α-ω1)-

sinα]/2tan(α-ω1)(Bsinα+Ccosα)

UE2=(dLOS1-dLOS2)[-cosαtan(α-ω1)-

sinα]/2tan(α-ω1)(Bsinα+Ccosα)(22)

將式(19),(20)代入式(15)得

UN1=(dLOS1-dLOS2)[sinαtan(α-ω1)+

cosα]/2tan(α-ω1)(Bsinα+Ccosα)

UN2=(dLOS1-dLOS2)[-sinαtan(α-ω1)+

cosα]/2tan(α-ω1)(Bsinα+Ccosα)

UE1=(dLOS1-dLOS2)[cosαtan(α-ω1)-

sinα]/2tan(α-ω1)(Bsinα+Ccosα)

UE2=(dLOS1-dLOS2)[-cosαtan(α-ω1)-

sinα]/2tan(α-ω1)(Bsinα+Ccosα)(23)

式(22),(23)即為工作面推進過程中地表三維位移的分解模型。式中的中間變量B,C與式(18)相同。

3.4.2 動態地表移動盆地中心o′的確定

為了求取地表移動過程中的水平位移方向角ω1，必須確定動態的地表移動盆地中心o′。由于地下開采影響向上傳播至地表是一個時間過程，使得動態的地表移動盆地中心o′滯后于對應的采空區中央位置，但滯后的具體位置因地質采礦條件不同而難以確定。根據開采沉陷的基本規律，近水平煤層開采地表移動過程中的動態移動盆地中心o′點在平行于推進方向并通過采空區中間的走向主斷面上，且具備以下特征:

(1) 在整個走向主斷上動態盆地中心o′點的下沉值應為最大值;

(2) 動態盆地中心o′點處的水平位移絕對值應趨近于零。

依據上述特征，基于走向主斷面線Z1Z2上InSAR-LOS位移的疊加計算，來搜索動態移動盆地中心o′。在式(18)雷達視線(LOS)方向與下沉、傾向水平位移、走向水平位移的幾何關系式中，由于垂直于走向主斷面上的橫向(即傾向)水平位移為零，則由式(18)可得走向主斷上任一點的三維位移分解公式:

dLOS=cosθW+(Bcosα-Csinα)UT(24)

由于盆地中心o′兩側對稱點的下沉值相等，水平位移則大小相等而方向相反(均指向中心o′)，則將中心點o′兩側對稱點(可以選點o′兩側的相鄰像元)的LOS位移進行疊加并取均值可得下沉值Wm，進行求差后可得水平位移值Um。由于該點處于動態盆地中心o′，則下沉值Wm應為整個走向主斷面上的最大下沉值，水平位移值Um應趨近為零。因此，在整個走向主斷面上針對任意兩個相鄰像元或一定窗口區域的LOS位移值或位移均值，進行上述疊加和求差計算，搜索出滿足下沉值達到最大且水平位移趨于零的像元對，取其中間位置即為動態盆地中心o′。具體步驟如下:

(1)將LOS位移圖中走向主斷面上的像元位移進行編號(LOS1,LOS2,LOS3,LOS4,…,LOSn-2，LOSn-1，LOSn)，由于盆地中心不可能偏到采空區以外，搜索起始點可以從開切眼正上方的地表位置開始。

(2)選取1,2號像元的位移值，用dLOSA,dLOSB表示，分別代入式(24)得:

dLOSA=cosθWA+(Bcosα-Csinα)UTA

dLOSB=cosθWB+(Bcosα-Csinα)UTB

WA=WB

UTA=-UTB(25)

將式(25)中的dLOSA與dLOSB分別進行疊加和差分，計算出像元1,2的下沉與水平位移值，記錄為第1組移動值(W1,UT1)。

(3)選取2,3號像元的位移值，分別賦值為dLOSA,dLOSB，由式(25)計算出像元2,3的下沉與水平位移值，記錄為第2組數據(W2,UT2)。

(4)以此類推，直至計算出第n-1,n號像元的下沉和水平位移值，記錄為第n-1組數據,形成位移值數組(W1,UT1),(W2,UT2),…,(Wn-1,UTn-1)。計算流程如圖7所示。

圖7 移動盆地中心搜索計算流程Fig.7 Flow diagram of mobile basin center search

(5)將上述全部(n-1)數組進行比較，分別搜索出下沉值最大的像元組合，以及水平位移絕對值最小的像元組合。當兩者位置一致時，該像元組的中間位置即為動態移動盆地的中心o′。當2個條件發生沖突時，優先將下沉值最大的位置選為盆地中心。

應該指出，當地表沉陷達到超充分采動后，動態移動盆地中央會出現“平底”現象，按照上述搜索方法將得到多組滿足條件的像元區域。此種情況下開切眼一側的半盆地內水平位移指向第一組滿足條件的像元位置，而推進邊界一側半盆地內水平位移指向最后一組滿足條件的像元位置。

4 三維位移反演算法實現與驗證

4.1 地表穩定后的三維位移反演計算流程

開采沉陷三維反演模型式(10)～(13)中涉及的參數包括:

(1)衛星入射方位角θ，可從雷達遙感數據的頭文件中獲取;

(2)衛星飛行方向與北方向的夾角αh，順時針方向為正。該參數可在雷達遙感數據的頭文件中獲取。

(3)地表點P1的水平位移方向與北方向的夾角ω1。在InSAR形變圖所處的坐標系統中，獲取任意地表P1及采空區中心O的平面坐標P1(N1，E1),O(NO，EO)，通過坐標反算得到兩點之間的距離Lp及兩點連線與北方向的夾角(即方位角)ω1。先用公式(15)求出P1O的象限角λ,再根據相應的象限條件將λ轉換為方位角ω1。

λ=arctan[(N0-N1)/(EO-E1)](26)

(4)地表點P1在移動盆地內的對稱點P3的坐標(N3，E3)根據對稱幾何關系由下式計算:

N3=NO+LPcosω1

E3=EO+LPsinω1(27)

地表任意點P1及其對稱點P3的LOS向位移dLOS1,dLOS3，可根據其坐標為索引，利用InSAR形變數據文件或LOS形變圖進行內插得到。因此，利用數學模型式(10)～(13)反演地表點的三維位移的算法流程如圖8所示。

圖8 穩定后的三維位移反演算法流程Fig.8 Algorithm flow of 3-D displacement inversion of Insar

根據上述算法在MATLAB中進行編程實現，將L向位移數據分解為下沉、水平南北向、東西向的三維位移數據，并繪制出開采沉陷三維位移分布圖。

4.2 動態地表三維位移反演計算流程

動態的地表三維位移反演模型中的衛星入射方位角θ、衛星飛行方向與北方向的夾角αh、地表點的水平位移方位角ω1、工作面推進方向的方位角α等參數含義與穩態的模型相同。此外，還要獲取對稱點P1，P2的坐標(N1,E1),(N2,E2)，并按照上述搜索算法確定動態的地表移動盆地中心O′的坐標(NO′，EO′)。其反演計算流程與圖8基本相同，不再贅述。應注意的是，對稱點P2坐標是利用與P1關于走向主斷面Z1Z2對稱的條件來求取的，與穩態情況下對稱點P3坐標的計算公式有所不同。

4.3 反演模型的驗證

4.3.1 地表移動穩定后的三維位移反演模型驗證

某礦開采近水平煤層，地表較為平緩。采用長壁工作面開采方式，工作面由正南向北推進。開采寬度L1=250 m，走向長度L3=500 m，開采厚度m=2 m，地表為平地，平均采深H=500 m;煤層傾角α=0°。工作面回采時間為2014-10-16—2015-04-18，平均每月推進83 m。在工作面正中央地表布設了走向和傾向觀測線。根據2015-12-03號最后一次實測數據繪制出地表穩定后的下沉與水平位移曲線，如圖9所示。

在2014-09-26—2015-12-10的時間段內，利用15景sentinel SAR影像做差分干涉處理并將LOS向位移進行疊加，繪出地表穩定后的LOS向位移等值線，經過平滑簡化后如圖10所示。LOS位移的最大值明顯偏離采空區中心位置。

圖9 實測數據繪制的下沉與水平位移曲線Fig.9 Subsidence and horizontal displacement curve

圖10 LOS向位移Fig.10 LOS displacement

圖11 忽略水平位移的下沉Fig.11 Subsidence diagram that neglected the horizontal displacement

首先采用常規忽略水平位移的方法，直接將圖10中的LOS位移投影為垂直下沉值。SAR影像的衛星入射角θ=30°，衛星飛行方向與北方向的夾角αh=345°，按式(2)計算下沉量，繪出下沉等值線分布如圖11所示。所得到的最大下沉值明顯偏離采空區中央上方，與圖9的實測數據相差較大。

按照圖8的計算流程對圖10的LOS位移進行三維分解。應用專門的MATLAB 程序計算地表格網點的下沉、南北和東西水平位移值，并生成對應的等值線圖，如圖12所示。圖12(a)中，最大下沉值位于采空區中心附近，地表下沉、南北(走向)水平位移、東西(傾向)向水平位移基本呈現以采空區中心(走向和傾向主斷面線)為對稱的分布特征。

圖12 LOS向位移三維分解Fig.12 LOS displacement to 3-D decomposition

4.3.2 動態地表三維位移反演模型驗證

選取2014-09-26—2015-02-20的前5 h段LOS向疊加位移，進行動態的三維位移分解。此時間段工作面推進距離為340 m，地表沉陷處在動態過程中。繪出疊加至2015-02-20的地表LOS向位移等值線圖，經過簡化和平滑處理后如圖13所示。地表LOS向位移明顯偏離動態的采空區中央。按照上述搜索計算方法確定地表動態的移動盆地中心位置，偏離采空區中心43 m。

按照動態的三維位移分解模型式(22),(23)對圖13的LOS向位移進行三維分解，生成地表下沉、南北向水平位移、東西向水平位移的等值線圖，如圖14所示。可以看出，工作面推進到340 m時，地表動態下沉、南北(即走向)、東西(即傾向)位移在走向主段面線的兩側都基本呈現對稱特征，而在開切眼和推進邊界上方地表的下沉和水平位移分布均為不對稱狀態。應該指出，本例中工作面走向為正南北向，傾向為東西向。如果工作面推進方向與北方向斜交，則應將計算出的南北向、東西向水平位移按照式(15)轉換為沿工作面走向和傾向的水平位移后，再繪制走向和傾向水平位移等值線圖。

圖13 動態的地表LOS向位移Fig.13 Dynamic surface LOS displacement

圖14 動態的地表三維位移Fig.14 Dynamic 3-D surface displacement

4.3.3 反演模型的精度驗證

以傾向線實測數據來驗證本文模型反演得到的下沉、水平位移值，對比曲線如圖15所示。統計表明，靜態模型反演的下沉值(圖15(a)曲線2)與實測值(圖15(a)曲線1)相比的均方根誤差為±20.1 mm，傾向水平位移(圖15(b)曲線2)與實測值(圖15(b)曲線1)相比的均方根誤差為±22.8 mm。

圖15 傾向線上反演曲線與實測曲線對比Fig.15 Comparison of displacement inversion curve on dip line with measured curve on dip line1—穩態實測曲線;2—穩態模型反演的曲線;3—推進至340 m時的實測曲線;4—動態模型反演曲線;5—直接忽略水平位移反演的穩態下沉曲線

如果忽略水平位移直接用LOS向位移來計算地表下沉值(圖15(a)曲線5)，則與實測下沉值相比的均方根誤差達±104 mm。工作推進至340 m時，基于動態模型反演的下沉值(圖15(a)曲線4)與實測值(圖15(a)曲線3)相比的均方根誤差為±17.5 mm，傾向水平位移(圖15(b)曲線4)與實測值(圖15(b)曲線3)相比的均方根誤差為±19.2 mm。考慮到InSAR監測本身的誤差影響，上述穩態和動態的模型反演三維位移的精度本身是可靠的。

5 結 論

(1)本文根據近水平煤層開采地表移動盆地中的下沉和水平位移所具備的對稱特征，建立了一種基于InSAR監測數據反演三維位移的新模型。現有的開采沉陷研究及大量實測資料已經證實，在近水平煤層長壁工作面開采及地形起伏較小的情況下，地表移動盆地在形態上確實存在明顯的對稱特征。在工作面推進中，地表動態下沉和水平位移也同樣存在以走向主斷面為軸線的對稱特征。因此，本文所建立的穩態和動態三維位移反演模型遵循了近水平煤層開采沉陷的基本原理。模型驗證表明，本文模型反演的結果符合實際情況，明顯優于忽略水平位移而直接分解地表下沉的方法。

(2)由于礦區地表下沉和水平向位移本身屬于獨立的形變分量，其相互關系難以用確定的函數模型來描述。因此，相對于現有在概率積分法基礎上推導的InSAR三維位移反演模型來說，本文方法具有更好的普適性和實用性。

(3)在穩定后的地表移動三維分解中，可以利用地表移動關于傾向和走向都存在對稱的特征(圖3中I,II,III,IV區均為對稱區)，將地表點的三維位移作為3個獨立未知量，列出四個區域對稱點的4個LOS向位移方程式，按照最小二乘原理進行三維位移解算。所以，在地表移動穩定后的三維位移反演中，也可以將動態的位移方程式(16)并入穩態反演模型中，求出三維位移的最小二乘解。

應該指出，由于礦區開采沉陷規律受多種復雜的地質采礦因素影響，在傾斜煤層、山區地形及不規則工作面開采條件下，地表移動變形可能喪失對稱性特征，筆者正在進一步研究非對稱沉陷情況下的InSAR三維位移反演模型。