基于SBAS-InSAR的西宁市滑坡识别和形变监测分析

何 丽,吴先谭,陈 欣,罗 芳,何政伟,薛东剑,白文倩,康桂川,张雨祥

(1.青海省青藏高原北部地质过程与矿产资源重点实验室,西宁 810000; 2.成都理工大学 a.地理与规划学院,b.成都理工大学 地球科学学院,成都 610059)

0 引言

中国黄土高原的面积约6.4×105km2,黄土厚度约为100 m～300 m[1]。黄土具有孔隙大、质地疏松、垂直节理发育等独特的构造特征,在降雨的条件下极易发生滑坡、崩塌、泥石流等地质灾害[2]。黄土滑坡是黄土高原分布最广的一类地质灾害[3]。近年来,由于黄土地区日益脆弱的地质环境,在强降雨条件下黄土斜坡极易发生倾倒、滑移,对人类生命财产威胁极大[4]。2010年榆林市发生大规模黄土崩塌,造成了严重的人员伤亡和财产损失[5](如2013年天水市受强降雨影响,诱发大规模黄土滑坡[6])。因此,大范围高效、精准地识别与监测具有潜在危险性的蠕变型黄土滑坡,对于滑坡防治、灾害预警以及保护人民生命财产安全具有重要意义。

随着遥感技术的发展,光学遥感[7]、无人机(UAV)摄影测量[8]、机载激光雷达(LiDAR)[9]和合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR[10])等技术逐渐成为滑坡识别和形变监测的主要手段。其中InSAR技术因具备高精度、高分辨率、覆盖范围大、综合成本低、全天时和连续跟踪微小形变等技术优势,使其在大范围滑坡识别和形变监测应用中得到青睐[11-13]。1996年Fruneau[14]首次证明了合成孔径雷达差分干涉测量技术(D-InSAR)在滑坡形变监测的有效性,同时也发现时空失相干和大气延迟对形变监测结果影响较大。随后,国际学者在D-InSAR技术的基础上提出了基于多景SAR影像数据的时间序列InSAR技术(TS-InSAR),包括永久散射体干涉测量(Persistent O Scatter InSAR,PS-InSAR)[15]、小基线集干涉测量(Small Baseline InSAR,SBAS-InSAR)[16],分布式散射体干涉测量SqueeSAR[17]等技术,这些技术能够获取精度更高、可靠性更强的滑坡形变结果,在大面积滑坡识别和形变监测中已得到广泛应用。目前,以SBAS为代表的TS-InSAR已被用于黄土滑坡的识别和形变监测[18-19]。如刘陈伟等[20]利用SBAS技术获取黑方台地区地面形变速率,并结合黄土滑坡的地形地貌特征,识别出黑方台地区6块存在较大形变的区域;Meng等[21]利用无人机(UAV)摄影测量与SBAS-InSAR技术表征甘肃省红壑岘黄土滑坡形变过程,对比分析发现SBAS-InSAR技术在蠕变变形监测方面更具有优势。

笔者以升降轨Sentinel-1A雷达影像作为数据源,在对时序InSAR技术分析的基础上,采用SBAS-InSAR技术,提取了西宁市的分布式散射体(Distributed Scatterers,DS),反演了研究区的形变速率和时间序列形变量,并结合Google Earth遥感影像圈定了活动滑坡边界,分析了西宁市活动滑坡的分布规律,以及滑坡累计形变与降雨的关系。本研究成果为同类型区域黄土滑坡识别和形变监测提供了借鉴。

1 研究区概况与数据来源

1.1 研究区概况

西宁市位于青藏高原东北部,隶属于湟水中游河谷盆地,地理位置:101°33′E～101°56'E;36°25'N～36°47N,区域面积为490 km2。研究区地势西北高东南低,地形切割强烈,沟壑纵横,区内发育一条近东西向的湟水河(图1)。研究区地层由老到新为长城系(Pta)、白垩系(K)、古近系(E)、新近系(N)、中更新统(Q2)、上更新统(Q3)、全新统(Q4)。岩土类型主要有黄土、泥岩砂岩和石膏岩互层、松散碎石土,均为易失稳变形地层。西宁市属高原半干旱大陆性气候,西宁地区全年降雨分配不均匀,年降雨量477.5 mm,主要集中在6月-9月,占全年总降雨量的67.4%,年降雨周期性变化明显。研究区地形主要以侵蚀构造低山丘陵与侵蚀堆积河谷平原为主,在低山丘陵前缘多发育高陡斜坡,加之区域降水集中,极易诱发滑坡、崩塌和泥石流等地质灾害。

