CN113343595A - 一种明渠输水系统事故的反演模型及事故流量、位置确定方法 - Google Patents

Abstract

本发明公开了一种明渠输水系统事故的反演模型及事故流量、位置确定方法，本发明在ID模型的基础上加入正常取水流量和事故流量，建立反演线性模型，反演线性模型包括水面面积和滞后时间两个参数，通过参数辨识识别水面面积和滞后时间，然后通过辨识得到的滞后时间计算事故影响滞后时间，结合上、下游水位变化的时间差确定事故位置，通过事故影响滞后时间和水面面积代入反演线性模型中计算事故流量。本发明能得到较为准确的事故源信息，对快速启动应急预案，降低事故不利影响，保障渠系或者引水工程本身结构安全有着十分重要的意义；建立的事故反演模型机理性强，所需算力较小，有利于提高工作效率，提高准确性、实时性。

Description

一种明渠输水系统事故的反演模型及事故流量、位置确定 方法

技术领域

本发明属于水利领域，涉及一种明渠输水系统事故识别技术，具体涉及一种明渠输水系统事故的反演模型及事故流量、位置确定方法。

背景技术

水资源是非常重要的自然资源，也是战略性经济资源，是一个国家综合国力的有机组成部分。我国的水资源问题是非常突出的，我国人均水资源量2200m³，目前有16个省(市、区)人均水资源低于严重缺水标准，有6个省、区人均水资源量低于500m³。同时我国水资源还呈现出①年内年际分配不均，旱涝灾害频发，供需矛盾突出；②地区分布不均，水土资源不匹配；③水资源利用率低，污染严重等特点。

针对我国水资源时空分布不均的问题，我国投入大量资金，众多水利工作者贡献、付出了诸多智慧与汗水，为我国设计、建设了诸多用以灌溉的灌溉渠系，和用以缓解地区水资源空间分配不均的长距离引水工程。我国灌区普查资料显示，全国设计灌溉面积30万亩及以上灌区456处，设计灌溉面积达2.80亿亩；设计灌溉面积1万～30万亩灌区7316处，设计灌溉面积2.33亿亩；50～1万亩灌区205.82万处，设计灌溉面积3.42亿亩，总体规模是十分庞大的。目前我国已经建成一大批以解决城市供水供需矛盾的大型引水工程，如南水北调工程、引滦入津、引青济秦、引黄入晋、东深供水等工程，其长度大都在100～200km，其中南水北调中线更是突破了1200km。

诸如此类灌溉渠系以及引水工程都是以明渠为主作为输水建筑物的，将其定义为明渠输水系统。此类明渠输水系统往往会由于运行管理不善，或者特殊地质条件或者不良气象条件，引发突发事故，对系统平稳运行带来巨大挑战，甚至引发较大的次生灾害，造成较大的经济损失。以“引江济汉"工程为例，推荐线-龙高Ⅰ线有70％的渠道经过膨胀土分布区，在渠道开挖和运行过程中，其边坡可能产生胀缩变形破坏，进而引发失稳事故，再如南水北调工程，总干渠焦作至沧河，长约100km的区段内及南水北调中线干线漕河段岗头隧洞曾发生过泥石流，南水北调工程中线总干渠高填方渠道累计48段，总长度53.3km，此类高填方渠道在极端工况或暴雨情况下较易发生决口事故。由于我国灌区分布极为广泛，因此由于部分灌区较差地质条件，管涌等渗透破坏事故也时有发生。因为闸门的操作不当或者降雨汇流入渠也可能引发渠道漫溢事故。

如上所述，突发事故对于明渠输水系统是不可避免的，而其发生之后往往会造成较大的人民生命财产损失，因此研究明渠输水系统在事故工况下调度是非常有必要的，目前对明渠事故的研究主要集中在事故应急调度及水力响应特性方面，以及污染物的溯源问题。事故应急调度预案主要是指对包括闸门、退水闸、溢流堰等建筑物的联合调度预案，以期望将损失降低到最低。由于针对明渠事故反演的定义尚不明确，因此这里对事故反演这一概念进行限定和说明。事故反演是一类典型的反问题，指的是基于明渠系统测点数据利用反演模型，对事故具体点位，事故发展过程进行较为准确的计算，难度更大，精度也更高。事故反演是针对性地启动应急预案的必要条件，因为只有知道事故类型、事故渠池、事故具体点位、事故严重程度才能对症下药地提出适当的应急预案，事故反演是事故应急调度的基础和前提条件，但是目前事故准确的具体信息的获取是非常困难的，往往是在事故持续发展一段时间后，引起系统较大的水位流量波动被管理人员观测到，然后进行现场巡渠才能确定，尤其是当事故发生在夜间，响应周期会更久。这样一来就浪费了大量的时间，错过了通过事故调度减小事故损失的最佳时机，增大了事故继续发展为其他更加恶性的事故或者引发次生灾害的风险。因此建立在线的明渠输水系统事故反演系统对快速启动应急预案，降低事故不利影响，保障渠系或者引水工程本身结构安全有着十分重要的意义。其次将事故反演纳入渠道控制系统，有利于提高工作效率，提高准确性、实时性，这也是我国智慧水网建设、水利数字化建设的必然要求。

发明内容

本发明要解决的技术问题在于针对明渠输水系统输水调度技术缺陷，提出一种明渠系统的事故反演模型。提出的事故反演模型将明渠系统实测得到的各个渠池下游水位、各个渠池的入流流量作为模型输入量，结合该反演模型可以反演求得明渠的事故流量过程和事故具体点位，其中事故具体点位是指反演得到的包含事故具体点的渠段范围，而事故流量过程指的是明渠系统因事故出现的流量出流或者入流过程，如发生管涌之后的流量泄露和损失过程、发生可能引发渠道漫溢事故的降雨汇流入渠流量过程。这一模型可以用于对明渠输水过程中泄露等水量损失进行监测，同时反演得到的事故点位和事故流量过程等较为准确的事故源信息将为应急调度提供依据，有利于快速启动应急预案，降低事故不利影响，防止事故进一步蔓延或引发次生灾害，有利于保障高效供水和渠系或者引水工程本身结构安全。

为了解决上述技术问题，本发明采用如下技术方案：

一种明渠输水系统事故的反演模型，该模型基于ID模型改进，ID模型将整个渠池概化为上游的均匀流区和下游的回水区，并认为均匀流区的存在将会使下游水位响应滞后于上游流量变化，而处于回水区段的渠池不存在滞后时间，模型包含两个基本参数：水面面积A_S和入流流量滞后时间τ，其特征在于：改进的方式为在ID模型的基础上加入正常取水流量q_off和事故流量q_ac，改进后的ID模型如下：

对公式(2)进行变形得到明渠事故的反演线性模型

上式中：A_S为回水区水面面积(m²)；y为下游水位相对于初值的增量(m)；x为回水区向下游方向的某处位置(m)；q_in为渠池上游入流流量相对于初值的变化量(m³/s)；q_out为渠池下游出流流量相对于初值的变化量(m³/s)；t为时间(s)；τ为入流流量滞后时间，即上游流量变化经均匀流区传递到下游回水区所用的时间(s)；

正常取水流量q_off为取水口的取水量的变化量，在事故反演计算过程中，时间相对较短，取水量近似认为无变化，因此，q_off在事故反演过程中认为为0，q_ac为事故流量相对于初始状态的增量，但初始状态一般认为明渠系统正常运行，故此该值就是事故流量本身；

对公式(3)进行离散化，得到离散后的反演线性模型如下：

DT为计算的时间间隔，k表示计算次数，即第k步，τ_d″为离散公式下的事故流量下游回水区水位影响的滞后时间序数，即事故影响滞后时间τ″含有DT的个数，由τ″/DT之后四舍五入得到；τ_d为离散公式下的上游流量下游回水区水位影响的滞后时间序数，即滞后时间含有DT的个数，由τ/DT之后四舍五入得到；i为渠池编号，q_i,in为渠池i的入流量变化量，通过实测获取数据；q_i,out为渠池i的出流量变化量，通过实测获取数据；q_i,ac为渠池i的事故流量，Z_i(k)为第k步，渠池i的下游水位，实测获取；

事故影响滞后时间滞后时间τ″通过公式

计算获得，l_ac为事故具体点位距离渠池入口的距离，(m)；l_tal为事故渠池的总长度，入流流量滞后时间τ和回水区水面面积A_S是模型参数，通过系统辨识的方法求得。

本发明还提供一种利用上述反演模型计算事故渠池的事故流量及事故位置的方法，其特征在于，包括以下步骤：

步骤1、通过圣维南方程组对需要辨识的渠池建立仿真模型；

步骤2、给定渠池渠池上游入流流量的变化量，通过步骤1中的仿真模型用直线拟合渠池下游水位上升过程，采用最小二乘法求取回水区水面面积A_S，通过渠池下游水位变化起始点获取该渠池的入流流量滞后时间τ；

步骤3、通过公式

计算得到事故具体点位置，l_ac为事故点位距离事故渠池进口的距离，Δτ表示上游开始监测到水位变化到下游开始监测到水位变化的时间差，通过对渠池实际观测得到，L为渠池实际总长，之后根据公式

计算得到事故影响滞后时间τ″；

步骤4、将步骤2和步骤3中得到的入流流量滞后时间τ、事故影响滞后时间和回水区水面面积A_S代入反演线性模型中即可得到事故渠池i的事故流量q_i,ac，完成事故流量及事故位置的确定。

进一步地，步骤2中，直线拟合渠池下游水位上升过程方法具体采用最小二乘法拟合，使目标函数

即最小二乘估计值，目标函数如下式所示：

上式中：y(k)为辨识渠池下游水位的实测值；u(k)为流量变化量u(k)＝q_in(k-τ)-q_out(k)，a为y(k-1)的线性系数，理论值为-1；b为u(k)前的线性系数，值为DT/A_S，通过对公式(5)的拟合求出b，通过b计算得到渠池下游回水区面积A_S。

本发明产生的有益效果：

(1)尚未有应用于明渠输水系统事故反演模型，提出的反演模型基于改进后的ID模型建立，对应用于明渠事故反演模型的建立具有重要参考意义；

(2)建立的事故反演模型能实现明渠事故具体点位和事故流量的反演，能得到较为准确的事故源信息，对快速启动应急预案，降低事故不利影响，保障渠系或者引水工程本身结构安全有着十分重要的意义；

(3)建立的反演模型也可用于明渠泄露监测，可以监测由于闸门泄露或结构破坏导致的流量损失，以及时进行监测维修，避免大量水资源的进一步浪费；

(4)建立的事故反演模型机理性强，所需算力较小(计算时长少于30s)，有利于提高工作效率，提高准确性、实时性。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为现有技术中ID模型基本原理图，是建立事故反演模型的基础，其中图1(a)为ID模型渠道概化示意图，图1(b)为渠池下游水位响应过程。

图2是改进的ID模型的原理图。

图3是确定事故点位的基本原理图。

图4为小流量类型事故的反演实施例。

图5是给出的小流量事故工况下反演模型的反演解效果图。

图6未修正模型反演解图。

图7决口工况2事故下游水位过程线。

图8大流量事故反演修正模型反演解图。

图9为ID模型参数辨识边界设置图。

图10为ID模型参数辨识下游水位响应过程。

具体实施方式

下面结合附图和实施例对本发明的实施方式作进一步详细描述。以下实施例用于说明本发明，但不能用来限制本发明的范围。

本发明一种明渠输水系统事故的反演模型，该模型基于ID模型改进，ID模型将整个渠池概化为上游的均匀流区和下游的回水区，并认为均匀流区的存在将会使下游水位响应滞后于上游流量变化，而处于回水区段的渠池不存在滞后时间，模型包含两个基本参数：水面面积AS和入流流量滞后时间τ，改进的方式为在ID模型的基础上加入正常取水流量q_off和事故流量q_ac。

一、反演模型的建立过程

(一)原模型

要建立的反演模型需要解决通过实测的渠池上下游水位和流量数据反推渠池内某点事故源基本信息(位置和流量过程)的问题，其内核在于建立可以描述渠池水力特性的水位流量关系。而ID(Integral Delay，积分滞后)模型忽略了次要的水力要素，建立了水位-流量线性模型。同时具有物理概念明确、结构简单、计算存储量小等优点，且目前已经在渠道控制器设计中得到了较为广泛的应用，故考虑基于该模型开发事故反演模型。

如图1所示的ID模型将整个渠池概化为上游的均匀流区和下游的回水区，其中均匀流区水深为正常水深，并认为均匀流区的存在将会使下游水位响应滞后于上游流量变化，而处于回水区段的渠池不存在滞后时间，模型包含两个基本参数：水面面积A_S和入流流量滞后时间τ，其表达式如公式(1)所示。如图1(b)所示，该模型反映了渠池下游水位对上游流量变化的响应过程。若入流流量从t₁时刻开始增大、而出流流量保持不变，当到达t₂时刻时，下游水位开始线性增大，斜率为ΔQ/A_S，入流流量滞后时间为τ＝t₂-t₁。

公式(1)中：A_S为回水区水面面积(m²)；y为下游水位相对于初值的增量(m)；x为回水区向下游方向的某处位置(m)；q_in为渠池上游入流流量相对于初值的变化量(m³/s)；q_out为渠池下游出流流量相对于初值的变化量(m³/s)；t为时间(s)；τ为入流流量滞后时间，即上游流量变化经均匀流区传递到下游回水区所用的时间(s)。

(二)改进后的ID模型

上述ID模型考虑了渠池进出口流量变化对回水区水位的影响，但无法解决事故流量发展过程对下游水位测点的影响，因此需要对ID模型进行改进，增加事故流量项和取水流量项。新模型的示意图如图2所示，其物理意义与原模型一致，但加入两项：正常取水流量(q_off)；事故流量(q_ac)。改进后模型则表征在上下游流量以及正常取水以及事故流量的综合影响下的回水区水位的动态响应规律，具体公式如公式(2)所示：

式中：A_S、y、t、τ、q_in、q_out表示的物理含义与公式(1)意义相同；q_off为取水口的取水量的变化量，为了简化问题，避免引入过多干扰，在事故反演过程中认为q_off为0，(单位m³/s)；τ′为取水口相对于回水区的滞后时间，(单位s)；q_ac为事故流量相对于初始状态的增量，但初始状态一般认为明渠系统正常运行，故此该值就是事故流量本身单位(m³/s)；τ″表示事故流量对下游回水区水位影响的滞后时间，该值的取值与事故具体点位有关，可以近似认为满足公式(3)，(单位s)。

式中：l_ac为事故具体点位距离渠池入口的距离，(m)；l_tal为事故渠池的总长度，(m)。

(三)基于改进的ID模型构建反演模型

改进后的ID模型中：q_in、q_out可以在得知渠池上下游进出流量的基础上求得；y可以通过下游水位测点的数据计算求得；A_S、τ是模型参数，可以通过系统辨识的方法求得。整个模型的未知量为事故流量过程q_ac以及反映事故具体点位的事故流量滞后时间τ″，这些也是模型的输出量。因此只需要对公式(2)进行变形，便可以得到明渠事故的反演线性模型，如公式(4)所示：

对公式(4)进行离散化，得到离散后的反演模型如公式(5)所示：

公式(5)中：i为渠池编号；k为计算步表示第k步；DT表示计算时间间隔，(s)；Z_i(k)为第k步的下游水位值，(m)；q_i,in为渠池i的入流量变化量，(m³/s)；q_i,out为渠池i的出流量变化量，(m³/s)；τ为入流流量滞后时间；τ″为事故流量滞后时间；τ_d″为离散公式下的事故流量下游回水区水位影响的滞后时间序数，即事故影响滞后时间τ″含有DT的个数，由τ″/DT之后四舍五入得到；τ_d为离散公式下的上游流量下游回水区水位影响的滞后时间序数，即滞后时间含有DT的个数，由τ/DT之后四舍五入得到。

本发明还提供一种利用上述反演模型求解计算事故流量和渠池中事故位置的方法

二、阐释反演模型的求解方法

由公式(5)可知反演线性系统待求量有两个(τ″，q_ac)，由于我们假设事故滞后时间τ″与事故点所在位置有关，详见公式(3)。τ″求解问题则转化为位置溯源问题。

该方法是基于亚临界流基本理论提出的，即认为水波向上游的传播速度与向下游的传播速度基本一致。因此当事故发生引起流量变化后，上下游初次监测到水位波动的时间差就反映了事故具体点位。如果渠池总长度用L表示，重力波波速用ω表示，则波从上游传递到下游的总时间为τ＝L/ω。用Δτ表示上游开始监测到水位变化到下游开始监测到水位变化的时间差，如果该值接近于0，则说明事故点靠近渠池中间点位，若Δτ＝±τ，则说明事故点位靠近渠池上下边界(为-Δτ时靠近渠池进口，为τ时靠近渠池出口)。基本原理示意图如图3所示，其具体位置的计算如公式(6)所示：

式中：l_ac为事故点位距离事故渠池进口的距离，m，Δτ表示上游开始监测到水位变化到下游开始监测到水位变化的时间差，通过对渠池实际观测得到，L为渠池实际总长。

求解方法具体如下：

步骤1、通过圣维南方程组对需要辨识的渠池建立辨识模型；

步骤2、给定渠池渠池上游入流流量的变化量，通过步骤1中的辨识模型用直线拟合渠池下游水位上升过程，采用最小二乘法求取回水区水面面积A_S，通过渠池下游水位变化起始点获取该渠池的滞后时间τ；

步骤3、通过公式(6)计算得到事故具体点位置，x_l为事故点位距离事故渠池进口的距离，Δτ表示上游开始监测到水位变化到下游开始监测到水位变化的时间差，通过对渠池实际观测得到，L为渠池实际总长，之后根据公式

计算得到事故影响滞后时间τ″；

步骤4、将步骤2和步骤3中得到的滞后时间τ、事故影响滞后时间τ″和回水区水面面积A_S代入反演线性模型中即可得到事故渠池i的事故流量q_i,ac，完成事故流量及事故位置的确定。

通过滞后时间τ和反应事故具体点位即可通过公式(3)计算事故滞后时间τ″，将事故滞后时间τ″和回水区水面面积A_S代入反演模型求解事故流量q_ac。若该值在0附近，则表示不存在流量的进入或者流出。若该值为正，则表明事故渠池流量的流出，这可能是管涌(人为偷水、泄漏)或者决口、结构破坏引发的。若数值为负，则表示事故渠池存在着流量入渠的情况，可能是降雨汇流导致的。

三、结合实施例说明针对管涌等小流量事故的反演模型精度

如图4所示的实施例表示一个由9个渠池构成的明渠系统，在渠池3中游位置发生了管涌事故，我们假定三种工况，然后通过该反演模型反演该事故点的位置以及流量过程，以展示模型反演精度。将反演结果汇总于下表及图5。

表1管涌事故反演解参数汇总表

由表1可知，事故发生具体点位的反演精度与实测数据的采样时间有关，同时也受到事故流量大小的影响，事故越轻微，这上下游波动差别越小，具体点位的反演难度也较高，如工况1、2的反演区间长达3072m。但总体来说，事故点位反演区间均包含了事故真实发生点位，且能将反演区间长度控制在1500m左右，该长度会随着采样时间间隔变短而变短，但对于工程实际而言，并不需要过分严格的反演精度，相反缩短采样时间会对设备的存储空间提出较高的要求。

对比多种工况发现建立的明渠事故反演线性模型针对小流量工况反演精度较高，其平均偏差小于5％，最大相对偏差不超过20％。同时也应当注意到对于事故流量的转折点，线性模型反演解偏差较大，这也是低阶模型模拟高阶变化不可避免的问题。而远离转折点时，模型反演解逐渐收敛于真实值。同时所有工况反演耗时不超过30s，其计算成本较小，有利于实现事故在线反演预警以及流量泄漏检测。

四、反演模型在大流量事故(决口或结构破坏事故)下的修正

由于大流量变幅情况下，整个明渠系统处于高度非线性，高耦合性的状态，依旧采用线性模型来反演事故流量过程，势必会产生较大误差。因此以某决口事故为例，对未修正的反演模型进行误差分析，之后进行模型修正，构建适用于大流量变幅情况下的明渠系统反演模型。

决口工况下仿真事故流量过程如图6所示，发现反演模型在真实事故流量到达平台期前都有较好反演性能，当时间到达t＝550min后反演解持续偏离真实值。这是线性模型的缺陷造成的，当事故流量维持不变时，整个事故渠池已经处于高度非线性，并由于渠池的蓄量之间的补充，处于耦合性较高的状态，此时用线性模型难以反映这种非线性和耦合性，因此反演解持续偏离真实值。但结合图7，发现t＝550min后的水位降速已经维持不变了，这也反映了事故流量已经和其他渠池的流量补充达到了平衡，在此之后事故流量势必是维持不变的。基于此对大流量情况下的事故反演模型进行修正，认为下游水位降速维持不变时，反演解将维持不变。即认为水位降速基本维持不变时，取后续事故流量值等于该时刻反演得到的事故流量。利用修正后的模型重新进行反演得到反演解如图8所示，发现大流量工况下反演解持续偏离真实值的问题得到了解决，且相对偏差仅为3％，因此可以在此类大流量事故反演时采用修正后的反演模型。

事故反演模型参数A_S和滞后时间τ辨识确定方法如下：

如图9和图10所示，本发明采用的参数确定方法是系统辨识方法，是通过实测数据或者仿真数据，求解满足精度要求的系统参数的一类方法的总称。钟乐等采用非恒定流的仿真计算获得数据，然后利用系统辨识的方法(最小二乘法)获得回水区面积A_S，与滞后时间τ。

具体方法如下：

(1)利用非恒定流仿真程序或者商业软件对需要辨识的渠池建模(采用圣维南方程组建模)。

(2)将渠池进出口流量均设为Q₀，即入流等于出流。获得较为稳定的初始状态。

(3)在t₀时刻在维持下游流量Q₀不变的情况下，给上游来流增加ΔQ。使得渠池水位上涨，持续运行一段时间，获取下游水位变化过程数据。通过计算下游水位起涨时间和上流流量阶跃增加时间的差值来获取滞后时间τ。

(4)通过仿真数据，用直线拟合水位上升过程，采用最小二乘法求取A_S。

使目标函数

这里的

也即最小二乘估计值。结合本发明用于辨识反演模型参数的目标函数如下式所示：

公式(7)中：y(k)为辨识渠池下游水位的实测值，(m)；u(k)为流量变化量

a为y(k-1)的线性系数，理论值为-1；b为u(k)前的线性系数，值为DT/A_S。由此可知通过系统辨识可获得回水区面积A_S。

本发明对ID模型进行了改进，并基于此构建了明去事故反演模型，能较为精确地反演事故具体点位和事故流量过程。在管涌和明渠泄露等小流量工况下事故流量反演精度较高，但事故点位反演精度有待进一步提高，在决口或结构破坏此类大流量变幅的事故工况下，反演模型需要以水位作为判别条件进行修正，修正后的模型有较高精度，满足工程实际需要。

该模型可用于输水监控系统的实时事故反演和流量监测。同时本文提出的模型方法对于我国大量已建、在建、规划中的大型引、调水工程具有重要的参考价值，事反演的实现有助于工程的长期安全、高效运行，同时对于我国大量的灌区输配水系统自动化、智能化升级亦具有一定的借鉴意义。

需要说明的是，本发明事故渠池确定方法采用本发明一起申请专利一种明渠输水系统事故识别方法中方法确定，具体方式如下：

步骤1、输入渠道建模参数，建立实际渠道的仿真模型，渠道建模参数包括渠池数量、闸门尺寸、渠道坡度、建筑物及糙率，模型以王长德、管光华2011年开发的输水渠道系统运行仿真与控制软件为基础建立。

步骤2、将用水计划作为输入量，通过仿真模型仿真计算正常运行工况下各渠池的闸门目标开度过程、下游目标水位过程以及目标过闸流量过程，得到不同时刻各渠池的闸门目标开度、下游目标水位以及目标过闸流量；

步骤3、读取实测得到各渠池闸门前后的水位、闸门开度以及过闸流量，结合步骤1得到的不同时刻各渠池的闸门目标开度、下游目标水位以及目标过闸流量计算每个渠池的四个判断指标，四个判断指标分别为绝对最大闸门开度偏差、最大水位偏速、末态流量增量和末态水位增量，具体如下：

(1)最大闸门开度偏差ΔABEmax

ΔABEmax指的是所有采样时刻下所有渠池实测的闸门开度与用水计划下的目标开度的差值的最大值的绝对值，反映了事故工况下闸门相对于目标开度的偏离程度，在不考虑反馈的情况下，非闸门故障情况下该值应该维持在零附近。若该值大于2cm则认为发生了闸门事故，最大闸门开度偏差计算公式如下：

ΔABEmax＝max(max(|E_rij-E_tij|))

上式中，E_rij为i渠池j时刻实测的闸门开度(m)；E_tij为用水计划下i渠池j时刻的目标闸门开度(m)；

(2)最大水位偏速ΔZspeedmax

ΔZspeedmax指的是所有采样时刻下所有渠池在单位采样时间步长下，下游测点水位相较于用水计划下游水位的最大偏移量，反映了事故工况下游水位相较于目标水位的偏离程度；最大水位偏速包括正值最大水位偏速ΔZspeedmax₁和负值最大水位偏速ΔZspeedmax₂，计算公式如下：

上式中，DT为计算的时间间隔，y_rij为i渠池在j时刻实测的下游水位(m)；y_tij为用水计划下i渠池在j时刻的下游目标水位(m)；

(3)末态流量增量ΔQend

ΔQend_i是指渠池i末时刻实测的入渠流量，或者通过实测水位，闸门开度间接计算的入渠流量与用水计划下的目标流量的差值，反映了事故工况入渠流量相较于目标流量的偏离程度，计算公式如下：

ΔQend_i＝Q_rend,i-Q_tend,i

上式中，Q_rend,i为渠池i末时刻实测的或利用其他实测数据间接计算出的入渠流量(m³/s)；Q_tend,i为用水计划下末时刻渠池的入渠流量(m³/s)；

(3)末态水位增量ΔZend

ΔZend_i是指渠池i末时刻实测的下游水位与用水计划下的目标水位之间的差值，反映了事故工况下渠池水位相较于目标水位的偏离程度，计算公式如下：

ΔZend_i＝Z_rend,i-Z_tend,i

上式中，Z_rend为渠池i末时刻实测的下游水位(m)，Z_tend,i为用水计划下末时刻的渠池i下游水位(m)；

步骤4、根据四个判断指标判断事故渠池和事故类型，具体方法如下：

步骤4.1、判断最大闸门开度偏差ΔABEmax的大小，如果最大闸门开度偏差大于阈值，则事故类型为闸门故障，通过闸门开度及开度变化，找出事故闸门，对应的渠池为事故渠池，判断结束，如果最大闸门开度偏差小于阈值，则执行下面步骤：

步骤4.2、设定水位波动速度阈值

对最大水位偏速进行判断，如果最大水位偏速满足

且

水位都出现下降，则为流量出流类型事故；同时通过末态流量增量指标判断事故渠池位置，具体的，计算各个渠池的末态流量增量，对每个渠池的末态流量增量，当n渠池上游末态流量增量为正(△Qend>0)，而下游渠池末态流量增量为负(△Qend<0)，则n渠池为事故渠池；否则执行下面步骤；

步骤4.3、如果最大水位偏速满足

且

水位都出现上涨，则为因降雨入流导致的渠道漫溢事故，同时通过末态流量增量指标判断事故渠池位置，具体的，计算各个渠池的末态流量为正(△Qend>0)，当n渠池上游末态流量增量为负(△Qend<0)，而下游渠池末态流量增量为正，则n渠池为事故渠池；否则整个明渠输水系统水位有升有降，(

且

)或者(

且

)，为渠道淤塞事故，计算各个渠池的末态水位增量Z_tend，当n渠池上游末态水位增量为正，上游出现壅水现象(△Zend>0)，而下游渠池末态水位增量为负，则n渠池为事故渠池；完成明渠输水系统事故识别。

