-
地表温度(land surface temperature，LST)是研究地表能量平衡、水文变化和气候变化的重要变量[1-2]。随着航天遥感技术的发展，很多传感器可提供空间连续的地表温度产品，但现有的地表温度产品都是对其遥感影像进行云检测和云覆盖像元标识，未估算云覆盖像元地表温度。云覆盖像元地表温度信息缺失在一定程度上破坏了地表温度时空分布的全局性和连续性，严重阻碍和限制了地表温度产品的应用和发展。因此，如何估算地表温度产品中云覆盖像元地表温度，形成无缝的卫星地表温度产品是热红外遥感的研究难点。
目前，LST数据重建方法主要有基于空间相关性的空间插值重建和以时间变化规律为前提的时间序列重建方法两种。基于地表温度的空间相关性的空间插值重建方法有很多[3-4]，例如距离反比法[5]，该方法假定相邻相似像元具有相同的地表温度，但是这类方法多适用于具有较高分辨率的极轨卫星，如中分辨率成像光谱仪(moderate resolution imaging spectroradiometer，MODIS)、Landsat，当空间分辨率降低时，地表温度的空间相关性会减弱，因此此类方法不适用于空间分辨率比较低的地球静止气象卫星。此外，此类方法针对小范围的缺失数据有较好的效果，但当缺失像元的区域较大时，重建后的LST精度会有明显的降低。以时间变化规律为前提的时间序列地表温度重建方法主要是利用温度昼夜变化模型(diurnal temperature cycle，DTC)进行空值修复[6]，目前DTC模型的研究已有很大的进展[7]，在参数增加的情况下，拟合效果会相应提高，但是增加参数会导致空值修复效率下降[8]。DTC模型在时间序列上空值数量较少的情况下，LST重建效果比较好，但在空值数量较多的情况下，修复效果会明显变差[9]。当DTC模型在有效值数目少于参数数目时，无法建立模型，故而无法实现对所有LST空值的重建。
因此，本文针对静止气象卫星风云四号A星(Fengyun 4A，FY-4A)提出一种顾及时空特征的重建模型，为获取时空连续的LST无缝数据集提供了一种解决思路。该模型利用LST时序变化特征和空间相似性特征对FY-4A的LST产品数据集进行初步填补，同时对时空填补后的LST时间序列采用Savitzky-Golay(S-G)滤波器进行平滑去噪，以提高FY-4A的LST产品空值像元的重建精度和质量。
-
本文研究区域位于中国西北部的黑河流域[10]，黑河流域地形复杂，包括平原、沙漠、高原和山地，图 1为该区域的土地利用类型图。
图 1 研究区域土地利用类型图
Figure 1. Land Cover Map of Study Area
本文考虑到FY-4A的空间分辨率(4 000 m)以及验证站点在该范围内地表覆盖类型的均一性问题，采用大满观测站的观测数据作为验证数据。
-
本文所使用的数据集包括以下4种：FY-4A的LST产品数据集(FY4A_AGRI_L2_LST)、FY-4A的快速大气订正云图产品(FY4A_AGRI_L2_ACI)、MODIS的每日陆地温度发射率产品数据集(MOD11B1/MYD11B1)和地面气象站点实测数据，数据的获取时间为2018-05，数据集具体参数见表 1[11]。
表 1 数据集参数
Table 1. Parameters of Data Sets
数据名称 空间分辨率 时间分辨率 数据用途 FY4A_AGRI_L2_LST 4 000 m 15 min、60 min、不定时 FY-4A的地表温度影像数据 FY4A_AGRI_L2_ACI 1 000 m 15 min、60 min、不定时 FY-4A的大气校正反射率数据，用于合成NDVI MOD11B1/MYD11B1 6 000 m 每天 MODIS的每日陆地温度发射率数据，用于计算宽波段发射率 气象站点数据 10 min 计算实测地表温度 -
FY-4A于2016-12-11发射成功[11]，携带的可见光红外成像仪(advanced geostationary radiation imager，AGRI) 有14个光谱通道，15 min能进行一次全圆盘探测[12]。本文获取了黑河地区的FY4A_AGRI _L2_LST数据集，一共998幅圆盘数据。由于该数据不具备投影和经纬度信息，因此根据风云卫星遥感数据服务网提供FY-4A数据行列号和经纬度的互换方法，依据研究区域对数据集进行了裁剪，获取了覆盖黑河流域998幅FY-4A地表温度数据集，并统计了每幅图像的有效LST数据占比，结果如图 2所示。每幅图像有效LST数据占总像元数目的比例均值为46.24%，LST有效率低于50%的图像数目占总图像数目的53.6%。由图 2可知，由于FY-4A的LST数据集的空间分辨率较低，其受云雨影响较大，导致FY-4A的LST数据集的空间有效率比较低。
图 2 每幅图像空间有效率的柱状图
Figure 2. Histogram of Space Efficiency Ratio of Each Image
获取的不同像元LST时序数据的有效率的空间分布如图 3所示。时相有效率平均值为46.24%，最大值是61.72%，最小值是5.41%。较高的时序缺失率主要分布在黑河流域上游(研究区域的右下角)，此区域主要是高山，受云雨影响较大。在其他区域尽管LST时序的有效率相对较高，但是大部分仍然低于50%。
图 3 像元的时序数据有效率空间分布图
Figure 3. Spatial Distribution of Time-Series Data Efficiency of Pixels
-
FY4A_AGRI_L2_ACI是经过大气校正后的红外、近红外、短波红外三色通道反射率数据集，用来计算归一化植被指数(normalized vegetation index，NDVI)，利用NDVI寻找相似像元。由于受云影响，NDVI数据存在大量空值，先采用最大值合成法得到每日合成的NDVI数据集[13]，如果遇到连续阴雨天，则采用最大值合成法得到每10天或者每月合成的NDVI数据集，获得较完整连续的NDVI数据。然后采用聚合平均的方法将计算得到的NDVI降采样到4 000 m，和FY-4A的LST数据集的空间分辨率保持一致。
-
本文使用数据集中的20、22、23、29、31、32波段的发射率合成3~14 μm宽波段发射率[14-15]，获得与地面气象观测站点对应像元的宽波段发射率数据，再结合地面气象站数据估算该气象站点的实际地表温度。
-
地面气象站点实测数据获取自寒区旱区科学数据中心，实测数据用于对FY4A_AGRI_L2_LST产品重建后的精度评估。使用的观测数据集是上下行长波辐射，气象站每10 min记录一次观测数据。其中上下行长波辐射以及MODIS的发射率数据集利用Stefan-Boltzmann定律来计算实测地表温度[16]。
-
本文提出顾及时空特征的LST重建模型，该模型从LST的时间角度和空间角度提出以下两个基本假设：(1)对于某一像元的LST时序，晴空天气下邻近天相同时段的地表温度变化是相似的[17]；(2)LST具有空间连续性，同一时刻，地表覆盖物类型相似的邻近像元具有相似的LST。
本文提出的LST重建算法的技术路线图如图 4所示。首先根据LST的时间变化规律，利用邻近时间LST相似的特征，基于时间域空值重建方法，实现LST产品空值的部分重建。其次，基于空间域空值重建方法，利用NDVI寻找待重建像元的相似像元，进行空域插补。当完成第一次时空重建后，由于连续阴雨天等极端天气情况导致大面积长时序缺值无法有效地完成所有像元的重建，因此需要再次进行时空重建，直至最终完成对所有无效像元的重建。然后对完成初步重建的LST数据集采用S-G滤波器方法进行去噪处理，最终实现LST产品空值像元的重建。最后结合模拟影像和验证站点数据，对重建的LST数据集进行精度评估，并将此模型应用于黑河流域，生成了998幅2018-05的LST无缝数据集。
图 4 本文方法的技术路线图
Figure 4. Technology Flowchart of Our Proposed Method
-
对于某一像元的LST时序，晴空天气下邻近天相同时段的地表温度变化是相似的，此时时间域空值重建的计算如下：
$$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right)-\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{k}\right)= \\ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{p}\right)-\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{k}\right) $$ (1) 式中，(xp，yp)为待重建像元的空间位置；$ {d}_{p} $为待重建日期；$ {d}_{k} $为邻近日期；$ {t}_{p} $为待重建时刻；$ {t}_{k} $为邻近时刻。由此可计算待重建时刻的地表温度为：
$$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right)=\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{k}\right)+ \\ [\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{p}\right)-\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{k}\right)] $$ (2) 一个无效的LST时序序列的开始前和结束后有两个有效的时刻t0和tend，因此可以借助t0和tend时刻的有效数据对待重建时刻tp的数据分别进行重建，在重建时采用多个邻近天(即时间窗口)的相同时段进行计算。
借助t0时刻的有效数据进行重建时，计算如下：
$$ \mathrm{L}\mathrm{S}{\mathrm{T}}_{0}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right)=\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{0}\right)+ \\ \frac{1}{{N}_{t}}\sum\limits_{{d}_{k}={d}_{p}-{n}_{t}}^{{d}_{p}+{n}_{t}}[\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{p}\right)- \\ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{0}\right)] $$ (3) 式中，$ \mathrm{L}\mathrm{S}{\mathrm{T}}_{0}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right) $表示利用t0时刻的有效数据对$ {d}_{p} $日期$ {t}_{p} $时刻的空值重建结果；$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{0}\right) $表示$ {d}_{p} $日期$ {t}_{0} $时刻的有效数据；$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{p}\right) $表示$ {d}_{k} $日期$ {t}_{p} $时刻的有效数据；$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{0}\right) $表示$ {d}_{k} $日期$ {t}_{0} $时刻的有效数据；nt表示时间窗口大小的一半；$ {N}_{t} $表示在2nt+1大小的时间窗口内LST有效值的数目。
借助tend时刻的有效数据进行重建时，计算如下：
$$ \begin{array}{c}\mathrm{L}\mathrm{S}{\mathrm{T}}_{\mathrm{e}\mathrm{n}\mathrm{d}}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right)=\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{\mathrm{e}\mathrm{n}\mathrm{d}}\right)+ \\ \frac{1}{{N}_{t}}\sum\limits_{{d}_{k}={d}_{p}-{n}_{t}}^{{d}_{p}+{n}_{t}}[\mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{p}\right)-\\ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{\mathrm{e}\mathrm{n}\mathrm{d}}\right)]\end{array} $$ (4) 式中，$ \mathrm{L}\mathrm{S}{\mathrm{T}}_{\mathrm{e}\mathrm{n}\mathrm{d}}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right) $表示利用tend时刻的有效数据对$ {d}_{p} $日期$ {t}_{p} $时刻的空值重建结果；$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{\mathrm{e}\mathrm{n}\mathrm{d}}\right) $表示$ {d}_{p} $日期$ {t}_{\mathrm{e}\mathrm{n}\mathrm{d}} $时刻的有效数据；$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{p}\right) $表示$ {d}_{k} $日期$ {t}_{p} $时刻的有效数据；$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{k}, {t}_{\mathrm{e}\mathrm{n}\mathrm{d}}\right) $表示$ {d}_{k} $日期tend时刻的有效数据。
对两次预测结果采用时相加权的方法进行组合，最终的重建结果可以表示为：
$$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right)={T}_{0}\times \mathrm{L}\mathrm{S}{\mathrm{T}}_{0}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right)+\\{T}_{\mathrm{e}\mathrm{n}\mathrm{d}}\times \mathrm{L}\mathrm{S}{\mathrm{T}}_{\mathrm{e}\mathrm{n}\mathrm{d}}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right) $$ (5) $$ {T}_{0}=\frac{\left|{t}_{p}-{t}_{0}\right|}{\left|{t}_{p}-{t}_{0}\right|+\left|{t}_{p}-{t}_{\mathrm{e}\mathrm{n}\mathrm{d}}\right|} $$ (6) $$ {T}_{\mathrm{e}\mathrm{n}\mathrm{d}}=\frac{\left|{t}_{p}-{t}_{\mathrm{e}\mathrm{n}\mathrm{d}}\right|}{\left|{t}_{p}-{t}_{0}\right|+\left|{t}_{p}-{t}_{\mathrm{e}\mathrm{n}\mathrm{d}}\right|} $$ (7) 式中，$ \mathrm{L}\mathrm{S}{\mathrm{T}}_{}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right) $表示利用对$ {d}_{p} $日期$ {t}_{p} $时刻的最终的空值重建结果；$ {T}_{0} $和$ {T}_{\mathrm{e}\mathrm{n}\mathrm{d}} $分别表示$ {t}_{0} $和$ {t}_{\mathrm{e}\mathrm{n}\mathrm{d}} $时刻对$ {t}_{p} $时刻像元重建数据的时间权重。
-
LST产品的空值完成时间域重建后，部分空值由于在选定的时间窗口内无有效值，故而没有完成重建，此时采用空间域重建的方法进行LST产品空值的重建。以目标像元为中心，选取一定大小的空间窗口进行LST空值的重建。在空间域进行LST空值重建时，先利用NDVI筛选得到目标像元空间窗口内的相似像元，然后通过引入相邻相似像元的LST信息，利用加权函数完成对无效LST像元的重建，计算如下：
$$ \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{p}, {y}_{p}, {d}_{p}, {t}_{p}\right)=\sum\limits_{i=1}^{{N}_{s}}{W}_{i}\times \mathrm{L}\mathrm{S}\mathrm{T}\left({x}_{i}, {y}_{i}, {d}_{p}, {t}_{p}\right) $$ (8) $$ {W}_{i}=\frac{1/\left({D}_{i}{S}_{i}\right)}{\sum\limits_{i=1}^{{N}_{s}}[1/({D}_{i}{S}_{i}\left)\right]} $$ (9) 式中，下标p为待重建像元的位置；i为相邻有效像元的位置；$ {t}_{p} $是待重建日期；Ns是在空间窗口内的有效LST数据的总数；$ {D}_{i} $为距离因子；$ {S}_{i} $为纹理相似性因子计算的加权函数。
根据地理学第一定律，两个像素越接近，其LST值应该越相似。因此，距离越近的像素越接近中心位置，距离因子$ {D}_{i} $的计算公式为：
$$ {D}_{i}=1+\frac{\sqrt{({x}_{i}-{x}_{p}{)}^{2}+({y}_{i}-{y}_{p}{)}^{2}}}{{n}_{s}} $$ (10) 式中，ns表示空间窗口大小的一半，完整的空间窗口大小为2ns+1。
为了筛选有效像元，首先确定空间窗口大小，计算空间窗口内所有像元的NDVI的标准差，然后计算各个像元与中心像元NDVI的差值，选取差值绝对值小于标准差的对应像元作为相似像元[18-19]。选取相似像元内具有有效LST值的像元作为有效像元，参与中心像元LST的重建。纹理相似性因子$ {S}_{i} $的计算如下：
$$ {S}_{i}=\left|\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}({x}_{i}, {y}_{i})-\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}({x}_{p}, {y}_{p})\right| $$ (11) -
文献[20]提出的S-G滤波器又称最小二乘法或数据平滑多项式滤波器，是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器最大的特点是在滤除噪声的同时可以确保信号的形状、宽度不变[21]。为了最大程度发挥S-G滤波器的优势，需要进行如下设定：(1)滑动窗口的大小，若设定的窗口偏小，则会生成冗余数据，若窗口偏大，则容易遗漏信息；(2)平滑多项式的次数，若平滑次数较高，能极大程度地去掉异常值，但因过度拟合容易产生新的异常，若平滑次数较低，则不能完全去除异常值。考虑到FY-4A_LST产品的时间分辨率为1 h，在白天温度上升阶段、白天温度下降阶段以及夜间3个阶段的LST数量为8~14个，因此将滤波窗口的大小设置为11，为保证平滑效果，将平滑次数设置为2。
-
为了验证本文模型的有效性，首先通过实验确定模型的时间窗口和空间窗口的大小，然后基于站点实测数据分情况对模型的重建效果进行评价，最后分析空值区域大小对LST重建效果的影响。采用均方根误差(root mean squared error，RMSE)和决定系数R2作为重建效果的评价指标。
-
时间窗口与空间窗口大小的选取对重建结果的精度有重要影响。实验中时间窗口大小设置为3~13个像元，步长为2个像元，空间窗口大小设置为3~11个像元，步长为2个像元。
将FY-4A的LST产品的重建结果与大满站实测LST进行对比验证，图 5展示了空间窗口大小和时间窗口大小变化对重建结果的RMSE和R2的影响。由图 5可知，当空间窗口大小固定时，随着时间窗口的增加，RMSE先降低后增加；当空间窗口大小固定时，随时间窗口的增加，R2先逐步增加，后趋于稳定；当时间窗口大小为11，空间窗口大小为5时，重建精度相对最好。
图 5 时间和空间窗口大小对重建效果的影响
Figure 5. Effects of Spatial and Temporal Window Sizes Reconstruction
-
本文以黑河流域大满站作为验证站点，验证本模型的重建效果。由于验证站点部分实测数据缺失，一共采用了931个时序数据进行精度评价。验证站点实测LST数据、FY-4A的LST产品的时序数据以及重建后的LST时序数据随时间的变化曲线如图 6所示。
图 6 大满站LST重建结果
Figure 6. Result of LST Reconstruction at Daman Station
由图 6可以看出，本文模型可以实现对无效LST值的重建。与原始的FY-4A的LST数据相比，使用本模型得到的LST数据集保留了绝大部分原始的LST数据，对于由于云雨等条件影响而导致的空值，本模型也实现了重建。对于部分离散值，利用S-G滤波进行了修正。
图 7为不同空间缺失率的图像重建前后对比图。由图 7可以发现，无论对于缺失像元少的图像还是缺失像元多的图像，本文模型均可以获取无缝的LST图像。在重建过程中，先使用时间域空值重建方法和空间域空值重建方法得到初步填补的图像，此时可以看到有一些异常值，经过S-G滤波平滑优化后，异常值明显减少，说明了本模型所采用方法的有效性。
图 7 不同空间缺失率图像重建前后对比图
Figure 7. Comparison of Results Before and After Image Reconstruction with Different Space Loss Rates
图 8展现了验证站点实测数据与重建前后的LST数据的对比情况。由图 8可知，与验证站点实测数据相比，重建前FY-4A的LST时序有效数据有427个，RMSE为4.805 K，R2为0.842。重建后FY-4A的LST时序有效数据有930个，RMSE为6.969 K，R2为0.605。重建精度相比FY-4A的LST产品有所降低，经分析，主要原因如下：(1)重建的LST时序是所有天气情况下的无效LST像元，考虑到本模型的理论基础，对于阴雨天气的重建能力不足，因此大量阴雨天气的重建LST时序参与精度评估导致重建精度有所降低。(2)FY-4A的LST产品在某些时刻其精度较低，这些时刻的邻近时刻受其影响，重建后LST精度较低。而验证站点需重建的数据有503个，故而导致重建后的LST时序与验证站点相比精度降低。与FY-4A的LST产品本身的精度RMSE相比较，重建后LST的RMSE增加了2.163 K。为了具体分析本文模型在不同天气情况下的重建效果，下面进行分情况讨论。
图 8 整体重建结果的散点图
Figure 8. A Scatter Chart of the Overall Reconstruction Results
-
为了更加准确地评估本模型的重建性能，本文依据FY-4A的LST产品时序的缺失情况和天气状况，分3种情况对大满站点的LST重建结果进行讨论。情况一：晴好天气，LST时序存在部分空值；情况二：多云天气，LST时序存在大量空值；情况三：阴雨天气，LST时序存在大量空值。图 9为不同情况站点实测值与重建前后的LST数据的对比散点图。从图 9中可以看到，对于情况一，以05-12 11：00—05-13 11：00的时序数据为例，时序数据有34个，有效数据有19个。分别将FY-4A的LST产品数据、重建的LST数据与实测LST数据进行比较，得到的RMSE分别为3.392 K、5.016 K，R2分别是0.948、0.874。对于误差比较大的散点，结合图 6可知，这些点的邻近时刻的有效值，其误差也较大，故而导致重建精度下降。对于情况二，以05-08 16：00—05-09 16：00的时序数据为例，时序数据有34个，有效数据有2个。将重建的LST数据与实测LST数据进行比较，RMSE为5.053 K，R2是0.726。在连续长时序无有效的LST数据的情况下，重建LST时序数据精度与情况一相比较，RMSE增加了0.037 K，R2下降了0.148。对于情况三，以05-15 11：00—05-16 11：00的时序数据为例，时序数据有32个，有效数据有3个。将重建的LST数据与实测LST数据进行比较，RMSE为7.872 K，R2是0.313，重建精度较差，这种情况下重建精度较差与本模型的理论基础有关，在阴雨天气情况下，受降雨以及地表热通量改变的影响，地物的LST变化速率与晴好天气和多云天气相比发生了显著的变化。故而本模型在阴雨天气等极端气象条件下，精度不理想。
图 9 不同情况重建结果的散点图
Figure 9. Scatter Plots of Reconstructed Results Under Different Conditions
排除FY-4A的LST产品本身的误差对模型重建精度的影响，对于情况一，重建后LST数据的RMSE为5.016 K，FY-4A的LST产品数据的RMSE为3.392 K，二者相比较可知，晴好天气下本模型的RMSE增加了1.624 K。对于情况二，假设此时FY-4A的LST产品本身的误差与情况一是一致的，RMSE也认为是3.392 K，此种情况重建后LST数据的RMSE为5.053 K，相比之下增加了1.661 K。可见，对于多云天气情况，FY-4A的LST产品由于受大气状况的影响，其精度相对于晴好天气的精度要低，故而此处的假设是合理的。对于情况三，与情况二相同，排除FY-4A的LST产品本身的误差，可以计算得到阴雨天气下模型的重建精度RMSE增加了4.480 K。可见，由于模型的理论假设限制，此种天气情况下重建结果并不理想。综上可得，本文模型不仅在晴好条件下具有较高的空值重建精度，而且在多云天气条件且数据缺失较多的情况下，也具有良好的重建精度，但是在阴雨天气条件下的重建精度不是很理想。
-
为了测试不同空值区空间大小对本文重建模型的影响，同时排除站点的空间代表性影响[22]，在实验区随机掩膜不同大小的空间范围，将重建LST数据与原始影像的LST数据进行比较。选取实验区2018-05-26 05：00的FY-4A的LST影像，模拟10×10像素、20×20像素、30×30像素、40×40像素、50×50像素、60×60像素等6种不同大小的空值区，计算其重建效果，结果如表 2所示。
表 2 不同空值区空间大小下的LST重建精度
Table 2. LST Reconstruction Accuracy Under Different Null Space Sizes
空值区空间大小/像素 RMSE/K R2 10×10 0.484 0.846 20×20 0.491 0.896 30×30 0.485 0.923 40×40 0.483 0.938 50×50 0.492 0.957 60×60 0.507 0.976 由表 2可知，对于10×10~60×60像素不同大小的空值区，与原始的FY-4A的LST数据相比较，重建后的RMSE为0.483~0.507 K，R2为0.846~0.976。RMSE变化值小于0.03 K，决定系数R2变化值也小于0.1。说明该重建模型精度受空值区空间大小的影响相对较小。由于本次实验模拟的是一个时刻的像元缺失，主要是采用时间域空值重建方法进行LST重建，考虑到该方法是以LST时间变化规律为前提的LST重建方法，故此结论是合理的。R2平均值在空值区较小的时候较差，原因是此时数据样本很少，LST的变化范围也较小，R2的大小受样本数量影响，此时值较小。随着空值区的增加，R2逐渐趋于稳定。
-
为了测试模型在时间域上的重建性能，以晴朗日期2018-05-02的LST数据集为例(这一天有32幅LST影像，分别对应不同的时刻)，确定大满站对应的像元，去除部分时刻的LST有效值且去除的时刻是连续的，去除掉的LST时刻数量分别是1、6、11、16、21、26、31个，然后对这些模拟的空值进行重建，将重建结果与原始影像的LST数据进行比较，重建效果如表 3所示。可以看出，RMSE为0.405~1.915 K，R2为0.952~0.989，证明了本模型在时间域上进行重建的有效性。随着连续无效LST时序的增加，其RMSE逐渐增加，且变化趋于稳定。而R2随着连续无效LST时序的增加最后趋于相对稳定。与文献[23]提出的DTC模型相比，本文提出的顾及时空特征的LST空值重建方法不受一天内有效LST数据数目和分布的影响。DTC模型对于LST有效值的分布时刻有要求，即为了求解DTC模型的参数要求LST有效值均匀分布在一天内，假如有效值只分布在上午时刻或者下午时刻，则无法进行模型的求解，此外DTC模型参数的求解对有效值的数目有要求，为了求解出模型参数，其有效值不得少于6个[9]。这就意味着在一天内LST有效值数目较少或者LST有效值只分布在上午或者下午或者晚上时无法实现对无效LST像元的重建。而本模型在进行时间域重建时的理论基础是邻近天相同时刻LST的变化是相似的，对LST有效值的时序分布和有效值的数目没有要求，由实验结果可以看出，即使在LST连续空值数目为31个也可以实现对LST的重建，且精度没有发生剧烈变化。由文献[9]可知，时序上空值数目为1~16个时，使用DTC模型对空值进行重建的精度RMSE为0.266 ~0.913 K，本模型RMSE为0.405~1.634 K，在空值数目较多时，使用DTC模型无法进行空值重建。与DTC模型相比较，本模型在LST无效时序较长时依然可以重建，且具有较好的精度，在时序上空值数目为21~31个时，RMSE为1.596~1.915 K，显示了模型的稳定性。这体现了本模型的优势，即对LST的有效值的时序分布和数目没有要求，在LST无效时序较长时依然具有较好的精度。
表 3 不同空值数目下LST重建效果评估
Table 3. Evaluation of LST Reconstruction Effect Under Different Null Values
连续无效LST数目/个 RMSE/K R2 1 0.405 6 1.234 0.952 11 1.388 0.984 16 1.634 0.986 21 1.915 0.981 26 1.438 0.989 31 1.596 0.985 -
针对FY-4A静止卫星的LST产品数据存在大量空值，且空值的分布具有时空连续性的特点，本文提出了顾及时空特征的LST重建模型，实现了无效LST像元的重建，为获取FY-4A的无缝数据集提供了一种解决思路。
本文提出的顾及时空特征的LST重建模型实现了FY-4A的LST重建，本模型所用的辅助数据均来自于FY-4A卫星的产品数据，数据源易于获取，模型易实现。与DTC模型相比，本模型重建LST不受像元有效值数目和分布时刻的影响，在有效数据很少时也可以实现对无效LST的重建；不仅可以完成对晴好天气的LST空值重建，而且对于多云天气情况依旧有较好的重建精度。
Land Surface Temperature Reconstruction Model of FY-4A Cloudy Pixels Considering Spatial and Temporal Characteristics
-
摘要: 热红外遥感为地表温度(land surface temperature, LST)时空全局快速获取提供了有效手段, 但目前已有的地表温度产品未估算云覆盖像元地表温度, 如何估算地表温度产品中空值像元的地表温度, 得到无缝的地表温度数据, 是热红外遥感的研究难点。针对该问题, 提出了一种顾及时空特征的LST重建模型, 该模型首先在时间域对LST空值进行重建, 然后在空间域对LST空值进行重建, 最后采用Savitzky-Golay滤波器对重建的LST数据进行平滑去噪, 实现LST的空值重建。以黑河流域为研究区域, 以风云四号A星(Fengyun 4A, FY-4A)数据为例, 计算了该模型在不同天气条件下的重建精度, 并分析了不同空值区域大小对重建结果的影响。结果表明, 所提方法能解决晴空和多云天气下有空值像元的LST重建问题, 一天内LST连续空值数目为1~31时, 重建的均方根误差为0.405~1.915 K, 决定系数R2为0.952~0.989;与传统的昼夜温度变化模型相比较, 该模型不受有效LST像元数量和LST分布时刻的影响。Abstract:
Objectives Resorting to remote sensing technology, land surface temperature (LST) data of high spatial and temporal resolution can be conveniently acquired. However, due to the influence of bad weather conditions such as cloud and rain, there are many null values in those LST data, being an obstacle to the further application of LST data. Methods A novel null reconstruction model of LST based on spatial and temporal characteristics is proposed.Firstly, the reconstruction of partial null values of LST data is realized by weighting the effective records of the same period within contiguous days in the time domain, followed the intrinsic assumption that the changes of LST data in adjacent time periods are similar. Secondly, the reconstruction of other null valued pixels are reconstructed in the spatial domain, based on the assumption that adjacent similar pixels have similar LST. The similar pixels are searched by the normalized vegetation index (NDVI). In addition, due to the influence of extreme weather conditions, the above two steps are repeated several times until all null pixels are entirely reconstructed. Finally, Savitzky-Golay (S-G) filter is employed to remove the noise of the reconstructed LST data. Results Compared to the in-situ data, the root mean squared errors (RMSE) of the LST data of all the FY-4A before and after reconstruction are 4.805 K and 6.969 K, respectively. Their determination coefficients (R2) are 0.842 and 0.605, respectively. Under clear sky conditions, RMSEs of FY-4A LST data before and after reconstruction are 3.392 K and 5.016 K, respectively, and their R2 are 0.948 and 0.874, respectively. Under cloudy weather conditions, RMSE and R2 of reconstructed LST data are 5.053 K and 0.726, respectively. Under rainy weather conditions, RMSE and R2 of reconstructed LST data are 7.872 K and 0.313, respectively. Compared with the original LST data of FY-4A, the RMSE of the reconstructed data ranges 0.483-0.507 K and the R2 ranges 0.846-0.976 for the null value regions of different sizes in spatial domain. RMSE ranges 0.405-1.915 K and R2 ranges 0.952-0.989 for invalid LST time series of different sizes in time domain. Conclusions The proposed null reconstruction model of LST can accomplish effective reconstruction not only for clear weather, but also for cloudy weather of long time series. Excluding the error of FY-4A LST data, RMSE of reconstruction results reaches 2.171 K. When the number of LST valid pixels is very small, the invalid LST can also be effectively reconstructed by the proposed model, which is tough for the diurnal temperature cycle (DTC) model. -
图 3 像元的时序数据有效率空间分布图
Figure 3. Spatial Distribution of Time-Series Data Efficiency of Pixels
图 5 时间和空间窗口大小对重建效果的影响
Figure 5. Effects of Spatial and Temporal Window Sizes Reconstruction
图 7 不同空间缺失率图像重建前后对比图
Figure 7. Comparison of Results Before and After Image Reconstruction with Different Space Loss Rates
图 9 不同情况重建结果的散点图
Figure 9. Scatter Plots of Reconstructed Results Under Different Conditions
表 1 数据集参数
Table 1. Parameters of Data Sets
数据名称 空间分辨率 时间分辨率 数据用途 FY4A_AGRI_L2_LST 4 000 m 15 min、60 min、不定时 FY-4A的地表温度影像数据 FY4A_AGRI_L2_ACI 1 000 m 15 min、60 min、不定时 FY-4A的大气校正反射率数据，用于合成NDVI MOD11B1/MYD11B1 6 000 m 每天 MODIS的每日陆地温度发射率数据，用于计算宽波段发射率 气象站点数据 10 min 计算实测地表温度 
下载: 导出CSV
表 2 不同空值区空间大小下的LST重建精度
Table 2. LST Reconstruction Accuracy Under Different Null Space Sizes
空值区空间大小/像素 RMSE/K R2 10×10 0.484 0.846 20×20 0.491 0.896 30×30 0.485 0.923 40×40 0.483 0.938 50×50 0.492 0.957 60×60 0.507 0.976
下载: 导出CSV
表 3 不同空值数目下LST重建效果评估
Table 3. Evaluation of LST Reconstruction Effect Under Different Null Values
连续无效LST数目/个 RMSE/K R2 1 0.405 6 1.234 0.952 11 1.388 0.984 16 1.634 0.986 21 1.915 0.981 26 1.438 0.989 31 1.596 0.985
下载: 导出CSV
-
[1] Li Z L, Tang B H, Wu H, et al. Satellite- Derived Land Surface Temperature: Current Status and Perspectives[J]. Remote Sensing of Environment, 2013, 131: 14-37 doi: 10.1016/j.rse.2012.12.008 [2] 向大享, 刘良明, 韩涛. FY-3A MERSI数据干旱监测能力评价[J]. 武汉大学学报·信息科学版, 2010, 35(3): 334-338 http://ch.whu.edu.cn/article/id/888 Xiang Daxiang, Liu Liangming, Han Tao. Estimation of Drought Monitoring Ability of FY-3A MERSI Data[J]. Geomatics and Information Science of Wuhan University, 2010, 35(3): 334-338 http://ch.whu.edu.cn/article/id/888 [3] Ke L, Ding X, Song C. Reconstruction of Time-Series MODIS LST in Central Qinghai-Tibet Plateau Using Geostatistical Approach[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(6): 1 602-1 606 doi: 10.1109/LGRS.2013.2263553 [4] 颜金彪, 段晓旗, 郑文武, 等. 顾及空间异质性的自适应IDW插值算法[J]. 武汉大学学报·信息科学版, 2020, 45(1): 97-104 doi: 10.13203/j.whugis20180213 Yan Jinbiao, Duan Xiaoqi, Zheng Wenwu, et al. Geomatics and Information Science of Wuhan University[J]. Geomatics and Information Science of Wuhan University, 2020, 45(1): 97-104 doi: 10.13203/j.whugis20180213 [5] 周义, 覃志豪, 包刚. GIDS空间插值法估算云下地表温度[J]. 遥感学报, 2012, 16(3): 492-504 https://www.cnki.com.cn/Article/CJFDTOTAL-YGXB201203006.htm Zhou Yi, Qin Zhihao, Bao Gang. Land Surface Temperature Estimation Under Cloud Cover with GIDS[J]. Journal of Remote Sensing, 2012, 16(3): 492-504 https://www.cnki.com.cn/Article/CJFDTOTAL-YGXB201203006.htm [6] Göttsche F M, Olesen F S. Modelling the Effect of Optical Thickness on Diurnal Cycles of Land Surface Temperature[J]. Remote Sensing of Environment, 2009, 113(11): 2 306-2 316 doi: 10.1016/j.rse.2009.06.006 [7] Duan S B, Li Z L, Wang N, et al. Evaluation of Six Land-Surface Diurnal Temperature Cycle Models Using Clear-Sky in Situ and Satellite Data[J]. Remote Sensing of Environment, 2012, 124: 15-25 doi: 10.1016/j.rse.2012.04.016 [8] 刘紫涵, 吴鹏海, 吴艳兰, 等. 风云静止卫星地表温度产品空值数据稳健修复[J]. 遥感学报, 2017, 21(1): 40-51 https://www.cnki.com.cn/Article/CJFDTOTAL-YGXB201701004.htm Liu Zihan, Wu Penghai, Wu Yanlan, et al. Robust Reconstruction of Missing Data in Feng Yun Geostationary Satellite Land Surface Temperature Products[J]. Journal of Remote Sensing, 2017, 21(1): 40-51 https://www.cnki.com.cn/Article/CJFDTOTAL-YGXB201701004.htm [9] 薛兴盛, 吴艳兰. 面向风云静止卫星地表温度产品的缺失数据修复方法对比[J]. 安徽农业大学学报, 2017, 44(2): 308-315 https://www.cnki.com.cn/Article/CJFDTOTAL-ANHU201702022.htm Xue Xingsheng, Wu Yanlan. A Comparison of Missing Data Reconstruction Methods for Feng Yun Geostationary Satellite Land Surface Temperature Products[J]. Journal of Anhui Agricultural University, 2017, 44(2): 308-315 https://www.cnki.com.cn/Article/CJFDTOTAL-ANHU201702022.htm [10] 李国亮. 2000—2012年黑河流域上游植被覆盖变化遥感监测与分析[D]. 兰州: 西北师范大学, 2015 Li Guoliang. Remote Sensing Monitoring and Analysis of Vegetation Cover Change in the Heihe River Basin from 2000—2012[D]. Lanzhou: Northwest Normal University, 2015 [11] 张志清, 董瑶海, 丁雷, 等. 我国首颗第二代静止气象卫星风云-4升空[J]. 国际太空, 2016, 12: 6-12 https://www.cnki.com.cn/Article/CJFDTOTAL-GJTK201612003.htm Zhang Zhiqing, Dong Yaohai, Ding Lei, et al. China's First Second-Generation FY-4 Meteorological Satellite Launched[J]. Space International, 2016, 12: 6-12 https://www.cnki.com.cn/Article/CJFDTOTAL-GJTK201612003.htm [12] 谭仲辉, 马烁, 韩丁, 等. 基于随机森林算法的FY-4A云底高度估计方法[J]. 红外与毫米波学报, 2019, 38(3): 381-388 https://www.cnki.com.cn/Article/CJFDTOTAL-HWYH201903020.htm Tan Zhonghui, Ma Shuo, Han Ding, et al. Estimation of Cloud Base Height for FY-4A Satellite Based on Random Forest Algorithm[J]. Journal of Infrared and Milmillimeter Waves, 2019, 38(3): 381-388 https://www.cnki.com.cn/Article/CJFDTOTAL-HWYH201903020.htm [13] Holben B N. Characteristics of Maximum-Value Composite Images from Temporal AVHRR Data[J]. International Journal of Remote Sensing, 1986, 7(11): 1 417-1 434 doi: 10.1080/01431168608948945 [14] Snyder W C, Wan Z, Zhang Y, et al. Thermal Infrared (3-14 μm) Bidirectional Reflectance Measurements of Sands and Soils[J]. Remote Sensing of Environment, 1997, 60(1): 101-109 doi: 10.1016/S0034-4257(96)00166-6 [15] Snyder W C, Wan Z. BRDF Models to Predict Spectral Reflectance and Emissivity in the Thermal Infrared[J]. IEEE Transactions on Geoscience and Remote Sensing, 1998, 36(1): 214-225 doi: 10.1109/36.655331 [16] Yang Y, Cao C, Pan X, et al. Downscaling Land Surface Temperature in an Arid Area by Using Multiple Remote Sensing Indices with Random Forest Regression[J]. Remote Sensing, 2017, 9(8): 789-795 doi: 10.3390/rs9080789 [17] Zhu X, Chen J, Gao F, et al. An Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model for Complex Heterogeneous Regions[J]. Remote Sensing of Environment, 2010, 114(11): 2 610-2 623 doi: 10.1016/j.rse.2010.05.032 [18] Kustas W P, Norman J M, Anderson M C, et al. Estimating Subpixel Surface Temperatures and Energy Fluxes from the Vegetation Index-Radiometric Temperature Relationship[J]. Remote Sensing of Environment, 2003, 85(4): 429-440 doi: 10.1016/S0034-4257(03)00036-1 [19] Gao F, Masek J, Schwaller M, et al. On the Blending of the Landsat and MODIS Surface Reflectance: Predicting Daily Landsat Surface Reflectance[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(8): 2 207-2 218 doi: 10.1109/TGRS.2006.872081 [20] Savitzky A, Golay M J E. Smoothing and Differentiation of Data by Simplified Least Squares Procedures[J]. Analytical Chemistry, 1964, 36(8): 1 627-1 639 doi: 10.1021/ac60214a047 [21] 贾明明, 王宗明, 张柏, 等. 综合环境卫星与MODIS数据的面向对象土地覆盖分类方法[J]. 武汉大学学报·信息科学版, 2014, 39(3): 305-310 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201403012.htm Jia Mingming, Wang Zongming, Zhang Bai, et al. Land Cover Classification of Compositing HJ-1 and MODIS Data Based on Object-based Method[J]. Geomatics and Information Science of Wuhan University, 2014, 39(3): 305-310 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201403012.htm [22] 马晋, 周纪, 刘绍民, 等. 卫星遥感地表温度的真实性检验研究进展[J]. 地球科学进展, 2017, 32(6): 615-629 https://www.cnki.com.cn/Article/CJFDTOTAL-DXJZ201706006.htm Ma Jin, Zhou Ji, Liu Shaomin, et al. Review on Validation of Remotely Sensed Land Surface Temperature[J]. Advances in Earth Science, 2017, 32(6): 615-629 https://www.cnki.com.cn/Article/CJFDTOTAL-DXJZ201706006.htm [23] Van den Bergh F, Van Wyk M A, Van Wyk B J, et al. Comparison of Data-Driven and Model-Driven Approaches to Brightness Temperature Diurnal Cycle Interpolation[C]//The 17th Annual Symposium of the Pattern Recognition Association of South Africa, Parys, South Africa, 2006 -
点击查看大图
图(9) / 表(3)
计量
- 文章访问数: 558
- HTML全文浏览量: 236
- PDF下载量: 53
- 被引次数: 0
