基于InSAR三維分解技術(shù)的貴州省貞豐某煤礦地下采空區(qū)反演*

姚佳明 姚 鑫 吳作啟 陳 劍 蔣永亨 劉星洪

（①活動構(gòu)造與地殼穩(wěn)定性評價重點實驗室，中國地質(zhì)科學(xué)院地質(zhì)力學(xué)研究所，北京100081，中國）

（②中國地質(zhì)大學(xué)（北京）工程技術(shù)學(xué)院，北京100083，中國）

（③煤炭科學(xué)技術(shù)研究院有限公司，北京100013，中國）

0 引 言

貴州貞豐縣煤礦采空區(qū)眾多，小煤窯的開采易造成地表塌陷、滑坡等一系列地質(zhì)災(zāi)害。而且在一些煤炭資源豐富的地區(qū)非法采礦事件時有發(fā)生，不僅嚴重破壞礦產(chǎn)資源和生態(tài)環(huán)境，易誘發(fā)斜坡變形開裂甚至失穩(wěn)（史文兵等，2016），造成礦山災(zāi)害和巨大的經(jīng)濟損失，威脅安全生產(chǎn)和社會穩(wěn)定。由于礦山分布廣泛，地下非法采礦活動隱蔽性強，僅依靠礦山執(zhí)法部門的調(diào)查，很難及時、準確地發(fā)現(xiàn)這些行為，需要耗費大量的時間、精力和財力（Xia et al.，2018）。了解地下開采情況可以對將來可能由采礦引起的地表破壞進行防護，減少生命財產(chǎn)損失，并且對煤礦越界開采起到一定監(jiān)測作用。我國煤礦采空區(qū)探測的技術(shù)方法較多，比如地球物理方法、鉆探方法、三維激光掃描等。地球物理方法利用波在不同地層和巖石的彈性、速度和運動方式差異來確定地下采空區(qū)（吳成平等，2007），然而物探類方法受環(huán)境噪音、適用性、分辨率、多解性、不確定性等多種條件的限制，無法做到準確定量預(yù)報（王升陽等，2015）；鉆探在工程施工、探測方面應(yīng)用較多，但施工難度大、成本高，探測成果受鉆孔數(shù)量、鉆孔密度、配套的監(jiān)測、測試方法等限制，分辨率較低，無法查清周邊采空區(qū)準確邊界（薛國強等，2018）；三維激光掃描技術(shù)可以高精度、高效率地測量出地下采空區(qū)，但性價比較低（薛國強等，2018）。對于已關(guān)閉或已廢棄煤礦，上述方法雖然可以在一定精度上測量出地下采空區(qū)范圍，但需要消耗大量人力和財力且有一定的作業(yè)危險性。

近年來很多學(xué)者將InSAR技術(shù)應(yīng)用于監(jiān)測與地球或人類活動有關(guān)的地表變形（Amelung et al.，1999；Du et al.，2017；Feng et al.，2015；Milillo et al.，2016a，2016b；Xu et al.，2016a，2016b；Zhao et al.，2016；Yang et al.，2017a，2017b；Blachowski et al.，2019），與傳統(tǒng)的地面水準測量技術(shù)、GPS技術(shù)相比，InSAR技術(shù)具有全天時、全天候、觀測范圍大、靈敏性高、高性價比、觀測時間可回溯等優(yōu)勢，并能穿透某些地物表面（吳立新等，2004）。煤礦開采產(chǎn)生的大區(qū)域變形可以使用InSAR技術(shù)獲取，在國內(nèi)外已有大量學(xué)者將InSAR技術(shù)應(yīng)用于礦區(qū)地面塌陷監(jiān)測，并通過實地GPS等地面測量驗證了InSAR計算結(jié)果的準確性。朱建軍等（2011）介紹了InSAR技術(shù)在礦區(qū)地表形變監(jiān)測中的應(yīng)用現(xiàn)狀及進展，分析了D-InSAR技術(shù)相比于傳統(tǒng)測量手段的優(yōu)勢；Xia et al.（2018）采用D-InSAR采礦沉陷監(jiān)測方法，建立了地表變形與以沉陷為特征的地下開采的時空關(guān)系模型；Fan et al.（2015）采用相疊加和亞像素偏移跟蹤相結(jié)合的方法，對采動引起的沉陷機理進行了研究；Wang et al.（2018）結(jié)合D-InSAR和概率積分法提出了一種基于單視距D-InSAR的概率積分參數(shù)模擬沉降的計算方法。

