CN101272168A - 一种信源数估计方法及其波达方向估计方法 - Google Patents

Abstract

本发明公开了一种波达方向估计中基于特征空间的信源数估计方法，包括如下步骤：(1)获得均匀线阵所接收的信号；(2)由接收信号计算协方差矩阵；(3)求复数的或者实数化后的协方差矩阵；(4)将协方差矩阵进行特征值分解，得到它的特征矢量矩阵；(5)由特征矢量求信源数判据；(6)求比例值；(7)信源数估计。本发明能够准确估计信源数并能够为波达方向估计信号处理过程省去大量的计算量，减少硬件开销。

Description

一种信源数估计方法及其波达方向估计方法

技术领域

本发明涉及波达方向估计中的信号处理方法，特别涉及一种基于特征空间的信源数估计方法，以及采用该信源数估计方法的波达方向估计方法。

背景技术

波达方向估计中的大多数处理方法是基于信号子空间算法，在这种处理方法中要求信号源数目(简称信源数)是已知的，则信号子空间与噪声子空间相互正交，由此得到正确的波达方向估计。但是在实际情形下，进行波达方向估计时往往不知道信源数，所以就需要首先估计信源数，如果信源数目估计错误，则造成信号子空间与噪声子空间不再正交，从而造成信源的漏警或虚警以及波达方向估计结果的偏差。因此信源数估计是波达方向估计中的关键技术。

王永宁在“空间谱估计理论与算法，北京：清华大学出版社，2004，11”一书中认为H.T.Wu等人提出的信源数的盖氏圆估计法是好方法。H.T.Wu等人的文章“Source number estimator using Gerschgorin disks，Proc.ICASSP，Adelaide，Australia，Apr.1994，261-264”，和“Source number estimation usingtransformed Gerschgorin radii，IEEE Trans.on SP，1995，43(6)：1325-1333”，以及“Gerschgorin radii based source number detection for closely spacedsignals，Proc.IEEE int.Conf.Acoust.Speech，Signal processing.Atlanta，May，1996，3054-3057”提出了信源数的盖氏圆盘估计法，该方法的基本技术方案为不需要具体知道信号协方差矩阵的特征值的数值，只是利用盖氏圆盘定理估计出各特征值的位置，信号的盖氏圆盘半径大，噪声的盖氏圆盘半径小，由半径较大的盖氏圆盘数目估计信源数。在上述基础上，H.T.Wu等人提出了一种启发式的判据，区分出信号和噪声的特征值，由此估计信源数。但是在信噪比较低时，常常会产生欠估计和过估计。而且在空间相关色噪声的相关系数较大时，性能退化。另外，在盖氏圆盘估计法中，要对(M-1)×(M-1)阶信号协方差矩阵进行特征值分解，其中M是传感器基元数，而在波达方向估计中要对M×M阶信号协方差矩阵进行特征值分解。因此，为了得到正确的波达方向，如果采用盖氏圆盘估计法进行信源数的估计，就不可避免地如刚才所述的要进行两次特征值分解，这耗费了大量的计算量，增加了硬件的开销。再次，它难以在实数空间中进行运算，不能获得实数空间中运算量减小的好处。

综上所述，鉴于现有的盖氏圆盘估计信源数的方法存在不足，因此就需要一种新的信源数估计方法，它的估计精度高、运算量小和鲁棒性好。

发明内容

本发明的目的是克服现有技术不足，提供适用于波达方向估计技术的一种具有估计精度高、运算量小和鲁棒性好的信源数估计方法。

为了达到上述目的，根据本发明的第一方面，提供一种信源数估计方法，包括如下步骤：

(1)获得均匀线阵所接收的信号X(t)，X(t)＝[x₁(t)，x(t)，…x_M(t)]^T为M×1的矢量，其中x₁，...x_M表示所述均匀线阵中基元1到M接收的信号，t表示时刻；

(2)由接收信号计算协方差矩阵

(3)将

进行特征值分解，得到它的特征矢量矩阵U_C＝[q₁，q₂…q_M]，其中q_i是

的特征矢量，q_i＝[q_1i，q_2i…q_mi…q_Mi]^T，i＝1，2，…M；

(4)根据q_i计算信源数m。

根据本发明的第二个方面，在所述方法的步骤(2)中，利用公式 ${\hat{R}}_X=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}X\left(t_l\right)X{\left(t_l\right)}^H,$ 进行协方差矩阵计算，其中L为快拍数。

