基于时序InSAR技术的中贵天然气管道天水市段沿线滑坡隐患识别与形变分析

方迎潮, 王小松, 蒋毅, 顿佳伟, 冯文凯*, 刘威, 丁治文, 张洋铭

(1.中国科学院水利部成都山地灾害与环境研究所山地灾害与地表过程重点实验室, 成都 610041; 2.中国科学院大学, 北京 100049; 3. 国家管网集团西南管道有限责任公司, 成都 610041; 4.地质灾害防治与地质环境保护国家重点实验室(成都理工大学), 成都 610051;5. 中国石油西南油气田公司川西北气矿, 绵阳 621000)

管道运输是随着石油与天然气工业快速发展而产生的一种特殊的货物运输方式,具有运量大、运费低、占地少、不受气候限制、易于自动控制及可持续作业的优点[1-3]。中国山地、高原及丘陵地势较为突出,管道线路在铺设过程中不可避免地会穿越地形起伏较大区域。天水市位于黄土高原地区,滑坡隐患地质灾害高发,对天然气管道布控及运行带来重大威胁与严峻挑战,因此亟需对管道沿线开展地表监测、滑坡隐患地质灾害识别。

近年来发展的新型遥感技术合成孔径雷达干涉测量(synthetic aperture radar interferometry,InSAR),可在不影响管道运行的情况下对管道沿线滑坡隐患进行早期识别,同时了解周边地质环境,掌握滑坡灾害的活动范围及运动特征,可初步判断受影响管道与灾害体的空间位置关系。运用InSAR技术对滑坡隐患地质灾害进行早期识别,黄锐等[4]对中缅油气管道采用了升降轨融合的永久散射体合成孔径雷达干涉测量(persistent scatterer interferometric synthetic aperture radar,PS-InSAR)技术进行监测,同时结合研究区多源信息数据建立滑坡风险评估达到了早期识别的目的。孙伟等[5]通过小基线集时序分析(small baseline subset InSAR,SBAS-InSAR)技术对输油管道重点段分析,有效获取了管道周边地表形变。张磊[6]通过时序InSAR技术全面识别晴隆县段管道沿线滑坡隐患点,并根据识别结果对研究区进行区域滑坡隐患灾害危险性评价。淦邦[7]分析了管道沿线滑坡灾害破坏特点,通过分析管道穿越方式如横穿、纵穿和斜穿,分析滑坡与管道之间的空间位置关系,进而做出安全监测及评价。白路遥等[8]分析卫星遥感技术在滑坡识别和应用的思路,分析国内外管道行业对管线周边地质灾害识别所运用的卫星遥感技术,对卫星遥感在滑坡隐患的识别应用的发展趋势做了总结。周定义等[9]运用SBAS-InSAR技术对高山峡谷地区滑坡进行早期识别,利用升降轨SAR数据互补的方式,精准全面地对高山峡谷滑坡隐患进行早期识别。邹永胜等[10]运用InSAR技术同时联合全球导航卫星系统(global navigation satellite system,GNSS)数据验证,提出线性区域管道地质灾害监测预警体系,结合GNSS现场数据,验证了体系的准确性和时效性。以上多为学者并未将SBAS-InSAR技术运用在管道沿线滑坡隐患灾害,同时结合降雨-形变进行归一分析,并分析滑坡隐患发育分布规律。

目前,中国使用时序InSAR技术对管道沿线滑坡灾害识别与分析研究较少,采用2019年1月—2021年12月的89景升轨、90景降轨Sentinel-1A数据,利用SBAS-InSAR技术,对中贵天然气管道滑坡隐患灾害进行早期识别研究,运用InSAR技术,同时结合降雨数据对黄土高原研究区中贵天然气管道重点段(天水段)滑坡隐患进行早期识别研究。并对管道沿线滑坡隐患发育特征与典型滑坡形变进行分析。联合升降轨的时序InSAR技术可以有效识别管道沿线滑坡隐患,为油气管线的安全运营及今后油气管道选线提供重要的科学依据和参考。

1 研究区概况

