A pan-sharpening method based on the ADMM algorithm

Abstract

Pan-sharpening is a method of integrating low-resolution multispectral images with corresponding high-resolution panchromatic images to obtain multispectral images with high spectral and spatial resolution. A novel variational model for pan-sharpening is proposed in this paper. The model is mainly based on three hypotheses: 1) the pan-sharpened image can be linearly represented by the corresponding panchromatic image; 2) the low-resolution multispectral image is down-sampled from the high-resolution multispectral image through the down-sampling operator; and 3) the satellite image has the low-rank property. Three energy components corresponding to these assumptions are integrated into a variational framework to obtain a total energy function. We adopt the alternating direction method of multipliers (ADMM) to optimize the total energy function. The experimental results show that the proposed method performs better than other mainstream methods in spectral and spatial information preserving aspect.

Keywords

pan-sharpening / multispectral image / panchromatic image / variational framework / energy function / ADMM

Introduction

The goal of pan-sharpening is to combine a high-resolution PANchromatic (PAN) image with a corresponding Low-resolution MultiSpectral (LMS) image with precise spectral information (^{Fang et al., 2013}). The pan-sharpened High-resolution MultiSpectral images (HMS) have more detailed information and the same spectral quality as the original LMS data, which is beneficial to target recognition and classification.

Several pan-sharpening methods have been proposed. For example, component-substitution (CS), including Gram-Schmidt (GS) (^{Laben and Brower, 2000}; ^{Lu and Zhang, 2014}), intensity-hue-saturation (IHS) (^{Rahmani et al., 2010}; ^{Mezouar et al., 2011}), Brovey (^{Zhang and Roy, 2016}), and principal component analysis (PCA) (^{Shah et al., 2007}; ^{Shah et al., 2008}), represent PAN image with a linear combination of bands in LMS images. Although these methods can provide pan-sharpening images with high spatial quality, they may suffer from spectral distortion. In addition, methods based on multiresolution analysis (MRA) (^{Park et al., 2001}), such as Wavelet and Contourlet (^{Metwalli et al., 2014}), inject the high-frequency details of PAN images into the spectral information of MS images. The MRA method can reduce spectral distortion to some extent, but will encounter a new problem of spatial degradation (^{Aiazzi et al., 2002}).

A model-based optimization (MBO) framework has emerged to overcome such degradation problems. All energy functional terms are added to a framework to obtain the optimal value of the total energy function, which corresponds to the best fusion quality (^{Möller et al., 2012}). The P+ XS method (^{Ballester et al., 2006}) was proposed in 2006, which is based on two assumptions, i.e., the spatial details are provided by the PAN image, and the up-sampled LMS image is obtained by blurring the desired HMS image. However, it is difficult to estimate an accurate convolution kernel. To overcome this defect, variational wavelet pan-sharpening (VWP) and alternate VWP (AVWP) (^{Moeller et al., 2008}; Möller et al., 2012) were presented by introducing spectral constraints and assuming that the coefficients of LMS/PAN and HMS are approximately equal. Then in 2014, ^{Aly and Sharma (2014)} determined the sampling interval to coincide with the frequency of the resulting image through approximate frequency-domain analysis. In 2015, Chen et al. proposed a convex optimization method based on the assumption that the down-sampled image of HMS is approximately the same as the LMS image.

Recently, machine learning has developed rapidly in the field of image processing, which gives birth to some learning-based pan-sharpening methods. ^{Ding et al. (2014)} characterized the sparse feature of latent HMS image by Bayesian learning model. ^{Masi et al. (2016)} introduced a three-layer network to tackle the pan-sharpening problem. ^{Wei et al. (2017)} proposed a deep residual pan-sharpening neural network (DRPNN) to perform high quality pan-sharpening. ^{Yuan et al. (2018)} proposed a pan-sharpening method based on a multi-scale and multi-depth convolutional neural network (CNN). ^{Yang et al. (2017)} proposed an end-to-end pan-sharpening system which incorporates problem-specific knowledge into the deep learning framework and is trained in the high-pass domain rather than the image domain.

The above methods based on variation and machine learning may have some shortcomings. Most of the models based on variational methods are relatively old, and the methods based on machine learning have some drawbacks, such as insufficient training samples and low time complexity. Therefore, in this paper, a novel model for pan-sharpening is proposed. The model is mainly composed of three energy terms, which correspond to three reasonable assumptions. The first term can preserve the spatial details of PAN images, the second term preserves spectral information as much as possible, and the third term introduces the low-rank property of HMS images. We solve this energy model with the alternating direction method of multipliers (ADMM) (^{Goldstein et al., 2014}). The time cost of obtaining the optimum is acceptable. Overall, we make three contributions:

• We design an efficient pan-sharpening model based on three reasonable assumptions, in which the low-rank term is newly introduced to pan-sharpening realm.

• We use the powerful ADMM algorithm to solve our non-differentiable optimization problem effectively and flexibility.

• We prove that the proposed method obtains better fusion quality than those obtained by other classic methods.

The rest of this paper is arranged as follows: The energy function of three terms are presented in Section 2. The ADMM algorithm is utilized to optimize the total energy function in Section 3. In Section 4, experimental results and analysis are given. Section 5 concludes the paper.

Methodology

Proposed variational model

In this section, we present three energy function terms based on reasonable assumptions, and then we build proposed total energy into a variational framework.

