This paper proposes a novel driving cycle construction method in consideration of velocity, road slope, and passenger load, based on a real-world bus route with a plug-in hybrid electric bus (PHEB). The main purpose is to address the disadvantage that an inaccurate reflection of the real-world driving characteristics for city buses will be caused when ignoring the passenger load in the course of a driving cycle synthesis. Two contributions are supplemented to distinguish from the previous research. Firstly, a novel station-based method is proposed aiming at developing a driving cycle with high accuracy. The kinematic segments are partitioned according to the distance of adjacent bus stops, while a two-dimensional Markov chain Monte Carlo method is employed to synthesize driving cycle between each interval of adjacent bus stops. Secondly, the random passenger load for different bus stops is treated as a discrete Markov chain model, according to the correlation analysis of the measured passenger data which are distinguished for off-peak and peak hours. Meanwhile, Monte Carlo simulation and maximum likelihood estimation are utilized to determine the most likely number of passengers for each bus stop. At last, the fuel consumption of the PHEB is simulated with the best-synthesized driving cycle and contrasted to the mean fuel consumption of the later measured data which is composed of the velocity, road slope, and the passenger load. The results demonstrate that the synthesized driving cycle has a higher accuracy on fuel consumption estimation.