CN115905447B - 基于广义三角帽不确定度模型的地表温度产品融合方法 - Google Patents

Abstract

本发明公开了一种基于广义三角帽不确定度模型的地表温度产品融合方法，包括以下步骤：步骤1：获取地表温度产品；步骤2：构建地表温度产品的时间序列数据集；步骤3：逐像元估算地表温度产品不确定度；步骤4：计算地表温度产品融合权重；本发明的有益效果是：在无需任何先验知识的情况下，利用广义三角帽方法逐像元估算每种月均地表温度产品的不确定度，计算得到最优的融合权重，通过融合多种地表温度产品提高地表温度产品精度。相比基于算术平均的简单融合方法，本发明充分地集成每种月均地表温度产品各自的优势，减少随机误差对地表温度产品融合的影响。

Description

基于广义三角帽不确定度模型的地表温度产品融合方法

技术领域

本发明涉及一种基于广义三角帽不确定度估算模型的地表温度产品融合方法，属于定量遥感技术领域。

背景技术

地表温度是表征地表过程变化的一个重要的特征物理量，是地表与大气能量交换的直接驱动因子，广泛地用于地表能量平衡、气候变化和资源环境监测等研究领域。准确地估算地表温度对于理解区域气候特征的变化至关重要，并且在农业、林业、生态、全球变化等领域中发挥着重要作用。国内外学者利用不同的地表温度反演算法生产了多种地表温度产品。由于数据源和反演算法的不同，导致每种地表温度产品各有优劣。为了充分利用每种地表温度产品的优势，对多种地表温度产品进行融合是提高地表温度产品精度的有效手段。

现有地表温度产品融合方法包括：(1)基于算术平均的简单融合方法。该方法没有考虑每种地表温度产品自身的不确定性，以算术平均值的形式进行等权重融合，很难充分地利用每种地表温度产品的有用信息。(2)基于卡尔曼滤波数据同化的融合方法。该方法考虑了产品自身的不确定性，然而其结果的最优求解需要产品不确定度的先验知识作为输入。该方法通常依赖于对同化观测值和模型参数的误差的主观统计描述，因此，利用模拟数据和观测数据的信息赋予的权重具有主观性，不能反映独立数据集的最优组合。产品的不确定度是产品融合的关键因素，权衡了产品的相对精度。现有技术或缺乏对产品自身不确定性的考虑，或需要人为的先验知识作为输入，缺乏客观性。综上所述，上述方法都不能最大限度地利用多种地表温度产品各自的优势进行地表温度产品融合。因此，现有技术存在缺陷，需要改进。

发明内容

本发明的目的是提出一种基于广义三角帽不确定度估算模型的地表温度产品融合方法。在无需任何先验知识且考虑地表温度产品自身不确定度的情况下，通过融合多种地表温度产品来提高地表温度产品精度。

本发明的技术方案如下：

一种基于广义三角帽不确定度模型的地表温度产品融合方法，包括以下步骤：

步骤1：获取地表温度产品

下载2003-2019年的四种地表温度产品；

步骤2：构建地表温度产品的时间序列数据集

所述四种地表温度产品具有不同的时间分辨率和空间分辨率以及不同的投影坐标系，将这四种地表温度产品统一到相同的时间分辨率和空间分辨率，以及相同的投影坐标系；

步骤3：逐像元估算地表温度产品不确定度

利用广义三角帽方法逐像元估算每种月均地表温度产品的不确定度；

步骤4：计算地表温度产品融合权重

假设每种月均地表温度产品的时间序列数据集呈正态分布，其概率密度函数表示为：

式中，P(·)为概率密度函数，LST_j为第j种月均地表温度产品的时间序列数据集，LST_t为月均地表温度产品的真值，ε_i和σ_i分别为第i种月均地表温度产品时间序列数据集的误差(零均值的白噪声)和标准差，L(·)为似然函数；

同样地，第j种月均地表温度产品时间序列数据集的概率密度函数表示为：

