CN116227142A - 基于modis遥感数据的逐日地表温度重建方法及装置 - Google Patents

Abstract

本发明提供一种基于MODIS遥感数据的逐日地表温度重建方法及装置。该方法能够在个人终端或云平台(如：GoogleEarthEngine)上独立实现假定晴空条件下任意区域的每日地表温度重建。该模型首先采用最邻近时相线性插值对以云污染为首的缺失值进行初步填充，再通过临年同期数据的加权融合缩小缺失范围，对于存在的噪声点和过拟合问题，引入信号学中的SG滤波进行平滑和校正，最终生成逐日1km的无缝地表温度产品。验证结果表明，该发明获得的产品能够较好地反映假定晴空下地表温度的真实状态，更能在排除天气因素干扰后有效反映地表温度的季节性特征，在长时序分析研究中更具有可靠性。

Description

基于MODIS遥感数据的逐日地表温度重建方法及装置

技术领域

本发明涉及遥感技术领域，提供一种基于MODIS遥感数据的逐日地表温度重建方法及装置。

背景技术

地表温度(Land Surface Temperature，LST)指的是太阳辐射到达地面后，经反射和吸收作用后对地面进行测量的温度。与平时提到的气温不同，地表温度代表陆地表面与空气交界处的温度，而前者代表距离地表1.5米左右高度的空气流动温度。尽管有所差异，大气温度与地表温度却都是反映大气圈、岩石圈、生物圈等泛在地球系统的晴雨表。除了全球温度变化，作为陆地表层和大气能量交换的直接驱动因子，地表温度还被广泛应用于生态环境监测、作物蒸散发估算和城市热环境分布等研究领域。作为获取区域性或全球性地表温度产品最可靠的途径，热红外遥感反演拥有地面定点观测等传统方法难以比拟的技术优势，长期以来都是获取高时空分辨率地表温度产品的首选方式。热红外遥感依靠搭载在卫星上的热红外传感器通过非接触测量的方式来获取地表温度。

1960年以来，大约有80种热红外辐射计搭载于240多颗卫星上，用来获取地表温度、云顶温度、大气水汽，降水量等一系列地球参量，基于辐射传输方程的各种算法已经相对成熟，如单通道算法、多通道算法、多角度算法、多时相算法、高光谱反演算法。虽然热红外反演到的地表温度空间分辨率和精度较高，但却容易受云污染的影响，数据有效率平均只有1/3。全覆盖地表温度重建产品主要由时空数据插值、多源数据融合和物理模型推演三种方法得到。但这些方法都有一些不足，比如在西部等气象站点较为稀疏的地区同化数据的产品精度较低，也导致局部融合效果不理想；时空插值方法，计算结果并非真实地表温度，在大面积或长时序数据缺失的情况下效果会严重失真等。

发明内容

为了解决上述技术问题，本发明提供了一种基于MODIS遥感数据的逐日地表温度重建方法及装置，达到MODIS地表温度数据时空无缝的效果。

根据本发明的一个方面，本发明提供的一种基于MODIS遥感数据的逐日地表温度重建方法，包括以下步骤：

S1：获取MOD11A1地表温度产品，并对其进行预处理后，筛选合格地表温度像元；

S2：对MOD11A1地表温度产品的单日全球缺失影像进行全域缺失时相弥补，得到连续的MOD11A1地表温度产品；

S3：以天为单位，对MOD11A1地表温度产品所缺影像的前后影像集使用最邻近法进行线性插值处理；

S4：结合邻年同期数据，利用地表温度具有正弦谐波函数的变化规律对重建像元进行加权融合，从而减少线性插值过程中的误差；

S5：获取MCD12Q1地表分类产品，计算同类地物的地表温度均值，进一步填充与之同类地物的地表温度缺失像元；

S6：使用SG滤波去除重建过程中的噪声，最终得到逐日地表温度重建结果。

进一步地，步骤S1中，所述预处理包括：条件筛选、影像镶嵌、研究区裁剪，像元质量审查等预处理。

进一步地，步骤S1具体包括：

S1.1：建立感兴趣区范围，获取MOD11A1地表温度产品，将所有影像裁剪至与感兴趣区同一范围；

S1.2：根据受云影响程度、数据质量、发射率误差和反演误差四个指标进一步对原始MOD11A1地表温度产品的质量进行筛选，获取合格地表温度像元。

优选地，步骤S1.2中，选取地表温度平均反演误差小于或等于2K和发射率平均误差小于或等于0.04这两个阈值作为筛选合格地表温度像元的依据。

进一步地，步骤S2中，所述全域缺失时相弥补的具体方法为：

对缺失的MOD11A1单日数据进行临近复刻，即通过搜索时序产品中距离缺失时相最近的影像作为该时期的初始值。

进一步地，步骤S3中，所述使用最邻近法进行线性插值处理，具体包括：

利用相邻前后时相晴空地表温度像元的观测值使用线性函数对缺失像元进行插值，从而快速弥补缺失像元的地表温度值；

结合地表温度短期变化趋向随机性的特点和GEE像元级单元运算的方式，线性插值过程逐像元进行，所引入的晴空地表温度像元都距离缺失像元在时相上是最接近的；

以预测时刻t_m为分界线，将距离该时刻间隔为n天以内的影像作为临近时相的影像集并划分为前后两个组合，其包含的时段分别为t_m-n-t_m-1和t_m+1-t_m+n；

对于预测时刻前的影像集采用末位压缩的方法，将该时段每个位置观测时间最新的晴空地表温度像元整合为一张影像；

对于预测时刻后的影像集采用首位压缩的方法，将该时段每个位置观测时间最新的晴空地表温度像元整合为一张影像；

对于压缩后的前后两张影像，每一个像元都附有不同的时间信息，依据像元代表的地表温度因变量和时间自变量进行一元线性插值，计算得到预测时刻像元的地表温度，其表达式为：

其中，LST_pre为预测时刻t_m的像元值，LST_up和LST_down分别为地表温度波段前后压缩的像元值，t_up和t_down为像元所包含的时间信息，b为一元线性方程的截距。

进一步地，步骤S4中，所述结合邻年同期数据，利用地表温度具有正弦谐波函数的变化规律对重建像元进行加权融合，具体包括：

对于不同年份的LST，采用反时间距离权重对预测时刻前后多个同期数据进行加权融合，其表达式为：

/>

其中，LST_fus为加权融合后的地表温度，LST_y和ω_y代表预测时刻t_y的地表温度和权重，y为年份，LST_i/j和ω_i/j代表预测时刻t_i/j的地表温度和权重，k为距离预测年份的最大年间隔；

权重约束条件为：

进一步地，步骤S5中，所述获取MCD12Q1地表分类产品，计算同类地物的地表温度均值，进一步填充与之同类地物的地表温度缺失像元，具体包括：

按照MCD12Q1土地覆盖利用产品的IGBP分类标准，将研究区地物分为17类，每类地物的地表温度均值计算方法为：

其中，LST_mean是某一类地表温度均值，n为该地类的像元总数，LST_i为该地类第i个像元的地表温度值；

使用LST_mean进一步填充与之同类地物的LST缺失像元。

进一步地，步骤S6中，所述使用SG滤波去除重建过程中的噪声，具体包括：

首先确定一个半径为m的移动窗口，将窗口内的所有数据纳入一个集合，窗口内所包含的时序数据时间节点为(-m，-m+1，…，0，1，…，m)，采用k-1次多项式对窗口内2m+1个时间点进行拟合，其多项式表达为：

y＝a₀+a₁x+a₂x²+…+a_k-1x^k-1

其中，y为滤波后的像元值，a_i(i＝1,2...k-1)为多项式参数；

窗口内每一个时间点对应一个多项式，共有2m+1个，当窗口内的像元数n大于多项式参数k时，求得方程方程组的参数，其表达式为：

用矩阵形式表示为：

Y_(2m+1)×1＝X_(2m+1)×k·A_K×1+E_(2m+1)×1

其中，Y_(2m+1)×1为滤波后LST值，X_(2m+1)×k为参与滤波的每年LST值，A_K×1为参数矩阵，E_(2m+1)×1为随机误差向量；

拟合参数A的最小二乘解

为：

Y的模型预测值或滤波值

为：

对于MOD11A1时序地表温度产品，2m+1的窗口代表预测时刻前后m天内的时段，X代表加权融合后并进行地类均值填充的地表温度值，

代表SG滤波去除噪声后的LST值。

根据本发明的另一方面，本发明提供了一种基于MODIS遥感数据的逐日地表温度重建装置，包括以下模块：

