基于模擬退火-粒子群算法的輸油管道泄漏定位技術(shù)

張 靜，徐春亮，馬彥濤，石 美，劉 勇，陸 亮

1.中國(guó)石油華北油田公司第一采油廠，河北任丘 062552

2.中國(guó)石油華北油田公司第三采油廠，河北河間 062400

3.中國(guó)石油華北油田公司第五采油廠，河北辛集 052360

隨著國(guó)民經(jīng)濟(jì)的發(fā)展，管道已成為油氣運(yùn)輸?shù)闹饕绞?，我?guó)有近80%的原油和90%的天然氣依靠管道輸送，截至2019年底，我國(guó)油氣管道總里程達(dá)15×104km。其中，很多管道的服役年限已超過(guò)20年，因第三方破壞、內(nèi)部腐蝕、土壤腐蝕、雜散電流腐蝕等原因造成的管道泄漏穿孔頻繁發(fā)生[1-2]。2013年11月22日發(fā)生的中石化東黃輸油管道泄漏爆炸，造成62人死亡，136人受傷，實(shí)際上在2時(shí)12分因出站壓力下降確認(rèn)管道泄漏，在3時(shí)40分才大致確定了泄漏位置，泄漏定位的不及時(shí)造成了應(yīng)急響應(yīng)和搶險(xiǎn)救援的滯后。因此，準(zhǔn)確、及時(shí)地對(duì)管道泄漏進(jìn)行檢測(cè)和定位，對(duì)確保管道的安全運(yùn)行具有重要意義。

基于模型的泄漏檢測(cè)方法，是通過(guò)建立流體動(dòng)力學(xué)模型對(duì)管道沿線的參數(shù)進(jìn)行在線預(yù)測(cè)，將預(yù)測(cè)值與儀表測(cè)量的實(shí)際值進(jìn)行對(duì)比，從而確定泄漏發(fā)生的位置[3-5]。Abduirallman[6]利用特征線法對(duì)不同時(shí)間層次上的管道參數(shù)進(jìn)行求解，得到仿真壓力曲線，但針對(duì)不同的泄漏位置，反應(yīng)時(shí)間不一致；劉恩斌等[7]通過(guò)對(duì)網(wǎng)格和時(shí)間步長(zhǎng)的改進(jìn)，擴(kuò)大了泄漏檢測(cè)的范圍和精度，但反應(yīng)時(shí)間較慢；賀寧等[8]通過(guò)設(shè)計(jì)高增益的卡爾濾波器對(duì)管道狀態(tài)進(jìn)行實(shí)時(shí)估計(jì)，使用粒子濾波和壓縮感知的目標(biāo)跟蹤算法對(duì)預(yù)測(cè)值進(jìn)行預(yù)處理，定位誤差較小，但各參數(shù)之間的殘差選擇需要經(jīng)驗(yàn)化處理。以上研究均在室內(nèi)環(huán)道試驗(yàn)中運(yùn)行且安裝了大量的在線傳感器用于數(shù)據(jù)傳輸，而實(shí)際工況中，往往只有首末站安裝了溫度、壓力和流量?jī)x表，因此需要利用有限的測(cè)量點(diǎn)對(duì)整條管道的狀態(tài)進(jìn)行評(píng)估。在此，通過(guò)龍格-庫(kù)塔法對(duì)管道參數(shù)進(jìn)行迭代處理，得到不同斷面的溫度、壓力、流量數(shù)據(jù)，利用模擬退火優(yōu)化-粒子群算法（SA-PSO）對(duì)最優(yōu)化問(wèn)題進(jìn)行求解，并與其他尋優(yōu)方式進(jìn)行對(duì)比，建立輸油管道泄漏定位模型。

1 流體動(dòng)力學(xué)模型

忽略流體在管道徑向上的參數(shù)變化，油品在管道中的流動(dòng)可看作一元穩(wěn)態(tài)流動(dòng)，滿足質(zhì)量守恒、能量守恒、動(dòng)量守恒、狀態(tài)方程等。為了簡(jiǎn)化模型，忽略dt項(xiàng)（流動(dòng)參數(shù)與時(shí)間無(wú)關(guān)），忽略管道自身的彈性形變，同時(shí)考慮管道傾角的影響，見(jiàn)公式（1）～（3）。

式中：ρ為流體密度，kg/m3；x為管道長(zhǎng)度，m；P為管道內(nèi)壓，Pa；v為流體流速，m/s；g為重力加速度，9.8 m/s2；θ為管道傾角，rad；λ為摩阻系數(shù)；d為管道內(nèi)徑，m；s為首末兩點(diǎn)的高程差，m；h為流體的焓值，J/kg；Q為流體對(duì)外放出的熱量，J/kg。

