多年冻土活动层厚度的遥感估算方法、装置及可读存储介质

摘要

本发明提供了一种多年冻土活动层厚度的遥感估算方法、装置及可读存储介质，所述方法包括：获取被估算区域的遥感数据及辅助数据；根据预设的多年冻土活动层厚度的计算公式、遥感数据及辅助数据估算所述被估算区域的多年冻土活动层厚度；其中，所述多年冻土活动层厚度的计算公式为：本发明实施例根据预设的多年冻土活动层厚度的计算公式及相关数据，如地表温度融化指数、叶面积指数、积雪日数等遥感数据以及土壤性质等，并获取多年冻土活动层厚度的测量数据，通过估计所述计算公式中的参数，估算多年冻土活动层厚度；通过遥感观测等获取的数据应用于多年冻土活动层厚度监测当中，可以在时间及空间上连续地监测多年冻土活动层厚度的变化。

说明书

技术领域

本发明涉及遥感技术领域，具体而言，涉及一种多年冻土活动层厚度的遥感估算方法、装置及可读存储介质。

背景技术

目前，多年冻土活动层厚度观测主要依赖于地面测量，主要测量手段有冻融管或钻孔观测，其空间代表性差，费时费力，成本高且时效性不强，且难以获取较大范围空间连续分布的活动层厚度信息。而在空间上估计活动层厚度多以斯蒂芬公式和库德里亚夫采夫经验公式为主，这些经验模型依赖于气象站的气温或土壤表面温度观测，由于多年冻土区气象站极其稀少，这严重限制了对多年冻土活动层厚度的监测。

发明内容

为了克服上述现有技术中的不足，本发明提供了一种多年冻土活动层厚度的遥感估算方法、装置及可读存储介质。

为了实现上述目的，本发明实施例所提供的技术方案如下所示：

第一方面，本发明实施例提供一种多年冻土活动层厚度的遥感估算方法，用于对多年冻土活动层的厚度进行估算，所述遥感估算方法包括：获取被估算区域的遥感数据和辅助数据，所述遥感数据包括：地表温度融化指数DDT、年平均叶面积指数LAI和积雪日数SnowD，辅助数据包括相关土壤数据，所述相关土壤数据的类别包括：土壤沙粒含量Sand、土壤粘粒含量Clay、土壤粉粒含量Silt、土壤有机碳含量Org和土壤密度Bld；

根据预设的多年冻土活动层厚度的计算公式、所述遥感数据及所述辅助数据估算所述被估算区域的多年冻土活动层厚度；其中，所述多年冻土活动层厚度的计算公式为：

其中，a、b、c、d、e、f、g、h和T为预先获取到的估算参数；ALT为所述被估算区域的多年冻土活动层厚度。

结合第一方面，在一些可能的实现方式中，在根据预设的多年冻土活动层厚度的计算公式和所述相关数据估算所述被估算区域的多年冻土活动层厚度之前，所述估算方法还包括：

预先获取所述被估算区域的所述遥感数据及辅助数据的多组数据以及所述多组数据对应的所述被估算区域的多年冻土活动层厚度、以及所述计算公式；

根据所述多组数据、所述多组数据对应的所述多年冻土活动层厚度及所述多年冻土活动层厚度的计算公式，通过参数估计方法得到所述估算参数。

结合第一方面，在一些可能的实现方式中，获取地表温度融化指数DDT的步骤包括：

对获取的一日四次的白天和夜晚地表温度，根据与所述白天和夜晚地表温度一起获取的相应质量控制符QC，选择质量最好的进行白天和夜晚地表温度的配对，对配对成功的像元直接进行平均，得到该像元的日平均地表温度；

提取配对成功像元的日地表温度变幅，利用克里金方法插值到每个像元，得到全年每日全部像元的地表温度变幅D；

对于没有形成配对的像元，根据当天全部像元的所述地表温度变幅D，计算该天没有形成配对像元的日平均地表温度；

根据前三步获取的日平均地表温度，利用时空增强方法估计完全无值像元的日平均地表温度值，得到时空完整的日平均地表温度；

根据时空完整的每日平均地表温度，确定所述一年中的大于0℃的每日平均地表温度的代数和为地表温度融化指数DDT。

结合第一方面，在一些可能的实现方式中，获取年平均叶面积指数LAI的步骤包括：

