复杂山区地形条件下ERA5再分析地表气温降尺度方法

图1 模型试验区和模型移植检验区

Fig. 1 The geographical location of the research area

2.2 研究数据来源

（1） ERA5再分析2 m高度气温数据

ERA5是ECMWF推出的第五代大气再分析产品,相比前一代产品ERA-Interim,时间和空间分辨率都有了较大提高。ERA5逐时再分析2 m高度气温是陆地、海洋或内水表面2 m高处空气的温度,从最低模式层到地表进行插值计算得到,空间分辨率是0.1 °,来源于欧洲气象中心（http://apps.ecmwf.int/datasets/）。

（2）自动气象站数据

逐时气温观测数据来自于2020年1月河北省张家口市的268个自动气象站和模型移植检验区的146个自动气象站,由河北省气象局提供,气象站位置如图1所示。数据内容包括站点经纬度位置、海拔和逐时气温。气象站观测数据主要用于降尺度模型构建和降尺度结果验证。

（3）Landsat8卫星影像数据

Landsat 8 是美国陆地卫星计划的第八颗卫星,于2013年2月发射升空,Landsat8搭载OLI和TIRS2个传感器,其中OLI陆地成像仪获取的多光谱影像分辨率是30 m,热红外传感器 TIRS包括 2个热红外波段,分辨率100 m。本文使用的数据下载于美国地质勘探局网站（https://glovis.usgs.gov/）,所用影像信息如表1所示。热红外影像用于计算地表温度（Land Surface Temperature,LST）、多光谱影像用于计算归一化植被指数（Normalized Difference Vegetation Index, NDVI）、归一化雪被指数（Normalized Difference Snow Index, NDSI）、地表光谱反照率（Albedo）等表征下垫面物理特性的地表参数,用于气温降尺度模型的输入分量。根据美国地质调查局（United States Geological Survey,USGS）公告,Landsat 8 TIRS传感器第11波段的定标参数存在较大不确定,参考胡德勇等^[33]、徐涵秋等^[34]的研究,选择第10波段数据利用单窗算法反演地表温度。

表1 本文使用的Landsat8遥感影像

Tab. 1 Landsat8 remote sensing image used in the research

覆盖范围	轨道号	成像日期
方法研究试验区	124_031	2020-01-31
	124_032	2020-01-31
	125_031	2020-01-22
	125_032	2020-01-22
模型移植检验区	122_032	2020-01-01

（4） SRTM DEM数据

SRTM（Shuttle Radar Topography Mission）DEM数据由美国航空航天局和国防部国家测绘局联合测量,空间分辨率为30 m,数据来源于美国地质勘探局（https://lpdaac.usgs.gov/）。DEM数据用于计算坡度（Slope）、坡向（Aspect）、天空开阔度(Sky View Factor, SVF)等地形因子,作为气温降尺度的输入分量。

3 研究方法

从气候层面讲,影响复杂地形条件下气温空间分布的因素主要包括3个方面：① 海拔、经纬度、大范围气候模式以及高大山脉的走向等宏观地理条件;② 坡度、坡向、地形遮蔽度等局部地形因子; ③ 土壤类型、植被状况等下垫面性质。再分析资料能够对宏观地理条件影响下的大尺度气温分布有较为准确的描述,而缺乏对地形因子和下垫面影响下的局部气温分布的描述;气象站能够获得其对应位置下垫面和局部地形影响下的气温,却因为其分布不均匀的特点,难以表达连续的气温空间分布特征,因此,准确获取山区气温空间分布需要综合使用再分析资料和气象站数据。对比再分析资料与站点观测气温数据发现,与逐时再分析气温数据相比,日均再分析气温具有更好的稳定性和可靠性,与实测日均气温更为接近,可基于气温在一天中的周期性变化规律建立日均气温与逐时气温的关系模型,实现逐时再分析气温降尺度。综合以上分析,本文提出了一种将局部薄盘样条函数插值与随机森林模型模拟相结合的降尺度方法（TPS-RF）,该方法利用再分析资料和站点数据,并辅以Landsat8影像和DEM数据,在0.1 °日均再分析气温插值的基础上,构建日均气温订正模型和逐时气温模拟模型,获取30 m分辨率的逐时气温空间分布,并开展降尺度方法的时空移植效果分析,研究路线如图2所示。

图2

图2 TPS-RF降尺度方法技术路线

Fig. 2 The flowchat of TPS-RF downscaling method

3.1 基于局部薄盘样条函数的日均再分析气温空间插值

局部薄盘光滑样条函数是样条函数的一个扩展^[35],它在引入普通样条函数自变量的基础上,允许加入线性协变量。局部薄盘样条函数模型如下：

（1）

z_{i} = f (x_{i}) + b^{T} y_{i} + e_{i} (i = 1,2, \dots, N)

式中： $z_{i}$ 是位于空间上第 $i$ 点的因变量; $x_{i}$ 为d维样条独立变量。 $f$ 是要估算的关于 $x_{i}$ 的未知光滑函数; $y_{i}$ 为p维独立协变量; $b$ 为 $y_{i}$ 的p维系数; $e_{i}$ 为具有期望值为0,方差为 $w_{i} σ^{2}$ 的自变量随机误差; $w_{i}$ 为已知的相对误差方差, $σ^{2}$ 为误差方差,在所有数据点上为常数,但通常未知。由式（1）可知,当式中不包含协变量 $y_{i}$ 时,模型简化为简单的薄盘样条函数;当式中缺少 $x_{i}$ 时,即没有独立自变量时,模型变化为多元线性回归。

模型中函数 $f$ 和系数 $b$ 可通过最小二乘估计使得下式中的 $β$ 最小来确定。

（2）

β = \overset{N}{\sum_{1}} {(\frac{Z_{i} - f (x_{i}) - b^{T} y_{i}}{w_{i}})}^{2} + ρ J_{m} (f)

