由于传染病可在大范围内传播，且传播速度快，传染病成为危及人类身体健康的主要因素之一，给人类带来了巨大的影响与危害. 根据监测到的传染病数据对其传播进程进行预测，从而采取应对措施，将大大降低传染病的危害，因此对传染病传播过程中的感染疾病人数进行预测，对传染病防控意义重大^[1-2].

目前，在预测传染病的模型和方法中，应用较为广泛的主要有传播动力学模型^[3-4]、马尔可夫模型^[5]以及时间序列方法^[6]等. 有学者针对传染病传播预测的问题提出了传播动力学模型，Kermack等^[3]提出的经典的易感-感染(susceptible-infected，SI)模型，可以合理地预测传染病在人与人之间的传播. 该模型后来被扩展到更复杂的情况，如易感-感染而未发病-感染-恢复-易感(susceptible-exposed-infected-recovered-susceptible，SEIRS)模型^[4]. Gomez等^[5]提出基于接触式传染病传播的离散马尔可夫链方法，将单个节点（个人）的接触过程推广到多个节点的复杂情况. 时间序列方法在传染病预测方面也有着很大的贡献，例如，Benvenuto等^[6]采用差分整合自回归移动平均(auto regressive integrated moving average，ARIMA)模型预测了新型冠状病毒疾病(coronavirus disease 2019，COVID-19)的流行病学趋势和发病率. 传染病传播数据是一种典型的时空数据，在时间和空间上都呈现出关联特性. 在时间上，每一时刻的感染疾病人数都与其之前时刻的感染人数有关；在空间上，通常邻近且相关因素相似的地区的发病水平和传播特征相关性较高，反之则相关性较低^[7]. 但上述方法主要关注传染病传播过程中的时间演变特性，而忽略了其空间相关特性. 此外，空间聚集性研究利用病例的发病时间、地理位置、人口数据以及传染病历史数据等信息，分析传染病传播过程中具有聚集性的时间区间和空间区域，从而对传染病是否爆发进行估计^[7]. 但在传染病传播过程中，地区间的空间相关性，即地区之间的人群流动也是一个非常关键的影响因素^[8]，空间聚集性研究则主要关注传染病的时空聚集性，忽略了地区间的空间相关性.

图信号处理(graph signal processing, GSP)利用图模型上的节点和边来刻画数据之间的关系, 将数据建模为图信号，例如社交网络图^[9]、传感器网络图^[10]，以对非规则网络数据进行处理. 因此，GSP逐渐受到研究者们的关注，也逐渐应用于传染病传播的相关研究^[11-13]. 基于GSP的相关理论，将传染病数据建模为时变图信号，将地区建模为图上节点，则可以刻画出数据间的相关性^[14]. 例如，Isufi等^[13]利用近似时间-顶点平稳性假设，采用基于图的时间序列预测模型对传染病的传播进行预测，从而反映出传染病传播过程中的时空关联性. 但是，Isufi等^[13]所提出的方法在对模型参数进行估计时仅考虑了预测时刻与其之前时刻信号的相关性，而未考虑到在传染病传播过程中，每一时刻的信号都与其之前时刻信号存在相关性，因此未能充分挖掘时序数据间的相关性，导致其方法不具有普适性，对于待预测的数据类型要求较高.

为了能够更好地描述传染病数据的空间相关性和时间相关性，从而更准确地刻画传染病的传播特性，采用图多项式-向量自回归(graph polynomial-vector autoregressive，GP-VAR)模型^[13]对传染病的传播进行预测，并提出了新的模型参数估计方法. 将传染病传播数据建模为时变图信号，将地区建模为图上节点，并采用GP-VAR模型预测时间序列在图上的演变过程. 考虑到预测时刻前的每一时刻与其之前所有时刻信号的相关性，提出基于最小二乘(least squares，LS)优化的模型参数估计方法. 实验部分使用SI模型、SEIRS模型生成的模拟数据^[15]以及德国COVID-19数据对感染疾病的人数进行预测，将基于LS优化的GP-VAR模型与原有的GP-VAR模型^[13]、图-向量自回归移动平均(graph-vector autoregressive moving average，G-VARMA)模型^[13]以及ARIMA模型^[6]进行比较，将归一化最小均方误差根(root normalized mean-square error，rNMSE)作为评价指标^[13]，以验证所提出方法的有效性.

图信号处理作为分析非规则域数据的有力工具，以图模型和图信号为基础，将代数信号处理理论^[16]、谱图理论^[17]与计算谐波分析相结合，并将一些经典信号处理技术扩展到图域，如图傅里叶变换(graph Fourier transform，GFT)^[17]、图滤波器^[18]，以解决非规则和复杂结构的信号处理问题.

给定一个无向加权图 $ G = \left\{ {V,E,{\boldsymbol{W}}} \right\} $. 其中， $ V $为图上节点的集合，共有 $ N $个节点， $ V = \left\{ {1,2, \cdots ,N} \right\} $； $ E $为图上节点之间边的集合； $ {\boldsymbol{W}} $为该图的加权邻接矩阵，表示节点间的连接关系，若节点 $ i $和节点 $j$之间有边相连，则邻接矩阵 $ {\boldsymbol{W}} $中的元素 $ {{W}}\left( {i,j} \right) = $ $ {{W}}\left( {j,i} \right) \ne 0 $，否则为0. 度矩阵 $ {\boldsymbol{D}} = {\text{diag}}\;\left( {{d_i}} \right) $， $ {d_i} = $ $ {\displaystyle \sum\limits}_{j = 1}^N {{{W}}(i,j)} $，其主对角线上的元素为邻接矩阵 $ {\boldsymbol{W}} $的行和. $ {\boldsymbol{L}} = {\boldsymbol{D}} - {\boldsymbol{W}} $为 $ G $的组合拉普拉斯矩阵，归一化拉普拉斯矩阵由 $ {{\boldsymbol{L}}_{{\rm{nor}}}} = {\boldsymbol{I}} - {{\boldsymbol{D}}^{ - 1/2}}{\boldsymbol{W}}{{\boldsymbol{D}}^{ - 1/2}} $表示，其中 $ {\boldsymbol{I}} $为单位矩阵. 对 $ {\boldsymbol{L}} $进行特征分解有 ${\boldsymbol{L}} = {\boldsymbol{U}}{{\textit{Λ}}}{{\boldsymbol{U}}^{\text{T}}}$, 其中 ${{\textit{Λ}}}$为 $ {\boldsymbol{L}} $的所有特征值组成的对角阵，特征值通常被称为图频率， $ {\lambda _1} \leqslant {\lambda _2} \leqslant \cdots \leqslant {\lambda _N} $， $ {\boldsymbol{U}} $为由特征值所对应的特征向量组成的矩阵.

图信号可以用向量 $ {\boldsymbol{x}} = {[{x_1}, \cdots ,{x_N}]^{\text{T}}} $表示，其中 $ {x_i} $为节点 $ i $的信号值. 如图1所示为社交网络图及图上信号，节点颜色的深浅表示该点信号值的大小. 在图信号处理领域中，图傅里叶变换^[17]作为图信号的频域分析工具，其定义与传统信号处理中的傅里叶变换相似，信号 $ {\boldsymbol{x}} $的GFT定义如下：

(1) $ {\text{GFT }}\left\{ {\boldsymbol{x}} \right\} = {\boldsymbol{\hat x}} = {{\boldsymbol{U}}^{\text{T}}}{\boldsymbol{x}}. $

同时，也可以得出逆GFT的表达式为 $ {\boldsymbol{x}} = {\boldsymbol{U\hat x}} $, 根据GFT的定义可以在图频域中对图信号进行滤波.

由于许多实际应用中的数据具有时变特性，例如传染病传播数据会随时间变化，可以将其建模为时变图信号. 由此将图信号处理扩展到时变图信号处理，以对具有时变特性的数据进行分析^[19]. 时变图信号是指随时间变化的图信号，可以表示为包含 $ P $个连续时刻信号的矩阵 ${\boldsymbol{F = }}\left[ {{\boldsymbol{x}}_{{0}}},{{\boldsymbol{x}}_{{1}}},{{{{\boldsymbol{x}}}}_{{2}}},\cdots , $ $ {{\boldsymbol{x}}_{P - 1}} \right]$. 如图2所示为有向周期时间序列图，其反映了周期时间序列的因果关系，每个节点对应一个时刻^[20].

图滤波器的定义与传统信号处理中的滤波器相似, 通过对图信号进行滤波从而得到期望的信号. 图滤波器^[18]可以定义为关于拉普拉斯矩阵 $ {\boldsymbol{L}} $的函数：

(2) $ h\left({\boldsymbol{L}}\right)={{\displaystyle \sum} _{z=0}^{Z}{a}_{z}{\boldsymbol{L}}^{z}}. $

式中: $ {a_z} $为系数, $ Z $为 $h\left( {\boldsymbol{L}} \right)$的最高阶数. 信号 $ {\boldsymbol{x}} $经图滤波器 $h\left( {\boldsymbol{L}} \right)$后的输出信号为

(3) $ {\boldsymbol{y}} = h\left( {\boldsymbol{L}} \right){\boldsymbol{x}} = {\boldsymbol{U}}h\left( {\boldsymbol{\varLambda }} \right){{\boldsymbol{U}}^{\text{T}}}{\boldsymbol{x}}. $

式中: $h\left( {{\textit{Λ}}} \right)$为关于特征值矩阵 ${{\textit{Λ}}}$的函数， $h\left( {{\textit{Λ}}} \right) \in {{\bf{R}}^{N \times N}}$. 可以看出，图滤波器 $h\left( {\boldsymbol{L}} \right)$与拉普拉斯矩阵 ${\boldsymbol{L}}$具有相同的特征向量, 其特征值为关于 ${\boldsymbol{L}}$的特征值的函数. 拉普拉斯矩阵 ${\boldsymbol{L}}$反映了节点与其 $1 {\text{-}} hop$邻居（1-hop是指一个节点到另一个节点的最短路径须经过的边数为1）之间的差异， ${{\boldsymbol{L}}^z}$则融合了每个节点的 $z {\text{-}} hop$邻居的信息，故 $h\left( {\boldsymbol{L}} \right)$作为核矩阵可以反映图节点之间的相关性.

已知历史时刻的传染病数据，可以对当前或未来时刻的病例人数进行预测. 基于图信号处理的相关理论，将传染病的发病地区建模为图上节点，节点间的边及其权重可以反映传染病传播过程中的空间关联性. 例如，在构建德国COVID-19数据的图模型时，图上的边可以根据地区间的位置信息确定，将距离较近的2个城市用边连接，若地区之间人群流动较大，则节点间边的权重设置为 $ {{W}}\left( {i,j} \right) = {{W}}\left( {j,i} \right) = 1 $. 由于病例人数会随时间变化，可以将每个地区不同时刻的病例人数建模为该节点的时变图信号. 由此，采用基于图的时间序列模型对传染病的传播进行建模，从而反映出传染病传播过程中的时间演变特性和地区间的空间相关性.

一般而言，节点的局部邻域内的图信号值通常是相互关联的，即 $t$时刻节点 $i$的信号值 $ {x_i}\left( t \right) $不仅与 $ t - k\left( {1 \leqslant k \leqslant t} \right) $时刻节点 $ i $的信号值 $ {x_i}\left( {t - k} \right) $相关，而且与其相邻节点 $ j $的信号值 $ {x_j}\left( {t - k} \right) $相关，这些相关性可以用于时变图信号处理中^[21-22].

GP-VAR模型是时间序列在图上的向量自回归(vector autoregressive, VAR)模型，其系数矩阵是图滤波器，由此可以在预测时融合更多邻居节点的信息^[14]. 若将传染病数据建模为时变图信号，采用GP-VAR模型进行预测，将T时刻的图信号 $ {{\boldsymbol{x}}_T} $表示为其之前 $K - 1\left( {2 \leqslant K \leqslant T} \right)$个时刻的图信号 $ {{\boldsymbol{x}}_{T - K + 1}},{{\boldsymbol{x}}_{T - K + 2}}, \cdots ,{{\boldsymbol{x}}_{T - 1}} $的函数，表达式^[13]如下：

(4) $ {{\boldsymbol{\tilde x}}_T} = {\sum\limits}_{k = 1}^{K - 1} {{{\boldsymbol{A}}_k}{{\boldsymbol{x}}_{T - k}}} + {{\boldsymbol{\varepsilon }}_T}. $

式中： $ {{\boldsymbol{\tilde x}}_T} $为 $ T $时刻的预测信号； $ {{\boldsymbol{\varepsilon }}_T} $为误差向量，是一个均值为0、协方差矩阵正定的随机向量； $ {{\boldsymbol{A}}_k} $为图滤波器， ${\boldsymbol{A}}_{k}={{\displaystyle \sum} _{q=0}^{Q}{a}_{q,k} {\boldsymbol{L}}^{q}}$为反映相关性的系数矩阵，将图拓扑嵌入到模型中， $ Q $为 $ {{\boldsymbol{A}}_k} $的最高阶数.

在采用GP-VAR模型进行预测时，首先，通过优化问题对模型中图滤波器的参数 $ {a_{q,k}} $进行估计；其次，利用估计出的参数 $ {a_{q,k}} $和已知的拉普拉斯矩阵 $ {\boldsymbol{L}} $设计出图滤波器 $ {{\boldsymbol{A}}_k} $；然后，通过图滤波器 $ {{\boldsymbol{A}}_k} $、 $ K - 1 $个时刻的已知图信号 $ {{\boldsymbol{x}}_{T - k}} $以及误差向量 $ {{\boldsymbol{\varepsilon }}_T} $对 $ T $时刻的图信号进行预测；最后，将式(4)作为一步预测结果表达式，通过递归的方式对其计算 $n$次，可以得到 $n$步预测结果，即 $ T + n - 1 $时刻的预测信号.

基于前面对预测模型的分析，考虑到连续 $ K $时刻的图信号中每一时刻与其之前所有时刻的相关性，将 $ {{\boldsymbol{x}}_t} $表示为 $ t $时刻的图信号，将GP-VAR模型中图滤波器的参数 $ {a_{q,k}} $的估计问题归结为如下最小二乘优化问题:

(5) $\begin{split} {{\mathop {\min }\limits_{{a_{q,k}}} }}\;& {\displaystyle \sum _{t=1}^{K-1}{\left\| {\boldsymbol{x}}_{t}-{\displaystyle \sum _{k=1}^{t}{\boldsymbol{A}}_{k}{\boldsymbol{x}}_{t-k}}\right\|}_{2}^{2}};\;\\ &{\rm{s.t}}. \;\;{\boldsymbol{A}}_{k}={\displaystyle \sum _{q=0}^{Q}{a}_{q,k}{\boldsymbol{L}}^{q}}. \end{split} $

式中： ${\left\| \cdot \right\|}_{2}^{2}$表示2范数的平方.

式(5)可以归结为无约束优化问题:

(6) $ {\mathop {\min }\limits_{{a_{q,k}}} }\;\;{\sum\limits_{t = 1}^{K - 1} {\left\| {{{\boldsymbol{x}}_t} - \sum\limits_{k = 1}^t {\sum\limits_{q = 0}^Q {{a_{q,k}}{{\boldsymbol{L}}^q}} {{\boldsymbol{x}}_{t - k}}} } \right\|_2^2} } . $

该问题中的输入为信号 $ {{\boldsymbol{x}}_t} $、 $ {{\boldsymbol{x}}_{t - k}} $与拉普拉斯矩阵 $ {\boldsymbol{L}} $, 可以解得 $ (Q+1) \left( {K - 1} \right) $个参数 $ {a_{q,k}} $. 由于 $ \left\| {\boldsymbol{x}} \right\|_2^2 = $ $ {\rm{tr}}\;\left( {{\boldsymbol{x}}{{\boldsymbol{x}}^{\text{T}}}} \right) $, 则有

(7) $ \begin{split} &\sum\limits_{t = 1}^{K - 1} \left\| {{{\boldsymbol{x}}_t} - \sum\limits_{k = 1}^t {\sum\limits_{q = 0}^Q {{a_{q,k}}{{\boldsymbol{L}}^q}} {{\boldsymbol{x}}_{t - k}}} } \right\|_2^2 = \hfill \\ &\qquad\sum\limits_{t = 1}^{K - 1} {\left[ {{\rm{tr}}\;\left( {{{\boldsymbol{x}}_t}{\boldsymbol{x}}_t^{\text{T}}} \right) - 2\sum\limits_{k = 1}^t {\sum\limits_{q = 0}^Q {{a_{q,k}}} {\rm{tr}}\;\left( {{{\boldsymbol{x}}_t}{\boldsymbol{x}}_{t - k}^{\text{T}}{{\boldsymbol{L}}^q}} \right)} } +\right.} \hfill \\ &\qquad\left. {\sum\limits_{{k_1} = 1}^t {\sum\limits_{{k_2} = 1}^t {\sum\limits_{{q_1} = 0}^Q {\sum\limits_{{q_2} = 0}^Q {{a_{{q_1},{k_1}}}{a_{{q_2},{k_2}}}} {\rm{tr}}\;\left( {{{\boldsymbol{L}}^{{q_1}}}{{\boldsymbol{x}}_{t - {k_1}}}{\boldsymbol{x}}_{t - {k_2}}^{\text{T}}{{\boldsymbol{L}}^{{q_2}}}} \right)} } } } \right]. \end{split} $

令 $ {{\boldsymbol{r}}_{t,t - k}} $、 $ {{\boldsymbol{v}}_k} $、 ${{{\textit{φ}}}}_{t}$分别为 $ Q+1 $、 $ Q+1 $、 $ {(Q+1) t} $维的向量：

(8) $ \left. \begin{aligned} {{\boldsymbol{r}}_{t,t - k}} =\;& \left[ {{\rm{tr}}\;\left( {{{\boldsymbol{x}}_t}{\boldsymbol{x}}_{t - k}^{\text{T}}{{\boldsymbol{L}}^0}} \right),} \right.{\left. { \cdots ,{\rm{tr}}\;\left( {{{\boldsymbol{x}}_t}{\boldsymbol{x}}_{t - k}^{\text{T}}{{\boldsymbol{L}}^Q}} \right)} \right]^{\text{T}}}, \\ {{\boldsymbol{v}}_k} =\;& {\left[ {{a_{0,k}}, \cdots ,{a_{q,k}}, \cdots ,{a_{Q,k}}} \right]^{\text{T}}}, \\ {{{{\textit{φ}}}}_t} =\;& {\left[ {{\boldsymbol{v}}_1^{\text{T}},{\boldsymbol{v}}_2^{\text{T}}, \cdots {\boldsymbol{v}}_k^{\text{T}}, \cdots ,{\boldsymbol{v}}_t^{\text{T}}} \right]^{\text{T}}}. \end{aligned}\right\}$

由此可以将式(7)中的求和形式写成 ${(Q+1) t} $维行向量与 $ {(Q+1) t} $维列向量相乘的形式：

(9) $\begin{split} &{\displaystyle \sum _{k=1}^{t}{\displaystyle \sum _{q=0}^{Q}{a}_{q,k}}{\rm{tr}}\;\left({\boldsymbol{x}}_{t}{\boldsymbol{x}}_{t}^{\text{T}}{\boldsymbol{L}}^{q}\right)}={\displaystyle \sum _{k=1}^{t}{\boldsymbol{r}}_{t,t-k}^{\text{T}}}{\boldsymbol{v}}_{k}=\\ &\qquad\left[{\boldsymbol{r}}_{t,t-1}^{\text{T}},{\boldsymbol{r}}_{t,t-2}^{\text{T}},\cdots ,{\boldsymbol{r}}_{t,0}^{\text{T}}\right]{\left[{\boldsymbol{v}}_{1},\;{\boldsymbol{v}}_{2},\;\cdots ,\;{\boldsymbol{v}}_{t}\right]}^{\text{T}}=\\ &\qquad\left[{\boldsymbol{r}}_{t,t-1}^{\text{T}},{\boldsymbol{r}}_{t,t-2}^{\text{T}},\cdots ,{\boldsymbol{r}}_{t,0}^{\text{T}}\right]{{{\textit{φ}}}}_{t}. \end{split} $

令

(10) $ \left. \begin{aligned} {\boldsymbol{b}} =\;& {\left[ {\sum\limits_{t = 1}^{K - 1} {{\boldsymbol{r}}_{t,t - 1}^{\text{T}}} ,\sum\limits_{t = 2}^{K - 1} {{\boldsymbol{r}}_{t,t - 2}^{\text{T}}} , \cdots ,{\boldsymbol{r}}_{K - 1,0}^{\text{T}}} \right]^{\text{T}}}, \hfill \\ {{{{\textit{φ}}}}_{K - 1}} =\;& {\left[ {{a_{0,1}}, \cdots ,{a_{Q,1}}, \cdots ,{a_{0,K - 1}}, \cdots {a_{Q,K - 1}}} \right]^{\text{T}}}. \end{aligned}\right\} $

式中： $ {\boldsymbol{b}} $与 $ {{{{\textit{φ}}}}_{K - 1}} $均为 $ (Q+1)\left( {K - 1} \right) $维的向量. 则有，

(11) $ \begin{split} &\sum\limits_{t = 1}^{K - 1} {\left[ {\sum\limits_{k = 1}^t {\sum\limits_{q = 0}^Q {{a_{q,k}}} {\rm{tr}}\;\left( {{{\boldsymbol{x}}_t}{\boldsymbol{x}}_{t - k}^{\text{T}}{{\boldsymbol{L}}^q}} \right)} } \right]} = \hfill \\ & \qquad \sum\limits_{t = 1}^{K - 1} {\left[ {{\boldsymbol{r}}_{t,t - 1}^{\text{T}},{\boldsymbol{r}}_{t,t - 2}^{\text{T}}, \cdots ,{\boldsymbol{r}}_{t,0}^{\text{T}}} \right]{{{{\textit{φ}}}}_t}} = \hfill \\ &\qquad{\boldsymbol{r}}_{1,0}^{\text{T}}{{{{\textit{φ}}}}_1} + \left[ {{\boldsymbol{r}}_{2,1}^{\text{T}},{\boldsymbol{r}}_{2,0}^{\text{T}}} \right]{{{{\textit{φ}}}}_2} + \cdots + \hfill \\ &\qquad\left[ {{\boldsymbol{r}}_{K - 1,K - 2}^{\text{T}},{\boldsymbol{r}}_{K - 1,K - 3}^{\text{T}}, \cdots ,{\boldsymbol{r}}_{K - 1,0}^{\text{T}}} \right]{{{{\textit{φ}}}}_{K - 1}} = {{\boldsymbol{b}}^{\text{T}}}{{{{\textit{φ}}}}_{K - 1}}. \end{split}$

同理，令

(12) $ \left.\begin{aligned} &{{\boldsymbol{M}}_{t - {k_1},t - {k_2}}} = \hfill \\ &\left[ {\begin{array}{*{20}{c}} { {\rm{tr}}\;\left( {{{\boldsymbol{L}}^0}{{\boldsymbol{x}}_{t - {k_1}}}{\boldsymbol{x}}_{t - {k_2}}^{\text{T}}{{\boldsymbol{L}}^0}} \right)}&{ \cdots}&{{\rm{tr}}\;\left( {{{\boldsymbol{L}}^0}{{\boldsymbol{x}}_{t - {k_1}}}{\boldsymbol{x}}_{t - {k_2}}^{\text{T}}{{\boldsymbol{L}}^Q}} \right)}\\ { \vdots}&{}&{ \vdots}\\ {{\rm{tr}}\;\left( {{{\boldsymbol{L}}^Q}{{\boldsymbol{x}}_{t - {k_1}}}{\boldsymbol{x}}_{t - {k_2}}^{\text{T}}{{\boldsymbol{L}}^0}} \right)}&{ \cdots}&{{\rm{tr}}\;\left( {{{\boldsymbol{L}}^Q}{{\boldsymbol{x}}_{t - {k_1}}}{\boldsymbol{x}}_{t - {k_2}}^{\text{T}}{{\boldsymbol{L}}^Q}} \right)} \end{array}} \right], \\ &{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} {\displaystyle \sum\limits_{t = 1}^{K - 1} {{{\boldsymbol{M}}_{t - 1,t - 1}}} }&{\displaystyle \sum\limits_{t = 2}^{K - 1} {{{\boldsymbol{M}}_{t - 1,t - 2}}} }& \cdots &{{{\boldsymbol{M}}_{t - 1,0}}} \\ {\displaystyle \sum\limits_{t = 2}^{K - 1} {{{\boldsymbol{M}}_{t - 2,t - 1}}} }&{\displaystyle \sum\limits_{t = 2}^{K - 1} {{{\boldsymbol{M}}_{t - 2,t - 2}}} }& \cdots &{{{\boldsymbol{M}}_{t - 2,0}}} \\ \vdots & \vdots & {} & \vdots \\ {{{\boldsymbol{M}}_{0,t - 1}}}&{{{\boldsymbol{M}}_{0,t - 2}}}& \cdots &{{{\boldsymbol{M}}_{0,0}}} \end{array}} \right]. \end{aligned}\right\} $

式中： ${{\boldsymbol{M}}_{t- {k_1},t - {k_2}}} \in {{\bf{R}}^{(Q+1)\times (Q+1)}}$, ${\boldsymbol{C}} \in {{\bf{R}}^{\left( {(Q+1)\left( {K - 1} \right)} \right) \times \left( {(Q+1)\left( {K - 1} \right)} \right)}}$，则有

(13) $\begin{split} \sum\limits_{t = 1}^{K - 1} &{\left[ {\sum\limits_{{k_1} = 1}^t {\sum\limits_{{k_2} = 1}^t {\sum\limits_{{q_1} = 0}^Q {\sum\limits_{{q_2} = 0}^Q {{a_{{q_1},{k_1}}}{a_{{q_2},{k_2}}}} {\rm{tr}}\;\left( {{{\boldsymbol{L}}^{{q_1}}}{{\boldsymbol{x}}_{t - {k_1}}}{\boldsymbol{x}}_{t - {k_2}}^{\text{T}}{{\boldsymbol{L}}^{{q_2}}}} \right)} } } } \right]} = \hfill \\ & \sum\limits_{t = 1}^{K - 1} {{{{\textit{φ}}}}_t^{\text{T}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_{t - 1,t - 1}}}&{{{\boldsymbol{M}}_{t - 1,t - 2}}}& \cdots &{{{\boldsymbol{M}}_{t - 1,0}}} \\ {{{\boldsymbol{M}}_{t - 2,t - 1}}}&{{{\boldsymbol{M}}_{t - 2,t - 2}}}& \cdots &{{{\boldsymbol{M}}_{t - 2,0}}} \\ \vdots & \vdots & {} & \vdots \\ {{{\boldsymbol{M}}_{0,t - 1}}}&{{{\boldsymbol{M}}_{0,t - 2}}}& \cdots &{{{\boldsymbol{M}}_{0,0}}} \end{array}} \right]{{{{\textit{φ}}}}_t}}= \hfill \\ & {{{\textit{φ}}}}_{K - 1}^{\text{T}}{\boldsymbol{C}}{{{{\textit{φ}}}}_{K - 1}}. \\[-13pt] \end{split} $

此时，无约束优化问题式(6)可以转化为关于图滤波器系数向量 $ {{{{\textit{φ}}}}_{K - 1}} $的优化问题：

(14) $ {\mathop {\min }\limits_{{{{{\textit{φ}}}}_{K - 1}}} }\;\;\;\;{\sum\limits_{t = 1}^{K - 1} {{\rm{tr}}\;\left( {{{\boldsymbol{x}}_t}{\boldsymbol{x}}_t^{\text{T}}} \right)} - 2{{\boldsymbol{b}}^{\text{T}}}{{{{\textit{φ}}}}_{K - 1}} + {{{\textit{φ}}}}_{K - 1}^{\text{T}}{\boldsymbol{C}}{{{{\textit{φ}}}}_{K - 1}}} . $

由于其目标函数为二次凸函数，令目标函数的梯度为0，则有

(15) $ {\boldsymbol{C}}{{{{\textit{φ}}}}_{K - 1}} - {\boldsymbol{b}} = {\bf{0}}. $

式中： $ {\boldsymbol{C}} $为Hessian矩阵，优化问题的最优解为 $ {{{{\textit{φ}}}}_{K - 1}} = {{\boldsymbol{C}}^{ - 1}}{\boldsymbol{b}} $.

基于上述分析, 基于LS优化的GP-VAR模型对传染病传播的预测算法如下.

算法: 基于LS优化的GP-VAR模型对传染病传播的预测算法

输入: 预测所需天数 $ K $、图滤波器阶数 $ Q $、图节点数 $ N $、邻接矩阵 $ {\boldsymbol{W}} $、预测步数 $ n $、预测时刻 $ T $及其之前连续K d的历史病例数据 $ {{\boldsymbol{x}}_{T - K}},{{\boldsymbol{x}}_{T - K + 1}}, \cdots ,{{\boldsymbol{x}}_{T - 1}} $；

步骤:

1. 计算归一化拉普拉斯矩阵 $ {{\boldsymbol{L}}_{{\rm{nor}}}} $；

2. 令 $ {{\boldsymbol{X}}_1} = \left[ {{{\boldsymbol{x}}_{T - K}},{{\boldsymbol{x}}_{T - K + 1}}, \cdots ,{{\boldsymbol{x}}_{T - 1}}} \right] $， ${\boldsymbol{Y}} \in {{\bf{R}}^{N \times n}}$为全零矩阵；

3. for $m = 1:n$

基于式(8)、(10)、(12)，通过 $ {{\boldsymbol{X}}_m} $及 $ {{\boldsymbol{L}}_{{\rm{nor}}}} $得到 $ {{\boldsymbol{b}}^{\left( m \right)}} $、 $ {{\boldsymbol{C}}^{\left( m \right)}} $；

求解优化问题式(14)，得到滤波器系数向量 $ {{{\textit{φ}}}}_{K - 1}^{\left( m \right)} $, 包含 $ (Q+1) \left( {K - 1} \right) $个参数 $ a_{q,k}^{\left( m \right)} $;

生成误差向量 $ {{\boldsymbol{\varepsilon }}_{T + m - 1}} $；

基于式(4)，通过参数 $ a_{q,k}^{\left( m \right)} $、误差向量 $ {{\boldsymbol{\varepsilon }}_{T + m - 1}} $及病例数据 $ {{\boldsymbol{X}}_m}(:,2),{{\boldsymbol{X}}_m}(:,3), \cdots ,{{\boldsymbol{X}}_m}(:,K) $得到 $ T + m - 1 $时刻的预测值 $ {{\boldsymbol{\tilde x}}_{T + m - 1}} $；

令 ${{\boldsymbol{X}}_{m + 1}}= \left[ {{{\boldsymbol{X}}_m}(:,2),\;{{\boldsymbol{X}}_m}(:,3), \cdots ,{{\boldsymbol{X}}_m}(:,K),\;{{{\boldsymbol{\tilde x}}}_{T + m - 1}}} \right]$， $ {\boldsymbol{Y}}(:,m) = {{\boldsymbol{\tilde x}}_{T + m - 1}} $； // $ {\boldsymbol{X}}(:,j) $为矩阵 $ {\boldsymbol{X}} $的第 $ j $列， $ {\boldsymbol{Y}}(:,m) $为矩阵 $ {\boldsymbol{Y}} $的第 $ m $列.

end；

输出: $n$步预测结果 $ {\boldsymbol{Y}} $.

采用模拟疾病数据(SI数据、SEIRS数据)^[15]和德国COVID-19数据^[23]来验证所提出方法的性能. SI数据 ${{\boldsymbol{S}}_1} \in {{\bf{R}}^{{N_1} \times {M_1}}}$与SEIRS数据 ${{\boldsymbol{S}}_2} \in {{\bf{R}}^{{N_2} \times {M_2}}}$是基于飞行网络生成的模拟疾病数据, 包含 ${N_1} = {N_2} = 125$个欧洲主要城市连续 ${M_1} = {M_2} = $ $ 487$ d的感染人数. 利用城市的位置信息和机场航班信息对数据 $ {{\boldsymbol{S}}_1} $、 $ {{\boldsymbol{S}}_2} $构建一个节点数为125的城市连接图，图中的每个节点代表一个城市，该城市的感染疾病人数作为相应节点的图信号. 此外，城市之间的关联性由2个图进行刻画^[21]，一个是基于城市位置信息的6近邻图，边的权重为 ${{{W}}_1}\left( {i,j} \right) = {{{W}}_1}\left( {j,i} \right) = 1$，另一个是基于城市机场的航班信息的连接图，若2个机场间有航班，则相应节点间边的权重为 $ {{{W}}_2}\left( {i,j} \right) = {{{W}}_2}\left( {j,i} \right) = 0.01 $. 如图3所示为欧洲125个主要城市的连接图，其邻接矩阵为 $ {\boldsymbol{W}} = {{\boldsymbol{W}}_1} + {{\boldsymbol{W}}_2} $. 德国COVID-19数据 ${{\boldsymbol{S}}_3} \in {{\bf{R}}^{{N_3} \times {M_3}}}$包含了德国 $ {N_3} = 412 $个地区连续 $ {M_3} = 507 $ d的实际感染新冠病毒的人数，基于地区的位置信息构建一个节点个数为412的6近邻图，其中，图上节点表示地区，节点的图信号表示相应地区的感染疾病人数，且对于有边相连的节点，其边的权重均设为1. 如图4所示为德国412个地区的连接图.

将所提出的基于LS优化的GP-VAR模型与G-VARMA模型^[13]、原有的GP-VAR模型^[13]以及ARIMA模型^[6]作比较. 对不同类型的传染病传播数据进行训练和测试，首先，利用所提出的方法对不同类型的训练数据在不同参数下进行预测，并对比其预测结果，从而选出对不同类型数据都能有较好预测效果的参数值 $ K $、 $ Q $，然后根据得到的 $ K $、 $ Q $对测试数据进行预测^[13]. 对于所提出的基于LS优化的GP-VAR模型、G-VARMA模型和原有的GP-VAR模型，预测所需的连续K个时刻的历史数据、图滤波器的最高阶数 $ Q $均相同，即选定 $ K = 3 $、 $ Q = 5 $，以便于对比不同方法的性能. 此外，由于ARIMA模型是单变量模型，根据文献[6]中的方法对每个地区分别选择其合适的参数进行预测. ${{\boldsymbol{\tilde x}}_t}$作为预测值， ${{\boldsymbol{x}}_t}$作为真实值，采用归一化最小均方误差根rNMSE作为预测误差评价指标^[13]：

(16) $ {\text{rNMSE}} = \left[{{{\sum\limits_{t = 1}^n {\left\| {{{{\boldsymbol{\tilde x}}}_t} - {{\boldsymbol{x}}_t}} \right\|_2^2} } \mathord{\left/ {\vphantom {{\sum\limits_{t = 1}^n {\left\| {{{{\boldsymbol{\tilde x}}}_t} - {{\boldsymbol{x}}_t}} \right\|_2^2} } {\sum\limits_{t = 1}^n {\left\| {{{\boldsymbol{x}}_t}} \right\|_2^2} }}} \right. } {\sum\limits_{t = 1}^n {\left\| {{{\boldsymbol{x}}_t}} \right\|_2^2} }}}\right]^{1/2} . $

如图5所示为不同类型的传染病数据下连续50 d的一步预测结果，展示了本研究提出的方法在不同数据下连续50 d进行一步预测的结果. 图中，D为天数，B为感染人数. 一方面，可以看出SI数据随时间变化剧烈，这是由于SI模型中患者恢复情况的突然转变会导致数据的突变. 而COVID-19数据随时间变化较平稳，感染疾病的人数一直随时间逐渐增加. 另一方面，可以看出本研究提出的方法对COVID-19数据的预测结果最接近其真实值，而对随时间变化剧烈的SI数据的预测结果与其真实值的差异最大，说明数据本身的特性会对预测效果产生重要影响.

如图6所示为 $ K = 3 $、 $ Q = 5 $时,不同方法在SI、SEIRS、COVID-19数据下的6步预测误差对比. 图中，n表示第n步预测. 相比其他方法，本研究所提出的方法在对不同特性的数据进行预测时效果最好. 在对SEIRS数据和COVID-19数据进行短期预测时，本研究所提出的方法的预测误差并没有明显低于G-VARMA模型的，但长期预测误差低于原有的GP-VAR模型和G-VARMA模型的. 这很大程度上取决于本研究提出的方法在估计滤波器系数时，更好地利用了信号间的相关性及数据的内在关联性，提升了滤波器整体性能，从而降低了预测误差. 在对COVID-19数据进行预测时，ARIMA模型的预测误差逐步下降，但本研究提出方法的短期预测效果更好. 这是由于COVID-19数据随时间变化较平稳，ARIMA模型可以较好地利用数据的平稳性，这也是其对于随时间变化剧烈的SI数据的预测效果最差的原因. 此外，GP-VAR模型的滤波器数目少于G-VARMA模型，设计自由度较低，导致原有的GP-VAR模型在变化较平稳的数据的短期预测精度方面明显落后于G-VARMA模型. 这也说明了原有的GP-VAR模型、G-VARMA模型以及ARIMA模型不具有普适性，对数据类型要求较高，而本研究所提出的基于LS优化的GP-VAR模型具有普适性，对不同类型的数据均能保持良好的预测效果.

如图7所示为不同参数下的第5步预测误差对比图. 当 $ K = 4 $时，不同 $ Q $下本研究提出的方法对SEIRS数据的第5步预测误差如图7(a)所示；当 $ Q = 4 $时，不同 $ K $下本研究提出的方法对SEIRS数据的第5步预测误差如图7(b)所示. 通过调整 $ Q $、 $ K $，可以观察到rNMSE整体上分别随 $ Q $、 $ K $的增加而增大. 这说明本研究提出的方法可以在已知的之前时刻数据较少、图滤波器的阶数较低时，对传染病传播进行较好的预测. 此外， $ Q $、 $ K $较小会减少须估计的模型参数的数量，极大地降低优化问题的复杂度，便于优化问题的求解.

基于GP-VAR模型对传染病传播过程中感染疾病的人数进行预测，重新设计了一种基于时间-空间相关的LS优化问题用于模型参数的估计. 相较于现有预测方法，本研究所提出的方法进一步考虑了数据在时间维和空间维的关联性，从而更准确地刻画了传染病的传播特性. 仿真实验表明，本研究所提出的方法能够更好地刻画传染病传播过程中的时间相关性以及空间相关性，有更优的预测性能.

本研究采用基于图的预测方法，以图上的边及权重来反映传染病传播过程中的空间关联性，图模型不随时间变化. 但在现实生活中，地区间的人群流动情况是随时间变化的，传染病传播过程中的空间关联性也会随时间变化. 未来改进工作如下：一是基于时变图进行预测，以进一步刻画传染病的传播特性；二是将模型扩展到其他时空数据预测场景，如交通流量预测、温度预测.