基于稀疏与低秩矩阵逼近的高光谱图像恢复方法与流程

本发明涉及图像处理，尤其是涉及基于稀疏与低秩矩阵逼近的高光谱图像恢复方法。
背景技术：
：高光谱成像技术是自20世纪发展起来的一种对地观测技术，在军事、地质勘测、农业检测等方面都实现了很高的实用价值。高光谱图像数据是一个三维的图像数据，除了二维图像本身外，还有波段维数，其主要特点是将传统的图像空间与光谱信息融合为一体。高光谱图像数据包含的丰富的地物光谱信息，可实现地物的目标识别等。虽然随着计算机技术，光谱成像技术的发展，高光谱图像数据在成像质量上有了很大的提高。但是不可避免的，由于成像仪物理缺陷、大气污染、传输损失和校准问题导致高光谱图像受到混合噪声影响。这些混合噪声包括高斯噪声、冲击噪声、条带噪声和死行。因此高光谱图像的去噪问题亟待解决，是高光谱图像数据在后续的信息分析中重要的预处理工作。目前，有许多图像去噪技术被应用于高光谱图像的去噪工作中。受到灰度图像去噪方法的启发，非局部均值方法、SVD分解法、Dabov等人在“Imagedenoisingbyspare3-Dtransform-domaincollaborativefiltering.IEEETrans.onImageProcess.”中提出的BM3D方法以及Yuan等人在“Hyperspectralimagedenoisingemployingaspectral-spatialadaptivetotalvariationmodel.IEEEGeosci.RemoteSens.”中提出SSATV方法，这些方法都可以逐波段地处理高光谱图像。然而这些方法的效果都不尽如人意，因为它们都忽略了高光谱图像不同波段之间的关联。为了利用不同波段图像之间的关联，一些研究者提出将大的三维数据分割成小的三维数据，再将其转换为二维数据后，利用其无噪声矩阵的低秩性来进行低秩矩阵恢复，从而达到去噪效果。例如Zhang等人在“HyperspectralImageRestorationUsingLow-RankMatrixRecovery.IEEETrans.onGeoscienceandRemoteSensing”中提出的LRMR的方法。技术实现要素：本发明的目的在于提供不仅去除绝大部分噪声，还保留丰富的图像细节的一种基于稀疏与低秩矩阵逼近的高光谱图像恢复方法。本发明包括以下步骤：(1)设一组多波段高光谱图像数据d，d的大小为M×N×B，其中M和N分别表示每个波段的高光谱图像的长和宽，B表示共有多少个波段，d的估计噪声等级为η；(2)初始化变量，去噪后数据有噪声数据高斯噪声等级η(0)＝η；(3)初始化迭代，令循环变量k＝1，设置最大外层循环次数K；(4)迭代正则化(5)以中心坐标位置为(i,j)，其中i和j分别表示中心点的横纵坐标，提取大小为w×h×B的数据其中w和h表示数据块的宽和高，B是高光谱数据波段数；(6)将三维数据转换成大小为wh×B的二维矩阵D；(7)对D构建加权的Schatten-p范式低秩矩阵近似模型；(8)使用扩展的拉格朗日乘数法对步骤(7)构建的加权的Schatten-p范式低秩矩阵近似模型进行求解得到去噪后的矩阵A；(9)更新噪声等级η(k)；(10)按照一定的步长取所有中心坐标(i,j)，重复步骤(5)～(9)，分别得到其对应的A(i,j)；(11)将所有大小为wh×B的A(i,j)转换回大小为w×h×B的三维数据(12)将所有拼合成完整的大小为M×N×B的三维数据(13)将k＞K作为循环终止条件，若k不满足大于K的条件，则在将k的值增加1后返回步骤(4)，否则直接执行步骤(14)；(14)将作为最后去除所有噪声后的高光谱图像数据。在步骤(5)中，所述提取大小为w×h×B的数据的方法是对于B个波段的M×N高光谱图像，分别取出其中以(i,j)为中心，w×h为大小的图像块后层叠在一起得到大小为w×h×B的三维数据块。在步骤(6)中，所述将三维数据转换成大小为wh×B的二维矩阵D是将B个w×h的图像块分别按顺序排列成wh×1的像素列后，再拼成一个wh×B的矩形D。在步骤(7)中，所述加权的Schatten-p范式低秩矩阵近似模型如下：s.tD＝A+E+N,||N||F≤η其中，C和λ是系数，A表示未受到噪声影响的干净图像数据，E表示冲击噪声、死行和条带噪声的混合，N表示高斯噪声，其噪声等级η人为设置；||·||1表示矩阵的1范数，||·||F表示矩阵的Frobenius范数；表示加权的Schatten-p范式，具体形式如下：其中，0＜p≤1，r＝min{wh,B}，σi是A的第i个奇异值，其中σi(D)表示D的第i个奇异值。在步骤(8)中，所述使用扩展的拉格朗日乘数法对步骤(7)构建的加权的Schatten-p范式低秩矩阵近似模型进行求解得到去噪后的矩阵A的具体步骤如下：(8.1)初始化变量：Γk＝D，Ak＝Ek＝0，βk＞0，ρ＞1，k＝0，Γ为拉格朗日乘子，βk是第k次迭代时约束项的惩罚系数，ρ是惩罚系数迭代因子，k为当前迭代数；(8.2)更新其中(8.3)更新(8.4)更新(8.5)更新Γk+1＝Γk-βk(Ak+1+Ek+1+Nk+1-D)；(8.6)令βk+1＝ρ*βk；(8.7)令k＝k+1；(8.8)若||Γk-Γk-1||1＜convergence，其中convergence为收敛参数，则执行步骤(8.9)，否则执行步骤(8.2)；(8.9)将Ak+1作为去噪后的矩阵A。在步骤(8.3)中，所述更新的方法如下：令矩阵Ek+1的每个元素：(Ek+1)ij＝sign(Xij)·max(|Xij|-λ,0)其中,(·)ij表示矩阵第i行第j列的元素，sign(·)为符号函数。在步骤(8.4)中，更新的具体步骤如下：(8.4.1)令对Y进行SVD分解得到Y＝UΣVΤ，其中Σ＝diag(σ1,...,σr)，也即σi为Y的奇异值；(8.4.2)初始化迭代，令循环变量i＝1，设置最大外层循环次数r；(8.4.3)利用σi,ωi,p，再通过广义软阈值法得到Ak+1的估计矩阵的第i个奇异值δi；(8.4.4)若i≥r，则直接执行步骤(8.4.5)，否则令i＝i+1，执行步骤(8.4.3)；(8.4.5)令变量Δ＝diag(δ1,...δr)，其中diag(·)为矩阵对角排列函数；(8.4.6)令Ak+1＝UΔVΤ；在步骤(8.4.3)中，所述利用σi,ωi,p，再通过广义软阈值法得到Ak+1的估计矩阵的第i个奇异值δi的具体步骤如下：(8.4.3.1)令变量σ＝σi，变量ω＝ωi，令这是一个关于ω和p的函数；(8.4.3.2)若|σ|≤τp(ω)，则令函数Sp(σ；ω)＝0，执行步骤(8.4.3.8)；否则执行步骤(8.4.3.3)；(8.4.3.3)初始化阈值法变量为δ(0)＝|σ|；(8.4.3.4)初始化迭代，令循环变量t＝0，设置最大外层循环次数J；(8.4.3.5)令δ(t+1)＝|σ|-ωp(δ(t))p-1；(8.4.3.6)若t≥J，则直接执行步骤(8.4.3.7)，否则令t＝t+1，执行步骤(8.4.3.5)；(8.4.3.7)设定函数Sp(σ；ω)＝sgn(σ)δ(t)，其中sgn为sign符号函数；(8.4.3.8)令δi＝Sp(σ；ω)。在步骤(9)中，所述更新噪声等级η(k)的方法为：其中，γ是比例因子。在步骤(12)中，所述将所有拼合成完整的大小为M×N×B的三维数据的具体步骤为:将放置回中以(i,j)为中心，w×h×B为大小的位置，各个之间重叠的部分求其平均值即可。本发明采用一种新颖的方法来去除高光谱图像在获取和传输过程中产生的混合噪声，这些噪声是混合的高斯噪声、冲击噪声、条带噪声和死行，针对高光谱图像的性质，本发明实现了一种稀疏与低秩矩阵逼近的高光谱图像恢复方法，不仅去除了绝大部分噪声，还保留了丰富的图像细节。本发明提出了一个基于低秩矩阵分解的迭代式的高光谱图像恢复方法，利用处理后的数据矩阵的低秩性，用加权的schattenp范式对其构建低秩矩阵近似模型，并给出了其求解方法。独特的加权的schatten-p范式的低秩正则项对低秩矩阵近似模型效果有明显的提升，其次在求解过程中自动地估计和更新对噪声等级的估计，达到了对不同的噪声能够自适应地处理。通过这些新颖的技术手段，不仅能很好地去除混合噪声，还极大程度地保留了有利用价值的图像信息，本发明的效果已领先于当前流行的高光谱图像恢复的方法。本发明具有如下技术效果：1.本发明先利用高光谱图像的经过处理后的低秩性构建扩展的鲁棒主成成分分析模型，这个模型很好地保留了低秩性和稀疏性。对于模型的低秩正则化项，本发明使用一种新颖的加权的Schatten-p范式，该范式不仅通对秩成分通过事假不同的权来控制，从而适应不同程度噪声，而且因为其对低秩假设有更好的实现，从而在理论和实际都更好地达到了去噪效果。2.针对上述难解的非凸优化去噪模型，本发明采用扩展的拉格朗日乘数法结合广义软阈值法对其进行求解。不仅如此，本发明还使用了迭代正则化框架，通过估计和调整噪声等级，使整个模型能自适应地应对不同的高光谱数据的恢复问题。附图说明图1是使用本发明和现有三种方法对烟囱的真实视频帧进行大气扰动去除后得到的图像；图2是使用本发明和现有三中方法对城市的模拟扰动视频进行大气扰动取出后得到的图像。具体实施方式本发明的具体实施步骤包括：步骤1，获取受到混合噪声影响的高光谱图像数据。利用高光谱成像仪获得一组多波段高光谱图像数据d，将它归一化到[0,1]。它的大小为M×N×B，其中M和N分别表示每个波段的高光谱图像的长和宽，B表示共有多少个波段。这组数据的估计的噪声等级为η＝20/255。步骤2，初始化迭代变量。(2a)令去噪后数据有噪声数据高斯噪声等级η(0)＝η。(2b)初始化迭代，令循环变量k＝1，中心横纵坐标分别为i＝10,j＝10。步骤3，迭代正则化步骤4，获取待建立低秩模型的二维数据。(4a)对于B个波段的M×N高光谱图像，分别取出其中以(i,j)为中心，20×20为大小的图像块后层叠在一起得到大小为20×20×B的三维数据块(4b)将B个20×20的图像块分别按顺序串联成400×1的像素列后，再拼成一个400×B的矩形D。步骤5，构建加权的Schatten-p范式低秩矩阵近似模型如下：s.tD＝A+E+N,||N||F≤η其中，C＝0.007，λ＝1.2，p＝0.7。p＝0.7，r＝min{400,B}，σi是A的第i个奇异值，其中σi(D)表示D的第i个奇异值。步骤6，使用扩展的拉格朗日乘数法对上述模型进行求解。(6a)初始化变量：Γh＝D，Ah＝Eh＝0，βh＞0，ρ＞1，h＝0，Γ为拉格朗日乘子，βh是第h次迭代时约束项的惩罚系数，ρ是惩罚系数迭代因子，h为当前迭代数。(6b)更新其中(6c)更新Eh+1，令矩阵Eh+1的每个元素(Eh+1)ij＝sign(Xij)·max(|Xij|-λ,0)，其中(6d)更新Ah+1如下：(6d.1)令对Y进行SVD分解得到Y＝UΣVΤ，Σ＝diag(σ1,...,σr)。(6d.2)初始化迭代，令循环变量q＝1，设置最大外层循环次数r,r如步骤5中所述。(6d.3)利用步骤5中的定义的ωq和p，步骤(6d.1)中的σq，通过如下广义软阈值法得到Ah+1的估计矩阵的第q个奇异值δq：(6d.3a)令变量σ＝σq，变量ω＝ωq，令函数(6d.3b)若|σ|≤τp(ω)，则令变量Sp(σ；ω)＝0，执行步骤(8.8)；否则执行步骤(8.3)。(6d.3c)令初始化阈值法变量δ(0)＝|σ|。(6d.3d)初始化迭代，令循环变量t＝0，设置最大外层循环次数J＝50。(6d.3e)令δ(t+1)＝|σ|-ωp(δ(t))p-1。(6d.3f)若t≥J，则直接执行步骤(6d.3g)，否则令t＝t+1，执行步骤(6d.3e)。(6d.3g)设定函数Sp(σ；ω)＝sgn(σ)δ(t)。(6d.3h)令δq＝Sp(σ；ω)。(6d.4)若q≥r，则直接执行步骤(6d.5)，否则令q＝q+1，执行步骤(6d.3)。(6d.5)令矩阵变量Δ＝diag(δ1,...δr)，其中diag(·)为矩阵对角排列函数。(6d.6)令Ah+1＝UΔVΤ。(6e)更新Γh+1＝Γh-βh(Ah+1+Eh+1+Nh+1-D)。(6f)令βh+1＝ρ*βh。(6g)令h＝h+1。(6h)若||Γh-Γh-1||1＜convergence，其中收敛参数convergence＝10-7，则执行步骤(5.9)，否则执行步骤(5.2)。(6i)将Ah+1作为去噪后的矩阵A(i,j)。步骤7，更新噪声等级步骤8，横纵坐标都按照步长为4，取遍所有中心坐标(i,j)，重复步骤(4)～(7)，分别得到其对应的A(i,j)。步骤9，得到本轮迭代的去噪高光谱数据。(9a)将所有大小为400×B的A(i,j)转换回大小为20×20×B的三维数据(9b)将所有放置回中以对应的(i,j)为中心，w×h×B为大小的位置，各个之间重叠的部分求其平均值得到大小为M×N×B的步骤10，将k＞6作为循环终止条件，若不满足，在将k的值增加1后返回步骤3，否则执行步骤11。步骤11，将作为最后去除所有噪声后的高光谱图像数据。本发明通过以下实验进行优点和有效性的证明。1.实验条件：实验室用台式机参数：CPU为Inter(R)Core(TM)i7-2600，主频为3.40GHz，内存为4G，操作系统为Win764位系统，实验平台为Matlab2014b。2.实验结果以及结果分析：实验一，用本发明和现有的方法进行真实高光谱图像数据恢复。从网络上下载一组真实的受混合噪声影响高光谱图像数据Urban，其大小为307×307×188，将其每个像素减去最小像素后，再除以最大与最小像素值之差，便得到归一化到[0,1]后的数据。再使用本发明提出的方法进行去噪处理，并用
背景技术：
中提到的BM3D、SSATV和LRMR三种方法做对比实验。实验结果如图2所示，其中图2(a)为其中一张原始的高光谱图像数据，图2(b)为BM3D方法效果，图2(c)为SSATV方法效果，图2(d)为LRMR方法效果，图2(e)为本发明效果。图2的实验结果说明：从直观视觉上看，图2(b)、图2(c)虽然去掉了大部分噪声，但同时也损失掉了许多细节信息，图2(b)还有明显的去噪痕迹。图2(d)的结果还有很多未去除的噪声。图2(e)中本发明的去噪效果明显，不仅很好地去除了噪声，还在很大程度上保留了丰富的细节信息。实验二，用本发明和现有的方法进行模拟高光谱图像数据恢复。从网络上下载一组未被噪声影响的高光谱图像，其大小为256×256×191。人为地对所有帧随机加上高斯噪声、冲击噪声、条带噪声和死行。再使用本发明提出的方法进行去噪处理，并用BM3D、SSATV和LRMR三种方法做对比实验。实验结果如图3所示，其中图3(a)为其中一张清晰的高光谱图像数据，图3(b)为加上模拟噪声后的高光谱图像，图3(c)为BM3D方法效果，图3(d)为SSATV方法效果，图3(e)为LRMR方法效果，图3(f)为本发明效果。表1给出了本发明与上述三种方法的平均峰值信噪比(MPSNR)和平均结构相似性(MSSIM)的对比结果。图3和表1的实验结果说明：图3(b)中受到混合噪声影响，视觉效果很差，信息干扰多；图3(c)损失大量细节信息，对死行噪声基本没有效果；图3(d)中所有噪声的去除都不尽如人意，模糊且有死行；图3(e)的结果相比前面的图都要好，但是在图中被绿色虚线圈出的部分，仍然有黑色印记残留，这是处理死行噪声不完整造成的。图3(f)中，本发明的方法不仅很好地去除了包括死行在内的各种噪声，也在在很大程度上保留了丰富的细节信息，与原数据图像相比，最为接近，效果最好。表1BM3DSSATVLRMR本发明MPSNR24.95817.34427.47430.303MSSIM0.61770.27520.81140.8711由表1表明，本发明的MPSNR和MSSIM值都高于现有的三种方法。当前第1页1 2 3