CN110596645B - 二维免求逆稀疏贝叶斯学习快速稀疏重构方法 - Google Patents

Abstract

本发明属于信号处理领域，具体涉及一种二维免求逆稀疏贝叶斯学习快速稀疏重构方法，包括以下步骤：S1：对二维稀疏重构问题进行稀疏表示建模；S2：对向量化稀疏信号x以及向量化噪声n进行统计建模；S3：求解向量化稀疏信号x、向量化方差倒数γ与噪声方差倒数α的后验概率；S4：更新辅助变量的矩阵形式Z。与IFSBL方法相比，本发明方法直接对二维信号进行处理，避免了二维信号向量化而产生大矩阵的问题，运算效率明显提升，且显著降低了对计算内存的需求；另一方面，本发明是在统计信号处理框架下实现稀疏重构，与非统计稀疏重构方法相比，具有更易获取全局最优解，对噪声鲁棒性更强，以及算法性能对参数初始化依赖程度不高等优势，工程实用性强。

Description

二维免求逆稀疏贝叶斯学习快速稀疏重构方法

技术领域

本发明属于信号处理领域，具体涉及一种二维免求逆稀疏贝叶斯学习(2Dimentional Inverse-Free Sparse Bayesian Learning，2D-IFSBL)快速稀疏重构方法。

背景技术

稀疏重构是压缩感知技术的核心，可从不完整观测数据中准确重构出稀疏信号。经过不断发展，稀疏重构技术已广泛应用于医学图像处理、计算机视觉、雷达成像等领域。经典稀疏重构算法包括l1正则化、基追踪(Basis Pursuit，BP)、正交匹配追踪(OrthogonalMatching Pursuit，OMP)、稀疏贝叶斯学习(Sparse Bayesian Learning，SBL)等。其中，SBL方法是在统计学理论框架下解决稀疏重构问题，对观测信号、稀疏信号与噪声信号的统计先验建模与后验求解，与l₁正则化、BP、OMP等方法相比，具有更易获取全局最优解、对噪声鲁棒性更强与参数自学习等优势。但SBL方法求解过程中涉及复杂的矩阵求逆操作，计算复杂度高。文献(Duan H,Yang L,Fang J,and Li H.Fast Inverse-Free Sparse BayesianLearning via Relaxed Evidence Lower Bound Maximization[J].IEEE SignalProcessing Letters,2017,24(6):774-778)提出了一种对传统SBL方法的改进方法，这种方法避免了SBL中的矩阵求逆过程，有效提升了运算效率，因而称为免求逆稀疏贝叶斯学习(Inverse-Free Sparse Bayesian Learning，IFSBL)方法。

在实际应用中，还存在二维信号稀疏重构问题，如在雷达成像中，需要对雷达回波进行二维傅里叶变换，而从不完整的雷达回波中重构雷达图像就是二维稀疏重构问题。在处理这类问题时，传统处理流程是首先将二维信号按列堆叠进行向量化，再利用传统稀疏重构方法进行处理，但此时引入了尺寸较大的矩阵，运算效率低，且对内存要求高，l₁正则化、BP、OMP以及SBL等传统方法均不满足工程实际需求。尽管IFSBL方法相对SBL方法提升了运算效率，但在处理二维信号时同样需要将二维信号向量化，无法避免大矩阵运算，运算效率仍然无法满足工程实际要求。

发明内容

本发明要解决的技术问题是如何从不完整观测数据中快速高效地重构二维稀疏信号，且降低对运算内存的要求，以提升工程实用性。

本发明的思路是针对传统稀疏重构方法处理二维信号运算效率低的问题，提出了一种二维免求逆稀疏贝叶斯学习方法。该方法首先对二维稀疏重构问题进行稀疏表示建模，并分别对待重构稀疏信号以及噪声进行统计先验建模；然后基于对待重构稀疏信号的先验以及观测信号的似然函数，通过贝叶斯公式计算稀疏信号的后验概率，从而实现稀疏信号重构。计算过程中直接对二维信号进行处理，而无需对其进行向量化，避免了大矩阵处理，从而有效提升二维稀疏重构的运算效率，并且降低内存需求。

本发明解决其技术问题所采取的技术方案是：一种二维免求逆稀疏贝叶斯学习快速稀疏重构方法，包括以下步骤：

S1对二维稀疏重构问题进行稀疏表示建模：

从不完整观测数据中重构二维稀疏信号，首先需要建立两者的关系表达式，如下式所示：

Y＝AXB^T+N (1)

其中

分别表示观测信号、左字典矩阵、稀疏信号、右字典矩阵以及噪声，

表示实数集，P、Q、M、N分别表示矩阵的维数。在二维稀疏重构问题中，观测信号维度小于稀疏信号，即P＜M且Q＜N。由于观测信号与稀疏信号均为二维信号，无法直接采用传统稀疏重构方法从观测信号Y中重构稀疏信号X，首先需要将两者进行向量化：

y＝Φx+n (2)

其中y＝vec(Y)、x＝vec(X)、n＝vec(N)分别表示观测信号、稀疏信号与噪声的向量化，vec(·)表示按列堆叠对矩阵进行向量化；

表示矩阵向量化对应的字典矩阵：

其中

表示矩阵的克罗内克积。字典矩阵Φ的维度较大，导致采用传统稀疏重构方法求解式(2)所示问题运算效率低，且对内存要求高，如当观测信号Y与稀疏信号X的尺寸分别为64×64与128×128时，字典矩阵Φ的尺寸为8192×8192，直接对其进行处理运算效率无法满足工程需求。

S2对向量化稀疏信号x以及向量化噪声n进行统计建模：

统计建模过程中，采用两层结构的高斯分布对向量化稀疏信号x进行统计建模。在第一层中，假设向量化稀疏信号x的各元素服从均值为零且相互独立的高斯分布：

其中γ表示向量化的方差倒数，x_n、γ_n分别表示向量化稀疏信号x和向量化方差倒数γ的第n个元素。

在第二层中，假设向量化方差倒数γ服从伽马分布：

其中a、b分别表示伽马分布的形状参数和尺寸参数。

同样采用两层结构的高斯分布对向量化噪声n进行统计建模。在第一层中，假设向量化噪声n服从零均值高斯分布，则向量化观测信号y的似然函数同样服从高斯分布：

其中

表示高斯分布，I表示单位矩阵，α表示噪声方差的倒数。

在第二层中，假设噪声方差的倒数α服从伽马分布：

其中c、d分别表示伽马分布的形状参数和尺度参数。

S3求解向量化稀疏信号x、向量化方差倒数γ与噪声方差倒数α的后验概率：

SBL方法实现稀疏重构的过程本质上是基于观测信号的似然函数与稀疏信号的先验模型，通过贝叶斯公式推导稀疏信号的后验概率的过程，而传统SBL方法求解向量化稀疏信号x后验概率时涉及复杂的矩阵求逆过程，运算效率低。免求逆稀疏贝叶斯学习方法(IFSBL)(Duan H,Yang L,Fang J,and Li H.Fast Inverse-Free Sparse BayesianLearning via Relaxed Evidence Lower Bound Maximization[J].IEEE SignalProcessing Letters,2017,24(6):774-778)有效避免了矩阵求逆，从而提升了运算效率。具体而言，IFSBL方法基于式(5)所示向量化观测信号y的似然函数，以及式(3)所示向量化稀疏信号x的先验概率，所得向量化稀疏信号x的后验概率为：

其中期望μ与协方差矩阵Σ如下式所示：

其中

