空时自适应处理STAP算法梳理

1. STAP算法简介

空时自适应处理(STAP)是由L.E. Brennan等人^[1]于1973年提出，利用机载雷达杂波信号存在空时耦合特性。

STAP技术充分利用多通道雷达提供的多个空域通道信息和相干脉冲串提供的时域信息，通过空域和时域二维自适应滤波的方式，实现杂波的有效抑制。STAP 技术的自适应体现在对外部杂波环境的准确感知及应对，其依赖于待检测距离单元( Rangecell Under Test, RUT) 杂波协方差矩阵(Clutter Covariance Matrix, CCM) 的实时获取，而CCM 在实际应用中通常是未知的，需要通过一定数量的独立同分布( Independent Identically Distributed, IID)训练样本最大似然估计得到。根据Reed-Mallett-Brennan准则^[2]，确保STAP输出信杂噪比损失小于$3\text{ dB}$以内所需IID样本数应至少为2倍系统自由度。但机载雷达通常工作在非均匀杂波环境中，难以获得足够的IID训练样本^[3]。

目前，国内外解决上述问题主要两类技术： 一类是将含奇异值样本剔除，使得样本均匀化；另一类是降低均匀样本需求，设计小样本条件下次最优STAP处理器。具体方法划分如图1.1。

• 降维STAP将样本需求由全局系统自由度降至局域系统自由度范畴。
• 降秩STAP基于子空间处理，摒除了由于噪声发散引起的自适应方向图畸变问题，将样本需求降至杂波秩量级，但性能严重依赖于杂波秩估计准确性，且运算量巨大。
• 平滑STAP可利用有限样本的空时平滑获取更多样本，但样本间的强相关性及其固有孔径损失使得该类算法性能不甚理想。
• 直接数据域STAP仅利用RUT数据，消除了非均匀杂波影响，但易受噪声影响且存在空时孔径损失导致性能无法达到次最优，同时空域平滑处理也决定了其仅适用于均匀线/面阵机载雷达。
• 参数化STAP本质为空-时最小二乘有限冲激响应滤波器，在理想条件下可显著降低样本需求，但在实际应用中性能受其模型准确性影响较大。
• 知识辅助STAP利用先验信息所估CCM对杂波进行预白化，以降低后续STAP 处理负担，但如何准确获取及有效利用先验信息仍是当前待解决难题。

S T A P 算 法 概 述	关键问题	涉及的相关方法	注释
	运算量和误差	降维STAP
		降秩STAP
	非均匀杂波	功率非均匀抑制法	机载雷达所照射的范围非常广，能达到几百公里，在这个范围内地面/地形的情况不可能是一成不变的(如陆海交界)。
		非均匀检测器
		直接数据域法
		模型参数化STAP法
		知识辅助的STAP法
		稀疏恢复STAP法
		混合STAP法
	非平稳杂波问题	一维补偿类法	主要是由雷达天线放置的形式和载机飞行方向之间的几何关系，例如若机载相控阵阵面于飞行方向不平行(存在夹角)，此时杂波回波谱会随着距离变化。
		二维补偿类法
		空时内插类法
		权值调整类法
		逆协方差矩阵预测类法
		基于俯仰维预滤波法
		3D-STAP法
	空时自适应检测	基于GLRT准则的STAD	传统STAP主要用于杂波抑制，STAD则是将杂波抑制与目标检测结合到一起考虑。
		基于Rao准则的STAD
		基于Wald准则的STAD

2. STAP信号模型

本部分主要参考文献[12]第二章，机载雷达的相关参数设置如下：

• 机载平台高度为$H$，运动速度为$v$且方向平行于$X$轴；
• 机载雷达的载频为$f_c$(波长为：$\lambda_c = \dfrac{c}{f_c}$)，初始相位为$\phi$；
• 机载平台的运动方向与阵面天线轴向的夹角为$\alpha$；
• 天线为含有$N$个阵元的水平均匀线阵，阵元间距为$d$；
• 一个相干处理间隔(CPI)内发射$M$个脉冲；
• 脉冲重复频率(PRF)为$f_r = \dfrac{1}{T_r}$，脉冲宽度为$T_p$；

此外，文献[15]和文献[16]关于STAP回波的建模写的也特别好，本博客的“2.3 干扰模型”和“2.4 噪声模型”主要参考这两篇文献。

并且文献[16]对于下文公式$(2-21)$中模糊函数部分也有解释，可以看一下。

2.1 窄带阵目标信号模型

假设此时存在一个目标，与雷达视线俯仰角为$\varphi$，相对阵面轴向的方位角为$\theta$，目标相对于天线轴向的入射锥角为$\psi$，目标与载机间的相对径向速度为$v_t$。具体示意图如下。

发射的时域完整波形为：

$$s^{\mathrm t}(t) = a_t u(t) \mathrm e^{\mathrm j 2 \pi f_c t} \mathrm e^{\mathrm j \phi} \tag{2-1}$$

式中上标用字母$t$表示该式为$s^{\mathrm t}(t)$的时域形式，$a_t$表示发射信号的复幅度，$u(t)$为脉冲调制的发射信号基带波形。

考虑到雷达多脉冲相参处理情况，$t = t_k + t_m$为全时间变量，$t_k$为快时间变量，其取值范围为$[0, T_r]$，$t_m = mT_r \quad (m =1,2, \cdots, M)$为慢时间变量。等效相位中心发射第$m$个脉冲经目标散射回到第$n$个阵元的回波信号可以表示为：

$$s^{\mathrm t}_{n, m}(t_k) = a_r u(t_k - \tau_{c, n, m}) \mathrm e^{\mathrm j 2 \pi f_c (t_k - \tau_{c, n, m})} \mathrm e^{\mathrm j \phi}, \quad n = 1,2, \cdots, N \tag{2-2}$$

注意，式$(2-2)$中时间变量为$t_k$ (PS：我之前一直认为这里用全时间变量$t$更合适，其实是错误的，这里$t_k$是表示已经减去了发射脉冲的慢时间间隔，根据下文$\tau_{c,n,m}$的计算可知，$\tau_{c,n,m}$表示的是第$nm$个脉冲从发射到接收的延迟，不是第$nm$个脉冲相对第$1~1$的延迟)，$a_r$表示回波信号幅度，主要由发射功率、系统损耗、距离、天线方向图和目标后向散射特性等决定。复包络信号$u(t)$的宽度为$T_p$ ，当$0 \leq t \leq T_p$时，$u(t)$为发射信号的复包络；当$T_p \leq t \leq T_r$时，$u(t) = 0$。

设第1号阵元为参考阵元，第1个脉冲为参考脉冲，则：

$$\tau_{c, n, m} = \dfrac{2R - (n-1)d \cos \psi - 2(m-1)T_r v_t}{c}$$

其中，$\tau_c = \dfrac{2R}{c}$是参考脉冲由等效相位中心发射经目标散射返回到参考阵元的双程时延，简称目标距离延时；$\tau_m = \dfrac{2(m-1)T_r v_t}{c}$为第$m$个脉冲相对参考脉冲的时延；$\tau_n = \dfrac{(n-1)d \cos \psi}{c}$为第$n$个阵元相对参考阵元的时延。

将接收到的回波信号进行混频，搬移到基带后，第$n$个阵元第$m$个脉冲的回波信号表示为：

$$s^{\mathrm t}_{n, m}(t_k) = a_r u(t_k - \tau_{c, n, m}) \mathrm e^{- \mathrm j 2 \pi f_c \tau_{c, n, m}} \tag{2-3}$$

当相控阵雷达系统为窄带系统时，式中回波信号的复包络在阵元间和脉冲间的相对延时引起的变化可以忽略不计，也即在包络延时中用$\tau_c$代替$\tau_{c, n, m}$，上式表示为：

$$s^{\mathrm t}_{n, m}(t_k) = a_r u(t_k - \tau_{c}) \mathrm e^{- \mathrm j 2 \pi f_c \tau_{c, n, m}} \tag{2-4}$$

接着，对混频后的回波信号进行离散时间采样，以$f_s$为采样频率均匀采样，脉冲间隔$T_r$时间内可得采样点数$L = T_r f_s$。离散采样后，$N$个阵元$M$个脉冲的回波信号离散形式为：

$$\boldsymbol S^{\mathrm t} = [\boldsymbol s^{\mathrm t}_{1,1} \quad \cdots \quad \boldsymbol s^{\mathrm t}_{N,1} \quad \boldsymbol s^{\mathrm t}_{1,2} \quad \cdots \quad \boldsymbol s^{\mathrm t}_{N,2} \quad \cdots \quad \boldsymbol s^{\mathrm t}_{1,M} \quad \cdots \quad \boldsymbol s^{\mathrm t}_{N,M} ] \tag{2-5}$$

其中，$\boldsymbol s^{\mathrm t}_{n,m}$为第$n$个阵元接收到第$m$个脉冲的$L$ 个离散回波信号组成的矢量：

$$\boldsymbol s^{\mathrm t}_{n,m} = \left[s^{\mathrm t}_{n,m}\left(\dfrac{0}{L}T_r\right) \quad \cdots \quad s^{\mathrm t}_{n,m}\left(\dfrac{l}{L}T_r\right) \quad \cdots \quad s^{\mathrm t}_{n,m}\left(\dfrac{L-1}{L}T_r\right) \right]^{\mathrm T} \tag{2-6}$$

其中，$s^{\mathrm t}_{n,m}\left(\dfrac{l}{L}T_r\right)$表示第$n$个阵元第$m$个脉冲回波的第$l$ 个采样点的信号：

$$s^{\mathrm t}_{n,m}\left(\dfrac{l}{L}T_r\right) = a_r u\left(\dfrac{l}{L}T_r - \tau_c\right)\mathrm e^{- \mathrm j 2 \pi f_c \tau_{c, n, m}} \tag{2-7}$$

时间离散后的发射复包络信号的矢量形式为：

$$\boldsymbol u = \left[u\left(\dfrac{0}{L}T_r\right) \quad \cdots \quad u\left(\dfrac{l}{L}T_r\right) \quad \cdots \quad u\left(\dfrac{L-1}{L}T_r\right)\right]^{\mathrm T} \tag{2-8}$$

考虑到雷达的脉冲压缩(脉冲的匹配滤波)通常在频域进行，我们对第$n$个阵元接收到第$m$个脉冲的离散回波信号$\boldsymbol s^{\mathrm t}_{n,m}$做快速傅里叶变换(FFT)处理，可以得到距离频域的信号形式：

$$\boldsymbol s^{\mathrm f}_{n,m} = \mathrm{FFT}(\boldsymbol s^{\mathrm t}_{n,m}) = \boldsymbol F^{\mathrm H}\boldsymbol s^{\mathrm t}_{n,m} \tag{2-9}$$

式中$\boldsymbol s^{\mathrm f}_{n,m}$为$L \times 1$维列向量，其上标用字母$f$ 表示该式为信号$\boldsymbol s^{\mathrm t}_{n,m}$的频域形式。$\boldsymbol F$为$L \times L$维的FFT矩阵，第$l$列矢量$\boldsymbol f_l$表示为：

$$\boldsymbol f_l=\left[\begin{array}{c} 1 \\ \mathrm e^{\mathrm j 2\pi 1 \times (l - L/2)/L} \\ \vdots \\ \mathrm e^{\mathrm j 2\pi (L-1) \times (l - L/2)/L} \end{array}\right], \quad \boldsymbol F=\left[\begin{array}{ccccc} \boldsymbol f_0 & \boldsymbol f_1 & \cdots & \boldsymbol f_{L-1} \end{array}\right] \tag{2-10}$$

第$l_f$号滤波器的中心频率为：$f_l = \dfrac{l_f - L_f/2}{L_f}f_s$，其中$l_f = 0,1,\cdots , L_f -1$表示距离频域滤波器序号，$L_f = L$表示距离频域滤波器个数(为了区分时域采样点数和频域滤波器数，在频域使用$L_f$）。由离散时域信号的DFT时移特性，距离频率为$f_l$的信号可以写成：

$$s^{\mathrm f}_{n,m} = a_r U(f_l) \exp\left(-\mathrm j 2\pi \dfrac{l_f - L_f/2}{L_f} f_s \tau_c\right) \exp\left(-\mathrm j 2\pi f_c \tau_{c,n,m}\right) \tag{2-11}$$

其中，$U(f_l)$是发射脉冲复包络的频谱，是由离散后的发射复包络$u\left(\dfrac{l}{L}T_r\right)$做$L$点FFT处理得到的。

然后在距离频域进行匹配滤波，使用基带信号匹配滤波函数：

$$h(t) = u^*(-t) \tag{2-12}$$

其时间离散后的矢量形式表示为：

$$\boldsymbol h = \left[h\left(\dfrac{0}{L}T_r\right) \quad \cdots \quad h\left(\dfrac{l}{L}T_r\right) \quad \cdots \quad h\left(\dfrac{L-1}{L}T_r\right)\right]^{\mathrm T} \tag{2-13}$$

上式做$L$点FFT处理，得到距离频域的匹配滤波响应：

$$\boldsymbol H(f) = \mathrm{FFT}(\boldsymbol h) = \boldsymbol F^{\mathrm H} \boldsymbol h \tag{2-14}$$

$\boldsymbol H(f)$是$L \times 1$维列向量，表示为：

$$\begin{aligned} \boldsymbol H & =\left[H(f_0), \cdots, H(f_l), \cdots, H(f_{L_f-1}) \right]^{\mathrm T} \\ & =\left[H\left(\dfrac{-L_f/2}{L_f}f_s\right) \quad \cdots \quad H\left(\dfrac{l_f - L_f/2}{L_f}f_s\right) \quad \cdots \quad H\left(\dfrac{L_f/2-1}{L_f}f_s\right)\right]^{\mathrm T} \end{aligned} \tag{2-15}$$

其中，$H(f_l)$为距离频率为$f_l = \dfrac{l_f - L_f/2}{L_f}f_s$的匹配滤波信号：

$$H(f_l) = U^*(f_l) \tag{2-16}$$

对窄带回波信号在频域进行匹配滤波，匹配滤波后的频域输出为：

$$\begin{aligned} x_{n, m, l_{f}}^{\mathrm{f}} & =s_{n, m, l_{f}}^{\mathrm{f}} H_{l_{f}}=s_{n, m}^{\mathrm{f}}\left(f_l\right) H\left(f_l\right) \\ & =a_r U(f_l) U^*(f_l) \exp\left(-\mathrm j 2 \pi\frac{l_{f}-L_{f} / 2}{L_{f}} f_s \tau_c\right) \exp\left(- \mathrm j 2 \pi f_c \tau_{c, n, m}\right) \end{aligned} \tag{2-17}$$

其中，$x_{n, m, l_{f}}^{\mathrm{f}}$的上标用字母$f $表示该式为信号的频域形式。再对$L$个距离频域输出进行逆傅里叶变换(IFFT)，即：

$$\boldsymbol x^{\mathrm t}_{n,m} = \mathrm{IFFT}(\boldsymbol x_{n, m}^{\mathrm{f}}) = \boldsymbol T^{\mathrm H} \boldsymbol x_{n, m}^{\mathrm{f}} \tag{2-18}$$

其中，$\boldsymbol x_{n, m}^{\mathrm{f}}​$为$L \times 1​$维距离频域输出列向量，$\boldsymbol T​$为$L \times L​$维的IFFT 矩阵，第$l​$列矢量$t_l​$表示为：

$$\boldsymbol t_l=\left[\begin{array}{c} 1 \\ \mathrm e^{-\mathrm j 2\pi 1 \times (l - L/2)/L} \\ \vdots \\ \mathrm e^{-\mathrm j 2\pi (L-1) \times (l - L/2)/L} \end{array}\right], \quad \boldsymbol T=\left[\begin{array}{ccccc} \boldsymbol t_0 & \boldsymbol t_1 & \cdots & \boldsymbol t_{L-1} \end{array}\right] = \boldsymbol F^{\mathrm H} \tag{2-19}$$

$\boldsymbol x^{\mathrm t}_{n,m}$为$L \times 1$维距离时域输出列向量，表示为：

$$\boldsymbol x^{\mathrm t}_{n,m} = \left[x^{\mathrm t}_{n,m, 0} \quad \cdots \quad x^{\mathrm t}_{n,m,l} \quad \cdots \quad x^{\mathrm t}_{n,m,L-1} \right]^{\mathrm T} \tag{2-20}$$

其中第$l$个距离时域的信号输出$x^{\mathrm t}_{n,m, l}$可以表示为：

$$\begin{aligned} & x_{n, m, l}^{\mathrm{t}}=\frac{1}{L} \sum_{l_{f}=0}^{L-1} x_{n, m, l_{f}}^{\mathrm{f}} \exp\left(\mathrm j 2 \pi \dfrac{l_{f}-L / 2}{L} l\right) \\ & =\exp\left(-\mathrm j 2 \pi f_c \tau_{c, n, m}\right) \cdot \frac{1}{L} \sum_{l_{f}=0}^{L-1} a_r U\left(f_l\right) U^*(f_l) \exp\left(-\mathrm j 2 \pi \dfrac{l_{f}-L / 2}{L_{f}} f_s \tau_{\mathrm{c}}\right) \exp \left({\mathrm j 2 \pi \frac{l_{f}-L / 2}{L} l}\right) \\ & \approx \exp\left(-\mathrm j 2 \pi f_c \tau_{c, n, m}\right) \cdot \exp \left(\mathrm j 2 \pi \frac{L}{2}\left(\frac{f_s \tau_{\mathrm{c}}}{L_{f}}-\frac{l}{L}\right)\right) \cdot \xi_0 \sum_{l_{f}=0}^{L-1} \exp\left(-\mathrm j 2 \pi \frac{l_{f}}{L_{f}} f_s \tau_{\mathrm{c}}\right) \exp\left(j 2 \pi \frac{l_{f}}{L} l\right) \\ & =\xi_0 \dfrac{\sin \pi L\left(\dfrac{f_{s} \tau_{c}}{L_{f}}-\dfrac{l}{L}\right)}{\sin \pi\left(\dfrac{f_{s} \tau_{c}}{L_{f}}-\dfrac{l}{L}\right)} \exp\left(\mathrm j 2 \pi \dfrac{L}{2}\left(\dfrac{f_{s} \tau_{c}}{L_{f}}-\dfrac{l}{L}\right)\right)\exp\left(-\mathrm j 2 \pi \frac{L-1}{2}\left(\dfrac{f_{s} \tau_{c}}{L_{f}}-\dfrac{l}{L}\right)\right) \exp\left(-j 2 \pi f_{c} \tau_{c, n, m}\right) \\ & \approx \xi_0 \dfrac{\sin \pi L\left(\dfrac{f_{s} \tau_{c}}{L_{f}}-\dfrac{l}{L}\right)}{\sin \pi\left(\dfrac{f_{s} \tau_{c}}{L_{f}}-\dfrac{l}{L}\right)} \exp\left(-\mathrm j 2 \pi f_c \tau_{\mathrm{c}, n, m}\right) \\ \end{aligned} \tag{2-21}$$

其中，$x_{n, m, l}^{\mathrm{t}}$的上标用字母$t$表示该式为信号的时域形式，$l = 0,1, \cdots, L -1$表示第$l$号距离门。式中指数项为阵元间的延时相位、脉冲间的延时相位和目标距离延时的相位。

$x_{n, m, l}^{\mathrm{t}}$的计算中由两个点需要注意：

① 复指数级数求和的计算，这一部分请参考博客文章记忆常用公式 - 1.4.3 常见其他数列求和 - 博客侦探 - 博客园；

② 式中$\xi_0$的意思，这个可是花费了我大量时间，目前我还是没有搞明白，先暂时放弃，但是目前查到的一些相关的资料，这一部分应该是与“匹配滤波器的输出(模糊函数)”有关，可能日后如果有需要可以再仔细研究 —— 文献[13]、文献[14]。

将目标距离延时的相位和信号输出幅度合并，得到窄带信号匹配滤波后的输出信号复幅度：

$$\eta_l = \xi_0 \dfrac{\sin \pi L \left(\dfrac{f_s \tau_c}{L_f} - \dfrac{l}{L}\right)}{\sin \pi \left(\dfrac{f_s \tau_c}{L_f} - \dfrac{l}{L}\right)} \exp(-\mathrm j 2 \pi f_c \tau_c) \tag{2-22}$$

代入式$(2-21)$整理得第$l$个距离单元回波信号的形式：

$$x_{n, m, l}^{\mathrm{t}} = \eta_l \exp\left(\mathrm j 2 \pi f_c \tau_{n, m}\right) = \eta_l \exp\left(\mathrm j 2 \pi f_c (\tau_{n} + \tau_m)\right) \tag{2-23}$$

观察上式可以发现，当目标出现在第$l$号距离单元时，$\eta_l$取最大值。

简要解释为什么目标出现在第$l$号距离单元时，$\eta_l$取最大值：

仅考虑$\eta_l$的模值项——$\xi_0 \dfrac{\sin \pi L \left(\dfrac{f_s \tau_c}{L_f} - \dfrac{l}{L}\right)}{\sin \pi \left(\dfrac{f_s \tau_c}{L_f} - \dfrac{l}{L}\right)}$在$\left(\dfrac{f_s \tau_c}{L_f} - \dfrac{l}{L}\right) = 0$时候取得最大值。因此若目标在第$l$号距离单元，则$\tau_c = \dfrac{l}{f_s}$，此时$\left(\dfrac{f_s \tau_c}{L_f} - \dfrac{l}{L}\right)$恰好等于0，$\eta_l$的模值取得最大。

对于同一个采样时刻，各阵元各脉冲的目标回波信号复幅度相同，目标峰值会出现在同一个距离单元中。

由式$(2-23)$可以得到窄带相控阵雷达系统的目标空域导向矢量：

$$\boldsymbol a^{(t)} = [1 \quad \mathrm e^{\mathrm j 2 \pi \vartheta} \quad \cdots \quad \mathrm e^{\mathrm j 2 \pi (N-1) \vartheta} ]^{\mathrm T} \tag{2-24}$$

其中，$\vartheta$表示目标的空域频率，与$\tau_n = \dfrac{d \cos \psi}{c} (n-1)$有关：

$$\vartheta=f_c \frac{d \cos \psi}{c} \tag{2-25}$$

同样地，根据式$(2-23)$可以得出目标时域导向矢量：

$$\boldsymbol b^{(t)} = [1 \quad \mathrm e^{\mathrm j 2 \pi \varpi} \quad \cdots \quad \mathrm e^{\mathrm j 2 \pi (M-1) \varpi} ]^{\mathrm T} \tag{2-26}$$

其中，$\varpi$表示目标的归一化多普勒频率，与$\tau_m = \dfrac{2 T_r v_t}{c} (m-1)$有关：

$$\varpi = f_c \frac{2 T_r v_{t}}{c} = f_c \frac{2 v_{t}}{c f_r} \tag{2-27}$$

由此我们可以推导出窄带目标的空时导向矢量：

$$\boldsymbol \hbar^{(t)} = \boldsymbol b^{(t)}(\varpi) \otimes \boldsymbol a^{(t)}(\vartheta) \tag{2-28}$$

其中，$\otimes$表示Kronecker积。则第$l$号距离单元的$NM \times 1$维目标回波数据$\boldsymbol{x}^{(t)}_l$可表示为：

$$\boldsymbol{x}^{(t)}_l=\eta_l \cdot \boldsymbol{\hbar}^{(t)} \tag{2-29}$$

统一在此说明一下，这里用上标$(t)$表示目标target的回波信号，同样下文中会使用上标$(c)$、$(j)$、$(n)$分别表示杂波、干扰、噪声。

2.2 窄带杂波模型

机载相控阵雷达与地面杂波的几何关系如图2.3所示。杂波的建模最初文献来自Ward撰写的文献[14]。

载机运动方向与阵面天线轴向的夹角为$\alpha$，杂波块相对雷达视线的俯仰角为$\varphi$，相对阵面轴向的方位角为$\theta$，此时杂波块相对于天线轴向的入射锥角为$\psi$，载机的飞行速度为$v$ 。将雷达所照射到的所有地面杂波区域，按照等距离环均匀划分，每个距离环内又按照等方位角均匀划分，得到众多的独立小杂波块，将各个小杂波块的回波叠加起来，就得到了完整的杂波回波。

PS：其实建立了前面的运动目标的信号模型，杂波模型基本就类比得到结论即可(因为杂波相当于运动速度为0的目标)。

在上述坐标系中，第$i$个距离环，第$j$个独立杂波块$(\theta_j, \varphi_i)$处杂波的空域频率$\vartheta_{i,j}$可表示为：

$$\vartheta_{i,j} = f_c \dfrac{d \cos \psi_{i,j}}{c} = f_c \dfrac{d \cos \theta_j \cos \varphi_i}{c} \tag{2-30}$$

则有，第$i$个距离环，第$j$个独立杂波块的空域导向矢量为：

$$\boldsymbol a^{(c)}_{i, j} = \left[1 \quad \mathrm e^{\mathrm j 2 \pi \vartheta_{i,j}} \quad \cdots \quad \mathrm e^{\mathrm j 2 \pi (N-1) \vartheta_{i,j}} \right]^{\mathrm T} \tag{2-31}$$

归一化多普勒频率$\varpi_{i,j}$可表示为：

$$\varpi_{i,j} = f_c \dfrac{2vT_r \cos(\theta_j + \alpha) \cos \varphi_i}{c} = f_c \dfrac{2v \cos(\theta_j + \alpha) \cos \varphi_i}{c f_r}\tag{2-32}$$

则第$i$个距离环，第$j$个独立杂波块的时域导向矢量为：

$$\boldsymbol b^{(c)}_{i, j} = \left[1 \quad \mathrm e^{\mathrm j 2 \pi \varpi_{i,j}} \quad \cdots \quad \mathrm e^{\mathrm j 2 \pi (M-1) \varpi_{i,j}} \right]^{\mathrm T} \tag{2-33}$$

由式$(2-31)$和式$(2-33)$可得，第$i$个距离环，第$j$个独立杂波块$(\theta_j, \varphi_i)$处杂波的空时导向矢量：

$$\boldsymbol \hbar^{(c)}_{i, j} = \boldsymbol b^{(c)}_{i, j}(\varpi) \otimes \boldsymbol a^{(c)}_{i, j}(\vartheta) \tag{2-34}$$

位于第$i$个距离门，第$j$个杂波块的杂波回波数据形式为：

$$\boldsymbol x^{(c)}_{i, j} = \kappa^{(c)}_{i, j} \cdot \boldsymbol \hbar^{(c)}_{i, j} \tag{2-35}$$

其中，$\kappa^{(c)}_{i, j}$代表该杂波块的复幅度，并且满足以下关系：

$$\mathbb{E}\left(\left|\kappa^{(c)}_{i, j}\right|^2\right) = \sigma^2_n \xi_{i,j } \tag{2-36}$$

其中，$\mathbb{E}(\cdot)$表示数学期望，$\xi_{i,j}$为第$i$个距离门，第$j$个杂波块的杂噪比(CNR)，$\sigma^2_n$表示噪声功率。

下面解释一个距离模糊次数的概念： 机载相控阵雷达信号带宽为$B$，载机高度为$H$，距离分辨率为$\Delta R = \dfrac{c}{2B}$，若地球曲率和载机高度决定的雷 达最大视距是$R_{\max}$，以距离分辨率为间隔将杂波区域等分成若干距离环，则地面上距离环的总数为： $$ L_{\max} = \dfrac{R_{\max} - H}{c/2B} \tag{2-37} $$ 而实际上雷达的最大不模糊距离为$R_u = \dfrac{c}{2f_r}$，以最大不模糊距离为间隔，将雷达视距范围$H \sim R_{\max}$的所 有距离区域混叠到最大不模糊距离的各个距离环上。雷达的距离模糊次数$N_r$为： $$ N_r = \left \lfloor \dfrac{R_{\max}}{R_u} \right \rfloor + 1 \tag{2-38} $$ 其中，$\left \lfloor \cdot \right \rfloor$表示向下取整。

综上所述，可以得到第$i$个距离单元的杂波回波数据：

$$\begin{aligned} \boldsymbol x^{(c)}_ i = &\sum_{k=1}^{N_r} \sum_{j = 1}^{N_c} \kappa^{(c)}_{k, i, j} \cdot \boldsymbol \hbar^{(c)}_{k, i, j} \\ = & \sum_{k=1}^{N_r} \sum_{j = 1}^{N_c} \kappa^{(c)}_{k,i, j} \cdot\boldsymbol b^{(c)}_{k, i, j}(\varpi) \otimes \boldsymbol a^{(c)}_{k, i, j}(\vartheta) \end{aligned} \tag{2-39}$$

其中，$\kappa^{(c)}_{k,i, j}$表示第$k$次模糊到第$i$个距离单元的第$j$个杂波块的复幅度，$N_c$为每个杂波环内的杂波块数，$\boldsymbol \hbar^{(c)}_{k, i, j}$表示相应的空时导向矢量。

2.3 干扰模型

本小节仅研究典型的噪声压制干扰信号，并对其进行建模，其他的干扰方式将在第三章和第四章作详细描述。噪声压制干扰机通常被置于距离敌方雷达较远的地基平台或机载平台，通过发射大功率、宽频带的噪声干扰信号来遮盖有用信号，尽可能地降低目标信噪比，从而实现对接收机的阻塞作用和对目标的压制作用。根据窄带假设，干扰信号在阵元间的传播时间相较于$1/B$要小得多，所以干扰信号在阵元之间是相关的；当雷达的PRF 明显小于接收机的瞬时带宽时，干扰信号在脉冲之间不相关。因此，干扰信号在空域上看起来像点目标或离散的杂波源，在时域的干扰效果则类似热噪声。

第$m$个脉冲重复周期(Pulse Repetition Interval, PRI)的空间快拍的干扰分量表示为：

$$\boldsymbol x^{(j)}_{m} = \gamma_m \boldsymbol a_j \tag{2-40}$$

其中，$\gamma_m$为第$m$个PRI 的干扰信号幅度，$\boldsymbol a_j$为$N \times 1$维干扰空间导向矢量。

设干扰信号俯仰角$\varsigma$，方位角$\epsilon$，可以得到干扰信号的空时快拍为：

$$\boldsymbol x^{(j)} = \boldsymbol \gamma_j \otimes \boldsymbol a_j \tag{2-41}$$

其中，$\boldsymbol \gamma_j$是服从复高斯分布的干扰幅度，具体形式如下：

$$\boldsymbol \gamma_j = [\gamma_0 \quad \gamma_1 \quad \cdots \quad \gamma_{M-1}]^{\mathrm T} \tag{2-42}$$ $$\boldsymbol a_j = \boldsymbol a_j(\epsilon, \varsigma) = \left[ 1 \quad \mathrm e^{\mathrm j 2 \pi (d \cos \varsigma \sin \epsilon)/\lambda_0} \quad \cdots \quad \mathrm e^{\mathrm j 2 \pi (N-1)(d \cos \varsigma \sin \epsilon)/\lambda_0} \right]^{\mathrm T} \tag{2-43}$$

已知来自不同脉冲的干扰信号之间互不相关，并进一步假设干扰信号在一个CPI内是平稳的，则有：

$$\mathbb{E}[\gamma_{m_1} \gamma_{m_2}^*] = \sigma_n^2 \xi_j \delta(m_1-m_2)\\ \Updownarrow \\ \mathbb{E}[\boldsymbol \gamma_j \boldsymbol \gamma_j^{\mathrm H}] = \sigma^2_n \xi_j \boldsymbol I_M \\ \tag{2-44}$$

可以计算得到其协方差矩阵如下式：

$$\begin{aligned} \mathbf{R}_j & = \mathbb{E}\left\{\boldsymbol{x}^{(j)} \boldsymbol{x}^{(j) \mathrm H}\right\}= \mathbb{E}[\boldsymbol \gamma_j \boldsymbol \gamma_j^{\mathrm H}] \otimes \mathbb{E}[ \boldsymbol a_j \boldsymbol a_j^{\mathrm H}]\\ & =\sigma^2_n \xi_j \boldsymbol{I}_M \otimes \boldsymbol{a}_j \boldsymbol{a}_j^{\mathrm H} \\ & =\boldsymbol{I}_M \otimes \sigma_n^2 \xi_j \boldsymbol{a}_j \boldsymbol{a}_j^{\mathrm H} \\ & =\boldsymbol{I}_M \otimes \mathbf{\Phi}_j \end{aligned} \tag{2-45}$$

其中，$\xi_j$为干扰噪声比(Jamming to Noise Ratio，JNR)。式$(2-45)$的空时协方差矩阵是一个块状对角矩阵，即在该矩阵的对角线上，有$M$个$N \times N$维矩阵，而其余部分均为$0$。

2.4噪声模型

雷达进行目标探测时，接收机内部噪声是主要的限制因素。首先，假设每个阵元的接收机噪声在空间上互不相关；此外，当瞬时带宽大于脉冲重复频率(Pulse Repetition Frequency, PRF)时，单个阵元的接收机噪声还满足时间上的不相关。

令$x_{n, n,m}$表示第$n$个阵元在第$m$个脉冲的噪声采样，则噪声的空间相关性和时间相关性分别表示为：

$$\mathbb{E}[x_{n, n_1,m} x^*_{n, n_2, m}] = \sigma_n^2 \delta(n_1 - n_2) \tag{2-46}$$ $$\mathbb{E}[x_{n,n,m_1} x^*_{n, n, m_2}] = \sigma_n^2 \delta(m_1 - m_2) \tag{2-47}$$

本文仅考虑接收机内部噪声。因为每个阵元有各自的接收机，所以可以认为不同阵元的噪声是不相关的。噪声呈随机分布，所以同一阵元的噪声在时间维也是不相关的。假设噪声信号服从均值为$0$，方差为$\sigma_n^2$的高斯分布。噪声空间的协方差矩阵可表示为：

$$\boldsymbol R_n =\mathbb{E}[\boldsymbol x^{(n)} \boldsymbol x^{(n)\mathrm H}] = \sigma_n^2 \boldsymbol I_M \times \boldsymbol I_N = \sigma_n^2 \boldsymbol I_{MN} \tag{2-48}$$

2.5 窄带阵信号模型适用假设讨论

在传统的窄带相控阵雷达系统中，点目标的方向$\psi$和与载机的相对径向速度$v_t$在CPI 内均是视为近似不变，目标空域导向矢量和时域导向矢量在CPI 内均不变。

3 稀疏STAP算法概述

3.1 SR STAP 基本原理及杂波稀疏性分析

(一) 基本原理

不考虑距离模糊次数，第$i$个距离单元机载雷达接收杂波(+噪声)数据可表示为：

$$\begin{aligned} \boldsymbol y =& \boldsymbol x^{(c)}_ i + \boldsymbol x^{(n)} \\ = & \sum_{j = 1}^{N_c} \kappa^{(c)}_{i, j} \cdot\boldsymbol b^{(c)}_{i, j}(\varpi) \otimes \boldsymbol a^{(c)}_{i, j}(\vartheta) + \boldsymbol x^{(n)} \end{aligned} \tag{3-1}$$

这里再回顾一下各个符号的含义：

• $\boldsymbol x^{(c)}_ i$为第$i$个距离单元的杂波回波数据
  • $\kappa^{(c)}_{i, j}$为该杂波块的回波复幅度；
  • $\boldsymbol b^{(c)}_{ i, j}(\varpi)$为第$i$距离环，第$j$个独立杂波块的时域导向矢量；
  • 归一化多普勒频率$\varpi_{i,j} = f_c \dfrac{2vT_r \cos(\theta_j + \alpha) \cos \varphi_i}{c}$
  • $a^{(c)}_{i, j}(\vartheta)$为第$i$距离环，第$j$个独立杂波块的空域导向矢量；
  • 空域频率$\vartheta_{i,j} = f_c \dfrac{d \cos \theta_j \cos \varphi_i}{c}$
• $ \boldsymbol x^{(n)}$表示随机噪声

由于我们现在只考虑第$i$个杂波距离环内的杂波抑制，因此下标$i$不再标出，并且将下标$j$改为$p$，将信号重写为：

$$\begin{aligned} \boldsymbol y =& \boldsymbol x^{(c)} + \boldsymbol x^{(n)} \\ = & \sum_{p = 1}^{N_c} \kappa^{(c)}_{p} \cdot \boldsymbol \nu^{(c)}_{p}(\varpi, \vartheta) + \boldsymbol x^{(n)} \\ = & \sum_{p = 1}^{N_c} \kappa^{(c)}_{p} \cdot \boldsymbol b^{(c)}_{p}(\varpi) \otimes \boldsymbol a^{(c)}_{p}(\vartheta) + \boldsymbol x^{(n)} \end{aligned} \tag{3-2}$$

其中，$\boldsymbol y \in \mathbb C^{NM \times 1}$，$ \kappa^{(c)}_{p}$第$p$个杂波块对应的复幅度，$\boldsymbol \nu^{(c)}_{p}(\varpi, \vartheta)$为空时二维导向矢量，时域导向矢量$\boldsymbol b^{(c)}_{p}(\varpi)$和空域导向矢量$\boldsymbol a^{(c)}_{p}(\vartheta)$分别为：

$$\boldsymbol b^{(c)}_{p}(\varpi) =\left[1 \quad \mathrm e^{\mathrm j 2 \pi \varpi_{p}} \quad \cdots \quad \mathrm e^{\mathrm j 2 \pi (M-1) \varpi_{p}} \right]^{\mathrm T} \tag{3-3}$$ $$\boldsymbol a^{(c)}_{p}(\vartheta) = \left[1 \quad \mathrm e^{\mathrm j 2 \pi \vartheta_{p}} \quad \cdots \quad \mathrm e^{\mathrm j 2 \pi (N-1) \vartheta_{p}} \right]^{\mathrm T} \tag{3-4}$$

由式$(3-2)$看出，机载雷达回波可由不同波束指向(空频)和不同多普勒频率的回波信号叠加而成。若将空频和多普勒频率分别离散化为$N_S = \rho_S N$ 和$N_D = \rho_D M$($\rho_S$和$\rho_D$分别为空频和多普勒频率离散化倍数，通常$\rho_S, \rho_D \gg 1$)，同时忽略量化误差影响，则待雷达回波信号还可表征为：

$$\boldsymbol y = \sum_{q = 1}^{N_S \cdot N_D} \alpha^{(c)}_q \boldsymbol v^{(c)}_{q}(\varpi, \vartheta) + \boldsymbol x^{(n)} = \boldsymbol{V}^{(c)}\boldsymbol{\alpha}^{(c)} + \boldsymbol x^{(n)} \tag{3-5}$$

其中，$\alpha^{(c)}_q$和$\boldsymbol v^{(c)}_q$分别为空时平面第$q$个网格点对应幅度和空时二维导向矢量。$\boldsymbol V^{(c)}$为空时二维平面所有网格点对应空时二维导向矢量集合，也称为字典：

$$\boldsymbol \alpha^{(c)} = [\alpha_1^{(c)}, \alpha_2^{(c)}, \cdots,\alpha_{N_S N_D}^{(c)}]^{\mathrm T} \in \mathbb C^{N_S N_D \times 1} \tag{3-6}$$ $$\boldsymbol V^{(c)} = [\boldsymbol v^{(c)}_{1}, \boldsymbol v^{(c)}_{2}, \cdots, \boldsymbol v^{(c)}_{N_SN_D}] \in \mathbb C^{NM \times N_S N_D} \tag{3-7}$$

由于$N_S × N_D \gg N × K$，式$(3-5)$为欠定方程． 根据SR理论，如果变量$\boldsymbol \alpha^{(c)}$足够稀疏，则可通过以下约束优化问题获得唯一解：

$$\begin{cases} \min & \| \boldsymbol{\alpha}^{(c)} \|_0 \\ \text { s.t. } & \| \boldsymbol y - \boldsymbol{V}^{(c)}\boldsymbol{\alpha}^{(c)} \|^2_2 \leq \varepsilon \end{cases} \tag{3-8}$$

其中，$\varepsilon$为噪声功率，此外，式$(3-8)$还可表述为如下形式：

$$\begin{cases} \min & \| \boldsymbol y - \boldsymbol{V}^{(c)}\boldsymbol{\alpha}^{(c)} \|^2_2\\ \text { s.t. } & \| \boldsymbol{\alpha} ^{(c)}\|_0 \leq r_s \end{cases} \tag{3-9}$$

其中，$r_s$表示杂波稀疏度。求解式$(3-8)$或式$(3-9)$得到稀疏系数$\boldsymbol \alpha$后，当前主要有两大类方法实现目标检测：

• ① 是利用训练样本求得稀疏系数来计算CCM，进而设计空时滤波器进行滤波处理后再检测；
• ② 是仅利用RUT数据求得稀疏系数，实现目标和杂波的同时超分辨谱估计，进而在该空时二维谱平面进行目标检测处理。

(二) 稀疏性分析

已知机载雷达杂波来向与其相应多普勒频率依从如下关系：

$$f_D = \dfrac{2v}{\lambda} \cos(\theta + \alpha)\cos\psi$$

由该式可看出杂波多普勒频率与空频存在依从关系，即一个空频对应一个多普勒频率( 不考虑距离模糊) ，因此杂波在空时二维平面仅分布于满足该依从关系的脊线及附近区域。也就是说，相比于整个空时二维平面，杂波仅分布于小部分特殊区域(一条直线上)，因此在该平面是稀疏的。

文献[17]分析和讨论了杂波在空时平面的稀疏性问题，并指出理想情况及正侧机载雷达杂波稀疏度$r_s$等于杂波秩$r$，即杂波子空间可由r个正交空时导向矢量完备表述。但该结论是否适用于非理想或非正侧情况，尚无定论。

上述讨论均是基于由空时导向矢量构成的过完备字典，其中字典中各导向矢量又称原子。在字典中包含上述$r$个正交原子情况下，通过设计合理SR算法即可求得相应原子位置以及幅度。

注意：对于非正侧或非理想情况，杂波秩要显著高于理想正侧情况，因此构成杂波子空间的正交基或正交空时二维导向矢量个数也相应增加。基于该判断，与理想正侧情况相比，非正侧或非理想情况下杂波稀疏度势必变大，同时杂波在空时平面的稀疏性将变差。

3.2 研究进展及相关问题

(一) 研究进展

(二) 相关问题讨论

(1) 空时谱估计还是杂波抑制?

目标由于自身尺寸限制，同杂波一样在空时平面也具有稀疏性。因此，不同于杂波抑制后检测目标的方式，若直接采用SR技术将RUT数据进行稀疏表征，得到目标和杂波的角-多普勒超分辨谱，则无需进行杂波抑制处理，而直接在空时二维谱检测目标，该方法称为SR空时谱估计方法。

综上所述，实际中采用SR杂波抑制的方式更为稳健，且有利于发现弱小目标。因此，基于SR类空时处理方法的研究应侧重于杂波抑制类方法。

(2) 单观测样本还是多观测样本?

单观测样本，即仅利用RUT数据进行恢复处理，但受噪声影响恢复性能不甚理想。在假定各观测数据具有相同稀疏结构，即不同观测数据中对应稀疏系数矢量的非零元素位置和个数均相同的前提下，式$(3-8)$在多观测样本条件下可进一步表述为：

$$\begin{cases} \min & \| \boldsymbol{A} \|_{2,0} \\ \text { s.t. } & \| \boldsymbol Y - \boldsymbol{V}^{(c)}\boldsymbol{A}^{(c)} \|_{\mathrm F} \leq \varepsilon \end{cases} \tag{3-10}$$

其中，$\boldsymbol A = \left[\boldsymbol{\alpha}^{(c)}_1, \boldsymbol{\alpha}^{(c)}_2, \cdots, \boldsymbol{\alpha}^{(c)}_L \right]$为稀疏系数矩阵，$\boldsymbol Y$为观测样本矩阵，$L$为观测样本个数。$|| \boldsymbol{A} ||_{2,0}$表示对矩阵各行取$l_2$范数，然后在列向取$l_0$范数；$|| \cdot ||_{\mathrm F}$表示对矩阵取Frobenius范数。

文献[40-44]将各类经典SR算法进行了拓展，并证明了如下结论：即单观测样本下SR求得唯一解的充要条件(不考虑噪声)：

$$r_s < \dfrac{\text{spark}(\boldsymbol V^{(c)})}{2} \tag{3-11}$$

在多观测样本条件下松弛为：

$$r_s = |\mathrm{supp}| < \dfrac{\boldsymbol V^{(c)} - 1 + \mathrm{rank}(\boldsymbol Y)}{2} \tag{3-12}$$

其中，$\mathrm{spark}(·)$表示矩阵的最小线性独立列个数，$\mathrm{rank}(·)$表示矩阵的秩。这也就是说多观测样本对稀疏度的要求显著降低，更有利于获得最优解。该结论同样适用于噪声条件下^[40]。

(3) 白化还是置零?

从杂波抑制方式来看，现有SR STAP方法主要分为两类，即杂波白化和杂波置零。

• 杂波白化主要是基于恢复得到原子及其相应幅度构造CCM，然后再进行自适应滤波处理^{[19],[23],[25]}。
• 杂波置零则利用SR得到杂波分量对应原子构造杂波子空间，再通过正交投影等方式进行杂波滤除^[31]。

与白化类方法相比，杂波置零类方法基于信号间正交特性进行杂波抑制，因此无需原子幅度信息，这部分降低了SR难度；同时可生成较深零深，将杂波处功率抑制为零，而白化类算法的零深取决于该处杂波功率强弱，且仅能将杂波抑制到噪声电平。 然而，杂波置零类算法性能严重依赖于杂波子空间维度信息的准确性。

通常来说，子空间维度的略微过估计并不会造成性能下降，但欠估计却可导致性能严重降低，甚至失效^[45]. 尽管文献[46]给出了理想情况各类配置情况杂波秩估计准则，但在实际应用中，受各类非理想因素影响，杂波秩很难准确得到，导致置零类SR算法的性能严重受损。尽管白化类SR算法依赖于原子幅度信息，但可通过合理设计算法而准确估计，因此更适用于机载雷达实际应用。

(4) 重构算法参数依赖还是不依赖?

通过拉格朗日乘子，式$(3-8)$和式$(3-9)$还可转化为如下优化问题：

$$\min _{\lambda, \boldsymbol \alpha^{(c)}} \dfrac{1}{2} \| \boldsymbol y - \boldsymbol{V}^{(c)} \boldsymbol \alpha^{(c)} \|^2_2 + \lambda \| \boldsymbol \alpha^{(c)} \|_0 \tag{3-13}$$

其中，第一项和第二项分别对应信号的保真度和稀疏度，而正则化参数$\lambda$则用于平衡两者间的相对重要性。由式$(3-8)$、$(3-9)$和$(3-10)$可知，参数$\varepsilon$、$r_s$和$\lambda$的设置对于上述优化问题的求解具有重要意义。实质上，$\lambda$可作为式$(3-8)$和$(3-9)$中约束优化问题的对偶形式中的对偶变量，如果已知噪声强度或稀疏度信息，则可对应求得其最优值。

在机载雷达中，准确噪声功率往往是未知的，特别是在低或中重频条件下。因此，参数设置问题成为影响其SR性能的重要因素。在实际应用中，由于空域幅相误差、杂波起伏和载机偏航等非理想因素的影响，杂波的空时谱展宽或泄漏到其他区域，无论杂波秩还是杂波稀疏度均无法准确获知，因此严重依赖该参数的各类算法在实际应用中往往是不稳健的。

将$l_0$拟范数松弛为$l_1$范数后，可通过典型凸优化方法求解上述优化问题，但运算量巨大，且性能严重依赖于正则化参数。从贝叶斯分析角度考虑，可将$l_1$范数优化方法转化为最大后验概率估计方法。在该框架下，SBL方法将稀疏系数各分量表征为不同方差且相互独立的零均值高斯随机变量，并且其参数化先验分布的期望形式为具有较强稀疏性的学生$t$分布，最终求取稀疏系数的过程转化为其超参数的迭代学习过程^[47]。

与$l_1$范数优化方法相比，SBL方法可获得更稀疏的全局最优解，同时不需要噪声强度先验信息，也无需设置正则化参数，因此也是当前SR STAP 中性能稳健的一类方法[48]。

(5) 非平稳杂波环境下是否适用?

除单基正侧放置均匀线(面)阵外，均存在杂波多普勒随距离变化而变化的情况，即杂波具有非平稳性^[49]。

与表述样本间功率或散射源差异性的非均匀不同，非平稳侧重于表述杂波在不同样本间空时二维分布的差异性。该差异性同样导致各观测样本稀疏结构不同，进而影响多样本SR性能。那么在非平稳杂波条件下，SR-STAP技术是否仍旧可行呢?

实质上，机载雷达非平稳杂波由天线俯仰下副瓣引入，主要成份为近程杂波^[50],[51]，因此可通过设计合理的俯仰维波束形成方式进行提前滤除，实现杂波的平稳化^{[50][51],[52]}。此时各观测样本间具有了相同稀疏结构。

如前所述，非均匀杂波条件下提升STAP性能的主要策略是降低样本需求，提升小样本条件下的算法性能；而非平稳杂波环境下则侧重于剔除非平稳杂波分量，将杂波平稳化。现实应用环境中往往是杂波非平稳性和非均匀性交织在一起，因此通过滤除非平稳杂波再级联SR STAP的方式可有效提升上述背景下的杂波抑制性能。

(6) 干扰条件下是否可行?

STAP技术基于空域和时域二维联合自适应处理，对于干扰抑制具有先天优势，因此利用空时自适应滤波同时抑制杂波和干扰是当前机载雷达的重要技术手段。尽管干扰多是宽频带的，但由于其仅来自有限几个方向，因此在空时二维域仍是稀疏的^[26],[29]

• 对于副瓣压制噪声干扰，直接采用传统STAP处理或者SR STAP处理即可实现对干扰和杂波的同时抑制；
• 对于主瓣压制噪声干扰，传统STAP处理由于主瓣畸变引起目标损失，导致算法性能失效，而如果仅利用RUT数据进行稀疏处理得到干扰、杂波和目标空时二维谱，则可期望实现对目标的有效检测和定位，但前提是目标功率相对较强；对于在时域功率变化剧烈的灵巧干扰，传统STAP处理由于无法获取有效样本往往性能较差，如果采用SR STAP进行处理，则可在极少样本条件下实现对干扰和杂波的有效抑制。

4 未来工作展望

4.1 网格失配问题

目前，SR STAP方法中字典由均匀离散化的空间锥角和多普勒频率所对应空时二维导向矢量构成。尽管这种构建字典的方法简单且易于处理，但不可避免存在网格失配问题，即真实杂波以较低概率落于空时平面网格点。网格失配导致杂波能量泄露到其他原子而引起杂波谱展宽，进而影响后续空时滤波器性能。加密网格可增大杂波落于网格点的概率，但过于密集的网格会造成基字典中相邻原子间相关性太强，从而降低恢复性能^[53]。

已有SR STAP方法的研究多集中于正侧视且设定杂波脊线与网格点划分吻合，并未考虑网格失配所带来的严重杂波估计损失问题。文献[39]研究了SR STAP中网格失配问题，并提出一种知识辅助的非均匀网格划分方法，以消除网格失配带来的不利影响，但方法依赖于载机速度和偏航角度等先验知识的准确性。

2012年Candés等^[54]提出利用全变分范数可从少量时域采样中重构无限精度的连续频率值，该方法直接在连续域上进行稀疏正则化，由于无需对连续参数空间进行离散化，因此能够提供信号的精确SR模型并从根本上消除网格失配问题；Tang等^[55]在此基础上提出无网格压缩感知(Gridless Compressive Sensing, GCS)，通过求解原子范数最小化问题，可从有限随机时域采
样中重构出无限精度的连续频率值，但对频率源间隔要求不低于系统自由度的四分之一。

文献[38]将GCS 理论引入到SR-STAP领域，初步验证了该理论对于消除网格失配的天然优势，但并未从理论上给出GCS适用于连续分布杂波的证明，且仅考虑了理想情况． 在后续的工作中，证明GCS理论适用于连续杂波SR，以及研究基于GCS理论的各类非理想条件下SR-STAP方法是该领域的重要研究方向。

4.2 空域误差问题

在实际工程中，由于通道内放大器的增益不一致导致的阵元幅相误差不可避免。在传统STAP方法中，误差分量隐含在接收回波数据中，因此基于回波数据所构造的CCM中蕴含了误差信息，只要有足够空域系统自由度参与自适应处理，则相应空时滤波器对误差具有较强的自适应补偿能力。而在SR-STAP方法中，由于误差未知，通常直接采用理想空时二维导向矢量来构建字典，因此各原子与真实阵列流形之间存在偏差，且误差越大，偏差越大。如果采用不准确的字典来进行SR处理，所估CCM 势必与真实情况存在较大偏差，从而造成后续滤波性能恶化甚至失效。因此，由阵元通道幅相误差导致的原子失配是SR-STAP应用于实际工程时必须要解决的问题。
在SR源定位应用中^[56]，可基于特定稀疏优化准则对信源方位和阵列幅相误差参数进行联合优化估计，以实现对误差参数的在线实时估计和校正，从而提升误差条件下源定位的稳健性。需要注意的是，当前源定位误差校正均基于离散信源背景。Ma 等^[28]首次将源定位中误差和信源来向联合优化技术移植到联合SR-STAP技术中，但其仿真实验仅基于有限强散射杂波源进行了验证。Sun等^[36]继而提出了误差条件下基于SBL的离散强杂波源抑制方法。Liu 等^[57]提出一种误差自校正子空间STAP算法，但严重依赖于对杂波秩的准确估计． 近期，文献^[58]提出的幅相误差与杂波联合优化估计方法，提升了误差条件下SR STAP方法的杂波抑制性能，但该算法性能严重依赖于正则化参数和惩罚因子的设置。因此，进一步验证连续杂波背景下误差自校正SR STAP技术的可行性，以及在此基础上设计稳健的误差与杂波联合稀疏优化方法是未来解决误差影响问题的重要研究方向。

5 基于DOA估计的On-Grid、Off-Grid、Gridless模型梳理^{[59], [61]}

注：目前看来文献[61]对其中几个问题有一些解释。

5.1 简要稀疏恢复模型

(一) 模型概述

稀疏恢复的模型为：

$$\boldsymbol y = \boldsymbol{Ax} + \boldsymbol e \tag{5-1}$$

其中，$\boldsymbol A \in \mathbb C^{M \times \bar{N}}$称为一个字典，$\boldsymbol A$中的列称为“原子”。$\boldsymbol y \in \mathbb C^{M \times 1}$是观测矢量，$\boldsymbol x \in \mathbb C^{\bar{N} \times 1}$是需要稀疏恢复的信号。$\boldsymbol x$中非零元素$K \ll \bar N$，也即$\boldsymbol y$仅需$\boldsymbol A$中$K$个原子就可线性表示。无噪声条件下，稀疏恢复问题为：

$$\begin{cases} \min\limits_{\boldsymbol x} & || \boldsymbol{x} ||_0 \\ \text { s.t. } & \boldsymbol y = \boldsymbol{Ax} \end{cases} \tag{5-2}$$

(二) 凸松弛

第一个实用稀疏恢复算法是基于凸松弛的——将$l_0$范数松弛$l_1$范数，此时问题为：

$$\begin{cases} \min\limits_{\boldsymbol x} & || \boldsymbol{x} ||_1 \\ \text { s.t. } & \boldsymbol y = \boldsymbol{Ax} \end{cases} \tag{5-3}$$

这个算法被称为基追踪(Basic Pursuit, BP)。由于$l_1$范数是凸的，$(5-3)$可以在多项式时间内求解。

当考虑噪声时，可以解决以下优化问题：

$$\min\limits_{\boldsymbol x} \lambda|| \boldsymbol{x} ||_1+ \dfrac{1}{2} ||\boldsymbol y - \boldsymbol{Ax}||_2^2 \tag{5-4}$$ $$\begin{cases} \min\limits_{\boldsymbol x} & || \boldsymbol{x} ||_1 \\ \text { s.t. } & ||\boldsymbol y - \boldsymbol{Ax}||_2 \leq \eta \end{cases} \tag{5-5}$$

其中，式$(5-4)$中的参数$\lambda>0$为正则化参数，式$(5-5)$中的$\eta \geq ||\boldsymbol e||_2$是噪声能量的上界。式$(5-4)$被称为LASSO(Least Absolute Shrinkage and Selection Operator)问题，式$(5-5)$被称为基追踪去噪问题(Basis Pursuit Denoising, BPDN)问题。当参数$\lambda,\eta \to 0$时，式$(5-4)$和式$(5-5)$都等价与式$(5-3)$。

此外，一种LASSO算法的改进版本为平方根LASSO(Square-Root LASSO, SR-LASSO)，其具体形式为：

$$\min\limits_{\boldsymbol x} \tau|| \boldsymbol{x} ||_1+ ||\boldsymbol y - \boldsymbol{Ax}||_2 \tag{5-6}$$

其中，$\tau>0$是正则化参数。对于LASSO，噪声通常被假设为高斯，并且正则化参数$\lambda$被选择为与噪声的标准偏差成比例，而SR-LASSO只需要对噪声分布的较弱假设，并且$\tau$可以被选择为独立于噪声水平的常数。

$l_1$范数虽然是凸优化，但$l_1$范数不是平滑的，高效求解很困难，到目前为止，在加速计算方面也有一些进展，比如$l_1$-magic算法、内点法、共轭梯度法、Nesterov’s带延拓平滑技术等。

(三) 非凸松弛

一个向量的$l_q(0<q<1)$范数定义为：

$$||\boldsymbol x_q|| = \left(\sum_n |x_n|^q\right)^{\frac{1}{q}} \tag{5-7}$$

可以将$l_0$问题非凸松弛为$l_q$范数，无噪声条件下，相应的稀疏恢复问题如下：

$$\begin{cases} \min\limits_{\boldsymbol x} & || \boldsymbol{x} ||_q^q \\ \text { s.t. } & \boldsymbol y = \boldsymbol{Ax} \end{cases} \tag{5-8}$$

其中，使用$|| \boldsymbol{x} ||_q^q$而不是$|| \boldsymbol{x} ||_q$是为了方便计算。与$l_1$范数相比，$l_q$范数更接近$l_0$范数，因此可以预期式$(5-8)$中的$l_q$优化会比BP产生更好的性能。

一个著名的$l_q$优化算法为FOCUSS(Focal Underdetermined System Solver)，FOCUS是一种迭代加权最小二乘法，无噪条件下，在每次迭代中FOCUS解决以下加权最小二乘问题：

$$\begin{cases} \min\limits_{\boldsymbol x} & \displaystyle\sum_n w_n |x_n|^2\\ \text { s.t. } & \boldsymbol y = \boldsymbol{Ax} \end{cases} \tag{5-9}$$

其中，权重系数$w_n = |x_n|^{q-2}$会根据当前迭代的最新解$\boldsymbol x$进行更新。注意，式$(5-9)$可以以closed form求解，因此可以通过适当的初始化来实现迭代算法。该这种算法可以解释为一种保证收敛到局部最小值的大化-最小化(Majorization-Minimization, MM)算法。

当考虑噪声时，可采用如下正则化方法：

$$\min\limits_{\boldsymbol x} \lambda|| \boldsymbol{x} ||_q^q+ \dfrac{1}{2}||\boldsymbol y - \boldsymbol{Ax}||_2^2 \tag{5-10}$$

其中，$\lambda>0$是正则化参数。与FOCUS中相同的主要思想开发了针对式$(5-10)$的正则化FOCUSS算法。但关于式$(5-10)$的一个难题是参数$\lambda$的选择。尽管一些文献中介绍了几种调整该参数的启发式方法，但据我们所知，在这方面还没有理论结果。

为了绕过参数调整问题，提出一种通过迭代最小化稀疏学习(Sparse Learning via Iterative Minimization, SLIM)的MAP估计方法。假设i.i.d的高斯噪声分布服从$\mathcal{CN}(0, \eta)$，$\boldsymbol x$具有如下先验分布：

$$f(\boldsymbol x) \propto \prod_{n} \exp\left[-\dfrac{2}{q}(|x_n|^q - 1)\right] \tag{5-11}$$

SLIM通过求解如下的$l_q$优化问题来计算MAP估计：

$$\min_{\boldsymbol x} M \log \eta + \eta^{-1} ||\boldsymbol{Ax} - \boldsymbol y||_2^2 + \dfrac{2}{q}||\boldsymbol x||_q^q \tag{5-12}$$

为局部求解式$(5-12)$，SLIM会像“正则化FOCUSS”算法一样，迭代更新$\boldsymbol x$，但与FOCUSS算法不同的是，SLIM算法也会同时根据最新的$\boldsymbol x$迭代更新$\eta$。

(四) 最大似然估计(MLE)

MLE是稀疏恢复的另一种常见方法，与凸松弛和OMP相比，MLE的一个优点是不需要了解噪声水平或稀疏程度。

假设$\boldsymbol x$服从均值为0，协方差为$\boldsymbol P = \mathrm{diag}(\boldsymbol p)$的多元高斯分布，其中$p_n \geq 0, n=1, \cdots, \bar N$，(这可视为$\boldsymbol x$的先验分布)，此外设噪声是方差为$\sigma$的i.i.d高斯分布，则从$\boldsymbol y = \boldsymbol{Ax} + \boldsymbol e$可得，$\boldsymbol y$服从均值为0，协方差为$\boldsymbol R = \boldsymbol{APA}^{\mathrm H} + \sigma \boldsymbol I$的高斯分布，可得$\boldsymbol y$的负对数似然函数为：

$$\mathcal L(\boldsymbol p, \sigma) = \log|\boldsymbol R| + \boldsymbol y^{\mathrm H}\boldsymbol R^{-1}\boldsymbol y \tag{5-13}$$

参数$(\boldsymbol p, \sigma)$可通过最小化$\mathcal L(\boldsymbol p, \sigma)$来估计：

$$\min_{\boldsymbol p, \sigma} \mathcal L(\boldsymbol p, \sigma) = \min_{\boldsymbol p, \sigma}\log|\boldsymbol R| + \boldsymbol y^{\mathrm H}\boldsymbol R^{-1}\boldsymbol y \tag{5-14}$$

一旦求解得到了参数$(\boldsymbol p, \sigma)$，则可以获得稀疏$\boldsymbol x$的后验分布，其均值和协方差为：

$$\begin{aligned} &\boldsymbol \mu = \boldsymbol{\Sigma A}^{\mathrm H}\boldsymbol y \\ \\ &\boldsymbol \Sigma = (\boldsymbol A^{\mathrm H}\boldsymbol A + \sigma \boldsymbol P^{-1})^{-1} \end{aligned} \tag{5-15}$$

稀疏向量$\boldsymbol x$可以估计为其后验均值$\boldsymbol \mu$。MLE一个主要的困难来自求解式$(5-14)$，该式的第一项$\log |\boldsymbol R|$是参数$(\boldsymbol p, \sigma)$的非凸(实际是凹的)函数，目前已提出了不同的方法，如重新加权优化法和稀疏贝叶斯学习法 (SBL)。

1. 重新加权优化法，给定最新的估计值$\boldsymbol R_j$后，采用MM算法在每次迭代时通过其切面$\mathrm{Tr}(\boldsymbol R_j^{-1} \boldsymbol R)$对$ \log |\boldsymbol R| $进行线性化。每次迭代产生的问题都是凸问题，使用一种称为基于协方差的稀疏迭代估计(SPICE)的算法求解。
2. 在稀疏贝叶斯SBL框架内，人们对MLE从不同的角度进行了阐释。特别是为了实现稀疏性，假设$\boldsymbol x$的先验分布是具有促进稀疏性的(稀疏先验)，比如以分层方式构建了$\boldsymbol x$的$\mathrm{Student-t}$分布。有趣的是，虽然方法不同，但是得到的目标函数却也是$\min\limits_{\boldsymbol p, \sigma} \mathcal{L}(\boldsymbol p, \sigma)$，为了优化这个目标函数，采用EM算法，在E-Step中计算$\boldsymbol x$的后验分布，在M-Step中对参数$(\boldsymbol p, \sigma)$进行更新。(注意：若$\boldsymbol x$设置不同的促进稀疏的先验分布，则SBL的目标函数会略有不同)

(五) 稀疏恢复与DOA估计的差异

稀疏表示中的字典通常包含有有限数量的原子，而DOA估计中的参数(角度)是连续值——即对应无限原子。

此外DOA一般有多个快拍，可以利用时间冗余性。

5.2 基于网格(On-Grid)的DOA估计

(一) 数据模型

为了填补连续DOA估计与离散稀疏表示之间的空白，网格稀疏方法简单地假定，连续DOA域可以由一组给定的网格点代替：

$$\bar{\boldsymbol \theta} = (\bar \theta_1, \bar \theta_2, \cdots, \bar \theta_{\bar N}) \tag{5-16}$$

其中，$\bar N \gg M$，从而得到如下的$M \times \bar N$的字典矩阵：

$$\boldsymbol A(\bar{\boldsymbol \theta}) = \left[\boldsymbol a(\bar \theta_1), \cdots, \boldsymbol a(\bar \theta_{\bar N})\right] \tag{5-17}$$

考虑$L$个快拍数据(MMV模型)，则有如下：

$$\boldsymbol Y = \boldsymbol{AX} + \boldsymbol E \tag{5-18}$$

其中，$\boldsymbol X = \left[\boldsymbol x(1), \cdots, \boldsymbol x(L)\right]$是行稀疏矩阵。

(二) $l_{2,0}$优化

矩阵$\boldsymbol X$的$l_{2,0}$范数为：

$$||\boldsymbol X||_{2,0} = \#\left\{n: ||\boldsymbol X_n||_2 > 0\right\} = \#\left\{n: X_n \neq \boldsymbol 0\right\} \tag{5-19}$$

其中，$\boldsymbol X_n$表示$\boldsymbol X$的第$n$行。在无噪情况下，$l_{2,0}$优化问题可以描述为：

$$\begin{cases} \min\limits_{\boldsymbol X} & || \boldsymbol{X} ||_{2,0} \\ \text { s.t. } & \boldsymbol Y = \boldsymbol{AX} \end{cases} \tag{5-20}$$

这是一个NP-Hard难题。

(三) 凸松弛

(1) $l_{2,1}$范数优化

矩阵$\boldsymbol X$的$l_{2,1}$范数为：

$$||\boldsymbol X||_{2,1} = \sum_n ||\boldsymbol X_n||_2 \tag{5-21}$$

在无噪情况下，$l_{2,0}$优化问题紧凸松弛为$l_{2,1}$问题可以描述为：

$$\begin{cases} \min\limits_{\boldsymbol X} & || \boldsymbol{X} ||_{2,1} \\ \text { s.t. } & \boldsymbol Y = \boldsymbol{AX} \end{cases} \tag{5-22}$$

在有噪情况下，可转换为LASSO问题或者BPDN问题：

$$\begin{aligned} & \mathrm{LASSO}: \min\limits_{\boldsymbol X} \lambda|| \boldsymbol{X} ||_{2,1}+ \dfrac{1}{2} ||\boldsymbol Y - \boldsymbol{AX}||_F^2 \\ \\ & \mathrm{BPDN}: \begin{cases} \min\limits_{\boldsymbol X} & || \boldsymbol{X} ||_{2,1} \\ \text { s.t. } & ||\boldsymbol Y - \boldsymbol{AX}||_F \leq \eta \end{cases} \end{aligned} \tag{5-23}$$

其中，$\eta \geq ||\boldsymbol E||_2$是噪声能量的上界。LASSO问题中的参数$\lambda$的选择是十分困难的。

(2) $l_{2,1}-SVD$算法优化

。。。。。。

(四) 非凸松弛

矩阵$\boldsymbol X$的$l_{2,q}$范数为：

$$||\boldsymbol X||_{2,q} = \left(\sum_n ||\boldsymbol X_n||^q_2 \right)^{\frac{1}{q}} \tag{5-24}$$

在无噪情况下，$l_{2,0}$优化问题紧凸松弛为$l_{2,q}$问题可以描述为：

$$\begin{cases} \min\limits_{\boldsymbol X} & || \boldsymbol{X} ||^q_{2,q} \\ \text { s.t. } & \boldsymbol Y = \boldsymbol{AX} \end{cases} \tag{5-25}$$

为了局部求解$(5-25)$问题，将FOCUSS算法拓展到多块拍的M-FOCUSS算法，与单快拍EOCUSS一样，M-FOCUSS算法每次迭代中解决如下加权最小二乘问题：

$$\begin{cases} \min\limits_{\boldsymbol X} & \displaystyle\sum_n w_n ||\boldsymbol X_n||_2^2\\ \text { s.t. } & \boldsymbol Y = \boldsymbol{AX} \end{cases} \tag{5-26}$$

其中，权重系数$w_n = ||\boldsymbol X_n||_2^{q-2}$会根据当前迭代的最新解$\boldsymbol X$进行更新。式$(5-26)$有封闭形式的解，因此可以利用迭代算法。

在有噪情况下，可转换为LASSO：

$$\min\limits_{\boldsymbol X} \lambda|| \boldsymbol{X} ||^q_{2,q}+ \dfrac{1}{2} ||\boldsymbol Y - \boldsymbol{AX}||_F^2 \tag{5-27}$$

为了避免式$(5-27)$中调整正则化参数$\lambda$的困难，将SLIM拓展到多快拍情况，也即噪声是方差为$\eta$的i.i.d高斯分布，$\boldsymbol X$服从如下的先验分布：

$$f(\boldsymbol X) \propto \prod_{n} \exp\left[-\dfrac{2}{q}(||\boldsymbol X_n||_2^q - 1)\right] \tag{5-28}$$

与单快拍类似，SLIM通过求解下述$l_{2,q}$优化问题得到$\boldsymbol X$的MAP：

$$\min_{\boldsymbol x} M \log \eta + \eta^{-1} ||\boldsymbol{AX} - \boldsymbol Y||_F^2 + \dfrac{2}{q}||\boldsymbol X||_{2,q}^q \tag{5-29}$$

使用与M-FOCUSS类似的重新加权技术，可以以封闭形式迭代更新$\boldsymbol X$和$\eta$，从而得到SLIM的多快拍版本。

(五) 稀疏迭代协方差估计SPICE^[62]

SPICE(Sparse Iterative Covariance-based Estimation)是一种用于稀疏信号处理和估计的算法，主要用于估计协方差矩阵。基本思想是利用协方差矩阵的稀疏性来降低算法的复杂度，并通过迭代方法逐步逼近协方差矩阵的真实值。在迭代过程中，SPICE算法利用了一种基于协方差矩阵估计的方法，即通过最小化原始数据和协方差矩阵之间的误差来估计协方差矩阵。在估计协方差矩阵的同时，SPICE算法还实现了稀疏性，即只有少数非零元素的矩阵表示。

(1) 广义最小二乘

为了推导广义最小二乘，做如下统计假设：

$$\left\{\boldsymbol x(1), \cdots, \boldsymbol x(L), \boldsymbol e(1), \cdots, \boldsymbol e(L)\right\} \tag{5-30}$$

是互不相关的，并且满足：

$$\begin{aligned} &\mathbb{E}\left[\boldsymbol e(t) \boldsymbol e^{\mathrm H}(t)\right] = \sigma \boldsymbol I \\ & \mathbb{E}\left[\boldsymbol x(t) \boldsymbol x^{\mathrm H}(t)\right] = \boldsymbol P = \mathrm{diag}(\boldsymbol p) \end{aligned} \tag{5-31}$$

其中，$\sigma \geq 0$，$p_n \geq 0,(n = 1,2,\cdots, \bar N)$是感兴趣的参数。由此可得，观测数据$\{\boldsymbol y(1), \cdots, \boldsymbol y(L)\}$是互不相关的，其协方差为：

$$\boldsymbol R = \mathbb{E}\left[\boldsymbol y(t) \boldsymbol y^{\mathrm H}(t)\right] = \boldsymbol{APA}^{\mathrm H} + \sigma \boldsymbol I = \boldsymbol A' \boldsymbol P' \boldsymbol A'^{\mathrm H} \tag{5-32}$$

其中，$\boldsymbol A’ = [\boldsymbol A, \boldsymbol I] \in \mathbb C^{M \times (\bar N + M)}$，$\boldsymbol P’ = \mathrm{diag}(\boldsymbol p, \sigma I) \in \mathbb C^{(\bar N + M)\times (\bar N + M)}$，该过程示意图如下。

据此，看到$\boldsymbol R$是参数$(\boldsymbol p, \sigma )$的线性函数。记样本协方差矩阵(Sample Covariance Matrix, SCM)为$\tilde{\boldsymbol R} = \dfrac{1}{L} \boldsymbol{YY}^{\mathrm H}$，给定$\tilde{\boldsymbol R}$，为了估计$\boldsymbol R$(实际上是估计参数$\boldsymbol p$和$\sigma$)，考虑广义最小二乘法。

首先向量化$\tilde{\boldsymbol R}$为$\tilde{\boldsymbol r} = \mathrm{vec}(\tilde{\boldsymbol R})$，向量化$\boldsymbol R$为$\boldsymbol r = \mathrm{vec}(\boldsymbol R)$，由于$\tilde{\boldsymbol R}$是${\boldsymbol R}$的无偏估计，因此有：

$$\mathbb E(\tilde{\boldsymbol r}) = \boldsymbol r \tag{5-33}$$

此外可以计算$\tilde{\boldsymbol r}$的协方差矩阵为(参考文档-加权协方差拟合准则推导)：

$$\mathrm{Cov}(\tilde{\boldsymbol r}) = \dfrac{1}{L} {\boldsymbol R}^{\mathrm T} \otimes \boldsymbol R \tag{5-34}$$

在广义最小二乘中最小化下述过程：

$$\begin{aligned} & \quad\dfrac{1}{L}(\widetilde{\boldsymbol{r}}-\mathbb{E} \widetilde{\boldsymbol{r}})^{\mathrm H} \mathrm{Cov}^{-1}(\widetilde{\boldsymbol{r}})(\widetilde{\boldsymbol{r}}-\mathbb{E} \widetilde{\boldsymbol{r}}) \\ =&\quad (\widetilde{\boldsymbol{r}}-\boldsymbol{r})^{\mathrm H}\left[\boldsymbol{R}^{-{\mathrm T}} \otimes \boldsymbol{R}^{-1}\right](\widetilde{\boldsymbol{r}}-\boldsymbol{r}) \\ =&\quad\operatorname{vec}^{\mathrm H}(\widetilde{\boldsymbol{R}}-\boldsymbol{R})\left[\boldsymbol{R}^{-{\mathrm T}} \otimes \boldsymbol{R}^{-1}\right] \operatorname{vec}(\widetilde{\boldsymbol{R}}-\boldsymbol{R}) \\ =&\quad\operatorname{vec}^{\mathrm H}(\widetilde{\boldsymbol{R}}-\boldsymbol{R}) \operatorname{vec}\left\{\boldsymbol{R}^{-1}(\widetilde{\boldsymbol{R}}-\boldsymbol{R}) \boldsymbol{R}^{-1}\right\} \\ =&\quad\operatorname{Tr}\left\{(\widetilde{\boldsymbol{R}}-\boldsymbol{R}) \boldsymbol{R}^{-1}(\widetilde{\boldsymbol{R}}-\boldsymbol{R}) \boldsymbol{R}^{-1}\right\} \\ =&\quad\left\|\boldsymbol{R}^{-\frac{1}{2}}(\widetilde{\boldsymbol{R}}-\boldsymbol{R}) \boldsymbol{R}^{-\frac{1}{2}}\right\|_{\mathrm{F}}^2 \end{aligned} \tag{5-35}$$

式$(5-35)$中的准则具有良好的统计特性；例如，在某些条件下，它提供了感兴趣的参数$(\boldsymbol p, \sigma )$的大快照最大似然ML估计器。不幸的是式$(5-35)$在${\boldsymbol R}$中是非凸的，因此在$(\boldsymbol p, \sigma )$中也是非凸的。因此不能保证它可以在全局范围内最小化。

根据式$(5-35)$的启发，提出一种凸准则：

$$\left|\left|\tilde{\boldsymbol R}^{-1/2}(\tilde{\boldsymbol R} - \boldsymbol R)\tilde{\boldsymbol R}^{-1/2}\right|\right|_F^2 \tag{5-36}$$

其中，式$(5-34)$中的$\mathrm{Cov}(\tilde{\boldsymbol r})$用其一致估计代替，也即$\dfrac{1}{L} \tilde{\boldsymbol R}^{\mathrm T} \otimes \tilde{\boldsymbol R}$。由此得到的估计是一个大样本ML估计值，但它只适用于$\tilde{\boldsymbol R}$是非奇异且$L \geq M$的情况。接下来介绍的SPICE算法依赖于式$(5-35)$或式$(5-36)$。

(2) SPICE算法

在SPICE中，当$\tilde{\boldsymbol R}$是非奇异且$L \geq M$时，采用如下协方差拟合准则：

$$h_1 = \left|\left|{\boldsymbol R}^{-1/2}(\tilde{\boldsymbol R} - \boldsymbol R)\tilde{\boldsymbol R}^{-1/2}\right|\right|_F^2 \tag{5-37}$$

在$L<M$，$\tilde{\boldsymbol R}$是奇异时，使用如下拟合标准：

$$h_2 = \left|\left|{\boldsymbol R}^{-1/2}(\tilde{\boldsymbol R} - \boldsymbol R)\right|\right|_F^2 \tag{5-38}$$

一个简单的计算显示？？？

$$\begin{aligned} h_1 &=\mathrm{Tr}\left(\boldsymbol R^{-1}\tilde{\boldsymbol R}\right) + \mathrm{Tr}\left(\tilde{\boldsymbol R}^{-1}{\boldsymbol R}\right)-2M \\ &=\mathrm{Tr}\left(\tilde{\boldsymbol R}^{1/2}\boldsymbol R^{-1}\tilde{\boldsymbol R}^{1/2}\right) + \sum_{n=1}^{\bar N} \left(\boldsymbol a_n^{\mathrm H} \tilde{\boldsymbol R}^{-1} \boldsymbol a_n\right) + \mathrm{Tr}\left(\tilde{\boldsymbol R}^{-1}\right) \sigma-2M \end{aligned} \tag{5-39}$$

因此基于$h_1$的SPICE的优化问题可以等价地表示为：

$$\min_{\boldsymbol p \succeq 0, \sigma>0} \mathrm{Tr}\left(\tilde{\boldsymbol R}^{1/2}\boldsymbol R^{-1}\tilde{\boldsymbol R}^{1/2}\right) + \sum_{n=1}^{\bar N} \left(\boldsymbol a_n^{\mathrm H} \tilde{\boldsymbol R}^{-1} \boldsymbol a_n\right)p_n + \mathrm{Tr}\left(\tilde{\boldsymbol R}^{-1}\right) \sigma \tag{5-40}$$

注意，上式中目标函数的第一项可以等价于：？？？

$$\mathrm{Tr}\left(\tilde{\boldsymbol R}^{1/2}\boldsymbol R^{-1}\tilde{\boldsymbol R}^{1/2}\right) = \begin{cases} \min & \mathrm{Tr}(\boldsymbol X)\\ \text { s.t. } & \left[\begin{array}{cc} \boldsymbol X& \tilde{\boldsymbol R}^{1/2} \\ \tilde{\boldsymbol R}^{1/2} & \boldsymbol R \end{array}\right] \succeq 0 \end{cases} \tag{5-41}$$

此时它在$\boldsymbol R$中是凸的，也即在参数$(\boldsymbol p, \sigma)$中是凸的，因此$h_1$在参数$(\boldsymbol p, \sigma)$上是凸的。

同样，对于$h_2$：

$$\begin{aligned} h_2 &=\mathrm{Tr}\left(\boldsymbol R^{-1}\tilde{\boldsymbol R^2}\right) + \mathrm{Tr}\left({\boldsymbol R}\right)-2\mathrm{Tr}\left(\tilde{\boldsymbol R}\right) \\ &=\mathrm{Tr}\left(\tilde{\boldsymbol R}\boldsymbol R^{-1}\tilde{\boldsymbol R}\right) + \sum_{n=1}^{\bar N} \left(\boldsymbol ||a_n||_2^2 p_n\right) -2 \mathrm{Tr}\left(\tilde{\boldsymbol R}\right) +2M\sigma \end{aligned} \tag{5-42}$$

此时，对应的优化问题：

$$\min_{\boldsymbol p \succeq 0, \sigma>0}\mathrm{Tr}\left(\tilde{\boldsymbol R}\boldsymbol R^{-1}\tilde{\boldsymbol R}\right) + \sum_{n=1}^{\bar N} \left(\boldsymbol ||a_n||_2^2 p_n\right) +2M\sigma \tag{5-43}$$

与式$(5-40)$类似，式$(5-43)$也是凸的，尽管式$(5-40)$、式$(5-43)$都可以被表示为二阶锥优化(SOCP)或半定规划(SDP)，并且有标准的求解器，但由于维度高，在实践中不能基于这些公式轻松求解。

我们现在介绍SPICE算法来处理上述计算问题。我们关注的是$L \geq M$的情况，但类似的结果也适用于$L<M$的情况。SPICE的主要基础是以下重新构造：

$$\mathrm{Tr}\left(\tilde{\boldsymbol R}^{1/2}\boldsymbol R^{-1}\tilde{\boldsymbol R}^{1/2}\right) = \begin{cases} \min\limits_{\boldsymbol C} & \mathrm{Tr}(\boldsymbol C^{\mathrm H}\boldsymbol P'^{-1} \boldsymbol C)\\ \text { s.t. } & \boldsymbol A' \boldsymbol C = \tilde{\boldsymbol R}^{1/2} \end{cases} \tag{5-44}$$

并表明$\boldsymbol C$的解由下式给出：

$$\boldsymbol C = \boldsymbol P' \boldsymbol A'^{\mathrm H}\boldsymbol R^{-1} \tilde{\boldsymbol R}^{1/2} \tag{5-45}$$

将式$(5-44)$带入到式$(5-40)$中，可得最小化$h_1$等价于：

$$\begin{cases} \min\limits_{\boldsymbol C, \boldsymbol p \succeq 0, \sigma>0} &\mathrm{Tr}(\boldsymbol C^{\mathrm H}\boldsymbol P'^{-1} \boldsymbol C) + \displaystyle\sum\limits_{n=1}^{\bar N} \left(\boldsymbol a_n^{\mathrm H} \tilde{\boldsymbol R}^{-1} \boldsymbol a_n\right)p_n + \mathrm{Tr}\left(\tilde{\boldsymbol R}^{-1}\right) \sigma\\ \text { s.t. } & \boldsymbol A' \boldsymbol C = \tilde{\boldsymbol R}^{1/2} \end{cases} \tag{5-46}$$

基于式$(5-46)$，通过迭代求解$\boldsymbol C$和$(\boldsymbol p, \sigma)$，得出SPICE算法为：

1. 首先使用传统波束形成器来初始化$(\boldsymbol p, \sigma)$；
2. 根据$(\boldsymbol p, \sigma)$的最新值带入式$(5-45)$更新$\boldsymbol C$；
3. 再固定$\boldsymbol C$来更新$(\boldsymbol p, \sigma)$；
4. 重复这一过程直到收敛。

注意：$(\boldsymbol p, \sigma)$是有解析解的在固定$\boldsymbol C$时，为进一步说明，考察：

$$\mathrm{Tr}(\boldsymbol C^{\mathrm H}\boldsymbol P'^{-1} \boldsymbol C) = \sum_{n}^{\bar N} \dfrac{||\boldsymbol C_n||_2^2}{p_n} + \sum_{n=\bar N + 1}^{\bar N +M}\dfrac{||\boldsymbol C_n||_2^2}{\sigma} \tag{5-47}$$

其中，$\boldsymbol C_n$表示$\boldsymbol C$的第$n$行。将式$(5-47)$带入式$(5-46)$可得：

$$\begin{aligned} p_n & =\dfrac{\left\|\boldsymbol{C}_n\right\|_2}{\sqrt{\boldsymbol{a}_n^{\mathrm H}\widetilde{\boldsymbol{R}}^{-1} \boldsymbol{a}_n}}, \quad n=1, \ldots, \bar{N}, \\ \\ \sigma & =\sqrt{\sum_{n=\bar{N}+1}^{\bar{N}+M}\dfrac{\left\|\boldsymbol{C}_n\right\|_2^2}{\operatorname{Tr}\left(\widetilde{\boldsymbol{R}}^{-1}\right)}} . \end{aligned} \tag{5-48}$$

由于问题是凸的，并且目标函数在迭代过程中单调递减，因此SPICE算法有望收敛到全局最小值。

接下来我们将讨论如何在SPICE中利用信号稀疏性和联合稀疏性。将式$(5-48)$插入$(5-46)$中，我们看到SPICE问题等价于：

$$\begin{cases} \min\limits _{\boldsymbol{C}} &\displaystyle\sum\limits_{n=1}^{\bar{N}} \sqrt{\boldsymbol{a}_n^{\mathrm H} \widetilde{\boldsymbol{R}}^{-1} \boldsymbol{a}_n}\left\|\boldsymbol{C}_n\right\|_2+\sqrt{\operatorname{Tr}\left(\widetilde{\boldsymbol{R}}^{-1}\right) \sum_{n=\bar{N}+1}^{\bar{N}+M}\left\|\boldsymbol{C}_n\right\|_2^2}\\ \text { s.t. } & \boldsymbol A^{\prime} \boldsymbol C=\widetilde{\boldsymbol{R}}^{1/2} \end{cases} \tag{5-49}$$

观察到，式$(5-49)$中第一项实际上是$\boldsymbol C$的前$\bar N$行的$l_2$范数的加权和(加权$l_{2,1}$范数)，从而促进了$\boldsymbol C$的行稀疏性，因此可以预期大部分$||\boldsymbol C_n||_2$将等于0，与式$(5-48)$一起看，则大部分$p_n$等于0，从而实现稀疏性。联合稀疏性是通过假设$\boldsymbol X$的每一行中的条目具有相同的方差$p_n$来实现的。

SPICE与单快照情况下的平方根LASSO有关。式$(5-43)$中的SPICE问题等效于：

$$\min_{\boldsymbol x} ||\boldsymbol x||_1 + ||\boldsymbol y - \boldsymbol{Ax}||_2 \tag{5-50}$$

也即相当于平方根LASSO中正则化参数$\tau = 1$的情况。

最后，请注意，式$(5-32)$中的分解在一般情况下并不是唯一的。这一观察结果的直接结果是SPICE算法通常不能提供$(\boldsymbol p, \sigma)$的唯一估计。这个问题将在SPICE的无网格版本中得到解决，

(六) 最大似然估计MLE

在MLE中也可以利用联合稀疏性，方法与SPICE类似。

设$\boldsymbol x(t), t=1,\cdots, L$是均值为0，协方差为$\boldsymbol P = \mathrm{diag}(\boldsymbol p)$的i.i.d的多元高斯分布。同时假设存在方差为$\sigma$的i.i.d的高斯噪声，且$\boldsymbol{X,E}$相互独立。那么观测数据$\boldsymbol y(t), t = 1, \cdots, L$则是i.i.d的0均值、协方差为$\boldsymbol R = \boldsymbol{APA}^{\mathrm H} + \sigma \boldsymbol I$的高斯分布，与$\boldsymbol Y$相关的负对数似然函数(参考文档-深入浅出多维高斯分布的最大似然估计矩阵推导)为：

$$\mathcal L(\boldsymbol p, \sigma) = \log |\boldsymbol R| + \mathrm{Tr}\left(\boldsymbol R^{-1} \tilde{\boldsymbol R}\right) \tag{5-51}$$

由此可通过该式进行参数$(\boldsymbol p, \sigma)$的估计：

$$\min_{\boldsymbol p, \sigma}\mathcal L(\boldsymbol p, \sigma) = \min_{\boldsymbol p, \sigma} \log |\boldsymbol R| + \mathrm{Tr}\left(\boldsymbol R^{-1} \tilde{\boldsymbol R}\right) \tag{5-52}$$

多快照MLE已在SBL或贝叶斯压缩传感框架内进行了研究。为了利用 $\boldsymbol x(t), t = 1, \cdots, L$的联合稀疏性，假定所有的$\boldsymbol x(t)$都有相同的稀疏先验。EM算法也可以通过最小化式$(5-52)$中的目标来进行参数估计。

在网格方法的框架内似乎很难获得对该问题的完全令人满意的解决方案，因为所采用的离散网格点与真正的连续DOA之间总是存在不匹配。

5.3 偏离网格(Off-Grid)的DOA估计

离散网格外的稀疏回复算法仍然需要进行网格划分，但与On-Grid的方法不同的是，DOA估计不在受限于网格上，Off-Grid的方法主要有两个方向：

1. 基于固定网格和网格偏移的联合估计方法；
2. 依赖动态网格的估计；

(一) 固定网格+网格偏移

对于固定网格$\bar{\boldsymbol \theta} = (\bar\theta_1, \cdots, \bar \theta_{\bar N})$，可引入如下的Off-Grid数据模型。假设$\bar{\boldsymbol \theta}$是由均匀间隔的网格点组成，网格间隔为$ r = \bar \theta_{k+1} - \bar \theta_k \propto \dfrac{1}{\bar N}$。则对于任意真实的DOA $\theta_k$，设离其最近的网格点为$\bar \theta_{nk}$，满足$|\theta_k - \bar \theta_{nk}|<\dfrac{r}{2}$，此时使用一阶泰勒展开来近似导向矢量$\boldsymbol a(\theta_k)$：

$$\boldsymbol a(\theta_k) \approx \boldsymbol a(\bar \theta_{nk}) + \boldsymbol b(\bar \theta_{nk})(\theta_k - \bar \theta_{nk}) \tag{5-53}$$

其中，$ \boldsymbol b(\bar \theta_{nk}) = \boldsymbol a’(\bar \theta_{nk})$，类比On-Grid数据模型式$(5-18)$，可得此时数据模型表示为：

$$\boldsymbol Y = \boldsymbol \Phi(\boldsymbol \beta) \boldsymbol X + \boldsymbol E \tag{5-54}$$

其中，$\boldsymbol \Phi(\boldsymbol \beta) = \boldsymbol A + \boldsymbol B \mathrm{diag}(\boldsymbol \beta)$，阵列流形矩阵$\boldsymbol A$的定义不变，具体各部分的定义如下：

$$\begin{aligned} & \boldsymbol A = \left[\boldsymbol a(\bar \theta_1), \cdots, \boldsymbol a(\bar \theta_{\bar N})\right] \\ &\boldsymbol B = \left[\boldsymbol b(\bar \theta_1), \cdots, \boldsymbol b(\bar \theta_{\bar N})\right] \\ & \boldsymbol \beta = \left[\beta_1, \cdots, \beta_{\bar N}\right] \in \left[-\dfrac{r}{2}, \dfrac{r}{2}\right]^{\bar N} \end{aligned} \tag{5-55}$$

且满足：

$$\begin{aligned} x_n(t) & = \begin{cases}s_k(t), & \text { if } \bar{\theta}_n=\bar{\theta}_{n_k} ; \\ 0, & \text { otherwise, }\end{cases} &\beta_n = \begin{cases}\theta_k-\bar{\theta}_{n_k}, & \text { if } \bar{\theta}_n=\bar{\theta}_{n_k} ; \\ 0, & \text { otherwise, }\end{cases} \end{aligned} \tag{5-56}$$

其中$n=1, \cdots, \bar{N}, t=1, \cdots, L$。

由式$(5-54)$可知，Off-Grid模型下，DOA估计问题可以表述为具有不确定参数的稀疏表示。 特别是，一旦可以根据$\boldsymbol Y$估计行稀疏矩阵$\boldsymbol X$和$\boldsymbol \beta$ ，则可以估计DOA。

与On-Grid模型式$(5-18)$相比，式$(5-54)$的Off-Grid模型引入了额外的网格偏移参数$\beta_n, n = 1, \cdots, \bar{N}$。注意，若$\boldsymbol \beta = \boldsymbol 0$，则式$(5-54)$的Off-Grid模型可化简为式$(5-18)$的On-Grid数据模型。式$(5-18)$的On-Grid数据模型可视为真实数据模型的0阶近似，而式$(5-54)$的Off-Grid模型可视为真实数据模型的1阶近似，网格失配可以通过联合估计网格偏移量得到补偿。

基于式$(5-54)$中的Off-Grid模型已经提出了几种联合估计$\boldsymbol X$和$\boldsymbol \beta$的算法，下面重点介绍基于$l_1$优化和SBL的2种算法。

(1) $l_1$优化算法

受标准稀疏恢复的启发，提出了几种$l_1$优化算法来处理Off-Grid中的DOA估计。首先是一种称为稀疏全最小二乘(Sparse Total Least-Squares, STLS)的算法，在单快拍情况下，可解决如下类似LASSO的问题：

$$\min_{\boldsymbol x , \boldsymbol \beta} \lambda_1 ||\boldsymbol x||_1 + \dfrac{1}{2}|| \boldsymbol y - [\boldsymbol A + \boldsymbol B \mathrm{diag}(\boldsymbol \beta)] \boldsymbol x ||_2^2 + \lambda_2 || \boldsymbol \beta||_2^2 \tag{5-57}$$

其中，$\lambda_1, \lambda_2$是正则化参数。但式$(5-57)$中没有使用$\boldsymbol \beta \in \left[-\dfrac{r}{2}, \dfrac{r}{2}\right]^{\bar N}$的先验信息。为了启发式控制$\boldsymbol \beta$的大小，它的功率也被最小化。注意，由于双线性项$\mathrm{diag}(\boldsymbol \beta) \boldsymbol x$，式$(5-57)$中的问题是非凸的。因此，为了求解式$(5-57)$采用一种交替算法，迭代求解$\boldsymbol x$和$\boldsymbol \beta$。此外式$(5-57)$也可以拓展到MMV情况来利用$\boldsymbol X$的联合稀疏性。但这些算法都存在一定的困难——参数$\lambda_1, \lambda_2$的选择。

为了利用先验信息$\boldsymbol \beta \in \left[-\dfrac{r}{2}, \dfrac{r}{2}\right]^{\bar N}$，提出了类似BPDN的算法用于单快拍情况：

$$\begin{cases} \min\limits_{\boldsymbol x, \boldsymbol \beta \in \left[-\dfrac{r}{2}, \dfrac{r}{2}\right]^{\bar N}} & ||\boldsymbol{x}||_1 \\ \text { s.t. } & \left|\left| \boldsymbol y - [\boldsymbol A \quad \boldsymbol B] \left[\begin{array}{c} \boldsymbol x \\ \boldsymbol \beta \odot \boldsymbol x \end{array}\right] \right|\right|_2 \leq \eta \end{cases} \tag{5-58}$$

其中，$\odot$表示Hadamard乘积，对应位置元素相乘(Hadamard乘积也可以使用$\circ$符号表示)，$\eta$与噪声水平有关。 类似式$(5-57)$，式$(5-58)$也是非凸的，因此也可以利用交替算法来单调地减小目标函数。

在上述算法基础上，另一种凸优化算法利用$\boldsymbol x$和$\boldsymbol v = \boldsymbol \beta \odot \boldsymbol x$地联合稀疏性：

$$\min_{\boldsymbol x, \boldsymbol v} \lambda ||[\boldsymbol x \quad \boldsymbol v]||_{2,1} + \dfrac{1}{2}\left|\left| \boldsymbol y - [\boldsymbol A \quad \boldsymbol B] \left[\begin{array}{c} \boldsymbol x \\ \boldsymbol v\end{array}\right] \right|\right|_2^2 \tag{5-59}$$

式$(5-59)$在参数选择合适的情况下，等价于下面的问题：

$$\begin{cases} \min\limits_{\boldsymbol x, \boldsymbol v} & ||[\boldsymbol{x \quad v}]||_{2,1} \\ \text { s.t. } & \left|\left| \boldsymbol y - [\boldsymbol A \quad \boldsymbol B] \left[\begin{array}{c} \boldsymbol x \\ \boldsymbol v \end{array}\right] \right|\right|_2 \leq \eta \end{cases} \tag{5-60}$$

这种方法的优势在于它是凸的，可以在多项式时间内求全局解。但该方法无法利用$\boldsymbol \beta$的先验信息，得到的解$\beta_n = \dfrac{v_n}{x_n}$甚至不正确，为解决这个问题，一种2阶段解决方案被提出：

1. 首先从式$(5-59)$中求解得到$\boldsymbol x$；
2. 再固定$\boldsymbol x$，通过最小化$|| \boldsymbol y - [\boldsymbol A + \boldsymbol B \mathrm{diag}(\boldsymbol \beta)] \boldsymbol x ||_2^2$来求解$\boldsymbol \beta$；

(2) SBL优化算法

相关文献在多快照下的SBL框架内提出一种Off-Grid DOA估计的系统方法，称为离网稀疏贝叶斯推理(OGSBI)。 为了估计附加参数$\boldsymbol \beta$，假设$\beta_n(n=1, \cdots, \bar N)$是i.i.d的均匀分布在区间$\left[-\dfrac{r}{2}, \dfrac{r}{2}\right]$上。 在所得的EM算法中，行稀疏信号$\boldsymbol X$的后验分布可以像标准SBL中一样在期望步骤中计算。 在最大化步骤中，除了更新行稀疏信号的功率$\boldsymbol p$和噪声方差$\sigma$之外，还更新$\boldsymbol \beta$。 与标准SBL 一样，似然度保证单调增加，因此可以获得算法的收敛性。

(二) 动态网格

数据模型被称为动态网格$\bar{\boldsymbol \theta}$，因为网格点$\theta_n$不是固定的：

$$\boldsymbol Y = \boldsymbol A(\bar{\boldsymbol \theta}) \boldsymbol X + \boldsymbol E \tag{5-61}$$

对于这个模型，我们需要联合估计行稀疏矩阵$\boldsymbol X$和网格$\bar{\boldsymbol \theta}$。 一旦获得它们，就使用与$\boldsymbol X$的非零行相对应的$\bar{\boldsymbol \theta}$网格点来估计DOA。由于$\theta_n$是根据数据估计的并且可以是连续 DOA 域中的任何值，因此这种离网数据模型是准确的，并且不会出现网格失配。 然而，由于映射$\boldsymbol a(θ)$的非线性，设计$\boldsymbol X$和$\bar{\boldsymbol \theta}$联合估计的算法是困难的。 请注意，我们将介绍的以下算法被指定为离网方法而不是无网格方法，因为它们仍然涉及网格选择(例如$\bar N$的选择和$\bar{\boldsymbol \theta}$的初始化），这会影响算法的计算速度和准确性。

基于式$(5-61)$中的数据模型，已经提出了几种算法：

• 基于SBL框架

第一类在SBL的框架内。但是，通常利用变分EM算法(或变分贝叶斯推断)来执行稀疏信号和参数估计，而不是像以前那样使用EM算法。原因是稀疏向量$\boldsymbol x$的后验分布通常不能在这里明确计算，并且EM需要该分布，而变分EM不需要。这些算法的主要困难是由于强非线性而更新$\boldsymbol \theta$。由于无法获得闭合形式的解，因此只能使用数值方法。

• 基于$l_1$优化

在单个快照的情况下，作为一个例子，论文[60]使用了一个小的$\bar N \geq K$，并试图通过迭代更新$\boldsymbol x$和$\boldsymbol \theta$来解决以下$l_1$优化问题：

$$\min_{\boldsymbol x,\bar{\boldsymbol \theta}} \lambda ||\boldsymbol x||_{1} + \dfrac{1}{2}\left|\left| \boldsymbol y - \boldsymbol A(\bar{\boldsymbol \theta}) \boldsymbol x \right|\right|_2^2 \tag{5-62}$$

为了避免某些$\bar \theta_n$收敛到相同的值，加入一个额外的(非凸)项$g(\bar{\boldsymbol \theta})$来惩罚位置紧邻的参数：

$$\min_{\boldsymbol x,\bar{\boldsymbol \theta}} \lambda_1 ||\boldsymbol x||_{1} + \dfrac{1}{2}\left|\left| \boldsymbol y - \boldsymbol A(\bar{\boldsymbol \theta}) \boldsymbol x \right|\right|_2^2 + \lambda_2 g(\bar{\boldsymbol \theta}) \tag{5-63}$$

其中，$\lambda_1, \lambda_2$是需要调整的正则化参数项。注意到，式$(5-62)$和式$(5-63)$均是非凸的，因此即使给定$\boldsymbol x$求解$\bar{\boldsymbol \theta}$也是困难的，且参数调整也很困难。此外文献[60]考虑了$l_q(0<q<1)$优化以增强稀疏性，但也存在同样的问题。

为增强稀疏性，相关文献又提出一种类似$l_1$优化算法：

$$\min_{\boldsymbol x,\bar{\boldsymbol \theta}} \sum_{n=1}^{\bar N} \lambda \log\left(|x_n|^2+\varepsilon\right)+ \left|\left| \boldsymbol y - \boldsymbol A(\bar{\boldsymbol \theta}) \boldsymbol x \right|\right|_2^2 \tag{5-64}$$

为了局部求解式$(5-64)$，要迭代更新$\boldsymbol x$和$\bar{\boldsymbol \theta}$。为了以封闭形式求解$\boldsymbol x$，式$(5-64)$中目标函数的第一项被一个二次代理函数代替，该函数保证了目标的下降。然后使用梯度下降法来解$\bar{\boldsymbol \theta}$，尽管难以选择$\lambda$，但被建议将$\lambda$设置为与噪声方差的倒数成比例。

总之，Off-Grid的方法引入了更多要估计的量，使算法复杂化了，且多数算法都涉及了非凸优化，只能保证局部收敛。此外很多算法没有理论保证。

5.4 无网格(Gridless)的DOA估计

(一) 数据模型

重新梳理数据模型，对于$N$阵元的ULA，数据模型为：

$$\boldsymbol Y = \boldsymbol A(\boldsymbol f)\boldsymbol S +\boldsymbol E \tag{5-65}$$

其中，$\boldsymbol f=(f_1, \cdots, f_K), f_k \in \mathbb T = \left(-\dfrac{1}{2}, \dfrac{1}{2}\right]$，$\boldsymbol A(\boldsymbol f) = [\boldsymbol a(f_1), \cdots, \boldsymbol a(f_K)] \in \mathbb C^{N \times K}$。

在单快照情况下，上述问题与线谱(频率)估计问题一致。 由于第一个无网格稀疏方法是针对线谱(频率)估计问题开发的，因此我们在单个快照情况下介绍它们，然后讨论如何通过利用快照的联合稀疏性将它们扩展到多个快照情况。 在此之前，下一小节将介绍一个重要的数学工具——Vandermonde分解。

(二) Toeplitz协方差矩阵的Vandermonde分解

在无噪情况下，数据的协方差矩阵为：

$$\boldsymbol R = \mathbb{E}\left(\boldsymbol y(t) \boldsymbol y^{\mathrm H}(t)\right) = \boldsymbol A(\boldsymbol f)\mathrm{diag}(\boldsymbol p)\boldsymbol A^{\mathrm H}(\boldsymbol f) \tag{5-66}$$

其中，$p_k >0 (k = 1, \cdots,K)$是源信号的功率。可以验证数据的协方差矩阵$\boldsymbol R$是一个Hemitian Toeplitz矩阵：

$$\boldsymbol R = \boldsymbol T(\boldsymbol u) = \left[\begin{array}{cccc} u_1 & u_2 & \cdots & u_N \\ u_2^* & u_1 & \cdots & u_{N-1}\\ \vdots & \vdots & \ddots & \vdots \\ u_N^* & u_{N-1}^* & \cdots & u_1 \end{array}\right] \tag{5-67}$$

其中，$\boldsymbol u = [u_1,u_2, \cdots, u_n] \in \mathbb C^N$。并且$\boldsymbol R$是一个半正定矩阵(PSD)，秩为$K$(在$K<N$的条件下)。

Vandermonde分解定理表明：任意非满秩且半正定的Toeplitz矩阵都可以唯一地分解为如式$(5-66)$所示的形式。这就意味着可以从数据协方差中精确地提取频率。形式上这个定理可以表述为如下形式：

证明请参考半正定Toeplitz矩阵的范德蒙德分解 - CSDN。

在同噪声方差的情况下，协方差$\boldsymbol R$仍是Toeplitz矩阵，此时自然可将$\boldsymbol R$分解为信号协方差与噪声协方差之和：

Toeplitz矩阵的Vandemonde分解是无网格稀疏方法中的一个重要工具，特别是这些方法将线谱(频率)估计问题转化为一个PSD的Toeplitz矩阵的估计问题，一旦矩阵被估计出来，则可以进行Vandemonde分解获得频率。

通过将Toeplitz矩阵释义为数据协方差矩阵(尽管可能因为某些统计假设不成立，导致并不是)，与传统的子空间方法直接从样本协方差矩阵估计频率不同，无网格方法利用更复杂的优化方法来估计数据协方差，利用其特殊结构/性质(Toeplitz性、低秩性、半正定性)有望实现更优的性能。

下面，在单快拍情况下讨论2种无网格稀疏优化方法：

1. 确定性优化算法：① 原子范数；② Hankle范数；
2. 协方差拟合方法：① SPICE的无网格版本；

注：“确定性”是指不对感兴趣的信号做出任何统计假设，相反信号是确定性的。我们在一组预定集合中寻找最稀疏的候选信号，候选信号的选取通过某种稀疏度来衡量。

(三) “确定性”优化算法的通用框架

单快拍情况下的数据：

$$\boldsymbol y = \boldsymbol z + \boldsymbol e = \boldsymbol A(\boldsymbol f) \boldsymbol s + \boldsymbol e \tag{5-70}$$

对于确定性稀疏方法，通常要求解以下形式的优化问题：

$$\begin{cases} \min\limits_{\boldsymbol z} & \mathcal M(\boldsymbol z) \\ \text { s.t. } & \left|\left| \boldsymbol z - \boldsymbol y \right|\right|_2 \leq \eta \end{cases} \tag{5-71}$$

其中，$\eta$与噪声功率有关，$\mathcal{M}(\boldsymbol z)$是稀疏度量，该稀疏度量被定义为使得通过最小化$\mathcal{M}(\boldsymbol z)$来减少用于表示$\boldsymbol z$的原子个数。

代替式$(5-71)$，我们可以解决一个正则化问题：

$$\min_{\boldsymbol z} \lambda \mathcal{M}(\boldsymbol z) + \dfrac{1}{2}\left|\left| \boldsymbol z - \boldsymbol y \right|\right|_2^2 \tag{5-72}$$

下面则是要重点考察如何构造$\mathcal{M}(z)$。

(1) 原子$l_0$范数

定义如下原子集：

$$\mathcal A = \left\{\boldsymbol a(f, \phi) = \boldsymbol a(f) \phi \mid f \in \mathbb T, \phi \in \mathbb C, |\phi| = 1 \right\} \tag{5-73}$$

它是无限精度的，因为$f$可以是任意实数。$\phi$是允许了一个初始相位的不同。 而根据式$(5-70)$，真实信号$\boldsymbol z$是原子集中$K$个原子的线性组合。 而原子$l_0$范数，就是指能组成$\boldsymbol z$的最少所需原子数， 即：

$$\begin{aligned} ||\boldsymbol z||_{\mathcal A, 0} =& \inf_{c_k, f_k, \phi_k} \left\{ \mathcal K \left| \boldsymbol z = \sum_{k=1}^{\mathcal K} c_k \boldsymbol a(f_k, \phi_k), \quad c_k >0, f_k \in \mathbb T, |\phi_k| = 1 \right.\right\} \\ = & \inf_{s_k, f_k} \left\{ \mathcal K \left| \boldsymbol z = \sum_{k=1}^{\mathcal K} s_k \boldsymbol a(f_k), \quad f_k \in \mathbb T \right. \right\} \end{aligned} \tag{5-74}$$

因为我们的目的就是为了恢复出$\boldsymbol A(\boldsymbol f)$，而$\boldsymbol z$可以写成无数种$\mathcal{A}$ 中原子线性组合的形式。 但只有对应所用原子数最少即对应最小原子$l_0$范数时，此时组成$\boldsymbol z$的原子才恰好对应待恢复的$\boldsymbol A(\boldsymbol f)$。因此，下面的任务就是最小化$\boldsymbol z$的原子$l_0$范数$||\boldsymbol z||_{\mathcal A, 0}$。

虽然模型建立出来了，但是容易看到，如何对$||\boldsymbol z||_{\mathcal A, 0}$进行最小化，可以说是完全摸不着头脑。这似乎比$l_0$范数最小化更为抽象。此时就是见证数学魅力的时刻，考虑如下形式的矩阵：

$$\left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right] \succeq 0 \tag{5-75}$$

其中，$x$是一个要被优化的自由变量(这里是标量)，$\boldsymbol T(\boldsymbol u)$的定义如之前，是一个由$\boldsymbol u$得到的Teoplitz矩阵。$x$则是一个待优化的变量。 这个约束隐含了如下的结论：

1. $\boldsymbol T(\boldsymbol u) \succeq 0$，否则无法保证对于任何向量$\boldsymbol y$都有$\boldsymbol y^{\mathrm H} \left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right]\boldsymbol y \geq 0$；
2. $\boldsymbol z$ 一定位于$\boldsymbol T(\boldsymbol u)$的列空间中，即$\boldsymbol z$能由$r=\mathrm{rank}(\boldsymbol T(\boldsymbol u))$个原子线性表示。

由于$\boldsymbol T(\boldsymbol u) \succeq 0$， 因此$\boldsymbol T(\boldsymbol u)$是一个半正定矩阵， 即存在范德蒙德分解：

$$\boldsymbol T(\boldsymbol u) = \sum_{k=1}^r p_k \boldsymbol a(f_k) \boldsymbol a^{\mathrm H}(f_k) = \boldsymbol A(\boldsymbol f) \mathrm{diag}(\boldsymbol p) \boldsymbol A^{\mathrm H}(\boldsymbol f) \tag{5-76}$$

其中，$r = \mathrm{rank}(\boldsymbol T(\boldsymbol u))$。由于$\boldsymbol z$一定位于$\boldsymbol T(\boldsymbol u)$的列空间中，那么$\boldsymbol z$必能写为$\boldsymbol a(f_k)$的线性组合！ 这一点至关重要，因为这引出了如下的结论：

最小化原子$l_0$范数等价于求解如下问题：

$$\begin{cases} \min\limits_{x, \boldsymbol u} & r = \mathrm{rank}(\boldsymbol T(\boldsymbol u)) \\ \text { s.t. } & \left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right] \succeq 0 \end{cases} \tag{5-77}$$

由于$\boldsymbol z$是$\boldsymbol T(\boldsymbol u)$范德蒙德分解所得的$r$个原子的线性组合。当找到秩最小的$r$时，也就找到了组成$\boldsymbol z$所需的最少原子数，也就对应$\boldsymbol z$的原子$l_0$范数最小化。 而对于这个转化后的问题，如果$\boldsymbol u$被解出，那么$\boldsymbol T(\boldsymbol u)$也能得到，那么对$\boldsymbol T(\boldsymbol u)$进行范德蒙德分解，也就获得了$\boldsymbol A(\boldsymbol f)$。但是美中不足的是，该目标函数并非凸函数，无法轻易求解。

将原子$l_0$范数进行凸松弛， 得到的就是原子$l_1$范数，简称原子范数。 其定义如下：

对于原子集合$\mathcal A$，若其凸包$\mathrm{conv}(\mathcal A)$相对于原点是一个中心对称的紧集，且包含原点作为内点，这意味着$\mathcal A$中的任一元素$\boldsymbol a \in \mathcal A$不会位于除$\boldsymbol a$以外的其他元素所构成的凸包$ \mathrm{conv}(\mathcal A \backslash \boldsymbol a)$内，即$\mathcal A$中的元素都是$\mathrm{conv}(\mathcal A)$的极值点，$\boldsymbol a \in \mathcal A$当且仅当$-\boldsymbol a \in \mathcal A$。此时由凸包$\mathrm{conv}(\mathcal A)$的尺度函数定义的范数称为原子范数，用$||\boldsymbol z||_{\mathcal A, 1}$或$||\boldsymbol z||_{\mathcal A}$表示(详细内容参考：如何解决基不匹配问题：从原子范数到无网格压缩感知或者文献[61])：

$$\begin{aligned} ||\boldsymbol z||_{\mathcal A, 1} =& \inf \left\{ t>0 \left| \boldsymbol z \in t\mathrm{conv}({\mathcal A}) \right.\right\} \\ = & \inf_{f_k, s_k} \left\{ {\sum_k ∣s_k∣ \left| \boldsymbol z = \sum_k \boldsymbol a(f_k) s_k , \quad f_k \in \mathbb T \right.} \right\} \end{aligned} \tag{5-78}$$

与原子$l_0$范数的定义进行比较，发现这和将传统的$\ell_0$-范数松弛为$\ell_1$-范数如出一辙。然而如何最小化$||\boldsymbol z||_{\mathcal A, 1}$看上去也非常困难。此时再度利用范德蒙德分解，有如下精彩的结论。

最小化原子范数(Atomic Norm Minimization, ANM)等价于求解如下问题：

$$\begin{cases} \min\limits_{x, \boldsymbol u} & \dfrac{1}{2}x + \dfrac{1}{2}u_1 \\ \text { s.t. } & \left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right] \succeq 0 \end{cases} \tag{5-79}$$

注意到，这是一个凸问题。首先目标函数显然是变量的仿射函数，因此为凸(既凸且凹)。限制条件也可以写为变量的仿射函数形式。 因此也满足凸问题的限制条件($f(x) \leq 0$，$f(x)$为凸函数)。($\boldsymbol X \succeq 0$可以等价为$\boldsymbol y^{\mathrm H} \boldsymbol X \boldsymbol y \geq 0, \quad \forall \boldsymbol y$，而对于每个$\boldsymbol y$， 都是关于$\boldsymbol X$的仿射变换）

$\boldsymbol T(\boldsymbol u)$的对角元素都是$u_1$， 那么事实上$u_1$就是$\mathrm{tr}(\boldsymbol T(u))$！而后者又被称为核函数， 也是$\mathrm{rank}(\boldsymbol T(u))$的经典凸松弛。 这从另一个角度解释了原子范数最小化是原子$l_0$范数最小化的凸松弛。

至此，压缩感知问题被转化为了一个可以由CVX进行直接求解的凸问题！还剩的最后一块拼图：即为何$\ell_1$-原子范数可以等效为这个凸问题(Ps：其实我已经推导过了，详细请参考相关博客)。

原子范数与原子$l_0$范数类似，频率也被编码在Toeplitz矩阵$\boldsymbol T(\boldsymbol u)$中，因此原子范数可以视为基于协方差矩阵的，但原子范数与原子$l_0$范数有一个区别在于强制$\boldsymbol T(\boldsymbol u)$的低秩性的方法：

相关文献研究了无噪声时，ANM的理论性能。在ULA情况下由$\min_{\boldsymbol z} \mathcal{M}(\boldsymbol z) \quad \text{s.t.} \quad \boldsymbol z = \boldsymbol y$导出的ANM问题实际允许平凡解$\boldsymbol{z=y}$。但是从式$(5-79)$导出的如下SDP问题的解仍有意义，可用于频率估计：

$$\begin{cases} \min\limits_{x, \boldsymbol u} & \dfrac{1}{2}x + \dfrac{1}{2}u_1 \\ \text { s.t. } & \left[\begin{array}{cc} x & \boldsymbol y^{\mathrm H} \\ \boldsymbol y & \boldsymbol T(\boldsymbol u) \end{array}\right] \succeq 0 \end{cases} \tag{5-80}$$

令$\mathcal{T}=\left\{f_1, \ldots, f_K\right\}$，定义$\mathcal T$的最小间隔为任意两个元素之间最近的环绕距离：

$$\Delta \mathcal{T} = \inf_{1 \leq j \neq k \leq K} \min\{|f_j - f_k|, 1-|f_j, f_k|\} \tag{5-81}$$

当考虑噪声时，由式$(5-72)$可得相应的SDP问题为：

$$\begin{cases} \min\limits_{x, \boldsymbol u, \boldsymbol z} & \dfrac{\lambda}{2}(x + \mu_1) + \dfrac{1}{2}||\boldsymbol z - \boldsymbol y||_2^2 \\ \text { s.t. } & \left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right] \succeq 0 \end{cases} \tag{5-82}$$

其中，正则化参数$\lambda$明显是用来平衡信号稀疏性和保真度的，但如何选择正则化参数却不清晰。在i.i.d高斯噪声假设下，相关文献进行了研究。

(2) 基于Hankel的核范数

另一种$\mathcal{M}(\boldsymbol z)$的选择是基于Hankel的核范数，这种$\mathcal{M}(\boldsymbol z)$度量是基于以下观察提出的：给定式$(5-70)$中的$\boldsymbol z$，构造如下Hankel矩阵：

$$\boldsymbol H(\boldsymbol z) = \left[\begin{array}{cccc} z_1 & z_2 & \cdots & z_n \\ z_2 & z_3 & \cdots & z_{n+1} \\ \vdots & \vdots & \ddots & \vdots \\ z_m & z_{m+1} & \cdots & z_N \end{array}\right] \tag{5-83}$$

其中，$m+n = N+1$，此时也就意味着：

$$\boldsymbol H(\boldsymbol z) = \sum_{k=1}^{K} s_k \left[\begin{array}{c} 1 \\ \mathrm{e}^{\mathrm i 2\pi f_k} \\ \vdots \\ \mathrm{e}^{\mathrm i 2\pi (m-1)f_k} \end{array}\right] \left[\begin{array}{c} 1 & \mathrm{e}^{\mathrm i 2\pi f_k} & \cdots & \mathrm{e}^{\mathrm i 2\pi (n-1)f_k} \end{array}\right] \tag{5-84}$$

如果$K < \min(m,n)$，则可以得到$\boldsymbol H(\boldsymbol z)$是一个低秩矩阵：$\mathrm{rank}(\boldsymbol H(\boldsymbol z)) = K$。

因此为了重构$\boldsymbol z$，可以选择$\mathrm{rank}(\boldsymbol H(\boldsymbol z))$为稀疏度量，如果可以确定$\boldsymbol z$，那么可以从$\boldsymbol z$中恢复频率。

由于$\mathrm{rank}(\boldsymbol H(\boldsymbol z))$很难解决，寻找$\mathrm{rank}(\boldsymbol H(\boldsymbol z))$的凸松弛，而核范数是一个$\mathrm{rank}$常用的凸松弛，此时有：

$$\mathcal{M}(\boldsymbol z) = ||\boldsymbol H(\boldsymbol z)||_\star \tag{5-85}$$

其中，核范数定义式为$||\boldsymbol A||_{\star} = \mathrm{Tr}\sqrt{\boldsymbol A^{\mathrm H}\boldsymbol A}$。将式$(5-71)$和式$(5-72)$中的$\mathcal{M}(\boldsymbol z)$替换为核范数，所导出的优化问题为增强矩阵补全(EMaC)问题。

来看一下核范数本身所导出的SDP问题：

$$\begin{cases} ||\boldsymbol H(\boldsymbol z)||_{\star} = &\min\limits_{\boldsymbol Q_1, \boldsymbol Q_2} & \dfrac{1}{2}\left[\mathrm{Tr}(\boldsymbol Q_1) + \mathrm{Tr}(\boldsymbol Q_1) \right] \\ &\text { s.t. } & \left[\begin{array}{cc} \boldsymbol Q_1 & \boldsymbol H^{\mathrm H}(\boldsymbol z) \\ \boldsymbol H(\boldsymbol z) & \boldsymbol Q_2 \end{array}\right] \succeq 0 \end{cases} \tag{5-86}$$

相关的证明以及矩阵补全问题请参考如下几个链接：
矩阵核范数优化和SDP等价性 - 吃饭吧唧吧唧嘴的文章 - 知乎
涉及矩阵范数的优化问题总结 - 博客园
矩阵补全中核范数优化问题如何求解？ - Orbitopes的回答 - 知乎
矩阵补全（matrix completion）的经典算法有哪些？目前比较流行的算法是什么？ - 知乎
研究领域总结（二）：稀疏——矩阵补全 - 博客园

与原子范数一样，EMaC问题也可以作为SDP，使用现有的优化器求解。

(3) ANM与EMaC的联系与区别

。。。。。。

(四) 协方差拟合：无网格SPICE(GLS)

GLS被引入作为SPICE算法的无网格版本的改进。由于SPICE是基于协方差的，并且数据协方差是感兴趣的DOA参数的高度非线性函数，因此在SPICE中执行网格化以基于零阶近似将问题线性化。GLS的关键思想是通过使用Toeplitz协方差矩阵的Vandermonde分解，使用在新参数向量$\boldsymbol u$中是线性的PSD Toeplitz矩阵$\boldsymbol T(\boldsymbol u)$来重新参数化数据协方差矩阵。为了推导GLS，我们自然会做出与SPICE相同的假设。

首先考虑ULA的情况。假设噪声是同方差的(注意，与SPICE一样，GLS可以扩展到异方差噪声的情况)。数据协方差矩阵$\boldsymbol R$是Toeplitz矩阵。因此，$\boldsymbol R$可以线性地重新参数化为：

$$\boldsymbol R = \boldsymbol T(\boldsymbol u), \quad \boldsymbol T(\boldsymbol u) \succeq 0 \tag{5-87}$$

在单快拍时，SPICE将最小化以下协方差拟合准则：

$$\left\|\boldsymbol{R}^{-\frac{1}{2}}\left(\boldsymbol{y} \boldsymbol{y}^{\mathrm{H}}-\boldsymbol{R}\right)\right\|_{\mathrm{F}}^2=\|\boldsymbol{y}\|_2^2 \cdot \boldsymbol{y}^{\mathrm{H}} \boldsymbol{R}^{-1} \boldsymbol{y}+\operatorname{Tr}(\boldsymbol{R})-2\|\boldsymbol{y}\|_2^2 \tag{5-88}$$

将式$(5-87)$带入式$(5-88)$，得到GLS的优化问题：

$$\begin{aligned} &\begin{cases} \min\limits_{\boldsymbol u} & \|\boldsymbol{y}\|_2^2 \cdot \boldsymbol{y}^{\mathrm{H}} \boldsymbol{T}^{-1}(\boldsymbol u) \boldsymbol{y}+\operatorname{Tr}(\boldsymbol{T}(\boldsymbol u)) \\ \text { s.t. } &\boldsymbol T(\boldsymbol u) \succeq 0 \end{cases} \\ & \qquad \qquad \qquad\Updownarrow \\ &\begin{cases} \min\limits_{x, \boldsymbol u} & \|\boldsymbol{y}\|_2^2 x+Mu_1 \\ \text { s.t. } &\boldsymbol T(\boldsymbol u) \succeq 0 \\ & x \geq \boldsymbol{y}^{\mathrm{H}} \boldsymbol{T}^{-1}(\boldsymbol u) \boldsymbol{y} \end{cases} \\ & \qquad \qquad \qquad\Updownarrow \\ &\begin{cases} \min\limits_{x, \boldsymbol u} & \|\boldsymbol{y}\|_2^2 x+Mu_1 \\ \text { s.t. } & \left[\begin{array}{cc} x & \boldsymbol y^{\mathrm H} \\ \boldsymbol y & \boldsymbol T( \boldsymbol u) \end{array}\right] \succeq 0 \end{cases} \end{aligned} \tag{5-89}$$

因此协方差拟合问题被视为可以在多项式时间内求解的SDP问题，问题求解后，就能得到数据协方差估计值$\hat{\boldsymbol R} = \boldsymbol T(\hat{\boldsymbol u})$，再利用推论(也即式$(5-69)$)，从$\hat{\boldsymbol R}$中估计出参数$(\hat{\boldsymbol f}, \hat{\boldsymbol p}, \sigma)$。

GLS保证产生一个具有最多$N-1$个源的稀疏解。这是频率检索步骤的直接结果(同样见推论)，一般来讲，在存在噪声时，GLS会高估真实源的数目$K$，这是合理的，因为GLS不假设对源数或者噪声方差有任何了解。

相关文献中根据数据协方差估计$\hat{\boldsymbol R}$的特征值提出了一种自动源数估计法(又称模型阶次选择)，其基本思想是——较大的特征值对应源，较小的特征值对应噪声。一种称为SORTE算法可以用来区分两组特征值。

(五) ANM与GLS的关联与区别

。。。。。。

参考文献

[1] Brennan L E, Reed L S. Theory of adaptive radar[J]. IEEE transactions on Aerospace and Electronic Systems, 1973 (2): 237-252.

[2] Reed I S, Mallett J D, Brennan L E. Rapid convergence rate in adaptive arrays[J]. IEEE Transactions on Aerospace and Electronic Systems, 1974 (6): 853-863.

[3] Melvin W L, Showman G A. An approach to knowledge-aided covariance estimation[J]. IEEE Transactions on Aerospace and Electronic Systems, 2006, 42(3): 1021-1042.

[4] Melvin W L, Guerci J R. Adaptive detection in dense target environments[C]//Proceedings of the 2001 IEEE Radar Conference (Cat. No. 01CH37200). IEEE, 2001: 187-192.

[5] Guerci J R, Baranoski E J. Knowledge-aided adaptive radar at DARPA: An overview[J]. IEEE Signal Processing Magazine, 2006, 23(1): 41-50.

[6] Brown R D, Wicks M C, Zhang Y, et al. A space-time adaptive processing approach for improved performance and affordability[C]//Proceedings of the 1996 IEEE National Radar Conference. IEEE, 1996: 321-326.

[7] Goldstein J S, Reed I S. Reduced-rank adaptive filtering[J]. IEEE Transactions on Signal Processing, 1997, 45(2): 492-496.

[8] Pillai S U, Lim Y L, Guerci J R. Generalized forward/backward subaperture smoothing techniques for sample starved STAP[J]. IEEE Transactions on Signal Processing, 2000, 48(12): 3569-3574.

[9] Sarkar T K, Wang H, Park S, et al. A deterministic least-squares approach to space-time adaptive processing (STAP)[J]. IEEE Transactions on Antennas and Propagation, 2001, 49(1): 91-103.

[10] Parker P, Swindlehurst A. Space-time autoregressive filtering for matched subspace STAP[J]. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39(2): 510-520.

[11] Melvin W L, Showman G A. An approach to knowledge-aided covariance estimation[J]. IEEE Transactions on Aerospace and Electronic Systems, 2006, 42(3): 1021-1042.

[12] 冯建婷. 宽带机载相控阵雷达空时自适应处理方法研究[D]. 西安电子科技大学, 2021.

[13] 孟庆统. 机载DDMA MIMO雷达杂波建模与杂波抑制[D].西安电子科技大学,2021.

[14] Ward J .Space-time adaptive processing for airborne radar[C]//Space-Time Adaptive Processing (Ref. No. 1998/241), IEE Colloquium on.1998.

[15] 秦兆锐. 对空时二维处理的干扰技术研究[D]. 陕西:西安电子科技大学,2020.

[16] 闵丛丛. 基于学习的空时自适应处理方法研究[D]. 四川:电子科技大学,2017.

[17] Yang Z, Li X, Wang H, et al. On clutter sparsity analysis in space–time adaptive processing airborne radar[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(5): 1214-1218.

[18] Maria S, Fuchs J J. Application of the global matched filter to STAP data an efficient algorithmic approach[C]//2006 IEEE International Conference on Acoustics Speech and Signal Processing Proceedings. IEEE, 2006, 4: IV-IV.

[19] Sun K, Zhang H, Li G, et al. A novel STAP algorithm in heterogeneous [A]. IEEE International Conference on Geoscience ＆ Remote Sensing Symposium [C]. Cape Town: IEEE, 2009. 336－339.

[20] Li J, Zhu X, Stoica P, et al. High resolution angle-Doppler imaging for MTI radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2010, 46(3): 1544-1556.

[21] Selesnick I W, Pillai S U, Li K Y, et al. Angle-Doppler processing using sparse regularization[C]//2010 IEEE International Conference on Acoustics, Speech and Signal Processing. IEEE, 2010: 2750-2753.

[22] Parker J T, Potter L C. A Bayesian perspective on sparse regularization for STAP post-processing[C]//2010 IEEE Radar Conference. IEEE, 2010: 1471-1475.

[23] Sun K, Meng H, Wang Y, et al. Direct data domain STAP using sparse representation of clutter spectrum[J]. Signal Processing, 2011, 91(9): 2222-2236.

[24] Yang Z, de Lamare R C, Li X. $ L_1 $-regularized STAP algorithms with a generalized sidelobe canceler architecture for airborne radar[J]. IEEE Transactions on Signal Processing, 2011, 60(2): 674-686.

[25] Ma Z, Liu Y, Meng H, et al. Jointly sparse recovery of multiple snapshots in STAP[C]//2013 IEEE Radar Conference (RadarCon13). IEEE, 2013: 1-4.

[26] Sen S. OFDM radar-space-time adaptive processing by exploiting spatio-temporal sparsity [J]. IEEE Transactions on Signal Processing, 2013, 61(1): 118 - 130.

[27] Yang Z, Li X, Wang H, et al. On clutter sparsity analysis in space–time adaptive processing airborne radar[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(5): 1214-1218.

[28] Ma Z, Liu Y, Meng H, et al. Sparse recovery‐based space‐time adaptive processing with array error self‐calibration[J]. Electronics letters, 2014, 50(13): 952-954.

[29] Feng W, Zhang Y, He X, et al. Cascaded clutter and jamming suppression method using sparse representation[J]. Electronics Letters, 2015, 51(19): 1524-1526.

[30] Wang Z, Wang Y, Duan K, et al. Subspace-augmented clutter suppression technique for STAP radar[J]. IEEE Geoscience and Remote Sensing Letters, 2016, 13(3): 462-466.

[31] Yang X, Sun Y, Zeng T, et al. Fast STAP method based on PAST with sparse constraint for airborne phased array radar[J]. IEEE Transactions on Signal Processing, 2016, 64(17): 4550-4561.

[32] Yang Z, Li X, Wang H, et al. Knowledge‐aided STAP with sparse‐recovery by exploiting spatio‐temporal sparsity[J]. IET Signal Processing, 2016, 10(2): 150-161.

[33] Guo Y, Liao G, Feng W. Sparse representation based algorithm for airborne radar in beam-space post-Doppler reduced-dimension space-time adaptive processing[J]. IEEE Access, 2017, 5: 5896-5903.

[34] Duan K, Wang Z, Xie W, et al. Sparsity‐based STAP algorithm with multiple measurement vectors via sparse Bayesian learning strategy for airborne radar[J]. IET Signal Processing, 2017, 11(5): 544-553.

[35] Wang Z, Xie W, Duan K, et al. Clutter suppression algorithm based on fast converging sparse Bayesian learning for airborne radar[J]. Signal Processing, 2017, 130: 159-168.

[36] Sun Y, Yang X, Long T, et al. Robust sparse Bayesian learning STAP method for discrete interference suppression in nonhomogeneous clutter[C]//2017 IEEE Radar Conference (RadarConf). IEEE, 2017: 1003-1008.

[37] Zetao W, Yongliang W, Fei G, et al. Clutter nulling space‐time adaptive processing algorithm based on sparse representation for airborne radar[J]. IET Radar, Sonar & Navigation, 2017, 11(1): 177-184.

[38] Feng W, Guo Y, Zhang Y, et al. Airborne radar space time adaptive processing based on atomic norm minimization[J]. Signal Processing, 2018, 148: 31-40.

[39] Duan K, Liu W, Duan G, et al. Off‐grid effects mitigation exploiting knowledge of the clutter ridge for sparse recovery STAP[J]. IET Radar, Sonar & Navigation, 2018, 12(5): 557-564.

[40] Cotter S F, Rao B D, Engan K, et al. Sparse solutions to linear inverse problems with multiple measurement vectors[J]. IEEE Transactions on signal processing, 2005, 53(7): 2477-2488.

[41] Tropp J A, Gilbert A C, Strauss M J. Algorithms for simultaneous sparse approximation. Part I: Greedy pursuit[J]. Signal processing, 2006, 86(3): 572-588.

[42] Tropp J A. Algorithms for simultaneous sparse approximation. Part II: Convex relaxation[J]. Signal Processing, 2006, 86(3): 589-602.

[43] Chen J, Huo X. Theoretical results on sparse representations of multiple-measurement vectors[J]. IEEE Transactions on Signal processing, 2006, 54(12): 4634-4643.

[44] Wipf D P, Rao B D. An empirical Bayesian strategy for solving the simultaneous sparse approximation problem[J]. IEEE Transactions on Signal Processing, 2007, 55(7): 3704-3716.

[45] Tsao T, Himed B, Michels J H. Effects of interference rank estimation on the detection performance of rank reduced STAP algorithms[C]//Proceedings of the 1998 IEEE Radar Conference, RADARCON’98. Challenges in Radar Systems and Solutions (Cat. No. 98CH36197). IEEE, 1998: 147-152.

[46] Goodman N A, Stiles J M. On clutter rank observed by arbitrary arrays[J]. IEEE Transactions on Signal processing, 2006, 55(1): 178-186.

[47] Tipping M E. Sparse Bayesian learning and the relevance vector machine[J]. Journal of machine learning research, 2001, 1(Jun): 211-244.

[48] Wipf D P, Rao B D. Sparse Bayesian learning for basis selection[J]. IEEE Transactions on Signal processing, 2004, 52(8): 2153-2164.

[49] Wang Y L, Peng Y N, Bao Z. Space–time adaptive processing for airborne radar with various array orientations[J]. IEE Proceedings-Radar, Sonar and Navigation, 1997, 144(6): 330-340.

[50] Meng X, Wang T, Wu J, et al. Short-range clutter suppression for airborne radar by utilizing prefiltering in elevation[J]. IEEE Geoscience and Remote Sensing Letters, 2009, 6(2): 268-272.

[51] 段克清, 谢文冲, 陈辉, 等. 基于俯仰维信息的机载雷达非均匀杂波抑制方法[J]. 电子学报, 2011, 39 (3): 585－590.

[52] Wu J, Wang T, Meng X, et al. Clutter suppression for airborne non-sidelooking radar using ERCB-STAP algorithm[J]. IET radar, sonar & navigation, 2010, 4(4): 497-506.

[53] Yang Z, Xie L. On gridless sparse methods for line spectral estimation from complete and incomplete data[J]. IEEE Transactions on Signal Processing, 2015, 63(12): 3139-3153.

[54] Candès E J, Fernandez‐Granda C. Towards a mathematical theory of super‐resolution[J]. Communications on pure and applied Mathematics, 2014, 67(6): 906-956.

[55] Tang G, Bhaskar B N, Shah P, et al. Compressed sensing off the grid[J]. IEEE transactions on information theory, 2013, 59(11): 7465-7490.

[56] Sun K, Liu Y, Meng H, et al. Adaptive sparse representation for source localization with gain/phase errors[J]. Sensors, 2011, 11: 4780-4793.

[57] Liu A, Sun H, Teh K C, et al. Robust space-time adaptive processing for nonhomogeneous clutter in the presence of model errors[J]. IEEE Transactions on Aerospace and Electronic Systems, 2016, 52(1): 155-168.

[58] Zhaocheng Y, de Lamare R C, Liu W. Sparsity-based STAP using alternating direction method with gain/phase errors[J]. IEEE Transactions on Aerospace and Electronic Systems, 2017, 53(6): 2756-2768.

[59] Yang Z, Li J, Stoica P, et al. Sparse methods for direction-of-arrival estimation[M]//Academic Press Library in Signal Processing, Volume 7. Academic Press, 2018: 509-581.

[60] C. Austin, J. Ash, R. Moses, Dynamic dictionary algorithms for model order and parameter estimation, IEEE Transactions on Signal Processing 61 (20) (2013) 5117–5130.

[61] 李慧启. 基于无网格压缩感知的跳频信号参数估计研究[D].战略支援部队信息工程大学,2019.

[62] Ottersten B, Stoica P, Roy R. Covariance matching estimation techniques for array signal processing applications[J]. Digital Signal Processing, 1998, 8(3): 185-210.