基于辅助粒子滤波的疲劳裂纹扩展预测研究

杨伟博， 袁慎芳， 邱 雷， 陈 健(南京航空航天大学 机械结构力学及控制国家重点实验室, 南京 210016)

疲劳裂纹是当前飞行器的主要损伤形式之一，实现疲劳裂纹的扩展预测，对飞行器延寿，提升飞行安全、降低维护成本具有重要意义。然而结构疲劳裂纹扩展是一个随机过程，具有非线性非高斯的复杂动态特性，它包含各种不确定因素，如材料参数的变化、环境参数的变化、测试数据的分散性、建模的局限性等，如何克服此类不确定性的因素的影响，实现疲劳裂纹这类非线性非高斯动态系统的预测，一直是学者们研究的重点。

目前，可用于实现裂纹扩展的寿命预测方法的主要有失效物理(Physic of Failure, PoF)方法、数据驱动(Data-Driven)方法和融合(Hybrid)算法[1-3]。其中，融合算法意在融合前两种算法的优点，使寿命预测算法具备更低的计算消耗和复杂度、需要更少的实验与训练数据，因此具备更好的发展前景。而作为融合算法的一种，基于粒子滤波(Particle Filter, PF)[4]的疲劳裂纹寿命预测方法结合了Paris规则的裂纹扩展规律，有限元建模和粒子滤波算法，与其他两种方法相比，其建模相对较容易，并不需要大量的实验与训练数据，且与Lamb波主动监测方法相结合，可实现在线预测，因此更适合于工程应用。

近年来，Orchard等[5-6]利用粒子滤波算法，建立了一种裂纹扩展在线预测系统，对UH-60行星齿轮进行了寿命预测，其预测结果相比传统的卡尔曼滤波等方法具有更高的预测精度，且能实时修正监测过程中的预测误差，但观测方程采用应变作为监测手段，观测方程的精度还有进一步的提高空间；Cadini等[7-8]结合文献中的仿真数据，利用粒子滤波算法实现了疲劳裂纹扩展的预测；袁慎芳等[9]以Paris规则作为状态方程，以Lamb波主动健康监测的在线监测技术建立观测方程，利用标准粒子滤波算法实现了中心裂纹铝板上的裂纹扩展在线预测，有效地消除了裂纹扩展参数及观测手段等多种不确定性的影响，获得了较高的预测精度。

从上述粒子滤波的裂纹损伤预测研究可以看出，目前相关领域的研究都是基于最简单的标准粒子滤波算法利用仿真数据或在简单结构中开展的。然而由于采用了次优估计[10]，其预测结果并未考虑当前观测更新的影响，因而依赖于建模的精确性。当模型误差较大时，将有大量粒子处于低相似区中而在重采样时被浪费，从而加剧了粒子多样性匮乏现象，即原有粒子群中很多粒子由于预测偏差大而权值太小没有“后代”，而少数几个权值较高的粒子则有很多相同的“后代”，重采样以后的粒子群由大量重复的粒子构成，使粒子群失去了多样性，粒子多样性匮乏现象的加剧，将导致预测精度下降。

解决上述问题的方法之一是增加粒子个数，然而这将影响计算效率。另一种方法是结合当前观测更新选择重要性函数。本文提出基于辅助粒子滤波(Auxiliary Particle Filter, APF)[11-12]的裂纹扩展寿命预测方法，通过Pairs规则建立状态方程，以基于压电主动Lamb波在线监测方法作为观测方程，构建描述孔边裂纹扩展的状态空间模型。最后，基于孔边裂纹的APF与PF寿命预测结果表明，APF算法充分结合裂纹在线观测值，使重采样后的粒子向似然函数的高值区移动，提高了先验估计与似然函数之间重叠性，能有效缓解多样性匮乏现象，克服PF算法过于依赖模型的缺点，具备更高的预测精度。

1 粒子滤波理论

结构疲劳裂纹扩展是一个动态变化的过程，它存在诸多的不确定性，其数学描述可用状态空间模型来表征。状态空间模型包含状态方程和观测方程

(1)

