基于时间序列形状分析的供水管网泄漏事件诊断方法

摘要

本公开提出一种基于时间序列形状分析的供水管网泄漏事件诊断方法，包括以下步骤：S1，利用供水管网监测数据构建时间序列，并确定时间序列的前向偏移量和后向偏移量；S2，利用所述监测数据，基于前向偏移量和后向偏移量构建历史子序列库；S3，利用历史子序列间的余弦夹角距离对历史子序列进行形状分析，删除历史子序列库中的形状异常子序列，构建参考子序列库；S4，构建新子序列，并将所述参考子序列库与新子序列在形状上进行比对，根据比对结果诊断泄漏事件。本公开基于时间序列形状分析的供水管网泄漏事件诊断方法，实现对供水管网中泄漏事件的准确判断。

权利要求

1.一种基于时间序列形状分析的供水管网泄漏事件诊断方法，包括以下步骤：S1，利用供水管网监测数据构建时间序列，并确定时间序列的前向偏移量nb和后向偏移量nf；

S2，利用所述监测数据，基于前向偏移量nb和后向偏移量nf构建历史子序列库；

S3，利用历史子序列间的余弦夹角距离对历史子序列进行形状分析，删除历史子序列库中的形状异常子序列，构建参考子序列库；

S4，构建新子序列，并将所述参考子序列库与新子序列在形状上进行比对，根据比对结果诊断泄漏事件。

2.根据权利要求1所述的方法，其中，在步骤S1中，利用供水管网监测数据构建时间序列，计算供水管网DMA中监测数据所构成的时间序列间的相关系数，确定序列的前向偏移量nb和后向偏移量nf。

3.根据权利要求1所述的方法，其中，所述步骤S1包括以下子步骤：S11，按日切分m天的流量或压力监测数据，构建m个时间序列，记每天的时间序列为Ti，i＝1，2，...，m；

S12，对时间序列Ti进行前、后向偏移，获得前、后向偏移的时间序列 和S13，计算时间序列Ti与前、后向偏移的时间序列 和 间的相关系数，获取前、后向偏移相关系数矩阵B和F；

S14，计算后向偏移相关系数矩阵F和前向偏移相关系数矩阵B每一列的均值，并计数均值中大于一预设值的个数，由此确定前向偏移量nb和后向偏移量nf。

4.根据权利要求3所述的方法，其中，

在子步骤S11中，若数据采集间隔为gmin，则时间序列Ti包含 个数据；

在子步骤S12中，所述后向偏移的时间序列 是指时间序列Ti向后滑动f个数据，即含有第(i+1)天的f个数据，同时删除第i天的f个数据，所述前向偏移是指时间序列Ti向前滑动b个数据，即 含有第(i-1)天的b个数据，同时删除第i天的b个数据；

在子步骤S13中，所述时间序列Ti与其后向偏移的时间序列 间后向偏移相关系数以及所述时间序列Ti与其前向偏移的时间序列 间的前向偏移相关系数 的计算公式为：其中，j＝1，2，...，m-1，cov表示时间序列间的协方差，σ表示时间序列中 个数据的方差；

在子步骤S14中，计算后向偏移相关系数矩阵F与前向偏移相关系数矩阵B每一列的均值，分别获得后向偏移相关系数均值向量sf与前向偏移相关系数均值向量sb，分别计数sf和sb中大于0.9的前nf和nb个元素，并将nf和nb作为最终的后向偏移量和前向偏移量。

5.根据权利要求4所述的方法，其中，所述步骤S2包括以下子步骤：S21，设定子序列长度l，利用m天的流量或压力监测数据，提取时刻t及其前nb+l-1个时刻和后nf个时刻的数据，构成第i天的时间序列片段，记为S22，将时间序列片段 切分为长度为l的历史子序列，构建历史子序列库Ht。

6.根据权利要求5所述的方法，其中，

在子步骤S21中，设定子序列长度l，在每个时刻t，提取一天中t时刻及其前(nb+l-1)个时刻和后nf个时刻的数据，每一天相应的数据所构成的时间序列片段记为 对于每天的前(nb+l-1)个时刻，仅提取后(m-1)天中相应的时间序列片段，即i＝2，3，...，m；对于每天的后nf个时刻，仅提取前(m-1)天中相应的时间序列片段，即i＝1，2，...，m-1，对于其他时刻，i＝1，2，...，m；l的取值范围为(nf+nb+1)到2(nf+nb)间的整数；

在子步骤S22中，设置长度为l的滑动窗口，则每个时间序列片段 可构成(nf+nb+1)个长度为l的子序列，使用时刻t所有时间序列片段构成的子序列构建时刻t的历史子序列库Ht；对于前(nb+l-1)个时刻和后nf个时刻，相应的历史子序列库Ht为(m-1)(nf+nb+1)行l列的矩阵，其他时刻的历史子序列库Ht则为m(nf+nb+1)行l列的矩阵。

7.根据权利要求6所述的方法，其中，所述步骤S3包括以下子步骤：S31，计算历史子序列库Ht中子序列间的余弦夹角距离，利用余弦夹角距离对历史子序列进行形状分析，获得距离矩阵Dt；

S32，删除距离矩阵Dt中的特定元素，逐行计算Dt中的最小值，构建最小距离向量dt；

S33，依据最小距离向量dt构建位置向量pt和天数向量p′t；

S34，设定初始化参数q，寻找p′t中前 个不重复的元素，每个元素记为S35，删除Ht中序列片段 所对应的形状异常子序列，构建参考子序列库Rt。

8.根据权利要求7所述的方法，其中，

在子步骤S31中，计算时刻t的历史子序列库Ht中子序列间的余弦夹角距离di，j，其计算公式为：式中， 和 分别表示时刻t历史子序列库中子序列 和 的模，其中j的取值与i相同，最终获得距离矩阵Dt，该矩阵为对称阵，行数和列数均等于Ht的行数，每一行或列表示第i个子序列与自身以及Ht中其他子序列间的余弦夹角距离；

在子步骤S32中，删除距离矩阵Dt中第i行的第 到第

个元素，同时计算剩余元素的最小值；对矩阵Dt的每行逐行进行

上述元素删除和剩余元素的最小值计算操作，并用所得到的最小距离构建最小距离向量dt；

在子步骤S33中，对最小距离向量dt进行降序排列，顺次记录排序后向量中每个元素在原向量dt中所在的位置并构建位置向量pt，然后将pt中的元素同时除以(nf+nb+1)，并对结果向上取整，最终得到天数向量p′t；

在子步骤S34中，设定初始化参数q，该参数表示m天中，有 天的数据可能存在形状上的异常 ；顺次寻找向量p ′t中前 个不重复的元素，每个元素记为在子步骤S35中，删除Ht中第 天的序列片段 所对应的子序列，构建时刻t的参考子序列库Rt。

9.根据权利要求8所述的方法，其中，所述步骤S4包括以下子步骤：S41，利用t时刻新采集的数据及该时刻t前l-1个时刻的数据构建新子序列Ct，将新子序列Ct与参考子序列库Rt合并，构建诊断矩阵Xt，获取Xt的距离矩阵D′t；

S42，删除D′t中的特定元素，然后逐步计算D′t中每一行的最小距离，构建最小距离向量d′t；

S43，判断d′t中的最后一个元素是否为最大值，若是，则t时刻发生泄漏事件；若否，则t时刻未发生泄漏事件。

10.根据权利要求9所述的方法，其中，

在子步骤S41中，将新子序列Ct作为矩阵的最后一行与参考子序列库Rt合并，构建诊断矩阵Xt，并获取Xt的距离矩阵D′t；对于前(nb+l-1)个时刻和后nf个时刻，D′t为行列均为的对称阵，其他时刻的D′t则为行列均为的对称阵；

在子步骤S42中，删除距离矩阵D′t中第i行的第 到第

个元素，然后计算D′t每一行的最小距离，以构成最小距离向量

d′t；

在子步骤S43中，所述d′t中的元素为d′i，i＝1，2，...，v，若 则判定时刻t发生泄漏事件，否则判定时刻t未发生泄漏事件；对于前(nb+l-1)个时刻和后nf个时刻，对于其他时刻，

说明书

基于时间序列形状分析的供水管网泄漏事件诊断方法

技术领域

[0001] 本公开属于城市供水管网领域与数据分析领域，具体涉及一种对供水管网泄漏事件的诊断方法，尤其涉及一种基于时间序列形状分析的供水管网泄漏事件诊断方法。

背景技术

[0002] 2016年，我国城镇供水管网漏损率约为18％，与《水污染防治行动计划》(“水十条”)要求的到2020年公共供水管网漏损率不超过10％的目标相比，仍有较大差距。供水管网的生命周期很长，服役期可达50年甚至100年。据统计，我国现有供水管网约50万km，这些管网将逐步进入高领服役期。预计到2030年，约55％的供水管网的管龄将超过30年，这将导致供水管网的泄漏更易发生，为供水管网的漏损控制造成极大障碍。

[0003] 大量工程经验表明，DMA(District Metering Area，独立计量分区)是控制漏损的有效手段。为更加有效地控制漏损，我国于2017年颁布了《城镇供水管网分区计量管理工作指南》，明确提出加快DMA地建设，加强对小区服务管线的漏损控制。

[0004] 快速发现泄漏点是有效减少管线泄漏水量的关键，同时也是降低供水管网漏损率的重要一环。由于材质、地面负荷、施工质量等因素影响，管线泄漏点在空间上分布不均匀。基于分布式光纤传感器的探漏技术(例如CN206504114U，CN106764460A)虽具有定位精度较高的特点，但应用于DMA中DN200以下小口径管线时，在经济性和工程适用性等方面不尽人意。以捕捉声学特征为主的技术(例如CN205919127U)，则存在时效性差的问题。此外，上述技术均需要专业的硬件设备，设备成本与管网改造费用高昂，不利于该种方法的大规模推广。

[0005] 管线泄漏会导致水量瞬增和压力突降，从理论上讲，可利用流量和压力监测设备读数的变化，判断是否存在泄漏。随着DMA的建设以及在线监测技术的广泛应用，供水管网逐渐积累了海量数据，如何通过数据挖掘，实现对泄漏的实时诊断，已成为供水管网领域的一大热点问题。在石油输配管道中，已经有基于人工神经网络的分类模型(例如CN106352244A)能够有效诊断管线是否发生泄漏。然而，分类模型的训练需要带有标签的数据(正常和泄漏两类)，属于有监督的学习方法。供水管网中的泄漏数据较少，而且由于管理等问题，自来水公司并未对所有泄漏事件进行记录，这都为分类模型的建立造成了障碍。因此，无需数据标签的非监督学习方法亟待开发。

发明内容

[0006] (一)要解决的技术问题

[0007] 本公开提供了一种基于时间序列形状分析的供水管网泄漏事件诊断方法，以至少部分解决以上所存在的技术问题。

[0008] (二)技术方案

[0009] 根据本公开的一个方面，提供了一种基于时间序列形状分析的供水管网泄漏事件诊断方法，包括以下步骤：

[0010] S1，利用供水管网监测数据构建时间序列，并确定时间序列的前向偏移量和后向偏移量；

[0011] S2，利用所述监测数据，基于前向偏移量和后向偏移量构建历史子序列库；

[0012] S3，利用历史子序列间的余弦夹角距离对历史子序列进行形状分析，删除历史子序列库中的形状异常子序列，构建参考子序列库；

[0013] S4，构建新子序列，并将所述参考子序列库与新子序列在形状上进行比对，根据比对结果诊断泄漏事件。

[0014] 在一些实施例中，在步骤S1中，利用供水管网监测数据构建时间序列，计算供水管网DMA中监测数据所构成的时间序列间的相关系数，确定序列的前向偏移量nb和后向偏移量nf。

[0015] 在一些实施例中，所述步骤S1包括以下子步骤：

[0016] S11，按日切分m天的流量或压力监测数据，构建m个时间序列，记每天的时间序列为Ti，i＝1，2，...，m；

[0017] S12，对Ti进行前、后向偏移，获得前、后向偏移的时间序列 和

[0018] S13，计算Ti与前、后向偏移的时间序列间的相关系数，获取前、后向偏移相关系数矩阵B和F；

[0019] S14，计算后向偏移相关系数矩阵F和前向偏移相关系数矩阵B每一列的均值，并计数均值中大于一预设值的个数，由此确定前向偏移量nb和后向偏移量nf。

[0020] 在一些实施例中，在子步骤S11中，若数据采集间隔为gmin，则时间序列Ti包含个数据；

[0021] 在子步骤S12中，所述后向偏移的时间序列 是指时间序列Ti向后滑动f个数据，即 含有第(i+1)天的f个数据，同时删除第i天的f个数据，所述前向偏移是指时间序列Ti向前滑动b个数据，即 含有第(i-1)天的b个数据，同时删除第i天的b个数据；

[0022] 在子步骤S13中，所述时间序列间后向偏移相关系数 前向偏移相关系数 的计算公式为：

[0023]

[0024]

[0025] 其中，j＝1，2，...，m-1，cov表示时间序列间的协方差，σ表示时间序列中 个数据的方差；

[0026] 在子步骤S14中，计算后向偏移相关系数矩阵F与前向偏移相关系数矩阵B每一列的均值，分别获得后向偏移相关系数均值向量sf与前向偏移相关系数均值向量sb，分别计数sf和sb中大于0.9的前nf和nb个元素，并将nf和nb作为最终的后向偏移量和前向偏移量。

