子空间方法——MUSIC算法

多重信号分类MUSIC算法

1. MUSIC算法原理

MUSIC算法，叫做多信号分类算法 (Multiple Signal Classification)，是一种基于特征结构的高分辨率DOA算法。该算法利用了信号子空间和噪声子空间正交性的特点，构造噪声空间然后通过谱峰搜索来检测信号的波达方向。需要注意的是，该算法有一个前提，即各个入射信号之间互不相关，这样才能保证入射信号的协方差矩阵是满秩的。

考虑以下这个线阵模型：

令 $X(t)$ 是在第 $t$ 个快拍观察到的数据向量。在阵列信号处理和空间谱估计中, $\mathbf{X}(t)=$ ${\left[X_1(t),X_2(t),\cdots ,X_N(t)\right]}^{\mathrm{T}}$ 由 $N$ 个阵元 (天线或传感器) 的观测数据组成。在时域谱估计中，向量 $\mathbf{x}(t)=[x(t),x(t-1),\cdots ,x(t-n-1){]}^{\mathrm{T}}$ 由连续的 $n$ 个观察数据样本组成。
假定数据向量 $\mathbf{X}(t)$ 是 $r$ 个窄带信号入射到 $N$ 个阵元组成的阵列的观察数据向量或者是 $r$ 个不相干的复谐波的叠加，即
$$\mathbf{X}(t)=\sum _{i=1}^r{s}_i(t)\mathbf{a}\left({\omega }_i\right)+\mathbf{v}(t)=\mathbf{A}\mathbf{s}(t)+\mathbf{n}(t)$$式中， $\mathbf{A}=\left[\mathbf{a}\left({\omega }_1\right),\cdots ,\mathbf{a}\left({\omega }_r\right)\right]$ 为 $N\times r$ 阵列流型矩阵，$\mathbf{a}\left({\omega }_i\right)={\left[1,{\mathrm{e}}^{\mathrm{j}{\omega }_i},\dots ,{\mathrm{e}}^{\mathrm{j}(N-1){\omega }_i}\right]}^{\mathrm{T}}$ 为方向向量; $\mathbf{s}(t)={\left[s_1(t),\cdots ,s_r(t)\right]}^{\mathrm{T}}$ 为随机信号向量，其均值为零向量，协方差矩阵为 ${\mathbf{R}}_{\mathbf{s}}=\mathrm{E}\left\{\mathbf{s}(t){\mathbf{s}}^{\mathrm{H}}(t)\right\};$ 而 $\mathbf{v}(t)={\left[v_1(t),\cdots ,v_N(t)\right]}^{\mathrm{T}}$ 为加性噪声向量，其各个分量为高斯白噪声，它们具有零均值和相同的方差 ${\sigma }^2$ 。在谐波恢复中，参数 ${\omega }_i$ 为复谐波的频率 ; 在阵列信号处理中, ${\omega }_i$ 是一空间参数
$${\omega }_i=2\pi \frac{d}{\lambda }\sin {\theta }_i$$式中, $d$ 为相邻两个阵元之间的距离 (假定阵元等间距排列成一直线)，$\lambda$ 为波长, 且 ${\theta }_i$ 表 示第 $i$ 个窄带信号达到阵元的入射方向，简称波达方向。
MUSIC算法要解决的是根据 $N$ 个快拍的观测数据向量 $\mathbf{X}(t)(t=1,2,\cdots ,N)$ 估计 $r$ 个参数 ${\omega }_{i\circ }$ 这相当于对 $r$ 个混合信号进行分类，简称多重信号分类。
假定噪声向量 $\mathbf{v}(t)$ 与信号向量 $\mathbf{s}(t)$ 统计不相关，并令观测数据向量的协方差矩阵 ${\mathbf{R}}_{xx}=\mathrm{E}\left\{\mathbf{x}(t){\mathbf{x}}^{\mathrm{H}}(t)\right\}$ 的特征值分解为
$${\mathbf{R}}_{xx}={\mathbf{AR}}_{ss}{\mathbf{A}}^{\mathrm{H}}+{\sigma }^2\mathbf{I}=\mathbf{U}\mathbf{\Sigma }{\mathbf{U}}^{\mathrm{H}}=[{\mathbf{U}}_{\mathbf{S}},{\mathbf{U}}_{\mathbf{N}}]\left[\begin{array}{cc}{\mathbf{\Sigma }}_{\mathbf{S}}& \mathbf{O}\\ \mathbf{O}& {\mathbf{\Sigma }}_{\mathbf{N}}\end{array}\right]\left[\begin{array}{l}{\mathbf{U}}_S^{\mathrm{H}}\\ {\mathbf{U}}_N^{\mathrm{H}}\end{array}\right]$$式中, ${\mathbf{R}}_{ss}=\mathrm{E}\left\{\mathbf{s}(t){\mathbf{s}}^{\mathrm{H}}(t)\right\},$ 且 ${\mathbf{\Sigma }}_{\mathbf{S}}$ 包含了 $r$ 个大特征值, 它们比 ${\sigma }^2$ 明显大很多，${\mathbf{\Sigma }}_{\mathbf{N}}={\sigma }^2{I}_{n-r}$。
因此有：
$${\mathbf{R}}_{xx}{\mathbf{U}}_{\mathbf{N}}=[{\mathbf{U}}_{\mathbf{S}},{\mathbf{U}}_{\mathbf{N}}]\left[\begin{array}{cc}{\mathbf{\Sigma }}_{\mathbf{S}}& \mathbf{O}\\ \mathbf{O}& {\mathbf{\Sigma }}_{\mathbf{N}}\end{array}\right]\left[\begin{array}{l}{\mathbf{U}}_S^{\mathrm{H}}\\ {\mathbf{U}}_N^{\mathrm{H}}\end{array}\right]{\mathbf{U}}_{\mathbf{N}}=[{\mathbf{U}}_{\mathbf{S}},{\mathbf{U}}_{\mathbf{N}}]\left[\begin{array}{cc}{\mathbf{\Sigma }}_{\mathbf{S}}& \mathbf{O}\\ \mathbf{O}& {\mathbf{\Sigma }}_{\mathbf{N}}\end{array}\right]\left[\begin{array}{l}\mathbf{O}\\ \mathbf{I}\end{array}\right]={\sigma }^2{\mathbf{U}}_{\mathbf{N}}$$
又由 ${\mathbf{R}}_{xx}={\mathbf{AR}}_{ss}{\mathbf{A}}^{\mathrm{H}}+{\sigma }^2\mathbf{I}$ 有 ${\mathbf{R}}_{xx}{\mathbf{U}}_{\mathbf{N}}={\mathbf{AR}}_{ss}{\mathbf{A}}^{\mathrm{H}}{\mathbf{U}}_{\mathbf{N}}+{\sigma }^2{\mathbf{U}}_{\mathbf{N}},$联立上式可以得到
$${\mathbf{AR}}_{ss}{\mathbf{A}}^{\mathrm{H}}{\mathbf{U}}_{\mathbf{N}}=\mathbf{O}$$进而有
$${\mathbf{U}}_N^{\mathrm{H}}{\mathbf{AR}}_{ss}{\mathbf{A}}^{\mathrm{H}}{\mathbf{U}}_{\mathbf{N}}=\mathbf{O}$$
众所周知, $\mathbf{Q}$ 非奇异时 ${\mathbf{t}}^{\mathrm{H}}\mathbf{Q}\mathbf{t}=0,$ 当且仅当 $\mathbf{t}=0,$ 故上式成立的充分必要条件是
$${\mathbf{A}}^{\mathrm{H}}{\mathbf{U}}_{\mathbf{N}}=\mathbf{O}$$因为 $R_{ss}=\mathrm{E}\left\{s(t)s^{\mathrm{H}}(t)\right\}$ 非奇异。将 $\mathbf{A}=\left[\mathbf{a}\left({\omega }_1\right),\cdots ,\mathbf{a}\left({\omega }_p\right)\right]$ 代入上式即有
$${\mathbf{a}}^{\mathrm{H}}(\omega )\mathbf{G}=0^{\mathrm{T}},\omega ={\omega }_1,{\omega }_2,\cdots ,{\omega }_p$$显然, 当 $\omega \ne {\omega }_1,{\omega }_2,\cdots ,{\omega }_p$ 时, ${\mathbf{a}}^{\mathrm{H}}(\omega )\mathbf{G}\ne 0^{\mathrm{T}}$ 。
将上式改写成标量形式, 可以定义一种类似于功率谱的函数
$$\mathbf{P}(\omega )=\frac{1}{{\mathbf{a}}^{\mathrm{H}}(\omega ){\mathbf{U}}_N{\mathbf{U}}_N^{\mathrm{H}}\mathbf{a}(\omega )}$$由于噪声的影响，${\mathbf{a}}^{\mathrm{H}}(\omega ){\mathbf{U}}_{\mathbf{N}}{\mathbf{U}}_N^{\mathrm{H}}\mathbf{a}(\omega )\ne 0$，而是一个很小的数，因此此时的$\mathbf{P}(\omega )$是一个极大值。对上式取峰值的 $r$ 个 $\omega$ 值 ${\omega }_1,{\omega }_2,\cdots ,{\omega }_r$ 给出 $r$ 个信号的波达方向 ${\theta }_1,{\theta }_2,\cdots ,{\theta }_r$。
同时在实际应用中，接收阵接收到的数据有限长，我们可以通过这些接收数据估计入射信号的协方差矩阵，达到极大似然估计的目的。假设有M个快拍snapshot，则有
$${\hat{R}}_{xx}=\frac{1}{M}\sum _{i=1}^M{x}_i{x}_i^H$$然后对协方差矩阵${\hat{R}}_{xx}$进行特征值分解：
$$\begin{array}{rl}{\hat{R}}_{xx}& =U\Sigma U^H\end{array}$$特征值分解后按照特征值的大小(从大到小)对特征向量进行重排，接下来就是要确定噪声子空间。假设排序后的特征值${\sigma }_1^2⩾{\sigma }_2^2⩾\cdots ⩾{\sigma }_r^2⩾{\sigma }_{r+1}^2⩾\cdots ⩾{\sigma }_N^2$，接收信号总能量定义为$$P={\sigma }_1^2+{\sigma }_2^2+\cdots +{\sigma }_N^2$$根据能量准则，如果J是满足下式的最小整数：$$\frac{{\sigma }_1^2+{\sigma }_2^2+\cdots +{\sigma }_J^2}{P}⩾0.95或0.9$$那么 { σ J + 1 2 , σ J + 2 2 ⋯   , σ N 2 } \{\sigma_{J+1}^{2},\sigma_{J+2}^{2}\cdots,\sigma_{N}^{2}\} {σJ+12​,σJ+22​⋯,σN2​}对应的特征向量组成的空间$\widehat{U_N}$即为噪声子空间。

