基于數據融合和非線性維納過程的埋地管道退化過程預測

宋正濤，穆化巍，趙紅衛，崔 潔，韓小妹，梁昌晶

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

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

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

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

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

1 KPCA

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

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

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

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

2 退化建模

2.1 WP模型

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

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

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

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

2.2 加速退化軌道模型

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

式中：α、β為Arrhenius模型參數；Si為第i個應力。

2.3 管道剩余壽命預測及在線更新

2.3.1 剩余壽命預測

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

圖1 管道退化過程

式中：inf為函數下界。

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

式中：Φ()為標準正態函數。

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

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

2.3.2 在線參數更新

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

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

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

3 實例分析

3.1 數據收集

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

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

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

表1 管道外腐蝕數據集（部分）

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

圖3 退化過程正態性概率分布

3.2 主成分選取

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

圖4 核主成分碎石圖

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

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

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

圖5 主成分的特征向量

3.3 參數估計

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

表2 主成分應力作用下的參數值

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

3.4 剩余壽命概率密度預測

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

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

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

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

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

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

3.5 模型對比

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

表4 不同模型的預測方法對比

4 結論

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

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

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