Hyperspectral Image Super-Resolution Algorithm Based on Graph Regular Tensor Ring Decomposition

Abstract

This paper introduces a novel hyperspectral image super-resolution algorithm based on graph-regularized tensor ring decomposition aimed at resolving the challenges of hyperspectral image super-resolution. This algorithm seamlessly integrates graph regularization and tensor ring decomposition, presenting an innovative fusion model that effectively leverages the spatial structure and spectral information inherent in hyperspectral images. At the core of the algorithm lies an iterative optimization process embedded within the objective function. This iterative process incrementally refines latent feature representations. It incorporates spatial smoothness constraints and graph regularization terms to enhance the quality of super-resolution reconstruction and preserve image features. Specifically, low-resolution hyperspectral images (HSIs) and high-resolution multispectral images (MSIs) are obtained through spatial and spectral downsampling, which are then treated as nodes in a constructed graph, efficiently fusing spatial and spectral information. By utilizing tensor ring decomposition, HSIs and MSIs undergo feature decomposition, and the objective function is formulated to merge reconstructed results with the original images. Through a multi-stage iterative optimization procedure, the algorithm progressively enhances latent feature representations, leading to super-resolution hyperspectral image reconstruction. The algorithm’s significant achievements are demonstrated through experiments, producing sharper, more detailed high-resolution hyperspectral images (HRIs) with an improved reconstruction quality and retained spectral information. By combining the advantages of graph regularization and tensor ring decomposition, the proposed algorithm showcases substantial potential and feasibility within the domain of hyperspectral image super-resolution.

1. Introduction

Hyperspectral remote sensing is a specialized remote sensing technique employed to acquire hyperspectral image data pertaining to the Earth’s surface. Compared with traditional remote sensing images, hyperspectral remote sensing can provide richer spectral information by observing and collecting data in more narrow wavelength bands. It has a wide range of applications in the fields of earth science, environmental monitoring, agriculture, forestry, urban planning, etc. Hyperspectral remote sensing uses sensors capable of collecting data in dozens or hundreds of consecutive spectral bands. These bands typically cover the visible, near-infrared, and short-wave infrared frequency bands. The spectral response of each band can provide a spectral signature, allowing more information to be obtained from remotely sensed images. By analyzing and processing these hyperspectral data, more detailed and accurate information, such as ground features, material composition, and environmental parameters, is obtained. However, the spatial resolution and spectral resolution of hyperspectral images often cannot reach the ideal state at the same time, considering the limitations of solar incident energy and sensors. Therefore, a popular approach is to improve the spatial resolution of hyperspectral images by fusing multispectral images (MSIs) with hyperspectral images (HSIs). This fusion method can make full use of the higher spatial resolution of multispectral images and combine it with the rich spectral information of hyperspectral images to obtain images with both a high spatial resolution and rich spectral information.

In the field of hyperspectral images, image fusion has been extensively and intensively studied by many researchers. According to the three dimensions in which the fused data are located, image fusion can be classified into pixel-level fusion, feature-level fusion, and decision-level fusion. Panchromatic sharpening [1] was one of the first spatial–spectral fusion methods developed. This method fuses a panchromatic image with a high spatial resolution with a multispectral image with a low spatial resolution but with multiple spectral bands. The panchromatic image has only one spectral band but a high spatial resolution, while the multispectral image contains rich spectral information, although it has a low spatial resolution. Via the panchromatic sharpening method, an image with both high spatial resolution containing multiple spectral bands can be obtained. This fusion method has an important role in the field of remote sensing and can provide more accurate and comprehensive information for various applications.

Pixel-level fusion methods are classified as component replacement (CS) [2,3,4], multi-resolution analysis (MRA) [5,6], hybrid methods [7,8], and model-based methods. These methods include intensity hue saturation (IHS) [9,10] methods, the Brovey transform [11], principal component analysis (PCA) [12], Gram–Schmidt (GS) [13] orthogonalization, the Laplace pyramid [14], the curvelet transform [15], and sparse matrix decomposition [16]. However, with the introduction of more algorithms and experiments, researchers have found some limitations of traditional methods to solve the problem of the hyper-resolution of hyperspectral images.

With the continuous development of model-based methods, researchers are increasingly focused on matrix-based decomposition [17,18,19,20,21] and tensor-based decomposition [22,23,24,25] methods. Common tensor decomposition methods include matrix decomposition, CP decomposition (CANDECOMP/PARAFAC) [26,27], Tucker decomposition [28,29], Tensor Train decomposition [30,31], and tensor ring decomposition [32,33,34]. These methods are receiving increasing attention and have more applications in the field of hyperspectral image super-resolution. A matrix can be viewed as a two-dimensional tensor, and therefore matrix decomposition can also be viewed as two-dimensional tensor decomposition.

Naoto Yokoya et al. [35] proposed a fusion method for hyperspectral and multispectral data called Coupled Non-negative Matrix Factorization (CNMF). This approach decomposes both types of data into endmember matrices and abundance matrices through unsupervised unmixing, resulting in a high-resolution hyperspectral image. However, CNMF is an ill-posed inverse problem and requires the addition of regularization terms to enhance its effectiveness. Some researchers have made improvements to CNMF. For instance, Yang et al. [36] introduced endmember vertex distance and iterative centroid proximity regularization terms, incorporating sparse and proximal regularization terms into CNMF. This approach reduces the computational complexity through proximal alternating optimization. Simões et al. [37] proposed the HySure method, which uses subspace regularization and vector total variation to effectively address the fusion of multispectral and panchromatic images. Dian et al. [38] utilized clustering structures to propose a low-dimensional spectral subspace representation model, which effectively captures the self-similarity of hyperspectral images and enhances image reconstruction. Xue et al. [21] developed a subspace clustering method based on structured sparse low-rank representation, where data samples are expressed as linear combinations of a dictionary and a basis matrix. This approach leverages learned spatial and spectral low-rank structures for hyperspectral image super-resolution.

However, in terms of computational efficiency and cost-effectiveness, tensor decomposition methods have shown unique advantages compared to matrix factorization methods. They comprehensively capture both spatial and spectral information in hyperspectral images. This comprehensiveness allows tensor decomposition to excel in preserving image details and features, thereby improving image resolution and information recovery capabilities.

Prévost et al. [39] proposed a new method for hyperspectral super-resolution (Coupled Tucker Approximation: Recoverability and SVD-Based, SCOTT), which is based on a coupled low-rank multicollinear (Tucker) model, using low-rank tensor approximation for hyperspectral super-resolution reconstruction of images. To implement the coupled tensor approximation, they proposed two SVD-based algorithms that are simple and efficient for dealing with unknown spatial degradation and generalized sharpening. By optimizing these algorithms, more accurate super-resolution reconstruction results of hyperspectral images can be obtained. Kanatsoulis et al. [40] proposed a new framework based on coupled tensor factorization (STEREO). This method aims to identify the super-resolution tensor and is suitable for dealing with cases where the spatial degradation model is unclear or inaccurately estimated. By optimizing this framework, the super-resolution reconstruction of hyperspectral images can be improved.

Dian et al. [41] introduced the non-local sparse tensor factorization (NLSTF) method to address hyperspectral image super-resolution. They treated hyperspectral data as a three-way tensor and applied Tucker decomposition techniques while incorporating the non-local similarity prior to hyperspectral images. This approach combines sparse tensor decomposition and non-local means into a unified framework, comprehensively addressing the issues of dictionary estimation and the sparse core tensor of each cube. This innovative method has brought new insights and solutions to the field of hyperspectral image super-resolution. Zeng et al. [25] investigated a multi-modal core tensor factorization (MCTF) method, which combines tensor low-rank modeling with a non-convex relaxation form (NC-MCTF). This model integrates the low-rank information provided by the Tucker and T-SVD methods, enabling simultaneous modeling of low-rank structures in multiple spectral directions and accurately recovering the intrinsic low-rank structural data using a small number of observed entries.

In summary, these methods make full use of tensor decomposition techniques to handle the multi-dimensional information in hyperspectral images, providing powerful tools and approaches for the field of hyperspectral image processing. In recent years, with the in-depth study of high-dimensional tensors, methods based on tensor chains and tensor rings have gradually emerged and gained popularity. These methods offer valuable insights and directions for future research in hyperspectral image processing.

In 2019, Dian [42] introduced a prior model based on the low tensor sequence rank, learning the correlations among spatial, spectral, and non-local modes in an HR-HSI cube with non-local similarity. This method transformed the super-resolution problem into a tensor chain rank regularization optimization problem. Jin et al. [23] employed a fully connected coupled tensor network to decompose the spatial structure of high-order tensors and used graph regularization techniques to preserve spectral information, thus improving the quality of the image.

Zhao et al. [43] first introduced the tensor ring decomposition (TR) model in 2016 based on the basic tensor decomposition model, which represents a high-dimensional tensor as a third-order tensor with cyclic interconnection properties by performing cyclic linear product operations on a low-dimensional core, and at the same time uses a graphical interpretation of it. The research focuses on how the properties of multilinear algebra can be effectively utilized by performing direct operations on TR representations (i.e., cores). This provides a potentially powerful framework for dealing with large-scale data. Additionally, it explores the connections with other tensor decomposition models, allowing for easy migration of potential representations from traditional models to TR models.

The coupled tensor ring model for hyperspectral image super-resolution problem was proposed by He et al. [44]. The model combines the advantages of coupling matrix and Tucker factorization and can better exploit the global spectral low-rank properties of recovered high-resolution hyperspectral images by using the kernel normalization tensor of the third core. The model has been shown to perform better in exploiting the low-rank properties of different hyperspectral image classes. By optimizing the CTRF model, accurate super-resolution reconstruction of hyperspectral images can be achieved, thus improving the spatial resolution and the quality of spectral information of the images.

Chen et al. [34] proposed an efficient FSTRD method to reconstruct HR-HSIs from an LR-HSI and HR-MSI pair of the same scene. The method introduces a smooth regularization term to constrain the factors in the tensor ring decomposition to capture the spatial–spectral continuity of HR-HSIs. Meanwhile, a PAM (Projected Alternating Minimization) algorithm was used to optimize each factor. Through the optimization process, the FSTRD method can effectively reconstruct the high-resolution hyperspectral images (HRIs) and improve the spatial resolution and spectral continuity of the images.

Despite the promising future of tensor ring decomposition methods in the field of hyperspectral image super-resolution, they still face several challenges and problems, which include difficulties in overfitting problems, real dataset validation, and noise interference. This paper presents a new algorithm for super-resolution of hyperspectral images based on graph-regular tensor ring decomposition. The contributions of the proposed method are as follows.

• (1). The graph regularization term is introduced, and the symmetric normalized Laplacian matrix is used to constrain the model, which provides a better structure representation and feature extraction ability for the model. This optimization method can not only effectively improve the structure and feature representation of the image, but also make full use of the spatial and spectral information in the process of image reconstruction, so as to preserve the key details and edge information.

• (2). By introducing a spatial smoothness constraint term, the paper successfully achieves spatial smoothness in the reconstructed images. This constraint not only preserves the textures and structural characteristics of the images but also effectively addresses overfitting and optimization oscillations. As a result, it ensures the stability and quality of the hyperspectral image super-resolution algorithm, providing a safeguard for its performance.

• (3). This method is based on tensor ring decomposition, fully utilizing the dual characteristics of hyperspectral images’ spatial structure and spectral information. By constructing a graph model and applying regularization techniques, it effectively captures both local and global features of the image, while also successfully preserving the spectral information. This enables achieving more comprehensive and accurate image reconstruction in the context of super-resolution.

• (4). The model proposed in this paper employs a multi-stage optimization framework, wherein the results of tensor ring decomposition are iteratively refined concurrently with the integration of regularization constraints. This multi-stage optimization strategy enhances the convergence and stability of the model, thereby yielding more precise and dependable super-resolution outcomes.

The algorithm in this paper has been experimentally demonstrated to achieve a significant improvement in the hyper-resolution task of hyperspectral images. Compared with traditional methods, our proposed algorithm not only produces clearer, detail-rich high-resolution images but also effectively retains the spectral information, particularly evident in real datasets. Therefore, the hyperspectral image super-resolution algorithm based on graph regular tensor ring decomposition has a broad potential and feasibility in practical applications.

The rest of the paper is organized according to the following structure. Section 2 provides the background information of the paper and details the hyperspectral image super-resolution algorithm (GRTR) based on graph regular tensor ring decomposition. Section 3 describes the optimization process of the GRTR algorithm in detail. Section 4 presents and analyzes the experimental results. Experimental evaluations are performed using multiple hyperspectral image datasets and compared with other classical methods. The performance of different algorithms in terms of spatial resolution enhancement and spectral fidelity is analyzed, and the quality of the reconstructed images is quantitatively evaluated. Finally, Section 5 summarizes the advantages and limitations of the GRTR algorithm and discusses possible future directions for improvement.

2. Background

2.1. Symbolic Representation

In this paper, the higher-order tensor is denoted using a calligraphic letter, $\mathcal{X}\in {\mathbb{R}}^{N_1\times N_2\times ...\times N_n}$ denotes a tensor of order n, and $N_n$ denotes the size of the nth dimension. The third-order tensor, matrix, and vector are represented by the cursive letter $\mathcal{X}\in {\mathbb{R}}^{N_1\times N_2\times N_3}$, the upper-case letter $X\in {\mathbb{R}}^{N_1\times N_2}$, and the lower-case letter $x\in {\mathbb{R}}^{N_1}$, respectively.

The high-resolution hyperspectral images is defined as the third-order tensor $\mathcal{X}$, $\mathcal{X}\in {\mathbb{R}}^{N_W\times N_H\times N_S}$, the low-resolution hyperspectral image (HSI) is defined as $\mathcal{Y}\in {\mathbb{R}}^{N_w\times N_h\times N_S}$, and the high-resolution multispectral image (MSI) is defined as $\mathcal{Z}\in {\mathbb{R}}^{N_W\times N_H\times N_s}$, where $N_W$ and $N_w$ denote the magnitude of the broad mode, $0<N_w<N_W$, $N_H$ and $N_h$ denote the magnitude of the high mode, $0<N_h<N_H$, $N_S$ and $N_s$ denote the magnitude of the spectral mode, $0<N_s<N_S$. ${\mathcal{X}}_{\times n}X$ represents the product of the tensor $\mathcal{X}$ and the matrix X in the nth mode. The three-state expansion of the third-order tensor is denoted as

(1)$X_{\left(1\right)}=[X{(:,:,1)}^T,X{(:,:,2)}^T,...,X{(:,:,N_S)}^T]\in {\mathbb{R}}^{N_W\times N_S{N}_H},$

(2)$X_{\left(2\right)}=[X{(:,:,1)}^T,X{(:,:,2)}^T,...,X{(:,:,N_S)}^T]\in {\mathbb{R}}^{N_H\times N_W{N}_S},$

(3)$X_{\left(3\right)}=[X{(:,1,:)}^T,X{(:,2,:)}^T,...,X{(:,N_H,:)}^T]\in {\mathbb{R}}^{N_S\times N_W{N}_H}.$

2.2. Tensor Decomposition

A tensor represents an abstraction of a multilinear mapping. It is an extension of a multidimensional array or multidimensional vector that can be used to represent a linear relationship between multiple vectors, matrices, or higher dimensional arrays. A tensor can be thought of as a multidimensional object that can have any number of dimensions, each of which can have an arbitrary size. The dimensions describe how the elements in the tensor are arranged, while the size of the tensor indicates the number of elements in each dimension. Common tensor decomposition methods include CP decomposition, Tucker decomposition, Tensor Train decomposition, and tensor ring decomposition.

CP decomposition [45] is a basic tensor decomposition method that decomposes a high-dimensional tensor into the product of multiple low-dimensional tensors. In hyperspectral image fusion, CP decomposition can decompose a hyperspectral image into three tensors that represent the relationship of spectra, space, and data, respectively, and each sub-tensor has rank 1 in a single dimension. The mathematical representation of CP decomposition is

(4)$\mathcal{X}=\sum _{r=1}^R{a}_r\circ b_r\circ c_r.$

where R is the magnitude of tensor rank, and the three sub-tensors are multiplied in the form of the outer product. $A=\left[a_1,a_2,...,a_r\right]$, $B=\left[b_1,b_2,...,b_r\right]$, and $C=\left[c_1,c_2,...,c_r\right]$ denote low-order potential factors, so CP decomposition can also be denoted as

(5)$\mathcal{X}=\left[A,B,C\right].$

Although CP decomposition is concisely expressed, determining its rank is an NP-hard problem, making it challenging to ascertain the magnitude of the tensor rank.

Tucker decomposition [46], also known as higher-order singular value decomposition (HOSVD), is a widely used higher-order tensor decomposition method, which decomposes a tensor into the product of a kernel tensor and multiple factor tensors. In hyperspectral image fusion, Tucker decomposition can decompose a hyperspectral image into a core tensor and multiple factor tensors. Each factor tensor represents different features, such as spectral, spatial, etc. The expression of Tucker decomposition is

(6)$\mathcal{X}=\mathcal{G}{\times }_1{A}^{\left(1\right)}{\times }_2{A}^{\left(2\right)}{\times }_3{A}^{\left(3\right)}.$

where $\mathcal{G}\in {\mathbb{R}}^{N_w\times N_h\times N_s}$ represents the core tensor, $A^{\left(1\right)}\in {\mathbb{R}}^{N_W\times N_w}$ and $A^{\left(2\right)}\in {\mathbb{R}}^{N_H\times N_h}$ represent the spatial feature tensor, and $A^{\left(3\right)}\in {\mathbb{R}}^{N_S\times N_s}$ represents the spectral feature tensor.

Tensor Train decomposition [47] processes higher-order tensors by decomposing them into a product of multiple small kernel tensors. In hyperspectral image fusion, Tensor Train decomposition can decompose a hyperspectral image into multiple small kernel tensors. Each kernel tensor represents features like spectra, space, and so on. The expression of Tensor Train decomposition is

(7)$\mathcal{X}(i_1,i_2,i_3)=G^{\left(1\right)}(:,i_1,:)G^{\left(2\right)}(:,i_2,:)G^{\left(3\right)}(:,i_3,:).$

where $G^{\left(1\right)}\in {\mathbb{R}}^{R_1\times N_W\times R_2}$, $G^{\left(2\right)}\in {\mathbb{R}}^{R_2\times N_H\times R_3}$, $G^{\left(3\right)}\in {\mathbb{R}}^{R_3\times N_S\times R_1}$. The rank of Tensor Train is defined as $ran{k}_{TT}\left(\mathcal{X}\right)=(R_1,R_2,R_3)$, $R_1=R_3=1$.

2.3. Tensor Ring Decomposition

Tensor ring decomposition is a method that represents higher-order tensors as a ring product of a series of lower-order tensors. It is an extension and improvement of traditional tensor decomposition methods such as CP decomposition and Tucker decomposition. In tensor ring decomposition, the original high-dimensional tensor is represented as a ring structure consisting of multiple low-dimensional tensors. Each low-dimensional tensor represents the projection of the original tensor onto a particular mode. By interconnecting these low-dimensional tensors, the original tensor can be recovered. The expression of tensor ring decomposition is

(8)$\mathcal{X}(i_1,i_2,i_3)=tr\left(G^{\left(1\right)}(:,i_1,:)G^{\left(2\right)}(:,i_2,:)G^{\left(3\right)}(:,i_3,:)\right).$

where $tr(·)$ denotes the matrix trace multiplication. The rank of the tensor ring is defined as $ran{k}_{TR}\left(\mathcal{X}\right)=(R_1,R_2,R_3)$, $R_1=R_3$. $G^{\left(n\right)}$ is the nth transverse slice matrix of the core tensor ${\mathcal{G}}^{\left(n\right)}$, so the HRI under the tensor ring decomposition can also be expressed as

(9)$\mathcal{X}=\Phi ({\mathcal{G}}^{\left(1\right)},{\mathcal{G}}^{\left(2\right)},{\mathcal{G}}^{\left(3\right)}).$

where ${\mathcal{G}}^{\left(1\right)}\in {}^{R_1\times N_W\times R_2}$, ${\mathcal{G}}^{\left(2\right)}\in {}^{R_2\times N_H\times R_3}$, and ${\mathcal{G}}^{\left(3\right)}\in {}^{R_2\times N_S\times R_3}$. The matrix representation of the TR can be expressed as $X_{<n>}=G_{{}^{\left(2\right)}}^n{\left(G_{<2>}^{\ne n}\right)}^T$. The main advantage of tensor ring decomposition is its ability to handle high-dimensional data efficiently and keep the storage and computational complexity low. It also has better local representation performance to capture the local structure and dependencies in high-dimensional data so that it can efficiently process high-dimensional data and extract useful information.

2.4. Graph Regularization

Graph regularization [48] is a graph-theory-based regularization method for processing data with a graph structure, such as social networks, images, videos, speech, etc. Its main aim is to improve the performance of a model by introducing structural information from the graph into the model. It is based on the concept of graph theory by establishing connection relations between data samples and using these connection relations to enhance the performance and robustness of the learning algorithm.

For a graph, $G=(V,E)$, where V is the set of vertices and E is the set of edges, let W be a symmetric weight matrix whose elements $W(i,j)$ denote the weights of the edges between vertices i and j. Define the degree matrix D to be a diagonal matrix whose element $D(i,i)$ denotes the degree of vertex i (the sum of the weights of the edges connected to it). The symmetric normalized Laplace matrix $L_{sym}$ is defined as:

(10)$L_{sym}=D^{-1/2}(D-W)D^{-1/2}=I-D^{-1/2}W{D}^{-1/2}.$

I is the unit matrix, and $D^{1/2}$ is the matrix obtained by calculating the positive square root of all elements of D. The symmetric normalized Laplace matrix $L_{sym}$ has the following characteristics: Symmetry—$L_{sym}$ is a symmetric matrix, i.e., $L_{sym}(i,j)=L_{sym}(j,i)$. Normalization—$L_{sym}$ is normalized by the normalization of the degree matrix, such that the weights corresponding to the degrees of each vertex are normalized.

In this paper, the spatial map of $\mathcal{X}$ is denoted as $G_{SPA}=(V_{SPA},E_{SPA})$, and the normalized spatial Laplace matrix with respect to the graph $G_{SPA}$ can be denoted as

(11)$L_{SPA}=D_{SPA}^{-1/2}(D_{SPA}-W_{SPA})D_{SPA}^{-1/2}=I_{SPA}-D_{SPA}^{-1/2}{W}_{SPA}{D}_{SPA}^{-1/2}.$

The spectral map of $\mathcal{X}$ is denoted as $G_{SPE}=(V_{SPE},E_{SPE})$, and the normalized spatial Laplace matrix with respect to the graph can be denoted as

(12)$L_{SPE}=D_{SPE}^{-1/2}(D_{SPE}-W_{SPE})D_{SPE}^{-1/2}=I_{SPE}-D_{SPE}^{-1/2}{W}_{SPE}{D}_{SPE}^{-1/2}.$

where $I_{SPA}$ and $I_{SPE}$ are unitary matrices with matching sizes.

$L_{SPA}$ focuses on the similarity between pixels in a local region, which can be defined based on the distance between pixels or other similarity metrics, and thus it is able to capture spatial features such as textures, details, and edges in an image. By constraining $L_{SPA}$, the super-resolution algorithm is able to better preserve and reconstruct the spatial features of the image, resulting in a clearer and more detailed synthetic high-resolution image. The $L_{SPE}$ can capture the spectral features in the image, i.e., the correlation and consistency between different bands. By constraining $L_{SPE}$, the super-resolution algorithm ensures that the synthesized high-resolution image is consistent with the original image in terms of spectral features, thus preserving the color accuracy and spectral detail of the image. Collectively, the spatial constraint term $L_{SPA}$ plays a role in preserving and reconstructing spatial attributes like texture, details, and edges within an image. Similarly, the spectral constraint term $L_{SPE}$ contributes to upholding the spectral precision and consistency of an image. The amalgamation of these two constraint terms synergistically enhances the reconstruction prowess of the hyperspectral image super-resolution algorithm, resulting in synthesized high-resolution images that are more faithful and accurate in representation.

Graph regularization involves constructing a graph for the hyperspectral image, which models the relationships between adjacent pixels. These relationships can encompass spatial positions, spectral similarities, and more. By introducing regularization terms defined on the graph, the output of the super-resolution algorithm can be constrained, promoting smoother and more coherent high-resolution images in both spatial and spectral aspects. Additionally, graph regularization helps maintain the image’s structural and textural features, preventing the generation of overly smooth or distorted images.

The objective of hyperspectral image super-resolution is to recover a high-resolution image from a low-resolution hyperspectral image, aiming to enhance image details and spatial resolution. Graph regularization, on the other hand, harnesses the spatial and spectral characteristics within hyperspectral images and leverages the connections within the constructed graph. By doing so, it introduces additional constraints and prior information to enhance the performance of super-resolution algorithms.

2.5. Spatial Smoothing Constraints

Spatial smoothing constraint is a regularisation constraint method commonly used in image processing and computer vision tasks. It achieves the smoothing of an image by constraining the neighboring pixels of the image and limiting the gradient change of the image. For an image $\mathcal{G}$, the spatial smoothing constraint on it is denoted as

(13)${\Psi }_d\left(\mathcal{G}\right)=\sqrt{\sum _{j=1}^{N_H}\sum _{\{m,n\}\in \Omega }{\left(G_j\left(m\right)-G_j\left(n\right)\right)}^2}$

In this context, $\Omega$ represents the set of all pairs of neighboring elements in the graph $\mathcal{G}$. The matrix G is derived by unfolding $\mathcal{G}$ in accordance with the spectral mode, where $G_j\left(m\right)$ signifies the mth element of the jth column. A diminished value of the smoothness constraint term suggests a decrease in differences between adjacent pixels in the image, producing smoother transitions and preserving the innate textures and structural attributes. Moreover, this results in a reduction in the potential noise and irregularities within the image. By incorporating the spatial smoothing constraint, a standard of smoothness is imposed on the output. This constraint compels the model to generate results that are more uniform and streamlined. Such a method not only diminishes the likelihood of overfitting due to noise or localized details in the training dataset but also stabilizes the optimization process. This stability enhances the possibility of the model arriving at a globally optimal solution or a more desirable local optimum, ultimately refining the generalization capabilities and optimization outcomes.

2.6. Proposed Method

In the context of hyperspectral image super-resolution, the relationship between HSRs and HSIs can be understood as arising from spatial downsampling. Specifically, a HRI undergoes downsampling operations along both spatial dimensions, resulting in the creation of an HSI. Thus, we can express the generation relationship of HSIs as follows:

(14)$\mathcal{Y}=\Phi ({\mathcal{G}}^{\left(1\right)}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right)}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right)})$

where $\mathcal{Y}$ represents the HSI, and $U_1\in {\mathbb{R}}^{N_w\times N_W}$ and $U_2\in {\mathbb{R}}^{N_H\times N_h}$ are spatial downsampling matrices. Similarly, MSIs can be seen as obtained through spectral downsampling.

Thus, MSIs can be represented as follows:

(15)$\mathcal{Z}=\Phi ({\mathcal{G}}^{\left(1\right)},{\mathcal{G}}^{\left(2\right)},{\mathcal{G}}^{\left(3\right)}{\times }_2{U}_3)$

where $\mathcal{Z}$ represents the MSI and $U_3\in {\mathbb{R}}^{N_S\times N_s}$ represents the spectral downsampling matrix.

The hyperspectral image super-resolution algorithm based on graph-regularized tensor ring decomposition integrates the concepts of graph regularization and tensor ring decomposition when addressing hyperspectral image super-resolution challenges. In this approach, graph regularization is employed to capture the structural information within the image, while tensor ring decomposition aids in extracting latent spectral and spatial features. Within the algorithm’s framework, spatial smoothness constraints are applied to spatial sub-tensors, further enhancing the reconstruction effectiveness of the images. The specific formulation of this algorithm is as follows:

(16)$\begin{array}{c}\underset{{\mathcal{G}}^{\left(1\right)},{\mathcal{G}}^{\left(2\right)},{\mathcal{G}}^{\left(3\right)}}{min}\frac{1}{2}{∥\mathcal{Y}-\Phi ({\mathcal{G}}^{\left(1\right)}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right)}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right)})∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-\Phi ({\mathcal{G}}^{\left(1\right)},{\mathcal{G}}^{\left(2\right)},{\mathcal{G}}^{\left(3\right)}{\times }_2{U}_3)∥}_F^2\hfill \\ +\alpha tr({{\mathcal{G}}_{\left(2\right)}^{\left(3\right)}}^T{L}_{SPE}{\mathcal{G}}_{\left(2\right)}^{\left(3\right)})+\beta tr(({\mathcal{G}}_{\left(2\right)}^{\left(2\right)}\otimes {\mathcal{G}}_{\left(2\right)}^{\left(1\right)}{)}^T{L}_{SPA}({\mathcal{G}}_{\left(2\right)}^{\left(2\right)}\otimes {\mathcal{G}}_{\left(2\right)}^{\left(1\right)}))+{\Psi }_d({\mathcal{G}}^{\left(1\right)})+{\Psi }_d({\mathcal{G}}^{\left(2\right)})\hfill \end{array}$

where ${∥·∥}_F^{}$ denotes the Frobenius norm and $\lambda$, $\alpha$, and $\beta$ are the regularization parameters. ⊗ denotes the Kronecker product. $L_{SPE}$ and $L_{SPA}$ denote the symmetric normalized Laplacian matrices, and ${\Psi }_d\left({\mathcal{G}}^{\left(1\right)}\right)$ and ${\Psi }_d\left({\mathcal{G}}^{\left(2\right)}\right)$ represent spatial smoothness constraints. Through a multi-stage optimization framework, the algorithm progressively refines the results of ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{G}}^{\left(2\right)}$, and ${\mathcal{G}}^{\left(3\right)}$ to achieve hyperspectral image super-resolution reconstruction. Within this algorithm, the combination of graph regularization and tensor ring decomposition, along with the introduction of spatial smoothness constraints, leads to better integration of the structure and features of hyperspectral images. Such comprehensive constraints not only enhance the quality of image reconstruction but also aid in preserving image details and features, resulting in outstanding performance in hyperspectral image super-resolution tasks.

Figure 1 illustrates the overall framework of this algorithm, encompassing key steps such as data preprocessing, graph construction, tensor ring decomposition, objective function formulation, and super-resolution reconstruction. The process begins with HRI and undergoes data preprocessing, including spatial and spectral downsampling operations, resulting in corresponding HSI and MSI outputs. Subsequently, graph construction is a pivotal step within the algorithm. In the process of graph construction, pixels in hyperspectral and multispectral images are considered as nodes of the graph, and the normalized Laplace matrix of the graph is built based on the similarity and distance relationship between pixels. This approach effectively integrates spatial and spectral information, providing a foundation for subsequent regularization and optimization processes. By introducing the tensor ring decomposition technique, the hyperspectral and multispectral images are decomposed to obtain latent feature representations ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{G}}^{\left(2\right)}$, and ${\mathcal{G}}^{\left(3\right)}$. These features have the capacity to capture underlying patterns and information within the data. To achieve the super-resolution reconstruction of hyperspectral images, we devised an objective function that integrates the features obtained from tensor ring decomposition with the relationships between the input hyperspectral image and the multispectral image. This objective function encompasses data terms, regularization terms, and spatial smoothness constraint terms, collectively aiming to enhance reconstruction quality while preserving image features.

The entire algorithm employs a multi-stage iterative optimization process to progressively refine ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{G}}^{\left(2\right)}$, and ${\mathcal{G}}^{\left(3\right)}$. Each stage of optimization thoroughly considers the constraints of graph regularization and tensor ring decomposition, aiming to iteratively improve the feature representations. Ultimately, based on the optimized features ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{G}}^{\left(2\right)}$, and ${\mathcal{G}}^{\left(3\right)}$, hyperspectral image super-resolution reconstruction is performed. Through reverse tensor ring decomposition and graph regularization, the algorithm recovers detailed information about the high-resolution image, achieving hyperspectral image super-resolution reconstruction.

3. Optimization

In this section, an efficient optimization algorithm is used to solve the above model. To achieve the optimization of the model, we continuously improve the performance of the algorithm by iteratively updating ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{G}}^{\left(2\right)}$, and ${\mathcal{G}}^{\left(3\right)}$. Specifically, the iterative update formulas for ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{G}}^{\left(2\right)}$, and ${\mathcal{G}}^{\left(3\right)}$ are shown below.

(17)$\left\{\begin{array}{c}{\mathcal{G}}^{\left(1\right),i+1}=\underset{{\mathcal{G}}^{\left(1\right)}}{argmin}\Gamma ({\mathcal{G}}^{\left(1\right)},{\mathcal{G}}^{\left(2\right),i},{\mathcal{G}}^{\left(3\right),i})+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(1\right)}-{\mathcal{G}}^{\left(1\right),i}∥}_F^2\hfill \\ {\mathcal{G}}^{\left(2\right),i+1}=\underset{{\mathcal{G}}^{\left(2\right)}}{argmin}\Gamma ({\mathcal{G}}^{\left(1\right),i+1},{\mathcal{G}}^{\left(2\right)},{\mathcal{G}}^{\left(3\right),i})+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(2\right)}-{{\mathcal{G}}^{\left(2\right)}}^{,i}∥}_F^2\hfill \\ {\mathcal{G}}^{\left(3\right),i+1}=\underset{{\mathcal{G}}^{\left(3\right)}}{argmin}\Gamma ({\mathcal{G}}^{\left(1\right),i+1},{\mathcal{G}}^{\left(2\right),i+1},{\mathcal{G}}^{\left(3\right)})+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(3\right)}-{{\mathcal{G}}^{\left(3\right)}}^{,i}∥}_F^2\hfill \end{array}\right.$

where the function $\Gamma \left(·\right)$ is implicitly defined, $\rho$ is a positive parameter, and ${\mathcal{G}}^{\left(n\right),i}$ denotes the value in the previous iteration. These iterative update formulas allow us to continuously adjust the values of ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{G}}^{\left(2\right)}$, and ${\mathcal{G}}^{\left(3\right)}$ at each iteration step. By iteratively updating, we can gradually approach the optimal solution and thus optimize the model.

The flow chart of the algorithm on the GRTR model is obtained by iteratively updating the formula, as shown in Algorithm 1.

Algorithm 1 GRTR for Hyperspectral Image Super-Resolution

Initialize spatial sub-tensor ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{G}}^{\left(2\right)}$, and spectral sub-tensor ${\mathcal{G}}^{\left(3\right)}$;  $\mathbf{While}$ $\mathbf{not}$ $\mathbf{do}$        Update sub-tensor ${\mathcal{G}}^{\left(1\right)}$ via solving Algorithm A1, Appendix A;        ${\mathcal{G}}^{\left(1\right),i+1}={\mathcal{G}}^{\left(1\right)}$;        Update sub-tensor ${\mathcal{G}}^{\left(2\right)}$ via solving Algorithm A2;        ${\mathcal{G}}^{\left(2\right),i+1}={\mathcal{G}}^{\left(2\right)}$;        Update sub-tensor ${\mathcal{G}}^{\left(3\right)}$ via solving Algorithm A3;        ${\mathcal{G}}^{\left(3\right),i+1}={\mathcal{G}}^{\left(3\right)}$;  $\mathbf{end}$ $\mathbf{While}$;  Estimating target HRI $\mathcal{X}$ via Equation (9)

3.1. Optimization Process for ${\mathcal{G}}^{\left(1\right)}$

After maintaining ${\mathcal{G}}^{\left(2\right)}$ and ${\mathcal{G}}^{\left(3\right)}$ constant, the optimization problem concerning ${\mathcal{G}}^{\left(1\right)}$ can be formulated as follows:

(18)$\begin{array}{c}\underset{{\mathcal{G}}^{\left(1\right)}}{min}\frac{1}{2}{∥\mathcal{Y}-\Phi ({\mathcal{G}}^{\left(1\right)}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right),i}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right),i})∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-\Phi ({\mathcal{G}}^{\left(1\right)},{\mathcal{G}}^{\left(2\right),i},{\mathcal{G}}^{\left(3\right),i}{\times }_2{U}_3)∥}_F^2\hfill \\ +\frac{\beta }{2}tr(({\mathcal{G}}_{\left(2\right)}^{\left(2\right),i}\otimes {\mathcal{G}}_{\left(2\right)}^{\left(1\right)}{)}^T{L}_{SPA}({\mathcal{G}}_{\left(2\right)}^{\left(2\right),i}\otimes {\mathcal{G}}_{\left(2\right)}^{\left(1\right)}))+{\Psi }_d({\mathcal{G}}^{\left(1\right)})+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(1\right)}-{{\mathcal{G}}^{\left(1\right)}}^{,i}∥}_F^2\hfill \end{array}$

The optimization problem regarding ${\mathcal{G}}^{\left(1\right)}$ can be solved using the Alternating Direction Method of Multipliers (ADMM) [49]. For computational convenience, let ${\Psi }_d\left({\mathcal{G}}^{\left(1\right)}\right)$ be represented in matrix form as ${\Psi }_d\left({\mathcal{G}}^{\left(1\right)}\right)={\parallel D_1{\mathcal{G}}^{\left(1\right)}\parallel }_2$, where $D_1$ is the spatial smoothness constraint matrix. By introducing the variable ${\mathcal{R}}_1={\mathcal{G}}^{\left(1\right)}$, the constrained problem is transformed into the optimization of the following augmented Lagrangian function:

(19)$\begin{array}{c}{L}_{{\mu }_1}({\mathcal{G}}^{\left(1\right)},{\mathcal{R}}_1,{\mathcal{O}}_1)=\hfill \\ \frac{1}{2}{∥\mathcal{Y}-\Phi ({\mathcal{G}}^{\left(1\right)}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right),i}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right),i})∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-\Phi ({\mathcal{G}}^{\left(1\right)},{\mathcal{G}}^{\left(2\right),i},{\mathcal{G}}^{\left(3\right),i}{\times }_2{U}_3)∥}_F^2\hfill \\ +\frac{\beta }{2}tr(({\mathcal{G}}^{\left(2\right),i}\otimes {\mathcal{R}}_1{)}^T{L}_{SPA}({\mathcal{G}}^{\left(2\right),i}\otimes {\mathcal{R}}_1))+{∥D_1{\mathcal{R}}_1∥}_2^2\hfill \\ +\frac{\rho }{2}{∥{\mathcal{G}}^{\left(1\right)}-{\mathcal{G}}^{\left(1\right),i}∥}_F^2+\frac{{\mu }_1}{2}{∥{\mathcal{G}}^{\left(1\right)}-{\mathcal{R}}_1-{\mathcal{O}}_1∥}_F^2\hfill \end{array}$

where ${\mu }_1$ represents the penalty parameter with ${\mu }_1>0$, and ${\mathcal{O}}_1$ denotes the Lagrange multiplier. By taking the derivative of the Lagrangian function and setting it equal to zero, the update rule for ${\mathcal{G}}^{\left(1\right)}$ can be obtained.

3.2. Optimization Process for ${\mathcal{G}}^{\left(2\right)}$

After fixing ${\mathcal{G}}^{\left(1\right)}$ and ${\mathcal{G}}^{\left(3\right)}$, the optimization problem for ${\mathcal{G}}^{\left(2\right)}$ can be expressed as

(20)$\begin{array}{c}\underset{{\mathcal{G}}^{\left(2\right)}}{min}\frac{1}{2}{∥\mathcal{Y}-\Phi ({\mathcal{G}}^{\left(1\right),i+1}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right)}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right),i})∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-\Phi ({\mathcal{G}}^{\left(1\right),i+1},{\mathcal{G}}^{\left(2\right)},{\mathcal{G}}^{\left(3\right),i}{\times }_2{U}_3)∥}_F^2\hfill \\ +\frac{\beta }{2}tr(({\mathcal{G}}_{\left(2\right)}^{\left(2\right),i}\otimes {\mathcal{G}}_{\left(2\right)}^{\left(1\right)}{)}^T{L}_{SPA}({\mathcal{G}}_{\left(2\right)}^{\left(2\right),i}\otimes {\mathcal{G}}_{\left(2\right)}^{\left(1\right)}))+{\Psi }_d({\mathcal{G}}^{\left(2\right)})+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(2\right)}-{{\mathcal{G}}^{\left(2\right)}}^{,i}∥}_F^2\hfill \end{array}$

By introducing the variable ${\mathcal{R}}_2={\mathcal{G}}^{\left(2\right)}$, the constrained problem is transformed into the optimization of the following augmented Lagrangian function:

(21)$\begin{array}{c}{L}_{{\mu }_2}({\mathcal{G}}^{\left(1\right)},{\mathcal{R}}_1,{\mathcal{O}}_1)=\hfill \\ \frac{1}{2}{∥\mathcal{Y}-\Phi ({\mathcal{G}}^{\left(1\right),i+1}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right)}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right),i})∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-\Phi ({\mathcal{G}}^{\left(1\right),i+1},{\mathcal{G}}^{\left(2\right)},{\mathcal{G}}^{\left(3\right),i}{\times }_2{U}_3)∥}_F^2\hfill \\ +\frac{\beta }{2}tr(({\mathcal{R}}_2\otimes {\mathcal{G}}^{\left(1\right),i+1}{)}^T{L}_{SPA}({\mathcal{R}}_2\otimes {\mathcal{G}}^{\left(1\right),i+1}))+{∥D_2{\mathcal{R}}_2∥}_2^2\hfill \\ +\frac{\rho }{2}{∥{\mathcal{G}}^{\left(2\right)}-{\mathcal{G}}^{\left(2\right),i}∥}_F^2+\frac{{\mu }_2}{2}{∥{\mathcal{G}}^{\left(2\right)}-{\mathcal{R}}_2-{\mathcal{O}}_2∥}_F^2\hfill \end{array}$

where ${\mu }_2$ represents the penalty parameter with ${\mu }_2>0$, and ${\mathcal{O}}_2$ denotes the Lagrange multiplier. $D_2$ is the spatial smoothness constraint matrix. By taking the derivative of the Lagrangian function and setting it equal to zero, the update rule for ${\mathcal{G}}^{\left(2\right)}$ can be obtained.

3.3. Optimization Process for ${\mathcal{G}}^{\left(3\right)}$

With ${\mathcal{G}}^{\left(1\right)}$ and ${\mathcal{G}}^{\left(2\right)}$ kept constant, the optimization problem concerning ${\mathcal{G}}^{\left(3\right)}$ can be represented as follows:

(22)$\begin{array}{c}\underset{{\mathcal{G}}^{\left(3\right)}}{min}\frac{1}{2}{∥\mathcal{Y}-\Phi ({\mathcal{G}}^{\left(1\right),i+1}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right),i+1}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right)})∥}_F^2\hfill \\ +\frac{\lambda }{2}{∥\mathcal{Z}-\Phi ({\mathcal{G}}^{\left(1\right),i+1},{\mathcal{G}}^{\left(2\right),i+1},{\mathcal{G}}^{\left(3\right)}{\times }_2{U}_3)∥}_F^2\hfill \\ +\frac{\alpha }{2}tr(\mathcal{G}{}_{\left(2\right)}^{\left(3\right)T}{L}_{SPE}{\mathcal{G}}_{\left(2\right)}^{\left(3\right)})+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(3\right)}-{{\mathcal{G}}^{\left(3\right)}}^{,i}∥}_F^2\hfill \end{array}$

By introducing the Lagrange multiplier ${\mathcal{R}}_3={\mathcal{G}}^{\left(3\right)}$, the constrained problem is transformed into the optimization of the following augmented Lagrangian function:

(23)$\begin{array}{c}{L}_{{\mu }_3}({\mathcal{G}}^{\left(3\right)},{\mathcal{R}}_3,{\mathcal{O}}_3)=\hfill \\ \frac{1}{2}{∥\mathcal{Y}-\Phi ({\mathcal{G}}^{\left(1\right),i+1}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right),i+1}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right)})∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-\Phi ({\mathcal{G}}^{\left(1\right),i+1},{\mathcal{G}}^{\left(2\right),i+1},{\mathcal{G}}^{\left(3\right)}{\times }_2{U}_3)∥}_F^2\hfill \\ +\frac{\alpha }{2}tr({{\mathcal{R}}_3}^T{L}_{SPE}{\mathcal{R}}_3)+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(3\right)}-{{\mathcal{G}}^{\left(3\right)}}^{,i}∥}_F^2+\frac{{\mu }_3}{2}{∥{\mathcal{G}}^{\left(3\right)}-{\mathcal{R}}_3-{\mathcal{O}}_3∥}_F^2\hfill \end{array}$

where ${\mu }_3$ represents the penalty parameter with ${\mu }_3>0$, and ${\mathcal{O}}_3$ denotes the Lagrange multiplier. By taking the derivative of the Lagrangian function and setting it equal to zero, the update rule for ${\mathcal{G}}^{\left(3\right)}$ can be obtained.

4. Experimental Results and Discussion

4.1. Datasets

In this section, four datasets are used to test the results of the experiments, including two simulated datasets and two real datasets.

4.1.1. Simulated Datasets

The Pavia University dataset is a selection of scenes acquired by the German airborne Reflection Optical Spectroscopic Imager (ROSIS) sensor during a flight campaign in Pavia, Italy, in 2003. The image size is 610 × 340, with 115 bands of continuous imaging in the wavelength range of 0.43–0.86 $\mu$m and a spatial resolution of 1.3 m. Ninety-three bands were selected for the experiment, and a spectral image of size 256 × 256 × 93 was retained as a reference image. Simulated hyperspectral images were obtained by spatial downsampling, and simulated multispectral images were obtained by spectral downsampling (https://paperswithcode.com/dataset/pavia-university (accessed on 1 March 2023)).

The Washington dataset was acquired by the Hyperspectral Digital Imaging Collection Experience (HYDICE) sensor in Washington. The data size is 1208 × 307 and it contains 191 bands in the 0.4–2.4 $\mu$m range. A spectral image of size 256 × 256 × 93 was selected as the reference image for the experiment. Furthermore, the hyperspectral and multispectral images were obtained by spatial downsampling and spectral downsampling, respectively, (https://engineering.purdue.edu/ biehl/MultiSpec/hyperspectral.html (accessed on 1 March 2023)).

4.1.2. Real Datasets

A multispectral dataset of Ningxia’s Sand Lake was acquired by the GF1 WFV sensor during a flight over Sand Lake in Shizuishan City, Ningxia, on 31 July 2019 [50]. The size of the acquired image is 12,000 × 13,400, the spatial resolution of the resulting image is 16 m, and the image has four spectral bands. The hyperspectral dataset was captured by the GF5-AHSI sensor during a flight over Sand Lake in Shizuishan City, Ningxia, on 27 July 2019 (Beijing time). The size of the acquired image is 2774 × 2554 with a spatial resolution of 30 m. There are 330 spectral bands in the image. After discarding unnecessary bands and selecting areas rich in feature information, we obtained a hyperspectral image of size 256 × 256 × 93 and a multispectral image of size 256 × 256 × 3. These were then used for practical data fusion experiments to validate the fusion algorithm’s feasibility.

The Qingtongxia Rice dataset was obtained through UAV (Unmanned Aerial Vehicle) flights equipped with the Pika-L hyperspectral sensor and the Micasense Altum multispectral sensor over experimental rice fields in Qingtongxia City, Ningxia [50]. The UAV operated at a speed of 3 m/s at an altitude of 120 m with a pitch angle of −90 degrees. The hyperspectral sensor captured images with dimensions of 5376 × 1492 which consisted of 150 bands, with a spatial resolution of 10 cm. The multispectral sensor acquired images with dimensions of 16103 × 9677 × 3 and a spatial resolution of 5.2 cm. During data processing, bands containing water vapor and those with significant noise were removed. This processing yielded hyperspectral images of size 510 × 510 with 60 bands and multispectral images of size 510 × 510 with 3 bands, underscoring the fusion algorithm’s viability.

4.2. Comparison Methodology

Our method was compared with several classical methods as well as currently popular methods, including algorithms based on HySure [37], CNMF [35], CTRF [44], SCOTT [39] hyperspectral super-resolution via coupled tensor factorization: identifiability and algorithms (STEREO) [40], and NLSTF [41]. This study used two simulated datasets and two real datasets in the same experimental environment and several experimental comparisons were performed. The algorithm proposed in this paper was compared and evaluated with the above methods in order to analyze its performance and superiority. These comparisons allow us to evaluate the effectiveness of this paper’s method in hyper-resolution tasks for hyperspectral images and to verify its advantages in improving spatial resolution and preserving spectral information.

4.3. Quantitative Indicators

In this paper, eight metrics were used to evaluate the quality of high-resolution hyperspectral images, namely, peak signal-to-noise ratio (PSNR), root mean square error (RMSE), spectral angle mapping (SAM), degree of distortion (DD), universal image quality index (UIQI), error relative global dimensionless synthesis (ERGAS), structural similarity (SSIM), and correlation coefficient (CC).

The PSNR is an objective metric used to measure the level of distortion or noise in an image, and its definition is based on the calculation of the mean square error (MSE). A larger PSNR value indicates that the difference between the fused image and the original image is smaller, suggesting that more detail has been preserved.

(24)$PSNR=10·{log}_{10}\left(\frac{MA{X}_i^2}{MSE}\right)$

(25)$MSE=\frac{1}{N_W{N}_H}\sum _{i=0}^{N_W-1}\sum _{j=0}^{N_H-1}{\left[Y\right(i,j)-X(i,j\left)\right]}^2$

where Y denotes the fused image and X denotes the reference image. $MA{X}_i$ indicates the maximum pixel value on the i-th spectral band.

ERGAS is a measure of spectral variability that reflects the global quality of the reconstructed image. A smaller ERGAS value means better fusion.

(26)$ERGAS(X,Y)=\frac{100}{S}\sqrt{\frac{1}{N_S}{\sum }_{i=1}^{N_S}\frac{MSE(X,Y)}{MEAN\left(Y\right)}}.$

where S denotes the spatial downsampling factor, $N_S$ denotes the number of bands, and $MEAN\left(Y\right)$ denotes the mean value of pixels in each band of the fused image.

SAM is spectral angle mapping, and is used to measure the similarity between the fused image and the spectral profile of the reference image. A smaller output value indicates that the two arrays are more matched and similar.

(27)$SAM(X,Y)=\frac{1}{N_W{N}_H}\sum _{i=1}^{N_W{N}_H}arccos\frac{(X,Y)}{{∥X_i∥}_2{∥Y_i∥}_2}$

The RMSE is a measure of the root mean square error between the image fusion result and each band of the reference image. The smaller the value obtained from the RMSE calculation, the better it is, indicating that the difference between the fusion result and the reference image is smaller and more details of the original two images are preserved.

(28)$RMSE(X,Y)=\sqrt{\frac{∥X,Y∥}{N_W{N}_H{N}_S}}.$

The CC denotes the spectral similarity between the fusion result of the image and the original image, and when the CC is approximately close to 1, this indicates a better spectral quality of the fused image.

(29)$CC(X,Y)=\frac{{\sum }_{i=1}^{N_W}{\sum }_{j=1}^{N_H}{[X(i,j)-V_X]}^2·{[Y(i,j)-V_Y]}^2}{\sqrt{{\sum }_{i=1}^{N_W}{\sum }_{j=1}^{N_H}{[X(i,j)-V_X]}^2·{\sum }_{i=1}^{N_W}{\sum }_{j=1}^{N_H}{[Y(i,j)-V_Y]}^2}}.$

where $V_X$ represents the pixel average in the image and $V_Y$ is the pixel average of the fused image Y.

The DD can indicate the size of the difference between the fusion result and the reference image in terms of spectral information, and a smaller value indicates a better effect of the fused image.

(30)$DD(X,Y)=\frac{1}{N_W{N}_H{N}_S}{∥vec\left(X\right)-vec\left(Y\right)∥}_1.$

SSIM is a measure of the similarity of the two images. SSIM takes values in [0,1]; the larger it is, the closer the image is to the original two images and more details of the original two images are preserved.

(31)$SSIM\left(X,Y\right)=\frac{\left(2{\mu }_X{\mu }_Y+C_1\right)\left({\sigma }_{XY}+C_2\right)}{\left({\mu }_X^2+{\mu }_Y^2+C_1\right)\left({\sigma }_X^2+{\sigma }_Y^2+C_2\right)}$

where ${{\sigma }_X}_Y$ and ${\mu }_Y$ represent the mean values of X and Y, respectively. ${\sigma }_X$ and ${\sigma }_Y$ represent the standard deviation of X and Y. ${{\sigma }_X}_Y$ represents the covariance of X and Y.

The UIQI is an index to detect the average correlation between the fused image and the reference image; the larger the UIQI value, the better the fusion effect.

(32)$UIQI\left(X,Y\right)=\frac{1}{M}\sum _{i=1}^M\frac{4{\sigma }_{XY}^2·X,Y}{\left({\sigma }_X^2+{\sigma }_Y^2\right)\left[{\left(X\right)}^2+{\left(Y\right)}^2\right]}.$

4.4. Parametric Analysis

The parameters (${\mu }_1$, ${\mu }_2$, ${\mu }_3$) control the balance between data fidelity and regularization terms during optimization. Parameter tuning is usually performed by trying different values and choosing to produce visually accurate results. Therefore, in experiments, tuning of the penalty parameters is usually required to find the optimal parameter values in order to obtain the best model performance. Figure 2a,b shows the PSNR values of the reconstructed HSI for different parameter fetches on the University of Pavia dataset, respectively. As shown in Figure 2a,b, the PSNR values remain relatively stable when ${\mu }_1$ and ${\mu }_2$ take values in the interval $[1\times {10}^{-8},1\times {10}^{-4}]$. Therefore, the parameter values were set to ${\mu }_1=1\times {10}^{-6}$ and ${\mu }_2=1\times {10}^{-6}$. Figure 2c shows the PSNR values of the reconstructed HSI at different fetches on the University of Pavia dataset. The horizontal axis in Figure 2c indicates the value taken by ${\mu }_3$. The PSNR value varies more when the parameter takes values in the interval $[1\times {10}^{-2},1\times {10}^{-1}]$, it varies less when it is in the interval $[1\times {10}^{-4},1\times {10}^{-3}]$, and it stabilizes when it is in the interval $[1\times {10}^{-8},1\times {10}^{-4}]$, so the experiment was set with ${\mu }_3=1\times {10}^{-6}$.

The regularization parameters $\alpha$ and $\beta$ govern the spatial and spectral regularization effects during the fusion process. Adjusting these values impacts both the spatial smoothness and spectral quality of the fused image. From Figure 3a,b, our analysis indicates that the experimental results are optimal when the parameter values lie within the interval $[1\times {10}^{-8},1\times {10}^{-4}]$. Consequently, we set the parameters as $\alpha =1\times {10}^{-6}$ and $\beta =1\times {10}^{-6}$.

To assess the algorithm’s convergence, we monitored its evaluation metrics across various iterations, producing the corresponding convergence curves shown in Figure 4a,b. Special attention was paid to the change in RMSE value with the number of iterations, where Figure 4b is a local enlargement of Figure 4 to show the convergence curves more clearly. The curves can be viewed to assess the convergence of the algorithm and to understand its speed and stability performance during the iteration process. The results indicate a stabilization of the RMSE value with an increase in iteration count, highlighting the algorithm’s continuous refinement of the objective function towards an optimal solution. These curves show a relatively smooth trend without obvious oscillations or sharp fluctuations, further proving the stability of the algorithm.

4.5. Experimental Results and Discussion

Figure 5 and Figure 6 display the fusion experiment results from various algorithms on the Pavia University dataset. For enhanced spatial detail visualization, the RGB band for the pseudo-color images was selected as $[61,25,13]$. Figure 5a,b shows the multispectral image after the spectral downsampling process and the hyperspectral image with spatial downsampling, respectively. Figure 5c–i shows super-resolution images after processing by different algorithms. Upon inspecting the comparison images, certain key differences and similarities are evident. Specifically, Figure 5f exhibits variations in image clarity compared to the MSI, evident from the residual noise present in the image. Furthermore, Figure 5g displays darker pixels relative to the MSI. It is evident that the various algorithms have achieved significant improvements in spatial resolution, which represents a substantial advancement compared to the HSI.

To further examine the differences between the fused images, we generated differential graphs. This involved calculating the pixel-by-pixel difference between the fused image and the original hyperspectral image to present the results of the various algorithms more intuitively.

In the differential graphs, if more spatial information is displayed, it indicates that less spatial information was lost during the fusion process. In Figure 6, we can clearly observe that the fused image still has a certain degree of noise, and some spatial details are lost. Specifically, we observe that Figure 6a,b,f exhibits a large loss of spatial information, resulting in large reconstruction errors and a failure to effectively recover the spatial structure of the multispectral images. However, these differential graphs retain more information in some edge parts, indicating that these algorithms have some advantages in edge preservation, although the overall effect needs to be further optimized. In contrast, the details of the houses in Figure 6d,e exhibit large areas of difference, indicating that the images still have more noise. Contrastingly, the differential graphs in Figure 6c,g almost obscure the finer details such as houses and roads, which indicates that the tensor-ring-based method is better at preserving the detail information of the image and provides better image quality. The differential graph in Figure 6g is even fainter compared to that in Figure 6c, which further highlights that the algorithm used in this paper has a greater advantage in dealing with the spatial details of the image.

Figure 7 and Figure 8 show the fusion results of the Washington dataset under different algorithms. To show the spatial details of the images, $[40,30,5]$ was chosen as the RGB band for the pseudo-color images. Figure 7a,b represents the multispectral image after the spectral downsampling process and the hyperspectral image with spatial downsampling, respectively. Figure 7c–h showcases the Washington images post-processed via different algorithms. Several of the algorithms can restore details well when observed in terms of the details of green spaces, roads, and buildings. Assessed in terms of image clarity and edge sharpness, Figure 7c–e, g–i contains edges and details with no blurring or blending, indicating that the algorithms are highly efficient in retaining spatial details.

The differential graphs in Figure 8 show the spatial detail and information lost during the fusion process. One can observe that in Figure 8a,b,e, a significant amount of spatial detail is not retained. By examining the variations in brightness and differences in detail within the differential graphs, it can be inferred that the fusion process resulted in a loss of spatial detail, as evidenced by the larger brightness variations and more noticeable detail differences observed in Figure 8a–c,f. The noise level in the difference image in Figure 8d,e is also clearly shown; the presence of significant noise or noise artefacts in the differential graphs indicates that noise or artefacts were introduced during the fusion process, resulting in a loss of spatial detail. The lack of significant edge blurring or texture distortion in Figure 8g suggests the algorithm’s efficacy in preserving the original image’s spatial structure.

Figure 9 presents the fusion results of the Ningxia Sand Lake dataset, obtained after employing various comparison algorithms, including the ones proposed in this paper. Figure 9a,b represents the MSI and the HSI obtained directly through the sensor and processed by spatial cropping and band cropping, respectively. Figure 9c–h shows the experimental results when processing this dataset using different algorithms. Directly observing the fused images reveals spectral distortion in Figure 9c. Moreover, observation of features like edges, textures, and fine structure in Figure 9f suggests that the algorithm lost significant spatial details during fusion. The fusion results in Figure 9e,g,h closely resemble the reference image, especially in Figure 9h, which appears free from distortion and blurring. Figure 10 demonstrates the fusion results of different algorithms on the Qingtongxia Rice dataset. Among them, Figure 10a,b represents the MSI and the HSI acquired by the sensor and that obtained after spatial cropping and band cropping, respectively. It is evident that Figure 10e,g is closer to the MSI, with the details of the houses and the vehicles retained more completely, and the resultant images from the proposed algorithm showcase spatial textures on roads and roofs that are more akin to the high-resolution images. This proves that the algorithm proposed in this paper performs equally well on real datasets.

Given that the human eye inherently perceives images with a certain bias, combining qualitative with quantitative results is a common method for evaluating image quality. In order to evaluate the effectiveness and accuracy of the algorithm more precisely, experiments were performed with precise comparisons of spatial features, spectral features, and signal-to-noise ratios, and the optimal values of the experimental results are shown in bold.

According to the experimental results shown in Table 1, the algorithm in this paper demonstrates excellent performance on the Pavia University dataset. Specifically, regarding the PSNR index, our algorithm suppresses noise more effectively than other algorithms, rendering the fused image closer to the original. The PSNR, an evaluation metric for overall image quality, underscores the superior performance of our algorithm in comparative experiments. Compared with other algorithms, the PSNR value of this paper’s algorithm is higher, which means that the fused image shows greater improvements in overall quality and in reducing image noise and improving image clarity. In addition, the evaluation metrics of spectral information, such as RMSE, ERGAS, SAM, DD, and CC, also approach optimal values. Although the SSIM value of NLSTF is slightly better than ours, our approach exhibits only a minor reduction in retaining spectral information. Additionally, our algorithm closely approaches optimal values in metrics related to spatial information, like SSIM and UIQI, which implies that the fused image has a significantly better preservation of spatial details and structures.

According to the detailed data analysis in Table 2, this paper’s algorithm shows superior results in describing the spectral information on the Washington dataset. In particular, key metrics such as SAM and CC are close to the optimal values, indicating that the algorithm in this paper shows excellent performance in maintaining the consistency and relevance of the spectral information during the image fusion process. Although our approach does not achieve optimal values in error metrics such as RMSE, ERGAS, and DD in the comparative experiments, which are mainly used to quantify the error and distortion level of image reconstruction, the algorithm in this paper still achieves satisfactory results in these aspects and successfully reduces the inter-band RMSE and the overall distortion level of the image. Furthermore, the PSNR, a classic metric in image quality evaluation, shows that our algorithm outperforms others in the comparative experiments, indicating that the fused image is closer to the original image with higher reconstruction fidelity compared to other algorithms. Meanwhile, the UIQI value also reaches the optimal level, which further confirms the excellent performance of this paper’s algorithm in maintaining image clarity and contrast. Although the algorithm is slightly inferior to CTRF in SSIM, it still achieves satisfactory results. SSIM is a comprehensive image quality assessment index that takes into account information on brightness, contrast, and structure. Although our algorithm does not achieve the optimal SSIM value, its performance remains within an acceptable range, indicating a commendable level of overall image quality and visual perception. In summary, the experimental results of this paper’s algorithm on the Washington dataset show that this paper’s method exhibits good performance in terms of both spectral information and image quality.

Based on the detailed quantitative metrics presented in Table 3, it can be concluded that the algorithm proposed in this paper results in less spatial information loss on the Ningxia Sand Lake dataset and demonstrates high spectral consistency with the source image. This indicates that the fused image is more consistent with the original image, as it maintains spectral information and reduces the degree of spectral distortion. Furthermore, the algorithm proposed in this study demonstrates the highest PSNR value in relation to the signal-to-noise ratio, thus evidencing its capability of efficiently attenuating noise in images and enhancing their quality and clarity. Moreover, to gain a deeper insight into the performance of our algorithm when processing larger datasets with an increased number of pixels, the data presented in Table 4 should be analyzed. Although the reference image comparison algorithm and the algorithm proposed in this paper share similarities in terms of spatial features for multispectral images, the latter exhibits superior performance regarding spectral information and consistency. Metrics such as RMSE and ERGAS show relatively lower values, highlighting the strength of this algorithm. This suggests that the algorithm presented in this paper accurately preserves spectral information and closely approximates the original dataset when applied to real dataset. This provides further evidence of the algorithm’s effectiveness.

The algorithm proposed in this paper exhibits a small loss in spatial information on real datasets with a high spectral coherence and signal-to-noise ratio. It is worth noting that in some comparison experiments, we observed that both the reference image comparison algorithm and the algorithm proposed in this paper were relatively close to the multispectral image, but there were some outlier cases. However, the algorithm in this paper is still able to maintain better results when processing data with more pixels, especially in terms of spectral characterization metrics such as RMSE, ERGAS, and SAM without outliers, which demonstrates that the algorithm in this paper exhibits robustness and stability when dealing with real datasets.

4.6. Ablation Experiments on Regularization Terms

In order to verify the validity of the graph regular term and the spatial smoothing regular term in the model proposed in this paper, ablation experiments were carried out on the Pavia University dataset. The effects of the regular terms and the changes in the algorithm performance were observed through the ablation experiments to quantify the contribution of the regular terms to the performance of the algorithm, so as to validate the effectiveness of the two regular terms separately and provide further experimental support for the algorithm. In Table 5, the comparison of quantitative metrics results between the NNTR (no spatial smoothing constraints, no graph regularity constraints) method and the NSTR (no spatial smoothing constraints, with graph regularity constraints) method proves that graph regularization plays a positive role in the algorithmic model and has a significant impact on the algorithmic performance. The experimental results show that after the introduction of graph regularization, the quantitative index results are significantly improved, limited not solely to spatial or spectral features, but reflecting a certain degree of improvement in both. This implies that graph regularization helps to better combine spatial and spectral information in hyperspectral image processing and improves the performance of the algorithm, which verifies its effectiveness. The comparison between the results of NSTR and GRTR verifies the effectiveness of spatial smoothing constraints, and GRTR with spatial smoothing constraints shows better performance in image reconstruction. The spatial smoothing constraints not only reduce the noise in the image to a certain extent and retain more structural information but also improve the peak signal-to-noise ratio of the image. This set of comparisons indicates that spatial smoothing constraints play an important role in improving the quality and performance of remote sensing image processing.

Collectively, the algorithm proposed in this paper shows good performance in terms of spectral and spatial information. The experimental results on both the simulated and real datasets show that, compared to other comparative algorithms, image fusion using the algorithm in this paper retains a large amount of spatial detail and maintains high-quality spectral information. The experimental results also demonstrate the excellent fusion performance of the algorithm on both the simulated and real datasets. The experimental results on the simulated dataset show that the algorithm is able to process the image data well and maintain its quality, while the experimental results on the real dataset also confirm the potential and feasibility of the algorithm for practical applications. In conclusion, the algorithm proposed in this paper has great potential and practicality in hyperspectral image fusion and can improve the quality of and detail retention in fused images while retaining spectral and spatial information.

5. Summary

In response to the limitations of hyperspectral and multispectral images, this paper proposes a fusion algorithm based on tensor ring decomposition and graph regularization. By transforming the image into a tensor representation and applying tensor ring decomposition, the process of decomposition and reconstruction is achieved. Additionally, graph regularization techniques and spatial smoothness constraints are introduced, enabling image smoothing and denoising. This enhances the quality and spectral fidelity of the fused image, resulting in a highly spatially resolved fused image. The experimental results validate the effectiveness and superiority of the proposed algorithm, showing significant improvements in fused image quality and spectral information retention. In summary, this research provides an effective method for fusing hyperspectral and multispectral images and has potential applications in the field of remote sensing image processing. Future research can further explore the performance and applicability of the algorithm and apply it to a wider range of remote sensing image processing tasks.

Despite the advantages of the tensor ring model in hyperspectral image fusion, there are some limitations regarding its computational complexity, parameter tuning, loss of spatial information, data dependency, and generality. When applying the model, these limitations need to be fully considered, evaluated, and adjusted according to actual needs. Future directions for improving the tensor ring model in the field of hyperspectral imagery include improving the model structure, introducing spatiotemporal information, combining deep learning and image generation techniques, incorporating cross-domain and cross-modal fusion, and applying reinforcement learning and adaptive fusion. These improvements will help improve the performance and applicability of the fusion model and drive research progress in the field of hyperspectral imagery.

Author Contributions

Data curation, X.Z. and X.M.; Writing—original draft preparation, S.S.; Writing—review and editing, W.B. and K.Q.; Funding acquisition, S.S., W.B. and W.F. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The dataset for pavia university can be downloaded from: https://engineering.purdue.edu/~biehl/MultiSpec/hyperspectral.html, (accessed on 1 March 2023). The dataset for Washington can be downloaded from: https://engineering.purdue.edu/~biehl/MultiSpec/hyperspectral.html, (accessed on 1 March 2023).

Acknowledgments

The authors would like to thank the Image and Intelligence Information Processing Innovation Team of the National Ethnic Affairs Commission of China for their support.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Illustration of the proposed method.

Figure 2. PSNR values for different parameters on the Pavia University dataset. (a) PSNR values for different [Forumla omitted. See PDF.] values. (b) PSNR values for different [Forumla omitted. See PDF.] values. (c) PSNR values for different [Forumla omitted. See PDF.] values.

Figure 3. PSNR values for different parameters on the Pavia University dataset. (a) PSNR values for different [Forumla omitted. See PDF.] values. (b) PSNR values for different [Forumla omitted. See PDF.] values.

Figure 4. The iteration curve. (a) The iteration curve with respect to the RMSE. (b) Local iteration curves with respect to the RMSE.

Figure 5. Comparative convergence results for the Pavia University dataset.

Figure 6. Comparison of differential graphs on the Pavia University dataset.

Figure 7. Comparative convergence results for the Washington dataset.

Figure 8. Comparison of differential graphs on the Washington dataset.

Figure 9. Fusion result comparison on the Ningxia Sand Lake dataset.

Figure 10. Fusion result comparison on the Qingtongxia Rice dataset.

Table 1

Results of quantitative metrics tested on the Pavia University dataset.

Method	PSNR	RMSE	ERGAS	SAM	UIQI	SSIM	DD	CC	TIME
$\mathbf{Best}$	$+\infty$	0	0	0	1	1	0	1	0
$\mathbf{HySure}$	39.7597	2.8047	1.5769	2.6782	0.9885	0.9821	1.8766	0.9943	31.7137
$\mathbf{CNMF}$	34.694	4.7239	2.8239	3.1308	0.9677	0.9592	3.1644	0.9879	14.5149
$\mathbf{CTRF}$	42.7362	2.0646	1.1527	2.1830	0.9922	0.9847	1.3576	0.9964	47.0982
$\mathbf{SCOTT}$	27.0018	11.6015	6.4980	6.9250	0.7565	0.6924	8.1506	0.8883	2.0046
$\mathbf{STEREO}$	36.7565	4.1604	2.3230	4.0671	0.9762	0.9596	2.7784	0.9864	9.4420
$\mathbf{NLSTF}$	42.0759	2.1348	1.2417	2.1593	0.9920	0.9865	1.3627	0.9963	83.4873
$\mathbf{Proposed}$	43.0236	1.9835	1.1139	2.0903	0.9928	0.9859	1.2977	0.9967	97.9190

The optimal value is in bold.

Table 2

Results of quantitative metrics tested on the Washington dataset.

Method	PSNR	RMSE	ERGAS	SAM	UIQI	SSIM	DD	CC	TIME
$\mathbf{Best}$	$+\infty$	0	0	0	1	1	0	1	0
$\mathbf{HySure}$	42.9681	2.0904	2.5829	2.8028	0.9863	0.9885	1.441	0.9918	31.4016
$\mathbf{CNMF}$	36.0655	5.0024	4.0038	3.8821	0.9301	0.9274	3.5911	0.9735	17.0852
$\mathbf{CTRF}$	47.4889	1.1991	3.1830	1.6988	0.9899	0.9937	0.7693	0.9929	46.6742
$\mathbf{SCOTT}$	31.4566	8.0958	6.1754	7.0211	0.8259	0.7922	5.7558	0.9111	2.4571
$\mathbf{STEREO}$	40.2923	2.7368	21.2904	3.9808	0.9577	0.9725	2.0252	0.9741	9.3929
$\mathbf{NLSTF}$	47.0276	1.2635	3.1164	1.7331	0.9896	0.9934	0.7903	0.9929	74.7197
$\mathbf{Proposed}$	47.4898	1.2167	2.6451	1.6970	0.9910	0.9934	0.7797	0.9941	98.2590

The optimal value is in bold.

Table 3

Results of quantitative metrics tested on the Ningxia Sand Lake dataset.

Method	PSNR	RMSE	ERGAS	SAM	UIQI	SSIM	DD	CC	TIME
$\mathbf{Best}$	$+\infty$	0	0	0	1	1	0	1	0
$\mathbf{HySure}$	21.7442	36.9608	54.6875	24.8293	0.6557	0.6538	22.9647	0.753	53.313
$\mathbf{CNMF}$	44.4076	4.5102	5.9954	2.2275	0.9848	0.9803	1.7758	0.9944	22.6497
$\mathbf{CTRF}$	48.1040	1.1031	2.9553	1.0954	0.9972	0.9973	0.7648	0.9996	106.4077
$\mathbf{SCOTT}$	29.129	14.4249	22.6336	9.0206	0.7852	0.6964	8.6452	0.9232	1.2054
$\mathbf{STEREO}$	26.6772	15.8444	43.2869	9.1986	0.9000	0.8935	11.8486	0.9944	19.2988
$\mathbf{Proposed}$	48.1946	1.0946	2.9362	1.0795	0.9972	0.9973	0.7569	0.9996	178.637

The optimal value is in bold.

Table 4

Results of quantitative metrics tested on the Qingtongxia Rice dataset.

Method	PSNR	RMSE	ERGAS	SAM	UIQI	SSIM	DD	CC	TIME
$\mathbf{Best}$	$+\infty$	0	0	0	1	1	0	1	0
$\mathbf{HySure}$	15.5155	44.1258	111.2612	36.6818	0.2143	0.3622	29.7187	0.9307	206.4664
$\mathbf{CNMF}$	24.5566	16.0158	31.8906	5.075	0.3737	0.7248	9.2112	0.943	62.4262
$\mathbf{CTRF}$	39.9188	2.9146	6.6945	2.1332	0.9562	0.9896	1.8246	0.9993	217.8644
$\mathbf{SCOTT}$	29.7305	8.5474	17.7987	3.9481	0.4565	0.825	5.4843	0.9833	1.8374
$\mathbf{STEREO}$	14.9351	48.082	123.9833	11.7194	0.3458	0.6618	30.4537	0.9496	57.5026
$\mathbf{Proposed}$	40.0091	2.9002	6.6672	2.0442	0.959	0.9902	1.7999	0.9993	356.521

The optimal value is in bold.

Table 5

Comparative results of ablation experiments on Pavia University dataset.

Method	Signal to Noise Ratio	Spectral Feature					Spatial Feature
	PSNR	RMSE	ERGAS	SAM	DD	CC	UIQI	SSIM
NNTR	40.6217	2.6205	1.4932	2.8833	1.8088	0.9941	0.9868	0.9722
NSTR	42.9886	1.9873	1.1163	2.1308	1.3188	0.9967	0.9927	0.9852
GRTR	43.0236	1.9835	1.1139	2.0903	1.2977	0.9967	0.9928	0.9859

Appendix A

Appendix A.1. Optimization Process for ${\mathcal{G}}^{\left(1\right)}$

The optimization of ${\mathcal{G}}^{\left(1\right)}$ can be achieved by iteratively updating ${\mathcal{G}}^{\left(1\right)}$, ${\mathcal{R}}_1$, and ${\mathcal{O}}_1$ through three sub-problems, as given in Equation (19):(A1)$$\left\{\begin{array}{c}{\mathcal{G}}^{\left(1\right),i+1}=\underset{{\mathcal{G}}^{\left(1\right)}}{argmin}\Gamma ({\mathcal{G}}^{\left(1\right)},{{\mathcal{R}}_1}^i,{{\mathcal{O}}_1}^i)\hfill \\ {\mathcal{R}}_1=\underset{{\mathcal{R}}_1}{argmin}\Gamma ({\mathcal{G}}^{\left(1\right),i+1},{\mathcal{R}}_1,{{\mathcal{O}}_1}^i)\hfill \\ {\mathcal{O}}_1=\underset{{\mathcal{O}}_1}{argmin}\Gamma ({\mathcal{G}}^{\left(1\right),i+1},{{\mathcal{R}}_1}^{i+1},{\mathcal{O}}_1)\hfill \end{array}\right.$$

Algorithm A1 Solving the ${\mathcal{G}}^{\left(1\right)}$ sub-problem

Input MSI $\mathcal{Z}$, HSI $\mathcal{Y}$, ${\mathcal{G}}^{\left(2\right)}$, ${\mathcal{G}}^{\left(3\right)}$, regularization parameters $\lambda$, $\alpha$, $\beta$, $\rho$ and ${\mu }_1$;  let ${\mathcal{R}}_1={\mathcal{G}}^{\left(1\right)}$;  $\mathbf{While}$  $\mathbf{not}\mathbf{coverged}$  $\mathbf{do}$         Updating the ${\mathcal{G}}^{\left(1\right)}$ via solving Equation (A4);         Updating the ${\mathcal{R}}_1$ via solving Equation (A8);         Updating the ${\mathcal{O}}_1$ via solving Equation (A9);  $\mathbf{end}$  $\mathbf{While}$ Output ${\mathcal{G}}^{\left(1\right)}$.

Appendix A.1.1. Optimization Process for ${\mathcal{G}}^{\left(1\right)}$

The subproblem on ${\mathcal{G}}^{\left(1\right)}$ can be formulated as (A2)$$\underset{{\mathcal{G}}^{\left(1\right)}}{min}\frac{1}{2}{∥\mathcal{Y}-{\mathcal{G}}^{\left(1\right)}{\times }_2{U}_1{\mathcal{P}}_1∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-{\mathcal{G}}^{\left(1\right)}{\mathcal{P}}_2∥}_F^2+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(1\right)}-{\mathcal{G}}^{\left(1\right),i}∥}_F^2+\frac{{\mu }_1}{2}{∥{\mathcal{G}}^{\left(1\right)}-{\mathcal{R}}_1-{\mathcal{O}}_1∥}_F^2.$$ where $\mathcal{P}1={({\mathcal{G}}^{\left(2\right)}{\times }_2{U}_2,{\mathcal{G}}^{\left(3\right)})}_{<2>}^T$, $\mathcal{P}2={({\mathcal{G}}^{\left(2\right)},{\mathcal{G}}^{\left(3\right)}{\times }_2{U}_3)}_{<2>}^T$. Utilizing the property of matrix unfolding, the sub-problem for ${\mathcal{G}}^{\left(1\right)}$ can be expressed as follows:(A3)$$\begin{array}{c}\underset{G_{<2>}^{\left(1\right)}}{min}\frac{1}{2}{∥Y_{<1>}-U_1{G}_{<2>}^{\left(1\right)}{P}_1∥}_F^2+\frac{\lambda }{2}{∥Z_{<1>}-G_{<2>}^{\left(1\right)}{P}_2∥}_F^2\hfill \\ +\frac{\rho }{2}{∥G_{<2>}^{\left(1\right)}-G_{<2>}^{\left(1\right),i}∥}_F^2+\frac{{\mu }_1}{2}{∥G_{<2>}^{\left(1\right)}-R_1-O_1∥}_F^2\hfill \end{array}$$

This problem is a quadratic optimization problem with a unique solution, equivalent to computing the following Sylvester matrix equation, which can be solved efficiently using the Conjugate Gradient (CG) method. (A4)$$\begin{array}{c}{U}_{{}_1}^T{U}_1{G}_{<2>}^{\left(1\right)}{P}_1{P}_{{}_1}^T+\lambda G_{<2>}^{\left(1\right)}{P}_2{P}_2^T+({\mu }_1+\rho )G_{<2>}^{\left(1\right)}\hfill \\ =\rho G_{<2>}^{\left(1\right),i}+U_1^T{Y}_{<1>}{P}_1^T+\lambda Z_{<1>}{P}_2^T+{\mu }_1{R}_1+{\mu }_1{O}_1\hfill \end{array}$$

Appendix A.1.2. Optimization Process for ${\mathcal{R}}_1$

The subproblem on ${\mathcal{R}}_1$ can be formulated as (A5)$$\underset{{\mathcal{R}}_1}{argmin}\frac{\beta }{2}tr\left(\right({\mathcal{G}}^{\left(2\right),i}\otimes {\mathcal{R}}_1{)}^T{L}_{SPA}({\mathcal{G}}^{\left(2\right),i}\otimes {\mathcal{R}}_1\left)\right)+{∥D_1{\mathcal{R}}_1∥}_2^2+\frac{{\mu }_1}{2}{∥{\mathcal{G}}^{\left(1\right)}-{\mathcal{R}}_1-{\mathcal{O}}_1∥}_F^2.$$

Using the properties of matrix expansion, the subproblem on ${\mathcal{R}}_1$ can be expressed as (A6)$$\underset{R_1}{argmin}\frac{\beta }{2}tr\left(\right(G_{<2>}^{\left(2\right)}\otimes R_1{)}^T{L}_{SPA}(G_{<2>}^{\left(2\right)}\otimes R_1\left)\right)+{∥D_1{R}_1∥}_2^2+\frac{{\mu }_1}{2}{∥G_{<2>}^{\left(1\right)}-R_1-O_1∥}_F^2$$

Cholesky decomposition [51] of $L_{SPA}$ leads to $L_{SPA}={\Sigma }^T\Sigma$, so the problem is transformed into the following equation for solving:(A7)$$\underset{R_1}{argmin}\frac{\beta }{2}{∥\Sigma (G_{<2>}^{\left(2\right)}\otimes R_1)∥}_F^2+{∥D_1{R}_1∥}_2^2+\frac{{\mu }_1}{2}{∥G_{<2>}^{\left(1\right)}-R_1-O_1∥}_F^2.$$

Taking the derivative of the above equation and making the other derivative 0, we obtain the solution with respect to $R_1$:(A8)$$R_1=(\beta {\Sigma }^{T}G{}_{<2>}^{\left(2\right)T}{G}_{<2>}^{\left(2\right)}\otimes I+{\mu }_{1}I+2D_1{)}^{-1}({\mu }_1{G}_{<2>}^{\left(1\right)}-{\mu }_1{O}_1).$$

Appendix A.1.3. Optimization Process for ${\mathcal{O}}_1$

The Lagrange multipliers ${\mathcal{O}}_1$ update is:(A9)$${\mathcal{O}}_1={\mathcal{O}}_1-({\mathcal{G}}^{\left(1\right)}-{\mathcal{R}}_1).$$

Appendix A.1.4. Optimization Process for ${\mathcal{G}}^{\left(2\right)}$

The optimization of ${\mathcal{G}}^{\left(2\right)}$ can be achieved by iteratively updating ${\mathcal{G}}^{\left(2\right)}$, ${\mathcal{R}}_2$, and ${\mathcal{O}}_2$ through three sub-problems, as given in Equation (21):(A10)$$\left\{\begin{array}{c}{\mathcal{G}}^{\left(2\right),i+1}=\underset{{\mathcal{G}}^{\left(2\right)}}{argmin}\Gamma ({\mathcal{G}}^{\left(2\right)},{{\mathcal{R}}_2}^i,{{\mathcal{O}}_2}^i)\hfill \\ {\mathcal{R}}_2=\underset{{\mathcal{R}}_2}{argmin}\Gamma ({\mathcal{G}}^{\left(2\right),i+1},{\mathcal{R}}_2,{{\mathcal{O}}_2}^i)\hfill \\ {\mathcal{O}}_2=\underset{{\mathcal{O}}_2}{argmin}\Gamma ({\mathcal{G}}^{\left(2\right),i+1},{{\mathcal{R}}_2}^{i+1},{\mathcal{O}}_2)\hfill \end{array}\right.$$

Algorithm A2 Solving the ${\mathcal{G}}^{\left(2\right)}$ sub-problem

Input MSI $\mathcal{Z}$, HSI $\mathcal{Y}$, ${\mathcal{G}}^{\left(1\right),i+1}$, ${\mathcal{G}}^{\left(3\right)}$, regularization parameters $\lambda$, $\alpha$, $\beta$, $\rho$ and ${\mu }_2$;  let ${\mathcal{R}}_2={\mathcal{G}}^{\left(2\right)}$;  $\mathbf{While}$  $\mathbf{not}\mathbf{coverged}$  $\mathbf{do}$         Updating the ${\mathcal{G}}^{\left(2\right)}$ via solving Equation (A13);         Updating the ${\mathcal{R}}_2$ via solving Equation (A17);         Updating the ${\mathcal{O}}_2$ via solving Equation (A18);  $\mathbf{end}$  $\mathbf{While}$ Output ${\mathcal{G}}^{\left(2\right)}$.

Appendix A.1.5. Optimization Process for ${\mathcal{G}}^{\left(2\right)}$

The subproblem on ${\mathcal{G}}^{\left(2\right)}$ can be formulated as (A11)$$\underset{{\mathcal{G}}^{\left(2\right)}}{min}\frac{1}{2}{∥\mathcal{Y}-{\mathcal{G}}^{\left(2\right)}{\times }_2{U}_2{\mathcal{Q}}_1∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-{\mathcal{G}}^{\left(2\right)}{\mathcal{Q}}_2∥}_F^2+\frac{\rho }{2}{∥{\mathcal{G}}^{\left(2\right)}-{\mathcal{G}}^{\left(2\right),i}∥}_F^2+\frac{{\mu }_2}{2}{∥{\mathcal{G}}^{\left(2\right)}-{\mathcal{R}}_2-{\mathcal{O}}_2∥}_F^2.$$

Using the properties of matrix expansion, the subproblem on ${\mathcal{G}}^{\left(2\right)}$ can be expressed as (A12)$$\begin{array}{c}\underset{G_{<2>}^{\left(2\right)}}{min}\frac{1}{2}{∥Y_{<1>}-U_2{G}_{<2>}^{\left(2\right)}{Q}_1∥}_F^2+\frac{\lambda }{2}{∥Z_{<1>}-G_{<2>}^{\left(2\right)}{Q}_2∥}_F^2\hfill \\ +\frac{\rho }{2}{∥G_{<2>}^{\left(2\right)}-G_{<2>}^{\left(2\right),i}∥}_F^2+\frac{{\mu }_2}{2}{∥G_{<2>}^{\left(2\right)}-R_2-O_2∥}_F^2.\hfill \end{array}$$

This problem is a quadratic optimization problem with a unique solution, equivalent to computing the following Sylvester matrix equation, which can be solved efficiently using the CG method. (A13)$$\begin{array}{c}{U}_2^T{U}_2{G}_{<2>}^{\left(2\right)}{Q}_1{Q}_{{}_1}^T+\lambda G_{<2>}^{\left(2\right)}{Q}_2{Q}_2^T+({\mu }_2+\rho )G_{<2>}^{\left(2\right)}\hfill \\ =\rho G_{<2>}^{\left(2\right),i}+U_2^T{Y}_{<1>}{Q}_1^T+\lambda Z_{<1>}{Q}_2^T+{\mu }_2{R}_2+{\mu }_2{O}_2\hfill \end{array}$$

Appendix A.1.6. Optimization Process for ${\mathcal{R}}_2$

The subproblem on ${\mathcal{R}}_2$ can be formulated as (A14)$$\underset{{\mathcal{R}}_2}{argmin}\frac{\beta }{2}tr(({\mathcal{R}}_2\otimes {\mathcal{G}}^{\left(1\right),i+1}{)}^T{L}_{SPA}({\mathcal{R}}_2\otimes {\mathcal{G}}^{\left(1\right),i+1}))+{∥D_2{\mathcal{R}}_2∥}_2^2+\frac{{\mu }_2}{2}{∥{\mathcal{G}}^{\left(2\right)}-{\mathcal{R}}_2-{\mathcal{O}}_2∥}_F^2.$$

Using the properties of matrix expansion, the subproblem on ${\mathcal{R}}_2$ can be expressed as (A15)$$\underset{R_2}{argmin}\frac{\beta }{2}tr((R_2\otimes G_{<2>}^{\left(1\right)}{)}^T{L}_{SPA}(R_2\otimes G_{<2>}^{\left(1\right)}))+{∥D_2{R}_2∥}_2^2+\frac{{\mu }_2}{2}{∥G_{<2>}^{\left(2\right)}-R_2-O_2∥}_F^2$$

According to the Cholesky decomposition, we can obtain:(A16)$$\underset{R_2}{argmin}\frac{\beta }{2}{∥\Sigma (R_2\otimes G_{<2>}^{\left(1\right)})∥}_F^2+{∥D_2{R}_2∥}_2^2+\frac{{\mu }_2}{2}{∥G_{<2>}^{\left(2\right)}-R_2-O_2∥}_F^2$$

Taking the derivative of the above equation and making the other derivative 0, we obtain the solution with respect to $R_2$:(A17)$$R_2={(\beta {\Sigma }^T{G}_{<2>}^{\left(1\right)T}{G}_{<2>}^{\left(1\right)}\otimes I+{\mu }_{2}I+2D_2)}^{-1}({\mu }_2{G}_{<2>}^{\left(2\right)}-{\mu }_2{O}_2)$$

Appendix A.1.7. Optimization Process for ${\mathcal{O}}_2$

The update on Lagrange multipliers ${\mathcal{O}}_2$ is:(A18)$${\mathcal{O}}_2={\mathcal{O}}_2-({\mathcal{G}}^{\left(2\right)}-{\mathcal{R}}_2)$$

Appendix A.1.8. Optimization Process for ${\mathcal{G}}^{\left(3\right)}$

The optimization of ${\mathcal{G}}^{\left(3\right)}$ can be achieved by iteratively updating ${\mathcal{G}}^{\left(3\right)}$, ${\mathcal{R}}_3$, and ${\mathcal{O}}_3$ through three sub-problems, as given in Equation (20):(A19)$$\left\{\begin{array}{c}{\mathcal{G}}^{\left(3\right),i+1}=\underset{{\mathcal{G}}^{\left(3\right)}}{argmin}\Gamma ({\mathcal{G}}^{\left(3\right)},{{\mathcal{R}}_3}^i,{{\mathcal{O}}_3}^i)\hfill \\ {\mathcal{R}}_3=\underset{{\mathcal{R}}_3}{argmin}\Gamma ({\mathcal{G}}^{\left(3\right),i+1},{\mathcal{R}}_3,{{\mathcal{O}}_3}^i)\hfill \\ {\mathcal{O}}_3=\underset{{\mathcal{O}}_3}{argmin}\Gamma ({\mathcal{G}}^{\left(3\right),i+1},{{\mathcal{R}}_3}^{i+1},{\mathcal{O}}_3)\hfill \end{array}\right.$$

Algorithm A3 Solving the ${\mathcal{G}}^{\left(3\right)}$ sub-problem

Input MSI $\mathcal{Z}$, HSI $\mathcal{Y}$, ${\mathcal{G}}^{\left(1\right),i+1}$, ${\mathcal{G}}^{\left(2\right),i+1}$, regularization parameters $\lambda$, $\alpha$, $\beta$, $\rho$ and ${\mu }_3$;  let ${\mathcal{R}}_3={\mathcal{G}}^{\left(3\right)}$;   $\mathbf{While}$  $\mathbf{not}\mathbf{coverged}$  $\mathbf{do}$          Updating the ${\mathcal{G}}^{\left(3\right)}$ via solving Equation (A22);          Updating the ${\mathcal{R}}_3$ via solving Equation (A25);          Updating the ${\mathcal{O}}_3$ via solving Equation (A26);   $\mathbf{end}$  $\mathbf{While}$   Output ${\mathcal{G}}^{\left(3\right)}$.

Appendix A.1.9. Optimization Process for ${\mathcal{G}}^{\left(3\right)}$

The subproblem on ${\mathcal{G}}^{\left(3\right)}$ can be formulated as (A20)$$\begin{array}{c}\underset{{\mathcal{G}}^{\left(3\right)}}{min}\frac{1}{2}{∥\mathcal{Y}-{\mathcal{G}}^{\left(3\right)}{\mathcal{T}}_1∥}_F^2+\frac{\lambda }{2}{∥\mathcal{Z}-{\mathcal{G}}^{\left(2\right)}{\times }_2{U}_3{\mathcal{T}}_2∥}_F^2\hfill \\ +\frac{\rho }{2}{∥{\mathcal{G}}^{\left(3\right)}-{\mathcal{G}}^{\left(3\right),i}∥}_F^2+\frac{{\mu }_3}{2}{∥{\mathcal{G}}^{\left(3\right)}-{\mathcal{R}}_3-{\mathcal{O}}_3∥}_F^2\hfill \end{array}$$ where ${\mathcal{T}}_1={({\mathcal{G}}^{\left(1\right),i+1}{\times }_2{U}_1,{\mathcal{G}}^{\left(2\right),i+1})}_{<2>}^T$, and ${\mathcal{T}}_2={({\mathcal{G}}^{\left(1\right),i+1},{\mathcal{G}}^{\left(2\right),i+1})}_{<2>}^T$. Using the properties of matrix expansion, the subproblem on ${\mathcal{G}}^{\left(3\right)}$ can be expressed as (A21)$$\begin{array}{c}\underset{G_{<2>}^{\left(3\right)}}{min}\frac{1}{2}{∥Y_{<1>}-G_{<2>}^{\left(3\right)}{T}_1∥}_F^2+\frac{\lambda }{2}{∥Z_{<1>}-U_3{G}_{<2>}^{\left(3\right)}{T}_2∥}_F^2\hfill \\ +\frac{\rho }{2}{∥G_{<2>}^{\left(3\right)}-G_{<2>}^{\left(3\right),i}∥}_F^2+\frac{{\mu }_3}{2}{∥G_{<2>}^{\left(3\right)}-R_3-O_3∥}_F^2\hfill \end{array}$$

This problem is a quadratic optimization problem with a unique solution, equivalent to computing the following Sylvester matrix equation, which can be solved efficiently using the CG method. (A22)$$\begin{array}{c}\lambda U_3^T{U}_3{G}_{<2>}^{\left(3\right)}{T}_1{T}_{{}_1}^T+G_{<2>}^{\left(3\right)}{T}_1{T}_1^T+({\mu }_3+\rho )G_{<2>}^{\left(3\right)}\hfill \\ =\rho G_{<2>}^{\left(3\right),i}+Y_{<1>}{T}_1^T+\lambda U_3^T{Z}_{<1>}{T}_2^T+{\mu }_3{R}_3+{\mu }_3{O}_3\hfill \end{array}$$

Appendix A.1.10. Optimization Process for ${\mathcal{R}}_3$

The subproblem on ${\mathcal{R}}_3$ can be formulated as (A23)$$\underset{{\mathcal{R}}_3}{argmin}\frac{\alpha }{2}tr({{\mathcal{R}}_3}^T{L}_{SPE}{\mathcal{R}}_3)+\frac{{\mu }_3}{2}{∥{\mathcal{G}}^{\left(3\right)}-{\mathcal{R}}_3-{\mathcal{O}}_3∥}_F^2$$

Using the properties of matrix expansion, the subproblem on ${\mathcal{R}}_3$ can be expressed as (A24)$$\underset{R_3}{argmin}\frac{\alpha }{2}tr({R_3}^T{P}_S{R}_3)+\frac{{\mu }_3}{2}{∥G_{<2>}^{\left(3\right)}-R_3-O_3∥}_F^2$$

Taking the derivative of the above equation and making the other derivative 0, we obtain the solution with respect to $R_3$:(A25)$$R_3={(\alpha P_S+{\mu }_{3}I)}^{-1}·{\mu }_3(G_{<2>}^{\left(3\right)T}-{O_3}^T)$$

Appendix A.1.11. Optimization Process for ${\mathcal{O}}_3$

The update on Lagrange multipliers ${\mathcal{O}}_3$ is:(A26)$${\mathcal{O}}_3={\mathcal{O}}_3-({\mathcal{G}}^{\left(3\right)}-{\mathcal{R}}_3)$$

References

1. Loncan, L.; de Almeida, L.B.; Bioucas-Dias, J.M.; Briottet, X.; Chanussot, J.; Dobigeon, N.; Fabre, S.; Liao, W.; Licciardi, G.A.; Simões, M. et al. Hyperspectral Pansharpening: A Review. IEEE Geosci. Remote Sens. Mag.; 2015; 3, pp. 27-46. [DOI: https://dx.doi.org/10.1109/MGRS.2015.2440094]

2. Vivone, G.; Restaino, R.; Licciardi, G.; Mura, M.D.; Chanussot, J. MultiResolution Analysis and Component Substitution techniques for hyperspectral Pansharpening. Proceedings of the 2014 IEEE Geoscience and Remote Sensing Symposium; Quebec City, QC, Canada, 13–18 July 2014; pp. 2649-2652. [DOI: https://dx.doi.org/10.1109/IGARSS.2014.6947018]

3. Palsson, F.; Sveinsson, J.R.; Ulfarsson, M.O. Optimal Component Substitution and Multi-Resolution Analysis Pansharpening Methods Using a Convolutional Neural Network. Proceedings of the IGARSS 2019 IEEE International Geoscience and Remote Sensing Symposium; Yokohama, Japan, 28 July–2 August 2019; pp. 3177-3180. [DOI: https://dx.doi.org/10.1109/IGARSS.2019.8899299]

4. Wang, W.; Jiao, L.; Yang, S. Novel Adaptive Component-Substitution-Based Pan-Sharpening Using Particle Swarm Optimization. IEEE Geosci. Remote Sens. Lett.; 2015; 12, pp. 781-785. [DOI: https://dx.doi.org/10.1109/LGRS.2014.2361834]

5. Wang, P.; Yao, H.; Huang, B.; Leung, H.; Liu, P. Multiresolution Analysis Pansharpening Based on Variation Factor for Multispectral and Panchromatic Images From Different Times. IEEE Trans. Geosci. Remote Sens.; 2023; 61, 5401217. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3252001]

6. Alparone, L.; Baronti, S.; Aiazzi, B.; Garzelli, A. Spatial Methods for Multispectral Pansharpening: Multiresolution Analysis Demystified. IEEE Trans. Geosci. Remote Sens.; 2016; 54, pp. 2563-2576. [DOI: https://dx.doi.org/10.1109/TGRS.2015.2503045]

7. Ma, W.; Shen, J.; Zhu, H.; Zhang, J.; Zhao, J.; Hou, B.; Jiao, L. A Novel Adaptive Hybrid Fusion Network for Multiresolution Remote Sensing Images Classification. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 5400617. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3062142]

8. Liang, W.; Wu, Y.; Li, M.; Cao, Y. High-Resolution SAR Image Classification Using Context-Aware Encoder Network and Hybrid Conditional Random Field Model. IEEE Trans. Geosci. Remote Sens.; 2020; 58, pp. 5317-5335. [DOI: https://dx.doi.org/10.1109/TGRS.2019.2963699]

9. Zhou, X.; Liu, J.; Liu, S.; Cao, L.; Zhou, Q.; Huang, H. A GIHS-based spectral preservation fusion method for remote sensing images using edge restored spectral modulation. ISPRS J. Photogramm. Remote Sens.; 2014; 88, pp. 16-27. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2013.11.011]

10. Gonzalez-Audicana, M.; Saleta, J.; Catalan, R.; Garcia, R. Fusion of multispectral and panchromatic images using improved IHS and PCA mergers based on wavelet decomposition. IEEE Trans. Geosci. Remote Sens.; 2004; 42, pp. 1291-1299. [DOI: https://dx.doi.org/10.1109/TGRS.2004.825593]

11. Alparone, L.; Aiazzi, B.; Baronti, S.; Garzelli, A. Fast classified pansharpening with spectral and spatial distortion optimization. Proceedings of the 2012 IEEE International Geoscience and Remote Sensing Symposium; Munich, Germany, 22–27 July 2012; pp. 154-157. [DOI: https://dx.doi.org/10.1109/IGARSS.2012.6351614]

12. Batur, E.; Maktav, D. Assessment of Surface Water Quality by Using Satellite Images Fusion Based on PCA Method in the Lake Gala, Turkey. IEEE Trans. Geosci. Remote Sens.; 2019; 57, pp. 2983-2989. [DOI: https://dx.doi.org/10.1109/TGRS.2018.2879024]

13. Kong, Y.; Hong, F.; Leung, H.; Peng, X. A Fusion Method of Optical Image and SAR Image Based on Dense-UGAN and Gram–Schmidt Transformation. Remote Sens.; 2021; 13, 4274. [DOI: https://dx.doi.org/10.3390/rs13214274]

14. Dong, W.; Zhang, T.; Qu, J.; Xiao, S.; Liang, J.; Li, Y. Laplacian Pyramid Dense Network for Hyperspectral Pansharpening. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 5507113. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3076768]

15. Zhang, Y.; Liu, F.; Jiao, L.; Yang, S.; Li, L.; Yang, M.; Wang, J.; Liu, X. Curvelet Adversarial Augmented Neural Network for SAR Image Classification. IEEE Trans. Geosci. Remote Sens.; 2023; 61, 4400717. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3239226]

16. Ghahremani, M.; Ghassemian, H. A Compressed-Sensing-Based Pan-Sharpening Method for Spectral Distortion Reduction. IEEE Trans. Geosci. Remote Sens.; 2016; 54, pp. 2194-2206. [DOI: https://dx.doi.org/10.1109/TGRS.2015.2497309]

17. Li, X.; Zhang, Y.; Ge, Z.; Cao, G.; Shi, H.; Fu, P. Adaptive Nonnegative Sparse Representation for Hyperspectral Image Super-Resolution. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2021; 14, pp. 4267-4283. [DOI: https://dx.doi.org/10.1109/JSTARS.2021.3072044]

18. Karoui, M.S.; Deville, Y.; Benhalouche, F.Z.; Boukerch, I. Hypersharpening by Joint-Criterion Nonnegative Matrix Factorization. IEEE Trans. Geosci. Remote Sens.; 2017; 55, pp. 1660-1670. [DOI: https://dx.doi.org/10.1109/TGRS.2016.2628889]

19. Shen, D.; Liu, J.; Wu, Z.; Yang, J.; Xiao, L. ADMM-HFNet: A Matrix Decomposition-Based Deep Approach for Hyperspectral Image Fusion. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 5513417. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3112181]

20. Fu, X.; Jia, S.; Xu, M.; Zhou, J.; Li, Q. Fusion of Hyperspectral and Multispectral Images Accounting for Localized Inter-Image Changes. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 5517218. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3124240]

21. Xue, J.; Zhao, Y.Q.; Bu, Y.; Liao, W.; Chan, J.C.W.; Philips, W. Spatial-Spectral Structured Sparse Low-Rank Representation for Hyperspectral Image Super-Resolution. IEEE Trans. Image Process.; 2021; 30, pp. 3084-3097. [DOI: https://dx.doi.org/10.1109/TIP.2021.3058590]

22. Xu, Y.; Wu, Z.; Chanussot, J.; Comon, P.; Wei, Z. Nonlocal Coupled Tensor CP Decomposition for Hyperspectral and Multispectral Image Fusion. IEEE Trans. Geosci. Remote Sens.; 2020; 58, pp. 348-362. [DOI: https://dx.doi.org/10.1109/TGRS.2019.2936486]

23. Jin, D.; Liu, J.; Yang, J.; Wu, Z. High-Order Coupled Fully Connected Tensor Network Decomposition for Hyperspectral Image Super-Resolution. IEEE Geosci. Remote Sens. Lett.; 2022; 19, 6014105. [DOI: https://dx.doi.org/10.1109/LGRS.2022.3207548]

24. Zhang, K.; Wang, M.; Yang, S.; Jiao, L. Spatial–Spectral-Graph-Regularized Low-Rank Tensor Decomposition for Multispectral and Hyperspectral Image Fusion. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2018; 11, pp. 1030-1040. [DOI: https://dx.doi.org/10.1109/JSTARS.2017.2785411]

25. Zeng, H.; Xue, J.; Luong, H.Q.; Philips, W. Multimodal Core Tensor Factorization and Its Applications to Low-Rank Tensor Completion. IEEE Trans. Multimed.; 2022; 3, pp. 1-15. [DOI: https://dx.doi.org/10.1109/TMM.2022.3216746]

26. Feng, J.; Li, X.; Liu, X.; Chen, C.; Chen, H. Seismic Data Denoising Based on Tensor Decomposition with Total Variation. IEEE Geosci. Remote Sens. Lett.; 2021; 18, pp. 1303-1307. [DOI: https://dx.doi.org/10.1109/LGRS.2021.3054765]

27. Qin, C.; Wu, Z.; Zhang, Y.; Sun, J.; Xu, Y.; Zheng, P.; Wei, Z. A Distributed Parallel Optimization of Remote Sensing Image Fusion Algorithm Based on Nonlocal Tensor CP Decomposition. Proceedings of the IGARSS 2022 IEEE International Geoscience and Remote Sensing Symposium; Kuala Lumpur, Malaysia, 17–22 July 2022; pp. 2370-2373. [DOI: https://dx.doi.org/10.1109/IGARSS46834.2022.9884571]

28. Borsoi, R.A.; Prévost, C.; Usevich, K.; Brie, D.; Bermudez, J.C.M.; Richard, C. Coupled Tensor Decomposition for Hyperspectral and Multispectral Image Fusion With Inter-Image Variability. IEEE J. Sel. Top. Signal Process.; 2021; 15, pp. 702-717. [DOI: https://dx.doi.org/10.1109/JSTSP.2021.3054338]

29. Bu, Y.; Zhao, Y.; Xue, J.; Chan, J.C.W.; Kong, S.G.; Yi, C.; Wen, J.; Wang, B. Hyperspectral and Multispectral Image Fusion via Graph Laplacian-Guided Coupled Tensor Decomposition. IEEE Trans. Geosci. Remote Sens.; 2021; 59, pp. 648-662. [DOI: https://dx.doi.org/10.1109/TGRS.2020.2992788]

30. Solgi, R.; Loaiciga, H.A.; Zhang, Z. Evolutionary Tensor Train Decomposition for Hyper-Spectral Remote Sensing Images. Proceedings of the IGARSS 2022–2022 IEEE International Geoscience and Remote Sensing Symposium; Kuala Lumpur, Malaysia, 17–22 July 2022; pp. 1145-1148. [DOI: https://dx.doi.org/10.1109/IGARSS46834.2022.9884813]

31. Liu, D.; Sacchi, M.D.; Chen, W. Efficient Tensor Completion Methods for 5-D Seismic Data Reconstruction: Low-Rank Tensor Train and Tensor Ring. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 5916317. [DOI: https://dx.doi.org/10.1109/TGRS.2022.3179275]

32. Zhang, H.; Huang, T.Z.; Zhao, X.L.; He, W.; Choi, J.K.; Zheng, Y.B. Hyperspectral Image Denoising: Reconciling Sparse and Low-Tensor-Ring-Rank Priors in the Transformed Domain. IEEE Trans. Geosci. Remote Sens.; 2023; 61, 5502313. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3237865]

33. Xu, H.; Qin, M.; Chen, S.; Zheng, Y.; Zheng, J. Hyperspectral-Multispectral Image Fusion via Tensor Ring and Subspace Decompositions. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2021; 14, pp. 8823-8837. [DOI: https://dx.doi.org/10.1109/JSTARS.2021.3108233]

34. Chen, Y.; Zeng, J.; He, W.; Zhao, X.L.; Huang, T.Z. Hyperspectral and Multispectral Image Fusion Using Factor Smoothed Tensor Ring Decomposition. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 5515417. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3114197]

35. Yokoya, N.; Yairi, T.; Iwasaki, A. Coupled Nonnegative Matrix Factorization Unmixing for Hyperspectral and Multispectral Data Fusion. IEEE Trans. Geosci. Remote Sens.; 2012; 50, pp. 528-537. [DOI: https://dx.doi.org/10.1109/TGRS.2011.2161320]

36. Yang, F.; Ping, Z.; Ma, F.; Wang, Y. Fusion of Hyperspectral and Multispectral Images With Sparse and Proximal Regularization. IEEE Access; 2019; 7, pp. 186352-186363. [DOI: https://dx.doi.org/10.1109/ACCESS.2019.2961240]

37. Simões, M.; Bioucas-Dias, J.; Almeida, L.B.; Chanussot, J. A Convex Formulation for Hyperspectral Image Superresolution via Subspace-Based Regularization. IEEE Trans. Geosci. Remote Sens.; 2015; 53, pp. 3373-3388. [DOI: https://dx.doi.org/10.1109/TGRS.2014.2375320]

38. Dian, R.; Li, S. Hyperspectral Image Super-Resolution via Subspace-Based Low Tensor Multi-Rank Regularization. IEEE Trans. Image Process.; 2019; 28, pp. 5135-5146. [DOI: https://dx.doi.org/10.1109/TIP.2019.2916734] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31107646]

39. Prévost, C.; Usevich, K.; Comon, P.; Brie, D. Hyperspectral Super-Resolution with Coupled Tucker Approximation: Recoverability and SVD-Based Algorithms. IEEE Trans. Signal Process.; 2020; 68, pp. 931-946. [DOI: https://dx.doi.org/10.1109/TSP.2020.2965305]

40. Kanatsoulis, C.I.; Fu, X.; Sidiropoulos, N.D.; Ma, W.K. Hyperspectral Super-Resolution Via Coupled Tensor Factorization: Identifiability and Algorithms. Proceedings of the 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP); Calgary, AB, Canada, 15–20 April 2018; pp. 3191-3195. [DOI: https://dx.doi.org/10.1109/ICASSP.2018.8462525]

41. Dian, R.; Fang, L.; Li, S. Hyperspectral Image Super-Resolution via Non-local Sparse Tensor Factorization. Proceedings of the 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR); London, UK, 21–26 July 2017; pp. 3862-3871. [DOI: https://dx.doi.org/10.1109/CVPR.2017.411]

42. Dian, R.; Li, S.; Fang, L. Learning a Low Tensor-Train Rank Representation for Hyperspectral Image Super-Resolution. IEEE Trans. Neural Networks Learn. Syst.; 2019; 30, pp. 2672-2683. [DOI: https://dx.doi.org/10.1109/TNNLS.2018.2885616]

43. Zhao, Q.; Zhou, G.; Xie, S.; Zhang, L.; Cichocki, A. Tensor Ring Decomposition. arXiv; 2016; arXiv: 1606.05535

44. He, W.; Chen, Y.; Yokoya, N.; Li, C.; Zhao, Q. Hyperspectral super-resolution via coupled tensor ring factorization. Pattern Recognit.; 2022; 122, 108280. [DOI: https://dx.doi.org/10.1016/j.patcog.2021.108280]

45. Zhao, Q.; Zhang, L.; Cichocki, A. Bayesian CP Factorization of Incomplete Tensors with Automatic Rank Determination. IEEE Trans. Pattern Anal. Mach. Intell.; 2015; 37, pp. 1751-1763. [DOI: https://dx.doi.org/10.1109/TPAMI.2015.2392756]

46. Li, S.; Dian, R.; Fang, L.; Bioucas-Dias, J.M. Fusing Hyperspectral and Multispectral Images via Coupled Sparse Tensor Factorization. IEEE Trans. Image Process.; 2018; 27, pp. 4118-4130. [DOI: https://dx.doi.org/10.1109/TIP.2018.2836307]

47. Zhang, H.; Zhao, X.L.; Jiang, T.X.; Ng, M.K.; Huang, T.Z. Multiscale Feature Tensor Train Rank Minimization for Multidimensional Image Recovery. IEEE Trans. Cybern.; 2022; 52, pp. 13395-13410. [DOI: https://dx.doi.org/10.1109/TCYB.2021.3108847]

48. Cao, M.; Bao, W.; Qu, K. Hyperspectral Super-Resolution Via Joint Regularization of Low-Rank Tensor Decomposition. Remote Sens.; 2021; 13, 4116. [DOI: https://dx.doi.org/10.3390/rs13204116]

49. Ghadimi, E.; Teixeira, A.; Shames, I.; Johansson, M. Optimal Parameter Selection for the Alternating Direction Method of Multipliers (ADMM): Quadratic Problems. IEEE Trans. Autom. Control; 2015; 60, pp. 644-658. [DOI: https://dx.doi.org/10.1109/TAC.2014.2354892]

50. Guo, H.; Bao, W.; Qu, K.; Ma, X.; Cao, M. Multispectral and Hyperspectral Image Fusion Based on Regularized Coupled Non-Negative Block-Term Tensor Decomposition. Remote Sens.; 2022; 14, 5306. [DOI: https://dx.doi.org/10.3390/rs14215306]

51. Ye, Q.; Amini, A.A.; Zhou, Q. Optimizing Regularized Cholesky Score for Order-Based Learning of Bayesian Networks. IEEE Trans. Pattern Anal. Mach. Intell.; 2021; 43, pp. 3555-3572. [DOI: https://dx.doi.org/10.1109/TPAMI.2020.2990820] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32340938]