CN101640541A - 一种稀疏信号的重构方法 - Google Patents

Abstract

本发明公开了一种稀疏信号的重构方法，主要解决从观测向量重构原始稀疏信号速率低的问题。该方法利用约束目标函数的可分解性，将优化约束目标函数问题分解为一系列小的约束目标函数来优化，提高重构速率。它包括获取原始数据部分、重构原始稀疏信号部分和重构准确度评价部分，其中：获取原始数据部分包括产生原始稀疏信号、观测矩阵和观测向量；重构原始稀疏信号部分主要包括设定无约束目标函数，推导约束目标函数、使用序列最小优化方法分解约束目标函数、计算重构信号和对重构信号进行除偏处理；重构准确度评价是比较均方误差的大小。本发明的重构方法可在保证重构准确率的前提下，提高重构速率，可应用于压缩感知等领域的稀疏信号重构问题。

Description

一种稀疏信号的重构方法

技术领域

本发明属于信号处理技术领域，特别是涉及稀疏信号重构，可用于压缩感知或者相反问题上。

背景技术

随着数字化信息时代的发展，信息来源的模拟化和信息处理工具的数字化越来越占有主导地位，信号的采样是两者联系的桥梁。傅立叶变换和奈奎斯特采样定理表明，欲无失真从离散信号中恢复信号，其采样速率必须是信号带宽的2倍。然而随着人们对信息需求量的增加，携带信息的信号的带宽越来越宽，据此理论为基础的信号处理框架，要求的采集速率和处理速度越来越高，对宽带信号处理的困难也越来越大，给信号采集、存储、传输和处理等带来了巨大的压力。

近几年来出现的一种新颖的理论——Compressed sensing，也称为Compressivesampling，得到了快速发展。该理论异于近代奈奎斯特采样定理，它指出：利用随机观测矩阵可以把一个稀疏或可压缩的高维信号投影到低维空间上，然后再利用这些少量的投影通过解一个优化问题就可以以高概率重构原始稀疏信号，并证明了这样的随机投影包含了原始稀疏信号的足够信息。于是一个可压缩信号能够用比奈奎斯特采样要求的信号长度小得多的信号来准确表示。该理论利用其它变换空间描述信号，建立新的信号描述和处理的理论框架，使得在保证信息不损失的条件下，用远低于奈奎斯特采样定理要求的速率采样信号而又可以完全恢复信，可极大地降低信号的采样频率、极大地减少数据存储和传输代价、极大地减少信号处理时间和计算成本。它带领信号处理进入一个新的革命时代。

在压缩感知理论中，由于观测数量K远小于信号长度N，因此不得不面对求解欠定方程组y＝Ax∈R^K，K＜N的问题。乍一看，求解欠定方程组似乎是无望的，但是，人们发现由于信号x∈R^N是稀疏的或可压缩的，也就是说，它在某个固定基ψ下可以精确地写成少量的基向量的叠加。这个前提从根本上改变了问题，使得求解压缩感知理论中欠定方程组问题可解，而观测矩阵具有RIP性质也为从K个观测值中精确恢复信号提供了理论保证。目前，存在很多压缩感知中重构原始稀疏信号的方法，例如：GPSR、state-of-the-art、IST、l1_ls和OMP等方法。

IST最初应用于反卷积方面，对于解决压缩感知问题，观测矩阵的行数一般小于列数，IST并不有效。l1_ls是基于共轭梯度下降算法提出的一种方法，每一次循环都需要寻找下降方向，浪费了时间。OMP方法虽然节约了大量时间，但是对于压缩感知中不一致情况，它对噪声的鲁棒性太低。GPSR和state-of-the-art当循环次数多时，浪费了时间。

发明内容

本发明的目的在于克服上述技术的不足，提出一种稀疏信号重构方法，以提高重构原始稀疏信号的速率和准确度。

为实现上述目的，本发明包括如下步骤：

(1)利用计算机得到原始稀疏信号f∈R^N、观测矩阵A∈R^K×N和观测向量y∈R^K；

(2)从含有K个元素的观测向量y中，按如下步骤重构含有N个元素的原始稀疏信号f：

2a)设定重构原始稀疏信号f的无约束目标函数：

$\begin{array}{cc}\min & \frac{1}{2}{\left|\right|y-\mathrm{Ax}\left|\right|}_2^2+\mathrm{\ϵ}{\left|\right|x\left|\right|}_1---\left(1A\right)\end{array}$

