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

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

0 引言

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

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

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

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

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

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

圖1 D-InSAR基本原理圖

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

(2)

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

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

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

2.1 研究區(qū)概況

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

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

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

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

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

dmax=λ/μ

(3)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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