基于同类地物地表温度日变化相关性的MODIS LST重建算法

崔建勇 张曼玉 宋冬梅 罗 升 单新建 王 斌

1)中国石油大学(华东)，海洋与空间信息学院，青岛 266580 2)中国石油大学(华东)研究生院，青岛 266580 3)海洋矿物资源实验室青岛海洋科学技术国家实验室，青岛 266071 4)四川省水利水电勘测设计研究院有限公司，成都 610072 5)中国地震局地质研究所，北京 100029

0 引言

地表温度(Land Surface Temperature，LST)是研究地表能量平衡和陆面过程的重要参数(Zhaoetal.，2017)，是区域以及全球尺度土壤水分、 蒸散发、 城市气候变化、 植被、 生态监测的关键因子(Chenetal.，2017)。准确获得不同空间尺度的地表温度对地震热辐射异常的研究亦十分重要(陈顺云等，2014)，有利于地震热异常信息的精确提取(刘文宝等，2020)。不仅如此，异常是相对于正常而言的，没有正常也就无所谓异常(陈顺云等，2004)，而得到地表温度年变基准场是热异常提取的基础与前提(陈顺云等，2009)，因此高精度的地表温度重建工作对于地震热异常的精确提取是十分必要的。利用MODIS数据反演地表温度具有速度快、 覆盖面积广、 观测周期长等优点，近年来已逐渐成为获取LST的重要手段(Zhaoetal.，2020)。然而无论何时，约65%的全球表面被云层覆盖，直接导致热红外遥感影像中存在大面积分布不均的缺值(Maoetal.，2019)，去除云及其他因素的干扰并实现对云覆盖像素的地表温度进行高精度补值的方法仍亟待完善(周义等，2014)。在全球和区域尺度进行的气候变化以及地表能量平衡模型等研究中均需要长时间序列、 空间连续的地表温度数据，有效的补值工作将为基于温度的地表过程研究提供坚实的数据基础。此外，大面积的数值缺失会对地表热辐射异常的识别和分析造成严重影响。因此，云下像元地表温度重建是利用热红外遥感识别发震信号不可缺少的重要环节。

迄今为止，关于地表温度的补值方法主要有3类： 回归法、 时空反距离加权法和时空克里金插值法。对于最常见的时空自回归模型，其主要思想是将一些加权回归方法集成到局部时空预测中，如在局部加权散点图平滑方法(William，2012)的基础上发展起来的时空异质性地理局部回归方法(樊子德等，2016)。目前，时空反距离加权方法也在不断创新，利用轻量级集成模型的补值精度有所提高，且该模型在时空突变中表现出较好的鲁棒性(Chengetal.，2020)。对于克里金插值法，起初将时间特性引入到普通克里金插值法中后，补值精度明显提高(Yangetal.，2015)，随后高精度的基于时空加权回归克里金法的地理统计模型(Duetal.，2018)、 时空约束条件下的非均匀时空克里金插值法(Yongetal.，2018)和时空协同克里金法(Huetal.，2019)也相继被提出，补值精度不断提高。

随着上述方法的不断优化，LST的可用性和补值精度均得到了很大提高，但现有方法必须依靠大量辅助数据或复杂模型，且模型的适用性较差。因此，本研究旨在开发一种基于少量辅助数据和非线性回归模型的MODISLST高效率重建方法，其理论依据为一天中同类地物不同时刻LST的变化具有规律性。该方法的优势在于使用简捷的回归模型即可为空值像元生成更精确的LST，且补值过程不依赖于大量辅助数据。

1 研究区与数据

1.1 研究区概况

和田地区隶属于新疆和田市，位于新疆维吾尔自治区南部，其中心位置为(37.12°N，79.92°E)，地势北低南高，南越昆仑山口直达藏北高原，北临塔里木盆地，且从西向东缓倾。和田地区自山麓向北，戈壁纵横，位于河流冲积扇平原上的绿洲分布较为广泛，扇缘与塔克拉玛干沙漠相连，一直延伸到塔里木盆地的中心。本次实验之所以将新疆和田地区作为研究区，是由于和田地区的地物覆盖种类相对简单，大部分为沙地，有少量植被覆盖，其地形构成主要为山脉和戈壁，四周有多条山脉，中间主要为戈壁滩以及少量绿洲。同时从气候角度看，和田属于明显的大陆性气候，降水量少，晴朗的天数较多，云量相对较低，且其白天温度高，水分蒸发速度快，卫星传感器在对该区域进行探测时受水汽及云层影响较小，获得的热红外辐射数据质量较高。

