CN111982294A - 集成热红外与再分析资料的全天候地表温度生成方法 - Google Patents

集成热红外与再分析资料的全天候地表温度生成方法 Download PDF

Info

Publication number
CN111982294A
CN111982294A CN202010707434.0A CN202010707434A CN111982294A CN 111982294 A CN111982294 A CN 111982294A CN 202010707434 A CN202010707434 A CN 202010707434A CN 111982294 A CN111982294 A CN 111982294A
Authority
CN
China
Prior art keywords
surface temperature
lst
data
pixel
weather
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
CN202010707434.0A
Other languages
English (en)
Other versions
CN111982294B (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.)
University of Electronic Science and Technology of China
Chengdu Univeristy of Technology
Original Assignee
University of Electronic Science and Technology of China
Chengdu Univeristy of Technology
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 University of Electronic Science and Technology of China, Chengdu Univeristy of Technology filed Critical University of Electronic Science and Technology of China
Priority to CN202010707434.0A priority Critical patent/CN111982294B/zh
Publication of CN111982294A publication Critical patent/CN111982294A/zh
Application granted granted Critical
Publication of CN111982294B publication Critical patent/CN111982294B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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/48—Thermography; Techniques using wholly visual means
    • G01J5/485—Temperature profile
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

本发明提供了一种集成热红外与再分析资料的全天候地表温度生成方法，该方法基于随机森林集成热红外遥感与再分析资料实现，具体包括以下步骤：数据预处理与时空匹配；时间算法的实现；空间算法的实现；以及生成最终的1km全天候地表温度。本发明引入再分析数据，以提供较长的地表温度及其相关观测数据，将其与热红外卫星遥感进行集成得到长时间序列、高时空分辨率的全天候地表温度，得到的地表温度同传统的热红外卫星遥感地表温度产品具有良好的一致性，且经过实测地表温度的检验，精度在可接受范围，图像质量较高，能展现更多地表温度的空间细节信息，总体质量和精度都要优于集成热红外卫星遥感数据和被动微波数据的方法。

Description

