1 引言

多年冻土是指持续三年或三年以上的冻结不融的土层，分为上部活动层与下部冻结层两部分，其中活动层冻融过程可以分为四个阶段：夏季融化过程、秋季冻结过程、冬季降温过程、夏季升温过程^[1].当冻土融化，土壤中水分逐渐由固态变为液态，地表表现为下沉，而当冻土冻结，土壤中水分由液态变为固态，地表表现为上胀.冻胀和融沉这两个过程的反复交替会造成地质环境的破坏，进而引发地基破裂及塌陷、地裂缝、滑坡及泥石流等地质灾害，对冻土区建筑的稳定性和人民的生命财产造成非常严重的威胁.青藏高原位于我国西南部，平均海拔在4000 m以上，其多年冻土占总面积的2/3，是目前世界上海拔最高、面积最大的多年冻土区.高原上自然环境极其复杂，主要土地类型为灌丛、湖泊、裸岩、积雪等，地形起伏大，海拔高，气压低，工作环境极其恶劣且危险，非常不利于进行人工作业.青藏铁路自格尔木至拉萨段位于其中，北起青海省格尔木市，经安多、那曲、当雄、羊八井至拉萨，全长1100余公里，有长达560余公里处于多年冻土区^[2].因此，为保证青藏铁路的正常运营及青藏高原居民区的经济社会可持续发展，对青藏铁路所处的冻土区进行形变监测具有非常重要的意义.就目前而言，对青藏高原地区进行形变监测的方法主要有水准测量、GPS等，然而这些方法无法在大面积无人区进行全天候无障碍监测，也难以获取研究区域内长时间的形变情况.

合成孔径雷达差分干涉测量技术（Differential Interferometric Synthetic Aperture Radar，D-InSAR）是近年来兴起的一种新型的地表形变监测手段，具有全天候、无接触、大面积（100km× 100km）、高空间分辨率（20m×20 m）、高精度（cm至mm级）等特点，相对于传统的形变监测方法具有不可代替的优越性^[3-4].2004年，李震等采用该技术对青藏公路沿线进行形变监测，所得结果与实地监测控制点形变数据非常吻合^[5]，2008年，谢酬等采用ALOS PALSAR数据利用D-InSAR技术监测青藏高原冻土形变，所得结果与水准测量结果有较好的一致性^[6]，证明了D-InSAR能对冻土冻融变化引起的地表形变进行精确的监测.但是，D-InSAR技术的监测精度受到时空失相关和大气延迟的影响^[7-8]，且只能得到单次形变结果，无法获取研究区域地表形变的时间演化情况.为了克服这一问题，近年来国际上一些学者在D-InSAR的基础上提出了时序InSAR技术（MultiTemporalInSAR），通过对大量数据中高相干点目标的时序分析，获取研究地区的地表形变速率和时间形变序列.其中，由Berardino等于2002年提出的小基线集技术（Small BAseline Subset，SBAS），选取空间基线、时间基线都较小的InSAR干涉对，能够有效地减缓由时空基线过长而引起的失相干情况，并且通过最小二乘法则提高了形变监测的时间分辨率^[9]，之后不断有学者将此技术应用于监测城市地表沉降、火山运动等方面，监测结果与实测结果非常符合^[10-11].

可见，目前在冻土监测方面，传统方法虽然时间联系性好，但只能对若干有限点的形变进行监测，不能从面维上监测大范围冻土区整体形变情况；新技术D-InSAR方法能从面维上研究冻土区形变，但是受到时间基线和空间基线的影响而不能获取形变时间演化情况.因此，本文提出将小基线集（SBAS）方法应用于监测冻土形变时间演化中.然而，传统的SBAS方法所采用的线性形变模型往往难以符合季节性变化的冻土形变规律，容易导致较为严重的地形残差和大气延迟误差的残留.基于此，本文提出利用周期形变模型来代替传统的线性形变模型，从而使得SBAS技术更加适用于冻土形变监测.

实验中，选用了ENVISAT卫星获取的21景覆盖青藏高原从羊八井站至当雄站铁路段冻土区的ASAR数据，首次利用SBAS-InSAR技术对冻土进行形变监测，并对所得到的研究区域冻土形变结果进行分析.鉴于目前对青藏高原冻土的研究大都局限于单次线性形变分析或者若干地物点时序形变分析，本文将重点研究青藏高原大面积冻土的季节性形变规律及其同温度变化的联系.

2 研究区概况和数据源 2.1 研究区概况

