CN115540788A - 多年冻土活动层厚度估算方法 - Google Patents

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

Info

Publication number
CN115540788A
CN115540788A CN202211392436.0A CN202211392436A CN115540788A CN 115540788 A CN115540788 A CN 115540788A CN 202211392436 A CN202211392436 A CN 202211392436A CN 115540788 A CN115540788 A CN 115540788A
Authority
CN
China
Prior art keywords
deformation
estimation
earth surface
image data
soil
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202211392436.0A
Other languages
English (en)
Other versions
CN115540788B (zh
Inventor
王琪洁
李佳晨
张亚
杨沙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202211392436.0A priority Critical patent/CN115540788B/zh
Publication of CN115540788A publication Critical patent/CN115540788A/zh
Application granted granted Critical
Publication of CN115540788B publication Critical patent/CN115540788B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B15/00—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons
    • G01B15/02—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons for measuring thickness
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9027—Pattern recognition for feature extraction

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本申请适用于多年冻土监测技术领域，提供了一种多年冻土活动层厚度估算方法，其中该方法包括：获取估算区域的多个不同轨道的SAR影像数据集；分别对多个SAR影像数据集进行时序InSAR处理，得到每个SAR影像数据集的LOS向地表时序形变；根据每个SAR影像数据集的LOS向地表时序形变，获取估算区域的地表真实二维形变结果；获取估算区域的遥感数据和土壤数据，并根据遥感数据和土壤数据对估算区域的未冻水含量进行预报；根据地表真实二维形变结果和未冻水含量，对估算区域的活动层厚度进行估算。本申请能提高活动层厚度估算的精度和可靠性。

Description

多年冻土活动层厚度估算方法
技术领域
本申请属于多年冻土监测技术领域，尤其涉及一种多年冻土活动层厚度估算方法。
背景技术
活动层(Active Layer)是位于多年冻土上限和地表之间的土壤层，极易受到温度变化的影响而发生冻胀融沉。活动层厚度(Active LayerThickness，ALT)指的是多年冻土每年融化的最大深度，也是监测多年冻土区域的不可或缺的气候变量之一。目前，多年冻土活动层厚度监测方法主要可分为传统的实地测量手段、基于遥感技术的方法和结合合成孔径雷达干涉测量(InSAR)技术的方法三类。
1)传统的实地测量手段通常依赖机械探针、冻融管、土温测量等方法来获取高质量的活动层厚度结果。这些方法费时费力成本高，仅能获取稀疏点尺度的活动层厚度结果，空间覆盖范围有限，无法体现多年冻土的空间异质性。
2)基于遥感技术的方法主要可分为探地雷达(GPR)技术和物理模型。相较于传统实地测量方法的稀疏点密度，GPR技术能够对活动层厚度进行几乎空间连续的有效估算，但该方法只适用于小范围的活动层厚度观测。物理模型的方法大多利用半经验或者统计模型建立活动层厚度与环境因子之间的关系，再结合环境因子遥感数据将活动层厚度外推插值至大范围区域，其中使用最为广泛的方法包括斯特藩(Stefan)和库德里亚夫采夫(Kudryavtsev)模型。然而，这类方法大多是基于区域的方法，难以应用于其它区域，同时部分输入参数难获取，且估算的活动层厚度结果分辨率较低，仅为一至几公里。
3)结合InSAR技术的方法，其首先利用具有高分辨率、高精度、大范围、长时段优势的InSAR技术准确获取冻土区地表形变，然后建立冻土区地表季节性形变与活动层变化之间的关系，从而达到大范围、高精度、高分辨率的活动层厚度估算。这类方法主要可分为两类。第一类方法是基于土壤一维热传导方程来反演活动层厚度，其假设InSAR获取的冻土时序形变和温度之间的时间延迟是温度从地表扩散至活动层底部的时间，然后依据土壤一维热传导方程建立该时间延迟和活动层厚度之间的函数关系。但是，该方法将整个区域假设为均匀体，采用了相同的土壤热扩散系数，因此仅适用于小区域的活动层厚度估算。另一类方法假设融化季冻土地表沉降完全是由活动层融化造成的，然后基于水质量守恒建立季节性形变和活动层厚度之间的积分关系，再通过逆运算估算活动层厚度。但该方法假设在活动层融化过程中，所有土壤水分都参与了水冰相变，而未考虑土壤负温状态下未冻水的影响。此外，仅利用InSAR一维形变观测值反演活动层厚度，会导致获取的形变结果出现误差，进而降低估算结果的精度和可靠性。
根据以上对现有多年冻土活动层厚度估算方法的分析，可以看出虽然InSAR技术已经成为大范围、高分辨率的活动层厚度估算的重要工具之一；但其对InSAR形变结果的依赖性高，且未全面考虑实际冻土变化的机理，会影响活动层厚度估算的精度和可靠性，从而影响该技术的应用和推广。
发明内容
本申请实施例提供了一种多年冻土活动层厚度估算方法，可以解决活动层厚度估算的精度和可靠性低的问题。
本申请实施例提供了一种多年冻土活动层厚度估算方法，包括：
获取估算区域的多个不同轨道的SAR影像数据集；
分别对多个SAR影像数据集进行时序InSAR处理，得到每个SAR影像数据集的LOS向地表时序形变；
根据每个SAR影像数据集的LOS向地表时序形变，获取估算区域的地表真实二维形变结果；
获取估算区域的遥感数据和土壤数据，并根据遥感数据和土壤数据对估算区域的未冻水含量进行预报；
根据地表真实二维形变结果和未冻水含量，对估算区域的活动层厚度进行估算。
可选的，根据每个SAR影像数据集的LOS向地表时序形变，获取估算区域的地表真实二维形变结果，包括：
分别将每个SAR影像数据集的LOS向地表时序形变分解为垂直、东西和南北三个方向的形变分量，得到每个轨道的观测方程；
利用多个轨道的观测方程获取估算区域的地表真实二维形变结果。
可选的，第i个轨道的观测方程为：
Figure BDA0003932525470000031
其中，
Figure BDA0003932525470000032
表示第i个轨道的SAR影像数据集的LOS向地表时序形变，dU表示
Figure BDA0003932525470000033
在垂直方向的形变分量，dE表示
Figure BDA0003932525470000034
在东西方向的形变分量，dN表示
Figure BDA0003932525470000035
在南北方向的形变分量，θi表示第i个轨道的SAR影像数据集的入射角，αi表示第i个轨道的SAR影像数据集的方位角，i＝1,…,j，j表示轨道的总数。
可选的，利用多个轨道的观测方程获取估算区域的地表真实二维形变结果，包括：
利用最小二乘法对
Figure BDA0003932525470000036
进行求解，得到估算区域的地表真实二维形变结果；地表真实二维形变结果包括dU和dE。
可选的，遥感数据包括多层土壤的温度数据，土壤数据包括每层土壤的土质类型以及每种土质类型对应的未冻水含量常数；
根据遥感数据和土壤数据对估算区域的未冻水含量进行预报，包括：
利用公式Wu＝aθ-b计算估算区域的未冻水含量；
其中，Wu表示估算区域的未冻水含量，θ表示温度数据的绝对值，a和b均表示未冻水含量常数。
可选的，在根据每个SAR影像数据集的LOS向地表时序形变，获取估算区域的地表真实二维形变结果之后，方法还包括：
利用傅里叶拟合方法对地表真实二维形变结果中的dU进行拟合；
从拟合结果中提取融化季时段对应的垂直方向的形变分量
Figure BDA0003932525470000041
可选的，根据地表真实二维形变结果和未冻水含量，对估算区域的活动层厚度进行估算，包括：
通过公式
Figure BDA0003932525470000042
计算估算区域从活动层融化开始t0至结束时刻tk的活动层厚度H；
其中，ΔHr表示估算区域在tr-1至tr时刻融化的活动层厚度，
Figure BDA0003932525470000049
Figure BDA00039325254700000410
Figure BDA0003932525470000044
表示
Figure BDA0003932525470000045
在tr时刻的值，
Figure BDA0003932525470000046
表示
Figure BDA0003932525470000047
在tr-1时刻的值，Wice表示估算区域在融化季时段的tr-1至tr时刻参与水冰相变过程的冰含量，
Figure BDA0003932525470000048
ρwater表示水的密度，ρice表示冰的密度，Wr表示估算区域在tr时刻的土壤水分含量，
Figure BDA00039325254700000411
表示估算区域在tr-1时刻的未冻水含量。
本申请的上述方案有如下的有益效果：
在本申请的一些实施例中，通过对多个不同轨道的SAR影像数据集分别进行时序InSAR处理得到多个LOS向地表时序形变，然后利用多个LOS向地表时序形变获取估算区域的地表真实二维形变结果，从而有效避免了单轨InSAR形变观测值难以反映地表真实形变特征的问题，同时在利用地表真实二维形变结果进行活动层厚度估算时，全面考虑了实际冻土变化过程中未冻水的影响，从而使得本申请的估算方法能大大提高活动层厚度估算的精度和可靠性。
本申请的其它有益效果将在随后的具体实施方式部分予以详细说明。
附图说明
为了更清楚地说明本申请实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本申请一实施例提供的多年冻土活动层厚度估算方法的流程图；
图2为本申请一实例中估算区域高相干点融化季开始时间统计直方图；
图3为本申请一实例中估算区域高相干点融化季结束时间统计直方图；
图4为本申请一实例中多年冻土活动层厚度实测值与估算结果对比图。
具体实施方式
以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本申请实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本申请。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本申请的描述。
应当理解，当在本申请说明书和所附权利要求书中使用时，术语“包括”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
还应当理解，在本申请说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
如在本申请说明书和所附权利要求书中所使用的那样，术语“如果”可以依据上下文被解释为“当...时”或“一旦”或“响应于确定”或“响应于检测到”。类似地，短语“如果确定”或“如果检测到[所描述条件或事件]”可以依据上下文被解释为意指“一旦确定”或“响应于确定”或“一旦检测到[所描述条件或事件]”或“响应于检测到[所描述条件或事件]”。
另外，在本申请说明书和所附权利要求书的描述中，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
在本申请说明书中描述的参考“一个实施例”或“一些实施例”等意味着在本申请的一个或多个实施例中包括结合该实施例描述的特定特征、结构或特点。由此，在本说明书中的不同之处出现的语句“在一个实施例中”、“在一些实施例中”、“在其他一些实施例中”、“在另外一些实施例中”等不是必然都参考相同的实施例，而是意味着“一个或多个但不是所有的实施例”，除非是以其他方式另外特别强调。术语“包括”、“包含”、“具有”及它们的变形都意味着“包括但不限于”，除非是以其他方式另外特别强调。
针对目前活动层厚度估算的精度和可靠性低的问题，本申请实施例提供了一种多年冻土活动层厚度估算方法，该方法通过对多个不同轨道的SAR影像数据集分别进行时序InSAR处理得到多个LOS向地表时序形变，然后利用多个LOS向地表时序形变获取估算区域的地表真实二维形变结果，从而有效避免了单轨InSAR形变观测值难以反映地表真实形变特征的问题，同时在利用地表真实二维形变结果进行活动层厚度估算时，全面考虑了实际冻土变化过程中未冻水的影响，从而使得本申请的估算方法能大大提高活动层厚度估算的精度和可靠性。
下面结合具体实施例对本申请提供的多年冻土活动层厚度估算方法进行示例性的说明。
本申请实施例提供了一种多年冻土活动层厚度估算方法，该方法可以由终端设备执行，也可以由应用于终端设备中的装置(比如芯片)来执行，下述实施例以该方法由终端设备执行为例。作为一种示例，该终端设备可以是平板，服务器或者笔记本电脑等，本申请实施例对此不做限定。
如图1所示，本申请实施例提供的多年冻土活动层厚度估算方法包括如下步骤：
步骤11，获取估算区域的多个不同轨道的SAR影像数据集。
上述估算区域为需进行多年冻土活动层厚度估算的区域；上述多个不同轨道的合成孔径雷达(SAR，SyntheticAperture Radar)影像数据集可以是合成孔径雷达卫星在多个不同轨道对估算区域进行观测得到的。
可以理解的是，为便于终端设备对估算区域进行多年冻土活动层厚度估算，合成孔径雷达卫星在得到上述SAR影像数据集后，可将得到的SAR影像数据集传输给终端设备。
步骤12，分别对多个SAR影像数据集进行时序InSAR处理，得到每个SAR影像数据集的LOS向地表时序形变。
在本申请的一些实施例中，计算第i个轨道的SAR影像数据集的LOS向地表时序形变的具体实现方式可以为：①对第i个轨道的SAR影像数据集的Ni(Ni第i个轨道的SAR影像数据集的总景数)景SAR影像进行配准；②选取合适的时空基线阈值，对满足阈值条件的干涉对采用外部数字高程模型(DEM)数据进行差分干涉处理得到Mi景差分干涉图；③对差分干涉图进行滤波、解缠，获得解缠后的差分干涉图；④通过合适的相干性阈值选出高相干点，并利用之前的Mi景解缠差分干涉图和短基线集合成孔径雷达干涉测量(SBAS-InSAR，smallbaseline subset InSAR)技术解算各高相干点的视线(LOS，Line OfSight)向地表时序形变
Figure BDA0003932525470000071
需要说明的是，由于时序InSAR处理为SAR影像数据集的常用处理方式，因此，在此不对时序InSAR处理的过程进行过多赘述。
可以理解的是，为方便后续的处理，可将不同轨道解算得到的高相干点LOS向地表时序形变结果地理编码、重采样至相同坐标系和分辨率。
步骤13，根据每个SAR影像数据集的LOS向地表时序形变，获取估算区域的地表真实二维形变结果。
在本申请的一些实施例中，为提高地表真实形变结果的准确性，可通过雷达成像的几何关系，结合多个不同轨道的LOS向地表时序形变，将一维LOS向形变分解为二维或三维地表真实形变，进而利用最小二乘法计算获取估算区域的地表真实二维形变结果。
步骤14，获取估算区域的遥感数据和土壤数据，并根据遥感数据和土壤数据对估算区域的未冻水含量进行预报。
在本申请的一些实施例中，可通过遥感数据采集设备(如卫星)采集得到上述遥感数据，该遥感数据包括估算区域的多层土壤的温度数据。
在本申请的一些实施例中，土壤数据包括上述多层土壤中每层土壤的土质类型以及每种土质类型对应的未冻水含量常数。
需要说明的是，上述土质类型可通过常用的土质识别方法对每层土壤的样本进行识别得到，上述未冻水含量常数为与未冻水含量相关的参数，具体可以通过实地钻孔测量的土壤温度和土壤水分数据结合“两点法”计算获取，更具的，可以通过构建土质类型与未冻水含量(即土壤水分数据)之间的关系式(该关系式包含未冻水含量常数)，然后在估算区域的多个位置钻孔，测量每个位置的土壤温度和土壤水分数据，再利用采集的数据代入关系式对未冻水含量常数进行求解，得到未冻水含量常数。
在本申请的一些实施例中，具体可利用公式Wu＝aθ-b计算估算区域的未冻水含量。其中，Wu表示估算区域的未冻水含量，θ表示温度数据的绝对值，a和b均表示未冻水含量常数。这里需要说明的是，由于反演过程都是逐像素反演的，也就是说在计算未冻水含量时，需针对每个像素点，利用该像素点对应的每一层土壤温度、未冻水含量常数进行计算，最终得到估算区域的未冻水含量。
步骤15，根据地表真实二维形变结果和未冻水含量，对估算区域的活动层厚度进行估算。
值得一提的是，在本申请的一些实施例中，通过对多个不同轨道的SAR影像数据集分别进行时序InSAR处理得到多个LOS向地表时序形变，然后利用多个LOS向地表时序形变获取估算区域的地表真实二维形变结果，从而有效避免了单轨InSAR形变观测值难以反映地表真实形变特征的问题，同时在利用地表真实二维形变结果进行活动层厚度估算时，全面考虑了实际冻土变化过程中未冻水的影响，从而使得本申请的估算方法能大大提高活动层厚度估算的精度和可靠性。
下面结合具体实施例对地表真实二维形变结果的获取进行示例性的说明。
在本申请的一些实施例中，上述步骤13，根据每个SAR影像数据集的LOS向地表时序形变，获取估算区域的地表真实二维形变结果的具体实现方式包括如下步骤：
步骤一，分别将每个SAR影像数据集的LOS向地表时序形变分解为垂直、东西和南北三个方向的形变分量，得到每个轨道的观测方程。
在本申请的一些实施例中，第i个轨道的观测方程为：
Figure BDA0003932525470000081
其中，
Figure BDA0003932525470000082
表示第i个轨道的SAR影像数据集的LOS向地表时序形变，dU表示
Figure BDA0003932525470000083
在垂直方向的形变分量，dE表示
Figure BDA0003932525470000084
在东西方向的形变分量，dN表示
Figure BDA0003932525470000085
在南北方向的形变分量，θi表示第i个轨道的SAR影像数据集的入射角，αi表示第i个轨道的SAR影像数据集的方位角，i＝1,…,j，j表示轨道的总数。
步骤二，利用多个轨道的观测方程获取估算区域的地表真实二维形变结果。
需要说明的是，要利用InSAR LOS向地表时序形变获取地表的真实三维形变结果，通常需要利用三个及以上独立的成像几何获取的多个测量值，从一维形变拓展至三维。考虑到差分干涉测量(D-InSAR)对南北向形变不敏感，为了保证最小二乘估算的正确性，通常忽略南北向形变，结合多个轨道的LOS向地表时序形变，解算估算区域地表真实二维形变结果。
具体的，可在忽略南北方向形变分量的基础上，利用最小二乘法对
Figure BDA0003932525470000091
Figure BDA0003932525470000092
进行求解，得到估算区域的地表真实二维形变结果，该地表真实二维形变结果包括dU和dE。
在本申请的一些实施例中，在得到估算区域的地表真实二维形变结果后，可利用傅里叶拟合方法对地表真实二维形变结果中的dU进行拟合，并从拟合结果中提取融化季时段对应的垂直方向的形变分量
Figure BDA0003932525470000093
以便对活动层厚度进行估算。
需要说明的是，傅里叶拟合方法对dU进行拟合，得到的拟合结果可以为一曲线，该曲线的横轴为时间，纵轴为dU的值，因此能从拟合结果中简单、快速提取出融化季时段对应的垂直方向的形变分量
Figure BDA0003932525470000094
在本申请的一些实施例中，在得到地表真实二维形变结果和未冻水含量后，可通过公式
Figure BDA0003932525470000095
计算估算区域从活动层融化开始t0至结束时刻tk的活动层厚度H。
其中，ΔHr表示估算区域在tr-1至tr时刻融化的活动层厚度，由于地表沉降完全是由活动层中冰水相变所造成的体积差引起的，可以建立ΔDr与ΔHr之间的关系为：
Figure BDA00039325254700000911
Figure BDA0003932525470000097
表示
Figure BDA0003932525470000098
在tr时刻的值，
Figure BDA0003932525470000099
表示
Figure BDA00039325254700000910
在tr-1时刻的值，Wice表示估算区域在融化季时段的tr-1至tr时刻参与水冰相变过程的冰含量。
假设融化季估算区域内多年冻土中的水处于饱和状态，且地表沉降完全是由于活动层中冰与水的相变所造成的体积差引起的。考虑到tr-1时刻冻结的活动层中仍然存在一定量的未冻水
Figure BDA0003932525470000101
结合水质量守恒，可以得到参与水冰相变的冰含量为：
Figure BDA0003932525470000102
ρwater表示水的密度，ρice表示冰的密度，Wr表示估算区域在tr时刻的土壤水分含量，可通过遥感数据采集设备(如卫星)采集得到，
Figure BDA0003932525470000103
表示估算区域在tr-1时刻的未冻水含量，即Wu在tr-1时刻的值。
由上可见，上述估算方法通过充分利用多轨道InSAR形变观测值，避免了单轨InSAR形变观测值难以反映地表真实形变特征的现象，并且在对冻土形变与活动层变化建模时全面考虑了实际冻土变化过程中未冻水的影响，同时方法易于实施，应用广泛，有助于提高大范围、高分辨率的InSAR活动层厚度估计的精度和可靠性。
为便于理解上述估算方法，在此以一具体实例对上述估算方法进行示例性说明。
在该实例中，采集的数据包括：覆盖了青藏高原北麓河多年冻土区(覆盖范围：92.7801°E-93.9707°E，34.5721°N-35.4318°N)的哨兵一号(Sentinel-1)卫星升轨和降轨影像数据和土壤水分主动-被动(SMAP)遥感土壤温湿度数据。
在该实例中，估算区域包含了楚玛尔河高平原、可可西里山和北麓河盆地；整个区域大多为富冰多年冻土，地下冰含量丰富，最大体积含冰量达到了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影像一致。
在该实例中，对估算区域的活动层厚度进行估算的具体步骤如下：
步骤1：对Sentinel-1升降和降轨影像分别进行配准，配准精度优于1/1000像素；选取时空基线阈值分别为150m和36天，对满足该阈值的干涉对采用航天飞机雷达地形测绘使命(SRTM)DEM数据进行差分干涉处理得到一系列的差分干涉图；采用Goldstein滤波算法(Goldstein滤波算法是一种干涉图滤波算法)对差分干涉图进行二维频率域滤波，并利用最小网络费用流方法进行二维相位解缠，为了避免低相干点对解缠的影响，在解缠过程中去除了相干性低于0.2的像素点，最后获取了解缠后的差分干涉图；选取最小相干性和平均相干性分别不低于0.35和0.85的像素点为高相干点，通过SBAS-InSAR技术逐像素求解出这些点的在升轨和降轨LOS向上的地表时序形变
Figure BDA0003932525470000111
和
Figure BDA0003932525470000112
将上述得到的地表时序形变结果编码、重采样到相同的地理坐标系和分辨率下。
步骤2：在忽略南北向形变的基础上，根据雷达成像的几何关系，将已获取的LOS向地表时序形变
Figure BDA0003932525470000113
和
Figure BDA0003932525470000114
视为观测值，可建立如下观测方程：
Figure BDA0003932525470000115
其中，θasc表示
Figure BDA0003932525470000116
对应的影像数据的入射角，αasc表示
Figure BDA0003932525470000117
对应的影像数据的方位角，θdsc表示
Figure BDA0003932525470000118
对应的影像数据的入射角，αdsc表示
Figure BDA0003932525470000119
对应的影像数据的方位角。根据最小二乘平差原理即可对上述观测方程进行求解得到地表真实二维时序形变dU和dE。考虑到观测方程可能存在不适定问题，此时可通过吉洪诺夫(Tikhonov)正则化的方法对其进行求解。
Figure BDA00039325254700001110
其中，B为公式中的系数矩阵；P为权阵，在本实例中取单位权阵；[t0 t1 … tN]为SAR影像的获取时刻。
获取各高相干点在垂直向上的时序形变dU之后，考虑到冻土形变一般由气候变暖引起的长期线性形变和冻融循环有关的周期性形变两部分组成，首先对各高相干点在垂直向上的时序形变dU进行去线性趋势，然后利用公式傅里叶拟合法对每个高相干点的垂直向时序形变进行参数拟合，确定每年融化季的时间(如图2和图3所示)并提取处于融化季之中的时序形变
Figure BDA00039325254700001111
其中，图2中的横轴为融化季开始时间，图3中的横轴为融化季结束时间，图2和图3中纵轴为确定融化季开始和结束时间的频率。
步骤3：根据估算区域周边的地质资料，可以获取估算区域的土壤类型和与未冻水含量相关的经验参数(即前文的未冻水含量常数a和b)，具体参数见
表1。
深度/m 土壤类型 未冻水含量常数a 未冻水含量常数b
0-2 砂土 0.07 0.17
2-10 黏土 0.12 0.15
>10 砂岩 0.01 0.10
表1
根据步骤2获取的融化季时间，可以从SMAP土壤温度数据中获取各土壤层在融化前的温度，将其与表1数据代入公式Wu＝aθ-b中，即可实现未冻水含量的预报，计算得到各土壤层在融化前的未冻水含量
Figure BDA0003932525470000123
表示Wu在ti时刻的值。
步骤4：将步骤2得到的融化季垂直向地表时序形变
Figure BDA0003932525470000124
步骤3得到的各土壤层融化前未冻水含量
Figure BDA0003932525470000125
和SMAP分层土壤水分数据代入公式
Figure BDA0003932525470000121
Figure BDA0003932525470000122
中，得到活动层估算结果。
由图4可以看出，通过与估算区域三个多年冻土监测站点实测的活动层厚度进行对比，本申请提供的估算方法均最接近与实测值，这说明本申请提供的估算方法估算出的活动层厚度更为准确，也验证了本申请提供的估算方法能够有效的提高大范围、高分辨率的InSAR活动层厚度估计的精度和可靠性。其中，图4中WD3、HR3和QT08分别表示3个不同的多年冻土监测站点，实测值指的是对应站点的活动层厚度的实际值，升轨和降轨分别表示利用单轨(升轨/降轨)SAR数据进行活动层厚度预测的方法。
以上所述是本申请的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本申请所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本申请的保护范围。

