白格滑坡处于中国四川省甘孜州白玉县与西藏昌都市江达县交界处，行政上隶属于江达县波罗乡白格村，属于青藏高原到四川盆地的过渡地带，滑坡区地势起伏较大，坡顶海拔约3 720 m，河谷高程2 880 m，垂直高差达840 m^[9]，易发生滑坡等地质灾害。

本文采用了从美国地质调查局（United States Geological Survey，USGS）官网获取的2015-11—2018-04期间的11景哨兵2号数据和2014-11—2018-09期间的16景Landsat 8陆地成像仪（operational land imager，OLI）数据，详细参数信息分别如表 1、表 2所示。两平台数据的最高空间分辨率分别为15 m和10 m。值得注意的是，哨兵2号和Landsat 8影像产品都已通过外部数字高程模型（digital elevation model，DEM）进行了正射校正^[11-12]。

本文对收集的哨兵2号和Landsat 8两平台

光学数据进行偏移量估计，获取其东西和南北向的形变，并对其中的一系列误差^[13-14]进行修正。与此同时，为了获取该滑坡时序形变，本文借鉴时序InSAR技术中短基线^[15]的思想来减少时空失相关、太阳高度角和方位角等对结果精度的影响。为了提高形变的时间分辨率，本文将两类光学数据联合解算，获取长时间序列形变，为滑坡的发育、发生和后续的预测提供丰富的监测信息。此外，通过统计稳定区域形变以及独立数据结果的交叉验证两种方式进行精度评定。本文的研究流程包含影像对选择、偏移量估计、误差后处理、时序偏移量联合解算、精度评定5个关键步骤。

光学影像匹配结果易受云雾或长时间地表植被变化等引起的失相关噪声、轨道误差、条带误差、地形阴影误差以及卫星姿态角误差影响^[13-14]。为了减少这些误差对偏移量估计的影响，本文选取月份数值差不大于3的影像对，以降低时间基线引起的噪声^[7，12]。哨兵2号和Landsat 8的空间基线阈值分别设置为0 m和1 000 m。辐射长度阈值设置为0.8h（h为像素地形高度）。其中，为了保证影像对的连通性，本文在上述选择基准上为Landsat 8添加了5个影像对，所有影像对的参数如图 1所示。

图  1  哨兵2号和Landsat 8影像对参数

Figure 1.  Parameters of Sentinel-2 and Landsat 8 Image Pairs

本文采用COSI-Corr（the Coregistration of Optical Sensed Images and Correlation）软件的基于亚像素频率域相关性匹配算法来获取其东西向、南北向形变，处理步骤及相关参数设置如下。

首先，为了计算多源数据联合的时序偏移量，本文将哨兵2号的近红外波段降采样至与Landsat 8全色波段一致的15 m空间分辨率，通过人为选取同名点结合最小二乘的方法配准^[16]。然后，利用COSI-Corr对滑动窗口大小（包括初始和最终相关搜索窗口）、步长、迭代次数和掩模阈值等参数进行设置。本文结合文献[6]的建议将两个搜索窗口都设置为64×64像素，步长设置为1像素，迭代次数设置为4，掩模阈值选择0.9。最后，采用软件批处理工具计算偏移量，每个影像对结果都包括东西向形变图、南北向形变图和信噪比（signal⁃to⁃noise ratio，SNR）图。正数代表向北和向东的形变，信噪比范围为0~1，其值越高代表越好的匹配质量。除了东西向和南北向，将欧氏距离作为总位移^[6，8]。

由§2.1可知，软件获取的偏移量包含失相干噪声、轨道误差、条带误差和卫星姿态角误差，为保证结果的精度，需要对各方向偏移量误差进行后处理。对于失相干噪声，本文通过设置SNR阈值（> 0.96）和手动掩模（针对由局部云量及其阴影产生的失相关）相结合的方法去除；对于轨道噪声，通过一次多项式曲面拟合模型方法去除。条带误差采用均值相减法去除。卫星姿态角误差采用改进的均值相减法^[14]去除。此外，还需要进行非局部均值滤波、坐标系转换、偏移量值重采样和二维水平形变场获取，本文结合文献[17]建议的默认参数进行非局部均值滤波。

短基线集时序方法最早由文献[15]提出，后被应用于光学影像时序形变监测。以东西向形变为例，假定有（N+1）幅同一区域的光学影像，获取时间按{t₀，t₁$\dots$t_n}排列，其中，t₀表示初始影像时间，M为互相关影像对。主从影像分别记为$m_i$和$s_i$。以形变速率作为未知参数构建观测方程为：

式中，系数阵$\mathit{B}$为M×N的矩阵，其元素由相邻影像的时间间隔组成，假设任意一个元素$B\left(i,j\right)$满足：

式中，$i=\mathrm{1,\; 2}\dots M$；d为形变观测值，即偏移量估计获取的形变值，是一个M×1的矩阵；$\mathit{v}$为待求参数速率矢量，$\mathit{v}={\left[v_{\mathrm{1}}{v}_{\mathrm{2}}\dots v_n\right]}^{\mathrm{T}}$。

