本发明属于水质异常检测领域，尤其涉及一种基于滑动窗多尺度主元分析的紫外光谱水质异常检测方法。
背景技术
水是生命之源，近年来，伴随着我国经济的高速发展，水体污染情况变得日益严峻。水污染事件会严重影响国民的生命健康安全及国家经济健康发展。水环境问题已经成为国民关注的焦点，开展水质监测技术的研究刻不容缓。
基于紫外光谱的水质异常检测方法，以其检测速度快、实时在线、无二次污染、可监测多种污染物的优点，成为了水质异常监测方法中的一个重要的研究方向。基于统计学的异常检测算法，以无监督的方式对水质进行监测，是紫外光谱水质异常检测中的一个重要组成部分，其中，又以主元分析方法提取水质指标特征值，结合统计量控制图的手段检测水质变化的方法最为常见。例如，kuhnert等人对光谱数据进行主元分析，运用t2统计量作为水质异常检测依据；hou等人运用概率主元分析和多变量监测图进行基于紫外光谱的水质异常检测。
传统的主元分析紫外水质异常检测方法仅限于模型固定的单尺度建模，需要进行过程没有序列相关性以及过程是线性的假设。但由于废水排放或天气等多种因素的影响，水质通常会呈现出一定的波动性及一定的时间相关性，同时进行在线光谱数据采集时，几乎所有实际过程的数据都是多尺度的，因此，需要对传统的主元分析紫外水质异常检测方法进行改进。研究表明，滑动窗口pca可针对传统pca模型固定的问题做出改进，mspca可实现多尺度的异常检测。
技术实现要素：
发明目的：针对现有技术存在的问题，本发明的内容在于改进传统主元分析紫外水质异常检测方法，使其可动态适应水质背景波动，并可对水质光谱数据进行多尺度的水质异常检测。
技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种基于滑动窗多尺度主元分析的紫外光谱水质异常检测方法，包括以下步骤：
(1)根据历史数据求取小波变换尺度l、窗口长度n、各尺度cusum控制限hl、小波重构数据cusum控制限h′；
(2)进行在线光谱数据采集，并等待数据填满窗口n；
(3)对窗口中的光谱数据进行基线校正及标准化预处理；
(4)对预处理后的光谱数据进行mspca计算，运用阈值法选取主元数目；
(5)应用cusum控制图进行各尺度上异常检测，若一尺度上的cusum统计量超出cusum控制限hl，即标记该尺度为异常尺度；
(6)组合异常小波尺度，并进行小波重构，对重构数据进行pca计算；
(7)计算重构数据的cusum统计量，若超出cusum控制限h′即标为水质异常，否则标记为正常，生成水质报告。
其中，步骤(1)中，运用白噪声检验的方法确定小波分解尺度l，包括以下步骤：
(1.1)对历史光谱数据进行i＝1层小波分解；
(1.2)对分解出的小波系数进行白噪声检验，若通过白噪声检验，进行阈值降噪处理，并进行i+1层分解，若不通过直接进入下一步；
(1.3)重复步骤(1.1)、(1.2)直至小波系数不能通过白噪声检验，放弃最后一次未通过检验的结果，最终确定小波分解尺度l。
其中，步骤(1)中，滑动窗口长度n的选择由小波分解尺度l决定：
n＝2l+1。
其中，步骤(1)中，各尺度cusum控制限hl，l＝1,2,…,l、小波重构数据cusum控制限h′由历史样本训练得到，其方法如下：
(1.4)选取正常水质样本，对光谱矩阵进行mspca计算，得到各尺度小波系数矩阵，并计算各尺度cusum统计量，剔除噪声点，选取各尺度cusum统计量最大边界量，作为对应cusum控制限，即hl，l＝1,2,…,l；
(1.5)选取重要小波系数组合重构，计算重构数据cusum统计量，剔除噪声点，选取重构数据cusum统计量最大边界量，作为对应cusum控制限，即h′；
(1.6)移动数据窗口，重复(1.4)-(1.6)，更新hl，h′数值，l＝1,2,…,l。
其中，步骤(4)中，对预处理后的光谱矩阵x进行滑动窗口mspca计算，计算方法如下：记光谱矩阵为x，x＝(x1′,x2′,…,xn′)t，xi′为基线校正及标准化处理后的光谱，选取正交小波算子w，对光谱矩阵x进行尺度l的分解；运用阈值消噪法对各个尺度的小波系数进行降噪处理，得到降噪后的粗尺度小波系数矩阵hl与细尺度小波系数矩阵gm，其中m＝1,2….l，对粗尺度小波系数矩阵hl、细尺度小波系数矩阵gm进行pca计算，采用阈值法动态确定主元数目。
其中，步骤(4)中，采用阈值法动态确定主元数目方法如下：计算各尺度上的小波系数矩阵的协方差矩阵c，对协方差矩阵c进行特征量分解，求得特征量λi，λ1≥λ2≥…≥λi，记特征量中非零个数为m，设定阈值为：
将特征量λi与阈值比较，选出所有比阈值大的特征量作为所提取的主元对应的特征值，选取出的特征量个数为主元数目k。
其中，步骤(5)中，引入平方预测误差qi，平方预测误差表征了小波系数矩阵g中小波行向量gi到残差空间的欧式距离，当某一尺度上小波发生异常时，平方预测误差会发生偏移；
小波行向量gi的平方预测误差计算公式为：
qi＝(git(1-pkpkt)gi)
pk为前k主元对应荷载向量pk组成的矩阵；
运用cusum控制图检测平方预测误差qi的偏移，cusum统计量ci的计算公式为：
ci＝max{0,qi-μ0+k+ci-1}
其中，ci-1为gi前一小波行向量gi-1的平方预测误差qi-1的cusum统计量，μ0为正常水样平方预测误差均值，k为检测过程中的允偏量；
当ci超出步骤(1)中计算出的cusum控制限hl时，即检测出平方预测误差qi发生偏移，控制图认为该尺度小波上发生了异常。
其中，步骤(6)中，异常小波系数进行重新组合重构，得到重构矩阵x′，重构矩阵计算公式如下：
其中
其中，步骤(7)中，对重构后的数据再次进行pca，根据重构数据的协方差矩阵重新计算cusum统计量，将cusum统计量与步骤(1)中计算出的cusum控制限h′作比较，若cusum统计量超出cusum控制限h′，则标记当前水质为异常水质，否则为正常水质，并生成最终的水质报告。
有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：
该方法改进了传统主元分析紫外水质异常检测方法，使其可动态适应水质背景波动，并可对水质光谱数据进行多尺度的水质异常检测，提高了针对动态变化水质检测的检测精度。
附图说明
图1为本发明框架结构图；
图2为本发明流程图；
图3为水质异常检测框架结构图。
具体实施方式
下面结合附图和实施例对本发明的技术方案作进一步的说明。
如图1及图2所示，本发明公开了一种基于滑动窗—多尺度主元分析的紫外光谱水质异常检测方法，包括以下步骤：
(1)根据历史数据求取小波变换尺度l、窗口长度n、各尺度cusum控制限hl，l＝1,2,…,l、小波重构数据cusum控制限h′等重要参数。其中，小波分解尺度l的选择与光谱图信噪比有关，运用白噪声检验的方式确定小波分解尺度，包括以下步骤：
(1.1)对原始光谱数据进行i＝1层小波分解；
(1.2)对分解出的小波系数进行白噪声检验，若通过白噪声检验，进行阈值降噪处理，并进行i+1层分解，若不通过直接进入下一步；
(1.3)重复步骤(1.1)、(1.2)直至小波系数不能通过白噪声检验，放弃最后一次未通过检验的结果，最终确定小波分解尺度l，并且，窗口长度n的选择由分解尺度l决定：n＝2l+1。
各尺度cusum控制限hl，l＝1,2,…,l、小波重构数据cusum控制限h′由历史样本训练得到，其方法为：
第一步：选取正常水质样本，对光谱矩阵进行mspca计算，得到各尺度小波系数矩阵，并计算各尺度cusum统计量，剔除噪声点，选取各尺度cusum统计量最大边界量，作为对应cusum控制限，即hl，l＝1,2,…,l；
第二步：选取重要小波系数组合重构，计算重构数据cusum统计量，剔除噪声点，选取重构数据cusum统计量最大边界量，作为对应cusum控制限，即h′。
第三步：移动数据窗口，重复第一步至第三步，更新hl，l＝1,2,…,l、h′数值。
(2)在线光谱数据采集，每条光谱采集m个波长点的吸光度，记为xi1,xi2…xim，m为每条光谱的波长点数，等待在线采集的光谱数据数量到达窗口长度n，n即为步骤(1)中求取的窗口长度。
(3)光谱数据基线校正及标准化预处理。应用广义whittaker平滑器对采集到的光谱数据xi进行基线校正预处理：
式中，xi为待分析的原始光谱，fi为基线校正后的光谱，，fi,j为fi中第j个波长点上的吸光度值，i＝1,…n，n为步骤(2)求取的窗口长度，j＝1,…,m，m为每条光谱的波长点数。引入参数μ，当xi,j>fi,j时，wj＝μ，当xi,j≤fi,j时，wj＝1-μ，μ＝0.01～0.1，ε＝102～109。
接着对基线校正后的光谱数据进行标准化处理，提高模型分析能力。
标准化处理后的光谱xi′计算公式为：
其中，fi(fi,1,fi,2,…,fi,m)为基线校正后的光谱，
其中，n为窗口中的光谱条数，k＝1,2….,m，m为每条光谱的波长点数。
(4)如图3所示，进行水质异常分析，首先对窗口中的光谱数据进行mspca计算，记光谱矩阵为x，x＝(x1′,x2′,…,xn′)t，xi′为标准化处理后的光谱，n为窗口尺度。选取正交小波算子w，对光谱矩阵x进行尺度l的分解，l为步骤(1)中求解的分解尺度。运用阈值消噪法对各个尺度的小波系数进行降噪处理，得到降噪后的粗尺度小波系数矩阵hl与细尺度小波系数矩阵gm，m＝1,2….l。应用阈值法动态确定主元数目，计算各尺度上的小波系数矩阵的协方差矩阵c，特征量λi(λ1≥λ2≥…≥λi)及对应的荷载向量pi，记特征量λi(λ1≥λ2≥…≥λi)中非零个数为m，设定阈值为：
其中，λi为小波系数矩阵的协方差矩阵特征量，λ1为第一大特征值，λ2为第二大特征值，m为特征值中非零个数。
将特征量λi与阈值比较，选出所有比阈值大的特征量作为所提取的主元对应的特征值，选取出的特征量个数为主元数目k。
接着应用cusum控制图进行各尺度上异常检测，引入平方预测误差qi，平方预测误差表征了小波系数矩阵g中小波行向量gi到残差空间的欧式距离，当某一尺度上小波发生异常时，平方预测误差会发生偏移。
小波行向量gi的平方预测误差计算公式为：
qi＝(git(1-pkpkt)gi)
pk为前k主元对应荷载向量pk组成的矩阵。
运用cusum控制图检测平方预测误差qi的偏移。cusum统计量ci的计算公式为：
ci＝max{0,qi-μ0+k+ci-1}
其中，ci-1为gi前一小波行向量gi-1的平方预测误差qi-1的cusum统计量，μ0为正常水样平方预测误差均值，k为检测过程中的允偏量。
当ci超出步骤(1)中历史数据训练出的cusum控制限hl时，控制图认为该尺度小波上发生了异常。
在对各尺度小波进行异常检测之后，可生成初步水质异常报告，但为降低虚警率，最后需要异常尺度小波系数进行重新组合，重构，应用cusum控制图对重构数据进行异常检测，生成最终的水质异常检测报告。
组合发生异常的小波系数，并进行小波重构，得到重构矩阵x′，重构矩阵计算公式为：
其中
对重构矩阵x′再次进行pca主元分解，运用阈值法确定主元数目，并计算重构矩阵x′中行向量平方预测误差，求取cusum统计量，将cusum统计量与步骤(1)中历史数据训练出的cusum控制限h′作比较，若cusum统计量超出cusum控制限h′，则标记当前水质为异常水质，否则为正常水质，并生成最终的水质报告。