煤礦地下采空與地表破壞存在一定變形規(guī)律，其變形量、變形范圍與地下開采方式、煤層埋深、開采煤層厚度、地層傾角、巖層組合、地表地形、坡面與地層的組合關(guān)系以及開采工作面的布置方向有直接關(guān)系。通過實地調(diào)查發(fā)現(xiàn)，研究區(qū)內(nèi)地質(zhì)構(gòu)造簡單，地下煤層展布特征基本一致，順層緩傾巖層以砂泥巖為主，通過周邊煤礦調(diào)查得知該區(qū)域開采方式一致皆采用下行式開采，且地下開采方向較為一致（自西向東開采），故在以上條件限制下可以根據(jù)該地區(qū)開采沉陷參數(shù)進行地下采空范圍反演，所選參數(shù)參照于《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采規(guī)范》（以下簡稱《三下采煤規(guī)范》）中的邊界角概念，建立InSAR短期動態(tài)監(jiān)測煤礦的影響角計算方法，本文目的是通過該方法計算地下采空地表變形趨于穩(wěn)定之后，在開采時間與穩(wěn)定時間段內(nèi)的短期地表變形與地下采空變形間的影響角規(guī)律，而后通過計算獲取研究區(qū)已知煤礦地下開采與地表沉降變形關(guān)系，來推算反演已廢棄煤礦的地下開采范圍及時間。

1 研究區(qū)概況

研究區(qū)位于貴州省黔西南布依族苗族自治州貞豐縣境內(nèi)，坐標范圍為北緯25°20′～25°21′，東經(jīng)105°31′～105°33′，變形區(qū)的坡度10°左右，坡向近北偏西30°（圖1）。

研究區(qū)地層從石炭系至第四系均有出露，其中以三疊系地層最為發(fā)育。貞豐縣主要含煤地層為三疊系上統(tǒng)火把沖組（T3h），由海陸交互相和沖積相含煤沉積組合形成，也是本研究區(qū)煤礦主要開采煤層（圖2）。T3h地層由石英砂巖、粉砂巖、黏土巖、碳質(zhì)黏土巖及數(shù)十層煤層（線）組成多個沉積旋回，與下伏地層把南組整合接觸，礦區(qū)內(nèi)僅出露至火把沖組第三段（T3h3）：淺灰、灰、黃灰綠色等色厚層，塊狀粗至細粒石英砂巖、細砂巖、黏土巖組成的不等厚韻律層，砂巖中具有斜層理。該區(qū)含有可采煤層兩層，編號為JK、K2。JK煤層厚0.70～0.80 m；K2煤層平均厚1.08 m，煤層傾角平均為10°，傾向290°～350°，礦區(qū)內(nèi)構(gòu)造比較簡單，巖層為單斜構(gòu)造。礦區(qū)斷層較發(fā)育，有柏枝樹斷層、大榜斷層（圖1紅實線），但變形區(qū)與斷層無接觸關(guān)系且相距較遠，故排除受到斷層影響。

經(jīng)過實地考察，受地下煤層和當?shù)貤l件限制，研究區(qū)內(nèi)礦井設(shè)計采用下行式開采，煤層開采順序為K2、JK聯(lián)合開采，平均層間距8～10 m，煤層傾角6°～15°，平均傾角12°，采深約250 m，采厚1.8 m，工作面走向長度420 m，傾向長度100 m，開采面積約42 000 m2。

2 研究與計算方法

2.1 開采沉陷參數(shù)計算與地下采空區(qū)反演

圖1 研究區(qū)位置與地質(zhì)圖Fig.1 Location and geological map of the study area

圖2 研究區(qū)地質(zhì)剖面圖（位置見圖1剖面線I－I′）Fig.2 Geological profile of the study area

本文研究區(qū)變形為已知煤礦開采和未知煤礦開采引起的兩部分變形。根據(jù)對本文研究區(qū)附近煤礦的調(diào)查，由于開采煤層與地形因素的影響，該地區(qū)煤礦開采方式皆為下行式開采，地層產(chǎn)狀、開采煤層、厚度等和鄰近煤礦開采條件較為一致。可以根據(jù)已知煤礦區(qū)計算的沉陷參數(shù)與變形時間滯后關(guān)系，反演未知區(qū)煤礦地下開采時間和范圍。

根據(jù)各時段地表LOS形變圖判斷地下掘進方向與地表變形時間以及在變形時段地表的累計變形，進行三維分解計算來獲取地表垂向變形信息。研究從兩個方面進行（圖3）：（1）在獲取已知煤礦地下開采資料和地質(zhì)條件的情況下，根據(jù)已知煤礦信息與地表變形信息計算出該區(qū)域開采沉陷參數(shù)（上、下山及走向影響角）和變形滯后時間；（2）根據(jù)已知煤礦開采計算獲得的參數(shù)和未知廢棄煤礦開采引起的地表變形反演其地下采空范圍和地下開采時間。

2.2 計算方法

本文使用D-InSAR方法的計算結(jié)果，對煤礦采空區(qū)三維地表變形進行反演，結(jié)合數(shù)值模擬技術(shù)對反演結(jié)果進行核查驗證。本文選取了日本JAXA宇宙航空研究中心研發(fā)的L波段ALOS-PALSAR2（以下簡稱P2）衛(wèi)星條帶模式數(shù)據(jù)，升、降軌各8景SAR影像。數(shù)據(jù)的斜距向和方位向分辨率分別為1.430 422 m和2.217 833 m，影像參數(shù)如表1。對比其他數(shù)據(jù)，本文所選的數(shù)據(jù)波長較長、分辨率較高、穿透力強，更有能力識別出礦區(qū)地表變形而不出現(xiàn)失相干等現(xiàn)象，適用于礦區(qū)地表變形監(jiān)測與分析。

使用D-InSAR方法分別對升、降軌影像做了差分干涉處理，為了排除地表大梯度變化帶來的失相干，時間基線就僅選取36對相鄰的影像對，來減少時間去相關(guān)造成的影響。處理過程中為了保持較高分辨率和信噪比，我們選擇了2×2的多視比，采用偏移多項式對影像進行配準，使用美國SRTM的30 m分辨率DEM去除地形相位。對干涉圖像使用最小費用流法（mcf）進行解纏，居民區(qū)較為穩(wěn)定，故選擇煤礦東北側(cè)居民區(qū)為解纏起算點，干涉時間基線圖如圖4所示。影像覆蓋區(qū)域相對于河谷區(qū)高差較小，且研究區(qū)范圍內(nèi)河流較少，水汽含量均一，大氣誤差較小，誤差源主要為軌道誤差。故本文在處理D-InSAR過程中，在兩次移除軌道誤差后，高程帶來的大氣誤差的去除和適應(yīng)性濾波去除湍流大氣均采用大窗口濾波，減少對礦區(qū)實際變形量的影響，最后投影到UTM坐標。

