CN106124044A - 实心球声源识别低旁瓣超高分辨率声学图像快速获取方法 - Google Patents

Abstract

本发明公开一种实心球声源识别低旁瓣超高分辨率声学图像快速获取方法。在识别三维空间内不相干声源时，本发明为快速获得清晰干净的超高分辨率声源成像图提供了有效途径，同时相对于实心球函数型延迟求和方法能更准确地量化声源，对声源识别结果的准确分析具有重要意义。

Description

实心球声源识别低旁瓣超高分辨率声学图像快速获取方法

技术领域

本发明涉及三维空间声源识别领域

背景技术

凭借旋转对称性好、声场记录全面、衍射作用强、记录信号信噪比高等优势，实心球形传声器阵列广泛流行于三维声源识别领域。球谐函数波束形成(Spherical HarmonicsBeamforming,SHB)是其普遍采用的声信号后处理方法。但是受球谐函数阶次截断、传声器离散采样等因素影响，该方法输出的声源成像图被大量高水平旁瓣污染，同时空间分辨率低，这些缺陷会使声源识别结果的分析具有严重不确定性。目前三维空间已有的反卷积(DAMAS)方法在SHB方法输出结果、阵列点传播函数、声源分布之间建立线性方程组，通过在反复迭代过程中引入正约束来定解该方程组，从而提取真实声源分布，可获得超高空间分辨率，但存在严重耗时问题；实心球的函数型延迟求和方法(FDAS)在实心球延迟求和(DAS)方法的基础上发展而来，能快速衰减DAS方法产生的旁瓣并能准确探测弱源，但该方法仅能轻微改善DAS方法空间分辨率，无法准确分辨彼此靠近的声源，并且实际应用中对声源的量化偏差也较大。

发明内容

本发明的目的是针对实心球阵列三维声源识别低旁瓣、超高分辨率声学图像的快速获取问题提供一种方法，同时克服FDAS方法对声源量化偏差大的不足。

为实现本发明目的而采用的技术方案如下，实心球声源识别低旁瓣超高分辨率声学图像快速获取方法：

1)：构建测量***：

参见图1，一个处于三维空间内的实心球阵列；所述实心球阵列包括一个半径为a的实心球，以及分布在所述实心球表面的Q个传声器；这些传声器的编号为q，q＝1,2,…,Q；所述实心球的球心作为原点，实心球阵列所处的三维空间内任意一点的位置坐标用(r,Ω)描述，r表示所描述位置与原点间的距离，表示所描述位置的方向，θ、分别为其仰角、方位角；(a,Ω_q)为q号传声器的位置坐标；

2)获取声压

各个传声器接收声压信号p(ka,Ω_q)；其中，处于所述三维空间内的单极子点声源，辐射声波为球面波；波数k＝2πf/c，声波频率为f，声速为c；

3)通过各个传声器接收声压信号P，得到互谱矩阵C：

$C=\stackrel{\‾}{{\mathrm{pp}}^H},$

P为各个传声器接收的复声压信号，

$p={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0),$

B表示所有声源位置坐标组成的集合，(r₀,Ω₀)和(r₀',Ω₀')均表示声源的位置坐标，s(kr₀,Ω₀)表示声源的强度，t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^T，t_q(kr₀,Ω₀)表示声源到q号传声器的声场传递函数，

$t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≡\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_0,ka)Y_n^{m*}\left({\mathrm{\Ω}}_0\right)Y_n^m\left({\mathrm{\Ω}}_q\right),$

分别为Ω₀、Ω_q方向的球谐函数，n、m为球谐函数的阶次，R_n(kr₀,ka)为声信号传播径向函数，

$R_n({\mathrm{kr}}_0,ka)=-4{\mathrm{\πih}}_n^{\left(2\right)}\left({\mathrm{kr}}_0\right)\left(j_n\right(ka)-\frac{j_n^{\′}\left(ka\right)}{h_n^{\left(2\right)\′}\left(ka\right)}{h}_n^{\left(2\right)}(ka\left)\right),$

