1 引言
在航天领域，姿态控制系统的设计、载荷的计算、结构的优化以及元器件的安装等都需要准确的结构模态数据。一般先通过地面模态试验对理论模型进行修正，再利用修正的模型计算各个状态下的模态特性[1]。然而地面模态试验大多局限于时不变结构，无法真实地模拟火箭飞行中的状态，因此开展火箭飞行过程中的时变模态参数辨识对整个火箭的安全性与稳定性具有重要意义。
传统的时变结构模态辨识通常采用“时间冻结”的理论，认为在冻结的时间窗口中结构参数保持不变，再利用成熟的时不变模态辨识方法对各个时间段进行处理最终得到时变模态参数，例如董严等[2]近似地认为火箭结构的模态在每个时间段内保持不变，利用从火箭发射至发动机分离时间段内各测点的测试数据，基于自回归滑动平均模型成功得到了各阶模态参数随时间的变化关系。王鹏辉等[3]采用基于自然激励法的组合方法对火箭模态参数进行识别，将不同时刻的参数进行拟合处理得到火箭的时变模态参数，获得了模态参数随加液氢注量增加的变化规律。但“时间冻结”法忽略了时间变化对结构参数的影响，模态辨识准确度并不高。
近年来随着信号处理技术的发展，将信号分解技术和时频变换相结合的时变模态辨识方法发展迅速。时频变换可以提供频率分辨率，信号分解可以提供时间精度，两者的结合提高了模态辨识的准确程度。例如王佐才等[4]将解析模态分解(analytical modal decomposition,AMD)与小波变换结合用于时变模态辨识，取小波脊线间的平均值作为AMD分解的时变截止频率，提高了瞬时频率的识别精度。王超等[5]利用变分模态分解(variational mode decomposition,VMD)加广义Morse小波的方法，对小车经过主梁时的加速度响应进行分析，成功识别了第一阶模态的时变频率。Wei等[6]将固有啁啾分量分解方法(intrinsic chirp component decomposition,ICCD)[7]与广义参数化时频变换(general parameterized time-frequency transform,GPTFT)[8]结合应用到时变系统自由响应下的模态辨识中，准确度远高于传统的希尔伯特-黄变换(hilbert-huang transform,HHT)。
但上述方法始终是单通道的模态辨识方法，一次仅能处理一个测点的数据，而火箭上的测点往往是成百上千个，逐点处理不仅效率低下，也没有充分利用各测量通道之间的相关性。鉴于目前鲜有学者开展该方面的研究，本文提出了基于多通道固有啁啾分量分解(multivariate ICCD,MICCD)[9]的时变结构模态辨识方法，该方法不仅利用GPTFT来提高瞬时频率(instantaneous frequency,IF)的分辨率，还通过扩展的最小二乘法求解多通道线性方程组来精确提取各时变分量，提高了时变模态参数辨识的准确性。
2 模态参数辨识
一个n自由度的时变结构振动微分方程可以表示为：
(1)
式中： M(t)、C(t)和K(t)分别为时变结构t时刻的质量矩阵、阻尼矩阵以及刚度矩阵；x(t)为结构的位移响应向量； F(t)为外部激励。根据模态叠加法，第v自由度的脉冲位移响应xv(t)可以表示为各阶模态响应的叠加，即：
(2)
其中
式中： K为频率分量总数； φk,v(t)、ak,v(t)分别为第v自由度第k阶位移振型向量值和瞬时振幅(instantaneous amplitude,IA)；qk(t)为第k阶归一化后的模态坐标；ck和φk代表qk(t)的振幅和初始相位； fnk(t)和ζk(t)代表t时刻第k阶的无阻尼固有频率和阻尼比； fdk(t)则为第k阶的有阻尼固有频率，由于结构的阻尼比通常很小，因此可近似认为fnk(t)≈fdk(t)。
值得注意的是，模态叠加法是基于时不变系统提出的。在时变系统领域，并没有关于模态响应的精确的理论推导，大多数研究采取了时变系统响应由多分量瞬时模态响应叠加而成的假设，幸运的是，这些论文的结果证明这种假设对于时变系统是适用的[10]，本文的仿真分析也证明了这一点。
若各通道的各阶的瞬时振幅已知，则可计算出时变结构的第k阶归一化瞬时振型向量φk(t)：
(3)
因此，瞬时频率(IF)和瞬时振幅(IA)的准确提取是模态参数辨识的关键。
3 多通道固有啁啾分量分解
3.1 多通道固有啁啾分量
在语音处理、雷达应用、机械故障诊断等各种应用中，信号通常是非平稳的，可以被建模为调幅和调频信号，也称为啁啾信号。对于多通道的啁啾信号而言，以结构响应为例，多个通道的响应信号都包括了同样的频率分量，故定义多通道固有啁啾分量(multivariate intrinsic chirp component,MICC)如下：
(4)
式中，m(t)即为多通道固有啁啾分量，是相同频率分量成分的集合。结合式(2)，多通道的测量信号X(t)可建模如下：
(5)
式中，αk,v(t)=ak,v(t)cosφk， βk,v(t)=-ak,v(t)sinφk； η(t)则代表分解误差和测量噪声。与单通道ICCD方法类似，可利用冗余傅里叶级数对瞬时振幅ak,v(t)与瞬时频率fk(t)进行拟合：
(6)
其中的待求系数为
其中，ω0=2πF0=2πfs/PN，N为信号样本个数，当P≥2时，上式变为冗余傅里叶级数模型，根据过完备字典下的信号稀疏分解理论[11]，这样做有利于减小拟合误差，提高系数求解和建模的准确性，本文取P=2；L和M分别为IA与IF的冗余傅里叶模型阶数，决定了拟合的复杂程度，模型阶数的确定方法可参考文献[7]。总的来看，式(5)和式(6)说明了瞬时振幅ak,v(t)与瞬时频率fk(t)之间的线性关系，即在瞬时频率fk(t)已知的情况下，通过求解一组线性方程组，可以很容易地得到瞬时振幅ak,v(t)。
3.2 瞬时频率提取
瞬时频率的准确提取依赖于能量高度集中的时频表示，能量越集中，时频脊线提取的准确度也就越高。传统的时频表示方法包括短时傅里叶变换、小波变换和Wigner-Ville分布等。虽然短时傅里叶变换和小波变换在变换基的使用上有所不同，但它们本质上都是使用水平线来近似表达给定信号在时频平面上的瞬时频率，Wigner-Ville分布虽然可以提供高浓度的时频表示，但在分析多分量信号时会存在交叉项干扰，频率分辨性不强，而由Yang等[8]提出的GPTFT在核函数表达式合理的情况下，能够准确表征各种非平稳信号的时频特征。为与式(6)对应，本文采用如下的核函数与变换框架：
(7)
其中
z(t)=s(t)+jH{s(t)}
式中，
和gσ分别为旋转算子、平移算子以及高斯窗函数，采用高斯窗函数的原因是其第一旁瓣衰减程度较海明窗、汉宁窗高，频带边界更加明显，更有利于时频脊线的提取；H{·}代表Hilbert变换；z(t)为多通道能量归一化累加信号s(t)的解析形式：
(8)
其中
式中: *表示归一化后的振幅幅值；θk(t)为第k个分量的相位。当旋转系数
等于IF的傅里叶系数
时，
旋转为平稳信号，平移算子
随后负责将旋转信号的能量在时频平面上移动到实际的瞬时频率处，因此即使面对频率快速变化的啁啾信号，GPTFT也能产生能量集中分布的时频表示。文献[8]提供了一种迭代拟合的方法以提高时频脊线提取的准确率，初始化旋转系数为零，通过式(6)对提取的时频脊线进行拟合来获得光滑的IF曲线和系数
下次迭代时令旋转系数等于上次拟合脊线得到的系数同时增大高斯窗长度来进一步提高频率分辨率。需要注意的是，因为GPTFT的旋转系数仅适用于某阶特定分量，因此每次变换只有一个频率分量的能量非常集中，故每个迭代过程只能提取一个分量成分。对于多频率分量信号，通常采用贪婪算法的思路，优先提取与最强分量相关联的时频脊线，之后将该脊线邻域内的能量分布置零，如此反复，便可获得所有分量的IF曲线及其系数
更多细节可参考文献[7-8]。
3.3 瞬时振幅估计
在单通道ICCD算法中，在已知IF信息的情况下，可建立具有l2范数约束的最小二乘模型求解瞬时振幅(IA)。然而，该模型不能应用于多通道信号，因为多通道信号处理需要同步分解所有通道信号，并不能用向量范数逐个描述，因此需要将最小二乘模型推广到多通道模型。为了更好地解释MICCD方法，将多通道输入信号模型X(t)重写为
X=ΦΘ+η
(9)
式中，X=[x1,x2,…,xV]； Φ=[Φ1,Φ2,…,ΦK]为包含模态信息的核函数，其中
Ck = diag[cos(θk(t0)) … cos(θk(tN-1))]
Sk = diag[sin(θk(t0)) … sin(θk(tN-1))]
θk(ti) = 2π
fk(τ)dτ
Θ为待求IA冗余傅里叶模型系数矩阵：
为使式(9)的分解误差η达到最小，同时避免求矩阵伪逆带来的“病态”问题，引入正则化因子λ：
(10)
式中，
代表Frobenius范数的平方；λ则用来平衡分解准确度与干扰敏感度[7]。上式的解析解为
(11)
式中，
为系数矩阵的计算值；I为大小与ΦTΦ一致的单位阵。瞬时振幅的估计值
则可根据求解出的系数矩阵计算得到：
(12)
对应的频率成分分量
可根据下式重建：
(13)
综上所述，本文算法分为两步：一是利用GPTFT估计多通道固有啁啾分量的瞬时频率信息；二是利用式(11)计算冗余傅里叶级数模型的系数，然后分别用式(12)和式(13)重建瞬时振幅和分量。
4 仿真校验
4.1 仿真设置
在Matlab/Simulink中对如图1所示的三自由度时变结构进行仿真模拟，给定第2个自由度以白噪声激励，使整个系统做随机振动，采样频率fs为100 Hz，采样时间为15 s，求解器选择Runge-Kutta算法[12]。
图1 三自由度弹簧阻尼系统
Fig.1 Three-degree-of-freedom spring damping system
图1中各项时变的物理参数为
m1(t)=0.5,m2(t)=1.5e-0.05t,m3(t)=2.5(kg)
c1(t)=0.2,c2(t)=0.2
c3(t)=0.3+0.03sin(0.5πt)
k2(t)=15 000,k4(t)=12 000
k1(t)=8 000+10 000
(arctan50-arctan(10(τ-8)2))dτ
为衡量瞬时固有频率辨识的精确度，定义如下的频率辨识误差指标：
(14)
其中:
为第k阶频率辨识误差，k值越大代表分量的频率越高； fk(n)和
分别为第n时刻第k分量的理论与辨识IF。
模态保证准则(modal assurance criterion,MAC)能够反映2个振型向量之间的相关程度，因此可以利用MAC值来表征某一时刻瞬时振型辨识的准确程度：
式中： φk为第k分量的理论振型向量；
为辨识出的振型向量；MACk(n)代表第n时刻第k分量的MAC值，MAC值越接近于1代表理论振型与辨识振型相关程度越强，辨识精度越高。
4.2 随机振动分析
利用GPTFT对多通道的累加信号s(t)进行分析来提取瞬时频率(IF)，令IF的傅里叶模型阶数M为15，每组频率分量的拟合迭代次数为3次，高斯窗长度按照128、256、512依次递增，3种频率成分迭代过程的频率辨识误差如表1所示，不难发现随着迭代次数的增加，频率辨识误差越来越小，提取的IF曲线越来越接近理论值，且迭代3次已基本稳定。
表1 各个频率成分迭代拟合的频率辨识误差
Table 1Frequency identification error of iterative fitting of each frequency component Hz
迭代次数第1阶第2阶第3阶10.220.320.3020.210.240.2730.210.230.28
以提取第3阶频率分量为例，GPTFT迭代效果如图2所示，图2中的曲线即为拟合IF曲线，可以随着迭代次数的增加，第3阶的能量分布越来越集中，证明了GPTFT在面对频率变化较快的分量时依然能够提供能量高度集中的时频分布。最终得到的各阶瞬时频率曲线如图3，IF拟合值与理论值基本吻合，说明本算法可以准确地提取瞬时频率。
将上述的IF信息代入式(11)～式(13)得各通道各分量的分解结果及残差如图4所示，图4中mk,v代表了第v通道第k阶啁啾分量；rv代表第v通道信号的分解误差，相比原始信号分解误差能量极小。
图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等[14]基于单通道VMD方法提出的多通道信号分解方法，通过建立约束变分模型，实现多通道信号非递归的自适应分解，将多通道信号在相同的频率尺度分解为相同数量的固有模态函数(intrinsic mode functions,IMFs)之和，保证了多通道信号分解时各阶分量频率的一致性。
在反复进行模态分解试验后，选择最优的分解参数，令分解的IMF数量为3，带宽约束参数为700，迭代收敛容差为10-7。经MVMD分解后，每个自由度的响应均被分解为3个频率一致的IMF分量，每一个IMF分量都反映了时变结构的某一阶固有频率。为获得瞬时固有频率，本文采用Hilbert变换求解每个IMF分量uk,v的瞬时频率和瞬时振幅：
(16)
式中: φk,v、fk,v和ak,v分别为第v通道第k个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自适应方法适用范围更广，分解误差更小，模态辨识精度更高。
[1] 王亮,张妍,蔡毅鹏,周国峰,等.基于滑动时域窗和ERA的模态阻尼比辨识技术研究[J].兵器装备工程学报,2018,39(09):20-25.
Wang L,Zhang Y,Cai Y P,et al.Study on the modal damping ratio identification based on sliding time window method and ERA[J].Journal of Ordnance Equipment Engineering,2018,39(09):20-25.
[2] 董严,付小燕,丁志伟.基于多测点数据的火箭飞行模态参数识别方法[J].固体火箭技术,2018,41(04):520-523.
Dong Y,Fu X Y,Ding Z W.Rocket flight modal identification method based on data of multi-measure points[J].Journal of Solid Rocket Technology,2018,41(04):520-523.
[3] 王鹏辉,黄佳,常洪振,等.自然激励下的运载火箭时变模态参数获取技术研究[J].强度与环境,2021,48(03):1-7.
Wang P H,Huang J,Chang H Z,et al.Time-varying modal acquisition technology of large-scale structure under natural excitation[J].Structure & Environment Engineering,2021,48(03):1-7.
[4] 王佐才,任伟新,邢云斐.基于解析模态分解的时变与弱非线性结构密集模态参数识别[J].振动与冲击,2014,33(19):1-7,16.
Wang Z C,Ren W X,Xing Y F.Analytical modal decomposition-based time-varying and weakly nonlinear structures’ modal parametric identification with closely-spaced modes[J].Journal of Vibration and Shock,2014,33(19):1-7,16.
[5] 王超,毛羚.基于VMD和广义Morse小波的结构瞬时频率识别[J].振动.测试与诊断,2020,40(05):957-962,1026.
Wang C,Mao L.Instantaneous frequency identification of structures based on VMD and generalized Morse wavelet[J].Vibration,Measurement& Diagnosis,2020,40(05):957-962,1026.
[6] Wei S,Chen S,Dong X,et al.Parametric identification of time-varying systems from free vibration using intrinsic chirp component decomposition[J].Acta Mechanica Sinica,2020,36(01):188-205.
[7] Chen S,Peng Z,Yang Y,et al.Intrinsic chirp component decomposition by using Fourier series representation[J].Signal Processing,2017,137:319-327.
[8] Yang Y,Peng Z K,Dong X J,et al.General parameterized time-frequency transform[J].IEEE Transactions on Signal Processing,2014,62(11):2751-2764.
[9] Chen Q,Lang X,Xie L,et al.Multivariate intrinsic chirp mode decomposition[J].Signal Processing,2021,183:108009.
[10] 张杰,史治宇,赵宗爽.基于线调频自适应分解的时变系统瞬时模态参数识别[J].振动与冲击,2020,39(22):103-109,118.
Zhang J,Shi Z Y,Zhao Z S.Instantaneous modal parameter identification of time-varying systems based on adaptive chirplet decomposition[J].Journal of Vibration and Shock,2020,39(22):103-109,118.
[11] Hou T Y,Shi Z Q.Sparse time-frequency representation of nonlinear and nonstationary data[J].Science China Mathematics,2013,56(12):2489-2506.
[12] 黄海阳,王成,李海波,等.基于滑动窗变步长等变自适应源分离的线性慢时变结构工作模态参数识别[J].计算机集成制造系统,2021,27(01):182-191.
Huang H Y,Wang C,Li H B,et al.Moving window variable step-size EASI based operational modal parameter identification for slow linear time-varying structure[J].Computer Integrated Manufacturing Systems,2021,27(01):182-191.
[13] 刘瀚超,王学智,李敏,等.基于多体动力学的发射装置托架拓扑优化设计[J].弹道学报,2021,33(01):11-15，22.
Liu H Z,Wang X A,Li M,et al.Topology optimization of launcher bracket based on multi-body dynamics[J].Journal of Ballistics,2021,33(01):11-15,22.
[14] Rehm An N U,Aftab H.Multivariate variational mode decomposition[J].IEEE Transactions on Signal Processing,2019,67(23):6039-6052.