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.