研究区位于甘肃省天水市中贵天然气管道沿线2 km范围内,管道总长达110 km,区内面积约275 km2,管道由北向南贯穿天水市秦州区、麦积区和秦安县(图1),南部与陇南市接壤,北部与定西市相交。研究区内地势西高东低,海拔为1 000～2 100 m,北部为黄土沟壑区,土壤以棕壤、黑垆土、冲击土等主要的耕层土壤,中部为渭河平原地貌,南部为山地地貌。管线所经区域发育的岩性有第四系、第三系、泥盆系、震旦系以及元古界系等多套地层。研究区地处副热带北部,年平均气温分布在8～12 ℃,年降雨量在425.0～553.2 mm,所处地段为中国地形、气候过渡带,气候复杂多变。同时区内人口较为密集,交通便利,修建有多条高速和国道,为地质灾害监测重点区。

2 数据和方法

2.1 数据获取

Sentinel-1A卫星从运行至今在地质灾害领域得到了广泛的应用[11-16],故研究区选取Sentinel-1A升降轨影像作为数据源,数据时间跨度均为2019年1月—2021年12月,升降轨影像入射角分别为43.95°、34.03°,SAR影像覆盖范围如图2所示。

图2 研究区Sentinel-1A数据覆盖范围Fig.2 Data coverage of Sentinel-1A in the study area

DEM采用美国地质调查区(SRTM)30mDEM数据,用作提供地理坐标系,来消除地形相位的影响[17]。

2.2 时序InSAR处理方法

常用的时序InSAR处理方法有永久散射体合成孔径雷达干涉测量技术(PS-InSAR)和小基线及干涉叠加测量技术(SBAS-InSAR),由于SBAS-InSAR技术有效解决了SAR数据集之间基线过长导致的时间不连续的问题,有效降低时空去相干和大气延迟误差等影响[9-10]。因此采用SBAS-InSAR技术对中贵天然气管道沿线进行解译识别。SBAS-InSAR技术最早由Berardino等[18]提出,其原理主要为利用小基线干涉测量,选择主影像与多景辅影像进行干涉,得到多幅差分干涉图,通过去除地形效应、相位解缠后利用最小二乘法得到干涉相位的形变时间序列,然后将干涉相位值用奇异值分解法联合求解,获得高分辨率时序形变监测结果。

基本原理:假设有N+1景单视复数影像,其成像时间为(t0,t1,…,tn),假设某一差分解缠图为j,某像元(x,r)对应的解缠相位可表示为

(1)

(2)

(3)

(4)

式中:λ为波长;θ为雷达视角;ΔZ为DEM误差;R为雷达距目标物体的斜距;d(tA,x,r)和d(tB,x,r)分别为不同参考时间tA、tB的视线向累计形变量;B⊥j为垂直基线长度;M为干涉图数量。

因此,设干涉区间内的形变速率为vk,k+1,则tA～tB的累计形变可表示为

(5)

式(5)中:δφdef(x,r)为干涉相位。

对M幅干涉图进行奇异值分解相位解缠,即得到不同时相SAR数据的形变速率。

使用ENVI软件中SARscape模块,用SBAS算法对Sentinel 1A数据进行干涉处理。根据研究区特点及数据优势,分析所获取数据的最大空间基线及平均时间基线,将阈值设置为250 m和72 d,以此为条件获取干涉数据降轨368对,升轨375对,其中降轨数据以获取日期2021年1月29日为主影像,升轨数据以日期2019年3月17日为主影像。对所获取的干涉数据进行相位解缠,结合30 m SRTM DEM数据,对地形相位进行修正,并使用POD精密轨道参数对卫星轨道信息进行修正。经过轨道精炼与重去平、数据反演及SAR数据地理编码后,得到所研究区域形变速率值。

3 InSAR识别结果与分析

3.1 地表形变监测结果

运用升降轨的SBAS-InSAR的方法,开展中贵油气管线天水段左右2 km范围内的地表形变监测。对监测结果的升降轨形变点进行统计分析,从升、降轨数据中共提取研究区内各56.5×104、76.1×104个SBAS形变点,分别如图3(a)、图3(b)所示。

图3 研究区升降轨年均形变速率统计图Fig.3 Annual average deformation rate statistics for the study area for ascending and descending rails

