基于SBAS-InSAR技術的M礦區沉降監測與建筑群安全性研究

夏 銳，李 鐵，張景發，田云鋒

（1.北京科技大學土木與資源工程學院，北京 100083；2.應急管理部國家自然災害防治研究院，北京 100085）

引言

在煤炭開采的過程中，由于地下煤層資源被采出，使其周圍巖體的應力平衡被打破，而當新的應力平衡形成時，就會導致巖體發生位移和形變，這種現象稱為“開采沉陷”［1］。礦井開采沉陷會破壞耕地、損壞房屋建筑及地面設施等［2］，給社會經濟的發展和居民的正常生產生活帶來了不可估量的損失［3］，探尋礦區地表的形變監測及其對地面建設的安全影響迫在眉睫。

傳統的地面沉降測量方法有水準測量、GNSS測量等。水準測量運用儀器設備，在野外選取水準點進行實地測量，實操方便，但耗費太多人力物力、監測范圍小［4－5］。GNSS測量運用衛星定位等，監測精度較高，但只是對離散點的監測，無法滿足礦區大范圍沉降區域監測的需要［6］。而近些年興起的合成孔徑雷達干涉測量技術（InSAR技術）能夠實現全天候、全天時、覆蓋廣、精度高的對地觀測，在地表形變的監測研究中得到了廣泛地應用和發展［7－8］。

早在1989年，Gabriel等［9］就基于Seasat衛星觀測數據首次利用D-InSAR技術獲取了農田灌溉區厘米級的地表形變，但對這門新興技術當時關注度還不大。直到1993年，Massonnet等［10］利用D-InSAR技術獲取了1992年Landers地區的地震形變場，并用模型驗證較好，引起了地震界的廣泛關注。后續學者對D-InSAR技術展開了大量研究，將其應用于各個相關的行業，并在此基礎上，不斷改進著D-InSAR技術存在的軌道誤差、時空失相干和大氣延遲誤差等局限性，發展起來了PS-InSAR技術和SBAS-InSAR技術等。PS-InSAR技術也叫永久散射體干涉技術，是一種單主影像時間序列技術，它基于具有穩定散射特性的PS點對研究區進行時間序列分析，可減弱失相干和大氣延遲因素對于干涉質量的影響［11－12］。但該技術對遙感影像數據量的要求比較高，處理過程耗時較長，而且需要研究區具有一定數量的PS點，也給后期解算帶來一定不便［13－14］。SBAS-InSAR技術也叫短基線集技術，它在PS-InSAR基礎上進行了改進，可選擇多幅影像作為主影像，利用較短時空基線生成的干涉對獲取差分相位信息，通過一定的模型算法來減弱誤差值對干涉結果的影響，最終獲得時間序列形變信息［15－16］。眾多學者基于SBAS技術展開對礦區沉降的大量研究［17－20］，也充分證明SBAS技術在礦區沉降監測研究的適用性。

文中利用SBAS-InSAR技術開展M礦區的沉降監測，并依據沉降結果開展礦區地表建筑群的安全性分析。通過處理2017年3月12日至2021年09月11日覆蓋M礦區的121景哨兵數據，確定了礦區近五年連續沉降區域以及年平均沉降速率，并對重點沉降區域選取監測點進行時序分析。另結合哨兵強度圖像及GF-1號影像解譯礦區地表建筑信息，聯合InSAR形變場對礦區地表主要建筑區的安全性進行分析，為地表建筑物的安全防護及安全預防提供指導性意見。

1 研究區概況