其中，被减数项为球面波在自由场中传播的径向函数，减数项为实心球引起的球面波散射径向函数，为虚数单位，j_n(ka)为n阶第一类球贝塞尔函数，为n阶第二类球汉克尔函数，j'_n(ka)、分别为j_n(ka)、的一阶导数。

上标“-”、“H”、“*”和“T”分别表示平均运算、转置共轭运算、共轭运算和转置运算；

4)确定各个频率对应的球谐函数阶次截断长度：

表示将数值向正无穷方向圆整到最近的整数；

5)将互谱矩阵C和球谐函数阶次截断长度N代入FDAS方法输出式：

$W_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\frac{1}{Q}{\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)C^{\frac{1}{\mathrm{\ξ}}}v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2}\right)}^{\mathrm{\ξ}}---\left(I\right)$

其中：(r_f,Ω_f)为聚焦点的位置坐标，v(kr_f,Ω_f)＝[v₁(kr_f,Ω_f),v₂(kr_f,Ω_f),…,v_Q(kr_f,Ω_f)]^T为聚焦列向量，v_q(kr_f,Ω_f)表示q号传声器的聚焦分量，R_n(kr_f,ka)为聚焦径向函数。|| ||₂表示2范数，指数参数ξ推荐取值为16。

6)对(Ⅰ)构造脊探测(RD)对象，

L(Ω_f)＝|W_F(kr_f,Ω_f)|

“||”表示取模运算，对每个聚焦方向，可按下列步骤检测其是否为极大值脊：

步骤1.构造L(Ω_f)的二阶偏导数矩阵，即黑塞矩阵：

$H=\left[\begin{array}{cc}\frac{{\∂}^{2}L}{\∂{\mathrm{\φ}}_f^2}& \frac{{\∂}^{2}L}{\∂{\mathrm{\φ}}_f\∂{\mathrm{\θ}}_f}\\ \frac{{\∂}^{2}L}{\∂{\mathrm{\θ}}_f\∂{\mathrm{\φ}}_f}& \frac{{\∂}^{2}L}{\∂{\mathrm{\θ}}_f^2}\end{array}\right]$

步骤2.特征值分解矩阵Η，找出绝对值最大的特征值λ及相应的特征向量u_λ，其中u_λ指示最大表面曲率方向，输出值在u_λ方向上若为极大值则λ为负，若为极小值则λ为正。

步骤3.计算标量ρ(Ω_f)为：

$\mathrm{\ρ}\left({\mathrm{\Ω}}_f\right)=-\frac{1}{2}sign\left(\mathrm{\λ}\right)|sign(\▿L({\mathrm{\Ω}}_f+{\mathrm{\χu}}_{\mathrm{\λ}})\·u_{\mathrm{\λ}})-sign(\▿L({\mathrm{\Ω}}_f-{\mathrm{\χu}}_{\mathrm{\λ}})\·u_{\mathrm{\λ}})|$

其中，“sign”为符号函数，为梯度算子，“．”表示求向量内积，χ为常数，用以控制检测空间精度，本发明取为0.2。当ρ(Ω_f)＝1时，聚焦方向为极大值脊，保留FDAS的输出值；反之，用零取代FDAS在该聚焦方向的输出值来剔除非脊点。

7)对DAS方法结果进行反卷积(DAMAS)求取各声源的平均声压贡献：

定义声源平均声压贡献P_AC(kr₀,Ω₀)为声源在阵列各传声器处产生声压信号自谱的均值，即为互谱矩阵C的迹，并带入DAS输出量表达式，

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈B}{P}_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}$

定义DAS的点传播函数为：

$psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))=\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}$

其表示(r₀,Ω₀)位置处的单位平均声压贡献点声源在聚焦点(r_f,Ω_f)处的波束形成贡献量。DAS输出量为各声源平均声压贡献与对应点传播函数的乘积的和。其中，

t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^T

