巨淑君, 高志良, 蔣明成, 高強

(1. 成都理工大學(xué)環(huán)境與土木工程學(xué)院, 成都 610059; 2.國能大渡河流域水電開發(fā)有限公司, 成都 610095; 3.大連理工大學(xué)建設(shè)工程學(xué)部, 大連 116081)

大渡河流域位于中國西南山區(qū),處在青藏高原邊緣向四川盆地的過渡地帶[1],多為高山峽谷地貌,自然地質(zhì)條件錯綜復(fù)雜,蘊含大量活躍的地震帶,地質(zhì)結(jié)構(gòu)極其脆弱,滑坡、泥石流、地震等災(zāi)害頻發(fā),近十多年來就發(fā)生過多起重大地質(zhì)災(zāi)害事故,較典型的有2008年的汶川地震、2017年九寨溝地震和2022年的瀘定地震[2],并伴隨有大量的滑坡等次生災(zāi)害[3-6]。據(jù)統(tǒng)計,大渡河地區(qū)已發(fā)現(xiàn)的滑坡、泥石流等地質(zhì)災(zāi)害隱患點超過了5 000多處[7-8],作為中國重要的水電資源開發(fā)基地,廣泛分布的地災(zāi)隱患區(qū)極大威脅了水電設(shè)施和人員安全[9],制約了流域內(nèi)社會經(jīng)濟發(fā)展。對大渡河地區(qū)開展大范圍地災(zāi)普查和早期識別,可以最大程度地減小災(zāi)害造成的損失,對于區(qū)域的防災(zāi)減災(zāi)工作具有重要意義。

目前,該地區(qū)的地質(zhì)災(zāi)害識別手段主要依賴于地面測量方法,采取人工排查、無人機攝影測量和GNSS (global navigation satellite system)基站監(jiān)測等技術(shù),對重點路段和點位進行監(jiān)視。然而,由于流域地質(zhì)環(huán)境惡劣,植被覆蓋茂密,這些隱患點具有隱蔽性強、分布范圍廣、調(diào)查難度大等特點,給地面監(jiān)測手段帶來極大挑戰(zhàn),當(dāng)前的技術(shù)手段無法滿足大范圍、長時間尺度的地質(zhì)災(zāi)害監(jiān)測識別。隨著空間遙感技術(shù)的不斷發(fā)展,合成孔徑雷達干涉測量(interferometric synthetic aperture radar, InSAR)技術(shù)以其具有高空間分辨率、高測量精度、長期監(jiān)測等特點,在地質(zhì)測量領(lǐng)域受到廣泛關(guān)注。InSAR技術(shù)通過分析一組多時相SAR影像,可以從中提取毫米級精度的地表形變信息,以較低成本實現(xiàn)對觀測區(qū)域的大規(guī)模、全天候、高頻次測量,特別是在大范圍地質(zhì)災(zāi)害監(jiān)測領(lǐng)域表現(xiàn)出極大的應(yīng)用潛力[10],被廣泛應(yīng)用在中國西南復(fù)雜地形區(qū)域的災(zāi)害監(jiān)測領(lǐng)域。文獻[11]利用DInSAR (differential InSAR)和時序InSAR技術(shù)分析了九寨溝地震后的區(qū)域滑坡情況;文獻[12]采用Stacking技術(shù)對金沙江地區(qū)滑坡群進行識別;文獻[13-17]采用了使用最為普遍的SBAS(small baseline subsets)技術(shù)對西南山區(qū)開展了大范圍的滑坡災(zāi)害研究,取得了較多成果。然而SBAS技術(shù)的多視處理會降低分辨率,尤其是在使用Sentinel-1這類中等分辨率數(shù)據(jù),可能會丟失許多局部滑坡區(qū)域,可以結(jié)合PSI(permanent scatterer InSAR)技術(shù)和SBAS技術(shù)特點,應(yīng)用于滑坡監(jiān)測中。

然而當(dāng)前對大渡河流域滑坡研究主要集中在丹巴、甲居和瀘定等個別隱患區(qū),缺少對大渡河全流域的滑坡災(zāi)害普查。基于此,現(xiàn)針對大渡河流域復(fù)雜地形和氣候條件下的地災(zāi)隱患監(jiān)測開展研究,探索利用時序InSAR技術(shù)開展大渡河流域地質(zhì)災(zāi)害全面、定期排查的可行性。收集了2018—2021年的Sentinel-1超過1 TB的數(shù)據(jù),采用PS-SBAS全流程自動化數(shù)據(jù)處理技術(shù)對海量SAR影像數(shù)據(jù)進行了深入分析,實現(xiàn)對大渡河流域丹巴至樂山段完整區(qū)域大范圍地表形變監(jiān)測,首次獲取了大渡河全流域歷史災(zāi)害隱患分布結(jié)果。

