AbstractIn the probabilistic seismic performance assessment of structures, intensity measures (IMs) represent seismic characteristics and variations. Traditional fragility analysis method based on the assumption of linear regression requires selecting an optimal IM as input variable. By introducing machine learning (ML) techniques, nonparametric fragility analysis theoretically allows for considering all potential IMs as inputs. Nevertheless, to reduce input dimensionality and improve training efficiency, the feature selection of IMs remains imperative. This paper proposes a method to select optimal ground motion IMs for data‐driven surrogate modeling of structures. Specifically, the elastic net algorithm is employed to select the optimal multiple IMs based on the coefficient of determination and regression coefficient, differing from the efficiency and practicality emphasized in the traditional method. Using the optimal multiple IMs as input variables, several ML techniques are employed to construct surrogate models for seismic damage assessment of structures, thereby developing fragility functions, that is, the conditional probability of exceeding a damaged state given seismic intensity. A 3‐span, 6‐storey, reinforced concrete frame is utilized to illustrate the proposed methodology. The predictive performance of all ML models with the optimal multiple IMs outperforms that of the models with the commonly used IM (e.g., peak ground acceleration, PGA) as sole input and all candidate IMs as inputs. Additionally, the surrogate models with the optimal multiple IMs enable a more comprehensive seismic fragility modeling of structures under two or more IMs simultaneously, such as the fragility surface under spectral acceleration at 1.0s (Sa‐1.0s) and velocity spectrum intensity (VSI).