EN中文

地球信息科学学报 >

2022 , Vol. 24 >Issue 2: 378 - 390

DOI: https://doi.org/10.12082/dqxxkx.2022.210257

遥感科学与应用技术

基于DCT-PLS算法的MODIS LST缺值填补方法研究

  • 刘恒孜 , 1, 2 ,
  • 吕宁 , 1, 3, * ,
  • 姜侯 1 ,
  • 姚凌 1, 3
展开
  • 1.中国科学院地理科学与资源研究所 资源与环境信息系统国家重点实验室,北京 100101
  • 2.中国科学院大学,北京 100049
  • 3.南方海洋科学与工程广东省实验室（广州）,广州 511458
*吕 宁（1981—）,男,河南洛阳人,副研究员,硕士生导师,主要从事大气遥感应用研究。E-mail:

刘恒孜（1997— ）,男,河南许昌人,硕士生,主要从事遥感地学分析与计算研究。E-mail:

收稿日期: 2021-05-10

  要求修回日期: 2021-07-05

  网络出版日期: 2022-04-25

基金资助

国家自然科学基金项目(41971312)

南方海洋科学与工程广东省实验室(广州)人才团队引进重大专项(GML2019ZD0301)

版权

版权所有,未经授权,不得转载、摘编本刊文章,不得使用本刊的版式设计。
收起

Research on Gaps Filling of MODIS LST based on DCT-PLS

  • LIU Hengzi , 1, 2 ,
  • Lü Ning , 1, 3, * ,
  • JIANG Hou 1 ,
  • YAO Ling 1, 3
Expand
  • 1. State Key Laboratory of Resources and Environmental Information Systems, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China
  • 2. University of Chinese Academy of Sciences, Beijing 100049, China
  • 3. Southern Marine Science and Engineering Guangdong Laboratory, Guangzhou 511458, China
*Lü Ning, E-mail:

Received date: 2021-05-10

  Request revised date: 2021-07-05

  Online published: 2022-04-25

Supported by

National Natural Science Foundation of China(41971312)

Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou)(GML2019ZD0301)

Copyright

Copyright reserved © 2022
Fold

摘要

受云层、传感器误差等因素影响,中分辨率成像光谱仪（MODIS）获取的地表温度产品（LST）在时间和空间上存在大量缺失数据,严重影响了基于时间序列数据的分析与应用。本文引进一种基于离散余弦变换与惩罚最小二乘的多维数据快速平滑方法（DCT-PLS）,利用数据集自身的时间和空间信息填补缺值。本文在粤港澳大湾区开展实证研究,将DCT-PLS算法用于填补该地区2001年1月—2017年12月的月值MODIS LST数据缺值,并引入人工模拟缺值对算法进行误差分析与精度验证。算法误差分析结果表明,填补误差主要来源于三维算法对数据集中有偏LST时间信息的使用,并因此产生显著高估或低估的填补结果。基于此,本文提出了利用MODIS LST数据集自身时间序列信息自动计算获取有效辅助LST信息的优化策略,从而实现填补算法计算效率和精度的提升：平均计算时间从12.0 s提高至1.7 s,平均R从0.94提高至0.97,平均RMSE从1.94 K提高至0.74 K（相较于三维算法）。在大湾区的填补结果表明（日间结果：R=0.98、RMSE=0.79 K;夜间结果：R=0.99、RMSE=0.56 K）,优化后的DCT-PLS算法可以快速鲁棒地填补MODIS LST月值数据产品中的缺值,并且具有稳定性强、不依赖外部数据集的计算特性,能够适应长时间序列MODIS LST缺值填补。

关键词： MODIS; 缺值填补; 地表温度; 长时序; 离散余弦变换; 惩罚最小二乘; 算法优化; 粤港澳大湾区

本文引用格式

刘恒孜 , 吕宁 , 姜侯 , 姚凌 . 基于DCT-PLS算法的MODIS LST缺值填补方法研究[J]. 地球信息科学学报, 2022 , 24(2) : 378 -390 . DOI: 10.12082/dqxxkx.2022.210257

Abstract

Affected by cloud layer, sensor error, and other factors, there are large area of data gaps in time and space of Land Surface Temperature (LST) products acquired by Moderate-resolution Imaging Spectroradiometer (MODIS), which seriously affects the analysis and application of LST data set. In this paper, a fully automatic smoothing algorithm for filling MODIS LST data sets is introduced. Based on a Penalized Least Squares (PLS) method, the algorithm allows fast smoothing of data in one and more dimensions by means of the Discrete Cosine Transform (DCT). By optimizing smoothing parameters, the DCT-PLS algorithm is a fast multidimensional data smoothing method, which uses the spatial and temporal information of LST data set to fill data gaps. This study carries out an empirical study in Guangdong-Hong Kong-Macao Greater Bay Area (GBA). The DCT-PLS algorithm is used to fill data gaps of monthly MODIS LST data set from January 2001 to December 2017 in this area. To analyze its performance, the error analysis and accuracy verification of the algorithm are performed by introducing artificial simulation data gaps. The results of error analysis show that errors in gap filling are mainly due to the use of biased LST time information in data set, which causes significantly overestimation or underestimation in filling results. In response to this problem, an optimized DCT-PLS algorithm is proposed to automatically obtain effective auxiliary LST layers from the time series of MODIS LST and provide unbiased temporal information. Meanwhile, the three-dimensional algorithm is changed to two-dimensional space to reduce computing consumption by transforming auxiliary temporal information into spatial information. The results of accuracy verification on the optimized DCT-PLS algorithm show that the computational efficiency and accuracy are optimized simultaneously. The average computing time decreases from 12.0s to 1.7s, average R increases from 0.94 to 0.97, and average RMSE decreases from 1.94 K to 0.74 K (compared with the three-dimensional algorithm). The optimized DCT-PLS algorithm is applied in GBA. The accuracy of gap filling of monthly MODIS LST data set from January 2001 to December 2017 is R=0.98 and RMSE=0.79 K in daytime, and R=0.99 and RMSE=0.56 K in nighttime. The results of accuracy verification in GBA show that the accuracy of gap filling is not affected by heterogeneous surface. The optimized DCT-PLS algorithm transforms the MODIS LST data in the time domain and frequency domain to two-dimensional space, and retains high-frequency information to fill data gaps. It is a fast and robust method for gap filling, which is completely independent of external auxiliary data sets, enabling gap filling on long-time sequence MODIS LST data set.