为辅助变量，<·>表示求期望算子，T为常量，其取值略大于利普希茨常数：T＝2λ_max(Φ^TΦ)+10^-10，其中λ_max(Φ^TΦ)表示Φ^TΦ的最大特征值。diag(·)表示利用向量生成对角矩阵。由式(9)可知，协方差矩阵Σ为对角矩阵，因此对其求逆可以通过对各元素分别取倒获得，运算效率明显提升。然而式(8)中涉及大矩阵Φ的计算，运算效率低，且对内存要求高。针对此问题，本发明对式(8)所示期望求解过程进行修正，使之可以直接处理二维数据，以避免矩阵向量化导致的大矩阵计算。具体而言，将式(9)代入式(8)可得：

其中

表示矩阵各元素分别相除，1_MN×1表示尺寸为MN×1，且各元素全为1的向量。进一步将

代入式(10)，并利用克罗内克积的相关性质，可得：

上式可以重新变换为矩阵形式：

其中

分别为期望μ、辅助变量

以及向量化方差倒数γ的矩阵形式，1_M×N表示1_MN×1的矩阵形式。如式(12)所示，U的计算过程只涉及到最大尺寸为M×N的矩阵运算，因此运算效率显著高于式(8)，且降低了对内存的要求。

进一步计算向量化方差倒数γ的矩阵形式Γ的后验概率，通过贝叶斯公式可得，其同样服从伽马分布：

其中形状参数

和尺度参数

如下式所示：

其中X_m，n表示稀疏信号X的第m行、n列元素，

的期望

的计算表达式为：

其中U_m，n、Γ_m，n分别表示U和Γ的第m行、n列元素。

最后计算噪声方差倒数α的后验概率，由贝叶斯公式可知，其同样服从伽马分布：

其中，α后验概率的形状参数

与尺度参数

分别为：

其中sum(·)、⊙分别表示矩阵内各元素相加以及两个矩阵每个元素分别相乘。

式(12)、(18)中涉及的期望<α>与<Γ>可分别由式(16)与(13)所示后验概率获得，如下式所示：

<Γ_m，n>表示向量化方差倒数γ期望的矩阵形式<Γ>中各元素的期望；

S4更新辅助变量的矩阵形式Z：

计算出向量化稀疏信号x、向量化方差倒数γ与噪声方差倒数α的后验概率后，进一步通过下式更新辅助变量的矩阵形式Z：

Z＝U (21)

综上所述，本发明的处理流程可概括为：循环迭代式(12)、式(19)、式(20)以及式(21)，直至相邻两次迭代所得U之间的相对误差达到额定门限(如：||U^(k+1)-U^(k)||₂/||U^(k)||₂＜10^-3，其中U^(k)、U^(k+1)分别表示第k次和第k+1次迭代所得U，||·||₂表示矩阵的2范数)，最后一步所得U即为重构的稀疏信号。此外，在迭代处理之前需要对参数进行初始化，具体而言：U初始化为

其中

表示矩阵的伪逆；噪声方差的倒数α初始化为α＝0.1var(Y)^-1，其中var(·)表示求方差算子；式(4)和式(6)中所示伽马分布中的参数a、b、c、d初始化为a＝b＝c＝d＝10^-6；向量化方差倒数矩阵形式Γ初始化为Γ＝1_M×N。

本发明取得的有益效果为：通过本发明可快速实现二维稀疏信号重构，与IFSBL方法相比，本发明方法直接对二维信号进行处理，避免了二维信号向量化而产生大矩阵的问题，运算效率明显提升，且显著降低了对计算内存的需求；另一方面，本发明所提2D-IFSBL方法是在统计信号处理框架下实现稀疏重构，与非统计稀疏重构方法相比，具有更易获取全局最优解，对噪声鲁棒性更强，以及算法性能对参数初始化依赖程度不高等优势，工程实用性强。

附图说明

图1本发明的实施流程图；

图2不同观测信号尺寸条件下的稀疏重构结果；

图3不同信噪比(SNR)条件下的稀疏重构结果；

图4不同观测信号尺寸下不同方法所重构稀疏信号的相对估计误差比较；

图5不同SNR条件下不同方法所重构稀疏信号的相对估计误差比较；

图6不同观测信号尺寸下不同方法计算时间比较。

