CN102121847A - 一种瞬态声场重建方法 - Google Patents

Abstract

本发明公开了一种瞬态声场重建方法，其特征是：在被测声场中设置全息面H，在全息面上均匀分布M个测量点；采用M个传声器同步测量全息面H上M个测量点的声压信号，获得各测量点的时域声压；在声源面S内设置虚源面S^*，在虚源面上分布N个时域等效源；采用线性插值函数对各测量点的时域声压进行插值处理，建立任意时刻全息面H上M个测量点的时域声压与M个时域等效源之间的传递关系；从初始时刻开始依次获得对应时刻时域等效源的源强；根据获得的各时刻时域等效源的源强以及时域等效源与任意场点之间的传递关系，重建任意场点的时域声压。本发明方法计算稳定性好、精度高、速度快、实施简单，特别适用于受到瞬态激励力的作用而产生的瞬态声场的重建。

Description

一种瞬态声场重建方法

技术领域

本发明涉及一种适用于噪声源产生的瞬态声场重建的瞬态声场重建方法。

背景技术

在计算三维辐射声场和散射声场时，最常用的方法是边界元法BEM，但BEM法存在着奇异积分处理困难和特征频率处解的非唯一性处理困难的问题。为寻找新的替代方法，最近20年来出现了各种方法，例如波叠加法、基本解法、声源模拟法、等效源法等，这些方法都有一个共同的思想，即通过在声源的内部放置满足给定边界条件的虚源来近似声源外部声场，同理也可以通过在声源的外部放置满足给定边界条件的虚源来近似声源内部声场。统称这些方法为等效源法。但等效源法在计算三维辐射声场和散射声场时主要被用于频域。Kropp和Svensson首次推导了等效源法的时域公式，并将该公式应用于箱体的辐射声场和散射声场的计算中；随后Lee提出了时域移动等效源法，该方法可实现在时域内对散射声场进行计算，他们运用该方法对旋翼飞机的噪声进行了分析。

等效源法不仅可以被用来计算辐射声场和散射声场，而且还可以与近场声全息技术相结合，实现声场的重建。Bi等提出了基于等效源法的近场声全息技术，并将其运用于多源声场和半自由声场；Sarkissian将等效源法运用到了Patch近场声全息中。但他们都是在频域中进行计算，并要求声源信号是稳态信号。Wu等运用Helmholtz最小二乘法在时域内重建了瞬态信号，但该方法要求声源面为球面。Hald等提出的时域全息技术，该技术需要对时域信号做短时傅立叶变换，对每一小段数据在频域进行全息计算，再变换到时域，因此计算时间较长、计算精度较低。Blais等提出采用拉普拉斯变换实现时域全息重建，但是该方法要求声源面仅为平面。

发明内容

本发明所解决的技术问题是避免上述现有技术所存在的不足之处，提供一种采用时域等效源法重建瞬态声场的方法，以全息面上时域声压为输入量，采用时域等效源法实现，直接在时域进行求解的计算稳定性好、计算精度高、计算速度快、实施简单的瞬态声场重建方法。

本发明解决技术问题所采用的技术方案是：

本发明瞬态声场重建方法的特点是：

在被测声场中设置全息面H，在全息面上均匀分布M个测量点；采用M个传声器同步测量全息面H上M个测量点的声压信号，获得各测量点的时域声压；在声源面S内设置虚源面S^*，在虚源面S^*上分布N个时域等效源；

采用线性插值函数对各测量点的时域声压进行插值处理，建立任意时刻全息面H上M个测量点的时域声压与M个时域等效源之间的传递关系；从初始时刻开始依次获得对应时刻时域等效源的源强；

根据获得的各时刻时域等效源的源强以及时域等效源与任意场点之间的传递关系，重建任意场点的时域声压。

本发明瞬态声场重建方法的特点也在于按如下步骤进行：

步骤a、测量全息面H上的时域声压

在由声源面S辐射的被测声场中，布置有全息面H；在全息面H上均匀分布M个测点；采用M个传声器同步测量全息面H上各测点处的时域声压；所述被测声场为瞬态声场、稳态声场或非稳态声场；

