王劍,董祥林,楊可明,姚樹一,石曉宇
(1.中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院,北京,100083;2.淮北礦業(yè)股份有限公司通防地測部,安徽淮北,235000)
地下采礦活動引起的地表大梯度快速沉降會導(dǎo)致上覆建筑物、管線以及道路等人工設(shè)施發(fā)生損壞[1],因此,常采取相應(yīng)科學(xué)有效的采礦工藝或開采措施以減弱地表下沉帶來的危害,確保礦區(qū)開采的可持續(xù)性和安全性。同時,也必須提供精確可靠的地表沉降監(jiān)測方法來獲取礦區(qū)地表下沉盆地范圍以及地表沉降量。差分干涉合成孔徑雷達測量(differential interferometry synthetic aperture radar,DInSAR)作為新興的空間對地觀測技術(shù),可以敏感地獲取地表微小變化,其監(jiān)測精度可達厘米級,甚至毫米級,并且具有全天候、全天時、大范圍便捷監(jiān)測的優(yōu)勢[2]。DInSAR 技術(shù)及由其發(fā)展而來的時序InSAR 技術(shù)已成為監(jiān)測礦區(qū)地表沉降的有效熱門手段之一[3-5]。目前,國內(nèi)外學(xué)者利用InSAR 技術(shù)監(jiān)測礦區(qū)地表沉降主要有3 種:PS(permanent scatterers,永久散射體)、SBAS(small baseline subsets,小基線集)等時序InSAR 技術(shù)[6-7],POT(Pixel offset-tracking,像素偏移量追蹤)技術(shù)[8-9]以及DInSAR 技術(shù)[10-11]。雖然時序InSAR 技術(shù)對形變梯度較小的地表形變具有很高的監(jiān)測精度[12],但得到的結(jié)果不連續(xù),只能獲取部分高相干點目標的形變信息[13],并且在對礦區(qū)地表大梯度沉降進行監(jiān)測的過程中存在結(jié)果明顯偏低的缺陷[14-15]。POT技術(shù)雖然可以監(jiān)測到地表較大沉降量[16],但數(shù)據(jù)處理時間較長,且結(jié)果精度低,僅為像元分辨率的1/10~1/30[14,17]。而DInSAR 技術(shù)處理時間短,所需影像少,得到的地表形變結(jié)果連續(xù),可監(jiān)測的地表沉降最大梯度適中,得到的結(jié)果精度較高,更適用于監(jiān)測礦區(qū)地表沉降[18]。此外,利用上述InSAR技術(shù)獲取煤礦區(qū)地表沉降信息研究大多數(shù)是針對開采的煤層上覆巖層在自然垮落情況下的礦區(qū)地表沉降監(jiān)測,對開采工作面注漿充填影響下的礦區(qū)地表動態(tài)沉降監(jiān)測研究很少[6,8-11];且大多采用的是商業(yè)L、X波段數(shù)據(jù)[3,17-18],若要得到地表長期的動態(tài)沉降信息則所需成本較高。因此,本文作者利用哨兵1 號A 星(Sentinel-1A)的C 波段影像數(shù)據(jù),以淮北礦業(yè)集團臨渙煤礦Ⅱ2 采區(qū)1021 和Ⅱ1021 工作面采動影響范圍為試驗區(qū),采用時序疊加DInSAR (time series stacking DInSAR, TSSDInSAR)方法并結(jié)合現(xiàn)場實測水準數(shù)據(jù)來監(jiān)測分析工作面開采時在注漿充填影響下的礦區(qū)地表動態(tài)沉降,以探尋注漿充填的地表減沉效果以及注漿充填影響下地表動態(tài)沉降規(guī)律,同時為礦區(qū)提供大范圍、低成本的地表長期動態(tài)沉降監(jiān)控分析新手段。
TSS-DInSAR 方法是建立在二軌法DInSAR 基礎(chǔ)上的。二軌法DInSAR 最早由MASSONNET等[19]提出,其主要思想為利用形變前后各一景SAR影像以及參考DEM數(shù)據(jù),通過差分處理生成包含不同組分相位信息的差分干涉圖,然后將除地表真實形變相位之外的其他相位信息剔除,即:


式中:λ為合成孔徑雷達波長。由于垂直沉降是礦區(qū)地表變形的主要貢獻分量[20],故忽略水平移動分量對雷達視線方向形變的貢獻,由三角函數(shù)關(guān)系得地表垂直方向的下沉量W為[6,20-21]

式中:θ為雷達成像幾何入射角。
對按時序獲取的同一位置的n 景影像A=[a1,a2,…,an](各影像對應(yīng)的獲取時間分別為t1,t2,…,tn,且滿足t1<t2<…<tn),每相鄰兩景重復(fù)進行二軌法DInSAR處理,得到每相鄰2個影像時段內(nèi)的地表沉降,即a1與a2得到w1~2,a2與a3得到w2~3,…,an-1與an得到wn-1~n,從而有n-1期地表下沉量:

然后,將這n-1期地表沉降結(jié)果與實地地物進行嚴格配準,并以第1 期的地表沉降結(jié)果w1~2為基準按時序依次疊加,即w1~3=w1~2+w2~3,w1~4=w1~2+w2~3+w3~4,…,w1~n=w1~2+w2~3+…+wn-1~n,從而得到從t1到tn時段內(nèi)的地表時序累積下沉量:

基于式(4)和式(5)結(jié)果分析在注漿充填影響下的礦區(qū)地表動態(tài)沉降規(guī)律。
選取中國安徽省北部淮北市境內(nèi)的淮北礦業(yè)集團臨渙煤礦Ⅱ2 采區(qū)作為研究區(qū),其地表覆蓋較多的植被、農(nóng)作物,地理位置如圖1所示,研究區(qū)內(nèi)布置有Ⅱ923 里、1021 和Ⅱ1021 共3 個采煤工作面,均使用綜合機械化長壁垮落法結(jié)合覆巖關(guān)鍵層注漿充填開采,詳細參數(shù)見表1。圖1 中,黑色多邊形區(qū)域表示工作面,在研究時段(2016-11-21—2017-02-25)內(nèi),有注漿活動的注漿孔共有10個,用同心環(huán)狀符號表示。

圖1 研究區(qū)內(nèi)的工作面與注漿孔位置Fig.1 Working-face and grouting-hole positions in study area
經(jīng)大量試驗發(fā)現(xiàn),受研究區(qū)地表濃密植被的影響,春、夏、秋季得到的C波段雷達影像差分干涉對相干性很低,很大程度上影響了DInSAR正確提取地表形變相位,導(dǎo)致結(jié)果精度較差[22],而且礦區(qū)地表大梯度沉降會導(dǎo)致相位缺失,最短重復(fù)周期的SAR 影像對可以使時間失相關(guān)減至最小,因此,本文在盡量滿足時間基線最短的前提下獲取研究區(qū)2016-11-21—2017-02-25 的8 景Sentinel-1A影像,影像數(shù)據(jù)主要參數(shù)見表2,組成的7個干涉對及其相關(guān)信息如表3 所示??梢姡撼?017-01-08—2017-02-01 的影像對時間基線為24 d 外,其余影像對的時間基線均為Sentinel-1A 最短重復(fù)周期12 d。用于與研究結(jié)果作對比、驗證的現(xiàn)場實測水準數(shù)據(jù)由研究區(qū)煤礦提供,其時間范圍為2016-11-16—2017-02-26,比研究時段多6 d,地面水準測量點數(shù)為40個。
利用式(1)~(4)對獲取的8 景Sentinel-1A 影像進行DInSAR處理,得到7期地表下沉速率,結(jié)果如圖2所示,對其統(tǒng)計得到單期的最大下沉量Wmax與最大下沉速率Vmax見表4。由圖2 可知:在研究時段內(nèi),研究區(qū)地表形成了明顯的下沉盆地,與工作面位置相吻合;隨著工作面的推進,每期地表下沉速率最大點也向工作面推進方向移動,并且符合在近水平煤層開采條件下,地表沉降具有超前影響以及最大下沉速率點滯后于開采活動中心的規(guī)律[1]。結(jié)合表4 可知:1)第5 期的地表最大下沉量Wmax最大,達47 mm,但其最大下沉速率Vmax在監(jiān)測時段內(nèi)最小,為1.96 mm/d,這是由于第5 期的DInSAR 時間基線為24 d,在此期間內(nèi)發(fā)生的地表沉降應(yīng)該遠大于此期監(jiān)測結(jié)果,但其實際下沉量超出了哨兵1號數(shù)據(jù)的最大監(jiān)測能力,因而只監(jiān)測到了47 mm的下沉量;2)第1期的地表最大下沉量最小,僅為30 mm;3)第3期的地表最大下沉速率最大,為3.50 mm/d。綜合各期數(shù)據(jù)分析可得,平均地表最大下沉量為38 mm,平均最大下沉速率為2.89 mm/d。