具体实施方式

下面结合附图对本发明进行进一步说明：

图1为本发明总处理流程。

本发明所述一种二维免求逆稀疏贝叶斯学习快速稀疏重构方法，包括以下四个步骤：

S1：对二维稀疏重构问题进行稀疏表示建模；

S2：对向量化稀疏信号x以及向量化噪声n进行统计建模；

S3：求解向量化稀疏信号x、向量化方差倒数γ与噪声方差倒数α的后验概率；

S4：更新辅助变量的矩阵形式Z。

采用仿真数据进行实验，比较本发明方法与IFSBL方法(Duan H,Yang L,Fang J,and Li H.Fast Inverse-Free Sparse Bayesian Learning via Relaxed EvidenceLower Bound Maximization[J].IEEE Signal Processing Letters,2017,24(6):774-778)以及二维平滑0范数方法(2D-SL0)(Aboozar G,Massoud B Z,Christian J.Sparsedecomposition of two dimensional signals[C].2009IEEE International Conferenceon Acoustics,Speech and Signal Processing,Taipei,2009:3157-3160)的性能。实验过程中，首先模拟产生一个尺寸为M×N的二维稀疏信号X，该信号包括K个非零点；左字典矩阵A和右字典矩阵B的尺寸分别为P×M、Q×N，两者的各元素均服从标准正态分布；噪声N为高斯白噪声；最后通过式(1)产生观测信号Y。

在实验一中，假设二维稀疏信号X包含K＝9个非零点，各非零点取值为1，且均匀分布于图像中央，以方便对比。二维稀疏信号X的尺寸设为80×80，即M＝N＝80。观测信号Y的信噪比(SNR)设为20dB。改变观测信号Y的尺寸：设P＝Q＝25、40、55，以比较2D-SL0、IFSBL以及本发明所提2D-IFSBL三种方法在不同观测信号尺寸条件下的性能。三种方法所重构的稀疏信号如图2所示，其中第一、二、三行分别对应P＝25、40、55条件下的稀疏重构结果。由图可知，本发明所提2D-IFSBL与IFSBL方法所得稀疏信号完全一致，当P＝25时，2D-SL0重构结果较优，而当P＝55时，本发明所提2D-IFSBL与IFSBL方法所得结果受噪声影响较小，重构的稀疏信号优于2D-SL0方法。

第二个实验比较三种方法在不同SNR条件下的性能。实验参数设置如下：P＝Q＝55、M＝N＝80、K＝9、SNR分别设为15dB、7.5dB以及0dB，三种方法在不同SNR条件下所得稀疏信号如图3所示。由图可知，本发明所提2D-IFSBL与IFSBL方法所重构的稀疏信号受噪声影响程度明显低于2D-SL0方法，尤其是在SNR为0dB的条件下，2D-SL0方法基本失效，而本发明所提2D-IFSBL与IFSBL方法仍然可以准确重构稀疏信号，表明两者对噪声的鲁棒性优于2D-SL0方法。

进一步通过蒙特卡洛实验比较算法性能。首先比较不同观测信号尺寸条件下三种方法所重构稀疏信号的相对误差，实验参数设置如下：M＝N＝80、K＝10、SNR＝20dB，改变P、Q的值，变化范围为10到60，变化步长为5。在每个不同观测信号尺寸条件下，重复进行100次蒙特卡洛实验，每次实验随机对A、B、X以及N进行取值，其中，稀疏信号X中的K＝10个非零点取值服从标准正态分布，且在X中的位置完全随机。100次蒙特卡洛实验中，三种方法所得平均稀疏重构相对误差如图4所示。由图可知，本发明所提2D-IFSBL方法与IFSBL方法性能完全一致，当P＜40时，2D-SL0方法较优，而当P＞40时，2D-IFSBL与IFSBL方法所得相对误差低于2D-SL0方法。

