CN116051620A - 基于InSAR技术冻土区活动层厚度估计方法和系统 - Google Patents

Abstract

本发明公开了一种基于InSAR技术冻土区活动层厚度估计方法和系统，该方法包括以下步骤：首先获取冻土实验区时间序列SAR图像、土壤含水量数据、土壤孔隙度数据并对数据预处理；然后对时间序列SAR图像进行差分干涉处理；再构建InSAR季节性形变模型；然后基于InSAR季节性形变模型的进行时序InSAR解算流程；最后构建InSAR季节性形变的活动层厚度反演模型，并根据冻土实验区的季节性形变获取大范围多年冻土区的活动层厚度。本发明可应用于大范围多年冻土区活动层厚度的估计，能够大范围、高精度、高分辨率的获取多年冻土区的活动层厚度结果，这对冻土区的环境、水文生态、工程建筑的运营具有重要科学与工程意义。

Description

基于InSAR技术冻土区活动层厚度估计方法和系统

技术领域

本发明涉及冰川冻土学领域，尤其涉及一种基于InSAR技术冻土区活动层厚度估计方法和系统。

背景技术

活动层厚度（Active Layer Thickness，ALT）在调节冻土区水文动态过程和气候变化起着至关重要的作用。冻土区近地表温度升高会导致多年冻土退化和活动层厚度增加，进而会引起局部塌陷和地表沉降。活动层厚度水热动态变化与土壤含水量变化，空气温度变化和地表、次地表热辐射传输过程密切相关，且这种变化对冻土区的水文和植物生态系统具有直接影响，因此实现多年冻土活动层厚度的大范围、高精度的估计对冻土区的环境、水文生态、工程建筑运营具有重要科学与工程意义。

目前冻土活动层厚度监测的方法主要分为实测法、经验/半经验模型法、基于遥感技术的活动层厚度估计。实测法包括金属探钎法、钻探法、测温法、探地雷达（GPR）法等，它通过实地测量冻土融化层和冻结层界限，是最直接的冻土活动层厚度的调查方法。经验/半经验模型法利用温度、积雪、土质等因素和冻融深度等因素建立经验关系，它适合空间异质性较强的冻土，在大尺度冻土区的活动层厚度监测得到了广泛应用。经验/半经验模型法分为Stefan模型、Kudryavtsev模型、Nelson模型、地球物理研究所多年冻土实验室2（GIPL2）模型等。基于遥感技术的活动层厚度估计包括基于光学和基于微波遥感两种手段来反演冻土活动层厚度。其中基于微波遥感的后向散射技术的活动层厚度反演主要利用SAR后向散射系数与实测活动层厚度建立的经验关系模型。近年来随着合成孔径雷达干涉（Interferometry Synthetic Aperture Radar，InSAR）技术的迅猛发展，多时相InSAR技术（Multi-Temporal InSAR，MT-InSAR）已广泛应用于监测冻土区的季节性地表形变，而MT-InSAR技术反演的季节性形变结果已经成功应用于估算冻土区地表活动层厚度，基于InSAR、MT-InSAR技术的冻土区活动层厚度反演成为目前青藏高原冻土领域研究的热点。

但是上述方法还存在以下问题：（1）实测法虽然精度高，但是仅能获取点尺度上的冻土活动层厚度，针对大区域尺度的冻土活动层厚度监测，由于环境及监测成本等因素，实测法难以满足应用需求。（2）经验/半经验模型法存在着模型复杂、计算量大、输入参数较多等问题，且生成的活动层厚度结果图分辨率较低，难以满足精细尺度下的应用要求。（3）基于InSAR技术的活动层厚度反演模型存在着地下土壤含水量时空分布无法准确表征，InSAR地表形变与活动层厚度关联不准确等问题。例如冻土中土壤含水量的影响往往被忽略或视为一个常数，土壤含水量的空间变化只在点尺度上进行判别，且没有考虑多层土壤和分层土壤孔隙度等条件。

为了获取大范围、高精度、高分辨率冻土区的活动层厚度信息，本发明提供了一种基于InSAR技术大范围多年冻土区活动层厚度估计方法和系统，基于InSAR技术获取的活动层厚度分辨率远优于基于经验、数值模型的结果，且获取的活动层厚度范围远远大于点尺度的实测方法，该方法在估计大范围多年冻土区活动层厚度方面具有广阔的应用前景。

发明内容

本发明的目的在于针对现有技术估计活动层厚度方法存在费时费力、测量密度低、模型参数较多等劣势，提供一种基于InSAR技术冻土区活动层厚度估计方法和系统。本发明能够实现大范围多年冻土区活动层厚度参数的估计。

