【学习笔记】【DOA子空间算法】5 SS-MUSIC 算法

5 SS-MUSIC 算法

  在学习 SS-MUSIC 算法之前需要理解相干信号的概念。首先定义两个平稳信号 $s_i(t)$ 和 $s_k(t)$ 的相关系数 ${\rho }_{ik}$ 如下：
$$\begin{array}{c}{\rho }_{ik}=\frac{\mathrm{E}[s_i(t)s_k^{\ast }(t)]}{\sqrt{\mathrm{E}[|s_i(t){|}^2]\mathrm{E}[|s_k(t){|}^2]}}\end{array}$$
若 ${\rho }_{ik}=0$，则称两个信号相互独立；若 $0\le |{\rho }_{ik}|\le 1$，则称两个信号相关；若 $|{\rho }_{ik}|=1$，则称两个信号相干。
  相干信号源 $s_i(t)$ 和 $s_k(t)$ 之间只差一个复常数，即 $s_i(t)=\alpha \ast s_k(t)$，其中 $\alpha$ 为复常数。

5.1 算法原理

  传统 MUSIC 算法通常假设信号间互不相关，如此便能保证信号协方差矩阵 ${\mathbf{R}}_S$ 是非奇异对角阵，从而推导出方向矢量矩阵 $\mathbf{A}$ 和噪声子空间矩阵 ${\mathbf{U}}_N$ 是正交的，最后再遍历 { a ( θ ) : θ ∈ [ − 9 0 ∘ , 9 0 ∘ ] } \{\mathbf{a}(\theta):\theta \in [-90^{\circ},90^{\circ}]\} {a(θ):θ∈[−90∘,90∘]} 找出和 ${\mathbf{U}}_N$ 存在正交关系的方向矢量。然而若存在相干信号，则原信号协方差矩阵 ${\mathbf{R}}_S\in {\mathbb{C}}^{K\times K}$ 是奇异的。
  假设 $K$ 个信号中存在两个信号是相干的，即 $s_2(t)=\alpha s_1(t)$，其中 $\alpha$ 为复常数，则之前所曾提及的信号模型可以被改写如下：
$$\begin{array}{c}\mathbf{x}(t)=\mathbf{As}(t)+\mathbf{n}(t)\in {\mathbb{C}}^{M\times 1}\end{array}$$
其中 $\mathbf{A}\in {\mathbb{C}}^{M\times (K-1)}$：
$$\begin{array}{c}\mathbf{A}=[\mathbf{a}({\theta }_1)+\alpha \mathbf{a}({\theta }_2),\mathbf{a}({\theta }_3),\cdots ,\mathbf{a}({\theta }_K)]\end{array}$$
和 $\mathbf{s}(t)\in {\mathbb{C}}^{(K-1)\times 1}$：
$$\begin{array}{c}\mathbf{s}(t)=[s_1(t),s_3(t),\cdots ,s_K(t)]\end{array}$$
  由前面的讨论不难看出，新的信号协方差矩阵 ${\mathbf{R}}_S\in {\mathbb{C}}^{(K-1)\times (K-1)}$ 仍然是一个非奇异矩阵，这意味着若此时对协方差矩阵 $\mathbf{R}$ 进行特征值分解，只能获得 $M-(K-1)$ 个较小特征值；并且从 $\mathbf{A}$ 的式子可以看出，$\mathbf{A}$ 的第一列为 $\mathbf{a}({\theta }_1)+\alpha \mathbf{a}({\theta }_2)$，该方向矢量无法从 { a ( θ ) : θ ∈ [ − 9 0 ∘ , 9 0 ∘ ] } \{\mathbf{a}(\theta):\theta \in [-90^{\circ},90^{\circ}]\} {a(θ):θ∈[−90∘,90∘]} 中获得，即该方向矢量不在原信号子空间中。
  对于相干信号的 DOA 估计，相干信号会导致信源协方差矩阵的秩亏缺，有部分能量就会散发到噪声子空间，使得传统 MUSIC 算法不能对其进行有效估计；因而去相干的思路就是把信号协方差矩阵的秩恢复到信号数，从而正确估计信号源方向。
  空间平滑 MUSIC 算法是一种降维处理，牺牲有效阵列孔径来实现信号源的去相干；该方法的降维体现在其将原阵列拆分成很多个子阵列，通过对子阵列的协方差矩阵作平均，从而重构接收数据协方差矩阵。

