一种基于压缩感知的无需超参数的测向算法

A Hyperparameter-Free Direction-Finding Algorithm in the Compressive Sensing Framework

摘要: 全自动的稀疏参数估计算法不需要用户对于超参数的数值做出任何艰难的决定(可能是试错而来)，因而更加实用。本文给出了现有方法的统一阐释，包括协方差矩阵拟合(CMF)、基于稀疏迭代协方差的估计(SPICE)以及基于似然的稀疏参数估计(LIKES)，认为它们全部都是不同统计距离下的基于协方差矩阵拟合的方法。在此基础上，本文提出了一种新的协方差矩阵拟合方案，将两种非对称Itakura-Saito距离的其中一种最小化。计算机仿真表明，本文提出方法相比上述方法的性能具有优势。

Abstract: The fully automatic sparsity-parameter estimation algorithms do not require the user to make any hard decision (possibly via trial-and-error) on the values of the hyperparameters, making them more pragmatic in practice. This paper provides a unified interpretation of the existing approaches including covariance matrix fitting (CMF), sparse iterative covariance based estimation (SPICE) and likelihood-based estimation of sparse parameters (LIKES). The point of view taken is that they are all covariance-fitting-based algorithms under different statistical distances. Following this, we present a new covariance-fitting scheme trying to minimize one of the two asymmetrical Itakura-Saito distances. Simulations show that the proposed method appears to be preferable as it outperforms the aforementioned algorithms in general.

1. 引言

频谱监测是频谱管理的重要工具，主要目的是完成对无线电频谱的监测和无线电信号参数的测量，其中，无线电信号的测向是其主要任务之一。现代测向技术包括旋转天线、多普勒、沃森–瓦特、相位干涉仪、相关干涉仪、高分辨率测向等等。本文讨论的是高分辨率测向技术，该项技术主要用来分辨相同频率的多个发射机的信号的方向 [1] 。

代表性的高分辨率测向技术有波束形成法 [2] 和子空间方法 [3] 。近年来，基于压缩感知的方法由于相较前两类方法进一步提升了估计精度而逐渐受到关注 [4] 。从数学上看，测向设备对空间谱的感知是一种欠定方程，即：以比较少的天线单元个数来分辨来自于各个方向的信号的能量。欠定方程的解通常有无穷多个，为了减少解的数量，提高对空间谱的分辨能力，需要对方程施加更多的约束。压缩感知(compressive sensing)，亦称为压缩传感，是一种利用未知矢量的稀疏性(sparsity)，对解进行约束的方法。稀疏性指的是一个矢量只有少数的几个元素不为零，大多数是零元素。空间谱符合稀疏性这一特征：对于某一个频点/信道而言，空间中只存在有限个方向的同频信号，大部分的方向是没有该频点/信道的信号的。该类方法的最新趋势是摆脱算法用户的人工干预，无需用户设定超参数(hyperparameter)从而实现全自动计算 [5] [6] [7] 。本文的目的即是提出一种基于压缩感知的无需超参数的测向算法，以实现高精度同频信号测向。

2. 数学模型

考虑一个N个天线单元和接收通道组成的测向设备接收到了M个同频窄带信号 ${\left\{s_m\left(t\right)\right\}}_{m=1}^M\in ℂ$ ，其导向矢量分别为 ${\left\{a_m\right\}}_{m=1}^M\in {ℂ}^{N\times 1}$ ，背景噪声为 $n\left(t\right)\in {ℂ}^{N\times 1}$ 。阵列测量矢量由下式给出：

$x\left(t\right)={\displaystyle {\sum }_{m=1}^M{a}_m{s}_m\left(t\right)}+n\left(t\right)$ (1)

令 ${\left\{{\theta }_k\right\}}_{k=1}^K$ 为可以覆盖所有可能入射区域的栅格，并覆盖了(至少是近似覆盖了)上述M个信号的入射角度，相应的导向矢量为，其中， $K\gg N$ 。于是，式(1)可以重新表述为：

$x\left(t\right)={\displaystyle {\sum }_{k=1}^K{d}_k{\bar{s}}_k\left(t\right)}+n\left(t\right)$ (2)

其中， ${\bar{s}}_k\left(t\right)$ 代表来自于 ${\theta }_k$ 方向的可能存在的辐射源的信号。基于压缩感知的到达角估计算法的目标即是从K个可能存在的辐射源中找出M个真实存在的辐射源。

