基于SBAS-InSAR時序分析技術的甘肅省紅會礦區地面沉降監測及其災害發育特征研究*

王 頌 張路青 周 劍 張奮翔 韓振華 尹凡龍

(①中國科學院地質與地球物理研究所,中國科學院頁巖氣與地質工程重點實驗室,北京 100029,中國)(②中國科學院大學,地球與行星科學學院,北京 100049,中國)(③中國科學院地球科學研究院,北京 100029,中國)(④北京工業大學,城市與工程安全減災教育部重點實驗室,北京 100124,中國)(⑤北京市勘察設計研究院有限公司,北京 100038,中國)(⑥北京中科地環工程技術有限公司,北京 101499,中國)

0 引 言

人類礦山工程與地質環境之間矛盾尖銳,大規模的采礦活動改變了區域能量交換和物質轉移過程,破壞了局部地質體的穩定性,繼而容易誘發地質災害(胡靜云等,2014; 周超等,2020)。礦山地質災害類型眾多,既包括發生于井下的冒頂、片幫、巖爆、突水、突泥等,也包括發生于地表的地面沉降、地面塌陷、地裂縫、崩塌、滑坡、泥石流等(寧宇,2009; 郭維君等,2010)。其中:地面沉降、地面塌陷和地裂縫是地下礦體開采后上覆巖層冒落,圍巖因內部原有應力平衡狀態被改變而不斷移動變形,最終造成地表下沉破壞的災害現象(何芳等,2008,2010; 黃平路等,2010)。高強度礦山開采所誘發的地面沉陷不僅導致了地表水和地下水循環紊亂、植被損傷、礦區沙漠化等一系列生態環境問題,還毀壞了采空區地表的房屋和農田,易使社會矛盾激化(王卓理等,2011)。

甘肅省紅會礦區是受高強度地下開采影響而形成的典型地表沉陷區,地表塌陷坑和地裂縫錯綜復雜,危害顯著。目前已有相關文獻針對礦區地面沉降進行了初步探索。李彥軍(2008)重點論述了礦區地裂縫和塌陷坑的發育特征、成因及其所造成的危害。然而,受當時科學技術條件所限,他所獲取的災害數量和空間分布數據十分有限。隨著遙感技術的飛速發展,大范圍內的地貌和地物得以快速便捷觀測和精確解譯。劉歡等(2010)、丁麗等(2010)和李瑩瑩等(2014)使用高分辨率的IKONOS衛星影像識別并統計了紅會礦區內的塌陷群和煤堆、煤矸石堆等土地占用情況,但沒有涉及到地裂縫的研究。王娟等(2012)和尹展等(2017)基于RS和GIS對紅會礦區地質環境進行了評價和分區,但評價過程中沒有考慮體現地面沉降動態過程的相關指標。因此,目前針對紅會礦區地面沉陷的研究還很不充分,特別是缺乏對災害發育特征的定量描述和礦區地表動態變形的監測。

廣泛使用的地表變形監測技術有精密水準測量、數字攝影測量、GPS實時監測、分布式光纖監測等(衛建東,2007; 董文文等,2016)。合成孔徑雷達干涉測量(InSAR)是一種基于主動式微波遙感區分地表屬性和記錄地表三維形變信息的新型監測手段。和傳統的變形監測技術相比,InSAR具有全天候、大面積、高精度、高分辨率、低成本等突出的優點,適用于火山活動、滑坡活動、冰川運動和地面沉降等諸多研究領域(李珊珊等,2013; 楊澤發等,2015; Rosi et al.,2018; 蘭恒星等,2019; 姚佳明等,2020)。目前,應用最多的InSAR變形監測方法包括D-InSAR、PS-InSAR、SBAS-InSAR、IPTA-InSAR、DS-InSAR和MAI技術等(Tao et al.,2014; 朱建軍等,2017;Xie et al.,2020; 張嚴等,2021)。例如,Chen et al.(2012)基于IPTA-InSAR技術研究了2006～2010年中國珠江三角洲地區地質、人類活動與地面沉降的相互作用,發現和第四紀沉積活動相比,人類城市化活動更容易誘發地面沉降。Zhang et al.(2018)使用D-InSAR技術提取了特拉華盆地鉀鹽礦地面沉降信息,發現地面沉降和鉀鹽礦的開采強度有密切的相關性,并且開采停止后地面仍然在持續下沉。

