基于SABS-InSAR的保德礦區地表形變監測與分析

樊昱初, 呂義清*, 楊娜

(1.太原理工大學礦業工程學院, 太原 030024; 2.航天宏圖信息技術股份有限公司, 北京 100094)

近年來由于煤炭能源的需求量逐漸增大,特別是資源整合后,礦山生產能力得到進一步的提高,保德礦區內形成了大量的采空區。在地下煤炭資源開采的過程中,引起了采場周圍巖體的應力重新分布[1],開采工作面及其周邊區域失去原始應力平衡。隨著采空區范圍的擴大,巖層移動逐步向上傳遞,最終導致地表發生形變,局部出現地面塌陷和地裂縫。破裂煤巖體在應力、地下水等多種因素影響下,其強度、承載力會有一定程度的降低,加大了地質災害發生的風險。為了保障礦區周邊居民的人身財產安全和可持續開采,對于礦區進行沉降監測分析具有重要的研究意義。利用技術手段掌握礦區地表形變信息,對潛在的災害點進行識別,為后續的防治工作提供科學的指導。

傳統的觀測方法有水準測量、三角高程測量、全站儀測量、靜態全球導航衛星系統(global navigation satellite system,GNSS)和動態全球導航系統等[2],這類型散點式觀測手段雖然有較高的精確度,但無法對整個礦區進行大面積,長時間的監測。合成孔徑雷達干涉測量(interferometry synthetic aperture radar,InSAR)作為一種新興的高精度的測量技術,能夠很好的解決宏觀監測的問題,其中合成孔徑雷達差分干涉測量技術(differential interferometry synthetic aperture radar,D-InSAR)可以全天時、全天候、近實時獲取大面積地表地形信息,且空間分辨率高,受大氣和季節影響較小。王志勇等[3]利用D-InSAR技術對濟寧礦區的沉降做了精細化的監測和分析,證明D-InSAR技術非常適合礦區大型變的沉降監測,監測精度可以達到厘米級。但是D-InSAR技術無法獲取研究區域長時間的時序變化結果,而且易受到時空失相干和大氣延遲效應影響,為了優化這一領域受限的問題,Berardino等[4]在2002年提出了小基線集技術(small baseline subset,SBAS)該技術不僅能夠有效的削弱時空失相干和大氣延遲效應對監測結果的影響,還能以更好的時空連續性來獲取相干點目標信息,使監測結果更加精確和直觀。朱建軍等[5]詳細介紹了InSAR技術目前的發展和主要應用領域,不僅可以對城市的沉降、礦區地表、地震、基礎設施等進行形變的監測,還能監測火山,冰川等的活動,同時對滑坡災害也有一定的識別和預警作用;劉志敏等[6]利用SBAS-InSAR技術對長治礦區的地表形變進行了監測和分析,并揭示了老采空區的殘余沉降趨勢;麻學飛等[7]先通過D-InSAR技術識別臨汾礦區沉陷區位置,后利用SBAS-InSAR技術進行時序監測,共獲取105處沉陷,為后續礦區沉陷災害防治提供了重要依據。

SBAS-InSAR在地表形變監測[8-9]以及地質災害預警[10]等方面已經有較為成熟的應用體系,但在礦區地表形變監測中,利用InSAR技術對形變驅動性因素的研究尚不充分,因此在此基礎上采用SBAS技術對研究區32景Sentinel-1A影像數據進行處理,獲取了保德礦區2年時間內的地表形變信息,結合野外勘查與相關資料,選取礦區內地質構造、地層巖性、采礦工程、降雨等因素,定性分析其對研究區內地表形變程度及范圍的影響,明確形變誘因,為煤礦安全生產風險性評價提供具體的評價指標,同時識別了潛在的地質災害風險點,為進一步的開采以及地質災害的防治提供了科學依據。

1 研究區概況與數據來源

1.1 研究區概況

