We introduce a robust mixing model to describe hyperspectral data resulting from the mixture of several pure spectral signatures. The new model extends the commonly used linear mixing model by introducing an additional term accounting for possible nonlinear effects, that are treated as sparsely distributed additive outliers. With the standard nonnegativity and sum-to-one constraints inherent to spectral unmixing, our model leads to a new form of robust nonnegative matrix factorization with a group-sparse outlier term. The factorization is posed as an optimization problem, which is addressed with a block-coordinate descent algorithm involving majorization-minimization updates. Simulation results obtained on synthetic and real data show that the proposed strategy competes with the state-of-the-art linear and nonlinear unmixing methods.