趙超英 劉曉杰高 楊馮曉松

1.長安大學地質(zhì)工程與測繪學院,陜西 西安 710054;

2.西部礦產(chǎn)資源與地質(zhì)工程教育部重點實驗室,陜西 西安 710054;

3.自然資源部生態(tài)地質(zhì)與災(zāi)害防控重點實驗室,陜西 西安 710054;

4.中國地質(zhì)科學院地質(zhì)力學研究所,北京 100081

0 引言

雅魯藏布江是中國最長的高原河流,位于西藏自治區(qū)境內(nèi),巨大的落差以及充沛的徑流量使其成為世界上水利資源最豐富的地區(qū)之一(江平等,2022;劉麗紅等,2022)。雅魯藏布江大峽谷位于雅魯藏布江下游,由于該區(qū)域是印度板塊和歐亞板塊碰撞的前緣地帶(李翔,2019),區(qū)域內(nèi)新構(gòu)造運動活躍,地層擠壓褶皺,斷裂切割,使得地勢起伏大,高山分布眾多,河谷深切現(xiàn)象普遍,山頂與河谷高差多大于2000 m,有的相對高差甚至大于3000 m,河流下切強烈,屬典型高山深切割區(qū)(李金城,2007;祝嵩,2012;辛聰聰,2019)。由于獨特的地質(zhì)構(gòu)造以及氣候變化的影響,區(qū)域內(nèi)崩塌、滑坡、泥石流等地質(zhì)災(zāi)害頻發(fā)(王軍等,2000;袁喆,2014;劉麗紅等,2022),例如1950年墨脫地震后引發(fā)冰川冰崩,導(dǎo)致雅魯藏布江下游派鎮(zhèn)100余人死亡(張沛全等,2008);2000年4月9日,林芝地區(qū)扎木弄溝發(fā)生大規(guī)模山體滑坡,截斷了易貢藏布河,之后堰塞湖潰決形成特大洪水,沖毀了下游道路以及沿岸村莊(殷躍平,2000;劉錚等,2020);2018年10月17日雅魯藏布江下游色東普溝發(fā)生特大型泥石流堵江事件,堰塞壩潰決后又造成洪水災(zāi)害,導(dǎo)致堰塞壩周邊鄉(xiāng)鎮(zhèn)受到威脅(劉傳正等,2019)。

雅魯藏布江下游復(fù)雜特殊的地質(zhì)環(huán)境導(dǎo)致了其地質(zhì)災(zāi)害問題的復(fù)雜性和特殊性(李濱等,2020)。為減少該區(qū)域地質(zhì)災(zāi)害給當?shù)丨h(huán)境及居民造成的影響,對區(qū)域內(nèi)高位地質(zhì)災(zāi)害進行早期識別是有必要的。合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar,InSAR)是多門學科交叉發(fā)展而來的新型技術(shù)(Moreira et al., 2013)。它可以全天候、全天時地實現(xiàn)對地表變形的厘米級至亞厘米級精度觀測。由于InSAR技術(shù)的覆蓋范圍廣及監(jiān)測精度高等特點,自從問世以來,就受到了地質(zhì)災(zāi)害研究人員的廣泛關(guān)注(Shi et al., 2016; Liu et al., 2021; Dun et al., 2021; Zhao et al., 2021; Wang et al., 2022; 朱怡飛等,2022)。但是由于米林-墨脫段地形起伏劇烈,且覆蓋有茂密的植被以及冰雪,使得InSAR監(jiān)測、識別地質(zhì)災(zāi)害的能力降低(Liu et al., 2021)。并且冰川快速移動所引起的失相干現(xiàn)象會導(dǎo)致InSAR結(jié)果的漏判(Sivalingam et al., 2022)。因此,文章以雅魯藏布江下游米林-墨脫段為研究區(qū),聯(lián)合時序InSAR技術(shù)以及SAR偏移量技術(shù)來進行高山峽谷區(qū)高位地質(zhì)災(zāi)害的識別(Gu et al., 2022)。根據(jù)Sentinel-1影像和ALOS/PALSAR-2影像獲得的形變速率圖,并結(jié)合數(shù)字高程模型(Digital Elevation Model,DEM)數(shù)據(jù)獲得的坡度圖和坡向圖以及多時相光學遙感影像,在研究區(qū)共探測到260處大小不等的地質(zhì)災(zāi)害形變區(qū)。此次研究為高山峽谷區(qū)高位災(zāi)害的識別和監(jiān)測提供了技術(shù)支持,給研究區(qū)內(nèi)的災(zāi)害防治提供了重要參考。

