本文研究区域位于西藏自治区贡觉县雄松乡（98.92°E，30.49°N）至沙东乡（98.9°E，30.61°N）的金沙江流域（见图 1），从谷歌影像目视解译发现可能存在滑坡。本文选取南北约15 km、东西约6 km的矩形区域进行滑坡探测实验，黄色矩形为研究区域，蓝色矩形为哨兵一号（Sentinel-1）升轨影像覆盖范围，绿色矩形为Sentinel-1降轨影像覆盖范围，红色矩形为陆地观测技术卫星2号（advanced land observing satellite 2，ALOS-2）升轨影像覆盖范围。研究区域海拔最高为4 608 m，最低为2 594 m，地形落差大，地形复杂，在河流附近极易孕育出滑坡群。

图  1  研究区域与SAR卫星影像覆盖范围

Figure 1.  Study Region with the Coverages of SAR Imagery

Sentine-l系列卫星是欧空局（European Space Agency，ESA）由哥白尼方案资助的对地观测卫星。Sentinel-1A卫星和Sentinel-1B卫星运行在同一轨道平面上，两颗卫星星座昼夜工作，并通过哥白尼访问中心提供免费、开放、连续的产品。该卫星工作在C波段，有条带模式（strip mapmode，SM）、干涉宽幅模式（interferometric wideswath，IW）、超宽幅模式（extra wideswath，EW）、波模式（wave mode，WV）4种工作模式。本文选择了IW模式下垂直发射垂直接收的同向极化方式（vertical vertical，VV）得到的单视复数影像（single look complex，SLC）。无论天气如何，Sentinel-1单颗卫星采用渐进式扫描地形观测模式（terrain observation by progressive scans，TOPS），每12 d就可以获得全球范围内的图像^[20]。本文选取2017-10-08—2019-10-10共61景Sentinel-1升轨影像和2017-07-23—2019-07-25共53景Sentinel-1降轨影像用于分析。

ALOS-2卫星携带L波段合成孔径雷达（synthetic aperture radar，SAR）传感器，穿透力强，具有聚光模式（Spotlight）、SM模式、扫描模式（ScanSAR）3种观测模式。由于观察模式和研究区域的不同，重访时间也从几个月到两周不等，与C波段相比，L波段由于雷达波长较长，能够更好地穿透植被获取地面信息，便于在高海拔、高植被覆盖地区保持较好的相干性，因此，本文选取SM3模式下水平发射水平接收的同向极化方式（horizontal horizontal，HH）得到的SLC，2017-07-24—2019-07-22 7景ALOS-2升轨影像用于研究。

本文采用Gamma处理1级单视复数格式的SAR影像^[21]，为了抑制斑点噪声，Sentinel-1影像在距离向和方位向采用10×2的多视因子。ALOS-2影像在距离向和方位向采用2×5的多视因子，与其他传统的SAR影像数据处理不同，Sentinel-1在TOPS模式下获取数据时，由于天线会发生旋转，需要更严格的图像配准，配准的方位角精度应优于0.001像素，以避免相邻Burst之间的相位跳跃^[22]。

利用30 m分辨率的航天飞机雷达地形测绘使命（shuttle radar topography mission，SRTM）数字高程模型（digital elevation model，DEM）（https://earthexplorer.usgs.gov）削弱地形的影响，Sentinel-1影像设置的时间基线为90 d，空间基线为150 m，由于ALOS-2数据偏少，因此本文对ALOS-2数据进行全干涉，得到ALOS-2数据干涉对组合21对（见表 1），Sentinel-1升轨影像干涉对组合212对，Sentinel-1降轨影像干涉对组合188对。

GACOS系统基于高分辨率欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts，ECMWF）$\mathrm{0.1}°\times \mathrm{0.1}°$格网的大气数值产品，采用迭代对流层分解（iterative tropospheric decomposition，ITD）模型从对流层延迟中分离出分层和湍流分量，分别进行空间内插，得到高空间分辨率（默认为90 m）的天顶对流层延迟（zenith troposphere delay，ZTD）改正影像^[23]。

InSAR Stacking的基本假设是研究区域内的地表形变是线性的，不存在非线性的变化，其优点是运算高效，可有效削弱随机误差的影响。引入GACOS的ZTD产品削弱了大气误差，残余的大气误差可视为随机误差，通过InSAR Stacking技术可提高信噪比^[24]。本文利用最小二乘对解缠相位求平均，干涉相位按时间间隔加权，得到研究区域LOS向年形变速率为^[25]：

$\mathrm{式}\mathrm{中},V_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$为平均相位速率；${\phi }_i$为解缠后的相位；$\mathrm{\Delta }T$为对于干涉图的时间基线；N为干涉图的个数。

