CN104767535A - 一种低复杂度块状稀疏信号重构方法 - Google Patents

Abstract

本发明属于稀疏信号恢复的技术领域，尤其涉及一种在信号接收中基于广义近似消息传递(Generalized Approximate Message Passing,GAMP)的低复杂度块状稀疏信号重构方法。一种低复杂度块状稀疏信号重构方法，包括如下步骤：通过压缩感知采样得到接收信号y＝A×x；初始化；进行循环迭代；输入恢复信号。与传统的基于贝叶斯的恢复算法相比，本发明方法利用广义近似消息传递方法，本发明可以在保持信号优越重构性能的基础上，有效地降低计算复杂度。将计算时间复杂度从O(n³)降为O(mn)，即从三次方复杂度降为线性时间复杂度，极大的缓解后端的信号处理压力，其中，m和n分别为观测值的维度和原始信号的维度。

Description

一种低复杂度块状稀疏信号重构方法

技术领域

本发明属于稀疏信号恢复的技术领域，尤其涉及一种在信号接收中基于广义近似消息传递(Generalized Approximate Message Passing,GAMP)的低复杂度块状稀疏信号重构方法。

背景技术

在过去的几年中，随着信号带宽的不断增加，在射频系统中，数字化信号对模数转换器(Analog to Digital Converter，AD)的要求越来越高。AD转换速率越高，功耗越大，有效位数也会降低。最近，压缩感知技术(Compressed Sensing)作为一种低速率采集稀疏信号的新技术在学术界和工业界被广泛研究与应用，它的应用背景是信号具有稀疏性，即信号只有很少的非零稀疏。例如，自然图像和通信信号在一些变换域上是稀疏分布的，或者有些图像本身就在时域上稀疏，这样的稀疏信号可以通过远低于奈奎斯特采样率直接进行采样。原始信号可以通过贪婪算法，线性规划或者基于贝叶斯推断等方法恢复信号。目前，基于贝叶斯的恢复算法具有相对优越的恢复性能，然而，基于贝叶斯的恢复算法一般包含矩阵求逆过程，在大规模的信号重构问题中并不具有实用价值。

发明内容

本发明为了解决现有技术中存在的不足，提供了一种低复杂度块状稀疏信号重构方法，该方法采用广义近似消息传递(GAMP)方法，将计算时间从三次方复杂度降为线性时间复杂度。

为了方便地描述本发明的内容，首先对本发明中使用的术语进行定义。

块状稀疏信号：信号中只有少数几个分量的值为非零，则称该信号是稀疏的，当这些分量呈块状聚集，则称该信号为块状稀疏信号。

期望最大算法(Expectation Maximization,EM)：不断建立似然函数的下界，并对下界进行优化，进而最大化似然函数，是一种解决存在隐含变量优化问题的有效方法。

一种低复杂度块状稀疏信号重构方法，具体步骤如下：

S1、通过压缩感知采样得到接收信号y＝A×x，其中，A为维度为m×n的投影矩阵，x为待恢复的块状稀疏信号，x＝[x₁,…,x_n]^T，y＝[y₁,…,y_m]^T，m＜＜n，m和n为大于1的自然数，[*]^T表示转置；

S2、构造函数z＝Ax，采用广义近似消息传递方法求解S1所述待恢复的块状稀疏信号x的后验均值，根据所述后验均值进行循环迭代求解，直到循环次数k等于系统预先设定的循环次数K，或前后两次迭代过程所得估计数据差值小于容许误差ε为止，其中，k为自然数，k≤K，K为大于等于1的自然数，具体步骤如下：

S21、初始化，

S1所述x的初始均值为S1所述x对应初始估计方差为τ^(x)(0)＝1，S1所述x的对偶变量s的估计初始值为稀疏控制参数α的初始值为α⁰＝1，噪声初始估计方差为σ²(0)＝var(y)/C，其中，参数α为对S1所述x进行稀疏控制的参数，var(y)表示S1所述接收信号y的方差，C为常数,j∈[1,...,m]；

S22、当循环次数k＝1时，利用S21所述α⁰计算出当前的后验估计均值，通过期望最大化(Expectation Maximuzation，EM)算法对所述α⁰和σ²(0)进行更新，得到α¹和σ²(1)，其中，α¹表示循环次数为1的稀疏控制参数α的值，σ²(1)表示循环次数为1的噪声估计方差，为循环次数k＝1时的块状稀疏信号；

S23、当循环次数k＝k+1时，利用α^k计算出的后验估计均值，通过EM算法对α^k和σ²(k)进行更新，得到

${\mathrm{\α}}_i^{k+1}=\frac{1}{{\left({\hat{x}}_i^{k+1}\right)}^2+{\mathrm{\τ}}_i^x(k+1)+\mathrm{\β}[{\left({\hat{x}}_{i-1}^{k+1}\right)}^2+{\mathrm{\τ}}_{i-1}^x(k+1)]+\mathrm{\β}[{\left({\hat{x}}_{i+1}^{k+1}\right)}^2+{\mathrm{\τ}}_{i+1}^x(k+1)]+2b}$ 和

其中，α^k表示循环次数k＝k时的稀疏控制参数α，α^k+1表示循环次数k＝k+1时的稀疏控制参数α，σ²(k)表示循环次数k＝k时的噪声估计方差，σ²(k+1)表示循环次数k＝k+1时的噪声估计方差，为循环次数k＝k+1时的待恢复块状稀疏信号，i∈[1,…,n]，b＝1e-8，β＝1，表示z当前估计值的第j个分量，表示z当前估计方差的第j个分量；

S24、若满足下述任意条件则停止循环迭代，条件如下：

条件1、k＝K，则此时为恢复信号，

条件2、则此时为恢复信号；

S3、输出恢复的信号

进一步地，S2所述采用广义近似消息传递方法求解后验均值的方法具体如下：

A、对于m维向量中任意的j∈[1,…,m]， ${\mathrm{\τ}}_j^p\left(k\right)=\underset{i}{\mathrm{\Σ}}{a}_{\mathrm{ji}}^2{\mathrm{\τ}}_i^x\left(k\right),{\hat{p}}_j^k=\underset{j}{\mathrm{\Σ}}{a}_{\mathrm{ji}}{\hat{x}}_i^k-{\mathrm{\τ}}_j^p\left(k\right){\hat{s}}_j^{k-1},$ 其中，为p^k的估计值，p^k是S2所述z加噪声后的变量，即p^k＝z+w₁，w₁为噪声方差等于τ^p(k)的零均值高斯白噪声；

B、 ${\hat{s}}_j^k=(y_j-{\hat{p}}_j^k)/({\mathrm{\τ}}_j^p\left(k\right)+{\mathrm{\σ}}^2\left(k\right)),{\mathrm{\τ}}_j^s\left(k\right)=1/({\mathrm{\τ}}_j^p\left(k\right)+{\mathrm{\σ}}^2\left(k\right)),$ 其中，s是x的对偶变量，为该变量的当前估计值，τ^s(k)为对应精确度；

C、对于n维向量中任意的i∈[1,…,n]， ${\mathrm{\τ}}_i^r\left(k\right)={\left(\underset{j}{\mathrm{\Σ}}{a}_{\mathrm{ji}}^2{\mathrm{\τ}}_j^s\left(k\right)\right)}^{-1},{\hat{r}}_i^k={\hat{x}}_i^k+{\mathrm{\τ}}_i^r\left(k\right)\underset{j}{\mathrm{\Σ}}{a}_{\mathrm{ji}}{\hat{s}}_j^k,$ 其中为r^k的估计值，r^k可以看作是x加噪声后的变量，即r^k＝x+w₂，w₂为零均值的高斯白噪声，噪声方差为τ^r(k)

