波束形成基于SOCP及QCQP的稳健自适应波束形成

量大龙神魔

已于 2024-03-17 22:03:45 修改

阅读量772

点赞数 5

文章标签：算法

于 2023-06-05 16:51:58 首次发布

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

版权

前言

新世纪以来，随着凸优化理论兴起，一些基于凸优化或者非凸优化凸近似的算法发展迅猛，比较著名的便是Sergiy A. Vorobyov老师提出的最坏情况性能最优波束形成算法以及黄永伟老师在TSP上发表的一系列的非凸问题凸近似的波束形成算法，而对于导向矢量估计，Capon功率谱的特性使得凸优化算法可以很好地应用，比较有代表性的便是QCQP与SOCP。读到了这篇文章，在此做一下简单的记录。

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

导向矢量估计

SOCP形式下的导向矢量估计问题可表示为：

$\min\limits_{\mathbf{\hat{a}}(\theta)}\mathbf{\hat{a}}^H(\theta)\mathbf{\hat{R}}^{-1}\mathbf{\hat{a}}(\theta),s.t.\left \| \mathbf{\hat{a}}(\theta)-\mathbf{{a}}(\theta) \right \|\leq \varepsilon$

其中 $\mathbf{\hat{R}}^{-1}$ 在理想情况下为一半正定阵，在带噪声情况下为一正定阵。该问题为凸，可通过CVX求解。则可重新定义在角度 $\theta$ 时，该阵列对应的导向矢量为 $\mathbf{\hat{a}}(\theta)$ 。(这种方法固然考虑了阵列失配对入射波导向矢量的影响，但是取该锥邻域内的Capon功率谱最大值会产生其他影响，若 $\varepsilon$ 过大，则会使得重建协方差阵的Capon功率谱谱峰变宽，若 $\varepsilon$ 过小，则会使得阵列失配问题考虑的不够充分)

对于第 $j$ 个入射波的导向矢量估计，作者采用在该先验角度领域内线积分重构协方差阵取其最大特征值对应的特征向量作为估计得到的导向矢量，有：

$\begin{aligned}\hat{\textbf{R}}_j&=\int\limits_{\bar{\Theta}_j}\hat{P}(\theta)\hat{\textbf{a}}(\theta)\hat{\textbf{a}}^H(\theta)d\theta=\int\limits_{\bar{\Theta}_j}\frac{\hat{\textbf{a}}(\theta)\hat{\textbf{a}}^H(\theta)}{\hat{\textbf{a}}^H(\theta)\hat{\textbf{R}}^{-1}\hat{\textbf{a}}(\theta)}d\theta\end{aligned}$

对该协方差阵进行特征分解，对其最大特征值进行幅值修正 $\mathbf{\hat{a}}_j = \sqrt{M} \mathbf{j}_1$ ，即为估计得到的导向矢量。这种方法较为不错地解决了模型中同时存在非相干散射与阵列失配的问题。当存在非相干散射时，上式协方差阵中的元素组成可表示为：

$\textbf{R}_j=\sigma_j^2\tilde{\textbf{a}}(\theta_j)\tilde{\textbf{a}}^H(\theta_j)+\sum_{p=1}^P\sigma_{j_p}^2\tilde{\textbf{a}}(\theta_{j_p})\tilde{\textbf{a}}^H(\theta_{j_p})+\sigma_n^2\mathbf{I}$

当阵列的分辨率足以分辨出散射角对应的不同导向矢量时，就需要通过入射波数量的先验信息或者一系列信源数估计算法进行导向矢量的选择，即是否需要接收信号的非相干散射分量。

对于期望信号的导向矢量，作者首先在期望信号先验角度邻域内进行协方差阵积分重建，并对其进行特征值分解，即：

$\hat{\textbf{R}}_\textsf{s}=\int\limits_{\Theta}\frac{\hat{\textbf{a}}(\theta)\hat{\textbf{a}}^H(\theta)}{\hat{\textbf{a}}^H(\theta)\hat{\textbf{R}}^{-1}\hat{\textbf{a}}(\theta)}d\theta=\sum_{m=1}^M v_m\textbf{d}_m\textbf{d}_m^H$

之后，求该矩阵的信号子空间的噪声子空间。当存在散射时，信号子空间的维数即为直射路径数加散射路径数，当不存在散射时，信号子空间的维数为1。定义信号子空间特征向量为：

$\textbf{D}\triangleq[\textbf{d}_1,\textbf{d}_2,...,\textbf{d}_R]$

则信号子空间的正交子空间(一说为噪声子空间)为： $\textbf{P}_{\textbf{s}}^{\perp}\triangleq\textbf{I}-\textbf{D}\textbf{D}^H$

实际导向矢量 $\tilde{\textbf{a}}_0$ 与名义导向矢量 $\textbf{a}_0$ 的关系为： $\tilde{\textbf{a}}_0=\textbf{a}_0+\textbf{e}$ ，其中 $\textbf{e}$ 为误差分量，而其中对Capon功率谱幅度起作用的则是误差分量中的正交分量，因此有以下QCQP问题可估计误差分量中的正交分量：

$\underset{\textbf{e}_\perp}{\text{min}}(\textbf{a}_0+\textbf{e}_\perp)^H\hat{\textbf{R}}^{-1}(\textbf{a}_0+\textbf{e}_\perp)\\ s.t.\textbf{P}_\textbf{s}^\perp(\textbf{a}_0+\textbf{e}_\perp)=\textbf{0},\quad\textbf{a}_0^H\textbf{e}_\perp=~0,\quad(\textbf{a}_0+~\textbf{e}_\perp)^H\hat{\textbf{R}}_{i+n}(\textbf{a}_0+~\textbf{e}_\perp)\leq~\textbf{a}_0^H\hat{\textbf{R}}_{i+n}\textbf{a}_0$

该优化问题中的三个约束项，第一个是保证估计得到的导向矢量属于信号子空间，第二个是保证正交分量与导向矢量正交，第三个是保证估计得到的导向矢量其Capon功率谱幅度更优。

则估计得到的导向矢量为： $\hat{\textbf{a}}_{\text{es}}=\textbf{a}_{0}+\textbf{e}_{\perp}$

噪声、信源功率估计

对于信源功率估计，作者采用Capon功率谱谱峰作为入射波功率估计值。

对于噪声功率估计，作者采用协方差阵噪声子空间对应的特征值均值作为噪声功率估计值。

协方差矩阵重建

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

$\hat{\mathbf{R}}_{i+n}=\hat{\mathbf{R}}_{\mathrm{int}}+\hat{\sigma}_{n}^{2}\mathbf{I}=\sum\limits_{j=1}^J(\hat{\sigma}_j^2\hat{\textbf{a}}_j\hat{\textbf{a}}_j^H+\sum\limits_{p=2}^{P+1}\hat{\sigma}_{jp}^2\hat{\textbf{a}}_{jp}\hat{\textbf{a}}_{jp}^H)$

仿真

部分仿真条件如表所示：

阵元数	10
信噪比	0
扰噪比	20
阵元间距	半波长
信源数	1
干扰数	2
信源到达角	5°
干扰到达角	-30°，45°
信源角度估计误差	[0°,3°]
角度采样步长	0.5°
噪声角度域采样点数	10
角度搜索范围	8°
SOCP问题约束项范数上确界	0.05
正交向量步长	0.001

代码如下所示：

%% 初始化及设定参数
clc;
clear;
jk = 0;
snr = 10;
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];%%信噪比、扰噪比
epsl = 0.05;
search_range = 8;
deg_grid = 0.5;
%% 导向矢量
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;
%% SOCP重建
inv_R = inv(R);
A_est = zeros(array_num,source_num);
pow_est = zeros(1,source_num);
R_source = zeros(array_num,array_num);
for ik = 1:source_num
    deg = source_dev(ik)-search_range:deg_grid:source_dev(ik)+deg_grid;
    sv_recons = zeros(array_num,length(deg));
    for jk = 1:length(deg)
        sv_nominal = exp(-1i*(0:array_num-1)'*pi*sind(deg(jk)));
        cvx_begin quiet
            variable a_est_grid(array_num,1) complex
            minimize(quad_form(a_est_grid,inv_R));
            subject to 
                norm(a_est_grid-sv_nominal) <= epsl;
        cvx_end
        sv_recons(:,jk) = a_est_grid;
    end
    R_j = zeros(array_num,array_num);
    for jk = 1:length(deg)
        R_j = R_j+sv_recons(:,jk)*sv_recons(:,jk)'/(sv_recons(:,jk)'*inv_R*sv_recons(:,jk));
    end
    if ik == 1
        R_source = R_j;
    end
    [v,d] = eig(R_j);
    [~,d_index] = sort(diag(d),"descend");
    v = v(:,d_index);
    A_est(:,ik) = sqrt(array_num)*v(:,1);
    pow_est(ik) = 1/(A_est(:,ik)'*inv_R*A_est(:,ik));
end
%% 协方差阵重建
R_recons = zeros(array_num,array_num);
for ik = 2:source_num
    R_recons = R_recons+pow_est(ik)*A_est(:,ik)*A_est(:,ik)';
end
[~,d] = eig(R);
diag_d = sort(diag(d),'descend');
pow_n = mean(diag_d(source_num+1:end));
R_recons = R_recons+pow_n*eye(array_num,array_num);
%% QCQP sv est
R_s = R_source;
[v,d,~] = svd(R_s);
[diag_d,d_index] = sort(diag(d),'descend');
v = v(:,d_index);
%%非相干散射情况下，这里我就直接拿噪声子空间上去了，因为信号子空间噪声子空间之和为单位阵
Ps = v(:,2:end)*v(:,2:end)';
a_0 = exp(-1i*(0:array_num-1)'*pi*sind(source_dev(1)));
cvx_begin
    variable e_prep(array_num,1) complex
    minimize(quad_form(a_0+e_prep,inv_R));
    subject to
        Ps*(a_0+e_prep) == 0;
        a_0'*e_prep == 0;
        quad_form(a_0+e_prep,R_recons) <= quad_form(a_0,R_recons);
cvx_end
a_est_0 = a_0+e_prep;
% a_est_0 = A_est(:,1);
w0 = inv(R_recons)*a_est_0/(a_est_0'*inv(R_recons)*a_est_0);
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

仿真结果运行如下所示：

作者的本意是想通过cvx进行导向矢量的精确求解，但是显然作者没考虑该算法的计算量，通过cvx的内点法及频繁的矩阵特征分解操作都大大增加了计算量，再一个就是CVX求解SOCP时，约束项过小就会使得寻找期望解失败。这篇文章的仿真让本人很好地回顾了下cvx的用法。

参考文献

Haoran Li, Jun Geng, Junhao Xie. Robust adaptive beamforming based on covariance matrix reconstruction with RCB principle. Digital Signal Processing, Volume 127, 2022, 103565.

量大龙神魔

关注

5
点赞
踩
12

收藏

觉得还不错? 一键收藏
2
评论
波束形成基于SOCP及QCQP的稳健自适应波束形成

文献 Robust adaptive beamforming based on covariance matrix reconstruction with RCB principle-DSP-2022
复制链接

扫一扫