Some response surface functions in complex engineering systems are usually highly nonlinear, unformed, and expensive to evaluate. To tackle this challenge, Bayesian optimization (BO), 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 (GP) surrogate model to be sporadic. In this article, 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 with 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 BO with GPs. Through extensive evaluations carried out on synthetic and real-world optimization problems, we show that AEN-RBF outperforms the 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. BO is one of the widely used strategies for learning the global optimum of black-box functions. BO has been applied to robotics, anomaly detection, automatic learning algorithm configuration, reinforcement learning, and deep learning. This article proposes one new kernel function, named after AEN-RBF. The new kernel function will make BO with GPs more robust to outliers and lower the data quality barrier of model training. This article was motivated by the hyperparameter 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.