CN105375931B - 一种基于卡尔曼滤波的复杂环境下信号重构方法 - Google Patents

Abstract

本发明涉及一种基于卡尔曼滤波的复杂环境下信号重构方法，与现有技术相比解决了计算复杂度高、效率低的缺陷。本发明包括以下步骤：对信号进行快速压缩，设计稀疏测量矩阵，在复杂环境下进行压缩测量获得测量值；建立信号的先验模型，输入信号的稀疏率，建立信号的先验模型；在二分图上进行置信传播计算；采用近似MMSE估计得到信号估计的初始值；使用卡尔曼滤波得到信号估计值。本发明采用了简单的稀疏测量矩阵，简化了测量矩阵的存储，信号重构时结合二分图和基于卡尔曼滤波的信号估计方法，进一步简化压缩感知的编码过程并提高重构精度。

Description

一种基于卡尔曼滤波的复杂环境下信号重构方法

技术领域

本发明涉及信号处理技术领域，具体来说是一种基于卡尔曼滤波的复杂环境下信号重构方法。

背景技术

压缩感知是一种新型信号获取技术，能够以低于香农采样定理的低采样率获得信号的无失真重建。对于一些实际应用环境，如超宽带通信、医学成像、无线传感器网络系统以及雷达等应用中，一方面由于大的信号带宽会导致高速采样从而产生海量数据，造成存储和通信的巨大压力；另一方面由于应用环境的复杂性导致采集的信号中含有大量的噪声，造成信号恢复的困难。压缩感知技术为解决上述问题提供了一个好的思路，即若信号在某个变换基或字典上稀疏，那么可利用一个测量矩阵将其映射到一个低维空间，这大大降低了采样频率，随后通过处理一个信号重构问题，就能够从这些少量的低维信号高概率地精确重构原始信号。

应用环境的复杂性导致噪声普遍存在，而在压缩感知中，抑制噪声影响的常用方法就是在信号重构时，将噪声的干扰考虑进去，求解一个约束限制的优化问题。传统方法有两类：一类采用l₁范数法求解噪声干扰下的信号重构；还有一类基于概率的稀疏信号重构算法，能够降低噪声对压缩感知的影响。如：稀疏贝叶斯学习(Sparse Bayesian Learning，SBL)算法、贝叶斯压缩感知(Bayesian Compressive Sensing，BCS)、近似消息传播(Approximate Message Passing，AMP)算法。这些方法可以在一定程度上抑制噪声干扰，但以上方法测量矩阵通常采用密集高斯矩阵，而在实际应用环境中，由于存储器的存储能力大多有限，导致这些方法实际应用价值不高、范围有限。如何开发出一种能够将测量矩阵简化，从而能够简单、高效地重构信号已经成为急需解决的技术问题。

发明内容

本发明的目的是为了解决现有技术中计算复杂度高、效率低的缺陷，提供一种基于卡尔曼滤波的复杂环境下信号重构方法来解决上述问题。

为了实现上述目的，本发明的技术方案如下：

一种基于卡尔曼滤波的复杂环境下信号重构方法，包括以下步骤：

对信号进行快速压缩，设计稀疏测量矩阵Φ_M×N，在复杂环境下进行压缩测量获得测量值y；

建立信号的先验模型，输入信号x的稀疏率q，建立信号x的先验模型f(x)；

在二分图上进行置信传播计算，定义变量节点b和校验节点c，建立二分图，以信号的先验为初始值进行迭代置信传播计算，得到信号的边缘分布f(v)；

采用近似MMSE估计得到信号估计的初始值

使用卡尔曼滤波得到信号估计值

所述的对信号进行快速压缩包括以下步骤：

定义稀疏信号x的维数N，压缩过以后的维数为M，计算压缩比ρ，其计算公式如下：

设定类LDPC矩阵的行权重r和列权重l，且

对类LDPC矩阵根据测量矩阵的密度或随机产生测量矩阵Φ_M×N中非零元素的位置向量Υ，令Φ_M×N中非零元素值交替定义为1和-1；

进行压缩测量，获得测量值y，其计算公式如下：

y＝Φx。

所述的建立信号的先验模型包括以下步骤：

定义信号x中的元素取x_i≠0时，使用高斯分布表示其概率分布，概率分布为

信号x中的元素取x_i＝0时，使用Dirac分布来近似概率分布，概率分布为δ(x)；

建立信号的先验模型f(x)，

其中：q:＝P_r{x_i≠0},i∈[1,2,...N]。

所述的在二分图上进行置信传播计算包括以下步骤：

设校验矩阵H的Tanner图G＝{(V,E)}，V为节点的集合，包含变量节点集合V_b和校验节点V_c，

V_b＝(b₁,b₂,……,b_n)，b_n为变量节点，其与校验矩阵H的各列相对应；V_c＝(c₁,c₂,……,c_m)，c_m表示校验节点，其与校验矩阵H的各行相对应；

E为不同类节点之间相连的边的集合，

建立测量矩阵Φ对应的二分图，二分图中每一条边连接变量节点X和测量值对应的校验节点Y，并且每一条边对应测量矩阵中的一个非零元素Φ_ij；

进行迭代置信传播计算；

设从变量节点到校验节点的消息编码为信号分量的后验概率概率密度，用μ_i→j表示，从校验节点到变量节点的消息编码为测量分量的概率密度，用μ_j→i表示；

从变量节点到校验节点的消息编码计算公式如下：

从校验节点到变量节点的消息编码计算公式如下：

对从变量节点到校验节点的消息编码计算公式和从校验节点到变量节点的消息编码计算公式进行迭代计算，直到消息值不再发生变化，迭代结束，输出信号值x的后验概率密度，如下所示：

所述的使用卡尔曼滤波得到信号估计值包括以下步骤：

令为输入卡尔曼滤波的初始值，

P(0)＝αI；

计算未经校正的变量估计误差的均方值P'(k)，其计算公式如下：

P'(k)＝AP(k-1)A^T，A是变量的增益矩阵，为常数，符号^T表示转置；

计算滤波增益矩阵H(k)，其计算公式如下：

其中为压缩测量时产生的噪声方差；

计算信号估计值其计算公式如下：

计算最小均方误差阵P(k)，其计算公式如下：

P(k)＝(I-H(k)Φ)P'(k)；

若||H(k)||₂＞ξ，ξ为常数，重复计算均方值P'(k)、滤波增益矩阵H(k)、信号估计值和最小均方误差阵P(k)；

若||H(k)||₂＜ξ，令k＝k+1，输出

选择中最大K个系数的位置作为支撑集Γ，其中K为稀疏信号中非零元素的个数，令：

还包括对信号估计值进行迭代更新，迭代执行在二分图上进行置信传播计算；采用近似MMSE估计得到信号估计的初始值使用卡尔曼滤波得到信号估计值三个步骤，直到重构精度满足误差要求或达到事先设定的最大迭代次数，输出更新后的信号估计值

所述的进行压缩测量，获得测量值y的计算公式如下：

y＝Φx+n。

有益效果

本发明的一种基于卡尔曼滤波的复杂环境下信号重构方法，与现有技术相比采用了简单的稀疏测量矩阵，简化了测量矩阵的存储，信号重构时结合二分图和基于卡尔曼滤波的信号估计方法，进一步简化压缩感知的编码过程并提高重构精度。其中测量矩阵为稀疏的类LDPC矩阵，编码速度快、节约压缩端的存储器；通过以置信传播计算获得的估计值作为初始值输入卡尔曼滤波估计信号，进一步消除噪声影响，且原始信号的分布状态不会影响重构结果，重构精度高同时具有普适性。

附图说明

图1为本发明的顺序流程图；

图2为本发明中校验矩阵H的Tanner图；

图3为本发明中测量矩阵Φ；

图4为图3的二分图。

具体实施方式

为使对本发明的结构特征及所达成的功效有更进一步的了解与认识，用以较佳的实施例及附图配合详细的说明，说明如下：

在日常复杂环境下信号测量值受噪声污染严重，如噪声为高斯噪声，服从均值为零，方差为本发明所述的一种基于卡尔曼滤波的复杂环境下信号重构方法，能够进一步简化压缩感知的编码过程并提高重构精度，如图1所示，其包括以下步骤：

