CN102089634A

CN102089634A - 重建声学场

Info

Publication number: CN102089634A
Application number: CN2009801264451A
Authority: CN
Inventors: J.哈德
Original assignee: Bruel and Kjaer Sound and Vibration Measurement AS
Current assignee: Hottinger Bruel and Kjaer AS
Priority date: 2008-07-08
Filing date: 2009-06-26
Publication date: 2011-06-08
Anticipated expiration: 2029-06-26
Also published as: WO2010003837A1; EP2297557A1; JP5220922B2; JP2011527422A; EP2297557B1; US8848481B2; US20110164466A1; CN102089634B

Abstract

公开了一种重建声场的方法。该方法包括接收在一组测量位置测量的第一声学量的测量值；根据平面波的叠加计算目标位置的第二声学量。该方法包括存储一组各个函数的内插表示，每个函数是两个或更少输入参数的函数；并且计算包括计算一组关联函数中的每一个，每个关联函数指示在所述测量位置中的第一位置处的平面波与第二位置处的平面波的关联，作为从该组内插表示获得的值的线性组合。

Description

重建声学场

技术领域

本发明涉及声学场的重建。

背景技术

近场声全息（NAH）是一种用于声幅射的3D可视化和基于声源附近的表面上的测量的精确噪声源定位的非常有用工具。它的还重建渐逝波分量的能力确保非常高的空间分辨率。

已知的近场声全息方法基于可分离坐标系中水平面上的规则网格测量，允许用空间离散傅立叶变换（DFT）执行NAH计算，参见例如E. G. Williams、J. D. Maynard和E. J. Skudrzyk等人在J. Acoust. Soc. Am. 68, 340-344 (1980)发表的“Sound source reconstruction using a microphone array”。因为DFT的使用，该处理非常快，但是使用DFT的负面效应包括严重的空间开窗效应，除非测量区域利用高声压完全覆盖这些区域。在一些情况下，不能实现对测量区域的这种需求，并且在许多情况下，必要的尺寸会变得过大。

已经提出一组技术来降低空间开窗效应，同时仍保持DFT空间处理，但是以增加复杂性和计算需求为代价，参见J.Hald在Inter-Noise 1994论文集中发表的“Reduction of spatial windowing effects in acoustical holography”。典型地，迭代过程首先被用于外推（extrapolate）被测区域之外所测量的声压，之后在扩展的数据窗口上应用基于DFT的全息方法。

已提出其他方法来设法避免使用空间DFT以及提供所需测量区域的减小。

一种这样的方法是Helmholtz方程最小二乘（HELS）方法，其使用根据球面波函数的声场的局部模型，参见例如Z.Wang和S.F.Wu的US 6,615,143的“Helmholtz equation-least-squares method for reconstructing the acoustic pressure field” J. Acoust. Soc. Am, 102(4), 2020-2032 (1997)或者S. F. Wu的“On reconstruction of acoustic fields using the Helmholtz equation-least-squares method”J. Acoust. Soc. Am, 107, 2511-2522 (2000)。然而，因为仅具有公共原点的球面波函数被用来表示声场，所以误差将被引入到关于源表面的声场重建中，除非源表面也是球面的并且以相同的原点为中心。该现有技术方法的另一缺点是估计问题的需求没有被充分确定，意味着测量位置的数目必须大于或等于在场表示中使用的球面波函数的数目。这可以需要大量的测量位置以获得足够精确的模型。第三个缺点是波函数的缩放（scaling）不以下述方式来应用，即在模型区中具有较强衰减的函数被缩放到相同区中的较低振幅。因为缺少该缩放，比如Tikhonov正则化之类的传统正则化方法不会适当地工作。相反，上述现有技术方法结合最小二乘解而无需正则化来应用从计算方面讲昂贵的迭代搜索来寻找球面波展开式的最优截断。

先前提出的另一种方法是在R.Steiner和J.Hald于Intern. J. Acoust.Vib.6,83-89(2001）发表的“Near-field Acoustical Holography without the errors and limitations caused by the use of spatial DFT” 中提出的统计最优近场声全息（SONAH）方法。该现有技术方法基于平面波函数，并且通过使用以下述方式定义的传递矩阵计算接近测量表面的映射表面上的声学量，所述方式为以最优平均精确度来投影（project）所有传播波和渐逝波的加权组。从1）所有测量位置对之间的自关联和互关联的自关联矩阵以及2）从每个测量点到每个重建位置的互关联的互关联矩阵获得传递矩阵。自关联和互关联处于平面传播和渐逝波的域中。

已在J. Hald发表于Inter-Noise 2003论文集中的“Patch near-field acoustical holography using a new statistically optimal method”和J. Hald发表于Euronoise 2006论文集中的“Patch holography in cabin environments using a two-layer handheld array and an extended SONAH algorithm”中描述了该框架在特定平面测量几何结构中的应用。

这些现有技术文档公开了垂直于测量平面的粒子速度分量和压力的计算。然而，该方法从计算方面来讲是昂贵的，因为在这些现有技术文档中公开的关联公式的数值计算从计算方面来讲是不高效的。

发明内容

本文公开了一种重建声场的方法，所述声场由至少一个声源生成。该方法的实施例包括：

－接收在一组测量位置中的各个位置处测量的至少第一声学量的一组测量值；

－计算一组关联函数中的每一个，每个关联函数指示在所述测量位置的第一位置处的至少一组平面传播和渐逝波和在第二位置处的至少一组平面波的关联；

－至少根据所计算的关联函数和该组测量值来计算指示在目标位置处的重建声场的第二声学量；

其中该方法还包括存储用于内插对于关联函数的各个贡献（contribution）的一组内插函数的表示，每个贡献是两个或更少输入参数的函数；以及其中计算该组关联函数中的每一个包括计算作为从所存储的表示获得的值的线性组合的每个关联函数。

发明人明白本文所述的方法可以以一致且计算高效的方式应用于一般的几何结构以及许多不同声学量的计算。进一步得到的结果是针对各种各样的测量几何结构和声学量的关联函数可以被表示为1维和/或2维函数的线性组合，所述1维和/或2维函数可以至少局部地被两个或更少输入值的各个内插函数近似。特别地，因为内插函数不多于二维，所以实现计算上的高效以及还有精确的内插，并且所存储的内插函数仅需要小的存储容量。因此，通过使用结合内插的表格或近似函数，数值计算所基于的公式具有非常适合于高效计算的形式。本文所描述的方法的实施例可以比现有技术方法执行快大约9倍。

声源可以是发射或反射声学辐射的任何对象。对象可以是任意形状的，并且声音可以是任何种类的声音，例如噪声、可听见的声音、不可听见的声音（例如超声或次声等等）或者其组合。包括一个或多个声源的区可以由将（一个或多个）声源从无源区分离的声表面定义。典型地，源表面可以是包括（一个或多个）声源的物理对象的表面。

该过程的实施例包括定义一个或多个虚拟源平面，每一个虚拟源平面表示一组平面传播和渐逝波函数；并且其中平面波具有相应波数并且在相对于虚拟源平面的方向上传播和/或渐逝。一般来说，对于虚拟源平面，可以定义指向远离源并且进入包含测量和重建位置的区（该区还将被称为测量区）的法线（normal）。对应的传播平面波通常可以覆盖来自法线方向以及偏离该方向高达90度的所有传播方向。渐逝波在该法线方向上即进入测量区以指数形式衰减。渐逝波可以形成不同衰减率和不同切线传播方向的连续谱。

例如，可以相对于包括声源的对象的表面的一部分（即相对于声表面）来定义虚拟声平面。例如，虚拟源平面可以被定义为从源表面间隔开短距离并且位于该对象内部的表面。

在测量位置处的声学量可以由任何适合的声学测量设备（例如传声器（microphone）、水听器（hydrophone）、压力梯度换能器、粒子速度换能器等等或其组合）测量。在一些实施例中，由声学测量设备阵列（例如以规则或不规则的网格布置的一组这种设备，诸如两维或三微网格）执行测量。所测量的第一声学量可以是声压、声压梯度、粒子速度和/或等等。

重建声场可以由适合的数据组来表示，所述数据组指示第二声学量（例如声压、声强度或粒子速度）的空间分布。第一和第二声学量可以是相同的量或不同的量。数据组可以被表示为映射，例如直接在任意表面的实际表面几何结构上的第二声学量的映射，所述任意表面例如是要被分析的声音/噪声发射对象的表面或者接近对象的表面。

因为硬表面的顶部上的粒子速度紧密对应于表面自身的实际振动，所以结果可以直接被用于与结构模型的关联。在一些实施例中，声场由保角映射表示。声场参数可以被重建并且被映射到表面上，通常被映射到源表面上或者源表面附近。可替换地或附加地，可以生成重建声场的特性的其他表示，例如波函数表示的有效性的3D区中的单个位置处的谱。

在此处所描述的方法实施例中，计算点处的粒子速度或压力的重建包括1）所有测量位置对之间的自关联和互关联的自关联矩阵的计算以及2）从每个测量点到目标位置的互关联的互关联矢量的计算。自关联矩阵可以针对所有目标位置被再次使用。在场表示中使用的平面传播和渐逝波函数的域中定义自关联和互关联，可以从目标位置和测量位置中选择第二位置。特别地，计算可以包括各个测量点之间的该组平面波的关联函数的自关联和互关联矩阵的计算以及各个测量点和目标点之间的该组平面波的关联函数的互关联矢量的计算。

术语内插函数（用于内插各个贡献）意图指局部近似于该贡献的函数组或任何函数。当内插函数的表示包括一组查找表时（每个查找表包括预先计算的相应函数的局部近似），提供特别高效的计算。每个查找表可以指示由两个输入参数编索引的函数值的二维阵列，并且每个函数值可以指示由两个输入值定义的输入点处的二维函数的值。各个函数中的每一个函数可以是包括由两个输入参数来参数化的积分的二维积分函数。

因此，在其他点处的函数值可以通过制成表格的函数值之间的内插（例如三次内插）而高效地获得。因为查找表是低维度的，即具有二维或更低的维度，所以内插过程是精确的并且计算高效的，并且查找表仅需要小存储并且可以利用很少的计算资源来生成。

两个输入参数中的第一输入参数可以指示在平行于虚拟源平面的平面中的第一测量位置和第二位置之间的由波数指示的因子缩放的投影距离。两个输入参数中的第二参数可以指示从虚拟源平面到第一测量位置和第二测量位置的相应距离的线性组合（通常为和或差），由波数指示的因子缩放。在一些实施例中，该参数化允许关联函数由一组预先计算的二维函数来表达。第二输入参数可以指示投影到虚拟源平面的法线（通常为z轴）上的第一测量位置和第二测量位置之间的并且由指示波数的因子缩放的投影距离。

在一些实施例中，根据预定的加权方案，每个平面波函数被相应加权因子缩放/加权。例如，预定的加权方案可以从均匀加权和方向性加权中选择。在均匀加权中，每个平面波被相同的常数因子缩放。在平面波的方向性加权中，例如通过向在传播方向上具有平行于源平面的较大分量的平面波提供较高的加权因子，每个平面波被相应加权因子缩放，所述相应加权因子取决于平面波的传播方向。这样的波还被称为离轴波，因为它们不平行于虚拟源平面的表面法线。已发现方向性加权来改进重建方法的精确度。

根据一些实施例，已发现对于给定的加权方案，可以基于六个二维函数的制成表格的表示的单个组来计算针对各种测量和重建几何结构的声压力和粒子速度二者。仅函数参数随着组结构而变化，必须计算在所述函数参数处的内插函数值。

在一些实施例中，选择虚拟源平面来定义：第一空间域（也被称为测量域），其包括测量位置和目标位置；以及第二空间域（也被称为源域），其不同于（相对于）第一空间域，所述第二空间域通常包括大多数声源，即测量位置和目标位置位于测量域中，并且在该测量域中重建该声场。尽管虚拟源平面通常被定义源平面内的短距离，但是为了重建过程，由虚拟源平面定义的测量域被假设成无源的。

在一些情况下，由单个虚拟源计划定义无源空间域是不可能或者不现实的。因此，在一些实施例中，该方法基于两组或更多组平面波，每一组都在相对于相应虚拟源平面的方向上传播和渐逝。例如，虚拟源平面可以被定义为一对或多对间隔开的平行虚拟源平面，其中无源区位于虚拟源平面的各自对的每一对之间。发明人已明白本方法和制成表格的内插函数的相同组也可以被应用于这种情况，从而提供效率和灵活的重建方法。

该组测量位置可以被布置在一个或多个测量平面中，例如在单个平面或两个或多个平行平面中。在每个平面内，测量位置可以以不规则的模式或者任何其他适合的方式布置在规则网格上。此外，本文所述的方法还可以被应用于非平面测量几何结构，即其中测量位置不位于一个或多个平行平面中但是例如位于弯曲表面上的布置。

在一些实施例中，当测量和计算点被布置在平行于所有虚拟源平面的各个平面中时可以获得进一步增加的效率。在这样的设置中，针对平面波的每个频率，仅出现对于关联函数的贡献的输入参数之一的一些值，因此使得1D查找表的生成具有吸引力，其中每个表格对应于所述输入参数的相应值。因此，对于某些几何结构，该方法甚至基于一维内插函数，并且因此在计算方面特别高效地执行。

当第二声学量是具有在平行于虚拟源平面的方向上的分量的粒子速度时，有助于对具有弯曲源表面的对象的改进分析。特别地，该实施例允许横跨弯曲重建表面分布（例如与源表面保角）的目标位置中的源场的改进重建。因为重建过程的该实施例提供与虚拟源平面平行的粒子速度的分量，所以可以横跨弯曲重建表面确定垂直于重建表面的粒子速度的分量（和/或强度）。

此处所公开的方法实施例被看作基于一组测量量在目标点处执行声场参数的线性预测，预测系数可以被确定以便在最小二乘意义下提供对所有平面传播波和渐逝波的加权组的最佳估计。最小二乘拟合可以被看作确定一组复系数的数值过程，该组复系数以适当方法解线性等式的过定（over-determined）组。

于是期望知道，平均来说在计算点处多好地预测这些平面波函数的值。因此，在一些实施例中，重建方法包括计算指示所计算的声学量的估计误差的量。可以使用先前描述的自关联和互关联（即利用少量的附加计算资源）来计算误差指示符。上述用于估计SONAH过程的预测精确度的方法已被证明在利用给定的测量几何结构来可视化并确定SONAH预测的有效性区的方面非常有用。

要注意，可以至少部分在软件或固件中实施上文和下文中描述的方法特征，并且在由诸如计算机可执行指令之类的程序代码装置的执行引起的数据处理设备或其他处理装置上实现上文和下文中描述的方法特征。在此处和下文中，术语处理装置包括被适合地适于执行上述功能的任何电路和/或设备。特别地，上述术语包括通用或专用可编程微处理器、数字信号处理器（DSP）、专用集成电路（ASIC）、可编程逻辑阵列（PLA）、现场可编程门阵列（FPGA）、专用电子电路等等，或其组合。

可以以不同方式来实施本发明的实施例，所述不同方式包括上文和下文中所述的方法、系统、设备和产品装置，每一个都产生结合首次提出的方法而描述的一个或多个优点和益处，并且每一个都具有对应于结合首次提出的方法而描述的和/或如在从属权利要求中所公开的实施例的一个或多个实施例。

特别地，用于重建声场的处理装置的实施例包括用于接收在一组测量位置的各个位置处测量的至少第一声学量的一组测量值的接口；以及被配置成执行下述内容的处理单元：

-计算一组关联函数中的每一个，每个关联函数指示在所述测量位置的第一位置处的至少一组平面波与在第二位置处的至少一组平面波的关联；以及

-至少根据所计算的关联函数和该组测量值来计算指示目标位置处的重建声场的第二声学量；

其中所述处理单元包括用于存储用于内插对于关联函数的各个贡献的一组内插函数的表示的存储介质，每个贡献是两个或更少输入参数的函数；并且其中所述处理单元适于计算作为从所存储的表示中获得的值的线性组合的该组关联函数中的每一个关联函数。

将会认识到，可以在用于执行重建过程（例如作为计算机程序的一部分或者在单独的文件或多个文件中）的计算机程序的软件编译时间生成所存储的表示，或者可以在运行时间由结果所得的程序代码生成所存储的表示，或者可以作为其组合。例如，该表示可以由计算机程序或安装程序生成并且在该计算机程序的第一次执行/安装期间被存储在计算机上。

用于重建声场的系统可以包括如上文和下文所公开的装置，以及用于在一组测量位置处测量第一声学量并且以通信连接的方式可连接到所述装置以将所测量的值转发到所述装置的一组换能器。

计算机程序可以包括程序代码装置，当在数据处理系统上执行程序代码装置时其适于使得所述数据处理系统执行上文公开的方法的步骤。该计算机程序可以被存储在存储装置上或者被体现为数据信号。存储装置可以包括用于存储数据的任何适合的电路或设备（例如RAM、ROM、EPROM、闪速存储器）、磁或光存储设备（例如CD ROM、DVD硬盘）、和/或等等。

附图说明

将参考附图从下文的描述的实施例中阐明并明白上述和其他方面，在附图中：

图1示出用于重建声场的系统的示意性框图。

图2示出计算重建声场的过程的流程图。

图3说明在具有仅一个虚拟源平面、具有限制匀质且无源区的两个平行虚拟源平面、以及具有盒形匀质且无源区的所有六个侧面上的虚拟源平面的情况下的测量几何结构。

图4说明在(k _x ,k _y)域和传播角度域中的平面波函数的密度。

图5示出对于两个不规则平面阵列的阵列平面中的声压预测的固有相对误差级E(r)。

图6示出在单层和双层8×8元素阵列的xz对称平面中的相对误差级。左列：单层阵列的固有误差级E(r)。中间列：由单层阵列估计的单极压力的相对误差。右列：双层阵列的固有误差级E(r)。

图7示出（a）垂直于8×8元素平面阵列的阵列平面的对称轴上的固有相对误差级E，并且图7b示出垂直于8×8元素平面阵列的阵列平面的对称轴上的振幅增益A。

图8示出在低频和高频下垂直于8×8元素平面阵列的阵列平面的对称轴上的固有相对误差级E(r)和振幅增益A(r)。曲线表示在重构过程中使用的不同动态范围D。

图9示出在低频和高频下垂直于8×8元素平面阵列的阵列平面的对称轴上的固有相对误差级E(r)和振幅增益A(r)。曲线表示在重建过程中使用的不同虚拟源平面距离d = - z ⁺，而定义渐逝波的实际内容的虚拟源距离在6cm处保持恒定。

图10说明在根据双层8×8阵列的中心平面中的单个单极点源来计算3D矢量强度的实施例中通过添加来自所有六个侧面而不是仅两个侧面的平面波函数来得到精确度的改进。所述点源位于阵列区域之外6cm。

遍及附图，凡是可能之处，相等的参考符号指的是相等或对应的元件、特征或部件。

具体实施方式

图1示出用于重建声场的系统的示意性框图。该系统包括一组声学接收机108和连接到所述声学接收机的分析单元103。

在下文中，声学接收机108还将被称为换能器。然而，将会认识到，该声学接收机可以是传声器、水听器、或用于测量诸如声压、声压梯度、粒子速度或其他线性量之类的声学特性的任何其他适合的设备。在图1的实施例中，换能器108被实施为换能器阵列102，其中以规则网格（例如1维、2维或3维网格）来布置该换能器108。尽管以规则网格的布置允许高效的数值实施，但是可以使用其他几何结构。在一些实施例中，甚至可以使用单个换能器来进行连续的不同点处的测量。这样的连续测量可以在例如其中声源是静止的，以及其中可以测量相对于一个或若干个参考信号的相位的情况下使用。可以根据要被分析的对象的几何结构复杂性和尺寸、感兴趣的频率范围以及期望的空间分辨率来选择换能器的数目和阵列的几何结构（例如换能器间间距）。

换能器阵列102被连接到分析单元103，使得该换能器108可以例如经由有线或无线信号连接将所测量的信号转发给分析单元。

分析单元103包括：用于接收和处理来自换能器阵列102的所测量的信号的接口电路104；与所述接口电路104进行数据通信的处理单元105；存储介质112；以及与所述处理单元105进行数据通信的输出单元106。即使在图1中被示出为单个单元，但是将会认识到分析单元103可以被物理地分成两个单独的设备（例如采集前端和计算机），或者甚至多于两个设备。类似地，将会认识到关于分析单元的不同子块而描述的功能可以被分成替代的或附加的功能或硬件单元/模块。

接口电路包括适合于从换能器108接收输出信号并且处理所接收的信号以用于处理单元105的后续分析的信号处理电路。接口电路104可以包括下述组件中的一个或多个：一个或多个用于放大所接收的信号的预放大器，一个或多个用于将所接收的信号转换成数字信号的模拟到数字（A/D)转换器，一个或多个滤波器（例如宽带滤波器）和/或等等。例如，接口电路可以提供作为输出数据的作为每一个换能器的频率的函数的振幅和相位。在一些实施例中，接口单元可以执行同时数据采集，并且于是由处理单元105来完成所有进一步处理，包括通常使用FFT的数据到频域的变换。

