文章目录

原子范数最小化(Atomic Norm Minimization)
- l0l_0l0-原子范数
- 最小化Toeplitz矩阵的秩
- - 引理1
  - 引理2
  - 问题的转化
- 转换为凸优化问题
- - 凸松弛(Relaxation)
  - l1l_1l1-原子范数最小化
  - 核范数最小化
  - 变换为矩阵的迹
- 稀疏阵列DOA估计
- 仿真结果
- 自己写的代码
- 相关文献与博客

原子范数最小化(Atomic Norm Minimization)

l0l_0l0-原子范数

原子集
A={a(θ,ϕ)=a(θ)ϕ:θ∈[−90∘,90∘],ϕ∈C,∣ϕ∣=1}\mathcal{A} = \{\boldsymbol{a}(\theta,\phi) = \boldsymbol{a}(\theta)\phi : \theta\in [-90^\circ,90^\circ], \phi \in \mathbb{C}, |\phi| = 1 \} A={a(θ,ϕ)=a(θ)ϕ:θ∈[−90∘,90∘],ϕ∈C,∣ϕ∣=1}
对于无噪声条件下的单快拍接收信号。
z=A(θ)sz = \boldsymbol{A}(\theta) \boldsymbol{s} z=A(θ)s
l0l_0l0范数的定义
∥z∥A,0=inf⁡sk,θk{K:z=∑k=1Ka(θk)sk,θk∈[−90∘,90∘]}\|\mathit{z}\|_{\mathcal{A},0} = \inf_{s_k,\theta_k}\{\mathcal{K}:z = \sum_{k = 1}^{\mathcal{K}} \boldsymbol{a} \left(\theta_k\right)s_k, \theta_k \in [-90^\circ,90^\circ] \} ∥z∥A,0=sk,θkinf{K:z=k=1∑Ka(θk)sk,θk∈[−90∘,90∘]}
l0l_0l0范数指的是，组成向量zzz的最小原子的个数。DOA估计反映在该问题中就是找出原子个数最小的情况下的方向向量的组合。

而对于稀疏阵列而言，信号zzz往往是残缺的，不完整的。信号恢复的过程，可以理解为：填充向量zzz，使其在原子集A\mathcal{A}A下，原子范数最小，即：
min⁡z∥z∥A,0\min_{z}\|\mathit{z}\|_{\mathcal{A},0} zmin∥z∥A,0

最小化Toeplitz矩阵的秩

引理1

对于任意的半正定矩阵T(z)∈CL×L\mathcal{T}\left(z\right) \in \mathbb{C}^{L \times L}T(z)∈CL×L，如果r=rank(T(z))≤L−1r=rank(\mathcal{T}(z)) \leq L-1r=rank(T(z))≤L−1，则Hermitian Toeplitz矩阵能够被唯一范德蒙分解为：
T(z)=∑k=1rpka(θk)aH(θk)=A(θ)PAH(θ)\mathcal{T}\left(z\right) = \sum_{k=1}^{r}p_k\boldsymbol{a}(\theta_k)\boldsymbol{a}^H(\theta_k) = \boldsymbol{A}(\theta)\boldsymbol{P}\boldsymbol{A}^H(\theta) T(z)=k=1∑rpka(θk)aH(θk)=A(θ)PAH(θ)
其中，A∈CL×r\boldsymbol{A}\in \mathbb{C}^{L\times r}A∈CL×r是范德蒙矩阵，它的列视为来自不同方向的信号的方向向量，P∈Cr×r\boldsymbol{P}\in \mathbb{C}^{r\times r}P∈Cr×r是对角矩阵。

引理2

考虑以下的式子：
[xzHzT(z)]⪰0\begin{bmatrix} x & \boldsymbol{z}^H \\ \boldsymbol{z} & \mathcal{T}\left(z\right) \end{bmatrix} \succeq 0 [xzzHT(z)]⪰0
其中，xxx是一个待优化的变量。这个式子对应了以下的结论：

T(z)⪰0\mathcal{T}\left(z\right)\succeq 0T(z)⪰0；
z\boldsymbol{z}z 一定位于Toeplitz矩阵 T(z)\mathcal{T}\left(z\right)T(z) 的列向量中。

问题的转化

由于上面的结论2，所以z\boldsymbol{z}z一定可以被A(θ)\boldsymbol{A}(\theta)A(θ)中的rrr个子向量线性表出。而这rrr个子向量就是矩阵T(z)\mathcal{T}\left(z\right)T(z)范德蒙分解后的rrr个特征向量。因此，向量z\boldsymbol{z}z的l0l_0l0范数最小化就是T(z)\mathcal{T}\left(z\right)T(z)矩阵秩的最小化。问题即可转换为：
min⁡x,zrank(T(z))s.t[xzHzT(z)]⪰0\min_{x,z} rank(\mathcal{T}\left(z\right)) \\ s.t \quad \begin{bmatrix} x & \boldsymbol{z}^H \\ \boldsymbol{z} & \mathcal{T}\left(z\right) \end{bmatrix} \succeq 0 x,zminrank(T(z))s.t[xzzHT(z)]⪰0
然而，该优化问题并不是一个凸优化，无法直接求解。

