基于Savitzky-Golay滤波算法的FY-2F地表温度产品时间序列重建

图1 研究区气象站点分布示意图

Fig.1 Distribution of meteorological stations in the study area

1.3 基于S-G滤波的LST重建方法

受研究区内云量、传感器等因素影响,LST遥感数据存在数值异常或缺值现象,FY-2F LST逐时数据即存在大量缺值,即使是每日合成产品仍存在很多无效像元,严重影响了数据的可用性,因此能否精确地重建LST数据就尤为重要。本文基于LST长时序变化特征,并结合滤波算法重建LST缺值。S-G滤波由Savitzky和Golay^[15]在1964年首次提出,是一种在时域内基于局部多项式最小二乘拟合的滤波方法。通过取点 $x_{i}$ 相邻时相(前后)固定个数点拟合多项式,多项式在 $x_{i}$ 的值即为S-G滤波的模拟值。该滤波器最大的特点为在滤除噪声的同时可以保持曲线的变化趋势。在对FY-2F LST数据重建的过程中,S-G滤波方法能够最大程度地保留原始数据,并且根据相邻数据模拟出缺值。S-G滤波过程可表示为

(1)

Y_{j}^{*} = \frac{\overset{m}{\sum_{i = - m}} C_{i} Y_{j + i}}{N}

式中: $Y_{j}^{*}$ 为重建后LST数据集; $Y_{j + i}$ 为LST原始值; $C_{i}$ 为滤波系数,即第i期原始LST数据 $Y_{j}$ 在平滑窗口中的权重; N为滑动窗口的数据个数(2m+1),其中m为窗口宽度。

当像元长时间连续缺失较多数据时,若单纯使用S-G滤波进行数据重建则可能与真实值相差较大。针对这一问题本文筛选出某像元一段不存在缺值的数据作为研究对象,从该段数据中分别提取15个相邻时相数据不参加S-G滤波,模拟出像元不同连续缺值情况,并将重建前、后LST值进行对比。图2中绘制出原始LST和5条不同连续缺值重建后曲线,方框代表模拟缺值,从图中可以清楚地看到连续缺失数据越多,模拟结果与原始LST差值越大,当连续缺失达到5个数据时误差大于3 K。

图2

图2 不同连续缺值重建效果

Fig.2 Reconstruction effects of different missing values

因此,本文首先遍历每一像元所有时相LST数据,并根据像元连续缺值情况将全年数据拆分为若干组,当连续缺失5个及以上时相数据时该段LST不进行重建。随后,根据实测数据统计出LST最大、最小值作为阈值,判断LST异常点,利用接口描述语言(interface description language,IDL)函数库提供的S-G滤波函数重建异常点及缺失数据,并且保留原始有效数据不变,通过反复迭代达到最佳拟合效果。S-G滤波函数需设定窗口宽度m及多项式拟合阶数Degree,通常情况下m值越大滤波结果越平滑; Degree一般设定在24范围内,阶数越低滤波结果平滑,但会带来误差,较高阶数可降低此项误差但“过拟合”会带来噪声影响滤波结果^[16-17]。为了确定S-G滤波LST重建最优参数,本文在连续30 d数据中随机选择5 d数据并赋为缺值,再按照窗口宽度m=3,5,7及阶数Degree=2,3,4组成9组参数组合进行S-G滤波,经实验平滑窗口宽度m取5,阶数Degree为4时模拟LST与原始值最为接近。

2 结果与分析

2.1 LST重建结果

本研究对2013年FY-2F LST时间序列数据进行了重建,通过IDL函数库提供的S-G滤波函数迭代拟合,最终达到最优拟合效果。S-G滤波能够很好地将FY-2F LST数据缺失值补充完整,图3展示了长三角地区每一像元点重建前、后时间序列缺值率(即表示每个像元点长时间序列上的缺值率,下文简称时相缺失率)。经统计,长三角地区重建前FY-2F LST平均时相缺失率为19.43%,其中单一像元最少缺失9.54%数据,最高时相缺失率则达到39.31%,整个长三角地区南部数据时相缺失率高于北部。经过重建后,FY-2F LST时相缺失率明显降低,长三角大部分区域LST补充完整,平均时相缺失率下降至1.69%,可见S-G滤波方法具有良好的重建效果。

图3

图3 LST时相缺失率空间分布

Fig.3 Distribution of LST vacancy rate in time-series

为了更直观地展现LST重建效果,图4选取不同空间缺值率的FY-2F LST原始数据与重建后数据进行对比。从图中可以看到,S-G滤波利用LST时序变化特征进行数据重建,重建效率较高,即便研究区内空间缺值率超过80%也能很好地完成重建工作,且能够较好地保证LST空间一致性。

图4

图4 FY-2F LST重建示例

Fig.4 Examples of FY-2F LST reconstruction result

为了从LST时间序列角度评价重建效果,本文进一步选取了南京站点作为验证点,对比站点实测值与重建前、后LST随时间变化曲线,如图5所示。

图5

图5 南京站点LST重建效果

Fig.5 LST reconstruction effect of Nanjing meteorological station

由南京站点2013年LST随时间变化总体曲线(图5(a))可以看出,本文重建方法可以有效地保留原始LST,并达到对缺值进行补充的目的。实测数据与重建前、后LST相比变化趋势一致,全年整体随时间变化先升高后降低。从图中也能看到,重建后的LST时间序列与实测数据仍然存在一定的误差,其中8月份的差值达到了10 K左右。为了更清楚地观察重建效果,分别提取4月及8月份数据如图5(b)和(c)。S-G滤波方法对LST时间序列进行滤波获取其长时间变化趋势,结合相邻时相点对缺值进行补充,使其符合LST变化规律。图5(b)中,实测数据与重建后FY-2F LST变化趋势较为一致,相邻2个时相FY-2F数据波动更大; 8月份(图5(c))存在连续缺失数据较多的情况,为了保证重建数据准确性该段缺值未进行补充,其余部分缺值经S-G滤波拟合效果较好,但与实测LST差值偏大,这种情况受数据本身精度的影响较大。

2.2 精度验证

为了进一步定量检验重建后LST结果的精度,根据经纬度提取了长三角地区46个气象站点重建前、后LST每日数据,并与相应日期实测LST数据进行对照。图6(a)和(b)分别展示了气象站点实测数据与LST重建前、后对比情况。从图中可以看到FY-2F LST原始产品与实测数据具有较好的拟合精度,判定系数R²为0.75,而平均绝对误差(mean absolute error,MAE)为4.88 K,原始LST产品精度较低。重建后样本数由11 201个增加至13 621个,R²及MAE与重建前基本保持一致,因此利用S-G滤波对FY-2F LST进行重建能够保证数据精度,重建结果较为可信。

图6

图6 实测值与重建前、后LST散点图

Fig.6 Scatter plots of the measured value, before and after reconstruction

经过对比发现,气象站点实测值与重建后LST值误差较大,主要有以下3点原因: ①FY-2F原始LST产品在某些地区误差较大,甚至超过10 K,由图5可以看出,高误差像元主要出现在夏季7—8月份,这可能受FY-2F LST产品本身反演算法的影响,也可能是该地区夏季云量较多,导致逐时数据缺失,进而影响了日平均LST数据的计算; ②气象站点实测值为“点数据”,而FY-2F LST数据为“面数据”,且空间分辨率为5 km,导致像元多为混合像元,由于两者尺度不同,在匹配时必然带来一定误差^[18,19]; ③本文是在FY-2F LST产品的基础上进行重建的,受原始产品质量影响,无法解决一些时相的高误差问题。因此,重建后产品与重建前产品相比更加合理。

