CN113624687A - 一种冻土层厚度反演方法、装置及计算机设备 - Google Patents

一种冻土层厚度反演方法、装置及计算机设备 Download PDF

Info

Publication number
CN113624687A
CN113624687A CN202110929930.5A CN202110929930A CN113624687A CN 113624687 A CN113624687 A CN 113624687A CN 202110929930 A CN202110929930 A CN 202110929930A CN 113624687 A CN113624687 A CN 113624687A
Authority
CN
China
Prior art keywords
surface temperature
earth surface
fitting
inversion
average
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.)
Pending
Application number
CN202110929930.5A
Other languages
English (en)
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.)
Meihang Remote Sensing Information Co ltd
Aerial Photogrammetry and Remote Sensing Co Ltd
Original Assignee
Meihang Remote Sensing Information Co ltd
Aerial Photogrammetry and Remote Sensing Co Ltd
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 Meihang Remote Sensing Information Co ltd, Aerial Photogrammetry and Remote Sensing Co Ltd filed Critical Meihang Remote Sensing Information Co ltd
Priority to CN202110929930.5A priority Critical patent/CN113624687A/zh
Publication of CN113624687A publication Critical patent/CN113624687A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
    • G01J5/007—Radiation pyrometry, e.g. infrared or optical thermometry for earth observation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/171—Systems in which incident light is modified in accordance with the properties of the material investigated with calorimetric detection, e.g. with thermal lens detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/31—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/31—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
    • G01N21/35—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light
    • G01N21/3563—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light for analysing solids; Preparation of samples therefor
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
    • G01J2005/0077—Imaging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N2021/178—Methods for obtaining spatial resolution of the property being measured
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N2021/1793—Remote sensing

Landscapes

  • Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Biochemistry (AREA)
  • Analytical Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Health & Medical Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geology (AREA)
  • Radiation Pyrometers (AREA)

Abstract

本申请提供一种冻土层厚度反演方法、装置及计算机设备。该方法包括：首先基于预设区域对应的ETM数据计算预设区域内各第一测量点的反演地表温度值。然后从全部第一测量点中，选取反演地表温度值与第一测量点的实测温度值处于预设偏差范围内的第二测量点。再将部分第二测量点的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。把全部反演地表温度值代入拟合方程式，得到各第一测量点对应的全部第二平均地表温度值。基于第二平均地表温度值反演各第一测量点对应的冻土层厚度。本申请通过精度高且时相跨度小的ETM数据计算各第一测量点的反演地表温度值，进一步计算对应的冻土层厚度，使得反演得到的冻土层厚度结果精度高。

Description