式中:xt表示t时刻裂纹状态值；f(·)为系统的状态方程，表征t时刻裂纹向t+1时刻扩展的规律；zt+1表示t+1时刻的观测值，此处指通过健康监测手段所获取的损伤特征参数；g(·)为系统的观测方程；ωt为t时刻的过程噪声，用于表征结构上的各种不确定影响因素；υt为t时刻的观测噪声，用于表征健康监测手段的不确定性监测误差。

粒子滤波具备平滑、滤波和预测的功能，在裂纹扩展预测的实现，是根据0到t时刻的裂纹状态值和1到t时刻的裂纹观测值，结合状态空间模型和t+1时刻的最新观测值，去预测t+1时刻及其以后的裂纹状态值，最后再根据裂纹参数计算结构是否可以使用，确定结构的寿命。上述过程实质是通过Bayes估计来实现裂纹扩展过程的无偏估计

(2)

式中状态值x的上标“^”表征后验估计，右上标“-”表征先验估计，p(xt+1|z1:t+1)表征观测值为z1,z2,…zt+1时，出现状态值xt+1的概率密度函数，由于它是通过t+1时刻的观测值zt+1修正的，因此常称为后验概率密度函数。p(x0:t+1|z1:t+1)为p(xt+1|z1:t+1)的边缘概率密度。

假设裂纹扩展满足一阶Markov过程，则后验概率密度函数p(xt+1|z1:t+1)可通过Bayes估计迭代计算获得

(3)

然而式(3)中在迭代过程中包含多重积分，一般情况下并无解析解，因此采用Monte Carlo仿真方法，从后验概率p(x0:t+1|z1:t+1)中抽取N个独立同分布的样本，当N～+∞时，样本的均值可逼近于数学期望，即

(4)

但在实际问题中，往往很难从p(x0:t+1|z1:t+1)中进行抽样，因此引入了一个容易采样的重要性函数q(x0:t+1|z1:t+1)，该参考分布已知且与p(x0:t+1|z1:t+1)同分布，则最终可将裂纹扩展的无偏估计式(2)化简为

(5)

式中：

(6)

以上过程即为序贯重要性采样(Sequential Importance Sampling, SIS)，最终通过裂纹先验估计与其权值的加权，可实现裂纹扩展过程的无偏估计。然而SIS算法存在严重的退化现象，因其最优解q(xt+1|xt,z1:t+1)=p(xt+1|xt,z1:t+1)在实际问题中常无法得到且对应积分不易求解。为此，Gordon提出使用次优解q(xt+1|xt,z1:t+1)=p(xt+1|xt)与重采样相结合以解决退化问题，即序贯重要性采样重采样 (Sequential Importance Resampling, SIR)，又称粒子滤波算法。

2 辅助粒子滤波

SIR算法解决了退化问题且抽样简单，易于实现，然而由于重要性函数使用了次优解p(xt+1|xt)，对比最优解p(xt+1|xt,z1:t+1)，次优解忽略了观测值z1:t+1的影响，更依赖于状态方程xt+1=f(xt,ωt)的建模精度，因此SIR(PF)算法存在过于依赖模型的特点。当重要性函数与真实分布非常接近时，SIR算法能达到很好的预测效果，但如果二者偏差较大，则预测效果将明显下降，甚至导致预测失败。

为克服上述缺陷，Pitty等提出了辅助粒子滤波算法，通过引入辅助变量l，将其与最新观测值的联合概率密度函数作为重要性函数，预先对原有粒子群进行抽样处理，使抽样后的粒子群向真实分布的高似然区移动，以此来克服算法过于依赖模型的缺点，提高预测精度。

假设t时刻p(xt|z1:t)已知，则t+1时刻的“经验”后验概率分布为

(7)

定义结合了辅助变量的联合概率密度函数

(8)

(9)

(10)

(11)

对式(9)进行边缘概率积分可得

(12)

辅助粒子滤波算法执行步骤如下：

3 疲劳裂纹扩展的状态空间模型

3.1 状态方程的建立

金属结构裂纹的扩展存在三个阶段：萌生阶段、稳定扩展阶段和快速扩展阶段。其中萌生阶段和快速扩展阶段的持续时间短，因此工程上进行结构寿命预测时，常采用稳定扩展阶段，此阶段又称为Paris规则区，裂纹扩展速率的Paris公式如式(13)所示

