CN104596636A - 声场分离方法 - Google Patents

Abstract

本发明旨在提供一种声场分离方法，包括以下步骤：A、测量测量面上的声压值B、根据振动体产生的声压值生成灰度图，根据实际声源点的位置计算出各点的灰度值，然后将各点灰度值作为各自的权值计算质心球虚拟源强球心坐标及非质心球虚拟源强球心坐标；C、对测量面上测定的声压值进行补零扩展，得到测量面扩展面上的声压值；D、建立测量面扩展面、测量面上声压值与多虚拟球等效源强声压之间的传递矩阵，获得测量面上的声压求解方程；E、结合、声压求解方程，获取两个声源面的声压值。该声场分离方法克服现有技术中波叠加法必须要求虚拟源强配置域与测量面共形的缺陷，同时还具有较高的计算效率和计算精度。

Description

声场分离方法

技术领域

本发明涉及噪声领域，具体涉及一种声场分离方法。

背景技术

在实际测量时，通常会遇到测量面两侧都有声源，或是测量面的一侧存在反射或散射的情况。这些干扰声源所产生的声场影响了对目标声源所辐射声场的准确测量。在实际工程中，为了更加准确的研究目标声源的声辐射特性或反射面的反射特性，需要将来自测量面两侧的辐射声分开。

现有分离方法包括：(1)基于Fourier变换的单面声场分离技术(《物理学报》(2009年58卷12期))。该方法能够去除奇异性，可以采用单面测量进行声场分离，该方法的缺点是受到傅里叶变换算法的影响，测量面孔径至少是重构面大小的数倍大时才能得到准确分离结果，因此测量量较大，无法对大型声源进行分离，即便能够得到声场分离结果，仍需要继续对声源面进行声场重构；(2)基于双面振速测量的声场分离技术(《声学学报》(2010年35卷6期))。该方法首先测量两个平行等间距测量面上的法向质点振速，再采用傅里叶法分离入射和辐射声场，该方法能够获得更高的法向质点振速精度，但却仍然受到测量孔径大小的限制；(3)基于声压和速度测量的统计最优声场分离方法。Jacobsen等人在《J.Acoust.Soc.Am》(2007年121卷3期)的文章中提出基于声压和速度测量的统计最优声场分离方法,该方法采用p-u声强探头，在单测量面上对局部声压和质点振速信息进行测量，再采用建立的联合求解公式实现来自测量面两侧的辐射声场分离和声源面的重建，该方法解决了测量空间大小的限制，但其缺点是计算量大、效率低、计算时间长。(4)基于球面波叠加的声场分离方法：1956年，J.Pachner采用球面波叠加法实现了任意波场中行波和驻波声场的分离，G.Weinreich等在1980年对J.Pachner提出的方法作了进一步改进，建立了基于双球面测量的声场分离理论。

现有声场分离技术仍然存在2个缺点：

(1)波叠加法必须要求虚拟源强配置域与测量面共形；

(2)计算量大、效率低、计算时间长，适应性不强。

发明内容

本发明旨在提供一种声场分离方法，该声场分离方法克服现有技术中波叠加法必须要求虚拟源强配置域与测量面共形的缺陷，同时还具有较高的计算效率和计算精度。

本发明的技术方案如下：一种声场分离方法，包括以下步骤：

A、在测量声场中两个声源面之间设置有测量面，测量面上呈网格设置声压振速传感器，相邻网格点之间的距离小于半个波长，声压振速传感器测量测量面上的声压值；

B、对振动体形状进行有限元仿真，建立多球形虚拟源，将多球形虚拟源分为质心球虚拟源强和非质心球虚拟源强两类，将声压图像的灰度值作为加权重心法的权重值，利用测量出的声压值生成灰度图，根据振动体上实际声源的位置坐标计算出各点的灰度值，然后将灰度值作为权值求系统重心，将系统重心作为质心球虚拟源强球心坐标，同时任选两个以上实际声源，利用各个实际声源的灰度值作为权值计算出与这些实际声源对应的非质心球虚拟源强球心坐标；

C、对测量面上测定的声压值进行补零扩展，得到测量面扩展面上的声压值；

D、建立测量面扩展面声压值与多虚拟球等效源强声压之间的传递矩阵，建立测量面上的声压值与多虚拟球等效源强声压之间的传递矩阵，对测量面扩展面声压值与多虚拟球等效源强声压之间的传递矩阵进行正则化处理，获得测量面上的声压求解方程，求解测量面上的声压值；

