联合主/被动微波遥感技术的活动层季节性冻融变化监测研究

寇 晓 康，刘 信 达，李 旭 东

(1.石家庄铁道大学土木工程学院，河北 石家庄 050043；2.道路与铁道工程安全保障省部共建教育部重点实验室(石家庄铁道大学)，河北 石家庄 050043；3.中核三维地理信息工程技术研究中心，河北 石家庄 050043)

0 引言

多年冻土[1]一般分布在高纬度或高海拔地区，我国青藏高原地区多年冻土覆盖区域占比近70%[2]。天然地面到多年冻土上限(多年冻土的顶板)之间是活动层(又称季节融化层)，活动层季节性变化引起的地表冻胀融沉是冻土区铁路、公路、管线、建筑物等产生形变的主要原因之一。近年来，随着全球不断升温[3]，多年冻土上限下移，活动层厚度及冻胀融沉幅度均随之发生改变，深刻影响着寒区基础设施的运营与安全。因此，开展活动层的冻融变化监测研究至关重要。

传统活动层冻融变化监测方法通过埋设传感器获取温湿度数据，辅助利用精密水准仪、GNSS等传统测量手段获取冻融引起的形变信息，得到的点位信息精确，数据准确度高，但点位布设密度低会产生较大误差，有些地方因很难抵达而无法监测，同时存在观测周期长、成本高等问题。相比而言，微波遥感穿透力强，且不受云、雨和太阳照度影响[4]，可用于大范围活动层的冻融变化形变监测，相关研究常采用被动微波遥感(Passive Microwave Remote Sensing,PMRS)技术，其基本探测原理为：水和冰的介电特性在微波波段差异巨大，当活动层内的土壤水发生冰/水相变之后，根据微波传感器接收的微波辐射亮温差异，可反推出土壤的冻结/融化状态。常用的基于PMRS技术的地表冻融监测方法有单/双指标算法[5-8]、决策树算法[9]、季节阈值算法[8,10]、判别式算法[11-13]等，由此生成了多种地表冻融状态数据集[14-16]。上述算法中用到的亮温数据多来源于SMMR、SSM/I、SSMIS、AMSR-E、AMSR-2、FY3B-MWRI、SMOS、SMAP等卫星传感器，时间分辨率较高(1～2 d)，但空间分辨率低，单个像元对应的地面单元边长约为25～40 km。因被动微波穿透能力与活动层内土壤水分含量及探测波段有关，而被动微波一般只能探测到土壤最表层的冻融状态信息[17]，对青藏高原地区而言，利用PMRS探测到的主要是活动层表层的冻融状态。在活动层季节性变化引起的冻胀融沉形变监测方面，一般采用主动微波遥感(Active Microwave Remote Sensing,AMRS)中的雷达差分干涉测量(D-InSAR)技术，利用合成孔径雷达(SAR)两次观测所获得的雷达复数影像数据中的相位信息，通过过滤观测区域的地形及平地相位反演地表形变，但该技术易受时(空)间失相干和大气延迟影响，精度难达要求，为此，学者们提出了以永久散射体法(PS)和小基线集(SBAS)方法为代表的时序差分雷达干涉测量 (MT-InSAR) 技术[18-20]，其克服了传统D-InSAR技术的缺点，理论精度能达到毫米级[18]。目前已将MT-InSAR技术应用于青藏高原的冻土形变监测中[21-29]，研究所需数据多来源于ERS、RADARSAT、ALOS、TerraSAR-X、ENVISAT、Sentinel-1、GF-3等卫星，空间分辨率在1 m～1 km之间，时间分辨率在11～46 d之间[30]。

活动层表层的冻融变化与整个活动层的冻融周期紧密相关，伴随活动层冻胀融沉的变化表层会产生形变,可见二者间有天然联系，但受制于PMRS的低空间分辨率及AMRS的低时间分辨率，现有研究鲜有从遥感视角分析二者间的关联。鉴于此，本文以青藏高原五道梁地区为研究区，拟联合主/被动微波遥感技术获取地表冻融与形变信息，并结合地面实测数据深入剖析活动层季节性冻融变化中冻融状态、土壤温度、地表形变、土壤湿度间的内在关联，从而明确主/被动微波遥感监测信息与活动层内部温湿度参量的协同变化机制，以期从遥感视角深入揭示活动层的季节性变化过程。

