InSAR技術在高寒礦區地災治理監測中的應用

吳啟紅 謝飛鴻 楊何 董建輝 魏占璽 毋遠召 黃程

摘 要：采用InSAR技術對大通煤礦沉陷區進行監測，得到了該礦區2014年10月至2021年11月的地面沉降變形情況和年均沉降速率等信息.結果表明，礦區主要沉降變形區集中在東北部和東南部2個片區，其沉降速率最大超過220 mm/yr；2017年后礦區沉降變形范圍和變形量值都呈顯著下降趨勢且沉降變形逐漸平穩，前期治理工程對控制礦區沉降變形起到了重要作用，可為高寒礦區防災減災及治理工程評價提供重要參考.

關鍵詞：InSAR；礦區治理；地面沉降；沉降速率

中圖分類號：P694

文獻標志碼：A

0 引 言

礦區地質災害主要包括崩塌、滑坡、泥石流、地面塌陷、地裂縫、礦井突水和矸石堆災害等；礦區環境災害主要包括水污染、土污染、水土流失、石漠沙化、地貌景觀破壞和植被破壞等.據不完全統計，在1970—2005年，各類礦山地質災害共計達到2 031次，而在2001—2010年僅10年時間就發生了礦山次生災害12 367起，其中尤以地面塌陷、滑坡和崩塌居多［1］.礦區內剝離大面積土壤植被及形成的煤矸石山，產生大量“黑色”污染和“灰色”污染，形成潛在不穩定斜坡，對礦區的生產設施和施工人員造成嚴重威脅.而高寒礦區由于高海拔、低氧和低溫，其地質環境尤為復雜及脆弱，如玉龍煤礦、木里煤礦和紅河溝煤礦等.因此，開展高寒礦區的變形監測對保障礦區安全至關重要.

傳統的礦區沉降監測方法主要采用精密水準測量和全球定位系統（GPS）測量等.這些監測方法存在監測工作量大、費時、費財和測點難以保存等缺陷，同時其測量精度并不高，從而導致從其測量數據獲得的結果會出現一定的偏差，對后續工作造成影響［2］.合成孔徑雷達干涉測量（InSAR）因在地表變形監測上具有全天時、全天候、高分辨率、高覆蓋率、直觀和形象等其他測量技術不具備的優點，可以在部分領域代替傳統的測量方式［3-4］.同時為去除大氣誤差的影響，加之合成孔徑雷達（SAR）影像數據的增加，在原有InSAR方法的基礎上，逐漸發展出Stacking-InSAR、D-InSAR、PS-InSAR和SBAS-InSAR等高級時序InSAR處理方法［5-6］.

InSAR技術最先被用于地球表面大范圍形變場的觀測，如地震位移測量和冰川漂移等.隨著理論發展與數據處理技術的進步，InSAR技術逐漸在國內外各行各業中得到廣泛應用，如煤田、石油、天然氣開產礦區以及地下水開采區地面沉降監測［7］.高海英等［8］通過小基線集InSAR技術識別黔西南州貞豐縣和安龍礦山地表形變；丁劉建等［9］基于SBAS-InSAR對滕州市附近礦區地面沉降進行監測；王志勇等［10］基于InSAR對濟寧礦區沉降進行精細化監測.這些研究主要針對礦區生產運營階段進行變形監測，而對礦區治理過程及后期的變形監測及應用研究不足.InSAR技術無論是在礦區生產運營階段、閉礦階段還是在治理過程中對災害的預防及評價方面都具有重要意義，同時對于指導我國在類似地區的防災減災工作也具有重要參考價值.因此，為了研究InSAR技術在高寒礦區治理中應用的有效性和實用性，本研究選取青海省大通煤礦區作為研究對象，采用InSAR技術對大通沉陷區進行處理，得到了該礦區治理過程及治理后的地面沉降變形情況和年均沉降速率等信息，分析變形規律及評價治理工程的有效性.

