In order to preserve both spectral and spatial information simultaneously in fused image, a method of remote sensing image fusion based on improved gradient projection Non-negative Matrix Factorization (NMF) and Complex Contourlet Transform (CCT) was proposed. Firstly, the panchromatic image was histogram matched to the intensity component of multispectral image to obtain a new panchromatic image. Then, the intensity component of multispectral image and the new panchromatic image were decomposed by CCT, respectively. Corresponding low frequency images and high frequency images were obtained. Next, the new low frequency image was obtained by fusion of two low frequency images using improved gradient projection NMF, and the new high frequency image was obtained by fusion of two high frequency images with absolute maximum criterion. Finally, the fused image was reconstructed by inverse CCT and inverse Hue-Saturation-Intensity (HSI) transform. A large number of experimental results show that, compared with other existing relative methods such as HSI method, the method combining NMF with nonsubsampled contourlet transform and the method combining lifting wavelet with HSI, the images fused by proposed method have higher spatial resolution and more spectral information.