基于MICCD的时变结构模态参数辨识

王 豪，蓝 鲲，夏国江，耿胜男

(北京宇航系统工程研究所， 北京 100076)

1 引言

在航天领域，姿态控制系统的设计、载荷的计算、结构的优化以及元器件的安装等都需要准确的结构模态数据。一般先通过地面模态试验对理论模型进行修正，再利用修正的模型计算各个状态下的模态特性。然而地面模态试验大多局限于时不变结构，无法真实地模拟火箭飞行中的状态，因此开展火箭飞行过程中的时变模态参数辨识对整个火箭的安全性与稳定性具有重要意义。

传统的时变结构模态辨识通常采用“时间冻结”的理论，认为在冻结的时间窗口中结构参数保持不变，再利用成熟的时不变模态辨识方法对各个时间段进行处理最终得到时变模态参数，例如董严等近似地认为火箭结构的模态在每个时间段内保持不变，利用从火箭发射至发动机分离时间段内各测点的测试数据，基于自回归滑动平均模型成功得到了各阶模态参数随时间的变化关系。王鹏辉等采用基于自然激励法的组合方法对火箭模态参数进行识别，将不同时刻的参数进行拟合处理得到火箭的时变模态参数，获得了模态参数随加液氢注量增加的变化规律。但“时间冻结”法忽略了时间变化对结构参数的影响，模态辨识准确度并不高。

近年来随着信号处理技术的发展，将信号分解技术和时频变换相结合的时变模态辨识方法发展迅速。时频变换可以提供频率分辨率，信号分解可以提供时间精度，两者的结合提高了模态辨识的准确程度。例如王佐才等将解析模态分解(analytical modal decomposition,AMD)与小波变换结合用于时变模态辨识，取小波脊线间的平均值作为AMD分解的时变截止频率，提高了瞬时频率的识别精度。王超等利用变分模态分解(variational mode decomposition,VMD)加广义Morse小波的方法，对小车经过主梁时的加速度响应进行分析，成功识别了第一阶模态的时变频率。Wei等将固有啁啾分量分解方法(intrinsic chirp component decomposition,ICCD)与广义参数化时频变换(general parameterized time-frequency transform,GPTFT)结合应用到时变系统自由响应下的模态辨识中，准确度远高于传统的希尔伯特-黄变换(hilbert-huang transform,HHT)。

但上述方法始终是单通道的模态辨识方法，一次仅能处理一个测点的数据，而火箭上的测点往往是成百上千个，逐点处理不仅效率低下，也没有充分利用各测量通道之间的相关性。鉴于目前鲜有学者开展该方面的研究，本文提出了基于多通道固有啁啾分量分解(multivariate ICCD,MICCD)的时变结构模态辨识方法，该方法不仅利用GPTFT来提高瞬时频率(instantaneous frequency,IF)的分辨率，还通过扩展的最小二乘法求解多通道线性方程组来精确提取各时变分量，提高了时变模态参数辨识的准确性。

2 模态参数辨识

一个自由度的时变结构振动微分方程可以表示为：

(1)

式中：()、()和()分别为时变结构时刻的质量矩阵、阻尼矩阵以及刚度矩阵；()为结构的位移响应向量；()为外部激励。根据模态叠加法，第自由度的脉冲位移响应()可以表示为各阶模态响应的叠加，即：

(2)

其中

式中：为频率分量总数；,()、,()分别为第自由度第阶位移振型向量值和瞬时振幅(instantaneous amplitude,IA)；()为第阶归一化后的模态坐标；和代表()的振幅和初始相位；()和()代表时刻第阶的无阻尼固有频率和阻尼比；()则为第阶的有阻尼固有频率，由于结构的阻尼比通常很小，因此可近似认为()≈()。

值得注意的是，模态叠加法是基于时不变系统提出的。在时变系统领域，并没有关于模态响应的精确的理论推导，大多数研究采取了时变系统响应由多分量瞬时模态响应叠加而成的假设，幸运的是，这些论文的结果证明这种假设对于时变系统是适用的，本文的仿真分析也证明了这一点。

若各通道的各阶的瞬时振幅已知，则可计算出时变结构的第阶归一化瞬时振型向量()：

(3)

