基于二元逆高斯過程的腐蝕輸油管道剩余壽命預(yù)測

王藝斐 蘇 春 謝明江

(東南大學(xué)機(jī)械工程學(xué)院，南京 211189)

管道是一種高效的石油和天然氣運輸方式.然而油氣管道長期在高溫、高壓下工作，并且油氣中的H2S、CO2、Cl-等腐蝕性介質(zhì)導(dǎo)致管道內(nèi)壁腐蝕.此外，管道外壁長期經(jīng)受不同性質(zhì)土壤、雜散電流的侵蝕.在多重因素的共同作用下，管壁會出現(xiàn)腐蝕、開裂、穿孔等現(xiàn)象[1].在突發(fā)的外力作用下或當(dāng)油氣的運行壓力大于管壁承載能力時，管線會發(fā)生破裂并導(dǎo)致油氣泄漏，造成環(huán)境污染和生態(tài)破壞[2].因此，開展管道剩余壽命預(yù)測、有效減少管道失效事件發(fā)生具有重要意義.

近年來，基于可靠性理論的油氣管道剩余壽命預(yù)測研究受到關(guān)注.Mosallam等[3]利用概率模型表征腐蝕尺寸和管道特性的不確定性.Gong等[4]考慮管段泄漏、爆裂失效模式和不同腐蝕缺陷下極限狀態(tài)函數(shù)的相關(guān)性，采用一階可靠度方法評價含腐蝕缺陷管道系統(tǒng)可靠性.Bouledroua等[5]采用二階矩方法評估受腐蝕管道因內(nèi)壓過大而失效的可靠性問題.但上述方法存在不足，主要表現(xiàn)在管道壽命預(yù)測精度不高、不適合于某些極端條件等[6].

管道腐蝕常表現(xiàn)為壁厚減薄、強(qiáng)度下降，是典型的性能退化過程.一些學(xué)者采用隨機(jī)過程描述管壁性能退化的不確定性.張新生等[7]基于逆高斯(IG)過程建立油氣管道腐蝕退化的狀態(tài)空間模型，開展管道壽命預(yù)測和維修決策優(yōu)化.Bazán等[8]提出非線性隨機(jī)過程的腐蝕增長模型，通過對比線性隨機(jī)變量腐蝕模型，驗證模型的擬合效果.Wang等[9]建立埋地管道結(jié)構(gòu)隨機(jī)腐蝕增長模型，采用幾何布朗橋過程模擬動態(tài)外部點蝕過程，用于地下管線結(jié)構(gòu)的完整性管理.Ossai等[10]采用非齊次線性增長純生馬爾可夫模型預(yù)測油氣管道內(nèi)腐蝕的深度分布.Amaya-Gómez等[11]基于在線檢測數(shù)據(jù)，集成應(yīng)用Lévy過程、伽馬過程和復(fù)合泊松過程構(gòu)建管道腐蝕退化模型.Zhang等[12]根據(jù)地下能源管道腐蝕缺陷的檢測數(shù)據(jù)，基于IG過程建立腐蝕深度增長模型，采用馬爾可夫鏈蒙特卡洛仿真得到預(yù)測結(jié)果.

長輸石油管道退化通常包含腐蝕穿孔、局部爆裂等多種失效模式，并且有不同的性能指標(biāo)[13].為更準(zhǔn)確地描述管道健康狀態(tài)，有必要采用2個或多個性能指標(biāo)表征管道健康狀態(tài)，并考慮性能指標(biāo)的相關(guān)性.Copula函數(shù)適合于描述多性能指標(biāo)間的相關(guān)性.Zhou等[14]利用Gaussian Copula、隨機(jī)過程理論描述管道缺陷生長之間的依賴關(guān)系.金曉航等[15]以Copula函數(shù)定義軸承2個性能指標(biāo)的相關(guān)性，得到軸承剩余壽命的聯(lián)合概率密度函數(shù).

本文以長輸石油管道為研究對象，考慮管道剩余壁厚、剩余強(qiáng)度2個退化性能指標(biāo)，提出一種基于二元IG過程的腐蝕輸油管道剩余壽命預(yù)測方法.首先采用Copula函數(shù)描述性能指標(biāo)之間的相關(guān)性，建立管道剩余壽命的聯(lián)合密度函數(shù)；采用期望值最大化(EM)算法估計失效概率密度函數(shù)中的未知參數(shù)；采用貝葉斯方法更新參數(shù)分布，完成石油管道剩余壽命的實時預(yù)測.最后通過數(shù)值案例驗證所提方法的可行性與有效性.

1 基于逆高斯過程的管道退化建模

1.1 逆高斯過程

逆高斯過程是一種具有單調(diào)性的隨機(jī)過程[16].假設(shè)產(chǎn)品退化量Y(t)服從逆高斯過程，它具有如下性質(zhì).

1) 當(dāng)t=0時，Y(0)=0.

2) 對任意t4>t3>t2>t1，Y(t4)-Y(t3)與Y(t2)-Y(t1)相互獨立.

其中，均值函數(shù)Λ(t)為一個單調(diào)增函數(shù)，ΔΛ=Λ(t+Δt)-Λ(t)，Λ(0)=0.逆高斯分布IG(βΛ(t),ηΛ2(t))的概率密度函數(shù)為

(1)

式中，β為退化率；η為尺寸參數(shù).

1.2 基于二元逆高斯過程的性能退化模型

采用多退化量預(yù)測設(shè)備剩余壽命時，需要考慮退化量之間的相關(guān)性.Copula函數(shù)是一種連接一元邊緣分布函數(shù)及其對應(yīng)多元聯(lián)合分布函數(shù)的特殊函數(shù)，通過構(gòu)建每個邊緣退化過程的聯(lián)合分布函數(shù)來定義多性能退化量之間的相關(guān)性.

1.2.1Copula函數(shù)定義

當(dāng)函數(shù)C(u,v)滿足如下性質(zhì)時，稱為二元Copula函數(shù)[17]:

1)C(u,v)的定義域為[0, 1]×[0, 1].

2)C(u,v)具有零基面，且二維遞增.

3) 對于任意的u,v∈[0, 1]，滿足C(u, 1)=u,C(1,v)=v.

當(dāng)已知各個變量的邊緣分布函數(shù)時，可根據(jù)Sklar定理得到它們的聯(lián)合分布函數(shù).

1.2.2Sklar定理

令F(x1,x2, …,xN)為具有邊緣分布F1(x1),F2(x2),…,FN(xN)的N元聯(lián)合分布函數(shù)，則存在一個Copula函數(shù)C(u1,u2, …,uN)滿足[18]

F(x1,x2,…,xN)=

C[F1(x1),F2(x2), …,FN(xN)]

(2)

若F1(x1),F2(x2),…,FN(xN)為連續(xù)函數(shù)，則C(u1,u2, …,uN)唯一確定；反之，若F1(x1),F2(x2), …,FN(xN)為一元分布函數(shù)，C(u1,u2, …,uN)是一個Copula函數(shù)，則由式(2)確定的F(x1,x2, …,xN)是具有邊緣分布F1(x1),F2(x2),…,FN(xN)的N元聯(lián)合分布函數(shù).

由Sklar定理可知，采用Copula函數(shù)描述2個退化性能指標(biāo)之間的相關(guān)特性可表示為[19]

F(x1,x2)=C[F1(x1),F2(x2)]

(3)

式中，F(xiàn)1(x1)、F2(x2)分別為2個性能指標(biāo)對應(yīng)的分布函數(shù).

2 管道剩余壽命預(yù)測模型

采用首達(dá)時間(first hitting time)表征退化量與壽命L之間的關(guān)系，并據(jù)此開展剩余壽命預(yù)測.不失一般性，假設(shè)管道退化量隨時間逐漸增加，其失效閾值為ω.當(dāng)退化量超過ω時即認(rèn)為發(fā)生失效，將管道剩余壽命預(yù)測問題轉(zhuǎn)化為預(yù)測退化量何時會到達(dá)失效閾值ω.

給定失效閾值ω，管道的壽命L可定義為[20]

L=inf{lk:Y(tk+lk)≥ω|Y(0)<ω}

(4)

式中，lk為tk時刻的剩余壽命.

設(shè)剩余壽命的PDF函數(shù)為fLk(lk|Yk)，則有

(5)

式中，Yk= {y0,y1, … ,yk}為tk時刻的退化數(shù)據(jù)集，yk為tk時刻的退化數(shù)據(jù)；μβ為均值；σβ為標(biāo)準(zhǔn)差；φ(x)為標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)；Φ(x)為標(biāo)準(zhǔn)正態(tài)分布累計分布函數(shù).

2.1 管道剩余壁厚的計算判據(jù)

當(dāng)管道剩余壁厚低于某臨界值時，油氣管道可能會發(fā)生泄漏失效.根據(jù)管道完整性管理規(guī)定，通常將腐蝕深度超過管壁厚度的80%作為管道需要維修的判據(jù)[21].通常，管道公司會不定期開展管道檢測以獲得腐蝕深度等基礎(chǔ)數(shù)據(jù).本文以管道剩余壁厚作為退化性能指標(biāo).考慮到不同性能指標(biāo)的量綱存在差異，為便于分析進(jìn)行歸一化處理，剩余壁厚退化量定義如下：

(6)

式中，δ為剩余壁厚退化量；d0為初始壁厚；d(t)為t時刻時的管道壁厚.

