MUSIC算法原理与信号DOA估计

一、平稳随机过程的自相关矩阵及其性质

1.1 自相关矩阵的定义

​ 对离散时间平稳随即构成，用$M$个时刻的随机变量$u(n),u(n-1),...,u(n-M+1)$构造随机向量

​ $$u(n)=[u(n),u(n-1),...,u(n-M+1){]}^T$$ $(\mathrm{1.1.1})$

​ 随机过程$u(n)$的自相关矩阵（correlation matrix）（简称相关矩阵）定义为

​ $$R=E[u(n)u^H(n)]$$ $(\mathrm{1.1.2})$

​ 将式（1.1.1）代入式（1.1.2），并考虑平稳条件，得到相关矩阵的展开式为

​ $$R=\left[\begin{array}{cccc}r(0)& r(1)& ...& r(M-1)\\ r(-1)& r(0)& ...& r(M-2)\\ ...& ...& ...& ...\\ r(-M+1)& r(-M+2)& ...& r(0)\end{array}\right]\in C^{M\times M}$$

​ 式中，$r(m)$是随机过程$u(n)$的自相关函数，为$r(m)=E[u(n)u^{\ast }(n-m)]$

​ 根据相关函数共轭对称性，即$r(-m)=r^{\ast }(m)$，上式可重写为

​ $$R=\left[\begin{array}{cccc}r(0)& r(1)& ...& r(M-1)\\ r^{\ast }(1)& r(0)& ...& r(M-2)\\ ...& ...& ...& ...\\ r^{\ast }(M-1)& r^{\ast }(M-2)& ...& r(0)\end{array}\right]$$

​ 因此，对于一个平稳随机过程，只需自相关函数$r(m)(m=0,1,...,M-1)$的$M$个值就可以完全确定相关矩阵$R$。

1.2 自相关矩阵的基本性质

​ 自相关矩阵在离散时间统计信号处理中具有极其重要的作用，由式$(\mathrm{1.1.2})$给出的定义，可以得到平稳离散时间随机过程相关矩阵的一些基本性质。

​ 性质1 平稳离散时间随机过程的相关矩阵是$Hermite$矩阵，即有$R^H=R$。

​ 注：$Hermite$矩阵又称作自共轭矩阵、埃尔米特矩阵，其含义是矩阵中每一个第$i$行第$j$列的元素都与第$j$行第$i$列元素的共轭相等，可推知$Hermite$矩阵的共轭转置矩阵等于其本身。

​ 性质2 平稳离散时间随机过程的相关矩阵是$Toeplitz$矩阵。

​ 性质3 平稳离散时间随机过程的相关矩阵$R$是非负定的，且几乎总是正定的。

​ 性质4 将观测向量$u(n)$元素倒排，重新定义向量

​ $$u_B(n)=[u(n-M+1),u(n-M+2),...,u(n){]}^T$$

​ 这里，下标$B$表示对向量$u(n)$内各分量做反序排列，则向量$u_B(n)$的相关矩阵可以表示如下式

​ $$R_B=E[u_B(n)u_B^T(n)]$$

​ $$=\left[\begin{array}{cccc}r(0)& r^{\ast }(1)& ...& r^{\ast }(M-1)\\ r(1)& r(0)& ...& r^{\ast }(M-2)\\ ...& ...& ...& ...\\ r(M-1)& r(M-2)& ...& r(0)\end{array}\right]$$

​ 性质5 平稳离散时间随机过程的自相关矩阵$R$从$M$维扩展为$M+1$维，有如下递推关系

​ $$R_{M+1}=\left[\begin{array}{cc}r(0)& r^H\\ r& R_M\end{array}\right]$$

​ 或等价地，有

​ $$R_{M+1}=\left[\begin{array}{cc}{R}_M& r_B^{\ast }\\ r_B^T& r(0)\end{array}\right]$$

​ 式中 $r^H=\left[\begin{array}{cccc}r(1)& r(2)& ...& r(M)\end{array}\right]$，$r_B^T=\left[\begin{array}{cccc}r(-M)& r(-M+1)& ...& r(-1)\end{array}\right]$

1.3 自相关矩阵的特征值与特征向量的性质

​ 对平稳随机过程的自相关矩阵$R$进行特征值分解，设向量$q_1,q_2,...,q_M$分别是特征值${\lambda }_1,{\lambda }_2,...,{\lambda }_M$所对应的特征向量，即