Key words： MODIS; gap filling; Land Surface Temperature (LST); long time sequence; discrete cosine transform; penalized least-squares; algorithm optimization; Guangdong-Hong Kong-Macao Greater Bay Area (GBA)

1 引言

陆地表面温度（Land Surface Temperature,LST）反映了地表大气相互作用和大气与地面之间的能量通量,是地表物理过程的重要参数之一,在水文学、气象学、生态学、流行病学和能源系统等各个研究领域有着广泛应用[1,2,3,4]。气象站点可以获得准确的LST数据,但存在站点稀疏、空间分布不规则、地面勘测难度较大以及点数据内插复杂性等问题[5]。基于中分辨率成像光谱仪仪器（Moderate-resolution Imaging Spectroradiometer,MODIS）卫星影像反演的LST产品克服了站点观测的不足,能够详细地捕获不同气象条件、不同局部气候和多种土地覆盖类型的温度差异,在生态环境以及农业研究中得到了广泛应用[6,7,8]。
然而,多变的大气条件、传感器故障和轨道间的扫描间隙经常将缺失数据引入到LST产品 中[9,10,11],这极大地降低了数据产品的可用性,并给后续研究与应用造成困难。目前对LST的填补方法可以分为对晴空条件下缺失像素的填补方法与对云像素的填补方法。其中,对晴空条件下缺失像素的填补主要通过数据集周围的无云LST值对缺值进行估计,按照所参考LST信息的差异可以分为空间填补方法[12,13,14,15]、时间填补方法[16,17,18,19]与时空填补方法[20,21,22]。空间填补方法一般不添加额外的辅助信息而直接利用自身数据进行缺值估计,并基于缺失数据和剩余有效数据拥有相同统计或几何结构的假设填补缺值,这些方法通常易于实现,但填补结果在异质性表面上表现出模糊的平滑效果。时间填补方法从同一空间位置的不同时间序列中获取互补信息,这些方法在规律变化的空间表面上是合适的,但由于忽略了邻近像素点的数据信息,填补能力受长时间序列中人工景观变化的影响较大。时空填补方法从额外的时空信息中同时提取互补信息,最大限度地利用所有可用的数据信息来提高填补精度,但填补模型对辅助数据集依赖程度较高,且模型参数较多。以上这些方法基于多种技术增强了晴空条件下MODIS LST数据产品的可用性,但其建模过程通常包含其他经验模型、外源LST数据集（例如,其他卫星观测、无线电探空测量和再分析产品等）或其他高分辨率辅助数据集（例如,太阳辐射、高程、植被指数和气象因子等）。外源信息的输入给填补算法引入了新的不确定性,同时,较多的参数与辅助信息也增加了填补算法在其他地区或长时间序列中迁移应用的难度。
针对以上这些填补方法的不足,找到一种快速鲁棒、独立于外部辅助数据集、且能够适应长时间序列MODIS LST的缺值填补算法十分重要。 Garcia等[23]基于离散余弦变换与惩罚最小二乘法（Penalized Least-Squares of the Discrete Cosine Transform,DCT-PLS）开发了一种多维数据平滑算法,该算法根据数据集本身在二维或三维上的有效信息填补缺失值,具有良好的计算效率,且不依赖于其他外部数据集;此外,其计算过程受模型内唯一的平滑参数控制,具有自动化迭代收敛的计算特点[24]。目前,这种方法已经在地表温度[25,26]、土壤湿度[27,28,29]、植被指数[30,31]、沿海海流[32,33]、臭氧[34]等地球物理数据集的缺值填补上进行了探索性的应用,但应用于填补长时间序列MODIS LST的研究较少,且填补能力与填补误差有待于进一步研究。本文以粤港澳大湾区作为核心研究区,应用基于DCT-PLS的平滑算法进行长时间序列MODIS LST月值数据产品的缺值填补,并分析该算法在填补LST缺值时的误差来源;在此基础上,结合二维与三维填补方法的优势对DCT-PLS算法进行优化,进一步提升算法的填补速度与填补精度。

2 研究方法

2.1 填补流程与验证方法

本文使用DCT-PLS算法填补MODIS LST月值数据集中缺值的流程如图1所示,包括：① LST数据集输入（包含原始缺值）;② “人工模拟缺值”处理（产生模拟缺值）;③ DCT-PLS算法模型填补缺值（原始缺值与模拟缺值被同时填补）;④ 填补后LST数据集输出;⑤ 填补结果精度验证。
图1 DCT-PLS算法填补流程

Fig. 1 DCT-PLS filling flow chart

由于整个填补过程只在MODIS LST月值数据集上进行,在“人工模拟缺值”步骤中,设计了包括“固定选取”与“随机选取” 2种缺值生成策略对算法填补精度进行验证分析。其中,“固定选取”生成策略在固定区域生成位置、大小、数量相同的人工模拟缺值,主要用来对算法进行误差分析与优化研究;而“随机选取”生成策略在待填补LST原始缺值周围随机生成一定数量的人工模拟缺值,主要用来对填补算法在目标LST的实际填补精度进行评价。2种策略在生成的人工模拟缺值时,均将非空的原始LST值替换为空值,当这些模拟空值被填补后,即形成了与其对应位置原始LST值的“填补LST-原始LST”像素对,并采用Pearson相关系数（R）、均方根误差（RMSE）与平均绝对误差（MAE）计算填补精度。
R = cov ( LST , LST ˆ ) σ LST σ LST ˆ
RMSE = ∑ i = 1 p LST - LST ˆ ) 2 p
MAE = ∑ i = 1 p ( | LST - LST ˆ | ) p
式中： LST代表MODIS地表温度原始值,其标准差为 σ LST; p为用于计算的样本总数; LST ˆ为地表温度的填补值,其标准差为 σ LST ˆ。