1 研究区及数据

1.1 研究区概况

五道梁地区位于青藏高原腹地青海省西南部曲麻莱县境内，海拔4 600～4 800 m，属高原山地气候，四季皆冬,1月平均气温为-16.2 ℃，7月平均气温为6.0 ℃，全年平均气温为-5.1 ℃，是全国最低值。研究区南部地表覆盖类型主要为高寒草甸，北部主要为裸地，该地区还分布着山体、热融湖塘等，青藏铁路贯穿其中(图1)，研究该地区地表冻融变化及形变对于冻土区铁路的维护具有重要意义。

图1 青藏高原五道梁地区地理位置及地表覆盖类型Fig.1 Geographical location and land cover types of Wudaoliang area in the Qinghai-Tibet Plateau

1.2 研究数据

(1)AMRS数据，来源于Sentinel-1卫星系统，下载网址https://search.asf.alaska.edu/#/，该系统由Sentinel-1A和Sentinel-1B两颗SAR卫星构成，单颗星重访周期为12 d，双星协同工作的重访周期为6 d。2017年10月至2018年12月研究区Sentinel-1B卫星数据缺失较严重，考虑到数据的一致性，本研究选取Sentinel-1A 卫星在干涉宽幅模式(IW)下的C波段VV极化方式升轨影像，共计36景，影像距离向分辨率为5 m，方位向分辨率为20 m，空间范围如图1所示。

(2)PMRS数据，来源于美国冰雪数据中心发布的增强型空间分辨率的亮度温度(CETB)数据集，该数据集包含SSM/I、SSMIS、AMSR-E、AMSR-2等传感器的长时序亮温数据，且各数据源之间进行了一致性校准，下载网址https://nsidc.org/data/NSIDC-0630/versions/1。选取DMSP F18卫星SSMIS传感器2017年10月至2018年12月的升轨数据，时间间隔为1 d,数据采用网格化存储，投影方式为EASE-Grid 2.0。不同通道的空间分辨率有区别，其中19 GHz亮温的分辨率为6.25 km，37 GHz亮温的分辨率为3.125 km。本文首先对19 GHz亮温进行降尺度，使之与37 GHz相同，再进行冻融判别。

(3)地面观测数据，来源于青藏高原多年冻土综合监测数据集[31,32]，可从国家青藏高原科学数据中心下载，选取位于研究区编号为QT08 2017年10月至2018年12月的站点数据，该站点对表层以下10 cm、40 cm、120 cm、200 cm、240 cm深度处的土壤温度和土壤湿度进行了长期观测。