式中，LST_j为第j种月均地表温度产品的时间序列数据集，ε_j和σ_j分别为第j种月均地表温度产品时间序列数据集的误差(零均值的白噪声)和标准差；

LST_t的最大似然值表示为其联合概率分布的最大值：

对公式(14)左右两边取对数，可得：

通过公式(15)左右两边对LST_t求偏导，可得：

令公式(16)等于0，可得：

将公式(17)表示为LST_t＝w_iLST_i+w_jLST_j，每种月均地表温度产品的融合权重表示为：

式中，w_i和w_j分别为第i和j种月均地表温度产品的融合权重；

对于N种月均地表温度产品，每种月均地表温度产品的融合权重表示为：

式中，N为月均地表温度产品的个数，w_k为第k种月均地表温度产品的融合权重；

利用每种月均地表温度产品的融合权重，逐像元计算融合后的月均地表温度产品：

式中，LST_m为融合后的月均地表温度产品，w_k为第k种月均地表温度产品的融合权值，LST_k为第k种月均地表温度产品。

所述的地表温度产品融合方法，步骤1中，四种地表温度产品为：(1)国家生态科学数据中心提供的全球日均地表温度产品(NESDC)，空间分辨率为1km，下载网址为http://www.nesdc.org.cn/sdo/detail？id＝60f4e35e7e28173cf0c8a771；(2)Zenodo数据共享平台提供的全球月均地表温度产品(Zenodo),空间分辨率为1km，下载网址为https://zenodo.org/record/6618442；(3)国家青藏高原科学数据中心提供的全球月均地表温度产品(TPDC)，空间分辨率为0.05°，下载网址为http://data.tpdc.ac.cn/zh-hans/data/212a9901-8a1f-43ef-bab6-c7d96e8c567a；(4)欧洲中期天气预报中心(ECMWF)提供的全球逐时月均地表温度产品(ERA5-Land)，空间分辨率为0.1°，下载网址为https://cds.climate.copernicus.eu/cdsapp#！/dataset/reanalysis-era5-land-monthly-means。

所述的地表温度产品融合方法，步骤2中，首先，将NESDC日均地表温度产品进行逐月平均，得到月均地表温度产品，空间分辨率由1km重采样为0.05°；其次，将Zenodo月均地表温度产品的空间分辨率由1km重采样为0.05°；然后，将ERA5-Land逐时月均地表温度产品进行逐时平均，得到月均地表温度产品，空间分辨率由0.1°重采样为0.05°；最后，将四种月均地表温度产品重新投影到统一的WGS84坐标系，完成空间匹配。将空间匹配后的四种月均地表温度产品按照月份排列进行时间匹配，构建四种月均地表温度产品的时间序列数据集。

所述的地表温度产品融合方法，步骤3中，对于N种月均地表温度产品，每种产品的时间序列数据集表示为：

LST_i＝LST_t+ε_i i＝1,2,......,N (1)

式中，N为月均地表温度产品的个数，LST_i为第i种月均地表温度产品的时间序列数据集，LST_t为月均地表温度产品的真值，ε_i为第i种月均地表温度产品时间序列数据集的误差；

任意选取一种月均地表温度产品作为参考，其他月均地表温度产品与参考月均地表温度产品的差值表示为：

y_i＝LST_i-LST_N＝ε_i-ε_N i＝1,2,...,N-1 (2)

式中，LST_N为任意选取的参考月均地表温度产品，ε_N为参考月均地表温度产品时间序列数据集的误差；

Y是由N-1个差值序列组成的矩阵：

式中，M为每种月均地表温度产品时间序列数据集的数据个数；

差值序列矩阵Y的协方差矩阵表示为：

式中，S为协方差矩阵，cov(·)为协方差算子，s_ij＝s_ji(i,j＝1,2,…,N-1)为y_i与y_j之间的方差(i＝j)或者协方差(i≠j)；

引入一个N×N协方差矩阵R，R为对称矩阵：

