基于时序InSAR的宁夏西吉县滑坡灾害隐患识别

2022-11-26弓永峰吴学华

安全与环境工程 2022年6期

弓永峰,王辉,吴学华,张佳,2*,刘君

(1.宁夏回族自治区国土资源调查监测院，宁夏银川 750004；2.中国地质大学(武汉)环境学院，湖北武汉 430078；3.长安大学水利与环境学院，陕西西安 710054)

早期识别是防治地质灾害的重要技术手段之一[1]。对于人工调查难以到达的陡峻山区以及植被覆盖度较高的区域，滑坡灾害早期识别可以解决隐患点底数不清的问题，具有良好的地质灾害预警和风险管理作用[2]。近年来，合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)技术以其全天候、大范围、高精度的优势不断得到推广，对于滑坡灾害风险评估的发展起着重要的推进作用[3-5]。传统的时序InSAR技术通过差分合成孔径雷达干涉测量(Differential Interferometric Synthetic Aperture Radar，D-InSAR)、短基线合成孔径雷达干涉测量(Small Baseline Subset Interferometric Synthetic Aperture Radar，SBAS-InSAR)、永久散射体合成孔径雷达干涉测量(Permanent Scatterer Interferometric Synthetic Aperture Radar，PS-InSAR)等对滑坡灾害进行识别。其中，D-InSAR技术是对两景或三景影像数据进行处理，可快速获取地表形变结果，但无法得到长时序的地表形变信息；SBAS-InSAR技术适用于分布式目标分析，可得到大范围的面状地表形变信息；PS-InSAR适用于点目标分析，可提取可靠性高的PS点[6-7]。利用InSAR技术可以高效地进行区域地质灾害的早期识别[8]。

InSAR、光学遥感等综合遥感测量技术的应用对大区域潜在滑坡灾害的早期识别及滑坡体风险的评估具有重大的作用[9-12]。如赵宝强[13]通过InSAR技术圈定出白龙江流域部分潜在的滑坡隐患，并验证了InSAR技术的可靠性；王战卫[14]基于InSAR和光学遥感技术，识别出青海省大通县144个潜在的滑坡灾害隐患点，表明InSAR技术对大范围滑坡监测的可行性；黎光艳[15]利用InSAR技术对白龙江中游滑坡进行监测，分析了该地区滑坡分布和形变规律；赵超英等[16]利用SAR数据对甘肃黑方台滑坡形变进行了监测，研究了滑坡失稳模式；许强[17]提出构建地质灾害隐患早期识别的“三查”体系，通过光学遥感、InSAR、无人机航拍等技术在滑坡监测中的综合应用可取得较好的效果；长安大学张勤教授研究团队和成都理工大学许强教授研究团队结合InSAR技术和GNSS监测系统提前6 h成功预警了甘肃永靖县突发性黄土滑坡，此次成功预警为后续滑坡监测奠定了基础[18]。尽管InSAR技术在滑坡灾害识别和监测方面的应用越来越广泛，但由于时空失相干的影响，在地形起伏较大且植被覆盖度高的山区，滑坡灾害的识别难度依然很大[19]。因此，InSAR技术针对山区潜在滑坡灾害的识别效果和适用性还需要进一步深入分析。本研究的创新在于综合了SBAS-InSAR技术和PS-InSAR技术在提取地表形变信息方面的优势，针对宁夏南部山区潜在滑坡灾害进行早期识别，有效克服了由于地形地貌所带来的时间和空间上失相干的局限性，更好地消除了平地效应，得到了精确度更高的地表形变结果，表明了InSAR技术在宁夏南部山区潜在滑坡灾害高精度识别和监测的可行性。

宁夏南部山区属西北黄土高原的一部分，区域内山地与平原错落分布，高差起伏大，地表沟壑纵横，是宁夏地质灾害的高易发区。西吉县位于宁夏南部山区，地貌以黄土丘陵为主，次为红土丘陵，地质环境复杂，地质灾害呈易发多发态势。滑坡占西吉县地质灾害总数的86%，是西吉县最发育的地质灾害。本文以宁夏西吉县为研究区，首先利用PS-InSAR和SBAS-InSAR融合技术对研究区2018—2021年Sentinel-1A升轨数据进行干涉处理，得到了该区域高精度的地表形变结果；然后通过目视判读光学影像在地表形变结果中圈定出滑坡灾害隐患区域，并通过野外核查，检验滑坡灾害隐患识别结果的可靠性和有效性；最后分析了典型老滑坡形变时空特征，对典型焦湾滑坡群进行实地勘查，以为宁夏南部山区地质灾害隐患早期识别、监测预警以及治理搬迁避让等地质灾害防治工作提供充分有效的数据支撑和方法参考。

1 研究区概况

西吉县位于宁夏回族自治区南部山区(105°20′E～106°04′E，35°35′N～36°14′N)，县内交通较便利，东西长约67 km，南北宽约74 km，总面积为3 144 km2，见图1。西吉县境内发育清水河、葫芦河和祖历河三大季节性河流，根据当地地貌特征可将该地区分为黄土丘陵沟壑区、土石山区和河谷川道3个二级地貌单元，其中黄土丘陵沟壑区占总面积的83%，地势北高南低、东高西低，海拔在1 656～2 606 m之间。西吉县生态环境脆弱，年平均气温为12.7℃，年平均降水量为570.2 mm，降水主要集中在7～9月份。西吉县境内土地类型复杂多样，土层深厚，土质疏松，受1920年8.5级海原大地震的影响，县境内地震后有大量的新滑坡形成，同时因遭受地震的影响和破坏，该区域局部地区岩土体结构松散，加剧了区域内崩塌、滑坡地质灾害发生的可能性[20]。

图1 西吉县地理位置图

2 研究数据与方法

2.1 研究数据

本次研究所需的数据有：Sentinel-1A卫星影像数据、数字高程模型(Digital Elevation Model,DEM)数据、哨兵卫星精密轨道文件数据、光学遥感影像数据。其中，哨兵一号卫星是欧洲航天局哥白尼计划(GMES)分别于2014和2016年升空的两颗对地监测卫星，双卫星以12 d为周期对地球全天候昼夜雷达成像，幅宽为250 km，主要应用于山体滑坡、城市地面沉降监测等方面，哨兵卫星影像数据下载网址为https://search.asf.alaska.edu/，数据中心免费提供同轨道双卫星影像；DEM(SRTM)数据可提供参考地形，空间分辨率为90 m，其下载网址为https://srtm.csi.cgiar.org/；哨兵卫星精密轨道文件数据能够对轨道误差信息进行修正，其下载网址为http://www.gscloud.cn/；光学遥感影像数据为天地图影像，影像拍摄时间为2021年3月27日，空间分辨率为0.48 m。考虑到数据运行的时间成本，在下载数据时筛选出成像质量较好的数据。本研究选取2018年12月至2021年6月Sentinel-1A卫星影像数据共28景，数据类型为单视复数影像(Single Look Complex，SLC)，数据模式为干涉宽带(Interferometric Wide，IW)，轨道方向为升轨，空间分辨率为5×20 m、入射角为39.07°。Sentinel-1A卫星影像数据具体成像时间，见表1。

表1 Sentinel-1A卫星影像数据成像时间

2.2 时序InSAR方法

小基线集干涉测量(Small Baseline Subset-Interferometric Synthetic Aperture Radar，SBAS-InSAR)技术是Berardino等[6]在2002年提出的适用于分布式目标的InSAR时序分析方法，主要应用于地表形变时序监测领域的研究。其基本原理是依据同一地区时序上的N+1景SAR影像数据获得的时空基线设置相应的合理阈值，选出公共主影像，其他影像以主影像为参照组成无数个小基线集合并生成差分干涉图，并利用最小二乘法及矩阵的奇异值分解法求解地表形变时间序列及形变速率。

假设从研究区域获取影像时间排列为(t0，t1,…,tn)的N+1幅SAR影像，经过处理可以得到M幅差分干涉图，则M满足：

(1)

假设对不同时刻ta和tb获取的两幅影像进行干涉处理得到差分干涉图j(tb>ta)，则差分干涉图j对应某一像素点(方位向x，距离向r)的相位值可以表示为

δφj(x,r)=φ(tb,x,r)-φ(ta,x,r)

(2)

式中:φ为干涉相位;j∈(1，2，…，M)；λ为雷达波长；d(tb,x,r)和d(ta,x,r)为时间ta和tb以初始时间t0为对照[d(t0,x,r)=0]沿雷达视线(LOS)方向累计的地表形变量。

为了获取符合规律且更为直观的地表沉降序列，生成第j幅差分干涉图的相关时段与平均相位速度进行求积，可计算出该幅差分干涉图的相位值，则平均相位速度与相位分别表示如下：

(3)

(4)

式中：v为平均相位速度；φ为干涉相位；t为相位获取时间；δφj为第j幅差分干涉图的干涉相位；k∈(1,2,…,M)。

对各个时间段地表沉降速度进行积分，即可获取各个时段的地表沉降量，具体表示为矩阵积分形式：

Bv=δφ

(5)

式中当系数矩阵B(M×N)满秩或秩亏时，研究区域地表形变速率可分别利用最小二乘法和奇异值分解法求解，再根据地表形变速率求出相应的累积地表形变量。

永久散射体干涉测量(Persistent Scatterers-InSAR，PS-InSAR)技术是由Ferretti等[21]提出的针对具有高相干性的点目标进行差分干涉测量，其基本原理是选取同一地区时间连续的多景影像数据，通过设定适当的基线阈值从影像中选取1景为主影像，其余为从影像，并将所有影像与主影像配准后生成时间序列干涉相对，从研究区域中获取稳定的具有高相干性的永久散射体点，再使用DEM对干涉相对进行差分干涉处理，得到每个永久散射体点的差分干涉相位，其又包含相位的4个分量：

φdint=φdem+φdef+φatm+φnoise

(6)

式中:φdint为每个永久散射体点的差分干涉相位;φdem为DEM误差造成的地形相位;φdef为雷达视线方向(LOS)地表形变引起的形变相位;φatm为大气延迟引起的大气相位;φnoise为噪声引起的噪声相位。

通过对永久散射体进行时间序列分析，根据各相位分量特性将差分干涉相位中大气、噪声、DEM等相位分离出来，获取到每个永久散射体点的地表形变相位，从而准确地监测出研究区域地表位移变化值。通过该方法识别出SAR影像中后向散射特征稳定的点，可以避免人为选取过程中存在不稳定GCP点带来的误差。

2.3 山区地表形变监测

山区地形复杂，沟壑纵横，高差起伏较大，植被生长繁茂，导致传统调查监测手段无法获取区域性的地表损伤，从而影响对地质灾害危险性和风险性的判别。虽然PS技术局限于线性形变，但可以获取相干性高的GCP点；而SBAS技术在分析分布式目标方面具有很好的应用效果。因此，本文基于两种方法的优势互补，对PS-InSAR、SBAS-InSAR技术进行融合用于获取地表形变监测结果，相比分别利用PS-InSAR或SBAS-InSAR技术可有效减小SAR数据在处理过程中造成的误差，增加结果的准确性。基于SBAS-InSAR技术处理数据过程中，在轨道精炼时导入了利用PS-InSAR技术获取的高相干性的GCP点，以此来去除残余的相位信息，更好地消除了平地效应，提高了SBAS-InSAR技术处理结果的精确性。其具体操作步骤如下：

(1) 数据预处理：由于原始数据量较大，通过对研究区进行裁剪，可减少数据处理时间，提高处理效率。

(2) 利用PS-InSAR技术提取高相干性的GCP点：设置形变采样频率、残余高度采样频率、时空相关滤波、相干系数阈值等参数，减小大气相位误差，获取高质量、稳定的GCP点。

(3) 利用SBAS-InSAR技术进行差分干涉处理：通过设置一定的时空基线阈值，生成干涉相对，再进行配准、去平、滤波、相位解缠及相干计算，得到去平和滤波后的相位图及解缠图。

(4) 轨道精炼和重去平：导入PS-InSAR技术获取的GCP点，消除平地效应。

(5) 时间/空间域形变计算：通过两次SBAS反演有效去除大气相位，得到形变速率及形变量。

融合PS技术的SBAS-InSAR技术流程如图2所示。

图2 融合PS技术的SBAS-InSAR技术流程图

3 研究结果与分析

3.1 潜在滑坡灾害隐患识别

利用SARScape 5.2软件对研究数据进行处理，将结果导入Arcgis 10.3软件中进行密度分割，并添加光学遥感影像底图，最终得到西吉县雷达视线(LOS)向年平均地表形变速率在空间上的变化图。根据时序InSAR技术监测结果显示，相干点目标(CTs)为824 557个，平均相干点密度约为262 CTs/km2，能够满足与滑坡有关方面研究的要求[11]。根据斜坡的岩性特征、数据反演精度和实地调查结果，以及在前人研究基础上[12]，设置±10 mm/a为雷达视线向地表形变速率阈值，借助光学影像对潜在的滑坡进行识别，得到西吉县雷达视线向年平均地表形变速率图和识别的潜在滑坡灾害隐患点，如图3所示。图中地表形变速率正值说明地表运动相对卫星方向靠近；地表形变速率负值说明地表运动相对卫星方向远离。

图3 西吉县雷达视线(LOS)向年平均地表形变速率图和识别的潜在滑坡灾害隐患点

由图3可知：研究区内共识别出11个特大潜在的滑坡灾害隐患点，其中有3个历史崩塌，均属于斜坡在自然和人工条件下遭受破坏变形，工程地质学上统称其为不稳定斜坡类地质灾害[22-23]。大气降水的频次以及地壳活动能够在一定程度上诱发滑坡、泥石流等斜坡类地质灾害[24-26]。据调查，1996—1998年间研究区强降雨的频次高，发生滑坡、崩塌的点多，完全符合此时段内降水频次高且强度大的特点。西吉县境内地层以第四系黄土为主，由于黄土的垂直节理非常发育，斜坡边缘常发育平行于斜坡走向的张裂缝带，宽数米，加之黄土结构松散、孔隙率高以及黄土特有的湿陷性等，尤其在强降雨时，地表水入渗速度很快，土体含水量较大，其抗拉、抗剪强度大大降低，连续性降水导致地表水持续对斜坡的不断渗透，最终造成滑坡灾害发生。

3.2 典型潜在滑坡形变的时空特征分析

3.2.1 向家村老滑坡

向家村老滑坡位于西吉县城西南方向约26 km处，其光学影像图和InSAR结果图如图4所示。图中红色界线区域为滑坡范围；黄色界线区域代表滑坡早期识别加剧范围。

图4 向家村老滑坡光学影像图和InSAR结果图

由图4可以看出，向家村老滑坡滑坡体呈不规则长舌状，周界清晰可见，后缘发生明显形变，滑坡体中下部有耕地分布，滑坡下部有道路横穿。

图5为向家村老滑坡的地表形变速率分布图。

图5 向家村老滑坡的地表形变速率分布图

由图5可知，向家村老滑坡最大地表形变速率位于滑坡体形变区后缘，其值达到-17 mm/a，滑坡范围超出原有滑坡边界。

为了对向家村老滑坡滑坡体的形变趋势进行深入研究，选取1#～5#5个样点(具体位置如图5标记所示)形变进行了时序分析，得到向家村老滑坡形变的时间序列图，见图6。其中，1#样点位于该滑坡体后缘，该处发生明显形变，累计地表形变量达-50 mm；2#样点位于该滑坡体后侧边缘，自2018年12月起该样点形变持续呈下降趋势，累计地表形变量达-30 mm；3#样点位于滑坡体中间位置，该样点形变趋于加速状态，累计地表形变量为10.21 mm；4#样点位于该滑坡治理区左侧，监测期间该样点累计地表形变量仅为-2 mm；5#样点位于该滑坡治理区右侧，该样点形变在监测期间内基本趋于稳定状态。这是因为：1#样点、2#样点位于该滑坡后缘，由于连续降雨造成松散物源下移，导致滑坡后缘形变明显；3#样点处缓慢抬升，主要是由于该滑坡后缘部分松散岩土下滑在该处堆积所致；4#样点、5#样点位于滑坡治理区附近，由于人为活动的影响，使该处地表形变趋于稳定状态。通过有效整治后该滑坡体前缘处地表形变基本保持稳定状态，而滑坡体中后缘地表形变趋势仍呈加速状态且累计地表形变量已较大，一旦平衡被打破，滑坡体下方居民的生命和财产安全将无法得到保证，且还可能会影响道路交通安全。

图6 向家村老滑坡的形变时间序列图

3.2.2 焦湾滑坡

焦湾滑坡为大型滑坡群，位于西吉县平峰镇焦湾村焦湾组，该处滑坡主要由1920年、1970年和2008年发生在滥泥河左岸黄土斜坡上的几次大地震诱发所致，其滑坡体较陡，上覆黄土，土体破碎，落水洞极其发育，滑坡前缘面临冲沟沟底有少量的新近系红泥岩出露。该地区为典型的黄土丘陵沟壑区，黄土梁整体走向西北，地势由东南向西北倾斜，峁梁起伏，沟壑纵横，滑坡发育于黄土丘陵斜坡上。根据InSAR结果显示，该大型滑坡群在监测时间内又产生了新的形变，具有老滑坡复活的可能性。结合光学遥感解译和InSAR技术对焦湾滑坡灾害隐患点进行识别，选取焦湾滑坡群中有形变的区域进行分析，其结果如图7所示。图中红色边界内代表滑坡范围。通过对InSAR得到的量化结果进行分级并用不同的颜色显示，图7(b)中蓝色部分表示滑坡形变区域，6#样点处地表形变速率为16 mm/a，7#样点处地表形变速率为13 mm/a。7#样点位于焦湾滑坡边界外部，随时有滑塌的危险。

图7 焦湾滑坡光学影像图和InSAR结果图

为了更加直观深入地了解焦湾滑坡群形变情况及滑坡成因，对该滑坡展开了野外实地勘查。野外实地勘查结果显示：该滑坡位于焦湾村西北部，地理位置为东经105°33′34″、北纬35°47′53″，滑坡全貌如图8(a)所示；该滑坡周界清晰，呈圆弧状，后壁较陡，高程为2 060 m，高约25 m，滑坡体较陡，坡度近65°，滑坡体上土体较为破碎，现已被人工改造为耕地[见图8(b)]；滑坡体遭流水侵蚀又出现新的侵蚀沟，使滑坡体稳定性下降，坡壁发育有落水洞[见图8(c)]；该滑坡周边路面塌陷现象严重[见图8(d)]；该滑坡滑体土主要为粉土，颜色为黄褐—黑褐色，均质结构，土层松散，抗剪强度低，垂直节理发育，含粉砂、暗色矿物，发育虫孔、根孔，岩土体含水率较低，处于可塑状态；该滑坡滑床主要岩性为晚更新世黄土，黑褐色，岩土体含水率较高，垂直节理发育。焦湾滑坡属于以地震为主要触发因素的滑坡，该滑坡前缘受河流冲刷作用，在雨季易发生蠕滑；滑坡区地层主要为更新统风积黄土和第四系滑坡堆积物，岩土体裂隙发育，物理力学性质较差，土体较破碎，降水易入渗，可能会引起滑坡整体失稳[27]。因此，对焦湾滑坡的防治应以监测为主，如有条件可在该滑坡体前缘辅以防护工程，以减少河流的冲蚀作用，一旦发现有形变迹象，立即采取防灾措施。

图8 焦湾滑坡的野外实地勘查结果

经分析，焦湾滑坡滑坡体裂缝、落水洞发育，应该及时填压，以防止雨水入渗；裸露的坡面易于地表水的下渗，应合理植树造林，改善坡面植被条件；滑坡体前缘存在临空面，应尽量避免可能扰动滑坡的各种人类工程活动，并严禁开挖坡脚；合理利用滑坡土地，不宜在滑坡上兴建有利于降水入渗的各种水利工程；建立滑坡灾害监测网络，加强雨季老滑坡动态的监测，发现形变迹象应立即上报主管部门，以便制定防灾应急对策。