式中，x∈R^N为重构信号，ε为惩罚因子，通过调节ε的大小，控制重构原始稀疏信号f的准确度；

2b)将无约束目标函数中重构信号x的N个元素用β_i-β_i ^*，i＝1，...，N代替，得到约束目标函数，令 $u=\left[\begin{array}{c}{\mathrm{\β}}_i\\ {\mathrm{\β}}_i^{*}\end{array}\right],$ i＝1，...，N，其中μ∈R^2N，初始化向量μ的所有元素都为零；

2c)求取约束目标函数的Lagrange对偶函数；

2d)由约束目标函数和Lagrange对偶函数得到满足Karush-Kuhn-Tucker，即KKT条件的最优解；

2e)以KKT条件作为重构稀疏信号的终止条件判断，判断当前向量μ中所有元素是否满足KKT条件，若有违反KKT条件的元素，执行步骤2f)，否则跳至步骤2h)；

2f)设定一个工作集work，使得工作集work中元素个数为nwork，并且选用当前向量μ中最大违反KKT条件的前nwork个元素作为工作集work中元素；

2g)计算工作集work中元素的最优解，更新向量μ中与工作集work中元素相对应的元素；

2h)计算重构信号的当前值x_i＝u_i-u_i+N，i＝1，...，N；

2i)对重构信号进行除偏处理，得到有除偏处理的重构信号x^*；

(3)计算重构原始稀疏信号f的准确度，用 $\mathrm{\δ}=(1/N){\left|\right|x-f\left|\right|}_2^2$ 评价重构信号的无除偏处理准确度，其中符号δ表示重构信号的无除偏处理均方误差，用重构信号的有除偏处理均方误差 $m=(1/N){\left|\right|x^{*}-f\left|\right|}_2^2$ 评价重构信号的有除偏处理重构准确度。

本发明与其他技术相比具有以下优点：

1、本发明在无约束目标函数的基础上，推导出一种约束目标函数，并且通过调整参数ε使重构信号的准确度和重构速率达到预想结果。

2、本发明考虑到约束目标函数的可分解性，将约束目标函数分解为一系列小的线性方程组来解决，节约了存储空间。

3、本发明将序列最小优化方法应用到重构原始稀疏信号的约束目标函数上，并且使用KKT条件作为重构稀疏信号的终止条件判断，具有很好的收敛性和快速性。

4、本发明在重构原始稀疏信号的过程中，通过人工设定当前工作集中元素的个数，控制重构原始稀疏信号过程的循环次数，提高了重构速率。

附图说明

图1是本发明对稀疏信号重构的过程图；

图2是用本发明方法和三种GPSR方法进行目标函数值随时间的仿真曲线图；

图3是用本发明方法和一种GPSR方法进行目标函数值和均方误差值随时间的仿真曲线图；

图4是用本发明方法、GPSR方法和两种OMP方法进行均方误差值和时间随原始稀疏信号中包含非零元素个数的仿真曲线图。

具体实施方式

参照图1，本发明包括原始数据获取过程、重构原始稀疏信号过程和准确度评价过程。具体实现如下：

步骤1：利用计算机得到原始稀疏信号f∈R⁴⁰⁹⁶、观测矩阵A∈R^1024×4096和观测向量y∈R¹⁰²⁴。

(1a)利用计算机产生包含θ个±1元素的原始稀疏信号f∈R⁴⁰⁹⁶，其他4096-θ个元素为零；

(1b)利用计算机产生满足均值为0、方差为1和独立同分布的1024×4096高斯矩阵AA，将矩阵AA按行进行正交化处理，得到观测矩阵A∈R^1024×4096；

(1c)根据原始稀疏信号f和观测矩阵A，利用计算机按照等式y＝Af+n计算观测向量y∈R¹⁰²⁴，其中n是方差为σ²＝10^-4的高斯白噪声。

步骤2：从含有1024个元素的观测向量y中精确重构含有4096个元素的原始稀疏信号f。

(2a)设定重构原始稀疏信号f的无约束目标函数；

$\begin{array}{cc}\min & \frac{1}{2}{\left|\right|y-\mathrm{Ax}\left|\right|}_2^2+\mathrm{\ϵ}{\left|\right|x\left|\right|}_1---a)\end{array}$

式中，x∈R⁴⁰⁹⁶为重构信号，ε为惩罚因子，通过调节ε的大小，控制重构原始稀疏信号f的准确度；

