Multi-sensor image fusion for pansharpening in remote sensing

Abstract

The main objective of this article is quality assessment of pansharpening fusion methods. Pansharpening is a fusion technique to combine a panchromatic image of high spatial resolution with multispectral image data of lower spatial resolution to obtain a high-resolution multispectral image. During this process, the significant spectral characteristics of the multispectral data should be preserved. For images acquired at the same time by the same sensor, most algorithms for pansharpening provide very good results, i.e. they retain the high spatial resolution of the panchromatic image and the spectral information from the multispectral image (single-sensor, single-date fusion). For multi-date, multi-sensor fusion, however, these techniques can still create spatially enhanced data sets, but usually at the expense of the spectral consistency. In this study, eight different methods are compared for image fusion to show their ability to fuse multitemporal and multi-sensor image data. A series of eight multitemporal multispectral remote sensing images is fused with a panchromatic Ikonos image and a TerraSAR-X radar image as a panchromatic substitute. The fused images are visually and quantitatively analysed for spectral characteristics preservation and spatial improvement. It can not only be proven that the Ehlers fusion is superior to all other tested algorithms, it is also the only method that guarantees excellent colour preservation for all dates and sensors used in this study.

Keywords:

• image fusion
• multitemporal
• multi-sensor
• pansharpening
• radar
• quality assessment

1. Introduction

Most of the operating earth observation satellites such as Landsat, Spot, Ikonos, Quickbird, Formosat, GeoEye or Orbview provide panchromatic images at a higher spatial resolution than in their multispectral mode. The difference in spatial resolution between the panchromatic and the multispectral mode can be measured by the resolution ratio (RR), i.e. the ratio of their respective ground sampling distances (GSD) that may vary between 1 : 2 and 1 : 5 (Ehlers et al. 2009). A large number of fusion or pansharpening techniques have been developed to combine panchromatic and multispectral images to produce an enhanced multispectral image of high spatial resolution. Most of these methods seem to work well with images that were acquired at the same time by one sensor (single-sensor, single-date fusion) (Price 1987, Pohl and van Genderen 1998, Alparone et al. 2007, Ehlers 2008).

However, a number of sensors acquire information about the earth in either only panchromatic (EROS-A, EROS-B, Worldview-1) or only multispectral (RapidEye) mode. It becomes therefore increasingly important to fuse image data from different sensors which are usually recorded at different dates. Thus, there is a need to investigate techniques that allow multi-sensor, multi-date image fusion.

With multi-sensor fusion, an additional problem has to be considered. If data from different satellites are used for fusion, the RR can get worse. For example, the RR between RapidEye (6.5 m GSD) and Worldview (0.5 m GSD) is 1 : 13, for Quickbird (0.61 m GSD) and SPOT-5 multispectral (10 m GSD) about 1 : 16. With the advent of high-resolution radar data from the TerraSAR-X or Radarsat 2 sensor, studies are underway to investigate the potential of an active sensor for pansharpening (Klonus 2008).

Generally, image fusion methods can be differentiated into three levels: pixel level (iconic), feature level (symbolic) and knowledge or decision level. Of highest relevance for remote sensing are techniques for iconic image fusion, as the data are least altered and most of the pansharpening techniques have been developed for this fusion (Wald et al. 1997, Pohl and van Genderen 1998, Zhang 2004a). Iconic image fusion techniques can be grouped into three classes: colour-related techniques, statistical methods and numerical methods (Pohl and van Genderen 1998). The first class includes colour compositions of three image bands in the RGB colour space as well as the more sophisticated colour transformations such as the intensity-hue-saturation (IHS) or hue-saturation-value (HSV) transform. The second class includes band statistics, such as the principal component (PC) or Gram–Schmidt (GS) transforms. The third class includes arithmetic operations such as image multiplication, summation and image rationing as well as sophisticated numerical approaches such as wavelets (Chipman et al. 1995), curvelets (Candès 1999) or ridgelets (Chen et al. 2005).

Wald (2002) proposed another categorisation of image fusion techniques: projection and substitution methods, relative spectral contribution and the ‘spatial improvement by injection of structures’ (améloration de la resolution spatial par injection de structures = ARSIS) concept. In the first class are those methods, which project the image into another coordinate system and substitute one component. The second class is comparable with the second class of Pohl and van Genderen (1998), with the exception that this category is restricted to band rationing and arithmetic combinations. The basis of the ARSIS concept is a multi-scale technique to inject the high spatial information into the multispectral images. Although this classification scheme bears some merits, we selected the one introduced by Pohl and van Genderen (1998) because of its mathematical precision.

Research on image fusion techniques in remote sensing started as early as the late 1980s and early 1990s (Cliche et al. 1985, Welch and Ehlers 1987, Chavez et al. 1991, Ehlers 1991) and concentrated on pixel level fusion (pansharpening). With the advent of the very high-resolution satellite programmes, the interest in pansharpening has significantly increased. Many publications focused on how to retain the spectral characteristics of the multispectral data after pansharpening (see, e.g. Alparone et al. 2007 or Thomas et al. 2008).

Currently, still uncomplex image fusion techniques such as IHS or Brovey are widely used for pansharpening. Often these methods do not retain the spectral characteristics of the multispectral datasets. More sophisticated techniques seem to work well for many applications when used for single-sensor, single-date fusion. When applied to multi-sensor data, they can still create spatially enhanced datasets but usually at the expense of the spectral consistency (Ehlers 2004, Klonus and Ehlers 2007, Klonus 2008).

