1 前言

温度是农业、生态系统研究非常重要的观测量, 农业作物生长的模拟、农业气象灾害监测、生态系统模拟中温度是必不可少的输入^{[1, 2]}. 随着农业、生态模拟的精细化, 要求温度数据具有更高的精度, 如农业气象灾害干热风监测、林块生态系统碳排放的模拟等^{[3, 4]}, 高精度的温度观测量是必不可少的输入参数. 温度观测数据一般通过野外气象观测站获取, 由于设备故障、恶劣环境或是认为操作失误等原因, 小气象观测难免会出现缺失^[5], 缺失数据插补或补全, 是温度观测数据运用前必不可少的预处理工作.

本文针对中国广州一个森林生态站长时间间隔温度观测数据缺失进行插补方法研究. 由于该森林生态站处在雷雨区, 小气象站夏季容易因恶劣天气损坏, 容易造成较长时间的数据缺失. 论文选择了该生态站同时有自动观测气象站, 作为比对观测, 该小气象场还有人工温度观测设施. 自动观测气象数据输出的观测产品时间频率为30分钟, 每天有48条观测记录数据; 人工观测分早、中、晚每天3次, 产生3条记录. 本文研究的实际应用问题: 如何运用不同数据插补方法, 通过每天低频的人工温度观测数据, 获取完整的高精度半小时频率温度观测数据.

数据插补是众多学科数据分析前必不可少的预处理工作. 目前以数据插值、统计分析和时间序列分析等为基础, 发展出了多种数据插补方法^[6,7], 但对高精度数据插补研究还处初步阶段^[5]. 高精度数据插补的要求是从已观测数据中学习到数据的规律或缺值模式, 从而实现对未观测数据的准确估计. 深度学习是机器学习领域一个新的研究方向, 是人工智能领域的一项颠覆性技术创新, 除了带来图像、语音和自然语言处理领域的突破, 也成功应用到了众多学科领域^[8,9]. 深度学习旨在获得样本数据的内在规律和表示^[10], 和数据插补的需求非常契合.

深度学习技术已经在交通、医疗、传感器网络等多个领域的数据插补中获得了成功应用^[11-13], 并发展出了GRU、LSTM、GAN等不同结构的数据插补深度学习神经网络^[14]. GRU (gate recurrent unit)和LSTM (long short-term memory)都是循环神经网络的不同形式, 可以解决RNN网络学习过程中的梯度消失或爆炸问题^[15]. 这两种结构神经网络在数据插补应用中, 不但可以从已观测数据之中学习规律, 也可从数据缺失值模式中进一步学习, 提高数据插补精度^[16,17]. GAN (generative adversarial networks)网络用于学习多变量时间序列的总体分布, 从而对观测数据中的缺失值进行插补^[18].

在众多深度学习时间序列数据插补模型中, 一种基于序列-序列(Seq2Seq)的深度学习模型在多个标准样本集数据插补都有很好的表现^[19], 该结构采用了双向循环LSTM网络, 在随后的实际应用也进一步验证了该结构适用于时间序列数据缺失插补问题^[20]. 另外, 基于Encoder-Decoder结构的深度学习神经网络在数据插补方面也获得成功应用^[21]. 这些不同结构的深度学习模型为本文研究提供了重要参考.

本文运用低频人工温度观测数据, 来插补高频次机器观测数据的长时间观测值缺失问题. 时间序列数据缺失值插补虽然已很丰富, 但针对这一特定应用场景的数据插补方法研究文献还是较少^[22]. 下面是一个具体的生态台站的小气象观测数据, 该小气象站同时具有半小时自动温度观测和每天3次人工温度观测, 但半小时自动温度观测数据存在较长时间的缺失值. 为了实现半小时温度观测数据的高精度插补, 文中详细给出了一个编码-解码结构的序列-序列深度学习温度插补模型的构建过程和数学公式, 并将其与其他插补方法进行了插补精度对比分析.

2 研究数据介绍

本文研究采用了我国广州鼎湖山森林生态系统国家野外科学观测研究站的气象温度观测数据. 该生态系统观测站开展有温度观测对比试验, 同时开展人工观测和气象机器自动观测活动, 有较长时间的温度观测数据记录, 表1是用于本文研究的温度人工观测数据和自动机器观测数据情况.

由于鼎湖山生态站位于中国南方的山区, 自动观测设备容易受雷雨季节影响而产生较长时间观测记录的缺失. 图1是某一机器自动气象观测数据的数据缺失情况分布图, 从图中可见在2020年7月有一次超过2个月的温度观测数据缺失.