一般情况下，当系数阵$\mathit{B}$满秩时，$\mathit{v}$通常采用最小二乘就能解算出来：

式中，权重矩阵$\mathit{P}$为一个M×M的矩阵，P=$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(w_{\mathrm{1}},w_{\mathrm{2}}\dots w_M)$，其中，

式中，$∆t$为影像时间间隔；w_i为每个干涉对的权重。最终通过每个时段形变速率乘以时间间隔来获得形变量。但当系数阵B秩亏时，即两两偏移量的影像组成的网络不连通时，此时需要用到加权奇异值分解（singular value decomposition，SVD）^[18]来求解未知参数$\mathit{v}$，即Landsat 8和哨兵2号两个独立的子集。

式中，svd为奇异值分解函数；U为$\mathit{Z}{\mathit{Z}}^{\mathrm{T}}$的特征向量，是一个M×M的正交矩阵；S为一个M×M的对角矩阵；V是一个N×M矩阵，为${\mathit{Z}}^{\mathrm{T}}\mathit{Z}$的特征向量；S⁺为S的逆；${\sigma }_i$为奇异值；L-1表示奇异值为0的数目。

由于缺乏实测数据，通常采用先验信息辅助进行精度评定，比如利用稳定区域的形变来验证。当用到两类数据源时，交叉验证也是常用的手段。本文在研究区随机选择了3个相同大小的稳定区域并结合文献[12]提出的误差评估方法，统计各影像对相比于初始影像对偏移量和时序偏移量结果。在交叉验证方面，选择哨兵2号和Landsat 8相近时段的偏移量结果进行对比分析。

本文时序形变的时间跨度约3.77 a。由已有成果可知白格滑坡滑移方向以东西向为主，南北向形变很小^[6]。因此，本文展示总形变的时序偏移量结果，包括哨兵2号（10个）和Landsat 8（15个）时序形变结果、双平台时序形变结果（26个）以及哨兵2号和Landsat 8平均形变速率结果，分别如图2~4所示。

图  2  哨兵2号（（a）~（j））和Landsat 8（（k）~（y）)时序偏移量

Figure 2.  Time-Series Deformation of Sentinel-2((a)—(j))and Landsat 8((k)—(y))

图  3  哨兵2号和Landsat 8联合时序偏移量

Figure 3.  Combined Time-Series Deformation of Sentinel-2 and Landsat 8 Results

图  4  哨兵2号和Landsat 8时序偏移量平均形变速率

Figure 4.  Time-Series Deformation Velocity for Sentinel-2 and Landsat 8

由图 2可知，哨兵2号计算了2015-11-23—2018-04-16白格滑坡时序形变，期间最大累计形变23.45 m；Landsat 8计算了2014-11-29—2018-09-05白格滑坡时序形变，期间最大累计形变45.43 m。

由图 3可知，双平台时序形变包括2014-11-29—2018-09-05期间的26个时序形变，截至2016-12-04、2017-11-07和2018-09-05，最大累计形变量分别为11.38 m、20.43 m和42.90 m。该时序形变结果包含了更多的形变细节，改善了滑坡监测的时间分辨率。

由图 4可知，哨兵2号和Landsat 8的东西向、南北向和总的平均形变速率的最大值分别为9.56、1.37、9.59 m/a以及9.16、1.07、9.18 m/a。两平台在3个方向上的平均形变速率差异都不大，只是滑坡边界略有不同，这与偏移量计算的影像对数量、质量和时间间隔有关。

本文从研究区域随机选择3个稳定区域（见图 5），并统计3个稳定区域在滤波后、单平台和双平台时序计算下的标准差来进行精度评定，结果如图 6所示。

图  5  选择的3个稳定区域位置

Figure 5.  Selected Three Stable Areas

图  6  哨兵2号与Landsat 8误差改进比例比较

Figure 6.  Deformation Improvement Comparison Between Sentinel-2 and Landsat 8

由图 6可知，大部分影像对通过单一影像对偏移量、单平台时序偏移量和双平台时序偏移量计算的东西向和南北向标准差都有所减小。其中，哨兵2号东西向分别平均提升21.47%、24.50%和27.52%，南北向分别平均提升24.59%、24.97%和30.34%；Landsat 8东西向分别平均提升33.53%、43.33%和46.94%，南北向分别平均提升36.60%、40.64%和41.04%。总体而言，Landsat 8不确定性改进效果优于哨兵2号。双平台时序偏移量精度相对于单平台都有所提高，哨兵2号东西向和南北向分别提升3.02%和5.37%；Landsat 8东西向和南北向分别提升3.61%和0.40%。

除了稳定区域统计值，本文还选择了哨兵2号和Landsat 8相近时间段的偏移量结果（见图 7）来进行交叉验证，影像对参数和结果统计分别如表 3和表 4所示。由图 7和表 4可知，哨兵2号无论平均值和标准方差都小于Landsat 8。前者标准方差在0.5 m左右，后者标准方差在0.75 m左右。根据COSI-Corr软件，互相关匹配精度最高可达1/20个像素^[19]，即哨兵2号和Landsat 8形变监测理论精度为0.5 m和0.75 m。两影像对参数基本一致，因此，就单影像对监测精度来说，哨兵2号优于Landsat 8，在单一数据下，应首选哨兵2号。