E、利用步骤D获得的声压求解方程，将两侧声源的声场分离，获取两个声源面的声压值。

所述步骤B具体步骤如下：

B1、对振动体形状进行有限元仿真，根据振动体产生的声压值生成灰度图，根据实际声源点的位置计算出各个声源的灰度值f(x_i,y_i,z_i)，(x_i,y_i,z_i)为各实际声源点的的坐标，设置的多球形虚拟源分为质心球虚拟源强和非质心球虚拟源强两类；

B2、利用式(1)算出整个系统的重心作为质心球虚拟源强球心坐标，

${\stackrel{\‾}{x}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\Σ}}}{x}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\Σ}}}f(x_i,y_i,z_i)};{\stackrel{\‾}{y}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\Σ}}}{y}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\Σ}}}f(x_i,y_i,z_i)};{\stackrel{\‾}{z}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\Σ}}}{z}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\Σ}}}f(x_i,y_i,z_i)}---\left(1\right);$

其中v为不规则振动体内部声源的个数；

B3、从振动体声源中任选w个声源，参照式(2)计算出一定数量的非质心球虚拟源强球心坐标其中2≤w<v；

${\stackrel{\&OverBar;}{x}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{x}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{y}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{y}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{z}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{z}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};1\&le;j<w---\left(2\right).$

所述步骤C具体步骤如下：

依照式(3)对测量面上的声压值进行补零扩展得到测量面扩展面上的声压值：

p_E(H⁺)＝D·p_E(H)    (3)；

其中，p_E(H)为测量面H上的声压测量值；p_E(H⁺)为测量面的扩展面上的声压值；D＝diag[D₁₁,D₂₂,…,D_NN]， $D_{\mathrm{ii}}=\left\{\begin{array}{cc}1& (y,z)\&Element;H\\ 0& (y,z)\&NotElement;H\end{array}\right.;$ D为采样算子，D_ii为矩阵D上对角线上的值；N为测量面扩展面上网格的点数；测量面H(x＝x_H)与(y,z)坐标面平行，测量面的法向为x方向。

所述步骤D具体步骤如下：

D1、建立测量面扩展面声压值与虚拟球等效源声压值之间的传递矩阵如式(4)：

p_E(H⁺)＝[T_H]Q    (4)；

其中p_E(H⁺)为测量面扩展面上的声压值，Q为虚拟球等效源强双向傅里叶分解后的系数矩阵；[T_H]为虚拟球等效源强与测量面扩展面上声压之间的传递矩阵；

其中[T_H]由以下步骤得到：

根据等效源强理论，假设S′是振动体内某一虚拟源强分布表面，则可将外场中某点r处的声压表示为

$p\left(r\right)=\underset{t=1}{\overset{t}{\mathrm{\&Sigma;}}}{\&Integral;}_{S^{\&prime;}}\mathrm{\&sigma;}\left(r_Q\right)K(r,r_Q){\mathrm{dS}}^{\&prime;}---\left(5\right);$

式中，r_Q是虚源面上的某一点，t为多虚拟球的个数，σ(r_Q)为待求的源强密度函数，K(r,r_Q)为积分核函数，K(r,r_Q)＝g(r,r_Q)＝(1/4πR)e^ikR，k＝ω/c为波数,R为两点间的距离；

将未知源强密度函数进行双向Fourier级数展开，并利用二维快速FFT计算积分，采用梯形公式离散格林函数，将与距离有关的矩阵规整成一个矩阵，即为[T_H]；

D2、建立测量面声压值与虚拟球等效源声压值之间的传递矩阵如下：

{p_E}＝[T_E]Q    (6)；

其中，{p_E}为测量面H上的声压值，[T_E]为虚拟球等效源强与测量面上声压之间的传递矩阵，[T_E]的计算参照[T_H]的计算进行；Q为虚拟球等效源强双向傅里叶分解后的系数矩阵；

D3、联立式(4)、(6)，求解测量面上的声压：

{p_E}＝[T_E][T_H]^-1p_E(H⁺)    (7)；

对T_H进行正则化处理，得到：

$\left\{p_E\right\}=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H\left\{p_{H^{+}}\right\}---\left(8\right);$

其中，α为正则化参数，I为单位对角矩阵，为T_H的共轭转置矩阵，为 的逆矩阵。

所述步骤E具体步骤如下：

测量面H上的声压和法向质点振速为：

p_E(H)＝p_1E(H)+p_2E(H)；v_E(H)＝v_1E(H)-v_2E(H)    (9)；

式中p_1E(H)为声源1在测量面H上所辐射的声压、p_2E(H)为声源2在测量面H上所辐射的声压、v_1E(H)为声源1在测量面H上所辐射的法向质点振速、v_2E(H)为声源2在测量面H上所辐射的法向质点振速；

对式(9)进行补零扩展得到：

p_E(H⁺)＝p_1E(H⁺)+p_2E(H⁺)；v_E(H⁺)＝v_1E(H⁺)-v_2E(H⁺)    (10)；

与式(8)相结合，得到

$p_1=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_{1E}\left(H^{+}\right)---\left(11\right);$

$p_2=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_{2E}\left(H^{+}\right)---\left(12\right);$

$p_1=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_{1E}\left(H^{+}\right)---\left(13\right);$

$p_2=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_{2E}\left(H^{+}\right)---\left(14\right);$

其中p₁为声源面S₁上的声压、p₂为声源面S₂上的声压,W_E为虚拟球等效源强与测量面上质点振速之间的传递矩阵，W_H为虚拟球等效源强与测量面扩展面上质点振速之间的传递矩阵；

将式(11)和式(12)相加，将式(13)和式(14)相减，可得：

$\begin{array}{c}{p}_1+p_2=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H(p_{1E}\left(H^{+}\right)+p_{2E}\left(H^{+}\right))\\ =T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)\end{array}---\left(15\right);$

