Sparsity-Based Spatiotemporal Fusion via Adaptive Multi-Band Constraints

Abstract

Remote sensing is an important means to monitor the dynamics of the earth surface. It is still challenging for single-sensor systems to provide spatially high resolution images with high revisit frequency because of the technological limitations. Spatiotemporal fusion is an effective approach to obtain remote sensing images high in both spatial and temporal resolutions. Though dictionary learning fusion methods appear to be promising for spatiotemporal fusion, they do not consider the structure similarity between spectral bands in the fusion task. To capitalize on the significance of this feature, a novel fusion model, named the adaptive multi-band constraints fusion model (AMCFM), is formulated to produce better fusion images in this paper. This model considers structure similarity between spectral bands and uses the edge information to improve the fusion results by adopting adaptive multi-band constraints. Moreover, to address the shortcomings of the ${\ell }_1$ norm which only considers the sparsity structure of dictionaries, our model uses the nuclear norm which balances sparsity and correlation by producing an appropriate coefficient in the reconstruction step. We perform experiments on real-life images to substantiate our conceptual augments. In the empirical study, the near-infrared (NIR), red and green bands of Landsat Enhanced Thematic Mapper Plus (ETM+) and Moderate Resolution Imaging Spectroradiometer (MODIS) are fused and the prediction accuracy is assessed by both metrics and visual effects. The experiments show that our proposed method performs better than state-of-the-art methods. It also sheds light on future research.

1. Introduction

Remote sensing satellites are important tools for the monitoring of processes such as vegetation and land cover changes on the earth surface [1,2,3]. Because of technological limitations in sensor designs [4], compromises have to be made between spatial and temporal resolutions. For example, Moderate Resolution Imaging Spectroradiometer (MODIS) can visit the earth once a day with 500-m spatial resolution. As a comparison, the spatial resolution of Landsat Enhanced Thematic Mapper Plus (ETM+) is 30 m, but its revisiting period is 16 days. Such a limitation restricts the application of remote sensing in problems that need images high in both spatial and temporal resolutions. Spatiotemporal reflectance fusion models [5] have thus been developed to fuse image data from different sensors to obtain high spatiotemporal resolution images.

The Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) [6] is a pioneering fusion model based on a weighting method. This model uses neighboring pixels to compute the center pixel at a point in time with a weighting function, and the weights are determined by spectral difference, temporal difference and location distance. Furthermore, Zhu et al. [7] proposed an Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model (ESTARFM) based on a STARFM algorithm to predict the surface reflectance of heterogeneous landscapes. Another improved STARFM method is a Spatial Temporal Adaptive Algorithm for mapping Reflectance Change (STAARCH) [8], which is designed to detect disturbance and changes in reflectance by using Tasseled Cap transformations. However, performance of the weighting methods are constrained because linear combination smooths out the changing terrestrial contents.

Another type of reflectance fusion method, known as dictionary learning methods, has been proposed to overcome the shortcoming of the weighting methods. Dictionary-based methods that use certain known dictionaries, such as wavelets and shearlets, have been proved to be efficient in multisensor and multiresolution image fusion [9,10,11]. In remote sensing data analysis, Moigne et al. [12] and Czaja et al. [13] proposed remote sensing image fusion methods based on wavelets and wavelet packet, respectively. Shearlet transform is also used in a fusion algorithm in [14] because shearlets can share the same optimality properties and enjoy similar geometrical properties. Using the capability of dictionary learning and sparsity-based methods in the super resolution analysis, Huang et al. [15] proposed a Sparse-representation-based Spatiotemporal Reflectance Fusion Model (SPSTFM) to integrate sparse representation and reflectance fusion by establishing correspondences between structures in high resolution images and their corresponding low resolution images through dictionary pair and sparse coding. SPSTFM assumes that high and low resolution images of the same area have the same sparse coefficients. Such assumption is, however, too restrictive [16]. Based on this idea, Wu et al. [17] proposed the Error-Bound-regularized Semi-Coupled Dictionary Learning (EBSCDL) model which assumes that the representation coefficients of the image pair have a stable mapping and coefficients of the dictionary pair have perturbations in the reconstruction step. Attempts have been made to improve the performance of the SCDL based models. For examples, Block Sparse Bayesian Learning for Semi-Coupled Dictionary Learning (BSBL-SCDL) [18] employs the structural sparsity of the sparse coefficients as a priori knowledge and Compressed Sensing for Spatiotemporal Fusion (CSSF) [19] considers explicitly the down-sampling process within the framework of compressed sensing for reconstruction. In comparison with the weighting methods, the advantage of the dictionary-learning-based methods is that they retrieve the hidden relationship between image pairs from the sparse coding space to better capture structure changes.

Besides the aforementioned methods, some researchers employed other approaches to fuse multi-source data. Unmixing techniques have been suggested for spatiotemporal fusion because of their ability to reconstruct images with high spectral fidelity [20,21,22,23,24]. Considering the mixed-class spectra within a coarse pixel, Xu et al. [25] proposed the Class Regularized Spatial Unmixing (CRSU) model. This method is based on the conventional spatial unmixing technique but is modified to include prior class spectra estimated by the known image pairs. To provide a formal statistical framework for fusion, Xue et al. [26] proposed Spatiotemporal Bayesian Data Fusion (STBDF) that makes use of the joint distribution to capture implicitly the temporal changes of images for the estimation of the high resolution image at a target point in time.

Considering structure similarity in spectral bands, structure information has been employed in pan-sharpening and image fusion. Shi et al. [27], for example, proposed a learning interpolation method for pan-sharpening by expanding sketch information of the high-resolution panchromatic (PAN) image which contains the structure features of the PAN image. Glasner et al. [28] verified that many structures in a natural image are similar at the same and different scales. Inspired by this, a self-learning approach was proposed by Khateri et al. [29] which uses similar structures at different levels to pan-sharpen the low resolution multi-spectral images. In multi-modality image fusion, Zhu et al. [30] proposed a method which decomposes images into cartoon and texture components, and preserves the structure information of two components based on spatial-based method and sparse representation, respectively.

