In order to enhance the spatial resolution of multispectral images and avoid severe spectral distortion, the amplified multispectral images and the original panchromatic image were decomposed by non-subsampled contourlet transform with the same layers. A 5×5 sliding window around the fusion pixel of the matrix in the low frequency sub-band of the decomposed images was selected to compute the local correlation coefficient and fourth-order correlation coefficient. If the calculated local correlation coefficient was greater than the value of fourth-order correlation coefficient, the center pixel of the matrix of the high frequency coefficients of multispectral images was replaced by high frequency coefficients of panchromatic image because the position of the ground objects were similar on the spectral characteristics. Otherwise, the high frequency coefficients of multispectral image were kept unchanged. At last, the fused image was reconstructed from the joint of approximation coefficients of multispectral image and detail coefficients obtained from the fusion procedure with inverse non-subsampled contourlet transform. The experimental simulation was conducted on Landsat 7 remote sensing images and the subjective and objective evaluation index of the fused results were given by comparing with the existing similar latest literature. The experimental results indicate that the proposed approach has good effects both on enhancing the original multispectral image and maintaining good spectral characteristics.