判断是否存在时间上缺失的被估算区域的叶面积指数；

根据时间序列重建方法对所述缺失的叶面积指数进行数据重建，得到全年完整的叶面积指数；

计算所述全年完整的叶面积指数的平均值作为年平均叶面积指数LAI。

结合第一方面，在一些可能的实现方式中，获取积雪日数SnowD的步骤包括：

判断所获取的积雪遥感产品是否存在时间上的缺失；

根据时间序列重建方法对所述缺失的积雪遥感产品进行数据重建，得到全年完整的积雪数据；

根据全年完整的积雪数据得到所述积雪日数SnowD。

结合第一方面，在一些可能的实现方式中，获取所述相关土壤数据包括：

根据活动层的深度进行分区，获取各种参数的多层土壤数据；

计算各种参数的多层土壤数据的平均值，得到各个土壤数据平均后的多层平均土壤数据；

确定所述多层平均土壤数据为所述相关土壤数据。

第二方面，本发明实施例还提供一种基于遥感的多年冻土活动层厚度的遥感估算装置，用于对多年冻土活动层的厚度进行估算；所述遥感估算装置包括：获取单元与处理单元；

获取单元用于获取被估算区域的遥感数据和辅助数据，所述遥感数据包括：地表温度融化指数DDT、年平均叶面积指数LAI和积雪日数SnowD，所述辅助数据包括相关土壤数据，所述相关土壤数据的类别包括：土壤沙粒含量Sand、土壤粘粒含量Clay、土壤粉粒含量Silt、土壤有机碳含量Org和土壤密度Bld；

处理单元用于根据预设的多年冻土活动层厚度的计算公式、所述遥感数据及所述辅助数据估算所述被估算区域的多年冻土活动层厚度；其中，所述多年冻土活动层厚度的计算公式为：

其中，a、b、c、d、e、f、g、h和T为预先获取到的估算参数；ALT为所述被估算区域的多年冻土活动层厚度。

结合第二方面，在一些可能的实现方式中，所述获取单元还用于预先获取所述被估算区域的所述遥感数据及辅助数据的多组数据以及所述多组数据对应的所述被估算区域的多年冻土活动层厚度、以及所述计算公式。

结合第二方面，在一些可能的实现方式中，所述处理单元还用于根据所述多组数据、所述多组数据对应的所述多年冻土活动层厚度及所述多年冻土活动层厚度的计算公式，通过参数估计方法得到所述估算参数。

第三方面，本发明实施例还提供一种可读存储介质，所述可读存储介质中存储有计算机程序，当所述计算机程序在计算机上运行时，使得所述计算机执行第一方面或第一方面中任一实施例所述的基于遥感的估算方法。

本发明的有益效果包括：

本发明实施例中的估算方法通过获取冻土活动层厚度的相关数据，利用多年冻土活动层厚度的计算公式(例如是一个线性估算模型)，估算得到被估算区域的多年冻土活动层厚度；直接通过遥感观测的方式获取地表温度融化指数DDT、年平均叶面积指数LAI和积雪日数SnowD，并结合一些易于获得的辅助数据，利用所述计算公式，便可以估算得到被估算区域的多年冻土活动层厚度，不再需要冻融管或钻孔观测等复杂且需人工操作的方式；并且可以在时间和空间上连续地监测被估算区域多年冻土活动层厚度。

为使本发明的上述目的、特征和优点能更明显易懂，下文举出本发明实施例，并配合所附附图，作详细说明。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明实施例1提供的多年冻土活动层厚度的遥感估算方法的流程示意图；

图2为本发明实施例1提供的计算得到估算参数的流程示意图；

图3为本发明实施例1提供的得到估算参数的举例示意图；

图4为本发明实施例2提供的一种多年冻土活动层厚度的遥感估算装置的功能框图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

在本发明实施例的描述中，需要说明的是，指示方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，或者是本领域技术人员惯常理解的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

实施例1

请参考图1，图1为本发明实施例1提供的基于遥感的多年冻土活动层厚度的估算方法的流程示意图。所述估算方法包括：

步骤S101：获取被估算区域的遥感数据和辅助数据；

步骤S102：根据预设的多年冻土活动层厚度的计算公式、所述遥感数据及所述辅助数据估算所述被估算区域的多年冻土活动层厚度。

