正交球面波源的近场声全息射线波函数波叠加法_石梓玉.pdf

第 36 卷第 2 期2023 年 4 月振 动 工 程 学 报Journal of Vibration EngineeringVol.36 No.2Apr.2023正交球面波源的近场声全息射线波函数波叠加法石梓玉1，2，向宇1，2，陆静1，2，王玉江1，2（1.广西科技大学广西汽车零部件与整车技术重点实验室，广西 柳州 545006；2.广西科技大学机械与汽车工程学院，广西 柳州 545006）摘要:在基于传统波叠加法的近场声全息技术中，多采用辐射球面波的单极子作为等效源，易导致传递矩阵病态，利用射线波函数替换球面波函数可有效改善传递矩阵病态性。然而，以往的射线波函数法采用格林函数的方向导数作为波函数，其解析表达式复杂，计算效率低。此外，以往方法的波函数指向设置对节点分布方式要求较高，限制了其应用的灵活性。针对上述问题，采用（n，0）阶的球面波源重新构造了一系列射线波函数，该射线波函数可利用球 Hankel函数和 Legendre多项式的递推形式方便地计算出其任意阶的表达式，大幅提高了效率。通过改进射线波函数的主指向设置，使其在实际使用中更加灵活，提出一种基于正交球面波源的射线波函数波叠加法。利用正四面体辐射体、两端带球帽的圆柱辐射体和简支矩形钢板声源 3个数值仿真，对比验证了传统方法和所提方法在声场重建中的效果。仿真结果表明：在声场细节信息较为丰富的高频下，即便采用正则化方法求解，传统波叠加法由于传递矩阵病态严重，在 3个仿真中均难以保证重建精度，误差为 20%左右；而建立的射线波函数法则有效降低了传递矩阵的条件数，改善了系统病态性，获得了更高的重建精度，误差在 5%10%之间，说明了所提方法更具优越性。关键词:近场声全息；波叠加法；正交球面波源；射线波函数中图分类号:O422.2；TB53 文献标志码:A 文章编号:1004-4523（2023）02-0534-11 DOI：10.16385/ki.issn.1004-4523.2023.02.025引言近场声全息技术（NAH）是一种有效的噪声源定位、识别和声场重建技术，经过近几十年的研究，在算法方面已相继提出了基于空间 Fourier 变换的NAH 方法1、基于 KirchhoffHelmholtz 边界积分方程的边界元方法24、源强模拟方法56、波叠加法7等。其中，波叠加法作为一种无奇异性、精度高且适用于任意形状结构的声场计算方法，自 1989年提出以来就已被广泛应用于各种声学问题的计算中811。其原理是在声源面内缩的一个封闭虚拟曲面上布置连续分布的等效源来表示声源向外辐射的声场，不仅克服了基于空间 Fourier 变换算法只能计算规则形状声源的缺点，又避免了边界元法所带来的复杂插值运算和奇异积分处理。但将波叠加法应用于声全息计算时，由于等效源面与全息测量面之间距离的影响，其离散后形成的传递矩阵通常是一个大条件数的病态矩阵，导致源强求解稳定性较差12。为了提高源强求解的稳定性，通常需采用正则化方法，目前最常用的正则化方法有截断奇异值方法（TSVD）、Tikhonov 正则化方法等13。但这些正则化方法本质上都是将传递矩阵中对测量误差非常敏感的小奇异值项进行截断或滤除，该过程必然会损失一部分声场细节信息。如果传递矩阵病态严重，那么正则化时就必须选取较大的截断点或正则化参数以过滤更多的小奇异值项，这会加剧声场细节信息的丢失。对于复杂程度不高的低频声场，小奇异值项对声场的贡献相对较小，在正则化后一般均可保证声场重建的分辨率。但对于声场细节信息丰富的高频声场，小奇异值项的滤除会对声场重建精度造成很大影响。因此，为保证波叠加法的重建精度，即便采用正则化方法，也应尽量改善传递矩阵的病态性。为改善传递矩阵的病态性，以往的研究大多侧重于优化等效源布置面的位置12，1418。但由于声源形状和性质、全息测量面形状、测点分布方式等的复杂性，等效源面最佳分布和位置的选择是一个非常复杂的问题，且至今尚无一个成熟有效的方法12，19。文献 2023 在对传统单极子波叠加法传递矩阵的病态性进行分析后发现，由于单极子向外辐射的波函数是以球面形式衰减的自由场格林函数，因此当等效源点之间或全息测量点之间的位置仅有微小改变时，格林函数的大小变化过于平缓，导致传递矩阵收稿日期:2021-08-27；修订日期:2021-11-04基金项目:国家自然科学基金资助项目（51775121）。第 2 期石梓玉，等：正交球面波源的近场声全息射线波函数波叠加法因不同行或列近似相等而病态。进而，文献 2223从优化波函数波阵面衰减速度的角度，提出了一种利用射线波函数替代传统球面波函数以改善传递矩阵病态性的射线波函数法，在改善传递矩阵病态性方面取得了一定效果。然而，该方法采用格林函数的方向导数作为射线波函数，这种类型射线波函数的解析表达式较为复杂，不仅计算效率低且难以计算到高阶导数。此外，该方法还将射线波函数的主指向设置为各等效源到其对应测点的方向，即要求等效源和测点无论是在数量还是分布方式上均必须一一对应，这无疑限制了射线波函数法在实际工程中的应用。为了提高射线波函数法的计算效率和灵活性，本文采用（n，0）阶正交球面波源构造了一种新型的射线波函数，该波函数无需求导运算，仅利用球Hankel函数和 Legrend多项式的递推形式即可方便地计算出其任意阶的表达式，不仅计算效率得到大幅提高，而且不会增加计算难度和复杂性。与此同时，本文通过分析射线波函数对传递矩阵病态性的改善机理，提出一种射线波函数的主指向设置方法，使等效源与测点的布置不再受到数量和分布均须一一对应的限制，大大提高了射线波函数法在应用中的灵活性。最后，通过 3 个不同声源的数值仿真算例对比验证了本文方法在声场重建中的有效性和准确性。1理论部分1.1传统单极子波叠加法的病态性及射线波函数法波叠加法的基本思想是：声源向外辐射的声场可由连续分布于其内部的等效源所辐射的声场叠加代替，该思想可用下式描述7：p(r)=(rE)(r，rE)d(rE)（1）式中p(r)为声源在空间r处辐射的声压；(rE)为位于rE处的等效源源强，下标表示其波函数为(r，rE)，根据等效源的类型不同，可为单极子等效源、偶极子等效源、单偶极子组合型等效源等；为等效源在声源内部的分布区域，SE为该区域的边界，如图 1（a）所示。为了便于计算，通常将等效源布置在虚拟边界SE上，并在积分离散时将每个单元内的等效源源强和波函数均视为常数，且配置在每个单元的中点，如图 1（b）所示。式（1）离散为如下形式56：p(r)=j=1Nq(rEj)(r，rEj)（2）式中q(rEj)为第j个离散等效源的源强；(r，rEj)为该等效源对应的波函数。由式（2）可知，只要确定了N个等效源的源强q(rE1)，q(rE2)，q(rEN)，即可计算出声源在任意场点r处辐射的声压。在基于波叠加法的近场声全息技术中，求解源强所需的方程组可以通过测量声源近场区域的声压或振速来建立。假设全息测量信息为声压，等效源为单极子，利用式（2）可得如下矩阵方程：pH=GHQG（3）式 中pH=p(rH1)p(rH2)p(rHM)T为M 1的测量声压列向量；GHij=G(rHi，rEj)为M N的全息测量数据与等效源强间的传递矩阵，G(rHi，rEj)为自 由 场 格 林 函 数，即 单 极 子 波 函 数；QG=qG(rE1)qG(rE2)qG(rEN)T为N 1的 单 极 子 等效源强列向量。求解方程（3），得源强向量QG为：QG=G+HpH（4）式中G+H表示矩阵GH的广义逆。为获得式（4）的最小二乘解，要求测点数M不小于等效源数目N，即M N。此外，在传统波叠加法中，一般采用单极子作为等效源，其波函数为自由图 1 波叠加法及等效源法的示意图Fig.1 The diagram for wave superposition method and equivalent source method535振 动 工 程 学 报第 36 卷场格林函数G(rHi，rEj)。文献 2023 在对由该函数所构成的传递矩阵GH的病态性进行分析后指出，由于格林函数G(rHi，rEj)是一个只与两点距离有关且以球面衰减的波函数，因此当不同测点或等效源点之间的位置仅有微小改变时，矩阵的不同行或列的元素将近似相等，导致传递矩阵GH因向量间的强线性相关性而病态。对于该病态问题，一般是借助正则化方法抑制测量误差的放大并以此稳定求解过程。常用的正则化方法包括截断奇异值方法（TSVD）、Tikhonov 正则化方法13等。但正则化方法本质上是通过类似滤波的方法将传递矩阵中对误差敏感的小奇异值项滤除。如果传递矩阵病态严重，在正则化时就需选取较大的截断点或正则化参数等以滤除更多的奇异值项，这将会加剧声场细节信息丢失，导致重建复杂声源或高频声场时精度下降。因而，为保证重建精度，应尽量改善波叠加法传递矩阵的病态性。为了改善上述问题，文献 2223 中提出了一种利用强指向性波函数替换球面形式波函数以提高声场重建稳定性的方法，并称之为射线波函数法。其基本原理是将传统波叠加法中单极子等效源辐射的球面波函数替换为满足 Helmholtz 方程和 Sommerfield辐射条件且主值指向等效源对应测点的射线波函数，如图 2所示。这样一来，等效源辐射的声波将仅在其对应测点处具有较大的声波激励，并生成较大的主对角元素，而在非对应测点处的声波激励则快速衰减，即生成较小的非对角元素，进而得到一个主对角元素占优的良态传递矩阵，以此提高声场重建的稳定性。1.2以往射线波函数法的缺陷及改进在对文献 2223 中提出的射线波函数法进行深入研究后发现，该方法虽然在改善传递矩阵病态性方面具有一定优势，但仍有以下两个缺陷：（1）波函数计算效率方面的缺陷。文献 2223均是利用格林函数的方向导数作为射线波函数，该类型波函数的高阶导数（超过 6阶后）解析表达式非常复杂，导致计算效率较低甚至无法计算。（2）射线波函数主指向设置方面的缺陷。文献2223 中均是将射线波函数的主指向设为等效源的对应测点方向，这样设置虽然可以形成主对角占优形态良好的传递矩阵，但要求等效源数量与测点数量相同，且在运算过程中须保证这两组节点的编号始终一一对应。而在实际应用中，为获得更高的计算精度，等效源最好均匀布置在与声源面共形的虚拟面上，同时为了便于制造和降低成本，全息测量面则通常为规则形状且测点规则分布，这必然难以保证全息测点与等效源点之间的一一对应。本文的主要内容则是针对以上两点提出如下改进办法：（1）对波函数计算效率方面的改进重新构造射线波函数。源模拟技术的研究表明，只要是满足 Helmholtz方程和 Sommerfeld辐射条件的解析函数(r，rE)均可作为等效源波函数56。因而射线波函数的选取并不局限于格林函数的导数，本文将利用 Helmholtz 方程在球坐标系下的基本解，即正交球面波源重新构造射线波函数。（2）对射线波函数主指向设置方面的改进。实际上，由 1.1 节中图 2 的分析可知，射线波函数能够降低传递矩阵线性相关性的主要原因在于它在非主指向的衰减速度远快于格林函数的球面波。因此，即便不采取等效源与主测点一一对应的设置方式，理论上仍应能显著降低传递矩阵因线性相关性过强导致的病态。在应用中，为了保证等效源辐射的声波在空间中分布均匀，只要各射线波函数的主指向均匀向外分散即可。例如，对于封闭的等效源面和全息面，可将射线波函数的主指向设置为从等效源面所包含空间的几何中心到各等效源连线的方向，如图 3（a）中的l1，l2，l3，lN所示；对于非封闭的等效源面和全息面，则可以设置为等效源面的外法向方向，如图 3（b）中的l1，l2，l3，lN所示。在后文的方法验证部分将通过仿真算例进行验证。1.3基于正交球面波源的射线波函数正交球面波源是由不同阶次的球谐函数和球Hankel函数构成，其具体表达式为24：pnm(r，)=hn(kr)Ymn(，)；n=0，1，2，3，；m=-nn（5）式中hn()为第一类或第二类 n 阶球 Hankel 函数，本文采用第二类 n 阶球 Hankel函数，即h(2)n()；Ymn()图 2 射线波和球面波在同一重建模型中的对比图Fig.2 Comparison of ray wave and spherical wave in the same reconstruction model536第 2 期石梓玉，等：正交球面波源的近场声全息射线波函数波叠加法为归一化的球谐函数；k为波数；(r，)表示场点在坐标系中的位置，如图 4所示。由式（5）可知，球面波源pnm(r，)的指向性取决于球谐函数项Ymn(，)。又由Ymn(，)的性质，当m=0时，其表达式中的角度变量将被消除，指向性仅由决定，此时该函数的指向形态必然关于 z轴回转对称且主值指向 z 轴。图 5 给出了m=0，n分别为 0，1，3，7，10，15 时球谐函数的指向形态图。由图 5 可见，当m=0时，球谐函数在 z 轴方向具有强指向性，且n越大，其指向性越强。因而，可取m=0时的球面波源pnm(r，)作为射线波函数。令式（5）中m=0，略去常系数后将其记为Dn(r，)：Dn(r，)=h(2)n(kr)Pn(cos )；n=0，1，2，3，（6）上式即为带有参数n的球面波源型射线波函数，其中Pn(cos )为n阶 Legendre 多项式。由于式（6）可利用球 Hankel函数和 Legendre多项式的递推形式方便地计算到任意阶，因此相较于格林函数导数型的射线波函数，不仅计算效率得到大幅提高，而且不会增加计算难度和复杂性。1.4射线波函数在波叠加法中的应用由 1.3 节的分析可知，将射线波函数应用于波叠加法时，需设置各等效源所辐射的射线波函数指向其对应主指向l1，l2，l3，lN，如图 3 中所示。下面假设第j个等效源在全局坐标系Oxyz中的位置为rEj，其对应的主指向为lj，如图 6 所示。由于射线波函数关于自身坐标的 z轴回转对称，因此，欲使其主瓣指向lj，可将位置rEj作为坐标原点、主指向lj为z轴作一局部坐标系Oxyz。此时，只要将射线波函 数Dn(r，)中 的 变 量和r替 换 为 局 部 坐 标 系Oxyz中的变量和r，即可使其指向z轴，即lj方向。利用式（6），可得rEj处单位强度等效源在场点r处辐射的声压为：Dn(r，)=h(2)n(kr)Pn(cos)（7）为方便起见，将上式中的局部坐标变量r和图 3 在采用封闭和非封闭全息采样面时射线波函数的主指向方向设置Fig.3 Setting of main direction of ray wave function in closed and nonclosed holographic surfaces图 4 球坐标系Fig.4 Spherical coordinate system图 5 球谐函数Y0n在三维空间中的指向形态Fig.5 Directional form of spherical harmonic function Y0n in threedimensional space537振 动 工 程 学 报第 36 卷用全局坐标变量rEj，r和lj表示如下：Dn(r，rEj，lj)=h(2)n(k|r-rEj|)Pnljr|lj|r（8）将上式中的场点位置r替换为全息面测点位置rHi，即可得到等效源与全息面间的传递矩阵DH，其元素为：DHij=Dn(rHi，rEj，lj)（9）进而式（3）和（4）的矩阵方程可分别改写为：pH=DHQD（10）QD=D+HpH（11）式中QD=qD(rE1)qD(rE2)qD(rEN)T为射线波函数所对应的等效源强向量。求解出等效源强后，再利用下式即可重建场点r处的声压：p(r)=j=1NqD(rEj)Dn(r，rEj，lj)（12）值得一提的是，在n=0和n=1时，式（6）分别对应单极子等效源和偶极子等效源25。此时，式（12）将变为传统的单层势和双层势波叠加法。2方法验证2.1仿真模型说明为验证本文所提方法在声场重建中的效果，该部分设计了 3 种不同类型声源的仿真模型，分别为正四面体辐射体、两端带球帽的圆柱辐射体和四边简支矩形钢板声源。需要指出的是，由于球面波源通常在单点多极法中被用于计算球形或近似球形的声源56，24，2627，而本文则是将其应用于波叠加法中，因此为了验证由球面波源构造的射线波函数在波叠加法中依然适用于任意形状声源的优势，上述3 种声源均为非球形声源。其中，正四面体辐射体主要考察等效源面和全息面均为封闭且共形，等效源数量与测点数量相同情况下的重建效果；两端带球帽的圆柱辐射体主要考察等效源面和全息面均为封闭，但仅近似共形，且等效源与测点并非一一对应情况下的重建效果；简支板声源则考察等效源面和全息面均为非封闭情况下对空间连续型结构声 源 的 重 建 效 果。3 个 仿 真 中 均 添 加 信 噪 比 为25 dB 的高斯白噪声，并使用 Tikhonov 正则化方法求解，正则化参数采用 L 曲线进行选择。同时，由前文所述，传递矩阵的病态性会导致正则化对声场细节信息的滤除增加，从而在重建高频复杂声场时精度下降，因此 3 个仿真均在 1000 Hz 以上的较高频下进行。2.2仿真 1：正四面体辐射体的声压重建如图 7（a）所示为一几何中心位于坐标原点O，底面的一边与y轴平行且边长为a=1 m的正四面体辐射体。假设该辐射体表面的声压由置于其中心O处的点声源S1和 z 轴正方向上与原点O距离为d=0.3 m的点声源S2叠加产生。为了增加声场的复杂度，将两个点声源设置为不同的多极子声源，其具体的表达式为：S1=h(2)0(krS1)Y00(S1，S1)（13）S2=h(2)4(krS2)Y04(S2，S2)（14）式中坐标变量rS1，S1，S1，rS2，S2，S2均为以点源所在位置为坐标原点、坐标轴方向与全局坐标系平行的局部坐标系中的位置标量。仿真中采用三角形单元均匀划分四面体表面，共计 340个节点。测量面设置为包裹辐射体且所有边长均为aH=1.2 m的共形四面体面，测点的数量和分布方式与声源表面节点相同，如图 7（b）所示。等效源布置在声源表面向内以 0.5 比例缩进的虚拟面上，等效源的数量和分布方式与声源表面节点相同。由于该模型中等效源与全息测点一一对应，因此设置射线波函数指向各等效源的对应测点。重建波数设置为k=35（1910.6 Hz）。图 8为正四面体声源的表面声压幅值分布。由图 8可见，即便已经采用了正则化方法进行求解，除了n=4阶波函数的重建结果与解析声压吻合较好以外，其余阶波函数的重建声压均产生了明显偏差。为了量化重建结果，图 9 给出了各阶波函数对应的重建声压相对误差曲线及传递矩阵条件数曲线，其中，相对误差由下式计算：Error=pn-p2p2 100%（15）图 6 指向性球面波源的局部坐标系Fig.6 Local coordinate system of directional spherical wave source538第 2 期石梓玉，等：正交球面波源的近场声全息射线波函数波叠加法式中pn表示第n阶波函数的重建声压向量；p表示解析声压向量。由图 9（a）可以发现，当波函数阶数n=0和n=1时（即 传 统 波 叠 加 法），重 建 误 差 较 大，超 过 了20%。而且当波函数的阶数在 04之间时，重建误差 基 本 呈 下 降 趋 势，并 在 4 阶 时 最 小，误 差 小 于10%。通过对比图 9（b）的条件数曲线可知，重建误差下降是由于传递矩阵条件数逐渐减少，因此正则化后损失的声场信息随之减少，重建结果更为精确和稳定。而当波函数的阶数大于 4 阶时，传递矩阵的条件数虽然仍保持下降趋势，但重建误差却反而逐渐上升。这是因为过高的阶数会使得射线波函数的指向性太强，并由此导致各等效源辐射的声波过于集中在其对应主指向的射线束内，而在非主指向上将因波函数衰减太快而产生较大误差，因此无法正确描述真实声场。由此可见，射线波函数的指向性并不是越强越好，而是应在合理的范围内选取。一个简单的方法是利用“辅助面法24”进行选取，即在声源面与全息面间设置若干辅助测量点，然后计算各阶射线波函数（通常只需要计算到前 1015阶）在辅助测点处的重建声压，最后将重建声压与该点的测量声压进行对比，取两者误差为最小时所对应的阶数即可。由于篇幅所限，本文不展开讨论波函数阶数的选取。实际上，由球面波源构造而成的射线波函数作为 Helmholtz 方程在球坐标系下的基本解，其阶数的选择需要综合考虑声源特性、NAH模型配置及噪声性质等28。2.3仿真 2：两端带球帽的圆柱辐射体的声压重建以坐标原点O为几何中心设置一两端带球帽的圆柱壳长条辐射体，其中圆柱部分长度为2a=1 m，圆柱和球帽的半径均为a=0.5 m，如图 10（a）所示。假设该辐射体表面声压由分别置于两个球帽球心处的点声源S1和S2叠加产生。与仿真 1 类似，将两个点声源设置为多极子，其中：S1=0.01 h(2)7(krS1)Y47(S1，S1)（16）图 7 正四面体辐射体的仿真模型示意图Fig.7 Simulation model of the radiator of regular tetrahedral图 8 正四面体辐射体表面的声压幅值分布Fig.8 Distribution of sound pressure amplitude on radiator of regular tetrahedral surface539振 动 工 程 学 报第 36 卷S2=h(2)2(krS2)Y02(S2，S2)（17）式中的坐标变量说明见仿真 1。仿真中仍采用三角形单元划分声源表面，共计 629个节点，如图 10（b）所示。测量面设置为包裹声源的长方体表面，其长宽高尺寸为2.1 m 1.1 m 1.1 m，在长方向等间隔分布 17 个测点，宽和高方向等间隔分布 9 个测点，测点数共计 642 个，如图 10（a）所示。等效源布置在声源表面向内以 0.5 比例缩进的虚拟面上，等效源的数量和分布方式与声源表面节点相同。注意，该模型中等效源与全息测点的数量和分布方式均不同，因此不存在一一对应关系，射线波函数的指向按照图 3（a）所示进行设置。重建波数设置为k=25（1364.8 Hz）。与仿真 1 类似，该仿真中分别给出了各阶波函数的声源表面重建声压幅值分布云图、声压相对误差曲线及传递矩阵条件数曲线，如图 11 和 12 所示。由图 12（a）可以发现，与仿真 1 不同，该重建模型中图 10 两端带球帽的圆柱辐射体仿真模型示意图Fig.10 Diagram of a radiator of cylinder with two spherical caps simulation model图 9 重建正四面体辐射体时的声压相对误差及传递矩阵条件数Fig.9 The relative error of sound pressure and the condition number of transfer matrix when reconstructing图 11 两端带球帽的圆柱辐射体表面声压幅值分布Fig.11 Distribution of sound pressure amplitude on radiator of cylinder with two spherical caps surface540第 2 期石梓玉，等：正交球面波源的近场声全息射线波函数波叠加法最小误差对应的波函数阶数为n=5。由图 12（b）的条件数曲线可以发现，虽然射线波函数并未指向等效源的对应测点，但传递矩阵的条件数仍然呈现良好的下降趋势，这表明本文在 1.3 节提出的射线波函数指向设置方法是有效的。2.4仿真 3：矩形简支板辐射声压的重建在波叠加法求解声场外问题的理论中，通常要求给定边界条件的全息面为一封闭曲面，但在实际工程中，重建由板、壳等结构振动产生的声场是一种常见的应用场景，此时一般采用非封闭的平面全息面测量声场信息。这一节的仿真将利用一个被无限大刚性障板围绕的简支板声源验证本文方法在全息面为非封闭平面时的声场重建效果。如图 13 所示为一左下角位于坐标原点、长方向和宽方向分别与y轴和x轴重合的长方形钢板。长宽厚尺寸设置为 1 m0.5 m0.003 m，杨氏模量为E=2.1 1011 Pa，泊松比为=0.23，密度为1=7.8 103 kg m3。边界条件设置为四边简支，并 在 板 的（0.25 m，0.75 m）位 置 处 施 加 一 幅 值 为1 N，波数为k=63（3439.2 Hz）的简谐激励力。全息测量面的大小与简支板相同并布置在其正上方0.1 m 处，测点在长方向和宽方向的扫描间隔均为0.02 m，共 1326 个测量点。等效源面位于简支板正下方 0.1 m 处，其分布方式与测点相同。射线波函数的指向按图 3（b）中的平面情形进行设置。仿真中重建简支板上方 0.05 m 处的声压，并与解析声压对比。简支板的解析声压是利用 Rayleigh 积分计算得到，详细的推导及表达式见文献 25。该仿真中除了给出与仿真 1 和仿真 2 相同的重建声压幅值分布云图、声压相对误差曲线及传递矩阵条件数曲线外，为了突出各阶波函数对简支板各区域声场的重建效果，还给出了重建面上的 相 对 误 差 分 布 云 图，上 述 结 果 如 图 1416 所示。从图 14 的声压幅值分布云图上看，各阶波函数的重建声压均与解析声压吻合得较好。但通过对比图 15 中的相对误差云图则可以发现，当波函数阶数为n=0和n=1时（即传统波叠加法），重建 误 差 云 图 中 有 大 面 积 误 差 15%的 深 红 色。而当波函数为n=28时，云图中深红色的大误差区域显著减少，并主要分布在解析声压接近于0 处和简支板边缘处，这是由于计算相对误差时分母接近 0 对误差产生的放大作用和声场信息泄露所导致的，其余大部分区域均为误差 0%5%的深蓝色和浅蓝色。由图 16（a）可以发现，该重建模型的最小误差对应的波函数阶数为n=4，最小误差在 5%左右，且阶数在n=28的范围内重建误差均相差不大，都在 5%10%之间。从图16（b）中可以看到，传递矩阵的条件数曲线仍然呈下降趋势。该仿真表明，对于非封闭的重建模型，本文方法仍具有一定优势。图 12 重建两端带球帽的圆柱辐射体时的声压相对误差及传递矩阵条件数Fig.12 The relative error of sound pressure and the condition number of transfer matrix when reconstructing图 13 四边简支矩形钢板的尺寸及等效源的分布Fig.13 Size of simply supported rectangular steel plate with four edges and distribution of equivalent source541振 动 工 程 学 报第 36 卷图 14 偏心激励下的四边简支矩形钢板在重建面上的声压幅值分布Fig.14 Amplitude distribution of sound pressure of simply supported rectangular steel plate with four edges on the reconstruction surface under eccentric excitation图 15 偏心激励下的四边简支矩形钢板在重建面上的重建声压相对误差分布Fig.15 Relative error distribution of reconstructed sound pressure of simply supported rectangular steel plate with four edges on the reconstruction surface under eccentric excitation图 16 重建偏心激励的简支矩形钢板声源时的声压相对误差及传递矩阵条件数Fig.16 The relative error of sound pressure and the condition number of transfer matrix when reconstructing the simply supported rectangular steel plate sound source542第 2 期石梓玉，等：正交球面波源的近场声全息射线波函数波叠加法3结 论针对波叠加法传递矩阵的病态问题，在射线波函数法的基础上，提出了一种由（n，0）阶正交球面波源构造而成的射线波函数，并改进了射线波函数的主指向设置。文中对所提方法进行了详细的推导和阐述，并通过数值仿真进行了验证，主要结论有：（1）由（n，0）阶正交球面波源构成的射线波函数，其指向性随着阶数 n的增大逐渐增强，并可利用Hankel 函数和 Legendre 多项式的递推形式方便地获得任意阶表达式，相较于以往格林函数导数型的射线波函数，其计算难度和复杂程度大大降低；（2）在射线波函数法中，射线波函数的主指向无需严格指向等效源的对应测点。只需各射线波函数的主指向在空间中均匀向外分散，即可有效地减少传递矩阵条件数，并保证声场重建精度，大大提高了射线波函数法在应用中的灵活性；（3）在 1000 Hz以上的含噪声高频声场中，由于声场细节信息丰富，即便采用 Tikhonov正则化方法求解，传统波叠加法仍无法保证重建精度，误差为20%左右，而射线波函数法则能获得更高的重建精度，误差在 5%10%之间。参考文献:1Williams E G，Maynard J D，Skudrzyk E.Sound source reconstructions using a microphone array J.The Journal of the Acoustical Society of America，1980，68（1）：340344.2 Veronesi W A，Maynard J D.Digital holographic reconstruction of sources with arbitrarily shaped surfaces J.The Journal of the Acoustical Society of America，1989，85（2）：588-598.3BAI M.Application of BEMbased acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometriesJ.The Journal of the Acoustical Society of America，1992，92（1）：533-549.4Cunefare K A，Koopmann G，Brod K.A boundary element method for acoustic radiation valid for all wavenumbers J.The Journal of the Acoustical Society of America，1989，85（1）：3948.5Ochmann M.The source simulation technique for acoustic radiation problemsJ.Acta Acustica United with Acustica，1995，81（6）：512527.6Ochmann M.The fullfield equations for acoustic radiation and scattering J.The Journal of the Acoustical Society of America，1999，105（5）：25742584.7Koopmann G H，Song L M，Fahnline J B.A method for computing acoustic fields based on the principle of wave superposition J.The Journal of the Acoustical Society of America，1989，86（6）：24332438.8Du B K，Zeng X Y，Vorlnder M.Multizone sound field reproduction based on equivalent source methodJ.Acoustics Australia，2021，49（2）：317-329.9Chen J，Liu C，Zhang X Z，et al.An approach for indoor prediction of the passby noise of a vehicle based on the timedomain equivalent source methodJ.Mechanical Systems and Signal Processing，2021，146：107037.10 刘延善，曾向阳，王海涛.封闭空间声场重构的多层等效源法 J.声学学报，2020，45（3）：367376.Liu Yanshan，Zeng Xiangyang，Wang Haitao.3D sound field reconstruction for the enclosed cavity using the multilayer equivalent sources methodJ.Acta Acustica，2020，45（3）：367376.11 Valdivia N P.Advanced equivalent source methodologies for nearfield acoustic holographyJ.Journal of Sound and Vibration，2019，438：6682.12 BAI M R，CHEN C C，Lin J H.On optimal retreat distance of equivalent source methodbased nearfield acoustical holography J.The Journal of the Acoustical Society of America，2011，129（3）：14071416.13 Williams E G.Regularization methods for nearfield acoustical holography J.The Journal of the Acoustical Society of America，2001，110（4）：19761988.14 Valdivia N P，Williams E G.Study of the comparison of the methods of equivalent sources and boundary element methods for nearfield acoustic holography J.The Journal of the Acoustical Society of America，2006，120（6）：36943705.15 Gounot Y J R，Musafir R E.On appropriate equivalent monopole sets for rigid body scattering problems J.The Journal of the Acoustical Society of America，2007，122（6）：31953205.16 Gounot Y J R，Mu