(2b)将无约束目标函数中重构信号x的4096个元素用β_i-β_i ^*，i＝1，...，4096代替，初始化向量μ的所有元素都为零，得到约束目标函数；

$\begin{array}{cc}\underset{\mathrm{\μ}}{\min }& F_M=\frac{1}{2}{\mathrm{\μ}}^T{H}_b\mathrm{\μ}+c^T\mathrm{\μ}+\frac{1}{2}{y}^{T}y\end{array}---b)$

st  u_i≥0，i＝1，...，8192

式中， $H_b=\left[\begin{array}{cc}H& -H\\ -H& H\end{array}\right],$ H＝A^TA， $c=\left[\begin{array}{c}-y^{T}A+\mathrm{\ϵ}\\ -y^{T}A+\mathrm{\ϵ}\end{array}\right],$ $\mathrm{\μ}=\left[\begin{array}{c}{\mathrm{\β}}_i\\ {\mathrm{\β}}_i^{*}\end{array}\right],$ i＝1，...，4096；

(2c)求约束目标函数b)的Lagrange对偶函数；

$F_D=\frac{1}{2}{\mathrm{\μ}}^T{H}_b\mathrm{\μ}+c^T\mathrm{\μ}+\frac{1}{2}{y}^{T}y-{\mathrm{\ω}}^T\mathrm{\μ}---c)$

st. u_i≥0，ω_i≥0，i＝1，...，8192

式中ω为Lagrange乘子；

(2d)由约束目标函数b)和Lagrange对偶函数c)得到最优解满足的KKT条件(Karush-Kuhn-Tucker)；

d-ω＝0

ω_i≥0

             d)

μ_i≥0

μ_iω_i＝0，i＝1，...，8192

式中d＝H_bμ+c，为约束目标函数b)的梯度，d＝ω；将式d)分解为两个条件：如果μ_i＝0，则要求d_i≥0；如果μ_i＞0，则要求μ_i＝0；

(2e)以KKT条件作为重构稀疏信号算法的终止条件判断，判断当前向量μ中所有元素是否满足KKT条件，若有违反KKT条件的元素，继续步骤(2f)，否则跳至步骤(2h)；

(2f)设定一个工作集work，使得工作集work中元素个数为nwork，并且选用当前向量μ中最大违反KKT条件的前nwork个元素作为工作集work中元素；

(2g)根据H值，找到当前工作集work中元素所对应的二次项系数矩阵Q和常数项系数向量b，得到一个nwork元线性方程组，求解该线性方程组，得到工作集work中元素的最优解，更新向量μ中与该工作集work中元素相对应的元素；

(2h)计算重构信号当前值x_i＝u_i-u_i+N，i＝1，...，4096；

(2i)对重构向量进行除偏处理，得到有除偏处理的重构信号x^*。

步骤3：计算重构原始稀疏信号f的准确度。

用 $\mathrm{\δ}=(1/4096){\left|\right|x-f\left|\right|}_2^2$ 评价重构信号的无除偏处理准确度，其中符号δ表示重构信号的无除偏处理均方误差，用重构信号的有除偏处理均方误差 $m=(1/4096){\left|\right|x^{*}-f\left|\right|}_2^2$ 评价重构信号的有除偏处理重构准确度。

本发明的效果通过以下仿真进一步说明：

1、仿真条件与内容

对压缩感知领域中稀疏信号的重构按照上述具体实施过程进行仿真。在我们的实验中，考虑一种比较典型的压缩感知问题。目的是从含有1024个元素的观测向量y中精确重构含有4096个元素的原始稀疏信号f。

仿真1，利用计算机产生包含160个值为±1的原始稀疏信号，并且在满足等式y＝Af+n的条件下，取ε＝0.08，从含有1024个元素的观测向量y中精确重构含有4096个元素的原始稀疏信号f。

用本发明方法进行循环次数、时间和均方误差值随原始稀疏信号中包含非零元素个数的10次仿真，记录下每次仿真的结果，计算10次的平均结果，结果如表1所示。

表1.用本发明方法进行10次仿真的平均结果

非零元素个数	循环次数	时间	均方误差值
				1	374	2.43	0.007
2	213	1.48	0.007
				3	152	1.13	0.007
4	117	0.97	0.007
				5	97	0.91	0.007
6	83	0.83	0.0068
				7	74	0.88	0.0068