对升轨、降轨SBAS点形变速率进行分析得到中贵管线周围地表形变情况。如图3所示,获得基于升降轨数据的SBAS-InSAR处理结果对其形变特征点进行统计分析:升轨降轨结果形变速率值大部分在-7.5～7.5 mm/a范围内,处于该范围的点位区段可视为基本稳定状态,因此设置研究区稳定阈值为±7.5 mm/a,在升轨结果中有38.7×104个形变特征点,占总个数的68.5%,降轨结果中56.3×104个形变特征点,占总数的74.0%处于稳定阈值中。图3(a)、图3(b)整体呈现出正态分布,其结果与预计情况较为符合。

3.2 滑坡隐患识别及野外验证

将地理编码后的处理结果(SBAS点)与DEM转山体阴影图相叠加(图4),得到监测时间段内中贵天然气管线天水段形变情况。从SBAS点分布情况(图4)可以看出,形变结果显示在研究区段内的绝大多数处呈现出绿色状态,少数区域呈现出深蓝色和红色。通过地表监测结果分析,设定形变阈值为±15 mm/a,即对管线周围形变值绝对值大于15 mm/a以上的形变区进行圈定分析,圈画结果如图5所示。经分析,管道地质灾害隐患区共17处。经现场调查,最终确定13处为滑坡灾害(表1),4处为人类工程活动导致的形变。

表1 黄土高原研究区InSAR形变区信息Table 1 InSAR deformation information of Loess Plateau study area

红色区域表示监测对象在视线方向上(LOS)远离卫星;深蓝色区域表示监测对象在视线方向上(LOS)靠近卫星;绿色区域表示监测段内无明显形变图4 SBAS点分布Fig.4 SBAS point distribution

HP01～HP13为灾害点名称图5 滑坡隐患灾害点分布Fig.5 Distribution of landslide hazard points

3.3 滑坡发育分布规律

对研究区内13处威胁管道的滑坡隐患灾害点进行发育分布规律分析。黄土滑坡灾害的分布一般与灾害区的地形地貌(坡度、坡向)、地层岩性等息息相关。由于本研究主要关注于管道沿线滑坡隐患灾害,因此在进行分析时加入了据管道距离因素,来进一步分析滑坡隐患发育分布规律。

3.3.1 距管道距离

从滑坡灾害分布[图6(a)]可以看出,研究区中13处滑坡中有10处主要分布在距离管道0～150 m区间范围内,占总灾害的76.92%,如图7(a)所示;2处滑坡距小于50 m,即HP03滑坡、HP13村滑坡,所占比例为15.38%;2处滑坡距管道距离为50～100 m,为小石沟地滑坡、杨家湾滑坡,所占比例为15.38%。距离管道距离大于150 m的灾害,整体规模较小,对管道整体影响较小。

图6 滑坡隐患与各项因素关系图Fig.6 Relationship between landslide hazards and various factors

图7 滑坡隐患与各项因素直方图Fig.7 Histogram of landslide hazards and factors

3.3.2 地形地貌

滑坡隐患的产生与坡度、高程有着紧密的联系,高程、坡度为滑坡灾害的发生提供了势能。经分析研究区内滑坡主要分布在50°以内,其中40～50°分布最多,占总灾害的61.54%[图6(c)],小于20°以及大于60°均不存在滑坡灾害,因为坡度过小时滑坡势能不足,无法为滑坡活动破坏提供条件。当坡度大于60°时,由于坡度太陡与SAR影像成像时易产生阴影、叠掩,无法形成有效SBAS点,对滑坡体的解译产生较大影响进而无法识别。

研究区内滑坡灾害前后缘高程要集中在400～600 m,共发育7处滑坡隐患[图6(b)],占总灾害的53.84%[图7(b)],前后缘高程即高差为滑坡隐患的发生提供了有利的条件,随着高差的增大滑坡隐患灾害数量逐渐减少。

研究区内滑坡灾害主要发育在坡向为270°～315°、270°、90°、45°～90°方向,其中45°～120°方向发育最多,占总灾害的69.23%[图6(d)、图7(d)],主要是卫星受自身飞行方向原因,使得对地表180°、360°(0°)方形变解译较为困难。

3.3.3 地层

4 典型滑坡形变分析

在解译出的13处滑坡隐患中,有2处滑坡隐患(HP03、HP13)距管道<50 m,且整体规模较大,一旦发生失稳,将对管道造成严重威胁。同时在对13处滑坡隐患进行调查的过程中发现,该研究区内地质灾害除了受自身条件(结构特征)以外,还受外界因素影响。近年来,研究区内降雨时间较为集中,这就导致各种地质灾害的发生概率大大增加,降雨一直以来就与各类地质灾害的稳定性息息相关[19-21],本次研究获取天水市天水气象站2019—2021年月降雨量数据,结合不同时间段时序InSAR技术,来分析HP03常沟村滑坡、HP13磨峪沟村滑坡不同阶段形变演化过程之间的关系。

