0 概述

在阵列信号处理过程中，空间谱表示信号在空间各个方向的能量分布^[1]，常用的空间谱估计算法有旋转不变子空间（Estimating Signal Parameters via Rotational Invariance Techniques，ESPRIT）算法^[2]、多重信号分类（Multiple Signal Classification，MUSIC）算法^[3]和相关改进算法^[4-5]，这些算法都是以事先已知信源的个数为前提条件。目前，信源个数估计多使用基于信息论准则的方法，包括最小信息准则（Akaike Information Criterion，AIC）^[6]、最小描述长度（Minimum Description Length，MDL）准则^[7]及相关改进算法^[8-9]，但这些算法仅适用于白噪声环境下的信源个数估计。

适用于色噪声环境的盖尔圆估计（Gerschgorin Disk Estimation，GDE）准则^[10]虽然能够弥补信息论准则在色噪声环境下信源个数估计失效的不足，但其在低信噪比和小快拍数条件下的性能会急剧下降甚至失效。针对该问题，文献[11]提出一种新的酉变换方法，并结合GDE准则进行信源个数估计。实验结果表明，该方法在快拍数为90时的信源个数检测准确率达到90%，但是要求信噪比在10 dB以上。文献[12]将GDE准则和AIC准则结合，克服了特征值的无序性导致估计错误的缺陷，使得检测准确率在信噪比为-5 dB时达到93%，但该方法要求较大的快拍数。文献[13]对盖尔圆半径进行压缩，并结合盖尔圆心值提出一种基于自适应调整因子的GDE准则。实验结果表明，该方法在信噪比为-4 dB、快拍数大于2 000的条件下检测准确率为80%，而在信噪比为15 dB、快拍数为50时检测准确率达到90%。由此可知，该方法无法同时保证低信噪比和小快拍数条件下的检测性能。文献[14]利用接收信号协方差矩阵的对角线平均值来构建新的协方差矩阵，并结合GDE准则估计信号源数。该算法在非平稳的色噪声环境中，信噪比为-8 dB、快拍数为100时的检测准确率就已经达到80%以上，虽然能同时保证低信噪比和小快拍数情况下的高检测准确率，但仅限于在10个阵元估计2个信源的条件下。文献[15]基于Khatri-Rao积先对接收信号的协方差矩阵做延迟处理后再对矩阵进行矢量化来构造新的矩阵，并结合盖尔圆准则估计信源数。该方法虽然能在$ M $个阵元下至多估计出$ (2M-1) $个信号源，但仅限于在白噪声条件下，具有一定的局限性。

上述GDE准则及改进方法对于含$ M $个阵元的阵列，在构造盖尔圆盘过程中都只用到了前（$ M-1 $）个阵元接收的信息，导致这些方法能估计的最大信源个数为$ M-2 $。文献[16]提出了加权盖尔圆估计（Weighted Gerschgorin Disk Estimation，WGDE）准则，由于该准则中的特征加权矩阵对增广盖尔圆矩阵中的盖尔圆半径做特征加权变换，因此进一步增大了信号盖尔圆和噪声盖尔圆之间半径的差异，使快拍数为64时的信源估计准确率达到90%，但是信噪比要求为13 dB。同时，WGDE准则在构造增广盖尔圆矩阵及其增广加权盖尔圆矩阵过程中都保留了$ M $个阵元的信息，这使得估计$ (M-1) $个信号源数成为可能。文献[17]将信源数估计问题转为模式识别问题，依据阵元在接收含有不同信源个数的入射信号时会产生不同相位差异这一特点，结合希尔伯特-黄变换（Hilbert-Huang Transformation，HHT）和SVM进行信源个数估计，从而利用$ M $个阵元估计（$ M-1 $）个信号源数，但该文并未对快拍数小于100时的信源数估计性能进行验证。

本文结合文献[17]方法，针对WGDE准则对增广加权盖尔圆矩阵信息利用不足的缺陷，提出一种在增广加权盖尔圆矩阵中获取多重特征并融合的信源个数估计方法。同时获取可用于描述信源个数的盖尔圆心值、盖尔圆半径和加权盖尔圆半径等多重特征进行融合，构建能够描述信源个数的高维特征向量，标定后代入支持向量机（Support Vector Machine，SVM）中训练可用于信源个数估计的分类器，并利用包含$ M $个天线的均匀圆阵（Uniform Circular Array，UCA）（下文简称M-UCA）进行仿真实验。

