正交球面波源的近场声全息射线波函数波叠加法

石梓玉 向宇 陆静 王玉江

摘要 在基于传统波叠加法的近场声全息技术中，多采用辐射球面波的单极子作为等效源，易导致传递矩阵病态，利用射线波函数替换球面波函数可有效改善传递矩阵病态性。然而，以往的射线波函数法采用格林函数的方向导数作为波函数，其解析表达式复杂，计算效率低。此外，以往方法的波函数指向设置对节点分布方式要求较高，限制了其应用的灵活性。针对上述问题，采用（n， 0）阶的球面波源重新构造了一系列射线波函数，该射线波函数可利用球Hankel函数和Legendre多项式的递推形式方便地计算出其任意阶的表达式，大幅提高了效率。通过改进射线波函数的主指向设置，使其在实际使用中更加灵活，提出一种基于正交球面波源的射线波函数波叠加法。利用正四面体辐射体、两端带球帽的圆柱辐射体和简支矩形钢板声源3个数值仿真，对比验证了传统方法和所提方法在声场重建中的效果。仿真结果表明：在声场细节信息较为丰富的高频下，即便采用正则化方法求解，传统波叠加法由于传递矩阵病态严重，在3个仿真中均难以保证重建精度，误差为20%左右；而建立的射线波函数法则有效降低了传递矩阵的条件数，改善了系统病态性，获得了更高的重建精度，误差在5%～10%之间，说明了所提方法更具优越性。

关键词 近场声全息; 波叠加法; 正交球面波源; 射线波函数

引 言

近场声全息技术（NAH）是一种有效的噪声源定位、识别和声场重建技术，经过近几十年的研究，在算法方面已相继提出了基于空间Fourier变换的NAH方法［1］、基于Kirchhoff?Helmholtz边界积分方程的边界元方法［2?4］、源强模拟方法［5?6］、波叠加法［7］等。其中，波叠加法作为一种无奇异性、精度高且适用于任意形状结构的声场计算方法，自1989年提出以来就已被广泛应用于各种声学问题的计算中［8?11］。其原理是在声源面内缩的一个封闭虚拟曲面上布置连续分布的等效源来表示声源向外辐射的声场，不仅克服了基于空间Fourier变换算法只能计算规则形状声源的缺点，又避免了边界元法所带来的复杂插值运算和奇异积分处理。但将波叠加法应用于声全息计算时，由于等效源面与全息测量面之间距离的影响，其离散后形成的传递矩阵通常是一个大条件数的病态矩阵，导致源强求解稳定性较差［12］。

为了提高源强求解的稳定性，通常需采用正则化方法，目前最常用的正则化方法有截断奇异值方法（TSVD）、Tikhonov正则化方法等［13］。但这些正则化方法本质上都是将传递矩阵中对测量误差非常敏感的小奇异值项进行截断或滤除，该过程必然会损失一部分声场细节信息。如果传递矩阵病态严重，那么正则化时就必须选取较大的截断点或正则化参数以过滤更多的小奇异值项，这会加剧声场细节信息的丢失。对于复杂程度不高的低频声场，小奇异值项对声场的贡献相对较小，在正则化后一般均可保证声场重建的分辨率。但对于声场细节信息丰富的高频声场，小奇异值项的滤除会对声场重建精度造成很大影响。因此，为保证波叠加法的重建精度，即便采用正则化方法，也应尽量改善传递矩阵的病态性。

為改善传递矩阵的病态性，以往的研究大多侧重于优化等效源布置面的位置［12，14?18］。但由于声源形状和性质、全息测量面形状、测点分布方式等的复杂性，等效源面最佳分布和位置的选择是一个非常复杂的问题，且至今尚无一个成熟有效的方法［12，19］。文献［20?23］在对传统单极子波叠加法传递矩阵的病态性进行分析后发现，由于单极子向外辐射的波函数是以球面形式衰减的自由场格林函数，因此当等效源点之间或全息测量点之间的位置仅有微小改变时，格林函数的大小变化过于平缓，导致传递矩阵因不同行或列近似相等而病态。进而，文献［22?23］从优化波函数波阵面衰减速度的角度，提出了一种利用射线波函数替代传统球面波函数以改善传递矩阵病态性的射线波函数法，在改善传递矩阵病态性方面取得了一定效果。然而，该方法采用格林函数的方向导数作为射线波函数，这种类型射线波函数的解析表达式较为复杂，不仅计算效率低且难以计算到高阶导数。此外，该方法还将射线波函数的主指向设置为各等效源到其对应测点的方向，即要求等效源和测点无论是在数量还是分布方式上均必须一一对应，这无疑限制了射线波函数法在实际工程中的应用。

