青藏高原是全球中低纬度海拔最高、面积最大的多年冻土分布区，现存多年冻土面积约130×10⁴ km²，约占高原总面积的56%(赵林等，2010).青藏高原冻土区是中国多数大河流域的发源地，为干旱和半干旱地区提供了主要水源(王庆锋等，2016).近年来，青藏高原变暖趋势明显，多年冻土出现活动层增厚、范围萎缩、厚度减薄等显著退化现象(Wang et al., 2000).多年冻土活动层内的水分随着温度的变化在固态和液态之间不断转换，土层发生冻胀融沉，导致多年冻土区大范围地面变形，严重破坏区域内重大工程基础设施(如青藏铁路、青藏公路)(Chen et al., 2012; Li et al., 2015).同时，活动层变化导致冻土区水文地质条件发生改变，影响寒区河流的地表径流和水资源调蓄功能，进而影响区域水文过程及生态环境(程国栋和金会军，2013).因此，为深入认识青藏高原冻土退化对人类工程活动和区域可持续发展的影响，对多年冻土区活动层厚度变化开展大范围、高分辨率监测具有十分重要的意义.

当前，野外原位观测地表形变的监测方法，如水准测量、全球导航卫星系统(GNSS)、测斜仪等(袁堃等，2013; Liu and Larson, 2017)，在冻土形变监测中具有精度高、布设方便等优势，但需要现场布设且监测点稀疏，难以有效开展大范围、长时间的冻土形变监测.作为近年来发展迅速的空间大地测量新技术，星载合成孔径雷达干涉测量(InSAR)具备全天候、高精度、大范围且无需布设地面控制点等技术优势，已逐渐应用于美国阿拉斯加、我国青藏高原和东北地区等多年冻土区的地表形变监测研究中.Wang和Li(1999)开创性地将差分干涉测量技术(DInSAR)应用到阿拉斯加冻土的冻胀形变研究中，但DInSAR方法易受大气扰动、地形误差和干涉失相干的影响，导致监测结果的精度和可靠性较低.为了降低大气延迟和失相干的影响，Xie等(2009)采用永久散射体雷达干涉测量(PS-InSAR)技术和Envisat/ASAR数据对青藏高原北麓河段青藏铁路进行监测，获得了与水准测量一致的结果.Daout等(2017)利用Envisat/ASAR数据进行小基线集(SBAS)时序InSAR分析，获取了西藏西北部超过60000 km²区域跨越8年的冻土区地表形变信息.近年来，国内学者先后利用DInSAR、PS-InSAR和SBAS等多种InSAR技术在我国多个冻土地区开展一些系列研究并取得了可靠的结果(李震等，2004; Liu et al., 2010; 胡波等，2010; 李珊珊等，2013; Chen et al., 2013; 王春娇，2015; Zhao R et al., 2016; Jia et al., 2017; 张正加，2017).

然而，PS-InSAR只针对不受时间、空间去相干的强反射点目标(如人工建筑物、裸露岩石等)进行干涉分析，在冻土地区的永久散射体(Persistent Scatterers，PS)十分有限，导致观测结果不能客观反映研究区整体变化(蒋弥等，2016).SBAS技术通过选取短时-空基线的干涉对减小时空失相干的影响，并将分布式雷达目标或称分布式散射体(Distributed Scatterers，DSs)作为时序分析对象以提高形变监测点的密度，但通常情况下冻土地区DSs目标的信噪比较低进而影响形变监测结果的可靠性(Fornaro et al., 2015).为进一步提高自然地表DSs点密度和时序形变估计质量，以SqueeSAR技术(Ferretti et al., 2011)及其变体技术(Wang et al., 2011; Goel et al., 2012; Verde et al., 2018)为代表的新一代DSs-InSAR方法近年来备受瞩目，该类方法可以很好地弥补PS-InSAR和SBAS方法的不足，但DSs点识别和相位优化运算十分耗时，且SAR影像数较少时(一般小于25景)假设检验显著性不高(蒋弥等，2016).此外，当前研究中常用的SAR影像(如ASAR、ERS、PALSAR等)，在我国大部分冻土区获取时间不均匀且数据量较少，难以获得完整的冻融形变过程，需提供先验知识对InSAR解算模型进行约束(如正弦模型和Stefan模型)，容易引入较大的模型误差(李珊珊等，2013).欧洲航空局(ESA)分别于2014年和2016年发射的Sentinel-1A/B SAR卫星，具有短重复周期(6/12/24天)、大幅宽、多极化等特点(Torres et al., 2012)，在一定程度上有利于减弱时间去相干的影响，并能够获取更丰富的冻土区域地表变化过程信息，对于多年冻土形变监测和活动层厚度反演具有很好的应用前景.