表1 研究區(qū)工作面參數(shù)Table 1 Parameters of working faces in study area

表2 影像數(shù)據(jù)主要參數(shù)Table 2 Main parameters of data

表3 干涉對組成及基線信息Table 3 Composition of interferometry pairs and baseline information

圖2 地表時序下沉速率Fig.2 Time series subsidence velocity of surface
結(jié)合實際注漿資料,第4期與第5期之間經(jīng)歷了較劇烈的注漿活動,因而相較于前幾期,第5期地表沉降中心區(qū)下沉速率明顯降低,并且各期最大沉降速率圍繞平均最大沉降速率波動較為劇烈,在第3期到達形變速率峰值后開始減小,至第5期減小為形變速率低谷,而第6期又到達形變速率次高峰,表明受注漿充填活動影響的采區(qū)地表沉降與無注漿充填也即巖層自然垮落引起的采區(qū)地表沉降規(guī)律有較大差異(正常情況下單點下沉速率呈現(xiàn)類似二次函數(shù)拋物線形狀)。另外,因受到注漿充填的影響,即使有些建筑物處于開采中心區(qū),相對于未受到注漿充填影響的地表,其中心及其附近地表下沉速率也要小很多,因而也使得DInSAR監(jiān)測的每期地表下沉速率圖以開采工作面走向線為中軸線表現(xiàn)出明顯的非對稱性(見圖2),這些結(jié)果均表明注漿減沉效果顯著。
為研究礦區(qū)地表長期動態(tài)沉降,將7幅地表沉降圖所示的階段沉降結(jié)果按時間先后順序依次進行疊加,得到地表時序累積沉降圖,結(jié)果如圖3所示,對其統(tǒng)計得到的時序累積最大下沉量與最大下沉速率見表4。從圖3 可以得出:在監(jiān)測期間,
地表下沉盆地的范圍隨時間推移有逐步擴大的跡象,主要表現(xiàn)在1021 與Ⅱ1021 工作面新開采的區(qū)域,并且在這2個工作面開采活動中心地表沉降異常明顯,形成了2個沉降中心,下沉量最大的沉降中心位于Ⅱ1021工作面采空區(qū)上覆地表,累積最大下沉量達到18 5 mm,另一個沉降中心位于1021工作面西北側(cè)地表,累積最大下沉量達到169 mm,以下沉量10 mm 為下沉邊界起算點,最終下沉盆地面積約為1.35 km2。由于Ⅱ923 里工作面已于2016-04-28 收作,距離開始時間(2016-11-21)已達7月,故II923里工作面上覆地表僅在靠近2個開采工作面的區(qū)域有較小沉降,最大下沉量為42 mm。研究區(qū)地表最終累積最大下沉速率為1.93 mm/d,平均累積最大下沉速率為2.18 mm/d,均大于地表移動持續(xù)時間三階段劃分閾值1.67 mm/d,表明該研究時段內(nèi)地表下沉盆地處于下沉活躍期(開始期為從下沉10 mm 開始至下沉速率達到1.67 mm/d;活躍期為下沉速率大于1.67 mm/d;衰退期為下沉速率小于1.67 mm/d 且6月內(nèi)下沉量小于30 mm)。

表4 地表時序沉降統(tǒng)計Table 4 Statistics of time series subsidence of surface

