波束形成基于协方差矩阵重建与对角加载的稳健自适应波束形成算法

量大龙神魔

已于 2023-06-01 17:11:04 修改

阅读量1.3k

点赞数 1

文章标签：算法矩阵信号处理信息与通信

于 2023-04-17 20:00:24 首次发布

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

版权

该文提出了一种新的稳健自适应波束形成算法，通过子空间投影法修正导向矢量，利用信号加噪声子空间估计信源、干扰和噪声功率，并通过在干扰源协方差矩阵上添加对角加载项来增强干扰抑制能力，尤其是在高信噪比情况下。该方法能有效抑制旁瓣，但可能造成波束旁瓣增益的不稳定性。

摘要由CSDN通过智能技术生成

前言

算法主要分为算法主要分为导向矢量估计、信源、干扰及噪声功率估计、协方差矩阵重构三部分。在导向矢量估计部分，作者采用子空间投影法对估计得到的角度构造的名义导向矢量进行进行进一步修正。在功率估计部分，作者提出了一种新颖的估计方法，通过构造信号加噪声子空间对各个信号、干扰及噪声的功率进行估计。在协方差矩阵重建部分，作者根据扰噪比，对每个信号及干扰的功率矩阵上增加对角加载项，对该稳健自适应波束形成算法的干扰抑制能力进行进一步调整。

导向矢量估计

作者采用的是在先验角度旁边的一定范围内寻找Capon功率谱最大处，其对应的角度的名义导向矢量即为估计得到的导向矢量。

信源、干扰及噪声功率估计

一般情况下，注意到

$\textbf{U}_s\textbf{U}_s^H+\textbf{U}_n\textbf{U}_n^H=\textbf{I}$

阵列协方差矩阵可被分解为：

$\begin{aligned}\hat{\textbf{R}}_x=\sum_{m=1}^M\lambda_m\textbf{u}_m\textbf{u}_m^H=\textbf{U}\Lambda\textbf{U}^H=\textbf{U}_s\Lambda_s\textbf{U}_s^H+\textbf{U}_n\Lambda_n\textbf{U}_n^H=\textbf{U}_s\Lambda_s\textbf{U}_s^H+\hat{\sigma}_n^2\left(\textbf{I}-\textbf{U}_s\textbf{U}_s^H\right)\end{aligned}$

在空域中随机选取L个导向矢量，则有：

$\mathbf{a}_l^H\hat{\textbf{R}}_x\mathbf{a}_l=\mathbf{a}_l^H[\textbf{U}_s\Lambda_s\textbf{U}_s^H+\hat{\sigma}_{nl}^2\left(\textbf{I}-\textbf{U}_s\textbf{U}_s^H\right)]\mathbf{a}_l =\mathbf{a}_l^{H}\textbf{U}_s\Lambda_s\textbf{U}_s^H\textbf{a}_l +\hat{\sigma}_{nl}^2\left(M-\|\textbf{U}_s^H\textbf{a}_l\|_2^2\right)$

因此根据第 $l$ 个角度的名义导向矢量获得的噪声功率估计值为：

$\hat{\sigma}^2_{nl}=\dfrac{\textbf{a}^H_I\hat{\textbf{R}}_X\textbf{a}_I-\textbf{a}^H_I\textbf{U}_s\Lambda_S\textbf{U}_s^H\textbf{a}_l }{M-\|\textbf{U}_s^H\textbf{a}_l \|^2_2}$

对所有获得的噪声功率估计值求平均即为估计得到的噪声功率，即：

$\hat{\sigma}_n^2=\frac{1}{L}\sum_{l=1}^L\hat{\sigma}_{nl}^2$

对于信源及干扰的功率估计，作者引入了以下方法：

一般情况下，接收信号的结构为：

$\textbf{x}(k)=\textbf{A}\textbf{s}(k)+\textbf{n}(k)$