將式（1）～（3）整理為常微分方程：

其中：

式中：T為流體溫度，K；Tm為環(huán)境溫度，K；m為質(zhì)量流量，kg/s。

由管道起點(diǎn)反復(fù)用式（6） 可計(jì)算不同截面的溫度、壓力、流量，同理從管道末點(diǎn)也可以反向計(jì)算起點(diǎn)的溫度、壓力、流量。這種求解常微分方程的方法對(duì)油氣管道的仿真準(zhǔn)確性很高。

2 泄漏定位原理

當(dāng)管道發(fā)生泄漏后，流體運(yùn)動(dòng)依然滿足動(dòng)力學(xué)方程，由于油品屬不可壓縮流體，密度基本不變，只有壓力、流量（流速）在泄漏點(diǎn)前后發(fā)生改變。由管道起點(diǎn)的T0、P0、m0計(jì)算t時(shí)刻管道末點(diǎn)的Tnt、Pnt、mnt，將Pnt與末點(diǎn)的實(shí)際壓力Pn對(duì)比，如|Pnt-Pn|＞σ，σ為報(bào)警閾值，則判定已經(jīng)泄漏；如Pnt＞Pn，則判定管道堵塞。發(fā)生泄漏后，假設(shè)泄漏點(diǎn)位置距離起點(diǎn)的軸向距離為l（且L為管道總長(zhǎng)度），根據(jù)起點(diǎn)T0、P0、m0計(jì)算泄漏點(diǎn)處的Pl，再結(jié)合末點(diǎn) Tn、mn計(jì)算末點(diǎn)壓力 P′n，如 P′n與Pn重合，則之前假定的距離l與實(shí)際的泄漏點(diǎn)位置一致。

在此，將管道泄漏定位問(wèn)題轉(zhuǎn)化為尋優(yōu)問(wèn)題，采用m組末點(diǎn)壓力計(jì)算的平均值，為了衡量個(gè)體之間的差異程度，引入偏差平方和：

3 SA-PSO算法

泄漏定位問(wèn)題轉(zhuǎn)為含不等式約束的尋優(yōu)問(wèn)題后，需要采用智能算法進(jìn)行求解。前人分別用遺傳算法（GA）、粒子群算法（PSO） 和L-M算法對(duì)管道泄漏問(wèn)題進(jìn)行了研究，但這幾種算法均容易陷入局部最優(yōu)，且收斂較慢，容易受到初始參數(shù)波動(dòng)的影響。在此，將模擬退火優(yōu)化算法（SA）與粒子群算法（PSO）相結(jié)合，利用模擬退火優(yōu)化-粒子群算法（SA-PSO）對(duì)最優(yōu)化問(wèn)題進(jìn)行求解。

式中：c1、c2為學(xué)習(xí)因子，兩者均為正數(shù)，c1為局部最優(yōu)解的步長(zhǎng)，c2為全局最優(yōu)解的步長(zhǎng)；ran（）為[0，1]的隨機(jī)數(shù)。

為了避免局部最優(yōu)，將式（8）進(jìn)行改進(jìn)，引入慣性因子ω，見(jiàn)式（10）：

全局搜索能力隨著ω的增加而增加，當(dāng)ω=1時(shí)，退化為基本的PSO算法。常規(guī)ω采用線性遞減方式，為了加快收斂速度，引入距離因子C（t）對(duì)慣性因子和學(xué)習(xí)因子進(jìn)行自適應(yīng)調(diào)整，可得式（11）：

式中：D（t）為全局最優(yōu)解的平均距離，max（D） 為全局最優(yōu)解的最大距離。

模擬退火算法是從鋼材熱處理退火過(guò)程中演化而來(lái)，通過(guò)概率接受準(zhǔn)則，控制退火中溫度的下降，反復(fù)迭代尋求近似的全局最優(yōu)解。針對(duì)SA算法可以跳出局部最優(yōu)解的限值，將SA和PSO融合，在PSO最優(yōu)粒子的選擇中加入Metropolis準(zhǔn)則，在接受最優(yōu)粒子的時(shí)候以一定概率接受較差粒子，隨著退火溫度的降低，逐漸達(dá)到全局收斂。

SA-PSO算法過(guò)程如下：

（1）對(duì)每個(gè)粒子的位置和速度進(jìn)行初始化，隨機(jī)產(chǎn)生混合矩陣作為初始粒子維度，計(jì)算初始退火溫度T0。

（2）根據(jù)模擬退火極大似然估計(jì)函數(shù)，計(jì)算當(dāng)前退火溫度下每個(gè)粒子的適應(yīng)度。

（5）采用衰減系數(shù)λ′進(jìn)行冷卻退火處理，滿足下式（Tk為第k次迭代時(shí)的溫度）：

