CN104330768A - 一种基于声矢量传感器的机动声源方位估计方法 - Google Patents

Abstract

本发明涉及一种基于声矢量传感器的机动声源方位估计方法，属于信号处理技术领域。本发明首先估计空域噪声协方差矩阵，通过噪声白化方法固化最大能量定位方法中的加权参数，以避免最优加权参数选择的一维搜索过程，提高了最大能量定向算法的估计精度，然后结合最大能量定向估计子输出和声源匀速运动的先验信息，在极坐标系下采用卡尔曼滤波技术进一步提高机动声源的方位估计精度。通过理论分析和仿真研究，本发明的基于声矢量传感器的机动声源方位估计与跟踪方法的估计精度优于原最大能量定向方法，并且由于采用了卡尔曼滤波，本发明的均方角度误差低于静态声源定位情况下的克拉美-罗下界。

Description

一种基于声矢量传感器的机动声源方位估计方法

技术领域

本发明涉及一种基于声矢量传感器的机动声源方位估计方法，属于信号处理技术领域。 

背景技术

在信号处理领域，波达方向(Direction of Arrival-DOA)估计是一个重要的研究课题，在导航、目标定位、波束形成方面都有着广泛的应用。1994年，Nehorai等将声学矢量传感器接收信号模型引入信号处理领域后，有关声学矢量传感器信号处理即成为研究的热点问题。与传统声压传感器仅感知声压信息相比，声学矢量传感器还可感知质点振速，增加了信息获取量，有望获得对声源状态更准确的认识。因而引起国内外研究者的关注，并进行了初步探索，取得系列研究成果。具体而言，Nehorai等研究了基于声学矢量传感器阵列DOA估计的CRLB；Hawkes等讨论了声学矢量传感器阵列的几何结构和传感器空间位置对参数估计性能的影响；顾陈、何劲等提出一种基于传播算子的声学矢量传感器阵列扩展孔径二维DOA估计算法。陈华伟，赵俊渭将宽带聚焦思想引入到了矢量传感器阵宽带处理，提出了基于矢量传感器阵的宽带相干信号子空间最优波束形成方法。此外，基于子空间的信号处理方法也均被应用到矢量传感器应用领域，如MUSIC算法以及ESPRIT算法。2011年Levin等提出基于一种基于单一矢量传感器的DOA估计方法，以期在空间非均匀高斯噪声环境中，通过加权最大能量梯度法估计声源方位，由于其最优权值需要通过一维搜索获取，该方法运算量较高。 

上述DOA估计的研究大多是以静止声源为研究对象，然而在实际工程应用中，对于DOA随着时间变化的机动声源目标，上述算法要重复对信号协方差矩阵进行特征值分解或奇异值分解，运算量极大，不能适用于实时性较高的DOA估计场合。因此，国内外部分学者尝试应用一些跟踪方向的滤波算法到声源定位与跟 踪领域，取得一定研究成果，但上述研究多局限于声压传感器阵，而基于矢量传感器阵列的声源定位与跟踪方法研究较少。 

发明内容

本发明的目的是提供一种基于声矢量传感器的机动声源方位估计方法，以解决目前声矢量传感器声源定位研究技术大多针对静止声源，在估计机动声源方位时估计精度低、运算量大、不适合实时处理的问题。 

本发明为解决上述技术问题而提供一种基于声矢量传感器的机动声源方位估计方法，该方法的步骤如下： 

1)利用改进最大能量梯度法获取初始时刻声源方向矢量和k时刻声源方向矢量

2)求出对应极坐标下的俯仰角与方位角初始状态向量 同时计算k时刻声源方向矢量对应极坐标下的俯仰角与方位角

3)根据卡尔曼滤波方法，获取k时刻的状态预测及其误差协方差矩阵P(k，k-1)， 

$\hat{X}(k,k-1)=\mathrm{\Φ}(k,k-1)X(k-1)$