处理单元105可以是适合的编程微处理器、计算机的中央处理单元、或用于处理从接口单元104接收的信号的任何其他适合的设备，例如ASIC、DSP和/或等等。处理单元适于处理经由接口单元104接收的换能器信号，以便如本文所描述的那样计算重建的声场。

存储介质112可以包括用于存储指示一组内插函数的表示的数据的任何适合的电路或设备（例如RAM、ROM、EPROM、EEPROM、闪速存储器），磁或光存储设备（例如CD ROM、DVD、硬盘）、和/或等等。在图1中，存储介质被示出为单独的，但是与处理单元通信连接。然而，将会认识到，存储介质112还可以被体现为处理单元105的一部分，例如内部存储器。

输出单元106可以包括显示器或用于提供重建声场的视觉表示的任何其他适合的设备或电路，例如用于提供打印表示的打印机和/或打印机接口。可替换地或附加地，输出单元106可以包括用于传达和/或存储指示重建的声场的数据的任何适合的电路或设备（例如RAM、ROM、EPROM、EEPROM、闪速存储器），磁或光存储设备（例如CD ROM、DVD、硬盘）、有线或无线数据通信接口（例如到计算机或电信网络诸如LAN、广域网和因特网的接口）和/或等等。

分析单元103可以被实施为适合编程的计算机，例如包括适合的信号采集板或电路的PC。

该系统还可以包括连接到分析单元103的适合的位置检测设备110，以用于检测换能器108例如相对于参考坐标系的位置。位置检测设备可以包括测试台、框架或其他结构，可以将换能器阵列以预定位置和定向安装到其上，以及/或者阵列的位置和定向可以由例如适合的传感器手动或自动测量之处。可替换地，换能器阵列102可以包括位置检测设备或至少其一部分。例如，换能器阵列可以包括一卷（a coil of）磁场位置检测系统、陀螺仪（gyroscope）、光学位置检测系统和/或等等。

在操作期间，换能器阵列102被定位在对象100的表面附件，所述对象100包括发射声学辐射并且针对其来重建声场的声源111。可以根据要被分析的对象的几何结构复杂性和尺寸、感兴趣的频率范围以及期望的空间分辨率来选择换能器的数目、阵列的几何结构（例如换能器间间距）以及阵列和源表面之间的距离。例如，对于汽车引擎的声场的重建，可以使用具有3cm元件间间距的传声器的8×8矩阵阵列。然而，将会认识到也可以使用其他类型和尺寸的阵列。

例如由位置检测设备110来确定阵列102的位置，并且将其馈送到分析单元103中。阵列102的换能器108测量确定位置处的声压或另一适合的声学量，并且结果所得的换能器信号可以被发送给分析单元103。

例如，换能器阵列可以是具有集成位置检测设备的手持阵列，因而允许在对象周围分布的不同的可访问位置处的测量，并且减少在用于测试的源的特殊准备上所花费的时间最小。例如，当将声源（例如压缩机、泵或其他机器）映射到测试台上时，改变测试设置（例如燃料管路、电子接线（wiring）等等的布线）以提供对完整面的不受妨碍的访问所需的时间被降低，如果没有消除的话。另一典型应用可以在汽车车厢内部，其中3D阵列网格可以被用于能够区别所有方向上的源，例如可以使用双层阵列（例如包括8×8×2传感器）。

分析单元103根据输入数据计算重建的声场，并且存储和/或输出重建声场的表示。

将参考图2并继续参考图1来描述用于生成重建声场的过程的实施例。

为了说明在此处所述的重建过程，现在将提供所谓的SONAH方法的最初概要。

为了本描述，我们假设已在由均质流体占用的无源区Ω中的一组位置 r _i , i= 1,2,…,I处测量复空间谐波声压p(r _i )。我们考虑到基于所测量的压力来在Ω中重建声场的问题。为此，我们最初引入基波函数的离散组Ψ_n , n = 1,…,N，实现均质波等式并且能够以足够高精确度来表示可以存在于Ω中的所有声压场。一般来说，这些基波函数可以是平面波、柱面波或球面波，或者它们可以是来自基源（例如分布在Ω之外的单极或偶极点源）的声场。

在所谓的HELS方法中，重建的第一步骤是确定声场表示中的每个基波函数的振幅和相位。这可以通过需要基波模型尽可能精确地表示所测量的压力来完成：

（1）

a_n是复展开系数。如果I ≥ N，并且如果已适合地选择测量点，则可以由（1）的最小二乘解来得到展开系数。否则，该问题是欠定的（under－determined），并且在对（1）的无限个解之中，可以采用具有最小2范数的一个解。该解将被称为亥姆霍兹方程最小范数（HELN）解。

为了得到矩阵形式，我们在测量位置处定义波函数值的矩阵B：

（2）

并且所测量的压力和展开系数分别是矢量p和a：

（3）

使用这些等式，等式（1）可以被重写成下面的形式：

（4）

其中I ≥ N，最小二乘（HELS）解是

并且其中I < N，最小范数（HELN）解是

符号H表示厄密特转置（Hermitian transpose）。特别地，必须在这些解中应用正则化。尽管可以使用其他正则化方案，但是为了本描述，我们将仅考虑Tikhonov正则化，并且在这种情况下HELS和HELN公式分别变成；

（5）

以及

（6）

ε是正的正则化参数，并且I是适当维度的单位对角二次矩阵。注意，实际上两个正则化公式（5）和（6）确切地提供与值I和N无关的相同解。然而，对于I < N，使用HELN公式（6）在计算上仍有优点，因为矩阵 BB ^H具有比 B ^H B 小的维度。

当已知展开系数时，可以通过波函数贡献（contribution）的总和来计算Ω中的任意点r处的声压：

（7）

T表示矩阵或矢量的转置。然后可以通过应用欧拉等式从声场表达式（7）估计方向χ上的粒子速度：

（8）

其中，ω是角频率，ρ₀是介质密度，并且隐式时间变量是

。随后可以使用的其他符号是声音在介质中的传播速度c，并且波数是k ≡ ω/c。

平面波构成波数域中的连续谱，所以在上述公式中使用这些必须包括通过在波数域中进行采样而获得的波函数的子组。本文所述的方法实施例通过利用SONAH方法以双向方式解决场估计问题来克服该采样问题，所述SONAH方法例如J.Hald发表于Inter-Noise 2003论文集中的“Patch near-field acoustical holography using a new statistically optimal method”以及J. Hald发表于Euronoise 2006论文集中的 “Patch holography in cabin environments using a two-layer handheld array and an extended SONAH algorithm” 。

根据SONAH方法，得到在Ω中的点r处的未知声压p(r)以作为已知声压值的线性组合：

（9）

其中矢量 c (r)包含复估计权重。然后通过需要预测公式（9）在最小二乘意义上为所有展开函数Ψ_n提供最佳平均预测来得到矢量c。参考等式（2）和（7），这意味着 c (r)应该被看作下述线性等式组的最小二乘解：

（10）

每当HELS等式（4）是过定（over-determined）的时候，该等式是欠定的，并且反之亦然。为了本描述，将会考虑Tikhonov正则化解，意味着根据最小二乘公式（10）的最小二乘和最小范数解可以表达为：

（11）

其中矩阵A是B的转置：

（12）

根据（7）和（12）的定义，之后通过下面的公式（13）和（14）分别给出 A ^H A 和 A ^H α(r)的元素，“*”表示复共轭：

（13）

（14）

矩阵 A ^H A 可以被看作基波函数（展开函数）的域中的测量点之间的互关联矩阵，并且然后矢量 A ^H α(r)包含测量点和估计点r之间的互关联（参见例如 R. Steiner和J. Hald的“Near-field Acoustical Holography without the errors and limitations caused by the use of spatial DFT,” Intern. J. Acoust. Vib. 6, 83-89 (2001))。还要注意， A ^H A 是具有非负本征值的厄密特对称矩阵，并且其维度等于测量点的数目I，其与展开函数的数目无关。

上述预测公式将被称为SONAH预测公式，并且它提供与HELS/HELN公式相同的解。实际上，HELN等式（7）和（6）产生SONAH公式：

（15）

在此处，我们已使用了这样的事实：对于任何非奇异二次矩阵，逆的转置等于转置的逆。

因为SONAH不明确地计算展开系数a，所以不能对一组展开系数执行加权/滤波。这样的滤波被集成在权重矢量 c (r)的计算中。这样做的工具是正则化，但是使它适当工作需要定义矩阵 A ^H A 的基波函数的适合的缩放。典型地，通过选择更大值的正则化参数ε，我们将有可能抑制还没有被抑制的这些渐逝基波，它们需要在源表面上的它们的重建中的最大放大。为了做到这一点，这些基波必须对应于仍在Tikhonov正则化的截止级之上的 A ^H A 的最小奇异值。通过确保在重建中需要较大放大的波函数具有横跨传声器位置的较小振幅，可以以基波函数的缩放来支持这样的行为。

在本实施例中，基波是平面传播和渐逝波，并且我们首先考虑其中如在图3a中说明的在自由场情形下重建声场的情况，即其中在源平面107之后z=z ₀ , z ₀ <0处的一个或多个源辐射到无源匀质半空间z ≥ z ₀（被表示为Ω）中。在这种情况下，可以由如下形式的平面波函数来表示辐射的声场

（16）

其中F是振幅缩放函数，并且z=z ⁺, z ⁺≤z ₀是所谓的虚拟源平面101，其定义波函数，例如在垂直于虚拟源平面的方向上衰减的渐逝分量。虚拟源平面通常被定义为位于真实源平面内部的小距离处。源平面107是平行于虚拟源平面的最近平面，所述虚拟源平面仍在无源区中，因此，平行平面仅仅碰到（即正切于）真实源表面。在仅具有单个虚拟源平面的自由场情况下，xy平面通常平行于虚拟源平面。在等式（16）中， r≡ (x,y,z)和 k≡ (k _x ,k _y ,k _z )，其中k_x和k_y是实波数，并且k_z是如下定义的复波数：

（17）

因此，展开函数的完整组包括两个波数变量k_x和k_y的连续谱。虚拟源平面z=z ⁺可以与真实源平面在z=z ₀处重合，但是通常它将在该平面之后，即z ⁺ <z ₀ 。虚拟源平面的坐标上的上标“+”指示在正z轴方向上传播的波被表示。

注意，在等式（16）中，缩放函数F(k _z )定义虚拟源平面z=z ⁺上的所有波函数的振幅。

为了使用上述离散HELS或SONAH公式，我们最初通过在2D波数域(k _x ,k _y )中的正则采样来选择离散子组：

（18）

当选择采样间距时，应该意识到波数域中的正则离散表示将表示空间上周期性的场，导致正如在基于空间DFT处理的传统NAH中的环绕误差。因此，采样间距应该比2π除以最小可接受周期长度更小。同时，采样区域应该覆盖在测量区中具有显著振幅的所有这些渐逝波以及所有传播波。这通常导致展开函数的数目N大于测量位置的数目I。因此，等同的HELS问题将是欠定的，并且解（6）将主要是最小范数解。

在一些实施例中，可以使用波函数的单位加权，其中F(k _z )=1。然后，根据等式（16），所有平面波函数具有在虚拟源平面

等于一的振幅。HELS展开系数矢量a因此仅仅是该平面中声场表示的波数谱，并且矩阵B可以被看作从该平面波谱到传声器位置处的声压力值的传递函数的矩阵。由等式（16）给出B的元素，其中插入传声器坐标。有趣的是，通过转换矩阵B，渐逝波分量（具有|k|>k的Φ_k）将得到随着波数以指数增加的衰减。因此，整个趋势是这些波数分量对应于矩阵B的越来越小的奇异值。该性质确保可以在波数谱a的反解中由Tikhonov正则化来执行渐逝波分量的滤波。传播波的一些组合也将对应于小的奇异值，例如表示来自远离测量位置的虚拟源平面中的区域的声辐射的组合。这样的组合也将被正则化抑制。

因为上面提到的HELS问题的欠定性质，解将主要是最小范数解。正则化将会引入解矢量a的范数的进一步减小。因此，在良好表示传声器处的压力的约束下，解过程将最小化波数谱中的能量，并且因此也最小化平行平面和虚拟源平面中的压力分布中的能量。为此，在测量位置周围的体积内实现好的声场表示，而远离测量位置将看到增加的低估。因为SONAH和HELS提供相同的解（15），所以对于SONAH来说结论将是相同的，如还将下面描述的结果说明的那样。

图4说明k _y =0处的沿k _x波数轴的等式（18）的正则采样。在辐射圆内部，即对于

在(k _x ,k _y )处的特定采样对应于在具有由等式（17）给出的k_z的矢量(k _x ,k _y ,k _z )的方向上传播的平面波。在图4中，可以通过将采样位置k _x投影到所示出的半径为k的半圆上来找到该方向。在(k _x ,k _y )平面中采样点具有恒定密度的情况下，平面波的角密度将在接近于z方向上最高，并且当偏轴角θ朝向90度增加时降低。因此，当使用常数加权函数F时，则在对于SONAH的预测加权的最小二乘解（11）中在近轴向方向传播的平面波将比远离轴波得到更高的权重。考虑到(k _x ,k _y )平面上的无限小区域分段，它将表示方向性半球上的区域；

其大到原来的1/cos(θ) = k/|k _z |倍。因此可以通过选择F来得到恒定的方向性功率谱密度，以使得

恒定，例如

（19）

其中F ₀ 是常数。等式（13）、（14）、（16）和（18）示出被缩放 |F| ²的波数贡献的加和而得到的矩阵 A ^H A 和 A ^Hα(r)的元素的贡献。

这是为什么选择对功率密度的均衡而不是对振幅密度的均衡是有优势的原因。已发现还应用于辐射圆之外的等式（19）的“全方向加权”以提供平均来说比常数加权F(k _z )=1稍微更好的精度。因此，在一些实施例中，使用全方向加权，并且它还将被用于说明遍及本描述的方法。

为了完全移除由波数域中的离散表示引起的环绕问题，我们现在使(k _x ,k _y )平面中的采样间距成为零，这意味着等式（13）和（14）中的加和变成(k _x ,k _y ) 平面上的积分。等式（11）示出，如果正则化因子被相应地调整，则应用于 A ^H A 和 A ^Hα(r)的元素的公共缩放因子不会改变SONAH结果。因此，我们可以选择以下述方式来缩放极限积分表达式表示

（20）

（21）

矩阵元素的极限值（20）和（21）对应于具有无限维度的矢量α(r)和矩阵A。不管怎样A和α(r)将被用于几个后续推导中，其中它们被应该理解为极限情况

。

通过等式（16）、（17）和（19）的应用，等式（20）和（21）中的二维表达式可以被降低成一维表达式：首先，对波数域中的极坐标(K,Ψ)执行变换(k _x ,k _y )= (K cos(Ψ), K sin(Ψ))，以及然后，可以分析得出极角Ψ的积分。结果，根据下面两个变量的两个函数，可以表示等式（20）和（21）的积分：

（22）

（23）

其中J ₀是零阶贝塞耳函数。函数I _p表示波数域中来自辐射圆内部的贡献，其中在积分表达式中用K=k sin(θ)替代，并且函数II _p表示来自辐射圆之外的贡献，其中在积分表达式中用代替。这些变量代替已被引入以使得（22）和（23）中的积分更适用于高斯类型的数值积分。使用函数I _p和II _p ，对于矩阵元素的公式（20）和（21）可以被重新公式化如下：

（24）

（25）

其中 r _i ≡ (x _i ,y _i ,z _i )和 r _j ≡ (x _j ,y _j ,z _j )是两个传声器位置，压力估计位置是 r≡ (x,y,z)，

是xy平面（虚拟源平面）上的传声器位置i和j的投影之间的距离，并且是xy平面上的压力计算位置和传声器位置i的投影之间的距离。现在使用等式（9）和（11）来估计位置r处的压力。J. Gomes 和P.C.Hansen在Acoustics’08, 2008的论文集中发表的“A study on regularization parameter choice in Near-field Acoustical Holography”中研究了用于确定等式（11）的正则参数ε的确定的不同方法，并且发现L曲线方法与SONAH结合提供整体大部分稳定结果。已在广义交叉校验（GCV）处指出早期的非关联测量误差的研究。

为了得到粒子速度估计，我们将欧拉等式（8）应用于由等式（9）和（11）给出的SONAH压力估计。结果，我们得到：

（26）

其中：

（27）

可以通过将等式(25)插入到（27）中，之后对函数Ip和IIp应用表达式（22）和（23）来得到 A ^H β _χ (r)的显式。对于z分量，我们得到：

（28）

并且对于x分量和y分量;

(29)

新的函数I _z , II _z , I _xy和II _xy连同I _p和II _p一起被列出在下面的表格1中。表格1还包括这些函数在利用常数加权(F(k _z )=1)的基波函数的实施例中所具有的形式。

表格1:在对于不同加权方案的关联函数的计算的实施例中使用的2维函数。

在上面的表格中，J ₀和J ₁分别是零阶和一阶贝塞尔函数。发明人已能够找到仅对于积分I _p和I _xy的实部的分析表达式，然而，积分通常可以由数值积分来计算，例如使用任何合适的数值积分方法，比如高斯类型的算法。作为对上述积分中的三个中的积分到无穷大的代替，积分上限可以在指数函数已达到最够低的等级的点处截断，或者可以基于频率和测量网格的间距来选择上限。

如上文所提到的那样，通过（例如由数值积分）生成并存储对于变量ρ和ζ的离散值的2D表格，可以相当大地加速大量的积分。为了计算对于ρ和ζ的特定值的积分，于是必须执行表的内插。可替换地或附加地，可以构造覆盖不同域的2D近似函数。

积分I _p , I _xy和I _z没有奇异点，所以它们可以被直接放在表格中并且被内插。积分II _p , II _xy 和 II _z有奇异点，其中两个参数等于零，即在(ρ,ζ)=(0,0)设置的参数的原点处。为了使内插方案在接近该原点处适当地工作，应该从表格化的函数中提取出该奇异点。

采用波函数的方向加权的情况，于是对于非常小的ρ和ζ值，我们具有下式：

对于，

对于

，

对于

，

以下面的方法在内插期间可以移除这些奇异点：对于II _p的情况，通过对于存储在表格中的参数组合组的数值积分算出积分，并且在将值存储在表格中之前将它们乘以。因此，可以对缩放的积分执行内插。在内插之后，然后用该结果除以

。可以类似地处理其他积分。利用单位加权，在截断它们之前以及在内插期间执行类似的奇异点移除。

有趣的是，当在平行于所有虚拟源平面的各个平面中布置测量和计算点时，可以获得进一步的效率增加。在这种设置中，对于每个频率，仅可能出现第二参数ζ的一些值，因此使得1D表格的生成很有吸引力，每个表格对应于ζ的相应值。例如，利用不同于测量平面的单个虚拟源平面、单个平行测量平面和单个平行计算平面，可以设置2组1D表格。在执行期间可以对每个频率高效地完成这些。类似地，利用一个虚拟源平面、两个平行测量平面以及单独的平行计算平面，可以对每个频率设置4组1D表格。例如，1D表格可以具有通过内插预先计算的2D表格而生成的优点。

一般来说，如上所述，当已在测量位置处测量声压时，可以得到在目标位置r处的压力如下：

其中p是给定位置处的测量/已知的声压数据的矢量，并且ε是正则化参数。可以从类似的表达式获得方向χ上的粒子速度：

上述公式使用Tikhonov类型的正则化，但是也可以使用采用类似正则化方案的其它公式，例如所谓的增强型Tikhonov正则化。

现在参考图2并继续参考图1，将描述基于上述方法的重建过程的实施例。

在最初的初始化步骤S1中，过程定义合适的坐标系113和虚拟源平面101。

一般来说，在测试设置期间，可以将测量位置选择为位于沿着要研究的对象的真实源表面。例如，假设自由场情形和平面换能器阵列，可以选择近似平行于一部分真实源表面的测量平面。类似地，在其中使用3D阵列的非自由场的情况下，阵列被假设成相当平坦，并且阵列的中间平面被定位为接近平行于一部分真实源表面。

在步骤S1中，过程可以接收关于测量点的位置的信息，即过程例如通过上述合适的位置检测系统接收一组位置 r _i , i = 1,2,…,I的坐标。该过程然后可以相对于测量几何结构选择坐标系。例如，对于平面换能器阵列，该过程可以定义xy平面位于测量平面（即位于z＝0平面中的测量平面）中，其中原点在平面阵列的中心，并且正z轴指向远离源表面，即进入到无源域中。类似地，对于3D阵列，z＝0平面可以被定义成阵列的中间平面。然而，将会认识到可以以另一种方式来选择坐标系。

在初始化步骤S1中，该过程可以进一步相对于所选的坐标系定义虚拟源平面101，例如平行于xy平面。可以自动或者至少部分基于用户输入来执行虚拟源平面的定位。

在一些应用中，对过程来说真实源几何结构不是已知的或者至少不会准确地已知。在这种情况下，该过程可以例如在与离阵列预定距离处平行于该阵列的平面上重建声场。在这种情况下，阵列和平行虚拟源平面之间的距离可以被指定为对于该过程的输入参数，例如指定为用户定义的输入参数。

在其他情况下，真实源几何结构可以是已知的并且可以被输入和对齐，或者可以利用还用于测量阵列位置的位置测量系统来测量。如果阵列位置和对象几何结构二者都是已知的，则处理软件可以自动确定二者之间的距离。并且基于此，该过程可以自动地将虚拟源平面定位于例如真实源表面内的指定距离且与阵列平行，若假设该阵列是平面的话。如果阵列不是平面的，则可以定义中间平面。

一般来说，可以在对象内部的短距离定义虚拟源平面，其中所选的距离可以取决于测量几何结构，例如距相邻测量点之间的典型或平均距离。例如，虚拟源平面可以被选择成大于（例如在大约1倍到大约3倍之间）源对象内部的远离源表面的测量点之间的平均最近相邻距离，即在与测量点的位置和重建声场的（一个或多个）目标位置的相对的真实源表面的侧面。已发现，通过将虚拟源平面定位于大约真实源表面之后的一个和半个测量网格间距（即比物理源表面进一步远离测量网格）来获得最优精度。常常发现稍稍更大的距离可以提供低频率处的更好结果。

为了重建过程，虚拟源平面的一侧被看作无源和匀质的，所以源场可以由与该虚拟源平面相关联的一组平面传播和渐逝波表示。

步骤S1的初始化还可以包括一个或多个参数的设置，（例如应用的加权方案、应用的正则化方案、所使用的测量系统的性质和类型（例如换能器数目）、要重建的声场的目标位置的坐标、要计算的一个或多个声学性质、频率范围、要被内插的函数表格的部分的初始化（可以用执行该过程的软件来提供其他部分）。例如，这些和附加或可替换的参数可以是用户经由分析单元103的适合的用户接口定义的。

在步骤S2中，该过程接收在该组位置 r _i测量的相应复时间谐波声压p(r _i ) 。 r _i表示该位置相对于上述坐标系的坐标矢量，例如根据从位置检测设备110接收的输入而确定的换能器位置。为了本描述，将假设位置 r _i位于由匀质流体（例如空气或水）占用的无源区Ω中。

在上文所讨论的自由场情况下，所有源都在z＝0平面的一侧，声场原则上可以是基于例如在该平面中的测量而完全重建的，但是此处所描述的方法还可以结合其他任意的3D测量网格来应用。此外，如将在下文所描述的那样，该过程还可以被应用于具有多个源平面的情况。

可以根据以各种各样的方式接收的换能器信号来确定复时间谐波声压p(r _i )。例如，可以利用所有的换能器来完成单个同时记录，之后是为每个FFT线提供复（时间谐波）声场的所有信号的FFT。可替换地，可以执行从参考信号到所有阵列换能器的交叉谱的测量，以及参考信号自谱的测量。被自谱的平方根除的交叉谱然后为所应用的谱的每个频率线提供复（时间谐波）声场。图2的过程基于所测量的声压力执行在Ω中的声场的重建。可替换的实施例可以测量一个或多个可替换的或附加的声学量，例如声压梯度、粒子速度等等，或者其组合。

在步骤S3中，该过程使用等式（24）计算自关联矩阵

，即作为分别在和处评估的二维函数I _p 和II _p的线性组合。为此，该过程分别将ρ和ζ的各个值用作进入I _p和II _p预先计算的函数值的查找表112的索引。查找表可以被存储在存储器或分析单元的其他适合的存储介质中。因此，该过程可以获得所存储的最接近ρ和ζ的函数值，并且在ρ和ζ中执行二维线性内插，或者另一适合的内插方案（例如三次内插）。

将会认识到，可替换地，该过程可以计算转置矩阵 BB ^H或者可替换的另一适合的变换矩阵。

在步骤S4中，该过程计算矢量q如下

或者其转置

在此处，ε是如上所述的正则化参数。根据指定信噪比（SNR）来计算正则化参数，或者可以使用广义交叉验证（GCV）、L曲线方法或类似的参数估计过程基于噪声测量的数据来自动计算所述正则化参数。当指定SNR时，则可以获得作为乘以10^(SNR/10)的矩阵 BB ^H中的对角元素的平均值或最大值的参数ε。在自由场的情况下，当仅使用单个虚拟源平面时，所有对角元素都相等，所以可以使用单个对角元素的值来代替平均或最大值。

可以通过任何用于求解矩阵等式的适合数值技术（例如对矩阵求逆）来执行上述矩阵操作。

根据所计算的q值，然后该过程可以计算声压。将会认识到，对声压来说可替换地或附加地，该过程可以计算不同的声学量，例如上述粒子速度。可以从所计算的压力和粒子速度值来获得声强度。

因此，该过程可以计算一组目标点r中的每一个处的声压和/或粒子速度，其中声场被重建。在图1中目标点的两个示例由109指定。在许多情况下，有趣的是在被研究的对象的源表面上重建声场。例如，该过程可以计算规则网格的一组网格点中的每一个网格点处的声压。为了本描述，将为其计算重建声场值的点还被称为计算点。

该过程可以通过所有计算点而循环，并且计算期望的声学量。具体地说，在步骤S5中，该过程计算在当前计算位置r处的期望声学量。

例如，该过程可以计算（多个）矢量

和/或

并且计算压力作为

和/或速度作为

该过程以与结合上面的步骤S3中的

的计算相似的方式（即作为在输入参数ρ和ζ的各个值处分别评估的二维函数I _p , II _p 和I _z , II _z 或者I _xy , II _xy的线性组合）分别使用等式（25）和（28）来计算关联矢量 A ^H α(r)和 A ^H β _χ (r)。如上那样，输入参数的值（在此处函数要被评估）被用作进入各个查找表112的索引以通过分别在表格化值I _p , II _p 和I _z , II _z或者I _xy ,II _xy之间执行三次内插或其他类型的内插获得2D函数的预先计算的函数值。

在步骤S6中，该过程确定是否留下未处理的计算点。如果情况如此，则处理返回步骤S5，否则，处理进入步骤S7，在此，该组计算出的声音压力和/或颗粒速度得以输出，例如，作为图形表示、作为数值阵列和/或以任何其他适合的形式。可替换地或附加地，从声压和粒子速度容易获得声强度。

在上文中，已结合如在图3a中说明的情况中的自由场情形描述了实施例。然而，此处所述的方法的实施例还可以用于其他情况，例如当源和反射对象不被限制成半空间时。在下文中，将描述上述实施例的修改，其为非自由场情形提供平面波公式，因此考虑了其中在等式（16）中定义的平面波函数不再构成对于测量区中的声场的完整基础的情况。

我们首先考虑在图3b中说明的情况，其中在位置z ⁺和z ^－（z ⁺ ≤ z ≤ z ^－）处的两个平行平面101a和101b之间的均质和无源介质中进行测量。对于z ≤ z ⁺的有源声场的部分可以由等式（16）的基本函数表示。为了也表示对于z ≥ z ^－的有源声场分量，我们添加对应的平面波函数组，产生组合的组：

对于z ≤ z ⁺的源（30）

对于z ≥ z ^－的源（31）

其中仍使用在等式（17）中k _z的定义，并且对这两个组应用全向加权：

（32）

因此，添加虚拟源平面对应于在SONAH（或HELS）方法中添加一组平面波函数。等式（13）、（14）、（20）和（21）示出考虑了所组合的基本函数组（30）和（31）的SONAH估计过程可以通过在矩阵

和 A ^H α(r)的元素中添加来自两个组的贡献而获得。这样做的结果是：

（33）

（34）

正如自由场的情况，估计粒子速度所需的矢量 A ^H β _χ (r)可以被导出，得到下面的结果：

（35）

（36）

上述公式对测量位置的选择没有约束，但是平面阵列可能不足以获得精确的预测。这可以通过在测量平面中图像对称的源然后不能被区分的事实来说明，因为它们在每个传声器处产生相等的压力。典型地，可以使用包括两个相等的平行平面矩形矩阵的双层阵列。可以通过常数系数F _0,z+ 和F _0,z-来调整阵列任一侧上的具有源的波的权重。典型地，两个参数将被选为等于一，但是如果在阵列的一侧不存在明显的源，则可以通过将它们中的一个设置成零、由此向该公式提供单侧源分布的信息来实现更好的估计。

通过在图3c中的截面图说明在无限平面段（例如板或片）的两侧上具有源的情况扩展到在盒形体积的所有六个侧面上具有源的情况。在这种情况下，我们仅将贡献添加到来自属于盒的每一侧面的平面波谱的矩阵元素。在下面的表格2中给出对应的公式。

表格2：对于所有六个侧面上的源的情况的公式

可以观察到对于仅两个侧面上的源的公式考虑来自所有方向的平面传播波，因为每个侧面包括来自具有相等权重的半球的所有平面波。但是两个侧面不包括所有类型的渐逝波：在x方向和y方向上衰减的波正在消失。

下面的示例给出模拟测量，其说明通过考虑在xy平面中衰减的波可以实现精确的改进。

因此，使用如在表格1中所示的相同2D函数，图2的方法可以被应用于在测量位置的多个侧面上具有源的情况。仍可以获得组成自关联和互关联矩阵的元素的关联函数作为这些2D函数的函数值的线性组合。

固有估计误差

在下文中，将描述用于至少确定平均来说在计算点处多好地预测此处所公开的重建方法的平面波函数的值的估计的方法的实施例。

考虑第一自由场情形，由等式（16）给出该组平面传播和渐逝波函数。根据等式（7）中的定义，矢量 α (r)包含在计算点r处的所有这些基波的真实值。等式（18）、（12）和（2）示出矩阵A的一行包含测量位置上的波函数之一的值，并且根据等式（9）和（11）因此之后矢量

被定义如下：

（37）

其包含在计算点处的波函数的估计。该矢量的2范数是预测的平均振幅的度量：

（38）

将由小于真实范数

的预测范数

来表征平均低估。为了找到真实范数

，我们从等式（18）、（12）和（2）注意到矩阵A的一个列包含在传声器位置处的所有波函数的真实值，正如 α (r) 包含在计算点处的所有函数值那样。因此，可以仅通过用计算位置r代替传声器位置 r _i 而使用用于矩阵 A ^H A 的对角元素 [A ^H A] _ii 的表达式（24）来计算 α ^H α 值。

（39）

在此处，可以使用定义（22）和（23）来容易地验证对于两个积分的分析表达式。然后，固有平均振幅增益被定义为：

（40）

可以通过应用等式（37）和（38）来计算估计误差矢量

的2范数：

（41）

并且最后，固有平均估计误差级被定义为：

（42）

通过应用对于矩阵

和 A ^H α (r)的元素的扩展公式（例如对于两个侧面上的源的情况的等式（33）和（34））可以简单地将上面的描述扩展到非自由场情形。在该情况下，对于

的表达式（39）必须被扩展，正如（24）被扩展成（33）那样：

（43）

因为下面两个原因，针对源的特定配置而观察到的相对预测误差将不同于上面的固有误差：

1）特定声场的平面波振幅谱将不同于该方法优化的那个平面波振幅谱。对于自由场情形，该优化振幅谱在所假设的虚拟源平面

上具有形式 |F(k _z )| 。

2）对于具体测量的声场，关于整个声压计算预测误差，而不是如平面波分量上的RMS那样。

下面将说明这两个差异中的每一个差异的影响。

示例：预测误差的数值研究

在下文中，将简述此处描述的方法对于几个典型阵列几何结构并作为频率的函数的预测的有效区的信息。

将不考虑测量误差，并且也将不研究对于正则化参数的自动确定的方法（GCV，L曲线）。为了本示例，使用下面的公式（参见例如R. Steiner 和J. Hald的“Near-field Acoustical Holography without the errors and limitations caused by the use of spatial DFT,” Intern. J. Acoust. Vib. 6, 83-89 (2001))）根据特定的动态范围D（信噪比）来确定在等式（11）和（26）中使用的正则化参数ε：

（44）

其中，[A ^H A] _ii是矩阵

的对角元素。导出对于平面阵列和自由场情形的该关系，在这种情况下所有对角元素是相等的。此处，我们通常使用对角元素平均值。D的典型值处于15dB和30dB的范围中。在下文中，我们使用D＝30dB，除非另外特别指出。

示例A：测量平面中的精确压力表示区

