CN101807045B - 基于数据的城市污水泵站系统建模方法 - Google Patents

Abstract

本发明涉及一种基于数据的城市污水泵站系统建模方法。现有排水管网泵站模型建模复杂度高。本发明方法首先建立排水管网泵站预测模型；其次SCADA系统采集的原始数据进行预处理；然后确定模型集类型，辨识模型结构和模型参数辨识；最后对模型参数进行在线动态修正。本发明方法突破传统建模需对机理特性完全了解的束缚，建模时间相对较短，获得模型结构较为简单。

Description

基于数据的城市污水泵站系统建模方法

技术领域

本发明属于自动化技术领域，具体涉及一种基于数据的城市污水泵站系统建模方法。

背景技术

随着城市的飞速发展，城市排水已成为制约城市快速发展的瓶颈之。城市排水主要针对城市现有排水设施，通过排水系统管网泵站模型和实时监测数据，整体优化排水系统运行，提高泵站网络的综合运行效率、节能降耗，提高城市的预洪、防洪能力。

排水管网泵站系统模型是城市排水泵站规划设计、运行控制、维护管理的核心。受天气、居民生活污水排放及管道淤堵、渗漏等复杂因素影响，泵站污水流量具有不确定性、非线性和滞后性。现有排水管网泵站模型(SWMM，STORM和DHI等)多采用基于圣维南方程的水力学建模，建模复杂度高，往往由于缺乏准确的边界条件信息而难以解析求解。同时排水管网系统存在的大量不确定因素，如降雨分布、生活污水排放、管道淤堵、渗漏等，都将增加建模复杂性。

发明内容

本发明的目标是针对现有排水泵站系统建模难题，提供一种基于数据的城市污水泵站系统建模方法，避免传统基于圣维南方程的水力学建模的不足，保证模型精度的同时，也保证模型结构简单并满足实际排水系统运行要求。

本发明首先利用污水泵站系统已知的先验知识，确定系统模型的类型结构和部分参数取值范围；然后根据泵站已有的大量运行数据认知模型参数，获得可行的数据模型。所建模型是实际系统的一个近似，结构简单，建模时间相对较短，适用于城市排水管网泵站实时优化运行和控制。

本发明利用数据采集、数据挖掘、系统辨识、预测机理、数据驱动和数值计算等手段，对难以建立机理模型的情况下，建立基于数据的城市排水管网系统污水泵站模型，以实现对排水运行过程的优化控制和调度。

本发明方法的具体步骤：

步骤(1)建立排水管网泵站预测模型。

现有城市排水系统数据采集与监控系统(SCADA系统)可采集泵站的水位信息，但因测量难题无法获得管道污水流量数据，为预测污水在排水管网系统中各泵站的分布情况，改进Gelormino和Ricker(1994)提出的管网模型，分别建立排水管网泵站污水储量(水位)和污水入流量的预测模型：

x(k)＝x(k-1)+T(z^-1)v(k-1)-s(k-1)+u(k-1)+ε(k-1)

和x′(k-1)＝x(k)-x(k-1)+s(k-1)＝T(z^-1)v(k-1)+u(k-1)+ε(k-1)

其中：k表示采样时刻，x(k)为k时刻该泵站污水储水量，可以通过泵站前池截面积乘以泵站前池水位求得；x′(k-1)为采样周期[k-1 k]内该泵站的污水入流量；s(k-1)为采样周期[k-1 k]内该泵站排出的污水量；u(k-1)和ε(k-1)分别为采样周期[k-1 k]内本地流入该泵站的污水量之和的均值和噪声；v(k-1)为采样周期[k-1 k]内上游泵站排出的污水量；由于上游泵站排出污水量经管网流入下游泵站的过程是非线性和滞后过程，采用T(z^-1)v(k-1)表示上游泵站排出的污水量经管网延迟流入该泵站的污水量，描述为：

$T\left(z^{-1}\right)v(k-1)=\underset{i=n_a}{\overset{n_b}{\mathrm{\Σ}}}{a}_i{z}^{-i}=(a_{n_a}{z}^{-n_a}+\·\·\·+a_{n_b}{z}^{-n_b})v(k-1)$

其中T(z^-1)为不确定延迟算子，z^-1为一步延迟算子，表示延时一步采样周期的相应量；n_a和n_b分别表示上游泵站排水量v对流入该泵站的污水量x′的最小和最大迟延，n_b＞n_a；系数a_i(i＝n_a，…，n_b)表示采样周期[k-i-1 k-i]内上游泵站排出量流入该泵站的污水量的比例。

由于污水在管网流动的延迟和不确定特性，不确定延迟算子T(z^-1)的结构(n_a和n_b)和参数a_i(i＝n_a，…，n_b)未知，后面将利用SCADA系统采集的运行数据辨识不确定算子T(z^-1)，建立数据模型。

步骤(2)数据预处理。

SCADA系统采集的原始数据可能存在噪声、数据不完整甚至是不一致，在利用这些数据进行分析建模前，需要对数据进行预处理。主要包括：

遗漏数据处理：SCADA采集的是时间序列数据，每隔10～30秒采样1次。针对可能存在的遗漏项，首先粗选时间序列数据，每分钟选1个数据；再利用忽略元组或历史数据补全方法处理遗漏数据。

噪声数据处理：采用分箱(按周围值的平均值平滑)方法处理噪声数据。

步骤(3)系统辨识。

由步骤(1)给出的排水管网泵站预测模型可知，由于上游泵站排水量大小对下游泵站的影响不同，不确定延迟算子T(z^-1)中模型结构(n_a和n_b)和模型参数(系数…，

)发生变化。辨识不确定延迟算子T(z^-1)具体步骤如下

a.确定模型集类型

根据先验信息分析，确定影响上游泵站排水量大小的两个具体因素：①开泵台数；②开泵时间。采用开泵率综合评估影响上游泵站排水量大小：开泵率(即计量时间内上游泵站所开泵运行总时间与计量时间内上游泵站可用泵总可运行时间的比值)。根据泵站大量历史运行数据统计分析，计量时间(如一天)内开泵率情况可聚类为3～5类模式，每一类对应建立一个预测模型。实际运行过程中，跟踪开泵率指标，动态判断当前运行时刻属于何种模式，从而切换至合适的模型。

b.模型结构辨识

对流入泵站的污水量x′与上游泵站排水量v进行互相关性分析，从而确定模型结构，即不确定延迟算子T(z^-1)中n_a和n_b值。

n组采样数据构成的时间序列[x′_t]和[v_t](t＝1，2，…，n)与时延q相关，其中1≤q≤n，互相关系数r_q(x′，v)为：

$r_q(x^{\′},v)=\frac{\underset{t=1}{\overset{n-q}{\mathrm{\Σ}}}(x_t^{\′}-{\stackrel{\‾}{x}}^{\′})(v_{t+q}-\stackrel{\‾}{v})}{{\left[\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{(x_t^{\′}-{\stackrel{\‾}{x}}^{\′})}^2\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{(v_t-\stackrel{\‾}{v})}^2\right]}^{1/2}},$ 式中： ${\stackrel{\‾}{x}}^{\′}=\frac{1}{n}\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{x}_t^{\′},$ $\stackrel{\‾}{v}=\frac{1}{n}\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{v}_t$

计算不同时延q下相关系数r_q(x′，v)，相关系数大于置信上限值的时刻就是模型结构n_a和n_b，从而确定各模式下预测模型的结构。

c.模型参数辨识

确定模型结构后，需对模型参数进行辨识，确定的值，采用最小二乘法进行参数辨识。

为便于使用最小二乘法，将步骤(1)得到的模型转化为：

z(k)＝x(k)-x(k-1)+s(k-1)-u(k-1)＝T(z^-1)v(k-1)+ε(k-1)

将上式化为一个标准的最小二乘形式：

z(k)＝h^T(k)θ+ε(k-1)

式中，系数

h(k)＝[v(k-n_a-1)，…，v(k-n_b-1)]^T

将上式可写成矩阵形式：

Z_L＝H_Lθ+∑_L

式中Z_L＝[z(1)，z(2)，…，z(L)]^T，∑_L＝[ε(1)，ε(2)，…，ε(L)]^T，L为辨识数据维数；

最小二乘法的准则函数J(θ)取为：

J(θ)＝(Z_L-H_Lθ)^TΛ_L(Z_L-H_Lθ)

其中，Λ_L为加权阵，一般取为正定对角阵。

极小化J(θ)，可求出系数θ的估计值

$\widehat{\mathrm{\θ}}={\left(H_L^T{H}_L\right)}^{-1}{H}_L^T{Z}_L$

步骤(4)模型参数在线动态修正。

随着观测数据的增多，常规最小二乘法计算量增加，不能及时辨识模型参数。本发明采用遗忘因子最小二乘法，克服数据饱和现象，减弱旧观测数据对未知参数估计的影响，动态修正模型参数。