研究區域所在的保德縣位于山西省西北部、地理范圍為38°39′00″N～39°6′56″N、110°56′30″E～111°19′40″E(圖1)。整體地形東高西低,北高南低,平均海拔為840 m,區域內地表主要被黃土覆蓋,水土流失、地形切割較為嚴重。土壤結構松散,溝壑縱橫,梁峁連綿,屬于典型的黃土丘陵溝壑地形地貌。大地構造位置處于鄂爾多斯東部大型聚煤盆地東緣河曲—保德斷隆以北,山西斷隆以西[11]。區內主要發育近南北走向的斷裂,整體構造比較簡單。

圖1 研究區位置Fig.1 Location of the study area

1.2 數據來源

選用覆蓋保德區2020-06-24—2022-09-13期間32景Sentinel-1A升軌單視復數(single look complex,SLC)影像作為基礎數據,其中path=11,frame=121有24景,path=113,frame=121有8景。衛星對地觀測模式為干涉寬幅模式(IW),極化方式為VV,重訪周期為12 d。獲取影像均為C波段,距離向分辨率為5 m,方位向分辨率為20 m。采用歐空局發布的精密定軌星歷數據(precise orbit ephe-merides,POD)對軌道信息進行修正,引入空間分辨率30 m的SRTM數據作為數字高程模型(digital elevation model,DEM)去除地形誤差相位。同時獲取的還有2020-06—2022-09期間保德縣境內五個氣象監測點的降水數據,研究區地質數據,以及重點監測區內各煤礦采掘工程信息。

2 數據處理方法

2.1 SBAS-InSAR技術原理

SBAS-InSAR是一種基于多主影像的InSAR時間序列技術,基于短基線原則,將所研究的SAR數據轉換組合為有多個主影像的干涉子集,并選擇其中一景影像作為公共主影像進行配準,利用時空基線較短的干涉對來提取地表形變信息。將在t1～tm時間內獲取的m幅SAR影像進行干涉組合,并對集合間影像進行差分干涉處理,形成N幅干涉條紋圖,滿足的條件為

(1)

對tA和tB兩個時間點的影像生成的第i幅干涉圖中,第x個像素的干涉相位表達式為

φint,x,i=φdef,x,i+φtopo,x,i+φaps,x,i+φorb,x,i+φnoise,x,i

(2)

式(2)中:φdef,x,i為斜距向變形即形變相位誤差;φtopo,x,i為地形相位誤差;φabs,x,i為大氣延遲相位引起的誤差;φorb,x,i為基線軌道引起的軌道相位誤差;φnoise,x,i為噪聲誤差。各誤差相位可在時間序列上采用濾波方法或多項式模型予以去除。

由最小二乘法求解得到監測區地表時間序列形變信息,假設不同干涉圖之間的變形速率vk,k+1,在tA～tB時間內的累積形變表達式為

(3)

將得到的N幅干涉圖進行解纏處理,可以進一步求出不同SAR影像獲取時間的形變速率,以上提及處理方式的具體計算方法可以參考文獻[4]。

2.2 數據處理流程

使用ENVI/SARscape軟件對數據進行計算處理,利用Arcgis[12]軟件對結果進一步分析,具體流程如下。

(1)為了提高影像數據的處理效率,首先根據礦區矢量范圍裁剪影像,通過設置時空基線閾值和多視比例值(距離向為5,方位向為1),生成連接圖和干涉對組合,選取2021-07-01影像為超級主影像,其他影像依次與主影像配對。

(2)干涉工作流處理,通過導入外部DEM來去除干涉圖中的地形相位,選用Goldstein方法對得到的差分干涉圖進行濾波處理,采用基于Delaunay三角網的最小費用流法(minimum cost flow,MCF)對差分干涉圖進行相位解纏,減小大氣延遲誤差、基線軌道誤差、相干斑噪聲誤差等對計算結果的影響,提高干涉條紋的清晰度。

(3)軌道精煉和重去平,在地形平坦遠離主形變區,或已知形變的區域設置36個控制點,通過這一步估算和去除殘余的恒定相位和解纏后還存在得相位坡道。

(4)采用奇異值分解(singular value decomposition, SVD)法估算平均形變速率和殘余地形相位,并且對合成的干涉圖去平,重新進行相位解纏和精煉,生成更加優化的結果。通過進一步反演,進行定制的大氣濾波,估算去除大氣相位,得到更精確的時間序列位移結果。

