1. 引言
冻土根据其存在时间及变化情况可以划分为短时冻土、季节冻土以及多年冻土。冻土作为自然环境的指示标,对于人类生存环境的变化具有很重大的意义。虽然永久冻土是一种地下现象,但其状态的变化会影响地表特征和工程基础设 [1]。多年冻土和季节性冻土在我国分布均非常广泛,分别占我国国土面积的25%和54% [2]。19世纪末,冻土学首次被提出为一门独立的学科,并成立了冻土研究学会。永久冻土区的研究主要集中在北半球或高纬度的一些国家和地区,我国西部青藏高原冻土分布较广泛 [3]。随着我国综合实力不断提升,国家实施一定的经济战略和政策,如西部大开发。但由于西部基础设施落后,环境险峻,建设实施受到多年冻土的影响,通过人工调研很难解决问题,且实施起来效率也不高,合成孔径雷达差分干涉测量(D-InSAR)技术广泛应用到冻土形变监测,山体滑坡,地面沉降等领域,已有学者利用此技术来监测青藏髙原冻土区的形变,但是,青藏高原地形起伏大,若釆用D-InSAR技术容易因空间基线过长而受到地形误差的影响,为解决相关问题和潜在的隐患,大批研究学者投入到我国西部的冻土研究中,包括地质构造、冻土监测、隐患排查等 [4] [5] [6]。为此,Xie [7] 等人利用PS-InSAR技术对北麓河区域多年冻土的情况进行了探测;2010年,Liu [8] 等人对阿拉斯加地区的冻土情况进行研究,获取该区域地表形变结果,并结合冻土物理形变规律进行分析;2016年,Zhao Rong [9] 等人利用SBAS-InSAR技术对青藏高原羊八井地区冻土特征进行探测,获得的形变结果能较好的符合冻土形变规律。但是由于传统SBAS-InSAR技术在冻土形变监测的实际应用中存在计算效率低与相位优化差的问题。本文通过ISCE开源平台添加相位优化及降噪算法进行高精度冻土形变监测的研究,改进的SBAS-InSAR通过干涉图反演进行解缠误差校正及不同的权重方式进行降噪处理,从而能够更加准确的得到形变相位。
研究区域选择盐湖一带(图1)由于该地位于青藏高原永久冻土区,地处中国青海可可西里国家级自然保护区,位于水系交汇地区,冻土区域发育较好,土壤含水量丰厚 [10]。受全球气候变暖及降水量变化的影响,近年来盐湖湖水面积急剧扩大,地表水的增加必然导致冻土发生变化,从而威胁到基础设施且存在地质灾害隐患 [11]。因此,对该区域进行冻土形变监测研究就显得格外重要。通过选取2014年10月22日至2019年8月21日的121景Sentinel雷达数据影像及地形气象等数据进行冻土形变的综合监测与分析,成功获取了青海西部盐湖及其周边冻土区的地表形变时间序列图,揭示了该冻土区从2014年到2019年的冻土引起的地表形变情况。研究内容主要分为三个方面:第一,探究冻土形变规律及影响因素;第二,研究冻土变化对于周边基础设施的影响;第三,发现冻土区域潜在的地质灾害隐患。通过与前人的研究成果进行比较发现所得到的地表形变规律非常吻合,证明了其在冻土形变监测领域具有一定的可行性。
Figure 1. Distribution and study area of permafrost in Qinghai Tibet Plateau
图1. 青藏高原冻土分布与研究区域图
2. 研究区概况及数据源
2.1. 研究区概况
研究区位于高海拔地区,属于高原大陆性气候,冬寒夏凉,雨热同季。既包含永久冻土也包含季节性冻土,区域内广泛发育着大面积的中、高纬度连续多年冻土,具有低温、高厚度的特点,对气温的变化非常敏感。冻土内部构造从上至下依次为季节冻土活动层与多年不融冻土层,地表形变受冻土变化的影响,热融湖、冻胀丘、石海、石河、斑土、冻胀草环等特征地貌形成。
2.2. 数据源
2.2.1. Sentinel-1数据
由于该地区植被较少,积雪相对较薄,相干性较好,很适合通过干涉合成孔径雷达(InSAR)去分析冻土的形变过程。Sentinel-1能够在较短的时间内获取长时间跨度的数据,这是长期挖掘缓慢变化的冻土的重要保证,极大地解决了冻土变形监测中雷达图像积累的瓶颈和困难,进一步降低了时间去相关的风险,使得产生具有良好空间覆盖的地表位移模型成为可能,使得其在长时间序列、短时间间隔的冻土变形监测中具有广阔的应用前景 [12]。Sentinel-1是一个全天时、全天候雷达成像系统,重访周期12天。Sentinel-1 基于C波段的成像系统采用4种成像模式来观测,工作模式如表1所示。可精确确定卫星位置和姿态角,可选用的轨道参数文件如表2所示,Sentinel-1具体的产品数据类型见表3。
本文研究的盐湖区域选用干涉宽幅模式的斜距单视复数产品,即通过IW模式获取Level-1的SLC产品来进行形变监测分析。本文选用的轨道文件为AUX_POEORB精确的轨道星历参数,分析得到的结果更加准确,该文件是在数据获取后21天内生产得到的,精度5 cm以内。通过获取121景Sentinel-1升轨影像数据对盐湖一带区域进行冻土形变研究,影像参数见表4所示。
2.2.2. SRTM DEM数据
SRTM数据能获取60˚N~60˚S覆盖全球面积80%以上的三维地表地形信息。SRTM数据包括SRTM1和SRTM3两种,两种数据的空间分辨率分别为30米和90米,本文采用的是30 m空间分辨率的SRTM1数据。SRTM的主要技术参数如表5所示,本文实验区域为青藏高原玉树市治多县附近,地理覆盖范围为(34~38N,92~96E),下载对应的SRTM1数据,共获得25个geotiff格式的DEM文件,影像数据三维显示如图2所示。
Table 5. Main technical parameters of SRTM
表5. SRTM主要技术参数
Figure 2. Three dimensional display of SRTM data in study area
图2. 研究区域SRTM数据三维显示
2.2.3. 轨道数据
利用SBAS-InSAR技术获取形变信息的过程中需要利用轨道姿态参数,由于雷达影像元数据中所包含的轨道姿态参数量较少,精度并不高,需要获取精度更高的轨道姿态量,通过结合外部卫星轨道数据数据来解决,对于Sentinel卫星轨道数据可在EARTHDATA平台上获取,本次研究共获取98个AUX_POEORB (精确轨道星历参数)轨道参数文件。
3. 理论及方法
3.1. 相位解缠误差校正
使用ISCE中的哨兵堆栈处理器来处理干涉图堆栈;将每景SAR图像与它在时间上最近邻的5景影像及进行配对(序列网络),干涉图采用最小成本流方法进行相位解缠。在相位解缠过程中,错误的整数周期(2π rad)会添加到干涉相位中,使其出现解缠误差,对于Sentinel-1数据集处理流程使用桥接和相位闭合方法纠正解缠误差,分别利用空间和时间域中的解缠误差。
3.1.1. 桥接相位
在空间域中,相位偏移会导致像素组出现解缠误差,绝大部分区域具有中到高的空间相干性,并且由于失相关或高形变相位梯度而相互分离。假设邻近的可靠区域之间的相位差异小于半个周期(π rad)。解缠误差校正的任务是要确定添加到每个可靠区域的整数周期相位偏移量,以便在区域之间对齐相位值,利用树搜索算法自动连接可靠区域的桥接方案。这类似于相位解缠中的区域集合。为实现相邻可靠区域间相位梯度平滑的假设,在相位解缠前将对流层、DEM误差、变形模型、坡道的贡献去除,校正后再将其加入。
对于每个干涉图,桥接方法可分为三步。首先是使用相位解缠算法如SNAPHU,使用所连接的组件信息来识别可靠区域。小于预选尺寸的区域将被丢弃。对于每个区域,在形态学图像处理中通过使用预先选择的形状和大小侵蚀,丢弃边界上的像素。第二步是使用最小生成树(MST)算法最小化桥的总长度,构造有向桥来连接所有可靠区域。使用广度优先算法来确定顺序和方向,从最大的可靠区域开始。第三步是估计两个区域之间的整周相位偏移量。为此,首先将相位差估计为以两个桥端点为中心的预选大小窗口内像素中值的差,整周相位偏移量是将相位差降到的整数周数。该算法可以选择基于最大可靠区域来估计线性或二次相位梯度,在偏移估计之前从整个干涉图中移除,然后在校正之后添加回来。
3.1.2. 相位闭合
在时域中,解缠误差可能会破坏三角干涉相位的一致性(s1s2-s2s3-s3s1)。闭合相位是解缠干涉相位的循环乘积 [13],公式(1)、(2)、(3)引用于参考文献 [13] :
(1)
,与是由在,,时刻获取的SAR图像所生成的三个解缠干涉相位,闭合相位的整数模糊度为:
(2)
其中wrap是将输入数值包裹为的数。在没有解缠误差的三角干涉相位中,在所有三联体中非零的三联体的数量为:
(3)
T为三联体的数目,可用于检测解缠误差。
使用自动检测和修正干涉图网络中解缠误差的算法,评估相位解缠并利用闭合相位信息修正解缠误差,在三维相位解缠过程中,使用闭合相位迭代调整成本,将整数周期相位偏移量添加到相位解缠误差的区域,直观地识别和纠正解缠误差 [14]。
对于冗余的干涉图网络,对于每个像素解缠干涉相位的整数模糊度,其时间一致性可以表示如下,公式(4)引用于参考文献 [14] :
(4)
其中C是所有可能的干涉图三联体的T × M设计矩阵,U为满足干涉相位一致性所需周期数的M × 1向量,并联合L1范数正则化最小二乘优化 [15] (也称为最小绝对收缩和选择算子(LASSO))来求解,公式(5)、(6)引用于参考文献 [15] :
(5)
其中为L1和L2范数项之间权衡的非负参数,取值基于不同值的拟合。校正后的解缠干涉相位为:
(6)
其中round是将输入数字四舍五入到最近整数的运算符。
3.2. 降噪
3.2.1. 根据残差相位判断噪声水平
在对于任意SAR干涉图的噪声水平,通过计算残差相位的均方根(RMS) [16] 来判断,公式(7)引用于参考文献 [16] :
(7)
其中,表示像素p的处的残差相位,是在网络反演过程中基于时间相干性选择的可靠像素集,总数量为。在计算RMS之前,对每幅干涉图的残差相位去除二次梯度的影响。在相位分量的噪声评估中,包括残余对流层波动、未校正的电离层波动,以及剩余的失相关噪声。将RMS
值大于预定义阈值干涉图标记为高噪声影像,进行地形残差和形变速率估计时不涉及这些影像,从而降低噪声的影响。
3.2.2. 基于相干性进行网络降噪
失相关噪声是主要的误差源,在通过时空基线阈值排除失相干影像的基础上,还通过使用基于相干性的网络校正来移除低相干性的干涉图,即排除高相干像素百分比较低的干涉图。
3.2.3. 通过相位方差倒数的权重函数抑制噪声
通过采取相应的权重函数来减少随机失相关噪声,给予高相干(低方差)的干涉图比低相干(高方差)的干涉图更大的权重,从而更准确地估计形变测量值和可靠性测量值 [17]。
干涉图由相位方差的倒数来加权,矩阵的形式如下,公式(8)引用于参考文献 [17] :
(8)
其中是通过相位概率分布函数(PDF)积分得到的第j个干涉图的相位方差。
4. 技术路线与数据处理
4.1. 技术路线
为研究盐湖周边冻土形变情况,本文从121景实验数据SAR影像生成475景干涉对影像,通过空间相干性来选取相干性高的干涉图288幅,通过差分干涉与相位优化获得研究区地表形变时间序列,影响分析效果的关键因素在于解缠误差的校正与降噪处理,通过干涉图反演及不同权重方式进行降噪,达到相位优化的效果,以获取更加准确的形变相位,图3是冻土形变分析的流程图。
Figure 3. Technical route of frozen soil deformation monitoring
图3. 冻土形变监测技术路线
4.2. 数据处理
4.2.1. 干涉对处理
根据获得的盐湖周边2014~2019年的Sentinel卫星121景SAR影像,根据时间空间基线、相关性等选择符合研究的影像,图4是时空基线图,用于SBAS方法中的小基线集的选取。
Figure 4. Sentinel-1 image time-space baseline
图4. 哨兵-1影像时空基线图
4.2.2. 干涉图网络校正
网络校正的目的是排除一些质量差的干涉图的影响,来获得更佳的时间序列分析结果。通过选择合适的SAR影像后对其进行差分干涉,生成的干涉图会有由于多种原因,如空间相干性不够,噪声、残差等影响最后形变速率的估计,通过对空间相关性设置一定的阈值,来选择最佳的干涉影像对,从而减小误差,也有利于形变速率的估计,为保证干涉影像一定的时间相干性,将阈值设置为0.65,经该阈值进行网络校正处理后的干涉图对如图5所示,根据空间相关性决定选取或丢弃的干涉图对,用蓝色的实线连通的是选用的干涉对,反之,虚线则表示相关性不够而舍弃的干涉图影像。
4.2.3. 选择参考点
由于研究区域位于高原冻土区,植被覆盖度较低,一部分相干性较差的河流、湖泊在进行掩膜时会被去除而不参与分析,整体相干性较好,而且受大气的影响基本一致,选取影像中受大气影响较小、与感兴趣区域海拔高度接近的高相干性参考点。
4.2.4. 相位解缠误差纠正
对于相位解缠误差的去除,首先基于空间相关性对相关性较低的干涉图进行去除,这时已经去除了几乎所有的解缠误差,还有一部分解缠误差可以通过权重方法来进行抑制,相位解缠误差通过桥接法与相位闭合法进行纠正。
4.2.5. 网络反演
该步骤对原始相位时间序列的干涉图网络进行反演。原始相位时间序列包括由DEM误差引起的地面变形、大气延迟和地形残差的贡献。通过相位方差的倒数对干涉图进行加权。反演出来的原始相位干涉图如图6所示。
4.2.6. 大气流程校正
研究区域使用全球大气模型Global Atmospheric Models (GAMs)对大气延迟误差进行校正,经过大气延迟校正后的盐湖一带区域干涉图时间序列如图7所示。
其他贡献相位校正的顺序可以互换,其运用的是INSAR数据的各个不同方面,如对流层延迟的校正通过去除与地形相关的信号;相位梯度的校正去除与空间坐标系相关的信号;地形残差的校正去除与垂直基线相关的信号等。
4.2.7. 相位梯度的去除
相位梯度是由剩余的对流层和电离层延迟引起的,在较小程度上是由轨道误差引起的。对于每景干涉图可靠像素,在位移时间序列中移除线性或二次梯度。相位梯度处理后的形变情况如图8所示。
Figure 7. Time series after correction of tropospheric delay
图7. 对流层延迟校正后的时间序列
4.2.8. 地形残差(DEM误差)校正
利用DEM与垂直基线时间序列的关系来校正由于DEM误差引起的相位残差,并得到形变时间序列图,本文实验区域为青藏高原玉树市治多县附近,地理覆盖范围为(34~38N,92~96E),下载对应的SRTM1数据,共获得25个geotiff格式的DEM文件。图9是去除地形残差影响后得到的最终的地表形变时间序列图。
Figure 8. Deformation after phase gradient treatment
图8. 相位梯度处理后的形变情况
Figure 9. Time series of deformation after terrain residual removal
图9. 地形残差去除后的形变时间序列
4.2.9. 用相位剩余均方根进行噪声评估
计算时间序列中每幅干涉图的残差相位均方根,然后选择RMS值最小的日期作为最优参考日期。检测RMS超出离群点检测阈值的噪声图像并将其去除。通过设定阈值,选择一副低RMS的影像作为参考,排除其他RMS较高的影像,以减少形变估计时的误差。RMS可以用来表现InSAR时间序列的噪声水平,通过计算残差相位的均方根(RMS)来判断噪声水平,最终选择2018年11~12月份生成的干涉图作为参考影像,能使形变时间序列图可视化效果更好,最终根据残差相位均方根(RMS)所获得的干涉图及参考影像如图10所示。
Figure 10. RMS estimation for reference image selection
图10. RMS估计用于参考影像选取
4.2.10. 形变速度估计
计算位移时间序列最佳拟合直线的斜率及其标准差,并以此估计平均速度。进行估计时不会使用残差相位均方根(RMS)较大的干涉图影像,即噪声较大的干涉图。然后对其进行地理编码,从雷达坐标重定向到地理坐标,最后得到平均形变速率图,图11表示的是盐湖一代区域的地表平均形变速率。
Figure 11. Surface average deformation rate map
图11. 地表平均形变速率图
5. 结果分析
5.1. 研究区域冻土形变情况
根据上图图11可知,研究区域2014年10月22日~2019年8月21日整体出现下沉。区域冻土的“热熔冻胀”现象呈现季节周期性变化,地表累积沉降量逐年增大,在经纬度35.740662˚N,93.073174˚E附近最大年均沉降量可达到−6.59 ± 0.10 [cm/yr],累积最大沉降量已达到−31.81 cm,这是因为在春夏季节,地面温度随着气温升高而升高,致使地表下冻土开始融化,融化深度也随之加深,表现为地表下沉;由于山体的岩石冰渍碎屑不断向下垮塌并堆积,经纬度位置为35.740662˚N,93.521896˚E附近最大年均形变量可达到11.69 ± 0.07 [cm/yr],累积最大形变量已达到55.13 cm,逐步形成岩石冰川地貌。
5.2. 冻土形变规律
5.2.1. 气候因素与冻土形变相关性
获取研究区域温度、降水量数据与形变数据进行联合分析,分析结果如图12所示。图中季节周期性变化存在延迟效应,温度高于0℃一段时间后开始沉降,低于0℃一段时间后开始抬升,主要是受到水的比热容的影响。该区域雨热同季,夏季温度升高时,伴随大量降雨,地表沉降速度加剧,反之,抬升加剧,是由于高温与强降雨共同作用使热量充分传递给冻土,水由于重力产生的渗透作用使其快速与冻土发生热熔,使得形变速率加快,多年冻土形变与年均气温的升高密切相关,相关性用Pearson指数表示如图13所示,相关性可以达到−0.689,呈现显著中等相关,这在一定程度上也是由于受到了季节周期性形变的延迟效应的影响。
Figure 13. Correlation between frozen soil deformation and temperature
图13. 冻土形变与温度相关性
5.2.2. 冻土季节周期性变化
冻土活动层变厚,永久冻土逐年衰减,多年冻土区活动层起始融化主要发生在4月初,起始冻结主要发生在9月中下旬,地表形变规律与冻土冻胀热融的规律相吻合,如图14所示,图中蓝色箭头处表示从该日期开始地表由于冻土季节性变化发生“冻胀”抬升,红色箭头处表示从该日期开始地表由于冻土季节性变化发生“热熔”沉降,且由于冻土活跃层的不断变深累积,季节周期性形变程度愈发显著。
Figure 14. Fitting of frozen soil deformation curve
图14. 冻土形变曲线拟合
5.3. 地质灾害隐患
图15表示岩石冰川地貌特征,是在高山环境中由松散的岩石碎片和冰混合而成的叶状或舌状地貌。图16表示岩石冰川区域形变情况,红色三角形区域表示特征区域,该区域年平均形变速率为11.69 ± 0.07 [cm/yr],累计形变量为55.13 cm,该区域在地貌形成过程中,岩石和冰渍不断下滑堆积。这是由于冻土发生变化导致熔融的岩石和冰向下爬移并黏合,不断堆积的过程中会伴随着掉落垮塌情形,会威胁到下坡地区的基础设施 [18],夏季的雨季和冬季的旱季对位移数据有显著的影响,前一个季节在许多斜坡上引起更明显的地面滑动,而后一个季节阻止了大多数滑动过程 [19],本研究中发现多处掉落垮塌隐患点,谷歌地球影像是由SPOT卫星获得的,空间分辨率小于5 m,因此我们将其叠加到谷歌地图进行显示。
Figure 16. Characteristic regional deformation of rock glacier
图16. 岩石冰川特征区域形变
5.4. 基础设施隐患
5.4.1. 青藏公路、铁路地表沉降东北路段较西南路段大
青藏公路、铁路路段年均形变量大,东北部路段灾害风险大于西南路段,如图17所示,极可能是由于东北路段受盐湖和周围河网的影响较大,使其土壤含水量较西南部高,使其发生沉降,垮塌等地质灾害的风险更高。
由于湖水面积的扩大,随时有可能再次发生溃堤,冲断青藏公路,因此研究青藏公路周边地表形变的变化就显得尤为重要;当青藏公路左侧沉降趋势明显,使得地形海拔降低时,盐湖的湖水对于公路的威胁就会增大。
Figure 17. Deformation of Qinghai Tibet railway and highway in Salt Lake section
图17. 青藏铁路、公路盐湖路段形变情况
5.4.2. 青藏公路、铁路地表形变情况
图18表示公路形变情况,红色三角形是谷歌影像上显示的特征公路区域,黑色的方形点是根据相关性选取的参考点,从离散点分布图中可以得到随时间变化的形变情况,公路年平均沉降速率可达−2.58 ± 0.09 [cm/yr],累计形变量为11.42 cm;图19表示铁路形变情况,从离散点分布图中可以得到随时间变化的形变情况,铁路年平均沉降速率可达−2.10 ± 0.08 [cm/yr],累计形变量为9.4 cm。
6. 总结
1) 盐湖地区冻土形变情况及规律:整体处于沉降趋势,地表累积沉降量逐年增大;冻土的“热熔冻胀”现象呈现季节周期性变化,雨热同季,夏季温度升高时,伴随大量降雨,地表沉降速度加剧,反之,抬升加剧,季节周期性变化存在延迟效应,主要是受到水的比热容的影响,多年冻土形变与气温变化密切相关,本文对冻土形变规律进行了较准确的拟合,之后研究可根据形变规律结合地质类型、土壤湿度等因素构建永久冻土形变模型,对同类地区冻土研究具有一定的指导意义。
2) 研究区基础设施影响情况:青藏公路、铁路路段年均形变量大,东北部路段灾害风险大于西南路段;公路年平均沉降速率可达−2.58 ± 0.09 [cm/yr],累计形变量为11.42 cm;铁路年平均沉降速率可达−2.10 ± 0.08 [cm/yr],累计形变量为9.4 cm,根据形变特征可对盐湖一带基础设施采取相应的保护措施,避免潜在性的风险发生,为其他地区永久冻土变化对基础设施影响的研究提供参考。
3) 研究区存在的地质灾害隐患:由于冻土发生变化,导致岩石和冰渍向下爬移并黏合堆积,会存在掉落垮塌的情形,威胁到下坡地区的基础设施,研究发现多处掉落垮塌隐患区域,年平均形变速率最严重区域可达11.69 ± 0.07 [cm/yr],累计形变量为55.13 cm;夏季的雨季会加速斜坡滑动,冬季的旱季则阻止了大多数的滑动过程,可采取措施对地质灾害的发生进行预防,避免严重性的损失,为探索冻土形变所致地质灾害的成因奠定了基础。