5.1.1 前向空间平滑 FSS-MUSIC 算法

  假设将 $M$ 个阵元均匀分成 $L$ 个相互重叠的子阵列，每个子阵列中的阵元数为 $P=M-L+1$，如取第 $1\sim P$ 个阵元作为第一个子阵列，取第 $2\sim (P+1)$ 个阵元作为第二个子阵列，以此类推。可以定义第 $l$ 个子阵列的接收数据向量 ${\mathbf{x}}_l(t)\in {\mathbb{C}}^{P\times 1}$ 为（为方便书写，忽略噪声部分）：
$$\begin{array}{c}\begin{array}{rl}{\mathbf{x}}_l(t)& =[x_l(t),x_{l+1}(t),\cdots ,x_{l+P-1}(t){]}^T\\ & ={\mathbf{A}}_P{\mathbf{D}}^{l-1}\mathbf{s}(t)\end{array}\end{array}$$
其中 ${\mathbf{A}}_P\in {\mathbb{C}}^{P\times K}$ 由 $\mathbf{A}$ 的前 $P$ 行组成，$x_i(t)$ 表示 $\mathbf{x}(t)$ 的第 $i$ 个元素，$\mathbf{D}\in {\mathbb{C}}^{K\times K}$ 表示如下：
D = diag ⁡ { e − j 2 π d sin ⁡ θ 1 / λ , ⋯   , e − j 2 π d sin ⁡ θ K / λ } \begin{equation*} \mathbf{D}=\operatorname{diag}\{e^{-j2\pi d\sin\theta_1/\lambda},\cdots,e^{-j2\pi d\sin\theta_K/\lambda}\} \end{equation*} D=diag{e−j2πdsinθ1​/λ,⋯,e−j2πdsinθK​/λ}​
  由此可以得到第 $l$ 个子阵的协方差矩阵 ${\mathbf{R}}_l\in {\mathbb{C}}^{P\times P}$：
$$\begin{array}{c}{\mathbf{R}}_l={\mathbf{A}}_P{\mathbf{D}}^{l-1}{\mathbf{R}}_S({\mathbf{D}}^{l-1}{)}^H{\mathbf{A}}_P^H\end{array}$$
  因此，通过前向空间平滑重构后的协方差矩阵 ${\mathbf{R}}_{\mathrm{FSS}}$ 可以通过对各个子阵的协方差矩阵求平均得到：
$$\begin{array}{c}\begin{array}{rl}{\mathbf{R}}_{\mathrm{FSS}}& =\frac{1}{L}\sum _{l=1}^L{\mathbf{R}}_l\\ & ={\mathbf{A}}_P\left(\frac{1}{L}\sum _{l=1}^L{\mathbf{D}}^{l-1}{\mathbf{R}}_S({\mathbf{D}}^{l-1}{)}^H\right){\mathbf{A}}_P^H\\ & ={\mathbf{A}}_P{\overline{\mathbf{R}}}_S{\mathbf{A}}_P^H\end{array}\end{array}$$
此时只需要证明 ${\overline{\mathbf{R}}}_S\in {\mathbb{C}}^{K\times K}$ 满秩即可。注意，要令 ${\overline{\mathbf{R}}}_S$ 满秩，需要满足 $P\ge K$ 和 $L\ge K$。
  首先假设信号模型完全相干，即 { s k ( t ) = α k s 0 ( t ) : k = 1 , ⋯   , K } \{s_k(t) = \alpha_k s_0(t):k = 1,\cdots,K\} {sk​(t)=αk​s0​(t):k=1,⋯,K}，则有 $\mathbf{S}=\mathbf{\alpha }{\mathbf{s}}_0^T$，其中 $\mathbf{\alpha }=[{\alpha }_1,\cdots ,{\alpha }_K{]}^T$；再根据 $\mathrm{E}[|s_0(t){|}^2]=1$，可以得到：