煤礦區(qū)地表變形為水平向與垂直向的三維變形，而D-InSAR計算出的結(jié)果為視線向（LOS）變形，若忽略水平運動而直接將LOS變形值進行轉(zhuǎn)化，求得的垂直向變形會存在很大誤差。國內(nèi)外很多學(xué)者將InSAR三維分解技術(shù)應(yīng)用于形變監(jiān)測與地表變形（王晨興等，2011；Li et al.，2015；Yang et al.，2017a，2017b；Huang et al.，2017；李凌婧等，2017；Liu et al.，2018；Yang et al.，2018；Zheng et al.，2019）。Hu et al.（2014）詳細介紹了目前國內(nèi)外學(xué)者在InSAR三維分解技術(shù)方面的發(fā)展情況，分別分析了各種方法的計算原理及其優(yōu)勢與劣勢。通過利用從至少3個成像幾何形狀獲得的多個InSAR測量，我們理論上可以將位移矢量擴展到3D（吳立新等，2004）。本文使用的方法是基于升、降兩軌計算方法，利用兩軌D-SAR影像干涉結(jié)果，根據(jù)升、降軌雷達成像幾何原理進行三維變形求解，因為只有升、降兩個方向的SAR影像，對于模糊度較高的N－S方向變形量不易準確計算，所以本文沒有N方向變形量分解，計算公式如下：

式中：du、de為垂直向、東西向變形量；Г為轉(zhuǎn)換系數(shù)；dlos，1、dlos，2為升、降軌兩干涉對變形值；?inc，i、αaz，i分別為入射角和方位角。

3 結(jié)果與分析

圖4 升、降軌基線圖（升－左、降－右）Fig.4 Baseline diagram of ascending and descending（ascend-left，descend-right）

通過兩軌干涉變形信息判斷，研究區(qū)變形可以分為兩部分，由兩個煤礦開采誘發(fā)：一部分為南側(cè)已知煤礦（已獲取地下開采資料）開采誘發(fā)的地表變形；另一部分為北側(cè)（已廢棄，無地下開采信息）未知煤礦開采帶來的地表塌陷變形（圖5）。通過南側(cè)煤礦開采信息（地下開采時間、開采范圍）和地表形變信息來計算開采沉陷參數(shù)，再根據(jù)未知煤礦地表變形信息反演地下開采時間和范圍。

3.1 已知區(qū)域開采沉陷參數(shù)和變形滯后時間計算

圖6為升、降軌在已知煤礦開采時間地表各時段地表變形信息，地下煤礦為自西向東開采，開采時間在2018年1月開始，2018年6月結(jié)束。兩個軌道的干涉影像時間不同，有助于形變區(qū)變形時間確定，在升、降軌共6個變形時間段內(nèi)，地表變形區(qū)隨時間的推移由西移動至東側(cè)，與地下開采方向一致，地表變形的時間與地下開采的時間存在一定延遲，滯后時間規(guī)律影響因素較多，在不同地區(qū)無法進行判斷與分析，但在地質(zhì)條件、煤層條件、開采方式與方向等相同的鄰近區(qū)域，滯后時間較為相近。

3.1.1 D-InSAR三維分解

根據(jù)煤礦開采及前期鉆探資料，已知煤礦區(qū)開采時間段為2018年1月至6月，為覆蓋變形時間段，本文的升、降軌數(shù)據(jù)分別選用了20171112～20180819和20171206～20180620共四景影像用于兩兩干涉處理。通過兩個干涉對進行三維分解解算，計算出地面垂向變形量。如圖7所示，在煤層開采的走向和傾向方向分別對垂向變形做變形剖面圖，開采面中心300 m長度范圍內(nèi)沉降量在20 cm以上，最大沉降30 cm，周邊區(qū)域變形量遞減。

圖6 煤礦開采過程多時間節(jié)點地表變形Fig.6 Known surface deformation diagram in each period of coal mining time

圖7 三維分解的垂向變形（a為已知煤礦地表變形垂向變形）Fig.7 Vertical deformation of three-dimensional decomposition（a shows the vertical deformation of known coal mine surface deformation）

3.1.2 開采沉陷參數(shù)計算