$\begin{array}{c}{p}_1-p_2=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H(v_{1E}\left(H^{+}\right)-v_{2E}\left(H^{+}\right))\\ =W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right)\end{array}---\left(16\right);$

将式17和式18相加，可得：

$p_1=\frac{1}{2}(T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)+W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right))---\left(17\right);$

将式17和式18相减，可得：

$p_2=\frac{1}{2}(T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)-W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right))---\left(18\right).$

本发明技术方案采用波叠加法进行声场分离，并采用灰度重心法对虚拟源强进行配置，灰度重心法通过将各实际声源点的灰度值作为权值来计算相应的虚拟源坐标，不仅克服了常规波叠加法必须要求虚拟源强配置域与测量面共形的缺陷，并且通过灰度权值进行配置的虚拟源强更贴近实际声源点的分布情况，提升声场的重建与分离的准确度；同时，通过将多虚拟球波叠加法与多虚拟域的二维快速Fourier变换算法相结合，极大地改善了计算精度、计算速度和适应性。

附图说明

图1为本发明声场分离方法的步骤流程图

图2为本发明声场分离原理图

图3为多球形虚拟源强配置面与测量面的布置示意图

图4为实施例1振动体内某截面声压幅值分布图的灰度图

图5为实施例1重建声场声压幅值图

图6为实施例1重建声场声压理论幅值图

图7为实施例1声压幅值相对误差图

图8为实施例1声源1声压幅值图

图9为实施例1声源1声压幅值相对误差图

图10为实施例1声源2声压幅值图

图11为实施例1声源2声压幅值相对误差图

图3中各部分名称及标号如下:

1为测量平面，2为传声器，3为虚拟球等效源。

具体实施方式

下面结合具体实施例说明本发明。

实施例1

如图1所示，本实施例声场分离方法包括以下步骤：

A、在测量声场中两个声源面之间设置有测量面，测量面上呈网格设置声压振速传感器，相邻网格点之间的距离小于半个波长，声压振速传感器测量测量面上的声压值；

B、对振动体形状进行有限元仿真，建立多球形虚拟源，将多球形虚拟源分为质心球虚拟源强和非质心球虚拟源强两类，将声压图像的灰度值作为加权重心法的权重值，利用测量出的声压值生成灰度图，根据振动体上实际声源的位置坐标计算出各点的灰度值，然后将灰度值作为权值求系统重心，将系统重心作为质心球虚拟源强球心坐标，同时任选两个以上实际声源，利用各个实际声源的灰度值作为权值计算出与这些实际声源对应的非质心球虚拟源强球心坐标；

具体步骤如下：

