^{1}

^{*}

^{2}

^{2}

^{2}

Due to the data acquired by most optical earth observation satellite such as IKONOS, QuickBird-2 and GF-1 consist of a panchromatic image with high spatial resolution and multiple multispectral images with low spatial resolution. Many image fusion techniques have been developed to produce high resolution multispectral image. Considering panchromatic image and multispectral images contain the same spatial information with different accuracy, using the least square theory could estimate optimal spatial information. Compared with previous spatial details injection mode, this mode is more accurate and robust. In this paper, an image fusion method using Bidimensional Empirical Mode Decomposition (BEMD) and the least square theory is proposed to merge multispectral images and panchromatic image. After multi-spectral images were transformed from RGB space into IHS space, next I component and Panchromatic are decomposed by BEMD, then using the least squares theory to evaluate optimal spatial information and inject spatial information, finally completing fusion through inverse BEMD and inverse intensity-hue-saturation transform. Two data sets are used to evaluate the proposed fusion method, GF-1 images and QuickBird-2 images. The fusion images were evaluated visually and statistically. The evaluation results show the method proposed in this paper achieves the best performance compared with the conventional method.

In the remote sensing satellite system, data storage amount and data transmission speed are limited, thus the spectral resolution of the image data usually is inversely proportional to the spatial resolution; the spectral resolution of the multispectral (MS) images is low; and the spatial resolution is high on the contrary; the spatial resolution of panchromatic (Pan) image is high; while the spectral resolution is low, for instance Landsat, IKONOS, and QuickBird-2. Many applications such as map updating, land-use classification, natural disaster and environment monitoring require remote sensing images with both high spectral resolution and spatial resolution. High spatial resolution information can allow accurate geometric analysis, while high spectral resolution information is essential for the objective interpretation of objects [

Many remote sensing image fusion methods have been developed in the last few years [

Empirical Mode Decomposition (EMD) [

Image fusion technique based on IHS transform is a representative method, which has been widely used to merge MS images and pan image. Due to its fast implementation and high efficiency, it has been integrated into some remote sensing image processing software, such as ENVI, ERDAS. It can convert a color image from RGB space to IHS space, the transformation models include sphere transformation, cylinder transformation, triangle transformation and hexagonal cone transformation and so on [

( I V 1 V 2 ) = ( 1 3 1 3 1 3 − 2 6 − 2 6 − 2 6 1 2 − 1 2 0 ) ( R G B ) (1)

H = tan − 1 ( V 2 V 1 ) S = V 1 2 + V 2 2 (2)

Inverse transformation formula is,

( R ′ G ′ B ′ ) = ( 1 − 1 2 1 2 1 − 1 2 − 1 2 1 2 0 ) ( I V 1 V 2 ) (3)

EMD is a new signal analysis method [

Bidimensional Empirical Mode Decomposition (BEMD) is an expansion of one dimensional EMD in two dimensional space [

1) Initializing. Letting r 0 , 0 = f ( x , y ) , r l , k represents the kth iteration process sifting for the lth residue.

2) Calculate IMFs. h l , k denotes the lth IMF sifted after k iteration process. Details of calculation IMFs are as follows:

a) Letting h l , k = r l , k .

b) Finding all local maxima and minima of h l , k , implements surface interpolation using maxima and minima, and obtain a maxima envelope surface u max ( x , y ) and minima envelope surface u min ( x , y ) .

c) Calculating average envelope surface.

m e a n l , k ( x , y ) = ( u max ( x , y ) + u min ( x , y ) ) / 2 (4)

d) Calculating h l , k − 1 .

h l , k − 1 ( x , y ) = h l , k − m e a n l , k ( x , y ) (5)

e) Calculating standard deviation (SD) of stopping iteration condition. If SD < ε ( ε is a predefined threshold (often 0.2 - 0.3)), then we can consider h l , k is the lth BIMF, otherwise, repeat steps (b) to (e) until the stop condition is satisfied.

S D = ∑ x = 0 M − 1 ∑ y = 0 N − 1 | h l , k ( x , y ) − h l , k − 1 ( x , y ) | h l , k 2 ( x , y ) (6)

3) Calculating Residue.

r i ( x , y ) = r i − 1 ( x , y ) − B I M F i ( x , y ) (7)

If r i ( x , y ) includes more than K extreme point, then go to step (1) to continue decompose until the number of r i ( x , y ) ’s extreme point is less than K.

4) Finally, the f ( x , y ) can be presented as follows:

f ( x , y ) = ∑ i = 1 n B I M F i ( x , y ) + r n ( x , y ) (8)

