CN109635309A - 一种地表温度空间降尺度方法 - Google Patents

Abstract

本发明公开了一种地表温度空间降尺度方法，本方法首先定量分析地表温度与地表参数包括不透水面盖度、植被盖度、土壤盖度、NDVI、NDBI、MNDWI、DEM和建筑密度间的相关性及其空间分布差异，然后利用机器学习的随机森林算法建立低空间分辨率地表温度产品与相关地表参数的回归模型，并结合高空间分辨率的地表参数进行预测得到高空间分辨率的地表温度预测值，然后通过地统计学理论的面到点克里格插值方法对随机森林回归模型的残差进行降尺度来提高回归模型残差的空间分辨率，最后将高空间分辨率的随机森林回归模型预测地表温度与面到点克里格插值的残差相加生成高分辨率、高精度地表温度产品，弥补已有地表温度产品空间分辨率不足的缺陷。

Description

一种地表温度空间降尺度方法

技术领域

本发明涉及遥感测量技术领域，具体涉及一种基于机器学习与地统计学理论的地表温度空间降尺度方法。

背景技术

地表温度是区域和全球尺度上陆地表层系统过程的关键参数，它综合了地表与大气的相互作用以及大气和陆地之间能量交换的综合结果。地表温度能够提供地表能量平衡状态的时空变化信息，在气象、地质、水文、生态、数值预报、以及区域气候模式等研究领域得到广泛的应用，现已成为城市热岛监测与热岛效应研究的主要参量。精确的地表温度不仅有助于评估地表能量与水文平衡、热惯量和土壤湿度，而且有助于获取全球表面温度及掌握其长期的变化，对区域及全球地气系统能量平衡和生态系统的研究具有重要意义。现在有卫星遥感反演地表温度产品存在着时间分辨率和空间分辨率的矛盾，单一卫星产品无法满足精细化的地表温度时空分辨监测与应用研究。ASTER和Landsat地表温度产品具有高空间分辨率，已被广泛应用于小面积监测，如城市热岛监测。但是这些地表温度产品时间分辨率低，并不适合时间序列监测。MODIS、FY3等地表温度产品具有较高时间分辨率，已广泛应用于全球尺度或大范围面积监测。由于空间分辨率较低，MODIS、FY3等地表温度产品并不适合小面积的高精度监测。因此，在较高时间分辨率地表温度数据基础上，增强、提高其空间分辨率，是当前热红外遥感应用研究中的热点与关键问题之一。

基于地表温度的降尺度转换方法可分为热红外锐化和温度分解2类。尽管针对地表温度降尺度的方法多种多样，但其实质都是建立热红外波段信息与各类地表参数间的关系模型，并假定这种关系模型不随空间尺度发生变化。上述简单的单因子、多因子回归关系并不能完全概述不同尺度因子与地表温度间的复杂关系。在物理机制尚不清楚的前提下，使用机器学习方法构建降尺度模型是较好的选择。与人工神经网络、支持向量机等机器学习方法相比，随机森林算法具有运算量小、容纳样本数量大等优点，适用于遥感降尺度研究。

由于现有星载热红外传感器技术的限制，现有的FY3C-VIRR每日地表温度产品空间分辨率较低(1km)，无法满足精细化的地表温度监测与应用研究，限制了高空间分辨率的地表温度信息在城市热环境研究、城市能效监测评估、建筑热能消耗检测、人工材料检测等研究中的应用。鉴于此，诸多学者利用高空间分辨率的地表覆盖信息进行地表温度的空间降尺度研究。这些地表温度降尺度方法大多是建立热红外波段信息与各类地表参数间的关系模型，并假定这种关系模型不随空间尺度发生变化。然而，这些简单的单因子、多因子回归关系并不能完全概述不同尺度因子与地表温度间的复杂关系。大范围或者地表覆盖复杂区域的地表温度空间降尺度具有较大的不确定性，有待进一步研究。

发明内容

针对现有技术的不足，本发明的目的旨在提供一种基于机器学习与地统计学理论的地表温度空间降尺度方法，以弥补已有地表温度产品空间分辨率不足的缺陷。

为实现上述目的，本发明采用如下技术方案：

一种地表温度空间降尺度方法，包括

获取地理空间影像数据，并对该影像数据进行预处理，以提取出高空间分辨率地表参数；

对高空间分辨率地表参数进行升尺度，以得到低空间分辨率地表参数；

选择与地表温度显著相关的低空间分辨率地表参数作为低空间分辨率地表参数指示因子；

选择与地表温度显著相关的高空间分辨率地表参数作为高空间分辨率地表参数指示因子；

利用机器学习的随机森林算法来建立低空间分辨率的地表温度与低空间分辨率地表参数指示因子的回归模型，以得到低空间分辨率的随机森林回归模型及其低空间分辨率的随机森林回归模型残差；

将低空间分辨率的随机森林回归模型应用于高空间分辨率地表参数指示因子，以获得高空间分辨率的随机森林回归模型的地表温度预测值；

对低空间分辨率的随机森林回归模型模型残差进行降尺度，以获得高空间分辨率的随机森林回归模型的地表温度残差估计值；

将所获得的高空间分辨率的随机森林回归模型地表温度预测值和高空间分辨率的随机森林回归模型的地表温度残差估计值相加，得到高空间分辨率的地表温度值。

所述高空间分辨率地表参数和低空间地表参数均包括不透水面盖度、植被盖度、土壤盖度、NDVI、NDBI、MNDWI、DEM和建筑密度。

提取出高空间分辨率地表参数的过程为：

分别利用地理空间影像数据的绿波段、红波段、近红外波段以及中红外波段的反射率计算NDVI、NDBI和MNDWI指数；

对地理空间影像数据首先利用改进的归一化水体指数去除水体，然后利用线性光谱解混模型对Landsat-8遥感影像数据的各波段反射率进行求解各端元比例；

线性光谱解混模型假设影像每个像元的反射率由该像元所有地物端元的反射率及其所占面积比例为权重系数的线性组合，其表达式为：

式中i＝1,2…,M，M为光谱波段数；n为地物端元数目；R_i是波段i的反射率，f_k是端元k在像元中所占的面积比例，R_ik波段i的第k个端元的光谱反射率，ER_i是波段i的残差；选择土壤、植被、高反照率和低反照率四个端元作为解混对象，采用最小二乘法求解各端元所占的比例，且线性光谱解混的求解必须满足以下的条件：

且f_k≥0 (5)

计算出高反照率、低反照率、植被和土壤盖度；

接下来，结合NDBI和NDVI指数对高反照率、低反照率、植被和土壤盖度进行修正，得到最终的不透水面、植被和土壤盖度。

所述地理空间影像数据为Landsat 8OLI遥感影像数据。

得到最终的不透水面、植被和土壤盖度的过程为：

首先利用Otsu算法对NDBI指数影像进行最优阈值提取，得到建筑物影像和非建筑物影像；如果当前像元的NDBI指数大于最优阈值时，即当前像元为建筑物像元时，此像元的低反照率盖度划分为低反照率不透水面盖度，否则，此像元的低反照率盖度划分为低反照率透水面盖度；如果当前像元的NDVI指数大于0.2时，此像元的低反照率透水面盖度属于低反照率植被盖度，否则，此像元的低反照率透水面盖度属于土壤盖度；最后将高反照率盖度和低反照率不透水面盖度相加得到不透水面盖度，将原来的土壤盖度与低反照率土壤盖度相加得到修正的土壤盖度，将原来的植被盖度与低反照率植被盖度相加得到修正的植被盖度。

对高空间分辨率地表参数采用平均聚合的方法来进行升尺度。

对低空间分辨率的随机森林回归模型模型残差采用地统计学理论的面到点克里格插值方法来进行进行降尺度处理，具体过程为：

首先计算低空间分辨率的随机森林回归模型残差在低分辨率面尺度上的变异函数值；然后通过去卷积方法估算低空间分辨率的随机森林回归模回归模型残差点尺度变异函数；在这个基础上，计算面与面、面与点在不同尺度上的低空间分辨率的随机森林回归模型残差间的变异函数值，从而构建用于估算高空间分辨率的随机森林回归模型残差的权重线性方程组，并进行求解，获得估算的权重值；结合估算的权重值以及对应的低空间分辨率的随机森林回归模型残差进行加权平均计算，得到高空间分辨率的随机森林回归模型的地表温度残差估计值。