综合以上分析,为了更加明确S-G滤波方法重建效果,本文随后采取人工模拟缺值的方法排除由于数据问题导致的影响,对比重建前、后数据验证S-G滤波重建数据的准确性(图7)。

图7

图7 重建前、后LST散点图

Fig.7 Scatter plots of LST before and after reconstruction

通过随机选取一些非缺值区域不参与S-G滤波,然后跟滤波之后的LST结果进行对比。选取验证点的时候,要求在长三角地区内均匀分布,且影像按照不同季节随机抽取,最大程度保证了验证结果的准确性。从图7中可见,重建前、后LST的MAE为1.35 K,R²为0.95,大部分样本都位于1∶1线附近,重建后LST与原始值吻合度较好。因此,除去数据质量等其他因素影响,利用S-G滤波方法进行数据重建效果良好。

为了进一步评价本文重建方法的鲁棒性,本文使用相同方法对2014年和2015年FY-2F LST日均值数据进行了重建,评价其R²和MAE,结果如表2所示。2014—2015年间站点实测值与重建前、后LST的R²和MAE基本持平,能够保证重建后LST数据精度不变; 人工模拟缺值点验证结果良好,与2013年结果大体相同,模型鲁棒性较好。

表2 2014—2015年LST重建精度评价

Tab.2 Evaluation of LST reconstruction accuracy from 2014 to 2015

类别	2014年		2015年
类别	R²	MAE/K	R²	MAE/K
实测值与重建前LST	0.65	4.52	0.69	4.31
实测值与重建后LST	0.63	4.53	0.68	4.32
人工模拟缺值点	0.95	1.37	0.96	1.58

新窗口打开| 下载CSV

3 结论

本文利用S-G滤波方法,结合LST时间序列变化特征,对长三角地区FY-2F LST产品进行了数据重建,并对重建结果进行了评价。结果表明:

1)长三角地区FY-2F LST原始产品平均时相缺失率为19.43%,S-G滤波方法能够有效地重建LST,重建后平均时相缺失率降低为1.69%。并且对于不同空间缺值率影像,该方法均能在保持LST空间一致性的基础上对其进行有效重建。

2)提取长三角地区46个气象站重建前、后LST值,结合站点LST进行精度验证与误差分析。重建后结果精度与重建前FY-2F LST产品基本持平,但与实测值间误差较大,其原因为: 产品本身质量影响了重建结果精度; 进行验证时,实测值为“点数据”,与遥感“面数据”因尺度问题会进一步扩大误差。

3)通过人为模拟验证区缺值,将重建后结果与重建前对比,发现重建精度较高,两者R²为0.95,MAE为1.35 K。通过重建2014—2015年间FY-2F LST日均值产品验证S-G滤波模型适用性,其结果与2013年大体相同,表明S-G滤波法可以用来进行LST重建。

本文在保持FY-2F LST原始产品精度的基础上进行数据重建,能够有效提高产品质量,为地表热环境研究提供数据基础。然而文中重建采用的S-G滤波方法仅考虑了LST时序变化规律,并未利用到其空间信息。在今后的研究中可结合多源遥感数据和下垫面特性等因子进行空间上的重建,进一步提高产品数据质量。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Dash

, Göttsche F

, Olesen F

, et al.

Land surface temperature and emissivity estimation from passive sensor data:Theory and practice-current trends

[J]. International Journal of Remote Sensing, 2002,23(13):2563-2594.

DOI:10.1080/01431160110115041 URL [本文引用: 1]

[2]

周芳成, 宋小宁, 李召良 .

地表温度的被动微波遥感反演研究进展

[J]. 国土资源遥感, 2014,26(1):1-7.doi: 10.6046/gtzyyg.2014.01.01.

<p>微波遥感能穿透云层，甚至可穿透一定程度的雨区，可以弥补热红外遥感的不足。发展基于被动微波遥感的地表温度反演算法可以全天候地为相关领域提供数据服务。根据前人研究，该文从方法论的角度将已有的反演方法分为统计模型法、物理模型法和神经网络算法3类，分析了每种方法的优缺点，并探讨了未来微波遥感反演地表温度的发展方向，以期为进一步研究提供参考。</p>