Claims (7)

1.一种多年冻土活动层厚度估算方法，其特征在于，包括：
获取估算区域的多个不同轨道的SAR影像数据集；
分别对多个所述SAR影像数据集进行时序InSAR处理，得到每个所述SAR影像数据集的LOS向地表时序形变；
根据每个所述SAR影像数据集的LOS向地表时序形变，获取所述估算区域的地表真实二维形变结果；
获取所述估算区域的遥感数据和土壤数据，并根据所述遥感数据和所述土壤数据对所述估算区域的未冻水含量进行预报；
根据所述地表真实二维形变结果和所述未冻水含量，对所述估算区域的活动层厚度进行估算。
2.根据权利要求1所述的方法，其特征在于，所述根据每个所述SAR影像数据集的LOS向地表时序形变，获取所述估算区域的地表真实二维形变结果，包括：
分别将每个所述SAR影像数据集的LOS向地表时序形变分解为垂直、东西和南北三个方向的形变分量，得到每个轨道的观测方程；
利用多个轨道的观测方程获取所述估算区域的地表真实二维形变结果。
3.根据权利要求2所述的方法，其特征在于，第i个轨道的观测方程为：
Figure FDA0003932525460000011
其中，
Figure FDA0003932525460000012
表示第i个轨道的SAR影像数据集的LOS向地表时序形变，dU表示
Figure FDA0003932525460000013
在垂直方向的形变分量，dE表示
Figure FDA0003932525460000014
在东西方向的形变分量，dN表示
Figure FDA0003932525460000015
在南北方向的形变分量，θi表示第i个轨道的SAR影像数据集的入射角，αi表示第i个轨道的SAR影像数据集的方位角，i＝1,…,j，j表示轨道的总数。
4.根据权利要求3所述的方法，其特征在于，所述利用多个轨道的观测方程获取所述估算区域的地表真实二维形变结果，包括：
利用最小二乘法对
Figure FDA0003932525460000021
进行求解，得到所述估算区域的地表真实二维形变结果；所述地表真实二维形变结果包括dU和dE。
5.根据权利要求1所述的方法，其特征在于，所述遥感数据包括多层土壤的温度数据，所述土壤数据包括每层所述土壤的土质类型以及每种土质类型对应的未冻水含量常数；
所述根据所述遥感数据和所述土壤数据对所述估算区域的未冻水含量进行预报，包括：
利用公式Wu＝aθ-b计算所述估算区域的未冻水含量；
其中，Wu表示所述估算区域的未冻水含量，θ表示所述温度数据的绝对值，a和b均表示所述未冻水含量常数。
6.根据权利要求4所述的方法，其特征在于，在所述根据每个所述SAR影像数据集的LOS向地表时序形变，获取所述估算区域的地表真实二维形变结果之后，所述方法还包括：
利用傅里叶拟合方法对所述地表真实二维形变结果中的dU进行拟合；
从拟合结果中提取融化季时段对应的垂直方向的形变分量
Figure FDA0003932525460000022
7.根据权利要求6所述的方法，其特征在于，所述根据所述地表真实二维形变结果和所述未冻水含量，对所述估算区域的活动层厚度进行估算，包括：
通过公式
Figure FDA0003932525460000023
计算所述估算区域从活动层融化开始t0至结束时刻tk的活动层厚度H；
其中，ΔHr表示所述估算区域在tr-1至tr时刻融化的活动层厚度，
Figure FDA0003932525460000024
Figure FDA0003932525460000025
Figure FDA0003932525460000026
表示
Figure FDA0003932525460000027
在tr时刻的值，
Figure FDA0003932525460000028
表示
Figure FDA0003932525460000029
在tr-1时刻的值，Wice表示所述估算区域在融化季时段的tr-1至tr时刻参与水冰相变过程的冰含量，
Figure FDA00039325254600000210
ρwater表示水的密度，ρice表示冰的密度，Wr表示所述估算区域在tr时刻的土壤水分含量，
Figure FDA0003932525460000031
表示所述估算区域在tr-1时刻的未冻水含量。
CN202211392436.0A 2022-11-08 2022-11-08 联合多轨道InSAR形变观测和未冻水含量的活动层厚度估算方法 Active CN115540788B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211392436.0A CN115540788B (zh) 2022-11-08 2022-11-08 联合多轨道InSAR形变观测和未冻水含量的活动层厚度估算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211392436.0A CN115540788B (zh) 2022-11-08 2022-11-08 联合多轨道InSAR形变观测和未冻水含量的活动层厚度估算方法