此外,在過去的研究基礎(chǔ)上,在該流域又識別出了10多處新的滑坡隱患區(qū)域,并對流域不同區(qū)域滑坡誘因進行了總結(jié)。本文成果可以為流域內(nèi)大壩和邊坡變形機理分析和災(zāi)害預(yù)警提供可靠數(shù)據(jù)和信息支撐,有效提升該地區(qū)災(zāi)害治理水平,同時方便后續(xù)對比分析瀘定地震前后的流域地質(zhì)災(zāi)害情況。

1 研究區(qū)與數(shù)據(jù)介紹

1.1 研究區(qū)概況

研究區(qū)位于四川盆地南部的大渡河流域,由北向南流經(jīng)金川、丹巴、瀘定等地,至石棉折向東,構(gòu)成“L”字形河流走向,整個流域狹長地帶共1 062 km。區(qū)域內(nèi)地形起伏極大,相對高差達到2 000～3 000 m,岸坡穩(wěn)定性差,遍布高陡邊坡與危巖體,是中國地質(zhì)災(zāi)害重災(zāi)區(qū)。隨著該地區(qū)經(jīng)濟快速發(fā)展,大規(guī)模工程建設(shè)活動的增加,進一步加劇災(zāi)害發(fā)生頻率。各類突發(fā)性的滑坡等災(zāi)害,不僅可能引發(fā)河流堵塞,同時也會阻斷流域內(nèi)唯一的通道S211省道,威脅水電站安全和區(qū)域施工建設(shè)活動,因此大渡河地區(qū)對于大范圍地質(zhì)災(zāi)害監(jiān)測需求較為迫切。本文選取了大渡河丹巴至樂山段作為研究區(qū),覆蓋了大渡河流域主要區(qū)域,如圖1(a)所示,其中藍線表示大渡河;圖1(b)為丹巴附近的光學(xué)影像圖,可看到區(qū)域谷深坡陡,是典型的高山峽谷地貌。

1.2 數(shù)據(jù)介紹

