地表温度(land surface temperature,LST)决定了陆地表层系统辐射能量平衡中的长波辐射,是地表辐射能量平衡和水循环物理过程中的关键参数^[1,2,3], 因此,准确、及时获取LST值,掌握其时空变化规律对农情信息监测^[4]、生态和环境研究^[5]以及气候变化^[6]等方面具有重要意义。

传统获取LST值的方法是通过气象站点或地面测量获得,在实际应用中,各种插值方法常用于将单点数据扩展获取到空间连续的面数据。然而海拔、纬度和下垫面性质等影响因子对LST值的影响时空变异性较大^[7]。插值后获得的空间栅格化的LST数据可靠性不高,因为大多数地区的地面站密度都远远不够。基于卫星热红外技术反演LST值的方法,较传统的地面定点观测具有能够获取大尺度、空间连续的LST值和近实时数据访问的优势,已成为获取LST值可靠的数据源^[8,9,10]。迄今为止, 已有许多的研究者利用热红外数据研究中小尺度水平城市地区LST时空变化的研究成果^[11,12,13], 但对于长时间序列,大范围,全局地区展开宏观研究较少。因此,本文基于重建的2003—2017年间无云条件下的MODIS LST数据,较为精准地获取了近15 a来中国LST值的时空分布格局及其变异规律,并对驱动因素进行了探析,以期对农作物生产、农业干旱监测及相关气候变化等的研究提供参考。

本研究选用了2种v006版本的 MODIS LST产品数据集: 逐日LST数据(MOD11C1和MYD11C1)和逐月LST数据(MOD11C3和MYD11C3)用于LST的变化特征研究,分别包括了昼夜2次的观测数据。空间分辨率均为0.05°。时间跨度为2003年1月—2017年12月,来自于美国国家航空航天局网站。原始的MODIS LST文件中还包含了相同时空分辨率的质量控制(quality control,QC)图层,为每个像元提供了可靠的质量控制信息。使用MODIS Reprojection Tool (MRT)工具为Terra和Aqua提取了昼夜2景LST数据及相应的QC数据,并进行了重投影、地理子集裁剪、格式和单位转换等预处理工作,最后得到tiff格式的Albers投影的LST数据。

地面观测站点的LST数据来源于中国气象局,涵盖了2003—2017年间中国2 399个气象地面站的LST记录。为确保数据重建和精度验证过程的独立性,卫星过境时刻对应的气象站点数据被随机分为2个完全独立的子集: ①地面站点总数的80%(1 919个站点)用于数据重建; ②20%(480个站点)的数据作为验证集用于重建结果的精度验证。其中,从子集①中提取出可用的地面站点的LST用于4个时刻的MODIS数据重建。子集②的月平均地面站点数据则用于LST数据的验证。所有站点数据均经过质量检查,对数据中的缺测数据和异常值进行了剔除。

LST重建过程包括基于QC低质量像元值的去除和LST值的重建2个步骤。为了避免不可靠像元对影像质量及其对重建结果的影响,必须进行低质量像元值的去除。此步骤是将LST数据的QC作为掩模文件将低质量像元设为空值,只保留高质量的LST像元。最后,QC中平均LST误差≤1 K、LST数据产生质量好和平均发射率误差≤0.01的像元被视为高质量数据,其余像元值被拒绝。本文使用QC去除受云影响的像元而非使用影像的直方图阈值,这是为了确保尽可能多的低质量的值被去除,否则采用较低的LST误差阈值也会去除掉一些可靠的LST值^[14]。经过QC数据掩模后,冬季(以1月份为例),Terra和Aqua不可用像元的覆盖率分别为23.45%和24.31%; 夏季(以7月份为例),Terra和Aqua不可用像元的覆盖率分别为19.68%和20.86%。

平均而言,约有65%的地球表面被云层覆盖 ^[15,16]。MODIS LST数据是在晴空条件下通过热红外波段反演的,由于受到云和其他大气干扰无法获取云覆盖下真实的LST值。同时,云会对到达地表的太阳辐射产生衰减作用,造成云覆盖下LST值往往低于同等条件未受云影响的晴空下的LST值^[17]。当前的研究大多基于时空信息及各种辅助数据恢复晴空下的LST,重建云覆盖下的LST值一直是LST估算研究中的难题,重建具有大面积缺失值的逐日LST数据更加困难。而月合成的LST产品,数据时空完整性和一致性得到很大的改善,很大程度上减少了在重建过程中引入的不确定性误差。因此在本研究中,MODIS逐月LST数据用于对LST时空变化情况的研究,并且所提出的重建模型基于逐日LST数据(MOD11C1/ MYD11C1),该模型有效地保留了原始逐日LST数据中的高精度像元,仅重建了缺失的和误差大的像元,最终将逐日的像元值合成平均填充到逐月LST数据中的空值像元。这种方法能够从逐日数据中继承高质量数据,并保证重建结果精度的稳定性,数据重建流程如图1所示。

由图1可知,重建过程可分为以下3步: ①确定逐月LST数据中空值像元位置,遍历该像元对应的所有逐日LST像元,对逐日LST数据中的空值像元重建; ②由于地面站点的LST数据为云覆盖下LST的可靠数据源,所以为了仅重建云覆盖下的LST,使用经纬度作为控制条件将自地面气象站的观测值填充相同时刻相同位置的缺失像元,并将其标记用于下一步骤中; ③利用临近非空像元采用反距离加权回归法重建缺失像元。为了避免参与计算的像元数据冗余降低执行效率,同时确保参与计算的像元值数量,经过多次试验后,选择空值像元邻域内39像元×39像元的滑动窗口。同时,由于经过地面站点赋值后的可用像元在代表云覆盖下真实的LST时更有代表性,对站点数据进行了相对多倍权重赋值,具体公式为:

式中: D为非空像元到空值像元t的距离; (x,y)和(x_t,y_t)分别为非空像元和空值像元t的位置; i和j为用于估算空值像元t的非空像元,i为未受云影响的非空像元,j为由地面站点LST赋值后的像元; W_i和W_j分别为像元i和j对目标像元的权重; D_i和D_j分别为像元i和j到空值像元t的距离; m为像元i的数量; n为像元j数量; M_c和M_g分别为非空像元i和地面站点LST赋值像元j的相对权重系数,为了更好地突出地面站点数据作为非空像元对云覆盖下像元值的代表性,本文中M_c和M_g分别设置为1和3。

为了验证MODIS LST数据重建后的精度,使用验证数据集②对重建结果进行验证,数据验证的精度使用决定系数R²、平均绝对误差MAE和均方根误差RMSE作为度量指标,计算公式分别为:

式中: N为样本个数; LST_i和GST_i分别为MODIS和地面站观测的LST值; LST¯和 GST¯分别为MODIS和地面站点的观测均值。

为了量化LST在时间和空间的变化模式,本文采用最小二乘线性回归的方法计算每个像元的时间变化模式,进而反映2003—2017年间中国LST的时空间变化模式。采用斜率k表示近15 a来LST的变化速率,相关系数R表示每个像元LST值和时间序列的相关性,即

式中: u为年序号; T_u为第u年对应的LST值; v为研究时段,v=15。当k为正时,表示LST随时间序列呈增加趋势; 当k为负时,表示LST随时间序列呈递减少趋势,k的绝对值越大表征变化趋势越明显。R的阈值范围为[-1,1],其中R值为正时,代表LST与时间序列呈正相关关系; R值为负时,代表LST与时间序列呈负相关关系,R的绝对值越大代表随时间序列变化的相关性越强。

本文利用地面站点的实测数据进行精度验证分析。为了提高重建结果的精度,应该尽可能多地利用地面站点信息,使得地面产品与卫星数据在验证时具有更高的一致性,这也是目前大尺度数据精度验证最好的方法。同时,为了进一步增加站点数据的代表性,在地表类型变化大、地形起伏大的地方未来需要进一步增加站点数量提高精度。

对未参与数据重建过程的地面站点数据集提取了其位置对应的MODIS LST值进行重建结果的对比验证,如图2所示。由于站点随机分布于研究区的不同位置代表了不同的下垫面情况和气候类型,因此可以用来表示重建结果的可靠性。重建后的MODIS LST数据与地面站点数据拟合度良好,线性拟合R²为0.92,MAE为1.36 ℃,RMSE为1.51 ℃。

由中国2003—2017年间年平均LST值分布(图3(a))可知,除青藏高原外,LST的空间分布呈现大致与纬度带平行由南向北逐渐降低的分布格局,这主要是受太阳辐射量递减的影响。而青藏高原地区LST显著低于同纬度带的东部地区,这主要受海拔高度的影响。年平均LST在不同地区差异显著,最大相差约36 ℃。高值中心主要位于华南地区南部海南岛附近,LST达到25 ℃,低值中心并未出现在高纬度地区,而是位于青藏高原的部分区域,年平均LST低于-11 ℃。图3(b)为2003—2017年间分别以年际和季节为尺度的LST计算的距平指数时间序列。对LST的时序变化规律进行分析,近15 a间LST距平曲线总体呈现高低温交错出现的阶段性变化特点。LST最低值出现在2012年,年平均LST仅为8.41 ℃,距平值达到-0.86 ℃。在进一步对季节和月度尺度分析LST之后发现2012年的降温主要发生在冬季,集中在1—2月份,并且南方降温非常显著,同时由图3(b)可知在2011年冬季也呈现出显著的偏冷现象,导致这种现象的原因是这2 a正值拉尼娜事件的发生,由于东亚冬季季风异常强,冬季南部出现异常降雨。LST最高温度的年份为2007年,距平值为0.61 ℃。而最冷的年份和最暖的年份最大的距平绝对值均出现在冬季,同时不同年份冬季LST距平值波动特征最为显著,对年平均LST的变冷/变暖的贡献最大。

为了更详细地了解不同区域LST值变化的空间格局和区域差异,计算了2003—2017年间每个像元的LST年平均变化率和相关系数(图4)。2003—2017年间,LST在不同地区显示出明显的地域差异。中国大部分地区呈变暖趋势(k>0 ℃·a^-1),占整个研究区的63.7%。另外,大约20.80%的区域经历了较为显著的变暖(k>0.05 ℃·a^-1,R>0.6),LST的变暖速率总体呈现北方大于南方,西部大于东部的分布特点。

在东北地区,15 a间表现出明显的降温趋势,对中国降温的贡献最大。除了大兴安岭东北部以外,在东北大多数地区都观察到了广泛降温(k在-0.12~-0.06 ℃·a^-1之间,R>0.55)。Sun等^[18] 也发现了气温降低的趋势。作为中国重要的粮食产区,东北地区的降温将对农业生产和区域经济的发展造成不利影响。寒冷的加剧将导致作物生长期积温不足,进一步造成大规模的农作物减产。

在华北地区,LST值的空间变化模式呈沿太行山两侧相反的变化趋势: 升温区主要位于太行山东部低海拔的黄淮海平原附近,而太行山脉西部的黄土高原呈现降温趋势。黄淮海平原及其南部的长江三角洲人口密集,都表现出了较明显的升温趋势。相对较高的变暖趋势可能归因于近年来人口的快速增加以及在快速城市化过程中排放了更多的温室气体。

在西北地区,天山山脉的大部分地区和内蒙古高原及其南部地区显著变暖(k>0.10 ℃·a^-1,R>0.71),是全国变暖最强的区域。15 a中,快速升温对西北地区有重大影响。已有研究表明,近年来西北地区呈现变湿的态势^[19,20,21]。西北地区降雨稀少,在大规模缺水地区,气候变暖变湿对当地水资源和生态环境的影响更为显著。在干旱环境下,植被对气候的正向变化极为敏感,变暖变湿促使西北地区的植被生长产生显著的正反馈。在青藏高原地区,生态环境复杂,高原地区独特的地形和周围地区的热力特性在调节周围大气环流系统中起着重要作用,也是全球气候变化的关键区域。如图4所示,在青藏地区中部唐古拉山附近及南部地区呈现明显的变暖趋势(k>0.08 ℃·a^-1,R>0.5),柴达木盆地也表现出显著的升温(k>0.1 ℃·a^-1,R>0.6)。

进一步分析年际平均LST的昼夜变化情况(图5)。由图5(a)可知,2003—2017年间,几乎所有地区的年度白天LST增温/降温趋势均显著高于夜间(k在-0.03~0.03 ℃·a^-1之间),因此,平均LST升温/降温趋势可归因于白天的变化。而人类活动、太阳辐射强度及日照时间的变化等是白天升温的主导因素。白天,青藏高原的南部是升温最大的区域(k>0.09 ℃·a^-1,R>0.6)。Duan等^[22]发现自1998年以来,青藏高原南部的白天云量迅速减少导致日照时间增加。白天地表接收的太阳辐射的增加对升温效果显著,是导致白天LST升高的重要因素。与白天相比,夜间的LST变化趋势平缓,可以认为是稳定的。