式中： $J_{m} (f)$ 是函数 $f$ 的粗糙测度函数,定义为函数 $f$ 的m阶偏导; $ρ$ 是正的光滑参数,在数据保真度和曲面粗糙度之间起平衡作用,通常由广义交叉验证GCV（Generalized Cross Validation）的最小化来确定。GCV计算可采用“one point move”方法,一次移去一个样点,用剩余样点在一定的光滑参数下进行曲面拟合得到该点的估计值,再计算观测值的方差。

本文使用逐时再分析气温数据计算日均气温,以经纬度为自变量,DEM为协变量,日均气温为因变量,应用局部薄盘样条函数训练插值模型,并将30 m分辨率DEM数据代入模型,得到30 m分辨率的日均气温。

3.2 基于随机森林算法的日均再分析气温插值结果订正

局部薄盘样条函数插值之后的日均再分析气温空间分布更为精细,但在局部细节纹理的描述和气温数值大小上缺乏准确性,其原因是用于插值的再分析气温数据的空间分辨率低且存在一定程度的误差。为了获取准确高空间分辨率气温数据,研究中采用随机森林方法构建日均气温订正模型,对局部薄盘样条函数空间插值结果进行订正。

随机森林方法最早由Breiman^[35]提出,是由多棵分类与回归决策树CART组合而成的一种新型机器学习算法^[35]。随机森林模型由众多CART决策树构成,模型决策能力取决于每一棵CART决策树,每一棵决策树由根节点、子节点和叶子节点组成,每个子节点包含一个评判规则,叶子节点则对应一个评判级别。随机森林算法的步骤大致分为以下4步：① 从原始数据集中抽取训练集;② 为每个训练集建立决策树,众多决策树就构成了森林;③ 设随机森林模型中共有M个风险指标,随机抽取m个（m≤M）作为节点指标,Breiman推荐 $m = \sqrt[]{M}$ ,在这m个节点指标中选择基尼系数最小作为该节点的分支标准;④ 最后根据所有决策树的预测结果,采用投票方式来决定新样本类别。每次抽样未被抽中的数据构成了袋外数据（Out-of-Bag,OOB）,利用这些袋外数据估计内部误差,称之为袋外误差（Error of Out-of-Bag, E_OOB）,计算公式如下：

（3）

E_{OOB} = \frac{1}{n} \overset{n}{\sum_{i = 1}} {[\hat{Y} (X_{i}) - Y_{i}]}^{2}

式中：n为OOB样本个数; $\hat{Y} (X_{i})$ 为根据给定样本通过随机森林算法计算的输出数据; $Y_{i}$ 为观测数据。

影响局地气温空间分布的主要因素是太阳直接辐射、地表热辐射和周围地物的反射,这3种辐射都是太阳、局地地形和下垫面相互作用的结果。为了描述这3种辐射,本文基于DEM计算了影响太阳辐射的坡度、坡向和天空开阔度,基于Landsat8遥感影像计算了反映下垫面发射及反射特性的归一化植被指数、归一化雪被指数、地表光谱反照率和与气温空间分布特征相似性较高的地表温度,各地形、地表参数计算方法如表2所示。本文对逐日平均气温与地表温度做了相关性检验,结果表明都达到了P<0.05的显著水平,二者之间的相关性平均为0.74,但为了进一步提高日均气温估算模型的精度,建立了气象站日均气温和对应位置地表温度之间的一元二次方程,对遥感反演的地表温度进行拟合得到LST',以提高日均温与地表温度的相关性。

（4）

T_{oi} = a \times LS {T_{i}}^{2} + b \times LS T_{i} + c

（5）

LST' = a \times LS T^{2} + b \times LST + c

式中： $T_{oi}$ 是第 $i$ 个站点的日均气温; $LS T_{i}$ 是第 $i$ 个站点位置对应的地表温度;LST是TIRS反演地表温度; $a 、 b 、 c$ 分别是一元二次方程的拟合系数。

表2 日均气温插值订正所用的地形地表参数

Tab. 2 Land surface parameters used to correct daily mean air temperature

参数名称	参数描述	参考文献
NDSI	$NDSI = (ρ_{Green} - ρ_{Swir 1}) / (ρ_{Green} + ρ_{Swir 1})$	[36]
NDVI	$NDVI = (ρ_{NIR} - ρ_{Red}) / (ρ_{NIR} + ρ_{Red})$	[37]
Albedo	$ρ_{short} = 0.356 \times ρ_{Coastal} + 0.13 \times ρ_{Green} + 0.373 \times ρ_{Red} + 0.085 \times ρ_{NIR} + 0.072 \times ρ_{Swir 2}$	[38]
LST'	大气与地表结合部的温度状况,和气温具有较高的相关性	[39-40]
DEM	通过有限的高程数据对地面地形的数字化模拟。气温随高程的增加而降低	[41-42]
Slope	坡面的垂直高度和水平方向的距离之比	[43]
Aspect	坡面在水平面上的投影方向	[43]
SVF	$SVF = 1 - \overset{n}{\sum_{i = 1}} si n^{2} β ∙ (\frac{α}{2 π})$	[44]

注： $ρ short$ 表示宽波段光谱反照率, $ρ Coastal$ 、 $ρ Green$ 、 $ρ Red$ 、 $ρ NIR$ 、 $ρ Swir 1$ 、 $ρ Swir 2$ 分别表示Landsat 8遥感影像的海岸波段、绿波段、红波段、近红外波段、短波红外1、短波红外2的光谱反射率, $α$ 表示方位角间距, $β$ 为在一个方位角间距范围内的最大山体高度角。

由于气象站空间分布不均且仅能观测对应位置的气温,在山区复杂地形条件下仅依靠气象站数据构建气温空间化模型,无法得到准确的气温空间分布,且该模型无法适用于不具备气象站点的区域。为解决这个问题,本文在日均再分析气温局部薄盘样条函数插值的基础上,计算气象站观测气温与其对应位置插值气温的比例系数 $k$ 。

（6）

k = \frac{T_{oi}}{T_{i}}

式中： $T_{oi}$ 是第 $i$ 个站点的气象站气温; $T_{i}$ 是第 $i$ 个站点位置的再分析插值气温。