​ $$R{q}_i={\lambda }_i{q}_i,i=1,...,M$$

​ 通过对自相关矩阵$R$进行特征值分解，可以得到随机过程$u(n)$的某些统计信息，这便是离散时间随机过程的特征值分析方法，是统计信号处理的基础。

​ 自相关矩阵$R$的特征值和特征向量的性质：

​ 性质1 特征值${\lambda }_1,{\lambda }_2,...,{\lambda }_M$都是实数，且是非负的

​ 性质2 对任意整数$k>0$，矩阵$R^k$的特征值为${\lambda }_1^k,{\lambda }_2^k,...,{\lambda }_M^k$

​ 性质3 若特征值${\lambda }_1,{\lambda }_2,...,{\lambda }_M$各不相同，则特征向量$q_1,q_2,...,q_M$相互正交

​ 性质4 若特征值${\lambda }_1,{\lambda }_2,...,{\lambda }_M$各不相同，$q_1,q_2,...,q_M$是相应的归一化特征向量，即

​ $$q_i^H{q}_j=\{\begin{array}{l}1,i=j\\ 0,i\ne j\end{array}$$

​ 定义矩阵$Q=\left[\begin{array}{c}{q}_1,q_2,...,q_M\end{array}\right]$，$\Lambda =diag[{\lambda }_1,{\lambda }_2,...,{\lambda }_M]$

​ 则矩阵$Q$是酉矩阵(unitary matrix)，且相关矩阵$R$可对角化为$Q^{H}RQ=$$\Lambda$

​ 性质5 特征值之和等于相关矩阵$R$的迹，即$tr(R)=Mr(0)={\sum }_{i=1}^M{\lambda }_i$

​ 性质6 $Karhunen-Loeve$展开：设零均值平稳随机过程$u(n)$构成的$M$维随机向量为$u(n)$，相应的相关矩阵为$R$，则向量$u(n)$可以表示为$R$的归一化特征向量$q_1,q_2,...,q_M$的线性组合，即

​ $$u(n)=\sum _{i=1}^M{c}_i{q}_i$$ $(\mathrm{1.3.1})$

​ 式中，展开式的系数$c_i$是由于内积$c_i=q_i^{H}u(n),i=1,2,...,M$定义的随机变量，且有

​ $$E[c_i]=0$$

​ $$E[c_i{c}_l^{\ast }]=\{\begin{array}{l}{\lambda }_i,i=1\\ 0,i\ne 1\end{array}$$

​ 式(1.3.1)称为$u(n)$的 $Karhunen-Loeve$展开式

二、MUSIC算法

​ 信号频率估计的多重信号分类(MUSIC，multiple signal classification)算法于1979年由R.O.Schmidt提出，该算法利用信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，估计信号频率。

​ 步骤1 根据$N$个观测样本值$x(0),x(1),...,x(N-1)$，估计自相关矩阵；

​ 步骤2 对自相关矩阵进行特征值分解，得到$M-K$个最小特征值对应的归一化特征向量，即得到噪声子空间的一组基向量，并构造噪声子空间矩阵$G=[u_{k+1},u_{k+2},...,u_M]\in C^{M\times (M-K)}$；

​ 步骤3 在$[-\pi ,\pi ]$内改变$w$，计算${\hat{P}}_{MUSIC(w)}=\frac{1}{a^H\hat{G}\widehat{G^H}a(w)}=\frac{1}{{\sum }_{i=K+1}^M|a^H(w)\widehat{u_i}{|}^2},w\in [-\pi ,\pi ]$，$P_{MUSIC}(w)$的峰值位置就是信号频率的估计值。

三、基于MUSIC算法的信号DOA估计方法

​ 将$K$个远场窄带信号从${\theta }_1,{\theta }_2,...,{\theta }_k$方向入射到$M$阵元的阵列时，阵列接收信号为

​ $$x(n)=As(n)+v(n)$$

​ 其中，$x(n)$为阵列接收数据向量，$A$为方向矩阵，$s(n)$是空间信号向量，$n(n)$是白噪声向量，将上式展开得到

​ $$\left[\begin{array}{c}{x}_0(n)\\ x_1(n)\\ ...\\ x_{M-1}(n)\end{array}\right]=\left[\begin{array}{cccc}1& 1& ...& 1\\ e^{-j{\varphi }_1}& e^{-j{\varphi }_2}& ...& e^{-j{\varphi }_k}\\ ...& ...& ...& ...\\ e^{-j(M-1){\varphi }_1}& e^{-j(M-1){\varphi }_2}& ...& e^{-j(M-1){\varphi }_k}\end{array}\right]\left[\begin{array}{c}{s}_1(n)\\ s_2(n)\\ ...\\ s_k(n)\end{array}\right]+\left[\begin{array}{c}{v}_0(n)\\ v_1(n)\\ ...\\ v_{M-1}(n)\end{array}\right]$$

​ 如果各信号源间相互统计独立，即

​ $$E[s_k(n)s_i^{\ast }(n)]=\{\begin{array}{l}{P}_k,k=i\\ 0,k\ne i\end{array}$$

​ 其中，$P_k$表示第$k$个信号的平均功率，则信号相关矩阵$P$是对角矩阵，即

​ $$P=E[s(n)s^k(n)]=diag[P_1,P_2,...,P_k]$$

​ 定义接收信号向量的空间相关矩阵为$R=E[x(n)x^H(n)]$，则

​ $$R=E[x(n)x^H(n)]=AP{A}^H+{\sigma }^{2}I$$

​ 其中，${\sigma }^2$为高斯白噪声的方差。为确保方向矩阵$A$的各列线性独立，应有$M>K$，即阵元数大于信号源数。因为矩阵$A$为Vandermonde（范德蒙德）矩阵，而且$P$为正定矩阵，则$AP{A}^H$的秩满足$rank(AP{A}^H)=K$，因此矩阵$AP{A}^H$存在$K$个正的特征值。

​ MUSIC算法DOA估计的原理和计算过程：

​ 对相关矩阵$R$进行特征值分解，并将这些特征值按单调非递增顺序排列，即

​ $${\lambda }_1\ge {\lambda }_2\ge ...\ge {\lambda }_K\ge {\lambda }_{K+1}={\lambda }_{K+2}=...={\lambda }_M={\sigma }^2$$

​ 这些特征值对应的归一化特征向量分别是$u_1,...,u_k,u_{k+1},...,u_M$，其中，$u_1,...,u_k$和$u_{k+1},...,u_M$分别张成信号子空间$E_s$和噪声子空间$E_N$，即

​ $$E_s=span[u_1,u_2,...,u_K]$$，$E_N=span[u_{K+1},u_{K+2},...,u_M]$

​ 定义矩阵$G=[u_{K+1},U_{K+2},...,u_M]\in C^{M\times (M-K)}$

​ 由于矩阵$A$是列满秩矩阵，$P$是满秩矩阵，可以证明$A^{H}G=0$

​ 或者等价地，有$G^{H}A=G^H\left[\begin{array}{cccc}a({\theta }_1)& a({\theta }_2)& ...& a({\theta }_K)\end{array}\right]=0$

​ 因此有$G^{H}a({\theta }_k)=0,k=1,2,...,K$

​ 其中，$a(\theta )$是阵列导向向量。

​ 实际应用中，根据$N$次快拍得到的接收数据$x(n),n=1,2,...,N$，用时间平均估计空间相关矩阵$R$，为

​ $$\hat{R}=\frac{1}{N}\sum _{n=1}^{N}x(n)x^H(n)$$

​ 可得MUSIC谱估计为$P_{MUSIC(\theta )}=\frac{1}{a^H(\theta )\hat{G}\widehat{G^H}a(\theta )},\theta \in (-\frac{\pi }{2},\frac{\pi }{2})$

​ 其中，矩阵$\hat{G}$是通过矩阵$\hat{R}$的特征值分解得到。MUSIC谱$P_{MUSIC}(\theta )$的$K$个峰值位置，就是信号波达方向${\theta }_k$的估计，其中$k=1,2,...,K$。由于MUSIC算法得到的谱并不是信号的空间功率谱，因此通常将式$P_{MUSIC(\theta )}=\frac{1}{a^H(\theta )\hat{G}\widehat{G^H}a(\theta )},\theta \in (-\frac{\pi }{2},\frac{\pi }{2})$称为伪谱。

​ 对于均匀矩形阵和均匀圆阵，其导向向量$a(\theta )$结构形式与均匀线阵一样，只要其方向矩阵$A$是列满秩的，仍可以利用MUSIC算法进行DOA估计。

（根据何子述《现代数字信号处理及其应用》整理而成。何子述 夏威 等.现代数字信号处理及其应用[M].北京：清华大学出版社, 2009年5月第1版：112~114）

​

​