許東麗， 王 濤， 達(dá)瓦桑杰， 楊愛平， 李金鎖， 朱建東

(1.四川地質(zhì)礦產(chǎn)勘查開發(fā)局207地質(zhì)隊(duì)， 樂山 614000； 2.四川嘉源蓉創(chuàng)地質(zhì)科技有限公司， 成都 610000；3.西藏自治區(qū)地質(zhì)礦產(chǎn)勘查開發(fā)局第二地質(zhì)大隊(duì), 拉薩 850003)

中國是滑坡災(zāi)害頻發(fā)國家，尤其是中國西部地區(qū)的滑坡，往往規(guī)模大、機(jī)制復(fù)雜，危害大[1]?；略谛纬砂l(fā)展過程中，往往伴隨著地表形變，通過對這些形變的監(jiān)測是識別及預(yù)測滑坡的重要手段。目前，滑坡監(jiān)測的主要方法有精密水準(zhǔn)測量、全球定位系統(tǒng)(global positioning system，GPS)、位移監(jiān)測等，但這些方法監(jiān)測覆蓋范圍小、工作量大，不適用于大型滑坡監(jiān)測。近年來，隨著合成孔徑雷達(dá)(interferometric synthetic aperture radar,InSAR)技術(shù)的發(fā)展，研究者們拓展了其應(yīng)用范圍，將原本主要用于地面沉降監(jiān)測的技術(shù)——星載InSAR用于滑坡監(jiān)測，取得了一系列良好的結(jié)果。

在InSAR變形監(jiān)測中，早期應(yīng)用的是D-InSAR[2-5](differential InSAR)技術(shù)，在D-InSAR基礎(chǔ)上發(fā)展了長時(shí)間序列InSAR技術(shù)并運(yùn)用在滑坡監(jiān)測中，主要包括永久散射技術(shù)PS-InSAR(permanent scatterer InSAR)和小基線集技術(shù)(SBAS-InSAR)。其中，PS-InSAR適用于城市等干涉條件和輻射比較穩(wěn)定的區(qū)域[6]，牛全福等[7]使用PS-InSAR技術(shù)對蘭州市地面形變進(jìn)行監(jiān)測；而SBAS-InSAR是一種基于多主影像的InSAR時(shí)序方法，只利用時(shí)空基線較短的干涉對來提取地表形變信息，對數(shù)據(jù)的數(shù)量要求較低，適用于山區(qū)等干涉條件不穩(wěn)定的地區(qū)。因此，目前主要將其運(yùn)用于非城區(qū)的線狀工程以及地質(zhì)災(zāi)害(滑坡為主)的形變監(jiān)測。文獻(xiàn)[8-11]使用SBAS-InSAR技術(shù)對川東地區(qū)、金沙江流域滑坡和鮮水河地質(zhì)災(zāi)害進(jìn)行識別和監(jiān)測;文獻(xiàn)[12-13]分別使用SBAS-InSAR技術(shù)對地鐵及高速公路進(jìn)行沉降監(jiān)測，以上研究均證明了SBAS-InSAR在低相干區(qū)也具有準(zhǔn)確估計(jì)形變的能力。

綜上所述，中外眾多學(xué)者采用InSAR技術(shù)對不同類型工程和災(zāi)害的形變特征進(jìn)行研究，取得了較滿意的研究成果。但將InSAR技術(shù)應(yīng)用于青藏高原特大型滑坡形變特征的研究相對較少。因此以規(guī)模大、植被覆蓋率低、處于變形加速階段的西藏自治區(qū)山南市龐村特大型滑坡為研究對象，采用SBAS-InSAR技術(shù)，選取2018年1月—2020年7月共75景Sentinel-1A數(shù)據(jù)，對滑坡區(qū)域的形變速率、累計(jì)形變量、時(shí)序形變特征及空間形變特征進(jìn)行分析，并與地面監(jiān)測站實(shí)測的位移數(shù)據(jù)進(jìn)行精度對比，驗(yàn)證了SBAS-InSAR技術(shù)應(yīng)用于青藏高原特大型滑坡形變監(jiān)測的可靠性，以期為青藏高原特大型滑坡的監(jiān)測工作提供有益參考。