上述各種方法中,SBAS-InSAR選取空間、時間基線均較小的干涉像對,能有效緩解數據處理過程中的失相干情況,提高了形變監測精度,更適合應用于城市及礦區地面沉降監測(李珊珊等,2013)。張金芝(2015)基于SBAS時序分析技術提取了1992～2000年黃河三角洲地面沉降速度,并重點分析了石油開采、地下水抽取和地表加載等典型影響因素對地面形變的作用。Zhao et al.(2011)監測到山西大同新城建設區、地下水開采區和采煤塌陷區的地面沉降速度分別是1cm/年、2cm/年和4cm/年。Zhang et al.(2019)將SBAS-InSAR形變分析結果與武漢市56個特定水準點測量結果進行對比,結果表明,SBAS-InSAR技術測量精度可達毫米級; 武漢市分布有4個主要沉降區,地下水開采、地鐵開挖和城市建設是造成城市沉降的主要因素。

綜上所述,本文以甘肅省白銀市平川區紅會礦區為研究對象,結合遙感解譯和現場調查定量分析地面變形災害的發育特征,并利用SBAS-InSAR技術系統研究該區域從2017年3月到2019年10月期間地面變形程度和隨時間的演化規律,為地面災害的預警和治理提供依據。

1 地質概況

1.1 區域地質背景

白銀市平川區位于甘肅省的中北部,區域海拔為1347～2858m,整體地形呈階梯狀,東南高、西北低,由東南向西北傾斜,多侵蝕-構造中山、低中山。平川區的東北和東南分別有六盤山和秦嶺作為屏障阻隔暖濕氣流,故區內氣候干燥,年平均降雨量不足200mm。北近騰格里沙漠,西伯利亞寒流易入侵,因此冬季寒冷,年平均氣溫約8.0℃(尹展,2012)。

大地構造上,平川區地處北祁連構造帶東段,秦祁褶皺強烈,經歷了加里東期、燕山期和喜山期多次大型構造運動,表現為強烈的沉降和局部的隆升。區內NNW向的斷裂構造發育,NNE向構造次之,主要有車輪口逆斷裂、響泉口逆斷裂、大水頭逆斷裂等(劉歡,2011)。前寒武系至第四系地層在區內均有出露,其中:第四系覆蓋面積最大,巖性主要為砂土和亞黏土; 白堊系地層巖性為紫紅色砂巖; 侏羅紀陸相含煤沉積形成該區域的含煤地層,可分為延安組和富縣組上、下兩部分,成礦條件好,煤炭資源總量位居甘肅省第二; 區內還出露小面積巖漿巖,主要是志留紀花崗閃長巖和斜長花崗巖(李瑩瑩等,2014)。

1.2 研究區工程地質條件

本文研究區主要為白銀市平川區中北部的紅會礦區及其周邊地區(圖1)。礦區中心經緯度為36°40′32″ N,105°04′47″ E,高程為1902m。煤層賦存以厚煤層為主,埋藏深度200～500m。研究區內有紅會一礦、紅會二礦、紅會三礦、紅會四礦、會通煤礦、長征煤礦、神峰煤礦等多家煤礦企業,均為地下開采。其中:紅會煤礦隸屬于靖遠煤業集團有限責任公司,含煤面積約18km2,是甘肅省重要的煤炭基地之一。礦區內主要涉及的行政村有西合村、牛拜村、玉灣村等,南側有多條省道和縣道穿過,兼有紅會鐵路直達礦山,交通便利。

圖1 研究區概況

