In order to modify the rail pre-grinding profile smoothly, non-uniform rational B-spline (NURBS) curve with weight factors is used to establish a parameterized model of the profile. A wheel-rail contact stochastic finite element model (FEM) is constructed by the Latin hypercube sampling method and 3D elasto-plastic FEM, in which the wheelset’s lateral displacement quantity is regarded as a random variable. The maximum values of nodal accumulated contact stress (NACS) and nodal mean contact stress (NMCS) in different pre-grinding profiles with differential weight factors are calculated and taken as the training samples to establish two Kriging models. A multi-objective optimization model of pre-grinding profile is established, in which the objective functions are the NACS and NMCS Kriging models. The optimum weight factors are sought using a non-dominated sorting genetic algorithm II (NSGA-II), and the corresponding optimum pre-grinding profile is obtained. The contact stress calculation before and after optimization indicates that the maximum values of NACS and NMCS decline significantly.