因此，本文提出一种基于分布式目标(DSs)的小基线集时序InSAR(DSs-SBAS)冻土形变监测方法，利用分布式目标时序InSAR方法和温度-形变约束关系，较好地克服了冻土区极易失相干和参考点选取的难题，可以提高冻土形变监测的时空分辨率和可靠性.以祁连山黑河流域上游西支源头野牛沟多年冻土区为研究区，利用该方法和27景Sentinel-1A数据，获取该区域大范围、高时空分辨率的冻土形变时间序列和年均形变速率；并使用Stefan模型联合地温数据估算其季节性形变幅度.此外，本文将分析研究区冻土形变的空间分布特性以及典型冻土地貌的冻融过程及其影响因素.

野牛沟位于青藏高原祁连山区黑河流域的西支源头(98.67°E—99.27°E, 38.55°N—39.02°N)，是我国西北第二大内陆河黑河流域的主要产流区，海拔3500~4700 m；属阿尔金山—祁连山高寒带山地多年冻土区，广泛发育着多年冻土和季节性冻土(周幼吾等，2000)；主要受西风环流影响，年平均气温为-4.1 ℃，年降水量约286 mm (Wang et al., 2015).研究区域为面积约为2000 km²纵向河谷，与山脉走向平行，沟底地形平坦开阔，越靠近两侧山脉坡度越大；二尕公路穿越其中，路面宽度为7 m，区域内长约50 km.主要植被类型为高寒草甸和高山沼泽，土壤类型主要为砂质土、亚黏土和黏土，土壤粗糙有机物含量少(图 1).区域内湿地广泛发育，夏季积雪融水和雨水下渗受下伏多年冻土层阻挡，聚集在冻胀草丘和洼地间.区域内多年冻土受坡面、植被、土壤、积雪等局部因素影响，广泛分布在海拔3650 m以上区域，该海拔以下区域发育着季节性冻土(Wang et al., 2013).研究区多年冻土融化/冻结始日分别在每年的4月中旬和10月下旬(王庆锋等，2016)，活动层厚度在1.07~4.86 m之间，平均值为2.72 m，且在不同植被类型间连续变化(Cao et al., 2017).

图 1

Fig. 1

图 1 野牛沟研究区域及土壤分类 Fig. 1 Study area and soil-type classification of the Yeniugou

本文选用欧空局Sentinel-1A卫星获取的27景升轨单视复数影像(表 1).数据获取时段为2014年10月19日至2016年10月20日，极化方式为VV极化，所处波段为C波段，方位向和距离向分辨率分别为20 m和5 m.原始数据的幅宽为250 km，选取了原始数据中长宽约47 km和43 km的感兴趣区域.为了提高Sentinel-1影像的轨道精度，采用了欧空局发布的精密轨道数据(PDGS).此外，使用美国宇航局发布的分辨率为30 m的SRTM数据作为外部DEM数据，用以去除干涉图的地形相位和地理编码；使用美国宇航局发布的SMAP和Sentinel-1数据反演得到的1 km分辨率的土壤湿度产品(图 4)辅助参考点选取，数据获取时间为2018年6月16日，均方根误差约为0.04 cm³·cm^-3 (Das et al., 2017, 2018).图 2为影像按照小基线原则组合而成干涉对的时空基线图.

