摘　要：基于图论理论、变流理论及贝叶斯算法建立环状热网泄漏工况水力模型、预测误差模型。在模拟管网泄漏工况时，将泄漏归结在节点处，将泄漏流量视为该节点流量的一部分，泄漏程度用泄漏面积反映。将由热网泄漏工况水力模型得到的热网泄漏数据库数据作为测量数据，将水力模型的模拟结果作为预测数据，采用预测误差模型计算预测数据的最优估计概率，指导确定泄漏位置。 

关键词：环状热网；  慢变流；  泄漏检测模型；  贝叶斯算法

Abstract：The hydraulic model and prediction error model for leakage condition of ring-shaped heatsupply network are built based on the graph theory，the variable flow theory and Bayesian algorithm．When the leakage accident conditions of pipe network is simulated，assumed that the leakage occurs at or near node，the leakage discharge is regarded as a part of the node discharge，the extent of leakage is expressed as the leakage area．The data in the heat-supply network leakage database obtained from the hydraulic model for leakage condition of ring-shaped heat-supply network are used as the measured data，the simulation results the hydraulic model are used as the predicted data，and the optimal prediction probability of the simulation results is calculated by the prediction error model to help to determine the leakage position．

Keywords：ring-shaped heat-supply network；slow transient flow；eakage detection model；Bayesian algorithm

黑龙江省地处高寒地区，冬季平均室外温度为-20～-15℃，最低室外温度达-40℃，是全国15个北方供暖区域中最寒冷的省份之一。近年来，热网泄漏的检测和定位技术得到了人们的重视^[1-3]。自20世纪70年代，许多国家开展了热网泄漏的检测和定位研究，并尝试了多种方法，取得了多项研究成果，并在工程中应用。

目前，国际上将检测和定位方法大体上分为基于硬件的方法、基于软件的方法。基于硬件的方法是指对泄漏对象直接进行检测，如直接观察法^[4-5]、检漏电缆法、放射性示踪法^[6-8]、红外线热成像检漏法^[9]、光纤检漏法^[10-11]等。基于软件的方法是指利用现代控制理论、信号处理和计算机技术等对泄漏造成的影响(如水头、流量、流速、比摩阻等参数的变化及泄漏引起的声波传输特性等)进行采集、处理和估计，从而对管道的非线性、不确定性、随机性等因素引起的误差进行补偿，进而提高泄漏检测的灵敏度和定位精度。

本文基于图论理论、慢变流理论，建立热网泄漏工况水力模型。在热网的参考节点(可选为定压点)处施加一个随时间改变的水头激励，并检测管段水头、流量变化响应。，将热网泄漏数据库数据(由热网水力模型计算得到)作为测量数据，将水力模型的模拟结果作为预测数据，采用贝叶斯算法，计算预测数据的最优估计概率，从而指导确定泄漏位置。由于施加的激励是慢变的过程，因此对管网的整体影响很小，但是对泄漏点的影响很明显，并且短时间内可采集大量数据，便于提高泄漏检测的效果，提高准确度^[12]。

在热网正常运行时，单纯的负荷变动以及调节可以抽象成非恒定流中的慢变流，通过对实际热网进行多种慢变工况的模拟，可以捕捉到丰富的信息。非恒定流动管网的基本方程为常微分方程，方程的求解相比恒定流管网更为复杂。在慢变流管网水力计算中需要引入惯性水头，其物理意义为以水头表示的由于时变惯性力做功引起的单位重量流体的能量。慢变流是指供热介质的流量、水头等参数随时间缓慢变化的流动，该流动状态下得到的参数介于恒定流与快变流之间。由于快变流对热网运行状态的影响较大，因此采用慢变流理论建立热网泄漏工况水力模型。

为了便于计算分析，我们将问题简化，设定热网中的节点流量并非单纯的供热流量，还包括与该节点相连管段上的泄漏流量。这样，我们可以将管网中的泄漏位置归结在节点处，其中泄漏流量依据泄漏面积按照孔口出流公式计算，这相当于在节点处增加一个额外的流量，而且泄漏量、泄漏位置都是未知的^[13]。节点处泄漏流量矩阵口。的计算式为：

qL=mAL

式中qL——节点处的泄漏流量矩阵，m³／s

m——孔Ll的流量系数，一般为0.60～0.62

AL——泄漏面积矩阵，m²

g——重力加速度，m／s²

hN——节点水头列向量，m

对某节点数为n+1，管段数为b的环状热网，任取一个节点作为参考节点，采用基本回路法(MKP)进行水力计算，基于慢变流的热网泄漏工况水力模型为^[13]：

Aq+qd+qL=0

A=(Ai,j)n×b

q=(qj)b×1

BDh=0

B=(Bp,j)(b-n)×b

Dh=(hj)b×1

Dh=hw+hG-Hp

hw=(hw,j)b×1

hG=C(dq/dt)

hp=(hp,j)b×1

hN=hCIN+h0

hO=(hO,i)n×1

hO=(AA^T)^-1ADh

式中A——基本关联矩阵

q——管段流量矩阵，m³／s

qd——用户的流量矩阵，m³／s

Ai,j——基本关联矩阵A的元素

qj——第j管段的流量，m³／s

B——基本回路矩阵

Dh——管段两端水头差列向量，m

Bp,j——基本回路矩阵B的元素

hj——第j管段两端水头差，m

hw——管段两端水头损失列向量，m

hG——管段惯性水头列向量，m

hp——循环泵扬程列向量，m

hw,j——第j管段水头损失，m

C——管段惯性特征矩阵，为b×b阶对角阵

t——时间，s

hp,j——第j管段循环泵扬程，m

hN——节点水头列向量，m

hC——参考节点水头，m

IN——b×1阶单位列向量，向量中每一个元素均为1

hO——节点与参考节点的相对水头列向量，m

hO,i——第节点与参考节点相对水头，

运用MATLAB程序建立热网泄漏工况水力模型，在参考节点处施加一个随时间变化的连续缓慢变化的激励(本文施加的是水头激励)，模拟计算各节点的水头与流量，建立热网节点泄漏数据库(数据作为测量数据)。设定管网中有U个水头测点、V个流量测点，并且随机选取肜个工况的水头与流量测量数据，得到U×W阶测量水头矩阵h^m与V×W阶测量流量矩阵q^m：

h^m=(h^mu,w)U×W

q^m=(q^mv,w)U×W

式中h^m——测量水头矩阵，Pa

h^mu,w——第／t个水头测点、第W个工况下的测量水头，Pa

q^m——测量流量矩阵，m³／s

q^mv,w——第v个流量测点、第w个工况下的测量流量，m³／s

当管网泄漏时(不考虑参考节点的泄漏)，可能发生泄漏的位置数量为L(L≤b)，这时泄漏位置和泄漏流量可以组成一列未知的参数集合?，参数集合?被分为两个子集，?={?L，?q}，集成到热网泄漏工况水力模型中，成为与泄漏时间密切相关的水力模型参数，用以描述管网的泄漏情况。子集?L为L维整数矢量，显示有泄漏的节点，也就是显示泄漏的位置。子集?q也是一个L维矢量，在求出泄漏流量的前提下，将泄漏流量量化在相应的三个泄漏位置。将参数集合?代入热网泄漏工况水力模型可计算得到预测水头矩阵h^w(?)、预测流量矩阵q^w(?)。根据测量结果，可计算得到预测结果的预测误差。

预测误差矩阵E(?)的计算式为：

 

式中E(?)——预测误差矩阵

h^w(?)——预测水头矩阵，m

h0——参考节点水头，m

q^w(?)——预测流量矩阵，m³／s

q⁰——参考节点流量，m³／s

h^wu,w——第u个水头测点、第w个工况下的预测水头，m

q^wv,w——第影个流量测点、第W个工况下的预测流量，m³／s

E(?)的出现是不可避免的，根据统计学的贝叶斯理论，将E(?)视为预测模型^[14]的一个随机数，并记预测模型的参数集合为s。则参数集合?与参数集合s可表示为参数集合{?，s}。在这里，通过测量值来更新参数集{?，s}，使用概率密度函数来量化它们的不确定性，进而确定预测结果的最优估计概率(即相对真实性)^[15]。

设定E(?)中的每一个元素都服从正态分布，则参数集合{?，s}的先验概率密度函数为：

 

 

式中f1——先验概率密度函数

h1^m、h2^m——第1、2个水头测点的测量水头，Pa

hU^m——第U个水头测点的测量水头，m

q1^m、q2^m——第1、2个流量测点的测量流量，m³／s

qV^m——第V个流量测点的测量流量，m³／s

h^wu,w(?)——第u个水头测点、第w个工况下的预测水头，m

q^wv,w(?)——第v个流量测点、第w个工况下的预测流量，m³／s

由贝叶斯公式可得参数集合{?，s}．的后验概率密度函数为：

 

式中f2——后验概率密度函数