基于声压测量的结构模态参数辨识

夏茂龙， 黎 胜,2

(1.大连理工大学 工业装备结构分析国家重点实验室 船舶工程学院,大连 116024；2.高新船舶与深海开发装备协同创新中心,上海 200240)

基于声压测量的结构模态参数辨识

夏茂龙1， 黎 胜1,2

(1.大连理工大学 工业装备结构分析国家重点实验室 船舶工程学院,大连 116024；2.高新船舶与深海开发装备协同创新中心,上海 200240)

通过建立辐射声压与激振力之间的声压频响函数矩阵，提出了一个基于测量声压识别振动结构模态参数的方法。该方法可以通过非接触测量声压来识别结构的固有频率、模态阻尼比与模态振型，避免了附加质量对结构的影响。声压频响函数矩阵是基于边界元和Rayleigh积分方法结合有限元结构动力学方程建立的，适用于任意结构且与结构模态参数有明确的关系；对于测量声压时多激励单输出与单激励多输出响应的不同试验模式，该方法都能识别结构的模态参数。以一平板结构为例，数值验证了该方法的准确性与适用性。

声压测量；模态参数；声压频响函数；振动声辐射

结构振动及辐射噪声可通过构建的模态分析模型来表述，若建立了准确的分析模型，可解决结构优化控制、结构设计、响应预测等一系列问题，因此如何获得表征模型的模态参数是一个重要的基础和前提[1]。试验模态分析作为一个传统的实验方法通过测试结构系统的激励与响应如位移、速度、加速度信号，建立系统的输入与输出之间的传递函数来识别结构的模态参数[2]。在测试过程中：①测量输出响应的传感器如加速度计等会对结构产生附加质量的影响，特别是对于对称结构或者密集模态结构，附加质量会对结构识别的模态参数产生较大影响；②测量输出响应的传感器须贴附在结构的表面，当结构表面处于高温等恶劣环境、结构表面无法接近或者无法贴附传感器时，则无法通过传感器来测量结构的响应。虽通过激光测振仪、高速摄像机等方法可通过非接触测量来获得结构的响应，但是仪器较为昂贵，也不能被广泛应用[3]。而基于声压的测量是一种非接触测量方法，其与结构的激励与响应相关，同时可以避免附加质量、结构表面不可测等问题，基于此本文探讨了利用测量的声压信号来识别结构模态参数的方法。

国内外学者针对声压测量来获得结构振动响应也做了很多研究，如近场声全息方法(Inverse Numerical Acoustics，NAH)[4-5]、声学逆问题求解方法(Noise,Vibration and Harshness,INA)[6]等可以通过声压测量阵列来实现结构表面振动响应的重建，但其需要辐射声场的近场声压阵列信息，且很少对结构的模态分析进行探讨。Prezelj等[7]利用表面近场的测量声压实现了振动结构模态重构，但是其要求所测量的声压须为距离结构振动表面非常近的近场声压信息。Zhu等[8-9]通过有限元方法建立声固耦合的动力学方程，并以此实现了结构的模态参数的辨识，然而此方法需建立声固耦合系统的有限单元模型，对于复杂的结构或者声场测量时，计算量非常大。国内学者基于声压测量的研究主要针对结构噪声控制，对模态参数识别方面的研究较少。

而利用边界元方法或者Rayleigh积分方法，可直接建立结构振动响应与辐射声场响应之间的关系[10-12]。由此，本文采用边界元方法和Rayleigh积分方法，无需重构结构表面振动信息，建立了结构激励与声压响应之间的声压频响函数，推导出辐射声场与结构模态参数之间的关系，进而实现了模态参数(固有频率、阻尼比、模态振型)的识别。该方法通过非接触测量辐射声场中的响应，可以避免加速度计等附加质量的影响；采用的边界元方法和Rayleigh积分方法可以仅离散表面的单元，计算简便且计算量小。若辐射声场中声压可测量得到，则可通过测量辐射声场中任一适当位置的声压信息来实现结构的模态参数辨识。对于不同的试验模式，该方法都可以识别结构的模态振型，对于单激励多输出系统来说，若识别了声压模态振型，不仅可以实现振动模态振型的识别，而且可通过声压模态振型表述振动产生的辐射声场。最后以一空气中的平板结构为例，通过数值分析测量结构的声压响应得到了结构的声压频响矩阵，分别计算了不同试验模式下识别的固有频率、模态阻尼比及模态振型，得到的结果与数值解吻合，并讨论了随着距离的变化，识别的声压振型随距离变化的影响。