2.2 管道剩余強(qiáng)度

文獻(xiàn)[11-13]多采用爆破壓力表征管道的承載能力，當(dāng)所計算的爆破壓力低于管道的運行壓力時管道將會失效.相關(guān)的國際標(biāo)準(zhǔn)包括B31G、DNV F101和PCORRC.根據(jù)上述標(biāo)準(zhǔn)可以計算管道爆破壓力[22].近年來，有限元方法在管道爆破壓力計算中受到重視.

但是，在某種條件下爆破壓力并不適用于作為管道剩余強(qiáng)度退化指標(biāo).在運行過程中，隨著腐蝕缺陷的擴(kuò)展，管道承壓能力不斷降低、剩余強(qiáng)度下降.本文綜合考慮等效應(yīng)力和屈服強(qiáng)度來描述剩余強(qiáng)度退化，采用有限元方法完成油氣管道建模，通過仿真得到不同腐蝕缺陷下管壁的等效應(yīng)力.當(dāng)最大等效應(yīng)力超過屈服強(qiáng)度時，管道將發(fā)生斷裂失效.管道剩余強(qiáng)度退化量的定義如下：

(7)

式中，ε為剩余強(qiáng)度退化量；G為屈服強(qiáng)度；τmax為最大等效應(yīng)力值.

2.3 基于二元逆高斯過程的管道剩余壽命預(yù)測

建立剩余壁厚和剩余強(qiáng)度2個性能退化量的剩余壽命邊緣密度函數(shù)，利用Copula函數(shù)分析2個性能退化量的相關(guān)性，可獲得管道剩余壽命的聯(lián)合概率密度函數(shù)[23]，即

(8)

式中，f(x1,x2)為聯(lián)合概率密度函數(shù)；f1(x1)、f2(x2)分別為2個性能指標(biāo)對應(yīng)的概率密度函數(shù)；c(x)為Copula函數(shù)的密度函數(shù).

通過計算聯(lián)合概率密度函數(shù)，可得到管道失效概率及其剩余壽命.Copula函數(shù)類型眾多，常用的Copula函數(shù)有Gaussian、T、Gumbel、Frank四種類型.不同Copula函數(shù)會導(dǎo)致不同的分析結(jié)果，因此需結(jié)合工程實際選擇合適的Copula函數(shù).

AIC(Akaike information criterion)信息準(zhǔn)則可用于評價模型擬合數(shù)據(jù)的效果，具有廣泛的適用性.根據(jù)AIC信息準(zhǔn)則計算得到的R值越小說明擬合效果越好.本文采用AIC信息準(zhǔn)則選擇合適的Copula函數(shù)[24]：

R=-2lnA+2m

(9)

式中，A為模型對應(yīng)的似然函數(shù);m為模型中參數(shù)的個數(shù).

3 模型參數(shù)估計

假設(shè)Δyi=yi-yi-1對應(yīng)從時刻ti-1到ti的退化增量，Δti=ti-ti-1.在給定參數(shù)β的條件下，退化數(shù)據(jù)Yk的抽樣分布函數(shù)可以表示為[25]

(10)

(11)

其中

(12)

(13)

(14)

本文采用EM算法估計模型中的未知參數(shù)，基本步驟如下：① 對隱變量計算對數(shù)似然函數(shù)的期望，記為E-step；② 對求過期望后的對數(shù)似然函數(shù)最大化，記為M-step；③ 重復(fù)步驟①和②，直到滿足指定的收斂準(zhǔn)則.

(15)

(16)

(17)

(18)

(19)

η的初值可由下式得到：

(20)

(21)

式中，xi為第i個觀測時刻的退化值；n為觀測次數(shù).

4 案例分析

本節(jié)采用文獻(xiàn)[7]中某管道24 a的腐蝕深度檢測數(shù)據(jù)開展管道剩余壽命預(yù)測.所選管道材料為APILX52，管道內(nèi)壓為10 MPa，外徑為340.8 mm，管壁厚度為9.9 mm，最小極限應(yīng)力為359 MPa.表1為該管道不同檢測時刻(tc)的管道壁厚退化量.

表1 不同檢測時刻管道壁厚退化量

考慮到腐蝕管道失效與其剩余強(qiáng)度不足有關(guān)，采用有限元方法仿真求解管道剩余強(qiáng)度.對腐蝕缺陷建模時需要將缺陷簡化為規(guī)則形狀，本文將腐蝕缺陷假設(shè)為長方體[26].缺陷的幾何尺寸包括最大腐蝕長度(lmax)、寬度(bmax)和深度(dmax).上述尺寸變化會影響等效應(yīng)力和管道剩余強(qiáng)度.參照文獻(xiàn)[7]設(shè)置管道的材料屬性參數(shù)，分別建立9種不同的缺陷模型.通過有限元仿真模擬不同缺陷下管道最大等效應(yīng)力(τmax)，結(jié)果如表2所示.