1 滑坡特征與數(shù)據(jù)源

1.1 滑坡基本特征

龐村特大型滑坡位于西藏自治區(qū)山南市隆子縣的加玉鄉(xiāng)東側(cè)，加玉河右岸(圖1)。該滑坡為一古滑坡，滑坡整體地形呈南高北低，上陡下緩，坡度20°～40°，局部地段可達(dá)50°以上，平均坡度約23°，整體面積約0.3 km2，整體體積約1.1×107m2，屬于特大型滑坡。滑坡平面形態(tài)呈圈椅狀，剖面形態(tài)呈“階梯”形，平面形態(tài)呈舌形，主滑方向約14°?；麦w后緣拉張裂縫已基本貫通并形成下錯陡坎，基巖出露，最大下錯達(dá)5.2 m，滑坡后壁可見明顯擦痕；滑坡前緣剪出口因G219改道開挖，出現(xiàn)不同程度垮塌(圖2)。當(dāng)?shù)卮迕裼?019年4月發(fā)現(xiàn)滑坡裂縫，隨著雨季和地震的影響，裂縫持續(xù)增大，滑坡整體處于不穩(wěn)定狀態(tài)。

圖1 研究區(qū)域衛(wèi)星影像圖Fig.1 Satellite imagery of the study area

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

1.2.1 雷達(dá)數(shù)據(jù)