上述问题可以采用约束l_p范数最小化方法解决，例如LASSO [8] 、FOCUSS [9] [10] 、l1-ACCV [11] 等等。然而，因为没有用来确定上述方法中涉及的拟合门限、正则化参数等超参数的规则，上述方法在实际应用中常常难以实施。最近，学者们提出了一系列基于协方差矩阵拟合 [12] 的稀疏参数估计方法，由于这些方法不需要设定超参数，因而可以全自动运行。

假设入射信号之间相互独立，背景噪声为空时白噪声，那么阵列协方差矩阵的稀疏表示为：

$R=E\left\{x\left(t\right)x^{\text{H}}\left(t\right)\right\}={\displaystyle {\sum }_{k=1}^K{\sigma }_k{d}_k{d}_k^{\text{H}}}+ϵI$ (3)

其中， ${\sigma }_k=E\left\{{\left|{\bar{s}}_k\left(t\right)\right|}^2\right\}$ 为第k个可能存在的辐射源信号的功率， $ϵI=E\left\{n\left(t\right)n^{\text{H}}\left(t\right)\right\}$ 为噪声协方差矩阵， ${(\cdot )}^{\text{H}}$ 代表共轭转置。为了考虑噪声功率不均匀的情况，可以定义 $\left[d_{K+1},d_{K+2},\cdots ,d_{K+N}\right]=I$ ， ${\sigma }_{K+q}=E\left\{{\left|n_q\left(t\right)\right|}^2\right\}$ ， $q=1,2,\cdots ,N$ ，所以有：

$R={\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k{d}_k{d}_k^{\text{H}}}$ (4)

同时，对于给定的L组阵列采样 ${\left\{x\left[l\right]\right\}}_{l=1}^L$ ，阵列采样协方差矩阵的估计公式为：

$\hat{R}=\frac{1}{L}{\displaystyle {\sum }_{l=1}^{L}x\left[l\right]x^{\text{H}}\left[l\right]}$ (5)

2.1. 协方差矩阵拟合算法(CMF)

协方差矩阵拟合算法(CMF)是在均匀噪声功率的假设下，将理论协方差矩阵 $R$ 与采样协方差矩阵 $\hat{R}$ 之间的欧式距离(SED)最小化 [5] ，即

$\begin{array}{l}\underset{{\left\{{\sigma }_k\right\}}_{k=1}^K,ϵ}{\min }\left\{d_{\text{SE}}\left(\hat{R},R\right)={\Vert \hat{R}-{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\sigma }_k{d}_k{d}_k^{\text{H}}}-ϵI\Vert }^2\right\}\\ \text{s}.\text{t}.{\left\{{\sigma }_k\right\}}_{k=1}^K\ge 0,{\displaystyle {\sum }_{k=1}^K{\sigma }_k}\le \lambda ,\lambda =-N{\gamma }_N+{\displaystyle {\sum }_{n=1}^N{\gamma }_n},ϵ\ge 0\end{array}$ (6)

其中， ${\gamma }_1\ge {\gamma }_2\ge \cdots \ge {\gamma }_N$ 为 $\hat{R}$ 的特征值。如果我们使用非均匀噪声功率模型，上述问题将会简化为：

$\begin{array}{l}\underset{{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}}{\min }\left\{d_{\text{SE}}\left(\hat{R},R\right)={\Vert \hat{R}-{\displaystyle \underset{k=1}{\overset{K+N}{\sum }}{\sigma }_k{d}_k{d}_k^{\text{H}}}\Vert }^2\right\}\\ \text{s}.\text{t}.{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}\ge 0\end{array}$ (7)

2.2. 基于稀疏迭代协方差的估计(SPICE)

SPICE算法 [6] 实际上是将两种高斯分布和 $$ 之间的Jeffreys距离(亦称为对称Kullback-Leibler距离) [13] 最小化：

$\begin{array}{l}\underset{{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}}{\min }\left\{{\Vert R^{-1/2}\left(\hat{R}-R\right){\hat{R}}^{-1/2}\Vert }^2=\text{tr}\left(R^{-1}\hat{R}\right)+\text{tr}\left({\hat{R}}^{-1}R\right)-2N=2d_{KL}\left(\hat{R},R\right)\right\}\\ \text{s}.\text{t}.R={\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k{d}_k{d}_k^{\text{H}}},{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}\ge 0\end{array}$ (8)

