张家慧，王英志，李新格，沈亮
（1.长春理工大学 电子信息工程学院，长春 130022；
2.北京控制工程研究所，北京 100190）
近几年，由 Candès、Donoho和 Tao等人提出了一种全新的信号采集和处理理论——压缩感知理论（Compressed Sensing，CS）［1-2］。该理论指出：只要信号是稀疏的或在某个变换域是稀疏的，就可以用一个与变换基不相关的观测矩阵将变换所得的高维信号投影到一个低维空间上，然后通过求解一个优化问题从这些少量但包含足够重构信息的投影中以高概率重构出原信号［3］。压缩感知理论突破了传统的采样定理，使用更低的采样率对信号进行随机采样和压缩，节约系统资源，在处理海量复杂信号时具有极大优势，广泛应用于各工程领域中。冲击波信号属于瞬态信号，具有时域持续时间短，有明显的开端和结束节点等特点，信息相对集中，即在整个采集过程中，信息密度较低，可以认为存在某稀疏域能对冲击波信号进行稀疏表示。因此，本文将冲击波信号应用于压缩感知理论中检测重构算法性能。
压缩感知理论的研究内容主要可分为：信号稀疏表示、观测矩阵建立以及重构算法设计。重构算法设计是整个信号处理过程的关键，所以本文围绕重构算法展开研究。目前，常用的重构算法主要有凸优化类算法、组合优化类算法以及匹配追踪类算法等［4］，其中匹配追踪类算法计算量较小，重建效果好，易于实现，应用范围最广。但在使用匹配追踪类算法重构信号时，由于其贪婪特性，虽然在局部搜索时能快速收敛，但缺少全局搜索机制，算法易陷入局部极值解。在冲击波测试领域，对冲击波信号的重构效果要求较高，需要提升重构算法的求解精度。考虑到信号重构过程是一个l0范数最小化问题，而智能优化算法是建立在生物智能或物理现象基础上的随机搜索算法，在没有集中控制且不提供全局模型的前提下，对计算数据的不确定性有很强的适应能力，可以更好地解决最优化问题。鉴于此，本文将搜索能力强的烟花算法引入到压缩感知技术中，提出一种改进型烟花重构算法，减小算法对信号的重构误差，优化算法在同类型算法中的运行效率，提升算法在冲击波测试领域中的实用性。
1.1 压缩感知数学模型
压缩感知理论的核心思想是对具有稀疏性的信号进行观测，降低观测过程中的冗余信息，从而实现以较少数据对原信号的高概率重构。
压缩感知理论应用的必要条件是信号稀疏，若信号x是稀疏的，可以对信号直接进行观测。设Φ∈RM×N，是观测矩阵，其观测过程为：
其中，y是M×1维观测信号，且M< 若是信号x非稀疏，需先对x进行稀疏表示，设Ψ∈RM×N，是标准正交基，θ∈RN×1，是稀疏度为K的稀疏信号，则信号x的稀疏表示过程为： 结合式（1）和式（2）得到非稀疏信号观测过程为： 其中，Α=ΦΨ，是传感矩阵。 当接收到y时，若A满足有限等距性质（Re⁃stricted Isometry Property，RIP），只要观测数满足M=O(Klog(N/K))，就能精准重构出原始稀疏信号。重构模型表示为： 除此之外，还有通过将l0范数松弛为lp（p≤1）范数或基于贝叶斯框架等方法实现信号重构。 在使用智能优化算法实现压缩感知重构过程时，可以使用直接对式（4）的最小化问题进行全局寻优的方法［5-7］：以被测信号长度N作为个体的维度，使用重构误差公式作为适应度函数，个体位置表示为0或1的二进制形式，通过算法的全局寻优能力获得信号非零元素位置信息，再利用最小二乘法获得非零元素位置的幅值。这种重构算法拥有全局搜索能力，可以求得算法的近似最优解，相较于只能寻求到次优解的重构算法，拥有更高的精确度，能以更大概率恢复出原始信号。 虽然这种方法理论上可以求出问题的近似最优解，但仅在N较小时可行，N增大时，算法寻优能力减弱，实际应用能力较差。本文基于此提出一种新型重构算法，相较于前述算法，以信号稀疏度K作为算法维度，［1，N］作为算法解空间，信号中非零位置索引集作为算法的解，算法的结构框图如图1所示。 图1 新型重构算法的结构框图 由图1知，本文的重构算法由前置算法和智能优化重构算法组合形成。信号重构时，观测信号y和传感矩阵A输入重构算法中。通过前置算法，得到从A中挑选出的最匹配y或y残差的列向量索引组成原子最终构造集F，将F传输至智能优化重构算法中，得到算法的维度和一个可行解，补全算法的初始参数后，通过算法全局寻优得到重构信号并输出。其中，因为信号稀疏度不能先验，无法确定算法维度，而稀疏度自适应匹配追踪（Sparsity Adaptive Matching Pur⁃suit，SAMP）算法虽然可以求解稀疏度，但算法步长固定，易产生过估计现象，影响重构效果。本文结合文献［8］中的双阈值变步长思想和文献［9］中的算法回溯思想得到一种改进型SAMP算法，减小对信号稀疏度的估计误差，在算法运行结束时将F提供给基于智能优化的重构算法。而重构算法则使用全局寻优能力求得近似最优解，提升重构算法精确度。重构算法的适应度函数仍为信号的重构误差公式，被测信号稀疏度K为算法维度，信号非零元素位置表示个体位置，F的大小表示算法维度K，F自身可表示一个可行解。 本文提出的新型重构算法，相较于文献［5-7］中的方法，算法维度K< 烟花算法（Firework Algorithm，FWA）是一种基于非生物种群的智能优化算法。算法模拟自然界中烟花爆炸过程进行数学建模，引入随机因素和选择策略，形成一种求解复杂优化问题最优解的并行式全局概率搜索算法。烟花算法一般由爆炸算子、变异算子、映射操作和选择策略四部分组成，其中，爆炸算子是FWA中的重要组成部分，包括对爆炸幅度、爆炸强度和位移操作的计算。虽然FWA提出的时间较晚，但由于算法具有良好的全局搜索方法和局部优化方法，在解决复杂优化问题上很有优势，因此本文使用FWA对信号重构过程进行优化。 为了方便表示，本文将基于FWA得到的重构算法命名为烟花重构算法（Firework Reconstruc⁃tion Algorithm，FWRA）。根据1.2节中本文提出的重构算法结构，FWRA的具体步骤如下所示： 输入：测量矩阵Ф，测量信号y，前置算法求解得到的原子最终构造集F。 初始化：解空间范围Bound∈[1,N]； 步骤一：解空间Bound内，随机生成（Seednum-1）个烟花。 步骤二：适应函数为： 根据式（6）计算烟花的适应度值，其中，适应度值fitness越小表示信号的重构效果越好； 步骤三：根据离散化的爆炸算子产生爆炸火花。 步骤四：为了保证种群的多样性，对任意数量火花进行离散高斯变异操作，产生高斯变异火花。 步骤五：计算适应度值，并找到种群中的适应度值最小的烟花或火花，其代表的空间位置为当代最优适应度值Pbest，与群体最优适应度值Xbest作比较，若Pbest更小则取代Xbest，替换相应可行解。 步骤六：判断是否达到算法最大迭代次数Maxiteration，若未达到则应用选择策略在当前种群内选择烟花进行下一次迭代，否则算法结束。 输出：Xbest和相应的可行解。 FWRA算法的流程图如图2所示。 图2 烟花重构算法流程图 为了进一步提升FWRA的算法性能，本文从两方面对算法进行改进：一方面是利用F中原子个数不变时，原子索引具有一致性的特点，将算法改为非数值优化的求解方式，有利于减小信号的重构误差； （1）非数值优化方法 在相同的实验条件下，不同的重构算法F中原子数固定时，选取的原子索引不变。即若已知测量信号y和原始信号的稀疏度K，在传感矩阵A中仅有唯一一组列向量可以成功匹配出原始信号，不存在其它列向量组合可以再重构出原始信号［10-11］。因而，在通过前置算法已知稀疏度K和一组原子索引的情况下，仅通过对固定解元素排列组合的方式进行求解。而在非数值优化下的FWRA，仍然保留了爆炸算子、变异算子、选择策略，但由于解元素不能超出解空间范围，需要重新设计爆炸算子和变异算子的计算公式并舍弃映射规则。 在爆炸算子中，规定每代所有烟花的爆炸幅度都相同，并使用两元素优化（2-opt）局部搜索方法进行位移操作，2-opt操作过程如图3所示。 图3 2-opt局部搜索操作 得到新可行解的适应度值f（Xg），与原解f（Xi）相比，若适应度值更小，则可以保留下来，若适应度值更大，也有可能被一定概率接受，设接受差解的概率为Pa，即： 其中，ρ为控制参数，由式（7）知，如果两个适应度值越相近，控制参数ρ越小，则接受的概率越高。这样的设定是为了使算法有一定的几率接受局部改变，同时减少一些无意义的改变。如果生成的火花与烟花相距太远，两个适应度值相差过大，那么该操作基本无效。参数ρ控制接受概率Pa的大小，防止其过大或过小。 在变异算子中，使用一种介于2-opt和3-opt之间的2h-opt局部搜索策略，优化解中的局部点，且变异算子中不考虑定义接收差解的概率，以强化算法的搜索方向，2h-opt操作过程如图4所示。 图4 2h-opt局部搜索操作 （2）精简爆炸强度 FWRA的爆炸强度公式如下： 其中，m0表示预设火花总量最大值； 使用FWRA解决信号重构问题时，算法的适应度值f（Xi）的数量级也在1e-15左右，此时爆炸强度公式中存在的极小常数ɛ不能忽略，致使每次烟花爆炸产生的火花数量都比预分配值多，最后分配到的火花总量远大于预设的火花总量最大值，严重影响算法的运行效率。为了避免ɛ的影响，本文提出一种精简爆炸强度公式如下： 其中，m表示更新后的火花总量最大值。 式（9）将适应度值差值分为零值与非零值两种情况分配火花，解决了公式无效化问题。同时，为了避免最优解烟花个数不唯一时火花总数远超预设值的情况，对m进行非实时更新，使生成的火花总量接近预设值。同时，相较于式（8）中最优烟花下分配到的火花数量过度，易产生重复解，浪费资源的情况。本式按比例合理分配火花数量，在保证算法搜索能力的前提下提升算法的运行效率，节约算法资源。 （3）自适应爆炸幅度 在算法离散优化过程中，由于局部函数内空间差异性较大，比较一次迭代是否产生更优解没有意义，因为当群体陷入局部收敛时，往往需要经多次变化才有机会跳出局部极值，这需要多次的迭代和计算代价。参考动态搜索烟花算法（dynFWA）中的爆炸幅度公式［12］，本文提出一种新的自适应爆炸幅度公式，用式（7）中控制参数ρ表示，设若连续15代不改变，则减小ρ，反之，则增大ρ，如式（10）所示： 其中，C+表示增大因子； 在群体最优适应度值Xbest比当代最优适应度值Pbest小时，爆炸产生更优值，扩大爆炸幅度加强全局搜索； （4）精英-轮盘赌选择策略 FWRA中基于欧式距离度量的选择策略，在每一代构建群体中任意两点间欧式距离矩阵，消耗算法的运行时间，且在离散化算法中保留优秀解的能力差，不适合使用。因此，除了继续应用精英选择策略，保留每代最优解，对其它个体采取轮盘赌的选择方式，提高优秀解被选择的概率，但不排除掉全部差解。定义烟花或火花Xi被选择的概率为： 由式（11）得到除最优个体的选择概率后，再根据轮盘赌策略选择子代。由定义可知，适应度值越小的个体有越大的概率被选择进入下一次迭代中，适应度值越大的个体性能虽然较差，但有利于发展种群的多样性，所以也存在着小概率的入选机会。 结合FWRA，可得IFWRA的具体步骤如下： 输入：测量矩阵Ф，测量信号y，前置算法求解得到的原子最终构造集F。 初始化：解元素固定为F中的原子索引； 步骤一：使用固定的解元素，随机生成（Seednum-1）个烟花。 步骤二：根据式（6）计算烟花的适应度值。 步骤三：首先使用式（9）计算每个烟花产生的火花数量，再使用式（10）自适应爆炸幅度，接着使用2-opt方法进行位移操作，最后结合式（7）得到爆炸火花。 步骤四：对任意数量火花使用2h-opt方法进行局部点变异操作，产生变异火花。 步骤五：计算适应度值，找到种群中当代最优适应度值Pbest，与群体最优适应度值Xbest作比较，若Pbest更小则取代Xbest，并替换相应可行解。 步骤六：判断是否达到算法最大迭代次数Maxiteration，若未达到则结合式（11）的精英-轮盘赌选择策略生成下一代烟花，否则算法结束。 输出：Xbest和相应的可行解。 IFWRA的算法流程图如图5所示。 图5 IFWRA算法流程图 IFWRA基于FWRA，使用非数值优化方法固定算法解空间，去除爆炸强度公式中极小常数导致的实际火花总量冗余影响，自适应烟花的爆炸幅度防止算法陷入局部最优，并使用耗时更少的精英-轮盘赌选择策略，有效提升了算法的重构性能。相较于FWRA，IFWRA对信号的重构误差更低，算法运行效率更高。 本文使用高斯稀疏信号作为测试信号，检验IFWRA的算法性能。IFWRA作为新型重构算法，在测试时，使用FWRA比较IFWRA中改进方法的有效性，蚁群重构算法（Ant Colony Reconstruc⁃tion Algorithm，ACRA）检验IFWRA的搜索能力，匹配追踪类算法中的SAMP算法对照IFWRA近似最优解的高精确性。其中，所有新型组合算法均使用相同前置算法，避免影响智能优化算法的搜索能力，使用均方根误差表示重构误差，设计实验一。 实验一：信号长度N=256的高斯稀疏信号作为被测信号进行实验，测量数M=128，稀疏度K=15∶5∶60，测量矩阵为伯努利随机矩阵，每个K值测试CNT=50次，分别使用IFWRA、FWRA、ACRA和SAMP（S=8）算法对信号进行重构，观察算法对信号的重构误差和运行效率。其中，FWRA和IFWRA中部分初始参数为：Seednum=5； 由图6可知，IFWRA的重构误差在稀疏度K从15到60的变化过程中始终最低，重构误差均值为5.982 4e-16，比次优的ACRA低约18.72%。IFWRA相较于FWRA，能避免算法陷入局部最优，拥有更高的搜索能力，重构误差比FWRA算法低约38.93%。而SAMP算法的重构误差在同类型匹配追踪类算法中较低，但相较于组合算法最高，误差均值为1.941 0e-15，IFWRA的重构误差相较于SAMP算法低约69.18%。 图6 算法重构误差与稀疏度K的关系 选取K=15和K=30时的 IFWRA、ACRA和FWRA每代最优解均值进行比较，观察智能优化算法的全局寻优能力。 从图7、图8中K=15和K=30显示的算法最优适应度值曲线的变化情况来看，IFWRA在三种智能优化算法中全局寻优能力最强。IFWRA和ACRA的曲线收敛速度较快，最优适应度值也一直远低于FWRA。而IFWRA的曲线相较于ACRA始终最低，证明IFWRA搜索能力最强。 图7 K=15时，算法最优适应度值与迭代次数的关系 图8 K=30时，算法最优适应度值与迭代次数的关系 由图6可知，SAMP算法的重构误差数量级在1e-15左右，而组合算法的重构误差数量级在1e-16左右，可见，SAMP算法无法满足对信号高精度重构的要求。且新型重构算法以牺牲部分算法运行时间为代价，寻求更小的重构误差，相较于SAMP算法，在运行效率上稍显劣势。所以，考虑到实际应用性与公平性，不再使用SAMP算法作对比。 由图9可知，在新型重构算法中，IFWRA的运行效率始终最高，平均运行时间为13.789 9 s，相较于平均运行时间为30.406 9 s的FWRA，运行效率提升约54.65%。ACRA耗时最长，平均运行时间为63.737 1 s，IFWRA相较于ACRA的运行效率提升约78.36%。ACRA算法复杂度约为O（D4），其中D表示算法维度，而IFWRA与IFWRA算法整体计算结构不变，算法复杂度均约为O（T2），其中T表示算法粒子数，所以ACRA的运行时间随K（K=D）增大呈四次幂递增，FWRA和IF⁃WRA的运行时间虽然也随K的增大而增加，但由于T基本无变化，算法的时间运行曲线增长速率较慢，所以ACRA仅在K=15和20时相较于FWRA稍具优势，当K>20后，算法运行时间迅速增大。而IFWRA相较于FWRA，算法运行时产生的火花总数更低，且使用了耗时更小的选择策略，算法计算量更少，运行效率更高。 图9 算法运行时间与稀疏度K的关系 在实际应用中，准确测量冲击波信号主要是为了合理规划武器装备与测试环境，降低武器测试风险［13］。所以，对冲击波信号进行信号重构时不要求实时性，反而更重视精确性。本文截取50 psi量程传感器实测冲击波信号，进行精确截取、降频等预处理，得到信号长度N=4 096的冲击波时域信号，再使用离散小波矩阵对信号进行表示，得到小波域冲击波信号［14］。冲击波的时域和小波域图形如图10所示。 图10 50 psi实测冲击波信号小波基稀疏结果 由图10可知，图10（a）中的时域冲击波信号是非稀疏信号，但图10（b）中的小波域冲击波信号波形具有稀疏性，表明离散小波变换可以实现冲击波信号的稀疏化。 本文为了验证IFWRA在冲击波测试领域的应用性，以冲击波信号为测试信号，比较IFWRA和其他同类型重构算法对信号的重构误差和运行效率，设计实验二。 实验二：截取一组50 psi量程传感器实测冲击波信号作为被测信号进行实验，信号长度N=4 096，测量数M=2 048，稀疏字典为离散小波变换矩阵，测量矩阵为伯努利随机矩阵，对信号进行CNT=50次测试，分别使用IFWRA、FWRA、ACRA对信号进行重构，观察算法的重构误差和运行效率。其中，FWRA和IFWRA部分初始参数：Seednum=5； 结合表1和图11—图13可知，IFWRA可应用于冲击波测试领域中实现对信号的重构，且在同类型组合算法中，IFWRA算法性能最优。IF⁃WRA在图11中的误差曲线相较于图12和图13中的算法误差曲线整体波动最小，表明重构误差最小，为1.148 8e-18，相较于次优的ACRA减小8.15%，相较于FWRA减小38.88%。IFWRA的算法运行效率也最高，算法运行时间为116.650 2 s，相较于次优的ACRA提升32.34%，相较于FWRA提升37.53%。分析图9知，ACRA仅在K较小时存在算法耗时低于FWRA的情况，表明本文选取的50 psi冲击波信号在小波域上稀疏性能好，稀疏度K极小。 表1 算法的运行时间和对信号的重构误差 图11 使用IFWRA对冲击波信号进行重构 图12 使用ACRA对冲击波信号进行重构 图13 使用FWRA对冲击波信号进行重构 IFWRA与其它重构算法在对比实验中算法性能依旧最佳，实验结果虽然与实验一相仿，但意义不同。实验一是为了验证IFWRA的算法性能，证明对IFWRA改进方法的有效性。而实验二主要是为了验证IFWRA的实用性。当IFWRA以实测冲击波信号为测试信号仍实现精确重构时，肯定了IFWRA的应用范围，证明IFWRA的研究技术较为成熟，已经可以应用于工程实践中。同时，实验二也为压缩感知技术在以冲击波为代表的瞬态信号领域的应用提供研究和参考价值，意义更重大。 本文提出了一种改进型烟花重构算法，在使用烟花算法的全局寻优能力优化信号重构过程的基础上，从算法的求解方式和算子缺陷两方面对算法进行改进：将算法改为非数值优化求解方式，固定解空间，增强算法搜索能力；
1.2 智能优化重构算法
2.1 烟花重构算法
烟花个体维度Dim=K；
初始烟花数量Seednum；
第一个烟花X1=F；
火花总数Sonnum；
算法最大迭代次数Maxiteration。
反之，表示信号的重构效果越差。
2.2 对烟花重构算法的改进
另一方面是对算法中的重要算子进行改进。算法的爆炸强度、爆炸幅度和选择策略公式在解决信号重构问题时存在缺陷，为避免计算资源的浪费，本文在保留爆炸强度特点的同时，消除公式中极小常数对爆炸火花数量的影响，自适应爆炸幅度平衡全局和局部搜索，并使用耗时更少的轮盘赌选择策略，提升算法的精确度和运行效率。为了将改进后的FWRA与FWRA作区分，将改进后的FWRA命名为IFWRA（Improved FWRA）。
fmax表示最大适应度值，ɛ表示数量级约为1e-15的极小常数。
C-表示缩小因子；
Num为记录烟花的解连续不改变的次数。
在Xbest连续15次比Pbest大时，烟花陷入局部最优，减小爆炸幅度，加强局部搜索，让烟花有更大可能跳出局部极值。
2.3 改进型烟花重构算法
烟花个体维度Dim=K；
初始烟花数量Seednum；
第一个烟花X1=F；
火花总数Sonnum；
算法最大迭代次数Maxiteration。
Sonnum=50；
Maxiteration=500；
X1=F。ACRA 中部分初始参数为：蚁群数量Antnum=50；
Maxiteration=500；
第一只蚂蚁位置A1=F。
Sonnum=50；
Maxiteration=500；
X1=F。ACRA中初始参数为：Antnum=50；
Maxiteration=500；
A1=F。
精简爆炸强度，自适应爆炸幅度和使用基于精英-轮盘赌的选择策略，有效节约算法计算资源。使用高斯稀疏信号和50 psi量程常感器实测冲击波信号，将本文算法与同类型和传统重构算法进行对比，测试结果表明：在算法对信号的重构误差比较上，本文算法始终优于其他对比算法；
在算法运行效率的比较上，本文算法的运行效率在同类型算法中始终最高。