(4)其他辅助数据，包括30 m空间分辨率的地表分类数据集(http://www.globallandcover.com/)以及30 m空间分辨率的SRTM DEM数据集(https://lpdaac.usgs.gov/)。

2 研究方法

本文技术流程(图2)为：首先利用Sentinel-1 SLC雷达影像数据，基于SBAS-InSAR监测方法获取研究区内的地表时序形变信息，并将其作为线性形变分量(年际形变)与周期性形变分量(季节形变)的组合，通过逐像元线性趋势拟合及趋势剔除，获取研究区地表季节形变信息，在此基础上通过空间聚合提取站点所对应微波像元(约3 km)的平均时序形变信息；其次利用SSMIS亮温数据，基于冻融判别式算法，通过对亮温重采样获取研究区内约3 km×3 km像元尺度的每日地表冻融状态信息；最后结合地面站点观测数据，分析活动层季节性冻融变化过程中的内部物理参量变化机制及其与外在遥感特征之间的关联。

图2 技术流程Fig.2 Technical flow chart

2.1 SBAS-InSAR方法

SBAS-InSAR方法是对相干目标进行相位分析以获取时序形变，通过选择合适的空间基线和时间基线阈值组成差分干涉对，然后选取相干目标点利用线性相位变化模型进行建模和解算，通过时空滤波去除大气延迟，在减少D-InSAR处理过程中的失相干影响及高程、大气误差的同时获取地表形变时间序列[20]。本文SBAS-InSAR方法实现流程为：根据研究区范围对卫星影像进行剪裁，将生成的数据进行干涉像对配对，加入外部DEM数据去除地形相位；依据小基线集技术原理，设置最大临界基线百分比为2%、时间基线为60 d，选择20～30个相干性较好的点作为相位解缠的参考点；将干涉图进行自适应滤波，使用最小费用流(MCF)方法进行相位解缠和基线精化估计；对残余的相位进行自适应滤波和第二次解缠，对输入的干涉图进行优化，在第一步得到的形变速率基础上进行大气滤波，从而估算和去除大气相位，最后将结果编码到WGS-84坐标系下，得到时序形变结果。

2.2 去趋势方法

多年冻土区地表形变分为年际尺度形变与季节尺度形变：在年际尺度上，夏季多年冻土区地表差异主要由多年冻土层上限的上下移动导致；在季节尺度上，地表形变主要由活动层的周期性冻胀融沉导致，即夏季活动层中的冰融化成水，地表发生下沉，冬季活动层水分重新凝结成冰，地表发生隆升，呈现周期性的形变特征[33]。本研究假设年际形变量在1年内为线性形变，可通过方程(Y=K×X+B，Y为年际形变量，X为时间(年)，K和B为线性拟合系数)进行拟合，基于该方程可获取各像元的年际变化趋势；假定季节形变量为原始形变量在年际变化趋势线上的投影转换分量D，计算公式为：

(1)

式中：X0、Y0为原始形变量在以时间和形变量组成的平面坐标系中的坐标。

2.3 重采样方法

本文采用双三次插值方法将SSMIS 19 GHz亮温数据的空间分辨率从6.25 km降至3.125 km，该方法用三阶多项式函数通过相邻的网格点拟合曲面，最后根据曲面拟合函数得到查询点的插值结果。具体计算参考文献[34]。

2.4 地表冻融判别式算法

冻融判别式算法比其他国际上主流算法的判对率高[14]，对裸地、植被或积雪覆盖下的土壤冻融情况均能进行有效区分，因此，本研究采用该算法进行活动层表层冻融状态的识别。该算法包括36.5 GHz V极化下的亮温(TB36.5V)和18.7 GHz H极化与36.5 GHz V极化下的亮温比(Qe18.7H/36.5V)两个主要指标，分别用于指示地表温度和土壤发射率的变化，具体判别公式为：

(2)

式中：DF、DT分别为冻土和融土的判别方程函数值,若DF-DT>0，则判别为冻土，否则判别为融土。

3 结果分析

3.1 季节形变量提取

3.1.1 假定合理性分析 本研究假定基于SBAS-InSAR方法获取的形变结果由线性形变分量(年际形变量)与周期性形变分量(季节形变量)两部分组成。由于研究期较短，为验证年际形变趋势是否存在，首先对时序形变结果进行逐像元线性拟合，统计各像元时序形变线性拟合方程的显著性检验结果(表1)并进行可视化(图3)，进而判定年际形变趋势的稳定性。由表1可知，显著相关区域面积之和占统计总面积的77.41%，表明研究区内绝大部分区域存在较明显的年际变化趋势，证明了本文假设的合理性。进一步对QT08站点3.125 km范围内的像元形变平均值进行统计(图4)，可以看出，其地表形变既表现出周期性特征，又表现出一定的年际沉降线性趋势。由于常年积雪等因素影响，部分区域失相干严重，造成形变结果存在少量的空白区(图3中的白色区域)，这部分并未参与计算。

图3 研究区逐像元时序形变线性拟合方程的显著性检验结果空间分布Fig.3 Spatial distribution of significance test results of linear fitting equation for time series deformation of pixels in the study area

图4 QT08站点对应的 3.125 km×3.125 km像元范围内形变平均值的时序变化Fig.4 Time series variation of average deformation in one 3.125 km×3.125 km pixel corresponding to QT08 station

表1 显著性检验标准及结果Table 1 Statistics of significance test criteria and results

3.1.2 季节形变量分布 利用2.2节的去趋势方法获取研究区2017年10月至2018年12月之间36景影像所对应的季节形变量的空间分布(图5)，可以看出，研究区大部分区域表现出较明显的季节性冻胀融沉趋势，其中北部尤为明显，南部稍弱，可能由于研究区北部以裸地为主，而南部以草地为主(图1)，受植被干扰产生失相干现象，进而对反演结果产生一定影响。进一步对QT08站点周围3.125 km的像元网格的季节形变量平均值进行可视化(图6)，可以看出，经去趋势处理后，该区域内的地表形变具有明显的冻胀融沉周期性变化趋势。

图5 研究区内季节形变量分布Fig.5 Distribution of seasonal deformation in the study area

图6 QT08站点对应的3.125 km×3.125 km像元范围内季节形变量平均值时序变化Fig.6 Time series variation of average seasonal deformation in one 3.125 km×3.125 km pixel corresponding to QT08 station

3.2 AMSR和PMSR监测结果对比分析

本研究利用SSMIS亮温数据，在对19 GHz亮温进行降尺度处理的基础上，基于冻融判别式算法最终获取研究区范围内3.125 km分辨率的每日地表冻融状态信息，然后提取QT08站点所对应的被动微波像元冻融状态信息，并与该地区的地表形变信息进行对比(图7)。可以看出，冻融状态及形变均表现出明显的周期特征，为更清晰地对比二者间的联系与差异，进行分时段对比分析。1)冻结期：自2017年10月18日开始，PMRS监测结果显示该区域地表经历一次冻融交替后持续处于冻结状态，说明自此进入冻结期；AMRS监测结果显示，随着地表冻结状态的持续，地表不断抬升，且有慢—快—慢—平稳的变化特征，而该特征在PMRS冻融状态监测信息中并未体现。2)冻融交替期：自2018年4月初开始，PMRS监测结果显示该区域地表进入冻融交替期；AMRS监测结果对应出现起伏特征，但变化量很小。3)融化期：自2018年5月中旬开始，PMRS监测结果显示除了偶尔的冻融交替现象之外(考虑可能为误判)，该区域地表始终保持融化状态，说明进入稳定的融化期；AMRS监测结果显示，随着地表融化状态的持续，地表不断沉降，沉降趋势同样有慢—快—慢—平稳的特征，但该特征在PMRS冻融状态监测信息中并未体现。4)再次冻融交替期：自2018年9月初开始，PMRS监测结果显示该区域再次进入冻融交替期；AMRS监测结果显示形变量再次出现波动，且波动幅度较小。5)再次冻结期：自2018年10月中旬之后，该区域地表持续处于冻结状态，地表不断抬升，说明再次进入冻结期。