通过相关性分析定量分析地表温度与高空间分辨率地表参数间的相关性及其空间分布差异，选择与地表温度显著相关的高空间分辨率地表参数作为高空间分辨率地表参数指示因子。

通过相关性分析定量分析地表温度与低空间分辨率地表参数间的相关性及其空间分布差异，选择与地表温度显著相关的低空间分辨率地表参数作为高空间分辨率地表参数指示因子。

本发明的有益效果在于：

现有FY3C-VIRR每日1km地表温度产品覆盖范围广，但空间分辨率较低，无法满足精细化的地表温度监测与应用研究，为解决此问题，本方法首先定量分析地表温度与地表参数包括不透水面盖度、植被盖度、土壤盖度、NDVI、NDBI、MNDWI、DEM和建筑密度间的相关性及其空间分布差异，然后利用机器学习的随机森林算法建立FY3C-VIRR地表温度产品与相关地表参数的回归模型，并结合高空间分辨率的地表参数进行预测得到高空间分辨率的地表温度预测值，然后通过地统计学理论的面到点克里格插值方法对随机森林回归模型的残差进行降尺度来提高回归模型残差的空间分辨率，最后将高空间分辨率的随机森林回归模型预测地表温度与面到点克里格插值的残差相加生成高分辨率、高精度地表温度产品，弥补已有FY3C-VIRR地表温度产品空间分辨率不足的缺陷。

附图说明

图1为本发明实施例提供的地表温度空间降尺度方法的流程图；

图2为改进的LSMA方法提取Landsat 8影像的不透水面、植被和土壤盖度流程；

图3为空间降尺度方法示意图；

图4为随机森林回归模型残差的空间降尺度过程图；

具体实施方式

下面，结合附图以及具体实施方式，对本发明做进一步描述：

针对FY3C-VIRR地表温度产品空间分辨率不高的问题，本实施例利用随机森林与面到点克里金方法，以NDVI、NDBI、MNDWI、不透水面盖度、DEM以及建筑密度等地表参数为指示因子，构建FY3C-VIRR地表温度产品的空间降尺度模型，将1km的地表温度产品降尺度到30m分辨率，也就是说，在本实施例中，低空间分辨率指的是1km分辨率，高空间分辨率指的的是30m分辨率率。参阅图1所示，本方法具体包括：

(1)影像数据获取及预处理