Zhou F

, Song X

, Li Z

Progress of land surface temperature retrieval based on passive microwave remote sensing

[J]. Remote Sensing for Land and Resources, 2014,26(1):1-7.doi: 10.6046/gtzyyg.2014.01.01.

[3]

李天祺, 朱秀芳, 潘耀忠 , 等.

MODIS陆地表面温度数据重构方法研究

[J]. 北京师范大学学报(自然科学版), 2015,51(s1):70-76.

Li T

, Zhu X

, Pan Y

, et al.

Study on reconstruction methods of MODIS LST products

[J]. Journal of Beijing Normal University (Natural Science), 2015,51(s1):70-76.

DOI:10.11867/j.issn.1001-8166.2011.07.0724 URL Magsci [本文引用: 1]

[4]

祝善友, 张桂欣 .

近地表气温遥感反演研究进展

[J]. 地球科学进展, 2011,26(7):724-730.

<p>高时间分辨率的近地表气温空间分布数据是许多陆面过程模型中非常重要的输入参数之一。在常规气象观测站点稀少或没有的情况下，利用遥感技术进行较高时空分辨率的近地表气温估算与反演，在理论方法与业务实践上都具有重要研究意义。根据地表能量平衡与辐射平衡原理，在气温遥感反演物理机制分析的基础上，总结了国内外近年来气温遥感反演的研究进展，主要方法可归纳为5类：单因子统计方法、多因子统计方法、神经网络方法、地表温度—植被指数方法和地表能量平衡方法，并从遥感反演气温的时空分辨率、反演模型中影响因子的考虑、模型的可移植性与实用性角度，讨论了已有研究方法中存在的困难与问题，最后对未来可能的研究方向做出了展望。</p>

Zhu S

, Zhang G

Progress in near surface air temperature retrieved by remote sensing technology

[J]. Advances in Earth Science, 2011,26(7):724-730.

[5]

Wan

, Zhang

, et al.

Validation of the land-surface temperature products retrieved from Terra moderate resolution imaging spectroradiometer data

[J]. Remote Sensing of Environment, 2002,83(1-2):163-180.

DOI:10.1016/S0034-4257(02)00093-7 URL [本文引用: 1]

[6]

涂丽丽, 覃志豪, 张军 , 等.

基于空间内插的云下地表温度估计及精度分析

[J].遥感信息, 2011(4):59-63.

Tu L

, Qin Z

, Zhang

, et al.

Estimation and error analysis of land surface temperature under the cloud based on spatial interpolation

[J].Remote Sensing Information, 2011(4):59-63.

[7]

李儒, 张霞, 刘波 , 等.

遥感时间序列数据滤波重建算法发展综述

[J]. 遥感学报, 2009,13(2):335-341.

DOI:10.3321/j.issn:1007-4619.2009.02.023 URL Magsci [本文引用: 1]

遥感时间序列数据(MODIS,NOAA/AVHRR,SPOT/VEGETATION等)在植被生长监测、物候信息提取、土地利用类型监测等诸多领域得到了广泛应用,是生产研究的重要数据源之一.由于传感器、云层大气等影响,遥感时间序列数据存在着严重的噪声,应用前必须进行序列滤波重建工作.综述现有各类滤波重建方法,对研究中广为采用的3类主要方法(基于最小二乘的非对称高斯函数拟合、SavitZky-Golay滤波、基于离散傅里叶的系列分析方法)集中阐述其理论基础、应用步骤和优缺点.总结当前遥感时间序列滤波重建方法需要进一步改进之处.

, Zhang

, Liu

, et al.

Review on methods of remote sensing time-series data reconstruction

[J]. Journal of Remote Sensing, 2009,13(2):335-341.

[8]