需要说明的是，所述遥感数据包括：地表温度融化指数DDT、年平均叶面积指数LAI和积雪日数SnowD，所述辅助数据包括相关土壤数据，所述相关土壤数据的类别包括：土壤沙粒含量Sand、土壤粘粒含量Clay、土壤粉粒含量Silt、土壤有机碳含量Org和土壤密度Bld；所述多年冻土活动层厚度的计算公式为：

其中，a、b、c、d、e、f、g、h和T为预先获取到的估算参数；ALT为所述被估算区域的多年冻土活动层厚度。

通过获取被估算区域的遥感数据、辅助数据及针对被估算区域的确定的估算参数，并通过多年冻土活动层厚度的计算公式，计算得到被估算区域的多年冻土活动层厚度。在所述某一估算区域估算参数确定的情况下，只需要通过遥感观测获取地表温度融化指数DDT、年平均叶面积指数LAI和积雪日数SnowD，结合一些可轻易得到的辅助数据，利用所述多年冻土活动层厚度的计算公式便可估算得到被估算区域的多年冻土活动层厚度，避免了现有需实地测量数据的复杂性及气象站不足的问题；还可以实现在时间和空间上连续地监测被估算区域多年冻土活动层厚度变化。

为了便于更清楚地描述本实施例中的估算方法的实施过程，下面将先介绍步骤S102中估算参数可以采用的获取过程。

请参考图2，图2为本发明实施例1提供的计算得到估算参数的流程示意图。

可选地，在根据预设的多年冻土活动层厚度的计算公式、所述遥感数据及所述辅助数据估算所述被估算区域的多年冻土活动层厚度之前，所述估算方法还包括：

步骤S103：预先获取所述被估算区域的所述遥感数据及所述副主数据的多组数据以及所述多组数据对应的所述被估算区域的多年冻土活动层厚度、以及所述计算公式；

步骤S104：根据所述多组数据、所述多组数据对应的所述多年冻土活动层厚度及所述多年冻土活动层厚度的计算公式，通过参数估计方法得到所述估算参数。

需要说明的是，本发明实施例采用的参数估计方法为最小二乘法，在其他实施例中，还可以采用其他参数估计方法，如最大似然法、最大后验估计等。

步骤S103负责获取所述计算公式中除了估算参数外的其他一整套参数，作为样本；进而通过步骤S104，根据现有的带有未知估算参数的计算公式，利用现有的多年冻土活动层厚度观测数据做最小二乘法估计，通过逼近的方式优化估算参数，进而使得计算公式无限靠近现有观测数据，进而得到尽可能接近真实厚度值的计算公式。接下来对上述估算参数的获取举例进行说明。

请参考图3，图3为本发明实施例1提供的参数估计样本得到估算参数的举例示意图。假设被估算区域为青藏高原沿线多年冻土区，本例根据公开的青藏高原沿线多年冻土区的活动层厚度观测数据，在选择51个采样点，相关系数达到0.98，显著水平P<0.0001时，通过参数估计方法得到图中所示结果，即得到针对青藏高原沿线多年冻土区的活动层厚度的估算公式为：

即获取到所述估算区域的估算参数，即a＝0.97853，b＝22.25079，c＝0.10733,d＝60.80734,e＝33.18032,f＝98.56106,g＝2.44596,h＝-13.93537,T＝-7512.3858。

接下来对步骤S101中提到的相关数据的可能的获取方式进行介绍。

可选地，获取地表温度融化指数DDT的步骤包括：

对获取的一日四次的白天和夜晚地表温度，根据与所述白天和夜晚地表温度一起获取的相应质量控制符QC，选择质量最好的进行白天和夜晚地表温度的配对，对配对成功的像元直接进行平均，得到该像元的日平均地表温度；

提取配对成功像元的日地表温度变幅，利用克里金方法插值到每个像元，得到全年每日全部像元的地表温度变幅D；

对于没有形成配对的像元，根据当天全部像元的所述地表温度变幅D，计算该天没有形成配对像元的日平均地表温度；

根据前三步获取的日平均地表温度，利用时空增强方法估计完全无值像元的日平均地表温度值，得到时空完整的日平均地表温度；

根据时空完整的每日平均地表温度，确定所述一年中的大于0℃的每日平均地表温度的代数和为地表温度融化指数DDT。

