Fusing Multiple Multiband Images

Abstract

High-resolution hyperspectral images are in great demand but hard to acquire due to several existing fundamental and technical limitations. A practical way around this is to fuse multiple multiband images of the same scene with complementary spatial and spectral resolutions. We propose an algorithm for fusing an arbitrary number of coregistered multiband, i.e., panchromatic, multispectral, or hyperspectral, images through estimating the endmember and their abundances in the fused image. To this end, we use the forward observation and linear mixture models and formulate an appropriate maximum-likelihood estimation problem. Then, we regularize the problem via a vector total-variation penalty and the non-negativity/sum-to-one constraints on the endmember abundances and solve it using the alternating direction method of multipliers. The regularization facilitates exploiting the prior knowledge that natural images are mostly composed of piecewise smooth regions with limited abrupt changes, i.e., edges, as well as coping with potential ill-posedness of the fusion problem. Experiments with multiband images constructed from real-world hyperspectral images reveal the superior performance of the proposed algorithm in comparison with the state-of-the-art algorithms, which need to be used in tandem to fuse more than two multiband images.

1. Introduction

The wealth of spectroscopic information provided by hyperspectral images containing hundreds or even thousands of contiguous bands can immensely benefit many remote sensing and computer vision applications, such as source/target detection [1,2,3], object recognition [4], change/anomaly detection [5,6], material classification [7], and spectral unmixing [8,9], commonly encountered in environmental monitoring, resource location, weather or natural disaster forecasting, etc. Therefore, finely-resolved hyperspectral images are in great demand [10,11,12,13,14]. However, limitations in light intensity as well as efficiency of the current sensors impose an inexorable trade-off between the spatial resolution, spectral sensitivity, and the signal-to-noise ratio (SNR) of existing spectral imagers [15]. As a results, typical spectral imaging systems can capture multiband images of high spatial resolution at a small number of spectral bands or multiband images of high spectral resolution with a reduced spatial resolution. For example, imaging devices onboard Pleiades or IKONOS satellites [16] provide single-band panchromatic images with spatial resolutions of less than a meter and multispectral images with a few bands and spatial resolutions of a few meters while NASA’s airborne visible/infrared imaging spectrometer (AVIRIS) [17] provides hyperspectral images with more than 200 bands but with a spatial resolution of several ten meters.

One way to surmount the abovementioned technological limitation of acquiring high-resolution hyperspectral images is to capture multiple multiband images of the same scene with practical spatial and spectral resolutions, then fuse them together in a synergistic manner. Fusing multiband images combines their complementary information obtained through multiple sensors that may have different spatial and spectral resolutions and cover different spectral ranges.

Initial multiband image fusion algorithms were developed to fuse a panchromatic image with a multispectral image and the associated inverse problem was dubbed pansharpening [18,19,20,21,22]. Many pansharpening algorithms are based on either of the two popular pansharpening strategies: component substitution (CS) and multiresolution analysis (MRA). The CS-based algorithms substitute a component of the multispectral image obtained through a suitable transformation by the panchromatic image. The MRA-based algorithms inject the spatial detail of the panchromatic image obtained by a multiscale decomposition, e.g., using wavelets [23], into the multispectral image. There also exist hybrid methods that use both CS and MRA. Some of the algorithms originally proposed for pansharpening have been successfully extended to be used for fusing a panchromatic image with a hyperspectral image, a problem that is called hyperspectral pansharpening [21].

Recently, significant research effort has been expended to solve the problem of fusing a multispectral image with a hyperspectral one. This inverse problem is essentially different from the pansharpening and hyperspectral pansharpening problems since a multispectral image has multiple bands that are intricately related to the bands of its corresponding hyperspectral image. Unlike a panchromatic image that contains only one band of reflectance data usually covering parts of the visible and near-infrared spectral ranges, a multispectral image contains multiple bands each covering a smaller spectral range, some being in the shortwave-infrared (SWIR) region. Therefore, extending the pansharpening techniques so that they can be used to inject the spatial details of a multispectral image into a hyperspectral image is not straightforward. Nonetheless, an effort towards this end has led to the development of a framework called hypersharpening, which is based on adapting the MRA-based pansharpening methods to multispectral–hyperspectral image fusion. The main idea is to synthesize a high-spatial-resolution image for each band of the hyperspectral image by linearly combining the bands of the multispectral image using linear regression [24].

In some works on multispectral–hyperspectral image fusion, it is assumed that each pixel on the hyperspectral image, which has a lower spatial resolution than the target image, is the average of the pixels of the same area on the target image [25,26,27,28,29]. Clearly, the size of this area depends on the downsampling ratio. Based on this pixel-aggregation assumption, one can divide the problem of fusing two multiband images into subproblems dealing with smaller blocks and hence significantly decrease the complexity of the overall process. However, it is more realistic to allow the area on the target image corresponding to a pixel of the hyperspectral image to span as many pixels as determined by the point-spread function of the sensor, which induces spatial blurring. The downsampling ratio generally depends on the physical and optical characteristics of a sensor and is usually fixed. Therefore, spatial blurring and downsampling can be expressed as two separate linear operations. The spectral degradation of a panchromatic or multispectral image with respect to the target image can also be modeled as a linear transformation. Articulating the spatial and spectral degradations in terms of linear operations forms a realistic and convenient forward observation model to relate the observed multiband images to the target image.

Hyperspectral image data is generally known to have a low-rank structure and reside in a subspace that usually has a dimension much smaller than the number of the spectral bands [8,30,31,32,33]. This is mainly due to correlations among the spectral bands and the fact that the spectrum of each pixel can often be represented as a linear combination of a relatively few spectral signatures. These signatures, called endmembers, may be the spectra of the material present at the scene. Consequently, a hyperspectral image can be linearly decomposed into its constituent endmembers and the fractional abundances of the endmembers for each pixel. This linear decomposition is called spectral unmixing and the corresponding data model is called the linear mixture model. Other linear decompositions that can be used to reduce the dimensionality of a hyperspectral image in the spectral domain are dictionary-learning-based sparse representation and principle-component analysis.

Many recent works on multiband image fusion, which mostly deal with fusing a multispectral image with a hyperspectral image of the same scene, employ the abovementioned forward observation model and a form of linear spectral decomposition. They mostly extract the endmembers or the spectral dictionary from the hyperspectral image. Some of the works use the extracted endmember or dictionary matrix to reconstruct the multispectral image via sparse regression and calculate the endmember abundances or the representation coefficients [34]. Others cast the multiband image fusion problem as reconstructing a high-spatial-resolution hyperspectral datacube from two datacubes degraded according to the mentioned forward observation model. When the number of spectral bands in the multispectral image is smaller than the number of endmembers or dictionary atoms, the linear inverse problem associated with the multispectral–hyperspectral fusion problem is ill-posed and needs be regularized to have a meaningful solution. Any prior knowledge about the target image can be used for regularization. Natural images are known to mostly consist of smooth segments with few abrupt changes corresponding to the edges and object boundaries [35,36,37]. Therefore, penalizing the total-variation [38,39,40] and sparse (low-rank) representation in the spatial domain [41,42,43,44] are two popular approaches to regularizing the multiband image fusion problems. Some algorithms, developed within the framework of the Bayesian estimation, incorporate the prior knowledge or conjecture about the probability distribution of the target image into the fusion problem [45,46,47]. The work of [48] obviates the need for regularization by dividing the observed multiband images into small spatial patches for spectral unmixing and fusion under the assumption that the target image is locally low-rank.