本文的研究区域位于青藏高原西北部，当雄县境内，区域中心经纬度为（北纬30°14'，东经90°56'），面积约1800km².图 1a显示的是青藏高原地区的地形图，其中的黑色虚线矩形框显示了研究区域的所在位置，图 1b是研究区域地形三维展示.该地区包括羊八井至当雄站约50km青藏铁路段，位于青藏高原念青唐古拉山南麓谷底，青藏公路及部分拉萨河穿越其中，铁路段海拔高度在4300~4640 m之间，铁路路基高2m，底宽10m.区域内以山体为主，地形起伏较大，铁路段沿线处于拉萨地块中部的当雄-羊八井地堑，经过宁中盆地和羊八井盆地，地形较为开阔平坦，西北侧为念青唐古拉山，平均海拔高度达6000m，东南侧为旁多山地，山顶海拔高度达5200m，地堑内部地势低洼，低阶地及河漫滩湿地发育，羊八井盆地与宁中盆地内发育大量中小型沼泽湿地^[12].所选区域自然区上划为青藏高原草原地带，界内主要以灌丛草原生态系统为主，覆盖度达10%~30%^[13].该地区气候属于高原亚干旱气候，冬寒夏冷，日温差非常大，根据拉萨气象站记录，区域内平均年降水量为407~468 mm，平均年温度已从2000年的8.4℃逐渐升高到2009年得10.3 ℃^[14].研究区域内含大面积多年连续冻土带^[15]，冻土活动层会发生季节性融沉和冻胀，引发冻胀丘、热融洼地等不良冻土现象.

2.2 实验数据

本文选用由欧洲空间局（ESA）发射的ENVISAT卫星所获取的覆盖研究区域的21景ASAR单视复数影像（single look complex，SLC）作为实验数据.表 1显示的是所用SAR影像的基本参数情况，其中时间跨度为2007年2月20日到2010年9月7日，影像数据的极化方式为VV，轨道方向为降轨，所处波段为C波段，影像中心入射角约为23°，方位向和距离向分辨率分别为4.053 m和7.083 m.为了提高ASAR影像的轨道数据精度，本研究采用了欧洲空间局发布的DORIS精轨数据.此外，本文实验所用数据还有从美国宇航局（NASA）官方网站下载的覆盖实验区域的SRTM数据，分辨率为90 m×90 m，将用作外部DEM数据以去除干涉图的地形相位.

3 小基线集（SBAS）技术

小基线集技术的基本原理就是将单次D-InSAR得到的形变结果作为观测值，再基于最小二乘法则来获取高精度的形变时间序列^{[9, 16-17]}.为了抑制由于时空基线过长所引起的时空去相干现象，小基线集技术从获取的SAR影像中选择时空基线均小于一定阈值的干涉影像对生成多视差分干涉图，这些干涉影像对会根据基线情况分成若干个集合，每个小集合的形变时间序列可以用最小二乘法进行解算求得，而利用奇异值分解法（SVD）将多个小基线集联合起来进行求解，得到整个时间段的形变时间序列.

首先假设获取了覆盖同一区域N+1幅SAR影像图，（t₀，…，t_N）代表影像的获取时间，并假设所有SAR影像已配准到同一坐标系下，从而可以得到M对时空基线均小于某一阈值的多视差分干涉对，且M满足以下条件：

(1)

随后对生成的多视差分干涉图进行滤波及解缠处理，并利用相干性信息在解缠后的差分干涉图中选取高相干点，然后对所选取的每个高相干点利用接下来所描述的算法模型进行计算.假设第i幅已解缠的多视差分干涉图是由时间t_A和时间t_B获取的两幅SAR影像进行干涉得到的，那么这幅差分干涉图中高相干点的相位值可以表示为^[9]：

(2)

式中，（x，r）为高相干点的影像坐标，ф（t_B，x，r）和ф（t_A，x，r）为两幅SAR影像上的相位值，ф_d（t_B，x，r）和ф_d（t_A，x，r）分别是时间t_B和t_A相对参考时间t₀在雷达视线方向（LOS）上的累积形变相位量，一般假设参考时间的形变相位值为0，即ф_d（t₀，x，r）=0.B_⊥为垂直基线长，Δz为高程误差，ф_atm（t_B，x，r）和ф_atm（t_A，x，r）是两幅影像中大气延迟引起的相位部分，Δn_i（x，r）是由失相关引起的噪声相位，r和θ分别为视线向斜距和入射角，λ为雷达波长.

为了获取准确的形变时间序列，需要首先去除地形误差相位及大气相位.对于本文研究所针对的冻土而言，地表形变的低频部分可以假设由线性沉降和季节周期形变两部分组成^[18]，可以由下式表示^[19-20]：

