基于巴特沃斯滤波器的近场声全息滤波方法.pdf

第45卷第2 4期2023年1 2 月舰船科学技术SHIP SCIENCEANDTECHNOLOGYVol.45,No.24Dec.,2023基于巴特沃斯滤波器的近场声全息滤波方法卢锦芳1，楼京俊，刘树勇1，毛为民3（1.海军工程大学动力工程学院，湖北武汉430 0 33；2.海军工程大学舰船与海洋学院，湖北武汉430033；3.中国人民解放军91 0 0 7 部队，浙江舟山31 6 0 0 0）摘要：针对近场声全息研究中声场重构时，指数滤波器无法兼顾通带衰减率和高频带衰减率的局限性，采用巴特沃斯滤波器作为滤波函数以消除噪声干扰，对2 个同向振动的脉动球源的重构声场进行仿真研究，给出巴特沃斯滤波器随参数变化的特性曲线。基于重构误差最小原则，给出该滤波器的最优参数，并与指数滤波器效果进行对比。结果表明，巴特沃斯滤波器可在声场重构过程中进行有效滤波，且在高频段，巴特沃斯滤波器重建误差明显小于指数滤波器重建误差，滤波效果更好，重建精度更高。关键词：近场声全息；波数域滤波；巴特沃斯滤波器；指数滤波器中图分类号：TB532文章编号：1 6 7 2-7 6 49(2 0 2 3)2 4-0 1 6 0-0 6Near-field acoustic holography filtering method based on Butterworth filter(1.College of Power Engineering,Naval University of Engineering,Wuhan 430033,China;2.College of Naval Architecture and Ocean Engineering,Naval University of Engineering,Wuhan 430033,China;3.No.91007 Unit of PLA Liberation Army,Zhoushan 316000,China)Abstract:In view of the limitation that the exponential filter cannot take into account both passband attenuation rateand high frequency attenuation rate when the sound field is reconstructed in the study of near-field acoustic holography,But-terworth filter is adopted as the filter function to eliminate the noise.The reconstructed sound field of two pulsating sphericalsources with the same direction of vibration is simulated.The characteristic curve of Butterworth filter changing with theparameters is given.Based on the principle of minimum reconstruction error,the optimal parameters of the filter are given,and the filtering effect is compared with that of the exponential filter.The results show that Butterworth filter can effectivelyfilter in the process of sound field reconstruction,and the reconstruction error of Butterworth filter is significantly smallerthan that of exponential filter in high frequency band and its filtering effect is better and the reconstruction accuracy is higher.Key words:near-field acoustic holography;k-space filtering;Butterworth filter;exponential filter0引言近场声全息技术（NearfieldAcousticHolography，NAH）由Maynard等在2 0 世纪8 0 年代提出 1-2 1，是声学研究中一种具有高分辨率重建精度的噪声源定位和声场可视化方法，在声场重构以及振动噪声控制等领域有着广泛应用。近场声全息通过多个传感器组成的测量面在被测声源物体表面的近场区域获取全息测量数据，建立重建面与测量面之间的空间变换关系，进收稿日期：2 0 2 2-1 1-30基金项目：湖北省自然科学基金面上项目（2 0 2 2 CFB405）；国家自然科学基金青年科学基金资助项目（52 2 0 1 38 9）作者简介：卢锦芳（1 999），女，博士研究生，研究方向为近场声全息。文献标识码：ALU Jin-fang,LOU Jing-jun,LIU Shu-yong,MAO Wei-mindoi:10.3404/j.issn.1672-7649.2023.24.029而重建三维空间的声压场、声强场、质点振速场或者预测远场指向性，从而为声源识别定位、有效控制结构振动噪声提供技术参考 3。在声场逆向重构过程中，位于高波数区的逝波成分将与噪声误差一起被逆传递因子按指数规律放大，导致重建结果产生很大误差，严重影响声场重建的准确性 4。针对这一问题，在重建过程中，必须对高波数区噪声误差进行滤波处理，以确保重建精度，为此国内外学者进行了广泛研究。Veronesi等 5 提出一第45卷种在波数域加低通滤波窗的方法来抑制高波数成分波对噪声的放大效应，该方法采取指数滤波器对波数域噪声信号进行处理，其算法简单，可操作性强，但该方法没有考虑声波频率、信噪比以及测量面与源面间距离等因素的影响，使用中常会出现因截止波数取值不当以致滤波效果不理想的情况。赵晨等 6 对指数滤波器中窗函数参数的合理设置进行分析，针对截止波数和窗函数陡度系数2 个关键参数，总结了在不同声源频率下取得最优重建效果的截止波数和陡度系数的最优取值，对窗函数的设计优化具有重要意义。张德俊7 在系统存在噪声情况下，对指数滤波器的截止波数进行修正，同时还考虑到指数滤波器仅在测量有限距离时有效，随着测量距离增大，滤波窗收效甚微，提出一种带约束条件的最小二乘法滤波函数。罗禹贡等 8 总结了候逝波随距离变化的规律并设计了相应的空间频域滤波器。由于指数滤波器仅对被逆传递函数按指数规律放大的噪声误差有效，对其他区域噪声误差过滤效果不佳，何其伟等 9 提出一种波数域循环加窗逼近法，通过循环加窗筛选最接近理论声压幅值的数据作为重建数据，从而确保声场重构的准确性。莫登沅等 1 0-1 2 将空间域滤波和波数域滤波相结合，提出一种空间域即波数域联合迭代滤波窗，该滤波窗可对重建面边缘误差进行有效消除，并且经过多次选代，明显提高了声场重建精度。在指数滤波器实际使用中发现，窗函数陡度系数的取值变化决定着窗函数的衰减速率，陡度系数取值越小，指数窗函数衰减越快，相较于截止波数的取值变化，陡度系数取值变化对反演效果影响较小，故其对指数滤波器衰减程度的调整也有限。此外，指数滤波器的陡度系数会同时影响通带内和高频带内滤波曲线的衰减速度，且不能时时兼顾两者之间的平衡，甚至会导致有价值的候逝波被剔除，影响滤波效果，进而影响声场重建精度。本文针对上述问题，选取2 个脉动球源作为研究对象，采用巴特沃斯数字滤波器（Butterworthfilter）作为滤波函数，对脉动球源声场反演过程进行数值仿真研究。基于重构误差最小原则，得到了使声场重建误差最小的巴特沃斯滤波器最优参数，并比较了该参数条件下，巴特沃斯滤波器与常用指数滤波器的效果。1基于Fourier空间声场变换的近场声全息在均匀理想流体介质中，小振幅声波的三维波动方程为：卢锦芳，等：基于巴特沃斯滤波器的近场声全息滤波方法卷积转化为波数域乘积可得：P(kx,ky,ZH)=P(kx,ky,zs)GD(kx,ky,zH-zs)。式中：P(kx,ky,zH)为测量面声压p(xH,yH,zH)的空间Fourier变换；P(kx,ky,zs)为重建面声压p(xs,ys,zs)的空间Fourier变换；Gp(k,ky,zH-zs)为Green函数gD(x,y,z)的空间Fourier变换。其解析表达式为：Gp(kr,ky,ZH-zs)=ei(cH-s)k.。式中：zZH为测量面z轴坐标，zs为声源面z轴坐标，kx、k y 分别为x轴和y轴方向波数，k,取值存在以下2种情况（+k=k2被称之为辐射圆，辐射圆区分了声波2 种不同的传播模式）1 3。当+k0z0的空间为自由声场时，由Green公式可得任意平面上的声压为：p(XH,H,ZH)=JJ,p(xs,ys,zs)gD.(XH-Xs,yH-ys,ZH-zs)dxsdyso式中：S为用于积分的无穷大平面；P(xH,YH,ZH)为测量面声压；p(xs,ys,zs)为重建面声压；gD(xH-xs,yH-y，zH-zs)为满足Dirichlet边界条件的格林函数。其表达式为：8D(a,2)=k2时：k,=i/(k2+k)-k2。当kz=-i/(k+磅)-k2时，其对应的解没有物理意义，在无穷远处趋向于无穷大，故此时对k,取值为kz=i/(+)-k2。在此情况下，k为纯虚数，GD(ks,k,ZH-zs)为波数空间的高频成分，在声波传播过程中，其幅值在传播方向随距离变化按指数规律衰减，称其为候逝波，为了获取到候逝波中包含的重要信息，需在声源近场对其进行测量。当已知测量面上声压分布p(xH,YH,zH)时，对式（5）进行空间二维Fourier逆变换，可得重建面上声压为：p(xs,ys,zs)=F-IP(kx,ky,zH)GB(k,ky,ZH-zs)。(9)式中，G(kx,ky,zH-zs)为声压重建的逆向传递因子。基于空间声场变换的平面近场声全息逆向重构流程，如图1 所示。测量面数据分析测量面声压二维空间测量面角谱数据采集傅里叶通P(k.k,zn)变换G(k,k,ZH-z)重建面声压二维空间加窗后重建面角谱加窗滤波p(a,y,z)博里叶过P(k,k,zm)f(c,M变换图1 基于空间声场变换的平面近场声全息逆向重构流程图Fig.1 Reverse reconstruction flow chart of planar near-field acous-tical holograph based on spatial sound field transformation2滤波器在实际声场反演过程中，位于高波数区的候逝波成分将与噪声误差一起被逆传递函数按指数规律放大，若不对噪声信号进行低通滤波，将导致重建结果产生很大误差，严重影响声场重建的准确性，因此必须在声场逆向重构过程中进行滤波处理2.1指数滤波器指数滤波器是声场重构过程中常用的滤波方法，其算法简单、可操作性强、应用广泛。截止波数和陡度系数是该方法的2 个关键系数，截止波数对参与声场重构过程的测量面声压角谱范围做出了限定，其取值受声波频率、信噪比以及测量面与源面间距离等因素的影响。截止波数取值过大，则达不到滤波效果，重建误差大；截止波数取值偏小，将导致过滤波，影响重建精度。陡度系数决定着窗函数的衰减速率，其取值越小，窗函数在截止波数处截止得越陡峭，且陡度系数会同时影响通带和高频带内滤波曲线的衰减速度，有时无法兼顾两者之间的平衡，甚至会导致有价舰船科学技术值的逝波被剔除，影响滤波效果，进而影响声场重建准确性。(8)指数滤波器计算公式为：exp(kr/kc-1)/a,k,kc,2V(Gex(1-/k)/al,k,ke.1式中：kr=V(k+)；k e 为截止波数；为窗函数陡度系数。2.2巴特沃斯滤波器巴特沃斯滤波器被称为最大平坦滤波器，其在通带的频率响应曲线符合最大限度平坦特性，是连续衰减的滤波器，该滤波器能够保持增益的平坦特性，尤为适合设计低通滤波器进行滤波处理，消除高频噪声干扰。巴特沃斯滤波器过渡带的陡峭程度与阶次成正比，随着阶次的增加，幅频特性曲线越陡峭，巴特沃斯滤波器更接近于理想低通状态，抑制噪声的能力变强。随着截止频率增加，系统响应速度变快，但系统稳定性也随之下降 4。巴特沃斯滤波器计算公式为：1重建面角谱f(C,M)=-P(k.,k,z)1+C2式中：C=(100.1p-1)为巴特沃斯滤波器截止频率；p为通带最大允许衰减系数；N为巴特沃斯滤波器阶数；;k=/&+)。3仿真研究与参数分析3.1两脉动球源声场模型如图2 所示，两脉动球源的中心坐标分别为C1=(0.175,0,0)和C2=(-0.175,0,0)，半径ro为0.0 1 m,两脉动球同相振动，表面振动速度Ua=1m/s，频率f-1000Hz。重建面zs=0.01m，测量面zH=0.05m，重建距离为0.0 4m，采样间距为0.0 5m，测点数量为2 52 5，在测量面上等间距分布。脉动球源辐射声压表达式：r式中：r为单个声源到测点距离；二绝对值为声压振幅。A=-i.kpcoQ0。4元式中：p为空气密度；co为空气中声速；Q=4元rUa第45卷(10)(11)2N。krAr(12)(13)第45卷1在频率为1 0 0 0 Hz时，测量面测得声压分布图如图3所示。在理论计算的基础上，添加高斯白噪声，模拟信噪比为2 0 dB时，测量面测得声压如图4所示。0.6160.414120.210u/00.20.4-0.6-0.5图3测量面声压分布图Fig.3 Sound pressure distribution on the measuring surface0.60.40.2u/00.2-0.4卢锦芳，等：基于巴特沃斯滤波器的近场声全息滤波方法X个重建面CZs图2 声场模型Fig.2Sound field model0 x/m163图5所示。可知，随着截止波数的增加，重建面误差全息面急剧下降，而在该阶段相较于截止波数的取值变化，陡度系数取值变化对重建误差影响较小。当指数滤波ZH器截止波数到达kc=2.5k附近时，滤波效果最好，重建误差最小，此时重建效果达到最优。当截止波数继续增大，参与重建过程的测量面角谱范围进一步扩大，被噪声误差淹没的高波数区逝波未被滤掉，滤波效果变差。当kc=2.5k，=0.2 时，重建误差达到最小，选取此时的指数滤波器参数为最优参数，该参数下指数滤波器如图6 所示。10090Pa80%/善重706050403082061040020.5Pa161412108642=0.01a-0.a-0:15252052.53.0510图5指数滤波器误差曲线图Fig.5Error curve of exponential filterk/k=2.5,a=0.21.00.80.60.40.2050k,/rad/m0-50k.=2.5ka-0.2L-14.07%1520k/k0-5025k./rad/m3050-0.6-0.5图4加噪后测量面声压分布图Fig.4 Sound pressure distribution image with noise on the measur-ing surface3.2指数滤波器参数设计为了能够定量描述滤波器滤波效果，定义重建误差为：E=Ps-pthe2100%。/pheP式中：ps为重建声压；pthe为理论声压。声源频率为1 0 0 0 Hz时，指数滤波器误差曲线如图6 最优参数下指数窗函数示意图00.5x/mFig.6 Diagram of exponential window function with3.3巴特沃斯滤波器参数设计声源频率为1 0 0 0 Hz时，巴特沃斯滤波器误差曲线如图7 所示。可以看出，在阶次N3时，随着通带允许最大衰减系数p的取值增大，滤波器截止频率增大(C=(100.1ap-1)，重建误差不断降低，误差曲线存在最低点，此时重建误差最小，重建效果达到最优。(14)随着阶次继续增加，误差曲线不存在拐点，随通带允许最大衰减系数p的取值单调递增，p越大，截止频率越大，重建过程中，过高波数的逝波未被滤掉，将连同噪声误差被逆传递函数按指数规律放大，造成optimal parameters:164 极大的重建误差。当N=2，p=0.6时，重建误差达到最小，选取此时的巴特沃斯滤波器参数为最优参数，该参数下巴特沃斯滤波器如图8 所示。100-0-N-190上-4N-280-N=370*N=4%/里重60上一N=5504030201000图7巴特沃斯滤波器误差曲线图Fig.7Error curve of Butterworth filter,=0.6,阶次N=21.00.80.60.40.2050k/rad/m0-50图8 最优参数下巴特沃斯窗函数示意图Fig.8Diagram of Butterworth window function withoptimal parameters3.4巴特沃斯滤波器与指数滤波器滤波效果分析当声源频率为1 0 0 0 Hz时，分析图5指数滤波器误差曲线可知，指数滤波器选择最优参数：截止波数kc=2.5k，窗函数陡度系数=0.2时，重建误差最小值为1 4.0 7%；分析图7 巴特沃斯滤波器误差曲线可知，巴特沃斯滤波器在选择最优参数阶次N=2，通带最大允许衰减系数p=0.6时，重建误差最小值为1 3.31%。经指数窗函数滤波后重建结果如图9 所示，经巴特沃斯窗函数滤波后重建结果如图1 0 所示。通过与理论值对比发现，这2 个滤波器均较好地实现声场重构，重构误差相对较小，其中巴特沃斯滤波器重建误差较指数滤波器更小。这一结果说明在声场模型相同时，选择最优参数的巴特沃斯滤波器滤波效果比选择最优参数的指数滤波器效果更好，巴特沃斯滤波器在声场重构过程中可进行有效滤波。不同频率条件下，巴特沃斯滤波器与指数滤波器误差曲线如图1 1 所示。重建距离为0.0 4m时，巴特沃舰船科学技术重建面声压理论值Pa1.530501.0200100.5-500-50050重建面声压计算值Pa重建面角谱计算值6001.52550201.0150100.55-504444/444/44444440.51.01.5a,-0.6-13.31%0.51.0通带允许最大衰减/dB0-50k/rad/m第45卷重建面角谱理论值500400300200100500400300200100-50050图9指数窗函数滤波后重建结果1.52.02.53.050Fig.9 Reconstructed result filtered by exponential重建面声压理论值Pa1.530501.0200100.5-5000.51.01.5重建面声压计算值Pa1.52550201.01500.5105-500.551.01.5图1 0 巴特沃斯窗函数滤波后重建结果Fig.10 The reconstructed result filtered by Butterworthwindowfunction斯滤波器全频段平均误差为1 4.3%，指数滤波器为14.5%，全频段滤波效果相当。从图中可看出，在高频段，巴特沃斯滤波器重建误差明显小于指数滤波器重建误差，滤波效果更好。这是由于指数滤波器的陡度系数会同时影响通带内和高频带内滤波曲线的衰减率，且不能兼顾两者之间的平衡，导致有价值的候逝波被滤出，从而影响滤波精度。后移测量面至zH=0.1m，重建距离为0.0 9 m时，不同频率条件下，巴特沃斯滤波器与指数滤波器误差曲线如图1 2 所示。此时，巴特沃斯滤波器全频段平均误差为30.6%，指数滤波器为31.6%，巴特沃斯滤波器滤波效果更好。从图中可看出，改变重建距离后，在高频段巴特沃斯滤波器重建误差仍明显小于指数滤波器重建误差，滤波效果更好。4结语针对指数滤波器无法兼顾通带衰减率和高频带的衰减率，可能导致有价值的候逝波被剔除，影响声场window function重建面角谱理论值500400300200100-500重建面角谱计算值600500400300200100-5005050第45卷16.516.015.5%15.014.514.013.513.012.512.00图1 1 重建距离0.0 4m时，巴特沃斯滤波器与指数滤波器误Fig.11The comparison of error curves of Butterworth filter andexponential filter with reconstruction distance 0.04 m403836343230282624220图1 2 重建距离0.0 9m时，巴特沃斯滤波器与指数滤波器误Fig.12The comparison of error curves of Butterworth filter andexponential filter with reconstruction distance 0.09 m重建精度问题，本文提出采用巴特沃斯数字滤波器作为滤波函数的近场声全息滤波方法，选取2 个同向振动的脉动球源作为研究对象，对两脉动球源声场反演过程进行数值仿真。通过仿真分析得到了使声场重建误差最小的巴特沃斯数字滤波器最优参数，并使用最优参数下该滤波器与最优参数下指数滤波器的滤波效果进行了对比。结果表明，巴特沃斯滤波器可替代指数滤波器在声场重构过程中，进行滤波处理，且在高频段巴特沃斯滤波器重建误差明显小于指数滤波器重建误差，滤波效果更好，重建精度更高。参考文献：1 MAYNARD J D,WILLIAMS E G,LEE Y.Nearfield acousticholography:I.Theory of generalized holography and the devel-opment of NAHJ.The Journal of the Acoustical Society ofAmerica,1985,78(4):1395413.2 WILLIAMS E G,MAYNARD J D,SKUDRZYK E.Soundsource reconstructions using a microphone arrayJ.The Jour-卢锦芳，等：基于巴特沃斯滤波器的近场声全息滤波方法一一巴特沃斯滤波器一一指数滤波器5001000Frequency/Hz差曲线一一巴特沃斯滤波器一一指数滤波器5001000Frequency/Hz差曲线165nal of the Acoustical Society of America,1980,68(1):340-344.3 WILLIAMS E G,FOURIER A.Sound radiation and nearfieldacoustical holographyJ.Academic Press,1999:590.【4于飞，陈剑，周广林，等.噪声源识别的近场声全息方法和数值仿真分析 .振动工程学报，2 0 0 3，1 6(3)：39-343.YU F,CHEN J,ZHOU G L,et al.A study on noise sourcesidentification using near field acoustic holography method andits numerical simulation analysisJ.Journal of Vibration Engi-1500200015002.000neering,2003,16(3):39-343.5 VERONESI W A,MAYNARD JD.Near-field acoustic holog-raphy(NAH)II.holographic reconstruction algorithms andcomputer implementationJj.The Journal of the Acoustical So-ciety of America,1987,81(5):;1307-22.【6】赵晨，冯丹平，杨明明，等.平面近场声全息中指数滤波器窗函数设计优化 .声学技术,2 0 2 1,40(6):7 2 3-7 2 7.ZHAO C,FENG D P,YANG M M,et al.Optimization of win-dow function design of exponential filter in planar near-fieldacoustic holography.Technical Acoustics,2021,40(6):723-727.【7】张德俊.近场声全息对振动体及其辐射场的成像 .物理学进展,1 996,1 6(3/4):6 1 4-6 2 3.ZHANG D J.Imaging for vibrating mode and radiation field ofvibrating object using NANJ.Progress in Physics,1996,16(3/4):614623.8罗禹贡，郑四发，杨殿阁，等.基于逝波衰减特性的空间频域滤波器研究 J.声学技术,2 0 0 4(1):54-56+6 6.LUO Y G,ZHENG S F,YANG D G,et al.K-space filterbased on decaying characteristics of evanescent wavesJ.Technical Acoustics,2004(1):54-56+66.【9何其伟，许国良，陈志敏，等.波数域循环加窗逼近法滤波方法研究 .现代电子技术,2 0 1 5,38(3):6 4-6 7.HE Q W,XU G L,CHEN Z M,et al.K-space filtering bycyclically windowing approximation methodJ.Modern Elec-tronics Technique,2015,38(3):64-67.10莫登沅,周其斗，谷高全,等.声全息重构中的滤波窗优化设计 .海军工程大学报,2 0 1 7,2 9(3):1 0 9-1 1 2.MO D Y,ZHOU Q D,GU G Q,et al.Design of filtering win-dow for acoustical holographyJ.Journal of Naval Universityof Engineering,2017,29(3):109-112.11】黄波.巴特沃斯数字滤波器的设计与仿真实现 .河南科技,2 0 2 1,40(36):1 0-1 2.HUANG B.Design and simulation of butterworth digitalfilterJJ.Henan Science and Technology,2021,40(36):10-12.12肖锦龙.巴特沃斯低通滤波器在图像增强中的应用研究 智能计算机与应用，2 0 2 1,1 1(5)：1 7 9-1 8 2.XIAO J L.Application of Butterworth low pass filter in imageenhancementJ.Intelligent Computer and Application,2021,11(5):179-182.13】陈心昭，毕传兴，等.近场声全息技术及其应用 M.北京：科学出版社,2 0 1 3.14】刘凯悦.复杂系统的解耦内模与事件触发故障补偿控制 D.北京：北京化工大学2 0 2 1.