When the endmembers or dictionary atoms are induced from an observed hyperspectral image, the problem of fusing the hyperspectral image with a multispectral image boils down to estimating the endmember abundances or representation coefficients of the target image, a problem that is often tractable (due to being a convex optimization problem) and has a manageable size and complexity. The estimate of the target image is then obtained by mixing the induced endmembers/dictionary and the estimated abundances/coefficients. It is also possible to jointly estimate the endmembers/dictionary and the abundances/coefficients from the available multiband data. This joint estimation problem is usually formulated as a non-convex optimization problem of non-negative matrix factorization, which can be solved approximately using block coordinate-descent iterations [49,50,51,52].

To the best of our knowledge, all existing multiband image fusion algorithms are designed to fuse a pair of multiband images with complementary spatial and spectral resolutions. Therefore, fusing more than two multiband images using the existing algorithms can only be realized by performing a hierarchical procedure that combines multiple fusion processes possibly implemented via different algorithms as, for example, in [53,54]. In addition, there are potentially various ways to arrange the pairings and often it is not possible to know beforehand which way will provide the best overall fusion result. For instance, in order to fuse a panchromatic, a multispectral, and a hyperspectral image of a scene, one can first fuse the panchromatic and multispectral images, then fuse the resultant pansharpened multispectral image with the hyperspectral image. Another way would be to first fuse the multispectral and hyperspectral images, then pansharpen the resultant hyperspectral image with the panchromatic image. Apart from the said ambiguity of choice, such combined pair-wise fusions can be slow and inaccurate since they may require several runs of different algorithms and may suffer from propagation and accumulation of errors. Therefore, the increasing availability of multiband images with complementary characteristics captured by modern spectral imaging devices has brought about the demand for efficient and accurate fusion techniques that can handle multiple multiband images simultaneously.

In this paper, we propose an algorithm that can simultaneously fuse an arbitrary number of multiband images. We utilize the forward observation and linear mixture models to effectively model the data and reduce the dimensionality of the problem. Assuming matrix normal distribution for the observation noise, we derive the likelihood function as well as the Fisher information matrix (FIM) associated with the problem of recovering the endmember abundance matrix of the target image from the observations. We study the properties of the FIM and the conditions for existence of a unique maximum-likelihood estimate and the associated Cramer–Rao lower bound. We regularize the problem of maximum-likelihood estimation of the endmember abundances by adding a vector total-variation penalty term to the cost function and constraining the abundances to be non-negative and add up to one for each pixel. The total-variation penalty serves two major purposes. First, it helps us cope with the likely ill-posedness of the maximum-likelihood estimation problem. Second, it allows us to take into account the spatial characteristics of natural images that is they mostly consist of piecewise plane regions with few sharp variations. Regularization with a vector total-variation penalty can effectively advocate this desired feature by promoting sparsity in the image gradient, i.e., local differences between adjacent pixels, while encourages the local differences to be spatially aligned across different bands [37]. The non-negativity and sum-to-one constraints on the endmember abundances ensure that the abundances have practical values. They also implicitly promote sparsity in the estimated endmember abundances.

We solve the resultant constrained optimization problem using the alternating direction method of multipliers (ADMM) [55,56,57,58,59,60]. Simulation results indicate that the proposed algorithm outperforms several combinations of the state-of-the-art algorithms, which need be cascaded to carry out fusion of multiple (more than two) multiband images.

2. Data Model

2.1. Forward Observation Model

Let us denote the target multiband image by $X\in {ℝ}^{L\times N}$ where $L$ is the number of spectral bands and $N$ is the number of pixels in the image. We wish to recover $X$ from $K$ observed multiband images $Y_k\in {ℝ}^{L_k\times N_k}$, $k=1,\dots ,K$, that are spatially or spectrally downgraded and degraded versions of $X$. We assume that these multiband images are geometrically coregistered and are related to $X$ via the following forward observation model

$$Y_k=R_{k}X{B}_k{S}_k+P_k$$

(1)

where

$L_k\le L$ and $N_k=N/D_k^2$ with $D_k$ being the spatial downsampling ratio of the $k$th image;

$R_k\in {ℝ}^{L_k\times N}$ is the spectral response of the sensor producing $Y_k$;

$B_k\in {ℝ}^{N\times N}$ is a band-independent spatial blurring matrix that represents a two-dimensional convolution with a blur kernel corresponding to the point-spread function of the sensor producing $Y_k$;

$S_k\in {ℝ}^{N\times N_k}$ is a sparse matrix with $N_k$ ones and zeros elsewhere that implements a two-dimensional uniform downsampling of ratio $D_k$ on both spatial dimensions and satisfies $S_k^{\top }{S}_k=I_N$;

$P_k\in {ℝ}^{L_k\times N_k}$ is an additive perturbation representing the noise or error associated with the observation of $Y_k$.

We assume that the perturbations $P_k$, $k=1,\dots ,K$, are independent of each other and have matrix normal distributions expressed by

$$P_k\sim \mathcal{M}{\mathcal{N}}_{L_k\times N_k}\left({\mathbf{0}}_{L_k\times N_k},{\Sigma }_k,I_{N_k}\right)$$

(2)

where ${\mathbf{0}}_{L_k\times N_k}$ is the $L_k\times N_k$ zero matrix, $I_{N_k}$ is the $N_k\times N_k$ identity matrix, and ${\Sigma }_k\in {ℝ}^{L_k\times L_k}$ is a diagonal matrix that represents the correlation among rows of $P_k$, which correspond to different spectral bands. Note that we consider the column-covariance matrices to be identity assuming that the perturbations are independent and identically-distributed in the spatial domain. However, by considering diagonal row-covariance matrices, we assume that the perturbations are independent in the spectral domain but may have nonidentical variances at different bands. Moreover, the instrument noise of an optoelectronic device can also have a multiplicative nature. A prominent example is the shot noise that is generally modeled using the Poisson distribution. By virtue of the central limit theorem and since a Poisson distribution with a reasonably large mean can be well approximated by a Gaussian distribution, our assumption of additive Gaussian perturbation for the acquisition noise/error is a sensible working hypothesis given that the SNRs are adequately high.

Note that $Y_k$, $k=1,\dots ,K$, in (1) contain the corrected (preprocessed) spectral values, not the raw measurements produced by the spectral imagers. The preprocessing usually involves several steps including radiometric calibration, geometric correction, and atmospheric compensation [61]. The radiometric calibration is generally performed to obtain radiance values at the sensor. It converts the sensor measurements in digital numbers into physical units of radiance. The reflected sunlight passing through the atmosphere is partially absorbed and scattered through a complex interaction between the light and various parts of the atmosphere. The atmospheric compensation counters these effects and converts the radiance values into ground-leaving radiance or surface reflectance values. To obtain accurate reflectance values, one additionally has to account for the effects of the viewing geometry and sun’s position as well as the surfaces structural and optical properties [10]. This preprocessing is particularly important when the multiband images to be fused are acquired via different instruments, from different viewpoints, or at different times. After the preprocessing, the images should also be coregistered.

2.2. Linear Mixture Model

Under some mild assumptions, multiband images of natural scenes can be suitably described by a linear mixture model [8]. Specifically, the spectrum of each pixel can often be written as a linear mixture of a few archetypal spectral signatures known as endmembers. The number of endmembers, denoted by $M$, is usually much smaller than the spectral dimension of a hyperspectral image, i.e., $M\ll L$. Therefore, if we arrange $M$ endmembers corresponding to $X$ as columns of the matrix $E\in {ℝ}^{L\times M}$, we can factorize $X$ as

$$X=EA+P$$

(3)

where $A\in {ℝ}^{M\times N}$ is the matrix of endmember abundances and $P\in {ℝ}^{L\times N}$ is a perturbation matrix that accounts for any possible inaccuracy or mismatch in the linear mixture mode. We assume that $P$ is independent of $P_k$, $k=1,\dots ,K$, and has a matrix normal distribution as

$$P\sim \mathcal{M}{\mathcal{N}}_{L\times N}\left({\mathbf{0}}_{L\times N},\Sigma ,I_N\right)$$

(4)

where $\Sigma \in {ℝ}^{L\times L}$ is its row-covariance matrix. Every column of $A$ contains the fractional abundances of the endmembers at a pixel. The fractional abundances are non-negative and often assumed to add up to one for each pixel.

The linear mixture model stated above has been widely used in various contexts and applications concerning multiband, particularly hyperspectral, images. Its popularity can mostly be attributed to its intuitiveness as well as relative simplicity and ease of implementation. However, remotely-sensed images of ground surface may suffer from strong nonlinear effects. These effects are generally due to ground characteristics such as non-planar surface and bidirectional reflectance, artefacts left by atmospheric removal procedures, and the presence of considerable atmospheric absorbance in the neighborhood of the bands of interest. The use of images of the same scene captured by different sensors, although coregistered, can also induce nonlinearly mainly owing to difference in observation geometry, lighting conditions, and miscalibration. Therefore, it should be taken into consideration that the mentioned nonlinear phenomena can impact the results of any procedure relying on the linear mixture model in any real-world application and the scale of the impact depends on the severity of the nonlinearities.

There are a few other caveats regarding the linear mixture model that should also be kept in mind. First, $X$ in model (3) corresponds to a matrix of corrected (preprocessed) values, not raw ones that would typically be captured by a spectral imager of the same spatial and spectral resolutions. However, whether these values are radiance or reflectance has no impact on the validity of the model, though it certainly matters for further processing of the data. Second, the model (3) does not necessarily require each endmember to be the spectral signature of only one (pure) material. An endmember may be composed of the spectral signatures of multiple materials or may be seen as the spectral signature of a composite material made of several constituent materials. Additionally, depending on the application, the endmembers may be purposely defined in particular subjective ways. Third, in practice, an endmember may have slightly different spectral manifestations at different parts of a scene due to variable illumination, environmental, atmospheric, or temporal conditions. This so-called endmember variability [62] along with possible nonlinearities in the actual underlying mixing process [63] may introduce inaccuracies or inconsistencies in the linear mixture model and consequently in the endmember extraction or spectral unmixing techniques that rely on this model. Lastly, the sum-to-one assumption on the abundances of each pixel may not always hold, especially, when the linear mixture model is not able to account for every material in a pixel possibly because of the effects of endmember variability or nonlinear mixing.

2.3. Fusion Model

Substituting (3) into (1) gives

$$Y_k=R_{k}EA{B}_k{S}_k+{\check{P}}_k$$

(5)

where the aggregate perturbation of the $k$th image is

$${\check{P}}_k=P_k+R_{k}P{B}_k{S}_k.$$

(6)

Instead of estimating the target multiband image $X$ directly, we consider estimating its abundance matrix $A$ from the observations $Y_k$, $k=1,..K$, given the endmember matrix $E$. We can then obtain an estimate of the target image by multiplying the estimated abundance matrix by the endmember matrix. This way, we reduce the dimensionality of the fusion problem and consequently the associated computational burden. In addition, by estimating $A$ first, we attain an unmixed fused image obviating the need to perform additional unmixing, if demanded by any application utilizing the fused image. However, this approach requires the prior knowledge of the endmember matrix $E$. The columns of this matrix can be selected from a library of known spectral signatures, such as the U.S. Geological Survey digital spectral library [64], or extracted from the observed multiband images that have the appropriate spectral dimension.

3. Problem

3.1. Maximum-Likelihood Estimation

In order to facilitate our analysis, we define the following vectorized variables

$$y_k=\mathrm{vec}\left\{Y_k\right\}\in {ℝ}^{L_k{N}_k\times 1}$$

(7a)

$$a=\mathrm{vec}\left\{A\right\}\in {ℝ}^{MN\times 1}$$

(7b)

$$p_k=\mathrm{vec}\left\{{\check{P}}_k\right\}\in {ℝ}^{L_k{N}_k\times 1}$$

(7c)

where $\mathrm{vec}\{·\}$ is the vectorization operator that stacks the columns of its matrix argument on top of each other. Applying $\mathrm{vec}\{·\}$ to both sides of (5) while using the property $\mathrm{vec}\left\{ABC\right\}=\left(C^{\top }⨂A\right)\mathrm{vec}\left\{B\right\}$ gives

$$y_k=\left(S_k^{\top }{B}_k^{\top }⨂R_{k}E\right)a+p_k$$

(8)

where $⨂$ denotes the Kronecker product.

Since $P_k$ and $P$ have independent matrix normal distributions (see (2) and (4)), $p_k$ has a multivariate normal distribution expressed as

$$p_k~{\mathcal{N}}_{L_k{N}_k}\left({\mathbf{0}}_{L_k{N}_k},I_{N_k}⨂{\Sigma }_k+S_k^{\top }{B}_k^{\top }{B}_k{S}_k⨂R_k\Sigma R_k^{\top }\right)$$

(9)

where ${\mathbf{0}}_{L_k{N}_k}$ stands for the $L_k{N}_k\times 1$ vector of zeroes. Using the approximation $S_k^{\top }{B}_k^{\top }{B}_k{S}_k\approx {\xi }_k{I}_{N_k}$ with ${\xi }_k>0$, we get

$$p_k ~ {\mathcal{N}}_{L_k{N}_k}\left({\mathbf{0}}_{L_k{N}_k},I_{N_k}⨂{\Lambda }_k\right)$$

(10)

where

$${\Lambda }_k={\Sigma }_k+{\xi }_k{R}_k\Sigma R_k^{\top }$$

(11)

In view of (9) and (10), we have

$$y_k ~ {\mathcal{N}}_{L_k{N}_k}\left(\left[S_k^{\top }{B}_k^{\top }⨂R_{k}E\right]a,I_{N_k}⨂{\Lambda }_k\right).$$

(12)

Hence, the probability density function of $y_k$ parametrized over the unknown $a$ can be written as

$$\begin{array}{ll}{f}_{y_k}\left(y_k;a\right)& ={\left|2\pi I_{N_k}⨂{\Lambda }_k\right|}^{-\frac{1}{2}}\\ & \times \exp \left\{-\frac{1}{2}{\left[y_k-\left(S_k^{\top }{B}_k^{\top }⨂R_{k}E\right)a\right]}^{\top }{\left(I_{N_k}⨂{\Lambda }_k\right)}^{-1}\left[y_k-\left(S_k^{\top }{B}_k^{\top }⨂R_{k}E\right)a\right]\right\}.\end{array}$$

(13)

Since the perturbations $p_k$, $k=1,\dots ,K$, are independent of each other, the joint probability density function of the observations is written as

$$\begin{array}{ll}{f}_{y_1,\dots ,y_K}\left(y_1,\dots ,y_K;a\right)& ={{\displaystyle \prod }}_{k=1}^K{f}_{y_k}\left(y_k;a\right)\\ & ={{\displaystyle \prod }}_{k=1}^K{\left|2\pi I_{N_k}⨂{\Lambda }_k\right|}^{-\frac{1}{2}}\exp \left\{-\frac{1}{2}{\displaystyle {\displaystyle \sum }_{k=1}^K}||\left(I_{N_k}⨂{\Lambda }_k^{-1/2}\right){\left[y_k-\left(S_k^{\top }{B}_k^{\top }⨂R_{k}E\right)a\right]||}^2\right\}\end{array}$$