2.2 DCT-PLS算法原理

DCT-PLS算法从一维简单向量的平滑衍生而来,目标是实现对一维向量离群值的平滑与去噪,同时保留数据集中主要信息不变。假设 y代表原始的LST数据。那么,一维数据的平滑模型如下：
y = y ˆ + ϵ
式中： ϵ为均值等于0的高斯噪声; y ˆ是平滑处理后的LST数据,通过求解 y ˆ的最优估计来实现对 y的平滑。DCT-PLS算法采用非参数回归惩罚最小二乘法进行平滑运算,其最小化的目标是：
F y ˆ = RSS + sP y ˆ = y ˆ - y 2 + sP y ˆ
式中： 表示欧几里得范数;参数 s是一个正实数,用于控制平滑程度; y ˆ的平滑程度随着 s的增大而增大,实现对估计结果平滑效果的参数化控制。DCT-PLS算法采用残差平方和RSS（Residual Sum-of Squares）来度量平滑数据的保真度,并采用粗糙度 P y ˆ表达数据的平滑程度,其中 P y ˆ定义为：
P y ˆ = D y ˆ 2
式中： D为预定义的三对角方阵,根据式（5）与式（6）得到原始LST数据与平滑后LST数据的线性关系：
y = ( I n + s D T D ) y ˆ
式中： I n为 n × n的单位矩阵; D T为 D的转置;式 ( I n + s D T D )是一个对称的五对角阵。MODIS LST为均匀间距的栅格数据, y ˆ i与 y ˆ i + 1的距离 h i恒等于1,矩阵 D可被重写为关于i的简单差分矩阵：
D = - 1 1 1 - 2 1 1 - 2 ⋱ - 2 1 1 - 2 1 1 - 1
对矩阵 D进行特征分解,得到 D = UΛ U - 1。 U为酉矩阵,存在 U - 1 = U T, U U T = I n, Λ是包含 D特征值的对角矩阵, Λ = diag λ 1 , … , λ n, λ i = - 2 + 2 cos i - 1 π / n。代入式（7）可得到：
y ˆ = U ( I n + s Λ 2 ) - 1 U T y = UΓ U T y
式中： U T为 n × n的离散余弦变换（DCT）; U为逆离散余弦变（IDCT）; Γ为包含平滑参数 s的对角矩阵。因此,平滑结果 y ˆ可以写为：
y ˆ = UΓ U T y = IDCT ( Γ DCT ( y ) )
在最终的平滑方程（式（10））中,原始数据先后经过DCT变换与IDCT变换实现平滑,且平滑程度只受 s值的影响与控制。因此,通过获取平滑参数 s的最优估计,就可以避免过度平滑或平滑不足带来的误差,从而得到与原始值最接近的平滑结果。在DCT-PLS算法中,Garcia提出利用广义交叉验证法（Generalized Cross Validation,GCV）进行优化,自动迭代计算最小化平滑误差,从而选择合适的平滑参数 s[23]。根据DCT和IDCT的运算性质,在数据空间分布均匀的情况下,GCV评分为：
GCV (s) = n ∑ i - 1 n 1 1 + s λ i 2 - 1 2 DC T i 2 ( y ) n - ∑ i - 1 n 1 1 + s λ i 2 2
通过计算GCV评分,DCT-PLS的计算复杂度仅为 O ( nlog ( n ) )。在进行平滑处理时,通过给缺值赋一个较小的权重,给有效数据赋一个较高的权重,即可达到平滑运算的同时完成对缺失数据的填补。此外,通过在一维向量平滑基础上进行多维扩展,将DCT与IDCT在数据的多个维度进行运算,DCT-PLS算法可以直接应用于多维数据的平滑与填补。定义 W为对角矩阵 diag ( I n ),根据 y i中缺失像素的位置, w ∈ [ 0,1 ]分别对应于缺失像素与有效像素,式（10）可以扩展到二维甚至N维数据中：
y ˆ { k + 1 } = IDCTN ( Γ T · DCTN ( W ( y - y ˆ { k } ) + y ˆ { k } ) )
MODIS LST数据的时空填补对应于N=2或 N=3的情形,输入含有缺值的单层或多层LST数据集 y后,DCT-PLS算法通过迭代计算合适的平滑参数 s,自动输出填补后完整连续的单层或多层LST数据集 y ˆ。

3 研究区概况与数据来源

3.1 研究区概况

本文的研究区域是粤港澳大湾区（Guangdong-Hong Kong-Macao Greater Bay Area,GBA）。粤港澳大湾区位于111.21°E—114.53°E,21.27°N—24.24°N,是世界级城市群,包括珠三角城市群的 9个市（广州、佛山、肇庆、深圳、东莞、惠州、珠海、中山、江门）,以及香港、澳门2个特别行政区。大湾区属亚热带气候,终年温暖湿润,年平均气温在 21~23 ℃,地势起伏较大,中部为平原,山地主要分布在肇庆市、惠州市及广州市从化区等北部边缘地带。大湾区位于珠江支流,区内主要河道有102条、总长度约1700 km,水网纵横交错,对大湾区内的生态环境起到调节作用。大湾区近10年来快速发展,城市化进程加快,城市人口的增加与下垫面类型的变化对LST产生了明显的影响;为推动大湾区的可持续发展,需要合理规划布局城市群建设。在此背景下,城市发展造成的局部LST变化及其对区域内生态环境产生的潜在影响成为了重要研究内容,获取长时间序列、空间连续的LST数据产品是开展相关研究的必要前提。本文研究基于DCT-PLS算法的MODIS LST数据缺值填补方法,旨在满足粤港澳大湾区LST相关研究中的数据需求。

3.2 研究数据及缺值统计

本文采用美国航空航天局（National Aeronautics and Space Administration,NASA）提供的MODIS LST数据产品,该产品具有高时间分辨率（每日在当地时间约01:30、10:30、13:30和22:30进行4次观测）、高空间分辨率（1 km分辨率）和全球覆盖范围的特点。本研究使用MOD11C3产品,包括每月日间LST与夜间LST,数据基于Terra卫星反演的每日LST数据产品（MOD11C1）合成。本文获取了大湾区（实际研究区为21°N-25°N、111°E-116°E范围内区域）2001年1月—2017年12月（共204个月）的日间与夜间LST数据集,以及用于统计分析的2001年—2017年土地利用数据（来自MODIS土地利用分类产品,MCD12Q1）,数据产品来自 https://ladsweb.modaps.eosdis.nasa.gov/search/。
在研究时段内,MODIS LST产品在大湾区的累计缺值数量空间分布如图2所示。由图可见,该产品日间LST与夜间LST在大湾区均存在大范围的缺值,占该地区陆地总面积的70%以上,空间上基本覆盖整个大湾区地区。日间缺值主要分布在大湾区中部与西北部地区,夜间缺值主要分布在大湾区中部与东部地区,且日间由于更大范围的云层遮挡而产生了更多的缺值。缺失值的大范围分布将大幅降低MODIS LST月值数据集在大湾区的可用度。在表1中,各种土地利用类型的累积缺失面积占比均达到了73%以上,其中,由于不透水层对水分的大量蒸发以及靠近海岸线带来的更多降水,城市区域更大范围与更频繁的云层覆盖导致了更多缺值的产生,日间LST产品在城市区域的缺值面积占比超过了90%,而林地与灌木林地中缺值比例相对较低,且所有土地类型的日间缺值面积比例均大于夜间缺值面积比例。
图2 粤港澳大湾区MODIS LST累计（2001年1月—2017年12月）缺值数量空间分布

注：蓝色线为叠加的粤港澳行政区划界限。该图基于自然资源部标准地图服务网站下载的审图号为GS(2019)4342号的标准 地图制作,底图无修改。

Fig. 2 Spatial distribution of cumulative missing pixels of MODIS LST in GBA from 2001.1 to 2017.12

表1 粤港澳大湾区不同土地利用类型MODIS LST累计（2001年1月—2017年12月）缺值面积占比

Tab. 1 Proportion of cumulative missing area of MODIS LST for different land covers in GBA from 2001.1 to 2017.12 (%)

土地利用类型 缺失面积占比（日间） 缺失面积占比（夜间）
林地 83.12 76.18
灌木林 83.65 73.08
草地 85.71 82.14
耕地 87.14 87.14
城市区域 90.95 87.65
如图3所示,日间缺值主要分布在2001、2005、2010和2015年,夜间缺值主要分布在2001、2005和2014年,日间与夜间缺值较多的年份基本一致,日间LST由于每月更大范围的云层覆盖在缺值数量上略高于夜间LST。此外,日间与夜间缺值数量在时间上没有表现出明显的分布规律,个别月份受到大范围云遮挡或特殊天气条件的影响产生明显较多的缺值;而在整个时间序列中,大部分月份仅存在较少数量的缺值。综合大湾区缺值的时空分布特征,日间LST与夜间LST数据集均存在不同位置、不同数量的缺失值,本文将基于DCT-PLS算法对这些缺值进行填补。
图3 粤港澳大湾区MODIS LST时间序列（2001年1月—2017年12月）缺值像素个数统计图

Fig. 3 Monthly quantity of missing pixels of MODIS LST in GBA from 2001.1 to 2017.12

4 误差分析与算法优化

4.1 误差分析

DCT-PLS算法是一种全三维的时空填补算法,根据其多维计算特性,DCT-PLS算法既可以对二维LST图像进行填补（2-D算法）,也可以对三维LST时空数据进行填补（3-D算法）。对于MODIS LST月值数据集,2-D算法利用每月LST的空间信息进行缺值填补,在2个维度上进行DCT与IDCT变换运算,具有更高的计算效率,但其缺乏对数据时间信息的使用;3-D算法通过在目标LST前后添加来自时间序列数据集中其他月份的辅助LST层,同时利用目标LST的空间信息与来自辅助LST层的时间信息进行填补,可以取得比2-D算法更高的填补精度,但3-D算法需要更大的计算时间成本,并可能会引入额外的误差。如图4所示,散点图显示了 3-D算法在2002年1月与2012年4月的填补结果,虽然目标LST填补结果的R值较高,但RMSE与MAE结果较差,表现出明显的高估与低估现象。
图4 3-D算法LST缺值填补结果

Fig. 4 LST filling results of 3-D algorithm

在3-D算法填补过程中,时间信息的添加可达到提供额外辅助信息从而提高填补精度的目的,但由辅助LST层提供的时间信息可能会是“有偏”的,这是导致填补结果出现高估与低估的主要原因。本文提取了3-D算法在填补2002年1月和2012年4月LST时使用的辅助LST层（图5）。在2002年1月的辅助LST层中,2001年11月、2002年2月和2002年3月的LST明显高于2002年1月的原始LST,这会给填补算法提供偏高的时间信息;在2012年4月的辅助LST层中,2012年2月和2012年3月的LST均明显低于2012年4月的原始LST,这会给填补算法提供偏低的时间信息。根据DCT-PLS的算法原理, 3-D算法会基于以上这些具有明显偏差的辅助LST信息进行平滑运算,运算结果将会使填补LST值接近辅助LST层,并最终产生如图4所示的高估与低估的填补结果。因此,本文提出以下假设：通过降低或避免有偏辅助LST层所带来的影响,将改善3-D算法产生的高估或低估现象并显著提高填补精度。
图5 3-D算法填补缺值使用的辅助LST层

注：该图基于自然资源部标准地图服务网站下载的审图号为GS(2019)4342号的标准地图制作,底图无修改。

Fig. 5 The auxiliary LST layers used by the 3-D algorithm to fill the missing values

为了验证上述猜想,本文对3-D算法填补2002年1月和2012年4月缺值的辅助LST层进行人工控制,在数据集中选取与待填补目标LST接近的LST作为辅助LST层：对2002年1月选取2007年1月与2015年12月LST;对2012年4月选取了2002年4月与2007年4月LST后重新进行填补。填补结果如 图6所示,相比图4填补结果,经人工控制辅助LST层后的填补精度明显提升,散点图分布在1:1直线附近,“有偏”辅助LST层导致的整体高估与低估现象得到很大改善。因此,3-D算法的填补误差主要来源于填补模型对于“有偏”辅助LST信息的使用,在原始时间序列中寻找与待填补目标LST接近的LST作为辅助LST层可有效改善整体高估或低估问题。
图6 3-D算法LST缺值填补结果（人工控制辅助LST层）

Fig. 6 LST filling results of 3-D algorithm (based on the manually controlled auxiliary LST layers)

4.2 算法优化

基于以上误差分析结果,本文提出一种DCT-PLS算法优化策略：通过对每一个待填补目标LST均在数据集中自动搜索与其接近的LST层来构建辅助LST层,从而达到在整个MODIS LST月值数据集上提高填补精度的目的。同时,将辅助LST层所提供的时间信息通过加权初始化的方式转化为目标LST可用的空间信息,对3-D算法进行“降维”处理,从而结合2-D算法与3-D算法的优势,实现填补速度与填补精度的优化改进,这种“优化算法”的填补流程如图7所示。首先,对每一层待填补目标LST,计算其与数据集中所有LST的平均绝对距离（Mean Absolute Distance,MAD）,MAD定义为目标LST与其他LST非空像素上LST值的平均绝对偏差;然后根据MAD的大小进行升序排列,得到除目标LST自身外,MAD最小的前4层LST作为辅助LST层,并按照目标LST缺值范围提取对应范围内有效的辅助LST信息;其次,将得到的辅助LST信息按照权重w进行加权平均产生用于初始化的空间信息,权重w定义为每一层辅助LST层MAD倒数值在4层辅助LST的MAD倒数值之和的占比;最后将初始化后的目标LST输入2-D算法完成填补。
图7 “优化算法”的缺值填补流程

Fig. 7 Process of optimized DCT-PLS to fill the missing values of MODIS LST

5 填补结果与分析

利用“优化算法”填补2002年1月和2012年4月LST缺值,结果如图8所示。通过在MODIS LST数据集中自动选取与待填补LST接近的有效辅助LST层,“优化算法”没有产生明显的高估与低估,相对于图6,数据点更为集中地分布在1:1直线附近,精度评价指标显示填补精度得到显著提升。
图8 “优化算法”LST缺值填补结果

Fig. 8 LST filling results of optimized DCT-PLS

本研究使用2-D、3-D以及“优化算法”对整个MODIS LST时间序列进行填补精度的对比验证,结果如图9所示。2-D算法对目标LST的填补效率更高,且在整个时间序列上的RMSE和MAE波动更小,体现了算法快速稳定的计算特征; 3-D算法虽然RMSE与MAE的平均值高于2-D算法,但是填补精度在整个时间序列上波动较大,平均填补精度受填补误差较大结果的影响,且在部分LST层上取得了比2-D算法更高的填补精度（R大于0.95,RMSE小于1.5 K,MAE小于1 K）。综合来看,2-D算法的计算优势主要表现为快速稳定的计算特性,但无法取得更高的填补精度,3-D算法虽然需要更多的计算时间且在个别LST上存在较大填补误差,但由于结合了辅助LST层提供的时间信息,在部分LST层中可取得显著优于2-D算法的填补结果。“优化算法”很好地结合了2-D与3-D算法的优势,R、RMSE、MAE相对于2-D与3-D算法均得到了提高,且计算速度达到了与2-D算法相当的水平,在整个MODIS LST数据集上表现出稳定、快速和高精度的填补性能。
图9 2-D、3-D、“优化算法”填补粤港澳大湾区2002年1月—2016年12月LST缺值的精度评价结果

Fig. 9 Accuracy of 2-D, 3-D and optimized DCT-PLS in filling missing values of LST in GBA from 2002.1 to 2016.12