预处理模块，用于获取MOD11A1地表温度产品，并对其进行预处理后，筛选合格地表温度像元；

全域缺失时相弥补模块，用于对MOD11A1地表温度产品的单日全球缺失影像进行全域缺失时相弥补，得到连续的MOD11A1地表温度产品；

最临近插值模块，用于以天为单位，对MOD11A1地表温度产品所缺影像的前后影像集使用最邻近法进行线性插值处理；

加权融合模块，用于结合邻年同期数据，利用地表温度具有正弦谐波函数的变化规律对重建像元进行加权融合，从而减少线性插值过程中的误差；

缺失像元填充模块，用于获取MCD12Q1地表分类产品，计算同类地物的地表温度均值，进一步填充与之同类地物的地表温度缺失像元；

SG滤波模块，用于使用SG滤波去除重建过程中的噪声，最终得到逐日地表温度重建结果。

本发明提供的技术方案具有以下有益效果：

本发明能够独立实现任意时相和区域的假定晴空条件下地表温度重建。本发明方法首先采用最邻近时相线性插值的方式对缺失值进行初步填充，再通过临年同期数据的加权融合再次缩小缺失范围，对于存在的噪声点和过拟合问题，引入SG滤波进行平滑和校正，最终生成逐日1km的无缝地表温度产品。经地表站点观测数据验证本发明提出的方法整体表现较好，将研究结果与同类产品GF(Gap-Filled)进行交叉验证，二者表现出高度相似性，但本发明提出的方法对大面积缺失区域的预测更为稳健。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例提供的一种基于MODIS遥感数据的逐日地表温度重建方法的流程图；

图2是本发明实施例中的最邻近线性插值示意图；

图3是本发明实施例中的2018-2020年DM站点MODIS晴空LST时间序列图；

图4是本发明实施例中的2020年7月15日时域重建局部过程图；

图5是本发明实施例提供的一种基于MODIS遥感数据的逐日地表温度重建装置的结构图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

参考图1，本发明实施例提供了一种基于MODIS遥感数据的逐日地表温度重建方法，该方法包括以下步骤：

S1：获取MOD11A1地表温度产品，并对其进行预处理后，筛选合格地表温度像元；

具体地，本实施例中选择Terra卫星携带的MODIS传感器所反演的MOD11A1地表温度产品，建立感兴趣区范围，将所有影像裁剪至与感兴趣区同一范围；根据受云影响程度、数据质量、发射率误差和反演误差四个指标进一步对原始数据的质量进行筛选，优选地表温度平均反演误差小于等于2K和发射率平均误差小于等于0.04这两个阈值作为筛选合格地表温度(LST)像元的依据，在其他实施例中，也可以根据实际情况和需求选择其他的阈值条件。

S2：对MOD11A1地表温度产品的单日全球缺失影像进行全域缺失时相弥补，得到连续的MOD11A1地表温度产品；

具体地，对MOD11A1地表温度产品的单日全球缺失影像进行临期复制，即通过搜索时序产品中距离缺失时相最近的影像作为该时期的初始值，以满足每日数据存在性的需求。

S3：以天为单位，对MOD11A1地表温度产品所缺影像的前后影像集使用最邻近法进行线性插值处理；

使用最邻近法进行线性插值处理的具体流程如图2所示，具体步骤为：

对于缺失像元恢复的第一步，是利用相邻前后时相晴空地表温度像元的观测值使用线性函数对缺失像元进行插值，从而快速弥补缺失像元的LST值；

地表温度虽然在长时序上存在正弦谐波函数的变价规律，但短期内的波动特征还是很难用合适的函数去拟合，线性插值往往采用多个时期的已知数据进行多元线性回归的方法来求解未知量，但该方法可能造成数据过拟合问题以至于在局部预测中出现明显偏差；

结合地表温度短期变化趋向随机性的特点和GEE像元级单元运算的方式，本发明的线性插值过程逐像元进行，所引入的晴空地表温度像元都距离缺失像元在时相上是最接近的。

以预测时刻t_m为分界线，将距离该时刻间隔为n天以内的影像作为临近时相的影像集并划分为前后两个组合，其包含的时段分别为t_m-n-t_m-1和t_m+1-t_m+n。对于预测时刻前的影像集采用末位压缩的方法，将该时段每个位置观测时间最新的晴空地表温度像元整合为一张影像，而预测时刻后的影像集则相反。对于压缩后的前后两张影像，每一个像元都附有不同的时间信息，依据像元代表的地表温度因变量和时间自变量进行一元线性插值，计算得到预测时刻像元的地表温度，其表达式为：

