孫建，王道順，姚文舉,3*，李進峰，甄洪帥,3，張方龍

(1.山東省魯南地質(zhì)工程勘察院(山東省地質(zhì)礦產(chǎn)勘查開發(fā)局第二地質(zhì)大隊)，山東 兗州 272100；2.高分辨率對地觀測系統(tǒng)濟寧數(shù)據(jù)與應(yīng)用中心，山東 兗州 272100；3.自然資源部采煤沉陷區(qū)綜合治理工程技術(shù)創(chuàng)新中心，山東 兗州 272100；4.梁山縣自然資源和規(guī)劃局，山東 梁山 272100)

0 引言

地下煤炭資源采出后，會導(dǎo)致開采區(qū)周圍巖體原始應(yīng)力平衡發(fā)生破壞，造成應(yīng)力重新分布和覆巖移動，引發(fā)地表沉陷，從而進一步破壞耕地、建(構(gòu))筑物和其他基礎(chǔ)設(shè)施，嚴重制約礦區(qū)和城市的可持續(xù)發(fā)展，而對開采沉陷區(qū)進行實時動態(tài)監(jiān)測則是研究礦區(qū)開采沉陷規(guī)律、預(yù)防災(zāi)害發(fā)生和解決礦區(qū)環(huán)境問題的重要手段，因此開展礦區(qū)地表沉陷監(jiān)測具有重要意義[1-2]。

傳統(tǒng)的礦區(qū)地表沉陷監(jiān)測方法主要通過布設(shè)巖移觀測站，利用GNSS、全站儀、水準儀等常規(guī)測量手段定期進行沉陷觀測，此類方法監(jiān)測精度雖然可觀，但監(jiān)測點位少、布設(shè)的監(jiān)測點容易遭到破環(huán)，難以整體反映礦區(qū)開采的地表沉陷特征，此外如果在觀測條件惡劣的地區(qū)還容易對測量人員構(gòu)成安全威脅[3]。合成孔徑雷達差分干涉測量技術(shù)(Differential Interferometric Synthetic Aperture Radar, D-InSAR)作為一種新興的空間對地觀測技術(shù)，通過非接觸測量的方式克服了傳統(tǒng)監(jiān)測方法的不足，可全天候連續(xù)獲取大范圍區(qū)域的微小形變信息。2004年，吳立新[4-5]在國內(nèi)首次分析了D-InSAR技術(shù)在煤礦開采沉陷變形監(jiān)測中的優(yōu)勢，并將該技術(shù)成功應(yīng)用于唐山市及開灤礦區(qū)。2008年，劉廣等人[6]利用ENVISAT和JERS-1雷達數(shù)據(jù)對河北峰峰礦區(qū)進行了地表沉陷監(jiān)測，分析了不同波段條件下干涉相位特性。2014年，王志勇等人[7]以濟寧礦區(qū)為研究對象，采用ALOS PALSAR雷達數(shù)據(jù)對地下采煤引起的地表沉降進行了精細化分析與解譯，證明D-InSAR在礦區(qū)地表開采沉陷監(jiān)測中的精度已達到cm級。2017年，郭山川等人[8]利用兩軌法D-InSAR技術(shù)對黃土高原的大同市南郊區(qū)進行了采煤沉陷監(jiān)測，并與開采沉陷預(yù)計結(jié)果進行了驗證，分析了結(jié)果差異原因，此后眾多學(xué)者不斷將D-InSAR技術(shù)應(yīng)用于礦區(qū)地表開采沉陷監(jiān)測研究當中[9-14]。

本文在分析雙軌D-InSAR獲取地表形變原理的基礎(chǔ)上，對兩景高分辨率RADARSAT-2超寬精細(Extra Fine,XF)模式雷達影像進行差分干涉處理，獲取了霄云礦區(qū)2016年12月18日—2017年2月25日開采沉陷區(qū)域和沉陷值，對監(jiān)測結(jié)果和精度進行了分析，證明了D-InSAR技術(shù)在礦區(qū)地表開采沉陷監(jiān)測中的有效性。

1 D-InSAR技術(shù)原理

D-InSAR技術(shù)是通過對同一地區(qū)前后跨越形變期獲取的SAR單視復(fù)數(shù)影像干涉得來的相位差，進而獲取地面形變信息的技術(shù)。根據(jù)在差分干涉處理過程中去除地形相位方法的不同，可以將D-InSAR測量形變分為二軌法、三軌法、四軌法等[15]。二軌法所需數(shù)據(jù)量少，干涉處理簡單，依靠外部DEM數(shù)據(jù)模擬去除地形相位，由此獲取地表形變信息。隨著高精度覆蓋全球陸地地表數(shù)字高程模型(Digital Elevation Model, DEM)數(shù)據(jù)的免費使用，二軌法已經(jīng)成為D-InSAR獲取地表形變的常用方法。

圖1為D-InSAR技術(shù)的基本原理圖。在雷達兩次成像前后，地物點P發(fā)生形變后的位置為P'，形變位移矢量為△d。在不考慮其他相位誤差的條件下，前后兩景影像的干涉相位可表示為：

圖1 D-InSAR基本原理圖

由于干涉條件的限制，D-InSAR的干涉相位φint中通常包含不同的相位成分，其中φref為參考相位，φtopo為地形相位，將這兩項從干涉相位去除后即可得到形變相位φdef，因此可以得到：

(2)

由式(2)可以看出，D-InSAR技術(shù)監(jiān)測得到的形變量是在雷達視距方向的投影，數(shù)據(jù)處理流程見圖2。

圖2 D-InSAR數(shù)據(jù)處理流程圖

2 研究區(qū)與數(shù)據(jù)源

2.1 研究區(qū)概況

霄云煤礦位于山東省濟寧市金鄉(xiāng)縣霄云鎮(zhèn)境內(nèi)，隸屬于濟寧礦業(yè)集團有限公司。礦井于2012年8月建成投產(chǎn)，可采儲量4791.5萬t，主采3煤賦存穩(wěn)定，可采范圍內(nèi)平均厚度3.49m，傾角一般為10°～21°，采用走向長壁采煤法，綜合機械化一次采全高采煤工藝。截至2017年3月，研究區(qū)共包括7個工作面，工作面開采情況見表1，區(qū)域概況見圖3。

表1 研究區(qū)工作面信息

圖3 研究區(qū)域概況圖

2.2 數(shù)據(jù)源

礦山開采沉陷具有大量級、非線性等特點，當?shù)乇硇巫兲荻容^大，超過InSAR技術(shù)所能測量的梯度范圍時，容易造成干涉圖干涉條紋雜亂和相位混疊，導(dǎo)致解纏錯誤，從而難以獲取準確的開采沉陷信息[16]。InSAR監(jiān)測的最大形變梯度可表示為[17]：

dmax=λ/μ

(3)

式中：λ—雷達波長，μ—影像像元分辨率?？梢钥闯?，InSAR可監(jiān)測的最大形變梯度與所用雷達數(shù)據(jù)參數(shù)相關(guān)，雷達波長越長、像元分辨率越高，能監(jiān)測到的形變梯度越大，越適用于礦區(qū)開采沉陷監(jiān)測，因此本次實驗選擇加拿大RADARSAT-2 XF的C波段5m高分辨率SAR影像作為監(jiān)測數(shù)據(jù)源，該影像在礦區(qū)開采沉陷監(jiān)測中具有顯著優(yōu)勢[18]。外部DEM數(shù)據(jù)選擇使用日本太空發(fā)展署(JAXA)免費分發(fā)的30m分辨率ALOS World 3D數(shù)據(jù)用以模擬并去除地形相位。所用影像對參數(shù)見表2。

表2 干涉SAR影像參數(shù)列表

3 數(shù)據(jù)處理與分析

3.1 數(shù)據(jù)處理過程

雙軌D-InSAR方法可以快速獲取監(jiān)測區(qū)域的形變信息，但需要保證監(jiān)測區(qū)域具有良好的時間和空間相干性，若所用數(shù)據(jù)的空間基線超過臨界基線值，會導(dǎo)致完全失相干，無法獲取有效監(jiān)測信息，因此在數(shù)據(jù)處理前一般需要對所用數(shù)據(jù)進行基線估計，以便判斷所用數(shù)據(jù)能否達到理想的監(jiān)測結(jié)果。由表2可知，兩景影像時間基線為72d，空間基線僅為-23.96m，遠遠小于其空間臨界基線值，滿足影像干涉要求。

實驗選擇2016年12月18日影像作為主影像，2017年2月25日獲取的影像為輔影像，根據(jù)圖2的數(shù)據(jù)處理流程進行差分干涉處理。首先對兩景影像進行配準處理，并保證配準精度小于1/8個像元，生成原始干涉圖(圖4a)；在距離向和方位向按照1∶2的視數(shù)比進行多視處理，去除平地相位和地形相位后生成差分干涉圖(圖4b)；由于數(shù)據(jù)成像和處理過程中不可避免會受到噪聲影響，實驗選擇自適應(yīng)的Goldstein濾波方法對差分干涉圖進行濾波處理(圖4c)，減弱噪聲相位，并生成相干系數(shù)圖(圖4d)；為了恢復(fù)真實相位值，利用最小費用流方法進行對干涉相位進行解纏處理，設(shè)置相干解纏閾值0.2，逐像元進行相位解纏。

圖4 差分干涉處理結(jié)果圖

在遠離形變區(qū)域、無相位解纏錯誤且相干性良好的區(qū)域選擇19個控制點，利用三次多項式模型去除殘余的恒定相位和相位坡道，軌道精煉后的相位中誤差為-0.0057rad，標準差為0.1614rad；去除殘余相位后將相位換算為形變，經(jīng)過地理編碼后生成監(jiān)測區(qū)域的沉陷結(jié)果(圖5)。

圖5 研究區(qū)沉陷量圖

3.2 監(jiān)測結(jié)果分析

從圖5可以看出，由于煤炭開采活動影響，研究區(qū)域在2016年12月18日—2017年2月25日內(nèi)發(fā)生明顯沉陷現(xiàn)象，共分布有3處沉陷盆地，沉陷盆地空間分布特征與采空區(qū)分布實際情況吻合，證明了D-InSAR技術(shù)在礦區(qū)地表沉陷監(jiān)測中的有效性。

為了進一步對監(jiān)測結(jié)果進行分析，實驗設(shè)置沉陷等值線初始值-0.005m，間距0.01m，繪制沉陷等值線見圖6。

圖6 研究區(qū)沉降等值線圖

可以看出，西側(cè)沉陷盆地沉陷范圍最大，中間位置的沉陷盆地沉陷量最大，東側(cè)沉陷盆地沉陷量最小。結(jié)合圖3采空區(qū)情況分析，a區(qū)域監(jiān)測時間段內(nèi)1311工作面自西向東開采，北側(cè)1309工作面為已停采區(qū)，1311工作面的采動引起了1309已采區(qū)上覆巖層“活化”而產(chǎn)生二次破裂，造成沉陷盆地中心位于兩工作面中間位置且東部沉陷量明顯大于西部，最大沉陷量為0.118m，而1309工作面自西向東開采停采5個月，處于沉陷衰退階段，東側(cè)區(qū)域沉陷量小于0.015m；b區(qū)域兩個沉陷盆地已接連成片，1304工作面的開采同樣引起了1302-1、1302-2和1306已采區(qū)的二次變形并使沉陷量和沉陷范圍擴大，最大沉陷量達到0.174m，1312工作面自東向西開采，僅停采2個月，最大沉陷值-0.064m，因此仍能呈現(xiàn)出沉陷盆地形態(tài)。

值得注意的是，根據(jù)研究區(qū)的煤層開采厚度和其它地質(zhì)采礦信息，本次研究監(jiān)測到的工作面開采下沉量顯著小于理論值。分析其原因認為，本次實驗監(jiān)測的時間區(qū)間為2016年12月18日至2017年2月25日(共69天)，未完整覆蓋1311和1304工作面的全部開采時段，監(jiān)測時段之后地表仍處于移動變形活躍期，因此監(jiān)測實驗得到的整體下沉值偏小。由于缺少實測數(shù)據(jù)，本次實驗采用內(nèi)符合評估的方式對監(jiān)測結(jié)果精度進行評價，通過統(tǒng)計形變量及其中誤差的大小和分布來驗證監(jiān)測結(jié)果的可靠性。

從表3、圖7可以看出，監(jiān)測結(jié)果的沉陷量峰值為-0.001m，沉陷分布接近于0m，符合大范圍穩(wěn)定的實際情況；沉陷精度峰值-0.026m，平均值0.011m，符合D-InSAR技術(shù)的cm級監(jiān)測精度要求。

表3 精度統(tǒng)計表

圖7 監(jiān)測沉降量統(tǒng)計分布圖

4 結(jié)論

本文利用D-InSAR技術(shù)對覆蓋研究區(qū)域的兩景RADARSAT-2雷達數(shù)據(jù)進行差分干涉處理，獲取了該區(qū)域2016年12月18日至2017年2月25日時間段內(nèi)的沉降信息，并對監(jiān)測結(jié)果和精度進行了分析。

(1)D-InSAR監(jiān)測到研究區(qū)域共分布有3處沉陷盆地，監(jiān)測到的最大沉陷量達到-0.174m。

(2)研究區(qū)域的監(jiān)測結(jié)果沉陷分布特征與工作面開采實際位置高度吻合，沉陷量大小和沉陷范圍與工作面開采情況相符。

(3)監(jiān)測區(qū)域沉陷量的統(tǒng)計分布峰值趨近于0m，內(nèi)符合精度誤差平均值0.011m，整體沉陷趨勢符合實際情況，證明D-InSAR技術(shù)在礦區(qū)地表沉陷監(jiān)測中具有良好的應(yīng)用效果。