圖3 地表時序累積沉降Fig.3 Time series cumulative subsidence of surface
基于圖3并結(jié)合實際注漿資料分析,由于1021工作面及其東南側(cè)儲灰場在研究時段內(nèi)持續(xù)實施注漿充填,故位于1021 工作面采空區(qū)上方的地面按無充填開采沉陷規(guī)律來說,應(yīng)該處于沉降中心區(qū),但事實上其下沉量比未受到注漿減沉影響的1021 工作面西北側(cè)區(qū)域下沉量小;另外,研究區(qū)地表下沉也無明顯向東南側(cè)擴張的跡象,下沉盆地?zé)o論是下沉范圍還是下沉量都表現(xiàn)出顯著的非對稱性,這些均有力表明注漿充填減緩地面沉降的效果顯著。
為進一步量化了解研究區(qū)地表動態(tài)沉降規(guī)律,對下沉盆地主斷面上的沉降特征進行分析。在研究區(qū)1021和Ⅱ1021工作面上選取了經(jīng)過2個工作面累積沉降最大值中心的線段AA1近似作為走向線,經(jīng)過Ⅱ1021 工作面累積下沉量最大值中心的線段BB1作為傾向線1,經(jīng)過1021 工作面累積下沉量最大值中心的線段CC1作為傾向線2,如圖4 所示。基于這3條線段提取出3組位于下沉盆地主斷面上的時序累積沉降曲線,結(jié)果如圖5 所示(其中坐標橫軸表示以A,B,C 為起點,主斷面上各點至起點的距離)。

圖4 所選主斷面位置示意圖Fig.4 Location map of selected main sections
1)在AA1的走向主斷面,以下沉量10 mm為下沉邊界起算點,則下沉盆地沉降范圍為距離A 點330~1 970 m 之間,記為A-330~1 970 m,長約1 640 m。Ⅱ1021 工作面的最大沉降中心和1021 工作面西北側(cè)的最大沉降中心分別位于850 m附近和1 350 m 附近,對應(yīng)的累積最大下沉量分別為185 mm和169 mm。另外,在800 m處的地表本應(yīng)處于最大下沉中心區(qū),卻存在凸起現(xiàn)象,結(jié)合實地注漿資料以及水準測量數(shù)據(jù)分析,這是由于該處緊鄰Ⅱ1021 工作面2 號注漿孔,而該注漿孔在DInSAR 第4 期與第5 期之間進行過劇烈的注漿活動,因此,地表出現(xiàn)了一定程度的抬升,進而使5~7 期時序累積沉降曲線在此處均存在局部峰值,與緊鄰此注漿孔的實地水準測點數(shù)據(jù)吻合。由于950~1 150 m 范圍內(nèi)的盆地距離2 個開采活動中心較遠,因而,其下沉量比最大沉降中心的下沉量小。1 500~1 900 m 范圍內(nèi)的地表升降波動也由第4~5期之間的地表注漿活動引起。
2)在BB1的傾向主斷面,以下沉量10 mm為下沉邊界起算點,則其下沉盆地范圍為B-470~1 570 m,長約1 100 m。在觀測期間,Ⅱ1021 工作面地表最大沉降中心從第1期的16 mm下沉量,以平均下沉速率2.01 mm/d 下沉到最后一期的185 mm,1 100 m 附近的5~7 期地表升降波動也由地表注漿充填引起。

圖5 下沉盆地主斷面時序累積沉降曲線Fig.5 Time series cumulative subsidence curves of main sections of subsidence basin
3)在CC1的傾向主斷面,以下沉量10 mm為下沉邊界起算點,則其下沉盆地范圍為C-500~1 620 m,長約1 120 m,與BB1主斷面的沉降長度接近。在觀測期間,1021 工作面西北側(cè)地表最大沉降中心從第1期的20 mm下沉量,以平均下沉速率1.76 mm/d 下沉到最后一期的169 mm。在5~7期,出現(xiàn)在1 300 m 附近的曲線峰值與1 500 m 附近的曲線低谷同樣是由于注漿活動引起的,并且,離最大沉降中心更近的1 300 m附近的地表下沉量小于1 500 m附近的地表下沉量,進一步證明了注漿減沉效果顯著。
由BB1傾向主斷面與CC1傾向主斷面對比可得,雖然Ⅱ1021工作面屬于重復(fù)采動并且傾向長度要比1021工作面的大,但因其為新開采的工作面,地表沉降有一定的延遲時間,即存在起動距,因而在前4期監(jiān)測結(jié)果中,Ⅱ1021工作面最大沉降中心的下沉量一直比1021 工作面西北側(cè)最大沉降中心的小,但隨著采掘工作面繼續(xù)向前推進,從第5期起,Ⅱ1021 工作面采空區(qū)走向長度達到155 m,Ⅱ1021 工作面最大沉降中心的下沉量大于1021 工作面西北側(cè)最大沉降中心的下沉量。
總體上看,這3組時序累積沉降曲線具有雙工作面同時開采的地表下沉特征,并且地表最大下沉量均顯示出增加的趨勢,表明研究區(qū)未達到充分采動狀態(tài)。另外,上述結(jié)果也間接反映了若研究區(qū)地表未進行注漿減沉,則下沉盆地范圍與下沉最大值皆會增大。
為了驗證利用TSS-DInSAR方法來監(jiān)測分析礦區(qū)地表動態(tài)沉降規(guī)律的可靠性,采用了地面實測水準數(shù)據(jù)進行對比驗證,水準測點位置見圖6。因本研究時段為2016-11-21—2017-02-25,而水準數(shù)據(jù)監(jiān)測時段為2016-11-16—2017-02-26,二者時間間隔并不完全相同,為了盡可能減少因時間基準不一致帶來的誤差,根據(jù)各位置時間段內(nèi)地表下沉速率等相關(guān)信息將水準數(shù)據(jù)起止日期內(nèi)插至與研究時段相同的日期,即水準數(shù)據(jù)與TSSDInSAR監(jiān)測結(jié)果起止時間一樣。

