地表水水质监测站点能够实时获取所处水域的水质状况. 基于获取的地表水水质数据并结合其他环境大数据可以实现高精度地表水水质预测，具有重大的实际应用意义. 例如监察人员能够提前感知流域的污染风险，向上下游站点发送预警信号，并回溯可能的污染来源. 对于预测值和实际值相差较大的点位，可以联动相关人员发起人工排查，提高线下执法效率.

实现高精度水质预测面临许多挑战. 水质数据在时间和空间2个维度上有复杂的依赖关系. 比如在时间维度上，站点某时刻的污染情况可能会与前一天的同一时刻相似；在空间上，站点的污染情况会受到周围水域的影响，传统的数据分析方法难以建模这种复杂的关联性.

许多研究者已经提出过一些方法解决地表水水质预测问题. 现有的预测方法可以粗略分为基于物理模型的方法和基于机器学习的方法2类.

基于物理模型的方法^[1-2]以领域内专业模型为基础，模拟水污染物从排放到扩散的全过程，从而预测各监测站点未来的水质情况. 该方法对污染源数据的完整性要求非常高，且模型参数需要专业人员根据不同的水域场景作调整适配，泛化能力较差.

基于机器学习的方法是一种数据驱动的方法，通过挖掘历史数据的隐含的复杂依赖关系来预测监测站点未来一段时间内的各项指标. 基于机器学习的方法，能够自动学习从多源数据输入到实际结果间的复杂映射关系. 这类方法可以进一步细分为时间序列分析模型、通用回归模型和深度学习模型. 时间序列分析模型^[3-4]通过发现历史序列中的线性规律来预测未来变化趋势，但无法接受序列数据以外的特征输入. 通用回归模型^[5-8]能够支持多源数据输入，如线性回归模型、支持向量机、随机森林等. 这类方法的缺点是为了达到可观的预测准确率，需要进行复杂的特征工程，提升了模型的落地应用难度.

深度学习模型^[9-11]通过堆叠多层神经网络来拟合从输入到输出的非线性复杂映射，可以实现自动化特征学习，解决传统回归模型的特征选择问题. 针对水质数据在时间和空间维度上的复杂依赖关系，本文提出基于图神经网络的地表水水质预测模型. 该模型根据站点在水网上的分布构建站点图，采用图神经网络建模地表水水质监测站点在空间上的复杂依赖关系，使用长短时记忆网络建模水质指标序列在时间上的复杂依赖关系，将时序编码结果输入到解码器中，得到最终的多步预测输出.

定义 $ {{X}}\in {\bf{R}}^{{\tau }_{\rm{in}}\times N\times D} $为模型输入，其中 $ {\tau }_{{\rm{in}}} $为输入时间窗口长度， $ N $为站点数量， $ D $为输入特征维度. 构建图 $ G=(V,E) $，用于建模站点间的复杂依赖关系. 其中 $ V $为站点集合， $ E $为图上边的集合. 水质预测的任务可以定义如下：建立地表水水质预测模型f(·)，利用历史水质指标监测数据和气象数据 $ {{X}} $，基于反映站点间依赖关系的图 $ G $，对未来一段时间内所有站点的水质指标 $ {{Y}}\in {\bf{R}}^{{\tau }_{\rm{out}}\times N\times d} $（ $ {\tau }_{\rm{out}} $为预测步长， $ d $为水质指标数量）进行预测：

(1) $ \left\{ {{{X}};G} \right\}\mathop \to \limits^{f\left( \cdot \right)} {{Y}}. $

模型的框架如图1所示. 预处理站点历史监测数据和气象等其他数据，根据站点在水网上的分布构建水网距离矩阵. 根据水网距离矩阵构建站点图，用于建模站点间的复杂依赖关系. 在每个时刻，将对应时刻的站点监测数据和气象数据作为节点输入，将站点间的距离权重作为边的输入. 利用图卷积神经网络（graph convolutional network，GCN）^[12]，实现节点间的信息交互. 采用基于LSTM^[13]的编码器，将站点信息编码为一个表示向量；将表示向量送入解码器，得到多步预测输出.

水质预测模型使用的数据输入包括地表水水质监测站点历史监测数据和气象数据.

1）站点监测数据. 站点监测数据的数据指标及其有效值范围如表1所示，包含5个反映地表水水质情况的监测指标，分别为PH值、氨氮质量浓度、总磷质量浓度、高锰酸盐指数和溶解氧质量浓度. 因为水质监测系统在运行的过程中会在中间链路上的环节出现问题，导致监测数据库中包含了一些异常值. 因此本文基于地表水水质国家标准^[14]，为不同的指标设定了不同的有效值范围（除PH值外，有效值上界为最差的V类水质标准上浮50%）以排除异常值. 所有不在上述有效值范围内的站点监测数据将被认为是异常值，在后续处理流程中视同缺失值.

2）气象数据. 气象因素对水污染物的产生与扩散有很大影响. 将气象因素引入模型计算中，能够基于深度学习技术，学习从气象因素到地表水水质的复杂影响关系.

气象数据的数据格式如表2所示，其中包含4个与水质相关性较大的气象指标，分别为气温、气压、湿度和降雨量. 本文为气象指标设定了对应的有效值范围（根据气象指标的实际含义）. 与站点监测数据类似，所有不在上述有效值范围内的站点监测数据将被认为是异常值，在后续处理流程中视同缺失值.

3）缺失值处理与归一化. 实际数据采集中出现的许多问题都会带来数据的缺失和异常，其中异常数据经过有效范围的判断后，可以视同缺失值，在实现时采用线性插值法实现缺失值填充. 线性插值法假设包含缺失值的一段连续序列满足线性变化关系，以线性函数去拟合数据变化，线性插值法的计算公式如下：

(2) $ p\left(x\right)=f\left({x}_{0}\right)+\frac{f\left({x}_{1}\right)-f\left({x}_{0}\right)}{{x}_{1}-{x}_{0}}(x-{x}_{0}). $

式中： $ {x}_{0} $和 $ {x}_{1} $为缺失值前、后的有效数据样本， $ p\left(x\right) $为拟合出的线性变化函数. 站点监测数据和气象数据都是在时间维度上平滑变化的，因此可以通过前、后有效观测值来填充.

模型输入数据的分布存在较大差异，在数据预处理中加入归一化操作，可以有效加速模型收敛，缓解潜在的梯度问题，归一化计算公式为

(3) $ {x_{{\rm{norm}}}} = \frac{{x - \mu }}{\sigma }. $

式中： $ \;\mu $为样本均值， $ \sigma $为样本标准差， $ {x}_{{\rm{norm}}} $为归一化操作后的输出数据.

1）站点图构建. 在输入时间窗口内的每个时刻，通过构建站点图 $ G=(V,E) $来建模站点间的复杂空间依赖关系. 在该场景中，图上的节点对应区域内的水质监测站点，边表示2个站点存在的交互关系. 构建图 $ G $的方式如图2所示. 图中， $ {{x}}_{a} $为 $ a $节点的信号， $ {{e}}_{a,b} $为连接 $ a $节点和 $ b $节点间边上的权重.

$ {{x}}_{a} $由对应时刻的水质指标监测数据和气象数据拼接组成， $ {{e}}_{a,b} $由站点地理位置计算得出，具体的计算公式如下：

(4) $ {{e}}_{a,b}=\left\{\begin{array}{l} {{{d}_{a,b}}^{-1}},\; 0 <{d}_{a,b} <{R}_{\rm{h}};\\ 0,\;\;\;\;\;\;\;\;\;{\text{其他}}.\end{array}\right. $

式中： $ {d}_{a,b} $为2个站点在水网上的距离， $ {R}_{\rm{h}} $为距离阈值. 在构图时，只计算距离小于 $ {R}_{\rm{h}} $的2个节点间边上的权重. 考虑到水污染物一般沿河道传播，采用水网距离来度量站点间的相关性，采用欧式距离度量会引入错误的先验知识.

2）图卷积计算. 在完成图的构建过程后，通过GCN实现节点间的信息交互，以建模站点间的空间依赖关系. 具体的计算公式为

(5) $ {{{H}}^{\left( {l + 1} \right)}} = \sigma \left( {{{{\tilde{ D}}}^{ - \frac{1}{2}}}{\tilde{ A}}{{{\tilde{ D}}}^{ - \frac{1}{2}}}{{{H}}^{\left( l \right)}}{{{W}}^{\left( l \right)}}} \right). $

式中： ${\tilde{ A}} = {{A}} + {{I}}$，即图的邻接矩阵加上self-loops，一般情况下邻接矩阵中的值只有1或者0以表示节点间的连接情况，该场景中的邻接矩阵计算时考虑了边权； $ {\tilde{ D}}^{-\frac{1}{2}} $用来对矩阵行列进行归一化处理； $ {{{H}}}^{\left(l\right)} $表示GCN第 $ l $层的节点信号， $ {{{H}}}^{\left(0\right)} $即对应时刻的节点输入信号 $ {{{X}}}^{t} $； $ {{{W}}}^{\left(l\right)} $表示第 $ l $层的变换参数，在训练过程中通过反向传播方法学习；σ为激活函数.

3）序列编解码. 在输入时间窗口内的所有图上完成图卷积操作后，每一个站点更新后的信号可以按时间顺序组成序列 $\{{{x}}^{1},{{x}}^{2},\cdots ,{{x}}^{{\tau }_{\rm{in}}}\}$. 序列将被送入编码LSTM中，得到每个站点的编码向量. LSTM是为了解决RNN网络梯度消失问题所提出的改进型序列建模神经网络，通过在模型内部引入多种功能的门机制，实现对序列长期依赖关系的建模. 具体的计算过程如下.

(6) $ {{{i}}^t} = {\rm{tanh}}\;\left( {{W_{\rm{i}}}{{{x}}^t} + {U_{\rm{i}}}{{{h}}^{t - 1}} + {{{b}}_{\rm{i}}}} \right), $

(7) $ {{{f}}^t} = {\rm{tanh}}\;\left( {{W_{\rm{f}}}{{{x}}^t} + {U_{\rm{f}}}{{{h}}^{t - 1}} + {{{b}}_{\rm{f}}}} \right), $

(8) $ {{{o}}^t} = {\rm{tanh}}\;\left( {{W_{\rm{o}}}{{{x}}^t} + {U_{\rm{o}}}{{{h}}^{t - 1}} + {{{b}}_{\rm{o}}}} \right), $

(9) $ {{{c}}^t} = {{{f}}^t} \circ {{{c}}^{t + 1}} + {{{i}}^t} \circ {\rm{tanh}}\;\left( {{W_{\rm{c}}}{{{x}}^t} + {U_{\rm{c}}}{{{h}}^{t - 1}} + {{{b}}_{\rm{c}}}} \right), $

(10) $ {{{h}}^t} = {{{o}}^t} \circ {\rm{tanh}}\;\left( {{{{c}}^t}} \right). $

式中： $ {W}_{\rm{i}} $、 $ {U}_{\rm{i}} $、 $ {W}_{\rm{f}} $、 $ {U}_{\rm{f}} $、 $ {W}_{\rm{o}} $、 $ {U}_{\rm{o}} $、 $ {W}_{\rm{c}} $和 $ {U}_{\rm{c}} $为可学习参数， $ {{b}}_{\rm{i}} $、 $ {{b}}_{\rm{f}} $、 $ {{b}}_{\rm{o}} $和 $ {{b}}_{\rm{c}} $为偏置向量， $ \circ $表示hadamard积.

编码LSTM得到的 $ {{h}}^{{\tau }_{\rm{in}}} $将作为站点的编码向量，用来初始化解码LSTM. 解码LSTM将每一步的隐状态 $ {{h}} $输入到多层感知机（multi-layer perception，MLP）中，以得到每一步的预测结果.

该任务的损失函数定义为均方误差（MSE），具体的计算方式如下：

(11) $ {\rm{Loss}} = \frac{1}{{N{\tau _{{\rm{out}}}}}}\mathop \sum \limits_{i = 1}^N \mathop \sum \limits_{j = 1}^{{\tau _{{\rm{out}}}}} \mathop \sum \limits_{k = 1}^d {\left| {{{Y}}_{i,k}^{t + j} - {\hat{ Y}}_{i,k}^{t + j}} \right|^2}. $

式中： $ {\tau }_{\rm{out}} $为预测步长， $ k $为地表水水质数据指标数， $ \hat{{{Y}}} $为地表水水质预测值.

采集杭嘉湖地区24个站点从2018年1月至2020年5月的历史水质指标监测数据和相应地区的气象数据. 地表水水质监测站点在地理上的分布如图3所示，地表水水质监测站点在地图上以圆点标注.

模型代码以Python语言编写，采用PyTorch框架^[15]构建神经网络模型，基于PyTorch Geometric库^[16]实现图神经网络计算. 模型的训练采用NVIDIA RTX2080Ti GPU.

模型的超参数设置如下：batch size设定为64，GCN的输出维度为16，LSTM的隐状态长度为32，GCN层数为1. 采用Adam优化器^[17]优化模型参数，学习率设定为0.002，训练epoch设定为150. 为了减少随机性，模型将独立训练5次，在实验部分取5次实验结果的平均值作为最终结果.

按照任务目标，预测任务定义为利用所有站点24个值（ $ {\tau }_{\rm{in}} $=24，采样周期为4 h）的历史水质指标监测数据及气象数据预测将来9个值（ $ {\tau }_{\rm{out}} $=9，采样周期为4 h）的水质情况. 将一个样本中包含长度为24的水质指标及气象数据作为模型输入，后面长度为9的水质指标作为地表水水质真实值. 用窗口大小为33的滑动窗口，以步长为1生成5 127个数据样本，按时间先后顺序以7∶1∶2的比例划分成训练集、评估集和测试集.

该实验以训练集的样本优化模型参数，以在评估集上的运行结果优化超参数，在测试集上验证该模型的性能.

性能评价指标包括平均绝对误差（MAE）和平均相对误差（MRE），计算方法如下：

(12) $ {{{\rm{MAE}}}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^N \left| {{{{Y}}_i} - {{{\hat{ Y}}}_i}} \right|, $

(13) $ {{{\rm{MRE}}}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^N \frac{{\left| {{{{Y}}_i} - {{{\hat{ Y}}}_i}} \right|}}{{{{{Y}}_i}}}. $

MAE直接计算真实值与预测值之间的绝对值. MRE在MAE的基础上除以真实值得到相对的误差率，结果更加直观，但不适用于数值范围在0附近的数据指标（相对误差会非常大）.

$ {R}_{\rm{h}} $直接影响了所构建图的边连接密度. 当 $ {R}_{\rm{h}} $很小时，有利于在图卷积计算时保留节点本地特征，但是不能充分地建模站点间的空间依赖关系. 当 $ {R}_{\rm{h}} $很大时，图上几乎所有的节点都有边的连接，在经图卷积计算后，图上所有节点的信号可能会趋同，丢失节点的本地特征，即产生过平滑现象. 模型损失函数MSE随 $ {R}_{\rm{h}} $的变化趋势如图4所示.

从实验结果可以看出，当 $ {R}_{{\rm{h}}} $=65 km时，MSE取得最小值0.375，优化效果最好. 在后续实验中， $ {R}_{{\rm{h}}} $设置为65 km.

除了历史水质指标监测数据外，模型还用到了气象数据和水网分布数据. 为了验证这2种数据的有效性，引入2个变体模型w/o wea和w/o water. 其中w/o wea去除了站点图节点信号的气象数据输入，w/o water以站点间的欧式空间距离代替水网距离来构建距离矩阵.

模型的性能对比结果如表3所示. 采用2个性能评价指标MAE和MRE. 5个水质指标分别为PH值、氨氮质量浓度、总磷质量浓度、高锰酸盐指数、溶解氧质量浓度. 其中因为氨氮质量浓度与总磷质量浓度的监测数据一直在0附近，MRE在数值上会趋近于无穷值，没有评价意义. 仅列出模型在PH值、高锰酸盐指数和溶解氧质量浓度3个指标上的实验结果. 表3中所有的结果为对应指标和模型所有步长结果的平均值.

从实验结果可以看出，该模型较2个变体模型具有性能上的优势，这说明模型中引入的气象数据和水网分布数据有助于实现水质指标预测. 水网分布数据比气象数据的影响更大，原因是基于空间欧式距离定义的相关性不仅不能充分表达站点间的依赖关系，还会引入错误的先验知识，造成模型无法正确地拟合真实的映射关系，而基于站点在水网上分布所构建的水网距离矩阵能够更加真实地反映站点间的相关性.

将提出的地表水水质预测模型与现有的基于机器学习方法的模型进行对比. 用于对比的基线模型如下.

1）支持向量回归机（SVR）^[5]. 该方法是传统支持向量机（SVM）的回归版本，目标是在特征空间内寻找最优的超平面，使得特征空间中的所有数据样本到该平面的距离最近. 原始版本的SVR只支持单变量输出，作为实验对比基线，本文为每个指标和预测步长建立SVR模型.

2）差分整合移动平均自回归模型（ARIMA）^[3]. 该方法将时间序列转化为平稳时间序列，将因变量对它的滞后值以及随机误差项进行回归拟合，是时间序列线性分析中最常见的模型之一.

3）长短期记忆网络（LSTM）^[10]. LSTM在序列建模过程中保持记忆性，较RNN与GRU有更加复杂的门控结构，设定遗忘门和输入门用以控制长期记忆的保留与更新，设定输出门用以得到每个时刻的输出结果. LSTM缓解了RNN模型中的梯度消失问题，是深度学习中建模时间序列最常用的神经网络.

为了评估各模型在短期预测及长期预测任务上的性能，列出了模型在不同预测步长下的误差. 如表4~6所示为提出模型与对比模型在PH监测值、高锰酸盐指数和溶解氧质量浓度3个指标上的预测结果. 表中，L为预测步长. 可见，与其他基线模型相比，该模型有明显的性能优势. 与时间序列分析方法ARIMA相比，该模型在PH监测值、高锰酸盐指数和溶解氧质量浓度3个指标上性能分别提升了23.3%、36.5%和37%. 与通用回归模型SVR相比，在3个指标上分别提升了26.6%、28.3%和36.8%. 与LSTM模型相比，在3个指标上分别提升了14.8%、20.0%和29.6%.

实验结果表明，与时间序列分析模型ARIMA和通用回归模型SVR相比，基于深度学习技术的LSTM和该模型在预测准确率上具有明显的优势. 这是因为在该场景中，输入数据的维度高达216维（24×9），开展特征工程的难度较大，而基于深度学习的模型支持多源异构数据输入，通过堆叠多层神经网络实现了自动化的特征学习. 该模型的预测性能比LSTM模型更好，这是因为该模型能够通过GCN建模站点在空间上的复杂依赖关系，例如水污染物在水网上的传播过程，相比深度单序列建模模型更完整地考虑了水质监测数据在时空两方面的复杂依赖关系.

本文提出新的地表水水质预测模型，使用图神经网络建模地表水水质监测站点在空间上的复杂依赖关系. 使用长短时记忆网络建模历史指标序列在时间上的复杂依赖关系，实现所有地表水水质站点水质指标的多步预测输出. 实验结果表明，与传统回归方法相比，利用该模型能够显著提高预测性能，与时间序列分析方法、通用回归方法和一般深度学习方法相比，该模型能够实现23.3%、26.6%和14.8%的性能提升.

本文提出的预测模型有进一步的拓展空间. 对于今后的工作，将重点探索以下3个方向. 1）引入更多影响地表水水质的因素，如站点周边污染源的实时排放信息、地表水文信息. 2）现有方法通过人为定义规则的方式实现图的构建，未来可以在图上实现节点间关系的自动推理，减少因预定义图带来的信息损失. 3）进一步扩大现有的建模范围，将区域从杭嘉湖局部区域扩展到整个浙江省水域或者是跨省域水网，以建模更加复杂的地表水水质监测站点间的空间依赖关系.