1 结构振动产生的辐射声场

对于封闭表面结构振动产生的辐射声场问题，可看做给定的结构表面声压或法向速度在空间中产生的声场，该声场可以通过Kirchoff-Helmholtz方程得到，若P为场点中任意一点，Q为结构封闭表面S上任意点，可建立结构的法向振速与声场中任意一点的声压的关系。即

(1)

式中：G(Q,P)=e-ikL/4πL为自由空间格林函数；p(Q)为点Q处的声压值；L=|Q-P|为声压场点P与结构表面法向振速点Q之间的距离；k=ω/c为波数；c为声速；n为表面S的外法线方向；C(P)为边界系数；取决于点P的位置。

对Kirchoff-Helmholtz方程进行离散，可得到结构表面单元的声压向量{PS}与结构表面法向振动速度{Vns}的公式

[Es]{Ps}=[Ds]{Vns}

(2)

式中：[Es]、[Ds]为系数矩阵。通过离散后的Kirchoff-Helmholtz方程，声场中任意一点的声压可用结构表面声压与结构表面法向振速表示

{Pf}=[Ef]{Ps}+[Df]{Vns}

(3)

式中，[Ef]、[Df]为系数矩阵。

将式(2)与式(3)联立，可以得到声场中任意一点声压P与结构表面振动法向速度之间的关系，即

{Pf}=([Ef][Es]-1[Ds]+[Df]){Vns}=[Dv]{Vns}

(4)

对于无限大平面障板结构，结构的振动表面可以离散成许多单元表面dS。则结构振动产生的辐射声场中的声压可以通过Rayleigh积分得到，即声场中的声压Pf可以表示为

(5)

式中：ρ0为介质的密度；ω为结构表面振动频率。将Rayleigh积分离散，可得到声场中任意一点声压Pf与离散的结构表面单元法向速度{Vns}之间的关系，即

{Pf}=[Dv]{Vns}

(6)

而对于结构表面法向速度{Vns}与结构表面速度{V}之间的可以通过转换矩阵G给出

{Vns}=[G]T{V}

(7)

式(4)、式(6)又可写作

{Pf}=[Dv]{Vns}=[Dv][G]T{V}=[R]{V}

(8)

式中，[R]可以定义为离散形式的辐射算子。

对于未封闭表面的三维结构，可通过间接边界元方法建立结构振动与产生的辐射声场的关系。综上所述，对于任意形状的结构而言，都可以由式(8)建立结构振动辐射声场中任一点的声压Pf与结构表面振动速度{V}之间的关系，据此可建立和分析下节中声压频响函数矩阵，得到结构振动激励与振动辐射声场中声压的关系。

2 结构振动的频响函数矩阵

若忽略结构与声场的耦合作用，当结构受到简谐激振力Feiωt作用时，结构的动力学方程可以表示为

(9)

则频响函数矩阵H

(-ω2M+iωC+K)=H-1(iω)

(10)

在上述方程的两端同时左乘ΦT，右乘Φ，由于振型Φ的正交性，则式(10)左端化为对角阵，即

(-ω2Mr+iωCr+Kr)=ΦTH-1(iω)Φ

(11)

两端取逆矩阵，得到结构振动的位移频响函数

H(iω)=Φ(-ω2Mr+jωCr+Kr)-1ΦT=Φ(S)-1ΦT(12)

式中，S为对角阵，其对角线元素Srr

(13)

当在j点加单位激振力，i点测量响应的幅值时，

(14)

式(14)可以写成更一般的形式

(16)