综上所述，新疆和田地区的地表温度数据受云等其他干扰因素相对于中国其他地区更少，数据质量相对较好，是进行地表温度补值实验的理想研究区。

1.2 遥感数据

本文所选取的数据产品来自MODIS(中分辨率成像光谱仪)传感器，搭载于Terra和Aqua卫星，其优点为光谱范围广，具有36个中等分辨率水平的光谱波段，每1～2d对地球表面进行1次观测，并提供陆海温度、 陆地表面覆盖、 云、 水汽、 火情等多种应用产品数据的免费服务(王慧等，2018)。

根据本实验需求，选取8d合成、 空间分辨率为1km的MODISLST数据产品： MYD11A2和MOD11A2，下载2003—2015年共计13a的数据(1)https： ∥search.earthdata.nasa.gov/。。有研究表明(Kangetal.，2018)，凌晨数据缺值少，数据重建保真性好。为了获得时空连续的地表温度数据，本文将MYD11A2的凌晨数据作为重建对象，将MOD11A2上午、 晚上和MYD11A2下午的地表温度作为补值数据，相关参数见表1。此外，本研究利用MODIS的MOD12Q1产品确定每个像元的地表覆盖类型，针对各类地物分别建立回归模型，将3种补值数据拟合回归至凌晨时的地表温度，为后续地表温度补值做准备。MOD12Q1包含5种不同的土地覆盖分类系统： 国际地圈生物圈计划(IGBP)全球植被分类方案、 马里兰大学(UMD)植被分类方案、 MODIS提取叶面积指数/光合有效辐射分量(LAI/fPAR)方案、 MODIS提取净第一生产力(NPP)方案和植被功能型(PFT)分类方案(王新闯等，2013)。

表1 MODIS地表温度数据的相关参数

1.3 LST质量评判和缺值率统计

MODIS数据集中包含LST的数据质量文件(QC)，该文件描述了每个像元的温度和发射率误差，其中LSTerror≤1K的数据点对应的QC均<63，因此将QC≥63的情况视为缺值，据此对2003年1月1日研究区4种数据源(不同时间段)的缺值情况进行了统计。图2 所示的灰度值表示地表温度，单位为K，其中白色部分表示缺值区域。

图1 新疆和田研究区示意图

图2 新疆和田地区MODIS地表温度不同时段的缺值情况

从研究区不同时间段的MODIS地表温度缺值情况对比结果来看，上午和下午的地表温度数据缺值情况比晚上和凌晨更严重。其中，MYD11A2凌晨1：30的地表温度数据缺值最少，MYD11A2下午1：30的数据缺值情况最差。为了分析不同季节LST的缺值率，对新疆和田地区13a共598期不同数据产品、 不同季节LST的缺值率进行了统计分析，结果如图3 所示。从表中可以看出，秋季的数据质量相对较好，夏季和冬季的数据质量较差。从数据类型角度来看，Aqua传感器MYD11A2凌晨的数据质量相对较好，因此将其作为LST重建对象。

图3 新疆和田地区(2003—2015年)四季地表温度数据的质量统计

2 研究方法

本文地表温度重建的基本思想是： MODISLST数据可以获取同类地物地表温度在一天中随时间变化的趋势，且同类地物地表温度随时间变化具有相似的规律，据此本文提出了基于同类地物不同时刻地表温度相关性的MODISLST重建算法。由于Aqua传感器的MYD11A2凌晨数据缺值少、 数据重建保真性好，因此本研究将其作为MODISLST重建对象。即当MYD11A2凌晨数据某像元缺值时，通过对同类像元不同时刻地表温度的相关性进行分析，分别建立凌晨时刻与MYD11A2下午、 MOD11A2上午和晚上同类地物地表温度的回归模型，按照回归模型的相关性高低选取最优值，分步实现凌晨地表温度的重建。