首先在美国联邦地质调查局地球资源观测与科学中心(https://espa.cr.usgs.gov/)下载Landsat 8OLI影像数据，并进行大气校正、几何校正等预处理；然后在国家气象卫星中心下载1km的每日FY3C-VIRR地表温度产品，并根据研究区域范围进行裁剪；在地理空间数据云上下载空间分辨率为30m的DEM(Digital Elevation Model-数字高程模型)数据；结合Wordview2影像，利用ArcGIS软件对研究区建筑物进行矢量化处理，得到建筑物的矢量图层，并转换为空间分辨率为0.5m的建筑物栅格数据，通过聚合计算得到空间分辨率为30m的建筑物密度栅格数据；最后将上述影像数据统一投影到UTM/WGS84中。

(2)高分辨率地表参数提取

利用30m空间分辨率的Landsat 8OLI遥感影像数据提取NDVI(NormalizedDifference Vegetation Index,归一化植被指数)、NDBI(Normalized DifferenceBuilding Index，归一化建筑指数)、MNDWI(Normalized Difference Water Index，归一化水指数)、不透水面盖度、植被盖度和土壤盖度等地表参数。

分别利用Landsat 8OLI遥感影像数据的绿波段、红波段、近红外波段以及中红外波段的反射率计算NDVI、NDBI和MNDWI指数：

式中，ρ_Green、ρ_R、ρ_NIR和ρ_MIR分别为Landsat 8遥感影像的绿波段、红波段、近红外和中红外波段反射率。

下面利用改进的线性光谱解混方法提取Landsat 8遥感影像的不透水面、植被和土壤盖度数据，如图2所示。对Landsat 8OLI遥感影像首先利用改进的归一化水体指数(MNDWI)去除水体，然后利用线性光谱解混模型(LSMA：Linear SpectralMixtureAnalysis)对Landsat-8的各波段反射率进行求解各端元比例。线性光谱解混模型(LSMA)假设影像每个像元的反射率由该像元所有地物端元的反射率及其所占面积比例为权重系数的线性组合，其表达式为：

式中i＝1,2…,M，M为光谱波段数；n为地物端元数目。R_i是波段i的反射率，f_k是端元k在像元中所占的面积比例，R_ik波段i的第k个端元的光谱反射率，ER_i是波段i的残差。这里选择土壤、植被、高反照率和低反照率四个端元作为解混对象，采用最小二乘法求解各端元所占的比例，且线性光谱解混的求解必须满足以下的条件：

且f_k≥0 (5)

通过传统的线性光谱解混方法计算出高反照率、低反照率、植被和土壤盖度。

接下来，结合NDBI和NDVI指数对高反照率、低反照率、植被和土壤盖度进行修正，得到最终的不透水面、植被和土壤盖度，以为后续的精确计算提供准确的参数。具体地：

首先利用Otsu算法对NDBI指数影像进行最优阈值提取，得到建筑物影像和非建筑物影像；如果当前像元的NDBI指数大于最优阈值时，即当前像元为建筑物像元时，此像元的低反照率盖度划分为低反照率不透水面盖度，否则，此像元的低反照率盖度划分为低反照率透水面盖度；如果当前像元的NDVI指数大于0.2时，此像元的低反照率透水面盖度属于低反照率植被盖度，否则，此像元的低反照率透水面盖度属于土壤盖度；最后将高反照率盖度和低反照率不透水面盖度相加得到不透水面盖度，将原来的土壤盖度与低反照率土壤盖度相加得到修正的土壤盖度，将原来的植被盖度与低反照率植被盖度相加得到修正的植被盖度。

(3)高分辨率地表参数升尺度处理

为了使得地表参数的空间分辨率与FY3C-VIRR地表温度产品的空间分辨率保持一致，将30m的不透水面盖度、植被盖度、土壤盖度、NDVI、NDBI、MNDWI、DEM以及建筑密度等地表参数，利用平均聚合的方法进行升尺度，得到1km的不透水面盖度、植被盖度、土壤盖度、NDVI、NDBI、MNDWI、DEM以及建筑密度。

(4)地表温度与地表参数的相关性分析

利用Pearson相关性分析法，定量探究地表温度与地表参数的相关性及其空间分布差异，最终选择与地表温度显著相关的地表参数作为地表参数指示因子。

(5)基于随机森林算法与地统计学理论的地表温度空间降尺度

利用机器学习算法和地统计学理论，以步骤(4)选择出的地表参数作为指示因子，构建FY3C-VIRR地表温度产品的空间降尺度模型，将1km的FY3C-VIRR地表温度产品降尺度到30m分辨率，其过程如图3所示。

首先，利用随机森林回归模型构建低空间分辨率(1km)下的FY3C-VIRR地表温度产品与相关地表环境指示因子的回归模型f(X)，得到低空间分辨率下的随机森林回归模型系数以及对应的回归模型残差。在此基础上，将低分辨率下的随机森林回归模型应用于高空间分辨率(30m)的地表环境指示因子数据，如30m的NDVI、NDBI、不透水面盖度、建筑密度等，从而得到具有高空间分辨率的FY3C地表温度预测数据f(X_i)。

在获取高空间分辨率的地表温度预测数据后，接下来还需要对低空间分辨率下的随机森林回归模型的残差进行处理。这里采用地统计学的面到点克里格方法对低空间分辨率的回归模型残差进行插值处理，得到高空间分辨率的回归模型残差，具体的计算过程如图4所示。

具体的计算过程为：首先计算随机森林回归模型残差在低分辨率面尺度上的变异函数值；然后通过去卷积方法估算回归模型残差点尺度变异函数；在这个基础上，计算面与面、面与点在不同尺度上的回归模型残差间的变异函数值，从而构建用于估算高空间分辨率回归模型残差的权重线性方程组，并进行求解；结合估算的权重值以及对应的低空间分辨率的随机森林回归模型残差进行加权平均计算得到高空间分辨率的随机森林回归模型残差估计值。

面到点的克里格插值基本原理如下：取空间分辨率为30m的遥感网格中心为规则分布待插值点，采用地统计学理论的面到点克里格插值(Area-to-pointKriging,ATPK)进行随机森林回归模型的残差进行降尺度。ATPK是用已知面对未知点进行插值的空间降尺度方法，其原理是利用z(x)所在面以及周边面ν_i(i＝1,…,K)数据进行线性加权求和得到z(x)估算值：

权重λ_x的求解方程组为：

式中x为待插值高空间分辨网格点，共有K个低空间分辨率的网格v_i(i＝1,…,k)，对随机森林回归模型残差z(v_i)进行加权平均估计得到x的估计值z(x)，对应的权重为λ_x，μ_x为拉格朗日算子。ATPK随机森林回归模型残差降尺度过程中最关键的是去卷积估算点与点协方差函数C，并以此分别构造面与点、面与面协方差函数和如下式所示：

N(v_i)和N(v_j)分别是网格v_i和v_j离散后包含点的数量，s即代表网格中的离散点。低空间分辨率的随机森林回归模型残差网格均被离散成点数据，面与点、面与面协方差通过被转化为众多点对协方差的均值来计算。

在获取高空间分辨率的随机森林回归模型地表温度预测值以及高空间分辨率的随机森林回归模型的温度残差值后，将两者相加，最终得到高空间分辨率(30m)的FY3C-VIRR地表温度数据，从而实现了对低空间分辨率的FY3C-VIRR地表温度产品空间降尺度。

由此可知，本首先利用改进的线性光谱解混方法提取Landsat 8OLI影像的不透水面、植被和土壤盖度，并通过相关性分析定量分析地表温度与地表参数包括不透水面盖度、植被盖度、土壤盖度、NDVI、NDBI、MNDWI、DEM和建筑密度间的相关性及其空间分布差异，选择与地表温度显著相关的地表参数作为指示因子；然后利用机器学习的随机森林算法建立FY3C-VIRR地表温度产品与相关地表参数的回归模型，得到低空间分辨率的随机森林回归模型及其模型误差，并利用该模型，获得高空间分辨率的地表温度数据；针对随机森林回归模型误差，采用地统计学理论的面到点克里格方法对其进行降尺度来提高回归模型残差的空间分辨率，最后将高空间分辨率的随机森林回归模型预测地表温度与面到点克里格插值的残差相加生成高分辨率、高精度地表温度产品，弥补已有FY3C-VIRR地表温度产品空间分辨率不足的缺陷。本方法适用于在大范围或者地表覆盖复杂的区域，能较好地提高FY3C-VIRR地表温度产品的降尺度精度。

对本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及形变，而所有的这些改变以及形变都应该属于本发明权利要求的保护范围之内。

Claims (9)

1.一种地表温度空间降尺度方法，其特征在于，所述方法包括

获取地理空间影像数据，并对该影像数据进行预处理，以提取出高空间分辨率地表参数；

对高空间分辨率地表参数进行升尺度，以得到低空间分辨率地表参数；

选择与地表温度显著相关的低空间分辨率地表参数作为低空间分辨率地表参数指示因子；

选择与地表温度显著相关的高空间分辨率地表参数作为高空间分辨率地表参数指示因子；

利用机器学习的随机森林算法来建立低空间分辨率的地表温度与低空间分辨率地表参数指示因子的回归模型，以得到低空间分辨率的随机森林回归模型及其低空间分辨率的随机森林回归模型残差；

将低空间分辨率的随机森林回归模型应用于高空间分辨率地表参数指示因子，以获得高空间分辨率的随机森林回归模型的地表温度预测值；

对低空间分辨率的随机森林回归模型模型残差进行降尺度，以获得高空间分辨率的随机森林回归模型的地表温度残差估计值；

将所获得的高空间分辨率的随机森林回归模型地表温度预测值和高空间分辨率的随机森林回归模型的地表温度残差估计值相加，得到高空间分辨率的地表温度值。

2.如权利要求1所述的地表温度空间降尺度方法，其特征在于，所述高空间分辨率地表参数和低空间地表参数均包括不透水面盖度、植被盖度、土壤盖度、NDVI、NDBI、MNDWI、DEM和建筑密度。

3.如权利要求2所述的地表温度空间降尺度方法，其特征在于，提取出高空间分辨率地表参数的过程为：

分别利用地理空间影像数据的绿波段、红波段、近红外波段以及中红外波段的反射率计算NDVI、NDBI和MNDWI指数；

对地理空间影像数据首先利用改进的归一化水体指数去除水体，然后利用线性光谱解混模型对Landsat-8遥感影像数据的各波段反射率进行求解各端元比例；

线性光谱解混模型假设影像每个像元的反射率由该像元所有地物端元的反射率及其所占面积比例为权重系数的线性组合，其表达式为：

式中i＝1,2…,M，M为光谱波段数；n为地物端元数目；R_i是波段i的反射率，f_k是端元k在像元中所占的面积比例，R_ik波段i的第k个端元的光谱反射率，ER_i是波段i的残差；选择土壤、植被、高反照率和低反照率四个端元作为解混对象，采用最小二乘法求解各端元所占的比例，且线性光谱解混的求解必须满足以下的条件：

计算出高反照率、低反照率、植被和土壤盖度；

接下来，结合NDBI和NDVI指数对高反照率、低反照率、植被和土壤盖度进行修正，得到最终的不透水面、植被和土壤盖度。

4.如权利要求1或3所述的地表温度空间降尺度方法，其特征在于，所述地理空间影像数据为Landsat 8 OLI遥感影像数据。

5.如权利要求3所述的地表温度空间降尺度方法，其特征在于，得到最终的不透水面、植被和土壤盖度的过程为：

首先利用Otsu算法对NDBI指数影像进行最优阈值提取，得到建筑物影像和非建筑物影像；如果当前像元的NDBI指数大于最优阈值时，即当前像元为建筑物像元时，此像元的低反照率盖度划分为低反照率不透水面盖度，否则，此像元的低反照率盖度划分为低反照率透水面盖度；如果当前像元的NDVI指数大于0.2时，此像元的低反照率透水面盖度属于低反照率植被盖度，否则，此像元的低反照率透水面盖度属于土壤盖度；最后将高反照率盖度和低反照率不透水面盖度相加得到不透水面盖度，将原来的土壤盖度与低反照率土壤盖度相加得到修正的土壤盖度，将原来的植被盖度与低反照率植被盖度相加得到修正的植被盖度。

6.如权利要求1所述的地表温度空间降尺度方法，其特征在于，对高空间分辨率地表参数采用平均聚合的方法来进行升尺度。

7.如权利要求1所述的地表温度空间降尺度方法，其特征在于，对低空间分辨率的随机森林回归模型模型残差采用地统计学理论的面到点克里格插值方法来进行进行降尺度处理，具体过程为：

首先计算低空间分辨率的随机森林回归模型残差在低分辨率面尺度上的变异函数值；然后通过去卷积方法估算低空间分辨率的随机森林回归模回归模型残差点尺度变异函数；在这个基础上，计算面与面、面与点在不同尺度上的低空间分辨率的随机森林回归模型残差间的变异函数值，从而构建用于估算高空间分辨率的随机森林回归模型残差的权重线性方程组，并进行求解，获得估算的权重值；结合估算的权重值以及对应的低空间分辨率的随机森林回归模型残差进行加权平均计算，得到高空间分辨率的随机森林回归模型的地表温度残差估计值。

8.如权利要求1所述的地表温度空间降尺度方法，其特征在于，通过相关性分析定量分析地表温度与高空间分辨率地表参数间的相关性及其空间分布差异，选择与地表温度显著相关的高空间分辨率地表参数作为高空间分辨率地表参数指示因子。

9.如权利要求1或8所述的地表温度空间降尺度方法，其特征在于，通过相关性分析定量分析地表温度与低空间分辨率地表参数间的相关性及其空间分布差异，选择与地表温度显著相关的低空间分辨率地表参数作为高空间分辨率地表参数指示因子。