B1、对振动体形状进行有限元仿真，根据振动体产生的声压值生成灰度图，根据实际声源点的位置计算出各个声源的灰度值f(x_i,y_i,z_i)，(x_i,y_i,z_i)为各实际声源点的的坐标，设置的多球形虚拟源分为质心球虚拟源强和非质心球虚拟源强两类；

B2、利用式(1)算出整个系统的重心作为质心球虚拟源强球心坐标，

${\stackrel{\&OverBar;}{x}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}{x}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{y}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}{y}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{z}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}{z}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)}---\left(1\right);$

其中v为不规则振动体内部声源的个数；

B3、从振动体声源中任选w个声源，参照式(2)计算出一定数量的非质心球虚拟源强球心坐标其中2≤w<v；

${\stackrel{\&OverBar;}{x}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{x}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{y}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{y}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{z}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{z}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};1\&le;j<w---\left(2\right);$

C、对测量面上测定的声压值进行补零扩展，得到测量面扩展面上的声压值；

所述步骤C具体步骤如下：

依照式(3)对测量面上的声压值进行补零扩展得到测量面扩展面上的声压值：

p_E(H⁺)＝D·p_E(H)    (3)；

其中，p_E(H)为测量面H上的声压测量值；p_E(H⁺)为测量面的扩展面上的声压值；D＝diag[D₁₁,D₂₂,…,D_NN]， $D_{\mathrm{ii}}=\left\{\begin{array}{cc}1& (y,z)\&Element;H\\ 0& (y,z)\&NotElement;H\end{array}\right.;$ D为采样算子，D_ii为矩阵D上对角线上的值；N为测量面扩展面上网格的点数；测量面H(x＝x_H)与(y,z)坐标面平行，测量面的法向为x方向；

D、建立测量面扩展面声压值与多虚拟球等效源强声压之间的传递矩阵，建立测量面上的声压值与多虚拟球等效源强声压之间的传递矩阵，对测量面扩展面声压值与多虚拟球等效源强声压之间的传递矩阵进行正则化处理，获得测量面上的声压求解方程，求解测量面上的声压值；

所述步骤D具体步骤如下：

D1、建立测量面扩展面声压值与虚拟球等效源声压值之间的传递矩阵如式(4)：

p_E(H⁺)＝[T_H]Q    (4)；

其中p_E(H⁺)为测量面扩展面上的声压值，Q为虚拟球等效源强双向傅里叶分解后的系数矩阵；[T_H]为虚拟球等效源强与测量面扩展面上声压之间的传递矩阵；

其中[T_H]由以下步骤得到：

根据等效源强理论，假设S′是振动体内某一虚拟源强分布表面，则可将外场中某点r处的声压表示为

$p\left(r\right)=\underset{t=1}{\overset{t}{\mathrm{\&Sigma;}}}{\&Integral;}_{S^{\&prime;}}\mathrm{\&sigma;}\left(r_Q\right)K(r,r_Q){\mathrm{dS}}^{\&prime;}---\left(5\right);$

式中，r_Q是虚源面上的某一点，t为多虚拟球的个数，σ(r_Q)为待求的源强密度函数，K(r,r_Q)为积分核函数，K(r,r_Q)＝g(r,r_Q)＝(1/4πR)e^ikR，k＝ω/c为波数,R为两点间的距离；

将未知源强密度函数进行双向Fourier级数展开，并利用二维快速FFT计算积分，采用梯形公式离散格林函数，将与距离有关的矩阵规整成一个矩阵，即为[T_H]；

D2、建立测量面声压值与虚拟球等效源声压值之间的传递矩阵如下：

{p_E}＝[T_E]Q    (6)；

其中，{p_E}为测量面H上的声压值，[T_E]为虚拟球等效源强与测量面上声压之间的传递矩阵，[T_E]的计算参照[T_H]的计算进行；Q为虚拟球等效源强双向傅里叶分解后的系数矩阵；

D3、联立式(4)、(6)，求解测量面上的声压：

{p_E}＝[T_E][T_H]^-1p_E(H⁺)    (7)；

对T_H进行正则化处理，得到：

$\left\{p_E\right\}=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H\left\{p_{H^{+}}\right\}---\left(8\right);$

其中，α为正则化参数，I为单位对角矩阵，为T_H的共轭转置矩阵，为 的逆矩阵；

E、利用步骤D获得的声压求解方程，将两侧声源的声场分离，获取两个声源面的声压值；

所述步骤E具体步骤如下：

测量面H上的声压和法向质点振速为：

p_E(H)＝p_1E(H)+p_2E(H)；v_E(H)＝v_1E(H)-v_2E(H)    (9)；

式中p_1E(H)为声源1在测量面H上所辐射的声压、p_2E(H)为声源2在测量面H上所辐射的声压、v_1E(H)为声源1在测量面H上所辐射的法向质点振速、v_2E(H)为声源2在测量面H上所辐射的法向质点振速；

对式(9)进行补零扩展得到：

p_E(H⁺)＝p_1E(H⁺)+p_2E(H⁺)；v_E(H⁺)＝v_1E(H⁺)-v_2E(H⁺)    (10)；

与式(8)相结合，得到

$p_1=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_{1E}\left(H^{+}\right)---\left(11\right);$

$p_2=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_{2E}\left(H^{+}\right)---\left(12\right);$

$p_1=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_{1E}\left(H^{+}\right)---\left(13\right);$

$p_2=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_{2E}\left(H^{+}\right)---\left(14\right);$

其中p₁为声源面S₁上的声压、p₂为声源面S₂上的声压,W_E为虚拟球等效源强与测量面上质点振速之间的传递矩阵，W_H为虚拟球等效源强与测量面扩展面上质点振速之间的传递矩阵；

将式(11)和式(12)相加，将式(13)和式(14)相减，可得：

$\begin{array}{c}{p}_1+p_2=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H(p_{1E}\left(H^{+}\right)+p_{2E}\left(H^{+}\right))\\ =T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)\end{array}---\left(15\right);$

$\begin{array}{c}{p}_1-p_2=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H(v_{1E}\left(H^{+}\right)-v_{2E}\left(H^{+}\right))\\ =W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right)\end{array}---\left(16\right);$

将式17和式18相加，可得：

$p_1=\frac{1}{2}(T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)+W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right))---\left(17\right);$

将式17和式18相减，可得：

$p_2=\frac{1}{2}(T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)-W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right))---\left(18\right).$

如图2所示，本实施例设置位于测量面的两侧为两个振动体S₁、S₂，其中S₁为36个点声源组成的振动体，S₂为9个点声源组成的振动体，都分布在0.4mх0.4mх0.2m的长方体内，以测量面的中心为坐标原点，测量面指向声源S₂方向为x轴正方向，多球形虚拟源强配置面与测量面的布置如图3所示；

本实施例采用灰度重心法布置多虚拟球的位置：根据测量出的声压值生成灰度图，振动体某截面声压幅值分布图的灰度图中灰度与该灰度点数量的关系如图4所示，根据实际声源点的位置计算出各点的灰度值，然后将各点灰度值作为各自的权值计算系统重心，进而作为质心球虚拟源强球心坐标；选择25个实际声源点，以其各自灰度值作为各自权值，计算出X个非质心球虚拟源强球心坐标，25个实际声源点灰度值如表1所示；

表1 选取的实际声源灰度值统计表

其中：对于质心位置的虚拟球的半径为r₀＝0.025m；非质心位置的多虚拟球的半径r₀＝0.05m；

本实施例重建的声场声压幅值如图5所示，进行声场分离后的声源1声压幅值、声源2声压幅值分别如图8、10所示；

在自由空间任意一场点P的理论声压表示为：

$p\left(P\right)=\mathrm{i\&rho;ck}\underset{n=1}{\overset{n}{\mathrm{\&Sigma;}}}{Q}_n\frac{\exp (-{\mathrm{ikr}}_n)}{{4\mathrm{\&pi;r}}_n}---\left(19\right);$

式中，r_n为场点P到第n个点源的距离，Q_n为第n个点源的强度；ρ为大气密度，c为声速，k＝2πf/c为波数；

声压幅值、幅值相对误差分别为：

$\left|p\left(r\right)\right|=\sqrt{{\left(\mathrm{Re}\left(p\left(r\right)\right)\right)}^2+{\left(\mathrm{Im}\left(p\left(r\right)\right)\right)}^2}---\left(20\right);$

${\mathrm{\&eta;}}_p=\frac{|p_T\left(r\right)-p_F\left(r\right)|}{\left|p_T\left(r\right)\right|}\&times;100\%---\left(21\right);$

式中p_T(r)为各点理论声压幅值，p_F(r)为本实施例重建的各点声压幅值；