转换为凸优化问题

凸松弛(Relaxation)

解释：通过对问题限制条件的松弛，将原非凸问题转换为凸优化问题。注：松弛原问题，只能得到一个可行域更大的问题，如果原问题是求最小，则松弛后的问题的最优值一定小于等于原问题的最优值，这也是一种给出下界的方法。松弛不仅仅用于整数约束，只要利于将定义域非凸变为凸集即可。

例如，l0l_0l0范数（非凸）的最紧凸松弛为l1l_1l1范数（凸），即：将向量中非零元素的个数转为元素绝对值的和。

那么，类似的，可以将l0l_0l0-原子范数松弛为l1l_1l1-原子范数。或者，将rank(T(z))rank(\mathcal{T}\left(z\right))rank(T(z))松弛为核范数∥T(z)∥∗\|\mathcal{T}\left(z\right)\|_*∥T(z)∥∗。

l1l_1l1-原子范数最小化

l1l_1l1-原子范数的定义如下：
∥z∥A,1=inf⁡sk,θ{∑k∣sk∣:z=∑k=1Ka(θk)sk,θk∈T}\|z\|_{\mathcal{A},1} = \inf_{s_k,\theta}\{\sum_{k}{|s_k|}:z=\sum_{k=1}^{\mathcal{K}}\boldsymbol{a}(\theta_k)s_k,\theta_k\in\mathbb{T}\} ∥z∥A,1=sk,θinf{k∑∣sk∣:z=k=1∑Ka(θk)sk,θk∈T}
从形式上看，问题已经转换成：恢复一个单快拍信号向量zzz，使得其l1l_1l1-原子范数最小，即：
min⁡z∥z∥A,1\min_{z}\|z\|_{\mathcal{A},1} zmin∥z∥A,1
此时，最小化l1l_1l1范数依然看起来十分困难，所以再次利用范德蒙分解，将该问题转化为：
min⁡x,z112(x+z1)s.t[xzHzT(z)]⪰0\min_{x,z_1} \frac{1}{2}(x+z_1) \\ s.t \quad \begin{bmatrix} x & \boldsymbol{z}^H \\ \boldsymbol{z} & \mathcal{T}\left(z\right) \end{bmatrix} \succeq 0 x,z1min21(x+z1)s.t[xzzHT(z)]⪰0
其中，z1z_1z1表示单快拍向量zzz中的第一个元素。至此，该问题可以用CVX工具箱来求解了。

核范数最小化

核范数（矩阵奇异值的和）的定义如下：
∥T(z)∥∗=Tr(TH(z)T(z))\|\mathcal{T}(z)\|_* = Tr\left(\sqrt{\mathcal{T}^H(z)\mathcal{T}(z)}\right) ∥T(z)∥∗=Tr(TH(z)T(z))
核范数∥T(z)∥∗\|\mathcal{T}(z)\|_*∥T(z)∥∗是一个凸函数，证明如下：

由于∥T(z)∥∗\|\mathcal{T}(z)\|_*∥T(z)∥∗和∥T(z)∥2\|\mathcal{T}(z)\|_2∥T(z)∥2为一对对偶范数，因此，若∥T(z)∥2\|\mathcal{T}(z)\|_2∥T(z)∥2为凸函数，那么∥T(z)∥∗\|\mathcal{T}(z)\|_*∥T(z)∥∗即为凸函数。
∥T∥2=sup⁡u,v{uHTv:∥u∥2=1,∥v∥2=1}\|\mathcal{T}\|_2 = \sup_{\boldsymbol{u},\boldsymbol{v}}\{\boldsymbol{u}^H\mathcal{T}\boldsymbol{v} : \|\mathcal{u}\|_2=1,\|\mathcal{v}\|_2=1\} ∥T∥2=u,vsup{uHTv:∥u∥2=1,∥v∥2=1}
向量u\boldsymbol{u}u和v\boldsymbol{v}v为各个奇异值对应的向量，uHTv\boldsymbol{u}^H\mathcal{T}\boldsymbol{v}uHTv是对T\mathcal{T}T的仿射变换，因此凸。由下面”逐点上确界“的性质，所以∥T∥2\|\mathcal{T}\|_2∥T∥2凸。
f(x)=sup⁡y,z{g(x,y,z):y∈Ay,z∈Az}f(x) = \sup_{y,z}\{g(x,y,z):y\in\mathcal{A}_y,z\in\mathcal{A}_z\} f(x)=y,zsup{g(x,y,z):y∈Ay,z∈Az}
ggg为凸函数，那么fff亦为凸函数，因而∥T(z)∥∗\|\mathcal{T}(z)\|_*∥T(z)∥∗的凸性得证。

而取矩阵秩的操作的最紧凸松弛为核范数，所以，可以将之前的非凸问题转换成：
min⁡x,z∥T(z)∥∗s.t[xzHzT(z)]⪰0\min_{x,z}\|\mathcal{T}(z)\|_* \\ s.t \quad \begin{bmatrix} x & \boldsymbol{z}^H \\ \boldsymbol{z} & \mathcal{T}\left(z\right) \end{bmatrix} \succeq 0 x,zmin∥T(z)∥∗s.t[xzzHT(z)]⪰0
该问题也可以通过CVX工具箱来求解。

变换为矩阵的迹

由于半正定的Hermitian Toeplitz矩阵∥T(z)∥\|\mathcal{T}(z)\|∥T(z)∥，所以其存在唯一的范德蒙分解：
T(z)=∑k=1rpka(θk)aH(θk)\mathcal{T}\left(z\right) = \sum_{k=1}^{r}p_k\boldsymbol{a}(\theta_k)\boldsymbol{a}^H(\theta_k) T(z)=k=1∑rpka(θk)aH(θk)
取该矩阵的迹：
Tr(T(z))=L∑k=1rpk∥z∥A,1=∑k=1rpk=Tr(T(z))/LTr(\mathcal{T}\left(z\right)) = L\sum_{k=1}^{r}p_k \\ \|z\|_{\mathcal{A},1} = \sum_{k=1}^{r}p_k = Tr(\mathcal{T}\left(z\right)) / L Tr(T(z))=Lk=1∑rpk∥z∥A,1=k=1∑rpk=Tr(T(z))/L
上式中，第一个等号可以理解为：由于方向向量a(θk),k=1,2,…\boldsymbol{a}(\theta_k),k=1,2,\ldotsa(θk),k=1,2,… 彼此相互独立，所以只要确定了信号向量z\boldsymbol{z}z存在由a(θk)\boldsymbol{a}(\theta_k)a(θk)张成的向量空间中，那么这一组系数pk,k=1,2,…,rp_k,k=1,2,\dots,rpk,k=1,2,…,r便可唯一确定。因此，取下界操作inf⁡\infinf就可无视。

因此，该凸问题就转换为：
min⁡zTr(T(z))s.tT(z)⪰0;\min_{z} \quad Tr(\mathcal{T}\left(z\right)) \\ s.t \quad \mathcal{T}\left(z\right) \succeq 0; zminTr(T(z))s.tT(z)⪰0;
该问题也可以通过CVX工具箱来求解。

稀疏阵列DOA估计

step1 通过图(a)中实际布置的天线阵列接收数据X\boldsymbol{X}X，计算接收阵列的协方差矩阵Rx\boldsymbol{R}_xRx;

step2 向量化Rx\boldsymbol{R}_xRx，并排序去冗余，得到图(b)中等效的单快拍信号向量zex\boldsymbol{z}_{ex}zex，包括的位置-12~12，空洞处用0元素代替;

step3 用zex\boldsymbol{z}_{ex}zex构造Toeplitz矩阵Rv\boldsymbol{R}_vRv，另外，定义与Rv\boldsymbol{R}_vRv同等大小的二进制矩阵G\boldsymbol{G}G，用来反映Rv\boldsymbol{R}_vRv有值的位置，z\boldsymbol{z}z代表待恢复的单快拍信号(对应的位置0~12);

step4 用CVX，求解以下任一优化问题，获取重构的Toeplitz矩阵T(z)\mathcal{T}\left(z\right)T(z);

仿真结果

信噪比20dB20dB20dB，快拍数200200200，来波方向−50∘-50^\circ−50∘至50∘50^\circ50∘，间隔5∘5^\circ5∘，共21个信源。

仿真使用的是增广展开的互质阵列（Augmented-Expanding Co-Prime Array），共10个阵元（M=3,N=5M=3,N=5M=3,N=5）。

原子范数最小化 ANM

核范数最小化 NNM

最小化迹的方法

相比之下，ANM算法的效果更好一些。

自己写的代码

如有不足，还望各位在评论区不吝赐教。

clc;clear all;close all;
%%%%%增广展开互质阵列（Augmented-Expanding Co-Prime Array）：填充内插阵元算法%%%%%
twpi = 2*pi;
derad = pi/180;
j = sqrt(-1);
c = 1500;%声速
f0 = 1e4;%载频
lambda = c/f0;%波长
dd = lambda/2;%阵元间距
M = 3;%子阵1的阵元数
N = 5;%子阵2的阵元数
kelm = 2*M+N-1;%实际阵元总数
d = [0:-N:(1-2*M)*N,M:M:(N-1)*M];%展开扩展互质阵
D = zeros(1,kelm*kelm);
for m=1:kelmD(1,(m-1)*kelm+1:m*kelm) = d-d(m);
end
[loc,index] = sort(D);
%%%%%构造排序去冗余矩阵FD%%%%%
indD = unique(D);%去除重复项
dof = length(indD);%阵列自由度
FD = zeros(dof,length(D));
for ii=1:length(D)count = length(find(D(ii)==D));temp = find(D(ii)==indD);FD(temp,ii) = 1/count;
end
%%%%%接收端%%%%%
d = d*dd;
% theta = [-0.1 0.1];%高精度到达角方向
theta = -45:6:45;%多到达角
iwave = length(theta);%信源数
snr = 30;%信噪比
n = 500;%采样数
fs = 2*f0;%采样频率
t = [1:n]/fs;
Amp = randn(iwave,n);
S = Amp.*(ones(iwave,1)*exp(j*twpi*(f0*t)));%独立信源
Ar = exp(j*twpi*d.'*sin(theta*derad)/lambda);
X = awgn(Ar*S,snr,'measured');
Rx = X*X'/n;
z = reshape(Rx,kelm*kelm,1);%转换为向量形式，成单快拍形式
%%%%%构造虚拟单快拍接收矢量：排序去冗余%%%%%
dm = unique(D);
zD = FD*z;
LD = (dof+1)/2;%%%%%构造RD矩阵%%%%%
RD = zeros(LD,LD);%重构的协方差矩阵
for mr = 1:LDtemp = LD+1-mr;RD(:,mr) = zD(temp:temp+LD-1);
end
vec_RD = reshape(RD,LD*LD,1);
%%%%%增加阵列自由度%%%%%
dm_ex = min(dm):1:max(dm);
zV = zeros(max(dm)-min(dm)+1,1);
dmax = max(dm);
for ii=1:dofzV(dm(ii)+dmax+1,1) = zD(ii,1);
end
%%%%%构造RV矩阵%%%%%%
LV = (length(zV)+1)/2;%Teoplitz矩阵大小
RV = zeros(LV,LV);%重构的协方差矩阵
for mr = 1:LVtemp = LV+1-mr;RV(:,mr) = zV(temp:temp+LV-1);
end%%%%%构造选择矩阵G%%%%%
G = ones(LV,LV);
posi0V = find(RV==0);
G(posi0V) = 0;%%%%%利用CVX工具箱求解凸优化问题%%%%%
mu = 0.25;
cvx_begin sdp quiet
% cvx_precision high
cvx_solver sdpt3variable T(LV,LV) hermitian toeplitz semidefinite;minimize( square_pos(norm(T .* G - RV,'fro')) + mu * trace(T) )    %目标函数
cvx_end%%%%%%%%%%构造MUSIC的谱函数%%%%%%%%%%%
[EV,Dv] = eig(T);%特征值分解
DD = diag(Dv);%将特征值变为向量形式
[DD,I] = sort(DD);%从小到大
DD = fliplr(DD');%翻转函数，从大到小
EV = fliplr(EV(:,I));
En = EV(:,iwave+1:end);%噪声子空间
dm_ss = dm_ex(1:LV)*dd;
for ii = 1:2001angle(ii) = (ii-1001)*90/1000;phim = derad*angle(ii);a = exp(j*twpi*dm_ss*sin(phim)/lambda).';Pmusic(ii) = 1/(a'*En*En'*a);
end
Pmusic = abs(Pmusic);
Pmax = max(Pmusic);
Pmusic_db = 10*log10(Pmusic/Pmax);
% [maxv_music,maxl_music]=findpeaks(Pmusic_db,'minpeakdistance',0.001,'minpeakheight',mean(Pmusic_db));
% [vl,ind] = sort(maxv_music,'descend');
% doa = angle(maxl_music(ind(1:iwave)));
h = plot(angle,Pmusic_db);
set(h,'Linewidth',0.2)
xlabel('angle(degree)')
ylabel('magnitude(dB)')
xlim([-90 90])
set(gca,'XTick',[-90:10:90])
grid on
hold on
m = ylim;
for i = 1:iwaveline([theta(i) theta(i)],[m(1) 0],'Color','k','LineWidth',0.1,'LineStyle','--');hold on
end
legend('ANM','Direction of Arrival');