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

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

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

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

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

1 方 法

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

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

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

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

(1)

(2)

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

(3)

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

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

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

2 實驗區(qū)及數(shù)據(jù)

2.1 實驗區(qū)工程地質(zhì)條件

撫順西露天煤礦位于遼寧省撫順市西南部,最早開采于1901年,1914年開始露天開采,至今形成了東西長6.6 km、南北寬2.2 km、深度418 m、面積約10.85 km2的亞洲第一大露天礦坑[1].2010年,西露天礦南幫受采礦活動及坡體軟弱夾層的綜合影響開始出現(xiàn)滑坡跡象,形成了東西長約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覆蓋于玄武巖上方,高陡邊坡處的玄武巖巖體相對穩(wěn)定;西側(cè)為特大型順層滑坡體,內(nèi)部分布有多條局部斷層,與F5斷層存在著伴生關(guān)系[26].

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

本文收集了兩組覆蓋西露天礦的SAR影像,共計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 雷達數(shù)據(jù)參數(shù)

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

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

3 實驗結(jié)果分析

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

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

圖6顯示,南幫滑坡體的滑移量從空間分布上看,呈現(xiàn)出西半部大、東半部較小的特點,其中最大處位于南幫千臺山及以北的邊坡處.該結(jié)果的空間分布特征與文獻[28]中的GPS監(jiān)測結(jié)果基本一致.從時間上來看,2014-03-16～2014-12-29期間的兩個監(jiān)測時段內(nèi),滑坡體大部分區(qū)域的累計滑移量都超過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期間的分時段坡向滑移量

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

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

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

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

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

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

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

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

4 監(jiān)測精度評價

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

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

5 結(jié) 論

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

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

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

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