1 研究區概況及數據

1.1 研究區位置

研究區位于青海省西寧市大通回族土族自治縣，地處青藏高原和黃土高原的過渡地帶，地理位置東經101°37′～101°56′，北緯36°55′～36°58′之間，海拔2 450～2 850 m，距省會西寧市35 km.

1.2 治理工程情況

2012—2015年共分4期對沉陷區內涉及橋頭鎮及良教鄉的7個村（橋頭鎮的元樹爾、小煤洞和大煤洞3個村，以及良教鄉的下甘溝、煤洞溝、白崖和上甘溝4個村）范圍內因煤炭開采造成地表土地資源壓占、挖損、植被破壞、粉塵污染和地下含水層破壞等，以及后期剩余煤炭資源開采可能引發的采空塌陷范圍進行治理，治理建設規模為1 385.02 hm2.雖然經過4期地質環境治理，但區內依然存在地質災害和生態環境破壞等問題，因此，2020—2022年，又對歷史遺留和新產生的塌陷區進行了綜合整治.

1.3 數據情況

1.3.1 Sentinel-1數據及獲取

Sentinel-1衛星數據是歐洲航天局為接替Envisat-ASAR衛星發射的科研SAR衛星，數據免費開放.該衛星成像質量高，軌道穩定，有多極化能力，單顆衛星重訪周期24 d，繼Sentinel-1A衛星后，Sentinel-1B衛星的發射使常用重訪周期變為12 d，遇到突發情況時可以2顆衛星合作將觀測重訪周期縮短至6 d.查詢可知，該地區Sentinel-1數據覆蓋良好，數據日期自2014年10月至今.Sentinel-1數據參數見表1.

Sentinel-1衛星執行全球對地觀測極化，以干涉寬幅成像為主，該模式有益于InSAR的干涉處理，且250 km×250 km的幅寬保證單景數據覆蓋觀測區.該衛星工作波段為C波段，雖然對植被的穿透能力較弱，但衛星數據量多、成像質量高，可以一定程度上彌補數據穿透能力弱的問題.本次計算使用2014年10月7日至2021年11月23日間共計173期降軌Sentinel-1數據，相對軌道號33，數據重訪周期為12/24 d.

1.3.2 DEM數據

本次計算使用ALOS World 3D DEM，其空間分辨率為30 m，該數據集是高精度全球數字地表模型數據，水平分辨率為 30 m，高程精度5 m，由高級陸地觀測衛星ALOS上搭載的全色遙感立體測繪儀PRISM獲取.是目前世界上最精確的3D地圖之一，覆蓋全球所有的土地尺度.

2 InSAR數據處理

2.1 InSAR技術原理

InSAR技術主要基于雷達影像的相位信息獲取目標點至雷達傳感器之間的距離，通過獲取2期甚至多期次影像并對其進行干涉即可丈量目標點在一定時間段內的位移變化量，其屬于定量遙感技術范疇［11-12］.目前，常用星載SAR平臺獲取的SAR影像進行差分干涉處理，其具有覆蓋面積大、飛行軌道穩定和獲取數據成本較低等優勢.星載SAR平臺的基本工作模式為重軌飛行，即1顆雷達衛星在相鄰軌道上對同一地區2次SAR成像，每幅SAR圖像中的每個像素均記錄了強度值和絕對相位值.針對同一地區獲取的2幅SAR影像進行差分干涉處理，可得到該地區地表形變信息.

2.2 數據處理

2.2.1 SAR數據配準

在進行干涉計算前，首先要對生成的單視復數數據（SLC）在外部DEM數據的參與下進行數據配準.Sentinel-1數據成像模式較為特殊，為多個子帶拼接構成的，為保證子帶間干涉處理時不發生相位的躍變，SLC數據配準精度要求較高，尤其是針對Sentinel-1數據而言，距離向（Azimuth）配準精度需要達到1/1 000像元級別.