根据式(20)算出测量面处的理论声压幅值如图6所示；

根据式(21)算出测量面处的理论声压幅值相对误差如图7所示；

根据式(20)、(21)结合声源1声压幅值算出声源1声压幅值相对误差如图9所示；

根据式(20)、(21)结合声源2声压幅值算出声源2声压幅值相对误差如图11所示。

由上述图示可以看出，本实施例的声场分离方法的分离出的声源声压幅值相对误差处于可接受的范围内，验证了本实施例声场分离方法的准确性。

Claims (5)

1.一种声场分离方法，其特征在于包括以下步骤：

A、在测量声场中两个声源面之间设置有测量面，测量面上呈网格设置声压振速传感器，相邻网格点之间的距离小于半个波长，声压振速传感器测量测量面上的声压值；

B、对振动体形状进行有限元仿真，建立多球形虚拟源，将多球形虚拟源分为质心球虚拟源强和非质心球虚拟源强两类，将声压图像的灰度值作为加权重心法的权重值，利用测量出的声压值生成灰度图，根据振动体上实际声源的位置坐标计算出各点的灰度值，然后将灰度值作为权值求系统重心，将系统重心作为质心球虚拟源强球心坐标，同时任选两个以上实际声源，利用各个实际声源的灰度值作为权值计算出与这些实际声源对应的非质心球虚拟源强球心坐标；

C、对测量面上测定的声压值进行补零扩展，得到测量面扩展面上的声压值；

D、建立测量面扩展面声压值与多虚拟球等效源强声压之间的传递矩阵，建立测量面上的声压值与多虚拟球等效源强声压之间的传递矩阵，对测量面扩展面声压值与多虚拟球等效源强声压之间的传递矩阵进行正则化处理，获得测量面上的声压求解方程，求解测量面上的声压值；

E、利用步骤D获得的声压求解方程，将两侧声源的声场分离，获取两个声源面的声压值。

2.如权利要求1所述的声场分离方法，其特征在于所述步骤B具体步骤如下：

B1、对振动体形状进行有限元仿真，根据振动体产生的声压值生成灰度图，根据实际声源点的位置计算出各个声源的灰度值f(x_i,y_i,z_i)，(x_i,y_i,z_i)为各实际声源点的的坐标，设置的多球形虚拟源分为质心球虚拟源强和非质心球虚拟源强两类；

B2、利用式(1)算出整个系统的重心作为质心球虚拟源强球心坐标，

${\stackrel{\&OverBar;}{x}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}{x}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{y}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}{y}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{z}}_0=\frac{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}{z}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{v}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)}---\left(1\right);$

其中v为不规则振动体内部声源的个数；

B3、从振动体声源中任选w个声源，参照式(2)计算出一定数量的非质心球虚拟源强球心坐标其中2≤w＜v；

${\stackrel{\&OverBar;}{x}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{x}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{y}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{y}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};{\stackrel{\&OverBar;}{z}}_j=\frac{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}{z}_{i}f(x_i,y_i,z_i)}{\underset{i=0}{\overset{w}{\mathrm{\&Sigma;}}}f(x_i,y_i,z_i)};1\&le;j<w---\left(2\right).$

3.如权利要求1所述的声场分离方法，其特征在于所述步骤C具体步骤如下：

所述步骤C具体步骤如下：

依照式(3)对测量面上的声压值进行补零扩展得到测量面扩展面上的声压值：

p_E(H⁺)＝D·p_E(H)    (3)；

其中，p_E(H)为测量面H上的声压测量值；p_E(H⁺)为测量面的扩展面上的声压值；D＝diag[D₁₁,D₂₂,...,D_NN]， $D_{\mathrm{ii}}=\left\{\begin{array}{cc}1& (y,z)\&Element;H\\ 0& (y,z)\&NotElement;H\end{array}\right.;$ D为采样算子，D_ii为矩阵D上对角线上的值；N为测量面扩展面上网格的点数；测量面H(x＝x_H)与(y,z)坐标面平行，测量面的法向为x方向。

4.如权利要求3所述的声场分离方法，所述步骤D具体步骤如下：

D1、建立测量面扩展面声压值与虚拟球等效源声压值之间的传递矩阵如式(4)：

p_E(H⁺)＝[T_H]Q    (4)；

其中p_E(H⁺)为测量面扩展面上的声压值，Q为虚拟球等效源强双向傅里叶分解后的系数矩阵；[T_H]为虚拟球等效源强与测量面扩展面上声压之间的传递矩阵；

其中[T_H]由以下步骤得到：

根据等效源强理论，假设S′是振动体内某一虚拟源强分布表面，则可将外场中某点r处的声压表示为

$p\left(r\right)=\underset{t=1}{\overset{t}{\mathrm{\&Sigma;}}}{\&Integral;}_{S^{\&prime;}}\mathrm{\&sigma;}\left(r_Q\right)K(r,r_Q){\mathrm{dS}}^{\&prime;}---\left(5\right);$

式中，r_Q是虚源面上的某一点，t为多虚拟球的个数，σ(r_Q)为待求的源强密度函数，K(r,r_Q)为积分核函数，K(r,r_Q)＝g(r,r_Q)＝(1/4πR)e^ikR，k＝ω/c为波数,R为两点间的距离；

将未知源强密度函数进行双向Fourier级数展开，并利用二维快速FFT计算积分，采用梯形公式离散格林函数，将与距离有关的矩阵规整成一个矩阵，即为[T_H]；

D2、建立测量面声压值与虚拟球等效源声压值之间的传递矩阵如下：

{p_E}＝[T_E]Q    (6)；

其中，{p_E}为测量面H上的声压值，[T_E]为虚拟球等效源强与测量面上声压之间的传递矩阵，[T_E]的计算参照[T_H]的计算进行；Q为虚拟球等效源强双向傅里叶分解后的系数矩阵；

D3、联立式(4)、(6)，求解测量面上的声压：

{p_E}＝[T_E][T_H]^-1p_E(H⁺)    (7)；

对T_H进行正则化处理，得到：

$\left\{p_E\right\}=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H\left\{p_{H^{+}}\right\}---\left(8\right);$

其中，α为正则化参数，I为单位对角矩阵，为T_H的共轭转置矩阵，为的逆矩阵。

5.如权利要求4所述的声场分离方法，其特征在于所述步骤E具体步骤如下：

测量面H上的声压和法向质点振速为：

p_E(H)＝p_1E(H)+p_2E(H)；v_E(H)＝v_1E(H)-v_2E(H)    (9)；

式中p_1E(H)为声源1在测量面H上所辐射的声压、p_2E(H)为声源2在测量面H上所辐射的声压、v_1E(H)为声源1在测量面H上所辐射的法向质点振速、v_2E(H)为声源2在测量面H上所辐射的法向质点振速；

对式(9)进行补零扩展得到：

p_E(H⁺)＝p_1E(H⁺)+p_2E(H⁺)；v_E(H⁺)＝v_1E(H⁺)-v_2E(H⁺)    (10)；

与式(8)相结合，得到

$p_1=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_{1E}\left(H^{+}\right)---\left(11\right);$

$p_2=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_{2E}\left(H^{+}\right)---\left(12\right);$

$p_1=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_{1E}\left(H^{+}\right)---\left(13\right);$

$p_2=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_{2E}\left(H^{+}\right)---\left(14\right);$

其中p₁为声源面S₁上的声压、p₂为声源面S₂上的声压,W_E为虚拟球等效源强与测量面上质点振速之间的传递矩阵，W_H为虚拟球等效源强与测量面扩展面上质点振速之间的传递矩阵；

将式(11)和式(12)相加，将式(13)和式(14)相减，可得：

$\begin{array}{c}{p}_1+p_2=T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H(p_{1E}\left(H^{+}\right)+p_{2E}\left(H^{+}\right))\\ =T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)\end{array}---\left(15\right);$

$\begin{array}{c}{p}_1-p_2=W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H(v_{1E}\left(H^{+}\right)-v_{2E}\left(H^{+}\right))\\ =W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right)\end{array}---\left(16\right);$

将式17和式18相加，可得：

$p_1=\frac{1}{2}(T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)+W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right))---\left(17\right);$

将式17和式18相减，可得：

$p_2=\frac{1}{2}(T_E{(\mathrm{\&alpha;I}+T_H^H{T}_H)}^{-1}{T}_H^H{p}_E\left(H^{+}\right)-W_E{(\mathrm{\&alpha;I}+W_H^H{W}_H)}^{-1}{W}_H^H{v}_E\left(H^{+}\right))---\left(18\right).$