In this paper, a new pansharpening method, which uses nonnegative matrix factorization, is proposed to enhance the spatial resolution of remote sensing multispectral images. This method, based on the linear spectral unmixing concept and called joint spatial-spectral variables nonnegative matrix factorization, optimizes, by new iterative and multiplicative update rules, a joint-variables criterion that exploits spatial and spectral degradation models between the considered images. This criterion considers only two unknown high spatial-spectral resolutions variables. The proposed method is tested on synthetic and real datasets and its effectiveness, in spatial and spectral domains, is evaluated with established performance criteria. Results show the good performances of the proposed approach in comparison with other standard literature ones.