波束形成基于原子范数最小化的稳健自适应波束形成

量大龙神魔

已于 2023-08-13 16:04:34 修改

阅读量2.1k

点赞数 14

文章标签：算法

于 2023-05-23 19:25:12 首次发布

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

版权

Compress Sensing 专栏收录该内容

2 篇文章 4 订阅

订阅专栏

前言

原子范数最小化作为压缩感知中比较热门的算法，因其通过在集合内寻找最优原子组合的思想，相较于其他On-grid、Off-grid的算法，Gridless的原子范数最小化算法具有最高的精度，但是随之而来的便是较高的复杂度与计算度，此外，在DOA及波束形成这一块由于其导向矢量结构特性与范德蒙德分解具有着极高的适配性，因此采用原子范数最小化这一方法具有更高的精度，这篇文章是20年由国防科大电科院的老师于TSP上进行发表。关于原子范数最小化及其在DOA中的应用场景，可以参考压缩感知的尽头: 原子范数最小化，该博客作者在其中已经提供了详细的公式推导。

模型提出及具体实现

问题提出

一般来讲，现有的大部分提出的稳健自适应波束形成算法仍然是基于MVDR，通过求解维纳最优解得到当前协方差矩阵下的最优权向量。该问题可以表示为如下形式：

$\min\limits_{\mathbf{w}}\mathbf{w}^H\mathbf{R}\mathbf{w},\text{s.t.}\mathbf{w}^H\mathbf{a}(\theta)=1$

在早期协方差矩阵通常以估计得到的数据协方差阵进行直接带入，这导致了在高信噪比下的信号自消的现象，近年来逐渐以干扰加噪声协方差阵代替，一定程度上提高了性能，在这篇文章中，作者通过将干扰矩阵带入，对协方差阵进行求解。即该问题可被表示为：

$\min\limits_{\mathbf{w}}\mathbf{w}^H\mathbf{R}_J\mathbf{w},\text{s.t.}\mathbf{w}^H\mathbf{a}(\theta)=1$

其中为干扰协方差阵，一说为干扰子空间，其理想情况下具有的特性为：半正定性，是一个典型的托普利茨结构的矩阵，并且可被范德蒙德分解。因此作者的核心思想便是通过原子范数最小化，得到以名义导向矢量集表示的与接收信号协方差矩阵近似的重构干扰协方差矩阵。

这篇文章中作者以表示阵元数，以表示快拍数，以表示干扰数。

原子范数最小化

定义原子集的表示形式为：

$\mathcal{A}=\{\textbf{A}(\theta,\textbf{b})|\theta\in[-\frac{\pi}{2},\frac{\pi}{2}],\|\textbf{b}\|_2=1\}$

其中的原子为：

$\textbf{A}(\theta,\textbf{b})=\textbf{a}(\theta)\textbf{b}^T$

其中 $\textbf{b}\in\mathcal{C}^{M\times1}$ ， $\|\mathbf{b}\|_2=1$ ，为阵元数，则接收信号中干扰信号分量可表示为：

$\textbf{J}=\sum_{i=1}^K c_i\textbf{A}(\theta_i,\textbf{b}_i)$

则可通过在接收信号的F范数邻域内，寻找能够被最少的原子数表示出来的干扰信号分量，即：

$\min\limits_{\mathbf{X},\theta_0,\mathbf{s}}\left\|\mathbf{X}-\mathbf{a}(\theta_0)\mathbf{s}^T\right\|_{\mathcal{A},0},s.t.\frac{1}{2}\left\|\mathbf{Y}-\mathbf{X}\right\|_{\mathcal{F}}^2\le\eta$

其中 $\|\cdot\|_{\mathcal{A},0}$ 便是原子范数，可表示为：

$\|\mathbf{X}\|_{\mathcal{A},0}=\inf\limits_r\bigg\{\mathbf{X}=\sum_i^rc_i\mathbf{A}(\theta_i,\mathbf{b}_i),c_i\ge0\bigg\}$

该原子范数显然是个NP-hard问题，但是原子范数有着一条非常优秀的性质，就是可以像传统的范数松弛成一个范数(这部分证明可参考B站-中科大凸优化-凸优化19)一样，从原子范数松弛成一个凸原子范数：

$\min\limits_{\mathbf{X},\theta_0,\mathbf{s}}\left\|\mathbf{X}-\mathbf{a}(\theta_0)\mathbf{s}^T\right\|_{\mathcal{A}}\quad\text{s.t.}\frac{1}{2}\left\|\mathbf{Y}-\mathbf{X}\right\|_{\mathcal{F}}^2\leq\eta$

其中 $||\cdot||_{\mathcal{A}}$ 为原子范数，可表示为：

$\|\mathbf{X}\|_\mathcal{A}=~\inf\left\{\sum_i c_i|\mathbf{X}=\sum_i c_i\mathbf{A}(\theta_i,\mathbf{b}_i),c_i\geq0\right\}$

可以发现，从原子范数松弛成一个原子范数，实际上就是把问题从最少分量表示转变为各分量的最小非负加权和。与此同时原子范数还满足范数的非负性、正定性、齐次性及三角不等式( 原子范数为凸，故根据凸函数性质，满足三角不等式)。然而原子范数中需要考虑的优化项过多，因此可以通过将期望信号分量 $\theta_0$ 以其估计值 $\hat {\theta}_0$ 进行代替，通过多次迭代进行 $\hat {\theta}_0$ 的更新，从而减少优化问题的复杂度。因此上面的干扰信号原子范数可表示为：

$\min\limits_{\mathbf{X},\mathbf{s}}\left\|\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\right\|_{\mathcal{A}}\quad\text{s.t.}\frac{1}{2}\left\|\mathbf{Y}-\mathbf{X}\right\|_{\mathcal{F}}^2\leq\eta$

其等价形式为：

$\min\limits_{\mathbf{X},\mathbf{s}~\mathbf{u}\in\mathcal{C}^M}\inf\limits_{,\mathbf{\Omega}\in\mathcal{C}^{T\times T}}\left(\frac{1}{2}\text{Tr}(\tau(\mathbf{u}))+\frac{1}{2}\text{Tr}(\mathbf{\Omega})\right),s.t.\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}\succeq~\mathbf{0},\frac{1}{2}\left\|\mathbf{Y}-\mathbf{X}\right\|_{\mathcal{F}}^2\le~\eta$

其中

$\tau(\mathbf{u})=\sum_{i=1}^K c_i^2\textbf{a}(\theta_i)\textbf{a}(\theta_i)^H$

为干扰协方差阵， $\mathbf{\Omega}$ 为一保证半正定约束的变量。该部分的思想主要就是通过舒尔补与范德蒙德分解，将一个原子范数约束问题等价为一个半正定约束规划问题，接下来对如何从第一个问题的问题的优化项转化为第二个问题的优化项与第一个问题的约束项作一点解释，首先引入如下两个概念:

范德蒙德分解

