In this paper, a multi-fidelity Bayesian optimization approach is presented to tackle computationally expensive constrained multi-objective optimization problems (MOPs). The proposed approach consists of a three-stage optimization framework designed to search for promising candidate points. In the first stage, an acquisition function is proposed to identify a feasible solution if none are available in the current set of sampling points. Subsequently, a new multi-fidelity weighted expected hypervolume improvement function is developed to find better solutions. In the third stage, a constrained weighted lower confidence bound acquisition function is presented to enhance the constraint predictions and refine the solutions near the constraint boundary. Additionally, a filter strategy is suggested to determine whether constraint updating is necessary, aiming to save computational resources and improve optimization efficiency. Moreover, to expedite the optimization process, a parallel optimization approach is further developed based on the suggested three-stage optimization framework. To achieve this, a multi-fidelity influence function is introduced, allowing the proposed approach to determine a desired number of candidate points within a single iteration. Lastly, the proposed approach is demonstrated through six numerical benchmark examples, which verifies its significant advantages in addressing expensive constrained MOPs. Besides, the proposed approach is applied to the multi-objective optimization of a metamaterial vibration isolator, resulting in the attainment of satisfactory solutions.