基于不确定集的波束形成算法原理介绍及MATLAB实现

原理介绍

该算法是通过将期望信号的导向矢量约束在一个不确定集上面，通过优化找到一个比名义导向矢量更加准确的导向矢量来提高波束形成器的稳健性。其中最具有代表性的是稳健Capon波束形成算法，其主要思想是在一定不确定区域内寻找使得阵列输出功率最大的那个导向矢量作为最佳估计值，以球形不确定集为例
S a = { a ∣ ∥ a ˉ 0 − a ∥ 2 < ε } {{S}_{a}}=\left\{ a|{{\left\| {{{\bar{a}}}_{0}}-a \right\|}_{2}}<\varepsilon \right\} Sa​={a∣∥aˉ0​−a∥2​<ε}
其中，$\epsilon$为大于0的常数，不确定集$S_a$描述的是以${\bar{a}}_0$为球心半径为$\epsilon$的球形区域，那么稳健的Capon波束形成可以描述为
$$\{\begin{array}{rl}& \underset{w}{\min }{w}^{H}Rw\\ & s.t.{\Vert {\bar{a}}_0-a\Vert }_2<\epsilon \end{array}$$
为了排除平凡解$a=0$， 需要$\epsilon <\sqrt{M}$
通过拉格朗日乘数法对其进行求解可以得到
$$\hat{a}={\left(\frac{R^{-1}}{\lambda }+I\right)}^{-1}{\bar{a}}_0={\bar{a}}_0-{\left(I+\lambda R\right)}^{-1}{\bar{a}}_0$$
带入约束函数取等可以得到
$$\Vert {\left(I+\lambda R\right)}^{\text{-}1}{\bar{a}}_0\Vert \text{=}\epsilon$$
最终求得的稳健Capon波束形成的权系数为
$$w=\frac{R^{-1}{a}_0}{a_0^H{R}^{-1}{a}_0}=\frac{{\left(R+\frac{1}{\lambda }I\right)}^{-1}{\bar{a}}_0}{{\bar{a}}_0^H{\left(R+\frac{1}{\lambda }I\right)}^{-1}R{\left(R+\frac{1}{\lambda }I\right)}^{-1}{\bar{a}}_0}$$
从形式上看，稳健Capon波束形成（Robust Capon Beamforming, RCB）算法属于对角加载算法范畴，只不过是将求解加载系数的问题转化为了确定不确定集系数问题，将通过不确定集约束间接估计出加载系数。

仿真参数设置

参数名称	参数值
阵元数	10
期望信号角度	$-5^{\circ }$
干扰信号角度	$-30^{\circ }$、$30^{\circ }$
SNR	10dB
INR	20
快拍数	60

基于上述仿真参数，可以得到基于不确定集的波束方向图为

从图中可以看出，形成的方向图在期望信号方向形成了最大增益，在干扰方向上形成了零陷，同时可以通过增加误差来进一步显示该方法的优越性，这里就不再叙述。
代码如下：

clear;
close all;
clc;
warning off

%% 初始化
M = 10;             %阵元数
fs = 5000;          % 采样频率
f = 1000;           % 信号频率
snap = 60;         % 快拍数
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; 

%% RCB算法
epsilon = 3; 
R_inv = inv(Rx); 
cvx_begin
    variable a(M) complex;
    minimize norm(sqrtm(R_inv)*a)
subject to
   norm(a0-a) <= epsilon;
cvx_end
w_rcb = R_inv*a/(a'*R_inv*a);
theta = -90:0.1:90; % scan angle
p = exp(-1j*2*pi*d*(0:M-1)'*sind(theta)/lamda);
y = w_rcb'*p;
yy = 20*log10(abs(y)/max(abs(y)));
%% 绘图
figure(1);
plot(theta,yy,'linewidth', 2);
xlabel('角度(\circ)');ylabel('归一化增益(dB)')
grid on;
xlim([-90 90])