基于时序InSAR的青海同仁市滑坡隐患早期识别与特征解析

李宗仁, 沙永莲, 辛荣芳*, 张兴, 隋嘉, 孙皓

(1.青海省地质调查院, 西宁 810012; 2.青海省遥感大数据工程技术研究中心, 西宁 810012; 3.青藏高原北部地质过程与矿产资源重点实验室, 西宁 810012; 4.青海省地质环境监测总站, 西宁 810001)

青海省地质灾害主要发育在山间谷地,具有点多面广的特点[1],面对严峻复杂的防灾形势,省内先后开展了系统全面的地质灾害详查和多轮次拉网式隐患排查,每年汛期不断巡查,并在此基础上建立了较为完善的群测群防体系。但近年来发生的很多地质灾害事件中有部分并不在已知隐患范围内,对于隐蔽性强、变形缓慢的隐患,传统的人工排查手段已显无力。当前防范地质灾害的核心需求是要搞清楚“隐患点在哪里”“什么时候发生”,为此,如何突破传统人工调查的局限,快速、有效地识别地质灾害隐患是青海省防灾工作急需解决的问题[2-3]。

目前,青海省地灾隐患主要靠人工地面摸排进行识别,这仅能从点上发现明显的灾害征兆,如地裂缝、鼓包、下错以及建筑物变形等迹象,面对大范围的整体位移且无明显变形迹象时难以察觉。近年来虽然开展了重点区域的光学遥感调查,实现了地质灾害的影像识别与变化监测,但也只能发现已发生或正在发生的灾害体形态变化,对于缓慢变形的隐患识别能力仍有不足。合成孔径雷达干涉测量(synthetic aperture radar interferometry,InSAR)技术具有全天候、全天时、覆盖范围广、空间分辨率高、非接触、综合成本低等优点,适宜于开展大范围地质灾害普查与长期持续观测,特别是其大范围连续跟踪微小形变的特性,对正在变形区具有独特的识别能力[4],但是,青海省地质灾害监管部门还未自行开展过地质灾害InSAR早期识别业务,导致防灾工作比较被动。现通过长时序SAR数据对同仁市进行持续形变观测,配合孕灾背景资料及高分辨率光学影像,识别圈定滑坡灾害隐患,并对隐患点形变特征进行时序分析,判断变形阶段,评估其危险性与风险,旨在为青海省地质灾害防治工作提供技术支撑。

1 研究区概况

同仁市属青海省黄南藏族自治州管辖,位于青海省东南部、黄南州东北部,市域东西宽75 km,南北长85 km,总面积3 275 km2,地理坐标为东经101°38′～102°27′,北纬35°01′～35°47′。研究区地处青藏高原与黄土高原相接的过渡地带,境内山峦起伏,河谷相间,地势南高北低,高差悬殊,隆务河纵贯全境南北,形成了东-西部山区和中部河谷的两峡-谷地地形特征(图1)。该区属典型的高原温带半干旱气候区,具有寒冷干燥,日照充足,降水集中等特点,年内降水分配极不均匀,多集中在5—9月,且降水量随地势增高垂直分带性明显。区内已造就的现代地形其发展演变过程是由地球内、外力共同作用的结果——地形破碎、坡陡沟深、沟壑纵横、地表裸露,地质环境较为脆弱,引发的环境地质问题主要有滑坡、崩塌、泥石流等地质现象以及修建道路引起的边坡失稳等[5]。

图1 研究区地理位置及地貌特征图Fig.1 Map of geographical location and geomorphic characteristics of the study area

2 研究思路与方法

地质灾害隐患早期识别是地质灾害防治工作的重要组成部分,相较于传统地面调查与光学遥感解译,时序InSAR技术作为一种新型的对地观测技术,不仅可以识别连续、缓慢的地表形变信息,同时还可以提升灾害隐患的识别效率[6]。研究区由于其特殊的地理环境,地质灾害发育具有数量多、类型全、规模大的特点,现针对区内的滑坡及潜在滑坡,采用基于Sentinel-1A升降轨数据的差分干涉测量短基线集时序分析技术(small baselines subset InSAR, SBAS-InSAR)进行地质灾害隐患的早期识别,技术流程如图2所示。

