Hyperspectral super-resolution via low rank tensor triple decomposition

Equipment for hyperspectral imaging uses an unusual sensor to produce images with great spectral resolution. The key factor in an image is the electromagnetic spectrum obtained by sampling at hundreds of continuous wavelength intervals. Because hyperspectral images have hundreds of bands, spectral information is abundant [32]. Hyperspectral images are now the core of a large and increasing number of remote sensing applications, such as environmental monitoring [4, 25], target classification [41, 8] and anomaly detection [7, 37].

With optical sensors, there is a fundamental compromise among spectral resolution, spatial resolution and signal-to-noise ratio. Hyperspectral images therefore have a high spectral resolution but a low spatial resolution. Contrarily, multispectral sensors, which only have a few spectral bands, can produce images with higher spatial resolution but poorer spectral resolution [5].

The approaches for addressing hyperspectral and multispectral images fusion issues can be categorized into four classes [42]. The pan-sharpening approach falls within the first category, which includes component substitution [27, 1] and multi-resolution analysis [6]. Selva [30] and Liu [22] proposed multi-resolution analysis method to extract spatial details that were injected into the interpolated hyperspectral image. Chen et al. [9] divided all bands of the images into several regions and then fused the images of each region by the pan-sharpening approach. Unfortunately, these methods may lead to spectral distortion. The second category is subspace-based formulas. The main idea is to explore the low-dimensional representation of hyperspectral images, which are regarded as a linear combination of a set of basis vectors or spectral features [18, 38, 34]. Hardie et al. [15] employed Bayesian formulas to solve the image fusion problem with the maximum posteriori probability. The method based on matrix factorization falls within the third category. Under these circumstances, HSI and MSI are regarded as the degradation versions of the reference SRI. Yokoya et al. [40] proposed a coupled nonnegative matrix factorization method, which estimated the basic vectors and its coefficients blindly and then calculated the endmember matrix and degree matrix by alternating optimization algorithm. In addition, there are methods combining sparse regularization and spatial regularization [19, 31], such as a low rank matrix decomposition model for decomposing a hyperspectral image into a super-resolution image and a difference image in [12, 42].

Tensor decomposition methods fall within the fourth category. Tensor, a higher-order extension of matrix, can capture the correlation between the spatial and spectral domains in hyperspectral images simultaneously. Representing hyperspectral and multispectral images as third-order tensors has been widely used in hyperspectral image denoising and super-resolution problems. The hyperspectral image is a three-dimensional data cube with two spatial dimensions (height and width) and one spectral dimension [17, 33, 23]. The spatial details and spectral structures of the SRI therefore are both more accurate.

By utilizing low rank tensor decomposition technology to explore the low rank properties of hyperspectral images, high spatial resolution and accurate restoration guarantees are obtained. Kanatsoulis et al. [20] proposed treating the super-resolution problem as a coupled tensor approximation, where the super-resolution image satisfied low rank Canonical Polynomial (CP) decomposition. Later, Li et al. [28, 24] considered a novel approach based on tucker decomposition. Xu et al. [39] first considered that the target SRI exhibited both the sparsity and the piecewise smoothness on three modes and then used tucker decomposition to construct model. Dian et al. [13] applied tensor kernel norm to explore low rank property of SRI. He et al. [16] proposed tensor ring decomposition model for addressing the image fusion problem.

In this paper, we propose a novel model based on low rank tensor triple decomposition to address the hyperspectral and multispectral images fusion problem. Triple decomposition breaks down a third-order tensor into three smaller tensors with predetermined dimensions. It is effectively applied in tensor completion and image restoration [29, 11]. The structures of image fusion and triple decomposition are presented in Figure 1. The optimization model is solved by the L-BFGS algorithm. The convergent results are given with the help of Kurdyka-$\L$ojasiewicz property of the objective function. Numerical experiments demonstrate that the proposed TTDSR method performs well and the convergent conclusion is also validated numerically.

The rest of this paper is organized as follows. In section 2, we first provide definitions used in the paper and present the basic algebraic operations of tensors. In section 3, we introduce our motivation and model for the hyperspectral and multispectral images fusion. In section 4, the L-BFGS algorithm and the gradient calculation are also presented in detail. Convergence analysis of the algorithm is available in section 5. In section 6, experimental results on two datasets are given. The conclusion is given in section 7.

Let $\mathbb{R}$ be the real field. We use lowercase letters and boldface lowercase letters to represent scalars and vectors respectively, while capital letters and calligraphic letters stand for matrices and tensors respectively.

An $m$th order $n$ dimensional tensor $\mathcal{A}$ is a multidimensional array with its entry being

Similarly, the $i$th entry in a vector a is symbolized by $(\mathbf{a})_{i}=a_{i}$, and the $(i,j)$th entry of a matrix $A$ is $(A)_{ij}=a_{ij}.$ Unless otherwise specified, the order of tensor is set to three hereafter in this paper. Without loss of generality, we denote $\mathcal{A}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}.$ A fiber of $\mathcal{A}$ is a vector produced by fixing two indices of $\mathcal{A},$ such as $\mathcal{A}(i,j,:)$ for any $i\in[1,\ldots,n_{1}],j\in[1,\ldots,n_{2}]$. A slice of $\mathcal{A}$ is a matrix by varying two of its indices while fixing another one, such as $\mathcal{A}(:,:,k)$ for any $k\in[1,\ldots,n_{3}].$ The geometric structures of the tensor expansion are shown in Figure 2.

The mode-$k$ product of the tensor $\mathcal{A}$ and a matrix $X$ is an extension of matrix product. Suppose matrices $F\in\mathbb{R}^{l\times n_{1}}$, $G\in\mathbb{R}^{m\times n_{2}}$, and $H\in\mathbb{R}^{n\times n_{3}}$. The mode-$1$ product of $\mathcal{A}$ and $F$ is denoted as $\mathcal{A}\times_{1}F$ with its elements being

Also we have

for $p=1,\ldots m,$ $q=1,\ldots,n.$ It is easy to verify that

The mode-$k$ matricization of $\mathcal{A}$ denoted by ${A}_{(k)}$ arranges its $n_{k}$ mode-$k$ fibers as the columns of $A_{(k)}$ in order. The Frobenius norm of tensor $\mathcal{A}$ is given by

The unfolding related no matter with tensors, matrices or vectors in the following part is carried out under the precondition that the left index varies more rapidly than the right one. Subscript represents dimension. We note that ${vec}\left(B_{m\times n}\right)$ creates a column vector $\mathbf{b}_{mn}$ by stacking the columns of $B_{m\times n}$ below one another, i.e., $\left(\mathbf{b}_{mn}\right)_{i+(j-1)m}=b_{ij}$. In addition, the kronecker product of matrices $B_{m\times n}$ and $C_{p\times q}$ is

for $s=1,\cdots,p,$ $t=1,\cdots,q.$ It generates a large matrix $D_{pm\times qn}=B_{m\times n}\otimes C_{p\times q}$, whose entries are ${d}_{gh}=b_{ij}c_{st}$ with $g=s+(i-1)p$ and $h=t+(j-1)q$. By this rule, unfolding matrices of tensor $\mathcal{A}$ could be represented by

Moreover, $A_{n_{1}n_{2}\times n_{3}}=\left(A_{n_{3}\times n_{1}n_{2}}\right)^{\top}.$

Tensors can describe high dimensional relationship. The hyperspectral and multispectral images are naturally third-order tensors. Methods based on low rank matrix decomposition reconstruct the 3-d tensor data into a 2-d matrix data at beginning, which destroy the original 3-d structure and may cause negative influence on the fusion consequences. Therefore, hyperspectral and multispectral image fusion problems have been investigated with the aid of various tensor decompositions in recent years. The tensor triple decomposition model has certain advantages over the classical CP and tucker decompositions. The triple rank is less than or equal to CP rank, which has been proved theoretically and practically. As described in theorem 2.1 of [29], low rank triple decomposition and triple ranks are well-defined. The proposed triple decomposition has the low rank properties and performs well in image restoration. When facing large-scale problems, the lower rank indicates that we transform the problem into a lower dimensional subspace, which is relatively easy to solve. Thus, we employ the tensor triple decomposition framework to establish the hyperspectral super-resolution model.

Define $\mathcal{Y}_{h}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, $\mathcal{Y}_{m}\in\mathbb{R}^{m_{1}\times m_{2}\times m_{3}}$ and $\mathcal{Z}\in\mathbb{R}^{m_{1}\times m_{2}\times n_{3}}$ as hyperspectral image (HSI), multispectral image (MSI) and super-resolution image (SRI) respectively. Here the first two indices $n_{1}$ and $n_{2}$ or $m_{1}$ and $m_{2}$ denote the spatial dimensions. The third index $n_{3}$ or $m_{3}$ indicate the number of spectral bands. Usually, the MSI contains more spatial information than the HSI, that is $m_{1}\gg n_{1},m_{2}\gg n_{2}.$ The HSI has more spectral bands than MSI, which means $n_{3}\gg m_{3}.$ Our target is to find the SRI that possesses the high spectral resolution of HSI and the spatial resolution of the MSI, i.e. a super-resolution image.

In this section, we focus on numerical approaches for computing a first-order stationary point of the optimization problem (12). Since the limited memory BFGS (L-BFGS) algorithm [21] is powerful for large scale nonlinear unconstrained optimization, we apply it to produce a search direction and consider using inexact line search techniques to update the iteration step size. In the computaion process, we either matricize or vectorize the tensor variable to get the gradient of the objective function. For convenience, we demonstrate the algorithm and its convergence analysis for