Close to saturation operation of high power amplifier (HPA) leads to strong nonlinear and dispersive characteristic of satellite channels. At the receiver, the observation signals are distorted by not only the linear inter-symbol interference (ISI), but also the nonlinear ones, which makes it challenge to perform optimal detection. In this paper, we study factor graph (FG)-based turbo equalization for nonlinear satellite channels characterized by Volterra series. Factor nodes on FG are classified into Belief propagation (BP) set and variational message passing (VMP) set to enable low complexity combined message passing implementation while with high performance. BP is used on the hard constraint nodes, such as demapping and decoding, while VMP is employed to update messages of the likelihood function node. It is shown that, without any approximation on the Volterra series channel model, messages can be expressed in a closed form via canonical parameters, and the extrinsic information from equalizer to decoder is derived in an explicit way. Simulation results demonstrate the superior performance of the proposed combined VMP-BP algorithm with low computational complexity.