本发明属于信号处理领域,可以应用于图像信息的稀疏重构,特别涉及基于多观测向量模型的低复杂度结构化重构方法。
背景技术:
根据压缩感知(compressivesensing,cs)理论,若信号在某一变换域满足一定的稀疏性,则能够利用远低于奈奎斯特采样率的速率对信号进行采样,并实现对信号的精确重构。目前,压缩感知技术在许多工程领域都有应用,例如欠奈奎斯特采样系统、压缩成像系统、压缩传感网络等。特别地,压缩感知技术还能够被应用于稀疏图像重构问题上,并已经在医学成像等重要领域得到相当广泛的应用。
然而,现有的压缩感知算法在解决图像重构问题时大多基于单一观测向量(smv)模型,即将整个稀疏图像视为一维向量进行处理;该处理方式所存在的问题是:当我们需要对高分辨率图像进行重构时,由大规模图像拉伸成的一维向量维度相当高,因此算法中涉及到的观测矩阵和协方差矩阵规模庞大,从而会对现有的计算平台带来相当大的存储压力和计算压力。比如,如果采用基于smv的压缩感知算法处理一张1024×1024维的图像,那么算法中涉及到的该图像矩阵展成一个列向量后维度将超过1.04×106维,存储这样规模的协方差矩阵需要超过64gb的计算平台内存;即使计算平台内存足够大,能够满足矩阵存储需求并使这些方法可以成功运行,但是由此耗费的计算时间和占用的计算资源也是难以承受的。因此,在高分辨率图像重构问题中,仅仅利用基于smv模型的压缩感知算法不具有现实意义。
为此,本发明充分考虑图像稀疏重构的特点,利用mmv模型对问题进行重新表征。mmv模型将原始信号的每一个列向量作为独立的子信号处理,并已经在阵列信号处理,认知无线电等领域已经取得了相应的研究成果,但并非被应用于图像的稀疏重构问题。因此,如果我们将多观测向量模型引入大规模图像重构问题,并且将图像矩阵的每一个列向量看作是mmv模型中的一个观测向量,则这样的模型表征方式不再涉及将图像转化为高维度单一向量的操作,相应的问题规模将会得到明显的降低,从而使得利用压缩感知的原理解决大规模图像的重构成为可以在大多数计算平台上实现的工程问题。但是,如果将现有的解决其他领域问题的mmv重构算法直接应用于解决图像重构问题时,则忽略了图像信息所具有的结构特性,从而造成信息的严重丢失,以及重构精度的恶化。
因此,为进一步克服直接利用mmv模型时可能出现的重构精度恶化问题,本发明所述重构方法在利用mmv模型对问题进行重新描述的基础之上,引入了图像信息中潜在的结构化先验约束,具体体现为两方面的结构化先验信息:一方面,图像的在小波变换域上的稀疏系数具有四叉树状结构的概率依赖关系,另一方面,图像矩阵的列向量之间并非独立,而是具有一定的空间相关性。
综上所述,本发明通过将mmv模型引入高分辨率图像重构问题,能够实现使其存储复杂度和计算复杂度大幅下降;进一步地,通过引入图像信息的结构化先验约束,保证了mmv模型下图像重构的精度。
技术实现要素:
(一)要解决的技术问题
本发明的目的是解决高分辨率图像的稀疏重构所面临的存储复杂度与时间复杂度较高,以及重构精度受限的问题。
(二)技术方案
为了解决上述技术问题,本发明提出一种低复杂度的基于多观测向量表征的结构化重构方法。本发明的目的是通过下述技术方案实现的:输入观测矩阵,多观测向量组成的矩阵,小波变换层数,每层小波系数数目;初始化系统模型参数,图像小波系数,最大迭代次数,精度性能指标;用基于马尔可夫链蒙特卡洛方法的吉布斯采样对图像小波系数进行重构,再对重构后小波系数矩阵进行小波逆变换,从而得到重构图像,
1、用结构化先验约束提高稀疏图像重构精度降低复杂度方法,其特征在于,是一种对分辨率为n×n=1024×1024的遥感图像x0利用多观测向量模型进行重新描述后,再利用x0在小波变换域上的稀疏系数具有四叉树状结构的概率依赖关系以及多观测向量矩阵在列向量之间的空间相关性这两方面的结构化先验约束来重构出分辨率为n×n=1024×1024的遥感图像x*的方法,是在计算机中依次按以下步骤实现的:
步骤(1),向计算机输入以下参数:
n×n=1024×1024的遥感图像矩阵x0,简称矩阵x0,
用matlab工具箱生成基于所述矩阵x0的四叉树状结构的二维离散小波系数变换矩阵θ,其分辨率为n×n=1024×1024,二维离散小波变换的层数s,基于所述矩阵x0的分辨率为n×n=1024×1024,所述层数的最大值smax=8,层号s=1,2,...,8,四叉树状结构的小波系数中每一层小波系数的个数为(nc)s,根结点在第一层s1,任取s=4,
用matlab工具箱生成的m×n的随机观测矩阵,m≤n,观测矩阵简称矩阵φ,下同,
基于所述矩阵θ生成的m×n的多观测向量组成的矩阵y,这是在以所述矩阵φ为观测矩阵的条件下对于所述矩阵x0的观测值,即y=φθ,
输入重构后得到的遥感图像x*的两个图像性能指标:峰值信噪比psnr的设定值,均方误差mse的设定值,
步骤(2),初始化系数模型参数,包括:
①不同层号的小波系数的正则化协方差矩阵c1,c2,...,c4,初始化时全为单位矩阵,简称cs={c1,c2,...,c4},
②控制各cs值的参数λ,令λ初始化为2,
③观测噪声所服从的高斯分布的精度参数,为观测噪声的方差的倒数αn,令其初始化为αn=10,
④不同层号的小波系数所服从的高斯分布的精度参数α1,α2,...,α4,为各层小波系数列向量之间方差的倒数,初始化时都等于1,
⑤设定:控制各层小波系数取零的概率:
对于第一层s1:控制小波系数矩阵θ(1)中的元素
对于第二层至第四层s2,s3,s4:
若小波系数矩阵θ(s),2≤s≤4的元素
若小波系数矩阵θ(s),2≤s≤4的元素
⑥系统超参数a0,b0,c0,d0,e0,f0,所述超参数是一种对步骤(2)中的参数所服从的先验概率分布参数化后得到的参数,其中:
a0,b0,c0,d0在初始化时全部等于10-6,用于控制概率分布函数的参数,为常数,
符号e0包括:(ec)1,以及(ecz)s,2≤s≤4:
(ec)1用于控制w1所服从的概率密度分布函数,取值如下:
(ec)1=0.9×(nc)1,
(ecz)2,(ecz)3,(ecz)4分别用于控制2≤s≤4时w02,w03,w04所服从的概率密度分布函数,取值如下:
(ecz)2=0.1×(nc)2,(ecz)3=0.1×(nc)3,(ecz)4=0.1×(nc)4,
(eco)2,(eco)3,(eco)4分别用于控制2≤s≤4时w12,w13,w14所服从的概率密度分布函数,取值如下:
(eco)2=0.5×(nc)2,(eco)3=0.5×(nc)3,(eco)4=0.5×(nc)4,
符号f0包括:(fc)1,以及(fcz)s,2≤s≤4:
(fc)1用于控制w1所服从的概率密度分布函数,取值如下:
(fc)1=0.1×(nc)1,
(fcz)2,(fcz)3,(fcz)4分别用于控制2≤s≤4时w02,w03,w04所服从的概率密度分布函数,取值如下:
(fcz)2=0.9×(nc)2,(fcz)3=0.9×(nc)3,(fcz)4=0.9×(nc)4,
(fco)2,(fco)3,(fco)4分别用于控制2≤s≤4时w12,w13,w14所服从的概率密度分布函数,取值如下:
(fco)2=0.5×(nc)2,(fco)3=0.5×(nc)3,(fco)4=0.5×(nc)4,
⑦初始化控制迭代过程的参数t,tmax,其中:
t为迭代次数的序号,tmax为最大迭代次数,故迭代次数的序号的设定值为t={1,2,...,tmax};
步骤(3),依次按以下过程进行计算小波系数矩阵θ的迭代过程:
步骤(3.1),依次按以下步骤进行第一次迭代,令t=1,
步骤(3.1.1),设定约束条件及步骤(2)中各种系统模型参数:
约束条件为:n×n=1024×1024的遥感图像矩阵x0,n×n=1024×1024的小波系数矩阵θ为四叉树结构,smax=8,任取s=4,
初始时,各层的所述小波系数矩阵θ为全零矩阵,
步骤(3.1.2),依次按下式计算根据系统模型参数和小波系数的初始值,计算中间变量
其中:
||·||2表示取向量或矩阵的2-范数,
αs和αn表示是初始时的系统模型参数值,
φ-i表示在φ中除去其第i列之外的其它所有列组成的矩阵,
φ·i表示φ矩阵的第i列的所有元素组成的向量,
θ(-i)j表示在向量θ·j中除去其第i个元素之外的其它所有元素组成的向量,
步骤(3.1.3),根据步骤(3.1.2)中得到的中间变量
其中:δ0代表原点处的单位冲击函数,
根据小波系数θ所服从的后验概率分布p(θij),采样得到矩阵θ的第i行、第j列元素θij的值,即得到矩阵θ的采样值,
步骤(3.1.4),根据系统模型参数的初始值和步骤(3.1.3)中得到的矩阵θ的采样值,计算第一次迭代中系统模型参数所服从的后验概率密度分布,以及矩阵cs的第一次迭代结果:
p(w1)=beta((ec)1+||θ(1)||0,(fc)1+n1-||θ(1)||0),
其中:
角标s的取值集合均为{1,...,4},注意,仅在输入的s≥2时,且当前考虑的层数s大于1时,才存在w0s和w1s,
||·||0表示矩阵(或向量)的0-范数,
θ(s)是第s层的小波变换系数按照自然顺序组合形成的矩阵,
rs是θs的行数,
θs,k表示第s层的第k个小波系数,θπ(s,k)是θs,k的父亲节点的系数值,
p(x)=gamma(a,b)表示随机变量x服从参数为(a,b)的伽马分布,其展开形式为
p(x)=beta(a,b)表示随机变量x服从参数为(a,b)的贝塔分布,其展开形式为
根据上述系统模型参数αn、αs、w1、w0s和w1s所服从的后验概率密度分布,对系统模型参数进行采样,并更新相应的系统模型参数取值,用于下一次的迭代过程,
步骤(3.1.5),存储步骤(3.1.3)中得到的矩阵θ的采样值,第一次迭代过程正式结束,
步骤(3.2),按照步骤(3.1.1)~步骤(3.1.5)的描述,进行第二次至第tmax次迭代的计算,
步骤(3.3),对所有迭代次数的序号,t={1,...,tmax},计算每次迭代中存储得到的所有矩阵θ的算术平均值,得到的平均值记作θ*,对平均值θ*进行s层的二维小波逆变换,得到x*,
步骤(3.4),统计出重构图像x*相对于原始图像x0的误差指数,psnr和mse,如果这两个指数满足输入的性能指标要求,则算法停止,否则调整步骤(1)中输入的小波系数层数s,且令s:=s+1,重新跳转至步骤(2),再次运行初始化和迭代计算过程,直到重构图像相对于原始图像的误差指数满足性能指标要求。
有益效果
本发明给出了一种低复杂度的高分辨率图像的稀疏重构方法。该重构方法基于mmv模型进行设计;为进一步提高重构精度,该方法引入了图像的结构化先验约束。该方法具有以下优点:(1)通过引入mmv模型使得问题的存储复杂度由o(n4)降低到o(n2),计算复杂度由o(m2n2)降低到o(mn),从而为高分辨率图像重构问题提供了可行的方案;(2)通过引入结构化先验概率分布约束,提高了mmv模型下图像重构的精确度;(3)利用该方法中推导得到的吉布斯采样方法,能够推断得到所述结构化mmv模型下图像系数与模型参数的理论最优解。
附图说明
图1是本发明实现稀疏图像重构的流程图。
图2刻画了图像多层离散小波变换的四叉树状概率依赖结构,其中s是多层小波变换的层数。
图3展示了本发明建立的基于结构化先验信息的层次化贝叶斯模型。其中a0,b0,c0,d0,e0,f0是刻画图像系数和模型参数先验分布的超参数,θ,y分别是图像的小波系数和基于多观测向量模型获得的稀疏观测矩阵,αs,cs,ωs是刻画θ的概率分布特性的参数,αn是刻画y的概率分布特性的参数。
图4是本发明提出的图像重构方法(cm-ssbl)与已有的基于mmv模型的稀疏重构算法(omp-mmv、m-focuss和msbl)在处理大规模图像稀疏重构问题时的效果对比。横轴是归一化采样率m/n,纵轴是相对重构误差;
图5是本发明提出的图像重构方法(cm-ssbl)与已有的基于mmv模型的稀疏重构算法(omp-mmv、m-focuss和msbl)在处理大规模图像稀疏重构问题时的效果对比。横轴是归一化采样率m/n,纵轴是峰值信噪比(psnr);
图6是采用分辨率为1024×1024的标准测试图像“cameraman”,在归一化采样率为0.4的情况下,对omp-mmv、m-focuss、msbl,以及本发明所提出的cm-ssbl四种方法进行图像重构实验的直观效果展示图。图6(a)是omp-mmv算法的重构结果及对应的psnr,图6(b)是m-focuss算法的重构结果及对应的psnr,图6(c)是msbl算法的重构结果及对应的psnr,图6(d)是本发明提出的cm-ssbl算法的重构结果及对应的psnr。
图7是优选与迭代控制相关的参数的过程中的参考曲线,横轴是峰值信噪比,纵轴是迭代次数,为了表现迭代次数的变化对重构误差性能的影响,图7截取了优选过程进行至最终,迭代次数相对较少时的峰值信噪比变化情况。
具体实施方式
下面结合附图和实施例,对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明,但不是限制本发明的范围。
如图1所示,本发明的目的是通过下述技术方案实现的:输入观测矩阵,多观测向量组成的矩阵,小波变换层数,每层小波系数数目;初始化系统模型参数,图像小波系数,最大迭代次数,精度性能指标;用基于马尔可夫链蒙特卡洛方法的吉布斯采样对图像小波系数进行重构,再对重构后小波系数矩阵进行小波逆变换,从而得到重构图像。本发明能对实际重构图像应用所需的精度性能指标,利用多观测向量模型降低复杂度,同时又通过自适应的方法调整迭代次数,使得重构图像的精度满足性能指标要求。
现在以对一幅维度为n×n(n=1024)的标准测试图像x进行稀疏重构为例,
首先确定稀疏观测矩阵,利用matlab随机数生成工具箱生成m×n(m=512)的观测矩阵φ,
根据图像稀疏表征的基本原理,在矩阵φ为表征基底的条件下对于原始遥感图像x0进行观测,得到观测向量组成的m×n的矩阵y,
设定二维离散小波变换的层数s=4,则每层小波系数的个数也随之确定:在n=1024,s=4的情况下,第0层小波系数个数为(nc)0=3×642=12288,第1层小波系数个数为(nc)1=3×1282=49152,第2层小波系数个数为(nc)2=3×2562=196608,第3层小波系数个数为(nc)3=3×5122=786432;
设定重构图像性能指标psnr=36(此处以mse=16.3设定性能指标,效果是等价的);
接下来,根据本发明提出的基于结构化先验信息的稀疏图像重构算法,依次进行如下求解步骤:
初始化n×n的小波系数θ,令θ等于全零矩阵;
初始化系统模型参数,包括:
不同层号的小波系数的正则化协方差矩阵c0,c1,...,cs-1,初始时全为单位矩阵,简称cs,cs={c0,c1,...,cs-1};
控制cs大小的参数λ,令λ=2;
观测噪声所服从的高斯分布的精度参数,为观测噪声的方差的倒数,记为αn,αn=10;
不同层号的小波系数所服从的高斯分布的精度参数α0,α1,...,αs-1,分别是第s层小波系数列向量之间协方差的倒数,初始化时都等于1;
w1控制第一层小波系数取零的概率,对于下标为层序号s,s=2,...,s,w02,...,w0s和w12,...,w1s都是s>1时控制小波系数取零的概率,
w1为第一层s=1时θij|s=1时的
w0s为s>1且θπ(ij)=0时的
w1s为s>1且θπ(ij)=1时的
初始化系统超参数a0,b0,c0,d0,e0,f0,其中:
a0,b0,c0,d0在初始化时全部等于10-6,
符号e0包括:(ec)1,以及(ecz)2,...,(ecz)s和(eco)2,...,(eco)s:
(ec)1用于控制w1所服从的概率密度分布函数,取值如下:
(ec)1=0.9×(nc)1,
(ecz)2,...,(ecz)s用于控制w02,...,w0s所服从的概率密度分布函数,取值如下:
(ecz)2=0.1×(nc)2,...,(ecz)s=0.1×(nc)s,
(eco)2,...,(eco)s用于控制w12,...,w1s所服从的概率密度分布函数,取值如下:
(eco)2=0.5×(nc)2,...,(eco)s=0.5×(nc)s,
符号f0包括:(fc)1,以及(fcz)2,...,(fcz)s和(fco)2,...,(fco)s:
(fc)1用于控制w1所服从的概率密度分布函数,取值如下:
(fc)1=0.1×(nc)1,
(fcz)2,...,(fcz)s用于控制w02,...,w0s-1所服从的概率密度分布函数,取值如下:
(fcz)2=0.9×(nc)2,...,(fcz)s=0.9×(nc)s,
(fco)2,...,(fco)s用于控制w12,...,w1s-1所服从的概率密度分布函数,取值如下:
(fco)2=0.5×(nc)2,...,(fco)s=0.5×(nc)s,
初始化控制迭代过程的参数t,tmax,其中:
t为迭代次数的序号,
设定最大迭代次数tmax=55。
依次按以下过程进行计算小波系数θ的迭代过程,
步骤(1),根据当前小波系数层数s确定系统模型参数的初始值;
步骤(2),tmax为当前设置好的常数,代表最多迭代计算的周期数,
步骤(3)依次按如下步骤实现第一次迭代计算,序号为t=t(1)=1,
步骤(3.1)以各系统模型参数和小波系数的初始值作为第一次迭代的初始值,
步骤(3.2),依次按下式计算根据系统模型参数和小波系数的初始值,计算中间变量
步骤(3.3),计算小波系数θ所服从的后验概率分布p(θij):
步骤(3.4),根据小波系数θ所服从的后验概率分布p(θij),采样得到矩阵θ的第i行、第j列元素θij的值,即得到矩阵θ的采样值,根据系统模型参数的初始值和矩阵θ的采样值,计算第一次迭代中系统模型参数所服从的后验概率密度分布,以及矩阵cs的第一次迭代结果,
步骤(4),按照步骤(3.1)~步骤(3.4)的描述,进行第二次至第tmax次迭代的计算,
步骤(5),对所有的θ(t),
步骤(6),统计出重构图像x*的psnr,此时仍未满足psnr=36的性能指标,则将tmax增大10,此时tmax=65,并将小波系数层数s增加1,重新跳转至步骤(1),再次运行迭代计算过程,此次得到的重构图像的psnr>36,则停止运算,输出重构图像。
在上述算法执行步骤的基础上,这里对本发明所述的引入结构化先验约束的mmv重构算法相较于已有mmv重构算法所能提供的重构性能增益进行说明。我们将重构得到的小波系数的相对重构误差和峰值信噪比作为评价标准,分别考察已有的3种mmv重构算法omp-mmv,m-focuss,msbl,以及本发明提出的cm-ssbl算法对分辨率为1024×1024的标准测试图像“cameraman”进行重构的效果。
图4展示了相对重构误差随归一化采样率m/n的变化情况,从图中可以看出,在归一化采样率从0.25到0.55的变化范围内,cm-ssbl算法与omp-mmv,m-focuss,msbl算法相比在误差性能方面有显著的提升。以归一化采样率为0.4时的实验结果为例,cm-ssbl算法的相对重构误差为0.03左右,而omp-mmv,m-focuss,msbl算法的相对重构误差分别为0.075,0.06,0.05左右,换言之,在归一化采样率为0.4,测试图像是分辨率为1024×1024的“cameraman”的情况下,cm-ssbl算法与omp-mmv算法相比,精度提升60%左右,与m-focuss算法相比,精度提升50%左右,而与msbl算法相比,精度提升40%左右。
图5展示了以峰值信噪比作为评价标准的情况下,omp-mmv,m-focuss,msbl,以及本发明提出的cm-ssbl算法对分辨率为1024×1024的“cameraman”测试图像进行重构的效果,从图中可以看出,在归一化采样率从0.25到0.55的变化范围内,cm-ssbl算法与omp-mmv,m-focuss,msbl算法相比在峰值信噪比的性能方面同样有显著的提升。
为了从直观效果上对比omp-mmv,m-focuss,msbl,以及本发明提出的cm-ssbl算法的性能,我们采用分辨率为1024×1024的标准测试图像“cameraman”对上述三种方法进行图像重构测试,归一化采样率为0.4。从图6可以看出,cm-ssbl算法的重构效果要明显优于omp-mmv,m-focuss,以及msbl算法。同时,从客观的psnr角度进行评价,cm-ssbl方法相较于omp-mmv,m-focuss,msbl算法也提供了5到8db的增益。这说明,本发明提出的引入结构化先验约束的想法是切实有效的。
为了检验cm-ssbl算法的鲁棒性,我们选取了不同的测试图像对omp-mmv,m-focuss,msbl,以及cm-ssbl算法进行测试,分别统计相对重构误差和峰值信噪比并进行对照,实验结果如表1所示。表1是采用不同测试图像(所有测试图像的分辨率均为1024×1024,归一化采样率均为0.4)时本发明提出的图像重构方法(cm-ssbl)与已有的基于mmv模型的稀疏重构算法(omp-mmv、m-focuss和msbl)在处理大规模图像的稀疏重构问题时的性能对比。
表1针对1024×1024维标准测试图像的不同算法重构性能对比
表1展示了4种标准测试图像的测试结果,4种图像分别是“lena”,“boats”,“barbara”,以及“goldhill”,所有测试图像分辨率是1024×1024,归一化采样率是0.4。从表1可以直观地看出,cm-ssbl算法在不同的测试图像下与omp-mmv,m-focuss以及msbl算法相比,均能够显著提升重构性能。
从图7可以看出,迭代次数越多,重构的峰值信噪比越高,进一步地考虑到时间性能,我们在实施方案中优选出tmax=65作为仿真例中的迭代总次数。
以上结合附图对本发明的具体实施方式作了说明,但这些说明不能被理解为限制了本发明的范围,本发明的保护范围由随附的权利要求书限定,任何在本发明权利要求基础上的改动都是本发明的保护范围。