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

摘要

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

主权项

一种基于封闭空间几何信息建模的单传声器声源定位方法，其特征在于包括以下步骤：步骤一、在封闭空间边界所包围的流体区域内布置n个节点，将此n个节点从1到n进行编号；n个节点一方面作为封闭空间的几何信息在后期进行建模计算，一方面用来表示未来声源定位时的位置；n个节点近似均匀分布于整个流体区域内，封闭空间的边界处需布有节点；当封闭空间内部存在物体时，物体内部不布置节点，且物体边界上应布置有节点，即节点应勾勒出物体形状；节点布置完毕后，封闭空间及其内部物体的边界用Γ表示，边界所包围的流体区域用Ω表示；步骤二、假定在封闭空间内声源的位置为r，它在单位时间内向单位体积内的空间提供了ρ0q(r,t)的媒质质量；根据质量守恒定律，媒质中声波的连续方程写为：∂ρ′∂t+div(ρ0v)=ρ0q---(1)]]>式中，ρ'为媒质密度增量，ρ0表示媒质静态密度，q为q(r,t)的简写，v为媒质质点速度，t表示时间，div为散度算子，在三维空间笛卡尔坐标系中，除了连续性方程之外，用来描述媒质声波的基本方程还有两个，它们不受声源的影响，分别为运动方程：ρ0∂v∂t=-gradp---(2)]]>和物态方程：p=c02ρ′---(3)]]>以上两式中，p代表声压，c0代表声速，grad为梯度算子，在三维空间笛卡尔坐标系中，由媒质中声波的三个基本方程得到有源情况下封闭空间中有关声压p的波动方程：▿2p-1c02∂2p∂t2=-ρ0∂q∂t---(4)]]>式中，为拉普拉斯算子，在三维空间笛卡尔坐标系中，在频域，声源强度q(r,t)表示为q(r,t)＝qω(r)ejωt                       (5)式中，ω为谐振频率，qω(r)为在位置r处频域内的声源强度；封闭空间声场各点的声压的频率与声源相同，声压表示为：p(r,t)＝pω(r)ejωt                       (6)式中，pω(r)为在位置r处频域内的声压；将式(6)及式(7)代入式(5)中，得到简谐声源激励下的声波波动方程为：▿2pω(r)ejωt+(ωc0)2pω(r)ejωt+jρ0ωqω(r)ejωt=0---(7)]]>令式中k＝ω/c0，称其为波数，并且消去ejωt，得到只依赖于空间坐标的那部分方程，即室内有源Helmholtz方程：▿2pω(r)+k2pω(r)+jρ0ωqω(r)=0---(8)]]>这样就将声压的时域问题转换为频域问题，式(9)即为封闭空间声场的控制方程；在封闭空间中，边界具有吸声能力，其声压梯度表示为：∂p∂n=-jkpζ---(9)]]>式中，n为封闭空间壁面外法线方向，ζ称为比声阻抗，满足下式：ζ=Zρ0c0---(10)]]>式中，Z为界面声阻抗；根据Galerkin型加权残量法，为了求解式(8)，首先设一试函数为代入有源Helmholtz方程及其边界条件，试函数产生残量R和R=▿2p‾+k2p‾+jρ0ωqω---(11)]]>R‾=∂p‾∂n+jkζp‾---(12)]]> 根据伽辽金法确定权函数，有∫Ωp‾·(▿2p‾+k2p‾+jρ0ωqω)dv-∫Γp‾·(∂p‾∂n+jkζp‾)ds=0---(13)]]> 由格林第一公式 式(13)简化为∫Ω(▿p‾·▿p‾-k2p‾·p‾-jρ0ωp‾qω)dv+∫Γjkζp‾·p‾ds=0---(15)]]> 在声场中任意一点的声压用各节点声压来表示，即p‾=Np=[N1,N2,...,Nn]p1p2...pn---(16)]]> 式中，Ni为节点i处的形函数，pi为节点i处的声压； 将式(16)代入式(15)，得∫Ω[pT(▿N)T(▿N)p-k2pTNTNp-jρ0ωpTNTqω]dΩ+∫Γjkζ(pTNTNp)dΓ=0---(17)]]>式中，为形函数的导数矩阵，其表达式为：▿N=∂N1∂x∂N1∂y∂N1∂z∂N2∂x∂N2∂y∂N2∂z.........∂Nn∂x∂Nn∂y∂Nn∂z=N1N2...Nn∂∂x∂∂y∂∂z---(18)]]> 整理式(17)，得到∫Ω[(▿N)T(▿N)]dv·p-∫Ω(k2NTN)dv·p-∫Ω(jρ0ωNTqω)dv+∫Γjkζ(NTN)ds·p=0---(19)]]> 令∫Ω(▿N)T(▿N)dΩ=K---(20)]]>1c02∫ΩNTNdΩ=M---(21)]]>1c0ζ∫ΓNTNdΓ=C---(22)]]>∫Ωjρ0ωNTqωdΩ=G---(23)]]> 其中，K称为刚度矩阵，M称为质量矩阵，C称为阻尼矩阵，G称为载荷矩阵；当声源位于位置r0(x0,y0,z0)处时，频域内的声源强度表示为： qω(r)＝qωδ(r‑r0)                      (24) 其中δ(r-r0)=0r≠r01r=r0---(25)]]> 将式(24)代入式(19)中，得 G＝∫Ω‑jρ0ωqωδ(r‑r0)NTdv＝‑jρ0ωqωNT            (26) 最后，将式(20)、(21)、(22)、(23)代入式(19)并整理得到(K+jωC-ω2M)N-1ωp=F---(27)]]>式中，F＝jρ0NTqω；K、C、M均为n×n阶的系数矩阵，各自的表达式分别为：M＝∫ΩNTNdΩ/c02，C＝∫ΓNTNdΓ/c0ζ，各式中的N为形函数，在实际求解中，K、C、M中的积分运算用求和运算代替其中m为积分点的数量，mb为边界上积分点的数量，ξi为积分系数；ω为圆频率；c0为空气中的声速；ζ称为比声阻抗，满足ζ＝Z/ρ0c0，ρ0为空气密度，Z为界面材料的声阻抗；p为封闭空间内任意位置处的声压，实际定位时为单传声器所测得的声信号；F为n×1阶的列向量，表示声源相关信息，其表达式为F＝jρ0NTqω，其中qω表示声源强度； 步骤三、根据节点坐标，利用移动最小二乘法，获得步骤二中所涉及的形函数N；利用移动最小二乘法构建形函数；一个场函数u(x)在一点的近似值表示为：uh(x,x‾)=Σi=1mpi(x‾)ai(x)=pT(x‾)a(x)---(28)]]>其中是计算点x的邻域范围内各节点的坐标，为基函数向量，m为基函数的个数，a(x)＝[a1(x),a2(x),…am(x)]为待定系数向量；使用单项式基函数做运算，在三维空间中常用的线性及二次单项式基函数分别为：p(x)=[1,x,y,z]T,m=4p(x)=[1,x,y,z,x2,xy,y2,yz,z2,xz]T,m=10---(29)]]>将求解域用节点离散后，在每个节点处定义一个权函数该函数只在支撑域内不为零，在支撑域之外为零，在三维情况下，权函数的支撑域为球形；选定权函数后，就求得近似函数在节点处的误差加权平方和：J=ΣI=1NwI(x)[Σi=1mpi(xI‾)ai(x)-uI]2---(30)]]> 令J取最小值，即∂J∂aj(x)=2ΣI=1NwI(x)[Σi=1mpi(xI‾)ai(x)-uI]pj(xI‾)=0,j=1,2,...,m---(31)]]> 经过整理后，得到下式： A(x)a(x)＝B(x)u                       (32) 式中，A(x)，B(x)的含义为：A(x)=ΣI=1NwI(x)p(xI‾)pT(xI‾)B(x)=[w1(x)p(x1‾),w2(x)p(x2‾),...,wN(x)p(xN‾)]---(33)]]> 由公式(32)得到a(x)，将其代入式(28)得：uh(x,x‾)=pT(x‾)A-1(x)B(x)u=N(x,x‾)u---(34)]]> 步骤四、在封闭空间内部任意位置处设置一传声器；当声源有声音发出时，传声器拾取到一段音频信号f(t)，将其进行傅里叶变换，得到此音频信号的频域信号F(ω)； 步骤五、将步骤四中得到的频域信号F(ω)作为式(29)中的p值代入式(29)，求解后得到n×1阶的列向量F，由于F中的j、ρ0、qω均为常数，因此，求得的列向量中值最大的元素所代表的节点的位置即为声源的位置。