P(k，k-1）＝Φ(k，k-1）P(k-1，k-1）Φ(k，k-1)^T+Γ(k-1)Q(k-1)Γ(k-1)^T

其中，为k时刻声源目标的状态向量，θ为声源信号入射向量的俯仰角，为声源的方位角，Φ(k，k-1)为k-1时刻至k时刻的状态转移矩阵，Γ(k-1)为系统噪声驱动阵，W(k-1)为系统激励噪声序列； 

4)利用k时刻先验信息对以上预测信息进行修正，从而得到k时刻波达方向估计及其误差协方差矩阵P(k，k)： 

$\hat{X}(k,k)=\hat{X}(k,k-1)+K\left(k\right)(Z\left(k\right)-H\hat{X}(k,k-1))$

P(k，k)＝P(k，k-1)-K(k)HP(k，k-1) 

其中 $H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],$ Z(k)为k时刻声源方向矢量对应的极坐标下的量测方程，K(k)为卡尔曼滤波器的增益。 

所述步骤1)中利用改进最大能量梯度法获取k时刻声源方向矢量的过程为： 

A.根据非均匀噪声协方差估计方法，估计入射声信号的噪声协方差矩阵 

$\hat{Q}=\left[\begin{array}{cc}{\mathrm{\δ}}_p^2& 0_{1\×3}\\ 0_{3\×1}& {\mathrm{\δ}}_v^2{I}_{3\×3}\end{array}\right]=\left(\begin{array}{ccc}{\hat{r}}_{11}-{\hat{x}}_2{\hat{x}}_4^{-1}{\hat{x}}_3& 0& 0_{1\×2}\\ 0& {\hat{r}}_{22}-{\hat{r}}_3{\hat{x}}_5^{-1}{\hat{x}}_6& 0_{1\×2}\\ 0_{2\×1}& 0_{2\×1}& \mathrm{diag}\{{\hat{r}}_4-{\hat{r}}_2{\hat{r}}_1^{-1}{\hat{r}}_3\}\end{array}\right);$

B.将声矢量传感器的接收信号进行预白化处理，得到预白化后的信号  ${\hat{Y}}_k\left(n\right),n=\mathrm{1,2},\·\·\·,N$

${\hat{Y}}_k\left(n\right)={\hat{Q}}^{-\frac{1}{2}}\left[\begin{array}{c}{Y}_k^p\left(n\right)\\ Y_k^v\left(n\right)\end{array}\right]=\left[\begin{array}{c}{\hat{Y}}_k^p\left(n\right)\\ {\hat{Y}}_k^v\left(n\right)\end{array}\right],n=\mathrm{1,2},\·\·\·,N;$

C.根据预白化后的信号计算预白化后振速信号与声压信号的协方差阵 和振速信号的协方差阵

$\begin{array}{c}{\hat{R}}_{\mathrm{vp}}=\frac{1}{N}\underset{n-1}{\overset{N}{\mathrm{\Σ}}}{\hat{Y}}_k^v\left(n\right){\hat{Y}}_k^p\left(n\right)\\ {\hat{R}}_{\mathrm{vv}}=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left[{\hat{Y}}_k^v\left(n\right)\right]{\left[{\hat{Y}}_k^p\left(n\right)\right]}^T\end{array};$

D.根据得到的振速信号与声压信号的协方差阵和振速信号的协方差阵 计算和

${\hat{u}}_0=q_0={\hat{R}}_{\mathrm{vp}}/\left|\right|{\hat{R}}_{\mathrm{vp}}\left|\right|$

其中μ为步长参数，为0.5。 

所述步骤3)中涉及的声源目标的运动状态方程可描述为： 

X(k)＝Φ(k，k-1)X(k-1)+Γ(k-1)W(k-1) 

其中，为k时刻声源目标的状态向量；Φ(k，k-1)为k-1时刻至k时刻的状态转移矩阵，Γ(k-1)为系统噪声驱动阵；W(k-1)为系统激励噪声序列。 

所述步骤1)得到的向量对应的极坐标下的量测方程可以描述为： 

其中

Z(k)的噪声为V(k)，其协方差矩阵 