2.2.2 差分干涉處理

如圖1所示，選擇2018年1月19日SAR影像為主影像進行圖像配準計算.其中，Azimuth方向配準精度小于1/1 000像元，滿足Sentinel-1數據配準要求.干涉像對最小時間間隔為12 d，最大時間間隔為120 d，最小垂直基線長度為0.3 m，最大垂直基線長度為216.8 m，完全滿足干涉計算需求，共有干涉像對513個.

2.2.3 自適應濾波處理

為平滑干涉圖噪聲，使用窗口為32的自適應濾波進行處理，增強干涉圖低質量區的相干性，方便下一步相位解纏處理.

2.2.4 相位解纏

使用最小費用流方法進行相位解纏處理，相干性閾值設置為0.2，解纏相位均過度平滑，未出現相位跳躍的情況，但部分解纏相位圖受大氣效應干擾嚴重.

2.2.5 大氣相位估計

估算高程相關的大氣延遲相位，并利用快速空間濾波的方式初步估算隨機大氣相位，而后去除大氣延遲相位，大氣相位干擾情況將被明顯改善.

2.2.6 Stacking-InSAR和SBAS-InSAR計算

利用去除大氣相位后的計算結果進行2014—2021年度的Stacking-InSAR計算，估算年均變形速率，并分析工作區逐年沉降變形情況.同時通過SBAS-InSAR方法計算典型形變點的累計變形曲線.

2.3 誤差源分析

InSAR計算中，會受到多種因素的影響，而且在不同地區，使用不同類型數據，同樣的影響因素造成的誤差程度不同，在計算過程中要結合實際情況加以分析，針對不同地區和不同數據來調節參數.誤差源主要有系統誤差、DEM誤差、失相干和大氣誤差.

3 InSAR結果分析

受制于Sentinel-1衛星數據本身分辨率，工作區InSAR觀測結果空間分辨率為近15 m.工作區觀測到的變形主要是采礦后沉降變形，所以將InSAR計算得到的雷達視線向（LOS）變形轉為垂直向變形.

3.1 2014—2021年逐年變化分析

利用Stacking-InSAR方法計算工作區2014—2021年間逐年地面沉降變形情況，2014年的InSAR數據是從2014年10月開始，經過運算后，圖中代表的是整個自然年的變形分布，如圖2所示.

由圖2可知，大通礦區地面沉降情況自2014—2021年呈現在變形范圍和沉降量值上均表現為減弱的趨勢.礦區的沉降區主要集中在東北部和東南部，可以分為2個大的沉降區，其他零星分布.沉降時間主要發生在2016年前，在2014—2015年間年均沉降速率大，最大可達220 mm/yr以上，自2017年后大通礦區沉降變形呈現為明顯減弱趨勢，2018年后明顯沉降區主要集中在礦區東南角.

3.2 2014—2021年多年期時間序列分析

采用SBAS-InSAR方法針對大通礦區獲取其2014—2021年的形變情況，SBAS-InSAR年均形變速率結果如圖3所示.

由圖3可知，大通礦區主要沉降區可分為1#沉降區和2#沉降區，整個礦區面上為零星分布點狀變形，高變形區主要分布在礦區東南角，未出現大面積高速沉降區，礦區2014—2021年間年均變形速率約為58 mm/yr.選取礦區6處典型變形點做時間序列分析，其累計變形曲線如圖4所示.SBAS1和SBAS2位于礦區東北部的1#沉降區.SBAS1位置累計變形量約340 mm，SBAS2位置累計變形量約320 mm，2條曲線都表明在2014—2017年初變形速率較大，后變形速率逐漸平穩.SBAS3和SBAS4位于礦區東南部的2#沉降區.SBAS3位置累計變形量約270 mm，曲線表明2014—2017年初變形速率較大，后變形速率逐漸平穩，到2020年變形速率又開始加大；SBAS4位置累計變形量已達400 mm，曲線表明該點2014—2017年初變形速率較大，后逐漸平穩，至2021年4月后短暫出現加速變形情況，后逐漸平穩.SBAS5位置累積變形量約160 mm，曲線表明2014年6月至2017年6月間變形速率較大，后逐漸平穩，并在2021年1月至2021年6月間出現加速變形情況.SBAS6位置累積變形量約40 mm，曲線表明該點2014—2019年變形較為平穩，而2020年開始出現加速變形的情況，但形變量值較小.

