设i时刻观测的光谱数据用X_{i, obs}, i=1, …, N表征, N为观测时长。 可以假设X_{i, obs}由两部分叠加组成, 纯净的光谱X_i, 以及很难估计的、 由噪声及基线漂移等为主引起的光谱差异部分 hi[6]。 如果令空间ε ^D表征光谱的差异性, 则可认为ε ^D是由{h_i}, i=1, …, N张成的。 根据正交投影定理, 校正后的光谱X_{i, corr}可以被式(1)表征, 其中, 矩阵P为空间ε ^D的一组基

Xi, corr=Xi, obsP⊥=Xi, obs[I-PT(PPT)-1P](1)

可以证明, 由{h_i}, i=1, …, N张成的空间ε ^D可通过如下方法估计。 首先, 随机在训练集中抽取一定数量的正常光谱, 由于可以保证检测开始时水质为正常, 故可采集相同数量的正常水质光谱, 依时序保存至测试集。 将两数据集中正常光谱数据合并成为一个数据表, 在数据表上做主成分分析(principal component analysis, PCA), 提取到的光谱主成分代表了上述的两数据集光谱的差异性, 故主成分分析得到的荷载矩阵 P˙可作为P的估计。 则由PCA求解过程中的约束条件可知, 式(1)可简化为式(2)

Xi, corr=Xi, obs(I-P˙TP˙)(2)

在实际应用中, 因为正常水质本身是具有一定波动差异的, 剔除的差异越多越容易改变光谱中与污染物相关的信息, 不可以将光谱所有的主成分去掉。 所以利用正交投影校正的方法, 需要确定一个自然数参数p, 只取主成分分析中贡献率大的前p维矩阵来代替上述方法中的$\hat{P}$。

偏最小二乘法(partial least squares, PLS)回归方法是化学计量学与光谱分析中的常用方法, 具体计算步骤和使用方法详见参考文献[7, 8], 这里仅说明PLS方法中的一个重要参数— — 成分选取个数f的确定。 在建模过程中, f选取过大容易过拟合, 选取过小泛化能力不足^[8], 本文通过K折交叉验证法选取f。 每次PLS迭代后, 可以获得一个新成分, 并更新了反映光谱特征的回归系数向量B。 记录不同回归系数向量B时测试所得RMSECV的值, 作RMSECV随f增加时的曲线。 利用F检验方法, 自动判断RMSECV首次下降不显著时, 取该处对应的f作为模型参数。

PLS方法本质上是一种基于特征变量的回归方法。 若将其中的因变量变为哑变量, 由0或1进行替代, 则可成为提取光谱信息并进行分类的方法, 通常称为PLS-DA^[9]。 PLS-DA在本文建模过程中0作为正常时刻, 1作为异常时刻, 利用上述迭代方法建模。 得到模型后, 根据训练数据ROC曲线确定最优阈值Threshold_best。 设PLS-DA模型用F_f(· )表示, 在i时刻输出的报警信号用y_i表示, 即y_i=F_f(X_{i, obs})。 通过y_i与Threshold_best对比, 得到最终的离群点序列Y_i。 对于i时刻输入到模型的光谱X_{i, obs}, i=1, …, N, 其离群点判别方法由式(3)描述

Yi=Outliers, yi=Ff(Xi, corr)> ThresholdbestNormal, yi=Ff(Xi, corr)≤Thresholdbest(3)

由于原始数据的噪声干扰和其他因素的影响, 在去掉基线漂移之后, 仍然会存在单个的离群点。 一般而言, 污染物入侵事件会持续多个检测周期, 故连续多个离群点有更大概率为水质异常。 利用序贯贝叶斯更新概率的方法, 可将连续多个水质异常点检出为事件。 在时序分析中, 当前点是否是异常点的概率是由历史的观测点决定的。 序贯贝叶斯方法得到的t+1时刻的状态由贝叶斯公式(4)决定

π1(t+1)=π1(t)P(yt+1|θ1)(1-π1(t))P(yt+1|θ0)+π1(t)P(yt+1|θ1)(4)

式(4)中, π ₁(t)为t时刻为异常的概率, 初始时刻可认为发生概率较小, 故π ₁(0)为一个很小的值, 实验中取10^-5; 通过式(4)得到一个概率来表征为异常的可能性, 传递给下一个时刻; 当概率π ₁(t+1)超过0.95时认为当前是异常点, 设定为报警, 低于0.95时认为是正常情况, 设定为不报警。 式(4)中θ ₁为异常, θ ₀为正常, P(y_t+1|θ )可根据历史数据或训练数据的最佳阈值点位置处的检出率、 误报率进行估计^[10]。