研究區地處西格拉灘盆地的西部,地形相對平坦,附近的地層巖性和地質構造如圖2示。該區域屬于寶積山向斜的北翼部分,北西向分布的壓性逆斷層發育(劉歡等,2010)。整個研究區內出露的巖性主要有第四系全新統洪積層(Qhpl)的砂、礫及砂土; 中侏羅統龍鳳山組(J2l)的含礫砂巖、細礫巖; 上三疊統西大溝組(T3x)的砂巖夾粉砂巖、頁巖; 志留紀花崗閃長巖(Sγδ)等。由于坐落在祁連山斷裂帶和西海固斷裂帶的交會處,研究區地震活動頻繁,根據中國地震動參數區劃圖(中華人民共和國國家質量監督檢驗檢疫總局等,2015),區內地震烈度為Ⅷ 度,地震動加速度反應譜特征周期為0.45s,地震動峰值加速度為0.3g。

圖2 紅會礦區周邊地質圖

礦區地表僅有數條因溝谷匯聚而形成的沙河,在暴雨時可形成洪流并最終匯入黃河。礦區地下水主要為賦存于白堊系、侏羅系砂礫巖中的碎屑巖類孔隙裂隙水和賦存于山前沖洪積物中的松散巖類孔隙水。研究區內主要人類活動是礦產資源開采,農耕等。受地形條件限制,區內崩塌、滑坡、泥石流等山地災害不發育。但由于大規模高強度的采煤活動,地面塌陷和地裂縫極其發育。

2 災害發育特征分析

紅會礦區歷史開采形勢復雜,尤其是20世紀中后期有近百家地方及私營企業同時進行近地表煤炭開采,不僅破壞了煤炭賦存的整體性,造成大量資源的浪費,而且使用落后無序的開采方式采煤后在井下形成了大范圍不規則采空區,誘發了地面沉降、塌陷坑和地裂縫等地質災害。盡管進入21世紀后,開采技術大幅提升,規?；?、機械化開采迅速推廣,但能源需求的增加迫使煤礦開采強度不斷加大,使得礦山地面地質災害加劇。

結合遙感影像,解譯和統計了截至2019年紅會礦區地面塌陷坑和地裂縫的空間分布,尺寸和數量,如圖3和圖4所示。在紅會礦區地表共發現地裂縫735條,塌陷坑559個,地裂縫平均長度約200m,塌陷坑平均面積約360m2。其中:地裂縫長度主要分布在100～300m范圍內,塌陷坑面積主要分布在<400m2的范圍內。從圖3a中還可以看出,地裂縫和塌陷坑密集分布在各個煤礦周圍,呈現出明顯的方向性,大致沿北東向分布,與該區域煤層的分布基本一致。

圖3 紅會礦區塌陷坑和地裂縫空間分布

圖4 紅會礦區塌陷坑和地裂縫數量及尺寸分布特征

此外,由于紅會一礦長年運用走向長壁綜采放頂煤采煤方法進行地下采煤活動,較薄的碳酸巖頂板大面積垮落,致使平川區黃嶠鄉牛拜村附近成為受地面塌陷影響最為嚴重的區域。調查發現,牛拜村地面沉陷區長約2.6km,寬約1.5km,影響面積約3.9km2,且隨著該區域地下開采的推進,礦山地面災害數量呈逐年遞增的趨勢。沉陷區內的部分塌陷坑深4m,半徑約15m,面積可達700m2; 地裂縫縱橫交錯,并穿過塌陷坑(圖5)。整個牛拜村周圍地面沉降大、塌陷群分布廣,造成周圍村鎮的房屋、耕地、道路破壞嚴重,牛拜村村民的生命和財產安全無法得到保障(圖6)。

圖5 牛拜村地面沉陷區內的塌陷坑與地裂縫照片

圖6 牛拜村地面塌陷產生的主要危害

針對紅會礦區一系列的礦山地質環境問題,白銀市常年開展環境恢復治理工程,其中:2012年和2013年開展的兩期重點治理工程(圖3b),主要包括紅會一礦周圍區域煤矸石堆、廢煤堆平整,挖土填方,土地復墾,種植果樹等,并且對周圍房屋損害嚴重的部分居民進行搬遷。通過對比治理工程前后的遙感影像并統計相關數據,發現紅會礦區內的地裂縫減少183個,塌陷坑減少108個,工程效果顯著(圖3b、圖3c和圖4)。然而,治理工程的范圍很小,約為整個沉陷區的十分之一。盡管平川區人民政府在2019年的《白銀市平川區采煤沉陷區綜合治理行動實施方案》中指出,計劃在未來5年實施兄弟村、中和村、西合村群眾1017戶4058人和牛拜村群眾95戶528人的搬遷安置任務,并繼續推進采煤沉陷區的生態修復和重大地質災害隱患治理(平川區人民政府,2019); 但隨著采礦活動的繼續進行,礦山地面地質災害仍有進一步發展的趨勢,因此,十分有必要對該區域的地表形變進行跟蹤監測,以支撐災害的預警和治理。