在礦區內做5條年均沉降速率剖面，剖面線位置如圖3所示，剖面結果如圖5所示.剖面A-A′表明沿該剖面線方向，地面沉降速率逐漸變大，在剖面線上1 200 m位置處出現轉折，變形速率開始減小；剖面B-B′表明沿該剖面方向變形速率較為平穩，局部出現較大沉降變形區；剖面C-C′表明沿該剖面線方向約1 300 m前變形速率平穩，1 300 m后出現較大速率變形區，在1 700m位置處變形速率最大；剖面D-D′表明沿該剖面線方向約1 000 m前變形速率平穩，1 000 m后出現較大速率變形區；剖面E-E′表明沿該剖面線方向整體變形平穩，局部單點存在較大變形，所以剖面呈離散狀.

3.3 2020年11月至2021年11月地面變形分析

2014—2021年間時間跨度較長，大通礦區地面變形在長時間跨度下可能會出現階段性的大幅度變化，所以使用2020年11月至2021年11月間33期影像對工作區做1年期的SBAS-InSAR沉降變形速率觀測，其結果如圖6所示.

由圖6可知，大通礦區2020年11月至2021年11月間最大年均變形速率不超過52 mm/yr，高變形區主要分布在礦區東南角，整個礦區面上為零星分布點狀變形，未出現大面積高速沉降區.同樣選取礦區6處典型變形點做時間序列分析，其累計變形曲線如圖7所示.SBAS1與SBAS2位置累積變形量較小，數據較離散，變形較為平穩，與圖4中趨勢一致；SBAS3和SBAS4位于礦區東南部，其變形曲線顯示該區域于2021年5月后開始加速變形，加速變形于2021年10月結束，后變形較為平穩；SBAS5點位于礦區中部位置，其變形曲線表明該位置于2021年5月后開始加速變形，并于2021年7月后加速變形停止，后變形逐漸平穩；SBAS6點位于礦區西部，其變形曲線表明該點一直處于勻速變形狀態.

礦區內沿坡向做5條年均變形速率剖面，剖面線位置如圖6所示，剖面線結果如圖8所示.與圖5年均速率變形分布圖一致，沿剖面線方向沉降變形表現為零星點狀分布，并未連接成為整體，所以剖面線方向上速率點較為離散.

3.4 對比分析

根據前期勘察設計資料和現場調查，2個集中沉降區主要位于煤礦采空區及影響區內.2012—2015年對礦山進行4期次的地質環境治理示范工程，主要治理措施包括煤矸石棄渣處置、廢棄礦井井口封堵、塌陷區土地整治、泵站及田間灌溉配套、溝道防洪、小煤洞渣堆景觀區、引水蓄水灌溉系統、排水系統、不穩定斜坡治理、田間道路、林網及礦山公園等工程.經過治理，沉降區范圍及沉降量逐漸減小，特別是2017年沉降基本趨于穩定，可見治理工程對控制沉降起到了重要作用，治理效果明顯.2020年礦區內部分沉降變形繼續增加，但變形速率較之前減小較多，減小約3/4.根據圖6和圖7變形特征圖，結合2021—2022年度綜合治理實施方案（見圖9）進行對比，可見實施方案中12個不穩定斜坡、3個滑坡和5個沉陷區在InSAR解譯圖中都有對應的變形解譯特征點.

