Enhanced Back-Projection as Postprocessing for Pansharpening

Abstract

Pansharpening is the process of integrating a high spatial resolution panchromatic image with a low spatial resolution multispectral image to obtain a multispectral image with high spatial and spectral resolution. Over the last decade, several algorithms have been developed for pansharpening. In this paper, a technique, called enhanced back-projection (EBP), is introduced and applied as postprocessing on the pansharpening. The proposed EBP first enhances the spatial details of the pansharpening results by histogram matching and high-pass modulation, followed by a back-projection process, which takes into account the modulation transfer function (MTF) of the satellite sensor such that the pansharpening results obey the consistency property. The EBP is validated on four datasets acquired by different satellites and several commonly used pansharpening methods. The pansharpening results achieve substantial improvements by this postprocessing technique, which is widely applicable and requires no modification of existing pansharpening methods.

1. Introduction

Due to physical and technical constraints [1], current optical earth observation satellites such as QuickBird, IKONOS, GeoEye-1, and WorldView-2, instead of providing images with both high spatial and spectral resolutions, can only produce a high spatial resolution panchromatic (PAN) image and a low spatial resolution multispectral (MS) image over the same area. The spatial and spectral resolutions allow identifying structures in the images through the electromagnetic spectrum and geometrical information, and thus are very important to many remote sensing tasks, such as change detection [2], digital soil mapping [3], and visual image analysis [4]. In order to benefit from these two kinds of information, the PAN and MS images can be integrated to get a MS image with both high spatial and spectral resolutions. This integrating process is classically referred to as pansharpening, which can be viewed as a special kind of image fusion [5] or super-resolution, Loncan15.

In the last two decades, many methods have been developed in the literature of pansharpening. They can be divided into two main groups: component substitution (CS) and multiresolution analysis (MRA) [6]. The CS group is usually called spectral class as the details that shall be injected into the upsampled MS bands are extracted from the PAN image by means of a spectral transformation of MS pixels, whereas the MRA group is also named spatial class as the injected details are extracted with respect to a spatial transformation of the PAN image, mostly relying on linear shift-invariant digital filters [7]. Classic methods of the CS or spectral family include the generalized Intensity-Hue-Saturation (GIHS) [8], the principal component analysis (PCA) [9], the Gram–Schmidt (GS) decomposition [10] and its adaptive version (GSA) [11], the partial replacement adaptive component substitution (PRACS) [12], the band-dependent spatial-detail (BDSD) [13], among many others [14,15,16]. Although these methods are easy to implement and can acquire outstanding spatial details, they may suffer from spectral distortion. On the contrary, the MRA or spatial class usually achieves a better spectral quality in the pansharpening result than that of the CS class. The representative methods belonging to the MRA class are high-pass modulation (HPM) [17], which is also named smoothing filter-based intensity modulation (SFIM) [18], generalized Laplacian pyramid (GLP) [19], morphological-based fusion (MF) [20], among many others [21,22]. However, the MRA methods may make the fused results produce ring artifacts, leading to spatial distortion.

There is in fact a tradeoff between the classes of CS and MRA methods in the improvement of spatial quality and the maintain of spectral information. In order to take advantage of the complimentary performances of the two classes, researchers have explored the hybrid methods of the two classes. For example, Núñez et al. [23] have proposed the additive wavelet luminance proportional (AWLP) method by combining the “à trous” wavelet transform with the GIHS method. Liao et al. [24] have proposed the guided filter PCA (GFPCA) method, which applies a guided filter in the PCA domain. Shah et al. [25] have proposed a combined adaptive PCA algorithm based on the discrete contourlet transform. In addition, the performances of the CS and/or MRA methods can also be improved by considering the modulation transfer function (MTF) (The MTF is a function of spatial resolution, and in theory the MTF decreases as the spatial resolution increase [26].) of the instruments [19] or by the elaborate design of injection gains in a context-adaptive manner [27].

In general, the restoration of ideal high resolution MS image from its degraded version is an ill-posed inverse problem. Different from the CS, MRA and hybrid classes, many regularization- and/or optimization-based methods have been proposed to resolve this ill-posed inverse problem [28], such as the total variation [29,30] or sparse regularization [31,32,33,34,35,36]. These methods usually achieve a better pansharpening quality, however, the high computational complexity makes them unsuitable for practical applications. There are still some works which rely on Bayesian framework [37] to solve this problem. Very recently, deep learning methods have been introduced into remote sensing problems [38,39], and several methods have been proposed to address pansharpening [40,41,42]. However, their performances would degrade if there are not enough training images. Detailed surveys of pansharpening methods can be found in [6]. Moreover, pansharpening methods have been recently expanded to sharp the hyperspectral images by either the high resolution PAN or MS images [43,44].

Although the recently developed methods provide practical advantages, it is not quite adequate to solve the spectral distortions originating from the spectral mismatch between the PAN and MS images. To reduce this problem, statistical matching (also known as histogram matching) technique is performed, explicitly or implicitly through the detail injection model. Several works [11,12] have utilized the histogram matching as a preprocessing or necessary step of pansharpening to boost the spectral accuracy, empirically indicating that the histogram matching can reduce the spectral distortions. Despite the effectiveness of histogram matching, it, until recently, has been analyzed theoretically and explicitly in [45,46]. On the other hand, the high-pass modulation (HPM) injection scheme [21,45,47,48] has recently been employed to improve the spatial details of the fusion products. However, these products will always be with a loss of spatial accuracy presented as a contrary to the spatial consistency property based on Wald’s protocol [49], which requires that the fused image should be as identical as possible to the original image once spatially degraded to its original scale. To satisfy this consistency property, the modulation transfer function (MTF)-matched filter [19], estimated by taking the appropriate cutoff frequency of the MS or PAN instruments into account, have been utilized to enhance the spatial resolution of the fused images while preserving the spectral characteristics of the original MS images in the most effective methods [6,13].

Motivated by the effectiveness of the histogram matching, HPM, and MTF-matched filter to improve the quality of the fused product, and the fact that the strategies of integrating these three techniques have not yet been investigated, a postprocessing algorithm, called enhanced back-projection (EBP), is proposed for pansharpening by integrating the three techniques into the back-projection process [50], which is well known and has been used to improve the quality of fused images in the field of traditional image super-resolution [51,52]. Our proposed EBP first enhances the fused results by both the histogram matching and HPM, then iteratively projects the reconstruction error back to tune the quality of fused images. To be spatial consistency, the MTF-matched filters are employed in the EBP to guide the back-projection process. Extensive experiments are carried out on four real datasets acquired by GeoEye-1, WorldView-1, QuickBird and IKONOS sensors at both the reduced and full scales, results of eight benchmark pansharpening methods with postprocessing by the our EBP indicate that the proposed EBP can further improve the fusion quality.