受试者工作特征曲线^[11](receiver operating characteristic curve, ROC curve)是以真阳性率(检出率)为纵坐标, 假阳性率(误报率)为横坐标绘制的, 直观判断检测模型性能的曲线^[12]。 在同一光谱数据集中测试, 利用ROC曲线下的面积(area under receiver operating characteristic, AUROC)对两种或两种以上检测模型进行性能比较, 一般0.5< AUROC< 1, 且AUROC越大, 则分类器效果越佳; 在同一种水质检测模型中, 一般认为检出率越高、 误报率越低, 检测结果越好, 故使用ROC曲线上最靠近左上角的点处所代表的阈值是一个平衡点, 可认为是该模型的最优阈值Threshold_best。

本文使用的性能指标包括如下三种: 查准率、 查全率及F-Score^[12]。 查准率(precision)是正确预报的异常光谱样本数占预测为异常的样本总数的比率; 查全率(recall)是正确预报的异常光谱样本数占样本中异常点总数的比率。 一般希望二分类算法的查准率和查全率都能接近1, 但查准率和查全率两个指标会出现此消彼长的矛盾情况, 因此需要综合考虑, 最常见的方法是使用F-Score参数, 即上述两参数的加权调和平均值, 如式(5)

Fα=(α2+1)Precision×Recallα2(Precision+Recall)(5)

当α =1时, 即为F₁参数。 F₁参数越接近1, 性能越好。

城市饮用水管网检测实验装置是由接入实际管网系统的测量主管道以及光谱传感器组成, 如图2所示。 实验装置中的开环调节流量为600 L· h^-1, 连接实际管网。 所配的污染物标准浓度溶液通过脉冲泵及橡胶管注入到静态混合管道, 注入的流速通过一个闭环反馈系统控制脉冲泵, 根据主管道流量进行微调, 保证加药流量稳定在0.2 L· min^-1(注入污染物的流量与主管道流量之比为1: 50)。 静态混合管道将主管道流量与污染物溶液进行充分混合, 以确保污染物入侵事件在管道横截面方向的浓度分布均匀。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 城市饮用水管网检测装置结构图Fig.2 Detection device for urban water distribution system

主管道外接下导管, 连接浸入式10 mm光程的UV/VIS光谱浸入式传感器(S: CAN analysis: : lyser)。 污染物从注入口到下导管的水力学时间大概是11.67 s。 光谱仪每隔30 s进行一次测量, 其中20 s为测量时间, 10 s为延迟时间。 光谱仪选取的波长范围为240~400 nm, 分辨率为2.5 nm。 光谱仪监测数据可通过USB接口导出分析。

本文选用苯酚(Phenol)作为注入试剂。 苯酚等芳香酚类物质是2012年国家饮用水健康标准中列出的重要污染物, 可溶于水, 不慎摄入一定量会出现急性中毒症状, 饮用含有苯酚的水会抑制中枢神经系统, 损害肝、 肾功。 它们被广泛使用在化工生产中, 近年来苯酚等芳香酚类物质在管网污染事浸中曾被检出。 因此, 对于管网水中出现的苯酚, 需做出快速、 及时、 有效的反应, 并积极采取措施减小事故危害。

本文数据采集历时2周, 包括正常和异常样本。 每次加药持续10 min模拟污染物入侵事件, 每个浓度的污染物入侵事件持续20个采样周期。 训练集的浓度由3个不同浓度梯度组成(50, 100和200 μ g· L^-1)。 为了进一步测试检测方法的效果, 测试集设计为5个浓度梯度(30, 40, 50, 100和200 μ g· L^-1)。 两个批次的水质数据在不同的时间段采集。

本文采用正交投影基线校正方法, 对采集到测试集光谱数据进行校正, 以减小用于建模的训练数据与新采集的测试数据之间的基线漂移。 用训练集随机选取100个正常水质的数据、 与依时序采集的前100个测试集光谱数据提取正交投影向量 P˙。 由于两个批次的数据在不同时间段采集, 受外部因素影响, 实验的训练集谱线较为分散, 而依时序取样的测试集谱线较为集中, 如图3(a)所示。 经正交投影基线方法校正之后的光谱如图3(b)所示(取正交投影校正的参数p=2), 可见减小了两个数据集之间的基线漂移。 图3(b)中校正后的曲线污染物峰的位置更加明显, 在270 nm特征峰处可见测试集、 训练集光谱数据分层, 校正方法对有效信号进行了放大。 将(b)图测试集的数据去均值, 可在时序上见5次事件, 如图3(c)所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 苯酚紫外可见光谱(特征波长270 nm) (a): 训练集和测试集原始数据; (b): 训练集和测试集校正后; (c): 测试集校正后去均值的时序展示Fig.3 UV/VIS Spectra of Phenol (Characteristic wavelength: 270 nm) (a): Original data of training and testing set; (b): After correction; (c): Corrected and normalized testing set

