基于小基线集雷达干涉测量的中巴公路盖孜河谷地质灾害早期识别

赵富萌，张 毅，孟兴民,3，苏晓军，石 伟

(1.兰州大学资源环境学院，甘肃 兰州 730000；2.甘肃省环境地质与灾害防治研究中心，甘肃 兰州 730000；3.兰州大学地质科学与矿产资源学院，甘肃 兰州 730000)

中巴喀喇昆仑公路(Karakoran Highway，KKH)北起中国新疆的喀什，穿过红旗拉普山口，南到巴基斯坦北部城市塔科特，全线长约1 036 km。公路位于帕米尔东北缘—西昆仑地区，属于帕米尔构造前缘地带，是中国大陆受板块动力作用和地震活动最强烈的地区之一[1]。在地震、降水和冰川活动下，公路两侧极易发生滑坡、泥石流、冰崩和雪崩等灾害，因此对公路沿线的灾害进行早期识别研究对公路正常通车和中巴经济走廊建设具有重要意义。

从20世纪50年代开始一直有学者对中巴公路沿线的地质构造、冰川活动及滑坡、泥石流发育等方面进行研究，但由于研究区的地质灾害具有灾种多、规模大、高频率分布的特点，使得研究多集中于对典型灾区的定性描述。自20世纪80年代开始对喀什至塔什库尔干路段冰川泥石流形成原因进行调查，得出区域泥石流暴发频繁的原因在于冰雪消融洪水和岩崩、沟岸坍塌等形成的松散堆积物为泥石流形成提供了水源和物源条件[2]；后来的学者又对该区域的地质条件、岩崩沉积对河流阻塞的影响、冰川灾害以及其他一些工程地质问题进行了讨论[3-6]。但由于公路地处高寒无人区，地质环境复杂，难以进行全面细致的调查，尤其是对区域地表变形监测和地质灾害早期识别缺乏深入研究。

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)技术由于不受天气制约和大范围、高精度的监测优势目前已被广泛应用于地质灾害的识别研究中[7-8]，尤其小基线集合成孔径雷达干涉测量(Small Baseline Subsets InSAR, SBAS-InSAR)技术更是克服了传统干涉测量方法中存在的时间、空间失相干以及大气效应的限制，使得形变结果在时间和空间上更为连续，适用于监测长时间序列的缓慢非线性形变，目前已被广泛应用于地表沉降[9]、地震[10]、火山[11]、活动断层[12]以及冻土形变和冰川运移[13-14]等领域，监测优势明显且潜力巨大，在识别不稳定斜坡如蠕移、滑坡[15]等也有良好结果。

本文以中巴公路地质灾害最为发育的盖孜河谷段10 km缓冲区为研究对象，利用2016—2017年共28景升轨Sentinel-1A卫星影像通过裁剪、干涉处理、轨道精炼、反演等得到研究区域雷达视线方向(Line of Sight, LOS)向形变速率图，并通过公式转换为斜坡向形变速率图，结合野外验证识别出区域不稳定斜坡，并结合坡度和降水数据对典型的地质灾害点进行分析，从而实现对该区域地质灾害的早期识别预警。

1 研究区概况

中巴喀喇昆仑公路(中国段)北起中国喀什的疏附县，沿G314经奥依塔克镇穿盖孜河谷到布伦口盆地，经苏巴士达坂和塔什库尔干到达中巴边境口岸红旗拉普山口。公路沿线地质环境复杂，构造活动强烈，地震频发，由地震引发的崩塌、泥石流等地质灾害对公路沿线人民的生命财产安全造成巨大威胁[16]。区域剧烈的构造活动及地震活动破坏岩土体结构，降低岩土体力学强度，使岩土体变形破碎，为地表渗流和地下水活动提供了良好通道，为崩塌、泥石流的形成、发育和发生提供了物质基础和水源条件。

盖孜河谷横切西昆仑山脉，是中巴公路的咽喉地段，河谷北依昆盖山，最高海拔6 670 m；南靠公格尔山，最高海拔7 649 m。河谷平均海拔约2 500 m，最宽处不到2 000 m，最窄处不足5 m。深狭的盖孜河谷两侧羽状分布着众多泥石流沟，沟源区至沟口的相对高差多在3 000 m以上，以公格尔山北侧的克拉牙依拉克冰川沟谷规模最大。河谷地区断层广泛发育，昆盖山南麓正断层、盖孜断层、主帕米尔断层和肯别尔特断层这四条断层在河谷走向均为NW—SE向，区内地层以石炭系和第四系为主，岩性包括石炭系的花岗岩、砂岩和灰岩，泥盆和二叠系的石英砂岩和片岩、大理岩等，三叠、侏罗和白垩系的红层泥岩和石膏层，新近系的泥岩和砂岩以及第四系复合成因的沉积物如砂土、黄土和砾石等(图1)。区域由于地处欧亚大陆腹地，典型的大陆性气候使得区域降水稀少，蒸发强烈，寒旱的气候特征和强烈的冰川及寒冻风化作用使得区域植被稀疏，山体裸露，仅在盖孜河的河道两侧生长有少量天然荒漠植被。盖孜河的支流包括奥依塔克河、木吉河、康西瓦尔河，径流补给以冰雪融水为主。通过盖孜河中游克勒克水文站1958—1996年的实测资料统计，盖孜河多年平均流量为31.3 m3/s，河水年内四季水量分配比例为春夏秋冬各占14.88%、60.75%、18.11%和6.26%，年内丰枯变化悬殊[17]。克勒克水文站点资料显示区域多年平均气温为8.18 ℃，最高月均温为20.7 ℃，出现在7月；多年平均降水量为123 mm，集中于每年的4—9月，占比达到76.3%；站点多年平均蒸发量约2 355 mm。区域地下水主要为第四系松散层孔隙水和基岩裂隙水，冰川融雪水是地下水的主要补给来源。受构造运动、河流侵蚀和冰川运动的影响，盖孜河谷地势险峻，地质灾害频发，滑坡、崩塌等斜坡地质灾害以及冰川区的冰川跃动和冰川泥石流等对公路安全造成巨大威胁，尤其在公格尔山北侧路段冰川泥石流风险性高，溜石坡也十分发育[18]。

图1 研究区地质图概况(源于中国地质调查局1∶50万地质图数据库)Fig.1 Survey of geological maps in the study area (from 1∶500 000 geological map database of China Geological Survey)

2 方法与数据

2.1 SBAS-InSAR技术原理及方法

InSAR相比于传统的形变监测方法如GPS等具有不受天气制约和大面积监测等优势而被广泛应用于地质灾害识别中[19]。而时序干涉测量技术则利用多期影像对地表形变进行反演，从而可以对区域地表形变的时间演化有很好的分析能力[20]。其中SBAS-InSAR的方法为最小化视角差异造成的误差和失相干，选取时间和空间基线都较小的SAR影像干涉对，对SAR影像干涉对进行差分多视处理，对差分干涉图的像元基于相干性筛选出高相干点，然后对这些高相干像元进行相位解缠，利用地面控制点进行轨道精炼，去除相干性较低和解缠错误的干涉像对，通过奇异值分解法和最小二乘法解得影像获取时刻的累计形变，最后利用时域的高通滤波和空间域的低通滤波去除大气影响得到最终的形变结果[9,21]。该方法在山区地质灾害早期识别中有较广泛应用。

SBAS方法获得沿雷达视线方向(LOS)的地表形变信息，而对于山区来说， LOS向的形变速率不足以真实反映斜坡的真正形变情况。通过公式将LOS向的形变速率转换为沿坡度方向的形变速率[22]：

Index=nlos·nslope

nlos=(-sinθcosαS,sinθsinαS,cosθ)

nslope=(-sinαcosφ,-cosαcosφ,sinφ)

式中：Vslope——沿坡度方向的形变速率；

Vlos——LOS向的形变速率；

α——斜坡的坡向/(°)；

φ——斜坡的坡度/(°)；

θ——雷达入射角/(°)；

αS——卫星轨道方向与正北方向的夹角，即雷达卫星飞行方向，升轨数据为负，降轨数据为正。

2.2 实验数据与处理

本研究选择中巴公路灾害发生频率最高的盖孜河谷段为研究对象，以公路线为中心分别向两侧延伸10 km的缓冲区作为实际研究范围，利用SBAS-InSAR方法得到研究区地表形变速率结果，数据选择覆盖研究范围从2016—2017年共28景的升轨Sentinel-1A卫星影像(表1)，通过设置最大空间基线阈值为110 m和最大时间基线阈值为90 d共生成162个干涉像对(图2)，借助美国地质调查局的SRTM(Shuttle RadarTopography Mission)模拟地形相位，在ENVI软件下的SARScape模块支持下对SAR影像进行裁剪、连接图生成、影像配准、差分干涉、自适应滤波、相位解缠、轨道精炼、形变相位估计、大气效应去除等处理，得到最终每个相干点的形变速率值和时间序列的非线性形变图。

图2 SBAS处理中的干涉对时间-空间基线连接图Fig.2 Time-space baseline connection diagram of the SBAS processing interference pair

表1 Sentinel-1A数据属性表

为分析时间序列的形变结果与降水之间的关系，需要获取研究区降水数据，但由于区域气象站点缺乏，我们选择了携带热带降雨测量任务(Tropical Rainfall Measuring Mission，TRMM)的遥感卫星的3B43数据，其空间分辨率为0.25°× 0.25°，基本满足研究区对降水数据空间分辨率的需求。TRMM数据在广泛应用的同时，也有许多学者将其与台站、雷达观测的其他降水数据进行了比较，对其应用于实际分析中的精度进行了验证。如最开始采用TRMM 3B41 RT数据在青藏高原周边区域与观测数据进行比对，证实了TRMM数据大致上可以反映降水的基本特征[23]，后又将TRMM 3B42 RT数据在日和月尺度上与气象观测数据进行一致性检测，表明TRMM降水数据与观测数据之间具有较高的线性相关特性[24]。而本研究使用的TRMM 3B43数据，经验证在月尺度上的局部地区的相关系数接近0.9[25]。

3 结果与分析

3.1 SBAS结果

SBAS-InSAR方法干涉处理得到区域LOS向形变值为-76～28 mm/a，依据形变速率值标准差选取-5～5 mm/a作为形变相对稳定值，负值表示形变点沿LOS向远离传感器的方向运动；正值表示形变点沿LOS向接近传感器的方向运动(图3)。由于采用的升轨数据，正值LOS向的形变在沿斜坡向的分量并不一定也是正值，另外雷达视线向形变结果受山区地形效应影响，与实际形变特征出现偏差，所以根据转换关系将LOS向形变速率值转为坡度向形变速率值。当公式中的Index趋于0，Vslope趋于无穷大。为剔除这一异常值，Herrera等以经验值Index=±0.3作为固定阈值，设定当-0.3

图3 盖孜河谷段2016—2017年LOS向形变速率Fig.3 Map showing the LOS deformation rate of the Gaizi valley from 2016 to 2017

图4 盖孜河谷段2016—2017年斜坡向形变速率Fig.4 Map showing the slope deformation rate of the Gaizi valley from 2016 to 2017

3.2 区域地质灾害变形分布特征

根据2016—2017年斜坡向滑移速率，结合前人相关研究成果和经验，选取数据标准差(σ=15 mm/a)为稳定点阈值，利用斜坡向滑移速率和野外验证圈定出449处地质灾害点，包括31处滑坡、416处不稳定斜坡和2处冰川运动(图5)。该区滑坡主要是表现在短时强降水条件下斜坡上部不稳定土体整体向下滑移的过程；不稳定斜坡空间分布频率最高、变形速率最大的集中于公路两侧的斜坡上，尤其是在公路北侧的东向斜坡存在大量不稳定点，在公路南侧冰川发育地区由冰碛物运动引起斜坡失稳；在公路南段公格尔山处冰川运动活跃，冰碛物广泛存在，地质环境很不稳定，斜坡向滑移量显示存在多处由冰碛物运动引起的滑坡和不稳定斜坡。由于InSAR结果不能识别泥石流，而盖孜河谷段由于降雨和冰川融雪引发的泥石流灾害尤为严重，所以本研究基于Google影像解译和现场调查圈定出23条泥石流沟，分布于公路两侧(图5)，在每年6—9月爆发的泥石流对公路正常通车造成极大阻碍。

图5 盖孜河谷段地质灾害类型分布图Fig.5 Map showing the distribution of geological hazard type in the Gaizi valley

3.2.1滑坡变形特征

研究区滑坡崩塌灾害的诱因包括两种：地震滑坡崩塌和降雨融雪型滑坡崩塌。在盖孜河谷段中部有多条断层穿过，地震活动频繁，再加上强降水的诱发，使得滑坡崩塌灾害严重[27]，以盖孜河谷段某处滑坡为例，分析其变形特征(图6)：该滑坡位于艾尔库然沟右侧斜坡上，地层上属于泥盆系的克孜勒陶组，上部为黄褐、灰黑色泥质粉砂岩、粉砂质泥岩夹生物碎屑灰岩；下部为碎屑岩，在降水作用下，极易发生失稳滑动。图6(a)中A1点斜坡向滑移速率为-33 mm/a，中间因形变速率超过InSAR监测能力而出现空值；图6(b)中的月降水数据来自于美国国家航空航天局(NASA)的0.25°空间分辨率的TRMM 3B43卫星数据，地面空间分辨率约2.75 km，基本能代表斜坡上的降水强度，另外TRMM数据能够保证不同海拔高度降水数据的均一性，所以基本满足地基台站降水数据极度缺乏的盖孜河谷段的降水分析[26]。图6(b)中A1点形变曲线显示在2017年1月和2月两月间的累计形变量为-32 mm，对应的两个月累计降水量为143 mm，强降水对应于强变形，但受水力作用于土壤的过程，所以斜坡滑移会稍滞后于强降水事件的发生，滑坡变形对降水响应较为敏感。图6 (c)野外验证中出现多条因降雨冲蚀出现的细沟，极易发生滑塌。

