Hyperspectral unmixing is an important technique for estimating fraction of different land covers from remote sensing imagery. In recent years, nonnegative matrix factorization (NMF) methods with various constraints have been introduced into hyperspectral unmixing. Among these methods, graph based constraint has been proved to be useful in capturing the latent manifold structure of the hyper-spectral data in the feature domain. However, due to the complexity of the data, only using single graph can not adequately reflect the intrinsic property of the data. In this paper, we propose a multiple graph regularized NMF method for hyperspectral unmixing, which approximates the manifold and consistency of data by a linear combination of several graphs constructed in different scales. Results on both synthetic and real data have validated the effectiveness of the proposed method, and shown that it has outperformed several state-of-the-arts hyperspectral unmixing methods.