巨淑君, 高志良, 蒋明成, 高强

(1. 成都理工大学环境与土木工程学院, 成都 610059; 2.国能大渡河流域水电开发有限公司, 成都 610095; 3.大连理工大学建设工程学部, 大连 116081)

大渡河流域位于中国西南山区,处在青藏高原边缘向四川盆地的过渡地带[1],多为高山峡谷地貌,自然地质条件错综复杂,蕴含大量活跃的地震带,地质结构极其脆弱,滑坡、泥石流、地震等灾害频发,近十多年来就发生过多起重大地质灾害事故,较典型的有2008年的汶川地震、2017年九寨沟地震和2022年的泸定地震[2],并伴随有大量的滑坡等次生灾害[3-6]。据统计,大渡河地区已发现的滑坡、泥石流等地质灾害隐患点超过了5 000多处[7-8],作为中国重要的水电资源开发基地,广泛分布的地灾隐患区极大威胁了水电设施和人员安全[9],制约了流域内社会经济发展。对大渡河地区开展大范围地灾普查和早期识别,可以最大程度地减小灾害造成的损失,对于区域的防灾减灾工作具有重要意义。

目前,该地区的地质灾害识别手段主要依赖于地面测量方法,采取人工排查、无人机摄影测量和GNSS (global navigation satellite system)基站监测等技术,对重点路段和点位进行监视。然而,由于流域地质环境恶劣,植被覆盖茂密,这些隐患点具有隐蔽性强、分布范围广、调查难度大等特点,给地面监测手段带来极大挑战,当前的技术手段无法满足大范围、长时间尺度的地质灾害监测识别。随着空间遥感技术的不断发展,合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)技术以其具有高空间分辨率、高测量精度、长期监测等特点,在地质测量领域受到广泛关注。InSAR技术通过分析一组多时相SAR影像,可以从中提取毫米级精度的地表形变信息,以较低成本实现对观测区域的大规模、全天候、高频次测量,特别是在大范围地质灾害监测领域表现出极大的应用潜力[10],被广泛应用在中国西南复杂地形区域的灾害监测领域。文献[11]利用DInSAR (differential InSAR)和时序InSAR技术分析了九寨沟地震后的区域滑坡情况;文献[12]采用Stacking技术对金沙江地区滑坡群进行识别;文献[13-17]采用了使用最为普遍的SBAS(small baseline subsets)技术对西南山区开展了大范围的滑坡灾害研究,取得了较多成果。然而SBAS技术的多视处理会降低分辨率,尤其是在使用Sentinel-1这类中等分辨率数据,可能会丢失许多局部滑坡区域,可以结合PSI(permanent scatterer InSAR)技术和SBAS技术特点,应用于滑坡监测中。

然而当前对大渡河流域滑坡研究主要集中在丹巴、甲居和泸定等个别隐患区,缺少对大渡河全流域的滑坡灾害普查。基于此,现针对大渡河流域复杂地形和气候条件下的地灾隐患监测开展研究,探索利用时序InSAR技术开展大渡河流域地质灾害全面、定期排查的可行性。收集了2018—2021年的Sentinel-1超过1 TB的数据,采用PS-SBAS全流程自动化数据处理技术对海量SAR影像数据进行了深入分析,实现对大渡河流域丹巴至乐山段完整区域大范围地表形变监测,首次获取了大渡河全流域历史灾害隐患分布结果。

此外,在过去的研究基础上,在该流域又识别出了10多处新的滑坡隐患区域,并对流域不同区域滑坡诱因进行了总结。本文成果可以为流域内大坝和边坡变形机理分析和灾害预警提供可靠数据和信息支撑,有效提升该地区灾害治理水平,同时方便后续对比分析泸定地震前后的流域地质灾害情况。

1 研究区与数据介绍

1.1 研究区概况

研究区位于四川盆地南部的大渡河流域,由北向南流经金川、丹巴、泸定等地,至石棉折向东,构成“L”字形河流走向,整个流域狭长地带共1 062 km。区域内地形起伏极大,相对高差达到2 000～3 000 m,岸坡稳定性差,遍布高陡边坡与危岩体,是中国地质灾害重灾区。随着该地区经济快速发展,大规模工程建设活动的增加,进一步加剧灾害发生频率。各类突发性的滑坡等灾害,不仅可能引发河流堵塞,同时也会阻断流域内唯一的通道S211省道,威胁水电站安全和区域施工建设活动,因此大渡河地区对于大范围地质灾害监测需求较为迫切。本文选取了大渡河丹巴至乐山段作为研究区,覆盖了大渡河流域主要区域,如图1(a)所示,其中蓝线表示大渡河;图1(b)为丹巴附近的光学影像图,可看到区域谷深坡陡,是典型的高山峡谷地貌。

1.2 数据介绍

目前能够满足大渡河流域形变监测应用的SAR影像主要是国外的Sentinel-1A数据,该卫星数据是由欧空局提供的C波段(波长5.6 cm)中等分辨率数据,数据来源(https://scihub.copernicus.eu/)。按轨道号排序为26-1、26-2和128-1,如图2红色框所示。每个片区获取了102景SAR影像,影像拍摄时间为2018年1月—2021年6月,极化方式为VV(vertical transmit, vertical receive)极化,TOPS(terrain observation with progressive scans)成像模式IW(interferometric wide swath)类型数据,超过300景的影像总数据量达到了1.25 TB,具体参数如表1所示。此外,还需要下载Sentinel-1数据对应影像时间的POD(precise orbit determination)精轨文件,获取地址为(https: //s1qc.asf.alaska.edu/aux_poeorb/),以修正轨道参数信息,以及时序InSAR处理过程需要使用的外部高程数据,采用SRTM(shuttle radar topography mission) 30 m DEM(digital elevation model),数据来源(http://srtm.csi.cgiar.org/srtmdata/)。

表1 研究区SAR影像数据详细参数Table 1 SAR image description used in this study

2 研究方法

传统的基于永久散射体PSI技术在大渡河区域,在时间基线较长时会出现严重的时间去相干误差,导致较大形变解算误差;而小基线集SBAS技术在复杂的西南山区,大量的阴影和叠掩区域的存在,会导致形变估计结果出现较大的分块跳变误差。因此,为了研究2018—2021年较长时间跨度的大渡河流域的地表形变信息,采用PSI和SBAS相结合的PS-SBAS方法,通过选取永久散射体PS,并基于小基线干涉对组合获取PS点对应的差分干涉相位,有效避免了长时间基线干涉对,同时不受阴影等低相干区域影响,最后采用SBAS处理思路完成数据处理,基本原理如下。

对一组获取时间为{t1,t2,…,tN}的SLC(single look complex)数据,基于幅度稳定性方法识别出M个具有稳定散射特性的PS点,{PS1,PS2,…,PSM}。根据区域干涉图质量,选择合适的时间基线和空间阈值,同时可以依据相干系数阈值,筛选得到一组具有K个干涉对的小基线集组合,并对PS点的复数据进行干涉和差分干涉处理,得到去除平地和地形相位的PS点差分干涉图,经过空域二维相位解缠后得到解缠干涉相位,此时由影像获取时间为{ti,tj}组成的解缠干涉图,第m个PS点解缠相位由以下几个分量组成,即

(1)

(2)

式(2)中:hm表示该点的残余高程;B⊥为空间垂直基线;λ、r、θ对应雷达参数中的波长、斜距和入射角。形变相位为时间t的线性函数,可以根据小基线集组合相位时间序列,得到相位变化速率,即

(3)

式(3)中:φk为对应时间的形变相位。组合每个干涉相位观测值,可以得到线性方程组为

[Δφ1,Δφ2,…,ΔφK]T

(4)

式(4)中:T为系数矩阵,每行只有1和-1两个非0值,位置一一对应于干涉对组合;Δφ为PS点差分解缠相位。对式(4)利用奇异值分解SVD(singular value decomposition)算法,可以解算得到形变速率和残余高程值的最小二乘解,最后叠加上经过大气滤波后得到的非线性形变分量,得到最终的形变结果。

由于完整覆盖流域的Sentinel-1影像数据量较大,传统完全依赖人工的数据处理方法,需要耗费大量的人力和时间代价。为了应对InSAR大数据处理的挑战,根据以往的数据处理经验,利用适合流域的最优InSAR处理参数构建了全自动化PS-SBAS数据处理流程,在干涉处理流程,能够根据相干性阈值选取合适数量的干涉对,同时在时序分析过程,能够自适应地完成多次迭代形变参数解算,从而极大提升了数据处理效率。具体数据处理过程如下:以26-1区域数据为例,选取获取时间为2019-12-04的SAR影像作为主参考影像数据,将其他辅图像配准到主影像参考坐标下,采用针对TOPS模式数据的配准算法,得到配准好的SLC数据。基于幅度离差法选取PS点后,以时间基线100 d、空间基线150 m作为阈值获取干涉对组合,并利用30 m SRTM数据生成差分干涉图。待利用MCF(minimum cost flow)(最小费用流)算法进行相位解缠后,对形变和高程参数进行估计,并根据残余相位质量迭代多次解算,得到线性形变速率和残余高程,最后利用SVD算法进行解算,并经过大气滤波后,得到最终的形变时间序列,同时依据时序相干系数阈值0.7精选PS点,地理编码后得到shapefile格式的形变结果文件。完整的数据处理流程如图3所示。

图3 PS-SBAS处理流程图Fig.3 The flowchart of PS-SBAS method

3 结果分析

3.1 大渡河流域总体结果分析

采用第2节所述的PS-SBAS处理流程,对区域26-1、26-2和128-1这三个片区2018—2021年Sentinel-1数据进行了处理分析,获取了大渡河流域丹巴至乐山段完整的大范围InSAR形变测量结果,由于关注重点是大渡河流域以及S211省道沿线两侧的区域,截取了大渡河流域沿岸1 km缓冲区的形变结果,其时空分布如图4所示。

图4中每个点都对应于SAR影像上选取的高相干散射点,形变解算结果均为雷达视线方向上的(line of sight, LOS)位移速率(mm/year),并根据其年形变速率进行颜色编码,其正负代表是远离或靠近卫星方向。时序InSAR技术获取到总的点密度分布在106个/km2,总体形变速率范围为-65.2～51.1 mm/y。根据形变速率大小以及时序变化曲线,从总体上看,大部分区域地表形变都比较小或接近0,仅有少部分形变较大的区域主要聚集在丹巴、石棉和汉源等地。为了方便解译,在图4中对这几个区域结果进行了局部放大,从图4中可以观测到,沿着大渡河存在多处明显沉降区域(黄色和红色点)。根据历史滑坡成因分析的经验[18-20],流域不同区域滑坡诱发因素存在差异,其中处于大渡河上游的丹巴地区,地层较为稳定,区域滑坡最容易受到陡峭地形和河流侵蚀的影响;位于中下游的石棉和汉源地区,地质构造活动和降雨较易诱发滑坡失稳。因此,可以看到大渡河流域滑坡与地质和气候条件紧密相关。此外,还有许多没有PS测量点覆盖的区域,主要是由于植被覆盖茂密,以及坡体朝向、地形起伏与SAR侧视成像几何问题的影响,缺少足够多的PS测量点覆盖,无法获取准确的形变信息。

3.2 隐患区识别与分析

对大渡河流域InSAR结果大范围解译,综合考虑形变速率大小、累积变形量、形变区域面积和坡向这几个因素,同时要保证测量区域在SAR影像上相干性较高、测量结果较为可靠。通过结合光学图像对比分析,最终确定了15处威胁较大、形变较严重的隐患区。其中多个隐患区的变形区域发育在古滑坡体附近,坡体失稳严重,对边坡下方的居民区和大渡河构成较大威胁,需要地面测量人员重点关注,具体的隐患区分布和形变信息见表2所示。在隐患点列表中,不仅包含了以往研究发现的甲居村和丹巴县等受到古滑坡威胁的区域,以及地面重点观测的郑家坪边坡,还新增加了10多处大型滑坡隐患点,为开展地质调查提供了数据支持,获得了泸定地震前的灾害分布情况。因此,可以看到基于卫星InSAR技术和光学影像辅助解译方法,极大提高了大渡河流域灾害监测效率,对于识别一些隐蔽滑坡隐患区具有无可比拟的优势。

为了进一步说明隐患区域情况,同时验证InSAR测量结果的合理性,下面针对其中三个典型形变区域进行深入分析,包括001乡道滑坡隐患区、瀑布沟水电站和郑家坪边坡。

(1)001乡道滑坡隐患区。001乡道位于丹巴县,距离丹巴水电站大约2.5 km,该隐患边坡紧邻大渡河,坡度较缓,从光学图上可看到,部分区域进行了边坡防护。提取该区域InSAR测量结果如图5所示,可看到该不稳定区域覆盖面积达到11.2×104m2,区域最大形变速率达到了43.4 mm/y,从2018年1月—2021年6月最大累积沉降量超过了110 mm。提取坡面两个点TS1、TS2的形变时间曲线序列如图5(b)、图5(c)所示,这两个典型点位的形变曲线上表现出匀速下沉趋势,沉降主要来源于区域内道路施工以及工程建设活动,加剧了边坡失稳现象,InSAR测量结果与实际情况相符。整个坡面每年沉降量平均在35～45 mm,对坡底的居民区和学校构成直接的安全隐患,同时对大渡河带来一定威胁,其危险程度较高,需要地面测量人员持续关注。

(2)瀑布沟水电站沉降隐患区。在大渡河流域部分水库、坝体也表现出较明显沉降,瀑布沟水电站位于汉源县和甘洛县交界处,是该流域重要的梯级水电站之一。InSAR测量结果如图6所示,该水电站坝体最大形变速率达-22.966 mm/y,整个大坝每年沉降量在20 mm左右,目前属于正常沉降范围。可以观察到该结果与其他水坝形变具有类似的变化特征[21]。位于水坝受力最大的TS1处,其下沉速率最大,达到了23 mm/y,此外该处变形表现出明显的年际变化,与水库蓄放水过程具有强关联性,证明了InSAR结果的合理性;随着坝体位置变化,从TS1到TS2和TS3处,坝体沉降速率由-26 mm/y逐渐减缓至-5.8 mm/y,由匀速沉降变为波动式下沉,坝体周边区域较为稳定,需要对水电站坝体进行持续监测。

(3)郑家坪边坡测量结果。郑家坪滑坡体位于泸定县,是流域的重点关注区,有地面测量人员长期监视。从图7所示的该区域时序InSAR处理结果中,可以看到该坡面底部较为稳定,在坡顶处存在明显沉降区,最大沉降速率达到20.8 mm/y,区域覆盖面积大不,紧邻大渡河,2018—2021年累积沉降量超过75 mm,由于该地区活跃的地质构造运动,使得该边坡存在一定的失稳风险,对坡底的公路和大渡河造成较大威胁。

表2 主要滑坡隐患区列表Table 2 The identified hidden landslide region list

图5 001乡道附近形变分布和特征点形变时间变化曲线Fig.5 The deformation pattern in Road 001, and the displacement of TS1 and TS2

为了对InSAR形变测量结果进行验证,利用附近的地面GNSS测量结果(图7中▲所示)进行对比分析。获取了离InSAR测量区域最近的两个GNSS站点ZP13和ZP14,将其GNSS测量结果投影到InSAR的卫星视线LOS方向,为了去除两者参考点差异,对InSAR与GNSS两点分别差分的结果如图8所示。

从图8可以看到,InSAR测量结果与GNSS在整体能够保持较好的一致性,其形变趋势是相同的,两者形变量曲线间最大误差低于10 mm,主要差异是由于InSAR测量点与GNSS站点并不完全相同所致,误差在可接受范围,证明了该区域InSAR测量结果的准确性。

图8 郑家坪区域InSAR与GNSS测量结果对比图Fig.8 The deformation results comparison between InSAR and GNSS in Zhenjiaping

4 结论

基于时序InSAR自动化处理技术完成了对大渡河流域丹巴至乐山段大范围滑坡灾害监测普查研究,通过对2018—2021年期间大规模的Sentinel-1数据处理,首次获取了大渡河流域完整的滑坡隐患分布情况,构建了流域地质灾害快速监测识别方法。通过综合考虑区域形变情况,同时重点针对给水电站和人员带来较大安全威胁,以及堵塞河流风险较大的隐患点,在前人研究基础上,从测量结果中又提取了10多处新的隐患区域,最大沉降速率在-65～-30 mm/y,对这些形变风险区进行了详细描述和危害评估。最后,对几个典型隐患区域进行了深入分析和验证,进一步核实了InSAR测量结果的合理性。

综上所述,本文的研究证明了InSAR技术可以被用于复杂地质环境下的大范围地质灾害普查,构建的自动化InSAR处理流程可以形成常态化的形变监测,并以较低成本快速获取灾害分布一张图,持续为相关部门进行防灾减灾工作提供丰富的数据支撑,也为中国其他流域的灾害监测提供了参考范例,该测量结果将用于指导进一步的地面测量工作,同时对于后续对泸定地震前后流域形变变化情况分析也具有重要参考意义。