该研究方法主要包括以下几部分工作： 首先，对所用数据集进行拼接、 重投影、 裁剪等预处理； 然后，以LAI/fPAR作为地表覆盖类型分类系统的依据，选取凌晨数据与其他3个时刻的同类地物像元，基于2阶多项式拟合分别建立凌晨与其他3个数据源各类地物地表温度的回归模型，并将三者拟合回归至凌晨时的地表温度，通过对所建模型进行分析发现晚上与凌晨的地表温度相关性最高，上午次之，下午最弱； 最后对凌晨地表温度进行两步重建： 先用晚上的拟合数据进行补值，然后取上午和下午拟合结果的最大值进行二次补值，完成对空值区域的地表温度重建。

图4 地表温度重建技术的流程图

2.1 数据预处理

本文所用的MODIS数据产品为HDF格式，借助NASA提供的MRT(MODIS Projection Tool)软件对原始影像进行拼接、 裁剪与重投影处理。首先对2003—2015年每期(8d)的3个Tile文件进行拼接，然后挑选出LST＿Day＿1km和LST＿Night＿1km，以新疆和田地区为中心，对10°×10°范围内的区域进行裁剪，将其重投影为Geographic Lat/Lon，最终得到4种地表温度数据： MOD＿DAY(上午10：30)、 MOD＿NIGHT(晚上10：30)、 MYD＿DAY(下午13：30)和MYD＿NIGHT(凌晨1：30)。

2.2 地表覆盖类型的选取

MOD12Q1地表覆盖类型数据产品提供5种分类系统，研究发现使用不同的分类系统对补值精度存在影响(图5)，因此在补值前需要选择一种适合的地表覆盖类型分类系统。随机选取1期MYD11A2凌晨地表温度数据(2015-43为2015年第43期)，采用本研究所提出的方法依据不同的分类系统分别对其空值区域进行补值，并按照不同的分类系统及精度区间统计补值结果的精度(图5)。从图中可以看出，以LAI/fPAR分类系统为依据得到的补值结果的精度均高于其他4种分类系统，其中误差≤1.5K的补值精度占比达80%以上。由此可知，使用LAI/fPAR分类系统得到的补值精度整体相对较高，该分类系统适用于后续的补值工作。

图5 5种地表覆盖类型补值精度的对比

2.3 LST回归模型的建立

地表温度重建就是根据大量已知数据求解温度随时间变化的规律和特征函数的过程，插值和拟合是实现这一过程的2种主要手段。拟合是基于大量有效值，找到一个连续曲面以最大限度地逼近这些点，而插值则是在建立连续曲面的过程中需顾及的待定点，多项式拟合是建立连续曲面用于回归的常用方法，对于不同类型的数据，要选择适合的拟合参数，否则会导致模型不稳定、 计算效率偏低。本研究建立了不同时刻、 空间相近的同类地物温度变化的时空回归模型，并分析讨论了不同拟合参数对模型精度的影响。

基于前文选取的2015年第43期凌晨地表温度数据分别对地表覆盖类型为水体、 稀树草原和裸地的MOD＿DAY和MOD＿NIGHT与MYD＿NIGHT数据进行回归拟合，地表覆盖类型为LAI/fPAR分类系统，结果如图6 所示。从图中可以看出，1、 2、 3阶多项式的拟合效果整体较好，当超过3阶时就会出现过拟合现象。表2 是地表覆盖类型分别为水体、 稀树草原和裸地时MOD＿DAY、 MOD＿NIGHT与MYD＿NIGHT进行回归拟合的效果评估，其中R2为相关系数，取值范围为[0，1]，该值越接近1，则表示模型对数据的拟合程度越好(裘伟等，2014)。随着拟合阶数的不断增大，R2也小幅增大，但复杂模型会存在精度稳定性差、 计算效率偏低等问题，故本研究选取2阶多项式进行拟合。从以上实验结果中还可以看出，3种地物的MOD＿NIGHT地表温度数据与MYD＿NIGHT地表温度数据的拟合效果相比MOD＿DAY更好，其2阶多项式拟合相关系数分别为0.940、 0.982、 0.965，而白天的相关系数为0.423、 0.890和0.590，说明晚上数据与凌晨数据的相关性最高。

