基于SBAS-InSAR技術(shù)的武安塌陷區(qū)形變監(jiān)測(cè)分析

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

(1.河北省科學(xué)院地理科學(xué)研究所/河北省地理信息開(kāi)發(fā)應(yīng)用工程技術(shù)研究中心，河北 石家莊 050011；2.河北交通職業(yè)技術(shù)學(xué)院，河北 石家莊 050035)

0 引言

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

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

1 數(shù)據(jù)與方法

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

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

圖1 研究區(qū)概況Fig.1 Survey of the study area

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

表1 研究區(qū)塌陷區(qū)與地裂縫災(zāi)害點(diǎn)和反演結(jié)果驗(yàn)證Table 1 Subsidence areas and ground fissures in the study area and verification of inversion results

1.2 SBAS-InSAR技術(shù)

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

2 結(jié)果分析

2.1 整體空間分布特征

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

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

利用ArcGIS將17處災(zāi)害點(diǎn)(表1)與累積形變量空間分布圖疊加(圖3)對(duì)比，發(fā)現(xiàn)12處災(zāi)害點(diǎn)分布在反演結(jié)果塌陷區(qū)內(nèi)，聚隆煤礦周?chē)亓芽p、豐里村及周?chē)亓芽p和西萬(wàn)年村塌陷區(qū)位于反演塌陷區(qū)的影響范圍內(nèi)，河西村西山坡塌陷區(qū)和龍洞山塌陷區(qū)兩處未顯示形變。由此可見(jiàn)，本文反演結(jié)果較準(zhǔn)確地反映了區(qū)內(nèi)災(zāi)害點(diǎn)的分布。

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

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

2.2 局部空間特征

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

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

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

在M2分區(qū)中，塌陷區(qū)6、9在兩個(gè)時(shí)段平均形變速率相近，規(guī)模變化不大，但形變中心均向南轉(zhuǎn)移；塌陷區(qū)1、4、5、8在后一時(shí)段平均形變速率減小，表明開(kāi)采強(qiáng)度有所減弱，其中塌陷區(qū)4形變速率最小，形變中心明顯消失，塌陷區(qū)1和8范圍有所增大，塌陷區(qū)5形變速率明顯減??；塌陷區(qū)2、3、7在后一時(shí)段平均形變速率增大，其中塌陷區(qū)2范圍增大，塌陷區(qū)3范圍減小且形變中心南移，塌陷區(qū)7規(guī)模和位置變化不大。

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

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

3 結(jié)論

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