$${\mathbf{R}}_S=\mathbf{S}{\mathbf{S}}^H=\mathbf{\alpha }{\mathbf{\alpha }}^H$$
  接下来将 ${\overline{\mathbf{R}}}_S$ 进一步化简：
$$\begin{array}{c}\begin{array}{rl}{\overline{\mathbf{R}}}_S& =\left[\begin{array}{cccc}\mathbf{I}& \mathbf{D}& \cdots & {\mathbf{D}}^{L-1}\end{array}\right]\left[\begin{array}{ccc}\frac{1}{L}{\mathbf{R}}_S& & \\ & \ddots & \\ & & \frac{1}{L}{\mathbf{R}}_S\end{array}\right]\left[\begin{array}{c}\mathbf{I}\\ {\mathbf{D}}^{-1}\\ ⋮\\ {\mathbf{D}}^{-L+1}\end{array}\right]\\ & =\frac{1}{L}\mathbf{G}{\mathbf{G}}^H\end{array}\end{array}$$
其中 $\mathbf{G}\in {\mathbb{C}}^{K\times L}$：
G = [ α , D α , ⋯   , D L − 1 α ] = diag ⁡ { α 1 , ⋯   , α K } ∗ A L T \begin{equation*} \begin{aligned} \mathbf{G} &= [\boldsymbol{\alpha}, \mathbf{D}\boldsymbol{\alpha},\cdots,\mathbf{D}^{L-1}\boldsymbol{\alpha}] \\ &=\operatorname{diag}\{\alpha_1,\cdots,\alpha_K\}*\mathbf{A}_L^T \end{aligned} \end{equation*} G​=[α,Dα,⋯,DL−1α]=diag{α1​,⋯,αK​}∗ALT​​​
  根据已知前提 $L\ge K$，则 $\mathrm{rank}({\mathbf{A}}_L^T)=K$，可得 $\mathrm{rank}({\overline{\mathbf{R}}}_S)=\mathrm{rank}(\mathbf{G})=K$，即 ${\overline{\mathbf{R}}}_S$ 满秩。
  由于子阵列的阵元数 $P=M-L+1$，则最多可以检测 $M-L$ 个相干信号，即 $K\le M-L$；同时又要满足 $K\le L$，因此最大情况下 $M-L=L$，表明 FSS-MUSIC 算法可以检测 $K=\frac{M}{2}$ 个相干信号。

5.1.2 后向空间平滑 BSS-MUSIC 算法

  后向空间平滑算法是对后向子阵列的共轭接收数据协方差矩阵进行平滑。同样假设将 $M$ 个阵元均匀分成 $L$ 个相互重叠的子阵列，每个子阵列中的阵元数为 $P=M-L+1$，但是和 FSS-MUSIC 算法不同，后向空间平滑算法取第 $(M-P+1)\sim M$ 个阵元的共轭信号作为第一个子阵列信号，取第 $(M-P)\sim (M-1)$ 个阵元的共轭信号作为第二个子阵列信号。
  容易得到后向空间平滑重构后的协方差矩阵 ${\mathbf{R}}_{\mathrm{BSS}}$ 与 ${\mathbf{R}}_{\mathrm{FSS}}$ 的关系：
$$\begin{array}{c}{\mathbf{R}}_{\mathrm{BSS}}=\mathbf{J}{\mathbf{R}}_{\mathrm{FSS}}^{\ast }\mathbf{J}\end{array}$$
其中 $\mathbf{J}$ 为翻转对角阵，它又被称为选择矩阵，不难看出这里就是将 ${\mathbf{R}}_{\mathrm{FSS}}$ 的共轭矩阵进行了对角翻转。同理 BSS-MUSIC 算法可以检测 $K=\frac{M}{2}$ 个相干信号。