First, we introduce the related symbol definitions. P:{\Omega }^{rw\times rh}\to \mathrm{R} represents a PAN image, \mathit{M}=(M_{\mathrm{1}},\mathrm{...},M_N):{\Omega }^{w\times h\times N}\to \mathrm{R} represents an LMS image, and \mathit{F}=(F_{\mathrm{1}},\mathrm{...},F_N):{\Omega }^{rw\times rh\times N}\to \mathrm{R} represents an HMS image. \Omega \subset \mathrm{R} denotes an bounded and open domain with a Lipschitz boundary, and r is the value of P/M_{\mathrm{1}}. w, h, and N represent width, height and the number of bands of the MS image, respectively. For a random pixel \mathit{x}\in \Omega, M_i(\mathit{x}) is the intensity of M_i at \mathit{x}, where i\in \{\mathrm{1},\mathrm{...},N\}. The Fourier transform is represented by ℱ, and the Fourier inversion transformation is represented by {ℱ}^{-\mathrm{1}}.

Spatial-information preserving term

Since a PAN image almost covers the entire space of the same scene, it contains nearly all the spatial information of the HMS image. Naturally, PAN P can be approximately expressed by a linear combination of all bands in HMS \mathit{F}. As is well known, we can use the gradient measure to characterize an image’s spatial information. Thus, the first assumption is described as

(1)\nabla P=\nabla {\displaystyle {\sum }_{i=\mathrm{1}}^N}{\alpha }_i{F}_i+\delta ,

where {\alpha }_i denotes the mixing coefficients such that {\displaystyle {\sum }_{i=\mathrm{1}}^N}{\alpha }_i=\mathrm{1}, \nabla is the gradient operator, and \delta is the additive noise with zero mean.

To preserve the spatial information of HMS image F and integrate such constraint into a variational framework, we give the data fidelity term as

(2){\mathrm{E}}_{\mathrm{1}}=\left|\right|\nabla P-\nabla {\displaystyle {\sum }_{i=\mathrm{1}}^N}{\alpha }_i{F}_i|{|}_{\mathrm{2}}^{\mathrm{2}},

where {\Vert \cdot \Vert }_{\mathrm{2}} is the L_{\mathrm{2}} norm.

Spectral-information preserving term

Another basic requirement of the pan-sharpening task is accurate spectrum preservation. To keep the spectral information of LMS image well, we assume that each M_i is formed from F_i after low-pass filtering and down-sampling. We formulate the relationship between M_i and F_i as:

(3){\mathit{M}}_i=DH{F}_i+\eta ,i=\mathrm{1},\mathrm{...},N,

where \eta is the additive noise with zero mean, H and D are low-pass filter and down-sampling operator, respectively.

Different from the P+ XS model (^{Ballester et al., 2006}), we abandon the assumption that the up-sampled LMS image is close to the blurred HMS image. In this paper, similar to the work (^{Chen et al., 2015}), we only assume that the blurred HMS image after down-sampling is close to the LMS image to preserve more accurate spectral information. To integrate Eq. (3) into a variational framework, the second hypothesis is proposed, and it is described as

(4)E_{\mathrm{2}}={\displaystyle {\sum }_{i=\mathrm{1}}^N}\left|\right|M_i-DH{F}_i|{|}_{\mathrm{2}}^{\mathrm{2}}=||\mathit{M}-DH\mathit{F}|{|}_{\mathrm{2}}^{\mathrm{2}}.

Low-rank regularization term

Some literature has shown that natural images have low-rank properties (^{Cai et al., 2010}; ^{Candès et al., 2011}; Wang and Zhang, 2012). How about satellite images?

We show the histogram of singular values of 800 collected HMS images in Fig. 1. One can see that most of the singular values tend to zero, which means the satellite images also have low-rank properties. This isn’t a coincidence, since satellite images usually have repeated textures.

Thus we have the following constraints:

(5)E_{\mathrm{3}}\mathrm{=}{\displaystyle {\sum }_{i=\mathrm{1}}^N}\mathrm{rank}(F_i),

Since \mathrm{rank}(\cdot ) is nonconvex and NP-hard, we use the nuclear norm instead. Then, it can be expressed as

(6)E_{\mathrm{3}}=\left|\right|\mathit{F}|{|}_{*}={\displaystyle {\sum }_{i=\mathrm{1}}^N}|\left|F_i\right|{|}_{*}.

Proposed pan-sharpening model

Taking the above analysis into account, we built the total energy function as

\mathrm{E}(\mathit{F})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{E}}_{\mathrm{1}}+{\displaystyle \frac{\sigma }{\mathrm{2}}}{\mathrm{E}}_{\mathrm{2}}+\beta {\mathrm{E}}_{\mathrm{3}}\mathrm{=}

(7){\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\parallel }\nabla P-\nabla {\displaystyle {\sum }_{i=\mathrm{1}}^N}{\alpha }_i{F}_i{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}+{\displaystyle \frac{\sigma }{\mathrm{2}}}\mathit{\parallel }\mathit{M}-DH\mathit{F}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}+\beta \mathit{\parallel }\mathit{F}{\mathit{\parallel }}_{*},

where \sigma and \beta are positive balance parameters. That is, to obtain HMS image with higher spectral quality, we should give \sigma a higher value and vice versa. Thus, the optimization problem of the proposed variational functional, i.e., Eq. (7), can be described as

(8)\mathit{F}=\underset{\mathit{F}}{\arg \min }\mathrm{E}(\mathit{F}).

Numerical algorithm based on ADMM

Here, we present the numerical solution of the proposed model. As introduced in ^{Boyd et al. (2010}), ADMM is a powerful algorithm that is able to solve our proposed optimization problem effectively. The wide applications in many image processing issues, such as image super-resolution and deblurring, have also shown the advantage of the ADMM algorithm.

The ADMM enables us to split the optimization problem into several subproblems. By replacing \mathit{F} by \mathit{C}, H\mathit{F} by {\mathit{B}}_{\mathrm{1}}, D{\mathit{B}}_{\mathrm{1}} by {\mathit{B}}_{\mathrm{2}}, and \mathit{F} by {\mathit{B}}_{\mathrm{3}}, Eq. (7) is rewritten as