根据发明的第三个方面，在所述方法的步骤(2)之后还包括步骤：

(a)由协方差矩阵

计算矩阵

其中为对整个阵列的接收数据作一次前后向空间平滑后得到的协方差矩阵；

(b)由

求得实数化处理后的协方差矩阵 ${\hat{R}}_T=Q_M^H{\hat{R}}_{\mathrm{fb}}{Q}_M;$

以及在步骤(3)中包括：使 ${\hat{R}}_X={\hat{R}}_T,$ 并将进行特征值分解，得到它的特征矢量矩阵 $U_R=Q_M^H{U}_C,$ U_R＝[q₁，q₂，…q_M]，其中q_i是

的特征矢量，q_i＝[q_1i，q_2i…q_mi…q_Mi]^T，i＝1，2，…M。

根据发明的第四个方面，所述方法的所述步骤(4)包括下列步骤：

(4.1)利用公式 $\stackrel{\‾}{\left|q_m\right|}=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\left|q_{\mathrm{im}}\right|$ 求判据值|q_m|，m＝1，2……M；

(4.2)利用判据值|q_m|确定信源数m。

根据发明的第五个方面，所述方法的步骤(4.2)还包括步骤：

(4.3)根据公式β_m＝|q_m|²/|q_m+1|求比例值β_m，其中m＝1，2……M；

(4.4)求β_m最大值对应的m，m＝1，2……M，m是信源数。

根据发明的第六个方面，所述方法的所述步骤(3)中q_i是按照其对应特征值大小排列的，i＝1，2，…M。

根据发明的第七个方面，提供一种波达方向估计方法，包括步骤：

(1)均匀线阵接收阵列信号X(t)，由上述第1至第6方面之一的信源数估计方法计算出信源数K；

(2)取得上述第1至第6方面之一的协方差矩阵估计值

，它是一个M×M阶的矩阵；

(3)采用上述

进行波达方向估计，获得波的入射方向。

根据本发明的第八个方面，本发明还提供一种信源数估计装置，包括：

(1)均匀线阵，包括基元1到M，用于获得均匀线阵所接收的信号X(t)，X(t)＝[x₁(t)，x(t)，…x_M(t)]^T为M×1的矢量，其中x₁，...x_M表示所述均匀线阵中基元1到M接收的信号，t表示时刻；

(2)协方差矩阵计算单元，用于由接收信号计算协方差矩阵

(3)特征值分解单元，用于将

进行特征值分解，得到它的特征矢量矩阵U_C＝[q₁，q₂…q_M]，其中q_i是

的特征矢量，q_i＝[q_1i，q_2i…q_mi…q_Mi]^T，i＝1，2，…M；

(4)信源数计算单元，根据q_i计算信源数m。

根据本发明的第九个方面，所述协方差矩阵计算单元利用公式 ${\hat{R}}_X=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}X\left(t_l\right)X{\left(t_l\right)}^H,$ 进行协方差矩阵计算，其中L为快拍数。

根据本发明的第十个方面，所述信源数估计装置还包括下列单元：

用于由协方差矩阵计算矩阵

的单元，其中

为对整个阵列的接收数据作一次前后向空间平滑后得到的协方差矩阵；

用于由

求得实数化处理后的防方差矩阵 ${\hat{R}}_T=Q_M^H{\hat{R}}_{\mathrm{fb}}{Q}_M$ 的单元；

以及所述特征值分解单元还进行下列计算：使 ${\hat{R}}_X={\hat{R}}_T,$ 并将

进行特征值分解，得到它的特征矢量矩阵 $U_R=Q_M^H{U}_C,$ U_R＝[q₁，q₂，…q_M]，其中q_i是的特征矢量，q_i＝[q_1i，q_2i…q_mi…q_Mi]^T，i＝1，2，…M。

根据本发明的第十一个方面，所述信源数计算单元还包括：

利用公式 $\stackrel{\‾}{\left|q_m\right|}=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\left|q_{\mathrm{im}}\right|$ 求判据值|q_m|的单元，m＝1，2……M；以及

利用判据值|q_m|确定信源数m的单元。

根据本发明的第十二个方面，所述利用判据值|q_m|确定信源数m的单元还包括：

根据公式β_m＝|q_m|²/|q_m+1|求比例值β_m的单元，其中m＝1，2……M；

