基于SBAS-InSAR的西寧市滑坡識別和形變監測分析

何 麗,吳先譚,陳 欣,羅 芳,何政偉,薛東劍,白文倩,康桂川,張雨祥

(1.青海省青藏高原北部地質過程與礦產資源重點實驗室,西寧 810000; 2.成都理工大學 a.地理與規劃學院,b.成都理工大學 地球科學學院,成都 610059)

0 引言

中國黃土高原的面積約6.4×105km2,黃土厚度約為100 m～300 m[1]。黃土具有孔隙大、質地疏松、垂直節理發育等獨特的構造特征,在降雨的條件下極易發生滑坡、崩塌、泥石流等地質災害[2]。黃土滑坡是黃土高原分布最廣的一類地質災害[3]。近年來,由于黃土地區日益脆弱的地質環境,在強降雨條件下黃土斜坡極易發生傾倒、滑移,對人類生命財產威脅極大[4]。2010年榆林市發生大規模黃土崩塌,造成了嚴重的人員傷亡和財產損失[5](如2013年天水市受強降雨影響,誘發大規模黃土滑坡[6])。因此,大范圍高效、精準地識別與監測具有潛在危險性的蠕變型黃土滑坡,對于滑坡防治、災害預警以及保護人民生命財產安全具有重要意義。

隨著遙感技術的發展,光學遙感[7]、無人機(UAV)攝影測量[8]、機載激光雷達(LiDAR)[9]和合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar,InSAR[10])等技術逐漸成為滑坡識別和形變監測的主要手段。其中InSAR技術因具備高精度、高分辨率、覆蓋范圍大、綜合成本低、全天時和連續跟蹤微小形變等技術優勢,使其在大范圍滑坡識別和形變監測應用中得到青睞[11-13]。1996年Fruneau[14]首次證明了合成孔徑雷達差分干涉測量技術(D-InSAR)在滑坡形變監測的有效性,同時也發現時空失相干和大氣延遲對形變監測結果影響較大。隨后,國際學者在D-InSAR技術的基礎上提出了基于多景SAR影像數據的時間序列InSAR技術(TS-InSAR),包括永久散射體干涉測量(Persistent O Scatter InSAR,PS-InSAR)[15]、小基線集干涉測量(Small Baseline InSAR,SBAS-InSAR)[16],分布式散射體干涉測量SqueeSAR[17]等技術,這些技術能夠獲取精度更高、可靠性更強的滑坡形變結果,在大面積滑坡識別和形變監測中已得到廣泛應用。目前,以SBAS為代表的TS-InSAR已被用于黃土滑坡的識別和形變監測[18-19]。如劉陳偉等[20]利用SBAS技術獲取黑方臺地區地面形變速率,并結合黃土滑坡的地形地貌特征,識別出黑方臺地區6塊存在較大形變的區域;Meng等[21]利用無人機(UAV)攝影測量與SBAS-InSAR技術表征甘肅省紅壑峴黃土滑坡形變過程,對比分析發現SBAS-InSAR技術在蠕變變形監測方面更具有優勢。

筆者以升降軌Sentinel-1A雷達影像作為數據源,在對時序InSAR技術分析的基礎上,采用SBAS-InSAR技術,提取了西寧市的分布式散射體(Distributed Scatterers,DS),反演了研究區的形變速率和時間序列形變量,并結合Google Earth遙感影像圈定了活動滑坡邊界,分析了西寧市活動滑坡的分布規律,以及滑坡累計形變與降雨的關系。本研究成果為同類型區域黃土滑坡識別和形變監測提供了借鑒。

1 研究區概況與數據來源

1.1 研究區概況

西寧市位于青藏高原東北部,隸屬于湟水中游河谷盆地,地理位置:101°33′E～101°56'E;36°25'N～36°47N,區域面積為490 km2。研究區地勢西北高東南低,地形切割強烈,溝壑縱橫,區內發育一條近東西向的湟水河(圖1)。研究區地層由老到新為長城系(Pta)、白堊系(K)、古近系(E)、新近系(N)、中更新統(Q2)、上更新統(Q3)、全新統(Q4)。巖土類型主要有黃土、泥巖砂巖和石膏巖互層、松散碎石土,均為易失穩變形地層。西寧市屬高原半干旱大陸性氣候,西寧地區全年降雨分配不均勻,年降雨量477.5 mm,主要集中在6月-9月,占全年總降雨量的67.4%,年降雨周期性變化明顯。研究區地形主要以侵蝕構造低山丘陵與侵蝕堆積河谷平原為主,在低山丘陵前緣多發育高陡斜坡,加之區域降水集中,極易誘發滑坡、崩塌和泥石流等地質災害。