图1 研究区位置

1.2 数据来源与数据处理

Sentinel-1A是由欧洲航天局(ESA)发射的C波段合成孔径雷达对地观测卫星,该卫星具有全天时、多种极化方式、时空分辨率高等特征。本研究采用Sentinel-1A卫星以干涉宽幅(IW)模式下VV极化方式获取的单视复数产品(SLC)进行时间序列形变反演,影像采集了2019年1月至2021年11月83幅升轨影像和87幅降轨影像。升轨影像与降轨影像在不同的入射角下,存在不同程度的几何畸变,形成不同的观测效果。结合升轨与降轨数据进行滑坡识别,能更有效地扩大监测范围。SAR数据的基本参数如表1所示,具体的技术流程如图2所示。

表1 Sentinel-1A卫星基础数据

图2 技术流程图

此外,笔者采用30 m分辨率的SRTM(shuttle radar topography mission) DEM(digital elevation model)去除干涉过程中的地形相位。选用POD Precise Orbit Ephemerides(POD精密定轨星历数据)(https://qc.sentinel1.eo.esa.int/)校正轨道信息,选用Generic Atmospheric Correction Online Service for InSAR(GACOS)(http://www.gacos.net/)大气延迟产品数据,该数据时间和空间分辨率较高,全球近实时,可有效去除干涉图中的大气延迟相位收集了西宁市气象站2019年-2021年日平均降水量数据,用于分析滑坡形变与降雨量之间的关系。

2 实验方法

Berardino等[16]提出的SBAS-InSAR技术是对分布式散射体(Distributed scatterer DS)进行相位分析获得时序形变。该技术对基于低分辨率、大尺度、时序形变监测能够取得很好的效果[22]。与传统的D-InSAR技术相比较而言,SBAS-InSAR技术能够较好地克服时空失相关和大气延迟效应等缺陷,实现更高精度的地表形变监测[23]。SBAS-InSAR技术实现过程为:

1)对原始的Sentinel-1A数据进行轨道信息校正以及裁剪到覆盖整个研究区范围。

2)综合考虑时间-空间基线,分别选取2019年3月10日和2019年11月10日的影像作为升轨和降轨数据的配准主影像。设置60 d时间基线,150 m空间基线,分别筛选出320个和304个干涉对(图3)。对干涉对进行多视、干涉图生成、去平地相位、Goldstein滤波、最小费用流(MCF)相位解缠等处理获得每个干涉图中任意像元(x,y)的差分干涉解缠相位。

图3 SBAS-InSAR时空基线图

3)形变速率和累计形变估算,形变分量可以进一步分解为线性形变分量和非线性形变分量。通过相位与时间基线间的最小二乘拟合,可以估算出线性位移。使用奇异值分解法(SVD)求出非线性形变最小范数意义上的解。

3 实验结果

3.1 SBAS-InSAR识别活动滑坡

通过SBAS-InSAR技术获取西宁市Sentinel-1A升轨和降轨数据雷达视线(Line of Sight,LOS)方向上的形变速率(图4、图5)。图4、图5中红色的点(负值)表示目标沿着LOS方向远离卫星移动,绿色的点表示在监测时间段内目标相对稳定,蓝色点(正值)表示目标沿着LOS方向靠近卫星移动。研究区升轨形变速率位于-98 mm/年～44 mm/年之间(图4(a)),降轨形变速率位于-95 mm/年～25 mm/年之间(图5(a))。由图4(b)、图5(b)可知,升轨和降轨形变速率均主要集中在-10 mm/年～10 mm/年的范围内,平均形变速率分别为-2.09 mm/年、-0.64 mm/年,且分布特征与正态分布相似,这表明研究区整体上处于稳定状态。基于形变速率信息,笔者将形变速率超过-15 mm/年区域确定为显著形变区。研究区共探测到74处显著形变区(升轨数据53处,降轨数据67处,升轨和降轨数据联合46处),总面积约9 km2。

图4 Sentinel-1A升轨形变速率图

图5 Sentinel-1A降轨形变速率图