然后就可以用$\mathbf{P}(\omega )=\frac{1}{{\mathbf{a}}^{\mathrm{H}}(\omega ){\hat{\mathbf{U}}}_{\mathbf{N}}{\hat{\mathbf{U}}}_{\mathbf{N}}^{\mathrm{H}}\mathbf{a}(\omega )}$进行谱峰搜索，寻找r个入射信号的来波方向。

2. MATLAB仿真代码

% Code For Music Algorithm
% Author:痒羊羊
% Date：2020/10/28

clc; clear all; close all;
%% -------------------------initialization-------------------------
f = 500;                                        % frequency
c = 1500;                                       % speed sound
lambda = c/f;                                   % wavelength
d = lambda/2;                                   % array element spacing
M = 10;                                         % number of array elements
N = 100;                                        % number of snapshot
K = 6;                                          % number of sources
doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals

%% generate signal
dd = (0:M-1)'*d;                                % distance between array elements and reference element
A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
S = sqrt(2)\(randn(K,N)+1i*randn(K,N));         % array of random signal, K*N
X = A*S;                                        % received data without noise, M*N
X = awgn(X,10,'measured');                      % received data with SNR 10dB

%% calculate the covariance matrix of received data and do eigenvalue decomposition
Rxx = X*X'/N;                                   % covariance matrix
[U,V] = eig(Rxx);                               % eigenvalue decomposition
V = diag(V);                                    % vectorize eigenvalue matrix
[V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
U = U(:,idx);                                   % reset the eigenvector
P = sum(V);                                     % power of received data
P_cum = cumsum(V);                              % cumsum of V

%% define the noise space
J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
J = J(1);                                       % number of principal component
Un = U(:,J+1:end);

%% music for doa; seek the peek
theta = -90:0.1:90;                             % steer theta
doa_a = exp(-1i*2*pi*dd*sind(theta)/lambda);    % manifold array for seeking peak
music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
music = 10*log10(music/max(music));             % normalize the result and convert it to dB

%% plot
figure;
plot(theta, music, 'linewidth', 2);
title('Music Algorithm For Doa', 'fontsize', 16);
xlabel('Theta(°)', 'fontsize', 16);
ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
grid on;

运行结果：

可以看到，入射信号互不相关时，传统MUSIC算法能够以高分辨率检测出6个源的大致波达方向，分别是-29.7°，0°，19.8°，39.8°，60.4°，74.7°，但仍存在估计精度的问题，有很多改进的MUSIC算法可以改善。
需要注意，阵元数目为M的半波长均匀线阵的自由度为DOF为M-1，表示该线阵最多可以分辨源的数目为M-1。同时若存在相干源，MUSIC算法的效果就会不理想，基于相干源的MUSIC算法会在后续文章详写。

部分MUSIC算法原理参考张贤达先生的《矩阵分析与应用》。
欢迎转载，表明出处。