稀疏DOA估计的经典算法——l1-SVD算法
On-Grid类DOA估计经典算法——
l
1
−
SVD
l_1-\text{SVD}
l1−SVD
文献”A Sparse Signal Reconstruction Perspective for Source Localization With Sensor Arrays”提出了一种稀疏表示的DOA定位算法，它属于On-Grid类算法的范畴。其核心要点有二：其一是，通过了奇异值(SVD)分解，把以大量快拍数衡量的信号模型，转换成以信源数衡量的低维信号模型；其二是，以二阶锥规划法替代通用的非线性优化方法来处理问题，使得算法更加高效。
目录
- On-Grid类DOA估计经典算法——
l
1
−
SVD
l_1-\text{SVD}
l1−SVD
- 过完备字典模型
l
1
l_1
l1范数最小化
- SVD降维
- 二阶锥(SOC)优化
- 二阶锥的定义与约束
- 二阶锥优化问题
- 理解问题的转化
- 与MUSIC谱估计对比
- 实验代码
- 文献
过完备字典模型
常用的参数标注：
M
M
M——均匀阵列的阵元数，
K
K
K——信源数，
T
T
T——时域快拍数，
N
θ
N_{\theta}
Nθ——过完备字典的尺寸
一般的，远场DOA估计模型为：
y
(
t
)
=
A
(
θ
)
s
(
t
)
+
n
(
t
)
y(t) = \boldsymbol{A}(\theta)s(t)+n(t)
y(t)=A(θ)s(t)+n(t)其中，
A
(
θ
)
∈
C
M
×
K
\boldsymbol{A}(\theta) \in \mathbb{C}^{M\times K}
A(θ)∈CM×K是阵列流形矩阵，
s
(
t
)
s(t)
s(t)和
n
(
t
)
n(t)
n(t)分别是
t
t
t时刻的信号向量和噪声向量，它们彼此独立。
同时，从稀疏的角度看待这个接收信号，我们可以把它重新写作：
Y
=
A
x
+
N
\boldsymbol{Y}=\boldsymbol{A}\boldsymbol{x}+\boldsymbol{N}
Y=Ax+N这时，
Y
∈
C
M
×
T
\boldsymbol{Y}\in \mathbb{C}^{M\times T}
Y∈CM×T仍为信号接收矩阵；而
A
∈
C
M
×
N
θ
\boldsymbol{A} \in \mathbb{C}^{M\times N_{\theta}}
A∈CM×Nθ为过完备集的流形矩阵，过完备集可以视作：由空域划分出来方向网格(Grid)组成的集合，如过完备集为：
S
θ
=
{
−
9
0
∘
:
1
∘
:
9
0
∘
}
\mathcal{S}_{\theta} = \{-90^{\circ}:1^{\circ}:90^{\circ}\}
Sθ={−90∘:1∘:90∘}其对应的集合尺寸为：
N
θ
=
181
N_{\theta} = 181
Nθ=181，对应的稀疏信号向量为：
x
∈
C
N
θ
×
1
\boldsymbol{x}\in \mathbb{C}^{N_{\theta}\times 1}
x∈CNθ×1，当然，它的绝大多数的元素都为0，所以才”稀疏“。
因此，恢复出信号向量
x
\boldsymbol{x}
x中不为0位置的索引，即代替了DOA估计的问题。
l
1
l_1
l1范数最小化
因为有了噪声项的干扰，信号
x
x
x出现了一些额外的非0项。因而，为了恢复这个信号向量，我们自然而然把最小化
l
0
l_0
l0范数作为目标函数，即使非0项数最小。
min
∥
x
∥
0
\min\|\boldsymbol{x}\|_0
min∥x∥0然而，
l
0
l_0
l0范数是非凸的，所以一般将其凸松弛为
l
1
l_1
l1范数来处理，即，向量元素绝对值的和。
min
∥
x
∥
1
\min \|\boldsymbol{x}\|_1
min∥x∥1另一方面，我们常用观测矩阵的”弗罗贝尼乌斯”范数(F)来衡量观测矩阵
Y
\boldsymbol{Y}
Y和理想接收矩阵
A
x
\boldsymbol{A}\boldsymbol{x}
Ax之间的差距，就是矩阵元素平方和再开根号：
∥
Y
−
A
x
∥
F
2
≤
β
2
\|\boldsymbol{Y} – \boldsymbol{A}\boldsymbol{x}\|_F^2 \leq \beta^2
∥Y−Ax∥F2≤β2这里的参数
β
2
\beta^2
β2指明了：我们能够允许的噪声参数的大小，一般可以从噪声功率处确定。
以上优化问题亦可以写成以下这种无约束的形式：
min
∥
Y
−
A
x
∥
F
2
+
λ
∥
x
∥
1
\min \|\boldsymbol{Y} – \boldsymbol{A}\boldsymbol{x}\|_F^2 + \lambda\|\boldsymbol{x}\|_1
min∥Y−Ax∥F2+λ∥x∥1这里的参数
λ
\lambda
λ用于控制”空间谱的稀疏性“和”范数偏移尺度“的平衡，约
0.1
−
10
0.1-10
0.1−10，它更像是一个经验值。
SVD降维
奇异值分解常常用于分析组成一个矩阵的主要成分，这很适合分析受噪声干扰的信号。将观测信号矩阵
Y
\boldsymbol{Y}
Y进行奇异值分解：
Y
=
U
L
V
H
\boldsymbol{Y} = \boldsymbol{U}\boldsymbol{L}\boldsymbol{V}^H
Y=ULVH其中，
U
∈
C
M
×
M
\boldsymbol{U}\in \mathbb{C}^{M\times M}
U∈CM×M和
V
∈
C
T
×
T
\boldsymbol{V}\in \mathbb{C}^{T\times T}
V∈CT×T都是酉矩阵，矩阵
L
∈
C
M
×
T
\boldsymbol{L}\in \mathbb{C}^{M\times T}
L∈CM×T是由一个
M
×
M
M\times M
M×M的对角矩阵和
M
×
(
T
−
M
)
M\times(T-M)
M×(T−M)的0矩阵构成，而对角矩阵中有
K
K
K个大元素和
(
M
−
K
)
(M-K)
(M−K)个小元素，它们之间存在数量级的差距。因此，很自然地，我们可以把这些主成分过滤出来。
定义一个选择矩阵
D
K
∈
C
T
×
K
\boldsymbol{D}_K\in \mathbb{C}^{T\times K}
DK∈CT×K，它由一个尺寸为
K
K
K的单位矩阵和
(
T
−
K
)
×
K
(T-K)\times K
(T−K)×K的0矩阵构成。
Y
S
V
=
Y
V
D
K
=
U
L
D
K
\boldsymbol{Y}_{SV} = \boldsymbol{Y}\boldsymbol{V}\boldsymbol{D}_K = \boldsymbol{U}\boldsymbol{L}\boldsymbol{D}_K
YSV=YVDK=ULDK得到一个新的接收矩阵
Y
S
V
∈
C
M
×
K
\boldsymbol{Y}_{SV}\in \mathbb{C}^{M\times K}
YSV∈CM×K，显然，它降低了信号的尺寸，使得算法更加高效，更适合高快拍的情景。
那么，此时，经过线性变换
V
D
K
\boldsymbol{V}\boldsymbol{D}_K
VDK，信号的模型变为：
Y
S
V
=
A
S
S
V
+
N
S
V
\boldsymbol{Y}_{SV} = \boldsymbol{A}\boldsymbol{S}_{SV} + \boldsymbol{N}_{SV}
YSV=ASSV+NSV对应的，
l
1
l_1
l1优化问题转变成了：
min
S
S
V
∥
Y
S
V
−
A
S
S
V
∥
F
2
+
λ
∥
s
~
l
2
∥
1
\min_{\boldsymbol{S}_{SV}} \|\boldsymbol{Y}_{SV} – \boldsymbol{A}\boldsymbol{S}_{SV}\|_F^2 + \lambda\|\boldsymbol{\tilde{s}}^{l_2}\|_1
SSVmin∥YSV−ASSV∥F2+λ∥s~l2∥1由于
S
S
V
\boldsymbol{S}_{SV}
SSV相比于
x
\boldsymbol{x}
x变成了矩阵，因而，
s
~
l
2
\boldsymbol{\tilde{s}}^{l_2}
s~l2为矩阵
S
S
V
\boldsymbol{S}_{SV}
SSV每一行向量的
l
2
l_2
l2范数，用以替代
x
\boldsymbol{x}
x的每一个元素。
二阶锥(SOC)优化
以上即为
l
1
−
SVD
l_1-\text{SVD}
l1−SVD算法的核心思想，以上优化问题是凸的，通过求解非线性优化即可。
而作者进一步将该问题转化为了二阶锥规划(SOCP)的框架，进而使用更高效的内点法求解。
二阶锥的定义与约束
定义二阶锥：
C
k
=
{
[
u
t
]
:
u
∈
R
K
−
1
,
t
∈
R
,
∥
u
∥
≤
t
}
\mathcal{C}_k = \left\{\begin{bmatrix} \boldsymbol{u} \\ t \end{bmatrix}: \boldsymbol{u} \in \mathbb{R}^{K-1}, t \in \mathbb{R},\|\boldsymbol{u}\| \leq t \right\}
Ck={[ut]:u∈RK−1,t∈R,∥u∥≤t}二阶锥约束可以表示为：
∥
A
x
+
b
∥
≤
c
T
x
+
d
\|\boldsymbol{A}\boldsymbol{x}+\boldsymbol{b}\| \leq \boldsymbol{c}^T\boldsymbol{x}+\boldsymbol{d}
∥Ax+b∥≤cTx+d等价于：
[
A
c
T
]
x
+
[
b
d
]
∈
C
k
\begin{bmatrix} \boldsymbol{A} \\ \boldsymbol{c}^T \end{bmatrix} \boldsymbol{x} + \begin{bmatrix} \boldsymbol{b} \\ \boldsymbol{d} \end{bmatrix} \in \mathcal{C}_k
[AcT]x+[bd]∈Ck
二阶锥优化问题
一般形式：
min
x
f
T
x
s
.
t
∥
A
i
x
+
b
i
∥
≤
c
i
T
x
+
d
i
F
x
=
g
\min_{\boldsymbol{x}} \quad \boldsymbol{f}^T\boldsymbol{x} \\ s.t \quad \|\boldsymbol{A}_i\boldsymbol{x}+\boldsymbol{b}_i\| \leq \boldsymbol{c}_i^T\boldsymbol{x}+\boldsymbol{d}_i \\ \quad \boldsymbol{F}\boldsymbol{x} = \boldsymbol{g}
xminfTxs.t∥Aix+bi∥≤ciTx+diFx=g其中，
f
∈
R
n
\boldsymbol{f} \in \mathbb{R}^n
f∈Rn，
A
i
∈
R
n
i
×
n
\boldsymbol{A}_i\in \mathbb{R}^{n_i \times n}
Ai∈Rni×n，
b
i
∈
R
n
i
\boldsymbol{b}_i\in \mathbb{R}^{n_i}
bi∈Rni，
c
i
∈
R
n
\boldsymbol{c}_i\in \mathbb{R}^{n}
ci∈Rn，
d
i
∈
R
\boldsymbol{d}_i\in \mathbb{R}
di∈R，
F
∈
R
p
×
n
i
\boldsymbol{F} \in \mathbb{R}^{p \times n_i}
F∈Rp×ni和
g
∈
R
p
\boldsymbol{g} \in \mathbb{R}^{p}
g∈Rp。另外，待优化变量
x
∈
R
n
\boldsymbol{x} \in \mathbb{R}^{n}
x∈Rn。
注意，二阶锥优化问题对目标函数没有绝对的要求。
理解问题的转化
作者Malioutov最终把非线性优化问题转化成了SOCP问题：
min
p
,
q
,
r
,
S
S
V
p
+
λ
q
s
.
t
∥
vec
{
Y
S
V
−
A
S
S
V
}
∥
2
2
≤
p
1
T
r
≤
q
∥
S
S
V
(
i
,
:
)
∥
2
≤
r
i
\min_{p,q,\boldsymbol{r},\boldsymbol{S}_{SV}} \quad p + \lambda q \\ s.t \quad \|\text{vec}\{\boldsymbol{Y}_{SV} – \boldsymbol{A}\boldsymbol{S}_{SV}\}\|_2^2 \leq p \\ \boldsymbol{1}^T\boldsymbol{r} \leq q \\ \|\boldsymbol{S}_{SV}(i,:)\|_2 \leq r_i
p,q,r,SSVminp+λqs.t∥vec{YSV−ASSV}∥22≤p1Tr≤q∥SSV(i,:)∥2≤ri二阶锥规划的约束条件，简单理解为：待优化变量的仿射变换的范数小于等于另一种仿射变换。
因此，对应关系为：
x
=
[
p
q
r
vec
{
S
S
V
}
]
\boldsymbol{x} = \begin{bmatrix} p \\ q \\ \boldsymbol{r} \\ \text{vec}\{\boldsymbol{S}_{SV}\} \end{bmatrix}
x=⎣
⎡pqrvec{SSV}⎦
⎤第一个约束条件转换为：
∥
vec
{
Y
S
V
}
−
(
I
⊗
A
)
vec
{
S
S
V
}
∥
≤
p
\|\text{vec}\{\boldsymbol{Y}_{SV}\} – \left(\boldsymbol{I}\otimes\boldsymbol{A}\right)\text{vec}\{\boldsymbol{S}_{SV}\}\| \leq p
∥vec{YSV}−(I⊗A)vec{SSV}∥≤p显然，这符合了二阶锥规划的约束条件的定义。
与MUSIC谱估计对比
两个信源的对比：
两个相近的信源的谱对比：
实验代码
clc;clear all;close all;
twpi = 2*pi;
derad = pi/180;
j = sqrt(-1);
%% 生成信号
M = 8;%阵元数
c = 1500;%声速1500m/s
f = 1e4;%载频10kHz
lambda0 = c/f;%载波波长
dd = 0.5*lambda0;%阵元间距
d = 0:dd:(M-1)*dd;
theta = [10 15];
K = length(theta);%信源数
fs = 2.5*f;%采样频率
L = 100;%快拍数
t = [1:L]/fs;%时域采样
A = exp(j*twpi*d'*sin(theta*derad)/lambda0);%方向向量
S = randn(K,L).*exp(j*twpi*ones(K,1)*f*t);
snr = 20;%信噪比
Y = awgn(A*S,snr,'measured');%加入高斯白噪声
%% 构造完备字典
Grid = -90:1:90;
AA = exp(j*twpi*d'*sin(Grid*derad)/lambda0);
Dk1 = eye(K);
Dk2 = zeros(K,L-K);
Dk = [Dk1,Dk2].';
[U,~,V] = svd(Y); %进行奇异值分解
Ysv = Y*V'*Dk;
%% 求解凸优化问题
cvx_begin quiet
variables p q
variables r(length(Grid))
variable SSV1(length(Grid),K) complex;
expression xsv(length(Grid),1)
expressions Zk(M,K) complex
minimize( p + 2*q );
subject to
Zk = Ysv - AA*SSV1;
Zvec = vec(Zk); %把矩阵转化为向量
norm(Zvec,2) <= p; %第一个不等式约束
sum(r) <= q; %第二个不等式约束
for i=1:length(Grid) %第三个不等式约束
norm(SSV1(i,:),2) <= r(i);
end
cvx_end
%% 绘制空间谱
power = 10*log10(abs(SSV1(:,1))/max(abs(SSV1(:,1))));
h1 = plot(Grid,power,'r');
hold on
xlabel('DOA/degree');
ylabel('PowerdB');
%% 绘制MUSIC空间谱
R = Y*Y'/L;
[EV,D] = eig(R);%特征值分解
DD = diag(D);%将特征值变为向量形式
[DD,I] = sort(DD);%从小到大
EV = fliplr(EV(:,I));
%%%%%%%%%%构造MUSIC的谱函数%%%%%%%%%%%
En = EV(:,K+1:M);%噪声子空间
Pmusic = zeros(1,length(Grid));
for ii = 1:length(Grid)
Pmusic(ii) = 1/(AA(:,ii)'*En*En'*AA(:,ii));
end
Pmusic=abs(Pmusic);
Pmax = max(Pmusic);
Pmusic_db = 10*log10(Pmusic/Pmax);
h2 = plot(Grid,Pmusic_db,'b');
xlim([-90 90])
set(gca,'XTick',[-90:30:90])
grid on
hold on
m = ylim;
for i = 1:K
line([theta(i) theta(i)],[m(1) 0],'Color','k','LineWidth',1,'LineStyle','--');
hold on
end
legend('l1-SVD','MUSIC','Direction of Arrival');
Tick',[-90:30:90])
grid on
hold on
m = ylim;
for i = 1:K
line([theta(i) theta(i)],[m(1) 0],'Color','k','LineWidth',1,'LineStyle','--');
hold on
end
legend('l1-SVD','MUSIC','Direction of Arrival');
文献
链接：文献获取
提取码：6666
本文来自网络，不代表协通编程立场，如若转载，请注明出处：https://net2asp.com/7dffe235a7.html
相关推荐
-
ArcNeural： AI 时代的多模数据库丨技术专栏
导读 本文根据 Fabarta 资深技术专家谭宇在“2023 中国软件技术大会”演讲实录整理而来。围绕以下四个方面进行介绍：首先简单介绍Fabarta 背景以及我们为什么要研发 A…
-
AI 神助攻，协同办公神器 —- ONLYOFFICE
人工智能不会取代人，只会淘汰那些不会使用人工智能的人。 – 鲁迅 一、人工智能重新定义办公新模式 随着GPT的横空出世，AI的应用场景已经无处不在，从智能客服、智能语音助手、智能家…
-
【人工智能课程】计算机科学博士作业二
使用TensorFlow1.x版本来实现手势识别任务中，并用图像增强的方式改进，基准训练准确率0.92，测试准确率0.77，改进后，训练准确率0.97，测试准确率0.88。 1 导…
-
小说推文的又一春天，AI漫画推文，高原创玩法详细教程，掌握流量密码【附全套工具软件】
小说推文项目很早就有了，传统的小说推文视频大都以混剪解压视频、游戏视频为主，虽然这类视频制作起来简单，但是原创度低。 随着短视频平台的规则变化，平台对非原创作品的打压力度越来越大，…
-
使用 Gradio 构建生成式 AI 应用程序(一): 图片内容读取app
今天我们来学习DeepLearning.AI的在线课程：Building Generative AI Applications with Gradio，该课程主要讲述利用gradi…
-
ISP-长短曝光融合生成HDR图像
1、高动态范围图像相关 图像的动态范围是指一幅图像中量化的最大亮度与最小噪声的比值。高动态范围HDR（high dynamic range）图像,能够完整表示真实场景中跨度很大的动…
-
20分钟，使用Amazon SageMaker快速搭建属于自己的AIGC应用
真火！ 作为最近一段时间人工智能领域内的顶流之一，AIGC（AI-Generated Content）早已火爆出圈，频登各大互联网平台热搜。 cite: 微软亚洲研究院官方微博 这…
-
【AI前沿】chatgpt还有哪些不足？
博客昵称：吴NDIR 个人座右铭：得之淡然，失之坦然 作者简介：喜欢轻音乐、象棋，爱好算法、刷题 其他推荐内容：计算机导论速记思维导图 其他内容推荐：五种排序算法 在这个愉快的周末…
-
三维点云处理之最近邻问题
前言：本系列文章是关于三维点云处理的常用算法，深入剖析pcl库中相关算法的实现原理，并以非调库的方式实现相应的demo。 1. 最近邻问题概述 （1）最邻近问题：对于点云中的点，怎…
-
第一章：AI大模型概述1.1 AI大模型的定义与特点1.1.1 什么是AI大模型
1.背景介绍 AI大模型是指具有极大规模、高度复杂结构和强大计算能力的人工智能模型。这类模型通常用于处理大规模、高维度的数据，并能够实现复杂的智能任务，如自然语言处理、图像识别、推…
