Compressive Two-Dimensional Harmonic Retrieval via Atomic Norm Minimization

摘要—— 本文关注从部分时间样本中估计二维 (2-D) 频率的问题，这一问题出现在许多应用中，如雷达、逆散射和超分辨率成像。假设研究对象是r rr个连续值二维正弦波的混合。

目标是在我们仅掌握n nn个等间距时间样本的随机子集信息时，识别出所有的频率分量。我们证明，在一定的温和谱分离条件下（under some mild spectral separation condition），只要样本复杂度超过r log ⁡ r log ⁡ n r \log r \log nrlogrlogn的量级，就有可能通过求解原子范数最小化规划来精确恢复所有频率。随后，我们提出通过半定规划来求解原子范数最小化，并提供数值算例以验证其实际能力。我们的工作将 Tang 等人针对线谱估计提出的框架扩展到了二维频率模型。

索引词—— 原子范数，基不匹配，连续值频率恢复，稀疏性。

II. PROBLEM FORMULATION AND RELATED WORK

A. Problem Formulation

不失一般性，考虑一个尺寸为n = ( 4 M + 1 ) × ( 4 M + 1 ) n = (4M + 1) \times (4M + 1)n=(4M+1)×(4M+1)的二维方形数据矩阵X ∗ \boldsymbol{X}^*X∗，其中J = { − 2 M , … , 0 , … , 2 M } × { − 2 M , … , 0 , … , 2 M } J = \{-2M, \dots, 0, \dots, 2M\} \times \{-2M, \dots, 0, \dots, 2M\}J={−2M,…,0,…,2M}×{−2M,…,0,…,2M}表示X ∗ \boldsymbol{X}^*X∗的索引的并集。施加此假设是为了简化理论保证的推导，只需稍作修改即可移除，类似处理方法见 [29]。X ∗ \boldsymbol{X}^*X∗的每个条目都可以表示为在时间索引k = [ k 1 , k 2 ] ∈ J \boldsymbol{k} = [k_1, k_2] \in Jk=[k1​,k2​]∈J处观测到的r rr个复正弦波的叠加，即

x k ∗ = x k 1 , k 2 ∗ = 1 ( 4 M + 1 ) ∑ i = 1 r d i e j 2 π f i ⊤ k , (1) x_{\boldsymbol{k}}^* = x_{k_1, k_2}^* = \frac{1}{(4M+1)} \sum_{i=1}^{r} d_i e^{j 2\pi \boldsymbol{f}_i^\top \boldsymbol{k}}, \tag{1}xk∗​=xk1​,k2​∗​=(4M+1)1​i=1∑r​di​ej2πfi⊤​k,(1)

其中d i d_idi​表示与每个1 ≤ i ≤ r 1 \le i \le r1≤i≤r相关的复振幅。令Ω = { f i = ( f 1 i , f 2 i ) ∈ [ 0 , 1 ) × [ 0 , 1 ) , 1 ≤ i ≤ r } \Omega = \{\boldsymbol{f}_i = (f_{1i}, f_{2i}) \in [0, 1) \times [0, 1), 1 \le i \le r\}Ω={fi​=(f1i​,f2i​)∈[0,1)×[0,1),1≤i≤r}为不同频率的集合。为了符号简洁，我们引入以下单位范数原子：

{ a ( f 1 i ) = 1 ( 4 M + 1 ) [ y i − 2 M , … , 1 , … , y i 2 M ] ⊤ , a ( f 2 i ) = 1 ( 4 M + 1 ) [ z i − 2 M , … , 1 , … , z i 2 M ] ⊤ , \begin{cases} \boldsymbol{a}(f_{1i}) = \frac{1}{\sqrt{(4M+1)}} \left[ y_i^{-2M}, \dots, 1, \dots, y_i^{2M} \right]^\top, \\ \boldsymbol{a}(f_{2i}) = \frac{1}{\sqrt{(4M+1)}} \left[ z_i^{-2M}, \dots, 1, \dots, z_i^{2M} \right]^\top, \end{cases}⎩⎨⎧​a(f1i​)=(4M+1)​1​[yi−2M​,…,1,…,yi2M​]⊤,a(f2i​)=(4M+1)​1​[zi−2M​,…,1,…,zi2M​]⊤,​

其中y i = e j 2 π f 1 i y_i = e^{j 2\pi f_{1i}}yi​=ej2πf1i​，且z i = e j 2 π f 2 i z_i = e^{j 2\pi f_{2i}}zi​=ej2πf2i​。这使得我们可以将X ∗ \boldsymbol{X}^*X∗写成如下矩阵形式

X ∗ = Y D Z ⊤ , (2) \boldsymbol{X}^* = \boldsymbol{Y} \boldsymbol{D} \boldsymbol{Z}^\top, \tag{2}X∗=YDZ⊤,(2)

其中Y \boldsymbol{Y}Y由下式给出

Y = [ a ( f 11 ) , … , a ( f 1 r ) ] ∈ C ( 4 M + 1 ) × r , (3) \boldsymbol{Y} = [\boldsymbol{a}(f_{11}), \dots, \boldsymbol{a}(f_{1r})] \in \mathbb{C}^{(4M+1)\times r}, \tag{3}Y=[a(f11​),…,a(f1r​)]∈C(4M+1)×r,(3)

Z = [ a ( f 21 ) , … , a ( f 2 r ) ] ∈ C ( 4 M + 1 ) × r , (4) \boldsymbol{Z} = [\boldsymbol{a}(f_{21}), \dots, \boldsymbol{a}(f_{2r})] \in \mathbb{C}^{(4M+1)\times r}, \tag{4}Z=[a(f21​),…,a(f2r​)]∈C(4M+1)×r,(4)

以及

D = diag ( [ d 1 , d 2 , ⋯ , d r ] ) = diag ( d ) ∈ C r × r . (5) \boldsymbol{D} = \text{diag}([d_1, d_2, \cdots, d_r]) = \text{diag}(\boldsymbol{d}) \in \mathbb{C}^{r \times r}. \tag{5}D=diag([d1​,d2​,⋯,dr​])=diag(d)∈Cr×r.(5)

令（Denote by）x ∗ = vec ( ( X ∗ ) ⊤ ) ∈ C ( 4 M + 1 ) 2 \boldsymbol{x}^* = \text{vec}((\boldsymbol{X}^*)^\top) \in \mathbb{C}^{(4M+1)^2}x∗=vec((X∗)⊤)∈C(4M+1)2表示向量化的数据矩阵，则有（then one has）

x ∗ = ( Y ⊗ Z ) d = ∑ i = 1 r d i a ( f 1 i ) ⊗ a ( f 2 i ) = ∑ i = 1 r d i c ( f i ) , (6) \begin{aligned} \boldsymbol{x}^* &= (\boldsymbol{Y} \otimes \boldsymbol{Z})\boldsymbol{d} = \sum_{i=1}^{r} d_i \boldsymbol{a}(f_{1i}) \otimes \boldsymbol{a}(f_{2i}) \\ &= \sum_{i=1}^{r} d_i \boldsymbol{c}(\boldsymbol{f}_i), \end{aligned} \tag{6}x∗​=(Y⊗Z)d=i=1∑r​di​a(f1i​)⊗a(f2i​)=i=1∑r​di​c(fi​),​(6)

其中⊗ \otimes⊗代表 Kronecker 积，且

c ( f i ) = c ( f 1 i , f 2 i ) : = a ( f 1 i ) ⊗ a ( f 2 i ) ∈ C ( 4 M + 1 ) 2 . \boldsymbol{c}(\boldsymbol{f}_i) = \boldsymbol{c}(f_{1i}, f_{2i}) := \boldsymbol{a}(f_{1i}) \otimes \boldsymbol{a}(f_{2i}) \in \mathbb{C}^{(4M+1)^2}.c(fi​)=c(f1i​,f2i​):=a(f1i​)⊗a(f2i​)∈C(4M+1)2.

满足∥ c ( f i ) ∥ 2 = 1 \|\boldsymbol{c}(\boldsymbol{f}_i)\|_2 = 1∥c(fi​)∥2​=1。

在本文中，我们假设X ∗ \boldsymbol{X}^*X∗的m mm个条目是均匀随机观测到的。具体而言，令T ⊂ J T \subset JT⊂J表示索引集，使得x k 1 , k 2 ∗ x^*_{k_1, k_2}xk1​,k2​∗​被观测到当且仅当( k 1 , k 2 ) ∈ T (k_1, k_2) \in T(k1​,k2​)∈T。定义算子P T \mathcal{P}_TPT​，使得P T ( M ) \mathcal{P}_T(\boldsymbol{M})PT​(M)表示M \boldsymbol{M}M在支撑于T TT上的矩阵子空间上的正交投影。我们将滥用符号（在不引起歧义的情况下），让T TT、J JJ和P T \mathcal{P}_TPT​同时也分别表示观测条目的集合、所有条目的集合以及关于向量化信号x ∗ \boldsymbol{x}^*x∗的观测算子。

本文的主要关注点是恢复原始数据矩阵X ∗ \boldsymbol{X}^*X∗中未被观测到的条目。我们注意到，一旦数据矩阵被恢复，频率Ω \OmegaΩ也可以使用诸如 MEMP 方法 [11] 等传统方法来恢复。

B. Conventional CS Approach

为了应用传统的压缩感知 (CS) 范式，我们将x ∗ \boldsymbol{x}^*x∗表示为预定基下的稀疏信号，方法是通过网格点t = [ t 1 , t 2 ] ∈ Ω d \boldsymbol{t} = [t_1, t_2] \in \Omega_dt=[t1​,t2​]∈Ωd​将二维平面[ 0 , 1 ) × [ 0 , 1 ) [0,1) \times [0,1)[0,1)×[0,1)离散化，其中t 1 , t 2 ∈ { 0 , … , 4 M 4 M + 1 } t_1, t_2 \in \{0, \dots, \frac{4M}{4M+1}\}t1​,t2​∈{0,…,4M+14M​}。将得到的 DFT 基写为

F = [ c ( t ) ∣ t ∈ Ω d ] = F 1 ⊗ F 1 ∈ C ( 4 M + 1 ) 2 × ( 4 M + 1 ) 2 , \boldsymbol{F} = [\boldsymbol{c}(\boldsymbol{t})|_{\boldsymbol{t} \in \Omega_d}] = \boldsymbol{F}_1 \otimes \boldsymbol{F}_1 \in \mathbb{C}^{(4M+1)^2 \times (4M+1)^2},F=[c(t)∣t∈Ωd​​]=F1​⊗F1​∈C(4M+1)2×(4M+1)2,

其中F 1 \boldsymbol{F}_1F1​是一个维数为( 4 M + 1 ) × ( 4 M + 1 ) (4M+1) \times (4M+1)(4M+1)×(4M+1)的 DFT 矩阵。向量化信号x ∗ \boldsymbol{x}^*x∗随后可以使用F \boldsymbol{F}F表示为

x ∗ = F d ~ , (7) \boldsymbol{x}^* = \boldsymbol{F} \tilde{\boldsymbol{d}}, \tag{7}x∗=Fd~,(7)

其中d ~ \tilde{\boldsymbol{d}}d~是近似稀疏的。CS 理论表明，我们可以通过如下ℓ 1 \ell_1ℓ1​最小化来恢复x ∗ \boldsymbol{x}^*x∗

min ⁡ d ~ ∥ d ~ ∥ 1 subject to P T ( F d ~ ) = P T ( x ∗ ) , \min_{\tilde{\boldsymbol{d}}} \|\tilde{\boldsymbol{d}}\|_1 \quad \text{subject to} \quad \mathcal{P}_T(\boldsymbol{F}\tilde{\boldsymbol{d}}) = \mathcal{P}_T(\boldsymbol{x}^*),d~min​∥d~∥1​subject toPT​(Fd~)=PT​(x∗),

其中最小化结果作为d ~ \tilde{\boldsymbol{d}}d~的估计值返回。上述方法的主要问题在于，频率f i \boldsymbol{f}_ifi​绝不会完美地落在网格Ω d \Omega_dΩd​上，从而导致真实频率与离散网格之间不可避免的不匹配问题。文献 [22] 已经证明，稀疏恢复算法的性能可能会显著退化。在本文中，我们将采用一种不同的方法，并尝试在不施加网格的情况下直接恢复频率。

III. ATOMIC NORM MINIMIZATION FOR 2-D HARMONIC RETRIEVAL

原子范数在 [30] 中被提出，作为一种设计用于模型选择的凸优化解的通用方法，其通过对简约模型的原子集合进行凸化来实现。信号模型的原子集合被定义为信号的最简单构建块，例如用于稀疏恢复的单位范数一稀疏向量，用于低秩矩阵补全的单位范数秩一矩阵等等。感兴趣的读者可参阅 [30] 以获取关于原子范数的详细讨论。在二维谐波检索的情况下，可以直接将原子集合定义为所有归一化二维复正弦波的集合：

A : = { c ( f ) ∣ f ∈ [ 0 , 1 ) × [ 0 , 1 ) } , \mathcal{A} := \{\boldsymbol{c}(\boldsymbol{f}) | \boldsymbol{f} \in [0, 1) \times [0, 1) \},A:={c(f)∣f∈[0,1)×[0,1)},

并将信号x \boldsymbol{x}x的原子范数定义为

∥ x ∥ A : = inf ⁡ f i ∈ [ 0 , 1 ) × [ 0 , 1 ) d i ∈ C { ∑ i ∣ d i ∣ ∣ x = ∑ i d i c ( f i ) } . (8) \|\boldsymbol{x}\|_{\mathcal{A}} := \inf_{\substack{f_i \in [0, 1) \times [0, 1) \\ d_i \in \mathbb{C}}} \left\{ \sum_i |d_i| \bigg| \boldsymbol{x} = \sum_i d_i \boldsymbol{c}(\boldsymbol{f}_i) \right\}. \tag{8}∥x∥A​:=fi​∈[0,1)×[0,1)di​∈C​inf​{i∑​∣di​∣​x=i∑​di​c(fi​)}.(8)

这是通过凸化x \boldsymbol{x}x的原子表示获得的，该表示使用最少数量的二维频率尖峰：

∥ x ∥ A , 0 = inf ⁡ f i ∈ [ 0 , 1 ) × [ 0 , 1 ) d i ∈ C { s ∣ x = ∑ i = 1 s d i c ( f i ) } . \|\boldsymbol{x}\|_{\mathcal{A}, 0} = \inf_{\substack{f_i \in [0, 1) \times [0, 1) \\ d_i \in \mathbb{C}}} \left\{ s \bigg| \boldsymbol{x} = \sum_{i=1}^s d_i \boldsymbol{c}(\boldsymbol{f}_i) \right\}.∥x∥A,0​=fi​∈[0,1)×[0,1)di​∈C​inf​{s​x=i=1∑s​di​c(fi​)}.

上述定义推广了 [29] 中一维谐波信号的原子范数，并允许适应更高维度。给定x ∗ \boldsymbol{x}^*x∗的部分观测值（或等价地P T ( x ∗ ) \mathcal{P}_T(\boldsymbol{x}^*)PT​(x∗)），我们尝试通过以下原子范数最小化规划进行恢复

x ^ = arg ⁡ min ⁡ x ∥ x ∥ A subject to P T ( x ) = P T ( x ∗ ) , (9) \hat{\boldsymbol{x}} = \arg\min_{\boldsymbol{x}} \|\boldsymbol{x}\|_{\mathcal{A}} \quad \text{subject to} \quad \mathcal{P}_T(\boldsymbol{x}) = \mathcal{P}_T(\boldsymbol{x}^*), \tag{9}x^=argxmin​∥x∥A​subject toPT​(x)=PT​(x∗),(9)

即寻找满足观测约束且具有最小原子范数的信号。该方法在 [29] 中被用于线谱估计，此时原子集合为A = { a ( f ) ∣ f ∈ [ 0 , 1 ) } \mathcal{A} = \{\boldsymbol{a}(f) | f \in [0, 1)\}A={a(f)∣f∈[0,1)}。在 [29] 中表明，在温和的频率分离条件下，包含O ( r log ⁡ r log ⁡ n ) \mathcal{O}(r \log r \log n)O(rlogrlogn)个样本的随机子集可以确保精确的频率恢复。

定理 1：令M ≥ 256 M \ge 256M≥256。假设我们在大小为∣ T ∣ = m |T| = m∣T∣=m的索引集T ⊂ J T \subset JT⊂J上均匀随机地观测 (1) 中数据矩阵X ∗ \boldsymbol{X}^*X∗的样本，其中f i ∈ [ 0 , 1 ) × [ 0 , 1 ) \boldsymbol{f}_i \in [0, 1) \times [0, 1)fi​∈[0,1)×[0,1)。假设d i d_idi​的符号是独立同分布 (i.i.d.) 的，且均匀抽取自{ + 1 , − 1 } \{+1, -1\}{+1,−1}，并且f i \boldsymbol{f}_ifi​之间的最小间隔满足
Δ min ⁡ ≜ min ⁡ i ≠ j ∥ f i − f j ∥ ∞ = min ⁡ i ≠ j max ⁡ { ∣ f 1 i − f 1 j ∣ , ∣ f 2 i − f 2 j ∣ } ≥ 1.19 M , (10) \Delta_{\min} \triangleq \min_{i \ne j} \|\boldsymbol{f}_i - \boldsymbol{f}_j\|_\infty = \min_{i \ne j} \max \{|f_{1i} - f_{1j}|, |f_{2i} - f_{2j}|\} \ge \frac{1.19}{M}, \tag{10}Δmin​≜i=jmin​∥fi​−fj​∥∞​=i=jmin​max{∣f1i​−f1j​∣,∣f2i​−f2j​∣}≥M1.19​,(10)
其中∣ f 1 i − f 1 j ∣ , ∣ f 2 i − f 2 j ∣ |f_{1i} - f_{1j}|, |f_{2i} - f_{2j}|∣f1i​−f1j​∣,∣f2i​−f2j​∣是单位圆上的环绕距离。那么存在一个数值常数C > 0 C > 0C>0，使得如果
m ≥ C max ⁡ { log ⁡ 2 M δ , r log ⁡ r δ log ⁡ M δ } , (11) m \ge C \max \left\{ \log^2 \frac{M}{\delta}, r \log \frac{r}{\delta} \log \frac{M}{\delta} \right\}, \tag{11}m≥Cmax{log2δM​,rlogδr​logδM​},(11)
则 (9) 的解在概率至少为1 − δ 1 - \delta1−δ下是精确且唯一的。当d i d_idi​的符号是在复单位圆上独立同分布均匀生成时，同样的结论成立，仅 (10) 中的常数不同。

证明可见附录 B。

定理 1 表明，只要频率如 (10) 中那样满足最小分离，一旦m mm达到max ⁡ { log ⁡ 2 n , r log ⁡ r log ⁡ n } \max\{\log^2 n, r \log r \log n\}max{log2n,rlogrlogn}的量级，通过原子范数最小化进行的恢复就是精确的。这一阶数界限与 [29] 中推导出的线谱估计性能保证一致。

我们将定理 1 与传统的子空间方法（如 ESPRIT）进行了比较。

• ESPRIT 能够从数据矩阵X ∗ \boldsymbol{X}^*X∗的Θ ( r ) \Theta(r)Θ(r)个连续样本中恢复潜在频率。精确恢复所需的样本数量仅取决于潜在的自由度，而与X ∗ \boldsymbol{X}^*X∗的环境维数无关。相比之下，所提出的算法 (9) 假设对数据矩阵X ∗ \boldsymbol{X}^*X∗进行随机子采样，并且需要略高的样本复杂度，大约为r log ⁡ r log ⁡ n r \log r \log nrlogrlogn的量级。
• 此外，在没有噪声的情况下，ESPRIT 允许在不施加如 (10) 所示的分离条件的情况下进行恢复。然而请注意，如 [28], [37] 中详述，当存在噪声时，分离条件是必须的。我们将通过数值算例证明，所提出的算法 (9) 在噪声观测下也是稳定的。

我们还将定理 1 与 CS [14] 中的标准结果进行了比较。当Ω \OmegaΩ中的频率确实位于 DFT 网格上时，CS 允许从数量为O ( r log ⁡ ( n / r ) ) \mathcal{O}(r \log(n/r))O(rlog(n/r))的样本中恢复r rr个复正弦波。所提出的算法 (9) 可以被视为针对偏离网格目标的 CS 方法的一种补救措施，其样本复杂度略大。

IV. APPROXIMATE SEMIDEFINITE PROGRAM TO SOLVE ATOMIC NORM MINIMIZATION

定理 1 表明，求解原子范数最小化问题 (9) 允许仅从少量时间样本中完美恢复数据矩阵。然而，一个自然的问题是如何以易于处理的方式求解 (9)。

• 不幸的是，文献 [29] 中提出的线谱情况下的原子范数最小化的精确半定规划表征，无法直接扩展到二维模型，这是由于将经典的 Caratheodory 定理 [34] 推广到更高维度存在根本性的困难。

• 尽管如此，在本节中，我们提出了一种半定规划来近似求解 (9)，其在第五节中表现出优异的经验性能。我们还提供了一个充分条件，在该条件下所提出的半定规划会返回 (9) 的解。

我们描述了当X ∗ \boldsymbol{X}^*X∗的维数为n 1 × n 2 n_1 \times n_2n1​×n2​（不一定是方阵）时的一般情况下的算法。我们仍然假设X ∗ \boldsymbol{X}^*X∗满足 (2)，但稍微滥用符号，令

Y = [ a ( f 11 ) , … , a ( f 1 r ) ] \boldsymbol{Y} = [\boldsymbol{a}(f_{11}), \dots, \boldsymbol{a}(f_{1r})]Y=[a(f11​),…,a(f1r​)]

其中a ( f 1 i ) = [ 1 , e j 2 π f 1 i , … , e j 2 π ( n 1 − 1 ) f 1 i ] ⊤ \boldsymbol{a}(f_{1i}) = [1, e^{j2\pi f_{1i}}, \dots, e^{j2\pi (n_1-1)f_{1i}}]^\topa(f1i​)=[1,ej2πf1i​,…,ej2π(n1​−1)f1i​]⊤，以及

Z = [ b ( f 21 ) , … , b ( f 2 r ) ] \boldsymbol{Z} = [\boldsymbol{b}(f_{21}), \dots, \boldsymbol{b}(f_{2r})]Z=[b(f21​),…,b(f2r​)]

其中b ( f 2 i ) = [ 1 , e j 2 π f 2 i , … , e j 2 π ( n 2 − 1 ) f 1 i ] ⊤ \boldsymbol{b}(f_{2i}) = [1, e^{j2\pi f_{2i}}, \dots, e^{j2\pi (n_2-1)f_{1i}}]^\topb(f2i​)=[1,ej2πf2i​,…,ej2π(n2​−1)f1i​]⊤，只要在上下文中清晰即可。