聯合LiDAR DEM與時序SAR數據的露天礦特大型滑坡監測

魏戀歡, 劉善軍, 楊天鴻, 趙 曄

(1. 東北大學 資源與土木工程學院, 遼寧 沈陽 110819; 2. 遼寧省自然資源事務服務中心, 遼寧 沈陽 110032)

撫順西露天礦是亞洲最大的露天煤礦,由于長期高強度開采,曾誘發了90余次滑坡災害,其中規模最大的為南幫的特大型滑坡[1].2010年至今,撫順西露天礦南幫受地質條件、采礦活動、降水等因素的綜合影響,形成了東西長約3.1 km、南北寬1.2～1.5 km、高差約102 m的特大型推移式順層巖質滑坡[1].針對該特大型滑坡,西露天礦已開展了GPS測量、現場邊坡雷達傳感器測量、深部巖移傳感器測量、應力測量等工作,但是由于存在監測點稀疏的缺陷,難以準確評估滑坡體大范圍、長期時空演化特征[2-4].

近年來,隨著合成孔徑雷達(synthetic aperture radar, SAR)衛星數據資源的不斷豐富,合成孔徑雷達干涉測量(interferometric SAR, InSAR)以其不受光照和氣象條件限制、高精度、低成本、廣覆蓋等優勢,在滑坡、地面沉降、地震、火山、開采沉陷、建/構筑物精細變形監測等領域取得了大量的成功應用[5].然而,在長期變形監測應用中,不同時段的位移速率往往差異較大,采用單一的InSAR方法難以滿足長期跟蹤監測的需求,因此需要多方法的協同監測.對于每年位移在幾十cm到十幾m之間的快速變形,可用SAR影像的像素偏移量跟蹤(pixel offset tracking, POT)技術進行監測.POT不受去相關因素的限制,根據不同時相SAR影像之間的精密配準參數,由同名像素點之間的偏移量推算雷達視線向(line of sight, LOS)和方位向(azimuth, Az)的位移參數,精度一般能達到dm級[6-7].針對每年位移在幾mm到幾十cm之間的緩慢變形,可采用多時相InSAR(multi-temporal InSAR，MTI)技術進行監測,理論上在雷達LOS向的測量精度可達到mm級[8].

在滑坡變形監測應用中,邊坡滑動位移通常是沿坡面向下的.受限于雷達衛星固有的側視成像特點,MTI和POT技術測量的僅為LOS或Az-LOS向的位移參數,可看作是邊坡坡向位移的投影分量.如何根據LOS或Az-LOS向位移參數求解真實的坡向位移是提高基于SAR數據的滑坡監測精度的關鍵[9].針對上述問題,本文提出一種聯合高精度激光雷達(light detection and ranging, LiDAR)數字高程模型(digital elevation model, DEM)與時序SAR影像的滑坡體坡向位移求解方法,可實現Az-LOS向變形參數到最大坡度方向的轉換.基于該方法,本文開展了撫順西露天礦南幫特大型滑坡長期變形監測,結合礦區地質資料、GPS監測數據、降水量等,分析了2014—2020年西露天礦南幫特大型滑坡變形的空間分布特征及時間變化規律.

1 方 法

1.1 POT方法及Az-LOS向到坡向變形轉換

傳統的基于光學影像的POT方法受到圖像噪聲、傳感器姿態等因素的限制,精度一般不高,且存在無法估算垂直位移、只能估算水平位移的問題[10-11].基于SAR影像的POT方法不受相位解纏及去相干問題的影響,通過使用歸一化互相關優化窗口,跟蹤不同時相SAR圖像之間后向散射特征點位置的變化(即偏移量),可以同時測量Az-LOS向的二維變形[12-15].由于SAR圖像之間的配準誤差通常為像素分辨率的1/10～1/20,當地表位移量小于配準誤差時,POT將難以提取精確的偏移量,因此POT適合監測快速變形[6].一般來說在其他條件相同時,SAR影像分辨率越高,POT對變形的敏感度及監測精度也相應越高.

