Some response surface functions in complex engineering systems are usually highly nonlinear, unformed, and expensive-to-evaluate. To tackle this challenge, Bayesian optimization, which conducts sequential design via a posterior distribution over the objective function, is a critical method used to find the global optimum of black-box functions. Kernel functions play an important role in shaping the posterior distribution of the estimated function. The widely used kernel function, e.g., radial basis function (RBF), is very vulnerable and susceptible to outliers; the existence of outliers is causing its Gaussian process surrogate model to be sporadic. In this paper, we propose a robust kernel function, Asymmetric Elastic Net Radial Basis Function (AEN-RBF). Its validity as a kernel function and computational complexity are evaluated. When compared to the baseline RBF kernel, we prove theoretically that AEN-RBF can realize smaller mean squared prediction error under mild conditions. The proposed AEN-RBF kernel function can also realize faster convergence to the global optimum. We also show that the AEN-RBF kernel function is less sensitive to outliers, and hence improves the robustness of the corresponding Bayesian optimization with Gaussian processes. Through extensive evaluations carried out on synthetic and real-world optimization problems, we show that AEN-RBF outperforms existing benchmark kernel functions.Note to Practitioners-Some industrial systems cannot be accurately represented by physical models. In this situation, data-driven black-box optimization is necessary for advancing the system automation and intelligence. Bayesian optimization is one widely used strategy for learning the global optimum of black-box functions. Bayesian optimization has been applied to robotics, anomaly detection, automatic learning algorithm configuration, reinforcement learning, and deep learning. This paper proposes one new kernel function, named after AEN-RBF. The new kernel function will make Bayesian optimization with Gaussian processes more robust to outliers and lower the data quality barrier of model training. This paper was motivated by the hyperparameters tuning problem of deep learning models for image defect detection in advanced manufacturing, but the method can be easily extended to other applications where kernel functions are needed. Our proposed method is verified by synthetic and real-world optimization problems.