\mathrm{E}(\mathit{F})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\parallel }\nabla P-\nabla {\displaystyle \sum _{i=\mathrm{1}}^N}{\alpha }_i{C}_i{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}+{\displaystyle \frac{\sigma }{\mathrm{2}}}\mathit{\parallel }\mathit{M}-{\mathit{B}}_{\mathrm{2}}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}+\beta \mathit{\parallel }{\mathit{B}}_{\mathrm{3}}{\mathit{\parallel }}_{*}

(9)\mathrm{subject}\mathrm{\hspace{0.17em}}\mathrm{to}\{\begin{array}{l}\mathit{C}=\mathit{F}\\ {\mathit{B}}_{\mathrm{1}}=H\mathit{F}\\ {\mathit{B}}_{\mathrm{2}}=D{\mathit{B}}_{\mathrm{1}}\\ {\mathit{B}}_{\mathrm{3}}=\mathit{F}\end{array}.

The augmented Lagrangian function of Eq. (9) is

\mathrm{L}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\parallel }\nabla P-\nabla {\displaystyle \sum _{i=\mathrm{1}}^N}{\alpha }_i{C}_i{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}+{\displaystyle \frac{\sigma }{\mathrm{2}}}\mathit{\parallel }\mathit{M}-{\mathit{B}}_{\mathrm{2}}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}

+\beta \mathit{\parallel }{\mathit{B}}_{\mathrm{3}}{\mathit{\parallel }}_{*}+\langle {\mathit{\Lambda }}_{\mathrm{1}},\mathit{C}-\mathit{F}\rangle +\langle {\mathit{\Lambda }}_{\mathrm{2}},{\mathit{B}}_{\mathrm{1}}-H\mathit{F}\rangle

+\langle {\mathit{\Lambda }}_{\mathrm{3}},{\mathit{B}}_{\mathrm{2}}-D{\mathit{B}}_{\mathrm{1}}\rangle +\langle {\mathit{\Lambda }}_{\mathrm{4}},{\mathit{B}}_{\mathrm{3}}-\mathrm{F}\rangle

+{\displaystyle \frac{\mu }{\mathrm{2}}}\mathit{\parallel }\mathit{C}-\mathit{F}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}+{\displaystyle \frac{\mu }{\mathrm{2}}}\mathit{\parallel }{\mathit{B}}_{\mathrm{1}}-H\mathit{F}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}

(10)+{\displaystyle \frac{\mu }{\mathrm{2}}}\mathit{\parallel }{\mathit{B}}_{\mathrm{2}}-D{\mathit{B}}_{\mathrm{1}}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}+{\displaystyle \frac{\mu }{\mathrm{2}}}\mathit{\parallel }{\mathit{B}}_{\mathrm{3}}-\mathit{F}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}},

where \mu is a weight parameter, {\mathit{\Lambda }}_i(i=\mathrm{1},\mathrm{...},\mathrm{4}) are Lagrange multipliers, and \mathit{C}=\{C_{\mathrm{1}},\mathrm{...},C_{\mathrm{4}}\}. We aim to minimize \mathrm{L} w.r.t. \mathit{F}, {\mathit{C}}_i, {\mathit{B}}_{\mathrm{1}}, {\mathit{B}}_{\mathrm{2}}, and {\mathit{B}}_{\mathrm{3}} and maximize {\mathit{\Lambda }}_i(i=\mathrm{1},\mathrm{...},\mathrm{4}). In the following, we give the iteration solution of each variable, and we denote t as the iteration number.

Updating \mathit{C}

We first fix {\mathit{F}}^t, {\mathit{\Lambda }}_i^t(i=\mathrm{1},\mathrm{...},\mathrm{4}) and C_i^t(j\ne i), and optimize C_i. By setting the partial derivative \partial \mathrm{L}/\partial C_i to zero, we obtain:

{\displaystyle \frac{\partial \mathrm{L}}{\partial C_i}}={\alpha }_i{\nabla }^{\mathrm{T}}(\nabla {\displaystyle \sum _{i=\mathrm{1}}^N}{\alpha }_i{C}_i-\nabla P)

(11)+{\mathit{\Lambda }}_{\mathrm{1}}+\mu (\mathit{C}-\mathit{F})=\mathrm{0},

where the superscript \mathrm{T} denotes transpose of the operator. Then, the solution of C_i is

(12)C_i^{t+\mathrm{1}}={ℱ}^{-\mathrm{1}}\left\{{\displaystyle \frac{{\alpha }_iℱ({\nabla }^{\mathrm{T}}\nabla (P-{\displaystyle {\sum }_{j\ne i}^N}{\alpha }_j{C}_j^t)+\mu F_i^t-{({\mathit{\Lambda }}_{\mathrm{1}})}_i^t)}{{\alpha }_i^{\mathrm{2}}ℱ({\nabla }^T\nabla )+\mu }}\right\},

where {\nabla }^{\mathrm{T}}\nabla is the negative Laplace operator.

Updating \mathit{F}

By fixing {\mathit{C}}^{t+\mathrm{1}}, {\mathit{{\rm B}}}_i^t(i=\mathrm{1},\mathrm{2},\mathrm{3}), and {\mathit{\Lambda }}_i^t(i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}), there are

{\displaystyle \frac{\partial \mathrm{L}}{\partial F}}=-{\mathit{\Lambda }}_{\mathrm{1}}-H^{\mathrm{T}}{\mathit{\Lambda }}_{\mathrm{2}}-{\mathit{\Lambda }}_{\mathrm{4}}-\mu (\mathit{C}-\mathit{F})