（6）滿足收斂條件，終止計(jì)算；反之，從上述第（2）步重新迭代計(jì)算。

為了對(duì)SA-PSO算法的優(yōu)化性能進(jìn)行分析，分別計(jì)算迭代1 000步后的適應(yīng)度，見(jiàn)圖1。

圖1 不同算法優(yōu)化性能

三種算法的適應(yīng)度在初期迭代過(guò)程中迅速降低，其中標(biāo)準(zhǔn)PSO算法的適應(yīng)度在迭代至700步后降為0.032 45；改進(jìn)后的PSO算法在迭代至300步后適應(yīng)度降為0.031 23，收斂精度上比標(biāo)準(zhǔn)PSO略微提升，但收斂速度明顯提高；SA-PSO在迭代至100步后適應(yīng)度降為0.032 47，與其余兩種算法的精度相差不多，繼續(xù)迭代至510步、550步和780步分別發(fā)生了階躍變化，最終適應(yīng)度降為0.020 11。證明SA-PSO算法可以有效跳出局部最優(yōu)解的限制，在保證迭代次數(shù)的前提下，得到全局最優(yōu)解。從精度和運(yùn)算速度來(lái)看，SA-PSO的效果最好。

4 實(shí)例計(jì)算

為了考慮安全性，選擇壓力等級(jí)較低的輸油管道試驗(yàn)，為此選擇華北油田采油一廠文三十到任三聯(lián)合站的輸油管道，長(zhǎng)度1.5 km，管道材質(zhì)為20鋼，規(guī)格為D159 mm×5 mm，采用環(huán)氧粉末防腐，外加聚氨酯泡沫保溫層。管道2001年投產(chǎn)，首站溫度61℃，出站壓力0.55 MPa，末站溫度50℃，進(jìn)站壓力0.38 MPa，管輸流量25 000 kg/h，綜合含水0.5%，中間經(jīng)過(guò)兩個(gè)閥組（距離起點(diǎn)分別為350 m和680 m），首末站均有體積流量計(jì)、溫度傳感器和壓力傳感器，儀表的可重復(fù)性均為0.25%。

4.1 算法實(shí)現(xiàn)

通過(guò)開(kāi)啟閥組處的泄壓閥向泄放罐中泄漏流體模擬管道泄漏，以某一時(shí)間節(jié)點(diǎn)為起點(diǎn)，考慮到壓力噪聲因素的影響，在120 s后緩慢打開(kāi)350 m閥組處的閥門，閥門開(kāi)度10%，將上位機(jī)數(shù)據(jù)和經(jīng)Matlab-toolbox編譯的泄漏定位模型連接。測(cè)試結(jié)果見(jiàn)圖2。

圖2 泄漏前后壓力分布情況

可見(jiàn)起點(diǎn)壓力和末點(diǎn)壓力都在120 s時(shí)迅速下降，30 s后趨于重新恢復(fù)穩(wěn)態(tài)，但穩(wěn)定值小于泄漏前的壓力?？倻y(cè)試時(shí)間300 s，泄漏后重新穩(wěn)態(tài)為150 s，上位機(jī)每秒接收一次數(shù)據(jù)，共得到150組數(shù)據(jù)，為了克服數(shù)據(jù)波動(dòng)造成的影響，將每10組的平均值合并為一組數(shù)據(jù)，即150～160 s為第1組數(shù)據(jù)，160～170 s為第2組數(shù)據(jù)，以此類推，290～300 s為第15組數(shù)據(jù)。

將數(shù)據(jù)代入SA-PSO模型，分別與遺傳算法（GA）[9]、人工蟻群算法（ACO）[10]、差分進(jìn)化算法（DE）[11]進(jìn)行比較，見(jiàn)圖3。

圖3 不同算法計(jì)算時(shí)間對(duì)比

由圖3可知SA-PSO算法的計(jì)算時(shí)間最短，平均5.5s，其余3種算法的計(jì)算時(shí)間較長(zhǎng)，并且均出現(xiàn)了計(jì)算時(shí)間的波動(dòng)，可能是現(xiàn)場(chǎng)需要參照壓力表手動(dòng)控制閥門開(kāi)度，因此對(duì)泄放速度的控制有滯后性，當(dāng)開(kāi)度過(guò)大時(shí)，需要進(jìn)行多次回調(diào)，導(dǎo)致捕捉到的壓力數(shù)據(jù)有波動(dòng)，這3種算法均需迭代計(jì)算到設(shè)定的殘差才能收斂，前期的數(shù)據(jù)波動(dòng)增加了計(jì)算時(shí)的累計(jì)殘差；但SA-PSO可以及時(shí)跳出局部最優(yōu)解，對(duì)泄漏量微小差值變化的適應(yīng)性較好。此外，該試驗(yàn)管道距離較短，但實(shí)際工況中，特別是針對(duì)長(zhǎng)輸管道或復(fù)雜集輸管網(wǎng)，需要更快的算法進(jìn)行準(zhǔn)確定位，因此SA-PSO算法的適應(yīng)性更好。