1 M-UCA的接收信号模型

假设$ k $个远场窄带信号入射到M-UCA上，其中$ k=\mathrm{1, 2}, \cdots , K $且$ K<M $，则M-UCA接收信号表示为：

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

(1)

其中，$ \mathit{\boldsymbol X}\left(t\right)=\left[{x}_{1}\right(t)\mathrm{ }, {x}_{2}(t)\mathrm{ }, \cdots , {x}_{M}{\left(t\right)]}^{\mathrm{T}} $为$ t $时刻的接收信号矢量，$ \mathit{\boldsymbol S}\left(t\right)=\left[{s}_{1}\right(t)\mathrm{ }, {s}_{2}(t)\mathrm{ }, \cdots , {s}_{k}{\left(t\right)]}^{\mathrm{T}} $为$ t $时刻信号的源矢量，$ \mathit{\boldsymbol N}\left(t\right)=\left[{n}_{1}\right(t)\mathrm{ }, {n}_{2}(t)\mathrm{ }, \cdots , {n}_{M}{\left(t\right)]}^{\mathrm{T}} $为$ t $时刻加性噪声矢量，$ \mathit{\boldsymbol A}=\left[\alpha \left({\mathit \Theta }_{1}\right)\mathrm{ }, \alpha \left({\mathit \Theta }_{2}\right)\mathrm{ }, \cdots , \alpha \left({\mathit \Theta }_{K}\right)\right] $为M-UCA的阵列流型矢量，$ {\mathit \Theta }_{k}=({\varphi }_{k}, {\vartheta }_{k}) $为第$ k $个信号源的入射角，$ {\varphi }_{k} $和$ {\vartheta }_{k} $分别为第$ k $个信号源的方位角和俯仰角，且$ {\varphi }_{k}\in \left(\mathrm{0, 2}\mathrm{\pi }\right) $，$ {\vartheta }_{k}\in \left(0, \frac{\mathrm{\pi }}{2}\right) $，$ \mathit{\boldsymbol \alpha }\left({\mathit \Theta }_{k}\right)=\left[{\mathrm{e}}^{\mathrm{j}\eta \mathrm{c}\mathrm{o}\mathrm{s}({\varphi }_{\mathit{k}}-{\beta }_{1})}, {\mathrm{e}}^{\mathrm{j}\eta \mathrm{c}\mathrm{o}\mathrm{s}({\varphi }_{\mathit{k}}-{\beta }_{2})}\right.{\left., \cdots , {\mathrm{e}}^{\mathrm{j}\eta \mathrm{c}\mathrm{o}\mathrm{s}({\varphi }_{\mathit{k}}-{\beta }_{\mathit{m}})}, \cdots , {\mathrm{e}}^{\mathrm{j}\eta \mathrm{c}\mathrm{o}\mathrm{s}({\varphi }_{\mathit{k}}-{\beta }_{\mathit{M}})}\right]}^{\mathrm{T}} $为第$ k $个信号源的导向矢量，其中，$ {\beta }_{m}=\frac{2\mathrm{\pi }}{M}\cdot (m-1) $为逆时针第$ m $个阵元与$ x $轴的夹角，$ \eta =\frac{2\mathrm{\pi }r}{\widehat{\lambda }}\mathrm{s}\mathrm{i}\mathrm{n}{\vartheta }_{k} $，$ m=\mathrm{1, 2}, \cdots , M $，$ \widehat{\lambda } $为载波波长。

2 增广加权盖尔圆矩阵与加权盖尔圆准则

在满足奈奎斯特采样定理^[18]的条件下，对$ \mathit{\boldsymbol X}\left(t\right) $做快拍数为$ L $的均匀采样，获得$ \mathit{\boldsymbol X}\left(t\right) $的观测信号数据集$ \mathit{\boldsymbol X}\left(l\right) $，$ l=\mathrm{1, 2}, \cdots , L $且$ \mathit{\boldsymbol X}\left(l\right)=\left[{x}_{1}\left(1\right)\mathrm{ }, {x}_{2}\left(2\right)\mathrm{ }, \cdots , \right.{\left.{x}_{M}\left(l\right)\right]}^{\mathrm{T}} $，则观测信号数据集的协方差矩阵为：