本文中數(shù)據(jù)源選用Sentinel-1A衛(wèi)星數(shù)據(jù)，Sentinel-1A衛(wèi)星是由歐空局發(fā)射的免費(fèi)對地觀測雷達(dá)衛(wèi)星，軌道高度693 km，多極化，影像采用TOPS(terrain observation with progressive scans in azimuth)模式成像，重訪周期為12 d，雷達(dá)波段為C波段，波長5.6 cm。數(shù)據(jù)空間分辨率30 m，本文中所采用的數(shù)據(jù)時(shí)間跨度為2018年1月3日—2020年7月3日，總計(jì)75景影像，數(shù)據(jù)來源(https://qc.sentinel1.eo.esa.int/)。數(shù)據(jù)的具體時(shí)序信息如表1所示，覆蓋范圍如圖1所示。

表1 數(shù)據(jù)時(shí)間序列列表Table 1 Time series list of data

同時(shí)，由于研究區(qū)域所處區(qū)域地形起伏較大，為消除地形相位對形變相位的影響，本文中采用30 m分辨率的SRTM(shuttle radar topography mission)數(shù)據(jù)作為輔助數(shù)據(jù)，計(jì)算地形相位分量。

1.2.2 地面監(jiān)測數(shù)據(jù)

相關(guān)管理部門于2019年8月初在滑坡體上安裝4臺一體化GNSS監(jiān)測站以及3臺拉線式裂縫計(jì)，對滑坡的位移狀況進(jìn)行監(jiān)測。本次搜集該地面監(jiān)測儀器從2019年8月1日—2020年7月3日的地面監(jiān)測數(shù)據(jù)，覆蓋InSAR數(shù)據(jù)的后半段，對SBAS-InSAR解譯結(jié)果進(jìn)行精度驗(yàn)證。地面監(jiān)測儀器安裝位置如圖3所示。

1.2.3 降雨量數(shù)據(jù)

本次實(shí)驗(yàn)還搜集了滑坡周邊日降雨量數(shù)據(jù)，該數(shù)據(jù)的采集地點(diǎn)位于滑坡威脅區(qū)外，由雨量計(jì)測得。數(shù)據(jù)時(shí)間為2019年8月1日—2020年7月3日，同樣覆蓋InSAR數(shù)據(jù)的后半段。

2 數(shù)據(jù)處理

2.1 時(shí)序InSAR形變監(jiān)測原理

InSAR技術(shù)能夠獲取大面積、高精度的數(shù)字高程模型，可用來監(jiān)測地震、火山、地面沉降等造成的微小地表形變信息。用于形變分析時(shí)，InSAR技術(shù)可獲得監(jiān)測區(qū)域厘米級甚至毫米級的形變量，因此許多研究者將其用于滑坡的緩慢形變監(jiān)測，分析滑坡的滑動趨勢，為滑坡有效治理提供可靠的技術(shù)依據(jù)。

在滑坡監(jiān)測領(lǐng)域，最早使用的是D-InSAR技術(shù)，它獲取同一研究區(qū)域不同時(shí)間獲取的三幅或以上的影像進(jìn)行干涉處理，獲得干涉像對，再對干涉像對進(jìn)行差分處理，消除地形起伏和幾何畸變影響，最終獲得該時(shí)段的滑坡形變信息。其原理如下：設(shè)形變前衛(wèi)星在A1、A2處獲取影像，斜距為r1、r2,θ為衛(wèi)星在A1位置成像的側(cè)視角；形變發(fā)生后在A3位置獲取影像，得到斜距r3、空間基線B2,δ1、δ2為基線B1、B2的基線傾角，B‖、B⊥和B′‖、B′⊥分別為空間基線B1、B2沿雷達(dá)視線和垂直于雷達(dá)視線方向上的分量，ΔRd為地面點(diǎn)的形變量(圖4)。則D-InSAR形變監(jiān)測原理可表達(dá)為

圖4 D-InSAR原理Fig.4 Principle of D-InSAR

r1r2?B1sin(θ-δ1)?B‖

(1)

由式(1)可得，干涉測量相位差與視線方向基線分量成正比。

將在位置A1和A2處獲得的地表未形變前的影像進(jìn)行干涉處理，相位僅包含地形分量，相位差為

(2)

式(2)中：λ為波長；P為重復(fù)軌道，P=2。

地表發(fā)生形變后在位置A3獲取影像，將其與未發(fā)生形變時(shí)在A1處獲取的影像進(jìn)行干涉處理，得到相位包含地形和形變分量，相位差為

(3)

設(shè)視線方向的形變量為ΔRLOS,則由式(1)～式(3)，可得形變量引起的干涉圖相位差為

(4)

式(4)中：ΔφLOS可根據(jù)相位解纏獲得，則根據(jù)式(1)～式(4)可計(jì)算沿雷達(dá)視線方向的形變量ΔRLOS。

由式(4)可得，當(dāng)ΔRLOS=λ/2時(shí)，ΔφLOS=2π，因此當(dāng)形變量為半波長時(shí)，即可引起一個周期的相位變化，說明差分相位對地表形變較為敏感，可用于探測微小形變。

由D-InSAR形變監(jiān)測原理可知，該技術(shù)僅僅使用三期數(shù)據(jù)進(jìn)行差分處理，適用于分析緩慢形變滑坡較短時(shí)間周期內(nèi)的形變，當(dāng)滑坡形變量較大和監(jiān)測周期較長時(shí)，使用D-InSAR技術(shù)會導(dǎo)致空間失相干和時(shí)間失相干，降低形變監(jiān)測精度。為了解決以上問題，研究者們提出了長時(shí)間序列InSAR技術(shù)，主要包括PS-InSAR技術(shù)和SBAS-InSAR技術(shù)，這兩種方法對時(shí)空失相干和大氣擾動兩方面進(jìn)行限制，顯著提高了形變監(jiān)測能力。其中PS-InSAR利用時(shí)間相干原理，對影像數(shù)量要求高，適用于城區(qū)等獲取時(shí)間規(guī)律、連續(xù)或干涉條件和輻射比較穩(wěn)定的區(qū)域；SBAS-InSAR技術(shù)利用空間相干原理，對數(shù)據(jù)數(shù)量要求低。

由于本文中研究區(qū)域面積較大、地形起伏大、滑坡持續(xù)時(shí)間長，干涉條件不穩(wěn)定，因此選用SBAS-InSAR進(jìn)行形變監(jiān)測和分析。

2.2 SBAS處理流程

SBAS技術(shù)是短基線差分干涉方法，將所有的SAR影像依據(jù)空間和時(shí)間基線分成不同的短基線子集，各子集內(nèi)進(jìn)行差分干涉，再對纏繞的干涉相位進(jìn)行解纏，根據(jù)各相干像元的相位與觀測時(shí)間的關(guān)系，利用奇異值分解方法鏈接各差分干涉圖，抑制DEM(digital elevation model)和大氣相位對形變相位的影響，獲得最小二乘解，其處理流程見圖5。

圖5 SBAS-InSAR處理流程Fig.5 Data processing flow of SBAS-InSAR

3 研究結(jié)果與分析

3.1 地表形變速率與累計(jì)形變量

3.1.1 龐村滑坡及其周邊區(qū)域形變情況

搜集Sentinel-1A降軌數(shù)據(jù)75景，基于GAMMA軟件處理平臺，使用SBAS-InSAR數(shù)據(jù)處理方法，獲得了龐村滑坡及其周邊區(qū)域的雷達(dá)視線(line of sight, LOS)方向形變速率圖(圖6)，其中負(fù)值表示形變方向?yàn)檠乩走_(dá)視線并遠(yuǎn)離衛(wèi)星方向。在 2018年1月3日—2020年7月3日期間，龐村及其周邊區(qū)域最大形變速率約為-128.09 mm/年。

