摘要——估计入射到传感器阵列上的信源数量是一个众所周知且被深入研究的问题。解决该问题的一个常用方法是使用信息论准则，例如最小描述长度（Minimum Description Length，MDL）或赤池信息准则（Akaike Information Criterion，AIC）。众所周知，MDL估计器是一种一致性估计器，它对于偏离高斯假设是鲁棒的，但对于偏离点源和/或时间或空间上的加性白噪声假设则不鲁棒。多年来，研究人员已经提出并测试了几种替代的估计算法。通常，这些算法相比MDL估计器具有更优的性能，并且对于偏离假设的空间模型也更具鲁棒性。然而，这些鲁棒算法具有很高的计算复杂度，需要进行多次多维搜索。

本文受现实生活问题的启发，采用一种系统性的方法来解决使用信息论准则鲁棒地估计信源数量的问题。文中研究了一种MDL类型的估计器，该估计器对于偏离“阵列上各处噪声水平相等”这一假设具有鲁棒性。该估计器的一致性得到了证明，即使在偏离相等噪声水平假设的情况下也是如此。文中还提出了一种新颖的低复杂度实现方法，该方法避免了多维搜索的需要，使得这种估计器成为实际应用中的一个有利选择。

I. INTRODUCTION

对信源数量的估计器，人们抱有相当大的兴趣。这些估计器不应需要先验知识，并且即使在偏离假设模型的情况下也应保持一致性。本文针对偏离假设模型的特定情况开发了此类估计器。我们特别考虑了传感器噪声水平在空间上非均匀的情况。文中将表明，尽管传统的信源数量估计方法在这种情况下倾向于过高估计信源数量，但我们提出的估计器没有这种倾向。

已有许多论文详细介绍了对MDL估计器的改进，即 [5], [6], [9], [14], [15], [16]，注意这只是此类工作的一部分。在许多这些工作中，采用了MDL方法，并通过利用某种类型的先验知识来实现性能提升[6], [10], [17]。通常所做的一个假设是加性噪声过程是一个空间白过程，而所提出的方法对偏离此假设的鲁棒性通常通过计算机仿真来评估[9]。总的来说，可以观察到，使用某种先验信息的方法是鲁棒的，而基于最小特征值重数的方法则不鲁棒。后一种估计器缺乏鲁棒性的原因是，当偏离假设模型时，最小特征值的重数等于1 [1]。因此，无法从接收信号相关矩阵的最小特征值的重数来推断信源的数量。另一方面，基于某些先验知识（例如，阵列导向矢量）的方法通常具有很高的计算复杂度，需要进行多次多维数值搜索[18]。此外，当发生某些对假设模型的偏离时，这些方法不一定是一致的，尽管它们在仿真中表现出良好的鲁棒性。

高效且鲁棒地估计信源数量在生物医学应用中非常重要（例如，参见[19]及其参考文献）。例如，在一种此类应用中，研究人员感兴趣的是估计对短暂刺激作出反应的神经元数量。这是通过在患者头部放置一个非常大的传感器阵列并记录这些传感器接收到的脑活动来完成的。在这些生物医学问题中，不存在先验知识（例如，关于导向矢量的知识）。此外，由于不同的传感器与患者皮肤的距离略有不同，传感器输出端的噪声水平差异很大。因此，生物医学应用是加性噪声不一定是空间均匀的一类重要问题的例子。

尽管存在鲁棒的信源数量估计器，但这些估计器需要一些通常无法获得的先验知识，并且如上所述，它们的计算复杂度很高。因此，计算上高效且鲁棒的信源数量估计器引起了人们相当大的兴趣。这些估计器不应需要先验知识，并且即使在偏离假设模型的情况下也应保持一致性。

本文针对偏离假设模型的特定类型开发了此类估计器。我们特别考虑了传感器噪声水平在空间上非均匀的情况。文中将表明，虽然传统的信源数量估计方法在这种情况下倾向于过高估计信源数量，但我们提出的估计器没有这种倾向。

B. Problem Formulation

考虑一个由 $p$ 个传感器组成的阵列，并用$\mathbf{x}(t)$表示在时刻 $t$ 接收到的$p$维信号矢量。用$q<p$表示入射到阵列上的信号数量。接收信号矢量的一个通用模型为 [11], [18]：

$$\mathbf{x}(t)=\mathbf{As}(t)+\mathbf{n}(t),t=1,2,\dots ,N\text{\hspace{1em}(1)}$$