1 研究區(qū)及影像數(shù)據(jù)概況

1.1 研究區(qū)概況

雅魯藏布江下游米林-墨脫段(圖1)位于西藏自治區(qū)東南部(藏東、藏南地區(qū)),區(qū)域內(nèi)水系為雅魯藏布江水系,主要包括帕隆藏布江、易貢藏布江、波都藏布、拉月曲等河流。研究區(qū)內(nèi)分布著廣泛冰川、積雪和凍土。主要的地形地貌為藏東高山峽谷區(qū),由于區(qū)域內(nèi)新構(gòu)造運動活躍,地層擠壓褶皺,斷裂切割,使得地勢起伏大,高山分布眾多,河谷深切現(xiàn)象普遍(邊紫璇,2021)。但受到青藏高原間歇性抬升和河谷下切的影響,在區(qū)域內(nèi)河谷區(qū)存在夷平面和階地等層狀地貌。研究區(qū)內(nèi)人類工程經(jīng)濟活動類型包括農(nóng)林牧活動、城鎮(zhèn)工程建設(shè)、交通工程、小型礦山工程、水利工程(水庫、堤防、渠道)等,人口密度小,少部分為無人區(qū)(高鵬,2010)。雅魯藏布江下游地區(qū)地質(zhì)災(zāi)害是地質(zhì)構(gòu)造、地層巖性、新構(gòu)造運動等內(nèi)部地質(zhì)特征和氣候,冰川水文、人類活動等外部誘發(fā)因素共同作用的結(jié)果;其形成機理、發(fā)生過程、方式和影響的范圍等方面與內(nèi)陸地區(qū)的地質(zhì)災(zāi)害均存在巨大差異。該區(qū)域是全球構(gòu)造應(yīng)力作用最強、隆升和剝蝕最快、新生代變質(zhì)和深熔作用最強的地區(qū),高位遠程地質(zhì)災(zāi)害鏈、高地應(yīng)力、高地溫等重大工程地質(zhì)問題突出(翟毅飛,2022)。

圖1 研究區(qū)位置及所用SAR影像空間覆蓋(位置圖據(jù)公開的SRTM DEM數(shù)據(jù)繪制)Fig.1 Location of the study area and coverage of SAR images (Location map based on publicly available SRTM DEM data)

1.2 影像數(shù)據(jù)

現(xiàn)有的SAR衛(wèi)星主要工作波段為X波段(波長3.1 cm)、C波段(5.6 cm)以及L波段(23.5 cm)。對于InSAR形變監(jiān)測而言,形變測量的靈敏度在一定程度上取決于雷達的波長,波長越短,靈敏度越高,反之波長越長,靈敏度越低(康亞,2020)。研究區(qū)處于地質(zhì)環(huán)境復(fù)雜的極高山區(qū),植被覆蓋極其茂密且大部分區(qū)域常年被冰雪覆蓋。在冰雪及植被覆蓋情況下,短波長的X波段和C波段SAR數(shù)據(jù)難以穿透植被冠層及干燥積雪,時間相干性保持能力弱(李凌婧等,2022);波長更長、穿透能力更強的L波段SAR數(shù)據(jù)往往能保持更好的相干性。在充分考慮每種波長SAR數(shù)據(jù)形變測量靈敏度及相干性保持能力的情況下,采用長短波長SAR數(shù)據(jù)相結(jié)合的方式進行研究區(qū)高位冰崩、崩塌及滑坡災(zāi)害的綜合判識與監(jiān)測研究。