(5)將生成的結果進行地理編碼,得到研究區雷達視線向(line of sight,LOS)的平均形變速率和累積形變量。

3 時序監測結果與分析

3.1 礦區地表形變結果分析

3.1.1 形變空間分布特征

在2020-06-24—2022-09-13監測時間段內,通過SBAS-InSAR技術處理得到研究區的年平均形變速率(圖2)和累積形變量(圖3),其中負值表示下沉,正值表示抬升。

圖2 保德礦區地表年平均形變速率Fig.2 Average annual surface deformation rate in Baode mining area

圖3 保德礦區地表累積形變Fig.3 Cumulative surface deformation in Baode mining area

結合相關資料分析,研究區主要變形區域與區內煤礦工作面位置空間分布一致,K1～K6重點監測點位于各沉降區域中心位置,監測時間段內中心點K5最大的沉降速率為-43.7 mm/a,最大累積形變量為-84.1 mm。從累積形變圖可以看出,研究區整體形變趨勢以下沉為主,采礦工程影響范圍內區域的平均沉降速率為-20.2～-5.5 mm/a,沉降量為-39～-12 mm;累積形變量為-23.1～-4.7 mm的區域在研究區內分布范圍最廣,該區域內分布有多個自然村落,地面沉降主要受各種人類活動影響;研究區部分邊緣地區在監測期內有較為明顯的抬升,平均形變速率為0.4～21.3 mm/a,最大累積抬升為38.1 mm。

3.1.2 礦區形變時空演變特征

以2020-08-23影像為基準,共選取9景影像數據來研究礦區地表形變的空間范圍演變特征(圖4)。

圖4 保德礦區地表形變時間序列圖Fig.4 Time series map of surface deformation in Baode mining area

以圖1重點研究區位置為基本劃分,各重點研究區內的特征點即重點監測點的累積形變曲線(圖5)為參考,分析形變時間序列圖,隨著時間的推移,礦區形成多處由中心向四周發展的漏斗狀下沉區域,整體沉降分布面積先增加后減少。在監測時間段內,形變量為2.0～38.1 mm的區域范圍有顯著增加,主要分布在監測區域邊緣;沉降主體外圍形變量為-4.7～2.0 mm的區域處于較為穩定的狀態,范圍基本保持不變;形變量為-23.1～-4.7 mm的區域范圍在前期有一定程度的增加,在2021-02-19～2021-09-11期間持續減小后逐漸趨于穩定;形變量為-84.1～-23.1 mm的沉降中心區域在監測時間段內面積有明顯的增加;隨著采礦工作面的推進,采掘工程對地表的影響范圍也進一步擴大,形變量為-39～-23.1 mm的區域面積也在持續增加。通過分析時間序列圖可以得出在未來一段時間內,各重點區域內沉降會持續進行,主形變區域的沉降范圍發展和沉降速率將會進一步加快。

圖5 重點監測點形變曲線圖Fig.5 Deformation curves of key monitoring points

從重點監測點K1～K6曲線(圖3)整體來看,累積形變量與時間基本呈線性關系,各點形變均以下沉趨勢為主,由于所處區域煤礦開采工程工況區別較大,所以各點的累積形變量大小和沉降速率都有明顯的差距,其中K1點累積沉降量為-38.36 mm,K4點累積沉降量為-81.13 mm,最大形變量絕對值差為42.77 mm。

3.1.3 重點研究區形變分析

選取煤礦分布集中,形變程度較大的6個重點研究區(圖1)進行詳細分析,將研究區年平均形變速率圖與Google衛星影像相結合(圖6)。其中有部分關停煤礦以及個人經營煤礦由于資料缺失,不計入分析。

圖6 重點研究區年平均形變速率圖Fig.6 The average annual deformation rate in key research areas

