压缩感知(compressed sensing, CS)在过去的十几年中，由于其在图像和音频处理^[1]、视觉跟踪^[2-4]、以及计算机视觉和水声信道估计^[5-8]等多个领域的巨大潜在应用价值，受到了越来越多的关注。

结合测量矩阵和压缩的信号对原始信号进行恢复是压缩感知最重要的部分之一，其可以对采样不足的信号进行充分的重构，突破了传统的香农奈奎斯特采样定理，实现了可接受的准确性和复杂度。常用的重构算法主要包括凸优化算法和贪婪迭代算法^[9]。其中，贪婪算法以其兼顾重构质量和运行效率的优势，受到越来越多的关注。匹配追踪算法(matching pursuit)是最典型的贪婪算法之一，其特点是根据传感矩阵与残差值相关性的强弱，挑选最为匹配的原子扩充支撑集，其重构速度快且有较好的重构质量。其改进算法还包括OMP (orthogonal MP)^[10]、ROMP (regularized OMP)^[11]、StOMP (stagewise OMP)^[12]等算法。然而，挑选的原子被选入支撑集后，将在整个迭代过程中保留，无法判断和更新错选的原子，影响了支撑集的准确性，且相关性判据有待改进。为此，出现的CoSaMP算法(compressive sampling MP)^[13]和SP算法(subspace pursuit)^[14]算法通过回溯步骤剔除相关性较低的原子，获得更精确的支持集，但是重构效果受到支撑集原子选择方式缺乏灵活性的约束。为了满足在未知稀疏度K的情况下获得更好的跟踪和重建性能，SAMP (sparsity adaptive MP)等算法对匹配跟踪算法进行了大量的改进^[15-18]。这类稀疏度自适应算法的核心是利用一个固定步长多次迭代不断逼近，达到事先不需要知道稀疏度就能准确估计真实的稀疏度K的目的。

在满足稀疏性条件的前提下，MMP算法^[19]从多个候选支撑集中选择最优方案，显然多路径搜索比单路径搜索有更好的性能。但由于MMP计算量大，在没有稀疏信息的情况下很难求解大规模问题。因此对MMP算法进行优化具有重要的意义。本文在MMP算法基础上提出了一种改进的快速稀疏度自适应多路径匹配追踪算法(improved multipath matching pursuit algorithm with sparsity adaptive, IMMP)。与传统算法相比，IMMP算法具有明显的重构精度和时间优势。

1 压缩感知与多路径匹配追踪算法 1.1 压缩感知理论

由CS理论可知，将一个K-稀疏或在某一稀疏域能稀疏表示的N维离散数字信号通过某一个测量矩阵进行有效压缩，我们就能够在线性变换下通过得到少量的观测值精确重构原始信号^[10]。信号x中除了只有K个非零值外其余N-K个值均为零，且K远小于N。采样和压缩过程可以描述为：

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} x}} $

(1)

式中：Φ是大小为M×N的测量矩阵，观测值y的长度为M且满足M≪N。压缩感知理论的目标是基于观测值y和测量矩阵Φ将重构问题转换为求具有稀疏约束的l₀范数的最优解问题。其求l₀范数最优解模型为：

$ \hat{\boldsymbol{x}}=\arg \min \limits_{\boldsymbol{x}}\|\boldsymbol{x}\|_{0}, \quad \text { s.t. } \boldsymbol{y}=\boldsymbol{\varPhi} \boldsymbol{x} $

(2)

式中‖x‖₀为x的l₀范数。如果在某一变换域中使信号x满足稀疏条件，称此变换域为x的稀疏域，在此稀疏域中，信号可以稀疏表示为：

$ \boldsymbol{x}=\boldsymbol{\varPsi} \boldsymbol{\theta} $

(3)

式中：Ψ=[ψ₁ ψ₂ … ψ_N]为稀疏基矩阵; θ为在稀疏基矩阵Ψ下的K稀疏信号矩阵。理论表明，求解最小l₀范数是一个NP-hard问题^[20]，而且求解最小l₀范数可以转化为求解最小l₁范数^[11]，即最小l₁范数优化问题：