根据2.2节中的计算原理,DCT-PLS算法在填补缺值时,算法基于待填补数据集自身信息,将带有缺值的原始数据在时域与频域中进行DCT与IDCT变换,并通过PLS方法计算合适的平滑程度,到达对原始数据主要信息的恢复,从而完成缺值填补。对2-D算法而言,其利用待填补LST中剩余的有效LST所提供的空间信息对缺值区域进行插值填补,这类似于一般的地统计学插值方法,可以完成对二维缺值图像的快速稳定地插值填充,但其填补精度在缺乏足够的空间信息时会受到限制。对于3-D算法而言,其在待填补LST自身空间信息基础上添加了来自辅助LST层的时间信息,并在3个维度上进行DCT与IDCT变换,三维运算虽然降低了模型的填补效率,但额外的时间信息起到了限制平滑范围的作用,使得填补结果接近辅助LST层提供的时间信息,因此,3-D算法的填补结果不具有稳定性：当时间信息与原始LST之间接近时,3-D算法可以在部分图层中取得比2-D算法更加精确的填补结果,但是当时间信息与原始LST之间偏差较大时,3-D算法将会产生有偏估计并带来较大的填补误差。而“优化算法”通过对每一层待填补LST寻找与其接近的辅助LST层,达到对辅助时间信息控制的目的,在“有效”时间信息的辅助下,“优化算法”可以完成对每一层LST的精确填补;此外,3-D算法在进行时间维的DCT与IDCT变换时,仅有少量辅助LST值参与运算,因此,本文将这些辅助LST值提供的时间信息在填补前转化为2-D算法可用的空间信息,使得填补算法仅在二维上即可同时利用目标LST自身空间信息与辅助LST有效的时间信息,达到了对填补速度与填补精度的优化。
通过“随机选取”策略对每一个目标LST随机生成人工模拟缺值,使用“优化算法”填补大湾区2001年1月—2017年12月月值MODIS LST日间产品与夜间产品的结果如图10所示。散点图显示在整个数据集上“优化算法”有效地填补了MODIS LST日间与夜间产品在大湾区的缺值,散点均匀的分布在1:1直线两侧,没有产生明显的离群值。在大湾区内,“优化算法”填补日间LST产品与夜间LST产品的R值大于0.98,RMSE小于0.79 K,MAE小于0.57 K,取得了全局高精度的填补结果。
图10 “优化算法”在粤港澳大湾区2001年1月—2017年12月LST的缺值填补结果

Fig. 10 Filling results of optimized DCT-PLS in GBA from 2001.1 to 2017.12

结合大湾区2001—2017年MODIS土地利用分类数据,日间LST与夜间LST填补结果在不同土地利用类型上没有表现出显著的精度差异（表2）。城市区域缺值比例更大,其RMSE与MAE值略高于其他土地利用类型。总之,“优化算法”在不同土地利用类型上均取得了高精度的填补结果,下垫面性质的变化没有对填补结果的精度带来明显影响。此外,本文将“优化算法”应用于填补2001—2017年MODIS LST日值数据（MOD11A1）中的缺值,取得了与月值数据集相似的高精度填补结果。
表2 粤港澳大湾区2001年1月—2017年12月日间和夜间LST缺值填补结果在不同土地利用类型中的精度

Tab. 2 Accuracy evaluation of day-time and night-time LST filling results for different land covers in GBA from 2001.1 to 2017.12

土地覆盖类型 RMSE/K R MAE/K
日间 夜间 日间 夜间 日间 夜间
林地 0.75 0.55 0.98 0.99 0.57 0.39
灌木林 0.78 0.58 0.97 0.99 0.58 0.41
草地 0.79 0.57 0.98 0.99 0.58 0.42
耕地 0.79 0.60 0.98 0.99 0.56 0.43
城市区域 0.85 0.67 0.98 0.99 0.59 0.45
基于以上分析,“优化算法”在MODIS LST数据集内有效辅助LST层的支持下,可完成对长时间序列数据缺值的有效填补。利用填补得到的2001年1月—2017年12月大湾区时空连续的LST数据集,本文计算了大湾区2001—2017年日间与夜间的LST均值空间分布（图11）。大湾区2001—2017年平均日间LST为292~304 K（18.85~30.85 ℃）,在广州、佛山、中山、东莞、深圳等城市群地区日间LST均值较高,在周边地区日间LST均值较低。大湾区2001—2017年平均夜间LST为284~294 K（10.85~20.85 ℃）,在广州、佛山、中山、东莞等城市群地区夜间LST均值较高,在周边地区LST均值较低。城市群地区存在大面积的不透水层,地表在日间吸收更多热量、在夜间释放热量,使得LST在日间与夜间均高于周边地区。在广州、深圳、东莞、中山、佛山等地区,日间与夜间的LST与周边地区差值在10 K左右,且日间LST具有更大面积的高温区域。
图11 粤港澳大湾区2001年1月—2017年12月平均LST空间分布

注：该图基于自然资源部标准地图服务网站下载的审图号为GS(2019)4342号的标准地图制作,底图无修改。

Fig. 11 Spatial distribution of average LST in GBA from 2001.1 to 2017.12

6 结论与展望

6.1 结论

针对卫星观测数据产品中普遍存在的数据缺失问题,本文介绍了一种快速填补长时间序列三维时空数据缺值的DCT-PLS优化算法,并在粤港澳大湾区进行了MODIS LST月值产品缺值填补的实证研究。2001年1月—2017年12月大湾区日间LST产品与夜间LST产品的缺值填补实验表明,DCT-PLS优化算法适用于填补长时间序列MODIS LST中的缺值。该优化算法结合了2-D算法与3-D算法的填补优势,可以达到对填补速度与填补精度的同步提升（在整个时间序列上,日间LST产品R均值达到0.97,分别提高0.06（2-D）与0.03（3-D）;RMSE均值达到0.74 K,分别降低1.14 K（2-D）与1.20 K （3-D）;MAE均值达到0.55 K,分别降低0.76 K（2-D）与1.02 K（3-D）;平均计算时间降低10.3 s（3-D））。
与其他MODIS LST缺值填补方法相比,DCT-PLS算法具有以下特点：① 多维填补：基于DCT与IDCT变换的多维计算特性使其可以很方便的扩展到高维数据中,从而应用于三维栅格时空数据集的缺值填补;② 独立于外部数据集：算法基于待填补数据集自身信息,将带有缺值的原始数据在时域与频域中进行DCT与IDCT变换,并通过PLS方法控制模型平滑程度保留数据集中的高频分量信息,完成对数据主要信息的恢复,降低了对高精度外源数据集的依赖,避免了外源数据集引入的额外误差与不确定性,且使算法模型更易于迁移应用;③ 稳定、快速的计算特性：算法对每一层目标LST,算法均通过模型内唯一的平滑参数控制填补结果,并使用GCV得分对平滑参数进行自动迭代计算,从而得到合适的平滑程度,达到对目标LST缺值快速稳定地填补输出,且能够适应长时间序列数据的缺值填补。

