波束形成基于体积分的协方差矩阵重建稳健自适应波束形成算法

量大龙神魔

已于 2023-09-11 21:20:35 修改

阅读量823

点赞数 6

文章标签：算法矩阵 matlab

于 2023-04-11 20:30:41 首次发布

本文链接：https://blog.csdn.net/qq_45471796/article/details/129760016

版权

前言

文中可能会出现错误，如果诸位大佬发现，请不吝指出谢谢。

基于协方差矩阵重建的稳健自适应波束形成算法，其能够比其他算法更为有效的减少甚至抑制协方差矩阵中的期望信号分量，从而在较高信(干)噪比下，相比于传统的信号协方差阵修正的波束形成算法(e.g. 对角加载、子空间投影等)，有效地防止了信号的自消现象的发生。早期的协方差矩阵重构的波束形成算法主要通过积分尽可能地还原干扰加噪声协方差阵在Capon功率谱中的特征，但是该方法很难在算法计算量与鲁棒性上均达到较好的水平，近年来协方差矩阵算法通过估计噪声功率与信号导向矢量对干扰加噪声协方差矩阵进行重构，这在大大降低了算法计算量的同时又较好地还原干扰加噪声协方差阵在Capon功率谱中的特征。

这篇文章是基于Y. Gu等人在2012年发表的文章Robust Adaptive Beamforming Based on Interference Covariance Matrix Reconstruction and Steering Vector Estimation中的工作，由L. Huang等人在2015年发表。文章的主要工作是针对线积分协方差矩阵重构法对导向矢量失配情况下鲁棒性不足的问题，将基于角度域上线积分的干扰加噪声协方差矩阵重建法拓展到基于Capon谱上高维体积分的协方差矩阵。

算法流程

干扰加噪声协方差矩阵重建

一般情况下，真实导向矢量与名义导向矢量一般满足球形不确定集的关系，即：

$\left \| \bar a(\theta) -a \right \|_{2} \leq \varepsilon,s.t.\left \| \bar a(\theta) \right \|_2 = \sqrt{M}$

其中 $\bar a(\theta)$ 为名义导向矢量， $a$ 为真实的导向矢量， $\varepsilon$ 为不确定系数。因此作者的思路主要是通过对每个干扰角度域上的名义导向矢量的球形不确定集中的所有导向矢量进行积分，即

