CN109145399A - 一种基于改进的粒子滤波算法的疲劳裂纹扩展预测方法 - Google Patents

Abstract

本发明公开了一种基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，包括：A、定义状态模型与观测模型；B、进行模型参数转移；C、进行裂纹状态转移；D、当有新的裂纹监测值时，将粒子值带入观测似然概率密度中进行计算，得到粒子的归一化权值；得到裂纹长度的后验分布和模型参数的后验分布；E、将状态模型参数作为裂纹长度的扩展，得到新的裂纹长度与模型参数的粒子集；F、将裂纹长度与模型参数粒子集带入状态转移方程，实现裂纹发展趋势的预测，得到任意时刻下裂纹长度的概率分布；对于给定的裂纹长度阈值，计算得出任意时刻剩余寿命的概率分布。采用本发明，能够通过构造的新模型参数转移过程，提高参数的收敛速度，并提高预测精度。

Description

一种基于改进的粒子滤波算法的疲劳裂纹扩展预测方法

技术领域

本发明涉及故障预测与健康管理领域，尤其涉及一种基于改进的粒子滤波算法的疲劳裂纹扩展预测方法。

背景技术

状态预测是结构健康状态评估的重要环节，包括对结构损伤发展趋势和剩余寿命的预测。其中疲劳裂纹损伤是金属结构中最为常见的损伤形式之一，疲劳裂纹扩展的预测包括对裂纹发展趋势与剩余寿命的预测。目前主要采用基于断裂力学的疲劳裂纹扩展模型对其进行预测，如Paris模型、NASGRO模型等，其模型参数通常采用实验室数据拟合或根据经验得到。然而，由于疲劳裂纹扩展会受到多种不确定性因素的影响，如材料特性、载荷、环境因素以及其本身固有的随机性，即使相同结构在相同应力环境下其裂纹扩展也会存在差异，因此很难精确的估计疲劳裂纹长度与裂纹扩展模型参数。贝叶斯滤波技术可将物理模型与损伤监测数据相结合，对裂纹状态与模型参数进行准确的估计。

在贝叶斯滤波方法中，粒子滤波由于不严格要求线性高斯的条件，对于裂纹扩展与剩余寿命预测问题具有良好的适用性。粒子滤波方法可将模型参数作为状态矢量的一部分，构建一个扩展状态矢量，从而实现对裂纹长度与模型参数的联合估计。对于模型参数的状态转移，为了缓解粒子贫瘠的问题，现有的方法通常采用加上一个人工动态噪声的方法，而这种方法造成模型参数收敛较慢。

发明内容

有鉴于此，本发明的主要目的在于提供一种基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，通过构造新的模型参数转移过程，实现对疲劳剩余寿命更准确的估计，以提高参数的收敛速度，并提高预测精度。

为达到上述目的，本发明的技术方案如下：

一种基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，包括：

A、定义状态模型与观测模型的步骤；

B、进行状态模型参数转移的步骤；

C、进行裂纹状态转移的步骤；

D、当有新的裂纹监测值时，将每一粒子值带入观测似然概率密度中进行计算，将其作为粒子权值，然后对所有粒子权值进行归一化处理，得到粒子的归一化权值；并通过裂纹长度粒子集及相应的归一化权值表示裂纹长度的后验分布，再通过模型参数粒子集及相应的归一化权值表示模型参数的后验分布；

E、将状态模型参数作为裂纹长度的扩展，裂纹长度与状态模型参数构成一个扩展状态变量，根据每个粒子的归一化权值，采用多项式重采样法进行重采样，得到新的裂纹长度与模型参数的粒子集，其中每一粒子的权值均变为粒子数之一；

F、将裂纹长度与模型参数粒子集带入状态转移方程，实现裂纹发展趋势的预测，得到任意时刻下裂纹长度的概率分布；

G、对于给定的裂纹长度阈值，计算得出任意时刻剩余寿命的概率分布。

其中，步骤D还包括：若无新的裂纹监测，则返回步骤C。

较佳地，进一步包括步骤H，即返回执行步骤B到步骤G。

其中，步骤A所述定义状态模型与观测模型，具体为：

其中，k为指第k个离散的时刻，x_k为系统的状态矢量，x_k＝f(x_k-1，ω_k)为状态模型，ω_k为状态转移噪声，该状态模型定义了状态转移概率p(x_k|x_k-1)，因此也称为状态转移方程；z_k为状态的观测矢量，z_k＝g(x_k,ν_k)为观测模型，ν_k为观测噪声，该观测模型定义了似然概率p(z_k|x_k)。

步骤A具体包括：

A1、选择Paris模型作为疲劳裂纹扩展模型，并据此定义裂纹的状态转移方程即状态转移概率密度：

其中：k表示离散的时刻x_k是k时刻的裂纹长度，x_k-1是k-1时刻的裂纹长度，lnC与m为与材料、环境因素相关的状态模型参数，ΔN表示两个相邻时刻之间的循环周次间隔，ΔK表示应力强度因子，ω_k-1为k-1时刻的状态模型过程噪声，其服从均值为方差为的正态分布，p(x_k|x_k-1)为状态转移概率密度；

A2、将模型参数lnC作为待估参数，与裂纹长度共同构成扩展的状态矢量：

由于在Paris模型中，参数lnC为定参数，不随时间变化，则lnC的转移状态方程表示为：

A3、定义观测模型为：

z_k＝x_k+ν_k，

其中，z_k为状态的观测矢量，x_k为系统的状态矢量，ν_k为观测噪声，N(0,δ_v ²)表示ν_k服从均值为0、方差为的正态分布。

其中，步骤B具体包括：

B1：求系统状态转移方程对模型参数的导数：

B2：计算k时刻模型参数lnC各个粒子的位置参数：

其中，α∈[0,1]为参数位置调整因子，且

B3、计算k时刻模型参数转移噪声的方差：

B4：生成新的模型参数粒子群：

其中，表示k时刻模型参数lnC各个粒子的位置参数，表示k时刻模型参数转移噪声的方差，表示一个由均值为0、方差为的正态分布生成的随机数。

步骤C具体包括：由步骤A定义的状态转移概率密度p(x_k|x_k-1)中采样，即将步骤B得到的新的参数粒子集与上一时刻的裂纹长度粒子集带入状态转移方程，得到下一时刻的裂纹状态粒子集：

其中，k表示离散的时刻，i表示粒子标号，表示k时刻第i个粒子的裂纹长度，是k-1时刻第i个粒子的裂纹长度，表示第i个粒子的状态模型参数，ΔN表示两个相邻时刻之间的循环周次间隔，ΔK表示应力强度因子。

步骤D具体包括：当有新的裂纹监测值时，根据下式对粒子权值进行更新：

然后对所有粒子权值进行归一化处理，得到粒子的归一化权值：

则裂纹长度与模型参数的后验分布分别为：

其中，δ(·)为狄拉克函数，为第i个粒子的归一化权值，与分别表示第i个粒子的裂纹长度与状态模型参数值。

步骤E具体包括：根据归一化权值对所有粒子进行重采样，重采样方法选用多项式重采样，得到新的粒子集及其相应权值：

其中，N_s表示粒子数目，与分别表示第i个粒子重采样后新生成的裂纹长度与状态模型参数值。

步骤F具体包括：将裂纹长度与模型参数的后验粒子集带入裂纹长度的状态转移方程中，对裂纹扩展趋势进行预测，得到未来任意时刻的裂纹长度的概率分布：

其中，x_k+p表示基于k时刻裂纹长度与状态模型参数的后验概率分布，预测的k+p时刻的裂纹长度。

步骤G具体包括：对于给定的裂纹阈值x_th，对剩余寿命进行预测，剩余寿命的概率分布为：

则，

剩余寿命的点估计值为：

其中，RUL_k表示k时刻的剩余寿命，表示k时刻剩余寿命的期望。

本发明基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，具有如下有益效果：

1)对于模型参数的状态转移，为了缓解粒子贫瘠的问题，标准粒子滤波算法通常采用加上一个人工动态噪声的方法，而这种参数转移方法引入了额外的不确定性，在参数初值与真实值相差较远的情况下，会导致模型参数收敛较慢。本发明通过在模型参数的转移过程中，引入观测值的信息，提出了一种改进的粒子滤波算法，与标准粒子滤波算法相比，提高了模型参数的收敛速度，可提供更加准确的早期预测。

2)由于在参数转移过程中引入了观测值的信息，该方法与标准粒子滤波算法相比，缓解了粒子退化程度，提高了剩余寿命的预测精度。

附图说明

图1为本发明基于改进的粒子滤波算法的疲劳裂纹扩展预测方法流程图；

图2为仿真生成的裂纹真实值与裂纹监测值随循环周次的变化示意图；

图3为模型参数lnC随滤波时间的变化示意图；

图4为裂纹长度随循环周次的变化示意图；

图5为剩余寿命随循环周次的变化示意图；

图6为剩余寿命估计精度随循环周次的变化示意图。

具体实施方式

下面结合附图及本发明的实施例对本发明的方法作进一步详细的说明。

图1为本发明基于改进的粒子滤波算法的疲劳裂纹扩展预测方法流程图。

如图1所示，该基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，具体包括如下步骤：

步骤11：定义状态模型与观测模型的步骤。

其中，k为指第k个离散的时刻，x_k为系统的状态矢量，x_k＝f(x_k-1，ω_k)为状态模型，ω_k为状态转移噪声，该状态模型定义了状态转移概率p(x_k|x_k-1)，因此也称为状态转移方程；z_k为状态的观测矢量，z_k＝g(x_k,ν_k)为观测模型，ν_k为观测噪声，该观测模型定义了似然概率p(z_k|x_k)。

在具体实施例中，具体包括：

步骤111：选择Paris模型作为疲劳裂纹扩展模型，并据此定义裂纹的状态转移方程(状态转移概率密度)：

上式中，k表示离散的时刻x_k是k时刻的裂纹长度，x_k-1是k-1时刻的裂纹长度，lnC与m为与材料、环境因素相关的状态模型参数，ΔN是两个相邻时刻之间的循环周次间隔，ΔK是应力强度因子，ω_k-1为k-1时刻的状态模型过程噪声，其服从均值为方差为的正态分布，p(x_k|x_k-1)为状态转移概率密度。

步骤112：将模型参数lnC作为待估参数，与裂纹长度共同构成扩展的状态矢量：

由于在Paris模型中，参数lnC为定参数，不随时间变化，则lnC的转移状态方程可写为：

步骤113：定义观测模型为：

z_k＝x_k+ν_k

其中，z_k为状态的观测矢量，x_k为系统的状态矢量，ν_k为观测噪声，N(0,δ_v ²)表示ν_k服从均值为0、方差为的正态分布。

步骤12：进行状态模型参数转移的步骤。

首先，求取状态转移方程对状态模型参数的导数；然后，分别计算状态模型参数的位置参数与转移噪声方差；最后，由以位置参数为均值、转移噪声方差为方差的正态分布中抽样，得到新的参数粒子集。

在具体实施例中，将状态模型参数转移，生成新的参数粒子集。其步骤包括：

步骤121：求状态模型对状态模型参数的导数：

步骤122：计算k时刻状态模型参数lnC各个粒子的位置参数：

其中，α∈[0,1]为参数位置调整因子，且

步骤123：计算k时刻状态模型参数转移噪声的方差：

步骤124：生成新的状态模型参数粒子群：

进一步地，上述状态模型参数的建议分布可写为：

其中，表示k时刻模型参数lnC各个粒子的位置参数，表示k时刻模型参数转移噪声的方差，表示一个由均值为0、方差为的正态分布生成的随机数。

在本发明的实施例中，步骤12中在计算状态模型参数的位置参数时，参数位置调整因子的最优值为0.96。

步骤13：进行裂纹状态转移的步骤。

由步骤11中定义的状态转移概率密度p(x_k|x_k-1)中采样，即将步骤12中得到的新的参数粒子集与上一时刻的裂纹长度粒子集，带入状态转移方程，得到下一时刻的裂纹状态粒子集。

在具体实施例中，具体为：由步骤11中定义的状态转移概率密度p(x_k|x_k-1)中采样，即将步骤12得到的新的参数粒子集与上一时刻的裂纹长度粒子集，带入状态转移方程，得到下一时刻的裂纹状态粒子集：

其中，k表示离散的时刻，i表示粒子标号，表示k时刻第i个粒子的裂纹长度，是k-1时刻第i个粒子的裂纹长度，表示第i个粒子的状态模型参数，ΔN表示两个相邻时刻之间的循环周次间隔，ΔK表示应力强度因子。

