周笑天 张茜茹 郭庆燕 陈益玲 周雪松 李长军 冯勇 张平
摘要 气温作为研究气候演变最基础的物理量,其日值序列的完整性和准确性对于气候分析与评估工作有着重要意义。近些年随着大量无人值守地面加密自动气象站的布设,不断出现随机站点和随机长度这种双随机特点的气象资料序列缺失,给气候分析和业务应用造成了不小的障碍。针对现有气象数据插补方案的不足,提出了一种全新的基于动态时间规整(dynamic time warping,DTW)的气温日值数据二次插补方法。该方法采用了一种实时的插补策略,主要技术内容包括:1)利用一元线性回归方程将原始气温观测时间序列分解出拟合直线和残差曲线,并将二者重构组成新的气温序列;2)给出了气温插补区的定义和插补条件;3)提出了利用动态时间规整方法计算站点间距离的新模式。利用山东省2021年的气温实况数据对该方法进行了双随机检验,检验结果表明:该方法可以满足日平均气温、日最高气温和日最低气温数据的插补需求;在插补流程中采用DTW距离测度和二次插补的组合方法,其插补效果优于目前常见的基于站点地理临近关系的组合方法;该方法对地形有一定的敏感性,平原或丘陵地区的插补效果要优于山地地区。
关键词气温日值;动态时间规整;重构;二次插补
气象资料序列的完整性是开展气候分析与评估的必要条件之一。随着我国精密气象监测系统建设的不断推进、地面气象观测站网布局逐步优化,区域观测盲区得到进一步消除,观测要素短板也实现了补足,为开展中小尺度灾害性天气监测预警业务和区域小气候特征分析创造了有利条件。与此同时,随着大量无人值守的地面加密自动气象站的布设,因仪器故障、通信中断、自然灾害等原因引起的观测中断或数据质量异常的发生概率大增,不断出现随机站点和随机长度这种双随机特点的气象资料序列缺失,给气候分析和业务应用造成了不小的障碍。
气温作为研究气候演变最基础的物理量,其日值序列的完整性和准确性对于气候统计与分析有着重要意义,也是业内关注的焦点之一(Della-Marta and Wanner,2006;Hansen et al.,2010;高庆九等,2018)。目前,国内外对气温日值缺失数据的插补,主要采用的方案是为数据缺失站点选择一个或者多个地理关系临近且气候相关性较高的数据参照站,并用参照站气温观测真值进行数据替代。王海军等(2008)采用距离最短原则,以最小绝对偏差(least absolute deviation,LAD)为目标函数求解模型参数,对处于平原地区的湖北蔡甸国家气象观测站日平均气温、日最高气温和日最低气温进行了插补试验和误差分析;余予等(2012)采用标准序列法(Steurer,1985)对1971—2000年我国2 000余个国家级地面气象观测站日平均气温进行插补试验,获得了较好的插补效果。以上插补方案虽然可以解决气温日值数据的插补,但在具体实现时,需要周边参考站累年历史同期气温平均值和标准差(余予等,2012),以便计算气候相关性,这种限定条件对于建站时间较短、迁建频繁、缺少长序列历史数据的无人值守的气象站来说并不适用,且时效上也较为滞后。
近些年,我国逐步建立并完善了包括气温要素在内的多源智能网格实况融合分析产品的小时级实时业务(李超等,2017;潘旸等,2018),客观上也促成了格点到站点间气温数据回插的实现(司鹏等,2022)。但是,由于气温实况格点是以ECMWF等数值预报产品为背景场,采用多重网格变分技术并融合地面站点的观测数据而生成的(张璐等,2017;师春香等,2019;刘莹等,2021),格点数据产品的质量依然依赖地面站点分布密度水平(龙柯吉等,2019;俞剑蔚等,2019)。当站点观测数据在时间和空间上出现随机性缺失时,关联时空范围内网格产品的稳定性也会受到极大影响(孙靖等,2021),从而极易造成数据回插异常。
因此,为了更科学和快捷地解决日益增多的双随机数据缺失问题,本文提出了一种全新的基于动态时间规整(dynamic time warping,DTW)的气温日值缺失数据二次插补方法。该方法采用了一种实时的插补策略,通过DTW距离计算、残拟分离和残拟重构等主要步骤,可以较准确地实现观测站点气温日值数据的插补。
1 资料与方法
1.1 资料来源
文中用于方法介绍、检验和分析的气温日值数据来源于山东省气象大数据云平台“天擎”数据库,且全部经过数据质量控制(周笑天等,2012;王海军等,2014;闵锦忠等,2018;叶小岭等,2019;邵宇行等,2022),时间范围涵盖2021年全年,包括日平均气温、日最高气温和日最低气温3种要素。
1.2 插补方法
1.2.1 气温时间序列的残拟分离和重构
气温日值缺失数据的插补,可以归结为气温时间序列中断点数据的预测与最优化求解问题,而一元线性回归作为天气气候业务中最基本也是最简单的预测分析方法之一(魏凤英,2007),可以将其原理和方法应用到插补过程中。
设Y=[y1,y2,…,yN]是由气温观测值构成的时间序列,长度为N,则可建立一元线性回归方程,记作:
将(2)式变换形式可得
1.2.2 气温插补区
B区(插补区)及其左邻A区和右邻C区的设定,扩展了《地面气象观测规范》(中国气象局,2003)中缺测数据内插仅限于单时次的限定,从而可以有效解决“双随机”中的长度随机问题。
1.2.3 DTW距离与参照站
动态时间规整(DTW)是时间规整与距离测度计算相结合的一种非线性规整计算方法,常被用在语音识别系统中(李正欣等,2014;闫宏宸和肖熙,2021)。该方法的最大特点就是可以通过路径规划来计算两个时间序列之间的最短累计距离,从而衡量两者的相似程度(Tormene et al.,2009;周笑天等,2022)。
本文将动态时间规整算法原理应用于气温序列的相似度计算。设X=[x1,x2,…,xN]和Y=[y1,y2,…,yN]为两个不同站点的气温时间序列,长度为N,如图2所示,上下两条实线分别代表序列X和序列Y,动态时间规整算法的路径规划过程即为从左向右刻画虚线的过程,该过程以(x1,y1)为起点,以(xN,yN)为终点,按照单调、连续的匹配原则,每向右行进一步画一条虚线,虚线两端分别匹配上下序列中的一个元素。
设行进至第k步时,虚线两端匹配的元素分别为xi和yj,则该步的步长距离为xi和yj气温之差的绝对值,记为d(wk)=|xi-yj|。
设L为起点(x1,y1)至终点(xN,yN)刻画虚线的总数量,则d(wk)最短累计步长即为序列X和序列Y的DTW距离,记作:
结合气温插补区的设定,设站点p为需要进行数据插补的插补站,气温时间序列Y(p)=[A(p),B(p),C(p)],B(p)为可插补的数据缺失序列,如果存在站点q,其气温时间序列Y(q)=[A(q),B(q),C(q)]无数据缺失,且Y(p)与Y(q)的计算区间一致,则插补站p和站点q的DTW距离记作:
DDTW(p,q)=RDTW(A(p),A(q))+RDTW(C(p),C(q))。 (6)
由式(6)可知,插补站p和站点q的DTW距离为两站A、C区序列的DTW距离之和。
将有限空间范围内的站点逐一与插补站p按照公式(6)计算DTW距离并排序,取DTW距离排序最小的站,称为参照站。
从以上定义可知,参照站是在DTW距离测度下,与插补站p气温序列最相似(DTW距离最近)的站,参照站的气温序列适合作为插补数据源,执行后续的插补操作。
DTW距离随着气温序列的变化而变化,因此是一种实时的、动态的距离计算方法。与传统的基于地理临近关系(如水平距离最近或者海拔高度最近)的参照站遴选方法不同,DTW距离不受站网分布密度的影响,可以使遴选过程更加灵活和精准,更适用于解决“双随机”中站点随机的问题。
1.2.4 插补流程
综合上述的概念和方法,设站点p为需要进行数据插补的站,B(p)为插补站p中连续数据缺失序列,则对B(p)的插补过程(流程如图3所示)描述如下:
第1步,输入插补站p的B(p)序列,内容为空值,长度为N;
第2步,确定B(p)序列的左邻序列A(p)和右邻序列C(p),长度都为N,合并得到序列Y(p)=[A(p),B(p),C(p)];
第3步,从插补站p周边临近站中,按照DTW测度得到参照站q,参照站q的序列Y(q)=[A(q),B(q),C(q)];
从图3所示的插补方法总体流程可以看出,该方法采用了数据嫁接的操作方式,共分为两个插补阶段,一次插补为B(q)对B(p)的直接替换,二次插补是在一次插补的基础上,对两站的气温序列进行残拟分离和残拟重构后,得到的插补结果。
需要说明的是,插补流程中第3步是以DTW距离测度为基准选定参照站,而当以地理距离测度为基准(如水平距离最近或者海拔高度最近等)选定参照站时,则应以相应距离计算方法代替。
2 方法检验与分析
2.1 检验设计
为了验证方法对实况数据插补的有效性,同时为了使检验过程更加完备,本文对检验条件和检验内容设计如下:
1)检验包括日平均气温、日最高气温和日最低气温3种要素在内的气温日值数据。
2)根据山东省气象地理一级区划规则,分别在鲁西北(以平原地形为主)、鲁中(以山地地形为主)、鲁南(以平原丘陵地形为主)和半岛(以丘陵和山地海岸地形为主)4个地区的每个地区中随机选择10个无人值守的地面气象观测站点作为插补测试站(图4),在兼顾地形特征的同时满足站点分布的随机性要求。
3)每个插补测试站的插补区起止时间随机产生,以满足插补长度随机性要求。
4)检验采用观测真值置空的方式模拟插补区的缺失数据,并在插补结束后计算插补值与观测真值之间的误差以评估插补效果。具体的误差评估指标包括均方根误差(RMSE):
和平均绝对误差(MAE):
5)在插补流程执行过程中,采用条件组合覆盖法,记录插补阶段(一次、二次插补)和距离测度(DTW距离、水平距离和海拔高度)的各选项组合所能产生的全部插补结果。
2.2 插补实例
本文挑选了具有代表性的插补实例进行插补过程检验。选择的插补站D0122(120.674 4°E,36.145 8°N)位于山东省青岛市,海拔高度为174 m,东、南两面临海,插补站D0122及其临近无人值守气象站站网分布如图5所示。下面以D0122站日平均气温在DTW距离测度下的插补为例,验证插补流程的可行性。
插补站D0122的日平均气温缺失日期区间为2021年7月1日至9月30日,长度为92 d,该时间段即为需要进行数据插补的B区(插补区)。由B区的日期区间范围可以确定B区的左邻A区的日期区间范围为2021年3月31日至6月30日,B区的右邻C区的日期区间范围为2021年10月1日至12月31日,A区和C区内的日平均气温数据完整,区间长度同为92 d,将A、B和C区按顺序合并,并将日期逐一映射为日序后,可得插补站D0122的日平均气温分区(图6)。
在DTW距离测度下,遍历插补站D0122的临近站(图5),按照式(6)对临近站逐一计算日平均气温的DTW距离(表1)。根据表1中的DTW距离排序,将DTW距离最小的站确定为参照站(站号D0181),其距离为218.7 ℃。
在选定D0181为参照站后,结合图7,这里对后续插补过程描述如下:
1)将参照站D0181的B区(图7a2中的B区),直接嫁接到插补站D0122的B区(图7a1中的B区红色曲线段),此为一次插补;
2)分别对一次插补序列(图7a1)和参照站D0181序列(图7a2)计算一元线性回归方程,残拟分离,获得相应的拟合线和残差线(图7b1和图7b2);
3)将一次插补的拟合线(图7b1中的B区取值部分)与参照站D0181残差线(图7b2的B区取值部分)相加重构,并再次嫁接于插补站D0122的B区(图7c1中B区紫色曲线段),完成二次插补。
图7c1中B区气温序列即为图6中B区缺失的日平均气温序列的最终插补结果,插补过程至此结束。
2.3 检验结果与分析
这里将日平均气温、日最高气温和日最低气温3种要素在山东省4个地区的双随机检验结果分别汇总在表2至表4中,并按地区分组对RMSE和MAE指标最小值进行了标注。
根据指标排序情况可知:
1)对于日平均气温(表2),在鲁西北、鲁中和半岛地区,DTW距离测度下的二次插补结果在RMSE指标和MAE指标上均表现最优;在鲁南地区,DTW距离测度下的二次插补结果虽然在RMSE指标上与一次插补表现持平(RMSE=0.406 ℃),但在MAE指标上表现更优(MAE=0.312 ℃)。
2)对于日最高气温(表3),在鲁西北、鲁中和半岛地区,DTW距离测度下的二次插补结果在RMSE指标和MAE指标上均表现最优;在鲁南地区,DTW距离测度下的二次插补结果虽然在RMSE指标上与水平距离测度下的二次插补表现持平(RMSE=0.859 ℃),但在MAE指标上表现更优(MAE=0.633 ℃)。
3)对于日最低气温(表4),在鲁西北、鲁南和半岛地区,DTW距离测度下的二次插补结果在RMSE指标和MAE指标上均表现最优;在鲁中地区,DTW距离测度下的二次插补结果虽然在RMSE指标上与海拔高度测度下的二次插补表现持平(RMSE=0.900 ℃),但在MAE指标上表现更优(MAE=0.693 ℃)。
综合检验结果可以看出,日平均气温、日最高气温和日最低气温3种要素均能成功完成双随机条件下的数据插补检验;3种要素在插补流程中采用DTW距离测度和二次插补的组合方法,其插补效果均优于基于水平距离、海拔高度的插补组合方法。
分别将表2至表4中4个地区的RMSE指标最小值作为比较对象,进一步分析不同地区的插补效果。由图8可见,日平均气温、日最高气温和日最低气温3种要素的插补误差具有近似的分布特征,均在鲁中地区最高、半岛地区其次、鲁西北和鲁南地区最低。考虑到鲁中地区和半岛地区多以山地地形为主,而鲁西北和鲁南地区多以平原和丘陵地形为主,因此可以认为,复杂地形是干扰和降低本文所提方法插补效果的一个重要因素。
3 结论
本文提出了一种实时的气温日值数据插补方法,该方法利用一元线性回归方程将气温观测时间序列分解出拟合直线和残差曲线,并通过将二者再次重构实现气温序列的重组;该方法给出了气温插补区的定义和插补区构成的充分条件,在满足条件的情况下,可以实现随机序列长度的插补需求;该方法提出了采用动态时间规整衡量气温序列相似性的新模式,使得站点间的距离计算更加科学和精准,可以满足站点随机分布的插补需求。
本文利用山东省2021年的气温实况数据,对该方法进行了双随机检验,检验结果表明:1)日平均气温、日最低气温和日最高气温3种气温日值要素均能够成功实现数据插补;2)在插补流程中采用DTW距离测度和二次插补的组合方法,其插补效果优于目前常见的基于水平距离或海拔高度等地理临近关系的组合方法;3)该方法对地形有一定的敏感性,平原或丘陵地区的插补效果要优于山地地区。
本文提出的基于时间序列的数据二次插补机制,对解决日益增加的双随机特点的气象资料缺失问题,有着较为广阔的应用前景,也为历史长序列气象数据的均一化订正提供了参考。
参考文献(References)
Della-Marta P,Wanner H,2006.A method of homogenizing the extremes and mean of daily temperature measurements[J].J Climate,19(17):4179-4197.doi:10.1175/JCLI3855.1.
高庆九,余汶樯,周小艳,2018.基于再分析资料与观测资料的中国低温阈值变化特征研究[J].大气科学学报,41(3):308-317. Gao Q J,Yu W Q,Zhou X Y,2018.Research on variation characteristics of low temperature threshold in China based on reanalysis data and observation data[J].Trans Atmos Sci,41(3):308-317.doi:10.13878/j.cnki.dqkxxb.20170903001.(in Chinese).
Hansen J,Ruedy R,Sato M,et al.,2010.Global surface temperature change[J].Rev Geophys,48(4):RG4004.doi:10.1029/2010rg000345.
李超,唐千红,陈宇,等,2017.多源数据融合系统LAPS的研究进展及其在实况数据服务中的应用[J].气象科技进展,7(2):32-38. Li C,Tang Q H,Chen Y,et al.,2017.An overview of progresses in LAPS and prospective applications in real time data service[J].Adv Meteor Sci Technol,7(2):32-38.doi:10.3969/j.issn.2095-1973.2017.02.005.(in Chinese).
李正欣,张凤鸣,李克武,等,2014.一种支持DTW距离的多元时间序列索引结构[J].软件学报,25(3):560-575. Li Z X,Zhang F M,Li K W,et al.,2014.Index structure for multivariate time series under DTW distance metric[J].J Softw,25(3):560-575.doi:10.13328/j.cnki.jos.004410.(in Chinese).
刘莹,师春香,王海军,等,2021.CLDAS气温数据在中国区域的适用性评估[J].大气科学学报,44(4):540-548. Liu Y,Shi C X,Wang H J,et al.,2021.Applicability assessment of CLDAS temperature data in China[J].Trans Atmos Sci,44(4):540-548.doi:10.13878/j.cnki.dqkxxb.20200819001.(in Chinese).
龙柯吉,师春香,韩帅,等,2019.中国区域高分辨率温度实况融合格点分析产品质量评估[J].高原山地气象研究,39(3):67-74. Long K J,Shi C X,Han S,et al.,2019.Quality assessment of high resolution temperature merged grid analysis product in China[J].Plateau Mt Meteor Res,39(3):67-74.doi:10.3969/j.issn.1674-2184.2019.03.011.(in Chinese).
闵锦忠,王晨珏,贾瑞怡,2018.苏皖地面自动站资料的质量控制及结果分析[J].大气科学学报,41(5):637-646. Min J Z,Wang C J,Jia R Y,2018.Quality control and result analysis for surface AWS data in Jiangsu and Anhui Provinces[J].Trans Atmos Sci,41(5):637-646.doi:10.13878/j.cnki.dqkxxb.20160417001.(in Chinese).
潘旸,谷军霞,宇婧婧,等,2018.中国区域高分辨率多源降水观测产品的融合方法试验[J].气象学报,76(5):755-766. Pan Y,Gu J X,Yu J J,et al.,2018.Test of merging methods for multi-source observed precipitation products at high resolution over China[J].Acta Meteorol Sin,76(5):755-766.doi:10.11676/qxxb2018.034.(in Chinese).
邵宇行,秦正坤,李昕,2022.基于EOF的高时空分辨率自动站温度观测资料质量控制[J].大气科学学报,45(4):603-615. Shao Y H,Qin Z K,Li X,2022.Quality control based on EOF for surface temperature observations from high temporal-spatial resolution automatic weather stations[J].Trans Atmos Sci,45(4):603-615.doi:10.13878/j.cnki.dqkxxb.20200516001.(in Chinese).
师春香,潘旸,谷军霞,等,2019.多源气象数据融合格点实况产品研制进展[J].气象学报,77(4):774-783. Shi C X,Pan Y,Gu J X,et al.,2019.A review of multi-source meteorological data fusion products[J].Acta Meteorol Sin,77(4):774-783.doi:10.11676/qxxb2019.043.(in Chinese).
司鹏,郭军,赵煜飞,等,2022.北京1841年以来均一化最高和最低气温日值序列的构建[J].气象学报,80(1):136-152. Si P,Guo J,Zhao Y F,et al.,2022.New series of daily maximum and minimum temperature observations for Beijing,China since 1841[J].Acta Meteorol Sin,80(1):136-152.doi:10.11676/qxxb2022.008.(in Chinese).
Steurer P,1985.Creation of a serially complete data base of high quality daily maximum and minimum temperature[M].Washington D C:National Climate Center:21.
孙靖,程光光,黄小玉,2021.中国地面气象要素格点融合业务产品检验[J].高原气象,40(1):178-188. Sun J,Cheng G G,Huang X Y,2021.The verification of gridded surface meteorological elements merging product in China[J].Plateau Meteor,40(1):178-188.doi:10.7522/j.issn.1000-0534.2019.00100.(in Chinese).
Tormene P,Giorgino T,Quaglini S,et al.,2009.Matching incomplete time series with dynamic time warping:an algorithm and an application to post-stroke rehabilitation[J].Artif Intell Med,45(1):11-34.doi:10.1016/j.artmed.2008.11.007.
王海军,涂诗玉,陈正洪,2008.日气温数据缺测的插补方法试验与误差分析[J].气象,34(7):83-91. Wang H J,Tu S Y,Chen Z H,2008.Interpolating method for missing data of daily air temperature and its error analysis[J].Meteor Mon,34(7):83-91.(in Chinese).
王海军,闫荞荞,向芬,等,2014.逐时气温质量控制中界限值检查算法的设计[J].高原气象,33(6):1722-1729. Wang H J,Yan Q Q,Xiang F,et al.,2014.Algorithm design of quality control for hourly air temperature[J].Plateau Meteor,33(6):1722-1729.doi:10.7522/j.issn.1000-0534.2014.00028.(in Chinese).
魏凤英,2007.现代气候统计诊断与预测技术[M].2版.北京:气象出版社:37. Wei F Y,2007.Modern climate statistical diagnosis and prediction technology[M].2nd ed.Beijing:China Meteorological Press:37.(in Chinese).
闫宏宸,肖熙,2021.概率线性判别分析在语音命令词置信度判决中的应用[J].计算机系统应用,30(1):54-62. Yan H C,Xiao X,2021.Application of probabilistic linear discriminant analysis in voice command confidence measures[J].Comput Syst Appl,30(1):54-62.doi:10.15888/j.cnki.csa.007732.(in Chinese).
叶小岭,陈洋,杨帅,等,2019.基于EEMD-CES的单站地面气温资料质量控制方法研究[J].大气科学学报,42(3):390-398. Ye X L,Chen Y,Yang S,et al.,2019.A quality control method of surface temperature observations based on the EEMD-CES algorithm for a single station[J].Trans Atmos Sci,42(3):390-398.doi:10.13878/j.cnki.dqkxxb.20171205001.(in Chinese).
俞剑蔚,李聪,蔡凝昊,等,2019.国家级格点实况分析产品在江苏地区的适用性评估分析[J].气象,45(9):1288-1298. Yu J W,Li C,Cai N H,et al.,2019.Applicability evaluation of the national gridded real-time observation datasets in Jiangsu Province[J].Meteor Mon,45(9):1288-1298.doi:10.7519/j.issn.1000-0526.2019.09.009.(in Chinese).
余予,李俊,任芝花,等,2012.标准序列法在日平均气温缺测数据插补中的应用[J].气象,38(9):1135-1139. Yu Y,Li J,Ren Z H,et al.,2012.Application of standardized method in estimating missing daily mean air temperature[J].Meteor Mon,38(9):1135-1139.(in Chinese).
张璐,田向军,刘宣飞,等,2017.基于多重网格策略的NLS-3DVar资料融合方法及其在气温数据融合中的应用[J].气候与环境研究,22(3):271-288. Zhang L,Tian X J,Liu X F,et al.,2017.NLS-3DVar data fusion method based on multigrid implementation strategy and its application in temperature data fusion[J].Clim Environ Res,22(3):271-288.doi:10.3878/j.issn.1006-9585.2016.16140.(in Chinese).
中国气象局,2003.地面气象观测规范[M].北京:气象出版社:121. China Meteorological Administration,2003.Specifications for surface meteorological observation[M].Beijing:China Meteorological Press:121.(in Chinese).
周笑天,褚希,姚志平,2012.一种基于k-means聚类的实时气温动态质量控制方法[J].气象,38(10):1295-1300. Zhou X T,Chu X,Yao Z P,2012.A dynamic method of quality control for real-time temperature measurements based on k-means clustering algorithm[J].Meteor Mon,38(10):1295-1300.doi:10.7519/j.issn.1000-0526.2012.10.016.(in Chinese).
周笑天,陈益玲,李芸,等,2022.一种基于特征曲线的自动土壤水分观测数据异常值检测方法[J].中国农业气象,43(3):229-239. Zhou X T,Chen Y L,Li Y,et al.,2022.An outliers detection method for automatic soil moisture observation data based on characteristic curve[J].Chin J Agrometeorol,43(3):229-239.doi:10.3969/j.issn.1000-6362.2022.03.006.(in Chinese).