表示声场传递函数列向量。6)中RD方法获取的脊中包含声源点，且声源点处的输出值相对很大，可默认一定动态范围外的脊和非脊处无声源，输出值为零，仅对动态范围内的脊进行DAMAS运算来获得超高空间分辨率，此时G为参与运算的脊的数目，该数值很小，可使反卷积具有极高的效率。W为G×1维列向量，对应元素为参与运算各脊处的DAS输出量。A为G×G维的点传播函数矩阵，P_AC为G×1维列向量，对应元素为参与运算各脊处声源的平均声压贡献，为未知非负数。如下线性方程组成立：

W＝AP_AC

该线性方程组采用高斯-赛德尔迭代方案求解P_AC，移除点传播函数的影响，能够显著衰减旁瓣并获得超高空间分辨率。初始化基于第l次迭代计算的结果进行第l+1次迭代的具体步骤为：

步骤1.计算残差r_e：

$\begin{array}{c}{r}_e^{\left(l\right)}({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈D}psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))P_{AC}^{(l+1)}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\\ +{\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈E}psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))P_{AC}^{\left(l\right)}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)-W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\end{array}$

其中，D为已完成l+1次迭代计算的声源点组成的集合，E为未完成l+1次迭代计算的脊点组成的集合，二者的并集包含所有参与运算脊点声源。

步骤2.确定

$\begin{array}{cc}{P}_{AC}^{(l+1)}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)=max\left(P_{AC}^{\left(l\right)}\right({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)-\frac{r_e^{\left(l\right)}({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))},0)& r_f=r_0,{\mathrm{\Ω}}_f={\mathrm{\Ω}}_0\end{array}$

按照该方案依次计算每个运算脊点，即可得出

8)根据迭代所得P_AC即能定位声源位置并且能绘制声源成像图实现三维声场可视化。

本发明的方法推导过程如下：

步骤1：获取各传声器接收复声压信号的互谱矩阵

图1为球阵列波束形成坐标系，原点位于阵列中心，三维空间内任意位置可用(r,Ω)描述，r表示所描述位置与原点间的距离，表示所描述位置的方向，θ、分别为其仰角、方位角。符号“●”代表嵌在实心球表面的传声器，记a为球半径，Q为传声器总数，q＝1,2,…,Q为传声器索引，(a,Ω_q)为q号传声器的位置坐标。假设声源为单极子点声源，辐射声波为球面波。声波频率为f，声速为c，波数k＝2πf/c，(r₀,Ω₀)表示声源的位置坐标，t_q(kr₀,Ω₀)表示声源到q号传声器的声场传递函数，根据散射理论，

$t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≡\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_0,ka)Y_n^{m*}\left({\mathrm{\Ω}}_0\right)Y_n^m\left({\mathrm{\Ω}}_q\right)---\left(1\right)$

其中，分别为Ω₀、Ω_q方向的球谐函数，n、m为球谐函数的阶次，上标“*”表示共轭运算，R_n(kr₀,ka)为声信号传播径向函数：

$R_n({\mathrm{kr}}_0,ka)=-4{\mathrm{\πih}}_n^{\left(2\right)}\left({\mathrm{kr}}_0\right)\left(j_n\right(ka)-\frac{j_n^{\′}\left(ka\right)}{h_n^{\left(2\right)\′}\left(ka\right)}{h}_n^{\left(2\right)}(ka\left)\right)---\left(2\right)$

其中，被减数项为球面波在自由场中传播的径向函数，减数项为实心球引起的球面波散射径向函数，为虚数单位，j_n(ka)为n阶第一类球贝塞尔函数，为n阶第二类球汉克尔函数，j'_n(ka)、分别为j_n(ka)、的一阶导数。

用s(kr₀,Ω₀)表示声源的强度，B表示所有声源位置坐标组成的集合，则q号传声器接收的声压信号p(ka,Ω_q)可写为：

$p(ka,{\mathrm{\Ω}}_q)={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)---\left(3\right)$