However, none of these spatiotemporal fusion methods consider the structure similarity between spectral bands in the fusion procedure. Although different bands have different reflectance ranges, the edge information is still similar [31]. Obviously, a reconstruction model can have a better performance if such information can be effectively used to predict the unknown high resolution image. Otherwise, the dictionary pair obtained by the training image pair are inefficient to predict the unknown images because of the lack of information for the target time. This can be explained from the experience in machine learning in which the ${\ell }_1$ norm is too restrictive in encoding the unknown data in the prediction process because it only uses the sparsity structure of the dictionary [32,33]. Therefore, the reconstruction model needs a replacement of the ${\ell }_1$ norm to reduce the impact of insufficient information and to improve the representation ability of the dictionary pair.

We propose a new model in this paper to enhance spatiotemporal fusion performance. Our model uses the edge information in different bands via adaptive multi-band constraints to improve the reconstruction performance. To overcome the disadvantage of the ${\ell }_1$ norm, the nuclear norm is adopted as the regularization term to increase the efficiency of the learnt dictionary pair. Nuclear norm considers not only the sparsity but also the coordination in producing a suitable coefficient that can harmonize the sparse and collaborative representations adaptively [32,33].

Overall, the main contributions of this work can be summarized as follows.

• The multi-band constraints are employed to reinforce the structure similarity of different bands in spatiotemporal fusion.
• Considering the different structure similarity between two bands, the adaptive regularization parameters are proposed to determine the importance of each multi-band constraint adaptively.
• The nuclear norm is employed to replace the ${\ell }_1$ norm in the reconstruction model because the nuclear norm considers both sparsity and correlation of the dictionaries and can overcome the disadvantage of the ${\ell }_1$ norm.

The remainder of this paper is organized as follows. Our method for spatiotemporal fusion, called adaptive multi-band constraints fusion model (AMCFM), is proposed in Section 2. Section 3 discusses the experiments carried out to assess the effectiveness of the AMCFM and four state-of-the-art methods in terms of statistics and visual effects. We then conclude the paper with a summary and direction for future research in Section 4.

2. Methodology

2.1. Problem Definition

In the following, MODIS images are selected as low resolution images and Landsat ETM+ images are selected as high resolution images. As shown in Figure 1, our spatiotemporal fusion model requires three low resolution images ${\mathit{M}}_1$, ${\mathit{M}}_2$ and ${\mathit{M}}_3$, and two high resolution images ${\mathit{L}}_1$ and ${\mathit{L}}_3$. The high resolution image ${\mathit{L}}_2$ is the target image that we want to predict. Let ${\mathit{L}}_{ij}$ (${\mathit{L}}_{ij}={\mathit{L}}_i-{\mathit{L}}_j$) and ${\mathit{M}}_{ij}$ (${\mathit{M}}_{ij}={\mathit{M}}_i-{\mathit{M}}_j$) be the high and low resolution difference images between $t_i$ and $t_j$ ($i,j\in \{1,2,3\}$), respectively. We assume that changes of remote sensing images between two points in time are linear. For effectiveness, the dictionary pair ${\mathit{D}}_l$ and ${\mathit{D}}_m$ is trained by the difference image pair ${\mathit{L}}_{31}$ and ${\mathit{M}}_{31}$ [15]. Then, the high resolution difference image ${\mathit{L}}_{21}$ can be produced by using the dictionary pair to encode the corresponding low resolution difference image ${\mathit{M}}_{21}$. ${\mathit{L}}_{32}$ can be obtained in the same way. Finally, the high resolution image at time $t_2$ can be predicted as follows:

$${\mathit{L}}_2=W_{21}\ast ({\mathit{L}}_1+{\mathit{L}}_{21})+W_{32}\ast ({\mathit{L}}_3-{\mathit{L}}_{32}).$$

(1)

The weights $W_{21}$ and $W_{32}$ we used are same as those in [19], which take the average of the two predicted difference images.

2.2. Dictionary Learning Fusion Model

As mentioned above, the conventional dictionary learning fusion models are usually comprised of two steps: the dictionary pair training step and the reconstruction step. The whole process is performed on each band separately. Here, we show the mathematical formulation of these two steps in SPSTFM [15], which is the initial dictionary learning model.

• Dictionary pair training step:
  In the training step, the difference image pair, ${\mathit{L}}_{31}$ and ${\mathit{M}}_{31}$, is used to train the high resolution dictionary ${\mathit{D}}_l$ and the corresponding low resolution dictionary ${\mathit{D}}_m$ as follows:

  $$[{\mathit{D}}_l,{\mathit{D}}_m]=arg\underset{{\mathit{D}}_l,{\mathit{D}}_m}{min}{∥\begin{array}{c}\mathit{Y}-{\mathit{D}}_l\mathit{A}\end{array}∥}_2^2+{∥\begin{array}{c}\mathit{X}-{\mathit{D}}_m\mathit{A}\end{array}∥}_2^2+\lambda {∥\begin{array}{c}\mathit{A}\end{array}∥}_1,$$

  (2)

  where $\mathit{Y}$ and $\mathit{X}$ are the column combination of the lexicographically stacked image patches, sampled randomly from ${\mathit{L}}_{31}$ and ${\mathit{M}}_{31}$, respectively. $\mathit{A}$ is the column combination of the representation coefficients corresponding to every column in $\mathit{Y}$ and $\mathit{X}$, and $\lambda$ is the regularization parameter. We adopt the K-SVD (K is the abbreviation for K-means and SVD is the abbreviation for Singular Value Decomposition) lgorithm [34] to solve for ${\mathit{D}}_l$ and ${\mathit{D}}_m$ in Equation (2).
• Reconstruction step:
  Then, ${\mathit{D}}_m$ is used to encode each patch of ${\mathit{M}}_{21}$ and the sparse coding coefficient $\mathit{\alpha }$ is obtained by solving the optimization problem:

  $${\mathit{\alpha }}^{\ast }=arg\underset{\mathit{\alpha }}{min}\frac{1}{2}{∥\begin{array}{c}{\mathit{m}}_{21}-{\mathit{D}}_m\mathit{\alpha }\end{array}∥}_2^2+\lambda {∥\begin{array}{c}\mathit{\alpha }\end{array}∥}_1,$$

  (3)

  where ${\mathit{m}}_{21}$ is a patch of ${\mathit{M}}_{21}$. The corresponding patch of the high resolution image can be produced by

  $${\mathit{l}}_{21}={\mathit{D}}_l{\mathit{\alpha }}^{\ast }.$$

  (4)
  Then, all patches of ${\mathit{L}}_{21}$ are merged to get the high resolution image ${\mathit{L}}_{21}$. ${\mathit{L}}_{32}$ can be obtained in the same way and the target image ${\mathit{L}}_2$ can be predicted through Equation (1).