本发明的目的是通过以下技术方案来实现的：本发明实施例第一方面提供了一种基于InSAR技术冻土区活动层厚度估计方法，包括以下步骤：

（1）获取冻土实验区的时间序列SAR图像数据和气象水文数据，并对时间序列SAR图像数据进行预处理，以获取配准后的时间序列SAR图像；

（2）对所述步骤（1）获得的配准后的时间序列SAR图像进行差分干涉处理，以获取差分干涉图；

（3）构建符合冻土物理形变过程的InSAR季节性形变模型；

（4）基于InSAR季节性形变模型对时间序列SAR图像和所述步骤（2）获得的差分干涉图进行时序InSAR解算，以获取季节性形变的结果；

（5）构建活动层厚度反演模型，并根据所述步骤（4）求解的冻土实验区的季节性形变的结果获取大范围多年冻土区的活动层厚度参数。

可选地，所述气象水文数据包括数字高程模型数据、日空气温度数据、土壤含水量数据、土壤孔隙度数据和土壤类型数据。

可选地，所述步骤（1）包括以下子步骤：

（1.1）获取数据：根据冻土实验区的经纬度地理位置获取时间序列SAR图像数据和数字高程模型数据，同时获取冻土实验区地表下不同深度的土壤含水量数据和土壤孔隙度数据，还需获取冻土实验区的日空气温度数据、体积含冰量数据、年平均地温数据、土壤类型数据和地形坡度数据；

（1.2）预处理：首先对时间序列SAR图像进行数据导入，以将原始SAR图像格式导入生成单视复数据格式，同时结合下载的精轨数据文件进行轨道参数的更新；然后选择一个主图像，对时间序列SAR图像进行图像配准，以将所有辅图像重采样到主图像的框架下，获取配准后的时间序列SAR图像。

可选地，所述步骤（2）包括以下子步骤：

（2.1）差分干涉：对所述步骤（1）获得的配准后的时间序列SAR图像进行差分干涉处理，对主图像和辅图像进行相位差分，以获取干涉图；

（2.2）地形和平地相位去除：利用数字高程模型数据根据反距离权重插值算法计算出每个像素的地形相位，同时计算出地形相位和平地相位并去除地形相位和平地相位，以获取差分干涉图；

（2.3）滤波：对所述步骤（2.2）获取的差分干涉图进行滤波处理，以获取滤波后的差分干涉图；

（2.4）相位解缠：对所述步骤（2.3）获取的滤波后的差分干涉图进行相位解缠处理，以获取相位解缠后的差分干涉图。

可选地，所述步骤（3）中构建符合冻土物理形变过程的InSAR季节性形变模型具体为：根据冻土区冻融过程与形变的关系，从冻土冻融物理过程出发，根据冻土实验区的日空气温度数据、体积含冰量数据、年平均地温数据、土壤类型数据、地形坡度数据、维度和地表高程构建符合冻土物理形变过程的InSAR季节性形变模型；

所述InSAR季节性形变模型的表达式为：

其中，D是InSAR观测期间内的累积形变；为冻土活动层引起的地表季节性的位移，为融合的冻胀和融沉过程的联合指数，主要与季节性融化沉降和冻结抬升的解冻和冻结累计天数的平方根以及冻结和融化n因子有关，S表示冻土在InSAR观测期间内的季节性形变速率；为位于多年冻土上限的蠕变形变，为体积含冰量；为多年冻土由于升温导致其物理力学性质改变引起的沉降，为年平均地表温度；和分别为体积含冰量和年平均地表温度因素的待求参数，e是残余形变项，t表示SAR图像获取的时刻。

可选地， 所述季节性形变包括活动层的冻融循环形变、多年冻土上限处的融沉形变以及多年冻土层的蠕变形变。

可选地，所述步骤（4）包括以下子步骤：

（4.1）将所述步骤（3）中构建的InSAR季节性形变模型的冻土形变项加入到InSAR干涉模型相位方程中，通过时序形变解算求取InSAR干涉模型相位方程中的形变信息；

（4.2）根据所述步骤（4.1）的InSAR干涉模型相位方程中的形变信息和所述步骤（2）获得的差分干涉图采用新小基线集方法对冻土形变参数的季节性形变量进行求解；

（4.3）根据季节性形变量利用最小二乘算法实现冻土实验区每个像素季节性形变参数的求解，以获取季节性形变的结果。