圖1 研究區位置

1.2 數據來源與數據處理

Sentinel-1A是由歐洲航天局(ESA)發射的C波段合成孔徑雷達對地觀測衛星,該衛星具有全天時、多種極化方式、時空分辨率高等特征。本研究采用Sentinel-1A衛星以干涉寬幅(IW)模式下VV極化方式獲取的單視復數產品(SLC)進行時間序列形變反演,影像采集了2019年1月至2021年11月83幅升軌影像和87幅降軌影像。升軌影像與降軌影像在不同的入射角下,存在不同程度的幾何畸變,形成不同的觀測效果。結合升軌與降軌數據進行滑坡識別,能更有效地擴大監測范圍。SAR數據的基本參數如表1所示,具體的技術流程如圖2所示。

表1 Sentinel-1A衛星基礎數據

圖2 技術流程圖

此外,筆者采用30 m分辨率的SRTM(shuttle radar topography mission) DEM(digital elevation model)去除干涉過程中的地形相位。選用POD Precise Orbit Ephemerides(POD精密定軌星歷數據)(https://qc.sentinel1.eo.esa.int/)校正軌道信息,選用Generic Atmospheric Correction Online Service for InSAR(GACOS)(http://www.gacos.net/)大氣延遲產品數據,該數據時間和空間分辨率較高,全球近實時,可有效去除干涉圖中的大氣延遲相位收集了西寧市氣象站2019年-2021年日平均降水量數據,用于分析滑坡形變與降雨量之間的關系。

2 實驗方法

Berardino等[16]提出的SBAS-InSAR技術是對分布式散射體(Distributed scatterer DS)進行相位分析獲得時序形變。該技術對基于低分辨率、大尺度、時序形變監測能夠取得很好的效果[22]。與傳統的D-InSAR技術相比較而言,SBAS-InSAR技術能夠較好地克服時空失相關和大氣延遲效應等缺陷,實現更高精度的地表形變監測[23]。SBAS-InSAR技術實現過程為:

1)對原始的Sentinel-1A數據進行軌道信息校正以及裁剪到覆蓋整個研究區范圍。

2)綜合考慮時間-空間基線,分別選取2019年3月10日和2019年11月10日的影像作為升軌和降軌數據的配準主影像。設置60 d時間基線,150 m空間基線,分別篩選出320個和304個干涉對(圖3)。對干涉對進行多視、干涉圖生成、去平地相位、Goldstein濾波、最小費用流(MCF)相位解纏等處理獲得每個干涉圖中任意像元(x,y)的差分干涉解纏相位。

圖3 SBAS-InSAR時空基線圖

3)形變速率和累計形變估算,形變分量可以進一步分解為線性形變分量和非線性形變分量。通過相位與時間基線間的最小二乘擬合,可以估算出線性位移。使用奇異值分解法(SVD)求出非線性形變最小范數意義上的解。

3 實驗結果

3.1 SBAS-InSAR識別活動滑坡

通過SBAS-InSAR技術獲取西寧市Sentinel-1A升軌和降軌數據雷達視線(Line of Sight,LOS)方向上的形變速率(圖4、圖5)。圖4、圖5中紅色的點(負值)表示目標沿著LOS方向遠離衛星移動,綠色的點表示在監測時間段內目標相對穩定,藍色點(正值)表示目標沿著LOS方向靠近衛星移動。研究區升軌形變速率位于-98 mm/年～44 mm/年之間(圖4(a)),降軌形變速率位于-95 mm/年～25 mm/年之間(圖5(a))。由圖4(b)、圖5(b)可知,升軌和降軌形變速率均主要集中在-10 mm/年～10 mm/年的范圍內,平均形變速率分別為-2.09 mm/年、-0.64 mm/年,且分布特征與正態分布相似,這表明研究區整體上處于穩定狀態?；谛巫兯俾市畔?筆者將形變速率超過-15 mm/年區域確定為顯著形變區。研究區共探測到74處顯著形變區(升軌數據53處,降軌數據67處,升軌和降軌數據聯合46處),總面積約9 km2。

圖4 Sentinel-1A升軌形變速率圖

圖5 Sentinel-1A降軌形變速率圖