LiCSBAS是一个与LiCSAR系统集成的开放源InSAR时间序列分析软件包^[16]，本文利用LiCSBAS解算InSAR雷达LOS位移时间序列。尽管Sentinel‐1影像的空间和时间基线较短，但由于研究区域的高植被、高积雪覆盖以及地形起伏出现严重的失相关和解缠误差，在进行时序分析之前需要对解缠后的干涉图进行质量检查，修改或删除^[26]。解算步骤如下：首先利用GACOS对流层延迟产品对干涉对进行有效的改正，剔除解缠误差，对整体干涉对进行质检。以SAR影像的有效像元率（每对干涉图上有效像元数/有效像元总数）和平均相干性作为阈值，对不满足条件的干涉对进行剔除，Sentinel-1升轨影像干涉对中22个干涉对被剔除，Sentinel-1降轨影像干涉对中10个干涉对被剔除。然后对干涉图进行相位环闭合的检验^[27]，逐一计算干涉对环闭合相位均方根（root mean square，RMS），为了防止RMS满足阈值的干涉对中仍然存在解缠误差的像元，该方法统计了每个像元环闭合相位误差大于阈值（1.5 rad）的个数，作为掩模像元的标准。最后在满足条件的情况下，得到研究区域的LOS向的年形变速率。

坡向是指坡面法线在水平面上投影的方向，坡度代表地表单元陡缓的程度，通常把坡面垂直高度与水平方向距离的比叫做坡度。坡度和坡向大致能反映一个地形的特征要素，往往滑坡的发生除了必要的外部环境条件，其位移特征与坡度和坡向有着很大的关系。

本文基于滑坡的运动规律，以滑坡的滑动面作为观测其运动的参考面，引入坡度和坡向两个地形因子作为滑动面坐标系的要素。在现实生活中滑动面是不断变化的，本文假设其为一个均匀的滑动面，设滑坡滑动面的坡向轴为$OA$，沿坡向$OA$向下为正，垂直于滑动面的法线方向$OC$向上为正，垂直于面AOC的$OB$轴向外为正构建以滑动面为参考面的独立坐标系。因为滑动面坐标系与站心坐标系（east/north/up，E/N/U）存在几何关系，见图 2。因此可将站心坐标系作为过渡，由LOS的年形变速率得到$OA$、$OB$、$OC$方向的年形变速率。在自然环境下，除了借助外力，一般沿滑动面的法线方向$OC$形变远远小于其他方向，所以本文将$OC$方向的结果视为零^[28]。

图  2  滑动面坐标系与站心坐标系的关系

Figure 2.  Relationship Between the Sliding-Surface-Defined Coordinate and Topocentric Coordinate Systems

由此得到三轨雷达影像的LOS向年形变速率，沿滑坡向（$OA$）、垂直滑坡向（$OB$）的年形变速率计算如下：

$\mathrm{式}\mathrm{中},V_{\mathrm{L}\mathrm{O}\mathrm{S}}$为SAR卫星LOS向的年形变速率；$\theta$为雷达卫星的入射角；$\alpha$为雷达卫星的方位角；$V_{\mathrm{N}}$、$V_{\mathrm{E}}$、$V_{\mathrm{U}}$分别为在站心坐标系下N、E、U方向的年形变速率；$\beta$为坡度；$\gamma$为坡向；$V_A$、$V_B$、$V_C$分别为滑坡向、垂直滑坡向、滑坡面法向（即垂直方向）的年形变速率。

本文将GACOS辅助下的InSAR Stacking和LiCSBAS技术分别应用于西藏自治区贡觉县雄松乡至沙东乡的金沙江流域进行滑坡群探测与监测，利用2017—2019年Sentinel-1和ALOS-2 3个轨道约2 a的影像对该滑坡群进行探测得到LOS向年形变速率，结果如图 3所示，负值表示远离雷达视线方向，正值表示靠近雷达视线方向。图 4为沿坡向和垂直滑坡向的年形变速率。

图  3  InSAR LOS向年形变速率图

Figure 3.  Annual Deformation Rates in the LOS Direction Revealed by InSAR

图  4  沿坡向$OA$和垂直滑坡向$OB$的年形变速率

Figure 4.  Annual Deformation Rate in the$OA$ and $OB$ Directions

由图3、4可以看出，该滑坡群有多处均发生明显形变，两种方法所得结果的量级和滑坡范围基本一致。将两种方法所得结果进行相关性分析，结果如图 5所示。

图  5  GACOS辅助下InSAR Stacking和LiCSBAS技术LOS向年形变速率相关性分析

Figure 5.  Correlation Analysis Between GACOS Assisted InSAR Stacking and LiCSBAS Derived LOS Annual Deformation Rates

由图 5可知，两者结果具有一致性，LOS向年形变速率的$R^{\mathrm{2}}$在0.92以上，沿坡向和沿垂直滑坡向年形变速率的$R^{\mathrm{2}}$在0.85以上。ALOS-2影像所得结果的标准差为1.3 mm/a，Sentinel-1A升轨影像所得结果的标准差为5.4 mm/a，Sentinel-1A降轨影像所得结果的标准差为6.8 mm/a，而由图 1可知，沿坡向年形变速率标准差为7.2 mm/a，沿垂直滑坡向年形变速率标准差为7.4 mm/a，这些均证明了GACOS辅助下InSAR Stacking结果的可靠性。