(13)

式中：a为表征裂纹长度；Nf为疲劳载荷周期数，C与m为材料参数，ΔK为应力强度因子幅，可通过有限元离线计算获得，此处采用奇异单元法计算获得。

假设两次扩展之间的时间间隔ΔNf足够小，则根据Paris规则可建立如式(14)所示的状态方程，其中为使过程噪声包含更大的预测空间，状态噪声采用加性的零均值高斯白噪声[13]，而非相乘的对数正态分布噪声。

at+1=at+Ct(ΔK)mΔNf+ω(t)

(14)

式中:at+1为t+1时刻的裂纹长度；Ct为t时刻的材料及环境相关的参数，针对特定的试件，C0可通过实验或查找数据库获得，此处采用实验手段，并结合曲先强验证所总结材料参数分布假设[14]，令lgC0～N(-14.307 5, 0.174 0)，m取值为4.038 3，用于初始化状态方程的不确定性；ω(t)～N(0,0.5)为高斯白噪声，用于表征状态方程预测过程中的不确定性。

式(14)存在两个不确定分布Ct与ω(t)，它们分别分于初始化表征Paris规则的不确定性以及预测过程中的不确定性，二者的结合充分考虑了疲劳裂纹扩展的随机因素，提升了状态方程的鲁棒性。

3.2 观测方程的建立

观测方程在裂纹扩展预测中的意义是将裂纹扩展尺寸与观测手段的某种特征关联起来，定量地描述裂纹扩展与观测手段之间的规律。Lamb波主动健康监测技术具有对小损伤敏感的特性，因此目前已被广泛应用于裂纹扩展的监测中[15-17]。该方法通过监测裂纹扩展时结构中传播的Lamb波的变化，量化分析裂纹尺寸与Lamb波响应信号间的规律，如图1所示为Lamb波响应信号随裂纹尺寸的变化。

图1 响应信号随裂纹尺寸的变化

由图1可知，金属疲劳裂纹扩展对Lamb波信号传播的影响主要表现为相位及幅值的变化，对应后续所开展的孔边裂纹验证实验，主要为裂纹尖端的扩展所导致的Lamb波直达波相位延迟及幅值衰减现象。其中相位延迟现象在本实验中占主导地位，因此根据无损检测中飞行时间(Time of Fight)的概念，定义如式(15)所示的时间延迟损伤因子，用于表征结构损伤程度与相位延迟之间的关系。

(15)

式中:T为裂纹扩展过程中直达波的飞行时间;T0为健康状态下直达波的飞行时间。实验准备初期，先通过几组实验获得DI与裂纹观测值，再通过多项式拟合，建立损伤因子与裂纹扩展之间的观测方程

DIt+1=g(at+1)+υ(t)

(16)

联立状态方程与观测方程，描述疲劳裂纹扩展的状态空间模型如式(17)所示

(17)

4 疲劳裂纹扩展预测的辅助粒子滤波实现

基于式(17)所建立的状态空间模型，为解决PF过于依赖模型的缺点，缓解多样性匮乏现象，将APF算法应用于裂纹扩展预测中。裂纹扩展的辅助粒子滤波实现流程如下，流程图如图2所示。

步骤1 初始化。由Lamb波在线监测获得裂纹的萌生尺寸初始化裂纹状态值a0及其对应循环周期数t0；由n组重复性疲劳实验通过割线法初始化材料参数C0与m的分布；初始化过程噪声ω(t)和观测噪声υ(t)分布；将a0和C0代入状态方程中，初始化裂纹预测粒子群{a0}N，权值{w0}N=1/N；

(18)

(19)

注：t与t+1的时间间隔为ΔNf。

图2 APF预测流程图

5 实验验证

选用孔边裂纹试验件，结合状态空间模型，对辅助粒子滤波改进算法的正确性进行验证。试件如图3所示，尺寸为240 mm×100 mm×2 mm，材料为LY12铝合金，共计6件，编号为T1-T6，传感器对称分布于裂纹两侧，间距80 mm，孔边预制3 mm裂纹。实验加载设备采用MTS810，试件T1先进行静载破坏实验，确定所能承受的最大拉力为45 kN，在此基础上，试件T2-T5采用正弦载荷等幅谱进行疲劳实验，载荷峰值为15 kN，应力比为0.1，载荷频率为10 Hz。疲劳实验过程中，裂纹长度每增加2 mm，即保持载荷至5 kN，采集Lamb波响应信号，并记录循环周期数，监测系统采用南京航空航天大学结构健康监测与预测研究中心自主研发的PHM扫查系统[18-19]，激励频率为290 kHz，采样率为10 MHz。