$ \hat{\boldsymbol{\theta}}=\underset{\boldsymbol{x}}{\arg \min }\|\boldsymbol{\theta}\|_{1}, \quad \text { s. t. } \boldsymbol{y}=\boldsymbol{\varPhi} \boldsymbol{\varPsi} \boldsymbol{\theta} $

(4)

式中$ \hat {\mathit{\pmb{θ}}} $是稀疏信号θ在变换域Ψ下的稀疏估计信号。

1.2 MMP算法

MMP算法通过将搜索最优路径建模成基于树的路径搜索问题，每次迭代依据测量矩阵与残差的相关性来更新一个支撑集原子，通过贪婪迭代搜索多个路径，每个路径经K次迭代后得到一个含K个原子的支撑集，最终达到选择最优路径即最优支撑集的目的。在多路径搜索中，当前第k层分支节点原子通过计算第k-1层分支节点的残差r^k-1与测量矩阵Φ的内积确定：

$ \varGamma=\arg \max \left|\left\langle\boldsymbol{\varPhi}, r^{k-1}\right\rangle\right| $

(5)

当前残差可以由观测向量y和当前重建的信号Φ_Ω^k$ \hat {\boldsymbol{x}}_k $做差得到:

$ r^{k}=\boldsymbol{y}-\boldsymbol{\varPhi}_{\varOmega^{k}} \hat{\boldsymbol{x}}_{k} $

(6)

式中Φ_Ω^k是在索引集Ω^k条件下由当前k个不相关的原子组成的支撑集，而重构信号$ \hat {\boldsymbol{x}} $通过采用最小二乘法求解，可以描述为：

$ \hat{\boldsymbol{x}}=\arg \min \limits_{\boldsymbol{x}}\|\boldsymbol{y}-\boldsymbol{\varPhi} \boldsymbol{x}\|_{2}, \quad \text { s. t. }\|\boldsymbol{x}\|_{0}=K $

(7)

式中: ‖·‖₂表示l₂范数，在最小二乘法求解过程中，重构信号$ \hat {\boldsymbol{x}}_k $由支撑集Φ_Ω^k的伪逆矩阵与观测向量进行矩阵运算得到：

$ \hat{\boldsymbol{x}}_{k}=\boldsymbol{\varPhi}_{\varOmega^{k}}^{\dagger} y $

(8)

伪逆矩阵为：

$ \boldsymbol{\varPhi}_{\varOmega^{k}}^{\dagger}=\left(\boldsymbol{\varPhi}_{\varOmega^{k}}^{\mathrm{T}} \boldsymbol{\varPhi}_{\varOmega^{k}}\right)^{-1} \boldsymbol{\varPhi}_{\varOmega^{k}}^{\mathrm{T}} $

(9)

式中Φ^T_Ω^k是支撑集Φ_Ω^k的转置矩阵。

在MMP算法中，分支原子选择的优先级由式(10)决定，其中C_l=[c₁  c₂ … c_K]，优先级函数layer_Order是其算法实现过程:

$ l=\sum\limits_{k=1}^{K}\left(c_{k}-1\right) d^{k-1}+1 $

(10)

式中：c_k与当前k层分支原子的选择相关; d表示当前原子节点的子路径数。C_l给予MMP搜索的第l条路径每个节点原子被选择的优先级，直至残差值满足一定停止迭代阈值ε₁或达到最大搜索路径数L_MAX时停止搜索。

layer_Order函数的实现步骤如下：

1) 输入当前路径次序当前路径次序l，子路径d，稀疏度K，初始值j=1, p=l-1并计算余数c_j=mod(p, d)。

2) 计算p除以d的商并赋值给p，如果j≤K, 更新迭代j值: j=j+1，并转步骤1)；否则停止迭代，输出当前路径节点搜索顺序C_l=[c₁  c₂ … c_K]。

MMP-DF算法步骤如下：

1) 输入观测向量y，测量矩阵Φ，稀疏度K，子路径d，最大路径迭代次数L_MAX和停止迭代阈值ε，初始值l=0, r_min=y, Ω⁰=Ø, r⁰=y。

2) 如果l < L_MAX或‖r_min‖₂>ε，则用当前参数l、d、k由优先级函数layer_Order算法步骤得到当前路径节点搜索顺序序列C_l，更新参数l=l+1，j=0；否则转到步骤6)结束。

3) 计算得到一个当前第c_j大值放入预选集Γ=arg max |〈Φ, r^j-1〉_cj |, 计算候选支撑集Ω^j=Ω^j-1∪Γ, 由最小二乘法求解当前阶段重建信号$ \hat{\boldsymbol{X}}_{j}={\mathit{\pmb{\Phi}}}_{{\mathit{Ω}^{j}}}^† \boldsymbol{y} $，并得到当前残差r^j=y-Φ_Ω^j$ \hat {\boldsymbol{X}}_j $。

4) 如果j < K，且数值长度length(Ω^j) < K，转到步骤3)，否则转到步骤5)。

5) 如果满足‖r^k‖₂ < ‖r_min‖₂, 更新重建信号$ \hat {\boldsymbol{X}} $=X^k和残差r_min=r^k，并转到步骤2)，否则直接转到步骤2)。

6) 输出重建的信号$ \hat {\boldsymbol{X}} $。

MMP算法改进了正交匹配追踪算法支撑集非最优的不足，将正交匹配追踪算法每个原子仅分支出一个原子改为分支出多个原子，其多样化的搜索路径使搜索范围更广，避免陷入局部最优解。理论上可以得到远多于OMP算法的支撑集，大大提高了重构精度。但当遇到大规模数据时，MMP算法需要搜索大规模路径，考虑实际应用需要对搜索路径进行有效优化降低计算复杂度。

2 改进的快速稀疏度自适应多路径匹配追踪算法

在实际应用中，重构信号的稀疏度一般是未知的，借鉴SAMP算法稀疏度自适应和SP算法回溯的思想，结合MMP算法的优点和不足，本文提出了IMMP算法。该算法对内积匹配准则进行改良，利用阈值法选择支撑集，并采用变步长思想来逼近信号稀疏度，并对MMP算法的搜索路径进行剪枝处理，在不影响重构质量的同时极大减少了运算时间，使之更具实用价值。下面主要从以下几个方面对IMMP算法进行描述。

2.1 改进的内积匹配准则

MMP等贪婪算法主要使用内积匹配准则度量测量矩阵与残差值间的相关性强弱。传统的内积匹配准则更多地是从方向上区分差异，而对绝对的数值不敏感，从而影响对相似原子的区分，使观测矩阵中任意2个相近似的原子可能影响残差信号的匹配过程。这种仅考虑向量维度方向上的相似性而没考虑到各个维度量纲差异的不合理性会降低信号的重构质量。假设有任意向量a=(a₁, a₂, …, a_n)和向量b=(b₁, b₂, …, b_n)。改良后的内积匹配准则公式为^[21]：

$ \Re(a, b)=\frac{\langle\boldsymbol{a}, \boldsymbol{b}\rangle}{\|\boldsymbol{a}-\bar{a}\| \cdot\|\boldsymbol{b}-\bar{b}\|} $

(11)

式中$ \bar{a}=\sum\limits_{i=1}^{n} a_{i} / n \text { 和 } \bar{b}=\sum\limits_{i=1}^{n} b_{i} / n $分别表示向量a和b的均值，其可以扩展到多个维度上。改良的内积匹配准则解决了在匹配过程中对绝对数值不敏感导致原始信号的部分信息被丢失的问题，能很好地保存每个向量的原始状态。改良的内积匹配准则相对于原有内积匹配准则更有效地实现了从冗余字典中选择最优的匹配原子。

2.2 剪枝优化搜索路径

MMP算法在基于树的多路径搜索的第k级，与残差值相关性最大的原子将成为新的子元素。假设树的根节点为S¹={S₁¹, S₂¹, S₃¹, …, S_d¹}，初始化的节点原子与初始残差y有最大绝对内积。其中:

$ S_{\mathrm{j}}^{1}=\arg \max \{|\langle\varPhi, \boldsymbol{y}\rangle|, j\}, j=1,2, \cdots, d $

(12)

很容易得到第1个候选支撑集包含的字典原子Ω¹={S₁¹, S₁², S₁³, …, S₁^K}，和相对应的当前残差值{y, r¹, r², …, r^K-1}，其余的候选支撑集按同样的方式选出。与单路径搜索相比，MMP没有只能在单一路径范围内寻找最优解的限制。

从根节点迭代到第K个分支节点, 其中每个节点代表不同的字典原子，当候选支撑集中选择的原子数量等于稀疏度K时，认为此路径搜索完成。根据原子选择优先级的不同，将会得到不同的候选支撑集。取当前起始路径l={1, 2, …, L_MAX}, 在后面每个路径迭代过程中, 一个最有希望的原子被添加到当前路径中。当搜索路径达到最大值L_MAX或者完全确定最优路径时终止迭代。图 1给出了该算法的树搜索过程，其中L_MAX=3，d=3, K=3。可以看到由C₁=(1, 1, 1)确定的原子构成第1个支撑集，第2个支撑集则由C₁=(2, 1, 1)确定的原子构成，直到选择最佳的候选支撑集为止。

在当前支撑集迭代完成后，将与上一个支撑集执行回溯步骤用来选择最佳的K个原子。在第2个候选支撑集{4, 1, 2}完成后，执行回溯细化没有替换任何原子，第3个候选支撑集{6, 3, 4}回溯重新选择原子后，支撑集更正为{6, 2, 4}，提高了路径搜索的精度。回溯虽然带来了额外的计算负担，但优化路径大大降低了计算复杂度，提高了算法的重构精度。

2.3 稀疏度自适应与指数变步长

稀疏度自适应算法设置的固定步长无法兼顾重构精度和重构效率，初始步长较小会降低算法重构效率，步长较大则会降低算法重构精度。采用兼顾两者的变步长是一个不错的方法。根据稀疏度估计的特点，改进算法稀疏度估计可以分为2个部分：1)初始阶段大步长快速估计，2)完成阶段小步长逐步逼近。

文献[22]研究表明, 重构信号的残差能量在2个相邻阶段中是递减的，且随着支撑集长度λ不断逼近真实稀疏度K，残差能量减小速度逐渐变缓趋于稳定。所以本文以相邻阶段重构信号的残差能量作为步长变换的条件。当残差能量小于等于‖r_min‖₂时，步长切换更新稀疏度到下一阶段，当残差能量小于ε₂时稀疏度估计过程完成。ε₂是迭代停止阈值，在这里取ε₂=10^-6。

利用a^x指数函数的变化规律来调整步长，初始阶段步长取一个合适的初值s₀，每次增加步长s=s₀a^x，支撑集长度λ=λ+s；其中a是固定常数，文献[23]证明当a∈[0.3, 0.5]时效果较好，x为当前候选支撑集迭代次数，$ \lceil\cdot\rceil $表示向上取整。当支撑集长度λ接近真实稀疏度K时，添加进支撑集的原子越来越少，直到支撑集长度λ=λ+1。采用此变步长方法可提高算法重构质量和计算效率。

2.4 算法步骤

根据上面对IMMP算法原理的描述，在下面分别详细地给出了阶段自适应估计稀疏度Search函数和IMMP算法的实现步骤：

Search函数的实现步骤如下：

1) 输入搜索序列C_l，残差值r⁰，稀疏度λ，input^j={r^j, Ω^j, j}。

2) 如果j < λ或λ>length(Ω^j)，更新参数j=j+1，否则转到步骤4)。

3) 计算得到一个当前第cj大值放入预选集Γ=arg max |〈Φ, rj-1〉cj |，计算得到候选支撑集Ωj=Ωj-1∪Γ，由最小二乘法求解当前阶段重建信号$ \hat {\boldsymbol{X}}_j $=ΦΩj†y，计算当前残差rj=y-ΦΩj$ \hat {\boldsymbol{X}}_j $，转到步骤2)。

4) 输出重建的信号$ \hat {\boldsymbol{X}} $_λ，残差r^λ，支撑集Ω^λ。

MMP-DF算法步骤如下：

1) 输入观测向量y，测量矩阵Φ，子路径d，初始步长s₀，最大路径迭代次数L_max，停止迭代阈值ε₁, ε₂，并且设定初始值：l=1, j=0, λ=1, x=1, r_min=y, Ω_l⁰=Ø input^j={y, Ø, j}。

2) 如果l < L_max或‖r_min‖₂>ε₁，用当前参数l, d, k由优先级函数layer_Order算法步骤得到当前路径节点搜索顺序序列C_l，用当前参数C_l、y、λ, input^j由阶段自适应估计稀疏度Search函数的算法步骤计算得到参数[X^j, r^j, Ω_l^j, j]，否则转到步骤6)结束。

3) 如果‖r^j‖₂ < ‖r_min‖₂，更新阶段自适应估计稀疏度Search函数输入参数input^j={r^j, Ω_l^j, j}，否则增加指数变量x=x+1, 计算迭代阶段步长s=s₀a^x，更新支撑集长度λ=λ+ $ \lceil s \rceil $，并更新当前最小残差r_min=r^j。

4) 如果l>1或‖r_min‖₂ < ε₂，开始搜索下一路径，路径参数增加l=l+1，置零j、x, 更新当前路径搜索的原子支撑集$ \hat {\boldsymbol{X}} $, 否则转到步骤5)。

5) 如果当搜索路径大于2条l≥2进入原子回溯阶段：更新预选集Ω′_l=Ω_l-1^λ∪Ω_l^λ，计算arg max(| < Φ, Ω′_l> |)，选前λ个最大值，将对应的索引值放入支撑集Γ′中，最小二乘法求解当前阶段重建信号$ \hat {\boldsymbol{X}} $=Φ_Γ′^†y, 计算并更新最小残差r_min=y-Φ_Γ′$ \hat {\boldsymbol{X}} $，并转到步骤2)，否则直接转到步骤2)。

6) 输出重建的信号$ \hat {\boldsymbol{X}} $。

3 仿真验证

本节使用一维信号模拟实验，以信号精确重构比(exact reconstruction ratio, ERR)和重构时间作为算法重构性能的评价标准，其中当重构信号的残差小于阈值ε₁时，认为信号重构是成功的，本文取ε₁=10^-6。信号精确重构率反映了重构精度，其中ERR定义为准确恢复次数p_n与总实验次数p_sum的比值：ERR=p_n/p_sum重构时间和迭代残差反映了重构效率。仿真实验使用Matlab(版本R2018b)在一台装配了Core i5-6400 @2.70GHz处理器的计算机上进行。实验中设置MMP和IMMP算法的最大迭代次数和子路径数(L_MAX, d)分别为(50, 6)和(10, 6)。选择高斯随机矩阵作为观测矩阵，为避免随机因素的干扰，每个算法都独立重复进行了10 000次试验且稀疏信号的非零位置都是随机选择的。

3.1 改进的内积匹配准则的有效性验证

以典型的单路径搜索算法OMP作为MMP算法中的一条搜索路径中为例，分别使用内积准则和改进的内积准则来度量观测矩阵原子和残差信号的相关性，以便挑选最为匹配的原子。随机生成稀疏度K=30，长度N=256的一维信号作为测试信号。比较在加入高斯噪声后(SNR=15 dB), 2种匹配准则下的残差值随不同迭代次数的变化情况，取对数表示的结果如图 2所示。

很明显，随着迭代次数增加，残差值逐渐减小。同等条件下，可以从图 2看出改进的内积匹配准则的残差值比内积准则更低，迭代过程中残差值逼近零的速度更快。这充分说明改进的内积匹配准则不仅更有利于挑选最优匹配原子，降低误选原子的概率，而且能够加速完成迭代逼近，提高重构效率。

表 1是当压缩比M/N=[0.234 4, 0.312 5, 0.396 0] 时，OMP算法应用2种不同原子匹配准则在不同高斯白噪声等级(SNR=[0, 5, 10, 15, 20, 25, 30] dB) 干扰下的平均重构信噪比，更进一步验证了改进的内积匹配准则的性能。其中，重构信噪比公式为：

$ \mathrm{SNR}=10 \lg \frac{\|\boldsymbol{x}\|_{2}}{\|\boldsymbol{x}-\hat{\boldsymbol{x}}\|_{2}} $

(13)

实验表明应用改进的内积匹配准则在稀疏度和测量数2个方面都优于改进前的内积匹配准则。在相同压缩比时，改进的内积匹配准则有更高的重构信噪比；在有相同的重构效果下，改进的内积匹配准则能适应更低的信噪比，抗干扰性能更优。改进的内积匹配准则从向量的夹角信息和长度数值信息2个方面更好地比较了向量间的相关性强弱，是一种更加合理有效的相关性度量准则。

3.2 算法性能分析

本部分比较了IMMP算法与OMP、StOMP、SP、CoSaMP和MMP等5种稀疏信号重构算法的平均重构性能，同时也分析了在相同迭代次数下不同初始步长对IMMP算法性能影响。

图 3(a)给出了固定测量数M=100时，重构信号的ERR性能随不同稀疏度K的变化曲线。从图 3(a)可以看出，本文提出的IMMP算法的ERR要高于传统的重构算法，尤其当K>30时，IMMP算法的ERR明显远高于传统算法。当K=40时，IMMP算法在初始步长s={1, 5, 8, 12}时的ERR为86.59%、87.51%、88.12%和88.79%，而除了CoSaMP算法ERR为0外，OMP、StOMP、SP和MMP算法的ERR分别仅为6.99%、34.91%、4.99%和40.85%。

图 3(b)显示了各种算法对应的CPU平均计算时间，客观上可以比较算法的计算效率。MMP多路径搜索算法计算时间随着稀疏度K的增加呈现近似指数式增长，而OMP、StOMP、SP和CoSaMP算法的计算时间增加幅度不明显。优化搜索路径后的IMMP算法的计算时间较OMP、StOMP、SP和CoSaMP等单路径搜索算法略有增加，但远低于MMP算法。

图 4给出了稀疏度K=30，随着测量数M增大时各种算法的重构成功率。容易由ERR曲线得到，除OMP算法外，其他算法在M>110时能取得极高的重构成功率，而IMMP算法在测量数M>90时也达到了相似的优良性能，所需测量数明显减少。

图 5说明了相应算法的平均计算时间，在M < 110时，IMMP算法显然比MMP算法花费更少的时间；而当M>110，IMMP算法由于要满足迭代停止条件计算时间比其他算法要长，所以IMMP算法在压缩比低的情况下更能体现出性能优势。

实验结果表明，初始步长的大小对IMMP算法的重构性能影响也较为明显。适当大的初始步长不仅能提高算法精确重构成功率，还能有效降低算法的计算时间，提高算法重构效率。IMMP算法具有良好的收敛速度和重构质量，表明优化后的精确的搜索策略可以加速搜索候选支撑集。

4 结论

1) 本文提出的IMMP算法可以盲稀疏度自适应估计，自适应选择支撑集原子个数；而且改进的原子选择方式能在高斯测量矩阵中更准确高效地挑选与残差信号匹配的原子，并结合剪枝技术优化了MMP算法的路径搜索，大大降低了算法计算时间。

2) MMP算法在迭代过程中采用支撑集原子回溯操作和指数变步长方法大大提高了信号重构精度。

3) 仿真结果表明，该算法在不同测量数下性能有所差异, 与传统算法相比，得到的测量数较大时，算法计算时间有所增加，可适当减少迭代次数减少计算时间，而不会降低算法重构精度；在测量数较少时，该算法的信号重建精度高，运行时间短，更贴合实际工程应用。