8	65	0.87	0.007
				9	59	0.91	0.0068
10	53	0.93	0.0067
				11	50	0.99	0.007
12	45	1.01	0.007
				13	42	1.07	0.007
14	41	1.15	0.007
				15	39	1.21	0.007
16	37	1.27	0.0067
				17	35	1.38	0.0067
18	33	1.43	0.007
				19	32	1.5	0.0067
20	30	1.55	0.0067

仿真2，利用计算机产生包含160个值为±1的原始稀疏信号，并且在满足等式y＝Af+n的条件下，取nwork＝6和ε＝0.08，从含有1024个元素的观测向量y中精确重构含有4096个元素的原始稀疏信号f。

用本发明方法和三种GPSR方法进行目标函数值随时间的仿真曲线图，如图2所示。

用本发明方法和一种GPSR方法进行目标函数值和均方误差值随时间的仿真曲线图，如图3所示。

用本发明方法、l1_ls、IST和三种GPSR方法进行10次仿真，计算平均结果，结果如表2所示。

表2.用本发明方法、l1_ls、IST和三种GPSR方法进行10次仿真的平均结果

方法	时间(秒)	均方误差值
			本发明方法	1.1	0.0068
GPSR_BB monotone	1.93	0.007
			GPSR_BB nonmonotone	1.86	0.0067
GPSR_Basic	2.05	0.007
			l1_ls	8.22	0.0067
IST	6.28	0.007

仿真3，利用计算机产生包含5-250个值为±1的原始稀疏信号，并且在满足等式y＝Af+n的条件下，从含有1024个元素的观测向量y中精确重构含有4096个元素的原始稀疏信号f。

用本发明方法、GPSR方法和两种OMP方法进行均方误差值和时间随原始稀疏信号中包含非零元素个数的仿真曲线图，如图4所示。

2、仿真结果分析

从表1所示的循环次数、时间和均方误差值随原始稀疏信号中包含非零元素个数的变化情况可以看出，当nwork＝6时，本发明方法所需时间最低。

从图2中可以看出，与三种GPSR方法相比，本发明方法的目标函数值随时间下降比较快。

图3(a)给出了用本发明方法和一种GPSR方法进行目标函数值随时间的仿真曲线图，图3(b)给出了用本发明方法和一种GPSR方法进行均方误差值随时间的仿真曲线图。从图3(a)中可以看出有除偏处理重构方法的目标函数值比无除偏处理重构方法的目标函数值大，是由有除偏处理过程与无除偏处理过程优化的目标函数不同所致。

从表2中可以看出，与l1_ls、IST和三种GPSR方法相比，在保证准确度的情况下，本发明方法的重构速率比较快。

图4(a)给出了用本发明方法、GPSR方法和两种OMP方法进行均方误差值随原始稀疏信号中包含非零元素个数的仿真曲线图。图4(b)给出了用本发明方法、GPSR方法和两种OMP方法进行时间随原始稀疏信号中包含非零元素个数的仿真曲线图。从图4可以看出，与GPSR方法和两种OMP方法相比，在保证准确度的情况下，本发明方法的重构速率比较快。

上述结果表明，本发明与三种GPSR方法、IST、l1_ls和两种OMP方法相比，具有明显的优点。首先，它只有一个参数，并且易于调节；其次，它的重构速率高于其他方法；最后，它的重构准确度并没有低于其他方法。

Claims (7)

1、一种稀疏信号的重构方法，包括如下步骤：

(1)利用计算机得到原始稀疏信号f∈R^N、观测矩阵A∈R^K×N和观测向量y∈R^K；

(2)从含有K个元素的观测向量y中，按如下步骤重构含有N个元素的原始稀疏信号f：

2a)设定重构原始稀疏信号f的无约束目标函数：

$\begin{array}{cc}\min & \frac{1}{2}{\left|\right|y-\mathrm{Ax}\left|\right|}_2^2+\mathrm{\ϵ}{\left|\right|x\left|\right|}_1\end{array}---\left(1A\right)$

式中，x∈R^N为重构信号，ε为惩罚因子，通过调节ε的大小，控制重构原始稀疏信号f的准确度；

2b)将无约束目标函数中重构信号x的N个元素用β_i-β_i ^*，i＝1，...，N代替，得到约束目标函数，令 $\mathrm{\μ}=\left[\begin{array}{c}{\mathrm{\β}}_i\\ {\mathrm{\β}}_i^{*}\end{array}\right],$ i＝1，...，N，其中μ∈R^2N，初始化向量μ的所有元素都为零；

2c)求取约束目标函数的Lagrange对偶函数；

2d)由约束目标函数和Lagrange对偶函数得到满足Karush-Kuhn-Tucker，即KKT条件的最优解；

2e)以KKT条件作为重构稀疏信号的终止条件判断，判断当前向量μ中所有元素是否满足KKT条件，若有违反KKT条件的元素，执行步骤2f)，否则跳至步骤2h)；

2f)设定一个工作集work，使得工作集work中元素个数为nwork，并且选用当前向量μ中最大违反KKT条件的前nwork个元素作为工作集work中元素；

2g)计算工作集work中元素的最优解，更新向量μ中与工作集work中元素相对应的元素；

2h)计算重构信号的当前值x_i＝u_i-u_i+N，i＝1，...，N；

2i)对重构信号进行除偏处理，得到有除偏处理的重构信号x^*；

(3)计算重构原始稀疏信号f的准确度，用 $\mathrm{\δ}=(1/N){\left|\right|x-f\left|\right|}_2^2$ 评价重构信号的无除偏处理准确度，其中符号δ表示重构信号的无除偏处理均方误差，用重构信号的有除偏处理均方误差 $m=(1/N){\left|\right|x^{*}-f\left|\right|}_2^2$ 评价重构信号的有除偏处理重构准确度。

2、根据权利要求1所述的重构方法，其中步骤(1)所述的利用计算机得到原始稀疏信号、观测矩阵和观测向量，包括如下步骤：

(a)利用计算机产生包含θ个非零元素的原始稀疏信号f，其他N-θ个元素值为零；

(b)利用计算机产生满足独立同分布、均值为0和方差为1的K×N高斯矩阵AA，将矩阵AA按行进行正交化处理，得到观测矩阵A∈R^K×N；

(c)根据原始稀疏信号f和观测矩阵A，利用计算机按照等式y＝Af，或者按照等式y＝Af+n计算观测向量y∈R^K，其中K＜N，n是方差为σ²的高斯白噪声。

3、根据权利要求1所述的重构方法，其中步骤2b)中的约束目标函数，用矩阵表示如下：

$\begin{array}{cc}\underset{\mathrm{\μ}}{\min }& F_M=\frac{1}{2}{\mathrm{\μ}}^T{H}_b\mathrm{\μ}+c^T\mathrm{\μ}+\frac{1}{2}{y}^{T}y\end{array}---\left(3A\right)$

st  u_i≥0，i＝1，...，2N

式中， $H_b=\left[\begin{array}{cc}H& -H\\ -H& H\end{array}\right],$ H＝A^TA， $c=\left[\begin{array}{c}-y^{T}A+\mathrm{\ϵ}\\ -y^{T}A+\mathrm{\ϵ}\end{array}\right].$

4、根据权利要求1所述的重构方法，其中步骤2c)中的Lagrange对偶函数，用矩阵表示如下：

$F_D=\frac{1}{2}{\mathrm{\μ}}^T{H}_b\mathrm{\μ}+c^T\mathrm{\μ}+\frac{1}{2}{y}^{T}y-{\mathrm{\ω}}^T\mathrm{\μ}---\left(4A\right)$

st.u_i≥0，ω_i≥0，i＝1，...，2N

式中ω为Lagrange乘子。

5、根据权利要求1所述的重构方法，其中步骤2d)中的KKT条件为：

d-ω＝0

ω_i≥0    (5A)

μ_i≥0

μ_iω_i＝0，i＝1，...，2N

式中d＝H_bμ+c，为目标函数(3A)的梯度，d＝ω；将式(5A)分解为两个条件：如果μ_i＝0，则要求d_i≥0；如果μ_i＞0，则要求μ_i＝0。

6、根据权利要求1所述的重构方法，其中步骤2g)中所述的计算该工作集work中元素的最优解，是根据H值，找到当前工作集work中元素所对应的二次项系数矩阵Q和常数项系数向量b，得到一个nwork元线性方程组，求解该线性方程组，得到工作集work中元素的最优解。

7、根据权利要求1所述的重构方法，其中步骤2i)中所述的对重构向量进行除偏处理，是采用共轭梯度下降算法对无约束目标函数(1A)的二次项||y-Ax||₂ ²进行优化，得到有除偏处理的重构信号x^*。

Images