5.1.3 前后向空间平滑 FBSS-MUSIC 算法

  前后向空间平滑算法将前向平滑、后向平滑的相关矩阵累加：
$$\begin{array}{c}{\mathbf{R}}_{\mathrm{FBSS}}=\frac{1}{2}({\mathbf{R}}_{\mathrm{FSS}}+{\mathbf{R}}_{\mathrm{BSS}})\end{array}$$
  不难看出，前后向空间平滑算法中利用了 $2L$ 个子阵列，即 $L$ 个前向子阵列和 $L$ 个后向共轭子阵列，因此可以将后向子阵列视为前向子阵列的镜像对称子阵列；但是要注意的是，FBSS-MUSIC 中 ${\mathbf{R}}_{\mathrm{FSS}}$ 和 ${\mathbf{R}}_{\mathrm{BSS}}$ 所对应的前向子阵列和后向子阵列必须中心对称。
  由于子阵列的阵元数 $P=M-L+1$，则最多可以检测 $M-L$ 个相干信号，即 $K\le M-L$；而前后向空间平滑算法进行了 $2L$ 次平滑，因而要满足 $K\le 2L$，因此最大情况下 $M-L=2L$，表明 FBSS-MUSIC 算法可以检测 $K=\frac{2M}{3}$ 个相干信号，所以说 FBSS-MUSIC 算法相比 FSS-MUSIC 算法更大地提高了阵列孔径。

5.2 算法步骤

  SS-MUSIC 算法步骤如下（输入为阵列接收矩阵 $\mathbf{X}$）：

1. 计算协方差矩阵 $\mathbf{R}=\frac{1}{T}\mathbf{X}{\mathbf{X}}^H$。
2. 采取相应的策略（FSS、BSS 或 FBSS）对 $\mathbf{R}$ 进行分块，从而累加求均值得到新的协方差矩阵 ${\mathbf{R}}_{\mathrm{SS}}\in {\mathbb{C}}^{P\times P}$，具体看代码实现。
3. 对 ${\mathbf{R}}_{\mathrm{SS}}$ 进行特征值分解，并对特征值进行排序，然后取得 $M-K$ 个较小特征值对应的特征向量来组成噪声子空间 ${\mathbf{U}}_N$。
4. 以下式遍历 $\theta =0^{\circ },1^{\circ }\cdots ,180^{\circ }$：
   $$\begin{array}{c}P(\theta )=\frac{1}{{\mathbf{a}}_P^H(\theta ){\mathbf{U}}_N{\mathbf{U}}_N^H{\mathbf{a}}_P(\theta )}\end{array}$$
   其中 ${\mathbf{a}}_P(\theta )$ 由 $\mathbf{a}(\theta )$ 的前 $P$ 行组成。此时得到一组 $P(\theta )$，$K$ 个最大值对应的 $\theta$ 就是需要返回的结果。

5.3 代码实现

% ss_music.m
clear;
clc;
close all;

%% 参数设定
c = 3e8;                                              % 光速
fc = 500e6;                                           % 载波频率
lambda = c/fc;                                        % 波长
d = lambda/2;                                         % 阵元间距，可设 2*d = lambda
twpi = 2.0*pi;                                        % 2pi
derad = pi/180;                                       % 角度转弧度
theta = [-30,-10,30,60]*derad;                        % 待估计角度
idx = 1:1:8; idx = idx';                              % 阵元位置索引

M = length(idx);                                      % 阵元数
K = length(theta);                                    % 信源数
T = 512;                                              % 快拍数
SNR = 10;                                             % 信噪比

%% 信号模型建立
coef = randn(K,1) + 1j*randn(K,1);
sig = randn(1,T) + 1j*randn(1,T);                     % 构造K个相干的信号
S = coef*sig;                                         % 复信号矩阵S，维度为K*T
% A = exp(-1j*twpi*d*idx*sin(theta)/lambda);            % 方向矢量矩阵A，维度为M*K
A = exp(-1j*pi*idx*sin(theta));                       % 2d = lambda，直接忽略不写
X = A*S;                                              % 接收矩阵X，维度为M*T
X = awgn(X,SNR,'measured');                           % 添加噪声