$R\left(k\right)=E\left\{\right[Z\left(k\right)-H\tilde{X}\left(k\right)\left]{[Z\left(k\right)-H\tilde{X}\left(k\right)]}^T\right\}.$

本发明的有益效果是：本发明首先估计空域噪声协方差矩阵，通过噪声白化方法固化最大能量定位方法中的加权参数，以避免最优加权参数选择的一维搜索过程，提高了最大能量定向算法的估计精度，结合最大能量定向估计子输出和声源匀速运动的先验信息，在极坐标系下采用卡尔曼跟踪进一步提高机动声源的方位估计精度，通过理论分析和仿真研究，本发明的声矢量传感器生源方位估计方法的估计精度优于原最大能量定向方法，并且由于采用了声源运动方程信息，新方法的均方角度误差低于静态声源定位的克拉美-罗下界。 

附图说明

图1是本发明所采用的声音矢量传感器模型示意图； 

图2是噪声协方差中声压与振速域元素的估计示意图； 

图3是本发明实施例在k＝121时的的MASE示意图； 

图4是本发明实施例在k＝121时的的MASE示意图； 

图5是本发明实施例中三维空间的跟踪效果图； 

图6本发明中的俯仰估计的RMSE量测时间示意图； 

图7是本发明中的方位估计的RMSR量测时间示意图。 

具体实施方式

下面结合附图对本发明的具体实施方式作进一步的说明。 

矢量传感器由三个正交的偶极子振速传感器与一个单极子声压传感器组成。利用声矢量传感器可以同时获得空间内某一点处声压和振速的信息，该传感器模型如图1所示。考虑在某一时刻k，载率为f的远场声源S(满足平面波条件)。定义u_k为与波达方向相反的声源方向向量： 

其中，θ_k∈[-π/2，π/2]为声信号入射向量的俯仰角，为信号源的方位角，p_k为声源声压；v_k＝[v_xk，v_yk，v_zk]^T表示声源振速且： 

v_k＝-p_ku_k/(ρ₀c)                 (2) 

其中，ρ₀为传播介质密度，c为传播介质中声信号的传播速度。 

在大多数水下声源定位及跟踪场景下，声源做速率较低的平稳运动，因此可以时间顺序分别在各个时刻进行数据采集并做出相应的DOA估计。假设每一时刻测得N组数据，则k时刻声源S_k∈R_1×N，噪声ε_k∈R_4×N以及矢量传感器接收数据Y_k∈R_4×N分别为 

S_k＝[S_k(1)，…，S_k(N)]                              (3) 

${\mathrm{\ϵ}}_k=\left[\begin{array}{c}{\mathrm{\ϵ}}_k^p\\ {\mathrm{\ϵ}}_k^v\end{array}\right]=[{\mathrm{\ϵ}}_k\left(1\right),\·\·\·,{\mathrm{\ϵ}}_k\left(N\right)]---\left(4\right)$

$Y_k=\left[\begin{array}{c}{Y}_k^p\\ Y_k^v\end{array}\right]=[Y_k\left(1\right),\·\·\·,Y_k\left(N\right)]---\left(5\right)$

其中，分别为声压域和振速域的加性噪声，矢量传感器声压域和振速域接收矢量可表示为 

$Y_k^p=a_p\left({\mathrm{\Θ}}_k\right)S_k+{\mathrm{\ϵ}}_k^p---\left(6\right)$

$Y_k^v=a_v\left({\mathrm{\Θ}}_k\right)S_k+{\mathrm{\ϵ}}_k^v---\left(7\right)$

即 

Y_k＝a(Θ_k)S_k+ε_k                   (8) 

其中， $a\left({\mathrm{\Θ}}_k\right)=\left[\begin{array}{c}{a}_p\left({\mathrm{\Θ}}_k\right)\\ a_v\left({\mathrm{\Θ}}_k\right)\end{array}\right]=e^{-j2\mathrm{\πf}{\mathrm{\τ}}_k}\left[\begin{array}{c}1\\ -u_k/{\mathrm{\ρ}}_{0}c\end{array}\right],$ 且a(Θ_k)∈R_4×N，τ_k为源S到传感器的波达时间。若声源信号S_k与噪声为广义平稳的零均值非相关随机过程，且方差分别为则噪声统计特性可表示为： 

