A transonic, high Reynolds number natural laminar flow airfoil is designed and studied. The γ-θ transition model is combined with the shear stress transport (SST) k-w turbulence model to predict the transition region for a laminar-turbulent boundary layer. The non-uniform free-form deformation (NFFD) method based on the non-uniform rational B-spline (NURBS) basis function is introduced to the airfoil parameterization. The non-dominated sorting genetic algorithm-II (NSGA-II) is used as the search algorithm, and the surrogate model based on the Kriging models is introduced to improve the efficiency of the optimization system. The optimization system is set up based on the above technologies, and the robust design about the uncertainty of the Mach number is carried out for NASA0412 airfoil. The optimized airfoil is analyzed and compared with the original airfoil. The results show that natural laminar flow can be achieved on a supercritical airfoil to improve the aerodynamic characteristic of airfoils.