表 1 (Table 1)

表 1 所用Sentinel-1数据获取时间 Table 1 Sentinel-1 SAR scenes used in this study No. Date 1 20141019 2 20141112 3 20141206 4 20150123 5 20150216 6 20150228 7 20150324 8 20150511 9 20150604 10 20150628 11 20150815 12 20150908 13 20151002 14 20151026 15 20151213 16 20160106 17 20160130 18 20160223 19 20160318 20 20160411 21 20160505 22 20160529 23 20160716 24 20160809 25 20160902 26 20160926 27 20161020 Note: The pixel spacing is 5 m in the slant range direction, 20 m in the azimuth direction.

表 1 所用Sentinel-1数据获取时间 Table 1 Sentinel-1 SAR scenes used in this study

图 2

Fig. 2

图 2 干涉对时空基线(蓝色节点表示SAR观测，连接线表示对应的干涉对) Fig. 2 Temporal-spatial baselines of interferometric pairs. Blue dots are SAR observations. Link lines are corresponding interferometric pairs

图 4

Fig. 4

图 4 土壤湿度及参考点候选区域 Fig. 4 Soil moisture and candidate areas of reference points

此外，图 1所示的PT3、PT4、PT5三个站点土壤温度数据，用于Stefan模型估算季节性形变幅度.这些站点的土壤温度观测采用Campbell公司生产的109土壤温度传感器，量程为-50~70 ℃，线性误差小于0.03 ℃(王庆锋等，2016).温度数据的获取时间与影像同步，为距地表以下5 cm处获取的数值.此外，为了实地验证InSAR结果，我们于2017年6月15—17日在野牛沟进行了为期3天的野外踏勘.

本文提出了一种基于分布式目标的小基线集时序InSAR(DSs-SBAS)冻土形变监测方法，方法流程如图 3所示，主要包括三部分：(1)取单视复数影像的幅度序列，进行DS点识别；(2)将幅度数据进行辐射定标和距离多普勒地形纠正；基于分布式目标进行自适应相干性估计，使用特征值分解算法优化分布式目标相位；根据土壤湿度、相干性结果和校正后幅度数据，进行参考点选取；(3)使用奇异值分解(SVD)模型计算年均形变速率和形变时间序列.此外，大气误差是InSAR数据处理中不可忽略的一个重要因素.本文主要通过建立与地形相关的模型采用参数估计的方法去除地形相关的大气误差；通过对多个干涉图滤波估计剩余大气误差并去除；最后通过使用特征值分解法，将残余的大气误差作为噪声的一部分做进一步去除.下文重点阐述DSs点选取、相位优化和参考点选取三个关键技术环节.

图 3

Fig. 3

图 3 DSs-SBAS时序InSAR冻土形变监测的技术流程 Fig. 3 Flow chart of time-series InSAR method for monitoring permafrost deformation

Jiang等(2015)提出了在高斯假设背景下将假设检验转换为置信区间估计的DSs点选取方法，大大提高了选点效率，并且在样本数较少的情况下，依然能保持较高的准确性.简言之，给定N个独立随机变量x_i, x_i服从瑞利分布根据中心极限定理，在N足够大的情况下，样本均值x= 服从均值为, 方差为δ_x²= 的高斯分布.因此X的区间估计为

(1)

为标准正态分布的1-α/2分位点，X表示参考点样本均值.根据SAR图像统计特征(Lee and Pottier, 2009)，单视幅度影像服从瑞利分布，其变异系数.因此，将(1)式转换为只包含μ_x的区间.

(2)