令p＝[p(ka,Ω₁),p(ka,Ω₂),…,p(ka,Ω_Q)]^T

t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^T分别为传声器声压信号列向量和声场传递函数列向量，上标“T”表示转置运算，对应式(3)，

$p={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)---\left(4\right)$

各传声器接收声压信号的互谱矩阵C写为：

$C=\stackrel{\‾}{{\mathrm{pp}}^H}={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈B}{\Σ}_{({r_0}^{\′},{{\mathrm{\Ω}}_0}^{\′})\∈B}\stackrel{\‾}{s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)s^{*}({{\mathrm{kr}}_0}^{\′},{{\mathrm{\Ω}}_0}^{\′})}t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({{\mathrm{kr}}_0}^{\′},{{\mathrm{\Ω}}_0}^{\′})---\left(5\right)$

其中，上标“-”、“H”分别表示平均运算和转置共轭运算，(r₀',Ω₀')也为声源位置坐标。

步骤2：获取实心球延迟求和(DAS)输出函数

过程201 定义实心球DAS输出函数

DAS输出函数为：

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\frac{1}{Q}\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)Cv({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2}---\left(6\right)$

其中，(r_f,Ω_f)为聚焦点的位置坐标，

v(kr_f,Ω_f)＝[v₁(kr_f,Ω_f),v₂(kr_f,Ω_f),…,v_Q(kr_f,Ω_f)]^T为聚焦列向量，|| ||₂表示2范数。v(kr_f,Ω_f)中，元素v_q(kr_f,Ω_f)表示q号传声器的聚焦分量，被定义为：

$v_q({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\underset{n=0}{\overset{N}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_f,ka)Y_n^{m*}\left({\mathrm{\Ω}}_f\right)Y_n^m\left({\mathrm{\Ω}}_q\right)---\left(7\right)$

其中，为Ω_f方向的球谐函数，R_n(kr_f,ka)为聚焦径向函数，其表达式与式(2)类同，只需将式(2)中的r₀替换为r_f即可，N为球谐函数阶次截断长度，需人为设定。式(6)也可写为：

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\underset{q,q^{\′}=1}{\overset{Q}{\Σ}}\frac{C_{{\mathrm{qq}}^{\′}}{v}_q^{*}({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)v_{q^{\′}}({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{Q{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2}---\left(8\right)$

其中，q'也为传声器索引，C_qq'为q号传声器接收声压信号相对于q'号传声器接收声压信号的互谱。

过程202 变形实心球DAS输出函数

定义平均声压贡献为声源在阵列各传声器处产生声压信号自谱的均值，记为P_AC(kr₀,Ω₀)，等于传声器接收声压信号互谱矩阵的迹平均。对于单声源，互谱矩阵C可写为：

C＝S(kr₀,Ω₀)t(kr₀,Ω₀)t^H(kr₀,Ω₀) (9)

其中，被定义为具有功率意义的声源强度。相应地，

$P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)=\frac{1}{Q}tr\left(C\right)=\frac{1}{Q}S({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2---\left(10\right)$

其中，“tr”表示求矩阵的迹。根据式(10)，用P_AC(kr₀,Ω₀)表示S(kr₀,Ω₀)后代入式(9)可得：

$C=\frac{P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)Q}{{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)---\left(11\right)$

将式(11)代入式(6)得：

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}---\left(12\right)$

psf被定义为实心球DAS声源识别方法的点传播函数，表示单位平均声压贡献单极子点声源的响应，根据式(12)该DAS方法的psf表达式为：

$psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))=\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}---\left(13\right)$

对单声源，该DAS方法的输出结果等于声源的平均声压贡献与相应点传播函数的乘积。

过程203 确定球谐函数阶次截断长度

球谐函数阶次截断长度按下式确定；

其中，符号表示将数值向正无穷方向圆整到最近的整数。

步骤3获取实心球函数型延迟求和(FDAS)输出函数

函数波束形成(FB)方法首先特征值分解阵列传声器接收声压信号的互谱矩阵C，得：

C＝UΣU^H (15)

其中，U＝[u₁,u₂,…,u_q,…,u_Q]为酉矩阵，Σ＝diag(σ₁,σ₂,…,σ_q,…,σ_Q)为对角矩阵，σ_q(q＝1,2,…,Q)为矩阵C的特征值，满足σ₁≤σ₂≤…≤σ_Q，u_q为σ_q对应的特征向量。然后，引入指数参数ξ，令：

$C^{\frac{1}{\mathrm{\ξ}}}\≡Udiag({\mathrm{\σ}}_1^{\frac{1}{\mathrm{\ξ}}},{\mathrm{\σ}}_2^{\frac{1}{\mathrm{\ξ}}},...,{\mathrm{\σ}}_q^{\frac{1}{\mathrm{\ξ}}},...,{\mathrm{\σ}}_Q^{\frac{1}{\mathrm{\ξ}}})U^H=\underset{q=1}{\overset{Q}{\Σ}}{\mathrm{\σ}}_q^{\frac{1}{\mathrm{\ξ}}}{u}_q{u}_q^H---\left(16\right)$

最后，采用已有声源识别方法后处理并计算所得结果的ξ次方。根据步骤2，构建的FDAS方法的输出函数为：

$\begin{array}{cc}{W}_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\frac{1}{Q}{\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)C^{\frac{1}{\mathrm{\ξ}}}v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2}\right)}^{\mathrm{\ξ}}& (\mathrm{\ξ}\≥2)\end{array}---\left(17\right)$

显然，当ξ＝1时，式(17)即为步骤2中DAS方法输出函数的表达式。假设三维空间内仅存在一个单极子点声源，平均声压贡献为P_AC(kr₀,Ω₀)，理论上，下列关系成立：

σ₁＝σ₂＝…＝σ_Q-1＝0 (18)

σ_Q＝tr(C)＝QP_AC(kr₀,Ω₀) (19)

$u_Q{u}_Q^H=\frac{t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)}{{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}---\left(20\right)$

联立式(16)～(20)可得：

$\begin{array}{c}{W}_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}\right)}^{\mathrm{\ξ}}\\ =P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\mathrm{psf}}^{\mathrm{\ξ}}\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))\end{array}---\left(21\right)$