对其进行自相关即得到自协方差阵 $\mathbf{R}_x$ 。对自协方差阵左乘阵列流形 $\mathbf{A}$ 的转置再右乘阵列流形 $\mathbf{A}$ ,即为：

$\textbf{A}^H\textbf{R}_x\textbf{A}=\textbf{A}^H\textbf{A}\textbf{C}_{s+i}\textbf{A}^H\textbf{A}+\sigma_n^2\textbf{A}^H\textbf{A}$

其中 $\textbf{C}_{s+i}$ 即为表示不同方向来向信号的功率值，不过一般通过对 $\mathbf{R}_x$ 进行特征分解，其特征值对角阵即为 $\textbf{C}_{s+i}+\sigma_n^2 \textbf{I}$ 。

因此作者通过以下等式求解 $\textbf{C}_{s+i}$ ：

$\hat{\textbf{C}}_{s+i}=(\overline{\textbf{A}}^H\overline{\textbf{A}})^{-1}\left(\overline{\textbf{A}}^H\hat{\textbf{R}}_x\overline{\textbf{A}}-\hat{\sigma}_n^2\overline{\textbf{A}}^H\overline{\textbf{A}}\right)(\overline{\textbf{A}}^H\overline{\textbf{A}})^{-1}$

协方差矩阵重建

作者提出的协方差矩阵重建法如下式所示：

$\widetilde{\textbf{R}}=\overline{\textbf{A}}\hat{\textbf{C}}_{s+i}\overline{\textbf{A}}^H+\hat{\sigma}_n^2\textbf{I}+\overline{\textbf{A}}_i\hat{\textbf{D}}_i\overline{\textbf{A}}_i^H=\overline{\textbf{A}}\tilde{\textbf{C}}_{s+i}\overline{\textbf{A}}^H+\hat{\sigma}_n^2\textbf{I}$

其中：

$\overline{\textbf{A}}_i=\left[\hat{\textbf{a}}_1,\hat{\textbf{a}}_2,\cdots,\hat{\textbf{a}}_Q\right]$

为各个干扰的估计导向矢量构成的阵列流形，

$\hat{\mathbf{D}}_i=\text{diag}\left\{\alpha_1\hat{\sigma}_1^2,\alpha_2\hat{\sigma}_2^2,\cdots,\alpha_Q\hat{\sigma}_Q^2\right\}$

为一对角加载阵，其中 $\alpha_{q},q=1,\text{}2,\cdots,Q$ 即为预先设定的对角加载量。

对角加载法原先通过增加协方差阵的特征值幅度，减小因噪声子空间最大特征值与最小特征值之比(一般称为白噪声特征值扩散程度)过大而带来的波束旁瓣提升的影响，从另一种方面来看，该种方法即通过间接降低信号信噪比，从而提高阵列的稳健性。

作者此处通过添加对角阵，使得协方差阵在经过对角加载后的干扰信号比：

$\text{ISR}=\dfrac{(1+\alpha_q)\hat{\sigma}_q^2}{\hat{\sigma}_0^2}$

始终高于一定的值，即：若第 $q$ 个干扰信号比低于于预先设定的阈值，则通过上式计算该干扰信号的对角比 $\alpha_{q}$ ，使得协方差阵中该信号的干扰信号比提高到一定的值。该种方法能够一定程度上解决高信噪比下信号的自消效应，同时也起到了加深波束图上零陷的深度的作用，对估计得到的导向矢量周围的Capon谱幅度进行了进一步抑制。(但是自消效应的前提是建立在信号因噪声与有限快拍数导致的协方差矩阵失真的情况下，干扰加噪声协方差阵中含有期望信号分量，且当估计得到的导向矢量与真实的导向矢量的二范数误差大于一定程度时，干扰信号比的提升将失去作用)

最后，得到的权向量可表示为：

$\mathbf{W}_{\text{p}}=\frac{\widetilde{\mathbf{R}}^{-1}\hat{\mathbf{a}}_0}{\hat{\mathbf{a}}_0^H\widetilde{\mathbf{R}}^{-1}\hat{\mathbf{a}}_0}$

