一种基于多时间尺度模型的遥感地表日均温计算方法与流程

本发明涉及一种基于多时间尺度模型的遥感地表日均温计算方法，属于遥感地表温度时间插值模型应用及遥感地表温度产品生产研究领域。

背景技术：

地表温度作为地表-大气能量交换过程中的关键参数，时刻处在动态变化当中。热红外遥感是快速大范围监测地表温度变化的重要手段。目前，热红外遥感动态监测地表温度面临以下问题需要解决：(1)受到云及大气因素的影响，地表温度观测存在大范围的缺值；(2)所观测到的地表温度均是瞬时值，无法反映地表温度的连续变化状况；(3)热红外遥感卫星产品的时间和空间分辨率通常无法兼顾。受限于上述问题，当前仍然缺乏空间尺度在公里级且时空连续的地表温度产品。目前，学术界通常采用8天平均或者月平均的方式来分析地表温度的动态变化特征。但是，这样一方面损失了地表温度逐日尺度的变化信息，另一方面几个瞬时过境时刻的简单平均与真实准确的日均温相比会存在系统性偏差。

为了生产空间尺度为公里级且真实准确的日均温产品，需要克服两个难题，一是如何重建云下瞬时温度，另一个是如何从简单几次的瞬时观测中构建出连续的地表温度日变化过程。云下地表温度重建方法可以大致分为统计法和模型法这两大类。统计法注重寻找地表温度与地面高程模型(digitalelevationmodel，dem)、归一化植被指数(normalizeddifferencevegetationindex，ndvi)、空气温度、微波温度观测等数据之间的相关性。模型法则多利用地表温度年变化模型(annualtemperaturecycle,atc)等时间插值模型，反距离权重插值等空间插值模型，以及时空插值模型来重建云下观测温度。从瞬时观测到连续的地表温度日变化过程构建的主要方法包括地表温度降尺度、图像融合及地表温度日变化模型(diurnaltemperaturecycle,dtc)。地表温度降尺度方法的原理为在低分辨率影像上建立地表温度与降尺度回归核(如ndvi、dem)之间的关系，再将这一关系映射到高分辨率影像上。图像融合方法通过各类融合算法将静止气象卫星数据高时间分辨率和极轨卫星高空间分辨率的特点结合起来。dtc模型则在晴空假设的基础上基于若干次瞬时地表温度观测来模拟地表温度日变化过程。

云下地表温度重建和从瞬时温度到连续的地表温度日变化过程构建，这两类方法之间的组合从理论上均可以实现地表日均温的生产。但是，在实际生产过程中，特别是大范围生产公里级别的产品，目前仍然缺乏一种能够基于以中分辨率成像光谱仪(moderateresolutionimagingspectroradiometer，modis)为代表的双极轨卫星地表温度数据快速生产日均温的方法，具体体现在：(1)在云下温度重建过程中，大多数云下重建方法仅适用于局地小区域，受限于计算速度，难以推广到全球等大范围产品生产；(2)在从瞬时温度到连续的地表温度日变化过程构建中，降尺度与图像融合方法受限于不同卫星数据之间的不完全匹配的系统性误差，从而影响到日均温产品的一致性。因此，如何在保证精度的前提下，用一种简单且快速的方法大范围生产日均温产品需要进一步的探讨。

技术实现要素：

本发明的目的在于：为实现快速准确生产地表温度日均温产品，在多时间尺度模型的框架下，提出一种耦合地表温度年变化与地表温度日变化模型的日均温产品生产方法。

本发明提出的基于多时间尺度模型的遥感地表日均温计算方法，包括以下步骤：

第一步、数据获取——获取modis地表温度产品，得到所有像元的每天四个过境时刻的地表温度观测值，当过境时刻像元上空有云遮挡时，则该像元在该过境时刻无对应的无地表温度观测值；

第二步、计算云下像元的地表温度初始拟合结果——针对每个过境时刻，将像元的地表温度观测值输入oatc模型，求解oatc模型三个自由参数，并利用参数确定oatc模型用于计算云下像元地表温度的初始拟合结果，oatc模型公式如下：