4.2 不同工況對(duì)泄漏定位的影響

4.2.1 壓力和流量變化對(duì)泄漏定位的影響

測(cè)試起點(diǎn)壓力 0.5～0.7 MPa、流量 25 000～28 000 kg/h、泄漏率10%條件下的泄漏定位結(jié)果，見(jiàn)表3。

表3 泄漏點(diǎn)350 m處4種算法的泄漏定位結(jié)果

不同的壓力和流量工況下，SA-PSO算法的定位誤差在1.3～2.9 m之間，遠(yuǎn)小于其余3種算法，GA算法雖適合求解離散問(wèn)題，但存在漢明懸崖的問(wèn)題，交叉和突變均難以跨越；ACO算法沒(méi)有進(jìn)行集中控制約束，個(gè)體行為會(huì)影響整個(gè)求解過(guò)程；DE算法容易發(fā)生收斂過(guò)慢現(xiàn)象，證明SA-PSO算法的融合對(duì)管道泄漏定位具有很好的效果。

4.2.2 泄漏率變化對(duì)泄漏定位的影響

在起點(diǎn)壓力0.51 MPa、末點(diǎn)壓力0.36 MPa、流量25 000 kg/h的條件下，對(duì)泄漏率小于5%的微小泄漏檢測(cè)結(jié)果見(jiàn)表4。

表4 泄漏點(diǎn)680 m處4種算法的泄漏定位結(jié)果

隨著泄漏率的逐漸減小，定位誤差逐漸增大，SA-PSO的最大定位誤差在5.2 m，遠(yuǎn)小于其余3種算法，證明SA-PSO算法對(duì)微小泄漏的適應(yīng)性較好。同時(shí)，隨著泄漏率的上升，SA-PSO算法的平均計(jì)算時(shí)間與泄漏率呈負(fù)相關(guān)，且泄漏量越大，所需的計(jì)算時(shí)間越相近，結(jié)果符合基于流體動(dòng)力學(xué)模型的泄漏檢測(cè)原理，見(jiàn)圖4。

圖4 不同泄漏率下的計(jì)算時(shí)間

4.2.3 t檢驗(yàn)

為了對(duì)不同工況條件下測(cè)試數(shù)據(jù)與標(biāo)準(zhǔn)數(shù)據(jù)之間的差異性進(jìn)行統(tǒng)計(jì)，檢驗(yàn)流量和泄漏率對(duì)定位精度的影響，檢查上述測(cè)試數(shù)據(jù)是否存在異常值，采用t檢驗(yàn)準(zhǔn)則，判定公式如下：

式中：t為t檢驗(yàn)的計(jì)算值；avgi和avgj分別為第i組和第j組測(cè)試數(shù)據(jù)的誤差平均值；σi和σj分別為第i組和第j組測(cè)試數(shù)據(jù)的誤差標(biāo)準(zhǔn)差；x和y分別為第i組和第j組測(cè)試數(shù)據(jù)的數(shù)據(jù)總數(shù)。

采用雙尾檢驗(yàn)，如t＜t（a，x+y-2），則測(cè)試數(shù)據(jù)與標(biāo)準(zhǔn)數(shù)據(jù)無(wú)顯著性差異；否則，有顯著性差異。其中，a為置信度水平，在此取0.05，結(jié)果見(jiàn)表5。

表5 流量和泄漏率對(duì)定位精度的影響分析結(jié)果

由表5可知，在流量和泄漏率這兩種因素的影響下定位精度依然可以得到有效保證，測(cè)試數(shù)據(jù)未出現(xiàn)異常值。

5 結(jié)束語(yǔ)

將泄漏定位問(wèn)題轉(zhuǎn)為尋優(yōu)問(wèn)題，利用模擬退火優(yōu)化粒子群算法（SA-PSO）對(duì)最優(yōu)化問(wèn)題進(jìn)行求解，其中SA-PSO算法的計(jì)算時(shí)間最短，平均1.5 s，不同工況下的定位誤差在1.3～2.9 m之間，同時(shí)該算法對(duì)微小泄漏的適應(yīng)性較好。在模型設(shè)置中未考慮當(dāng)量管徑和摩阻系數(shù)的影響，實(shí)際中這兩個(gè)參數(shù)均為變化參量，今后需改進(jìn)泄漏模型和求解過(guò)程。