以 $k$ 为因变量,以同一位置的各地表地形参数为自变量,构建基于随机森林算法的日均气温订正模型,并将模型应用到整个研究区,得到订正后的日均气温降尺度结果 $T_{day}$ 。

为了获取最优降尺度结果,分别对随机森林算法中的决策树数量和每个决策树随机选择的特征数目进行参数调优。一般来讲,决策树数量越多结果越好,但当决策树到达一定数量后,再增加决策树数量对结果的影响不大,并且会增加计算量。本文通过MATLAB计算了误差随决策树数量和特征数目的变化情况,考虑误差大小和计算量两方面因素,最终确定决策树数量为100,特征数目为4。

3.3 基于随机森林算法的逐时气温估算

在没有突发天气状况的条件下,各位置气温在一天中的变化呈现出一定的周期性规律,但山区地形复杂,太阳辐射角度和坡度、坡向的相互作用使得各位置气温在周期变化的基础上又表现出独特性特点。本研究认为逐时气温是在日均气温的基础上加上或者减去一定数值 $Δ T (t)$ 得到,但由于不同空间位置的地表地形要素存在着差异,所以不同位置的 $Δ T (t)$ 不同。为了描述气温随时间变化的数值大小和空间差异,并提高模型的时空移植性能,论文将 $Δ T (t)$ 分解为气温平均变化幅度 $Δ T_{era} (t)$ 和气温变化的空间分布差异 $Δ T_{oe} (t)$ 两部分。研究中通过再分析数据计算 $Δ T_{era} (t)$ 来表达研究区域t时刻气温相对于日均气温的平均变化幅度,通过气象站观测数据计算 $Δ T_{oi} (t)$ 来表示站点位置t时刻气温相对于日均气温的变化大小,进而计算 $Δ T_{oi} (t)$ 与 $Δ T_{era} (t)$ 的差异并将其空间化得到 $Δ T_{oe} (t)$ ,来表达 $Δ T (t)$ 的空间分布差异。具体计算过程如下：

（1）首先计算 $t$ 时刻再分析气温与日均再分析气温之差在整个区域上的平均值 $Δ T_{era} (t)$ ;再计算 $t$ 时刻站点观测气温与观测日均气温之间的差值 $Δ T_{oi} (t)$ ,并计算 $Δ T_{oi} (t)$ 与 $Δ T_{era} (t)$ 的差值 $∆ T_{oei} (t)$ ,计算公式如（7）所示。

（7）

\begin{array}{l} ∆ T_{oei} (t) = Δ T_{oi} (t) - Δ T_{era} (t) = \\ [T_{oi} (t) - \frac{\overset{24}{\sum_{t = 1}} T_{oi} (t)}{24}] - [{\bar{T}}_{era} (t) - \frac{\overset{24}{\sum_{t = 1}} {\bar{T}}_{era} (t)}{24}] \\ (t = 1, 2, 3, \dots, 24) \end{array}

式中： $T_{oi} (t)$ 是第 $i$ 个站点 $t$ 时刻的观测气温; ${\bar{T}}_{era} (t)$ 是研究区 $t$ 时刻再分析气温的平均值。

（2）以各地形地表参数为自变量,以 $∆ T_{oei} (t)$ 为因变量,构建随机森林模型,并将模型应用至整个研究区域,得到空间分布的 $Δ T_{oe} (t)$ 。

（3）假设订正后的日均再分析气温与实测日均气温非常接近,可计算得到t时刻30 m分辨率的气温空间分布,公式如下：

（8）

T_{t} = T_{day} + Δ T (t) = T_{day} + Δ T_{era} (t) + Δ T_{oe} (t)

式中： $T_{day}$ 是订正后的日均再分析气温降尺度结果。

3.4 降尺度结果精度验证

选择平均绝对误差（Mean Absolute Error,MAE）、均方根误差（Root Mean Square Error,RMSE）、平均偏差（Mean Bias Error,MBE）和相关系数r,对再分析气温降尺度方法及其时空移植效果进行精度评定。平均绝对误差反映误差的大小,均方根误差可以反映样点数据的估值灵敏度和极值效应,这2个误差值越小,表明模型模拟精度越高;平均偏差反映预测值与实际值的总体相对大小关系,其绝对值越小,预测结果就越接近实际值;相关系数r表示估计值与实测值相关程度的度量,其值越接近1表示结果越好。

4 结果及分析

4.1 ERA5再分析气温数据精度分析

根据逐时再分析数据计算日均气温和月均气温,随后基于局部薄盘样条函数分别将3种时间尺度的0.1 °气温降尺度到30 m。选取张家口市268个气象站气温数据和对应位置再分析资料气温降尺度结果,分别计算了二者在月均、日均和逐时时间尺度上气温的相关系数r、平均绝对误差MAE和均方根误差RMSE。

从逐时、日均到月均的时间尺度变化上,再分析资料与气象站气温的相关性逐渐升高,误差逐渐减小（表3）。在月均时间尺度上,再分析资料与气象站气温高度相关,相关性达到了0.88;在日均时间尺度上,二者在1月各天相关性的平均值是0.80,除了1月5日和21日这两天相关性较低之外,再分析资料与气象站气温都表现出了较高的相关性;在逐时时间尺度上,再分析资料与气象站气温相关性的平均值能达到0.70以上,但是随着时刻的变化,相关性大小有较大幅度的波动（图3）,表现出了不稳定的特点。再分析资料与气象站气温在月均时间尺度上虽然相关性高,但月均气温向逐时气温转换的时间尺度跨度大,难以建立二者之间的稳定关系。因此,结合再分析资料气温在日均时间尺度上的可靠性和气温在一天中的周期性变化规律,本文在再分析日均气温降尺度基础上,构建逐时气温模型估算30 m分辨率气温数据。

表3 不同时间尺度上再分析气温精度分析

Tab. 3 Accuracy analysis of reanalysis temperature on different time scales

指标	月均气温	日均气温	逐时气温
r	0.88	0.80	0.71
MAE/℃	2.40	2.70	3.10
RMSE/℃	3.00	3.30	3.80