可选地，所述新小基线集方法具体为：对于M幅SAR图像和所述步骤（2）获得的N个相位解缠后的差分干涉图中的第 l个像素，其相位解缠后的干涉相位观测集合等效于每个像素单独SAR图像相位值的线性组合。

可选地，所述步骤（5）包括以下子步骤：

（5.1）根据所述步骤（1）获取的土壤含水量数据和土壤孔隙度数据以及所述步骤（4）求解的冻土区季节性形变结果数据，构建活动层厚度反演模型；

所述活动层厚度反演模型的表达式如下：

其中，H是冻土融沉期的最大活动层厚度，为时间间隔，i为SAR图像序列号，即，为第i景SAR图像与第j景SAR图像的时间间隔内的融化深度，为第i景SAR图像与第j景SAR图像的时间间隔内的季节性沉降量，为纯水的密度，为纯冰的密度，q为不同地下深度的第q层土壤含水量图层，为相邻两个时相的第k层的土壤含水量均值，为第q层的土壤孔隙度。

（5.2）根据所述步骤（4）求解的冻土实验区的季节性形变的结果和所述步骤（5.1）构建的活动层厚度反演模型获取大范围多年冻土区的活动层厚度参数。

本发明实施例第二方面提供了一种基于InSAR技术冻土区活动层厚度估计系统，用于实现上述的基于InSAR技术冻土区活动层厚度估计方法，包括：

数据获取和预处理模块：用于获取冻土实验区的时间序列SAR图像数据、数字高程模型数据、日空气温度数据、土壤含水量数据、土壤孔隙度数据、土壤类型数据、体积含冰量数据、年平均地温数据和地形坡度数据，并对时间序列SAR图像数据进行预处理；

时间序列SAR图像差分干涉流程模块：用于根据配准后的时间序列SAR图像进行差分干涉流程；

时序InSAR解算模块：用于构建冻土实验区的InSAR季节性形变模型，并采用新小基线集方法对冻土形变参数进行求解，获取冻土实验区的季节性形变量；和

活动层厚度估计模块：用于构建基于InSAR季节性形变的活动层厚度估计模型，计算多年冻土区大范围活动层厚度参数。

本发明的有益效果是，本发明基于InSAR技术的活动层厚度进行反演，根据冻土融化期的实际形变量并基于反演模型估计活动层厚度，符合活动层厚度形变规律，为冻土区活动层厚度反演提供了新的技术手段，该技术估计的活动层厚度具有空间分辨率高、覆盖范围广等优势，这种方法在冰川冻土学领域估计冻土区大范围活动层厚度具有极大的应用潜力，也为研究冻土区的环境、水文、生态和寒区工程设施的灾害防治和监测具有重要科学意义；本发明能够实现大范围多年冻土区活动层厚度参数的估计。

附图说明

图1是本发明的基于InSAR技术冻土区活动层厚度估计方法流程图；

图2是本发明实施例GPR实测活动层厚度结果和Sentinel-1与TerraSAR-X传感器反演活动层厚度结果对比图；其中，图2中的（a）为高寒草甸活动层厚度结果；图2中的（b）为高寒荒漠活动层厚度结果；

图3是本发明的基于InSAR技术冻土区活动层厚度估计系统的结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的基于InSAR技术冻土区活动层厚度估计方法，如图1所示，具体包括以下步骤：

（1）获取冻土实验区的时间序列SAR（合成孔径雷达，Synthetic Aperture Radar）图像数据和气象水文数据，并对时间序列SAR图像数据进行预处理，以获取配准后的时间序列SAR图像。

本实施例中，气象水文数据包括DEM数据、日空气温度数据、土壤含水量数据、土壤孔隙度数据和土壤类型数据。

（1.1）获取数据：根据冻土实验区的经纬度地理位置获取时间序列SAR图像数据和DEM（数字高程模型，Digital Elevation Model）数据，同时获取冻土实验区地表下不同深度的土壤含水量数据和土壤孔隙度数据，为构建符合冻土物理形变过程的InSAR季节性形变模型，还需获取冻土实验区的日空气温度数据、土壤类型数据。

具体地，根据冻土实验区的经纬度地理位置获取DEM数据和时间序列SAR图像，其中搜集的SAR图像数据集见表1。本实施例中，此次冻土实验区为青海省西南部北麓河地区，该区域地处青藏高原永久冻土区，气候为寒冷半干旱气候，平均气温约-3.8℃，冻土的体积冰含量大于25%，年平均地温大于-1.0℃，年平均降水量为300 ~ 400mm。本次实验获取重返周期12天、250公里幅宽、VV极化方式、N.150降轨模式2017年8月7日至2019年10月25日期间的全部Sentinel-1数据以更精确反演北麓河地区的冻土活动层厚度参数。同时采集TerraSAR-X数据进行对比实验，TerraSAR-X影像为聚束模式，降轨轨道号为N.385，HH极化方式，覆盖面积约为2.8×7.5 km²，时间跨度为2018年12月15日至2019年10月8日。