求β_m最大值对应的m的单元，m＝1，2……M，m是信源数。

与现有技术相比，本发明的有益效果在于：

1.本发明应用了信号子空间与噪声子空间的正交性，将协方差矩阵

降噪后分别直接投影到信号子空间和噪声子空间，以利于准确估计信源数。

2.本发明是对M×M阶协方差矩阵进行特征值分解，它可以直接用于以后的波达方向估计中，无需再进行一次特征值分解。由于在波达方向估计中，特征值分解的运算量占总运算量至少一半，因此利用本发明能够为波达方向估计信号处理过程省去大量的计算量，减少硬件开销。

3.本发明的方法可以在实数空间中运算，进一步减少运算量。

附图说明

图1是本发明的信源数估计方法流程图；

图2是本发明的一个应用实施例与现有技术的盖氏圆盘估计法性能比较，其中图2(a)是采用本发明的结果，它是实数特征空间法。图2(b)是采用本发明的结果，它是复数特征空间法，简称为特征空间法。图2(c)是采用现有技术的盖氏圆盘估计法的结果。

图3是信源功率相等情况下，实数特征空间法、复数特征空间法和盖氏圆信源数估计的正确概率随信噪比的变化。

图4是信源功率不等情况下，实数特征空间法、复数特征空间法和盖氏圆信源数估计的正确概率随信噪比的变化。

图5是空间相关色噪声情况下，实数特征空间法、复数特征空间法和盖氏圆信源数估计的正确概率比较。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细描述：

为了便于理解本发明，下面首先从原理上来阐述本发明。

一.特征空间的信源数估计方法

1.1理论

设一均匀线阵由M个基元组成，用来接收信号源发出的信号，设接收到K个远场窄带信号(即信号源的个数为K)。设各个信号源发出的信号对于均匀线阵的入射角为(θ₁，θ₂，…θ_K)，则线阵接收的阵列信号为：

X(t)＝AS(t)+N(t)                  (1)

其中t是时间，各项为：

X(t)＝[x₁(t)，x₂(t)，…x_M(t)]^T，是M×1的矢量       (2a)

A＝[a(θ₁)，a(θ₂)，…a(θ_K)]，为M×K的阵列导向矩阵(2b)

$a\left({\mathrm{\θ}}_k\right)=[1,e^{-j\frac{\mathrm{\ω}d\sin {\mathrm{\θ}}_k}{c}},...e^{-j(M-1)\mathrm{\ω}\frac{d\sin {\mathrm{\θ}}_k}{c}}],$ 称为导向矢量，k＝1，…K(2c)S(t)＝[s₁(t)，s₂(t)，s₃(t)…s_K(t)]^T，是K个窄带信号    (2d)

N(t)是M×1的噪声矢量                              (2e)

其中 $j=\sqrt{-1},c$ 表示介质中的声速，ω是信号的中心频率，d是相邻基元间距，上标T表示矩阵转置。A还可以写成另一个形式：

A＝[b₁，b₂，…b_M]^T                                (3a)

$b_i={[e^{-j(i-1)\mathrm{\ω}\frac{d\sin {\mathrm{\θ}}_1}{c}},...e^{-j(i-1)\mathrm{\ω}\frac{d\sin {\mathrm{\θ}}_K}{c}}]}^T,i=\mathrm{1,2},...M---\left(3b\right)$

其中b_i是A的行矢量。设信号S与噪声不相关，当噪声为空间白噪声的情况下，阵列信号协方差矩阵可写成

R_X＝E[x(t)x^H(t)]

(4)

＝AR_SA^H+σ²I

其中R_S为信号S的协方差矩阵，σ²是噪声的方差，I是单位矩阵。实际情况下，往往得不到R_X，常用下面的估计值近似：

${\hat{R}}_X=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}X\left(t_l\right)X{\left(t_l\right)}^H---\left(5\right)$

其中L为快拍数，上标H表示矩阵共轭转置，由此得

${\hat{R}}_X=[A{\hat{R}}_S{b}_1^H,A{\hat{R}}_S{b}_2^H,...A{\hat{R}}_S{b}_M^H]+{\mathrm{\σ}}^{2}I---\left(6\right)$

其中是信号的协方差矩阵估计值，其表式为：

${\hat{R}}_S=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}S\left(t_l\right)S{\left(t_l\right)}^H---\left(7\right)$