为了对结果进行更加详细和直观的分析，在形变量较大的滑坡上随机选取3个剖面线a、b、c，并根据滑坡探测所得结果，将该滑坡群分为$A$、$B$、$C$、$D$、$E$、$F$、$G$ 7个区域，见图 6。图 7为两种方法沿剖面线a、b、c获得的LOS向年形变速率对比。由图 7可知，两种方法所得结果相关性都在0.9以上，尤其是Sentinel-1升轨和ALOS-2影像的结果相关性都在0.95以上。对于ALOS-2影像，两种方法计算结果的RMS均在2 mm/a以下，Sentinel-1升轨影像的两种方法解算结果的RMS均在2~4 mm/a之间；相较于其他两轨道的对比结果，Sentinel-1降轨影像的两种方法解算结果的RMS较大，剖线a的RMS为8.5 mm/a，剖线b的RMS为4 mm/a，剖线c的RMS为3.7 mm/a，这可能是由于降轨影像采集的时间在白天，而升轨影像采集的时间在晚上，白天大气水汽时空变化可能比较大，扰动型误差占比较大，而GACOS改正后的残余也随之变大。通过以上对比分析更加进一步验证该区域滑坡识别的可靠性。

图  6  研究区域剖面线a、b、c和滑坡区域A~G

Figure 6.  Profiles a, b, c in the Study Area and Landslide Areas A-G

图  7  沿剖面线a、b、c的GACOS辅助下InSAR Stacking和LiCSBAS技术获得的LOS向年形变速率对比

Figure 7.  Comparison Between GACOS Assisted InSAR Stacking and LiCSBAS Derived LOS Annual Deformation Rates Along Profiles a, b, c

图 8为滑坡区域A~G的光学遥感影像，通过谷歌地球将滑坡区域标记，目视解译发现7个区域均有滑坡迹象，该滑坡群在相同的监测时间段内表现出不同的区分性。

图  8  滑坡区域A~G的光学遥感影像

Figure 8.  Landslide Areas A-G Revealed by Optical Imagery

由图 8可知，区域A、B、C、D为较活跃滑坡，具有整体失稳的可能。区域A滑坡覆盖范围约为8.5$\mathrm{k}{\mathrm{m}}^{\mathrm{2}}$，通过光学影像可以看出该区域为积层滑坡，在滑坡的前后缘均存在明显的形变，尤其是在滑坡前缘区域沿坡向最大年形变速率为$-$142 mm/a，可能原因是常年受河水的侵蚀，导致局部发生垮塌，可判定为老滑坡；区域B滑坡覆盖范围约为3$\mathrm{k}{\mathrm{m}}^{\mathrm{2}}$，沿坡向最大年形变速率为$-$82 mm/a，通过光学遥感影像可以判断该滑坡为岩质滑坡，该滑坡后缘裂缝明显，前缘上部多级滑动交错，中部滑动明显；区域C和区域D滑坡覆盖范围分别约为12$\mathrm{k}{\mathrm{m}}^{\mathrm{2}}$、2.4$\mathrm{k}{\mathrm{m}}^{\mathrm{2}}$，同属于积层滑坡，尤其是区域C前后缘均存在明显的形变，该区域属于迎风坡，常年受强烈的太阳光照，坡脚被江水掏蚀严重，沿坡向最大年形变速率为$-$163 mm/a。区域F属于岩质滑坡，前后缘均形变缓慢，区域E、G为积层滑坡，目前区域E、F、G滑坡仍处于缓慢变形阶段。区域A、B、C、D形变明显，在强降雨江水的侵蚀以及地震等诱发下，可能会整体失稳造成堵江后果。后续应继续对该区域进行InSAR时间序列分析，在条件允许情况下进行现场设备监测，避免不必要的损失。

本文对西藏自治区贡觉县雄松乡至沙东乡的金沙江流域范围进行滑坡探测与监测，获取该区域的LOS向、坡向和垂直滑坡向年形变速率，沿坡向最大年形变速率为-163 mm/a，探测结果较可靠。并通过谷歌地球光学遥感影像对滑坡群进行判别，发现了一个包含7处滑坡隐患点（A~G）的滑坡群。该滑坡群形成的原因包括：（1）研究区域海拔为4 000 m左右，距离冰川近，近年来气温上升，冰川融水增加，导致以冰川为驱动力的滑坡数量增加；（2）在2018-10-11和2018-11-03，研究区域上游白格村附近连续两次发生大规模高位滑坡形成堰塞湖，在人工泄洪时，对该滑坡群前缘带来较大影响。探测到7个隐患点具有较明显的形变，区域A、B、C、D为较活跃滑坡，可能会整体失稳造成堵江后果，目前区域E、F、G滑坡仍处于缓慢变形阶段。GACOS辅助下InSAR Stacking技术可有效去除长波段及与地形相关大气延迟误差，具有简单、有效、快速、易于推广应用等优点。