式中，σ_ij＝σ_ji(i,j＝1,…,N)为ε_i与ε_j之间的方差(i＝j)或者协方差(i≠j)；矩阵中的对角线元素σ_ii(i＝1,…,N)是待求解的未知数；

R与S的关系表示为：

式中，R_s为R的子矩阵，矩阵维数为(N-1)×(N-1)；σ_s＝[σ_1N σ_2N K σ_(N-1)N]^T是维数为(N-1)的向量；σ_NN是ε_N与ε_N之间的方差；J＝[I -U]，I是维数为(N-1)×(N-1)的单位矩阵，U＝[1 1 K 1]^T是维数为(N-1)的向量；

由于有N×(N+1)/2个未知数(R中不同元素的个数)，而只有N×(N-1)/2个方程(S中不同元素的个数)，使得公式(6)不能得到确定解；剩余的N个自由参数(σ_s，σ_NN)通过定义一个目标函数，使其最小化来确定这N个自由参数的唯一解：

式中，

为了满足det(R)＞0的条件，目标函数(公式(7))的约束条件表示为：

满足约束条件的迭代初值设置为：

通过使目标函数公式(7)在约束条件公式(8)下最小化，得到N个自由参数(σ_s，σ_NN)的唯一解；当这N个自由参数确定后，其他未知数可以通过公式(10)求解得到：

逐像元计算每种月均地表温度产品的不确定度：

式中，σ_i为第i种月均地表温度产品时间序列数据集的标准差。

本发明的有益效果是：在无需任何先验知识的情况下，利用广义三角帽方法逐像元估算每种月均地表温度产品的不确定度，计算得到最优的融合权重，通过融合多种地表温度产品提高地表温度产品精度。相比基于算术平均的简单融合方法，本发明充分地集成每种月均地表温度产品各自的优势，减少随机误差对地表温度产品融合的影响。

附图说明

图1基于广义三角帽不确定度估算模型的地表温度产品融合方法流程图；

图2利用地面实测地表温度评估四种月均地表温度产品和融合后的月均地表温度产品的精度。(a)NESDC月均地表温度产品，(b)TPDC月均地表温度产品，(c)ERA5-Land月均地表温度产品，(d)Zenodo月均地表温度产品，(e)融合后的月均地表温度产品；

具体实施方式

以下结合具体实施例，对本发明进行详细说明。

步骤1：获取地表温度产品

下载2003-2019年的四种地表温度产品：(1)国家生态科学数据中心提供的全球日均地表温度产品(NESDC)，空间分辨率为1km，下载网址为http://www.nesdc.org.cn/sdo/detail？id＝60f4e35e7e28173cf0c8a771；(2)Zenodo数据共享平台提供的全球月均地表温度产品(Zenodo),空间分辨率为1km，下载网址为https://zenodo.org/record/6618442；(3)国家青藏高原科学数据中心提供的全球月均地表温度产品(TPDC)，空间分辨率为0.05°，下载网址为http://data.tpdc.ac.cn/zh-hans/data/212a9901-8a1f-43ef-bab6-c7d96e8c567a；(4)欧洲中期天气预报中心(ECMWF)提供的全球逐时月均地表温度产品(ERA5-Land)，空间分辨率为0.1°，下载网址为https://cds.climate.copernicus.eu/cdsapp#！/dataset/reanalysis-era5-land-monthly-means。

步骤2：构建地表温度产品的时间序列数据集

上述四种地表温度产品具有不同的时间分辨率和空间分辨率以及不同的投影坐标系，需要将这四种地表温度产品统一到相同的时间分辨率和空间分辨率，以及相同的投影坐标系。首先，将NESDC日均地表温度产品进行逐月平均，得到月均地表温度产品，空间分辨率由1km重采样为0.05°；其次，将Zenodo月均地表温度产品的空间分辨率由1km重采样为0.05°；然后，将ERA5-Land逐时月均地表温度产品进行逐时平均，得到月均地表温度产品，空间分辨率由0.1°重采样为0.05°；最后，将四种月均地表温度产品重新投影到统一的WGS84坐标系，完成空间匹配。将空间匹配后的四种月均地表温度产品按照月份排列进行时间匹配，构建四种月均地表温度产品的时间序列数据集。