对于非高斯分布而言，上述拟合公式仍然是一种度量两个协方差矩阵差异的优秀标准。论文 [6] 中有关随机相位信号的仿真结果证明了这一点。

2.3. 基于似然的稀疏参数估计(LIKES)

LIKES算法 [7] 使用的是高斯最大似然(GML)标准：

$\begin{array}{l}\underset{{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}}{\min }\left\{X^{\text{H}}{R}^{-1}X+\text{ln}\left|R\right|\right\}\\ \text{s}.\text{t}\text{.}R={\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k{d}_k{d}_k^{\text{H}}},{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}\ge 0\end{array}$ (9)

其中， $X=\left[x\left[1\right],x\left[2\right],\cdots ,x\left[L\right]\right]$ 。论文 [14] 证明了对于参数估计而言，GML代价函数最小化等价于两种高斯分布 ${\mathcal{G}}_1\left(0,R\right)$ 和 ${\mathcal{G}}_2\left(0,\hat{R}\right)$ 之间的两种非对称Itakura-Saito距离(ISD)的其中一种最小化，由下式给出：

$\begin{array}{l}\underset{{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}}{\min }\left\{d_{\text{IS}}\left(\hat{R},R\right)=\text{tr}\left(\hat{R}{R}^{-1}-I\right)-\text{ln}\left|\hat{R}{R}^{-1}\right|\right\}\\ \text{s}.\text{t}\text{.}R={\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k{d}_k{d}_k^{\text{H}}},{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}\ge 0\end{array}$ (10)

和式(8)一样，式(9)中的拟合标准对于非高斯分布来说仍然有效。然而，式(9)中的优化问题不是凸优化问题，需要利用SPICE算法进行迭代优化 [7] 。

3. 本文提出的协方差矩阵拟合方案

本文提出利用另外一种非对称Itakura-Saito距离(ISD)用作拟合标准，即

$\begin{array}{l}\underset{{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}}{\min }\left\{d_{\text{IS}}\left(R,\hat{R}\right)=\text{tr}\left(R{\hat{R}}^{-1}-I\right)-\text{ln}\left|R{\hat{R}}^{-1}\right|\right\}\\ \text{s}.\text{t}\text{.}R={\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k{d}_k{d}_k^{\text{H}}},{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}\ge 0\end{array}$ (11)

经过简单的计算可知，

$d_{\text{IS}}\left(R,\hat{R}\right)=\text{tr}\left(R{\hat{R}}^{-1}\right)-\text{ln}\left|R\right|+\text{const}={\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k{d}_k^{\text{H}}{\hat{R}}^{-1}{d}_k}-\text{ln}\left|{\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k{d}_k{d}_k^{\text{H}}}\right|+\text{const}$ (12)

容易验证，对于 ${\left\{{\sigma }_k\right\}}_{k=1}^{K+N}\ge 0$ ， ${\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k}>0$ ，有 ${\displaystyle {\sum }_{k=1}^{K+N}{\sigma }_k{d}_k{d}_k^{\text{H}}}\succ 0$ 。因此， $d_{\text{IS}}\left(R,\hat{R}\right)$ 的最小化问题属于“max-det”凸优化问题 [15] 。如同LIKES方法，该方法也属于最大似然方法 [16] 。

有关鲁棒自适应波束形成的研究表明，直接使用采样协方差矩阵的逆矩阵 ${\hat{R}}^{-1}$ 可能会导致算法性能下降，尤其是在采样点数较小的情况下 [14] 。为了解决这一问题，我们将 $\hat{R}$ 替换为基于广义线性组合(GLC) [17] 的收缩型增强估计：

$\bar{R}=\hat{\alpha }I+\hat{\beta }\hat{R}$ (13)

其中，

$\hat{\rho }=\frac{1}{L^2}\underset{l=1}{\overset{L}{{\displaystyle \sum }}}{\Vert x\left[l\right]\Vert }^4-\frac{1}{L}{\Vert \hat{R}\Vert }^2$

$\hat{\nu }=\frac{\text{tr}\left(\hat{R}\right)}{N}$

$\hat{\beta }=1-\frac{\hat{\alpha }}{\hat{\nu }}$ (14)

因此，本文提出的方法可以表述为：

