摘　要：基于圖論理論、變流理論及貝葉斯算法建立環(huán)狀熱網(wǎng)泄漏工況水力模型、預(yù)測誤差模型。在模擬管網(wǎng)泄漏工況時，將泄漏歸結(jié)在節(jié)點處，將泄漏流量視為該節(jié)點流量的一部分，泄漏程度用泄漏面積反映。將由熱網(wǎng)泄漏工況水力模型得到的熱網(wǎng)泄漏數(shù)據(jù)庫數(shù)據(jù)作為測量數(shù)據(jù)，將水力模型的模擬結(jié)果作為預(yù)測數(shù)據(jù)，采用預(yù)測誤差模型計算預(yù)測數(shù)據(jù)的最優(yōu)估計概率，指導(dǎo)確定泄漏位置。?

關(guān)鍵詞：環(huán)狀熱網(wǎng)；??慢變流；??泄漏檢測模型；??貝葉斯算法

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

黑龍江省地處高寒地區(qū)，冬季平均室外溫度為-20～-15℃，最低室外溫度達(dá)-40℃，是全國15個北方供暖區(qū)域中最寒冷的省份之一。近年來，熱網(wǎng)泄漏的檢測和定位技術(shù)得到了人們的重視^[1-3]。自20世紀(jì)70年代，許多國家開展了熱網(wǎng)泄漏的檢測和定位研究，并嘗試了多種方法，取得了多項研究成果，并在工程中應(yīng)用。

目前，國際上將檢測和定位方法大體上分為基于硬件的方法、基于軟件的方法。基于硬件的方法是指對泄漏對象直接進(jìn)行檢測，如直接觀察法^[4-5]、檢漏電纜法、放射性示蹤法^[6-8]、紅外線熱成像檢漏法^[9]、光纖檢漏法^[10-11]等?；谲浖姆椒ㄊ侵咐矛F(xiàn)代控制理論、信號處理和計算機技術(shù)等對泄漏造成的影響(如水頭、流量、流速、比摩阻等參數(shù)的變化及泄漏引起的聲波傳輸特性等)進(jìn)行采集、處理和估計，從而對管道的非線性、不確定性、隨機性等因素引起的誤差進(jìn)行補償，進(jìn)而提高泄漏檢測的靈敏度和定位精度。

本文基于圖論理論、慢變流理論，建立熱網(wǎng)泄漏工況水力模型。在熱網(wǎng)的參考節(jié)點(可選為定壓點)處施加一個隨時間改變的水頭激勵，并檢測管段水頭、流量變化響應(yīng)。，將熱網(wǎng)泄漏數(shù)據(jù)庫數(shù)據(jù)(由熱網(wǎng)水力模型計算得到)作為測量數(shù)據(jù)，將水力模型的模擬結(jié)果作為預(yù)測數(shù)據(jù)，采用貝葉斯算法，計算預(yù)測數(shù)據(jù)的最優(yōu)估計概率，從而指導(dǎo)確定泄漏位置。由于施加的激勵是慢變的過程，因此對管網(wǎng)的整體影響很小，但是對泄漏點的影響很明顯，并且短時間內(nèi)可采集大量數(shù)據(jù)，便于提高泄漏檢測的效果，提高準(zhǔn)確度^[12]。

在熱網(wǎng)正常運行時，單純的負(fù)荷變動以及調(diào)節(jié)可以抽象成非恒定流中的慢變流，通過對實際熱網(wǎng)進(jìn)行多種慢變工況的模擬，可以捕捉到豐富的信息。非恒定流動管網(wǎng)的基本方程為常微分方程，方程的求解相比恒定流管網(wǎng)更為復(fù)雜。在慢變流管網(wǎng)水力計算中需要引入慣性水頭，其物理意義為以水頭表示的由于時變慣性力做功引起的單位重量流體的能量。慢變流是指供熱介質(zhì)的流量、水頭等參數(shù)隨時間緩慢變化的流動，該流動狀態(tài)下得到的參數(shù)介于恒定流與快變流之間。由于快變流對熱網(wǎng)運行狀態(tài)的影響較大，因此采用慢變流理論建立熱網(wǎng)泄漏工況水力模型。

為了便于計算分析，我們將問題簡化，設(shè)定熱網(wǎng)中的節(jié)點流量并非單純的供熱流量，還包括與該節(jié)點相連管段上的泄漏流量。這樣，我們可以將管網(wǎng)中的泄漏位置歸結(jié)在節(jié)點處，其中泄漏流量依據(jù)泄漏面積按照孔口出流公式計算，這相當(dāng)于在節(jié)點處增加一個額外的流量，而且泄漏量、泄漏位置都是未知的^[13]。節(jié)點處泄漏流量矩陣口。的計算式為：

qL=mAL

式中qL——節(jié)點處的泄漏流量矩陣，m³／s

m——孔Ll的流量系數(shù)，一般為0.60～0.62

AL——泄漏面積矩陣，m²

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

hN——節(jié)點水頭列向量，m

對某節(jié)點數(shù)為n+1，管段數(shù)為b的環(huán)狀熱網(wǎng)，任取一個節(jié)點作為參考節(jié)點，采用基本回路法(MKP)進(jìn)行水力計算，基于慢變流的熱網(wǎng)泄漏工況水力模型為^[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——基本關(guān)聯(lián)矩陣

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

qd——用戶的流量矩陣，m³／s

Ai,j——基本關(guān)聯(lián)矩陣A的元素

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

B——基本回路矩陣

Dh——管段兩端水頭差列向量，m

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

hj——第j管段兩端水頭差，m

hw——管段兩端水頭損失列向量，m

hG——管段慣性水頭列向量，m

hp——循環(huán)泵揚程列向量，m

hw,j——第j管段水頭損失，m

C——管段慣性特征矩陣，為b×b階對角陣

t——時間，s

hp,j——第j管段循環(huán)泵揚程，m

hN——節(jié)點水頭列向量，m

hC——參考節(jié)點水頭，m

IN——b×1階單位列向量，向量中每一個元素均為1

hO——節(jié)點與參考節(jié)點的相對水頭列向量，m

hO,i——第節(jié)點與參考節(jié)點相對水頭，

運用MATLAB程序建立熱網(wǎng)泄漏工況水力模型，在參考節(jié)點處施加一個隨時間變化的連續(xù)緩慢變化的激勵(本文施加的是水頭激勵)，模擬計算各節(jié)點的水頭與流量，建立熱網(wǎng)節(jié)點泄漏數(shù)據(jù)庫(數(shù)據(jù)作為測量數(shù)據(jù))。設(shè)定管網(wǎng)中有U個水頭測點、V個流量測點，并且隨機選取肜個工況的水頭與流量測量數(shù)據(jù)，得到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

當(dāng)管網(wǎng)泄漏時(不考慮參考節(jié)點的泄漏)，可能發(fā)生泄漏的位置數(shù)量為L(L≤b)，這時泄漏位置和泄漏流量可以組成一列未知的參數(shù)集合?，參數(shù)集合?被分為兩個子集，?={?L，?q}，集成到熱網(wǎng)泄漏工況水力模型中，成為與泄漏時間密切相關(guān)的水力模型參數(shù)，用以描述管網(wǎng)的泄漏情況。子集?L為L維整數(shù)矢量，顯示有泄漏的節(jié)點，也就是顯示泄漏的位置。子集?q也是一個L維矢量，在求出泄漏流量的前提下，將泄漏流量量化在相應(yīng)的三個泄漏位置。將參數(shù)集合?代入熱網(wǎng)泄漏工況水力模型可計算得到預(yù)測水頭矩陣h^w(?)、預(yù)測流量矩陣q^w(?)。根據(jù)測量結(jié)果，可計算得到預(yù)測結(jié)果的預(yù)測誤差。

預(yù)測誤差矩陣E(?)的計算式為：

?

式中E(?)——預(yù)測誤差矩陣

h^w(?)——預(yù)測水頭矩陣，m

h0——參考節(jié)點水頭，m

q^w(?)——預(yù)測流量矩陣，m³／s

q⁰——參考節(jié)點流量，m³／s

h^wu,w——第u個水頭測點、第w個工況下的預(yù)測水頭，m

q^wv,w——第影個流量測點、第W個工況下的預(yù)測流量，m³／s

E(?)的出現(xiàn)是不可避免的，根據(jù)統(tǒng)計學(xué)的貝葉斯理論，將E(?)視為預(yù)測模型^[14]的一個隨機數(shù)，并記預(yù)測模型的參數(shù)集合為s。則參數(shù)集合?與參數(shù)集合s可表示為參數(shù)集合{?，s}。在這里，通過測量值來更新參數(shù)集{?，s}，使用概率密度函數(shù)來量化它們的不確定性，進(jìn)而確定預(yù)測結(jié)果的最優(yōu)估計概率(即相對真實性)^[15]。

設(shè)定E(?)中的每一個元素都服從正態(tài)分布，則參數(shù)集合{?，s}的先驗概率密度函數(shù)為：

?

?

式中f1——先驗概率密度函數(shù)

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個工況下的預(yù)測水頭，m

q^wv,w(?)——第v個流量測點、第w個工況下的預(yù)測流量，m³／s

由貝葉斯公式可得參數(shù)集合{?，s}．的后驗概率密度函數(shù)為：

?

式中f2——后驗概率密度函數(shù)