CMIF方法中模态参数不确定性的计算

孙鑫晖，郝木明，李振涛，张令弥，王 彤

（1．中国石油大学（华东）化学工程学院，山东 青岛 266580；2．南京航空航天大学振动工程研究所，江苏 南京 210016）

引 言

模态 参 数 识 别 （Modal Parameter Identification，MPI）是试验模态分析（Experimental Modal Analysis，EMA）的重要内容。由于实验过程中不可避免的受到环境噪声、仪器测量误差等不确定因素的影响，从而导致模态参数识别结果中同样存在不确定性。从目前的模态参数识别方法来看，绝大多数方法将模态参数作为一个确定性参数进行识别，并未考虑其不确定性的影响。模态参数的不确定性信息可以有多种应用，例如用于区分结构模态与计算模态［1］，用于模型修正中作为加权信息使用［2］，因此一种模态参数识别方法在获得识别结果的同时，有必要给出其不确定性信息。当测量噪声满足正态分布时，模态参数的不确定性可以通过其统计特性进行描述（如方差）［3］。在模态参数不确定性的研究中，时域方法：Longman对ERA方法进行摄动分析，确定出模态参数对测量误差的灵敏度［4］；Paez使用Bootstrap方法，对模态参数统计特性进行了估计［5］；Reynders，Carden讨论了随机子空间方法中模态参数的不确定性计算［3，6］。频域方法：Guillaume分析了频响函数有理分式模型中系数摄动对极点的影响［7］；Pintelon给出了运行模态分析中模态参数不确定性的计算方法［8］；Troyer对模态参数置信区间的快速计算方法进行了研究［9，10］。

本文以频域中CMIF模态参数识别方法（Complex Modal Indicator Function）为研究对象［11］，通过Kronecker代数以及矩阵的灵敏度分析，详细推导了由测试不确定性（频响函数中的噪声方差）获得模态参数不确定性（模态参数方差）的计算公式。

1 CMIF模态参数识别方法理论基础

CMIF方法最初作为一种模态指示函数来使用，用于指示模态的分布情况，进而发展为一种模态参数识别方法，用于EMA情况下模态参数识别。经过发展，该方法已经推广到运行模态分析（Operational Modal Analysis，OMA），用于仅有响应数据情况下的参数识别［12］。CMIF识别方法的思想是利用复模态指示函数曲线，对固有频率处的频响函数进行奇异值分解，获得增强频响函数（Enhanced FRF）。增强频响函数在固有频率附近近似为单自由度系统，可运用单自由度拟合方法进行识别。由于其拟合频段范围很窄，属于一种窄带识别算法［13］，其识别过程可以按照以下两步进行。

Step1：计算固有频率处的增强频响函数

对频响函数矩阵H（jω）∈CNo×Ni（No，Ni为输出、输入数目）在固有频率附近处进行奇异值分解，可以得到

Step2：对增强频响函数进行拟合获得模态参数

式中λ为极点，待识别参数b1，b0，λ可以通过以下最小二乘方程得到

其中X，θ，Y的表达式为：

式中Nf为拟合频段内包含的谱线数。固有频率、阻尼比可以通过式（6）计算

2 CMIF方法中模态参数不确定性的计算

图1 模态参数方差的计算流程Fig．1 The calculation flow of modal parameter variance

根据CMIF方法的识别过程，模态参数方差的计算流程按照图1分4步完成。

2．1 频响函数中噪声方差的估计

噪声的方差信息可以在频响函数估计过程中获得（如H1，H2估计），对于频响函数矩阵H（jω）中的每一频响函数Ho，i（jω）（o＝1，2，…，No，i＝1，2，…，Ni），方差的估计方法如下［14］

式中 var表示方差，E表示期望，M为频响函数估计过程中的平均次数，γo，i（jω）为相干函数。

2．2 增强频响函数协方差的计算

通过线性代数知识可知［15］，对于矩阵A，B，C存在如下关系

为了以后的推导过程中表示方便，这里引入2个下标xre和xRe，二者的定义［8］

式中x可以为标量，也可以为矩阵；Re（x），Im（x）表示x的实部与虚部。容易证明以上两个符号满足如下的关系

其中cov（［vec（H（jω））］re）∈R2NoNi×2NoNi为频响函数矩阵按列展开后实部、虚部的协方差矩阵。

2．3 极点协方差的计算

对式（4）两端左乘XH得到如下正规方程

式中M＝XHX，N＝XHY，对式（17）进行一阶灵敏度分析可以得到