图6 不同地物的MOD＿DAY和MOD＿NIGHT与MYD＿NIGHT的温度拟合图

表2 温度拟合的效果评估

2.4 两步重建方法

本研究发现，在实验区内仅使用晚上的LST对凌晨数据进行补值后仍然存在空值像元，因此需借助上午和下午的数据进行2次补值，以提高LST的补值率。这种利用夜晚数据进行第1次补值，利用上午、 下午数据进行第2次补值的方法称为地表温度的两步重建法。由于3种数据源的组合方式不同，补值精度存在较大差异，为此我们对比分析了5种补值方案，最终确定两步重建的具体方案为： 首先，使用晚上的数据对凌晨LST进行补值，当存在缺值时再用上午和下午数据的最大值对其进行二次补值。地表温度重建算法的具体流程如图7 所示。

图7 MODIS地表温度两步重建的算法流程图

通过2.3节拟合分析可知，上午(10：30)、 下午(13：30)、 晚上(22：30)3种数据与凌晨(1：30)数据的相关性之间存在较大差异。随机选取9期数据，建立以下5种补值方案。方案1： 单独使用晚上数据补值； 方案2： 取上午、 下午和晚上数据的最大值补值； 方案3： 取上午和下午数据的最大值补值； 方案4： 单独使用上午数据补值； 方案5： 单独使用下午数据补值。分别从补值精度和补值率2个方面评价5种方案的补值效果，其结果如表3 所示。

表3 5种补值方案的效果对比

由表3 可知，单独使用晚上数据(方案1)对凌晨地表温度进行补值具有最高的精度，平均误差为1.06K。由于白天数据受云层影响而缺值多，故当加入上午和下午数据联合补值时(方案2)精度有所下降。但晚上数据也会存在缺值的情况，此时若单独使用晚上数据进行补值，则无法对缺值区域的凌晨LST开展补值工作。而当三者联合补值时补值率可达100%，因此本文考虑对凌晨地表温度进行两步重建。根据回归模型相关性的高低和本节补值精度的对比分析确定两步重建的具体方案为： 首先使用与凌晨数据相关性最高的晚上数据对空值区域的LST进行补值，对于剩余的小部分未补值的空值像元，取同一像元的上午、 下午与凌晨数据拟合结果中的最大值进行二次补值。此方法在确保补值精度的前提下提高了补值率，从而保证了重建数据的时空连续性。

3 研究结果

为了检验本文所提方法的有效性，选择2003—2015年间的9期影像，随机选取窗口大小分别为131×154(图8a)和144×171(图8b)的区域进行零值化处理，然后按照本文方法对该区域的地表温度进行补值，将补值结果与该窗口的原有地表温度数据进行对比分析，地表温度重建前后的精度对比如表4 所示。从表中可以看出，该方法在2个不同区域的重建精度均较高，其中区域1地表温度的重建精度为0.92K，区域2地表温度的重建精度为1.05K。

图8 零值化区域示意图(白色区域)

表4 重建前后的精度对比

与此同时，为了验证本文所提方法对其他3个时间段地表温度数据重建的普适性，按照本文所述方法流程，分别对区域1和区域2的上午、 下午、 晚上数据也进行了补值，结果如表4 所示。从表中可以看出，基于该方法可以实现一天中4个时间段的温度重建，且凌晨和晚上的地表温度重建精度比上午和下午高。为了更加直观地展示基于本文方法的重建效果，以2015年43期为例，对2个子区域的地表温度重建前后的效果进行对比分析。从图9 中可以看出，重建后的地表温度数据与原始数据基本一致，且能够较好地保留原始影像的纹理信息。因此，在以后的地震热异常提取实验中，可选择缺值较少、 重建精度较高的凌晨和晚上数据作为异常提取对象。