本次研究的礦區位于遼寧省沈陽市南20 km，41.45°～41.65°N，123.16°～123.39°E，面積約為750 km2，如圖1所示。自1970年起，礦區上先后新建了三個大型礦井，一礦、二礦、三礦。一礦礦井始于1970年10月份破土興建，1982年11月25日生產，二礦始建于1969年，三礦礦井于1991年12月份破土興建，1997年5月1日開始試生產，2000年正式投產。到2000年三礦建成并正式投產后，礦區的整體建設規模基本形成。礦區之上，雖然沒有密集的城市建筑、設施，但也零散分布著較多以村為單位聚集的建筑群，自北向南一礦上有紅星村、韓城村、英窩村、林盛村等，二礦上有吉祥村、南紅村、西羊角灣村等，三礦上有北青墩子村、尤家甸子村、楊甸子村、樹碑臺村等。這些村落在礦區上也構成了大大小小的建筑區，生活著較多的居民，矗立著眾多的建筑體。隨著礦井的開采，開采深度最深處已經超過了1 000 m，采空區面積也大幅增加，地表沉陷帶來的問題也日益突出，對地面的建筑區以及居民的生活都將帶來一定的影響。

圖1 礦區范圍圖Fig.1 Mining area map

2 SBAS-InSAR技術原理

短基線集（Small Baseline Subsets，SBAS）是Berardino等［21］于2002年提出的一種新的差分干涉測量技術，主要用于提取低分辨率、大范圍、形變速率緩慢的地表形變信息。SBAS-InSAR技術選用較短的時空基線形成干涉對，在干涉處理中組成干涉對的兩幅影像能夠保持較強的相干性，保證了相位處理結果的穩定性。SBAS-InSAR技術對低相干地區的形變仍有較強的適應性，因此也適合用于礦區形變的監測。

假設在t0～tN時間段內共有N+1影像，這些影像自由組成干涉對，可生成M幅干涉圖，且有：

將每組干涉對中兩幅影像的像素值進行差分，理論上將第i幅差分干涉圖的干涉相位可以表示為：

式中：φ(t a,r,x)和φ(t b,r,x)為組成干涉對的兩幅影像分別在ta、tb時刻像元所對應的相位值，ta＜tb；λ表示雷達波長；(r,x)表示像元點坐標；d(t a,r,x)和d(t b,r,x)分別表示像元在ta、tb時刻，相對于t0時刻發生在雷達視線方向上的地表形變。

在實際計算長期連續變化的形變量時，需求解相鄰時相的相位平均變化速率，可以表示為：

而第i副干涉圖的相位又可表示為：

將上式寫成矩陣形式：

利用奇異值分解法求解上式，可以算得不同時段上的相位平均變化速率，最后對各時段上速率進行積分從而獲取各時間段的形變量。

3 SAR數據與處理

文中采用的雷達影像數據是2017年3月12日至2021年9月11日期間獲取的121景C波段哨兵（Sentinel-1 A/B）升軌數據。數據采用GMTSAR開源軟件（https://topex.ucsd.edu/gmtsar/）進行處理。DEM數據采用90 m空間分辨率的SRTM高程數據。Sentinel-1數據的基本信息見表1。用于SBAS分析的干涉對時空基線組合如圖2所示。

表1 Sentinel-1數據參數Table 1 Sentinel-1 data parameters

3.1 配準和干涉

獲取的121景哨兵數據利用GMTSAR軟件進行處理，選取第一景圖像作為主圖像，基于歐空局提供的精密軌道將其他圖像與主圖像進行精確匹配。為了保證圖像之間的相干性，選擇時間基線為120 d，共生成了1 096副干涉對（圖2）。基于DEM生成地形相位，從干涉相對中減去，從而得到主要包含形變信息的干涉相位。

圖2 時空基線Fig.2 Temporal-spatial baselines

3.2 濾波和解纏

干涉相位中含有一定的斑點噪聲，選用高斯濾波器對干涉相位進行濾波處理，去除干涉處理中留下的噪聲信息，使得相位更連續、更加清晰。最后利用Snaphu程序進行相位解纏，恢復相位主值中被模糊的整周相位，從而計算出正確的高程或形變信息。

3.3 小基線集分析