表明：FDAS方法对单极子点声源的响应等于该声源的平均声压贡献与其点传播函数ξ次方的乘积。

步骤4对函数型延迟求和(FDAS)结果进行脊探测(RD)

若二元函数在某位置的输出值在其最大表面曲率方向上为极值，则该位置被称为脊。特定频率及聚焦距离下，FDAS算法的输出函数W_F(kr_f,Ω_f)是关于聚焦方向Ω_f＝(θ_f,φ_f)的二元函数，其在声源方向的输出值远大于在其它聚焦方向的输出值，声源方向为明显的极大值脊，可通过脊检测剔除主瓣内的非脊点，以缩减主瓣宽度、提高空间分辨率。根据定义，脊也是最大表面曲率方向上一阶导数符号改变的转折点。记L(Ω_f)＝|W_F(kr_f,Ω_f)|为脊检测对象，“||”表示取模运算，对每个聚焦方向，可按下列过程检测其是否为极大值脊：

过程401.构造L(Ω_f)的二阶偏导数矩阵，即黑塞矩阵：

$H=\left[\begin{array}{cc}\frac{{\∂}^{2}L}{\∂{\mathrm{\φ}}_f^2}& \frac{{\∂}^{2}L}{\∂{\mathrm{\φ}}_f\∂{\mathrm{\θ}}_f}\\ \frac{{\∂}^{2}L}{\∂{\mathrm{\θ}}_f\∂{\mathrm{\φ}}_f}& \frac{{\∂}^{2}L}{\∂{\mathrm{\θ}}_f^2}\end{array}\right]---\left(22\right)$

过程402.特征值分解矩阵Η，找出绝对值最大的特征值λ及相应的特征向量u_λ，其中u_λ指示最大表面曲率方向，输出值在u_λ方向上若为极大值则λ为负，若为极小值则λ为正。

过程403.计算标量ρ(Ω_f)为：

$\mathrm{\ρ}\left({\mathrm{\Ω}}_f\right)=-\frac{1}{2}sign\left(\mathrm{\λ}\right)|sign(\▿L({\mathrm{\Ω}}_f+{\mathrm{\χu}}_{\mathrm{\λ}})\·u_{\mathrm{\λ}})-sign(\▿L({\mathrm{\Ω}}_f-{\mathrm{\χu}}_{\mathrm{\λ}})\·u_{\mathrm{\λ}})|---\left(23\right)$

其中，“sign”为符号函数，为梯度算子，“．”表示求向量内积，χ为常数，用以控制检测空间精度，本发明取为0.2。当ρ(Ω_f)＝1时，聚焦方向为极大值脊，保留FDAS的输出值；反之，用零取代FDAS在该聚焦方向的输出值来剔除非脊点。

步骤5对延迟求和(DAS)结果进行反卷积(DAMAS)

记G为聚焦点总数，W为G×1维列向量，对应元素为各聚焦点的延迟求和输出量，由式(12)计算得出，A为G×G维的点传播函数矩阵，由式(13)计算得出，P_AC为G×1维列向量，对应元素为各聚焦点处声源的平均声压贡献，为未知非负数，对于不相干声源如下线性方程组成立：

W＝AP_AC (24)

(24)式采用高斯-赛德尔迭代方案求解P_AC，移除点传播函数的影响，能够显著衰减旁瓣并获得超高空间分辨率。初始化基于第l次迭代计算的结果进行第l+1次迭代的具体过程为：

过程501.计算残差r_e：

$\begin{array}{c}{r}_e^{\left(l\right)}({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈D}psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))P_{AC}^{(l+1)}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\\ +{\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈E}psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))P_{AC}^{\left(l\right)}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)-W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\end{array}---\left(25\right)$

其中，D为已完成l+1次迭代计算的声源点组成的集合，E为未完成l+1次迭代计算的声源点组成的集合，二者的并集包含所有声源。

过程502.确定

$\begin{array}{cc}{P}_{AC}^{(l+1)}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)=max\left(P_{AC}^{\left(l\right)}\right({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)-\frac{r_e^{\left(l\right)}({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))},0)& r_f=r_0,{\mathrm{\Ω}}_f={\mathrm{\Ω}}_0\end{array}$

按照该方案依次计算每个聚焦点，即可得出

直接用DAMAS清晰化DAS的成像结果需要考虑所有聚焦点，庞大的聚焦点数会使反卷积承受严重耗时问题。步骤4中对FDAS结果进行脊探测获取的脊中包含声源点，且声源点处的输出值相对很大，可默认一定动态范围外的脊和非脊处无声源，输出值为零，仅对动态范围内的脊进行DAMAS运算来获得超高空间分辨率，此时，G不再为聚焦点总数，而是参与运算的脊的数目，该数值很小，可使反卷积具有极高的效率。

步骤6定位声源位置并且能绘制声源成像图实现三维声场可视化。

根据步骤5求得的P_AC即能定位声源并能绘制声源成像图实现三维声场可视化。

附图说明

图1是本发明球阵列波束形成坐标系图；

图2是本发明测量布置示意图；

图3是DAMAS方法3000Hz单声源识别声源成像图；

图4是FDAS方法3000Hz单声源识别声源成像图；

图5是本发明方法3000Hz单声源识别成像图；

图6是DAMAS方法3000Hz不相干双声源识别声源成像图；

图7是FDAS方法3000Hz不相干双声源识别声源成像图；

图8是本发明方法3000Hz不相干双声源识别成像图；

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

实心球声源识别低旁瓣超高分辨率声学图像快速获取方法:

1)：构建测量***：