${\hat R}_{int}^{'} = \int_{S_a(\theta \in \Theta_{int})}^{}\frac{aa^H}{a^H {\hat R}^{-1}a}dv$

如下图所示：

其中集合 $S_a(\theta \in \Theta_{int})$ 定义如下：

$S_a(\theta \in \Theta_{int}) = \{ a| \left \| \bar a (\theta) - a \right \|_2 \leq \varepsilon,\theta \in \Theta_{int} \}$

然而该集合中的很多导向矢量都是冗余的，一些导向矢量之间的关系仅为倍数。因此需要对积分操作进行进一步改进。已知导向矢量在协方差矩阵中的体现与导向矢量的幅度无关，即共线的导向矢量构造的协方差矩阵元素 $\frac{aa^H}{a^H {\hat R}^{-1}a}$ 是完全相同的，因此需要将该体积分转化为面积分减少冗余导向矢量对协方差矩阵的作用，即：

${\hat R}_{int}^{} = \frac{1}{2}\int_{\partial S_a(\theta \in \Theta_{int})}^{}\frac{aa^H}{a^H {\hat R}^{-1}a}d\sigma$

式中 $\partial S_a(\theta \in \Theta_{int})=\{ a| \left \| \bar a (\theta) - a \right \|_2 = \varepsilon,\theta \in \Theta_{int} \}$ 表示图中高维环的表面，事实上对该环表面进行面积分仍然有一倍的导向矢量冗余，即该高维环的下平面中的导向矢量总有一个对应的导向矢量位于上平面。因此该面积分需要乘以 $\frac{1}{2}$ 以防止重复计算。

然而在理论上对该环表面进行面积分难以达到，因此会对该面积分进行以下近似：

${\hat R}_{int}^{} \approx \frac{1}{2} \int_{\Theta_{int}}^{}\int_{\partial S_a(\theta)}^{}\frac{aa^H}{a^H {\hat R}^{-1}a}d\sigma d\theta=\frac{1}{2} \int_{\Theta_{int}}^{}r(\theta)d\theta$

其中 $r(\theta) = \int_{\partial S_a(\theta)}^{}\frac{aa^H}{a^H {\hat R}^{-1}a}d\sigma$ 即为协方差阵的矩阵元素。

实际计算过程中，连续积分难以实现，因此作者对该二重积分及协方差阵的矩阵元素进一步离散近似：

$r(\theta_i) \approx \sum_{l=1}^{L} \frac{a_{il}a_{il}^H}{a_{il}^H {\hat R}^{-1}a_{il}}$

${\hat R}_{int}^{} \approx \frac{1}{2} \sum_{i=1}^{I}r(\theta_i) \approx \frac{1}{2} \sum_{i=1}^{I}\sum_{l=1}^{L} \frac{a_{il}a_{il}^H}{a_{il}^H {\hat R}^{-1}a_{il}},\theta_i \in \Theta_{int}$

其中：

$a_{il} = \bar a(\theta_i)+ \frac{\varepsilon}{\sqrt{M}}\left [ e^{j\phi_0^l} , e^{j\phi_1^l} ,..., e^{j\phi_{M-1}^l} \right ], \phi_m^l \in~ [0,2 \pi), m=~0,1,...,M-~1$

在实际仿真过程中为方便计算，作者在每个干扰角度域上均采样 $2^M$ 个点，M为阵元数。

噪声功率估计

作者此处采用的是将原协方差矩阵的最小的特征值作为噪声功率 $\sigma_n^2$ 。

导向矢量估计

作者此处采用构造一个QCQP问题对导向矢量进行求解。该QCQP问题可以表示为以下形式：

$\min_{e_\perp } (\bar a_{0}+e_\perp)^H \hat R^{-1} (\bar a_{0}+e_\perp) ,s.t. \bar a_{0}^He_\perp =~0,(\bar a_{0}+~e_\perp)^H \hat R^{-1}_{i+n} (\bar a_{0}+~e_\perp) \leq~ \bar a_{0}^HR^{-1}_{i+n} \bar a_{0}$

通过CVX求解该优化问题很容易得到期望信号的估计导向矢量。

则估计得到的干扰加噪声协方差矩阵为：

$\hat R_{i+n} = \hat R_{int}+ \hat {\sigma}_n^2 I$

阵列权向量为：

$w = \frac{\hat R^{-1}_{i+n} \hat a_0}{\hat a_0^H \hat R^{-1}_{i+n} \hat a_0}$

仿真

部分仿真条件如下所示：

阵元数	10
信噪比	20
扰噪比	20
阵元间距	半波长
信源数	1
干扰数	2
信源到达角	5°
干扰到达角	-30°，45°
信源角度估计误差	[0°,3°]
不确定系数 $\varepsilon$	$\sqrt{0.1}$

仿真代码如下所示(运行该matlab代码需要安装cvx工具包)：

clc;
clear;
jk = 0;
loop = 50;
SINRO = zeros(1,length(-20:5:40));
snr = 20;
deg_dev_predf = 3;%%预定义的角度最大偏移量
ang_mismatch1 = randi(deg_dev_predf*2)-deg_dev_predf;
ang_mismatch2 = randi(deg_dev_predf*2)-deg_dev_predf;
ang_mismatch3 = randi(deg_dev_predf*2)-deg_dev_predf; %%三个不同方向的DOA误差
%% 初始化及设定参数
array_num = 10;%%阵元数
snapshot_num = 30;%%快拍数
source_aoa = [5,-30,45];%%信源到达角
tgt_ang = 1;
c = 340;%%波速
f = 1000;%%频率
lambda = c/f;%%波长
d = 0.5*lambda;
inr1 = 20;
inr2 = 20;
source_num = length(source_aoa);%%信源数
sig_nr = [snr,inr1,inr2];%%信号功率dBw
deg_deviate = 3;%%误差角度偏离
source_dev = source_aoa+[ang_mismatch1,ang_mismatch2,ang_mismatch3];
i_dev = sort(source_dev(1,2:length(source_dev)),'ascend');
sample_p = 40;
reco_size = 8;
i_region = [];
for ik = 1:length(i_dev)
    i_region = [i_region linspace(i_dev(ik)-reco_size,i_dev(ik)+reco_size,floor(sample_p/length(i_dev)))];%%确定干扰域
end
epsl = sqrt(0.1);
%% 转向矢量
a_bar = exp(-1i*(0:array_num-1)'*pi*sind((source_dev(tgt_ang))));%%实际带有误差的阵列响应矢量
A = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(source_aoa));%%阵列响应矩阵
tar_sig = wgn(1,snapshot_num,sig_nr(tgt_ang));
inf1 = wgn(1,snapshot_num,sig_nr(2));
inf2 = wgn(1,snapshot_num,sig_nr(2));
n = (randn(array_num,snapshot_num)+randn(array_num,snapshot_num)*1i)/sqrt(2);
rec_sig = A(:,1)*tar_sig+A(:,2)*inf1+A(:,3)*inf2+n;
R = rec_sig*rec_sig'/snapshot_num;%%接收阵的协方差矩阵
[U,Gamma,V] = svd(R);
Rs = A(:,1)*tar_sig*(A(:,1)*tar_sig)'/snapshot_num;
Ri = (A(:,2)*inf1+A(:,3)*inf2)*(A(:,2)*inf1+A(:,3)*inf2)'/snapshot_num;
Rn = n*n'/snapshot_num;
%% 算法
pow_n = min(diag(Gamma));%%噪声功率
e_base = exp(1i*pi)*ones(array_num,1);
R_recons = zeros(array_num,array_num);
sample_m = 2^array_num;%%单个弧上的采样点数
for mk = 1:length(i_region)
    a_bar_i = exp(-1i*(0:array_num-1)'*pi*sind((i_region(mk))));
    for kk = 1:sample_m
    arr_out = numbertoarray(kk-1,array_num);
    arr_out = (arr_out-0.5)*2;
    a_il = a_bar_i+epsl*(arr_out.*e_base)/sqrt(array_num);
    R_recons = R_recons+a_il*a_il'/(a_il'*inv(R)*a_il);
    end
end
R_recons = R_recons+pow_n*eye(array_num);
R_recons = R_recons/2;
R_inv_recons = inv(R_recons);
cvx_begin quiet
    variable a_e(array_num) complex
    minimize norm(sqrtm(R_inv_recons)*(a_bar+a_e));
subject to
    a_bar'*a_e == 0;
    norm(sqrtm(R_recons)*(a_bar+a_e)) <= norm(sqrtm(R_recons)*a_bar);
cvx_end
a_est = a_bar+a_e;
w0 = inv(R_recons)*a_est/(a_est'*inv(R)*a_est);
beam_plot(w0);

function arr_out = numbertoarray(num_in,arr_size)
    arr_out = zeros(arr_size,1);
    char_out = dec2bin(num_in);
    for zk = 1:length(char_out)
        arr_out(arr_size+1-zk,1) = str2double(char_out(length(char_out)+1-zk));
    end
end
function [] = beam_plot(input_weight_vector)
array_num = length(input_weight_vector);
theta = -90:0.01:90;
p = exp(-1j*2*pi*(0:array_num-1)'*sind(theta)/2);
y = input_weight_vector'*p;
outputval = 20*log10(abs(y)/max(abs(y)));
figure()
plot(theta,outputval);

仿真运行结果如下图所示：

其中红线为OPT SINR，横轴为输入SNR。

~~个人仿真出的该算法的结果并不是很理想。~~

感谢评论区指正，改正后仿真结果与原文章成功对应。

部分参考文献

Y. Gu and A. Leshem. Robust Adaptive Beamforming Based on Interference Covariance Matrix Reconstruction and Steering Vector Estimation. IEEE Transactions on Signal Processing. vol. 60, no. 7, pp. 3881-3885, July 2012.

L. Huang, J. Zhang, X. Xu and Z. Ye. Robust Adaptive Beamforming With a Novel Interference-Plus-Noise Covariance Matrix Reconstruction Method. IEEE Transactions on Signal Processing. vol. 63, no. 7, pp. 1643-1650, April1, 2015.

量大龙神魔

关注

6
点赞
踩
12

收藏

觉得还不错? 一键收藏
15
评论
波束形成基于体积分的协方差矩阵重建稳健自适应波束形成算法

文献 Robust Adaptive Beamforming With a Novel Interference-Plus-Noise Covariance Matrix Reconstruction Method-TSP-2015
复制链接

扫一扫