因此，瞬时频率(IF)和瞬时振幅(IA)的准确提取是模态参数辨识的关键。

3 多通道固有啁啾分量分解

3.1 多通道固有啁啾分量

在语音处理、雷达应用、机械故障诊断等各种应用中，信号通常是非平稳的，可以被建模为调幅和调频信号，也称为啁啾信号。对于多通道的啁啾信号而言，以结构响应为例，多个通道的响应信号都包括了同样的频率分量，故定义多通道固有啁啾分量(multivariate intrinsic chirp component,MICC)如下：

(4)

式中，()即为多通道固有啁啾分量，是相同频率分量成分的集合。结合式(2)，多通道的测量信号()可建模如下：

(5)

式中，,()=,()cos，,()=-,()sin；()则代表分解误差和测量噪声。与单通道ICCD方法类似，可利用冗余傅里叶级数对瞬时振幅,()与瞬时频率()进行拟合：

(6)

其中的待求系数为

其中，=2π=2π，为信号样本个数，当≥2时，上式变为冗余傅里叶级数模型，根据过完备字典下的信号稀疏分解理论，这样做有利于减小拟合误差，提高系数求解和建模的准确性，本文取=2；和分别为IA与IF的冗余傅里叶模型阶数，决定了拟合的复杂程度，模型阶数的确定方法可参考文献[7]。总的来看，式(5)和式(6)说明了瞬时振幅,()与瞬时频率()之间的线性关系，即在瞬时频率()已知的情况下，通过求解一组线性方程组，可以很容易地得到瞬时振幅,()。

3.2 瞬时频率提取

瞬时频率的准确提取依赖于能量高度集中的时频表示，能量越集中，时频脊线提取的准确度也就越高。传统的时频表示方法包括短时傅里叶变换、小波变换和Wigner-Ville分布等。虽然短时傅里叶变换和小波变换在变换基的使用上有所不同，但它们本质上都是使用水平线来近似表达给定信号在时频平面上的瞬时频率，Wigner-Ville分布虽然可以提供高浓度的时频表示，但在分析多分量信号时会存在交叉项干扰，频率分辨性不强，而由Yang等提出的GPTFT在核函数表达式合理的情况下，能够准确表征各种非平稳信号的时频特征。为与式(6)对应，本文采用如下的核函数与变换框架：

(7)

其中

()=()+{()}

(8)

其中

3.3 瞬时振幅估计

在单通道ICCD算法中，在已知IF信息的情况下，可建立具有范数约束的最小二乘模型求解瞬时振幅(IA)。然而，该模型不能应用于多通道信号，因为多通道信号处理需要同步分解所有通道信号，并不能用向量范数逐个描述，因此需要将最小二乘模型推广到多通道模型。为了更好地解释MICCD方法，将多通道输入信号模型()重写为

=+

(9)

式中，=[,,…,]；=[,,…,]为包含模态信息的核函数，其中

= diag[cos(()) … cos((-1))]

= diag[sin(()) … sin((-1))]

为待求IA冗余傅里叶模型系数矩阵：

为使式(9)的分解误差达到最小，同时避免求矩阵伪逆带来的“病态”问题，引入正则化因子：

(10)

(11)

(12)

(13)

综上所述，本文算法分为两步：一是利用GPTFT估计多通道固有啁啾分量的瞬时频率信息；二是利用式(11)计算冗余傅里叶级数模型的系数，然后分别用式(12)和式(13)重建瞬时振幅和分量。

4 仿真校验

4.1 仿真设置

在Matlab/Simulink中对如图1所示的三自由度时变结构进行仿真模拟，给定第2个自由度以白噪声激励，使整个系统做随机振动，采样频率为100 Hz，采样时间为15 s，求解器选择Runge-Kutta算法。

图1 三自由度弹簧阻尼系统Fig.1 Three-degree-of-freedom spring damping system

图1中各项时变的物理参数为

()=05,()=15e-005,()=25(kg)

()=02,()=02

()=03+003sin(05π)

()=15 000,()=12 000

为衡量瞬时固有频率辨识的精确度，定义如下的频率辨识误差指标：

(14)

模态保证准则(modal assurance criterion,MAC)能够反映2个振型向量之间的相关程度，因此可以利用MAC值来表征某一时刻瞬时振型辨识的准确程度：

4.2 随机振动分析

表1 各个频率成分迭代拟合的频率辨识误差Table 1Frequency identification error of iterative fitting of each frequency component Hz

以提取第3阶频率分量为例，GPTFT迭代效果如图2所示，图2中的曲线即为拟合IF曲线，可以随着迭代次数的增加，第3阶的能量分布越来越集中，证明了GPTFT在面对频率变化较快的分量时依然能够提供能量高度集中的时频分布。最终得到的各阶瞬时频率曲线如图3，IF拟合值与理论值基本吻合，说明本算法可以准确地提取瞬时频率。

将上述的IF信息代入式(11)～式(13)得各通道各分量的分解结果及残差如图4所示，图4中,代表了第通道第阶啁啾分量；代表第通道信号的分解误差，相比原始信号分解误差能量极小。

图2 第3阶频率分量的拟合迭代过程示意图Fig.2 Fitting iterative process of the 3rd frequency component

图3 基于MICCD的瞬时频率曲线Fig.3 Instantaneous frequency identification result based on MICCD

图4 MICCD分解结果及线差图Fig.4 MICCD decomposition results

图4证明了MICCD方法不仅同时处理多个通道的信号，还可以实现同一种频率分量的对齐，且每个分量的重构精度极高。将分解得到的瞬时振幅(IA)代入式(3)可以得到瞬时振型，以第2阶振型为例，辨识振型与理论振型随时间变化趋势如图5所示，可以在误差允许的范围内，辨识振型的变化趋势与理论振型保持一致，说明本文方法可以有效地提取结构瞬时振型。各阶振型的MAC曲线如图6所示，任意时刻的MAC值都十分接近1，说明瞬时振型的辨识精度很高，与理论振型有很高的相关性。

图5 第2阶瞬时振型辨识结果图Fig.5 2nd instantaneous mode shape identification result

图6 基于MICCD的瞬时振型辨识结果曲线Fig.6 Instantaneous mode shape identification result based on MICCD

4.3 对比分析

为了展示本文方法的优越性，又引入MVMD方法对时变结构模态参数进行辨识。MVMD是在2019年由Rehman等基于单通道VMD方法提出的多通道信号分解方法，通过建立约束变分模型，实现多通道信号非递归的自适应分解，将多通道信号在相同的频率尺度分解为相同数量的固有模态函数(intrinsic mode functions,IMFs)之和，保证了多通道信号分解时各阶分量频率的一致性。

在反复进行模态分解试验后，选择最优的分解参数，令分解的IMF数量为3，带宽约束参数为700，迭代收敛容差为10。经MVMD分解后，每个自由度的响应均被分解为3个频率一致的IMF分量，每一个IMF分量都反映了时变结构的某一阶固有频率。为获得瞬时固有频率，本文采用Hilbert变换求解每个IMF分量,的瞬时频率和瞬时振幅：

(16)

式中:,、,和,分别为第通道第个IMF的瞬时相位、瞬时频率与瞬时振幅。图7展示了基于MVMD方法得到的瞬时频率曲线，由于Hilbert方法直接求得的IF曲线呈现反复振荡的特点，并不符合真实结构的频率变化规律，利用式(6)对其拟合获得拟合IF曲线，结果表明前两阶的瞬时频率比较准确，但第3阶出现了明显的模态混叠现象，瞬时频率曲线后半段跳跃到第2阶。对比图3发现基于MICCD方法的瞬时频率辨识精度更高，瞬时振型亦是如此(对比图6与图8)。

图7 基于MVMD的瞬时频率辨识结果曲线Fig.7 Instantaneous frequency identification result based on MVMD

图8 基于MVMD的瞬时振型辨识结果曲线Fig.8 Instantaneous mode shape identification result based on MVMD

5 结论

1) 引入广义参数化时频变换(GPTFT)，提高了时频面内的能量集中程度，再迭代优化提高了瞬时频率辨识的准确率；

2) 扩展最小二乘法求解多元线性方程组，充分利用了各通道间的关系，使得啁啾分量重构更加精确，获得了精确的振幅包络和瞬时振型；

3) MICCD方法需要对各阶的瞬时频率进行准确提取作为输入信息，比MVMD自适应方法适用范围更广，分解误差更小，模态辨识精度更高。