To summarize, the main contribution of this paper are in three-fold:

• A simple yet effective postprocessing framework is proposed for pansharpening. The proposed method is widely applicable and requires no modification of existing pansharpening methods.
• It has been shown that the back-projection process with the MTF-matched filter as the filter kernel can refine high-frequency texture details and make the pansharpening results satisfying the spatial consistency to some extent.
• Extensive experiments on four kinds of datasets have been conducted and show that the pansharpening results achieve substantial improvements by the proposed method.

The remainder of this paper is organized as follows. Section 2 introduces some background about the histogram matching, HPF, and MTF-matched filter. The proposed EBP method and experiments are presented in Section 3 and Section 4, respectively. And conclusions have been drawn in Section 5.

2. Background

The mathematical notation used hereafter is defined in the following. Images are indicated in a vectorial form by bold lowercase, e.g., $\mathit{z}$, with the ith element indicated as $\mathit{z}\left(i\right)$. In this paper, let ${\mathit{x}}_k^{\mathrm{ms}}$ indicate the ith band of the high resolution MS image, ${\mathit{y}}_k^{\mathrm{ms}}$ be its low resolution MS band obtained by spatially blurred and downsampled of ${\mathit{x}}_k^{\mathrm{ms}}$, and ${\tilde{\mathit{y}}}_k^{\mathrm{ms}}$ be the interpolated MS band from ${\mathit{y}}_k^{\mathrm{ms}}$, i.e., the low resolution MS image at the scale of its high resolution one, where the superscript “ms” denotes the multispectral image. The PAN image with high and low resolutions are indicated by ${\mathit{x}}^{\mathrm{pan}}$ and ${\mathit{y}}^{\mathrm{pan}}$, respectively.

Given the observed low-resolution MS band ${\mathit{y}}_k^{\mathrm{ms}}$ and high-resolution PAN image ${\mathit{x}}^{\mathrm{pan}}$, to obtain the high-resolution MS band ${\mathit{x}}_k^{\mathrm{ms}}$, a general formulation is given by

$${\mathit{x}}_k^{\mathrm{ms}}={\tilde{\mathit{y}}}_k^{\mathrm{ms}}+{\mathit{d}}_k^{\mathrm{ms}},k=1,2,\cdots ,n,$$

(1)

where n is the number of multispectral bands and ${\mathit{d}}_k^{\mathrm{ms}}$ are the missing high-resolution spatial details of the kth MS band ${\mathit{x}}_k^{\mathrm{ms}}$. Usually, one assumes the missed spatial details can be extracted from the weighed high-frequency component of the PAN image as

$${\mathit{d}}_k^{\mathrm{ms}}=g_k({\mathit{x}}^{\mathrm{pan}}-{\tilde{\mathit{y}}}^{\mathrm{pan}}),$$

(2)

where $g_k$ is the injection gain corresponding to the kth band and ${\tilde{\mathit{y}}}^{\mathrm{pan}}$ is a low-pass-filtered version of the PAN image ${\mathit{x}}^{\mathrm{pan}}$ and/or a combination of the MS bands (i.e., intensity component $\mathit{I}$). Different ways used to generate the low resolution image ${\tilde{\mathit{y}}}^{\mathrm{pan}}$ (or intensity component $\mathit{I}$) yield different pansharpening methods [6]. Due to the difference between the MS sensors and the PAN sensor, the extracted spatial details ${\mathit{d}}_k^{\mathrm{ms}}$ from the PAN image are easy to be redundant resulting spectral distortion or deficient resulting in insufficient spatial quality of the fusion products. To overcome these drawbacks, the following three techniques are usually adopted in the filed of pansharpening.

2.1. Histogram Matching

Histogram matching (HM) is a process where a time series, image or higher dimension scalar data is adjusted such that its histogram matches that of another reference dataset. A common application of HM is to match the images from two sensors to balance their different responses. In the field of pansharpening, due to the difference between MS and PAN sensor, HM is often used to compensate for the radiometry differences between the MS and PAN images as a relative sensor calibration technique. For example, the high resolution PAN image ${\mathit{x}}^{\mathrm{pan}}$ is preliminarily histogram matched to the intensity component as for the case of CS class. To ensure a global preservation of radiometry and be band-dependent, the PAN image ${\mathit{x}}^{\mathrm{pan}}$ is histogram-matched to each interpolated MS band ${\tilde{\mathit{y}}}_k^{\mathrm{ms}}$ as for the case of MRA class in the following way

$${\mathit{x}}_{h,k}^{\mathrm{pan}}=({\mathit{x}}^{\mathrm{pan}}-{\mu }_{{\mathit{x}}^{\mathrm{pan}}})·\frac{{\sigma }_{{\tilde{\mathit{y}}}_k^{\mathrm{ms}}}}{{\sigma }_{{\tilde{\mathit{y}}}^{\mathrm{pan}}}}+{\mu }_{{\tilde{\mathit{y}}}_k^{\mathrm{ms}}},$$

(3)

where ${\mathit{x}}_{h,k}^{\mathrm{pan}}$ is the histogram-matched PAN image to the kth interpolated MS band ${\tilde{\mathit{y}}}_k^{\mathrm{ms}}$, ${\mu }_{(·)}$ and ${\sigma }_{(·)}$ denote the mean and standard deviation of an image, respectively. By doing this, the histogram-matched PAN ${\mathit{x}}_{h,k}^{\mathrm{pan}}$, once degraded to the low resolution as that of the kth interpolated MS band ${\tilde{\mathit{y}}}_k^{\mathrm{ms}}$, is expected to exhibit the same mean and variance as ${\tilde{\mathit{y}}}_k^{\mathrm{ms}}$. Although HM is not required, many studies [11,12,45,46] have shown that a correct HM is the key to achieve extra performance from established methods and can help to reduce the spectral distortion.

2.2. High-Pass Modulation

The principle of high-pass modulation (HPM) is to transfer the high-frequency components of the high-resolution PAN image to the MS bands with modulation coefficients being the ratio between each low-resolution MS band and the low-resolution PAN [21,47,53], formulated as