為了進一步判斷形變區是否為滑坡及圈定滑坡邊界范圍,結合顯著形變區的分布情況,利用多期次Google Earth光學遙感影像和高分2號影像數據(圖6),可以清晰地觀察到地表物體形態。結合滑坡在光學遙感色調、紋理、平面形態、微地貌形態以及變形特征等解譯標志,特別是在高分辨率影像中表現出明顯的裂縫及滑塌現象,對SBAS-InSAR技術探測到的顯著形變區進一步劃分,最終劃定了升軌滑坡、降軌滑坡以及其他形變區的分布范圍(圖7)。筆者將35處具有明顯滑坡光學影像解譯特征的顯著形變區識別為滑坡,并在Google Earth影像中勾畫出滑坡邊界,39處具有明顯人類活動特征的顯著形變區識別為其他形變區(圖7(a))。

圖6 西寧市不同時期遙感影像圖

圖7 識別區結果分布及滑坡主要發育區局部放大圖

經過分析,InSAR技術識別的35處滑坡中,升軌數據識別出14處、降軌數據識別出28處、升軌和降軌數據聯合識別出7處(圖7(a))。識別出的滑坡總面積約2.4 km2,其中面積最小為5 945.64 m2,最大為0.56 km2,規模以中大型為主?；鲁蕳l帶狀集中分布于湟水河及其一級支流河谷平原區向低山丘陵區過渡的斜坡地帶,且湟水河北岸滑坡發育密度大于南岸。圖7(b)、圖7(c)為研究區滑坡典型發育區,分別記為A區和B區。

圖8(a)、圖8(b)分別表示A和B典型區InSAR降軌數據形變速率空間三維分布,分析可知,A區和B區域分別發育4處滑坡(H1、H2、H3、H4)和6處滑坡(H5、H6、H7、H8、H9、H10),其基本特征如表2所示。

表2 重點區InSAR識別滑坡的基本特征

圖8 研究區典型滑坡形變速率及空間分布

由表2及圖8(c)、圖8(d)可知,A區域中H2滑坡最大形變速率達48.87 mm/年,坡體形變集中在中部和前緣,在滑坡體前緣和中部發育3處次級滑坡。后緣下錯約10 m并且發育多條走向近南北裂縫,長約180 m。中部有一系列橫向延伸長度70 m～90 m的不連續圓弧形拉裂縫,且呈逐步貫通趨勢,在自重作用下,上部巖土不斷向下擠壓。由表2及圖8(e)、圖8(f)可知,B區域中H8滑坡最大形變速率達81.26 mm/年,坡體形變集中在中部區域,滑體后部發育有三條裂縫,均呈東西向展布,平面形狀呈弧形,最長一條達到170 m,寬0.2 m～0.6 m,滑坡體北側分布呈帶狀陡崖,坡度近90°?，F陡崖上發育多處危巖,危巖與母巖間裂隙寬10 cm～30 cm,上下貫通,呈孤立狀(圖8(d))。

識別到的39處其他形變區在影像可清晰識別為挖填區(治溝造地、土地整平、垃圾填埋、流域治理)、建筑施工區、道路沉降區、礦區和工業區等。經過現場調查核實,最終確定13處為挖填區、16處為建筑施工區、2處為道路沉降區、4處為礦區和4處為工業區。圖9(a)～圖9(c)分別顯示了建筑施工、工業區和治溝造地引起的地表形變。

圖9 研究區其他形變區展示圖

4 討論

4.1 滑坡時間序列形變與降雨關系

降雨是誘發黃土滑坡的重要因素之一,這與黃土孔隙比、濕陷性等獨特性質密切相關。本研究以張家灣滑坡和付家寨滑坡為例,分析降雨與滑坡時間序列累計形變(SBAS-InSAR)的演化關系。

4.1.1 張家灣滑坡形變與降雨關系