参见图1，一个处于三维空间内的实心球阵列；所述实心球阵列包括一个半径为a的实心球，以及分布在所述实心球表面的Q个传声器；这些传声器的编号为q，q＝1,2,…,Q；实施例中，Q＝36；所述实心球的球心作为原点，实心球阵列所处的三维空间内任意一点的位置坐标用(r,Ω)描述，r表示所描述位置与原点间的距离，表示所描述位置的方向，θ、分别为其仰角、方位角；(a,Ω_q)为q号传声器的位置坐标；作为验证，本实施例中，声源为稳态白噪声信号激励的扬声器，采用公司、半径为0.0975m、集成4958型传声器的36通道实心球阵列采样声压信号，扬声器声源到阵列中心的距离均为1m。

参见图2，为本实施例测量布置示意图。

2)获取声压

通过硬件采集***采集各个传声器接收声压信号p(ka,Ω_q)；其中，波数k＝2πf/c，声波频率为f，声速为c，(a,Ω_q)为q号传声器的位置坐标；

3)通过各个传声器接收声压信号，得到互谱矩阵C(可采用B&K公司的pulse软件获得)：

$C=\stackrel{\‾}{{\mathrm{pp}}^H},$

$p={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0),$

其中：B表示所有声源位置坐标组成的集合，(r₀,Ω₀)表示声源的位置坐标，s(kr₀,Ω₀)表示声源的强度，

t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^T，t_q(kr₀,Ω₀)表示声源到q号传声器的声场传递函数，

$t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≡\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_0,ka)Y_n^{m*}\left({\mathrm{\Ω}}_0\right)Y_n^m\left({\mathrm{\Ω}}_q\right),$

分别为Ω₀、Ω_q方向的球谐函数，n、m为球谐函数的阶次，

上标“-”、“H”、“*”和“T”分别表示平均运算、转置共轭运算、共轭运算和转置运算；

4)确定各个频率对应的球谐函数阶次截断长度：

符号表示将数值向正无穷方向圆整到最近的整数。

5)将互谱矩阵C和球谐函数阶次截断长度N分别代入FDAS和DAS输出式：

$W_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}\right)}^{\mathrm{\ξ}}$

和

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\frac{1}{Q}\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)Cv({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2}$

其中：P_AC(kr₀,Ω₀)为互谱矩阵C的迹，(r_f,Ω_f)为聚焦点的位置坐标，t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^T为传声器声压信号声场传递函数列向量，v(kr_f,Ω_f)＝[v₁(kr_f,Ω_f),v₂(kr_f,Ω_f),…,v_Q(kr_f,Ω_f)]^T，v_q(kr_f,Ω_f)表示q号传声器的聚焦分量，ξ为指数参数，推荐取值为16，|| ||₂表示2范数，上标“T”表示转置运算；

6)对W_F(kr_f,Ω_f)进行脊探测(RD)并剔除非脊点和30dB动态范围之外的脊点，然后针对剩余脊点对延迟求和结果W(kr_f,Ω_f)进行反卷积(DAMAS)求得P_AC。

7)根据P_AC绘制声源成像图，实现声场可视化，准确定位声源位置。

用DAMAS方法、FDAS方法和本发明方法分别识别3000Hz单声源，成像动态范围均为30dB，聚焦声源面被设定为与阵列同心的球面，仰角从0°到180°，方位角从0°到360°，间隔均为5°，共37×73个聚焦点，聚焦距离为1m，反卷积的迭代次数均取为5000。该扬声器声源平均声压贡献的实测值为49.1649dB。

参见图3，为DAMAS方法3000Hz单声源识别声源成像图，输出主瓣各聚焦点处输出值的线性叠加为49.1866dB，与实测值相差0.0217dB，具有超高分辨率并且旁瓣很小；

参见图4，为FDAS方法3000Hz单声源识别声源成像图,输出主瓣峰值为48.2983dB,与实测值相差较大为0.8666dB且分辨率较低；