(3)

式中，t为影像获取时间相对于所选参考影像累积时间量；T为季节性形变周期，一般为1年；α₁，α₂，α₃为模型参数，则（2）式可表示如下：

(4)

式中Δф_i，res（x，r）为第i幅干涉图中的残余相位.对（4）式进行转换，所有干涉图相位可以表示为如下矩阵：

(5)

其中P^T=[α₁，α₂，α₃，Δz]为待求参数组成的向量矩阵；δф为由所有差分干涉图中相位值所组成的向量矩阵；A为系数矩阵，比如若有δф₁=ф₄ -ф₂及δφ₂=φ₃ -φ₀，则：

(6)

（5）式中有M个观测量和4个未知数，若M≥4，就可以利用最小二乘法来进行求解：

(7)

由式（7）计算可以得到高程误差Δz和残差部分Δф_res，残差部分中包括高频形变、大气延迟和噪声影响，根据大气延迟误差和形变以及噪声的不同时空特性，可以将残余相位中的大气延迟误差分离出来.一般大气延迟相位在空间域上表现为低频信号，而在时间域上是一个随机过程，没有任何相关性，为高频信号；形变相位在空间域和时间域上都表现为低频信号；而噪声相位在空间域和时间域上都表现为高频信号.因此可以通过时间高通滤波和空间低通滤波来得到大气延迟引起的相位^[21-22]，本文采用三角滤波法进行时间滤波，具体计算过程如下：

首先将差分干涉图中残差作为观测值采用SVD法获取SAR影像获取时间相位图像序列，得到的图像序列m（x，r，t）可以认为是由非线性形变相位def（x，r，t）和大气相位及噪声atm_noi（x，r，t）叠加而成：

(8)

式中t=（t₀，…，t_N）为SAR影像获取时间序列，（x，r）为像元位置.因为大气和噪声相位atm_noi（x，r，t）在时间上是不相关的，根据时间不相关特点选用式（9）来求取像元大气延迟及噪声相位之和的时间序列，

(9)

逐点计算之后得到所有高相干点包括大气延迟及噪声之和的相位时间序列.再通过均值法进行空间滤波，滤波窗口为5×5，从而得到高相干点大气延迟时间序列，再进行差分获取对应差分干涉图中大气延迟相位值.

通过上述计算处理过程获取了高程误差和大气延迟，将它们从差分干涉相位图中减去，得到只包括形变及噪声的差分干涉图，从而可以将该相位值作为观测值求取时间序列形变值：

(10)

其中，δφ'_i（x，r）为去除了高程误差和大气延迟的相位；Δn'_i（x，r）为噪声相位以及一些地形和大气的残余相位.同样，所有干涉图相位可以表示为如下矩阵：

(11)

式中，系数矩阵B[M ×N]每一行对应一个干涉图，每一列对应一景SAR影像图^[9]，式中有M个方程，N个未知量，若M≥N，并且A的秩为N，则可利用最小二乘来求解每幅影像形变相位值：

(12)

由于M个干涉影像对可能由于时空基线的分布情况不处于同一个小基线集组合，矩阵B的秩小于N，那么相应的法方程系数矩阵B^TB秩亏，因此根据最小二乘法所得到的结果并不唯一.为了解决法方程秩亏问题，可以采用奇异值分解法（Singular ValueDecomposition，SVD）来进行求解^[9].

如果利用SVD法直接对相位进行求解，解算得到的相位值在时间上会出现很大的跳跃性，不符合形变的物理规律^[9]，那么为了获取具有物理意义的形变序列，可将两个获取时间之间的平均相位速度作为未知量来进行求解

(13)

则式（10）可以表示为：

(14)

代入（11）式可得：

(15)

式中，B[M ×N]为系数矩阵，对矩阵B使用SVD法得到速度矢量v的广义逆解：

(16)

其中U为M行M列的正交矩阵，其列向量为BB^T的特征向量；V为N行M列正交矩阵，其列向量为B^TB的特征向量；S为M行M列对角矩阵，对角元素为B的特征值σ，当小基线集个数为L时，B的秩为N-L+1，S中有N-L+1个空特征值，那么S可以表示如下：

(17)

求出v之后，根据每个时间区间形变速率对各时间段的速度在时间域上积分得到每个时间段的形变相位值，最后通过乘以相位形变转换系数（λ/4π）得到整个监测时间内地表在雷达视线方向（Line-Of-Sight，LOS）上的形变时间序列.整个SBAS方法的技术流程图如图 2所示，本文主要对流程图中基于高相干点建立观测方程这一步骤进行了改进，用（3）式所示的周期形变模型来代替传统SBAS过程中所用的线性模型^[9].