边金虎, 李爱农, 宋孟强 , 等.

MODIS植被指数时间序列Savitzky-Golay滤波算法重构

[J]. 遥感学报, 2010,14(4):725-741.

利用Savitzky-Golay(S-G)滤波方法对若尔盖高原湿地区2000-2009年MODIS 16d最大值合成的NDVI时间序列数据进行了重构,并与中值迭代滤波法、傅里叶变换法进行了比较.结果表明,基于S-G滤波的时间序列重构方法重构后的NDVI时间序列在直观及像元的时间序列曲线上均取得了较好的效果,对提高该数据产品质量有很大帮助,通过该方法重构后的高质量的NDVI时间序列对利用该数据源对若尔盖湿地生态系统监测提供了良好的基础.

Bian J

, Li A

, Song M

, et al.

Reconstruction of NDVI time-series datasets of MODIS based on Savitzky-Golay filter

[J]. Journal of Remote Sensing, 2010,14(4):725-741.

[9]

Nguyen O

, Kawamura

, Trong D

, et al.

Temporal change and its spatial variety on land surface temperature and land use changes in the Red River Delta,Vietnam,using MODIS time-series imagery

[J]. Environmental Monitoring and Assessment, 2015,187(7):1-11.

DOI:10.1007/s10661-014-4167-x URL [本文引用: 1]

[10]

Neteler

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 URL [本文引用: 1]

[11]

臧琳, 宋冬梅, 甘宇亮 , 等.

统计模型与滤波算法的地表温度重建方法探讨

[J]. 测绘科学, 2016,41(7):11-17.

Zang

, Song D

, Gan Y

, et al.

Study on reconstruction of LST based on the statistical model and filtering algorithm

[J]. Science of Surveying and Mapping, 2016,41(7):11-17.

[12]

韩冬锐, 徐新良, 李静 , 等.

长江三角洲城市群热环境安全格局及土地利用变化影响研究

[J]. 地球信息科学学报, 2017,19(1):39-49.

Han D

, Xu X

, Li

, et al.

Study on the security pattern of the heat environment and the influence of land use change in the Yangtze River Delta urban agglomeration

[J]. Journal of Geo-Information Science, 2017,19(1):39-49.

[13]

Ooka

Recent development of assessment tools for urban climate and heat-island investigation especially based on experiences in Japan

[J]. International Journal of Climatology, 2007,27(14):1919-1930.

DOI:10.1002/(ISSN)1097-0088 URL [本文引用: 1]

[14]

姚远, 陈曦, 钱静 .

城市地表热环境研究进展

[J]. 生态学报, 2018,38(3):1134-1147.

Yao

, Chen

, Qian

Research progress on the thermal environment of the urban surfaces

[J]. Acta Ecologica Sinica, 2018,38(3):1134-1147.

DOI:10.3969/j.issn.1000-6362.2015.01.012 Magsci [本文引用: 1]

[15]

Savitzky

, Golay

M J E

Smoothing and differentiation of data by simplified least squares procedures

[J]. Analytical Chemistry, 1964,36(8):1627-1639.

DOI:10.1021/ac60214a047 URL [本文引用: 1]

[16]

Chen

, Jonsson

, Tamura

, et al.

A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter

[J]. Remote Sensing of Environment, 2004,91(3-4):332-334.

DOI:10.1016/j.rse.2004.03.014 URL [本文引用: 1]

[17]

权文婷, 周辉, 李红梅 , 等.

基于S-G滤波的陕西关中地区冬小麦生育期遥感识别和长势监测

[J]. 中国农业气象, 2015,36(1):93-99.