以上实施方式仅用于说明本发明，而非对本发明的限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行各种组合、修改或者等同替换，都不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围当中。

Claims (3)

1.一种明渠输水系统事故的反演模型，该模型基于ID模型改进，ID模型将整个渠池概化为上游的均匀流区和下游的回水区，并认为均匀流区的存在将会使下游水位响应滞后于上游流量变化，而处于回水区段的渠池不存在滞后时间，模型包含两个基本参数：水面面积A_S和入流流量滞后时间τ，其特征在于：改进的方式为在ID模型的基础上加入正常取水流量q_off和事故流量q_ac，改进后的ID模型如下：

对公式(2)进行变形得到明渠事故的反演线性模型

上式中：A_S为回水区水面面积(m²)；y为下游水位相对于初值的增量(m)；x为回水区向下游方向的某处位置(m)；q_in为渠池上游入流流量相对于初值的变化量(m³/s)；q_out为渠池下游出流流量相对于初值的变化量(m³/s)；t为时间(s)；τ为入流流量滞后时间，即上游流量变化经均匀流区传递到下游回水区所用的时间(s)；

正常取水流量q_off为取水口的取水量的变化量，在事故反演计算过程中，时间相对较短，取水量近似认为无变化，因此，q_off在事故反演过程中认为为0，q_ac为事故流量相对于初始状态的增量，但初始状态一般认为明渠系统正常运行，故此该值就是事故流量本身；

对公式(3)进行离散化，得到离散后的反演线性模型如下：

DT为计算的时间间隔，k表示计算次数，即第k步，τ″_d为离散公式下的事故流量下游回水区水位影响的滞后时间序数，即事故影响滞后时间τ″含有DT的个数，由τ″/DT之后四舍五入得到；τ_d为离散公式下的上游流量下游回水区水位影响的滞后时间序数，即滞后时间含有DT的个数，由τ/DT之后四舍五入得到；i为渠池编号，q_i,in为渠池i的入流量变化量，通过实测获取数据；q_i,out为渠池i的出流量变化量，通过实测获取数据；q_i,ac为渠池i的事故流量，Z_i(k)为第k步，渠池i的下游水位，实测获取；

事故影响滞后时间τ″通过公式

计算获得，l_ac为事故具体点位距离渠池入口的距离，(m)；l_tal为事故渠池的总长度，入流流量滞后时间τ和回水区水面面积A_S是模型参数，通过系统辨识的方法求得。

2.一种利用权利要求1所述反演模型计算事故渠池的事故流量及事故位置的方法，其特征在于，包括以下步骤：

步骤1、通过圣维南方程组对需要辨识的渠池建立仿真模型；

步骤2、给定渠池渠池上游入流流量的变化量，通过步骤1中的仿真模型用直线拟合渠池下游水位上升过程，采用最小二乘法求取回水区水面面积A_S，通过渠池下游水位变化起始点获取该渠池的入流流量滞后时间τ；

步骤3、通过公式

计算得到事故具体点位置，l_ac为事故点位距离事故渠池进口的距离，Δτ表示上游开始监测到水位变化到下游开始监测到水位变化的时间差，通过对渠池实际观测得到，L为渠池实际总长，之后根据公式

计算得到事故影响滞后时间τ″；

步骤4、将步骤2和步骤3中得到的入流流量滞后时间τ、事故影响滞后时间τ″和回水区水面面积A_S代入反演线性模型中即可得到事故渠池i的事故流量q_i,ac，完成事故流量及事故位置的确定。

3.如权利要求2所述的方法，其特征在于：步骤2中，直线拟合渠池下游水位上升过程方法具体采用最小二乘法拟合，使目标函数

即最小二乘估计值，目标函数如下式所示：

上式中：y(k)为辨识渠池下游水位的实测值；u(k)为流量变化量u(k)＝q_in(k-τ)-q_out(k)，a为y(k-1)的线性系数，理论值为-1；b为u(k)前的线性系数，值为DT/A_S，通过对公式(5)的拟合求出b，通过b计算得到渠池下游回水区面积A_S。

Images