多年冻土活动层厚度估算方法

1.本技术属于多年冻土监测技术领域，尤其涉及一种多年冻土活动层厚度估算方法。


背景技术：

2.活动层(active layer)是位于多年冻土上限和地表之间的土壤层，极易受到温度变化的影响而发生冻胀融沉。活动层厚度(active layerthickness，alt)指的是多年冻土每年融化的最大深度，也是监测多年冻土区域的不可或缺的气候变量之一。目前，多年冻土活动层厚度监测方法主要可分为传统的实地测量手段、基于遥感技术的方法和结合合成孔径雷达干涉测量(insar)技术的方法三类。
3.1)传统的实地测量手段通常依赖机械探针、冻融管、土温测量等方法来获取高质量的活动层厚度结果。这些方法费时费力成本高，仅能获取稀疏点尺度的活动层厚度结果，空间覆盖范围有限，无法体现多年冻土的空间异质性。
4.2)基于遥感技术的方法主要可分为探地雷达(gpr)技术和物理模型。相较于传统实地测量方法的稀疏点密度，gpr技术能够对活动层厚度进行几乎空间连续的有效估算，但该方法只适用于小范围的活动层厚度观测。物理模型的方法大多利用半经验或者统计模型建立活动层厚度与环境因子之间的关系，再结合环境因子遥感数据将活动层厚度外推插值至大范围区域，其中使用最为广泛的方法包括斯特藩(stefan)和库德里亚夫采夫(kudryavtsev)模型。然而，这类方法大多是基于区域的方法，难以应用于其它区域，同时部分输入参数难获取，且估算的活动层厚度结果分辨率较低，仅为一至几公里。
5.3)结合insar技术的方法，其首先利用具有高分辨率、高精度、大范围、长时段优势的insar技术准确获取冻土区地表形变，然后建立冻土区地表季节性形变与活动层变化之间的关系，从而达到大范围、高精度、高分辨率的活动层厚度估算。这类方法主要可分为两类。第一类方法是基于土壤一维热传导方程来反演活动层厚度，其假设insar获取的冻土时序形变和温度之间的时间延迟是温度从地表扩散至活动层底部的时间，然后依据土壤一维热传导方程建立该时间延迟和活动层厚度之间的函数关系。但是，该方法将整个区域假设为均匀体，采用了相同的土壤热扩散系数，因此仅适用于小区域的活动层厚度估算。另一类方法假设融化季冻土地表沉降完全是由活动层融化造成的，然后基于水质量守恒建立季节性形变和活动层厚度之间的积分关系，再通过逆运算估算活动层厚度。但该方法假设在活动层融化过程中，所有土壤水分都参与了水冰相变，而未考虑土壤负温状态下未冻水的影响。此外，仅利用insar一维形变观测值反演活动层厚度，会导致获取的形变结果出现误差，进而降低估算结果的精度和可靠性。
6.根据以上对现有多年冻土活动层厚度估算方法的分析，可以看出虽然insar技术已经成为大范围、高分辨率的活动层厚度估算的重要工具之一；但其对insar形变结果的依赖性高，且未全面考虑实际冻土变化的机理，会影响活动层厚度估算的精度和可靠性，从而影响该技术的应用和推广。


技术实现要素：

