CN103712684A - 声场重建方法 - Google Patents

Abstract

本发明旨在提供一种声场重建方法，包括以下步骤A、测量测量声场中测量面上的声压值；B、根据声源的实际位置和声强的大小，利用重心法对于不同的声源赋以一定的权重，计算出多球虚拟源强的位置、个数以及每个虚拟球的位置坐标；C、对测量面上测定的声压值进行补零扩展，得到测量面扩展面上的声压值；D、计算测量面扩展面声压值与虚拟球声压之间的传递矩阵；E、建立重建面上的声压值与多虚拟球等效源强声压值之间的传递矩阵；F、联立步骤D和E获得的传递矩阵，求解重建面上的声压值。本发明声场重建方法通过引入多球面域二维快速Fourier变换算法，不需虚拟源强配置域与测量面共形，具有很高计算效率和计算精度。

Description

声场重建方法

技术领域

本发明涉及噪声源识别方法领域，具体涉及一种采用声压测量和基于多虚拟球的patchNAH声场重建方法。

背景技术

Patch近场声全息(Patch near-field acoustic holography，PNAH)是在近场声全息理论基础上发展形成的一种新型近场声全息技术。PNAH允许测量面小于实际源面,重建区域可以是局部源面，突破了常规NAH对全息孔径尺寸的限制，更便于NAH的实际应用。

现有局部近场声全息技术包括：（1）数据外推方法。Yoshikawa在《美国声学学会》（2001年110卷4期）的文章中对数据外推技术进行了研究，提出时域和波数域的两种数据外推方法，该方法的缺点是需要进行迭代计算，计算过程复杂计算量大；（2）Patch近场声全息技术。Williams在《美国声学学会》（2003年113卷3期）的文章中通过对近场概念的进一步研究，提出了一种Patch NAH的概念，这种方法允许测量的传感器阵列比声源小，克服了测量孔径的有限大小引起的重建误差，降低了测量时间，但在计算时间方面该技术仍然受到迭代计算的限制，计算过程复杂计算，效率低；（3）2005年，Lee等将PNAH推广应用到柱面坐标系。SAIJYOU等采用边界元法对Patch原理的重建思想进行推广，提出了基于BEM的PNAH。SARKISSIAN等利用波叠加法(Wave superposition approach，WSA)实现了全息声压的近场外推，建立了基于WSA的PNAH。

等效源法是由Koopmann等最先提出,该方法又被称为波叠加方法、源模拟法等。它的基本原理是用一系列定位的基本源（即等效源），来代替实际声源。通过测量声压、声速等参数，在一系列等效源产生的合成值之间联系的相容性条件上，得到声源的理论分布情况。该方法避免了边界元法存在的复杂插值运算、奇异积分处理和特征波数处解的非唯一性等问题,具有运算速度快、计算精度高、求解稳定性好等优点。

然而，现有波叠加声全息技术仍然存在2个缺点：（1）虚拟源强的分布域（面）需与测量面（全息面）共形；（2）采用离散源强或基于边界元的分布源强，计算效率和精度低，适应性不强。

发明内容

本发明旨在提供一种声场重建方法，该声场重建方法克服现有技术虚拟源强的分布域需与测量面共形、计算效率与精度低的缺陷，通过引入多球面域二维快速Fourier变换算法，不需要虚拟源强配置域与测量面共形，具有很高计算效率和计算精度。

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

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

B、对振动体形状进行有限元仿真，根据声源的实际位置和声强的大小，利用重心法对于不同的声源赋以一定的权重，计算出多球虚拟源强的位置和个数，并生成每个虚拟球的位置坐标；

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

为了减少小全息孔径造成的误差，对在较小孔径内测得的声压数据进行外推，获得较大测量孔径内的声压数据的近似值。依照式（1）对测量面上的声压测量数据进行补零扩展，得到测量面扩展面上的声压值。

D、计算测量面扩展面声压值与虚拟球声压之间的传递矩阵；

E、建立重建面上的声压值与多虚拟球等效源强产生的声压值之间的传递矩阵；

F、联立步骤D和E获得的传递矩阵，同时对步骤D获得的传递矩阵进行正则化方法处理，获得重建面上的声压求解方程，求解重建面上的声压值。

所述步骤B中，将虚拟球区分为质心位置的虚拟球和非质心位置的虚拟球，其中对于质心位置的虚拟球的半径为r₀=0.2；非质心位置的各虚拟球的半径为r₀=0.1。

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

依照式（1）对测量面上的声压值进行补零扩展得到测量面扩展面上的声压值：

p_E(H⁺)=D·p_E(H)   （1）；

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

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

测量面扩展面声压值与虚拟球等效源声压值之间的传递矩阵如式（2）：

p_E(H⁺)=[T_H]Q   （2）；

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

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

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

$p\left(r\right)=\underset{t=1}{\overset{t}{\mathrm{\Σ}}}{\∫}_{S^{\′}}\mathrm{\σ}\left(r_Q\right)K(r,r_Q){\mathrm{dS}}^{\′}---\left(3\right);$

式中，r_Q是虚源面上的某一点，t为多虚拟球的个数，σ(r_Q)为待求的源强密度函数，K(r,r_Q)为积分核函数，K(r,r_Q)=g(r,r_Q)=(14πR)e^ikR将未知源强密度函数进行双向Fourier级数展开，并利用二维快速FFT计算积分，采用梯形公式离散格林函数，将与距离有关的矩阵规整成一个矩阵，即为[T_H]。

所述步骤E中重建面声压值与虚拟球等效源声压值之间的传递矩阵如下：

{p_R}=[T_R]Q   （4）；

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

所述步骤F中重建面声压值的求解具体步骤如下：

联立式（2）、（4），求解重建面上的声压：

{p_R}=[T_R][T_H]^-1p_E(H⁺)   （5）；

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

$\left\{p_R\right\}=T_R{(\mathrm{\αI}+T_H^H{T}_H)}^{-1}{T}_H^H\left\{p_{H^{+}}\right\}---\left(6\right);$

其中，α为正则化参数，I为单位对角矩阵，

为T_H的共轭转置矩阵，

为

的逆矩阵。

本发明声场重建方法的优点如下：

1.采用多球形虚拟源强，通过合理配置多虚拟球的位置和个数，克服了虚拟源强的分布域需与测量面共形问题，同时对测量面进行了补零扩展，对在较小孔径内测得的声压数据进行外推，获得较大测量孔径内的声压数据的近似值，减少小全息孔径造成的误差；

2.通过建立多虚拟球波叠加法，并采用多虚拟域的二维快速Fourier变换算法，两者相结合，极大地改善了计算精度、计算速度和适应性。

附图说明

图1为本发明声场重建方法的步骤流程图

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

图3为测量面补零扩展后的声压幅值图

图4为测量面补零扩展后的初相位图

图5为本实施例方法重建声压幅值图

图6为本实施例方法重建初相位图

图7为声压幅值的理论值图

图8为初相位理论值图

图9为幅值误差图

图10为初相位误差图

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

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

具体实施方式

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

实施例1

如图1所示，本实施例声场重建方法的步骤如下：

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

B、对振动体形状进行有限元仿真，根据声源的实际位置和声强的大小，利用重心法对于不同的声源赋以一定的权重，计算出多球虚拟源强的位置和个数，并生成每个虚拟球的位置坐标；

其中，将虚拟球区分为质心位置的虚拟球和非质心位置的虚拟球，其中对于质心位置的虚拟球的半径为r₀=0.2；非质心位置的各虚拟球的半径为r₀=0.1；

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

依照式（1）对测量面上的声压值进行补零扩展得到测量面扩展面上的声压值：

p_E(H⁺)=D·p_E(H)   （1）；

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

D、计算测量面扩展面声压值与虚拟球声压之间的传递矩阵，具体如下：

测量面扩展面声压值与虚拟球等效源声压值之间的传递矩阵如式（2）：

p_E(H⁺)=[T_H]Q   （2）；

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

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

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

$p\left(r\right)=\underset{t=1}{\overset{t}{\mathrm{\Σ}}}{\∫}_{S^{\′}}\mathrm{\σ}\left(r_Q\right)K(r,r_Q){\mathrm{dS}}^{\′}---\left(3\right);$

式中，r_Q是虚源面上的某一点，t为多虚拟球的个数，σ(r_Q)为待求的源强密度函数，K(r,r_Q)为积分核函数，K(r,r_Q)=g(r,r_Q)=(14πR)e^ikR将未知源强密度函数进行双向Fourier级数展开，并利用二维快速FFT计算积分，采用梯形公式离散格林函数，将与距离有关的矩阵规整成一个矩阵，即为[T_H]；

E、建立重建面上的声压值与多虚拟球等效源强产生的声压值之间的传递矩阵，具体如下：

所述步骤E中重建面声压值与虚拟球等效源声压值之间的传递矩阵如下：

{p_R}=[T_R]Q   （4）；

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

F、联立步骤D和E获得的传递矩阵，同时对步骤D获得的传递矩阵进行正则化方法处理，获得重建面上的声压求解方程，求解重建面上的声压值，具体如下:

联立式（2）、（4），求解重建面上的声压：

{p_R}=[T_R][T_H]^-1p_E(H⁺)   （5）；

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

$\left\{p_R\right\}=T_R{(\mathrm{\αI}+T_H^H{T}_H)}^{-1}{T}_H^H\left\{p_{H^{+}}\right\}---\left(6\right);$

其中，α为正则化参数，I为单位对角矩阵，

为T_H的共轭转置矩阵，

为

的逆矩阵。

本实施例假设一振动体的声源为六个点源，分布在0.4mх0.4mх0.2m的长方体的六个面的中心，以测量面的中心为坐标原点，测量面指向声源方向为x轴正方向，长方体内腔内任意位置分布着16个点源；

本实施例采用重心法布置多虚拟球的位置和个数，利用有限元软件计算出它们的重心坐标，本实施例虚拟球个数为9个；各虚拟球半径设置如下：

