基于SCLSTM模型的MODIS地表温度产品重建方法

宋冬梅 张曼玉 单新建 王 斌

1)中国石油大学(华东)海洋与空间信息学院,青岛 266580 2)海洋矿物资源实验室,青岛海洋科学技术国家实验室,青岛 266071 3)中国地震局地质研究所,北京 100029

0 引言

地表温度是区域及全球尺度上地表能量平衡和水循环物理过程的关键参数(Zhao Betal.,2020),可用于干旱监测(Xieetal.,2021)、城市热岛效应(闫章美等,2021)与地震前热异常信息提取(Songetal.,2018; Chenetal.,2020)等。准确获得时空完整的LST数据有助于精确提取地震前兆热异常信息(陈顺云等,2014),不仅如此,异常是相对于正常而言的(陈顺云等,2004),而地表温度背景场是热异常信息提取的基础与前提(陈顺云等,2009),因此重建LST对地震热异常的相关研究十分必要。MODIS传感器具有观测范围广、时间分辨率高的优点,是获取LST数据的重要来源。然而,在任何时刻65%的全球表面都会被云层覆盖(Maoetal.,2019),这将直接导致热红外遥感影像中存在大面积缺值,且各地的缺值区范围不等。其中,中国大陆地区MODIS LST数据的年均缺值率最高可达20%(臧琳,2017),特别是中国西南地区缺值最为严重,导致该区域MODIS地表温度数据的可用性降低。因此,MODIS地表温度数据的重建工作具有十分重要的意义。

迄今为止,关于地表温度的重建方法主要包括以下4种:1)基于空间域信息的地表温度重建方法。该类方法主要利用克里金插值(Neteler,2010)、样条函数(Kilibardaetal.,2014)、回归树分析(Fanetal.,2014)、反距离加权(张军等,2011)等算法,基于缺失像元与其相邻晴空像元之间的空间相关性实现地表温度的重建。此类方法易于实现,无需其他辅助信息。然而,由于该类方法受到地表覆盖变化和地表高程变化的影响较大(Wuetal.,2021),导致其LST重建精度不够理想。2)基于时间域信息的地表温度重建方法。该类方法主要利用日温度循环建模(Zhangetal.,2021)、谐波分析(Nguyenetal.,2015)、多时相鲁棒回归(Zengetal.,2018)、奇异谱分析(Malamirietal.,2018)、物理建模(Fuetal.,2019)等算法,基于缺失像元与其时间域上相邻的晴空像元之间的相关性实现地表温度的重建。尽管该类方法在通常情况下重建效果较好,但当相邻期的缺值数据较多时,其重建效果将变差。3)综合时空信息的地表温度重建方法。该类方法主要基于RSDAST(Sunetal.,2017)、多元回归(Zhao Wetal.,2020)、地表能量平衡(Tanetal.,2021)等算法,一般先在时间域上进行补值,然后再在空间域上进行补值,从而实现地表温度数据的重建。该类方法虽然在重建过程中综合考虑了时间和空间信息,但并没有实现真正意义上的时空联合建模,因此制约了其重建精度的进一步提升。此外,由于该类方法既顾及了时间信息又顾及了空间信息,使得重建过程较为复杂。4)基于多传感器数据的LST重建。该方法的主要思想是将PMW(Passive Microwave,被动微波)LST数据与TIR(Thermal Infrared,热红外)LST数据相结合,以重建LST数据(Kouetal.,2016; Prigentetal.,2016; Shwethaetal.,2016)。PMW(苏扬等,2022)可穿透云层获得全天候LST数据(Zhang Qetal.,2020),为TIR LST数据的重建提供有用的辅助信息。因此,通常使用PMW LST数据估算云下LST数据。其中,结合TIR LST(如MODIS)和PMW LST的LST重建方法不断被提出,此类方法在生成高精度、时间连续的LST数据方面具有很大的潜力。然而,由于PMW特定的扫描方式导致相邻2个轨道之间存在扫描间隙(Zhang Xetal.,2020),在中低纬度地区缺失宽度可达数百千米,使PMW获取的LST数据空间不连续。此外,由于PMW LST和TIR LST数据的物理意义不一致,基于PMW观测数据反演的LST不能代表实测的地表温度,且2种数据的时空分辨率也不一致(Zhouetal.,2017),通常需要在转换后使用(Moncetetal.,2011)。因此,虽然PMW可以有效地辅助云下TIR LST数据的重建,但其数据本身的局限性阻碍了该方法的广泛应用。