图3

图3 不同时间尺度上再分析资料与气象站气温相关性

注：横坐标表示1月2日0:00至1月31日23:00的720个时刻。

Fig. 3 Correlation between reanalysis data and temperature of meteorological stations on different time scales

4.2 1月14日试验区再分析逐时气温降尺度结果分析

为了检验降尺度方法在试验区2020年1月逐时气温估算方面的效果,选择了1月14日的数据构建模型,开展降尺度方法研究。首先利用0.1°分辨率的再分析资料计算日均气温,并基于局部薄盘样条函数将其降尺度到30 m,随后随机选取试验区内40个自动气象站气温数据构建日均气温订正模型和逐时气温估算模型,获得30 m分辨率的逐时气温空间分布。分别使用用于训练的40个站点数据和剩余的228个站点数据计算降尺度前后的平均绝对误差MAE和平均偏差MBE（表4）,并以2、8、14、20时（北京时间,以下相同）为例,图4给出了气温空间分布降尺度前后对比结果。

表4 再分析资料降尺度前后精度对比

Tab. 4 Error comparison of reanalysis air temperature data before and after downscaling (℃)

数据	MAE	MBE
0.1°分辨率再分析气温	3.88	3.53
降尺度后30 m分辨率气温（训练集）	1.22	-0.09
降尺度后30 m分辨率气温（验证集）	2.13	0.47

图4

图4 试验区降尺度前后2020年1月14日气温空间分布

Fig. 4 Spatial distribution of air temperature on Jan 14, 2020 before and after downscaling in the test area

根据图4,通过降尺度,气温空间分布的精细化程度得到了显著提高,细节纹理更加丰富,尤其是地形起伏较大的中部山区,气温空间分布的纹理与地形纹理有很高的相似性。随着时间推移,山区气温相对高低会发生较大变化,这既是太阳辐射与地形、下垫面相互作用的结果,也与冷热气流运动密切相关。在地形起伏较小的南部城乡建设用地和北部坝上高原,气温空间分布与下垫面属性更为相关,下垫面接收太阳辐射后的反射和发射决定了气温大小和空间分布上的相对高低。根据表4,降尺度之后气温的MAE和MBE相较于原始ERA5再分析气温都有了较大幅度的下降,表明降尺度之后的气温精度有了显著提高。再分析资料的MBE达到了3.53 ℃,说明再分析气温数据在该地区该时间段的气温高于实测气温,降尺度之后的气温MBE下降到了0.47 ℃,与气象站气温更加接近,这说明降尺度在一定程度上解决了再分析气温在该区域偏高的问题。

为了验证本文降尺度方法（TPS-RF）的结果精度,将其与基于薄盘样条函数的气象站气温插值（TPS-Sta）、基于薄盘样条函数的再分析逐时气温插值并根据气象站数据进行一元线性模型订正（TPS-Line）的方法进行对比,部分时刻的结果如图5所示,利用未参与建模的气象站数据计算各降尺度结果的MAE、RMSE和相关系数r,结果如表5所示。其中TPS-Line的一元线性模型如式（9）所示。

（9）

T_{oi} (t) = u \times T_{ei} (t) + v

式中： $T_{oi} (t)$ 是 $t$ 时刻第 $i$ 个气象站点获取的气温, $T_{ei} (t)$ 是 $t$ 时刻第 $i$ 个气象站点位置处再分析资料插值得到的气温, $u$ 和 $v$ 分别是一元线性模型的拟合系数。

图5

图5 试验区2020年1月14日部分时刻3种方法降尺度结果对比

Fig. 5 Comparison of downscaling results of different methods in the test area on January 14, 2020

表5 试验区2020年1月14日不同降尺度方法精度

Tab. 5 Error comparison of different downscaling methods in the test area on January 14, 2020

指标	TPS-RF	TPS-Sta	TPS-Line
RMSE/℃	2.68	4.17	2.92
MAE/℃	2.13	3.61	2.33
r	0.82	0.84	0.79

根据图5可知,3种降尺度方法都能描述气温随海拔和纬度升高而逐渐降低的特点,但是3种方法得到的局部区域空间分布结果差异较大。根据表5可知,由于气象站分布及其空间代表性问题,TPS-Sta降尺度方法的结果误差最大,研究中提出的TPS-RF方法精度最高。对比TPS-RF降尺度结果和TPS-Line降尺度结果可以发现,二者最大区别在细节纹理信息的描述上,TPS-RF降尺度结果有更好的细节纹理信息（图5(a)、图5(d)、图5(g)、图5(j)）,而TPS-Line方法只能得到该区域的大致气温分布,无法给出气温空间分布细节信息（图5(c)、图5(f)、图5(i)、图5(l)）。对比TPS-RF降尺度结果和TPS-Sta插值结果,二者在表达气温空间分布相对高低和细节纹理上都有较大差异,其中TPS-Sta插值方法无法准确表达北部坝上高原和南部城乡建设用地的气温空间分布细节(图5(e))。根据图5,在3种降尺度方法中,TPS-RF降尺度得到的气温空间分布值域范围更大,与实际情况更为吻合,在复杂地形条件下的山区具有更高精度,而其他两种方法分别受局部薄盘样条函数、再分析资料分辨率低和站点数量少的限制,不能准确表达高温和低温区域。

利用228个气象站点气温作为验证数据,计算1月14日TPS-RF降尺度结果的均方根误差,并将其插值到整个研究区,得到误差的空间分布（图6）,并通过比较各时刻降尺度结果误差、再分析气温与实测气温之间的相关性（图7）,分析误差产生的可能原因。

图6

图6 2020年1月14日试验区TPS-RF降尺度结果误差空间分布

Fig. 6 Error distribution for downscaling result of TPS-RF in the test area on January 14, 2020

图7

图7 逐时降尺度结果误差分析

Fig. 7 Error analysis diagram of hourly downscaling results