4 青藏高原冻土实例研究

本文选取青藏高原冻土进行实例研究，实验具体区域及实验数据本文已在前面进行了详细介绍，下面将对实验流程主要处理环节进行描述，并对所得到的青藏高原冻土形变结果进行分析.

（1）实验的首要步骤是选取符合时空基线条件的干涉对.本文研究采用的空间基线阈值为150m，时间基线阈值为2年，最后共得到65对干涉对.如图 3所示，干涉对被分为两个小基线集.随后将所有的干涉对配准到同一坐标系上，便于后续的计算.

（2）本文选取二轨差分干涉测量方法对选取的干涉对进行处理，并采用SRTM数据作为外部DEM来消除地形相位信息^[23-24].为了减少计算量并抑制噪声，对干涉图采取距离向2视、方位向10视的多视操作.随后利用基于最小二乘的滤波方法进一步消除多视差分干涉图中噪声的影响^[25]，并利用最小费用流方法进行相位解缠^[26-27].

（3）虽然实验中已经采用了精密的轨道数据进行定轨，但是在解缠后的差分干涉图中仍然存在一些由轨道误差引起的趋势相位，所以本文采用了迭代最小二乘法对其进行拟合：

(18)

其中x，y为SAR坐标系下像素在距离向和方位向的坐标值；b₁~b₅为待定模型参数；R_orbit为在差分干涉图中均匀提取的相位值，本文所用的采样间隔为20×20，并在选取过程中舍去了位于低相干区域和形变区域的像素点.利用迭代最小二乘法对（18）式进行解算，基本上进行两次迭代就可以得到可靠的结果.最后利用所求得的模型参数计算出所有像素由轨道误差引起的趋势相位，并从差分干涉图中将其去除.

（4）为了避免低相干点带来的计算误差，在差分干涉图中选取高相干点进行计算.本文利用所有SAR影像平均相干值和最小相干值来选取目标，即选取的高相干点的平均相干值大于0.5，且在所有差分干涉图的相干值的最小值大于0.3（小于此值的差分图中相应像素点被视为失相干）.图 4显示的就是研究区域的平均相干图，可以看出，相干性较好的地区为地形较为平坦的地堑区和盆地，两旁为念青唐古拉山和旁多山区，相干性非常低，铁路沿线的相干性也相对较低.

（5）对选取的高相干点逐点利用所描述SBAS方法进行求解，利用（4）-（9）式分离出高程残差和大气延迟相位，再结合（10）-（17）式获取地表在LOS向的形变时间序列结果.由于ENVISAT卫星的雷达入射角较小（约23°），因此LOS向的形变结果基本上反映的就是地表的抬升或下沉.图 5为实验所获取的高程残差图（图内黑色实线代表青藏铁路所在位置），从图中结果可以看出实验数据高程误差绝对值大面积介于0~6 m的范围，局部地区达到了8~10m的误差，与SRTM数据10m的相对高程精度非常吻合^{[23, 28]}.而10 m高程误差在垂直基线为100m的情况下引起差分干涉结果中形变误差基本都在1cm之内.

图 6a为研究区域的季节性形变幅度图，选取的形变参考点为念青唐古拉山山坡上裸露于地面的岩石（图中白色方框所示），假设该点稳定，其形变值为0.可以看出研究区域整体上季节性形变较为明显，且东南方向形变幅度相对较大，并且从山顶至山底方向形变程度越来越大，区域内形变幅度最大可达到2cm，部分出现在铁路沿线.研究区域的冻土形变除了有明显的季节性变化特征之外，还存在着长期逐年沉降的趋势，图 6b为实验得到的平均形变速率图，从图中可以看出青藏铁路附近区域存在每年0~1mm的沉降，且从西南至东北方向有逐渐增大的趋势，越接近山脊处，沉降量越大.这与山脊处积雪量较多有关，随着温度的升高，积雪不断融化，积雪表面高程降低，同时由于积雪覆盖，地表下层热量积累，导致冻土逐年融化下沉，因此地表下沉量更大.