图9 不同区域、 不同时间地表温度重建前后的效果对比(以2015年第43期为例)

对比图10a 和图10b 可见，基于本文方法的地表温度重建结果与原图像差异不大，纹理信息一致性较好，证实了本方法的有效性。随后，又分别对补值前后共9期数据的补值精度和补值贡献率进行了统计分析，如表5 所示。在补值精度方面，补值后的最小误差为0.57K，9期的平均误差为0.92K，且均<1.2K； 在补值贡献方面，晚上数据的补值贡献占比平均为0.996，上午和下午的数据补值贡献极小。综上所述，本文方法可基于少量辅助数据获得较高精度的补值结果，大大提升了地表温度补值效率。

图10 地表温度重建前后对比

表5 补值精度对比

为了进一步验证本方法的有效性，选取与本文所用数据源相似的MODISLST补值方法(李楠等，2018)作为对比。该方法使用了较为经典的空间域加权平均插值算法，充分利用了邻近像素属性的空间相关性和同类地物LST变化的相似性，其补值率和补值精度较高。具体方法为： 若数据A在(x，y)处存在空值，则用数据B(x，y)的值乘上A、B数据所在子区域的均值之比，结果见图10c。通过对比可以明显看出，对比方法的补值结果比原图稍微偏高(颜色越亮表示值越大)，其中2004-43期和2015-25期影像的重建效果相对较好，2009-38期LST重建后影像子区域的分界线处存在显著的阶跃式变化，且与原始影像差别明显。由此可知，此方法补值精度的稳定性受各子区域因素的影响较大。

表6 现有方法的重建精度对比

使用对比方法补值后的最高精度为0.89K，平均精度为1.13K，最大误差达1.6K，精度低于本文所提方法。综上所述，本文所提方法的补值精度较高，算法稳定性较强，避免了在对比方法中出现的“偏差”。此外，还将本文方法与近5年提出的地表温度补值方法的补值精度进行了对比分析。 由表6 可知，本方法在使用少量辅助数据的前提下达到了良好的补值效果，再次证明了本文所提方法的优越性。

4 结论

MODISLST产品因存在大量空值导致其应用受限。为了获得时空连续的地表温度产品，本文以一天中同类地物不同时刻地表温度的变化规律为理论依据，提出了一种基于同类地物地表温度日变化相关性的MODISLST重建算法。该算法主要包括以下2个步骤： 首先，建立上午、 下午和晚上分别与凌晨地表温度的回归模型，将3个时刻的地表温度拟合回归至凌晨，然后使用两步重建法实现对凌晨LST的重建。本研究的主要结论如下：

1)对不同地表覆盖类型的MOD11A2上午、 晚上和MYD11A2下午与MYD11A2凌晨地表温度建立回归模型分析可知，2阶多项式的拟合效果更优，且模型稳定性较好； 2)对于该研究区，使用夜晚数据对凌晨地表温度进行补值的精度最高，但不能完全保证重建地表温度的时空连续性，若加入上午和下午地表温度数据，则补值率可进一步得到提高； 3)基于本方法的地表温度重建数据的区域连贯性更好，不会出现阶跃式变化。通过与近5年提出的LST补值方法的对比发现，本方法所用辅助数据较少，且补值精度相对较高，补值后误差均在1.2K以下，平均误差为0.92K，最小误差为0.57K，可满足地表温度补值要求，为后续地震前兆热异常精细提取研究提供数据支撑。

虽然本文方法可以实现地表温度的高效补值，但仍存在以下不足： 由于用来补值的辅助数据自身也存在空值，以致单独使用晚上数据进行补值后的数据产品中仍存在缺值，尽管加入上午和下午数据进行第2步补值可提高补值率，但不可避免会引入误差，未来可考虑使用更优质的数据源替代白天地表温度数据。此外，该方法只对新疆和田区域进行了地表温度补值验证，未验证算法的普适性。未来，我们将在算法的区域适用性方面不断开展深入研究，在实验中考虑加入S-G(Savitzky-Golay)滤波或深度学习模型以进一步完善地表温度重建，从而提升地表温度补值方法的普适性。