进一步通过蒙特卡洛实验比较不同SNR条件下三种方法性能。实验参数设置如下：M＝N＝80、P＝Q＝55、K＝10，SNR变化范围为-12dB到16dB，变化步长为2dB。在每个SNR条件下，分别进行100次蒙特卡洛实验，并记录三种方法的平均稀疏重构相对误差，结果如图5所示。由图可知，本发明所提2D-IFSBL方法与IFSBL方法全面优于2D-SL0方法，在任何SNR条件下均获得了较低的稀疏重构相对误差，验证了其较强的鲁棒性。

最后比较三种方法的运算效率。同样采用蒙特卡洛实验，实验参数设置如下：

SNR＝20dB、M＝N。改变N的值，取值范围为50到250，步长为10，每个N条件分别进行100次蒙特卡洛实验，记录三种算法的平均计算时间(计算平台：Intel(R)Corei7-8550U@1.8GHz)，结果如图6所示。由图可知，2D-SL0方法运算效率最高，本发明所提2D-IFSBL方法其次，IFSBL方法最低。IFSBL方法由于涉及大矩阵计算，对计算内存要求较高，当N＞150时，计算内存超限，算法失效。相比而言，本发明所提2D-IFSBL方法在保持与IFSBL性能一致的前提下，运算效率最高提升100余倍，且降低了对计算内存的要求，当N＝250时，仍然只需约5秒的时间重构稀疏信号，满足工程实际需求。

以上实验结果表明，本发明对传统IFSBL进行修正，使之可以直接处理二维数据。在二维稀疏重构问题中，本发明方法在保持传统IFSBL方法性能的前提下，运算效率提升100余倍，并且显著降低了对计算内存的要求。与2D-SL0方法相比，本发明方法属于统计稀疏重构方法，具有鲁棒性较强的优势。尽管在观测数据尺寸较小的条件下，算法性能低于2D-SL0方法，在低SNR条件下算法性能明显优于2D-SL0方法，适用于SNR较低的应用场景。

Claims (4)

1.一种二维免求逆稀疏贝叶斯学习快速稀疏重构方法，其特征在于，该方法包括以下步骤：

S1对二维稀疏重构问题进行稀疏表示建模：

从不完整观测数据中重构二维稀疏信号，首先需要建立两者的关系表达式，如下式所示：

Y＝AXB^T+N (1)

其中

分别表示观测信号、左字典矩阵、稀疏信号、右字典矩阵以及噪声，

表示实数集，P、Q、M、N分别表示矩阵的维数，在二维稀疏重构问题中，观测信号维度小于稀疏信号，即P＜M且Q＜N；由于观测信号与稀疏信号均为二维信号，无法直接采用传统稀疏重构方法从观测信号Y中重构稀疏信号X，首先需要将两者进行向量化：

y＝Φx+n (2)

其中y＝vec(Y)、x＝vec(X)、n＝vec(N)分别表示观测信号、稀疏信号与噪声的向量化，vec(·)表示按列堆叠对矩阵进行向量化；

表示矩阵向量化对应的字典矩阵：

其中

表示矩阵的克罗内克积；

S2对向量化稀疏信号x以及向量化噪声n进行统计建模：

统计建模过程中，采用两层结构的高斯分布对向量化稀疏信号x进行统计建模：在第一层中，假设向量化稀疏信号x的各元素服从均值为零且相互独立的高斯分布：

其中γ表示向量化的方差倒数，x_n、γ_n分别表示向量化稀疏信号x和向量化方差倒数γ的第n个元素；

在第二层中，假设向量化方差倒数γ服从伽马分布：

其中a、b分别表示伽马分布的形状参数和尺寸参数；

同样采用两层结构的高斯分布对向量化噪声n进行统计建模：在第一层中，假设向量化噪声n服从零均值高斯分布，则向量化观测信号y的似然函数同样服从高斯分布：

其中

表示高斯分布，I表示单位矩阵，α表示噪声方差的倒数；

在第二层中，假设噪声方差的倒数α服从伽马分布：

其中c、d分别表示伽马分布的形状参数和尺度参数；

S3求解向量化稀疏信号x、向量化方差倒数γ与噪声方差倒数α的后验概率：

IFSBL方法基于式(5)所示向量化观测信号y的似然函数，以及式(3)所示向量化稀疏信号x的先验概率，所得向量化稀疏信号x的后验概率为：

其中期望μ与协方差矩阵Σ如下式所示：

其中

为辅助变量，<·>表示求期望算子，T为常量，其取值略大于利普希茨常数；diag(·)表示利用向量生成对角矩阵；由式(9)可知，协方差矩阵Σ为对角矩阵，因此对其求逆可以通过对各元素分别取倒获得，运算效率明显提升；然而式(8)中涉及大矩阵Φ的计算，运算效率低，且对内存要求高；针对此问题，本发明对式(8)所示期望求解过程进行修正，使之可以直接处理二维数据，以避免矩阵向量化导致的大矩阵计算；具体而言，将式(9)代入式(8)可得：

其中

表示矩阵各元素分别相除，1_MN×1表示尺寸为MN×1，且各元素全为1的向量；进一步将

代入式(10)，并利用克罗内克积的相关性质，可得：

上式可以重新变换为矩阵形式：

其中

分别为期望μ、辅助变量z以及向量化方差倒数γ的矩阵形式，1_M×N表示1_MN×1的矩阵形式；

进一步计算向量化方差倒数γ的矩阵形式Γ的后验概率，通过贝叶斯公式可得，其同样服从伽马分布：

其中形状参数

和尺度参数

如下式所示：

其中X_m,n表示稀疏信号X的第m行、n列元素，

的期望

的计算表达式为：

其中U_m,n、Γ_m,n分别表示U和Γ的第m行、n列元素；

最后计算噪声方差倒数α的后验概率，由贝叶斯公式可知，其同样服从伽马分布：

其中，α后验概率的形状参数

与尺度参数

分别为：

其中sum(·)、⊙分别表示矩阵内各元素相加以及两个矩阵每个元素分别相乘；

式(12)、(18)中涉及的期望<α>与<Γ>可分别由式(16)与(13)所示后验概率获得，如下式所示：

<Γ_m,n>表示向量化方差倒数γ期望的矩阵形式<Γ>中各元素的期望；

S4更新辅助变量的矩阵形式Z：

计算出向量化稀疏信号x、向量化方差倒数γ与噪声方差倒数α的后验概率后，进一步通过下式更新辅助变量的矩阵形式Z：

Z＝U (21)

综上所述，本发明的处理流程可概括为：循环迭代式(12)、式(19)、式(20)以及式(21)，直至相邻两次迭代所得U之间的相对误差达到额定门限，最后一步所得U即为重构的稀疏信号。

2.一种根据权利要求1所述二维免求逆稀疏贝叶斯学习快速稀疏重构方法，其特征在于：在迭代处理之前对参数进行初始化，具体而言：U初始化为

其中

表示矩阵的伪逆；噪声方差的倒数α初始化为α＝0.1var(Y)^-1，其中var(·)表示求方差算子；式(4)和式(6)中所示伽马分布中的参数a、b、c、d初始化为a＝b＝c＝d＝10^-6；向量化方差倒数矩阵形式Γ初始化为Γ＝1_M×N。

3.一种根据权利要求1或2所述二维免求逆稀疏贝叶斯学习快速稀疏重构方法，其特征在于：S3中常数T的具体取值为：T＝2λ_max(Φ^TΦ)+10^-10，其中λ_max(Φ^TΦ)表示Φ^TΦ的最大特征值。

4.一种根据权利要求1所述二维免求逆稀疏贝叶斯学习快速稀疏重构方法，其特征在于：S4中的额定门限为||U^(k+1)-U^(k)||₂/||U^(k)||₂＜10^-3，其中U^(k)、U^(k+1)分别表示第k次和第k+1次迭代所得U，||·||₂表示矩阵的2范数。

Images