$\begin{array}{l}\underset{{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}}{\min }\left\{\underset{k=1}{\overset{K+N}{{\displaystyle \sum }}}{\sigma }_k{d}_k^{\text{H}}{\hat{R}}^{-1}{d}_k-\text{ln}\left|\underset{k=1}{\overset{K+N}{{\displaystyle \sum }}}{\sigma }_k{d}_k{d}_k^{\text{H}}\right|\right\}\\ \text{s}.\text{t}.{\left\{{\sigma }_k\right\}}_{k=1}^{K+N}\ge 0\end{array}$ (15)

式(15)可以利用“max-det”问题专用数学工具MAXDET [18] 、通用凸优化数据工具CVX [19] 等来实现优化计算。例如，当使用CVX数学工具时，可以定义：

$c={\left[d_1^{\text{H}}{\hat{R}}^{-1}{d}_1,d_2^{\text{H}}{\hat{R}}^{-1}{d}_2,\cdots ,d_{K+N}^{\text{H}}{\hat{R}}^{-1}{d}_{K+N}\right]}^{\text{T}}$ (16)

$p={\left[{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_{K+N}\right]}^{\text{T}}$ (17)

$D=\left[d_1,d_2,\cdots ,d_{K+N}\right]$ (18)

那么，式(15)的CVX优化代码为：

cvx_begin quiet

cvx_precision best

cvx_expert true

variable p(K+N,1)

minimize (c'*p-log_det(D*diag(p)*D'))

subject to

p >= zeros(K+N,1);

cvx_end

对矢量 $p$ 进行局部极值搜索，即可获得测向结果。

4. 计算机仿真

本节给出计算仿真示例来验证本文提出算法的性能。仿真场景为一个10阵元的均匀线性阵列(阵列间距为半个波长)，接收三个同频率的远场窄带信号(分别来自−5˚，15˚和25˚)，背景噪声为加性高斯白噪声。使用步进为0.1˚的均匀栅格来覆盖可能的入射角度(−90˚,90˚]，即 $K=1800$ 。图1展示了本文提出的算法与CMF、SPICE、LIKES方法的均方根误差(RMSE)随着信噪比(SNR)变化的曲线，其中采样点数固定为200。图2则展示了上述算法的RMSE随着采样点数变化的曲线，其中信噪比固定为10 dB。图1和图2的结果均为300次蒙特卡洛随机实验的平均结果，RMSE利用下式进行计算：

$\text{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{n=1}^{300}\left[{\left({\hat{\theta }}_{1,n}-{\theta }_1\right)}^2+{\left({\hat{\theta }}_{2,n}-{\theta }_2\right)}^2+{\left({\hat{\theta }}_{3,n}-{\theta }_3\right)}^2\right]}{900}}$ (19)

其中， ${\hat{\theta }}_{1,n}$ 、 $$ 、 ${\hat{\theta }}_{3,n}$ 分别是对 $$ 、 ${\theta }_2=-{15}^{\circ }$ 、 ${\theta }_3={25}^{\circ }$ 的第n次估计结果。

Figure 1. RMSE versus varied SNR

图1. 均方根误差随信噪比的变化曲线

Figure 2. RMSE versus varied numbers of samples

图2. 均方根误差随采样数的变化曲线

从图1和图2可以看出，本文提出算法在参与比较的算法中是最准确的。它的性能与LIKES算法相仿，略有提升。这两种算法是基于非对称Itakura-Saito距离的两个相反的方向：一个衡量的是 ${\mathcal{G}}_1$ 到 ${\mathcal{G}}_2$ 的收敛度，另一个衡量的是 ${\mathcal{G}}_2$ 到 ${\mathcal{G}}_1$ 的收敛度。需要指出的是，LIKES算法的优化是一个非凸优化问题，在某些随机实验中难以做到全局收敛，而本文算法的优化是一个凸优化问题，其全局收敛性得到了保证。这最终造成了两者在角度估计精度上的差别。

5. 结论

本文针对数种基于协方差拟合的稀疏系数估计算法，给出了一种统一的基于收敛距离的解释，并在此基础上提出了一种利用Itakura-Saito距离的新方法，该方法的估计性能优于参与比较的方法。鉴于收敛距离的种类繁多，在将来还有许多研究需要做。举例来说，Bhattacharyya距离 [20] 和Wasserstein度量 [21] 不需要用户指定参数，具有作为拟合标准的潜力。