为了提高射线波函数法的计算效率和灵活性，本文采用（n， 0）阶正交球面波源构造了一种新型的射线波函数，该波函数无需求导运算，仅利用球Hankel函数和Legrend多项式的递推形式即可方便地计算出其任意阶的表达式，不仅计算效率得到大幅提高，而且不会增加计算难度和复杂性。与此同时，本文通过分析射线波函数对传递矩阵病态性的改善机理，提出一种射线波函数的主指向设置方法，使等效源与测点的布置不再受到数量和分布均须一一对应的限制，大大提高了射线波函数法在应用中的灵活性。最后，通过3个不同声源的数值仿真算例对比验证了本文方法在声场重建中的有效性和准确性。

1 理论部分

1.1 传统单极子波叠加法的病态性及射线波函数法

波叠加法的基本思想是：声源向外辐射的声场可由连续分布于其内部的等效源所辐射的声场叠加代替，该思想可用下式描述［7］：

式中 p（r）为声源在空间r处辐射的声压；σΦ（rE）为位于rE处的等效源源强，下标Φ表示其波函数为Φ（r，rE），根据等效源的类型不同，可为单极子等效源、偶极子等效源、单?偶极子组合型等效源等；Ω为等效源在声源内部的分布区域，SE为该区域的边界，如图1（a）所示。

为了便于计算，通常将等效源布置在虚拟边界SE上，并在积分离散时将每个单元内的等效源源强和波函数均视为常数，且配置在每个单元的中点，如图1（b）所示。式（1）离散为如下形式［5?6］：

式中 qΦ（rEj）为第j个离散等效源的源强；Φ（r，rEj）为该等效源对应的波函数。

由式（2）可知，只要确定了N个等效源的源强qΦ（rE1），qΦ（rE2），…，qΦ（rEN），即可计算出声源在任意场点r处辐射的声压。在基于波叠加法的近场声全息技术中，求解源强所需的方程组可以通过测量声源近场区域的声压或振速来建立。假设全息测量信息为声压，等效源为单极子，利用式（2）可得如下矩阵方程：

式中 pH=[p（rH1）p（rH2）…p（rHM）]T为M×1的测量声压列向量；[GH]ij=G（rHi，rEj）为M×N的全息测量数据与等效源强间的传递矩阵，G（rHi，rEj）为自由场格林函数，即单极子波函数；QG=[qG（rE1）qG（rE2）…qG（rEN）]T为N×1的单极子等效源强列向量。

求解方程（3），得源强向量QG为：

式中 G+H表示矩阵GH的广义逆。

为获得式（4）的最小二乘解，要求测点数M不小于等效源数目N，即M≥N。此外，在传统波叠加法中，一般采用单极子作为等效源，其波函数为自由场格林函数G（rHi，rEj）。文献［20?23］在对由该函数所构成的传递矩阵GH的病态性进行分析后指出，由于格林函數G（rHi，rEj）是一个只与两点距离有关且以球面衰减的波函数，因此当不同测点或等效源点之间的位置仅有微小改变时，矩阵的不同行或列的元素将近似相等，导致传递矩阵GH因向量间的强线性相关性而病态。对于该病态问题，一般是借助正则化方法抑制测量误差的放大并以此稳定求解过程。常用的正则化方法包括截断奇异值方法（TSVD）、Tikhonov正则化方法［13］等。但正则化方法本质上是通过类似滤波的方法将传递矩阵中对误差敏感的小奇异值项滤除。如果传递矩阵病态严重，在正则化时就需选取较大的截断点或正则化参数等以滤除更多的奇异值项，这将会加剧声场细节信息丢失，导致重建复杂声源或高频声场时精度下降。因而，为保证重建精度，应尽量改善波叠加法传递矩阵的病态性。

为了改善上述问题，文献［22?23］中提出了一种利用强指向性波函数替换球面形式波函数以提高声场重建稳定性的方法，并称之为射线波函数法。其基本原理是将传统波叠加法中单极子等效源辐射的球面波函数替换为满足Helmholtz方程和Sommerfield辐射条件且主值指向等效源对应测点的射线波函数，如图2所示。这样一来，等效源辐射的声波将仅在其对应测点处具有较大的声波激励，并生成较大的主对角元素，而在非对应测点处的声波激励则快速衰减，即生成较小的非对角元素，进而得到一个主对角元素占优的良态传递矩阵，以此提高声场重建的稳定性。

1.2 以往射线波函数法的缺陷及改进

在对文献［22?23］中提出的射线波函数法进行深入研究后发现，该方法虽然在改善传递矩阵病态性方面具有一定优势，但仍有以下两个缺陷：

（1）波函数计算效率方面的缺陷。文献［22?23］均是利用格林函数的方向导数作为射线波函数，该类型波函数的高阶导数（超过6阶后）解析表达式非常复杂，导致计算效率较低甚至无法计算。

（2）射线波函数主指向设置方面的缺陷。文献［22?23］中均是将射线波函数的主指向设为等效源的对应测点方向，这样设置虽然可以形成主对角占优形态良好的传递矩阵，但要求等效源数量与测点数量相同，且在运算过程中须保证这两组节点的编号始终一一对应。而在实际应用中，为获得更高的计算精度，等效源最好均匀布置在与声源面共形的虚拟面上，同时为了便于制造和降低成本，全息测量面则通常为规则形状且测点规则分布，这必然难以保证全息测点与等效源点之间的一一对应。

本文的主要内容则是针对以上两点提出如下改进办法：

（1）对波函数计算效率方面的改进——重新构造射线波函数。源模拟技术的研究表明，只要是满足Helmholtz方程和Sommerfeld辐射条件的解析函数Φ（r，rE）均可作为等效源波函数［5?6］。因而射线波函数的选取并不局限于格林函数的导数，本文将利用Helmholtz方程在球坐标系下的基本解，即正交球面波源重新构造射线波函数。

（2）对射线波函数主指向设置方面的改进。实际上，由1.1节中图2的分析可知，射线波函数能够降低传递矩阵线性相关性的主要原因在于它在非主指向的衰减速度远快于格林函数的球面波。因此，即便不采取等效源与主测点一一对应的设置方式，理论上仍应能显著降低传递矩阵因线性相关性过强导致的病态。在应用中，为了保证等效源辐射的声波在空间中分布均匀，只要各射线波函数的主指向均匀向外分散即可。例如，对于封闭的等效源面和全息面，可将射线波函数的主指向设置为从等效源面所包含空间的几何中心到各等效源连线的方向，如图3（a）中的l1，l2，l3，…，lN所示；对于非封闭的等效源面和全息面，则可以设置为等效源面的外法向方向，如图3（b）中的l1，l2，l3，…，lN所示。在后文的方法验证部分将通过仿真算例进行验证。

1.3 基于正交球面波源的射线波函数

正交球面波源是由不同阶次的球谐函数和球Hankel函数构成，其具体表达式为［24］：

式中 hn（?）为第一类或第二类n阶球Hankel函数，本文采用第二类n阶球Hankel函数，即h（2）n（?）；Ymn（?）为归一化的球谐函数；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，θ）：

上式即为带有参数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轴作一局部坐标系Ox'y'z'。此时，只要将射线波函数Dn（r，θ）中的变量θ和r替换为局部坐标系Ox'y'z'中的变量θ'和r'，即可使其指向z'轴，即lj方向。

利用式（6），可得rEj处单位强度等效源在场点r处辐射的声压为：

为方便起见，将上式中的局部坐标变量r'和θ'用全局坐标变量rEj，r和lj表示如下：

将上式中的场点位置r替换为全息面测点位置rHi，即可得到等效源与全息面间的传递矩阵DH，其元素为：

进而式（3）和（4）的矩阵方程可分别改写为：

式中 QD=[qD（rE1）qD（rE2）…qD（rEN）]T为射线波函数所对应的等效源强向量。

求解出等效源强后，再利用下式即可重建场点r处的声压：