參考《三下采煤規(guī)范》中的邊界角概念，建立以地表InSAR監(jiān)測10 mm變形邊界與地下開采邊界連線與水平線的角度為InSAR監(jiān)測煤礦開采變形影響角，本文根據(jù)三維變形中的垂向變形量，結(jié)合地下煤礦開采邊界對已知煤礦開采區(qū)的走向，上、下山影響角做了計算（圖8），計算得出上山影響角75°，下山影響角80°，走向影響角83°。

3.1.3 變形時間滯后關(guān)系

地下采空反映到地表形成塌陷需要一定的時間，開采條件、地質(zhì)條件、開采方向等因素皆會影響變形滯后時間，而在本研究區(qū)上述條件較為類似，根據(jù)SAR影像干涉對時間節(jié)點的地表變形和地下開采時間，在本研究區(qū)做如圖9地表變形時間與地下采空時間的關(guān)系曲線，兩條曲線在橫軸（時間軸）方向的差值或偏移即為地表變形滯后時間，根據(jù)曲線走向趨勢可以判斷出在研究區(qū)的地質(zhì)和開采條件下，地下采空至地表變形處1個月的時間滯后關(guān)系。

3.2 未知區(qū)地下采空區(qū)反演

圖10為升、降軌不同時間段的干涉影像，地表變形移動與煤礦采空引起的變形一致，兩軌影像在近乎相同時間內(nèi)的地表形變量、形變區(qū)大小基本一致也驗證D-InSAR結(jié)果的準確性。通過兩軌干涉圖地表變形信息分析得出，該地區(qū)煤礦自西向東開采，變形時間在2017年4月26日至2017年10月25日之間。隨時間的推移，地表變形區(qū)也隨之向東移動，該地區(qū)屬干旱地區(qū)，年均降雨量少、植被稀疏，地下水位變化不大；在歷史SAR干涉影像中未開采時段該地區(qū)無變形產(chǎn)生，且該地區(qū)地表坡度、地下巖層較緩，不易發(fā)生滑坡，故也可以排除因地下水位變化、滑坡等引起的地表變形。

從圖10中8景干涉對中得出變形時間大致可以確定為：2017年4月末至2017年10月末。在考慮影像獲取時間與誤差允許范圍內(nèi)，我們對該區(qū)升軌20170514～20171112和降軌20170426～20171025兩個相近時間段做三維分解解算。

3.2.1 D-InSAR三維分解

圖8 走向、上山、下山影響角計算示意圖（橫縱軸1:1）Fig.8 Schematic diagram for calculating the ascending and descending influence angles（horizontal and vertical axis 1:1）

圖9 地下采空與地上沉陷時域?qū)Ρ葓DFig.9 Time-domain comparison of underground mining and overland subsidence

圖10 a～d為升軌干涉圖，e～h為降軌干涉圖Fig.10 a～d shows the ascending orbital interference pattern，e～h shows the descending orbital interference pattern

獲取的升軌影像入射角 ?inc，i、方位角 αaz，i分別為39°、－10°，降軌影像入射角 ?inc，i、方位角 αaz，i分別為43°、－170°，經(jīng)過2×2重采樣并地理編碼后的P2 SAR影像距離向和方位向分辨率分別為4.5 m和4.4 m，故三維分解的垂直向和東西向分辨率也為4.5 m和4.4 m。對升、降軌兩個干涉對進行等尺寸裁剪獲得1290×1290的干涉對，根據(jù)如上參數(shù)和公式計算得出了垂直和東西向變形量。

在圖11a、圖11b以及5條剖面線結(jié)果中看出，圖11a為垂直向形變結(jié)果，整體變形較為集中，最大形變量在區(qū)域中心為0.5 m，周圍大部分區(qū)域變形為0.2～0.4 m左右；圖11b為東西向形變結(jié)果，紅色值和藍色值分別代表運動向東和向西，從EE’剖面線中可見明顯分界線。中心左側(cè)整體為向東運動，符合采礦引起的地表塌陷變形特征。最大位移值0.5 m，中心右側(cè)整體向西運動，最大位移值0.36 m。