其中，LST_pre为预测时刻t_m的像元值，LST_up和LST_down分别为地表温度波段前后压缩的像元值，t_up和t_down为像元所包含的时间信息，b为一元线性方程的截距。

S4：结合邻年同期数据，利用地表温度具有正弦谐波函数的变化规律对重建像元进行加权融合，从而减少线性插值过程中的误差；

由于地物接收太阳辐射强度的周期性变化特征，所代表的地表温度整体呈正弦曲线的变化规律。如图3所示，以黑河流域地表过程综合观测网DM站点为例，其下垫面的玉米地在2018-2020年的LST移动平均拟合曲线不仅短期内具有良好的线性特征，并且在多年同期LST的变化区间也高度吻合。S4的具体步骤为：

对于不同年份的LST，采用反时间距离权重对预测时刻前后多个同期数据进行加权融合，其表达式为：

其中，LST_fus为加权融合后的地表温度，LST_y和ω_y代表预测时刻t_y的地表温度和权重，y为年份，LST_i/j和ω_i/j代表预测时刻t_i/j的地表温度和权重，k为距离预测年份的最大年间隔；

权重约束条件为：

S5：获取MCD12Q1地表分类产品，计算同类地物的地表温度均值，进一步填充与之同类地物的地表温度缺失像元；

具体地，S5包括：

按照MCD12Q1土地覆盖利用产品的IGBP分类标准，将研究区地物分为17类，每类地物的地表温度均值计算方法为：

其中，LST_mean是某一类地表温度均值，n为该地类的像元总数，LST_i为该地类第i个像元的地表温度值，使用LST_mean对缺失值进行进一步填充。

S6：使用SG滤波去除重建过程中的噪声以保证时序的稳定性，最终得到逐日地表温度重建结果；

步骤S6具体包括：

首先确定一个半径为m的移动窗口，将窗口内的所有数据纳入一个集合，窗口内所包含的时序数据时间节点为(-m，-m+1，…，0，1，…，m)，采用k-1次多项式对窗口内2m+1个时间点进行拟合，其多项式表达为：

y＝a₀+a₁x+a₂x²+…+a_k-1x^k-1

其中，y为滤波后的像元值，a_i(i＝1,2...k-1)为多项式参数；

窗口内每一个时间点对应一个多项式，共有2m+1个，当窗口内的像元数n大于多项式参数k时，求得方程方程组的参数，其表达式为：

用矩阵形式表示为：

Y_(2m+1)×1＝X_(2m+1)×k·A_K×l+E_(2m+1)×1

其中，Y_(2m+1)×1为滤波后LST值，X_(2m+1)×k为参与滤波的每年LST值，A_K×1为参数矩阵，E_(2m+1)×1为随机误差向量；

拟合参数A的最小二乘解

为：

Y的模型预测值或滤波值

为：

对于MOD11A1时序地表温度产品，2m+1的窗口代表预测时刻前后m天内的时段，X代表加权融合后并进行地类均值填充的地表温度值，

代表SG滤波去除噪声后的LST值。

本发明提供的一种基于MODIS遥感数据的逐日地表温度重建方法的重建结果如图4所示，a为中国行政区内部分地区2020年7月15日的重建结果，b为图中矩形融合前的地表温度由于云遮档等一些原因其大部分面积地表温度为空值，c为通过S3提出的相邻时相线性插值处理后的地表温度影像，可以看出大部分缺失像元得到填补，d为通过S4提出的邻年同期数据加权融合后的地表温度影像，经过处理后所有缺失像元均得到填补，e为SG滤波处理后的地表温度，去除了重建过程中的噪声以保证时序的稳定性。

下面对本发明提供的一种基于MODIS遥感数据的逐日地表温度重建装置进行描述，下文描述的一种基于MODIS遥感数据的逐日地表温度重建装置与上文描述的一种基于MODIS遥感数据的逐日地表温度重建方法可相互对应参照。

参考图5，本发明提供的一种基于MODIS遥感数据的逐日地表温度重建装置，具体包括以下模块：

预处理模块01，用于获取MOD11A1地表温度产品，并对其进行预处理后，筛选合格地表温度像元；