根据图6和图7可知,逐时降尺度结果误差随时间变化具有一定规律性,其中下午时刻降尺度结果误差最小,主要原因可能是降尺度方法中考虑了地形、下垫面对太阳辐射的影响,而对于没有太阳辐射的夜间时刻,气温空间分布及变化主要受地表长波辐射以及水平气流的影响,研究中提出的降尺度方法未能充分考虑这种影响,因此误差相对较大。根据图7可知,再分析资料与实测气温的相关性与降尺度结果误差存在着一定关联性,在一段时间范围内,相关性越高则再分析资料越可靠,降尺度结果误差越小。根据图6可知,误差在空间分布上的特征与变化规律性不明显,因此,再分析资料的准确程度是降尺度误差大小的最主要影响因素。

气象站数据、Landsat遥感影像、DEM的数据本身误差,地表地形参数的计算过程产生的误差,日均气温订正模型和气温逐时化模型构建过程产生的误差最终都会转化为结果误差。另外,站点气温数据和降尺度结果分别代表气象站所在点的气温和30 m分辨率栅格区域的气温,在空间尺度上有不同的代表性,这也是误差的来源之一。最后,本文降尺度方法缺乏对地表长波辐射和水平气流的表达也在一定程度上增加了夜间降尺度结果误差。

4.3 TPS-RF降尺度方法的时空移植精度分析

将TPS-RF方法移植到其他时间和区域,以检验该方法在气象站点较少、甚至没有气象站点时用于ERA5再分析气温降尺度的效果。

4.3.1 方法的时间移植精度分析

将TPS-RF降尺度方法移植用于其他日期时,设置了有气象站观测数据和没有观测数据两种对照情况。当假设有气象站点时,对比本文提出的降尺度方法（TPS-RF）、基于局部薄盘样条函数的气象站气温插值（TSP-Sta）和基于局部薄盘样条函数的再分析资料插值及其线性订正3种降尺度方法（TSP-Line）的降尺度结果;当假设没有气象站观测数据时,将1月14日构建的日均气温订正模型和逐时气温模型应用到1月的其他日期（TPS-RF-time）,实现降尺度模型的时间移植,并与基于局部薄盘样条函数的再分析资料插值（TPS-era5）结果相比较。

以2020年1月20日为例,表6比较了不同降尺度方法得到的气温空间分布结果,表7给出了精度分析结果。根据表6和表7可知,当假设有足够的气象站数据时,本文提出的降尺度模型（TPS-RF）和再分析资料插值后订正方法（TPS-Line）可获得更加精细的气温空间分布,但气象站气温数据插值(TPS-Sta)结果仅与DEM空间分布相似,而未考虑坡度、坡向等地形要素和下垫面属性的影响,与实际气温空间分布差异较大,且结果误差最大,从误差大小和空间分布来看,TPS-RF方法能够获得最优的降尺度结果。当假设没有气象站观测数据时,基于1月14日观测数据训练的模型（TPS-RF-time）表现出了较好的时间移植性,能够获得更为精细的气温空间分布,RMSE较再分析资料的直接插值（TPS-era5）有1.04 ℃左右的降低。对比无气象站数据的模型时间移植结果和有气象站数据的降尺度结果,TPS-RF降尺度模型迁移误差增加了0.52 ℃左右（表7）,主要原因是其他日期与建模日期（1月14日）的天气状况不完全相同,气温随时间的变化规律也不完全一致;另外,随着日期的变化,太阳高度角也随之变化,太阳辐射与复杂地形下的下垫面之间的相互作用也会发生细微变化,也在一定程度增大了误差。模型降尺度结果在细节纹理上与观测数据具有一定差异,但降尺度结果与站点观测数据的空间分布总体相似,误差增大在可被接受范围,表明TPS-RF模型具有较好的时间移植性,可用于无气象站观测数据情况下的再分析气温降尺度。

表6 2020年1月20日试验区各方法降尺度空间分布结果对比

Tab. 6 Comparison of downscaling spatial distribution results of different methods in the test area on January 20, 2020

模型	2时	8时	14时	20时

表7 试验区2020年1月各降尺度方法精度对比

Tab. 7 Comparison of different downscaling methods in the test area in January 2020

是否使用站点观测数据	方法	RMSE/℃	MAE/℃	r
否	TPS-era5	3.93	3.14	0.74
否	TPS-RF-time	2.89	2.35	0.77
是	TPS-RF	2.37	1.88	0.80
	TPS-Line	2.63	2.10	0.73
	TPS-Sta	3.50	3.10	0.87

4.3.2 方法的空间移植精度分析

在验证模型的空间移植效果时,首先将模型移植检验区的再分析日均气温降尺度到30 m,然后使用试验区对应时刻的气象站数据训练日均气温订正模型和逐时气温模型,并将模型迁移到检验区域,实现检验区的再分析气温降尺度。与模型时间移植分析一样,设置检验区有、无气象站观测数据2组对照组,对比不同方法的降尺度效果,结果如表8和表9所示。

表8 模型空间移植检验区域的多种降尺度结果对比

Tab. 8 Comparison of different downscaling results by transfering the model from the test area to the checking region

模型	2时	8时	14时	20时

表9 模型空间移植检验区域1月各降尺度方法精度对比

Tab. 9 Comparison of different downscaling methods for the checking region in January

是否使用站点观测数据	方法	RMSE/℃	MAE/℃	r
否	TPS-era5	2.69	2.12	0.72
否	TPS-RF-area	2.41	1.93	0.73
是	TPS-RF	2.00	1.53	0.76
	TPS-Line	2.00	1.56	0.74
	TPS-Sta	2.25	2.00	0.86