仿真

部分仿真条件如表所示：

阵元数	10
信噪比	0
扰噪比	20
阵元间距	半波长
信源数	1
干扰数	2
信源到达角	5°
干扰到达角	-30°，45°
信源角度估计误差	[0°,3°]
角度采样步长	0.1°
随机采样点数L	10
角度搜索范围	4°
ISR门限	8dB

代码如下所示：

clc;
clear;
jk = 0;
loop = 50;
SINRO = zeros(1,length(-20:5:40));
snr = 0;
%% 初始化及设定参数
array_num = 10;%%阵元数
snapshot_num = 100;%%快拍数
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_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误差
source_dev = source_aoa+[ang_mismatch1,ang_mismatch2,ang_mismatch3];%%实际获得的先验角度
sample_acc = 0.1;%%信号域及干扰域离散采样步长
search_range = 4;
%% 转向矢量
A_bar = exp(-1i*(0:array_num-1)'*pi*sind((source_dev)));%%实际带有误差的阵列响应矢量
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;%%接收阵的协方差矩阵
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;
%% 算法
A_est = zeros(array_num,source_num);
for ik = 1:length(source_dev)
    max_capon_spec_val = 0;
    for jk = source_dev(ik)-search_range:sample_acc:source_dev(ik)+search_range
        current_sv = exp(-1i*(0:array_num-1)'*pi*sind(jk));
        if 1/(current_sv'*inv(R)*current_sv) > max_capon_spec_val
            max_capon_spec_val = 1/(current_sv'*inv(R)*current_sv);
            max_deg = jk;
            A_est(:,ik) = current_sv;
        end
    end
    disp(max_deg);
end
[U,D,V] = svd(R);
Us = U(:,1:source_num);
Lambda_s = D(1:source_num,1:source_num);
sample_L_dot = 10;
L = -90:180/sample_L_dot:90-180/sample_L_dot;
noise_est = zeros(1,sample_L_dot);
for ik = 1:length(L)
    A_l = exp(-1i*(0:array_num-1)'*pi*sind(L(ik)));
    noise_est(ik) = (A_l'*R*A_l-A_l'*Us*Lambda_s*Us'*A_l)/(array_num-norm(Us'*A_l,2)^2);
end
noise_est = mean(abs(noise_est));
C_si = inv(A_est'*A_est)*A_est'*(R-noise_est*eye(array_num,array_num))*A_est*inv(A_est'*A_est);
R_est = A_est*C_si*A_est'+noise_est*eye(array_num,array_num);
ISR = 8;
D = zeros(source_num-1,source_num-1);
for ik = 2:source_num
    if 10*log10(C_si(ik,ik)/C_si(1,1)) < ISR
        D(ik-1,ik-1) = C_si(1,1)/C_si(ik,ik)*10^(ISR/10)-1;
    end
end
R_est = R_est+A_est(:,2:end)*D*A_est(:,2:end)';
w0 = inv(R)*A_est(:,1)/(A_est(:,1)'*inv(R)*A_est(:,1));
beam_plot(w0);

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,'LineWidth',2);
end

仿真运行结果如下所示：

可以发现，对干扰使用对角加载法，虽然有效地进一步抑制了噪声，其缺点便是波束图旁瓣增益不稳定。

总的来说，这篇文章对噪声的估计方法是比较新颖的。

参考文献

Yue Yang, Xu Xu, Huichao Yang, Wangjie Li. Robust adaptive beamforming via covariance matrix reconstruction with diagonal loading on interference sources covariance matrix, Digital Signal Processing. Volume 136, 2023, 103977.

量大龙神魔

关注

1
点赞
踩
23

收藏

觉得还不错? 一键收藏
1
评论
波束形成基于协方差矩阵重建与对角加载的稳健自适应波束形成算法

文献 Robust adaptive beamforming via covariance matrix reconstruction with diagonal loading on interference sources covariance matrix-DSP-2023
复制链接

扫一扫