Selecting reasonable laser beam welding (LBW) process parameters is very helpful for obtaining a good welding bead profile and hence a high quality of the welding joint. Existing process parameter optimization approaches for LBW either based on cost-expensive physical experiments or low-fidelity (LF) computer simulations. This paper proposes a multi-fidelity (MF) metamodel based LBW process parameter optimization approach, in which different levels fidelity information, both from LF computer simulations and high-fidelity (HF) physical experiments can be integrated and fully exploited. In the proposed approach, a three-dimensional thermal finite element model is developed as the LF model, which is fitted with a LF metamodel firstly. Then, by taking the LF metamodel as a base model and scaling it using the HF physical experiments, a MF metamodel is constructed to approximate the relationships between the LBW process parameters and the bead profile. Two metrics are adopted to compare the prediction accuracy of the MF metamodel with the single-fidelity metamodels solely constructed with physical experiments or computer simulations. Results illustrate that the MF metamodel outperforms the single-fidelity metamodels both in global and local accuracy. Finally, the fast elitist non-dominated sorting genetic algorithm (NSGA-II) is used to facilitate LBW process parameter space exploration and multi-objective Pareto optima search. LBW verification experiments verify the effectiveness and reliability of the obtained optimal process parameters.