基于高斯权值‑混合粒子滤波的疲劳裂纹扩展预测方法与流程

本发明涉及一种基于高斯权值-混合粒子滤波的疲劳裂纹扩展预测方法，属于故障预测与健康管理技术领域。

背景技术：

随着现代工程系统的复杂程度和综合化程度越来越高，传统的事后维修和定期维护策略已经无法满足其保障维护的要求。近年来，故障预测与健康管理(Prognostics and Health Management，PHM)技术可以通过传感器在线获取系统的实际健康状态，并预测其退化情况，从而制定最佳的运行和维护策略，以实现视情维护。因此PHM技术受到越来越多的关注。工程结构作为系统的核心组成部分，其服役环境复杂并承受各种交变载荷，容易产生疲劳裂纹。疲劳裂纹的存在和扩展将严重削弱结构承载能力，甚至导致灾难性的事故的发生。因此对结构的疲劳裂纹进行在线准确地预测具有重要的理论意义和工程应用价值。

然而疲劳裂纹扩展过程是一个包含各种不确定性的过程，比如材料性质的不确定性，服役环境的不确定性，以及载荷的不确定性等。同时基于结构健康监测方法的在线裂纹监测也受到各种不确定性因素的影响。这些不确定性因素严重影响了疲劳裂纹扩展预测的准确性。近年来，基于贝叶斯理论的概率方法通过结合裂纹扩展模型与实际的结构健康监测数据以消除这些不确定性的影响，受到越来越多的关注。在结构健康监测方法中，基于导波的结构健康监测方法具有小损伤敏感，监测范围广等优点，被认为是非常有前景的方法之一。

在基于贝叶斯理论的方法中，粒子滤波方法由于不需要满足线性和高斯过程假设，被认为非常适用于解决非线性非高斯的疲劳裂纹扩展预测问题。但将粒子滤波方法应用于疲劳裂纹扩展预测时，常常难以准确地定义先验裂纹扩展状态方程，也就是描述当前结构疲劳裂纹扩展的模型，这会加剧粒子滤波算法的粒子退化和多样性匮乏问题，使得裂纹扩展预测结果存在较大的误差。

技术实现要素：

本发明针对先验定义的疲劳裂纹扩展状态方程通常与结构实际的裂纹扩展存在较大偏差的问题，提出了一种基于高斯权值-混合粒子滤波的疲劳裂纹扩展预测方法，结合基于主动导波的裂纹监测方法实现结构疲劳裂纹扩展的准确预测。

本发明为解决其技术问题采用如下技术方案：

一种基于高斯权值-混合粒子滤波的疲劳裂纹扩展预测方法，包括如下步骤：

(1)在结构处于服役条件下，采用基于主动导波的结构健康监测方法对结构的疲劳裂纹进行在线监测，通过获取的导波信号计算损伤因子，并以损伤因子作为裂纹长度的观测值；同时根据结构的实际服役情况定义结构的疲劳裂纹扩展状态空间模型；

(2)在结构处于服役条件下，每当获得一个新的裂纹长度观测值，按照定义的混合系数λ，从均匀分布U(0,1)中随机采样得到随机数，如果该随机数大于λ，从观测概率密度中随机采样得到一个表征裂纹长度的粒子，否则从先验转移概率密度中随机采样得到粒子，重复该过程的得到确定数目粒子组成的粒子集；

(3)对于从先验转移概率密度中采样得到的粒子，通过观测似然概率密度加权；而对于从观测概率密度中采样得到的粒子，根据先验估计定义一个高斯权值分布，通过粒子与先验估计之间的高斯权值作为粒子的权值，然后对所有粒子的权值进行归一化，得到粒子的归一化权值；

(4)通过粒子集和相应的归一化权值表征裂纹长度的后验概率密度，并计算裂纹长度的后验估计；在获得的后验概率密度基础上，将粒子向未来时刻投影得到当前时刻的疲劳裂纹扩展预测结果；

(5)根据每个粒子的归一化权值，使用系统重采样算法进行重采样；

(6)重复上述步骤(2)，(3)，(4)，(5)，(6)实现在线的疲劳裂纹扩展预测。

步骤(2)所述混合系数λ取值为0.5；此外，观测概率密度定义为p(x_k|y_k)，式中y_k为k时刻获得的裂纹长度观测值，x_k为k时刻的裂纹长度；观测概率密度的概率密度形式与表征观测不确定性的概率密度一致，以观测值y_k对应的裂纹长度为均值，k时刻先验粒子集的标准差作为观测概率密度的标准差。

步骤(3)中定义了一个高斯权值分布，对从观测概率密度中采样得到的粒子加权，权值计算方式如下，

式中，是从观测概率密度p(x_k|y_k)中采样得到的第i个粒子，为第i个粒子在k-1时刻的权值；为第i个粒子在k时刻的权值；σ_v为表征观测不确定性的概率密度标准差；为定义的高斯权值分布，其均值为裂纹长度的先验估计标准差σ_p为经验设定值。