目前能夠滿足大渡河流域形變監(jiān)測應(yīng)用的SAR影像主要是國外的Sentinel-1A數(shù)據(jù),該衛(wèi)星數(shù)據(jù)是由歐空局提供的C波段(波長5.6 cm)中等分辨率數(shù)據(jù),數(shù)據(jù)來源(https://scihub.copernicus.eu/)。按軌道號排序為26-1、26-2和128-1,如圖2紅色框所示。每個片區(qū)獲取了102景SAR影像,影像拍攝時間為2018年1月—2021年6月,極化方式為VV(vertical transmit, vertical receive)極化,TOPS(terrain observation with progressive scans)成像模式IW(interferometric wide swath)類型數(shù)據(jù),超過300景的影像總數(shù)據(jù)量達到了1.25 TB,具體參數(shù)如表1所示。此外,還需要下載Sentinel-1數(shù)據(jù)對應(yīng)影像時間的POD(precise orbit determination)精軌文件,獲取地址為(https: //s1qc.asf.alaska.edu/aux_poeorb/),以修正軌道參數(shù)信息,以及時序InSAR處理過程需要使用的外部高程數(shù)據(jù),采用SRTM(shuttle radar topography mission) 30 m DEM(digital elevation model),數(shù)據(jù)來源(http://srtm.csi.cgiar.org/srtmdata/)。

表1 研究區(qū)SAR影像數(shù)據(jù)詳細參數(shù)Table 1 SAR image description used in this study

2 研究方法

傳統(tǒng)的基于永久散射體PSI技術(shù)在大渡河區(qū)域,在時間基線較長時會出現(xiàn)嚴重的時間去相干誤差,導(dǎo)致較大形變解算誤差;而小基線集SBAS技術(shù)在復(fù)雜的西南山區(qū),大量的陰影和疊掩區(qū)域的存在,會導(dǎo)致形變估計結(jié)果出現(xiàn)較大的分塊跳變誤差。因此,為了研究2018—2021年較長時間跨度的大渡河流域的地表形變信息,采用PSI和SBAS相結(jié)合的PS-SBAS方法,通過選取永久散射體PS,并基于小基線干涉對組合獲取PS點對應(yīng)的差分干涉相位,有效避免了長時間基線干涉對,同時不受陰影等低相干區(qū)域影響,最后采用SBAS處理思路完成數(shù)據(jù)處理,基本原理如下。

對一組獲取時間為{t1,t2,…,tN}的SLC(single look complex)數(shù)據(jù),基于幅度穩(wěn)定性方法識別出M個具有穩(wěn)定散射特性的PS點,{PS1,PS2,…,PSM}。根據(jù)區(qū)域干涉圖質(zhì)量,選擇合適的時間基線和空間閾值,同時可以依據(jù)相干系數(shù)閾值,篩選得到一組具有K個干涉對的小基線集組合,并對PS點的復(fù)數(shù)據(jù)進行干涉和差分干涉處理,得到去除平地和地形相位的PS點差分干涉圖,經(jīng)過空域二維相位解纏后得到解纏干涉相位,此時由影像獲取時間為{ti,tj}組成的解纏干涉圖,第m個PS點解纏相位由以下幾個分量組成,即

(1)

(2)

式(2)中:hm表示該點的殘余高程;B⊥為空間垂直基線;λ、r、θ對應(yīng)雷達參數(shù)中的波長、斜距和入射角。形變相位為時間t的線性函數(shù),可以根據(jù)小基線集組合相位時間序列,得到相位變化速率,即

(3)

式(3)中:φk為對應(yīng)時間的形變相位。組合每個干涉相位觀測值,可以得到線性方程組為

[Δφ1,Δφ2,…,ΔφK]T

(4)

式(4)中:T為系數(shù)矩陣,每行只有1和-1兩個非0值,位置一一對應(yīng)于干涉對組合;Δφ為PS點差分解纏相位。對式(4)利用奇異值分解SVD(singular value decomposition)算法,可以解算得到形變速率和殘余高程值的最小二乘解,最后疊加上經(jīng)過大氣濾波后得到的非線性形變分量,得到最終的形變結(jié)果。

由于完整覆蓋流域的Sentinel-1影像數(shù)據(jù)量較大,傳統(tǒng)完全依賴人工的數(shù)據(jù)處理方法,需要耗費大量的人力和時間代價。為了應(yīng)對InSAR大數(shù)據(jù)處理的挑戰(zhàn),根據(jù)以往的數(shù)據(jù)處理經(jīng)驗,利用適合流域的最優(yōu)InSAR處理參數(shù)構(gòu)建了全自動化PS-SBAS數(shù)據(jù)處理流程,在干涉處理流程,能夠根據(jù)相干性閾值選取合適數(shù)量的干涉對,同時在時序分析過程,能夠自適應(yīng)地完成多次迭代形變參數(shù)解算,從而極大提升了數(shù)據(jù)處理效率。具體數(shù)據(jù)處理過程如下:以26-1區(qū)域數(shù)據(jù)為例,選取獲取時間為2019-12-04的SAR影像作為主參考影像數(shù)據(jù),將其他輔圖像配準到主影像參考坐標(biāo)下,采用針對TOPS模式數(shù)據(jù)的配準算法,得到配準好的SLC數(shù)據(jù)?；诜入x差法選取PS點后,以時間基線100 d、空間基線150 m作為閾值獲取干涉對組合,并利用30 m SRTM數(shù)據(jù)生成差分干涉圖。待利用MCF(minimum cost flow)(最小費用流)算法進行相位解纏后,對形變和高程參數(shù)進行估計,并根據(jù)殘余相位質(zhì)量迭代多次解算,得到線性形變速率和殘余高程,最后利用SVD算法進行解算,并經(jīng)過大氣濾波后,得到最終的形變時間序列,同時依據(jù)時序相干系數(shù)閾值0.7精選PS點,地理編碼后得到shapefile格式的形變結(jié)果文件。完整的數(shù)據(jù)處理流程如圖3所示。

圖3 PS-SBAS處理流程圖Fig.3 The flowchart of PS-SBAS method

3 結(jié)果分析

3.1 大渡河流域總體結(jié)果分析

采用第2節(jié)所述的PS-SBAS處理流程,對區(qū)域26-1、26-2和128-1這三個片區(qū)2018—2021年Sentinel-1數(shù)據(jù)進行了處理分析,獲取了大渡河流域丹巴至樂山段完整的大范圍InSAR形變測量結(jié)果,由于關(guān)注重點是大渡河流域以及S211省道沿線兩側(cè)的區(qū)域,截取了大渡河流域沿岸1 km緩沖區(qū)的形變結(jié)果,其時空分布如圖4所示。

圖4中每個點都對應(yīng)于SAR影像上選取的高相干散射點,形變解算結(jié)果均為雷達視線方向上的(line of sight, LOS)位移速率(mm/year),并根據(jù)其年形變速率進行顏色編碼,其正負代表是遠離或靠近衛(wèi)星方向。時序InSAR技術(shù)獲取到總的點密度分布在106個/km2,總體形變速率范圍為-65.2～51.1 mm/y。根據(jù)形變速率大小以及時序變化曲線,從總體上看,大部分區(qū)域地表形變都比較小或接近0,僅有少部分形變較大的區(qū)域主要聚集在丹巴、石棉和漢源等地。為了方便解譯,在圖4中對這幾個區(qū)域結(jié)果進行了局部放大,從圖4中可以觀測到,沿著大渡河存在多處明顯沉降區(qū)域(黃色和紅色點)。根據(jù)歷史滑坡成因分析的經(jīng)驗[18-20],流域不同區(qū)域滑坡誘發(fā)因素存在差異,其中處于大渡河上游的丹巴地區(qū),地層較為穩(wěn)定,區(qū)域滑坡最容易受到陡峭地形和河流侵蝕的影響;位于中下游的石棉和漢源地區(qū),地質(zhì)構(gòu)造活動和降雨較易誘發(fā)滑坡失穩(wěn)。因此,可以看到大渡河流域滑坡與地質(zhì)和氣候條件緊密相關(guān)。此外,還有許多沒有PS測量點覆蓋的區(qū)域,主要是由于植被覆蓋茂密,以及坡體朝向、地形起伏與SAR側(cè)視成像幾何問題的影響,缺少足夠多的PS測量點覆蓋,無法獲取準確的形變信息。

3.2 隱患區(qū)識別與分析

對大渡河流域InSAR結(jié)果大范圍解譯,綜合考慮形變速率大小、累積變形量、形變區(qū)域面積和坡向這幾個因素,同時要保證測量區(qū)域在SAR影像上相干性較高、測量結(jié)果較為可靠。通過結(jié)合光學(xué)圖像對比分析,最終確定了15處威脅較大、形變較嚴重的隱患區(qū)。其中多個隱患區(qū)的變形區(qū)域發(fā)育在古滑坡體附近,坡體失穩(wěn)嚴重,對邊坡下方的居民區(qū)和大渡河構(gòu)成較大威脅,需要地面測量人員重點關(guān)注,具體的隱患區(qū)分布和形變信息見表2所示。在隱患點列表中,不僅包含了以往研究發(fā)現(xiàn)的甲居村和丹巴縣等受到古滑坡威脅的區(qū)域,以及地面重點觀測的鄭家坪邊坡,還新增加了10多處大型滑坡隱患點,為開展地質(zhì)調(diào)查提供了數(shù)據(jù)支持,獲得了瀘定地震前的災(zāi)害分布情況。因此,可以看到基于衛(wèi)星InSAR技術(shù)和光學(xué)影像輔助解譯方法,極大提高了大渡河流域災(zāi)害監(jiān)測效率,對于識別一些隱蔽滑坡隱患區(qū)具有無可比擬的優(yōu)勢。

為了進一步說明隱患區(qū)域情況,同時驗證InSAR測量結(jié)果的合理性,下面針對其中三個典型形變區(qū)域進行深入分析,包括001鄉(xiāng)道滑坡隱患區(qū)、瀑布溝水電站和鄭家坪邊坡。

(1)001鄉(xiāng)道滑坡隱患區(qū)。001鄉(xiāng)道位于丹巴縣,距離丹巴水電站大約2.5 km,該隱患邊坡緊鄰大渡河,坡度較緩,從光學(xué)圖上可看到,部分區(qū)域進行了邊坡防護。提取該區(qū)域InSAR測量結(jié)果如圖5所示,可看到該不穩(wěn)定區(qū)域覆蓋面積達到11.2×104m2,區(qū)域最大形變速率達到了43.4 mm/y,從2018年1月—2021年6月最大累積沉降量超過了110 mm。提取坡面兩個點TS1、TS2的形變時間曲線序列如圖5(b)、圖5(c)所示,這兩個典型點位的形變曲線上表現(xiàn)出勻速下沉趨勢,沉降主要來源于區(qū)域內(nèi)道路施工以及工程建設(shè)活動,加劇了邊坡失穩(wěn)現(xiàn)象,InSAR測量結(jié)果與實際情況相符。整個坡面每年沉降量平均在35～45 mm,對坡底的居民區(qū)和學(xué)校構(gòu)成直接的安全隱患,同時對大渡河帶來一定威脅,其危險程度較高,需要地面測量人員持續(xù)關(guān)注。

(2)瀑布溝水電站沉降隱患區(qū)。在大渡河流域部分水庫、壩體也表現(xiàn)出較明顯沉降,瀑布溝水電站位于漢源縣和甘洛縣交界處,是該流域重要的梯級水電站之一。InSAR測量結(jié)果如圖6所示,該水電站壩體最大形變速率達-22.966 mm/y,整個大壩每年沉降量在20 mm左右,目前屬于正常沉降范圍。可以觀察到該結(jié)果與其他水壩形變具有類似的變化特征[21]。位于水壩受力最大的TS1處,其下沉速率最大,達到了23 mm/y,此外該處變形表現(xiàn)出明顯的年際變化,與水庫蓄放水過程具有強關(guān)聯(lián)性,證明了InSAR結(jié)果的合理性;隨著壩體位置變化,從TS1到TS2和TS3處,壩體沉降速率由-26 mm/y逐漸減緩至-5.8 mm/y,由勻速沉降變?yōu)椴▌邮较鲁?壩體周邊區(qū)域較為穩(wěn)定,需要對水電站壩體進行持續(xù)監(jiān)測。

(3)鄭家坪邊坡測量結(jié)果。鄭家坪滑坡體位于瀘定縣,是流域的重點關(guān)注區(qū),有地面測量人員長期監(jiān)視。從圖7所示的該區(qū)域時序InSAR處理結(jié)果中,可以看到該坡面底部較為穩(wěn)定,在坡頂處存在明顯沉降區(qū),最大沉降速率達到20.8 mm/y,區(qū)域覆蓋面積大不,緊鄰大渡河,2018—2021年累積沉降量超過75 mm,由于該地區(qū)活躍的地質(zhì)構(gòu)造運動,使得該邊坡存在一定的失穩(wěn)風(fēng)險,對坡底的公路和大渡河造成較大威脅。

表2 主要滑坡隱患區(qū)列表Table 2 The identified hidden landslide region list

圖5 001鄉(xiāng)道附近形變分布和特征點形變時間變化曲線Fig.5 The deformation pattern in Road 001, and the displacement of TS1 and TS2

為了對InSAR形變測量結(jié)果進行驗證,利用附近的地面GNSS測量結(jié)果(圖7中▲所示)進行對比分析。獲取了離InSAR測量區(qū)域最近的兩個GNSS站點ZP13和ZP14,將其GNSS測量結(jié)果投影到InSAR的衛(wèi)星視線LOS方向,為了去除兩者參考點差異,對InSAR與GNSS兩點分別差分的結(jié)果如圖8所示。