Previous studies [

1) Assume that the vector Y does not vary with time, multiple observations are made to obtain its true values,

Y = C X + V (9)

where Y is a vector which size is m × 1 , C is an observation matrix which size is m × n , X is an observation value and V is an observation error which it’s mean is zero.

2) In the case that the number of observation equation is more than the number of unknown variables, the optimal estimate of the true value X can be expressed as follows:

{ X ¯ L S W = ( C T R − 1 C ) − 1 C T R − 1 Y R = V a r [ V ] (10)

where V a r [ V ] is the observation error matrix.

3) Assume three MS images and a Pan image which spatial resolution ratio are n : 1 are used to produce fusion images, the observation error of panchromatic spatial information is r , and the corresponding multispectral image observation error is n r . The parameter of Formula (9) is as follows,

Y = C X + V

Y = [ Y 1 Y 2 Y 3 Y 4 ] , C = [ 1 1 1 1 ] T

R = [ r 2 0 0 0 0 ( n r ) 2 0 0 0 0 ( n r ) 2 0 0 0 0 ( n r ) 2 ]

Then the optimal evaluation of true X is,

X ¯ L S W = ( C T R − 1 C ) − 1 C T R − 1 Y = n 2 n 2 + 3 Y 1 + 1 n 2 + 3 Y 2 + 1 n 2 + 3 Y 3 + 1 n 2 + 3 Y 4 (11)

The detailed fusion steps are as follows, and

1) Geometric registration of the images to be fused.

2) Transform MS images into Intensity, hue and saturation component by IHS transform.

3) Apply histogram matching between intensity component and pan image, and produce a new pan image.

4) Decompose intensity component and new pan image using BEMD method, and obtain IMFs and reduce component corresponding.

5) Take IMFs into Formula (11) to calculate optimal IMFs which is optimal spatial information evaluation.

6) Use Formula (8) to reconstruct new I component.

7) Perform Inverse IHS transform to complete fusion.