将式（16）中的高阶项略去，得到

式中

将式（18）代入式（17）得到Δθ的表达式

根据式（5）可知λ位于列向量θ的第3行，因此令P（jωk）为列向量（jωk－λ）M－1X（jωk）H的第3行，则Δλ可以表示为

将式（12）应用于式（20），（Δλ）re可以表示为

之所以要将Δλ表示（Δλ）re的形式，原因是式（24）中计算模态参数方差时需要使用λre的协方差矩阵，而不是λ的方差。由式（21）得到λre的协方差为

假设频响函数中噪声在不同频率分量之间相互独立［9］，式（22）可以简化为如下形式

2．4 模态参数方差的计算

由于CMIF方法在s域中识别，s域中极点的cov（λre）与模态参数方差之间关系可通过极点的灵敏度分析得到，计算公式如下，详细推导见文献［8］。

3 仿真算例

图2 5自由度仿真算例Fig．2 Five DOF simulation case

采用图2中的一个5自由度数值算例进行仿真检验。首先由表1数据计算出模态参数，如表2所示。在各阶模态中添加1%阻尼比，利用模态叠加法得到理论频响函数，频率范围0～3Hz，谱线数为1 024。在理论频响函数的基础上生成200个频响函数样本，在每个样本中添加5%的高斯白噪声，图3为一个频响函数样本以及噪声的标准差。通过Monte Carlo模拟（简称MC）对每个样本进行识别，得到200组模态参数，统计出模态参数的方差，并以此作为真值。然后由本文方法计算出模态参数的方差，将此结果与MC结果作比较。

表1 模型的物理参数Tab．1 Physical parameters

表2 模态参数的理论值Tab．2 Theoretical modal parameters

图3 一个FRF样本及噪声标准差（5%噪声）Fig．3 A sample FRF and noise standard deviation（5%noise）

图4为5%噪声下一个频响函数样本的识别结果。图中曲线从上向下依次为频响函数奇异值分解得到的奇异值曲线，用于指示模态。第一条奇异值曲线（最上面）的每个峰值对应一阶模态，如果存在密集模态，则第一、第二两条奇异值曲线会同时出现峰值。采用MC模拟对200组频响函数样本进行参数识别，结果如图5所示。

图4 一个频响函数样本的CMIF方法识别结果Fig．4 CMIF curve and identification results

图5 Monte Carlo模拟Fig．5 Monte Carlo simulation

表3为5%噪声情况下，频率、阻尼方差的MC模拟结果与本文的计算结果，其中频率标准差最大误差为6．6%，阻尼标准差最大误差为8．3%。另外在频响函数中添加3%，10%的高斯白噪声，这3种情况的MC模拟结果与计算结果如图6所示，随着噪声量级的增加，模态参数的方差也在增加。在本算例中，各阶模态随着阶次增加方差也在增大，原因是各阶模态的强度依次减弱（图4中各阶模态的峰值高度），因此受到同样量级噪声干扰的情况下，弱模态的抗噪声干扰能力低，导致其模态参数方差增加。

按照文献［9］计算了5%噪声情况下ployLSCF算法中模态参数的方差，其中识别为10，结果如表4所示。由于采用两种不同的模态参数方法，其识别结果的方差除了与噪声大小有关之外，还与识别算法自身有关，因此二者之间数值并不相同，但是结果的数量级以及各阶模态方差的相对大小基本一致。

表3 计算结果与MC模拟结果（CMIF）Tab．3 Estimation and MC simulation results（CMIF）

表4 计算结果与MC模拟结果（polyLSCF）Tab．4 Estimation and MC simulation results（polyLSCF）

图6 不同量级噪声下计算结果与MC模拟结果Fig．6 Estimation and MC results for different noise level

4 实测算例

采用图7中T型构件进行实验验证。实验由锤击法完成，数据通过法国OROS公司的OR34动态分析仪采集。由分析仪所获得数据为UFF58格式［16］，该数据格式广泛应用于结构动力学分析中。本文中的程序在MATLAB2007下完成，为了读取UFF格式文件，需要参照文档进行数据转换［16］。

图7 实验装置与测量仪器Fig．7 Experimental setup and instruments

如图8所示，测试重复进行40次，每一次敲击可以估计一组频响函数；图9为一个频响函数样本及其噪声分布。噪声的方差通过40组频响函数统计得到，模态参数方差的计算过程与仿真算例相同。

图8 一个测点多组激励与响应Fig．8 Multi－group input and output

