本发明公开了一种基于VARX模型的供水管网异常事件侦测方法。本发明首先进行测点分组，确定输入样本。其次按分组建立VARX模型。然后进行压力预测并计算差异，并计算差异的平均值和标准差。最后按照异常事件判定规则侦测异常事件。本发明基于VARX模型采用差异分析的供水管网异常事件侦测方法，具有抗干扰性强，侦测能力强等特点。
1.一种基于VARX模型的供水管网异常事件侦测方法，其特征在于该方法包括以下步骤：步骤1.测点分组，确定输入样本，具体是：建立供水管网预测数据库；其中输入数据包括：测量点压力、流量值；输出数据包括：预测压力、流量值；
将监测点通过相关性分组，相关性计算公式如下：其中x，y是两组时间序列数据；Ex，Ey是对应的两组时间序列数据的期望值；Dx，Dy是对应的两组时间序列的方差；corr(x,y)则是对应的两组时间序列数据的相关系数，用来表征相关性；相关系数的取值范围在0～1，0表示不相关，反之1表示相关性明显；
根据监测点分组确定输入样本为{xt,yt}i，其中xt为d维外生变量，yt为k维内生变量，即d个外生变量的监测点k个内生变量的监测点，t为时间单位，i表示第i组监测点；
为每一组监测点的每一次预测建立VARX模型，根据输入样本数据{xt,yt}i进行模型定阶并建立模型；
其中 为残差，T为样本周期，p为模型滞后项，AIC值取最小值时的p值为模型最佳选择值；
其中：xt为d维外生变量；Θj为k×d维待估计的系数矩阵；yt为k维的内生变量；c为k维的常数项；Φi为k×k维待估计的系数矩阵；p1,q1为滞后阶数；{εt}为k维白噪声序列；
VARX模型的参数估计使用极大似然法得到，设定：β＝vec([c Φ1...Φp1 Θ1...Θq1])， 则上式可表示为：yt＝Wt-1β+εt
根据样本数据和建立好的模型得到每个监测点的预测值,并计算监测点n天内每一时刻预测值和观测值的差异；其中VARX模型预测值为(Y1'j,Y2'j,...,Yij')，观测值为(x1j,x2j,...,xij)，i＝1...k，j＝1...n，k表示监测点压力数据个数，则差异可表示为：Δyij＝(x1j-Y′1j,x2j-Y′2j,...,xij-Y′ij)步骤4.计算差异的平均值和标准差
计算每个监测点在M天中差异的平均值和标准差；假定压力时间序列的采样频率为u min/次，M天的差异可表示如步骤3，那么一天差异的平均值 可由下面的公式计算得到，其中q＝1...u...1440/u；
标准差(σ1,σ2,...,σq)可通过下式计算得到：步骤5.侦测异常事件
按照下述异常事件判定规则，实时侦测供水管网的异常情况：a.任一时间点压力、流量观测值低于-4σ边界；
b.连续2个时间点压力、流量观测值低于-3σ边界；
c.连续4个时间点压力、流量观测值低于-2σ边界；
d.连续8个时间点压力、流量观测值低于-1σ边界；
一旦观测值触动判定规则，则认定管网异常事件发生，记录发生时刻，并触发相应预警。
一种基于VARX模型的供水管网异常事件侦测方法
技术领域
[0001] 本发明属于城市供水领域，具体是一种基于VARX模型的供水管网异常事件侦测方法。
背景技术
[0002] 暗漏和爆管是城市供水管网常见的两种异常事件，前者具有时间持续长、规模小、不易察觉的特点，而后者则正好相反，具有突发性、不确定性，规模大的特点。及时发现供水管网中存在的漏损事件，准确定位事件发生地点，迅速采取措施防止事态恶化，对供水管网安全运行来说显得尤为重要。
[0003] 目前供水管网异常事件侦测主要基于微观水力学模型和宏观数据模型两种。其中微观模型需要已知管网系统的详细资料，如管网的拓扑结构、管段的材质、长度、直径、摩阻等具体信息，国内基于微观模型的应用研究较少。基于宏观模型的研究主要分为边界分析和差异分析两种方法，边界分析方法易受噪声的干扰而差异分析方法抗干扰能力较强。
[0004] 有鉴于此，使用差异分析方法可以较好地对供水管网异常事件进行侦测。本发明基于VARX模型高精度预测进行差异分析法侦测供水管网异常事件。
发明内容
[0005] 本发明针对现有技术的不足，提供了一种基于VARX模型的供水管网异常事件侦测方法。其中VARX(a Vector Auto-Regressive with eXogenous variables)模型将单向因果关系的变量作为外生变量加入VAR模型中，可以有效地增强模型的抗扰动性，同时该模型具有较高的精度。
[0006] 为实现以上目的，本发明采取以下步骤：
[0007] 步骤1.测点分组，确定输入样本
[0008] 建立供水管网预测数据库。其中输入数据包括：测量点压力、流量值等；输出数据包括：预测压力、流量值等。
[0009] (1)监测点分组
[0010] 将监测点通过相关性分组，相关性计算公式如下：
[0011]
[0012] 说明：x，y是两组时间序列数据；Ex，Ey是对应的两组时间序列数据的期望值；Dx，Dy是对应的两组时间序列的方差；corr(x,y)则是对应的两组时间序列数据的相关系数，用来表征相关性。相关系数的取值范围在0～1，0表示不相关，反之1表示相关性明显，将相关性明显的监测点归为一类。
[0013] (2)确定输入样本
[0014] 根据监测点分组确定输入样本为{xt,yt}i，其中xt为d维外生变量，yt为k维内生变量，即d个外生变量的监测点k个内生变量的监测点，t为时间单位，i表示第i组监测点。
[0015] 步骤2.按分组建立VARX模型
[0016] 为每一组监测点的每一次预测建立VARX模型，根据输入样本数据{xt,yt}i进行模型定阶并建立模型。
[0017] (1)VARX模型定阶
[0018] 采用AIC(Akaike Information Criterion)定阶准则。其定义如下：
[0019]
[0020] 其中 为残差，T为样本周期，p为模型滞后项，AIC值取最小值时的p值为模型最佳选择值。
[0021] (2)建立VARX模型
[0022] VARX模型的表达式为：
[0023]
[0024] 其中：xt为d维外生变量；Θj为k×d维待估计的系数矩阵。yt为k维的内生变量；c为k维的常数项；Φi为k×k维待估计的系数矩阵；p1,q1为滞后阶数；{εt}为k维白噪声序列。
[0025] VARX模型的参数估计使用极大似然法得到，设定：β＝vec([c Φ1...Φp1Θ1...Θq1])， 则上式可表示为：
[0026] yt＝Wt-1β+εt
[0027] 随后根据极大似然算法有：
[0028]
[0029] 从β和 方面我们可以衍生模型参数如下：
[0030]
[0031]
[0032] 步骤3.压力预测并计算差异
[0033] 根据样本数据和建立好的模型得到每个监测点的预测值,并计算监测点n天内每一时刻预测值和观测值的差异。其中VARX模型预测值为(Y′1j,Y′2j,...,Y′ij)，观测值为(x1j,x2j,...,xij)，i＝1...k(k表示监测点压力数据个数)，j＝1...n，则差异可表示为：
[0034] Δyij＝(x1j-Y′1j,x2j-Y′2j,...,xij-Y′ij)
[0035] 步骤4.计算差异的平均值和标准差
[0036] 计算每个监测点在M天中差异的平均值和标准差。假定压力时间序列的采样频率为u min/次，M天的差异可表示如步骤3，那么一天差异的平均值 可由下面的公式计算得到，其中q＝1...u...1440/u。
[0037]
[0038] 标准差(σ1,σ2,...,σq)可通过下式计算得到：
[0039]
[0040] 步骤5.侦测异常事件
[0041] 按照下述异常事件判定规则，实时侦测供水管网的异常情况：
[0042] a.任一时间点压力、流量观测值低于-4σ边界；
[0043] b.连续2个时间点压力、流量观测值低于-3σ边界；
[0044] c.连续4个时间点压力、流量观测值低于-2σ边界；
[0045] d.连续8个时间点压力、流量观测值低于-1σ边界。
[0046] 一旦观测值触动判定规则，则认定管网异常事件发生，记录发生时刻，并触发相应预警。
[0047] 本发明的有益效果：本发明基于VARX模型采用差异分析的供水管网异常事件侦测方法，具有抗干扰性强，侦测能力强等特点。
附图说明
[0048] 图1：DMA区10个监测点分布图；
具体实施方式
[0049] 为使本发明实现的技术手段与创作特征易于表达，下面结合附图和实例，对本发明的实时方式进一步详述。
[0050] 本实例中考虑某DMA(District Metering Area)区，有10个有效压力测点，其地理位置如图1。同时引入三个入水口压力点(分别以P1，P2，P3表示)和DMA区用水量作为外生变量。压力数据时间范围从X年3月20日至X年4月3日，其中使用3月20日至4月2日数据建立模型对4月3日“5次模拟爆管试验”异常事件进行侦测。这里，仅用压力值进行示例说明，但本发明方法同样适用于流量值。
[0051] 步骤1.确定输入数据
[0052] 建立供水管网预测数据库。其中输入数据包括：测量点压力值、流量值等；输出数据包括：预测压力值等。
[0053] (1)监测点分组
[0054] 将压力监测点和入水口压力点通过相关性分组，相关性计算公式如下：
[0055]
[0056] 说明：x，y是两组时间序列数据；Ex，Ey是对应的两组时间序列数据的期望值；Dx，Dy是对应的两组时间序列的方差；corr(x,y)则是对应的两组时间序列数据的相关系数。得到压力监测点分组如下表1。
[0057] 表1：监测点分组情况表
[0058]组号 监测点成员
第一组 No.3、No.5、No.6、No.7、No.10、P2
第二组 No.1、No.8、No.9、P3
第三组 No.2、No.4、P1
[0059] (2)确定输入样本
[0060] 根据监测点分组确定输入样本为{xt,yt}i，其中xt为d维外生变量，yt为k维内生变量，t为时间单位，i表示第i组监测点。第一组监测点的输入样本为2维外生变量，包括P2入水口压力测点和DMA区用水量，5维的内生变量即为5个压力监测点，同理可得其余两组输入样本情况，采样频率为5min/次。
[0061] 步骤2.建立VARX模型
[0062] 根据输入样本确定模型阶数和模型，采用五天的数据预测一天压力监测点的数据，所以需要建立3×10个模型，即三组监测点别需要建立预测3.25～4.3日10天的VARX压力预测模型。
[0063] (1)VARX模型定阶
[0064] 采用AIC(Akaike Information Criterion)定阶准则。其定义如下：
[0065]
[0066] 其中 为残差，T为样本周期，p为模型滞后项，AIC值取最小值时的p值为模型最佳选择值。计算得到模型阶数为p＝4。
[0067] (2)建立VARX模型
[0068] VARX模型的表达式为：
[0069]
[0070] 其中：xt为d维外生变量；Θj为k×d维待估计的系数矩阵。yt为k维的内生变量；c为k维的常数项；Φi为k×k维待估计的系数矩阵；p1,q1为滞后阶数；{εt}为k维白噪声序列。
[0071] VARX模型的参数估计使用极大似然法得到，设定：β＝vec([c Φ1...Φp1Θ1...Θq1])， 则上式可表示为：
[0072] yt＝Wt-1β+εt
[0073] 随后根据极大似然算法有：
[0074]
[0075] 从β和 方面我们可以衍生模型参数如下：
[0076]
[0077]
[0078] 步骤3.压力预测并计算差异
[0079] 根据样本数据和建立好的模型得到每个监测点的预测值,并计算监测点n天内每一时刻预测值和观测值的差异。得到预测数据VARX模型预测值(Y′1j,Y′2j,...,Y′ij)，观测值为(x1j,x2j,...,xij)，其中i＝1...k(k表示监测点压力数据个数，第一组k＝5，第二组k＝3，第三组k＝2)，j＝1:10，则差异可表示为：
[0080] Δyij＝(x1j-Y′1j,x2j-Y′2j,...,xij-Y′ij)
[0081] 步骤4.计算差异的平均值和标准差
[0082] 计算每个监测点在9天中每一采样点的平均值和标准差。采样频率为5min/次，则9天的压力数据可表示为(x1j,x2j,...,xij)，其中i＝1:k(第一组k＝5，第二组k＝3，第三组k＝2)，j＝1:9，那么一天差异的平均值 可由下面的公式计算得到，其中q＝1...5...288，M＝9。
[0083]
[0084] 标准差(σ1,σ2,...,σq)可通过下式计算得到：
[0085]
[0086] 步骤5.异常侦测
[0087] 对4月3日10个监测点的压力观察值，根据下述异常事件判定规则：
[0088] a.任一时间点压力观测值低于-4σ边界；
[0089] b.连续2个时间点压力观测值低于-3σ边界；
[0090] c.连续4个时间点压力观测值低于-2σ边界；
[0091] d.连续8个时间点压力观测值低于-1σ边界。
[0092] 此处以2号监测点为例，侦测到的异常事件如表2，吻合5次模拟爆管试验情况。
[0093] 表2：2号监测点所侦测的异常事件
[0094]序号 异常时间段 持续时间(min) 最大误差(σ)
1 9:25 <5 -5.946
2 9:40 <5 -4.242
3 10:05 <5 -5.339
4 12:25 <5 -9.326
5 19:50 <5 -4.493