对于任意秩 $r \leq N$ 的半正定Hermitian-Toeplitz矩阵 $T(u) \in \mathcal{C}^{N \times N}$ ,可以进行如下范德蒙德分解：

$T = \sum _{k=1}^r p_k a(f_k)a^H(f_k)$

当且仅当时该分解唯一，其中为第k个“频率分量”，其形式就像阵列信号处理中的导向矢量一样，其中的各个元素在相位上的变化趋势时线性的，而为每个频率分量对应的幅度以及相位。

舒尔补条件

假设 $C \in \mathcal{S}^m_{++}$ ，即为一正定矩阵， $A \in \mathcal{S}^n$ ,即为一对称矩阵，则有以下不等式成立：

$S \triangleq~ \left[\begin{array}{cc}A&B\\ B^T&C\\ \end{array}\right] \succeq~ 0 \Leftrightarrow~ S_C\triangleq~ A-~BC^{-1}B^T \succeq~0$

在获知以上两个概念后，首先观察以下矩阵：

$\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}$

则若该矩阵半正定，有：

1、 $\tau(\mathbf{u})$ 为半正定矩阵，存在唯一范德蒙德分解。

2、 $\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T$ 一定位于 $\tau(\mathbf{u})$ 的列空间中

在满足作者给定的条件下该矩阵是半正定的，因此首先对该矩阵尝试作范德蒙德分解，即得到：

$\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}=~\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&0\end{bmatrix}+~\begin{bmatrix}0&0\\ 0&\mathbf{\Omega}\end{bmatrix}\\=~\sum_{i=1}^K c_i^2\begin{bmatrix} \textbf{a}(\theta_i) \\ 0 \end{bmatrix}\begin{bmatrix} \textbf{a}(\theta_i) \\ 0 \end{bmatrix}^H+~\sum_{i=1}^K c_i^2\begin{bmatrix} 0 \\ \phi_i \end{bmatrix}\begin{bmatrix} 0 \\ \phi_i \end{bmatrix}^H=~\sum_{i=1}^K c_i^2\begin{bmatrix} \textbf{a}(\theta_i) \\ \phi_i \end{bmatrix}\begin{bmatrix} \textbf{a}(\theta_i) \\ \phi_i \end{bmatrix}^H$

其中 $\mathbf{a}(\theta_i) \in \mathcal{C}^{M\times 1}$ 为 $\tau(\mathbf{u})$ 范德蒙德分解得到的频率分量，其意义为协方差矩阵中的导向矢量， $\phi_i \in \mathcal{C}^{T \times 1}$ 为范德蒙德分解得到的频率分量，其意义为第个入射信号在每个快拍上的相位信息所构成的向量。则若找到了满足该矩阵半正定条件的秩最小的 $\tau(\mathbf{u})$ ，其秩便是组成 $\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T$ 所需的最少原子数，即实现了相应的原子范数最小化，该问题可被表示为：

$\min _{\mathbf{\Omega},\mathbf{u}} rank(\tau(\mathbf{u})),s.t.\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}\succeq ~0$

该问题非凸，可凸松弛至

$\min _{\mathbf{\Omega},\mathbf{u}} Tr(\tau(\mathbf{u}))+~Tr(\mathbf{\Omega})),s.t.\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}\succeq ~0$

该松弛实际上的意义是范数松弛至范数的相似的表示形式，即最少原子数松弛至最小非负加权和。该问题显然为凸：首先矩阵的迹为一仿射，对于约束条件，显然有以下等价条件成立：

$\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}\succeq~ 0 \Leftrightarrow~ \begin{bmatrix} x \\ y \end{bmatrix}^H\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}\begin{bmatrix} x \\ y \end{bmatrix}\\=~x^T \tau(\mathbf{u})x+~y^H(\mathbf{X}-~\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H x+~x^H(\mathbf{X}-~\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)y+~y^H \mathbf{\Omega} y \geq~ 0,\forall x \in ~\mathcal{C}^{M \times 1}, y \in~ \mathcal{C}^{T \times 1}$

该约束为凸。

至此已经可以通过CVX对该凸问题进行求解。作者随后提供了该凸问题的ADMM的迭代解法，对于SDP类的优化问题，通过ADMM能够降低优化算法的复杂度。关于ADMM在原子范数最小化中的应用可以参考：ADMM算法的应用：降低SDP算法复杂度

首先对原问题

进行等价变换：

$\min\limits_{\mathbf{X},\mathbf{s},\mathbf{u},\mathbf{\Omega}}\left(\frac{1}{2}\text{Tr}(\tau(\mathbf{u}))+\frac{1}{2}\text{Tr}(\mathbf{\Omega})\right)+~\frac{\lambda }{2}\left\|\mathbf{Y}-\mathbf{X}\right\|_{\mathcal{F}}^2,s.t.\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}\succeq~\mathbf{0}$

该问题也等价为作者在文中提出的形式：

$\min\limits_{\mathbf{X},\mathbf{s},\mathbf{u},\mathbf{\Omega}}\frac{\varepsilon}{2}\left(\text{Tr}(\tau(\mathbf{u}))+\text{Tr}(\mathbf{\Omega})\right)+~\frac{1}{2}\left\|\mathbf{Y}-\mathbf{X}\right\|_{\mathcal{F}}^2,s.t.\textbf{Z}=~\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix},\mathbf{Z}\succeq~\mathbf{0}$

从原问题，到原问题的等价问题是以罚函数的形式，将约束条件放进优化条件，相当于对原问题进行松弛。同时，引入罚函数法得到的最优值总比原问题的最优值更小(这部分可参考B站-中科大凸优化-凸优化40、41)。

从原问题的等价问题到作者在文中提出的约束问题，作者首先是将惩罚因子移到了矩阵迹的部分，这就使得该问题的目标从最稀疏原子表示出于接收信号相似的信号到对接收信号以最稀疏的信号进行表示。对于约束函数部分，通过引入新的变量 $\mathbf{Z}$ ，使得半正定矩阵

$\begin{bmatrix}\tau(\mathbf{u})&\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T\\ (\mathbf{X}-\mathbf{a}(\hat{\theta}_0)\mathbf{s}^T)^H&\mathbf{\Omega}\end{bmatrix}$

变为无约束，对于该半正定矩阵的优化问题变为了无约束优化问题。

ADMM求解

对作者在文中提出的约束问题形式进行增广拉格朗日函数求解，增广拉格朗日函数如下所示：

$\mathcal{L}(\mathbf{X},\mathbf{s},\mathbf{\Omega},\mathbf{\Lambda},\mathbf{Z})=~\frac{1}{2}\left\|\mathbf{X}-\mathbf{Y}\right\|_{\mathcal{F}}^2+~\frac{\varepsilon}{2}\left(\text{Tr}(\tau(\mathbf{u}))+\text{Tr}(\mathbf{\Omega})\right)\\+~\left<\mathbf{\Lambda},\textbf{Z}-\left[\begin{matrix}\tau(\textbf{u})&\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T\\ (\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T)^H&\mathbf{\Omega}\end{matrix}\right]\right>+~\frac{\rho}{2}\left\|\textbf{Z}-\begin{bmatrix}\tau(\textbf{u})&\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T\\ (\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T)^H&\mathbf{\Omega}\\ \end{bmatrix}\right\|_{\mathcal{F}}^{2}$

