朱旭程

(海军航空工程学院飞行器工程系 烟台，264001)



旋翼振动监测信号的演化分析与损伤跟踪方法*

朱旭程

(海军航空工程学院飞行器工程系 烟台，264001)

为获得一种对损伤敏感而对正常状态扰动不敏感的旋翼损伤观测信号、解决直升机旋翼在飞行状态下的损伤监测问题，研究了旋翼在重构相空间中的损伤演化特性并提出一种新的损伤跟踪方法。该方法首先应用嵌入技术将旋翼气弹有限元模型产生的振动监测信号重构到维数更高的相空间中，采用非线性Volterra级数建立旋翼基准状态预测模型，以基准状态预测结果与实测相轨迹之件的差异作为状态预测残差，在多个邻域内对其统计平均形成一个损伤观测特征向量；然后,应用奇异谱分解方法从损伤观测特征向量的时间序列中提取出桨叶损伤演化的维度和趋势信息，采用双指数平滑方法建立损伤演化趋势预测模型并估计出桨叶损伤故障剩余寿命。采用桨叶损伤模型和旋翼气弹模型仿真数据验证了所提方法的可行性和有效性。结果表明:该方法能充分利用监测信号的非线性特性在高维相空间中重建系统动力学本质，以不同的时间尺度来观测系统的演化特性；具有损伤模式自动识别能力，可用于难以事先确定系统损伤演化模型或维数信息的场合。

直升机；旋翼；健康监测；故障诊断；相空间

引言

旋翼是直升机飞行的关键动力部件，在桨叶周期变距和气动交变载荷作用下桨叶、桨毂、减摆器和变距拉杆等部件会出现磨损和疲劳失效故障。能够准确监测旋翼的损伤状况，对于提高直升机飞行安全性极具现实意义。文献[1]对旋翼损伤监测问题开展了研究，提出了通过桨叶相对位移和机体振动识别旋翼故障的方案。Ganguli等[2-3]分析了典型损伤对旋翼气弹响应信号的影响，运用神经网络、模糊推理等人工智能技术建立了一些旋翼故障诊断算法。目前，旋翼在飞行条件下的损伤监测技术还存在瓶颈[4]：在高速旋转的旋翼上加装各种数据采集系统尚存在一些技术上的限制；此外，由于桨距时变、气弹耦合和非定常尾流等因素影响，传统损伤跟踪方法难以获得一种对损伤敏感而对正常状态扰动不敏感的损伤观测信号，有效的旋翼损伤跟踪方法依然缺乏。

笔者应用混沌信号嵌入分析理论的相空间重构技术，将旋翼状态监测数据变换到更高维数的状态空间内、以不同的时间尺度来研究旋翼系统的损伤演化特性。采用非线性Volterra级数建立旋翼基准状态预测模型，用基准状态预测结果与实测相轨迹之间的差异作为状态预测残差，发展一种新的旋翼损伤跟踪方法。

1 旋翼在重构相空间的损伤演化特性

1.1 旋翼气弹有限元模型

为研究损伤对桨叶位移、桨毂力等旋翼状态监测信号的影响，需要建立旋翼气弹有限元模型。由Hamilton能量变分原理

(1)

其中：Π为系统总能；T为动能；U为应变能；W为外力虚功。

气动载荷采用叶素理论建立，旋翼诱导速度采用White-Blake入流模型[5]计算。

采用有限元方法离散上述能量积分，进一步采用分部积分法建立旋翼气弹有限元方程

(2)

其中：M，C，K，F分别表示惯性矩阵、阻尼矩阵、刚度矩阵和载荷向量，在定常平飞状态下它们为方位角的周期函数。

旋翼气弹有限元方程可采用时域有限元方法化为代数方程求解，或者将其转换为非线性状态方程采用数值积分方法求解。

旋翼的故障模式种类多样，不同类型的故障对系统参数会有不同影响。采用复合材料桨叶刚度疲劳损伤演化模型说明本研究方法的有效性。

(3)

桨叶损伤变量随使用寿命的扩展如图1所示。

图1 桨叶刚度损伤变量扩展曲线Fig.1 Damage accumulation curve of blade stiffness

目前，直升机健康和使用监测系统(health and usage monitor system,简称HUMS)能够提供的旋翼状态监测信号主要有二类：a. 桨叶的挥舞、摆动和扭转位移，可通过在机身上安装激光多普勒测振仪实时监测[6]；b. 桨毂力和力矩，可通过机身上的应变计或振动传感器结合飞行参数得到。下面以桨尖位移、桨毂力和力矩作为旋翼典型状态监测信号。

将旋翼气弹有限元模型和桨叶损伤演化模型结合，分析桨叶损伤对桨尖位移和桨毂载荷的影响。笔者以BO-105直升机旋翼为参考设置实例旋翼的基本参数：桨叶数为4，旋翼半径R为4.94 m，旋翼转速Ω为40 rad/s，桨叶密度m为6.46 kg/m，桨叶弦长为0.055R，弦向刚度为0.026 8mΩ2R4，挥向刚度为0.010 8mΩ2R4，扭转刚度为0.006 15mΩ2R4，洛克数为5.2，桨叶升力线斜率为6.0，型阻系数为0.002，诱导阻力系数为0.2，俯仰力矩系数为0.00。飞机在前进比0.3，旋翼拉力系数0.004 9的条件下定常平飞。旋翼气弹模型计算程序采用Matlab软件建立(如图2)，桨叶在健康(t/tF=0)和损伤(t/tF=1/3)状态下仿真产生的桨尖位移、桨根力和力矩如图3所示。图中:w为桨尖挥舞位移；v为桨尖摆动位移；Fz为桨毂力；Mx为桨毂力矩；n为样本长度。

图2 旋翼气弹模型计算流程Fig.2 Flowchart of rotor simulation procedure

图3 健康和损伤旋翼桨尖位移和桨毂载荷Fig.3 Health and damaged blade tip deflections and Hub loads

1.2 旋翼状态监测信号的相空间重构

旋翼损伤监测的一个基本问题是找到一种对损伤敏感而对系统状态变化、干扰噪声不敏感的观测信号。然而桨叶损伤与气弹耦合现象会引起气动载荷沿桨叶结构的分布发生变化[7]，随着桨叶损伤程度的扩展，在挥舞、摆动、扭转方向的一些固有频率互相穿越与交叉，使得桨尖位移和桨毂载荷的峰值等特征信号在单调性关系上出现突变，因而难以由这些特征信号辨识出损伤的变化趋势。源自混沌理论的相空间重构技术可以在不同维数的空间内、以不同的时间尺度来研究演化系统的动力学特性，为非线性振动系统的损伤跟踪提供了一条新途径。只要时间尺度和空间维数选择合适，可通过相轨迹的移动变化定量跟踪系统损伤的演化过程[8]。将某监测信号的时序数据{y(i):i=1,2,…,M} 进行时间坐标延迟，构造一状态点列 {x(n):n=1,2,…,N}。其中：x(n) = [y1(n),y2(n), …,ym(n)]T；yk(n) =y(n+(k-1)τ)；m为嵌入维数；τ为延迟时间。该点列描述了系统在m维状态空间中的相轨迹，由它构成的空间称为相空间，这个过程称为相空间重构。根据嵌入理论[9]，如果重构系统与原动力学系统的动态响应可限制在相同维数的吸引子流形上，重构系统与原系统在拓扑关系或微分意义上属于同一等价(微分同胚)，具有相同的动力学性质和几何性质，系统的许多特性如分维、K熵、李指数等均相同。

图4 平均互信息法和假近邻法求重构参数流程Fig.4 Flowchart for time delay and embedding dimension based on average mutual information and FNN percentage

图5 相轨迹在三维空间中的投影Fig.5 3D views of reconstructed phase space trajectories

为直观表明旋翼在相空间中的损伤演化过程，将带桨叶损伤BO-105旋翼的桨尖位移和桨毂载荷数据进行相空间重构。在相空间重构时，延迟时间和嵌入维数是两个主要参数，采用平均互信息法和假近邻法来确定[10]，计算流程如图4所示。求得延迟时间为30，嵌入维数为4，信号相轨迹在三维空间中的投影如图5所示。可见，旋翼系统的稳态相轨迹(吸引子)在重构空间中完全展开，通过相轨迹随桨叶损伤扩展而出现的漂移和变形等几何变化可直观地了解和跟踪旋翼损伤演化过程。

2 在重构相空间跟踪旋翼损伤的方法