图7 地表冻融状态与季节形变量的时序变化对比Fig.7 Temporal variation comparison of surface freezing-thawing state and seasonal deformation

综上，在活动层季节性变化的各个周期内，PMRS地表冻融状态监测结果与AMRS季节形变量提取结果在时间上均具有较强的关联性与一致性，但值得注意的是，在冻结期和融化期内，AMRS所显示的地表形变速率特征并未在PMRS监测结果中体现。

3.3 PMRS监测结果与活动层土壤温度对比分析

有学者依据活动层冻融过程不同阶段的水热特征差异，将青藏高原活动层的年际变化过程划分为夏季融化、秋季冻结、冬季降温和春季升温4个阶段[35]。为明晰主/被动微波遥感技术在活动层季节性冻融变化监测中的作用，本研究将年际冻融变化过程重新划分为稳定冻结期、稳定融化期和两次冻融交替期。由于被动微波的穿透能力有限，PMRS的监测结果主要反映近地表约0～5 cm的土壤冻融状态，现结合QT08站点在10 cm、40 cm、120 cm、200 cm和240 cm深度处的实测温度进行对比分析(图8)。从图8可以看出，基于PMRS获取的冻融状态变化情况与各层温度之间既存在明显的相关性，又存在一定的差异性，具体到各时期总结如下：1)冻结期：PMRS监测结果显示，两次地表冻结期分别自2017年10月18日和2018年10月18日左右开始，而从QT08站点各层土壤温度看，浅层土壤(40 cm以内)冻结开始时间与其大体相当，伴随外界低温向下传导，冻结起始时间随活动层深度的增加有不断向后延迟趋势。2)两次冻融交替期：PMRS监测结果显示，2018年4月初及9月下旬的冻融交替现象主要出现在活动层浅层，120 cm及以下温度振幅较小，几乎无明显的冻融交替现象。由于活动层浅层厚度较小，推断其为冻融交替期季节形变幅度较小的原因。3)融化期：PMRS监测结果显示，自2018年5月初表层融化状态基本保持稳定，活动层浅层土壤的融化期基本与PMRS监测结果一致，随着外界热量不断向下传导，虽然各层土壤温度均有抬升，但中层(120 cm)及以下土壤温度仍维持在0℃以下，其冻结期并未结束，随后经过约1个月的热量向下传导，这些层次才相继进入融化期。