3 小基線集(SBAS)技術

3.1 基本原理

在SBAS技術中通過使用每次D-InSAR計算的形變結果即可獲取高精度的形變時間序列(Berardino et al.,2002; 李凌婧等,2014)。首先,選取N+1景覆蓋同一區域、時間分布為t0,t1,…,tN的SAR數據影像圖,進行基線估算后,選取其中1景作為超級主影像,其余N景為輔影像,將時空基線均小于某一閾值的影像進行配對,共形成M對干涉像對,其中M滿足:

(1)

將M對干涉像對采用雙通道D-InSAR進行差分干涉處理,去除平地相位和地形相位后,假設在tA和tB時刻獲取的兩景SAR影像生成第i幅解纏的多視差分干涉圖,那么SAR坐標系下任意坐標(x,r)處相干點的相位值可以表示為(Guo et al.,2017):

Δφi(x,r)=φ(tB,x,r)-φ(tA,x,r)

≈Δφi,def(x,r)+Δφi,topo(x,r)+

Δφi,atm(x,r)+Δφi,θ(x,r)+

Δφi,nosie(x,r)

datm(tA,x,r)]+Δn

(2)

式中:λ為雷達波長; Δφi,def(x,r)是累積形變引起的相位; d(tA,x,r)和d(tB,x,r)分別是tA和tB時刻相對于參考時間t0時刻在雷達視線方向(LOS)上的累積形變相位量; Δφi,topo(x,r)是由地形誤差所導致的相位,與垂直基線B⊥和高程誤差Δz成正比,與入射角θ的正弦成反比; Δφi,atm(x,r)是兩幅影像中大氣延遲引起的相位datm(tA,x,r)和datm(tB,x,r)所造成的相位誤差;Δφi,noise(x,r)為噪聲殘余相位; Δn表示失相關引起的噪聲相位。其中:大氣相位誤差和噪聲誤差可以通過時域高通濾波和空域低通濾波與干涉相位分離。

在去除上述各種相位誤差分量后,由式(2)可以得到M個方程,方程中含有N個未知量,用矩陣表示為:

Aφ=Δφ

(3)

式中:A對應一個M×N的系數矩陣,每一行和每一列分別對應一個干涉圖和一景SAR影像,因此M和N分別為干涉圖和SAR影像的數量。

直接求解相位將違反形變的物理規律,故將兩景SAR影像獲取時間之間的平均相位速度作為未知量來進行求解(李珊珊等,2013),

(4)

式中:φ(t0,x,r)=0。

結合式(3),(4),可以得到一個新的矩陣方程:

Bv=Δφ

(5)

式中:B為M×N的系數矩陣。采用奇異值分解法(SVD)和最小二乘法進行求解方程(4),即可得到形變速率,進一步在時間域上積分得到整個時間序列的形變相位值。

3.2 數據介紹

本文地表形變監測選用由歐洲空間局(ESA)發射的Sentinel-1A衛星獲取的覆蓋研究區域的32景單視復數影像(SLC)作為實驗數據,數據范圍如圖1。表1顯示的是所用SAR影像的基本參數情況,其中時間跨度為2017年3月27日到2019年10月19日,影像數據的極化方式為VV,軌道方向為升軌,分辨率為5m×20m。為了提高影像的軌道數據精度,本研究還采用了歐洲空間局發布的AUX＿POEORB精密軌道參數文件進行軌道校準。此外,所用數字高程模型數據為美國宇航局(NASA)公開的SRTM數據(圖7),分辨率為30m×30m,用于消除地形相位的影響。