$E\left\{\left[\begin{array}{c}{\mathrm{\ϵ}}_k^p\\ {\mathrm{\ϵ}}_k^v\end{array}\right]\right\}=0_{4\×N}---\left(9\right)$

$E\left\{\left[\begin{array}{c}{\mathrm{\ϵ}}_k^p\\ {\mathrm{\ϵ}}_k^v\end{array}\right]{\left[\begin{array}{c}{\mathrm{\ϵ}}_k^p\\ {\mathrm{\ϵ}}_k^v\end{array}\right]}^T\right\}=\left[\begin{array}{cc}{\mathrm{\δ}}_p^2& 0_{1\×3}\\ 0_{3\×1}& {\mathrm{\δ}}_v^2{I}_{3\×3}\end{array}\right]---\left(10\right)$

本发明的声矢量传感器声源方位估计方法通过白化单极子传感器和偶极子传感器接收噪声，然后结合卡尔曼滤波方法，利用声源运动方程先验信息，获取声源状态(方位、俯仰及速度)更精确的估计，该方法的具体过程如下： 

1.利用最大能量梯度算法，针对k时刻噪声分量高斯非均匀的情况，得到声源导向矢量的估计为： 

$\begin{array}{c}{\hat{u}}_k=\underset{q}{\arg \max }\left\{\frac{1}{N}\underset{n-1}{\overset{N-1}{\mathrm{\Σ}}}{[{\mathrm{\α}}_p{Y}_k^p\left(n\right)+{\mathrm{\α}}_v{q}^T{Y}_k^v\left(n\right)]}^2\right\}\\ \mathrm{subject\; to}{q}^{T}q=1\end{array}---\left(11\right)$

式(11)中，参数α_p、α_v与声压功率、振速功率在各个坐标轴上的分量有关，且满足关系α_p+α_v＝1。将公式(11)中多项式展开，求q的梯度，并对做适当变换得 

$\begin{array}{c}{\▿}_{q}T\left(q\right)={\mathrm{\α}}_p{R}_{\mathrm{vp}}+{\mathrm{\α}}_v{R}_{\mathrm{vv}}q\\ ={\mathrm{\α}}_p{R}_{\mathrm{vp}}+(1-{\mathrm{\α}}_p)R_{\mathrm{vv}}q\end{array}---\left(12\right)$

其中，R_vp为振速信号与声压信号的协方差阵，R_vv为振速信号的协方差阵。由式子(11)得出，与单极子元素部分的波束形成输出相关，与各个偶极子元素的波束形成输出相关，在条件下，即当α_p/α_v的值与单极子和各个偶极子噪声功率呈反比时，DOA估计均方误差接近Cramer-Rao下界。 

2.将传感器接收信号的协方差矩阵由若干次快拍采样估计，表示为 

$R=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left[\begin{array}{c}{Y}_k^p\left(n\right)\\ Y_k^v\left(n\right)\end{array}\right]{\left[\begin{array}{c}{Y}_k^p\left(n\right)\\ Y_k^v\left(n\right)\end{array}\right]}^T---\left(13\right)$

定义过渡矩阵D＝[I_3×3 0_3×1]，得到 

$R_D={\mathrm{DRD}}^T=\left(\begin{array}{ccc}{\hat{r}}_{11}& {\hat{x}}_1& {\hat{x}}_2\\ {\hat{x}}_3& {\hat{r}}_{22}& {\hat{x}}_4\\ {\hat{x}}_5& {\hat{x}}_6& {\hat{r}}_{33}\end{array}\right)---\left(14\right)$

将协方差阵R划分为如下形式 

$R=\left(\begin{array}{ccc}\mathrm{zz}& \mathrm{zz}& {\mathrm{zz}}_{1\×2}\\ {\hat{r}}_1& \mathrm{zz}& {\hat{r}}_{3_{1\×2}}\\ {\hat{r}}_{2_{2\×1}}& {\mathrm{zz}}_{2\×1}& {\hat{r}}_{4_{2\×2}}\end{array}\right)---\left(15\right)$