(13)-\mu H^{\mathrm{T}}({\mathit{B}}_{\mathrm{1}}-H\mathit{F})-\mu ({\mathit{B}}_{\mathrm{3}}-\mathit{F})=\mathrm{0.}

We use the fast Fourier transform to solve Eq. (13), obtaining:

(14){\mathit{F}}^{t+\mathrm{1}}={ℱ}^{-\mathrm{1}}\left\{{\displaystyle \frac{ℱ({\mathit{\Lambda }}_{\mathrm{1}}^t+H^{\mathrm{T}}{\mathit{\Lambda }}_{\mathrm{2}}^t+{\mathit{\Lambda }}_{\mathrm{4}}^t+\mu ({\mathit{C}}^{t+\mathrm{1}}+H^{\mathrm{T}}{\mathit{B}}_{\mathrm{1}}^t+{\mathit{B}}_{\mathrm{3}}^t))}{\overline{ℱ(H)}ℱ(H)+\mathrm{2}\mu }}\right\}\mathrm{.}

Updating {\mathit{B}}_{\mathrm{1}}

Next, to optimize {\mathit{B}}_{\mathrm{1}}, we fix {\mathit{C}}^{t+\mathrm{1}}, {\mathit{F}}^{t+\mathrm{1}}, {\mathit{{\rm B}}}_i^t(i=\mathrm{1},\mathrm{2},\mathrm{3}), and {\mathit{\Lambda }}_i^t(i=\mathrm{1},\mathrm{...},\mathrm{4}). Similarly, we obtain:

{\displaystyle \frac{\partial \mathrm{L}}{\partial {\mathit{B}}_{\mathrm{1}}}}={\mathit{\Lambda }}_{\mathrm{2}}^t-D^{\mathrm{T}}{\mathit{\Lambda }}_{\mathrm{3}}^t+\mu ({\mathit{B}}_{\mathrm{1}}-H{\mathit{F}}^{t+\mathrm{1}})

(15)-\mu D^{\mathrm{T}}({\mathit{B}}_{\mathrm{2}}^t-D{\mathit{B}}_{\mathrm{1}})=\mathrm{0.}

Then the solution for B_{\mathrm{1}}is

(16){\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}}={\displaystyle \frac{-{\Lambda }_{\mathrm{2}}^t+D^{\mathrm{T}}{\Lambda }_{\mathrm{3}}^t+\mu (H{\mathit{F}}^{t+\mathrm{1}}+D^{\mathrm{T}}{\mathit{B}}_{\mathrm{2}}^t)}{\mu (\mathrm{1}+D^{\mathrm{T}}D)}}.

Updating {\mathit{B}}_{\mathrm{2}}

With {\mathit{C}}^{t+\mathrm{1}}, {\mathit{F}}^{t+\mathrm{1}}, {\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}}, {\mathit{B}}_{\mathrm{3}}^t, and {\mathit{\Lambda }}_i^t(i=\mathrm{1},\mathrm{...},\mathrm{4}) fixed, we optimize {\mathit{B}}_{\mathrm{2}}:

{\displaystyle \frac{\partial \mathrm{L}}{\partial {\mathit{B}}_{\mathrm{2}}}}=-\sigma (\mathit{M}-{\mathit{B}}_{\mathrm{2}})+{\mathit{\Lambda }}_{\mathrm{3}}^t

(17)+\mu ({\mathit{B}}_{\mathrm{2}}-D{\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}})=\mathrm{0.}

Then the solution for {\mathit{B}}_{\mathrm{2}} is

(18){\mathit{B}}_{\mathrm{2}}^{t+\mathrm{1}}={\displaystyle \frac{\sigma \mathit{M}-{\mathit{\Lambda }}_{\mathrm{3}}^t+\mu D{\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}}}{\mu +\sigma }}.

Updating {\mathit{B}}_{\mathrm{3}}

{\mathit{B}}_{\mathrm{3}} has a different solution from above variables. By utilizing the soft thresholding (SVT) (^{Bredies and Lorenz, 2008}), we have

{\mathit{B}}_{\mathrm{3}}^{t+\mathrm{1}}=\min \beta \mathit{\parallel }{\mathit{B}}_{\mathrm{3}}{\mathit{\parallel }}_{*}+\langle {\mathit{\Lambda }}_{\mathrm{4}}^t,{\mathit{B}}_{\mathrm{3}}-{\mathit{F}}^{t+\mathrm{1}}\rangle +{\displaystyle \frac{\mu }{\mathrm{2}}}\mathit{\parallel }{\mathit{B}}_{\mathrm{3}}-{\mathit{F}}^{t+\mathrm{1}}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}

=\min {\displaystyle \frac{\beta }{\mu }}\mathit{\parallel }{\mathit{B}}_{\mathrm{3}}{\mathit{\parallel }}_{*}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\parallel }{\mathit{B}}_{\mathrm{3}}-{\mathit{F}}^{t+\mathrm{1}}+{\displaystyle \frac{{\mathit{\Lambda }}_{\mathrm{4}}^t}{\mu }}{\mathit{\parallel }}_{\mathrm{2}}^{\mathrm{2}}

(19)=\mathrm{SVT}\left({\mathit{F}}^{t+\mathrm{1}}-{\displaystyle \frac{{\Lambda }_{\mathrm{4}}^t}{\mu }},{\displaystyle \frac{\beta }{\mu }}\right).

Updating {\mathit{\Lambda }}_i(i=\mathrm{1},\mathrm{...},\mathrm{4})

We readily obtain the maximum of {\mathit{\Lambda }}_i(i=\mathrm{1},\mathrm{...},\mathrm{4}):