图 7显示的是2007年5月1日至2010年9月7日研究区域的冻土形变时间序列图.各形变图显示的均为以2007年6月5日为参考时间的累积形变量，每幅图左上角所标的数值为SAR影像的获取时间（如20070501代表 2007年5月1日）.从图中可以看出，研究区域的地表形变存在明显的以年为周期的季节性变化规律，即表现为融沉与冻胀过程的不断交替，且在7-8月份间，由融沉引起的地表下沉量达到最大，达到了-29mm，在1-2月份间，由冻胀引起的地表抬升量达到最大值，可达到34mm.此外，从图中还发现从念青唐古拉山山坡至山底形变幅度越来越大，这是因为山底部分接近拉萨河，水分补给充足，土壤湿度含量较高，冻土发育比较充分，活动层相对较厚，冻土形变量相对也要大一些^[18]，除此之外，山底接近铁路及公路沿线，青藏铁路和青藏公路的建成，改变了冻土与外界大气间的热交换平衡状态^[29]，由于路面强烈的吸热作用导致过剩的热量积累^[30]，使得山底部分出现相对较大的融沉变形.同时还可以看出从东北至西南走向形变程度有逐渐减小的趋势，这是因为东北部为宁中盆地，盆地内发育大量中小型沼泽湿地，相对而言盆地地质更为湿润，冻土发育更充分，活动层也更厚^[12]，所以形变量更大.

为了对研究区域进一步分析，本文分别选取研究区域内四个典型地物点（具体位置如图 6中的红色五角星所示），获取它们的时间形变规律.图 8所示为所选取的地物点的时间形变图，（a）、（b）、（c）、（d）分别对应图 6中a、b、c、d点的时间形变序列图，形变值均以2007年6月5日作为参考，图中位于上方的虚线为拉萨市月平均温度拟合曲线，位于下方的虚线为各地物点形变值拟合曲线.a、b两点分别位于当雄盆地和羊八井盆地，且都远离青藏铁路及公路沿线等人工建筑活动区，接近山脊处，从图 8（a、b）中可以看出a、b点所在地理位置地表形变以周期形变为主，其中周期形变以年为周期，形变幅度达到15mm，伴随着周期形变的同时每年还存在2~3mm的线性沉降.c点位于宁中盆地、d点位于宁中盆地中心铁路沿线，从图中可以看出，c点和d点同样是以周期形变占主导，同时伴随着一个长期的沉降趋势，形变幅度达到25 mm，每年同时间发生了0~1mm的线性沉降.可以看出越接近盆地低洼处的区域冻土周期形变的幅度越大，因为盆地中心发育着大小不等的沼泽，水含量较为丰富，冻土活动层较厚^[12].同时从图中还可以看出冻土形变同冻土区温度变化存在明显的联系，即所选地物点形变情况与温度成一定的反向变化，在春夏季，随着温度升高冻土融化，地表下沉，在8月份左右融沉量达到最大，进入秋冬季，随着温度降低，冻土冻胀，地表上升，在2月份左右冻胀量达到最大.由此可见，冻土区域的形变受到冻土的气温、地质和地貌等的综合影响.

5 结论

近年来，不少专家学者一直在探讨对青藏铁路沿线区域冻土形变进行监测的方法，但是目前能大面积对冻土形变规律进行监测的手段还很局限.在此基础上，本文对InSAR技术中的小基线集（SBAS）方法进行改进，利用周期形变模型来代替传统的线性形变模型，从而使其可以更好地应用于冻土形变监测中.利用ENVISAT卫星在青藏高原羊八井至当雄青藏铁路段区域获取的3年共21景ASAR影像作为实验数据，改进的算法模型很好地除了地形残差及大气延迟引起的误差，最终获取的冻土形变时间序列结果同其物理变化规律非常符合，证明了利用SBAS-InSAR技术来获取冻土点、线、面多维形变趋势及规律具有良好的应用前景.

从实验结果可以看出，羊八井至当雄青藏铁路段区域的冻土形变主要呈现周期性变化，以1年为周期，其变化规律同温度变化存在联系，即当温度升高时冻土融化，地表形变现象表现为下沉，当温度不断降低时冻土冻胀，地表形变表现为上升，且从西南往东北走向形变幅度逐渐增大.除此之外，随着全球气温的逐年不断持续上升以及人类工程活动的影响，使得部分冻土永久消融^[31-32]，导致研究区域每年同期还有一个毫米级的沉降量.此次监测结果成功地识别出来羊八井至当雄青藏铁路段区域的冻土形变的时空分布规律，选择在形变较小的地方进行人工建设活动，避开形变很大、超出可承受范围的地区，可以减少冻土形变对建筑物及人们生活的影响.此外，研究结果还发现，铁路沿线区域也存在一定程度的形变，因此为了保证青藏铁路正常运营，必须在冻土区铁路段采取一些保护措施，如通风路堤、热桩、片石器冷路基和热棒路基等，以减少热量进入路基同时加强散热效果来达到保护冻土的效果^[32-33].

致谢

感谢欧空局提供的ENVISAT ASAR数据（AO-4458，4914，9235）.