共計收集了覆蓋研究區(qū)2017年3月—2020年7月284景Sentinel-1影像以及2014年9月—2020年5月111景ALOS/PALSAR-2影像,SAR數(shù) 據(jù) 基本參數(shù)見表1。采用30 m空間分辨率的SRTM DEM及日本AW3D30 DSM進行地形相位的去除及結(jié)果分析。

表1 所用SAR數(shù)據(jù)基本參數(shù)Table 1 Fundamental parameters for SAR images used in this study

2 高位地質(zhì)災(zāi)害早期識別方法概述

采用時序InSAR技術(shù)以及SAR偏移量技術(shù)來進行高山—極高山區(qū)大范圍地質(zhì)災(zāi)害識別。由于研究區(qū)域海拔高,地形起伏劇烈,監(jiān)測環(huán)境的差異性較大,并且不同的多時相InSAR/SAR技術(shù)具有各自不同的適用特點(Osmanogˇlu et al., 2016)。因此針對每個區(qū)域的環(huán)境以及存檔SAR數(shù)據(jù)使用了多種InSAR/SAR數(shù)據(jù)處理流程。在滑坡識別中所應(yīng)用的InSAR/SAR技術(shù)包括:Stacking-InSAR(劉 國 祥,2019)、SBAS-InSAR(Tizzani et al., 2007)以及SAR偏移量技術(shù)(Strozzi et al., 2007)。使用InSAR/SAR技術(shù)的主要目的在于恢復(fù)出覆蓋研究區(qū)域的大范圍地表形變。在恢復(fù)出地表形變之后,結(jié)合形變信息、DEM以及多時相光學遙感影像進行滑坡的圈定工作。具體的操作流程如圖2所示,主要分為以下幾步。

圖2 高山峽谷區(qū)高位地質(zhì)災(zāi)害SAR/InSAR識別技術(shù)流程圖Fig.2 Flow chart of the SAR/InSAR identification technology for high-elevation geohazard in a mountain-valley area

(1)SAR數(shù)據(jù)的合理獲取。根據(jù)研究區(qū)域的地形起伏程度、降雨、植被覆蓋、地質(zhì)環(huán)境等以及結(jié)合每種SAR衛(wèi)星的波長、入射角大小、飛行方向等進行SAR數(shù)據(jù)的合理選取,最大程度增加SAR數(shù)據(jù)的可視區(qū)域,減少地質(zhì)災(zāi)害的漏判與誤判(Dai et al., 2022)。

(2)設(shè)置合理的時空基線閾值生成所有可能的干涉圖,抑制時空去相干的影響(Li et al., 2022)。

(3)對所生成的所有干涉圖進行濾波、相位解纏以及誤差改正。濾波的主要目的是削弱噪聲的影響,使相位變得更加連續(xù),最終達到正確相位解纏與地質(zhì)災(zāi)害精準解譯的目的。相位解纏的主要目的是恢復(fù)出相位主值中被模糊掉的整周相位,從而計算出正確的形變信息。誤差改正包括大氣延遲誤差改正和DEM誤差改正。大氣延遲誤差包括垂直分層大氣延遲和湍流混合效應(yīng)。對于垂直分層大氣延遲,可采用數(shù)學模型以及外部大氣產(chǎn)品(例如GACOS、ECMWF等)進行改正;對于湍流混合效應(yīng),可采用時空域濾波的方式進行改正。ALOS/APLSAR-2部分影像大氣誤差改正前后解纏圖如圖3所示(Kang et al., 2021)。對于DEM誤差,首先選取具有短時間基線、長空間基線的干涉對,并對其采用最小二乘法求解DEM誤差(李偉華,2020),然后對原始DEM進行更新。