A.对于质心位置的虚拟球的半径为r₀=0.2；

B.对于非质心位置的多虚拟球的半径r₀=0.1；

本实施例的测量面采用的是6*6的声阵列布置方式，测量面区域H为-0.2≤y≤0.2，-0.2≤z≤0.2；补零扩展面H⁺为-0.4≤y≤0.4，-0.4≤z≤0.4；

图2为多球形虚拟源强配置面与测量面的布置示意图，图3和图4分别为测量面补零扩展后的声压幅值和初相位分布图；

本实施例的重建面与测量面的距离x为0.25m，重建面的区域为-0.2≤y≤0.2，-0.2≤z≤0.2；基于上述方法步骤重建出来的声压幅值和初相位分别如图5，图6所示；

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

$p\left(P\right)=\mathrm{i\ρck}\underset{n=1}{\overset{16}{\mathrm{\Σ}}}{Q}_n\frac{\exp (-{\mathrm{ikr}}_n)}{{4\mathrm{\πr}}_n}$

式中，r_n为场点P到第n个点源的距离，Q_n为第n个点源的强度；根据上述公式算出x=0.25m处的理论声压幅值和理论初相位分别如图7和图8所示；

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

$\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},{\mathrm{\η}}_p=\frac{|p_T\left(r\right)-p_F\left(r\right)|}{\left|P_T\left(r\right)\right|}\×100\%,$

$\mathrm{\Φ}\left(r\right)={\tan }^{-1}(\left(\mathrm{Re}\left(p\left(r\right)\right)\right)/\mathrm{Im}\left(p\left(r\right)\right)),{\mathrm{\η}}_{\mathrm{\Φ}}=\frac{|{\mathrm{\Φ}}_T\left(r\right)-{\mathrm{\Φ}}_F\left(r\right)|}{\left|{\mathrm{\Φ}}_T\left(r\right)\right|}100\%,$

p_T(r)、Φ_T(r)分别为各点理论声压幅值、初相位，p_F(r)、Φ_F(r)分别为本实施例重建的各点声压幅值、初相位，声压幅值相对误差与相位相对误差分别如图9、10所示。

Claims (6)

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

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

B、对振动体形状进行有限元仿真，根据声源的实际位置和声强的大小，利用重心法对于不同的声源赋以一定的权重，计算出多球虚拟源强的位置和个数，并生成每个虚拟球的位置坐标；

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

D、计算测量面扩展面声压值与虚拟球声压之间的传递矩阵；

E、建立重建面上的声压值与多虚拟球等效源强产生的声压值之间的传递矩阵；

F、联立步骤D和E获得的传递矩阵，同时对步骤D获得的传递矩阵进行正则化方法处理，获得重建面上的声压求解方程，求解重建面上的声压值。

2.如权利要求1所述的声场重建方法，其特征在于：

所述步骤B中，将虚拟球区分为质心位置的虚拟球和非质心位置的虚拟球，其中对于质心位置的虚拟球的半径为r₀=0.2；非质心位置的各虚拟球的半径为r₀=0.1。

3.如权利要求1所述的声场重建方法，其特征在于：

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

依照式（1）对测量面上的声压值进行补零扩展得到测量面扩展面上的声压值：

p_E(H⁺)=D·p_E(H)   （1）；

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

4.如权利要求1所述的声场重建方法，其特征在于：

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

测量面扩展面声压值与虚拟球等效源声压值之间的传递矩阵如式（2）：

p_E(H⁺)=[T_H]Q   （2）；

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

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

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

$p\left(r\right)=\underset{t=1}{\overset{t}{\mathrm{\Σ}}}{\∫}_{S^{\′}}\mathrm{\σ}\left(r_Q\right)K(r,r_Q){\mathrm{dS}}^{\′}---\left(3\right);$

式中，r_Q是虚源面上的某一点，t为多虚拟球的个数，σ(r_Q)为待求的源强密度函数，K(r,r_Q)为积分核函数，K(r,r_Q)=g(r,r_Q)=(14πR)e^ikR将未知源强密度函数进行双向Fourier级数展开，并利用二维快速FFT计算积分，采用梯形公式离散格林函数，将与距离有关的矩阵规整成一个矩阵，即为[T_H]。

5.如权利要求1所述的声场重建方法，其特征在于：

所述步骤E中重建面声压值与虚拟球等效源声压值之间的传递矩阵如下：

{p_R}=[T_R]Q   （4）；

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

6.如权利要求1所述的声场重建方法，其特征在于：所述步骤F中重建面声压值的求解具体步骤如下：

联立式（2）、（4），求解重建面上的声压：

{p_R}=[T_R][T_H]^-1p_E(H⁺)   （5）；

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

$\left\{p_R\right\}=T_R{(\mathrm{\αI}+T_H^H{T}_H)}^{-1}{T}_H^H\left\{p_{H^{+}}\right\}---\left(6\right);$

其中，α为正则化参数，I为单位对角矩阵，

为T_H的共轭转置矩阵，

为

的逆矩阵。

Images