步骤3：逐像元估算地表温度产品不确定度

利用广义三角帽方法逐像元估算每种月均地表温度产品的不确定度。对于N种月均地表温度产品，每种产品的时间序列数据集表示为：

LST_i＝LST_t+ε_i i＝1,2,......,N (1)

式中，N为月均地表温度产品的个数，LST_i为第i种月均地表温度产品的时间序列数据集，LST_t为月均地表温度产品的真值，ε_i为第i种月均地表温度产品时间序列数据集的误差(零均值的白噪声)。

由于无法获得月均地表温度产品的真值，因此任意选取一种月均地表温度产品作为参考，其他月均地表温度产品与参考月均地表温度产品的差值表示为：

y_i＝LST_i-LST_N＝ε_i-ε_N i＝1,2,...,N-1 (2)

式中，LST_N为任意选取的参考月均地表温度产品，ε_N为参考月均地表温度产品时间序列数据集的误差(零均值的白噪声)。需要说明的是，由于广义三角帽方法在理论上对参考月均地表温度产品的选择不敏感，所以每种月均地表温度产品的不确定度不会因参考月均地表温度产品的选取的不同而改变。

Y是由N-1个差值序列组成的矩阵：

式中，M为每种月均地表温度产品时间序列数据集的数据个数。

差值序列矩阵Y的协方差矩阵表示为：

式中，S为协方差矩阵，cov(·)为协方差算子，s_ij＝s_ji(i,j＝1,2,…,N-1)为y_i与y_j之间的方差(i＝j)或者协方差(i≠j)。

引入一个N×N协方差矩阵R(R为对称矩阵)：

式中，σ_ij＝σ_ji(i,j＝1,…,N)为ε_i与ε_j之间的方差(i＝j)或者协方差(i≠j)；矩阵中的对角线元素σ_ii(i＝1,…,N)是待求解的未知数。

R与S的关系表示为：

式中，R_s为R的子矩阵，矩阵维数为(N-1)×(N-1)；σ_s＝[σ_1N σ_2N K σ_(N-1)N]^T是维数为(N-1)的向量；σ_NN是ε_N与ε_N之间的方差；J＝[I -U]，I是维数为(N-1)×(N-1)的单位矩阵，U＝[1 1 K 1]^T是维数为(N-1)的向量。

由于有N×(N+1)/2个未知数(R中不同元素的个数)，而只有N×(N-1)/2个方程(S中不同元素的个数)，使得公式(6)不能得到确定解。剩余的N个自由参数(σ_s，σ_NN)通过定义一个目标函数，使其最小化来确定这N个自由参数的唯一解：

式中，

为了满足det(R)＞0的条件，目标函数(公式(7))的约束条件表示为：

满足约束条件的迭代初值设置为：

通过使目标函数(公式(7))在约束条件(公式(8))下最小化，得到N个自由参数(σ_s，σ_NN)的唯一解。当这N个自由参数确定后，其他未知数可以通过公式(10)求解得到：

逐像元计算每种月均地表温度产品的不确定度：

式中，σ_i为第i种月均地表温度产品时间序列数据集的标准差。

步骤4：计算地表温度产品融合权重

假设每种月均地表温度产品的时间序列数据集呈正态分布，其概率密度函数表示为：

式中，P(·)为概率密度函数，LST_j为第j种月均地表温度产品的时间序列数据集，LST_t为月均地表温度产品的真值，ε_i和σ_i分别为第i种月均地表温度产品时间序列数据集的误差(零均值的白噪声)和标准差，L(·)为似然函数。

同样地，第j种月均地表温度产品时间序列数据集的概率密度函数表示为：

