CN106772220A - 基于封闭空间几何信息建模的单传声器声源定位方法 - Google Patents

Abstract

本发明公开了一种基于封闭空间几何信息建模的单传声器声源定位方法，用于解决现有封闭空间声源定位方法实用性差的技术问题。技术方案是在封闭空间边界所包围的流体区域内布置n个节点，将此n个节点从1到n进行编号。利用此n个节点构建出封闭空间及内部物体的几何形状，并根据此n个节点坐标信息建立描述封闭空间数学性质的系统矩阵；在实际声源定位时，利用单传声器测得多个频率下的声信号，将此声信号与系统矩阵共同运算，得到声信号的位置信息，从而实现声源定位，实用性好。

Description

基于封闭空间几何信息建模的单传声器声源定位方法

技术领域

本发明涉及一种封闭空间声源定位方法，特别涉及一种基于封闭空间几何信息建模的单传声器声源定位方法。

背景技术

封闭空间内的声源定位问题是声学研究中的一项重要内容，在人们生活中的很多方面都有较广的应用。例如飞机舱室的噪声源定位问题，会议室的智能话筒传音控制系统，以及地下车库的车位识别系统等。这类空间的主要特点是环境复杂且混响程度较高，某些小尺度的空间还有严重的声波干涉、衍射效应，对声源定位技术的精度、稳定性要求较高。

文献“基于传声器阵列的声源定位，信息技术，2016(6):136-138”采用了基于声达时间差的定位技术，这种声源定位方法一般分为2个步骤进行：先进行声达时间差估计，并从中获取传声器阵列中阵元间的声延迟；再利用获取的声达时间差，结合已知的传声器阵列的空间位置进一步定出声源的方向。这种方法相对于可控波束形成技术和基于高分辨率的谱估计技术的精度较好，计算量较小。

从算法和实际应用来看，这种声源定位技术还存在着一些缺陷。首先，要获得较高的定位精度，需要大规模的传声器阵列，这在经济上会造成很大的负担；对于基于声达时间差的定位技术，声源定位实际上只是定位出了声源的方向，具体的距离则无法给出，另外，这种方法受环境混响影响很大，当混响程度较高时，精度很低。同时该方法在定位阶段用的参数已经是对过去时间的估计，因此不是最优估计。

发明内容

为了克服现有封闭空间声源定位方法实用性差的不足，本发明提供一种基于封闭空间几何信息建模的单传声器声源定位方法。该方法在封闭空间边界所包围的流体区域内布置n个节点，将此n个节点从1到n进行编号。利用此n个节点构建出封闭空间及内部物体的几何形状，并根据此n个节点坐标信息建立描述封闭空间数学性质的系统矩阵；在实际声源定位时，利用单传声器测得多个频率下的声信号，将此声信号与系统矩阵共同运算，得到声信号的位置信息，从而实现声源定位，实用性好。

本发明解决其技术问题所采用的技术方案：一种基于封闭空间几何信息建模的单传声器声源定位方法，其特点是包括以下步骤：

步骤一、在封闭空间边界所包围的流体区域内布置n个节点，将此n个节点从1到n进行编号。n个节点一方面作为封闭空间的几何信息在后期进行建模计算，一方面用来表示未来声源定位时的位置。

n个节点近似均匀分布于整个流体区域内，封闭空间的边界处需布有节点。当封闭空间内部存在物体时，物体内部不布置节点，且物体边界上应布置有节点，即节点应勾勒出物体形状。节点布置完毕后，封闭空间及其内部物体的边界用Γ表示，边界所包围的流体区域用Ω表示。

步骤二、假定在封闭空间内声源的位置为r，它在单位时间内向单位体积内的空间提供了ρ₀q(r,t)的媒质质量。根据质量守恒定律，媒质中声波的连续方程写为：

式中，ρ'为媒质密度增量，ρ₀表示媒质静态密度，q为q(r,t)的简写，v为媒质质点速度，t表示时间，div为散度算子，在三维空间笛卡尔坐标系中，

除了连续性方程之外，用来描述媒质声波的基本方程还有两个，它们不受声源的影响，分别为运动方程：

和物态方程：

p＝c₀ ²ρ' (3)以上两式中，p代表声压，c₀代表声速，grad为梯度算子，在三维空间笛卡尔坐标系中，

由媒质中声波的三个基本方程得到有源情况下封闭空间中有关声压p的波动方程：

式中，为拉普拉斯算子，在三维空间笛卡尔坐标系中，

在频域，声源强度q(r,t)表示为

q(r,t)＝q_ω(r)e^jωt (5)

式中，ω为谐振频率，q_ω(r)为在位置r处频域内的声源强度。

封闭空间声场各点的声压的频率与声源相同，声压表示为：

p(r,t)＝p_ω(r)e^jωt (6)

式中，p_ω(r)为在位置r处频域内的声压。

将式(6)及式(7)代入式(5)中，得到简谐声源激励下的声波波动方程为：

令式中k＝ω/c₀，称其为波数，并且消去e^jωt，得到只依赖于空间坐标的那部分方程，即室内有源Helmholtz方程：

这样就将声压的时域问题转换为频域问题，式(9)即为封闭空间声场的控制方程。

在封闭空间中，边界具有吸声能力，其声压梯度表示为：

式中，n为封闭空间壁面外法线方向，ζ称为比声阻抗，满足下式：

式中，Z为界面声阻抗。

根据Galerkin型加权残量法，为了求解式(8)，首先设一试函数为代入有源Helmholtz方程及其边界条件，试函数产生残量R和

根据伽辽金法确定权函数，有

由格林第一公式

式(13)简化为

在声场中任意一点的声压用各节点声压来表示，即

式中，N_i为节点i处的形函数，p_i为节点i处的声压。

将式(16)代入式(15)，得

式中，为形函数的导数矩阵，其表达式为：

整理式(17)，得到

令

其中，K称为刚度矩阵，M称为质量矩阵，C称为阻尼矩阵，G称为载荷矩阵。当声源位于位置r₀(x₀,y₀,z₀)处时，频域内的声源强度表示为：

q_ω(r)＝q_ωδ(r-r₀) (24)

其中

将式(24)代入式(19)中，得

最后，将式(20)、(21)、(22)、(23)代入式(19)并整理得到

式中，F＝jρ₀N^Tq_ω。K、C、M均为n×n阶的系数矩阵，各自的表达式分别为：M＝∫_ΩN^TNdΩ/c₀ ²，C＝∫_ΓN^TNdΓ/c₀ζ，各式中的N为形函数，在实际求解中，K、C、M中的积分运算用求和运算代替 其中m为积分点的数量，m_b为边界上积分点的数量，ξ_i为积分系数；ω为圆频率；c₀为空气中的声速；ζ称为比声阻抗，满足ζ＝Z/ρ₀c₀，ρ₀为空气密度，Z为界面材料的声阻抗；p为封闭空间内任意位置处的声压，实际定位时为单传声器所测得的声信号；F为n×1阶的列向量，表示声源相关信息，其表达式为F＝jρ₀N^Tq_ω，其中q_ω表示声源强度。

步骤三、根据节点坐标，利用移动最小二乘法，获得步骤二中所涉及的形函数N。利用移动最小二乘法构建形函数。一个场函数u(x)在一点的近似值表示为：

其中是计算点x的邻域范围内各节点的坐标，为基函数向量，m为基函数的个数，a(x)＝[a₁(x),a₂(x),…a_m(x)]为待定系数向量。使用单项式基函数做运算，在三维空间中常用的线性及二次单项式基函数分别为：

将求解域用节点离散后，在每个节点处定义一个权函数该函数只在支撑域内不为零，在支撑域之外为零，在三维情况下，权函数的支撑域为球形。选定权函数后，就求得近似函数在节点处的误差加权平方和：

令J取最小值，即

经过整理后，得到下式：

A(x)a(x)＝B(x)u (32)

式中，A(x)，B(x)的含义为：

由公式(32)得到a(x)，将其代入式(28)得：

步骤四、在封闭空间内部任意位置处设置一传声器。当声源有声音发出时，传声器拾取到一段音频信号f(t)，将其进行傅里叶变换，得到此音频信号的频域信号F(ω)。

步骤五、将步骤四中得到的频域信号F(ω)作为式(29)中的p值代入式(29)，求解后得到n×1阶的列向量F，由于F中的j、ρ₀、q_ω均为常数，因此，求得的列向量中值最大的元素所代表的节点的位置即为声源的位置。

本发明的有益效果是：该方法在封闭空间边界所包围的流体区域内布置n个节点，将此n个节点从1到n进行编号。利用此n个节点构建出封闭空间及内部物体的几何形状，并根据此n个节点坐标信息建立描述封闭空间数学性质的系统矩阵；在实际声源定位时，利用单传声器测得多个频率下的声信号，将此声信号与系统矩阵共同运算，得到声信号的位置信息，从而实现声源定位，实用性好。本发明将几何信息建模引入到封闭空间的声源定位之中，仅需单传声器即可完成精度较高的定位。解决了背景技术声源定位方法需要使用传声器阵列进行声源定位的技术问题，仅仅利用单传声器测得多个频率下的声信号，将此声信号与系统矩阵共同运算，得到声信号的位置信息，从而实现快速实时声源定位。

下面结合附图和具体实施方式对本发明作详细说明。

附图说明

图1是本发明基于封闭空间几何信息建模的单传声器声源定位方法的流程图。

图2是本发明方法所涉及的封闭空间及传声器、声源示意图。

图3是本发明方法所涉及的声源点及接收点设置示意图。

具体实施方式

参照图1-3。本发明基于封闭空间几何信息建模的单传声器声源定位方法具体步骤如下：

实际封闭空间为一矩形空间，其长宽高分别为0.9m、1.0m、1.2m，其边界为亚克力玻璃结构。坐标原点设置在封闭空间的底角处，传声器吊装在封闭空间顶部，其坐标为(0.5m，0.5m，1m)，声源位置坐标为(0.7m，0.8m，0.3m)。

步骤1、建立包含封闭空间几何信息的节点模型，由于此空间形状规则，因此节点为均匀分布，每个坐标轴向上有9个节点，共有729个节点。

步骤2、获取计算声源位置的系统方程。

步骤3、计算形函数N。

步骤4、在(0.5m，0.5m，1m)处吊装单传声器。当声源有声音发出时，传声器拾取到一段音频信号f(t)，将其进行傅里叶变换，得到此音频信号的频域转换信号，为F(ω)。

步骤5、将步骤4中得到的频域信号F(ω)作为式(29)中的p值代入式(29)，求解后得到n×1阶的列向量F，由于F中的j、ρ₀、q_ω均为常数，列向量F实际由N控制。因此，求得的列向量F中值最大的元素所代表的节点的位置即为声源的位置。

在本实施例中，对10个声源位置分别利用不同的声音作为测试信号进行了定位测试，每个声源位置测试10次。最终测试正确率为96％，有效地证明了本发明的有效性。

更具体的方法步骤如下：

步骤一、设有一封闭空间，在其内部要实现声源定位。在封闭空间边界所包围的流体区域内布置n个节点，将此n个节点从1到n进行编号。n个节点一方面作为封闭空间的几何信息在后期进行建模计算，一方面用来表示未来声源定位时的位置。

节点布置方式没有严格规定，但是应近似均匀分布于整个区域内，封闭空间的边界处需布有节点。当封闭空间内部存在物体时，物体内部不布置节点，且物体边界上应布置有节点，即节点应勾勒出物体形状。节点布置完毕后，封闭空间及其内部物体的边界用Γ表示，边界所包围的流体区域用Ω表示。

步骤二、获取计算声源位置的系统方程。假定在封闭空间内声源的位置为r，它在单位时间内向单位体积内的空间提供了ρ₀q(r,t)的媒质质量。根据质量守恒定律，媒质中声波的连续方程写为：

式中ρ'为媒质密度增量，ρ₀表示媒质静态密度，q为q(r,t)的简写，v为媒质质点速度，t表示时间，div为散度算子，在三维空间笛卡尔坐标系中，

除了连续性方程之外，用来描述媒质声波的基本方程还有两个，它们不受声源的影响，分别为运动方程：

和物态方程：

p＝c₀ ²ρ' (3)

上两式中，p代表声压，c₀代表声速，grad为梯度算子，在三维空间笛卡尔坐标系中，

利用与推导无源波动方程相类似的方法，由媒质中声波的三个基本方程得到有源情况下封闭空间中有关声压p的波动方程:

式中为拉普拉斯算子，在三维空间笛卡尔坐标系中，

在频域，声源强度q(r,t)表示为

q(r,t)＝q_ω(r)e^jωt (5)

式中ω为谐振频率，q_ω(r)为在位置r处频域内的声源强度。

由于通常将封闭空间声场作为线性系统考虑，因此空间内部各点的声压的频率与声源相同，声压表示为：

p(r,t)＝p_ω(r)e^jωt (6)

式中p_ω(r)为在位置r处频域内的声压。

将式(6)及式(7)代入式(5)中，得到简谐声源激励下的声波波动方程为：

令式中k＝ω/c₀，称其为波数，并且消去e^jωt，得到只依赖于空间坐标的那部分方程，即室内有源Helmholtz方程：

这样就将声压的时域问题转换为频域问题，式(9)即为封闭空间声场的控制方程。

在封闭空间中，边界通常具有一定的吸声能力，其声压梯度表示为：

式中n为封闭空间壁面外法线方向，ζ称为比声阻抗，满足下式：

式中Z为界面声阻抗，查表得到。

根据Galerkin型加权残量法，为了求解式(8)，首先设一试函数为代入有源Helmholtz方程及其边界条件，由于试函数通常不是精确解，因此将产生残量R和

根据伽辽金法确定权函数，有

由格林第一公式

式(13)简化为

在声场中任意一点的声压用各节点声压来表示，即

式中N_i为节点i处的形函数，p_i为节点i处的声压。

将式(16)代入式(15)，得

式中为形函数的导数矩阵，其表达式为：

整理式(17)，得到

令

∫_Ωjρ₀ωN^Tq_ωdΩ＝G (23)

其中K称为刚度矩阵，M称为质量矩阵，C称为阻尼矩阵，G称为载荷矩阵。当声源位于某一特定位置r₀(x₀,y₀,z₀)处时，频域内的声源强度表示为：

q_ω(r)＝q_ωδ(r-r₀) (24)

其中

将式(24)代入式(19)中，得

G＝∫_Ω-jρ₀ωq_ωδ(r-r₀)N^Tdv＝-jρ₀ωq_ωN^T (26)

最后，将式(20)、(21)、(22)、(23)代入式(19)并整理得到

式中F＝jρ₀N^Tq_ω。

式中K、C、M均为n×n阶的系数矩阵，各自的表达式分别为：M＝∫_ΩN^TNdΩ/c₀ ²，C＝∫_ΓN^TNdΓ/c₀ζ，各式中的N为形函数，在实际求解中，K、C、M中的积分运算用求和运算代替 其中m为积分点的数量，m_b为边界上积分点的数量，ξ_i为积分系数；ω为圆频率；c₀为空气中的声速；ζ称为比声阻抗，满足ζ＝Z/ρ₀c₀，ρ₀为空气密度，Z为界面材料的声阻抗，查表得到；p为封闭空间内任意位置处的声压，实际定位时为单传声器所测得的声信号；F为n×1阶的列向量，表示声源相关信息，其表达式为F＝jρ₀N^Tq_ω，式中q_ω表示声源强度。

步骤三、计算形函数。根据节点坐标，利用移动最小二乘法，获得步骤2中所涉及的形函数N。利用移动最小二乘法构建形函数。一个场函数u(x)在一点的近似值表示为：

其中是计算点x的邻域范围内各节点的坐标，为基函数向量，m为基函数的个数，a(x)＝[a₁(x),a₂(x),…a_m(x)]为待定系数向量。通常使用单项式基函数做运算，在三维空间中常用的线性及二次单项式基函数分别为：

将求解域用节点离散后，在每个节点处定义一个权函数该函数只在一个有限区域(支撑域)内不为零，在区域之外为零，在三维情况下，权函数的支撑域通常为球形。常用的权函数有高斯函数、样条函数等。选定权函数后，就求得近似函数在节点处的误差加权平方和：

令J取最小值，即

经过整理后，得到下式：

A(x)a(x)＝B(x)u (32)

式中A(x)，B(x)的含义为：

由公式(32)得到a(x)，将其代入式(28)得：

步骤四、在封闭空间内部任意位置处设置一传声器。当声源有声音发出时，传声器拾取到一段音频信号f(t)，将其进行傅里叶变换，得到此音频信号的频域转换信号，为F(ω)。

步骤五、将步骤四中得到的频域信号F(ω)作为式(29)中的p值代入式(29)，求解后得到n×1阶的列向量F，由于F中的j、ρ₀、q_ω均为常数，因此，求得的列向量中值最大的元素所代表的节点的位置即为声源的位置。

Claims (1)

1.一种基于封闭空间几何信息建模的单传声器声源定位方法，其特征在于包括以下步骤：

步骤一、在封闭空间边界所包围的流体区域内布置n个节点，将此n个节点从1到n进行编号；n个节点一方面作为封闭空间的几何信息在后期进行建模计算，一方面用来表示未来声源定位时的位置；

n个节点近似均匀分布于整个流体区域内，封闭空间的边界处需布有节点；当封闭空间内部存在物体时，物体内部不布置节点，且物体边界上应布置有节点，即节点应勾勒出物体形状；节点布置完毕后，封闭空间及其内部物体的边界用Γ表示，边界所包围的流体区域用Ω表示；

步骤二、假定在封闭空间内声源的位置为r，它在单位时间内向单位体积内的空间提供了ρ₀q(r,t)的媒质质量；根据质量守恒定律，媒质中声波的连续方程写为：

$\frac{\∂{\mathrm{\ρ}}^{\′}}{\∂t}+div\left({\mathrm{\ρ}}_{0}v\right)={\mathrm{\ρ}}_{0}q---\left(1\right)$

式中，ρ'为媒质密度增量，ρ₀表示媒质静态密度，q为q(r,t)的简写，v为媒质质点速度，t表示时间，div为散度算子，在三维空间笛卡尔坐标系中，

除了连续性方程之外，用来描述媒质声波的基本方程还有两个，它们不受声源的影响，分别为运动方程：

${\mathrm{\ρ}}_0\frac{\∂v}{\∂t}=-gradp---\left(2\right)$

和物态方程：

$p={c_0}^2{\mathrm{\ρ}}^{\′}---\left(3\right)$

以上两式中，p代表声压，c₀代表声速，grad为梯度算子，在三维空间笛卡尔坐标系中，

由媒质中声波的三个基本方程得到有源情况下封闭空间中有关声压p的波动方程：

${\▿}^{2}p-\frac{1}{c_0^2}\frac{{\∂}^{2}p}{\∂t^2}=-{\mathrm{\ρ}}_0\frac{\∂q}{\∂t}---\left(4\right)$

式中，为拉普拉斯算子，在三维空间笛卡尔坐标系中，

在频域，声源强度q(r,t)表示为

q(r,t)＝q_ω(r)e^jωt (5)

式中，ω为谐振频率，q_ω(r)为在位置r处频域内的声源强度；

封闭空间声场各点的声压的频率与声源相同，声压表示为：

p(r,t)＝p_ω(r)e^jωt (6)

式中，p_ω(r)为在位置r处频域内的声压；

将式(6)及式(7)代入式(5)中，得到简谐声源激励下的声波波动方程为：

${\▿}^2{p}_{\mathrm{\ω}}\left(r\right)e^{j\mathrm{\ω}t}+{\left(\frac{\mathrm{\ω}}{c_0}\right)}^2{p}_{\mathrm{\ω}}\left(r\right)e^{j\mathrm{\ω}t}+{\mathrm{j\ρ}}_0{\mathrm{\ωq}}_{\mathrm{\ω}}\left(r\right)e^{j\mathrm{\ω}t}=0---\left(7\right)$

令式中k＝ω/c₀，称其为波数，并且消去e^jωt，得到只依赖于空间坐标的那部分方程，即室内有源Helmholtz方程：

${\▿}^2{p}_{\mathrm{\ω}}\left(r\right)+k^2{p}_{\mathrm{\ω}}\left(r\right)+j{\mathrm{\ρ}}_0\mathrm{\ω}{q}_{\mathrm{\ω}}\left(r\right)=0---\left(8\right)$

这样就将声压的时域问题转换为频域问题，式(9)即为封闭空间声场的控制方程；

在封闭空间中，边界具有吸声能力，其声压梯度表示为：

$\frac{\∂p}{\∂n}=-\frac{jkp}{\mathrm{\ζ}}---\left(9\right)$

式中，n为封闭空间壁面外法线方向，ζ称为比声阻抗，满足下式：

$\mathrm{\ζ}=\frac{Z}{{\mathrm{\ρ}}_0{c}_0}---\left(10\right)$

式中，Z为界面声阻抗；

根据Galerkin型加权残量法，为了求解式(8)，首先设一试函数为代入有源Helmholtz方程及其边界条件，试函数产生残量R和

$R={\▿}^2\stackrel{\‾}{p}+k^2\stackrel{\‾}{p}+{\mathrm{j\ρ}}_0{\mathrm{\ωq}}_{\mathrm{\ω}}---\left(11\right)$

$\stackrel{\‾}{R}=\frac{\∂\stackrel{\‾}{p}}{\∂n}+\frac{jk}{\mathrm{\ζ}}\stackrel{\‾}{p}---\left(12\right)$

根据伽辽金法确定权函数，有

${\∫}_{\mathrm{\Ω}}\stackrel{\‾}{p}\·({\▿}^2\stackrel{\‾}{p}+k^2\stackrel{\‾}{p}+{\mathrm{j\ρ}}_0{\mathrm{\ωq}}_{\mathrm{\ω}})dv-{\∫}_{\mathrm{\Γ}}\stackrel{\‾}{p}\·(\frac{\∂\stackrel{\‾}{p}}{\∂n}+\frac{jk}{\mathrm{\ζ}}\stackrel{\‾}{p})ds=0---\left(13\right)$

由格林第一公式

式(13)简化为

${\∫}_{\mathrm{\Ω}}(\▿\stackrel{\‾}{p}\·\▿\stackrel{\‾}{p}-k^2\stackrel{\‾}{p}\·\stackrel{\‾}{p}-{\mathrm{j\ρ}}_0\mathrm{\ω}\stackrel{\‾}{p}{q}_{\mathrm{\ω}})dv+{\∫}_{\mathrm{\Γ}}\frac{jk}{\mathrm{\ζ}}\stackrel{\‾}{p}\·\stackrel{\‾}{p}ds=0---\left(15\right)$

在声场中任意一点的声压用各节点声压来表示，即

$\stackrel{\‾}{p}=Np=\[N_1,N_2,...,N_n\]\left\{\begin{array}{c}{p}_1\\ p_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ p_n\end{array}\right\}---\left(16\right)$

式中，N_i为节点i处的形函数，p_i为节点i处的声压；

将式(16)代入式(15)，得

${\∫}_{\mathrm{\Ω}}\[p^T{(\▿N)}^T(\▿N)p-k^2{p}^T{N}^{T}Np-{\mathrm{j\ρ}}_0{\mathrm{\ωp}}^T{N}^T{q}_{\mathrm{\ω}}\]d\mathrm{\Ω}+{\∫}_{\mathrm{\Γ}}\frac{jk}{\mathrm{\ζ}}\left(p^T{N}^{T}Np\right)d\mathrm{\Γ}=0---\left(17\right)$

式中，为形函数的导数矩阵，其表达式为：

$\▿N=\left[\begin{array}{ccc}\frac{\∂N_1}{\∂x}& \frac{\∂N_1}{\∂y}& \frac{\∂N_1}{\∂z}\\ \frac{\∂N_2}{\∂x}& \frac{\∂N_2}{\∂y}& \frac{\∂N_2}{\∂z}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{\∂N_n}{\∂x}& \frac{\∂N_n}{\∂y}& \frac{\∂N_n}{\∂z}\end{array}\right]=\left\{\begin{array}{cccc}{N}_1& N_2& ...& N_n\end{array}\right\}\left\{\begin{array}{ccc}\frac{\∂}{\∂x}& \frac{\∂}{\∂y}& \frac{\∂}{\∂z}\end{array}\right\}---\left(18\right)$

整理式(17)，得到

${\∫}_{\mathrm{\Ω}}\[{\left(\▿N\right)}^T\left(\▿N\right)\]dv\·p-{\∫}_{\mathrm{\Ω}}\left(k^2{N}^{T}N\right)dv\·p-{\∫}_{\mathrm{\Ω}}\left({\mathrm{j\ρ}}_0{\mathrm{\ωN}}^T{q}_{\mathrm{\ω}}\right)dv+{\∫}_{\mathrm{\Γ}}\frac{jk}{\mathrm{\ζ}}\left(N^{T}N\right)ds\·p=0---\left(19\right)$

令

${\∫}_{\mathrm{\Ω}}{(\▿N)}^T(\▿N)\mathrm{d\Ω}=K---\left(20\right)$

$\frac{1}{{c_0}^2}{\∫}_{\mathrm{\Ω}}{N}^{T}Nd\mathrm{\Ω}=M---\left(21\right)$

$\frac{1}{c_0\mathrm{\ζ}}{\∫}_{\mathrm{\Γ}}{N}^{T}Nd\mathrm{\Γ}=C---\left(22\right)$

${\∫}_{\mathrm{\Ω}}{\mathrm{j\ρ}}_0{\mathrm{\ωN}}^T{q}_{\mathrm{\ω}}d\mathrm{\Ω}=G---\left(23\right)$

其中，K称为刚度矩阵，M称为质量矩阵，C称为阻尼矩阵，G称为载荷矩阵；当声源位于位置r₀(x₀,y₀,z₀)处时，频域内的声源强度表示为：

q_ω(r)＝q_ωδ(r-r₀) (24)

其中

$\mathrm{\δ}(r-r_0)=\left\{\begin{array}{cc}0& r\≠r_0\\ 1& r=r_0\end{array}\right.---\left(25\right)$

将式(24)代入式(19)中，得

G＝∫_Ω-jρ₀ωq_ωδ(r-r₀)N^Tdv＝-jρ₀ωq_ωN^T (26)

最后，将式(20)、(21)、(22)、(23)代入式(19)并整理得到

$\frac{(K+j\mathrm{\ω}C-{\mathrm{\ω}}^{2}M)N^{-1}}{\mathrm{\ω}}p=F---\left(27\right)$

式中，F＝jρ₀N^Tq_ω；K、C、M均为n×n阶的系数矩阵，各自的表达式分别为：M＝∫_ΩN^TNdΩ/c₀ ²，C＝∫_ΓN^TNdΓ/c₀ζ，各式中的N为形函数，在实际求解中，K、C、M中的积分运算用求和运算代替 其中m为积分点的数量，m_b为边界上积分点的数量，ξ_i为积分系数；ω为圆频率；c₀为空气中的声速；ζ称为比声阻抗，满足ζ＝Z/ρ₀c₀，ρ₀为空气密度，Z为界面材料的声阻抗；p为封闭空间内任意位置处的声压，实际定位时为单传声器所测得的声信号；F为n×1阶的列向量，表示声源相关信息，其表达式为F＝jρ₀N^Tq_ω，其中q_ω表示声源强度；

步骤三、根据节点坐标，利用移动最小二乘法，获得步骤二中所涉及的形函数N；利用移动最小二乘法构建形函数；一个场函数u(x)在一点的近似值表示为：

$u^h(x,\stackrel{\‾}{x})=\underset{i=1}{\overset{m}{\Σ}}{p}_i\left(\stackrel{\‾}{x}\right)a_i\left(x\right)=p^T\left(\stackrel{\‾}{x}\right)a\left(x\right)---\left(28\right)$

其中是计算点x的邻域范围内各节点的坐标，为基函数向量，m为基函数的个数，a(x)＝[a₁(x),a₂(x),…a_m(x)]为待定系数向量；使用单项式基函数做运算，在三维空间中常用的线性及二次单项式基函数分别为：

$\begin{array}{c}p\left(x\right)={\[1,x,y,z\]}^T,m=4\\ p\left(x\right)={\[1,x,y,z,x^2,xy,y^2,yz,z^2,xz\]}^T,m=10\end{array}---\left(29\right)$

将求解域用节点离散后，在每个节点处定义一个权函数该函数只在支撑域内不为零，在支撑域之外为零，在三维情况下，权函数的支撑域为球形；选定权函数后，就求得近似函数在节点处的误差加权平方和：

$J=\underset{I=1}{\overset{N}{\Σ}}{w}_I\left(x\right){\[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{p}_i\left(\stackrel{\‾}{x_I}\right)a_i\left(x\right)-u_I\]}^2---\left(30\right)$

令J取最小值，即

$\begin{array}{c}\frac{\∂J}{\∂a_j\left(x\right)}=2\underset{I=1}{\overset{N}{\mathrm{\Σ}}}{w}_I\left(x\right)\[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{p}_i\left(\stackrel{\‾}{x_I}\right)a_i\left(x\right)-u_I\]p_j\left(\stackrel{\‾}{x_I}\right)\\ =0,j=1,2,\mathrm{...},m\end{array}---\left(31\right)$

经过整理后，得到下式：

A(x)a(x)＝B(x)u (32)

式中，A(x)，B(x)的含义为：

$\begin{array}{c}A\left(x\right)=\underset{I=1}{\overset{N}{\mathrm{\Σ}}}{w}_I\left(x\right)p\left(\stackrel{\‾}{x_I}\right)p^T\left(\stackrel{\‾}{x_I}\right)\\ B\left(x\right)=\[w_1\left(x\right)p\left(\stackrel{\‾}{x_1}\right),w_2\left(x\right)p\left(\stackrel{\‾}{x_2}\right),...,w_N\left(x\right)p\left(\stackrel{\‾}{x_N}\right)\]\end{array}---\left(33\right)$

由公式(32)得到a(x)，将其代入式(28)得：

$u^h(x,\stackrel{\‾}{x})=p^T\left(\stackrel{\‾}{x}\right)A^{-1}\left(x\right)B\left(x\right)u=N(x,\stackrel{\‾}{x})u---\left(34\right)$

步骤四、在封闭空间内部任意位置处设置一传声器；当声源有声音发出时，传声器拾取到一段音频信号f(t)，将其进行傅里叶变换，得到此音频信号的频域信号F(ω)；

步骤五、将步骤四中得到的频域信号F(ω)作为式(29)中的p值代入式(29)，求解后得到n×1阶的列向量F，由于F中的j、ρ₀、q_ω均为常数，因此，求得的列向量中值最大的元素所代表的节点的位置即为声源的位置。