研究區內主要涉及12個產能為百萬噸級以上的煤礦,其中1號研究區主體為南茆露天煤礦,其年平均形變速率最小,為17.64 mm/a,區域內累積最大形變量為-68.93 mm;5號研究區內年平均形變速率為-27.36 mm/a,沉降程度最為嚴重,區域內累積最大形變量為-84.87 mm,形變主體為泰安煤礦和孫家溝煤礦。通過分析研究區的地表形變時間序列圖和重點區域的形變速率圖可以得出,煤礦的采掘工程和相應的人類工程活動是礦區出現局部沉降,形成沉降中心的主要因素。

3.2 研究區地表形變驅動因素分析

為了進一步研究控制和影響保德礦區地表形變的因素,下文將從礦區的地層巖性和地質構造,煤礦的開采工程以及研究區內降水量這4個方面進行具體的分析。

3.2.1 開采工程

以5號研究區中的泰安煤礦生產井田為例(圖7),井田內含煤地層主要為石炭系上統太原組(C3t)和二疊系下統山西組(P1s),其中有5層全區可采煤層,自上而下為山西組的8號煤層和太原組的11、12、13、15號煤層,目前煤礦主要開采8、11、12號煤層。

圖7 泰安煤礦采掘煤層工作面位置Fig.7 Location of coal seam working face in Tai’an Coal Mine

(1)8號煤層。該煤層井田內西南部大面積采空,煤層厚度為3.31～6.55 m,平均4.25 m,煤層賦存標高為1 100～720 m,其中8101工作面于2013年采空,8102工作面于2014年采空,8103、8104工作面于2015年采空,8105、8106工作面于2016年采空,8107工作面于2017年采空。

(2)11號煤層。該煤層在井田中南部部分采空,上距8號煤層23.9～44.48 m,平均為35.72 m,煤層厚度為1.30～1.95 m,平均為1.53 m,煤層賦存標高為1 060～680 m,其中11101工作面于2017年采空,11102工作面于2018年采空,11103工作面Ⅴ區在2019年回采347 m,Ⅵ區在2018年回采588 m。

(3)12號煤層。該煤層井田內東部大面積采空,上距11號煤層2.43～7.84 m,平均為3.63 m,煤層厚度1.00～6.50 m,平均為3.21 m,煤層賦存標高1 050～690 m,其中東部大面積采空區于2012年采空,12101工作面Ⅶ區在2019年采空,Ⅷ區在2018年采空,12102工作面Ⅴ區在2020年采空,Ⅵ區在2019年采空。

將各煤層采礦工作面與圖6中泰安煤礦礦區內地表年平均形變速率圖對比可以看出,地表發生變形的區域與采掘工作面空間分布相對應,其中11號煤層中的11104、11105工作面以及12號煤層中12102～12106工作面均在2020年后有采掘工程,2020-06—2022-09,泰安礦區內局部最大累積形變量達為-84.87 mm,12103、12104工作面對應空間的地表平均形變速率可達-33.25 mm/a;8號煤層中的8106及8107工作面分別于2016年和2017年采空,在監測時間內,其累積最大形變量為-30.91 mm,對應區域的平均形變速率為-7.29 mm/a;12號煤層東部采空區于2012年已完成開采,但在監測時間內,仍能觀測到較為明顯的地面沉降,其中最大累計形變量為-23.11 mm,采空區域內平均形變速率為-4.35 mm/a。根據相關文獻資料,采空區在形成十年后,殘余沉降基本消失,完成總變形量90%以上[13-14],地表一般已經進入相對穩定的狀態,造成東部采空區異常沉降的原因推測是因為采空區與12號煤層活躍開采工作面距離較近,且受到地下水抽取以及地下礦井巷道工作影響,導致其仍存在一定程度的下沉。通過對比正在開采、采空三到四年以及開采完成十年以上不同工況工作面對應區域的形變速率可以得出:人為的采礦工程是造成和加劇礦區地表形變的主要原因。

3.2.2 地層巖性與地質構造

地層巖性與地質構造是地面沉降發生的重要地質背景,研究區內斷層構造的分布以及地面可壓縮層特性等因素對沉降的發展都有一定程度的控制[15-16]。因此研究區域構造背景是分析礦區地表形變驅動性因素的重要一環。保德地區位于鄂爾多斯盆地東緣的北部,以升降運動為主體,整體表現為面狀隆起與局部相對下沉。