$ \boldsymbol{R}_{\mathrm{xx}}=\frac{1}{L} \sum\limits_{l=1}^{L} \boldsymbol{X}(l) \boldsymbol{X}(l)^{\mathrm{H}}$

(2)

其中，$ {(\mathrm{ }\cdot \mathrm{ })}^{\mathrm{H}} $表示共轭转置。对$ {\mathit{\mathit{\boldsymbol R}}}_{\mathrm{x}\mathrm{x}} $做特征分解得到：

$ {\mathit{\boldsymbol R}}_{\mathrm{x}\mathrm{x}}=\sum\limits_{m=1}^{M}{\lambda }_{m}{\mathit{\boldsymbol u}}_{m}{\mathit{\boldsymbol u}}_{m}^{\mathrm{H}} $

(3)

其中，$ {\lambda }_{m} $为$ {\mathit{\boldsymbol R}}_{\mathrm{x}\mathrm{x}} $的第$ m $个特征值，$ {\lambda }_{1}>{\lambda }_{2}>, \cdots , >{\lambda }_{k-1}>{\lambda }_{k}>>, \cdots , >{\lambda }_{M} $，当$ m\le k $时，$ {\lambda }_{m} $为信号特征值，反之为噪声特征值，$ {\mathit{\boldsymbol u}}_{m} $为$ {\lambda }_{m} $所对应的特征向量。定义如式（4）所示的矩阵：

$ {\mathit{\boldsymbol R}}_{M+1}=\left[\begin{array}{cc}{\mathit{\boldsymbol R}}_{\mathrm{x}\mathrm{x}}& \mathit{\boldsymbol r}\\ {\mathit{\boldsymbol r}}^{\mathrm{H}}& \boldsymbol 1\end{array}\right] $

(4)

其中，$ \mathit{\boldsymbol r} $为$ {\mathit{\boldsymbol R}}_{\mathrm{x}\mathrm{x}} $的最后一列向量。定义酉变换矩阵T如式（5）所示：

$ \mathit{\boldsymbol T}=\left[\begin{array}{cc}\mathit{\boldsymbol u}& \boldsymbol 0\\ \boldsymbol 0& \boldsymbol 1\end{array}\right] $

(5)

其中，$ \mathit{\boldsymbol u}={\left[{u}_{1}, {u}_{2}, \cdots , {u}_{M}\right]}^{\mathrm{T}} $。经过如式（6）所示的酉变换计算得到增广盖尔圆矩阵$ {\mathit{\boldsymbol R}}_{\mathrm{G}} $：

$ {\mathit{\boldsymbol R}}_{\mathrm{G}}={\mathit{\boldsymbol T}}^{\mathrm{H}}{\mathit{\boldsymbol R}}_{M+1}\mathit{\boldsymbol T}=\left[\begin{array}{cc}\mathit{\pmb{\Sigma }}& {\mathit{\boldsymbol u}}^{\mathrm{H}}\mathit{\boldsymbol r}\\ {\mathit{\boldsymbol r}}^{\mathrm{H}}\mathit{\boldsymbol u}& \boldsymbol 1\end{array}\right]= \\ \;\;\;\; \left[\begin{array}{ccccc}{\lambda }_{1}& 0& \cdots & 0& {p}_{1}\\ 0& {\lambda }_{2}& \cdots & 0& {p}_{2}\\ 0& 0& & ⋮& ⋮\\ ⋮& ⋮& & {\lambda }_{M}& {p}_{M}\\ {p}_{1}^{\mathrm{*}}& {p}_{2}^{\mathrm{*}}& \cdots & {p}_{M}^{\mathrm{*}}& 1\end{array}\right] $

(6)

其中，$ \mathit{\pmb{\Sigma }}$为所有$ {\lambda }_{m} $组成的对角阵。定义特征加权矩阵W如式（7）所示：