在正交投影校正中, p参数的选择是关键。 如果p选择过大, 那么被去掉的成分越多, 越可能将含有污染物信息的成分去掉; 如果p参数选择过小, 那么可能会导致基线校正不足。 本文选择p参数秉承两个原则, 其一是消除不同批次水质光谱的背景差异, 其二是尽可能多的保持原光谱的信息, 由于不同数据间的差异形式和种类不同, p值会有所不同。 针对本文实验, 一种经验的选取参数p的方法是: 令p从0开始自增, 观察被去掉p维主成分后的正常水质光谱主成分。 当去掉p维主成分后的两个数据集中正常水质光谱主成分刚好接近重合时, 选择此时的p。 图4显示了p从0自增到3时, 两个数据集随机挑选的正常水质光谱的前两维主成分。 可见p=2时已经重合, 故选取的最佳p值不应大于2。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同参数p下的正常光谱主成分PC1, PC2 (p从0变化到3)Fig.4 First and second principle component among different p in normal spectrum(p Ranged from 0 to 3)

2.2.1 无监督异常检测方法

在无监督异常检测的研究中, 利用PCA提取新采集的光谱数据主成分, 通过统计量查看报警信号是否超过阈值, 超过阈值则被判为离群点。 图5所示为基于无监督的水质异常检测效果图。 图5(a)所示事件强度分别为30, 40, 50, 100和200 μ g· L^-1的苯酚污染事件。 利用PCA提取每条谱线的主成分, 根据主成分分析原理, 对按列服从高斯分布的原始矩阵进行主成分变换, 则可以根据Chi-square分布构建统计量^[4]进行统计检验, 若统计量超过阈值, 则认为当前监测数据为离群点。 阈值T_α 需要根据历史数据及经验确定。 图5(b)为计算的统计量指标作为报警信号, 统计量指标通过与阈值T_α 对比, 得到离群点序列如图5(c)。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 非监督方法检测效果 (a): 入侵事件实际浓度; (b): 统计量作为报警信号, 阈值Tα =0.65; (c): 离群点; (d): 异常概率; (e): 水质报警序列; (f): 与真实情况对比Fig.5 Detecting results of unsupervised method (a): Concentrations of intrusion events; (b): Statistic as alarm signal Tα =0.65; (c): Outliers; (d): Anomaly probability; (e): Alarm sequences; (f): Compared results with ground truth

由于背景水质波动、 工况变化或噪声干扰等情况也会产生离群点, 故对离群点序列采用了序贯贝叶斯方法, 进行异常概率更新, 得出该离群点是否属于水质报警序列的概率。 选取滑动窗长为200, 考虑当前时刻以前200个点的历史数据, 在滑动窗内根据序贯贝叶斯公式, 可计算得到当前时刻的异常概率为图5(d)所示, 其中的检出率、 误报率参数, 是根据历史数据通过经验得到的。 根据小概率事件原理, 设置概率阈值为0.95, 当超过该阈值时认为污染入侵事件一定发生, 将结果记录到水质预警序列, 如图5(e)所示。 图5(f)为报警序列图5(e)与实际情况图5(a)的对比, 得到的误报漏报结果。

综上可见, 在当前T_α =0.65情况下, 前两个浓度较低的30和40 μ g· L^-1事件漏报。 50 μ g· L^-1以上浓度的污染, 在水质基线稳定的情况下能够检出。 如果减小阈值T_α , 会增加误报FP, 增加阈值T_α , 会增加漏报FN, 很难使用现有方法在低浓度事件检出时获得一个高检出低误报的效果。 表1统计了在该数据集下的检测效果指标。

表1

Table 1

表1(Table 1)

表1 无监督方法检出效果 Table 1 Detecting results based on unsupervised method 检测效果指标 Precision Recall F1 非监督方法 0.45 0.37 0.40

表1 无监督方法检出效果 Table 1 Detecting results based on unsupervised method

2.2.2 监督学习异常检测模型建立

基于2.2.1提出的异常检测研究思路, 针对训练集的先验信息增加有监督学习环节。 利用本文提出的基于PLS-DA监督学习异常检测方法, 按图1中的建模流程, 首先利用全部训练集数据进行建模。 建模基于K=5的交叉验证, 通过多次迭代, 获取不同f下的RMSECV。 如图6所示, 为迭代25次RMSECV随f的变化, 由F检验得出f变化显著的位置为f=3。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 RMSECV变化(f从1到25)Fig.6 RMSECV curve changes (parameter f ranged from 1 to 25)

选取f=3, 可得PLS-DA的回归系数向量B, 向量长度同光谱的波数点一致, 如图7所示。 通过观察可见, 在苯酚特征波长270 nm附近出现了较大的数值。 由于使用PLS-DA方法在提取自变量矩阵和因变量矩阵主成分的同时, 能够使提取到的自变量主成分对因变量主成分有很强的解释能力, 所以在光谱数据分层的波长上回归系数的数值会偏大。 由此可见, 选用PLS-DA有监督方法对光谱进行特征提取, 具有一定的物理解释。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 PLS-DA回归系数(标出苯酚特征波长270 nm)Fig.7 PLS-DA coefficient of regression (labelled the characteristic wavelength 270 nm)

获取系数B之后, 可以得到训练集的测试结果。 在表2列了训练集的AUROC结果。 可见结果接近于1, 说明该模型对自身的解释能力较强。 通过训练集的ROC结果, 可以选取最接近左上角的阈值, 作为模型的最优阈值Threshold_best。

表2

Table 2

表2(Table 2)

表2 PLS-DA建模参数 Table 2 PLS-DA modelling parameters 污染物 f RMSECV 训练集AUROC Thresholdbest 苯酚 3 0.137 7 0.985 3 0.216 8

表2 PLS-DA建模参数 Table 2 PLS-DA modelling parameters

2.2.3 监督学习异常检测方法

针对本文提出的方法, 按图1中的检测流程, 将与2.2.1中同样的数据集作为观测所得的测试集, 利用正交投影校正方法进行基线校正(取p=2), 输入到2.2.2得到的PLS-DA模型中, 输出报警信号, 如图8(b)所示, 虚线部分为Threshold_best, 通过Threshold_best及式(3)判断报警信号。 得到的水质离群点如图8(c)所示。 同样利用了2.2.1中用于滚动更新异常概率的序贯贝叶斯方法, 取滑动窗长为200, 得到每个时刻离群点为异常的概率和最终的水质报警序列, 分别如图8(d)和图8(e)。 对检测效果指标进行统计, 结果如表3。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 监督学习方法检测效果 (a): 入侵事件浓度; (b): 报警信号, 阈值为0.216 8; (c): 离群点; (d): 异常概率; (e): 水质报警序列; (f): 与真实情况对比Fig.8 Detecting results of unsupervised method (a): Concentrations of intrusion events; (b): Alarm signal when threshold=0.216 8; (c): Outliers; (d): Anomaly probability; (e): Alarm sequences; (f): Compared results with ground truth

表3

Table 3

表3(Table 3)

表3 有监督方法检出效果 Table 3 Detecting results based on supervised method 检测效果指标 Precision Recall F1 有监督方法 0.78 0.76 0.77

表3 有监督方法检出效果 Table 3 Detecting results based on supervised method

根据测试数据集上的结果, 得到ROC曲线如图9。 得到基于PLS-DA的监督学习方法和基于PCA的非监督方法, AUROC分别为0.95和0.76。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 有监督PLS-DA方法, 和无监督的方法, ROC比较及AUROCFig.9 ROC curve and value of AUROC from PLS-DA method and unsupervised method

由表1、 表3及图5、 图8、 图9所得结果比较可知, 有监督方法不仅可有效检出浓度较低的30和40 μ g· L^-1事件, 也可以检出无监督方法有效的高浓度事件。 在当前含有低浓度情况的数据集下, 有监督方法的查准率、 查全率、 F₁以及AUROC均优于无监督方法。 由此可见, 有监督的PLS-DA方法能够在具有一定物理意义的情况下对预处理后的特征污染物光谱进行针对性的特征提取, 能够利用特征污染物的先验信息降低检出限。 在已知待测污染物类型的情况下, 有监督方法可作为非监督的方法的补充, 提升了在低浓度情况下的水质异常事件检测效果。

2.2.4 关于正交投影校正预处理对检出结果的影响

考虑实际使用情况中, 训练数据和测试数据之间的基线漂移现象, 2.2.3中采用正交投影方法校正基线。 为了验证正交投影方法可以有效校正基线, 提高监督学习方法对水质基线变化的适应性, 在训练集(左侧蓝色)和测试集(右侧红色)上的未经基线校正处理(p=0)的报警信号展示如图10(a)。 由于水质基线的漂移, 可以看出未经基线校正得到的Threshold_best测试集报警信号, 会发生一定程度的漂移, 从而使训练集中得到的在测试集中失效, 当前算法中无法有效检出测试集的低浓度事件。 作为对比, 图10(b)补充了图8(b)中未显示的训练集报警信号, 可见经过基线校正后的测试集报警曲线可用最优阈值判定低浓度情况下的异常。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 基线校正对训练集(左)和测试集(右)报警信号的影响 (a): 未基线校正处理的报警信号, p=0; (b): 基线校正后的报警信号, p=2Fig.10 Baseline correction impact on training set (left) and testing set (right) (a): Alarm signal without baseline correction, p=0; (b): Alarm signal with baseline correction, p=2