近年来,随着深度学习方法的应用不断深入,越来越多的学者尝试在LST数据重建工作中使用该类方法。其中,LSTM模型由于能够充分学习时序数据中的长期依赖关系,逐渐被应用于地表温度重建工作中(Cuietal.,2022)。但该模型对于时序数据中极值点信息的表达能力较差,导致部分缺值点的重建效果不理想。为此,有学者提出了基于BiLSTM(Bi-directional Long Short-Term Memory,双向长短时记忆)模型的LST重建方法(Arslanetal.,2019),试图解决由于存在极值点导致重建效果欠佳的问题(Cuietal.,2022)。但由于该模型对时间序列进行双向学习,导致模型的训练时间长、效率低。时间序列预测模型中除了最常用的LSTM和BiLSTM外,一维卷积神经网络由于能够很好地挖掘出数据本身的潜在特征,因此被逐渐应用于地表温度数据的重建工作中。例如,Wu等(2019)提出了一种基于多尺度特征卷积神经网络的地表温度重建方法,该方法根据在不同时间获取的同一区域的地表温度之间的非线性关系实现数据重建,虽然原理较为简单,但它要求辅助的LST图像不能有缺值,其苛刻的数据要求限制了算法的适用性。

LST数据重建本质上是一种数据预测工作(Cuietal.,2022)。而随着数据预测研究的不断深入,学者们发现将数据分解模型和预测模型相结合的混合模型比单一模型具有更好的预测效果(Yangetal.,2021)。其中,SSA相比于经验模态分解等其他数据分解方法(Cuietal.,2022),能够更加充分地挖掘数据潜在的周期及趋势特征,所获得的数据特征更加丰富。时间序列预测模型中的LSTM具有学习数据长期依赖关系的能力,而CNN能够充分提取一维时间序列数据的局部特征,因此,将CNN与LSTM结合则可以更好地发挥模型的预测优势。为此,本文提出一种基于混合模型的地表温度两步重建方法——SCLSTM(SSA-CLSTM)。该方法包括2个主要步骤:首先,使用SSA模型提取地表温度时间序列的趋势特征并填补缺值像元,实现地表温度的初步重建。然后,再利用SCLSTM(即1DCNN-3层堆叠LSTM)模型学习数据的局部时序特征和长期依赖关系。其中,1DCNN负责局部时序特征的挖掘与学习,而3层堆叠LSTM负责数据的长期规律的学习,并通过LSTM模型的3层串联对缺失像元的地表温度进行迭代预测,从而实现缺值像元LST数据的精细重建。

1 研究区与数据

1.1 研究区

本文的研究区主要包括干旱少雨、植被覆盖稀少的新疆和田地区和湿润多雨、植被生长茂盛的四川汶川地区,研究区的概况如图1所示。和田地区位于中国西北,地处新疆南部,其中心位置为(37.1°N,79.9°E),地势北低南高。该地区属于典型的大陆性气候,全年热量充足、降水稀少,晴空天数较多。汶川地区位于四川省中部,其中心位置为(31.0°N,103.4°E),地势西北高东南低。该区域在夏季受到东南季风和西南季风的共同影响,导致气温高、雨水多,夏季降水占全年降水量的55%以上(倪楠,2020)。因此,较多的云雨天气严重影响了该地区MODIS地表温度产品的数据质量。据统计(臧琳,2017),中国大陆的地表温度年均缺值率最高可达20%。其中,中国西南区域的缺值情况最为严重,空间分布特征为东南高、西北低。与新疆和田地区相比,四川汶川地区的地表温度缺值情况更加严重。本文选择云量不同的区域进行重建实验,以便更好地评估新方法的可靠性及区域适用性。

图1 研究区概况

1.2 地表温度数据

本文采用新疆和田地区2003—2008年以及汶川地区2015—2020年空间分辨率为1km、8d合成的MODIS MYD11A2 LST凌晨(1:30am)数据产品为实验数据。所用MODIS数据产品为HDF格式,需借助NASA提供的MRT软件对原始影像进行拼接、裁剪与重投影(Geographic Lat/Lon)处理。其中,2个研究区的范围分别是:新疆和田(30°～40°N,75°～85°E); 四川汶川(26°～36°N,98.5°～108.5°E)。