3.2.2 采空區(qū)邊界反演及開采時間推算

本文在廢棄煤礦采空沉降中心區(qū)域建立5條剖面線來反演地下開采范圍。在開采傾向（近自西向東）方向的4條剖面線最大沉降、10 mm沉降邊界位置以及沉降趨勢大致一致，故取4條剖面線平均值來反演開采傾向邊界；主斷面上地表移動盆地的范圍最大，移動最充分、移動形變量大，故在開采走向方向最大區(qū)域變形處拉取一條剖面線作為走向邊界的反演，根據(jù)已知影響角、地下開采深度、煤層產(chǎn)狀以及地表高程和形變信息通過三角函數(shù)幾何計算出地下開采邊界（圖12，圖13）。

綜合走向、傾向確定的開采范圍，可以確定開采面位置如圖14。

根據(jù)地表變形時間和滯后關(guān)系，本文對地下開采時間進行了推測（圖15）。在鄰近已知礦區(qū)的地下開采與地表破壞30 d左右的時間滯后關(guān)系判斷，地表破壞對應(yīng)的地下采空時間應(yīng)為地表變形時間的1個月之前，在兩軌D-InSAR各時段節(jié)點變形結(jié)果中做出地表變形在煤層開采走向上的時間－距離曲線（藍、黃虛線），兩曲線變形時間－距離曲線基本吻合，再根據(jù)滯后時間推算出地下煤礦的開采時間曲線（紅色虛線）。

3.2.3 采掘范圍數(shù)值模擬驗證

根據(jù)該區(qū)地質(zhì)巖層剖面與煤層信息，在MIDAS中建立簡化FLAC3D計算模型（圖16），計算模型范圍包含了未知采空區(qū)，尺寸為680 m×400 m×300 m，模型上表面為該地區(qū)30 m分辨率DEM雙線性重采樣至5 m分辨率構(gòu)成，依據(jù)當?shù)孛旱V鉆孔勘察報告該地區(qū)煤層上覆地層主要為泥質(zhì)砂巖與泥質(zhì)粉砂巖，為方便計算地層劃分過程中忽略了較薄的泥巖、粉砂等地層；在煤層位置設(shè)置東西向長380 m、南北向長150 m、傾向330°、傾角10°的開采煤層，具體分層見圖16。

圖11 基于D-InSAR的升、降軌三維分解Fig.11 D-InSAR based three-dimensional decomposition of ascending and descending orbits，where Fig.11a represents vertical deformation and red integrity represents settlement；Fig.11b is the east-west deformation，the red negative value represents the east，and the blue represents the west；A－A′equilinear lines are shaped variable profile lines；the NaN white areas are incoherent and without data areas

圖12 開采煤層傾向方向的地下開采區(qū)反演Fig.12 Inversion of underground mining area in mining seam inclination direction

圖13 開采煤層走向方向的地下開采區(qū)反演Fig.13 Inversion of underground mining area in strike direction of coal seam

計算模型左、右邊界x＝0im和x＝680 m兩個垂直面僅約束邊界面法向位移，平面內(nèi)無約束；前、后邊界y＝0 m和y＝400 m兩個垂直面僅約束邊界面法向位移，平面內(nèi)無約束；模型底部z＝－300 m水平邊界采用固定約束；水平地表面為自由面，收斂標準為體積最大不平衡力與典型內(nèi)里的比率小于10－5。巖體采用六面體八節(jié)點單元模擬，本構(gòu)模型采用基于彈塑性理論分析的莫爾－庫侖模型。研究區(qū)降雨較少，開采面上無地下水賦存，未考慮孔隙水、靜水壓力，巖體物理力學(xué)參數(shù)選取見表2。

圖14 開采區(qū)位置圖Fig.14 Location of mining area

圖15 地下開采時間推測Fig.15 Estimation of underground mining time

表2 巖體物理力學(xué)參數(shù)表Table 2 Physical and mechanical parameters of rock mass

計算過程如下：

