毫米波雷达系列 | 基于前后向空间平滑的MUSIC算法详解

忙了一阵子的中期和专利，基本上告一段落，简单的写一个比较常见的解相干MUSIC角度估计算法。

算法是基于MIMO体制雷达的基础上实现的，MIMO体制形成虚拟接收阵元可以来改善毫米波雷达的成像效果，但是角度分辨率从仍受限于阵列孔径大小，为了提高MIMO毫米波雷达的成像分辨率只能以增加系如成本为代价继续增加天线阵列。

DOA阵列模型

假设信号均为窄带信号，前方有k个目标则接受阵列信号可以表示为：
$$x(t)=A(\theta )S(t)+n(t)$$
其中，$S(t)=[s1(t),s2(t),...sk(t){]}^T$是入射信号，$A(\theta )=[\alpha 1,\alpha 2,...\alpha k]$为接收阵列的导向矢量，$n(t)=[n1(t),n2(t),...nN(t){]}^T$是阵列接收到的噪声矢量，噪声功率为σ² ,展开可得

MUSIC算法

MUSIC算法的核心思想就是将阵列接收数据的协方差矩阵进行特征值分解得到特征值和特征向量，将特征向量划分为信号子空间和噪声子空间，这两个空间是正交且不相关的，利用该性质可以实现空间信号的超分辨。

接受阵列数据的协方差矩阵可以表示为：
$$R=E[x{x}^H]=AE[S{S}^H]+{\sigma }^{2}I=ARsA++{\sigma }^{2}I$$
其中，$Rs=E[S{S}^H]$表示信号的自相关矩阵，ARsA为信号相关部分，σ²I为噪声相关部分
对R进行特征值分解：
$$R=U\Sigma U^H$$
其中U是特征向量组成的矩阵，Σ是所有特征值所构成的对角矩阵，Σ可以表示为：

Σ中的对角线上的各个特征值满足以下条件：

这里假设有两个对角矩阵：

前一个为K个较大的特征值所构成的对角矩阵, 为剩下的(N-K)个小特征值构成的对角矩阵。这样可以将特征向量矩阵分为两个部分：第一部分为K个大特征值对应的信号子空间 ，第二部分为(N-K)个小特征值对应的噪声子空间 ，可以重写为：

在信号源相互独立的条件下，可以知道导向向量与信号噪声子空间正交，且和也是正交的，所以可以得到：

在实际工程应用当中，我们得到的接收数据都是在有限时间内得到的所以一般用协方差矩阵的最大似然估计来代替数据协方差矩阵：

所以谱空间函数可以表示为：

空间平滑算法

空间平滑技术的基本思想是：如果存在远场窄带且相干的电磁信号照射线性天线阵列，那么这些阵列对信号进行协方差处理后的矢量有可能是满秩的。基于此可以将该阵列分割成很多子阵，计算时将所有的协方差矢量进行和相加，这样可能得到一个矩阵，它的秩的大的数量值相等于相干信号源的个数，对该矩阵矢量实施子空间的算法可算出每个信源的波达方向角。一般应用中，空间平滑算法只适合于阵元等间隔的线性阵(ULA)。

假设存在一个均匀阵列，阵元个数为N，阵元间距为d，空间中有p(p<N)个窄带信号,首先使用前向空间平滑把阵列划分为M个相互重叠的前向子阵，子阵的阵元数相等均为k，可得知k=N-M+1。

则第m个前向子阵的输出为：
$$x_m^f(t)=[x_m(t),x_{m+1}(t),...,x_{m+k-1}(t){]}^T=A_k{D}^{(m-1)}s(t)+n_m(t)$$
其中D为：

作协方差运算可以得到：
R m f = E { x m f ( t ) x m f ( t ) H } R_m^f=E\lbrace x_m^f(t)x_m^f(t)^H\rbrace Rmf​=E{xmf​(t)xmf​(t)H}
定义前向空间平滑矩阵为：
$$R^f=1/M\sum _{m-1}^M{R}_m^f$$
后向与前向的滑动方向相反，可得：
$$R^b=1/M\sum _{m-1}^M{R}_m^b$$
$R^b$是$R^f$的共轭倒序矩阵，两者之间的关系满足共轭倒叙不变性，可以表示为：
$$R^f=J(R^f{)}^{\ast }J$$
矩阵$J$为副对角线元素均为1，因此前后向空间平滑协方差矩阵的定义为：
$$R^{fb}=1/2(R^f+R^b)$$
可以利用该协方差矩阵去求解特征值继续MUSIC算法来完成角度估计。

整体流程

仿真代码

clear all
close all
format long
snr = 0;
N = 256;
source_number=2;%信元数
sensor_number=8;%阵元数
source_doa =[20 25];

%% 前后向解相干
w=[pi/4 pi/4].';%信号频率
l=((2*pi*3e8)/w(1)+(2*pi*3e8)/w(2))/2;%信号波长  
d=0.5*l;%阵元间距
m=6;%每个子阵阵元数
p=3;%相互交错的子阵数
snr=10;%信噪比
A=[exp(-j*(0:sensor_number-1)*d*2*pi*sin(source_doa(1)*pi/180)/l);exp(-j*(0:sensor_number-1)*d*2*pi*sin(source_doa(2)*pi/180)/l)].';%阵列流型

s=sqrt(10.^(snr/10))*exp(j*w*[0:N-1]);%仿真信号
x = A*s;
x=awgn(x,snr);

x1=x([1:6],:);R1=x1*x1'/1024;
x2=x([2:7],:);R2=x2*x2'/1024;
x3=x([3:8],:);R3=x3*x3'/1024;

Rf=(R1+R2+R3)/3;

[U,S,V]=svd(Rf);
Un=U(:,source_number+1:m);
Gn=Un*Un';

searching_doa=-90:0.1:90;%线阵的搜索范围为-90~90度
 for i=1:length(searching_doa)
   a_theta=exp(-j*(0:m-1)'*2*pi*d*sin(pi*searching_doa(i)/180)/l);
   Pmusic(i)=1./abs((a_theta)'*Gn*a_theta);
 end
plot(searching_doa,10*log(Pmusic));
%axis([-90 90 -90 90]);
hold on;
grid on;
title('前后向空间平滑MUSIC')
xlabel('入射角/度');
ylabel('谱峰/dB');
grid on;