DI=-6.234×10-5a3+4.854×10-3a2-

4.069×10-2a+9.231×10-2+υ(t)

(20)

作为验证的T6试件，其裂纹扩展过程如表1所示，图4所示为T2-T5的观测方程拟合曲线与T6实验损伤因子的对比，可见T6的实验值均落于拟合曲线的95%的置信区间内，其最大误差为0.144 3，验证了观测手段与观测拟合方程的正确性。

表1 疲劳裂纹扩展过程

图4 T6实验损伤因子与观测方程拟合曲线的对比

选取粒子数N为500，载荷循环迭代间隔ΔNf为40，分别利用PF与APF算法，在实时获取T6试件的损伤因子更新时，实现T6试件的实时监测与疲劳寿命预测。APF与PF的预测结果如图5所示，以预测值与真实值的均方根误差(Root Mean Square Error, RMSE)作为预测精度对比指标，RMSE越小，则精度越高，从图5中可以看出，APF的RMSE为0.438 5，低于PF的0.712 9，因此结合了最新Lamb波DI观测值的APF具备比PF更优秀的预测效果。

图5 PF与APF的预测结果

剃除重采样后的重复粒子(材料参数Ct)，对比APF与PF算法的粒子多样性匮乏程度，以重采样后具有不同数值的材料参数Ct粒子个数作为多样性指标，则APF与PF在迭代过程中粒子多样性变化如图6所示。

从图 6可以看出，相比PF算法，APF算法的粒子多样性匮乏现象得到了缓解，预测精度得到提升。此外，由于进行了二次重采样，运算时间由0.276 1 s增至0.299 3 s(计算采用MATLAB R2013a的64位4G内存Win7系统的笔记本电脑)，虽然增加了部分计算量，但预测精度获得了提升，且这样的计算消耗，对在线监测的实现影响不大。

图6 APF与PF多样性的对比

将PF算法的中粒子数N由500增加至2 500，则APF与PF算法在不同粒子数下的预测结果如图7所示。从图7中可以看出，随着粒子数的增加PF算法的预测精度增加，RMSE由0.712 9降至0.415 2，相比于粒子数为500的APF算法，其RMSE=0.438 5，PF算法的预测精度稍显优势。然而PF算法由于粒子数的增加，运算时间由0.276 1 s增至0.804 6 s。且若将此部分运算消耗提供给APF算法，APF将获得更高的预测精度，因此APF算法的寿命预测性能要优于PF算法。

图7 APF与PF不同粒子数的预测结果

6 结论

APF算法结合了Lamb波主动健康监测技术获得的最新裂纹观测值，克服了传统PF算法过于依赖模型的特点，在孔边裂纹结构的寿命预测中，有效地缓解了多样性匮乏现象，虽然算法计算时间存在小幅度的增加，但并不影响在线预测的实现，预测精度相比PF得到了整体提升，其RMSE由0.712 9降低至0.438 5，更适用于结合了Lamb波在线监测的裂纹扩展实时寿命预测。

[1] ZHANG H, KANG R, PECHT M. A hybrid prognostics and health management approach for condition-based maintenance[C]∥2009 IEEE International Conference on Industrial Engineering and Engineering Management. Hongkong, China: IEEE, 2009.

[2] PECHT M, GU J. Physics-of-failure-based prognostics for electronic products[J]. Transactions of the Institute of Measurement and Control, 2009, 31(3/4): 309-322.

[3] SHEPPARD J W, KAUFMAN M A, WILMER T J. IEEE standards for prognostics and health management[J]. Aerospace and Electronic Systems Magazine, 2009, 24(9): 34-41.

[4] GORDON N J, SALMOND D J, SMITH A F M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation[J].IEE Proceedings F (Radar and Signal Processing). IET Digital Library, 1993, 140(2): 107-113.

[5] ORCHARD M, VACHTSEVANOS G. A particle filtering approach for on-line failure prognosis in a planetary carrier plate[J]. International Journal of Fuzzy Logic and Intelligent Systems, 2007, 7(4): 221-227.

[6] ORCHARD M E, VACHTSEVANOS G J. A particle-filtering approach for on-line fault diagnosis and failure prognosis[J]. Transactions of the Institute of Measurement and Control, 2009, 31(3-4): 221-246.

[7] CADINI F, ZIO E, AVRAM D. Monte Carlo-based filtering for fatigue crack growth estimation[J]. Probabilistic Engineering Mechanics, 2009, 24(3): 367-373.

[8] ZIO E, PALONI G. Particle filtering prognostic estimation of the remaining useful life of nonlinear components[J]. Reliability Engineering & System Safety, 2011, 96(3): 403-409.

[9] 袁慎芳, 张华, 邱雷, 等．基于粒子滤波算法的疲劳裂纹扩展预测方法[J]．航空学报, 2013, 34(12): 2740-2747.

YUAN Shenfang, ZHANG Hua, QIU Lei, et al. A fatigue crack growth prediction method based on particle filter[J]. Acta Aeronautica et Astronautic Sinica, 2013, 34(12): 2740-2747.

[10] DOUCET A, GODSILL S, ANDRIEU C. On sequential Monte Carlo sampling methods for Bayesian filtering[J]. Statistics and Computing, 2000, 10(3): 197-208.

[11] PITT M K, SHEPHARD N. Filtering via simulation: auxiliary particle filters[J]. Journal of the American Statistical Association, 1999, 94(446): 590-599.

[12] PITT M K, SHEPHARD N. Auxiliary variable based particle filters[M]∥Sequential Monte Carlo Methods in Practice. New York: Springer, 2001: 273-293.

[13] SUN J, ZUO H, WANG W, et al. Prognostics uncertainty reduction by fusing on-line monitoring data based on a state-space-based degradation model[J]. Mechanical Systems and Signal Processing, 2014, 45(2): 396-407.

[14] 曲先强, 马永亮, 崔洪斌, 等．Paris 公式中材料参数的统计特性分析[C]//2008 全国MTS断裂测试研讨会论文集, 2008: 26-31.

[15] 曹俊, 袁慎芳, 蔡建, 等．疲劳裂纹扩展的实时健康监测[J]．压电与声光, 2008, 30(6): 776-778.

CAO Jun, YUAN Shenfang, CAI Jian, et al. Real-time structural health monitoring for fatigue crack growth[J]. Piezoelectrics and Acoustooptics, 2008, 30(6): 776-778.

[16] 袁慎芳, 邱雷, 王强, 等．压电-光纤综合结构健康监测系统的研究及验证[J]．航空学报, 2009, 30(2): 348-356.

YUAN Shenfang, QIU Lei, WANG Qiang, et al. Application research of a hybrid piezoelectric-optic fiber integrated structural health monitoring system[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(2): 348-356.

[17] 陆希,孟光,李富才．基于Lamb波的薄壁槽状结构损伤检测研究[J]．振动与冲击, 2012, 31(12):63-67.

LU Xi, MENG Guang, LI Fucai, et al. Lamb wave-based damage detection for a channel-like thin-wall structure[J]. Journal of Vibration and Shock, 2012, 31(12):63-67.

[18] 邱雷, 袁慎芳．集成压电健康监测扫查系统的研制及其应用[J]．压电与声光, 2008, 30(1): 39-41.

QIU Lei, YUAN Shenfang. Research and application of integrated health monitoring scanning system based on PZT sensor array[J]. Piezoelectrics & Acoustooptics, 2008,30(1): 39-41.

[19] QIU L, YUAN S F, WANG Q, et al. Design and experiment of PZT network-based structural health monitoring scanning system[J]. Chinese Journal of Aeronautics, 2009, 22(5): 505-512.

[20] DONG M, WANG N. Adaptive network-based fuzzy inference system with leave-one-out cross-validation approach for prediction of surface roughness[J]. Applied Mathematical Modelling, 2011, 35(3): 1024-1035.