基于广义线性组合的Capon波束形成原理介绍及MATLAB实现

原理介绍

在小快拍的情况下，采样协方差矩阵${\hat{\mathbf{R}}}_x$和理想协方差矩阵之间有较大的差距。为了能够提高采样协方差矩阵估计的准确性，研究学者提出了一种基于协方差矩阵的广义线性组合（General-Linear-Combination，GLC）方法。
考虑组合后的协方差矩阵${\tilde{\mathbf{R}}}_x$
$${\tilde{\mathbf{R}}}_x\text{=}\alpha \mathbf{I}+\beta {\hat{\mathbf{R}}}_x$$
其中${\tilde{\mathbf{R}}}_x\ge 0$是半正定的，$\mathbf{I}$为单位阵，$\alpha \ge 0$，$\beta \ge 0$为组合系数，也叫收缩参数，可以通过最小化与理想协方差矩阵的均方误差来确定
$$\begin{array}{rl}& f=E\left\{\Vert {\tilde{\mathbf{R}}}_x-{\mathbf{R}}_x\Vert \right\}\\ & =E\left\{{\Vert \alpha \mathbf{I}-(1-\beta ){\mathbf{R}}_x+\beta ({\hat{\mathbf{R}}}_x-{\mathbf{R}}_x)\Vert }^2\right\}\\ & ={\Vert \alpha \mathbf{I}-(1-\beta ){\mathbf{R}}_x\Vert }^2+{\beta }^{2}E\left\{{\Vert {\hat{\mathbf{R}}}_x-{\mathbf{R}}_x\Vert }^2\right\}\\ & ={\alpha }^{2}M-2\alpha (1-\beta )\text{tr}\left({\mathbf{R}}_x\right)+{\left(1-\beta \right)}^2{\Vert {\mathbf{R}}_x\Vert }^2+{\beta }^{2}E\left\{{\Vert {\hat{\mathbf{R}}}_x-{\mathbf{R}}_x\Vert }^2\right\}\end{array}$$
通过求解上述表达式，可以得到最优的$\alpha$和$\beta$
$$\{\begin{array}{rl}& {\beta }_0=\frac{\gamma }{\rho +\gamma }\\ & {\alpha }_0=\upsilon \left(1-{\beta }_0\right)=\upsilon \frac{\rho }{\gamma +\rho }\end{array}$$
其中，$\rho =E\left\{{\Vert {\hat{\mathbf{R}}}_x-{\mathbf{R}}_x\Vert }^2\right\}$，$\upsilon \text{=tr}\left({\mathbf{R}}_x\right)/M$，$\gamma ={\Vert \upsilon \mathbf{I}-{\mathbf{R}}_x\Vert }^2$。同时，我们可以看出$\beta \in \left[\text{0}\text{1}\right]$以及$\alpha \ge 0$，但是无法求解具体的值，因为上述表达式包含理想协方差矩阵。只有得到$\rho$的值后，才可以去估计$\alpha$和$\beta$，因此首先要对$\rho$进行估计，可以得到以下的表达式
$$\hat{\rho }=\frac{1}{N^2}\sum _{n=1}^N{\Vert \mathbf{x}\left(n\right)\Vert }^4-\frac{1}{N}{\Vert {\hat{\mathbf{R}}}_x\Vert }^2$$
同时，注意到
$$\gamma +\rho ={\Vert \upsilon \mathbf{I}-{\mathbf{R}}_x\Vert }^2+E\left\{{\Vert {\hat{\mathbf{R}}}_x-{\mathbf{R}}_x\Vert }^2\right\}=E\left\{{\Vert {\hat{\mathbf{R}}}_x-\upsilon \mathbf{I}\Vert }^2\right\}$$
因此，可以通过采样协方差矩阵${\hat{\mathbf{R}}}_x$来估计$\hat{\upsilon }$
基于上述过程，可以得到$\alpha$和$\beta$的估计值
$$\{\begin{array}{rl}& {\hat{\beta }}_0=\text{1}-\frac{\hat{\rho }}{{\Vert {\hat{\mathbf{R}}}_x-\hat{\upsilon }\mathbf{I}\Vert }^2}\\ & {\hat{\alpha }}_0=\hat{\upsilon }\left(1-{\beta }_0\right)=\hat{\upsilon }\frac{\hat{\rho }}{{\Vert {\hat{\mathbf{R}}}_x-\hat{\upsilon }\mathbf{I}\Vert }^2}\end{array}$$
但是由于$\alpha \ge 0$，$\beta \ge 0$，而上面的求解方法得到的结果可能是负值，所以可以采用下列表达式进行计算
$$\{\begin{array}{rl}& {\hat{\alpha }}_0=\min \left[\hat{\upsilon }\frac{\hat{\rho }}{{\Vert {\hat{\mathbf{R}}}_x-\hat{\upsilon }\mathbf{I}\Vert }^2},\hat{\upsilon }\right]\\ & {\hat{\beta }}_0=1-\frac{{\hat{\alpha }}_0}{\hat{\upsilon }}\end{array}$$
基于估计的系数，可以得到组合后的协方差矩阵
$${\tilde{\mathbf{R}}}_x\text{=}{\hat{\alpha }}_0\mathbf{I}+{\hat{\beta }}_0{\hat{\mathbf{R}}}_x$$
那么，GLC波束形成器的加权系数可以表示为
$$\mathbf{w}=\frac{{\tilde{\mathbf{R}}}_x^{-1}{\mathbf{a}}_0}{{\mathbf{a}}_0^H{\tilde{\mathbf{R}}}_x^{-1}{\mathbf{a}}_0}$$
如果${\hat{\beta }}_0\ne 0$，那么加权系数还可以表示为
$$\mathbf{w}=\frac{{\left[\frac{{\hat{\alpha }}_0}{{\hat{\beta }}_0}\mathbf{I}+{\hat{\mathbf{R}}}_x\right]}^{-1}{\mathbf{a}}_0}{{\mathbf{a}}_0^H{\left[\frac{{\hat{\alpha }}_0}{{\hat{\beta }}_0}\mathbf{I}+{\hat{\mathbf{R}}}_x\right]}^{-1}{\mathbf{a}}_0}$$
可以看出，其等效为对角加载，而且可以自动计算对角加载的大小，避免了如何确定对角加载量大小的问题，所以该方法又称作为自动对角加载技术。

仿真参数设置

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

%% GLC算法
v=trace(Rx)/M;
p_sum=norm(rec_sig,'fro')^4;

rho=p_sum/(snap^2)-norm(Rx,'fro')^2/snap;
alpha=min(v*rho/(norm(Rx-v*eye(size(Rx)))^2),v);
beta=1-alpha/v;
R_glc=alpha*eye(M)+beta*Rx;
w_glc =inv(R_glc)*a0/(a0'*inv(R_glc)*a0);
theta = -90:0.1:90; % scan angle
p = exp(-1j*2*pi*d*(0:M-1)'*sind(theta)/lamda);
y = w_glc'*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] L. Du, J. Li and P. Stoica, “Fully Automatic Computation of Diagonal Loading Levels for Robust Adaptive Beamforming,” in IEEE Transactions on Aerospace and Electronic Systems, vol. 46, no. 1, pp. 449-458, Jan. 2010, doi: 10.1109/TAES.2010.5417174.
[2]Ledoit, Olivier & Wolf, Michael, 2004. “A well-conditioned estimator for large-dimensional covariance matrices,” Journal of Multivariate Analysis, Elsevier, vol. 88(2), pages 365-411, February.
[3]P. Stoica, J. Li, X. Zhu and J. R. Guerci, “On Using a priori Knowledge in Space-Time Adaptive Processing,” in IEEE Transactions on Signal Processing, vol. 56, no. 6, pp. 2598-2602, June 2008, doi: 10.1109/TSP.2007.914347.