SMI波束形成和LSMI波束形成算法原理介绍及MATLAB实现

SMI波束形成

为了最大化输出信噪比，Capon等人提出了在保证对期望信号响应一定的条件下最小化阵列输出功率的方法，可以表示为
$$\{\begin{array}{rl}& \underset{w}{\min }{w}^{H}Rw\\ & s.t.w^H{a}_0=1\end{array}$$
其中，$a_0$表示期望信号的导向矢量。
利用拉格朗日乘数法对上述优化问题进行求解，定义函数
$$J\left(w,\lambda \right)=w^{H}Rw+\lambda \left(w^H{a}_0-1\right)$$
对上式关于$w^{\text{*}}$求导并令导数为零可以得到
$$Rw+\lambda a_0=0$$
即
$$w=-\lambda R^{-1}{a}_0$$
将其带入$w^H{a}_0=1$，可以得到
$$\lambda =-\frac{1}{a_0^H{R}^{-1}{a}_0}$$
因此最优解为
$$w=\frac{R^{-1}{a}_0}{a_0^H{R}^{-1}{a}_0}$$
这就是著名的Capon波束形成算法，在理想的条件下可以最大化输出信干噪比。将加权系数带入优化函数中，可以得到阵列输出的功率为
$$P=w^{H}Rw=\frac{1}{a_0^H{R}^{-1}{a}_0}$$
该功率又称Capon功率。当导向矢量早空间变化时，上述功率变化为整个空间的功率，又称空间谱，常常被用来进行DOA估计。
在实际中，由于$R$无法获得，常用样本协方差矩阵$\hat{R}$代替，此时的波束形成器称为样本协方差矩阵求逆（Sample Matrix Inversion，SMI）波束形成器。其加权系数为
$$w=\frac{{\hat{R}}^{-1}{a}_0}{a_0^H{\hat{R}}^{-1}{a}_0}$$

LSMI波束形成

为了提高SMI的稳健性，LSMI（Loading Sample Matrix Inversion，LSMI）被提了出来。对角加载就是在样本协方差矩阵的对角元素加一个常数，其对应的权系数为
$$w=\frac{{\left(\hat{R}\text{+}\lambda I\right)}^{-1}{a}_0}{a_0^H{\left(\hat{R}\text{+}\lambda I\right)}^{-1}{a}_0}$$

方向图

阵列方向图指的是阵列权矢量与波达方向之间的关系，它反映了该阵列对各个方向信号的接收响应情况，对于$M$阵元的均匀线阵来说，阵列响应可以表示为
$$F\left(\theta \right)=w^{H}a\left(\theta \right)=\sum _{m=0}^{M-1}{w}_m^{\ast }{e}^{-j2\pi d\sin \theta /{\lambda }_0}$$
其中，$w=\left[w_1,w_2,\cdots ,w_{M-1}\right]$为权系数。
对其进行归一化可以得到
$$G\left(\theta \right)=20\lg \frac{\left|F\left(\theta \right)\right|}{\max \left|F\left(\theta \right)\right|}$$
该式刻画了阵列对不同方向信号的幅度增益情况，被称作阵列方向图

仿真参数设置

参数名称	参数值
阵元数	10
期望信号角度	$-5^{\circ }$
干扰信号角度	$-30^{\circ }$、$30^{\circ }$
SNR	10dB
INR	20dB
快拍数	60
在该条件下得到SMI和LSMI的方向图如下
从图中可以看出，LSMI的方向图较SMI有了较大的改善。SMI算法性能不好的原因是快拍数太少，在上述条件不变的情况下，增加快拍数至600，可以得到如下结果

从图中可以看出，两者在期望信号方向形成最大增益接收，在干扰方向上均能形成零陷。此外，还可以看出LSMI有更低的旁瓣，也体现了LSMI较SMI的有更好的效果。
代码如下：

clear;
close all;
clc;
warning off
%% 初始化
M = 10;             %阵元数
fs = 5000;          % 采样频率
f = 1000;           % 信号频率
snap = 600;         % 快拍数
T = 0.5;           %采样时间
t = 1/fs:1/fs:T;
c = 340;
lamda = c/f;              %波长
d = 0.5*lamda;          %阵元间距
theta0 =-5;                %期望信号角度
theta1 =-30;                %干扰角度
theta2 = 30;                %干扰角度
snr=10;                     %信噪比
inr1 =20;                   %干噪比
inr2 = 20;                   %干噪比
snr_noise = 0;              %噪声功率1,为0dBW

%% 导向矢量
a0 = exp(-1j*2*pi*d*sind(theta0)*(0:M-1)'/lamda);
a1 = exp(-1j*2*pi*d*sind(theta1)*(0:M-1)'/lamda);
a2 = exp(-1j*2*pi*d*sind(theta2)*(0:M-1)'/lamda);

%% 信号、干扰和噪声
tar_sig = wgn(1,length(t), snr);
inf1 = wgn(1,length(t),inr1);
inf2 = wgn(1,length(t),inr2);
noise = wgn(M,length(t),snr_noise);

%% 阵列接收信号
rec_sig = a0*tar_sig + a1*inf1 + a2*inf2 + noise;
interference = a1*inf1 + a2*inf2;
sig = a0 * tar_sig;

%% 协方差矩阵
Rx = rec_sig(:,1:snap)*rec_sig(:,1:snap)'/snap;
Rs = sig(:,1:snap)*sig(:,1:snap)'/snap;
Ri = interference(:,1:snap)*interference(:,1:snap)'/snap;
Rn = noise(:,1:snap)*noise(:,1:snap)'/snap;

%% SMI算法

w_smi = inv(Rx)*a0*inv(a0'*inv(Rx)*a0);
theta = -90:0.1:90; % scan angle
p = exp(-1j*2*pi*d*(0:M-1)'*sind(theta)/lamda);
y_smi = w_smi'*p;
yy_smi = 20*log10(abs(y_smi)/max(abs(y_smi)));

%% LSMI算法
LNR = 10;
Loading = LNR * eye(M);
w_lsmi = inv(Rx+Loading)*a0*inv(a0'*inv(Rx+Loading)*a0);
theta = -90:0.1:90; % scan angle
p = exp(-1j*2*pi*d*(0:M-1)'*sind(theta)/lamda);
y_lsmi = w_lsmi'*p;
yy_lsmi = 20*log10(abs(y_lsmi)/max(abs(y_lsmi)));
%% 绘图
figure(1);
plot(theta,yy_smi,'linewidth', 2);
hold on;
plot(theta,yy_lsmi,'linewidth', 2);
xlabel('角度(\circ)');ylabel('归一化增益(dB)');
legend('SMI','LSMI');
grid on;
xlim([-90 90])