Spectral unmixing of hyperspectral images is an important issue in the fields of remotesensing. Jointly exploring the spectral and spatial information embedded in the data is helpful toenhance the consistency between mixing/unmixing models and real scenarios. This paper proposesa graph regularized nonlinear unmixing method based on the recent multilinear mixing model(MLM). The MLM takes account of all orders of interactions between endmembers, and indicates thepixel-wise nonlinearity with a single probability parameter. By incorporating the Laplacian graphregularizers, the proposed method exploits the underlying manifold structure of the pixels’ spectra,in order to augment the estimations of both abundances and nonlinear probability parameters.Besides the spectrum-based regularizations, the sparsity of abundances is also incorporated for theproposed model. The resulting optimization problem is addressed by using the alternating directionmethod of multipliers (ADMM), yielding the so-called graph regularized MLM (G-MLM) algorithm.To implement the proposed method on large hypersepectral images in real world, we proposeto utilize a superpixel construction approach before unmixing, and then apply G-MLM on eachsuperpixel. The proposed methods achieve superior unmixing performances to state-of-the-artstrategies in terms of both abundances and probability parameters, on both synthetic and real datasets.