值得一提的是，在n=0和n=1时，式（6）分别对应单极子等效源和偶极子等效源［25］。此时，式（12）将变为传统的单层势和双层势波叠加法。

2 方法验证

2.1 仿真模型说明

为验证本文所提方法在声场重建中的效果，该部分设计了3种不同类型声源的仿真模型，分别为正四面体辐射体、两端带球帽的圆柱辐射体和四边简支矩形钢板声源。需要指出的是，由于球面波源通常在单点多极法中被用于计算球形或近似球形的声源［5?6，24，26?27］，而本文则是将其应用于波叠加法中，因此为了验证由球面波源构造的射线波函数在波叠加法中依然适用于任意形状声源的优势，上述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叠加产生。为了增加声场的复杂度，将两个点声源设置为不同的多极子声源，其具体的表达式为：

式中 坐标变量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给出了各阶波函数对应的重建声压相对误差曲线及传递矩阵条件数曲线，其中，相对误差由下式计算：

式中 pn表示第n阶波函数的重建声压向量；p表示解析声压向量。

由图9（a）可以发现，当波函数阶数n=0和n=1时（即传统波叠加法），重建误差较大，超过了20%。而且当波函数的阶数在0～4之间时，重建误差基本呈下降趋势，并在4阶时最小，误差小于10%。通过对比图9（b）的条件数曲线可知，重建误差下降是由于传递矩阵条件数逐渐减少，因此正则化后损失的声场信息随之减少，重建结果更為精确和稳定。而当波函数的阶数大于4阶时，传递矩阵的条件数虽然仍保持下降趋势，但重建误差却反而逐渐上升。这是因为过高的阶数会使得射线波函数的指向性太强，并由此导致各等效源辐射的声波过于集中在其对应主指向的射线束内，而在非主指向上将因波函数衰减太快而产生较大误差，因此无法正确描述真实声场。由此可见，射线波函数的指向性并不是越强越好，而是应在合理的范围内选取。一个简单的方法是利用“辅助面法［24］”进行选取，即在声源面与全息面间设置若干辅助测量点，然后计算各阶射线波函数（通常只需要计算到前10～15阶）在辅助测点处的重建声压，最后将重建声压与该点的测量声压进行对比，取两者误差为最小时所对应的阶数即可。由于篇幅所限，本文不展开讨论波函数阶数的选取。实际上，由球面波源构造而成的射线波函数作为Helmholtz方程在球坐标系下的基本解，其阶数的选择需要综合考虑声源特性、NAH模型配置及噪声性质等［28］。

2.3 仿真2：两端带球帽的圆柱辐射体的声压重建

以坐标原点O为几何中心设置一两端带球帽的圆柱壳长条辐射体，其中圆柱部分长度为2a=1 m，圆柱和球帽的半径均为a=0.5 m，如图10（a）所示。假设该辐射体表面声压由分别置于两个球帽球心处的点声源S1和S2叠加产生。与仿真1类似，将两个点声源设置为多极子，其中：

式中的坐标变量说明见仿真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不同，该重建模型中最小误差对应的波函数阶数为n=5。由图12（b）的条件数曲线可以发现，虽然射线波函数并未指向等效源的对应测点，但传递矩阵的条件数仍然呈现良好的下降趋势，这表明本文在1.3节提出的射线波函数指向设置方法是有效的。

2.4 仿真3：矩形简支板辐射声压的重建

在波叠加法求解声场外问题的理论中，通常要求给定边界条件的全息面为一封闭曲面，但在实际工程中，重建由板、壳等结构振动产生的声场是一种常见的应用场景，此时一般采用非封闭的平面全息面测量声场信息。这一节的仿真将利用一个被无限大刚性障板围绕的简支板声源验证本文方法在全息面为非封闭平面时的声场重建效果。

