基于CUSUM 方法的集中供熱管網泄漏故障診斷

周守軍張國正劉臣林旭東曹長紅

(1.山東建筑大學熱能工程學院，山東 濟南 250101；2.華能山東發電有限公司白楊河發電廠，山東 淄博 255200；3.淄博市博山區熱力公司，山東 淄博 255200)

0 引言

城鎮集中供熱管道是城市地下管線的重要組成部分，相比于自來水管道，其應用時間短，但發展卻較為迅速。 據統計，截至2019 年，我國運營的熱水管網總長達48.8×104km，其中一次網長度為12.4×104km、庭院管網長為36.4×104km[1]。 由于管網的運行狀況愈加復雜，保障這些管網的運行安全就顯得尤為重要。 其中，泄漏是影響管網運行安全的主要問題，其主要原因包括腐蝕、設備老化、機械沖擊等[2-3]，大量介質的泄漏會導致經濟損失和環境破壞。 由于管道性能會隨著時間的推移而降低，不可避免地產生局部泄漏，因此及時有效的泄漏故障診斷在實踐中非常重要。

最初，我國大部分供熱企業處于“先漏后檢”的階段，一般是由工作人員“聽聲判漏”，或結合管網壓力變化試挖泄漏位置，這種方式的效率低、耗時長、成本高且定位不夠精準[4]。 隨著研究深入，管道泄漏診斷技術根據其作用的位置及實施手段主要可分為管道內和管道外診斷。 管道外診斷方法是將探測儀器置于在地表或管道外監測周圍環境因素，獲取數據并進行分析，從而及時地診斷泄漏的發生。楊漢瑞等[5]提出一種基于分布式光纖傳感技術的膨脹節膨脹量的泄漏診斷方法，建立了有效反映膨脹量與溫度變化關系的檢測模型，從而診斷了泄漏位置，其診斷結果與膨脹量的測量誤差相關。 吳晉湘等[6]利用紅外熱像儀和土壤濕度速測儀等記錄異常管段沿線土壤的溫度、濕度信息，分析管網泄漏和保溫層破損時傳熱傳質對管道上方土壤溫濕度的影響規律，判定泄漏點位置。 PERPAR 等[7]利用土壤的溫度梯度診斷供熱管道的泄漏狀況，其精度較高。 FRIMAN 等[8]通過無人機搭載遙感熱相機，將檢測到具有異常高溫度的區域認定為泄漏區域，此方法可以自動分割影響檢測結果的建筑物，從而減小誤報率。 上述管道外診斷方法有一定的效果，但受外部環境因素影響較大，且大部分存在排查周期長、成本過高、實施過程復雜以及確定泄漏位置不夠及時等缺點。

管道內診斷法發展源于20 世紀末，主要利用監測設備置于管內獲取管網運行數據進行直接診斷，或結合其他學科開展數據分析及泄漏診斷。FILIPPINI 等[9]基于壓力波和質量流量數據分析，直接診斷泄漏時間和泄漏位置，但其準確度較差。PéREZ-PéREZ 等[10]提出人工神經網絡結合在線監測數據的管網泄漏診斷方法，利用模擬數據來訓練人工神經網絡，并結合數學模型估計管網阻力系數，得到較好的診斷結果，但所需數據量較多，對模型的精度要求也較高，且管網運行中的壓力波動對泄漏定位的精度有較大影響。 石光輝等[11]利用高頻壓力傳感器，通過捕捉管網泄漏瞬間形成的負壓水擊波，并根據小波分析方法得到負壓波傳播到傳感器的時間，再以此定位泄漏位置，但供熱管網距離長、分支多，負壓水擊波會產生較大的衰減，難以對其準確定位。

供熱管網泄漏初期，補水流量及壓力變化相對較小，泄漏情況難以確定。 累積和(Cumulative Sum，CUSUM)方法基于似然比推出，通過分段累積監測數據序列信息，將過程的較小偏移疊加，以達到“質變”的效果，從而實現對較小偏移的檢測[12]。CUSUM 方法已成功應用于金融分析[13]、醫學診斷[14]、地震探測[15]及給水管網爆管檢測[16]等領域。文章提出利用CUSUM 方法和集中供熱管網泄漏仿真模型，采用仿真模擬與實驗研究相結合的方法，通過監測管網補水流量和壓力開展泄漏故障診斷，旨在泄漏發生初期就能對供熱管網進行泄漏故障報警和準確定位。

1 泄漏故障診斷方法及管網仿真

1.1 集中供熱管網仿真模型

文章提出的泄漏故障診斷方法的基礎是集中供熱管網仿真模型[17-18]。 設任意一個供熱管網，其支路數為m、節點數為n＋1，由網絡圖論理論、基爾霍夫定律、伯努利方程和管路特性方程，得到其水力工況仿真模型方程組，由式(1)表示為

式中A為管網關聯矩陣(為n×m階矩陣)；G為管段流量向量；B為管網的基本回路矩陣，為(m-n)×m階矩陣；ΔH為管段阻力壓降；S為管段阻力特性系數矩陣(m階對角矩陣)；｜G｜為管段流量G的絕對值m階對角矩陣；Z為管段支路中兩節點的位能差向量；DH為管段的水泵向量，當管段不含水泵時，該管段DH＝0，當有水泵時，DH為水泵揚程；Q為管網各節點的凈出流量向量，入流為正，出流為負。 當Q為零向量時，表示供熱管網正常運行，管網各節點的出流量與入流量代數和為0，即管網為一個閉合環路；當Q中某元素為負值，表示該管網發生泄漏。

而S可由式(2)表示為

式中K為管壁的當量絕對粗糙度，對于供熱管道，一般為0.000 5 m；d為管道內徑，m；l、ld分別為管網計算管段的長度和局部阻力當量長度，m；ρ為管內流體介質的平均密度，kg/m3。

式(1) 對應的求解模型由式(3)表示為

式中M為馬克斯威(Max Well)矩陣，是以B為基礎的(m-n)×(m-n)階的對稱正定矩陣，M矩陣對應于一定的樹，不同的樹相對應的M矩陣也不同；ΔHk為基本回路管段壓降代數和，當Gk為方程組的解時，其值為0。

上述管網仿真模型方程組包括正常工況模型、節點泄漏模型和管段泄漏模型。 其中節點泄漏是指泄漏位置在供回水干管到用戶(換熱站)的分支節點上，原管網的拓撲結構與正常工況一致，在節點泄漏仿真時，只需對該節點設定泄漏流量值即可。 管段泄漏是指泄漏位置在供回水干管管段上，在管段泄漏仿真時，泄漏位置處設置流量接近零的“虛擬用戶”，將管段泄漏轉化成“虛擬用戶”的分支節點泄漏，管網拓撲結構發生改變，在其分支節點處給定泄漏流量值，建立虛擬節點泄漏仿真模型，即為管段泄漏模型。

1.2 基于CUSUM 方法的泄漏故障診斷方法

1.2.1 CUSUM 方法

假設x1，x2，…，xζ，…，xn為按順序記錄的一組過程監測樣本數據，對這組數據做以下假設，由式(4)[12]表示為

式中H0為原假設，表示過程無變點；H1為備擇假設，表示過程有變點；μ0、μ1分別為變點出現前后樣本均值，變點即樣本數據突變的點，對應供熱管網中泄漏發生的點；σ為樣本標準差，假設為定值；n為樣本數據組別上限；ζ為變點發生處數據組別，ζ＜n。

H1對H0的似然比可由式(5)表示為

式中Ln，ζ為備擇假設H1對原假設H0的似然比統計量；L(n，μ0)、L(n，μ1)分別為H0與H1的似然函數；f0(xi)、f1(xi)分別為H0與H1的密度函數。

上、下偏移累積和由式(6)表示為

1.2.2 泄漏故障診斷方法

文章提出的泄漏故障診斷方法是利用CUSUM方法監測管網補水流量變化，當管網補水流量累積和變化超過預定閾值時，即可發出報警信號，并確定泄漏發生時間及泄漏量，然后將泄漏流量帶入管網泄漏仿真模型，對每個節點進行泄漏模擬計算，并根據給定的判定準則，比較分析采樣壓力表的實驗運行值與仿真模擬值，確定泄漏位置。

(1) 診斷是否發生泄漏

①獲取一定時間間隔(5 min)的管網補水流量監測數據序列qi，計算其均值μ0與方差σ0，并將數據序列標準化為yi，由式(7)～(9)表示為

②選取CUSUM 累積和報警閾值h和漂移參數k，計算上偏移累積和，由式(10)表示為

(2) 診斷何時發生泄漏以及泄漏程度

②根據泄漏開始時刻補水流量值q(is)和泄漏最大時刻補水流量值q(if)，計算最大泄漏流量Δqm＝q(if)-q(is)。

(3) 診斷泄漏發生的準確位置

①計算泄漏開始時刻is(正常工況)和泄漏最大時刻if(泄漏工況)各個采集點實際壓力變化值ΔHMj和仿真壓力變化值ΔHsimMj，由式(11)和(12)表示為

式中Mj為壓力采集點序號，j＝1，2，…，r，其中r為壓力采集點總個數；HMj(is)、HMj(if)分別為is、if時刻該采集點實際壓力值，由實際數據得到；分別為is、if時刻該采集點仿真壓力值，將最大泄漏流量代入泄漏仿真模型可以得到。

假設管網拓撲結構中所有分支節點都為泄漏點，則目標函數OFx由式(13)[19]表示為

式中x為節點編號；N為管網總節點數；OFx即為假設第x分支節點為泄漏點所得目標函數值；Mp(p＝j＋1，…，r)為序號大于Mj的采樣壓力表節點。

②比較所有的目標函數值OFx，找到其最小值OFmin，則對應節點為一級泄漏點

③利用二級管段泄漏模型，重復(3)中①和②，只計算一級泄漏點及其兩側的“虛擬節點”，找到最小目標函數值對應的節點，即為二級泄漏點，再根據與的具體距離L，準確定位泄漏位置。

具體泄漏故障診斷流程圖如圖1 所示。

圖1 泄漏故障診斷流程圖

2 管網泄漏實驗與模型驗證

2.1 實驗系統

實驗系統由集中供熱實驗管網和在線數據監測系統兩部分組成(如圖2 所示)。 實驗管網包含18 個用戶和3 個熱源，熱源由循環水泵和電加熱水箱組成。 管道布置為上供下回，供回水管段對稱布置，且二者管徑相同。 定壓方式為變頻水泵補水定壓，定壓點為循環水泵入口處，在泄漏實驗運行的過程中，可實現自動補水。

圖2 集中供熱實驗管網示意圖

實驗管網各用戶管段上依次安裝有進口壓力表、溫度表、電磁流量計、調節閥、截止閥和出口壓力表。 其中，截止閥用來模擬各用戶阻力，調節閥用來進行實驗管網初調節，使用戶達到設計流量。 循環水泵的前后位置依次安裝進口溫度表、進口壓力表、出口壓力表、電磁流量計和出口溫度表，補水泵處亦有進口壓力表、電磁流量計和出口壓力表。 在線數據監測系統即對上述各儀表進行實時監測并將實驗數據采集、存儲到MySQL 數據庫中。

2.2 實驗方案設計

文章進行單熱源枝狀供熱管網泄漏工況模擬實驗。 所用單熱源枝狀管網拓撲結構如圖3 所示，其中n、s 分別表示節點和管段。 該管網包括1 號熱源和18 個用戶，共分為36 個管段、38 個節點。 實驗管網設計總流量為18 m3/h，而每個用戶設計流量為1 m3/h。