基于上述干涉對進行SBAS分析。即使針對低相干性地表，SBAS也能得到位移結果。研究區內部分地區因水體或植被覆蓋出現了失相干現象，經過統計分析，以0.05～0.1為相干性閾值剔除了此類地區。經過SBAS分析最終得到影像覆蓋區域各象元點的位移時間序列，基于時序可以分析不同時間段內的沉降速率及5年來的平均沉降情況（圖3）。如圖3，從整個圖幅來看，沉降區域呈離散狀分布，最大年平均沉降速率超過190 mm/a。經過與光學遙感影像進行對比，發現沉降地點大多位于礦區，與采礦等人類活動有關，這從一方面也驗證了文中所得InSAR場結果的準確性。在部分地區，InSAR年平均形變場顯示出抬升趨勢，可能與當地的地下水回灌等人類活動有關。

圖3 年平均沉降速率圖Fig.3 Annual average settlement rate map

4 結果分析

4.1 時序監測

圖4為疊加GF-1號背景影像的礦區沉降年平均速率圖，圖中標示了時序監測點的位置。下載的GF-1號影像為2020年4月11日的數據，影像經過正射、融合等處理。從圖4可看出，礦區沉降區域主要集中于三礦和二礦，在一礦也有少量分布。這也與一礦開采年份較早，如今主要集中于二礦、三礦開采有關，近些年三礦的加劇開采導致的沉降最為嚴重，年平均沉降速率最大達93 mm/a。如圖4所示，將沉降區域由南及北依次命名A～G區。A區位于三礦內，包含了南李大人屯村、蒿子屯村、樹碑臺村、楊甸子村等，范圍最大，統計出沉降超過10 mm/a的沉降面積約為15.88 km2，中心沉降速率最大約為93 mm/a，計算沉降體積約190.94萬立方米。B區位于三礦與二礦的交界處，主要分布于二礦內，包含了尤家甸子村、北青堆子村、西大堡村等，沉降中心速率最大約為87 mm/a，但沉降范圍沒有A區大，統計出沉降超過10 mm/a的沉降面積約為4.93 km2，計算沉降體積約73.58萬立方米。圖5為A區、B區的沉降三維示意圖。一礦由于早期已經開采，近些年并未出現較大范圍及較大速率的沉降。一礦沉降區域僅包括若干小的沉降區域，C區～G區，沉降面積都不是很大，分別約為0.59 km2、0.46 km2、0.6 km2、1.15 km2、1.22 km2，中心沉降速率也不是很高，約為20～30 mm/a，合計累計沉降體積約29.43萬m3。

圖4 礦區沉降年平均速率圖Fig.4 Map of the annual average rate of subsidence in the mining area

圖5 A區、B區沉降三維示意圖Fig.5 3D schematic diagram of settlement of zone A and zone B

文中選取沉降區域A區、B區的沉降中心a、b點作為監測點進行時序分析，圖6展示了2個點的時序監測結果。沉降中心點的時序監測結果初步反映了各沉降區域長時間的沉降趨勢。由圖6可看出，自2017年3月12日開始，一直到2021年9月11日，近5年的時間，礦井沉降區A區、B區隨時間推移沉降量一直在增大。且在2018年12月份左右，A區、B區均出現沉降速率先后加大的跡象。B區在2020年7月底至2021年5月初，沉降有一定放緩趨勢，但之后又繼續保持增大。至2021年9月份，A區沉降中心區域點的累計沉降量達384 mm，而B區沉降中心區域點的累計沉降量達351 mm。根據GF-1號影像查找沉降中心位置，目前A區、B區沉降中心地表位置均是農田，并無直接地表建筑物，但按照目前沉降監測點的時序沉降效果來看，A區、B區都有沉降繼續加深、影響范圍繼續擴大的趨勢。

圖6 a、b監測點時間序列圖Fig.6 Time series diagram of monitoring points a and b

4.2 精度驗證