集成热红外与再分析资料的全天候地表温度生成方法
技术领域
本发明涉及卫星遥感地表温度产品生成领域，具体涉及一种基于随机森林集成热红外遥感与再分析资料的1km全天候地表温度生成方法。
背景技术
地表温度(LST)在地球表面和大气界面的物质循环和能量交换中一个非常重要的物理量。地表温度的准确估算对能量循环、气候预报、农业生产等方面具有推动作用。及时掌握区域和全球尺度上的地表温度时空分布，尤其是全面、完整和连续的地表温度时空分布信息，对地气系统能量平衡和生态系统研究具有重要意义。随着遥感技术的大力发展，越来越丰富的遥感数据被获取，使得精确估算全球范围的地表温度成为可能。过去几十年来，利用星载热红外遥感数据反演地表温度得到了极大发展，其算法较为成熟，空间分辨率较高，反演的精度基本能够满足实际工作与公开应用的要求。然而，在类似青藏高原这样具有频繁云覆盖特征的地区，传统的地表温度产品无法达到全天候监测的应用要求。使用集成热红外和被动微波数据的方法能够生成具有较高分辨率的全天候地表温度，并有效地填补数据缺失部分，但是由于该方法使用涉及的数据源时间跨度短，时间分辨率不够，无法满足较长时间序列、高时空分辨率地表温度数据生产需要。再分析数据，如全球土地数据同化系统(GLDAS)，同样可以提供时空连续的全天候地表温度，与传统的地表温度产品相比，再分析数据能够提供具有更高时间分辨率的地表温度及其相关参数，但受制于其空间分辨率较低，使其无法在具体的区域研究中提供具有足够空间变化细节的地表温度。
发明内容
针对上述问题，本发明的目的是提供一种将热红外遥感地表温度产品(TIR LST)和再分析数据集成起来的方法，并结合两种数据各自的优点，以产生具有高时空分辨率的长时间序列全天候地表温度。一方面，该方法生成的全天候1km地表温度与传统的TIR LST产品具有很好的一致性，有效地填补了因云层污染导致的数据空洞；另一方面，同集成热红外和被动微波遥感数据生成的1km分辨率的全天候地表温度产品相比，使用本发明方案生成的产品具有更好的精度和图像质量。另外，该方法在不同时间跨度下都表现出了良好的适应性，使得将其进一步应用，以生成长时间序列的全天候地表温度数据，从而对促进区域和全球范围相关研究具有重大意义。
本发明首先从众多机器学习算法中选择随机森林回归(random forestregression)用于构建TIR LST和再分析资料之间的映射(函数)关系，并基于时间和空间维度提出两种算法用以刻画和描述地表温度的特征。一种是时间算法(temporalalgorithm)，即描绘时间维度上地表温度变化的算法，该算法利用晴空条件下的低分辨率再分析数据生成全天候高分辨率的地表温度，最终生成的地表温度称为LSTT-RF。同时，考虑到地表温度在时空维度上变化的差异性，针对性地提出一种空间算法(spatialalgorithm)，用以描述地表温度与其相关的空间因子之间的映射关系，该算法利用全天候条件下高空间分辨率数据，结合较低空间分辨率的再分析资料，生成全天候高空间分辨率地表温度，通过该算法得到的地表温度称为LSTS-RF。最终，通过随机森林回归的方法，结合所有与地表温度相关的时空因子，集成时间算法和空间算法的结果，生成精度最好的1km全天候地表温度(LSTT-RF-S)。其中，
(1)时间算法的公式如下：
Figure BDA0002594835400000021
式中，LSTTIR-Q表示高空间分辨率的卫星地表温度产品中像元Q对应的时间序列，单位为K，本发明中默认时长为1年；GLsd代表包含像元Q的再分析数据相关因子的时间序列；RFQ表示LSTTIR-Q与GLsd之间建立的随机森林回归关系；PT-m-sd表示第m个地表温度相关因子的时间序列，本发明中设置m＝12。考虑云层覆盖的影响，将利用公式(1)求取初始1km全天候地表温度的过程具体到晴空和非晴空的条件下：
Figure BDA0002594835400000022
式中，
Figure BDA0002594835400000023
和
Figure BDA0002594835400000024
分别表示在晴空和非晴空条件下，像元所对应的1km地表温度估算值；td-cd-n表示非晴空条件下的第n个像元。
进一步，像元Q对应的1km地表温度估算值和热红外遥感地表温度产品之间的系统误差可以表示为：
Figure BDA0002594835400000031
式中，
Figure BDA0002594835400000032
和
Figure BDA0002594835400000033
分别代表像元Q对应的1km地表温度估算值和热红外遥感地表温度产品的平均值，ΔδT-Q表示像元Q对应的系统误差。
因此，通过分别优化晴空和非晴空条件下地表温度估计值，可以获得最终的全天候地表温度(LSTT-RF)，如下所示：
Figure BDA0002594835400000034
式中，LSTT-RF-sd-Q和LSTT-RF-cd-Q分别代表晴空和非晴空条件下像元Q对应的1km地表温度值，LSTT-RF-Q表示像元Q对应的最终的LSTT-RF值。
(2)在空间算法中，对于包含像元P的特定区域，在再分析数据的像元尺度下，地表温度可以与其相关的空间因子的关系可以表示为：
Figure BDA0002594835400000035
式中，其中LSTS-low表示从再分析资料中得到的低分辨率地表温度数据，单位为K；NS-k-low表示对应LSTS-low分辨率的第k个与地表温度相关的空间因子；RFs表示像元P对应的地表温度与多个与其相关的空间因子之间的基于随机森林回归的函数模型。
将公式(4)应用到较高空间分辨率情况下(例如1km)，则可通过公式(6)估算像元P的初始1km地表温度：
Figure BDA0002594835400000036
式中，
Figure BDA0002594835400000037
在表示较高分辨率情况下像元P对应的地表温度估计值，单位K；NS-k-high-P表示像元P对应的第k个较高空间分辨率的地表温度空间因子；RFs表示较高分辨率下像元P对应的地表温度与多个相关空间因子之间随机森林回归关系。
进一步，像元P对应的1km地表温度估算值和热红外遥感地表温度产品之间的系统误差可以表示为：
Figure BDA0002594835400000041
式中，
Figure BDA0002594835400000042
代表像元P对应的1km地表温度估算值的平均值，单位为K；
Figure BDA0002594835400000043
代表像元P对应的1km地表温度产品的平均值，单位为K。
因此，通过分别优化晴空和非晴空条件下地表温度估计值，可以获得最终的全天候地表温度(LSTS-RF)：
Figure BDA0002594835400000044
式中，LSTS-RF-sd-P和LSTS-RF-cd-P分别代表晴空和非晴空条件下像元P对应的1km地表温度值，LSTS-RF-P表示像元P对应的最终LSTS-RF值。
为了提高该方法在复杂条件下的适应能力，应当充分表征地表温度的时空变化特征。因此，本发明将时间和空间算法得到的结果进行融合，通过以下表达式生成最终的1km地表温度产品：
LSTT-RF-S-P＝RFM(LSTT-RF-P,LSTS-RF-P) (9)
其中RFM表示基于随机森林的回归函数，用于衡量像元P对应的从时间和空间算法生成的地表温度对最终地表温度的贡献度。
本发明所采用方法的优点在于：集成热红外卫星遥感数据和被动微波数据生产全天候地表温度的方法，受制于时间分辨率和时间跨度的影响，无法满足生产高时空分辨率的长时间序列全天候地表温度的要求。本发明引入具有更高时间分辨率的再分析数据，以提供较长的地表温度及其相关观测数据，将其与热红外卫星遥感进行集成得到长时间序列、高时空分辨率的全天候地表温度。本发明得到的地表温度同传统的热红外卫星遥感地表温度产品具有良好的一致性，且经过实测地表温度的检验，精度在可接受范围，图像质量较高，能展现更多地表温度的空间细节信息，总体质量和精度都要优于集成热红外卫星遥感数据和被动微波数据的方法。
附图说明
图1是本发明的全天候地表温度生成方法流程图
图2是包含研究区在内的地表温度图
图3是2014年MODIS地表温度和本发明提出方法得到的地表温度之间的空间分布评价指数图
图4是时间算法和空间算法分别在2014年1月1日，4月1日，7月1日和10月1日得到的地表温度图
图5是时间算法对最终结果的年内平均贡献度
图6是2010年本发明温度结果和TIR LST之间的散点图
图7是3种地表温度产品在2010年4月1日的温度图
图8-10是3种地表温度与地面实测站点之间的散点图
具体实施方式
下面结合附图和具体实施例进行说明。
本发明提出的一种集成热红外与再分析资料的全天候地表温度生成方法首先从众多机器学习算法中选择随机森林回归(random forest regression)用于构建TIR LST和再分析资料之间的映射(函数)关系，并基于时间和空间维度提出两种算法用以刻画和描述地表温度的特征。一种是时间算法(temporal algorithm)，即描绘时间维度上地表温度变化的算法，该算法利用晴空条件下的低分辨率再分析数据生成全天候高分辨率的地表温度，最终生成的地表温度称为LSTT-RF。同时，考虑到地表温度在时空维度上变化的差异性，针对性地提出一种空间算法(spatial algorithm)，用以描述地表温度与其相关的空间因子之间的映射关系，该算法利用全天候条件下高空间分辨率数据，结合较低空间分辨率的再分析资料，生成全天候高空间分辨率地表温度，通过该算法得到的地表温度称为LSTS-RF。最终，通过随机森林回归的方法，结合所有与地表温度相关的时空因子，集成时间算法和空间算法的结果，生成精度最好的1km全天候地表温度(LSTT-RF-S)。其中，
(1)时间算法的公式如下：
Figure BDA0002594835400000061
式中，LSTTIR-Q表示高空间分辨率的卫星地表温度产品中像元Q对应的时间序列，单位为K，本发明中默认时长为1年；GLsd代表包含像元Q的再分析数据相关因子的时间序列；RFQ表示LSTTIR-Q与GLsd之间建立的随机森林回归关系；PT-m-sd表示第m个地表温度相关因子的时间序列，本发明中设置m＝12。考虑云层覆盖的影响，将利用公式(1)求取初始1km全天候地表温度的过程具体到晴空和非晴空的条件下：
Figure BDA0002594835400000062
式中，
Figure BDA0002594835400000063
和
Figure BDA0002594835400000064
分别表示在晴空和非晴空条件下，像元所对应的1km地表温度估算值；td-cd-n表示非晴空条件下的第n个像元。
进一步，像元Q对应的1km地表温度估算值和热红外遥感地表温度产品之间的系统误差可以表示为：
Figure BDA0002594835400000065
式中，
Figure BDA0002594835400000066
和
Figure BDA0002594835400000067
分别代表像元Q对应的1km地表温度估算值和热红外遥感地表温度产品的平均值，ΔδT-Q表示像元Q对应的系统误差。
因此，通过分别优化晴空和非晴空条件下地表温度估计值，可以获得最终的全天候地表温度(LSTT-RF)，如下所示：
Figure BDA0002594835400000071
式中，LSTT-RF-sd-Q和LSTT-RF-cd-Q分别代表晴空和非晴空条件下像元Q对应的1km地表温度值，LSTT-RF-Q表示像元Q对应的最终的LSTT-RF值。
(2)在空间算法中，对于包含像元P的特定区域，在再分析数据的像元尺度下，地表温度可以与其相关的空间因子的关系可以表示为：
Figure BDA0002594835400000072
式中，其中LSTS-low表示从再分析资料中得到的低分辨率地表温度数据，单位为K；NS-k-low表示对应LSTS-low分辨率的第k个与地表温度相关的空间因子；RFs表示像元P对应的地表温度与多个与其相关的空间因子之间的基于随机森林回归的函数模型。
将公式(4)应用到较高空间分辨率情况下(例如1km)，则可通过公式(6)估算像元P的初始1km地表温度：
Figure BDA0002594835400000073
式中，
Figure BDA0002594835400000074
在表示较高分辨率情况下像元P对应的地表温度估计值，单位K；NS-k-high-P表示像元P对应的第k个较高空间分辨率的地表温度空间因子；RFs表示较高分辨率下像元P对应的地表温度与多个相关空间因子之间随机森林回归关系。
进一步，像元P对应的1km地表温度估算值和热红外遥感地表温度产品之间的系统误差可以表示为：
Figure BDA0002594835400000075
式中，
Figure BDA0002594835400000076
代表像元P对应的1km地表温度估算值的平均值，单位为K；
Figure BDA0002594835400000077
代表像元P对应的1km地表温度产品的平均值，单位为K。
因此，通过分别优化晴空和非晴空条件下地表温度估计值，可以获得最终的全天候地表温度(LSTS-RF)：
Figure BDA0002594835400000081
式中，LSTS-RF-sd-P和LSTS-RF-cd-P分别代表晴空和非晴空条件下像元P对应的1km地表温度值，LSTS-RF-P表示像元P对应的最终LSTS-RF值。
为了提高该方法在复杂条件下的适应能力，应当充分表征地表温度的时空变化特征。因此，本发明将时间和空间算法得到的结果进行融合，通过以下表达式生成最终的1km地表温度产品：
LSTT-RF-S-P＝RFM(LSTT-RF-P,LSTS-RF-P) (9)
其中RFM表示基于随机森林的回归函数，用于衡量像元P对应的从时间和空间算法生成的地表温度对最终地表温度的贡献度。
本实施例选择区域为青藏高原和其周边地区(19.5°N-44°N,73°E-104°E)，青藏高原是世界上海拔最高的高原，平均海拔超过4000米，俗称"世界屋脊"。由于其独特的地理位置和地形，这一地区具有非常复杂的气候条件和各种土地覆盖类型。青藏高原的动力、热力效应对亚洲季风的形成具有关键性的作用，同时也是全球气候变化的敏感区域之一。地表温度是地表能量平衡中的主要参数，对于进一步研究该地区的陆地-大气相互作用具有重要意义，因此迫切需要一种长时间序列、高时空分辨率的地表温度数据。
本实施例中具体采用的数据集包括：1)遥感数据集：MODIS逐日地表温度产品(MOD21A1/MYD21A1,1km)，16天合成归一化植被指数(NDVI)产品(MOD12A2,1km)，逐日地表反照率产品(MCD43A3，1km)，逐日归一化雪指数(NDSI)产品(MOD10A1，500m)，DEM(250m)；2)再分析资料：全球陆面数据同化系统(GLDAS)，时间分辨率3小时，空间分辨率为0.25°(GLDAS_NOAH025_3H)；3)验证数据集：9个科研站点的地表实测温度数据。
1.数据预处理与时空匹配
首先，利用IDL语言和MRT工具对MODIS数据产品进行批量处理(格式转换、投影变换、拼接、裁剪)；对于MOD21A1/MYD21A1数据，利用相应的质量控制数据，提取高质量的MODIS LST像元。
其次，对GLDAS数据的参数时间进行三次样条函数内插，以统一到MODIS LST的观测时间。然后根据读取的GLDAS数据和MODIS数据的经纬度信息，在MATLAB软件平台中进行空间匹配，生成对应位置关系的查找表，最后以此为基础逐年输出相关参数数据。
第三，分别将250m的DEM和500m的MODIS土地覆盖数据重采样到1km分辨率，以匹配MODIS LST数据。16天1km分辨率的NDVI数据则进行时间插值，以符合逐日的时间分辨率要求。
最后，利用基于统计的时间滤波器填补由云层覆盖导致缺失的1km反照率和NDSI数据。
2.时间算法的实现
1)根据公式(1)，在Matlab平台上对1年的数据逐像元训练MODIS LST和GLDAS数据的回归关系。
2)根据(1)的回归结果，利用公式(2)估算初始的全天候1km地表温度。
3)通过公式(3)估算初始的全天候地表温度和MODIS LST之间的系统误差。
4)通过方程式(4)估算优化的全天候1km LSTT-RF。
3.空间算法的实现
1)对于一年的数据，逐日使用1km NDVI和NDSI将研究区域划分为若干子区域，包括人口密集区、稀疏区、荒地面积、雪冰区和水体面积。
2)对于每个子区域，将空间因子统一重采样到0.25°分辨率，以匹配GLDAS LST。然后，通过公式(5)分别训练GLDAS LST和晴空和非晴空区域的空间因子之间的回归关系。在此步骤中，仅选择包含80％及以上晴空条件子像元的GLDAS像元进行训练。
3)通过公式(6)估计初始全天候1km地表温度。
4)通过公式(7)估计初始全天候地表温度和MODIS LST之间的系统误差。
5)通过公式(8)估计优化的全天候1km LSTS-RF。
4.生成最终的1km全天候地表温度
1)根据公式(9)，对于一年的数据，在晴空条件下使用LSTT-RF、LSTS-RF和MODIS LST逐像元训练基于RF的回归函数贡献度。
2)利用LSTT-RF和LSTS-RF，分别将训练的函数应用于晴空和非晴空条件下，最终得到1km全天候LSTT-RF-S。
进一步地，具体的结果分析如下：
以2014年的结果为例，使用本发明的方法集成Aqua MODIS和GLDAS数据得到高分辨率的全天候地表温度。根据图2-3的结果，其中图2是包含研究区在内的地表温度图，其中(a)-(d)是2014年1月1日，4月1日，7月1日和10月1日的再分析资料(GLDAS)0.25°地表温度，(e)-(h)是MODIS地表温度，(i)-(l)是本发明提出方法得到的地表温度；图3是2014年MODIS地表温度和本发明提出方法得到的地表温度之间的空间分布评价指数图，其中(a)-(b)是0.25°再分析资料(GLDAS)的R2(相关系数)；(c)-(d)是本发明方法得到地表温度(LSTT-RF-S)的R2；(e)-(f)是0.25°再分析资料(GLDAS)的标准误差；(g)-(h)是本发明方法得到地表温度的标准误差。LSTT-RF-S与MODIS LST在空间分布上非常接近，而且没有明显的系统误差，可以很好地弥补MODIS LST的不足。此外，LSTT-RF-S图像中没有出现斑块效应，说明本方法在从低分辨率再分析数据中获取中/高分辨率LST方面具有很好的效果。与GLDAS LST相比，LSTT-RF-S相关系数值在0.89～0.99之间，比GLDAS LST总体上高0.06～0.36，而标准误差值在0.61～1.69K之间，比GLDAS LST低1.71～1.19K，与MODIS LST的一致性较好。此外，在夜间的LSTT-RF-S和MODIS LST的一致性要优于白天，其原因是根据MODIS观测时间进行时间内插的地表参数在白天的变化幅度更大，使结果具有一定的不确定性。另外，比较青藏高原以外的地区，LSTT-RF-S同MODIS LST之间表现出比在青藏高原更好的一致性。
根据图4，图4是时间算法和空间算法分别在2014年1月1日，4月1日，7月1日和10月1日得到的地表温度图，其中(a)-(d)是时间算法得到温度，(e)-(h)是空间算法得到温度。两种算法中，LSTT-RF的空间格局更接近MODIS LST，但两种算法的结果总体上十分接近。其中，LSTT-RF的相关系数略高于LSTS-RF，标准误差略低于LSTS-RF-。LSTT-RF-S与MODIS LST的一致性优于两种算法单独比较。这表明本方法能够根据RF回归机制动态确定每个LSTT-RF和LSTS-RF的最优贡献度。此外，时间算法在像元尺度对LSTT-RF-S的年内平均贡献如图5所示。图5是时间算法对最终结果的年内平均贡献度，其中(a)代表白天，(b)代表夜间。
综合在白天和夜间的情况，LSTT-RF对最终LSTT-RF-S的贡献超过71％，说明时间算法在这两种算法中占据主导地位。
图6展示的是2010年LSTT-RF-S和TIR LST之间的散点图，包括MODIS LST和AATSRLST。即图6是2010年本发明温度结果和TIR LST之间的散点图，其中(a)-(b)是晴空条件下本发明温度结果与MODIS LST之间的散点图；(c)-(d)是晴空下本发明温度结果与AATSRLST之间的散点图。LSTT-RF-S的平均误差为0.08K～0.16K，在1.12K～1.46K范围内。LSTT-RF-S和MODIS LST之间的一致性较好，表明了前一年建立的RF回归模型的良好适应性。与AATSRLST相比，平均误差为-0.21K～0.25K，标准误差为1.27K～1.36K，说明该模型在不同的星载热红外传感器上具有较好的精度。
图7比较了MODIS LST、集成热红外和微波遥感数据的LST以及LSTT-RF-S。图7是3种地表温度产品在2010年4月1日的温度图，其中(a)、(b)、(c)分别为MODIS LST、集成热红外与微波数据得到的地表温度以及本发明得到地表温度。
LSTT-RF-S和集成热红外和微波遥感数据的LST在大部分地区都表现出了空间分布的一致性，但后者由于微波数据热采样深度的影响，导致了MWLST和TIR LST之间存在较大的不一致性。另外，集成热红外和微波遥感数据的方法会对AMSR-E的条带进行填充，填补的区域存在明显的空间不连续性。与之相比本方法在视觉上具有更好的图像质量。
图8-10分别展示了3种地表温度产品与9个地面站点实测地表温度之间的散点图。图8中，(a)代表是阿柔站，(b)是大满站，(c)是DSL。图9，中(a)代表花寨子站，(b)代表D66，(c)是D105。图10中，(a)是GZ，(b)是玛曲，(c)是PRD。
在青藏高原包含的所有的站点中，MODIS LST和LSTT-RF-S的精度接近于MBE/RMSE，误差不超过0.10K/0.14K，这表明LSTT-RF-S在不同时间跨度上保持了良好的鲁棒性。所有站点的检验结果均表明LSTT-RF-S的精度最高，平均误差为-0.06K，均方根误差为2.14K。
本发明涉及前人关注较少的领域，提出了一种方法来产生高分辨率的全天候LST。本发明的核心是从时间维的再分析数据中构造地表温度的时间序列与其多个因子之间的回归关系，即时间算法。考虑到该算法对LSTT-RF-S具有较大的贡献度，经过地面站点实测温度数据的检验可以表明：在没有高分辨率输入数据的条件下，依靠较低分辨率的多时相地表温度相关参数，通过回归映射就能得到精度和质量可靠的全天候高时空分辨率的地表温度。在这一结论的启发下，可以进一步构造多个类似的模型函数，得到在云层覆盖下无法获得的全天候高分辨率地球物理参数，如反照率、LAI、NDVI等，从而有利于提高现有陆面过程模型的时空分辨率。
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围。凡采用等同替换或等效替换，这些变化是显而易见，一切利用本发明构思的发明创造均在保护之列。