表1：冻土实验区SAR数据集详细信息

同时获取冻土实验区不同深度的分层土壤含水量和分层土壤孔隙度的数据，另外，为构建符合冻土物理形变过程的InSAR季节性形变模型，还需获取冻土实验区的日空气温度数据、体积含冰量数据、年平均地温数据、土壤类型数据、地形坡度数据。土壤类型数据来源于中国科学院南京土壤研究所提供的《1:100 万中华人民共和国土壤图》，并根据土体导热系数对土质进行赋值。使用欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts，ECMWF）综合预测系统模型的第五代大气再分析全球气候数据（ERA5-Interim）的日空气温度数据、分层土壤含水量数据，其中分层土壤含水量数据为四个土壤层（0~7 cm、7~21 cm、21~72 cm、72~189 cm）且分辨率为的数据产品，该数据可在欧洲中期天气预报中心https://apps.ecmwf.int.网站获得。分层土壤孔隙度的数据见表2。从国家青藏高原科学数据中心获取北麓河地区的体积含冰量数据、年平均地温数据、土壤类型数据，从北麓河地区的DEM数据利用ArcGIS10.6软件提取冻土实验区的地形坡度、纬度和高程数据。

表2：北麓河地区土壤孔隙度数据

（1.2）预处理：首先对时间序列SAR图像进行数据导入，以将原始SAR图像格式导入生成单视复数据格式，同时结合下载的精轨数据文件进行轨道参数的更新；然后选择一个主图像，对时间序列SAR图像进行图像配准，以将所有辅图像重采样到主图像的框架下，获取配准后的时间序列SAR图像。

具体地，采用开源软件GMTSAR进行处理，该软件可通过https://topex.ucsd.edu/gmtsar/ 网站获取。首先对时间序列SAR图像进行数据导入，即将原始SAR图像tiff格式导入生成单视复数据格式，同时结合下载的精轨数据文件进行轨道参数的更新用于后续的图像配准。对于TerraSAR-X数据采用基于精轨数据的几何配准进行图像配准，对于Sentinel-1数据采用基于精轨数据的几何配准和增强谱分集方法进行图像配准，选择主图像对时间序列SAR图像进行配准，即将所有辅图像重采样到主图像的框架下，从而得到配准后的时间序列SAR图像。

（2）对步骤（1）获得的配准后的时间序列SAR图像进行差分干涉处理，以获取差分干涉图。应当理解的是，步骤（2）的过程中均可采用开源软件GMTSAR进行处理，此乃InSAR的差分干涉处理流程。

（2.1）差分干涉：对步骤（1）获得的配准后的时间序列SAR图像进行差分干涉处理，对主图像和辅图像进行相位差分，以获取干涉图。应当理解的是，本实施例中采用开源软件GMTSAR对配准后的时间序列SAR图像进行差分干涉处理。

（2.2）地形和平地相位去除：利用DEM数据根据反距离权重插值算法计算出每个像素的地形相位，同时计算出地形相位和平地相位并去除地形相位和平地相位，以获取差分干涉图。

（2.3）滤波：为抑制差分干涉图中的噪声影响，通过滤波操作对步骤（2.2）获取的差分干涉图进行滤波处理，以获取滤波后的差分干涉图。

具体地，可以选择非线性自适应Goldstein空域滤波对生成的差分干涉图进行滤波处理，有利于抑制差分干涉图中的噪声影响。

（2.4）相位解缠：对步骤（2.3）获取的滤波后的差分干涉图进行相位解缠处理，以获取相位解缠后的差分干涉图。

具体地，可以对滤波后的差分干涉图采用最小费用流算法进行相位解缠处理。根据相位值可以得出两次成像中微波的路程差，从而计算出目标地区的地形、地貌以及表面的微小变化，可用于数字高程模型建立、地壳形变探测等，即可从干涉条纹中获取地形高程数据等。

（3）构建符合冻土物理形变过程的InSAR季节性形变模型。

