Hyperspectral unmixing (HU) is an important task in hyperspectral image (HSI) processing, which estimates endmembers and their corresponding abundances. Generally, the unmixing process of HSI can be approximated by a linear mixing model. Since each type of material only appears in a few pixels in HSI, sparse non-negative matrix factorization (NMF) and its extensions are regarded as the main unmixing methods. However, due to the lack of consideration of the spatial distribution in the local domain, most HU methods lose the inducing effect of spatial correlation in the decomposition process, and ignore the correlation of the real distribution of different materials in different images. In this paper, a novel NMF unmixing model is proposed, called SMRNMF, which learns multiple subspace structures from the original hyperspectral images and combines them into a sparse NMF framework to improve the performance of the model. Firstly, subspace clustering is embedded into the sparse NMF model. According to the spatial correlation of the original data, two similarity matrices are learned, which make full use of the local correlations between the pixels of the original data. Secondly, based on the self-expression characteristics of the data in the subspace, the global similarity pixel graph matrix is embedded into the model to construct a self-expression regularizer to improve the unmixing performance. Finally, a smoothing matrix is cleverly constructed and embedded in the model to overcome the adverse effects of noisy information in the abundance images. Experiments on several simulated and real HSI data sets show that our method has superior performance compared with the existing methods.