基于InSAR技術的大同市云岡礦區地面沉降監測

黃德華

(福建省地質測繪院,福建 福州 350011)

開采地下礦產資源容易破壞巖層結構,可能引發地面沉降等一系列地質問題,從而對沉陷區的地表建筑、公共設施、耕地、水環境等人類生存環境、經濟財產安全造成極大的影響[1-4]。為了區域的可持續發展,科學合理地進行開發、生態治理,對礦區長時間序列的動態監測十分必要。

由于礦區形變周期性長、范圍廣,而常規的監測手段具有監測點離散、周期長、成本高的特點,難以進行時空連續性監測。合成孔徑雷達差分干涉測量技術是通過短周期內的影像對干涉相位的差分計算,可大范圍觀測地表形變,其監測成本低,精度可達厘米至毫米級[5]。2000年CARNEC C等[6]首次將D-InSAR技術用于煤礦區地表形變場監測;2005年吳立新等[7]對中國東部典型工礦區進行了礦區地表沉陷D-InSAR監測研究。在D-InSAR技術的基礎上,2001年FERRETTI A等[8]提出了PS-InSAR技術;2002年BERARDINO P等[9]隨后提出了SBAS-InSAR技術,此后更是發展了眾多時序監測算法。近年來,多國學者利用InSAR技術對礦區地面開展了形變監測,例如:2015年GRZOVIC M等[10]使用PS-InSAR與SBAS-InSAR方法對美國伊利諾伊州斯普林菲爾德地下煤礦采區進行了地表形變監測,監測結果顯示地下礦井坍塌導致了該區域的地表形變;2023年虎小強等[11]使用SBAS-InSAR方法對新疆阿希礦區地面沉降時空變化特征進行監測,分析了沉降主要驅動因素,等等?，F有研究表明,SBAS-InSAR與PS-InSAR是目前比較有代表性的時序InSAR算法,可有效監測礦區地面沉陷。

本文闡述了InSAR技術原理,運用D-InSAR、時序InSAR技術對云岡礦區開展地面沉降監測,其中在2019年1月—2020年12月時段采用24景哨兵-1A衛星影像數據進行SBAS-InSAR處理獲取形變信息,在2020年1月—2020年9月時段采用25景Sentinel影像數據進行PS-InSAR處理獲取形變信息,并選取2020年1月—2020年9月兩種技術方法不同形變信息進行對比分析,驗證兩種時序InSAR技術在礦區地表沉降監測應用中的可靠性。

1 研究區及數據概況

1.1 研究區概況

研究區位于山西省北部、大同市西部云岡區內,面積約137.16 km2,西北部為煤礦開采區,79%以上為地下開采區,其余為露天開采區。東南部為大同市城市區域(圖1)。

(a).大同市范圍;(b).研究區范圍圖1 研究區概況圖Fig. 1 The overview of the study area

1.2 實驗數據

本文采用2019年1月—2020年12月的哨兵-1A衛星數據及SRTM 30 m分辨率數字高程數據作為數據支撐,其中SBAS-InSAR和PS-InSAR采用數據情況詳見表1。

表1 哨兵-1A衛星影像主要信息表Table 1 Main information of Sentinel-1A satellite image

2 數據處理原理及流程

InSAR技術中相位信息由平地相位、軌道誤差、地形相位、形變相位、大氣相位以及噪聲相位共6部分相位信息組成[12],公式如下:

φdiff=φtopo+φdef+φorb+φatm+φnoise,

式中:φtopo為參考數字高程模型(Digital Elevationg Model, DEM)殘余地形相位;φdef為形變相位,由線性形變及非線性形變組成;φorb為去平引入的軌道誤差相位,空間上表征為一定的線性特性,可包含于φatm中;φatm為大氣延遲相位;φnoise為噪聲相位。

對形變前后影像進行干涉處理,引入參考DEM模擬地形相位,去平、降噪、減弱大氣延遲后,剔除模擬地形相位,即可獲取形變相位從而生成差分干涉圖,反演一維地表視向形變信息[13]。

相比D-InSAR技術的短期形變場獲取,多時相合成孔徑雷達干涉測量技術(Multi-temporal Synthetic Aperture Radar Interferometry, MT-InSAR)增大影像數量,通過空間域、時間域濾波,緩解大氣延遲效應等誤差,能對長時間序列下的影像集進行綜合分析,得到地表長期演變趨勢[14]。永久散射體技術(PS-InSAR)及小基線集技術(SBAS-InSAR)是2個目前比較有代表性的MT-InSAR算法[15]。