步骤b、在声源面S内设置虚源面S^*，在所述虚源面S^*上分布N个时域等效源；所述时域等效源的个数不大于对应全息面H上分布的测点数；所述时域等效源为标准点源、面源或体源；

步骤c、采用线性插值函数对各测点的时域声压进行插值处理，建立任意时刻全息面H上M个测点的时域声压与时域等效源之间的传递关系：

$P_H^i={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^1+{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^2+{\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^j+\·\·\·+{\mathrm{\Φ}}_H^i\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^i,$ 其中

$P_H^i={\left[\begin{array}{cccc}{p}_1\left(t_i\right)& p_2\left(t_i\right)& \·\·\·& p_M\left(t_i\right)\end{array}\right]}^T,$

${\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)=\left[\begin{array}{cccc}\frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{12}^i\right)}{R_{12}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{12}^i\right)}{R_{12}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{1N}^i\right)}{R_{1N}}\\ \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{22}^i\right)}{R_{22}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{22}^i\right)}{R_{22}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{2N}^i\right)}{R_{2N}}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{M2}^i\right)}{R_{M2}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{M2}^i\right)}{R_{M2}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{\mathrm{MN}}^i\right)}{R_{\mathrm{MN}}}\end{array}\right],$

$Q^j={\left[\begin{array}{cccc}{q}_1^j& q_2^j& \·\·\·& q_N^j\end{array}\right]}^T,$

${\mathrm{\τ}}_{\mathrm{mn}}^i=t_i-R_{\mathrm{mn}}/c;$ τ^t＝t_i-R_min/c；j＝1，2，…i；

上标T表示矩阵的转置；下标H表示全息面；

上标i表示第i个时间步阶，i＝1，2，…l，I为总采样点数；

Q^j为t_j时刻等效源源强列向量；

p_m(t_i)为t_i时刻第m个测点处的时域声压；

为线性插值函数；

R_mn为第m个测点与第n个等效源点之间的距离；

R_min为测点与等效源点之间的距离的最小值；

为虚源面上t_j时刻第n个等效源点处时域等效源的源强；

c为声音传播速度；

步骤d、从t₁时刻开始依次获得对应时刻时域等效源的源强

当i＝1时， $P_H^1={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^1\right)Q^1;$

由于

和

均为已知量，因此可得：

其中：为的广义逆矩阵；