2 研究方法

利用单一预测模型进行LST重建时可能会引入与地表温度毫无关系的噪声,导致重建效果较差。为此,本文提出了一种基于混合模型的地表温度两步重建方法——SCLSTM(SSA-CLSTM),该方法的具体流程如图2所示。

图2 重建方法的流程

2.1 基于SSA方法的LST初步重建

为了避免时间序列中连续缺失数据对地表温度精细重建的影响,在精细重建前需要对原始数据中的缺值进行初步重建。SSA是近年来兴起的一种用于提取非线性时间序列数据特征的强大方法(Zengetal.,2018)。该方法通常首先根据所观测到的时间序列构造出轨迹矩阵,并对轨迹矩阵进行分解,提取出原时间序列中不同成分的信号,如长期趋势信号、周期信号和噪声信号等,然后依据分解后所得到的不同成分的数据系列对原始数据贡献率的大小进行数据重构(陈浩文,2020)。需要说明的是,取过大或过小的贡献率都不利于数据的有效重建,这是由于过大的贡献率会使得重构后的数据仍然包含很多噪声,从而无法很好地得到时间序列数据的趋势特征,而过小的贡献率虽然能够很好地滤除噪声,但却使得时间序列数据中的趋势特征过于平滑。本研究所用的地表温度数据产品在一年中有46个值,因此本文将SSA的窗口长度设置为46。使用SSA对每个像元的时间序列进行分解后可以得到46个长度为46的子序列数据。

本文拟采用上文提到的SSA方法提取出地表温度时间序列中的趋势值用以填补缺值像元,实现数据的初步重建。基于SSA方法的LST初步重建过程如图3所示,其具体步骤为:首先,利用SSA方法对缺值像元的地表温度时间序列进行分解,得到包含不同特征的子序列数据。然后,计算分解后各子序列的贡献率,选择贡献率之和 >85% 的前r个子序列,其中,子序列贡献率的计算过程如式(1)所示。最后,将所选出的各子序列对应时间段的数值进行相加,并将所得的计算结果填充至缺值像元,完成地表温度时间序列的初步重建。

(1)

图3 基于SSA方法的LST初步重建过程

式中,Cj表示前j项子序列的贡献率。

2.2 基于1DCNN-3层堆叠LSTM方法的LST精细重建

CNN是一种深度神经网络,它可以从复杂的数据中挖掘其潜在特征。近年来,CNN已越来越多地被应用于目标识别和图像分类等领域。其模型一般由输入层、卷积层、池化层、全连接层和输出层5个部分构成(Wuetal.,2020)。其中,卷积层的作用是通过对输入数据进行特征提取从而得到更具表现力的特征序列,池化层的作用是通过对输入特征进行降采样处理以强化部分显著特征(陈赐阳等,2021)。不同维度的CNN具有类似的结构和处理步骤,主要差异在于输入数据的维数、卷积核的大小及卷积核在数据之间的滑动方式。由于本文所提出的地表温度数据重建方法是基于像素的一维时序数据进行的,因此采用一维CNN(1DCNN)模型。其中,一维卷积层的结构及卷积核的滑动方式如图4所示。

图4 一维卷积层结构及卷积核滑动方式示意图

为了克服RNN网络的梯度消失问题,Hochreiter等(1997)提出了LSTM模型,是一种时间循环神经网络。LSTM模型通过引入遗忘门、输入门和输出门,并通过“门(gate)”的控制决定网络中哪些信息被保留、哪些信息被遗忘(Volkanetal.,2020)。而网络中的“细胞状态”通过有选择性地接收遗忘门和输入门的信息得以更新。更新后的“细胞状态”和输出门的信息共同决定最终的输出结果,并基于迭代运算得到训练好的优化网络。本研究所用LSTM模型的预测原理如图5所示,其原理为利用前M个时刻的地表温度值预测第L时刻(缺值时刻)的地表温度,并不断迭代,直至完成时间序列上所有缺值像元的重建。

图5 LSTM模型预测示意图

基于1DCNN-3层堆叠LSTM模型的地表温度精细重建过程如图2所示。本文所提出的精细重建模型主要包括1DCNN和3层堆叠LSTM 2个部分。其中,1DCNN由5个卷积层组成,用于提取每个时间步中时间序列的主要特征,以减少输入到LSTM网络中的冗余信息。3层堆叠LSTM是由3个LSTM串联组成,主要用于下一时间步的LST数值预测。LSTM层数是时间序列数据预测中的关键参数,本研究通过逐渐增加LSTM的层数,并观察地表温度重建误差的变化,从而选取出合适的LSTM层数,结果如图6所示。从图中可以看出,当LSTM的层数为3时,模型的误差最小,重建后数据与原始数据的相关性最高,因此本研究将LSTM的层数取为3。将3层LSTM进行堆叠,既能增加LSTM循环层及网络容量,也可进一步提升模型的预测能力。本研究中CNN和LSTM模型的详细参数信息如表1所示。其中,SCLSTM模型的输入内容为带有数据缺失的一维时间序列数据,输出内容为重建后完整的时间序列数据。其中,训练与测试数据的比为0.85,迭代次数为250,批数量为125。

表1 CNN与LSTM模型的参数信息

图6 SCLSTM模型中不同LSTM层数的数据重建精度

基于1DCNN-3层堆叠LSTM模型的地表温度精细重建过程可概括为以下2个主要步骤:

(1)基于1DCNN的LST局部时序特征提取

首先,将基于SSA方法获得的LST时间序列的初步粗重建结果输入到1DCNN模型中。然后,通过1DCNN的卷积运算提取出LST时间序列的初步特征。之后,将该初步特征作为第2个卷积层的输入内容,再次执行卷积运算,由此得到一个更高级的特征,以此类推,直至完成所有卷积运算,其操作的具体过程如图7所示。进行卷积计算后,池化层把最后一个卷积层计算得到的特征序列作为输入内容,并使用池化窗口在该序列上进行滑动,每滑动一次取窗口的平均值进行池化,以输出更具表现力的特征序列。

图7 卷积核权重与LST数据局部序列段之间的卷积运算

(2)基于3层堆叠LSTM模型的LST精细化重建

为使LSTM模型具有更强的表达能力,本研究采用3层堆叠LSTM模型的策略以实现LST的精细重建,过程如图8所示。首先,将1DCNN模型提取出的时序特征作为第1层LSTM模型的输入内容,并将其输出结果作为第2层LSTM的输入内容,以此类推,完成3层LSTM的预测,从而实现对LST时间序列中第1个缺值时刻数据的精细化重建。为了完成其他所有缺值时刻数据的精细化重建,迭代实施1DCNN-3层堆叠LSTM模型,过程如图9所示。综上所述,由于该模型同时使用了CNN和LSTM,保证了地表温度时间序列中的突变信息和趋势特征均能够被准确提取,这使得地表温度数据的精确重建成为可能。

图8 3层堆叠LSTM精细化重建过程

图9 LST时序数据的迭代精细重建过程

2.3 对比方法介绍

本文所采用的2种对比方法分别为基于SSA-LSTM和基于SSA-BiLSTM(Cuietal.,2022)的地表温度重建方法。2种方法的共同之处均是首先利用SSA算法对缺值像元的地表温度的时间序列进行分解,得到包含不同特征的子序列数据,然后将不同子序列数据重构的LST值作为预测模型的输入内容。2种方法的不同之处在于所使用的预测模型不同,前者使用LSTM模型,后者使用BiLSTM模型。

为了评价模型的重建效果,本研究使用“去除-重建-对比”的策略分别对所提方法和2种对比方法的重建效果进行了定性和定量分析。其中,“去除-重建-对比”的具体步骤为:首先,在完整的LST图像中去除部分已有数据并用零值填充; 然后,分别使用3种方法对缺值像元的LST进行数据重建; 最后,将重建后的LST值与原始数据进行对比分析,并采用RMSE(Root Mean Square Error,均方根误差)、R2(相关系数)和AD(平均偏差)进行精度评定,计算公式如式(2)—(4)所示。此外,为了进一步地比较3种方法的重建效果,本文还依据重建数据的误差大小,按照<1K、1～1.5K、1.5～2K和>2K 4个等级进行详细分析。

(2)

(3)

(4)

3 研究结果

3.1 不同重建方法的精度对比

为了验证SCLSTM的重建效果,比较3种重建方法的重建精度,在新疆和田地区使用不同方法进行LST重建对比实验。首先,对新疆和田研究区2008年第29期和第35期的数据进行了抠值处理。其中,第29期图像抠除的范围为21×21,总计441个像元; 第35期图像抠除的范围为26×26,总计676个像元。然后分别使用3种方法对缺值区域的LST进行重建,结果如图10和图11所示。

图10 CNN-LSTM与对比方法重建效果的对比

图11 新疆和田地区LST重建精度与效率分析

观察第29期和第35期的重建结果发现,本文方法重建后的图像与原始图像的一致性更高,纹理信息更相近,而其他2种方法的重建结果相对较差。此外,为了定量比较3种方法的重建效果,本文还对重建数据与原始数据进行了比较,结果如图11所示。从图11a 和11b中可以看出,本文方法的重建误差整体较低。其中,重建误差<1K的像元数量占总缺值像元数量的一半以上,而误差>2K的像元数量明显少于其他2种方法。图11c 展示了3种方法的重建误差RMSE,从图中可以发现,本文方法的RMSE分别为0.916K和1.114K,均低于其他2种方法。从图11d 中可以看出,本文方法的重建结果与原始数据的相关性最高。图11e 展示了3种方法重建结果的平均偏差,从图中可以看出基于SCLSTM方法重建后的平均偏差最小。此外,图11f 展示了3种方法重建效率的统计分析结果,从中可以发现SSA-LSTM方法所用时间最少,本文方法次之,而SSA-BiLSTM方法所用时间最多。对于新疆和田地区,尽管本文方法所耗费的时间比SSA-LSTM略多,但在重建精度上具有明显优势。

3.2 区域适用性验证

为了证明新方法的区域适用性,我们在四川汶川地区进行了实验。首先,将2020年第8期和第44期四川汶川地区空间上部分LST数据进行去除,生成数据缺失的LST数据。第8期图像抠除的范围为21×21,总计441个像元; 第44期图像抠除的范围为26×26,总计676个像元。然后分别使用3种方法对缺值区域的LST进行重建,其重建结果如图12和图13所示。

图12 CNN-LSTM与对比方法重建效果的对比

图13 四川汶川地区LST重建精度与效率分析

从图12中观察2期图像的重建结果发现,SSA-LSTM方法和SSA-BiLSTM方法重建后的图像分别存在大量的高值和低值区,而本文方法重建后的图像与原始图像的一致性更高。此外,从图13a 和图13b 中可以看出本文方法的重建误差整体较低,其中<1K的像元数量占总缺值像元数量的一半以上,而误差>2K的像元数量明显少于其他2种方法。且本文方法2期图像的RMSE分别为0.724K和0.712K,均低于其他2种对比方法。此外,重建后的数据与原始LST的相关系数为0.95。从图13e 中可以看出,基于本文方法重建后的平均偏差最小。最后,对比3种方法的重建效率分析发现,SSA-LSTM方法所用时间最少,本文方法次之,SSA-BiLSTM方法所用时间最多(图13f),这与新疆和田的实验结果一致。通过四川汶川地区的重建实验可以看出,本文的重建方法即使是在云雾覆盖较多、天气状况较差、地表覆盖类型复杂的区域也能取得较好的数据重建效果,这证明了本文方法的可靠性和区域普适性。

3.3 地面实测数据验证

为了进一步验证新方法的可靠性,利用新疆和田6个气象站0cm地表温度的实测数据对新方法的准确性进行了评价。6个气象站分别为皮山、墨玉、和田、洛浦、策勒和于田,所有站点的空间位置信息如表2所示。其中,0cm地表温度数据来自气象台站的观测,它是指疏松、平整的土壤表面的观测温度,主要利用玻璃液体地温表测量获取。此外,气象站0cm实测地表温度数据只代表观测点局部区域的地表温度状况。本研究所用气象站实测数据的时间分辨率为1h,测量误差为0.1K。在本研究中,基于2015—2019年的MODIS LST数据的时序变化特征,重建2020年的LST缺失值,并将重建结果与实测数据进行对比分析。然而,由于MODIS 8d合成数据和气象站实测数据的时间分辨率不一致,无法直接将实测数据和MODISLST进行对比。因此,本研究采用MODIS 8d合成数据的合成方式(对8d中晴好天气下凌晨 1︰30 的LST取平均值),对气象站实测地表温度数据进行相同的处理,以保证2种数据的时间一致性。此外,为了保证验证结果的准确性,本研究只选取气象站点对应MODIS像元上的LST值与实测数据进行对比分析。

表2 气象站的空间位置

具体验证过程包括3个步骤:首先,选取6个气象站所在的6个对应像素上的MODIS LST数据,并分别去除所选像素时间域上20%的数据,生成6组数据缺失的时间序列LST。然后,基于新方法对缺失的LST数据进行重建。最后,基于MODIS LST原始数据和实测LST数据对重建精度进行验证(图14)。其中,MODIS LST原始数据与气象站实测数据的相关性如图14a 所示,R2为0.932,AD为4.003K,将其作为补值效果的判别标准。本文通过对比新方法与其他方法重建后的MODIS LST数据与实测数据的R2和AD判断补值效果的优劣。由图14b—d可知,基于SCLSTM方法重建后的LST数据和实测数据的相关系数和平均偏差(R2为0.923,AD为4.153K)与原始数据和实测数据的相关系数和平均偏差更接近。此外,图14e—g分别显示了MODIS LST原始数据及3种重建方法所得结果之间的R2、RMSE和AD。通过对比可以看出,本文所述方法的重建精度最高,R2为0.987,RMSE为1.143K,AD为1.061K。综上可知,基于本文方法重建后的LST与原始数据能够保持较好的一致性。

图14 基于气象站实测LST的可靠性分析

3.4 本文方法的实际应用

为了更直观地显示新方法在空间上的重建效果,本文还对2008年新疆和田地区和2020年四川汶川地区的LST缺值区域进行了数据重建,得到了完整的地表温度数据。图15和图16分别展示了2个研究区图像重建前后的对比效果。仔细观察重建图像不难发现,补值后的图像纹理细致自然,没有明显的边界效应,由此可以证明本文所述方法能够较好地实现大面积LST缺值区域的数据重建。

图15 新疆和田重建前后对比

图16 四川汶川重建前后对比

4 结论

MODIS LST数据产品因存在大量缺值而导致其应用受限。为获得时空连续的地表温度数据,本研究提出了一种基于混合模型(SCLSTM)的地表温度2步重建方法。该方法主要包括基于SSA方法的初步重建和基于CNN-3层堆叠LSTM方法的精细重建两大部分。利用MODIS MYD11A2遥感数据,基于本文所提新方法,采用“去除-重建-对比”的策略对新疆和田地区和四川汶川地区进行了LST重建实验,并与其他2种混合模型方法进行比较。此外,基于气象站0cm实测地表温度数据对所提新方法的LST重建精度进行了可靠性验证。主要结论如下:

(1)与现有的其他2种混合模型重建方法相比,新方法能够更好地捕捉到LST数据的时间序列特征,使得重建图像不仅能够较好地保持原始图像的纹理特征,同时还可以提升数据重建的准确性,在3种方法中重建误差最小,RMSE可降至0.712K,重建后LST与原始数据的相关系数可达0.95以上。此外,通过气象站实测数据验证进一步证明了该方法的可靠性。

(2)新方法不仅适用于缺值量少的地区,而且在云雾量大、地表覆盖类型复杂的地区也能取得较好的重建效果,证明该方法具有良好的区域普适性。综上可知,基于该方法对LST数据进行重建,可为后续地震前兆热异常信息的精确提取提供数据支撑。

虽然本文方法可以实现LST的有效重建,但仍存在以下不足:本文在时间序列预测模型LSTM的基础上加入了一维卷积神经网络1DCNN,虽然该模型取得了较好的重建效果,但由于重建模型的复杂度有所增加,因此模型的运行时间比单一的LSTM模型更长。其次,本文方法缺少对空间信息的利用,在未来的研究中可以考虑将二维空间卷积网络与时间预测模型相结合,建立一种能够同时获得地表温度时空联合特征的重建模型,以期进一步提升LST数据的重建精度。