2.1 SBAS-InSAR

SBAS-InSAR技術將所有覆蓋同一地區的SAR影像組成若干個子集,子集內的影像基線距(包括時間基線距和空間基線距)小,子集間的影像基線距大,運用奇異值分解方法獲取最小形變速率,基于最小形變速率標準獲得小基線干涉圖。SBAS-InSAR方法限制了長基線導致的幾何去相干,而且使更多的SAR圖像參與到形變計算,增加了時間上的采樣。提出SBAS方法是為解決影像數量不足、基線失相干無法進行相位分離的問題,基于一定的空間時間基線的若干干涉對集進行后續相位解纏等操作,從而產生解纏相位圖用以分析地表形變情況,以減少時間和空間去相干因素的影響,可以得到監測區域的平均形變速率[16]。

SBAS-InSAR主要處理流程(圖2)包括:估算影像對基線,選擇配準的主影像,完成時序影像配準;按照短空間基線組合原則,生成短基線干涉對數據集,進行常規差分干涉步驟,去除平地、地形相位,獲取差分干涉圖集;利用相干系數圖,選取高相干點,進行相位解纏;通過相干性設定閾值進行分布式目標分析提取;采用奇異值分解法模型求解形變參數和高程誤差值即殘余地形的最小范數解;進行高通、低通濾波進行第二次反演,計算非線性形變和大氣相位并進行分離,從而反演形變總和,將結果進行地理編碼生成時序形變信息[17]。

圖2 SBAS-InSAR方法基本流程圖Fig. 2 Basic flowchart of SBAS-InSAR

2.2 PS-InSAR

PS-InSAR技術選取影像集采用長時間序列下各個分辨單元內部保持穩定、散射特性強烈的永久散射體,利用相位信息可靠的特征,對這些離散點進行時間序列分析,從而提取目標形變量,探測地表形變信息[18]。

PS-InSAR技術(圖3)利用K+1幅SAR單視復數影像集,結合外部DEM,計算得到K幅差分干涉圖,并通過幅度離差系數法等算法提取穩定的N個PS候選點[19]。

圖3 PS-InSAR方法基本流程圖Fig. 3 Basic flowchart of PS-InSAR

根據大氣相位時序隨機、空間相關的特點,構建差分網絡,對相鄰的兩PS點干涉相位進行3次差分,以降低大氣噪聲相位影響。利用二維周期圖估計高程以及形變參數,根據求解結果進行相位解纏,得到的相位殘差中包括殘余大氣信號、非線性形變信號。由于大氣相位空間域上低頻信號、時間域上高頻信號,進行空間域低通濾波、時間域高通濾波,將大氣相位與非線性形變相位分離,將線性形變相位與非線性形變相位相加即為視向形變相位,最終反演形變場[20]。

3 結果分析

3.1 D-InSAR沉降結果分析

對2020年6月2日—2020年7月8日的2景哨兵-1A影像進行D-InSAR處理,生成差分干涉圖(圖4)。結果表明,研究區西北部榮華皂村、栗莊村、永定莊村、里南溝村等地可見明顯沉降漏斗,最大視向沉降量達4.6 cm/36 d,沉降區域普遍形變量為2 ～3 cm/36 d;東南部城市區域地表穩定,無明顯形變。

圖4 D-InSAR監測研究區沉降結果圖Fig. 4 Subsident results of D-InSAR monitoring area

3.2 SBAS-InSAR沉降結果分析

2019年1月—2020年12月的24景哨兵-1A影像進行SBAS-InSAR處理,獲得平均形變速率圖(圖5)。研究區西北部存在12處明顯沉降,分布范圍基本與D-InSAR結果重合,且由外緣至中心形變量逐漸增大,沉降區最大平均沉降速率達180 mm/y,其余區域平均沉降速率為5～0 cm/y,2019—2020年最大累積沉降量為333 mm,存在沉降中心超出形變測量梯度的情況;東南部城市區域地表穩定,無明顯形變。

圖5 SBAS-InSAR技術監測研究區平均形變速率圖Fig. 5 Average deformation rate of SBAS-InSAR monitoring area

根據所屬地理位置將沉降漏斗劃分為榮華皂村沉降區、晉華宮街道沉降區、忻州窯村沉降區、曹家窯村沉降區、里南溝村沉降區以及永定莊村沉降區6個沉降區域,其中里南溝村沉降區沉降現象最為嚴重,沉降區域面積達1.27 km2。