式中，zz表示在本讨论中无关的块矩阵，通过非均匀噪声协方差估计方法^[12]，利用(14)、(15)两式、估计入射声信号的噪声协方差矩阵 

$\hat{Q}=\left[\begin{array}{cc}{\mathrm{\δ}}_p^2& 0_{1\×3}\\ 0_{3\×1}& {\mathrm{\δ}}_v^2{I}_{3\×3}\end{array}\right]=\left(\begin{array}{ccc}{\hat{r}}_{11}-{\hat{x}}_2{\hat{x}}_4^{-1}{\hat{x}}_3& 0& 0_{1\×2}\\ 0& {\hat{r}}_{22}-{\hat{r}}_3{\hat{x}}_5^{-1}{\hat{x}}_6& 0_{1\×2}\\ 0_{2\×1}& 0_{2\×1}& \mathrm{diag}\{{\hat{r}}_4-{\hat{r}}_2{\hat{r}}_1^{-1}{\hat{r}}_3\}\end{array}\right)---\left(16\right)$

将传感器接收信号预白化得到 

${\hat{Y}}_k\left(n\right)={\hat{Q}}^{-\frac{1}{2}}\left[\begin{array}{c}{Y}_k^p\left(n\right)\\ Y_k^v\left(n\right)\end{array}\right]=\left[\begin{array}{c}{\hat{Y}}_k^p\left(n\right)\\ {\hat{Y}}_k^v\left(n\right)\end{array}\right],n=\mathrm{1,2},\·\·\·,N---\left(17\right)$

从而，预白化后振速信号与声压信号的协方差阵、振速信号的协方差阵分别为 

$\begin{array}{c}{\hat{R}}_{\mathrm{vp}}=\frac{1}{N}\underset{n-1}{\overset{N}{\mathrm{\Σ}}}{\hat{Y}}_k^v\left(n\right){\hat{Y}}_k^p\left(n\right)\\ {\hat{R}}_{\mathrm{vv}}=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left[{\hat{Y}}_k^v\left(n\right)\right]{\left[{\hat{Y}}_k^p\left(n\right)\right]}^T\end{array}---\left(18\right)$

3.改进的最大能量梯度DOA估计算法为静态声源定位方法，如将其直接应用于动态目标跟踪场景，由于没有利用声源运动方程的先验信息，所以计算量大且精度有待提高。考虑声源跟踪的水域舰船探测背景，声源状态常见为匀速，而运动模型则为线性模型，因此可引入卡尔曼滤波算法到运动声源的DOA估计中，以提高声源位置估计精度，卡尔曼滤波算法表示如下 

$\hat{X}(k,k-1)=\mathrm{\Φ}(k,k-1)X(k-1)---\left(19\right)$

P(k，k-1)＝Φ(k，k-1)P(k-1，k-1)Φ(k，k-1)^T+Γ(k-1)Q(k-1)Γ(k-1)^T    (20) 

K(k)＝P(k，k-1)H^T(HP(k，k-1)H^T+R(k))^-1               (21) 

$\hat{X}(k,k)=\hat{X}(k,k-1)+K\left(k\right)(Z\left(k\right)-H\hat{X}(k,k-1))---\left(22\right)$

P(k，k)＝P(k，k-1)-K(k)HP(k，k-1)           (23) 

其中 $H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],$ Z(k)为k时刻声源方向矢量对应的极坐标下的量测方程，K(k)为卡尔曼滤波器的增益，为k时刻的状态更新值。 

声矢量传感器DOA估计的研究中，传感器观测模型为极坐标模型，而在目标跟踪算法中如果采用空间直角坐标系则会造成目标状态与观测之间的非线性关系，不利于数据处理。因此需要针对机动声源目标建立极坐标系下的运动模型，设声源初始位置为并以的角速度做匀速运动。 

则声源目标的运动状态方程可描述为 

X(k)＝Φ(k，k-1)X(k-1)+Γ(k-1)W(k-1)            (24) 