其中 $\mathbf{\Lambda}$ 为拉格朗日乘子， $\rho$ 为增广拉格朗日乘子。接下来便是以优化变量-约束变量-拉格朗日乘子的顺序进行交替优化，如下所示：

$(\mathbf{X}^{t+1},\mathbf{s}^{t+1},\mathbf{u}^{t+1},\mathbf{\Omega}^{t+1})=~\underset{\mathbf{X},\mathbf{s},\mathbf{u},\mathbf{\Omega}}{\text{argmin}}\mathcal{L}(\mathbf{X},\mathbf{s},\mathbf{u},\mathbf{\Omega},\mathbf{\Lambda}^t,\mathbf{Z}^t)$

$\mathbf{Z}^{t+1}=\underset{\mathbf{Z}\succeq\mathbf{0}}{\operatorname{argmin}}\mathcal{L}(\mathbf{X}^{t+1},\mathbf{s}^{t+1},\mathbf{u}^{t+1},\mathbf{\Omega}^{t+1},\mathbf{\Lambda}^t,\mathbf{Z})$

$\begin{aligned}\mathbf{\Lambda}^{t+1}=\mathbf{\Lambda}^t+\rho\left(\mathbf{Z}^{t+1}-\begin{bmatrix}\tau(\mathbf{u}^{t+1})&\mathbf{X}^{t+1}-\mathbf{a}(\theta_0)(\mathbf{s}^{t+1})^T\\ (\mathbf{X}^{t+1}-\mathbf{a}(\theta_0)(\mathbf{s}^{t+1})^T)^H&\mathbf{\Omega}^{t+1}\end{bmatrix}\right)\end{aligned}$

其中上标为其迭代次数。式1与式2的要求均为取到全局最优解，式3则是通过取其一阶偏导为0得到其闭式解。至此为止迭代形式仍有些许复杂，因此可通过对增广拉格朗日函数求导对得出一些变量的闭式更新解，首先有以下公式成立：

$\frac{\partial}{\partial A}({Tr}(AX))=\frac{\partial}{\partial A}({Tr}(XA))=X^T$

$\mathrm{Tr}(A^HB)=\sum_{i j}A_{i j}B_{i j} = \left \langle A,B \right \rangle$

$\frac{\partial \left \| f(\mathbf{X}) \right \|_\mathcal{F}^2}{\partial \mathbf{X}}=2f(\mathbf{X})\frac{\partial f(\mathbf{X})}{\mathbf{X}}$

$\frac{\partial x^H}{\partial x} = 0$

对于一些参与求导的分块矩阵，作者进行了如下处理：

$\mathbf{\Lambda}=\begin{bmatrix}\mathbf{\Lambda}_{M\times M}&\mathbf{\Lambda}_{M\times T}\\ \mathbf{\Lambda}_{T\times M}&\mathbf{\Lambda}_{T\times T}\end{bmatrix},\mathbf{Z}=\begin{bmatrix}\mathbf{Z}_{M\times M}&\mathbf{Z}_{M\times T}\\ \mathbf{Z}_{T\times M}&\mathbf{Z}_{T\times T}\end{bmatrix}$

则上面的全局最优解操作即可重新表示为：

$\mathbf{\Omega}^{t+1} =~ \mathbf{\Omega}\big|_{\frac{\partial \mathcal{L}(\mathbf{X}^{t},\mathbf{s}^{t},\mathbf{u}^{t},\mathbf{\Omega},\mathbf{\Lambda}^t,\mathbf{Z}^t) }{\partial \mathbf{\Omega}}=0}=~\mathbf{Z}^t_{T\times T}+~\frac{1}{\rho}\left(\boldsymbol{\Lambda}^t_{T\times T}-\frac{\varepsilon}{2}\mathbf{I}\right)$

$\textbf{u}^{t+1}=~\textbf{u}\big|_{\frac{\partial \mathcal{L}(\mathbf{X}^{t},\mathbf{s}^{t},\mathbf{u},\mathbf{\Omega}^t,\mathbf{\Lambda}^t,\mathbf{Z}^t) }{\partial \mathbf{u}}=0}=~\frac{1}{\rho}\mathbf{\Upsilon}\left(g(\mathbf{\Lambda}_{M\times M}^t)+\rho g(\mathbf{Z}_{M\times M}^t)-\frac{\varepsilon}{2}\mathbf{e}_1\right)$

其中 $\mathbf{\Upsilon}$ 为一对角阵，其对角线上的元素值满足：

$\mathbf{\Upsilon}_{i,i}=\frac{1}{M-i+1}$

$g\left ( \cdot \right )$ 为一矩阵到向量的线性映射，该映射输出的第 $i$ 个元素表示输入矩阵中所有的行数 $p$ 与列数 $q$ 满足 $p-q+1=i$ 的元素的和。 $\mathbf{e}_1$ 则表示首项为1，其余为0的列向量。至于 $\textbf{u}^{t+1}$ 是如何得到的，将增广拉格朗日函数中关于 $\textbf{u}$ 的项进行展开，有

${\frac{\partial \mathcal{L}(\mathbf{X}^{t},\mathbf{s}^{t},\mathbf{u},\mathbf{\Omega}^t,\mathbf{\Lambda}^t,\mathbf{Z}^t) }{\partial \mathbf{u}}} =~ \frac{\varepsilon }{2}M\mathbf{u}_1-~\sum_{i=1}^{M}\sum_{j=1}^{M}(\mathbf{\Lambda}_{M\times M})_{ij}T(\mathbf{u})_{ij}+~\frac{\rho}{2}\sum_{i=1}^{M}\sum_{j=1}^{M}((\mathbf{Z}_{M\times M})_{ij}-~T(\mathbf{u})_{ij})^2$

结合复Toeplitz阵的特性，对 $\mathbf{u}$ 进行逐项求导，即可得到 $\textbf{u}^{t+1}$ 。

对于 $\mathbf{X}$ 与 $\mathbf{s}$ ，这里作者用了一个很巧妙的代换运算，将两个全局最优解以拼接矩阵的方式进行联合求解。首先提取出增广拉格朗日函数中与 $\mathbf{X}$ 和 $\mathbf{s}$ 相关的部分，有：

$\mathcal{L}(\mathbf{X},\mathbf{s})=~\frac{1}{2}\left \| \mathbf{X}-\mathbf{Y} \right \|_{\mathcal{F}}^2-~2\left \langle \mathbf{\Lambda}_{M \times T},\mathbf{X}-\mathbf{a}(\theta_0)\mathbf{s}^T\right \rangle+~\rho\left \| \mathbf{Z}_{M \times T}-\mathbf{X}+\mathbf{a}(\theta_0)\mathbf{s}^T \right \|_{\mathcal{F}}^2$