式中，LST_j为第j种月均地表温度产品的时间序列数据集，ε_j和σ_j分别为第j种月均地表温度产品时间序列数据集的误差(零均值的白噪声)和标准差。

LST_t的最大似然值表示为其联合概率分布的最大值：

对公式(14)左右两边取对数，可得：

通过公式(15)左右两边对LST_t求偏导，可得：

令公式(16)等于0，可得：

将公式(17)表示为LST_t＝w_iLST_i+w_jLST_j，每种月均地表温度产品的融合权重表示为：

式中，w_i和w_j分别为第i和j种月均地表温度产品的融合权重。

对于N种月均地表温度产品，每种月均地表温度产品的融合权重表示为：

式中，N为月均地表温度产品的个数，w_k为第k种月均地表温度产品的融合权重。

利用每种月均地表温度产品的融合权重，逐像元计算融合后的月均地表温度产品：

式中，LST_m为融合后的月均地表温度产品，w_k为第k种月均地表温度产品的融合权值，LST_k为第k种月均地表温度产品。

图2显示了利用地面实测地表温度评估四种月均地表温度产品和融合后的月均地表温度产品的精度。相比单一的月均地表温度产品，融合后的地表温度产品的精度提高了0.2-1.0K。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (4)

1.一种基于广义三角帽不确定度模型的地表温度产品融合方法，其特征在于，包括以下步骤：

步骤1：获取地表温度产品

下载四种地表温度产品；

步骤2：构建地表温度产品的时间序列数据集

所述四种地表温度产品具有不同的时间分辨率和空间分辨率以及不同的投影坐标系，将这四种地表温度产品统一到相同的时间分辨率和空间分辨率，以及相同的投影坐标系；

步骤3：逐像元估算地表温度产品不确定度

利用广义三角帽方法逐像元估算每种月均地表温度产品的不确定度；

步骤4：计算地表温度产品融合权重

假设每种月均地表温度产品的时间序列数据集呈正态分布，其概率密度函数表示为：

式中，P(·)为概率密度函数，LST_j为第j种月均地表温度产品的时间序列数据集，LST_t为月均地表温度产品的真值，ε_i和σ_i分别为第i种月均地表温度产品时间序列数据集的误差(零均值的白噪声)和标准差，L(·)为似然函数；

同样地，第j种月均地表温度产品时间序列数据集的概率密度函数表示为：

式中，LST_j为第j种月均地表温度产品的时间序列数据集，ε_j和σ_j分别为第j种月均地表温度产品时间序列数据集的误差(零均值的白噪声)和标准差；

LST_t的最大似然值表示为其联合概率分布的最大值：

对公式(14)左右两边取对数，可得：

通过公式(15)左右两边对LST_t求偏导，可得：

令公式(16)等于0，可得：

将公式(17)表示为LST_t＝w_iLST_i+w_jLST_j，每种月均地表温度产品的融合权重表示为：

式中，w_i和w_j分别为第i和j种月均地表温度产品的融合权重；

对于N种月均地表温度产品，每种月均地表温度产品的融合权重表示为：

式中，N为月均地表温度产品的个数，w_k为第k种月均地表温度产品的融合权重；

利用每种月均地表温度产品的融合权重，逐像元计算融合后的月均地表温度产品：

式中，LST_m为融合后的月均地表温度产品，w_k为第k种月均地表温度产品的融合权值，LST_k为第k种月均地表温度产品。

2.根据权利要求1所述的地表温度产品融合方法，其特征在于，步骤1中，下载2003-2019年的四种地表温度产品，四种地表温度产品为：(1)国家生态科学数据中心提供的全球日均地表温度产品(NESDC)，空间分辨率为1km，下载网址为http://www.nesdc.org.cn/sdo/detail？id＝60f4e35e7e28173cf0c8a771；(2)Zenodo数据共享平台提供的全球月均地表温度产品(Zenodo),空间分辨率为1km，下载网址为https://zenodo.org/record/6618442；(3)国家青藏高原科学数据中心提供的全球月均地表温度产品(TPDC)，空间分辨率为0.05°，下载网址为http://data.tpdc.ac.cn/zh-hans/data/212a9901-8a1f-43ef-bab6-c7d96e8c567a；(4)欧洲中期天气预报中心(ECMWF)提供的全球逐时月均地表温度产品(ERA5-Land)，空间分辨率为0.1°，下载网址为https://cds.climate.copernicus.eu/cdsapp#！/dataset/reanalysis-era5-land-monthly-means。