为了进一步判断形变区是否为滑坡及圈定滑坡边界范围,结合显著形变区的分布情况,利用多期次Google Earth光学遥感影像和高分2号影像数据(图6),可以清晰地观察到地表物体形态。结合滑坡在光学遥感色调、纹理、平面形态、微地貌形态以及变形特征等解译标志,特别是在高分辨率影像中表现出明显的裂缝及滑塌现象,对SBAS-InSAR技术探测到的显著形变区进一步划分,最终划定了升轨滑坡、降轨滑坡以及其他形变区的分布范围(图7)。笔者将35处具有明显滑坡光学影像解译特征的显著形变区识别为滑坡,并在Google Earth影像中勾画出滑坡边界,39处具有明显人类活动特征的显著形变区识别为其他形变区(图7(a))。

图6 西宁市不同时期遥感影像图

图7 识别区结果分布及滑坡主要发育区局部放大图

经过分析,InSAR技术识别的35处滑坡中,升轨数据识别出14处、降轨数据识别出28处、升轨和降轨数据联合识别出7处(图7(a))。识别出的滑坡总面积约2.4 km2,其中面积最小为5 945.64 m2,最大为0.56 km2,规模以中大型为主。滑坡呈条带状集中分布于湟水河及其一级支流河谷平原区向低山丘陵区过渡的斜坡地带,且湟水河北岸滑坡发育密度大于南岸。图7(b)、图7(c)为研究区滑坡典型发育区,分别记为A区和B区。

图8(a)、图8(b)分别表示A和B典型区InSAR降轨数据形变速率空间三维分布,分析可知,A区和B区域分别发育4处滑坡(H1、H2、H3、H4)和6处滑坡(H5、H6、H7、H8、H9、H10),其基本特征如表2所示。

表2 重点区InSAR识别滑坡的基本特征

图8 研究区典型滑坡形变速率及空间分布

由表2及图8(c)、图8(d)可知,A区域中H2滑坡最大形变速率达48.87 mm/年,坡体形变集中在中部和前缘,在滑坡体前缘和中部发育3处次级滑坡。后缘下错约10 m并且发育多条走向近南北裂缝,长约180 m。中部有一系列横向延伸长度70 m～90 m的不连续圆弧形拉裂缝,且呈逐步贯通趋势,在自重作用下,上部岩土不断向下挤压。由表2及图8(e)、图8(f)可知,B区域中H8滑坡最大形变速率达81.26 mm/年,坡体形变集中在中部区域,滑体后部发育有三条裂缝,均呈东西向展布,平面形状呈弧形,最长一条达到170 m,宽0.2 m～0.6 m,滑坡体北侧分布呈带状陡崖,坡度近90°。现陡崖上发育多处危岩,危岩与母岩间裂隙宽10 cm～30 cm,上下贯通,呈孤立状(图8(d))。

识别到的39处其他形变区在影像可清晰识别为挖填区(治沟造地、土地整平、垃圾填埋、流域治理)、建筑施工区、道路沉降区、矿区和工业区等。经过现场调查核实,最终确定13处为挖填区、16处为建筑施工区、2处为道路沉降区、4处为矿区和4处为工业区。图9(a)～图9(c)分别显示了建筑施工、工业区和治沟造地引起的地表形变。

图9 研究区其他形变区展示图

4 讨论

4.1 滑坡时间序列形变与降雨关系

降雨是诱发黄土滑坡的重要因素之一,这与黄土孔隙比、湿陷性等独特性质密切相关。本研究以张家湾滑坡和付家寨滑坡为例,分析降雨与滑坡时间序列累计形变(SBAS-InSAR)的演化关系。

4.1.1 张家湾滑坡形变与降雨关系

该滑坡地理坐标为101°39′04″E,36°38′25″N。滑坡高差约325 m,坡度30°,面积5.6×106m2,总体积2.2×106m3,表层披覆黄土,属平缓层状岩质斜坡,为特大型推移式滑坡。滑坡呈不规则状,坡体有多处浅表层滑塌迹象,后缘裂缝明显(图10(a))。图10(b)为该滑坡Sentinel-1A升轨影像LOS方向形变速率图,形变速率范围为-25 mm/年～10 mm/年,形变部位集中在滑坡后缘及中部具有明显拉裂缝的区域。根据形变点分布规律,可将滑坡分为四个显著性形变区(Ⅰ、Ⅱ、Ⅲ、Ⅳ)。经野外调查可得,滑坡Ⅰ区后缘贯通,有土质垮塌的现象,前缘的挡墙有开裂(图10(c))。Ⅱ区因人类工程活动导致部分地带有失稳的可能。Ⅲ区发育走向130°的张性拉裂缝,下错明显约80 cm,形成新的不稳定体,目前该区已进行了滑坡治理。Ⅳ区后缘可见局部的垮塌,且发育多条裂缝(图10(d))。滑坡体上季节性降雨导致裂缝、落水洞、座落等现象时有发生。分别在四个区域取一点(P1、P2、P3、P4)进行时间序列累计形变量与日平均降雨量关系分析(图10(e)),四个区域的累计形变量均呈现波动增加,其中P1点累计形变量最大为-78.6 mm,形变曲线反映该区域处于快速形变阶段。滑坡在2019年6月-2020年2月和2020年8月-2021年2月(浅蓝色区域),由于雨季和冻融沉降的到来,形变量呈快速增加趋势,可见滑坡形变与日平均降雨量密切相关,具有明显的季节性特征,尤其是在雨季变形明显加速,并且变形加速时间与雨季开始时间存在明显的滞后效应。

图10 张家湾滑坡InSAR识别结果

4.1.2 付家寨滑坡形变与降雨关系

该滑坡于2019年11月18日发生在城东区付家寨北部,地理坐标:E:101°54′10.54″、N:36°34′52.03″。滑坡长270 m,横310 m,滑向295°,平均坡度25°,平面形状近似正方形,滑坡体较为破碎,发育多个次级滑坡,多条裂缝贯穿,后缘下错明显(图11(a)～图11(b))。

图11 付家寨滑坡InSAR探测结果

图11(c)为该滑坡Sentinel-1A降轨影像LOS方向形变速率图,坡体形变集中在中前部区域,最大形变速率达40 mm/年,结合地形条件和形变空间分布特征可将滑坡分为两处显著形变区(Ⅰ、Ⅱ)。在两处形变显著区域各选1个特征点P1、P2,其中P1平均形变速率为31 mm/年,累计形变量为92.7 mm。该形变区前缘堆积体堆覆沟道内,挤压沟道,堆积高度达3 m,灾害造成4间房屋损毁(图11(d))。P2平均形变速率为38 mm/年,累计形变量为105.6 mm。该形变区为一中型牵引式次级滑坡,滑体中部树木歪斜,绿化管网断裂,道路垮塌,后缘下错约5 m,前缘堆积体堆覆沟道内,挤压沟道,堆积高度达3 m(图11(e))。

图11(f)为P1特征点、P2特征点累计形变量与日平均降雨量关系图。可知累计形变曲线可分为两个阶段,2019年1月-2020年10月时间段形变量快速增加,P1点、P2点累计形变量分别增加81 mm、89 mm。2020年10月之后形变趋势变缓慢,整体趋于稳定。形变曲线在2019年5月和2020年5月雨季之后有明显的加速变形趋势(图中浅蓝色区域),说明该滑坡变形与日降雨量密切相关。

5 结论

研究区位于黄土地区,地表黄土物质广泛覆盖,植被稀少,区域内相干性良好,采用升降轨Sentinel-1A数据,基于时序InSAR技术对研究区活动滑坡进行识别和形变监测。结论如下:

1)通过对原始Sentinel-1A数据进行裁剪,综合考虑时间和空间基线,筛选出若干组干涉对,随后进行干涉处理(影像配准、生成干涉图、去地形相位、相位解缠等)以及轨道精炼和重去平,最后对形变速率进行估算,识别出区域的显著形变区域。同时结合高分辨率卫星影像,进一步判定及划定滑坡边界范围。共探测到74处显著形变区,经过分析,35处为活动滑坡,其中升轨识别出14处、降轨识别出28处、升轨和降轨联合识别出7处。升轨和降轨数据相结合,能够有效地减少因几何畸变导致的监测“盲区”,从而提高卫星数据的有效探测率。

2)以张家湾滑坡和付家寨滑坡两个典型滑坡为例,分析了降雨与滑坡时间序列累计形变与降水的关系。分析发现滑坡受降雨入渗的影响,造成黄土、部分泥岩、石膏块软化与崩解,从而使土体物理力学性质发生改变,相对集中且强烈的降雨最终导致滑坡稳定性降低,内在(工程岩组、断层)和外在(降雨)条件的共同作用导致已识别的滑坡在研究时段处于活动状态。