需要说明的是，本发明实施例所述获取一日四次白天和夜晚地表温度的方式是通过搭载在Terra和aqua卫星上的Modis传感器返回的数据直接获取，具体而言，返回的数据指即Modis传感器返回的产品：MOD11A1和MYD11A1中的数据。在其他实施例中，也可以通过搭载在NOAA系列气象卫星上的AVHRR传感器获取，本发明对此不作限制。

由于卫星传感器时常会受到云的影响，时常缺失部分日期的观测数据，或是某些日期的观测数据失真，通过上述方法能够获取到时空完整的日平均地表温度，进而获取优质的地表温度融化指数DDT，最终提高了估算精度。

通过充分利用MODIS一天四次地表温度观测，区分白天晚上数据的配对、未配对和完全无值三种情况以年为单位实现了全年日平均地表温度的恢复和重建。

需要说明的是，上述获取地表温度融化指数DDT的方式仅是一种实现获取地表温度融化指数DDT的可能的方式，在其他实施例中，为实现本发明的目的采用其他获取地表温度融化指数DDT的方式不应成为对本发明的限制。

可选地，获取年平均叶面积指数LAI的步骤包括：

判断是否存在时间上缺失的被估算区域的叶面积指数；

根据时间序列重建方法对所述缺失的叶面积指数进行数据重建，得到全年完整的叶面积指数；

计算所述全年完整的叶面积指数的平均值作为年平均叶面积指数LAI。

需要说明的是，叶面积指数LAI指植物叶片总面积与土地面积的比值；本发明实施例中叶面积指数的获取采用卫星遥感方法，具体为统计模型法，如归一化植被指数NDVI(英文：Normalized Difference Vegetation Index，简称NDVI)方法等，本发明对此不做限制。

同样由于存在云的影响，单纯获取的叶面积指数会存在缺失的问题，影响最终计算结果的准确度，采用时间序列重建方法能减轻影响准确度的程度。

所述时间序列重建方法采用常用的时间序列重建方法即可，例如阈值法、曲线拟合法及最大值合成法等。

需要说明的是，上述获取叶面积指数LAI的方式仅是一种获取叶面积指数LAI的可能的方式，在其他实施例中，为实现本发明的目的采用其他叶面积指数LAI的方式不应成为对本发明的限制。

可选地，获取积雪日数SnowD的步骤包括：

判断所获取的积雪遥感产品是否存在时间上的缺失；

根据时间序列重建方法对所述缺失的积雪遥感产品进行数据重建，得到全年完整的积雪数据；

根据全年完整的积雪数据得到所述积雪日数SnowD。

需要说明的是，积雪日数指一年中有积雪的日子数。本发明实施例采用的获取积雪遥感产品为结合Terra和Aqua卫星的积雪产品，进而获取逐日积雪覆盖数据；在其他实施例中，也可以采用其他积雪数据产品，本发明对此不做限制。所述时间序列重建方法同样采用常用的时间序列重建方法即可，例如阈值法、曲线拟合法及最大值合成法等。

可选地，获取相关土壤数据包括：

根据活动层的深度进行分区，获取各种参数的多层土壤数据；

计算各种参数的多层土壤数据的平均值，得到各个土壤数据平均后的多层平均土壤数据；

确定所述多层平均土壤数据为所述相关土壤数据。

需要说明的是，本发明实施例中所述相关土壤数据中各类土壤数据取自全球土壤信息设施(GSIF：Global Soil Information Facilities)项目制备的全球SoilGrid土壤属性数据。在其他实施例中，也可以通过如中国土壤数据库等获取，本发明对此不做限制。

以下将举例对获取相关土壤数据进行说明：在本实施例中，将所述冻土活动层以5cm为间距分为5cm深度层、10cm深度层及15cm深度层。根据上述的全球SoilGrid土壤属性数据，将三个层中的各种土壤数据分别平均，得到了所述的多层平均土壤数据，提高了估算的精度。

以土壤沙粒Sand含量为例，通过获取5cm深度层的土壤沙粒含量Sand₁、10cm深度层的土壤沙粒含量Sand₂及15cm深度层的土壤沙粒含量Sand₃，将三者做平均得到多层平均土壤沙粒含量Sand，即Sand＝(Sand₁+Sand₂+Sand₃)/3。