当i＝2时， $P_H^2={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^1+{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^2;$

其中

和

均为已知量，Q¹已求得，因此可得

$Q^2={\left({\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)\right)}^{+}(P_H^2-{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^1);$

根据依次获得的等效源源强，可获得t_i时刻等效源源强列向量为：

$Q^j={\left[{\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)\right]}^{+}[P_H^j-{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)Q^1-{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)Q^2-\·\·\·-{\mathrm{\Φ}}_H^{j-1}\left({\mathrm{\τ}}_{\mathrm{mn}}^{\mathrm{ij}}\right)Q^{j-1}];$

步骤e、根据获得的各时刻时域等效源的源强以及时域等效源与任意场点之间的传递关系，重建任意场点的时域声压；

根据步骤d获得各个时刻的等效源源强，得任意t_i时刻任意场点f处时域声压为：

$p_f\left(t_i\right)={\mathrm{\Φ}}_f^1\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^1+{\mathrm{\Φ}}_f^2\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^2+\·\·\·+{\mathrm{\Φ}}_f^j\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^j+\·\·\·+{\mathrm{\Φ}}_f^1\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^1,$ 其中：

${\mathrm{\τ}}_{\mathrm{fn}}^i=t_i-R_{\mathrm{fn}}/c.$

所述各测点上的时域声压的测量是采用传声器同步测量获得。

所述全息面H和声源面S为封闭面或非封闭面，为平面或任意曲面。

理论模型：

1、时域等效源法的公式，参见图1，假设空间中存在一个任意形状的声源，声源面为S。若在声源的内部布置N个虚拟的等效源，在此等效源为标准点源，则声源在声场中任意一点任意时刻所产生的时域声压可由N个等效源的源强与脉冲响应函数的卷积叠加得到，即：

$p_m\left(t\right)={\mathrm{\Σ}}_{n=1}^N[q_n\left(t\right)*g_{\mathrm{mn}}(R_{\mathrm{mn}},t)]---\left(1\right)$

式中，符号“*”表示两个函数之间的卷积；m表示声场中第m个测点；n表示第n个时域等效源；q_n(t)表示t时刻第n个时域等效源的源强；R_mn表示点m与第n个虚拟的等效源之间的距离；g_mn(R_mn，t)表示脉冲响应函数，其表达式为：

$g_{\mathrm{mn}}(R_{\mathrm{mn}},t)=\frac{1}{R_{\mathrm{mn}}}\mathrm{\δ}(t-R_{\mathrm{mn}}/c)---\left(2\right)$

式中，δ(x)表示Dirac函数。将公式(2)代入公式(1)可得

$p_m\left(t\right)={\mathrm{\Σ}}_{n=1}^N\frac{1}{R_{\mathrm{mn}}}[q_n\left(t\right)*\mathrm{\δ}(t-R_{\mathrm{mn}}/c)]={\mathrm{\Σ}}_{n=1}^N\frac{1}{R_{\mathrm{mn}}}{q}_n(t-R_{\mathrm{mn}}/c)---\left(3\right)$

若令τ_mn＝t-R_mn/c，则公式(3)可写成

$p_m\left(t\right)={\mathrm{\Σ}}_{n=1}^N\frac{1}{R_{\mathrm{mn}}}{q}_n\left({\mathrm{\τ}}_{\mathrm{mn}}\right)---\left(4\right)$

要想在计算机上实现对公式(4)的计算，则必须对时间t和时间τ分别进行离散化处理。设

t_i＝t₁+(i-1)Δt      (5)

τ^j＝τ¹+(j-1)Δτ   (6)

式中，t₁和τ¹均表示初始时刻，即当t≤t₁时声源不发声，当τ≤τ¹时，等效源强度为0；i＝1，2，…I，j＝1，2，…；Δt和Δτ均表示时间步长。由于R_mn会随着等效源和场点位置的变化而变化，所以无法保证t_i-R_mn/c全部落在所设定的τ的时间步阶上，这就为计算等效源的强度带来了困难。

2、插值处理，为解决上述问题，对等效源强度的时间函数采取了插值处理。令：

$q_n\left(\mathrm{\τ}\right)={\mathrm{\Σ}}_{j=1}^J{\mathrm{\Φ}}^j\left(\mathrm{\τ}\right)q_n^j---\left(7\right)$

式中，Φ^j(τ)为线性插值函数，其表达式为

将公式(7)代入公式(4)中可得：

$p_m\left(t\right)={\mathrm{\Σ}}_{n=1}^N{\mathrm{\Σ}}_{j=1}^J\frac{1}{R_{\mathrm{mn}}}{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{\mathrm{mn}}\right)q_n^j---\left(9\right)$

由公式(9)可知，通过插值保证了用于计算的等效源的强度所对应的时刻是所设定的τ的时间步阶。

3、声场重建，设Δt＝Δτ，τ¹＝t₁-R_min/c，I＝j，其中R_min表示取所有R_mn中的最小值。则当t取某一时间步阶t₁时，

因此：

$p_m\left(t_i\right)={\mathrm{\Σ}}_{n=1}^N{\mathrm{\Σ}}_{j=1}^i\frac{1}{R_{\mathrm{mn}}}{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)q_n^j={\mathrm{\Σ}}_{j=1}^i{\mathrm{\Σ}}_{n=1}^N\frac{1}{R_{\mathrm{mn}}}{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)q_n^j---\left(10\right)$

将公式(10)写成向量形式为：

$p_m\left(t_i\right)=\left[\begin{array}{cccc}\frac{{\mathrm{\Φ}}^1\left({\mathrm{\τ}}_{m1}^i\right)}{R_{m1}}& \frac{{\mathrm{\Φ}}^2\left({\mathrm{\τ}}_{m2}^i\right)}{R_{m2}}& \·\·\·& \frac{{\mathrm{\Φ}}^i\left({\mathrm{\τ}}_{\mathrm{mN}}^i\right)}{R_{\mathrm{mN}}}\end{array}\right]\left[\begin{array}{c}{q}_1^1\\ q_2^1\\ \·\\ \·\\ \·\\ q_N^1\end{array}\right]+$ (11)

$\·\·\·+\left[\begin{array}{cccc}\frac{{\mathrm{\Φ}}^i\left({\mathrm{\τ}}_{m1}^i\right)}{R_{m1}}& \frac{{\mathrm{\Φ}}^i\left({\mathrm{\τ}}_{m2}^i\right)}{R_{m2}}& \·\·\·& \frac{{\mathrm{\Φ}}^i\left({\mathrm{\τ}}_{\mathrm{mN}}^i\right)}{R_{\mathrm{mN}}}\end{array}\right]\left[\begin{array}{c}{q}_1^i\\ q_2^i\\ \·\\ \·\\ \·\\ q_N^i\end{array}\right]$

假设声场中测点的个数为M个，将这M个测点在t_i时刻的时域声压写成向量形式，即：

$P_H^i={\left[\begin{array}{cccc}{p}_1\left(t_i\right)& p_2\left(t_i\right)& \·\·\·& p_M\left(t_i\right)\end{array}\right]}^T---\left(12\right)$

式中，下标H表示全息面，上标i表示第i个时间步阶。由公式(11)可推得：

$P_H^i={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^1+{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^2+{\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^j+\·\·\·+{\mathrm{\Φ}}_H^i\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^i---\left(13\right)$

式中，j＝1，2，…i；

为t_i时刻全息面时域声压与等效源之间的传递矩阵，其形式为：

${\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)=\left[\begin{array}{cccc}\frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{12}^i\right)}{R_{12}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{12}^i\right)}{R_{12}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{1N}^i\right)}{R_{1N}}\\ \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{22}^i\right)}{R_{22}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{22}^i\right)}{R_{22}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{2N}^i\right)}{R_{2N}}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{M2}^i\right)}{R_{M2}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{M2}^i\right)}{R_{M2}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{\mathrm{MN}}^i\right)}{R_{\mathrm{MN}}}\end{array}\right]---\left(14\right)$

Q^j为t_j时刻等效源源强列向量，其形式为：

$Q^j={\left[\begin{array}{cccc}{q}_1^j& q_2^j& \·\·\·& q_N^j\end{array}\right]}^T---\left(15\right)$

由公式(13)可知，当i＝1时，

$P_H^1={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^1\right)Q^1---\left(16\right)$

在公式(16)中，和

均为已知量，因此可求得：

$Q^1={\left[{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^1\right)\right]}^{+}{P}_H^1---\left(17\right)$

式中，

为的广义逆矩阵。当i＝2时，

$P_H^2={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^1+{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^2---\left(18\right)$

在公式(18)中，

和

均为已知量，Q¹已由公式(17)求得，因此可求得：

$Q^2={\left({\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)\right)}^{+}(P_H^2-{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^1)---\left(19\right)$

同理，依次类推可获得：

$Q^j={\left[{\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)\right]}^{+}[P_H^j-{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)Q^1-{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)Q^2-\·\·\·-{\mathrm{\Φ}}_H^{j-1}\left({\mathrm{\τ}}_{\mathrm{mn}}^{\mathrm{ij}}\right)Q^{j-1}]---\left(20\right)$

一旦获得各个时刻的等效源源强，则可求得t_i时刻任意场点f处的时域声压为：

$p_f\left(t_i\right)={\mathrm{\Φ}}_f^1\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^1+{\mathrm{\Φ}}_f^2\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^2+\·\·\·+{\mathrm{\Φ}}_f^j\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^j+\·\·\·+{\mathrm{\Φ}}_f^1\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^1---\left(21\right)$

式中 ${\mathrm{\τ}}_{\mathrm{fn}}^i=t_i-R_{\mathrm{fn}}/c.$

在上述求得等效源源强的过程中为保证解的唯一性，须满足M≥N，同时求得等效源源强的过程也属于逆问题，因此可采用Tikhonov正则化方法稳定求解过程，正则化参数用GCV法来进行选取。

通过上述方法，实现了任意场点任意时刻瞬态声场的重建。

与已有技术相比，本发明的有益效果：

1、本发明采用时域等效源法作为重建算法，与传统的方法相比，本发明方法具有计算稳定性好、计算精度高、计算速度快、实施简单等优点。

2、本发明采用采用时域等效源法建立声场与声源之间的关系，直接在时域进行计算，不需要进行任何频域变换。

3、本发明所采用的全息面和声源面，可以是任意形状面，解决了传统方法只能适用于平面、柱面或球面等规则形状的缺陷。

4、本发明方法可以用于瞬态声场、稳态声场、非稳态声场等分析。

附图说明

图1为声源、时域等效源以及全息面位置示意图；

图2爆炸性膨胀球在点A的时域声压理论值和重建值的比较；

图3爆炸性膨胀球在点B的时域声压理论值和重建值的比较；

图4爆炸性膨胀球在t＝0.04ms时刻重建面上所有点的时域声压理论值和重建值的比较；

图5爆炸性膨胀球在t＝0.14ms时刻重建面上所有点的时域声压理论值和重建值的比较；

图6脉冲加速度刚性球在点A的时域声压理论值和重建值的比较；

图7脉冲加速度刚性球在点B的时域声压理论值和重建值的比较；

图8脉冲加速度刚性球在t＝0.04ms时刻重建面上所有点的时域声压理论值和重建值的比较；

图9脉冲加速度刚性球在t＝0.14ms时刻重建面上所有点的时域声压理论值和重建值的比较；

具体实施方式

参见图1，本实施例中，在由声源面S构成的被测声场中，布置有全息面H；在全息面H上均匀分布M个测点；采用M个传声器同步测量全息面H上各测点处的时域声压；在声源面S内部有虚源面S^*，在虚源面S^*上布置有N个时域等效源；全息面H上测点数M不小于虚源面S^*上时域等效源数N。

瞬态声场重建是按如下步骤进行：

步骤a、测量全息面H上的时域声压

在由声源面S辐射的被测声场中，布置有全息面H；在全息面H上均匀分布M个测点；采用M个传声器同步测量全息面H上各测点处的时域声压；被测声场为瞬态声场、稳态声场或非稳态声场；全息面H和声源面S为封闭面或非封闭面，为平面或任意曲面。

步骤b、在声源面S内设置虚源面S^*，在虚源面S^*上分布N个时域等效源；时域等效源的个数不大于对应全息面H上分布的测点数；时域等效源为标准点源、面源或体源；

步骤c、采用线性插值函数对各测点的时域声压进行插值处理，建立任意时刻全息面H上M个测点的时域声压与时域等效源之间的传递关系：

$P_H^i={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^1+{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^2+{\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^j+\·\·\·+{\mathrm{\Φ}}_H^i\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^i,$ 其中

$P_H^i={\left[\begin{array}{cccc}{p}_1\left(t_i\right)& p_2\left(t_i\right)& \·\·\·& p_M\left(t_i\right)\end{array}\right]}^T,$

${\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)=\left[\begin{array}{cccc}\frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{12}^i\right)}{R_{12}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{12}^i\right)}{R_{12}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{1N}^i\right)}{R_{1N}}\\ \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{22}^i\right)}{R_{22}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{22}^i\right)}{R_{22}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{2N}^i\right)}{R_{2N}}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{M2}^i\right)}{R_{M2}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{M2}^i\right)}{R_{M2}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{\mathrm{MN}}^i\right)}{R_{\mathrm{MN}}}\end{array}\right],$

$Q^j={\left[\begin{array}{cccc}{q}_1^j& q_2^j& \·\·\·& q_N^j\end{array}\right]}^T,$

${\mathrm{\τ}}_{\mathrm{mn}}^i=t_i-R_{\mathrm{mn}}/c;$ τⁱ＝t_i-R_min/c；j＝1，2，…i；

上标T表示矩阵的转置；下标H表示全息面；

上标i表示第i个时间步阶，i＝1，2，…I，I为总采样点数；

Q^j为t_j时刻等效源源强列向量；

p_m(t_i)为t_i时刻第m个测点处的时域声压；

为线性插值函数；

R_mn为第m个测点与第n个等效源点之间的距离；

R_min为测点与等效源点之间的距离的最小值；

为虚源面上t_j时刻第n个等效源点处时域等效源的源强；

c为声音传播速度；

步骤d、从t₁时刻开始依次获得对应时刻时域等效源的源强

当i＝1时， $P_H^1={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^1\right)Q^1;$

由于

和

均为已知量，因此可得：

其中：

为

的广义逆矩阵；

当i＝2时， $P_H^2={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^1+{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^2;$

其中

和均为已知量，Q¹已求得，因此可得

$Q^2={\left({\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)\right)}^{+}(P_H^2-{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^1);$

根据依次获得的等效源源强，可获得t_i时刻等效源源强列向量为：

$Q^j={\left[{\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)\right]}^{+}[P_H^j-{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)Q^1-{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)Q^2-\·\·\·-{\mathrm{\Φ}}_H^{j-1}\left({\mathrm{\τ}}_{\mathrm{mn}}^{\mathrm{ij}}\right)Q^{j-1}];$

步骤e、根据获得的各时刻时域等效源的源强以及时域等效源与任意场点之间的传递关系，重建任意场点的时域声压；

根据步骤d获得各个时刻的等效源源强，得任意t_i时刻任意场点f处时域声压为：

$p_f\left(t_i\right)={\mathrm{\Φ}}_f^1\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^1+{\mathrm{\Φ}}_f^2\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^2+\·\·\·+{\mathrm{\Φ}}_f^j\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^j+\·\·\·+{\mathrm{\Φ}}_f^1\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^1,$ 其中：

${\mathrm{\τ}}_{\mathrm{fn}}^i=t_i-R_{\mathrm{fn}}/c.$

方法的检验：

仿真1

声源为一个处于自由空间中的爆炸性膨胀球，以球心为球面坐标系的原点，则该球的表面振速可表示为：

v(r₀，θ，φ，t)＝v_s·E(t-t₁)     (22)

式中，r₀表示球的半径，在此取r₀＝00138m；v_s表表示振速幅值，在此取v_s＝1m/s，；E表示Heaviside函数；t₁初始的时间步阶，在此取t₁＝0s。则该声源所产生的时域声压的解析解表达式为：

$p(r,\mathrm{\θ},\mathrm{\φ},t)={\mathrm{\ρ}}_0{\mathrm{cv}}_s{r}_0\frac{e^{-(\mathrm{ct}-r+r_0)/r_0}}{r}E(t-\frac{r-r_0}{c})---\left(23\right)$

式中，r表示声场中任意场点到坐标系原点的距离，θ为r与z轴正方向之间的夹角，φ为r与x轴正方向之间沿逆时针方向的夹角；ρ₀表示声场中的介质密度。

设全息面是半径r_h＝00344m的球面，将该球面沿θ方向分成7等份，沿φ方向分成14等份，共有86个结点。重建面是半径r_s＝0.0275m的球面，其球面网格划分与全息面相同。虚源面是半径r_v＝0.0034m的球面，将该球面沿θ方向分成4等份，沿φ方向分成8等份，共有26个结点，将等效源布置在这些结点上。设时间采样频率为10kHz，时间采样点数为101。为了定量地表示重建结果和理论结果之间的吻合程度，定义声压理论值和声压重建值之间的均方根误差计算公式为：

$\mathrm{\μ}=\frac{\left|\right|P_r-P_t\left|\right|}{\left|\right|P_t\left|\right|}\×100\%---\left(24\right)$

式中，P_r表示声压重建值，P_t表示声压理论值，||P_t||表示矩阵P_t的二范数。

为检验重建面上各点的时域波形的重建效果，在重建面选取两点A和B，点A和点B的坐标分另为A(r，θ，φ)＝A(0.0275m，π/7，π/7)，B(r，θ，φ)＝B(0.0275m，3π/7，-6π/7)。点A和点B的时域声压理论值和重建值的比较参见图2和图3。利用公式(24)分别计算两点的均方根误差，点A的均方根误差为0.04％，点B的均方根误差为0.05％。

为检验任意时刻整个重建面的重建效果，选取两个时刻t＝0.04ms和0.14ms，在这两个时刻重建面上所有点的时域声压理论值和重建值的比较参见图4和图5。利用公式(24)分别计算了两个时刻的均方根误差，t＝0.04ms时刻的均方根误差为0.01％，t＝0.14ms时刻的均方根误差为0.14％。

上述仿真结果表明利用本发明方法所获得的时域声压重建值与理论值很好地吻合。

仿真2

声源为一个处于自由空间中的脉冲加速度刚性球，以球心为球面坐标系的原点，则该球的表面振速可表示为：

v(r₀，θ，φ，t)＝v_scosθ·E(t-t₁)    (25)

式中，0≤θ≤π，其它参数取值同仿真1。则该声源所产生的时域声压的解析解表达式为：

$p(r,\mathrm{\θ},\mathrm{\φ},t)={\mathrm{\ρ}}_0{\mathrm{cv}}_s{r}_0\frac{e^{-(\mathrm{ct}-r+r_0)/r_0}}{r}E(t-\frac{r-r_0}{c})$ (26)

$\×[\cos \left(\frac{\mathrm{ct}-r+r_0}{r_0}\right)-(1-\frac{r_0}{r})\sin \left(\frac{\mathrm{ct}-r+r_0}{r_0}\right)]$

全息面位置和网格划分、重建面位置和网格划分、等效源布置位置、时间采样频率和时间采样点数均同仿真1。

为检验重建面上各点的时域波形的重建效果，在重建面选取两点A和B，点A和点B的坐标分另为A(r，θ，φ)＝A(0.0275m，π/7，π/7)，B(r，θ，φ)＝＝B(0.0275m，3π/7，-6π/7)。点A和点B的时域声压理论值和重建值的比较参见图6和图7。利用公式(24)分别计算了两点的均方根误差，点A的均方根误差为3.98％，点B的均方根误差为4.44％

为检验任意时刻上整个重建面的重建效果，选取两个时刻t＝0.04ms和0.14ms，在这两个时刻重建面上所有点的时域声压理论值和重建值的比较参见图8和图9。利用公式(24)分别计算了两个时刻的均方根误差，t＝0.04ms时刻的均方根误差为2.26％，t＝0.14ms时刻的均方根误差为10.92％。

该仿真结果进一步表明利用本发明方法所获得的时域声压重建值与理论值具有很好地吻合效果。

Claims (4)

1.一种瞬态声场重建方法，其特征是：

在被测声场中设置全息面H，在全息面上均匀分布M个测量点；采用M个传声器同步测量全息面H上M个测量点的声压信号，获得各测量点的时域声压；在声源面S内设置虚源面S^*，在虚源面S^*上分布N个时域等效源；

采用线性插值函数对各测量点的时域声压进行插值处理，建立任意时刻全息面H上M个测量点的时域声压与M个时域等效源之间的传递关系；从初始时刻开始依次获得对应时刻时域等效源的源强；

根据获得的各时刻时域等效源的源强以及时域等效源与任意场点之间的传递关系，重建任意场点的时域声压。

2.根据权利要求1所述的瞬态声场重建方法，其特征是按如下步骤进行：

步骤a、测量全息面H上的时域声压

在由声源面S辐射的被测声场中，布置有全息面H；在全息面H上均匀分布M个测点；采用M个传声器同步测量全息面H上各测点处的时域声压；所述被测声场为瞬态声场、稳态声场或非稳态声场；

步骤b、在声源面S内设置虚源面S^*，在所述虚源面S^*上分布N个时域等效源；所述时域等效源的个数不大于对应全息面H上分布的测点数；所述时域等效源为标准点源、面源或体源；

步骤c、采用线性插值函数对各测点的时域声压进行插值处理，建立任意时刻全息面H上M个测点的时域声压与时域等效源之间的传递关系：

$P_H^i={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^1+{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^2+{\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^j+\·\·\·+{\mathrm{\Φ}}_H^i\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)Q^i,$ 其中

$P_H^i={\left[\begin{array}{cccc}{p}_1\left(t_i\right)& p_2\left(t_i\right)& \·\·\·& p_M\left(t_i\right)\end{array}\right]}^T,$

${\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^i\right)=\left[\begin{array}{cccc}\frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{12}^i\right)}{R_{12}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{12}^i\right)}{R_{12}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{1N}^i\right)}{R_{1N}}\\ \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{22}^i\right)}{R_{22}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{22}^i\right)}{R_{22}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{2N}^i\right)}{R_{2N}}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{M2}^i\right)}{R_{M2}}& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{M2}^i\right)}{R_{M2}}& \·\·\·& \frac{{\mathrm{\Φ}}^j\left({\mathrm{\τ}}_{\mathrm{MN}}^i\right)}{R_{\mathrm{MN}}}\end{array}\right],$

