Hyperspectral unmixing is one of the most important procedures for remote sensing image processing. The multilayer non-negative matrix factorization (MLNMF)-based method has been widely used for hyperspectral unmixing due to its good performance for highly mixed data with multiple-decomposition structure. However, few works consider the spatial information in the image, which may enhance the performance. In order to solve this issue, we propose a homogeneous region regularized multilayer non-negative matrix factorization (HR-MLNMF) method for hyperspectral unmixing. In HR-MLNMF, the spatial information, depicted by the homogeneous region, is applied to regularize MLNMF, which could enhance the smoothness of each homogeneous spatial field to achieve better performance. Experiments on both synthetic and real datasets have validated the effectiveness of our method and shown that it has outperformed several state-of-the-art approaches of hyperspectral unmixing.