一文读懂MUSIC算法DOA估计的数学原理并仿真

一文读懂MUSIC算法DOA估计的数学原理并仿真

---

前言

MUSIC（Multiple Signal Classification）算法于1979年由R.O.Schmidt提出，是阵列信号处理中广泛应用的经典DOA（Direction of Arrival）估计算法，凭借其超分辨的估计性能受到广泛关注。本文将从数学公式推导的角度出发系统阐述MUSIC算法的基本原理，包括为什么噪声子空间与信号子空间是正交的以及如何根据正交性来构建空间谱函数实现对信号方向的高精度估计。

---

提示：以下是本篇文章正文内容，转载请附上链接！

一、DOA估计基本原理

当$K$ 个远场窄带信号$s_1(n),s_2(n),...,s_K(n)$从${\theta }_1,{\theta }_2,\cdot \cdot \cdot ,{\theta }_K$ 方向人射到$M$ 阵元的阵列时，阵列接收信号为
$$\mathbf{x}(n)=\mathbf{As}(n)+\mathbf{v}(n)$$ 其中$,\mathbf{x}(n)\in {\mathbb{C}}^{M\times 1}$为阵列接收数据向量$,\mathbf{A}\in {\mathbb{C}}^{M\times K}$为方向矩阵$,\mathbf{s}(n)\in {\mathbb{C}}^{K\times 1}$为空间信号向量，$\mathbf{v}(n)\in {\mathbb{C}}^{M\times 1}$是白噪声向量。对于均匀线阵，展开即有
$$\mathbf{x}(n)=[x_0(n)x_1(n)...x_{M-1}(n){]}^{\mathrm{T}}$$ $$\mathbf{A}=[\mathbf{a}({\theta }_1)\mathbf{a}({\theta }_2)...\mathbf{a}({\theta }_K)]$$ $$\mathbf{s}(n)=[s_1(n)s_2(n)...s_K(n){]}^{\mathrm{T}}$$ $$\mathbf{v}(n)=[v_0(n)v_1(n)...v_{M-1}(n){]}^{\mathrm{T}}$$如果各信号源间相互统计独立，即 E ⁡ { s k ( n ) s i ∗ ( n ) } = { σ k 2 , k = i 0 , k ≠ i \operatorname{E}\{s_k(n)s_i^*(n)\}= \begin{cases} \sigma_k^2, & k=i \\ 0, & k\neq i & & \end{cases} E{sk​(n)si∗​(n)}={σk2​,0,​k=ik=i​​​其中，${\sigma }_k^2$表示第$k$个信号的平均功率，则信号相关矩阵$\mathbf{P}$是对角矩阵，即
$$\mathbf{P}=\mathrm{E}\left[\mathbf{s}\left(n\right){\mathbf{s}}^{\mathrm{H}}\left(n\right)\right]=\mathrm{diag}\left\{{\sigma }_1^2,{\sigma }_2^2,\cdot \cdot \cdot ,{\sigma }_k^2\right\}$$因为是零均值、方差为${\sigma }^2$的白噪声，所以有
E ⁡ { v ( n ) v H ( n ) } = diag ⁡ { σ 2 , ⋅ ⋅ ⋅ , σ 2 } = σ 2 I \operatorname{E}\{\boldsymbol{v}(n)\boldsymbol{v}^\mathrm{H}(n)\}=\operatorname{diag}\{{\sigma}^2,\cdotp\cdotp\cdotp,{\sigma}^2\}={\sigma}^2\boldsymbol{I} E{v(n)vH(n)}=diag{σ2,⋅⋅⋅,σ2}=σ2I定义接收信号向量的空间相关矩阵为 R = E { x ( n ) x H ( n ) } = E { [ A s ( n ) + v ( n ) ] [ s H ( n ) A H + v H ( n ) ] } = A P A H + σ 2 I \begin{aligned} \mathbf{R}&=\mathbb{E}\left\{\mathbf{x}(n)\mathbf{x}^{\mathrm{H}}(n)\right\}\\&=\mathbb{E}\{\begin{bmatrix} \boldsymbol{A}\boldsymbol{s}(n)+\boldsymbol{v}(n)\end{bmatrix}\begin{bmatrix} \boldsymbol{s}^\mathrm{H}(n)\boldsymbol{A}^\mathrm{H}+\boldsymbol{v}^\mathrm{H}(n) \end{bmatrix}\}\\&=\boldsymbol{A}\boldsymbol{P}\boldsymbol{A}^\mathrm{H}+\sigma^2\boldsymbol{I} \end{aligned} R​=E{x(n)xH(n)}=E{[As(n)+v(n)​][sH(n)AH+vH(n)​]}=APAH+σ2I​当然，$\mathbf{R}$也可以表示为$$\begin{array}{rl}\mathbf{R}& =[\mathbf{a}({\theta }_1)\mathbf{a}({\theta }_2)...\mathbf{a}({\theta }_K)]\mathrm{diag}\left\{{\sigma }_1^2,{\sigma }_2^2,\cdot \cdot \cdot ,{\sigma }_k^2\right\}[\mathbf{a}({\theta }_1)\mathbf{a}({\theta }_2)...\mathbf{a}({\theta }_K){]}^{\mathrm{H}}+{\sigma }^2\mathbf{I}\\ & =\sum _{k=1}^K{\sigma }_k^2\mathbf{a}({\theta }_k)\mathbf{a}({\theta }_k{)}^{\mathrm{H}}+{\sigma }^2\mathbf{I}\end{array}$$为确保方向矩阵$\mathbf{A}$的各列线性独立，应有$M>K$，即阵元数大于信号源数。
方向矩阵$\mathbf{A}$是列满秩的，因此，${\mathbf{A}}^{\mathrm{H}}$是行满秩矩阵，即$$\mathrm{rank}(\mathbf{A})=\mathrm{rank}({\mathbf{A}}^{\mathrm{H}})=K$$又因为$\mathbf{P}$是秩为$K$的满秩对角阵。因此，矩阵乘积$\mathbf{A}\mathbf{P}$是对矩阵$\mathbf{A}$作满秩变换，变换后的秩保持不变，即$$\mathrm{rank}(\mathbf{AP})=K$$同样，${\mathbf{APA}}^{\mathrm{H}}$也 为 对 矩 阵 ${\mathbf{A}}^{\mathrm{H}}$作 满 秩 变 换 , 即
$$\mathrm{rank}({\mathbf{APA}}^{\mathrm{H}})=K$$这里用到线性代数的一个知识是左乘列满秩，右乘行满秩，不改变矩阵的秩。
因此，矩阵 ${\mathbf{APA}}^{\mathrm{H}}$ 共有 $K$ 个非零特征值。对 ${\mathbf{APA}}^{\mathrm{H}}$ 进行特征值分解，设 ${\tilde{\lambda }}_1,{\tilde{\lambda }}_2,\cdot \cdot \cdot ,{\tilde{\lambda }}_M$ 为其特征值；${\mathbf{u}}_1,{\mathbf{u}}_2,\cdots ,{\mathbf{u}}_M$为对应的正交归一化特征向量。不妨将其$K$ 个非零特征值设为${\tilde{\lambda }}_1,{\tilde{\lambda }}_2,\cdot \cdot \cdot ,{\tilde{\lambda }}_K\ne 0；$其余特征值 ${\tilde{\lambda }}_{K+1}={\tilde{\lambda }}_{K+2}=\cdot \cdot \cdot {\tilde{\lambda }}_M=0$，所以有
$$\begin{array}{rl}& ({\mathbf{APA}}^{\mathrm{H}}){\mathbf{u}}_i={\tilde{\lambda }}_i{\mathbf{u}}_i,i=1,2,\cdot \cdot \cdot ,K\\ & ({\mathbf{APA}}^{\mathrm{H}}){\mathbf{u}}_i={\tilde{\lambda }}_i{\mathbf{u}}_i=0,i=K+1,K+2,\cdot \cdot \cdot ,M\end{array}$$对上面两个式子右乘${\mathbf{u}}_i^{\mathrm{H}}$,有
$$({\mathbf{APA}}^{\mathrm{H}}){\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}={\tilde{\lambda }}_i{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}},i=1,\cdot \cdot \cdot ,M$$上式中分别取$i=1,\cdot \cdot \cdot ,M$，可得$M$个等式，将各等式两边分别相加，得
$$({\mathbf{APA}}^{\mathrm{H}})\sum _{i=1}^M{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}=\sum _{i=1}^M{\tilde{\lambda }}_i{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}$$考虑到${\mathbf{u}}_1,{\mathbf{u}}_2,\cdots ,{\mathbf{u}}_M$是正交的归一化特征向量，所以有
$$\sum _{i=1}^M{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}=\mathbf{I}$$所以有
$${\mathbf{APA}}^{\mathrm{H}}=\sum _{i=1}^M{\tilde{\lambda }}_i{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}=\sum _{i=1}^K{\tilde{\lambda }}_i{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}$$因此，自相关矩阵$\mathbf{R}$可以表示为$$\begin{array}{rl}\mathbf{R}& =\sum _{i=1}^K{\tilde{\lambda }}_i{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}+{\sigma }^2\sum _{i=1}^M{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}\\ & =\sum _{i=1}^K({\tilde{\lambda }}_i+{\sigma }^2){\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}+{\sigma }^2\sum _{i=K+1}^M{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}\\ & =\sum _{i=1}^M{\lambda }_i{\mathbf{u}}_i{\mathbf{u}}_i^{\mathrm{H}}\end{array}$$其中$$\begin{gathered} {\lambda }_i={\tilde{\lambda }}_i+{\sigma }^2,i=1,\cdot \cdot \cdot ,K \\ {\lambda }_i={\sigma }^2,i=K+1,\cdot \cdot \cdot ,M \end{gathered}$$ $\mathbf{R}$的$M$个特征值中仅有$K$个特征值${\lambda }_1,\cdot \cdot \cdot ,{\lambda }_K$与信号有关，其余$M-K$个特征值${\lambda }_{K+1},\cdot \cdot \cdot ,{\lambda }_M$ 仅与噪声有关。根据以上事实，有以下重要定义：
（1）信号子空间${\mathbf{E}}_{\mathrm{s}}($signal subspace)：${\mathbf{E}}_{\mathrm{s}}$是由${\lambda }_1,\cdot \cdot \cdot ,{\lambda }_K$对应的特征向量${\mathbf{u}}_1,\cdot \cdot \cdot ,{\mathbf{u}}_K$生成的子空间，记为
E S = s p a n { u 1 , u 2 , ⋅ ⋅ ⋅ , u K } \boldsymbol{E}_\mathrm{S}=\mathrm{span}\{\mathbf{u}_1,\mathbf{u}_2,\cdotp\cdotp\cdotp,\mathbf{u}_K\} ES​=span{u1​,u2​,⋅⋅⋅,uK​}（2）噪声子空间${\mathbf{E}}_N$(noise subspace)：${\mathbf{E}}_N$是由${\lambda }_{K+1},\cdot \cdot \cdot ,{\lambda }_M$对应的特征向量${\mathbf{u}}_{K+1},\cdot \cdot \cdot ,{\mathbf{u}}_M$ 生成的子空间，记为
E N = span ⁡ { u K + 1 , u K + 2 , ⋅ ⋅ ⋅ , u M } \boldsymbol{E}_N=\operatorname{span}\{\boldsymbol{u}_{K+1},\boldsymbol{u}_{K+2},\cdotp\cdotp\cdotp,\boldsymbol{u}_M\} EN​=span{uK+1​,uK+2​,⋅⋅⋅,uM​} ${\mathbf{E}}_{\mathrm{S}}$ 既与信号有关，又与噪声有关；而${\mathbf{E}}_{\mathrm{N}}$ 仅与噪声有关。
由前面的推导有：
$$({\mathbf{APA}}^{\mathrm{H}}){\mathbf{u}}_i={\tilde{\lambda }}_i{\mathbf{u}}_i=0,i=K+1,K+2,\cdot \cdot \cdot ,M$$对上式两边同时左乘 ${\mathbf{A}}^{\mathrm{H}}$，得
$${\mathbf{A}}^{\mathrm{H}}{\mathbf{APA}}^{\mathrm{H}}{\mathbf{u}}_i=0,i=K+1,\cdot \cdot \cdot ,M$$由前面的推导易知：rank(${\mathbf{A}}^{\mathrm{H}}\mathbf{A})=K$，即 ${\mathbf{A}}^{\mathrm{H}}\mathbf{A}$ 可逆。上式等号两边同时左乘$({\mathbf{A}}^{\mathrm{H}}\mathbf{A}{)}^{-1}$，有
$$({\mathbf{A}}^{\mathrm{H}}\mathbf{A}{)}^{-1}{\mathbf{A}}^{\mathrm{H}}{\mathbf{APA}}^{\mathrm{H}}{\mathbf{u}}_i={\mathbf{PA}}^{\mathrm{H}}{\mathbf{u}}_i=0,i=K+1,\cdot \cdot \cdot ,M$$又因为矩阵$\mathbf{P}$为正定的对角矩阵，上式两边可再同时左乘${\mathbf{P}}^{-1}$,有
$${\mathbf{P}}^{-1}{\mathbf{PA}}^{\mathrm{H}}{\mathbf{u}}_i={\mathbf{A}}^{\mathrm{H}}{\mathbf{u}}_i=0,i=K+1,\cdot \cdot \cdot ,M$$定义矩阵$\mathbf{G}$是由噪声子空间的向量构成的矩阵，即
$$\mathbf{G}=\left[\begin{array}{cccc}{\mathbf{u}}_{K+1}& {\mathbf{u}}_{K+2}& \cdots & {\mathbf{u}}_M\end{array}\right]\in {\mathbb{C}}^{M\times (M-K)}$$所以有
$${\mathbf{A}}^{\mathrm{H}}\mathbf{G}=0$$或者等价地，有
$${\mathbf{G}}^{\mathrm{H}}\mathbf{A}={\mathbf{G}}^{\mathrm{H}}\left[\mathbf{a}\left({\theta }_1\right)\mathbf{a}\left({\theta }_2\right)\cdots \mathbf{a}\left({\theta }_K\right)\right]=0$$即每一个方向矢量与噪声子空间都是正交的。
实际应用中，根据$N$次快拍得到的接收数据 $\mathbf{x}(n),n=1,2,\cdot \cdot \cdot ,N$，用时间平均估计空间相关矩阵 $\hat{\mathbf{R}}$ 为
$$\hat{\mathbf{R}}=\frac{1}{N}\sum _{n=1}^N\mathbf{x}(n){\mathbf{x}}^{\mathrm{H}}(n)$$在实际工程中，由于用相关矩阵的估计$\hat{\mathbf{R}}$代替 $\mathbf{R}$ 进行特征分解，因此，在给定的角度${\theta }_k$，信号方向向量$\mathbf{a}({\theta }_k)$与噪声子空间并不严格地满足正交条件方程式。所以左边并不严格等于零，而是一个很小的值。于是，可以构造如下的空间谱扫描函数：
$$\begin{array}{rl}{P}_{\mathrm{MUSIC}}(\theta )& =\frac{1}{{\mathbf{a}}^{\mathrm{H}}(\theta ){\mathbf{GG}}^{\mathrm{H}}\mathbf{a}(\theta )}\\ & =\frac{1}{{\sum }_{i=K+1}^M\mid {\mathbf{a}}^{\mathrm{H}}(\theta ){\mathbf{u}}_i{\mid }^2},\theta \in \left[\begin{array}{c}-\pi /2,\pi /2\end{array}\right]\end{array}$$信号到达角度的估计可以由函数$P_{\mathrm{MUSIC}}(\theta )$的$K$个峰值位置确定。
谱函数 $P_{\mathrm{MUSIC}}(\omega )$的峰值位置反映了信号的到达角度，但它并不是信号的功率谱，通常将 $P_{\mathrm{MUSIC}}(\theta )$称为伪谱(pseudo spectrum),或 MUSIC 谱。

二、MATLAB仿真

考虑一个8阵元的均匀线阵，采样快拍数为1024，信号入射角度分别为-40°，-10°，0°，40°，用MUSIC算法估计信号到达角度如下图所示

需要注意的点：
1、对矩阵$\hat{\mathbf{R}}$分解得到的大的$K$个特征值并不是信号的功率，较小的$M-K$个特征值一定是噪声的功率。
2、信噪比较低时，MUSIC算法失效，即MUSIC算法适用于信噪比较大的场景。
3、当信号相干时，MUSIC算法也会失效。

---

总结

以上就是本文所有的内容，本文详细介绍了谱估计公式推导的数学过程，本文展示了MUSIC算法如何利用噪声子空间和信号子空间的正交性实现信号的到达方向估计，同时还给出了MATLAB仿真的结果。