6.2 展望

基于本文研究,未来可在如下方面开展工作：① 目前所提出的“优化算法”已经在晴空条件下在大湾区内取得了全局高精度的填补结果,但针对云下LST的重建与还原能力有待于进一步探索。因此,可以对DCT-PLS算法加入适当的外源辅助数据集,探索在其他数据的支持下,DCT-PLS算法重建云下LST的能力与适应性。② 本文基于MODIS LST数据对DCT-PLS算法填补缺值的能力进行了研究,算法填补其他大型地球物理数据集（如气溶胶、地表辐射）中缺值的适应性与误差来源有待于进一步研究。③ 本文填补算法是否适用于其他复杂地形地区（如云贵高原）或高海拔地区的LST缺值填补有待于进一步探索。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
闫李月, 李洪忠, 韩宇, 等. MODIS与Landsat 8地表温度融合拼接——以粤港澳大湾区为例[J]. 热带地理, 2019, 39(5):689-700.

[ Yan L Y, Li H Z, Han Y, et al. Surface temperature splicing study fusing MODIS and landsat 8: A case study in the Guangdong-Hong Kong-Macao greater bay area[J]. Tropical Geography, 2019, 39(5):689-700. ] DOI: 10.13284/j.cnki.rddl.003186

DOI

[2]
战川, 唐伯惠, 李召良. 近地表大气逆温条件下的地表温度遥感反演与验证[J]. 遥感学报, 2018, 22(1):28-37.

[ Zhan C, Tang B H, Li Z L. Retrieval and validation of land surface temperature for atmospheres with air temperature inversion[J]. Journal of Remote Sensing, 2018, 22(1):28-37. ]

[3]
李召良, 段四波, 唐伯惠, 等. 热红外地表温度遥感反演方法研究进展[J]. 遥感学报, 2016, 20(5):899-920.

[ Li Z L, Duan S B, Tang B H, et al. Review of methods for land surface temperature derived from thermal infrared remotely sensed data[J]. Journal of Remote Sensing, 2016, 20(5):899-920. ]

[4]
管延龙, 王让会, 李成, 等. 基于MODIS数据的天山区域地表温度时空特征[J]. 应用生态学报, 2015, 26(3):681-688.

[ Guan Y L, Wang R H, Li C, et al. Spatial-temporal characteristics of land surface temperature in Tianshan Mountains area based on MODIS data[J]. Chinese Journal of Applied Ecology, 2015, 26(3):681-688. ] DOI: 10.13287/j.1001-9332.20141223.027

DOI

[5]
刘志武, 党安荣, 雷志栋, 等. 利用ASTER遥感数据反演陆面温度的算法及应用研究[J]. 地理科学进展, 2003, 22(5):507-514,544.

[ Liu Z W, Dang A R, Lei Z D, et al. A retrieval model of land surface temperature with ASTER data and its application study[J]. Progress in Geography, 2003, 22(5):507-514,544. ]

[6]
魏然, 单杰. 城市地表温度影像时空融合方法研究[J]. 武汉大学学报·信息科学版, 2018, 43(3):428-435.

[ Wei R, Shan J E. Spatial and temporal fusion for urban land surface temperature image mapping[J]. Geomatics and Information Science of Wuhan University, 2018, 43(3):428-435. ] DOI: 10.13203/j.whugis2015048

DOI

[7]
Xu Y M, Shen Y. Reconstruction of the land surface temperature time series using harmonic analysis[J]. Computers & Geosciences, 2013, 61:126-132. DOI: 10.1016/j.cageo.2013.08.009

DOI

[8]
于晓静, 刘焕军, 曲长祥, 等. 基于MODIS数据的黑龙江省黑土带地温空间格局分析[J]. 水土保持研究, 2012, 19(2):66-70,287.

[ Yu X J, Liu H J, Qu C X, et al. Study on the spatial pattern of land surface temperature in black-soil zone, Heilongjiang Province[J]. Research of Soil and Water Conservation, 2012, 19(2):66-70,287. ]

[9]
Julien Y, Sobrino J A. Comparison of cloud-reconstruction methods for time series of composite NDVI data[J]. Remote Sensing of Environment, 2010, 114(3):618-625. DOI: 10.1016/j.rse.2009.11.001

DOI

[10]
Zhou J E, Jia L, Menenti M. Reconstruction of global modis ndvi time series: Performance of harmonic analysis of time series (hants)[J]. Remote Sensing of Environment, 2015, 163:217-228. DOI: 10.1016/j.rse.2015.03.018

DOI

[11]
周芳成, 宋小宁, 李召良. 地表温度的被动微波遥感反演研究进展[J]. 国土资源遥感, 2014, 26(1):1-7.

[ Zhou F C, Song X N, Li Z L. Progress of land surface temperature retrieval based on passive microwave remote sensing[J]. Remote Sensing for Land & Resources, 2014, 26(1):1-7. ]

[12]
Guillemot C, Le Meur O. Image inpainting: Overview and recent advances[J]. IEEE Signal Processing Magazine, 2014, 31(1):127-144. DOI: 10.1109/MSP.2013.2273004

DOI

[13]
柯灵红, 王正兴, 宋春桥, 等. 青藏高原东北部MODIS LST时间序列重建及与台站地温比较[J]. 地理科学进展, 2011, 30(7):819-826.

[ Ke L H, Wang Z X, Song C Q, et al. Reconstruction of MODIS LST time series and comparison with land surface temperature (T) among observation stations in the northeast Qinghai-Tibet plateau[J]. Progress in Geography, 2011, 30(7):819-826. ]

[14]
Neteler M. Estimating daily land surface temperatures in mountainous environments by reconstructed MODIS LST data[J]. Remote Sensing, 2010, 2(1):333-351. DOI: 10.3390/rs1020333

DOI