其中，为k时刻声源目标的状态向量；Φ(k，k-1)为k-1时刻至k时刻的状态转移矩阵，Γ(k-1)为系统噪声驱动阵；W(k-1)为系统激励噪声序列。以改进的最大能量梯度法的输出向量作为卡尔曼滤波方程的量测信息。 

则向量对应的极坐标下的量测方程可以描述为 

其中

Z(k)的噪声为V(k)，其协方差矩阵 

$R\left(k\right)=E\left\{\right[Z\left(k\right)-H\tilde{X}\left(k\right)\left]{[Z\left(k\right)-H\tilde{X}\left(k\right)]}^T\right\}---\left(26\right)$

式中， $H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right].$ 然而在滤波过程中，k时刻的状态向量的实际值未知，因而无法采用式(26)获得R(k)，当量测噪声ε_p(k)、ε_v(k)是方差为目互不相关的高斯白噪声时，采用改进的最大能量梯度DOA估计算法获得DOA估计误差基本处于稳定，因而可以采用统计学的方法，通过蒙特卡洛算法，基于每一组快拍数据得到DOA估计与其整体均值的误差，进而得出的k时刻的Z(k)的噪声 协方差矩阵。 

$R\left(k\right)=\left\{\underset{l=1}{\overset{N}{\mathrm{\Σ}}}\right[Z_l\left(k\right)-\underset{m=1}{\overset{N}{\mathrm{\Σ}}}{Z}_m\left(k\right)/N]/N\}{\left\{\underset{l=1}{\overset{N}{\mathrm{\Σ}}}\right[Z_l\left(k\right)-\underset{m=1}{\overset{N}{\mathrm{\Σ}}}{Z}_m\left(k\right)/N]/N\}}^T---\left(27\right)$

综上所述，算法将DOA估计的俯仰和方位信息作为卡尔曼滤波算法的量测信息，自适应地跟踪移动声源的波达方向，具体算法流程如下表所示： 

仿真分析 

本发明结合卡尔曼滤波算法和改进的最大能量算法(IMP)在三维空间定位跟踪运动声源，以矢量传感器的位置为极坐标原点，设定声源为匀速，假定初始时刻，声源与传感器的距离d＝1km，俯仰角θ₀＝25°，方位角在俯仰方向和方位方向的声源速度分别为0.25°/second和0.2°/second。在噪声功率时变的情况下，生成具有不同信噪比的三组信号和噪声，在第k跟踪时间步可以将它们区分出来。信号噪声方差分别设置如下：(1)时间步k＝1～60， ${\mathrm{\δ}}_{v,k}^2=1.3;$ (2)时间步k＝61～120， ${\mathrm{\δ}}_{s,k}^2=10,{\mathrm{\δ}}_{p,k}^2=0.9,{\mathrm{\δ}}_{v,k}^2=1.1;$ (3)时间步k＝121～180，运动声源的运动状态方程由公式(3)给出。 

在第k步，设定过程噪声W_k-1为零均值的高斯白噪声，其正定协方差矩阵  $Q_{k-1}=\left[\begin{array}{cc}{\mathrm{\σ}}^2& 0\\ 0& {\mathrm{\σ}}^2\end{array}\right],$ 其中σ＝0.01°。过程噪声分布矩阵 ${\mathrm{\Γ}}_{k-1}=\left[\begin{array}{cc}F& 0\\ 0& F\end{array}\right],$ 其中  $F=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right],$ 采样间隔T＝1。 

根据声压与振速域信号的预白化协方差矩阵，采用IMP算法进行DOA估计，所以得到噪声协方差Ω_k的精确估计是非常重要的，在以上仿真条件下，设定总的跟踪时间步L＝180，每个跟踪时间步k采取N＝8000次快拍采样，通过公式(11)～(14)估计主对角阵Ω_k的元素。噪声协方差矩阵的单极子和偶极子噪声功率估计如图2所示。 