Claims (4)

1.一种集成热红外与再分析资料的全天候地表温度生成方法，其特征在于，该方法首先从众多机器学习算法中选择随机森林回归(random forest regression)用于构建热红外遥感地表温度产品TIR LST和再分析资料之间的映射关系，并基于时间和空间维度提出两种算法用以刻画和描述地表温度的特征：一种是时间算法(temporal algorithm)，即描绘时间维度上地表温度变化的算法，该算法利用晴空条件下的低分辨率再分析数据生成全天候高分辨率的地表温度，最终生成的地表温度称为LSTT-RF；同时，考虑到地表温度在时空维度上变化的差异性，针对性地提出一种空间算法(spatial algorithm)，用以描述地表温度与其相关的空间因子之间的映射关系，该算法利用全天候条件下高空间分辨率数据，结合较低空间分辨率的再分析资料，生成全天候高空间分辨率地表温度，通过该算法得到的地表温度称为LSTS-RF，最终，通过随机森林回归的方法，结合所有与地表温度相关的时空因子，集成时间算法和空间算法的结果，生成精度最好的1km全天候地表温度LSTT-RF-S，该方法具体包括以下步骤：
1)数据预处理与时空匹配：
选择区域为青藏高原和其周边地区，即19.5°N-44°N,73°E-104°E，青藏高原是世界上海拔最高的高原，平均海拔超过4000米，由于其独特的地理位置和地形，这一地区具有非常复杂的气候条件和各种土地覆盖类型，青藏高原的动力、热力效应对亚洲季风的形成具有关键性的作用，同时也是全球气候变化的敏感区域之一；地表温度是地表能量平衡中的主要参数，对于进一步研究该地区的陆地-大气相互作用具有重要意义，因此迫切需要一种长时间序列、高时空分辨率的地表温度数据；具体采用的数据集包括：1)遥感数据集：MODIS逐日地表温度产品MOD21A1/MYD21A1，1km，16天合成归一化植被指数NDVI产品MOD12A2，1km，逐日地表反照率产品MCD43A3，1km，逐日归一化雪指数NDSI产品MOD10A1，500m，DEM 250m；2)再分析资料：全球陆面数据同化系统GLDAS，时间分辨率3小时，空间分辨率为0.25°GLDAS_NOAH025_3H；3)验证数据集：9个科研站点的地表实测温度数据；
首先，利用IDL语言和MRT工具对MODIS数据产品进行批量处理；对于MOD21A1/MYD21A1数据，利用相应的质量控制数据，提取高质量的MODIS LST像元；
其次，对GLDAS数据的参数时间进行三次样条函数内插，以统一到MODIS LST的观测时间，然后根据读取的GLDAS数据和MODIS数据的经纬度信息，在MATLAB软件平台中进行空间匹配，生成对应位置关系的查找表，最后以此为基础逐年输出相关参数数据；
第三步，分别将250m的DEM和500m的MODIS土地覆盖数据重采样到1km分辨率，以匹配MODIS LST数据；16天1km分辨率的NDVI数据则进行时间插值，以符合逐日的时间分辨率要求；
最后，利用基于统计的时间滤波器填补由云层覆盖导致缺失的1km反照率和NDSI数据；
2)时间算法的实现：
第一步，根据公式(1)，将预处理好的1年再分析资料GLDAS和MODIS数据在Matla b平台上逐像元训练MODIS LST和GLDAS数据的回归关系：
Figure FDA0002594835390000021
式中，LSTTIR-Q表示高空间分辨率的卫星地表温度产品中像元Q对应的时间序列，单位为K，预设时长为1年；GLsd代表包含像元Q的再分析数据相关因子的时间序列；RFQ表示LSTTIR-Q与GLsd之间建立的随机森林回归关系；PT-m-sd表示第m个地表温度相关因子的时间序列；
第二步，根据第一步的回归结果，利用公式(2)估算初始的全天候1km地表温度，考虑云层覆盖的影响，将利用公式(1)求取初始1km全天候地表温度的过程具体到晴空和非晴空的条件下：
Figure FDA0002594835390000022
式中，
Figure FDA0002594835390000023
和
Figure FDA0002594835390000024
分别表示在晴空和非晴空条件下，像元所对应的1km地表温度估算值；td-cd-n表示非晴空条件下的第n个像元；
第三步，通过公式(3)估算初始的全天候地表温度和MODIS LST之间的系统误差，具体为像元Q对应的1km地表温度估算值和热红外遥感地表温度产品之间的系统误差可以表示为：
Figure FDA0002594835390000031
式中，
Figure FDA0002594835390000032
和
Figure FDA0002594835390000033
分别代表像元Q对应的1km地表温度估算值和热红外遥感地表温度产品的平均值，ΔδT-Q表示像元Q对应的系统误差；
第四步，通过方程式(4)估算优化的全天候1km LSTT-RF，具体为通过分别优化晴空和非晴空条件下地表温度估计值，可以获得最终的全天候地表温度LSTT-RF：
Figure FDA0002594835390000034
式中，LSTT-RF-sd-Q和LSTT-RF-cd-Q分别代表晴空和非晴空条件下像元Q对应的1km地表温度值，LSTT-RF-Q表示像元Q对应的最终的LSTT-RF值；
3)空间算法的实现：
第一步，对于1年的数据，逐日使用1km NDVI和NDSI将研究区域划分为若干子区域；
第二步，对于每个子区域，将空间因子统一重采样到0.25°分辨率，以匹配GLDAS L ST；然后，通过公式(5)分别训练GLDAS LST和晴空和非晴空区域的空间因子之间的回归关系，其中，非晴空区域的空间因子包括NDVI，NDSI，DEM，反照率；在此步骤中，仅选择包含80％及以上晴空条件子像元的GLDAS像元进行训练，具体为对于包含像元P的特定区域，在再分析数据的像元尺度下，地表温度可以与其相关的空间因子的关系可以表示为公式(5)：
Figure FDA0002594835390000035
式中，其中LSTS-low表示从再分析资料中得到的低分辨率地表温度数据，单位为K；NS-k-low表示对应LSTS-low分辨率的第k个与地表温度相关的空间因子；RFs表示像元P对应的地表温度与多个与其相关的空间因子之间的基于随机森林回归的函数模型；
第三步，通过公式(6)估计初始全天候1km地表温度，将公式(4)应用到1km的高空间分辨率情况下，则可通过公式(6)估算像元P的初始1km地表温度：
Figure FDA0002594835390000041
式中，
Figure FDA0002594835390000042
在表示较高分辨率情况下像元P对应的地表温度估计值，单位K；NS-k-high-P表示像元P对应的第k个较高空间分辨率的地表温度空间因子；RFs表示较高分辨率下像元P对应的地表温度与多个相关空间因子之间随机森林回归关系；
第四步，通过公式(7)估计初始全天候地表温度和MODIS LST之间的系统误差，具体为像元P对应的1km地表温度估算值和热红外遥感地表温度产品之间的系统误差可以表示为公式(7)：
Figure FDA0002594835390000043
式中，
Figure FDA0002594835390000044
代表像元P对应的1km地表温度估算值的平均值，单位为K；
Figure FDA0002594835390000046
代表像元P对应的1km地表温度产品的平均值，单位为K；
第五步，通过公式(8)估计优化的全天候1km LSTS-RF，具体为通过分别优化晴空和非晴空条件下地表温度估计值，可以获得最终的全天候地表温度LSTS-RF：
Figure FDA0002594835390000045
式中，LSTS-RF-sd-P和LSTS-RF-cd-P分别代表晴空和非晴空条件下像元P对应的1km地表温度值，LSTS-RF-P表示像元P对应的最终LSTS-RF值；
4)生成最终的1km全天候地表温度：
第一步，根据公式(9)，对于1年的数据，在晴空条件下使用LSTT-RF、LSTS-RF和MO DIS LST逐像元训练基于RF的回归函数贡献度：
LSTT-RF-S-P＝RFM(LSTT-RF-P,LSTS-RF-P) (9)
其中RFM表示基于随机森林的回归函数，用于衡量像元P对应的从时间和空间算法生成的地表温度对最终地表温度的贡献度；
第二步，利用LSTT-RF和LSTS-RF，分别将训练的函数应用于晴空和非晴空条件下，最终得到1km全天候LSTT-RF-S。
2.根据权利要求1所述的集成热红外与再分析资料的全天候地表温度生成方法，其特征在于，所述步骤1)中所述批量处理包括格式转换、投影变换、拼接和裁剪。
3.根据权利要求2所述的集成热红外与再分析资料的全天候地表温度生成方法，其特征在于，所述步骤2)中所述m取值为12。
4.根据权利要求1-3任意一项所述的集成热红外与再分析资料的全天候地表温度生成方法，其特征在于，所述步骤3)中所述若干子区域包括人口密集区、稀疏区、荒地面积、雪冰区和水体面积。
CN202010707434.0A 2020-07-21 2020-07-21 集成热红外与再分析资料的全天候地表温度生成方法 Active CN111982294B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010707434.0A CN111982294B (zh) 2020-07-21 2020-07-21 集成热红外与再分析资料的全天候地表温度生成方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010707434.0A CN111982294B (zh) 2020-07-21 2020-07-21 集成热红外与再分析资料的全天候地表温度生成方法