為了進一步研究斷層構造對礦區地表形變的影響,將構造圖(圖8)與礦區地表形變速率圖(圖9)進行疊加分析,可以直觀看出礦區整體范圍內沉降發展趨勢與斷層方向一致,斷層兩側的地面沉降速率受斷層的上、下盤滑動及相對升降運動的影響有較為明顯的差異,另一方面斷層兩側地層巖性即可壓縮層厚度的不同,也是導致沉降速率變化的原因之一。根據形變速率圖和監測點數據計算顯示,斷層兩側年沉降速率的平均差異約為3.1 mm/a;保德縣地處晉西北黃土高原,其黃土分布范圍較廣,第四系堆積物占總面積的69.98%。根據巖層的成因類型、結構和工程地質特征,區內的工程地質巖組劃分為3個大類,9個亞類(表1)。

圖8 保德礦區地質構造要圖Fig.8 Geological structure map of Baode mining area

圖9 保德礦區地表形變與斷層分布Fig.9 Surface deformation and fault distribution in Baode mining area

研究區內砂土由第四系全新統(Q4)沖、洪積物組成,巖性為砂卵礫石和砂層,屬于低壓縮性土。黏性土主要為新近系保德組紅色黏土(N2),底部為礫巖,壓縮系數為0.02～0.082 MPa-1。特殊性巖土則按照形成時間分為新黃土和老黃土,據土工試驗,老黃土顆粒成份中黏粒含量高于新黃土,其壓縮系數為0.005～0.132 MPa-1,新黃土則為0.002～0.076 MPa-1。相關研究表明,在相同開采工況下,以粉土、黏性土層為主的區域,單位沉降量是以砂層、砂礫石層為主要巖性區域的2～3倍[17],以碳酸鹽巖帶為基巖的區域的平均形變速率要快于非碳酸鹽巖帶區域。在研究區內第四系堆積物分布廣泛的地質背景下,相關的人類地面活動,地下工程的開發都會導致軟土層的壓實。此外,人工開采引起的周邊區域內地下水水位持續波動,不僅會導致含水砂層壓實還會進一步潛蝕碳酸鹽巖[18],這也是開采區附近區域沉降速率與其余地區有較大差異的重要原因之一。

3.2.3 降雨影響

據保德縣氣象局1980—2020年共計40年統計資料顯示,區域多年平均降水量為475.9 mm,且年際變化較大,豐水年平均降水量為861.4 m,枯水年平均降水量為194.8 m,其中70%的降水量集中在7、8、9三個月,并且多以暴雨出現。為了研究降雨量對礦區地表形變的影響,向相關氣象部門收集了研究區域內自2020年6月起5個地面氣象監測點逐月降水資料,涉及義門鎮、橋頭鎮、孫家溝鄉、南河溝鄉四個區域,由于鄉鎮之間的降水量有一定的差異,因此分別進行了單獨統計,將重點研究區內監測點(圖6)的地面累積形變量與監測時間內區域月降水量進行疊加分析(圖10)。

圖10 保德礦區降水量與累積地表形變Fig.10 Precipitation and cumulative surface deformation in Baode mining area

2～6號重點監測區內,雖然各煤礦采礦的工況有所不同,但從整體形變量走勢分析,降水量的增加會減緩地表形變的發展。礦區地表主要被黃土覆蓋,黃土由粉土、粉質黏土組成,透水性一般較差,因此在降雨強度較小的4、5、10月時間段內,一定量的雨水滲入土壤,淺層土壤中水分含量會有顯著增加[19],此時體現在地表形變圖中,監測點下沉趨勢減緩,部分測點由于土壤吸水膨脹,地表會發生小幅度抬升。而在此后的11月到來年2月期間,因水分蒸發等原因,表層土壤的含水量會進一步下降,此時地表會發生季節性沉降。而在7、8、9月,由于多暴雨且降雨時間集中,降水會在短時間內匯集,形成的地表水流會在黃土構造節理、風化裂隙、陷穴等發育部位的空隙下滲或灌入,在相對隔水的部位形成上層滯水或飽水帶,增大巖土體重力[20],體現在形變圖中,強降雨會加速地表下沉進程。圖6中C2、B4、A5、C6測點的形變狀態對降雨量大小的響應相對及時。以B4測點為例,該點位于王家嶺煤礦2019年采空區附近,因此采礦工程擾動對地表形變影響較小,在2021-05—2021-07降雨期間,其地表累積抬升7.32 mm,在2021-11—2022-02期間,地表累積下沉21.07 mm,B4測點在監測時間內地表形變趨勢以下降為主,但從折線圖可以看出降雨對該點的形變發展起到了一定的控制所用。