其中 $\mathbf{A}=[\mathbf{a}({\psi }_1),\mathbf{a}({\psi }_2),\dots ,\mathbf{a}({\psi }_q)]$是一个由$q$个$p$维矢量组成的$p\times q$矩阵，且$\mathbf{a}(\psi )$位于阵列流形 { A = a ( ψ ) ∣ ψ ∈ Ψ } \{\mathcal{A}=\mathbf{a}(\psi) \mid \psi \in \mathbf \Psi\} {A=a(ψ)∣ψ∈Ψ}上，这里$\mathbf{\Psi }]$表示描述阵列响应的一组参数。$\mathbf{a}(\psi )$被称为阵列响应矢量或导向矢量，$\mathbf{A}$被称为导向矩阵，而${\psi }_i$是与第$i$个信源相关的未知参数向量。$\mathbf{s}(t)=[s_1(t)\dots s_q(t){]}^T$是一个零均值、具有正定相关矩阵${\mathbf{R}}_s$的白色复平稳高斯随机过程；$\mathbf{n}(t)$是一个时间白色的复高斯随机矢量过程，与信号无关，其均值为零，相关矩阵为$\text{diag}([{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_p^2])$，其中$\text{diag}([{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_p^2])$表示一个对角矩阵，其对角线上的元素为矢量$[{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_p^2]$。这个相关矩阵代表了每个传感器可能面临不同噪声水平的场景。

面临不同的噪声水平。定义${\sigma }^2=(1/p){\sum }_{i=1}^p{\sigma }_i^2$，$w_i={\sigma }_i^2-{\sigma }^2$，以及$\mathbf{w}=[w_1,w_2,\dots ,w_p]$。借助${\sigma }^2$和$\mathbf{w}$，加性噪声的相关矩阵可以描述如下： E { n ( t ) n H ( t ) } = σ 2 I + diag ( w ) E\{\mathbf{n}(t)\mathbf{n}^H(t)\} = \sigma^2\mathbf{I} + \text{diag}(\mathbf{w}) E{n(t)nH(t)}=σ2I+diag(w)。这种替代表示简化了后文中的一些证明和推导。注意，矢量$\mathbf{w}$代表了对“噪声水平在整个阵列上恒定”这一假设的偏离。最后，假设导向矩阵$\mathbf{A}$的所有元素都是未知的 [1]，唯一的限制是$\mathbf{A}$为满秩矩阵。在后文中，高斯假设将会被放宽。

我们用 ${\mathbf{\theta }}_q$ 表示假设有 $q$ 个信源时的未知参数集，即 ${\mathbf{\theta }}_q=[{\mathbf{R}}_{s,q},{\mathbf{A}}_q,{\sigma }_{n,q}^2,{\mathbf{w}}_q]$，其中${\mathbf{R}}_{s,q}$是假设有$q$个信源时的发射信号相关矩阵；${\mathbf{A}}_q$是假设有$q$个信源时的导向矩阵；${\sigma }_{n,q}^2$是白噪声水平；而${\mathbf{w}}_q$是包含了表示与空间白噪声假设偏离程度的参数的矢量。给定$q$个信源时，未知参数的参数空间（parameter space）用 ${\mathbf{\Theta }}_q$ 表示。问题是基于阵列输出的 $N$ 个独立快照${\mathbf{x}}_1=\mathbf{x}(t_1),\dots ,{\mathbf{x}}_N=\mathbf{x}(t_N)$来估计 $q$ [1]。

C. Information Theoretic Criteria and MDL Estimators

信息论准则（Information Theoretic Criterion，ITC）是一种用于在多个竞争的参数化模型中进行选择的估计准则[3]。

给定一个参数化的概率密度族 $f_{\mathbf{X}}(\mathbf{X}|{\mathbf{\theta }}_q)$，其中${\mathbf{\theta }}_q\in {\mathbf{\Theta }}_q$、$\mathbf{X}=[{\mathbf{x}}_1,\dots ,{\mathbf{x}}_N]$，对于不同的$q$，一个ITC估计器选择$\hat{q}$以满足[3]：

$${\hat{q}}_{\mathrm{ITC}}=\arg \underset{q}{\min }\mathrm{ITC}(q)=\arg \underset{q}{\min }\left\{-L(\hat{\mathbf{\theta }}{)}_q+\mathrm{penalty}(q)\right\}\text{\hspace{1em}(2)}$$

其中$L({\theta }_q)\triangleq \log f_{\mathbf{X}}(\mathbf{X}|{\theta }_q)$是测量值的对数似然，${\hat{\mathbf{\theta }}}_q=\arg {\max }_{{\mathbf{\theta }}_q\in {\mathbf{\Theta }}_q}{f}_{\mathbf{X}}(\mathbf{X}|{\mathbf{\theta }}_q)$是在给定第$q$个分布族下未知参数的最大似然（ML）估计，而 $\text{penalty}(q)$是与所使用的特定ITC相关联的某个通用惩罚函数。

MDL和AIC估计器分别由 $\mathrm{penalty}(q)=\left|{\mathbf{\Theta }}_q\right|\log (N)/2$ 和 $\text{penalty}(q)=|{\mathbf{\Theta }}_q|$给出，其中$|{\theta }_q|$是${\theta }_q$中自由参数的数量[1], [20], [21]。

众所周知，在渐近情况下并满足某些正则性条件时，MDL估计器会最小化测量值$\mathbf{X}$和模型${\hat{\theta }}_q$两者的描述长度（以比特为单位）[22]，而AIC准则则最小化不同模型与真实模型之间的库尔贝克-莱布勒散度。在本文的其余部分，我们将只考虑MDL准则，尽管其他的惩罚函数也可以用类似的方式处理。

尽管在许多与阵列处理相关的问题中，例如波达方向（DOA）估计，当估计信源数量时人们对阵列结构有一定的先验知识，但这些先验知识通常被忽略了[1], [18]。其原因是，通过忽略阵列结构并假设信号和噪声为高斯分布，且 $\mathbf{w}\equiv 0$，所得到的MDL估计器（2）（在此称为高斯-MDL（Gaussian-MDL，GMDL）估计器[11]）具有一个由[1]给出的简单闭式表达式：

$$\begin{array}{rl}{\hat{q}|}_{\mathrm{GMDL}}=\arg \underset{q=0,\dots ,p-1}{\min }[-& N\log \frac{{\prod }_{i=q+1}^p{l}_i}{{\left(\frac{1}{p-q}{\sum }_{i=q+1}^p{l}_i\right)}^{p-q}}\\ & +\frac{1}{2}(q(2p-q)+1)\log N]\end{array}\text{\hspace{1em}(3)}$$

其中 $l_1\ge l_2\ge \cdots \ge l_p$是经验接收信号相关矩阵（empirical received signal’s correlation matrix）$\hat{\mathbf{R}}=(1/N)\sum {\mathbf{x}}_i{\mathbf{x}}_i^H$的特征值。众所周知，

• 当 $\mathbf{w}\equiv 0$ 时，GMDL估计器是信源数量的一致性估计器；
• 而当$\mathbf{w}\ne 0$时，GMDL估计器（3）则不具有一致性，

实际上，在这种情况下，随着快照数量趋于无穷大，GMDL估计器所产生的错误概率会趋近于1 [11]。

用 ${\mathcal{R}}_q$ 表示所有正定、秩为$q$的$p\times p$厄米特矩阵的集合，用$\mathcal{W}$表示所有零均值、$p$维向量的集合。在问题陈述所作的假设下，用于估计信源数量的MDL估计器（下文称作鲁棒MDL（RMDL）估计器）由显示在下一页底部的公式（4）给出，其中${\hat{\mathbf{A}}}_q,{\hat{\mathbf{R}}}_{s,q},{\hat{\sigma }}_{n,q}^2$和${\hat{\mathbf{w}}}_q$是在假设有$q$个信源时未知参数的最大似然（ML）估计，即，我们有显示在下一页底部的公式（5）。

注意，由于 ${\mathbf{A}}_q{\mathbf{R}}_{s,q}{\mathbf{A}}_q^H\in {\mathcal{R}}_q$，通过使用特征分解，我们可以写出${\mathbf{A}}_q{\mathbf{R}}_{s,q}{\mathbf{A}}_q^H={\sum }_{i=1}^q{\lambda }_i{\mathbf{v}}_i{\mathbf{v}}_i^H$，其中 { v i } \{\mathbf{v}_i\} {vi​}是一组标准正交向量。因此，假设有$q$个信源时的未知参数向量也可以由下式给出：
$${\mathbf{\theta }}_q=[{\lambda }_1,\dots ,{\lambda }_q,{\mathbf{v}}_1^T,\dots ,{\mathbf{v}}_q^T,{\mathbf{w}}^T,{\sigma }_n^2].\text{\hspace{1em}(6)}$$

II. IDENTIFIABILITY AND CONSISTENC OF THE RMDL ESTIMATOR

A. Identifiability

考虑一个参数化的概率密度函数（pdf）族 $f_{\mathbf{X}}(\mathbf{x}|\mathbf{\theta }),\mathbf{\theta }\in \mathbf{\Theta }$。如果对于任意的 $\mathbf{\theta }\ne {\mathbf{\theta }}^{\prime }$，$f_{\mathbf{X}}(\mathbf{x}|\mathbf{\theta })$ 和 $f_{\mathbf{X}}(\mathbf{x}|{\mathbf{\theta }}^{\prime })$ 之间的库尔贝克-莱布勒散度（Kullback–Leibler divergence）都大于零，那么该密度族就被称为是可辨识的（identifiable）。也就是说，$D(f_{\mathbf{X}}(\mathbf{x}|\mathbf{\theta })||f_{\mathbf{X}}(\mathbf{x}|{\mathbf{\theta }}^{\prime }))>0$，其中 $D(f(\mathbf{x})||g(\mathbf{x}))=\int f\log (f/g)$ 是 $f(\mathbf{x})$ 和 $g(\mathbf{x})$ 之间的库尔贝克-莱布勒散度 [23]。这个条件确保了参数空间与测量值的统计特性之间存在一一对应的关系。

I-B节中讨论的问题是一个模型阶数选择问题 [22]。如果对于某个 $k\ne l$，能够在参数空间中找到两个点 ${\mathbf{\theta }}_k\in {\mathbf{\Theta }}_k$ 和 ${\mathbf{\theta }}_l\in {\mathbf{\Theta }}_l$，使得 $f(\cdot |{\mathbf{\theta }}_k)=f(\cdot |{\mathbf{\theta }}_l)$，那么这个问题就是不可辨识的。不幸的是，我们实际上可以找到这样两个点，从而得出结论：I-B节中讨论的估计问题是不可辨识的。接收信号的pdf完全由接收信号的相关矩阵来表征。因此，为了证明该问题是不可辨识的，只需找到两个不同的参数值，在这两个参数值下，对应的接收信号相关矩阵是相等的。例如，考虑以下接收信号相关矩阵：$\text{diag}([11,10.5,9.5,10])$。

• 该相关矩阵可能由一个只有噪声的场景产生，其中 ${\sigma }_n^2=10.25$ 并且$\mathbf{w}=[0.75,0.25,-0.75,-0.25]$，
• 也可能由一个单信源场景产生，其中 ${\sigma }_n^2=10,\mathbf{w}=[0,0.5,-0.5,0],\mathbf{a}=[1,0,0,0{]}^T$，并且${\mathbf{R}}_s=1$。

因此，我们找到了两种场景：第一种对应于纯噪声场景，另一种对应于单信源场景，在这两种场景下接收信号矢量的分布是相同的。因此，这个例子表明所构建的估计问题是不可辨识的。

为了使估计问题变得可辨识，所有具有相同接收信号pdf的点，除了一个之外，都必须从参数空间中移除。按照模型阶数选择问题的惯例，在所有具有相同接收信号pdf的点中，将信源数量最少的那个点，即未知参数数量最少的点，保留在参数空间中，其余的点则被删除。由此产生的主要问题是，参数空间中的大部分点是否是可辨识的。幸运的是，这个问题的答案是肯定的；也就是说，参数空间中的大部分点是可辨识的。下面的引理描述了参数空间中所有不可辨识的点。

引理1： 假设 $q<p$。那么，${\mathbf{\theta }}_q$ 是参数空间中的一个不可辨识点，当且仅当对于某个 $j\in [1,\dots ,p]$，矩阵 ${\sum }_{i=1}^q{\lambda }_i{\mathbf{v}}_i{\mathbf{v}}_i^H$ 包含 $\alpha {\mathbf{e}}_j=[0_{j-1},1,0_{p-j}]$ 作为其第j行，其中 ${\mathbf{v}}_i$ 在(6)中定义。

引理1的证明： 见附录A。

引理1的证明为不可辨识点提供了一个有趣的物理解释。特别地，从引理的证明中可以看出，所有不可辨识点都与上述用于证明问题不可辨识的例子中的点相似。也就是说，一个不可辨识点对应于这样一种场景：存在（比方说）$q$个信源，其中一个信源仅被其中一个传感器接收到。由于这个信源无法与相应元件上噪声水平某个标称值的偏差区分开来，这种情况可能会与另一个信源数量更少且相应元件上噪声水平增加的场景相混淆。从实践的角度来看，这种情况是罕见的。