气象资料序列的完整性是开展气候分析与评估的必要条件之一。随着我国精密气象监测系统建设的不断推进、地面气象观测站网布局逐步优化，区域观测盲区得到进一步消除，观测要素短板也实现了补足，为开展中小尺度灾害性天气监测预警业务和区域小气候特征分析创造了有利条件。与此同时，随着大量无人值守的地面加密自动气象站的布设，因仪器故障、通信中断、自然灾害等原因引起的观测中断或数据质量异常的发生概率大增，不断出现随机站点和随机长度这种双随机特点的气象资料序列缺失，给气候分析和业务应用造成了不小的障碍。

气温作为研究气候演变最基础的物理量，其日值序列的完整性和准确性对于气候统计与分析有着重要意义，也是业内关注的焦点之一（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=[y₁，y₂，···，y_N]是由气温观测值构成的时间序列，长度为N，则可建立一元线性回归方程，记作:

式中:t_i为y_i所对应的时间，在日时间尺度下，时间自变量t_i为气温观测日期的日序（t₁=1，t₂=2，···，t_N=N），日期与日序一一对应; a为回归常数，b为回归系数，a和b用最小二乘法估计。$\widehat{Y}=\left[{\widehat{y}}_1，{\widehat{y}}_2，\cdots ，\right.\left.{\widehat{y}}_N\right]$为拟合时间序列，其形态为一条直线，代表了气温的整体变化趋势。

观测值y_i到拟合值${\widehat{y}}_i$之间的残差记为:

将（2）式变换形式可得

由（3）式可看出，观测值y_i可以由残差值e_i和拟合值${\widehat{y}}_i$相加构成，即序列Y=[y₁，y₂，···，y_N]可以分解出残差序列E=[e₁，e₂，···，e_N]和拟合序列$\widehat{Y}=\left[{\widehat{y}}_1，{\widehat{y}}_2，\cdots ，\right.\left.{\widehat{y}}_N\right]$，实现“残拟分离”; 将残差序列E=[e₁，e₂，···，e_N]和拟合序列$\widehat{Y}=\left[{\widehat{y}}_1，{\widehat{y}}_2，\cdots ，\right.\left.{\widehat{y}}_N\right]$的再次结合称为“残拟重构”。

1.2.2 气温插补区

设Y=[A，B，C]=[a₁，···，a_N，b_N+1，···，b_2N，c_2N+1，···，c_3N]为气温观测时间序列，长度为3N（N≥1），序列B=[b_N+1，···，b_2N]（长度为N）为序列Y中气温连续缺失部分，元素为空值。如果序列B的左邻序列A=[a₁，···，a_N]和右邻序列C=[c_2N+1，···，c_3N]存在且元素完整，则称序列B是可插补的。插补后的序列B记为$\tilde{B}=\left[{\tilde{b}}_{N+1}，\cdots ，{\tilde{b}}_{2N}\right]$，插补后的序列Y记为:

序列A、B和C所在的时间区间分别称为A、B和C区，B区也被称为插补区，序列Y和$\tilde{Y}$所在的时间区间被称作计算区间，如图1所示。

B区（插补区）及其左邻A区和右邻C区的设定，扩展了《地面气象观测规范》（中国气象局，2003）中缺测数据内插仅限于单时次的限定，从而可以有效解决“双随机”中的长度随机问题。

1.2.3 DTW距离与参照站

动态时间规整（DTW）是时间规整与距离测度计算相结合的一种非线性规整计算方法，常被用在语音识别系统中（李正欣等，2014; 闫宏宸和肖熙，2021）。该方法的最大特点就是可以通过路径规划来计算两个时间序列之间的最短累计距离，从而衡量两者的相似程度（Tormene et al.，2009; 周笑天等，2022）。

本文将动态时间规整算法原理应用于气温序列的相似度计算。设X=[x₁，x₂，···，x_N]和Y=[y₁，y₂，···，y_N]为两个不同站点的气温时间序列，长度为N，如图2所示，上下两条实线分别代表序列X和序列Y，动态时间规整算法的路径规划过程即为从左向右刻画虚线的过程，该过程以（x₁，y₁）为起点，以（x_N，y_N）为终点，按照单调、连续的匹配原则，每向右行进一步画一条虚线，虚线两端分别匹配上下序列中的一个元素。