本发明的有益效果如下：

本发明提出了一种基于高斯权值-混合粒子滤波的疲劳裂纹扩展预测方法，以观测概率密度和先验转移概率密度的线性组合作为粒子滤波的重要性密度函数。从重要性密度函数中采样得到的粒子集同时考虑了观测信息和先验信息，降低了对准确定义疲劳裂纹扩展状态方程的依赖，同时结合了基于主动导波的裂纹监测方法，能有效地用于结构的疲劳裂纹扩展预测。

附图说明

图1为本发明提出的疲劳裂纹扩展预测方法流程图。

图2为实施例中的结构尺寸图。

图3为实施例中的结构疲劳裂纹扩展图。

图4为实施例中结构健康监测方法获得的导波损伤因子图。

图5为实施例中裂纹长度后验估计对比图。

图6为实施例中预测的失效循环载荷数的相对误差对比图。

具体实施方式

下面结合附图对本发明创造的技术方案进行详细说明。

图1所示，为一种基于高斯权值-混合粒子滤波的疲劳裂纹扩展预测方法的流程图，方法步骤如下，

(1)在结构处于服役条件下，采用基于主动导波的结构健康监测方法对结构的疲劳裂纹进行在线监测。在结构的关键部位粘贴压电传感器，通过压电传感器激励和采集结构中的导波信号。利用结构健康监测系统实现导波信号的在线采集和调理，对获取的导波信号，计算其损伤因子作为裂纹长度的观测值；

(2)根据结构的实际情况定义疲劳裂纹扩展状态空间模型，包括状态方程和观测方程两个部分，其中状态方程由结构的形状及其服役条件确定，观测方程由基于导波的结构健康监测方法确定，如下所示。

状态方程：x_k＝x_k-1+Δx_k-1(x_k-1)·exp(ω)

观测方程：y_k＝g(x_k)+v

式中，k是离散的时间；x_k是k时刻的裂纹长度；x_k-1是k-1时刻的裂纹长度；Δx_k-1(x_k-1)是裂纹扩展增量，为x_k-1的函数，由疲劳裂纹扩展模型定义，比如Paris模型，NASGRO模型；ω是一个随机变量，服从高斯分布表征疲劳裂纹扩展的不确定性，其均值选择为以保证exp(ω)的期望为1，指数项exp(ω)保证裂纹扩展增量非负；y_k为k时刻获得的导波损伤因子；g(·)是观测映射，表征裂纹长度和损伤因子之间的关系；v是一个随机变量，表征观测的不确定性。

(3)定义状态空间模型后，通过结构健康监测方法在线获取结构当前时刻k的裂纹长度观测值y_k，每当获得一个新的裂纹长度观测值，从如下所示的混合建议分布中采样，

式中，k是离散的时间；x_k是k时刻的裂纹长度；x_k-1是k-1时刻的裂纹长度，y_k是k时刻的裂纹长度观测值；i为粒子序号，i＝1,...,N_s，N_s为粒子总数；为重要性密度函数；为先验转移概率密度；λ为混合系数取为λ＝0.5；p(x_k|y_k)为观测概率密度，其形式与表征观测不确定性的概率密度一致。

对于观测概率密度，其均值为观测值y_k对应的裂纹长度g^-1(y_k)；由上一时刻粒子集通过状态方程一步转移得到先验粒子集以的标准差作为观测概率密度的标准差。

从混合建议分布中采样得到粒子集的方式如下：

(a)令i＝1

(b)从均匀分布U(0,1)中随机采样得到u⁽ⁱ⁾，

如果u⁽ⁱ⁾＞λ，从p(x_k|y_k)中随机采样得到

如果u⁽ⁱ⁾≤λ，从中随机采样得到

(c)i＝i+1，并且重复步骤(a)(b)，直到i＝N_s

(4)采样得到粒子集后，分别计算每个粒子的权值。从先验概率分布中采样得到的样本，其权值为，

式中，为第i个粒子在k-1时刻的权值；为第i个粒子在k时刻的权值；为观测似然概率密度。

如果该粒子是从观测概率密度p(x_k|y_k)中采样得到的，权值计算如下，

式中，为第i个粒子在k-1时刻的权值；为第i个粒子在k时刻的权值；为定义的高斯权值分布，其均值为裂纹长度的先验估计标准差为σ_p，该条件概率密度表征粒子与先验估计之间的高斯权值。先验估计由先验粒子集加权求和得到，比例项(σ_p/σ_v)保证两种不同加权方式得到的粒子权值在同一数量级，式中σ_v为表征观测不确定性的概率密度标准差。

计算得到所有粒子的权值后，对所有粒子的权值进行归一化处理，如下式所示

当前时刻裂纹长度的后验概率密度由粒子集以及对应的归一化权值表征如下，

式中，p(x_k|y_1:k)为当前时刻裂纹长度的后验概率密度；δ是狄利克雷函数，其表达式如下，

当前时刻裂纹长度的后验估计由粒子加权求和得到，

(5)在获得的粒子集和归一化权值的基础上，通过裂纹扩展状态空间模型预测每个粒子在未来时刻的疲劳裂纹扩展，得到粒子集里面所有粒子的裂纹扩展预测结果，其表达式如下，

式中为第i个粒子在未来时刻k+d的裂纹长度，为第i个粒子在未来时刻k+d-1的裂纹长度(d＝1,2,3,......)。最后得到预测的裂纹长度概率密度以及裂纹长度预测值如下，

式中，为k时刻的粒子权值；p(x_k+d|y_1:k)为预测的k+d时刻裂纹长度概率密度，粒子为预测的k+d时刻裂纹长度

(6)根据归一化权值对所有的粒子进行重采样，采用系统重采样算法，步骤如下：

(a)计算归一化权值的累积数列C，共N_s个元素。

(b)从均匀分布U(0,1)中随机采样得到u，

(c)令i＝1

计算u⁽ⁱ⁾＝u+(i-1)/N_s

顺序查找数列C中第一个大于u⁽ⁱ⁾的元素下标j，

令重采样得到的粒子并且权值

(d)i＝i+1，并重复步骤(c)，直到i＝N_s

(7)系统重采样算法完成后，将重采样得到的粒子和相应的权值作为当前值，继续通过结构健康监测方法获得新的损伤因子，k＝k+1，算法循环迭代。

本实施例中以仿真的AL2024-T6铝板的中心裂纹扩展为例来具体说明本发明方法的具体实施过程。

如图2所示的中心裂纹铝板，其厚度为3mm，尺寸为b＝200mm,h＝400mm。所承受疲劳载荷为常幅正弦载荷，载荷峰值S_max＝60Mpa，应力比为R＝0。通过下式所述的NASGRO模型仿真结构真实的疲劳裂纹扩展，

式中，k为离散的时刻，x_k为k时刻的裂纹长度，x_k-1为k-1时刻的裂纹长度。ΔK为应力强度因子幅；K_max为应力强度因子最大值；ΔN为载荷循环的步长，C，m，p，q为经验常数；K_C为断裂韧度；ΔK_th为应力强度因子幅阈值；f为裂纹张开系数；应力强度因子式中形状参数b为矩形板宽度，S为矩形板承受的拉伸应力。

取ΔN＝50个载荷循环；f＝0.3958；C＝exp(-26.76)；m＝3.2；p＝0.25；q＝1；初始裂纹长度x₀＝3mm，得到疲劳裂纹扩展轨迹如图3所示。

假设导波损伤因子与实际的裂纹长度存在如下关系，

y＝-5×10^-5(x-x₀)³+2×10^-3(x-x₀)²+1.5×10^-2(x-x₀)+v

式中，x为裂纹长度；y为导波损伤因子；x₀为初始裂纹长度；观测过程受到高斯白噪声v的影响，其服从高斯分布N(0,0.03²)；并且每隔15000个载荷循环进行裂纹长度观测，得到的表征裂纹长度的损伤因子如图4所示。

首先根据具体情况定义该结构的状态空间模型，包括状态方程和观测方程。其中状态方程定义如下，

式中，x_k为k时刻的裂纹长度；x_k-1为k-1时刻的裂纹长度；参数C选取为C＝exp(-26.6)以体现先验定义的裂纹扩展状态方程不准确的情形，其余参数保持不变；表征裂纹扩展不确定性的随机变量ω标准差为σ_ω＝1.3。

观测方程定义根据观测映射定义，如下式所示。

y_k＝-5×10^-5(x_k-x₀)³+2×10^-3(x_k-x₀)²+1.5×10^-2(x_k-x₀)+v

式中，y_k为导波损伤因子；x₀为初始裂纹长度；观测不确定性v服从高斯分布N(0,0.03²)，其标准差σ_v＝0.03。

接下来对高斯权值-混合粒子滤波算法进行初始化，选取粒子数N_s＝500，高斯权值分布的标准差σ_p＝0.4。在初始时刻k＝0所有的粒子的裂纹长度都初始化为x₀。

在结构处于服役过程中，通过基于主动导波的结构健康监测方法对结构的裂纹长度进行在线监测，序贯获得导波损伤因子。在当前时刻，一旦获得新的损伤因子，粒子滤波方法首先计算裂纹长度的后验概率密度，以及裂纹长度的后验估计，然后在该后验概率密度的基础上对未来的疲劳裂纹扩展进行预测。图5所示为结合了所有观测值后的裂纹长度后验估计结果，其中先验模型表征的是裂纹扩展状态方程的裂纹扩展期望值。可以看到，先验模型的误差很大。而本发明的方法在裂纹长度观测点的后验估计均方根误差仅为0.28mm。另一方面，假设以裂纹长度18.69mm为临界裂纹长度，超过该裂纹长度的循环载荷数被认为是失效循环载荷数。本发明提出的方法预测失效循环载荷数的相对误差如图6所示，可以看到，本文提出的方法可以较为准确的预测结构疲劳裂纹扩展。