在实际处理中，我们首先假设μ_x≈X，取显著性水平为0.5以减小分布尾部的估计偏差，且不考虑像元之间的连接性，将所有落入区间的样本点都视为分布式目标点集Ω_init.将Ω_init中的像元值平均之后再重新以0.05的显著性水平估计DSs候选点(DSC)，获得新的DSC点集Ω.再对集合中所有点进行连通性评估，只保留符合八连通的点.在获取DSC点集之后，对每一个点计算其自适应相干系数

(3)

Ω表示DSC点集合，Z_iZ_j^*表示对两复数影像共轭相乘.由于样本数量较少，得到的相干系数是有偏的，因此使用第二类统计估计得到一个无偏的结果(Abdelfattah and Nicolas, 2006).

(4)

这里〈·〉_k∈Ω代表取平均.若某一点在80%的干涉对中的自适应相干系数都大于0.5，则将此点保留为DSs点.

青藏高原冻土区大多覆盖着低矮的高原草甸或发育着冻胀草丘等冻土地貌，并且每年的冻融循环及春季降雪极易造成SAR干涉对失相干.在这种低相干区域，Fornaro等(2015)提出了一种基于特征值分解的方法来改善时序InSAR处理结果.这种方法的关键在于对干涉数据集进行特征值分解，提取相位观测值中相互独立的主要成分，并对每一个像元干涉组合的协方差矩阵进行特征分解.本文采用该方法从干涉相位协方差矩阵中估计相位观测值，每一个特征值对应着特定的形变特性以及目标特性.协方差矩阵为：

(5)

Z=[z₁, z₂, z₃, …, z_N]^T表示N景影像归一化后的观测值，E[|z_i|²]=1，H代表共轭转置.

(6)

(7)

这里，e^jϕ_mn表示第m景影像和第n景影像的干涉相位，为无偏的自适应相干系数.根据特征值分解理论，我们可以将协方差矩阵C表示为C=UΛU^-1=UΛU^T.这里U=[u₁, u₂, …, u_N]为特征值对应的特征向量，Λ=diag(λ₁, λ₂, …, λ_N)是N个特征值组成的对角阵，其中λ₁≥λ₂≥…≥λ_N.因此，我们可以将协方差矩阵写为N个独立散射体的和，

(8)

特征值代表每个散射体的权重，假设这里存在n(0＜n＜N)个较大的非零特征值(λ₁, λ₂, …, λ_n)对应着n个不同的散射机制，和N-n个较小的特征值(λ_n+1, λ_n+2, …, λ_N)对应着由于失相干造成的噪声.因此可以将(8)式中协方差矩阵进一步写为

(9)

在实际应用过程中，我们并不知道主散射体的个数.由于本文研究的是极易失相干的冻土区域，第二个特征值对应的信号很可能只是噪声.因此，本文仅采用最大特征值对应的特征向量进行分析.值得注意的是最大特征值对应的特征向量并不唯一且还存在一个固定项，因此我们将最终的结果表示为干涉相位的形式.

将InSAR的相对观测结果转换为绝对值，需要选定一个或若干个绝对形变已知的点(Hanssen, 2001).但在野牛沟冻土区难以找到稳定的裸岩或已知形变点作为参考点，因此我们提出一种地温-形变约束下的冻土区参考点选取新策略.冻土地表形变和该区域温度呈现强烈的负相关(Zhao R et al., 2016)，而稳定区域的形变序列和温度序列的相关性基本为零.据此，可以根据形变和温度的相关性来判定区域地表的稳定性.另外，冻土地表形变与土壤湿度具有一定的相关性，冻土区域干燥地表保持稳定的可能性较高.因此，本文使用NASA发布的覆盖研究区域的土壤湿度产品作为参考，数据获取时间为2018年6月16日.冻胀期土壤冻结地表土壤湿度均保持在较低水平，因此融沉期的土壤湿度可以很好地代表区域内土壤湿度分布状况.土壤湿度分布如图 4，参考点的选取需要进行以下步骤：首先根据土壤湿度数据选择土壤湿度小于0.2 cm³·cm^-3的区域；然后，根据相干性、校正后的幅度数据和区域地形地貌，将图 4所示的四个区域(R₁、R₂、R₃、R₄)作为待选区域.值得注意的是区域R₄非大面积干燥区域，但此处有人工建筑存在，在干涉对中呈现高相干性，故将其选为参考点候选区域.最后，我们依次选择每个区域中相干性高、地形平坦的点(r₁、r₂、r₃、r₄)作为InSAR解算参考点，剩下的区域作为验证集，计算基于某一个参考点下其他三个区域所有点的形变序列与温度的相关性.图 5中a图为当参考点选r₁时，R₂、R₃和R₄区域所有点的形变序列与温度的相关系数直方图.r₁的相关系数基本位于零值附近(图 5a)，表明此参考点较为稳定；而当参考点选r₂、r₃和r₄时，其相关系数的统计分布较大程度偏离零值(图 4中b、c、d)，说明这些参考点存在较大形变.