2.3. Adaptive Multi-Band Constraints Fusion Model

Our model uses the same strategy for dictionary pair training and focuses on the improvement of the reconstruction step. We propose the following model for spatiotemporal fusion by replacing the ${\ell }_1$ norm with the nuclear norm ${\parallel ·\parallel }_{\ast }$ and incorporating the multi-band constraints. The reconstruction formulation is given as follows:

$$\begin{array}{cc}\hfill [{\mathit{\alpha }}_N^{\ast },{\mathit{\alpha }}_R^{\ast },{\mathit{\alpha }}_G^{\ast }]=argmin\sum _{c\in \{N,R,G\}}\frac{1}{2}& {∥\begin{array}{c}{\mathit{m}}_c-{\mathit{D}}_{m_c}{\mathit{\alpha }}_c\end{array}∥}_2^2+\lambda {∥\begin{array}{c}{\mathit{D}}_{m_c}\mathrm{Diag}({\mathit{\alpha }}_c)\end{array}∥}_{\ast }\hfill \\ \hfill +{\tau }_{NR}& {∥\begin{array}{c}S{\mathit{D}}_{l_N}{\mathit{\alpha }}_N-S{\mathit{D}}_{l_R}{\mathit{\alpha }}_R\end{array}∥}_2^2\hfill \\ \hfill +{\tau }_{RG}& {∥\begin{array}{c}S{\mathit{D}}_{l_R}{\mathit{\alpha }}_R-S{\mathit{D}}_{l_G}{\mathit{\alpha }}_G\end{array}∥}_2^2\hfill \\ \hfill +{\tau }_{GN}& {∥\begin{array}{c}S{\mathit{D}}_{l_G}{\mathit{\alpha }}_G-S{\mathit{D}}_{l_N}{\mathit{\alpha }}_N\end{array}∥}_2^2,\hfill \end{array}$$

(5)

where $\lambda$, ${\tau }_{NR}$, ${\tau }_{RG}$ and ${\tau }_{GN}$ are the regularization parameters. ${\parallel M\parallel }_{\ast }$ is the nuclear norm of M that is the sum of all the singular values of matrix M. For a vector $\mathit{v}$, $\mathrm{Diag}(\mathit{v})$ represents a diagonal matrix whose diagonal elements are the corresponding elements of the vector $\mathit{v}$. S is a high-pass detector filter. Here, we choose the two-dimensional Laplacian operator. The subscripts N, R and G mean the near-infrared (NIR), red and green band, respectively.

The dictionary pair ${\mathit{D}}_l$ and ${\mathit{D}}_m$ is trained by the difference images ${\mathit{L}}_{31}$ and ${\mathit{M}}_{31}$ which do not contain sufficient information of the images at time $t_2$. When reconstructing ${\mathit{L}}_{21}$ or ${\mathit{L}}_{32}$, if the model only uses the ${\ell }_1$ norm regularization, then the performance are unsatisfactory. It is more reasonable to integrate sparsity and correlation of the dictionaries. The nuclear norm term is just the kind of regularization that can adaptively balance sparsity and correlation via a suitable coefficient. As the property shown in [33,35],

$$\begin{array}{c}\hfill {∥\begin{array}{c}{\mathit{\alpha }}_c\end{array}∥}_2\le {∥\begin{array}{c}{\mathit{D}}_{m_c}\mathrm{Diag}({\mathit{\alpha }}_c)\end{array}∥}_{\ast }\le {∥\begin{array}{c}{\mathit{\alpha }}_c\end{array}∥}_1,\end{array}$$

(6)

where all columns of ${\mathit{D}}_{m_c}$ have unit norm. When the column vectors of ${\mathit{D}}_{m_c}$ are orthogonal, ${∥\begin{array}{c}{\mathit{D}}_{m_c}\mathrm{Diag}({\mathit{\alpha }}_c)\end{array}∥}_{\ast }$ is equal to ${∥\begin{array}{c}{\mathit{\alpha }}_c\end{array}∥}_1$. When the column vectors of ${\mathit{D}}_{m_c}$ are highly correlated, ${∥\begin{array}{c}{\mathit{D}}_{m_c}\mathrm{Diag}({\mathit{\alpha }}_c)\end{array}∥}_{\ast }$ will be close to ${∥\begin{array}{c}{\mathit{\alpha }}_c\end{array}∥}_2$ [35]. Generally, remote sensing images in the dictionary ${\mathit{D}}_{m_c}$ are neither too independent nor too correlated because the test images and training images can contain high correlative information (i.e., stable land-cover) and independent information (i.e., land-cover change). Therefore, as shown in Equation (6), our model can benefit from both the ${\ell }_1$ norm and the ${\ell }_2$ norm. The advantage of the nuclear norm is that the nuclear norm can capture the correlation structure of the training images which the ${\ell }_1$ norm cannot.

The last three terms in the model are the multi-band regularization terms. Taking the NIR band as an example, ${\mathit{D}}_{l_N}{\mathit{\alpha }}_N$ denotes a high resolution patch of the NIR band and $S{\mathit{D}}_{l_N}{\mathit{\alpha }}_N$ stands for the edge in the patch. These terms make the sparse codes (The codes may not be sparse based on the nuclear norm regularization, but, for convenience and without confusion, we still call them sparse codes.) in different bands no longer independent and reinforce the structure similarity of different bands.

Nevertheless, the nuclear norm regularization and multi-band regularization make it more complicated to solve the model. In Section 2.5, we propose the method to get a solution efficiently.

2.4. Adaptive Parameters between Bands

The ranges of reflectance in different bands are different in remote sensing images. In natural images, the range of the three channels is $[0,255]$.

Table 1. Mean and standard deviation in three bands of a multi-band image.

Index	Image 1			Image 2			Image 3
	NIR	Red	Green	NIR	Red	Green	NIR	Red	Green
Mean Value	0.2212	0.0378	0.0483	0.2113	0.0345	0.0476	0.2320	0.0355	0.0487
Standard Deviation	0.0546	0.0105	0.0097	0.0383	0.0141	0.0121	0.0413	0.0108	0.0084

implicitly shows the range differences of three bands in terms of mean and standard deviation. Obviously, the structures are more similar when the means and standard deviations are closer. Based on this rationale, we propose an adaptive regularization parameter as follows:

$$\begin{array}{c}\hfill {\tau }_{NR}=\gamma ·{10}^{-\frac{min({\overline{m}}_N+{\sigma }_N,{\overline{m}}_R+{\sigma }_R)}{max({\overline{m}}_N+{\sigma }_N,{\overline{m}}_R+{\sigma }_R)}},\end{array}$$

(7)

where ${\overline{m}}_c$ is the mean value of band c and ${\sigma }_c$ is the standard deviation of band c; $\gamma$ is a parameter to control the magnitude; and ${\tau }_{RG}$ and ${\tau }_{GN}$ can be obtained by the definition as well.

This parameter estimates the distribution of the reflectance values between two bands and produces a suitable parameter adaptively. When two bands have similar reflectance values, the parameter is close to $\gamma$. When the difference between two bands increases, the parameter decreases exponentially. The more similar two bands are, the more important a role the corresponding term plays in the model. This property fits such intuitive perception.

2.5. Optimization of the AMCFM

In this section, we solve the reconstruction model. For the optimization, a simplification is made first. We introduce the following vectors and matrices:

$$\begin{array}{cc}\hfill \mathit{\alpha }& =\left[\begin{array}{c}{\mathit{\alpha }}_N\\ {\mathit{\alpha }}_R\\ {\mathit{\alpha }}_G\end{array}\right],\mathit{m}=\left[\begin{array}{c}{\mathit{m}}_N\\ {\mathit{m}}_R\\ {\mathit{m}}_G\end{array}\right],{\mathit{D}}_m=\left[\begin{array}{ccc}{\mathit{D}}_{m_N}& 0& 0\\ 0& {\mathit{D}}_{m_R}& 0\\ 0& 0& {\mathit{D}}_{m_G}\end{array}\right],\hfill \\ \hfill \mathit{A}& =\left[\begin{array}{ccc}({\tau }_{GN}+{\tau }_{NR}){\mathit{D}}_{l_N}^T{S}^{T}S{\mathit{D}}_{l_N}& -{\tau }_{NR}{\mathit{D}}_{l_N}^T{S}^{T}S{\mathit{D}}_{l_R}& 0\\ 0& ({\tau }_{NR}+{\tau }_{RG}){\mathit{D}}_{l_R}^T{S}^{T}S{\mathit{D}}_{l_R}& -{\tau }_{RG}{\mathit{D}}_{l_R}^T{S}^{T}S{\mathit{D}}_{l_G}\\ -{\tau }_{GN}{\mathit{D}}_{l_G}^T{S}^{T}S{\mathit{D}}_{l_N}& 0& ({\tau }_{RG}+{\tau }_{GN}){\mathit{D}}_{l_G}^T{S}^{T}S{\mathit{D}}_{l_G}\end{array}\right],\hfill \end{array}$$

where $\mathit{\alpha }$ and $\mathit{m}$ are concatenations of the sparse coefficients and low resolution image patches, respectively. ${\mathit{D}}_m$ is a dictionary that contains low resolution dictionaries of three bands in its diagonals. We also define $\mathit{B}$ as:

$$\mathit{B}=\frac{1}{2}{\mathit{D}}_m^T{\mathit{D}}_m+\mathit{A}.$$

(8)

Then, Equation (5) can be simplified as follows:

$${\mathit{\alpha }}^{\ast }=arg\underset{\mathit{\alpha }}{min}{\mathit{\alpha }}^T\mathit{B}\mathit{\alpha }-{\mathit{m}}^T{\mathit{D}}_m\mathit{\alpha }+\lambda {∥\begin{array}{c}{\mathit{D}}_m\mathrm{Diag}(\mathit{\alpha })\end{array}∥}_{\ast }.$$

(9)

Here, we use the alternating direction method of multipliers (ADMM) [36,37,38] algorithm to approximate the optimal solution of Equation (9). The optimization problem can be written as follows:

$$\begin{array}{c}\underset{\mathit{\alpha }}{min}{\mathit{\alpha }}^T\mathit{B}\mathit{\alpha }-{\mathit{m}}^T{\mathit{D}}_m\mathit{\alpha }+\lambda {∥\begin{array}{c}\mathit{Z}\end{array}∥}_{\ast }\hfill \\ \mathrm{s}.\mathrm{t}.{\mathit{D}}_m\mathrm{Diag}(\mathit{\alpha })-\mathit{Z}=0,\hfill \end{array}$$

(10)

where $\mathit{Z}$ is the dual variable in the ADMM algorithm. The augmented Lagrangian function $\mathit{L}$ of the optimization problem is given as

$$\begin{array}{c}\hfill \mathit{L}={\mathit{\alpha }}^T\mathit{B}\mathit{\alpha }-m^T{\mathit{D}}_m\mathit{\alpha }+\lambda {∥\begin{array}{c}\mathit{Z}\end{array}∥}_{\ast }+\frac{\rho }{2}{∥\begin{array}{c}{\mathit{D}}_m\mathrm{Diag}(\mathit{\alpha })-\mathit{Z}+\mathit{\mu }\end{array}∥}_F^2,\end{array}$$

(11)

where $\rho$ is a positive scalar and $\mathit{\mu }$ is a scaled variable. The ADMM consists of the following iterations:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\mathit{\alpha }}_{k+1}=arg\underset{\mathit{\alpha }}{min}{\mathit{\alpha }}^T\mathit{B}\mathit{\alpha }-{\mathit{m}}^T{\mathit{D}}_m\mathit{\alpha }+\frac{\rho }{2}{∥\begin{array}{c}{\mathit{D}}_m\mathrm{Diag}(\mathit{\alpha })-{\mathit{Z}}_k+{\mathit{\mu }}_k\end{array}∥}_F^2\hfill \\ {\mathit{Z}}_{k+1}=arg\underset{\mathit{Z}}{min}\lambda {∥\begin{array}{c}\mathit{Z}\end{array}∥}_{\ast }+\frac{\rho }{2}{∥\begin{array}{c}{\mathit{D}}_m\mathrm{Diag}({\mathit{\alpha }}_{k+1})-\mathit{Z}+{\mathit{\mu }}_k\end{array}∥}_F^2\hfill \\ {\mathit{\mu }}_{k+1}={\mathit{\mu }}_k+{\mathit{D}}_m\mathrm{Diag}({\mathit{\alpha }}_{k+1})-{\mathit{Z}}_{k+1}.\hfill \end{array}\right.\end{array}$$

(12)

To minimize the augmented Lagrangian function, we solve each of the subproblems in Equation (12) by fixing the other two variables alternatively. For the step of updating ${\mathit{\alpha }}_{k+1}$, ${\mathit{\alpha }}_{k+1}$ can be deduced as follows:

$$\begin{array}{cc}\hfill {\mathit{\alpha }}_{k+1}& =arg\underset{\mathit{\alpha }}{min}\frac{\rho }{2}({\mathit{\alpha }}^T\mathrm{Diag}(\mathrm{diag}({\mathit{D}}_m^T{\mathit{D}}_m))\mathit{\alpha }-2{(\mathrm{diag}({\mathit{Z}}^T{\mathit{D}}_m+{\mathit{\mu }}^T{\mathit{D}}_m))}^T\mathit{\alpha })\hfill \\ & +{\mathit{\alpha }}^T\mathit{B}\mathit{\alpha }-{\mathit{m}}^T{\mathit{D}}_m\mathit{\alpha }\hfill \\ & =arg\underset{\mathit{\alpha }}{min}{\mathit{\alpha }}^T(\frac{\rho }{2}\mathrm{Diag}(\mathrm{diag}({\mathit{D}}_m^T{\mathit{D}}_m))+\mathit{B})\mathit{\alpha }\hfill \\ & -(\rho {(\mathrm{diag}({\mathit{Z}}^T{\mathit{D}}_m+{\mathit{\mu }}^T{\mathit{D}}_m))}^T+{\mathit{m}}^T{\mathit{D}}_m)\mathit{\alpha }\hfill \\ & ={(\mathit{C}+{\mathit{C}}^T)}^{-1}(\rho {(\mathrm{diag}({\mathit{Z}}^T{\mathit{D}}_m+{\mathit{\mu }}^T{\mathit{D}}_m))}^T+{\mathit{m}}^T{\mathit{D}}_m),\hfill \end{array}$$

(13)

where $\mathit{C}=\frac{\rho }{2}\mathrm{Diag}(\mathrm{diag}({\mathit{D}}_m^T{\mathit{D}}_m))+\mathit{B}$. For a matrix $\mathit{M}$, $\mathrm{diag}(\mathit{M})$ represents a vector whose ith element is the ith diagonal element of the matrix $\mathit{M}$.

For the step of updating ${\mathit{Z}}_{k+1}$, ${\mathit{Z}}_{k+1}$ can be calculated by the singular value thresholding operator [39] as follows:

$$\begin{array}{cc}\hfill {\mathit{Z}}_{k+1}& =arg\underset{\mathit{Z}}{min}\frac{1}{2}{∥\begin{array}{c}\mathit{Z}-({\mathit{D}}_m\mathrm{Diag}({\mathit{\alpha }}_{k+1})+{\mathit{\mu }}_k)\end{array}∥}_F^2+\frac{\lambda }{\rho }{∥\begin{array}{c}\mathit{Z}\end{array}∥}_{\ast }\hfill \\ & ={\mathcal{D}}_{\frac{\lambda }{\rho }}({\mathit{D}}_m\mathrm{Diag}({\mathit{\alpha }}_{k+1})+{\mathit{\mu }}_k),\hfill \end{array}$$

(14)

where ${\mathcal{D}}_{\frac{\lambda }{\rho }}$ is the singular value shrinkage operator, which is defined as follows:

$$\begin{array}{c}\hfill {\mathcal{D}}_{\frac{\lambda }{\rho }}(\mathit{X}):=\mathit{U}{\mathcal{D}}_{\frac{\lambda }{\rho }}(\Sigma ){\mathit{V}}^{\ast },{\mathcal{D}}_{\frac{\lambda }{\rho }}(\Sigma )=\mathrm{diag}(\{max({\sigma }_i-\frac{\lambda }{\rho },0)\}),\end{array}$$

(15)

where $\frac{\lambda }{\rho }$ is a positive scalar, $\mathit{U}\Sigma {\mathit{V}}^{\ast }$ is the singular value decomposition of a matrix $\mathit{X}$, ${\sigma }_i$ is the ith positive singular value of $\mathit{X}$, and $max({\sigma }_i-\frac{\lambda }{\rho },0)=\left\{\begin{array}{cc}{\sigma }_i-\frac{\lambda }{\rho },\hfill & \mathrm{if}{\sigma }_i-\frac{\lambda }{\rho }\ge 0,\hfill \\ 0,\hfill & \mathrm{if}{\sigma }_i-\frac{\lambda }{\rho }<0.\hfill \end{array}\right.$ Now, we use $\mathit{U}\Sigma {\mathit{V}}^{\ast }$ to denote the singular value decomposition of $({\mathit{D}}_m\mathrm{Diag}({\mathit{\alpha }}_{k+1})+{\mathit{\mu }}_k)$ and use ${\sigma }_i$ to denote the ith positive singular value of $({\mathit{D}}_m\mathrm{Diag}({\mathit{\alpha }}_{k+1})+{\mathit{\mu }}_k)$. Then,

$$\begin{array}{c}\hfill {\mathit{Z}}_{k+1}=\mathit{U}{\mathcal{D}}_{\frac{\lambda }{\rho }}(\Sigma ){\mathit{V}}^{\ast }=\mathit{U}\mathrm{diag}(\{max({\sigma }_i-\frac{\lambda }{\rho },0)\}){\mathit{V}}^{\ast }.\end{array}$$

(16)

The implementation details of the whole reconstruction procedure based on the ADMM algorithm can be summarized in Algorithm 1.

Algorithm 1 Reconstruction Procedure of the Proposed Method

1: Input: The regularization parameter $\lambda$, learnt dictionary pair ${\mathit{D}}_l$ and ${\mathit{D}}_m$, low resolution difference image and initial parameters ${\mathit{\alpha }}_0$, ${\mathit{Z}}_0$, ${\mathit{\mu }}_0$, $\rho$. 2: Preprocessing: Normalize the low resolution difference image, then segment the low resolution difference image into patches $\mathit{M}={\{{\mathit{m}}_i\}}_{i=1}^N$ with a 7 × 7 patch size and a four-pixel overlap in each direction. 3: Calculate: The structure similarity parameters ${\tau }_{NR}$, ${\tau }_{RG}$ and ${\tau }_{GN}$. 4: Repeat: (1) Update the sparse coefficient ${\alpha }_{k+1}$ as: ${\mathit{\alpha }}_{k+1}={(\mathit{C}+{\mathit{C}}^T)}^{-1}(\rho {(\mathrm{diag}({\mathit{Z}}^T{\mathit{D}}_m+{\mathit{\mu }}^T{\mathit{D}}_m))}^T+{\mathit{m}}^T{\mathit{D}}_m)$. (2) Update $Z_{k+1}$ as: ${\mathit{Z}}_{k+1}=\mathit{U}\mathrm{diag}(\{max({\sigma }_i-\frac{\lambda }{\rho },0)\}){\mathit{V}}^{\ast }$. (3) Update ${\mathit{\mu }}_{k+1}$ as: ${\mathit{\mu }}_{k+1}={\mathit{\mu }}_k+{\mathit{D}}_m\mathrm{Diag}({\mathit{\alpha }}_{k+1})-{\mathit{Z}}_{k+1}$. Repeat the above procedure until the convergence criterion ${∥\begin{array}{c}{\mathit{D}}_m\mathrm{Diag}(\mathit{\alpha })-\mathit{Z}\end{array}∥}_F\le \epsilon$ is met or the pre-specified number of iterations is reached and get the desired sparse coefficient ${\mathit{\alpha }}^{\ast }$. 5: Output: The corresponding patch of the high resolution image can be reconstructed as $\mathit{l}={\mathit{D}}_l{\mathit{\alpha }}^{\ast }$ and the predicted image $\mathit{L}$ can be obtained by merging all patches.

2.6. Strategy for More Bands

The proposed method considers the structure similarity of different bands and uses pairwise comparisons of the NIR, red and green bands. It should be noted that the relationship between n and m is quadratic ($m=\frac{n^2-n}{2}$), where n and m represent the number of bands and the number of multi-band regularization terms, respectively. When n increases, the model will be much more complicated and difficult to solve.

Table 2. Bandwidth of Landsat and MODIS.

Band Name	Landsat		MODIS
	Band Number	Bandwidth (nm)	Band Number	Bandwidth (nm)
Blue	1	450–520	3	459–479
Green	2	530–610	4	545–565
Red	3	630–690	1	620–670
Near-Infrared	4	780–900	2	841–876
Middle-Infrared	5	1550–1750	6	1628–1652
Middle-Infrared	7	2090–2350	7	2105–2155

shows that adjacent bands have consistent bandwidths. This property indicates that structures of adjacent bands are more similar than those of the other pairs. It is thus reasonable to use adjacent bands constraints instead of pairwise comparisons of all bands. Otherwise, the number of combinations of the adjacent bands will get to be smaller and m will become linear to n ($m=n-1$). Therefore, to efficiently extend it to more bands, we use the strategy that only considers structure similarity of two adjacent bands. This smaller model (AMCFM-s) can be reformulated as follows:

$$\begin{array}{cc}\hfill [{\mathit{\alpha }}_N,{\mathit{\alpha }}_R,{\mathit{\alpha }}_G]=argmin\sum _{c\in \{N,R,G\}}\frac{1}{2}& {∥\begin{array}{c}{\mathit{m}}_c-{\mathit{D}}_{m_c}{\mathit{\alpha }}_c\end{array}∥}_2^2+\lambda {∥\begin{array}{c}{\mathit{D}}_{m_c}{\mathit{\alpha }}_c\end{array}∥}_{\ast }\hfill \\ \hfill +{\tau }_{NR}& {∥\begin{array}{c}S{\mathit{D}}_{l_N}{\mathit{\alpha }}_N-S{\mathit{D}}_{l_R}{\mathit{\alpha }}_R\end{array}∥}_2^2\hfill \\ \hfill +{\tau }_{RG}& {∥\begin{array}{c}S{\mathit{D}}_{l_R}{\mathit{\alpha }}_R-S{\mathit{D}}_{l_G}{\mathit{\alpha }}_G\end{array}∥}_2^2.\hfill \end{array}$$

(17)

The procedure of solving AMCFM-s is the same as that of AMCFM. Details can be found in Section 2.5.

3. Experiments

The performance of our proposed method is compared to those of the four state-of-the-art methods for evaluation. ESTARFM [7] is a weighting method and CRSU [25] is an unmixing-based method. The other two are dictionary learning methods, named SPSTFM [15] and EBSCDL [17].

All programs are run in Windows 10 system (Microsoft, Redmond, Washington, DC, USA) and the processor is Intel Core i7-6700 3.40 GHz (Intel, Santa Clara, CA, USA). All of these fusion algorithms are coded in Matlab 2015a (MathWorks, Natick, MA, USA) except the ESTARFM, which is in IDL 8.5 (Harris Geospatial Solutions, Broomfield, CO, USA).

3.1. Experimental Scheme

In this experiment, we use the data acquired from the Boreal Ecosystem-Atmosphere Study (BOREAS) southern study area on 24 May, 11 July and 12 August in 2001, respectively. The products from Landsat ETM+ and MODIS (MOD09GHK) are selected as the source data for fusion. The Landsat image on 11 July 2001 is set as the target image for prediction. All the data are registered for fine geographic calibration.

In the fusion process, we focus on three bands: NIR, red and green. The size of the test images is 300 × 300. Before the test, we up-sample the MODIS images to the same resolution as the Landsat images via bi-linear interpolation because the spatial resolutions of these two source images are different.

3.2. Parameter Settings and Normalization

The parameters of AMCFM are set as follows. The dictionary size is 256, the patch size is 7 × 7, the overlap of patches is 4, the number of training patches is 2000, $\lambda$ is 0.15, ${\mathit{\alpha }}_0$ is 0, ${\mathit{Z}}_0$ and ${\mathit{\mu }}_0$ are both 0, and $\rho$ is 0.1. All the comparative methods keep their original parameter settings.

Normalization can speed up the computation time and has an effect on the fusion results. As a preprocessing step, the high and low resolution images are normalized as follows:

$$\mathit{L}=\frac{\mathit{L}-\overline{\mathit{L}}}{{\sigma }_{\mathit{L}}},\mathit{M}=\frac{\mathit{M}-\overline{\mathit{M}}}{{\sigma }_{\mathit{M}}},$$

(18)

where $\overline{\mathit{L}}$ is the mean value of image $\mathit{L}$ and ${\sigma }_{\mathit{L}}$ is the standard deviation of image $\mathit{L}$.

3.3. Quality Measurement of the Fusion Results

Several metrics have been used to evaluate the fusion results by different methods. These metrics can be classified into two types, namely the band quality metrics and the global quality metrics.

We employ three assessment metrics, namely the root mean square error (RMSE), average absolute difference (AAD) and correlation coefficient (CC) to assess the performance of the algorithms in each band. The ideal result is 0 for RMSE and AAD, while it is 1 for CC.

Three other metrics are adopted to evaluate the global performance, including relative average spectral error (RASE) [40], Erreur Relative Globale Adimensionnelle de Synthèse (ERGAS) [41] and Q4 [42]. The mean RMSE (mRMSE) of three bands is also used as a global index. The ideal result is 0 for mRMSE, RASE and ERGAS, while it is 1 for Q4. It should be noted that Q4 is defined for four spectral bands. For our comparisons, the real part of a quaternion is set to 0.

3.4. Results

Table 3. Band quality metrics in image 1.

Approaches	RMSE			AAD			CC
	NIR	Red	Green	NIR	Red	Green	NIR	Red	Green
ESTARFM	0.0177	0.0049	0.0043	0.0133	0.0039	0.0033	0.9553	0.9102	0.8951
CRSU	0.0286	0.0099	0.0071	0.0217	0.0084	0.0059	0.9050	0.7670	0.7975
SPSTFM	0.0157	0.0051	0.0040	0.0121	0.0040	0.0031	0.9717	0.8935	0.9148
EBSCDL	0.0157	0.0051	0.0040	0.0120	0.0040	0.0031	0.9721	0.8938	0.9152
AMCFM	0.0152	0.0050	0.0039	0.0118	0.0039	0.0030	0.9755	0.8979	0.9163
AMCFM-s	0.0154	0.0050	0.0039	0.0119	0.0039	0.0030	0.9743	0.8978	0.9158

RMSE, root mean square error; AAD, average absolute difference; CC, correlation coefficient.

,

Table 4. Band quality metrics in image 2.

Approaches	RMSE			AAD			CC
	NIR	Red	Green	NIR	Red	Green	NIR	Red	Green
ESTARFM	0.0201	0.0109	0.0102	0.0133	0.0046	0.0045	0.8678	0.6488	0.5766
CRSU	0.0284	0.0120	0.0108	0.0215	0.0065	0.0054	0.7330	0.5436	0.4933
SPSTFM	0.0193	0.0110	0.0102	0.0125	0.0051	0.0044	0.8696	0.6508	0.5736
EBSCDL	0.0193	0.0110	0.0102	0.0125	0.0051	0.0044	0.8704	0.6515	0.5746
AMCFM	0.0189	0.0110	0.0102	0.0122	0.0051	0.0044	0.8756	0.6525	0.5730
AMCFM-s	0.0188	0.0110	0.0102	0.0121	0.0051	0.0044	0.8789	0.6509	0.5706

and

Table 5. Band quality metrics in image 3.

Approaches	RMSE			AAD			CC
	NIR	Red	Green	NIR	Red	Green	NIR	Red	Green
ESTARFM	0.0194	0.0048	0.0049	0.0136	0.0036	0.0035	0.9136	0.8998	0.8366
CRSU	0.0285	0.0093	0.0064	0.0224	0.0073	0.0044	0.7965	0.6968	0.7251
SPSTFM	0.0165	0.0059	0.0040	0.0122	0.0044	0.0030	0.9276	0.8480	0.8859
EBSCDL	0.0164	0.0059	0.0040	0.0121	0.0044	0.0030	0.9280	0.8483	0.8861
AMCFM	0.0161	0.0057	0.0039	0.0120	0.0042	0.0029	0.9330	0.8583	0.8883
AMCFM-s	0.0162	0.0057	0.0040	0.0121	0.0042	0.0029	0.9324	0.8584	0.8878

show the digital values of these methods in each band. All these methods can reconstruct the target high resolution image. ESTARFM has a good performance in the red band. CRSU performs well in the red band and green band of image 2, but, in most cases, this method has undesirable results. SPSTFM and EBSCDL have similar results and EBSCDL produces slightly higher quality in these three images. AMCFM and AMCFM-s produce the best results for NIR band. Moreover, AMCFM has the best or the second best results in almost all metrics, showing the stability and efficiency in its performance.

The global metrics of different methods are shown in

Table 6. Global quality metrics in image 1.

Approaches	mRMSE	RASE	ERGAS	Q4
ESTARFM	0.0090	10.6579	0.6155	0.9515
CRSU	0.0152	17.5449	1.1343	0.8670
SPSTFM	0.0083	9.5958	0.6013	0.9628
EBSCDL	0.0083	9.5662	0.6002	0.9631
AMCFM	0.0080	9.2880	0.5866	0.9656
AMCFM-s	0.0081	9.4056	0.5884	0.9641

RASE, relative average spectral error; ERGAS, Erreur Relative Globale Adimensionnelle de Synthèse.

,

Table 7. Global quality metrics in image 2.

Approaches	mRMSE	RASE	ERGAS	Q4
ESTARFM	0.0137	14.7738	1.3645	0.8251
CRSU	0.0171	19.3090	1.5131	0.6519
SPSTFM	0.0135	14.4394	1.3702	0.8074
EBSCDL	0.0135	14.4177	1.3686	0.8080
AMCFM	0.0134	14.2524	1.3644	0.8113
AMCFM-s	0.0133	14.1951	1.3671	0.8116

and

Table 8. Global quality metrics in image 3.

Approaches	mRMSE	RASE	ERGAS	Q4
ESTARFM	0.0097	11.2841	0.6511	0.9025
CRSU	0.0147	16.8017	1.0984	0.7556
SPSTFM	0.0088	9.8227	0.6849	0.9117
EBSCDL	0.0088	9.8013	0.6840	0.9120
AMCFM	0.0086	9.5761	0.6652	0.9159
AMCFM-s	0.0086	9.6288	0.6658	0.9148

. AMCFM has the best global performance in all three images, except for Q4 in image 2 and ERGAS in image 3. Image 1 is best captured by our proposed model with a noticeable performance in all four metrics. The outstanding performance of AMCFM is attributed to its improved performance in the NIR band.

Figure 2 and Figure 3 compare the target (true) Landsat images with the images predicted by ESTARFM, CRSU, SPSTFM, EBSCDL, AMCFM and AMCFM-s. We use NIR-red-green band as the red-green-blue-band composite to show the images. These images are displayed with an ENVI 5.3 (Harris Geospatial Solutions, Broomfield, Colorado, United States) 2% linear enhancement.

All these fusion algorithms have the capability to reconstruct the main structure and details of the target image. It appears that the colors of the dictionary learning methods are visually more similar to the true Landsat image than the weighting method and unmixing-based method. The details captured by AMCFM are more prominent than those captured by SPSTFM and EBSCDL, which can be observed in the two-times enlarged red box in the images. Overall, our proposed method has the best performance in visualization.

Figure 4, Figure 5 and Figure 6 display the 2D scatter plots of NIR, red and green band of image 1. ESTARFM performs slightly better than the other methods in the red band. This result is consistent with the statistics in Table 3. However, in the NIR and green band, it is obvious that dictionary learning methods outperform the weighting method and unmixing-based method because scatter plots of ESTARFM and CRSU are more dispersed. The scatter plots of our proposed methods, AMCFM and AMCFM-s, are closer to the 1-1 line than the other methods, indicating that using the edge information can actually improve fusion performance, especially in the NIR band. In general, Figure 4, Figure 5 and Figure 6 show that our proposed methods reconstruct images closest to the true Landsat image.

4. Discussion

Although the model performs well in the experiments, there still exists some questions to be discussed. Therefore, more experiments are performed to answer these questions.

4.1. Which Condition Is Better for AMCFM

Table 6, Table 7 and Table 8 show that image 1 best fits our model. This can be explained by the level of details in

Table 9. Level of details in a target image.

	Image 1			Image 2			Image 3
	NIR	Red	Green	NIR	Red	Green	NIR	Red	Green
$\frac{StandardDeviation}{Mean}$	0.2468	0.2777	0.2008	0.1812	0.4087	0.2542	0.1780	0.3042	0.1725

. We employ “$\frac{StandardDeviation}{Mean}$” to represent the level of details of a target image as in [19]. It is clear that image 1 has the highest level of details in NIR band and the most similar levels of details in the three bands. Therefore, more structure similarity information can be captured to improve the fusion results. When there is a large divergence in a certain band, such as the red band in image 2, the results of the dictionary learning methods in this band are unsatisfactory. Under this situation, the ESTARFM performs better in red band.

4.2. Computational Cost

Computational cost is an important factor in practical application.

Table 10. Computational cost (in seconds).

Time	ESTARFM	CRSU	SPSTFM	EBSCDL	AMCFM	AMCFM-s
training	-	-	29	731	29	29
predicting	103	52	8	101	948	945
total	103	52	37	842	977	974

records the running time of all algorithms in image 1. It shows that SPSTFM has the fastest running speed. EBSCDL is time-consuming because the algorithm models the relationship between high and low resolution patches by a mapping function. AMCFM is a little slower than EBSCDL because of the complexity of the ADMM algorithm. However, for the improvement in results obtained, the slightly increased running time is acceptable. To accelerate the computation, an alternative approach can be designed to solve the reconstruction model efficiently, or the program can be coded with Graphics Processing Unit (GPU) support for parallel running in the future work.

4.3. Parameters

The parameters of the multi-band constraints determine the importance of the corresponding terms in the fusion model and the scalar $\gamma$ affects the value of the parameter $\tau$ directly. In order to find a suitable $\gamma$, Figure 7 depicts how Q4 changes with respect to $\gamma$. Q4 is an index which encapsulates both spectral and radiometric measurements of the fusion result. Thus, we choose it to reflect the fusion results. A larger value of Q4 means a better fusion performance. When $\gamma$ is smaller than 10, the performance of AMCFM evidently improves with the increase of $\gamma$. However, Q4 hardly increases when $\gamma$ is larger than 10. Therefore, we set $\gamma$ to 10.

5. Conclusions and Future Work

In this paper, we have proposed a novel dictionary learning fusion model, called AMCFM. This model considers the structure similarity between bands via adaptive multi-band constraints. These constraints essentially enforce the similarities of the edge information across bands in high resolution patches to improve the fusion performance. Moreover, different from existing dictionary learning models which only emphasize on sparsity, we use the nuclear norm as the regularization term to represent both sparsity and correlation. Therefore, our model can reduce the impact of inefficient dictionary pair and improve the representation ability of the dictionary pair. Comparing with four state-of-the-art fusion methods in metrics and visual effects, the experimental results support our proposed model in the improvements of image fusion. Although our model is slower than the other two dictionary learning methods in this empirical analysis because of the complexity of the optimization algorithm, the fusion results obtained from our model are improved indeed. One may wonder whether it is justifiable to achieve a slight improvement on the expense of an increase in computational time. Our argument is that, on a theoretical basis, our model is more reasonable and appealing than SPSTFM and EBSCDL because it capitalizes on the structure information and correlation of dictionaries for image fusion. Such advantages will be more evident when structure similarity increases.

However, there remains some room for improvement. Firstly, the ${\ell }_2$ norm loss term assumes that noise is an i.i.d. Gaussian. We can consider the use of other noise hypotheses, such as i.i.d. Gaussian mixture and non-i.i.d noise structure, to improve the fusion results. Secondly, the computation cost of the proposed method is high because of the complexity of the ADMM algorithm. To reduce the computation time, an alternative approach can be designed to solve the reconstruction model efficiently for practical applications. To analyze hyperspectral data efficiently, dimension reduction methods might need to be incorporated into the fusion process.