信号间常常是不相关的，则

是满秩的，且为对角矩阵。矢量

都不为零，i＝1，2，…M。为了减小噪声的贡献，执行降噪运算，取下式

${{\hat{R}}^{\′}}_X={\hat{R}}_X-{\mathrm{\γ}}^{2}I=[A{\hat{R}}_S{b}_1^H,A{\hat{R}}_S{b}_2^H,...A{\hat{R}}_S{b}_M^H]+{\mathrm{\δ}}^{2}I---\left(8\right)$

其中γ²为阵列信号协方差矩阵估计最小的两个特征值的平均值，即 ${\mathrm{\γ}}^2=({{\mathrm{\σ}}^2}_{M-1}+{\mathrm{\σ}}_M^2)/2.$ 当噪声为空间白噪声时γ²＝σ²，δ²＝0。实际情况下，噪声与空间白噪声有差异，此时δ²不为零，是一个小值。

对于本领域的普通技术人员来说，上述实施例中给出了计算协方差矩阵和执行降噪运算的一种方式，但是，本领域的普通技术人员根据本发明的教导，可以设计和实现其他方式的协方差矩阵计算和降噪运算而不脱离本发明的精神。

阵列信号协方差矩阵表示为

${\hat{R}}_X=U_C\mathrm{\Σ}{U}_C^H=U_S{\mathrm{\Σ}}_S{U}_S^H+U_N{\mathrm{\Σ}}_N{U}_N^H---\left(9\right)$

其中∑，∑_S和∑_N分别是阵列信号、信号和噪声协方差矩阵特征值组成的对角矩阵，U_C、U_S和U_N分别是阵列信号、信号和噪声协方差矩阵的特征矢量矩阵。可以证明，经过降噪处理的阵列信号协方差矩阵

为

${{\hat{R}}^{\′}}_X=U_C{\mathrm{\Σ}}^{\′}{U}_C^H=U_S{{\mathrm{\Σ}}^{\′}}_S{U}_S^H+U_N{{\mathrm{\Σ}}^{\′}}_N{U}_N^H---\left(10\right)$

其中∑′、∑′_S和∑′_N分别为阵列信号、信号和噪声的特征值减去噪声功率估计γ²后组成的对角阵，U_C、U_S和U_N与降噪前相同。U_C的表示为U_C＝[q₁，q₂，…q_M]，其中q_i是

的特征矢量，i＝1，2，…M。

的特征值为λ_i，i＝1，2，…M，将它从大到小排列。q_i按与其对应的特征值大小排列。求

在U_C中的投影

其i列为

P_i＝[p_i1，p_i2，…p_im…p_iM]^T    i＝1，2，…M    (11)

则

$P_i=[U_C^{H}A{\hat{R}}_S{b}_i^H+{\mathrm{\δ}}^2],i=\mathrm{1,2},...M---\left(12\right)$

由于源信号特征子空间与导向矢量在同一空间内，对应的投影p_im的模应该是一个相对较大的值。由于噪声的特征子空间与导向矢量正交，对应的投影p_im的模应该是一个相对较小的值。这是本发明的理论基础。下面介绍具体方法。

1.2信号功率归一化

各个信号的功率常常是不等的，为了提高信源数估计的性能，要对信号的功率进行归一化。

还可以表示为下式

$U_C^H{{\hat{R}}^{\′}}_X={\left[\begin{array}{cc}{U}_S& U_N\end{array}\right]}^H[U_S{{\mathrm{\Σ}}^{\′}}_S{U}_S^H+U_N{{\mathrm{\Σ}}^{\′}}_N{U}_N^H]=\left[\begin{array}{cc}{{\mathrm{\Σ}}^{\′}}_S& U_S^H\\ {{\mathrm{\Σ}}^{\′}}_N& U_N^H\end{array}\right]---\left(13\right)$

容易证明，∑′_SU_S ^H和∑′_NU_N ^H的各行分别正比于λ₁-γ²，λ₂-γ²，……λ_M-γ²。所以分别用λ₁-γ²，λ₂-γ²，……λ_M-γ²归一化，可以消除各个信号功率不同对信源数估计产生的不利影响。

由(13)式可知，功率归一化后，在U_C中的投影就是U_C ^H，因此我们可以用U_C＝[q₁，q₂，…q_M]估计目标数。q_i表示为下式