本实施例中，根据冻土区冻融过程与形变的关系，从冻土冻融物理过程出发，建立拟合冻土物理形变过程的形变模型以提高InSAR技术形变反演精度，如考虑日空气温度数据、体积含冰量、年平均地表温度、土壤类型数据、地形坡度数据、维度和地表高程等因素。将多年冻土区地表总形变分为活动层的冻融循环形变（与多年冻土活动层冻融过程有关）、多年冻土上限处的融沉形变（与多年冻土上限富含冰层有关）以及多年冻土层的蠕变形变（由多年冻土升温导致土力学性质改变而生），从而构建符合冻土物理形变过程的InSAR季节性形变模型，其表达式如下：

其中，D是InSAR观测期间内的累积形变；为冻土活动层引起的地表季节性的位移，为融合的冻胀和融沉过程的联合指数，主要与季节性融化沉降和冻结抬升的解冻和冻结累计天数的平方根以及冻结和融化n因子有关，S表示冻土在InSAR观测期间内的季节性形变速率；为位于多年冻土上限的蠕变形变，为体积含冰量；为多年冻土由于升温导致其物理力学性质改变引起的沉降，为年平均地表温度；和分别为体积含冰量和年平均地表温度因素的待求参数，e是残余形变项，t表示SAR图像获取的时刻。

进一步地，的计算公式为：

其中，和分别表示融化沉降和冻胀抬升的累积日数，该参数由日空气温度数据获得；为尺度因子，和分别为冻土区冻结期和融化期的土壤导热率，和分别为冻结n因子和融化n因子（n-factor），用于表征冻土与大气圈之间能量的交换过程，该过程是冻土形成的主要原因，该参数可由地表温度数据和日空气温度数据求解。

需要说明的是，和两个参数可由步骤（1）获取ERA5-Interim再分析日空气温度数据求解。示例性地，根据冻土实验区的日空气温度数据可得到冻土的融化和冻结日期（此次实验中取日空气温度数据大于0摄氏度确定为融化日期，小于或等于0摄氏度确定为冻结日期），同时结合SAR图像观测时期，计算出与冻土的融化和冻结日期相差的天数，即可得到和两个参数。

进一步地，的计算公式为：

其中，为体积含冰量（%），ST表示土壤类型，NDVI为归一化植被指数，该参数可由欧空局的Sentine-2影像的近红外（NIR）与红光（Red）的波段信息计算，即NDVI = (NIR -Red) / (NIR + Red)，SD表示地形坡度信息。

需要说明的是，土壤类型、地形坡度信息为步骤（1）所获取的数据。

进一步地，的计算公式如下：

其中，表示年平均地温（℃），N和H分别表示纬度(°)和地表高程（m）。

需要说明的是，纬度、地表高程为步骤（1）所获取的数据。

（4）基于InSAR季节性形变模型对时间序列SAR图像和步骤（2）最终获得的差分干涉图进行时序InSAR解算，以获取季节性形变的结果。

（4.1）将步骤（3）中构建的InSAR季节性形变模型的冻土形变项加入到InSAR干涉模型相位方程中，通过时序形变解算即可求得InSAR干涉模型相位方程中的形变相位项。其中，构建的InSAR干涉模型相位方程见下式：

其中，表示InSAR干涉模型相位方程，n为时序InSAR干涉对的第n个干涉图，x为步骤（2.1）求得的干涉图中的某一像素点；表示主图像m相位图的第x个像素的干涉相位值，表示辅图像s相位图的第x个像素的干涉相位值，为T时刻的主图像m，为时刻T的辅图像s；代表InSAR季节性形变模型；为第n个干涉图第x个像素的形变相位值，def表示形变相位项；为第n个干涉图第x像素的大气相位值，atmos表示大气相位项；为第n个干涉图第x像素的地形相位值，topo表示地形相位项；为第n个干涉图第x像素的轨道误差相位值，orb表示轨道误差相位项；为第n个干涉图第x像素的残余噪声相位值，res表示残余相位项。

（4.2）根据步骤（4.1）的InSAR干涉模型相位方程中的形变信息和步骤（2）获得的差分干涉图采用NSBAS（新小基线集，New Small BAseline Subset）方法对冻土形变参数的季节性形变量进行求解。

本实施例中，NSBAS（新小基线集）方法具体为：对于M幅SAR图像和步骤（2.4）生成的N个相位解缠后的差分干涉图中的第 l个像素，步骤（2.4）中相位解缠后的干涉相位观测集合可以写成每个像素单独SAR图像相位值的线性组合，具体见下式：