如图13所示为一左下角位于坐标原点、长方向和宽方向分别与y轴和x轴重合的长方形钢板。长宽厚尺寸设置为1 m×0.5 m×0.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相同的重建声压幅值分布云图、声压相对误差曲线及传递矩阵条件数曲线外，为了突出各阶波函数对简支板各区域声场的重建效果，还给出了重建面上的相对误差分布云图，上述结果如图14～16所示。从图14的声压幅值分布云图上看，各阶波函数的重建声压均与解析声压吻合得较好。但通过对比图15中的相对误差云图则可以发现，当波函数阶数为n=0和n=1时（即传统波叠加法），重建误差云图中有大面积误差≥15%的深红色。而当波函数为n=2～8时，云图中深红色的大误差区域显著减少，并主要分布在解析声压接近于0处和简支板边缘处，这是由于计算相对误差时分母接近0对误差产生的放大作用和声场信息泄露所导致的，其余大部分区域均为误差0%～5%的深蓝色和浅蓝色。由图16（a）可以发现，该重建模型的最小误差对应的波函数阶数为n=4，最小误差在5%左右，且阶数在n=2～8的范围内重建误差均相差不大，都在5%～10%之间。从图16（b）中可以看到，传递矩阵的条件数曲线仍然呈下降趋势。该仿真表明，对于非封闭的重建模型，本文方法仍具有一定优势。

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）： 340?344.

2Veronesi 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 BEM?based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries［J］. 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）： 39?48.

5Ochmann M. The source simulation technique for acoustic radiation problems［J］. Acta Acustica United with Acustica， 1995， 81（6）： 512?527.

6Ochmann M. The full?field equations for acoustic radiation and scattering［J］. The Journal of the Acoustical Society of America， 1999， 105（5）： 2574?2584.

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）： 2433?2438.

8Du B K， Zeng X Y， Vorl?nder M. Multizone sound field reproduction based on equivalent source method［J］. Acoustics Australia， 2021， 49（2）： 317-329.

9Chen J， Liu C， Zhang X Z， et al. An approach for indoor prediction of the pass?by noise of a vehicle based on the time?domain equivalent source method［J］. Mechanical Systems and Signal Processing， 2021， 146：107037.

10刘延善， 曾向阳， 王海涛. 封闭空间声场重构的多层等效源法［J］. 声学学报， 2020， 45（3）： 367?376.

Liu Yanshan， Zeng Xiangyang， Wang Haitao. 3D sound field reconstruction for the enclosed cavity using the multilayer equivalent sources method［J］. Acta Acustica， 2020， 45（3）： 367?376.

11Valdivia N P. Advanced equivalent source methodologies for near?field acoustic holography［J］. Journal of Sound and Vibration， 2019， 438： 66?82.

12BAI M R， CHEN C C， Lin J H. On optimal retreat distance of equivalent source method?based near?field acoustical holography［J］. The Journal of the Acoustical Society of America， 2011， 129（3）： 1407?1416.

13Williams E G. Regularization methods for near?field acoustical holography［J］. The Journal of the Acoustical Society of America， 2001， 110（4）： 1976?1988.

14Valdivia N P， Williams E G. Study of the comparison of the methods of equivalent sources and boundary element methods for near?field acoustic holography［J］. The Journal of the Acoustical Society of America， 2006， 120（6）： 3694?3705.

15Gounot 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）： 3195?3205.

16Gounot Y J R， Musafir R E. Genetic algorithms： a global search tool to find optimal equivalent source sets ［J］. Journal of Sound and Vibration， 2009， 322（1?2）： 282?298.

17Gounot Y J R， Musafir R E. Simulation of scattered fields： some guidelines for the equivalent source method ［J］. Journal of Sound and Vibration， 2011， 330（15）： 3698?3709.

18WU S W， XIANG Y. Location optimization of monopole equivalent sources in wave superposition method ［J］. International Journal of Acoustics and Vibration， 2018， 23（2）： 254?263.

19LEE S. Review： the use of equivalent source method in computational acoustics［J］. Journal of Computational Acoustics， 2017， 25（1）： 1630001.

20张永斌. 基于等效源法和质点振速测量的近场声全息技术［D］. 合肥： 合肥工业大学， 2010.

Zhang Yongbin. Research on nearfield acoustic holography based on the equivalent source method and the measurement of particle velocity［D］. Hefei： Hefei University of Technology， 2010.

21Fernandez?Grande E， Xenaki A， Gerstoft P. A sparse equivalent source method for near?field acoustic holography［J］. The Journal of the Acoustical Society of America， 2017， 141（1）：532?542.

22向宇， 石梓玉， 陆静， 等. 基于波叠加法的非共形近场声全息波函数的构造与选择［J］. 振动与冲击， 2020， 39（15）： 183?192.