為了精確監測露天礦滑坡沿坡面的滑移參數,本文提出一種高精度DEM輔助的Az-LOS向到坡向位移的轉換方法,原理如圖1所示.假設A點為某坡面上的滑坡點,坡向變形dA的水平分量dh=dA·cosβ,垂直分量dv=dA·sinβ,其中β為A點的坡度.通過圖1的幾何投影關系可建立dA與dAz和dLOS之間的關系式,即

dAz=dh·cos(α0-αA)=dA·cosβ·cos(α0-αA) ,

(1)

(2)

式中:dAz為方位向變形量;dLOS為LOS向變形量;β為坡度;α0為航向角(順時針為正)；αA為坡向角;θ為入射角.將式(1),式(2)聯立消去β可得坡向變形量dA的計算公式:

(3)

1.2 MTI方法及LOS到坡向變形轉換

首個MTI方法是意大利米蘭理工大學于2000年提出的永久散射體干涉測量(permanent scatterers InSAR, PSI)方法,它克服了傳統差分干涉測量方法中存在的時間失相干、幾何失相干及大氣延遲的影響,通常可以獲得mm級精度的變形[16-17].在PSI基礎上,雷達遙感領域的學者相繼提出了包括小基線子集法(small baseline subset, SBAS)在內的多種類似方法,統稱MTI方法[18-23].這些方法原理與PSI基本相同,主要區別在于干涉圖組合方式、提點策略、相位解纏算法、是否采用變形模型等方面[8].由于本文主題不是為了比較不同MTI方法的性能,考慮到SBAS法能夠充分利用所有可能的高質量干涉圖，提高觀測時間采樣率，對于時空失相干具有很好的魯棒性,因此本文使用SBAS法進行MTI處理.

針對基于MTI的滑坡監測中面臨的LOS向變形到坡度方向的轉換問題,作者曾提出一種基于LiDAR DEM輔助的LOS向到坡向變形的轉換方法,已在鞍鋼集團各大鐵礦高陡邊坡監測中取得了成功應用[24-25].本文繼續采用這一方法實現LOS向到坡向的轉換.轉換過程中,需要根據LiDAR DEM計算礦坑內精確的坡度、坡向參數,結合雷達衛星數據的成像幾何求解坡向變形.

2 實驗區及數據

2.1 實驗區工程地質條件

撫順西露天煤礦位于遼寧省撫順市西南部,最早開采于1901年,1914年開始露天開采,至今形成了東西長6.6 km、南北寬2.2 km、深度418 m、面積約10.85 km2的亞洲第一大露天礦坑[1].2010年,西露天礦南幫受采礦活動及坡體軟弱夾層的綜合影響開始出現滑坡跡象,形成了東西長約3.1 km、南北寬1.2～1.5 km、高差約102 m的特大型推移式順層巖質滑坡,滑坡區如圖2中綠色虛線區域所示[26-27].礦坑E400剖面工程地質分布如圖3所示.

圖2 研究區Google Earth影像

撫順西露天礦南幫為典型的順層邊坡,自下而上依次分布有片麻巖、玄武巖、凝灰巖等,并夾雜有少量薄煤層,發育有東西向的F2斷層和北西—南東向的F5斷層.其中,F5斷層將南幫一分為二,東側以玄武巖為主,坑底區域有少量凝灰巖覆蓋于玄武巖上方,高陡邊坡處的玄武巖巖體相對穩定;西側為特大型順層滑坡體,內部分布有多條局部斷層,與F5斷層存在著伴生關系[26].

2.2 數據源

本文收集了兩組覆蓋西露天礦的SAR影像,共計109幅,其中包括:2014-03-16～2016-12-18期間的7幅COSMO-SkyMed影像;2017-01-11～2020-05-31期間的102幅Sentinel-1B影像.兩組數據的空間覆蓋范圍見圖4,詳細參數如表1所示.

圖3 撫順西露天礦坑E400剖面典型工程地質分布圖

圖4 SAR影像覆蓋范圍

表1 雷達數據參數

