Hyperspectral Image Super-Resolution via Nonlocal Low-Rank Tensor Approximation and Total Variation Regularization

Abstract

Hyperspectral image (HSI) possesses three intrinsic characteristics: the global correlation across spectral domain, the nonlocal self-similarity across spatial domain, and the local smooth structure across both spatial and spectral domains. This paper proposes a novel tensor based approach to handle the problem of HSI spatial super-resolution by modeling such three underlying characteristics. Specifically, a noncovex tensor penalty is used to exploit the former two intrinsic characteristics hidden in several 4D tensors formed by nonlocal similar patches within the 3D HSI. In addition, the local smoothness in both spatial and spectral modes of the HSI cube is characterized by a 3D total variation (TV) term. Then, we develop an effective algorithm for solving the resulting optimization by using the local linear approximation (LLA) strategy and the alternative direction method of multipliers (ADMM). A series of experiments are carried out to illustrate the superiority of the proposed approach over some state-of-the-art approaches.

1. Introduction

Hyperspectral images (HSIs) are recordings of reflectance of light of some real world scenes or objects including hundreds of spectral bands ranging from ultraviolet to infrared wavelength [1,2]. The abundant spectral bands of HSIs could possess fine spectral feature differences between various materials of interest and enable many computer vision and image processing tasks to be more successfully achievable. However, due to various factors, e.g., the constraints of imaging sensor and time constraints, the acquired HSIs unfortunately are of low spatial resolution, which cannot provide any further help for high precision processing requirements in many fields including mineralogy, manufacturing, environmental monitoring, and surveillance. Therefore, the task of enhancing the spatial resolution of HSIs is a valuable research issue and has received much attention in recent years.

Super-resolution is a widely-used signal post-processing technique for image spatial resolution enhancement, which produces the high resolution (HR) image from the observed low resolution (LR) image or sequence that are noisy, blurred and downsampled [3,4]. Tsai and Huang [5] presented the first algorithm for the reconstruction of the HR image from a set of undersampled, aliased and noise-free LR images. Then, Kim et al. [6] extended the work to consider observation noise as well as the spatial blurring effect. Along this line, several frequency domain methods have been proposed. See, e.g., [7,8]. However, such methods cannot utilize the prior information in the spatial domain. Therefore, to overcome their weaknesses, many kinds of spatial domain approaches have been developed. For instance, Ur et al. [9] proposed a non-uniform interpolation by exploiting the structural information in the spatial domain; Irani et al. [10] proposed an iterative back projection method, which can be considered as solving a simple least square problem.

In the last few years, incorporating the prior information of HR auxiliary images into the process of HSI super-resolution. In particular, multispectral and hyperspectral data fusion have been becoming more and more popular. See, e.g., [11,12,13,14] and references therein. Due to the limitations of remote sensing systems, however, it is not easy to get such HR images. Therefore, single HSI super-resolution has still attracted much interest in many real world scenarios. In [2], Akgun et al. proposed a novel hyperspectral image acquisition model, with which a projection onto convex sets based super-resolution method was proposed to enhance the resolution of HSIs. Guo et al. [15] used the unmixing information and total variation (TV) minimization to produce a higher resolution HSI. By modeling the sparse prior underlying HSIs, a sparse HSI super-resolution model was proposed in [16]. Zhang et al. [17] proposed a maximum a posteriori based HSI super-resolution reconstruction algorithm, in which PCA is employed to reduce computational load and simultaneously remove noise. Huang et al. [18] presented a novel super-resolution approach of HSIs by joint low-rank and group-sparse modeling. Their approach can also deal with the situation that the system blurring is unknown. In a recent study, Li et al. [19] explored sparse properties in both spectral and spatial domains for HSI super-resolution. Recently, sparse representation and compressed sensing have successfully been used in various real image super-resolution applications. See, e.g., [20,21]. Actually, such methods can be naturally used to enhance the spatial resolution of one single HSI cube in a band-by-band manner.

In this paper, following the ideas of our preliminary work [22] and the work [23] for MRI super-resolution, we consider the HSI cube as the tensor with three modes, namely, width, height, and band, and then exploit the underlying structures in both spatial and spectral domain by using direct tensor modeling techniques to achieve the spatial resolution enhancement. Precisely, the spectral bands of an HSI commonly have global strong correlations, and for each local fullband patch of an HSI (which is stacked by patches at the same location of HSI over all bands), there are many same size fullband patches similar to it; this spatial-and-spectral correlation is modeled by a nonconvex low-rank tensor penalty. In addition, for each voxel, its intensity seems to be almost equal to those in its neighbourhood from both the spatial and the spectral viewpoints; this local spatial-and-spectral smoothness is exploited by using the 3D TV penalty. As such, the HSI super-resolution task resorts to solving a nonconvex optimization problem, which could be effectively solved by combing the local linear approximation (LLA) strategy and the alternative direction method of multipliers (ADMM) that we shall show later.

This paper is organized as follows. Section 2 introduces some notations and preliminaries of tensors, which will be used for presenting our procedure. In Section 3, the proposed tensor model and its motivations are introduced. Then, an effective algorithm is developed for solving the resulting nonconvex optimization problem in Section 4. Extensive experiments on three real datasets are carried out in Section 5 to illustrate the outperformance of our model over several state-of-the-art ones. Finally, we conclude this paper with some discussions on future research in Section 6.

2. Notation and Preliminaries

It is known that a tensor can be regarded as a multi-index numerical array, and its order is defined as the number of its modes or dimensions. A real-valued tensor of order K is denoted by $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \dots \times I_K}$ and its entries by $x_{i_1,i_2,\dots ,i_K}$. We then can consider an $n\times 1$ vector x as a tensor of order one, and an $n\times m$ matrix X as a tensor of order two. Following [24], we shall provide a brief introduction on tensor algebra.

The inner product of two same-sized tensors $\mathcal{X},\mathcal{Y}$ is defined as

$$\langle \mathcal{X},\mathcal{Y}\rangle :=\sum _{i_1,i_2,\dots ,i_K}{x}_{i_1,i_2,\dots ,i_K}·y_{i_1,i_2,\dots ,i_K}.$$

(1)

Then, the corresponding Frobenius norm can be defined as ${\parallel \mathcal{X}\parallel }_F=\sqrt{\langle \mathcal{X},\mathcal{X}\rangle }$. The so-called mode-n matricization or unfolding of the tensor $\mathcal{X}$ is denoted as ${\mathcal{X}}_{\left(n\right)}$, where the tensor element $(i_1,i_2,\dots ,i_K)$ is mapped to the matrix element $(i_n,j)$ satisfying $j=1+{\sum }_{k=1,k\ne n}^K(i_k-1)J_k$ with $J_k={\prod }_{m=1,m\ne n}^{k-1}{I}_m$. Then, the operation “fold" is defined as ${\mathrm{fold}}_n\left({\mathcal{X}}_{\left(n\right)}\right):=\mathcal{X}$. The mode-n multiplication of $\mathcal{X}$ with the matrix U is denoted by $\mathcal{X}{\times }_{n}U$ and, elementwise, we have ${\left(\mathcal{X}{\times }_{n}U\right)}_{i_1,\dots ,i_{n-1}j{i}_{n+1}\dots ,i_K}={\sum }_{i_n}{x}_{i_1,i_2,\dots ,i_N}·u_{j,i_n}.$ The multi-linear rank is defined as an array $(r_1,r_2,\dots ,r_K)$, where $r_n=\mathrm{rank}\left({\mathcal{X}}_{\left(n\right)}\right)$, $n=1,2,\dots ,K$.

Following [25], for a given matrix X, its folded-concave norm is defined as ${\parallel X\parallel }_{P_{\lambda }}:={\sum }_{j=1}^r{P}_{\lambda }\left({\sigma }_j\left(X\right)\right)$, where ${\sigma }_j\left(X\right)$ is the j-th singular value of X and r is its rank, and $P_{\lambda }$ is the so-called folded-concave penalty function defined by [26]. In this work, we use one popular folded-concave penalty i.e., the minmax concave plus (MCP) penalty, which has the form:

$$P_{\lambda }\left(t\right)=\left\{\begin{array}{cc}a{\lambda }^2/2,\hfill & \mathrm{if}|t|\ge a\lambda ,\hfill \\ \lambda |t|-\frac{t^2}{2a},\hfill & \mathrm{otherwise}..\hfill \end{array}\right.$$

(2)

Here, a is a fixed constant. Then, similar to the tensor nuclear norm defined in [27], we can define the tensor MCP penalty by applying the MCP penalty function to each matricization ${\mathcal{X}}_{\left(i\right)}$:

$${\parallel \mathcal{X}\parallel }_{P_{\lambda }}=\sum _i^N{\alpha }_i{\parallel {\mathcal{X}}_{\left(i\right)}\parallel }_{P_{\lambda }},$$

(3)

where ${\alpha }_i\ge 0$ and satisfies ${\sum }_i^N{\alpha }_i=1$.

3. HSI Super-Resolution via Nonlocal Low-Rank Tensor Approximation and TV Regularization

In this section, we firstly propose the observation model. Then, we use the 3D TV regularization to describe local smoothness of one HSI, and adopt a path-based tensor folded-concave penalty to characterize both the global correlation and the nonlocal self-similarity of the HSI. Finally, we propose a novel regularization model for dealing with the HSI super-resolution task. The detailed framework of our procedure is presented in Figure 1.

3.1. Observation Model

The low spatial resolution HSI is expressed by the following observation model:

$$\mathcal{Y}=DS\mathcal{X}+\mathcal{E},$$

(4)

where the tensor $\mathcal{Y}\in {\mathbb{R}}^{w_l\times h_l\times B}$ denotes the observed low spatial resolution HSI, D is the downsampling operator, S is the blurring operator, $\mathcal{X}\in {\mathbb{R}}^{w_h\times h_h\times B}$ is the HR image to be reconstructed and $\mathcal{E}$ represents the observation noise with the same size of $\mathcal{Y}$, and can be considered as the additive Gaussian noise with zero-mean and variance ${\sigma }^2$. It is easy to see that this is an ill-posed problem ($w_l{h}_l<w_h{h}_h$), and some regularization terms of $\mathcal{X}$ based on its prior knowledge, denoted by $R\left(\mathcal{X}\right)$, can be introduced to make this problem to be well-posed, leading to:

$$\widehat{\mathcal{X}}=arg\underset{\mathcal{X}}{min}{\{\parallel DS\mathcal{X}-\mathcal{Y}\parallel }^2+\lambda R\left(\mathcal{X}\right)\},$$

(5)

where $\lambda$ is a regularization parameter to balance the fidelity term and the regularization term. Then, the main objective is to design suitable regularization terms to explore the underlying priors of $\mathcal{X}$.

3.2. 3D TV Regularization

Total variation (TV) has been widely used to explore the spatial piecewise smooth structure for tackling various HSI processing tasks [15,28,29]. It has the ability of preserving local spatial consistency and suppressing observed noise. Considering that the high spatial resolution HSI to be reconstructed is treated as a 3rd-tensor, its local spatial-and-spectral smoothness can be exploited by a weighted 3D TV term, which has the form:

$$\mathrm{TV}\left(\mathcal{X}\right)=\sum _{ijk}{w}_1|x_{ijk}-x_{ij,k-1}|+w_2|x_{ijk}-x_{i,j-1,k}|+w_3|x_{ijk}-x_{i-1,j,k}|,$$

(6)

where $x_{ijk}$ is the $(i,j,k)$-th element of $\mathcal{X}$, and $w_j(j=1,2,3)$ is the weight along the j-th mode of $\mathcal{X}$ that controls its regularization strength. This weighted 3D TV models the piecewise smooth structure in both spatial and spectral domains. For simplicity, we fix the value of $w_j(j=1,2)$ to 1 since the spatial dimensions have similar effect, and tune the spectral weight $w_3$ between 0 to 1. We find that the reconstruction performance is stable when $0.4\le w_3\le 0.8$. As such, we fix the value of $w_3$ to 0.6 in the experimental study presented in this paper.

3.3. Nonlocal Low-Rank Tensor Approximation

Nonlocal self-similarity [30,31,32] is a patch-based useful prior, which means that, for a given local patch in one image, there are many patches similar to it. Motivated by [33], separating $\mathcal{X}$ into a set of image patches $\mathsf{\Omega }={\{{\mathcal{X}}_n\in {\mathbb{R}}^{b\times b\times B}\}}_{p=1}^P$ (where b is the patch size, P is the number of 3D patches with overlap), and by performing block matching [34], a group of patches that is most similar to each patch ${\mathcal{X}}_p$ can be extracted. By stacking all these patches together, we can get a clustered 4th-order tensor ${\mathcal{T}}_k$ with size $b\times b\times B\times d$, where d is the number of 3D patches in the kth cluster. The global spectral correlation among all the bands implies that ${\mathcal{T}}_{k\left(3\right)}$ is low-rank. In addition, the observation that the patches in each cluster possess very similar structures implies that ${\mathcal{T}}_{k\left(4\right)}$ is also low-rank. Combining these two points with the certain correlations in both spatial modes, we can measure the low-rank structure of each 4th-order tensor ${\mathcal{T}}_k$ by a weighted sum of the rank of each unfolding, that is,

$$\mathrm{Rank}\left({\mathcal{T}}_k\right)=\sum _{i=1}^4{\alpha }_i\mathrm{rank}\left({\mathcal{T}}_{k\left(i\right)}\right),$$

(7)

where ${\alpha }_i\ge 0$ and satisfies ${\sum }_{i=1}^4{\alpha }_i=1$. Since minimizing the rank function (7) is generally intractable, and the matrix nuclear norm is usually used as a tight convex surrogate of the matrix rank [35], the rank function (7) can be replaced with the tensor nuclear norm [27]:

$$\parallel {\mathcal{T}}_k{\parallel }_{*}=\sum _{i=1}^4{\alpha }_i{\parallel {\mathcal{T}}_{k\left(i\right)}\parallel }_{*},$$

(8)

where ${\parallel Z\parallel }_{*}:={\sum }_{r=1}^{min(m,n)}{\sigma }_r\left(Z\right)$ is so-called the nuclear norm of matrix Z with size $m\times n$.

Although a convex tensor nuclear norm, such as (8), could provide satisfactory results in various tensor recovery problems, studies like [35] demonstrated that the matrix nuclear norm over-penalizes large singular values, and thus gives a biased estimator in low-rank structure learning. Fortunately, a folded-concave penalty [26] can be considered to remedy such modeling bias [25,26]. Thus, we use the tensor MCP penalty defined by (3) to exploit the low-rank structure of ${\mathcal{T}}_k$, namely,

$$\parallel {\mathcal{T}}_k{\parallel }_{P_{\lambda }}=\sum _{i=1}^4{\alpha }_i{\parallel {\mathcal{T}}_{k\left(i\right)}\parallel }_{P_{\lambda }},$$

(9)

where $P_{\lambda }$ is the so-called MCP penalty function defined in (2).

Thus far, we have only considered how to model the low-rank characteristic of one cluster $T_k$. Basically, it should be necessary to exploit the spatial-and-spectral correlation of $\mathcal{X}$ by the following form:

$${\parallel \mathcal{X}\parallel }_{Nonlocal-P_{\lambda }}=\sum _{k=1}^N{\parallel {\mathcal{T}}_k\parallel }_{P_{\lambda }},$$

(10)

where N is the number of clustered groups.

3.4. Proposed Model

With the above discussions, we now propose the following regularization model for the HSI super-resolution task:

$$\underset{\mathcal{X}}{min}{\parallel DS\mathcal{X}-\mathcal{Y}\parallel }_F^2+{\parallel \mathcal{X}\parallel }_{Nonlocal-P_{{\lambda }_1}}+{\lambda }_2\mathrm{TV}\left(\mathcal{X}\right),$$

(11)

where the scalars ${\lambda }_1$ and ${\lambda }_2$ are regularization parameters. It is not hard to see that (11) is a nonconvex optimization problem due to the nonconvexity of MCP penalty function. Thus, we only wish to find a good local solution. Specifically, as is shown in Figure 1, by choosing the direct upsampled version of observered LR tensor $\mathcal{Y}$ as the intilization point, we would like to estimate a good HR tensor $\mathcal{X}$ via the ADMM.

4. Optimization Procedure

The model (11) can be rewritten as the following equivalent form:

$$\underset{\mathcal{X}}{min}{\parallel \mathcal{Y}-DS\mathcal{X}\parallel }^2+\sum _k\sum _{i=1}^4{\alpha }_i{\parallel {\mathcal{T}}_{k\left(i\right)}\parallel }_{P_{{\lambda }_1}}+{\lambda }_2{\parallel G_w\left(\mathcal{X}\right)\parallel }_1,$$

(12)

where $G_w(·)=\left[w_1\times G_h(·);w_2\times G_v(·);w_3\times G_t(·)\right]$ is the so-called weighted three-dimensional difference operator, and $G_h$, $G_v$, $G_t$ are the first-order difference operators with respect to three different directions of HSI cube. Define the operator $\mathcal{P}$ that first extracts patches within a 3rd-order HSI cube and then applies nonlocal block matching on these 3rd-order patches to form a 4th-order tensor. Thus, we can denote that ${\mathcal{P}}_{\mathcal{X}}^k={\mathcal{T}}_k$. By introducing some auxiliary variables, i.e., ${\left\{\mathcal{M}\right\}}_{i=1}^4$, $\mathcal{F}$ into (12) leads to

$$\begin{array}{c}\underset{\mathcal{X},{{\{\mathcal{M}}_n\}}_{i=1}^4,\mathcal{F}}{min}{∥\mathcal{Y}-DS\mathcal{X}∥}^2+\sum _k\sum _{i=1}^4{\alpha }_i{\parallel {{\mathcal{P}}_{{\mathcal{M}}_i}^k}_{\left(i\right)}\parallel }_{P_{{\lambda }_1}}+{\lambda }_2{\parallel \mathcal{F}\parallel }_1\hfill \\ \mathrm{subject}\mathrm{to}\mathcal{X}={\mathcal{M}}_i,\mathcal{F}=G_w\left(\mathcal{X}\right),i=1,\dots ,4.\hfill \end{array}$$

(13)

Based on the methodology of ADMM [36], the augmented Lagrangian function of the above problem is written as follows:

$$\begin{array}{c}{L}_{\mu ,\nu }(\mathcal{X},{\mathcal{M}}_1,\dots ,{\mathcal{M}}_4,{\mathcal{U}}_1,\dots ,{\mathcal{U}}_4,\mathcal{F},\mathcal{V})\hfill \\ ={\parallel \mathcal{Y}-DS\mathcal{X}\parallel }_F^2+\sum _{i=1}^4\left(\sum _k{\alpha }_i{\parallel {{\mathcal{P}}_{{\mathcal{M}}_i}^k}_{\left(i\right)}\parallel }_{P_{{\lambda }_1}}+\langle \mathcal{X}-{\mathcal{M}}_i,{\mathcal{U}}_i\rangle +\frac{\mu }{2}{∥\mathcal{X}-{\mathcal{M}}_i∥}_F^2\right)\hfill \\ +{\lambda }_2{\parallel \mathcal{F}\parallel }_1+\langle \mathcal{F}-G_w\mathcal{X},\mathcal{V}\rangle +\frac{\nu }{2}{\parallel \mathcal{F}-G_w\mathcal{X}\parallel }_F^2,\hfill \end{array}$$

(14)

where ${\left\{{\mathcal{U}}_i\right\}}_{i=1}^4$ and $\mathcal{V}$ are the Lagrangian parameters. We will transform (14) into three subproblems and iteratively update each variable by fixing the other ones. Let l denote the lth iteration step; then, in the $(l+1)$th iteration, variables involved in model (13) can be updated as follows:

Update $\mathcal{X}$. Extracting all items related to $\mathcal{X}$ from (14), it can be easily deduced that:

$$\begin{array}{cc}\hfill {\mathcal{X}}^{(l+1)}=\underset{\mathcal{X}}{argmin}{∥\mathcal{Y}-DS\mathcal{X}∥}_F^2& +\sum _{i=1}^4\left(\langle \mathcal{X}-{\mathcal{M}}_i^{\left(l\right)},{\mathcal{U}}_i^{\left(l\right)}\rangle +\frac{\mu }{2}{∥\mathcal{X}-{\mathcal{M}}_i^{\left(l\right)}∥}_F^2\right)\hfill \\ & +\langle {\mathcal{F}}^{\left(l\right)}-G_w\mathcal{X},{\mathcal{V}}^{\left(l\right)}\rangle +\frac{\nu }{2}{\parallel {\mathcal{F}}^{\left(l\right)}-G_w\mathcal{X}\parallel }_F^2,\hfill \end{array}$$

(15)

which equals solving the following linear equation:

$$2\tilde{S}\tilde{D}DS\mathcal{X}+4\mu \mathcal{X}+\nu G_w^{\prime }{G}_w\mathcal{X}=2\tilde{S}\tilde{D}\mathcal{Y}+\sum _{i=1}^4(\mu {\mathcal{M}}_i^{\left(l\right)}-{\mathcal{U}}_i^{\left(l\right)})+G_w^{\prime }{\mathcal{V}}_n^l+\nu G_w^{\prime }{\mathcal{F}}_n^l,$$

(16)

where $\tilde{S}$, $\tilde{D}$ are the transposes of S and D, respectively, and $G_w^{\prime }$ indicates the adjoint of $G_w$. This linear system can be efficiently solved by using the well-known preconditioned conjugate gradient method.

Update ${\left\{{\mathcal{M}}_i\right\}}_{i=1}^4$. Similar to the works [22,25], we adopt the LLA algorithm to transform the MCP penalization problem into a weighted nuclear norm penalization problem. Specifically, we need to solve

$$\underset{{\left\{{\mathcal{M}}_i\right\}}_{i=1}^4}{min}\sum _{i=1}^4\left(\sum _k{\alpha }_i{Q}_{P_{{\lambda }_1}}\left(\left.\sigma \left({{\mathcal{P}}_{{\mathcal{M}}_i}^k}_{\left(i\right)}\right)\right|\sigma \left({{\mathcal{P}}_{{\mathcal{X}}^{(l+1)}}^k}_{\left(i\right)}\right)\right)+\langle {\mathcal{X}}^{(l+1)}-{\mathcal{M}}_i,{\mathcal{U}}_i^{\left(l\right)}\rangle +\frac{\mu }{2}\parallel {\mathcal{X}}^{(l+1)}-{\mathcal{M}}_i{\parallel }_F^2\right),$$

(17)

where $Q_{P_{{\lambda }_1}}\left(\sigma (X|X^{\left(l\right)})\right)$ is the locally linear approximation of ${\parallel X\parallel }_{P_{\lambda }}$ for a given matrix X when $X^{\left(l\right)}$ is given. For solving (17), we can first consider the following problem for each kth patch set:

$$\underset{{\left\{{\mathcal{P}}_{{\mathcal{M}}_i}^k\right\}}_{i=1}^4}{min}\sum _{i=1}^4\left({\alpha }_i{Q}_{P_{{\lambda }_1}}\left(\left.\sigma \left({{\mathcal{P}}_{{\mathcal{M}}_i}^k}_{\left(i\right)}\right)\right|\sigma \left({{\mathcal{P}}_{{\mathcal{X}}^{(l+1)}}^k}_{\left(i\right)}\right)\right)+\langle {\mathcal{P}}_{{\mathcal{X}}^{(l+1)}}^k-{\mathcal{P}}_{{\mathcal{M}}_i}^k,{\mathcal{P}}_{{\mathcal{U}}_i^{\left(l\right)}}^k\rangle +\frac{\mu }{2}{∥{\mathcal{P}}_{{\mathcal{X}}^{(l+1)}}^k-{\mathcal{P}}_{{\mathcal{M}}_i}^k∥}_F^2\right).$$

(18)

Then, according to [27], we can update each ${\mathcal{P}}_{{\mathcal{M}}_i}^k$ by solving

$$\underset{{\left\{{\mathcal{P}}_{{\mathcal{M}}_i}^k\right\}}_{i=1}^4}{min}{\alpha }_i{Q}_{P_{{\lambda }_1}}\left(\left.\sigma \left({{\mathcal{P}}_{{\mathcal{M}}_i}^k}_{\left(i\right)}\right)\right|\sigma \left({{\mathcal{P}}_{{\mathcal{X}}^{(l+1)}}^k}_{\left(i\right)}\right)\right)+\langle {\mathcal{P}}_{{\mathcal{X}}^{(l+1)}}^k-{\mathcal{P}}_{{\mathcal{M}}_i}^k,{\mathcal{P}}_{{\mathcal{U}}_i^{\left(l\right)}}^k\rangle +\frac{\mu }{2}{∥{\mathcal{P}}_{{\mathcal{X}}^{(l+1)}}^k-{\mathcal{P}}_{{\mathcal{M}}_i}^k∥}_F^2.$$

(19)

Thus, by applying the Theorem 3.1 in [25], we can get the following closed-form solution:

$${\mathcal{P}}_{{\mathcal{M}}_i^{(l+1)}}^k={\mathrm{fold}}_i\left[S_{{\lambda }_1{\alpha }_i/\mu ,W_i}\left({{\mathcal{P}}_{{\mathcal{X}}^{(l+1)}}^k}_{\left(i\right)}+\frac{{{\mathcal{P}}_{{\mathcal{U}}_i^{\left(l\right)}}^k}_{\left(i\right)}}{\mu }\right)\right],i=1,2,3,4.$$

(20)

Here, the singular value shrinkage operator $S_{\tau }\left(X\right)$ is defined by $S_{\tau }\left(X\right):=U_X{D}_{\tau }\left({\mathsf{\Sigma }}_X\right)V_X^T$, where $X=U_X{\mathsf{\Sigma }}_X{V}_X^T$ is the singular value decomposition of the matrix X and for a given matrix A, ${\left[D_{\tau }\left(A\right)\right]}_{ij}=sgn\left(A_{ij}\right)(|A_{ij}{|-\tau )}_{+}$. Additionally, the weight matrix $W_i$ is defined by $W_i=Diag\left({(\lambda -(\sigma \left({\mathcal{X}}_{\left(i\right)}^{(l+1)}\right)/a))}_{+}\right)$ for some fixed $a>1$. Then, as shown in Figure 1, by aggregating all such patch sets ${\{{\mathcal{P}}_{{\mathcal{M}}_i^{(l+1)}}^k\}}_k$, we can update each ${\mathcal{M}}_i$ as

$${\mathcal{M}}_i^{(l+1)}=\sum _k{\mathcal{P}}_{{\mathcal{M}}_i^{(l+1)}}^k.$$

(21)

Update $\mathcal{F}$. Extracting all terms containing $\mathcal{F}$ from (14), we can get

$$\begin{array}{cc}\hfill {\mathcal{F}}^{(l+1)}& =\underset{\mathcal{F}}{argmin}\sum _{i=1}^4\left({\lambda }_2{\beta }_i{∥{\mathcal{F}}_{\left(i\right)}∥}_1+\langle \mathcal{F}-G_w{\mathcal{X}}^{(l+1)},{\mathcal{V}}^{\left(l\right)}\rangle +\frac{\nu }{2}{∥\mathcal{F}-G_w{\mathcal{X}}^{(l+1)}∥}_F^2\right)\hfill \\ & =\underset{\mathcal{F}}{argmin}\sum _{i=1}^4\left({\lambda }_2{\beta }_i{∥{\mathcal{F}}_{\left(i\right)}∥}_1+\frac{\nu }{2}\parallel \mathcal{F}-G_w{\mathcal{X}}^{(l+1)}+\frac{{\mathcal{V}}_{\left(i\right)}^{\left(l\right)}}{\nu }{\parallel }_F^2\right).\hfill \end{array}$$

(22)

By introducing the so-called soft-thresholding operator:

$${Soft}_{\Delta }\left(x\right)=\left\{\begin{array}{cc}\hfill x-\Delta ,& ifx>\Delta ,\\ \hfill x+\Delta ,& ifx<\Delta ,\\ \hfill 0,& otherwise,\end{array}\right.$$

where $x\in \mathbb{R}$ and $\Delta >0$, and then we update ${\mathcal{F}}^{(l+1)}$ as

$${\mathcal{F}}^{(l+1)}={\mathrm{fold}}_i\left[{Soft}_{{\lambda }_2{\beta }_i/\nu }\left(G_w{\mathcal{X}}_{\left(i\right)}^{{}^{(l+1)}}-\frac{{\mathcal{V}}_{\left(i\right)}^{\left(l\right)}}{\nu }\right)\right].$$

(23)

Update ${\left\{{\mathcal{U}}_i\right\}}_{i=1}^4$ and $\mathcal{V}$. The multipliers are updated by

$$\begin{array}{c}{\mathcal{U}}_i^{(l+1)}={\mathcal{U}}_i^{\left(l\right)}+\gamma ·\mu ({\mathcal{X}}^{(l+1)}-{\mathcal{M}}_i^{(l+1)}),i=1,\dots ,4,\hfill \\ {\mathcal{V}}^{(l+1)}={\mathcal{V}}^{\left(l\right)}+\gamma ·\nu ({\mathcal{F}}^{(l+1)}-G_w{\mathcal{X}}^{(l+1)}),\hfill \end{array}$$

(24)

where $\gamma$ is the parameter with a fixed value, namely, 1.05, and the penalty parameters $\mu$ and $\nu$ follow a certain adaptive updating scheme, which can facilitate the convergence of the proposed optimization procedure. Take $\mu$ as an example. We first initialize $\mu$ as a small value, i.e., $\mu ={10}^{-3}$, and then update it by the scheme:

$$\mu \leftarrow c_1·\mu \mathrm{if}Res>c_2·Re{s}_{pre},$$

(25)

where $Res=\parallel \mathcal{Y}-DS{\mathcal{X}}^{(l+1)}{\parallel }_F$ and $Re{s}_{pre}={\parallel \mathcal{Y}-DS{\mathcal{X}}^{\left(l\right)}\parallel }_F$, and $c_1$ and $c_2$ can be taken as 1.15 and 0.95, respectively.

5. Experimental Study

Three popular real-world HSI data sets were used in our experiments, i.e., the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) Moffett Field data set [37], the Hyperspectral Digital Imagery Collection Experiment (HYDICE) Urban data set [38] and the HYDICE Washington DC Mall data set [39]. After removing seriously polluted bands and cropping images for each data set, the HSI cube used for the experiments are of $256\times 256\times 146$, $256\times 256\times 140$ and $256\times 256\times 140$, respectively. To thoroughly evaluate the performance of the proposed approach, we considered three popular super-resolution methods for comparison, that is, the nonlocal autoregressive model (NARM) proposed by [21], the spatial–spectral group sparsity method (SSGS) proposed by [19], the low-rank and total variation regulariztions (LRTV) method proposed by [23]. We also considered the nearest neighbor interpolation (NN) method that is used to achieve the upsampled HSI for comparison. For brevity, our proposed approach is dubbed as NLRTATV. In this set of experiments, the blurring kernel is chosen as the popular Gaussian kernel and all the LR HSIs are obtained by downsampling the original HR HSIs with a factor of 2 or 3, i.e., the LR HSIs are of spatial size $128\times 128$ or $85\times 85$. In the following experiments, similar to [19], the gray values of each band of HSI were normalized to [0, 255] to facilitate the numerical calculation, though this operation may change the relative spectral properties of the HSI bands. In addition, the parameters a and ${\alpha }_i$ in (2) and (3) were fixed to 5 and $\frac{1}{4}$, respectively. For another two regularization parameters ${\lambda }_1$ and ${\lambda }_2$ of the model (12), though they should be carefully tuned, based on the analysis presented in Section 5.3, we set ${\lambda }_1$ as 0.3 and ${\lambda }_2$ as 0.04. For the ease of interested readers to reproduce the results of NLRTATV, we have released the codes at http://vision.sia.cn/our%20team/Hanzhi-homepage/vision-ZhiHan%28English%29.html.

5.1. Quantitative Comparison

Table 1. Super-resolution reconstruction results by different methods under noise free setting.

		NN	NARM	SSGS	LRTV	NLRTATV
	MPSNR (dB)	19.6357	20.1774	21.2733	27.7744	34.3935
DC Mall (factor 2)	MSSIM	0.6771	0.7575	0.8142	0.8380	0.9665
	SAM (rad)	0.0998	0.0913	0.0875	0.0805	0.0487
	MPSNR (dB)	20.1281	21.5613	22.7087	28.2954	35.2754
Urban (factor 2)	MSSIM	0.6924	0.7770	0.8092	0.8485	0.9663
	SAM (rad)	0.0971	0.0837	0.0789	0.0710	0.0470
	MPSNR (dB)	19.9535	20.3799	21.4042	29.8499	35.9816
Moffett Field (factor 2)	MSSIM	0.6836	0.7621	0.8042	0.8349	0.9545
	SAM (rad)	0.0958	0.0820	0.0754	0.0630	0.0365
	MPSNR (dB)	24.9190	25.1116	26.0086	27.4270	33.1308
Moffett Field (factor 3)	MSSIM	0.6888	0.7130	0.7139	0.7955	0.9301
	SAM (rad)	0.0963	0.0915	0.0901	0.0810	0.0476

and

Table 2. Super-resolution reconstruction results by different methods under noisy setting ($\sigma =5$).

		NN	NARM	SSGS	LRTV	NLRTATV
	MPSNR (dB)	19.4578	20.1558	21.2191	27.2824	32.8701
DC Mall (factor 2)	MSSIM	0.6346	0.7575	0.7901	0.8126	0.9329
	SAM (rad)	0.1034	0.0913	0.0889	0.0845	0.0473
	MPSNR (dB)	20.8761	21.4544	21.6906	27.7523	31.8675
Urban (factor 2)	MSSIM	0.6263	0.7770	0.8090	0.8148	0.9042
	SAM (rad)	0.0928	0.0897	0.0835	0.0745	0.0628
	MPSNR (dB)	19.7640	20.3251	20.4034	29.1580	32.5487
Moffett Field (factor 2)	MSSIM	0.6107	0.7621	0.7982	0.8158	0.8834
	SAM (rad)	0.0905	0.0720	0.0695	0.0635	0.0470
	MPSNR (dB)	24.4130	25.0118	25.0782	26.8280	29.4416
Moffett Field (factor 3)	MSSIM	0.6223	0.7130	0.7247	0.7458	0.8441
	SAM (rad)	0.0929	0.0915	0.0864	0.0839	0.0660

give the super-resolution reconstruction results by all the compared methods under noise free setting and noisy setting with $\sigma =5$, respectively, in terms of three image quality measures, including the mean peak signal-to noise ratio (MPSNR) and the mean feature similarity (MSSIM), as used in [40], and spectral angle mapper (SAM) [41]. The PSNR (in dB) [42] and SSIM [43] are two traditional image quality measures in image processing and computer vision, which are used to evaluate the similarity between the target image and the reference image based on mean squared error and structural consistency, respectively. Note that, despite its popularity, PSNR may not be a very good image quality measure for HSI based problems. SAM (in rad) calculates the average angle between spectrum vectors of the target HSI and the reference one across all spatial positions. It is well-known that the higher the value of PSNR or SSIM is, the better the image quality is; the lower the value of SAM is, the smaller the spectral distortion.

As shown in Table 1 and Table 2, firstly, the proposed NLRTATV method significantly outperforms other competing ones in terms of the aforementioned three quality measures. This clearly shows that more fine priors underlying the HSI cube were exploited by NLRTATV. Secondly, LRTV illustrates a superior performance over NN, SSGS and NARM because of the use of a direct tensor sparse representation of the hyperspectral cube. Both of these findings demonstrate the power of tensor modeling techniques for tackling HSI super-resolution tasks.

To provide more detailed quantitative comparison, we also show in Figure 2 and Figure 3 the PSNR and SSIM values of each band for all four of the cases listed in Table 1 and Table 2, respectively. As can be seen, the tensor methods including LRTV and the proposed NLRTATV can get much higher values of SSIM and PSNR than other ones for almost every band. In addition, NLRTATV can provide significantly improved reconstruction results over LRTV, especially for noise free cases. All of these show that the proposed NLRTATV can discover more intrinsic structures of HSIs than other competing ones, and thus it possesses a better reconstruction ability.

5.2. Visual Quality Comparison

In terms of visual quality, one representative band of recovered HSIs in three typical cases obtained by all of the compared methods are shown in Figure 4, Figure 5 and Figure 6. It is evident that the images shown by (f) in the Figure 4, Figure 5 and Figure 6 recovered by our NLRTATV method preserve clearer and sharper edges compared with the other images shown by (b–e) in Figure 4, Figure 5 and Figure 6. One can zoom in on each figure for more visual comparison details.

5.3. Analysis of the Parameters

In the proposed NLRTATV model, there exist two regularization parameters ${\lambda }_1$ and ${\lambda }_2$ that need to be carefully identified. Thus, we provide some experimental analysis to show how we can choose them for practical use. (a) and (b) in Figure 7 describe the relationship between MPSNR and the regularization parameters ${\lambda }_1$ and ${\lambda }_2$ when NLRTATV performed on DC Mall data under noise free setting, respectively, with the other parameters fixed at optimal values. It can be observed that, as ${\lambda }_1$ and ${\lambda }_2$ vary in a wide range, the MPSNR value is rather stable and relatively stable correspondingly. Actually, we performed NLRTATV on several other experimental cases and found a similar behaviour as shown in Figure 7. Therefore, ${\lambda }_1$ can be chosen as a value in $[0.1,1]$ and ${\lambda }_2$ can be chosen as a value in $[0.005,0.05]$ in real implementation. For simplicity, in all the experiments conducted in this work, we set ${\lambda }_1$ and ${\lambda }_2$ as 0.3 and 0.04, respectively.

Since the nonlocal block matching operation is used in NLRTATV, the number of similar patches in each group is another important parameter that needs to be carefully considered. Figure 8 depicts the MPSNR and MSSIM gains versus the number of similar patches when NLRTATV performed on DC Mall data under a noise free setting. Here, we can observe that, as the number of similar patches increases to a value larger than 10, the curves of MPSNR and MSSIM both became stable. For simplicity, in all of the experiments, the number of similar patches was fixed to 30.

6. Conclusions

In this paper, we have proposed a novel method named NLRTATV for dealing with the problem of HSI super-resolution by using tensor structural modeling. The proposed method has considered the global correlation across spectral domain, the nonlocal self-similarity across spatial domain, and the local smooth structure across both spatial and spectral domains of the HSI cube by combining the low-rank tensor and 3D total variation regularizations. Extensive experimental study on three HSI datasets have demonstrated the superior performance of the proposed NLRTATV over several popular methods, which clearly shows the advantages of NLRTATV in exploring the intrinsic characteristics of the HSI cube.

Along this line, several desirable research directions can be conducted for future study. Firstly, though the proposed NLRTATV has only tested on the experiments with Gaussian blurring kernel, which is the situation in which no prior information on the blurring operator is known, a similar idea used in [18] could be incorporated into our model for tackling HSI super-resolution tasks with unknown blurring kernels. Secondly, it could be interesting to extend the deep learning ideas of [44] and [45] to design a deep tensor architecture to learn the more intrinsic characteristics of the HSI cube, and, as a result, to provide much better super-resolution reconstruction results. Thirdly, the problem of hyperspectral and multispectral data fusion has received much attention in recent years. It would be also interesting to extend the presented tensor modeling technique for dealing with such problems, which could give more effective methods for enhancing the spatial resolution of hyperspectral data.