圖6 龐村滑坡及其周邊地區(qū)地表形變速率圖Fig.6 The surface deformation rate of Pangcun landslide and the surrounding area

3.1.2 地表形變速率與累計(jì)形變量

選取龐村滑坡處的解譯結(jié)果進(jìn)行詳細(xì)分析，結(jié)果如圖7所示，在2018年1月3日—2020年7月3日期間，龐村滑坡區(qū)域的最大形變速率為-38.182 mm/年。從解譯結(jié)果中可以看出滑坡形變主要發(fā)生在右側(cè)；在形變區(qū)域存在明顯的形變分區(qū)，形變速率較大的地區(qū)主要位于滑坡右側(cè)的中下部，形變速率大于23.55 mm/年，其中形變速率最大的區(qū)域主要有三處，形變速率均大于33.340 mm/年(圖7)。

依據(jù)形變速率，對滑坡區(qū)域2018年1月3日—2020年7月3日的累計(jì)形變量進(jìn)行統(tǒng)計(jì)分析，可得滑坡區(qū)域內(nèi)最大累計(jì)形變量為110.393 mm(圖8)。

圖8 龐材滑坡地表累計(jì)形變量圖Fig.8 The cumulative displacement of Pangcun landslide

3.2 滑坡形變空間分布特征

根據(jù)形變速率以及累計(jì)形變量圖，可將滑坡體劃分為3個不同形變程度的分區(qū)，如圖9所示。

圖9 龐材滑坡地表滑坡形變分區(qū)Fig.9 The deformation zone of Pangcun landslide

為進(jìn)一步研究滑坡體形變的空間分布特點(diǎn)，根據(jù)初步劃定的形變分區(qū)，沿滑坡主滑方向取兩條剖面線A—A′和B—B′(圖10),長度分別為895 m和876 m，從起點(diǎn)開始，每隔25 m取剖面點(diǎn)進(jìn)行形變速率和累計(jì)形變量分析，結(jié)果如圖11所示。

圖10 剖面線分布圖Fig.10 Section line distribution

圖11 剖面各點(diǎn)形變速率與累計(jì)形變量Fig.11 Deformation rate and cumulative displacement at each point of the section line

在剖面線A—A′上，形變速率和累計(jì)形變量最大的點(diǎn)距離起點(diǎn)600 m，形變速率達(dá)26.036 mm/年，累計(jì)形變量73.54 mm；剖面線B—B′上，形變速率和累計(jì)形變量最大的點(diǎn)位于距起點(diǎn)500 m處，形變速率和累計(jì)形變量分別為28.786 mm/年，累計(jì)形變量80.086 mm。

由剖面點(diǎn)形變速率和累計(jì)形變量分析結(jié)果可知，剖面上形變最大的點(diǎn)位于滑坡前緣區(qū)域，與實(shí)地調(diào)查結(jié)果一致，該區(qū)域主要為G219國道改造，人工修路開挖造成邊坡失穩(wěn)，形成陡坎并下滑。

3.3 滑坡形變時(shí)間分布特征

根據(jù)滑坡形變分區(qū)結(jié)果，在重點(diǎn)形變區(qū)(Ⅳ、Ⅴ、Ⅵ表示)內(nèi)選取兩個特征點(diǎn)P1、P2(圖12)進(jìn)行時(shí)序分析，研究滑坡形變的時(shí)間變化特點(diǎn)，結(jié)果如圖13所示。

圖12 重點(diǎn)形變區(qū)及特征點(diǎn)分布Fig.12 Distribution of the main deformation areas and feature points

由圖13可以看出，滑坡在兩個時(shí)段存在明顯的形變加速，分別為2018年11月29日—2019年1月22日(Ⅰ期)和2019年10月1日—2019年11月18日(Ⅱ期)，其余時(shí)期為均勻沉降。

根據(jù)實(shí)地調(diào)查結(jié)果以及降雨量數(shù)據(jù)：Ⅰ期的形變加速主要是受2018年10月對G219國道的改造工作以及隆子縣2018年12月1日地震(28.12°N， 92.62°E，震源深度5 km，據(jù)龐村滑坡23 km)影響，國道G219改造工作對滑坡前緣邊坡開挖，造成邊坡失穩(wěn)，并在地震的作用下加速下滑。Ⅱ期形變加速主要受2019年7月19日錯那5.6級地震(27.67°N， 92.89°E，震源深度10 km，據(jù)龐村滑坡72 km)以及研究區(qū)域連續(xù)降雨(圖13)，降低了巖土體的穩(wěn)定性，滑坡體加速形變。

圖13 特征點(diǎn)累計(jì)形變量時(shí)序圖Fig.13 Time series diagram of cumulative displacement of feature points

3.4 精度驗(yàn)證

精度驗(yàn)證部分，將本文中SBAS-InSAR解譯結(jié)果與地面監(jiān)測結(jié)果進(jìn)行對比，得到如下結(jié)果。

(1)SBAS與GNSS監(jiān)測結(jié)果。如圖14所示，將地面GNSS監(jiān)測站監(jiān)測結(jié)果與SBAS在相應(yīng)區(qū)域的解譯結(jié)果進(jìn)行對比發(fā)現(xiàn)：四個區(qū)域均處于持續(xù)形變階段，并且儀器GNSS-4處的累計(jì)形變量比其他三處小，SBAS解譯結(jié)果可得到相同結(jié)論。

圖14 SBAS-InSAR與GNSS結(jié)果對比Fig.14 Comparison of SBAS-InSAR and GNSS results

(2)SBAS與拉線式位移計(jì)監(jiān)測結(jié)果。如圖15所示，將地面拉線式位移計(jì)監(jiān)測結(jié)果與SBAS在相應(yīng)區(qū)域的解譯結(jié)果進(jìn)行對比發(fā)現(xiàn)：LF-1與LF-2兩個區(qū)域處于持續(xù)形變階段，LF-3在前期形變后處于穩(wěn)定階段，SBAS解譯結(jié)果可得到相同結(jié)論；同時(shí)結(jié)合監(jiān)測點(diǎn)位分布圖(圖3：LF-1LF-2處于滑坡后緣右側(cè)，LF-3處于滑坡后緣左側(cè))，驗(yàn)證了SBAS解譯結(jié)果以及根據(jù)解譯結(jié)果劃分的形變分區(qū)與地面拉線式位移監(jiān)測結(jié)果的一致性，即滑坡形變主要集中在滑坡右側(cè)。

圖15 SBAS-InSAR與拉線式位移計(jì)結(jié)果對比Fig.15 Comparison of SBAS-InSAR and linear displacement results

