## 1 Introduction  
Complicated mechanical assets have served an increasingly significant role in industrial systems, which tend to accomplish diverse functions and operate in comprehensive and dynamic environment load, particularly restricted by multiple performance indicators, resulting in various probable failure modes. In addition, they are prone to induce downtime, trigger large economic losses, and even threaten personnel safety once failure occurs [1; 2; 3]. Notably, these multiple performance parameters, when operating, exhibit random uncertainties, such as variations in material properties, loads and structural geometry, which can produce inevitable effects on the reliability and safety analysis of complicated mechanical assets [4; 5; 6; 7]. Furthermore, there exists more than one limit state function, and basic component variables vary with the different limit state functions induced by multiple failure modes. Consequently, it is extraordinarily significant to investigate the performance reliability analysis oriented to complicated mechanical systems, promoting operation efficiency in the economy, safety and other terms [8; 9; 10]. As aforementioned, the practical limit state function is generally characterized by being high-dimension, nonlinear and implicit, which generates a challenge in calculation. To this end, the Monte Carlo simulation earns a wide utilization [11; 12], whose advantages encompass the execution convenience and excellent accuracy based on sufficient samples. However, the calculation efficiency of the Monte Carlo method will be seriously restricted when the limit state function is nonlinear, complex or implicit, especially for the system with low failure probability, indicating enormous computation cost and time [13; 14; 15]. Therefore,how to promote the calculation efficiency of the limit state function has become a widely concerned point in the reliability analysis of complicated mechanical assets. In recent decades, the surrogate model has been widely employed to address the relevant issues, whose core is constructing a surrogate model to substitute the original limit state function during calculation. Guided by the mentioned principle, quantities of approaches arise, such as response surface method (RSM) [16; 17; 18], Kriging [19; 20; 21], polynomial chaos expansion (PCE) [22; 23], radial basis function (RBF) [24; 25], support vector machine (SVM) [26; 27], artificial neural network (ANN) [28; 29; 30; 31; 32; 33], etc., all of which are relatively convenient to be enforced and to promote the computation efficiency to a great extent. More specifically, the Kriging model possesses superior local characteristics, enabling it to concurrently predict sample response value and corresponding errors. Nonetheless, the intricate nature of its construction poses challenges, particularly in scenarios involving vast sample sizes, where a significant computational load will be incurred. Additionally, it behaves poorly in the fitting procedure for high-dimensional problems. The PCE model possesses excellent global characteristics, but its computation efficiency decreases exponentially with the increase in the dimension upon input variables. The RBF model has a distinct advantage in simple structure, concise training, fast and convergence aspects, and can approximate any nonlinear function, overcoming local minimum issues. Nevertheless, it falls short in interpretation, determining the number of hidden layer nodes, node center and width, and data pathology occurs in the optimization process. By contrast, the SVM model exhibits superior capabilities in handling high-dimensional data and possesses robot generalization abilities, but the computational complexity of the SVM model is high, especially for large-scale and high-dimensional data sets, leading to substantial demands on computation time and space [34]. Moreover, the training of the SVM algorithm necessitates numerous iterations, escalating computational complexity.