由于滑坡通常位于地形起伏較大的區域,高空間分辨率及高精度的DEM數據對提高基于時序SAR數據的滑坡監測精度至關重要,具體原因包括2點:1)MTI處理過程中,需要使用DEM去除地形相位分量.目前大多數研究成果采用的SRTM DEM或TanDEM-X DEM(只有90 m×90 m)空間分辨率和高程精度都較低,導致地形起伏較大區域的MTI結果中存在明顯的地形殘差；2)在基于時序SAR的滑坡監測中,由于地形條件復雜,高分辨率高精度的DEM有利于確定滑坡體的滑動方向與Az-LOS向之間的空間幾何關系,進而實現高精度坡向變形參數解算.因此,本文使用了空間分辨率為1.3 m、高程精度為cm級的LiDAR DEM,如圖5所示.

圖5 激光雷達掃描生成的高精度DEM

3 實驗結果分析

3.1 2014—2016年快速變形期POT監測結果

2014—2016年間,南幫特大型滑坡處于快速滑移狀態.為了監測該時段的大尺度滑坡變形特征,本文對2014-03-16～2016-12-18期間的7幅高分辨率COSMO-SkyMed數據開展了POT處理.為減小時間失相干的影響,對每間隔半年左右的影像進行分時段跟蹤,分別測量各時段內的Az-LOS變形量,然后根據本文提出的方法轉換到坡向,結果如圖6所示.

圖6顯示,南幫滑坡體的滑移量從空間分布上看,呈現出西半部大、東半部較小的特點,其中最大處位于南幫千臺山及以北的邊坡處.該結果的空間分布特征與文獻[28]中的GPS監測結果基本一致.從時間上來看,2014-03-16～2014-12-29期間的兩個監測時段內,滑坡體大部分區域的累計滑移量都超過12.5 m,最大處超過20 m.2014-12-29～2015-07-09期間,南幫滑移量減小,大部分區域未超過7.5 m;2015-07-09～2015-12-16期間,滑移量再次增大,但仍明顯小于2014年同期,最大值未超過12.5 m;2015-12-16～2016-06-09期間,南幫滑移量再次減小,大部分區域未超過7.5m;2016-06-09～2016-12-18期間,滑移量再次增大,但大部分區域在12.5 m以內.

圖6 2014-03-16～2016-12-18期間的分時段坡向滑移量

為了更直觀展示南幫滑坡體2014—2016年間的時間變化規律,繪制了滑坡軸線附近某剖面的累計滑移量.該剖面經度為123.88°E,位置如圖6f中的黑色虛線所示,剖面滑移歷史如圖7所示,其中黑色箭頭為滑坡體后緣次級地裂縫所在位置,兩側滑移量差異較大.地裂縫以北的紅色虛線框內,滑移量明顯大于其他區域,為該次級地裂縫前緣的拉張斷陷帶.坑底部的藍色虛線框內,滑移量由南到北逐漸減小,為礦坑底部滑坡前緣的隆起帶.從該滑坡剖面的滑動歷史上來看, 2014-03-16～2014-07-22期間累計滑移量最大值為24.588 m,至2014-12-19期間累計滑移量最大值為39.137 m,至2015-07-09期間累計滑移量最大值為45.770 m,至2015-12-16期間累計滑移量最大值為55.879 m,至2016-06-09期間累計滑移量最大值為61.531 m,至2016-12-18期間累計滑移量最大值為73.309 m.總體來說,該滑坡體呈現出2014年快速滑移、2015—2016年速度稍微減慢、秋冬季減速、雨季加速的特征.

3.2 2017—2020年緩慢變形期MTI監測結果

2017年以來西露天礦將E1000以西區域作為東露天礦的排土場,可對南幫的大型滑坡體實現內排壓腳治理[29].經過一段時間的排放后,截止2020年7月已回填至-200 m標高左右.同時,西露天礦對南幫的特大型滑坡體實施了注漿、抗滑樁、防治水等綜合治理方案,最終使該特大型滑坡體基本回歸穩定[29].本文使用MTI方法處理了2017-01-11～2020-05-31期間的Sentinel-1B影像,然后將LOS向結果轉換到坡向,所得的坡向累計滑移量如圖8所示.