由兩種對比結(jié)果可以看出，地面監(jiān)測結(jié)果在數(shù)值上比SBAS解譯結(jié)果大，這是由于地面監(jiān)測的主要為滑坡滑動方向的形變量，而SBAS解譯結(jié)果為沿著雷達(dá)視線方向的形變量；并且由于Sentinel-1A數(shù)據(jù)的地面分辨率較低，區(qū)域內(nèi)形變結(jié)果代表了該對應(yīng)像素分辨率內(nèi)地物的平均形變量值。但二者在總體形變趨勢、形變的空間分布以及時(shí)間分布特征上基本一致。

因此運(yùn)用SBAS-InSAR對免費(fèi)的Sentinel-1數(shù)據(jù)進(jìn)行解譯處理，所得的滑坡形變結(jié)果可為滑坡形變監(jiān)測與治理提供可靠的技術(shù)依據(jù)；并且相對于地面監(jiān)測站來說，成本更低，結(jié)果較為準(zhǔn)確，可在部分無地面監(jiān)測站點(diǎn)的大型滑坡監(jiān)測中，使用SBAS解譯結(jié)果代替地面監(jiān)測結(jié)果。

4 結(jié)論

以西藏自治區(qū)龐村滑坡為研究對象，采用SBAS-InSAR技術(shù)進(jìn)行特大型滑坡形變監(jiān)測，獲得了研究區(qū)域2018年1月3日—2020年7月3日時(shí)間段內(nèi)的地表形變速率與累計(jì)形變量圖，并對結(jié)果進(jìn)行分析，得出了以下結(jié)論。

(1)在2018年1月3日—2020年7月3日期間，滑坡區(qū)域的最大形變速率為-38.182 mm/年?；滦巫冎饕l(fā)生在滑坡右側(cè)中下部。

(2)在2018年1月3日—2020年7月3日期間，滑坡區(qū)域內(nèi)最大累計(jì)形變量為110.393 mm，滑坡區(qū)域存在明顯的分界線，可分為三個重點(diǎn)形變區(qū)。

(3)空間分布上，剖面線A—A′和B—B′的起點(diǎn)和終點(diǎn)均不存在形變，形變最大的區(qū)域分別位于距離剖面A—A′和B—B′起點(diǎn)600 m和500 m處。

(4)時(shí)間分布上，在2018年1月3日—2020年7月3日期間，滑坡表面存在兩個明顯的形變加速時(shí)段，分別為2018年12月29日—2019年1月22日和2019年10月1日—11月18日，其余時(shí)期為均勻沉降。

(5)精度驗(yàn)證上，本文中SBAS-InSAR解譯結(jié)果與地面監(jiān)測所獲得的形變結(jié)果基本一致，驗(yàn)證了使用免費(fèi)的Sentinel-1A數(shù)據(jù)，基于SBAS-InSAR技術(shù)對大型滑坡進(jìn)行形變監(jiān)測的有效性。

因此，使用免費(fèi)的Sentinel-1A數(shù)據(jù)，基于SBAS-InSAR技術(shù)對西藏自治區(qū)龐村特大型滑坡進(jìn)行形變監(jiān)測與分析，充分展示了該滑坡形變的時(shí)間和空間分布特征，為探究滑坡的滑動機(jī)理以及后續(xù)滑坡治理工程提供了可靠的技術(shù)依據(jù)；將形變監(jiān)測結(jié)果與地面調(diào)查以及地面監(jiān)測結(jié)果進(jìn)行對比，驗(yàn)證了解譯結(jié)果的可靠性，可將監(jiān)測結(jié)果進(jìn)行由單點(diǎn)到面域的轉(zhuǎn)換；并且在部分無地面監(jiān)測站點(diǎn)的滑坡災(zāi)害中，可將SBAS-InSAR監(jiān)測結(jié)果在一定程度上代替地面監(jiān)測結(jié)果，為后續(xù)類似的特大型滑坡監(jiān)測提供新的思路和參考。