云下遥感地表温度重构方法研究

周芳成, 唐世浩, 韩秀珍, 宋小宁, 曹广真

(1.国家卫星气象中心，北京 100081; 2.中国科学院大学，北京 100049)

0 引言

地表温度是地表-大气之间水热平衡中一个重要的参数，不仅可用于判断火灾和地震带、找矿和地热资源、研究城市热岛效应[1-2]； 还直接影响着蒸散发、土壤含水量、植被和土壤生化特性、大气可降水量、区域CO2含量等重要参数[3-8]。在大气、生态、水文和生物地球化学等诸多领域的研究模型中，地表温度是常见的输入变量，其精度直接影响着模型的输出精度[9]。因此，获取高精度和连续性的地表温度数据具有非常重要的意义。

虽然地面测站可以获取高精度和全天候的地表温度观测数据，但是稀疏分布的地面测站无法覆盖地表温度强烈的空间异质性，使得站点观测可能会遗漏大尺度地表温度的空间变化信息。随着对地观测卫星的发展，遥感技术被认为是唯一可以在全球尺度保证高时间分辨率和空间全覆盖的地表温度观测手段[10]。基于遥感技术的地表温度反演成为大尺度地表温度获取的有效手段。当前，基于热红外遥感的地表温度反演理论与方法已经较为成熟，可以获得较高精度的晴空地表温度遥感产品[11-12]，但是受制于热红外遥感无法穿透云层，有云地区地表温度遥感数据存在缺失，极大影响了地表温度数据的全天候应用需求。对云下缺失地表温度遥感数据进行重构，获得全天候地表温度遥感数据、满足其全天候应用的需求是当前研究的热点。

当前，国内外已发展出很多缺失数据重构方法，如最优插值法(optimal interpolation, OI)、经验正交函数分解法(empirical orthogonal function, EOF)、期望最大化法(expectation maximization, EM)、奇异谱方法(singular spectrum analysis, SSA)、卡尔曼滤波(Kalman filter, KF)、本征模态分解法(proper orthogonal decomposition, POD)、变分资料同化(variational data assimilation, VDA)等[13]。在国内，赵冰等[14]利用邻近非空像元采用反距离加权回归法重建缺失像元； 吴迪等[15]以FY-2F地表温度日均值产品为数据源，利用地表温度时间序列特征，开展了基于Savitzky-Golay滤波算法的地表温度长时间序列的重建研究。上述方法通常受到输入数据时间序列长度、地表覆盖类型等因素影响。本文发展的第1种重构方法是借助地表温度同化数据集发展的一种时空匹配的数据融合方法； 第2种重构方法是已经在海表数据重构方面获得了广泛应用的基于EOF法的经验正交函数插值法(data interpolating empirical orthogonal function, DINEOF)方法。DINEOF方法在陆地参数的重构方面研究较少，本文尝试将其应用于地表温度重构的研究中并评价其精度。2种方法基于不同的输入数据，可为未来地表温度的全天候获取研究提供有益的帮助。

1 数据与方法

1.1 数据

相比白天，夜间的地表温度空间覆盖数据具有更小的空间异质性，因此本文用到的原始地表温度遥感数据是来自Aqua卫星搭载的中分辨率成像光谱仪(Moderate-resolution Imaging Spectroradiometer, MODIS)的夜间晴空地表温度数据(MYD11C1)，0.05°空间分辨率，过境时间约为当地时间01∶30，下载自美国国家航空航天局网站。为了满足方法1对输入数据的要求，本文使用了中国气象局陆面数据同化系统(China Meteorological Administration Land Data Assimilation System Version 2.0, CLDAS-V2.0)的地表温度同化数据集。CLDAS的地表温度数据集是中国气象局推出的覆盖亚洲区域(N0°～65°，E60°～160°)、全天候、具有0.062 5°空间分辨率与1 h时间分辨率的等经纬度网格融合分析产品，该数据集的优点是利用了多源的地面、卫星等观测资料，在中国区域质量和时空分辨率都优于国际同类产品。为了满足方法2对输入数据的要求，下载了2017年全年的每日MYD11C1数据以构成时间维。

为了对比2种方法的适应性和精度，本文分别使用遥感数据交叉验证和地面实测数据直接验证2种方法。遥感数据验证选择了位于不同经纬度、不同土地利用类型、不同气候条件的3个验证区(图1)，分别位于塔克拉玛干沙漠地区(验证区1，裸土、少云)、山东半岛地区(验证区2，农田、多云)、南方某地(验证区3，林地和草地、多云)。重构前，首先对3个验证区产生人造云(即将原有的有效值改为空值)，重构后对比原始值与重构值的精度。地面验证数据是来自中国地面气象站逐小时观测资料的地表温度数据，用于验证有云条件下重构值的精度。精度评价指标为均方根误差(root mean square error, RMSE)和偏差(Bias)。在2017年全年数据中随机选取1月、4月、7月、10月的一天数据分别代表冬季、春季、夏季、秋季，进行时间上的验证。

