1. 引言
青藏高原(简称“高原”)被称为地球的“第三级”,是世界上面积最大、海拔最高的高原,其特殊的地形和下垫面使得高原地–气相互作用过程对东亚大气环流甚至全球变化都有重要影响,被认为是我国气候变化的启动区[1]-[4]。地表温度是陆面研究的重要参数,表征着地面热源特征,其值准确与否直接关系到地–气相互作用过程有关问题研究结果的正确性,在高原气候变化研究中起着十分重要作用[5]-[8]。因此,高原地表温度的时空变化特征对于进一步研究高原地–气相互作用具有重要意义。
目前,高原地表温度数据主要有三个来源:气象观测站现场测量、卫星产品和再分析数据[9]-[11]。气象观测站可以直接测量地表温度,但地域连续性较差,不能够较好地表现出整个高原的地表温度分布状况。大气再分析资料可以提供空间覆盖完整、时间均一稳定的长序列历史天气资料,但由于再分析资料在数据同化过程中使用的卫星资料在地形起伏较大的山区具有明显的局限性[12],同时,再分析资料作为非独立观测资料,不能完全取代观测资料来描述大气的真实三维状态、客观反映气候的长期变化趋势。因此在使用前需对其在高原地区的可信程度、适用范围进行验证,从而提高气候变化研究的可靠性和再分析资料的质量[13]-[15]。
许多学者已经对多种再分析资料在高原地区的适用性进行了评估,例如:除多等[16]利用高原60个气象站点地面观测气温数据,表明MERRA-2再分析产品相比其他全球再分析资料在高原地面气温表达方面具有一定的优势。刘婷婷等[17]利用全国站点日温度数据评价ERA5再分析地面气温数据在中国区域的适用性,发现ERA5数据低海拔区域较高海拔区域精度高,在高原气候区精度最低。胡梦玲等[18]发现JRA-55能很好地反映出位势高度的变化趋势,与位势高度资料和探空资料相关性最好,但在气候变化趋势方面存在时空差异。赵洪宇等[19]研究发现JRA-55再分析水汽资料在模式及观测资料同化方面均有诸多改进,但与高原观测资料的偏差较大。基于以上研究,我们发现ERA5、JRA-55、MERRA-2等新一代大气再分析资料在高原地区的适用性研究主要针对气温、风场、水汽等要素,但针对地表温度要素的适用性对比评估分析的研究较少。
影响青藏高原地区地表温度的因子有海拔、归一化植被指数(NDVI)和积雪覆盖率等,各种因素相互作用相互联系,共同影响地表温度。海拔、NDVI以及积雪对地表温度的影响均呈现负相关关系[20]-[24]。其中,海拔是地形因子中影响LST的主要因素,地表温度随海拔升高呈自然对数形式降低[20]。同时海拔还会影响积雪的分布,青藏高原积雪分布以高海拔特征为主,有明显的垂直地带性[25]。植被的分布同时受积雪及海拔影响,同时随着海拔的升高,植被与积雪之间转化的可能性随之增加[23] [26] [27]。因此,青藏高原复杂的地形地势地貌、高海拔和冬季的积雪对于地表温度有重要影响。
基于以上论述,本研究依据高原地区2001~2016年站点观测数据,定量评估ERA5、JRA-55和MERRA-2的三种再分析地表温度数据,分析三种再分析资料在年、季尺度上高原地表温度的时间和空间特征,并且分析海拔、归一化植被指数(NDVI)、积雪覆盖率三个因子对高原地表温度的影响。
2. 数据来源与方法介绍
2.1. 站点介绍
青藏高原(Qinghai-Xizang Platea, TP)位于欧亚大陆的中心,范围为24˚~40˚N,67˚~105˚E。本文所用观测资料为2001~2016年青藏高原地区130个基准气象台站(图1)的逐月平均地表温度数据。
Figure 1. Distribution of meteorological station in the Qinghai-Xizang Platea
图1. 青藏高原观测气象站点空间分布
2.2. 数据来源和仪器介绍
2.2.1. 再分析资料
ERA5是欧洲中期天气预报中心生产的全球气候的第五代大气再分析数据集。其时空分辨率大幅提升,水平分辨率为0.25˚ × 0.25˚,产品要素包括大气、海洋和陆地等共240种[28],本文选取2001~2016年ERA5月尺度的地表温度数据(https://cds.climate.copernicus.eu/)。JRA-55是日本气象厅的第2套全球大气再分析资料,采用更先进的数据同化方案(4Dvar与3Dvar)以及新的辐射方案和变分偏差校正。本文选取2001~2016年JRA-55月尺度的地表温度数据(https://rda.ucar.edu/)。MERRA-2再分析数据提供了从1980年开始的数据,空间分辨率为0.5˚ × 0.625˚,取代了原来的MERRA数据集。与MERRA数据集相比MERRA-2同化系统取得了进步,能够同化现代高光谱辐射和微波观测,以及GPS-无线电掩星数据集。本文选取2001~2016年MERRA-2月尺度的地表温度数据(https://rda.ucar.edu/)。
2.2.2. 卫星遥感产品
本研究选用的归一化差异植被指数(NDVI)和积雪覆盖率(Snow Cover)来自在Terra航天器上运行的中分辨率成像分光仪(MODIS)仪器,能够提供空间分辨率为0.05˚的MODIS气候模拟网格单元(CMG)的月平均数据,具有较高的空间分辨率和光谱分辨率,并且有良好的动态变化观测能力[21] [29]-[31]。NDVI 作为表征植被生长状况的指示因子,用来反映绿色植被的相对丰度和发育状况;Snow Cover产品提供每个CMG的月平均积雪覆盖百分比[32]。这两个数据产品分别来自于美国地质调查局EROS数据中心(EDC) (https://lpdaac.usgs.gov/)的陆地过程DAAC的MOD13C2产品和国家冰雪数据中心(NSIDC) (https://nsidc.org/home)的MOD10CM产品。本研究使用了如下MODIS数据产品(V061):
1) MODIS/Terra Vegetation Indices Monthly L3 Global 0.05Deg CMG V061;
2) MODIS/Terra Snow Cover Monthly L3 Global 0.05Deg CMG V061。
2.2.3. 统计指标
本文选取偏差(bias)、均方根误差(RMSE)和相关系数(r)对ERA5、JRA-55和MERRA-2地表温度数据进行定量评估,以上统计指标的具体计算公式如下:
偏差:
(1)
均方根误差:
(2)
相关系数:
(3)
其中,
和
分别为再分析地表温度和站点观测地表温度。偏差(公式(1))表示再分析与站点观测地表温度之间的差异;均方根误差(公式(2))是再分析与站点观测地表温度之间偏差的度量;相关系数(公式(3))反映了再分析与站点观测地表温度之间的线性相关程度[33]。
站点数据与网格数据在空间尺度上的不匹配是需要解决的一个问题,本研究采用了最近邻插值法,最邻近插值方法是在经纬度基础上从网格点中提取距离站点最近的网格数据作为缺失的站点数据,而且相较于其他插值方法该种方法产生的误差较小[34] [35] [36]。
2.2.4. 路径分析法
本研究选取的海拔(X1)、NDVI (X2)和积雪覆盖率(X3)三个因子不仅影响地表温度(Y),且三个因子之间也有一定的相互影响,为了更准确地评估三个因子对地表温度的影响,本研究使用路径分析方法。路径分析是一种多元统计分析方法,常用于研究自变量与因变量之间的相关关系,它将因子对因变量的影响分解为直接和间接两部分[37] [38]。使用路径系数表征变量之间因果联系的强度,并用标准化回归系数进行分析[39] [40]。本文使用SPSS Statistics 26.0 (IBM Corporation, Armonk, NY, USA)计算路径系数和显著性检验进行路径分析。除Bonferroni校正方法外,低于0.05的P值被视为具有统计学意义[40]。
本研究将路径分析中的各个变量分为自变量、中介变量、应变量和残差变量。根据以往研究者的研究结果,海拔是影响地表温度的主要地形因子,且地表温度随海拔升高呈自然对数形势降低[20],显然海拔是不受其他因素影响的因子,因此本研究中将海拔作为自变量,即只能为因的变量;中介变量为既是因又是果的变量,而NDVI和积雪覆盖率对地表温度都存在负相关关系[21] [23],海拔对NDVI空间分布的影响作用显著[41],积雪覆盖率以高海拔为主要分布区域,海拔越高,积雪覆盖率越大[25] [42],本研究中将NDVI和积雪覆盖率作为中介变量;应变量为只能是果的变量,即本研究中的地表温度;残差变量是因果模型之外的影响因变量的所有变量的总称,即图2中的e1、e2和e3。
Figure 2. A model of multivariable path analysis
图2. 路径分析模型
本研究的路径分析模型如图2所示,其中P12、P13、P1Y、P2Y和P3Y为路径系数。同时,由于NDVI与积雪覆盖率之间存在显著的负相关性[43],本研究考虑了NDVI和积雪覆盖率两个中介变量之间的相关性,即图2中的r23。此外,由于积雪覆盖率因子具有强烈的季节性响应,本研究仅在冬季分析了积雪覆盖率因子的影响,而在全年范围内仅分析海拔和NDVI对地表温度的影响。路径分析的结构方程组为:
(4)
3. 结果
3.1. 多种大气再分析地表温度定量评估
分析2001~2016年观测资料与ERA5、JRA-55、MERRA-2再分析资料的年、月平均地表温度变化。由图3可知,三种再分析资料地表温度年际变化与观测值基本一致,但显著低于观测值,整体呈现增温趋势。地表温度月际变化基本一致,1月平均地表温度最低,2月以后开始升高,到7月达到峰值,8月之后开始下降,均表现为夏季高、冬季低的单峰型特征。
Figure 3. Surface temperature time series of ERA5, JRA-55 and MERRA 2 observed surface temperatures from 2001 to 2016 in the Qinghai-Xizang Platea
图3. 2001~2016年青藏高原地表温度时间序列
首先比较ERA5、JRA-55和MERRA-2三种再分析资料与观测平均地表温度偏差的空间分布(图4)。三种再分析平均地表温度与实际观测地表温度相比均明显偏低,所有站点的平均偏差均为负值,对实际地表温度过冷估计,在空间尺度上均呈现南北差异,具有边缘高,中部低的分布特征具有边缘高,中部低的分布特征。ERA5在高原东南部、东北祁连山地区偏差较大,JRA-55和MERRA-2在青藏高原南部、柴达木盆地地区偏差较大,MERRA-2全年偏差变化不明显。三种再分析资料均存在季节变化。ERA5在春季和冬季偏差大,而JRA-55在秋季和冬季偏差较小,MERRA-2整体偏差较小,对地表温度的模拟能力较强,能较好的表现出青藏高原地表温度的空间分布特征。
Figure 4. Deviation spatial distribution map of ERA5, JRA-55 and MERRA-2 observed surface temperatures from 2001 to 2016 in the Qinghai-Xizang Platea (unit: ˚C)
图4. ERA-5、JRA-55和MERRA-2年、季平均地表温度偏差的空间分布(单位:℃)
进一步考察ERA5、JRA-55和MERRA-2三种再分析资料对青藏高原地表温度变化情况(图5)。ERA5数据与观测资料相比,整体呈现南高北低的分布态势,分布不均,高原南部地表温度均方差较大存在高值中心,高原北部和边缘均方差较小,季节尺度来看,在春季、冬季有明显的高值中心为22℃~27℃,其余地区均方差达到8℃~16℃。JRA-55数据在春季、夏季,高原北部和南部的均方差较大,高原中部较低,在秋季、冬季,高原平均地表温度呈现南高北低的分布态势,整体均方差均在9℃~13℃。MERRA-2数据与观测资料相比,逐月平均地表温度均方差整体较小,均方差均在4℃~12℃,没有明显的高值中心,四季变化不明显。
Figure 5. Mean square error spatial distribution map of ERA5, JRA-55 and MERRA-2 observed surface temperatures from 2001 to 2016 in the Qinghai-Xizang Platea (unit: ˚C)
图5. ERA-5、JRA-55和MERRA-2年、季平均地表温度均方差的空间分布(单位:℃)
Figure 6. Correlation coefficient spatial distribution map of ERA5, JRA-55 and MERRA-2 observed surface temperatures from 2001 to 2016 in the Tibetan Platea
图6. ERA-5、JRA-55和MERRA-2年、季平均地表温度相关分布
考察ERA5、JRA-55和MERRA-2三种再分析资料与观测平均地表温度的相关情况(图6)。MERRA-2地表温度与观测地表温度呈现正相关,最高相关系数达到0.9,冬季相关系数最大,相关性不显著的地区集中在高原西南部、南部、中部的局部地区,JRA-55和ERA5与观测地表温度的相关分布较相似,在东南部、东部边缘地区、柴达木盆地地区相关性最显著,均在夏季相关性最显著,秋季最弱,ERA5与高原大部分地区相关性最差甚至出现负相关的情况。对比而言,MERRA-2地表温度与观测平均地表温度的相关性最好。其次是JRA-55,而ERA-5与观测平均地表温度的相关性最差,仅在高原东南部和北部相关性显著。
逐月计算ERA-5、JRA-55和MERRA-2三种再分析地表温度平均值以考察它们的时间变化特征(图7)。从偏差看,ERA-59月平均地表温度偏差最小,3月最大;JRA-55月平均地表温度偏差以1、11、12月最小,3、8月最大;MERRA-2月平均地表温度偏差8月最大,1、12月最小。从均方差看,ERA-5月平均地表温度均方差以2、3月最大,8、9月最小;JRA-55和MERRA-2月平均地表温度均方差均在1、11、12月最小。各月三种再分析平均地表温度均为正相关,MERRA-2与观测月平均地表温度的相关性最好,4月、7月的相关系数最高,超过了0.5;JRA-55以6月相关系数最大,12月最低;ERA5各月平均地表温度与观测平均地表温度相关系数均低于相应月份MERRA-2和JRA-55,相关性最好的是6月,也仅有0.4。
Figure 7. Surface temperature time series of ERA5, JRA-55 and MERRA-2 observed surface temperatures from 2001 to 2016 in the Qinghai-Xizang Platea (unit: ˚C)
图7. ERA-5、JRA-55和MERRA-2年、月平均地表温度时间变化序列(单位:℃)
3.2. 多种大气再分析地表温度区域变化特征分析
从前2种模态的解释方差(表1)来看,观测资料的解释方差与再分析资料相差较大,分别为95.46%和1.48%。三种再分析资料第一模态的解释方差在30%~40%之间,其中JRA-55资料最大(38.61%),ERA5资料最小(30.61%);第二模态的解释方差在12%~21%之间,JRA-55最大(20.89%),ERA5最小(12.90%)。
Table 1. Variance contribution of EOF analysis in the first two modes of EOF data for various land temperature data (%)
表1. 多种地表温度资料EOF分析前两个模态的方差贡献(%)
|
OBS |
MERRA-2 |
ERA5 |
JRA-55 |
第一模态 |
95.46 |
37.79 |
30.61 |
38.61 |
第二模态 |
1.48 |
19.50 |
12.90 |
20.89 |
累计方差贡献 |
96.94 |
57.29 |
43.51 |
59.5 |
3.2.1. 站点资料地表温度区域变化特征分析
本文对青藏高原站点资料的年平均地表温度经过标准化处理后做经验正交分解(EOF) (图8(a1)~(d1))。第一模态的空间分布显示出出高原主体为一致的负值区,表现为高原地表温度变化的整体一致性,结合时间系数可以看出青藏高原地表温度呈现逐年上升的趋势,但升温幅度不大。从第二模态的空间分布中,可以看出高原从中部到东部呈现出“-+”的变化特征,高原中部负值的绝对值略大于东部正值,结合时间系数得出2012年青藏高原中部和东部地表温度的反向变化最明显。
3.2.2. 三种再分析资料地表温度区域变化特征分析
对青藏高原三种再分析资料的年平均地表温度经过标准化处理后做经验正交分解(图8(a2)~(d4))。MERRA-2第一模态的空间分布显示出高原大部分区域为一致的负值区,只有东南部和西北部有小部分正值区,结合时间系数可以看出仅有东北部与西南部小部分区域显示出与观测资料相一致的地表温度升高的趋势,2012年正值区与负值区温度反向变化最明显。分析MERRA-2第二模态的空间分布和时间系数可以得出,青藏高原除西南小部分区域外其余地区呈现地表温度上升的趋势,这与观测资料的分析结果在大部分区域中保持一致。
ERA5第一模态的空间分布和时间系数以及JRA-55第一模态的空间分布和时间系数均没有显示出明显的温度变化趋势。ERA5第二模态的空间分布结合时间系数显示出西部降温东部升温的趋势;JRA-55第二模态的空间分布结合时间系数显示出JRA-55模拟的地表温度仅在2001~2006年中大部分区域呈现升温趋势,其余年份趋势并不明显。
Figure 8. EOF analysis of various surface temperature data
图8. 多种地表温度资料EOF分析
通过站点资料EOF分析,我们得出青藏高原全域地表温度均呈现出随时间升高的趋势。再综合MERRA-2、ERA5和JRA-55的EOF分析与站点资料EOF分析的比较情况,三种再分析资料中仅有MERRA-2的第二模态结果在大部分区域里(除西南小部分区域外)与站点资料得出的结果相一致,进一步说明MERRA-2在模拟青藏高原地表温度及其年际变化趋势方面的适用性比ERA5和JRA-55的适用性好。
3.3. 多种大气再分析地表温度与气象因子的关系
3.3.1. 各变量的相关关系
所有变量的双变量相关系数(r)见表2。双变量相关矩阵表示地表温度与三个因子都有显著相关,其中与海拔因子的相关性最大(r = −0.363),与NDVI呈正相关关系,与积雪覆盖率呈负相关。在三个因子之间的相关关系中,海拔与NDVI的相关性最大(r = 0.478),而积雪覆盖率与海拔和NDVI的相关性均呈负相关关系。
Table 2. Correlation (r) matrix of the variables
表2. 各变量之间的相关系数(r)
|
X1 |
X2 |
X3 |
Y |
X1海拔 |
1.000 |
|
|
|
X2 NDVI |
0.478** |
1.000 |
|
|
X3积雪覆盖率 |
−0.105** |
−0.155** |
1.000 |
|
Y地表温度 |
−0.363** |
0.135** |
−0.180** |
1.000 |
**P < 0.01.
3.3.2. 地表温度影响因子的路径分析
通过路径分析,表3列出了三个因子分别在全年和冬季尺度下对站点观测和再分析地表温度总效果系数,其中LST-O,LST-E,LST-J和LST-M分别表示站点观测、ERA5、JRA-55和MERRA-2的地表温度。可以看出,海拔对LST的影响效果最为显著。对所有资料地表温度来说,在全年范围内,海拔对地表温度有负影响,而NDVI对地表温度有正影响效果;但是在冬季,由于积雪因子的引入,NDVI对地表温度的影响变为负影响,海拔对地表温度的负影响也有所加强。考虑到站点观测地表温度相对于再分析资料更为准确,在后面具体讨论再分析地表温度对三个因子的响应时,以三个因子对站点观测地表温度的影响程度作为参考真值进行比较分析。
Table 3. The total effect coefficients during winter and all year
表3. 冬季和全年尺度下的总效果系数
|
LST-O |
LST-E |
LST-J |
LST-M |
X1海拔 |
冬季 |
−0.584 |
−0.485 |
−0.622 |
−0.541 |
全年 |
−0.394 |
−0.230 |
−0.426 |
−0.291 |
X2 NDVI |
冬季 |
−0.117 |
−0.140 |
−0.246 |
−0.160 |
全年 |
0.134 |
0.146 |
0.075 |
0.164 |
X3积雪覆盖率 |
冬季 |
−0.180 |
−0.073 |
−0.114 |
−0.153 |
**P < 0.01.
由表3可知海拔对地表温度有负影响效应,并且随着海拔的升高地表温度变冷。在全年尺度上,海拔对站点观测地表温度(LST-O)的总效果系数为−0.394;海拔对ERA5产品(LST-E)和MERRA-2产品(LST-M)的总效果系数略小于海拔对LST-O的影响,但对LST-J的影响略大。在冬季,海拔对地表温度的负影响增大到−0.584,海拔对再分析地表温度与对LST-O的总效果系数的相对大小与全年尺度相比保持不变。在这两个时间尺度上,LST-J对海拔因子的响应程度与LST-O最接近,LST-M次之且相差不大。
同时可得到,NDVI对地表温度的影响在全年尺度下为正效应,即植被长势越好,地表温度越高;而在冬季正相反,NDVI对地表温度的影响为负效应。全年范围内NDVI对LST-O的总效果系数为0.134,在冬季为−0.117。综合两个时间尺度,相比之下LST-E对NDVI因子的响应效果最接近于LST-O,其次是LST-M。
本文仅在冬季尺度上研究了积雪覆盖率因子对LST的影响,积雪覆盖率对地表温度呈负影响效应,这意味着积雪覆盖率越高,地表温度越冷。积雪覆盖率对LST-O的总效果系数为−0.180,对LST-M的影响(−0.153)与LST-O最接近。
4. 讨论
本文通过ERA5、JRA-55和MERRA-2三种再分析地表温度模拟资料,分析了近10年来青藏高原地表温度变化趋势,并对地表温度模拟资料的适用性进行了评估,高原变化趋势与朱智等[44]的研究结果较为一致,本文认为MERRA-2在高原的适用性较好,但JRA-55、MERRA-2再分析资料存在较大偏差的原因及其在青藏高原其他区域适用性还有待进一步检验。
本文通过站点资料EOF分析得出的青藏高原地表温度变化趋势与朱智等[44] [45]研究得到的青藏高原地表温度随年份呈上升趋势结论相一致,但是三种再分析资料EOF分析得出的结果均与站点资料相比较均产生了不同程度的差异,这可能与本文选取研究的年份较少对于趋势的呈现不够明显以及不同再分析资料对于植被等地表温度影响因素的同化情况有关,也可能与青藏高原的观测站点的分布情况不均匀有关。
本文中通过路径分析方法得到NDVI与地表温度呈正相关,仅在冬季引入积雪覆盖率因子后NDVI的影响为负效应,但是有研究表示NDVI与地表温度呈负相关关系[21],这可能与积雪和地表温度之间有显著的反馈作用有关,而本文在研究三个因子对LST的影响时,在全年范围内并没有使用积雪覆盖率因子,忽略了积雪对LST的影响,导致全年范围内结果出现偏差。以上导致研究结果不一致的归因还需要进一步具体研究说明。
5. 结论
本研究利用2001~2016年站点观测资料评估三套地表温度再分析产品(ERA5、JRA-55和MERRA-2)在高原的适用性,并讨论了三个气象因子(海拔、NDVI和积雪覆盖率)对高原地表温度的影响程度,得到以下主要结论:
1) 2001~2016年青藏高原地表温度观测数据,三种地表温度模拟资料均与实际观测值出现不同程度的低估现象,呈现南北差异,具有季节性变化特征,相比于其他两种地表温度模拟资料,MERRA-2再分析地表温度资料与实际观测资料较吻合。
2) 2001~2016年青藏高原全域地表温度均呈现出随时间升高的趋势,三种再分析资料中仅有MERRA-2第二模态分析得出的趋势与青藏高原实际情况在大部分区域上(除西南部分区域)相接近,在模拟青藏高原地表温度时空变化情况上的表现优于其他两种资料。
3) 三个因子对站点观测地表温度的影响程度与它们之间的相关性差别不大。海拔和积雪覆盖率因子对地表温度为负影响效果,而NDVI因子对地表温度的影响随季节变化,冬季为负影响,全年平均为正影响。三个再分析产品的地表温度对这三个因子的响应与站点观测地表温度相差不大,综合来看,MERRA-2地表温度的响应效果优于另外两个再分析资料,更接近于作为参考真值的站点观测地表温度。
基金项目
成都信息工程大学大学生创新创业训练项目(校级)202310621021资助。