$${\mathit{x}}_k^{\mathrm{ms}}={\tilde{\mathit{y}}}_k^{\mathrm{ms}}·\frac{{\mathit{x}}^{\mathrm{pan}}}{{\tilde{\mathit{y}}}^{\mathrm{pan}}},k=1,2,\cdots ,n.$$

(4)

According to above formulation, the HPM assumes that each pixel of the fused MS band is simply proportional to the corresponding high-resolution image at each pixel. This proportionality is a spatially variable injection gain. The HPM is a successful method in image fusion, and recent works in [6,21,47] have confirmed its capability for improving the quality of the fused product. Although the HPM can help to get very appealing visual features, some possible numerical issues could appear due to division operator, e.g., creating some fused pixels with very high values, resulting in spatial distortion.

2.3. MTF-Matched Filter

Characteristics of the imaging instrument carried on a satellite are important to describe the geometry of the landscape. Each imaging instrument features its own spatial and spectral model embodied by the modulation transfer function (MTF), which is the amplitude spectrum of the system point spread function (PSF). The MTF determines the upper limits of image quality. In practice, the MTFs of the imaging instrument are different from each other for both MS and PAN imaging sensors. Thus, a special structure of the MTF should be taken into account when sharpening the MS image to the resolution of the PAN image. However, the exact filter responses of the MTF for each sensor are not provided by the instrument manufactures. Fortunately, the filter gains at Nyquist cutoff frequency of the MTFs are provided by the manufacturers. Using this information, the MTF filters for MS and Pan sensors can be estimated based on the assumption that the frequency response of each filter is approximately Gaussian-shaped and matched to the gain at Nyquist cutoff frequency of the exact MTF filter [19]. By taking use of the MTF-matched filters, several improved methods, such as the BDSD [13], GLP with MTF-matched filter (MTF-GLP) [19], MTF-GLP-HPM [21], and recently proposed methods [22,54,55], have been developed, and proving its effectiveness [47].

3. The Enhanced Back-Projection Applied to Pansharpening

The strategy used for our proposed post-processing method is based on the well-known back-projection algorithm [50,56]. Thus, in this section, the back-projection method is first introduced and then the proposed EBP algorithm is described in details.

3.1. Back-Projection

The generation process of an low-resolution image $\mathit{y}$ can be modeled by a combination of a blur effect and a down-sampling operator to an ideal high-resolution image $\mathit{x}$, i.e., formulated as follows

$${\mathit{y}=(\mathit{x}\ast \mathcal{F})\downarrow }_r,$$

(5)

where $\mathcal{F}$ is the blurring filter caused by the optical imaging system, * is the convolution operator, and ${\downarrow }_r$ is the down-sampling operator with a scaling factor of r (typically 4 in the literature of pansharpening). Based on the above image acquisition model, back-projection (BP) was firstly proposed by Iranni and Peleg in [50,57] to reconstruct the ideal high resolution image $\mathit{x}$. It is a popular method and applied successfully in the field of image super-resolution [58,59,60].

The BP is originally designed for the image super-resolution reconstruction problem with multiple low resolution inputs [57], and then specialized for the case with single input [56]. The main idea of BP is to iteratively compute the reconstruction error $\mathit{e}$ and then fuse it back to tune the estimated high-resolution image $\widehat{\mathit{x}}$ of $\mathit{x}$. By starting with an initial estimation of high-resolution image ${\widehat{\mathit{x}}}^0$, the updating procedure of BP can be summarized as repeatedly doing the following two steps:

$$\begin{array}{cc}\hfill {\mathit{e}}_t& =\mathit{y}-({\widehat{\mathit{x}}}_t\ast \mathcal{F}){\downarrow }_r,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\widehat{\mathit{x}}}_{t+1}& ={\widehat{\mathit{x}}}_t+{\mathit{e}}_t{\uparrow }_r\ast \mathcal{G},\hfill \end{array}$$

(7)

where ${\widehat{\mathit{x}}}_t$ is the estimated high resolution image at the tth iteration, and $\mathcal{G}$ is a back-projection filter, ↑ is the upsampling operator, until the differences between the simulated low resolution image and the input low resolution image $\mathit{y}$ are small enough. The combination of the above two step is referred to as BP algorithm.

It has been shown in the following theorem [50,56] that, under certain conditions, the BP algorithm can converge to the desired deblurred image, which satisfies the image acquisition model Equation (5). Although the BP has been proven to improve the image quality [52,61], because of the ill-conditional nature of the generation process model, it is very sensitivity to (inappropriate) selection of initial high resolution image ${\widehat{\mathit{x}}}_0$ and the filters $\mathcal{F}$ and $\mathcal{G}$, thus leading to failure to converge and variability in results.

Theorem 1 

(see [50,56]). The updating rules of Equations (6) and (7) will converge to a desired image $\widehat{\mathit{x}}$, which satisfies Equation (5), with an exponential rate for all $r\ge 1$, if the following condition holds:

$${\parallel \delta -\mathcal{F}\ast \mathcal{G}\downarrow }_r{\parallel }_2<1,$$

where δ denotes the unity pulse function centered at $(0,0)$.

3.2. EBP as Postprocessing for Pansharpening

Motivation. According to Theorem 1, the result of BP will converge to an image that satisfies the generation model Equation (5), which can also be seen as the model of how the low-resolution MS band ${\mathit{y}}_k^{\mathrm{ms}}$ be generated from the corresponding high-resolution MS band ${\mathit{x}}_k^{\mathrm{ms}}$. This just fulfill the spatial consistency property of Wald’s protocol, which requires that the fused MS band, once spatially degraded to its original scale, should be as identical as possible to the original MS band ${\mathit{y}}_k^{\mathrm{ms}}$ [49,62]. However, the BP is sensitivity to the initial choice of image ${\widehat{\mathit{x}}}_0$ and the blur operators $\mathcal{F}$ and $\mathcal{G}$, usually leading to variability in results and failure in convergence. In the literature of pansharpening, the blur operators $\mathcal{F}$ and $\mathcal{G}$ can be estimated by the MTF-matched filters for each MS band as stated previously, and a natural choice for the initial image ${\widehat{\mathit{x}}}_0$ is the result obtained by other pansharpening methods since it is more closer to the ideal high resolution MS image.

On the other hand, HM can reduce the spectral distortion and HPM can improve the spatial quality, especially when the image to be matched or modulated is as close as possible to the reference. All of these motivate us to exploit them to further improve the fusion results. Therefore, by considering the characteristics of the imaging instrument (i.e., the MTF-matched filter), careful design of using the above three techniques of HM, HPM, and BP as postprocessing steps will significantly improve the quality of the fusion results.