Publications (2)

Publication Number Publication Date
CN115540788A true CN115540788A (zh) 2022-12-30
CN115540788B CN115540788B (zh) 2023-08-29

Family

ID=84720947

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211392436.0A Active CN115540788B (zh) 2022-11-08 2022-11-08 联合多轨道InSAR形变观测和未冻水含量的活动层厚度估算方法

Country Status (1)

Country Link
CN (1) CN115540788B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116051620A (zh) * 2023-04-03 2023-05-02 之江实验室 基于InSAR技术冻土区活动层厚度估计方法和系统

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140062764A1 (en) * 2012-09-04 2014-03-06 Fugro Earthdata, Inc. Method and apparatus for mapping and characterizing sea ice from airborne simultaneous dual frequency interferometric synthetic aperture radar (ifsar) measurements
CN104792818A (zh) * 2015-04-07 2015-07-22 河南大学 将土中水相变潜热进行能量替代的黏土冻结阶段比热计算方法
CN109165463A (zh) * 2018-09-12 2019-01-08 中国科学院寒区旱区环境与工程研究所 多年冻土活动层厚度的遥感估算方法、装置及可读存储介质
CN112986981A (zh) * 2021-02-18 2021-06-18 中国科学院西北生态环境资源研究院 多年冻土区地表冻融形变监测方法、装置和电子设备
CN113111531A (zh) * 2021-04-23 2021-07-13 中国水利水电科学研究院 面向分布式水文模型的季节性冻土区冻土层厚度模拟方法
CN113419044A (zh) * 2021-06-02 2021-09-21 中国科学院西北生态环境资源研究院 基于黏土扩散层离子浓度梯度的冻土未冻水含量计算方法
CN114563787A (zh) * 2022-02-25 2022-05-31 中山大学 一种冻土形变监测方法、装置、设备及存储介质
CN114707344A (zh) * 2022-04-15 2022-07-05 西南交通大学 一种基于系统动力学的多年冻土活动层厚度计算方法
CN114966692A (zh) * 2022-07-19 2022-08-30 之江实验室 基于Transformer的InSAR技术冻土区多变量时序形变预测方法及装置

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140062764A1 (en) * 2012-09-04 2014-03-06 Fugro Earthdata, Inc. Method and apparatus for mapping and characterizing sea ice from airborne simultaneous dual frequency interferometric synthetic aperture radar (ifsar) measurements
CN104792818A (zh) * 2015-04-07 2015-07-22 河南大学 将土中水相变潜热进行能量替代的黏土冻结阶段比热计算方法
CN109165463A (zh) * 2018-09-12 2019-01-08 中国科学院寒区旱区环境与工程研究所 多年冻土活动层厚度的遥感估算方法、装置及可读存储介质
CN112986981A (zh) * 2021-02-18 2021-06-18 中国科学院西北生态环境资源研究院 多年冻土区地表冻融形变监测方法、装置和电子设备
CN113111531A (zh) * 2021-04-23 2021-07-13 中国水利水电科学研究院 面向分布式水文模型的季节性冻土区冻土层厚度模拟方法
CN113419044A (zh) * 2021-06-02 2021-09-21 中国科学院西北生态环境资源研究院 基于黏土扩散层离子浓度梯度的冻土未冻水含量计算方法
CN114563787A (zh) * 2022-02-25 2022-05-31 中山大学 一种冻土形变监测方法、装置、设备及存储介质
CN114707344A (zh) * 2022-04-15 2022-07-05 西南交通大学 一种基于系统动力学的多年冻土活动层厚度计算方法
CN114966692A (zh) * 2022-07-19 2022-08-30 之江实验室 基于Transformer的InSAR技术冻土区多变量时序形变预测方法及装置

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
黄兆欢: "基于真实二维InSAR技术的北麓河多年冻土区地表形变时空特征分析与活动层厚度反演", 《中国优秀硕士学位论文全文数据库基础科学辑》, no. 09, pages 7 - 8 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116051620A (zh) * 2023-04-03 2023-05-02 之江实验室 基于InSAR技术冻土区活动层厚度估计方法和系统

