高分辨率時序InSAR在礦區開采沉陷監測中的應用
——以霄云礦區為例

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

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

0 引言

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

1 SAR影像選擇

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

表1 SAR數據集基本參數

2 時序InSAR分析技術

2.1 差分相位模型

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

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

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

2.2 監測方法

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

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

式中：

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

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

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

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

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

3 礦區開采沉陷提取

3.1 干涉組合選取

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

表2 研究區圖幅干涉組合時空基線

3.2 差分干涉處理

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

3.3 研究區沉陷監測結果

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

圖1 霄云礦區時序形變演化過程

圖2 霄云礦區累計沉降量及其等值線

4 結果分析

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

4.1 研究區地質采礦條件

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

4.2 空間影響分析

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

4.3 開采沉陷規律分析

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

圖3 礦區開采面空間沉降影響范圍

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

5 結論

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

(1)高分辨率SAR影像與時序InSAR技術的結合可以實現高精度的礦區開采形變信息的準確獲取，具有監測效率高，經濟代價小、可延續性好等突出優勢，在大范圍采煤沉陷區沉陷監測方面極具應用潛力。

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

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