选择准则函数为：

J_L+1＝αJ_L+[z(L+1)-h^T(L+1)θ_L]²

给出递推公式求解模型系数θ的估计值

$\left\{\begin{array}{c}{\widehat{\mathrm{\θ}}}_{L+1}={\widehat{\mathrm{\θ}}}_L+K(L+1)[z(L+1)-h^T(L+1){\widehat{\mathrm{\θ}}}_L]\\ K(L+1)=\frac{P\left(L\right)h(L+1)}{\mathrm{\α}+h^T(L+1)P\left(L\right)h(L+1)}\\ P(L+1)=\frac{1}{\mathrm{\α}}[1-K(L+1)h^T(L+1)]P\left(L\right)\end{array}\right.$

其中α为遗忘因子，一般取0.95＜α＜1；K(L+1)为权因子；P(L+1)为正定的协方差阵。

本发明方法的有益效果：

(1)该方法突破传统建模需对机理特性完全了解的束缚，建模时间相对较短，获得模型结构较为简单，是实际系统的一个合理近似。

(2)该方法避免求解传统水力学数学模型-圣维南方程组，同时也避免计算机模拟模型对边界参数等基础数据的需求，只需利用大量的泵站历史运行数据，实用更强。

(3)该方法通过模型结构聚类，建立不同模式下预测模型，可更好地模拟预测各种条件下的排水情况。

(4)该方法通过遗忘因子最小二乘法对模型参数进行在线辨识，克服了数据饱和现象，减弱旧观测数据对未知参数估计的影响，能及时反映模型参数的变化，从而可以减少降雨分布、生活污水排放、管道淤堵、渗漏诸多不确定因素对系统性能的影响，提高系统的鲁棒性和可靠性，弥补了传统建模方法的不足。

(5)该方法建立的基于数据的城市污水泵站系统模型，可以很好的预测污水泵站系统的多个参数(污水流入量、污水储量、污水水位)，用于指导泵站的实际运行：(a)合理设置泵站的控制水位，以避免不必要的频繁开关泵，从而降低能耗，延长水泵的寿命；(b)实现预警判断，为保证安全裕度提前开泵，甚至排空泵站污水，以避免污水溢出，实现安全生产。

具体实施方式

以杭州留下和磁带厂两个上下游泵站为例，对本发明作进一步说明。

(1)建立排水管网泵站预测模型；排水管网泵站预测模型包括排水管网泵站污水储量预测模型和污水入流量的预测模型，分别为：

x(k)＝x(k-1)+T(z^-1)v(k-1)-s(k-1)+u(k-1)+ε(k-1)

和x′(k-1)＝x(k)-x(k-1)+s(k-1)＝T(z^-1)v(k-1)+u(k-1)+ε(k-1)

其中：k为采样时刻，x(k)为k时刻磁带厂泵站污水储水量，x′(k-1)为采样周期[k-1 k]内磁带厂泵站的污水入流量；s(k-1)为采样周期[k-1 k]内磁带厂泵站排出的污水量；u(k-1)为采样周期[k-1 k]内本地流入磁带厂泵站的污水量之和的均值；ε(k-1)为采样周期[k-1 k]内本地流入磁带厂泵站的污水量之和的噪声；v(k-1)为采样周期[k-1 k]内留下泵站排出的污水量，T(z^-1)v(k-1)为留下泵站排出的污水量经管网延迟流入该泵站的污水量：

$T\left(z^{-1}\right)v(k-1)=\underset{i=n_a}{\overset{n_b}{\mathrm{\Σ}}}{a}_i{z}^{-i}=(a_{n_a}{z}^{-n_a}+\·\·\·+a_{n_b}{z}^{-n_b})v(k-1)$

其中T(z^-1)为不确定延迟算子，z^-1为一步延迟算子，表示延时一步采样周期的相应量；n_a和n_b分别表示留下泵站排水量v对流入磁带厂泵站的污水量x′的最小迟延和最大迟延，n_b＞n_a；系数a_i(i＝n_a，…，n_b)为采样周期[k-i-1 k-i]内上游泵站排出量流入该泵站污水量的比例；

(2)数据预处理：对SCADA采集的数据进行数据遗漏处理和噪声处理。

遗漏数据处理：泵站SCADA每隔20秒采样1次数据，首先从采集的时间序列数据中进行粗选，每分钟选1个数据，根据前后数据补全，关联数据补全方法进行遗漏数据处理；对于连续缺失较多的数据，通过对比前一天的数据，进行类比补全。

噪声数据处理：对开关泵时水面波动造成的水位测量误差，通过多传感器测量值平均和前后数据平滑校正进行处理；对于明显的偏大偏小的数据，通过去除毛刺方法进行校正。

(3)数据聚类：根据预处理后的历史数据计算统计指标值-开泵率，通过统计分析将开泵率聚类为3～4类模式，每一类对应建立一个预测模型。实际运行过程中，跟踪开泵率指标，动态判断当前运行模式，切换至合适的预测模型。

(4)模型结构辨识：分别对每种模式的杭州磁带厂泵站污水流入量x′与其上游留下泵站的排水量v进行相关性分析，获得每类模式下留下泵站和磁带厂泵站互相关系数值r_q(x′，v)，污水泵站液位控制的运行机理决定了互相关系数的周期性。

根据经验公式计算得到的流速和上下游泵站的距离计算污水流经管网时间，确定正相关性对应的周期，即上游留下泵站的排水量影响下游磁带厂泵站流入量的时间段。该周期内的相关系数大于置信上限值时所对应的最小时刻和最大时刻就是模型中不确定延迟算子T(z^-1)中n_a和n_b值，从而确定各模式下预测模型的结构。

(5)模型参数辨识：采用最小二乘法对各种模式下预测模型进行参数辨识，即确定

的值。

将预测模型转化为如下形式：

z(k)＝x(k)-x(k-1)+s(k-1)-u(k-1)＝T(z^-1)v(k-1)+ε(k-1)

将上式转化为一个标准的最小二乘形式：

z(k)＝h^T(k)θ+ε(k-1)

式中，系数

h(k)＝[v(k-n_a-1)，…，v(k-n_b-1)]^T。

将上式可写成矩阵形式：

Z_L＝H_Lθ+∑_L

式中Z_L＝[z(1)，z(2)，…，z(L)]^T，∑_L＝[ε(1)，ε(2)，…，ε(L)]^T，L为辨识数据维数；

根据从而得到模型参数的估计值

(6)模型参数在线动态修正。

通过步骤(2)至步骤(4)得到不同模式下的预测模型，在实际运行时，为了能及时反映不确定因素对模型参数的影响，采用遗忘因子最小二乘法对模型参数进行在线修正。

取遗忘因子α值为0.98。初值θ和P取步骤(4)中参数辨识的结果，即

误差ε＝0.01。按如下递归方程进行在线参数修正：

$\left\{\begin{array}{c}{\widehat{\mathrm{\θ}}}_{L+1}={\widehat{\mathrm{\θ}}}_L+K(L+1)[z(L+1)-h^T(L+1){\widehat{\mathrm{\θ}}}_L]\\ K(L+1)=\frac{P\left(L\right)h(L+1)}{\mathrm{\α}+h^T(L+1)P\left(L\right)h(L+1)}\\ P(L+1)=\frac{1}{\mathrm{\α}}[1-K(L+1)h^T(L+1)]P\left(L\right)\end{array}\right.$

(7)排水管网泵站污水分布预测。

由步骤(3)至步骤(5)确定各模式下泵站预测模型，利用该模型可以预测下游泵站的污水流入量和污水储量，以及泵站前池水位。

Claims (1)

1.基于数据的城市污水泵站系统建模方法，其特征在于该方法包括如下步骤：

步骤(1)建立排水管网泵站预测模型；排水管网泵站预测模型包括排水管网泵站污水储量预测模型和污水入流量的预测模型，分别为：

x(k)＝x(k-1)+T(z^-1)v(k-1)-s(k-1)+u(k-1)+ε(k-1)

和x′(k-1)＝x(k)-x(k-1)+s(k-1)＝T(z^-1)v(k-1)+u(k-1)+ε(k-1)

其中：k为采样时刻，x(k)为k时刻该泵站污水储水量，x′(k-1)为采样周期[k-1 k]内该泵站的污水入流量；s(k-1)为采样周期[k-1 k]内该泵站排出的污水量；u(k-1)为采样周期[k-1 k]内本地流入该泵站的污水量之和的均值；ε(k-1)为采样周期[k-1 k]内本地流入该泵站的污水量之和的噪声；v(k-1)为采样周期[k-1 k]内上游泵站排出的污水量，T(z^-1)v(k-1)为上游泵站排出的污水量经管网延迟流入该泵站的污水量：

$T\left(z^{-1}\right)v(k-1)=\underset{i=n_a}{\overset{n_b}{\mathrm{\Σ}}}{a}_i{z}^{-i}=({a_n}_a{z}^{{-n}_a}+...+{a_n}_b{z}^{{-n}_b})v(k-1)$

其中T(z^-1)为不确定延迟算子，z^-1为一步延迟算子，表示延时一步采样周期的相应量；n_a和n_b分别表示上游泵站排水量v对流入该泵站的污水量x′的最小迟延和最大迟延，n_b＞n_a；系数a_i(i＝n_a，…，n_b)为采样周期[k-i-1 k-i]内上游泵站排出量流入该泵站污水量的比例；

步骤(2)对数据进行预处理，具体包括处理遗漏数据和处理噪声数据；

处理遗漏数据的方法为：对数据采集与监控系统采集的时间序列数据首先进行数据粗选，即每分钟选1个数据；然后再利用忽略元组或历史数据补全方法处理遗漏数据；

处理噪声数据的方法为：采用分箱方法处理噪声数据；

步骤(3)辨识不确定延迟算子T(z^-1)，具体方法是：

a、确定模型集类型；

根据泵站大量历史运行数据将上游泵站开泵率聚类为3～5类模式，每一类对应建立一个预测模型，跟踪当前开泵率指标，动态判断当前运行时刻属于何种模式，从而切换至与该模式对应的预测模型；

所述的开泵率为计量时间内上游泵站所开泵运行总时间与计量时间内上游泵站可用泵总可运行时间的比值；

b、对流入该泵站的污水量x′与上游泵站排水量v进行互相关性分析，确定模型结构n_a和n_b；具体方法是：

n组采样数据构成时间序列[x′_t]和[v_t](t＝1，2，…，n)与时延q相关，其中1≤q≤n，互相关系数r_q(x′，v)为：

$r_q(x^{\′},v)=\frac{\underset{t=1}{\overset{n-q}{\mathrm{\Σ}}}(x_t^{\′}-{\stackrel{\‾}{x}}^{\′})(v_{t+q}-\stackrel{\‾}{v})}{{\left[\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{(x_t^{\′}-{\stackrel{\‾}{x}}^{\′})}^2\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{(v_t-\stackrel{\‾}{v})}^2\right]}^{1/2}},$ 式中： ${\stackrel{\‾}{x}}^{\′}=\frac{1}{n}\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{x}_t^{\′},$ $\stackrel{\‾}{v}=\frac{1}{n}\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{v}_t$

计算不同时延q的相关系数r_q(x′，v)，当r_q(x′，v)大于设定的置信上限值时，r_q(x′，v)所对应的最小时刻和最大时刻就是模型结构n_a和n_b；

c、采用最小二乘法确定预测模型参数

将步骤(1)得到的预测模型转化为：

z(k)＝x(k)-x(k-1)+s(k-1)-u(k-1)＝T(z^-1)v(k-1)+ε(k-1)

将上式化为一个标准的最小二乘形式：

z(k)＝h^T(k)θ+ε(k-1)

式中，系数

h(k)＝[v(k-n_a-1)，…，v(k-n_b-1)]^T；

将上式写成矩阵形式：

Z_L＝H_Lθ+∑_L

式中Z_L＝[z(1)，z(2)，…，z(L)]^T，∑_L＝[ε(1)，ε(2)，…，ε(L)]^T，L为辨识数据维数；

最小二乘法的准则函数J(θ)取为：

J(θ)＝(Z_L-H_Lθ)^TΛ_L(Z_L-H_Lθ)

其中，Λ_L为加权阵，取正定对角阵；

极小化J(θ)，得到系数θ的估计值

d、在线动态修正模型参数

选择准则函数J_L+1，J_L+1＝αJ_L+[z(L+1)-h^T(L+1)θ_L]²；

根据递推公式求解系数θ的估计值

$\left\{\begin{array}{c}{\widehat{\mathrm{\θ}}}_{L+1}={\widehat{\mathrm{\θ}}}_L+K(L+1)[z(L+1)-h^T(L+1){\widehat{\mathrm{\θ}}}_L]\\ K(L+1)=\frac{P\left(L\right)h(L+1)}{\mathrm{\α}+h^T(L+1)P\left(L\right)h(L+1)}\\ P(L+1)=\frac{1}{\mathrm{\α}}[1-K(L+1)h^T(L+1)]P\left(L\right)\end{array}\right.$

其中α为遗忘因子，一般取0.95＜α＜1；K(L+1)为权因子；P(L+1)为正定的协方差阵。