圖3 大氣誤差改正前后解纏圖Fig.3 Unwrapping images before and after the correction of atmospheric error

(4)基于干涉圖相干性以及誤差改正后每個解纏圖的殘余大氣噪聲水平來自動選取優(yōu)質(zhì)的干涉圖。對于一些觀測環(huán)境極為復(fù)雜的地區(qū),例如有劇烈的地形起伏以及極其茂密的植被覆蓋,往往會造成相位的不連續(xù)以及不一致,從而引起較為嚴重的相位解纏誤差。因此,在自動選取干涉圖基礎(chǔ)上通過目視方式進行解纏誤差的檢查與改正,并進一步選取高質(zhì)量的解纏圖進行形變的反演。對于解纏誤差的改正,可采用相位閉合或者相位補償?shù)姆绞竭M行改正(王霞迎,2018)。ALOS/PALSAR-2部分影像相位解纏誤差改正前后解纏圖如圖4所示。

圖4 相位解纏誤差改正前后解纏圖Fig.4 Unwrapping images before and after the correction of phase error

(5)采用多種SAR/InSAR技術(shù)獲取研究區(qū)域地表形變圖。Stacking-InSAR、SBAS-InSAR與偏移量等SAR/InSAR技術(shù)被用來反演研究區(qū)域的地表形變圖。

3 雅魯藏布江下游高位地質(zhì)災(zāi)害早期識別

3.1 形變結(jié)果分析

采用多種SAR/InSAR技術(shù)包括Stacking-InSAR、SBAS-InSAR與Offset-tracking以及光學偏移量技術(shù)對所獲得的SAR數(shù)據(jù)及光學影像進行了處 理。研 究 區(qū)2017年3月 至2020年7月Sentinel-1數(shù)據(jù)地表形變速率如圖5所示,研究區(qū)2014年9月至2020年5月ALOS/PALSAR-2數(shù) 據(jù)地表形變速率如圖6所示。對于Sentinel-1數(shù)據(jù)得到的地表形變速率,大部分點目標的形變速率分布在-10～20 mm/a;對于ALOS/PALSAR-2數(shù) 據(jù)得到的地表形變速率,大部分點目標的形變速率分布在-10～20 mm/a,說明研究區(qū)整體穩(wěn)定性較好。在整個研究區(qū)中,在海拔較高的山谷及溝道探測到大量的形變聚集區(qū),主要是由于冰川運動、冰崩、滑坡以及崩塌引起的變形。林芝市巴宜區(qū)附近部分區(qū)域形變信號主要是由于誤差所引起的(圖5),該地區(qū)植被覆蓋極其茂密加之冰雪的覆蓋導(dǎo)致C波段的Sentinel-1數(shù)據(jù)相干性極低,受失相干噪聲影響較大。此外,海拔較高的部分山峰由于快速的冰川運動、冰崩以及巖崩引起了相位測量失相干,未獲得有效的觀測點目標(圖5,圖6)。對于相位測量失相干的區(qū)域,采用SAR偏移量技術(shù)及光學偏移量技術(shù)進行形變的探測。

圖5 米林-墨脫段2017年3月至2020年7月Sentinel-1數(shù)據(jù)地表形變速率Fig.5 Surface deformation rate of the Sentinel-1 data from March 2017 to July 2020 for the Milin-Motuo section

圖6 米林-墨脫段2014年9月至2020年5月ALOS/PALSAR-2數(shù)據(jù)地表形變速率Fig.6 Surface deformation rate of the ALOS/PALSAR-2 data from September 2014 to May 2020 for the Milin-Motuo section

3.2 地質(zhì)災(zāi)害編目與分析

基于Sentinel-1數(shù)據(jù)及ALOS/PALSAR-2數(shù)據(jù)相位測量獲得的地表形速率及時間序列、ALOS/PALSAR-2數(shù)據(jù)SAR偏移量測量獲得的方位向與視線向形變以及光學影像偏移量測量獲得的南北向與東西向形變,并結(jié)合DEM數(shù)據(jù)獲得的坡度圖和坡向圖以及多時相光學遙感影像,在米林-墨脫段共探測到260處大小不等地質(zhì)災(zāi)害變形區(qū),其中崩滑災(zāi)害111處,冰川、冰崩變形區(qū)149處。探測到的崩滑災(zāi)害以及冰川、冰崩變形區(qū)編目圖見圖7。

圖7 米林-墨脫段地質(zhì)災(zāi)害編目圖Fig.7 Catalogued map of geohazards in the Milin-Metuo section

為進一步分析探測到的崩滑及冰川、冰崩變形區(qū)空間上分布規(guī)律,分別繪制并統(tǒng)計了距離雅魯藏布江5 km、10 km、15 km和20 km崩滑及冰川、冰崩變形區(qū)分布。距離雅魯藏布江不同距離崩滑及冰川、冰崩變形區(qū)分布和統(tǒng)計結(jié)果分別如圖8和圖9所示。結(jié)果表明:距離雅魯藏布江5 km范圍內(nèi),分布有崩滑災(zāi)害34處,冰川、冰崩變形區(qū)15處;距離雅魯藏布江10 km范圍內(nèi),分布有崩滑災(zāi)害70處,冰川、冰崩變形區(qū)97處;距離雅魯藏布江15 km范圍內(nèi),分布有崩滑災(zāi)害91處,冰川、冰崩變形區(qū)114處;距離雅魯藏布江20 km范圍內(nèi),分布有崩滑災(zāi)害100處,冰川、冰崩變形區(qū)120處。崩滑以及冰川、冰崩變形區(qū)數(shù)量隨著離雅魯藏布江距離的增加而增加(圖9),說明探測到的大多數(shù)崩滑及冰川、冰崩變形區(qū)分布在距離雅魯藏布江較遠、海拔較高的溝道及山峰。

圖8 距離雅魯藏布江不同距離地質(zhì)災(zāi)害分布Fig.8 Distribution of geohazards at different distances from the Yarlung Zangbo River

圖9 距離雅魯藏布江不同距離地質(zhì)災(zāi)害分布統(tǒng)計Fig.9 Distribution statistics of geohazards at different distances from the Yarlung Zangbo River

在探測到的149處冰川、冰崩變形區(qū)中,有13處為大冰川的運動,其位置分布如圖10所示。在這些大冰川中,流速最快的為培龍貢支大冰川,該冰川2019年8月29日至2019年9月26日雷達方位向與視線向二維地表形變見圖11,28天時間里最大位移在雷達方位向超過10 m,在視線向超過14 m。最長的冰川為果登冰川,其長度達到16 km,該冰川2020年1月至4月南北向與東西向形變見圖12。在澤巴隆巴冰川溝中探測到一危險巖崩隱患點,其2016年6月15日至2018年3月7日雷達方位向與視線向二維地表形變見圖13。在不到兩年的時間里,最大累積形變在方位向達到-13 m,在視線向達到7 m。多時相光學遙感影像表明該變形體已形成多條大型拉張裂縫,一旦發(fā)生崩落可能會堵塞雅魯藏布江形成堰塞湖,危及雅魯藏布江上下游。

圖10 探測大冰川運動位置分布Fig.10 Location distribution of detected large glacier movements

圖11 培龍貢支大冰川2019年8月29日至2019年9月29日方位向與視線向二維地表形變Fig.11 Two dimensional azimuthal and LOS deformation of the Peilonggongzhi Glacier from August 29, 2019 to September 29, 2019

圖12 果登冰川2020年1月至2020年4月南北向與東西向形變Fig.12 North-south and east-west deformation of the Guodeng Glacier from January 2020 to April 2020

圖13 澤巴隆巴溝巖崩危險隱患點2016年6月15日至2018年3月7日方位向與視線向二維形變Fig.13 Two-dimensional azimuthal and LOS deformation from June 15, 2016 to March 7, 2018 at the rock avalanche potential sites in the Zebalongba gorge