第一步，对信号进行快速压缩。设计稀疏测量矩阵Φ_M×N，在复杂环境下进行压缩测量获得测量值y。首先需要设计一个合适的压缩测量矩阵，在这里采用一种类LDPC码，通过一个生成矩阵将信号序列映射为码字序列，生成矩阵采用一个奇偶校验阵H_M×N，压缩测量矩阵Φ_M×N等效为类LDPC码的校验矩阵H_M×N，原始稀疏信号x等效为信息码字C，整个压缩测量过程可等效为一个类LDPC的编码过程。由于校验矩阵H_M×N具有稀疏性，因此测量矩阵Φ_M×N为一个稀疏矩阵，相比于传统的高斯密集测量矩阵则大大节省了存储空间。其具体包括以下步骤：

(1)定义稀疏信号x的维数N，压缩过以后的维数为M，计算压缩比ρ，压缩比ρ可根据人为需要设定，根据压缩比ρ可设定类LDPC矩阵的行权重r和列权重l，压缩比ρ的计算公式如下：

(2)设定类LDPC矩阵的行权重r和列权重l，且

(3)对类LDPC矩阵根据测量矩阵的密度或随机产生测量矩阵Φ_M×N中非零元素的位置向量Υ，令Φ_M×N中非零元素值交替定义为1和-1。在此，对于需要压缩的信号x设计稀疏测量矩阵Φ_M×N，采用非规则码，此外为提高压缩感知系统的性能，我们在应用时允许负值元素出现在矩阵H中，即基于类LDPC编码的压缩测量矩阵元素为{0,1,-1}，且1和-1交替出现。

(4)进行压缩测量，获得测量值y，其计算公式如下：

y＝Φx。

由于压缩测量矩阵Φ_M×N等效为类LDPC码的校验矩阵H_M×N，原始稀疏信号x等效为信息码字C，整个压缩测量过程可等效为一个LDPC的编码过程，即y＝Φx＝H_M×NC。考虑到复杂环境下测量过程不可避免的引入噪声，因此在测量值y的计算公式可以为y＝Φx+n＝H_{M× N}C+n。

第二步，建立信号的先验模型。输入信号x的稀疏率q，建立信号x的先验模型f(x)。通过先验模型，可以建立信号变量x到其对应元素{x_i}的概率密度函数f(x_i)的映射。相对于密集高斯矩阵，稀疏测量矩阵Φ在进行压缩测量时能够获得的关于原始稀疏信号信号变量x的信息要更少，只要满足稀疏率那么测量次数取M＝Ο(K log(N))足够恢复出原始信号。其具体包括以下步骤：

(1)无论信号x的真实分布如何，定义信号x中的元素取x_i≠0时，使用高斯分布表示其概率分布，概率分布为

信号x中的元素取x_i＝0时，使用Dirac分布来近似概率分布，概率分布为δ(x)。

(2)建立信号的先验模型f(x)，

其中：q:＝P_r{x_i≠0},i∈[1,2,...N]。

第三步，在二分图上进行置信传播计算，定义变量节点b和校验节点c，建立二分图，以信号的先验为初始值进行迭代置信传播计算，得到信号的边缘分布f(v)。信号重构时在二分图上进行置信传播计算，可以得到一个基于最小均方误差的近似估计量，以此作为初始值进行卡尔曼滤波进一步消减初始估计的误差值，消除了噪声的影响，提高了重构精度。其具体包括以下步骤：

(1)如图2所示，按现有技术的方法设校验矩阵H的Tanner图G＝{(V,E)}，V为节点的集合，包含变量节点集合V_b和校验节点V_c，

V_b＝(b₁,b₂,……,b_n)，b_n为变量节点，其与校验矩阵H的各列相对应；V_c＝(c₁,c₂,……,c_m)，c_m表示校验节点，其与校验矩阵H的各行相对应；

E为不同类节点之间相连的边的集合，即两个变量节点或者两个校验节点之间不存在相连的边。也就是说在校验矩阵H中的任一个非零元素h_ij，在图2的Tanner图的映射中，第i个校验节点与第j个变量节点存在一条相连的边，则可以定义节点的度为与某个节点相连的边的总数。

(2)建立测量矩阵Φ对应的二分图，二分图中每一条边连接变量节点X和测量值对应的校验节点Y，并且每一条边对应测量矩阵中的一个非零元素Φ_ij。如图3所示，在此以测量矩阵Φ为例，按图2的Tanner图规则演变出图4所示的图3的二分图。在图3中可以看到，测量矩阵Φ的第一列，在第1位和第3位存在非零元素Φ_ij，在图4中，X的第一个点与Y的第一个点、第三个点连线；同理，测量矩阵Φ的第三列，在第2位和第4位存在非零元素Φ_ij，在图4中，X的第三个点与Y的第二个点、第四个点连线，以此产生测量矩阵Φ对应的二分图。

(3)进行迭代置信传播计算。由于需要在信号分量x对应的变量节点和测量值分量y对应的校验节点间双向传递消息，这里我们设从变量节点到校验节点的消息编码为信号分量的后验概率概率密度，用μ_i→j表示，从校验节点到变量节点的消息编码为测量分量的概率密度，用μ_j→i表示；

从变量节点到校验节点的消息编码计算公式如下：

符号η[]表示对归一化操作，因为置信传播迭代传递的消息是概率密度，因此对于迭代过程中传递的每一个消息都需要进行归一化处理。

从校验节点到变量节点的消息编码计算公式如下：

(4)对从变量节点到校验节点的消息编码计算公式和从校验节点到变量节点的消息编码计算公式进行迭代计算，直到消息值不再发生变化，迭代结束，输出信号值x的后验概率密度，如下所示：

第四步，按现有技术方法采用近似MMSE估计得到信号估计的初始值求的运算等效为求x的条件均值。

第五步，使用卡尔曼滤波得到信号估计值以作为输入的初始值进行卡尔曼滤波得到信号估计值其具体步骤如下：

(1)令为输入卡尔曼滤波的初始值，

P(0)＝αI，其中α可取一个很大的正数，I为一单位阵。

(2)计算未经校正的变量估计误差的均方值P'(k)，其计算公式如下：

P'(k)＝AP(k-1)A^T，A是变量的增益矩阵，为常数，符号^T表示转置。

(3)计算滤波增益矩阵H(k)，其计算公式如下：

其中为压缩测量时产生的噪声方差。

(4)计算信号估计值其计算公式如下：

(5)计算最小均方误差阵P(k)，其计算公式如下：

P(k)＝(I-H(k)Φ)P'(k)。

(6)若||H(k)||₂＞ξ，ξ为常数,即迭代次数，令k＝k+1，重复计算均方值P'(k)、滤波增益矩阵H(k)、信号估计值和最小均方误差阵P(k)；

若||H(k)||₂＜ξ，输出

(7)选择中最大K个系数的位置作为支撑集Γ，其中K为稀疏信号中非零元素的个数，令：

在此就可对信号x进行重构，但为了达到重构精度、满足误差要求，还可以对信号估计值进行迭代更新。

第六步，迭代执行在二分图上进行置信传播计算；采用近似MMSE估计得到信号估计的初始值使用卡尔曼滤波得到信号估计值三个步骤，直到重构精度满足误差要求或达到事先设定的最大迭代次数，输出更新后的信号估计值

在实际应用中，假设有一个一维的高斯信号x，其稀疏率为q，非零系数服从高斯分布，其应用面向复杂环境下的压缩和重构，建立测量值如y＝Φx+n，要从y中重构出按照传统的重构方法，的输出即为重构结果，可以看到其中对噪声的抑制并不完全，直接利用边缘分布采用条件均值进行估计信号误差较大，精度较差。本发明中利用的输入作为卡尔曼滤波的初始值，通过卡尔曼滤波进一步抑制噪声，提高重构精度。

再如，假设有一个二维的稀疏图像信号，设定小系数阈值为s，计算图像信号的稀疏率q，建立测量值y＝Φx+n，要从y中重构出建立先验模型，通过迭代置信传播计算得到以此作为输入卡尔曼滤波的初始值，通过卡尔曼滤波进一步抑制噪声，提高重构精度。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明的范围内。本发明要求的保护范围由所附的权利要求书及其等同物界定。

Claims (6)

1.一种基于卡尔曼滤波的复杂环境下信号重构方法，其特征在于，包括以下步骤：

11)对信号进行快速压缩，设计稀疏测量矩阵Φ_M×N，在复杂环境下进行压缩测量获得测量值y；

12)建立信号的先验模型，输入信号x的稀疏率q，建立信号x的先验模型f(x)；

13)在二分图上进行置信传播计算，定义变量节点b和校验节点c，建立二分图，以信号的先验为初始值进行迭代置信传播计算，得到信号的边缘分布f(v)；

14)采用近似MMSE估计得到信号估计的初始值

15)使用卡尔曼滤波得到信号估计值所述的使用卡尔曼滤波得到信号估计值包括以下步骤：

151)令为输入卡尔曼滤波的初始值，

152)计算未经校正的变量估计误差的均方值P'(k)，其计算公式如下：

P'(k)＝AP(k-1)A^T，A是变量的增益矩阵，为常数，符号T表示转置；

153)计算滤波增益矩阵H(k)，其计算公式如下：

其中为压缩测量时产生的噪声方差；

154)计算信号估计值其计算公式如下：

155)计算最小均方误差阵P(k)，其计算公式如下：

P(k)＝(I-H(k)Φ)P'(k)；

156)若||H(k)||₂＞ξ，ξ为常数，令k＝k+1，重复计算均方值P'(k)、滤波增益矩阵H(k)、信号估计值和最小均方误差阵P(k)；

若||H(k)||₂＜ξ，输出

157)选择中最大K个系数的位置作为支撑集Γ，其中K为稀疏信号中非零元素的个数，令：

2.根据权利要求1所述的一种基于卡尔曼滤波的复杂环境下信号重构方法，其特征在于，所述的对信号进行快速压缩包括以下步骤：

21)定义稀疏信号x的维数N，压缩过以后的维数为M，计算压缩比ρ，其计算公式如下：

22)设定类LDPC矩阵的行权重r和列权重l，且

23)对类LDPC矩阵根据测量矩阵的密度或随机产生测量矩阵Φ_M×N中非零元素的位置向量Υ，令Φ_M×N中非零元素值交替定义为1和-1；

24)进行压缩测量，获得测量值y，其计算公式如下：

y＝Φx。

3.根据权利要求1所述的一种基于卡尔曼滤波的复杂环境下信号重构方法，其特征在于，所述的建立信号的先验模型包括以下步骤：

31)定义信号x中的元素取x_i≠0时，使用高斯分布表示其概率分布，概率分布为

信号x中的元素取x_i＝0时，使用Dirac分布来近似概率分布，概率分布为δ(x)；

32)建立信号的先验模型f(x)，

其中：q:＝P_r{x_i≠0},i∈[1,2,...N]。

4.根据权利要求1所述的一种基于卡尔曼滤波的复杂环境下信号重构方法，其特征在于，所述的在二分图上进行置信传播计算包括以下步骤：

41)设校验矩阵H的Tanner图G＝{(V,E)}，V为节点的集合，包含变量节点集合V_b和校验节点V_c，

V_b＝(b₁,b₂,……,b_n)，b_n为变量节点，其与校验矩阵H的各列相对应；V_c＝(c₁,c₂,......,c_m)，c_m表示校验节点，其与校验矩阵H的各行相对应；

E为不同类节点之间相连的边的集合，

42)建立测量矩阵Φ对应的二分图，二分图中每一条边连接变量节点X和测量值对应的校验节点Y，并且每一条边对应测量矩阵中的一个非零元素Φ_ij；

43)进行迭代置信传播计算；

设从变量节点到校验节点的消息编码为信号分量的后验概率概率密度，用μ_i→j表示，从校验节点到变量节点的消息编码为测量分量的概率密度，用μ_j→i表示；

从变量节点到校验节点的消息编码计算公式如下：

从校验节点到变量节点的消息编码计算公式如下：

44)对从变量节点到校验节点的消息编码计算公式和从校验节点到变量节点的消息编码计算公式进行迭代计算，直到消息值不再发生变化，迭代结束，输出信号值x的后验概率密度，如下所示：

5.根据权利要求1所述的一种基于卡尔曼滤波的复杂环境下信号重构方法，其特征在于：还包括对信号估计值进行迭代更新，迭代执行在二分图上进行置信传播计算；采用近似MMSE估计得到信号估计的初始值使用卡尔曼滤波得到信号估计值三个步骤，直到重构精度满足误差要求或达到事先设定的最大迭代次数，输出更新后的信号估计值

6.根据权利要求2所述的一种基于卡尔曼滤波的复杂环境下信号重构方法，其特征在于，所述的进行压缩测量，获得测量值y的计算公式如下：

y＝Φx+n。