$Q^j={\left[\begin{array}{cccc}{q}_1^j& q_2^j& \·\·\·& q_N^j\end{array}\right]}^T,$

${\mathrm{\τ}}_{\mathrm{mn}}^i=t_i-R_{\mathrm{mn}}/c;$ τ^t＝t_i-R_min/c；j＝1，2，…i；

上标T表示矩阵的转置；下标H表示全息面；

上标i表示第i个时间步阶，i＝1，2，…I，I为总采样点数；

Q^j为t_j时刻等效源源强列向量；

p_m(t_i)为t_i时刻第m个测点处的时域声压；

为线性插值函数；

R_mn为第m个测点与第n个等效源点之间的距离；

R_min为测点与等效源点之间的距离的最小值；

为虚源面上t_j时刻第n个等效源点处时域等效源的源强；

c为声音传播速度；

步骤d、从t₁时刻开始依次获得对应时刻时域等效源的源强

当i＝1时， $P_H^1={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^1\right)Q^1;$

由于

和

均为已知量，因此可得：其中：

为

的广义逆矩阵；

当i＝2时， $P_H^2={\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^1+{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^2;$

其中和

均为已知量，Q¹已求得，因此可得

$Q^2={\left({\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)\right)}^{+}(P_H^2-{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^2\right)Q^1);$

根据依次获得的等效源源强，可获得t_i时刻等效源源强列向量为：

$Q^j={\left[{\mathrm{\Φ}}_H^j\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)\right]}^{+}[P_H^j-{\mathrm{\Φ}}_H^1\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)Q^1-{\mathrm{\Φ}}_H^2\left({\mathrm{\τ}}_{\mathrm{mn}}^j\right)Q^2-\·\·\·-{\mathrm{\Φ}}_H^{j-1}\left({\mathrm{\τ}}_{\mathrm{mn}}^{\mathrm{ij}}\right)Q^{j-1}];$

步骤e、根据获得的各时刻时域等效源的源强以及时域等效源与任意场点之间的传递关系，重建任意场点的时域声压；

根据步骤d获得各个时刻的等效源源强，得任意t_i时刻任意场点f处时域声压为：

$p_f\left(t_i\right)={\mathrm{\Φ}}_f^1\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^1+{\mathrm{\Φ}}_f^2\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^2+\·\·\·+{\mathrm{\Φ}}_f^j\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^j+\·\·\·+{\mathrm{\Φ}}_f^1\left({\mathrm{\τ}}_{\mathrm{fn}}^i\right)Q^1,$ 其中：

${\mathrm{\τ}}_{\mathrm{fn}}^i=t_i-R_{\mathrm{fn}}/c.$

3.根据权利要求1所述的瞬态声场重建方法，其特征是所述各测点上的时域声压的测量是采用传声器同步测量获得。

4.根据权利要求1所述的瞬态声场重建方法，其特征是所述全息面H和声源面S为封闭面或非封闭面，为平面或任意曲面。

Images