3.根据权利要求1所述的地表温度产品融合方法，其特征在于，步骤2中，首先，将NESDC日均地表温度产品进行逐月平均，得到月均地表温度产品，空间分辨率由1km重采样为0.05°；其次，将Zenodo月均地表温度产品的空间分辨率由1km重采样为0.05°；然后，将ERA5-Land逐时月均地表温度产品进行逐时平均，得到月均地表温度产品，空间分辨率由0.1°重采样为0.05°；最后，将四种月均地表温度产品重新投影到统一的WGS84坐标系，完成空间匹配；将空间匹配后的四种月均地表温度产品按照月份排列进行时间匹配，构建四种月均地表温度产品的时间序列数据集。

4.根据权利要求1所述的地表温度产品融合方法，其特征在于，步骤3中，对于N种月均地表温度产品，每种产品的时间序列数据集表示为：

LST_i＝LST_t+ε_ii＝1,2,......,N (1)

式中，N为月均地表温度产品的个数，LST_i为第i种月均地表温度产品的时间序列数据集，LST_t为月均地表温度产品的真值，ε_i为第i种月均地表温度产品时间序列数据集的误差；

任意选取一种月均地表温度产品作为参考，其他月均地表温度产品与参考月均地表温度产品的差值表示为：

y_i＝LST_i-LST_N＝ε_i-ε_Ni＝1,2,...,N-1 (2)

式中，LST_N为任意选取的参考月均地表温度产品，ε_N为参考月均地表温度产品时间序列数据集的误差；

Y是由N-1个差值序列组成的矩阵：

式中，M为每种月均地表温度产品时间序列数据集的数据个数；

差值序列矩阵Y的协方差矩阵表示为：

式中，S为协方差矩阵，cov(·)为协方差算子，s_ij＝s_ji(i,j＝1,,2,…,N-1)为y_i与y_j之间的方差(i＝j)或者协方差(i≠j)；

引入一个N×N协方差矩阵R，R为对称矩阵：

式中，σ_ij＝σ_ji(i,j＝1,…,N)为ε_i与ε_j之间的方差(i＝j)或者协方差(i≠j)；矩阵中的对角线元素σ_ii(i＝1,…,N)是待求解的未知数；

R与S的关系表示为：

式中，R_s为R的子矩阵，矩阵维数为(N-1)×(N-1)；σ_s＝[σ_1N σ_2N K σ_(N-1)N]^T是维数为(N-1)的向量；σ_NN是ε_N与ε_N之间的方差；J＝[I -U]，I是维数为(N-1)×(N-1)的单位矩阵，U＝[11 K 1]^T是维数为(N-1)的向量；

由于有N×(N+1)/2个未知数(R中不同元素的个数)，而只有N×(N-1)/2个方程(S中不同元素的个数)，使得公式(6)不能得到确定解；剩余的N个自由参数(σ_s，σ_NN)通过定义一个目标函数，使其最小化来确定这N个自由参数的唯一解：

式中，

为了满足det(R)＞0的条件，目标函数(公式(7))的约束条件表示为：

满足约束条件的迭代初值设置为：

通过使目标函数公式(7)在约束条件公式(8)下最小化，得到N个自由参数(σ_s，σ_NN)的唯一解；当这N个自由参数确定后，其他未知数可以通过公式(10)求解得到：

逐像元计算每种月均地表温度产品的不确定度：

式中，σ_i为第i种月均地表温度产品时间序列数据集的标准差。

Images