(14)

and the log-likelihood function of $a$ given the observed data as

$$\begin{array}{ll}l\left(a|y_1,\dots ,y_K\right)& =\ln f_{y_1,\dots ,y_K}\left(y_1,\dots ,y_K;a\right)\\ & =-\frac{1}{2}\ln \left({{\displaystyle \prod }}_{k=1}^K\left|2\pi I_{N_k}⨂{\Lambda }_k\right|\right)-\frac{1}{2}{\displaystyle {\displaystyle \sum }_{k=1}^K}||\left(I_{N_k}⨂{\Lambda }_k^{-1/2}\right){\left[y_k-\left(S_k^{\top }{B}_k^{\top }⨂R_{k}E\right)a\right]||}^2.\end{array}$$

(15)

Accordingly, the maximum-likelihood estimate of $a$ is found by solving the following optimization problem

$$\begin{array}{ll}\widehat{a}& =\arg \underset{a}{\max } l\left(a|y_1,\dots ,y_K\right)\\ & =\underset{a}{\mathrm{argmin}} \frac{1}{2}{\displaystyle {\displaystyle \sum }_{k=1}^K}||\left(I_{N_k}⨂{\Lambda }_k^{-1/2}\right){\left[y_k-\left(S_k^{\top }{B}_k^{\top }⨂R_{k}E\right)a\right]||}^2.\end{array}$$

(16)

This problem can be stated in terms of $A={\mathrm{vec}}^{-1}\left\{a\right\}$ as

$$\widehat{A}=\underset{A}{\mathrm{argmin}} \frac{1}{2}{\displaystyle \sum }_{k=1}^K||{\Lambda }_k^{-1/2}{\left(Y_k-R_{k}EA{B}_k{S}_k\right)||}_{\mathrm{F}}^2.$$

(17)

The Fisher information matrix (FIM) of the maximum-likelihood estimator $\widehat{a}$ in (16) is calculated as

$$\mathcal{F}=-\mathrm{E}\left[{\mathcal{H}}_l\left(a\right)\right]$$

(18)

where ${\mathcal{H}}_l\left(a\right)$ denotes the Hessian, i.e., the Jacobian of the gradient, of the log-likelihood function $l\left(a|y_1,\dots ,y_K\right)$. The entry on the $i$th row and the $j$th column of ${\mathcal{H}}_l\left(a\right)$ is computed as

$$\frac{{\partial }^2}{\partial a_i\partial a_j}l\left(a|y_1,\dots ,y_K\right)$$

(19)

where $a_i$ and $a_j$ denote the $i$th and $j$th entries of $a$, respectively. Accordingly, we can show that

$$\mathcal{F}={\displaystyle \sum }_{k=1}^K\left(B_k{S}_k{S}_k^{\top }{B}_k^{\top }⨂E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E\right).$$

(20)

If $\mathcal{F}$ is invertible, the optimization problem (16) has a unique solution given by

$$\widehat{a}={\left[{\displaystyle \sum }_{k=1}^K\left(B_k{S}_k{S}_k^{\top }{B}_k^{\top }⨂E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E\right)\right]}^{-1}{\displaystyle \sum }_{k=1}^K\left(B_k{S}_k⨂E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}\right)y_k$$

(21)

and the Cramer–Rao lower bound for the estimator $\widehat{a}$, which is a lower bound on the covariance of $\widehat{a}$, is the inverse of $\mathcal{F}$. The FIM $\mathcal{F}$ is guaranteed to be invertible when, for at least one image, the matrix $B_k{S}_k{S}_k^{\top }{B}_k^{\top }⨂E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E$ is full-rank.

The matrix $S_k{S}_k^{\top }$ has a rank of $N_k$ hence for $D_k>1$ is rank-deficient. The blurring matrix $B_k$ does not change the rank of the matrix that it multiplies from the right. In addition, as ${\Lambda }_k^{-1}$ is full-rank, $E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E$ has a full rank of $M$ when the rows of $R_{k}E$ are at least as many as its columns, i.e., $L_k\ge M$. Therefore, $A$ and consequently $X$ is guaranteed to be uniquely identifiable given $Y_k$, $k=1,\dots ,K$, only when at least one observed image, say the $q$th image, has full spatial resolution, i.e., $N_q=N$, with the number of its spectral bands being equal to or larger than the number of endmembers, i.e., $L_q\ge M$, so that, at least for the $q$th image, $B_q{S}_q{S}_q^{\top }{B}_q^{\top }⨂E^{\top }{R}_q^{\top }{\Lambda }_q^{-1}{R}_{q}E$ is full-rank.

In practice, it is rarely possible to satisfy the abovementioned requirements as multiband images with high spectral resolution are generally spatially downsampled and the number of bands of the ones with full spatial resolution, such as panchromatic or multispectral images, is often less than the number of endmembers. Hence, the inverse problem of recovering $A$ from $Y_k$, $k=1,\dots ,K$, is usually ill-posed or ill-conditioned. Thus, some prior knowledge need be injected into the estimation process to produce a unique and reliable estimate. The prior knowledge is intended to partially compensate for the information lost in spectral and spatial downsampling and usually stems from experimental evidence or common facts that may induce certain analytical properties or constraints. The prior information is commonly incorporated into the problem in the form of imposed constraints or additive regularization terms. Examples of prior knowledge about $A$ that are regularly used in the literature are non-negativity and sum-to-one constraints, matrix normal distribution with known or estimated parameters [45], sparse representation with a learned or known dictionary or basis [41], and minimal total variation [38].

3.2. Regularization

To develop an algorithm for effective fusion of multiple multiband images with arbitrary spatial and spectral resolutions, we employ two mechanisms to regularize the maximum-likelihood cost function in (17).

As the first regularization mechanism, we impose a constraint on $A$ such that its entries are non-negative and sum to one in all columns. We express this constraint as $A\ge 0$ and ${\mathbf{1}}_M^{\top }A={\mathbf{1}}_N^{\top }$ where $A\ge 0$ means all the entries of $A$ are greater than or equal to zero. As the second regularization mechanism, we add an isotropic vector total-variation penalty term, denoted by $||\nabla {A||}_{2,1}$, to the cost function. Here, ${||·||}_{2,1}$ is the l${}_{2,1}$-norm operator that returns the sum of l${}_2$-norms of all the columns of its matrix argument. In addition, we define

$$\nabla A=\left[\begin{array}{c}\mathbf{A}{D}_h\\ \mathbf{A}{D}_v\end{array}\right]\in {ℝ}^{2M\times N}$$

(22)

where $D_h$ and $D_v$ are discrete differential matrix operators that, respectively, yield the horizontal and vertical first-order backward differences (gradients) of the row-vectorized image that they multiply from the right. Consequently, we formulate our regularized optimization problem for estimating $A$ as

$$\begin{gathered} \underset{A}{\min }\frac{1}{2}{\displaystyle \sum }_{k=1}^K||{\Lambda }_k^{-1/2}{\left(Y_k-R_{k}EA{B}_k{S}_k\right)||}_F^2+\alpha ||\nabla {A||}_{2,1} \\ \mathrm{subject}\mathrm{to}:A\ge 0\mathrm{and}{\mathbf{1}}_M^{\top }A={\mathbf{1}}_N^{\top } \end{gathered}$$

(23)

where $\alpha \ge 0$ is the regularization parameter.

