基于数据融合和非线性维纳过程的埋地管道退化过程预测

2023-08-30宋正涛穆化巍赵红卫韩小妹梁昌晶

石油工程建设 2023年4期

宋正涛，穆化巍，赵红卫，崔洁，韩小妹，梁昌晶

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模型忽略了管道腐蚀的随机性和不确定性，认为管道腐蚀深度随时间延长而递增，这与实际腐蚀深度在波动中递增的结论不符。综上所述，本文提出的模型对腐蚀管道剩余寿命的预测更接近实际值。