CN111951183A - 基于近端交替惩罚算法的低秩全变分高光谱图像修复方法 - Google Patents

Abstract

一种基于近端交替惩罚算法的低秩全变分模型求解及其在高光谱图像修复领域的应用，包括以下步骤：步骤1)获取待去噪的高光谱图像数据

其中m、n、p分别表示高光谱图像的宽度、高度，及谱带的数量；步骤2)构建LRTV修复模型模型；步骤3)使用近端交替惩罚算法求解步骤4)输出修复后的高光谱图像。本发明可以有效的修复损失大量像素点的高光谱图像，并且其求解效率高。

Description

基于近端交替惩罚算法的低秩全变分高光谱图像修复方法

技术领域

本发明涉及一种基于近端交替惩罚算法的低秩全变分高光谱图像修复方法。

背景技术

高光谱遥感图像(Hyperspectral Imagery，HSI)由于其丰富的光谱信息，有着广泛的应用。特别是在精准农业、植被生态、地址矿产、大气环境检测等领域中有较大的应用价值。但是，高光谱图像在经由传感器采集过程中，不可避免地会损失部分像素点，这会严重降低图像质量，并影响后续处理。因此高光谱图像的修复是一个很重要的数据预处理过程。

全变分(Total Variation，TV)正则项是图像处理中广泛应用的修复模型，其能很好地保护图像边缘信息和保留平滑空间分段。近几年，基于TV的图像修复方法已被扩展至HSI修复领域。Kuiteing等人提出了用于HSI重构的迭代TV架构，Yuan等人提出了一个采用频谱空间自适应TV的HSI修复算法。

低秩矩阵分解(Low Rank Matrix Factorization，LRMF)也已被广泛用作图像分析，网页搜索和计算机视觉等领域。低秩模型描述了在高维数据中发现和利用低维结构的问题。对于HSI，光谱上的相邻波段通常表现出很强的相关性，并且HSI中的相邻像素通常也具有高度空间相关性，这两者都揭示了高光谱图像的低秩结构。基于此，LRMF已被广泛应用于HSI修复领域。

He等人提出了低秩全变分(Low Rank and Total Varaition，LRTV)模型，提出将TV正则项与LRMF结合起来，利用低秩模型获取HSI数据的频谱相关性，使用TV正则化保留HSI数据平滑空间分段。LRTV模型在高光谱图像修复领域获得了很好的效果。

求解LRTV模型，作为一个凸优化问题，可以用许多经典的凸优化算法对其求解，如交替方向乘子法(Alternating Direction Method of Multipliers，ADMM),原始对偶分裂法(Primal-dual Splitting Algorithm，PDS)等。本发明使用了一个新的算法，即近端交替惩罚算法(Proximal Alternating Penalty Algorithm，PAPA)来对LRTV模型进行求解并将其应用于高光谱图像修复。PAPA算法主要使用经典二次惩罚法，交替最小化方法，内斯特罗夫的加速方法和参数自适应策略方法的新颖组合，来有效求解不光滑的约束凸优化问题，并且在非遍历意义上获得了著名的

收敛，k是迭代次数。二次惩罚法是处理约束问题的经典优化框架，此方法单独使用时通常效率不高，但当其与交替策略结合起来，却可以获得很好的效果。应用交替最小化方法解耦联合变量z＝(x,y)，在x和y之间进行交替求解，此方法已被广泛应用于非约束问题的求解。PAPA算法的关键之处在于内斯特罗夫的加速方案中的参数自适应策略，这个方法可以加快PAPA的收敛速度，并且可以自动更新惩罚项参数和其他参数，无需手动调整。因此，在求解效率方面，PAPA算法显著高于ADMM，PDS等算法。

发明内容

本发明针对ADMM，PDS等算法求解LRTV模型计算效率低下，修复效果差等问题，提出了使用一个新的算法PAPA来有效求解LRTV模型，并将其应用到高光谱图像修复领域。实验结果表明，其修复效果与计算效率都得到了显著提升。

本发明解决其技术问题所采用的技术方案是：

基于近端交替惩罚算法的低秩全变分模型求解并将其应用于高光谱图像修复领域，包括以下步骤：

步骤1，获取待修复的高光谱图像数据

其中m、n、p分别表示高光谱图像的宽度、高度、及谱带数量；

步骤2，构建LRTV修复模型；

步骤3，使用近端交替惩罚算法求解；

步骤4，输出修复后的高光谱图像。

本发明的效果表现在：高光谱图像的修复具有一定的挑战性，通过LRTV模型可以有效的对高光谱图像进行修复，但用一般的基于原始-对偶分裂的算法如ADMM，PDS等，求解效率过低，修复效果较差，因此本发明提出了用近端交替惩罚算法来求解LRTV模型，并将其应用于高光谱图像的修复。实验结果表明，本发明提出的方法可以有效的修复损失大量像素点的高光谱图像，并且其求解效率明显高于ADMM，PDS，ALM等算法。

附图说明

图1是损失85％像素点的高光谱图像。

图2是本发明的流程图。

图3是使用ADMM算法求解LRTV模型修复后的高光谱图像及其所消耗时间。

图4是使用PDS算法求解LRTV模型修复后的高光谱图像及其所消耗时间。

图5是使用PAPA算法求解LRTV模型修复后的高光谱图像及其所消耗时间。

具体实施方式

下面结合附图，进一步说明本发明的技术方案。

基于近端交替惩罚算法的低秩全变分高光谱图像修复方法，包括以下步骤：

步骤1)输入丢失像素点的高光谱图像数据

其中m、n、p分别表示高光谱图像的宽度、高度、及谱带数量；具体如图1所示；

步骤2)针对高光谱图像修复的LRTV模型定义如下：

是输入的丢失像素点的图像，

表示输出的HSI数据矩阵。||Y||_LR表示低秩项。根据在信号和图像处理中对于使用低秩正则项的研究，本发明使用核范数来表示低秩项，其定义为：||Y||_LR≈||Y||_*。||Y||_tv表示全变分正则项，其定义为：

||·||_F为Frobenius范数。其定义为：

是一个索引矩阵，0表示输入数据的丢失元素，1表示输入数据T的可观测元素。支持集Ω定义为{(i,j)|q_i,j＝1}。α与β则分别表示低秩项与TV正则项的参数。而δ的值取决于T中观察到的元素丢失率。

步骤3)使用PAPA算法求解LRTV模型

31)优化低秩全变分模型：

其中||X||₁表示l₁范数，其定义为：

D是矩阵的差分算子。

32)令f(X)＝||X||₁，g(Y)＝||Y||_*，h(Y)＝||p_Ω(Y-T)||² _F；

33)初始化X⁰，

Y⁰，

α，β，ρ₀；

34)开始迭代，其中k表示迭代次数，初始值为0；

341)更新X^k+1：

342)更新Y^k+1：

343)更新

344)更新ρ_k+1：

ρ_k+1:＝(k+2)ρ₀ (6)

其中

这一步使用了经典的二次惩罚方法，

是一个二次惩罚函数，Ax+By+c＝0是目标函数的约束项，projκ(·)表示凸集上的投影算子，ρ是惩罚参数，

prox表示函数的近似算子。其定义为:

表示对函数

中的变量y求梯度，其值为

其中

在开始迭代前，首先给X与Y赋初值，同时使

低秩项参数α与全变分参数β分别设置为0.009与1.2。ρ₀的初值设为0.5。

根据近似算子的定义，可以将式(3)表示为：

用软阈值求解式(7)得到：

对于式(4)，可以通过引入以下引理1来求解：

引理1：令

是一个给定矩阵，则对秩为r的矩阵W的奇异值分解定义如下：

W＝UE_rV^T,E_r＝diag({σ_i}_1≤i≤r) (9)

diag({σ_i}_1≤i≤r)表示对角矩阵，其i行对应的对角元素为σ_i，σ_i是矩阵W的第i个奇异值。U为左奇异矩阵，V^T为右奇异矩阵。

然后，奇异值收缩运算符将服从

其中

J_ξ(W)＝UJ_ξ(W)V^T，J_ξ(E_r)＝diag{max((σ_i-ξ),0)} (11)

ξ为输入的阈值，通过引理1，可以得到式7的解：

算法1：

1.输入：含有损失像素点的高光谱图像T

2.初始化：α，β，ρ₀，

3.根据(3)计算X^k+1；

4.根据(4)计算Y^k+1；

5.根据(5)更新

6.根据(6)更新ρ_k+1；

7.收敛之后输出Y^k+1

步骤4)输出修复后的高光谱图像，如图3所示。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (2)

1.基于近端交替惩罚算法的低秩全变分高光谱图像修复方法，包括以下步骤：

步骤1)输入丢失像素点的高光谱图像

其中m、n、p分别表示高光谱图像的宽度，高度，及谱带的数量；

步骤2)针对高光谱图像修复的低秩全变分模型定义如下：

其中

是输入的丢失像素点的图像，

表示输出的HSI数据矩阵；||Y||_LR表示低秩项；使用核范数来表示低秩项，其定义为：||Y||_LR≈||Y||_*，||Y||_tv表示全变分正则项，其定义为：

||·||_F为Frobenius范数，其定义为：

X_i,j为X的元素值，

θ∈{0,1}^mn×p是一个索引矩阵，0表示输入数据的丢失元素，1表示输入数据T的可观测元素；支持集Ω定义为{(i,j)|q_i,j＝1}，表示未丢失的像素集；α与β分别表示低秩项与TV正则项的参数；而δ的值取决于T中观察到的元素丢失率；

步骤3)使用近端交替惩罚算法求解低秩全变分模型；

步骤4)输出修复后的高光谱图像。

2.如权利要求1所述的一种基于近端交替惩罚算法的低秩全变分高光谱图像修复方法，其特征在于：步骤3)具体包括：

31)优化低秩全变分模型：

其中||X||₁表示l₁范数，其定义为：

||Y||_*表示Y的核范数，D是矩阵的差分算子；

32)令f(X)＝||X||₁，g(Y)＝||Y||_*，h(Y)＝||p_Ω(Y-T)||² _F；

33)初始化X⁰，

Y⁰，

α，β，ρ₀；

34)开始迭代，其中k表示迭代次数，初始值为0；

341)更新X^k+1:

其中

ρ_k是惩罚参数；

AX+BY+C＝0是目标函数的约束项，proj_κ(·)表示凸集上的投影算子；近似算子的定义如下：

根据近似算子的定义，可将

改写为：

对于l₁范数的近似算子，可用软阈值求解，得到：

342)更新Y^k+1：

其中

表示对函数

中的变量Y求梯度，其值为

其中

对于式(7)，通过引入引理1来求解：

引理1：令

是一个给定矩阵，则对秩为r的矩阵W的奇异值分解定义如下：

W＝UE_rV^T,E_r＝diag({σ_i}_1≤i≤r) (8)

diag({σ_i}_1≤i≤r)表示对角矩阵，其i行对应的对角元素为σ_i，σ_i表示矩阵W的第i个奇异值，U为矩阵W的左奇异矩阵，V^T为矩阵W的右奇异矩阵；

然后，奇异值收缩运算符将服从

其中

J_ξ(W)＝UJ_ξ(W)V^T，J_ξ(E_r)＝diag{max((σ_i-ξ),0)} (10)

ξ为输入的阈值，通过引理1，可以得到式7的解：

343)更新

344)更新ρ_k+1：

ρ_k+1:＝(k+2)ρ₀ (13)。

Images