大渡河流域時序InSAR地質災害監測

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

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

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

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

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

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

1 研究區與數據介紹

1.1 研究區概況

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

1.2 數據介紹

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

表1 研究區SAR影像數據詳細參數Table 1 SAR image description used in this study

2 研究方法

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

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

(1)

(2)

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

(3)

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

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

(4)

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

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

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

3 結果分析

3.1 大渡河流域總體結果分析

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

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

3.2 隱患區識別與分析

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

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

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

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

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

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

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

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

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

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

4 結論

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

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