黄河流域刘家峡—兰州段滑坡灾害的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软件绘制,在此一并表示感谢！