選取里南溝北西向(C—D)、北東向(A—B)剖面做沉降漏斗分析,得到其空間累積沉降分布圖(圖6)。A—B剖面(圖6(b))相對完整,即累積沉降剖面線由北部邊緣-中心-南部邊緣呈淺-深-淺漏斗狀,該剖面沉降中心隨時間變化基本維持不變,各處沉降量隨時間變化而增加,總體不斷下沉;C—D剖面(圖6(c))不完整,邊緣至中心沉降逐步加深,由累積沉降剖面可知,2019年沉降中心尚可監測,隨著時間推移,該剖面沉降中心逐步向東轉移,預計未來沉降中心沉降值將超出探測梯度。

(a).里南溝村沉降區剖面位置圖;(b). A—B剖面示意圖;(c).C—D剖面示意圖圖6 2019—2020年研究區累積沉降剖面示意圖Fig. 6 Schematic diagram of cumulative subsidence profile from 2019 to 2020 in the study area

3.3 PS-InSAR沉降結果分析

對2020年1月—9月的25景哨兵-1A影像進行PS-InSAR處理,獲得平均形變速率圖(圖7)。研究區西北部存在12處明顯沉降,分布范圍基本與SBAS-InSAR結果重合,從形變分布上可以看出,兩者所獲取的形變區域空間分布特征一致,沉降漏斗位置分布吻合,且形變趨勢及累積形變量較為一致。

圖7 PS-InSAR技術監測研究區平均形變速率圖Fig. 7 Average deformation rate of PS-InSAR monitoring area

3.4 SBAS與PS技術結果交叉驗證

在3個明顯沉降區,選取里南溝村沉降區邊緣至沉降中心的特征點進行時序對比分析,采用2020年的SBAS-InSAR監測結果,并提取PS-InSAR中對應時相的監測結果,將2020年1月設為起始月份,繪制特征點2020年1月—2020年9月的累積形變曲線(圖8)。

(a). 里南溝村沉降區;(b). A-1點形變量;(c). A-2點形變量;(d). A-3點形變量圖8 里南溝村沉降區域3個特征點累積形變圖Fig. 8 Cumulative deformation of three characteristic points in area of Linangou Village

由里南溝村沉降區域3個特征點累積形變曲線圖(圖8)可看出,SBAS-InSAR與PS-InSAR監測結果總體趨勢一致,表明2020年以來,該區域持續沉陷,其中A-1特征點SBAS-InSAR技術監測累積形變量最大達35.93 mm,PS-InSAR技術監測累積形變量最大達41.99 mm;A-2特征點SBAS-InSAR技術監測累積形變量最大達49.88 mm,PS-InSAR技術監測累積形變量最大達69.17 mm;A-3特征點SBAS-InSAR技術監測累積形變量最大達57.71 mm,PS-InSAR技術監測累積形變量最大達70.99 mm。用均方根誤差(RMSE)對A-1、A-2、A-3點的形變量觀測值xsbas,t和xps,t進行比較:

計算得A-1點RMSE值為2.89,A-2點RMSE值為7.48,A-3點RMSE值為5.31。

對兩組形變量觀測數據做線性擬合(圖9),A-1觀測點的決定系數R2為0.947 2,A-2觀測點的決定系數R2為0.960 1,A-3觀測點的決定系數R2為0.962 6。

圖9 里南溝村沉降區域3個特征點兩組形變量觀測結果擬合曲線圖Fig. 9 Fitted curve of the two observation groups of shape variables at three characteristic points in area of Linangou Village

由上可知,SBAS-InSAR形變監測與PS-InSAR形變監測結果相近,形變區域吻合,形變趨勢基本一致,決定系數均>0.9,驗證了SBAS與PS兩種時序InSAR技術在區域形變監測具有一致性。

4 結論

大同市云岡區作為全國最大的產煤區之一,其地下煤礦開采區引起地面沉降災害。本文應用InSAR技術對區域的多景哨兵-1A雷達影像進行處理,獲取了地面沉降特征,詳細分析了時空演變規律并進行了初步探討,主要結論如下:

(1) 沉降主要分布于云岡區西北部煤礦開采區,沉降面積大,沉降速率達180 mm/y,最大累積沉降量為333 mm;東南部城市區域沒有明顯形變跡象。

(2) 采用SBAS-InSAR與PS-InSAR技術對其地表沉降情況進行長時間序列監測,監測結果表明形變空間分布吻合,時序沉降趨勢一致,形變速率具有較高的相關性(R2決定系數均>0.9),均能有效監測地表沉降變化。