[0027] 在一些实施例中，所述步骤S2包括以下子步骤：

[0028] S21，设定子序列长度l，利用m天的流量或压力监测数据，提取时刻t及其前nb+l-1个时刻和后nf个时刻的数据，构成第i天的时间序列片段，记为

[0029] S22，将时间序列片段 切分为长度为l的历史子序列，构建历史子序列库Ht。

[0030] 在一些实施例中，在子步骤S21中，设定子序列长度l，在每个时刻t，提取一天中t时刻及其前(nb+l-1)个时刻和后nf个时刻的数据，每一天相应的数据所构成的时间序列片段记为 对于每天的前(nb+l-1)个时刻，仅提取后(m-1)天中相应的时间序列片段，即i＝2，3，...，m；对于每天的后nf个时刻，仅提取前(m-1)天中相应的时间序列片段，即i＝1，
2，...，m-1，对于其他时刻，i＝1，2，...，m；l的取值范围为(nf+nb+1)到2(nf+nb)间的整数；

[0031] 在子步骤S22中，设置长度为l的滑动窗口，则每个时间序列片段 可构成(nf+nb+1)个长度为l的子序列，使用时刻t所有时间序列片段构成的子序列构建时刻t的历史子序列库Ht；对于前(nb+l-1)个时刻和后nf个时刻，相应的历史子序列库Ht为(m-1)(nf+nb+1)行l列的矩阵，其他时刻的历史子序列库Ht则为m(nf+nb+1)行l列的矩阵。

[0032] 在一些实施例中，所述步骤S3包括以下子步骤：

[0033] S31，计算历史子序列库Ht中子序列间的余弦夹角距离，利用余弦夹角距离对历史子序列进行形状分析，获得距离矩阵Dt；

[0034] S32，删除距离矩阵Dt中的特定元素，逐行计算Dt中的最小值，构建最小距离向量dt；

[0035] S33，依据最小距离向量dt构建位置向量pt和天数向量p′t；

[0036] S34，设定初始化参数q，寻找p′t中前 个不重复的元素，每个元素记为

[0037] S35，删除Ht中序列片段 所对应的形状异常子序列，构建参考子序列库Rt。

[0038] 在一些实施例中，在子步骤S31中，计算时刻t的历史子序列库Ht中子序列间的余弦夹角距离di，j，其计算公式为：

[0039]

[0040] 式中， 和 分别表示时刻t历史子序列库中子序列 和 的模，其中j的取值与i相同，最终获得距离矩阵Dt，该矩阵为对称阵，行数和列数均等于Ht的行数，每一行或列表示第i个子序列与自身以及Ht中其他子序列间的余弦夹角距离；

[0041] 在子步骤S32中，删除距离矩阵Dt中第i行的第 到第 个元素，同时计算剩余元素的最小值；对Dt的每行逐行进行上
述元素删除和剩余元素的最小值计算操作，并用所得到的最小距离构建最小距离向量dt；

[0042] 在子步骤S33中，对最小距离向量dt进行降序排列，顺次记录排序后向量中每个元素在原向量dt中所在的位置并构建位置向量pt，然后将pt中的元素同时除以(nf+nb+1)，并对结果向上取整，最终得到天数向量p′t；

[0043] 在子步骤S34中，设定初始化参数q，该参数表示m天中，有 天的数据可能存在形状上的异常；顺次寻找向量p′t中前 个不重复的元素，每个元素记为

[0044] 在子步骤S35中，删除Ht中第 天的序列片段 所对应的子序列，构建时刻t的参考子序列库Rt。

[0045] 在一些实施例中，所述步骤S4包括以下子步骤：

[0046] S41，利用t时刻新采集的数据及该时刻t前l-1个时刻的数据构建新子序列Ct，将新子序列Ct与参考子序列库Rt合并，构建诊断矩阵Xt，获取Xt的距离矩阵D′t；

[0047] S42，删除D′t中的特定元素，然后逐步计算D′t中每一行的最小距离，构建最小距离向量d′t；

[0048] S43，判断d′t中的最后一个元素是否为最大值，若是，则t时刻发生泄漏事件；若否，则t时刻未发生泄漏事件。

[0049] 在一些实施例中，在子步骤S41中，将新子序列Ct作为矩阵的最后一行与参考子序列库Rt合并，构建诊断矩阵Xt，并获取Xt的距离矩阵D′t；对于前(nb+l-1)个时刻和后nf个时刻，D′t为行列均为 的对称阵，其他时刻的D′t则为行列均为 的对称阵；

