郑万基, 孙 倩, 刘小鸽
(1. 中南大学 地球科学与信息物理学院,湖南 长沙 410083;2. 湖南师范大学 资源与环境科学学院,湖南 长沙 410081;3. 湖南师范大学 地理空间大数据挖掘与应用湖南省重点实验室,湖南 长沙 410081)
金沙江位于长江上游,流经中国西藏、云南及四川等地,流域内坡陡水急,落差高达5100 m,蕴含丰富的水能资源,是中国具有重要战略地位的水电基地。然而,金沙江流域内地形情况极为复杂,气候多变,易发生滑坡及泥石流等自然灾害。伴随着流域内水库蓄水,库区内的地下水会发生相应的变化,使得稳定的堆积体有再次发生滑动的可能。库区内发生滑坡,可能导致建筑物发生损毁、水库造成淤积,形成涌流,给经济生产建设带来极大的安全隐患,因此需要进行长期的地表形变监测。目前常用的监测技术主要包括GPS、水准等地面测量手段,然而这些传统手段已经难以满足大范围、高精度、高时空分辨率的监测需求,并且需要投入大量人力物力,作业人员的人身安全难以得到保障。近些年发展起来的时序合成孔径雷达干涉测量(InSAR)技术,在地震[1]、火山[2]、滑坡[3]及城市基础设施监测[4]等方面已经得到较为广泛的应用,该技术弥补了传统监测手段空间分辨率低的不足,无需投入大量的人力物力,兼具较高的监测精度,其cm到mm级的监测精度已得到广泛的验证[5-6]。
目前使用较为广泛的时序InSAR监测方法主要是SBAS(小基线集)技术[7]及PSI(永久散射体InSAR)技术[8]。然而SBAS技术难以解决大气延迟等与空间相关的误差,且需要进行相位解缠,在面对一些“孤岛”区域时,易产生解缠误差;而PSI技术则需要监测对象在时间序列上保持较高的相干性,否则将难以选取较多的形变监测点,降低InSAR监测的空间分辨率。TCP-InSAR技术是近年发展起来的一种新型时序InSAR技术[9],该技术融合PS技术与SBAS技术的特点,在无需解缠的情况下,可以有效地克服与空间相关的误差,并确保可以选取较多的监测点,保证InSAR监测的空间分辨率,其监测精度在滑坡[3]、沉降监测[10]及火山活动监测[11]等领域已得到验证。
本文基于欧空局的新一代中高分辨率C波段SAR卫星——哨兵一号A星(Sentinel-1A),使用TCP-InSAR技术,获取金沙江中游兴培当至草可都段的滑坡结果。最后根据监测结果,结合从NOAA获取的气象资料及当地地质情况,分析6处滑坡的影响因素。
本文的研究区域如图1所示,地处云南省西北部,距离丽江市以北约69 km处的金沙江中游地区,该河段距离下游的梨园水电站约28 km。如图1中(c)所示,金沙江连通区域南北,格吉河与白水河于此处汇入金沙江,右岸有草可都、兴培当等居民点,左岸则为下只恩水电站和白水河水电站,即从北纬27°27′37.2″至27°30′15.4″,东经100°8′45.1″至100°12′22.7″的范围。该地区处于滇西北高原地区,地势起伏较大,山川走向以北北西和北北东为主。研究区内气候分为干、湿两季,湿季主要集中于5—10月,年降水量约为900 mm,在雨季集中降雨时,常发生山洪及泥石流等灾害。该地区地层岩性以玄武岩及灰岩为主,存在有多处活动中的大型堆积体滑坡,滑坡体上分布有少量的耕地和居民,滑坡体物质以碎块石土、残坡积层、粘土及砂土为主[12]。本文通过TCP-InSAR技术对该区域的滑坡体进行监测,而InSAR技术通常会受地表覆盖植被的影响,影响大小与植被高度、SAR信号频率有关[13]。研究区海拔均处于2 500 m以上,没有高大的树林遮挡,而由图1(c)可知,在第一景SAR影像时间点前后,该区域植被覆盖较少,个别有植被覆盖的区域植被也较为稀疏,不会对C波段的哨兵一号卫星造成失相关。此外,本文采用的TCP-InSAR技术,将根据时间序列上干涉图的相干性进行选点,若存在低相干性的植被观测点,将不会被选中用于形变估计。
实验采用的数据为欧洲空间局的新一代中高分辨率卫星哨兵一号A星(Sentinel-1A)数据。该卫星为C波段的合成孔径雷达卫星,采用全新的TOPSAR成像模式,实现在中等分辨率(距离向3.7 m,方位向14 m)下,具备12 d的重返周期及250 km×250 km的大范围观测能力。在A星与B星串联模式下,更可实现6 d的重返观测。此外,基于欧盟的哥白尼计划,该卫星目前可免费提供全球最新的实时观测数据,这为实现高时空分辨率的滑坡监测提供重要的数据保障。实验从欧洲空间局的哨兵科学数据中心处获取17景Sentinel-1A数据,时间跨度为2015-12-18—2017-02-10,覆盖区域如图1(b)中红色方框所示。由于哨兵卫星具备较好的轨道控制能力,因此此次获取的哨兵数据其垂直基线皆控制在150 m以内,具体数据参数如表1所示。
注:图(a)(b)底图为当地DEM数据,图(c)底图为google获取的2015-11-22影像图1 研究区位置及数据覆盖情况
序号时间升/降轨入射角/(°)垂直基线/m时间基线/d12015-12-18升轨33.90022016-01-11升轨33.926.92432016-02-04升轨33.94.84842016-02-28升轨33.9-67.57252016-04-16升轨33.9123.712062016-05-10升轨33.9-71.114472016-06-03升轨33.9-18.216882016-07-21升轨33.950.821692016-08-14升轨33.9-44.4240102016-09-07升轨33.9-10.5264112016-10-01升轨33.913.6288122016-10-25升轨33.9-50.3312132016-11-18升轨33.9-42.3336142016-12-12升轨33.951.8260152017-01-05升轨33.98.3384162017-01-29升轨33.9-14.7408172017-02-10升轨33.916.6420
TCP-InSAR技术是近几年发展起来的一种新兴的时序InSAR技术[9],该技术是在早期PSI和SBAS等常用时序InSAR技术的基础上,充分结合二者的优势发展出来的。在干涉图组合上,TCP-InSAR技术使用SBAS技术中的多主影像、短时空基线的干涉数据组合方法,有效减少干涉图中的时空失相干噪声;而在观测量的选取上,该技术则使用类似于PSI技术的以构网弧段作为观测量的方法,通过相邻点互相差分的方式,有效地抑制了大气等空间相关的噪声。而与二者不同的是,该技术在参数估计时,除估计线性形变速率及高程残余外,还同时估计由于基线不准确导致的轨道趋势项[14];并且,该技术选取的相干点,并不需要在整个时间序列上保持高相干特征,该方法根据数据配准时的偏移量信息或相干性信息,选取出在整个时间序列及在部分序列上保持高相干的点,极大扩展观测点的数量[15]。此外,TCP-InSAR技术的另一大优势是,该方法不需要对差分干涉图进行解缠处理[16]。基于大部分相邻TCP点在短时空基线的差分干涉对上不包含整周模糊度的思想,TCP-InSAR技术在进行第一次参数估计之后,将根据观测残差判断出具有整周模糊度的弧段,并将其剔除。TCP-InSAR技术的技术流程如图2所示。
图2 TCP-InSAR数据处理流程
实验首先使用GAMMA软件将所有的SAR数据配准于同一空间坐标系下,为了避免TOPSAR模式的Sentinel-1A数据在配准时在不同的burst上产生的相位跳变现象,实验过程中利用USGS提供的30 m分辨率的SRTM数据用于辅助SAR数据配准。在完成配准后,将试验区域单独裁出,进行差分干涉处理,以减少数据的运算量,在处理过程中使用5×1的多视视数(最终影像分辨率约为14 m×14 m)。差分过程使用USGS提供的30 m分辨率的SRTM数据作为外部DEM去除地形相位。为了避免产生孤立的小基线集,降低解算结果的稳健性,同时还需保持较高的相干性以抑制时空失相干噪声,以80 m作为空间基线阈值;此外,对于滑坡而言,当滑坡的累积形变超出InSAR的监测梯度后,会造成失相关,影响InSAR监测结果精度,因此以90 d作为时间基线的阈值,由此最终选出如图3所示的24对干涉基线构成短基线集。此外,在滑坡地区,短时间内滑坡可能会因为崩塌而对InSAR信号造成严重的失相关,影响TCP-InSAR的监测质量,为此,以相干性0.4作为阈值选取出滑坡上未崩塌的形变点,最终在该地区沿金沙江两岸共选出33 987个TCP观测点构成观测网络。由于该试验区域选取的TCP点位在金沙江两岸,中间有金沙江隔断,为了使结果具有较高的稳健性,将TCP点分成东岸TCP点及西岸TCP点单独进行处理。在西岸,通过第一次参数估计,获得每个弧段的线性形变速率及高程残差,基于第一次参数估计的结果,以1.5作为TCP-InSAR模糊度探测的阈值[15],剔除具有模糊度的观测弧段,同时将轨道趋势误差从每个观测弧段中去除,并进行第二次的参数再估计,最终得到线性形变速率及高程残余,将线性形变速率加到每个SLC的残余相位中,得到时间序列上的形变结果。求得每个弧段的线性形变速率及时序形变序列后,选取如图1(c)所示金沙江西岸的下之恩水电站作为稳定的解算起始点,进行网络平差,求得每个TCP点的年均形变速率及时序形变结果。在金沙江东岸,以同样的参数进行处理,即得到东岸TCP点的线性形变速率、高程残差及形变时间序列,在东岸选取的解算起始点为远离形变区、地形较高且地势平坦的高相干点,如图1(c)金沙江东岸所示。
图3 TCP-InSAR构建的小基线集
图4 研究区雷达LOS向年均形变速率场,底图为google earth 2015-11-22影像
如图4所示,通过TCP-InSAR技术最终获取金沙江中游兴培当至草可都段的雷达视线向(LOS向)的年均地表形变信息。LOS向形变为地表真实形变在雷达实现方向上的投影,靠近卫星为正值,远离卫星为负值。由于目前多数星载SAR卫星均以极地轨道的方式运行,飞行方向近似南北向,因此LOS向形变中以东西向和垂直向投影结果的贡献最多。从图4中可以发现,该地区有多处滑坡区域,其中对金沙江会产生直接影响的滑坡则主要有4处,如图4红色虚线区域所示,形变方向主要以由两岸向金沙江滑动为主,按照形变量的大小,分别为位于格吉河口与白水河口之间的滑坡A、草可都滑坡B、白水河口右岸滑坡C、兴培当滑坡D。由于本次实验采用的是升轨SAR数据,因此,在金沙江西侧形变为负值,在金沙江东侧形变为正值。位于金沙江东岸的滑坡体为兴培当滑坡B和草可都滑坡D,草可都滑坡B的形变量约为25~50 mm/a,由于形变较大的缘故,在滑坡的中后部出现形变超出监测梯度的情况,因此,在草可都滑坡的中后部无法选出相干性较高的点;兴培当滑坡D的形变量约为15~27 mm/a,二者均为厚层大型堆积体滑坡。在金沙江西岸的主要滑坡体则为白水河口滑坡C,及位于格吉河口与白水河口之间的滑坡A。其中,滑坡A的形变量级为-40~-65 mm/a,而滑坡C为厚层中型堆积体滑坡,形变量级在-35~-47 mm/a之间。
根据Brabb等人的实验结果,地质、土壤及坡度是影响滑坡稳定性的主要原因[17]。草可都滑坡B的前缘坡度为50°~65°,中后部约为10°~30°,而兴培当滑坡D的前缘坡度为45°~75°,中后部为5°~25°,因此滑坡坡度较大的草可都滑坡的形变量比兴培当滑坡的形变量要大。此外,根据现场考察资料[12],在草可都滑坡体的两侧有发育小冲沟,前缘有泉水点出露,从而导致该地区土壤含水量增加,孔隙压力变大,有效正压力因此减小,使滑坡速度加快,也是使得草可都滑坡速度比兴培当滑坡速度大的原因。并且,草可都滑坡的前缘出现鼓胀裂缝及崩塌的现象[12],因此在实验结果中,草可都滑坡的前缘出现失相干现象。对于滑坡A所在区域,其坡度前缘为5°~10°,中后部约为10°~25°;而对于滑坡C所在区域,其坡度前缘为15°~35°,中部为5°~15°,后部为25°~45°。虽然滑坡A的地形坡度小于滑坡C,然而滑坡A的地表滑动速度却大于滑坡C,根据现场资料显示[12],滑坡A的堆积物以砾石与少量粉土组成,而滑坡C的堆积物为碎块石土,粉土塑性指数更大,因此滑坡A的形变更大。
为了进一步了解4个形变区域的动态变化过程,通过TCP-InSAR技术得出不同时间序列上的LOS向形变场。图5显示从2015-12-18—2017-02-10共16个时刻的形变。由时序结果中可看出,在滑坡A,B,C和D虽然在不同的时间形变的速度变化不一,然而主要形变仍以线性形变为主,随着时间的推移形变量在不断地增加。此外,在个别区域出现明显的季节性形变,如图中的P5点和P6点。
图5 时序形变场,图中P1~P4点为形变区域A~D中的一点,P5~P6为出现季节性形变区域中的一点
其中,P6点位于草可都滑坡与兴培当滑坡之间,其形变在2016年7月达到约25 mm后开始下降,而后在2016年9月再次达到25 mm,且形变范围更大;而P5点从2016年2月开始形变逐渐增大,在2016年8月达到最大35 mm的形变后,开始逐渐下降。两点的最大形变集中于2016年的夏季,而此时属于该地区的湿季,是降水最多的时期,推测这两个地区的形变与降水有关。
为了验证P1-P6点与降水之间存在的关系,从NOAA(美国国家海洋与大气管理局)获得该地区以南69 km的丽江市在对应时间段内的降水与气温数据作进一步分析,其结果如图6所示。气温为每日气温均值,降水量为SAR时间序列间隔内的总降水量,形变量为LOS向累计形变量的绝对值。该地区在夏季主要受海洋性西南季风和东南季风的影响,降水充沛,而冬季则受西风带气流影响,天气干燥。
图6 不同形变点累计形变与降水/气温之关系
如图6所示,位于白水河口右岸的滑坡P3点和位于兴培当滑坡的P4点的形变主要以线性形变为主,随时间迁移形变量不断递增,仅在降水较为集中的6—8月出现较为轻微的滑坡速度加快的现象,整体滑动速度则较为稳定。而对位于白水河口与格吉河口间的不稳定斜坡P1点和位于草可都的P2点而言,其形变主要是持续的线性形变与季节性形变共同作用的结果。两点在大体上处于形变不断递增的趋势,而在2016年6月至9月期间,即降水最集中的时间,P1点和P2点的形变量都出现较大增长。滑坡的发生一般与土壤渗透量有直接相关,渗透量一般由土壤入渗率与降水获得,因此此处结合相关文献中给出的地质土壤资料,结合降水对其进行分析。根据相关资料显示[12],P1点和P2点地区的滑坡物质以粘土、粉土为主,相比起P3点的碎块石土和P4点的碎块石砂土而言,P1点和P2点的土壤环境更为松散,由大气降水带来的地表水在这种土壤环境下更容易渗透,从而导致孔隙压力增大,致使有效正压力减小,从而使P1点和P2点比P3及P4点更易受到降水的影响。
相较于前4个点,P5点和P6点虽然在年均形变速率上较小,然而却表现出十分明显的季节性形变现象。P5点位于格吉河口与白水河口之间,是P1点所在山峰的另一侧。该点的形变主要由降水所引起,在降水集中的6月至10月期间该点的形变值在不断增加,随后开始递减,在12月出现少许降水时,该点又开始活跃,而后迅速递减,与降水表现出较强的相关性。P6点是一个堆积体,位于含可地区。该点的形变虽然也呈现出明显的季节性趋势,然而与P5点不同的是,该点的趋势与气温变化的相关性更大,虽然该点海拔只有1700 m,然而在该点附近有海拔4000 m以上的雪山,因此,推测该点的季节性主要由于气温的升高或降低,会导致雪山融雪或结冰,引起地下水位变化所导致的,据此,推断该点的形变主要是气温与大气降水共同作用的结果。在P1-P6点的形变时间序列中,除了P5点外,其它点在降水极少的2017-01-05—2017-01-29期间,出现形变速度较快的异常现象。由于这5个点都处于沿江地带,初步推测有可能与库区水位变化有关,然而由于缺乏当地水文数据的缘故,此处不作进一步讨论。
本文以云南省丽江市北部的金沙江中游兴培当至草可都段的滑坡作为研究对象,基于新一代中高分辨率卫星哨兵一号A星,采用TCP-InSAR技术,获取该地区2015年12月至2017年2月的LOS向形变信息。形变结果表明,该地区山体总体处于不稳定状态,多处山体出现滑坡的现象,在沿江的草可都地区,LOS向年平均形变速率最大达到了约55 mm/a,在形变大的区域更出现形变超出监测梯度的现象。然后,根据地形及当地地质情况,分析该区域内不同滑坡形变速率不同的原因。并由NOAA提供的气象数据,分析该地区共计6个主要滑坡体的影响因素,其中,草可都滑坡与位于格吉河口与白水河口之间的不稳定斜坡的形变速率均由降水与线性形变共同影响,较为稳定的兴培当滑坡和白水河口右岸的滑坡主要受线性形变的影响,而位于格吉河口与白水河口之间的不稳定斜坡另一侧的滑坡对于降水则较为敏感,位于含可的堆积体形变则主要与气温有关。
然而,该地区处于梨园水库的上游地区,滑坡极有可能与水库水位存在一定的联系,由于水文数据的缺乏,本文暂时不能对其进行分析,而真实地表湿度数据缺失,也会使得本文的分析可能存在一定的不足。此外,本文只获取了该地区的一维LOS向的形变结果,仅一维的形变结果,难以对滑坡结果有充分详实的理解[18]。因此,未来的工作将可以考虑融合多孔径InSAR( MAI)技术或偏移量跟踪技术[19],获取更加全面的滑坡信息。