D、 ${\hat{x}}_i^{k+1}=\frac{{\hat{r}}_i^k/{\mathrm{\τ}}_i^r\left(k\right)}{({\mathrm{\α}}_i^k+\mathrm{\β}{\mathrm{\α}}_{i-1}^k+\mathrm{\β}{\mathrm{\α}}_{i+1}^k)+1/{\mathrm{\τ}}_i^r\left(k\right)},{\mathrm{\τ}}_i^x(k+1)=1/(({\mathrm{\α}}_i^k+\mathrm{\β}{\mathrm{\α}}_{i-1}^k+\mathrm{\β}{\mathrm{\α}}_{i+1}^k)+1/{\mathrm{\τ}}_i^r\left(k\right)),$ 其中是x的当前估计值，τ^x(k+1)为对应方差。

进一步地，S2所述K＝300，ε＝1e-6。

本发明的有益效果是：

与传统的基于贝叶斯的恢复算法相比，本发明方法利用广义近似消息传递方法，可以在保持信号优越重构性能的基础上，有效地降低计算复杂度。将计算时间复杂度从O(n³)降为O(mn)，即从三次方复杂度降为线性时间复杂度，极大的缓解后端的信号处理压力，其中，m和n分别为观测值的维度和原始信号的维度。

本发明在恢复块状稀疏信号时，相对于高计算复杂度的贝叶斯重构方法几乎没有性能损失，可以保持与高计算复杂度的贝叶斯重构方法同等的恢复性能。

本发明计算复杂度低，更适用于实际信号恢复(比如常见的图片一般包含上万个像素点)。

在能量以及计算资源高度受限制的无线传感器网络中可以应用本发明作信号的采集和重构算法，有效地节省能量。同时，本发明可以满足对实时性要求高的视频监控等应用需求。

附图说明

图1是使用本发明方法对基于压缩感知的块状稀疏信号数据采集与重构方法流程图。

图2是在稀疏度不变的情况下重构成功率随观测值变化的性能比较图。

图3是在稀疏度不变的情况下各个恢复算法的平均运行时间比较图。

图4是各算法恢复性能随观测值变化的比较图(卫星图像恢复(256×256))。

具体实施方式

下面结合实施例和附图，详细说明本发明的技术方案。

本发明采用的模型：

在信号采集过程中基于压缩感知技术对信号进行采样，该场景可以是无线宽带信号接收端，或者是视频监控等场景中视频采集端等，模型可以描述为y＝Ax，y为在信号接收端一次采样时隙得到的数据(在视频监控等图像处理场景中x为矩阵，

可以通过向量化操作转换为向量形式，同时y和A作相应变形)，y是一个m维向量，x是一个n维向量，并且m＜＜n，其中x＝[x₁,…,x_n]^T，y＝[y₁,…,y_m]^T，矩阵A在硬件实现上可以通过PN伪随机码实现。

如图1所示，一种低复杂度块状稀疏信号重构方法，包括如下步骤：

步骤1、在无线宽带信号接收端，通过压缩感知采样得到接收信号y＝A×x，其中，A为投影矩阵，所述A的维度为m×n，x为将要恢复的块状稀疏信号，x＝[x₁,…,x_n]^T，y＝[y₁,…,y_m]^T，

矩阵A有PN伪随机码产生，测量信号y可以由量化后的数字信号与PN伪随机码矩阵A相乘实现，这一过程我们实现了原始信号的采样，由于m＜＜n，即测量量远小于原始信号长度。

步骤2、进行初始化，给定变量x的初始估计均值为x对应初始估计方差为τ^(x)(0)＝1，变量x的对偶变量s的估计初始值为参数α(α为对x进行稀疏控制的参数)的初始值为α⁰＝1，噪声初始估计方差为σ²(0)＝var(y)/C，其中，var(y)表示接收信号y的方差，C为一个常数。

步骤3、首先求解x的后验均值，在硬件实现上，可以采用数位讯号处理器(Digital Signal Processor，DSP)或者现场可编程门阵列(Field Programmable GateArray，FPGA)。

待恢复信号x进行循环迭代求解，主要分为两步：

第一步、在广义近似消息传递步骤，利用当前的稀疏控制参数估计值α^k计算出x，z＝Ax的后验估计均值；第二步、基于近似计算得到的x，z＝Ax后验估计均值，通过EM方法对超参数进行更新并获得α^k+1。首先采用广义近似消息传递方法，求解x的后验均值。

在广义近似消息传递方法中，

对于m维向量中任意的j∈[1,…,m]：

${\mathrm{\τ}}_j^p\left(k\right)=\underset{i}{\mathrm{\Σ}}{a}_{\mathrm{ji}}^2{\mathrm{\τ}}_i^x\left(k\right),{\hat{p}}_j^k=\underset{j}{\mathrm{\Σ}}{a}_{\mathrm{ji}}{\hat{x}}_i^k-{\mathrm{\τ}}_j^p\left(k\right){\hat{s}}_j^{k-1},$ 其中为p^k的估计值，p^k可以看作是z加噪声后的变量，即p^k＝z+w₁，w₁为零均值的高斯白噪声，噪声方差为τ^p(k)；

${\hat{s}}_j^k=(y_j-{\hat{p}}_j^k)/({\mathrm{\τ}}_j^p\left(k\right)+{\mathrm{\σ}}^2\left(k\right)),{\mathrm{\τ}}_j^s\left(k\right)=1/({\mathrm{\τ}}_j^p\left(k\right)+{\mathrm{\σ}}^2\left(k\right)),$ 其中s是x的对偶变量，为该变量的当前估计值，τ^s(k)为对应灵敏度；

于n维向量中任意的i∈[1,…,n]：

${\mathrm{\τ}}_i^r\left(k\right)={\left(\underset{j}{\mathrm{\Σ}}{a}_{\mathrm{ji}}^2{\mathrm{\τ}}_j^s\left(k\right)\right)}^{-1},{\hat{r}}_i^k={\hat{x}}_i^k+{\mathrm{\τ}}_i^r\left(k\right)\underset{j}{\mathrm{\Σ}}{a}_{\mathrm{ji}}{\hat{s}}_j^k,$ 其中为r^k的估计值，r^k可以看作是x加噪声后的变量，即r^k＝x+w₂，w₂为零均值的高斯白噪声，噪声方差为τ^r(k)；

${\hat{x}}_i^{k+1}=\frac{{\hat{r}}_i^k/{\mathrm{\τ}}_i^r\left(k\right)}{({\mathrm{\α}}_i^k+\mathrm{\β}{\mathrm{\α}}_{i-1}^k+\mathrm{\β}{\mathrm{\α}}_{i+1}^k)+1/{\mathrm{\τ}}_i^r\left(k\right)},{\mathrm{\τ}}_i^x(k+1)=1/(({\mathrm{\α}}_i^k+\mathrm{\β}{\mathrm{\α}}_{i-1}^k+\mathrm{\β}{\mathrm{\α}}_{i+1}^k)+1/{\mathrm{\τ}}_i^r\left(k\right)),$ 其中是x的当前估计值，τ^x(k+1)为对应方差，其中β＝1。

步骤4、然后利用期望最大算法(EM)算法更新噪声方差估计值σ²和x稀疏控制参数向量α的值，对于n维向量α中任意的i∈[1,…,n]： ${\mathrm{\α}}_i^{k+1}=\frac{1}{{\left({\hat{x}}_i^{k+1}\right)}^2+{\mathrm{\τ}}_i^x(k+1)+\mathrm{\β}[{\left({\hat{x}}_{i-1}^{k+1}\right)}^2+{\mathrm{\τ}}_{i-1}^x(k+1)]+\mathrm{\β}[{\left({\hat{x}}_{i+1}^{k+1}\right)}^2+{\mathrm{\τ}}_{i+1}^x(k+1)]+2b}$ 其中b＝1e-8。

步骤5、重复步骤3和步骤4，直到循环次数k等于系统预先设定的循环次数K＝300，或前后两次迭代过程所得估计数据差值小于容许误差ε为止，即ε＝1e-6。

步骤6、输出恢复的信号

经过上述操作，完成了对信号的采集与重构的过程。

下面对经典重构算法及本发明方法的算法性能对比分析，以进一步验证本发明的性能。

图2与图3表征了采用传统基追踪算法(Basis Pursuit)、SBL算法(SparseBayesian Learning)和PCSBL算法(Pattern Coupled Sparse Bayesian Learning)以及使用本发明采用的低复杂度信号恢复算法的性能比较，图2可以看出本发明的提出的算法在稀疏度情况相同下与高复杂度的PCSBL方法具有同等的恢复性能，即几乎相同的重构成功率(当小于10-6则表示该次信号重构实验成功，x为真实信号，为恢复信号)，同时由图3表明在相同条件下，本发明提出的算法具有最小的计算复杂度，单次恢复平均计算时间远小于其他算法。

图4表示了在计算复杂度相当的几类算法中，本发明的算法利用不同的观测值进行恢复时在m/n＜0.25时具有远小于其他算法的NMSE(Normalized Mean SquaredError)，也就是说只要利用很少的观测值，本发明的算法便可以实现对原始信号的可靠恢复。

Claims (3)

1.一种低复杂度块状稀疏信号重构方法，其特征在于，包括如下步骤：

S1、通过压缩感知采样得到接收信号y＝A×x，其中，A为维度为m×n的投影矩阵，x为待恢复的块状稀疏信号，x＝[x₁,...,x_n]^T，y＝[y₁,...,y_m]^T，m＜＜n，m和n为大于1的自然数，[*]^T表示转置；

S2、构造函数z＝Ax，采用广义近似消息传递方法求解S1所述待恢复的块状稀疏信号x的后验均值，根据所述后验均值进行循环迭代求解，直到循环次数k等于系统预先设定的循环次数K，或前后两次迭代过程所得估计数据差值小于容许误差ε为止，其中，k为自然数，k≤K，K为大于等于1的自然数，具体步骤如下：

S21、初始化，

S1所述x的初始均值为S1所述x对应初始估计方差为τ^(x)(0)＝1，S1所述x的对偶变量s的估计初始值为稀疏控制参数α的初始值为α⁰＝1，噪声初始估计方差为σ²(0)＝var(y)/C，其中，参数α为对S1所述x进行稀疏控制的参数，var(y)表示S1所述接收信号y的方差，C为常数,j∈[1,...,m]；

S22、当循环次数k＝1时，利用S21所述α⁰计算出当前的后验估计均值，通过期望最大化(Expectation Maximuzation，EM)算法对所述α⁰和σ²(0)进行更新，得到α¹和σ²(1)，其中，α¹表示循环次数为1的稀疏控制参数α的值，σ²(1)表示循环次数为1的噪声估计方差，为循环次数k＝1时的块状稀疏信号；

S23、当循环次数k＝k+1时，利用α^k计算出的后验估计均值，通过EM算法对α^k和σ²(k)进行更新，得到

${\mathrm{\α}}_i^{k+1}=\frac{1}{{\left({\hat{x}}_i^{k+1}\right)}^2{\mathrm{\τ}}_i^x(k+1)+\mathrm{\β}[{\left({\hat{x}}_{i-1}^{k+1}\right)}^2+{\mathrm{\τ}}_{i-1}^x(k+1)]+\mathrm{\β}[{\left({\hat{x}}_{i-1}^{k+1}\right)}^2+{\mathrm{\τ}}_{i+1}^x(k+1)]+2b}$ 和

其中，α^k表示循环次数k＝k时的稀疏控制参数α，α^k+1表示循环次数k＝k+1时的稀疏控制参数α，σ²(k)表示循环次数k＝k时的噪声估计方差，σ²(k+1)表示循环次数k＝k+1时的噪声估计方差，为循环次数k＝k+1时的待恢复块状稀疏信号，i∈[1,...,n]，b＝1e-8，β＝1，表示z当前估计值的第j个分量，表示z当前估计方差的第j个分量；

S24、若满足下述任意条件则停止循环迭代，条件如下：

条件1、k＝K，则此时为恢复信号，

条件2、 $\mathrm{norm}{\left|\right|{\hat{x}}^{k+1}{-\hat{x}}^k\left|\right|}_2<\mathrm{\&epsiv;},$ 则此时为恢复信号；

S3、输出恢复的信号

2.根据权利要求1所述一种低复杂度块状稀疏信号重构方法，其特征在于：S2所述采用广义近似消息传递方法求解后验均值的方法具体如下：

A、对于m维向量中任意的j∈[1,...,m]， ${\mathrm{\&tau;}}_j^p\left(k\right)=\underset{i}{\mathrm{\&Sigma;}}{a}_{\mathrm{ji}}^2{\mathrm{\&tau;}}_i^x\left(k\right),{\hat{p}}_j^k=\underset{i}{\mathrm{\&Sigma;}}{a}_{\mathrm{ji}}{\hat{x}}_i^k-{\mathrm{\&tau;}}_j^p\left(k\right){\hat{x}}_j^{k-1},$ 其中，为p^k的估计值，p^k是S2所述z加噪声后的变量，即p^k＝z+w₁，w₁为噪声方差等于τ^p(k)的零均值高斯白噪声；

B、 ${\hat{x}}_j^k=(y_j-{\hat{p}}_j^k)/({\mathrm{\&tau;}}_j^p\left(k\right)+{\mathrm{\&sigma;}}^2\left(k\right)),{\mathrm{\&tau;}}_j^s\left(k\right)=1/({\mathrm{\&tau;}}_j^p\left(k\right)+{\mathrm{\&sigma;}}^2\left(k\right)),$ 其中，s是x的对偶变量，为该变量的当前估计值，τ^s(k)为对应精确度；

C、对于n维向量中任意的i∈[1,...,n]， ${\mathrm{\&tau;}}_i^r\left(k\right)={\left(\underset{j}{\mathrm{\&Sigma;}}{a}_{\mathrm{ji}}^2{\mathrm{\&tau;}}_j^s\left(k\right)\right)}^{-1},{\hat{r}}_i^k={\hat{x}}_i^k+{\mathrm{\&tau;}}_i^r\left(k\right)\underset{j}{\mathrm{\&Sigma;}}{a}_{\mathrm{ji}}{\hat{s}}_j^k,$ 其中为r^k的估计值，r^k可以看作是x加噪声后的变量，即r^k＝x+w₂，w₂为零均值的高斯白噪声，噪声方差为τ^r(k)

D、 ${\hat{x}}_i^{k+1}=\frac{{\hat{r}}_i^k/{\mathrm{\&tau;}}_i^r\left(k\right)}{({\mathrm{\&alpha;}}_i^k+\mathrm{\&beta;}{\mathrm{\&alpha;}}_{i-1}^k+{\mathrm{\&beta;\&alpha;}}_{i+1}^k)+1/{\mathrm{\&tau;}}_i^r\left(k\right)},{\mathrm{\&tau;}}_i^x(k+1)=1/(({\mathrm{\&alpha;}}_i^k+{\mathrm{\&beta;\&alpha;}}_{i-1}^k+{\mathrm{\&beta;\&alpha;}}_{i+1}^k)+1/{\mathrm{\&tau;}}_i^r\left(k\right)),$ 其中是x的当前估计值，τ^x(k+1)为对应方差。

3.根据权利要求1所述一种低复杂度块状稀疏信号重构方法，其特征在于：S2所述K＝300，ε＝1e-6。