The profile is a crucial characteristic of a ballast particle. It is important to characterize the uncertainties of three-dimensional ballast particle profiles and to determine the effects of these profiles on the dynamic performance of ballast bed. In this study, the spherical harmonics method was employed to achieve efficient reconstruction of the three-dimensional profiles of 300 ballast particles, which were obtained by performing indoor scanning. The relationships between elongated or flaky particles and the various orders of spherical harmonics functions were proven. The findings demonstrate that, when the spherical harmonics order exceeded 15, the errors in the volume, surface area and sphericity of the reconstructed particles were less than 5% compared to real particles. Furthermore, the spherical harmonics spectra approximately conformed to a gamma distribution, for which the average coefficient of determination (({R}^{2})) was 0.95. The second order relative spherical harmonics frequency significantly influenced the major ballast shape, which was considered to be cubic when this frequency was less than 0.17. Then, by utilizing the Nataf transformation, a joint probability density function of the spherical harmonics frequencies was established, which enabled probabilistic characterization and rapid random generation of the 3D profiles of ballast particles. Direct shear tests revealed that the proposed model exhibited an error of only 3%; in contrast, models that employed spherical particles or single type particles showed substantially higher errors that could reach 66%. Finally, a ballast bed was developed to analyze the effects of elongated and flaky particles on the lateral resistance of ballast bed and the sleeper acceleration receptance. The results indicated that the lateral resistance could increase by up to 30% as the proportions of elongated and flaky particles increased, while the sleeper acceleration receptance below 500 Hz not changed.