There are a huge number of harmonics in the railway power supply system. Accurately estimating the harmonic impedance of the system is the key to evaluating the harmonic emission level of the power supply system. A harmonic impedance estimation method is proposed in this paper, which takes the Gaussian mixture regression (GMR) as the main idea, and is dedicated to calculating the harmonic impedance when the load changes or the background harmonic changes in the traction power supply system. First, the harmonic voltages and currents are measured at the point of common coupling (PCC); secondly, a Gaussian mixture model (GMM) is established and optimized parameters are obtained through the EM algorithm; finally, a Gaussian mixture regression is performed to obtain the utility side harmonic impedance. In the simulation study, different harmonic impedance estimation models with uniform distribution and Gaussian distribution are established, respectively, and the harmonic impedance changes caused by different system structures in the railway power supply system are simulated. At the same time, the error is compared with the existing method to judge the accuracy and robustness of this method. In the case analysis, the average value, average error, standard deviation and other indicators are used to evaluate this method. Among them, the average error and standard deviation of this method are about one-fifth to one-third of those of the binary linear regression (BLR) method and the independent random vector (IRV) method. At the same time, its index is slightly better than that of the support vector machine (SVM) method.