参见图5，为本发明方法3000Hz单声源识别成像图,输出主瓣各聚焦点处输出值的线性叠加为49.1649dB,与实测值相差很小为0.0275dB，同时也具有超高分辨率并且旁瓣衰减能力极好；

参见图6，为DAMAS方法3000Hz不相干双声源识别声源成像图；

参见图7，为FDAS方法3000Hz不相干双声源识别声源成像图；

参见图8，为本发明方法3000Hz不相干双声源识别成像图；

显见，本发明方法对单声源及不相干声源均能快速准确地定位声源，成像图清晰干净。相对于直接对延迟求和结果进行反卷积(DAMAS)，本发明方法的计算时间仅为前者用时的千分之四，极大地提高了效率，同时与DAMAS方法具有同等的超高空间分辨及更优的旁瓣衰减能力；相对于函数型延迟求和(DAS)方法，本发明延续了前者的强旁瓣衰减能力同时分辨率获得极大提高，并且能更准确地量化声源。

Claims (1)

1.实心球声源识别低旁瓣超高分辨率声学图像快速获取方法，其特征在于：

1)：构建测量***：

一个处于三维空间内的所述实心球阵列；所述实心球阵列包括一个半径为a的实心球，以及分布在所述实心球表面的Q个所述传声器；这些传声器的编号为q，q＝1，2，…，Q；所述实心球的球心作为原点，实心球阵列所处的三维空间内任意一点的位置坐标用(r，Ω)描述，r表示所描述位置与原点间的距离，表示所描述位置的方向，θ、分别为其仰角、方位角；(a，Ω_q)为q号传声器的位置坐标；

2)获取传声器声压

各个传声器接收声压信号p(ka，Ω_q)；其中，波数k＝2πf/c，声波频率为f，声速为c，(a，Ω_q)为q号传声器的位置坐标；

3)通过各个传声器接收声压信号。得到互谱矩阵C：

$C=\stackrel{\‾}{{\mathrm{pp}}^H},$

$p={\mathrm{\Σ}}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0),$

其中：B表示所有声源位置坐标组成的集合，(r₀，Ω₀)表示声源的位置坐标，

s(kr₀，Ω₀)表示声源的强度，

t(kr₀，Ω₀)＝[t₁(kr₀，Ω₀)，t₂(kr₀，Ω₀)，…，t_Q(kr₀，Ω₀)]^T，t_q(kr₀，Ω₀)表示声源到q号传声器的声场传递函数，

$t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≡\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_0,ka)Y_n^{m*}\left({\mathrm{\Ω}}_0\right)Y_n^m\left({\mathrm{\Ω}}_q\right),$

分别为Ω₀、Ω_q方向的球谐函数，n、m为球谐函数的阶次，R_n(kr₀，ka)为声信号传播径向函数；

上标“-”、“H”、“*”和“T”分别表示平均运算、转置共轭运算、共轭运算和转置运算；

4)确定各个频率对应的球谐函数阶次截断长度：

表示将数值向正无穷方向圆整到最近的整数；

5)将互谱矩阵C和球谐函数阶次截断长度N分别代入：

$W_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2{||t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)||}_2^2}\right)}^{\mathrm{\ξ}}$

和

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\frac{1}{Q}\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)Cv({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{{||v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)||}_2^2}$

其中：P_AC(kr₀，Ω₀)为互谱矩阵C的的迹，(r_f，Ω_f)为聚焦点的位置坐标，v(kr_f，Ω_f)＝[v₁(kr_f，Ω_f)，v₂(kr_f，Ω_f)，…，v_Q(kr_f，Ω_f)]^T，v_q(kr_f，Ω_f)表示q号传声器的聚焦分量， 和分别为聚焦点方向和各麦克风方向的球谐函数，ξ为指数参数，||||₂表示2范数；

6)对W_F(kr_f，Ω_f)进行脊探测并剔除非脊点和设定的动态范围之外的脊点，然后针对剩余脊点对延迟求和结果W(kr_f，Ω_f)进行反卷积求得P_AC°

7)根据P_AC绘制声源成像图。