根据表8和表9可知,模型空间移植表现出了和模型时间移植较为相似的结果。假设有气象站观测数据时,TPS-RF降尺度结果的MAE和RMSE较TPS-Sta分别有0.47 ℃和0.25 ℃的减小,空间细节纹理较TPS-Line方法更为精细,综合气温空间分布和降尺度结果误差来看,TPS-RF模型能够获得最优的降尺度结果。当假设没有气象站观测数据时,模型迁移（TPS-RF-area）效果在精度和精细化程度上都优于基于局部薄盘样条函数的再分析资料空间插值（TPS-era5）。模型空间迁移的降尺度结果（RMSE）较TPS-RF方法有0.41 ℃左右的增加（表9）,主要原因是模型的训练和应用区域不同,下垫面属性和地形起伏都相差较远;另外,两块区域与海洋之间的距离差异较大,这也是误差增大的原因之一。相比于模型的时间移植（表7）,模型的空间移植结果误差更小,主要原因是模型迁移区的海拔高度和地形起伏与试验区相比更小、再分析气温误差相对更小,因此降尺度模型也更为稳定。

5 结论与讨论

利用2020年1月份模型试验区和移植检验区的再分析资料、气象站气温数据、DEM和Landsat8影像,本文构建了局部薄盘样条函数和随机森林算法相结合的TPS-RF降尺度模型,开展复杂地形条件下的ERA5再分析气温降尺度方法研究,以获取30 m分辨率的逐时气温数据。为了检验TPS-RF降尺度模型的时空移植效果,设置了对照试验对比了在有、无气象站观测数据情况下的多种降尺度方法精度。

论文研究结果表明,TPS-RF模型可用于山区复杂地形条件下的ERA5再分析气温降尺度,降尺度后的再分析气温具有准确的空间分布与丰富的细节纹理信息,降尺度结果MAE和MBE相比原始再分析数据分别有1.75 ℃和3.06 ℃明显下降,有效解决了再分析气温偏高的问题;本文降尺度模型有较好的时空移植性,模型的时空移植结果在精度和空间分布上都优于再分析数据的直接插值,可用于气象站数据较少甚至没有气象站数据情况下的再分析数据降尺度。

本文创新性的提出了局部薄盘样条函数与随机森林算法相结合的TPS-RF统计降尺度方法,该方法将表示大尺度气温变化的再分析资料和精准表示局部气温的气象站数据相结合,并辅以高分辨率DEM、Landsat8影像来获取高时空分辨率的气温空间分布。该方法所用数据易于获取,且具有明确的物理意义,表明TPS-RF具有较高的适用性和可靠性。针对山区地形复杂,气温在时间和空间上变化剧烈的特点,本文以逐时刻30 m分辨率的气温空间分布作为降尺度目标,在时间和空间尺度上都优于以往的气温降尺度研究。

对于不同季节的气温降尺度和模型的时间移植,由于各地表参数空间分布、太阳辐射与地形的相互作用都发生了较大变化,需要对特征变量、模型参数等做相应调整。此外,本文降尺度模型只考虑了宏观地理条件、局部地形因子、地形下垫面属性,却没有考虑风向、风速、气流等因素对气温降尺度误差带来的影响,需要后续进一步开展相关研究工作。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Gao

, Bernhardt

, Schulz

Elevation correction of ERA-Interim temperature data in complex terrain[J]. Hydrology and Earth Syst Gao L, Bernhardt M, Schulz K. Elevation correction of ERA-Interim temperature data in complex terrain

[J]. Hydrology and Earth System Sciences, 2012,16(12):4661-4673. DOI: 10.5194/hess-16-4661-2012

[2]

Suggitt A

, Wilson R

, Isaac

N J B

, et al.

Extinction risk from climate change is reduced by microclimatic buffering

[J]. Nature Climate Change, 2018,8(8):713-717. DOI: 10.1038/s41558-018-0231-9

Protecting biodiversity against the impacts of climate change requires effective conservation strategies that safeguard species at risk of extinction(1). Microrefugia allowed populations to survive adverse climatic conditions in the past(2,3), but their potential to reduce extinction risk from anthropogenic warming is poorly understood(3-5), hindering our capacity to develop robust in situ measures to adapt conservation to climate change(6). Here, we show that microclimatic heterogeneity has strongly buffered species against regional extirpations linked to recent climate change. Using more than five million distribution records for 430 climate-threatened and range-declining species, population losses across England are found to be reduced in areas where topography generated greater variation in the microclimate. The buffering effect of topographic microclimates was strongest for those species adversely affected by warming and in areas that experienced the highest levels of warming: in such conditions, extirpation risk was reduced by 22% for plants and by 9% for insects. Our results indicate the critical role of topographic variation in creating microrefugia, and provide empirical evidence that microclimatic heterogeneity can substantially reduce extinction risk from climate change.

[3]

Gerlitz

Using fuzzified regression trees for statistical downscaling and regionalization of near surface temperatures in complex terrain: A case study from Khumbu Himal

[J]. Theoretical and Applied Climatology, 2015,122(1-2):337-352. DOI: 10.1007/s00704-014-1285-x

[4]

Di-Vittorio A

, Miller N

Evaluating a modified point-based method to downscale cell-based climate variable data to high-resolution grids

[J]. Theoretical and Applied Climatology, 2013,112(3-4):495-519. DOI: 10.1007/s00704-012-0740-9

[5]

Dirksen

, Knap W

, Steeneveld G

, et al.

Downscaling daily air-temperature measurements in the Netherlands

[J]. Theoretical and Applied Climatology, 2020(1-2):1-17. DOI: 10.1007/s00704-020-03313-1

[6]

Tang

, Ali

, Ma

, et al.

Downscaling of ERA-Interim temperature in the contiguous United States and its implications for Rain-Snow partitioning

[J]. Journal of Hydrometeorology, 2018,19(7):1215-1233. DOI: 10.1175/JHM-D-18-0041.1

[7]

段平, 盛业华, 李佳, 等.

自适应的IDW插值方法及其在气温场中的应用

[J]. 地理研究, 2014,33(8):31-38.

[ Duan

, Sheng Y

, Li

, et al.

Adaptive IDW interpolation method and its application in the temperature field

[J]. Geographical Research, 2014,33(8):31-38. ]

[8]

Merwade V

, Maidment D

, Goff J

Anisotropic considerations while interpolating river channel bathymetry

[J]. Journal of Hydrology, 2006,331(3):731-741. DOI: 10.1016/j.jhydrol.2006.06.018