实施例2

请参考图4，图4为本发明实施例2提供的一种多年冻土活动层厚度的遥感估算装置的功能框图。所述遥感估算装置用于对多年冻土活动层的厚度进行估算，包括：获取单元201及处理单元202。

获取单元201用于获取被估算区域的遥感数据和辅助数据，所述遥感数据包括：地表温度融化指数DDT、年平均叶面积指数LAI和积雪日数SnowD，所述辅助数据包括相关土壤数据，所述相关土壤数据的类别包括：土壤沙粒含量Sand、土壤粘粒含量Clay、土壤粉粒含量Silt、土壤有机碳含量Org和土壤密度Bld；

处理单元202用于根据预设的多年冻土活动层厚度的计算公式、所述遥感数据及所述辅助数据估算所述被估算区域的多年冻土活动层厚度；其中，所述多年冻土活动层厚度的计算公式为：

其中，a、b、c、d、e、f、g、h和T为预先获取到的估算参数；ALT为所述被估算区域的多年冻土活动层厚度。

可选地，所述获取单元还用于预先获取所述被估算区域的所述遥感数据及辅助数据的多组数据以及所述多组数据对应的所述被估算区域的多年冻土活动层厚度、以及所述计算公式。

可选地，所述处理单元还用于根据所述多组数据、所述多组数据对应的所述多年冻土活动层厚度及所述多年冻土活动层厚度的计算公式，通过参数估计方法得到所述估算参数。

获取单元201获取上述相关数据的可以实现的方式如下：

可选地，获取单元201获取地表温度融化指数DDT的步骤包括：

对获取的一日四次的白天和夜晚地表温度，根据与所述白天和夜晚地表温度一起获取的相应质量控制符QC，选择质量最好的进行白天和夜晚地表温度的配对，对配对成功的像元直接进行平均，得到该像元的日平均地表温度；

提取配对成功像元的日地表温度变幅，利用克里金方法插值到每个像元，得到全年每日全部像元的地表温度变幅D；

对于没有形成配对的像元，根据当天全部像元的所述地表温度变幅D，计算该天没有形成配对像元的日平均地表温度；

根据前三步获取的日平均地表温度，利用时空增强方法估计完全无值像元的日平均地表温度值，得到时空完整的日平均地表温度；

根据时空完整的每日平均地表温度，确定所述一年中的大于0℃的每日平均地表温度的代数和为地表温度融化指数DDT。

可选地，获取单元201获取年平均叶面积指数LAI的步骤包括：

判断是否存在时间上缺失的被估算区域的叶面积指数；

根据时间序列重建方法对所述缺失的叶面积指数进行数据重建，得到全年完整的叶面积指数；

计算所述全年完整的叶面积指数的平均值作为年平均叶面积指数LAI。

可选地，获取单元201获取积雪日数SnowD的步骤包括：

判断所获取的积雪遥感产品是否存在时间上的缺失；

根据时间序列重建方法对所述缺失的积雪遥感产品进行数据重建，得到全年完整的积雪数据；

根据全年完整的积雪数据得到所述积雪日数SnowD。

可选地，获取单元201获取相关土壤数据包括：

根据活动层的深度进行分区，获取各种土壤数据的多层土壤数据；

计算各种参数的多层土壤数据的平均值，得到各个土壤数据平均后的多层平均土壤数据；

确定所述多层平均土壤数据为所述相关土壤数据。

本实施例中的遥感估算装置与实施例1所述的估算方法是基于同一构思下的发明，通过前述对估算方法及其各种变化形式的详细描述，本领域技术人员可以清楚的了解本实施例中遥感估算装置的实施过程，所以为了说明书的简洁，在此不再赘述。

实施例3

本发明实施例3提供一种可读存储介质，所述可读存储介质中存储有计算机程序，当所述计算机程序在计算机上运行时，使得所述计算机执行实施例1中所述的基于遥感的估算方法。本发明对此不再赘述。

在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，也可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，附图中的流程图和框图显示了根据本发明的多个实施例的装置、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

另外，在本发明各个实施例中的各功能模块可以集成在一起形成一个独立的部分，也可以是各个模块单独存在，也可以两个或两个以上模块集成形成一个独立的部分。

所述功能如果以软件功能模块的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。

以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。