波束形成基于导向矢量估计与干扰协方差阵重建的稳健自适应波束形成算法

量大龙神魔

已于 2023-06-15 15:18:30 修改

阅读量766

点赞数 5

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

于 2023-06-03 21:12:41 首次发布

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

版权

前言

这篇文章是我接触的第一篇的关于协方差阵重建波束形成算法的文章，是中科大的同一个课题组基于2015年TSP的体积分的工作，提出的导向矢量估计重建协方差矩阵的思想，在之前仿真过程中碰到了不少问题，导致波束形成器实际性能不达预期，在此进行简单记录。

整个算法将分成导向矢量估计、协方差阵重建与噪声、信源功率估计三部分。

P.S.这篇文章的作者本人(孙思聪)的硕士论文第三章即为这篇文章，如果对这篇博客或原英文文献不清楚可前往知网自行查找。

导向矢量估计

作者提出的导向矢量估计法从某一角度来说可算作是网格法求解最优导向矢量，只不过该网格是由名义导向矢量中的角度以及导向矢量的超平面的正交向量的模联合确定的。空域功率谱，在DOA方面可通过角度-功率进行表示，可视为一个受角度影响的一维频谱。在波束形成方面，又可通过导向矢量-功率进行表示，此时导向矢量中的每个元素的变化皆可影响其幅度，因此又可视为一个受导向矢量影响的高维频谱。

对于正交向量的确定，作者采用在第 $l$ 个入射波的先验角度周围一定范围的邻域 $[\hat{\theta}_l-\frac{\Delta\theta_l}{2},\hat{\theta}_l+\frac{\Delta\theta_l}{2}]$ 内进行均匀离散选点，若有 $M$ 个阵元则选取 $M-1$ 个点，保证该阵列流形矩阵的列欠秩性。则所得到的阵列流形矩阵可表示为：

$\textbf{A}_l=[\textbf{a}(\hat{\theta}_{l1}),\textbf{a}(\hat{\theta}_{l2}),...,\textbf{a}(\hat{\theta}_{l(M-1)})]$

根据SVD，有 $\textbf{A}_l = UDV$ ，则若想得到 $\textbf{A}_l$ 列空间的正交向量，只需 $\textbf{A}_l\textbf{A}_l^H = UDVV^HD^HU^H$ 。则可对矩阵 $\textbf{A}_l\textbf{A}_l^H$ 作特征分解，有：

$\textbf{A}_l\textbf{A}_l^H = \textbf{A}_l\textbf{A}_l^H=\sum\limits_{m=1}^M\lambda_{lm}\textbf{e}_{lm}\textbf{e}_{lm}^H=\textbf{E}_{l\text{h}}\Lambda_{l\text{h}}\textbf{E}_{l\text{h}}^H+\lambda_{lM}\textbf{e}_{lM}\textbf{e}_{l\text{M}}^H$

最小特征值对应的特征向量即为正交向量，在得到正交向量后，最优导向矢量的求解问题变为：

$\underset{\tilde{\textbf{a}}}{\text{max}}P(\tilde{\textbf{a}}_l)=\frac{1}{\tilde{\textbf{a}}_l^H\hat{\textbf{R}}^{-1}\tilde{\textbf{a}}_l},s.t.\tilde{\textbf{a}}_l=\sqrt{M}\frac{\hat{\textbf{a}}_l+\eta\textbf{e}_{lM}}{||\hat{\textbf{a}}_l+\eta\textbf{e}_{lM}||_2},-\tilde{\eta} \le \eta \le \tilde{\eta}$

该问题非凸，因此只能采用网格法遍历求解。

或者，该问题可改写为

$\underset{\tilde{\textbf{a}}}{\text{min}}~\tilde{\textbf{a}}_l^H\hat{\textbf{R}}^{-1}\tilde{\textbf{a}}_l,s.t.\left \| \tilde{\textbf{a}}_l-\bar{\textbf{a}}_l \right \|_2 \le \left \| \bar{\textbf{a}}_l-{\textbf{a}}(\hat {\theta}_k+\Delta\theta/2) \right \|_2$

该问题为一服从球形约束的QCQP凸问题，可通过CVX求解。

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

考虑单信号的协方差阵的组成结构，为：

$\textbf{R}=\sigma_l^2\textbf{a}_l\textbf{a}_l^H+\sigma_n^2\textbf{I}$