[9]

Lu G

, Wong D

An adaptive inverse-distance weighting spatial interpolation technique

[J]. Computers & Geosciences, 2008,34(9):1044-1055. DOI: 10.1016/j.cageo.2007.07.010

[10]

Zhang

, Liu

, Zheng

, et al.

Regional precipitation model based on geographically and temporally weighted Regression Kriging

[J]. Remote Sensing, 2020,12(16):2547. DOI: 10.3390/rs12162547

[11]

Cort

, Willmott, et al.

Climatologically aided interpolation (CAI) of terrestrial air temperature

[J]. International Journal of Climatology, 1995,15(2):221-229. DOI: 10.1002/joc.3370150207

[12]

李莎, 舒红, 徐正全.

利用时空Kriging进行气温插值研究

[J]. 武汉大学学报·信息科学版, 2012,37(2):237-241.

[ Li

, Shu

, Xu

Interpolation of temperature based on spatial-temporal Kriging

[J]. Geomatics and Information Science of Wuhan University, 2012,37(2):237-241. ]

[13]

Zhang

, Shao J

, Luo

Spatial interpolation of air temperature with ANUSPLIN in Three Gorges Reservoir Area. International Conference on Remote Sensing

[C], Environment and Transportation Engineering, IEEE. 2011(147):3465-3468. DOI: 10.1109/RSETE.2011.5965057

[14]

Liu

, Mcvicar T

, Li L

, et al.

Modeling spatial distribution of pan evaporation based on quinat-variate thin plate spline function

[J]. Science of Soil and Water Conservation, 2006,4(6):23-30.

[15]

Wahba

Spline models for observational data

[J]. Mathematics of Computation. 1991. 195(57):444-445. DOI: 10.2307/2938687

[16]

徐翔, 许瑶, 孙青青, 等.

复杂山地环境下气候要素的空间插值方法比较研究

[J]. 华中师范大学学报(自然科学版), 2018,52(1):122-129.

[ Xu

, Xu

, Sun Q

, et al.

Comparision study on meteorological spatial interpolation approaches in Kangdian region of China

[J]. Journal of Huazhong Normal University(Natural Sciences), 2018,52(1):122-129. ]

[17]

易桂花, 张廷斌, 何奕萱, 等.

四种气温空间插值方法适用性分析

[J]. 成都理工大学学报(自然科学版), 2020,47(1):115-128.

[ Yi G

, Zhang T

, He Y

, et al.

Applicability analysis of four spatial interpolation methods for air temperature

[J]. Journal of Chengdu University of Technology(Social Sciences), 2020,47(1):115-128. ]

DOI:10.7522/j.issn.1000-0534.2015.00037 [本文引用: 1]

[18]

谭剑波, 李爱农, 雷光斌.

青藏高原东南缘气象要素Anusplin和Cokriging空间插值对比分析

[J]. 高原气象, 2016,35(4):875-886.

选取地形起伏度巨大的青藏高原东南缘为研究区,利用该研究区96个气象站点,结合高程数据,分别采用Cokriging和Anusplin空间插值方法,获取2010年250 m分辨率的年均温度和年累计降水插值曲面。并采用交叉验证方法对比Anusplin与Cokriging插值精度,分析了误差的空间分布特征,重点对比两种插值曲面差异较大的区域精度优劣,评价两种方法在复杂地区的适用性。结果表明,Anusplin在复杂地表的插值表现优于Cokriging,其中Anusplin气温插值的均方差仅为0.82℃,而Cokriging的均方差为1.45℃；两者的降水插值精度基本一致,但Anusplin在气象要素空间异质性大的区域优于Cokriging。因此,与Cokriging相比,Anusplin更适合青藏高原东南缘复杂地表气象要素空间插值。

[ Tan J

, Li A

, Lei G

Contrast on Anusplin and Cokriging meteorological spatial interpolation in southeastern margin of Qinghai-Xizang Plateau

[J]. Plateau Meteorology, 2016,35(4):875-886. ]

[19]

郑小波, 罗宇翔, 于飞, 等.

西南复杂山地农业气候要素空间插值方法比较

[J]. 中国农业气象, 2008,29(4):458-462.

[ Zheng X

, Luo Y

, Yu

, et al.

comparisons of spatial interpolation methods for agroclimate factors in complex mountain areas of southwest China

[J]. Chinese Journal of agrometeorology, 2008,29(4):458-462. ]

[20]

Liang

, Zeng

, Li

, et al.

Consistent comparison of remotely sensed sea ice concentration products with ERA-Interim reanalysis data in Polar Regions

[J]. Remote Sensing, 2020,12(18):2880. DOI: 10.3390/rs12182880

[21]

Johnston B

, Randel W

, Sjoberg J

Evaluation of tropospheric moisture characteristics among COSMIC-2, ERA5, and MERRA-2 in the tropics and subtropics

[J]. Remote Sensing, 2021,13(5):880. DOI: 10.3390/rs13050880

[22]

Wang

, Xu

, Nie

, et al.

Evaluation of precipitable water vapor from five reanalysis products with ground-Bbased GNSS observations

[J]. Remote Sensing, 2020,12(11):1817. DOI: 10.3390/rs12111817

[23]

Kumar

, Yadav

, Mohan

Seasonal sea-ice variability and its trend in the Weddell Sea sector of West Antarctica

[J]. Environmental Research Letters, 2021,16(2):24-46. DOI: 10.1088/1748-9326/abdc88

DOI:10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2 URL [本文引用: 1]

[24]

Kalnay

, Kanamitsu

, Kistler

, et al.

The NCEP/NCAR 40-year reanalysis project

[J]. Bulletin of American Meteorological Society, 1996,77(3):437-471.

[25]

Kistler

, Kalnay

, Collins

, et al.

The NCEP-NCAR 50-year reanalysis: Monthly means CD-ROM and documentation

[J]. Bulletin of the American Meteorological Society, 2001,82(2):247-268. DOI: 10.1175/1520-0477(2001)082

[26]

Kazutoshi

, Junichi

, Hiroshi

, et al.

The JRA-25 reanalysis

[J]. Journal of the Meteorological Society of Japan, 2007,85(3):369-432. DOI: 10.2151/jmsj.85.369

[27]

Hersbach

, Bell

, Berrisford

, et al.

The ERA5 global reanalysis

[J]. Quarterly Journal of the Royal Meteorological Society, 2020,146(730):1999-2049. DOI: 10.1002/qj.3803

[28]

Mearns L

, Gutowski

, Jones

, et al.

A regional climate change assessment program for North America

[J]. Eos Transactions American Geophysical Union, 2009,90(36):311-311. DOI: 10.1029/2009EO360002

[29]

高谦, 江志红, 李肇新.

多模式动力降尺度对中国中东部地区极端气温指数的模拟评估

[J]. 气象学报, 2017,75(6):917-933.

[ Gao

, Jiang Z

, Li Z

Simulation and elaluation of multi-model dynamical downscaling of temperature extreme indices over the Middle and East China

[J]. Acta Meteorologica Sinica, 2017,75(6):917-933. ]

DOI:10.7522/j.issn.1000-0534.2016.00039 [本文引用: 1]

[30]

肖林鸿, 高艳红, Chen

Fei

, 等.

青藏高原极端气温的动力降尺度模拟

[J]. 高原气象, 2016,35(3):574-589.

利用ERA-Interim再分析资料作为边界条件,基于耦合陆面模式Noah-MP的区域气候模式WRF在东亚区域进行了动力降尺度模拟(简称WRF2),对比格点观测资料,评估了动力降尺度对青藏高原极端气温指数的模拟能力,在此控制试验基础上,分别将WRF的陆面模式替换为Noah LSM,边界条件替换为CCSM4,进行了两组敏感性试验(分别是WRF1 和WRF3),通过与控制试验的比较,分析了边界条件和陆面模式对极端气温指数模拟的影响。结果表明,WRF2 能较好地模拟青藏高原极端气温指数气候态的空间分布,但存在一定的冷偏差;受边界条件影响WRF3 模拟的极端气温指数的气候倾向率存在负偏差。同时,尽管采用不同的边界条件,耦合相同陆面过程的两次数值试验对极端气温空间分布的模拟能力相似,相比WRF2,WRF1 表现出更强的冷偏差;但边界条件对极端气温指数气候倾向率的影响大于陆面模式,WRF3 模拟的极端气温指数气候倾向率与观测结果更为接近。

[ Xiao L

, Gao Y

, Chen

, et al.

Dynamic downscaling simulation of extreme temperature indices over the Qinghai-Xizang Plateau

[J]. Plateau Meteorology, 2016,35(3):574-589. ]

[31]

钞振华

中国西部地区气温资料的统计降尺度研究

[J]. 干旱区研究, 2011,28(5):879-879.

[ Chao Z

Study on Statistical Downscaling of Air Temperature in West China

[J]. Arid Zone Research, 2011,28(5):879-879. ]

[32]

Oshima

, Kato

, Kadokura

An application of statistical downscaling to estimate surface air temperature in Japan

[J]. Journal of Geophysical Research, 2002,107(D10):ACL-14.

[33]

胡德勇, 乔琨, 王兴玲

等.

利用单窗算法反演Landsat 8 TIRS数据地表温度

[J]. 武汉大学学报·信息科学版, 2017,42(7):869-876.

[ Hu D

, Qiao

, Wang X

, et al.

Comparion of three single-window algorithms for retrieving land-surface temperature with landsat 8 TIRS data

[J]. Geomatics and Information Science of Wuhan University, 2017,42(7):869-876. ]

[34]

徐涵秋

新型Landsat8卫星影像的反射率和地表温度反演

[J]. 地球物理学报, 2015,58(3):741-747.

[ Xu

Hanqiu

Retrieval of the reflectance and land surface temperature of newly-launched landsat 8 satellite

[J]. Chinese Journal of Geophysics, 2015,58(3):741-747. ]

[35]

Breiman

Random Forest

[J]. Machine Learning, 2001,45(1):5-32.

DOI:10.1023/A:1010933404324 URL [本文引用: 3]

[36]

Salomonson

, Appel

Estimating fractional snow cover from MODIS using the normalized difference snow index

[J]. Remote Sensing of Environment, 2004,89(3):351-360. DOI: 10.1016/j.rse.2003.10.016

[37]

Dirksen

, Knap W

, Steeneveld G

, et al.

Downscaling daily air-temperature measurements in the Netherlands

[J]. Theoretical and Applied Climatology, 2020,142:1-17. DOI: 10.1007/s00704-020-03313-1

[38]

Liang

Narrowband to Broadband Conversions of Land Surface Albedo I Algorithms

[J]. Remote Sensing of Environment, 2001,76(2):213-238. DOI: 10.1016/S0034-4257(00)00205-4

[39]

Mutiibwa

, Strachan

, Albright

Land Surface Temperature and Surface Air Temperature in Complex Terrain

[J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2016,8(10):4762-4774. DOI: 10.1109/JSTARS.2015.2468594

[40]

Phan

, Martin

, Jan

Estimating Daily Maximum and Minimum Land Air Surface Temperature Using MODIS Land Surface Temperature Data and Ground Truth Data in Northern Vietnam

[J]. Remote Sensing, 2016,8(12):1002-1025. DOI: 10.3390/rs8121002

[41]

Minder J

, Mote P

, Lundquist J

Surface temperature lapse rates over complex terrain: Lessons from the Cascade Mountains

[J]. Journal of Geophysical Research Atmospheres, 2010,115(D14):333-351. DOI: 10.1029/2009jd013493

[42]

Weng D

, Sun Z

A preliminary study of the lapse rate of surface air temperature over mountainous regions of china

[J]. Geographical Research, 1984,3(2):24-34.

[43]

Mahrt

Variation of surface air temperature in complex terrain

[J]. J Appl Meteorol Climatol. 2006. 45(11):1481-1493. DOI: 10.1175/JAM2419.1