笔者所提的旋翼损伤跟踪方法首先在多个邻域内对状态预测残差统计平均形成一个损伤观测特征向量，然后应用奇异谱分解方法从损伤观测特征向量的时间序列中提取出桨叶损伤演化的维度和趋势信息，采用双指数平滑(double exponential smoothing，简称DES)方法[11]建立损伤演化趋势预测模型并估计出桨叶损伤故障剩余寿命。程序主要包括损伤观测信号生成模块和损伤状态优化估计模块(如图6)。损伤观测信号生成模块事先根据健康旋翼在相空间内的运动轨迹建立基准状态预测模型，运行时利用现场采集数据和预测模型产生短时状态预测残差并统计生产损伤观测特征向量。损伤分析和剩余寿命预测模块采用奇异谱分解方法获得旋翼桨叶损伤维数和演化趋势信息，利用损伤演化趋势预测模型估计桨叶剩余寿命。

图6 旋翼损伤监测程序原理图Fig.6 Rotor health monitoring system framework

2.1 采用Volterra预测模型生成损伤观测信号

在旋翼系统状态空间中，用损伤旋翼和健康旋翼轨迹之间的短时偏差作为系统的损伤观测信号。设旋翼在某任意时刻的状态为xn，损伤为φ，经过短时间τ后旋翼状态变为x(xn,τ,φ)。若按健康系统(φ=φR)估计其状态，则应为x(xn,τ,φR)，两状态之间的残差为

S(xn,τ)(φ-φR)+O(ε)

(4)

短时状态残差利用系统的动力学模型去除信号中因运行时间和系统状态的正常转移而造成的变化，保留了损伤造成的变化，是一种更好的损伤观测信号。若能忽略高阶小量的影响，式(4)可近似为

(5)

(6)

其中:k5(x0;φ)=(k1+2k2+2k3+k4)/6;k1=F(x0;φ)；k2=F(x0+0.5τk1;φ)；k3=F(x0+0.5τk2;φ)；k4=F(x0+τk3;φ)。

应用式(6)计算得到短时状态残差后，就可研究灵敏度S(xn,τ)随系统初始状态和损伤大小的变化情况，具体分析高阶小量的影响。下面分析实例旋翼桨尖挥舞位移的短时残差灵敏度随初始状态和桨叶刚度损伤的变化情况。首先，在健康旋翼的状态空间中随机选择16个点作为计算初始状态，在桨叶寿命周期(0

本文提出了(N,m)维修模型,在充分考虑不同部件的维修优先级不同的条件下，对高维修优先级部件，采用修复性维修以及预防性维修；对低维修优先级部件，采用成组维修。不同于传统的依赖于单一团队的集中式维修模式，本文提出了异构多部件系统的竞争式维修模式。以全周期的最低维修费率为目标，通过计算不同部件的失效次序统计量，对参数进行优化，进而得到最优维修策略。通过数值算例验证了维修策略的有效性，并对参数敏感度进行了分析。

图7 短时状态残差与寿命和损伤的关系Fig.7 Curves of short time state error vs. life time and damage variable

(7)

其中：y为传感器输出序列；x为相点；xR为预测器输出的参考点；ex为状态预测残差；ek为损伤观测信号。Pha为相空间重构算子；m，τ分别为嵌入维数和时延参数；Bi为相空间邻域(i=1,2，…,M)；Ni为邻域内相点数目；Vol为Volterra预测器；p为预测器阶数，Volterra预测器的滤波系可根据健康旋翼的相轨迹采用自适应算法确定[12]。

该方法的计算流程如图8所示。根据健康旋翼在相空间内的运动轨迹，采用非线性Volterra级数建立基准状态短时预测模型。监测时将旋翼的状态采集数据运用延迟嵌入方法重构出实测相轨迹，采用Volterra预测模型计算出参考轨迹，将参考轨迹和实测相轨迹之差统计平均生成损伤观测信号。由于Volterra预测模型在不同区域内的精度变化，单个短时状态残差灵敏度具有较大的不确定性。为提高损伤检测的鲁棒性并增加损伤模式的跟踪数目，将参考相空间划分成灵敏度大体相等且互不重叠的几个区域Bi(i=1,2，…,M)，状态预测残差对多种损伤类型的非线性特性可通过其在不同相空间区域内的灵敏度差异来体现，用各个区域内所有残差的平均值来消除灵敏度矩阵的随机误差，将一条记录中各区域内残差均值组合形成一个维数足够高的损伤观测特征向量ek。

图8 损伤观测信号的计算流程Fig.8 Flow chart of damage measurement prediction

应用该方法对实例旋翼在不同损伤状态下的桨叶位移和桨毂载荷进行分析。在桨叶寿命周期内(0≤t/tF<1)，等时间间隔采集60个损伤程度(0≤φ<1)下的仿真记录，将参考相空间按照坐标分量x1大小划分成8个等间隔的邻域，在各邻域内统计每条记录的短时状态残差，计算所得每个残差特征向量ek包含32个列，各残差特征向量随使用寿命的变化曲线如图9所示。可见，各监测信号的残差特征向量可大致反映桨叶损伤曲线的变化趋势，但数据依然有较大起伏，难以判断其中含有几种损伤发展模态。

图9 残差特征向量随时间的变化Fig.9 Curves of residual tracking metric(Scaled)

2.2 桨叶损伤信息的提取

理论上可以利用残差特征向量ek与损伤变量之间的函数关系φk=H(ek)来估计桨叶损伤。然而，旋翼的实际损伤情况比较复杂，损伤产生的原因可能来自疲劳、腐蚀和照射等不确定性因素，损伤变量的实际维数往往不能事先知道；有些损伤类型不但损伤扩展规律不确定，甚至没有事先可直接测量的损伤变量，难以建立损伤与观测之间映射关系。因此，需要直接应用ek序列分析损伤变量的维数和损伤演化的模态信息。奇异谱分析方法[13]能够研究多维时序数据在模态空间中的能量分布情况，可用于分析旋翼的损伤维数和演化趋势信息。

(8)

式(8)是一种Rayleigh商的优化问题[13]，可化为广义特征值问题

(9)

其中：

Δt为序列的时间间隔；D为微分算子矩阵，大小为(n-1)n；λ为特征值；q为特征向量。

在Matlab软件中广义特征值方程(8)可采用奇异谱分解法求解：(Q,λ) =eig(YTY, (DY)TDY)。矩阵Y在模态坐标系中的投影分量为

φk=Yqk

(10)

其中：qk为第k个广义模态向量；φk为第k个模态投影。

矩阵Y中既包含损伤变量缓慢演化信号，又包含旋翼系统的状态变量快速运动的信号。如果残差矩阵中有k个特征值比其它特征值高出一个数量级，则表明系统存在着k个缓慢的损伤演化模式，这些特征值对应的投影可指示损伤演化过程。

采用奇异谱分解算法对BO-105旋翼各监测信号的残差矩阵Y进行分析，得到的32个广义特征值如图10(a)所示，相应的，各残差矩阵在32个模态坐标系中的投影分量如图10(b)所示。可见，残差信号的能量按模态频率重新分配后，有一个特征值远远大于其他特征值，说明残差矩阵中有一种确定的趋势分量。与之对应模态的投影分量较为清晰地反映桨叶损伤扩展的趋势信息，容易验证桨叶损伤与其之间有较为显著的线性关系。将该模态投影分量作为一种名义损伤变量，可对旋翼的损伤状态进行定量跟踪。

图10 广义特征值和模态坐标投影Fig.10 Eigenvalues and modal projections

2.3 桨叶剩余寿命的预测

为估计桨叶剩余寿命需要建立合理的损伤趋势预测模型。桨叶刚度损伤过程是一种非线性随机过程，双指数平滑方法[11]既能灵活适应回归变量的非线性动态特性，又能快速计算多变量回归预测模型的统计分析参数，是一种优选的自适应预测技术。DES方法的基本预测公式为

(11)

当给定失效条件φF后，根据式(10)求得桨叶剩余寿命。采用DES方法对实例旋翼各残差矩阵的模态投影分量进行分析，得到桨叶剩余寿命估计结果如图11所示。在使用时间t≥0.2tF区间内，采用各特征模态投影分量对桨叶剩余寿命跟踪的误差如表1所示。可见，这种方法能较为有效地跟踪桨叶寿命过程。

图11 应用不同观测数据对剩余寿命估计Fig.11 Prediction of remaining time to failure using different measurements

表1 桨叶剩余寿命的跟踪误差

3 结束语

对结构损伤机理和材料失效理论的研究通常都是在微观尺度上进行的。这种研究是必要的，但这类研究不能阐明这些损伤过程在整个系统宏观行为中是如何表现的，也不能提供一种系统级的损伤跟踪手段。源自混沌信号分析方法的相空间重构技术可以从不同的时空角度观察系统的损伤行为特性。笔者利用信号重构技术和Volterra预测模型发展了一种新的损伤观测残差生成方法，对损伤观测残差序列分析可提取损伤变量的维度和趋势信息。所提方法能充分利用被监测信号的非线性特性在高维相空间中重建系统动力学本质，以不同的时间尺度来观测系统的演化特性。该方法具有损伤模式的自动识别能力，可用于损伤演化模型或维数信息难以事先确定的场合。

[1] Azzam H, Andrew M J.The use of math-dynamic model to aid the development of integrated health and usage monitoring[J].Journal of Aerospace Engineering,1992,206(1):71-96.

[2] Ganguli R,Chopra I,Haas D J.Detection of helicopter rotor system simulated faults using neural networks[J].Journal of American Helicopter Society,1997,42(2):161-171.

[3] Pawar P M, Ganguli R.Fuzzy logic based health monitoring and residual life prediction for composite rotor[J].Journal of Aircraft,2007,44(3): 981-995.

[4] Pawar P M, Ganguli R. Helicopter rotor health monitoring-a review[J]. Journal of Aerospace Engineering,2007,221(1): 631-647.

[5] Johnson W. Comprehensive analytical model for rotorcraft aerodynamics and dynamics volumes ii:components theory [M].Palo Alto:Johnson Aeronautics,2007:280-283.

[6] Miles T J, Lucas M.Torsional and bending vibration measurements on rotors using laser technology[J]. Journal of Sound Vibration, 1999,266(3):1441-1467.

[7] Pawar P M,Ganguli R. On the effect of progressive damage on composite helicopter rotor system behavior[J]. Composite Structures, 2007,78:410-423.

[8] Chelidze D, Liu M.Dynamical systems approach to fatigue damage identification [J]. Journal of Sound Vibration,2005, 281(3):887-904.

[9] Hilborn R C.Chaos and nonlinear dynamics[M].Oxford:Oxford University Press,2006:375-389.

[10]Cao Liangyue.Practical method for determining the minimum embedding dimension of a scalar time series[J]. Physica D:Nonlinear Phenomena,1997,110(1-2):43-50.

[11]Guo Zhihao, Shahdi M, Shaya S. Multi-objective OLSR for proactive routing in manet with delay,energy, and link lifetime predictions[J].Applied Mathematical Modeling,2011,35:1413-1426.

[12]张家树.混沌信号的非线性自适应预测技术及其应用研究[D].成都：中国电子科技大学, 2001.

[13]钱征文,程礼,李应红.利用奇异值分解的信号降噪方法[J].振动、测试与诊断,2011，31(4):459-463.

Qian Zhengwen,Cheng Li,Li Yinghong. Noise reduction method based on singular value decomposition[J].Journal of Vibration,Measurement & Diagnosis, 2011,31(4):459-463.(in Chinese)

[14]Marquardt D. An algorithm for least-squares estimation of nonlinear parameters [J]. SIAM Journal on Applied Mathematics, 1963, 11: 431-441．

欢迎订阅《振动、测试与诊断》

《振动、测试与诊断》由工业和信息化部主管，南京航空航天大学和全国高校机械工程测试技术研究会联合主办，是反映振动、动态测试及故障诊断学科领域的科研成果及其应用情况的技术性刊物。主要刊登国内外以振动测试与故障诊断为中心的动态测试理论、方法和手段的研究及应用方面的技术文献，包括实验测试技术、测试仪器的研制、方法和系统组成、信号分析、数据处理、参数识别与故障诊断以及有关装置的设计、使用、控制、标定和校准等，不拘泥于行业和测试项目。

本刊为EICompendex数据库收录期刊和中文核心期刊，双月刊，每逢双月末出版，每本定价20元，全年120元。欢迎订阅和投稿，欢迎在本刊刊登各类广告和科技信息。

编辑部地址：南京市御道街29号 邮政编码：210016 电话：(025)84893332

传真：(025)84893332 E-mail:qchen@nuaa.edu.cn 网址：http:∥zdcs.nuaa.edu.cn

*国家自然科学基金资助项目(51375490)

2013-01-07；

2013-02-24

10.16450/j.cnki.issn.1004-6801.2015.01.020

V275.1; TH113.1

朱旭程,男，1973年11月，博士、副教授。主要研究方向为航空机械维修、综合保障。曾发表《直升机旋翼健康监测系统地面试验台方案设计与开发》(《兵工学报》2008年第29卷第4期)等论文。 E-mail：zhu133@189.cn