7.本技术实施例提供了一种多年冻土活动层厚度估算方法，可以解决活动层厚度估算的精度和可靠性低的问题。
8.本技术实施例提供了一种多年冻土活动层厚度估算方法，包括：
9.获取估算区域的多个不同轨道的sar影像数据集；
10.分别对多个sar影像数据集进行时序insar处理，得到每个sar影像数据集的los向地表时序形变；
11.根据每个sar影像数据集的los向地表时序形变，获取估算区域的地表真实二维形变结果；
12.获取估算区域的遥感数据和土壤数据，并根据遥感数据和土壤数据对估算区域的未冻水含量进行预报；
13.根据地表真实二维形变结果和未冻水含量，对估算区域的活动层厚度进行估算。
14.可选的，根据每个sar影像数据集的los向地表时序形变，获取估算区域的地表真实二维形变结果，包括：
15.分别将每个sar影像数据集的los向地表时序形变分解为垂直、东西和南北三个方向的形变分量，得到每个轨道的观测方程；
16.利用多个轨道的观测方程获取估算区域的地表真实二维形变结果。
17.可选的，第i个轨道的观测方程为：
[0018][0019]
其中，表示第i个轨道的sar影像数据集的los向地表时序形变，du表示在垂直方向的形变分量，de表示在东西方向的形变分量，dn表示在南北方向的形变分量，θi表示第i个轨道的sar影像数据集的入射角，αi表示第i个轨道的sar影像数据集的方位角，i＝1,
…
,j，j表示轨道的总数。
[0020]
可选的，利用多个轨道的观测方程获取估算区域的地表真实二维形变结果，包括：
[0021]
利用最小二乘法对进行求解，得到估算区域的地表真实二维形变结果；地表真实二维形变结果包括du和de。
[0022]
可选的，遥感数据包括多层土壤的温度数据，土壤数据包括每层土壤的土质类型以及每种土质类型对应的未冻水含量常数；
[0023]
根据遥感数据和土壤数据对估算区域的未冻水含量进行预报，包括：
[0024]
利用公式wu＝aθ-b
计算估算区域的未冻水含量；
[0025]
其中，wu表示估算区域的未冻水含量，θ表示温度数据的绝对值，a和b均表示未冻
水含量常数。
[0026]
可选的，在根据每个sar影像数据集的los向地表时序形变，获取估算区域的地表真实二维形变结果之后，方法还包括：
[0027]
利用傅里叶拟合方法对地表真实二维形变结果中的du进行拟合；
[0028]
从拟合结果中提取融化季时段对应的垂直方向的形变分量
[0029]
可选的，根据地表真实二维形变结果和未冻水含量，对估算区域的活动层厚度进行估算，包括：
[0030]
通过公式计算估算区域从活动层融化开始t0至结束时刻tk的活动层厚度h；
[0031]
其中，δhr表示估算区域在t
r-1
至tr时刻融化的活动层厚度，时刻融化的活动层厚度，时刻融化的活动层厚度，表示在tr时刻的值，表示在t
r-1
时刻的值，w
ice
表示估算区域在融化季时段的t
r-1
至tr时刻参与水冰相变过程的冰含量，ρ
water
表示水的密度，ρ
ice
表示冰的密度，wr表示估算区域在tr时刻的土壤水分含量，表示估算区域在t
r-1
时刻的未冻水含量。
[0032]
本技术的上述方案有如下的有益效果：
[0033]
在本技术的一些实施例中，通过对多个不同轨道的sar影像数据集分别进行时序insar处理得到多个los向地表时序形变，然后利用多个los向地表时序形变获取估算区域的地表真实二维形变结果，从而有效避免了单轨insar形变观测值难以反映地表真实形变特征的问题，同时在利用地表真实二维形变结果进行活动层厚度估算时，全面考虑了实际冻土变化过程中未冻水的影响，从而使得本技术的估算方法能大大提高活动层厚度估算的精度和可靠性。
[0034]
本技术的其它有益效果将在随后的具体实施方式部分予以详细说明。
附图说明
[0035]
为了更清楚地说明本技术实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本技术的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0036]
图1为本技术一实施例提供的多年冻土活动层厚度估算方法的流程图；
[0037]
图2为本技术一实例中估算区域高相干点融化季开始时间统计直方图；
[0038]
图3为本技术一实例中估算区域高相干点融化季结束时间统计直方图；
[0039]
图4为本技术一实例中多年冻土活动层厚度实测值与估算结果对比图。
具体实施方式
[0040]
以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本技术实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本技术。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本技术的描述。
[0041]
应当理解，当在本技术说明书和所附权利要求书中使用时，术语“包括”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
[0042]
还应当理解，在本技术说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
[0043]
如在本技术说明书和所附权利要求书中所使用的那样，术语“如果”可以依据上下文被解释为“当...时”或“一旦”或“响应于确定”或“响应于检测到”。类似地，短语“如果确定”或“如果检测到[所描述条件或事件]”可以依据上下文被解释为意指“一旦确定”或“响应于确定”或“一旦检测到[所描述条件或事件]”或“响应于检测到[所描述条件或事件]”。
[0044]
另外，在本技术说明书和所附权利要求书的描述中，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
[0045]
在本技术说明书中描述的参考“一个实施例”或“一些实施例”等意味着在本技术的一个或多个实施例中包括结合该实施例描述的特定特征、结构或特点。由此，在本说明书中的不同之处出现的语句“在一个实施例中”、“在一些实施例中”、“在其他一些实施例中”、“在另外一些实施例中”等不是必然都参考相同的实施例，而是意味着“一个或多个但不是所有的实施例”，除非是以其他方式另外特别强调。术语“包括”、“包含”、“具有”及它们的变形都意味着“包括但不限于”，除非是以其他方式另外特别强调。
[0046]
针对目前活动层厚度估算的精度和可靠性低的问题，本技术实施例提供了一种多年冻土活动层厚度估算方法，该方法通过对多个不同轨道的sar影像数据集分别进行时序insar处理得到多个los向地表时序形变，然后利用多个los向地表时序形变获取估算区域的地表真实二维形变结果，从而有效避免了单轨insar形变观测值难以反映地表真实形变特征的问题，同时在利用地表真实二维形变结果进行活动层厚度估算时，全面考虑了实际冻土变化过程中未冻水的影响，从而使得本技术的估算方法能大大提高活动层厚度估算的精度和可靠性。
[0047]
下面结合具体实施例对本技术提供的多年冻土活动层厚度估算方法进行示例性的说明。
[0048]
本技术实施例提供了一种多年冻土活动层厚度估算方法，该方法可以由终端设备执行，也可以由应用于终端设备中的装置(比如芯片)来执行，下述实施例以该方法由终端设备执行为例。作为一种示例，该终端设备可以是平板，服务器或者笔记本电脑等，本技术实施例对此不做限定。
[0049]
如图1所示，本技术实施例提供的多年冻土活动层厚度估算方法包括如下步骤：
[0050]
步骤11，获取估算区域的多个不同轨道的sar影像数据集。
[0051]
上述估算区域为需进行多年冻土活动层厚度估算的区域；上述多个不同轨道的合成孔径雷达(sar，syntheticaperture radar)影像数据集可以是合成孔径雷达卫星在多个
不同轨道对估算区域进行观测得到的。
[0052]
可以理解的是，为便于终端设备对估算区域进行多年冻土活动层厚度估算，合成孔径雷达卫星在得到上述sar影像数据集后，可将得到的sar影像数据集传输给终端设备。
[0053]
步骤12，分别对多个sar影像数据集进行时序insar处理，得到每个sar影像数据集的los向地表时序形变。
[0054]
在本技术的一些实施例中，计算第i个轨道的sar影像数据集的los向地表时序形变的具体实现方式可以为：
①
对第i个轨道的sar影像数据集的ni(ni第i个轨道的sar影像数据集的总景数)景sar影像进行配准；
②
选取合适的时空基线阈值，对满足阈值条件的干涉对采用外部数字高程模型(dem)数据进行差分干涉处理得到mi景差分干涉图；
③
对差分干涉图进行滤波、解缠，获得解缠后的差分干涉图；
④
通过合适的相干性阈值选出高相干点，并利用之前的mi景解缠差分干涉图和短基线集合成孔径雷达干涉测量(sbas-insar，small baseline subset insar)技术解算各高相干点的视线(los，line ofsight)向地表时序形变需要说明的是，由于时序insar处理为sar影像数据集的常用处理方式，因此，在此不对时序insar处理的过程进行过多赘述。
[0055]
可以理解的是，为方便后续的处理，可将不同轨道解算得到的高相干点los向地表时序形变结果地理编码、重采样至相同坐标系和分辨率。
[0056]
步骤13，根据每个sar影像数据集的los向地表时序形变，获取估算区域的地表真实二维形变结果。
[0057]
在本技术的一些实施例中，为提高地表真实形变结果的准确性，可通过雷达成像的几何关系，结合多个不同轨道的los向地表时序形变，将一维los向形变分解为二维或三维地表真实形变，进而利用最小二乘法计算获取估算区域的地表真实二维形变结果。
[0058]
步骤14，获取估算区域的遥感数据和土壤数据，并根据遥感数据和土壤数据对估算区域的未冻水含量进行预报。
[0059]
在本技术的一些实施例中，可通过遥感数据采集设备(如卫星)采集得到上述遥感数据，该遥感数据包括估算区域的多层土壤的温度数据。
[0060]
在本技术的一些实施例中，土壤数据包括上述多层土壤中每层土壤的土质类型以及每种土质类型对应的未冻水含量常数。
[0061]
需要说明的是，上述土质类型可通过常用的土质识别方法对每层土壤的样本进行识别得到，上述未冻水含量常数为与未冻水含量相关的参数，具体可以通过实地钻孔测量的土壤温度和土壤水分数据结合“两点法”计算获取，更具的，可以通过构建土质类型与未冻水含量(即土壤水分数据)之间的关系式(该关系式包含未冻水含量常数)，然后在估算区域的多个位置钻孔，测量每个位置的土壤温度和土壤水分数据，再利用采集的数据代入关系式对未冻水含量常数进行求解，得到未冻水含量常数。
[0062]
在本技术的一些实施例中，具体可利用公式wu＝aθ-b
计算估算区域的未冻水含量。其中，wu表示估算区域的未冻水含量，θ表示温度数据的绝对值，a和b均表示未冻水含量常数。这里需要说明的是，由于反演过程都是逐像素反演的，也就是说在计算未冻水含量时，需针对每个像素点，利用该像素点对应的每一层土壤温度、未冻水含量常数进行计算，最终得到估算区域的未冻水含量。
[0063]
步骤15，根据地表真实二维形变结果和未冻水含量，对估算区域的活动层厚度进
行估算。
[0064]
值得一提的是，在本技术的一些实施例中，通过对多个不同轨道的sar影像数据集分别进行时序insar处理得到多个los向地表时序形变，然后利用多个los向地表时序形变获取估算区域的地表真实二维形变结果，从而有效避免了单轨insar形变观测值难以反映地表真实形变特征的问题，同时在利用地表真实二维形变结果进行活动层厚度估算时，全面考虑了实际冻土变化过程中未冻水的影响，从而使得本技术的估算方法能大大提高活动层厚度估算的精度和可靠性。
[0065]
下面结合具体实施例对地表真实二维形变结果的获取进行示例性的说明。
[0066]
在本技术的一些实施例中，上述步骤13，根据每个sar影像数据集的los向地表时序形变，获取估算区域的地表真实二维形变结果的具体实现方式包括如下步骤：
[0067]
步骤一，分别将每个sar影像数据集的los向地表时序形变分解为垂直、东西和南北三个方向的形变分量，得到每个轨道的观测方程。
[0068]
在本技术的一些实施例中，第i个轨道的观测方程为：
[0069][0070]
其中，表示第i个轨道的sar影像数据集的los向地表时序形变，du表示在垂直方向的形变分量，de表示在东西方向的形变分量，dn表示在南北方向的形变分量，θi表示第i个轨道的sar影像数据集的入射角，αi表示第i个轨道的sar影像数据集的方位角，i＝1,
…
,j，j表示轨道的总数。
[0071]
步骤二，利用多个轨道的观测方程获取估算区域的地表真实二维形变结果。
[0072]
需要说明的是，要利用insar los向地表时序形变获取地表的真实三维形变结果，通常需要利用三个及以上独立的成像几何获取的多个测量值，从一维形变拓展至三维。考虑到差分干涉测量(d-insar)对南北向形变不敏感，为了保证最小二乘估算的正确性，通常忽略南北向形变，结合多个轨道的los向地表时序形变，解算估算区域地表真实二维形变结果。
[0073]
具体的，可在忽略南北方向形变分量的基础上，利用最小二乘法对具体的，可在忽略南北方向形变分量的基础上，利用最小二乘法对进行求解，得到估算区域的地表真实二维形变结果，该地表真实二维形变结果包括du和de。
[0074]
在本技术的一些实施例中，在得到估算区域的地表真实二维形变结果后，可利用
傅里叶拟合方法对地表真实二维形变结果中的du进行拟合，并从拟合结果中提取融化季时段对应的垂直方向的形变分量以便对活动层厚度进行估算。
[0075]
需要说明的是，傅里叶拟合方法对du进行拟合，得到的拟合结果可以为一曲线，该曲线的横轴为时间，纵轴为du的值，因此能从拟合结果中简单、快速提取出融化季时段对应的垂直方向的形变分量
[0076]
在本技术的一些实施例中，在得到地表真实二维形变结果和未冻水含量后，可通过公式计算估算区域从活动层融化开始t0至结束时刻tk的活动层厚度h。
[0077]
其中，δhr表示估算区域在t
r-1
至tr时刻融化的活动层厚度，由于地表沉降完全是由活动层中冰水相变所造成的体积差引起的，可以建立δdr与δhr之间的关系为：之间的关系为：表示在tr时刻的值，表示在t
r-1
时刻的值，w
ice
表示估算区域在融化季时段的t
r-1
至tr时刻参与水冰相变过程的冰含量。
[0078]
假设融化季估算区域内多年冻土中的水处于饱和状态，且地表沉降完全是由于活动层中冰与水的相变所造成的体积差引起的。考虑到t
r-1
时刻冻结的活动层中仍然存在一定量的未冻水结合水质量守恒，可以得到参与水冰相变的冰含量为：ρ
water
表示水的密度，ρ
ice
表示冰的密度，wr表示估算区域在tr时刻的土壤水分含量，可通过遥感数据采集设备(如卫星)采集得到，表示估算区域在t
r-1
时刻的未冻水含量，即wu在t
r-1
时刻的值。
[0079]
由上可见，上述估算方法通过充分利用多轨道insar形变观测值，避免了单轨insar形变观测值难以反映地表真实形变特征的现象，并且在对冻土形变与活动层变化建模时全面考虑了实际冻土变化过程中未冻水的影响，同时方法易于实施，应用广泛，有助于提高大范围、高分辨率的insar活动层厚度估计的精度和可靠性。
[0080]
为便于理解上述估算方法，在此以一具体实例对上述估算方法进行示例性说明。
[0081]
在该实例中，采集的数据包括：覆盖了青藏高原北麓河多年冻土区(覆盖范围：92.7801
°
e-93.9707
°
e，34.5721
°
n-35.4318
°
n)的哨兵一号(sentinel-1)卫星升轨和降轨影像数据和土壤水分主动-被动(smap)遥感土壤温湿度数据。
[0082]
在该实例中，估算区域包含了楚玛尔河高平原、可可西里山和北麓河盆地；整个区域大多为富冰多年冻土，地下冰含量丰富，最大体积含冰量达到了70％，平均过剩冰含量可达19％。多年冻土年平均温度介于-1.5～0℃之间，活动层厚度介于1.4～3.4m之间。insar数据采用sentinel-1升轨和降轨影像，其中升轨影像轨道和帧(frame)编号分别为41和90，降轨影像的轨道和frame编号分别为150和495，时间范围为2018年5月10日至2020年7月22日。smap遥感土壤温湿度数据提供土壤表层5cm处和土壤1m深处的土壤温度和湿度数据，其时间覆盖范围与sar影像一致。
[0083]
在该实例中，对估算区域的活动层厚度进行估算的具体步骤如下：
[0084]
步骤1：对sentinel-1升降和降轨影像分别进行配准，配准精度优于1/1000像素；选取时空基线阈值分别为150m和36天，对满足该阈值的干涉对采用航天飞机雷达地形测绘使命(srtm)dem数据进行差分干涉处理得到一系列的差分干涉图；采用goldstein滤波算法(goldstein滤波算法是一种干涉图滤波算法)对差分干涉图进行二维频率域滤波，并利用最小网络费用流方法进行二维相位解缠，为了避免低相干点对解缠的影响，在解缠过程中去除了相干性低于0.2的像素点，最后获取了解缠后的差分干涉图；选取最小相干性和平均相干性分别不低于0.35和0.85的像素点为高相干点，通过sbas-insar技术逐像素求解出这些点的在升轨和降轨los向上的地表时序形变和将上述得到的地表时序形变结果编码、重采样到相同的地理坐标系和分辨率下。
[0085]
步骤2：在忽略南北向形变的基础上，根据雷达成像的几何关系，将已获取的los向地表时序形变和视为观测值，可建立如下观测方程：
[0086][0087]
其中，θ
asc
表示对应的影像数据的入射角，α
asc
表示对应的影像数据的方位角，θ
dsc
表示对应的影像数据的入射角，α
dsc
表示对应的影像数据的方位角。根据最小二乘平差原理即可对上述观测方程进行求解得到地表真实二维时序形变du和de。考虑到观测方程可能存在不适定问题，此时可通过吉洪诺夫(tikhonov)正则化的方法对其进行求解。
[0088][0089]
其中，b为公式中的系数矩阵；p为权阵，在本实例中取单位权阵；[t
0 t1ꢀ…ꢀ
tn]为sar影像的获取时刻。
[0090]
获取各高相干点在垂直向上的时序形变du之后，考虑到冻土形变一般由气候变暖引起的长期线性形变和冻融循环有关的周期性形变两部分组成，首先对各高相干点在垂直向上的时序形变du进行去线性趋势，然后利用公式傅里叶拟合法对每个高相干点的垂直向时序形变进行参数拟合，确定每年融化季的时间(如图2和图3所示)并提取处于融化季之中的时序形变其中，图2中的横轴为融化季开始时间，图3中的横轴为融化季结束时间，图2和图3中纵轴为确定融化季开始和结束时间的频率。
[0091]
步骤3：根据估算区域周边的地质资料，可以获取估算区域的土壤类型和与未冻水含量相关的经验参数(即前文的未冻水含量常数a和b)，具体参数见
[0092]
表1。
[0093]
深度/m土壤类型未冻水含量常数a未冻水含量常数b0-2砂土0.070.172-10黏土0.120.15》10砂岩0.010.10
[0094]
表1
[0095]
根据步骤2获取的融化季时间，可以从smap土壤温度数据中获取各土壤层在融化前的温度，将其与表1数据代入公式wu＝aθ-b
中，即可实现未冻水含量的预报，计算得到各土壤层在融化前的未冻水含量表示wu在ti时刻的值。
[0096]
步骤4：将步骤2得到的融化季垂直向地表时序形变步骤3得到的各土壤层融化前未冻水含量和smap分层土壤水分数据代入公式和smap分层土壤水分数据代入公式中，得到活动层估算结果。
[0097]
由图4可以看出，通过与估算区域三个多年冻土监测站点实测的活动层厚度进行对比，本技术提供的估算方法均最接近与实测值，这说明本技术提供的估算方法估算出的活动层厚度更为准确，也验证了本技术提供的估算方法能够有效的提高大范围、高分辨率的insar活动层厚度估计的精度和可靠性。其中，图4中wd3、hr3和qt08分别表示3个不同的多年冻土监测站点，实测值指的是对应站点的活动层厚度的实际值，升轨和降轨分别表示利用单轨(升轨/降轨)sar数据进行活动层厚度预测的方法。
[0098]
以上所述是本技术的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本技术所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本技术的保护范围。