则该协方差阵的Capon功率谱可表示为：

$P(\theta)=\frac{1}{\textbf{a}^H(\theta)(\sigma_l^2\textbf{a}_l\textbf{a}_l^H+\sigma_n^2\textbf{I})^{-1}\textbf{a}(\theta)}$

此时首先通过WoodBury矩阵公式处理分母中的逆矩阵，有：

$(\sigma_l^2\textbf{a}_l\textbf{a}_l^H+\sigma_n^2\textbf{I})^{-1} = \sigma_n^{-2}\textbf{I}-\sigma_n^{-4}\textbf{a}_l(\sigma_l^{-2}+\sigma_n^{-2}\textbf{a}_l^H\textbf{a}_l)^{-1}\textbf{a}_l^H\\~~~~~~~~~~~~~~~~~~~~~~~=\sigma_n^{-2}\textbf{I}-\frac{\sigma_n^{-4}}{\sigma_l^{-2}+\sigma_n^{-2}M_l}\textbf{a}_l\textbf{a}_l^H$

且已知两个复向量的内积为：

$\textbf{u}^H\textbf{v}=conj(\textbf{v}^H\textbf{u})=||\textbf{u}||_2||\textbf{v}||_2\cos\theta_H(\textbf{u},\textbf{v})e^{j\varphi},\textbf{0}\le~\theta_H(\textbf{u},\textbf{v})\le~\frac{\pi}{2},-\pi\le~\varphi\le~\pi$

其中 $\theta_H(\textbf{u},\textbf{v})$ 为厄密特角， $e^{j\varphi}$ 为伪角度。

将展开的逆矩阵代回Capon功率谱表达式中，可得：

$P(\theta)=\frac{1}{\textbf{a}^H(\theta)(\sigma_n^{-2}\textbf{I}-\frac{\sigma_n^{-4}}{\sigma_l^{-2}+\sigma_n^{-2}M_l}\textbf{a}_l\textbf{a}_l^H)\textbf{a}(\theta)}=~\frac{1}{\sigma_n^{-2}M-\frac{\sigma_n^{-4}}{\sigma_l^{-2}+\sigma_n^{-2}M_l}\cos^2\theta_H[\textbf{a}(\theta),\textbf{a}_l]MM_l}\\~~~~~~~=\frac{\sigma_l^{-2}+\sigma_n^{-2}M_l}{\sigma_n^{-2}\sigma_l^{-2}M+\sigma_n^{-4}MM_l-\sigma_n^{-4}\cos^2\theta_H[\textbf{a}(\theta),\textbf{a}_l]MM_l}=\frac{\sigma_\text{n}^2+\sigma_l^2M_l}{M+MM_lr(1-\rho^2)}$

其中 $\rho=\cos\theta_H[\textbf{a}(\theta),\textbf{a}_l]$ ， $r=\sigma_l^2/\sigma_n^2$ ， $\sqrt{M_l}=\left \| \mathbf{a}_l \right \|_2$ ， $\sqrt{M}=\left \| \mathbf{a}(\theta ) \right \|_2$ 。当两导向矢量近乎不相关时，可得 $\rho^2\approx 0$ ，则有： $P(\theta) = \sigma_n^2/M$ ，即可据此估计噪声。

在获得噪声功率与第 $l$ 个入射波的导向矢量 $\mathbf{\hat {\tilde {a}}}_l$ 后，即可得到其功率为：

$\sigma_l^2 =P\left (\mathbf{\hat {\tilde {a}}}_l \right )-\sigma_0^2$

干扰加协方差矩阵重建

重建得到的干扰加噪声协方差阵可表示为：

$\hat{\mathbf{R}}_{\text{i+n}}=\sum\limits_{l=1}^L\sigma_l^2\hat{\tilde{\mathbf{a}}}_l\hat{\tilde{\mathbf{a}}}_l^H+\hat{\sigma}_{\text{n}}^2\mathbf{I}$

仿真

部分仿真条件如表所示：

阵元数	10
信噪比	0
扰噪比	20
阵元间距	半波长
信源数	1
干扰数	2
信源到达角	5°
干扰到达角	-30°，45°
信源角度估计误差	[0°,3°]
角度采样步长	0.1°
噪声角度域采样点数	10
角度搜索范围	4°
正交向量模边界	0.1
正交向量步长	0.001

代码如下所示：

%% 初始化及设定参数
clc;
clear;
jk = 0;
loop = 50;
SINRO = zeros(1,length(-20:5:40));
snr = 10;
% for kkk = 1:100
%     disp(kkk);
array_num = 10;%%阵元数
deg_dev_predf = 4;
sample_acc = 0.1;
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误差
snapshot_num = 30;%%快拍数
source_aoa = [5,-30,45];%%信源到达角
tgt_ang = 1;
source_dev = source_aoa + [ang_mismatch1,ang_mismatch2,ang_mismatch3];%%DOA得到的带有误差的
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];%%信噪比、扰噪比
search_range = 4;
points = array_num-1;%%%%%%%%
eta_predf = 0.1;
eta_acc = 1e-3;
%% 导向矢量
A_bar = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*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(3));
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;
%% 噪声功率估计
source_dev1 = sort(source_dev,'ascend');
noise_deg = -90:sample_acc:source_dev1(1)-search_range;
for ik = 2:source_num
noise_deg = [noise_deg source_dev1(ik-1)+search_range:sample_acc:source_dev1(ik)-search_range];
end
noise_deg = [noise_deg source_dev1(end):sample_acc:90];
flag = floor(linspace(1,length(noise_deg),array_num));
pow_n = 0;
for ik = 1:length(flag)
    aa = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(noise_deg(ik)));
    pow_n = pow_n+1/(aa'*inv(R)*aa);
end
A_est = zeros(array_num,source_num);
pow_est = zeros(1,source_num);
for ik = 1:source_num
%%正交向量求解
    sample_dot = linspace(source_dev(ik)-search_range,source_dev(ik)+search_range,points);
    A = exp(-1i*(0:array_num-1)'*pi*sind(sample_dot));
    [V,D] = eig(A*A');
    [diag_d,d_index] = sort(D,'descend');
    V = V(:,d_index);
    e_otho = V(:,end);
    axis_deg = source_dev(ik)-search_range:sample_acc:source_dev(ik)+search_range;
    axis_e_vec = -1*eta_predf:eta_acc:eta_predf;
    capon_value= zeros(length(axis_deg),length(axis_e_vec));
    for jk = 1:length(axis_deg)
        for kk = 1:length(axis_e_vec)
            normalization_sv = (exp(-1i*(0:array_num-1)'*pi*sind(axis_deg(jk)))+axis_e_vec(kk)*e_otho)/norm((exp(-1i*(0:array_num-1)'*pi*sind(axis_deg(jk)))+axis_e_vec(kk)*e_otho),2);
            capon_value(jk,kk) = 1/(normalization_sv'*inv(R)*normalization_sv);
        end
    end
    max_val = max(max(capon_value));
    [row,col] = find(max_val == capon_value);
    A_est(:,ik) = (exp(-1i*(0:array_num-1)'*pi*sind(axis_deg(row)))+axis_e_vec(col)*e_otho)/norm((exp(-1i*(0:array_num-1)'*pi*sind(axis_deg(row)))+axis_e_vec(col)*e_otho),2);
    pow_est(ik) = max_val-pow_n;
end
R_recons = zeros(array_num,array_num);
for ik = 1:source_num
    R_recons = R_recons+pow_est(ik)*A_est(:,ik)*A_est(:,ik)';
end
R_recons = R_recons+pow_n*eye(array_num);
w0 = inv(R_recons)*A_est(:,1)/(A_est(:,1)'*inv(R_recons)*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

仿真结果运行如下所示：

对导向矢量的非凸估计是一个典型的On-Grid的做法，精度受网格步长限制，可以发现，在小步长的前提下算法有着不错的性能。

参考文献

Sicong Sun, Zhongfu Ye. Robust adaptive beamforming based on a method for steering vector estimation and interference covariance matrix reconstruction. Signal Processing. Volume 182, 2021, 107939.

量大龙神魔

关注

5
点赞
踩
14

收藏

觉得还不错? 一键收藏
2
评论
波束形成基于导向矢量估计与干扰协方差阵重建的稳健自适应波束形成算法

文献 Robust adaptive beamforming based on a method for steering vector estimation and interference covariance matrix reconstruction-SP-2021
复制链接

扫一扫