即通过式(15)得到多自由度结构的频响函数矩阵H(iω)。

3 结构辐射声场中声压频响函数矩阵与模态参数识别

若测得在施加结构j点的激振Fj，结构振动辐射声场中i点的声压Pi，则可定义声压传递矩阵

(17)

由式(8)可知，结构振动辐射声场中的任意一点的声压可由结构表面振动的速度获得，在简谐激励作用下

(18)

若结构自由度为N，结构的单元数为n，则式(18)可写作

(19)

式中，Rim为矩阵R第i行第m个元素，继续化简可得到

(21)

式中，ψir为第r阶模态振型振动时辐射声压响应，可定义为第r阶声压模态振型，其与结构的激振频率ω、声压位置i对应的辐射算子R，第r阶模态振型φr相关。

为进一步确定各个参数之间的关系，以离散的Rayleigh积分为例，化简可得到

即结构的声压模态振型的幅值与结构的振动频率ω2，介质密度ρ0成正比，当结构表面单元确定后，声压模态振型与振动模态振型φmr，声压点与结构表面单元距离L相关。

则结构振动产生的辐射声场中i点的声压也可表示为

(23)

式中，qr为第r阶的结构模态坐标；通过式(23)可以看出辐射声场中任一点的声压可以通过声压模态振型与结构的模态坐标求得。

因此n个点声压与m个激振力之间的声压频响函数矩阵为

式中，HAi,j由式(20)、式(21)表示为

(25)

进一步的，通过表达式(24)、式(25)可以得到在特定频率下，结构振动产生的辐射声压与结构的激振力成正比，这也与Zhu等的结论相同，验证了文献[13]的假设，另外，文献[14]假设在非常近场的情况下，声压的振幅与结构振动的振幅成比例。在此，在固定的频率下，若声压的位置i与结构振动表面单元j的位置非常近，由式(22)可得，

(26)

由式(26)得到结构的声压振幅与结构振动幅值成线性关系[14]。

基于频响函数矩阵分析结构模态参数的方法有很多，如峰值提取法、有理分式多项式法、正交多项式法、PolyMax法等[15]。式(25)与式(15)对比可知，声压频响函数矩阵与振动频响函数矩阵具有相同的形式，两式都包含的相同极点sr。因此，可以通过模态分析法来获得结构的固有频率、模态阻尼比、模态振型。

值得注意的一点是，振型矢量通常可以通过对一系列响应测点求出的留数得到。而声压频响函数虽然与振动频响函数具有相同形式的留数矩阵，但是其包含的声压模态振型与振动模态振型不同，识别的结果也有所不同。

声压频响函数的留数定义为

式中：ψmr为声压模态振型第r阶第m个元素；φnr为振动模态振型第r阶第n个元素。

(1)当结构进行任意n点激励，测得任意一点j的声压响应输出时，即通过声压频响函数矩阵获得声压留数矩阵Br的第j行，由式(27)通过计算得到振动模态振型φr。

(2)当结构受任意一点i激励，m个声压响应值输出时，即通过声压频响函数矩阵获得留数矩阵Br的第i列，由式(27)通过计算得到声压模态振型ψr。

因此，根据不同的试验模式，我们可以分别得到结构的振动模态振型φr与声压模态振型ψr。

进一步，若无法通过(1)直接获得结构振动模态振型φr，而只能通过(2)得到结构的声压模态振型ψr，则式(21)写成矩阵的形式

ψr=Rpφr,Rp=iωR

(28)

所以

(29)

但是在声学逆问题求解过程中，仍存在很多问题[16]。当Rp的条件数很大即为病态矩阵时，若ψr引入很小的噪声干扰，就会引起φr的扰动。而这在试验模态分析中测量是难以避免的，当Rp为病态矩阵时直接由式(29)计算的模态振型φr是不准确的。

对此很多学者也做了研究。若Rp可随试验条件改动，则可通过增加声场测点的位置，来丰富声场信息从而降低矩阵Rp的条件数，再通过式(29)来直接计算结构的模态振型。若Rp不可改动，则可通过截断奇异值分解(Truncated Singular Value Decomposition,TSVD)[17]或者正则化方法[18]来降低结构的“病态性”，从而计算得出结构的模态振型。