其中，为个相位解缠后的差分干涉图像素 l的相位向量，为由0、-1、1组成的维系数矩阵，为像素 l组成的相位增量，即待求解的解缠相位值，为第i景和第j景SAR图像像素 l组成的解缠相位值，为SAR图像采集时间n与n + 1之间的像素 l的相位增量。具体地，该公式的求解方法可用最小二乘算法进行解算，当矩阵为奇异矩阵时，可采用奇异值分解算法（Singular Value Decomposition，SVD）进行求解。然而，当矩阵存在秩亏问题时，即在差分干涉图网络连接中缺少一个关键链接时，在这种情况下，SAR图像构成的网络连接被分割成两个或多个独立的图像组，即有个别图像或者图像组没有被连接上，SVD算法代替最小二乘算法来反演时间序列位移时，它将连续图像组之间的相位延迟增量设为零。因此会使时序形变解算产生偏差。为了克服这一缺陷，NSBAS在反演中添加约束条件来进行优化形变解算过程。像素 l在时刻的累积相位延迟定义为：

其中，k为像素 l的第k个相位增量，为InSAR相位中像素 l第k个相位增量的模型相位，为k个相位增量的序号，为像素 l两个连续SAR图像相位延迟增量。

然后将约束添加到形变反演中，表示为：

其中，为像素 l在t时刻的线性形变速率，为像素 l第一个相位增量的时刻，为像素 l第k个相位增量的时刻，为像素 l在t时刻的线性形变加速度，为像素 l第k个相位增量跟垂直基线相关的DEM误差，表示像素 l跟垂直基线相关的参数，表示干涉图像素 l第k个相位增量时刻的垂直基线，为干涉图像素 l在t时刻的残余形变。因此针对没有连接的小基线干涉对的情况，通过添加权重来进行计算，即见下式：

基于上式，将步骤（3）获得的InSAR季节性形变模型应用到NSBAS方法中，对冻土实验区的冻融形变进行监测，可见下式：

其中，S表示季节性形变速率。

（4.3）根据季节性形变量利用最小二乘算法实现冻土实验区每个像素季节性形变参数的求解，以获取季节性形变的结果，例如，基于时序Sentinel-1图像求解的北麓河地区季节性形变的结果包括高寒荒漠和高寒草甸相对应的结果，该地区基于Sentinel-1数据季节性形变范围为- 54.50 mm~4.50mm。

（5）构建活动层厚度反演模型，并根据步骤（4）求解的冻土实验区的季节性形变的结果获取大范围多年冻土区的活动层厚度参数。

（5.1）根据步骤（1）获取的土壤含水量数据和土壤孔隙度数据以及步骤（4）求解的冻土区季节性形变结果数据，构建活动层厚度反演模型，其表达式为：

其中，H为冻土融沉期的最大活动层厚度，为时间间隔，i为第i景SAR图像，即，为第i景SAR图像与第j景SAR图像的时间间隔内的融化深度，为第i景SAR图像与第j景SAR图像的时间间隔内的季节性沉降量，由步骤（4）求解得到；为纯水的密度，为纯冰的密度，q为不同地下深度的第q层土壤含水量图层，为相邻两个时相的第k层的土壤含水量均值，该参数可由ERA5-Interim再分析分层土壤含水量计算；为第q层的土壤孔隙度，该参数可由步骤（1）获取的分层土壤含水量数据求解。

应当理解的是，由于步骤（1）中获取了冻土实验区地表下不同深度的土壤含水量数据和土壤孔隙度数据，根据分层土壤含水量数据建立对应层的土壤孔隙度数据，即可很容易得知相邻两个时相的第q层的土壤含水量均值和第q层的土壤孔隙度。

本实施例中，该活动层厚度反演模型考虑了垂直方向的土壤含水量和土壤孔隙度的因素，因此估计的活动层厚度更精确。

（5.2）根据步骤（4）求解的冻土实验区的季节性形变的结果和步骤（5.1）构建的活动层厚度反演模型获取大范围多年冻土区的活动层厚度参数。

具体地，将步骤（4）求解的冻土实验区的季节性形变的结果代入步骤（5.1）构建的活动层厚度反演模型中，即可求得大范围多年冻土区的活动层厚度参数。应当理解的是，根据步骤（4）求解的冻土实验区的季节性形变的结果可以进一步求解得到第i景SAR图像与第j景SAR图像的时间间隔内的季节性沉降量。