采用SPOT VGT S10产品数据，通过S G滤波处理，重构NDVI时间序列，分析陕西关中地区冬小麦返青期和抽穗期在遥感影像上NDVI时序曲线变化特征，对冬小麦的返青期和抽穗期进行遥感识别。通过与站点观测的生育期数据比较，遥感识别的均方根误差大部分在±20d之内，返青期均方根误差（RMSE）为14.29d，抽穗期RMSE为9.16d，表明基于SPOT VGT数据对关中地区冬小麦返青期和抽穗期进行遥感识别是可行的。采用差值模型基于环境减灾小卫星（HJ星）数据进行渭南市2012年冬小麦抽穗期长势监测研究，结果显示大部分地区2012年冬小麦长势稍好于2011年，经与实地调查资料比较，监测结果与实际情况基本一致。采用HJ星数据进行冬小麦长势监测，结果纹理清晰，且以同源数据提取的冬小麦种植面积作为掩膜剔除非麦区域，降低了重采样带来的误差，用其进行长势监测研究具有明显优势。

Quan W

, Zhou

, Li H

, et al.

Remote recognition and growth monitoring of winter wheat in key stages based on S-G filter in Guanzhong region

[J]. Chinese Journal of Agrometeorology, 2015,36(1):93-99.

DOI:10.3969/j.issn.0001-5733.2011.03.014 Magsci [本文引用: 1]

[18]

陈顺云, 刘培洵, 刘力强 , 等.

遥感与实测地表温度的对比分析及在地震研究中的意义

[J]. 地球物理学报, 2011,54(3):747-755.

伴随着空间观测技术的发展,卫星热红外遥感在地震领域受到越来越多的关注,同时存在许多基础性工作亟待完善.本文以由卫星遥感影像与实际测量两种不同方法获取的地表温度为基础,选取2006年3月~2008年2月近2年的数据,进行遥感与实测地表温度之间的对比研究.分析结果表明遥感与实测地表温度之间:夜间差值比白天要小,白天的相关性比夜间相关性差;无论是白天还是夜间,长周期信号相关性比短周期要好,对于几十天以上的变化,两者相关程度均极高.意味着分析地表温度的长周期变化时,选用白天或者夜间数据产品没有实质性区别.通过本文分析,初步明确了遥感地表温度与实测地表温度之间的相关性,有助于选取合适的遥感资料开展地震研究.

Chen S

, Liu P

, Liu L

, et al.

Comparative analysis between land surface temperatures obtained by field measurement and satellite remote sensing and its implication in earthquake research

[J]. Chinese Journal of Geophysics, 2011,54(3):747-755.

[19]

杨贵军, 孙晨红, 历华 .

黑河流域ASTER与MODIS融合生成高分辨率地表温度的验证

[J]. 农业工程学报, 2015,31(6):193-200.

融合多源遥感数据生成高时空分辨率数据具有重要的应用价值。为了解决高空间分辨率数据重访周期长及云雨天气带来的数据短缺问题，该文基于增强自适应的遥感图像时空融合方法（enhanced spatial and temporal adaptive reflectance fusion model, ESTARFM），使用多时相MODIS数据提供地物时间变化信息，结合ASTER影像提供的空间细节信息，选择多波段数据（可见光近红外数据和地表温度数据）共同作为输入变量融合生成高时空地表温度。融合结果分别与地表红外辐射计观测温度和ASTER温度产品进行了验证。验证结果表明：基于ESTARFM方法降尺度地表温度影像清晰，融合结果与地表红外辐射计观测温度呈显著的线性正相关关系，相关系数均高于为0.71，预测得到的地表温度与真实测得的数据的平均绝对偏差均低于2.00 K，均方根误差均低于2.60 K。与ASTER地表温度产品的验证中，整体验证结果的R2均在0.95以上。此外，ESTARFM方法在各个地类中的融合效果较好，均表现出非植被区域的相关性高于植被和水体，尤其在2012年8月27日非植被的R2达到0.91。

Yang G

, Sun C

, Li

Verification of high-resolution land surface temperature by blending ASTER and MODIS data in Heihe River Basin

[J]. Transactions of the Chinese Society of Agricultural Engineering, 2015,31(6):193-200.