甄洪帥，張紅月，姬宗皓，殷幼松

(1.自然資源部采煤沉陷區(qū)綜合治理工程技術(shù)創(chuàng)新中心，山東 濟(jì)寧 272000；2.山東省魯南地質(zhì)工程勘察院，山東 兗州 272100；3.濟(jì)寧市采煤沉陷地治理中心，山東 濟(jì)寧 272000；4.環(huán)球星云遙感科技有限公司，浙江 湖州 313200)

0 引言

霄云礦區(qū)位于濟(jì)寧市金鄉(xiāng)縣霄云鎮(zhèn)，于2008年4月開工建設(shè)，區(qū)內(nèi)地勢平坦，地面標(biāo)高為35.8～44.9m，總體趨勢為西高東低。土地肥沃，地面村莊較多。區(qū)內(nèi)地層為走向近EW、傾向北的單斜構(gòu)造，井田北至劉吳莊斷層，西為曹馬集斷層，東距魚臺縣城約20km，南起各煤層露頭。該礦區(qū)自2011年底到2015年中由水準(zhǔn)測量觀測到的最大沉降量超過300mm。為了獲取霄云礦區(qū)地表沉降的現(xiàn)狀資料以服務(wù)于后期礦區(qū)治理方案的規(guī)劃設(shè)計(jì)，研究該礦區(qū)采煤沉陷的空間分布及時序演化特征具有重要意義。然而，傳統(tǒng)觀測手段雖然具有精度高、監(jiān)測量級大等特點(diǎn)，但在礦區(qū)面域形變監(jiān)測應(yīng)用中存在著經(jīng)濟(jì)代價高、觀測周期長、采樣點(diǎn)密度難以反映礦區(qū)整體形變特征等突出問題[1-7]。Lidar技術(shù)可以通過高程變化從而提取礦區(qū)形變特征，但費(fèi)用昂貴，在大規(guī)模工程化應(yīng)用過程中存在一定局限性[8-10]。利用合成孔徑雷達(dá)干涉測量(InSAR)技術(shù)方法進(jìn)行采煤沉陷區(qū)的沉降監(jiān)測在德國、澳大利亞和波蘭等國研究較早[11-13]，近年我國相關(guān)學(xué)者應(yīng)用該技術(shù)也對礦區(qū)沉陷開采做了相應(yīng)的研究及應(yīng)用工作[14-17]。而InSAR技術(shù)監(jiān)測礦區(qū)沉陷的應(yīng)用效果取決于諸多因素，主要包括：SAR影像空間分辨率、影像數(shù)量及時間采樣頻率、雷達(dá)波長、礦區(qū)工作面地物特征及其開采特點(diǎn)等。就空間分辨率這一要素而言，中等分辨率SAR影像在解譯過程中容易造成相位混疊、相干點(diǎn)目標(biāo)密度較小等問題，對整個礦區(qū)形變場空間形態(tài)及量級程度的刻畫存在不足，繼而影響礦區(qū)沉陷范圍線的提取精度[18-22]。鑒此，該文以加拿大RADARSAT-2超精細(xì)寬(ExtraFine，XF)模式5m高分辨率時序影像為數(shù)據(jù)源，利用時序InSAR技術(shù)對濟(jì)寧市霄云礦區(qū)開展沉陷調(diào)查監(jiān)測，為后期礦區(qū)綜合治理及生態(tài)修復(fù)工作的順利開展提供切實(shí)可靠的現(xiàn)狀資料。

1 SAR影像選擇

影像源的選擇上主要根據(jù)當(dāng)前主流SAR衛(wèi)星存檔數(shù)據(jù)可利用情況，并結(jié)合采煤沉陷區(qū)地面沉降監(jiān)測對SAR數(shù)據(jù)源幅寬大、分辨率高、重訪周期穩(wěn)定及數(shù)據(jù)獲取可靠等方面的參數(shù)需求，最終確定使用RADARSAT-2的超精細(xì)寬高分辨率SAR存檔數(shù)據(jù)(ExtraFine，XF)，實(shí)現(xiàn)工作區(qū)全覆蓋、高精度、分年度的監(jiān)測，確保監(jiān)測區(qū)因時間差異而形成的塌陷形變可以在InSAR結(jié)果中有很好的體現(xiàn)。同時，為了平衡影像成本和監(jiān)測效果，選用每年度時間上分布較為均勻的5期影像，2016--2018年共30景XF時間序列數(shù)據(jù)進(jìn)行觀測，SAR數(shù)據(jù)集基本參數(shù)見表1。

表1 SAR數(shù)據(jù)集基本參數(shù)

2 時序InSAR分析技術(shù)

2.1 差分相位模型