全域缺失时相弥补模块02，用于对MOD11A1地表温度产品的单日全球缺失影像进行全域缺失时相弥补，得到连续的MOD11A1地表温度产品；

最临近插值模块03，用于以天为单位，对MOD11A1地表温度产品所缺影像的前后影像集使用最邻近法进行线性插值处理；

加权融合模块04，用于结合邻年同期数据，利用地表温度具有正弦谐波函数的变化规律对重建像元进行加权融合，从而减少线性插值过程中的误差；

缺失像元填充模块05，用于获取MCD12Q1地表分类产品，计算同类地物的LST均值，进一步填充与之同类地物的地表温度缺失像元；

SG滤波模块06，用于使用SG滤波去除重建过程中的噪声，最终得到逐日地表温度重建结果。

所述基于MODIS遥感数据的逐日地表温度重建装置的各个模块与所述基于MODIS遥感数据的逐日地表温度重建方法的各个步骤一一对应，分别用于实现上述方法的各步骤，在此不再赘述。

本发明实施例提供一种基于MODIS遥感数据的逐日地表温度重建方法及装置。该方法能够在个人终端或云平台(如：Google Earth Engine)上独立实现假定晴空条件下任意区域的每日地表温度重建。该模型首先采用最邻近时相线性插值对以云污染为首的缺失值进行初步填充，再通过临年同期数据的加权融合缩小缺失范围，对于存在的噪声点和过拟合问题，引入信号学中的SG滤波进行平滑和校正，最终生成逐日1km的无缝地表温度产品。验证结果表明，该发明获得的产品能够较好地反映假定晴空下地表温度的真实状态，更能在排除天气因素干扰后有效反映地表温度的季节性特征，在长时序分析研究中更具有可靠性。

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。

上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序，可将这些词语解释为标识。

以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (10)

1.一种基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，包括以下步骤：

S1：获取MOD11A1地表温度产品，并对其进行预处理后，筛选合格地表温度像元；

S2：对MOD11A1地表温度产品的单日全球缺失影像进行全域缺失时相弥补，得到连续的MOD11A1地表温度产品；

S3：以天为单位，对MOD11A1地表温度产品所缺影像的前后影像集使用最邻近法进行线性插值处理；

S4：结合邻年同期数据，利用地表温度具有正弦谐波函数的变化规律对重建像元进行加权融合，从而减少线性插值过程中的误差；

S5：获取MCD12Q1地表分类产品，计算同类地物的地表温度均值，进一步填充与之同类地物地表温度缺失像元；

S6：使用SG滤波去除重建过程中的噪声，最终得到逐日地表温度重建结果。

2.根据权利要求1所述的基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，步骤S1中，所述预处理包括：条件筛选、影像镶嵌、研究区裁剪，像元质量审查。

3.根据权利要求1所述的基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，步骤S1具体包括：

S1.1：建立感兴趣区范围，获取MOD11A1地表温度产品，将所有影像裁剪至与感兴趣区同一范围；

S1.2：根据受云影响程度、数据质量、发射率误差和反演误差四个指标进一步对原始MOD11A1地表温度产品的质量进行筛选，获取合格地表温度像元。

4.根据权利要求3所述的基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，步骤S1.2中，选取地表温度平均反演误差小于或等于2K和发射率平均误差小于或等于0.04这两个阈值作为筛选合格地表温度像元的依据，其中，K为温度量纲，开尔文。

5.根据权利要求1所述的基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，步骤S2中，所述全域缺失时相弥补的具体方法为：

对缺失的MOD11A1单日数据进行临近复刻，即通过搜索时序产品中距离缺失时相最近的影像作为该时期的初始值。

6.根据权利要求1所述的基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，步骤S3中，所述使用最邻近法进行线性插值处理，具体包括：

利用相邻前后时相晴空地表温度像元的观测值使用线性函数对缺失像元进行插值，从而快速弥补缺失像元的地表温度值；

结合地表温度短期变化趋向随机性的特点和GEE像元级单元运算的方式，线性插值过程逐像元进行，所引入的晴空地表温度像元在时相上都最接近缺失像元所对应的时刻；

以预测时刻t_m为分界线，将距离该时刻间隔为n天以内的影像作为临近时相的影像集并划分为前后两个组合，其包含的时段分别为t_m-n-t_m-1和t_m+1-t_m+n；

