Multiband Image Fusion via Regularization on a Riemannian Submanifold

Abstract

Multiband image fusion aims to generate high spatial resolution hyperspectral images by combining hyperspectral, multispectral or panchromatic images. However, fusing multiband images remains a challenge due to the identifiability and tracking of the underlying subspace across varying modalities and resolutions. In this paper, an efficient multiband image fusion model is proposed to investigate the latent structures and intrinsic physical properties of a multiband image, which is characterized by the Riemannian submanifold regularization method, nonnegativity and sum-to-one constraints. An alternating minimization scheme is proposed to recover the latent structures of the subspace via the manifold alternating direction method of multipliers (MADMM). The subproblem with Riemannian submanifold regularization is tackled by the projected Riemannian trust-region method with guaranteed convergence. The effectiveness of the proposed method is demonstrated on two multiband image fusion problems: (1) hyperspectral and panchromatic image fusion and (2) hyperspectral, multispectral and panchromatic image fusion. The experimental results confirm that our method demonstrates superior fusion performance with respect to competitive state-of-the-art fusion methods.

1. Introduction

With the development of Earth observation satellites, multiband images (e.g., hyperspectral (HS) images, multispectral (MS) images and panchromatic (PAN) images), can be obtained. These data are widely applied in various applications [1,2], such as monitoring land use [3] and ice sheets [4]. However, these images have critical trade-offs between the spatial resolution and the spectral resolution due to hardware limitations. For example, an HS image has a high spectral resolution with a reduced spatial resolution. An MS image offers a moderate spatial resolution with only a few bands. PAN images provide a much better spatial resolution.

Due to the increasing availability of optical imaging systems, they are of interest to multiband image fusion, including the fusion of HS and PAN images (HS-PAN image fusion), hyperspectral and multispectral (HS-MS) image fusion or simultaneous hyperspectral, multispectral and panchromatic (HS-MS-PAN) image fusion. Specifically, the multiband image fusion problem refers to the process of recovering a 3D data cube from two or three degraded data cubes. Thus, proper modeling of these data plays an essential role in achieving these goals [5].

There has been a remarkable effort in the community toward multiband image fusion [6,7,8,9]. For example, HS-PAN image fusion (also known as hyperspectral pan sharpening) [6,8,9] can be broadly divided into six classes: component substitution (CS), multi-resolution analysis (MRA), Bayesian, matrix factorization, deep learning and hybrid methods. For improvement of the fusion performance, jointly fusing HS, MS and PAN images has been investigated [10]. However, this method did not take into account the low-dimensional structure of multiband image, which has recently gained much interest [11,12]. More importantly, it is a common trend to identify and represent the intrinsic structure of a multiband image.

A promising method for multiband imaging is subspace representation [12,13,14,15,16], which assumes the latent subspace to be of a low rank. A first successful attempt was presented in [13,14] by exploiting the local low-rank subspace of the spectral signatures in each patch. In the work of Zhang et al. [15], a similar local geometry represented by the group spectral embedding regularization method was proposed for HS-MS image fusion. Furthermore, Kanatsoulis et al. [16] exploited the multi-dimensional structure of HS imaging and approximated the spectral images via canonical polyadic decomposition. It has been demonstrated that these methods give competitive fusion performance. However, there are several challenges to identifying and tracking the underlying subspace across varying modalities. It is difficult to assert the uniqueness of the latent subspace to preserve the significant structure in the data, which is sensitive to the existence of outliers.

To address the aforementioned issues, a Riemannian submanifold regularization method is proposed to identify and capture the local geometry of the estimated HS image by exploiting the geometry structure in an abundance matrix [14]. Specifically, this regularization method is implicitly imposed on the linear spectral mixture model (i.e., abundance matrix). The advantages of the Riemannian submanifold regularization method are that it can represent the abundance correlation and recover the subspace embedded in the HS image. These advancements motivated us to investigate the smooth structure of the underlying subspace for integrating the spatial and spectral information.

1.1. Scope and Contributions

In this work, this paper proposes an efficient multiband image fusion method with Riemannian submanifold regularization and additional constraints. The key idea is to establish a new multiband image fusion model to estimate the endmember and the abundance simultaneously. It overcomes the two problems in the following ways. First, the problem related to the identifiability of the latent subspace is addressed by the Riemannian submanifold regularization method, nonnegativity and sum-to-one constraints. Second, the subproblem of tracking the underlying subspace is reformulated as an optimization problem on a Riemannian submanifold.

The contributions can be summarized as follows:

• An efficient multiband image fusion model utilizing the Riemannian submanifold regularization method is proposed. This model is characterized by rank equality constraints with matrix manifold, nonnegativity and sum-to-one constraints. This is a new problem formulation for investigating the latent structures across varying modalities and resolutions.
• An alternating minimization scheme is proposed to recover the latent structures of the subspace using the framework of the manifold alternating direction method of multipliers. An efficient projected Riemannian trust region method with guaranteed convergence is adopted to track the latent subspace.
• The proposed method is validated in two applications: (1) hyperspectral and panchromatic image fusion and (2) the fusion of hyperspectral, multispectral and panchromatic images. The experimental results show that the proposed method is more effective than the competitive state-of-the-art fusion methods.

1.2. Related Work

1.2.1. HS-PAN Image Fusion

HS-PAN image fusion aims to perform pan sharpening with HS imaging [6]. It is a special case of HS-MS image fusion [17]. The early HS-PAN image fusion methods involve combining an MS image with a PAN one. Traditional pan sharpening methods are extended for HS-PAN image fusion. Two representative methods include CS and MRA. Some well-known CS methods are principal components analysis (PCA) [18], Brovey transform (BT) [19] and the intensity-hue-saturation (IHS) method [20]. MRA methods are characterized by a multi-resolution decomposition, such as wavelet transform methods [21,22], intensity modulation with a smoothing filter [23] and the generalized Laplacian pyramid [24]. Recently, many MS-PAN image fusion methods have been extended to fusing HS images with PAN ones. Vivone et al. [25] extended and investigated MRA- and CS-based methods for fusing HS images. A guided filter in the PCA domain was proposed [18] for sharpening HS images, and it provided good performance. Recently, deep learning methods within the framework of a convolutional neural network [8,26] (CNN) were investigated. These methods showed promising fusion performance because of the ability to extract high-level features. However, deep learning methods require extremely large datasets for training.

1.2.2. HS-MS Image Fusion

HS-MS image fusion has attracted increasing interest in recent years [7,27]. Among the existing state-of-the-art HS-MS image fusion methods, three dominant frameworks are matrix factorization (MF) approaches, tensor factorization (TF) approaches and deep learning approaches.

The key idea of MF approaches is to factorize the HS image into the spectral basis and the corresponding coefficients. The advantages of MF methods derive from the spectral signature represented by sparse representation [28,29] or a low rank [5,13,30]. The sparse representation methods model the HS image via sparse dictionary learning. Recently, Dong et al. [31] promoted the nonlocal self-similarities in an HS image via the promising alternative framework of nonlocal sparse representation. The low-rank methods regard the spectral signatures represented by a low-dimensional subspace. For example, Simões et al. [30] made use of the total variation regularization method to promote sparsity in the gradient domain. However, these methods rely on accurate estimation of the sensor parameters, such as the spectral response function.

TF approaches have been an active topic. These methods view an HS image datacube as a three-dimensional tensor, where the HS image is factored into different low-dimensional subspaces. Dian et al. [32] proposed a sparse tensor factorization-based fusion model by utilizing the properties of nonlocal self-similarity in the spatial domain and global correlation of an HS image in the spectral domain. Coupling the sparse tensor factorization method with a core tensor regularizer [11] was proposed for the fusion of HS and MS images, which provides promising fusion performance. Furthermore, the fusion framework with the spatial spectral graph-regularized low-rank tensor decomposition method was developed for HS-MS image fusion [33]. In the work of Dian et al. [12], the prior low tensor train rank of the grouped 4D tensors was exploited. However, the TF methods involve the operation of singular value decomposition (SVD), which may lead to high computational cost.

Deep learning approaches have been successfully applied to this task [34,35,36]. For example, Yang et al. [34] developed a deep two-branche convolutional neural network for image fusion, which takes the spectral correlation of HS and MS images into consideration. Huang et al. [35] proposed a deep hyperspectral image fusion network (DHIF-Net) via an iterative spatio-spectral regularization scheme. Specifically, they exploited the spatio-spectral regularization and physical imaging models simultaneously. Recently, Xie et al. [36] introduced an interpretable deep network, utilizing the low-rank knowledge of an HS image and the imaging models of HS and MS images. However, deep learning-based approaches still have some challenges, such as estimating the architecture and complexity of the deep network and the high inference cost.

Furthermore, most existing MF-based methods for HS-MS image fusion share two common limitations: (1) identifiability of the underlying subspace across modalities and (2) tracking the underlying subspace due to the non-uniqueness of SVD. Although, these methods construct the subspace representation for a pair of multiband images with varied modalities or resolutions [30], few studies established the relationship between the identifiability and the tracking of the latent subspace for HS-MS image fusion. Furthermore, fusing multiband images more than two times is still an open issue. Meanwhile, the HS image is assumed to be of a low rank [30,37], which indicates that it lies in a common subspace across modalities. This motivated us to propose a new multiband image fusion model for recovering and tracking the latent subspace. It provides an alternative way to investigate the underlying structure of these high-dimensional datasets.

This paper is organized as follows. In Section 2, some preliminaries about the manifold optimization method are presented. In Section 3, the fusion model and the notations are presented. In Section 4, the details of the proposed alternating minimization scheme are presented. In Section 5, the experimental results are illustrated. A detailed discussion is presented in Section 6. At last, the conclusion is provided in Section 7.

2. Preliminaries for Riemannian Manifold Optimization

In this section, some ingredients of Riemannian manifold optimization are provided for finding the optimal solution to the problem [38,39]. The Riemannian manifold is referred to as a manifold due to its tangent spaces endowed with a smoothly varying metric. This paper focuses on the manifold of low-rank matrices of rank r embedded in ${\mathbb{R}}^{n\times m}$ (i.e., ${\mathcal{M}}_r$), where $r<m,n$. Embedded manifolds have many inherited properties from ${\mathbb{R}}^{n\times m}$.

We focus on the optimization of functions f defined on Riemannian manfolds $({\mathcal{M}}_r,g)$:

$$\underset{x}{min}f\left(x\right)s.t.x\in {\mathcal{M}}_r,$$

where g is a Riemannian metric defined by an inner product on the tangent spaces. However, there exists a challenge in this problem.

Recently, a class of operations called retractions was introduced to deal with this problem [38]. In the context of Euclidean optimization, a gradient step update $x_t-\nabla f\left(x_t\right)$ takes the model problem outside of the manifold ${\mathcal{M}}_r$ at every step t of the optimization algorithm. Thus, it has to be pulled back onto the manifold. An ideal operation is called exponential mapping, which maps the tangent vector to a point along a geodesic curve. Unfortunately, it may be computationally expensive to calculate the geodesic curve.

2.1. Riemannian Gradient and Tangent Space

A Riemannian optimization algorithm typically conducts a line search or solves a model problem in a tangent space [38]. Therefore, the gradient and the Hessian of an objective function on a Riemannian manifold are two basic concepts.

Let it be given that $\mathcal{M}$ is a smooth submanifold of a Euclidean space, and let $x\in \mathcal{M}$. The tangent space of $\mathcal{M}$ at x, denoted by $T_x\mathcal{M}$, is a group of derivatives of all the smooth curves passing through x, and $T_x\mathcal{M}=\left\{\gamma \left(0\right)\right|\gamma \left(t\right)\}$ is a curve in $\mathcal{M}$ with $\gamma \left(0\right)=x$. The tangent space is a vector space, and its element in $T_x\mathcal{M}$ corresponds to a linear mapping from the set of smooth, real-valued functions in the neighborhood from x to $\mathbb{R}$. $T_x\mathcal{M}$ is equipped with an inner product (or metric) $g_x(·,·):T_x\mathcal{M}\times T_x\mathcal{M}\to \mathbb{R}$.

Consider a smooth function $f:\mathcal{M}\to \mathbb{R}$ defined on a Riemannian submanifold. The Riemannian gradient of f at x (i.e., $\mathrm{grad}f\left(x\right)$) is the unique tangent vector

$$g_x(\mathrm{grad}f\left(x\right),{\xi }_x)=\mathrm{D}f\left(x\right)\left[{\xi }_x\right],\forall {\xi }_x\in T_x\mathcal{M},$$

where $\mathrm{D}f\left(x\right)\left[{\xi }_x\right]$ is the directional derivative of f along the direction ${\xi }_x$, specifically when ${\mathcal{M}}_r$ is a Riemannian submanifold embedded in ${\mathbb{R}}^{m\times n}$. A useful property of embedded manifolds is that their Riemannian gradients are defined as the orthogonal projection on the tangent space of the gradient of f; in other words, we have

$$\mathrm{grad}f\left(x\right)={\mathrm{P}}_{T_x\mathcal{M}}(\nabla f\left(x\right)),{\xi }_x\in T_x\mathcal{M},$$

where $\nabla f\left(x\right)$ is the Euclidean gradient of f at x. It should be noted that ${\mathrm{P}}_{T_x\mathcal{M}}$ is viewed as a retraction [38,40].

The Hessian of f at x (i.e., $\mathrm{Hess}f\left(x\right)$) is a mapping from $T_x\mathcal{M}$ to $T_x\mathcal{M}$. Furthermore, when $\mathcal{M}$ is a Riemannian submanifold of a Euclidean space ${\mathbb{R}}^{m\times n}$, the Riemannian Hessian of f is expressed as follows:

$$\mathrm{Hess}f\left(x\right)\left[{\xi }_x\right]={\mathrm{P}}_{T_x\mathcal{M}}\left(\mathrm{D}\mathrm{grad}f\left(x\right)\left[{\xi }_x\right]\right),{\xi }_x\in T_x\mathcal{M},$$

where $\mathrm{D}\mathrm{grad}f\left(x\right)\left[{\xi }_x\right]$ denotes a classical directional derivative of $\mathrm{grad}f\left(x\right)$ along the direction ${\xi }_x$. Similarly, an orthogonal projection is performed.

2.2. Retractions

Retractions can find a point on the manifold that is in the direction of the gradient. They can approximate the geodesic gradient flow on the manifold [38]. In other words, retractions allow us to move on a manifold (i.e., moving in the direction of a tangent vector) while staying on the manifold. It provides an alternative to exponential mapping when designing optimization algorithms on manifolds. For the manifold ${\mathcal{M}}_r$, a retraction is defined as follows:

$$x_{t+1}=R_{x_t}(x_t+\xi )=\mathrm{P}(x_t+\alpha ).$$

Mathematically, the retractions are defined as a mapping, where $R_x$ is a smooth mapping from $T_x\mathcal{M}$ to $\mathcal{M}$ such that $R_x\left(0\right)=x$. The corresponding differential of retraction at zero is an identity mapping (i.e., $\frac{d}{dt}{R}_x{\left(t\xi \right)|}_{t=0}=\xi ,\forall \xi \in T_x\mathcal{M}$), which is referred to as the local rigidity condition [38].

3. Problem Formulation

3.1. Degradation Model for Multiband Imaging

Let $P\in {\mathbb{R}}^{n_s\times n}$ denote the matricization of the latent multiband image. We denote $n_s$ to be the number of spectral bands, while $n=n_x\times n_y$ is the vectorization of the HS image at the band $n_x$. The goal is to recover an HS image of the high spatial resolution P from K observations (i.e., $Y_k\in {\mathbb{R}}^{n_{s_k}\times n_k},k=1,\cdots ,K$, where $n_{s_k}$ and $n_k$ represent the number of bands and pixels, respectively). Moreover, we assume that $n_{s_k}\le n_s$, and $n_k=n/d_k^2$. $d_k$ denotes the scale factor.

Usually, the responses of the imaging sensors are treated as a number of linear transformations. The corresponding spectral and spatial degradations are denoted as the operators. The input multiband images are assumed to be geometrically registered. The degradation model can be expressed as

$$Y_k=R_{k}P{B}_k{S}_k+N_k,$$

(1)

where $R_k\in {\mathbb{R}}^{n_k\times n_s}$ is the spectral response of the imaging sensor, $B_k$ denotes the point spread function of the kth imaging sensor, $S_k$ represents a downsampling operator for the spatial dimensions subject to $S_k^T{S}_k=I_{n_k}$ and $N_k$ indicates a spectral independent perturbation matrix. In particular, $Y_k$ involves some preprocessing steps [41], such as radiometric calibration and geometric correction.

One commonly used approach for the spectral representation of an observed multiband image is the linear mixture model [42] (LMM), which reveals the relationship between mixtures and endmembers. This model assumes that every pixel of an HS image can be decomposed into a convex combination of a small number $m\ll n_s$ of endmembers [43]. This indicates that the number of endmembers is usually sparse (i.e., the number of nonzero entries is very small compared with the size of the abundances). The latent multiband image is the product of an endmember matrix E and the abundance matrix X. Thus, the LMM model can be expressed as follows:

$$P=EX+N,$$

(2)

where $E\in {\mathbb{R}}^{n_s\times m}$ represents the spectral signatures matrix with respect to the endmembers, $X\in {\mathbb{R}}^{m\times n}$ is the matrix of endmember abundances stored in columns and $N\in {\mathbb{R}}^{n_s\times n}$ denotes the measurement noise or model error [42]. It can be seen that each row of the matrix P consists of all the pixels in a given spectral band.

By substituting Equation (2) into Equation (1), we can obtain a multiple multiband image formation model:

$$Y_k=R_{k}EX{B}_k{S}_k+{\widehat{N}}_k,$$

(3)

where ${\widehat{N}}_k$ denotes the kth multiband image’s additive perturbation matrix, which can be is expressed as follows:

$${\widehat{N}}_k=N_k+R_{k}N{B}_k{S}_k.$$

Given the observation $Y_k$ and its corresponding endmember matrix E, the measurement P can be recovered from the the abundance matrix X. However, additional constraints on E or X should be used to improve the quality of the solutions. Thus, additional information is introduced (i.e., the Riemannian submanifold regularization method, nonnegativity and sum-to-one constraints).

3.2. Proposed Fusion Model

For the multi-source information modeling, naturally, the prior information about X is required. To reconstruct P, the proposed fusion model can be formulated as an optimization problem:

$$\begin{array}{c}\underset{X\in {\mathcal{M}}_r}{min}\frac{1}{2}{\mathsf{\Sigma }}_{k=1}^K\parallel Y_k-R_{k}EX{B}_k{S}_k{\parallel }_F^2+\frac{\alpha }{2}{\parallel X\parallel }_F^2+{\mathcal{I}}_{\Omega }\left(X\right),\hfill \end{array}$$

(4)

where ${\mathcal{M}}_r=\{X\in \mathcal{M}|\mathrm{rank}\left(X\right)=r\}$ represents the rank equality constraints with the matrix manifold, the symbol $\mathrm{rank}(·)$ denotes the rank of the matrix, r is a positive integer, the term $\frac{\alpha }{2}{\parallel X\parallel }_F^2$ is the Riemannian submanifold regularization (i.e., the constraint set ${\mathcal{M}}_r$), $\alpha$ is a positive regularization parameter, ${\mathcal{I}}_{\Omega }(·)$ denotes an indicator function and the expression $\Omega =\left\{X\right|X\ge 0,1_m^{T}X=1_n^T\}$ indicates that all the elements of X are greater than or equal to zero and sum to one. The proposed framework for multiple multiband image fusion is presented in Figure 1.

3.3. Riemannian Submanifold Regularization

We observe that multiband image lies on a low-dimensional manifold $\mathcal{M}$. Specifically, we take a different approach which exploits the problem structure of rank equality constraints with a matrix manifold (i.e., the latent structures of the abundance matrix). More specifically, the structural information of the abundance matrix is captured and represented by the geometric structure of the matrix manifold ${\mathcal{M}}_r$. This paper proposes a regularization method to extract a common subspace from multiple observations indirectly, which provides more flexibility to handling data with nonlinear structures. It can be seen that this regularization process is implemented through rank equality constraints with a matrix manifold.

Definition 1 

(Riemannian submanifold regularization). Let $X\in {\mathcal{M}}_r$ be a smooth manifold [38] with the constraint set of matrices of rank being r at most, where $r<min(m,n)$ is a positive integer. To stabilize the fusion model in Equation (3), additional a priori knowledge about the abundance matrix X is defined as follows:

$$\frac{\alpha }{2}{\parallel X\parallel }_F^2,s.t.X\in {\mathcal{M}}_r=\{X\in {\mathbb{R}}^{m\times n}|\mathit{rank}\left(X\right)=r\},$$

(5)

where ${\mathcal{M}}_r$ is viewed as a submanifold of a dimension $(m+n-r)r$ embedded in the Euclidean space ${\mathbb{R}}^{m\times n}$. Here, the constraint $\mathit{rank}\left(X\right)=r$ denotes the rank of matrix X as equal to r, which can be viewed as a rank equality constraint with a matrix manifold. The Tikhonov regularization term [44] (i.e., $\frac{\alpha }{2}{\parallel X\parallel }_F^2$) ensures the available solution.

However, it is challenging to deal with the rank equality constraints with a matrix manifold, where the underlying subspace is spanned by endmembers indirectly. The corresponding subspace of the target image P lies in a low-dimensional manifold. Moreover, these constraints are nonconvex, which may be numerically expensive.

In this paper, the manifold geometry [45] of Riemannian submanifold regularization (i.e., the set of rank equality constraints) is utilized to address this issue. Thus, both the constraint search space and the geometry of the smooth embedded submanifold should be considered. Extensive studies of the optimization method on a Riemannian manifold led to the retraction-based framework using the differential geometry structure of the underlying manifold.

3.4. Nonnegativity and Sum-to-One Constraints

In addition to the geometric structure given by Riemannian submanifold regularization, it is reasonable to assume that the coefficients should satisfy the following constraints [42,43]:

$$X\ge 0,1_m^{T}X=1_n^T,$$

where $1_m^T$ means an $m\times 1$ vector with all ones. The combined constraints preserve the inherent characteristics of the solutions. In other words, it makes a solution more feasible. The nonnegativity and sum-to-one constraints can be rewritten as follows:

$$\Omega =\left\{X\right|X\ge 0,1_m^{T}X=1_n^T\}.$$

After introducing the indicator function $\mathcal{I}\left(X\right)$, we can obtain the following formulation:

$${\mathcal{I}}_{\Omega }\left(X\right)=\left\{\begin{array}{cc}0\hfill & X\in \Omega \hfill \\ +\infty \hfill & X\notin \Omega .\hfill \end{array}\right.$$

4. Proposed Method

4.1. Alternating Minimization Scheme

Given the fusion problem in Equation (4), we propose a numerical solution using the framework of the manifold alternating direction method of multipliers. The constrained search space is equipped with the Riemannian structure of the submanifold, and the nonnegativity and sum-to-one constraints are considered. A key challenge of this problem is the geometry of the underlying manifold.

After applying the variable splitting strategy, we solve the problem in an alternating minimization way. Several auxiliary variables ${\wedge }_{1,\cdots ,k},W,Z$ are introduced. We have the following the minimization problem:

$$\begin{array}{c}\underset{X\in {\mathcal{M}}_r,{\wedge }_1,\cdots ,{\wedge }_k,U,V}{min}\frac{1}{2}{\mathsf{\Sigma }}_{k=1}^K{\parallel Y_k-R_{k}E{\wedge }_k{S}_k\parallel }_F^2\hfill \\ +\frac{\alpha }{2}{\parallel W\parallel }_F^2+{\mathcal{I}}_{\Omega }\left(Z\right),\hfill \\ s.t.{\wedge }_k=X{B}_k,W=X,Z=X.\hfill \end{array}$$

(6)

Then, the augmented Lagrangian formulation can be expressed as

$$\begin{array}{c}\mathcal{L}(X,{\wedge }_1,\cdots ,{\wedge }_k,W,Z,F_1,\cdots ,F_K,G,H)=\hfill \\ \frac{1}{2}{\mathsf{\Sigma }}_{k=1}^K\parallel Y_k-R_{k}E{\wedge }_k{S}_k{\parallel }_F^2+\frac{\alpha }{2}{\parallel W\parallel }_F^2+{\mathcal{I}}_{\Omega }\left(Z\right)\hfill \\ +\frac{\mu }{2}{\mathsf{\Sigma }}_{k=1}^K\parallel X{B}_k-{\wedge }_k-F_k{\parallel }_F^2+\frac{\mu }{2}{\parallel X-W-G\parallel }_F^2\hfill \\ +\frac{\mu }{2}{\parallel X-Z-H\parallel }_F^2,\hfill \end{array}$$

(7)

where $F_k\in {\mathbb{R}}^{m\times n},k=1,\cdots ,K,G\in {\mathcal{M}}_r$, $H\in {\mathbb{R}}^{m\times n}$ are Lagrangian multipliers and $\mu \ge 0$ is a penalty parameter. Then, the problem has the following formulations:

$$\begin{array}{c}{X}^n=arg\underset{X}{min}\mathcal{L}(X,{\wedge }_1^{n-1},\cdots ,{\wedge }_k^{n-1},W^{n-1},\hfill \\ Z^{n-1},F_1^{n-1},\cdots ,F_K^{n-1},G^{n-1},H^{n-1}),\hfill \end{array}$$

(8)

$$\begin{array}{c}\{{\wedge }_{1,\cdots ,K}^n,U^n,V^n\}=arg\underset{{\wedge }_{1,\cdots ,K},U,V}{min}\mathcal{L}(X^n,{\wedge }_1,\hfill \\ \cdots ,{\wedge }_k,W,Z,F_1^{n-1},\cdots ,F_K^{n-1},G^{n-1},H^{n-1}),\hfill \end{array}$$

(9)

$$\begin{array}{c}{F}_k^n=F_k^{n-1}-\left(X^n{B}_k-{\wedge }_k^n\right),k=1,\cdots ,K,\hfill \end{array}$$

(10)

$$\begin{array}{c}{G}^n=G^{n-1}-\left(X^n-W^n\right),\hfill \end{array}$$

(11)

$$\begin{array}{c}{H}^n=H^{n-1}-\left(X^n-Z^n\right).\hfill \end{array}$$

(12)

Thus, the problem (i.e., Equation (7)) breaks down into three subproblems. The details are summarized in Algorithm 1. The computational complexity of the proposed method is $max(\mathcal{O}\left(r^3\right),\mathcal{O}(mnlog\left(mn\right)))$, where $r<min(m,n)$:

$$\begin{array}{c}{\wedge }_k^n=arg\underset{{\wedge }_k}{min}\frac{1}{2}{\parallel Y_k-R_{k}E{\wedge }_k{S}_k\parallel }_F^2+\hfill \\ \frac{\mu }{2}{\parallel X^n{B}_k-{\wedge }_k-F_k^{n-1}\parallel }_F^2,k=1,\cdots ,K,\hfill \end{array}$$

(13)

$$W^n=arg\underset{W}{min}\frac{\alpha }{2}{\parallel W\parallel }_F^2+\frac{\mu }{2}{\parallel X^n-W-G^{n-1}\parallel }_F^2,$$

(14)

$$Z^n=arg\underset{Z}{min}{\mathcal{I}}_{\Omega }\left(Z\right)+\frac{\mu }{2}{\parallel X^n-Z-H^{n-1}\parallel }_F^2.$$

(15)

Note that the subproblem in Equation (14) can be cast as an optimization problem on a Riemannian submanifold. This problem can be solved by the retraction-based optimization method on the manifold. In addition, the solution to the subproblem Z has a closed-form expression. The value of $\mu$ is updated at each iteration. The optimization of $\mathcal{L}$ with respect to each of the variables $X,\wedge$ is described in Appendix A. For the variables W and Z, the solutions are described in detail in Appendix B.

Algorithm 1 Optimization procedures for the problem in Equation (6)

1: Initializing:$E=\mathrm{VCA}\left(Y_l\right)$;// estimate E by the method in [46] 2: $X^0=\mathrm{Interpolate}(Y_l,E)$; // interpolate the outcome of SUnSAL [47] 3: for k = 1:K do 4:    ${\wedge }_k^0=X^0;F_k^n=0_{m\times n};$ 5:    $W^0=X^0;Z^0=X^0;$ 6:    $G^0=0_{m\times n};H^0=0_{m\times n};$ 7: end for 8: for n = 1:N do 9:     Computing $X^n$; // provided in Equation (A2) 10:    for k = 1:K do 11:      Computing ${\wedge }_k^n$; // provided in Equation (A5) 12:      Computing $X^n-G^{n-1}$; 13:    end for 14:    Computing $W^n$; // refer to Algorithm A1 15:    Computing $Z^n$; // refer to Equation (A10) and [48]. 16:    for k = 1:K do 17:      $F_k^n=F_k^{n-1}-(X^n{B}_k-{\wedge }_k^n)$ 18:      $G^n=G^{n-1}-(X^n-W^n)$ 19:      $H^n=H^{n-1}-(X^n-Z^n)$ 20:    end for 21: end for 22: return $X=E{X}^n$

4.2. Convergence Analysis

This subsection studies the behavior of the convergence of the proposed method. Note that the convergence behavior and results of the ADMM were widely studied in [49]. Thus, we analyzed the convergence of Algorithm A1 for the subproblem with Riemannian submanifold regularization, which makes use of the retraction-based manifold optimization theory [38]:

Theorem 1.

If there exists a global minimizer $W^{*}$ to the problem in Equation (A6), then the first-order necessary optimality condition holds:

$$\langle \zeta ,(1+\frac{\alpha }{\mu })W^{*}-(X^n-G^{n-1})\rangle \ge 0,\forall \zeta \in T_W{\mathcal{M}}_r,$$

where $T_W{\mathcal{M}}_r$ denotes the tangent cone of the constraint set ${\mathcal{M}}_r$ at $W^{*}$.

Proof. 

Suppose that $W^n\in {\mathcal{M}}_r$ is a sequence generated by Equation (A6); in other words, we have

$$\underset{n\to \infty }{lim}{f}_W^n\left(W^n\right)=\underset{W\in {\mathcal{M}}_r^{*}}{inf}{f}^n\left(W\right).$$

Since $f^n$ is bounded and coercive, it gives

$$f^n\left(W\right)\to \infty ,\mathrm{if}\parallel W\parallel \to \infty ,$$

Thus, this indicates that the sequence $W^n$ is uniformly bounded. Furthermore, all points belong to the set W, which is compact. Then, it follows that the sequence is bounded.

It is of interest to note that $\left\{W^n\right\}$ has an accumulation point $W^{*}$. As the feasible set ${\mathcal{M}}_r$ is closed, and $f^n:{\mathcal{M}}_r\to \mathbb{R}$ is continuous, therefore, we conclude that $W^{*}$ has a global minimizer.

Then, the first-order necessary optimality condition of the global minimizer is presented directly:

$$\langle \zeta ,(1+\frac{\alpha }{\mu })W^{*}-(X^n-G^{n-1})\rangle \ge 0,\forall \zeta \in T_W{\mathcal{M}}_r,$$

It is noticeable that $\mathrm{rank}\left(W^{*}\right)=r$. Thus, $T_W{\mathcal{M}}_r$ turns into a linear subspace in ${\mathbb{R}}^{m\times n}$ (i.e., the tangent space at $W^{*}$). The inequality condition can be

$$P_{T_{W^{*}}{\mathcal{M}}_r}((1+\frac{\alpha }{\mu })W^{*}-(X^n-G^{n-1}))=0,$$

where $P_{T_{W^{*}}{\mathcal{M}}_r}$ denotes the orthogonal projection onto the linear subspace $T_{W^{*}}{\mathcal{M}}_r$.    □

Theorem 2.

Assume that κ is a positive constant and $\left\{{ϵ}^n\right\}$ is a sequence of nonnegative scalars. If $\left\{W^n\right\}\subset {\mathcal{M}}_r$ is a sequence generated by Algorithm A1, then we have

$$f^n\left(W^{n+1}\right)\le f^n\left(W^n\right)-\kappa {\parallel W^{n+1}-W^n\parallel }^2,$$

(16)

$$\langle \xi ,(1+\frac{\alpha }{\mu })W^{n+1}-(X^n-G^{n-1})\rangle \ge -{ϵ}^k\parallel \xi \parallel ,\forall \xi \in T_W{\mathcal{M}}_r,$$

(17)

If $\left\{W^{n^s}\right\}$ is a convergent subsequence of $\left\{W^n\right\}$ with the limit point $W^{*}\in {\mathcal{M}}_r$, then $\mathit{rank}\left(W^{*}\right)=r$, and ${lim}_{s\to \infty }{ϵ}^{n^s}=0$. Finally, $W^{*}$ satisfies the first-order optimality condition for the problem in Equation (A6).

Proof. 

It can be seen that $f\left(W^{n+1}\right)\le f\left(W^n\right),\forall n$. Since f is bounded, we have

$$\underset{n\to \infty }{lim}f\left(W^{n+1}\right)-f\left(W^n\right)=\underset{n\to \infty }{lim}f\left(W^{n+1}\right)-f\left(W^{n+1}\right)=0,$$

where the conditions in Equations (16) and (17) indicate that ${lim}_{n\to \infty }\parallel W^n-W^{n+1}\parallel =0$.

If $\left\{W^{n^s}\right\}$ is a subsequence that converges to $\left\{W^{*}\right\}$ with a rank r, then we have

$$\mathrm{rank}\left(W^{n^s}\right)=r,\forall s.$$

Since ${\mathcal{M}}_r$ is a smooth manifold in a neighborhood of $W^{*}$, the condition in Equation (17) indicates that

$$\parallel P_{T_W{\mathcal{M}}_r}((1+\frac{\alpha }{\mu })W^{n+1}-(X^n-G^{n-1}))\parallel \le {ϵ}^{n^s},$$

for all sufficiently large values of s. The continuity of the smooth mappings $\left(W\right)\in {\mathcal{M}}_r\rightarrowtail P_{T_W{\mathcal{M}}_r}$ holds true by passing $s\to \infty$.    □

For studying the convergence of $W^{n+1}$, we have the following direct lemma (Lemma 1). In particular, it can be concluded that

$$\underset{n\to \infty }{lim}\parallel P_{T_W{\mathcal{M}}_r}((1+\frac{\alpha }{\mu })W^{n+1}-(X^n-G^{n-1}))\parallel =0.$$

Lemma 1. 

Let $\left\{W^n\right\}\subset {\mathcal{M}}_r$ be a sequence generated by Algorithm A2. Then, the following statements hold true:

$$\underset{n\to \infty }{lim}\parallel W^n-W^{n+1}\parallel =0.$$

Proof. 

It can be seen that the statement is the first part of the proof for Theorem 2.    □

5. Performance Evaluation

5.1. Experimental Settings

We conducted two experiments to assess the performance of the proposed method: (1) hyperspectral and panchromatic image fusion and (2) hyperspectral, multispectral and panchromatic image fusion. All experiments were performed on a laptop with Windows 10 64 bit and implemented and tested in MATLAB 2019b.

For HS-PAN image fusion, the proposed method was compared with 14 related HS-PAN image fusion methods, including the smoothing filter-based intensity modulation (SFIM) method [23], the method of a modulation transfer function-generalized Laplacian pyramid with high-pass modulation (MT) [50], the MT with high-pass modulation (MT-HPM) method [25], the Gram–Schmidt (GS) spectral sharpening method [51], the GS adaptive (GSA) method [52], PCA [53], PCA with a guided filter (PCA-GF) [18], robust fast (RF) fusion of multiband images [54], coupled nonnegative matrix factorization [5] (CF), the hyperspectral image super-resolution method with subspace regularization [30] (HY), MT with regression [55] (MT-R), MT with full-scale regression [56] (MT-FSR), MT with dual-scale regression [57] (MT-DSR) and the method of context-aware detail injection fidelity (CDIF) with adaptive coefficient estimation [58].

For HS-MS-PAN image fusion, the proposed method was compared with some state-of-art methods [5,30,54,59,60]. These methods contain two procedures. First, two popular methods for pan sharpening were used separately (i.e., the MT method and band-dependent spatial detail (BD)) [59]. Second, RF-, CF- and HY-based methods were used for HS-PAN image fusion. In addition, the endmembers presented in the HS image were extracted with a vertex component analysis (VCA) algorithm [46]. Meanwhile, sparse unmixing with the variable splitting and augmented Lagrangian (SUnSAL) method was utilized to extract the abundance matrix [47].

5.2. Datasets and Quality Measures

Four real-world datasets were used [61]: the Botswana, Washington DC Mall, Indian Pines and Kennedy Space Center datasets. The tests in our experiments were performed in a semi-synthetic way. The intensities of these images were arranged from 0 to 1. The observed three input images (PAN, MS and HS images) from the referred HS image (ground truth) were generated according to Wald’s protocol [62]. Gaussian white noise with zero mean was added to each band of the three images such that the band-specific signal-to-noise ratio (SNR) was 30 dB for the MS and HS images and 35 dB for the PAN image. The bands related to the noisy measurements were removed from the reference image. A brief description of the datasets is provided:

• Botswana dataset: The HS image was collected by a Hyperion sensor over Okavango Delta, Botswana in 2001–2004. The number of bands in our experiment was 145. The spatial resolution is $400\times 240$. The observed scene contains the land cover type information.
• Indian Pines dataset: The imaging sensor is the airborne visible infrared image spectrometer (AVIRS) airborne hyperspectral instrument. Images were captured over northern and western Indiana in the USA. The number of bands was 200. The dimensions of the HS images are $400\times 400$. The scenery of this dataset includes housing, built structures, and forests.
• Kennedy Space Center dataset: This dataset was captured at Kennedy Space Center in Florida, United States by an AVIRIS. This dataset comprises 176 bands with an image size of $500\times 400$. The content of the HS image contains various land cover types.
• Washington DC Mall dataset: This was collected by the hyperspectral digital image collection experiment (HYDICE) over the Washington DC Mall in the United States. A portion of the original data was used. The resolution of the HS image is $400\times 300$. The number of spectral bands was 191.

Four widely used quality measures [27] were used in our experiments. These measures included the relative dimensionless global error in synthesis (ERGAS), spectral angle mapper (SAM) and universal image quality index (UIQI).

5.3. Results

5.3.1. Hyperspectral and Panchromatic Image Fusion

Given an HS image, the goal of HS-PAN image fusion is to fuse a PAN image with the same spectral information while keeping some innovative content. All the methods were applied to the three datasets (i.e., Botswana, Indian Pines and Washington DC Mall). The blurring kernel was Gaussian with a size of $3\times 3$ and standard deviation of 0.8. The scaling factor was four.

Table 1. Numerical results of hyperspectral and panchromatic image fusion for three datasets.

Datasets	Botswana			Indian Pines			Washington DC Mall
Fusion Methods	ERGAS↓	SAM↓	UIQI (%)↑	ERGAS↓	SAM↓	UIQI (%)↑	ERGAS↓	SAM↓	UIQI (%)↑
SFIM	24.9670	2.6758	86.49	0.9583	1.2023	87.81	3.9752	2.6530	94.11
MT	3.0794	2.5482	90.02	0.9249	1.1800	88.83	3.3584	2.8864	95.76
MT-HPM	47.0230	2.6108	87.52	0.9722	1.2086	88.48	3.4854	2.5888	95.51
GS	2.8263	2.4166	91.25	1.3886	1.4620	81.38	5.7251	5.2996	84.39
GSA	3.1798	2.5806	90.06	0.8252	1.0897	89.76	3.5158	2.9546	95.87
PCA	2.9592	2.5074	90.42	1.8095	1.9336	76.61	3.8862	3.4006	91.50
PCA-GF	3.1833	3.1687	77.98	1.1607	1.6353	79.49	5.6117	2.8216	79.55
RF	1.8364	2.4047	94.65	1.0243	0.9509	91.41	3.5888	3.6285	93.21
CF	1.2706	2.5176	93.89	0.9441	1.7088	83.50	1.8495	3.5331	96.14
HY	1.7612	2.1741	91.97	1.0495	0.9555	90.94	2.7327	3.0719	96.33
MT-R	2.7732	3.1392	91.81	1.5811	1.4481	90.24	4.7134	2.1903	92.02
MT-FSR	2.7487	3.1581	91.90	1.5811	1.4594	90.21	4.7096	2.4774	91.85
MT-DSR	2.7728	3.1394	91.80	1.5822	1.4480	90.19	4.7108	2.1918	92.04
CDIF	2.2615	2.5503	93.15	1.2430	1.0518	93.88	4.0160	1.8530	95.85
Proposed Method	1.4321	1.9561	95.64	0.8025	0.9487	95.67	1.4122	2.3475	98.81

gives a quantitative evaluation of all the competing methods. The best fusion results are marked in bold. We can see that the proposed method delivered better ERGAS, SAM and UIQI scores than its competitors in most cases. Meanwhile, the SFIM- and MT-HPM-based methods introduced significant spectral distortions. Remarkably, the proposed method performed relatively satisfactory on all three datasets, demonstrating the effectiveness of the proposed method.

The visual comparison results are illustrated in Figure 2. The ground truth image is provided in Figure 2a. The visual results of the referred fusion methods are demonstrated. The selected part is enlarged and provided in the bottom right corner. The PCA-GF method seems to be insufficient because the fused image lost some local details. The visual result of the proposed method is presented in Figure 2p. As can be seen from Figure 2, the proposed method gave better fused images than the other methods in terms of smooth areas and texture regions. In other words, the proposed method provided the fused results closest to the ground truth in Figure 2a.

5.3.2. Hyperspectral, Multispectral and Panchromatic Image Fusion

We conducted the experiments to jointly fuse HS, MS and PAN images. In our experiments, these images were fused in four ways: PAN + HS, (MS + HS) + PAN, (PAN + MS) + HS and PAN + MS + HS. We compared the proposed method with two state-of-art fusion methods: the RF method [54], and HY method [30]. The experiments were performed under the same experimental setting with 100 iterations.

The quantitative results of all the referred methods are presented in

Table 2. Quantitative evaluation of fusing hyperspectral, multispectral and panchromatic images on four datasets.

Images	Fusion Methods	Botswana				Indian Pines
		ERGAS↓	SAM↓	UIQI (%)↑	Time (s)↓	ERGAS↓	SAM↓	UIQI (%)↑	Time (s)↓
PAN + HS	HY	1.8450	2.4035	94.49	48.34	0.8234	1.0494	87.17	86.45
	RF	1.8364	2.4047	94.65	49.18	0.8191	1.0529	87.12	90.37
(PAN + MS) + HS	BD + HY	3.0323	3.5445	93.61	49.17	0.7491	1.0801	87.16	86.26
	BD + RF	3.0192	3.5324	93.69	51.54	0.7434	1.0846	87.20	86.40
	MT + HY	2.9706	3.5923	93.88	51.83	0.9252	1.1439	84.86	84.92
	MT + RF	2.9573	3.5786	93.95	51.06	0.9170	1.1519	84.89	86.25
	CF	1.8310	2.7891	89.96	38.93	0.9213	1.7485	86.83	73.48
PAN + (MS + HS)	HY + HY	2.0442	2.1535	93.52	49.51	0.7789	0.9373	85.24	86.71
	RF + RF	2.9769	3.6799	93.57	51.77	0.8110	0.9481	84.66	90.99
PAN + MS + HS	Our Method	1.6317	1.6408	98.52	38.01	0.5812	0.8850	95.18	65.51
Images	Fusion Methods	Washington DC Mall				Kennedy Space Center
		ERGAS↓	SAM↓	UIQI (%)↑	Time (s)↓	ERGAS↓	SAM↓	UIQI (%)↑	Time (s)↓
PAN + HS	HY	3.9132	4.6076	92.63	83.29	3.9600	4.5045	94.55	111.47
	RF	3.9207	4.6075	92.64	84.08	3.8971	4.2528	94.94	127.33
(PAN + MS) + HS	BD + HY	4.0388	4.7666	92.32	83.91	3.6127	4.7748	95.42	110.40
	BD + RF	4.1409	4.7867	92.06	87.09	3.4127	4.5649	95.94	110.49
	MT + HY	4.2401	4.8087	91.63	83.42	8.3050	5.3050	86.94	128.82
	MT + RF	4.3543	4.8267	91.29	84.47	8.2069	5.2489	87.24	121.36
	CF	4.0730	5.2249	82.12	57.11	3.6472	5.3522	86.37	88.72
PAN + (MS + HS)	HY + HY	3.0955	3.6234	95.28	94.44	2.8052	5.6430	96.67	111.39
	RF + RF	3.6510	3.8968	94.26	88.87	3.6730	6.1409	95.80	116.21
PAN + MS + HS	Our Method	2.5503	2.8112	97.77	62.39	2.5044	4.2221	97.23	82.49

. The best results are indicated in bold. The referenced fusion methods were combined for fusing HS, MS and PAN images. Note that the proposed approach outperformed the other methods in almost all the quality measures. Table 2 also presents the computational time of the fusion methods across the complete datasets. We can see that the proposed method generally demonstrated competitive performance. These results suggest that our method could recover the detailed structures of the underlying HS image well.

The visual fusion results are illustrated in Figure 3 and Figure 4. Examples of two datasets (i.e., the Indian Pines and Washington DC Mall datasets) are shown separately in Figure 3a and Figure 4a, respectively. The visual fusion results of our method are presented individually in Figure 3b and Figure 4b. The reconstructed images of the HY method are presented in Figure 3c and Figure 4c one by one. In Figure 3d and Figure 4d, the per pixel residuals of these two datasets are presented visually, where all the pixels of each dataset were sorted by magnitude and plotted. We can observe that our method was better than the other methods in terms of the normalized root mean square error (NRMSE).

The reconstruction errors of the referred methods with respect to the original data are illustrated in Figure 3e and Figure 3f for the Indian Pines dataset and Figure 4e and Figure 4f for the Washington DC Mall dataset, respectively. These figures indicate that our method can achieve better fusion performance compared with the competing methods, which supports the effectiveness of the proposed method.

We have shown that the underlying subspace was captured and tracked by the proposed method effectively. The numerical experiments show that the proposed fusion model with hybrid constraints can reconstruct more local details in comparison with the referred fusion methods. The outcomes of the proposed method (e.g., the quality measures in Table 1 and Table 2 and visual comparisons from Figure 2, Figure 3 and Figure 4 provide evidence that the involvement of a Riemannian submanifold regularization method led to the improvement in fusion performance. The fusion experiments indicate that the use of a modeling strategy and the alternating minimization scheme is critical for the success of the proposed method.

6. Discussion

We showed that the proposed fusion model and its alternating minimization scheme can track the underlying subspace effectively. This difference might be related to the method used. It is possible that the geometric representation of the latent subspace across modalities provides an efficient way to describe the geometry of the subspace. One limitation of the proposed method is that the input multiband images should be registered spatially.

In summary, the proposed method with the Riemannian submanifold regularization method is useful for estimating high-fidelity hyperspectral images from multiband images. Although the local structure of the abundance matrix was selectively exploited, the geometry of the constraint manifold provided a more intuitive understanding of the geometric structure of the underlying subspace. The proposed framework provides new insights into existing multiband image fusion methods.

7. Conclusions

In this paper, an efficient multiband image fusion method was proposed to obtain a high-fidelity fused image using Riemannian submanifold regularization, nonnegativity and sum-to-one constraints. Furthermore, the identifiability and tracking of the underlying subspace was completed by an alternating minimization scheme. To exploit the internal structures of the abundance matrix, the constrained search subspace was reformulated as an optimization problem on a smooth Riemannian submanifold. An efficient projected Riemannian trust region method was developed. The experiments for two multiband image fusion problems showed that the proposed method outperformed the state-of-the-art competing fusion methods. It should be noted that there also remains room for further study. For example, some developed image fusion models account for spectral variability.