一种冻土层厚度反演方法、装置及计算机设备
技术领域
本申请涉及遥感数据处理领域，尤其涉及一种冻土层厚度反演方法、装置及计算机设备。
背景技术
随着社会的发展，对于天然气水合物形成和保存的研究也在不断地深入。陆域天然气水合物形成于低温高压的冻土分布带中，因此，冻土层厚度是研究天然气水合物的一个非常重要的参量。同时，冻土层分布厚度与陆地表面温度之间存在着内在的联系。因此，利用遥感技术进行温度信息的反演是分析天然气水合物形成条件和保存条件中最关键的一步。
目前，用于进行陆地表面温度信息提取的遥感数据类型主要包括MODIS数据和ASTER数据。MODIS数据虽然易于获取，但其空间分辨率主要适用于大比例尺的宏观工作，难以满足高精度的研究场景；而ASTER数据的时相跨度太大，会间接导致计算结果不准确，难于满足地表温度信息的提取。
所以，现有技术存在冻土层厚度反演结果精度低的技术问题。
发明内容
为了解决上述技术问题，本发明提供了一种冻土层厚度反演方法、装置及计算机设备，具体方案如下：
第一方面，本申请实施例提供了一种冻土层厚度反演方法，所述方法包括：
基于预设区域对应的ETM数据计算所述预设区域内各第一测量点的反演地表温度值；
从全部所述第一测量点中，选取所述反演地表温度值与所述第一测量点的实测温度值处于预设偏差范围内的第二测量点；
将部分所述第二测量点的所述反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式，其中，所述第一平均地表温度值用于表征预设时间段内所述第一测量点对应的平均地表温度值；
把全部所述反演地表温度值代入所述拟合方程式，得到各第一测量点对应的全部第二平均地表温度值；
基于所述第二平均地表温度值反演各所述第一测量点对应的冻土层厚度。
根据本申请公开的一种具体实施方式，基于预设区域对应的ETM数据计算所述预设区域内各第一测量点的反演地表温度值的步骤，包括：
基于预设区域对应的ETM数据计算辐射亮度值；
根据公式TS＝K2/ln(K1/B(TS)+1)计算出预设区域内各第一测量点的反演地表温度值，其中，TS为反演地表温度值，B(TS)为温度为T的黑体在热红外波段的辐射亮度，K1＝666.09W/(m2·sr·μm)，K2＝1282.71K。
根据本申请公开的一种具体实施方式，所述ETM数据包括ETM影像，基于预设区域对应的ETM数据计算辐射亮度值的步骤，包括：
从ETM影像中获取所述预设区域内各第一测量点对应的目标像元；
基于公式FV＝(NDVI-NDVIS)/(NDVIV-NDVIS)计算全部所述目标像元对应的植被覆盖度，其中，NDVI为归一化差异植被指数，NDVIV＝0.43和NDVIS＝0，若所述目标像元的NDVI大于0.43，FV取值为1，若所述目标像元的NDVI小于0，FV取值为0；
基于各所述目标像元对应的所述植被覆盖度确定各所述目标像元的类型，所述目标像元的类型包括水体像元、建筑物像元和自然表面像元；
根据所述目标像元对应的地表辐射率计算公式计算所述目标像元的地表辐射率，其中，自然表面像元的地表辐射率计算公式为
Figure BDA0003210207040000031
建筑物像元的地表辐射率计算公式为
Figure BDA0003210207040000032
水体像元的地表辐射率为定值0.995；
根据公式B(TS)＝[Lλ-L↑-τ(1-ε)L↓]/(τ·ε)计算得到辐射亮度值，其中，B(TS)为所述辐射亮度，τ为0.6，ε为所述目标像元对应的地表辐射率，L↑为大气向上辐射亮度＝3.39W/(m2·sr·μm)，L↓为大气向下辐射亮度＝5.12W/(m2·sr·μm)。
根据本申请公开的一种具体实施方式，将部分所述第二测量点的所述反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式的步骤，包括：
按照预设选取规则选取部分所述第二测量点为拟合点；
根据公式tε(i)＝65.71321-1.00181N(i)-0.07273E(i)-0.00553H(i)计算各所述拟合点对应的第一平均地表温度值，其中，tε(i)为第i个拟合点对应的第一平均地表温度值，N(i)为第i个拟合点对应的纬度，E(i)为第i个拟合点对应的经度，H(i)为第i个拟合点对应的海拔高度；
将各拟合点的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。
根据本申请公开的一种具体实施方式，基于所述第二平均地表温度值反演各所述第一测量点对应的冻土层厚度的步骤，包括：
获取所述预设区域对应的地温梯度，并根据公式q＝gλ计算得到所述预设区域对应的地中热流，其中，q为所述预设区域对应的地中热流，g为所述预设区域对应的地温梯度，λ为定值；
基于公式
Figure BDA0003210207040000041
计算各所述第一测量点对应的冻土层厚度，其中，Hf(i)为第i个第一测量点对应的冻土层厚度，
Figure BDA0003210207040000042
为与第一平均地表温度tε(i)对应的第二平均地表温度，h为定值。
第二方面，本申请实施例提供了一种冻土层厚度反演装置，所述装置包括：
计算模块，用于基于预设区域对应的ETM数据计算所述预设区域内各第一测量点的反演地表温度值；
选取模块，用于从全部所述第一测量点中，选取所述反演地表温度值与所述第一测量点的实测温度值处于预设偏差范围内的第二测量点；
拟合模块，用于将部分所述第二测量点的所述反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式，其中，所述第一平均地表温度值用于表征预设时间段内所述第一测量点对应的平均地表温度值；
代入模块，用于把全部所述反演地表温度值代入所述拟合方程式，得到各第一测量点对应的全部第二平均地表温度值；
反演模块，用于基于所述第二平均地表温度值反演各所述第一测量点对应的冻土层厚度。
根据本申请公开的一种具体实施方式，所述选取模块具体用于按照预设选取规则选取部分所述第二测量点为拟合点；
根据公式tε(i)＝65.71321-1.00181N(i)-0.07273E(i)-0.00553H(i)计算各所述拟合点对应的第一平均地表温度值，其中，tε(i)为第i个拟合点对应的第一平均地表温度值，N(i)为第i个拟合点对应的纬度，E(i)为第i个拟合点对应的经度，H(i)为第i个拟合点对应的海拔高度；
将各拟合点的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。
根据本申请公开的一种具体实施方式，所述反演模块具体用于获取所述预设区域对应的地温梯度，并根据公式q＝gλ计算得到所述预设区域对应的地中热流，其中，q为所述预设区域对应的地中热流，g所述预设区域对应的地温梯度，λ为定值；
基于公式
Figure BDA0003210207040000051
计算各所述第一测量点对应的冻土层厚度，其中，Hf(i)为第i个第一测量点对应的冻土层厚度，
Figure BDA0003210207040000052
为与第一平均地表温度tε(i)对应的第二平均地表温度，h为定值。
第三方面，本申请实施例提供了一种计算机设备，所述计算机设备包括处理器和存储器，所述存储器存储有计算机程序，所述计算机程序在所述处理器上执行时实现第一方面中任一项实施例所述的冻土层厚度反演方法。
第四方面，本申请实施例提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序在处理器上执行时实现第一方面中任一项实施例所述的冻土层厚度反演方法。
相对于现有技术而言，本申请具有以下有益效果：
本申请基于预设区域对应的ETM数据计算预设区域内各第一测量点的反演地表温度值。然后从全部第一测量点中，选取反演地表温度值与第一测量点的实测温度值处于预设偏差范围内的第二测量点。再将部分第二测量点的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。把全部反演地表温度值代入拟合方程式，得到各第一测量点对应的全部第二平均地表温度值。基于第二平均地表温度值反演各第一测量点对应的冻土层厚度。本申请通过精度高且时相跨度小的ETM数据计算各第一测量点的反演地表温度值，进一步计算对应的冻土层厚度，使得反演得到的冻土层厚度结果精度高。
附图说明
为了更清楚地说明本发明的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对本发明保护范围的限定。在各个附图中，类似的构成部分采用类似的编号。
图1为本申请实施例提供的一种冻土层厚度反演方法的流程示意图；
图2为本申请实施例提供的一种冻土层厚度反演装置的模块框图。
具体实施方式
下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。
通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
在下文中，可在本发明的各种实施例中使用的术语“包括”、“具有”及其同源词仅意在表示特定特征、数字、步骤、操作、元件、组件或前述项的组合，并且不应被理解为首先排除一个或更多个其它特征、数字、步骤、操作、元件、组件或前述项的组合的存在或增加一个或更多个特征、数字、步骤、操作、元件、组件或前述项的组合的可能性。
此外，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
除非另有限定，否则在这里使用的所有术语(包括技术术语和科学术语)具有与本发明的各种实施例所属领域普通技术人员通常理解的含义相同的含义。所述术语(诸如在一般使用的词典中限定的术语)将被解释为具有与在相关技术领域中的语境含义相同的含义并且将不被解释为具有理想化的含义或过于正式的含义，除非在本发明的各种实施例中被清楚地限定。
下面结合附图，对本申请的一些实施方式作详细说明。在不冲突的情况下，下述的实施例及实施例中的特征可以相互结合。
参见图1，为本申请实施例提供的一种冻土层厚度反演方法的流程示意图，如图1所示，所述方法主要包括。
步骤S101，基于预设区域对应的ETM数据计算所述预设区域内各第一测量点的反演地表温度值。
目前，用于进行地表温度值反演的遥感数据类型主要包括中分辨率成像光谱仪(Moderate-resolution Imaging Spectroradiometer，简称MODIS)数据、先进星载热发射和反射辐射仪(Advanced Spaceborne Thermal Emission and Reflection Radiometer，简称ASTER)数据以及增强型专题绘图仪(Enhanced Thematic Mapper，简称ETM)数据。其中，MODIS数据虽然具有每日两次过境，易于获取的优势，但其1000m的空间分辨率主要适用于大比例尺的宏观工作，不适用于小面积的预设区域，其较低的空间分辨率难以满足工作精度。ASTER数据，包含14个波段，第十至第十四波段为热红外波段，分辨率为90m，是用于地表温度信息反演的主要波段。但是，ASTER数据的时相对应的时间跨度太大，若使用ASTER数据进行地表温度值的反演，会导致反演结果精度较低，无法用于后续的冻土厚度的反演工作。ETM数据包括1个全色波段空间分辨率为15m及7个多光谱波段，其中1-5、7波段空间分辨率为30m，第6波段为红外波段，空间分辨率60m，时相上的跨度比较小，同时数据质量好，有利于对预设区域的地表温度值的反演工作。
第一测量点是在预设区域内预先选取的用于反演对应的地表温度值，并实际测量对应地表温度值的地理位置点。
首先，获取预设区域对应的ETM影像，可以将ETM影像划分为面积大小和形状都相同的若干个像元。预设区域内的各第一测量点必然处于某个像元内，可以将各第一测量点所属的像元确定为各第一测量点对应的目标像元。具体实施时，由于第一测量点选取规则的不同，不同的第一测量点可能同时处于同一个目标像元。然后计算各目标像元的预设参数，得到各目标像元对应的第一测量点的反演温度值，其中，预设参数包括植被覆盖度、地表辐射率和辐射亮度值，具体计算步骤如下：
A.首先，采用混合像元分解法计算预设区域对应的ETM影像中各目标像元的植被覆盖度。具体地，基于公式FV＝(NDVI-NDVIS)/(NDVIV-NDVIS)计算全部目标像元对应的植被覆盖度，其中，NDVI为归一化差异植被指数，NDVIV＝0.43和NDVIS＝0，若目标像元的NDVI大于0.43，FV取值为1，若目标像元的NDVI小于0，FV取值为0。然后基于各目标像元对应的植被覆盖度以及预设的植被覆盖度与各目标像元类型的对应关系，确定各目标像元的类型，以此将ETM影像的地类大致分为水体、植被和裸地。其中，目标像元的类型包括水体像元、建筑物像元和自然表面像元，水体像元对应水体，建筑物像元对应裸地，自然表面像元对应植被。
B.根据不同目标像元的类型对应的地表辐射率计算公式计算目标像元的地表辐射率。其中，自然表面像元的地表辐射率计算公式为
Figure BDA0003210207040000101
建筑物像元的地表辐射率计算公式为
Figure BDA0003210207040000102
水体像元的地表辐射率可以设置为定值，优选地，水体像元的地表辐射率可以为0.995。
C.根据公式B(TS)＝[Lλ-L↑-τ(1-ε)L↓]/(τ·ε)计算得到辐射亮度值，其中，B(TS)是通过普朗克定律推得的黑体在地表温度值为TS时的辐射亮度，τ为大气在热红外波段的透过率，优选地，τ可以设置为定值0.6，ε为不同目标像元类型对应的地表辐射率，L↑为大气向上辐射亮度＝3.39W/(m2·sr·μm)，L↓为大气向下辐射亮度＝5.12W/(m2·sr·μm)。
D.在获取温度为TS的黑体在热红外波段的辐射亮度后，根据普朗克公式的反函数以及公式TS＝K2/ln(K1/B(TS)+1)计算出预设区域内各第一测量点的反演地表温度值，其中，TS为反演地表温度值，K1＝666.09W/(m2·sr·μm)，K2＝1282.71K，K为热力学温度单位开尔文。
上述步骤通过预设区域对应的ETM影像计算各预设参数，得到各第一测量点对应的反演地表温度值。选取ETM这类遥感数据，时相跨度小且数据质量好，有利于提高地表温度值的反演。
步骤S102，从全部所述第一测量点中，选取所述反演地表温度值与所述第一测量点的实测温度值处于预设偏差范围内的第二测量点。
具体实施时，各第一测量点的反演地表温度值可以通过数据统计表和灰度图的形式进行显示。若显示形式为灰度图，可以进一步选取目标像元并读取该目标像元点的反演地表温度值。在此基础上，可以按照不同的温度范围级别将灰度图和各目标像元点的反演地表温度值进行阈值分割，得到预设区域的反演地表温度分布图，以使基于反演地表温度图进行分析使用。例如，可以按照<15℃、15-20℃、20-22℃、22-25℃、25-30℃、30-35℃、35-40℃七个级别进行阈值分割，这里不作进一步限定。
然后，可以根据不同的精度需求或实际需求，选择由不同时期的ETM数据计算得到的反演地表温度值与实测温度值进行比较。例如，对各第一测量的地表温度值的实际测量时间是2021年8月，若需要较高的反演精度，可以选择2021年8月对应时相的ETM数据得到各反演地表温度值；也可以选择历史年份中8月对应的ETM数据得到各反演地表温度值。具体地，从全部第一测量点中，选取反演地表温度值与第一测量点的实测温度值处于预设偏差范围内的第二测量点。由于实际影响地表温度值的因素较多，例如水、风、积雪覆盖、降水、云层、空气湿度、雾及悬浮物质等，以及反演时采用的ETM数据并非瞬时数据，无法与现时温度情况完全吻合。因此，通过设置“预设偏差范围”从全部第一测量点中选取与对应的实测温度值大致相同的第二测量点即可。
步骤S103，将部分所述第二测量点的所述反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式，其中，所述第一平均地表温度值用于表征预设时间段内所述第一测量点对应的平均地表温度值。
由于ETM数据反演的是瞬间的地表温度值，如某天或者某月，而计算冻土层厚度时用的是年平均地表温度值，因此，需要利用计算得到的反演地表温度值和部分实际的年平均地表温度值进行拟合，进一步求算整个预设区域的年平均地表温度值。
到达陆地表面的太阳辐射的纬向变化特征是冻土分带性的原因，因而也决定了年平均地表温度值在空间上的分带性规律。首先，按照预设选取规则选取部分第二测量点为拟合点。可以按照用户的输入进行选取，也可以是根据预设选取规则选择均匀分布于预设区域的第二测量点来进行拟合。
具体实施时，可以根据下述公式计算各拟合点对应的第一平均地表温度值:tε(i)＝65.71321-1.00181N(i)-0.07273E(i)-0.00553H(i),其中，tε(i)为第i个拟合点对应的第一平均地表温度值，N(i)为第i个拟合点对应的纬度，E(i)为第i个拟合点对应的经度，H(i)为第i个拟合点对应的海拔高度。上式表明，tε(i)与纬度N(i)、经度E(i)和海拔H(i)均呈负相关，即tε(i)随纬度、经度和海拔增大而降低。具体地，可以通过对应预设区域的数字高程模型(Digital Elevation Model,简称DEM)数据进行程序编写，完成各第一测量点的纬度、经度和海拔高度信息提取工作。然后将各拟合点的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。
步骤S104，把全部所述反演地表温度值代入所述拟合方程式，得到各第一测量点对应的全部第二平均地表温度值。
拟合方程式表征各拟合点的反演地表温度值与对应的第一平均地表温度值间具有较高的相关系数，说明通过ETM数据计算得到的反演地表温度值与对应的第一平均地表温度值间存在着紧密的联系，从而为进一步反演冻土层厚度奠定基础。
步骤S105，基于所述第二平均地表温度值反演各所述第一测量点对应的冻土层厚度。
在求得预设区域中各第一测量点对应的第二平均地表温度值的基础上，可以进一步求取各第一测量点对应的冻土层厚度。首先获取预设区域对应的地温梯度。在地形变化较大的区域，地温梯度变化较大，因此要通过实测资料的研究取值或者取经验值。然后根据公式q＝gλ计算得到预设区域对应的地中热流。其中，q为预设区域对应的地中热流，g为预设区域对应的地温梯度，λ为冻土层的导热系数，对于匀质地层，λ可以取定值。再基于公式
Figure BDA0003210207040000131
计算各第一测量点对应的冻土层厚度，其中，Hf(i)为第i个第一测量点对应的冻土层厚度，
Figure BDA0003210207040000132
为与第一平均地表温度tε(i)对应的第二平均地表温度，h为地温年变化深度，通常可以取经验值或者定值。
本申请提供的冻土层厚度反演方法，通过时相跨度小且空间分辨率高的ETM数据计算得到反演地表温度值，然后将反演地表温度值与实测温度进行比较，将相似性高的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。通过拟合方程式进一步得到全部第一测量点对应的第二平均地表温度值，再根据第二平均地表温度值反演得到对应的冻土层厚度。本申请通过上述技术手段可以反演得到高精度的冻土层厚度结果，提高了冻土层厚度反演的准确性。
与上述方法实施例相对应，参见图2，本发明还提供一种冻土层厚度反演装置200，所述冻土层厚度反演装置200包括：
计算模块201，用于基于预设区域对应的ETM数据计算所述预设区域内各第一测量点的反演地表温度值；
选取模块202，用于从全部所述第一测量点中，选取所述反演地表温度值与所述第一测量点的实测温度值处于预设偏差范围内的第二测量点；
拟合模块203，用于将部分所述第二测量点的所述反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式，其中，所述第一平均地表温度值用于表征预设时间段内所述第一测量点对应的平均地表温度值；
代入模块204，用于把全部所述反演地表温度值代入所述拟合方程式，得到各第一测量点对应的全部第二平均地表温度值；
反演模块205，用于基于所述第二平均地表温度值反演各所述第一测量点对应的冻土层厚度。
具体实施时，所述选取模块具体用于按照预设选取规则选取部分所述第二测量点为拟合点；
根据公式tε(i)＝65.71321-1.00181N(i)-0.07273E(i)-0.00553H(i)计算各所述拟合点对应的第一平均地表温度值，其中，tε(i)为第i个拟合点对应的第一平均地表温度值，N(i)为第i个拟合点对应的纬度，E(i)为第i个拟合点对应的经度，H(i)为第i个拟合点对应的海拔高度；
将各拟合点的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。
具体实施时，所述反演模块具体用于获取所述预设区域对应的地温梯度，并根据公式q＝gλ计算得到所述预设区域对应的地中热流，其中，q为所述预设区域对应的地中热流，g所述预设区域对应的地温梯度，λ为定值；
基于公式
Figure BDA0003210207040000151
计算各所述第一测量点对应的冻土层厚度，其中，Hf(i)为第i个第一测量点对应的冻土层厚度，
Figure BDA0003210207040000152
为与第一平均地表温度tε(i)对应的第二平均地表温度，h为定值。
此外，还提供一种计算机设备，计算机设备包括处理器和存储器，存储器存储有计算机程序，计算机程序在处理器上执行时实现上述冻土层厚度反演方法。
此外，还提供一种计算机可读存储介质，计算机可读存储介质存储有计算机程序，计算机程序在处理器上执行时实现上述冻土层厚度反演方法。
本申请提供的冻土层厚度反演装置、计算机设备和计算机可读存储介质，通过时相跨度小且空间分辨率高的ETM数据计算得到反演地表温度值，然后将反演地表温度值与实测温度进行比较，将相似性高的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。通过拟合方程式进一步得到全部第一测量点对应的第二平均地表温度值，再根据第二平均地表温度值反演得到对应的冻土层厚度。通过上述技术手段可以反演得到高精度的冻土层厚度结果，提高了冻土层厚度反演的准确性。
所提供的计算机设备和计算机可读存储介质的具体实施过程，可以参见上述实施例提供的冻土层厚度反演方法的具体实施过程，在此不再一一赘述。
在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，也可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，附图中的流程图和结构图显示了根据本发明的多个实施例的装置、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，结构图和/或流程图中的每个方框、以及结构图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
另外，在本发明各个实施例中的各功能模块或单元可以集成在一起形成一个独立的部分，也可以是各个模块单独存在，也可以两个或更多个模块集成形成一个独立的部分。
所述功能如果以软件功能模块的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是智能手机、个人计算机、服务器、或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种冻土层厚度反演方法，其特征在于，所述方法包括：
基于预设区域对应的ETM数据计算所述预设区域内各第一测量点的反演地表温度值；
从全部所述第一测量点中，选取所述反演地表温度值与所述第一测量点的实测温度值处于预设偏差范围内的第二测量点；
将部分所述第二测量点的所述反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式，其中，所述第一平均地表温度值用于表征预设时间段内所述第一测量点对应的平均地表温度值；
把全部所述反演地表温度值代入所述拟合方程式，得到各第一测量点对应的全部第二平均地表温度值；
基于所述第二平均地表温度值反演各所述第一测量点对应的冻土层厚度。
2.根据权利要求1所述的方法，其特征在于，基于预设区域对应的ETM数据计算所述预设区域内各第一测量点的反演地表温度值的步骤，包括：
基于预设区域对应的ETM数据计算辐射亮度值；
根据公式TS＝K2/ln(K1/B(TS)+1)计算出预设区域内各第一测量点的反演地表温度值，其中，TS为反演地表温度值，B(TS)为温度为T的黑体在热红外波段的辐射亮度，K1＝666.09W/(m2·sr·μm)，K2＝1282.71K，K为热力学温度单位开尔文。
3.根据权利要求2所述的方法，其特征在于，所述ETM数据包括ETM影像，基于预设区域对应的ETM数据计算辐射亮度值的步骤，包括：
从ETM影像中获取所述预设区域内各第一测量点对应的目标像元；
基于公式FV＝(NDVI-NDVIS)/(NDVIV-NDVIS)计算全部所述目标像元对应的植被覆盖度，其中，NDVI为归一化差异植被指数，NDVIV＝0.43和NDVIS＝0，若所述目标像元的NDVI大于0.43，FV取值为1，若所述目标像元的NDVI小于0，FV取值为0；
基于各所述目标像元对应的所述植被覆盖度确定各所述目标像元的类型，所述目标像元的类型包括水体像元、建筑物像元和自然表面像元；
根据所述目标像元对应的地表辐射率计算公式计算所述目标像元的地表辐射率，其中，自然表面像元的地表辐射率计算公式为
Figure FDA0003210207030000021
建筑物像元的地表辐射率计算公式为
Figure FDA0003210207030000022
水体像元的地表辐射率为定值0.995；
根据公式B(TS)＝[Lλ-L↑-τ(1-ε)L↓]/(τ·ε)计算得到辐射亮度值，其中，B(TS)为所述辐射亮度，τ为0.6，ε为所述目标像元对应的地表辐射率，L↑为大气向上辐射亮度＝3.39W/(m2·sr·μm)，L↓为大气向下辐射亮度＝5.12W/(m2·sr·μm)。
4.根据权利要求1所述的方法，其特征在于，将部分所述第二测量点的所述反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式的步骤，包括：
按照预设选取规则选取部分所述第二测量点为拟合点；
根据公式tε(i)＝65.71321-1.00181N(i)-0.07273E(i)-0.00553H(i)计算各所述拟合点对应的第一平均地表温度值，其中，tε(i)为第i个拟合点对应的第一平均地表温度值，N(i)为第i个拟合点对应的纬度，E(i)为第i个拟合点对应的经度，H(i)为第i个拟合点对应的海拔高度；
将各拟合点的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。
5.根据权利要求1所述的方法，其特征在于，基于所述第二平均地表温度值反演各所述第一测量点对应的冻土层厚度的步骤，包括：
获取所述预设区域对应的地温梯度，并根据公式q＝gλ计算得到所述预设区域对应的地中热流，其中，q为所述预设区域对应的地中热流，g为所述预设区域对应的地温梯度，λ为定值；
基于公式
Figure FDA0003210207030000031
计算各所述第一测量点对应的冻土层厚度，其中，Hf(i)为第i个第一测量点对应的冻土层厚度，
Figure FDA0003210207030000032
为与第一平均地表温度tε(i)对应的第二平均地表温度，h为定值。
6.一种冻土层厚度反演装置，其特征在于，所述装置包括：
计算模块，用于基于预设区域对应的ETM数据计算所述预设区域内各第一测量点的反演地表温度值；
选取模块，用于从全部所述第一测量点中，选取所述反演地表温度值与所述第一测量点的实测温度值处于预设偏差范围内的第二测量点；
拟合模块，用于将部分所述第二测量点的所述反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式，其中，所述第一平均地表温度值用于表征预设时间段内所述第一测量点对应的平均地表温度值；
代入模块，用于把全部所述反演地表温度值代入所述拟合方程式，得到各第一测量点对应的全部第二平均地表温度值；
反演模块，用于基于所述第二平均地表温度值反演各所述第一测量点对应的冻土层厚度。
7.根据权利要求6所述的装置，其特征在于，所述选取模块具体用于按照预设选取规则选取部分所述第二测量点为拟合点；
根据公式tε(i)＝65.71321-1.00181N(i)-0.07273E(i)-0.00553H(i)计算各所述拟合点对应的第一平均地表温度值，其中，tε(i)为第i个拟合点对应的第一平均地表温度值，N(i)为第i个拟合点对应的纬度，E(i)为第i个拟合点对应的经度，H(i)为第i个拟合点对应的海拔高度；
将各拟合点的反演地表温度值与对应的第一平均地表温度值进行拟合，得到拟合方程式。
8.根据权利要求6所述的装置，其特征在于，所述反演模块具体用于获取所述预设区域对应的地温梯度，并根据公式q＝gλ计算得到所述预设区域对应的地中热流，其中，q为所述预设区域对应的地中热流，g所述预设区域对应的地温梯度，λ为定值；
基于公式
Figure FDA0003210207030000041
计算各所述第一测量点对应的冻土层厚度，其中，Hf(i)为第i个第一测量点对应的冻土层厚度，
Figure FDA0003210207030000051
为与第一平均地表温度tε(i)对应的第二平均地表温度，h为定值。
9.一种计算机设备，其特征在于，所述计算机设备包括处理器和存储器，所述存储器存储有计算机程序，所述计算机程序在所述处理器上执行时实现权利要求1至5中任一项所述的冻土层厚度反演方法。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序在处理器上执行时实现权利要求1至5中任一项所述的冻土层厚度反演方法。
CN202110929930.5A 2021-08-13 2021-08-13 一种冻土层厚度反演方法、装置及计算机设备 Pending CN113624687A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110929930.5A CN113624687A (zh) 2021-08-13 2021-08-13 一种冻土层厚度反演方法、装置及计算机设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110929930.5A CN113624687A (zh) 2021-08-13 2021-08-13 一种冻土层厚度反演方法、装置及计算机设备