步骤14：当有新的裂纹监测值时，将每一粒子值带入观测似然概率密度中进行计算，将其作为粒子权值，之后对所有粒子权值进行归一化处理，得到粒子的归一化权值。通过裂纹长度粒子集及相应的归一化权值表示裂纹长度的后验分布，通过模型参数粒子集及相应的归一化权值表示模型参数的后验分布。若无新的裂纹监测，则返回步骤13。

在具体实施例中，当有新的裂纹监测值时，根据下式对粒子权值进行更新：

之后，对所有粒子权值进行归一化处理，得到粒子的归一化权值：

则裂纹长度与模型参数的后验分布分别为：

若无新的裂纹监测值，返回步骤13。

其中，δ(·)为狄拉克函数，为第i个粒子的归一化权值，与分别表示第i个粒子的裂纹长度与状态模型参数值。

步骤15：将模型参数作为裂纹长度的扩展，两者构成一个扩展状态变量，根据每个粒子的归一化权值，采用多项式重采样方法进行重采样，得到新的裂纹长度与模型参数的粒子集，其中每一粒子的权值均变为粒子数之一。

在具体实施例中，根据归一化权值对所有粒子进行重采样，重采样方法选用多项式重采样，得到新的粒子集及其相应权值：

其中，N_s表示粒子数目，与分别表示第i个粒子重采样后新生成的裂纹长度与状态模型参数值。

步骤16：将裂纹长度与模型参数粒子集带入状态转移方程可实现裂纹发展趋势的预测，得到任意时刻下裂纹长度的概率分布。

在具体实施例中，将裂纹长度与模型参数的后验粒子集带入裂纹长度的状态转移方程中，对裂纹扩展趋势进行预测，得到未来任意时刻的裂纹长度的概率分布：

其中，x_k+p表示基于k时刻裂纹长度与状态模型参数的后验概率分布，预测的k+p时刻的裂纹长度。

步骤17：对于给定的裂纹长度阈值，可得到任意时刻剩余寿命的概率分布。

在具体实施例中，对于给定的裂纹阈值x_th，可对剩余寿命进行预测，剩余寿命的概率分布为：

剩余寿命的点估计值为：

其中，RUL_k表示k时刻的剩余寿命，表示k时刻剩余寿命的期望。

进一步的，还可以包括：

步骤18：返回执行步骤12到步骤17。

其中，本说明书上下文中所列公式中的参数的物理意义或含义，前后均相同或一致，故不再赘述。

下面以仿真的铝合金裂纹扩展为例来具体说明本发明上述方法的具体实验数据和有益效果。

图2为仿真生成的裂纹真实值与裂纹监测值随循环周次的变化示意图。

如图2所示为仿真的裂纹扩展数据与裂纹长度监测数据，相关参数设置如表1所示。

表1：仿真参数设置

之后，对改进的粒子滤波算法的进行初始化，相关参数如表2所示。

图3为模型参数lnC随滤波时间的变化示意图；图4为裂纹长度随循环周次的变化示意图；图5为剩余寿命随循环周次的变化示意图。

滤波结果如图3、图4、图5所示，从图3～图5中可以看出，经过少量迭代，模型参数与剩余寿命就迅速收敛于真实值。

图6为剩余寿命估计精度随循环周次的变化示意图。

采用如下精度指标评价剩余寿命的估计精度：

其中，RUL_k表示k时刻的剩余寿命，表示k时刻剩余寿命的期望。

如图6所示，显示了剩余寿命估计精度随循环周次的变化过程。

由此可见，采用上述的仿真实验验证了本发明方法的有效性，该方法能够适用于对含有裂纹损伤的结构进行实时的损伤发展趋势和剩余寿命预测。

以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (10)

1.一种基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，包括：

A、定义状态模型与观测模型的步骤；

B、进行状态模型参数转移的步骤；

C、进行裂纹状态转移的步骤；

D、当有新的裂纹监测值时，将每一粒子值带入观测似然概率密度中进行计算，将其作为粒子权值，然后对所有粒子权值进行归一化处理，得到粒子的归一化权值；并通过裂纹长度粒子集及相应的归一化权值表示裂纹长度的后验分布，再通过模型参数粒子集及相应的归一化权值表示模型参数的后验分布；

E、将状态模型参数作为裂纹长度的扩展，裂纹长度与状态模型参数构成一个扩展状态变量，根据每个粒子的归一化权值，采用多项式重采样法进行重采样，得到新的裂纹长度与模型参数的粒子集，其中每一粒子的权值均变为粒子数之一；