[0050] 在子步骤S42中，删除距离矩阵D′t中第i行的第 到第 个元素，然后计算D′t每一行的最小距离，以构成最小距离向量
d′t；

[0051] 在子步骤S43中，所述d′t中的元素为d′i，i＝1，2，...，v，若 则判定时刻t发生泄漏事件，否则判定时刻t未发生泄漏事件；对于前(nb+l-1)个时刻和后nf个时刻， 对于其他时刻，

[0052] (三)有益效果

[0053] 从上述技术方案可以看出，本公开基于时间序列形状分析的供水管网泄漏事件诊断方法至少具有以下有益效果其中之一：

[0054] (1)本公开的核心之一在于充分分析了DMA时间序列数据间的形状相似性。由于DMA内用户较为固定，每日相同时段内的用水行为相对一致，即使因为天气、季节等因素造成用水总量上的增减变化，但每日相同时段内的水量增减趋势(即形状)较为相似。因此，截取每日相同时段内的监测数据而获取的子序列间具有较高的相似性。

[0055] (2)利用余弦夹角距离对子序列的形状进行分析，寻找形状相似度低的子序列作为异常序列，以达到泄漏事件诊断的目的，无需获取数据标签(正常或泄漏)，即可自动识别出表征泄漏事件的监测数据，属于无监督的学习方法。

[0056] (3)本公开中的泄漏诊断方法具有较强的鲁棒性，能够保证较少的误报数量。在构建历史子序列库时，对每个时刻的子序列加入前、后向偏移量，构成多个前、后平移的子序列，使子序列库更为丰富。这使得流量(或压力)时间序列发生前、后平移(由用户用水行为造成)时，该诊断方法也不会产生误报。

附图说明

[0057] 图1示出了一种基于时间序列形状分析的供水管网泄漏事件诊断方法的流程图；

[0058] 图2示出了确定时间序列前、后向偏移量的流程图；

[0059] 图3示出了构建历史子序列库的流程图；

[0060] 图4示出了构建参考子序列库的流程图；

[0061] 图5示出了在线诊断泄漏事件的流程图。

具体实施方式

[0062] 为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。

[0063] 针对上述问题，本公开的目的是提供一种基于时间序列形状分析的供水管网泄漏事件诊断方法，实现以无监督学习的方法对供水管网中是否发生泄漏事件进行准确的判断。

[0064] 所述基于时间序列形状分析的供水管网泄漏事件诊断方法，包括以下步骤：利用供水管网监测数据构建时间序列，并确定时间序列的前向偏移量和后向偏移量；利用所述监测数据，基于前向偏移量和后向偏移量构建历史子序列库；利用历史子序列间的余弦夹角距离对历史子序列进行形状分析，删除历史子序列库中的形状异常子序列，构建参考子序列库；构建新子序列，并将所述参考子序列库与新子序列在形状上进行比对，根据比对结果在线诊断泄漏事件。本公开的核心在于充分分析了DMA时间序列数据间的形状相似性，属于无监督的学习方法。

[0065] 具体的，本公开基于时间序列形状分析的供水管网泄漏事件诊断方法，包括以下步骤：

[0066] S1，计算供水管网DMA中流量(或压力)监测数据所构成的时间序列间的相关系数，确定序列的前向偏移量nb和后向偏移量nf。此步骤为序列前、后向偏移量的确定，具体包括：

[0067] S11，将m天的流量(或压力)历史监测数据按日切分，构成m个时间序列；记每日的时间序列为Ti，i＝1，2，...，m；若数据采集间隔为gmin，则时间序列Ti包含 个数据；其中m的取值范围为30到90之间整数，g一般为5或15；

[0068] S12，对Ti进行前、后向偏移，并将进行前、后向偏移的时间序列分别记为 和特别地，对于T1，仅进行后向偏移，对于Tm，仅进行前向偏移。其中，后向偏移是指时间序列Ti向后滑动f个数据， 将含有第(i+1)天的f个数据，同时删除第i天的f个数据，i＝1，2，...，m-1， 前向偏移是指时间序列Ti向前滑动b个数据， 将含有
第(i-1)天的b个数据，同时删除第i天的b个数据，i＝2，3，...，m，

[0069] S13，计算Ti与 之间以及Ti与 之间的相关系数，j＝1，2，...，m-1，k＝2，3，...，m；最终获取后向偏移相关系数矩阵F与前向偏移相关系数矩阵B；F与B均有m(m-1)行， 列，每一列按照偏移量f(或b)的大小升序排列；时间序列间后向偏移相关系数 (或前向偏移相关系数 )的计算公式为：

