InSAR技術在高寒礦區(qū)地災治理監(jiān)測中的應用

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

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

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

中圖分類號：P694

文獻標志碼：A

0 引 言

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

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

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

1 研究區(qū)概況及數(shù)據(jù)

1.1 研究區(qū)位置

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

1.2 治理工程情況

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

1.3 數(shù)據(jù)情況

1.3.1 Sentinel-1數(shù)據(jù)及獲取

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

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

1.3.2 DEM數(shù)據(jù)

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

2 InSAR數(shù)據(jù)處理

2.1 InSAR技術原理

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

2.2 數(shù)據(jù)處理

2.2.1 SAR數(shù)據(jù)配準

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

2.2.2 差分干涉處理

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

2.2.3 自適應濾波處理

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

2.2.4 相位解纏

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

2.2.5 大氣相位估計

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

2.2.6 Stacking-InSAR和SBAS-InSAR計算

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

2.3 誤差源分析

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

3 InSAR結果分析

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

3.1 2014—2021年逐年變化分析

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

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

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

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

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

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

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

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

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

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

3.4 對比分析

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

4 結 論

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

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

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

參考文獻：

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

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

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

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

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

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

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

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

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

［10］王志勇，張繼賢，黃國滿.基于InSAR的濟寧礦區(qū)沉降精細化監(jiān)測與分析［J］.中國礦業(yè)大學學報，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