表2 不同缺陷所對應(yīng)的等效應(yīng)力

由表2可知，當(dāng)腐蝕長度不變時，等效應(yīng)力隨著腐蝕缺陷深度的增加而增加；當(dāng)腐蝕深度不變時，等效應(yīng)力將隨缺陷長度或?qū)挾鹊脑黾佣黾?因此，管道剩余強(qiáng)度受不同方向腐蝕缺陷大小的綜合作用，不能簡單地忽略某個方向的缺陷尺寸.

建立管道腐蝕缺陷有限元模型，劃分網(wǎng)格并設(shè)定邊界條件，可以計算出管道退化過程中不同時刻的等效應(yīng)力.圖1所示為缺陷模型對應(yīng)的應(yīng)力云圖.根據(jù)式(7)，可計算得到不同檢測時刻的剩余強(qiáng)度退化量，如圖2所示.

圖1 管道應(yīng)力云圖

由式(8)可得到2個性能指標(biāo)剩余壽命的邊緣概率密度函數(shù)，并根據(jù)AIC信息準(zhǔn)則選擇合適的Copula函數(shù)管道剩余壽命聯(lián)合概率密度函數(shù).

圖2 管道剩余強(qiáng)度性能指標(biāo)退化趨勢

按照AIC信息準(zhǔn)則得到4種常用Copula函數(shù)的計算結(jié)果如表3所示.根據(jù)AIC準(zhǔn)則計算所得R值越小，說明所對應(yīng)的Copula函數(shù)適應(yīng)性越好.由表3可知，Gaussian Copula函數(shù)所對應(yīng)的R值最小.因此，選擇Gaussian Copula函數(shù)來表征2個性能指標(biāo)之間的相關(guān)性，得到管道剩余壽命的聯(lián)合概率密度函數(shù).

表3 4種Copula函數(shù)的R值

圖3為剩余壁厚指標(biāo)所對應(yīng)的不同檢測時刻對參數(shù)μ、σ2和η的估計值.管道剩余強(qiáng)度參數(shù)所對應(yīng)的參數(shù)估計與剩余壁厚退化指標(biāo)趨勢相似，此處不再贅述.由參數(shù)估計結(jié)果可知，隨著有效退化數(shù)據(jù)量的增加，各參數(shù)估計結(jié)果逐漸逼近真實參數(shù)值.原因在于，貝葉斯更新與EM算法結(jié)合能夠充分利用新獲得的退化數(shù)據(jù)，實時更新估計結(jié)果，有助于提高預(yù)測精度.上述結(jié)果驗證了本文所提出的參數(shù)估計方法的有效性.

(a) μ

(b) σ2

(c) η

表4 不同檢測時刻模型的參數(shù)估計值

圖4 不同時刻管道剩余壽命預(yù)測的概率密度函數(shù)

由圖4可知，隨著管道運行時間的增加，失效概率密度函數(shù)峰值呈增加趨勢，管道發(fā)生失效的概率不斷增大，同時管道的剩余壽命減小.上述結(jié)論與工程實際相符.

為客觀評價所提方法剩余壽命預(yù)測效果，定義相對誤差指標(biāo)(E)：

(22)

式中，S為剩余壽命預(yù)測值；T為實際的剩余壽命.

管道實際壽命采用文獻(xiàn)[7]中的腐蝕油氣管道的剩余壽命值，即設(shè)計壽命減去當(dāng)前已服役的時間.剩余壽命預(yù)測誤差結(jié)果如表5所示.基于二元IG過程的剩余壽命預(yù)測方法的最大預(yù)測誤差為10.7%，最小預(yù)測誤差為2.2%，總體上，采用IG過程的預(yù)測誤差相對較小，具有較高的預(yù)測精度.

表5 剩余壽命預(yù)測值與實際值比較

5 結(jié)論

1) 考慮到單一性能退化指標(biāo)難以全面反映輸油管道的健康狀態(tài)，文中提出同時利用2個性能指標(biāo)評估管道健康狀態(tài)的方法，并構(gòu)建基于二元IG過程的管道性能退化模型.

2) 結(jié)合數(shù)值算例，采用Gaussian Copula函數(shù)評估2個性能指標(biāo)的相關(guān)性；通過更新IG過程確定管道失效概率密度函數(shù)，分析管道健康狀態(tài)，完成管道剩余壽命預(yù)測.案例分析表明，基于二元IG過程的剩余壽命預(yù)測最大預(yù)測誤差為10.7%，最小預(yù)測誤差為2.2%，所提出的預(yù)測方法具有較高的預(yù)測精度.