CN110532620A - 一种基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法 - Google Patents

Abstract

本发明公开了一种基于递推最小二乘‑核平滑粒子滤波的疲劳裂纹扩展预测方法，包括：A、定义状态模型与观测模型；B、采用递推最小二乘‑核平滑方法进行状态模型参数转移；C、裂纹状态转移；D、当有新的裂纹监测值时，计算粒子权值；并通过裂纹长度粒子集、模型参数及相应的归一化权值表示两者的后验分布；E、将状态模型参数作为裂纹长度的扩展，得到新的裂纹长度与模型参数的粒子集；F、将裂纹长度与模型参数粒子集带入状态转移方程，实现裂纹发展趋势的预测，得到任意时刻下裂纹长度的概率分布；G、对于给定的裂纹长度阈值，计算得出任意时刻剩余寿命的概率分布。本发明的方法提高了对裂纹扩展剩余寿命的预测性能。

Description

一种基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预 测方法

技术领域

本发明涉及故障预测与健康管理领域，具体而言涉及一种基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法。

背景技术

疲劳裂纹是金属结构中最为常见的损伤形式之一，它将会造成结构退化甚至失效。疲劳裂纹扩展的准确预测对保证结构安全、指定合理的维修计划至关重要。现已有多种疲劳裂纹扩展模型可用于预测，当选定模型后，当前裂纹状态和模型参数直接决定了预测结果的准确性。随着目前结构健康监测技术的发展，粒子滤波方法逐渐应用于疲劳裂纹扩展的预测。粒子滤波基于贝叶斯理论和重要性采样方法，将物理模型与监测数据相结合，可实现对裂纹状态与模型参数的联合估计，由于其不存在对非线性非高斯的严格假设，对疲劳裂纹扩展预测问题具有良好的适用性。

在对裂纹状态和模型参数进行联合估计时，需要分别定义裂纹状态和模型参数的转移过程。裂纹状态的转移通常根据经典裂纹扩展模型Paris公式得到。Paris公式中存在两个需要估计的模型参数ln C和m，其转移通常是基于人工演化(Artificial Evolution,AE)或核平滑(Kernel Smoothing,KS)的方法。然而我们发现，这种参数转移方法在参数先验与参数真实值较为接近的情况下可以获得准确的参数估计和裂纹扩展预测结果，但当参数先与真实值相差较远时，则会出现一个“协调变化”现象，即裂纹扩展可得到正确的预测结果，但参数却没有收敛至其真实值。理论上应该只有唯一的一对ln C和m，与材料和所处环境有关，相反我们发现似乎存在多个组合均能正确的预测裂纹扩展轨迹，在滤波过程中，ln C和m自动向距其较近的一对组合移动直至收敛。在这种情况下，估计所得到的ln C和m不具备任何的实际意义，应用它们具备很大的风险，难以指导进行后续的维修和改进。

发明内容

有鉴于此，本发明的主要目的在于提供一种基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法。在粒子滤波框架内，提出一种递推最小二乘-核平滑方法，重点在于构造Paris公式模型参数的转移过程，获得对模型参数的正确估计，进而实现对疲劳裂纹扩展的准确预测。

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

一种基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，包括以下步骤：

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

B、采用递推最小二乘-核平滑方法进行状态模型参数转移；

C、进行裂纹状态转移；

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

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

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

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

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

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

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

x_k＝f(x_k-1，ω_k) (1)

z_k＝g(x_k,ν_k) (2)

其中，k指离散的时刻，x_k为k时刻系统的状态矢量，x_k＝f(x_k-1，ω_k)为状态模型，ω_k为状态转移噪声；z_k为状态的观测矢量，z_k＝g(x_k,ν_k)为观测模型，ν_k为观测噪声。

其中，步骤A具体包括：

A1、定义状态模型：以Paris模型作为疲劳裂纹扩展模型定义状态模型，其参数lnC和m为待估参数，与裂纹长度共同构成扩展的状态矢量：

其中：x_k是k时刻的裂纹长度，θ_k为k时刻的参数矢量，lnC_k与m_k为k时刻的Paris模型参数。据此定义状态模型为：

其中：ΔN表示两个相邻时刻之间的循环周次间隔，F表示形状因子，Δσ表示应力范围，ω_k-1为k-1时刻的状态模型过程噪声，表示ω_k-1服从均值为方差为的正态分布；

A2、定义观测模型为：

其中，z_k为裂纹观测值，ν_k为观测噪声，表示ν_k服从均值为0、方差为的正态分布。

其中，步骤B具体包括：

B1：对状态模型左右两侧分别取期望值，并进行双对数变换：

其中：为k时刻裂纹长度粒子的期望值，和分别为k-1时刻裂纹长度和模型参数粒子的期望值；

假设在k时刻的观测误差为零，则有等于z_k，可将上式写为下面的形式：

式中：θ_l,k-1为经过双对数变换后的新的参数矢量，为输入矢量，y_k为输出变量，三者表达式如下：

B2：采用带遗忘因子的最小二乘方法对参数矢量θ_l,k进行估计：

式中：G_k为增益；P_k-1为计算中的过程矩阵，其初始值设置为P₀＝10⁶I，I为二阶单位矩阵；λ为遗忘因子，其值范围通常设置为[0.95,1]，得到θ_l,k，可通过反算得到k时刻经过带遗忘因子的递推最小二乘法估计得到的参数矢量

式中：θ_l,k(1)和θ_l,k(2)分别表示矢量θ_l,k的第一个元素和第二个元素；

B3：对参数粒子进行收缩，计算每个参数粒子的核位置矢量：

式中，i为粒子标号，为k时刻第i个参数粒子的核位置矢量，为k-1时刻第i个参数粒子矢量，α_k为核参数，计算如下：

式中，α₀为核参数的初始值，一般设置为[0.95,1]；

B4：对参数粒子进行高斯噪声抖动，

首先计算高斯噪声的方差：

式中：为高斯噪声的方差矩阵，表示的方差矩阵；

然后计算转移后新的模型参数粒子：

其中，表示一个由均值为0、方差为的二维正态分布生成的随机数。

步骤C具体包括：根据步骤A定义的状态模型，将步骤B得到的新的参数粒子集与上一时刻的裂纹长度粒子集带入状态方程，得到下一时刻的裂纹状态粒子集：

其中，表示k时刻第i个粒子的裂纹长度，是k-1时刻第i个粒子的裂纹长度，和表示第i个粒子的模型参数，ΔN表示两个相邻时刻之间的循环周次间隔，F表示形状因子，Δσ表示应力范围。

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

式中：为k时刻第i个粒子的权值，p(|)表示条件概率。然后对所有粒子权值进行归一化处理，得到粒子的归一化权值：

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

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

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

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

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

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

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

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

本发明基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，具有如下有益效果：

1)针对现有方法中出现的“协调变化”问题，本发明通过将递推最小二乘和核平滑的方法相结合，设计出一种新的状态模型参数转移方法。该方法解决了当参数先验与真实值相距较远情况下，ln C和m之间出现的“协调变换”问题，实现了对模型参数的准确估计。

2)不论是在参数先验与真实值相距较近还是相距较远的情况，与现有方法相比，所提出的方法均提高了对裂纹扩展剩余寿命的预测性能。

附图说明

图1为本发明基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法流程图；

图2为试样尺寸图和实际图；

图3为实际裂纹扩展轨迹和相应的裂纹监测结果图；

图4为比较例一中的典型裂纹扩展预测结果；

图5为比较例一中的典型参数估计结果；

图6为比较例二中的典型裂纹扩展预测结果；

图7为比较例二中的典型参数估计结果。

具体实施方式

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

图1为本发明基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法流程图。

如图1所示，该基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，具体包括如下步骤：

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

x_k＝f(x_k-1，ω_k) (1)

z_k＝g(x_k,ν_k) (2)

其中，k指离散的时刻，x_k为k时刻系统的状态矢量，x_k＝f(x_k-1，ω_k)为状态模型，ω_k为状态转移噪声；z_k为状态的观测矢量，z_k＝g(x_k,ν_k)为观测模型，ν_k为观测噪声。

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

步骤111：定义状态模型：以Paris模型作为疲劳裂纹扩展模型定义状态模型，其参数ln C和m为待估参数，与裂纹长度共同构成扩展的状态矢量：

其中：x_k是k时刻的裂纹长度，θ_k为k时刻的参数矢量，lnC_k与m_k为k时刻的Paris模型参数，据此定义状态模型为：

其中：ΔN表示两个相邻时刻之间的循环周次间隔，F表示形状因子，Δσ表示应力范围，ω_k-1为k-1时刻的状态模型过程噪声，表示ω_k-1服从均值为方差为的正态分布。

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

其中，z_k为裂纹观测值，ν_k为观测噪声，表示ν_k服从均值为0、方差为的正态分布。

步骤12：采用递推最小二乘-核平滑方法进行状态模型参数转移的步骤。

首先，对状态模型进行双对数变换为线性形式；然后采用带遗忘因子的最小二乘方法对参数进行估计；之后，对参数粒子进行收缩，计算每个参数粒子的核位置矢量；最后，对参数粒子进行高斯噪声抖动，得到转移后新的参数粒子集。

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

步骤121：对状态模型左右两侧分别取期望值，并进行双对数变换：

其中：为k时刻裂纹长度粒子的期望值，和分别为k-1时刻裂纹长度和模型参数粒子的期望值。假设在k时刻的观测误差为零，则有等于z_k，可将上式写为下面的形式：

式中：θ_l,k-1为经过双对数变换后的新的参数矢量，为输入矢量，y_k为输出变量，三者表达式如下：

步骤122：采用带遗忘因子的最小二乘方法对参数矢量θ_l,k进行估计：

式中：G_k为增益；P_k-1为计算中的过程矩阵，其初始值设置为P₀＝10⁶I，I为二阶单位矩阵；λ为遗忘因子，其值范围通常设置为[0.95,1]，得到θ_l,k，可通过反算得到k时刻经过带遗忘因子的递推最小二乘法估计得到的参数矢量

式中：θ_l,k(1)和θ_l,k(2)分别表示矢量θ_l,k的第一个元素和第二个元素。

步骤123：对参数粒子进行收缩，计算每个参数粒子的核位置矢量：

式中，i为粒子标号，为k时刻第i个参数粒子的核位置矢量，为k-1时刻第i个参数粒子矢量，α_k为核参数，计算如下：

式中，α₀为核参数的初始值，一般设置为[0.95,1]。

步骤124：对参数粒子进行高斯噪声抖动，

首先计算高斯噪声的方差：

式中：为高斯噪声的方差矩阵，表示的方差矩阵；

然后计算转移后新的模型参数粒子：

其中，表示一个由均值为0、方差为的二维正态分布生成的随机数。

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

由步骤11定义的状态模型，将步骤12得到的新的参数粒子集与上一时刻的裂纹长度粒子集带入状态转移方程，得到下一时刻的裂纹状态粒子集。

在具体实施例中，具体为：由步骤11定义的状态模型，将步骤12得到的新的参数粒子集与上一时刻的裂纹长度粒子集带入状态方程，得到下一时刻的裂纹状态粒子集：

其中，表示k时刻第i个粒子的裂纹长度，是k-1时刻第i个粒子的裂纹长度，和表示第i个粒子的模型参数，ΔN表示两个相邻时刻之间的循环周次间隔，F表示形状因子，Δσ表示应力范围。

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

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

式中：为k时刻第i个粒子的权值，p(|)表示条件概率。然后对所有粒子权值进行归一化处理，得到粒子的归一化权值：

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

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

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

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

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

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

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

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

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

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

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

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

进一步的，还可以包括：

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

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

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

图2为试样尺寸图和实际图，采用压电传感器对裂纹进行了监测，图3为实际裂纹扩展轨迹和相应的裂纹监测结果。

试验验证采用两个比较例进行。在比较例一中参数先验与其真实值相距较远，在比较例二中参数先验与其真实值相距较近，其他参数设置相同。在每一个比较例中，均分别采用现有的AE方法、KS方法和本发明所提出的递推最小二乘-核平滑(recursive leastsquares-kernel smoothing,RLS-KS)方法(下文称为RLS-KS方法)进行状态模型参数的转移，将三种方法的结果对比分析，从而来证明本发明方法的创新性。粒子滤波算法初始化设置如表1所示。与状态模型参数转移方法相关的参数设置单独说明如下：在AE方法中，logC和m的人工噪声方差分别为0.005和0.001；在KS方法中，核参数为0.1；在本发明的RLS-KS方法中，遗忘因子和核参数初始值分别为0.998和0.99。为了避免抽样误差的影响，每种情况下的滤波均进行100次。

表1粒子滤波算法初始化设置

为了定量评价方法的性能，进行了指标的定义。好的性能通常意味着具有以下几个特点：1)真实值可落入预测的95％置信区间内；2)均值估计与真实值接近；3)95％置信区间的宽度较小。据此，我们定义了以下三个指标：

(1)命中率：

式中：N为总滤波次数，在比较例中为100；n为真实值落入95％置信区间内的次数。

(2)相对误差

式中：为均值，X_r为真实值

(3)相对置信区间宽度

式中：X_u为95％置信区间上限，X_l为95％置信区间下限。

显然，对于以上三个指标，HR越大，RE和RIW越小，意味着性能越好。

图4和图5分别为比较例一中三种方法的典型裂纹扩展预测和参数估计结果。其中，图4中左侧图为裂纹扩展轨迹的预测结果，右侧为在给定阈值下的剩余寿命概率分布。由图4可以看出，三种方法均可提供准确的裂纹扩展预测结果，真实值均落入预测的95％置信区间内。与AE方法和KS方法相比，本发明的RLS-KS方法具有集中的后验分布，且裂纹扩展轨迹和剩余寿命的均值与真实值更为接近。因此，本发明的RLS-KS方法具有相对更好的裂纹扩展预测性能。由图5可以看出，三种方法均可正确的预测模型参数log C和m，大约在70000循环周次左右，log C和m均收敛于其真实值并在之后保持稳定。表2给出了比较例一中三种方法的指标对比。由表2可以看出，三种方法中剩余寿命的HR均为100％，即预测的95％置信区间均将真实值包含在内。参数log C和m的HR在95％-100％之间，说明有几次遗漏，但这是由于抽样误差的存在因此是可以接受的。此外，对于参数log C和m，本发明的RLS-KS方法的RE和RIW均稍大于AE方法和KS方法，但对于剩余寿命，本发明的RLS-KS方法的RE和RIW却明显小于AE方法和KS方法。也就是说，虽然本发明的RLS-KS方法对于参数估计的性能稍差，但预测性能却明显提高。

表2比较例一中三种方法的指标对比

图6和图7分别为比较例二中的典型裂纹扩展预测和参数估计结果。其中，图6中左侧图为裂纹扩展轨迹的预测结果，右侧为在给定阈值下的剩余寿命概率分布。从图6可以看出，与比较例一中相似，三种方法中本发明的RLS-KS方法具有相对更好的预测性能。由图7可以看出，AE方法和KS方法中出现了严重的“协调变化”，在迭代初始，参数log C和m分别向其各自的真实值靠近，但在约70000循环左右停止移动，最终收敛于错误的参数组合值，而错误的参数组合却均提供了正确的裂纹扩展预测，如图6(a)和(b)所示。相反，在本发明的RLS-KS方法中，随着迭代的进行参数log C和m不断向其各自的真实值靠近，并最终于15000循环左右收敛于其真实值。表3给出了比较例二中三种方法的指标对比。从HR的结果来看，AE方法和KS方法均能获得正确的剩余寿命预测结果，但却不能正确估计模型参数；而本发明的RLS-KS方法可同时提供正确的剩余寿命预测和模型参数估计。观察其他两种指标可以看出，本发明的RLS-KS方法的RE均小于另外两种方法，参数的RIW稍大但剩余寿命的RIW却明显较小。因此，本发明的RLS-KS方法在裂纹扩展预测和参数估计方法均要好于AE方法和KS方法。

表3比较例二中三种方法的指标对比

由以上分析可以得出，当参数先验与其真实值相距较远的情况下，本发明所提出方法有效解决了“协调变化”的问题，可实现对Paris模型参数log C和m的准确估计。同时，不论是在参数先验与真实值相距较近还是相距较远的情况，本发明所提出的方法均提高了对裂纹扩展剩余寿命的预测性能。

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

Claims (11)

1.一种基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于，包括以下步骤：

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

B、采用递推最小二乘-核平滑方法进行状态模型参数转移；

C、进行裂纹状态转移；

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

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

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

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

2.根据权利要求1所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤D还包括若无新的裂纹监测，则返回步骤C。

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

4.根据权利要求1所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤A所述定义状态模型与观测模型具体为：

x_k＝f(x_k-1，ω_k) (1)

z_k＝g(x_k,ν_k) (2)

其中，k指离散的时刻，x_k为k时刻系统的状态矢量，x_k＝f(x_k-1，ω_k)为状态模型，ω_k为状态转移噪声；z_k为状态的观测矢量，z_k＝g(x_k,ν_k)为观测模型，ν_k为观测噪声。

5.根据权利要求4所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤A具体包括：

A1、定义状态模型：以Paris模型作为疲劳裂纹扩展模型定义状态模型，其参数ln C和m为待估参数，与裂纹长度共同构成扩展的状态矢量：

其中：x_k是k时刻的裂纹长度，θ_k为k时刻的参数矢量，lnC_k与m_k为k时刻的Paris模型参数，据此定义状态模型为：

其中：ΔN表示两个相邻时刻之间的循环周次间隔，F表示形状因子，Δσ表示应力范围，ω_k-1为k-1时刻的状态模型过程噪声，表示ω_k-1服从均值为方差为的正态分布；

A2、定义观测模型为：

其中，z_k为裂纹观测值，ν_k为观测噪声，表示ν_k服从均值为0、方差为的正态分布。

6.根据权利要求5所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤B具体包括：

B1：对状态模型左右两侧分别取期望值，并进行双对数变换：

其中：为k时刻裂纹长度粒子的期望值，和分别为k-1时刻裂纹长度和模型参数粒子的期望值，假设在k时刻的观测误差为零，则有等于z_k，可将上式写为下面的形式：

式中：θ_l,k-1为经过双对数变换后的新的参数矢量，为输入矢量，y_k为输出变量，三者表达式如下：

B2：采用带遗忘因子的最小二乘方法对参数矢量θ_l,k进行估计：

式中：G_k为增益；P_k-1为计算中的过程矩阵，其初始值设置为P₀＝10⁶I，I为二阶单位矩阵；λ为遗忘因子，其值范围通常设置为[0.95,1]，得到θ_l,k，可通过反算得到k时刻经过带遗忘因子的递推最小二乘法估计得到的参数矢量

式中：θ_l,k(1)和θ_l,k(2)分别表示矢量θ_l,k的第一个元素和第二个元素；

B3：对参数粒子进行收缩，计算每个参数粒子的核位置矢量：

式中，i为粒子标号，为k时刻第i个参数粒子的核位置矢量，为k-1时刻第i个参数粒子矢量，α_k为核参数，计算如下：

式中，α₀为核参数的初始值，设置为[0.95,1]；

B4：对参数粒子进行高斯噪声抖动，

首先计算高斯噪声的方差：

式中：为高斯噪声的方差矩阵，表示的方差矩阵；

然后计算转移后新的模型参数粒子：

其中，表示一个由均值为0、方差为的二维正态分布生成的随机数。

7.根据权利要求6所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤C具体包括：

根据步骤A定义的状态模型，将步骤B得到的新的参数粒子集与上一时刻的裂纹长度粒子集带入状态方程，得到下一时刻的裂纹状态粒子集：

其中，表示k时刻第i个粒子的裂纹长度，是k-1时刻第i个粒子的裂纹长度，和表示第i个粒子的模型参数，ΔN表示两个相邻时刻之间的循环周次间隔，F表示形状因子，Δσ表示应力范围。

8.根据权利要求7所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤D具体包括：

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

式中：为k时刻第i个粒子的权值，p(|)表示条件概率；

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

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

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

9.根据权利要求8所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤E具体包括：

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

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

10.根据权利要求9所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤F具体包括：

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

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

11.根据权利要求10所述基于递推最小二乘-核平滑粒子滤波的疲劳裂纹扩展预测方法，其特征在于步骤G具体包括：对于给定的裂纹阈值x_th，对剩余寿命进行预测，剩余寿命的概率分布为：

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