(20){\mathit{\Lambda }}_{\mathrm{1}}^{t+\mathrm{1}}={\mathit{\Lambda }}_{\mathrm{1}}^t+\mu ({\mathit{C}}^{t+\mathrm{1}}-{\mathit{F}}^{t+\mathrm{1}}),

(21){\mathit{\Lambda }}_{\mathrm{2}}^{t+\mathrm{1}}={\mathit{\Lambda }}_{\mathrm{2}}^t+\mu ({\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}}-H{\mathit{F}}^{t+\mathrm{1}}),

(22){\mathit{\Lambda }}_{\mathrm{3}}^{t+\mathrm{1}}={\mathit{\Lambda }}_{\mathrm{3}}^t+\mu ({\mathit{B}}_{\mathrm{2}}^{t+\mathrm{1}}-D{\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}}),

(23){\mathit{\Lambda }}_{\mathrm{4}}^{t+\mathrm{1}}={\mathit{\Lambda }}_{\mathrm{4}}^t+\mu ({\mathit{B}}_{\mathrm{3}}^{t+\mathrm{1}}-{\mathit{F}}^{t+\mathrm{1}}).

Using the above solvers, we summarize the proposed method in Algorithm 1.

Algorithm 1. Proposed Pan-sharpening Algorithm

Input: LMS image \mathit{M}, high resolution PAN image \mathit{P},

Initialize:

{\mathit{F}}^{\mathrm{0}}={\mathit{B}}_{\mathrm{1}}^{\mathrm{0}}={\mathit{B}}_{\mathrm{2}}^{\mathrm{0}}={\mathit{B}}_{\mathrm{3}}^{\mathrm{0}}=\mathrm{0},

{\mathit{\Lambda }}_{\mathrm{1}}^{\mathrm{0}}={\mathit{\Lambda }}_{\mathrm{2}}^{\mathrm{0}}={\mathit{\Lambda }}_{\mathrm{3}}^{\mathrm{0}}={\mathit{\Lambda }}_{\mathrm{4}}^{\mathrm{0}}=\mathrm{1},

\mu >\mathrm{0},\sigma >\mathrm{0},\beta >\mathrm{0},tol>\mathrm{0},

while {\displaystyle \frac{\mathit{\parallel }{F}^{t-\mathrm{1}}-F^t\mathit{\parallel }}{F^t}}\ge tol

For i=\mathrm{1},\mathrm{...},N

C_i^{t+\mathrm{1}}={ℱ}^{-\mathrm{1}}\left\{{\displaystyle \frac{{\alpha }_iℱ({\nabla }^{\mathrm{T}}\nabla (P-{\displaystyle {\sum }_{j\ne i}^N}{\alpha }_j{C}_j^t)+\mu F_i^t-{({\mathrm{\Lambda }}_{\mathrm{1}})}_i^t)}{{\alpha }_i^{\mathrm{2}}ℱ({\nabla }^{\mathrm{T}}\nabla )+\mu }}\right\},

end-for

{\mathit{F}}^{t+\mathrm{1}}={ℱ}^{-\mathrm{1}}\left\{{\displaystyle \frac{ℱ({\mathit{\Lambda }}_{\mathrm{1}}^t+H^{\mathrm{T}}{\mathit{\Lambda }}_{\mathrm{2}}^t+{\mathit{\Lambda }}_{\mathrm{4}}^t+\mu ({\mathit{C}}^{t+\mathrm{1}}+H^{\mathrm{T}}{\mathit{B}}_{\mathrm{1}}^t+{\mathit{B}}_{\mathrm{3}}^t))}{\overline{ℱ(H)}ℱ(H)+\mathrm{2}\mu }}\right\},

{\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}}={\displaystyle \frac{-{\mathit{\Lambda }}_{\mathrm{2}}^t+D^{\mathrm{T}}{\mathit{\Lambda }}_{\mathrm{3}}^t+\mu (H{\mathit{F}}^{t+\mathrm{1}}+D^{\mathrm{T}}{\mathit{B}}_{\mathrm{2}}^t)}{\mu (\mathrm{1}+D^{\mathrm{T}}D)}},

{\mathit{B}}_{\mathrm{2}}^{t+\mathrm{1}}={\displaystyle \frac{(\sigma \mathit{M}-{\mathit{\Lambda }}_{\mathrm{3}}^t+\mu D{\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}})}{\mu +\sigma }},

{\mathit{B}}_{\mathrm{3}}^{t+\mathrm{1}}=\mathrm{SVT}\left({\mathit{F}}^{t+\mathrm{1}}-{\displaystyle \frac{{\mathit{\Lambda }}_{\mathrm{4}}^t}{\mu }},{\displaystyle \frac{\beta }{\mu }}\right),

{\mathit{\Lambda }}_{\mathrm{1}}^{t+\mathrm{1}}={\mathit{\Lambda }}_{\mathrm{1}}^t+\mu ({\mathit{C}}^{t+\mathrm{1}}-{\mathit{F}}^{t+\mathrm{1}}),

{\mathit{\Lambda }}_{\mathrm{2}}^{t+\mathrm{1}}={\mathit{\Lambda }}_{\mathrm{2}}^t+\mu ({\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}}-H{\mathit{F}}^{t+\mathrm{1}}),

{\mathit{\Lambda }}_{\mathrm{3}}^{t+\mathrm{1}}={\mathit{\Lambda }}_{\mathrm{3}}^t+\mu ({\mathit{B}}_{\mathrm{2}}^{t+\mathrm{1}}-D{\mathit{B}}_{\mathrm{1}}^{t+\mathrm{1}}),

{\mathit{\Lambda }}_{\mathrm{4}}^{t+\mathrm{1}}={\mathit{\Lambda }}_{\mathrm{4}}^t+\mu ({\mathit{B}}_{\mathrm{3}}^{t+\mathrm{1}}-{\mathit{F}}^{t+\mathrm{1}}),

