叢康林, 岳建平, 董 超, 方云波

(1. 山東農(nóng)業(yè)大學(xué) 信息科學(xué)與工程學(xué)院，山東 泰安 271018；2. 河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 211100；3. 上海衛(wèi)宸測(cè)繪科技有限公司，上海 200001)

采礦區(qū)地面沉降是一種變化緩慢且不可逆轉(zhuǎn)的地表形變現(xiàn)象，會(huì)造成建筑物地基下沉、房屋開裂、地下管道破損等一系列問題，成因機(jī)制錯(cuò)綜復(fù)雜，致災(zāi)過程緩慢，一旦形成便難以恢復(fù)。新汶礦區(qū)作為百年老礦區(qū)，存在采煤時(shí)間久、采空區(qū)范圍大、地表沉降嚴(yán)重等問題。常規(guī)地面大地測(cè)量技術(shù)和地質(zhì)監(jiān)測(cè)手段雖能夠準(zhǔn)確獲取地面沉降結(jié)果，但該類方法是通過離散的監(jiān)測(cè)點(diǎn)獲取形變信息，無法得到大范圍連續(xù)面狀區(qū)域的地表形變信息，存在一定的局限性。

合成孔徑雷達(dá)干涉測(cè)量(Interferometric Synthetic Aperture Radar，InSAR)作為一種影像測(cè)地技術(shù)，以其全天候、全天時(shí)、高精度和區(qū)域測(cè)量的特點(diǎn)，可獲取大范圍面狀連續(xù)區(qū)域的地表形變信息，被廣泛應(yīng)用于城市地面沉降監(jiān)測(cè)[1-2]。由于礦區(qū)多被農(nóng)作物和林木等植被覆蓋，導(dǎo)致DInSAR(Differential InSAR)技術(shù)的相干性弱，使得監(jiān)測(cè)結(jié)果可靠性低[3]。時(shí)序合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)能夠克服常規(guī)DInSAR技術(shù)失相干的缺陷[4]。為了提高時(shí)序InSAR在低相干區(qū)的監(jiān)測(cè)能力，Andy Hopper等人提出了StaMPS/MTI(Stanford Method for Persistent Scatterers/Multi-Temporal InSAR)方法[5-6]，采用三維時(shí)空解纏算法獲取目標(biāo)的時(shí)序形變信息。本文采用StaMPS/MTI方法對(duì)覆蓋新汶礦區(qū)的48景Sentinel1-A影像進(jìn)行時(shí)序InSAR數(shù)據(jù)處理，獲取2016～2020年間的新汶礦區(qū)地表形變時(shí)空分布信息，對(duì)近4年的地面沉降情況進(jìn)行分析，揭示“后煤炭時(shí)代”老礦采煤沉陷區(qū)的形變規(guī)律及現(xiàn)狀。

1 研究區(qū)域概況及數(shù)據(jù)

1.1 研究區(qū)域概況

新汶礦區(qū)地處山東省新泰市，東起澇坡、員外哨，西至南楊家莊、大霞霧，東西均以煤系地層露頭為界，北至蓮花山斷裂，南至蒙山背斜北坡奧陶系地層出露處，東西長(zhǎng)約30 km，南北寬約8～10 km，面積約305 km2，地層呈北西走向。地理位置如圖1所示，礦區(qū)主要包括山東能源新礦集團(tuán)國有煤礦和新泰市地方煤礦，是華東地區(qū)重要的能源基地之一。

圖1 新汶礦區(qū)地理位置圖

1.2 數(shù)據(jù)

SAR數(shù)據(jù)選用2016年1月～2020年1月共48期Sentinel-1A升軌影像，如表1所示。在數(shù)據(jù)處理過程中，使用歐空局發(fā)布的精密軌道數(shù)據(jù)和SRTM的90 m分辨率的外部DEM。

表1 Sentinel-1A影像軌道和日期

2 SBAS-StaMPS方法