In this paper, we choose two sets of data that comes from different sensor and different region to verify the effectiveness of the fusion method, and the fusion result of the proposed method is compared with the fusion result of DWT- and BEMD-based fusion method [

1 16 | 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 | (12)

[

As show in

Due to the spatial resolution of quickbird-2 image in

In terms of objective evaluation, we selects the correlation coefficient (CC) and the distortion degree (DD) to evaluate the spectral properties of the fusion image relative to the reference image, selects the high frequency correlation coefficient (HFCC) [

Correlation coefficient (CC) reflects the similarity of spectral properties of two images. If the correlation coefficients tend to be 1, the spectral properties of the two images are close to each other, otherwise the opposite. The Correlation coefficients can be calculated as follows:

CC = ∑ i = 1 M ∑ j = 1 N [ X ( x i , y i ) − X ¯ ] [ Y ( x i , y i ) − Y ¯ ] ∑ i = 1 M ∑ j = 1 N [ X ( x i , y i ) − X ¯ ] ∑ i = 1 M ∑ j = 1 N [ Y ( x i , y i ) − Y ¯ ] (13)

where X ¯ and Y ¯ are means of image X and image Y.

The distortion degree (DD) indicates the deviation of the fused image relative to the reference image. The smaller the distortion degree, the closer the spectral information between two images is. The distortion degree can be calculated as follows:

DD = 1 M N ∑ i = 1 M ∑ j = 1 N | X ( x i , y j ) − Y ( x i , y j ) | (14)

The high frequency correlation coefficient (HFCC) presents the similarity of the high frequency information of the fusion image and that of the reference image. The main idea is based on the fact the most spatial information in image is concentrated in the high frequency information. If the high frequency correlation coefficients tend to be 1, the spatial structure of the images is close to each other. The filter which we used to extract high frequency of image is Laplace mask as follows:

[ − 1 − 1 − 1 − 1 8 − 1 − 1 − 1 − 1 ] (15)

The structural similarity Index Measurement (SSIM) is a comprehensive index which takes into account the gray correlation, the gray deviation and the image contrast of the image. The closer the value is to 1, the closer the two images are. The SSIM can be calculated as follows:

SSIM = 4 σ X Y μ X μ Y ( σ X 2 + σ Y 2 ) [ ( μ X ) 2 + ( μ Y ) 2 ] (16)

where σ X Y is the covariance between image X and image Y, μ X and μ Y are the means, σ X and σ Y are variance of image X and image Y.

From

Fusion method | Band | CC | DD | HFCC | SSIM |
---|---|---|---|---|---|

DWT + IHS | R | 0.693 | 0.076 | 0.198 | 0.464 |

G | 0.718 | 0.068 | 0.219 | 0.068 | |

B | 0.636 | 0.077 | 0.193 | 0.402 | |

mean | 0.682 | 0.073 | 0.203 | 0.311 | |

DWT + HIS + LS | R | 0.761 | 0.066 | 0.208 | 0.501 |

G | 0.718 | 0.068 | 0.219 | 0.508 | |

B | 0.706 | 0.068 | 0.201 | 0.436 | |

mean | 0.728 | 0.067 | 0.209 | 0.481 | |

BEMD | R | 0.874 | 0.049 | 0.207 | 0.535 |

G | 0.841 | 0.051 | 0.215 | 0.537 | |

B | 0.814 | 0.056 | 0.197 | 0.453 | |

mean | 0.843 | 0.052 | 0.206 | 0.508 | |

The method proposed in this paper | R | 0.899 | 0.043 | 0.258 | 0.566 |

G | 0.866 | 0.047 | 0.259 | 0.555 | |

B | 0.864 | 0.045 | 0.244 | 0.508 | |

mean | 0.876 | 0.045 | 0.253 | 0.543 |

Fusion method | Band | CC | DD | HFCC | SSIM |
---|---|---|---|---|---|

DWT+IHS | R | 0.839 | 0.091 | 0.390 | 0.387 |

G | 0.825 | 0.098 | 0.408 | 0.414 | |

B | 0.819 | 0.102 | 0.391 | 0.335 | |

mean | 0.828 | 0.097 | 0.396 | 0.379 | |

DWT+IHS+LS | R | 0.873 | 0.082 | 0.412 | 0.412 |

G | 0.864 | 0.087 | 0.430 | 0.443 | |

B | 0.861 | 0.090 | 0.415 | 0.366 | |

mean | 0.866 | 0.086 | 0.419 | 0.407 | |

BEMD+IHS | R | 0.918 | 0.062 | 0.267 | 0.409 |

G | 0.927 | 0.058 | 0.280 | 0.456 | |

B | 0.929 | 0.059 | 0.271 | 0.394 | |

mean | 0.925 | 0.060 | 0.273 | 0.420 | |

The method proposed in this paper | R | 0.935 | 0.054 | 0.396 | 0.455 |

G | 0.942 | 0.051 | 0.413 | 0.506 | |

B | 0.944 | 0.051 | 0.398 | 0.440 | |

mean | 0.940 | 0.052 | 0.402 | 0.467 |

component’s high frequency information, their fusion images’ correlation coefficient is more larger, and the distortion is more smaller. It shows the least square evaluation is a more efficient model than the model which the high frequency information is replaced directly in spectral preservation.

The correlation coefficients of high frequency of the fusion images based on the method proposed in this paper are slightly higher than the corresponding index of the image using other three methods, the index of image based on other three methods are basically close, just the method based on wavelet transform, it is not efficient for improving high frequency correlation coefficient to inject spatial details using least square method. The ranking of SSIM indexes is the same as the ranking of correlation coefficients and distortions. The order of the two kinds of indexes shows that the fusion image based on the method proposed in this paper is closer to the reference image in spectral information and spatial information.

The index’s distribution in

In order to make fusion image contain the same spectral information and spatial information as the image obtained using the same parameter sensor, this paper proposes a fusion method based on BEMD + IHS to merge MS images and Pan image. Simultaneously, considering spatial information in MSI and Pan forms a set of observation values of unequal precision, an optimal evaluation of spatial information can be calculated via the least square method. Two sets of data, GF-1 images and QuickBird-2 images, are used to verify the fusion method. After visual analysis and objective evaluation, experimental result shows that the fusion results based on BEMD proposed in this paper are much closer to the reference image than fusion results based on wavelet transform, and the fused images using least square method to integrate spatial detail are much closer to the reference image than fusion image using panchromatic image’s high frequency information to substitute I component’s high frequency information. In addition, the fused image using substitution model to inject spatial details exists the phenomenon that spatial information is injected overly.

This research is supported by Natural Science Foundation of Hunan Province (14JJ7039) and Educational Commission of Hunan Province (14C1071), China.

Huang, D.S., Yang, P., Li, J. and Ma, C.H. (2017) Remote Sensing Image Fusion Using Bidimensional Empirical Mode Decomposition and the Least Squares Theory. Journal of Computer and Communications, 5, 35-48. https://doi.org/10.4236/jcc.2017.512004