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

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

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

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

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

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

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

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

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

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

其中：

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

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

2 泄漏定位原理

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

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

3 SA-PSO算法

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

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

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

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

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

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

SA-PSO算法過程如下：

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

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

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

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

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

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

三種算法的適應(yīng)度在初期迭代過程中迅速降低，其中標(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)算速度來看，SA-PSO的效果最好。

4 實例計算

為了考慮安全性，選擇壓力等級較低的輸油管道試驗，為此選擇華北油田采油一廠文三十到任三聯(lián)合站的輸油管道，長度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)過兩個閥組（距離起點分別為350 m和680 m），首末站均有體積流量計、溫度傳感器和壓力傳感器，儀表的可重復(fù)性均為0.25%。

4.1 算法實現(xiàn)

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

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

可見起點壓力和末點壓力都在120 s時迅速下降，30 s后趨于重新恢復(fù)穩(wěn)態(tài)，但穩(wěn)定值小于泄漏前的壓力。總測試時間300 s，泄漏后重新穩(wěn)態(tài)為150 s，上位機(jī)每秒接收一次數(shù)據(jù)，共得到150組數(shù)據(jù)，為了克服數(shù)據(jù)波動造成的影響，將每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)行比較，見圖3。

圖3 不同算法計算時間對比

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

4.2 不同工況對泄漏定位的影響

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

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

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

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

4.2.2 泄漏率變化對泄漏定位的影響

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

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

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

圖4 不同泄漏率下的計算時間

4.2.3 t檢驗

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

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

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

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

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

5 結(jié)束語

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