图1 验证区的位置及土地利用类型(底图为MODIS的2017年土地利用类型产品MCD12C1)Fig.1 Locations and land use types of verification zones

1.2 两种云下缺失数据重构方法

1.2.1 同化数据重构法

同化数据通常融合了多源地面和卫星观测资料，结合同化、插值、反演、校正等技术，一般具有较高的精度，且不受气象条件影响，可以实现数据的全天候获得。当前常见的同化数据集有全球陆面数据同化系统、北美陆面数据同化系统、欧洲陆面数据同化系统、加拿大陆面数据同化系统、中国陆面数据同化系统等。本文以中国气象局的CLDAS地表温度数据集为例，首先将其与MODIS的地表温度产品进行时空尺度匹配，然后将时空匹配后的地表温度值赋给MODIS云下缺失像元。CLDAS地表温度数据集具有全天候、高时空分辨率、高精度的优势，利用中国区域业务的质量控制后的地表温度自动站观测资料对CLDAS地表温度数据集进行评估，全国区域平均相关系数为0.98，RMSE为1.8 K，偏差为1.4 K。但是CLDAS是1 h分辨率的全范围地表温度，与极轨卫星的过境时间地表温度具有差异性，特别是中国区域横跨多个时区，如果用单一时刻的CLDAS地表温度去重构会产生较大的误差。本文发展了一种时空数据匹配的融合方法，使CLDAS的地表温度尽量符合极轨卫星地表温度，其工作流程为：

1)根据MYD11C1携带的晴空像元卫星过境时间信息，对有云像元的过境时间进行插值。

2)将各个像元的卫星过境时间进行四舍五入取整，以符合CLDAS数据1 h时间分辨率的要求，同时也保证了重构的地表温度因为时间差异而造成的误差控制在0.5 h以内。

3)对CLDAS的地表温度数据集进行空间分辨率的降尺度操作，将原始的0.062 5°降尺度到与MYD11C1相同的0.05°，由于空间分辨率相差不大，忽略由于降尺度可能产生的误差。

4)根据有云像元的过境时间和位置信息寻找与之对应时刻和位置的CLDAS地表温度数据，并将该值赋给MYD11C1缺失像元。

5)由于CLDAS和MODIS在晴空像元处的地表温度值也有差异，因此可以假设缺失点处也存在同样的误差。因此，本文对缺失点处的初始猜测值又提出了二次订正的思路： 用非缺失点处的重构误差帮助订正缺失点处的猜测值，对猜测结果进行二次优化。计算公式为：

ER=Vrecon-Vorig，

(1)

Vrecon,after=Vrecon,before-ER，

(2)

式中：Vorig为非缺失点处的原始值；Vrecon为非缺失点处的猜测值；ER为重构误差；Vrecon,before和Vrecon,after分别为缺失点处二次订正前和二次订正后的猜测值。

1.2.2 DINEOF重构法

DINEOF法是由Beckers等[16]提出的一种基于经验正交函数分解法来重构时空场中缺失点的方法，该方法需要大量的历史数据构成时间维，通过对缺失点的历史数据分析提高重构精度。相比于经典的OI法[17]，DINEOF法具有无需先验值、自适应、效率高的优点[13,18]，目前已经在海表温度、叶绿素a浓度、海面风场数据的重构应用中发挥了重要作用[19-22]。当前国内文献对DINEOF方法在地表温度重构上的研究并不多，国际上，Zhou等[23]将其应用于在了青藏高原阿里地区的地表温度重构。DINEOF法的基本原理与工作流程可参考相关文献[16,18]。值得注意的是，DINEOF模型在数据重构的赋值过程中，在缺失点处赋以猜测值，在非缺失点处则同时存在猜测值和原始值，虽然最终只保留了原始值，但在非缺失点处其实存在一个潜在的重构误差(猜测值与原始值的差异)，相对应的，可以假设缺失点处也存在同样的重构误差。此处采用与同化数据重构法相同的二次订正对缺失点重构值进行纠正。

2 结果与分析

2.1 分布趋势

图2—5分别是2017年的冬季(1月9日)、春季(4月27日)、夏季(7月12日)和秋季(10月12日)共4天的地表温度分布图。

从图中看出，MYD11C1晴空地表温度分布图由于有云遮挡，地表温度分布呈破碎状，特别是南方部分地区，云覆盖范围大且覆盖时间长，极大影响了对地表温度全天候应用的需求； 而经过2种重构方法重构后的地表温度分布图修复了地表温度分布的破碎性，更好地表现出随季节的更替(冬-春-夏-秋)，北部的低温区域(蓝色)消失，南部的高温区域(红色)从南向北推进又回到南部的趋势。从视觉效果上看，重构图不仅更加美观，而且具有与地势起伏和时间变化相符合的丰富的地表温度纹理特征，体现了地表温度分布的空间异质性和时间差异性，其分布趋势较为合理，表明2种重构方法均有较高的可信度。

2.2 遥感数据验证

遥感数据交叉检验可以在同一个空间尺度上对重构值进行评价，原始的MODIS地表温度作为真值，产生人造云对其遮挡，通过对比原始值与重构值来评价其精度。本文对3个验证区所属的土地利用类型进行了评价(4季4天)。图6— 8分别是晴空条件下2种重构方法在验证区1(沙漠)、验证区2(农田)和验证区3(林草混合地)的重构精度的对比散点图。验证区1塔克拉玛干沙漠地区的数据点分布在3个区域，主要是因为所选验证区在10月12日当天无有效数据； 验证区2山东农田地区的数据点也分布在3个区域，分别代表了冬季(低温)、夏季(高温)和春秋(中温)，表明山东农田地区的地表温度有较为明显的四季变化，春秋季节地表温度则较为接近； 验证区3林草混合地区数据点主要集中在285 K和293 K附近的2个区域，说明南方林草混合地区全年较为温暖。从图中看出，不同地类条件下2种方法的精度均较为接近，其中农田和沙漠地区的RMSE都小于2 K，Bias都小于1 K，精度较好； 林草混合地的RMSE约为2.7～2.9 K，精度略差的原因体现在2方面： 一是南方地区有云天数较多，有效值较少，因此DINEOF法条件下可参考的时间维信息较少； 二是同化数据重构法当中CLDAS主要表现的是土壤表层温度信息，而MODIS在高植被覆盖区则主要体现植被冠层温度，两者存在一定误差。从结果看，2种方法在少云、地类单一的地区重构精度较好； 在多云、地类复杂地区重构精度略有降低。

2.3 地面实测数据验证

重构是为了解决云下地表温度数据缺失的问题，因此本文更关注2种方法在有云地区的重构精度。有云地区由于无遥感地表温度产品，只能依赖于地面站点的实测值进行验证，站点实测值是基于点尺度，而本文重构的地表温度是基于面尺度(0.05°空间分辨率)，2种尺度的地表温度存在一定的误差，因此本文在利用地面实测数据验证前，首先建立了地面实测数据和晴空遥感数据之间的尺度匹配关系：

LSTRS=aLSTground+b，

(3)

式中：LSTRS为地面站点所处像元在晴空条件的MODIS地表温度值；LSTground为地面站点测量的地表温度值；a和b为系数。将式(3)应用于有云地区，计算出像元尺度地面实测值，缩小尺度差异对结果的影响。

地面实测数据的验证首先根据中国地面气象站所处经纬度判断每一个站点是否处于有云地区，判断为“是”的站点地表温度参与验证。图9是云下地面站点数据和2种重构地表温度的对比散点图。从图9中看出，冬季整体精度略差于其他3个季节，2种方法的RMSE分别为3.21 K和3.20 K； 秋季精度最好，2种方法RMSE均在2.5 K左右； 春季2种方法精度相差较大，同化数据重构法的RMSE为3.05 K，DINEOF重构法只有2.54 K； 夏季的散点数量最多，因为夏季云量多，表明大部分站点地区都被云遮挡，2种方法在夏季的RMSE分别为3.08 K和3.02 K，相差不大。总体来看，2种方法在春季精度相差较大，达到0.5 K，DINEOF重构法较优，在其他季节2种方法精度比较接近； 2种方法在有云地区全年精度可以达到2.5～3.5 K之间。

3 结论

本文发展了2种方法来重构全天候地表温度遥感产品，并对其进行了验证和评价： 方法1是基于同化地表温度数据集的重构方法，本文发展了一套时空匹配融合方法，减小了同化数据获取时间与极轨卫星过境时间下地表温度的差异，并根据晴空条件下2种地表温度的系统误差对云下数据二次订正以消除重构误差； 方法2是将当前在海表参数重构研究中较为流行的DINEOF法应用于地表温度重构研究，并对其增加了二次订正来消除重构误差。

通过对重构数据的地表温度分布图观察，2种方法都能较好地重构云下缺失地表温度，并能体现出地形差异和时间差异导致的地表温度纹理特征，从空间分布趋势来看，2种方法的重构结果都较为可信。从对遥感数据的验证结果来看，2种方法在沙漠、农田、林草混合地3个验证区的RMSE均小于3 K，Bias小于1 K； 通过对地面实测数据的验证结果来看，在预先消除验证数据与重构数据的尺度差异之后，2种方法在有云地区全年精度可以达到2.5～3.5 K之间； 2种方法精度在夏秋冬比较接近，春季DINEOF重构法优于同化数据重构法。