（1）依據(jù)彈性求解法，利用莫爾－庫侖模型計算初始地應(yīng)力，以不平衡力與典型應(yīng)力的比值小于1×10－5為收斂條件判斷初始應(yīng)力平衡，以模型角點（5，5，－250）與除頂面所有面為固定邊界，并固定底面x、y方向速度。

（2）在走向、傾向兩個主要變形線設(shè)置34個監(jiān)測點用于變形量驗證。

（3）煤層開挖，計算步長設(shè)置為2000步，保證不平衡力收斂時足夠小，使模型內(nèi)部應(yīng)力狀態(tài)充分調(diào)整。

計算結(jié)果如下，數(shù)值模擬結(jié)果顯示變形區(qū)主要集中在開采區(qū)上方地表附近，其他地區(qū)明顯變形，開采中心上方形變最大，垂向沉降達30 cm左右，變形特征與煤礦開采沉陷特征相似，中心沉降最大，四周呈梯度減?。▓D17）。

在數(shù)值模擬過程中對走向、傾向剖面線設(shè)置了34個監(jiān)測點用于變形量的檢驗。與D-InSAR結(jié)果對比變形區(qū)內(nèi)南北、東西主要兩條主要剖面如圖18所示，變形量總體變形趨勢一致，但有較小差別，在傾向、走向沉降剖面的差值均方根都在6 cm左右（圖19），在礦區(qū)大變形的條件下6 cm誤差可以接受，根據(jù)數(shù)值模擬對未知采礦區(qū)變形的計算結(jié)果來看，兩者大致趨勢保持一致，變形值相差不大，標準方差較小，計算結(jié)果較為可靠。

3.2.4 野外驗證

圖16 FLAC3D計算模型Fig.16 FLAC3D calculation model

圖17 數(shù)值模擬垂向形變量Fig.17 Numerical simulation of vertical shape variables

我們圍繞采空區(qū)及附近山體進行了野外現(xiàn)場調(diào)查，調(diào)查發(fā)現(xiàn)采空區(qū)上方地表存在多處擠壓和拉張裂縫以及大規(guī)模的巖土體變形塌陷現(xiàn)象，塌陷變形區(qū)主要在開采區(qū)上部，與InSAR計算沉降區(qū)位置基本吻合；在開采區(qū)四周，有多條地面水平運動產(chǎn)生的拉張裂縫（圖20），由地面不均勻水平運動產(chǎn)生，證實了計算結(jié)果的合理性和可靠性。

4 討 論

本文只針對研究區(qū)InSAR觀測周期內(nèi)進行分析討論，煤礦采空造成的地表塌陷影響因素很多，包括地下開采方式、煤層埋深、開采煤層厚度、地層傾角、巖層組合、地表地形、坡面與地層的組合關(guān)系以及開采工作面的布置方向等，且開采沉陷參數(shù)在各個地區(qū)各不相同，故本文只針對貴州貞豐縣地區(qū)對煤層采空區(qū)反演。反演不僅需要該區(qū)地表垂向形變，還需要地層信息、煤層深度、煤層產(chǎn)狀等信息，且開采沉陷參數(shù)獲取來源為該區(qū)附近煤礦，僅對研究區(qū)煤礦有一定應(yīng)用意義。由于該采區(qū)煤礦已經(jīng)廢棄，無法獲得實際地下采空資料，該區(qū)地表也無水準設(shè)備做GPS測量，故本文應(yīng)用數(shù)值模擬技術(shù)進行了驗證，驗證結(jié)果與反演結(jié)果大致相同，一定程度驗證了反演的準確性。另外，由于研究區(qū)植被覆蓋、地物反射特性、季節(jié)因素等，不同SAR數(shù)據(jù)的波段選擇、分辨率、極化方式等皆會對計算結(jié)果產(chǎn)生影響。本文研究區(qū)地表以裸露砂巖或草植被覆蓋砂巖為主，而采用的L波段3 m分辨率數(shù)據(jù)對植被穿透力較強，且地物反射信號較高，計算數(shù)據(jù)結(jié)果較為可靠。