1號研究區主體為南茆露天煤礦,A1、B1測點分別在兩個黃土開挖邊坡工作面附近,C1測點位于煤礦排土場東側,排土場是露天煤礦開采過程中產生的廢棄土石的堆放場所,成份主要由黃土、泥巖、砂質泥巖、粉砂質泥巖等組成,其中還混有前期開采排棄的碎石,砂質土[21-22]。雖然排土場經過一定程度的壓實處理,但結構仍然為松散多孔,整體強度較低,在豐沛降雨時間段內,大量雨水會入滲排土場中,導致黃土等巖土體的抗剪強度降低,密度變大,內部的力學性質發生改變,致使其整體有一定量的下沉。而在集中暴雨情況下,降水會在排土場的表面和坡面形成徑流,在雨水的侵蝕和沖刷作用下,排土場局部會發生塌陷,在極端降雨條件下,排土場邊坡可能會發生失穩,滑坡等災害。從C1測點的形變折線圖分析,在2021-08—2021-09時間內,累積下沉達11.35 mm,2022-06—2022-09期間,累積下沉17.78 mm,總體來看,降雨加速了該點的沉降速率。A1、B1測點整體形變發展受降雨影響較小,但在少量降雨的月份,地表會發生小幅度抬升,無雨的時間段內也發生了一定的季節性沉降。需要說明的是,由于露天開采礦區部分區域有較大程度的形變(圖11),因此在數據處理過程中會出現形變失相干的現象,導致這部分區域形變信息的缺失,如大規模巖土體的轉移,局部發生嚴重的塌陷,滑坡等,但通過分析已知測點的形變發展,仍然可以得出降雨量作為誘發地面沉降的重要自然因素之一,與礦區的地表形變有很大程度的相關性。

圖11 1號重點監測區衛星影像Fig.11 Satellite image of No. 1 key monitoring area

4 結論

利用SBAS技術對覆蓋研究區的32景Sentinel-1A影像數據進行處理,獲取保德礦區2020—2022期間形變信息,分析了礦區的形變特征以及引起形變的驅動性因素,得出如下主要結論。

(1)研究區存在6個較大的形變區域,區內地表形變以下沉為主且局部形成多處由中心向四周發展的漏斗狀沉降中心,主形變區域監測點累積形變量與時間基本呈線性關系。

(2)監測時間段內,區域的平均地表形變速率分布在-44～24 mm/a,最大累積沉降為85 mm,最大累積抬升為38 mm。

(3)采礦工程是研究區發生地表形變的主要驅動因素,沉降中心空間分布范圍與各煤礦工作面位置有良好的一致性,斷層分布和地表可壓縮層厚度差異對研究區整體的形變起到了一定的控制作用。從監測點累積形變量整體來看,大多數情況下降水量的增加會減緩礦區地表形變的發展,而在極端降雨條件和某些特殊區域中,降雨則會加快地表形變速率,甚至會引發崩塌,滑坡等地質災害。

(4)通過對形變信息的分析,明確研究區內地表形變誘因,在后續工作中以上述驅動因素為評價指標,結合研究區內人口密度、建筑密度、道路密度、人類工程活動等數據對礦區地質災害發生的風險性以及易損性做出合理的評價,為煤礦安全生產和形變區的綜合治理提供了科學依據,但SBAS-InSAR方法對礦區最大形變量探測能力有限,需要結合其他勘測方法才能更加準確、全面地反映礦區實際情況。