一种基于稀疏表示的RAW图像去噪方法与流程

本发明属于图像处理技术领域，特别涉及一种基于稀疏表示的RAW图像去噪方法。

背景技术：

微小卫星相对于大卫星，具有质量轻，成本低，设计与制造周期短等特点。对地观测作为一项空间应用的重要任务，有着广泛的应用，如植被监测和灾难监测。在该领域，微小卫星已经成功得到应用。图像处理系统是对地观测任务的关键组成部分。

但是，由于复杂的太空环境，相机拍摄的图像受到了严重的噪声影响。例如光线强度较弱将会带来很大的噪声，相机通过增加输入信号的感光度(增加ISO)来放大感光元件的电压，会造成噪声增强。数字相机中的图像前处理技术在原本噪声源的基础上还会带来特定的方法噪声，例如白平衡、色彩校正、伽马校正和边缘增强会改变噪声源的性质或增加噪声的强度。如果直接压缩图像，传输到地面处理，那么方法噪声、量化噪声和传输过程的冲击噪声会使得噪声的模型更加复杂。

目前传统的信号稀疏表示去噪方法对已知噪声模型(如加性高斯白噪声)的灰度或者彩色图像去噪。自然信号在冗余字典中稀疏分解，是一种有效的信号表示方法。在图像去噪领域，稀疏表示很好地将信号分解成稀疏向量，通过稀疏向量恢复无噪图像。

但实际的问题是，经过白平衡、去马赛克、颜色变换三种主要的图像前处理技术，数字相机输出YUV或者RGB图像噪声模型是未知的。而且图像具有信号相关、颜色相关和空间相关的噪声特性。如果简单地利用传统稀疏表示方法去噪，那么就会出现去噪效果差的问题。另一方面，原始图像(或称为：RAW图像)的噪声模型是已知的，即加性高斯白噪声。RAW图像的格式一般采用Bayer颜色滤波阵列，即图像中像素总数的一半是绿色(G)，而总数的四分之一分别为红色(R)和蓝色(B)，像素首先按GRGR排列，下一行像素按BGBG排列。如果直接处理RAW图像，相应的传统稀疏表示图像处理方法有限制，因为RAW图像数据不符合一般的图像假设，如：平滑性、分段不变性等。

技术实现要素：

本发明的目的是提供一种基于稀疏表示的RAW图像去噪方法，以解决在处理未去噪RAW图像过程中，不能有效地估计出噪声的参数，计算复杂的问题。

本发明采用以下技术方案：一种基于稀疏表示的RAW图像去噪方法，其特征在于，包括以下步骤：

步骤一、将待去噪RAW图像分解成为固定大小相互重叠的RAW矩形图像块；

步骤二、将步骤一中得到的RAW矩形图像块按照不同颜色通道，重新排列成G1RBG2颜色层；

步骤三、将步骤二中得到的G1RBG2颜色层按照G1、R、B、G2的顺序依次拉伸后得到列向量形式的图像块向量，将每个图像块向量从左至右拼接成G1RBG2矩阵，G1RBG2矩阵的每一列为一个训练样本；

步骤四、从步骤三中的G1RBG2矩阵中随机选择部分训练样本作为训练样本集T，采用K-SVD方法对训练样本集T进行学习得出字典D；

步骤五、采用改进后的OMP算法对步骤四中得出的字典D进行图像重建后，形成无噪RAW图像。

进一步地，步骤四中采用K-SVD方法对训练样本集T进行学习得出字典D的具体步骤如下：

阶段一：初始化阶段，即从训练样本集T中随机选取i个n维样本t_i用来组成字典D，并始字典D进行J’次迭代，J’为迭代总次数；

阶段二：固定字典D阶段，即采用OMP算法找出训练样本集T中的每个训练样本的稀疏表示α_i；

阶段三：更新字典D阶段，即将固定字典D阶段得到的每个训练样本的稀疏表示组成稀疏矩阵A，通过稀疏矩阵A对阶段二中得出的字典D进行处理，依次更新字典D中的每一列，得出训练完毕的字典D。

进一步地，步骤五采用如下步骤具体实现：

5.1)设置各个训练样本的初始参数，并计算残差与训练完毕的字典D中每列d_l的内积，l为字典D的列数，找出最大内积所对应的下标p；

5.2)更新索引集，更新残差，并将迭代变量的数值增加1；

5.3)当更新后的残差<重构误差阈值时，结束迭代，且当i个训练样本t_i均都已找到稀疏表示α_i时，采用加权求和将重叠的图像块组合，得出无噪RAW图像，当i个训练样本t_i未全部找到稀疏表示α_i，则返回步骤5.1)，继续执行；

当更新后的残差≥重构误差阈值时，返回步骤5.2)，继续执行。

进一步地，步骤5.1)采用如下步骤具体实现：

找到训练样本集T中每个训练样本t_i的稀释表示，i＝1，2，…，i，针对其中的一个样本，设置每个训练样本t_i的初始化残差r⁽⁰⁾，设置重构误差阈值为n(Cσ)²，设置索引集Ω⁰＝0；

计算残差r^(L-1)与训练完毕的字典D中每列d_l的改进的内积，找出最大的内积所对应的下标p，p∈{1,2,…,k}，即：

式中，d_l为字典D的每一列，l＝1，2，…，k；r^(L-1)为OMP算法中迭代(L-1)次的残差；B_n为n×n的全1矩阵；γ用于调节颜色相关性；a为γ的分解因子；n表示图像块向量的大小，I为单位矩阵。

进一步地，步骤5.2)采用如下步骤具体实现：

更新索引集Ω^(L)(L)＝p，更新列d_p构成集D^(L)＝D(:,Ω^(L)(1:L)),D^(L)为迭代L次时的字典D；

通过最小二乘法计算得到稀疏逼近β_i^(L)＝argmin||t_i-D^(L)α_i||₂，

更新残差r^(L)＝t_i-D^(L)β_i^(L)，更新残差之后，迭代变量L的数值加1。

进一步地，步骤5.3)采用如下步骤具体实现：

5.3.1)当r^(L)＜n(Cσ)²，得到样本t_i的稀疏表示：α_i＝Sβ_i^(L)，结束迭代；并判断i个训练样本t_i是否都找到稀疏表示α_i；

当i个训练样本t_i都已找到稀疏表示α_i，则通过公式

采用加权求和将重叠的图像块组合，得出无噪RAW图像，式中，取λ＝30/σ；

当未全部找到i个样本t_i的稀疏表示，则返回步骤5.1)，继续执行；

5.3.2)当r^(L)≥n(Cσ)²时，则返回步骤5.2)，继续执行。

本发明的有益效果是：通过对OMP算法进行改进，使其适应于RAW图像的重建，通过该方法能够很好地估计出来噪声的参数，而且图像去噪计算复杂度减少至少1/3，本方法主要是针对在太空中低照度下，低信噪比成像系统去噪方法，不仅在太空中噪声较大的恶劣环境中具有较好的去噪性能，而且在成像过程中的RAW图像阶段去噪，能有效地在源头上抑制噪声，为后续图像处理和压缩提供必要的准备。

【附图说明】

图1为本发明一种基于稀疏表示的RAW图像去噪方法的流程图；

图2为本发明一种基于稀疏表示的RAW图像去噪方法中将图像块分离成G1RBG2颜色层步骤的示意图；

图3为本发明一种基于稀疏表示的RAW图像去噪方法将G1RBG2颜色层重组成G1RBG2矩阵步骤的示意图；

图4为本发明一种基于稀疏表示的RAW图像去噪方法中通过G1RBG2矩阵训练字典并重组图像的示意图。

图中，1.RAW矩形图像块，2.G1RBG2颜色层，3.G1RBG2矩阵，4.字典D，5.无噪RAW图像。

【具体实施方式】

下面结合附图和具体实施方式对本发明进行详细说明。

本发明公开了一种基于稀疏表示的RAW图像去噪方法，包括以下步骤：

步骤一、将待去噪RAW图像以给定步长分解成为固定大小相互重叠的RAW矩形图像块1；

步骤二、将步骤一中得到的RAW矩形图像块1按照不同颜色通道，重新排列成G1RBG2颜色层2；

步骤三、将步骤二中得到的G1RBG2颜色层2按照G1、R、B、G2的顺序依次拉伸后得到列向量形式的图像块向量，将每个图像块向量从左至右拼接成G1RBG2矩阵3，G1RBG2矩阵3的每一列为一个训练样本；

步骤四、从步骤三中的G1RBG2矩阵3中随机选择部分训练样本作为训练样本集T，采用K-SVD方法对训练样本集T进行学习得出字典D4。

步骤五、采用改进后的OMP算法对步骤四中得出的字典D4进行图像重建后，形成无噪RAW图像。

步骤一中的RAW矩形图像块1的大小为b*b，即RAW矩形图像块1中相邻两边的长度为均b个像素点，一般的b取值为8～16之间的偶数，优选的b取值为8；为了使图像块重叠，步长取小于b的偶数，优选的步长取值为2；且每个RAW图像块的左上角均是同一颜色分量。

步骤二中的G1RBG2颜色层2的颜色格式为：G1、R、B、G2；

步骤四中的采用K-SVD方法得出字典D4，其具体通过公式：

式中：n表示图像块向量的大小，n＝4*4*4＝64，k表示字典D4的大小，即列向量的个数，k＝256，C表示噪声增益，取值为1.15，α_ij为稀疏表示向量；R_ijX为图像块向量；σ为噪声的标准差，D为字典D4，i为字典D4中的列数，j为字典D4的列数。

步骤四中，采用K-SVD方法对训练样本集T进行学习得到字典D4，目的是找出更好的字典稀释来表示样本数据。

此步骤分为三个阶段进行实施：

阶段一：初始化阶段，具体步骤如下：

从训练样本集T中随机选区i个n维样本t_i，i＝1，2，…，i；用t_i来初始化初始字典D4的每一列d_l，l表示字典D4的列数，l＝1，2，…，l，即初始字典D4的每一列设置为i个训练样本的前l个；

初始字典D＝[d₁，d₂，…d_l]＝[t₁，t₂，…t_l]；

然后，归一化初始字典D4，即将字典D4进行迭代，

设置迭代次数变量J，初始状态下J＝1，每迭代一次后，次数变量J的值增加1，总共迭代J’次，J’为迭代总次数；

当迭代次数增加，去躁效果会提升，迭代次数越多，噪声标准差越小，而且根据零均值高斯噪声模型，采用经典PCA主成分分析方法估计标准差参数；

所以本方法设置：当初始标准差小于等于5时，J’＝10；当标准差大于5时，J’＝15；

阶段二：稀疏编码阶段，具体步骤如下：

当初始化阶段完成后，进入稀疏编码阶段，即固定字典D4，此步骤采用OMP算法来实现，目的是为了找到训练样本集T中每个训练样本的稀释表示。

4.1)依次针对训练样本集T的每一个样本：

初始化残差r⁽⁰⁾等于训练样本t_i，设置稀疏表示迭代次数L＝1，总共迭代L’＝3次，即稀疏度为3，设置索引集Ω⁰＝0；

其中，r⁽⁰⁾就是没有迭代的初始残差，r⁽¹⁾为第一次迭代计算生成的更新生成的残差，直到更新完r⁽³⁾，停止迭代；

4.2)计算前一次迭代中的残差r^(L-1)与字典D4中每列d_l的内积，通过r^(L-1)和d_l的内积值，找出最大的内积所对应的下标p，p∈{1,2,…,k}，即：p＝argmax_l＝1,...,k|<r^(L-1),d_l>|；

将索引集更新为Ω^(L)(L)＝p，更新所选择的列d_p构成集合D^(L)＝D(:,Ω^(L)(1:L))，D^(L)为迭代L次后的字典D4；

通过最小二乘法计算得到稀疏逼近α_i^(L)＝argmin||t_i-D^(L)α_i||₂，α_i为t_i的稀疏表示；

更新本次迭代中生成的残差r^(L)＝t_i-D^(L)α_i^(L)之后，再次进行迭代，将迭代变量的次数更新为L+1次；

4.3)将L+1与L’进行比较；

4.3.1)如果L+1>L’，则结束迭代，并判断i个训练样本t_i是否全部找到稀疏表示α_i；

当i个训练样本t_i全部找到稀疏表示α_i，则进入字典更新阶段；

当i个训练样本t_i未全部找到稀疏表示α_i，返回步骤4.2)；

4.3.2)如果L+1≤L’，继续返回步骤4.2)；

阶段三：字典更新阶段，将上述稀疏编码阶段得到的稀疏表示α_i，组成k行i列的稀疏矩阵A，其元素为α_i，通过稀疏矩阵A处理字典D4，依次更新字典D4中的每一个原子列d_l，l＝1，2，…，l，具体步骤如下：

4.4)找到使用了当前需更新原子列d_l的稀疏矩阵A的列元素集合{(l，i)|α_li≠0}，即对应于d_l原子列的是稀疏矩阵A中的l行，找到l行中不为零元素所在的列序号i构成的集合(l,i)∈w_l；

对于该集合w_l中所有的元素i，选择稀疏矩阵A和训练样本集T的第i列分别组合成新的稀疏矩阵A’和样本矩阵T’，其中i为小于g的非0整数，再去除新的稀疏矩阵和样本矩阵中的第l行后，得出稀疏矩阵A_l和样本矩阵T_l；

然后，将字典D4中去除第l行得到新的字典D4，通过公式E_l＝T_l-D_lA_l得出误差矩阵E_l；

然后，应用奇异值分解SVD分解误差矩阵E_l＝UΔV^T，得出U矩阵和V矩阵，更新列d_l为U矩阵的首列，更新系数值为V矩阵的首列，然后乘以Δ(1,1)；

4.5)判断字典D4的k列原子是否依次更新完；

4.5.1)如果字典D4的k列原子没有更新完，则返回步骤4.4)，继续执行；

4.5.2)如果字典D4的k列原子已更新完，判断迭代次数变量J的大小；

若J>J’，完成字典D4的更新，得出更新后的字典D4，并结束迭代；否则，返回步骤4.2)，继续执行；

步骤五中，OMP算法是基于l₀惩罚算法，实际就是寻找信号在字典D4的稀疏表示，即为下列优化问题：

该式的优化目标只有α_ij，因为字典D4在步骤四中已经训练完成。

通过优化完α_ij之后，再通过公式得出无噪图像。

在步骤五中，由于经典的OMP算法导致RAW图像空间的稀疏趋于R＝B＝G，使得彩色图像趋于灰度图像，无法保证重建图像块保持与原图像块具有相同的平均色彩，因此，采用对OMP算法进行改进，使其适应于RAW图像的重建。

通过改进后的OMP算法，计算所有n维图像块的样本集T_n×m的稀疏表示，固定步骤四中训练好的字典D4，其具体步骤如下：

5.1)找到每个训练样本t_i的稀释表示，i＝1，2，…，i，针对其中的一个样本，首先，初始化残差r⁽⁰⁾等于训练样本t_i，设置重构误差阈值为n(Cσ)²，设置索引集Ω⁰＝0；

然后，计算残差r^(L-1)与字典D4中每列d_l的改进的内积，找出最大的内积所对应的下标p，p∈{1,2,…,k}，即：

式中，d_l为字典D4的每一列，l为字典D4的列数，l＝1，2，…，l；r^(L-1)为OMP算法中迭代(L-1)次的残差，B_n为n×n的全1矩阵，γ用于调节颜色相关性；a为γ的分解因子，n表示图像块向量的大小，I为单位矩阵。。

即每个图像块和字典D4中的原子列都乘以尺度矩阵(I+(a/n)K)，然后在计算它们的内积。

值得注意的是，由于字典D4乘以尺度矩阵(I+(a/n)K)后，各个列向量需要重新归一化，即用(I+(a/n)K)DS处理字典D4，得到归一化之后的字典D4，式中，S代表对角矩阵，矩阵S的对角线元素对应为(I+(a/n)K)D各列2范数的转置：

式中，diag表示把向量元素对角化，sum表示对矩阵各列求和，“./”和“.*”分别表示矩阵点除和点乘；

5.2)更新索引集Ω^(L)(L)＝p，更新所选择的列d_p构成集D^(L)＝D(:,Ω^(L)(1:L))；

通过最小二乘法计算得到稀疏逼近β_i^(L)＝argmin||t_i-D^(L)α_i||₂；更新残差r^(L)＝t_i-D^(L)β_i^(L),更新残差之后，迭代变量L的数值加1；

5.3)如果r^(L)＜n(Cσ)²，得到样本t_i的稀疏表示：α_i＝Sβ_i^(L)，结束迭代；并判断i个训练样本t_i是否都找到稀疏表示α_i；

5.3.1)如果i个训练样本t_i都找到稀疏表示α_i，则通过公式

采用加权求和将重叠的图像块组合，估计去噪后的RAW图像；式中，λ表示拉格朗日乘子，它与噪声水平有关，噪声越大λ越小，因为噪声小的图像对重建图像影响大，噪声大的图像对重建图像影响应该小，取λ＝30/σ。

5.3.2)如果未全部找到i个样本t_i的稀疏表示，则返回步骤5.1)，继续执行；

5.4)若r^(L)≥n(Cσ)²，则返回步骤5.2)，继续执行。