本文誤差來源主要有：D-InSAR計算過程中對較大變形區(qū)的解纏可能會丟失相位信息造成變形量結(jié)果計算的誤差；兩軌D-InSAR三維分解使用的兩對干涉對的時間有一定差別，造成煤礦采空變形周期中某軌形變量的缺失，對三維變形量分解的結(jié)果產(chǎn)生影響；本文采用的開采沉陷參數(shù)為定義的InSAR監(jiān)測地表10imm沉降變形曲線計算出的影響角，未使用以0 mm沉降量邊界計算的移動角，而開采沉陷邊界參數(shù)的選取可能對地下采空區(qū)的反演結(jié)果有一定差異；因煤礦變形一般2 a左右才會達到沉降的最大值，從SAR數(shù)據(jù)量積累和InSAR變形計算所允許時間基線的角度來看，都無法實現(xiàn)礦區(qū)最大沉降量的觀測計算，直接使用《三下采煤規(guī)范》中邊界角概念可能會影響計算準確度，故本文僅能在礦區(qū)地表短期沉降穩(wěn)定后獲取該時段沉降，并計算新定義的煤礦開采沉陷影響角；反演結(jié)果的驗證使用數(shù)值模擬技術(shù)，該區(qū)地層條件比較復(fù)雜，故在數(shù)值模擬中進行了模型簡化，不同地層的物理力學(xué)性質(zhì)不同，可能造成一定的誤差。

圖18 數(shù)值模擬與D-InSAR三維分解對比圖Fig.18 Numerical simulation and D-InSAR three-dimensional decomposition comparison

圖19 數(shù)值模擬與D-InSAR三維分解誤差柱狀圖Fig.19 Numerical simulation and D-InSAR 3D decomposition error histogram

圖20 野外驗證地表變形破壞Fig.20 Field verification of surface deformation and failure

Offset-Tracking技術(shù)可以應(yīng)用于礦區(qū)變形測量中，其應(yīng)用于變形量大、變形范圍大的地區(qū)有較大優(yōu)勢，其計算結(jié)果與窗口選取、濾波技術(shù)的使用有較大關(guān)系，誤差較大，但是本文研究區(qū)變形較小，變形范圍較小，而且本文為L波段3 m高分辨率數(shù)據(jù)測量精度較大D-InSAR計算結(jié)果會更加精確，故未使用Offset-Tracking技術(shù)進行變形量計算。

5 結(jié) 論

本文利用升、降軌觀測16期3 m空間分辨率的L波段PALSAR-2數(shù)據(jù)，采用InSAR三維分解技術(shù)，以貴州省貞豐縣煤礦采空區(qū)為研究區(qū)，開展了多期地表形變監(jiān)測，計算出煤礦采區(qū)地表變形區(qū)分布及沉降值，并得到以下結(jié)論。

（1）利用InSAR三維分解技術(shù)與煤礦地下開采資料對該區(qū)域開采沉陷參數(shù)與變形時間滯后規(guī)律進行計算，得出該地區(qū)開采沉陷參數(shù)（走向影響角83°、上山影響角75°、下山影響角80°）以及約1個月變形滯后時間。

（2）利用InSAR三維分解技術(shù)結(jié)合已知區(qū)計算的開采沉陷參數(shù)與變形滯后時間，對已廢棄的未知煤礦采空區(qū)進行了反演，計算得出廢棄煤礦380 m×150 m的地下采空范圍，與地下開采的大致時間，野外調(diào)查和FLAC3D數(shù)值模擬計算結(jié)果與D-InSAR變形結(jié)果總體較為一致，證實了該方法的合理性和可靠性。

該方法可以在資料足備的情況下，對未知煤礦采空區(qū)進行反演，為計算地下采空范圍及礦山越界開采提供了新思路，并對越界開采起到一定監(jiān)測作用。