[发明专利]基于封闭空间几何信息建模的单传声器声源定位方法

权利要求书

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_ω均为常数，因此，求得的列向量中值最大的元素所代表的节点的位置即为声源的位置。