我们以两个不同不规则平面阵列设计的阵列平面z＝0中的普通固有误差分布的研究开始。这两个都具有等于0.75m的直径，并且针对高达至少20kHz的波束形成应用来优化二者。一个是包括如轮子中的辐条那样布置的13个相等的6元素线阵列的辐轮阵列，参见图5。另一个是包括具有12个不规则但均匀分布在每个扇区中的传声器的7个相等角度扇区的扇形轮阵列（参见例如J. Hald, “Array designs optimized for both low-frequency NAH and high-frequency beamforming”, Proceedings of Inter-Noise 2004）。在阵列几何结构优化期间保持扇形轮阵列的均匀分布，以便结合SONAH应用获得良好的性能。因为阵列是平面的，所以我们假设自由场情形，并且仅使用单个虚拟源平面。通过计算每个元素的平均阵列面积来定义平均元素间距，并且对于两个阵列几何结构该平均元素间距接近于7cm。对于规则的矩阵阵列，该间距将支持高达2kHz左右的频率范围，并且为了避免高达该频率的空间混叠，应该选择不小于阵列网格间距的相隔距离。因为虚拟源平面应该被定位于大约在真实源表面之后的z ₀ =-7cm处的一个传声器间距，所以我们选择离阵列平面14cm的虚拟源平面，即在z ⁺ =-14cm 处。

在图5中的等高线图示出阵列平面中的声压预测的固有相对误差级 E(r) ，由圆圈指示传声器位置。白色区域表示低于－20dB的固有误差级E(r) ，并且灰色等高线示出以2dB步长的精确度降级到大于由黑色表示的－6dB的误差级。在低频率处，两个阵列可以提供跨越整个阵列区域的良好声压表示，但是在1kHz处辐轮阵列的传声器分布中的间隔已变得太大以致不能在间隔中提供良好的表示。扇形轮阵列提供全阵列区域上一直都高于2kHz的可接受的压力表示，对于辐轮阵列几何结构的情况显然并非如此。辐轮阵列对于波束成型应用来说非常有利，不过其支持结构更容易以大尺寸构建。

示例B：测量平面之外的精确压力表示的区

我们现在转向垂直于阵列平面的平面中的精确声压预测的区的研究。作为对设计用于波束成型的不规则阵列的代替，我们考虑规则的8×8元素单层（8×8×1）和双层（8×8×2）阵列，其中在x方向和y方向上元素间距是3cm，从而覆盖区域 0 ≤ x ≤ 21 cm, 0 ≤ y ≤ 21 cm。这两种类型的阵列在z＝0的xy平面上具有传声器层，并且双层阵列在z＝3cm处具有附加的相等层。我们查看y＝10.5cm的xz平面（其是阵列的对称平面）中的相对声压预测误差。正如对于不规则平面阵列那样，我们为单层规则阵列选择离阵列等于两倍阵列元素间距的距离处的虚拟源平面，即z ⁺ = -6 cm。对于双层阵列，我们选择到前虚拟源平面的相同距离和到后虚拟源平面等于19倍元素间距的距离，即z ⁺ = -6 cm, z ^- = 60 cm。

图6在左列中示出单层阵列（具有单个虚拟源平面）的固有误差级E(r) ，并且在右列中示出双层阵列（具有两个虚拟源平面）的固有误差级E(r) 。中间的列示出单层阵列对于特定声场（在(x,y,z) = (7.5,7.5,-6.0) cm处的单极点源的场）的声压预测的相对误差。映射区域在沿着x轴的两个方向上延伸超过阵列区域6cm，并且在z方向上从z = -3 cm延伸到等于4倍正z方向上的阵列元素间距的距离。清楚地，特定声场预测的误差分布更不规则，并且一般来说比普通固有误差更低。在特定区域获得比由固有误差级预测的显著小的误差。这可以是因为来自不同平面波分量的误差贡献彼此消除。固有误差级避免这样的消除效应，并且提供该方法的估计精确度的更普通的图片。与单层阵列相比，看到双层阵列明显延伸精确表示的区，尽管已通过添加额外的虚拟源平面来增加声场中的自由度。具有低预测误差的区在z方向上的范围相当有限的原因是SONAH假设与xy方向上的源分布无关。为了提供良好的远场估计，可以将关于源区域的大小的一些假设集成到该方法中。这例如通过使用基于单极源的有限分布而非平面传播和渐逝波的SONAH公式来完成。

示例C：对改变参数的预测精确度的敏感性

为了限制对参数改变敏感性的研究的数据量，我们查看仅沿着垂直于阵列的中心的预测误差。将仅考虑示例B的单层8×8×1阵列。图7a示出与在图6的左列中表示的情形相同的情形的固有误差级，所以我们看到的仅是通过图6的左列中的等高线图的中间的一组水平切片。几个频率被添加到图7a中，尽管例如以3cm间隔的规则阵列网格的较高极限频率是5kHz。清楚地，在5kHz处，固有误差级已开始快速增加，甚至在测量平面中，但是仍可以实现可接受的预测。图7b示出振幅增益 A(r)的对应曲线。甚至在5kH处，在离阵列平面一个阵列网格间距内仍可以精确地预测振幅，但是在该间隔之外，观察到增加的低估。

在图7中，根据等式（44）确定正则化参数，其中固定的动态范围D等于30dB。下一步骤是看当动态范围D改变时预测误差如何改变。图8中的左列表示在低频（500Hz）和高频（4kHz）处固有误差级E对D中的改变的敏感性。对于每个频率，存在一组用对应的D值标记的曲线，所述D从10dB变化到40dB。一般来说，对比测量平面更远离源区的压力的预测不敏感，而更接近源区的预测用增加的D值（对应于更小的正则化参数值）得到了改进，尤其在低频处。这可以被预测，因为较小的D值将限制渐逝波的重建，并且因为渐逝波在低频且接近源区处构成更大部分的声场。图8中的右列给出在SONAH估计中固有振幅增益A的对应曲线。在500Hz处，除了在非常接近源区之外，振幅预测对D的改变不是非常敏感，而在4kHz处，当D的值变成比20dB左右小时在所有位置（甚至更远离源区）处都观察到增加的低估。解释必须是，在高频处传播平面波构成较大部分的平面波谱，并且因此构成较大部分的SONAH自关联矩阵

的奇异值谱。当D降低超过某一点时（正则化参数增加），则甚至部分传播波分量将被正则化明显衰减。

如上文所定义的固有振幅增益A和误差级E表示在所有平面传播和渐逝波上平均的估计位置的误差，其中如在SONAH估计算法的导出中所使用的那样将相同的权重应用于渐逝波分量。在这两种情况下，在（一个或多个）相同的虚拟源平面（对于自由场情形的

）中缩放波函数。现在要研究的问题是如果在SONAH计算中使用的虚拟源距离不代表在测量应用中渐逝波的实际内容，则误差级如何改变。应该考虑仅自由场情形和单个虚拟源平面的使用。为了完成该研究，我们注意到在等式（38）和（41）中关联矩阵

和

表示我们在其上平均估计误差的平面波场，而估计权重矢量 c (r)表示SONAH估计。因此，我们可以通过使用不同虚拟源距离来研究SONAH估计对1）测量中的渐逝波的实际内容和2）他们在SONAH计算中的权重之间的失配程度的变化的敏感性以用于1）在（38）和（41）中

和

的计算，和2）SONAH估计权重矢量 c (r)的计算。对于在图9中给出的结果，我们保留先前针对实际声场而使用的虚拟源平面z = z ⁺ = -6cm，即在等式（38）和（41）中使用的

和

的计算，同时变化的虚拟源平面距离d已被用于SONAH处理。大的d值将在SONAH估计中将低权重施加于渐逝波。图9中的左例给出在低频（500Hz）和高频（4kHz）处固有误差级E对d的改变的灵敏性。对于每个频率，存在一组标记有对应d值的曲线，所述d值从3cm变到12cm。右列包含振幅增益的对应曲线。在所有情况下，通过等于30dB的固定动态范围D指定正则化参数。

在500Hz处，误差级对太小的虚拟源距离（即重建中的渐逝波的太高权重）相当不敏感，但是太大的距离最终会导致比测量平面更远离源区的等级的稍稍过估计（A>1)。解释是波数泄漏（leakage），使得小的渐逝波分量如传播分量那样被处理，并且因此有助于场进一步远离：因为更小数目的测量点，传播和渐逝波分量之间的区别不是非常明显。当（由权重矢量 c (r) 表示的）SONAH处理将非常低的权重施加于渐逝波时，则渐逝波分量到传播分量的泄露会增加。查看z<0的误差增加（还在d为大时出现），解释是不同的：因为在SONAH计算矩阵中施加于渐逝波的权重太低，所选择的正则化参数将防止这些波分量中的某些被充分地重建。

在4kHz处，渐逝波通常构成非常小部分声场，所以对在SONAH计算中将或多或少的权重施加于他们的敏感性更小。最大的影响是当使用太小的虚拟源距离时对振幅增益的影响，从而在重建中施加于渐逝波的权重比实际平面波谱所建议的更大。靠近源区观察到高估，而远离时看到比正常情况更严重的低估。解释与针对500Hz处的高估给出的相同，泄漏现在仅以相反的方式进行：一些小的传播波分量将如渐逝波那样被处理，并且因此远离源区域时他们将衰减并且朝向源区域时将被放大。

示例D：阵列平面中的源-使用6个虚拟源平面

