使用多类型预测变量的随机森林回归模型对陆地表面温度进行降尺度分析

本文利用多类型预测变量的随机森林回归模型(MTVRF)，以地表反射率、光谱指数、地形因子和土地覆盖类型等四种输入变量建立地表温度(LSTs)与其他地表参数的非线性关系。本文的主要目的是分析MTVRF模型在多变量回归中的优越性，无论在简单或复杂的底层地表，并进一步证明在一个研究区域训练的随机森林（RF）回归降尺度模型在应用于另一个区域时的稳健性。中分辨率成像光谱仪LST产品的空间分辨率被MTVRF降尺度从990米降至90米，并与其他两种降尺度方法（如基本RF模型和热锐化算法）进行了比较。通过计算降尺度和参考LST之间的平均误差、决定系数( )和均方根误差(RMSE)，MTVRF模型取得了满意的效果。在对不同土地植被的LST进行降尺度时, MTVRF也取得了令人满意的结果, 并对不同地区的训练模型进行了评估。在研究区 B 训练的 MTVRF 模型的 RMSE 为 3.13k, 而在研究区 A 训练的 RMSE 为 2.11k; 这表明在特定区域训练的 MTVRF 模型被认为对其他不同地表条件下的 LST 的降尺度具有足够的鲁棒性。

引言

地表温度(LST)被描述为区域和全球范围内陆地生态系统中大气和地表相互作用、物质循环和能量交换的最重要参数之一，在地表能量平衡和生物物理参数的求解中起着至关重要的作用[1]-[6]。分辨率较高的LST已被广泛用于水文平衡评估、全球变暖研究、城市热岛效应评估和地表蒸发量计算[7]-[12]。

热红外遥感数据通常来自热红外传感器，如热红外传感器(TIRS)、高级星载热发射和反射辐射计(ASTER)和中分辨率成像光谱仪(MODIS)，已被持续用于检索陆地表面温度。然而，遥感的LST产品在高时间分辨率和高空间分辨率之间产生了权衡[13]。例如，Landsat ETM plus在波段6的分辨率为60米，但周期为16天，不能满足同一地区的连续时间序列观测。高时间分辨率的传感器，如MODIS，每天获得两次观测，但只能提供1公里的LST。低空间分辨率的热红外图像由于每个像素内的子组件，会受到空间异质性的影响。这个问题阻碍了LST的应用。因此，需要一种有效的降尺度方法来提高LST的分辨率[14]-[16]。

在过去的几年中, 人们提出了各种降尺度算法, 将 LST 从粗略的分辨率降到更细的分辨率[13], [17]-[25], 其中基于尺度因子的降尺度方法被认为是一种流行的算法[15]-[16], [26].基于尺度因子的降尺度方法通常分为三类：基于统计回归的降尺度方法[13], [19], [27], [28], 基于调制分布的降尺度方法和基于线性谱系混合模型的降尺度方法[29]-[35]。统计回归方法，如辐射表面温度分解程序（the disaggregation procedure for radiometric surface temperature, DisTrad）方法[13]和热锐化（Thermal Sharpening, TsHARP）算法[19]，在这些方法中被普遍认为是容易操作和令人满意的准确度。早期的回归算法倾向于关注植被指数与 LST 之间的统计关系, 如归一化差异植被指数(NDVI)[13]、分位植被覆盖率[17]-[19]、土壤调整植被指数(SAVI)[24], 其依据是地表温度的空间差异主要由植被覆盖率控制[36], [37].然而, 这些算法的主要限制是, 对于有不透水表面、不毛之地和湿地的地区, 植被指数和LST之间的相关性似乎不够充分[29], [38], [39].因此, 许多其他光谱指数被用来拟合 LST 和土地参数之间的非线性关系, 包括归一化差异建筑指数 (NDBI)、修正的归一化水指数(MNDWI)[28]和归一化差异干旱指数(NDDI) [22].其他尺度因子，如地表反照率、数字评价模型和土地覆盖类型，也不断被引入以改进统计回归算法[27]，[41]。

除此方法外, 近期的研究还涉及到 LST 与所述参数之间的线性或非线性回归算法问题[42]-[45], 例如, Duan 和 Li 提出了一种独创的地理加权回归降尺度模型, 并取得了良好的效果[45].基于贝叶斯的模型[42]、人工神经网络(ANN)[46]-[48]、支持向量机(SVM)[49]、结合全局窗和移动窗的回归树以及随机森林(RF)[25]、[28]、[50]、[51]等机器学习算法在拟合LST与其他变量之间的非线性关系时获得了很高的精度.

最近, 随机森林模型被用于植被和干旱地区的 LST 降尺度研究[28], [39].此外, 将地表反射率、地形因素(DEM, 坡度)和土地利用图引入随机森林降尺度模型, 在模拟Landsat LST和MODIS LST产品中取得了很好的效果[51].但仍有一些细节问题值得讨论, 如随机森林降尺度模型在不同地区的应用, 以及对预测因子的选择应更全面和系统地考虑.因此, 本研究的主要目的是验证多类型预测变量的随机森林回归模型(MTVRF)在不同底层地面的LST降尺度上具有满意的表现, 并对MTVRF模型的稳健性进行了综合评价。本研究首先应用MTVRF模型将MODIS(990 m)重采样的LST产品降尺度到90 m，使用了四种预测变量，包括地表反射率、光谱指数、地形因子和土地分类图；然后用ASTER(90 m)的LST产品来评价其结果。为了进一步比较MTVRF模型的降尺度性能，我们选择了两个具有不同土地覆盖和地形因素的区分区域。本文的其余部分组织如下。第2节介绍了研究区域、数据和方法。第3节对降尺度的结果进行了详细分析。第4节得出了进一步的讨论，第5节得出了一个结论。

材料和方法

研究区

本研究选择了两个具有不同地形和土地覆盖的典型地区，以充分了解这两个地区的基本表面特征。图1显示了由Landsat-8 OLI反射率数据得出的这两个地区的假彩色图像。

利用Landsat 8数据生成的假彩色图像（R：波段5；G：波段4；B：波段3）对研究区进行地理位置定位。(a)西班牙塞戈维亚的研究区A，(b)中国北京的研究区B，(c)水区的子区1，(d)植被区的子区2，以及(e)不透水表面区的子区3。

研究区域A包括位于西班牙塞戈维亚的Peñarora山区。该地区的纬度和经度分别为北纬40度至41度和西经4度至3度。该地区80%以上的面积被植被覆盖，主要由森林和耕地组成。气候为地中海亚热带气候，温度在-17.0℃至39.7℃之间。这个地区的年平均温度为11.5℃，年降水量约为464毫米。

研究区B位于中国北京的西北部，包括海淀和昌平区。这个地区的坐标范围是北纬39度到40度，东经115度到116度。除西部山区外，该地区大部分地区地势平坦。气候是典型的半湿润大陆性季风气候，年平均温度为10℃至12℃，平均降水量为450毫米至550毫米。本研究区包含四种土地覆盖：植被、耕地、不透水表面（包括建筑物和道路）和水；该地区的主要特点是西部地区的植被和东部地区的不透水表面。

为了更详细地讨论LST降尺度，我们还选择了三个分区：研究区A的分区1，即Santillana水库，被标记为水区，最大地表覆盖面积为1052公顷；研究区B的分区2，以森林为主，被视为植被区；研究区B（即西城区）的分区3，有高层建筑和纵横交错的道路，可以说是不透水区域。

数据预处理

Landsat数据

Landsat 8 Operational Lad Imager（OLI）和TIRS图像是在USGS Earth Resources Observation and Science（EROS）中心科学处理架构下获取的，并由USGS Earth Explorer（http://earthexplorer.usgs.gov/）进行检索，分辨率分别为30米和100米。在这项研究中，两个Landsat瓦片（WRS-2 Path 201 / Row 32和Path 123 / Row 32）被应用于覆盖塞戈维亚和北京的选定研究区域。带有2-7个波段的Landsat地表反射率产品和经过处理的光谱指数被作为输入变量之一，用于将MODIS LST产品从990米降级到90米。

ASTER数据

更精细分辨率的LST来自ASTER，它是1999年12月18日发射的Terra卫星上的一个传感器。ASTER LST产品(AST08)可从NASA Earthdata Search(https://search.earthdata.nasa.gov/)获得, 空间分辨率为90米, 由温度/辐射分离(TES)算法生成, 精度约为1.5K[52].ASTER LST将作为参考数据来验证MTVRF模型在更精细尺度上的性能。

MODIS数据

研究区A的MODIS产品是在2016年6月22日获得的，研究区B的MODIS产品是在2014年7月24日获得的，它们来自Terra卫星上的另一个传感器。由于ASTER和MODIS传感器使用的是同一个卫星平台，LST产品中由几何观测偏差引起的差异往往被忽略。同时的数据也消除了ASTER和MODIS LST产品之间的观察时间差异。虽然MODIS的LST可以通过广义分割窗（GSW）算法[53]和TES算法[54]产生，但由于科学数据质量问题，无法获得基于TES的LST，我们只能使用基于GSW的传感器LST。收集的6个MOD11A1 LSTs是基于瓦片的全球产品，提供每像素的温度和发射率值，分辨率为1公里，精度约为1K[55], [56]。这些基于GSW的LST产品在研究区A登记到WGS 84/UTM 30 N区，在研究区B登记到WGS 84/UTM 50 N区，重采样间隔为990米。值得注意的是，由于生成LST产品的算法不同，ASTER TES的LST和MODIS GSW的LST之间存在偏差。幸运的是，一个简单的线性回归已经被用来消除不同传感器之间的系统性LST差异[57]。这个线性回归模型首先为 ASTER 和 MODIS LST 产品建立了较粗的分辨率。随后，用这种简单的线性回归方法对 MODIS LST 产品进行转换，然后应用于 MTVRF 模型训练，而不是原始 MODIS LST 产品。图2显示了ASTER和MODIS LSTs在990米分辨率下的散点图和线性回归关系。在95%的置信区间内，研究区A的决定系数（R2）为0.93，研究区B的决定系数为0.90，这表明ASTER和MODIS LST产品之间存在明显的相关性。研究区A的RMSE为1.82K，研究区B为1.07K。有60-70%的点的RMSE值在-1到1K之间, 这表明这种转换可能适合于评估降尺度LST的准确性.

图2.研究区域A和B的MODIS和ASTER的LST产品的散点图和线性回归关系。

SRTM DEM数据

数字高程模型（DEM）数据来自NASA的航天飞机雷达地形任务（SRTM），这是由NASA、国家国防局（NIMA）、德国和意大利航天局共同完成的。空间分辨率为90米的SRTM DEM数据也被重投影到WGS 84/UTM 30N区和WGS 84/UTM 50N区。随后，通过ARCGIS 10.2的空间分析模块计算出坡向、坡度和山体阴影。这些地形因素在空间上被重采样到990米的分辨率，通过空间平均化到MODIS LST的分辨率来进行RF模型训练。

方法

MTVRF

RF是一种综合机器学习算法，由bagging算法演变而来[58]-[60]。作为一种非线性统计集合回归方法，RF是由一组不相关的分类树和决策回归树构成。每一个自举样本都是从训练集中选取的，在训练树的每个模式时，使用的特征都是从所有特征中随机抽取一定比例的特征[61]。RF训练的结果变成了所有决策树的投票输出。RF有望获得几乎所有期望的结果，对数据和使用的参数具有高度的适应性，而且不需要像SVM那样繁琐地调整参数。RF被认为是当今最好的机器学习算法之一，已被广泛用于涉及遥感图像处理的各个领域，如遥感图像分类、特征识别和空间降尺度等[62]-[65]。

随机森林降尺度算法中的预测变量数据集是为了有效地反映不同地区LST的空间变化.在 TsHARP 和基本 RF 模型中, 只有植被覆盖率等与 LST 关联度最高的少数变量才有资格被引入模型中.然而，由于土地覆盖的复杂性，这些变量的性能可能较弱，在其他地区进行测试时，会进一步阻碍其效果。在本文中，MTVRF将在对多重共线性不敏感的前提下，对算法的复杂性和输入变量的数量进行权衡，这保证了结果对缺失和非平衡数据的稳健性，并对成千上万的输入有一个令人满意的预测。此外，MTVRF还可以避免过度拟合，从而提高了模型的通用性，因此它可以描述复杂的地表状态，并适用于其他地区，结果可以接受。

根据 LST 与生物物理参数的统计相关性, 在 MTVRF 降尺度模型中主要引入了 4 种预测变量以提高泛化能力.输入变量列表如下:

可见、近红外和短波红外波段的表面反射率，其中包含大量的植被覆盖和土壤湿度条件。
典型的光谱指数，可能对特定的土地覆盖类型敏感。例如，NDVI、植被分数和SAVI令人满意地表明了植被密度和生物量，而NDBI的训练是为了更精确地识别不透水表面，其他光谱指数如NDDI、NMDI和MNDWI被用来拟合地表温度和土壤水分之间的关系，这对LST的变化有很大影响。
地形因素，包括DEM及其衍生物，如研究区域的坡向、坡度和山影，被认为与山区的LST有显著的相关性。
土地分类图，被认为是促进认识不同地区土地覆盖对LST影响的预测器。

降尺度模型过程

开发的MTVRF LST降尺度模式的具体步骤如图3所示，并总结如下。

输入的变量应该进行汇总，以配合要降尺度的产品的粗略分辨率。在本研究中，首先将分辨率分别为30 m和90 m的Landsat 8 OLI和SRTM图像数据用空间平均法汇总到990 m（即与MODIS的LST数据相同），这样可以在粗放的水平上建立解释变量与LST之间的统计关系。
在990米的水平上，LST与解释变量之间的MTVRF回归模型可以表示为：。其中，下标c代表具有较粗分辨率的变量，下标i代表第i个变量，ρ代表反射率，S代表光谱指数，TF代表地形因子，LC代表土地覆盖类型。函数F(⋅)表示转换后的MODIS LST与这些变量之间的非线性关系。
值得注意的是, 由于 RF 回归无法解释 LST 分布的所有时空变化, 因此由更粗分辨率的原始 LST 与模拟 LST 之间的差异计算出的残余温度 (LSTc )是为了修正预测模型误差[51] . 其中下标o表示经ASTER转换的原始MODIS LST。
考虑到LST和其他变量之间关系的尺度不变性，训练后的模型随后被用于更细的尺度（90米）。重新取样后的模型误差被分配给更精细分辨率的每个像素点。因此，在 90 米分辨率下，最终降尺度的 LST 可按以下方式计算。其中下标f指的是具有更精细分辨率的变量。

随机森林模型最重要的参数之一是袋外误差，其中袋外（out-of-bag，OOB）样本是指没有参与模型训练的大约三分之一的样本。每个树中的袋外样本作为测试数据集，以确保对误差的无偏估计。因此，不需要执行交叉验证或使用单一的测试数据集来获得无偏估计的误差，因为无偏估计可以在模型的生成过程中建立。平均而言，每个研究地区的训练数据集的大小约为10000个。

根据OOB误差估计，对带来最佳回归结果的RF模型的关键参数，即决策树的数量（n_estimators）和要分割的特征的最大数量（max_features）进行遍历，最佳结果是n估计器为200，max_features为9。

MTVRF 降尺度模型的性能与以下模型进行了比较: (1) Hutengs 和 Vohland[51] 提出的基本 RF 模型, 该模型选择近红外和红色波段的反射率值、DEM 和土地利用图作为输入变量; (2) TsHARP 方法, 该方法基于 LST 和 NDVI 的线性关系.

误差评估

降尺度的结果是由R2、ME和RMSE来评价的，它们通常被用作拟合问题的可测量指标。其表达式如下。

其中LSTe指的是来自MTVRF的估计LST，LSTr代表参考ASTER LST，n代表参与估计的样本总数。

结果和分析

选择变量重要性

图4显示了MTVRF在两个研究区计算的输入变量的重要性排名。值得一提的是，这个重要性得分给出了一个关于输入变量贡献的相对排名，但它不等同于相关系数。地形因素的贡献，特别是DEM和坡度的贡献，在两个地区都显示为大于其他变量的贡献。这意味着沿山地表面的太阳入射辐射和长波表面冷却存在着巨大的地形影响。此外，以地形为主的研究区A，在DEM和其他种类的变量之间显示出明显的不一致，而研究区B的重要性得分由于地势平坦而显示出同质的变化。

对于研究区A，主要由海拔约616-2409米的丘陵覆盖，LST的空间分布几乎完全由这些地形因素和植被分量的相关因素控制，如红波段和NMDI。对于研究区B，在市中心上空获取的场景，地貌复杂，地形相对平坦，其DEM的重要性变得不如研究区A，分类图的权重更高。被视为城市中建筑物和道路指标的NDBI也显示出较高的重要性。

因为这些重要性分数随着输入变量的数量而变化，也就是说，当输入变量被添加或删除时，这些重要性分数会发生变化，它们可能为选择具有较高相关性的输入变量提供证据。为了进一步分析，应进行量化分析以评估MTVRF的降尺度效果。

不同地表情况下的降尺度表现

表 1 显示了 RF 模型从 990m 到 90m 的 LST 降尺度表现, 三种模型都提高了 LST 从 990m 到 90m 的分辨率, 且精度令人满意.平均而言, MTVRF 模型获得的降尺度效果最好, 其次是基本 RF 模型和 TsHARP 方法.