由于本研究缺少礦區的實地水準測量數據，僅獲取了影像范圍內的少量GPS監測資料，以此初步驗證SBAS反演結果的可靠性。在影像區域內總共獲取了2個流動站GPS數據，分別為E317、E320點。其在影像的位置如圖3所示。SBAS監測這兩點的沉降圖如圖7、圖8所示。由于獲取的GPS點監測數據的時間只有2009年到2017年，與SBAS反演所用的時間不同，但鑒于GPS站點具有一定的穩定性，在已知GPS點垂向年沉降速率基礎上，計算其在雷達視線向的年沉降速率，并與SBAS反演結果計算出的年平均沉降速率進行對比，初步驗證結果的可靠性。SBAS反演結果計算的年平均沉降速率即為該站點時間序列上沉降量擬合出來的直線斜率。GPS監測在雷達視線向的年沉降速率值與SBAS結果擬合的年平均沉降速率值對比如表2所示。可以看到，E317點的誤差在1.52～3.78 mm/a，E320點的誤差在1.24～3.08 mm/a。由于InSAR形變監測過程中本身就存在一定的系統誤差，如噪聲、大氣延遲等都會帶來一定誤差，相對于礦區大量級的形變來說，這個誤差是相對較小的，不會對監測結果及建筑區安全評定產生較大的影響，因此這個誤差結果是可以允許的，也表明SBAS-InSAR反演的結果較為可靠。

圖7 E317點沉降圖Fig.7 E317 point settlement chart

圖8 E320點沉降圖Fig.8 E320 point settlement chart

表2 SBAS-InSAR監測結果與GPS測量結果對比Table 2 Comparison of SBAS-InSAR monitoring results and GPS measurement results

5 建筑區安全性研究

礦井開采沉陷會對地面環境與建筑房屋等造成破壞，造成嚴重的安全問題。文中通過覆蓋礦區的GF-1號影像，解譯出礦區地面的建筑物區域，結合礦區的沉降結果，分析礦區沉降對地面建筑物的影響。并通過處理初時刻2017年3月12日與末時刻2021年9月11日的哨兵影像數據，獲取兩個時刻的強度圖，根據強度數據提取地面建筑物信息，判斷地面建筑群在近5年期間是否發生大的變化，以免對評定結果產生影響。

圖9、圖10為礦區范圍內建筑群的解譯結果，結果顯示近5年礦區上的建筑群分布格局及位置基本上是保持不變，在各區域內部有的較少建筑物數目的增加，對本次安全性評價影響極小。通過強度圖像的解譯初步確定了重點研究的建筑群。

圖9 2017.03.12哨兵強度數據建筑解譯Fig.9 Architectural interpretation of sentinel strength data for March 12，2017

圖10 2021.09.11哨兵強度數據建筑解譯Fig.10 Architectural interpretation of sentinel strength data for September 11，2021

文中參考《民用建筑可靠性鑒定標準》（GB 50292－2015）、《建筑地基基礎設計規范》（GB 50007－2011）等眾多標準規范，在鑒定建筑區沉降安全性時，判定標準主要依據是：建筑房屋處于自然狀態時，連續2個月地基沉降量是否大于每月2 mm［22－23］。文中將年沉降速率不超過24 mm/a的礦區地表建筑區定義為相對安全，而超過這個沉降速率值的區域定義為安全性較差，需采取相應措施。

圖11為通過GF-1號數據精細解譯建筑信息并疊加了礦區沉降的示意圖，由圖11可看出，礦區范圍內存在著較多的地面建筑區，但需重點關注的主要是位于沉降區A區、B區附近的建筑群。礦區建筑沉降細節圖如圖12所示。根據建筑群的沉降情況分別設定不同的安全程度等級，如表3依次劃分為：相對安全、值得關注、現場勘探。

圖11 礦區GF-1號建筑解譯及沉降圖Fig.11 Mining building interpretation and settlement map based on Gaofen No.1

表3 建筑區安全分級Table 3 Safety classification of construction area

圖12 建筑區細節圖Fig.12 Detail map of building area

