CN106125047A - 一种气流环境下的声源定位方法 - Google Patents

Abstract

本发明提供一种气流环境下的声源定位方法，属于声源定位技术领域。该方法步骤如下：搭建传声器阵列，采集气流环境下的声信号；根据阵列传声器坐标等几何参数，借助于Amiet模型，求解各传声器到扫描点的声传播路径，得到气流环境下的阵列流形矩阵A；预估声源个数并构造降噪后的互谱矩阵，在扫描点声压矩阵、阵列流形矩阵及降噪后的互谱矩阵间建立价值函数；通过凸优化迭代求解声压矩阵，最终绘制出声压重构云图，获得声源位置信息。本发明能够准确定位气流环境下的声源位置，修正气流造成的定位偏差，定位分辨率高，能够实现多相干声源定位，降低了凸优化迭代过程中的计算量，减少了重构云图中旁瓣数，声源定位效果良好。

Description

一种气流环境下的声源定位方法

技术领域

本发明涉及声源定位技术领域，特别是指一种气流环境下的声源定位方法。

背景技术

近年来，民用、军用飞行器等航空运输快速发展，航空气动噪声所引发的一系列问题开始逐渐引起人们的关注。机身噪声成为严重制约飞行器设计的因素。机身噪声主要来自于升降装置，起落架，轮舱以及机翼后缘等，严重威胁人体健康，影响附近居民日常生活。同时，高速的飞行器在高声强的热噪声激励下，飞行器构件材料内部会产生快速交变应力，局部结构会出现热噪声疲劳破坏，严重影响飞行器服役性能。为提高飞行器服役安全性，降低航空噪声对周围环境造成的不良影响，需开展相应的气动噪声测量和评估工作以优化设计、降低噪声。因此，确定噪声源位置的高分辨率定位方法显得极为重要。

由于飞行器噪声实测成本较高，测试过程中各种环境参数不易控制，因此，在风洞内进行飞行器缩比模型气动噪声测试获得了广泛应用。在风洞内进行噪声测试时，缩比模型表面产生的气动热噪声，经过气流及其边界层后被气流外的传声器阵列接收到，声波的传播路径受气流流速及剪切层折射作用而发生改变。常规的波束形成算法在定位声源的过程中，未考虑气流环境对阵列流形向量的影响，造成定位结果偏差较大，且存在分辨率低，旁瓣值多，计算量大等问题。

发明内容

本发明要解决的技术问题是提供一种气流环境下的声源定位方法，实现了气动噪声源的准确、高效定位。该方法主要用于汽车、飞机等风洞试验中，用于确定噪声源位置。

为实现上述目的，本发明公开以下技术方案：

该方法通过建立三维剪切流模型，计算气流影响下的阵列流形矩阵A，进而构造扫描点声压矩阵、阵列流形矩阵及降噪后互谱矩阵间的价值函数，设置自适应更新的稀疏化约束条件，凸优化求解扫描点声压矩阵P并绘制声压重构云图，进而获得声源位置信息。

所述方法主要包括以下步骤：

(一)搭建传声器阵列，采集声信号；

(二)建立三维模型，计算剪切流作用下的阵列流形矩阵A；

(三)预估声源个数，构造降噪后的互谱矩阵并建立价值函数；

(四)设置稀疏化约束条件，凸优化求解扫描点声压矩阵并绘制重构云图。

其中，传声器阵列为7条臂组成的阿基米德螺旋阵，其上分布28个驻极体自由场传声器，阵列布置于气流外安全区域，利用基于NI-PXIe总线的数据采集系统完成传声器阵列的多通道同步采集，采样频率44100Hz。

阵列流形矩阵A考虑了气流影响造成的声线偏折问题，具体计算过程如下：

由Snell折射定律和速度矢量合成定理，建立气流环境下的声传播路径模型：

其中，M₀表示气流的马赫数，c为常温下声传播速度，i₁为气流外声波传入气流内部时初始入射方向，i₂声波进入气流后折射方向，i₃为气流内矢量合成后的速度方向，i₅为声波传出气流时折射方向，通过离散化和迭代方法，对式(1)进行求解，进而求得扫描点与传声 器之间的近似声程r和声传播时延τ为：

