CN113222860A - 基于噪声结构多重正则化的图像恢复方法及系统 - Google Patents
基于噪声结构多重正则化的图像恢复方法及系统 Download PDFInfo
- Publication number
- CN113222860A CN113222860A CN202110614392.0A CN202110614392A CN113222860A CN 113222860 A CN113222860 A CN 113222860A CN 202110614392 A CN202110614392 A CN 202110614392A CN 113222860 A CN113222860 A CN 113222860A
- Authority
- CN
- China
- Prior art keywords
- image
- noise
- sparse
- parameters
- mean value
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 89
- 238000011084 recovery Methods 0.000 title claims abstract description 34
- 238000009826 distribution Methods 0.000 claims abstract description 34
- 230000006870 function Effects 0.000 claims abstract description 34
- 230000009466 transformation Effects 0.000 claims abstract description 22
- 238000002939 conjugate gradient method Methods 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 41
- 238000005259 measurement Methods 0.000 claims description 28
- 239000013598 vector Substances 0.000 claims description 28
- 238000005457 optimization Methods 0.000 claims description 10
- 238000004422 calculation algorithm Methods 0.000 abstract description 13
- 238000000844 transformation Methods 0.000 abstract description 3
- 238000004364 calculation method Methods 0.000 description 9
- 238000002474 experimental method Methods 0.000 description 6
- 230000014509 gene expression Effects 0.000 description 6
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 239000000126 substance Substances 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000011426 transformation method Methods 0.000 description 2
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 230000003190 augmentative effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000000701 chemical imaging Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000009629 microbiological culture Methods 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000001308 synthesis method Methods 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2136—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on sparsity criteria, e.g. with an overcomplete basis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N7/00—Computing arrangements based on specific mathematical models
- G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Artificial Intelligence (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Mathematical Analysis (AREA)
- Bioinformatics & Computational Biology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Probability & Statistics with Applications (AREA)
- Algebra (AREA)
- Computational Mathematics (AREA)
- Evolutionary Biology (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Computing Systems (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Image Processing (AREA)
Abstract
本发明公开了基于噪声结构多重正则化的图像恢复方法及系统，选择多个稀疏变换并加入相互独立的噪声，建立相互独立的稀疏表示似然函数；基于稀疏表示的似然函数，在各个稀疏域内，假设稀疏变换系数是相互独立的，各自建立满足共轭匹配关系的条件先验概率密度，定义服从Gamma分布的超参数；确定噪声不确定度和超参数的均值估计模型；最后通过变分期望最大化方法结合共轭梯度法迭代实现目标图像的快速重建，不但克服了综合法贝叶斯压缩感知对稀疏变换方法的单重性和可逆性限制，而且联合稀疏域能够有效提高算法的收敛速度。
Description
技术领域
本发明属于图像恢复领域，尤其涉及基于噪声结构多重正则化的图像恢复方法及系统。
背景技术
本部分的陈述仅仅是提供了与本发明相关的背景技术信息，不必然构成在先技术。
压缩感知(compressive sensing,CS)技术的主要原理是，对满足稀疏性要求的信号，以固定方式测量或采样，能够以远小于Nyquist-Shannon采样定律要求的采样数据量高精度地重建原信号。由于在重建稀疏信号的同时，能够实现超越常规方法的性能，被广泛应用于相关科学领域，如声呐和雷达成像、信号和图像处理等。根据CS原理，原问题相当于求解有约束条件的l0范数优化问题，该问题是NP-hard问题，无法通过凸优化方法求解。因此，常常采用贪婪算法近似计算，如正交匹配追踪(orthogonalmatchingpursuit，OMP)算法。另外，也可以采用lp(0＜p＜1)范数或l1范数降低求解难度，但是求解精度和重建成功率也会有所下降。
从贝叶斯理论角度看，CS方法通过引入目标信号的先验分布假设，结合Type-I求解方法，实现最大后验概率估计器，常用的岭回归、LASSO、重加权l1等CS方法都可以归纳和找到对应的先验分布概率密度函数，因此常规的CS方法也视为确定性CS，与统计型CS相对应。与Type-I方法对应，为了实现Type-Ⅱ最大似然求解，采用贝叶斯推算方法学习或训练多层次模型中的超参数。将条件后验概率密度函数作为真实后验的近似值，在已知超参数的前提下估计信号的均值和协方差。该计算过程要求模型中，相邻层次间满足共轭匹配关系，这样可以避免采用马尔可夫链蒙特卡罗(MCMC)进行低效率的迭代学习参数。采用分层的先验模型求解CS问题的方法——Type-Ⅱ方法，也称为贝叶斯压缩感知(Bayesiancompressive sensing,BCS)方法，在重建成功率和鲁棒性等方面优于确定性CS。BCS方法优势在于，除了提高稀疏信号的重建性能外，还不需要用户定义重要的正则化参数，并给出不确定性估计结果。根据报道，BCS相对于CS，在低采样率、信号高稀疏度和低信噪比的恶劣条件下具有更高的重建成功率和更低的归一化均方误差(Normalized Mean Square Error，NMSE)。因此，BCS被广泛应用于雷达成像、到达方向(DOA)估计、高光谱成像、阵列诊断、无线通信和图像恢复应用。
为提高待恢复图像的稀疏度，常采用稀疏表示(sparse representation,SR)方法或正则化方法来提升算法的输出图像质量和恢复成功率。常用的SR方法如全变分(TV)，非局部TV(NLTV)，离散余弦变换(DCT)，离散小波变换(DWT)，contourlet变换和shearlet变换。不同类型的图像，对SR方法的响应也不同。例如，TV在小尺度范围内提取图像的局部平滑性，而NLTV在大尺度上获取全局的相似度。contourlet和shearlet是图像处理中的高级SR方法，用于覆盖各向异性特征和纹理，特别适合于于二维(2D)和三维(3D)数据恢复。尽管这些SR方法是针对特定情况精心设计的，但很难在固定的单一稀疏变换域内对复杂信号进行稀疏化，例如自然图像和雷达/医学成像的重建目标。
针对上述问题，可以采用多重正则化结合CS(MRCS)来提高SR和信号重建的性能，MRCS方法大多基于交替方向乘子法(ADMM)，包括增广拉格朗日法、分裂Bregman迭代或优化-最小化(MM)方法，MRCS常常采用具有互补特征的SR组合来提高恢复性能，例如，用于图像处理和磁共振图像(MRI)的TV和NLTV组合，用于图像恢复的基于群的稀疏表示(GSR)和NLTV，加权核范数和TV用于高光谱图像去噪。
然而，BCS的多重正则化扩展却不如CS那样直接，因为必须在考虑先验概率时，同时建立共轭匹配关系，综合法能够避免上述考虑，但是对稀疏变换有单重性和可逆性要求，这里无法使用。因此，目前尚未有理论基础完整、适用性广泛的多重正则化BCS方法出现。
发明内容
本发明为了解决上述问题，提出了基于噪声结构多重正则化的图像恢复方法及系统，通过高斯白噪声结构引入多个稀疏变换，即通过引入相互独立的多个满足高斯分布的概率密度函数，根据贝叶斯原理，建立多重正则化多层次贝叶斯模型，进一步对各重变换对应的变换系数给出稀疏性假设，即满足均值为零的高斯分布，从而建立各层次之间共轭匹配关系，最后通过变分期望最大化方法结合共轭梯度法迭代实现目标图像的快速重建，不但克服了综合法贝叶斯压缩感知对稀疏变换方法的单重性和可逆性限制，而且联合稀疏域能够有效提高算法的收敛速度。
根据一些实施例，本发明采用如下技术方案：
基于噪声结构多重正则化的图像恢复方法，包括：
获取待恢复图像，并将图像转化为一维信号后，用测量矩阵投影到观测向量；
将测量矩阵和观测向量输入基于噪声结构多重正则化的图像恢复模型，得到目标图像均值；其中，基于噪声结构多重正则化的图像恢复模型通过超参数、稀疏变换系数和噪声的不确定度的迭代更新实现目标图像均值优化。
基于噪声结构多重正则化的图像恢复系统，包括：
图像获取模块，其用于获取待恢复图像，并将图像转化为一维信号后，用测量矩阵投影到观测向量；
图像恢复模块，其用于将测量矩阵和观测向量输入基于噪声结构多重正则化的图像恢复模型，得到目标图像均值；其中，基于噪声结构多重正则化的图像恢复模型通过超参数、稀疏变换系数和噪声的不确定度的迭代更新实现目标图像均值优化。
一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成第一方面所述方法的步骤。
一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成第一方面所述方法的步骤。
本发明的有益效果是：
1、本发明通过联合多重正则化方法来提高对目标图像的稀疏表示能力，采用贝叶斯压缩感知方法获得比传统压缩感知方法更强的稀疏重建性能，并有效提高重建成功率。
2、本发明提出基于噪声结构的贝叶斯多层次模型，不但克服了综合法贝叶斯压缩感知对稀疏变换方法的单重性和可逆性限制，而且在联合稀疏域能够有效提高算法的收敛速度。
3、本发明为提高算法对大尺寸图像的重建效率，提出采用共轭梯度法迭代估计目标图像均值，并利用迭代过程中间参数估计的超参数和噪声方差，避免了协方差矩阵的求逆运算，有效降低了内存需求和计算量。
4、本发明的基于噪声结构多重正则化的图像恢复方法适用于低信噪比图像的恢复，由于采用了高斯噪声结构，在低信噪比情况下，具有更高的重建精度和稳定性。
附图说明
构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
图1是本发明基于噪声结构的多重正则化多层次模型结构图；
图2是分段稀疏信号举例；
图3是分段稀疏信号恢复不确定性性能曲线；
图4是分段稀疏信号恢复稀疏度性能曲线；
图5是分段稀疏信号恢复信噪比性能曲线；
图6是用于图像恢复的原始图像：生物细胞；
图7是用于图像恢复的原始图像：电子线路图；
图8是用于图像恢复的原始图像：月亮；
图9是用于图像恢复的原始图像：汽车轮胎；
图10(a)是SNR＝20dB时，采用GN-MRSBL生物细胞图像恢复结果；
图10(b)是SNR＝20dB时，采用ADMM生物细胞图像恢复结果；
图10(c)是SNR＝20dB时，采用proximal方法生物细胞图像恢复结果；
图11(a)是SNR＝20dB时，采用GN-MRSBL电子线路图像恢复结果；
图11(b)是SNR＝20dB时，采用ADMM电子线路图像恢复结果；
图11(c)是SNR＝20dB时，采用proximal方法电子线路图像恢复结果；
图12(a)是SNR＝20dB时，采用GN-MRSBL月亮图像恢复结果；
图12(b)是SNR＝20dB时，采用ADMM月亮图像恢复结果；
图12(c)是SNR＝20dB时，采用proximal方法月亮图像恢复结果；
图13(a)是SNR＝20dB时，采用GN-MRSBL汽车轮胎图像恢复结果；
图13(b)是SNR＝20dB时，采用ADMM汽车轮胎图像恢复结果；
图13(c)是SNR＝20dB时，采用proximal方法汽车轮胎图像恢复结果；
图14(a)是SNR＝40dB时，采用GN-MRSBL生物细胞图像恢复结果；
图14(b)是SNR＝40dB时，采用ADMM生物细胞图像恢复结果；
图14(c)是SNR＝40dB时，采用proximal方法生物细胞图像恢复结果；
图15(a)是SNR＝40dB时，采用GN-MRSBL电子线路图像恢复结果；
图15(b)是SNR＝40dB时，采用ADMM电子线路图像恢复结果；
图15(c)是SNR＝40dB时，采用proximal方法电子线路图像恢复结果；
图16(a)是SNR＝40dB时，采用GN-MRSBL月亮图像恢复结果；
图16(b)是SNR＝40dB时，采用ADMM月亮图像恢复结果；
图16(c)是SNR＝40dB时，采用proximal方法月亮图像恢复结果；
图17(a)是SNR＝40dB时，采用GN-MRSBL汽车轮胎图像恢复结果；
图17(b)是SNR＝40dB时，采用ADMM汽车轮胎图像恢复结果；
图17(c)是SNR＝40dB时，采用proximal方法汽车轮胎图像恢复结果。
具体实施方式
下面结合附图与实施例对本发明作进一步说明。
应该指出，以下详细说明都是例示性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
实施例1
如图1所示，本实施例提供的基于噪声结构多重正则化的图像恢复方法，通过高斯白噪声结构引入多个稀疏变换，即通过引入相互独立的多个满足高斯分布的概率密度函数，根据贝叶斯原理，建立多重正则化多层次贝叶斯模型，进一步对各重变换对应的变换系数给出稀疏性假设，即满足均值为零的高斯分布，从而建立各层次之间共轭匹配关系，最后通过变分期望最大化(variational expectation maximization,VEM)方法结合共轭梯度法(conjugate gradient method,CGM)迭代实现目标图像的快速重建。图像恢复是利用图像中已知区域信息对破损区域进行信息填充，以弥补信息的损失，从而获得目标图像。
基于噪声结构多重正则化的图像恢复方法，包括如下步骤：
步骤(1)：通过测量矩阵(或隐函数形式的测量函数)和观测向量，建立图像恢复的线性方程；即，将测量方程写成显式的线性方程或者隐式函数方程，前者的形式为，观测向量等于测量矩阵与一维向量化待重建目标图像的乘积，后者直接给出测量函数和测量函数的共轭转置函数，加入高斯噪声，构建满足高斯分布的测量似然函数。在本发明中，能够采用矩阵形式表达的函数称为显式函数，例如傅里叶变换和全变差变换；由于维度限制或者非线性关系，难以采用矩阵形式表达，而采用函数输入-输出形式表达的，称为隐函数，例如剪切波变换等。
具体为：假设待恢复目标图像表示为二维矩阵X，通过将X的所有列向量拼接为一个列向量x，假设
(即待恢复图像向量的维度为N)，可以建立观测向量
(即观测向量的维度为M)与目标图像x的线性关系
y＝Ax+n (1)
其中，
为测量矩阵，
为复高斯白噪声，一般满足M＜N。
或者给出测量矩阵的隐式表达式
y＝A(X)+n (2)
并同时给出测量函数的共轭转置函数AH(·)，因此，y的似然函数表示为
其中，噪声的不确定度β＝σ-2，其中σ-2为噪声n的方差。
步骤(2)：通过高斯白噪声结构引入多个稀疏变换Dj，建立相互独立的稀疏表示似然函数。
根据先验信息，选择合理的稀疏表示方法组合，建立相互独立的稀疏表示似然函数方程，满足高斯分布；即，根据所选择的多个稀疏变换方法，给出稀疏变换的显式变换矩阵或隐式变换函数及对应变换函数的共轭转置变换函数，加入相互独立的高斯白噪声，建立稀疏表示的似然函数方程。
具体的：选择J个复数域稀疏变换以Dj表示，其中
Nj为第j个变换域稀疏系数的维度，j＝1，2，...，J，因此变换系数满足噪声方程
nj＝sj-Djx (4)
其中
nj为零均值高斯白噪声，即
Θj＝diag(θj)，其中diag(θj)表示以向量θj为对角阵元素的矩阵。因此得到稀疏表示似然函数：
步骤(3)：基于稀疏表示的似然函数，在各个稀疏域内，假设稀疏变换系数是相互独立的，各自建立满足共轭匹配关系的先验条件概率密度函数，定义服从Gamma分布的超参数αj、θj，并定义服从Gamma分布的噪声不确定度β。
定义超参数的先验分布参数，构建贝叶斯模型各层次的条件先验概率密度；即，在各个稀疏域内，假设稀疏表示系数是相互独立的，并且各自建立满足共轭匹配关系的条件先验概率密度函数，定义其中Gamma分布的参数，鉴于更深层次的参数分布对推算结果影响有限，分布参数定义为常数即可；
具体的：
首先，定义公式(3)中噪声不确定度β满足Gamma分布，即
p(β)＝Gamma(β；c，d) (6)
其中，c和d为预先定义的常数参数。
根据公式(5)，条件先验概率密度p(x|sG；θG)表示为
其中，下标G表示稀疏的标号集，即G＝{1，2，...，J}。
在先验概率密度p(x)服从无信息假设(即均匀分布)时，得到p(x|sG；θG)也服从高斯分布，即
为了给稀疏系数施加稀疏性先验条件并保留层次之间的共轭匹配性，假设
其中，Λj＝diag(αj)，超参数αj是独立同分布的，而且服从Gamma分布
其中，aj和bj为较小的常数参量。
同理，超参数θj也满足独立同分布的Gamma分布
其中，ej和fj为较小的常数参量。
步骤(4)：确定噪声不确定度和超参数的均值表达式，基于噪声不确定度和超参数的均值表达式，对待修复图像进行目标图像均值估计。
通过VEM算法完成所有超参数的条件后验概率密度近似估计，并根据估计均值得到算法的迭代公式，结合CGM实现目标图像均值估计。即，根据VEM原理，所有超参数满足条件后验概率密度函数相互独立假设，以此推算各超参数满足的近似后验概率密度分布，根据超参数的均值表达式，利用CGM算法实现图像均值的迭代求解。
根据VEM理论，采用变分概率分布q(x，sG，θG，αG，β)来近似后验概率分布p(x，sG，θG，αG，β|y)，而且可以分解为如下相互独立的后验概率分布
在M-step中，为实现全概率密度最大化，即
只需分别估计qx(x)，qs(sG)，qθ(θG)，qα(αG)和qβ(β)的分布。
步骤(4.1)：推算目标图像满足的近似后验概率分布，采用CGM迭代计算待重建目标图像的均值，并提取CGM迭代中的向量参数；
步骤(4.2)：根据CGM向量参数与图像协方差之间的关系，计算中间参数，用于后续的参数重建，而不必完整记录图像协方差矩阵；
步骤(4.1)和步骤(4.2)的步骤为：
首先，给出目标图像x的近似条件后验概率密度满足
其中，
表示关于x的概率密度函数qx(x)的均值。
显然，qx(x)满足高斯分布，其均值和协方差矩阵为
其中，均值μx即为估计均值。由于在后续超参数更新中，多次用到协方差矩阵∑x，因此一般采用直接计算公式(15)中的矩阵求逆，但是对于大尺度信号，如512×512像素图像，直接求逆不但在计算量上无法接受，而且常规计算机内存无法实现。
本发明采用CGM方法迭代计算均值μx和与协方差矩阵∑x相关的中间参数，方法如下：
引入矩阵
公式(14)可以写成线性厄米特方程
根据CGM原理，将公式(16)写成线性方程形式Bx＝b，其中B＝Lx，x＝μx，
则求解步骤为
(1)初始化：x(0)＝BHb，r(0)＝b-Bx (0)和w(0)＝r(0)；
(2)对于第k次迭代，完成如下步骤：
x(k+1)＝x(k)+ζ(k)w(k)
r(k+1)＝r(k)-ζ(k)Bw(k)
w(k+1)＝r(k+1)+β(k)w(k)
其中，上角标k表示迭代次数，每次迭代的列向量组成矩阵Wk＝[Wk-1w(k)]满足Lx共轭关系，即
因此，协方差矩阵
可以表示为
∑x＝∑kw(k)(w(k))H (17)
于是，∑x的对角线元素可以表示为
diag(∑x)＝∑k|w(k)|2 (18)
引入两个辅助变量
h(k)＝Aw(k) (20)
得到
和
diag(A∑xAH)＝∑k|h(k)|2 (22)
最后，计算中间参数，即中间参数
tr(A∑xAH)＝sum(diag(A∑xAH)) (23)
步骤(4.3)：利用步骤(4.2)计算得到的中间参数，计算所有超参数和噪声方差的估计均值；具体的步骤为：
根据稀疏变换系数sj的相互独立假设，其近似条件后验概率密度满足
可见，qs(sj)满足高斯分布，稀疏变换系数sj的均值和协方差可以表示为
其中，<·>表示估计均值。
下面推算θj的更新公式：θj满足Gamma先验分布，即
其中，θj(m)表示矢量θj的第m个元素，因此θj的条件后验概率密度qθ(θj)满足
由上述结果可见，qθ(θj)也是Gamma分布
其中，
表示矩阵
的第m行第m列元素，对所有元素m，相当于计算
由于
为对角阵，因此这里的计算容易实现；当计算
时，需要用到CGM算法中的中间参数，如公式(21)所示可以实现快速计算。
下面推算αj的更新公式：αj满足如下先验概率密度
于是近似条件后验概率密度函数qα(αj)满足
类似于θj，超参数αj也服从Gamma分布
最后，推算β的更新公式：条件后验概率密度函数qβ(β)满足
因此
其中，tr(A∑xA)也需要采用CGM中的中间参数，如公式(23)所示为相应的快速计算方法。
总结上述过程，需要迭代求解如下参数的均值
步骤(4.4)：判断迭代收敛性，若满足收敛条件中的一种，结束迭代输出图像均值的估计结果；若不满足，继续迭代。
其中，收敛条件包括：
(1)
其中tol为容差；
(2)
容差tol可与条件(1)设相同值，其中
(3)迭代次数k＞预设最大迭代次数Imax。
本申请的基于噪声结构多重正则化的图像恢复方法，接收待恢复图像后，将图像转化为一维信号后，用测量矩阵A投影到观测向量y；将测量矩阵和观测向量输入基于噪声结构多重正则化的图像恢复模型，得到目标图像均值；基于噪声结构多重正则化的图像恢复模型通过超参数、稀疏变换系数和噪声的不确定度的迭代更新实现目标图像均值优化。
基于噪声结构多重正则化的图像恢复模型通过超参数、稀疏变换系数和噪声的不确定度的迭代更新实现目标图像均值优化的具体步骤包括：
(1)输入：稀疏变换Dj，矢量形式y、测量矩阵A；
(2)初始化：最大迭代次数Imax，容差tol，定义aj＝bj＝ej＝fj＝c＝d为较小的数值(如10-9)，j＝1，2，...，J，超参数初始值α、θ，噪声不确定度β，稀疏变换稀疏均值<sj>；
(3)基于测量矩阵、观测向量、超参数的均值、稀疏变换系数的均值、噪声不确定度的均值，采用CGM计算目标图像的均值μx和与协方差矩阵∑x，并根据公式(21)-(23)得到相关的中间参数，具体的：
目标图像x的均值和协方差矩阵的表达式为
引入矩阵
目标图像x的均值的表达式可以写成线性厄米特方程
根据CGM原理，将公式线性厄米特方程写成线性方程形式Bx＝b，其中B＝Lx，x＝μx，
则求解步骤为
(a)初始化：x(0)＝BHb，r(0)＝b-Bx(0)和w(0)＝r(0)；
(b)对于第k次迭代，完成如下步骤：
x(k+1)＝x(k)+ζ(k)w(k)
r(k+1)＝r(k)-ζ(k)Bw(k)
w(k+1)＝r(k+1)+β(k)w(k)
其中，上角标k表示迭代次数，每次迭代的列向量组成矩阵Wk＝[Wk-1w(k)]满足Lx共轭关系，即
协方差矩阵
表示为
∑x＝∑kw(k)(w(k))H
∑x的对角线元素表示为
diag(∑x)＝∑k|w(k)|2
引入两个辅助变量
h(k)＝Aw(k)
得到
和
diag(A∑xAH)＝∑k|h(k)|2
最后，计算中间参数
tr(A∑xAH)＝sum(diag(A∑xAH))。
(4)基于中间参数、目标图像x的均值、稀疏变换系数sj的均值和协方差，根据公式(27)、(28)、(36)计算超参数θj的均值，具体的：
基于稀疏变换系数sj的均值和协方差，根据公式(30)、(31)、(37)计算超参数αj的均值，具体的：
(5)基于超参数θj的均值、超参数αj的均值、目标图像x的均值，结合公式(24)、(25)，计算稀疏变换系数sj的均值和协方差，具体的：
(6)基于中间参数、目标图像x的均值，根据公式(34)、(35)、(38)计算噪声不确定度β的均值<β>，具体的：
(7)判断迭代收敛性，即判断是否满足收敛条件之一，如不满足，返回步骤(3)；若满足，结束迭代，输出目标图像的均值，得到最终的目标图像，即恢复后的图像。
在实验中，待恢复的分段稀疏信号由三部分组成，分别对应恒等变换、一阶TV和Haar小波变换稀疏性，一个信号实例如图2所示。特别地，我们将分段稀疏信号的稀疏性定义为所有本征稀疏信号的稀疏度和，用K表示。测量矩阵A的列服从标准正态分布，各列被归一化为单位范数。每个实验重复200次，在无噪声实验中，确定度定义为M/N，稀疏度定义为K/N。信号的恢复成功定义为，当归一化均方误差(NMSE)满足预设准则NMSE＜10-3时，信号被认为是成功恢复。NMSE也是噪声实验的测量性能指标，定义为
其中，xgen是待恢复的原始信号，
是方法恢复的估计结果。
采用的用于比对的方法包括两种多重正则化CS方法：ADMM和proximal方法；附图中，本申请的基于噪声结构多重正则化的图像恢复方法被简称为GN-MRSBL。从图3和图4的结果可以看出，GN-MRSBL在不确定性和稀疏度方面的成功率最高，显著优于ADMM和proximal方法。在多约束正则化的情况下，采用proximal方法分段稀疏信号的NMSE普遍大于0.6，这在本实验中被认为是恢复失败。为了清晰地展示信噪比的恢复性能，在图5中采用NMSE的对数形式(Log-NMSE)。从结果看，在中等信噪比条件下，即20-35dB，GN-MRSBL似乎优势较弱，而在低信噪比条件下有明显优势。
为了验证多重正则化的有效性，采用多重正则化方法验证图像数据的恢复性能，分别采用恒等变换、一阶TV和非局部TV的联合稀疏表示。如图6-9所示，采用四幅图像作为原始数据，包括生物细胞、电子电路图、月球和汽车轮胎。图像的大小是[256，256]。分别采用GN-MRSBL、ADMM和proximal方法3种方法用于图像恢复。测量矩阵A采用置乱块Hadamard集合生成，测量的不确定度比为0.625。实验中，分别在信噪比(SNR)为20dB和40dB条件下比较图像恢复性能，恢复结果如图10-17所示。采用NMSE和结构相似度(SSIM)指标定量衡量图像恢复性能，在信噪比为20dB和40dB条件下的结果分别总结在表1和2中。对于每一张图像，性能最好的用粗体突出显示。从恢复结果看，在信噪比为20～40dB之间，ADMM的结果过于平滑，噪声对其影响较小。从两种性能指标NMSE和SSIM衡量，GN-MRSBL在信噪比为20dB和40dB时均能达到最佳性能。
表1信噪比为20dB时NMSE/SSIM估计结果
表2信噪比为40dB时NMSE/SSIM估计结果
本发明通过联合多重正则化方法来提高对目标图像的稀疏表示能力，采用贝叶斯压缩感知方法获得比传统压缩感知方法更强的稀疏重建性能，并有效提高重建成功率；基于噪声结构的贝叶斯多层次模型，不但克服了综合法贝叶斯压缩感知对稀疏变换方法的单重性和可逆性限制，而且在联合稀疏域能够有效提高算法的收敛速度；最后，采用共轭梯度法迭代估计目标图像均值，并利用迭代过程中间参数估计的超参数和噪声方差，避免了协方差矩阵的求逆运算，有效降低了内存需求和计算量。
实施例2
本实施例提供基于噪声结构多重正则化的图像恢复系统，包括：
图像获取模块，其用于获取待恢复图像，并将图像转化为一维信号后，用测量矩阵投影到观测向量；
图像恢复模块，其用于将测量矩阵和观测向量输入基于噪声结构多重正则化的图像恢复模型，得到目标图像均值；其中，基于噪声结构多重正则化的图像恢复模型通过超参数、稀疏变换系数和噪声的不确定度的迭代更新实现目标图像均值优化。
实施例3
本实施例还提供了一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成实施例1所述方法的步骤。
实施例4
本实施例还提供了一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成实施例1所述方法的步骤。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.基于噪声结构多重正则化的图像恢复方法，其特征在于，包括：
获取待恢复图像，并将图像转化为一维信号后，用测量矩阵投影到观测向量；
将测量矩阵和观测向量输入基于噪声结构多重正则化的图像恢复模型，得到目标图像均值；其中，基于噪声结构多重正则化的图像恢复模型通过超参数、稀疏变换系数和噪声的不确定度的迭代更新实现目标图像均值优化。
2.如权利要求1所述的基于噪声结构多重正则化的图像恢复方法，其特征在于，所述基于噪声结构多重正则化的图像恢复模型通过超参数、稀疏变换系数和噪声的不确定度的迭代更新实现目标图像均值优化的具体步骤为：
基于测量矩阵和观测向量，采用共轭梯度法计算目标图像均值和中间参数；
基于中间参数，更新超参数、稀疏变换系数和噪声的不确定度；
判断是否满足收敛条件，若满足，输出目标图像均值；否则，采用更新后的超参数、稀疏变换系数和噪声的不确定度，返回重新计算目标图像均值和中间参数。
3.如权利要求2所述的基于噪声结构多重正则化的图像恢复方法，其特征在于，所述基于中间参数，更新超参数、稀疏变换系数和噪声的不确定度的具体步骤为：
基于中间参数、稀疏变换系数的均值和协方差，确定超参数的均值；
基于超参数的均值，更新稀疏变换系数的均值和协方差；
基于中间参数，计算噪声不确定度的均值。
4.如权利要求1所述的基于噪声结构多重正则化的图像恢复方法，其特征在于，所述超参数的更新基于超参数满足条件后验概率密度函数相互独立假设。
5.如权利要求1所述的基于噪声结构多重正则化的图像恢复方法，其特征在于，所述噪声的不确定度服从Gamma分布。
6.如权利要求1所述的基于噪声结构多重正则化的图像恢复方法，其特征在于，所述变换系数满足均值为零的高斯分布。
7.如权利要求1所述的基于噪声结构多重正则化的图像恢复方法，其特征在于，所述超参数是独立同分布的，且服从Gamma分布。
8.基于噪声结构多重正则化的图像恢复系统，其特征在于，包括：
图像获取模块，其用于获取待恢复图像，并将图像转化为一维信号后，用测量矩阵投影到观测向量；
图像恢复模块，其用于将测量矩阵和观测向量输入基于噪声结构多重正则化的图像恢复模型，得到目标图像均值；其中，基于噪声结构多重正则化的图像恢复模型通过超参数、稀疏变换系数和噪声的不确定度的迭代更新实现目标图像均值优化。
9.一种电子设备，其特征是，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成权利要求1-6任一项方法的步骤。
10.一种计算机可读存储介质，其特征是，用于存储计算机指令，所述计算机指令被处理器执行时，完成权利要求1-6任一项方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110614392.0A CN113222860B (zh) | 2021-06-02 | 2021-06-02 | 基于噪声结构多重正则化的图像恢复方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110614392.0A CN113222860B (zh) | 2021-06-02 | 2021-06-02 | 基于噪声结构多重正则化的图像恢复方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113222860A true CN113222860A (zh) | 2021-08-06 |
| CN113222860B CN113222860B (zh) | 2022-09-20 |
Family
ID=77082329
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110614392.0A Active CN113222860B (zh) | 2021-06-02 | 2021-06-02 | 基于噪声结构多重正则化的图像恢复方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113222860B (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113421208A (zh) * | 2021-08-24 | 2021-09-21 | 北京与光科技有限公司 | 高分辨率光谱恢复方法 |
| CN114296087A (zh) * | 2021-12-13 | 2022-04-08 | 哈尔滨工程大学 | 一种线上贝叶斯压缩水下成像方法及系统、设备和介质 |
| CN116452425A (zh) * | 2023-06-08 | 2023-07-18 | 常州星宇车灯股份有限公司 | 图像超分辨率重建方法、设备及介质 |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100246920A1 (en) * | 2009-03-31 | 2010-09-30 | Iowa State University Research Foundation, Inc. | Recursive sparse reconstruction |
| CN105550993A (zh) * | 2016-01-18 | 2016-05-04 | 中国空间技术研究院 | 基于多重变换域的超分辨率重建方法 |
| CN205642634U (zh) * | 2016-02-04 | 2016-10-12 | 北京理工大学 | 基于压缩感知的双波长温度场成像设备及系统 |
| WO2018099321A1 (zh) * | 2016-11-30 | 2018-06-07 | 华南理工大学 | 一种基于广义树稀疏的权重核范数磁共振成像重建方法 |
| US20180198992A1 (en) * | 2016-08-30 | 2018-07-12 | Shenzhen University | Infrared image reconstruction method based on block-sparse compressive sensing and system thereof |
| CN109064406A (zh) * | 2018-08-26 | 2018-12-21 | 东南大学 | 一种正则化参数自适应的稀疏表示图像重建方法 |
| CN109447921A (zh) * | 2018-12-05 | 2019-03-08 | 重庆邮电大学 | 一种基于重构误差的图像测量矩阵优化方法 |
| CN110852973A (zh) * | 2019-11-12 | 2020-02-28 | 华中科技大学 | 一种脉冲噪声模糊图像非线性复原方法及系统 |
| CN111640080A (zh) * | 2020-06-05 | 2020-09-08 | 郑州轻工业大学 | 一种基于高光谱全变差的cs图像去噪重建方法 |
| CN112529826A (zh) * | 2020-12-10 | 2021-03-19 | 南京航空航天大学 | 截断式张量贝叶斯多光谱图像压缩感知重构方法 |
| CN112541965A (zh) * | 2020-12-02 | 2021-03-23 | 国网重庆市电力公司电力科学研究院 | 基于张量逼近和时空相关性的压缩感知图像和视频恢复 |
-
2021
- 2021-06-02 CN CN202110614392.0A patent/CN113222860B/zh active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100246920A1 (en) * | 2009-03-31 | 2010-09-30 | Iowa State University Research Foundation, Inc. | Recursive sparse reconstruction |
| CN105550993A (zh) * | 2016-01-18 | 2016-05-04 | 中国空间技术研究院 | 基于多重变换域的超分辨率重建方法 |
| CN205642634U (zh) * | 2016-02-04 | 2016-10-12 | 北京理工大学 | 基于压缩感知的双波长温度场成像设备及系统 |
| US20180198992A1 (en) * | 2016-08-30 | 2018-07-12 | Shenzhen University | Infrared image reconstruction method based on block-sparse compressive sensing and system thereof |
| WO2018099321A1 (zh) * | 2016-11-30 | 2018-06-07 | 华南理工大学 | 一种基于广义树稀疏的权重核范数磁共振成像重建方法 |
| CN109064406A (zh) * | 2018-08-26 | 2018-12-21 | 东南大学 | 一种正则化参数自适应的稀疏表示图像重建方法 |
| CN109447921A (zh) * | 2018-12-05 | 2019-03-08 | 重庆邮电大学 | 一种基于重构误差的图像测量矩阵优化方法 |
| CN110852973A (zh) * | 2019-11-12 | 2020-02-28 | 华中科技大学 | 一种脉冲噪声模糊图像非线性复原方法及系统 |
| CN111640080A (zh) * | 2020-06-05 | 2020-09-08 | 郑州轻工业大学 | 一种基于高光谱全变差的cs图像去噪重建方法 |
| CN112541965A (zh) * | 2020-12-02 | 2021-03-23 | 国网重庆市电力公司电力科学研究院 | 基于张量逼近和时空相关性的压缩感知图像和视频恢复 |
| CN112529826A (zh) * | 2020-12-10 | 2021-03-19 | 南京航空航天大学 | 截断式张量贝叶斯多光谱图像压缩感知重构方法 |
Non-Patent Citations (6)
| Title |
|---|
| CHAO LI 等: "Compressive Beamforming Based on Multiconstraint Bayesian Framework", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 * |
| NICOLAS DOBIGEON 等: "REGULARIZED BAYESIAN COMPRESSED SENSING IN ULTRASOUND IMAGING", 《20TH EUROPEAN SIGNAL PROCESSING CONFERENCE (EUSIPCO 2012)》 * |
| SANJAY VISWANATH, MUTHUVEL ARIGOVINDAN: "Structurally Adaptive Multi-Derivative Regularization for Image Recovery from Sparse Fourier Samples", 《HTTPS://ARXIV.ORG/ABS/2105.12775V1》 * |
| 党宏社，张娜: "一种基于贝叶斯压缩感知的图像修复方法", 《河南大学学报（自然科学版）》 * |
| 刘梅 等: "一种基于神经网络的超分辨图像重构方法", 《哈尔滨工业大学学报》 * |
| 郭从洲 等: "空间目标图像的非凸稀疏正则化波后复原", 《光学精密工程》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113421208A (zh) * | 2021-08-24 | 2021-09-21 | 北京与光科技有限公司 | 高分辨率光谱恢复方法 |
| CN113421208B (zh) * | 2021-08-24 | 2021-11-02 | 北京与光科技有限公司 | 高分辨率光谱恢复方法 |
| CN114296087A (zh) * | 2021-12-13 | 2022-04-08 | 哈尔滨工程大学 | 一种线上贝叶斯压缩水下成像方法及系统、设备和介质 |
| CN114296087B (zh) * | 2021-12-13 | 2024-07-12 | 哈尔滨工程大学 | 一种线上贝叶斯压缩水下成像方法及系统、设备和介质 |
| CN116452425A (zh) * | 2023-06-08 | 2023-07-18 | 常州星宇车灯股份有限公司 | 图像超分辨率重建方法、设备及介质 |
| CN116452425B (zh) * | 2023-06-08 | 2023-09-22 | 常州星宇车灯股份有限公司 | 图像超分辨率重建方法、设备及介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113222860B (zh) | 2022-09-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113222860B (zh) | 基于噪声结构多重正则化的图像恢复方法及系统 | |
| CN109116293B (zh) | 一种基于离格稀疏贝叶斯的波达方向估计方法 | |
| US8760572B2 (en) | Method for exploiting structure in sparse domain for magnetic resonance image reconstruction | |
| CN105741333B (zh) | 一种Video-SAR图像实时压缩重构方法 | |
| CN111047661B (zh) | 一种基于稀疏流形联合约束的cs-mri图像重构方法 | |
| CN109887050B (zh) | 一种基于自适应字典学习的编码孔径光谱成像方法 | |
| CN108765313B (zh) | 基于类内低秩结构表示的高光谱图像去噪方法 | |
| CN111754598B (zh) | 基于变换学习的局部空间邻域并行磁共振成像重构方法 | |
| Cao et al. | CS-MRI reconstruction based on analysis dictionary learning and manifold structure regularization | |
| CN109584330A (zh) | 一种基于压缩感知逼近l0范数的梯度投影图像重建方法 | |
| CN114624646B (zh) | 一种基于模型驱动复数神经网络的doa估计方法 | |
| Zhang et al. | LR-CSNet: low-rank deep unfolding network for image compressive sensing | |
| CN109920017B (zh) | 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法 | |
| Chaabouni et al. | Using ROI with ISOM compression to medical image | |
| CN111161184A (zh) | 一种基于mcp稀疏约束的快速mr图像去噪方法 | |
| Wen et al. | The power of complementary regularizers: Image recovery via transform learning and low-rank modeling | |
| CN109188327B (zh) | 基于张量积复小波紧框架的磁共振图像快速重构方法 | |
| Shastri et al. | Denoising generalized expectation-consistent approximation for MR image recovery | |
| CN113222861B (zh) | 基于等式结构多重正则化的图像恢复方法及系统 | |
| CN115471580A (zh) | 一种物理智能高清磁共振扩散成像方法 | |
| CN111397733B (zh) | 一种单/多帧快照式光谱成像方法、系统及介质 | |
| CN113487507A (zh) | 一种基于多模块特征聚合的双域递归网络mr重建方法 | |
| Chen et al. | Robust image compressive sensing based on m-estimator and nonlocal low-rank regularization | |
| Gao et al. | High-quality image interpolation via local autoregressive and nonlocal 3-D sparse regularization | |
| Xie et al. | A probabilistic model-based method with nonlocal filtering for robust magnetic resonance imaging reconstruction |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |