基于线性积分的波束形成方法原理介绍及MATLAB实现

原理介绍

基于线性积分重构干扰加噪声协方差矩阵是首次利用重构的协方差矩阵进行波束形成的，其性能较之前的波束形成方法有了很大的改进。
线性积分重构的原理就是在不包含期望信号的角度内进行积分重构干扰加噪声协方差矩阵，利用重构后的协方差矩阵进行加权系数的计算。

协方差矩阵重构

Capon功率谱为
$$P\left(\theta \right)=\frac{1}{{\bar{\mathbf{a}}}^H\left(\theta \right){\hat{\mathbf{R}}}_x^{-1}\bar{\mathbf{a}}\left(\theta \right)}$$
其中，$\bar{\mathbf{a}}\left(\theta \right)$表示角度$\theta$的名义导向矢量
在除去信号区域内利用Capon功率谱积分重构干扰加噪声协方差矩阵
$$\begin{array}{rl}& {\hat{\mathbf{R}}}_{i+n}^{-1}={\int }_{\bar{\Theta }}^{}P\left(\theta \right)\bar{\mathbf{a}}\left(\theta \right){\bar{\mathbf{a}}}^H\left(\theta \right)d\theta \\ & ={\int }_{\bar{\Theta }}^{}\frac{1}{{\bar{\mathbf{a}}}^H\left(\theta \right){\hat{\mathbf{R}}}_x^{-1}\bar{\mathbf{a}}\left(\theta \right)}\bar{\mathbf{a}}\left(\theta \right){\bar{\mathbf{a}}}^H\left(\theta \right)d\theta \end{array}$$
其中，$\bar{\Theta }$表示除去期望信号的区域，下图为积分示意图

通过在去除期望信号的区域进行积分获得干扰加噪声协方差矩阵。

期望信号导向矢量估计

最初的期望信号导向矢量估计，都是将其约束在一个不确定集内，最大化阵列输出功率，也就可以表示为
$$\begin{array}{rl}& \underset{\mathbf{e}}{\min }{\left({\bar{\mathbf{a}}}_0+\mathbf{e}\right)}^H{\hat{\mathbf{R}}}_{i+n}^{-1}\left({\bar{\mathbf{a}}}_0+\mathbf{e}\right)\\ & s.t.{\left({\bar{\mathbf{a}}}_0+\mathbf{e}\right)}^H{\hat{\mathbf{R}}}_{i+n}\left({\bar{\mathbf{a}}}_0+\mathbf{e}\right)\le {\bar{\mathbf{a}}}_0^H{\hat{\mathbf{R}}}_{i+n}{\bar{\mathbf{a}}}_0\end{array}$$
其中，$\mathbf{e}$是导向矢量的误差不等式的约束时为了防止期望信号导向矢量收敛到非期望信号区域内。
由于真实的导向矢量${\mathbf{a}}_0$可以表示为
$${\mathbf{a}}_0={\bar{\mathbf{a}}}_0+\mathbf{e}={\bar{\mathbf{a}}}_0+{\mathbf{e}}_{\perp }+{\mathbf{e}}_{\parallel }$$
而${\mathbf{e}}_{\parallel }$是不影响波束形成器的性能的，因此，上述优化问题可以表示为
$$\begin{array}{rl}& \underset{{\mathbf{e}}_{\perp }}{\min }{\left({\bar{\mathbf{a}}}_0+{\mathbf{e}}_{\perp }\right)}^H{\hat{\mathbf{R}}}_{i+n}^{-1}\left({\bar{\mathbf{a}}}_0+{\mathbf{e}}_{\perp }\right)\\ & s.t.{{\bar{\mathbf{a}}}_0}^H{\mathbf{e}}_{\perp }=0\\ & {\left({\bar{\mathbf{a}}}_0+{\mathbf{e}}_{\perp }\right)}^H{\hat{\mathbf{R}}}_{i+n}\left({\bar{\mathbf{a}}}_0+{\mathbf{e}}_{\perp }\right)\le {{\bar{\mathbf{a}}}_0}^H{\hat{\mathbf{R}}}_{i+n}{\bar{\mathbf{a}}}_0\end{array}$$
最终的期望信号导向矢量估计可以表示为
$${\hat{\mathbf{a}}}_0={\bar{\mathbf{a}}}_0+{\mathbf{e}}_{\perp }$$
波束形成器的加权系数可以表示为
$$\mathbf{w}=\frac{{\hat{\mathbf{R}}}_{i+n}^{-1}{\hat{\mathbf{a}}}_0}{{\hat{\mathbf{a}}}_0^H{\hat{\mathbf{R}}}_{i+n}^{-1}{\hat{\mathbf{a}}}_0}$$
后续有研究线性积分推广到环积分，也就是将其推广到高维空间，能够获得更好的效果，但是运算量较大，感兴趣的可以参考Robust Adaptive Beamforming With a Novel Interference-Plus-Noise Covariance Matrix Reconstruction Method这篇文章，因为原理较为类似，这里就不再介绍。

仿真参数设置

参数名称	参数值
阵元数	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;

%% 积分重构
theta_range1 =-90:0.1:-10;
theta_range2 = 0:0.1:90;
C1 = zeros(M,M);
for iter1 = 1:length(theta_range1)
    a1 = exp(-1j*2*pi*.5*sind(theta_range1(iter1))*(0:M-1)'); 
    C1 = C1 + inv(a1'*inv(Rx)*a1)*a1*a1';
end
C2 = zeros(M,M);
for iter2 = 1:length(theta_range2)
    a2 = exp(-1j*2*pi*.5*sind(theta_range2(iter2))*(0:M-1)'); 
    C2 = C2 + inv(a2'*inv(Rx)*a2)*a2*a2';
end
C_recon = C1 + C2; %重构协方差矩阵
R_inv = inv(C_recon);

cvx_begin
    variable e(M) complex
    minimize norm(sqrtm(R_inv)*(a0+e))
subject to
    a0' * e == 0;
    norm(sqrtm(C_recon)*(a0+e)) <= norm(sqrtm(C_recon)*a0);
cvx_end
a0_recon = a0 + e;
w = R_inv*a0_recon*inv(a0_recon'*R_inv*a0_recon);
theta = -90:0.1:90; 
p = exp(-1j*2*pi*d*(0:M-1)'*sind(theta)/lamda);
y = w'*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]Y. Gu and A. Leshem, “Robust Adaptive Beamforming Based on Interference Covariance Matrix Reconstruction and Steering Vector Estimation,” in IEEE Transactions on Signal Processing, vol. 60, no. 7, pp. 3881-3885, July 2012, doi: 10.1109/TSP.2012.2194289.