同时对IMP算法的DOA估计性能进行了评估，该仿真试验采用均方角度误差(MSAE)为标准评估DOA估计算法性能，采用蒙特卡洛算法提高仿真结果的代表性和说服力，MASE及其克拉美-罗界(CRLB)具体表达式如下： 

$\mathrm{MSAE}=\left[\underset{m=1}{\overset{\mathrm{MC}}{\mathrm{\Σ}}}\underset{N\→\∞}{\lim }(N\×E\{{{\mathrm{AE}}_m}^2\left\}\right)\right]/\mathrm{MC},{\mathrm{AE}}_m={2\sin }^{-1}\left(\frac{\left|\right|{\hat{u}}_{m,k}-u_k\left|\right|}{2}\right)---\left(28\right)$

${\mathrm{MSAE}}_{\mathrm{CRLB}}=\frac{{\mathrm{\δ}}_{v,k}^2}{{\mathrm{\δ}}_{s,k}^2}(1+\frac{\left({\mathrm{\δ}}_{p\left|\right|v,k}^2\right)}{{\mathrm{\δ}}_{s,k}^2}),{\mathrm{\δ}}_{p\left|\right|v,k}^2={({\mathrm{\δ}}_{p,k}^{-2}+{\mathrm{\δ}}_{v,k}^{-2})}^{-1}---\left(29\right)$

其中蒙特卡洛仿真次数MC＝1000，在第m次蒙特卡洛仿真中，与u_k分别表 示第k步由矢量传感器指向声源的估计和真实矢量。 

该仿真实验中，在相同的仿真环境下利用MP算法和IMP算法进行移动声源DOA估计，分别设定k＝1和k＝121作为对照研究。在图3和图4中，将MP算法的均方角度误差(MSAE)表示为加权参数的函数，其中取值范围为0-1，且MP算法的最小均方角度误差已在图中标出。IMP算法的均方角度误差(MSAE)和克拉美-罗界(CRLB)如图3和图4所示，从这两个图中可以发现，在不同信噪比条件下，IMP算法的DOA估计均方角度误差接近MP算法的最小均方角度误差，并接近克拉美-罗界，这说明在未知噪声环境下IMP算法可以达到MP算法的最高估计精度。 

利用IMP算法和卡尔曼滤波的移动声源DOA估计 

以IMP算法得到的DOA估计结果作为量测，正如公式(20)所提到的， 卡尔曼滤波量测Z_k的噪声协方差矩阵R_k可利用蒙特卡洛仿真实验通过公式(22)得到。信息模型，仿真环境和相关参数与上文相同，图5展示了两种算法的跟踪效果。图6和图7分别展示了方位估计和俯仰估计的RMSE。 

从图5-7中可以发现，利用定位跟踪联合算法可以有效地跟踪运动声源，与MP算法相比，其位置估计RMSE较低。为了进一步说明在不同环境下定位跟踪联合算法的优越性，表2对两种算法的RMSE时间均值进行了对比，我们发现新算法的估计精度比MP DOA估计算法有所提高，且其MSAE低于静态CRLB，这是因为定位跟踪联合算法有效运用先验信息对初始估计进行校正。 

Claims (4)

1.一种基于声矢量传感器的机动声源方位估计方法，其特征在于，该方法的步骤如下：

1)利用改进最大能量梯度法获取初始时刻声源方向矢量和k时刻声源方向矢量

2)求出对应极坐标下的俯仰角与方位角初始状态向量同时计算k时刻声源方向矢量对应极坐标下的俯仰角与方位角

3)根据卡尔曼滤波方法，获取k时刻的状态预测及其误差协方差矩阵P(k,k-1)，

$\hat{X}(k,k-1)=\mathrm{\Φ}(k,k-1)X(k-1)$

P(k,k-1)＝Φ(k,k-1)P(k-1,k-1)Φ(k,k-1)^Τ+Γ(k-1)Q(k-1)Γ(k-1)^Τ

其中，为k时刻声源目标的状态向量，θ为声源信号入射向量的俯仰角，为声源的方位角，Φ(k,k-1)为k-1时刻至k时刻的状态转移矩阵，Γ(k-1)为系统噪声驱动阵，W(k-1)为系统激励噪声序列；