end-while

Output: the pan-sharpened image F.

Experiments and analyses

To verify the feasibility and superiority of the proposed method, we use QuikBird satellite datato analyze the experimental result. There are four bands, including red (R), green (G), blue (B), and near-infrared (NIR), contained in the MS images of this satellite. We first discuss the effect of the three parameters in our augmented Lagrangian function (10). Then we compare the proposed method with the current mainstream methods to evaluate the quality of pan-sharpened images. These benchmark methods included wavelet (^{Zhou et al., 1998}; ^{Yocky, 1995}), the P+ XS (^{Ballester et al., 2006}), VWP (^{Moeller et al., 2008}), SIRF (^{Chen et al., 2014}, ²⁰¹⁵), AVWP (^{Moeller et al., 2008}), and PanNet (^{Yang et al., 2017}. We set the parameter tol in the stop criterion to 10⁻³. All our experiments are carried out with MATLAB 2015b running on an Intel CORE™ 2.33 GHz computer with 8 GB RAM.

Evaluation metrics for image fusion

The effectiveness evaluation of pan-sharpening of remote-sensing images is complex. A subjective assessment method has one sidedness and is greatly influenced by the observer’s own factors. Therefore, the objective evaluation of image quality should be carried out at the same time. The objective quantitative evaluation metrics used for reduced-resolution verification in this paper include UIQI (^{Alparone et al., 2007}; Pushparaj and Hegde, 2017), ERGAS (^{Alparone et al., 2007}; Pushparaj and Hegde, 2017), SCC (Pushparaj and Hegde, 2017), Q⁴ (^{Alparone et al., 2004}, ²⁰⁰⁷), RMSE (Pushparaj and Hegde, 2017), RASE (^{Choi, 2006}), QAVE (^{Wang and Bovik, 2002}; ^{Toet and Hogervorst, 2003}), CC (Pushparaj and Hegde, 2017), and PSNR (^{Wang et al., 2009}).

The effect of the three parameters

There are three parameters in our augmented Lagrangian function (10), i.e., \sigma, \beta, and \mu. To find out how to get best pan-sharpening performance by tuning these parameters, we conduct three groups of experiments, each with one parameter changed and the other two fixed to suitable values. We compute the PSNR value of each result and then draw the curve about PSNR versus the three parameters, respectively, which is shown in Fig. 2. From the three figures, one can clearly see that our algorithm isn't very sensitive to the two weighting parameters \sigma and \beta, but a little sensitive to parameter \mu. To get the best pan-sharpening performance, we set \sigma =\mathrm{0.01}, \beta =\mathrm{0.0001}, and \mu =\mathrm{1} in our experiments.

Reduced-resolution verification

We first show results from experiments where we generate LMS images from HMS. We regard the original MS image as HMS, and down-sample the PAN and MS images using the method illustrated in (^{Vivone et al., 2015}) to obtain the input PAN and IMS images.

Visual comparison

From Figs. 3–5, we see that the results of P+ XS (^{Ballester et al., 2006}), VWP, and AVWP (^{Moeller et al., 2008}) are inferior because the whole image are blurred. For instance, the buildings are blurred in Fig. 4, and the car on the road is so blurred that it cannot be recognized in Fig. 5. Many artifacts exist in the results of wavelet (^{Metwalli et al., 2014}; ^{Zhang, 2010}). Although no obvious artifacts can be found in the images generated by SIRF (^{Chen et al., 2014}, ²⁰¹⁵), a closer look shows that the brightness and color of the images tend to change, which is a spectral distortion phenomenon. In addition, white lines and the trees turn dim in Fig. 3, but the color of the buildings is brighter than that of the ground truth in Fig. 4. In Fig. 5, the lower-right corner of the edge turns bright while the edge looks blurry. Besides, the object details of the fused images obtained by PanNet (^{Yang et al., 2017}) are delicate, and the ground objects are clear, but there are also distortion phenomena, such as slight brightening of the image in Fig. 3, some artifacts in the forest area in Fig. 4, and obvious blurring at the building in Fig. 5. In contrast, our method performs better without obvious artifacts and spectral distortion.

Overall, by comparing the pan-sharpened images in Figs. 3–5, one can see that all methods improve the visual effects of input LMS image. However, results produced by benchmarks show a relative lack of features or details. By comparison, the proposed method generates pan-sharpened results with more accurate spatial details and higher-quality spectral information. Therefore, we can conclude that our method is superior to other methods in a visual comparison.

Quantitative analysis

The data in Tables 1–3 correspond to Figs. 3–5, respectively. We highlight the best result in bold and indicate the reference values in the last row of each table.

In Table 1, compared with the other methods, our method is the best. The UIQI and ERGAS values of our method are slightly worse than that of the SIRF and PanNet, but the other index values are better than those of other methods. The second is the SIRF method, which is better than wavelet, VWP, P+ XS, AVWP, and PanNet. The index values of wavelet, VWP, and P+ XS are average, and there are no outstanding values. The index values of VWP are poor, which explains why the images obtained by this method look more ambiguous in Fig. 3. Although PanNet is the best in the ERGAS metric, it is not ideal on other indicators. In Tables 2 and 3, the index values of AVWP and PanNet are relatively poor, which is consistent with the blurriness and color change of the pan-sharpened images in Figs. 4 and 5. The values of wavelet, P+ XS, and VWP are average. Our method has better numerical results than the other seven methods for most indexes. Therefore, our method is superior to other methods in quantitative metrics, which echoes the previous visual analysis.

Efficiency comparison

To evaluate the efficiency of our method among these methods, i.e., P+ XS (^{Ballester et al., 2006}), VWP, AVWP (^{Moeller et al., 2008}) and SIRF (^{Chen et al., 2014}, ²⁰¹⁵), we make a comparison in terms of computational cost. For different sizes of PAN images, the average computational costs of these five methods are listed in Table 4. The number of iterations of all methods is set to 100, which is enough to achieve the convergence. From this table, our method is much faster than other variational pan-sharpening methods, which is enough to prove the efficiency of this method.

Full-resolution verification

Finally, we evaluate the proposed method at full resolution, i.e., use no ground-truth images as references in this part of the verification. Also, because there are no unified and reliable metrics to measure the quality of pan-sharpened images with full resolution, we make only visual comparisons.

In this section, we select only several methods with better results in the above experiment for comparison, including wavelet, P+ XS, AVWP, VWP, PanNet, and SIRF. The pan-sharpened images obtained by these methods are shown in Fig. 6. The images obtained by VWP, AVWP, P+ XS, and PanNet are so blurry that it is impossible to distinguish the details of ground objects, such as the trees in the images. wavelet, which is confirmed with the visual analysis in Fig. 4(d) and Fig. 5(d). SIRF shows details clearly, but it also has some artifacts, especially the road marked in the red square, which is enlarged in Fig. 8. We then take a closer look at the pan-sharpened images from VWP, AVWP, P+ XS, and PanNet, which are shown in Fig. 7. The details of the building can be clearly identified from the pan-sharpened images by VWP and AVWP, but the trees are still very blurry. The pan-sharpened images obtained by P+ XS and PanNet methods are blurrier than those from VWP and AVWP. Although wavelet preserves spatial information well, there are many artifacts. We show the zoomed-in area circled in red for SIRF and our method in Fig. 8, where the spectral distortion and obvious artifacts of results obtained by SIRF are clearly seen.

The performance at full resolution is consistent with that at reduced resolution. We can conclude that our method is the best among these methods because it has no obvious artifacts or spectral distortion. Again, these results demonstrate the merits of the proposed method, i.e., better preservation of the spectral and spatial information.

Discussion

In this study, a novel variational pan-sharpening model based on three reasonable assumptions is proposed. Experiments at reduced- and full-resolution indicate the efficiency and robustness of the proposed method. Generally, there are two advantages of the proposed method:

1) The widely used prior knowledge on the up-sampled LMS image is abandoned owing to the difficulty of estimating the convolution kernel. Instead, we assume that the down-sampled HMS image is close to the LMS image, which can avoid spectral distortion to some degree.