图 5

Fig. 5

图 5 不同参考点下形变序列与温度序列相关系数直方图 Fig. 5 Correlation coefficient histogram between the series of deformation and temperature under different reference points

目前Stefan方程的应用主要有两种方式，一种以时间为自变量，此模型在冻融季交替处会出现跳变，违背了冻土实际形变规律；一种以温度为自变量，此模型更加接近Stefan方程的原理和冻土形变的物理机制，在冻融季交替阶段更加贴合事实，结果较为可靠(Liu et al., 2012; Daout et al., 2017).本文使用基于温度的Stefan模型，从两年(2014-10-17—2016-10-20)的形变时间序列中估算平均冻融形变幅度.

(10)

(11)

这里t_f和t_t分别表示冻胀季和融沉季的开始时刻，A为年均形变速率，B为季节性形变幅度，c为常数.DDF和DDT分别是全年日均温度大于零和小于零两个区间的日均温度累积值的归一化结果.由于研究区内实测地温数据有限，本文使用仅有的3个站点数据(位置见图 1)参与模型解算.另外，考虑到每个站点的数据都存在不同程度的缺失，我们对3个站点的数据取平均值作为最后的温度数据.

将27个InSAR形变观测结果代入上述模型使用最小二乘求解出研究区季节性形变幅度B.InSAR观测值的中误差由InSAR观测值和模型拟合值决定，可以表达为：

(12)

这里e是由27景影像解得的InSAR观测值与通过最小二乘确定的模型拟合值的差值，N为观测值个数.根据误差传播定律可知，季节性形变幅度中误差与观测值中误差相同，年均形变速率中误差为观测值中误差的0.5倍.计算结果显示，研究区内大部分区域季节性形变幅度中误差小于3 mm；季节性形变较大区域的季节性形变幅度中误差达到5~7 mm，考虑到这些区域土壤含水量大，相干性较其他区域低.总体来说，本文的InSAR处理结果具有较高的信噪比，可以满足冻土形变测量的要求.另外，由于研究区域冻土的真实形变很难与模型完全相同，因此本文所用方法评定的中误差可能偏大.

分别以r₁、r₂、r₃和r₄为参考点解算得到4个形变时间序列(图 6).其中以r₁为参考点解算的形变幅度为1.9 cm，形变速率为-1.1 cm·a^-1；以r₂为参考点解算的形变幅度为2.2 cm，形变速率为-1.1 cm·a^-1，以r₃为参考点解算的形变幅度为1.7 cm，形变速率为-0.1 cm·a^-1，以r₄为参考点解算的形变幅度为1.6 cm，形变速率为-0.1 cm·a^-1.可以明显看出，当r₃、r₄为参考点时和r₁、r₂为参考点时年均形变速率差一个量级，参考点选取不当会使监测到的形变结果大小发生变化，从而导致对年均形变速率的错误估计.

图 6

Fig. 6

图 6 不同参考点下解算的形变序列 Fig. 6 Time series of deformation under different reference points

