宋正濤，穆化巍，趙紅衛(wèi)，崔 潔，韓小妹，梁昌晶

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

2.中國石油華北油田公司，河北任丘 062552

3.中國石油華北油田公司勘探開發(fā)研究院，河北任丘 062552

油氣管道是油氣資源開發(fā)和運輸?shù)闹匾浇椤ｋS著管道運行時間的增加以及復(fù)雜多變的環(huán)境和外界因素的作用，管道腐蝕、開裂、泄漏等事故時有發(fā)生，對人員生命安全和環(huán)境產(chǎn)生了一定影響[1?3]。對管道的剩余壽命進行預(yù)測是完整性管理的重要組成部分，可以有效減少事故發(fā)生，為管道維護、維修提供決策意見。

目前，對管道剩余壽命的研究主要從力學(xué)角度[4?5]、機器學(xué)習(xí)[6?7]和數(shù)理統(tǒng)計[8?10]等方面開展。力學(xué)角度：通過比較管道載荷和結(jié)構(gòu)抗力之間的關(guān)系，計算管道可靠性，但該方法是基于結(jié)構(gòu)本身特性進行計算，未考慮外部環(huán)境和不確定因素的影響，結(jié)果相對保守；機器學(xué)習(xí)：通過灰色模型、神經(jīng)網(wǎng)絡(luò)和支持向量機等實現(xiàn)已有壽命的參數(shù)擬合和回歸，但對于影響剩余壽命因素間的相關(guān)性考慮不足，對于管道長壽命、高可靠的產(chǎn)品預(yù)測效果并不好；數(shù)理統(tǒng)計：模型從單一的腐蝕深度或腐蝕速率數(shù)據(jù)出發(fā)，根據(jù)數(shù)據(jù)分布情況，預(yù)測管道最可能發(fā)生失效的年限，但未考慮數(shù)據(jù)波動和隨機退化對腐蝕趨勢的影響。此外，上述模型對于剩余壽命的概率密度和可靠度函數(shù)并未描述，無法獲得剩余壽命的眾數(shù)。維納過程（Wiener Process，WP）可以描述連續(xù)的性能退化過程，對于非線性、個體差異性和測量隨機性具有很好的適應(yīng)性[11]?；诖?，利用核主成分分析（Kernel Prin?cipal Component Analysis，KPCA）對影響管道外腐蝕的因素進行數(shù)據(jù)融合，結(jié)合非線性WP 和加速退化軌道模型建立管道剩余壽命的概率密度函數(shù)和可靠度函數(shù)；采用蒙特卡洛馬爾科夫鏈（Markov Chain Monte Carlo，MCMC）方法中的Gibbs 抽樣求解，進行模型參數(shù)估計，得到不同加速應(yīng)力下的參數(shù)估計值；通過實例分析，驗證對比模型結(jié)果和預(yù)測精度。

1 KPCA

KPCA 是Sch?lkoph 團隊在對主成分分析（PCA）進行非線性擴展的基礎(chǔ)上得到的[12]，通過特征值累計貢獻率的大小，衡量影響管道剩余壽命的關(guān)鍵因素。設(shè)有m個影響因素、n個樣本條目，則映射函數(shù)為：

式中：i為樣本編號；Rm、Rk分別為m維、k維向量；φ()為映射函數(shù)。

特征空間的協(xié)方差矩陣可以表示為：

通過核函數(shù)將式（2）化簡，得到特征值λ和特征向量，當(dāng)前q個特征值的累計貢獻率c大于85%時，則滿足要求，將數(shù)據(jù)重構(gòu)并進行反歸一化處理，得到降維樣本。

2 退化建模

2.1 WP模型

設(shè)X(t)為管道在t時刻的性能退化量，即腐蝕深度，當(dāng)隨機過程X(t)滿足增量正態(tài)性、增量獨立性和路徑連續(xù)性時，可以用WP 模型表示性能退化軌跡[13]，公式如下：

式中：X(0)為管道在初始時刻的腐蝕深度，即X(0)=0；W(t)為布朗運動，表示退化過程的時變性，且W(t)～N(0，t)；μ為漂移系數(shù)，用于表征管道的退化速率，即腐蝕速率；σ為擴散參數(shù)，用于表征管道退化軌跡的波動隨機過程。

為充分考慮管道的非線性退化過程，對式（4）進行改進，得：

式中：Λ(t∣b)為管道隨時間t的連續(xù)非線性退化函數(shù)；b為待定參數(shù)；ε為測量誤差?？紤]到不同管段位置的腐蝕差異性，令μ和ε相互獨立。

2.2 加速退化軌道模型

對于管道這類長壽命產(chǎn)品，在可見的壽命周期內(nèi)有可能出現(xiàn)零失效的現(xiàn)象，即從投產(chǎn)至腐蝕穿孔可能經(jīng)歷數(shù)年時間，無法在投產(chǎn)前的壽命試驗中進行可靠性評價，因此必須引入加速退化軌道模型用于描述失效機理不變條件下，管道壽命與應(yīng)力的關(guān)系。應(yīng)用較為廣泛的有Arrhenius 模型、Eyring 模型、Inverse Power 模型和Exponential 模型等[14]，考慮到漂移系數(shù)和擴散參數(shù)均與應(yīng)力相關(guān)，且管道腐蝕深度的變化多與時間呈冪次變化，則采用Arrhenius模型，公式如下：

式中：α、β為Arrhenius模型參數(shù)；Si為第i個應(yīng)力。

2.3 管道剩余壽命預(yù)測及在線更新

2.3.1 剩余壽命預(yù)測

由布朗運動特征可知，假設(shè)管道從t到t+Δt內(nèi)的壁厚減薄量ΔX是由無數(shù)個微小隨機退化量組成，且這些退化量彼此屬于獨立同分布，與時間呈正比，則ΔX服從正態(tài)分布。根據(jù)GB 50251—2015 中工藝計算的相關(guān)要求，在滿足腐蝕裕量的前提下，外徑與壁厚的比值不應(yīng)超過100，故以此為最大腐蝕深度的閾值l，X(t)首次達到l的時間即為首達時刻T（T為管道壽命），退化過程見圖1，公式如下：

圖1 管道退化過程

式中：inf為函數(shù)下界。

結(jié)合非線性WP 模型的性質(zhì)，T在進行函數(shù)變換后服從逆高斯分布，且漂移系數(shù)μ服從正態(tài)分布，通過推導(dǎo)得到管道剩余壽命的分布函數(shù)F(t)和可靠度函數(shù)R(t)分別為：

式中：Φ()為標(biāo)準(zhǔn)正態(tài)函數(shù)。

此時，管道的完全壽命期望值為E(T)=l/μ，概率密度函數(shù)f(t)為F(t)的積分，公式如下：

將KPCA 篩選后的數(shù)據(jù)用于描述管道腐蝕規(guī)律，建立腐蝕因素與腐蝕速率之間的加速關(guān)系，采用蒙特卡洛馬爾科夫鏈（MCMC） 方法中的Gibbs 抽樣求解式（4）～式（9）的模型參數(shù)，求解b、α、β和σ的估計值，再代入式（6）求μ值，最后通過式（10）求解剩余壽命的概率密度函數(shù)。

2.3.2 在線參數(shù)更新

根據(jù)管道現(xiàn)場監(jiān)測數(shù)據(jù)，實時更新漂移系數(shù)和擴散參數(shù)等模型參數(shù)，為突出不同管段退化過程的差異性，采用卡爾曼濾波原理對其進行實時更新。設(shè)X(t1)，X(t2)，…，X(tj)為t1，t2，…，tj時刻的退化量，退化過程如式（5）所示，則下一時刻的狀態(tài)轉(zhuǎn)移方程為：

式中：S0為初始時刻的應(yīng)力值；μj?1為上一時刻的漂移系數(shù)，μj為j時刻的漂移系數(shù)，兩者相等表示退化速率保持恒定，即為恒定應(yīng)力加速退化；k為漂移系數(shù)與擴散參數(shù)平方的比值。

在式（11）的基礎(chǔ)上，通過真實狀態(tài)預(yù)測、協(xié)方差預(yù)測、濾波增益、狀態(tài)更新、協(xié)方差更新等步驟，利用MCMC 方法更新參數(shù)估值，得到不同時間下的剩余壽命概率密度。

3 實例分析

3.1 數(shù)據(jù)收集

某輸氣管道跨越2 省4 市，全長521 km，管材為X80，管徑1 016 mm，壁厚10～12 mm，運行壓力（8±1）MPa。2015 年投產(chǎn)使用，沿線地質(zhì)環(huán)境復(fù)雜，因此管道的外腐蝕較嚴(yán)重。技術(shù)人員在2017年采用實地埋片的方式對同埋深的管道腐蝕速率進行了測定，試片沒有采取任何涂層、鍍膜、氧化等保護措施，將試片與管道用導(dǎo)線等電位連接，表示兩者在相同的陰極保護電位下。實驗時間30 d，腐蝕前后的部分試片的表面形貌見圖2。此外，對試片處的土壤成分進行實驗測定，包括電阻率、管地電位、含水量、硫化物含量等指標(biāo)。

圖2 腐蝕前后試片的表面形貌

考慮到土壤中的微生物類型主要為硫酸鹽還原菌、鐵細菌和腐生菌，他們均會在土壤中代謝產(chǎn)生硫化物，進而改變土壤氧含量、pH值等參數(shù)；土壤中的有機質(zhì)為含碳化合物，改善土壤通透性和透光性、提高保溫效果等均可改變土壤容重、孔隙度和質(zhì)地，進而影響含水量、含鹽量等參數(shù)，因此不再考慮以上間接因素對腐蝕的影響。為了便于處理，將Na+、K+、Ca2+、Mg2+、CO32?、HCO3?等離子的濃度統(tǒng)一作為含鹽量考慮，鑒于SO42?、Cl?會加快對管材的腐蝕作用，故單獨考慮。最終確定影響管道外腐蝕的12 項因素，形成的檢測數(shù)據(jù)集見表1。

表1 管道外腐蝕數(shù)據(jù)集（部分）

為了驗證退化模型的有效性，構(gòu)建退化數(shù)據(jù)的正態(tài)性概率分布，見圖3。其中，大部分數(shù)據(jù)均在對角線附近，通過AD、K?S 等統(tǒng)計量的假設(shè)檢驗，證明了管道壁厚的退化過程服從WP模型。

圖3 退化過程正態(tài)性概率分布

3.2 主成分選取

對表1 中的數(shù)據(jù)應(yīng)用式（1）～式（3），選擇徑向基核函數(shù)為KPCA 的核函數(shù)，得到特征值和特征向量，見圖4。

圖4 核主成分碎石圖

前5 個主成分的累計貢獻率已達到85.39%，且特征值均大于1，故可以選擇前5 個主成分代表之前的12個變量，重構(gòu)后的數(shù)據(jù)如式（12）：

式中：F1～F5為主成分，x1，x2，…，x12分別對應(yīng)表1表頭中的電阻率、氧化還原電位、…、陰極保護率等。

分析特征向量，絕對值越大，其因素越可以代表管道的外腐蝕特征，見圖5。以±0.5 為閾值，主成分F1上電阻率、含水量、pH值的代表性較強，F(xiàn)2上氯離子和硫酸根離子含量的代表性較強，F(xiàn)3上雜散電流的代表性較強，F(xiàn)4上陰極保護率的代表性較強，F(xiàn)5上氧化還原電位和自然電位的代表性較強。這5個主成分分別與土壤理化性質(zhì)、土壤加速腐蝕特性、雜散電流、陰極保護和管材理化特性等有關(guān)，基本涵蓋了埋地管道的外腐蝕影響因素，驗證了KPCA算法的科學(xué)性。

圖5 主成分的特征向量

3.3 參數(shù)估計

將KPCA 算法得到的數(shù)據(jù)進行非線性特征提取，采用處理后的數(shù)據(jù)作為應(yīng)力代入加速方程，隨后采用MCMC 方法進行參數(shù)估計。為確保參數(shù)收斂至全局最優(yōu)，有效提升參數(shù)估計的準(zhǔn)確性，預(yù)先設(shè)置b、α、β和σ的先驗分布均為Gamma 函數(shù)，位置、形狀、尺度分別為0.1、0.1 和3。以主成分F1對參數(shù)α的樣本路徑和累計均值為例，見圖6。在抽樣過程中，馬爾科夫鏈只在小范圍內(nèi)波動，說明建立的馬爾科夫鏈已通過多次轉(zhuǎn)移達到穩(wěn)態(tài)，同時累計均值也下降至平穩(wěn)狀態(tài)，95%的置信區(qū)間可以覆蓋對應(yīng)的參數(shù)取值范圍，說明抽樣誤差較小且迭代收斂良好。其他應(yīng)力下參數(shù)的迭代軌跡類同，結(jié)果見表2。

表2 主成分應(yīng)力作用下的參數(shù)值

圖6 參數(shù)α在主成分F1下的迭代軌跡

3.4 剩余壽命概率密度預(yù)測

取表2 中各參數(shù)的平均值作為終值，此時b=1.116，α=1.088，β=75.056，σ=0.446；將20 個樣本的數(shù)據(jù)代入式（12）中，取平均值作為終值，此時F1=60.135，F(xiàn)2=45.214，F(xiàn)3=92.452，F(xiàn)4=152.742，F(xiàn)5=26.11。將上述終值代入式（6）得到5 個主成分應(yīng)力作用下的μ1=0.312，μ2=0.207，μ3=0.483，μ4=0.665，μ5=0.061，均值μ=0.346。設(shè)初始腐蝕深度為0 mm，該管道的最大腐蝕深度閾值l=10 mm，得到管道的完全壽命期望值為E(T)=l/μ=10/0.346=28.9 a，與常規(guī)的設(shè)計年限相符（25～30 a）。代入式（10）得到第2 a時剩余壽命的概率密度函數(shù)：

剩余壽命的概率密度分布見圖7，峰值對應(yīng)的時刻即為檢測時間對應(yīng)的剩余壽命眾數(shù)，則剩余壽命為26.5 a，考慮到檢測時間為第2 a，則管道設(shè)計壽命為28.5 a，與之前計算的壽命期望值基本一致，說明了采用非線性WP 過程構(gòu)建的剩余壽命概率密度函數(shù)是正確的。

圖7 第2 a的管道剩余壽命概率密度分布

同理，根據(jù)卡爾曼濾波原理中的狀態(tài)轉(zhuǎn)移方程，得到不同檢測時間下的剩余壽命概率密度分布，見圖8。隨著管道運行時間的延長，剩余壽命的眾數(shù)不斷減小，分布形狀從“矮胖”向“高瘦”轉(zhuǎn)移，概率密度分布越來越集中，失效概率不斷增加，符合工程實際運行方式。對應(yīng)的設(shè)計壽命與壽命期望值見表3，兩者的相對誤差在[1.35%，3.43%]之間，且誤差在波動中逐漸減小，說明非線性WP過程的預(yù)測具有穩(wěn)定性和準(zhǔn)確性。

表3 基于非線性WP過程的壽命預(yù)測值

圖8 不同檢測時間的剩余壽命概率密度分布

3.5 模型對比

為了驗證主成分和非線性WP 模型的預(yù)測效果，將二元逆高斯模型[15]和線性WP 模型進行對比，采用有限元分析，結(jié)合Von Mises 等效應(yīng)力準(zhǔn)則，確定管道的實際剩余壽命，對比結(jié)果見表4。其中，二元逆高斯分布未考慮退化量之間的相關(guān)性，同時Copula 函數(shù)是對應(yīng)多元聯(lián)合分布的特殊函數(shù)，再用最大期望值法進行參數(shù)估計時容易陷入局部最優(yōu)解，導(dǎo)致誤差較大；而線性WP模型忽略了管道腐蝕的隨機性和不確定性，認為管道腐蝕深度隨時間延長而遞增，這與實際腐蝕深度在波動中遞增的結(jié)論不符。綜上所述，本文提出的模型對腐蝕管道剩余壽命的預(yù)測更接近實際值。

表4 不同模型的預(yù)測方法對比

4 結(jié)論

1）從管道腐蝕退化數(shù)據(jù)入手，基于主成分分析和非線性WP 過程，結(jié)合加速退化軌道模型，建立了用于描述腐蝕管道剩余壽命的概率函數(shù)，相較于其余模型，剩余壽命的預(yù)測誤差更小，說明本文所建模型用于預(yù)測外腐蝕影響下的剩余壽命是可行的。

2）經(jīng)KPCA 算法后，原有的12 個影響因素在非線性提取的過程中被降低至5個主成分，數(shù)據(jù)融合降低了因素間的多重共線性，無需求解線性優(yōu)化問題。

3）隨著管道運行時間的延長，剩余壽命的概率密度分布越來越集中，其剩余壽命不斷減小，只需不斷更新歷史數(shù)據(jù)即可實現(xiàn)長壽命、高可靠管道的可靠性評價。