图6 滑坡识别分布特征图Fig.6 Map showing the landslide identification distribution

3.2.2不稳定斜坡变形特征

图7 不稳定斜坡数量与地层岩性和坡度的关系Fig.7 Relationship between the numbers of unstable slopes and stratigraphic lithology and slope

中巴公路横穿帕米尔西构造结，破碎的构造带、活跃的构造活动、短时间的强降暴雨和盖孜河对斜坡坡脚不断的深切侵蚀是斜坡运动的主要控制因素和诱发因子。通过斜坡向滑移速率识别出区域不稳定斜坡416处，占所有识别的地质灾害类型的85.6%。图5显示不稳定斜坡集中于坡度较大、植被稀疏的高海拔山区地带。区域地质条件复杂，寒冻风化作用强烈，使得崩塌滑落发育。影响斜坡灾害发生的因素主要包括地层岩性和坡度两个方面，地层岩性条件表现为不稳定斜坡点的51%分布于石炭、二叠系的粉细砂岩、片岩和石英片岩区域，18%分布于泥盆系泥质粉砂岩区域，在三叠系的砂岩、粉砂岩，白垩系的泥岩，新近系的泥岩、砂岩以及全新世的复合成因堆积物也分布有规模不一的不稳定斜坡(图7a)。在区域强烈的冻融风化作用下，岩石内部的强度降低，节理裂隙发育，一方面气温升高时裂隙内的积雪融化会降低岩体的抗剪强度，使得岩体极易破裂而发生崩塌滑落；另一方面冰雪融水沿地裂缝渗入岩体，聚集于相对隔水层之上，软化侵蚀岩体结构，使岩体失稳而发生滑落。统计识别出的416处不稳定斜坡发育与坡度值对应关系如图7(b)所示，可以看出不稳定斜坡集中发育于坡度在20°～50°之间，占总的不稳定斜坡数量的75%，不稳定斜坡发育数量与坡度并没有明显的正相关关系，不稳定斜坡的发育受多种因素控制，除岩性和斜坡坡度外，区域断层活动、地震、降水和植被盖度对斜坡稳定性也会产生影响。

图8 不稳定斜坡识别分布特征图Fig.8 Map showing the unstable slope identification distribution

以盖孜检查站附近的斜坡灾害为例分析其发生与降水的关系(图8)。该斜坡海拔约4 300 m，由斜坡向形变速率图可以看出，斜坡处于极不稳定状态，形变剧烈且分布广泛(图8 a)。在斜坡下部附近出现空白区域是因为实际形变值超出InSAR监测能力而出现空值。图8 (a)中点A1和A2的时间序列形变图见图8(b)，A1点斜坡向滑移量从2016年1月6日到2017年12月14日累计达-81 mm，且在2017年1月的位移量超过-11 mm，通过TRMM降水数据可以看到2017年1月的降水量为64 mm，占全年降水量的18%。这一降水过程也是全球气候异常的一个表现，因强降水事件的发生而导致斜坡变形，所以对区域强降水的发生需给以重视，预防灾害的发生；A2点从2016年1月6日到2017年12月14日的累积位移量为-54 mm。图8(c)、(d)是发生斜坡坡面碎屑流的野外照片，坡度较陡，为上部碎屑岩体滑落提供了动力条件，表面几乎无植被覆盖，碎屑物在雨水作用下沿斜坡冲刷出多条细沟，在坡脚处堆积。盖孜河谷属于干热河谷，降水稀少且集中，每年的降水量以历时短暂的暴雨形式出现。降水是引发斜坡失稳的重要原因，主要通过以下方式影响斜坡的稳定性：一方面降水入渗使土体达到饱和，孔隙水压力增加，斜坡的抗剪强度降低，易失稳下滑；另一方面降水入渗到岩土体中，增加土体重量，从而降低软弱结构面的抗剪强度，也会使滑带土得到软化，极易发生土体滑动，同时也说明斜坡变形会稍滞后于强降水事件。另外地震也是斜坡失稳的重要原因。地震导致坡体结构破坏，岩土体疏松，强降雨沿着疏松滑坡坡体裂缝、落水洞等优势通道快速入渗，从而导致岩土体含水量上升、强度降低，地下水位也随之上升，诱发坡体变形[28]。