(1)数据选择。以Sentinel-1A卫星SAR数据为主要数据源,选用2019年1月—2020年10月覆盖研究区的52景升轨影像和45景降轨影像进行区域形变信息提取(表1、图3)。除SAR数据外,还采用了ALOS卫星30 m分辨率的DEM数据用于消除干涉过程中的地形相位,降低地形误差,同时还采用了精密定轨星历数据(POD)作为轨道数据对轨道信息进行修正,消除因轨道误差产生的系统性误差;采用了2020年5—9月覆盖研究区的14景空间分辨率为1 m的国产高分二号(GF-2)光学影像对隐患点进行筛选与确认(图3)。高分二号卫星是中国自主研制的首颗空间分辨率优于1 m的民用光学遥感卫星,具有分辨率高、定位精度准、姿态机动快等特点,其亚米级空间分辨率能够有效识别地质灾害形态要素。

图2 本次研究技术流程图Fig.2 Technical flow chart of the study

表1 研究所用Sentinel-1A数据参数表Table 1 Main parameters of Sentinel-1A data

图3 本次研究所用卫星数据覆盖分布图Fig.3 Satellite data coverage map used in this study

(2)SAR数据处理。对覆盖研究区的总计97景升降轨数据进行SBAS-InSAR技术处理,包括连接图生成、差分干涉、相位解缠、轨道精炼与重去平、SBAS反演、地理编码等一系列工作流[7-8]。首先进行干涉像对的配对,对升、降轨数据分别设置最大时间基线为60 d和48 d,通过参数设置及经验分析,剔除质量差的干涉对,最终得到有效干涉对183组和105组,影像配对情况良好;其次进行滤波与干涉处理,对升降轨数据均设置多视比为10∶2,提高数据信噪比,降低方位向与距离向分辨率之间的差异。同时采用高斯滤波(Goldstein)方法进行滤波处理,提高干涉条纹的清晰度和减少时空基线引起的失相干,从而提升差分干涉的稳定与可靠性;利用狄洛尼最小费用流(Delaunay minimum cost flow, DMCF)方法进行相位解缠处理[9-10],实现大范围的全局最优解缠结果,处理过程中对相干性低或存在明显大气相位的像对进行了移除,升降轨数据分别移除干涉对18组和35组;通过选取地面控制点(GCP)实现干涉结果的轨道精炼与重去平,利用相位解缠图和相干性图辅助完成GCP点的选取,选择的控制点均匀分布在高相干区域;SBAS形变反演分两次进行,第一次反演是在忽略大气效应的基础上进行地表形变速率和高程改正数估算,同时通过二次解缠进行干涉图优化。第二次反演是使用定制滤波去除低频形变信息、高程误差以及大气误差,得到最终的时间序列位移结果;最后将SBAS反演结果转化成统一的地理坐标,实现SAR数据从斜距标坐标到地理坐标系的转换。

(3)滑坡隐患识别与特征解析。在SBAS-InSAR反演的地表形变信息基础上,以“青海省同仁县1∶5万地质灾害详细调查(2018)”成果资料中地质灾害分布情况为参考,根据GF-2影像中判识的地形地貌、植被覆盖等特征,圈定具有孕灾条件且有威胁对象斜坡处的形变集中区,根据光学影像特征确定滑坡隐患的边界范围,经与收集的“青海省地质灾害隐患数据库(2020年)”套合比对,完成“是否新增隐患”属性的判识。并选取典型滑坡隐患从地表形变的时空分布、光学影像特征、实地变形特征等方面进行特征分析。

3 基于升降轨数据的滑坡隐患识别

3.1 区域地表形变SBAS反演结果

基于SBAS-InSAR技术获得了同仁市在2019年1月—2020年10月的区域地表形变信息。从升轨数据反演结果中可以看出,监测周期内大部分区域处于稳定状态,其中远离卫星视线方向的最大累积形变量为130 mm,靠近卫星视线方向的最大累积形变量为56 mm,平均累积形变量为5 mm,强形变区域多分布于隆务河西岸的陡峭斜坡上[图4(a)]。从中识别出8个形变特征点A01～A08,经时序特征分析,以2019年1月9日为基准,各特征点均表现为随着时间的变化其累积形变量不断增加。从量级上看,监测时段内A01与A02的累积形变量最大,达105 mm。根据研究区降雨特征,在2020年汛期,各特征点形变量呈明显加速趋势(图5)。

基于降轨数据的反演结果显示,研究区在监测周期内远离卫星视线方向的最大累积形变量为95 mm,靠近卫星视线方向的最大累积形变量为73 mm,平均累积形变量为1 mm[图4(b)],共识别出11个形变特征点,全部特征点均呈持续变形趋势,不同特征点的累积形变量级与形变特征存在差异,其中,D08的累积形变量最大,达85 mm;D06、D07、D10、D11的累积形变量较小,且变化相对较缓。受降雨影响,雨季时(曲线图中橙色部分)各特征点均出现形变加速或波动现象,此外,D10和D11在2020年3月出现波动是因冰雪消融产生的影响(图6)。

图4 基于升降轨数据的累积形变量反演结果Fig.4 Inversion results of cumulative deformation based on ascending and descending data

图6 降轨数据中形变特征点时序曲线图Fig.6 Time series curve of deformation characteristic points in descending data

3.2 滑坡隐患早期识别

通过时序SBAS-InSAR技术反演的区域地表形变结果是进行滑坡隐患早期识别的重要依据[11],在形变速率、累计形变量分析基础上,以形变特征点的形变数量级、坡体是否具备孕灾条件、形变速率是否突变、是否具有威胁对象等综合因素为参考,结合GF-2光学影像进行筛查,剔除明显受地形或大气影响、不满足滑坡发育条件、长时间处于冰雪覆盖、农田耕作区等形变特征区域(点),在区内共识别出19处滑坡隐患,集中分布于隆务河两侧的斜坡上,其中升轨数据中8处,降轨数据中11处,主要涉及黄乃亥乡、年都乎乡、兰采乡、隆务镇、加吾乡、曲库乎乡、扎毛乡等7个乡镇,配合GF-2光学影像明确了各滑坡隐患坡体的位置及其边界范围,滑坡隐患早期识别详情见表2,各隐患分布情况如图7所示。

基于升轨数据识别出的8处隐患中A-02、A-03、A-04、A-05和A-08与历史地质灾害面重合,A-01、A-06和A-07为新识别隐患,降轨数据中D-02、D-05、D-09与已知灾害点重合,各隐患点形变特征详如图8所示。

3.3 识别结果有效性分析

3.3.1 形变速率精度分析

经形变速率平均精度计算公式统计,升轨数据中共获取有效观测值16 346 546个,其中有效形变点的年均形变速率集中于-10～10 mm/a,形变速率精度范围为0～8 mm,平均形变速率精度2 mm,标准差为4 mm;降轨数据中共获取有效观测值15 748 139个,其中有效形变点的年均形变速率集中于-20～20 mm/a,形变速率精度范围为0～11 mm,平均形变速率精度为0,标准差为5 mm(图9)。总体而言,研究区内相干性较高,基于升、降数据反演的形变速率精度可达毫米级,反演结果精度优于8～11 mm,升轨形变速率精度整体优于降轨,监测结果有效可靠。

表2 基于升降轨数据的滑坡隐患早期识别详情表Table 2 Details of Early Identification of Landslide Hazards Based on ascending and descending data

图8 基于升降轨数据的形变速率及隐患点形变特征图Fig.8 Deformation feature map of deformation rate and hidden danger points based on ascending and descending data

3.3.2 结合几何畸变的识别结果差异性分析

星载SAR数据的几何畸变与雷达成像侧视角、地形坡度坡向有着直接的联系,受几何畸变影响而导致的相位丢失或重叠是雷达卫星侧视成像的先天不足。从所用SAR数据几何畸变分析结果可知,畸变类型主要分为主动叠掩、被动叠掩,主动阴影、被动阴影以及叠掩阴影,主、被动根据观测地物本身产生畸变或因周围地物影响产生畸变而划分。升降轨数据中几何畸变主要发生在隆务河东西两岸,区域内受几何畸变影响较降轨数据更为明显。经统计,升降轨数据主要受到叠掩畸变的影响,且升轨受影响范围比降轨更广(图10)。

图9 形变速率精度反演结果与统计图Fig.9 Inversion results and statistical chart of deformation rate accuracy

图10 升降轨数据几何畸变分布及统计图Fig.10 Geometric distortion distribution and statistical chart of ascending and descending data

为直观反映升降轨几何畸变与坡度坡向的关系,以黄乃亥乡为例,升轨数据中叠掩区域主要分布于隆务河东西两侧东南朝向的斜坡上,阴影区集中分布于坡度较大的区域;降轨数据中叠掩区域主要分布在靠近隆务河西侧的东北朝向斜坡上,阴影区则分布在背向雷达信号的陡峭斜坡上。此外,叠掩区主要为坡向复杂区域,阴影区则为坡度较大区域,难以获取有效的监测信号,几何畸变的类型及其分布在升降轨数据中表现差异较大(图11)。

为进一步说明几何畸变对升降轨识别结果的影响,以D-04号隐患为例,该坡体在升轨数据中受到叠掩影响,形变速率反演结果中存在明显的信号丢失,形变速率密度及精度均较低;而在降轨数据中虽受到轻微阴影影响,但得到了有效且完整的形变速率结果,反映出了明显的形变特征(图12),这表明升降轨数据在形变解算时具有互补性,联合不同入射角的SAR影像获取不同成像几何下目标地物的形变特征,在一定程度上弥补了单一成像方式带来的监测盲区,识别结果更加全面,从而提高了滑坡隐患早期识别结果的准确性和可靠性。

4 典型滑坡隐患特征解析

在滑坡监测中,地表形变的时间和空间变化是评价滑坡阶段的重要指标,通过SBAS-InSAR时序方法获取到滑坡隐患不同时期的形变量级与速率参数,通过形变特征分析为滑坡监测与预警提供可借鉴的参考依据。现选取西山滑坡群、羊智滑坡群、尕沙日滑坡作为典型滑坡隐患进行特征解析。

图11 黄乃亥乡升降轨数据几何畸变分布与坡度坡向关系图Fig.11 Mapping of geometric distortion distribution and slope aspect of ascending and descending data in Huangnaihai Township

图12 D-04号隐患升降轨数据几何畸变与形变速率关系图Fig.12 Relationship between geometric distortion and deformation rate of ascending and descending data at D-04 hidden danger

4.1 西山滑坡群特征分析

西山滑坡群位于同仁市隆务镇夏琼南路西侧斜坡地带,二郎寺西侧,隆务河西岸,中心经纬度坐标为102°00′E,35°31′N,由Ⅰ～Ⅻ号12个滑坡体组成,为典型的大型土质滑坡群。该滑坡群发育于临夏组地层上部、残坡积层及隆务河四级阶地上,以西山南段古滑坡为主,严重威胁坡下居民点及学校等市政设施,隐患风险等级较高[12]。

由覆盖西山滑坡群的斜坡向形变速率分布图可以看出,整个坡面中存在2处强形变区域,呈环形、漏斗状形态分布于Ⅱ号和Ⅴ号滑坡体,表现出持续变形趋势。在两处强形变分布区沿斜坡方向(西北向东南)绘制了形变速率剖面,其中,Ⅱ号滑坡中部的形变速率最大,达-13 m/a;Ⅴ号滑坡的最大形变速率为-12 mm/a,位于坡体后缘,该滑坡后缘和中部的形变速率明显大于前缘,稳定性较差(图13)。经实地调查发现,Ⅴ号滑坡体后缘有3条雁形排列的拉张裂缝,中部发育6、7条裂缝,呈横向断续或连续发育,局部有下错现象(图14),这与InSAR形变监测结果相一致。

监测周期内西山滑坡群最大形变量为24 mm,最大形变中心位于Ⅴ号滑坡。在Ⅱ号和Ⅴ号滑坡分别从上部、中部、下部选取特征点进行时序形变分析,结果显示,在2019年6月23日之前形变呈线性匀速下降趋势,至2019年12月5日,整个坡体形变出现线性波动,表现出多次极快的小幅度下降和抬升,这与夏季降雨集中导致降雨入渗裂缝,表层土壤扰动有关,坡体稳定性变差。而在2020年1月—2020年5月18日,形变恢复平稳,呈缓慢下降,随着6月雨季的到来,形变明显加快,并在2020年9月之后进入新一轮的波动变化(图15)。

图13 西山滑坡群形变速率分布及剖面图Fig.13 Distribution and profile of deformation rate of Xishan landslide group

图14 Ⅴ号滑坡中部(左)和后缘(右)拉张裂缝实地照片Fig.14 Field photos of tension cracks in the middle (left) and rear edge (right) of No.V landslide

图15 西山滑坡群累积形变量分布及特征点时序曲线图Fig.15 Time series curve of cumulative deformation variables and characteristic points of Xishan Landslide Group

4.2 羊智滑坡群特征分析

羊智滑坡群位于黄乃亥乡东北部,地处隆务河西岸,海拔2 275～2 864 m,平均坡度42.5°,中心地理坐标为102°12′E,35°40′N,为典型土质滑坡区,所在坡体植被较少,黄土覆盖,南侧存在历史泥石流冲沟,坡下为谷地平原,分布有居民点和农田,严重威胁坡顶的羊智村,若发生滑动也会对隆务河产生影响,隐患风险等级较高[13]。

由形变速率反演结果可知,该斜坡范围内存在4个明显的强形变中心,自北向南编号为YZ_1#、YZ_2#、YZ_3#、YZ_4#,最大形变速率为-64 mm/yr,最大形变中心位于YZ_1#,平均形变速率为-13 mm/a。选取4条剖面线沿斜坡向下方向绘制了形变速率变化曲线,结果显示,各滑坡强形变区均位于斜坡中部或后缘位置,其中YZ_1#滑坡形变核心位于坡顶往下450 m处,形变速率达-64 mm/a,在滑坡体前缘古滑坡面范围内形变速率较小,变化范围在-20～-10 mm/a;YZ_2#与YZ_1#滑坡形变速率变化相似;YZ_3#滑坡形变速率相对较小,最大形变速率为-27 mm/a,形变核心位于坡体中部;YZ_4#滑坡形变核心位于坡体后缘。其中YZ_3#靠近村庄居民点,人类活动迹象明显(图16)。

图16 羊智滑坡群形变速率分布及剖面图Fig.16 Distribution and profile of deformation rate of Yangzhi landslide group

结合GF-2光学影像圈定了该滑坡群内各滑坡的边界范围,明确了强形变中心对应的滑坡体,各形变中心均位于滑坡体的中上部区域,坡体陡峭,植被稀疏,黄土覆盖,为滑坡灾害的发生提供了环境条件及物质来源,滑坡群区域内无断层发育,断裂运动影响较小。各滑坡在GF-2影像中特征明显,可清晰识别出YZ_1#滑坡前缘及中部有次级滑坡痕迹;YZ_2#中部发育有次级滑坡;YZ_3#滑坡后缘及中部有次级滑坡,沿斜坡方向自上而下发育有历史泥石流冲沟;YZ_4#后缘呈典型的“圈椅状”地貌,且发育有次级小型滑坡(图17)。

羊智滑坡群范围内最大累积形变量达-122 mm,位于YZ_1#滑坡。整体上累积形变量的分布特征与形变速率分布具有很好的一致性。针对4个滑坡体共选取了14个特征点绘制时间序列曲线图,其中,YZ_1#滑坡最大累积形变量接近120 mm,在监测时段内滑坡前缘形变较小,后缘则处于快速形变阶段,在2020年5月之后,因降雨增多形变量呈加速递增趋势,而斜坡中部处于缓慢变形阶段;YZ_2#滑坡整体形变特征与YZ_1#相似,形变集中于斜坡中部偏后缘,其形变量随着时间的变化逐渐增加;YZ_3#滑坡中部累积形变量明显大于前缘和后缘,整体形变量随着时间变化不断递增;YZ_4#滑坡后缘形变明显,在2020年之后其形变量变化趋势明显加快。整体而言,羊智滑坡群受夏季降雨影响存在较大隐患(图18)。

4.3 尕沙日滑坡特征分析

沙日村位于同仁盆地中的河谷阶地区,处于隆务河西岸,地理坐标为102°01′33″E,35°36′6″N。村庄后侧所靠山体为南北走向,海拔2381～3005 m,相对高差约为620 m,其山体中下部由水平产状的第三系红色砂泥岩构成,尕沙日滑坡上部覆盖有风积黄土,纵向坡度约65°。该滑坡发生于2009年9月13日,滑体前缘在河谷阶地上形成碎屑流冲击,直接损坏公路、房屋、引水渠等,并有人员伤亡,造成直接经济损失1 256.4万元。

尕沙日滑坡范围内最大形变速率达-41 mm/a,形变核心集中于坡体中部及后缘位置,监测时段内后壁处于不稳定状态。由于受到几何畸变的影响,后壁靠近滑坡边界位置失相干,存在监测空值。沿斜坡向下方向(G-G′)绘制的形变速率剖面图中显示,坡体中后部的形变速率明显大于斜坡前缘,即斜坡体后缘处于不稳定状态(图19)。

经GF-2影像解译发现,滑坡体后壁“圈椅状”特征明显,且发育有多处横纵向拉裂缝,最长裂隙约300 m,坡体后壁失稳存在滑动的可能性,这与时序InSAR的形变监测结果相吻合(图20)。

图18 羊智滑坡群累积形变量分布及特征点时序曲线图Fig.18 Time series curve of cumulative deformation variables and characteristic points of Yangzhi Landslide Group

图19 尕沙日滑坡形变速率分布及剖面图Fig.19 Distribution and profile of deformation rate of Gashari landslide

图20 尕沙日滑坡光学影像解译图Fig.20 Optical image interpretation map of Gasari landslide

该滑坡在监测周期内最大累积形变量达-79 mm,整体处于蠕变状态。沿斜坡方向选取5个特征点绘制了时序曲线,可以看出,除了靠近斜坡前缘的G4和G5两个点从2019年7月—2020年1月存在一个缓慢抬升的趋势,位于斜坡中部和后缘部位的G1、G2和G3的形变量随着时间的变化而不断增加。整体而言,斜坡中部处于不稳定的缓慢蠕动变形阶段,而后缘则处于剧烈变形状态,由于滑坡体后壁拉张裂缝的发育,在降雨的影响下,极易再次发生滑动(图21)。

图21 尕沙日滑坡累积形变量及时序位移曲线图Fig.21 Cumulative deformation and time sequence displacement curve of Gasari Landslide

5 结论

基于时序SBAS-InSAR技术,联合升降轨SAR影像从不同成像几何下反演获取了研究区地表形变信息,结合孕灾背景资料,辅以高分辨率光学影像筛选,对区内滑坡隐患进行了早期识别,得到如下结论。

(1)在2019年1月—2020年8月监测周期内,区内共识别出滑坡隐患19处,其中历史滑坡8处,新发现隐患11处。所有隐患集中分布于隆务河两岸和隆乌古曲河畔,涉及黄乃亥乡、年都乎乡、兰采乡、隆务镇、加吾乡、曲库乎乡、扎毛乡等7个乡镇。

(2)经形变速率精度分析,本次形变监测结果精度优于8～11 mm,监测结果有效;结合区内升降轨SAR影像几何畸变类型及分布,对识别出的19处隐患进行了相干性和几何畸变复核检验,所有隐患均处在非几何畸变区,形变监测结果可靠。同时结合几何畸变对升降轨识别结果进行了差异性分析,19处隐患中,升轨识别8处,降轨识别11处,联合升降轨数据从不同成像几何下对目标地物的形变特征进行提取,在一定程度上弥补了单一成像几何带来的监测盲区,识别结果更加全面,提高了滑坡隐患早期识别的准确性。

(3)通过典型滑坡特征分析发现,本次研究识别出的隐患均处于坡度陡、高差大的斜坡上,冲沟发育、溯源侵蚀现象严重,为其滑动提供了良好的地形条件;滑坡类型均为土质滑坡,所处黄土区,坡体结构松散、空隙大,具有较强的透水性,独特的岩土体条件极易发育形成滑坡灾害;坡体植被稀疏,多分布有地表裂缝,为地表水的入渗提供了通道,致使坡体稳定性变差;从典型隐患形变的时间序列分析结果可知,各隐患坡体均呈缓慢变形趋势,稳定性较差,在研究区5—10月,尤其是6—8月集中降雨期间,形变呈加速趋势,降雨是该区诱发滑坡灾害的主要因素。