Method. With the aforementioned motivation in mind, by integrating these techniques previously described into the BP method, an enhanced back-projection (EBP) method as postprocessing for pansharpening is proposed. The proposed EBP consists of two stages: enhancement stage and back-projection stage, which are described in detail below.

1.

Enhancement stage-This stage is to further adjust the spectral accuracy by HM and improve the spatial quality by HPM. Given the initial HR MS band ${\widehat{\mathit{x}}}_{k,0}^{\mathrm{ms}},k=1,2,\cdots ,n,$ obtained by another existing method, panchromatic image ${\mathit{x}}^{\mathrm{pan}}$, and the observed low resolution MS band ${\mathit{y}}_k^{\mathrm{ms}}$, $k=1,2,\cdots ,n$, the enhancement stage is formed by the following two consecutive steps:

1.1

Histogram matching: The panchromatic image ${\mathit{x}}^{\mathrm{pan}}$ is histogram-matched to each interpolated low-resolution MS band ${\tilde{\mathit{y}}}_k^{\mathrm{ms}}$ by the following equation

$$\begin{array}{cc}\hfill {\mathit{x}}_{h,k}^{\mathrm{pan}}& =({\mathit{x}}^{\mathrm{pan}}-{\mu }_{{\mathit{x}}^{\mathrm{pan}}})·\frac{{\sigma }_{{\tilde{\mathit{y}}}_k^{\mathrm{ms}}}}{{\sigma }_{{\mathit{x}}^{\mathrm{pan}}}}+{\mu }_{{\tilde{\mathit{y}}}_k^{\mathrm{ms}}},k=1,2,\cdots ,n.\hfill \end{array}$$

(8)

1.2

High-pass modulation: Each initial HR MS band ${\widehat{\mathit{x}}}_{k,0}^{\mathrm{ms}}$ is modulated by the corresponding histogram-matched PAN image ${\mathit{x}}_{h,k}^{\mathrm{pan}}$ as follows

$${\widehat{\mathit{x}}}_{k,1}^{\mathrm{ms}}={\widehat{\mathit{x}}}_{k,0}^{\mathrm{ms}}·\frac{{\mathit{x}}_{h,k}^{\mathrm{pan}}}{{\tilde{\mathit{y}}}_{h,k}^{\mathrm{pan}}},$$

(9)

where ${\tilde{\mathit{y}}}_{h,k}^{\mathrm{pan}}$ is an low-resolution version of panchromatic image obtained by low-pass filtering the ${\mathit{x}}_{h,k}^{\mathrm{pan}}$.

2.

Back-projection stage-This stage is to tune the spatial details injected into the estimated high resolution MS band ${\widehat{\mathit{x}}}_{k,1}^{\mathrm{ms}}$. Although the HPM can significantly improve the injected high-resolution spatial details, the fusion results usually sharpen too much to comply with the spatial consistency property. Therefore, in order to assure the fusion results to be consistent with the low resolution MS band, for each MS band ${\widehat{\mathit{x}}}_{k,1}^{\mathrm{ms}}$ iteratively doing the following two steps:

2.1

Compute the tth iteration reconstruction error ${\mathit{e}}_{k,t}^{\mathrm{ms}}$ between the original low-resolution MS band ${\mathit{y}}_k^{\mathrm{ms}}$ and the the back-projected low-resolution MS band as follows

$${\mathit{e}}_{k,t}^{\mathrm{ms}}={\mathit{y}}_k^{\mathrm{ms}}-({\widehat{\mathit{x}}}_{k,t}^{\mathrm{ms}}\ast {\mathcal{F}}_k^{\mathrm{MTF}}){\downarrow }_r$$

(10)

where ${\mathcal{F}}_k^{\mathrm{MTF}}$ is a corresponding MTF-matched filter to the kth MS band.

2.2

Back-project the error ${\mathit{e}}_{k,t}^{\mathrm{ms}}$ to adjust the fused MS band ${\widehat{\mathit{x}}}_{k,t}^{\mathrm{ms}}$ as follows

$${\widehat{\mathit{x}}}_{k,t+1}^{\mathrm{ms}}={\widehat{\mathit{x}}}_{k,t}^{\mathrm{ms}}+{\mathit{e}}_{k,t}^{\mathrm{ms}}{\uparrow }_r\ast {\mathcal{G}}_k^{\mathrm{MTF}}$$

(11)

where ${\mathcal{G}}_k^{\mathrm{MTF}}$ is a back-projection filter corresponding to the MTF-matched filter ${\mathcal{F}}_k^{\mathrm{MTF}}$.

With this approach, the following two main aspects should be pointed out.

• Currently, pansharpening postprocessing has not received sufficient attention. To the best of our knowledge, the integration of the above three techniques for pansharpening postprocessing is lacking. It will be shown that the proposed EBP as a refinement of the fusion results obtained by existing methods can enhance the quality of fusion product.
• It should be pointed out that the aim of a HM procedure (i.e., Equation (8)) is to obtain a PAN with the same mean of a MS band, but not with the same standard deviation. This is because the PAN is a high spatial resolution image with a large quality of details at high resolution, whereas the MS bands are low-resolution images without details at high resolution. Therefore, the standard deviation is higher than the one of MS band, and the normalization of HM usually has to be made with respect to the standard deviation of a low-pass version of the PAN, not the one of the PAN itself. Here the high resolution PAN is used in the PM procedure because this paper aims at the postprocessing of the high resolution MS bands and this will results in better scores in the experiments.
• The proposed EBP has been designed as a postprocessing approach, and hence, it does not require any modifications for the existing pansharpening methods. Additionally, promising results can be obtained by the proposed EBP efficiently.

4. Experimental Results and Analysis

In this section, several experiments on four kinds of real datasets at both reduced and full scales are conducted to demonstrate the effectiveness of our proposed postprocessing method (i.e., EBP) in enhancing the performances of other pansharpening methods. And the following eight benchmark pansharpening methods are selected to evaluate the proposed EBP:

• BDSD [13], which obtains the optimal minimum mean square error (MMSE) joint estimation of the spectral weights and injection gains at a reduced resolution by using the MTF of the MS sensor.
• GFPCA [24], which is a hybrid method of the CS and MRA classes by applying the guided filter in the PCA domain.
• GSA [11], which is an improved version of GS [10] to capture the spectral responses of sensors by optimizing the mean square error (MSE) with respect to the PAN image.
• MF [20], which is based on the nonlinear decomposition scheme of morphological operators.
• Nonlinear IHS (NLIHS) [63], which estimates the intensity component via local and global synthesis approaches.
• PRACS [12], which generates high resolution details by partial replacement and uses statistical ratio-based injection gains.
• SFIM [18], which is base on the idea of using the ratio between the high resolution PAN image and its low resolution version obtained by low-pass filtering.
• CNN-based Pansharpening (PNN) [41], which adapts a three-layer convolutional neural networks (CNN) to perform the pansharpening task. Note that its results for QuickBird dataset are not reported since the trained model for QuickBird sensor is not provided.

All of the parameters for the above eight methods are set in accordance with the authors’ statements in their papers. It should be pointed out that the results of these eight methods reported in the following experiments have been obtained by the public available toolbox or the source codes kindly provided by the original authors. For example, the implementations of BDSD, GSA, PRACS and SFIM is from the pansharpening Matlab toolbox (Available online: https://openremotesensing.net/knowledgebase/a-critical-comparison-among-pansharpening-algorithms/ (accessed on 8 May 2012), Vivone15), and the source codes of NIHS and MF can also be found at https://openremotesensing.net/, the codes of GFPCA and PNN are provided by the original authors.

4.1. Datasets

The fusion results are evaluated on several datasets acquired by four different satellites. These datasets cover a variety of scenes. The parameters of the four satellites (IKONOS, QuickBird, WorldView-2, GeoEye-1) are reported in

Table 1. Parameters of different satellites. Note that the gains at Nyquist cutoff frequency of the MTFs (MTF gain) are reported in parentheses and in blue.

Parameters		IKONOS	QuickBird	WorldView-2	GeoEye-1
Launch date		24 September 1999	18 October 2001	8 October 2009	6 September 2008
Temporal resolution		<3 days	1–5 days	1.1 day	<3 days
Radiometric resolution		11	11	11	11
Spatial resolution	MS	4 m	2.4 m	1.84 m	1.84 m
	PAN	1 m	0.6 m	0.46 m	0.46 m
Spectral range (MTF gain)	Blue	445–516 nm (0.27)	450–520 nm (0.34)	450–510 nm (0.35)	450–510 nm (0.23)
	Green	506–595 nm (0.28)	520–600 nm (0.32)	510–580 nm (0.35)	510–580 nm (0.23)
	Red	632–698 nm (0.29)	630–690 nm (0.30)	630–690 nm (0.35)	655–690 nm (0.23)
	NIR1	757–900 nm (0.28)	760–900 nm (0.22)	770–895 nm (0.35)	780–920 nm (0.23)
	Red edge			705–745 nm (0.35)
	Coastal			400–450 nm (0.35)
	Yellow			585–625 nm (0.27)
	NIR2			860–1040 nm (0.35)
	PAN	450–900 nm (0.17)	450–900 nm (0.15)	450–800 nm (0.11)	450–800 nm (0.16)

, where the MTF gain for each band is also shown in parentheses and in blue. And the details of the datasets are described below.

• IKONOS Data Set: This dataset is composed of a pair of MS and PAN images, which are acquired by the IKONOS satellite in Sichuan, China, on 16 November 2000 and can be downloaded from http://glcf.umiacs.umd.edu. The IKONOS satellite can produce a PAN image with 1-m spatial resolution and MS images with 4-m spatial resolution. Each MS image has four different bands, namely blue, green, red, and nearinfrared (NIR). This test site contains abundant objects such as mountainous and farmland, roads, and some houses after an earthquake.
• QuickBird Data Set: The second dataset is collected by the QuickBird satellite on an areas of Shatin, Hong Kong, on 7 January 2007. Similar to the IKONOS dataset, the QuickBird dataset also has the MS image with four bands (blue, green, red and NIR) and a PAN image, and with the spatial resolution of 0.6-m for the MS images and of 2.4-m for the PAN image. The test scene covers a number of large buildings such as skyscrapers, commercial and industrial structures, and a number of small objects such as cars, small housings, a playground and so on.
• WorldView-2 Data Set: This dataset has been acquired by the WorldView-2 satellite on 3 April 2011 and can be downloaded from http://cms.mapmart.com/Samples.aspx. The WorldView-2 satellite was launched on 8 October 2009. Different from the above two satellites, the WorldView-2 satellite can provide the MS images with 8 bands, including 4 standard bands (red, green, blue, and NIR1) and 4 new bands (coastal, yellow, red edge, and NIR2), refer to Table 1 for more details. And it produces 0.46-m spatial resolution for the PAN image and 1.84-m spatial resolution for the MS images. The land cover types of the test PAN and MS images for this dataset are mainly some buildings with shadows and some trees.
• GeoEye-1 Data Set: The last dataset is provided by the GeoEye-1 satellite, which is capable of acquiring data at 0.41-m for PAN and 1.65-m for the MS images. Similar to the IKONOS and QuickBird imagery, the GeoEye-1 imagery is also composed of four bands covering visible and near-infrared for the MS images. This test site contains both homogeneous and heterogeneous areas with a lot of fine spatial details.

In our experiments, the PAN images are of size $512\times 512$ and the MS images are of size $128\times 128$.

4.2. Quality Evaluation

The quality of the pansharpening results can be evaluated subjectively and/or objectively. For subjective evaluation, a visual analysis on a color display of the fused image is often performed to see whether the fused objects are more clear than the original MS image and their colors are natural and similar to those of their low resolution ones. While the objective evaluation is a challenging problem since the reference image (i.e., the ideal high-resolution MS images ${\mathit{x}}_k^{\mathrm{ms}}$, $k=1,2,\cdots ,n$) are not available in practice. Currently, there are two ways to quantitatively measure the fusion results. One is Wald’s protocol [49], which is based on the assumption of scale invariance, i.e., the quantitative evaluation is performed on a reduced scale such that the original MS image can be used to as a reference and to compare with the pansharpenend MS image. Another is to make the quantitative evaluation without a reference at a full scale. In this case, the quantitative metrics are designed by exploiting both the relationship between the pansharpened MS images and the original MS images and the relationship between the pansharpened MS image and the original PAN image [64].

Over the past decays, a lot of quantitative metrics were developed to measure the results based on Wald’s protocol. In our experiments, four metrics are used to quantitatively evaluate the performance of the proposed EBP on the above four datasets. They are defined as follows:

• Correlation Coefficient (CC) [5]

  $$\mathrm{CC}=\frac{{\sum }_{i=1}^p\left({\mathit{x}}^{\mathrm{ms}}\left(i\right)-{\mu }_{{\mathit{x}}^{\mathrm{ms}}}\right)\left({\widehat{\mathit{x}}}^{\mathrm{ms}}\left(i\right)-{\mu }_{{\widehat{\mathit{x}}}^{\mathrm{ms}}}\right)}{\sqrt{{\sum }_{i=1}^p{\left({\mathit{x}}^{\mathrm{ms}}\left(i\right)-{\mu }_{{\mathit{x}}^{\mathrm{ms}}}\right)}^2{\sum }_{i=1}^p{\left({\widehat{\mathit{x}}}^{\mathrm{ms}}\left(i\right)-{\mu }_{{\widehat{\mathit{x}}}^{\mathrm{ms}}}\right)}^2}},$$

  where ${\mathit{x}}^{\mathrm{ms}}$ and ${\widehat{\mathit{x}}}^{\mathrm{ms}}$ are the reference and fused MS images with size of p pixels.
• Root Mean Square Error (RMSE) [65] is calculated for the kth MS band as follows

  $$\mathrm{RMSE}\left(k\right)=\sqrt{\frac{1}{p}\sum _{i=1}^p{({\mathit{x}}_k^{\mathrm{ms}}\left(i\right)-{\widehat{\mathit{x}}}_k^{\mathrm{ms}}\left(i\right))}^2}.$$

  (12)
• Erreur Relative Globale Adimensionnelle de Synthése (ERGAS, or relative dimensionless global error in synthesis) [66] is defined as

  $$\mathrm{ERGAS}=\frac{100}{\beta }\sqrt{\frac{1}{n}\sum _{k=1}^n{\left(\frac{\mathrm{RMSE}\left(k\right)}{{\mu }_{{\mathit{x}}_k^{\mathrm{ms}}}}\right)}^2},$$

  (13)

  where n is the number of bands, $\beta$ is the scale ratio between the PAN and the original MS images and ${\mu }_{{\mathit{x}}_k^{\mathrm{ms}}}$ is the mean of the kth reference MS band ${\mathit{x}}_k^{\mathrm{ms}}$.
• Spectral Angle Mapper (SAM) [67] between two spectral vectors $\mathit{x}$ and $\widehat{\mathit{x}}$ is defined as

  $$\mathrm{SAM}(\mathit{x},\widehat{\mathit{x}})=arccos\left(\frac{\langle \mathit{x},\widehat{\mathit{x}}\rangle }{{\parallel \mathit{x}\parallel }_2{\parallel \widehat{\mathit{x}}\parallel }_2}\right),$$

  (14)

  where $\langle ·,·\rangle$ denotes the inner product and ${\parallel ·\parallel }_2$ denotes the $l_2$-norm.
• Q4/Q8 [68,69] is an extension of the universal image quality index (UIQI) [70], and Q4 is given by

  $$Q4=\frac{{\sigma }_{{\mathit{z}}_1{\mathit{z}}_2}}{{\sigma }_{{\mathit{z}}_1}{\sigma }_{{\mathit{z}}_2}}·\frac{2|{\mu }_{{\mathit{z}}_1}\left|\right|{\mu }_{{\mathit{z}}_1}|}{|{\mu }_{{\mathit{z}}_1}{|}^2+\left|{\mu }_{{\mathit{z}}_2}\right|}·\frac{2{\sigma }_{{\mathit{z}}_1}{\sigma }_{{\mathit{z}}_2}}{{\sigma }_{{\mathit{z}}_1}^2+{\sigma }_{{\mathit{z}}_2}^2},$$

  (15)

  where ${\mathit{z}}_1={\mathit{x}}_1^{\mathrm{ms}}+\mathrm{i}{\mathit{x}}_2^{\mathrm{ms}}+\mathrm{j}{\mathit{x}}_3^{\mathrm{ms}}+\mathrm{k}{\mathit{x}}_4^{\mathrm{ms}}$, ${\mathit{z}}_2={\widehat{\mathit{x}}}_1^{\mathrm{ms}}+\mathrm{i}{\widehat{\mathit{x}}}_2^{\mathrm{ms}}+\mathrm{j}{\widehat{\mathit{x}}}_3^{\mathrm{ms}}+\mathrm{k}{\widehat{\mathit{x}}}_4^{\mathrm{ms}}$, ${\mathit{x}}_k^{\mathrm{ms}}$ and ${\overline{\mathit{x}}}_k^{\mathrm{ms}}$ are the kth band of the reference and fused MS images, respectively. Here, $\mathrm{i},\mathrm{j}$ and k are imaginary units, ${\mu }_{\mathit{z}}$ and ${\sigma }_{\mathit{z}}$ are the mean and variance of variable $\mathit{z}$, and ${\sigma }_{{\mathit{z}}_1{\mathit{z}}_2}$ is the covariance between ${\mathit{z}}_1$ and ${\mathit{z}}_2$. $Q4$ is usually calculated using a sliding window $r\times r$ (typically $16\times 16$ or $32\times 32$) and averaged on the entire image. Q4 has been extended to Q8 index such that it is suitable for images whose number of bands is any power of two, refer to [69] for more details.

Among them, CC, RMSE, and SAM are usually averaged over all of the MS bands to yield an overall score. The closer to 1 the values of CC and Q4 are, the better are the quality of the pansharpened MS images. For SAM, RMSE and ERGAS, the ideal values are 0.

Additionally, the Quality with No Reference (QNR) metric proposed by Alparone et al. in [64] is applied to perform the quantitative assessment at a full scale. The QNR metric is defined as follows

$$\mathrm{QNR}=(1-D_{\lambda })(1-D_s),$$

(16)

where $D_{\lambda }$ is a spectral distortion metric, given by

$$D_{\lambda }=\frac{1}{n(n-1)}\sum _{k=1}^n\sum _{l=1}^n|Q({\mathit{y}}_k^{\mathrm{ms}},{\mathit{y}}_l^{\mathrm{ms}})-Q({\widehat{\mathit{x}}}_k^{\mathrm{ms}},{\widehat{\mathit{x}}}_l^{\mathrm{ms}})|,$$

(17)

and $D_s$ is a spatial quality metric, defined by

$$D_s=\frac{1}{n}\sum _{k=1}^n|Q({\widehat{\mathit{x}}}_k^{\mathrm{ms}},{\mathit{x}}^{\mathrm{pan}})-Q({\mathit{y}}_k^{\mathrm{ms}},{\mathit{y}}^{\mathrm{pan}})|.$$

(18)

Here Q denotes the universal image quality index (UIQI) [70]. The optimal values of $D_{\lambda }$ and $D_s$ are 0, and thus the closer to 1 the value of QNR is, the better is the quality of the fused product. Note that the implementations of the quality indices used in the experiments are also from the above pansharpening Matlab toolbox.

4.3. Result Analysis

Experiments were carried out at both reduced and full scales. As for the reduced scale experiments, one can follow Wald’s protocol [49] where the original MS and PAN images are degraded by a low-pass Gaussian filters with the gains at the Nyquist frequency (as reported in Table 1) being the same as those of the sensors. As stated above, the results are evaluated by visual analysis and quantitative measures, respectively.

Visual Analysis. The pansharpened images on all the four datasets for the eight methods with and without our proposed EBP as postprocessing are shown in Figure 1, Figure 2, Figure 3 and Figure 4, where the RGB bands are displayed for visual comparison. It is worth noting that the conclusions are similar for both reduced and full scales. For the sake of brevity and the limitation of space, we don’t show the visual results at the reduced scale.

All the eight methods can effectively sharp the expanded MS image by bicubic interpolation as in their paper [11,12,13,18,19,20,24,35,63]. However, as one can see from these figures, the pansharpened images obtained by the eight methods with our proposed EBP as postprocessing, as shown on the right side of each subfigure, have better spatial and color quality than those of methods without EBP, as shown on the left side of each sugfigure. The proposed EBP can significantly improve the spatial details of the fused images obtained by other methods, this is clearly visible on the edges of the buildings, especially for GFPCA and NLIHS, which produce fused images with very strong blurring for the mountains in Figure 1 and the buildings in Figure 2 and Figure 3, as well as on the houses and ground surface in Figure 4. Although some methods such as BDSD can produce the better results with less blurring, they sometimes suffer from strong spectral distortions, as shown in Figure 4, the false RGB colors indicate that the results after postprocessing by our proposed EBP can preserve well the spectral information, as also confirmed by the following quantitative comparisons. As for the GSA, MF, PNN methods, although they have a good balance between the injected spatial details and the maintain of original spectral information, their results with the EBP postprocessing can obtain more spatial details and more pleasant colors that those of without EBP postprocessing, see the displays in Figure 1, Figure 2, Figure 3 and Figure 4 for comparisons.

As a summary of the visual analysis, the proposed EBP algorithm can significantly enhance the performances of other methods for improving the spatial details while presenting natural and pleasant visual results.

Quantitative Analysis. Quantitative evaluation further helps explaining the conclusions drawn from the visual analysis. The quantitative results by the considered metrics are provided as an objective evaluation. Specially,

Table 2. Quantitative comparison for the eight methods without (denoted by ✗ ) and with (denoted by ✓) postprocessing by our proposed EBP on four kinds of satellite datasets in the full scale. Note that, for better comparison, the results improved by EBP are highlighted in blue, otherwise the original results are colored in green.

Datasets	Metrics	EBP	BDSD	GFPCA	GSA	MF	NLIHS	PRACS	SFIM	PNN
IKONOS	$D_{\lambda }$	✗	0.010	0.020	0.054	0.041	0.028	0.021	0.034	0.082
		✓	0.010	0.041	0.030	0.022	0.012	0.028	0.021	0.020
	$D_s$	✗	0.066	0.129	0.454	0.244	0.103	0.236	0.143	0.242
		✓	0.062	0.056	0.171	0.152	0.116	0.119	0.141	0.123
	QNR	✗	0.925	0.854	0.517	0.724	0.872	0.748	0.828	0.696
		✓	0.928	0.906	0.804	0.830	0.873	0.856	0.841	0.859
QuickBird	$D_{\lambda }$	✗	0.039	0.035	0.031	0.035	0.006	0.037	0.021	-
		✓	0.019	0.020	0.013	0.010	0.011	0.025	0.009	-
	$D_s$	✗	0.022	0.079	0.065	0.036	0.060	0.061	0.025	-
		✓	0.020	0.024	0.018	0.027	0.022	0.019	0.020	-
	QNR	✗	0.940	0.889	0.907	0.930	0.935	0.904	0.954	-
		✓	0.962	0.957	0.969	0.962	0.967	0.956	0.972	-
WorldView-2	$D_{\lambda }$	✗	0.033	0.037	0.009	0.015	0.009	0.013	0.009	0.028
		✓	0.006	0.010	0.008	0.007	0.009	0.018	0.005	0.013
	$D_s$	✗	0.092	0.048	0.044	0.017	0.038	0.039	0.026	0.030
		✓	0.045	0.012	0.009	0.015	0.004	0.009	0.006	0.011
	QNR	✗	0.878	0.916	0.947	0.969	0.953	0.948	0.965	0.943
		✓	0.949	0.978	0.982	0.979	0.987	0.973	0.989	0.976
GeoEye-1	$D_{\lambda }$	✗	0.200	0.126	0.111	0.163	0.004	0.035	0.177	0.037
		✓	0.228	0.102	0.154	0.174	0.007	0.025	0.175	0.036
	$D_s$	✗	0.056	0.125	0.172	0.155	0.055	0.116	0.139	0.051
		✓	0.031	0.080	0.034	0.035	0.051	0.081	0.037	0.038
	QNR	✗	0.755	0.765	0.736	0.708	0.941	0.853	0.709	0.914
		✓	0.748	0.827	0.817	0.797	0.942	0.896	0.794	0.927

reports the results for the four datasets obtained by eight benchmark methods with and without postprocessing by our proposed EBP at full scale. Besides,

Table 3. Quantitative comparison for the eight methods without (denoted by ✗) and with (denoted by ✓) postprocessing by our proposed EBP on four kinds of satellite datasets in the reduced scale. Note that, for better comparison, the results improved by EBP are highlighted in blue, otherwise the original results are colored in green.

Datasets	Metrics	EBP	BDSD	GFPCA	MF	GSA	NLIHS	PRACS	SFIM	PNN
IKONOS	CC	✗	0.923	0.916	0.937	0.933	0.922	0.946	0.942	0.931
		✓	0.940	0.948	0.941	0.939	0.930	0.943	0.942	0.937
	ERGAS	✗	3.234	3.620	3.051	3.016	3.372	2.887	3.019	2.906
		✓	2.878	2.592	2.730	2.773	2.733	2.704	2.703	3.038
	SAM	✗	3.989	4.429	3.599	3.773	4.323	3.523	3.663	3.572
		✓	3.155	3.082	3.017	3.125	3.244	3.111	3.071	3.429
	RMSE	✗	21.43	23.93	20.19	19.99	22.45	19.05	19.89	19.21
		✓	19.16	17.32	18.22	18.48	18.65	18.13	18.04	19.99
	Q4	✗	0.854	0.768	0.860	0.832	0.815	0.873	0.862	0.843
		✓	0.876	0.884	0.878	0.858	0.854	0.879	0.878	0.860
QuickBird	CC	✗	0.846	0.917	0.929	0.938	0.926	0.945	0.931	-
		✓	0.928	0.954	0.949	0.945	0.949	0.940	0.947	-
	ERGAS	✗	3.945	3.923	2.912	2.816	3.517	2.626	3.080	-
		✓	3.214	2.309	2.509	2.615	2.385	2.773	2.471	-
	SAM	✗	4.391	3.444	2.514	2.942	2.981	2.796	2.649	-
		✓	2.269	2.089	2.090	2.187	2.183	2.269	2.144	-
	RMSE	✗	55.55	53.29	39.64	38.47	47.48	36.20	41.70	-
		✓	44.14	31.71	34.43	35.90	33.06	38.78	33.94	-
	Q4	✗	0.838	0.755	0.904	0.891	0.834	0.912	0.882	-
		✓	0.902	0.933	0.932	0.924	0.930	0.917	0.929	-
WorldView-2	CC	✗	0.868	0.947	0.964	0.972	0.961	0.975	0.965	0.974
		✓	0.954	0.985	0.984	0.985	0.984	0.985	0.984	0.983
	ERGAS	✗	11.481	7.534	4.863	4.939	6.383	4.625	5.297	4.332
		✓	6.601	3.138	3.348	3.057	2.990	3.201	3.213	3.375
	SAM	✗	12.179	5.547	3.476	3.879	4.692	3.881	3.736	4.542
		✓	5.152	2.702	2.899	2.763	3.039	3.161	2.909	3.163
	RMSE	✗	121.64	79.62	51.80	52.12	66.92	45.97	56.01	46.36
		✓	70.03	33.16	35.37	32.37	32.40	35.16	33.99	35.43
	Q8	✗	0.860	0.847	0.953	0.947	0.908	0.958	0.937	0.965
		✓	0.939	0.979	0.979	0.981	0.979	0.980	0.979	0.978
GeoEye-1	CC	✗	0.852	0.795	0.894	0.863	0.875	0.890	0.892	0.908
		✓	0.903	0.922	0.912	0.914	0.917	0.921	0.910	0.918
	ERGAS	✗	3.724	4.309	2.982	3.278	3.586	3.177	3.157	2.783
		✓	3.002	2.539	2.706	2.676	2.592	2.515	2.707	2.623
	SAM	✗	5.276	5.152	3.920	4.661	4.018	3.848	3.943	3.591
		✓	3.549	3.389	3.411	3.517	3.444	3.421	3.429	3.422
	RMSE	✗	34.82	39.86	27.74	30.59	32.81	29.86	29.21	25.95
		✓	27.97	23.73	25.33	25.06	24.27	23.63	25.30	24.41
	Q4	✗	0.833	0.573	0.823	0.796	0.722	0.781	0.790	0.842
		✓	0.875	0.882	0.875	0.884	0.871	0.883	0.868	0.880

presents the corresponding results by the aforementioned five numeric metrics at a reduced scale based on Wald’s protocol [49]. Note that the best result for each method with and without our proposed EBP as postprocessing is shown in boldface blue. According to the two tables, one can see that the eight benchmark methods produce fusion products that resulted in better quantitative scores for those five metrics at reduced scale and the three metrics at full scale when our proposed EBP algorithm was used for postprocessing, compared by those obtained without our proposed EBP as postprocessing, except that the original results of GFPCA and PRACS have a relatively better spectral quality (corresponding to the metrics $D_{\lambda }$) for the IKONOS datasets and the NLIHS for the QuickBird dataset in the full scale experiments. As for the GeoEye-1 dataset in the full scale experiments, although the results of GSA, MF and NLIHS methods after postprocessing by our proposed EBP have a little bit precision loss in spectral quality, they have much better performances on the whole as shown by the value of the index QNR in Table 2. Another interesting observation, however, can be made for BDSD method on GeoEye-1 dataset in the full scale experiments. In this case, we find that the BDSD algorithm did not seem to benefit from the postprocessing, resulting in higher $D_{\lambda }$ while lower QNR scores. This effect may be due to the severer spectral distortion of the fused products of BDSD as shown in Figure 4a. This reveals that additional attention must be taken when using the proposed method to process significantly spectral distortion images.

As for the reduced scale experiments, it is shown in Table 3 that the results of eight benchmark methods with our proposed EBP as postprocessing can almost always obtain the improved performances with the five metrics except for the PRACS method for the IKONOS and QuickBird datasets and the PNN method for the IKONOS dataset.

Overall, the results in Table 2 and Table 3 suggest that the use of our proposed EBP as postprocessing can be beneficial since a clear pattern of improvement based on visual inspection (as shown in Figure 1, Figure 2, Figure 3 and Figure 4) and quantitative comparison on five metrics at the reduced scale and three metrics at the full scale for four kinds of satellite datasets were observed for eight benchmark algorithms.

5. Conclusions

Pansharpening is an important preprocessing step for a variety of applications based on high resolution multispectral images. Over the past decades, lots of techniques, e.g., HM, HPM, the MTF-matched filter, and so on, have been proposed in this literature. However, the postprocessing for pansharpening has not received sufficient attention. Specifically, the integration of the HM, HPM and MTF-matched filter for pansharpening postprocessing is rarely investigated. In this paper, a post-processing method, called enhanced back-projection (EBP), which strives to integrate the HM, HPM and MTF-matched filters into the back-projection algorithm, is present and applied to pansharpening. Eight most benchmark pansharpening methods were selected and several experiments on four different kinds of satellite datasets were carried out to verify the effectiveness of the proposed EBP. The experimental results have shown that although our proposed EBP as postprocessing does not always guarantee better results, but it does frequently produce results that are equal to or better than the results without using our proposed EBP as postprocessing. Additionally, it should be noted that the proposed EBP as postprocessing works as an independent but complementary postprocessing module for pansharpening and thus does not need any modification of the existing pansharpening methods, which is a highly desirable feature. As is well-known, radiometric calibration (except the spectral and spatial) information is relevant to classification accuracy of remote sensing data. Thus, the studies on whether the radiometric calibration is maintained or lost after applying the post-processing pansharpening algorithm may be the focus of future studies.