[0070]

[0071]

[0072] 其中，cov表示时间序列间的协方差，σ表示时间序列中 个数据的方差；

[0073] S14，计算F与B每一列的均值，分别获得后向偏移相关系数均值向量sf与前向偏移相关系数均值向量sb；分别计数sf和sb中大于0.9的前nf和nb个元素，并将nf和nb作为最终的后向偏移量和前向偏移量；由此，可确定与原时间序列的变化趋势(形状)具有高度相似性的偏移时间序列。当然，本领域技术人员也可以对此数值0.9进行适当调整。

[0074] S2，设定子序列长度l，利用m天的流量(或压力)历史监测数据，按照时刻(共个时刻)，构建 个历史子序列库，每个时刻的历史子序列库记为Ht。此步骤为历史子序列库的构建，具体包括：

[0075] S21，在每个时刻t，提取一天中该时刻及其前(nb+l-1)个时刻和后nf个时刻的数据，每一天相应的数据所构成的时间序列片段记为 上标i代表天数；特别地，对于每天的前(nb+l-1)个时刻，仅提取后(m-1)天中符合条件的数据构成时间序列片段，i＝2，3，...，m，对于每天的后nf个时刻，仅提取前(m-1)天中符合条件的数据构成时间序列片段，i＝1，2，...，m-1；对于其他时刻，i＝1，2，...，m；l的取值范围为(nf+nb+1)到2(nf+nb)间的整数；
其中，符合条件的数据是指，在每个时刻，所提取的该时刻及其前(nb+l-1)个时刻和后nf个时刻的数据。

[0076] S22，设置长度为l的滑动窗口，则每个时间序列片段 可构成(nf+nb+1)个长度为l的子序列，使用时刻t所有时间序列片段构成的子序列构建时刻t的历史子序列库Ht；对于前(nb+l-1)个时刻和后nf个时刻，相应的历史子序列库Ht为(m-1)(nf+nb+1)行l列的矩阵，其他时刻的历史子序列库Ht则为m(nf+nb+1)行l列的矩阵。

[0077] S3，使用余弦夹角距离对历史子序列进行形状分析，其中，形状相似的历史子序列间距离较小，反之，距离较大；若某一子序列与其他子序列间距离的最小值相对较大(例如大于一阈值)，则该子序列为形状异常的子序列；删除历史子序列库Ht中异常的子序列，构成 个参考子序列库，每个时刻的参考子序列库记为Rt。此步骤为参考子序列库的构建，具体包括：

[0078] S31，记时刻t的历史子序列库中的子序列为 对于前(nb+l-1)个时刻和后nf个时刻，i＝1，2，...，(m-1)(nf+nb+1)，对于其他时刻，i＝1，2，...，m(nf+nb+1)；计算时刻t的历史子序列库中子序列间的余弦夹角距离di，j，其计算公式为：

[0079]

[0080] 式中， 和 分别表示时刻t历史子序列库中子序列 和 的模，其中j的取值与i相同；最终获得距离矩阵Dt，该矩阵为对称阵，行数和列数均等于Ht的行数，每一行(或列)表示第i个子序列与自身以及Ht中其他子序列间的余弦夹角距离；

[0081] S32，删除距离矩阵Dt中第i行中的第 到第个元素(这些元素代表第i个子序列与自身的距离值，以及第i个
子序列与跟其来自同一天的子序间的距离值)，然后计算剩余元素的最小值；对Dt的每行逐行进行上述元素删除和最小值计算操作，并用所得到的最小距离构建最小距离向量dt。对于前(nb+l-1)个时刻和后nf个时刻，i＝1，2，...，(m-1)(nf+nb+1)，相应的dt包含(m-1)(nf+nb+1)个最小距离；对于其他时刻，i＝1，2，...，m(nf+nb+1)，相应的dt则包含m(nf+nb+1)个最小距离；

[0082] S33，对最小距离向量dt进行降序排列，顺次记录排序后向量中每个元素在原向量dt中所在的位置并构建位置向量pt，然后将pt中的元素同时除以(nf+nb+1)，并对结果向上取整，最终得到天数向量p′t；p′t中的元素表示排序后的子序列对应的天数，排在越靠前的元素，表明来自该天的子序列越异常；对于前(nb+l-1)个时刻和后nf个时刻，向量p′t中各个元素的取值范围为1到(m-1)间的整数，对于其他时刻，p′t中各个元素的取值范围为1到m间的整数；