圖6 水準驗證數(shù)據(jù)點位Fig.6 Point positions on verification data of leveling
TSS-DInSAR與水準測量的地表累積下沉量對比如圖7 所示。由圖7 可知:除了5~8 點差值較大(差值分別為-109,-96,-121 和-75 mm)以及38號,39 號點曲線趨勢不一致外,其余點位差值都比較小,擬合度較好,均方根誤差RMSE 為35.8 mm。這是由于5~8 號點位于下沉盆地中下沉梯度最大的位置,以4號和5號點為例,其下沉梯度約為2.0 mm/m,這遠超出C 波段Sentinel-1A 數(shù)據(jù)可監(jiān)測到的最大下沉梯度1.4 mm/m[21],導(dǎo)致這些點的DInSAR 監(jiān)測結(jié)果與水準測量結(jié)果差值較大。此外,對于39號點,其位置比38號點更靠近最大沉降中心,但水準測量數(shù)據(jù)卻出現(xiàn)抬升異常,綜合研究區(qū)的多次DInSAR處理試驗結(jié)果,該異??赡苁怯伤疁蕼y量中人為誤差傳遞導(dǎo)致。于是,以2.5倍RMSE為限差,剔除觀測相位失相干點5~7 后,得到新的RMSE 為19.7 mm,精度提高了45%。
為進一步驗證本研究精度,將去除失相干點5~7 后的TSS-DInSAR 累積下沉量與水準數(shù)據(jù)累積下沉量進行回歸分析,并進行了F檢驗,結(jié)果如圖8所示。從圖8可知:決定系數(shù)R2達到0.910 4,二者之間的標準化殘差值95%分布于-2~2,且服從正態(tài)分布;F′=355.53,F(xiàn)=F0.05(1,35)=4.13,F(xiàn)′>F,且p<0.001。這些結(jié)果表明:去除失相干點后的最終DInSAR累積下沉量與水準測量累積下沉量之間具有很強的線性關(guān)系,回歸模型具有高擬合優(yōu)度與置信度,利用TSS-DInSAR方法來監(jiān)測分析礦區(qū)地表動態(tài)沉降規(guī)律是切實可靠的,滿足實際應(yīng)用需求。

圖7 TSS-DInSAR與水準測量的累積下沉量對比Fig.7 Comparison of TSS-DInSAR cumulative subsidence and leveling cumulative subsidence
沉陷中心區(qū)域TSS-DInSAR監(jiān)測結(jié)果與水準測量結(jié)果差值較大并且研究中存在誤差,其產(chǎn)生原因主要有:
1) 由于Sentinel-1A 數(shù)據(jù)波長較短、分辨率較低,能夠監(jiān)測到的最大下沉梯度遠小于研究區(qū)地表實際的最大下沉梯度,因此在下沉梯度較大的地區(qū)出現(xiàn)相位丟失現(xiàn)象,難以探測出準確的下沉量。