Publications (2)

Publication Number Publication Date
CN111982294A true CN111982294A (zh) 2020-11-24
CN111982294B CN111982294B (zh) 2022-06-03

Family

ID=73438368

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010707434.0A Active CN111982294B (zh) 2020-07-21 2020-07-21 集成热红外与再分析资料的全天候地表温度生成方法

Country Status (1)

Country Link
CN (1) CN111982294B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113158570A (zh) * 2021-04-26 2021-07-23 电子科技大学 一种融合多源卫星遥感的全天候地表温度近实时反演方法
CN113255148A (zh) * 2021-06-04 2021-08-13 中国科学院地理科学与资源研究所 基于modis产品数据估算全天气气温及其时空分布方法
CN113743000A (zh) * 2021-08-13 2021-12-03 电子科技大学 一种生成高时间分辨率全天候地表温度的方法
CN113834572A (zh) * 2021-08-26 2021-12-24 电子科技大学 一种无人机非制冷型热像仪测温结果漂移去除方法
CN114047563A (zh) * 2021-10-18 2022-02-15 中国人民解放军国防科技大学 一种红外高光谱的全天候同化方法
CN114417728A (zh) * 2022-01-27 2022-04-29 中国农业科学院农业资源与农业区划研究所 基于温度和发射率及深度学习的近地表空气温度反演方法
CN115859211A (zh) * 2022-11-17 2023-03-28 中国农业科学院农业资源与农业区划研究所 一种基于三温不确定度估算模型的地表温度产品融合方法
CN116933664A (zh) * 2023-09-15 2023-10-24 航天宏图信息技术股份有限公司 地表温度数据重构模型构建方法及地表温度数据重构方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103983360A (zh) * 2014-05-30 2014-08-13 中国科学院遥感与数字地球研究所 一种基于hj-1b irs卫星数据反演lst的方法
CN105204024A (zh) * 2015-10-12 2015-12-30 电子科技大学 微波遥感表面温度向热红外遥感肤面温度转换的方法
CN107576417A (zh) * 2017-09-04 2018-01-12 电子科技大学 一种全天候地表温度生成方法
WO2018145229A1 (zh) * 2017-02-10 2018-08-16 广西壮族自治区气象减灾研究所 一种近地面气温的大面积精确反演方法
CN109580003A (zh) * 2018-12-18 2019-04-05 成都信息工程大学 一种静止气象卫星热红外数据估算近地面大气温度方法
CN110516816A (zh) * 2019-08-30 2019-11-29 中国科学院、水利部成都山地灾害与环境研究所 基于机器学习的全天候地表温度生成方法及装置
US20200118012A1 (en) * 2017-04-18 2020-04-16 Hewlett-Packard Development Company, L.P. Monitoring the Thermal Health of an Electronic Device
US20200167640A1 (en) * 2018-11-27 2020-05-28 The Boeing Company System and method for generating an aircraft fault prediction classifier

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103983360A (zh) * 2014-05-30 2014-08-13 中国科学院遥感与数字地球研究所 一种基于hj-1b irs卫星数据反演lst的方法
CN105204024A (zh) * 2015-10-12 2015-12-30 电子科技大学 微波遥感表面温度向热红外遥感肤面温度转换的方法
WO2018145229A1 (zh) * 2017-02-10 2018-08-16 广西壮族自治区气象减灾研究所 一种近地面气温的大面积精确反演方法
US20200118012A1 (en) * 2017-04-18 2020-04-16 Hewlett-Packard Development Company, L.P. Monitoring the Thermal Health of an Electronic Device
CN107576417A (zh) * 2017-09-04 2018-01-12 电子科技大学 一种全天候地表温度生成方法
US20200167640A1 (en) * 2018-11-27 2020-05-28 The Boeing Company System and method for generating an aircraft fault prediction classifier
CN109580003A (zh) * 2018-12-18 2019-04-05 成都信息工程大学 一种静止气象卫星热红外数据估算近地面大气温度方法
CN110516816A (zh) * 2019-08-30 2019-11-29 中国科学院、水利部成都山地灾害与环境研究所 基于机器学习的全天候地表温度生成方法及装置

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
WEI ZHAO 等: "Normalization of the temporal effect on the MODIS land surface temperature product using random forest regression", 《ISPRS JOURNAL OF PHOTOGRAMMETRY AND REMOTE SENSING》 *
华俊玮 等: "基于随机森林算法的地表温度降尺度研究", 《国土资源遥感》 *
权凌 等: "基于时间序列建模的城市热岛时间尺度成分分离方法与应用", 《地球科学进展》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113158570A (zh) * 2021-04-26 2021-07-23 电子科技大学 一种融合多源卫星遥感的全天候地表温度近实时反演方法
CN113255148A (zh) * 2021-06-04 2021-08-13 中国科学院地理科学与资源研究所 基于modis产品数据估算全天气气温及其时空分布方法
CN113743000A (zh) * 2021-08-13 2021-12-03 电子科技大学 一种生成高时间分辨率全天候地表温度的方法
CN113743000B (zh) * 2021-08-13 2023-04-18 电子科技大学 一种生成高时间分辨率全天候地表温度的方法
CN113834572A (zh) * 2021-08-26 2021-12-24 电子科技大学 一种无人机非制冷型热像仪测温结果漂移去除方法
CN113834572B (zh) * 2021-08-26 2023-05-12 电子科技大学 一种无人机非制冷型热像仪测温结果漂移去除方法
CN114047563A (zh) * 2021-10-18 2022-02-15 中国人民解放军国防科技大学 一种红外高光谱的全天候同化方法
CN114417728A (zh) * 2022-01-27 2022-04-29 中国农业科学院农业资源与农业区划研究所 基于温度和发射率及深度学习的近地表空气温度反演方法
CN115859211A (zh) * 2022-11-17 2023-03-28 中国农业科学院农业资源与农业区划研究所 一种基于三温不确定度估算模型的地表温度产品融合方法
CN115859211B (zh) * 2022-11-17 2023-07-21 中国农业科学院农业资源与农业区划研究所 一种基于三温不确定度估算模型的地表温度产品融合方法
CN116933664A (zh) * 2023-09-15 2023-10-24 航天宏图信息技术股份有限公司 地表温度数据重构模型构建方法及地表温度数据重构方法
CN116933664B (zh) * 2023-09-15 2024-01-05 航天宏图信息技术股份有限公司 地表温度数据重构模型构建方法及地表温度数据重构方法

