最差性能最优波束形成算法原理介绍及MATLAB实现

原理介绍

标准的Capon波束形成是最小化阵列输出功率，并且对期望信号无失真接收来达到最大化输出信干噪比的目的。但是当真实导向矢量未知且存在一个不确定集中时，为了对真实导向矢量设置无失真响应约束，最差情况性能最优算法通过对不确定集合中所有可能的导向矢量进行无失真响应约束来确保对期望信号的无失真接收，可以表示为
$$\{\begin{array}{rl}& \underset{w}{\min }{w}^H\hat{R}w\\ & s.t.\left|w^{H}a\right|\ge 1,a\in S_a\end{array}$$
其中，$S_a$为球形不确定集。对于每一个$a\in S_a$来说，$\left|w^{H}a\right|\ge 1$表示一个非凸的二次约束函数，也就是说，上述问题相当于含有无数个非凸的二次约束问题，是一个NP-hard问题，然而可以通过转化，将其转化为一个凸问题，进而利用内点法进行求解。
该问题等价于
$$\underset{a\in S_a}{\min }\left|w^{H}a\right|\ge 1$$
令 $a=\bar{a}+e$，那么${\Vert e\Vert }_2\le \epsilon$。根据柯西许瓦兹不等式可以得到
$$\left|w^{H}a\right|=\left|w^H\bar{a}+w^{H}e\right|\ge \left|w^H\bar{a}\right|-\left|w^{H}e\right|\ge \left|w^H\bar{a}\right|-\epsilon {\Vert w\Vert }_2$$
因此上述问题可以转化为
$$\left|w^H\bar{a}\right|-\epsilon {\Vert w\Vert }_2\ge 1$$
该约束函数由于左边绝对值运算仍然是一个非凸问题。但是最初的优化问题的目标函数对$w$的相位不敏感，那么可以调节$w$的相位使得$w^H\bar{a}$是一个实数，即
$$\{\begin{array}{rl}& \mathrm{Re}\left\{w^H\bar{a}\right\}\ge 0\\ & \mathrm{Im}\left\{w^H\bar{a}\right\}=0\end{array}$$
联合上述问题，原来的非凸问题可以转化为
$$\{\begin{array}{rl}& \underset{w}{\min }{w}^H\hat{R}w\\ & s.t.\left|w^H\bar{a}\right|-\epsilon {\Vert w\Vert }_2\ge 1\mathrm{Im}\left\{w^H\bar{a}\right\}=0\end{array}$$
该问题是一个很难获取闭式解的凸优化问题，可进一步转化为二阶锥规划（Second-Order Cone Programming, SOCP）问题，通过内点法求解。

仿真参数设置

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

仿真结果如下

从结果来看，WCP能够早期望信号方向形成最大增益，在干扰方向上形成零陷，从而达到了增强期望信号抑制干扰的效果。
代码如下：

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 = 1;              %噪声功率

%% 导向矢量
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; 

%% WCP算法
epsilon = 3; 
Cov = Rx; % 采样协方差矩阵
cvx_begin quiet
    variable w(M) complex;
    minimize norm(sqrt(Cov)*w)
subject to
    epsilon * norm(w,2) <= real(w'*a0) - 1;
    imag(w'*a0) == 0;
cvx_end
optimal_weight = w;
theta = -90:0.1:90; % scan angle
p = exp(-1j*2*pi*d*(0:M-1)'*sind(theta)/lamda);
y = optimal_weight'*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])

参考文献：
[1]S. A. Vorobyov, A. B. Gershman and Zhi-Quan Luo, “Robust adaptive beamforming using worst-case performance optimization: a solution to the signal mismatch problem,” in IEEE Transactions on Signal Processing, vol. 51, no. 2, pp. 313-324, Feb. 2003, doi: 10.1109/TSP.2002.806865.