圖8 TSS-DInSAR與水準累積下沉量回歸分析Fig.8 Regression analysis on TSS-DInSAR cumulative subsidence and leveling cumulative subsidence
2)地面實測水準數(shù)據(jù)起止時段與DInSAR監(jiān)測起止時段有些差異,在一定程度上會影響數(shù)據(jù)擬合驗證分析。
3)單元范圍不一致。DInSAR 技術(shù)獲取的是1個面的結(jié)果,而水準測量針對的則是單點。具體而言,DInSAR技術(shù)測量的結(jié)果是影像上一個分辨單元的沉降量,受影像分辨率的限制,地表上某點的沉降情況會受到分辨單元范圍內(nèi)所有點的影響。因此,用1個“點”來評價1個“面”本身就存在著誤差[2]。
4)緊鄰研究區(qū)西邊的沉陷積水區(qū),并且所在區(qū)域覆蓋大量植被,這降低了DInSAR干涉對的相干性,進而影響解纏結(jié)果與最終提取的地表沉降結(jié)果。
5)未考慮水平移動分量對雷達視線方向形變的貢獻,直接通過投影關(guān)系將視線向形變轉(zhuǎn)換至垂直向形變,這必然會帶來誤差。
6)其他因素。包括所采用的DInSAR處理算法本身存在誤差、大氣延遲效應(yīng)等因素。
綜合考慮以上誤差來源,在仍使用高性價比的Sentinel-1 衛(wèi)星影像作為數(shù)據(jù)源的前提下,減小上述誤差的方式主要有:
1) 引入Sentinel-1B 衛(wèi)星影像數(shù)據(jù),將時間分辨率提高1倍,從而可以顯著減少因礦區(qū)地表形變梯度過大而導(dǎo)致的時空失相干,并且能更精細的提取礦區(qū)地表形變。
2)優(yōu)化DInSAR處理算法,比如增大圖像配準窗口大小以及配準窗口數(shù)量來提高配準精度,對得到的干涉圖進行更為精細的濾波以盡量保證下沉盆地中的干涉相位連續(xù)且下沉盆地外的噪聲較少,相位解纏可使用基于不規(guī)則格網(wǎng)的方法來減少MCF解纏結(jié)果存在較多的相位不連續(xù)現(xiàn)象。
3)在干涉時盡量避免雨雪天氣所獲取的影像數(shù)據(jù),以提高得到的干涉圖的質(zhì)量并減少大氣延遲相位的組分,另外,可以借助外部GPS 數(shù)據(jù)來減弱大氣延遲誤差。
4)反演出地表三維形變,針對某個方向上的真實形變值進行對比。
1)利用TSS-DInSAR方法監(jiān)測到的研究區(qū)地表下沉盆地在研究時段內(nèi)處于下沉活躍期,其平均最大下沉量為38 mm,平均最大下沉速率為2.89 mm/d,最終累積最大下沉速率為1.93 mm/d,平均累積最大下沉速率為2.18 mm/d,最終下沉盆地面積約為1.35 km2,且在Ⅱ1021工作面和1021工作面西北側(cè)形成了2個沉降中心,二者的最終累積最大下沉量分別為185 mm和169 mm。
2)受注漿充填影響的采區(qū)建筑物及其附近地表的下沉速率、下沉量以及下沉范圍均顯著減小。下沉盆地以工作面走向線為中軸線表現(xiàn)出明顯的非對稱性,表明結(jié)合注漿充填的采煤工藝可以有效減緩控制區(qū)域內(nèi)的地表沉降,有注漿充填活動的采區(qū)地表沉降與巖層自然垮落引起的采區(qū)地表沉降規(guī)律存在較大差異。
3)TSS-DInSAR監(jiān)測結(jié)果與地面實測水準數(shù)據(jù)就最終累積下沉量進行驗證精度,剔除3個相位失相干點后,RMSE 為19.7 mm,回歸分析R2為0.910 4,F(xiàn) 檢驗p<0.001,表明基于TSS-DInSAR方法的礦區(qū)地表沉降監(jiān)測結(jié)果可靠性較高,具有實際應(yīng)用可行性。同時也反映出Sentinel-1A 數(shù)據(jù)能監(jiān)測到地表微小沉降且測量精度很高,而對下沉梯度大于1.4 mm/m 的地表沉降,其監(jiān)測值會偏低。