Xiang Yu， Shi Ziyu， Lu Jing， et al. Construction and selection of nonconformal near?field acoustic holography wave function based on wave superposition method［J］. Journal of Vibration and Shock， 2020， 39（15）： 183?192.

23石梓玉， 向宇， 陆静， 等. 一种提高声场重构稳定性的射线等效源法［J］. 广西科技大学学报， 2019， 30（3）： 1?7.

Shi Ziyu， Xiang Yu， Lu Jing， et al. A ray equivalent source method for improving the stability of acoustic field reconstruction［J］. Journal of Guangxi University of Science and Technology， 2019， 30（3）： 1?7.

24毕传兴， 陈心昭， 陈剑， 等. 正交球面波源边界点法及其在声全息中的应用［J］. 科学通报， 2004， 49（13）： 1322-1331.

Bi Chuanxing， Chen Xinzhao， Chen Jian， et al. Orthogonal spherical wave source boundary point method and its application to acoustic holography［J］. Chinese Science Bulletin， 2004， 49（13）：1322-1331.

25Williams E G. Fourier Acoustics： Sound Radiation and Nearfield Acoustical Holography［M］. London： Cambridge University Press， 1999.

26Gomes J， Hald J， Juhl P， et al. On the applicability of the spherical wave expansion with a single origin for near?field acoustical holography［J］. The Journal of the Acoustical Society of America， 2009， 125（3）： 1529?1537.

27毕传兴， 胡定玉， 徐亮， 等. 采用球面波叠加法还原自由声场的方法［J］. 声学学报， 2014， 39（3）： 339?346.

Bi Chuanxing， Hu Dingyu， Xu Liang， et al. Recovery of the free field in a noisy environment by using the spherical wave superposition method［J］. Acta Acustica， 2014， 39（3）： 339?346.

28Wu H J， Jiang W K， Zhang H B. A mapping relationship based near?field acoustic holography with spherical fundamental solutions for Helmholtz equation［J］. Journal of Sound and Vibration， 2016， 373： 66?88.

Wave superposition method with ray wave functions for near?field acoustic holography based on orthogonal spherical wave sources

SHI Zi?yu 1，2 ?XIANG Yu 1，2 ?LU Jing 1，2WANG Yu?jiang 1，2

1. Guangxi Key Laboratory of Automobile Components and Vehicle Technology， Guangxi University of Science and Technology， Liuzhou 545006， China；

2. School of Mechanical and Automotive Engineering， Guangxi University of Science and Technology， Liuzhou 545006， China

Abstract In the near-field acoustic holography based on the traditional wave superposition method， the monopole radiating the spherical wave is often used as the equivalent source， which is easy to cause the ill conditioned transfer matrix. Replacing the spherical wave function with the ray wave function can effectively improve the ill conditioned transfer matrix. However， the previous ray wave function method uses the directional derivative of Green's function as the wave function， which has complex analytical expression and low computational efficiency. In addition， the wave function direction setting of the previous methods has high requirements for the node distribution， which limits the flexibility of its application. In order to solve the above problems， a series of ray wave functions are reconstructed by using （n， 0） order spherical wave source， this ray wave function can be easily calculated to any order by using the recursive form of spherical Hankel function and Legendre polynomial， and the efficiency is greatly improved. By improving the main direction setting of ray wave function to make it more flexible in practical use， a ray wave function wave superposition method based on orthogonal spherical wave source is proposed. Three numerical simulations， i.e.， radiator of regular tetrahedral， radiator of cylinder with two spherical caps and simply supported rectangular steel plate acoustic source， are performed to verify the sound field reconstruction effect in both traditional method and proposed method. The simulation results show that in the high frequency with abundant sound field details， due to the serious ill conditioned transfer matrix， the traditional wave superposition method cannot guarantee the reconstruction accuracy in the three simulations even with the help of regularization method， and the error is about 20%. The established ray wave function method effectively reduces the condition number of the transfer matrix and improves the ill condition of the system. Therefore， higher reconstruction accuracy is obtained， and the error is between 5% and 10%， which shows that the method in this paper is more superior.

Keywords near?field acoustic holography; wave superposition method; orthogonal spherical wave source; ray wave function