Hyperspectral unmixing is one of the most important procedures for remote sensing image processing. The non-negative matrix factorization (NMF)-based method has been widely used for hyperspectral unmixing since it can get endmember and abundance matrices simultaneously. However, the inherent single-decomposition structure of NMF may not achieve good performance for highly mixed data. To solve this issue, we propose a spectral and spatial total-variation-regularized multilayer non-negative matrix factorization (SSTV-MLNMF) for hyperspectral unmixing. In SSTV-MLNMF, we designed an effective multilayer factorization process and combined spectral and spatial total variation as extra regularization. These could enhance the smoothness for spectral signatures and spatial fields, which could achieve better performance. Experiments on both synthetic and real datasets have validated the effectiveness of our method and have shown that it has outperformed several state-of-the-art approaches of hyperspectral unmixing.