Over the last few years, a number of improved algorithms have been developed with the promise to minimise colour distortion while retaining the spatial improvement of the standard data fusion algorithms. Especially, wavelet techniques led to a number of new fusion methods (Otazu et al. 2005, Lillo-Saavedra and Gonzalo 2006, Yunhao et al. 2006). We selected six advanced fusion methods for comparison. These algorithms are colour normalisation spectral sharpening (CN), GS fusion, modified IHS fusion, Ehlers fusion, University of New Brunswick (UNB) fusion and a wavelet-based fusion, the proportional additive wavelet fusion (AWLP). For comparison, we also selected the ‘traditional’ Brovey and PC fusion methods.

2. Study area and datasets

The study area is located in the North of Spain representing the region around Santo Domingo de la Calzada (Figure 1). This area was used as a control site of the Joint Research Centre of the European Commission (JRC) in the project ‘Control with Remote Sensing of Area-Based Subsidies’ (CwRS) (http://agrifish.jrc.it/marspac/DCM/). A series of eight multitemporal, multispectral remote sensing images (seven SPOT scenes and one Formosat scene) was used (Table 1). These multispectral images cover a time frame of almost two years and virtually all seasons and thus pose an excellent challenge for a spectral characteristics preserving data fusion. All images were fused with a panchromatic Ikonos image from 30 May 2005 (1 m GSD) and in a second investigation with a TerraSAR-X spotlight image with 1 m spatial resolution as a panchromatic substitute. The TerraSAR-X radar image was recorded on 3 May 2008 and was provided by the German Aerospace Centre DLR. The image was despeckled with a proprietary algorithm developed at the University of Würzburg in Germany.

Figure 1. Study site near Santo Domingo de la Calzada, Spain.

Table 1. Multispectral remote sensing datasets for the study site.

Satellite sensor	Recording date	GSD (m)	Abbreviation
SPOT 4	3 November 2003	20	SPOT 4 2003
SPOT 4	24 April 2004	20	SPOT 4 2004a
SPOT 4	15 May 2004	20	SPOT 4 2004b
SPOT 5	24 July 2004	10 (20 m SWIR)	SPOT 5 2004
SPOT 2	10 December 2004	20	SPOT 2 2004
SPOT 5	10 April 2005	10 (20 m SWIR)	SPOT 5 2005
SPOT 4	20 July 2005	20	SPOT 4 2005
Formosat 2	12 August 2005	8	Formosat

This agricultural area was chosen, because most of the colour changes appear in agricultural fields depending on crop phenology and harvesting dates. These spectral changes pose a serious problem for most fusion methods and consequently prevent a correct visual or automated agricultural assessment.

In urban areas, the colour changes are less severe. Here, the dominant challenge is the spatial improvement due to man-made structures (e.g. houses, roads). It is therefore mandatory to investigate not only the spectral fidelity of the multi-date, multi-sensor fusion algorithms but also the actual spatial improvement. Consequently, our test sites also included houses and roads as artificial structures.

3. Methodology

3.1 Preprocessing

The creation of coregistered datasets is the first and an extremely important part of the fusion process. Especially for radar data with its geometric properties such as foreshortening, layover, and shadowing, it is necessary to ensure that the images are in perfect registration. Otherwise, the quality of the fused image will decrease significantly. Even small misalignments can cause colour artefacts in the fusion process. Fortunately, the new German radar satellite TerraSAR-X provides image data in a perfect quality and a spatial resolution of 1 m in its spotlight mode (Schubert et al. 2008). This provided an excellent basis for a correct rectification. The other images were orthocorrected and accuracy checked by different European Union (EU) data contractors as part of the standard CwRS routine (JRC 2008).

After the orthocorrection of the original multispectral images, they were resampled using a cubic convolution function to the ground resolution of the high-resolution images (1 m). All bands of the original images were used for the fusion process.

3.2 Image fusion algorithms

We used eight different algorithms for fusion of the images. All eight images were fused with the panchromatic Ikonos image resulting in a dataset of 64 fused images. In a second investigation, four of these images (Formosat, SPOT 5 2005 and SPOT 4 2004a, b) were fused with the TerraSAR-X image. All images were visually and statistically evaluated for colour preservation and spatial improvement. Use was made of standard techniques such as Brovey and PC fusion and a number of recently developed algorithms: CN spectral sharpening fusion, GS fusion, modified IHS fusion, Ehlers fusion as well as two methods that participated in the IGARSS image fusion contest, the UNB fusion and the AWLP. The AWLP method ranked first in the fusion contest (Alparone et al. 2007). We are aware that other methods have been proposed but most are either in an experimental stage or have not been fully published. We concentrated on those that were readily available in remote sensing image analysis packages and/or were described in detail so that they could be implemented. The algorithms are described in the following sections.

3.2.1 Brovey transform

The Brovey transform is a combination of arithmetic operations. Each spectral band is first divided by the sum of the three chosen bands and then multiplied with the panchromatic image (Hallada and Cox 1983). For our study, this method was extended to four bands because the multispectral datasets were all four-band images.

3.2.2 Principal component transform

PC transform is a statistical technique that transforms a multivariate dataset of correlated variables into a dataset of uncorrelated linear combinations of the original variables. For images, it creates an uncorrelated feature space that can be used for further analysis instead of the original multispectral feature space. The PC is applied to the multispectral bands. The panchromatic image is histogram matched to the first PC. It then replaces the selected component and an inverse PC transform takes the fused dataset back into the original multispectral feature space (Chavez et al. 1991). The advantage of the PC fusion is that the number of bands is not restricted (such as the original IHS or Brovey fusion). It is, however, a statistical procedure which means that it is sensitive to the area to be sharpened and produces fusion results that may vary depending on the selected image subsets (Ehlers and Klonus 2004).

3.2.3 Modified IHS transform

The IHS transform separates spatial (intensity) and spectral (hue and saturation) information from a standard RGB image. The intensity refers to the total brightness of the image, hue to the dominant or average wavelength of the light contributing to the colour and saturation to the purity of colour. This is comparable with the human perception. To fuse the images, three bands of a multispectral image are transformed from the RGB domain into the IHS colour space. The panchromatic component is matched to the intensity of the IHS image and then replaces the intensity component. We make use of the modified IHS fusion from Siddiqui (2003) which was developed for a better fit of the fused multispectral bands to the original data. After the matching, the panchromatic image replaces the intensity in the original IHS image and the fused image is transformed back into the RGB colour space (Ehlers 1991). For more than three bands, the IHS transform has to be used more than once depending on the number of bands to be fused.

3.2.4 Additive wavelet proportional fusion

The AWL method (Nunez et al. 1999) is one of the existing multi-resolution wavelet-based image fusion techniques. It was originally designed for a three-band RGB multispectral image. In this method, the spectral signature is preserved because the high-resolution panchromatic structure is integrated into the luminance L-band of the original low-resolution multispectral image. Therefore, this method is only defined for three bands. It was extended to an arbitrary number of bands by Otazu et al. (2005). It maintains the spectral signature of an n-band image in the same way as AWL does with RGB images. This generalised method is called AWLP. The algorithm was implemented in Matlab for this study by the authors.

3.2.5 Colour normalisation spectral sharpening

CN spectral sharpening is an extension of the Brovey transform. The input image bands are grouped into spectral segments defined by the spectral range of the panchromatic image. The corresponding band segments are processed in the following manner: each input band is multiplied by the sharpening band and then normalised by dividing it by the sum of the input bands in the segment. Bands outside the spectral range of the panchromatic image are not sharpened (Vrabel et al. 2002).

3.2.6 Gram Schmidt fusion

The GS process transforms a set of vectors into a new set of orthogonal and linear independent vectors. By averaging the multispectral bands, the GS fusion simulates a low-resolution panchromatic band. As the next step, a GS transform is performed for the simulated panchromatic band and the multispectral bands with the simulated panchromatic band applied as the first band. Then the high spatial resolution panchromatic band replaces the first GS component. Finally, an inverse GS transform is applied to create the pansharpened multispectral bands (Laben et al. 2000).

3.2.7 Ehlers fusion

The Ehlers fusion was developed specifically for a spectral characteristics preserving image merging (Ehlers 2004). It is based on an IHS transform coupled with a Fourier domain filtering. The principal idea behind a spectral characteristics preserving image fusion is that the high-resolution image has to sharpen the multispectral image without adding new grey level information to its spectral components. An ideal fusion algorithm would enhance high-frequency changes such as edges and grey level discontinuities in an image without altering the multispectral components in homogeneous regions. To facilitate these demands, two prerequisites have to be addressed. First, colour and spatial information have to be separated. Second, the spatial information content has to be manipulated in a way that allows an adaptive enhancement of the images. This is achieved by a combination of colour and Fourier transforms.

For optimal colour separation, use is made of an IHS transform. This technique is extended to include more than three bands by using multiple IHS transforms until the number of bands is exhausted. If the assumption of spectral characteristics preservation holds true, there is no dependency on the selection or order of bands for the IHS transform. Subsequent Fourier transforms of the intensity component and the panchromatic image allow an adaptive filter design in the frequency domain. Using fast Fourier transform (FFT) techniques, the spatial components to be enhanced or suppressed can be directly accessed. The intensity spectrum is filtered with a low-pass filter (LP) whereas the spectrum of the high-resolution image is filtered with an inverse high-pass filter (HP). After filtering, the images are transformed back into the spatial domain with an inverse FFT and added together to form a fused intensity component with the low-frequency information from the low-resolution multispectral image and the high-frequency information from the high-resolution panchromatic image. This new intensity component and the original hue and saturation components of the multispectral image form a new IHS image. As the last step, an inverse IHS transformation produces a fused RGB image that contains the spatial resolution of the panchromatic image and the spectral characteristics of the multispectral image. These steps can be repeated with successive three band selections until all bands are fused with the panchromatic image. The order of bands and the inclusion of spectral bands for more than one IHS transform are not critical because of the colour preservation of the procedure (Klonus and Ehlers 2007).

3.2.8 University of New Brunswick fusion

The UNB fusion algorithm (Zhang 2004b) applies first histogram standardisation for the multispectral and panchromatic bands of the input images. The multispectral bands in the spectral range of the panchromatic image are selected and a regression analysis is calculated using a least square algorithm. The results are used as weight factors for the multispectral bands. Via multiplication with the corresponding bands and a following addition, a new synthesised image is produced. To create the fused image, each standardised multispectral image is multiplied with the standardised panchromatic image and divided by the synthesised image. This method was primarily designed for single-sensor, single-date images and is used as the standard method for Quickbird pansharpening. For this study, the implementation in PCI Geomatica was used.

3.3 Evaluation methods

3.3.1 Visual analysis

For visual analysis, colour preservation and spatial improvement were examined. For colour preservation, each band was evaluated for grey level changes and different band combinations were checked for spectral consistency. Although a number of different band combinations were tested, in this article we concentrate on the standard false colour infrared display because this band combination is widely used for many remote-sensing applications. For spatial improvement analysis, we checked for each band, if more details could be recognised.

These visual analysis techniques, however, have the disadvantage that they are subjective and depend heavily on the experience of the respective interpreter. Therefore, we also applied a number of quantitative-statistical tests.

3.3.2 Spectral evaluation

A number of statistical evaluation methods are used to measure after-fusion colour fidelity. These methods should be objective, reproducible, and of quantitative nature. Therefore, we selected the following statistical measures.

The correlation coefficient (CC) between the original multispectral bands and the equivalent fused bands is the most frequently used method to evaluate the spectral value preservation (see, e.g. González-Audícana et al. 2005 or Klonus and Ehlers 2007). The values range from −1 to 1. The best correspondence between fused and original image data shows the highest correlation value and should be close to 1.

For a per-pixel deviation (PD; Wald 2002), the fused image is first degraded to the spatial resolution of the original image and is then subtracted from the original image on a per-pixel basis. In a final step, the average deviation per pixel is calculated using the number of pixel in the scene. Low differences between the original and the fused image are indicated by low values of the PD. The lowest (and best) value is 0. A perfect 0, however, means that original and fused image are identical and no pansharpening is performed. A good correspondence is therefore indicated by values close to 0.

As an additional measurement, the root mean square error (RMSE) was also proposed by Wald (2002). It is computed by the difference of the standard deviation and the mean of the fused and the original image. The best possible value is again 0.

The structure similarity (SSIM) index was proposed by Wang et al. (2004). The SSIM is a further development of the Q index, which is widely used to evaluate the spectral value preservation of pansharped images (Wang and Bovik 2002, Alparone et al. 2004). This method combines a comparison of luminance, contrast and structure and is applied locally in an 8 × 8 window. This window is moved pixel-by-pixel over the entire image. At each pixel, the local statistics and the SSIM index are calculated within the window. The values vary between 0 and 1. Values close to 1 show the highest correspondence with the original images.

3.3.3 Spatial improvement evaluation

The objective is to find the fused image with the optimal combination of spectral characteristics preservation and spatial improvement. If the image shows no spatial improvement, no pansharpening is performed which makes the fusion obsolete. The ‘fused image’ would score best for spectral evaluation with no spatial improvement, however. Most evaluations have concentrated on spectral consistency with little or no emphasis on spatial improvement. To avoid this shortcoming, we selected two different quantitative methods to measure the quality of spatial improvement.

For the calculation of the high-pass correlation (HCC) coefficient, a HP filter with a 3 × 3 Laplacian kernel is first applied to the panchromatic image and to each band of the fused image. Then the correlation coefficients between the HP filtered bands and the HP filtered panchromatic image are calculated. This analysis was proposed by Zhou et al. (1998). According to Pradhan et al. (2006), there is no reason for the preference of a specific filter, any HP filter can be used for this purpose.

For the second measure, an edge detection (ED) approach was used. A Canny edge detector (Canny 1986) is applied to the panchromatic image and each band of the fused multispectral image. The detected edges are then compared to the panchromatic image edges for each individual band. ED correspondence is measured in per cent; 100% means that all the edges in the panchromatic image are detected in the fused image. Only fusion techniques which produce spectral fidelity and spatial improvement should be used, especially for multi-date and multi-sensor fusion.

4. Fusion results and quality evaluation

For the investigation of the electro-optical (E–O) fusion (i.e. pansharpening with the Ikonos image) a multi-date, multi-sensor dataset with 64 fused images was produced from eight original images with eight different fusion methods. For the investigation using TerraSAR-X radar data as panchromatic substitute (radar fusion), 32 fused images were produced from four original images and again the eight different fusion methods. All images were visually and statistically analysed for all fusion methods. Our comprehensive evaluation is based on an individual analysis of all 96 fusion results.

To illustrate our findings in the following section, we use a subset of the SPOT 5 2004 image, because it is very representative for the encountered fusion effects. Figure 2 shows the original multispectral SPOT 5 2004 (band 3, near infrared; band 2, red; and band 1, green) as a standard false colour infrared display and the panchromatic Ikonos image. For the investigation of the radar fusion, a subset of the SPOT 5 2005 image is used to illustrate the differences between the fusion methods. It was not possible to use the same subsets for both fusion tests, because Ikonos and TerraSAR-X had no overlapping ground coverage (Figure 3).

Figure 2. Original multispectral SPOT 5 2004 image (a) recorded on 24 July 2004 (© CNES 2004) in the band combination 3 (near infrared), 2 (red), 1 (green) resampled to 1m and the original panchromatic Ikonos image and (b) recorded on 30 May 2005 (satellite image by GeoEye). Note the inverse brightness in the images for the fields on the left side.

Figure 3. Original multispectral SPOT 5 2005 (© CNES 2005) image (a) in the band combination 3 (near infrared), 2 (red), 1 (green) resampled to 1m and the original panchromatic TerraSAR-X image and (b) recorded on 3 May 2008(©DLR 2008). Note the differences in the field structure for the agricultural areas in the upper right part.

The results for the panchromatic fusion are shown in Figures 4–10 (even numbers); the ones for the radar fusion are presented in Figures 5–11 (odd numbers). All multispectral images are again shown in the band combination near infrared, red, green.

Figure 4. Fused multispectral SPOT 5 2004 image using Brovey (a) and PC (b) in the band combination 3 (near infrared), 2 (red), and 1 (green).

Figure 5. Fused multispectral SPOT 5 2005 image using Brovey (a) and PC (b) in the band combination 3 (near infrared), 2 (red), and 1 (green).

Figure 6. Fused multispectral SPOT 5 2004 image using IHS (a) and AWLP (b) in the band combination 3 (near infrared), 2 (red), and 1 (green). Note the artefacts in homogeneous areas for the ALWP image.

Figure 7. Fused multispectral SPOT 5 2005 image using IHS (a) and AWLP (b) in the band combination 3 (near infrared), 2 (red), and 1 (green).

Figure 8. Fused multispectral SPOT 5 2004 image using CN (a) and GS (b) in the band combination 3 (near infrared), 2 (red), and 1 (green).

Figure 9. Fused multispectral SPOT 5 2005 image using CN (a) and GS (b) in the band combination 3 (near infrared), 2 (red), and 1 (green).

Figure 10. Fused multispectral SPOT 5 2004 image using Ehlers (a) and UNB (b) in the band combination 3 (near infrared), 2 (red), and 1 (green).

4.1 Visual analysis

In the fused Brovey images, the improved structure is clearly visible. For the E–O fusion (Figure 4a) colour distortions are visible, especially for the fields in the west. The reason for this is the grey value distribution of the panchromatic Ikonos image. The brightness structure of the panchromatic image is directly injected into the multispectral image and causes the alteration in colour. The best result for this fusion technique was achieved for the SPOT 5 2005 image which was recorded in April, i.e. the closest to the acquisition date of the panchromatic Ikonos image. In the fused radar image (Figure 5a), the colour changes are striking; in some parts of the image, colours appear that seem ‘unnatural’ in a false colour infrared display.

The fused PC images show some spatial improvement but not as detailed as for the Brovey fusion. In both images (Figures 4b and 5b), dramatic changes in the spectral values appear. The reason for these colour changes is again directly related to the panchromatic/radar information. For the E–O fusion, the fused image falsely suggests that the bare field in the southwest is actually vegetated which renders it inappropriate for any photointerpretation or automated classification process. For the radar fusion, the near infrared band with high grey values in the original image is represented by low grey values in the fused image. These low values come from the radar image.

In general, the spatial structures in the IHS fused images are improved, however, not as much as for the Brovey fusion. The IHS fusion shows sometimes acceptable results but never without serious colour changes (Figure 6a). The best result for the E–O fusion image was achieved for the SPOT 5 2005 image. This image was recorded in April, the shortest time difference to the acquisition date of the panchromatic Ikonos image. In the fused radar images (Figure 7a), the IHS methods show a better colour preservation than most of the other methods. However, there are still significant colour changes.

The worst results for the spatial enhancement are associated with the AWLP fusion (Figures 6b and 7b). The spatial structures are only slightly improved. Especially along the edges, artefacts appear which disrupt and blur the information in the whole image. Particularly in homogeneous areas, no spatial improvement could be found at all. This is true for all images of the E–O fusion process. The spectral values are relatively well preserved making it the second best method for colour fidelity.

For the E–O fusion, the CN spectral sharpening shows good to very good spatial improvement in all images, but often associated with serious spectral distortions. The results were often similar to those of the Brovey fusion (Figure 8a). This comes as no surprise as the CN fusion is an extension of the Brovey method. For the radar fusion, the TerraSAR-X header had to be manipulated so that the wavelength of the multitemporal images could be adapted to the radar image. With this, it was possible to use the TerraSAR-X image as a substitute for a high-resolution panchromatic image. This incompatibility was probably the main reason for the enormous change of the spectral values in the fused radar images (Figure 9a). This is especially true for regions where the TerraSAR-X image exhibits extremely high digital numbers (DNs), such as in edges of buildings. These colour artefacts distort also the spatial resolution.

The GS fusion shows some improvement in spatial resolution but never as much as Brovey, PC, CN or IHS. In some cases, the fused images look still blurred (Figure 8b). The correspondence of the spectral values to the original is low. For the E–O fusion, it displays similar colour deficiencies as the PC fusion. This holds also true for the radar fusion (Figure 9b). In the near infrared band, high grey values in the original images are replaced with low values in the fused images. This leads to spectral changes, especially in this band.

It is clearly visible that only the Ehlers fusion retains the spectral characteristics of the original image (Figures 10a and 11a). This is true for all images, as there are almost no visible colour changes in comparison to the original. The spatial improvement is evident and provides good to very good results, particularly for the E–O fusion. For the radar fusion, the general colour preservation is satisfactory. Only for the edges of buildings, which act as corner reflectors, high DNs are produced which do not exist in the multispectral data. As the Ehlers fusion injects these edges into the multispectral images, this effect cannot be avoided. On the other hand, the buildings become more distinct, a fact that can support a subsequent photointerpretation.

Figure 11. Fused multispectral SPOT 5 2005 image using Ehlers (a) and UNB (b) in the band combination 3 (near infrared), 2 (red), and 1 (green).

The worst results are produced by the UNB fusion (Figures 10b and 11b). The spectral information is replaced by a more or less homogeneous structure for all bands which makes the fused image look like a black-and-white display. In addition, the edges are deformed and many artefacts are visible in the image. Therefore, the spatial structure is also not much improved. It has to be noted that this dramatic effect appears only in two cases for the E–O fusion, the SPOT 5 2005 and the SPOT 2 2004. The reason for this is probably that this method uses reference images for the regression, which are in the spectral range of the panchromatic image. However, all the other fused images (including the radar fusion) also show severe colour distortions (Figure 11b). This function was developed for single-sensor, single-date fusion. The reason, why it does not work for multi-sensor fusion, however, is not clear. It may have something to do with its implementation in the used software package. Consequently, it should not be used for multitemporal and/or multi-sensor fusion.

4.2 Quantitative spectral analysis

Although a visual analysis has to be a significant part of all quality evaluation procedures, it has the disadvantage that it is very subjective and depends on the interpreter (Klonus and Ehlers 2009). Statistical analysis techniques have to accompany the visual interpretation for objective and reproducible quantitative results. For this purpose, we make use of the quantitative-statistical methods presented in Section 3.3. Comprehensive results of this evaluation are presented in Tables 2–13. Best values in each column (scene) in the tables are marked with bold and the worst in grey. In general, the quantitative analyses for colour preservation did confirm the visual inspection findings. The Ehlers fusion scored the highest marks for colour preservation in nearly 90% of the images.

Table 2. CCs for the E–O fusion.

	SPOT 4 2003	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2004	SPOT 2 2004	SPOT 5 2005	SPOT 4 2005	Formosat
Brovey	0.715	0.730	0.666	0.625	0.418	0.945	0.564	0.472
PC	0.721	0.418	0.440	0.686	0.443	0.747	0.397	0.294
IHS	0.682	0.843	0.773	0.792	0.314	0.942	0.662	0.511
AWLP	0.926	0.933	0.939	0.986	0.723	0.997	0.936	0.999
CN	0.777	0.830	0.795	0.735	0.399	0.965	0.725	0.471
GS	0.735	0.387	0.372	0.709	0.415	0.853	0.341	0.417
Ehlers	0.996	0.998	0.998	0.985	0.971	0.999	0.999	0.994
UNB	0.228	0.434	0.395	0.554	0.063	0.173	0.363	0.236

Table 3. CCs for the radar fusion.

	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2005	Formosat
Brovey	0.480	0.504	0.448	0.055
PC	0.508	0.602	0.382	0.129
IHS	0.844	0.837	0.833	0.285
AWLP	0.904	0.952	0.927	0.942
CN	0.525	0.517	0.555	0.105
GS	0.616	0.618	0.612	0.459
Ehlers	0.984	0.980	0.990	0.987
UNB	0.323	0.418	0.244	0.144

Table 4. Deviation per pixel for the E–O fusion.

	SPOT 4 2003	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2004	SPOT 2 2004	SPOT 5 2005	SPOT 4 2005	Formosat
Brovey	57.8	64.2	95.3	65.1	8.22	73.5	129	84.5
PC	8.00	10.2	16.7	2.43	8.11	10.4	44.4	17.4
IHS	5.74	5.39	9.13	3.37	5.70	3.29	10.9	6.72
AWLP	2.29	3.24	4.43	0.62	2.80	0.65	5.33	0.19
CN	0.92	8.10	23.7	2.56	64.8	7.55	38.4	21.0
GS	4.97	10.3	14.6	3.28	4.27	5.95	16.1	6.86
Ehlers	0.84	0.74	0.71	0.80	0.60	0.38	0.72	0.70
UNB	16.7	10.1	14.0	4.43	5.05	18.8	14.6	7.59

Table 5. Deviation per pixel for the radar fusion.

	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2005	Formosat
Brovey	89.95	97.46	67.64	101.82
PC	31.61	42.14	64.31	72.24
IHS	33.29	7.83	5.76	7.30
AWLP	32.79	4.19	4.01	1.82
CN	30.28	27.62	12.62	19.25
GS	34.65	11.75	9.26	5.91
Ehlers	1.54	2.32	1.28	0.57
UNB	38.56	14.06	12.54	6.68

Table 6. RMSE for the for the E–O fusion.

	SPOT 4 2003	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2004	SPOT 2 2004	SPOT 5 2005	SPOT 4 2005	Formosat
Brovey	59.68	68.32	101.25	65.78	8.30	77.82	134.40	85.54
PC	7.63	12.61	20.95	4.45	7.96	14.71	45.76	16.11
IHS	3.78	3.26	8.31	3.56	3.45	2.36	6.30	2.13
AWLP	1.12	1.33	0.88	0.46	1.44	2.33	1.21	0.66
CN	9.90	10.50	26.33	4.92	67.16	8.45	40.65	20.70
GS	0.72	1.34	1.84	0.24	0.30	1.88	3.67	1.64
Ehlers	0.38	0.42	0.45	0.24	0.16	0.32	0.32	0.46
UNB	32.61	0.93	2.41	0.01	0.10	25.24	8.68	2.93

Table 7. RMSE for the radar fusion.

	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2005	Formosat
Brovey	93.09	101.00	55.53	103.12
PC	29.96	55.30	81.84	105.48
IHS	34.02	6.86	1.84	3.51
AWLP	33.57	11.13	9.85	18.41
CN	33.33	27.25	14.14	26.14
GS	34.84	6.94	2.12	0.52
Ehlers	0.20	0.27	0.24	0.23
UNB	45.21	13.42	6.87	6.46

Table 8. SSIM for the E–O fusion.

	SPOT 4 2003	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2004	SPOT 2 2004	SPOT 5 2005	SPOT 4 2005	Formosat
Brovey	0.326	0.267	0.162	0.398	0.837	0.248	0.099	0.281
PC	0.737	0.491	0.448	0.862	0.830	0.631	0.438	0.716
IHS	0.666	0.779	0.673	0.855	0.833	0.899	0.509	0.682
AWLP	0.845	0.806	0.780	0.957	0.830	0.991	0.763	0.997
CN	0.762	0.770	0.698	0.812	0.304	0.920	0.566	0.693
GS	0.723	0.477	0.419	0.844	0.878	0.757	0.350	0.714
Ehlers	0.991	0.994	0.995	0.969	0.987	0.999	0.996	0.991
UNB	0.355	0.446	0.373	0.773	0.838	0.277	0.249	0.638

Table 9. SSIM for the radar fusion.

	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2005	Formosat
Brovey	0.172	0.120	0.221	0.024
PC	0.294	0.395	0.121	0.131
IHS	0.629	0.667	0.706	0.490
AWLP	0.714	0.820	0.744	0.871
CN	0.319	0.328	0.376	0.201
GS	0.406	0.441	0.584	0.635
Ehlers	0.961	0.949	0.967	0.987
UNB	0.182	0.221	0.312	0.554

Table 10. HCC coefficients for the E–O fusion.

	SPOT 4 2003	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2004	SPOT 2 2004	SPOT 5 2005	SPOT 4 2005	Formosat
Brovey	0.940	0.937	0.939	0.952	0.934	0.938	0.942	0.948
PC	0.547	0.619	0.680	0.594	0.948	0.592	0.799	0.729
IHS	0.933	0.938	0.882	0.869	0.864	0.893	0.859	0.839
AWLP	−0.059	−0.069	−0.060	0.019	−0.076	−0.100	−0.049	−0.016
CN	0.765	0.763	0.760	0.791	0.871	0.771	0.697	0.970
GS	0.985	0.584	0.589	0.873	0.940	0.525	0.941	0.976
UNB	0.972	0.958	0.938	0.983	0.625	0.979	0.053	0.960
Ehlers	0.966	0.955	0.936	0.933	0.894	0.957	0.925	0.963

Table 11. HCC coefficients for the radar fusion.

	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2005	Formosat
Brovey	0.856	0.983	0.859	0.644
PC	0.367	0.621	0.481	0.480
IHS	0.459	0.562	0.394	0.298
AWLP	0.001	−0.025	−0.068	−0.112
CN	0.101	0.111	0.105	0.125
GS	0.509	0.353	0.472	0.856
Ehlers	0.582	0.570	0.527	0.568
UNB	−0.726	0.298	0.231	−0.756

Table 12. ED results for the E–O fusion.

	SPOT 4 2003	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2004	SPOT 2 2004	SPOT 5 2005	SPOT 4 2005	Formosat
Brovey	96.1	95.7	96.2	96.7	95.3	95.2	96.1	95.8
PC	95.9	96.1	95.7	94.9	95.6	94.6	95.3	95.2
IHS	94.9	94.9	94.0	93.8	93.9	92.6	93.5	93.1
AWLP	86.7	86.3	86.3	86.3	87.5	85.7	86.0	86.6
CN	95.2	94.8	95.1	95.0	95.6	93.8	95.3	96.3
GS	96.4	96.6	96.5	95.0	95.4	96.6	96.5	95.5
Ehlers	95.0	94.6	94.2	94.8	93.6	92.0	94.1	92.7
UNB	97.1	96.2	96.2	96.8	95.3	95.5	92.3	96.4

Table 13. ED results for the radar fusion.

	SPOT 4 2004a	SPOT 4 2004b	SPOT 5 2005	Formosat
Brovey	95.5	97.1	95.7	93.9
PC	94.8	95.4	93.8	93.9
IHS	89.3	90.0	88.3	87.0
AWLP	84.0	83.4	81.8	78.0
CN	90.8	90.6	90.4	90.2
GS	95.3	94.8	93.9	93.2
Ehlers	88.7	88.6	90.1	89.8
UNB	93.7	92.6	89.9	94.9

4.2.1 Correlation coefficient

The calculation of the CC has the benefit that a (constant) bias between the fused and the original image does not alter the results. Best results were obtained by the Ehlers fusion for 10 images and by the AWLP method for two images. Brovey, IHS, PC, and GS fusion have their best results for the SPOT 5 2005 image, which is the closest in time to the acquisition date of the panchromatic Ikonos image. In general, however, CN, IHS, Brovey, PC, and GS exhibit low CCs (with the exception of SPOT 5 2005). The lowest CCs were calculated for the UNB fusion. Highest value was only 0.554 for the SPOT 5 2004 scene. Only the Ehlers fusion exceeded 0.9 for all images. Table 2 presents the correlation results for the E–O fusion. The coefficient values are calculated as average value for all bands.

CCs are lower for most methods for the radar fusion (Table 3). One exception is the IHS function where the values are higher on average than for the E–O case. Again, the Ehlers fusion has the highest correlation with coefficients above 0.9 for all images. Second is the AWLP fusion. All other methods produced low correlation values between 0.05 and 0.6. CCs are the lowest for the UNB and PC fusion, thus confirming the visual analysis.

4.2.2 Per-pixel deviation

The calculated PD values are presented in Tables 4 and 5. The lowest (and best) values are associated with the Ehlers fusion again confirming the results of the visual analysis. Results for AWLP and the modified IHS are in most cases acceptable. For the SPOT 4 2004b image, however, also these two methods fail to produce satisfactory results. The GS fusion produces acceptable results for only three images with deviation values below 5.0. PC, CN, UNB and Brovey fusion produce insufficient results, especially for the radar fusion with a PD above 5.0 for all methods and exceeding 10.0 for most of the images (Table 5). The worst results are associated the Brovey fusion with PD values above 50 for all images with the exception of SPOT 2 2004. It should be mentioned that the result of the Brovey method was adapted to the 8-bit range, but no histogram matching was applied. Nevertheless, the visual inspection and the low CCs prove that the spectral quality of the Brovey method is not sufficient for image fusion.

4.2.3 Root mean square error

The best values for the RMSE (Table 6) are again associated with the Ehlers fusion. Not in correspondence with the visual interpretation are the results for the GS fusion. For the E–O test case, the GS fusion shows the second best results whereas in the visual analysis it was ranked very low. However, all RSME values are still a factor of 4–20 higher than those of the Ehlers fusion. A close third is the AWLP fusion. Unacceptable results were produced by the PC, CN and Brovey methods. Especially, Brovey shows again the highest (worst) values, with the exception of the winter scene (SPOT 2 2004). The UNB fusion shows very different results ranging from 0.01 to 32.6, where 0.01 is the best value for one band of the SPOT 5 2004 scene and 32.6 the worst value for one band of the SPOT 4 2003 scene.

For the radar fusion (Table 7), all methods with the exception of the Ehlers fusion produce very poor results for the scene SPOT 4 2004a. For the three other scenes, the two simplest methods, Brovey and PC, show the worst results. AWLP, IHS, and GS results are somewhat better but still a factor of 2–10 worse than the Ehlers fusion. Especially, the values for the GS fusion could not be confirmed by the visual analysis. CN and UNB did not produce acceptable results.

4.2.4 Structure similarity index

Again, best results for the E–O and the radar fusion were achieved by the Ehlers algorithm followed by the AWLP (Tables 8 and 9). The Ehlers fusion was the only one with SSIM values of 0.9 and better for all images.

For the radar case (Table 9), no other fusion algorithm scored above 0.9. For the E–O fusion, PC, IHS, CN and GS show some acceptable results with values between 0.4 and 0.9. For the radar fusion, the values are generally lower (with the exception of Ehlers). UNB and Brovey produce unsatisfactory result for both test cases.

4.3 Quantitative spatial analysis

In most fusion papers, the quality assessment is only concerned with spectral quality. The idea of pansharpening, however, is to inject the spatial structure of the high-resolution panchromatic image into the multispectral bands. If only the spectral fidelity is considered, the unaltered multispectral image would reach the best results which would make the fusion obsolete. Consequently, we assessed also the spatial quality of the fusion algorithms. Here, the results are very different from the spectral analysis. From the visual analysis, we could already conclude that the spatial enhancement in the fused images proved to be sufficient for most of the methods.

4.3.1 High-pass correlation

The CCs after HP filtering are to some degree influenced by the different fusion types. The reason for this should be investigated, because visual inspection showed sometimes different results. For the E–O fusion, most of the methods could improve the structure in the image (Table 10). Brovey and Ehlers were very consistent and produced average values above 0.9. IHS and CN were lower but also consistent, whereas PC, CN, and UNB were very inconsistent. High values for SPOT 2 2004 (PC and GS) or SPOT 4 2004a (UNB) are contradicted by low values for other dates such as SPOT 5 2005 (PC and GS) or SPOT 4 2005 (UNB). These results make their use for image fusion unreliable. The lowest values are produced by the AWLP with CC around zero meaning that there is no correspondence between the fused bands and the high-resolution image. We suspect that the observed artefacts are the reason for the insufficient results of the HP correlation analysis (Table 10).

In general, the radar fusion showed lower CCs than the E–O fusion which can be explained by the special characteristics of radar imagery (Table 11). Best results were obtained by the Brovey transform with the only value above 0.9 (SPOT 4 2004b). CCs were above 0.6 for all Brovey fused images. The only other function that was consistently above 0.5 was the Ehlers fusion. All other methods were inconsistent (PC, IHS and GS) or even insufficient (AWLP, CN and UNB).

4.3.2 Edge detection

For the ED evaluation, results were very similar for most of the fusion functions (Tables 12 and 13). With the exception of the AWLP method, all of them show very good results with an edge correspondence above 90%. The only fusion function with results below 90% for all images was the AWLP. This confirms the visual interpretation results that the ALWP does not really improve the spatial resolution. Best results were obtained by UNB and CN, which failed completely in the spectral quality assessment.

For radar fusion, the results differ somewhat from the E–O case. Here, Brovey is the best overall with values between 94 and 97%. The other functions above 90% are PC, CN, and GS followed by UNB, Ehlers and IHS. The latter functions had values between 87 and 95%, which means that they still produce sufficient results. AWLP is the only fusion technique which has values below 90% for all images with a minimum of 78%, making it the only technique that does not produce adequate spatial improvement.

5. Conclusion

In this study, we could prove that standard and most of the advanced fusion methods cannot cope with the demands that are placed on them by multi-sensor/multi-date fusion. The spectral distortions are manifold: brightness reversions, a complete change of spectral characteristics, artificial artefacts or unnatural and artificial colours are the results of many operational fusion techniques that are readily available in many commercial image processing systems. Fusion methods such as PC, CN, GS or UNB should only be used for single-sensor, single-date images. Wavelet-based fusions can retain most of the spectral characteristics which comes unfortunately at the expense of spatial improvement. The AWLP wavelet method produced additional spatial artefacts instead of spatial improvements. This is probably caused by the wavelet characteristics.

For the radar fusion, we can conclude that simple replacement of the panchromatic image by a high-resolution radar image has to be treated with caution. It could be shown that pansharpening techniques could be used to improve the spatial resolution of multispectral images with TerraSAR-X data as high-resolution input. It is also evident that only the Ehlers fusion was capable of preserving the spectral characteristics of the multispectral image making it suitable for further image analysis (photointerpretation or classification). The only technique of those that were applied in our comparative study that delivers pansharpened images with almost no spectral change is the Ehlers fusion.

As discussed in Section 3.3, a combined method for a quantitative assessment of spatial improvement and spectral preservation is needed because, otherwise, the best colour preservation is observed if no pansharpening is performed. This would, however, defeat the purpose of the fusion process. A combined assessment of spatial improvement and spectral characteristics preservation can be used as a comprehensive quality measure for image fusion. So far, we have used a narrative combination of results for different quantitative analyses. A mathematical equation that integrates spectral- and spatial-quality evaluation would most likely be a better and more consistent quality indicator. The development of such an integrative quality indicator should therefore be high on the list for future fusion research.

Acknowledgements

Thanks to the JRC for providing the data for this research. The TerraSAR-X data and project is supported by the German BMWI through the DLR with the contract number 50EE0704.