圖3 單熱源枝狀實驗管網拓撲結構圖

實驗管網共13 個泄漏點，根據距離熱源由近到遠、先節點后管段、先供水后回水的原則編號。 其中節點泄漏點為3 個，在用戶的進口或出口節點，編號為1～3；管段泄漏點為10 個，在供、回水管路上各5 個，編號為4～13。 具體泄漏點位置及編號如圖3所示(括號內編號表示泄漏點在回水管路)。

城市供熱管網正常工況下補水流量不宜大于總循環水量的1%[20]，而補水流量＞5% 時，表明管網發生了嚴重泄漏，此時泄漏位置也相對容易被確定。為了驗證文章泄漏診斷方法對不同泄漏工況的有效性，并使實驗工況研究更豐富，此次實驗根據不同泄漏率(泄漏流量占總流量的百分比)分為1.1%(工況1)、2.5%(工況2)、4.0%(工況3)和5.5%(工況4)4 種泄漏工況，13 個泄漏點共52 組實驗。 每組實驗分為正常工況和泄漏工況兩部分，各運行1 870 s，數據采集頻率為2.2 s/次。 采用自適應遞推最小二乘濾波算法[21]對管網運行數據(主要包括壓力、補水流量)進行預處理，消除噪聲。

2.3 泄漏仿真模型實驗驗證

針對供熱管網的正常工況，龔璞等[18]已經開展了模擬研究，而文章主要對一級節點泄漏模型和二級管段泄漏模型進行驗證。

節點泄漏工況選擇1 號泄漏點工況2，管段泄漏工況選擇11 號泄漏點工況1。 將兩組實驗中泄漏工況數據各取300 組，計算其平均值，得到各節點泄漏工況下的運行壓力以及循環水泵揚程、定壓點壓力、泄漏流量等數據；再將循環水泵揚程、定壓點壓力、泄漏流量分別代入已編譯的節點和管段泄漏仿真模型MATLAB 程序中，得到管網各節點仿真與實驗壓力數據。 如圖4(a)所示，節點泄漏工況下各節點壓力實驗值與仿真值的相對誤差最大、最小值分別為4.71%和0.86%；如圖4 (b)所示，管段泄漏工況下各節點壓力實驗值與仿真值的相對誤差最大、最小值分別為4.64%、0.06%。 兩級泄漏仿真模型最大相對誤差均不超過5%，且從圖4 中可以看出，仿真模型的供回水壓力曲線與實驗曲線的變化基本相符，說明仿真模型的準確性較高。

圖4 不同泄漏工況實驗與仿真數據對比

3 泄漏故障診斷方法實驗研究

3.1 CUSUM 方法參數選取

在泄漏故障診斷方法中，閾值h和漂移參數k的選取至關重要。 如果參數選取過小，會使得泄漏報警靈敏度過高，可能會將管網正常狀態顯示為泄漏狀態，增加誤報率；而參數選取過大，又將降低靈敏性，可能會將管網泄漏狀態誤認為是正常狀態，增加漏報率。 主要采用3 種方法選取h和k。

方法1 為基于公式的推導[22]，由式(14)表示為

式中α為誤報率；β為漏報率；δ＝(μ1-μ0)/σ，為變點的標準化移位量。 此方法期望值為α≤10%、β≤5%。

方法2 為基于平均運行鏈長(Average Run Length， ARL)的選取[23]。 ARL 指從檢測開始到檢測出發生異常所需要的平均樣本個數，取ARL 為150，得到h＝3.5、k＝0.5。

方法3 為基于經驗值的選取[16]，h＝5σ，k＝1.015。

以3 號泄漏點的4 組工況為例，將上述3 組h、k參數值代入上述泄漏診斷方法，計算泄漏報警時間及最大泄漏流量。 泄漏報警時間結果見表1，最大泄漏流量如圖5 所示。

表1 不同方法所得泄漏報警時間單位:s

從表1 可以看出，3 組不同h、k值得到的泄漏報警時間一致，都是在泄漏發生2.2 s 后即發出報警信號，泄漏報警效果良好。 由圖5 可以看出:工況1 和3中，方法3 得到的最大泄漏流量與實驗值基本一致，另外兩種方法計算值則比實驗值要低；工況2 中，方法1 的計算值比實驗值小，另外兩種方法計算結果與實驗值相等；工況4 中，方法3 計算值與實驗值相等，另外兩種方法計算值都比實驗值稍低。 綜上所述，方法3 的整體計算結果更接近實驗值，故選擇參數h＝5σ、k＝1.015 進行泄漏故障診斷方法的應用。

3.2 實驗結果分析

3.2.1 泄漏報警計算

以3 號泄漏點4 組工況為例，按照上述方法對實驗管網補水流量數據進行泄漏故障診斷，實驗管網補水流量運行數據及CUSUM 累積和曲線如圖6所示。 可以看出，當管網發生泄漏時，管網補水流量明顯增大。 經過累積和計算，4 組工況在實驗運行第1 872.2 s 時均有，并在此之后累積和不斷增大，因此可以判定此時發生泄漏，與實際情況相符。

圖6 補水流量運行監測數據及累積和曲線圖

利用此泄漏診斷方法對13 個泄漏點共52 組泄漏實驗計算了泄漏報警時間，發現所有工況實驗均在實驗運行第1 870.0 s 時發生泄漏，而計算結果顯示泄漏報警信號發出的時間均為實驗運行第1 872.2 s，沒有誤報、漏報或者延遲報警。 由于現有實驗臺的壓力表采集頻率最高只能達到2.2 s/次，若提高采集頻率，報警時間還能縮短。

泄漏報警之后，計算管網泄漏流量。 從圖6 可以看出，工況1 ～4，計算得到泄漏開始時刻is與實際時刻一致，泄漏最大時刻if也是泄漏工況運行穩定后的時間，所以兩個時刻對應補水流量差值即為最大泄漏流量。

所有節點及工況最大泄漏流量診斷結果如圖7所示，其中實驗值是實驗過程中管網最大泄漏流量值，計算值即為泄漏診斷得到的管網最大泄漏流量值。 4 種工況下計算值與實驗值的誤差均較小，其中相對誤差最大值為1.6%，發生在2 號泄漏點工況1 和10 號泄漏點工況1，其余工況誤差＜1.0% 。 泄漏流量計算值與實驗值誤差越小，說明最大泄漏時間if的計算也就越準確，這也為泄漏點位置準確診斷提供了良好的基礎。

圖7 管網最大泄漏流量診斷結果圖

3.2.2 泄漏點位置診斷

根據泄漏位置診斷準則可以看出，選取合適的壓力采集點，是準確定位的關鍵。 由于供熱管網泄漏位置的出現沒有規律可循，所以壓力采集點的分布應當能反映泄漏發生后管網整體的壓力變化，各個壓力采集點距離過近或過遠都有一定的局限性。 利用距管網熱源近、中、遠的采集點選取原則，詳細對比分析不同的壓力采集點分布方案，以期得到一個合理的分布。 具體分布方案見表2。

表2 不同壓力采集點分布方案設計

將泄漏診斷所得位置與實際泄漏位置的誤差距離分為l≤2 m、l≤4 m 以及l≤6 m，以方案1 為例，得到診斷結果見表3、4 及6(√表示診斷結果在誤差范圍內，×表示診斷結果在誤差范圍外)。

表3 方案1 泄漏位置診斷結果( l ≤2 m)

表4 方案1 泄漏位置診斷結果( l ≤4 m)

表5 方案1 泄漏位置診斷結果( l ≤6 m)

由上述泄漏工況實驗條件可知，每組工況的泄漏位置診斷準確與否可視為隨機事件，因此定義泄漏位置診斷準確率PA為診斷準確的工況數占總工況數的百分比[24]。 令PAi表示第i個方案的泄漏位置診斷準確率，各方案泄漏位置診斷準確率PA如圖8 所示。 誤差距離l≤2 m 時，PA1＝PA2＝PA5＞PA4＞PA3；誤差距離l≤4 m 時，PA1＞PA2＞PA5＞PA4＞PA3；誤差距離l≤6 m 時，PA1＞PA2＞PA3＝PA5＞PA4。

圖8 各方案泄漏位置診斷準確率圖

分析不同壓力采集點方案診斷結果可以得到:

(1) 方案2 相比于方案1，減少了2 個壓力采集點，但PA降低不多，而方案3 比方案1 增加了4 個壓力采集點，反而PA降低較多，說明PA跟壓力采集點個數有直接的關系，壓力采集點個數過多或過少都會影響位置診斷效果；

(2) 方案4 與1 壓力采集點個數相同，但位置分布不同，導致PA整體下降較多，說明當壓力采集點個數固定時，分布位置太過靠近熱源，也會使PA降低；

(3)方案5 與1 相比，只選取了供水管段上的3個壓力采集點，采集點個數減少了3 個，但位置分布與方案1 一致，PA卻發生下降，說明只在供水管段或回水管段設置壓力采集點，都不能準確反映出泄漏工況發生時管網的整體變化，會影響泄漏位置診斷效果。

通過以上壓力采集點分布方案的泄漏位置診斷準確率PA對比發現，壓力采集點個數并不是越多越好，當管網規模較小，若壓力采集點個數太多，則兩點間距離近，因泄漏產生的壓力變化并不明顯，反而會對泄漏位置診斷結果產生較大影響，應該合理分布，使各壓力采集點正確反映出管網泄漏時整體壓力變化，這樣才能使泄漏位置診斷結果更準確。

最終確定方案1 為最優方案，其誤差距離l≤2 m時，PA＝57.7%；l≤4 m 時，PA＝67.3%；l≤6 m時，PA＝73.1%。 隨著誤差距離的增大，PA也增大，誤差范圍＞6 m 時，對診斷結果基本無影響。

泄漏位置判定的影響因素包括實驗運行壓力和仿真模擬壓力。 其中，由于循環水泵工作點波動、儀表誤差等，難以完全消除實驗運行壓力對泄漏位置診斷的影響；但仿真模型可以通過提高其精度、減小仿真模擬壓力對泄漏位置診斷的影響，進而提高位置診斷準確度。

不論連續還是間斷補水，只要能獲取到補水流量監測數據序列，即可利用文章方法診斷泄漏故障。雖未驗證環狀供熱管網，但環狀管網與枝狀管網在泄漏工況下的壓力和補水流量變化規律一致，故文章方法對環狀管網同樣適用。

4 結論

基于變點理論中的CUSUM 方法，結合集中供熱管網泄漏仿真模型，提出了一種包括泄漏報警、泄漏流量計算以及泄漏位置診斷的集中供熱管網泄漏故障診斷方法，并利用MATLAB 軟件進行了算法開發。 結合集中供熱實驗管網設計了實驗方案，得到單熱源枝狀管網不同泄漏位置及不同泄漏流量下的壓力與流量數據，并驗證了該方法，得到結論如下:

(1) 在利用合適的CUSUM 參數值進行泄漏情況報警時，泄漏工況發生后2.2 s 即可報警，且沒有誤報和漏報，說明文章方法可以及時對管網泄露報警；

(2) 文章方法診斷所得泄漏開始時間準確，泄漏最大時間也是泄漏工況運行穩定的時間，對應計算得到最大泄漏流量與實驗值相對誤差最大僅為1.6%，說明方法對泄露時間和最大泄漏流量的診斷較為準確；

(3) 壓力采集點的分布距離應該適中，采集點個數也并非越多越好，采用其中的最優方案診斷泄漏位置，其準確率可達73.1%，說明文章方法對泄漏位置的診斷準確率較高。