表1 Sentinel-1A 影像技術參數

3.3 數據處理流程

利用SBAS-InSAR技術進行紅會礦區地表形變監測的數據處理流程如圖8所示。32景數據裁剪后依次以每一景數據作為超級主影像,進行基線估算,其中最大的絕對空間基線為254.409m,最小的絕對空間基線為1.665m,最大的絕對時間間隔為936d,最小的絕對時間間隔為24d。設置最大空間基線為臨界基線的45%,最大時間基線為550d,計算各組像對的空間基線和時間基線如圖9。

圖8 SBAS-InSAR數據處理流程

圖9 時空基線分布圖

根據時空基線圖設置2018-05-21的影像為超級主影像,其他影像與其配準,經差分干涉生成、軌道精煉和重去平、形變速率和殘余地形估算、地形殘余相位和大氣相位去除、形變序列估算、地理編碼等多項處理步驟后獲得最終的形變結果。

4 紅會礦區地表形變分析

紅會礦區各時間段地表變形的累積量及變形速率如圖10a。從圖中可以看出,計算區域內累積沉降量較大的位置集中在劃定的研究區內,2017年3月27日～2019年10月19日近兩年半的時間里,研究區內累積最大沉降量可達177mm,最大平均年沉降速度達70mm·a-1?？焖俪两祬^主要集中在紅會一礦、紅會二礦和紅會四礦周圍,沉降量從沉降中心向周圍遞減,共和鎮、牛拜村、西合村和玉灣村均受地面沉降災害影響。將計算范圍的結果圖進行裁剪得到研究區內的結果圖(圖10b),并與遙感解譯的地裂縫和塌陷坑分布圖相疊加,可見,紅會礦區現有的地裂縫幾乎全部分布在快速沉降區,隨著煤礦采空區的進一步增大,極有可能繼續擴展。多數塌陷坑則分布在緩慢沉降區,變形相對穩定,進一步災變的可能性很小。盡管環境恢復治理區沒有新的地裂縫和塌陷坑產生,但依舊在遭受采礦活動的劇烈擾動,沉降較為明顯。

圖10 SBAS-InSAR時序分析獲取的2017/03/27～2019/10/19期間的累積沉降量和沉降速度(雷達視線方向)

以2個月為時差展示地面沉降過程,圖11所示的研究區形變時間序列。以2017年3月27日為起點,假設此時地面沉降為0,則從圖中可以看出,隨著礦山開采活動的進行,地面沉降逐漸嚴重,最先呈現紅色的即為沉降中心,并隨著時間開始向周圍發展。沉降最為嚴重的為紅會一礦和紅會二礦周圍,這也符合礦山生產的實際情況。

圖11 研究區形變時間序列圖

為了進一步研究地面沉降速度及其對周圍村鎮、公路的影響,如圖3a所示,選取典型的監測點(P1～P22),給出其整個時間范圍內的沉降曲線和距離研究區約30km的靖遠縣雨量站日降雨量數據(圖12)。其中:P1和P2位于2012～2013年的環境恢復治理區,解譯結果表明,該區域仍以每年40～50mm的速度沉降,治理效果難以維持。P3和P9代表紅會一礦和二礦礦區內沉降比較嚴重的位置點,在兩年半的時間內,幾乎以線性速度沉降,最終累積沉降量分別達130mm和120mm。P4～P7所在區域為歷史塌陷嚴重區,該區域還在不斷沉降,累積沉降量達80mm,但是其不是目前沉降的中心區域,沉降速度呈降低趨勢。P8位于快速沉降區的邊緣地帶,沉降速度大幅降低,年均沉降速度小于10mm。結合P10～P12空間位置和時序變形數據,表明了距離開采中心越遠,地面形變強度逐漸減小。P1～P12所在區域原本為農田,由于地裂縫和塌陷坑導致無法灌溉,已基本荒廢。

圖12 紅會礦區典型監測點時序變形(雷達視線方向)