Also Published As

Publication number Publication date
CN111982294B (zh) 2022-06-03

Similar Documents

Publication Publication Date Title
CN111982294B (zh) 集成热红外与再分析资料的全天候地表温度生成方法
Vander Jagt et al. The effect of spatial variability on the sensitivity of passive microwave measurements to snow water equivalent
CN111553245A (zh) 基于机器学习算法和多源遥感数据融合的植被分类方法
CN111678866A (zh) 一种多模型集成学习的土壤含水量反演方法
CN107014753B (zh) 作物长势监测方法和系统
CN106021868A (zh) 一种基于多规则算法的遥感数据降尺度方法
Liu et al. Estimation of surface and near-surface air temperatures in arid northwest china using landsat satellite images
CN114120101A (zh) 一种土壤水分多尺度综合感知方法
Sun et al. Microwave and meteorological fusion: A method of spatial downscaling of remotely sensed soil moisture
CN113447137A (zh) 一种面向无人机宽波段热像仪的地表温度反演方法
Holmlund et al. Constraining 135 years of mass balance with historic structure-from-motion photogrammetry on Storglaciären, Sweden
Wang et al. Estimation of 30 m land surface temperatures over the entire Tibetan Plateau based on Landsat-7 ETM+ data and machine learning methods
Petropoulos et al. Exploring the use of Unmanned Aerial Vehicles (UAVs) with the simplified ‘triangle’technique for soil water content and evaporative fraction retrievals in a Mediterranean setting
CN111401644A (zh) 一种基于神经网络的降水降尺度空间预测方法
Bartolacci et al. Patterns of co-variability between physical and biological parameters in the Arabian Sea
El Garouani et al. Water balance assessment using remote sensing, Wet-Spass model, CN-SCS, and GIS for water resources management in Saïss Plain (Morocco)
CN107688712B (zh) 一种基于dem和ndvi的气温降尺度方法
KR102608357B1 (ko) 위성 기반 해수면 염분 데이터와 해양 수치 모델 기반 해수면 염분 데이터를 머신러닝 모델에 적용하여 해수면 염분을 예측하기 위한 방법 및 장치
CN109357765B (zh) 一种面向土壤属性预测与制图的亮温或陆面温度协同变量构建选择方法
Georgiou et al. Evaluation of MODIS‐derived LST products with air temperature measurements in Cyprus
CN113762383B (zh) 一种基于多源数据的植被指数融合方法
CN115993668B (zh) 一种基于多项式改正和神经网络的pwv重建方法及系统
Huss et al. Temporal and spatial changes of Laika Glacier, Canadian Arctic, since 1959, inferred from satellite remote sensing and mass-balance modelling
Kiyofuji et al. A visualization of the variability of spawning ground distribution of Japanese common squid (Todarades pacificus) using Marine-GIS and satellite data sets
Wu et al. Global statistical assessment of Haiyang-2B scanning microwave radiometer precipitable water vapor

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