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

原理介绍

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

仿真参数设置

参数名称	参数值
阵元数	10
期望信号角度	−5∘-5^{\circ}−5∘
干扰信号角度	−30∘-30^{\circ}−30∘、30∘30^{\circ}30∘
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.