3.2.3冰川运动变形特征

图9 冰川运动点分布特征图Fig.9 Features of the glacier movement identification distribution

公路南段的公格尔山平均海拔超过4 000 m，冰川活动强烈，尤其是公格尔山，最高峰海拔7 649 m，是冰川源区，由于变形相对较大超出InSAR识别的范围导致严重的失相干。图9显示了公格尔山北侧冰川携带冰碛物运动情况。在InSAR识别的初始阶段，冰川运动强烈，A1点的冰碛物沿斜坡向下游接收区输送冰体和碎屑物质，大量物质堆积造成该处冰面隆起，变现为曲线的上升(图9 b)，当部分碎屑流继续向下运动， A1点失去源区冰碛物补给而不再发生正的变形，但由于该地区的土质多为季节性冻土，极不稳定，在冲击作用下使得下部土体下沉从而表现为负的变形，如此循环，对斜坡土体造成更大破坏，加剧了不稳定斜坡的发生，甚至引发冰碛物滑坡。

3.2.4泥石流沟运动特征

中巴公路盖孜河谷段两侧泥石流广泛发育。据调查发现区域泥石流按照形态可分为沟谷型泥石流和坡面型泥石流。沟谷型泥石流多发育于公路北侧的山谷内，流通区和发育区的分区不明显，多分布在山谷和山前坡地；沟谷两侧的坡脚处发育坡积物，沟床内的堆积物厚度较大，松散固体物质较为丰富。坡面型泥石流的形成区一般分布在公路两侧的山体上。山体较为破碎，堆积物较多，颗粒松散，粒径较大，新发生的泥石流堆积于老泥石流扇上，坡面水流汇集，表层松散堆积物直接转化为泥石流。

区域泥石流按照水源条件可分为融雪型泥石流和洪水型泥石流两类。融雪型泥石流集中于公路两侧的昆盖山和公格尔山冰川分布区，一般为沟谷型泥石流，海拔多在4 000 m以上。冰川泥石流沟由于物源区被冰川积雪所覆盖，状态极不稳定：一方面夏季的持续高温使冰川积雪急剧消融，突发水源大量增加，另一方面冰川融水在向下流动时会携带大量的冰碛或雪崩岩屑等固体物质。洪水型泥石流分布广泛，在年降水量极少的干旱期由于强烈的物理风化作用使得地表形成丰富的松散碎屑物质，在公路北侧普遍存在着有利于汇水的围谷状地貌，且形成区面积大、坡面多、坡度陡、沟壑密度大，暴雨作用下造成泥石流发生[29]。

4 结论

(1)SBAS的监测结果显示2016—2017年研究区LOS向形变速率为-76～28 mm/a；结合研究区的坡度、坡向以及卫星采集数据的几何姿态等信息将LOS向形变速率转换为斜坡向形变速率最大达到-157 mm/a；基于斜坡向滑移速率，得到覆盖研究区的449处灾害点，包括31处滑坡，416处不稳定斜坡和2处冰川运动，形变点集中分布于道路北侧；不稳定斜坡集中的沟道为泥石流发育提供了充足的物源物质，使得盖孜河谷段的泥石流沟广泛存在。遥感解译结合野外调查数据统计出研究区存在23个泥石流沟。

(2)结合区域TRMM降水数据对典型地质灾害点的时间序列变形结果进行分析，证实了斜坡运动与强降水存在较强相关性，但因水对岩体的入渗侵蚀需要一个过程，所以斜坡滑移滞后于强降水事件的发生。在当前全球气候变化背景下，会发生越来越多气候异常事件，这也迫切需要对地质灾害隐患点进行监测和早期识别。

(3)该研究证明了SBAS-InSAR方法在地质灾害早期识别和监测中的应用优势，但受限于方法和数据的不足，对植被覆盖区和冰川区识别结果较少，受卫星重访周期影响，对短时间发生的崩塌等地质灾害没有有效识别，在今后研究中应合理选择雷达数据和处理方法，以有效快速地对地质灾害进行早期识别。