在探測到的111處崩滑災(zāi)害中,有17處滑坡靠近雅魯藏布江分布,若坡體發(fā)生崩落,其巖土體會直接沖入雅魯藏布江,在一定程度上可能會威脅沿線居民群眾生命財產(chǎn)安全及交通設(shè)施安全。17個滑坡的詳細信息如表2所示,部分滑坡光學遙感影像見圖14,可以看到滑坡形態(tài)特征已非常明顯,大多數(shù)滑坡體表面存在崩塌現(xiàn)象以及曾經(jīng)發(fā)生滑動的痕跡。在這17個滑坡中,最危險的滑坡體為達波滑坡。達波滑坡體2016年6月15日至2018年3月7日雷達方位向與視線向二維形變見圖15。該滑坡體是一巨型古滑坡體,長度約1054 m,寬度約1089 m,其在2017年11月之前處于穩(wěn)定狀態(tài),未觀測到任何形變信息。2017年11月18日米林縣6.9級地震誘發(fā)了該滑坡體的復(fù)活,在視線向產(chǎn)生約-10 m及方位向產(chǎn)生約3 m的變形,2020年該滑坡又勢于穩(wěn)定狀態(tài)。SAR強度圖上看到該滑坡體后緣已完全脫離,左右兩側(cè)裂縫已完全貫通,在極端事件下(地震)一旦發(fā)生失穩(wěn)會完全堵塞雅魯藏布江,危及上下游居民生命財產(chǎn)安全,該滑坡體應(yīng)受到高度重視,布設(shè)地面監(jiān)測設(shè)備開展連續(xù)跟蹤監(jiān)測。

圖14 靠近雅魯藏布江部分典型滑坡光學遙感影像Fig.14 Optical remote sensing images of typical landslides near the Yarlung Zangbo RiverRed curves indicate the landslide boundary

圖15 達波滑坡2016年6月15日至2018年3月7日方位向與視線向二維形變Fig.15 Two-dimensional azimuthal and LOS deformation of the Dapo landslide from June 15, 2016 to March 7, 2018

表2 17處滑坡詳細信息Table 2 Detailed information of 17 landslides

4 結(jié)論

文章通過284景Sentinel-1影像與111景ALOS/PALSAR-2影像聯(lián)合Stacking-InSAR、SBASInSAR以及SAR偏移量技術(shù)識別了雅魯藏布江下游米林-墨脫段大小不等的地質(zhì)災(zāi)害形變區(qū)共260處。災(zāi)害編目圖顯示探測到的大多數(shù)崩滑以及冰川、冰崩變形區(qū)分布在距離雅魯藏布江較遠、海拔較高的溝道及山峰。

在探測到的13處大型冰川中,澤巴隆巴冰川溝中的危險巖崩隱患點在不到兩年的時間里方位向累積形變已經(jīng)達到-13 m,并且該變形體已經(jīng)形成多條大型拉張裂縫,一旦發(fā)生崩落可能會形成堰塞湖,給雅魯藏布江上下游環(huán)境及居民帶來危害。

探測到的17處沿岸滑坡中,達波古滑坡受米林地震影響而復(fù)活,視線向形變量達-10 m,并且該滑坡體后緣已經(jīng)完全脫離,左右兩側(cè)裂縫完全貫通。該滑坡一旦失穩(wěn)會完全堵塞雅魯藏布江,危及上下游居民的生命財產(chǎn)安全,應(yīng)對該滑坡體提高重視,布設(shè)地面監(jiān)測設(shè)備來預(yù)防災(zāi)害的發(fā)生。

致謝:本文Sentinel-1數(shù)據(jù)由歐空局免費提供,30m分辨率SRTM DEM數(shù)據(jù)由美國USGS免費提供,再次一并致謝。感謝各位評審專家、編輯對文章提出的寶貴意見。