图8 地表冻融状态与各层温度时序变化对比Fig.8 Temporal variation comparison of surface freezing-thawing state and soil temperature for various layers

综上可知，PMRS监测结果显示的冻结期/融化期的开始和结束时间与活动层浅层冻结期/融化期的开始和结束时间基本一致，活动层中层、深层冻结期/融化期的开始和结束时间相对滞后，其原因在于浅层土壤与外界环境交互较剧烈，而热量在土壤中的传递需要时间；此外，PMRS监测结果中的冻融交替现象主要出现在活动层的浅层，中层和深层土壤温度振幅较小，无明显冻融交替现象，因此冻融交替期内AMRS获取的季节形变幅度一般也较小。

3.4 AMRS季节形变监测结果与活动层土壤湿度对比分析

AMRS监测到的形变实际为活动层各层冰/水相变的积累，本文结合QT08站点在不同层次的实测土壤湿度进行对比分析。从图9可以看出，基于AMRS获取的季节形变量与各层湿度之间有明显的相关性，具体到各时期总结如下：1)冻结期：AMRS监测结果显示，自2017年10月中旬至2018年2月中旬期间，地面整体出现抬升趋势，对比各层土壤湿度变化可知，地面抬升主要与各层水分持续下降有关，冰/水相变(相变水量)速率决定了抬升趋势的变化速率，自2018年2月下旬开始，当水分全部冻结后，地表保持稳定状态。2)两次冻融交替期：AMRS监测结果显示，在2018年4月初及9月中旬，地表形态出现微小波动，这主要与浅层土壤的冰/水相变有关，在3.3节中已做出解释。3)融化期：AMRS监测结果显示,在2018年4月中旬至9月初期间，地面整体出现沉降趋势，对比各层土壤湿度变化可知，地面沉降主要与各层水分的回升有关，水分的回升速率与地表形变速率呈正相关，个别层次土壤湿度突然增高可能是降雨所致。综上可知：AMRS监测得到的季节形变量与活动层各层冰/水相变总量有关，季节形变速率与冰/水相变速率呈正相关，当活动层内部冰/水不再发生变化时，季节形变趋于稳定。

图9 季节形变反演结果与各层土壤水分时序变化情况对比Fig.9 Comparison between inversion results of seasonal deformation and temporal variation of soil moisture in each layer

4 结论

本研究以青藏高原五道梁地区为研究区，分别基于主/被动微波遥感技术获取了3 km×3 km分辨率的每日活动层表层冻融状态信息及5 m×20 m分辨率、时间间隔约12 d的活动层季节形变信息；在此基础上，通过空间重采样及时序拟合技术，结合地面实测土壤温度、湿度数据，对比分析了冻结期、融化期和冻融交替期内活动层各冻融参量之间的内在联系，得出以下结论：1)PMRS冻融状态监测结果可用于提取活动层表层冻融周期，该周期与活动层浅层冻融周期基本一致，活动层中层及深层的冻融起始时间相对滞后；2)AMRS季节形变监测结果与活动层内部总的相变水量紧密相关，形变速率与冰/水相变速率呈正相关关系；3)冻融交替现象主要发生在活动层的浅层，深层不明显，此为AMRS监测结果波动幅度较小的原因；4)PMRS监测结果与AMRS监测结果在时间上有较强的关联性和一致性，在冻融周期上的统一是季节形变监测中主/被动微波遥感协同应用的基础。本文关于AMRS和PMRS的联合分析对于大尺度范围内的活动层季节性变化监测研究具有重要意义。由于研究区内可用的站点数据较少，通过综合分析本文只得出一些定性结论，未来可结合更多实测数据，将AMRS与PMRS技术深入结合，应用于活动层冻融参量的定量监测中。