[0083] S34，设定初始化参数q，该参数表示m天中，有 天的数据可能存在形状上的异常，q的取值范围为[0.05，0.15]；顺次寻找向量p′t中前 个不重复的元素，每个元素记为 是指对q和m的乘积向上取整；

[0084] S35，删除Ht中第 天的序列片段 所对应的子序列(相对于历史子序列库中的其他子序列而言，所删除的这些子序列为在形状上是异常的子序列，即异常子序列)，构建时刻t的参考子序列库Rt；对于前(nb+l-1)个时刻和后nf个时刻，相应的参考子序列库Rt为行l列的矩阵，其他时刻的参考子序列库Rt则为行l列的矩阵。

[0085] S4，利用时刻t的参考子序列库，与时刻t的新子序列在形状上进行比对，从而判断时刻t是否发生泄漏事件；时刻t的新子序列由时刻t新采集的数据及该时刻前(l-1)个时刻的数据(监测数据)构成，记为Ct。此步骤为泄漏事件的在线诊断，具体包括：

[0086] S41，将Ct作为矩阵的最后一行与Rt合并，构建诊断矩阵Xt；使用步骤S31中的计算方法，获取Xt的距离矩阵D′t，对于前(nb+l-1)个时刻和后nf个时刻，D′t为行列均为的对称阵，其他时刻的D′t则为行列均为的对称阵；

[0087] S42，使用步骤S32中的计算方法逐行删除D′t中的元素，然后计算每行的最小距离，以构成最小距离向量d′t；特别地，对于D′t的最后一行，不做上述删除操作，直接计算该行除零以外的最小距离；

[0088] S43，记d′t中的元素为d′t，i＝1，2，...，v，若 则判定时刻t发生泄漏事件，否则判定时刻t未发生泄漏事件；对于前(nb+l-1)个时刻和后nf个时刻，对于其他时刻，

[0089] 在以下实施方式中，利用Matlab 9.1作为泄漏事件诊断方法的开发平台，通过在M文件编辑器中编写相应的代码即可实现对DMA监测数据的读取、存储和分析。本实施例利用Matlab 9.1强大的计算能力和丰富的内置函数，提高了开发效率，缩短了开发时间。

[0090] 图1示出了本公开基于时间序列形状分析的供水管网泄漏事件诊断方法的流程图，为了实现泄漏事件的诊断，首先应通过计算时间序列的相关系数确定序列的前、后向偏移量；然后使用历史监测数据构建历史子序列库；随后，基于子序列间的余弦夹角距离，删除历史子序列库中与其他子序列形状相似性较低的子序列，构建参考子序列库；最后，在形状上比对新子序列与参考子序列库中的子序列，进行泄漏事件的在线诊断。

[0091] 具体地，作为示例，选取一真实的DMA进行实施。所选DMA仅有1个入水口，使用电磁流量计进行计量，流量数据每15min采集并传输一次，即g＝15，

[0092] 一天中共有96个时刻。收集该管网4个月的历史数据(2016年4月1日到2016年7月31日)用以实施本公开。使用4-5月的数据构建历史子序列库，6-7月的数据用于在线诊断。
在6-7月间，发生消火栓损坏事故两次，管道冲洗一次。这些事件可以视为泄漏事件，从而验证本公开的诊断效果。

[0093] 本实施例中所用3次事件的具体信息如表1所示：

[0094] 表1事件描述

[0095]

[0096] (1)计算供水管网DMA中流量(或压力)监测数据所构成的时间序列间的相关系数，确定序列的前向偏移量nb和后向偏移量nf：

[0097] 由于使用4-5月(共61天)的历史监测数据，则可切分为61个时间序列，即m＝61；对61个时间序列进行前、后向偏移，表2示出了对第2天的时间序列进行前、后向偏移的部分结果。

[0098] 表2时间序列前后偏移部分结果示例

[0099]

[0100]

[0101] 计算Ti与 之间以及Ti与 之间的相关系数，获得的后向偏移相关系数矩阵F与前向偏移相关系数矩阵B的形式相同，表.3示出了后向偏移相关系数矩阵F的部分结果。

[0102] 表3后向偏移相关系数矩阵F的部分结果示例

[0103]

[0104] 根据后向偏移相关系数矩阵F与前向偏移相关系数矩阵B计算后向偏移相关系数均值向量sf与前向偏移相关系数均值向量sb，部分结果如表4所示。

[0105] 表4前后偏移相关系数均值向量部分结果示例

[0106]f/b 1 2 3 4 5 6 7 8 9 10 11 12 13 …
sf 0.95 0.93 0.89 0.84 0.79 0.73 0.67 0.60 0.54 0.47 0.41 0.34 0.27 …
sb 0.95 0.93 0.89 0.85 0.79 0.73 0.67 0.60 0.53 0.47 0.40 0.33 0.27 …

[0107] 基于上述结果，sf和sb中的前两个元素大于0.9，即nf＝nb＝2。

[0108] (2)设定子序列长度l，利用m天的流量(或压力)历史监测数据，按照时刻(共个时刻)，构建 个历史子序列库：

[0109] 本实施例中选取子序列长度为6，即l＝6，则每天时刻t的时间序列片段 可构成5个长度为6的子序列；对于前7个时刻(0:00到1:30)和后2个时刻(23:30和23:45)，相应的历史子序列库Ht为300行6列的矩阵，其他时刻的历史子序列库Ht则为305行6列的矩阵；表5示出了12:00(第49个时刻)历史子序列库H49的一部分；可以发现，与其他子序列相比，第20天的子序列中有元素异常增大。

[0110] 表5历史子序列库H49部分示例

[0111]

[0112] (3)设定初始化参数q，使用余弦夹角距离评估子序列间的形状相似性，定义与其他子序列间距离的最小值较大的子序列为形状相似性低的子序列(即异常子序列)，删除历史子序列库Ht中异常的子序列，构成参考子序列库Rt：

[0113] 本实施例以H49，示例参考子序列库的构建。表6示出了由H49得到的距离矩阵D49的一部分以及相应的最小距离向量d49；可以发现，来自第20天的子序列(96-100行)与其他子序列间的距离明显较大。由于余弦夹角距离的值普遍偏小，表6中的距离全部放大了一万倍。

[0114] 表6距离矩阵D49及其相应的最小距离向量d49部分示例

[0115]

[0116] 表7示出了对d49进行降序排列后获取的位置向量pt以及经过进一步处理后得到的向量p′t的部分结果。

[0117] 表7位置向量pt以及向量p′t的部分结果示例

[0118]pt 96 97 98 99 100 63 225 62 64 65 224 223 92 …
p′t 20 20 20 20 20 13 45 13 13 13 45 45 19 …

[0119] 由向量p′t的结果可知，H49中第20天的子序列形状相似性与其他子序列相比最低(最异常)，表5以及表6中现实的结果相一致；本实施例中，设置初始化参数为0.15，则需要删除10天 的子序列；顺次寻找向量p′t中前10个不重复的元素，分别为：20，13，45，19，53，9，61，35，12和4。因此，删除这些天对应的共50个子序列，构成12:00的参考子序列库R49。基于上述结果，表5中示出的H49的第96-100行以及301-305行均被删除。
最终，R49仅包含255个子序列，其他时刻的参考子序列库均按照上述过程构建。

[0120] (4)利用时刻t的参考子序列库，与时刻t的新子序列在形状上进行比对，从而判断时刻t是否发生泄漏事件：

[0121] 本实施例以7月14日的消火栓损坏为例，可知该天12:00的子序列应该被识别为异常，记相应的序列C49＝(273.2，282.5，286.4，337.0，409.5，368.3)。

[0122] 将C49作为矩阵的最后一行与R49合并，构成诊断矩阵X49，对X49中的所有子序列进行距离计算，所得到的距离矩阵D′49与D49形式相同。与步骤S3类似，逐行删除D′49中的元素，然后计算每行的最小距离，以构成最小距离向量d′49；特别地，对于D′49的最后一行(即第256行)，不做上述删除操作，直接计算该行除零以外的最小距离。表8示出了最小距离向量d′49的 一部 分 ，与 表 6 相同 ，表 8 中的 距离 全 部 放 大 了 一万 倍 。经 计 算 ，因此，判断确定7月14日12:00发生了泄漏事件。

[0123] 表8最小距离向量d′49部分示例

[0124]

[0125] 在对6-7月所有流量数据的5856次诊断结果中，除了准确诊断了表1中所列的3次泄漏事件，本公开中的诊断方法还进行了29次报警，其中有12次均是由于信号传输问题所导致的，其他17次报警的原因未知，视为误报。相对于5856次诊断，仅出现17次误报，误报数量较少。

[0126] 至此，已经结合附图对本公开实施例进行了详细描述。依据以上描述，本领域技术人员应当对本公开有了清楚的认识。

[0127] 需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。

[0128] 当然，根据实际需要，本公开方法还包含其他的步骤，由于同本公开的创新之处无关，此处不再赘述。

[0129] 以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。