從圖8可以看到,InSAR測量結(jié)果與GNSS在整體能夠保持較好的一致性,其形變趨勢是相同的,兩者形變量曲線間最大誤差低于10 mm,主要差異是由于InSAR測量點與GNSS站點并不完全相同所致,誤差在可接受范圍,證明了該區(qū)域InSAR測量結(jié)果的準確性。

圖8 鄭家坪區(qū)域InSAR與GNSS測量結(jié)果對比圖Fig.8 The deformation results comparison between InSAR and GNSS in Zhenjiaping

4 結(jié)論

基于時序InSAR自動化處理技術(shù)完成了對大渡河流域丹巴至樂山段大范圍滑坡災(zāi)害監(jiān)測普查研究,通過對2018—2021年期間大規(guī)模的Sentinel-1數(shù)據(jù)處理,首次獲取了大渡河流域完整的滑坡隱患分布情況,構(gòu)建了流域地質(zhì)災(zāi)害快速監(jiān)測識別方法。通過綜合考慮區(qū)域形變情況,同時重點針對給水電站和人員帶來較大安全威脅,以及堵塞河流風(fēng)險較大的隱患點,在前人研究基礎(chǔ)上,從測量結(jié)果中又提取了10多處新的隱患區(qū)域,最大沉降速率在-65～-30 mm/y,對這些形變風(fēng)險區(qū)進行了詳細描述和危害評估。最后,對幾個典型隱患區(qū)域進行了深入分析和驗證,進一步核實了InSAR測量結(jié)果的合理性。

綜上所述,本文的研究證明了InSAR技術(shù)可以被用于復(fù)雜地質(zhì)環(huán)境下的大范圍地質(zhì)災(zāi)害普查,構(gòu)建的自動化InSAR處理流程可以形成常態(tài)化的形變監(jiān)測,并以較低成本快速獲取災(zāi)害分布一張圖,持續(xù)為相關(guān)部門進行防災(zāi)減災(zāi)工作提供豐富的數(shù)據(jù)支撐,也為中國其他流域的災(zāi)害監(jiān)測提供了參考范例,該測量結(jié)果將用于指導(dǎo)進一步的地面測量工作,同時對于后續(xù)對瀘定地震前后流域形變變化情況分析也具有重要參考意義。