2. 绍兴市自来水有限公司, 绍兴 312000;
3. 浙江和达科技股份有限公司, 嘉兴 314006
2. Shaoxing Tap-water Co., Ltd, Shaoxing 312000, China;
3. Zhejiang HEDA Technology Co., Ltd, Jiaxing 314006, China
爆管往往伴随着短时的大水量漏失,是一个困扰供水行业的典型问题,不仅会造成水资源的浪费,同时也会影响正常供水,城市中发生的爆管事故甚至会引发道路积水、交通中断以及地面塌陷等次生灾害。因此,及时发现供水管网中的爆管事故,对于保障供水安全以及公共安全具有重要意义。
近年来,基于数据挖掘的爆管识别方法开始得到关注与研究。Mounce等[1-2]利用人工神经网络,建立基于规则的神经网络系统,分析观测值和预测值之间的偏差以判断管网数据的异常。同时,还将模糊推理系统结合神经网络应用于管网异常检测。Romano等[3]协同使用人工神经网络、统计过程控制技术和Bayes推理系统方法,进行爆管事故的检测。Ye和Fenner[4]建立自适应Kalman滤波模型,并通过估计值与观测值的残差与爆管的相关性进行事故判断,与神经网络方法相比,这种方法不需要大量的训练数据。此外,Ye和Fenner[5]还提出了一种基于加权最小二乘的方法,辅以期望最大化算法,这种检测方法可自动剔除历史异常数据,提高了模型的准确性。上述方法均在真实管网的计量分区(district metering area, DMA)中进行了验证,表明基于日常水力监测数据的爆管识别技术具有较强的爆管识别能力。但上述文献中的DMA大多是指互不连通、结构简单的独立计量区域[6],数据分析方法也仅仅基于单个流量计的数据。而国内的DMA情况有所不同,特别是由极其复杂的环状管网改造而成的DMA,可能存在多个出口和入口,与临近的DMA相互连通,这种多出入口的DMA在相关文献中少有研究。另外,前述方法均需建立预测模型,不仅过程复杂,预测的精度也会影响识别效果。
本文以Wu等[7]的研究为基础,将爆管识别分为3步:1) 通过数据转换重构DMA出入口流量数据并减小数据的波动;2) 使用聚类算法检测所有异常流量数据;3) 在所有异常数据中识别爆管事件。最后,将该方法应用到2个真实的多出入口DMA中,验证其有效性并初步探讨该方法的适用条件。
1 研究方法 1.1 研究地点与监测数据选取中国南方某市的2个DMA,并收集相应的流量数据进行分析。DMA1的面积约为6.5 km2,共有2个入口和3个出口;DMA2的面积约为DMA1的2倍,共有6个入口和8个出口。其中,2个DMA的出入口处都有相应的流量计进行数据采集,由于DMA2的一个出口流量计长期故障,将使用其他13个流量计的数据进行分析。
DMA1包含2015年5月31日到2015年12月16日共计200天的流量数据,并分别于11月17日和12月16日在4个不同排放口进行了7次爆管模拟实验,具体情况如表 1所示。DMA2包含2015年5月31日到2015年11月16日共计170天的流量数据,没有进行爆管模拟实验。其中,所有数据的采集间隔为5 min。
排放口 | 日期 | 实验时间 | 爆管流量/(m3·h-1) | 爆管流量占入流总量比例/% |
P1 | 2:17—2:27 | 约200 | 约13.3 | |
P2 | 2015年11月17日 | 2:30—2:40 | 约300 | 约17.2 |
P3 | 2:46—2:57 | 约400 | 约23.1 | |
P4 | 2015年12月16日 | 1:03—1:18 | 约200 | 约11.8 |
1:23—1:33 | 约100 | 约6.3 | ||
1:33—1:43 | 约150 | 约9.4 | ||
1:43—1:53 | 约200 | 约12.6 |
1.2 数据转换与数据标准化
对原始的时间序列按照时刻进行切分,构成切分序列,每个序列都代表某一时刻的流量数据随天数的变化[5, 7]。若数据采集间隔为5 min,每个流量计的数据都可以转化为288段(0:00,0:05,…,23:55) 切分序列。为了利用同一DMA所有流量计的数据关联性(指管网中发生事件时,DMA出入口流量计会随之有所响应,在数值上出现增减现象),将来自不同流量计相同时刻的切分序列进行组合,构成288个矩阵。矩阵形式如下:
$ \begin{matrix} \text{da}{{\text{y}}_{1}} \\ \vdots \\ \text{da}{{\text{y}}_{p}} \\ \end{matrix}\left( \begin{matrix} {{q}_{11}}&\cdots &{{q}_{1n}} \\ \vdots &{}&\vdots \\ {{q}_{p1}}&\cdots &{{q}_{pn}} \\ \end{matrix} \right). $ |
其中: q代表流量计读数;p为监测天数;n为DMA中流量计的个数,矩阵中的行向量代表每一天所有流量计在某一时刻的数值。每天新采集的流量数据将以行向量的形式传输进入相应时刻的矩阵,进行爆管事故的识别,本文将这些数据命名为测试向量。
各流量计监测得到的数据值域不同(即矩阵中的每列数据具有不同值域),为了避免数值较大的流量数据影响后续的计算结果,需要进行标准化处理,计算方法如下:
$ {x}'=\frac{{{x}_{i}}-\mu }{\sum\limits_{i=1}^{p}{\left| {{x}_{i}}-\mu \right|}}. $ | (1) |
其中:xi表示矩阵中某列数据的第i个值, p为监测天数(矩阵的行数), μ为该列流量数据的中值, 分母为该列数据的绝对标准差。
1.3 基于聚类算法的爆管识别方法 1.3.1 聚类算法在没有发生异常情况时,同一矩阵中的向量具有较高的相似度; 当出现异常情况时,相应的向量会与其他向量产生较大差异。基于此,本文选择了广泛用于异常值检测的聚类算法[8]来评估矩阵中向量间的相似度,从而发现异常的流量数据。
在对某个矩阵进行聚类分析时,需要计算每个向量的局部密度ρ和距离δ[9]。向量i的局部密度ρi的计算方法如下:
$ {{\rho }_{i}}=\sum\limits_{j}{\chi ({{d}_{ij}}-{{d}_{\text{c}}})}. $ | (2) |
其中:如果x < 0,则χ(x)=1,否则χ(x)=0;dij表示向量i与向量j间使用标准化数据计算的欧氏距离;dc为截断距离,dc使得向量的平均局部密度值为矩阵中所有向量数的1%~2%[9]。向量i的距离δi的计算方法如下:
$ {{\delta }_{i}}=\underset{j:{{\rho }_{j}}>{{\rho }_{i}}}{\mathop{\text{min}}}\, {{d}_{ij}}. $ | (3) |
其中,如果向量i具有最大的局部密度,其距离定义为
$ {{\delta }_{i}}=\text{ma}{{\text{x}}_{j}}{{d}_{ij}}. $ | (4) |
根据Wu等[7]的研究,本文对具有最低密度的向量按其距离由大到小排序,距离越大,向量为异常的可能性越大。同时,使用相应的参数来量化异常向量的距离,如果参数取值为0.1,则排在前10%的向量为异常向量。
1.3.2 异常向量的检测1) 初始化过程。
历史监测数据中可能含有因为各类事件导致的异常向量,如果历史异常向量与测试向量有较高相似度,该向量就不能被成功检测为异常。因此,应先使用聚类算法去除所有历史异常向量。在此过程中,量化异常距离的参数为初始化参数,它决定了一个历史数据矩阵中异常向量的个数。本文对2个DMA各使用其前150天(5月31日—10月27日)的数据建立历史数据矩阵并进行初始化,最终生成288个初始化矩阵,剩余数据用于生成测试向量。
2) 检测过程。
异常向量的检测是一个反复使用聚类算法并对矩阵进行更新的过程。完成初始化后,测试向量会依次进入相应时刻的矩阵并进行聚类分析,此阶段量化异常距离的参数为检测参数,与初始化参数略有不同,它仅决定了测试向量是否异常。一旦被识别为异常,测试向量则被剔除出矩阵,否则将继续保留在矩阵中。当第2个测试向量进入相应矩阵时,进行相同操作直到所有测试向量被检测完毕。在真实管网中,异常流量的出现频率较低,因此新向量会不断进入矩阵并保留下来,这将造成数据量增多并降低聚类效率。为了避免上述情况的发生,每对矩阵进行一次聚类分析后,都会将矩阵中的第一个向量删除。例如,对0:00时刻矩阵进行第一次聚类后,5月31日对应的向量就会被删除。所有检测出的异常向量将进入下一阶段,识别是否有爆管事故发生。
1.3.3 爆管事故的识别根据Wu等[7]对DMA流量波动特征的分析与总结,可以依据入口流量与出口流量的数值变化特征对爆管事故与其他异常情况进行区分。记异常的测试向量为x=(x1, x2, …, xn),并计算其相应矩阵的中值向量m=(μ1, μ2, …, μn),此中值向量的各个元素都由矩阵中相应列的中值构成。由于切分后的流量数据并不符合正态分布[10],相较于Wu等[7]采用的均值向量,中值向量更具鲁棒性[8, 10]。分析x中的各个元素,以其中较大元素(相较于中值而言)的个数为标准,将其分为3类,如表 2所述。爆管识别的逻辑判断过程如图 1所示。
类别 | 原因 | 流量变化特征 | 较大元素的数量 | 数学表达式 | |
入口流量 | 出口流量 | ||||
突增 | 晴热天气、节假日 | 数值增加 | 数值增加 | n | x-m>0 |
突减 | 突然降温 | 数值减小 | 数值减小 | 0 | m-x>0 |
爆管 | 爆管事故、 区域内大用户用水、 区域内管道冲洗 | 数值增加 | 变化不明显或出现 数值减小的情况 | (0, n) | 其他 |
注:n代表流量计的个数, x为异常的测试向量, m为相应的中值向量。 |
1.4 评价指标
本文的爆管识别方法本质上是一种二元分类器,即把所有测试向量分为爆管和非爆管两类。因此,该方法的评估指标可以采用真阳性率(true positive rate, TPR)和假阳性率(false positive rate, FPR),FPR亦可以称为误报率,计算公式如下所示:
$ \text{TPR=}\frac{\text{TP}}{\text{TP+FN}}, $ | (5) |
$ \text{FPR=}\frac{\text{FP}}{\text{FP+TN}}\text{.} $ | (6) |
其中: TP代表真阳性,即正确被识别为爆管的次数;FP代表假阳性,即本应该是非爆管却被识别为爆管的次数;TN为真阴性,即被正确识别为非爆管的次数;FN为假阴性,即本应是爆管却被识别为非爆管的次数。TP和FP出现时,均会生成报警。
管网中爆管发生频率较低,因此一定要控制报警次数,即在保证TPR较高的前提下,误报率一定要处于较低水平[11]。否则,过多的报警次数会加重供水企业工作人员的工作任务,反而带来不利影响。
2 结果与讨论 2.1 聚类过程为了阐述聚类算法的有效性,本文以DMA1的数据为例,将2015年11月17日得到的模拟爆管数据随机替换到历史数据中(本文中替换了7月27日的数据,即在7月27日应检测出爆管事故)。图 2显示了对2:50时刻相应矩阵进行聚类分析后的结果。图 2a是进行聚类后得到的决策图,描述了所有向量的局部密度和距离,图中最左侧的点都具有最低的密度,而方形点与三角形点具有较大距离,因此相应的向量被识别为异常。图 2b描绘了DMA1的5支流量计150 d中在2:50时刻的流量变化曲线,其中方形点代表的异常向量是由于10月9日区域下游的管道冲洗造成的,可以发现,每支流量计的数值都明显增大;而三角形点代表的异常向量是由于7月27日的爆管造成的,仅有入口流量计的数值增大,而出口流量计的数值变化不明显甚至减小。通过上述分析可知,聚类算法能够准确识别异常的流量波动。
2.2 参数的取值与识别效果
DMA1剩余的50天(10月28日—12月16日)的流量数据可生成14 400(50×288) 个测试向量,在第21天和第50天共进行了7次爆管模拟实验,根据这些实验的持续时间,共应检测出16个表征爆管事故的测试向量。
对所有数据进行聚类分析时,共包括初始化参数和检测参数2个参数,这2个参数的值共同决定了该爆管识别方法的TPR和FPR。采用枚举法,以0.01为步长,在0.01~0.15的范围内搜索2个参数的最优值,即在所得225组结果中寻找最优结果(高TPR与低FPR)。最终,初始化参数和检测参数的值分别被确定为0.13和0.10。如表 3的识别结果所示,在最优参数下,该方法共识别出13个表征爆管的向量,TPR和FPR分别为81.25%和1.51%。Wu等[7]使用相同的历史数据进行初始化并对21天的数据进行测试时,FPR仅为0.60%,这是因为本文加入了爆管流量更小的实验,为了能够成功识别,2个参数的值有所增加。这使得该方法对爆管事故更加灵敏,但误报率也相应升高。
日期 | 爆管模拟实验 | 爆管向量 对应时间 | 识别结果 |
2015年 11月17日 | 实验1 | 2:20 | FN |
2:25 | TP | ||
实验2 | 2:30 | FN | |
2:35 | TP | ||
2:40 | TP | ||
实验3 | 2:50 | TP | |
2:55 | TP | ||
2015年 12月16日 | 实验4 | 1:05 | TP |
1:10 | TP | ||
1:15 | TP | ||
实验5 | 1:25 | TP | |
1:30 | FN | ||
实验6 | 1:35 | TP | |
1:40 | TP | ||
实验7 | 1:45 | TP | |
1:50 | TP |
2.3 异常向量与天气的关系
本文使用最优参数的情况下,共检测出了1 673个异常向量,除了被识别出的13次真实爆管,有217次误报,根据表 1所述,这些极有可能由DMA内大用户用水造成。其他两类异常向量中,突增向量共862个,突减向量共581个,它们并非均匀分布在50 d中,其出现的数量与天气的变化有很强的关联性。如图 3所示,在11月5日天气由雨转晴,突增向量开始增多,随着温度的大幅上升,一直到11月7日,突增向量大量出现;而在11月8日,温度突然下降,被识别为突减的异常向量显著增加,突增向量随之减少。
2.4 DMA出入口数量与位置的影响
DMA2剩余的20 d的流量数据可生成5 760个测试向量,由于没有进行爆管实验,默认所有测试向量都为正常向量。在本研究中,设置初始化参数和检测参数分别为0.10和0.05,得到FPR高达10.80%,尽管2个参数的值小于在DMA1中使用的参数值,FPR却高出数倍。在DMA2中,没有TPR来验证该方法对爆管的敏感性,但是仅从FPR来看,其应用效果并不理想。经分析,DMA2的面积大、出入口多,而增设某些出入口的目的是防止水质恶化,其位置并不在供水干管上,流量小且波动大;另外,有些出入口可能靠近大用户,流量很容易因为大用户非常规的用水行为产生较大波动。上述情况都会导致相应的流量数据缺乏较强的规律性,图 4展示了该区域某只流量计连续3天流量数值的变化情况,其中负数表示管道中的水反向流动。图中信息表明,该流量计的数据变化无规律,即使进行数据转换也无法减小其波动性。当将此流量计的数据删除,仅使用其他12支数据规律较强的流量计时,在相同的参数值(0.10和0.05) 下,FPR降低为4.43%,但误报次数依然相对较多。因此,DMA出入口数量及其位置会显著影响该方法的识别效果。
3 结论
本文基于流量监测数据,应用聚类算法识别了多出入口DMA的爆管事故,研究结论如下:
1) 聚类算法避免了建立复杂预测模型的过程,依靠对流量数据进行相似度的比较,即可发现各类异常数据。
2) 多类事件均可引起DMA流量数据的异常波动,但是对出入口流量数据的影响各不相同。因此,利用所有出入口流量数据的关联性进行分析能够将异常数据分类,进一步降低爆管识别的误报率。结果表明突增和突减两类异常向量的出现与天气的变化相符,证实了对异常向量分类的正确性。
3) 当DMA的出入口较多时,某些流量计的数据可能没有较强的变化规律,该方法的效果也因此显著下降。这与建设DMA时边界的划分、出入口的设置等因素有关,但此情况在我国真实管网中普遍存在,所以应注重对弱规律数据的研究以提升识别效果,而更为科学合理的DMA划分方法也值得进一步研究。
4) 在聚类过程中,参数值越大,对爆管事件越灵敏,但误报次数也会增加, 在DMA出入口相对较少且位置适宜的前提下,通过选定合适的参数,可以保证该方法对爆管事件较灵敏的同时获得较低误报率。
随着城市供水管网监测系统的进一步完善以及分区计量在我国的迅速发展,该技术将有广阔的应用前景。另外,可使用该技术对临近的多个独立DMA的用水量数据进行同步分析(爆管会造成某DMA用水量突然增加),将基于聚类的爆管识别方法应用在独立计量区域上。这都将有助于管网的漏损控制,同时降低爆管带来的其他负面影响。
[1] | Mounce S R, Khan A, Wood A S, et al. Sensor-fusion of hydraulic data for burst detection and location in a treated water distribution system[J]. Information Fusion, 2003, 4(3): 217–229. DOI:10.1016/S1566-2535(03)00034-4 |
[2] | Mounce S R, Boxall J B, Machell J. Development and verification of an online artificial intelligence system for detection of bursts and other abnormal flows[J]. Journal of Water Resources Planning and Management, 2010, 136(3): 309–318. DOI:10.1061/(ASCE)WR.1943-5452.0000030 |
[3] | Romano M, Kapelan Z, Savic D A. Automated detection of pipe bursts andother events in water distribution systems[J]. Journal of Water Resources Planning and Management, 2014, 140(4): 457–467. DOI:10.1061/(ASCE)WR.1943-5452.0000339 |
[4] | Ye G L, Fenner R A. Kalman filtering of hydraulic measurements for burst detection in water distribution systems[J]. Journal of Pipeline Systems Engineering and Practice, 2011, 2(1): 14–22. DOI:10.1061/(ASCE)PS.1949-1204.0000070 |
[5] | Ye G L, Fenner R A. Weighted least squares with expectation-maximization algorithm for burst detection in U.K. water distribution systems[J]. Journal of Water Resources Planning and Management, 2014, 140(4): 417–424. DOI:10.1061/(ASCE)WR.1943-5452.0000344 |
[6] | 凌文翠, 张涛, 强志民, 等. 城市供水管网独立计量区域的研究与应用进展[J]. 中国给水排水, 2011, 27(13): 46–50. LING Wencui, ZHANG Tao, QIANG Zhimin, et al. Research and application of district metering area for urban water distribution network[J]. China Water & Wastewater, 2011, 27(13): 46–50. (in Chinese) |
[7] | Wu Y, Liu S, Wu X, et al. Burst detection in district metering areas using a data driven clustering algorithm[J]. Water Research, 2016, 100: 28–37. DOI:10.1016/j.watres.2016.05.016 |
[8] | Tan P N, Steinbach M, Kumar V. Introduction to Data Mining[M]. New Jersey: Addison-Wesley, 2005. |
[9] | Rodriguez A, Laio A. Clustering by fast search and find of density peaks[J]. Science, 2014, 344(6191): 1492–1496. DOI:10.1126/science.1242072 |
[10] | Loureiro D, Amado C, Martins A, et al. Water distribution systems flow monitoring and anomalous event detection:A practical approach[J]. Urban Water Journal, 2016, 13(3): 242–252. DOI:10.1080/1573062X.2014.988733 |
[11] | Metz C E. Basic principles of ROC analysis[J]. Seminars in Nuclear Medicine, 1978, 8(4): 283–298. DOI:10.1016/S0001-2998(78)80014-2 |