時序InSAR分析是D-InSAR監(jiān)測技術(shù)的應(yīng)用拓展，通過對不同時刻獲取的同一區(qū)域SAR影像上受失相干現(xiàn)象影響較小的像元(高相干點(diǎn)目標(biāo))開展相位穩(wěn)定性分析，以實(shí)現(xiàn)長時間跨度地表沉降信息的高精度獲取。在經(jīng)過時間序列SAR影像同名點(diǎn)位的配準(zhǔn)、從雷達(dá)坐標(biāo)系向地理坐標(biāo)系的編碼轉(zhuǎn)換、不同時刻影像多基線組合選取及利用外部DEM對干涉相位去除地形影響后，時序差分相位模型可表示為：

φx,i=φdef,x,i+Δφε,x,i+φα,x,i+Δφorb,x,i+φn,x,i

式中：φx,i代表第i幅差分干涉對上第x個像元的差分干涉相位；φdef,x,i為雷達(dá)視向上的形變相位；Δφε,x,i為外部參考DEM不準(zhǔn)確引入的殘余地形相位誤差；φα,x,i代表大氣延遲誤差；Δφorb,x,i為衛(wèi)星軌道數(shù)據(jù)不準(zhǔn)確引入的相位誤差；φn,x,i為其他噪聲相位，如：熱噪聲、體散射變化誤差及數(shù)據(jù)處理中引入的誤差等。差分相位模型的構(gòu)建是各種時序InSAR分析方法的解算基礎(chǔ)。

2.2 監(jiān)測方法

當(dāng)前，根據(jù)研究對象的不同，從差分干涉相位中獲取形變成分主要有2種處理思路：一種是假定形變區(qū)符合某種時序演化規(guī)律，先解算出主體形變成分，然后從殘余相位中分離出剩余部分，兩者之和便是形變信息；另一種則是通過從差分干涉相位中逐步剝離其他干擾相位成分進(jìn)而獲取總的形變信息。與城市緩慢地表沉降可采用前一種時間線性形變模型描述不同，由煤炭開采引起的形變發(fā)育過程較為復(fù)雜、沉陷量級大，且各工作面開采時段并不一致，故應(yīng)該采用第二種處理思路，即：不假定時間形變模型的空間相關(guān)性InSAR分析方法[23-24]。其處理步驟包括：時序SAR影像配準(zhǔn)、干涉圖生成、地形相位去除、初選相干點(diǎn)目標(biāo)、空間相關(guān)和非相關(guān)視角誤差估算、終選相干點(diǎn)目標(biāo)、相位解纏和大氣相位獲取等，形變分析過程如下。

假設(shè)φdef,φα.Δφorb在一定空間尺度上具有相關(guān)性，Δφε.φn在此距離上非空間相關(guān)，且在x以為中心、L為半徑的圓內(nèi)其均值為零。以第x像素半徑L內(nèi)所有高相干點(diǎn)目標(biāo)的相位平均值可表示為：

式中：

由DEM誤差造成的相位差與垂直基線成正比例：

φε，x,i=B⊥,x,i·Kε,x

式中：Kε,x為比例常量，改寫為：

Kε,x是唯一與基線相關(guān)的因子，可通過最小二乘法估算。相位穩(wěn)定性評定因子表述如下：

式中：φ＇ε,x,i為由Kε,x估計(jì)不準(zhǔn)確而引起的殘余誤差，由于相鄰相干點(diǎn)目標(biāo)間空間相關(guān)DEM誤差相差并不大，故該項(xiàng)對解纏誤差的影響較小。對上式進(jìn)行三維迭代相位解纏處理，并利用大氣、軌道和形變信號三者之間時空頻譜特性的差異，實(shí)現(xiàn)形變信號和噪聲信號的有效分離。

3 礦區(qū)開采沉陷提取

3.1 干涉組合選取

對于時序InSAR分析而言，可通過兩種組合模式構(gòu)建干涉對，一種是單一主影像的放射狀網(wǎng)絡(luò)；另一種則為短時間、空間基線的小基線集組合。前一種模式一般在影像數(shù)量大于25期以上時采用；而后一種模式能夠以更高的空間密度提取到相干點(diǎn)目標(biāo)，對礦區(qū)空間形變特征的反應(yīng)更加明顯、直觀。該文受限于數(shù)據(jù)量少且影像間隔時間較長，采用小基線集干涉組合方式，將空間垂直基線300m、時間基線200d作為閾值，所選組合情況如表2所示。

表2 研究區(qū)圖幅干涉組合時空基線

3.2 差分干涉處理

RADARSAT-2XF模式影像單視情況下方位向分辨率為2.5m左右，而距離向則為5m。為抑制相干斑噪聲且同時保持高分辨影像的精細(xì)化監(jiān)測優(yōu)勢，采用距離向和方位向1∶2的窗口對配準(zhǔn)后影像進(jìn)行多視處理，在提高相位信噪比的同時，使得方位向和距離向分辨率接近一致，改善影像判讀效果。此外，由于煤炭開采區(qū)多位于農(nóng)田等植被覆蓋度高、相干性低的區(qū)域，需要對去除地形貢獻(xiàn)后的相位進(jìn)行濾波。

3.3 研究區(qū)沉陷監(jiān)測結(jié)果

在對影像解譯處理過程中需要注意兩個方面：一方面是大氣條件和解算傳播誤差的影響。大氣一般在空間1～3km呈現(xiàn)出較強(qiáng)的相關(guān)性，且隨著距離的增大，相關(guān)性逐漸減弱；而解算傳播誤差同樣和距離有關(guān)，空間范圍越大，誤差積累越明顯。鑒此，為有效控制大區(qū)域時序InSAR分析過程中大氣條件及解算傳播誤差的影響，以霄云礦區(qū)邊界外擴(kuò)1～3km為界，對獲取的時序差分干涉相位進(jìn)行解算分析。另外一方面則是需要排除地下水開采這類由居民生產(chǎn)活動產(chǎn)生的地面沉降現(xiàn)象。與礦采形變不同，地下水開采不會在短期內(nèi)發(fā)生劇烈的地表形變，因此其相位條紋不會過于明顯；此外該類現(xiàn)象多發(fā)生在人口聚集地，可通過形變場與村莊空間分布的疊加分析予以判別。綜合上述考慮后，得到2016—2018年研究區(qū)范圍內(nèi)的沉陷演化過程，如圖1和2所示。

圖1 霄云礦區(qū)時序形變演化過程

圖2 霄云礦區(qū)累計(jì)沉降量及其等值線

4 結(jié)果分析

時序InSAR技術(shù)獲取的監(jiān)測結(jié)果既可以從空間上分析礦區(qū)開采引起的沉陷對于周邊環(huán)境的影響，尤其在“三下壓煤”開采情況下地面沉陷對于建、構(gòu)筑物的安全穩(wěn)定構(gòu)成的威脅，同時也可以與開采面的工作時段相結(jié)合，從而綜合礦區(qū)開采工藝、地質(zhì)條件等因素對沉陷發(fā)育規(guī)律進(jìn)行研究，進(jìn)一步加深礦區(qū)沉陷過程的認(rèn)知，為開采沉陷控制方案的制定提供輔助資料。

4.1 研究區(qū)地質(zhì)采礦條件

霄云煤礦采用立井開拓，井底水平為-790m，采煤方法為綜采放頂煤采煤法，開采標(biāo)高為-430m～-1500m[25]。其中，首采區(qū)地面影響范圍內(nèi)有共有11個村莊，分布較為零散；二采區(qū)范圍內(nèi)有9個村莊。該區(qū)為全隱蔽式華北型石炭二疊紀(jì)含煤地層，地層自上而下分別為第四系、古近系、二疊系、石炭系、奧陶系。區(qū)內(nèi)總體為向北傾伏的單斜構(gòu)造，地層傾角10°～21°，淺部陡深部緩；斷裂構(gòu)造較發(fā)育，共發(fā)現(xiàn)斷層19條，其中落差大于100m的2條；50～100m的2條；20～50m的11條，斷層走向多為NE向和EW向，且主要斷層控制程度均已達(dá)到可靠或較可靠。礦區(qū)煤系地層平均總厚度260m，含煤19層。可采、局部可采煤層3層(3、12下、16上)，平均總厚度5.37m。3煤層厚2.25～7.69m，平均3.98m，穩(wěn)定全區(qū)可采；12下煤層厚0.34～1.28m，平均0.71m，局部可采；16上煤層厚0.60～0.74m，平均0.68m，較穩(wěn)定局部可采。

4.2 空間影響分析

霄云礦區(qū)井田存在較厚的沖積層，在進(jìn)行工作面開采時，地表下沉量相比于基巖情況下要大，由于厚沖積層的移動傳遞，會使得淺層土體固結(jié)壓縮時形變較為均勻，進(jìn)而反映到地表時會相對平緩，不會引起較大量級的集中形變；此外，開采煤層埋藏較深，礦區(qū)內(nèi)煤層頂板以砂巖為主且較為穩(wěn)定，石炭-二疊系硬巖超過150m，條帶開采條件下上覆巖層與硬巖巖梁將整體緩慢下沉，使得開采引起的地表下沉、變形值會比較均衡[26]，這些從監(jiān)測結(jié)果中得到了很好的體現(xiàn)。霄云煤礦主采3煤層，采深約-790m，采用條帶式或隔面開采，因此，開采沉陷并不存在非常集中的大量級變形。從圖3中可以看出，在監(jiān)測時段內(nèi)，地面沉陷主要集中在2015—2018年度的開采工作面周邊，而2012—2014年的開采面已經(jīng)基本上達(dá)到了穩(wěn)沉狀態(tài)，不再有明顯的下沉趨勢。

4.3 開采沉陷規(guī)律分析

霄云礦區(qū)共劃分為6個采區(qū)，2012—2018年的開采工作面位置分布如圖3中藍(lán)色框體。對礦區(qū)內(nèi)不同開采時段工作面選擇點(diǎn)位繪制沉陷形變曲線以探究地面沉陷發(fā)育規(guī)律，所選取的6處具有代表性點(diǎn)位空間位置如圖3中A—F所示，形變序列演化曲線見圖4。

圖3 礦區(qū)開采面空間沉降影響范圍

經(jīng)分析，圖4中的6處形變曲線主要體現(xiàn)出：①InSAR形變演化過程與開采時段引起的沉陷吻合度高。圖4a位置沉陷由2015年工作面的開采所導(dǎo)致，從曲線發(fā)育趨勢可以看出，存在明顯的收斂趨勢，而結(jié)合空間沉陷范圍可以看出，大部分在2012—2014年間的開采面周邊已經(jīng)達(dá)到穩(wěn)沉狀態(tài)，不存在突出形變跡象；而在圖4b中可以看到，在2017年工作面尚未處于開采狀態(tài)時，地面點(diǎn)位周邊較為穩(wěn)定，而隨著開采工作的進(jìn)行，隨之發(fā)生快速下沉過程，這些點(diǎn)位的沉降跡象都與礦區(qū)沉陷發(fā)育規(guī)律較為一致。②距離開采工作面的遠(yuǎn)近不同，地表沉陷曲線演化形態(tài)非常相似，這一點(diǎn)如圖4b和圖4c所示。兩處點(diǎn)位雖然在時序沉降數(shù)值以及曲線的緩和程度上存在一定差異，但不難看出，在不受鄰近開采工作面影響情況下，兩者曲線形態(tài)保持一致；同時，距離開采工作面越近，曲線的變化程度越為劇烈，沉陷量值越大。③多個鄰近開采工作面共同作用對于沉陷曲線的形態(tài)存在顯著影響。為了便于對比分析不同工作面、不同開采時間對開采沉陷的影響，將D，E和F點(diǎn)繪制于一張圖上，如圖4d所示。從圖4d中可以看出，D點(diǎn)2016年開采引起的形變?nèi)匀惶幱谥泻笃诎l(fā)育過程中。雖然E和F點(diǎn)位都主要受2018年開采工作面影響，但可以明顯看出這兩處點(diǎn)位變形過程并非遵循先穩(wěn)定而后快速沉陷的演化規(guī)律，這主要是因?yàn)檫@兩處點(diǎn)位在2018年工作面進(jìn)行開采之前，周邊還存在2016年的開采面，由于2018年開采時該處尚未達(dá)到穩(wěn)沉狀態(tài)，因此出現(xiàn)了2016年中期到2017年初的緩慢沉陷過程，而到了2017年基本處于沉陷發(fā)育的末端。在2018年隨著新的工作面進(jìn)行開采，導(dǎo)致地面沉陷現(xiàn)象又重復(fù)發(fā)生，便出現(xiàn)了上述的形變發(fā)育混疊過程。

5 結(jié)論

該文利用15期RADARSAT-2XF模式高分辨率數(shù)據(jù)，以濟(jì)寧市霄云礦區(qū)為例，采用基于空間相關(guān)性分析的時間序列InSAR方法對其采空區(qū)沉陷現(xiàn)狀進(jìn)行了解譯獲取，結(jié)果表明：

(1)高分辨率SAR影像與時序InSAR技術(shù)的結(jié)合可以實(shí)現(xiàn)高精度的礦區(qū)開采形變信息的準(zhǔn)確獲取，具有監(jiān)測效率高，經(jīng)濟(jì)代價小、可延續(xù)性好等突出優(yōu)勢，在大范圍采煤沉陷區(qū)沉陷監(jiān)測方面極具應(yīng)用潛力。

(2)高分辨率時序InSAR分析方法可以很好捕捉到開采沉陷區(qū)的時空形變特征，通過對時序相位成分的分析，在獲取沉陷區(qū)累計(jì)沉降量的同時，還可以得到每景影像獲取時刻的變形信息，為礦區(qū)地面沉陷演化過程的分析提供了資料。

(3)受當(dāng)前雷達(dá)衛(wèi)星運(yùn)轉(zhuǎn)條件以及傳感器自身成像要求的限制，InSAR技術(shù)在礦區(qū)沉陷監(jiān)測中的應(yīng)用效果受到了一定程度的制約，但可以預(yù)見，隨著SAR衛(wèi)星及傳感器的迭代更新及InSAR技術(shù)方法的不斷發(fā)展，其在礦區(qū)沉陷監(jiān)測方面將會發(fā)揮更加突出的作用。