示例性地，基于时序Sentinel-1图像北麓河地区反演的活动层厚度结果与基于TerraSAR-X图像北麓河地区反演的活动层厚度结果分别为0.3 ~ 4.23m和0.3 ~ 4.04m。为了更详细分析北麓河地区不同地貌景观活动层厚度的时空变化和验证本发明反演的活动层厚度的精度，使用2018年8月28日和2020年9月24日~2020年10月7日在北麓河进行两次野外探地雷达（GPR）实验数据进行对比分析。GPR数据采集的具体地理位置包括了两种地貌，分别是采集到的高寒荒漠的现场照片和高寒草甸的现场照片。GPR实测活动层厚度剖面结果与这两种传感器反演的活动层厚度结果进行对比分析见图2，其中，图2中的（a）为高寒草甸活动层厚度结果，图2中的（b）为高寒荒漠活动层厚度结果。结果分析由于TerraSAR-X和Sentinel-1空间分辨率较高，反演的活动层厚度结果分辨率较高。两种传感器反演的活动层厚度变化基本符合GPR数据的探测剖面变化，证明了本发明提供的多年冻土区活动层厚度估计方法可靠性。

值得一提的是，本发明实施例还提供了一种基于InSAR技术冻土区活动层厚度估计系统，用于实现上述的基于InSAR技术冻土区活动层厚度估计方法。

本实施例中，该系统包括数据获取和预处理模块、时间序列SAR图像差分干涉流程模块、时序InSAR解算模块和活动层厚度估计模块，如图3所示。

本实施例中，数据获取和预处理模块用于获取冻土实验区的时间序列SAR图像数据、DEM数据、日空气温度数据、土壤含水量数据、土壤孔隙度数据、土壤类型数据、体积含冰量数据、年平均地温数据和地形坡度数据，并对时间序列SAR图像数据进行预处理。需要说明的是，用户可输入冻土监测区的经纬度信息，系统会自动下载对应区域的SAR数据。

本实施例中，时间序列SAR图像差分干涉流程模块用于根据配准后的时间序列SAR图像进行差分干涉流程。其中，该模块包括差分干涉、地形和平地相位去除、滤波及相位解缠流程。需要说明的是，用户可自行配置不同流程的算法参数。

本实施例中，时序InSAR解算模块用于构建冻土实验区的InSAR季节性形变模型，并采用NSBAS方法对冻土形变参数进行求解，获取冻土实验区的季节性形变量。需要说明的是，用户可根据需求选择不同的形变模型和时序InSAR解算算法。

本实施例中，活动层厚度估计模块用于构建基于InSAR季节性形变的活动层厚度估计模型，计算多年冻土区大范围活动层厚度参数。需要说明的是，用户可根据需求选择输出的活动层厚度结果的范围和输出格式。

以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (10)

1.一种基于InSAR技术冻土区活动层厚度估计方法，其特征在于，包括以下步骤：

（1）获取冻土实验区的时间序列SAR图像数据和气象水文数据，并对时间序列SAR图像数据进行预处理，以获取配准后的时间序列SAR图像；

（2）对所述步骤（1）获得的配准后的时间序列SAR图像进行差分干涉处理，以获取差分干涉图；

（3）构建符合冻土物理形变过程的InSAR季节性形变模型；

（4）基于InSAR季节性形变模型对时间序列SAR图像和所述步骤（2）获得的差分干涉图进行时序InSAR解算，以获取季节性形变的结果；

（5）构建活动层厚度反演模型，并根据所述步骤（4）求解的冻土实验区的季节性形变的结果获取大范围多年冻土区的活动层厚度参数。

2.根据权利要求1所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于，所述气象水文数据包括数字高程模型数据、日空气温度数据、土壤含水量数据、土壤孔隙度数据和土壤类型数据。

3.根据权利要求2所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于，所述步骤（1）包括以下子步骤：

（1.1）获取数据：根据冻土实验区的经纬度地理位置获取时间序列SAR图像数据和数字高程模型数据，同时获取冻土实验区地表下不同深度的土壤含水量数据和土壤孔隙度数据，还需获取冻土实验区的日空气温度数据、体积含冰量数据、年平均地温数据、土壤类型数据和地形坡度数据；

（1.2）预处理：首先对时间序列SAR图像进行数据导入，以将原始SAR图像格式导入生成单视复数据格式，同时结合下载的精轨数据文件进行轨道参数的更新；然后选择一个主图像，对时间序列SAR图像进行图像配准，以将所有辅图像重采样到主图像的框架下，获取配准后的时间序列SAR图像。

4.根据权利要求1所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于，所述步骤（2）包括以下子步骤：

（2.1）差分干涉：对所述步骤（1）获得的配准后的时间序列SAR图像进行差分干涉处理，对主图像和辅图像进行相位差分，以获取干涉图；

（2.2）地形和平地相位去除：利用数字高程模型数据根据反距离权重插值算法计算出每个像素的地形相位，同时计算出地形相位和平地相位并去除地形相位和平地相位，以获取差分干涉图；

（2.3）滤波：对所述步骤（2.2）获取的差分干涉图进行滤波处理，以获取滤波后的差分干涉图；

（2.4）相位解缠：对所述步骤（2.3）获取的滤波后的差分干涉图进行相位解缠处理，以获取相位解缠后的差分干涉图。

5.根据权利要求1所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于，所述步骤（3）中构建符合冻土物理形变过程的InSAR季节性形变模型具体为：根据冻土区冻融过程与形变的关系，从冻土冻融物理过程出发，根据冻土实验区的日空气温度数据、体积含冰量数据、年平均地温数据、土壤类型数据、地形坡度数据、维度和地表高程构建符合冻土物理形变过程的InSAR季节性形变模型；

所述InSAR季节性形变模型的表达式为：

与S的乘积加上与的乘积，再加上与的乘积，再加上e得到D；

其中，D是InSAR观测期间内的累积形变；为冻土活动层引起的地表季节性的位移，为融合的冻胀和融沉过程的联合指数，主要与季节性融化沉降和冻结抬升的解冻和冻结累计天数的平方根以及冻结和融化n因子有关，S表示冻土在InSAR观测期间内的季节性形变速率；为位于多年冻土上限的蠕变形变，为体积含冰量；为多年冻土由于升温导致其物理力学性质改变引起的沉降，为年平均地表温度；和分别为体积含冰量和年平均地表温度因素的待求参数，e是残余形变项，t表示SAR图像获取的时刻。

6.根据权利要求5所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于， 所述季节性形变包括活动层的冻融循环形变、多年冻土上限处的融沉形变以及多年冻土层的蠕变形变。

7.根据权利要求1所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于，所述步骤（4）包括以下子步骤：

（4.1）将所述步骤（3）中构建的InSAR季节性形变模型的冻土形变项加入到InSAR干涉模型相位方程中，通过时序形变解算求取InSAR干涉模型相位方程中的形变信息；

（4.2）根据所述步骤（4.1）的InSAR干涉模型相位方程中的形变信息和所述步骤（2）获得的差分干涉图采用新小基线集方法对冻土形变参数的季节性形变量进行求解；

（4.3）根据季节性形变量利用最小二乘算法实现冻土实验区每个像素季节性形变参数的求解，以获取季节性形变的结果。

8.根据权利要求7所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于，所述新小基线集方法具体为：对于M幅SAR图像和所述步骤（2）获得的N个相位解缠后的差分干涉图中的第l个像素，其相位解缠后的干涉相位观测集合等效于每个像素单独SAR图像相位值的线性组合。

9.根据权利要求1所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于，所述步骤（5）包括以下子步骤：

（5.1）根据所述步骤（1）获取的土壤含水量数据和土壤孔隙度数据以及所述步骤（4）求解的冻土区季节性形变结果数据，构建活动层厚度反演模型；

所述活动层厚度反演模型的表达式如下：

个累加得到H，与的乘积除以与的差，再除以与乘积的平均值，再累加次；

其中，H是冻土融沉期的最大活动层厚度，为时间间隔，i为SAR图像序列号，即，为第i景SAR图像与第j景SAR图像的时间间隔内的融化深度，为第i景SAR图像与第j景SAR图像的时间间隔内的季节性沉降量，为纯水的密度，为纯冰的密度，q为不同地下深度的第q层土壤含水量图层，为相邻两个时相的第k层的土壤含水量均值，为第q层的土壤孔隙度；

（5.2）根据所述步骤（4）求解的冻土实验区的季节性形变的结果和所述步骤（5.1）构建的活动层厚度反演模型获取大范围多年冻土区的活动层厚度参数。

10.一种基于InSAR技术冻土区活动层厚度估计系统，用于实现权利要求1-9中任一项所述的基于InSAR技术冻土区活动层厚度估计方法，其特征在于，包括：

数据获取和预处理模块：用于获取冻土实验区的时间序列SAR图像数据、数字高程模型数据、日空气温度数据、土壤含水量数据、土壤孔隙度数据、土壤类型数据、体积含冰量数据、年平均地温数据和地形坡度数据，并对时间序列SAR图像数据进行预处理；

时间序列SAR图像差分干涉流程模块：用于根据配准后的时间序列SAR图像进行差分干涉流程；

时序InSAR解算模块：用于构建冻土实验区的InSAR季节性形变模型，并采用新小基线集方法对冻土形变参数进行求解，获取冻土实验区的季节性形变量；和

活动层厚度估计模块：用于构建基于InSAR季节性形变的活动层厚度估计模型，计算多年冻土区大范围活动层厚度参数。

Images