The non-negativity and sum-to-one constraints on $A$, which force the columns of $A$ to reside on the unit $\left(M-1\right)$-simplex, are naturally expected and help find a solution that is physically plausible. In addition, they implicitly induce sparseness in the solution. The total-variation penalty promotes solutions with a sparse gradient, a property that is known to be possessed by images of most natural scenes as they are usually made of piecewise homogeneous regions with few sudden changes at object boundaries or edges. Note that the subspace spanned by the endmembers is the one that the target image $X$ lives in. Therefore, through the total-variation regularization of the abundance matrix $A$, we regularize $X$ indirectly.

4. Algorithm

Defining the set of values for $A$ that satisfy the non-negativity and sum-to-one constraints as

$$\mathcal{S}=\left\{A|A\ge 0,{\mathbf{1}}_M^{\top }A={\mathbf{1}}_N^{\top }\right\}$$

(24)

and making use of the indicator function ${ı}_{\mathcal{S}}\left(A\right)$ defined as

$${ı}_{\mathcal{S}}\left(A\right)=\{\begin{array}{cc}0& A\in \mathcal{S}\\ +\infty & A\notin \mathcal{S},\end{array}$$

(25)

we rewrite (23) as

$$\underset{A}{\min } \frac{1}{2}{\displaystyle \sum }_{k=1}^K||{\Lambda }_k^{-1/2}{\left(Y_k-R_{k}EA{B}_k{S}_k\right)||}_{\mathrm{F}}^2+\mathsf{\alpha }||\nabla {A||}_{2,1}+{ı}_{\mathcal{S}}\left(A\right).$$

(26)

4.1. Iterations

We use the alternating direction method of multipliers (ADMM), also known as the split-Bregman method, to solve the convex but nonsmooth optimization problem of (26). We split the problem to smaller and more manageable pieces by defining the auxiliary variables, $U_k\in {ℝ}^{M\times N}$, $k=1,\dots ,K$, $V\in {ℝ}^{2M\times N}$, and $W\in {ℝ}^{M\times N}$, and changing (26) into

$$\begin{gathered} \underset{A,U_1,\dots ,U_K,V,W}{\min } \frac{1}{2}{\displaystyle \sum }_{k=1}^K||{\Lambda }_k^{-1/2}{\left(Y_k-R_{k}E{U}_k{S}_k\right)||}_{\mathrm{F}}^2+\mathsf{\alpha }||{V||}_{2,1}+{ı}_{\mathcal{S}}\left(W\right) \\ \mathrm{subject}\mathrm{to}:U_k=A{B}_k,V=\nabla A,W=A. \end{gathered}$$

(27)

Then, we write the augmented Lagrangian function associated with (27) as

$$\begin{gathered} \mathcal{L}\left(A,U_1,\dots ,U_K,V,W,F_1,\dots ,F_K,G,H\right)=\frac{1}{2}{\displaystyle \sum }_{k=1}^K||{\Lambda }_k^{-1/2}{\left(Y_k-R_{k}E{U}_k{S}_k\right)||}_{\mathrm{F}}^2+\mathsf{\alpha }||{V||}_{2,1}+{ı}_{\mathcal{S}}\left(W\right) \\ +\frac{\mu }{2}{\displaystyle \sum }_{k=1}^K||A{B}_k-U_k-F_k{||}_{\mathrm{F}}^2+\frac{\mu }{2}||\nabla A-V-{G||}_{\mathrm{F}}^2+\frac{\mu }{2}||A-W-{H||}_{\mathrm{F}}^2 \end{gathered}$$

(28)

where $F_k\in {ℝ}^{M\times N}$, $k=1,\dots ,K$, $G\in {ℝ}^{2M\times N}$, and $H\in {ℝ}^{M\times N}$ are the scaled Lagrange multipliers and $\mu \ge 0$ is the penalty parameter.

Using the ADMM, we minimize the augmented Lagrangian function in an iterative fashion. At each iteration, we alternate the minimization with respect to the main unknown variable $A$ and the auxiliary variables; then, we update the scaled Lagrange multipliers. Hence, we compute the iterates as

$$A^{\left(n\right)}=\underset{A}{\mathrm{argmin}} \mathcal{L}\left(A,U_1^{\left(n-1\right)},\dots ,U_K^{\left(n-1\right)},V^{\left(n-1\right)},W^{\left(n-1\right)},F_1^{\left(n-1\right)},\dots ,F_K^{\left(n-1\right)},G^{\left(n-1\right)},H^{\left(n-1\right)}\right)$$

(29)

$$\begin{gathered} \left\{U_1^{\left(n\right)},\dots ,U_K^{\left(n\right)},V^{\left(n\right)},W^{\left(n\right)}\right\} \\ =\underset{U_1,\dots ,U_K,V,W}{\mathrm{argmin}}\mathcal{L}\left(A^{\left(n\right)},U_1,\dots ,U_K,V,W,F_1^{\left(n-1\right)},\dots ,F_K^{\left(n-1\right)},G^{\left(n-1\right)},H^{\left(n-1\right)}\right) \end{gathered}$$

(30)

$$F_k^{\left(n\right)}=F_k^{\left(n-1\right)}-\left(A^{\left(n\right)}{B}_k-U_k^{\left(n\right)}\right), k=1,\dots ,K$$

(31a)

$$G^{\left(n\right)}=G^{\left(n-1\right)}-\left(\nabla A^{\left(n\right)}-V^{\left(n\right)}\right)$$

(31b)

$$H^{\left(n\right)}=H^{\left(n-1\right)}-\left(A^{\left(n\right)}-W^{\left(n\right)}\right)$$

(31c)

where superscript $\left(n\right)$ denotes the value of an iterate at iteration number $n\ge 0$. We repeat the iterations until convergence is reached up to a maximum allowed number of iterations.

Since we define the auxiliary variables independent of each other, the minimization of the augmented Lagrangian function (28) with respect to the auxiliary variables can be realized separately. Thus, (30) is equivalent to

$$U_k^{\left(n\right)}=\underset{U_k}{\mathrm{argmin}} \frac{1}{2}||{\Lambda }_k^{-1/2}{\left(Y_k-R_{k}E{U}_k{S}_k\right)||}_{\mathrm{F}}^2+\frac{\mu }{2}||A^{\left(n\right)}{B}_k-U_k-F_k^{\left(n-1\right)}||{}_{\mathrm{F}}^2, k=1,\dots ,K$$

(32)

$$V^{\left(n\right)}=\underset{V}{\mathrm{argmin}} \mathsf{\alpha }||{V||}_{2,1}+\frac{\mu }{2}||\nabla A^{\left(n\right)}-V-G^{\left(n-1\right)}||{}_{\mathrm{F}}^2$$

(33)

$$W^{\left(n\right)}=\underset{W}{\mathrm{argmin}} {ı}_{\mathcal{S}}\left(W\right)+\frac{\mu }{2}||A^{\left(n\right)}-W-H^{\left(n-1\right)}||{}_{\mathrm{F}}^2.$$

(34)

4.2. Solutions of Subproblems

Considering (28), (29) can be written as

$$A^{\left(n\right)}=\underset{A}{\mathrm{argmin}} {\displaystyle \sum }_{k=1}^K||A{B}_k-U_k^{\left(n-1\right)}-F_k^{\left(n-1\right)}||{}_{\mathrm{F}}^2+||\nabla A-V^{\left(n-1\right)}-G^{\left(n-1\right)}||{}_{\mathrm{F}}^2+||A-W^{\left(n-1\right)}-H^{\left(n-1\right)}||{}_{\mathrm{F}}^2 .$$

(35)

Calculating the gradient of the cost function in (35) with respect to $A$ and setting it to zero gives

$$\begin{gathered} A^{\left(n\right)}=\left[{\displaystyle \sum }_{k=1}^K\left(U_k^{\left(n-1\right)}+F_k^{\left(n-1\right)}\right)B_k^{\top } +Q_1^{\left(n-1\right)}{D}_h^{\top }+Q_2^{\left(n-1\right)}{D}_v^{\top }+W^{\left(n-1\right)}+H^{\left(n-1\right)}\right] \\ \times {\left({\displaystyle \sum }_{k=1}^K{B}_k{B}_k^{\top }+D_h{D}_h^{\top }+D_v{D}_v^{\top }+I_N\right)}^{-1} \end{gathered}$$

(36)

where, for the convenience of presentation, we define $Q_1^{\left(n-1\right)}$ and $Q_2^{\left(n-1\right)}$ as

$$\left[\begin{array}{c}{Q}_1^{\left(n-1\right)}\\ Q_2^{\left(n-1\right)}\end{array}\right]=V^{\left(n-1\right)}+G^{\left(n-1\right)}.$$

(37)

To make the computation of $A^{\left(n\right)}$ in (36) more efficient, we assume that the two-dimensional convolutions represented by $B_k$, $k=1,\dots ,K$, are cyclic. In addition, we assume that the differential matrix operators $D_h$ and $D_v$ apply with periodic boundaries. Consequently, multiplications by $B_k^{\top }$, $D_h^{\top }$, and $D_v^{\top }$ as well as by ${\left({{\displaystyle \sum }}_{k=1}^K{B}_k{B}_k^{\top }+D_h{D}_h^{\top }+D_v{D}_v^{\top }+I_N\right)}^{-1}$ can be performed through the use of the fast Fourier transform (FFT) algorithm and the circular convolution theorem. This theorem states that the Fourier transform of a circular convolution is the pointwise product of the Fourier transforms, i.e., a circular convolution can be expressed as the inverse Fourier transform of the product of the individual spectra [65].

Equating the gradient of the cost function in (32) with respect to $U_k$ to zero results in

$$E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E{U}_k^{\left(n\right)}{S}_k{S}_k^{\top }+\mu U_k^{\left(n\right)}=E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{Y}_k{S}_k^{\top }+\mu \left(A^{\left(n\right)}{B}_k-F_k^{\left(n-1\right)}\right).$$

(38)

Multiplying both sides of (38) from the right by the masking matrix $M_k=S_k{S}_k^{\top }$ and its complement $I_N-M_k$ yields

$$U_k^{\left(n\right)}{M}_k={\left(E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E+\mu I_N\right)}^{-1}\left[E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{Y}_k{S}_k^{\top }+\mu \left(A^{\left(n\right)}{B}_k-F_k^{\left(n-1\right)}\right)M_k\right]$$

(39)

and

$$U_k^{\left(n\right)}\left(I_N-M_k\right)=\left(A^{\left(n\right)}{B}_k-F_k^{\left(n-1\right)}\right)\left(I_N-M_k\right),$$

(40)

respectively. Note that we have $S_k^{\top }{S}_k=I_N$ and $M_k$ is idempotent, i.e., $M_k{M}_k=M_k$. Summing both sides of (39) and (40) gives the solution of (32) for $k=1,\dots ,K$ as

$$\begin{array}{ll}{U}_k^{\left(n\right)}& =U_k^{\left(n\right)}{M}_k+U_k^{\left(n\right)}\left(I_N-M_k\right)\\ & ={\left(E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E+\mu I_N\right)}^{-1}\left[E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{Y}_k{S}_k^{\top }+\mu \left(A^{\left(n\right)}{B}_k-F_k^{\left(n-1\right)}\right)M_k\right]+\left(A^{\left(n\right)}{B}_k-F_k^{\left(n-1\right)}\right)\left(I_N-M_k\right).\end{array}$$

(41)

The terms ${\left(E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E+\mu I_N\right)}^{-1}$ and $E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{Y}_k{S}_k^{\top }$ do not change during the iterations and can be precomputed.

The subproblem (33) can be decomposed pixelwise and its solution is linked to the so-called Moreau proximity operator of the ${\mathcal{l}}_{2,1}$-norm given by column-wise vector-soft-thresholding [66,67]. If we define

$$Z^{\left(n\right)}=\nabla A^{\left(n\right)}-G^{\left(n-1\right)},$$

(42)

the $j$th column of $V^{\left(n\right)}$, denoted by $v_j^{\left(n\right)}$, is given in terms of the $j$th column of $Z^{\left(n\right)}$, denoted by $z_j^{\left(n\right)}$, as

$$v_j^{\left(n\right)}=\frac{\max \{||z_j^{\left(n\right)}||{}_2-\frac{\mathsf{\alpha }}{\mu },0\}}{||z_j^{\left(n\right)}||{}_2}{z}_j^{\left(n\right)}.$$

(43)

The solution of (34) is the value of the proximity operator of the indicator function ${ı}_{\mathcal{S}}\left(W\right)$ at the point $A^{\left(n\right)}-H^{\left(n-1\right)}$, which is the projection of $A^{\left(n\right)}-H^{\left(n-1\right)}$ onto the set $\mathcal{S}$ defined by (24). Therefore, we have

$$\begin{array}{ll}{W}^{\left(n\right)}& =\underset{W\in \mathcal{S}}{\mathrm{argmin}}|| A^{\left(n\right)}-H^{\left(n-1\right)}-{W||}_{\mathrm{F}}^2\\ & ={\mathsf{\Pi }}_{\mathcal{S}}\left\{A^{\left(n\right)}-H^{\left(n-1\right)}\right\}\end{array}$$

(44)

where ${\mathsf{\Pi }}_{\mathcal{S}}\{·\}$ denotes the projection onto $\mathcal{S}$. We implement this projection onto the unit $\left(M-1\right)$-simplex employing the algorithm proposed in [68].

Algorithm 1 presents a summary of the proposed algorithm.

Algorithm 1 The proposed algorithm

1: initialize 2:    $E\leftarrow \mathrm{VCA}\left(Y_l\right)$ % if $E$ is not known and $Y_l$ has full spectral resolution 3:    $A^{\left(0\right)}\leftarrow$ upscale and interpolate the output of $\mathrm{SUnSAL}\left(Y_l,E\right)$ 4:    for $k=1,\dots ,K$ 5:      $U_k^{\left(0\right)}=A^{\left(0\right)}$ 6:      $F_k^{\left(n\right)}={\mathbf{0}}_{M\times N}$ 7:    $V^{\left(0\right)}=A^{\left(0\right)}$, $W^{\left(0\right)}=A^{\left(0\right)}$ 8:    $G^{\left(0\right)}={\mathbf{0}}_{M\times N}$, $H^{\left(0\right)}={\mathbf{0}}_{M\times N}$ 9: for $n=1,2,\dots$ % until a convergence criterion is met or a given maximum number of iterations is reached 10:    $\left[\begin{array}{c}{Q}_1^{\left(n-1\right)}\\ Q_2^{\left(n-1\right)}\end{array}\right]=V^{\left(n-1\right)}+G^{\left(n-1\right)}$ 11:    $A^{\left(n\right)}=\left[{{\displaystyle \sum }}_{k=1}^K\left(U_k^{\left(n-1\right)}+F_k^{\left(n-1\right)}\right)B_k^{\top } +Q_1^{\left(n-1\right)}{D}_h^{\top }+Q_2^{\left(n-1\right)}{D}_v^{\top }+W^{\left(n-1\right)}+H^{\left(n-1\right)}\right]$        $\times {\left({{\displaystyle \sum }}_{k=1}^K{B}_k{B}_k^{\top }+D_h{D}_h^{\top }+D_v{D}_v^{\top }+I_N\right)}^{-1}$ 12:    for $k=1,\dots ,K$ 13:       $U_k^{\left(n\right)}={\left(E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{R}_{k}E+\mu I_N\right)}^{-1}\left[E^{\top }{R}_k^{\top }{\Lambda }_k^{-1}{Y}_k{S}_k^{\top }+\mu \left(A^{\left(n\right)}{B}_k-F_k^{\left(n-1\right)}\right)M_k\right]$          $+\left(A^{\left(n\right)}{B}_k-F_k^{\left(n-1\right)}\right)\left(I_N-M_k\right)$ 14:    $Z^{\left(n\right)}=\nabla A^{\left(n\right)}-G^{\left(n-1\right)}$ 15:    for $j=1,\dots ,N$ 16:       $v_j^{\left(n\right)}=\frac{\max \{||z_j^{\left(n\right)}{||}_2-\frac{\mathsf{\alpha }}{\mu },0\}}{||z_j^{\left(n\right)}{||}_2}{z}_j^{\left(n\right)}$ 17:    $W^{\left(n\right)}={\mathsf{\Pi }}_{\mathcal{S}}\left\{A^{\left(n\right)}-H^{\left(n-1\right)}\right\}$ 18:    for $k=1,\dots ,K$ 19:       $F_k^{\left(n\right)}=F_k^{\left(n-1\right)}-\left(A^{\left(n\right)}{B}_k-U_k^{\left(n\right)}\right)$ 20:    $G^{\left(n\right)}=G^{\left(n-1\right)}-\left(\nabla A^{\left(n\right)}-V^{\left(n\right)}\right)$ 21:    $H^{\left(n\right)}=H^{\left(n-1\right)}-\left(A^{\left(n\right)}-W^{\left(n\right)}\right)$ 22: calculate the fused image 23:    $\widehat{X}=E{A}^{\left(n\right)}$

4.3. Convergence

By defining

$$\mathcal{U}={\left[U_1,\cdots ,U_K,V,W\right]}^{\top }$$

(45)

and

$$\mathcal{C}={\left[B_1,\cdots ,B_K,D_h,D_v,I_N\right]}^{\top },$$

(46)

(27) can be expressed as

$$\begin{gathered} \underset{A}{\min } f\left(\mathcal{U}\right) \\ \mathrm{subject}\mathrm{to}\mathcal{U}=\mathcal{C}{A}^{\top } \end{gathered}$$

(47)

where

$$f\left(\mathcal{U}\right)=\frac{1}{2}{\displaystyle \sum }_{k=1}^K||{\Lambda }_k^{-1/2}{\left(Y_k-R_{k}E{U}_k{S}_k\right)||}_{\mathrm{F}}^2+\mathsf{\alpha }||{V||}_{2,1}+{ı}_{\mathcal{S}}\left(W\right).$$

(48)

The function $f\left(\mathcal{U}\right)$ is closed, proper, and convex as it is a sum of closed, proper, and convex functions and $\mathcal{C}$ has full column rank. Therefore, according to Theorem 8 of [56], if (47) has a solution, the proposed algorithm converges to this solution, regardless of the initial values as long as the penalty parameter $\mu$ is positive. If no solution exists, at least one of $A^{\left(n\right)}$ and ${\mathcal{U}}^{\left(n\right)}$ will diverge.

5. Simulations

To examine the performance of the proposed algorithm in comparison with the state-of-the-art, we simulate the fusion of three multiband images, viz. a panchromatic image, a multispectral image, and a hyperspectral image. To this end, we adopt the popular practice known as the Wald’s protocol [69], which is to use a reference image with high spatial and spectral resolutions to generate the lower-resolution images that are fused and evaluate the fusion performance by comparing the fused image with the reference image.

We obtain the reference images of our experiments by cropping five publicly available hyperspectral images to the spatial resolutions given in

Table 1. The spatial and spectral dimensions of the considered hyperspectral datasets (reference images) and the value of the regularization parameter used in the proposed algorithm with each dataset.

Image	No. of Rows	No. of Columns	No. of Bands	$\mathit{\alpha }$
Botswana	$400$	$240$	$145$	$5$
Indian Pines	$400$	$400$	$200$	$7$
Washington DC Mall	$400$	$300$	$191$	$5$
Moffett Field	$480$	$320$	$176$	$22$
Kennedy Space Center	$500$	$400$	$176$	$28$

. These images are called Botswana [70], Indian Pines [71], Washington DC Mall [71], Moffett Field [17], and Kennedy Space Center [70]. The Botswana image has been captured by the Hyperion sensor aboard the Earth Observing 1 (EO-1) satellite, the Washington DC Mall image by the airborne-mounted Hyperspectral Digital Imagery Collection Experiment (HYDICE), and the Indian Pines, Moffett Filed, and Kennedy Space Center images by the NASA Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) instrument. All images cover the visible near-infrared (VNIR) and short-wavelength-infrared (SWIR) ranges with uncalibrated, excessively noisy, and water-absorbed bands removed. The spectral resolution of each image is also given in Table 1. The data as well as the MATLAB code used to produce the results of this paper can be found at [72].

We generate three multiband images (panchromatic, multispectral, and hyperspectral) using each reference image. We obtain the hyperspectral images by applying a rotationally-symmetric 2D Gaussian blur filter with a kernel size of $13\times 13$ and a standard deviation of $2.12$ to each reference image followed by downsampling with a ratio of $4$ in both spatial dimensions for all bands. For the multispectral images, we use a Gaussian blur filter with a kernel size of $7\times 7$ and a standard deviation of $1.06$ and downsampling with a ratio of $2$ in both spatial dimensions for all bands of each reference image. Afterwards, we downgrade the resultant images spectrally by applying the spectral responses of the Landsat 8 multispectral sensor. This sensor has eight multispectral bands and one panchromatic band. Figure 1 depicts the spectral responses of all the bands of this sensor [73]. We create the panchromatic images from the reference images using the panchromatic band of the Landsat 8 sensor without applying any spatial blurring or downsampling. We add zero-mean Gaussian white noise to each band of the produced multiband images such that the band-specific signal-to-noise ratio (SNR) is $30$ dB for the multispectral and hyperspectral images and $40$ dB for the panchromatic image. In practice, SNR may vary along the bands of a multiband sensor and the noise may be non-zero-mean or non-Gaussian. Our use of the same SNR for all bands and zero-mean Gaussian noise is a simplification adopted for the purpose of evaluating the proposed algorithm and comparing it with the considered benchmarks.

Note that we have selected the standard deviations of the abovementioned 2D Gaussian blur filters such that the normalized magnitude of the modulation transfer function (MTF) of both filters is approximately $0.25$ at the Nyquist frequency in both spatial dimensions [74] as shown in Figure 2. We have also selected the filter kernel sizes in accordance with the downsampling ratios and the selected standard deviations. In our simulations, we use symmetric 2D Gaussian blur filters for simplicity and ease of computations. Gaussian blur filters are known to well approximate the real acquisition MTFs, which may be affected by a number of physical processes. With a pushbroom acquisition, the MTF is different in the across-track and along-track directions. This is because, in the along-track direction, the blurring is due to the effects of the spatial resolution and the apparent motion of the scene. On the other hand, in the across-track direction, the blurring is mainly attributable to the resolution of the instrument, fixed by both detector and optics. There is also the contribution to MTF by the propagation trough a scattering atmosphere.

The current multiband image fusion algorithms published in the literature are designed to fuse two images at a time. In order to compare the performance of the proposed algorithm with the state-of-the-art, we consider fusing the abovementioned three multiband images in three different ways, which we refer to as Pan + HS, Pan + (MS + HS), and (Pan + MS) + HS, using the existing algorithms for pansharpening, hyperspectral pansharpening, and hyperspectral-multispectral fusion. In Pan + HS, we only fuse the panchromatic and hyperspectral images. In Pan + (MS + HS), and (Pan + MS) + HS, we fuse the given images in two cascading stages. In Pan + (MS + HS), first, we fuse the multispectral and hyperspectral images. Then, we fuse the resultant hyperspectral image with the panchromatic image. We use the same algorithm at both stages, albeit with different parameter values. In (Pan + MS) + HS, we first fuse the panchromatic image with the multispectral one. Then, we fuse the pansharpened multispectral image with the hyperspectral image. We use two different algorithms at each of the two stages resulting in four combined solutions.

For pansharpening, which is the fusion of a panchromatic image with a multispectral one, we use two algorithms called the band-dependent spatial detail (BDSD) [75] and the modulation-transfer-function generalized Laplacian pyramid with high-pass modulation (MTF-GLP-HPM) [76,77,78]. The BDSD algorithm belongs to the class of component substitution methods and the MTF-GLP-HPM algorithm falls into the category of multiresolution analysis. In [18], where several pansharpening algorithms are studied, it is shown that the BDSD and MTF-GLP-HPM algorithms exhibit the best performance among all the considered ones.

For fusing a panchromatic or multispectral image with a hyperspectral image, we use two algorithms proposed in [38,79,80], which are called HySure and R-FUSE-TV, respectively. These algorithms are based on total-variation regularization and are among the best performing and most efficient hyperspectral pansharpening and multispectral–hyperspectral fusion algorithms currently available [21,81].

We use three performance metrics for assessing the quality of a fused image with respect to its reference image. The metrics are the relative dimensionless global error in synthesis (ERGAS) [82], spectral angle mapper (SAM) [83], and $Q{2}^n$ [84]. The metric $Q{2}^n$ is a generalization of the universal image quality index (UIQI) proposed in [85] and an extension of the $Q4$ index [86] to hyperspectral images based on hypercomplex numbers.

We extract the endmembers (columns of $E$) from each hyperspectral image using the vertex component analysis (VCA) algorithm [87]. The VCA is a fast unsupervised unmixing algorithm that assumes the endmembers as the vertices of a simplex encompassing the hyperspectral data cloud. We utilize the SUnSAL algorithm [88] together with the extracted endmembers to unmix each hyperspectral image and obtain its abundance matrix. Then, we upscale the resulting matrix by a factor of four and apply two-dimensional spline interpolation on each of its rows (abundance bands) to generate the initial estimate for the abundance matrix $A^{\left(0\right)}$. We initialize the proposed algorithm as well as the HySure and R-FUSE-TV algorithms by this matrix.

To make our comparisons fair, we tune the values of the parameters in the HySure and R-FUSE-TV algorithms to yield the best possible performance in all experiments. In addition, in order to use the BDSD and MTF-GLP-HPM algorithms to their best potential, we provide these algorithms with the true point-spread function, i.e., the blurring kernel, used to generate the multispectral images.

Apart from the number of endmembers, which can be estimated using, for example, the HySime algorithm [31], the proposed algorithm has two tunable parameters, the total-variation regularization parameter $\mathsf{\alpha }$ and the ADMM penalty parameter $\mu$. The automatic tuning of the values of these parameters is an interesting and challenging subject. There are a number of strategies that can be employed such as those proposed in [66,89]. We found through experimentations that although the value of $\mu$ impacts the convergence speed of the proposed algorithm, as long as it is within an appropriate range, it has little influence on the accuracy of the proposed algorithm. Therefore, we set it to $\mu =1.5\times {10}^3$ in all experiments. The value of $\mathsf{\alpha }$ affects the performance of the proposed algorithm in subtle ways as shown in Figure 3 where we plot the performance metrics, ERGAS, SAM, and $Q{2}^n$, against $\mathsf{\alpha }$ for the Botswana and Washington DC Mall images. The results in Figure 3 suggest that, for different values of $\mathsf{\alpha }$, there is a trade-off between the performance metrics, specifically, ERGAS and $Q{2}^n$ on one side and SAM on the other. Therefore, we tune the value of $\mathsf{\alpha }$ for each experiment only roughly to obtain a reasonable set of values for all three performance metrics. We give the values of $\mathsf{\alpha }$ used in the proposed algorithm in Table 1.

In Table 2, we give the values of the performance metrics to assess the quality of the images fused using the proposed algorithm and the considered benchmarks. We provide the performance metrics for the case of considering only the bands within the spectrum of the panchromatic image as well as the case of considering all bands, i.e., the entire spectrum of the reference image. We also give the time taken by each algorithm to produce the fused images. We used MATLAB (The MathWorks, Natick, MA, USA) with a 2.9-GHz Core-i7 CPU and 24 GB of DDR3 RAM and ran each of the proposed, HySure, and R-FUSE-TV algorithms for $200$ iterations as they always converged sufficiently after this number of iterations. According to the results in Table 2, the proposed algorithm significantly outperforms the considered benchmarks. It is also evident from the required processing times that the computational (time) complexity of the proposed algorithm is lower than those of its contenders.

In Figure 4 and Figure 5, we plot the sorted per-pixel normalized root mean-square error (NRMSE) values of the proposed algorithm and the best performing algorithms from each of the Pan + HS, Pan + (MS + HS), and (Pan + MS) + HS categories. Figure 4 corresponds to the case of considering only the spectrum of the panchromatic image and Figure 5 to the case of considering the entire spectrum. We define the per-pixel NRMSE as $||x_j-{\widehat{x}}_j{||}_2 /||x_j{||}_2$ where $x_j$ and ${\widehat{x}}_j$ are the $j$th column of the reference image $X$ and the fused image $\widehat{X}$, respectively. We sort the NRMSE values in the ascending order.

In Figure 6, we show RGB renderings of the reference images together with the panchromatic, multispectral, and hyperspectral images generated from them and used for the fusion. We also show the fused images yielded by the proposed algorithm and Pan + (MS + HS) fusion using the HySure algorithm, which generally performs better than the other considered benchmarks. The multispectral images are depicted using their red, green, and blue bands. The RGB representations of the hyperspectral images are rendered through transforming the spectral data to the CIE XYZ color space and then transforming the XYZ values to the sRGB color space. From visual inspection of the reference and fused images shown in Figure 6, it is observed that the images fused by the proposed algorithm match their corresponding reference images better than the ones produced by the Pan + (MS + HS) fusion using the HySure algorithm do.

6. Conclusions

We proposed a new image fusion algorithm that can simultaneously fuse multiple multiband images. We utilized the well-known forward observation model together with the linear mixture model to cast the fusion problem as a reduced-dimension linear inverse problem. We used a vector total-variation penalty as well as non-negativity and sum-to-one constraints on the endmember abundances to regularize the associated maximum-likelihood estimation problem. The regularization encourages the estimated fused image to have low rank with a sparse representation in the spectral domain while preserving the edges and discontinuities in the spatial domain. We solved the regularized problem using the alternating direction method of multipliers. We demonstrated the advantages of the proposed algorithm in comparison with the state-of-the-art via experiments with five real hyperspectral images that were done following the Wald’s protocol.