式中，toatc(d)表示第d天的地表温度初始拟合值；n表示一年中的总天数，闰年为366，平年为365；ta表示地表温度年均温；a表示地表温度一年中的变化幅度；β表示地表温度年变化相位角，ta、a和β是oatc模型待求解的三个自由参数；

第三步、计算云下像元地表温度的最终拟合结果——计算对应的空气温度观测值与地表温度初始拟合结果的差δtair(d)，针对每个境时刻，将像元的地表温度观测值和δtair(d)输入eatc模型，求解eatc模型四个自由参数，并利用参数确定eatc模型计算的云下像元地表温度的最终拟合结果，eatc模型公式如下：

式中，teatc(d)表示第d天的地表温度最终拟合值；n表示一年中的总天数，闰年为366，平年为365；ta表示地表温度年均温；a表示地表温度一年中的变化幅度；β表示地表温度年变化相位角；k是空气温度与地表温度之间调节参数，ta,a,β和k是四个待求解的自由参数；

第四步、计算日均地表温度——无云遮挡时，地表温度为地表温度观测值，有云遮挡时，地表温度为eatc模型最终拟合结果，将四个过境时刻及其对应的地表温度输入到四参数dtc模型中，逐像元拟合一天中地表温度曲线，再根据拟合得到的地表温度曲线计算日均地表温度。

本发明可以重建云遮挡造成的地表温度缺值，同时能够有效纠正传统四次瞬时观测平均较真实日均温偏高的系统性偏差，获得时间和空间上均连续的日均温产品。本发明方法适用性强，能用于大范围遥感地表日均温产品的快速生产，能够满足实际生产的需要。

附图说明

下面结合附图对本发明作进一步的说明。

图1为本发明方法流程示意图。

图2为本发明基于modis数据，在编号为h11v04，2017年第180天(2017年6月29号)的结果示意图。

具体实施方式

如图1所示，本发明基于modis数据生产日均温产品的过程如下：

第一步：数据获取。

获取modis空间分辨率为1千米的地表温度产品，产品编号为mod11a1和myd11a1。如图1a部分所示，modis一天中有四个过境时刻(terra白天、aqua白天、terra夜晚、aqua夜晚)。

第二步、计算云下像元的地表温度初始拟合结果。

针对每个境时刻，将像元的地表温度观测值输入oatc模型，求解oatc模型三个自由参数，并利用参数确定oatc模型用于计算云下像元地表温度的初始拟合结果，oatc模型公式如下：

式中，toatc(d)表示第d天的地表温度初始拟合值；n表示一年中的总天数，闰年为366，平年为365；ta表示地表温度年均温；a表示地表温度一年中的变化幅度；β表示地表温度年变化相位角，ta、a和β是oatc模型待求解的三个自由参数。

以terra白天过境时刻为例，将像元terra白天过境时刻的地表温度观测数据输入到oatc模型中，使用最小二乘法拟合求解oatc模型的三个自由参数，各参数的初始值及拟合的上下界如表1所示：

表1oatc模型参数拟合时的初始值和上下界

*tmax表示输入地表温度数据中最大值，tmin表示输入地表温度数据中最小值。

然后利用确定后的模型计算像元的地表温度初始拟合值。

第三步：计算云下像元的地表温度最终拟合结果。

