The increasing traction station load of high-speed railway with single-phase, high-power, spatiotemporal mobility, and uncertainty has brought major impacts and challenges to the power system. The non-negligible characteristics are the spatial distribution of traction stations on high-speed railway lines, the time sequence of traction load moving with EMUs, as well as the nonparametric probabilistic features and correlation of traction load. A method of time sequence three-phase probabilistic power flow (3PPF) calculation for power system including traction station load of high-speed railway is proposed in this paper. Firstly, through the port transformation of traction power supply system (TPSS), the three-phase equivalent model of TPSS is established to construct an unbalanced three-phase power flow calculation model of grid with traction station load. Secondly, considering the spatial correlation, time sequence and non-parametric probabilistic feature of traction load, the time sequence probabilistic models of multiple traction station loads are built based on diffusion-based kernel density estimator (DKDE). Then, the algorithm of time sequence 3PPF calculation based on three-point estimation method (3PEM) is applied to deal with the computational efficiency and accuracy problems of 3PPF considering correlation and timing. Finally, the two modified three-phase power systems are employed for case studies in combination with the measured data of traction station of a high-speed railway. The validity and computing efficiency of the proposed method are verified.