图 7为基于4个不同参考点(r₁、r₂、r₃、r₄)解得的形变序列与温度序列的相关系数.结果表明，图 7(r₁、r₂)中a区域形变与温度为强负相关，b区域形变与温度为弱相关；通过实地踏勘和土壤湿度结果(图 4)发现研究区域内，a区域广泛发育着冻胀草丘且土壤为亚黏土对应较大的季节性形变幅度，b区域地表干燥有稀疏植被对应形变很小或无形变与图 7(r₁)中的结果十分吻合.图 7(r₃、r₄)中a区域形变与温度为负相关，b区域形变与温度为强正相关，而且其他一些相对土壤湿度小于0.4的区域(图 4)形变与温度也呈现出强正相关.这是因为，r₃和r₄本身可能具有较大形变，以其为参考点解算的结果使得形变较小的区域出现违反冻融机理的形变结果，从而和温度序列呈现强正相关.但较大形变区域的形变序列依然可以和温度保持负相关.以上分析说明冻土区InSAR解算参考点选择不当会影响年均速率和形变空间分布.通过分析4个不同参考点的形变序列与温度的相关系数的空间分布并结合图 5的结果，发现r₁参考点较r₂、r₃和r₄更为可靠，证明了本文使用的冻土区InSAR解算参考点选择策略可以很好地选择出合格的参考点.另外，研究区域地壳形变活跃，增加了冻土形变研究的难度.根据现有研究资料可知，研究区域位于两逆断层之间，水平形变约为6 mm·a^-1，垂直形变小于3 mm·a^-1，形变无明显空间差异(Zhao Q et al., 2016; Li et al., 2018).本文通过选择合适的参考点以及差分计算很好地去除了地壳形变信号(背景信号).在后续大区域的研究中，青藏高原活跃的地壳形变会严重地干扰冻土形变信号，因此需要选取多个分布均匀的基岩点作为参考点联合解算，方能较好地去除地壳形变信号(Daout et al., 2017).

图 7

Fig. 7

图 7 不同参考点下形变序列与温度序列相关系数 (图r1、r2、r3、r4是以图 4中r1、r2、r3、r4为参考点的结果) Fig. 7 Correlation coefficients between the series of deformation and temperature under different reference points (subgraphs r1, r2, r3 and r4 correspond to reference points r1, r2, r3 and r4, respectively)

利用DSs-SBAS方法和27景Sentinel-1A数据，获取了研究区跨时2年的形变时间序列和年均形变速率，主要由位于活动层下的多年冻土上界融化导致.年均形变速率结果(图 8)表明，研究区域大部分处于稳定状态(-1.0~+1.0 cm·a^-1)，形变主要集中在土壤性质为亚黏土的野牛沟沟底区域，同时南坡有零星分布.由于InSAR观测结果对视线向形变十分敏感，南坡较为陡峭的土坡边缘显示出较大(-3.0~-2.5 cm·a^-1)的年均形变速率.南坡靠近沟底部分(区域D及其周边)为一个缓坡，进入融沉期时冻土中孔隙冰融化流失，活动层不断增厚，导致年均形变速率达-2.5~-1.5 cm·a^-1.总体来看，研究区域年均形变速率主要受地形和土壤性质的影响，形变集中的沟底土壤性质为亚黏土(图 1)，土壤湿度较大且该区域地形便于排水固结.

图 8

Fig. 8

图 8 年平均形变速率. A、B、C和D矩形框为4个典型区域的位置和范围 Fig. 8 Annual deformation rates of the study area. Rectangles A, B, C and D indicate typical areas

