基于SBAS-InSAR的西寧市滑坡識(shí)別和形變監(jiān)測(cè)分析

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

(1.青海省青藏高原北部地質(zhì)過(guò)程與礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室,西寧 810000; 2.成都理工大學(xué) a.地理與規(guī)劃學(xué)院,b.成都理工大學(xué) 地球科學(xué)學(xué)院,成都 610059)

0 引言

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

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

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

1 研究區(qū)概況與數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

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

圖1 研究區(qū)位置

1.2 數(shù)據(jù)來(lái)源與數(shù)據(jù)處理

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

表1 Sentinel-1A衛(wèi)星基礎(chǔ)數(shù)據(jù)

圖2 技術(shù)流程圖

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

2 實(shí)驗(yàn)方法

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

1)對(duì)原始的Sentinel-1A數(shù)據(jù)進(jìn)行軌道信息校正以及裁剪到覆蓋整個(gè)研究區(qū)范圍。

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

圖3 SBAS-InSAR時(shí)空基線圖

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

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

3.1 SBAS-InSAR識(shí)別活動(dòng)滑坡

通過(guò)SBAS-InSAR技術(shù)獲取西寧市Sentinel-1A升軌和降軌數(shù)據(jù)雷達(dá)視線(Line of Sight,LOS)方向上的形變速率(圖4、圖5)。圖4、圖5中紅色的點(diǎn)(負(fù)值)表示目標(biāo)沿著LOS方向遠(yuǎn)離衛(wèi)星移動(dòng),綠色的點(diǎn)表示在監(jiān)測(cè)時(shí)間段內(nèi)目標(biāo)相對(duì)穩(wěn)定,藍(lán)色點(diǎn)(正值)表示目標(biāo)沿著LOS方向靠近衛(wèi)星移動(dòng)。研究區(qū)升軌形變速率位于-98 mm/年～44 mm/年之間(圖4(a)),降軌形變速率位于-95 mm/年～25 mm/年之間(圖5(a))。由圖4(b)、圖5(b)可知,升軌和降軌形變速率均主要集中在-10 mm/年～10 mm/年的范圍內(nèi),平均形變速率分別為-2.09 mm/年、-0.64 mm/年,且分布特征與正態(tài)分布相似,這表明研究區(qū)整體上處于穩(wěn)定狀態(tài)。基于形變速率信息,筆者將形變速率超過(guò)-15 mm/年區(qū)域確定為顯著形變區(qū)。研究區(qū)共探測(cè)到74處顯著形變區(qū)(升軌數(shù)據(jù)53處,降軌數(shù)據(jù)67處,升軌和降軌數(shù)據(jù)聯(lián)合46處),總面積約9 km2。

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

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

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

圖6 西寧市不同時(shí)期遙感影像圖

圖7 識(shí)別區(qū)結(jié)果分布及滑坡主要發(fā)育區(qū)局部放大圖

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

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

表2 重點(diǎn)區(qū)InSAR識(shí)別滑坡的基本特征

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

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

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

圖9 研究區(qū)其他形變區(qū)展示圖

4 討論

4.1 滑坡時(shí)間序列形變與降雨關(guān)系

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

4.1.1 張家灣滑坡形變與降雨關(guān)系

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

圖10 張家灣滑坡InSAR識(shí)別結(jié)果

4.1.2 付家寨滑坡形變與降雨關(guān)系

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

圖11 付家寨滑坡InSAR探測(cè)結(jié)果

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

圖11(f)為P1特征點(diǎn)、P2特征點(diǎn)累計(jì)形變量與日平均降雨量關(guān)系圖。可知累計(jì)形變曲線可分為兩個(gè)階段,2019年1月-2020年10月時(shí)間段形變量快速增加,P1點(diǎn)、P2點(diǎn)累計(jì)形變量分別增加81 mm、89 mm。2020年10月之后形變趨勢(shì)變緩慢,整體趨于穩(wěn)定。形變曲線在2019年5月和2020年5月雨季之后有明顯的加速變形趨勢(shì)(圖中淺藍(lán)色區(qū)域),說(shuō)明該滑坡變形與日降雨量密切相關(guān)。

5 結(jié)論

研究區(qū)位于黃土地區(qū),地表黃土物質(zhì)廣泛覆蓋,植被稀少,區(qū)域內(nèi)相干性良好,采用升降軌Sentinel-1A數(shù)據(jù),基于時(shí)序InSAR技術(shù)對(duì)研究區(qū)活動(dòng)滑坡進(jìn)行識(shí)別和形變監(jiān)測(cè)。結(jié)論如下:

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

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