本文主要是介绍波束形成 基于不等式惩罚项约束与自由度限制的稳健自适应波束形成算法,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
前言
这篇文章是一篇典型的基于凸优化的稳健自适应波束形成算法。作者的核心思想主要是通过惩罚约束自由度与不等式项约束不确定域从而达到阵列的稳健性,同时考虑了导向矢量失配、阵列自由度限制与协方差矩阵失真三个因素。虽然很早就打算尝试对这篇文章进行仿真,但是囿于能力问题还是没能成功,现在回过头来打算再仔细看看,对作者的内容进行补充。
阵列稳健性之导向矢量适配
针对导向矢量失配问题,作者提出了以下约束:
定义1
给定先验角度,与该角度对应的真实导向矢量,则有以下约束成立:
其中为人为定义的误差阈值,为针对导向矢量失配的稳健性系数。已知在导向矢量失配问题存在的情况下,真实导向矢量可分解为其名义分量与阵列误差分量,即:
其中为其名义导向矢量分量,为阵列误差分量,且有成立。
MVDR的等价优化形式可表示为:
对于该最优解,与邻域内的导向矢量,有:
其中为一预先定义的不确定角度域。
阵列稳健性之协方差矩阵失真
针对小快拍数下协方差矩阵的失真问题,常见的方法是根据有限的快拍数进行协方差矩阵重建,尽量得到无失真情况下的协方差矩阵;另一种方法则是根据先验角度信息对干扰域邻域范围内的导向矢量都施加约束以进行抑制,避开协方差矩阵的进一步处理过程;此处作者提出的方法明显是属于后一种。
定义2
对干扰域角度的约束可表示为:
其中为一根据第个角度的先验信息而定义的离散角度集合,同时表示着其抑制程度,同样的,对于该阵列权向量的最优解,与邻域内的导向矢量,有:
可发现,作者采取的是对波束形成器的空域响应施加约束以提高阵列稳健性。
阵列稳健性之有限自由度分配
针对阵列有限的自由度及干扰的抑制,作者在干扰域的约束问题上进一步添加了惩罚项,通过降低干扰域空域响应的最大值进行整体抑制,即:
其中为一预先指定的干扰抑制参数,更大的表明该干扰源更应优先抑制,故应与干扰源功率呈正相关的关系;为一控制空域响应的参数,在的控制与max的约束下,大的总会有小的与之对应,随之而来的便是对应方向上更为严格的空域响应约束;当足够大时,便总能存在一最优解。(私以为该部分与其说是自由度分配不如说是空域谱控制)
P-ICMV的优化问题约束形式
考虑了以上三种情况下不同的优化问题与约束问题,作者提出的波束形成器可表示为:
定义3
令为一由期望信号邻域(该邻域应足够小)内的离散角度点的集合中的元素对应的名义导向矢量拼接而成的阵列流型,假定的列秩为该集合中的元素个数,(此处可得采样点数显然小于阵元数)则若,P-ICMV的优化问题总是可行的。
证明如下:
考虑对于一权向量,其相对于阵列流型而得到空域响应复向量,则可通过该空域响应复向量,用最小二乘法得到的最小二乘解,有:
带入约束问题中的第一个约束项,有
令,即得
该波束形成器还具有以下性质:
- 有限空域响应
针对的情况,即误差在允许的范围内时,针对P-ICMV波束形成器的可行解,其期望信号邻域内角度的阵列响应存在一个上确界与下确界,即:
同样的,干扰信号邻域内角度的阵列响应也存在一个上确界:
若进一步考虑导向矢量失配的情况,则有,或
- 自由度分配(过载情况下信号抑制的优先级)
P-ICMV波束形成器的优化项中,在的控制与max的约束下使得波束形成器能够根据提供的对对应方向上的信号进行进一步约束。
P-ICMV的优化问题约束形式对应的ADMM解
因P-ICMV的优化问题中存在一个min-max的优化项难以处理,因此作者在此进行了一个等价变换:
显然此为一松弛操作,考虑:
该优化问题的最优值显然由确定的前提下,最小化序列最大值的过程中最先到达下确界的决定,当某一个达到下确界时最优值即给出。(e.g.)
因此,以下优化问题的可行域显然比min-max问题的可行域更为严格。
同时,进一步引入辅助变量与,有:
则原问题可进一步转化为:
则可得到该问题的增广拉格朗日式为:
需要注意的是,复值优化问题,乘子与约束的内积是需要取实部的。(第一遍看这个问题的ADMM解法的时候本人也感到疑惑,到后面才发现作者是将约束15.a与约束15.b带到后面的拉格朗日乘子最小的环节,使之从一个无约束问题变化为了一个带约束问题,个人感觉这块处理得非常精妙。(可能是本人知识浅薄))
该增广拉格朗日乘子中各变量的迭代式为:
其中,,,
定义4
在定义3成立的前提下,ADMM各变量的迭代式会收敛至最优解。
ADMM中第一个变量迭代式的解
除去增广拉格朗日函数中无关的分量后,求解关于的最优值的优化问题如下所示:
其中:
则该优化问题的最优解为:
其中,。
过程如下:
从的等式可得,该矩阵为一正定阵,因此首先对该矩阵作特征分解,有:;根据正交矩阵的性质,有:,则令,故原问题可表示为:
显然该问题的最优解与关于的最优值的优化问题的最优解可相互转换,则其拉格朗日乘子如下:
又因为,注意到该式在处不可导,则需对可行域进行进一步划分。
定义6
关于的最优值的优化问题的等价问题的最优解为,当且仅当。
证明如下:
已知,令,则有:
假定,定义一可行解为,有:
则。
因此为了避开最优解为的情况,首先考虑的条件。在下该问题的KKT条件为:
当时,根据4.可得:,根据5.可得:,当时1.成立。
当时,根据3.可得,根据5.可得,带入4.即得:
从上式可得:
且为单调递减。且进一步考虑与,有:
且有:
进一步验证了。同样地,考虑,则有:
即确定了,因此最优解即可得到。
ADMM中第二个变量迭代式的解
该优化问题是可分别对每对分别求解得到,该优化问题可表示为:
其中,,。该系列的复SOCP优化问题可通过以下引理进行求解:
引理
考虑一个复SOCP优化问题为:
其中,,,,令,,其最优解为:
证明如下:
由约束项得:,首先求解,对于任何在约束项内的,有:
带入,显然约束项成立,因此slater条件成立,注意到为复数,且该问题优化函数为强凸函数,对其求偏导有点难度,因此可根据本人在前一文章中的KKT条件部分,得该凸问题的KKT条件为:
由此可发现若,则,不可确定,因此首先考虑时的情况。在的前提下,根据4.,有:
当时,根据3.,有:,其中为转向角,带入,有:
根据2.,有:,因此可得。
当时,,带入1.,有:,因此的最优解为:
将的情况包括进去,令(复数在复坐标轴上的角度),,则有:
即得 。
对于的求解,将带回原式,可得:
其中:
显然连续且强凸,因此存在两种情况:1.,2.。若,则在时单调递减,即得。若,则。若非以上两种情况,则:
因此有:
其中。注意到这三种情况下都是其最小值,因此有:
即得。
ADMM中第三个变量迭代式的解
该优化问题可表示为:
其中
将该优化问题拆解为两部分,第一部分为:
根据上一引理,其最优解为:
第二部分为一无约束问题:
其中,且令,有:
则该问题的最优解可表示为:
其过程如下:
首先对进行求导,有:
显然在上是单调递增且连续的,且在上有,根据一阶最优性条件:
因此:
若,则在该可行域内的最优点为,有:
若,则最优解一定在可行域内,因此将与中的元素重新以递增的形式拼接成一个集合,当且时即得到最优解,因此有:
其中,。
仿真
囿于时间原因,打算通过对文中的优化问题12.进行CVX求解权向量,部分仿真条件如表所示:
阵元数 | 10 |
信噪比 | 0 |
快拍数 | 100 |
扰噪比 | 20 |
阵元间距 | 半波长 |
信源数 | 1 |
干扰数 | 2 |
信源到达角 | 5° |
干扰到达角 | -30°,45° |
信号域宽度 | 4° |
角度误差上限 | 2° |
角度采样步长 | 1° |
期望信号误差控制量 | {6,4,2,4,6}*0.1 |
导向矢量误差分量上确界 | 0.24 |
干扰信号误差控制量 | , |
干扰抑制优先权重 | , |
噪声和干扰抑制的权衡参数 |
代码如下所示:
clc;
clear;
jk = 0;
snr = 0;
deg_dev_predf = 2;%%预定义的角度最大偏移量
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 = 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_deviate = 3;%%误差角度偏离
source_dev = source_aoa+[ang_mismatch1,ang_mismatch2,ang_mismatch3];
reco_size = 4;
reso = 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;
Us = U(:,1:source_num);
Un = U(:,source_num+1:end);
%% 算法
%%mu
para_mu = max(abs(diag(Gamma)))*10;
%%delta
para_delta = 0.24;
ds_ang_set = source_dev(1)-reco_size:reso:source_dev(1)+reco_size;
ds_sv = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(ds_ang_set));
%%c_theta
c_theta = (abs(-2*floor(length(ds_ang_set)/2):2:2*floor(length(ds_ang_set)/2))+2)*0.1;
intf_ang_set = [];
%%gamma_k and c_phi
for ik = 2:source_numintf_ang_set = [intf_ang_set,source_dev(ik)-reco_size:reso:source_dev(ik)+reco_size];
end
intf_sv = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(intf_ang_set));
c_phi = zeros(1,length(intf_ang_set));
for ik = 1:length(intf_ang_set)sv_current = intf_sv(:,ik);c_phi(ik) = sqrt(sv_current'*inv(R)*sv_current);
end
c_phi = c_phi/max(c_phi);
beta_k = zeros(1,source_num-1);
for ik = 2:source_numphi_k = source_dev(ik)-reco_size:reso:source_dev(ik)+reco_size;for jk = 1:length(phi_k)sv_current = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(phi_k(jk)));beta_k(ik-1) = beta_k(ik-1)+1/(sv_current'*inv(R)*sv_current);end
end
beta_k = beta_k/max(beta_k);
gamma_k = [];
for ik = 2:source_numgamma_k = [gamma_k,beta_k(ik-1)*ones(1,length(phi_k))];
end
c_theta_length = length(c_theta);
c_phi_length = length(c_phi);
c_phi_gamma_k = c_phi./gamma_k;
cvx_begin quietvariable w0(array_num,1) complexvariable t minimize(w0'*R*w0+para_mu*t);subject to%%constrainsts1for ik = 1:c_theta_lengthabs(w0'*ds_sv(:,ik)-1)+para_delta*norm(w0,2) <= c_theta(ik);end%%constrainsts2for ik = 1:c_phi_length(abs(c_phi_gamma_k(ik)*w0'*intf_sv(:,ik))+para_delta*norm(w0*sqrt(c_phi_gamma_k(ik)),2)) <= t;end
cvx_end
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)));
for ik = 1:length(outputval)if outputval(ik) <= -100outputval(ik) = -100;end
end
figure()
plot(theta,outputval,'LineWidth',2);
end
仿真结果如下所示:
由于本次仿真只采用了10个阵元,故从波束图上不能很好地看出惩罚项的添加对干扰域信号的抑制作用(实际上可从上图模糊看出干扰角度那块的旁瓣低了不少),该类算法还有一个痛点个人觉得便是直接通过凸优化求解权向量只能大致地对干扰域邻域内的信号进行抑制,而不能像导向矢量约束类算法那样对干扰进行精确抑制。该篇文献中对min-max问题的松弛及对优化问题15的ADMM求解是本人觉得的亮点,后者通过将难处理的约束项拉到增广拉格朗日式的最小化问题中当作约束项,一定程度上减轻了绝对值操作的处理难度。
代码可能有错,大佬们若发现问题可及时指出,谢谢。
参考文献
Wenqiang Pu, Jinjun Xiao, Tao Zhang, Zhi-Quan Luo, A penalized inequality-constrained approach for robust beamforming with DoF limitation, Signal Processing, Volume 202, 2023, 108746.
这篇关于波束形成 基于不等式惩罚项约束与自由度限制的稳健自适应波束形成算法的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!