黃河流域劉家峽—蘭州段滑坡災害的InSAR識別及成因分析

吳東霖,葛偉鵬,2,魏聰敏,張波,2

(1.中國地震局蘭州地震研究所,甘肅 蘭州 730000;2.蘭州地球物理國家野外科學觀測研究站,甘肅 蘭州 730000)

0 引言

中國地質(zhì)災害分布廣泛且種類繁多,其中滑坡所占比例超過70%[1]。滑坡災害在我國分布廣泛、發(fā)生頻率高,每年都會造成巨大的人員傷亡、財產(chǎn)損失和環(huán)境破壞,尤其是在我國西南、西北地區(qū)。如2017年6月24日四川茂縣新磨村發(fā)生山體高位滑坡,造成83人失聯(lián)或死亡,64戶農(nóng)房被掩埋,岷江支流堵塞近2 km[2]。因此滑坡災害的形變監(jiān)測和隱患識別意義重大。

合成孔徑雷達干涉測量技術(shù)(InSAR)由于其全天候、全天時、范圍廣、高精度、低成本等特點,近年來逐漸被國內(nèi)外學者投入滑坡研究中并取得了較好的效果。1996年Achache等[3]將DInSAR技術(shù)引入滑坡領(lǐng)域,觀測研究Saint-Etienne-de-Tinee滑坡的形變,得到的結(jié)果與地表測量結(jié)果具有高度一致性,證明InSAR技術(shù)可用于滑坡形變監(jiān)測;2000年Vietmeier等[4]通過觀測La Valette滑坡的形變速率,進一步證明了InSAR技術(shù)在滑坡變形監(jiān)測中的可行性和準確性;2001年Ferretti等[5]利用PS-InSAR技術(shù)監(jiān)測Ancona地區(qū)滑坡,證實了PS點形變監(jiān)測精度可達1 mm;2004年Hilley等[6]利用PS-InSAR技術(shù)監(jiān)測美國Hayward斷層附近區(qū)域,并指出降水可能引發(fā)滑坡;2007年程滔[7]在研究區(qū)布設(shè)角反射器(CR),結(jié)合PS-InSAR對陜西子長縣滑坡進行了監(jiān)測研究;2010年王騰[8]提出一種準永久散射技術(shù)(QPS-InSAR),并以此監(jiān)測了巴東新城區(qū)的兩個滑坡;2016年康亞[9]采用InSAR Stacking、IPTA、SBAS-InSAR和PS-InSAR技術(shù)對西南山區(qū)進行形變監(jiān)測,成功探測出多處滑坡災害點;2020年張毅等[10]綜合利用InSAR和無人機測繪技術(shù)提出一種預測潛在滑坡位置及范圍的方法,并得到了很好的驗證。

基于上述研究,本文提取2014年10月—2019年12月覆蓋黃河流域劉家峽—蘭州段的升軌Sentinel-1數(shù)據(jù),利用PS-InSAR技術(shù)結(jié)合GPS形變率對數(shù)據(jù)進行改正,得到研究區(qū)LOS向形變速率;將研究結(jié)果與GPS結(jié)果及前人研究進行對比,驗證研究結(jié)果的有效性,并進一步分析形變典型區(qū)的形變特征和潛在危險性。

1 研究區(qū)與數(shù)據(jù)概況

1.1 研究區(qū)概況

本文研究區(qū)是黃河流域劉家峽—蘭州段(圖1),主要位于甘肅省臨夏州北部永靖縣,西接青海省民和縣,南與積石山縣、東鄉(xiāng)縣隔黃河相望,東北與蘭州市接壤,東南與定西市接壤。該區(qū)域鄰近青藏高原和黃土高原過渡帶,屬隴西中新生界盆地構(gòu)造,是地震多發(fā)地帶。境內(nèi)平均海拔約2 000 m,地形復雜,山川交錯,河谷縱橫,屬盆地邊沿的次高山群及高原淺山丘陵區(qū),大部分地面被黃土覆蓋,地勢東西高、中部低[11]。區(qū)內(nèi)水資源較為豐富,黃河干流從積石山縣由西南進入劉家峽水庫后轉(zhuǎn)向東北,在劉家峽鎮(zhèn)與支流洮河交匯,后又轉(zhuǎn)向西北,經(jīng)太極鎮(zhèn)、小川—大川盆地后再轉(zhuǎn)向東北,經(jīng)鹽鍋峽鎮(zhèn)與支流湟水匯合,后一路向東進入蘭州市區(qū)。研究區(qū)屬大陸型半干旱氣候,陽光充足,年日照2 200～2 800 h,氣候溫和,年平均氣溫8～9 ℃,年平均降水280 mm,60%～70%的降水集中在5—8月[11]。

圖1 研究區(qū)活動斷裂與地震分布圖(圖中紅緣黑色箭頭表示本文所用7個GPS站點的水平速度)Fig.1 Distribution of active faults and earthquakes in the study area

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

本文選取歐空局2014年發(fā)射的Sentinel-1A衛(wèi)星的干涉寬幅(interferometric wide-swath,IW)模式、Level-1 SLC (single look complex)級別、升軌的C波段SAR數(shù)據(jù)為數(shù)據(jù)源,其空間分辨率為5 m×20 m,時間分辨率為12 d。所用SAR數(shù)據(jù)共111景,圖1中紫色矩形框表示其覆蓋范圍(軌道號為128),監(jiān)測時段選取2014年10月—2019年12月。同時選取歐空局提供的精密軌道星歷數(shù)據(jù)(POD Precise Orbit Ephemerides)和精密軌道輔助文件(Precise Orbit Auxiliary File)。DEM采用由美國宇航局(NASA)提供的SRTM數(shù)據(jù),其分辨率為30 m×30 m。

隨后引入研究區(qū)GPS速度場來校正InSAR初始結(jié)果,使用的GPS觀測數(shù)據(jù)來源于中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(CMONOC)二期工程以及國際GNSS服務提供的部分IGS觀測數(shù)據(jù)、星歷數(shù)據(jù)和表文件數(shù)據(jù)。利用GAMIT[12]進行單日松弛解處理,再利用GLOBK[13]軟件將其與SOPAC (Scripps Orbital and Position Analysis Center)提供的ITRF2014全球框架解合并,并通過卡爾曼濾波得到站點在ITRF2014參考框架下的位置和速度,最后再將其轉(zhuǎn)化到鄂爾多斯塊體區(qū)域參考框架下。

2 研究方法

2.1 PS-InSAR基本原理

在SAR影像中,每個像素的回波是其相關(guān)地面分辨單元中所有散射體的相干和,而數(shù)據(jù)獲取過程中散射體間相對運動和雷達視角的改變會使相應散射體返回能量以不同方式求和,從而引發(fā)影像去相關(guān)[14]。當分辨單元內(nèi)某一散射體返回能量遠大于其他散射體時,對應像素主要受這一個散射體影響,則去相關(guān)程度會大大減小,這種散射體即被定義為永久散射體(Persistent scatterer/Permanent scatterer,PS)[14]。人居環(huán)境下常見的PS點有建筑物屋頂、道路等人造地物,人跡稀少環(huán)境下則多為特定方向的裸露巖石。這些目標的散射特性較為穩(wěn)定,雷達回波信號反射能力很強,具有較高信噪比,幾乎不受時間和空間基線影響,能在較長時間內(nèi)保持較高相干性[15]。

PS-InSAR技術(shù)本質(zhì)上仍然屬于D-InSAR,其基本原理是:假設(shè)研究區(qū)內(nèi)有N幅SAR影像,基于空間基線最小、時間基線最小、多普勒質(zhì)心頻率基線最小等原則選定其中一幅為主影像,其余為輔影像[16],配準、干涉后可形成N-1幅差分干涉圖,再利用振幅離差閾值、相位離差閾值、小波相位分析等方法從差分干涉序列SAR影像中選取PS點,對這些PS點進行時間序列分析得到其形變量,進而得到整個研究區(qū)的地表形變場。

2.2 趨勢改正原理

PS-InSAR得到的初始視距向形變速率存在一個斜坡趨勢[17],引入測區(qū)內(nèi)GPS站點速率對其進行改正,同時也可以給InSAR結(jié)果一個參考框架。

首先,計算每個GPS站周圍圓形(半徑約400 m)范圍內(nèi)InSAR的反距離加權(quán)LOS向速率和入射角。其次,通過式(1)將GPS站速率(VN、VE、VU)轉(zhuǎn)化合成到雷達視線向(Vlos)[18],并求出其與對應InSAR視距向加權(quán)速率的差值:

Vlos=VNsinθsinα-VEsinθcosα+VUcosθ

(1)

式中:θ和α分別代表雷達入射角和衛(wèi)星飛行方位角,對于Sentinel-1升軌,取α=-15°[18],θ即取上步算得加權(quán)入射角。

接著,使用最小二乘法令式(2)中的R2最小以求出趨勢參數(shù)m1、m2、m3,變形后的m1、m2、m3則由式(3)求出:

(2)

(3)

式中:loni、lati和ΔVi分別代表第i個GPS站點的經(jīng)度、緯度及其GPS與InSAR視距向速率差值。

最后,移除InSAR視距形變率初始值的線性趨勢。

3 數(shù)據(jù)處理與結(jié)果分析

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

根據(jù)時間基線和垂直基線最優(yōu)原則選取2017年10月22日的SAR影像為主影像,其余為輔影像,組成110組干涉對。干涉對信息如圖2所示,干涉對空間基線最大為139 m,符合要求。

圖2 時空基線分布Fig.2 Distribution of temporal-spatial baselines

基于ISCE平臺[19]進行預處理,再基于StaMPS平臺[5]采用StaMPS方法來提取PS點,最終得到2 884 074個PS點(振幅離差閾值為0.4),其在城鎮(zhèn)地區(qū)分布較密集,野外山區(qū)分布較稀疏(圖3)。基于選定的PS點進行形變反演,依次去除殘余地形和大氣延遲等影響,解算出各點的形變速率[圖3(a)]。以區(qū)域內(nèi)7個GPS點三維速率為約束,去除InSAR初始結(jié)果中的線性趨勢后,得到修正的鄂爾多斯塊體區(qū)域參考框架下的PS點LOS向形變速率[圖3(b)]。

圖3 InSAR LOS向形變速率改正前后對比圖(圖中暗紅色橢圓為袁道陽等[20]推測的1125年蘭州M7.0地震極震區(qū))Fig.3 InSAR LOS deformation rate before and after calibration

3.2 結(jié)果分析

將改正后的InSAR結(jié)果裁剪出研究區(qū)范圍,得到圖4(a)。由圖4(a)可知,相對于鄂爾多斯地塊,整個區(qū)域基本處于LOS沉降狀態(tài),沉降嚴重區(qū)主要分布在黃河流域及其支流和公路附近,沉降最顯著的區(qū)域位于三條峴鄉(xiāng)下莊村附近(圖4中橢圓區(qū)域),而非研究熱點黑方臺。將沉降最顯著區(qū)域——下莊沉降區(qū)放大,得到圖(4b)。由圖(4b)可知,下莊沉降區(qū)表現(xiàn)為橢圓形,其南側(cè)是洮河,西、北、東側(cè)均為居民聚居地,如紅峴子、大臺子、下莊、青和等村落,LOS沉降速率>4 mm/a,面積約25 km2,體量遠大于黑方臺;沉降最大值點距下莊、青和等村約2 km,距洮河約4 km,距永靖縣城約12 km,距劉家峽水庫約20 km,距黑方臺約20 km,距蘭州市區(qū)約40 km。

由圖4(b)可見,沉降速率>10 mm/a區(qū)域也表現(xiàn)為橢圓狀。為了進一步分析下莊沉降區(qū)的變形特征,我們沿該橢圓的長軸(BB′)和短軸(AA′)各作一條剖面(圖5)。如圖5所示,灰色圓點為LOS向形變速率剖面,綠色實線為對應的地表高程剖面。根據(jù)形變速率剖面可知,沉降速率>8 mm/a的長軸約4 km,形成一個較大的沉降漏斗[圖5(b)]。中心沉降點P1的LOS向形變速率為-14.4 mm/a,約位于橢圓右焦點。結(jié)合高程剖面,山頂點P2和P3的沉降速率約為0 mm/a和-1 mm/a,分別向西和向南延伸至山下達到最大沉降值,說明若滑坡被觸發(fā),可能會向南方滑動,對東南的居民聚集地造成巨大危害,甚至可能堵塞南側(cè)的洮河。

圖4 研究區(qū)LOS向形變速率Fig.4 LOS deformation rate of the study area

圖5 沉降剖面Fig.5 Subsidence profile

圖3和圖4(a)中疊加了袁道陽等[20]推斷的1125年蘭州M7.0地震極震區(qū),可看出下莊沉降區(qū)處于極震區(qū)內(nèi),因此推斷該滑坡隱患點可能與1125年蘭州地震有關(guān)。結(jié)合圖6(a)和圖6(b)中山底點P1與山頂點P2、山頂點P3間的時間序列差值可看出,2017年7月30日—8月23日沉降差值分別累計近20 mm和30 mm,之后未回到原位置,產(chǎn)生瞬態(tài)滑移,引起永久變形。2018年8月雖然也存在這樣的趨勢,但差值并非持續(xù)增大,后續(xù)基本回到原位置,故推測每年8月沉降差值的增大或與降水有關(guān)。2017年8月的差值突增則可能是受到降水和2017年8月8日九寨溝MS7.0地震的雙重影響。

圖6 時間序列Fig.6 Time series of some special points

通過實地現(xiàn)場調(diào)查,沉降最大點P1處未見到明顯滑坡形態(tài),僅有一些邊坡失穩(wěn)現(xiàn)象。通過對沉降區(qū)東南方下莊村居民進行走訪,得知不少居民家房屋有裂縫產(chǎn)生,最為明顯的是P4點某居民家墻上產(chǎn)生了約1 cm寬的裂縫(圖7)。房屋建設(shè)已近十年,裂縫約產(chǎn)生于2017年,與時間序列顯示的形變特征相符[圖6(c)],證明了InSAR結(jié)果的可信度。

圖7 P4點下莊村某居民家墻上的裂縫Fig.7 A crack on the wall of one resident’s house in Xiazhuang village (Point P4)

3.3 精度評定

為驗證InSAR結(jié)果的精度,將GPS臺站的N、E、U分量合成LOS向與InSAR結(jié)果進行對比,將結(jié)果列于表1。由表1可知其較差基本較小,證明InSAR結(jié)果精度較好。其中G095和G360兩點GPS和InSAR結(jié)果的差值較大,可能是因為這兩點接近圖幅邊緣,受到的約束較差。

表1 GPS與InSAR視距向速率對比Table 1 Comparison between LOS velocity of GPS and InSAR

最后,對比本研究得到的黑方臺地區(qū)的LOS向速率和前人結(jié)果以驗證InSAR結(jié)果的可靠性。前人對黑方臺地區(qū)形變的研究較為豐富,其中王晨興等[21]得到黑方臺LOS向速率為-4～-13 mm/a,曾珠[22]得出LOS向區(qū)域均值為-2.7 mm/a,史緒國等[23]、朱文峰等[24]和張毅等[10]得到的LOS向結(jié)果最大值分別為-60 mm/a、-70 mm/a和-64 mm/a。經(jīng)ArcGIS統(tǒng)計,本研究得到的黑方臺地區(qū)視線向形變速率為-0.05～-10.09 mm/a,均值為-4.19 mm/a[圖4(a)],與王晨興和曾珠的結(jié)果較為相近,與朱文峰和張毅的結(jié)果相差較大。這可能是因為前兩者和本研究一樣采用PS-InSAR方法,而后三者采用SBAS-InSAR方法;在研究期間發(fā)生滑坡的區(qū)域會失相干,長時間序列PS-InSAR方法無法在這些區(qū)域提取到PS點,這可能導致較大滑動速率觀測值的丟失[25]。

4 結(jié)論

本文利用PS-InSAR技術(shù)處理111景覆蓋黃河流域劉家峽—蘭州段的Sentinel-1 C-SAR數(shù)據(jù),采用GPS數(shù)據(jù)對其進行趨勢改正后獲得研究區(qū)2014年10月—2019年12月間的年平均地表LOS向形變場;通過對比發(fā)現(xiàn)InSAR結(jié)果與GPS的較差小于2.3 mm/a,且與同期前人研究較為符合,驗證了研究結(jié)果的有效性。

研究區(qū)沉降主要分布在黃河及其支流沿線和交通干線附近,最大沉降區(qū)位于永靖縣三條峴鄉(xiāng)下莊村附近,其規(guī)模和量級都大于研究熱點黑方臺。如果被外力(如強降雨、周邊強震等)觸發(fā)滑坡,對周圍村落和附近河流、交通都將造成嚴重損失,應予以重視。通過歷史地震和時間序列分析,此沉降區(qū)可能與1125年蘭州M7.0地震有關(guān),2017年九寨溝M7.0地震可能對其產(chǎn)生了觸發(fā)作用。

本文通過長時間序列PS-InSAR監(jiān)測獲取了大范圍的滑坡監(jiān)測,但會丟失一些較大的快速形變,后續(xù)可以考慮和短時間序列PS-InSAR、SBAS-InSAR或機載/地基SAR相結(jié)合。本文的結(jié)果是LOS向形變場,只能較為粗糙地分析形變,后續(xù)考慮利用升降軌和GPS聯(lián)合解算出該區(qū)域的三維速率場。

致謝:感謝ASF DAAC(https://asf.alaska.edu/)為本文提供的歐空局哥白尼Sentinel-1數(shù)據(jù)[201410—201912];感謝國家重大科學工程“中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡”(http://www.neiscn.org)為本文提供的GPS觀測資料;感謝IGS數(shù)據(jù)中心(http://www.igs.org/)為本文提供的全球IGS站資料;感謝中國地震臺網(wǎng)中心國家地震科學數(shù)據(jù)中心(http://data.earthquake.cn)提供的數(shù)據(jù)支撐。本文采用ISCE、StaMPS、GAMIT/GLOBK、ArcGIS等軟件進行數(shù)據(jù)處理,部分圖件由GMT軟件繪制,在此一并表示感謝！