4)利用k时刻先验信息对以上预测信息进行修正，从而得到k时刻波达方向估计及其误差协方差矩阵P(k,k)：

$\hat{X}(k,k)=\hat{X}(k,k-1)+K\left(k\right)(Z\left(k\right)-H\hat{X}(k,k-1))$

P(k,k)＝P(k,k-1)-K(k)HP(k,k-1)

其中 $H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],$ Z(k)为k时刻声源方向矢量对应的极坐标下的量测方程，K(k)为卡尔曼滤波器的增益。

2.根据权利要求1所述的基于声矢量传感器的机动声源方位估计方法，其特征在于，所述步骤1)中利用改进最大能量梯度法获取k时刻声源方向矢量的过程为：

A.根据非均匀噪声协方差估计方法，估计入射声信号的噪声协方差矩阵

$\hat{Q}=\left[\begin{array}{cc}{\mathrm{\δ}}_p^2& 0_{1\×3}\\ 0_{3\×1}& {\mathrm{\δ}}_v^2{I}_{3\×3}\end{array}\right]=\left(\begin{array}{ccc}{\hat{r}}_{11}-{\hat{x}}_2{\hat{x}}_4^{-1}{\hat{x}}_3& 0& 0_{1\×2}\\ 0& {\hat{r}}_{22}-{\hat{x}}_3{\hat{x}}_5^{-1}{\hat{x}}_6& 0_{1\×2}\\ 0_{2\×1}& 0_{2\×1}& \mathrm{diag}\{{\hat{r}}_4-{\hat{r}}_2{\hat{r}}_1^{-1}{\hat{r}}_3\}\end{array}\right);$

B.将声矢量传感器的接收信号进行预白化处理，得到预白化后的信号 ${\hat{Y}}_k\left(n\right),n=\mathrm{1,2},...,N$

${\hat{Y}}_k\left(n\right)={\hat{Q}}^{-\frac{1}{2}}\left[\begin{array}{c}{Y}_k^p\left(n\right)\\ Y_k^v\left(n\right)\end{array}\right]=\left[\begin{array}{c}{\hat{Y}}_k^p\left(n\right)\\ {\hat{Y}}_k^v\left(n\right)\end{array}\right]n=\mathrm{1,2},...,N;$

C.根据预白化后的信号计算预白化后振速信号与声压信号的协方差阵和振速信号的协方差阵

$\begin{array}{c}{\hat{R}}_{\mathrm{vp}}=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\hat{Y}}_k^v\left(n\right){\hat{Y}}_k^p\left(n\right)\\ {\hat{R}}_{\mathrm{vv}}=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left[{\hat{Y}}_k^v\left(n\right)\right]{\left[{\hat{Y}}_k^p\left(n\right)\right]}^T\end{array};$

D.根据得到的振速信号与声压信号的协方差阵和振速信号的协方差阵计算和

${\hat{u}}_0=q_0={\hat{R}}_{\mathrm{vp}}/\left|\right|{\hat{R}}_{\mathrm{vp}}\left|\right|$

其中μ为步长参数，为0.5。

3.根据权利要求1所述的基于声矢量传感器的机动声源方位估计方法，其特征在于，所述步骤3)中涉及的声源目标的运动状态方程可描述为：

X(k)＝Φ(k,k-1)X(k-1)+Γ(k-1)W(k-1)

其中，为k时刻声源目标的状态向量；Φ(k,k-1)为k-1时刻至k时刻的状态转移矩阵，Γ(k-1)为系统噪声驱动阵；W(k-1)为系统激励噪声序列。

4.根据权利要求1所述的基于声矢量传感器的机动声源方位估计方法，其特征在于，所述步骤1)得到的向量对应的极坐标下的量测方程可以描述为：

其中

Z(k)的噪声为V(k)，其协方差矩阵

$R\left(k\right)=E\left\{\right[z\left(k\right)-H\tilde{X}\left(k\right)\left]{[z\left(k\right)-H\tilde{X}\left(k\right)]}^T\right\}.$