F、将裂纹长度与模型参数粒子集带入状态转移方程，实现裂纹发展趋势的预测，得到任意时刻下裂纹长度的概率分布；

G、对于给定的裂纹长度阈值，计算得出任意时刻剩余寿命的概率分布。

2.根据权利要求1所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，步骤D还包括：若无新的裂纹监测，则返回步骤C。

3.根据权利要求1所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，进一步包括步骤H，即返回执行步骤B到步骤G。

4.根据权利要求1所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，步骤A所述定义状态模型与观测模型，具体为：

其中，k为指第k个离散的时刻，x_k为系统的状态矢量，x_k＝f(x_k-1，ω_k)为状态模型，ω_k为状态转移噪声，该状态模型定义了状态转移概率p(x_k|x_k-1)即状态转移方程；z_k为状态的观测矢量，z_k＝g(x_k,ν_k)为观测模型，ν_k为观测噪声，该观测模型定义了似然概率p(z_k|x_k)。

5.根据权利要求4所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，具体包括：

A1、选择Paris模型作为疲劳裂纹扩展模型，并据此定义裂纹的状态转移方程即状态转移概率密度：

其中：k表示离散的时刻x_k是k时刻的裂纹长度，x_k-1是k-1时刻的裂纹长度，lnC与m为与材料、环境因素相关的状态模型参数，ΔN表示两个相邻时刻之间的循环周次间隔，ΔK表示应力强度因子，ω_k-1为k-1时刻的状态模型过程噪声，其服从均值为方差为的正态分布，p(x_k|x_k-1)为状态转移概率密度；

A2、将模型参数lnC作为待估参数，与裂纹长度共同构成扩展的状态矢量：

由于在Paris模型中，参数lnC为定参数，不随时间变化，则lnC的转移状态方程表示为：

A3、定义观测模型为：

z_k＝x_k+ν_k，

其中，z_k为状态的观测矢量，x_k为系统的状态矢量，ν_k为观测噪声，N(0,δ_v ²)表示ν_k服从均值为0、方差为的正态分布。

6.根据权利要求5所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，步骤B具体包括：

B1：求系统状态转移方程对模型参数的导数：

B2：计算k时刻模型参数lnC各个粒子的位置参数：

其中，α∈[0,1]为参数位置调整因子，且

B3、计算k时刻模型参数转移噪声的方差：

B4：生成新的模型参数粒子群：

其中，表示k时刻模型参数lnC各个粒子的位置参数，表示k时刻模型参数转移噪声的方差，表示一个由均值为0、方差为的正态分布生成的随机数。

7.根据权利要求6所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，步骤C具体包括：

由步骤A定义的状态转移概率密度p(x_k|x_k-1)中采样，即将步骤B得到的新的参数粒子集与上一时刻的裂纹长度粒子集带入状态转移方程，得到下一时刻的裂纹状态粒子集：

其中，k表示离散的时刻，i表示粒子标号，表示k时刻第i个粒子的裂纹长度，是k-1时刻第i个粒子的裂纹长度，表示第i个粒子的状态模型参数，ΔN表示两个相邻时刻之间的循环周次间隔，ΔK表示应力强度因子。

8.根据权利要求7所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，步骤D具体包括：

当有新的裂纹监测值时，根据下式对粒子权值进行更新：

然后对所有粒子权值进行归一化处理，得到粒子的归一化权值：

则裂纹长度与模型参数的后验分布分别为：

其中，δ(·)为狄拉克函数，为第i个粒子的归一化权值，与分别表示第i个粒子的裂纹长度与状态模型参数值。

9.根据权利要求8所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，步骤E具体包括：

根据归一化权值对所有粒子进行重采样，重采样方法选用多项式重采样，得到新的粒子集及其相应权值：

其中，N_s表示粒子数目，与分别表示第i个粒子重采样后新生成的裂纹长度与状态模型参数值。

10.根据权利要求9所述基于改进的粒子滤波算法的疲劳裂纹扩展预测方法，其特征在于，步骤F具体包括：

将裂纹长度与模型参数的后验粒子集带入裂纹长度的状态转移方程中，对裂纹扩展趋势进行预测，得到未来任意时刻的裂纹长度的概率分布：

其中，x_k+p表示基于k时刻裂纹长度与状态模型参数的后验概率分布，预测的k+p时刻的裂纹长度。