根据公式(10)和(11)估算得到研究区季节性形变幅度(图 9)，主要由活动层的冻结隆起和融化下沉交替所致.从图 9中可以看出，季节性形变幅度的大小主要受土壤性质和土壤湿度的影响，较大的季节性形变值基本分布在含水量丰富(图 4)且为亚黏土土质(图 1)的野牛沟河两侧.黏土和亚黏土相对砂质土含有更多的水分，在冻融循环中能产生更大的地表形变.另外，其大小和空间分布与高程呈现负相关，随着高程不断增加，年均温度降低、活动层厚度变小、形变幅度也相应变小.研究区域内北坡季节性形变幅度较南坡小，北坡形变集中在靠近沟底的一块洼地(区域B)中，季节性形变幅度为3~5 cm；南坡形变分散(区域C和区域D及其周边区域)季节性形变幅度在2~3 cm之间，其中一些分布在陡峭的土坡边缘.形成这样分布的原因可能是北坡坡度大且较南坡干燥(图 4)，南坡地势平缓且远处有大量雪山融水补给.特别需要指出的是，研究区域东部大部分区域季节性形变幅度均小于2 cm，而东部沟底(区域A)的季节冻融形变幅度集中在3~6 cm.这是因为，东部整体海拔较低，多年冻土分布稀少；而沟底常年有大量积水，土壤含水量大，因此季节冻融形变幅度也最大.

图 9

Fig. 9

图 9 季节性形变幅度. A、B、C和D矩形框为4个典型区域的位置和范围 Fig. 9 Magnitude of averaged seasonal surface subsidence for the study area based on InSAR measurements and Stefan model. Rectangles A, B, C and D indicate typical areas

本文根据野牛沟区域冻土地貌特征以及实地踏勘的结果，选取了4个典型区域(图 10)，计算每个典型区域的年均形变速率和季节性形变幅度，并分析其形变时间序列(图 11).结果表明，多年冻土活动层冻结/融化始日和冻土形变存在明显的空间差异，主要与冻土地貌和土壤类型相关.其中，(1)区域A位于野牛沟河下游河滩，海拔约为3550 m，由于河流搬运作用将上游风化形成的石屑带到下游并分布在河滩上，土壤为黏土且水分补给充足.该区域约从每年8月中下旬开始冻结，次年3月上旬开始融化，季节性形变幅度达3.34 cm，年均形变速率为-0.36 cm·a^-1(图 11a)；(2)区域B位于野牛沟河中游靠近沟底区域，海拔3700 m，土壤为亚黏土，二尕公路贯穿其中.区域内普遍发育着直径约2 m的小热融坑，在公路两旁已经发育成小型热融湖塘，表明此区域地下富含冰层，含水量充足且处于退化状态.该区域约从每年10月初开始冻结，次年3月上旬开始融化，季节性形变幅度达3.20 cm，年均形变速率为-0.67 cm·a^-1(图 11b)；(3)区域C位于PT3观测点处，海拔3850 m，土壤为砂质土，覆盖着低矮的高寒草甸，发育着少量的冻胀草丘且在融沉期基本无地表水.该区域约从每年9月中旬开始冻结，次年4月中旬开始融化，季节性形变幅度为1.89 cm，年均形变速率为-0.27 cm·a^-1(图 11c)；(4)区域D位于野牛沟南坡靠近沟底区域，坡度大且土壤为亚黏土，海拔3730 m，广泛发育着冻胀草丘，融沉期可见草丘间有少量地表积水，坡上发育着数条流水沟.该区域约从每年10月初开始冻结，次年3月上旬开始融化，季节性形变幅度达2.50 cm，年均形变速率为-1.30 cm·a^-1(图 11d)，处于显著退化状态.

图 10

Fig. 10

图 10 典型冻土区实地照片 (a)区域A河滩; (b)区域B热融湖; (c)区域C高寒草甸; (d)区域D沼泽化草甸.典型冻土区位置标注于图 8、图 9. Fig. 10 Photographs of typical frozen-earth areas (a) River shoal in area A; (b) Thermal melt lake in area B; (c) Meadow on highland in area C; (d) Swamping meadow in area D. Their locations are shown in Figs. 8 and 9.

图 11

Fig. 11