对于预测时刻前的影像集采用末位压缩的方法，将该时段每个位置观测时间最新的晴空地表温度像元整合为一张影像；

对于预测时刻后的影像集采用首位压缩的方法，将该时段每个位置观测时间最新的晴空地表温度像元整合为一张影像；

对于压缩后的前后两张影像，每一个像元都附有不同的时间信息，依据像元代表的地表温度因变量和时间自变量进行一元线性插值，计算得到预测时刻像元的地表温度，其表达式为：

其中，LST_pre为预测时刻t_m的像元值，LST_up和LST_down分别为地表温度波段前后压缩的像元值，t_up和t_down为像元前后所对应的时间，b为一元线性方程的截距。

7.根据权利要求1所述的基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，步骤S4中，所述结合邻年同期数据，利用地表温度具有正弦谐波函数的变化规律对重建像元进行加权融合，具体包括：

对于不同年份的地表温度数据，采用反时间距离权重对预测时刻前后多个同期数据进行加权融合，其表达式为：

其中，LST_fus为加权融合后的地表温度，LST_y和ω_y代表预测时刻t_y的地表温度和权重，y为年份，LST_i/j和ω_i/j代表预测时刻t_i/j的地表温度和权重，k为距离预测年份的最大年间隔；

权重约束条件为：

8.根据权利要求1所述的基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，步骤S5中，获取MCD12Q1地表分类产品，计算同类地物的地表温度均值，进一步填充与之同类地物的LST缺失像元，具体包括：

按照MCD12Q1土地覆盖利用产品的IGBP分类标准，将研究区地物分为17类，每类地物的地表温度均值计算方法为：

其中，LST_mean是某一类地表温度均值，n为该地类的像元总数，LST_i为该地类第i个像元的地表温度值；

使用LST_mean进一步填充与之同类地物的地表温度缺失像元。

9.根据权利要求1所述的基于MODIS遥感数据的逐日地表温度重建方法，其特征在于，步骤S6中，所述使用SG滤波去除重建过程中的噪声，具体包括：

首先确定一个半径为m的移动窗口，将窗口内的所有数据纳入一个集合，窗口内所包含的时序数据时间节点为(-m,-m+1,…,0,1,…,m)，采用k-1次多项式对窗口内2m+1个时间点进行拟合，其多项式表达为：

y＝a₀+a₁x+a₂x²+…+a_k-1x^k-1

其中，y为滤波后的像元值，a_i(i＝1,2…k-1)为多项式参数，x为窗口内2m+1个时间点所对应的地表温度像元值；

窗口内每一个时间点对应一个多项式，共有2m+1个，当窗口内的像元数n大于多项式参数k时，求得方程方程组的参数，其表达式为：

用矩阵形式表示为：

Y_(2m+1)×1＝X_(2m+1)×k·A_K×1+E_(2m+1)×1

其中，Y_(2m+1)×1为滤波后地表温度像元值，X_(2m+1)×k为参与滤波的每年地表温度像元值，A_K×1为参数矩阵，E_(2m+1)×1为随机误差向量；

拟合参数A的最小二乘解

为：

Y的模型预测值或滤波值

为：

对于MOD11A1时序地表温度产品，2m+1的窗口代表预测时刻前后m天内的时段，X代表加权融合后并进行地类均值填充的地表温度值，

代表SG滤波去除噪声后的地表温度像元值。

10.一种基于MODIS遥感数据的逐日地表温度重建装置，其特征在于，包括以下模块：

预处理模块，用于获取MOD11A1地表温度产品，并对其进行预处理后，筛选合格地表温度像元；

全域缺失时相弥补模块，用于对MOD11A1地表温度产品的单日全球缺失影像进行全域缺失时相弥补，得到连续的MOD11A1地表温度产品；

最临近插值模块，用于以天为单位，对MOD11A1地表温度产品所缺影像的前后影像集使用最邻近法进行插值处理；

加权融合模块，用于结合邻年同期数据，利用地表温度具有正弦谐波函数的变化规律对重建像元进行加权融合，从而减少线性插值过程中的误差；

缺失像元填充模块，用于获取MCD12Q1地表分类产品，计算同类地物的地表温度均值，进一步填充与之同类地物的地表温度缺失像元；

SG滤波模块，用于使用SG滤波去除重建过程中的噪声，最终得到逐日地表温度重建结果。

Images