图  7  哨兵2号和Landsat 8影像对形变结果

Figure 7.  Deformation Results for Sentinel-2 and Landsat 8 Pairs

为了分析白格滑坡的演化发育过程，本文在坡顶、坡中部和坡底分别选择1个（P₁）、3个（P₂、P₃和P₄）和1个（P₅）特征点分析形变时间序列，5个特征点位置及其累计形变如图 8所示。由图 8可知，特征点P₁和P₅分别处于坡顶和坡底，总的累计形变分别为36.70 m和32.14 m。特征点P₂、P₃、P₄都处于滑坡中部，总的累计形变量分别为39.27、42.90、38.13 m。值得注意的是，2018-04-16—2018-09-05期间5个特征点形变速率差异较大，分别为93.63、71.8、108.6、168.1、48.94 mm/d。坡底特征点形变速率最小，坡顶次之，坡中部平均最大，该时段形变特征基本符合推移式滑坡自上而下的特点^[20]。

图  8  5个特征点位置及其累计形变

Figure 8.  Time-Series Deformation of Five Selected Points

根据蠕变变形理论^[21]，滑坡发生前通常经历初始变形阶段、等速变形阶段和加速变形阶段（细分为初加速、中加速和加加速3个亚阶段）。文献[22]提出将切线角α_i作为3阶段划分依据，当切线角远小于45°认为是滑坡初始变形阶段；当切线角约等于45°认为是滑坡等速变形阶段；当切线角远大于45°认为是滑坡加速变形阶段（80°和85°分别为3个亚阶段判定依据）。为了分析此滑坡的形变阶段，本文选择累计形变量最大的特征点P₃结合切线角分析累计形变与时间关系（见图 9）。由图 9和滑坡3阶段特征，本文将滑坡划分为初始变形（2014-11-29—2015-09-29）、等速变形（2015-09-29—2015-12-23）和加速变形3个阶段，其中加速变形又可细分为初加速（2015-12-23—2018-02-05）、中加速（2018-02-05—2018-09-05）2个亚阶段。整个时段形变速率逐渐增加，其中，初始变形阶段的平均形变速率约7.05 mm/d，切线角约36.29°；等速变形阶段的平均形变速率约10.27 mm/d，切线角约46.91°；初加速阶段的平均形变速率约26.80 mm/d，切线角约70.28°；中加速阶段的平均形变速率约65.79 mm/d，切线角约81.69°。

图  9  白格滑坡发生前3.77 a形变阶段划分

Figure 9.  Deformation Stage Division of Baige Landslide Before the Occurrence

为了分析形变与降水的关系，本文收集了来自全球气象站距离白格滑坡最近的德格县站点（98.567°E，31.733°N）的逐年、逐月和逐日的降雨数据。特征点P₃形变速率、累计形变与降雨数据如图 10所示。

图  10  白格滑坡2014-11至2018-09降雨量、形变速率与累计形变

Figure 10.  Rainfall, Deformation Rate, and Cumulative Deformation of the Baige Landslide from November 2014 toSeptember 2018

由图 10可知，2015—2017年累计形变量较小，2018年累计形变量较大，形变速率逐渐增大。此外，本文对雨量和累计形变（图 10中的形变点）进行了相关性分析，其相关系数R²为0.88，表明形变与降雨高度相关，但该滑坡可能还受其他内外部因素^[23]作用。

本文联合哨兵2号和Landsat 8的光学数据计算了白格滑坡2014-11—2018-09二维水平时序形变，比较了不同平台的偏移量形变精度差异，根据蠕动变形理论、切线角对时序偏移量进行阶段划分，分析了形变与降雨的相关性。

结果表明：（1）多源光学影像时序偏移量跟踪技术可改善滑坡监测的时间分辨率，获取更多形变细节。（2）双平台时序偏移量精度高于单平台，哨兵2号东西向和南北向分别提升3.02%和5.37%；Landsat 8东西向和南北向分别提升3.61%和0.40%。此外，相同观测条件下，哨兵2号监测精度高于Landsat 8。（3）白格滑坡在2014-11—2018-09期间的时序形变可划分为初始变形、等速变形和初加速变形和中加速变形阶段，平均形变速率分别为7.05、10.27、26.80、65.79 mm/d。（4）白格滑坡形变与降雨因素高度相关，可能受内外因素共同作用。总体来说，本文证实了双平台时序偏移量在滑坡监测中具有提高时间分辨率和精度的优势，未来多平台时序偏移量估计可能是该方法的趋势。此外，多平台时序偏移量有助于滑坡阶段划分，把握其演化过程。

此外，由于不同平台光学数据传感器对地表形变的敏感度存在差异，多源光学数据联合解算时序偏移量方法存在交叉时段形变波动等噪声，后续研究可以考虑影像监测能力并引入方差分量估计加权来减小此误差。