图 11 典型冻土区形变时间序列图, 蓝色虚线为冻融期分界 a)区域A; (b)区域B; (c)区域C; (d)区域D. Fig. 11 Time series of deformation of typical frozen-earth areas. Blue dashed line is the boundary between freezing and thawing periods (a) Area A; (b) Area B; (c) Area C; (d) Area D.

由图 11可见，使用Stefan模型可以准确地估算出冻土季节性形变幅度，其中偏离模型的值可能是因为噪声的影响或真实形变本身与模型有着微弱的差异.研究区融沉期基本开始于3月上旬至4月中旬，冻胀期开始于8月中下旬至10月初，冻融开始时刻均略微早于土壤温度的零度时刻且该时间差在不同区域又有所差异.这种时间差的产生可能是因为所用温度数据有限难以代表每个区域并且在土壤温度接近零度时，地表土壤含水量迅速变化，影响雷达信号的传播路径，也会造成这种提前“冻胀”和“融沉”的现象(Zwieback et al., 2015).另外，活动层的冻结融化过程不仅受温度的制约, 还受地形地貌、岩性、地表植被特征以及水文状况等综合因素的影响.因此，冻土的冻融过程具有很大的时空分布差异，不同地区活动层冻融时刻和冻融天数具有很大的差异(赵林等，2000).冻土季节性形变周期也受高程和地下水位变化的影响，随着高程不断增加地下水位相应降低，冻土融化始日和冻结始日也发生了相应的推迟(Daout et al., 2017).区域A、B、D在3月上旬开始融化，而区域C的融化始日要比它们晚15~30天左右，这是因为区域C的地势较高，地下水水位处于较低位置，对温度变化反应相对迟缓.区域A于8月中下旬开始冻结，而区域B、C和D于9月中旬至10月初开始冻结，较区域A推迟30~45天左右，可能是因为区域A富含黏土，地势较低导致含水量充足，地下水水位高且为季节性冻土，当温度开始降低至零度时，表层水迅速冻结发生抬升.

此外，图 11结果表明，4个典型区域的冻土形变时间序列存在季节性周期形变(a，c)和季节性波动下沉(b，d)两种特征.这两种时间上的形变特征广泛存在研究区域内.冻土地貌和水热条件的不同以及人类活动的影响可能是导致形变特征差异的原因，需要今后进一步收集该区域活动层厚度、土壤含水量和地温等信息开展深入研究.

本文提出一种基于分布式目标的小基线集时序InSAR(DSs-SBAS)的冻土形变监测方法，结合分布式目标提取和特征值分解算法提高了冻土活动层形变的空间分辨率和可靠性，并利用地温-形变约束关系提出了一种冻土区InSAR解算参考点选取的新策略，较好地解决了冻土区稳定参考点选取问题.利用2014—2016年27景Sentinel-1 SAR数据对黑河西支源头野牛沟区域的冻土进行形变监测，获取了该冻土区形变时间序列和年均形变速率，并采用Stefan模型估算季节性形变幅度.结果表明：(1)研究区内大部分多年冻土处于稳定状态(-1.0~+1.0 cm·a^-1)，在地形陡峭的南坡边缘及含水量丰富的野牛沟河上游两侧沟底部分区域存在较大形变(区域B、C、D)，年均形变速率最大为-3.0 cm·a^-1；区域内较大的季节性形变幅度主要集中在中下游土质为亚黏土的沟底区域(区域A、B、D)，最大可达6.0 cm.(2)区域内冻土形变时间序列主要呈现以一年为周期的变化，约从每年3月上旬左右开始融化，到9月底达到最大融深；从10月初开始冻胀，到1月中旬达到最大值并持续到3月上旬.受不同冻土地貌、水热条件以及人类活动的影响，存在季节性周期形变和季节性波动下沉两种形变特征.(3)不同区域的活动层冻结/融化始日和冻土形变存在明显空间差异，主要与冻土地貌、土壤类型以及活动层厚度密切相关.