This paper presents a nonlinear mixing model for hyperspectral image unmixing. The proposed model assumes that the pixel reflectances are post-nonlinear functions of unknown pure spectral components (referred to as endmembers) contaminated by an additive white Gaussian noise. The nonlinear effects affecting endmembers are approximated by polynomial functions leading to a polynomial post-nonlinear mixing model. A Bayesian strategy is used to estimate the parameters of this model yielding an unsupervised nonlinear unmixing algorithm. Due to the large number of parameters to be estimated, an efficient constrained HamiltonianMarkov chain Monte Carlo method is developed to sample according to the posterior of the Bayesian model. The performance of the resulting unmixing strategy is evaluated on synthetic data.