Gaussian Processes (GP) is an effective surrogate model for globally accurate emulation of noisy computer simulations. With the goal of building surrogate model for expensive simulations with heterogeneous noise, we utilize a Bayesian framework paired with GP, from which a novel empirical integrated mean squared error-based sequential design strategy is proposed to approximate the mean response surface with a small fixed simulation budget. Through different synthetic examples, we show that the proposed strategy has the potential to achieve high predictive accuracy under a small budget compared to existing state-of-the-art methods. We also demonstrate the performance of our strategy on a real-data simulation of finding the reliable region in podium building seismic design.