LST变化具有显著的季节差异,因此通过分析不同季节的变化斜率(图6(a)—(d))和相关系数(图6(e)—(h)),更加精细地揭示LST的时空变化格局。由图6可知,四季LST的变暖趋势为冬季最强,春季其次,且夏秋季总体变化不大。冬季时高达69.4%的地区呈升温趋势(k>0 ℃·a^-1),这种显著的变暖趋势来自西北地区和青藏高原地区。这可能与该地区西伯利亚高压系统的减弱有关,西伯利亚高压系统是起源于蒙古—西伯利亚寒冷地区的强大高压系统,其带动大量冷且干燥的气团向南移动,以促进亚洲东部地区的快速降温。 Li等^[23]计算冬季温度与西伯利亚高压指数之间的相关性为0.715,表明中国西北地区的冬季温度与西伯利亚高压指数显著相关。此外,在天山山脉附近及塔里木盆地的西部边缘地区都形成了许多变暖区,具有沿绿洲分布的特点,这可能与城市化进程中人类活动产生的更多温室气体的促进作用有关。具体来说,春季的LST变化呈现北方升温南方微弱降温的空间格局,南北差异显著,变化斜率的均值为0.02 ℃·a^-1。在北方地区的黄淮海平原、大兴安岭的北部地区以及西北地区的内蒙古高原和新疆的北部地区,均出现了较集中的快速升温区域,特别是在内蒙古高原和大兴安岭的西北部的一些地区升温特别显著(k>0.18 ℃·a^-1,R>0.6)。虽然春季中国北方和西北地区的变暖有利于生长季的开始,但在黄淮海平原作为中国重要的作物种植区,温度升高使正值返青期和拔节期的冬小麦遭受干旱的频率增加且干旱程度加剧。而夏秋季节LST变化格局相似,总体变暖趋势均非常微弱,几乎没有变化。夏秋两季最显著的变暖均来自青藏高原(变化显著的地区均达到0.08以上,而且表现出强的正相关)。秋季东北地区呈现最显著的降温趋势,特别是松嫩平原附近以及三江平原的西部地区(k<-0.1 ℃·a^-1)。

1)近年来LST总体呈微弱上升趋势且空间分布不均,总体具有“北方大于南方,西部大于东部”的分布特点。其中变暖的区域占总体的63.7%,变冷的区域36.3%。最强的变暖趋势出现在西北地区,其次是青藏高原和黄淮海平原。显著的降温趋势也有体现,主要集中在东北地区。

2)白天LST变化显著高于夜间,这可能与太阳辐射的变化以及人类活动释放的大量温室气体密切相关。同时,LST的趋势在不同季节显著变化。冬季变暖趋势最为明显,变暖面积所占比例最大,特别是在干旱、半干旱的西北地区及青藏高原地区。

3)通常,温度变化影响的原因可分为以下2个方面: 一方面,气候变化受自然因素的影响,如厄尔尼诺/南方涛动等海洋变化、火山活动、太阳辐射强迫等; 另一方面人为因素的作用,包括温室气体的强迫作用、气溶胶的增加和地表覆盖的变化。人类活动释放的二氧化碳(例如化石能源消耗)对地表变暖的影响是不容忽视的。温室气体是影响温度变化的重要因素,它主要通过吸收地面反射的红外长波辐射,对变暖产生积极的反馈作用。此外,温室气体还与季风产生耦合作用导致季风区能源收支的变化,进而影响季风环流的强度。同时,温度短期变化还会受到大气污染物排放的气溶胶(例如硫酸盐气溶胶)增加的影响。气溶胶污染将加热空气流动并降低到达表面的太阳辐射产生冷却效果,从而影响局部地区甚至全球气候变化。特别是在我国工业和交通运输相对发达的东部地区,例如华北地区这种反应将更加明显。此外,在当前城市扩张过程中,建设用地、耕地和林地等土地利用的变化也会对温度变化产生持续的影响。

热红外数据反演的LST数据受云的影响很大,在云覆盖区域会造成LST数据信息的缺失。与以往利用时空信息重建晴空下的LST不同,本研究将地面站点数据应用相对多倍权重与非空像元一起用于重建云覆盖下的LST数据的过程,更能代表被云污染区域的LST情况,在提高有云情况下区域LST产品精度时更具优势。然而重建结果的精度会受到地面站点数量的影响,同时也会受到地形和地表类型的影响。在今后研究中, 可以考虑这2方面信息,并进一步增加地面观测站点的数量以进一步提高重建后LST数据的精度。