计算对应的空气温度观测值(在小尺度生产中，可以采用站点数据观测结果，而在大尺度生产中，空气温度观测值可以利用再分析数据模拟得到的空气温度结果。计算空气温度观测值与地表温度初始拟合值的差δtair(d)，针对每个过境时刻，将像元的地表温度观测值和δtair(d)输入eatc模型，求解eatc模型四个自由参数，并利用参数确定eatc模型计算的云下像元地表温度的最终结果，eatc模型公式如下：

式中，teatc(d)表示第d天的地表温度最终拟合值；n表示一年中的总天数，闰年为366，平年为365；ta表示地表温度年均温；a表示地表温度一年中的变化幅度；β表示地表温度年变化相位角；k是空气温度与地表温度之间调节参数，ta,a,β和k是四个待求解的自由参数。四个待求解的自由参数，通过最小二乘法拟合确定，将求解得到的参数带入到公式(2)中即可得到eatc模型拟合得到的地表温度年变化结果(如图1b部分所示)。eatc模型的参数初值及拟合上下界见下表：

表2eatc模型参数拟合时的初始值和上下界

*tmax表示输入地表温度数据中最大值，tmin表示输入地表温度数据中最小值。

第四步：对四个过境时刻分别重复第二步和第三步，这样利用eatc模型能够逐像元重建出四个过境时刻一年的云下地表温度。其中，eatc模型重建的晴空观测部分保留modis原始的晴空观测结果，云下观测部分使用eatc模型的重建结果，旨在保证modis各过境时刻晴空观测部分的精度。eatc模型云下重建这一步使得modis每个像元每天四个过境时刻均有地表温度观测值，这为下一步dtc模型的使用奠定基础(如图1c部分所示)。

第五步：计算日均地表温度。

无云遮挡时，地表温度为地表温度观测值，有云遮挡时，地表温度为eatc模型最终拟合结果，将四个过境时刻及其对应的地表温度输入到四参数dtc模型中，逐像元拟合一天中地表温度曲线，再根据拟合得到的地表温度曲线计算日均地表温度。

本实施例中，通过线性插值计算出云下过境时刻的时间。第三步中eatc模型重建出了各过境时刻的云下地表温度结果，但mod11a1/myd11a1产品中并没有云下缺值像元的过境时刻，由于modis过境时刻波动范围并不是很大，在此处对某个像元一年内缺失的过境时间采用线性插值的方式来估算各过境时刻云下的过境时间。

在eatc模型与过境时间线性插值的基础上，利用四参数dtc模型逐像元拟合地表温度的日内变化情况。将modis四个过境时刻的温度及其过境时间输入到本发明使用的四参数dtc模型(命名为got09-dt-001)中。需要注意的是，dtc模型的拟合过程是从当天的日出阶段开始，到第二天的日出时刻。通常情况下modisaqua卫星的夜间过境时刻(大约在当地时间01：30附近)在当天日出之前，因此在输入过程中通常将当天aqua卫星夜间过境时刻数据作为第二天aqua卫星夜间过境结果用于拟合dtc模型(如图1d部分所示)。本发明使用的dtc模型如公式(3)所示。

式中，t0表示地表温度残差项；ta表示白天温度余弦变化的振幅；tm表示一天中地表温度达到最大的时刻；ts表示地表开始自然冷却的时刻；θ表示热时角；热时角对应的“太阳天顶角”为θz；m(θz)为大气质量；re表示地球半径，取为6371千米；h表示大气标高，估算后定为8.43千米。θ、θz、m(θz)、k值的计算方法参考andolesen(2009)即可。当t＝tm或者t＝ts时，θz分别取值为θz,min,θzs，m(θz)分别取值为mmin,m(θzs)。got09-dt-001模型包含t0,ta,tm,ts四个参数，四个参数通过最小二乘法拟合确定，将计算得到的参数结果代入公式(3)中即可拟合出地表温度的日变化状况(如图1d部分所示)。got09-dt-001模型各参数的初始值及拟合的上下界如下表所示：

表3got09-dt-001模型参数拟合时的初始值和上下界

*tmin,tmax,tss分别表示地表温度输入数据中的最大值、最小值、以及日落时刻。

在got09-dt-001模型拟合出地表温度的日变化状况的基础上，从日出阶段开始，通过积分的手段计算出地表温度的日均温。本例中在一天中每隔5分钟取一个温度然后平均作为最终计算得到日均温，也可以是每隔60分钟取一个温度然后平均，或者30分钟取一个温度然后平均。得到日均温产品的同时也可以计算在eatc模型重建基础上四次过境时刻的平均温度减去日均温的系统性误差(如图1e部分所示)。

图2展示了编号为h11v04的modis文件在2017年第180天的生产结果。图2a～d表示modis四个过境的原始观测，图2e～h分别表示对应过境时刻eatc模型的云下重建结果，图2i表示日均温生产结果，图2j表示四次过境时刻的平均温度减去日均温的系统性误差结果，图2k表示示例区的地理位置及地表覆盖类型状况。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。