[15]
柯灵红, 王正兴, 宋春桥, 等. 青藏高原东北部MODIS地表温度重建及其与气温对比分析[J]. 高原气象, 2011, 30(2):277-287.

[ Ke L H, Wang Z X, Song C Q, et al. Reconstruction of MODIS land surface temperature in northeast Qinghai-Xizang plateau and its comparison with air temperature[J]. Plateau Meteorology, 2011, 30(2):277-287. ]

[16]
Shen H F, Li X H, Cheng Q, et al. Missing information reconstruction of remote sensing data: A technical review[J]. IEEE Geoscience and Remote Sensing Magazine, 2015, 3(3):61-85. DOI: 10.1109/MGRS.2015.2441912

DOI

[17]
韩晓勇, 韩玲, 戚鹏程. 秦巴山区MODISLST时序数据重建及特征分析[J]. 地理与地理信息科学, 2016, 32(1):71-77.

[ Han X Y, Han L, Qi P C. Reconstruction of MODIS LST time series and features analysis in Qinling-Bashan mountain area[J]. Geography and Geo-information Science, 2016, 32(1):71-77. ]

[18]
Kang J, Tan J L, Jin R, et al. Reconstruction of MODIS land surface temperature products based on multi-temporal information[J]. Remote Sensing, 2018, 10(7):1112. DOI: 10.3390/rs10071112

DOI

[19]
Zeng C, Shen H F, Zhong M L, et al. Reconstructing MODIS LST based on multitemporal classification and robust regression[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(3):512-516. DOI: 10.1109/LGRS.2014.2348651

DOI

[20]
Liu X, Deng C W, Chanussot J, et al. StfNet: A two-stream convolutional neural network for spatiotemporal image fusion[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(9):6552-6564. DOI: 10.1109/TGRS.2019.2907310

DOI

[21]
臧琳, 宋冬梅, 单新建, 等. 基于被动微波与时空联合算法的云下像元LST重建[J]. 遥感技术与应用, 2016, 31(4):764-772.

[ Zang L, Song D M, Shan X J, et al. Reconstruction of LST under the cloud based on passive microwave remote sensing and spatio-temporal domain algorithm[J]. Remote Sensing Technology and Application, 2016, 31(4):764-772. ]

[22]
Sun L A, Chen Z X, Gao F, et al. Reconstructing daily clear-sky land surface temperature for cloudy regions from MODIS data[J]. Computers & Geosciences, 2017, 105:10-20. DOI: 10.1016/j.cageo.2017.04.007

DOI

[23]
Garcia D. Robust smoothing of gridded data in one and higher dimensions with missing values[J]. Computational Statistics & Data Analysis, 2010, 54(4):1167-1178. DOI: 10.1016/j.csda.2009.09.020

DOI

[24]
Garcia D. A fast all-in-one method for automated post-processing of PIV data[J]. Experiments in Fluids, 2011, 50(5):1247-1259. DOI: 10.1007/s00348-010-0985-y

DOI

[25]
Pham H T, Kim S, Marshall L, et al. Using 3D robust smoothing to fill land surface temperature gaps at the continental scale[J]. International Journal of Applied Earth Observation and Geoinformation, 2019, 82:101879. DOI: 10.1016/j.jag.2019.05.012

DOI

[26]
Liu H Z, Lu N, Jiang H, et al. Filling gaps of monthly terra/MODIS daytime land surface temperature using discrete cosine transform method[J]. Remote Sensing, 2020, 12(3):361. DOI: 10.3390/rs12030361

DOI

[27]
Kim S, Paik K, Johnson F M, et al. Building a flood-warning framework for ungauged locations using low resolution, open-access remotely sensed surface soil moisture, precipitation, soil, and topographic information[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2018, 11(2):375-387. DOI: 10.1109/JSTARS.2018.2790409

DOI

[28]
Wang G J, Garcia D, Liu Y, et al. A three-dimensional gap filling method for large geophysical datasets: Application to global satellite soil moisture observations[J]. Environmental Modelling & Software, 2012, 30:139-142. DOI: 10.1016/j.envsoft.2011.10.015

DOI

[29]
张桂欣, 郝振纯, 祝善友, 等. AMSR2缺失数据重建及其土壤湿度反演精度评价[J]. 农业工程学报, 2016, 32(20):137-143.

[ Zhang G X, Hao Z C, Zhu S Y, et al. Missing data reconstruction and evaluation of retrieval precision for AMSR2 soil moisture[J]. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(20):137-143. ]

[30]
Lei F N, Crow W T, Shen H F, et al. Assessment of the impact of spatial heterogeneity on microwave satellite soil moisture periodic error[J]. Remote Sensing of Environment, 2018, 205:85-99. DOI: 10.1016/j.rse.2017.11.002

DOI

[31]
张慧娴. 中国农田植被时空变化特征及其气候影响因子分析[D]. 南京:南京信息工程大学, 2018.

[ Zhang H X. Spatial and temporal characteristics of cropland vegetation variations in China and related climate factor analysis[D]. Nanjing: Nanjing University of Information Science & Technology, 2018. ]

[32]
Fredj E, Kohut J, Roarty H, et al. Fast gap filling of the coastal ocean surface current in the seas around Taiwan[C]// OCEANS 2016 - Shanghai. IEEE, 2016:1-4.

[33]
肖江洪. 高频地波雷达海流数据处理方法研究[D]. 武汉:武汉大学, 2017.

[ Xiao J H. Research on current data processing method of high frequency surface wave radar[D]. Wuhan: Wuhan University, 2017. ]

[34]
彭晓琳. 全球卫星臭氧产品的时空重建研究[D]. 武汉:武汉大学, 2017.

[ Peng X L. Spatio-temporal reconstruction for globally remotely sensed total ozone production[D]. Wuhan: Wuhan University, 2017. ]

Options
文章导航

/

〈 〉