聯(lián)合LiDAR DEM與時(shí)序SAR數(shù)據(jù)的露天礦特大型滑坡監(jiān)測(cè)

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

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

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

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

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

1 方 法

1.1 POT方法及Az-LOS向到坡向變形轉(zhuǎn)換

傳統(tǒng)的基于光學(xué)影像的POT方法受到圖像噪聲、傳感器姿態(tài)等因素的限制,精度一般不高,且存在無法估算垂直位移、只能估算水平位移的問題[10-11].基于SAR影像的POT方法不受相位解纏及去相干問題的影響,通過使用歸一化互相關(guān)優(yōu)化窗口,跟蹤不同時(shí)相SAR圖像之間后向散射特征點(diǎn)位置的變化(即偏移量),可以同時(shí)測(cè)量Az-LOS向的二維變形[12-15].由于SAR圖像之間的配準(zhǔn)誤差通常為像素分辨率的1/10～1/20,當(dāng)?shù)乇砦灰屏啃∮谂錅?zhǔn)誤差時(shí),POT將難以提取精確的偏移量,因此POT適合監(jiān)測(cè)快速變形[6].一般來說在其他條件相同時(shí),SAR影像分辨率越高,POT對(duì)變形的敏感度及監(jiān)測(cè)精度也相應(yīng)越高.

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

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

(1)

(2)

式中:dAz為方位向變形量;dLOS為L(zhǎng)OS向變形量;β為坡度;α0為航向角(順時(shí)針為正)；αA為坡向角;θ為入射角.將式(1),式(2)聯(lián)立消去β可得坡向變形量dA的計(jì)算公式:

(3)

1.2 MTI方法及LOS到坡向變形轉(zhuǎn)換

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

針對(duì)基于MTI的滑坡監(jiān)測(cè)中面臨的LOS向變形到坡度方向的轉(zhuǎn)換問題,作者曾提出一種基于LiDAR DEM輔助的LOS向到坡向變形的轉(zhuǎn)換方法,已在鞍鋼集團(tuán)各大鐵礦高陡邊坡監(jiān)測(cè)中取得了成功應(yīng)用[24-25].本文繼續(xù)采用這一方法實(shí)現(xiàn)LOS向到坡向的轉(zhuǎn)換.轉(zhuǎn)換過程中,需要根據(jù)LiDAR DEM計(jì)算礦坑內(nèi)精確的坡度、坡向參數(shù),結(jié)合雷達(dá)衛(wèi)星數(shù)據(jù)的成像幾何求解坡向變形.

2 實(shí)驗(yàn)區(qū)及數(shù)據(jù)

2.1 實(shí)驗(yàn)區(qū)工程地質(zhì)條件

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

圖2 研究區(qū)Google Earth影像

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

2.2 數(shù)據(jù)源

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

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

圖4 SAR影像覆蓋范圍

表1 雷達(dá)數(shù)據(jù)參數(shù)

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

圖5 激光雷達(dá)掃描生成的高精度DEM

3 實(shí)驗(yàn)結(jié)果分析

3.1 2014—2016年快速變形期POT監(jiān)測(cè)結(jié)果

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

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

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

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

3.2 2017—2020年緩慢變形期MTI監(jiān)測(cè)結(jié)果

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

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

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

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

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

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

4 監(jiān)測(cè)精度評(píng)價(jià)

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

表2 南幫3個(gè)點(diǎn)累計(jì)滑移量精度評(píng)價(jià)

5 結(jié) 論

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

2) 滑坡易發(fā)區(qū)地形及地質(zhì)條件復(fù)雜,開展基于時(shí)序SAR數(shù)據(jù)的滑坡監(jiān)測(cè)時(shí),應(yīng)充分考慮觀測(cè)值方向與滑坡變形的空間一致性問題,實(shí)現(xiàn)Az-LOS向到坡向變形的轉(zhuǎn)換,提取沿滑坡坡面的變形參數(shù),提高監(jiān)測(cè)精度.

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

4) 通過對(duì)西露天礦南幫滑坡成因進(jìn)行分析,發(fā)現(xiàn)斷層、軟弱地層、地下水賦存等復(fù)雜地質(zhì)條件及降水量是影響變形的主要因素,在開展露天礦滑坡等地質(zhì)災(zāi)害監(jiān)測(cè)工作時(shí)不能忽視其影響.