StaMPS是一種MT-InSAR算法[4]，為近年來解決長(zhǎng)時(shí)序、小位移監(jiān)測(cè)的有效方法，主要用于探測(cè)非城市地表區(qū)域的形變信息，采用相位時(shí)序分析算法，根據(jù)時(shí)序相干系數(shù)判斷相位的時(shí)域穩(wěn)定性[7]。根據(jù)干涉相位中雷達(dá)信號(hào)和噪聲的不同時(shí)空屬性特征，通過時(shí)空濾波方法加以分離，通過三維相位解纏來獲取地表形變信息[8]。

SBAS-StaMPS方法是在給定的多幅SAR影像中，利用相干測(cè)度定義像元信噪比，由滿足信噪比閾值的永久散射體候選點(diǎn)構(gòu)成StaMPS-MTI技術(shù)的高相干點(diǎn)集，按照SBAS(Small Baseline Subsets)法形成時(shí)序差分干涉影像序列，從而提取出時(shí)序形變信息和年平均形變速率。該方法一方面提供了在高空間分辨率下操作的可能性，并能夠利用干涉相位的空間相關(guān)性識(shí)別慢去相關(guān)濾波相位(SDFP，Slowly-decorrelating Filtered Phase)像元[9]；另一方面，較好地揭示了多主時(shí)間序列的三維小基線干涉圖[10]。為了最大限度地提高干涉圖的相干性，SBAS-StaMPS通過考慮空間基線、時(shí)間間隔和圖像對(duì)的多普勒頻率質(zhì)心差生成小基線干涉圖，在保證圖像對(duì)的網(wǎng)絡(luò)中不存在孤立簇的情況下，選擇空間、時(shí)間和多普勒基線低于給定閾值的圖像對(duì)形成干涉圖。再用頻譜濾波器對(duì)干涉圖進(jìn)行濾波，消除方位向上的非重疊多普勒譜，減少距離向上的幾何去相關(guān)[8]。數(shù)據(jù)處理流程如圖 2所示，基本步驟包括干涉圖生成、相位穩(wěn)定性估計(jì)、PS點(diǎn)選擇和時(shí)序位移估計(jì)。采用SAR影像生成差分干涉圖，利用外部DEM去除地形相位；基于振幅濾波方法，采用振幅閾值和振幅離差指數(shù)選出永久散射體候選點(diǎn)集(PSCs)；對(duì)PSCs進(jìn)行相位穩(wěn)定性分析，基于時(shí)間相干系數(shù)閾值篩選PS點(diǎn)；運(yùn)用3D相位解纏，通過時(shí)空濾波提取PS點(diǎn)的相位時(shí)間序列，獲取視線向形變信息。

圖2 StaMPS數(shù)據(jù)處理流程

SBAS-StaMPS方法在短時(shí)間基線中進(jìn)行空間上高相干性的濾波，能夠克服PSInSAR方法的長(zhǎng)基線低相干問題，對(duì)于植被覆蓋較多而散射體密度較低的礦區(qū)，通過PS點(diǎn)和SDFP點(diǎn)的共同選取，能夠獲取更多的PS點(diǎn)，且提高了點(diǎn)位的可信度。

3 數(shù)據(jù)處理結(jié)果

根據(jù)三基線和最小原則[11]選取2017年12月10日的影像作為主影像，其余影像與其配準(zhǔn)，形成201個(gè)SBAS干涉對(duì)，最大時(shí)間基線為192d，最大垂直基線為99.27 m，得到新泰市2016～2020年間地表形變的雷達(dá)視線向(LOS)年平均速率圖，如圖3所示，正值(藍(lán)色)表示上升，負(fù)值(紅色)表示下沉。選取覆蓋新汶礦區(qū)的區(qū)域作為研究區(qū)，如圖4所示，提取沉降明顯的A～H共8處區(qū)域進(jìn)行分析，其形變時(shí)間序列如圖5所示，藍(lán)色表示相對(duì)于主影像的形變量，紅色直線表示年平均形變速率。

圖3 新泰市地表形變年平均速率圖(2016～2020)

圖4 新汶礦區(qū)地表形變年平均速率圖(2016～2020)

圖5 礦區(qū)沉降時(shí)間序列

4 實(shí)驗(yàn)結(jié)果分析