該滑坡地理坐標為101°39′04″E,36°38′25″N?；赂卟罴s325 m,坡度30°,面積5.6×106m2,總體積2.2×106m3,表層披覆黃土,屬平緩層狀巖質斜坡,為特大型推移式滑坡?；鲁什灰巹t狀,坡體有多處淺表層滑塌跡象,后緣裂縫明顯(圖10(a))。圖10(b)為該滑坡Sentinel-1A升軌影像LOS方向形變速率圖,形變速率范圍為-25 mm/年～10 mm/年,形變部位集中在滑坡后緣及中部具有明顯拉裂縫的區域。根據形變點分布規律,可將滑坡分為四個顯著性形變區(Ⅰ、Ⅱ、Ⅲ、Ⅳ)。經野外調查可得,滑坡Ⅰ區后緣貫通,有土質垮塌的現象,前緣的擋墻有開裂(圖10(c))。Ⅱ區因人類工程活動導致部分地帶有失穩的可能。Ⅲ區發育走向130°的張性拉裂縫,下錯明顯約80 cm,形成新的不穩定體,目前該區已進行了滑坡治理。Ⅳ區后緣可見局部的垮塌,且發育多條裂縫(圖10(d))?；麦w上季節性降雨導致裂縫、落水洞、座落等現象時有發生。分別在四個區域取一點(P1、P2、P3、P4)進行時間序列累計形變量與日平均降雨量關系分析(圖10(e)),四個區域的累計形變量均呈現波動增加,其中P1點累計形變量最大為-78.6 mm,形變曲線反映該區域處于快速形變階段?；略?019年6月-2020年2月和2020年8月-2021年2月(淺藍色區域),由于雨季和凍融沉降的到來,形變量呈快速增加趨勢,可見滑坡形變與日平均降雨量密切相關,具有明顯的季節性特征,尤其是在雨季變形明顯加速,并且變形加速時間與雨季開始時間存在明顯的滯后效應。

圖10 張家灣滑坡InSAR識別結果

4.1.2 付家寨滑坡形變與降雨關系

該滑坡于2019年11月18日發生在城東區付家寨北部,地理坐標:E:101°54′10.54″、N:36°34′52.03″?；麻L270 m,橫310 m,滑向295°,平均坡度25°,平面形狀近似正方形,滑坡體較為破碎,發育多個次級滑坡,多條裂縫貫穿,后緣下錯明顯(圖11(a)～圖11(b))。

圖11 付家寨滑坡InSAR探測結果

圖11(c)為該滑坡Sentinel-1A降軌影像LOS方向形變速率圖,坡體形變集中在中前部區域,最大形變速率達40 mm/年,結合地形條件和形變空間分布特征可將滑坡分為兩處顯著形變區(Ⅰ、Ⅱ)。在兩處形變顯著區域各選1個特征點P1、P2,其中P1平均形變速率為31 mm/年,累計形變量為92.7 mm。該形變區前緣堆積體堆覆溝道內,擠壓溝道,堆積高度達3 m,災害造成4間房屋損毀(圖11(d))。P2平均形變速率為38 mm/年,累計形變量為105.6 mm。該形變區為一中型牽引式次級滑坡,滑體中部樹木歪斜,綠化管網斷裂,道路垮塌,后緣下錯約5 m,前緣堆積體堆覆溝道內,擠壓溝道,堆積高度達3 m(圖11(e))。

圖11(f)為P1特征點、P2特征點累計形變量與日平均降雨量關系圖?？芍塾嬓巫兦€可分為兩個階段,2019年1月-2020年10月時間段形變量快速增加,P1點、P2點累計形變量分別增加81 mm、89 mm。2020年10月之后形變趨勢變緩慢,整體趨于穩定。形變曲線在2019年5月和2020年5月雨季之后有明顯的加速變形趨勢(圖中淺藍色區域),說明該滑坡變形與日降雨量密切相關。

5 結論

研究區位于黃土地區,地表黃土物質廣泛覆蓋,植被稀少,區域內相干性良好,采用升降軌Sentinel-1A數據,基于時序InSAR技術對研究區活動滑坡進行識別和形變監測。結論如下:

1)通過對原始Sentinel-1A數據進行裁剪,綜合考慮時間和空間基線,篩選出若干組干涉對,隨后進行干涉處理(影像配準、生成干涉圖、去地形相位、相位解纏等)以及軌道精煉和重去平,最后對形變速率進行估算,識別出區域的顯著形變區域。同時結合高分辨率衛星影像,進一步判定及劃定滑坡邊界范圍。共探測到74處顯著形變區,經過分析,35處為活動滑坡,其中升軌識別出14處、降軌識別出28處、升軌和降軌聯合識別出7處。升軌和降軌數據相結合,能夠有效地減少因幾何畸變導致的監測“盲區”,從而提高衛星數據的有效探測率。

2)以張家灣滑坡和付家寨滑坡兩個典型滑坡為例,分析了降雨與滑坡時間序列累計形變與降水的關系。分析發現滑坡受降雨入滲的影響,造成黃土、部分泥巖、石膏塊軟化與崩解,從而使土體物理力學性質發生改變,相對集中且強烈的降雨最終導致滑坡穩定性降低,內在(工程巖組、斷層)和外在(降雨)條件的共同作用導致已識別的滑坡在研究時段處于活動狀態。