2) We assume that a linear combination of all bands of HMS image is close to the corresponding PAN image in spatial domain because PAN image contains the entire spatial range of the same scene. The gradient operator is exploited to extract the spatial information of images, which ensures the clear details and sharp edges of the pan-sharpened image.

Although the proposed method can obtain high quality pan-sharpened images, there are some shortcomings. For example, the parameters need to be tuned for each data set, and promising results are accompanied by high computational complexity. In our future work, we will try to find some automatic method for tuning parameters. More efficient algorithms for finding the optimal solution will also be studied.

Conclusions

In this paper, a novel pan-sharpening model is proposed, which is based on three assumptions about the relationship between HMS and available LMS and PAN images. The desired HMS image is obtained by optimizing our variational model using ADMM algorithm. To verify the reliability and efficiency of the proposed method, we compare it with five popular pan-sharpening methods by visual judgment and quantitative analysis. Comparisons at full- and reduced-resolution indicate that our method works better than other methods. From experiments, we can conclude that the proposed method is able to effectively construct HMS images with high resolution.

Pan-sharpening is still a growing research direction. Many problems need to be further studied. For future research, we will try to combine our variational approach with hyperspectral analysis techniques, such as demixing, to further improve the quality of the reconstructed image. Remote-sensing images captured by other satellites may be also exploited.

References

[1] Alparone L, Baronti S, Garzelli A, Nencini F (2004). A global quality measurement of pan-sharpened multispectral imagery. IEEE Geosci Remote Sens Lett, 1(4): 313–317

[2] Alparone L, Wald L, Chanussot J, Thomas C, Gamba P, Bruce L M (2007). Comparison of pansharpening algorithms: outcome of the 2006 GRS-S data-fusion contest. IEEE Trans Geosci Remote Sens, 45(10): 3012–3021

[3] Aiazzi B, Alparone L, Baronti S, Garzelli A (2002). Context-driven fusion of high spatial and spectral resolution images based on oversampled multiresolution analysis. IEEE Trans Geosci Remote Sens, 40(10): 2300–2312

[4] Aly H, Sharma J (2014). A regularized model-based optimization framework for pan-sharpening. IEEE Trans Image Process, 23(6): 2596–2608

[5] Ballester C, Caselles V, Igual L, Verdera J, Roug B (2006). A variational model for P+XS image fusion. Int J Comput Vis, 69(1): 43–58

[6] Boyd S, Parikh N, Chu E, Peleato B, Eckstein J (2010). Distributed optimization and statistical learning via the alternating direction method of multipliers. Found Trends Mach Learn, 3(1): 1–122

[7] Bredies K, Lorenz D A (2008). Soft-thresholding. J Fourier Anal Appl, 14(5–6): 813–837

[8] Cai J, Candes E, Shen Z (2010). A singular value thresholding algorithm for matrix completion. SIAM J Optim, 20(4): 1956–1982

[9] Candes E, Li X, Ma Y, Wright J (2011). Robust principal component analysis. J Assoc Comput Mach, 58(3): 1–37

[10] Chen C, Li Y, Liu W, Huang J (2015). SIRF: simultaneous satellite image registration and fusion in a unified framework. IEEE Transactions on Image Processing, A Publication of the IEEE Signal Processing Society, 24(11): 4213

