基于Matlab的倍频程分析.docx
《基于Matlab的倍频程分析.docx》由会员分享，可在线阅读，更多相关《基于Matlab的倍频程分析.docx（16页珍藏版）》请在冰豆网上搜索。
基于Matlab的倍频程分析
用Matlab语言实现噪声的1/3倍频程分析
摘要:
在声学测量研究中，1/3倍频程谱反映了声源的能量分布情况.本文基于Matlab
软件开发平台，实现了对高斯白噪声的1/3倍频程分析，验证了该算法的正确性，具有精度高,性能稳定的特点。
关键词:
故障诊断；传感器优化布置；高斯白噪声；功率谱；1/3倍频程
0引言
“设备故障诊断（ConditionMonitoringandFaultsDiagnosis）”是近十几年发展起来的一门新兴技术，包含两方面的内容：
一是对设备的现场运行状态进行监测；二是在出现故障情况时对故障进行分析与诊断。
这二者是密不可分和相互关联的。
掌握设备现在的状况及信息，预知和预测有关故障或异常的程度，分析故障产生的原因，判断故障发展趋势及其对将来的影响，从而找出必要的对策或解决方法，是设备故障诊断的功能。
运用设备诊断技术所取得的经济效益是明显的。
据日本资料报道，采用诊断技术后，事故率减少75%，维修费用可降低25-50%。
英国对两千个工厂进行的调查表明，采用诊断技术后，维修费用每年可节约3亿英磅。
目前我国的机器设备总值约为8000多亿元，每年用于设备大修、小修和处理故障的费用一般占固定资产原值的3%—5%，采用诊断技术改进维修方式和方法后，一年取得的经济效益可达数百亿元。
因而减少停机时间而创造的社会效益将非常巨大。
显然，设备故障诊断与监测技术对企业的正常生产经营是必不可少的，必须把它作为企业管理与发展的一个重要内容。
设备故障诊断一般分两个阶段四个步骤实施。
两个阶段为状态监测和故障诊断。
故障诊断的四个步骤为：
信号检测、特征提取（信号处理）、状态识别和诊断决策。
其具体内容为：
（l）信号检测：
按不同诊断目的选择最能表征工作状态的信号。
这种工作状态信号称为初始模式。
（2）特征提取（信号处理）：
将初始模式向量进行信号处理、变换，去掉冗余信息，提取故障特征、形成待检模式；
（3）状态分析：
将待检模式与样式模式（故障档案）对比和状态分类，判断出故障类型。
（4）诊断决策：
根据判别结果采取相应对策，对机械设备及工作进行必要的预测及修正。
1传感器优化布置
结构健康监测检测是近年来发展起来的结构无损检测技术。
它利用结构的某些信息,运用一定的数学方法,来判定结构是否损伤以及损伤的位置和程度。
近10年来,国内外学者普遍认同的损伤评估方法是试验模态分析法。
进行模态实验的第1步就是获得被测结构激励和响应的时域信号,而传感器的配置方案是首先要确定的。
不适当的传感器配置将影响识别参数的精度,而且传感器本身需要一定的成本,与其配套使用的数据采集和处理设备的代价也都较高。
从经济方面考虑,希望采用尽可能少的传感器。
因此,确定传感器的最佳数目,并将它们配置在最优位置,具有重要的实用价值。
要进行传感器的优化配置,首先要确定合理的、能反映设计要求的优化配置准则。
目前发展起来的优化准则很多,其中,基于识别误差最小准则的方法是使用较多的1种方法;如:
Kammer提出的有效独立法;模型缩减准则也是1种常用的准则,但这种方法只能保证低阶模态的精度;另外还有可控可观度准则,模态应变能准则等。
Came等认为模态置信度MAC（ModalAssuranceCriterion）矩阵是评价模态向量交角的1个很好的工具。
其次,传感器的优化配置还必须选用适当的优化方法。
近几年发展起来的随机类方法主要有模拟退火算法和遗传算法,这种方法不易陷入局部最优解,但现有的随机搜索技术可靠性并不高。
目前使用最多的是序列法中的逐步削减法,它每次从剩余传感器的可选位置中去掉1个或多个对目标函数贡献最小或较小的可选位置,直到剩余最优可选位置为止。
逐步累积法与削减法相反,它是不断的从剩余可选位置中选取1个最优的加入到优化配置中,直至达到最优的数目为止。
本文在此基础上提出了修正的逐步累计法,算例表明,修正的逐步累计法具有较高的计算精度和优化结果。
1.1基本理论
1.1.1基于QR分解的传感器配置原理
根据模态叠加原理,系统的响应可以表示为:
（1-1）
其中其中,{u}为物理坐标,{u}∈
为第i阶模态向量,
∈
为模态向量矩阵,
为第i阶模态坐标,{q}∈
。
s代表传感器数量,m为所需识别的模态个数。
式
（1）的最小二乘解为
（1-2）
如果考虑测量噪声,式
（1）应改写为
（1-3）
其中,v代表方差为
的高斯分布白噪声,这里假设测量噪声相互独立并且对各个传感器测量信号的统计特性相同,则{
}与{q}的协方差[2]为:
（1-4）
其中,[Q]称为Fisher信息矩阵。
当[Q]取极大值时,协方差[P]最小,就能够得到较好的估计。
所以必须使[Q]的某1种范
数最大,这里选取常用的2-范数‖Q‖2。
由于
（1-5）
所以
（1-6）
因此,以上对[Q]的要求可通过
的选择来实现。
根据矩阵理论,列主元QR分解是选取
矩阵列向量组具有较大范数子集的1种简捷有效的方法。
设有限元模型所得的模态向量矩阵对应于可测自由度的子集为5,5∈
。
一般有m<
n,并且r（
）=m,即矩阵
列满秩。
由于列主元QR分解选择的是列向量组的子集,所以,进
行
的列主元QR分解:
（1-7）
其中,E为置换矩阵,则
中对应于{
1},{
2},…,{
n}的行（即自由度）就是
的行向量组中具有较大范数的子集（其中{
}代表R矩阵的第i列）。
1.1.2模态置信度MAC矩阵
由结构动力学原理可知,结构各阶模态向量在节点上的值形成了1组正交向量,但由于量测自由度远小于结构模型的自由度并且受到测试精度和测量噪声的影响,测得的模态向量已不可能保证其正交性,在极端的情况下甚至会由于向量间的空间交角过小而丢失重要的模态。
因此,在选择测点时有必要使量测的各模态向量保持较大的空间交角,从而尽可能地把原来模型的特性保留下来。
Kammer的EI法事实上从另1个角度起到了这个作用[。
Came等认为模态置信度MAC矩阵是评价模态向量空间交角的1个很好的工具,其公式表达如下:
（1-8）
其中
和
分别为第i阶和第j阶模态向量。
由式（1-8）可以看出,MAC矩阵考虑的是模态向量矩阵列空间的度量特性。
而根据1.1.1中的讨论,利用QR分解的传感器配置实际上讨论的是模态向量矩阵行空间的特性,所以虽然这种传感器配置可以保证q估计值的质量,但并不一定能够得到良好的MAC,必须采取新的措施来提高MAC以满足振型匹配的要求。
由式（1-8）可以看出,MAC矩阵的非对角元
（i≠j）代表了相应2模态向量的交角状况。
换句话说,当MAC阵的某一非对角元
（i≠j）=1时,表明了第i向量与第j向量交角为0,2向量不可分辨;而当
（i≠j）=0时,则表明第i向量与第j向量相互正交,2向量较易识别,故测点的布置应力求使MAC的非对角元向最小化发展。
2噪声的物理量度
2.1.1噪声的基本概念
各种频率和声强杂乱无序组合的声音叫做噪声。
声音是由物体振动引起的。
物体振动通过媒质中传播引起人耳或其它接收器的反应，就是声音。
振动的物体是声音的声源，产生噪声的物体或机械设备等成为噪声源。
振动在弹性介质中以波的形式传播，这种弹性波叫做声波。
当振源频率在20～20000Hz之间时，人的耳朵可以感受到它。
当振源频率低于20Hz或高于20000Hz时，人耳无法听到。
低于20Hz的波动叫次声波，高于20000Hz的波动叫超声波。
人耳日常听到的声音，通常来自空气所传播频率在20～20000Hz之间的声波。
多年来噪声信号都是作为有害信号在系统设计或改进时加以抑制。
在多数以控制为目的的测量中，传感器的设计都是尽量排除噪声信号。
其实，噪声信号中带有大量设备运行状态信息，通过对噪声信号进行分离和分析，可以了解设备运行状态，对设备进行状态检测和故障诊断。
近场噪声分析方法是一种非接触测量方法，它是利用声音传感器接近被测部件来测取噪声信号，通过分析噪声信号进行频谱分析，比较被测信号频谱特性与正常运行状态下的频谱特性来诊断机械故障。
2.1.2声压和声压级
当没有声波时，空气处于静止状态时，其压强为大气压强P0，当有声压存在时，局部空气产生压缩或膨胀，在压缩的地方压强增加，在膨胀的地方压强减少，这样就在原来大气压上又迭加了一个压强的变化量。
这个迭加上去的压强变化与静压强的差值称为声压，用Pe表示。
声压可正可负是围绕钧值摆动的量，其均值恒为零，通常仪器检测的声压为声压的均方根值Pe，称为有效声压。
（Pa）（2-1）
其中T为周期，根据简谐波
，Pm为声压幅值。
声压的大小表示了声波的强弱。
由于正常人耳能听到的最弱声音的声压和能使人耳感到疼痛的声音的声压大小之间相差一百万倍，表达和应用起来很不方便。
同时，实际上人耳对声音的感受也不是线性的，它不是正比与声压绝对值的大小，而是同它的对数近似成正比。
因此如果将两个声音的声压之比用对数的标度来表示，那么不仅应用简单，而且也接近于人耳的听觉特性。
这种用对数标度来表示的声压称为声音级，它用dB表示。
即
（2-2）
其中，Lp为声压级，单位为分贝（dB），P0是参考声压级，国际上规定P0＝20x10
Pa，这就是人耳刚能听到的1KHz的纯音的声压值。
当采用声压级的概念后，听阈与痛阈的声压比从100万倍的变化范围变成0～120dB的变化，所以“级”的大小能衡量声音的相对强弱，声压级被定义成声压平方比值的对数值。
由于声压平方的比值也与声功率成比例，因此声压级与声功率级就联系起来了。
2.1.3声强与声强级
在声波传播方向上单位时间内垂直通过单位面积的声能量，成为声音的强度或简称声强，用I表示，单位W/
。
声强的大小可用来衡量声音的强弱，声强越大，我们听到的声音越强；声强越小，我们听到的声音越轻。
声强与离开声源的距离有关，距离越远，声强越小。
根据定义，瞬时声强I可表示为
I（t）＝P（t）*u（t）（2-3）
式中P（t）为瞬时声压，u（t）为质点的振速，由于u（t）是矢量，所以I（t）也是矢量，即代表声能的传播方向。
与声压一样，声强用“级”来表示，即声强级
它的单位也是分贝（dB）,即
（2-4）
其中
为基准声强，
＝
W/
，它相当于人耳能听到的最弱的声音的强度。
声强级与声压级的关系是：
式中ρ为空气密度，c为空气中声速。
媒质的ρc随媒质的温度和气压而改变。
如果在测量条件时恰好ρc＝400，则
＝
。
对一般情况，声压级与声强级相差一个修正项
，此值是比较小的。
3幅度谱分析方法
噪声控制中的数字化技术是指对噪声信号的检测,变换,记录,传输和分析处理采用数字编码信号,通过数字计算方法对信号进行滤波,检测,快速傅立叶变换等,用以提高噪声信号的处理速度和精度。
3.1傅立叶变换
随机振动,可以认为是一个具有零均值的各态历经的平稳高斯过程1,"因此现场测试所得的振动时程曲线可以视为许多不同频率的正弦波叠加的结果,即可以用傅立叶级数的形式表示"有限时间内非周期连续函数的傅立叶变换为
正变换:
逆变换:
实际的采样信号总是离散的并且信号长度是有限的,对等时距的离散点进行
傅立叶分析称为离散点的傅立叶变换:
（n=0,1,2,…,N-1）
式中,N为采样数据量。
在噪声测量中，采用1/3倍频程谱分析能详细的反映出噪声源的频谱特性，便于较全面的了解声源产生激励和提出最佳的降噪对策。
本文就是基于上述背景和目的，应用Matlab语言实现了1/3倍频程分析功能，并将结果以柱状图的形式可视化，给人以直观的视觉感受。
3.21/3倍频程的实现
一般噪声的频率分布宽阔，在实际的频谱分析中，不需要也不可能对每个频率成分进行具体分析，为了便于观察噪声信号宏观上的能量分布，忽略信号频率或相位信息微小变化对观察结果的影响，把20～20000Hz的声频范围分为几个段落，每个段落称为频带或频程，频带的中心频率为：
，
和
为该频带的上限和下限频率。
频带的划分也有规定，一般规定集中m倍频程，m由下式确定：
当m＝1时，称为倍频程，类推，没＝1/3时，称为1/3倍频程。
在此频带（
～
）的频谱，称为1/3倍频程谱。
按照1/3倍频程的方法，可将声频范围分为更多的频带，便于较仔细地研究。
对于噪声来讲，由于1/3倍频程能够很好地体现噪声带宽地能量分布情况，为噪声控制提供参数，采取合理地措施，从而达到降噪地目的，所以在噪声分析过程1/3倍频程地分析显得尤为重要。
我国对1/3倍频程地频率范围及其中心频率都做了规定，程序中按照上、下限的规定将功谱中的频率划分成多个频带,在每个频带中将所有的功率谱幅值相加、平均,这样就借助功率谱实现了1/3倍频程功能,最终结果以柱状图的形式显示。
这里需要注意的是，在求和的过程中要用原始的测量值（噪声为Pa），如果要使用声压级dB（A）,需要转化为Pa，噪声dB（A）与Pa的关系：
dB=20xlog（Pa/20-5）（在空气中），所以在叠加时，一定要采取对数加法，否则，就会导致算法错误。
1/3倍频程的实现流程如图1所示（请见下页）：
三分之一倍频程的中心频率和频率范围如表1所示：
图11/3倍频程的实现流程
表11/3倍频程的中心频率和频率范围（Hz）
中心频率
频率范围
中心频率
频率范围
25
22.4~28
800
710~900
31.5
28~35.5
1000
900~1120
40
35.5~45
1250
1120~1400
50
45~56
1600
1400~1800
63
56~71
2000
1800~2240
80
71~90
2500
2240~2800
100
90~112
3150
2800~3550
125
112~140
4000
3550~4500
160
140~180
5000
4500~5600
200
180~224
6300
5600~7100
250
224~280
8000
7100~9000
310
280~355
10000
9000~11200
400
355~450
12500
11200~14000
500
450~560
16000
14000~18000
630
560~710
4程序结果
此次程序验证信号为高斯白噪声,
结果如下所示:
1功率谱函数
2柱状图显示
5实验结论
由白噪声的频率特性可知,在相等带宽内能量是相同的。
用倍频程分析，则中心频率增加一倍，能量就增加一倍，即3dB；而用1/3倍频程进行分析时，中心频率增加1/3倍，能量就增加1/3倍，即1dB；由上图可知，该算法的正确性（由于低频部分的fft的线速很少，与理论有一些偏差，而中高频部分很好的反映了这一特性）。
与传统的方法相比，不但计算速度显著提高了，而且计算精度也提高了。
6实验体会
此次科研训练为期32个学时，由于对专业知识的匮乏，几经波折，终于完成。
通过这次科研训练，学到了不少东西：
首先，最大的体会是，自己的基础很薄弱，对很多东西都不了解，碰到了很多问题。
因此，在以后的学习中，需要加强对专业知识的学习。
其次，我对科研的过程有了一个大体的了解，并且对科学精神有了一个明确的认识，做学问要一丝不苟，不可敷衍了事，不可弄虚作假。
第三，学会了查找资料的方法，我现在明白了如果需要学习知识，除了图书馆里藏书和期刊之外，网络数据库也是一个不错的选择，比如中国期刊网，WebofScience，等；并且网络的作用也很大，比如一些专业的论坛，可以去论坛浏览，寻找自己需要的东西，还可以发帖问问题，等等。
第四，学会了一些软件的使用。
学会了一些Matlab编制一些简单的程序，使用Ansys对一个模型进行应力，应变分析和模态分析。
第五，对实验内容有一个深入的了解。
了解了声压级，功率谱分析和1/3倍频程
第六，学会了独立解决问题，在程序编制过程中遇到不少问题，出现了几次的反复和错误，通过查阅图书，上网搜集资料，逐一解决。
附：
Matlab程序带注释
%生成高斯白噪声
n=8192;
y=randn（1,n）;
y=y-mean（y）;
y=y/std（y）;
a=0;
b=sqrt
（1）;
y=a+b*y;
%采样频率
fs=51200;
%频率范围
f=[22.4,28,35.5,45,56,71,90,112,140,180,224,280,355,450,560,710,900,1120,1400,1800,2240,2800,3550,4500,5600,7100,9000,11200,14000,18000];
%中心频率
fm=[25,31.5,40,50,63,80,100,125,160,200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150,4000,5000,6300,8000,10000,12500,16000];
n3=length（fm）;
deltf=fs/n;
ff=0:
fs/n:
fs*（n/2-1）/n;
yf=do_fft（y）;
figure,plot（ff,yf）;
title（'高斯白噪声的功率谱'）;
ylabel（'功率/Pa'）;
xlabel（'频率/Hz'）;
yn=yf.*yf;
fori=1:
29
f1=f（i）;
f2=f（i+1）;
n1=floor（f1/deltf）;
n2=floor（f2/deltf）;
temp=0;
forj=n1:
1:
n2
temp1=yn（j）/20-5;
yn（j）=20*10^（temp1）;
temp=temp+yn（j）;
end
energy（i）=temp;
temp2=log10（energy（i）/（20*10^（-5）））;
ydb（i）=temp2*20/deltf;
end
figure
bar（ydb）;
title（'1/3倍频程分析结果'）;
ylabel（'分贝/dB'）;
xlabel（'频率/Hz'）;
参考文献
[]
【1】马大猷，噪声与振动手册【M】，机械工业出版社，2002。
9
【2】张沛上，噪声防治【M】，电子工业出版社，1988.10
【3】王济，胡晓，Matlab在振动信号处理中的应用，中国水利出版社，2006.1
【4】陈花玲，机械工程测试技术，机械工业出版社，2002.1
【5】王建春等，噪声测量中1/3倍频程与倍频程频谱的关系与计算【J】，噪声与振动控制，1996（5）：
39-41
【6】成峰等，基于Matlab的1/3倍频程FIR数字滤波器设计【J】，福州大学学报（自然科学版），2003，31
（2）：
160-165
【7】杨永红，一种基于MeasurementStudio计算动态信号1/3倍频程谱的方法【J】，北京工商大学学报（自然科学版），2003，21（3）：
37-39
【8】黄晶晶，雷勇，基于VC++的1/3倍频程设计与实现，电子测量技术，2006，29（6）135-136
【9】刘福强，张令弥，作动器传感器优化配置的研究进展【J】，力学进展，2001，3
（2）：
380
【10】KammerDC.Sensorplacementforon-orbitmodalidentificationandCorrelationoflarge
Structures[J],ControlandDynamics,1991,14
（2）:
251-259
【11】AmirFijanyandFarrokhVatan，ANewEfficientAlgorithmforAnalyzingand
OptimizingtheSystemofSensors，IEEEACpaper#1568,Version2,UpdatedDec,202005.
【12】M.Meo∗,G.Zumpano，Ontheoptimalsensorplacementtechniquesforabridgestructure
EngineeringStructures27（2005）1488–1497
【13】DanielC.Kammera,*,MichaelL.Tinker，Optimalplacementoftriaxialaccelerometersformodalvibrationtests，MechanicalSystemsandSignalProcessing18（2004）29–41
【14】AN-PANCHERNG，OPTIMALSENSORPLACEMENTFORMODALPARAMETERIDENTIFICATIONUSINGSIGNALSUBSPACECORRELATIONTECHNIQUES，MechanicalSystemsandSignalProcessing（2003）17
（2）,361d378
【15】FrankY.S.LinandP.L.Chiu,ANear-OptimalSensorPlacementAlgorithmto
AchieveCompleteCoverage/DiscriminationinSensorNetworks，IEEECOMMUNICATIONSLETTERS,VOL.9,NO.1,JANUARY2005
【16】EmanueleTrucco,Model-BasedPlanningofOptimalSensorPlacementsforInspection，IEEETRANSACTIONSONROBOTICSANDAUTOMATION,VOL.13,NO.2,APRIL1997
升级会员 