圖7 南幫滑坡剖面(123.88°E)2014-03-16～2016-12-18期間累計滑移量

圖8 2017-01-11～2020-05-31期間南幫累計 滑移量

圖8表明,南幫滑坡體2017年年初至2020-05-31期間大部分區域的累計滑移量在300 mm以內,滑移量最大的區域位于E1300附近(原劉山舊河道位置).

圖8中所標4個點的累計坡向滑移曲線如圖9所示.其中P1～P3點在2014—2016年間布設了GPS監測設備,P4點位于滑移量最大的區域.從圖9可看出,P1點與P3點在該時段內累計滑移量均未超過200 mm,而在2014-03-16～2016-12-18期間POT測得的滑移量分別為33.099 m和60.810 m.P2點在該時段的滑移量約為300 mm,而2014-03-16～2016-12-18期間POT測得的滑移量為50.732 m.這三個點在兩個時段的滑移量差異巨大,表明南幫滑坡體已基本穩定.P4點由于位于劉山舊河道附近,該時段內的累計滑移量接近600 mm,且與月降水量存在明顯的相關性.尤其是在2019年8月的強降雨后,P4點的滑移速率明顯增大,但2019年年底速率已經減小.這是由于舊河道附近基巖下切,存在河卵石層,滲透性良好,使得邊坡裂縫容易滲入雨水及地表水,增加了巖土層的流動性,降低了抗拉強度.

圖9 2017-01-11～2020-05-31期間P1～P4點的 滑移曲線

4 監測精度評價

為了驗證本文Az-LOS向至坡向變形轉換方法的精度,將2014-03-16～2016-06-09期間P1～P3點的GPS監測數據與本文監測結果進行分時段的對比分析,結果如表2所示.其中,P1點的POT監測結果在各時段內相對于GPS監測結果的誤差均小于1 m,誤差均值為0.130 m,RMSE為0.489 m.POT測得的累計滑移量為33.099 m,GPS測得的累計滑移量為32.447 m,誤差為0.652 m.P2點POT監測結果在2014-07-22～2014-12-29期間與GPS結果相比誤差較大,為2.78 m,其余時段誤差均在1 m左右或少于1 m.P2點的POT監測結果與GPS結果相比誤差均值為-0.662 m,RMSE為1.349 m;累計偏移量分別為50.732和54.041 m,誤差為-3.309 m,誤差主要來源于2014-07-22～2014-12-29期間的配準誤差.P3點在各時段的POT監測結果與GPS結果相比誤差均為dm級,誤差均值為-0.067 m,RMSE為0.459 m.P3點POT方法和GPS測得的累計滑移量分別為60.810和61.144 m,誤差為-0.334 m.總體來說,3個點的POT監測結果與GPS監測結果的誤差均值均在dm級,表明本文提出的轉換方法在基于POT的滑坡快速變形監測應用中精度較高.

表2 南幫3個點累計滑移量精度評價

5 結 論

1) 撫順西露天礦南幫滑坡是一個空間分布范圍廣、發育時間長的特大型滑坡,其變形具有階段性特征,不同時段內的滑移速率差異巨大,通過協同使用POT和MTI方法,獲得了滑坡體長時序、高精度的時空發育特征.

2) 滑坡易發區地形及地質條件復雜,開展基于時序SAR數據的滑坡監測時,應充分考慮觀測值方向與滑坡變形的空間一致性問題,實現Az-LOS向到坡向變形的轉換,提取沿滑坡坡面的變形參數,提高監測精度.

3) 西露天礦南幫特大型滑坡體2014年處于快速滑移階段;2015—2016年間滑移速度減慢,呈現出秋冬季減速、雨季加速的特點;2017年初以后,南幫滑坡體回歸穩定,僅在E1300附近存在小范圍蠕滑.將監測結果與現場GPS監測數據及月降水量數據進行對比分析,精度較高,且與文獻中的結論一致.

4) 通過對西露天礦南幫滑坡成因進行分析,發現斷層、軟弱地層、地下水賦存等復雜地質條件及降水量是影響變形的主要因素,在開展露天礦滑坡等地質災害監測工作時不能忽視其影響.