%% SS-MUSIC 算法
% 计算协方差矩阵
R = X*X'/T;

% 重构协方差矩阵使得秩恢复，注意P和L均需要大于或等于K
P = M-3;                                              % 单个子阵列中的阵元数，自行设置
L = M-P+1;                                            % 子阵列数目

% 方法一，前向空间平滑 FSS
Rf = zeros(P, P);
for i = 1:L
    Rf = Rf + R(i:i+P-1, i:i+P-1);
end

% 方法二，后向空间平滑 BSS
Rb = zeros(P, P);
J = fliplr(eye(M));                                   % 翻转对角阵
Rj = J*(conj(R))*J;
for i = 1:L
    Rb = Rb + Rj(i:i+P-1, i:i+P-1);
end

% 方法三，前后向空间平滑 FBSS
Rfb = zeros(P, P);
Rjj = (R+Rj)/2;                                       % 利用前面求的Rj
for i = 1:L
    Rfb = Rfb + Rjj(i:i+P-1, i:i+P-1);
end

% 重新构造协方差之后需要重新设定一些参数，如果想查看传统MUSIC算法对相干信号模型的结果，只需注释下面三行即可
R = Rfb/L;                                            % 该行用作选取具体的重构方法，假若选用方法一，则R = Rf/L
M = P;                                                % 此时需要设置构造后的阵元数为P
idx = 1:1:M; idx = idx';                              % 并且设置构造后的索引

% 特征值分解并取得噪声子空间
[U,D] = eig(R);                                       % 特征值分解
D = diag(D)';                                         % 提取特征值
[D,I] = sort(D);                                      % 将特征值排序从小到大
U = fliplr(U(:, I));                                  % 对应特征矢量排序，fliplr 之后，较大特征值对应的特征矢量在前面
Un = U(:, K+1:M);                                     % 噪声子空间
Gn = Un*Un';
% 空间谱搜索
searchGrids = 0.1;                                    % 搜索间隔
ang = (-90:searchGrids:90)*derad;
Pmusic = zeros(1, length(ang));                       % 空间谱
for i = 1:length(ang)
    a = exp(-1j*pi*idx*sin(ang(i)));
    Pmusic(:, i) = 1/(a'*Gn*a);
end
% 归一化处理，单位化为 db
Pmusic = abs(Pmusic);
Pmusic = 10*log10(Pmusic/max(Pmusic));
% 作图
figure;
plot(ang/derad, Pmusic, '-', 'LineWidth', 1.5);
set(gca, 'XTick',(-90:30:90));
xlabel('\theta(\circ)', 'FontSize', 12, 'FontName', '微软雅黑');
ylabel('空间谱(dB)', 'FontSize', 12, 'FontName', '微软雅黑');
% 找出空间谱Pmusic中的峰值并得到其对应角度
[pks, locs] = findpeaks(Pmusic);                      % 找极大值
[pks, id] = sort(pks);
locs = locs(id(end-K+1:end))-1;
Theta = locs*searchGrids - 90;
Theta = sort(Theta);
disp('估计结果：');
disp(Theta);

5.4 参考内容

1. Shan T J, Wax M, Kailath T. On spatial smoothing for direction-of-arrival estimation of coherent signals[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1985, 33(4): 806-811.
2. Pillai S U, Kwon B H. Forward/backward spatial smoothing techniques for coherent signal identification[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1989, 37(1): 8-15.
3. 【CSDN】基于空间平滑MUSIC算法的相干信号DOA估计(1)
4. 【CSDN】基于空间平滑MUSIC算法的相干信号DOA估计(2)
5. 【CSDN】关于相干信号的阵列处理
6. 【博客园】空间谱专题11：子阵平滑与秩亏缺