Due to the changeable internal mechanisms or external working conditions, the degradation trend of products usually presents multiple phase characteristics. However, most existing multi-phase degradation methods treat each phase as linear and rely on artificial designation to directly specify the forms of drift models, which may result in an inaccurate description of degradation characteristics for complex equipment. To this end, we formulate a general nonlinear multi-phase degradation model with three-source variability based on the Wiener process. Meanwhile, a stage division method is developed to automatically determine the degradation phase number, change-point locations, and the forms of drift models. Then, we obtain the expressions of remaining useful life (RUL) by considering the uncertainty of change-point degradation observations. Especially, we derive the approximate analytical solution of RUL based on the linear model. Furthermore, to fully consider the unit-to-unit heterogeneity and utilize the degradation observations of the in-service unit and prior information simultaneously, we propose a parameter estimation method based on variational Bayesian approach, which adaptively updates all parameters as random variables. Finally, two numerical examples and three practical examples are provided to verify the effectiveness of the proposed method.