通过以下两个等价代换，即可把该双变量增广拉格朗日函数变换成单变量增广拉格朗日函数，即：

$\mathbf{Y}-\mathbf{X} = \mathbf{Y}-{\mathbf{F}}_2\mathbf{\bar{X}}$

$\mathbf{X}-\mathbf{a}\left ( \hat{\theta}_0 \right )=\mathbf{F}_1\mathbf{\bar {X}}$

等价的单变量增广拉格朗日函数为：

$\mathcal{L}(\mathbf{\bar {X}})=~\frac{1}{2}\left \| \mathbf{Y}-{\mathbf{F}}_2\mathbf{\bar{X}} \right \|_{\mathcal{F}}^2-~2\left \langle \mathbf{\Lambda}_{M \times T},\mathbf{F}_1\mathbf{\bar {X}}\right \rangle+~\rho\left \| \mathbf{Z}_{M \times T}-\mathbf{F}_1\mathbf{\bar {X}} \right \|_{\mathcal{F}}^2$

对应的全局最优解为：

$\mathbf{\bar{X}}^{t+1} =~\mathbf{\bar{X}}\big|_{\frac{\partial \mathcal{L}(\mathbf{\bar {X}})}{\partial \mathbf{\bar {X}}}=0}=~\left(2\rho\mathbf{F}_1^H\mathbf{F}_1+\mathbf{F}_2^H\mathbf{F}_2\right)^{-1}\left(\mathbf{F}_2^H\mathbf{Y}+2\mathbf{F}_1^H\mathbf{\Lambda}_{M\times T}^t \right)$

其中 $\textbf{F}_1=[\textbf{I}|-\textbf{a}(\theta_0)]$ ， $\textbf{F}_2=[\textbf{I}\mid\textbf{0}_{M\times1}]$ 。

最后为半正定矩阵 $\mathbf{Z}$ 的更新，增广拉格朗日函数中与 $\mathbf{Z}$ 相关的分量如下所示：

$\mathcal{L}(\mathbf{Z})=~\left \langle \mathbf{\Lambda},\mathbf{Z} \right \rangle+~\frac{\rho}{2}\left\|\textbf{Z}-\begin{bmatrix}\tau(\textbf{u})&\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T\\ (\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T)^H&\mathbf{\Omega}\\ \end{bmatrix}\right\|_{\mathcal{F}}^{2}$

则根据ADMM的Scaled Form(该部分可参考ADMM算法简介)，该函数可被表示为：

$\mathcal{L}(\mathbf{Z})=~\frac{\rho}{2}\left\|\textbf{Z}-\begin{bmatrix}\tau(\textbf{u})&\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T\\ (\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T)^H&\mathbf{\Omega}\\ \end{bmatrix}+ \rho^{-1}\mathbf{\Lambda}\right\|_{\mathcal{F}}^{2}-~\frac{1}{2 \rho}\left \| \mathbf{\Lambda} \right \|_{\mathcal{F}}^2$

则该函数的目标即为寻找使第一项F范数的平方最小化的半正定矩阵，则有：

$\left[\begin{matrix}\tau(\textbf{u}^{t+1})&\textbf{X}^{t+1}-\textbf{a}(\theta_0)(\textbf{s}^{t+1})^T\\ (\textbf{X}^{t+1}-\textbf{a}(\theta_0)(\textbf{s}^{t+1})^T)^H&\textbf{a}(\theta_0)(\textbf{s}^{t+1})^T\\ \end{matrix}\right]-~\frac{1}{\rho}\mathbf{\Lambda}^t=~\sum_{i \in \mathcal{D}}\sigma_i^t\textbf{U}_i^t(\textbf{U}_i^t)^H$

$\textbf{Z}^{t+1}=\sum_{i \in \mathcal{D}}\sigma_i^t\textbf{U}_i^t(\textbf{U}_i^t)^H$

其中 $\mathcal{D}=\{i|\sigma_i^t\geq0\}$ 。至此ADMM算法的全部闭式更新解已获得。

收敛性及误差分析

在无噪声环境下，最小化噪声信号的原子范数总是可行的，因此作者考虑在噪声环境下最小化噪声信号的原子范数可行性。

引理1： $\mathbf{\hat {X}}$ 与 $\mathbf{\hat {s}}$ 为上文中的原子范数最小化问题(文中(14))的最优解，当且仅当以下条件满足：

$\left\|\textbf{Y}-\hat{\textbf{X}}\right\|_{\mathcal{A}}^*\leq\varepsilon$

$\left\langle\mathbf{Y}-\mathbf{\hat{X}},\mathbf{\hat{X}}\right\rangle\ge\varepsilon\left\lVert\mathbf{\hat{X}}-\mathbf{a}(\theta_0)\mathbf{\hat{s}}^T\right\rVert_{\mathcal{A}}$

$(\mathbf{Y}-\mathbf{\hat{X}})^H\textbf{a}(\theta_0)=0$

其中 $\left \langle \cdot \right \rangle$ 为两个矩阵的复内积，定义为 $\left\langle\bar{\mathbf{X}},\mathbf{X}\right\rangle=\mathrm{Re}(\mathrm{Tr}(\mathbf{X}^H\bar{\mathbf{X}}))$ ， $||\mathbf{X}||_{\mathcal{A}}^*$ 为 $||\mathbf{X}||_{\mathcal{A}}$ ，即原子范数的对偶范数，满足范数四要素的范数可视为一凸函数，因此对偶范数的定义如下：

$\text{sup}\left\{z^Tx\mid||z||\leq1\right\}=\sup\limits_{||z||\leq1}z^Tx$

即对于范数 $||x||$ ，其对偶范数是找到一个向量或矩阵 $z^T$ ，使得 $z^T$ 和 $x$ 的内积达到最大，这个最大的内积就是 $||x||$ 的对偶范数。因此原子范数的对偶范数定义为(第二步到第三步，从 $\left \| a \right \|\left \| b \right \| \cos(\theta) = \left \langle a,b \right \rangle$ 的角度上看待)：

$\begin{aligned}\|\mathbf{X}\|_A^*&=\sup\limits_{\left\|\mathbf{\bar{X}}\right\|_A\leq1}\left\langle\mathbf{X},\mathbf{\bar{X}}\right\rangle =\sup\limits_{\theta\in[-\frac{\pi}{2},\frac{\pi}{2}],\|\mathbf{b}\|_2=1}\left\langle\mathbf{X},\mathbf{a}(\theta)\mathbf{b}^T\right\rangle =\sup\limits_{\theta\in[-\frac{\pi}{2},\frac{\pi}{2}]}\left\|\mathbf{X}^H\mathbf{a}(\theta)\right\|_2\\ \end{aligned}$

谈到范数的对偶函数，就不得不提范数的共轭函数。根据共轭函数的定义，范数的共轭函数有以下定义成立(以原子范数为例)：

$f(\mathbf{X}) =~ \left \| \mathbf{X} \right \|_{\mathcal{A}}\rightarrow~ f^{*} (\mathbf{Y})=~\sup\limits_{\mathbf{X} \in domf}\left ( \left \langle \mathbf{X,\mathbf{Y}} \right \rangle - \left \| \mathbf{X} \right \|_{\mathcal{A}}\right )$

范数的共轭函数为指示函数，即若 $f(\mathbf{X}) =~ \left \| \mathbf{X} \right \|_{\mathcal{A}}$ ，其共轭函数为：

$f^*(\mathbf{Y})=\begin{cases}0&||\mathbf{Y}||^*\leq1\\ \infty&\text{otherwise}\end{cases}=I_{\{\mathbf{Z}:||\mathbf{Z}||^*\leq1\}}(\mathbf{Y})$

且对偶范数有以下性质成立：

$\left\langle\mathbf{A},\mathbf{B}\right\rangle\leq\left\|\mathbf{A}\right\|_{\mathcal{A}}\left\|\mathbf{B}\right\|_{\mathcal{A}}^*$

对于该引理的证明，首先考虑引理1，其条件1为(14)的约束项，因此显然成立。随后作者将最小化噪声信号分量的原子范数优化问题转化为其等价的LASSO形式(个人感觉更像是罚函数，赋约束以权值将原问题变成一个无约束问题)，即：

$f(\bar{\mathbf{X}})=\frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_2\bar{\mathbf{X}}\right\|_\mathcal{F}^2+\varepsilon\left\|\mathbf{F}_1\bar{\mathbf{X}}\right\|_\mathcal{A}$

该问题显然为凸((14)的等价形式(15)为凸，因此14为凸)，因此存在最小值 $\mathbf{\hat{\bar{X}}}$ ，使得：

$f(\mathbf{\hat{\bar{X}}}+\alpha(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}}))\geq f(\mathbf{\hat{\bar{X}}})$

展开上式，有：

$\frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_2(\mathbf{\hat{\bar{X}}}+\alpha(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}}))\right\|_\mathcal{F}^2+~\varepsilon\left\|\mathbf{F}_1(\mathbf{\hat{\bar{X}}}+\alpha(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}}))\right\|_\mathcal{A} \geq ~\frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}}\right\|_\mathcal{F}^2+~\varepsilon\left\|\mathbf{F}_1\mathbf{\hat{\bar{X}}}\right\|_\mathcal{A}$

$\frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}}-\alpha \mathbf{F}_2(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}}))\right\|_\mathcal{F}^2+~\varepsilon\left\|\mathbf{F}_1(\mathbf{\hat{\bar{X}}}+\alpha(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}}))\right\|_\mathcal{A} \geq ~\frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}}\right\|_\mathcal{F}^2+~\varepsilon\left\|\mathbf{F}_1\mathbf{\hat{\bar{X}}}\right\|_\mathcal{A}$

注意到：

$\left \| \mathbf{A-B} \right \|_{\mathcal{F}}^2 = \left \| \mathbf{A} \right \|_{\mathcal{F}}^2-2\left \langle \mathbf{A}, \mathbf{B}\right \rangle+\left \| \mathbf{B} \right \|_{\mathcal{F}}^2$

则：

$\frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}}\right\|_\mathcal{F}^2-~\alpha\left \langle\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}} ,\mathbf{F}_2(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}}) \right \rangle+~\frac{1}{2}\left\|\alpha \mathbf{F}_2(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}})\right\|_\mathcal{F}^2+~\varepsilon\left\|\mathbf{F}_1(\mathbf{\hat{\bar{X}}}+\alpha(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}}))\right\|_\mathcal{A} \geq ~\frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}}\right\|_\mathcal{F}^2+~\varepsilon\left\|\mathbf{F}_1\mathbf{\hat{\bar{X}}}\right\|_\mathcal{A}$

$\left\langle\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{X}},\mathbf{F}_2(\mathbf{\bar{X}}-\mathbf{\hat{X}})\right\rangle-~\frac{1}{2}\alpha\left\lVert\mathbf{F}_2(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}})\right\rVert_{\mathcal{F}}^2 \le~\alpha^{-1}\varepsilon\left(\left\|\mathbf{F}_1\left(\hat{\bar{\mathbf{X}}}+\alpha(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}})\right)\right\|_\mathcal{A}-\left\|\mathbf{F}_1\mathbf{\hat{\bar{X}}}\right\|_\mathcal{A}\right)$

又因为原子范数 $\left \| \cdot \right \|_{\mathcal{A}}$ 为凸， $\left \| \mathbf{F}_i \left ( \cdot \right )\right \|_{\mathcal{A}}$ 为仿射与凸函数的嵌套，因此也为凸，则根据凸函数的定义，有：

$\left\|\textbf{F}_1\left(\hat{\bar{\textbf{X}}}+\alpha(\bar{\textbf{X}}-\mathbf{\hat{\bar{X}}})\right)\right\|_\mathcal{A}=~\left\|\textbf{F}_1\left((1-\alpha)\hat{\bar{\textbf{X}}}+\alpha\bar{\textbf{X}}\right)\right\|_\mathcal{A}\leq ~(1-~\alpha)\left\|\textbf{F}_1 \hat{\bar{\textbf{X}}}\right\|_\mathcal{A}+~\alpha\left\|\textbf{F}_1 \bar{\textbf{X}}\right\|_\mathcal{A}$

即得：

$\alpha^{-1}\left(\left\|\mathbf{F}_1\left(\hat{\bar{\mathbf{X}}}+\alpha(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}})\right)\right\|_{\mathcal{A}}-\left\|\mathbf{F}_1\mathbf{\hat{\bar{X}}}\right\|_{\mathcal{A}}\right)\leq~\left\|\textbf{F}_1\bar{\textbf{X}}\right\|_{\mathcal{A}}-~\left\|\textbf{F}_1\hat{\bar{\textbf{X}}}\right\|_{\mathcal{A}} ,\forall \mathbf{\bar{X}},\forall \alpha \in~ (0,1)$

将(39)带入(38)对不等式进行进一步放缩：

$\left\langle\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}},\mathbf{F}_2(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}})\right\rangle-~\frac{1}{2}\alpha\left\lVert\mathbf{F}_2(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}})\right\rVert_{\mathcal{F}}^2 \le ~\varepsilon \left\|\textbf{F}_1\bar{\textbf{X}}\right\|_{\mathcal{A}}-~\varepsilon \left\|\textbf{F}_1\hat{\bar{\textbf{X}}}\right\|_{\mathcal{A}}$

同时令 $\alpha\rightarrow 0$ ，则有：

$\inf_{\mathbf{\bar{X}}} \left ( \varepsilon \left\|\textbf{F}_1\bar{\textbf{X}}\right\|_{\mathcal{A}} -\left \langle \mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}},\mathbf{F}_2\mathbf{\bar{X}}\right \rangle \right )\geq~\varepsilon \left\|\textbf{F}_1\mathbf{\hat{\bar{X}}}\right\|_{\mathcal{A}} -~\left \langle \mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}},\mathbf{F}_2\mathbf{\hat{\bar{X}}}\right \rangle$ 则可推出：

$\varepsilon\left|\left|\textbf{F}_1\mathbf{\bar{X}}\right|\right|_\mathcal{A}\geq ~\left\langle\mathbf{Y}-\mathbf{F}_2\mathbf{\hat{\bar{X}}},\mathbf{F}_2(\mathbf{\bar{X}}-\mathbf{\hat{\bar{X}}})\right\rangle +~\varepsilon \left\|\textbf{F}_1\hat{\bar{\textbf{X}}}\right\|_{\mathcal{A}}$

进一步得到：

$\frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_{2}\mathbf{\bar{X}}\right\|_{\mathcal{F}}^{2}+~\varepsilon\left\|\mathbf{F}_{1}\mathbf{\bar{X}}\right\|_{\mathcal{A}} \geq~ \frac{1}{2}\left\|\mathbf{Y}-\mathbf{F}_2\hat{\bar{\textbf{X}}}+\mathbf{F}_2(\hat{\bar{\textbf{X}}}-\mathbf{\bar{X}})\right\|^2_{\mathcal{F}}+~\left<\mathbf{Y}-\mathbf{F}_2\hat{\bar{\textbf{X}}},\mathbf{F}_2(\mathbf{\bar{X}}-\hat{\bar{\textbf{X}}})\right> +~\varepsilon\left\|\mathbf{F}_{1}\hat{\bar{\textbf{X}}}\right\|_{\mathcal{A}}$

上式右边部分以 $\hat{\bar{\textbf{X}}}$ 代替 $\mathbf{\bar{X}}$ 即证 $\hat{\bar{\textbf{X}}}$ 为 $f(\mathbf{\bar{X}})$ 最优解：

$f(\mathbf{\bar{X}})\geq f(\hat{\bar{\textbf{X}}})$

令 $\textbf{X}_c=\textbf{X}-\textbf{a}(\theta_0)\textbf{s}^T$ ， $\mathbf{\hat{\bar{X}}}=[\mathbf{\hat{X}}^T\mid\mathbf{\hat{s}}]^T$ ，则式(40)可表示为：

$\inf_{\mathbf{\bar{X}}}\left(\varepsilon\|\mathbf{X}_{c}\|_{\mathcal{A}}-\left\langle\mathbf{Y}-\mathbf{\hat{X}},\mathbf{X}_{c}+\mathbf{a}(\theta_{0})\mathbf{s}^{T}\right\rangle\right) \geq ~\varepsilon\left\|\mathbf{\hat{X}}-\mathbf{a}(\theta_0)\hat{\mathbf{s}}^T\right\|_{\mathcal{A}}-~\left\langle\mathbf{Y}-\mathbf{\hat{X}},\mathbf{\hat{X}}\right\rangle$

等价于：

$\begin{aligned}-\sup\limits_{\mathbf{X}_c}\left(\left\langle\mathbf{Y}-\hat{\mathbf{X}},\mathbf{X}_c\right\rangle-\varepsilon\left\lVert\mathbf{X}_c\right\rVert_A\right)-\sup\limits_{\mathbf{s}}\left\langle\mathbf{Y}-\hat{\mathbf{X}},\mathbf{a}(\theta_0)\mathbf{s}^T\right\rangle\end{aligned} \geq~\varepsilon\left\|\hat{\mathbf{X}}-\mathbf{a}(\theta_0)\hat{\mathbf{s}}^T\right\|_{\mathcal{A}}-~\left\langle\mathbf{Y}-\hat{\mathbf{X}},\hat{\mathbf{X}}\right\rangle$

因此根据原子范数的对偶范数的定义及共轭函数的性质(可参考【凸优化笔记3】-凸优化问题、示性函数、共轭函数、对偶范数)，有

$\sup\limits_{\mathbf{X}_c}\left(\left\langle\mathbf{Y}-\mathbf{\hat{X}},\mathbf{X}_c\right\rangle-\left\|\mathbf{X}_c\right\|_{\mathcal{A}}\right)=~I_{\left\{\mathbf{M}\left\|\mathbf{M}\right\|_{\mathbf{A}}\leq1\right\}}(\mathbf{Y}-~\mathbf{\hat{X}}) =~\left\{\begin{matrix}0,\textrm{if}\left\|\mathbf{Y}-\mathbf{\hat{X}}\right\|_{\mathcal{A}}^{*}{\leq1}\\ \infty,~~~~\textrm{otherwise}\end{matrix}\right.$

式17条件1即证。

$\sup\limits_{\mathbf{s}}\Big\langle\mathbf{Y}-~\mathbf{\hat{X}},\mathbf{a}(\theta_0)\mathbf{s}^T\Big\rangle=~\begin{cases}0,\mathrm{if}(\mathbf{Y}-\mathbf{\hat{X}})^H\mathbf{a}(\theta_0)=\mathbf{0}\\ \infty,~~\mathrm{otherwise}\end{cases}$

式(17)条件3即证，将(44)与(45)代入(43)，即得：

$\varepsilon\left\|\mathbf{\hat{X}}-\mathbf{a}(\theta_0)\hat{\mathbf{s}}^T\right\|_{\mathcal{A}}-\left\langle\mathbf{Y}-\mathbf{\hat{X}},{\mathbf{\hat{X}}}\right\rangle\leq 0$

式(17)条件2即证，因此满足式(17)三条件的 $\mathbf{\hat{X}}$ 与 $\mathbf{\hat{s}}$ 即为问题(16)的一个解。

定理1 ：对于优化问题(16)，若噪声的 l_1 原子范数的对偶范数存在上确界 $\mathbb{E}\|\mathbf{N}\|_\mathcal{A}^*\leq\varepsilon$ ，则最优解 $\mathbf{\hat{X}}$ 的均方误差可表示为：

$\frac{1}{T}\mathbb{E}\left\|\mathbf{\hat{X}}-\mathbf{X}^{\star}\right\|_\mathcal{F}^2\leq\frac{1}{T}(2\varepsilon\left\|\mathbf{X}^{\star}\right\|_{\mathcal{A}}+\varepsilon\sqrt{\eta})$

其中 $\mathbf{X}^\star=\mathbf{S}+\mathbf{J}$ 为信号的真实值。

证明：根据式(17)，(19)及定理1的前提 $\left\|\mathbf{\hat{X}}-\mathbf{X}^{\star}\right\|_\mathcal{F}^2$ 可表示为：

$\left\|\hat{\mathbf{X}}-\mathbf{X}^{\star}\right\|_\mathcal{F}^2=\left\langle\mathbf{\hat{X}}-\mathbf{X}^{\star},\mathbf{\hat{X}}-(\mathbf{Y}-\mathbf{N})\right\rangle\begin{aligned}&=\Big\langle\mathbf{\hat{X}},\mathbf{\hat{X}}-\mathbf{Y}\Big\rangle+\Big\langle\mathbf{\hat{X}},\mathbf{N}\Big\rangle-\Big\langle\mathbf{X}^{\star},\mathbf{\hat{X}}-\mathbf{Y}\Big\rangle-\langle\mathbf{X}^{\star},\mathbf{N}\rangle\end{aligned}\\~~~~~~~~~~~~~~~~~\leq \begin{aligned}-\varepsilon\left\|\hat{\mathbf{X}}-\mathbf{a}(\theta_0)\hat{\mathbf{s}}^T\right\|_{\mathcal{A}}+\varepsilon\left\|\hat{\mathbf{X}}\right\|_{\mathcal{A}}+2\varepsilon\left\|\mathbf{X}^{\star}\right\|_{\mathcal{A}}\end{aligned}$

同时最优解须满足式14中的条件：

$\left\|\mathbf{Y}-\mathbf{\hat{X}}\right\|_{\mathcal{F}}^2=\left\|\mathbf{Y}-\mathbf{\hat{X}}_c-\mathbf{a}(\theta_0)\mathbf{\hat{s}}^T\right\|_{\mathcal{F}}^2\le\eta$

根据式(17)条件3，有：

$(\textbf{Y}-\mathbf{\hat{X}})^H\textbf{a}(\theta_0)=(\textbf{Y}-\mathbf{\hat{X}}_c)^H\textbf{a}(\theta_0)-\mathbf{\hat{s}}^*=\textbf{0}$

其中 $\mathbf{\hat{s}}^*$ 为 $\mathbf{\hat{s}}$ 的共轭带入式(48)，有：

$\eta\geq\left\|\mathbf{Y}-\mathbf{\hat{X}}_c-\mathbf{a}(\theta_0)\mathbf{\hat{s}}\right\|_\mathcal{F}^2=\left\|\mathbf{Y}-\mathbf{\hat{X}}_c\right\|_\mathcal{F}^2 - 2\left\|(\mathbf{Y}-\mathbf{\hat{X}}_c)^H\mathbf{a}(\theta_0)\mathbf{\hat{s}}\right\|_\mathcal{F}+||\hat{\mathbf{s}}||_{2}^{2}=\left\|\textbf{Y}-\hat{\textbf{X}}_c\right\|^2_{\mathcal{F}}-\left\|\mathbf{\hat{s}}\right\|^2_2\geq0$

即证 $\|\mathbf{\hat{s}}\|_2^2\le\eta$ 。因此有：

$\|\mathbf{a}(\theta_0)\mathbf{\hat{s}}^T\|_{\mathcal{A}}=\|\mathbf{\hat{s}}\|_2\leq\sqrt{\eta}\Rightarrow \|\hat{\mathbf{X}}-\mathbf{X}^{\star}\|_{\mathcal{F}}^2\leq\varepsilon\sqrt{\eta}+2\varepsilon\|\mathbf{X}^{\star}\|_{\mathcal{A}}$

定理2 ：对于优化问题(16)，若噪声的 l_1 原子范数的对偶范数存在上确界，且真实期望信号满足 $\|\mathbf{s}^{\star}\|_2\le\xi$ ，则干扰信号的最优解满足：

$\frac{1}{T}\mathbb{E}\left\|\hat{\mathbf{X}}_{c}-\mathbf{X}_{c}^{\star}\right\|_{\mathcal{F}}^{2}\leq \frac{1}{T}(\left\|\hat{\mathbf{X}}-\mathbf{X}^{\star}\right\|_\mathcal{F}^2+\left\|\mathbf{\hat{s}}-\mathbf{s}^{\star}\right\|_2^2)\leq \frac{1}{T}(\left\|\hat{\mathbf{X}}-\mathbf{X}^{\star}\right\|_\mathcal{F}^2+\left\|\mathbf{\hat{s}}\right\|_2^2+\left\|\mathbf{s}^{\star}\right\|_2^2)\\~~~~~~~~~~~~~~~~~~~~~~~\leq \frac{1}{T}(2\varepsilon\left\|\mathbf{X}^{\star}\right\|_{\mathcal{A}}+\varepsilon\sqrt{\eta}+\eta+\xi)$

其中 $\mathbf{X}_c^\star=\mathbf{J}=\mathbf{X}^\star-\mathbf{a}(\theta_0)(\mathbf{s}^\star)^T$ 为估计值 $\mathbf{\hat{X}}_c$ 的真实值。则可得到干扰信号与干扰子空间均一致收敛。

期望信号方向再估计

已知上一循环所得到的期望信号方向为 $\hat{\theta}_0$ ，则方向误差可表示为 $\Delta\theta=\theta_0-\hat{\theta}_0$ ，求解该方向误差可表示为以下 l_1 原子范数最小化问题：

$\min\limits_{\mathbf{X},\Delta\theta}\left\|\mathbf{X}-\left(\mathbf{a}(\hat{\theta}_0)+\Delta\theta\mathbf{a}'(\hat{\theta}_0)\right)\mathbf{\hat{s}}^T\right\|_{\mathcal{A}}\text{s.t.}\frac{1}{2}\left\|\mathbf{Y}-\mathbf{X}\right\|_{\mathcal{F}}^2\le\eta$

其中 $\mathbf{a}'(\hat{\theta}_0)$ 为 $\mathbf{a}(\theta)$ 在 $\hat{\theta}_0$ 处的导数(与On-grid到Off-Grid有异曲同工之妙)，其SDP形式为

$\min\limits_{\mathbf{X},\Delta\theta,\mathbf{u},\mathbf{\Omega}2}\left(\text{Tr}(\tau(\mathbf{u}))+\text{Tr}(\mathbf{\Omega})\right)+\frac{1}{2}\left\|\mathbf{Y}-\mathbf{X}\right\|_{\mathcal{F}}^2,\\s.t.\mathbf{Z} =\begin{bmatrix}\tau(\textbf{u})&\textbf{X}-(\textbf{a}(\hat{\theta}_0)+\Delta\theta\textbf{a}'(\hat{\theta}_0))\hat{\textbf{s}}^T\\ \textbf{X}^H-\hat{\textbf{s}}^*(\textbf{a}(\hat{\theta}_0)^H+\Delta\theta\textbf{a}'(\hat{\theta}_0)^H) & \mathbf{\Omega}\end{bmatrix},\textbf{Z}\succeq 0$ 该问题同样可以通过CVX或者ADMM进行求解。

自适应权向量估计

$\mathbf{w}=\frac{\tau(\mathbf{u})^{-1}\mathbf{a}(\hat{\theta}_0)}{\mathbf{a}(\hat{\theta}_0)^H\tau(\mathbf{u})^{-1}\mathbf{a}(\hat{\theta}_0)}$

需要注意的是， $\tau(\mathbf{u})$ 为一欠秩阵，因此直接取逆不可行，可采用以该矩阵的M-P广义逆代替之，有：

$\tau(\mathbf{u})^H\tau(\mathbf{u}) = \mathbf{P}\bar{\mathbf{D}}\mathbf{P}^H$

$\tau(\textbf{u})^{-1}=\textbf{PD}^{-1}\textbf{P}^H\tau(\textbf{u})^H$

在这之后作者又将该算法拓展到了多目标场景，即构造多个期望信号的约束，其思想与MVDR拓展到LCMV类似，因此不再赘述与进行仿真。

算法步骤总结

该迭代原子范数最小化自适应波束形成算法步骤如下所示：

输入：量测数据 $\mathbf{Y}$ ，期望信号估计值 $\hat{\theta}_0^0$

初始化：设定噪声相关量 $\varepsilon$ ， $\eta$

设定迭代变化值 $\delta$ ，进行迭代：

迭代步骤1：对问题(14)进行求解，以获得期望信号的估计值，求解得到的信号估计值 $\mathbf{X}$ 记为 $\mathbf{X}_1$ 。

迭代步骤2：对问题(23)进行求解，以获得期望信号方向估计偏差 $\Delta \theta$ ，求解得到的信号估计值 $\mathbf{X}$ 记为 $\mathbf{X}_2$ 。

迭代步骤3：若 $\left\|\hat{\mathbf{X}}_{1}-\hat{\mathbf{X}}_{2}\right\|_{\mathcal{F}}^2\le\delta$ 满足，跳出循环并输出噪声子空间估计值 $\tau(\mathbf{u})$ 与期望信号方向估计偏差 $\Delta \theta$ ，否则返回迭代步骤1。

获得自适应权向量 $\mathbf{w}$ 。

仿真

部分仿真条件如表所示：

阵元数	10
信噪比	10
扰噪比	10
阵元间距	半波长
信源数	1
干扰数	2
信源到达角	5°
干扰到达角	-30°，45°
信源角度估计误差	[0°,3°]
信号F范数容忍度 $\delta$	0.1
噪声阵F范数约束上限 $\eta$	100

代码如下所示(采用的是CVX求解器，若想要提高求解速度请用ADMM对变量进行逐个迭代求解)：

运行该段代码需要安装MATLAB CVX组件与SDPT3组件，若未安装可能运行报错。

运行需求时间较长，主要原因是对角度估计线性逼近过程中一阶泰勒展开导致的逼近速度过慢。

clc;
clear;
jk = 0;
loop = 50;
SINRO = zeros(1,length(-20:5:40));
snr = 10;
deg_dev_predf = 3;%%预定义的角度最大偏移量
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误差
% ang_mismatch1 = 0;
% ang_mismatch2 = 0;
% ang_mismatch3 = 0; %%三个不同方向的DOA误差
%% 初始化及设定参数
array_num = 10;%%阵元数
snapshot_num = 10;%%快拍数
source_aoa = [5,-30,45];%%信源到达角
tgt_ang = 1;
c = 340;%%波速
f = 1000;%%频率
lambda = c/f;%%波长
d = 0.5*lambda;
inr1 = 10;
inr2 = 10;
source_num = length(source_aoa);%%信源数
sig_nr = [snr,inr1,inr2];%%信号功率dBw
deg_deviate = 3;%%误差角度偏离
source_dev = source_aoa+[ang_mismatch1,ang_mismatch2,ang_mismatch3];
%% 导向矢量
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;
%% 算法
theta_est = [];
tar_ang_est = source_dev(1);
theta_est = [theta_est tar_ang_est];
judge_fro_tol = 0.1;
noise_var_tol = 100;
current_fro = norm(rec_sig,"fro")^2;
while current_fro >= judge_fro_tol
%% iteration step 1
    current_ang = tar_ang_est(end);
    current_sv_0 = exp(-1i*(0:array_num-1)'*pi*sind(current_ang));
    cvx_begin quiet
        cvx_solver sdpt3
        variable Tau1(array_num,array_num) hermitian toeplitz
        variable omega_matrix1(snapshot_num,snapshot_num) hermitian toeplitz
        variable source_sig1(snapshot_num,1) complex
        variable X1(array_num,snapshot_num) complex
        minimize(array_num*Tau1(1,1)+snapshot_num*omega_matrix1(1,1))
        subject to 
            [Tau1 X1-current_sv_0*conj(source_sig1');(X1-current_sv_0*conj(source_sig1'))' omega_matrix1] == hermitian_semidefinite(array_num+snapshot_num);
            norm(rec_sig-X1,"fro") <= sqrt(2*noise_var_tol);
    cvx_end
%% iteration step 2
    der_sv_0 = -1i*pi*(0:array_num-1)'*cosd(current_ang).*current_sv_0;
    cvx_begin quiet
        cvx_solver sdpt3
        variable Tau2(array_num,array_num) hermitian toeplitz
        variable omega_matrix2(snapshot_num,snapshot_num) hermitian toeplitz
        variable delta_theta 
        variable X2(array_num,snapshot_num) complex
        minimize(array_num*Tau2(1,1)+snapshot_num*omega_matrix2(1,1))
        subject to 
            [Tau2 X2-(current_sv_0+delta_theta*der_sv_0)*conj(source_sig1');(X2-(current_sv_0+delta_theta*der_sv_0)*conj(source_sig1'))' omega_matrix2] == hermitian_semidefinite(array_num+snapshot_num);
            norm(rec_sig-X2,"fro") <= sqrt(2*noise_var_tol);
    cvx_end
%% judgment
    current_fro = norm(X1-X2,'fro');
    disp(current_fro);
    current_ang = current_ang+delta_theta;
    disp(delta_theta);
    disp(current_ang);
    tar_ang_est = [tar_ang_est current_ang];
end
[P,D] = eig(Tau2'*Tau2);
inv_Tau = P*inv(D)*P'*Tau2';
a_est = exp(-1i*(0:array_num-1)'*pi*sind(current_ang));
w0 = inv_Tau*a_est/(a_est'*inv_Tau*a_est);
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

仿真运行结果如下所示：

部分参考文献

Stephen Boyd; Neal Parikh; Eric Chu; Borja Peleato; Jonathan Eckstein, Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers

Zai Yang, Jian Li, Petre Stoica, and Lihua Xie. Sparse Methods for Direction-of-Arrival Estimation

B. N. Bhaskar, G. Tang and B. Recht. Atomic Norm Denoising With Applications to Line Spectral Estimation. IEEE Transactions on Signal Processing, vol. 61, no. 23, pp. 5987-5999, Dec.1, 2013

X. Zhang, W. Jiang, K. Huo, Y. Liu and X. Li. Robust Adaptive Beamforming Based on Linearly Modified Atomic-Norm Minimization With Target Contaminated Data. IEEE Transactions on Signal Processing, vol. 68, pp. 5138-5151, 2020

量大龙神魔

关注

14
点赞
踩
42

收藏

觉得还不错? 一键收藏
2
评论
波束形成基于原子范数最小化的稳健自适应波束形成

文献 Robust Adaptive Beamforming Based on Linearly Modified Atomic-Norm Minimization With Target Contaminated Data-TSP-2020
复制链接

扫一扫