(1)在沉降區(qū)域整體空間分布上，從圖 3可以看出，整個(gè)新泰市轄區(qū)的主要地面沉降區(qū)域位于中東部的新汶地區(qū)，該區(qū)域是山東能源新汶礦業(yè)集團(tuán)和新泰市地方煤礦的主要采煤區(qū)。圖 4顯示了新汶礦區(qū)A～H的8個(gè)明顯沉降區(qū)域，均分布在柴汶河以北，經(jīng)過實(shí)地調(diào)繪及同相關(guān)煤礦生產(chǎn)企業(yè)調(diào)研，沉降區(qū)域所屬煤礦如表2所示，除兩處地方煤礦外，其余6處均為山東能源新汶礦業(yè)集團(tuán)的權(quán)屬企業(yè)采煤區(qū)。說明SBAS-StaMPS方法提取的形變區(qū)域與實(shí)際相符，能夠反映礦區(qū)沉降區(qū)域。

表2 沉降區(qū)域關(guān)系對(duì)應(yīng)表

(2)在形變速率上，提取各煤礦沉降中心區(qū)域按照沉降速率分類，如表3所示。年平均速率超過40 mm/a的有1處，為孫村煤礦，年平均速率在35～40 mm/a的有1處煤礦，年平均速率在30～35 mm/a的有3處煤礦，其余小于30 mm/a。

表3 沉降區(qū)分類表

(3)從沉降量上看，自2016年以來，各煤礦累計(jì)沉降量及年平均速率，如表 4所示。其中，4處礦區(qū)沉降量超過150 mm，分別是協(xié)莊煤礦、翟鎮(zhèn)煤礦、孫村煤礦和碗窯頭地方礦；其余煤礦沉降量都在100～150 mm，良莊礦業(yè)沉降量最小。

表4 沉降量及平均沉降速率

(4)從沉降趨勢(shì)上看(如圖6所示)，A、C、F沉降速率平穩(wěn)，為線性勻速沉降趨勢(shì)。B、E、H在2019年后有明顯的沉降加速趨勢(shì)，沉降速率分別為-70.86 mm/a、-74.84 mm/a和-52.12 mm/a，表明存在沉降速度增大現(xiàn)象，不均勻沉降具有災(zāi)害隱患，需要加強(qiáng)監(jiān)測(cè)，關(guān)注沉降趨勢(shì)。D、G表現(xiàn)為階梯型趨勢(shì)，特別地，D在2019年后沉降量趨于平穩(wěn)，沉降趨勢(shì)明顯減緩，該區(qū)域?qū)儆谑⑷V業(yè)采煤區(qū)，由于該企業(yè)資源枯竭，按照國家退出產(chǎn)能政策于2018年11月后開始實(shí)施關(guān)井閉坑，停止開采，故沉降量減緩，監(jiān)測(cè)結(jié)果與實(shí)際情況吻合。

圖6 典型沉降區(qū)的沉降趨勢(shì)圖

5 結(jié) 論

SBAS-StaMPS方法以其獨(dú)特的PS和SDFP選點(diǎn)算法，有效克服了D-InSAR技術(shù)的時(shí)空和幾何失相干問題，對(duì)林地和農(nóng)田覆蓋較多而散射體密度較低的礦區(qū)能夠提取可靠的PS點(diǎn)，且以mm級(jí)的精度獲取年沉降速率。本文采用該方法對(duì)2016～2020年間的C波段Sentinel-1A數(shù)據(jù)進(jìn)行時(shí)序分析，有效識(shí)別出新汶礦區(qū)8處采煤沉陷區(qū)，提取了4年的大范圍連續(xù)沉降結(jié)果，累計(jì)沉降量均在100 mm以上，最大平均沉降速率為44.49 mm/a，揭示了各采煤區(qū)地面沉降的時(shí)空分布和形變趨勢(shì)。研究結(jié)果表明SBAS-StaMPS是一種礦區(qū)長(zhǎng)時(shí)序地表形變監(jiān)測(cè)的有效技術(shù)，能夠?yàn)楸O(jiān)測(cè)礦區(qū)地表形變，指導(dǎo)煤礦合理開采和可持續(xù)發(fā)展提供一定的技術(shù)支持和理論借鑒。