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

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

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

0 引言

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

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

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

1 研究區與數據概況

1.1 研究區概況

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

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

1.2 數據簡介

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

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

2 研究方法

2.1 PS-InSAR基本原理

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

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

2.2 趨勢改正原理

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

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

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

(1)

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

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

(2)

(3)

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

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

3 數據處理與結果分析

3.1 數據處理

根據時間基線和垂直基線最優原則選取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),其在城鎮地區分布較密集,野外山區分布較稀疏(圖3)。基于選定的PS點進行形變反演,依次去除殘余地形和大氣延遲等影響,解算出各點的形變速率[圖3(a)]。以區域內7個GPS點三維速率為約束,去除InSAR初始結果中的線性趨勢后,得到修正的鄂爾多斯塊體區域參考框架下的PS點LOS向形變速率[圖3(b)]。

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

3.2 結果分析

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

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

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

圖5 沉降剖面Fig.5 Subsidence profile

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

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

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

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

3.3 精度評定

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

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

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

4 結論

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

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

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

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