4 結 論

1）利用SBAS-InSAR觀測大通礦區2014—2021年間長時間序列變化，表明礦區主要沉降變形區集中在東北部和東南部2個片區，大通礦區2014—2016年間沉降速率最大超過220 mm/yr，自2017年后礦區沉降變形范圍和變形量值都呈顯著下降趨勢且沉降變形逐漸平穩，但是存在局部或某時間段內發生突然變形加速的現象.

2）2014—2021年的5條沉降速率剖面表明主體沉降位于礦區中前部，后部（西部）沉降速率很小；2020年11月至2021年11月間的InSAR結果表明，高變形區主要分布在礦區東南角，整個礦區面上為零星分布點狀變形，未出現大面積高速沉降區.

3）通過對比分析，表明前期治理工程對控制礦區沉降變形起到了重要作用，近期變形有所增加，但變形速率減小較大；可見InSAR技術對于隱患點的識別及變形監測具有大范圍和高精度的特點，可以為防災減災及治理工程評價提供依據.

參考文獻：

［1］陳家彪，沈冰，劉飛燕.西南地區礦山環境地質問題研究［J］.礦產綜合利用，2007，146（4）：43-46.

［２］劉國林，張連蓬，成樞，等．合成孔徑雷達干涉測量與全球定位系統數據融合監測礦區地表沉降的可行性分析［J］.測繪通報，2005，51（11）：13-16．

［3］蒲川豪，許強，趙寬耀，等.利用小基線集InSAR技術的延安新區地面抬升監測與分析［J］.武漢大學學報（信息科學版），2021，46（7）：983-993.

［4］陸會燕，李為樂，許強，等.光學遙感與InSAR結合的金沙江白格滑坡上下游滑坡隱患早期識別［J］.武漢大學學報（信息科學版），2019，44（9）：1342-1354.

［5］李承航，張文春.基于InSAR技術的監測應用及研究進展［J］.北方建筑，2020，5（4）：15-19.

［6］范洪冬，鄧喀中，承達瑜.InSAR拓展和融合技術在礦山開采監測中的應用［J］.金屬礦山，2008，382（4）：7-10.

［7］張訓虎，章磊，郝樹賓，等.合成孔徑雷達干涉（InSAR）測量技術應用及展望［J］.北京測繪，2014，115（2）：28-31.

［8］高海英，趙爭，章彭.時序InSAR的貴州地質災害監測［J］.測繪科學，2020，45（7）：91-99.

［9］丁劉建，陶秋香，高騰飛，等.SBAS InSAR技術在礦區地面沉降監測中的應用［J］.中國科技論文，2019，14（3）：320-325.

［10］王志勇，張繼賢，黃國滿.基于InSAR的濟寧礦區沉降精細化監測與分析［J］.中國礦業大學學報，2014，43（1）：169-174.

［11］Bamler R，Hartl P.Synthetic aperature radar interferometry［J］.Inverse Probl，1998，14（4）：1-54.

［12］Rosen P A，Hensley S，Joughin I R，et al.Synthetic aperture radar interferometry［J］.P IEEE，2000，88（3）：333-382.

（實習編輯：姚運秀）

Abstract：

InSAR technology is used to monitor the subsidence area of Datong coal mine，and the information of land subsidence deformation and average annual subsidence rate in the mining area from October 2014 to November 2021 are obtained.The results show that the main subsidence deformation areas in the mining area are concentrated in the northeast and southeast，and the maximum subsidence rate is more than 220 mm / yr.After 2017，the settlement deformation range and deformation value of the mining area has shown a significant downward trend，and the settlement deformation has been gradually stable.The preliminary treatment work has played an important role in controlling the settlement and deformation of the mining area.The study can provide an important reference for disaster prevention and reduction and treatment engineering evaluation in alpine mining areas.

Key words：

InSAR;mining area treatment;land subsidence;sedimentation rate