q_i＝[q_1i，q_2i，…q_mi，…q_Mi]^T        i＝1，2，…M    (14)

U_C ^H第m行元素模的均值为

$\stackrel{\‾}{\left|q_m\right|}=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\left|q_{\mathrm{im}}\right|,m=\mathrm{1,2}......M---\left(15\right)$

|q_m|称为信源数估计的判据。基于上面的分析可见，|q_m|的大小分为两组，|q_m|大的一组对应于信源，此组中m的个数就是信源数。

1.3信源数估计方法

基于上述研究和大量的数值分析，提出下面的信源数估计方法，它在模拟计算和实际应用中均证明是有效的。方法简述如下：

(1)由求U_C。

(2)按(15)式求|q_m|，m＝1，2，……M。

(3)求比例值β_m＝|q_m|²/|q_m+1|，m＝1，2，……M。

(4)β_m最大值对应的m即为信源数。

二、实数空间中特征空间的信源数估计方法

本节与上一节一样假设信号间互不相关，信号与噪声不相关，噪声为空间白噪声。定义如下稀疏矩阵Q_M，当M分别为2n和2n+1时

$Q_{2n}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_n& j{I}_n\\ J_n& -{\mathrm{jJ}}_n\end{array}\right]$ $Q_{2n+1}=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_n& 0& {\mathrm{jI}}_n\\ 0& \sqrt{2}& 0\\ J_n& 0& {-\mathrm{jJ}}_n\end{array}\right]---\left(16\right)$

其中J_n为n×n维置换矩阵，I_n为n×n维单位矩阵。

采用王永宁所著“空间谱估计理论与算法”一书中的方法，在实数空间中经过实数化处理后的阵列协方差矩阵可以表示为

${\hat{R}}_T=Q_M^H{\hat{R}}_{\mathrm{fb}}{Q}_M---\left(17\right)$

其中

为对整个阵列的接收数据作一次前后向空间平滑后得到的协方差矩阵。可以证明

的特征向量组成的矩阵为 $U_R=Q_M^H\left[\begin{array}{cc}{U}_S& U_N\end{array}\right],$ 又有 $Q_M^H{Q}_M=I,$ 则得下式

${\hat{R}}_T=Q_M^{H}A{\hat{R}}_S{A}^H{Q}_M+{\mathrm{\σ}}^{2}I---\left(18\right)$

下面的叙述与上一节类似，本节作一简单介绍。执行降噪运算

${{\hat{R}}^{\′}}_T={\hat{R}}_T-{\mathrm{\η}}^{2}I=Q_M^{H}A{\hat{R}}_S{A}^H{Q}_M+{\mathrm{\Δ}}^{2}I---\left(19\right)$

其中η²和Δ²与γ²和δ²定义类似，只是η²和Δ²在实数空间中。求

的特征矢量矩阵为U_R＝[q₁，q₂，…q_M]，其中

q_i＝[q_i1，q_i2…q_im…q_iM]^T  i＝1，2，……M    (20)

为了简化符号，在这里讨论的实数空间中，仍用了q_i，符号在形式上与复数空间中(14)式相同。在实际计算中注意它们的区别即可。功率归一化后，

就是U_R ^H，因此就可以用U_R估计目标数。按(15)式求出信源数估计的判据|q_m|，|q_m|大的一组中m的个数就是信源数。实数特征空间法的运算量是复数特征空间法的四分之一，如果计上变为实数的运算量开销，则约为四分之三。

三、本发明的波达方向估计中基于特征空间的信源数估计方法

下面总结出本发明的波达方向估计中基于特征空间的信源数估计方法，包括如下步骤：

(1)获得均匀线阵所接收的信号X(t)，X(t)＝[x₁(t)，x₂(t)，…x_M(t)]^T为M×1的矢量，其中x₁，...x_M表示基元1到M接收的信号，t表示时刻；M为均匀线阵的基元数目；

(2)计算协方差矩阵的估计值为：

${\hat{R}}_X=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}X\left(t_l\right)X{\left(t_l\right)}^H$

其中L为快拍数；

方法1，特征空间信源数估计方法，X(t)常为复数，复数空间特征空间信源数估计方法常简称为特征空间信源数估计方法，包括如下步骤：

(1-3)将

进行特征值分解，得到的特征矢量矩阵U_C＝[q₁，q₂…q_M]，其中q_i是