Publications (1)

Publication Number Publication Date
CN113624687A true CN113624687A (zh) 2021-11-09

Family

ID=78385282

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110929930.5A Pending CN113624687A (zh) 2021-08-13 2021-08-13 一种冻土层厚度反演方法、装置及计算机设备

Country Status (1)

Country Link
CN (1) CN113624687A (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107290503A (zh) * 2017-06-20 2017-10-24 武汉大学 冻土中水非饱和流运动和溶质迁移通量及等效参数监测的方法
CN111337550A (zh) * 2020-03-13 2020-06-26 山东航向电子科技有限公司 基于gnss干涉信号的冻土土壤湿度监测方法
CN111812645A (zh) * 2020-06-10 2020-10-23 西南交通大学 冻土季节形变的卫星干涉测量方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107290503A (zh) * 2017-06-20 2017-10-24 武汉大学 冻土中水非饱和流运动和溶质迁移通量及等效参数监测的方法
CN111337550A (zh) * 2020-03-13 2020-06-26 山东航向电子科技有限公司 基于gnss干涉信号的冻土土壤湿度监测方法
CN111812645A (zh) * 2020-06-10 2020-10-23 西南交通大学 冻土季节形变的卫星干涉测量方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
周强: "陆域天然气水合物遥感探测研究", 《中国博士学位论文全文数据库 基础科学辑》, no. 8, 15 August 2006 (2006-08-15), pages 58 - 68 *
王春娇: "基于D-InSAR的东北多年冻土退化区地表形变研究", 《中国博士学位论文全文数据库 基础科学辑》, no. 1, 15 January 2016 (2016-01-15), pages 29 - 36 *

Similar Documents

Publication Publication Date Title
Zhang et al. Lakes’ state and abundance across the Tibetan Plateau
Azizian et al. Determining the best remotely sensed DEM for flood inundation mapping in data sparse regions
Feyisa et al. Locally optimized separability enhancement indices for urban land cover mapping: Exploring thermal environmental consequences of rapid urbanization in Addis Ababa, Ethiopia
Wang et al. Application of a normalized difference impervious index (NDII) to extract urban impervious surface features based on Landsat TM images
Tang et al. An improved NDVI-based threshold method for estimating land surface emissivity using MODIS satellite data
US20210209803A1 (en) Computer-based method and system for geo-spatial analysis
Fu et al. Estimating landscape net ecosystem exchange at high spatial–temporal resolution based on Landsat data, an improved upscaling model framework, and eddy covariance flux measurements
Zhang et al. Glacier change and glacier runoff variation in the Tuotuo River basin, the source region of Yangtze River in western China
Zhao et al. Estimating fractional vegetation cover from leaf area index and clumping index based on the gap probability theory
Wu et al. Evaluating and comparing performances of topographic correction methods based on multi-source DEMs and Landsat-8 OLI data
Zhang et al. A stepwise downscaling method for generating high-resolution land surface temperature from AMSR-E data
Maeda Downscaling MODIS LST in the East African mountains using elevation gradient and land-cover information
Liu et al. Comparative analysis of fractional vegetation cover estimation based on multi-sensor data in a semi-arid sandy area
Ramdani et al. Spatio-temporal analysis of urban temperature in Bandung City, Indonesia
Amani et al. A parametric method for estimation of leaf area index using landsat ETM+ data
Hassanpour et al. Modification on optical trapezoid model for accurate estimation of soil moisture content in a maize growing field
Li et al. Estimation of aboveground vegetation biomass based on Landsat-8 OLI satellite images in the Guanzhong Basin, China
Alqasemi et al. Retrieval of monthly maximum and minimum air temperature using MODIS aqua land surface temperature data over the United Arab Emirates
Luo et al. Spatial downscaling of MODIS land surface temperature based on a geographically and temporally weighted autoregressive model
Natsagdorj et al. An integrated methodology for soil moisture analysis using multispectral data in Mongolia
Yeom et al. Mapping rice area and yield in northeastern asia by incorporating a crop model with dense vegetation index profiles from a geostationary satellite
Guo et al. Using VIIRS-DNB and landsat data for impervious surface area mapping in an arid/semiarid region
Yang et al. Estimating evapotranspiration from terrestrial groundwater-dependent ecosystems using Landsat images
Zhang et al. Comparative study of remote sensing estimation methods for grassland fractional vegetation coverage–a grassland case study performed in Ili prefecture, Xinjiang, China
Jin-Ming et al. Temporal-spatial variations and influencing factors of Lakes in inland arid areas from 2000 to 2017: a case study in Xinjiang

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