$ \mathit{\boldsymbol W}=\left[\begin{array}{cc}\mathit{\pmb{\Sigma }}& \boldsymbol 0\\ 0& \boldsymbol 1\end{array}\right] $

(7)

对$ {\mathit{\boldsymbol R}}_{\mathrm{G}} $做特征加权变换，得到增广加权盖尔圆矩阵$ {\mathit{\boldsymbol R}}_{\mathrm{F}} $，如式（8）所示：

$ {\mathit{\boldsymbol R}}_{\mathrm{F}}=\mathit{\boldsymbol W}{\mathit{\boldsymbol R}}_{\mathrm{G}}{\mathit{\boldsymbol W}}^{-1}= \\ \;\;\;\;\;\; \left[\begin{array}{ccccc}{\lambda }_{1}& 0& \cdots & 0& {\lambda }_{1}{p}_{1}\\ 0& {\lambda }_{2}& \cdots & 0& {\lambda }_{2}{p}_{2}\\ 0& 0& & ⋮& ⋮\\ ⋮& ⋮& & {\lambda }_{M}& {\lambda }_{M}{p}_{M}\\ \frac{1}{{\lambda }_{1}}{p}_{1}^{\mathrm{*}}& \frac{1}{{\lambda }_{2}}{p}_{2}^{\mathrm{*}}& \cdots & \frac{1}{{\lambda }_{M}}{p}_{M}^{\mathrm{*}}& 1\end{array}\right] $

(8)

令$ {r}_{m}={\lambda }_{m}{p}_{m} $，得到判别准则为：

$ \begin{array}{l}{k}^{*}=\underset{1<m<M}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}m\\ \mathrm{s}.\mathrm{t}.\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{r}_{m}-\frac{1}{M}\sum\limits_{i=1}^{M}\sqrt{{r}_{i}}<0\end{array} $

(9)

则利用WGDE准则得到的最终信源估计个数为：

$ {k}_{\mathrm{W}\mathrm{G}\mathrm{D}\mathrm{E}}={k}^{*}-1 $

(10)

3 增广加权盖尔圆矩阵多重特征融合及建模

在增广加权盖尔圆矩阵$ {\mathit{\boldsymbol R}}_{\mathrm{F}} $中包含了$ {p}_{m} $、$ {\lambda }_{m} $、$ {r}_{m} $3种用于信源个数估计的特征，分别表示盖尔圆心、盖尔圆半径和加权盖尔圆半径，其中，AIC准则利用$ {\lambda }_{m} $对信源个数进行估计^[6]，GDE准则利用$ {p}_{m} $对信源个数进行估计^[10]，而WGDE准则利用$ {r}_{m} $对信源个数进行估计^[16]，显然WGDE准则对$ {\mathit{\boldsymbol R}}_{\mathrm{W}} $中的特征应用并不充分。因此，本文将上述都能用于信源个数估计的3种特征进行融合，构建$ 3M $维的特征列向量$ \mathit{\boldsymbol \theta } $，如式（11）所示：

$ \mathit{\boldsymbol \theta }={\left[{\lambda }_{1}, {\lambda }_{2}, \cdots , {\lambda }_{M}, {p}_{1}, {p}_{2}, \cdots , {p}_{M}, {r}_{1}, {r}_{2}, \cdots , {r}_{M}\right]}^{\mathrm{T}} $

(11)

分别对式（11）中的各特征做归一化处理，如式（12）~式（14）所示：

$ {\dot{\lambda }}_{m}=\frac{{\lambda }_{m}}{\sum\limits_{i=1}^{M}\lambda {}_{{}_{i}}} $

(12)

$ {\dot{p}}_{m}=\frac{{p}_{m}}{\sum\limits_{i=1}^{M}{p}_{i}} $

(13)

$ {\dot{r}}_{m}=\frac{{r}_{m}}{\sum\limits_{i=1}^{M}{r}_{i}} $

(14)

其中，$ m=\mathrm{1, 2}, \cdots , M $。由此得到归一化后的特征向量$ \dot{\mathit{\boldsymbol \theta }}\in {\mathbb{R}}^{3M} $用于描述信源的个数，如式（15）所示：

$ \dot{\mathit{\boldsymbol \theta }}={\left[{\dot{\lambda }}_{1}, {\dot{\lambda }}_{2}, \cdots , {\dot{\lambda }}_{M}, {\dot{p}}_{1}, {\dot{p}}_{2}, \cdots , {\dot{p}}_{M}, {\dot{r}}_{1}, {\dot{r}}_{2}, \cdots , {\dot{r}}_{M}\right]}^{\mathrm{T}} $

(15)

给定一个包含$ N $个已标定的训练样本的数据集$ \left\{({\dot{\mathit{\boldsymbol \theta }}}_{1}, {y}_{1})\mathrm{ }, ({\dot{\mathit{\boldsymbol \theta }}}_{2}, {y}_{2})\mathrm{ }, \cdots , \right.\left.({\dot{\mathit{\boldsymbol \theta }}}_{n}, {y}_{n}), \cdots , ({\dot{\mathit{\boldsymbol \theta }}}_{N}, {y}_{N})\right\} $，其中，$ {\dot{\mathit{\boldsymbol \theta }}}_{n} $为根据式（11）~式（15）计算所得的第n个样本的特征向量，$ {y}_{n} $为对$ \dot \theta_{n} $的标定，且$ {y}_{n}\in \left\{\mathrm{1, 2}, \cdots , M-1\right\} $，$ n=\mathrm{1, 2}, \cdots , N $。由于本文中的模式分类为非线性，因此选用式（16）所示的高斯径向基函数作为核函数对原模式$ \dot \theta_{n} $做特征提取：

$ K\left(\phi \left({\dot{\theta }}_{n}\right)\mathrm{ }, \phi \left({\dot{\theta }}_{j}\right)\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-g{‖\phi \left({\dot{\theta }}_{n}\right)-\phi \left({\dot{\theta }}_{j}\right)‖}^{2}\right) $

(16)

在式（16）中，$ g $为高斯径向基函数的核参数，$ \phi :{\mathbb{R}}^{3M}\to {\mathbb{R}}^{\mathrm{S}} $，且有：

$ {\dot{\theta }}_{j}=\frac{\sum\limits_{n=1}^{N}{\dot{\theta }}_{n}}{N} $

(17)

本文模式分类的SVM设计通过寻找以下优化问题的解来实现：

$ \begin{array}{l}\underset{{\omega }_{m}\mathrm{ }, {b}_{m}}{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{ }\sum\limits_{n=1}^{N}\sum\limits_{m=1}^{M-1}\frac{1}{2}{‖{\omega }_{m}‖}^{2}+c\underset{m}{\mathrm{m}\mathrm{i}\mathrm{n}}{\epsilon }_{n, m}\\ \mathrm{s}.\mathrm{t}.\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\left|{\omega }_{m}^{\mathrm{T}}\phi \left({\dot{\theta }}_{n}\right)+{b}_{m}-{y}_{n}\right|\le {\epsilon }_{n, m}\end{array} $

(18)

其中，$ {\omega }_{m}, {b}_{m}\in {\mathbb{R}}^{\mathrm{S}} $，$ c $为惩罚因子，$ {\epsilon }_{n, m} $为松弛变量。由于式（18）中的最优解为$ ({\omega }_{m}^{*}, {b}_{m}^{*}) $，因此本文利用Libsvm工具包在寻优的过程中同时求解最优的参数对$ ({c}^{*}, {g}^{*}) $，从而提高对训练数据集的识别率。

4 仿真实验及分析 4.1 数学模型训练

实验采用4-UCA作为接收信号阵列，阵元间距为载波波长的1/2，3个远场窄带信号源的入射角随机设置，分别为$ ({0}^{\circ }, {25}^{\circ }) $、$ ({29}^{\circ }, {69}^{\circ }) $和$ ({60}^{\circ }, {106}^{\circ }) $，信噪比为$ -7\mathrm{ }\mathrm{d}\mathrm{B} $，采样快拍数$ L=120 $，在色噪声条件下^[19]，随机选取信源数分别为1、2、3时的各1 000个样本生成特征向量并标定，放入Libsvm中训练，令$ \mathrm{l}{\mathrm{b}}_{}g\in [-\mathrm{10, 10}] $、$ \mathrm{l}{\mathrm{b}}_{}c\in [-\mathrm{10, 20}] $，搜索步长为1，采用10倍交叉验证法^[20]寻找$ ({c}^{*}, {g}^{*}) $，寻优过程的二维等高线图和三维等高线图分别如图 1（a）和图 1（b）所示，其中，$ \mathrm{l}{\mathrm{b}}_{}{c}^{\mathrm{*}}=18 $，$ \mathrm{l}{\mathrm{b}}_{}{g}^{\mathrm{*}}=-4 $，此时的估计准确率约为$ 92.2\mathrm{\%} $，从而获得能够进行信源个数估计的4-UCA数学模型。

4.2 仿真实验

实验1  不同快拍数下信源个数估计的对比实验

选取4-UCA接收的包含1个、2个和3个信号源的测试样本各200个，且测试样本与训练样本的入射角不同，信噪比为$ {0}_{}\mathrm{d}\mathrm{B} $，快拍数$ L $从40增加到200，步长为20，将本文方法与文献[10]提出的GDE方法、文献[12]提出的改进GDE算法（此处记为GDE_AIC方法）、文献[16]提出的WGDE方法以及文献[17]方法（此处记为HHT_SVM方法）做对比实验，针对不同拍数下采集的测试样本，每个方法重复Monte Carlo实验^[21]200次。分别对含1个、2个和3个信号源的阵列信号做测试，得到5种方法在不同快拍数下的检测准确率。

实验2  不同信噪比下信源个数估计的对比实验

实验中的条件除了$ L=100 $，SNR范围为-20 dB~20 dB，步长为5 dB外，其他条件与实验1相同，得到5种方法对不同信噪比阵列信号的检测准确率。

4.3 结果分析

实验1的检测准确率对比如图 2所示，可以看出：当信号源个数为1和2时，在快拍数$ L\le 200 $条件下，GDE方法、MGDE方法、WGDE方法失效，GDE_AIC有一定的检测准确率但最高只有72%，而本文方法则表现出优异的检测性能，当$ L>90 $时，对信号源个数的检测准确率就已经超过了90%；当信号源为3时，GDE方法、MGDE方法、WGDE方法、GDE_AIC方法完全失效，本文方法在快拍数$ L=60 $时检测准确率就达到了80%；HHT_SVM方法只是在$ L\ge 100 $时才表现出检测性能优于本文方法。

实验2的检测准确率对比如图 3所示，可以看出：当信号源个数为1和2时，GDE方法和WGDE方法在信噪比大于等于$ {7}_{}\mathrm{d}\mathrm{B} $时才能获得80%以上的检测准确率，GDE_AIC方法在信噪比大于等于$ {5}_{}\mathrm{d}\mathrm{B} $时的检测准确率达到90%以上，而本文方法在信噪比为$ -{20}_{}\mathrm{d}\mathrm{B} $时检测准确率就已经达到80%，且检测准确率随着信噪比的提高而稳步提高；当信源个数为3时，GDE方法、WGDE方法、GDE_AIC方法失效，本文方法则始终保持较高的检测准确率；HHT_SVM方法的检测性能优越一直保持在80%以上，在信噪比大于等于$ {10}_{}\mathrm{d}\mathrm{B} $时其检测性能略优于本文方法。

5 结束语

GDE准则及现有改进方法大多只能估计比阵列天线阵元个数少2个的信源个数。为弥补这一缺陷，本文提出一种基于多重特征融合的信源个数估计方法，从WGDE准则得到的增广加权盖尔圆盘矩阵中同时获取盖尔圆心、盖尔圆半径和加权盖尔圆半径等多重特征进行融合，构建可描述信源个数的特征向量，并利用SVM训练分类器数学模型用于信源个数估计。基于4-UCA的仿真结果表明，本文方法不仅能够准确估计仅比阵元数少1的信源个数，而且在低信噪比和小快拍数条件下也具有良好的估计性能。本文未考虑信号源数目大于或等于阵元数目的情况，后续将从这一角度出发对信源估计方法做进一步探索。