4.1 HP03

HP03常沟村滑坡位于甘肃省天水市秦州区常沟村,此滑坡位于常沟村南部,滑坡边界较为清晰整体呈矩形,坡体斜长为180 m,宽为97 m,前缘至坡脚高程约为1 676 m,后缘至坡体中部平台陡缓交界处,高程约为1 728 m,高差为52 m,滑坡主滑方向为260°,坡度约为30°。如图8(a)所示在光学遥感影像上可见后缘及左侧边界明显下错,滑坡体较陡,坡体整体植较少,覆盖层较厚,以第四系黄土堆积为主,坡体透水性较差。该地区为典型黄土丘陵沟壑区,滑坡发育在黄土丘陵斜坡上,影响下易产生滑动。通过对时序InSAR结果的量化分析,从图8(b)中可以看出,滑坡体右侧形变较为明显,滑坡体整体形变速率在7～-19 mm/a,为了更加准确地分析滑坡体的形变特征,在滑坡体中部及靠近管线部位选取2个形变特征点,如图8(b)白色点标记所示。如图9所示,特征点P1、P2时序形变速率波动趋势大致相同,在2020年7月出现明显的形变加速(图9),累计形变量值最大达到-30 mm。

图8 HP03滑坡光学遥感影像图和InSAR结果Fig.8 Optical remote sensing image map and InSAR result of HP03 landslide

图9 HP03滑坡特征点形变速率图Fig.9 Characteristic point shape deformation rate diagram of HP03 landslide

绘制P1、P2特征点的时间序列曲线与常沟村的月降雨量变化图(图10),由于P1点、P2点形变趋势相似,故只分析P1点。如图10所示,P1点在2019年7月—2019年9月期间丰富的降雨下,由于雨水的渗透,导致产生季节性沉降,数值达到13 mm,2019年10月—2020年4月非雨季型,形变较为稳定;2020年5月—2020年8月,月降雨量150 mm达到近年之最,土壤吸水膨胀导致地表发生抬升至12.3 mm。2020年9月—2020年11月,累计形变较大达到了-25 mm,推测在此区域发生了滑坡的蠕变;在2021年1月—2021年7月非雨季,P1点有沉降趋势但整体形变较为稳定,在2021年8月雨季过后有些许形变,累计型变量达到-37.5 mm,整体呈现季节性波动。

图10 滑坡特征点形变与降雨变化图Fig.10 Deformation of landslide characteristic points and rainfall variation diagram

对HP03滑坡展开了野外实地调查(图11)。调查结果显示:滑坡体在平面形态上呈矩形,坡体土体较为破碎且大部分区域已被人工改造为耕地。坡体后缘的下错陡坎清晰可见,在坡体中部靠近左边界处存在剪切裂隙,延伸6～8 m,滑坡体上植被覆盖较为稀少,覆盖层较厚。该滑坡体主要为粉土,主要为浅黄色粉土,均质结构,抗剪强度低,发育虫洞、根洞,岩土体盐水率较低,处于可塑状态;该滑坡地层主要为第四系上更新统风积黄土层(Qeol3),为浅黄色粉质黏土,土体孔隙较大,结构较为疏松,节理垂直发育,土体松散破碎,降水易入渗,可能会引起滑坡整体性失稳。管道从滑坡体前缘穿过,若滑坡体发生失稳会对管道造成一定影响。经分析常沟村滑坡裂缝发育,应及时处理,以防雨水入渗使得滑坡整体失稳对管线造成影响。滑坡后缘存在下挫陡坎,应尽可能避免人类工程活动。

图11 HP03滑坡现场调查图Fig.11 Site survey drawing of HP03 landslide

4.2 HP13

HP13磨峪沟村滑坡位于甘肃省天水市秦州区磨峪沟村附近。如光学遥感图所示[图12(a)],可见滑坡后缘及左右边界有明显下错,滑坡前缘堆积清晰可见,滑坡剖面形态为凹形。滑坡前缘高程为1 460 m,后缘至陡缓交界处,高程为1 513 m,坡体斜长为120 m,宽为110 m,坡度为45°～55°。InSAR结果形变异常区域主要集中在滑坡体中部,最大形变速率达-19 mm/a[图12(b)]。

图12 HP13滑坡光学遥感影像图和InSAR结果Fig.12 Optical remote sensing image map and InSAR result of HP13 landslide

对HP13滑坡选取P1、P2点2个样点形变进行时序分析。从图13中可以看出,2019—2021年每年降雨量最大集中在6—9月份,其他月份均有降雨,但降雨量较小,相对来说不是很集中。从累计形变量变化趋势来看,P1、P2两处特征点的累计形变量曲线的变化趋势大致相同,不同点位之间的形变量有所不同。两处特征点的历史形变量分别经历三次与降雨变化的响应阶段。以P1点最大累计形变量(35 mm)曲线为例进行降雨分析,P2特征点的累计形变量曲线与降雨的相关性分析与P1点基本一致。

图13 滑坡特征点形变与降雨变化图Fig.13 Deformation of landslide characteristic points and rainfall variation diagram

从P1点的时序累计形变量曲线可以看出,在2019年1月—2019年4月时间段,该区域存在一定的降雨量,最大降雨量为45.2 mm,该时间段内形变量呈缓慢增长的趋势,累计形变达到-10.3 mm。2019年7月降雨量加大,P1点的累计形变随之加大。随后在2019年10月—2020年3月时间段降雨量减少,形变增长速度减小。

在2020年7月降雨量与同年7月相比急剧增加,P1点在2020年6月之前的形变量曲线斜率小于2020年6月之后的形变量曲线斜率。说明6月之后形变量增长速度大于6月之前的形变增长速度。从图14中可以看出,P1点累计形变量曲线2020年7月处存在明显的沉降漏斗。此时P1点、P2点的行遍知分别达到-18.35、-19.23 mm。

图14 HP13滑坡现场调查图Fig.14 Site survey drawing of HP13 landslide

2021年9月降雨量再次增加,P1点累计形变量曲线斜率同2020年9月之前曲线斜率相比再次有所增加。P1点、P2点形变值达到-35.38、-33.48 mm。

通过上述降雨分析,结合特征点时序曲线的变化情况,P1、P2两处特征点的时序形变加速趋势与降雨量变化高度吻合,每年6—9月时间段内,形变量显著加速。由此可见HP13滑坡体上局部形变与降雨呈正相关。目前滑坡处于蠕滑阶段,极有可能在强降水作用下再次发生复活,并对管线造成一定影响。

经现场调查,发现滑坡后缘及侧壁陡坎明显[图14(a)],后缘下错陡坎约10 m,可见表面风化壳与其下紫红色泥页岩明显分界,后缘拉张裂隙发育;滑坡体中部裂缝较为发育前缘为季节性河道,坡体前缘多有石块堆积,堆积体受风化作用影响较大,局部可见强变形迹象[图14(b)]。滑坡目前处于破坏阶段。坡表变形迹象强烈。在降雨作用下易进一步扩展,中贵管线位于滑坡前缘12 m处,滑坡对管线造成一定威胁。

5 结论

将SBAS-InSAR技术运用于中贵天然气管道天水段滑坡隐患识别与分析,并分析典型滑坡的形变特征,得出如下主要结论。

(1)研究区内采用SBAS-InSAR技术,结合形变阈值(15 mm/a)一共识别出17处威胁管道隐患点,通过野外调查,确认13处为滑坡地质灾害,4处为人类工程活动造成,野外形变迹象与InSAR解译结果较为吻合,即运用SBAS-InSAR技术可有效识别管道沿线滑坡隐患。

(2)对研究区内滑坡灾害发育分布规律进行了分析,滑坡多发生的坡度40°～50°,前后缘高差多为400～600 m,坡向45°～120°,与管道距离100～150 m范围内;研究区内地层复杂多样,岩体结构较为破碎,滑坡多发育在第四系全新统冲洪积层。

(3)将典型滑坡特征点形变值与降雨数据进行分析,得出HP03常沟村滑坡以及HP13磨峪沟村滑坡的形变与降雨成正相关。两处滑坡处于蠕滑阶段,在降雨作用下易复活,对管线造成一定影响。