设行进至第k步时，虚线两端匹配的元素分别为x_i和y_j，则该步的步长距离为x_i和y_j气温之差的绝对值，记为d（w_k）=|x_i-y_j|。

设L为起点（x₁，y₁）至终点（x_N，y_N）刻画虚线的总数量，则d（w_k）最短累计步长即为序列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距离记作:

由式（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）];

第4步，将参照站q的B^（q）直接嫁接于插补站p的A^（p）和C^（p）之间，用以替换B^（p），形成插补站p第一次插补序列${\tilde{Y}}^{（p）}=\left[A^{（p）}，{\tilde{B}}^{（p）}，C^{（p）}\right]，{\tilde{B}}^{（p）}=B^{（q）}$，完成一次插补;

第5步，分别对参照站q的Y^（q）=[A^（q），B^（q），C^（q）]和插补站p的一次插补序列${\tilde{Y}}^{（p）}=\left[A^{（p）}，{\tilde{B}}^{（p）}，C^{（p）}\right]$建立时间为自变量的一元线性回归方程，Y^（q）对应拟合序列${\widehat{Y}}^{（q）}=\left[{\widehat{A}}^{（q）}，{\widehat{B}}^{（q）}，{\widehat{C}}^{（q）}\right]$和残差序列$E^{（q）}=\left[E_A^{（q）}，E_B^{（q）}，E_C^{（q）}\right]$，${\tilde{Y}}^{（p）}$对应拟合序列${\widehat{\tilde{Y}}}^{（p）}=\left[{\widehat{A}}^{（p）}，{\widehat{\tilde{B}}}^{（p）}，{\widehat{C}}^{（p）}\right]$和残差序列${\tilde{E}}^{（p）}=\left[E_A^{（p）}，{\tilde{E}}_B^{（p）}，\right.\left.E_C^{（p）}\right]$，实现残拟分离;

最后，将序列$E_B^{（q）}$和序列${\widehat{\tilde{B}}}^{（p）}$相加，残拟重构得到序列$B_r^{（p）}$，$B_r^{（p）}$即为插补站p数据缺失序列B^（p）的二次插补结果。

从图3所示的插补方法总体流程可以看出，该方法采用了数据嫁接的操作方式，共分为两个插补阶段，一次插补为B^（q）对B^（p）的直接替换，二次插补是在一次插补的基础上，对两站的气温序列进行残拟分离和残拟重构后，得到的插补结果。

需要说明的是，插补流程中第3步是以DTW距离测度为基准选定参照站，而当以地理距离测度为基准（如水平距离最近或者海拔高度最近等）选定参照站时，则应以相应距离计算方法代替。

2 方法检验与分析

2.1 检验设计

为了验证方法对实况数据插补的有效性，同时为了使检验过程更加完备，本文对检验条件和检验内容设计如下:

1）检验包括日平均气温、日最高气温和日最低气温3种要素在内的气温日值数据。

2）根据山东省气象地理一级区划规则，分别在鲁西北（以平原地形为主）、鲁中（以山地地形为主）、鲁南（以平原丘陵地形为主）和半岛（以丘陵和山地海岸地形为主）4个地区的每个地区中随机选择10个无人值守的地面气象观测站点作为插补测试站（图4），在兼顾地形特征的同时满足站点分布的随机性要求。

3）每个插补测试站的插补区起止时间随机产生，以满足插补长度随机性要求。

4）检验采用观测真值置空的方式模拟插补区的缺失数据，并在插补结束后计算插补值与观测真值之间的误差以评估插补效果。具体的误差评估指标包括均方根误差（RMSE）:

和平均绝对误差（MAE）:

式中: ${\tilde{y}}_i$为插补值; y_i为观测真值; n为插补总样本数（站次数）。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℃）。

注：1）标注同一地区中RMSE指标最小值，2）标注同一地区中MAE指标最小值.

注：1）标注同一地区中RMSE指标最小值，2）标注同一地区中MAE指标最小值.

注：1）标注同一地区中RMSE指标最小值，2）标注同一地区中MAE指标最小值.

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）该方法对地形有一定的敏感性，平原或丘陵地区的插补效果要优于山地地区。

本文提出的基于时间序列的数据二次插补机制，对解决日益增加的双随机特点的气象资料缺失问题，有着较为广阔的应用前景，也为历史长序列气象数据的均一化订正提供了参考。

参考文献