基于SBAS-InSAR技術的武安塌陷區形變監測分析

邸 明 慧，張 麗 云，賀 月 娟，徐 寧，喬 良*，李 曉 婧

(1.河北省科學院地理科學研究所/河北省地理信息開發應用工程技術研究中心，河北 石家莊 050011；2.河北交通職業技術學院，河北 石家莊 050035)

0 引言

地下礦產資源開采易導致地面塌陷，不僅嚴重破壞土地(特別是耕地)和地面建筑物，還會影響交通運輸和水資源等[1],因此，礦區地表形變監測對于礦區安全生產、地質環境治理和災害防控等具有重要作用[2]。InSAR(Interferometric Synthetic Apeture Radar)技術具有不受氣候制約和大范圍、高精度的監測優勢，已被廣泛應用于各種地質災害的監測研究中[3，4]，尤其是SBAS-InSAR(Small Basedline Subsets InSAR)技術克服了傳統測量方法存在的時(空)間失相干和大氣效應限制，適用于長時序緩慢非線性形變監測[5]，并且在采空塌陷區應用較好。

太行山地區是河北省重要的礦產資源分布區之一，其北段為Fe-Mo、Cu-Pb、Zn-Au成礦區，南段為Fe-煤-鋁土礦-石膏成礦區,礦產資源的大量開采導致該區出現嚴重的地質災害和生態環境問題。武安地區是河北太行山南段礦產資源集中分布區，因采礦活動造成多處不同程度的地面塌陷[6]，尤其在中部丘陵和盆地區域地面塌陷更普遍[7]。張振生等[8]基于D-InSAR技術發現1993-2004年武安礦山存在11處以煤礦為主的塌陷區，沉降范圍由小變大再趨于平穩；張景發等[6，9]利用InSAR技術提取武安市慧蘭村礦區塌陷區，發現塌陷區范圍和沉降量均隨開挖時間推移而變化；周洪月等[10，11]基于時序InSAR處理發現武安市存在一個沉降中心。但以上研究均未針對武安市域內的塌陷區開展區域性監測和分析。1996-2014年武安市大量小礦被關停整合，2009年啟動沉陷區綜合治理工程，2017年推出礦區生態修復政策，許多礦區得到整治。在此背景下，停采和修復礦區的地面塌陷是否存在繼續擴大以及是否有新的塌陷區形成，均直接關系到地方政府相關政策的制定，亟須從全區出發開展地面塌陷區的形變監測研究。因此，本文以武安市為研究區，利用2016年10月-2017年3月和2017年10月-2018年3月兩時段各15景降軌Sentinel-1A衛星影像，通過SBAS-InSAR技術反演得到研究區雷達視線方向(Line Of Sight，LOS)的時序累積形變量，結合塌陷區和地裂縫災害點數據對反演結果進行驗證，分析塌陷區的空間分布和形變特征，為地方政府防災減災提供科學依據。

1 數據與方法

1.1 研究區與數據

武安市位于河北省邯鄲市西北(113°45′～114°22′E，36°28′～37°1′N)，處于太行山隆起向華北平原沉降帶的過渡區，地勢西北高(最高約為1 900 m)、東南低，連續下降地形被紫山—鼓山打斷，形成寬緩的武安盆地(圖1)。武安市是全國58個重點產煤縣(市)和全國四大富鐵礦基地之一，“邯邢式鐵礦”因鐵礦資源在邯鄲—邢臺一帶集中分布而得名。邯邢式鐵礦集中區可以劃分為“三帶兩盆”，其中“三帶”自西向東依次為符山成礦帶、固鎮—武安—礦山—綦村成礦帶和鼓山—洪山—新城成礦帶，“兩盆”指“三帶”間的陽邑盆地和武安盆地[12]。在《河北省地質災害防治“十三五”規劃》中，武安—沙河鐵、煤礦區屬于地面塌陷高易發區。本文將研究區劃分為M1、M2和M3共3個區域(圖1)，以便進行更細尺度的分析。

圖1 研究區概況Fig.1 Survey of the study area

為減少地表散射機制改變對差分干涉的影響[13]，本研究選取2016年10月-2017年3月和2017年10月-2018年3月各15景降軌Sentinel-1A衛星影像(https://search.asf.alaska.edu/)作為數據源，以30 m分辨率的SRTM1(http://www.gscloud.cn/)地形數據作為外部參考DEM；利用武安市已知災害點(表1)對反演結果進行驗證，災害點數據來源于《邯鄲市2016-2019年地質災害防治方案的通知》[14]，包括14處塌陷區和3處地裂縫，其中中型規模11處，小型規模6處，均處于欠穩定狀態。

表1 研究區塌陷區與地裂縫災害點和反演結果驗證Table 1 Subsidence areas and ground fissures in the study area and verification of inversion results

1.2 SBAS-InSAR技術

相比D-InSAR技術，SBAS-InSAR技術可獲取時間序列的干涉結果并削弱大氣的影響；相比PS-InSAR技術，其數據量要求少，從而在一定程度上可實現同一季節內地物反射條件基本不變，保證像對干涉質量。該技術實現流程為：假設共有N+1幅SAR影像，根據一定的時(空)間基線選取原則，選擇影像組合，生成M((N+1)/2≤M≤N(N+1)/2)景短基線距的干涉圖；基于數據基線長度和研究區實際情況，對像對進行短基線干涉處理，再利用外部DEM模擬地形相位去除地形相位的干擾；依據干涉圖、相位離散度和相位標準差，選擇高相干的點進行相位解纏和殘余地形相位估計，從而建立分時段的平均變形速率、高程誤差和差分相位的模型方程組；最后利用奇異值分解法(SVD)計算最小二乘解，估計時序非線性形變和DEM誤差，利用高斯濾波器進行空間濾波和時間濾波以去除大氣相位。

2 結果分析

2.1 整體空間分布特征

本研究應用ENVI軟件的SARscape模塊，對影像依次進行裁剪、配準、干涉(最大空間基線閾值為5%、最大時間基線閾值為120 d)、濾波、解纏、軌道精煉、形變估計和去大氣效應處理，生成兩時段的時序形變分布圖，結合Google Earth影像去除非塌陷區，得到研究區主要塌陷區LOS向累積形變量空間分布(圖2)。結果顯示，2016年10月-2017年3月塌陷區最大累積形變量為197.84 mm，平均累積形變量為16.44 mm；2017年10月-2018年3月塌陷區最大累積形變量為172.16 mm，平均累積形變量為14.22 mm，均小于前者，說明研究區塌陷有向好發展趨勢。塌陷區的總體空間分布趨勢基本一致，主要分布在研究區東部，在西部地區分布比較零散且規模較小，但塌陷區主要形變中心位置和規模都隨時間有所變化。

圖2 研究區LOS向累積形變量空間分布Fig.2 Spatial distribution of the line of sight accumulative deformation in the study area

利用ArcGIS將17處災害點(表1)與累積形變量空間分布圖疊加(圖3)對比，發現12處災害點分布在反演結果塌陷區內，聚隆煤礦周圍地裂縫、豐里村及周圍地裂縫和西萬年村塌陷區位于反演塌陷區的影響范圍內，河西村西山坡塌陷區和龍洞山塌陷區兩處未顯示形變。由此可見，本文反演結果較準確地反映了區內災害點的分布。

圖3 研究區塌陷區和地裂縫分布Fig.3 Distribution of subsidence areas and ground fissures in the study area

結合研究區主要鐵、煤礦分布帶可以看出，本研究反演出的塌陷區主要位于鐵礦分布帶(固鎮—武安—礦山—綦村成礦帶)和煤礦分布帶(鼓山—洪山—新城成礦帶)內，與武安地區鐵、煤礦資源分布范圍高度一致，這也印證了本文反演結果的準確性。固鎮—武安—礦山—綦村成礦帶下伏基巖地層傾向東南，地表高度從西往東雖然逐漸變小，但盆地內的松散沉積物和巖層厚度自西向東遞增[15]，這可能導致西側礦產資源分布深度比東側淺，發現和開采時間更早。因此該成礦帶西側塌陷區規模比東側小，推測西側礦區多為老礦區或棄采礦區，開采量小或礦區恢復使塌陷區累積形變量逐漸變小；而成礦帶東側為新礦區，正處于開采階段，塌陷區累積形變量相對較大。

2.2 局部空間特征

由圖4和表2可知：在M1分區中，塌陷區1、3、4、9在兩時段形變速率變化不大，形變中心轉移方向均不相同，塌陷區1向北，塌陷區3向西南，塌陷區4、9向東；塌陷區2、6、8在后一時段形變速率減小，其中塌陷區2和8的規模和累積形變量顯著減小，平均形變速率分別由7.62 mm/d、8.61 mm/d減小為2.57 mm/d、1.34 mm/d，塌陷區6規模變化不大，但累積形變量明顯變小；塌陷區5、7在后一時段規模變化不大，但形變速率增大，增幅分別為5.67 mm/d和3.19 mm/d，形變中心發生轉移。

表2 兩時段不同分區塌陷區平均形變速率Table 2 Average deformation rates of different subsidence areas in three sub-regions in two periods mm/d

圖4 M1、M2和M3分區LOS向累積形變量Fig.4 Maps of the line of sight accumulative deformation in M1，M2 and M3 sub-regions

在M2分區中，塌陷區6、9在兩個時段平均形變速率相近，規模變化不大，但形變中心均向南轉移；塌陷區1、4、5、8在后一時段平均形變速率減小，表明開采強度有所減弱，其中塌陷區4形變速率最小，形變中心明顯消失，塌陷區1和8范圍有所增大，塌陷區5形變速率明顯減小；塌陷區2、3、7在后一時段平均形變速率增大，其中塌陷區2范圍增大，塌陷區3范圍減小且形變中心南移，塌陷區7規模和位置變化不大。

在M3分區中，塌陷區2和5在后一時段形變速率降幅較大。塌陷區7為北洺河鐵礦塌陷區，是研究區比較典型的老塌陷區，位于武安市上團城村東北約1 km處，礦區面積20.95萬m2，2002年建成投產，2003年2月出現塌陷坑，2004年2月主要采空區上部約1 000 m范圍內整體下沉2 m左右，塌陷坑周邊地表裂縫增多，最寬超過400 mm，地表傾斜明顯，在北洺河(改道)北岸形成1.24 km2的塌陷區[16]。北洺河鐵礦塌陷區開挖近20年，但本研究顯示該塌陷區目前依然存在形變，兩時段的平均形變速率分別為2.60 mm/d和3.47 mm/d。由此可見礦產開挖所形成的塌陷區具有形變時間長、破壞性大和難治理等特點。

由表2可知，M1塌陷區的形變速率為1.15～11.48 mm/d、平均形變速率為5.96 mm/d，M2塌陷區的形變速率為0.28～9.50 mm/d、平均形變速率為5.78 mm/d，M3塌陷區的形變速率為0.76～8.87 mm/d、平均形變速率為3.02 mm/d，可見M1、M2塌陷區的形變量和形變速率明顯大于M3塌陷區；3個分區中各塌陷區的形變中心位置變化各異，主要與各礦區實際開挖方向有關，在單時段形變累積圖中最大形變中心可視為階段開挖的起始點，形變中心位置不變、形變加快，這可能是垂向開挖造成的。

3 結論

本文采用小基線集合成孔徑雷達干涉測量(SBAS-InSAR)技術，對2016年10月-2017年3月和2017年10月-2018年3月武安市Sentinel-1A影像進行干涉與反演，分別得到兩時段的時序地表累積形變量。結果分析表明：反演塌陷區主要分布在固鎮—武安—礦山—綦村成礦帶和鼓山—洪山—新城成礦帶內，且受礦山斷裂和鼓山—紫山斷裂控制；各塌陷區的范圍和形變中心都隨時間有所變化，并出現個別塌陷區形變中心的消失和新增；2017-2018年最大形變量和平均形變量均比2016-2017年有所減小，這在一定程度上表明研究區塌陷問題整體有所緩解；M1、M2分區塌陷區范圍、累積形變量和形變速率均比M3分區大，建議對M1、M2分區建立長期穩定的InSAR和原位監測體系，加強對塌陷區突發災害問題的預警。由于缺少塌陷區原位沉降監測數據，本文未對反演速率進行精度驗證，未來將結合原位監測手段開展精度驗證及數據融合研究。