运用每天3次的人工观测数据对2020年7月超过2个月的机器温度观测数据缺失进行插补, 是本文方法研究的应用背景. 由于人工观测数据和机器观测数据之间很容易建立线性相关关系, 本文数据插补问题研究核心是如何运用低频的人工温度观测数据对高频的机器自动温度观测数据进行插补.

3 温度观测数据深度学习插补 3.1 基本定义

气温具有天的周期性, 很自然的将长时间序列温度观测数据按照天为单位进行划分, 变为每天48个观测值的分段序列. 为了更好的把研究集中到长时间间隔缺失值的插补, 对时间序列中偶尔或短时间数据缺先采用上文的卡尔曼平滑的方法进行插补. 这样以天为分段的温度时间序列包含两种, 即没有缺失值的每天分段, 记为 ${d}_{\rm full}^{j}$ , 和只包含早中晚3次观测值的每天分段, 记为 ${d}_{\rm miss}^{j}$ , 这样以天为分段的时间分段序列可以表示为:

$ {\text{\{ }}d_{\rm full}^1, \cdots, d_{\rm full}^i, d_{\rm miss}^{i + 1}, \cdots, d_{\rm miss}^{i + m}, d_{\rm full}^{i + m + 1}, \cdots, d_{\rm miss}^n{\text{\} }} $

(1)

序列(1)表示长度为n天温度时间序列中, 缺失值窗口宽度为m天. 半小时温度观测序列(1)为长度为48n, 存在缺失值48m的半小时温度观测数据序列, 以天分段后的表达形式.

为了表示时间序列(1)中缺失值的位置, 对长度为L(48n)的半小时采样温度时间序列{ $ {T}_{t}^{i} $ }, 构建相应长度为L的掩码时间序列{ $ {m}_{t}^{i} $ }, 其中:

$ m_{t}^{i}=\left\{\begin{array}{l}0 , \;{\text{如}} T_{t}^{i} {\text{未被观测}} \\ 1,\;{\text{否则}} \end{array}\right.$

现在以天为单位, 对长度为L的半小时掩码序列进行分段, 没有缺失值的掩码每天分段记为 ${M}_{\rm full}^{j}$ , 和只包含早中晚3次观测值的掩码每天分段, 记为 ${M}_{\rm miss}^{j}$ , 这样就可以建立与式(1)对应的以天为分段的掩码序列:

$ \{ M_{\rm full}^1, \cdots, M_{\rm full}^i, M_{\rm miss}^{i + 1}, \cdots, M_{\rm miss}^{i + m}, M_{\rm full}^{i + m + }, \cdots, M_{\rm full}^n\} $

(2)

3.2 滚动窗口采样

采用滚动窗的方法, 基于以天为分段的时间序列为深度学习模型训练构建样本集. 对长度为m(天)的缺失值进行插补, 需构建样本滚动窗口的长度大于m, 并且在m的两端各保留长度为s(天)的观测数据, 这样滚动窗口长度w为m+2×s天. 训练样本为适应序列-序列(Seq2Seq)的训练方法来构建, 对长度为w的训练输入样本中温度观测序列为:

$ {\text{\{ }}d_{\rm full}^1, \cdots, d_{\rm full}^s, d_{\rm miss}^{s + 1}, \cdots, d_{\rm miss}^{s + m}, d_{\rm full}^{s + m + }, \cdots, d_{\rm miss}^w{\text{\} }} $

(3)

可通过训练形成如下的时间序列结果输出:

$ {\text{\{ }}d_{\rm full}^1, \cdots, \hat d_{\rm full}^s, \hat d_{\rm impt}^{s + 1}, \cdots, \hat d_{\rm impt}^{s + m}, d_{\rm full}^{s + m + }, \cdots, d_{\rm miss}^w{\text{\} }} $

(4)

序列(4)中, ${\hat d}_{\rm impt}^{j}$ 为对缺失值插补后, 某一天每半小时温度观测数据值的完整分段. 在构建训练样本时, 按照观测值与掩码序列以天为顺序对应关系, 构建长度为w(天)的掩码序列, 作为训练样本的另一个输入.

训练样本需要在没有缺失值的温度观测序列基础上构建, 样本中观测值缺失的模式同实际情况一致, 即每天只有早中晚3次观测值. 表2是训练样本中存在缺失值的某一天温度数据及其对应的掩码示例.

3.3 深度学习模型的设计

典型的基于Seq2Seq的时间序列数据插补深度学习模型有SSIM, BRTS-I等^[19,21]. 本文吸收了这些模型的优点, 将Seq2Seq和Encoder-Decoder深度学习架构结合起来, 所设计的深度学习模型结构如下文所述. 上面输入序列(1)被记为x={x₁, x₂, …, x_n}, 输出序列(4)被记为y={y₁, y₂, …, y_n}, 掩码序列(2)被记为m={m₁, m₂, …, m_n}.

(1)编码

从图2中可见, 深度学习结构中的编码部分的基本结构是LSTM-I, 该结构与BRTS-I结构中的RTS-I结构相似, 其中的循环神经网络单元直接采用了长短期记忆单元; 另外, 本文温度观测缺失值部分, 每天48个半小时温度值, 只有3个观测值, 有效值比较稀疏, 所有没有采用RTS-I中的缺失值时间间隔的变量和相应的训练公式. 下面定义中LSTM被简化为一个简单算子的形式, 将LSTM-I单元过程数学描述为:

$ {\tilde x_t} = {W_x}{h_{t - 1}} + {b_x}$

(5)

$ x_t^c = {x_t} \odot {m_t} + (1 - {m_t}) \odot {\tilde x_t} $

(6)

$  {h}_{t}=LSTM({x}_{t}^{c}, {h}_{t-1}) $

(7)

$ l_{t}=\left\langle m_{t}, £\left(x_{t}, \tilde{x}_{t}\right)\right\rangle $

(8)

式(5)将上一个LSTM单元的隐状态h_t–1转化为估计向量 $\tilde{x}_{t}$ , 其中W_x、b_x为模型参数; 式(6)通过运用掩码向量m_t, 把输入向量x_t中的缺失值替换为估计向量 $\tilde{x}_{t}$ 对应的值; 式(7)通过LSTM网络单元把 $ {x}_{t}^{c} $ 和隐状态h_t–1产生预测状态h_t; 式(8)是LSTM-I单元的估计误差, 为缺失值位置上观测值与估计值绝对差的累计量.

图2中神经网络编码部分由双向的LSTM-I神经网络构成: 一个是从时间序列的开始到结束读取输入, 产生前向隐状态向量序列 $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over h} = \{ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over h} _1},{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over h} _2}, \cdots ,{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over h} _n}\}$ ; 另一个是从时间序列的结束到开始反向读取输入, 产生后向隐状态序列 $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftharpoonup$}} \over h} = \{ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftharpoonup$}} \over h} _1},{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftharpoonup$}} \over h} _2}, \cdots ,{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftharpoonup$}} \over h} _n}\}$ . 将前向和后向隐状态序列拼接到一起, 构成编码层的编码输出h={h₁, h₂, …, h_n}, 其中向量h_i为:

$ {h}_{i}=\text{ }\{{\overrightarrow{h}}_{i}, {\overleftarrow{h}}_{i}\}   $

(9)

双向编码LSTM-I编码网络误差包括正向和逆向估计误差两部分.

(2)解码

解码层接收编码输出序列h, 产生插补后的时间序列结果序列y. 本文通过优选, 采用了2种解码结构: 一种是全连接层和一层线性层组合, 如图3所示; 另一种是LSTM和一层线性层的组合, 如图4所示.

① 基于两层线性层组合解码过程数学描述如下:

$ {\tilde h_t} = Dropout({h_t}) $

(10)

$ {s_t} = g({W_s}{\tilde h_t} + {b_s}) $

(11)

$ {y_t} = {W_y}{s_t} + {b_y} $

(12)

$ l_{y}=\left\langle{m}_{{t}}, £\left({x}_{{t}}, {y}_{{t}}\right)\right\rangle $

(13)

式(10)中的Dropout通过对解码输入随机丢弃部分神经元, 能够起到预防过拟合的作用; 式(11)为全连接层, $ g $ 为激活函数, 全连接层产生输出状态序列s={s₁, s₂, …, s_n}; 由于温度值是连续值, 式(12)为最上层即线性变换层, 输出插补结果序列y; 式(13)是解码的插补结果误差, 为缺失值位置上观测值与插补值绝对差的累计量结果.

② LSTM层和线性层组合的解码结构

LSTM层和线性层组合的解码结构如图4所示. 解码层接收编码层的输出序列h, 产生插补后的时间序列结果序列y.

LSTM层和线性层组合的解码过程数学描述为:

$ {s_t} = LSTM\left( {{h_t}, {s_{t - 1}}} \right) $

(14)

$ {y}_{t}={W}_{y}{s}_{t}+{b}_{y}  $

(15)

$l_{y}=\left\langle m_{t}, £\left(x_{t}, y_{t}\right)\right\rangle$

(16)

如式(14), 解码层底部是一个标准的LSTM网络, 该网络综合编码输出序列h, 产生包含更丰富信息输出状态序列s={s₁, s₂, …, s_n}; 如式(15), 由于温度值是连续值, 解码层顶部采用了线性全连接层, 输出插补结果序列y. 同式(13), 式(16)是解码层的插补结果误差.

上述两种解码机制的数据插补深度学习网络的误差构成是相同的, 神经网络的误差包括3部分, 即:

$ {l_{\text{t}}} = l_t^f + l_t^b + {l_y} $

(17)

式(17)中, $ {l}_{t}^{f} $ 前向LSTM-I编码层的估计误差, $ {l}_{t}^{b} $ 为后向LSTM-I编码层的估计误差.

4 插值效果评估方法

本文采用多个指标评价不同数据插补方法的性能, 评价指标的数值在测试样本集上计算. 包括均方根误差(RMSE), 平均绝对误差(MAE), 平均相对误差(MRE)和皮尔逊相关系数(PCC), 定义如下:

$ \textit{RMSE} = \sqrt {\frac{1}{n}\sum\nolimits_{i = 1}^n {{{({x_i} - {y_i})}^2}} } $

(18)

$ MAE = \frac{1}{n}\sum\limits_{i = 1}^n {|{x_i} - {y_i}|} $

(19)

$ MRE = \frac{1}{n}\sum\limits_{i = 1}^n {|\frac{{{x_i} - {y_i}}}{{{x_i}}}|} $

(20)

$ PCC = \frac{{\displaystyle\sum\limits_{i = 1}^{{n}} {({x_i} - \bar x)({y_i} - \bar y)} }}{{\displaystyle\sum\limits_{i = 1}^n {{{({x_i} - \bar x)}^2}\displaystyle\sum\limits_{i = 1}^n {{{({y_i} - \bar y)}^2}} } }} $

(21)

上面指标计算中, x_i为所构造测试本中缺失值位置的实际观测值, y_i为缺失值位置对应的插补结果值. PCC评价指标计算如式(21), 其中 $\bar x$ 为样本中缺失值位置实际观测值的总体平均, $\bar y$ 为缺失值位置插补结果的总体平均.

5 方法结果

作为对比, 本文选用了简单的总体平均插补方法(Mean); 基于时间序列分解的卡尔曼插补方法(Kalman-struct); 基于深度学习的BRTS-I时间序列插补方法; 本文设计的两种编码-解码结构的双向LSTM网络插补方法, 解码层为全连接(BiLSTM-FC-I)和解码层为LSTM (BiLSTM-LSTM-I).

上述方法中, 总体平均插补方法可以在整个数据集上计算获取, 而插补方法BRTS-I、BiLSTM-FC-I、BiLSTM-LSTM-I、Kalman-struct均需要先把整个数据集分为训练集和测试集, 然后在相同的训练集上进行训练, 在同一测试集上进行精度分析. 基于深度学习的BRTS-I、BiLSTM-FC-I、BiLSTM-LSTM-I虽然模型结构有所不同, 但其关键参数LSTM状态的维度均相同, 均取值为108; 另外模型的训练参数也相同, 采用了相同的mini batch参数和优化方法, 优化方法均选用Adam, 初始学习率取值为0.001; 训练终止策略均为连续10轮测试精度均不超过训练过程中测试精度的最优值. 深度模型实现是以PyTorch深度学习框架为基础, 以2020年07月缺失值窗口左侧的观测数据构建训练集, 右侧的观测数据构建测试集. 深度学习插补方法构建了两种训练样本, 一种缺失值时间窗口设定为30天, 另一种缺失值时间窗口设定为60天, 这两种训练样本缺失值窗口前后连续观测值均设定为14天. 为了区别这两种训练样本, 下面将缺失值时间窗口长度作为了深度学习插补方法的后缀, 将各种插补方法的结果总结如表3.

从表3中可见, 深度学习方法要明显优于简单的总体平均方法. 深度学习方法之间的精度也存在较大的差别, 图5是卡尔曼插值方法, 以及各种深度学习方法插值RMSE精度的比较图.

从图5中比较BRTS-I、卡尔曼方法、BiLSTM-I三种方法, 本文设计的两种BiLSTM-I深度学习气温插补方法要优于其他两种方法; 基于LSTM解码的BiLSTM-I深度学习方法精度优于基于全连接解码的BiLSTM-I深度学习方法; BRTS-I深度学习时间序列插补方法精度最低.

卡尔曼平滑方法时间序列插补方法精度取决于状态方程是否准确表达了时间序列, Kalman-struct假定时间序列的趋势性和季节性的成分可通过基本的线性方程进行拟合. 深度学习方法, 没有对时间序列的表达形式进行任何假设, 而是通过训练数据集, 自动学习时间序列的准确表达形式. 从测试结果看, BiLSTM-I方法比Kalman-struct卡尔曼方法, 以及BRTS-I深度学习方法更有可能获取了时间序列的准确表达, 从而获得了更高的数据插补精度.

本文设计的两种解码结构的BiLSTM-I模型插补精度都高于BRTS-I模型. BiLSTM-I模型与BRTS-I模型的区别主要有两点: 首先从模型结构上, BiLSTM-I采用了Encoder-Decoder结构, 而BRTS-I只相当于BiLSTM-I模型的Encoder部分, 模型结构上BiLSTM-I有利于充分学习数据的潜在分布规律, 从而可以提高数据插补精度; 其次模型代价函数有区别, BiLSTM-I和BRTS-I的代价函数均由3部分构成^[19], 前两部分是相同的, BRTS-I模型代价函数的第3部分是前向和后向LSTM网络对缺失估计值的差; BiLSTM-I模型代价函数的第3部分则为最后的估计值和真实观测值的差, BiLSTM-I模型误差函数对插补结果的评价更直接, 模型收敛误差和插补精度直接对应, 从而确保模型收敛时插补误差也能达到最小.

两种解码结构的BiLSTM-I模型, 基于LSTM解码的模型插补精度优于全连接解码结构模型插补精度, 这主要是LSTM解码时, 不但可以利用当前的编码输出信息, 还可以利用之前的编码输出信息; 而全连接编码则只能利用当前的编码信息, LSTM解码方法对编码信息的利用更为充分.

图5中两种解码结构的BiLSTM-I模型, 缺值窗口分别为30天和60天的测试精度基本一致. 深度学习方法应用中模型的泛化能力非常重要, 本文问题体现在模型对不同宽度缺失值窗口插补精度是否一致. 为了对这一点进行检验, 我们用缺失值窗口为30天的模型对缺失值为60天时间温度观测时间序列进行插补, 然后用缺失值窗口为60天的模型对缺失值为30天时间温度观测时间序列进行插补, 表4是这两种情况下, 两种解码结构的插补方法结果的精度统计表.

从表4可见, 无论是缺失值窗口为60天的模型应用到缺失值情况为30天, 还是缺失值窗口为30天模型应用到缺失值情况为60天, 两种解码结构的模型精度的各项指标都非常稳定, 这表明这两种解码结构的BiLSTM-I深度学习模型都对不同缺失值窗口有较好的泛化能力.

6 总结

本文运用不同深度学习数据插补方法, 通过每天低频的人工温度观测数据, 获取完整的高精度半小时频率温度观测数据. 本文采用序列-序列的时间序列插补方法, 基于编码-解码结构的深度学习模型(BiLSTM-I), 编码层采用双向LSTM-I网络, 解码层分别采用LSTM解码结构与全连接解码结构, 设计了两种解码结构的深度学习数据插补模型.

试验分析结果表明, 本文设计的BiLSTM-I深度学习气温插补方法要优于其他方法. 基于LSTM解码结构的深度学习模型, 缺失值时间为30天的测试集, 测试结果精度RMSE为0.47℃; 缺失值时间为60天的测试集, 测试结果精度RMSE为0.49℃. 基于全连接解码结构的深度学习模型, 缺失值时间为30天的测试集, 测试结果精度RMSE为0.76℃; 缺失值时间为60天的测试集, 测试结果精度RMSE为0.70℃.

最后, 文章还分析了BiLSTM-I深度学习插补方法对不同时间温度缺失长度的适应能力. 分别用缺失值时间长度为30天的训练模型, 对缺失值为60天的测试集进行插补; 以缺失值时间长度为60天的训练模型, 对缺失值为30天的测试集进行插补, 结果表明两种解码结构的深度学习训练模型对不同的温度缺失时间长度具有泛化能力.