P13為紅會三礦周圍的位置點,沉降量很小,不到10mm。P14點位于紅會四礦礦區,累積沉降量達70mm,但在2018年7月之前沉降速度較大,之后沉降緩慢,這主要因為四礦資源枯竭,生產減慢,即將面臨關停。P15神峰煤礦沉降量約20mm,P16西合村距離紅會四礦較近,近兩年沉降量很小,不到10mm。P17、P18和P19分別為牛拜村、玉灣村和共和鎮位置點,累積沉降量分別為30mm、15mm和25mm。P20、P21和P22為省道308上面的3個位置點,沉降量相差不大,為20～25mm,調查時路面未見明顯裂縫,但已立警示牌告知注意路面塌陷。

觀察P17～P22監測點時序變形曲線,可見6條曲線的形態大致相同,幾乎都在每年7月累積沉降量達到峰值,之后沉降變形減弱,地面甚至開始抬升?？紤]到降雨主要集中在每年的6、7、8月份,這說明這些人口聚集區的地面沉降變形,一定程度上還與地下水過量開采有關,雨季時雨水補給地下水,使得原本沉降的地面又有所回升。但從整體沉降趨勢來看,這些區域仍處于緩慢下沉階段。值得注意的是,本文所使用的SAR影像數據僅包含了近兩年半的地表變形信息,結合遙感歷史數據和現場調查發現,牛拜村、西合村、兄弟村等村莊及農田區域地面沉降嚴重,伴有大量的塌陷坑和地裂縫,這是煤礦數十年高強度地下開采導致的結果。因此,十分有必要盡快推進沉陷區群眾搬遷安置以及生態修復、環境整治和重大地質災害隱患治理等工作。

圖A1 SBAS-InSAR時序分析獲取的2017/03/15～2019/11/12期間的累積沉降量和沉降速度(雷達視線方向)

5 結 論

(1)依托遙感解譯和現場調查,獲取了甘肅省紅會礦區地面塌陷坑和地裂縫的空間分布,尺寸和數量,結果表明,紅會礦區數十年高強度的地下開采導致至少700條地裂縫和500個塌陷坑沿北東向密集分布于礦區地表,造成牛拜村、西合村、兄弟村等村莊的房屋嚴重受損,農田無法灌溉。

(2)基于2017年3月至2019年10月期間的32景Sentinel-1A數據,利用SBAS-InSAR技術獲取了紅會礦區的地面平均形變速率與累積沉降量,并針對典型位置點進行了時序變形分析。結果發現,快速沉降區集中在紅會一、二、四礦周圍,與地裂縫空間分布位置相吻合,累積最大沉降量達170mm; 環境恢復治理區仍受采煤活動的劇烈擾動,以年均40～50mm的速度沉降。

(3)牛拜村、西合村、兄弟村等人類聚集區的地表變形整體呈緩慢下沉趨勢,累積沉降量約15～30mm; 這些區域的地面沉降除受地下采煤影響之外,還可能與地下水過量開采有關,雨季時原本沉降的地面有明顯回升。

(4)整個紅會礦區地面變形災害嚴重且處于發展態勢,應規范礦業行為,采用合理的采礦技術方法,進一步控制災害的惡性發生,同時盡快推進沉陷區群眾搬遷安置、生態修復和地質災害隱患治理等工作。

附 錄A

在未能獲取同時間階段內現場實測地表形變數據的條件下,使用SBAS-InSAR技術和位于不同軌道時的數據再次提取了研究區域內的地表變形信息,以對前文所述的監測結果進行交叉驗證。其中:所采用的數據為2017年3月15日到2019年11月12日期間33景降軌的Sentinel-1A數據,影像數據的極化方式為VV,其他參數與3.2節中所述的一致。最終,獲得的紅會礦區各時間段地表變形的累積量及變形速率如圖A1所示。

對比圖A1和圖10a,在兩種不同類型觀測數據模式下,兩者的累積最大沉降量和年均沉降速度有小幅差距,這可能是因為最終獲取結果為雷達視線方向形變而不是真實三維形變,故而在一定程度上受雷達觀測方向的影響。但整體上,研究區域內兩者的地表形變基本吻合,快速沉降區的空間分布保持一致,集中在紅會一、二、四礦周圍,證明了監測結果可靠性較高。