將A區沉降范圍內的建筑群分別命名為A1、A2、A3、A4，B區沉降范圍內的建筑群分別命名為B1、B2、B3。建筑區A1、A2分別位于A區沉降中心的北側和西南側，A1南側處于礦區以內的一小部分位于南李大人屯村，該區已經產生了一定沉降，A2區位于蒿子屯村，雖然并未處于沉降區域正中心，但較大部分已經產生了沉降，這2個地區處于礦區部分沉降速率都在30～60 mm/a之間，遠遠超過了24 mm/a的安全沉降標準，后續需前往現場進行建筑房屋的勘探，查驗房屋變形及損壞情況，確定安全防護措施。建筑區A3、A4處于沉降區域的邊沿，分別位于樹碑臺村、楊甸子村，沉降速率在0～30 mm/a，絕大部分處于相對安全狀態，但小部分已經超過了安全沉降標準，因此值得后期關注。在B區沉降范圍的西北側有建筑群B1和B2，同樣處于沉降區域的邊沿，分別位于北青堆子村、尤家甸子村，B1沉降速率在0～20 mm/a，相對處于安全狀態，而B2沉降速率在0～30 mm/a，且位于A區、B區的中間地帶，隨著A區、B區的擴張，沉降速率可能會繼續增大，也值得后期關注。在B區的南側，有著較大的建筑群B3，位于西大堡村，B區的沉降中心區域剛好緊貼B3建筑區的邊緣，B3范圍內已經發生了較大沉降，沉降速率在30～60 mm/a，需要進行現場勘探。而位于一礦的沉降區C～G區，基本避開了建筑群，而且一礦沉降量較小，幾乎對建筑區域不會產生較大的影響，一礦內的建筑群相對來說較為安全。

綜上，沉降區中的建筑群B1以及C～～G區內的建筑均為相對安全，而A3、A4、B2建筑群，因部分區域沉降速率超過了安全標準，且隨A、B沉降區擴大區域沉降速率可能繼續增大，因此后期值得關注。而A1建筑群南部以及A2、B3建筑群，沉降速率達到了30～60 mm/a，已經極大地超過了建筑安全沉降標準，后續將針對這些區域進行現場勘探，確認房屋變形及損壞情況，并根據實際情況制定安全防護措施。

6 結論

文中基于2017.3～2021.9獲取的121景哨兵數據，利用SBAS-InSAR技術獲取了M礦區近5年的沉降情況，確定了沉降區域及年平均沉降速率，針對重點區域開展了沉降時序分析。結合哨兵數據強度圖像及GF-1光學影像對礦區內的建筑群進行了精細解譯，聯合InSAR形變場對礦區地表主要建筑區的安全性進行分析。

結果表明，礦區以三礦和二礦的地表沉降最為嚴重，一礦的沉降相對較小。這也與礦井的開采歷史及先后開采順序有關。礦井以三礦的A區和二礦的B區沉降最為嚴重。位于三礦和二礦的A區、B區沉降中心累計沉降量分別達384 mm、351 mm，平均年沉降速率分別為93 mm/a、87 mm/a，統計沉降面積分別達15.88 km2、4.93 km2，計算沉降體積分別約為190.94萬m3、73.58萬m3。而一礦有多個沉降區域，但沉降范圍及沉降量相對來說最小，中心沉降速率在20～30 mm/a之間，沉降面積在0.59～1.22 km2之間，累計沉降體積約為29.43萬m3。另結合影像覆蓋區域的GPS監測資料，驗證了SBAS-InSAR形變監測結果的可靠性。礦區地表建筑群中，A1南部進入礦區的部分、A2以及B3建筑群安全性較差，分別位于南李大人屯村、蒿子屯村、西大堡村，沉降速率已經在30～60 mm/a之間，后續還需進行現場勘探，進一步做好建筑物安全防護及預防工作。而A3、A4、B2建筑群分別位于樹碑臺村、楊甸子村、尤家甸子村，沉降速率在0～30 mm/a之間，小部分區域已經超過了安全沉降標準，而且這些建筑群也處于A區、B區之間，隨礦井開采推進沉降速率也會進一步加大，值得后期關注。B1及C～G區建筑群沉降速率為0～20 mm/a，安全性相對較好。

以上研究，一方面給礦區的沉降治理提供參考依據，便于礦區今后進行實地考察，做好相關監測以及安全預防工作。另一方面，根據沉降結果對礦區地表建筑群進行的安全分析，有利于礦區今后合理的開采布置，為礦區地表建筑物的安全防護及安全預防提供指導性意見。