图9 实测频响函数及其噪声标准差Fig．9 Measurement FRF and noise standard deviation

图10 T型构件一个FRF样本的识别结果Fig．10 Identification result of T shape structure

图11 计算结果与实验结果的比较Fig．11 Estimation and experiment results

图10为一个频响函数样本的识别结果，在分析范围内包含7阶模态。模态参数方差的计算结果与实验结果如图11所示，二者吻合程度相对于仿真算例有所降低，但也能够接受。误差主要来自两方面：（1）频响函数矩阵中的所有频响函数互不相关假设不能完全满足，因此计算噪声的协方差矩阵有一定误差；（2）式（23）中假设频响函数中噪声在不同频率分量之间相互独立，由图9可知，噪声与频响函数之间具有一定的相关性。

5 结 论

（1）利用矩阵的灵敏度分析与Kronecker代数理论，通过测量频响函数中噪声的方差信息，给出了CMIF识别算法中模态参数不确定性的计算方法，在识别结果的基础上增加了统计信息。

（2）仿真结果表明，推导过程中采用的一阶灵敏度分析方法已经具有较高的精度。实测算例中由于一些假设不能完全满足，计算精度相对于仿真结果有所降低。

（3）在工程应用中，如果原始数据只有频响函数，没有噪声的方差信息，这时模态参数方差的计算将受到限制。

［1］ Verboven P，Cauberghe B．User－assisting tools for a fast frequency－domain modal parameter estimation method［J］．MSSP，2004，18（4）：759—780．

［2］ Cauberghe B．Applied frequency－domain system identification in the field of experimental and operational modal analysis［D］．Belgium：Vrije Universiteit Brussel，2004．

［3］ Reynders E，Pintelon R，Roecka G．Uncertainty bounds on modal parameters obtained from stochastic subspace identification［J］．Mechanical Systems and Signal Processing，2008，22（4）：948—969．

［4］ Longman R W，Juang J N．Variance Based Confidence Criterion for ERA Identification Modal Parameters［R］．AAS 87－454：581—601．

［5］ Paez T L，Hunter N F．Statistcal series：part25：fundamental concepts of the bootstrap for statistical analysis of mechanical systems［J］．Experimental Techniques，1998，22（3）：3 523．

［6］ Carden E P，Mita A．Challenges in developing confidence intervals on modal parameters estimated for large civil infrastructure with stochastic subspace identification［J］．Structural Control and Health Monitoring，2009，23（1）：217—225．

［7］ Guillaume P，Schoukens J，Pintelon R．Sensitivity of roots to errors in the coefficients of polynomials ob－tained by frequency domain estimation methods［J］．IEEE Transactions on Instrumentation and Measurement，1989，38（6）：1 050—1 056．

［8］ Pintelon R，Guillaume P，Schoukens J．Uncertainty calculation in（operational）modal analysis［J］．MSSP，2007，21（6）：2 359—2 373．

［9］ Troyer T，Guillaume P，Steenackers G．Fast variance calculation of polyreference least－squares frequency－domain estimates［J］．MSSP，2009，23（5）：1 423—1 433．

［10］Troyer T，Guillaume P，Pintelon R，et al．Fast calculation of confidence intervals on parameter estimates of least－squares frequency－domain estimators［J］．MSSP，2009，23（2）：261—273．

［11］Shih C Y，Tsuei Y G．Complex mode indication function and its applications to spatial domain parameter estimation［J］．MSSP，1988，2（4）：367—377．

［12］Zhang L M，Wang T，Tamura Y．A frequency－spatial domain decomposition（FSDD）method for operational modal analysis［J］．MSSP，2010，24（5）：1 227—1 239．

［13］Zhang L M，Sun X H，Wang T．Narrow－band，selectband vs broadband modal identification：their features and comparisons［A］．Proc of the 26th International Modal Analysis Conference［C］．USA，February 2009．

［14］贝达特，皮尔索，著．相关分析和谱分析的工程应用［M］．令福根，译．北京：国防工业出版社，1983．

Bendat J S，Piersol A G．Engineering Applications of Correlation and Spectral Analysis［M］．Beijing：National Defence Industry Press，1983（in Chinese）．

［15］戴华．矩阵论［M］．北京：科学出版社，2001．Dai H．Matrix Theory［M］．Beijing：Science Press，2001．

［16］UC－SDRL．Universal File Formats for Modal Analysis Testing［EB／OL］．http：／／www．sdrl．uc．edu／universal－file－formats－for－modal－analysis－testing－1，2008．