Also Published As

Publication number Publication date
CN115540788B (zh) 2023-08-29

Similar Documents

Publication Publication Date Title
Necsoiu et al. Rock glacier dynamics in Southern Carpathian Mountains from high-resolution optical and multi-temporal SAR satellite imagery
Cao et al. Spatial variability of active layer thickness detected by ground‐penetrating radar in the Qilian Mountains, Western China
Liu et al. Surface motion of active rock glaciers in the Sierra Nevada, California, USA: inventory and a case study using InSAR
Minet et al. Mapping shallow soil moisture profiles at the field scale using full-waveform inversion of ground penetrating radar data
Strozzi et al. Detecting and quantifying mountain permafrost creep from in situ inventory, space-borne radar interferometry and airborne digital photogrammetry
Belart et al. Winter mass balance of Drangajökull ice cap (NW Iceland) derived from satellite sub-meter stereo images
Li et al. Investigating mountain glacier motion with the method of SAR intensity-tracking: Removal of topographic effects and analysis of the dynamic patterns
Samsonov et al. SAR-derived flow velocity and its link to glacier surface elevation change and mass balance
Li et al. Anomalous glacier changes in the southeast of Tuomuer‐khan Tengri Mountain ranges, central Tianshan
CN116051620B (zh) 基于InSAR技术冻土区活动层厚度估计方法和系统
Buckel et al. Insights into a remote cryosphere: a multi-method approach to assess permafrost occurrence at the Qugaqie basin, western Nyainqêntanglha Range, Tibetan Plateau
CN115540788B (zh) 联合多轨道InSAR形变观测和未冻水含量的活动层厚度估算方法
Chen et al. Spatial variability in melting on Himalayan debris-covered glaciers from 2000 to 2013
Boniface et al. Potential of shipborne GPS atmospheric delay data for prediction of Mediterranean intense weather events
Carrara et al. Post-emplacement dynamics of andesitic lava flows at Volcán de Colima, Mexico, revealed by radar and optical remote sensing data
Wu et al. Surface-deformation monitoring in the permafrost regions over the Tibetan Plateau, using Sentinel-1 data
Frey et al. Tomographic profiling with snowscat within the ESA snowlab campaign: Time series of snow profiles over three snow seasons
Zhao et al. Investigation, Monitoring, and Simulation of Permafrost on the Qinghai‐Tibet Plateau: A Review
Zhang et al. Glacier elevation change in the Western Qilian mountains as observed by TerraSAR-X/TanDEM-X images
Evans et al. Determination of snow depth using elevation differences determined by interferometric SAR (InSAR)
CN115061112A (zh) 面状雪水当量获取方法、装置及电子设备
Li et al. Monitoring snow depth and its change using repeat-pass interferometric SAR in Manas River Basin
Ho et al. Measuring ground subsidence in Hanoi city by radar interferometry
Xu et al. High-precision measurements of the inter-annual evolution for Urumqi Glacier No. 1 in eastern Tien Shan, China
Thomas Remote sensing reveals shrinking Greenland ice sheet

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant