CN103985093B - 基于多随机测量迭代像素判决的压缩感知稳健重构方法 - Google Patents

Abstract

本发明属于噪声条件下的压缩感知重构处理领域，具体涉及一种基于多随机测量迭代像素判决的压缩感知稳健重构方法，具体包括，获取成像区域图像，对成像区域图像进行成像网格划分，根据成像模型构造稀疏字典；生成初始测量矩阵，并从原始数据中获取初始观测数据；通过随机抽取初始测量矩阵若干行组成测量子矩阵，并计算获得观测子集，根据测量子矩阵、稀疏字典和观测子集计算正则化参数，进一步求解稀疏系数，根据TSC算法，得到稀疏系数的估计值，当稀疏系数满足初始化判决门限值后，即得到了压缩感知稳健重构的图像结果，本发明利用噪声的随机性通过多次随机测量及TSC处理，采用自适应迭代处理，有效改善了CS对含噪信号的重构质量。

Description

基于多随机测量迭代像素判决的压缩感知稳健重构方法

技术领域

本发明属于噪声条件下的压缩感知(Compressed Sensing,CS)重构处理领域，具体涉及一种基于多随机测量迭代像素判决的压缩感知稳健重构方法。

背景技术

CS是近几年发展起来建立在信号稀疏表示和逼近理论基础上的新的研究领域，它与传统奈奎斯特采样定理不同，理论指出：只要信号是可压缩的或在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号，可以证明这样的投影包含了重构信号的足够信息。CS理论的本质是一种非适应性的、非线性的信号(图像)重构方法，对于一个N×1维离散信号f在正交基Ψ＝[ψ₁,ψ₂,…,ψ_n]下可以展开为式中s为f的展开系数，s_i＝<f,ψ_i>。假设s中非零元素的个数K(K＜＜N)，那么称s为K稀疏的。假设利用一个J×N维测量矩阵Φ采集数据，那么获取的数据为y＝Φf＝ΦΨs。通常情况下，为了从y中完全恢复s，J应该大于或等于N。但是如果s稀疏并且ΦΨ满足有限等距特性(Restricted Isometry Property,RIP)时，f则可以从J＝cK log(NK)次测量值以高概率精确重建，只需求解下列凸规划问题

min||s||₀   s.t.y＝ΦΨs          (1)

其中，c代表一个常数，s.t.表示满足约束条件的数学符号(本文下同)，由于求解式(1)的数值计算极不稳定而且是NP(Nonconvex Programming)难题，而该最优化问题与信号的稀疏分解十分类似，所以有学者从信号稀疏分解的相关理论中寻找更有效的求解途径。由于||s||₁最小范数下在一定条件下和||s||₀最小范数具有等价性，那么就可以将上式转化为||s||₁最小范数下的最优化问题：

min||s||₁  s.t.y＝ΦΨs       (2)

但在实际环境中，重构信号通常会受到噪声的污染而不再稀疏，如果对含噪信号仍然采用单一的稀疏性约束条件则无法有效恢复原始信号。但CS理论仍然可以采用其他的信号恢复方法，不同之处在于恢复过程所使用的优化目标函数的形式及参数的设置不同。

为解决这一问题，现阶段研究较多的有CSDN算法、LASSO算法、算法、FOCUSS(FOCal Underdetermined System Solver)算法等。

(1)当噪声分布已知的情况，沿用BP方法对噪声的抑制方法，即CSDN(Compressed Sensing De-Noising)算法，修改式(2)约束条件为

min||s||₁  s.t.||ΦΨs-y||₂＜σ       (3)

其中σ是与噪声方差及数据长度N相关的常量。

(2)当已知信号的稀疏度(1-范数大小)时，则可以采用LASSO(Least AbsoluteShrinkage and Selection Operator)算法来对信号进行有效恢复

min||ΦΨs-y||₂  s.t.||s||₁＜τ        (4)

其中τ为稀疏度。

(3)当对信号和噪声都未知时，可把寻找稀疏解问题归结为带约束二次规划(Bound Constrained Quadratic programming,BCQP)问题，可以利用梯度投影(Gradient Projection,GP)算法来有效求解

min(||ΦΨs-y||₂+λ||s||₁)         (5)

其中，λ称为正则化参数(regularization parameter)，λ的选择是正则化处理的基本问题，可以利用测量数据通过一些算法进行有效估计估计，如Stein无偏风险估计(Stein’s Unbiased Risk Estimator,SURE)算法、广义交叉验证(GeneralizedCross-Validation,GCV)算法以及用于图像处理的L曲线(L-curve)算法(O.Batu,M.Cetin,Parameter selection in sparsity-driven SAR imaging,IEEE Transactions onAerospace and electronic systems,October2011,47(4):3040-3050.)。

尽管上述重构方法采取了不同的约束形式，但它们都采用2-范数对噪声进行抑制。由于无法体现信号的稀疏性，2-范数易于造成重构系数的幅度损失。幅度损失是大部分CS重构算法都很难解决的问题，近年来，为方便研究大部分研究人员甚至有意回避了该问题。显然，上述各式都是对含噪信号稀疏重构的有效形式，但由于噪声破坏了信号在基矩阵Ψ上的稀疏度，当信噪比较低且未知时，仅仅利用单次测量，即使通过式(5)也很难实现对原始信号的精确重构。

发明内容

本发明解决的技术问题是针对低信噪比条件下基于单次随机测量的CS算法重构含噪稀疏信号稳健性差的不足，提出了基于多随机测量迭代像素判决(Multiple Measurements Iterative Pixel Discrimination,MMIPD)的CS稳健重构方法；

本发明采用的技术方案是：

一种基于多随机测量迭代像素判决的压缩感知稳健重构方法，包括以下步骤：

(S1)获取成像区域图像，对成像区域图像进行成像网格划分，根据成像模型构造稀疏字典Ψ；生成初始测量矩阵Φ₀，并从原始数据x中获取初始观测数据y₀，即：y₀＝Φ₀x；定义l表示重构序号，初始化l＝1，l取值为自然数；初始化判决门限ε_S与ε_SD的值；

(S2)随机抽取Φ₀若干行组成测量子矩阵Φ_l，计算获得观测子集y_l，即

y_l＝Φ_ly₀；

(S3)若l等于1，则根据Φ_l、Ψ及y_l的值和GCV算法计算正则化参数λ，若l值不等于1，则进入第(S4)步骤；

(S4)将Φ_l、Ψ、y_l、λ代入min(||Φ_lΨs_l-y_l||₂+λ||s_l||₁)公式，求解得到稀疏系数s_l；

(S5)若l等于1时，则l值增加1后，返回步骤(S2)进行迭代重构；

若l不等于1，根据时间序列聚类算法，得到稀疏系数s的估计值即：

$\hat{s}=\mathrm{TSC}[s_1,s_2,\&CenterDot;\&CenterDot;\&CenterDot;,s_l];$

计算当前图像稀疏度与前一次图像稀疏度两者差的对数值ΔI_S，

其中，表示第l次重构得到的中所有元素值超过均值的个数；

计算当前图像标准差与前一次图像标准差两者之差的对数值ΔI_SD，

其中，其中，表示第l次重构得到的的标准差；

如果ΔI_S≤ε_S且ΔI_SD≤ε_SD，迭代终止，所获取的即为压缩感知稳健重构的图像结果；

否则，则l值增加1后，返回步骤(S2)进行迭代重构。

优选地，所述判决门限ε_S与ε_SD的取值范围为-30dB～0dB。

本发明在单次观测数据的基础上，通过多次随机抽取和CS处理获得重构序列，再利用时间序列聚类(Time Series Clustering,TSC)算法完成对原信号的稳健重构，有效地抑制了重构杂波，提高了含噪信号的重构稳健性。

采用本发明可以达到以下技术效果：

本发明利用噪声的随机性通过多次随机测量及TSC处理，在遵循CS基本原理的基础上采用自适应迭代处理，有效改善CS对含噪信号的重构质量，为噪声条件下的稳健稀疏重构拓展了新思路。本发明的应用并不需要完备的原始数据集，具有良好的适用性，可用于实际工程中的信号处理。需要指出的是，本发明并非一种特定的含噪信号稳健重构方法，而是一个稳健重构框架，其所采用的CS重构算法和TSC算法均可利用相似的算法进行替代，因此具有更优的适用性和通用性。

附图说明

图1是本发明方法流程图；

图2是本发明仿真试验场景图；

图3是部分无噪声回波数据的CS重构处理结果图；

图4是部分含噪声回波数据的CS重构处理结果图；

图5是全部含噪声回波数据的CS重构处理结果图；

图6是不同迭代次数下的多次测量及相干处理结果图；

图7是不同迭代次数下的多次测量及时间序列聚类处理结果图；

图8是两种方法对实测数据处理后的结果图，图8-a是传统的CS方法处理后的结果图，图8-b是本发明处理后的结果图。

具体实施方式

为更加清楚阐明本发明原理过程，下面结合附图和实施例作进一步详细说明。

图1为本发明流程图，第一步：对成像区域图像进行成像网格划分，根据成像模型构造稀疏字典Ψ；生成初始测量矩阵Φ₀，并依此从雷达获取的原始回波数据x(即满足Nyquist采样定理的完备数据集)中获取初始观测数据y。

具体步骤为：

步骤①：根据预先设定的成像分辨率(ρ_r与ρ_a分别为距离向与方位向分辨率)对长L_r米、宽L_a米的成像区域进行稀疏成像网格划分，得到N_r×N_a成像空间矩阵，其中N_r＝L_r/ρ_r，N_a＝L_a/ρ_a，然后根据信号体制、成像模型构造稀疏字典Ψ(GurbuzAli Cafer,McClellan James H.,Scott Waymond R.A compressive sensing dataacquisition and imaging method for stepped frequency GPRs[J].IEEE Transactionson Signal Processing,2009,57(7):2640-2650.)；

步骤②：生成0、1随机矩阵作为初始测量矩阵Φ₀(Gurbuz Ali Cafer,McClellan James H.,Scott Waymond R.A compressive sensing data acquisition andimaging method for stepped frequency GPRs[J].IEEE Transactions on SignalProcessing,2009,57(7):2640-2650.)，并依此对原始数据x进行观测得到初始观测数据y₀，如下式所示：

y₀＝Φ₀x              (6)

定义l表示重构序号，初始化l＝1，l取值为自然数；

第二步：随机抽取初始测量矩阵Φ₀的若干行构造成测量子矩阵Φ_l，然后对初始测量数据y₀进行观测，获得观测子集y_l，再利用GCV算法估计正则化参数λ。

具体步骤为：

随机抽取初始测量矩阵Φ₀的M行(其中M要小于矩阵Φ₀的最大行数)作为测量子矩阵Φ_l，其中l表示重构序号，l取值为自然数。然后利用Φ_l对y₀进行子测量以获得测量子集y_l，具体如下式所示：

y_l＝Φ_ly₀        (7)

第三步：若l等于1，根据Φ_l、Ψ及y_l利用GCV算法估计正则化参数λ；若l值不等于1，则进入第第四步；

第四步：将Φ_l、Ψ、y_l及λ代入公式min(||Φ_lΨs_l-y_l||₂+λ||s_l||₁)公式，求解得到稀疏系数s_l；

第五步，若l等于1时，则l值增加1后，返回步骤(2)进行迭代重构；若l不等于1，利用TSC算法(Keogh E,Lin J,Truppel W,Clustering of time seriessubsequences is meaningless:implications for previous and future research[C],Proceedings of the IEEE International Conference on Data Mining,Melbourne,FL,2003:115-122.)对重构序列[s₁,s₂,…,s_l]中相同位置的元素进行聚类处理，得到稀疏系数s的估计值即

$\hat{s}=\mathrm{TSC}[s_1,s_2,\&CenterDot;\&CenterDot;\&CenterDot;,s_l]\&DoubleRightArrow;\left[\begin{array}{c}\hat{s}\left(1\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \hat{s}\left(N_r{N}_a\right)\end{array}\right]=\left[\begin{array}{c}\mathrm{TSC}[s_1\left(1\right)\&CenterDot;\&CenterDot;\&CenterDot;s_l\left(1\right)]\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{TSC}[s_1\left(N_r{N}_a\right)\&CenterDot;\&CenterDot;\&CenterDot;s_l\left(N_r{N}_a\right)]\end{array}\right]---\left(8\right)$

本发明定义图像稀疏度表示第l次重构得到的中所有元素的值超过均值的个数；定义重构图像标准差表示第l次重构得到的的标准差。

为提高MMIPD方法的处理效率，在实际处理中通常采用变量I_S与I_SD各自当前值与前一次所得值的差值对数形式，即和并设定迭代终止判决门限ε_S与ε_SD(ε_S与ε_SD的取值可根据实际需求选取范围为-30dB～0dB之间值，通常选为-10dB)。

若满足条件：ΔI_S≤ε_S且ΔI_SD≤ε_SD，则MMIPD方法处理结束，迭代终止，稀疏系数s的估计输出为：

$\hat{s}=\mathrm{TSC}[s_1,\&CenterDot;\&CenterDot;\&CenterDot;,s_l]---\left(9\right)$

所获取的即为压缩感知稳健重构的图像结果；

否则，重构序号l值增加1，返回至第二步，继续计算稀疏系数s_l。

具体实施例结合附图进行说明如下：

图2为采用宽带步进频率体制的合成孔径雷达成像的场景图：以雷达航迹方向为方位向，设为X轴，垂直于航迹方向为距离向，设为Y轴，雷达孔径中点O为坐标原点，以包含目标和雷达孔径的平面为坐标平面建立直角坐标系。雷达孔径长度10米，方位采样间隔为0.1米，孔径中心距离成像区域中心径向距离为8米，成像区域为[-1～1]米×[7.5～9.5]米，共包含5个理想点目标，坐标分别为(0.3,9.2)、(-0.5,8.6)、(-0.2,7.8)、(0.6,8.0)、(0.4,8.2)，设雷达发射的步进频率脉冲信号，其频率步进增量Δf＝4MHz，频点总数N_f＝300，利用Matlab仿真平台仿真全频带、全通道原始回波数据x；

根据图2设定的成像区域可知，成像区域长度L_r＝2米、成像宽度L_a＝2米，设定稀疏成像距离向与方位向的分辨率分别为ρ_r＝0.1米、ρ_a＝0.1米，那么成像空间矩阵大小为N_r×N_a＝201×201＝40401，然后根据图2中步进频率信号的相关参数构造稀疏字典Ψ；

图3是利用随机测量的10％数据通过式(5)求解的重构结果图，具体过程是通过Matlab仿真软件来实现的，见参考文献“M Grant and S Boyd,CVX:MatlabSoftware for Disciplined Convex Programming(Web Page and Software)2008[online].Available:http://stanford.edu/boyd/cvx.”。

图4是对有噪数据(信噪比为-15dB)随机抽取10％利用式(5)所得的重构结果图，其中正则化参数λ是利用观测数据通过GCV算法估计得到的。可以得出：虽然5个目标点的稀疏系数均已重构出来(虚线圆圈所包围的点)，但均存在不同程度的幅度损失(目标点的色度偏灰)，而且重构图像中存在较多的杂波，其原因是噪声破坏了信号的稀疏性，增加了稀疏系数中非零值的个数。

图5是利用全部含噪数据通过求解式(5)所得的重构结果图，其中正则化参数λ同样是通过GCV算法估计所得，但利用了全部观测数据。可以得出：用以重构的输入数据量越大，目标信息量也越多，就越有利于重构质量的改善。与图4对比可知：基于全部数据的稀疏重构，其杂波数量较图4明显减少，并且目标稀疏系数幅度损失也相应减少(目标点的色度较亮)。

图6是采用多次测量及相干处理后的重构结果图，其中每次迭代重构均采用采用同分布的不同随机测量矩阵(用于重构的输入数据量为随机抽取10％)，通过求解式(5)得图6所示结果，图6-a、图6-b、图6-c、图6-d分别对应的迭代次数为5、10、15、20。可以看出：相干处理虽然在常规雷达信号处理中具有较好处理效果，但是其平均加权的特性使得该方法很难有效提取目标或杂波的稳定序列聚类中心，虽然可以保证目标高概率的重构并减少幅度损失，但其杂波抑制能力有限(输出图像信噪比不高)。另外，在设定的判决门限下，利用相干处理方法即使迭代了20次，重构图像质量指标依然无法满足迭代终止条件，因此相干处理需要更多的迭代次数，增大了存储与计算的负担。

图7是采用多次测量及时间序列聚类处理后的重构结果，其中每次重构均采用与图6相同的测量矩阵并通过求解式(5)获得，图7-a、图7-b、图7-c、图7-d分别对应的迭代次数为5、10、15、20。对比各图可以看出：当迭代次数达到15时，目标稀疏系数基本重构处理，且幅度损失较小，杂波抑制效果较好。重构过程中，在设定的判决门限下，当迭代次数为17时满足终止条件，迭代重构结束。对比图6与图7可知：本发明在测量矩阵、正则化参数及迭代次数相同的情况下，更有效的提取了目标或杂波稳定的幅度聚类中心，改善了稀疏重构质量，提高了输出图像信噪比。

图8为实测浅埋金属圆柱体的处理结果图，金属圆柱体的双散射中心特性可参考文献“JIN Tian,LOU Jun,ZHOU Zhimin,Extraction of Landmine FeaturesUsing a Forward-Looking Ground Penetrating radar with MIMO Array,IEEETransactions on Geoscience and Remote Sensing,2012.10,50(10):4135-4144.”利用20％数据，图8-a给出了常规CS算法的重构结果，图8-b为本发明方法处理结果图，对比发现，本发明重构图像中的杂波较少，图像信杂比较高。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (2)

1.一种基于多随机测量迭代像素判决的压缩感知稳健重构方法，其特征在于包括以下步骤：

(S1)获取成像区域图像，对成像区域图像进行成像网格划分，根据成像模型构造稀疏字典Ψ；生成初始测量矩阵Φ₀，并从原始数据x中获取初始观测数据y₀，即：y₀＝Φ₀x；定义l表示重构序号，初始化l＝1，l取值为自然数；初始化判决门限ε_S与ε_SD的值；

(S2)随机抽取Φ₀若干行组成测量子矩阵Φ_l，计算获得观测子集y_l，即y_l＝Φ_ly₀；

(S3)若l等于1，则根据Φ_l、Ψ及y_l的值和GCV算法计算正则化参数λ，若l值不等于1，则进入第(S4)步骤；

(S4)将Φ_l、Ψ、y_l、λ代入min(||Φ_lΨs_l-y_l||₂+λ||s_l||₁)公式，求解得到稀疏系数s_l；

(S5)若l等于1时，则l值增加1后，返回步骤(S2)进行迭代重构；

若l不等于1，根据时间序列聚类算法，得到稀疏系数s的估计值，即：

其中，TSC表示时间序列聚类算法；

计算当前图像稀疏度与前一次图像稀疏度两者差的对数值△I_S，

其中，表示第l次重构得到的中所有元素的值超过均值的个数；

计算当前图像标准差与前一次图像标准差两者之差的对数值△I_SD，

其中，表示第l次重构得到的的标准差；

如果△I_S≤ε_S且△I_SD≤ε_SD，迭代终止，所获取的即为压缩感知稳健重构的图像结果；

否则，则l值增加1后，返回步骤(S2)。

2.如权利要求1所述的基于多随机测量迭代像素判决的压缩感知稳健重构方法，其特征在于：所述判决门限ε_S与ε_SD的取值范围均为-30dB～0dB。