本发明涉及一种室内噪声源定位方法，尤其涉及一种利用球面近场声全息算法重构声源附近表面的无功声强分布来进行室内噪声源定位方法。
背景技术：
小型房间内部噪声过大会危害工作人员的身体健康，因此对室内噪声源进行有效治理是至关重要的工作，对噪声源进行有效控制的必要前提是对噪声源进行准确定位。
在噪声源定位领域，采用传声器阵列进行声场测量是最常用的手段。比较常用的是平面传声器阵列，但平面阵的局限在于不能区分阵列前方或后方的声音。小型房间内部由于壁面和天花板的反射作用，使得内部声场形成混响场，在如此复杂的声场环境中，平面阵很难对噪声源进行有效定位。球形阵列具有经纬角全指向性和三维对称性，能一次性测量三维空间各个方向的声场信息。球形阵特殊的结构使其具有很大的灵活性，测量速度快，因此可用于潜艇、飞机、汽车等内部空间的声场测量。
近场声全息技术是基于传声器阵列的噪声源定位识别算法，对于低频近场的噪声源定位问题具有很好的适用性。在自由场条件下，采用传声器阵列测量的三维空间声场数据可作为近场声全息技术的全息面数据，通过该算法重构出声源附近表面的声压分布来确定噪声源的位置。但在小型房间等内部空间中，由于天花板和壁面的反射等作用，使得内部声场会产生多阶强模态叠加成为混响场，此时适用于自由场的近场声全息算法重构声源附近表面的声压往往不能确定噪声源位置。
技术实现要素：
本发明要解决的技术问题是：提供一种利用目标声源附近表面的复声强的无功声强部分来对室内噪声源进行准确定位的方法，以克服现有技术的不足。
本发明的技术方案是：一种基于球面近场声全息重构无功声强的室内噪声源定位方法，包括以下步骤：
(1)建模仿真计算规则小型房间内部由机械声源辐射的三维空间声场数据；
(2)采用阵元数为q的球面传声器阵列对目标声源辐射到三维空间的声压信号进行采集p(r,t)；
(3)对球形阵的q个传声器数据分别和第1个传声器数据进行互谱分析，获取整个球形阵的频域复声压p(r,θ,φ)；
(4)基于球形阵的频域复声压，利用球面近场声全息算法重构声源面的声压p(rs,θ,φ)和质点振速vn(rs,θ,φ)；
(5)利用复声压和质点振速重构声源附近表面复声强ic(r,w)，提取其无功声强部分ir(rs,w)；
(6)通过声源附近表面的无功声强分布ir(rs,w)确定室内噪声源位置。
所述步骤(2)中由q个传声器组成球形阵列整个球形阵测量数据的复声压数据为：
其中：以随机球形阵中心为坐标原点(0,0,0)；
球形阵阵元个数q，阵元编号为1,…,q,…q，第q个传声器的空间坐标为(xq,yq,zq)；
p0为距离点声源1m处的声压幅值，单位为pa；
ω为角频率，rad/s；t为时间点，s；k为波数。假设空间有k个点源目标入射到球形阵列上，则每个传声器接收的数据为k个点源目标的叠加。
所述步骤(3)中对球形阵的q个传声器数据分别和第1个传声器数据进行互谱分析，获取整个球形阵的频域复声压为：
所述步骤(4)中采用球面近场声全息算法重构声源附近表面的声压分布的计算公式如下：
式中：rs为重构面半径，rh为全息面半径；
谐函数，重构球面声压球面波谱pnm(rs)的计算公式如下：
式中：pnm(rs)为全息面的球面波谱；
jn(kr)为第一类球bessel函数，jn(krs)/jn(krh)为全息面球面波谱pnm(rs)与重构面球面波谱pnm(rs)之间的传递函数，全息面的球面波谱的计算公式如下：
式中：(rq,θq,φq)为第q个传声器位置的坐标，pq(rq,θq,φq)为第q个传声器采样获取的声压；αq＝4πa2/q为权重系数，是每个传声器位置对应的球面网格面积；n为需要计算的球面波谱的最高阶数；
所述步骤(4)中采用球面近场声全息算法重构声源附近表面的质点振速分布为：
式中：i2＝-1；ρ0c0为空气的介质特性阻抗；
j'n(kr)为n阶第一类球贝塞尔函数的导数。
所述步骤(5)声源附近表面的复声强用复声压分布和复质点振速分布获得，其计算计算公式如下：
ic(r,w)＝p(rs,w)·vn(rs,w)
式中：p(rs,w)为声压的fourier变换；u(rs,w)为振速的fourier变换；
所述复声强分解为有功声强部分和无功声强部分，用下式表示：
ic(r,w)＝ia(rs,w)+i·ir(rs,w)
式中：ia(rs,w)为有功声强部分，可用来表示声波向远处传播的声能；
ir(rs,w)为无功声强部分，表示不传播的声能。
本发明的工作原理是：有功声强与无功声强在近场的强弱有区别，无功声强在近场较强，在远场较弱，无功声强可用来反映近场声场的复杂性。
与现有技术相比，本发明具有以下优点：
1)球面阵经纬角全指向性和三维对称性，能一次性测量三维空间各个方向的声场信息。
2)利用无功声强在声源近场区域较强的特点，能够排除室内混响环境的影响，对室内噪声源进行有效定位。
附图说明
图1为本发明的球阵数据接收模型；
图2为本发明的室内噪声源定位模型；
图3为本发明的流程图；
图4为本发明的100hz声场计算结果；
图5为本发明的400hz声场计算结果；
图6为本发明的100hz噪声源定位结果；(a)f＝100hz声源面理论声压分布；(b)重构球面无功声强分布；(c)重构球面无功声强分布展开；
图7为本发明的400hz噪声源定位结果；(a)声源面理论声压分布；(b)重构球面无功声强分布；(c)重构球面无功声强分布展开。
具体实施方式
下面结合附图及具体的实施例对发明进行进一步介绍：
本发明的流程如图3所示，首先计算声源在给定尺寸房间的辐射声场，提取球面阵个采集点的数据，将q个通道数据分别与第1通道数据做互谱分析，获取球面阵的频域复声压，采用球面近场声全息算法重构声源附近表面的声压分布与质点振速分布，通过复声强计算公式得到声源附近表面的复声强，提取其虚部为其无功部分，利用声源附近表面的无功声强分布确定噪声源位置。具体方法如下：
(1)建模仿真计算规则小型房间内部由机械声源辐射的三维空间声场数据。
利用多物理场耦合软件comsol建模仿真，球阵数据接收模型如图1所示，在三维空间中以球面阵球心为坐标原点(0,0,0)，点声源位于(x0,y0,z0)。室内噪声源定位模型如图2所示，在一个长lx＝4m,ly＝3m,lz＝3m的房间内，四面墙均为硬边界条件，球面阵放置于房间内，以球面阵球心为坐标原点(0,0,0),在距离球面阵球心放为0.6m处放置一点声源，正对球阵上第25号传声器位置，其位置信息在直角坐标系下为(x0,y0,z0)＝(-0.5619,0.2100,0.0075)，对应球坐标下为(r0,θ0,φ0)＝(0.6m,89°,160°)。
(2)采用阵元数为q的球面传声器阵列对目标声源辐射到三维空间的声压信号进行采集p(r,t)。
由q个传声器组成球形阵列整个球形阵测量数据的复声压数据为：
其中：以随机球形阵中心为坐标原点(0,0,0)；
球形阵阵元个数q，阵元编号为1,…,q,…q，第q个传声器的空间坐标为(xq,yq,zq)；
p0为距离点声源1m处的声压幅值，pa；
ω为角频率，rad/s；t为时间点，s；k为波数。
假设空间有k个点源目标入射到球形阵列上，则每个传声器接收的数据为k个点源目标的叠加。
(3)将球形阵的q个传声器数据分别和第1个传声器数据进行互谱分析，获取整个球形阵的频域复声压p(r,θ,φ)。公式如下：
(4)基于球形阵的频域复声压，利用球面近场声全息算法重构声源面的声压p(rs,θ,φ)和质点振速vn(rs,θ,φ)。
其中采用球面近场声全息算法重构声源附近表面的声压分布的计算公式如下：
式中：rs为重构面半径，rh为全息面半径；
谐函数，重构球面声压球面波谱pnm(rs)的计算公式如下：
式中：pnm(rs)为全息面的球面波谱；
jn(kr)为第一类球bessel函数，jn(krs)/jn(krh)为全息面球面波谱pnm(rs)与重构面球面波谱pnm(rs)之间的传递函数，全息面的球面波谱的计算公式如下：
式中：(rq,θq,φq)为第q个传声器位置的坐标，pq(rq,θq,φq)为第q个传声器采样获取的声压；αq＝4πa2/q为权重系数，是每个传声器位置对应的球面网格面积；n为需要计算的球面波谱的最高阶数；
采用球面近场声全息算法重构声源附近表面的质点振速分布为：
式中：i2＝-1；ρ0c0为空气的介质特性阻抗；
j'n(kr)为n阶第一类球贝塞尔函数的导数。
(5)利用复声压和质点振速重构声源附近表面复声强ic(r,w)，提取其无功声强部分ir(rs,w)。
声源附近表面的复声强用复声压分布和复质点振速分布获得，其计算计算公式如下：
ic(r,w)＝p(rs,w)·vn(rs,w)
式中：p(rs,w)为声压的fourier变换；u(rs,w)为振速的fourier变换；
所述复声强分解为有功声强部分和无功声强部分，用下式表示：
ic(r,w)＝ia(rs,w)+i·ir(rs,w)
式中：ia(rs,w)为有功声强部分，可用来表示声波向远处传播的声能；
ir(rs,w)为无功声强部分，表示不传播的声能。
(6)通过声源附近表面的无功声强分布ir(rs,w)确定室内噪声源位置。
房间模型三维声场计算结果如图4和图5所示，首先设置计算参数和边界条件，将球阵球心作为坐标原点，在距离球心的0.6m，正对球形阵25号传声器位置放置一个激励点声源，声源辐射频率f分别为100hz和400hz，计算出房间内各点声场信息，提取出球阵表面采集点的声压数据和声源半径所在球面的声压数据。
声源面理论声压分布如图6所示，通过球面近场声全息算法重构出的声源附近表面的无功声强分布如图7所示。将球面阵表面的声压数据作为全息数据，重构出声源附近表面得到声压和质点振速分布，进而求出声源附近表面的无功声强分布，实现噪声源的准确定位。
图6和图7分别是频率为100hz和400hz的噪声源定位结果，图中“+”表示理论声源位置，直角坐标系下(x0,y0,z0)＝(-0.5619,0.2100,0.0075)，球坐标系下(r0,θ0,φ0)＝(0.6m,89°,160°)。图中(a)对应的是从声场计算结果中提取出声源半径所在球面的声压分布，图中(b)(c)为利用球面阵测量数据重构出的无功声强分布对噪声源进行定位的结果，其中(b)为直角坐标系下的定位结果，(c)为球坐标系下的定位结果，分别对比图6和图7的(a)和(c)可知采用声源附近表面的无功声强分布可对噪声源进行准确定位。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。