其中，L₁为扫描声源面到气流一侧边缘距离，W_i为气流区域宽度，L₂为传声器阵列到气流另一侧边缘距离，c₁为矢量合成后的速度；

为求解三维问题，将各传声器单独取出，在其与扫描点单独组成的平面内进行传播时延的计算，以此来完成三维声传播时延的求解，获得阵列流行矩阵A：

步骤(三)具体为：

对采集到的声信号数据y(b)求解其互谱矩阵G＝y(b)y^H(b)，然后对阵列信号对阵列信号互谱矩阵G进行特征分解G＝UΛU^H，根据分解所得特征值预估声源个数取其前个较大的特征值和对应的特征向量，构造降噪后的互谱矩阵利用二阶锥规划模型建立扫描点声压矩阵P、阵列流形矩阵A及降噪后互谱矩阵间的价值函数(4)，通过最小化价值函数迭代求解声压矩阵P：

式中，Q为一个任意酉矩阵，XX^H＝P，x为矩阵X中所有元素组成的向量，各变量的脚标n表示迭代次数，ξ_n为稀疏化约束参数。

其中，稀疏化约束条件||x_n||₁≤ξ_n约束了扫描声源面信号的总功率大小，稀疏化约束参数ξ_n的值越小，就会有尽可能多的非声源位置的自相关函数值为零，稀疏化约束参数设置遵循以下原则：

当n＝1时，约束参数ξ_n的值可通过柯西-施瓦茨不等式(5)获得：

式中，x_l为矩阵X中的元素，β＝tr(Λ)；

当n＞1时，约束参数ξ_n可由式(6)实现自适应更新，达到更好的稀 疏化效果:

ξ_n＝||x_n-1||₁(6)

该方法步骤(四)中通过凸优化循环方式最小化价值函数获得扫描点声压矩阵P，直至连续两次迭代求解结果小于设定的阈值时，迭代计算停止，利用所获得的矩阵P绘制扫描点声压重构云图，对比各重构点声压值大小确定声源位置。

本发明的上述技术方案的有益效果如下：

利用Amiet方法对三维空间下的声传播路径进行了修正，提高了阵列流形向量A的准确度；与原始定位算法相比，改进后的算法在每一步迭代过程中能自适应地更改稀疏化约束参数，加强了算法的稀疏化程度，提高了计算效率，减小了旁瓣数量；在气动噪声源定位中，该方法具有较高的分辨率，且对于多相干声源具有良好的辨识能力，有利于气动噪声源的精准定位。

传统方法未考虑声源定位时气流造成的偏差，虽已有研究过气流修正问题，但修正模型简化为圆柱型气流，成轴对称形状，而汽车、飞机等风洞测试气流多为方形截面，较为复杂，圆柱型模型不再适用，本方法是基于方形截面气流，同时更适用于简化的圆柱形气流。

附图说明

图1为本发明的气流环境下的声源定位方法的定位流程示意图；

图2为本发明的气流环境下的声源定位方法的定位设备分布示意图；

图3为本发明的剪切流作用下声波传播路径示意图；

图4为本发明气流环境下三维空间声线折射示意图；

图5为本发明双声源定位声压重构云图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明提供一种气流环境下的声源定位方法，流程如下：

首先根据所需测量的对象及其所处流速环境，确定能满足声压级和频率响应范围的传声器；然后搭建阿基米德螺旋形传声器阵列，通过BNC线缆连接NI-PXIe的多通道数据采集设备，通过labview软件实现多通道声信号的数据采集并存储。

其次对数据信号进行后续处理，将环境参数和采集到的声信号数据导入matlab软件进行编程计算，求解阵列流行向量以及声压矩阵P，根据求解的声压矩阵P绘制出重构云图，以此确定噪声源位置。同时也可在labview采集软件中添加MathScript程序模块实现对采集数据的实时计算，达到声源位置的动态定位。

如图1所示，为该方法的流程图，具体步骤如下：

(一)搭建传声器阵列，采集声信号；

在本实施例中，对消音风洞内的两个直径为9.5cm的球形扬声器进行定位，气流截面为正方形，截面边长20cm，将两扬声器置于距离气流中心线一侧的15cm处，两声源间距25cm，传声器阵列置于距离气流中心线另一侧25cm处，各装置位置分布如图2所示。

所搭建传声器阵列为7条臂组成的阿基米德螺旋阵，其上分布28个驻极体自由场传声器，阵列置于气流外安全区域，利用基于NI-PXIe总线的数据采集系统完成传声器阵列的多通道同步采集，采样频率44100Hz，在流速为80m/s时，采集两扬声器所发出的相干纯音信号，并存储数据。

(二)建立三维模型，计算剪切流作用下阵列流形矩阵A；

利用Amiet方法对三维空间下的声传播路径进行修正，以提高阵列流形向量的准确度；

在气流作用下，声波传播路径示意图如图3所示，设扫描点与一侧剪切层距离为L₁，传声器阵列与另一侧气流剪切层的距离为L₂，扫描点F与第m个传声器的水平距离为St；

从扫描点F处发出一条与水平方向夹角为i₁的声线，即FB段，该声线在穿越气流剪切层时发生了折射，其折射角i₂为线段BC与剪 切层的夹角；由于气流的存在，声线不会沿着BC所指方向传播，而是会往气流下游方向偏移，即沿着线段BD传播至气流另一侧剪切层上的D点处，设线段BD与水平方向夹角为i₃；最后，声波从D点穿越剪切层进而发生第二次折射，沿直线DE方向传播至麦克风处，其折射角i₅为线段DE与水平方向夹角。显然，为了准确计算出气流环境下的声程，需要求出图3中所示的FB，BD和DE这三条线段长度之和r；

根据Snell折射定律和速度矢量合成定理，建立气流环境下声传播路径模型：

其中，M₀表示气流的马赫数，c为常温下声传播速度，通过离散化和迭代方法，对式(1)进行求解，进而求得扫描点与传声器之间的近似声程r和声传播时延τ为：

τ＝(|FB|+|BC|+|DE|)/c＝(L₁/sin i₁+W_i/sin i₂+L₂/sini₅)/c(8)

需要注意的是，在式(2)中计算声波在气流中的传播路径时所用的线段是BD，而在计算声波在气流中的传播时间所用的线段是BC，这是因为声波在气流中沿BD方向的运动是声波以声速c沿BC方向和以流速V₀沿水平方向两种运动的合运动，因此声波在气流内的运动时间就等于声波以声速c沿BC方向的运动时长。在图3所示的模型中，传声器、扫描点都处在同一水平面。当要计算的扫描点同传声器不在同一水平面时，如图4所示，各传声器到扫描点的声传播路径计算便成为一个三维的声线路径求解问题。为了求解三维问题，可将各传声器单独取出，在其与扫描点单独组成的平面内进行传播时延的计算，以此来完成三维声传播时延的求解，具体操作步骤如下所述：

如图4所示，设在声源面内有一个扫描点B，发出一条与x轴，y轴夹角分别为(σ,β)的声线，该声线在穿越气流所在区域后被传声器阵列某一阵元Q接收到。设气流区域内马赫数均为M₀，气流区域的宽度为W，且传声器Q和扫描点B并不在同一水平面上。此时，在阵列面上，过传声器Q做一条水平直线l₁。在声源所在平面上，过扫描点B，做另外一条水平线直线l₂。由于阵列面和声源面平行，所以水平直线l₁和l₂一定平行，由l₁和l₂构成的平面与矩形气流的交平面即为图4所示的矩形CDEF，其宽度为W′。之后，可根据图3的模型和式(1)、式(2)的计算步骤进行计算，获得在气流环境下声波传播路径r和声传播时延τ，进而可以计算出经过修正的阵列流形矩阵A：

(三)预估声源个数，构造降噪后的互谱矩阵并建立价值函数；

对采集到的声信号数据y(b)求解其互谱矩阵G＝y(b)y^H(b)，然后对互谱矩阵G进行特征分解G＝UΛU^H，将分解后得到的M个特征值从大到小排列，即λ₁＞λ₂＞...＞λ_M，令 组成的序列反映了互谱矩阵特征值之间的突变程度，预估声源个数为 取值应使得其次，取其前个较大的特征值和对应的特征向量，构造降噪后的互谱矩阵利用二阶锥规划模型建立扫描点声压矩阵P、阵列流形矩阵A及降噪后互谱矩阵间的价值函数(4)，通过最小化价值函数迭代求解声压矩阵P：

式中，Q为一个任意酉矩阵，XX^H＝P，x为矩阵X中所有元素组成的向量，各变量的脚标n表示迭代次数，ξ_n为稀疏化约束参数。

由于无需求解信号互谱矩阵G的逆矩阵，而是对G进行特征分解，再通过凸优化方式求解声压矩阵P，因此所述方法可以用于矩阵G在秩亏损情况下的相干声源的定位；同时，降噪后的互谱矩阵具有较高的稀疏性，使得计算效率提高。

(四)设置稀疏化约束条件，凸优化求解扫描点声压矩阵并绘制重构云图；

稀疏化约束条件||x_n||₁≤ξ_n约束了扫描声源面信号的总功率大小，稀疏化约束参数ξ_n的值越小，就会有尽可能多的非声源位置的自相关函数值为零，对MACS算法进行改进，在每一步迭代过程中自适应地修正稀疏化约束参数ξ_n，增强算法的稀疏化程度，提高计算效率，减少旁瓣个数，进一步提高气动噪声源的定位精度，稀疏化约束参数设置遵循以下原则：

当n＝1时，约束参数ξ_n的值可通过柯西-施瓦茨不等式(5)获得：

式中，x_l为矩阵X中的元素，β＝tr(Λ)；

当n＞1时，约束参数ξ_n可由式(6)实现自适应更新，达到更好的稀疏化效果:

ξ_n＝||x_n-1||₁(12)

所述方法通过凸优化循环方式，最小化价值函数获得扫描点声压矩阵P，具体计算过程如下：

首先，令n＝1，根据柯西-施瓦茨不等式，初始化稀疏算子并预先设置维的单位矩阵Q₁，求解式(4)得到X_n；

其次，将求解得到的X_n代入式(7)得到更新后的Q_n+1：

其中，将矩阵X_n中所有元素组成的向量x_n代入式(6)得到更新后的稀疏化约束参数ξ_n+1；

将更新后的Q_n+1和ξ_n+1代入式(4)，开始下一次的凸优化过程，求解得到新的X_n+1，按照上述过程反复迭代，根据可以求得声压矩阵P_n，直到小于设定的阈值(如10^-10)为止，此时迭代终止，利用所获得的矩阵P绘制扫描点声压重构云图，对比各重构点声压值大小确定声源位置。

如图5所示，为所述方法对气流影响下的双声源定位结果，其中 两圆环位置为扬声器实际所在的物理位置，通过对比修正前后的定位结果可以看出，阵列流行矩阵未作修正时，由于剪切流的影响会出现定位偏差，重构声源位置在圆环外区域；对阵列流行矩阵修正后，能够较为准确定位声源位置，且具有较高的分辨率，表明所述方法在气动噪声源定位中，具有较好的辨识能力和适用性。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种气流环境下的声源定位方法，其特征在于：具体步骤如下：

(一)搭建传声器阵列，采集声信号；

(二)建立三维模型，计算剪切流作用下的阵列流形矩阵A；

(三)预估声源个数，构造降噪后的互谱矩阵并建立价值函数；

(四)设置稀疏化约束条件，凸优化求解扫描点声压矩阵并绘制重构云图。

2.根据权利要求1所述的气流环境下的声源定位方法，其特征在于：所述步骤(一)中的传声器阵列为7条臂组成的阿基米德螺旋阵，传声器阵列上分布28个驻极体自由场传声器，传声器阵列布置于气流外安全区域，利用基于NI-PXIe总线的数据采集系统完成传声器阵列的多通道同步采集，采样频率44100Hz。

3.根据权利要求1所述的气流环境下的声源定位方法，其特征在于：所述步骤(二)中的阵列流形矩阵A考虑了气流影响造成的声线偏折问题，具体计算过程如下：

由Snell折射定律和速度矢量合成定理，建立气流环境下的声传播路径模型：

$\left\{\begin{array}{c}\frac{c}{cos\left(i_2\right)}+M_0=\frac{c}{cos\left(i_1\right)}\\ tan\left(i_3\right)=\frac{\sin \left(i_2\right)}{M_0+cos\left(i_2\right)}\\ \frac{c}{\cos \left(i_2\right)}+M_0=\frac{c}{\cos \left(i_5\right)}\end{array}\right.$

\*MERGEFORMAT(1)，

其中，M₀表示气流的马赫数，c为常温下声传播速度，i₁为气流外声波传入气流内部时初始入射方向，i₂为声波进入气流后折射方向，i₃为气流内矢量合成后的速度方向，i₅为声波传出气流时折射方向；

通过离散化和迭代方法，对式(1)进行求解，进而求得扫描点与传声器之间的近似声程r和声传播时延τ为：

$\left\{\begin{array}{c}r=L_2/\sin i_5+W_i/\sin i_3+L_1/\sin i_1\\ \mathrm{\τ}=(L_2/\sin i_5+L_1/\sin i_1)/c+(W_i/\sin i_3)/c_1\end{array}\right.$

\*MERGEFORMAT(2)，

其中，L₁为扫描声源面到气流一侧边缘距离，W_i为气流区域宽度，L₂为传声器阵列到气流另一侧边缘距离，c₁为矢量合成后的速度；

为求解三维问题，将各传声器单独取出，在传声器与扫描点单独组成的平面内进行传播时延的计算，以此来完成三维声传播时延的求解，获得阵列流行矩阵A：

$A=\[{\hat{a}}_1,...{\hat{a}}_l...,{\hat{a}}_L\],{\hat{a}}_l=\[e^{-{\mathrm{j\τ}}_{l,1}}/r_{l,1},e^{-{\mathrm{j\τ}}_{l,2}}/r_{l,2},...,e^{-{\mathrm{j\τ}}_{l,M}}/r_{l,M}\]$

\*MERGEFORMAT(3)。

4.根据权利要求1所述的气流环境下的声源定位方法，其特征在于：所述步骤(三)具体为：

对采集到的声信号数据y(b)求解其互谱矩阵G＝y(b)y^H(b)，然后对阵列信号互谱矩阵G进行特征分解G＝UΛU^H，将分解后得到的M个特征值从大到小排列，即λ₁＞λ₂＞...＞λ_M，令 组成的序列反映了互谱矩阵特征值之间的突变程度，预估声源个数为 取值应使得取前个中最大的特征值和对应的特征向量，构造降噪后的互谱矩阵利用二阶锥规划模型建立扫描点声压矩阵P、阵列流形矩阵A及降噪后的互谱矩阵间的价值函数(4)，通过最小化价值函数迭代求解声压矩阵P：

$\begin{array}{ccc}\underset{Q,X}{min}\left|\right|\stackrel{\‾}{G}{Q_n}^H-{\mathrm{AX}}_n|{|}_F^2,& s.t.& \left|\right|x_n|{|}_1\≤{\mathrm{\ξ}}_n,{Q_n}^H{Q}_n=I,n>1\end{array}$

\*MERGEFORMAT(4)，

式中，Q为一个任意酉矩阵，XX^H＝P，x为矩阵X中所有元素组成的向量，各变量的脚标n表示迭代次数，ξ_n为稀疏化约束参数。

5.根据权利要求1或4所述的气流环境下的声源定位方法，其特征在于：所述稀疏化约束条件||x_n||₁≤ξ_n约束了扫描声源面信号的总功率大小，稀疏化约束参数设置遵循以下原则：

当n＝1时，约束参数ξ_n的值可通过柯西-施瓦茨不等式(5)获得：

$\underset{l=1}{\overset{L\hat{L}}{\Σ}}\left|x_l\right|\≤\sqrt{\underset{l=1}{\overset{L\hat{L}}{\Σ}}|x_l{|}^{2}L\hat{L}}\≤\sqrt{\mathrm{\β}L\hat{L}}={\mathrm{\ξ}}_1$

\*MERGEFORMAT(5)，

式中，x_l为矩阵X中的元素，β＝tr(Λ)；

当n＞1时，约束参数ξ_n由式(6)实现自适应更新，达到稀疏化效果:

ξ_n+1＝||x_n||₁

\*MERGEFORMAT(6)。

6.根据权利要求1所述的气流环境下的声源定位方法，其特征在于：所述步骤(四)具体为通过凸优化循环方式最小化价值函数获得扫描点声压矩阵P，直至连续两次迭代求解结果小于设定的阈值时，迭代计算停止，利用所获得的矩阵P绘制扫描点声压重构云图，对比各重构点声压值大小确定声源位置。