[11] Chen C, Li Y, Liu W, Huang J (2014). Image fusion with local spectral consistency and dynamic gradient sparsity. IEEE Conference on Computer Vision & Pattern Recognition, 2760–2765

[12] Choi M (2006). A new intensity-hue-saturation fusion approach to image fusion with a tradeoff parameter. IEEE Trans Geosci Remote Sens, 44(6): 1672–1682

[13] Ding X, Jiang Y, Huang Y, Paisley J (2014). Pan-sharpening with a Bayesian nonparametric dictionary learning model. In: Proceeding of 17th International Conference of Artificial Intelligence and Statistics, 176–184

[14] Dong W, Shi G, Li X (2013). Nonlocal image restoration with bilateral variance estimation: a low-rank approach. IEEE Trans Image Process, 22(2): 700–711

[15] Fang F, Li F, Shen C, Zhang G (2013). A variational approach for pan-sharpening. IEEE Trans Image Process, 22(7): 2822–2834

[16] Goldstein T, O’Donoghue B, Setzer S, Baraniuk R (2014). Fast alternating direction optimization methods. SIAM J Imaging Sci, 7(3): 225–231

[17] Laben C, Brower B (2000). Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening. Websterny Uspenfieldny. USA: Eastman Kodak Company, 275–299

[18] Lu X, Zhang J (2014). Panchromatic and multispectral images fusion based on modified GS-SWT. Geoscience & Remote Sensing Symposium (IGARSS), 2530–2533

[19] Masi G, Cozzolino D, Verdoliva L, Scarpa G (2016). Pansharpening by convolutional neural networks. Remote Sens, 8(7): 594

[20] Metwalli M, Nasr A, Faragallah O, Rabaie E E, Abbas A (2014). Efficient pan-sharpening of satellite images with the contourlet transform. Int J Remote Sens, 35(5): 1979–2002

[21] Mezouar M E, Taleb N, Kpalma K, Ronsin J (2011). An IHS-based fusion for color distortion reduction and vegetation enhancement in IKONOS imagery. IEEE Trans Geosci Remote Sens, 49(5): 1590–1602

[22] Moller M, Wittman T, Bertozzi A, Burger M (2013). A variational approach for sharpening high dimensional images. SIAM J Imaging Sci, 5(1): 150–178

[23] Moller M, Wittman T, Bertozzi A (2008). Variational wavelet pan-sharpening. IEEE Trans Geosci Remote Sens, 1–9

[24] Park J, Kim K, Yang K (2001). Image fusion using multiresolution analysis. IEEE Int Geosci Remote Se, 2(2): 864–866

[25] Pushparaj J, Hegde A (2016). Evaluation of pan-sharpening methods for spatial and spectral quality. Appl Geomatics, 9(1): 1–12

[26] Rahmani S, Strait M, Merkurjev D, Moeller M, Wittman T (2010). An adaptive HIS pan-sharpening method. IEEE Geosci Remote Sens Lett, 7(4): 746–750

[27] Shah V, Younan N, King R (2008). An efficient pan-sharpening method via a combined adaptive PCA approach and contourlets. IEEE Trans Geosci Remote Sens, 46(5): 1323–1335

[28] Shah V, Younan N, King R (2007). An adaptive PCA-based approach to pan-sharpening. Remote Sens, 6748: 1–9

[29] Toet A, Hogervorst M (2003). Performance comparison of different gray level image fusion schemes through a universal image quality index. In: Kadar I (Ed.), Signal Processing, Sensor Fusion, and Target Recognition XII, SPIE-5096, 552–561

[30] Vivone G, Alparone L, Chanussot J, Mura M, Garzelli A (2015). A critical comparison among pansharpening algorithms. IEEE Trans Geosci Remote Sens, 53(5): 2265–2586

[31] Wang Q, Yu D, Shen Y (2009). An overview of image fusion metrics. IEEE Instrumentation & Measurement Technology Conference, Singapore, 918–923

[32] Wang S, Zhang Z (2012). Colorization by Matrix Completion. In: Proceedings of the Twenty-Sixth AAAI Conference on Artificial Intelligence (AAAI'12), AAAI Press, 1169–1175

[33] Wang Z, Bovik A (2002). A universal image quality index. IEEE Signal Process Lett, 9(3): 81–84

[34] Wei Y, Yuan Q, Shen H, Zhang L (2017). Boosting the accuracy of multispectral image pansharpening by learning a deep residual network. IEEE Geosci Remote Sens Lett, 14(10): 1795–1799

[35] Yang J, Fu X, Hu Y, Huang Y, Ding X, Paisley J (2017). PanNet: A deep network architecture for pan-sharpening. Proceedings of IEEE International Conference on Computer Vision.,1753–1761

[36] Yocky D (1995). Image merging and data fusion by means of the discrete two-dimensional wavelet transform. J Opt Soc Am A Opt Image Sci Vis, 12(9): 1834–1841

[37] Yuan Q, Wei Y, Meng X, Shen H, Zhang L (2018). A multiscale and multidepth convolutional neural network for remote sensing imagery pan-sharpening. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 11(3): 978–989

[38] Zhang B (2010). Study on image fusion based on different fusion rules of wavelet transform. International Conference on Advanced Computer Theory & Engineering, 3(7): 649–653

[39] Zhang H, Roy D (2016). Computationally inexpensive Landsat 8 operational land imager (OLI) pansharpening. Remote Sens, 8(3): 180

[40] Zhou J, Civco D, Silander J (1998). A wavelet transform method to merge Landsat TM and SPOT panchromatic data. Int J Remote Sens, 19(4): 743–757

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature