波达方向(direction-of-arrival, DOA)^[1-3]估计作为阵列信号处理的一个研究重点，该技术一直被广泛应用于雷达、声呐探测等方面。传统的高分辨，例如最大似然估计(maximum likelihood, ML)^[4]、旋转不变子空间算法(estimating signal parameter via rotational invariance techniques, ESPRIT)^[5]、多重子空间分类算法(multiple signal classification, MUSIC)^[6]以及加权子空间类算法(weighted subspace fitting, WSF)^[7]等，均需要将信源个数作为先验知识，才能进行DOA估计^[8]。然而，在实际工程中，不能直接获取目标个数，并且对目标数目进行准确估计在工程上也具有一定难度和挑战。现有的对于目标数进行估计的算法，包括最小长度描述法(minimum description length, MDL)^[9-10]以及信息论准则(Akaike information criterion, AIC)^[11]。这些经典算法最大的弊端在于当环境噪声不满足高斯白噪声时，很难获得准确的信源数目。虽然，有一些算法针对此情况进行了改进，然而，在信噪比(signal-to-noise ratio, SNR)较低、快拍数不足时，依然难以获得准确的信源数目。

针对上述问题，有学者提出了无需信源数目的DOA估计方法，例如众所传统波束形成(conventional beamforming, CBF)以及最小方差无畸变响应(minimum variance distortionless response, MVDR)^[12]算法，此类算法虽然无需已知信源个数，但是却不能获得超分辨。为此，Qian等^[13]提出了一种无需信源个数的DOA算法，此算法主要利用了Toeplitz矩阵重构，并利用特定矩阵的最大特征值进行DOA估计，实验证明该算法可以获得较好的估计性能，但是该算法需要进行Toeplitz矩阵重构，因而使得算法所能探测的目标数目减少至阵元数的一半。文献[14]提出一种利用神经网络的来进行DOA估计的方法，虽然，此方法也无需已知信源个数，并且不会损失阵列的孔径，但是为了对神经网络进行训练，此方法需要提供大量的有效数据进行学习。

现有的DOA算法，大多从特征向量入手进行分析，而忽略了对于特征值的运用。虽然文献[15]提出了一种NSP(noise subspace to power)算法，该算法从特征值入手，利用噪声功率的不变性，构建了空间谱，从而使得算法性能优于MUSIC算法，但是此算法依旧没有摆脱子空间类算法的弊端，需要已知信源个数。为此，本文在NSP算法的框架上进行分析，提出了一种特殊的特征值排布方式，限定了其最小特征值对应的特征向量并利用其在不同角度上的正交特性，构建空间谱，从而使得本文算法无需已知信源个数，且拥有着超分辨效果。

1 基于特征值排序的DOA估计算法 1.1 数据模型

假设K个远场且相互独立的窄带信号入射到一个由M个阵元组成的均匀线阵中，均匀线阵的阵间距d为半波长。因此，阵列接收的信号可以表示为：

$ \mathit{\boldsymbol{X}}(t) = \mathit{\boldsymbol{A}}(\mathit{\boldsymbol{ \boldsymbol{\varTheta} }})\mathit{\boldsymbol{S}}(t) + \mathit{\boldsymbol{N}}(t) $

(1)

式中：S(t)∈C^K×1表示信号；N(t)∈C^M×1表示传感器接收到的加性噪声；A(Θ)=[a(θ₁), a(θ₂), …, a(θ_K)]∈C^M×K代表阵列导向矢量，Θ=[θ₁, θ₂, …, θ_K]代表入射角度的集合，对于角度θ的导向矢量可以表示为：

$ \mathit{\boldsymbol{a}}(\theta ) = {[1,{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}d\sin (\theta )/\lambda }}, \cdots ,{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}(M - 1)d\sin (\theta )/\lambda }}]^{\rm{T}}} $

(2)

式中λ代表信号的波长，因此接收数据的协方差矩阵可以表示为：

$ \begin{array}{*{20}{l}} {R = {\rm{E}}[\mathit{\boldsymbol{X}}(t){\mathit{\boldsymbol{X}}^{\rm{H}}}(t)] = \mathit{\boldsymbol{A}}(\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}){\mathit{\boldsymbol{R}}_S}{\mathit{\boldsymbol{A}}^{\rm{H}}}(\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{E}}[\mathit{\boldsymbol{N}}(t){\mathit{\boldsymbol{N}}^{\rm{H}}}(t)]} \end{array} $

(3)

其中R_S=E[S(t)S^H(t)]代表信号协方差矩阵。

1.2 算法思路及空间谱构造

MUSIC算法主要利用噪声子空间与导向矢量的正交性进行测向，即：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{U}}_N^{\rm{H}}{\mathit{\boldsymbol{U}}_N}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ) \to 0,}&{\forall \theta \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}\\ {\mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{U}}_N^{\rm{H}}{\mathit{\boldsymbol{U}}_N}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ) \ne 0,}&{\forall \theta \notin \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \end{array}} \right. $

(4)

式中U_N代表由式(3)分解得到的噪声子空间。由式(4)可知，导向矢量必与噪声子空间中的每一列正交，因此，可以得到：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{e}}_M^{\rm{H}}{\mathit{\boldsymbol{e}}_M}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ) \to 0,}&{\forall \theta \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}\\ {\mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{e}}_M^{\rm{H}}{\mathit{\boldsymbol{e}}_M}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ) \ne 0,}&{\forall \theta \notin \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \end{array}} \right. $

(5)

式中e_M代表第M个特征值对应的特征向量。利用式(5)构建空间谱，可得：

$ f(\theta ) = \frac{1}{{\mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{e}}_M^{\rm{H}}{\mathit{\boldsymbol{e}}_M}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta )}} $

(6)

利用式(6)进行仿真，设置SNR为0 dB，快拍数为500，阵元数为6，入射角度为0°和10°进行空间谱仿真实验，其结果如图 1所示，从图 1中可以看到，算法虽然在0°和10°处形成了谱峰，但是在非目标处也形成谱峰，产生了虚假目标，从而导致式(6)算法失效。

为此，可以按照文献[15]的方法构建一个协方差矩阵：

$ \mathit{\boldsymbol{\bar R}} = \mathit{\boldsymbol{R}} + \beta \mathit{\boldsymbol{a}}(\theta ){\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ) $

(7)

式中参数β为一个待确定的变量。此时，希望R的最小特征值对应的特征向量e_M满足：

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\bar e}}}_M} = {m_1}{\mathit{\boldsymbol{e}}_M},}&{\forall \theta \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}\\ {{{\mathit{\boldsymbol{\bar e}}}_M} = m{\mathit{\boldsymbol{e}}_M} + n\mathit{\boldsymbol{e}},}&{\forall \theta \notin \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \end{array}} \right. $

(8)

式中：e代表a(θ)a^H(θ)的特征向量；m₁、m、n为相应的系数。将式(8)获得的e_M代入式(6)中，可得：

$ f(\theta ) = \frac{1}{{\mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{\bar e}}_M^{\rm{H}}{{\mathit{\boldsymbol{\bar e}}}_M}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta )}} $

(9)

因此，根据正交性可得：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{\bar e}}_M^{\rm{H}}{{\mathit{\boldsymbol{\bar e}}}_M}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ) \to 0,}&{\forall \theta \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}\\ {\mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{\bar e}}_M^{\rm{H}}{{\mathit{\boldsymbol{\bar e}}}_M}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ) \ne 0,}&{\forall \theta \notin \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \end{array}} \right. $

(10)

并且存在：

$ \frac{{\{ \mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{\bar e}}_M^{\rm{H}}{{\mathit{\boldsymbol{\bar e}}}_M}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta )|\forall \theta \notin \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\} }}{{\{ \mathit{\boldsymbol{a}}(\theta )\mathit{\boldsymbol{\bar e}}_M^{\rm{H}}{{\mathit{\boldsymbol{\bar e}}}_M}{\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta )|\forall \theta \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\} }} \to \infty $

(11)

由式(11)可以看出，目标处与非目标处数值相差较大，因此式(9)将不会在非目标处形成谱峰，从而克服了式(6)的弊端。

1.3 β边界条件分析

为后续讨论方便，假设N(t)为高斯白噪声，从而可得E[N(t)N^H(t)]=σ_n²I，σ_n²为噪声的功率。为了满足式(8)，因此R的特征向量必须满足：

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \lambda (\mathit{\boldsymbol{\bar R}}|\theta = {\theta _k},{\theta _k} \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}) = \\ {\rm{diag}} (\underbrace {{\lambda _1}, \cdots ,{\lambda _k} + \varDelta , \cdots ,{\lambda _K}}_K,\underbrace {\sigma _n^2, \cdots ,\sigma _n^2}_{M - K}) \end{array} $

(12)

其中，Δ=βM，令S={1, 2，…, K}，{λ_i|i∈S}代表R的特征值。式(8)另外一种情况下的特征值表示为：

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \lambda (\mathit{\boldsymbol{\bar R}}|\theta = {\theta _k},{\theta _k} \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}) = \\ {\rm{diag}} (\underbrace {{\lambda _1}, \cdots ,{\lambda _K}}_K,\underbrace {\sigma _n^2, \cdots ,\varDelta + \sigma _n^2}_{M - K}) \end{array} $

(13)

对于式(12)、式(13)的特征值排布模式进行推导：

$ \mathit{\boldsymbol{\bar R}}(\theta = {\theta _k},{\theta _k} \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}) = \sum\limits_{i = 1}^M {{\lambda _i}} {\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{e}}_i^{\rm{H}} + \beta M\mathit{\boldsymbol{e}}{\mathit{\boldsymbol{e}}^{\rm{H}}} $

(14)

当θ=θ_k时，必存在ee^H=e_ke_k^H，因此式(14)变换为：

$ \begin{array}{l} \mathit{\boldsymbol{\overline R}} \left( {\theta = {\theta _k},{\theta _k} \in \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \right) = \sum\limits_{i = 1,i \ne k}^M {{\lambda _i}} {\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{e}}_i^{\rm{H}} + \left( {\beta M + {\lambda _k}} \right){\mathit{\boldsymbol{e}}_k}\mathit{\boldsymbol{e}}_k^{\rm{H}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{i = 1,i \ne k}^K {{\lambda _i}} {\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{e}}_i^{\rm{H}} + \left( {\beta M + {\lambda _k}} \right){\mathit{\boldsymbol{e}}_{\bf{k}}}e_k^{\rm{H}} + \sum\limits_{i = K + 1}^M {\sigma _n^2} {\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{e}}_i^H \end{array} $

(15)

因此只要选取合适的β值便可满足式(12)的特征值排布情况，同理也可得到式(13)的特征值排布情况。式(12)说明，最小特征值对应的特征向量是噪声所形成的；而式(13)说明，最小特征值对应的特征向量是由信号和噪声组合而成，满足式(8)。为此文章讨论β的取值范围，根据式(12)和式(13)中特征值的排布情况，可以得到关于β的限制条件：

$ \left\{ {\begin{array}{*{20}{l}} {\beta M + \sigma _n^2 < \sigma _n^2}\\ {\sigma _n^2 \le {\lambda _k} + \beta M \le {\lambda _1}} \end{array}} \right. $

(16)

为了能够准确获得β取值范围，首先必须给出信号功率与信号特征值之间的关联，将式(3)进行改写可得：

$ \mathit{\boldsymbol{R}} = \sum\limits_{i = 1}^K {\sigma _i^2} \mathit{\boldsymbol{a}}({\theta _i}){\mathit{\boldsymbol{a}}^{\rm{H}}}({\theta _i}) + \sigma _n^2\mathit{\boldsymbol{I}} $

(17)

式中：{σ_i²|i∈S}代表第i个信号的功率。根据rank(a(θ)a^H(θ))=1^[16]，可知：

$ \mathit{\boldsymbol{a}}(\theta ){\mathit{\boldsymbol{a}}^{\rm{H}}}(\theta ) = M\mathit{\boldsymbol{e}}{\mathit{\boldsymbol{e}}^{\rm{H}}} $

(18)

由于rank(a(θ)a^H(θ))=1，所以a(θ)a^H(θ)特征值可以表示为tr(a(θ)a^H(θ))=tr(a^H(θ)a(θ))=M。根据式(18)，式(17)可以改写为：

$ \mathit{\boldsymbol{R}} = \sum\limits_{i = 1}^K M \sigma _i^2{\mathit{\boldsymbol{e}}_i}\mathit{\boldsymbol{e}}_i^{\rm{H}} + \sigma _n^2\mathit{\boldsymbol{I}} $

(19)

根据文献[17]中特征值的运算性质：

$ { {\rm{eig}} (\mathit{\boldsymbol{A}} - c\mathit{\boldsymbol{I}}) = {\lambda _{\mathit{\boldsymbol{A}},i}} - c} $

$ { {\rm{eig}} (\mathit{\boldsymbol{I}} + c\mathit{\boldsymbol{A}}) = c{\lambda _{\mathit{\boldsymbol{A}},i}} + 1} $

式中：A为对称矩阵；c为常数；λ_{A, i}为A的特征值。因此理想情况下信号特征值与信号功率之间的关系可以表示为：

$ \{ {\lambda _i} - \sigma _n^2 = M\sigma _i^2|i \in S\} $

(20)

根据式(20)和式(16)可以得到β的边界条件：

$ - \sigma _k^2 \le \beta < 0 $

(21)

根据SNR的定义，并将β进行适量的缩放，选取边界内的一个合适的β，其选取方式可以表示为：

$ \beta = - {\lambda _{\min ,R}}\frac{{{{10}^{{\lambda _{\max ,R}}}}/10M}}{M} $

(22)

2 计算机仿真实验及结果分析 2.1 算法可行性分析

假设2个等功率且不相关的远场窄带信号，以0°和10°的波达角入射到一个由6个声压传感器组成的均匀线阵上，阵元间距d=0.5λ，快拍数为300，设置SNR为0 dB。进行MVDR、NSP以及本文算法的对比实验。其结果如图 2所示，从图 2中，可以看到当2个目标的角度差较小时，MVDR算法无法进行准确识别，而本文算法以及NSP算法却能准确识别2个目标，从而说明本文算法与NSP算法拥有着高于MVDR算法的分辨率。而且此时本文算法，旁瓣比较平坦，没有形成错误的谱峰，说明式(11)的正确性，进一步说明所提算法理论的正确性可行性。

2.2 算法统计性能分析

为了进一步验证算法的性能，进行200次蒙特卡洛独立实验。保持其他条件不变，设置SNR从-5 dB开始，每次增加2 dB，一直到15 dB，进行3种算法的均方根误差(root mean square error，RMSE)对比分析实验，其结果如图 3所示，从图 3中可以看到随着信噪比的增加，本文算法以及NSP算法的均方根误差随之减小，但是MVDR算法在SNR小于7 dB时，其均方根误差减小较为缓慢，而在SNR大于7 dB时，算法均方根误差急剧减小，MVDR算法在小于7 dB时，无法区分0°和10°的2个目标，从而说明本文算法以及NSP算法在低信噪比时分辨率优于MVDR算法，所提算法拥有着超分辨效果。

保持其他条件不变，设置一个目标的波达角为0°，另外一个角度为(0+ϕ)°，ϕ可表示为一个由5开始每次递增2，直到25为止的等差数列。因此，ϕ可以看作2个目标之间的角度差；设置SNR为0 dB，进行200次蒙特卡洛实验，得到不同角度差下算法的均方根误差，其结果如图 4所示，从图 4中，可以看到，当目标角度差小于10°时，本文算法性能是最好的，说明本文算法更加适用于小角度差的环境下，随着角度差的增加，3种算法的均方根误差趋近于一致，从整体上看本文算法的性能优于MVDR算法，说明本文算法可以在无需已知信源个数的前提下获得超分辨效果。

为分析快拍数对于算法的影响，保持其他条件不变，设置目标入射角度为0°和15°，快拍数为100开始，每次递增200，直到900为止，进行200次蒙特卡洛实验，其均方根误差如图 5所示。从图 5中，可以看到随着快拍数的增加，MVDR算法的均方根误差也随之减小，但此时其性能在3种算法中是最差的，而NSP算法是受快拍数影响较小的算法，本文算法在快拍数小于300时，受快拍影响较大，当快拍数大于300时，本文算法的均方根误差较为稳定且与NSP算法较为接近。

3 结论

1) 本文通过推导信号特征值与信号功率之间的关联，给出了β的边界条件以及取值方法，从而实现了本文算法，并克服了噪声最小特征值对应的特征向量与扫描导向矢量正交时，扫描谱上出现虚假谱峰的弊端。

2) 本文给出一种无需已知信源数目的DOA算法，其性能优于MVDR算法，且在小角度差时，其分辨率高于NSP算法。

3) 本文算法虽然拥有较好的估计精度以及分辨率，但算法需要在每个搜索角度进行特征值分解，以获得所需特征向量，计算量较大。

本文的给出了信号特征值与信号功率之间的数学表达式，揭示了两者之间的内在关联。但本文算法的重点在于β的选择，如何利用β提高算法性能、减小信号检测门限以及减小算法计算量，是后期的研究重点。