图10示出通过添加来自所有六个侧面而不是来自仅两个侧面的平面波函数以允许对沿着xy平面中的方向衰减的渐逝波分量正确建模而获得的精确度改进的示例。利用与结合示例B所描述的阵列几何结构相同的8×8×2双层阵列几何结构来执行模拟测量。单个单极点源被置于(x,y,z) = (-6.0,7.5,1.5) cm 处，即在两个阵列层之间、阵列区域之外的6cm的中心平面中。清楚地，该单极将创建沿着xy平面的方向衰减的渐逝波类型。SONAH被用于基于128个测量的压力值来预测在前阵列平面z = 0 cm中的测量位置处的强度矢量。然后使用下面的公式来计算相对平均误差级：

（45）

其中

是在计算点数目i处的真实强度矢量，

是相同矢量的SONAH估计，并且加和遍及i。使用30dB动态范围（D＝30dB）来执行SONAH计算，并且虚拟源平面被选成在所有方向上远离阵列6cm：

,

,

,

,,

。

因此，单极位于由6个虚拟源平面限定的盒子的边界上。查看图10，看出当存在所有6个虚拟源平面时低频强度估计的精确度比当仅存在与xy平面平行的虚拟源平面时明显更精确。然而，改进会随着增加的频率而消失，因为渐逝波在较高频率处构成较小部分的声场。

尽管已详细描述并示出了一些实施例，但是本发明不限于此，而是在下述权利要求限定的主题的范围内本发明可以以其他方式体现。例如，主要结合其中测量的声学量是声压的实施例来描述此处所公开的方法实施例。然而，将会认识到在其他实施例中，可以测量不同声学量。例如，作为对声压的代替，可以测量测量位置r _i处χ方向上的粒子速度u _χ ( r _i )。在这样的实施例中，作为对等式（1）的代替，当用于求解展开系数矢量a的线性等式的系统被设置时可以使用针对位置r _i的下述等式：

其中我们已引入展开函数的空间导数

这意味着在矩阵A和B中，函数Ψ_n(r _i)被Γ_χ,n(r _i) 代替（n=1,2,…,N），并且在等式（3）的矢量p中，压力p(r _i) 被所测量的粒子速度u _χ(r _i)代替。然后相应地必须改变自关联矩阵所包括的元素。

在此描述的方法和装置可以用于例如当分析机器、电机、引擎、交通工具（例如车辆）和/或等的声学特性时重建各种声音/噪声源（例如振动对象）的声音场。

在此描述的方法的实施例可以通过包括若干独特元件的硬件、并且/或者至少部分地通过合适编程的微处理器而得以实现。

在列举了若干部件的装置权利要求中，这些部件中的几个可以由一个以及相同元件、组件或硬件项来实施。在相互不同的从属权利要求中陈述或者在不同实施例中描述特定措施的仅仅事实并非表示这些措施的组合不能被有利地使用。

应强调，当在该说明书中使用术语“包括/包含”时，其指定所声明的特征、元件、步骤或组件的存在，但不排除一个或多个其它特征、元件、步骤、组件或其群组的存在或增加。

Claims

1.一种重建声场的方法，所述声场由至少一个声源生成，该方法包括：

－接收在一组测量位置中的各个位置测量的至少第一声学量的一组测量值；

－计算一组关联函数中的每一个，每个关联函数指示在所述测量位置中的第一位置处的至少一组平面波与第二位置处的至少一组平面波的关联；

－至少根据所计算的关联函数和该组测量值来计算指示目标位置处的重建声场的第二声学量；

其中该方法还包括存储用于内插对于关联函数的各个贡献的一组内插函数的表示，每个贡献是两个或更少输入参数的函数；并且其中计算该组关联函数中的每一个包括计算作为从所存储的表示获得的值的线性组合的每个关联函数。

2.根据权利要求1所述的方法，其中该组平面波是一组平面传播和渐逝波。

3.根据权利要求2所述的方法，其中该组平面波是平面传播和渐逝波的连续谱。

4.根据前述任一项权利要求所述的方法，其中第二位置是从目标位置和测量位置选择的。

5.根据前述任一项权利要求所述的方法，其中每个内插函数的表示包括查找表，每个查找表包括各个贡献的预先计算的值。

6.根据权利要求5所述的方法，其中每个查找表指示由两个输入参数编索引的函数值的二维阵列，并且其中每个函数值指示由所述两个输入值定义的输入点处的二维函数的值。

7.根据前述任一项权利要求所述的方法，其中至少一组平面波中的每一个平面波具有相对于预定虚拟源平面而限定的传播方向。

8.根据权利要求7所述的方法，其中包括至少一个声源的区由将至少一个声源从无源区分开的源表面限定，其中虚拟源平面具有指向无源区的第一侧，并且其中平面波点的传播方向远离虚拟声源并指向无源区。

9.根据权利要求7或8所述的方法，其中虚拟源平面被限定成位于接近至少一个声源的源表面的一侧上。

10.根据权利要求9所述的方法，其中虚拟源平面被限定为具有距源表面的、比该组测量位置的最近相邻测量位置之间的平均距离更大的距离。

11.根据权利要求7至10中任一项所述的方法，其中对于关联函数的各个贡献中的每一个是包括被两个输入参数参数化的积分的二维积分函数。

12.根据权利要求11所述的方法，其中两个输入参数中的第一个输入参数指示投影到虚拟源平面中的第一测量位置和第二位置之间的投影距离，所述投影距离由指示波数的因子缩放。

13.根据权利要求11至12中任一项所述的方法，其中两个输入参数中的第二个输入参数是指示从虚拟源平面到第一测量位置和到第二位置处的相应距离的线性组合的参数，其被指示波数的因子缩放。

14.根据权利要求7至13中任一项所述的方法，其中由根据预定加权方案的相应加权因子加权每个平面波。

15.根据权利要求14所述的方法，其中从均匀加权和方向性加权中选择预定加权方案，平面波的方向性加权是平面波的传播方向的函数。

16.根据权利要求14或15所述的方法，其中一组内插函数的表示包括用于关于每个加权方案内插相应六个二维函数的六个内插函数的表示。

17.根据权利要求7至16中任一项所述的方法，其中该组测量位置被布置在一个或多个平行测量平面中，所述测量平面与一个或两个虚拟源平面平行。

18.根据权利要求17所述的方法，其中对于关联函数的各个贡献中的每一个是被第一和第二输入参数参数化的二维函数；其中测量平面与虚拟源平面平行；其中目标位置被布置在与测量平面平行的重建平面中；并且其中存储一组内插函数的表示包括针对第二输入参数的各个值生成查找表，每个查找表指示由第一输入参数编索引的二维函数的函数值的一维阵列，所述阵列包括在第二输入参数的对应值以及第一输入参数的各个值处的二维函数的函数值。

19.根据权利要求18所述的方法，其中第一输入参数指示用于为其计算关联函数的两个位置的投影之间的距离，所述投影是虚拟源平面上的投影，并且所述距离乘以波数。

20.根据权利要求7至19中任一项所述的方法，其中第二声学量是在平行于虚拟源平面的方向上具有分量的粒子速度。

21.根据权利要求7至20中任一项所述的方法，包括限定多于一个虚拟源平面，每个对应于相应组的平面波；并且其中每一个关联函数指示在所述测量位置的第一个测量位置处的该组平面波与第二位置处的该组平面波的关联。

22.根据前述任一项权利要求所述的方法，包括使用该组关联函数来计算指示所计算的第一声学量的估计误差的量。

23.根据前述任一项权利要求所述的方法，其中以不规则的平面阵列布置测量位置。

24.根据前述任一项权利要求所述的方法，其中从声压、声压梯度和粒子速度选择第一声学量。

25.根据前述任一项权利要求所述的方法，其中从声压、粒子速度和声强度选择第二声学量。

26.根据前述任一项权利要求所述的方法，其中至少一个内插函数的表示包括查找表，所述查找表包括各个贡献的预先计算的函数值；其中所述贡献是包括对于至少一对输入值的奇异点的函数；并且其中计算作为从所存储的表示获得的值的线性组合的关联函数包括：

－获得多个乘以相应缩放因子的所存储的函数值以补偿奇异点，所述缩放因子取决于对于贡献的各个输入值；

－执行所缩放的函数值之间的内插以获得缩放的内插值，并且内插值除以对应的缩放因子。

27.一种用于重建声场的处理装置，处理设备包括用于接收在一组测量位置中的各个位置测量的至少第一声学量的一组测量值的接口；并且处理单元被配置成：

其中所述处理单元包括用于存储用于内插对于关联函数的各个贡献的一组内插函数的表示的存储介质，每个贡献是两个或更少输入参数的函数；并且其中所述处理单元适于计算作为从所存储的表示获得的值的线性组合的该组关联函数中的每一个关联函数。

28.一种用于重建声场的系统，该系统包括：在权利要求27中限定的处理装置；以及用于在一组测量位置处测量第一声学量并且能够以通信连接的方式连接到所述装置以便将所测量的第一声学量转发到所述处理装置的一组换能器。

29.将根据权利要求27的装置用于重建由交通工具的噪声源生成的声场。

30.包括程序代码装置的计算机程序，当在数据处理系统上执行所述程序代码装置时，适于使得数据处理系统执行根据权利要求1至26中任一项所述的方法的步骤。