4 数值计算与讨论

4.1 模型的的建立

以四边简支的矩形钢板在空气中的振动声辐射为例进行数值计算，平板板长Lx=0.455 m，板宽Ly=0.379 m，板厚h=0.003 m，板密度ρs=7 850 kg/m3，弹性模量E=2.1×1011N/m2，泊松比υ=0.3，空气密度ρ0=1.21 kg/m3，声速c=343 m/s，平板结构设为瑞利阻尼，阻尼系数基于结构的前两阶频率ω1、ω2和阻尼比0.01得到：α=7.362，β=1.177×10-5。

计算中采用四边形四节点等参单元将平板单元离散为16×16的单元，单元节点数为289，因四边简支，边界上的每个节点有2个转动自由度，中间的节点每个有1个垂向平动自由度和2个转动自由度，故自由度数为803。有限元是基于Mindlin板弯曲理论，板单元Rayleigh积分的离散与有限元计算同样的网格。可通过有限元数值计算得到结构前6阶的固有频率与模态振型，模态阻尼比可由瑞利阻尼系数与固有频率得到。

在试验模态数值计算过程中，主要在两种试验模式下计算声压频响函数，分别为单激励输入多声压输出系统与多激励输入单声压输出系统，具体如下：①激振力作用于编号40(x=0.284 m,y=0.071 m)的位置，大小为1 N;测量距离平板z=0 m，z=0.01 m，z=0.05 m，z=0.08 m，z=0.1 m，z=0.15 m，z=0.2 m的所有节点位置处的声压幅值，通过式(17)、式(24)获得随距离、频率变化的声压频响函数矩阵。②对板表面除边界外的编号位置进行激励，测量编号40，距离平板z=0 m,z=0.01 m,z=0.05 m,z=0.08 m,z=0.1 m,z=0.15 m,z=0.2 m位置处的声压幅值，通过式(17)、式(24)获得获得随距离、频率变化的声压频响函数矩阵。

采用有理分式多项式法对获得的声压频响函数矩阵进行拟合，获得结构的频率和模态阻尼比。对于多激励输入单声压输出系统的试验模式来说，可计算出结构振动模态振型，对于在节点或者节点线上测点的结果因其无法准确获得其对应阶的频响函数，故所识别的对应阶的频率、模态阻尼比、模态振型系数不予采用。对于单激励输入多声压输出的的试验模式来说，可由式(25)、式(26)获得结构的固有频率、模态阻尼比与声压模态振型。同样的若测点在节点或者节点线上，则所识别的对应阶的频率、模态阻尼比、声压模态振型系数不予采用。进一步的，若通过结构表面的离散单元计算求出Rp，当z=0.08 m时，Rp的条件数为1.58×107，此时无法通过式(29)通过声压模态振型准确计算出振动模态振型。因此当z≥0.08 m时采用TSVD法来改善Rp矩阵的性态。 平板结构及测点位置见图1。

图1 平板结构及测点位置Fig.1 The plate structure and positions of measurement

4.2 计算结果与讨论

通过数值计算分别为单激励输入多声压输出系统与多激励输入单声压输出系统的声压频响函数矩阵，再通过拟合进行结构的模态参数识别。由于结构的固有频率和阻尼比仅通过适当一点的声压频响函数即可识别，本节采用了编号40(x=0.284 m,y=0.071 m)、49(x=0.144 m,y=0.095 m)、152(x=0.057 m,y=0.284 m)位置的测量数据求出结构的固有频率与阻尼比并取平均值。以距离表面z=0.01 m,z=0.08 m,z=0.2 m为例，识别的结果与采用有限元方法计算的结果如表1～表3，得到的结构模态振型以z=0.08 m为例如表4。将求得的结构模态振型与有限元分析得到的模态振型进行模态置信度准则(Modal Assurance Criterion，MAC)分析，以此来判断识别的振型与理论振型的准确性，得到MAC矩阵如表5。

由MAC分析可知，矩阵的对角元素大于99%，非对角元素远远小于主对角元素，验证了识别振型的准确性。

表1 z=0.01 m模态识别结果与有限元计算结果对比

表2 z=0.08 m模态识别结果与有限元计算结果对比

表3 z=0.2 m模态识别结果与有限元计算结果对比

表4 模态振型结果与有限元计算结果(z=0.08 m)

表5 识别模态振型与有限元分析模态振型的MAC矩阵(z=0.08 m)

测量距离对于模态参数识别的结果有一定的影响。且随着距离的增加，识别的声压模态振型也会有所不同，在此同样以平面钢板为例，探讨了随距离变化，对声压模态振型的识别结果如表6。

通过表6可知，随着距离的增加，所识别的模态声压振型也有所不同，即对于单模态振动而言，当距离表面较近时，结构的声压模态振型与结构的模态振型成比例关系；当距离增加时，结构的声压模态振型随着波阵面扩展，所识别的声压模态振型结果同振动模态振型也会有较大变化；当距离gt;0.5 m时，对于奇-奇振动模态而言，在观测面上不再有明显的振型区别，对于偶-偶模态而言，在观测面上只有边角处有较大数值。因此，在传声器布置位置时，需考虑距离的变化对声压模态振型的影响，选择适当的传声器测量结果来进行模态参数的识别。

表6 不同距离下前6阶声压模态振型识别结果

以一平板结构的数值计算结果分析，在多激励输入单声压输出与单激励多输出的试验分析中，表1、表2、表3的结果表明，识别的前6阶固有频率、模态阻尼比与有限元分析的模态固有频率与阻尼比吻合，且随着与平板距离的增加，传声器在z=0.01 m,z=0.08 m,z=0.2 m上也同样能识别准确的结构固有频率与模态阻尼比。之所以能准确地识别前6阶模态参数在于：所选的观测平面位于平板结构的正上方，且与平板的激振方向一致；所选测点的位置避开了振动模态振型或声压模态振型的节点与节点线，可以同时得到多阶模态振动响应产生的辐射声压；同样随着距离的增加，所测声压虽小，但没有其他噪声的干扰，因此能得到各阶模态振动声辐射的声压响应。若所选的测点位置不当或考虑噪声的影响而无法获准确得对应阶模态振动辐射声场的响应，则影响此方法识别的精度。

以z=0.08 m为例，由表5、表6可知，在多激励输入单声压输出与单激励多单声输出的试验分析中，所识别的振动模态振型与有限元分析结果吻合，经MAC矩阵的验证，矩阵的主对角元素大于99%，这也验证了所识别振动模态的准确性。当z≤0.08 m，所识别的声压振型与所识别的模态振型相似，这也数值验证了Elwali等的假设。随着距离的增大，在观测面的声压模态振型有了较大变化，因此在传声器布置时同样要考虑距离的影响。

5 结 论

本文基于边界元和Rayleigh积分方法结合有限元结构动力学方程，推导出结构的声压频响函数矩阵，建立了声压频响函数矩阵与结构模态的关系，并基于此提出了基于振动产生的辐射声压来识别振动结构的模态参数(固有频率、阻尼比、结构模态振型)的方法，该方法可以通过非接触测量振动产生的辐射声压来识别，避免了附加质量对结构的影响。对于声压测量的多激励单输出与单激励多输出响应的不同试验模式，该方法都能识别结构的固有频率、模态阻尼比与模态振型。最后以一平板结构为例，数值验证了该方法的准确性与适用性。数值结果表明，该方法能准确地识别平板结构的固有频率、模态阻尼比与模态振型，验证了基于声压测量的模态参数识别方法的准确性，具有一定的工程实用价值。

[ 1 ] 傅志方, 华宏星. 模态分析理论与应用[M]. 上海：上海交通大学出版社, 2000.

[ 2 ] EWINS D J. Modal testing: theory and practice[M]. Letchworth: Research Studies Press, 1995.

[ 4 ] WILLIAMS E G. Regularization methods for near-field acoustical holography[J]. The Journal of the Acoustical Society of America, 2001, 110(4): 1976-1988.

[ 5 ] MAYNARD J D, WILLIAMS E G, LEE Y. Near field acoustic holography: I. Theory of generalized holography and the development of NAH[J]. The Journal of the Acoustical Society of America, 1985, 78(4): 1395-1413.

[ 6 ] TOURNOUR M, CREMERS L, GUISSET P, et al. Inverse numerical acoustics based on acoustic transfer vectors[C]//7th International Congress on Sound and Vibration. Germany:Garmisch-Partenkirchen,2000: 2069-1076.

[ 8 ] ZHU W D, LIU J M, XU Y F, et al. A modal test method using sound pressure transducers based on vibro-acoustic reciprocity[J]. Journal of Sound and Vibration, 2014, 333(13): 2728-2742.

[ 9 ] PIERRO E, MUCCHI E, SORIA L, et al. On the vibro-acoustical operational modal analysis of a helicopter cabin[J]. Mechanical Systems and Signal Processing, 2009, 23(4): 1205-1217.

[10] Boundary element acoustics: Fundamentals and computer codes[M].[S.l.]：Wit Pr/Computational Mechanics, 2000.

[11] FAHY F J, GARDONIO P. Sound and structural vibration: radiation, transmission and response[M]. Academic Press, 2007.

[12] HAMDI M A, VILLE J M. Development of a sound radiation model for a finite-length duct of arbitrary shape[J]. AIAA Journal, 1982, 20(12): 1687-1692.

[13] ELWALI W, SATAKOPAN H, SHAUCHE V, et al. Modal parameter estimation using acoustic modal analysis[M].New York:Springer, 2011.

[14] XU Y F, ZHU W D. Operational modal analysis of a rectangular plate using non-contact excitation and measurement[J]. Journal of Sound and Vibration, 2013, 332(20): 4927-4939.

[15] PINTELON R, SCHOUKENS J. System identification: a frequency domain approach[M].Hoboken: John Wiley amp; Sons, 2012.

[16] NELSON P A, YOON S H. Estimation of acoustic source strength by inverse methods: Part I, conditioning of the inverse problem[J]. Journal of Sound and Vibration, 2000, 233(4): 639-664.

[17] HANSEN P C. Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numerical rank[J]. SIAM Journal on Scientific and Statistical Computing, 1990, 11(3): 503-518.

[18] LIU W, WU C. A predictor-corrector iterated Tikhonov regularization for linear ill-posed inverse problems[J]. Applied Mathematics and Computation, 2013, 221(1): 802-818.

Identificationofstructuralmodalparametersbasedonsoundpressuremeasurement

XIA Maolong1, LI Sheng1,2

(1.State Key Laboratory of Structural Analysis for Industrial Equipment, School of Naval Architecture,Dalian University of Technology, Dalian 116024, China;2. Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, Shanghai 200240, China)

Based on the frequency response function of sound pressure, a method of identifying structural modal parameters from measured sound pressure was proposed. It can accurately identify the structural natural frequencies, damping ratios and mode shapes from the measured sound pressure, and can avoid mass loading due to the use of vibration sensors on the structure in traditional measurements. Combined with structural dynamic equations, the sound pressure frequency response function was established based on the boundary element method and the Rayleigh integral formulation for sound radiation from an arbitrary body. The method was also able to identify the structural modal parameters with Single-Input-Multiple-Output testing and Multiple-Input-Single-Output testing. Finally the accuracy and the validity of the method were verified through numerical simulations.

sound pressure measurement; identification of modal parameters; frequency response function of sound pressure; vibration and sound radiation

辽宁省教育厅重点实验室基础研究项目(LZ2014004)

2016-04-11 修改稿收到日期：2016-07-18

夏茂龙 男，博士生，1988年生

黎胜 男，博士，教授，博士生导师，1973年生

TU311.3

A

10.13465/j.cnki.jvs.2017.22.036