的特征矢量，i＝1，2，…M，q_i是按与其对应的特征值从大到小排列；此处进行特征值分解的方法可以采用比如G.Xu，T.Kailath的文章“A fast algorithm forsignal subspace decomposition and its performance analysis，Proc.IEEE mt.Conf.Acouct.Speech.Sig.Proc.1991：3069～3072”中提出的子空间迭代法，或者李有明、保铮的文章“基于子空间迭代的快速高分辨技术，电子学报，1994.2”，提出的预处理子空间迭代法；

(1-4)按(15)式求判据|q_m|，m＝1，2，……M；

(1-5)求比例值β_m＝|q_m|²/|q_m+1|，m＝1，2，……M；

(1-6)求β_m最大值对应的m，m＝1，2，……M，它是信源数。

方法2，实数特征空间信源数估计方法

实数特征空间法执行下列步骤

(2-3)按照王永宁所著“空间谱估计理论和方法”一书的方法由

求得

(2-4)由(17)式求得

(2-5)将

进行特征值分解，得到

的特征矢量矩阵 $U_R=Q_M^H\left[\begin{array}{cc}{U}_S& U_N\end{array}\right]=[q_1,q_2...q_M],$ 其中q_i是

的特征矢量，i＝1，2…M，q_i是按与其对应的特征值从大到小排列；

(2-6)按(15)式求|q_m|，m＝1，2，……M；

(2-7)求比例值β_m＝|q_m|²/|q_m+1|，m＝1，2，……M；

(2-8)求β_m最大值对应的m，m＝1，2，……M，它是信源数。

四、.采用本发明的信源数估计方法的波达方向估计方法

设一均匀线阵由M个基元组成，用来接收信号源发出的信号，设接收到K个远场窄带信号(即信号源的个数为K)。设各个信号源发出的信号对于均匀线阵的入射角为(θ₁，θ₂，…θ_K)。采用本发明的信源数估计方法的波达方向估计方法步骤如下：

(1)线阵接收的阵列信号为X(t)，它由(1)式和(2)式表示，其中信源数K由本发明的信源数估计方法给出。

(2)求协方差矩阵估计值

它与(5)式相同，是一个M×M阶的矩阵。上面在进行信源数估计时已经求出

和它的特征矢量矩阵，将它直接供给波达方向估计方法使用。

(3)采用王永宁所著“空间谱估计理论与算法”一书中介绍的多种波达方向估计方法获得波的入射方向。

尽管在上述实施例中结合附图以及处理流程对于本发明的构思、精神进行了描述，但是对于本领域的普通技术人员来说，可以认识到也可以以实现上述功能或者效果的装置或者产品的形式实现本发明，而不脱离本发明的构思的范围。

下面将采用本发明的方法进行模拟计算并与现有技术的方法进行比较，以体现本发明的上述有益效果。下面模拟计算用的基元接收到的信号X(t)是根据公式(1)和公式(2)计算出来的，并且S(t)是窄带高斯信号。

例1.

本例进行了模拟计算，模拟计算采用的主要参数与H.T.Wu等人文章中的例子类似，H.T.Wu等人采用其文章中的例子是表示其盖氏圆盘估计法的优点，其后的计算结果将表明本发明明显优于现有技术。X(t)模拟计算的参数为：均匀线阵相邻基元间距d为半波长，基元数M＝6个，目标数K＝2个，入射角(θ₁，θ₂)分别为-10°和10°，信噪比|S|²/2σ_n ²为5db，快拍数为100，入射信号频率为150kHz。本发明的实数特征空间法、复数特征空间法与盖氏圆盘法的结果均示于图2，前两者的判据值用o表示，后者的判据值，也即盖氏圆盘半径，用×表示，黑点·为噪声的判据值；图2(a)中横坐标为实验次数，每一次实验产生三个判据值，纵坐标为本发明的实数特征空间法信源数估计判据值，即前面所述的|q_m|的大小；图2(b)中的横坐标为实验的次数，每次试验产生三个判据值。纵坐标为本发明的复数特征空间法信源数估计判据值，即前面所述的|q_m|的大小；图2(c)中横坐标为实验次数，每一次实验产生三个判据值，纵坐标为盖氏圆盘半径判据值。图2(a)和(b)中，判据值较大的一组是位于上方的两个数值，它是两个目标的估计判据值。本发明的方法的两个目标的估计判据值离噪声的判据很远，约为2；而如图(c)所示，盖氏圆盘估计法的一个信号源目标的判据值与噪声的判据值很接近，很易产生欠估计。由于本发明的信源数估计中，信号源判据值和噪声判决值差别很大，所以本领域技术人员很容易判断信源数，提高了信源数估计精度。

例2.

模拟计算信源等功率情况下，本发明的实数特征空间法和复数特征空间法信源数估计的成功概率随信噪比的变化，并与盖氏圆法进行比较。X(t)模拟计算的参数为：8个基元的均匀线阵，基元间距为半波长。三个独立的等功率的窄带信号分别以10°、20°和30°的入射角入射到声阵上，信号频率为150kHz，噪声为空间高斯白噪声。快拍数为20，信噪比由-5db变到20db。实数特征空间法、复数特征空间法和盖氏圆法分别独立估计信源数200次，用蒙特卡罗法得到它们各自的成功概率，见图3。由图可见，在低信噪比区域，实数特征空间法和复数特征空间法的估计成功概率明显好于盖氏圆法，以实数特征空间法最好。此外，盖氏圆法要人工协助选取估计判据的阈值，本发明的实数特征空间法和复数特征空间法是自动设定的。

例3.

模拟计算信源功率不等情况下，本发明的实数特征空间法和复数特征空间法信源数估计的成功概率随信噪比的变化，并与盖氏圆法进行比较。X(t)模拟计算的参数为：三个独立的不等功率的信号入射到声阵上，入射角分别为40°、45°和60°，信号功率比为1∶2∶1，信噪比以最小信号功率为基准，快拍数为50。其它条件同例2，计算结果见图4。由图可见，在信源功率不等的情况下，实数特征空间法和复数特征空间法依然好于盖氏圆法，以实数特征空间法最好。此外，盖氏圆法要人工协助选取估计判据的阈值，实数特征空间法和复数特征空间法是自动设定的。

例4.

模拟计算空间相关色噪声下，本发明的实数特征空间法和复数特征空间法信源数估计成功概率，并与盖氏圆法进行比较。空间相关色噪声的模型为

$R_{i,k}=E\left[n_i^{*}\left(t\right)n_k\left(t\right)\right]={\mathrm{\β}}^{|i-k|}\exp \left(j2\mathrm{\π\φ}(i-k)\right),i,k=\mathrm{1,2}......M---\left(20\right)$

其中β∈[0，1]是阵列基元上色噪声的空间相关系数，φ为0.1弧度，它控制噪声功率谱锋值的位置。X(t)模拟计算的参数为：3个独立的不等功率窄带信号分别以入射角20°、40°和60°入射到均匀线阵上，信号功率比1∶2∶1，信噪比为10db，以最小信号功率为基准。其它条件同例2，计算结果见图5，由图可见，在相关系数大时，实数特征空间法和复数特征空间法明显好于盖氏圆法，以实数特征空间法性能最好。

本发明是对M×M阶协方差矩阵进行特征值分解，它可以直接用于以后的波达方向估计中，无需再进行一次特征值分解。由于在波达方向估计中，特征值分解的运算量占总运算量至少一半，因此利用本发明能够为波达方向估计信号处理过程省去大量的计算量，减少硬件开销。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，本发明所述的技术方案既可以通过软件的手段实现，也可以通过硬件或软硬结合的方式实现，同时，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (12)

1.一种信源数估计方法，包括如下步骤：

(1)获得均匀线阵所接收的信号X(t)，X(t)＝[x₁(t)，x(t)，…x_M(t)]^T为M×1的矢量，其中x₁，...x_M表示所述均匀线阵中基元1到M接收的信号，t表示时刻；

(2)由接收信号计算协方差矩阵

(3)将

进行特征值分解，得到它的特征矢量矩阵U_C＝[q₁，q₂…q_M]，其中q_i是的特征矢量，q_i＝[q_1i，q_2i…q_mi…q_Mi]^T，i＝1，2，…M；

(4)根据q_i计算信源数m。

2、根据权利要求1的方法，其特征在于在步骤(2)中，利用公式 ${\hat{R}}_X=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}X\left(t_l\right)X{\left(t_l\right)}^H,$ 进行协方差矩阵计算，其中L为快拍数。

3、根据权利要求2的方法，其特征在于在步骤(2)之后还包括步骤：

(a)由协方差矩阵

计算矩阵

其中为对整个阵列的接收数据作一次前后向空间平滑后得到的协方差矩阵；

(b)由

求得实数化处理后的协方差矩阵 ${\hat{R}}_T=Q_M^H{\hat{R}}_{\mathrm{fb}}{Q}_M;$

以及在步骤(3)中包括：使 ${\hat{R}}_X={\hat{R}}_T,$ 并将

进行特征值分解，得到它的特征矢量矩阵 $U_R=Q_M^H{U}_C,$ U_R＝[q₁，q₂，…q_M]，其中q_i是

的特征矢量，q_i＝[q_1i，q_2i…q_mi…q_Mi]^T，i＝1，2，…M。

4、根据权利要求1-3之一的方法，其特征在于步骤(4)包括下列步骤：

(4.1)利用公式 $\stackrel{\‾}{\left|q_m\right|}=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\left|q_{\mathrm{im}}\right|$ 求判据值|q_m|，m＝1，2……M；

(4.2)利用判据值|q_m|确定信源数m。

5、根据权利要求4的方法，其特征在于步骤(4.2)还包括步骤：

(4.3)根据公式β_m＝|q_m|²/|q_m+1|求比例值β_m，其中m＝1，2……M；

(4.4)求β_m最大值对应的m，m＝1，2……M，m是信源数。

6、根据权利要求1的方法，其特征在于所述步骤(3)中q_i是按照其对应特征值大小排列的，i＝1，2，…M。

7、一种波达方向估计方法，包括步骤：

(1)均匀线阵接收阵列信号X(t)，由权利要求1-6之一的信源数估计方法计算出信源数K；

(2)取得权利要求1-6之一的协方差矩阵估计值它是一个M×M阶的矩阵；

(3)采用上述

进行波达方向估计，获得波的入射方向。

8、一种信源数估计装置，包括：

(1)均匀线阵，包括基元1到M，用于获得均匀线阵所接收的信号X(t)，X(t)＝[x₁(t)，x(t)，…x_M(t)]^T为M×1的矢量，其中x₁，...x_M表示所述均匀线阵中基元1到M接收的信号，t表示时刻；

(2)协方差矩阵计算单元，用于由接收信号计算协方差矩阵

(3)特征值分解单元，用于将进行特征值分解，得到它的特征矢量矩阵U_C＝[q₁，q₂…q_M]，其中q_i是

的特征矢量，q_i＝[q_1i，q_2i…q_mi…q_Mi]^T，i＝1，2，…M；

(4)信源数计算单元，用于根据q_i计算信源数m。

9、根据权利要求8的信源数估计装置，其特征在于，所述协方差矩阵计算单元利用公式 ${\hat{R}}_X=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}X\left(t_l\right)X{\left(t_l\right)}^H,$ 进行协方差矩阵计算，其中L为快拍数。

10、根据权利要求9的信源数估计装置，其特征在于还包括下列单元：

用于由协方差矩阵

计算矩阵的单元，其中

为对整个阵列的接收数据作一次前后向空间平滑后得到的协方差矩阵；

用于由

求得实数化处理后的协方差矩阵 ${\hat{R}}_T=Q_M^H{\hat{R}}_{\mathrm{fb}}{Q}_M$ 的单元；

以及所述特征值分解单元还进行下列计算：使 ${\hat{R}}_X={\hat{R}}_T,$ 并将

进行特征值分解，得到它的特征矢量矩阵 $U_R=Q_M^H{U}_c,$ U_R＝[q₁，q₂，…q_M]，其中q_i是

的特征矢量，q_i＝[q_1i，q_2i…q_mi…q_Mi]^T，i＝1，2，…M。

11、根据权利要求8-10之一的信源数估计装置，其特征在于信源数计算单元还包括：

利用公式 $\stackrel{\‾}{\left|q_m\right|}=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\left|q_{\mathrm{im}}\right|$ 求判据值|q_m|的单元，m＝1，2……M；以及利用判据值|q_m|确定信源数m的单元。

12、根据权利要求11的信源数估计装置，其特征在于利用判据值|q_m|确定信源数m的单元还包括：

根据公式β_m＝|q_m|²/|q_m+1|求比例值β_m的单元，其中m＝1，2……M；求β_m最大值对应的m的单元，m＝1，2……M，m是信源数。

Images