Abstract: Quantification of structural damage states in the fragility analysis and seismic performance evaluation of underground structures usually adopts the deterministic threshold values, which is a limitation in the framework of performance-based seismic engineering. A new seismic fragility analysis method considering the fuzzy probability of damage states of the underground structures is proposed based on the fragility analysis framework of the traditional underground structures. A two-dimensional nonlinear finite element model for a three-story and three-span subway station structure is established fully considering the soil-structure interaction (SSI). The input ground motions used in the nonlinear incremental dynamic analysis of the SSI model are back-calculated from an ensemble of 21 actual earthquake records on the ground surfaces. The triangle and quasi-normal membership functions are introduced to compute the membership degrees of the maximum interlayer drift ratios for different damage states so as to derive the fuzzy failure probability under different intensity levels of ground motions. The maximum likelihood estimation method is used to fit the numerical results and establish the seismic fragility curves of subway station structures considering the fuzzy probability of damage states. The results of the numerical study indicate that the seismic fragility curves of the underground structures derived from the proposed method generally give a safer estimation of the seismic performance of the structures. The dispersions of the seismic fragility curves increase with the consideration of the fuzzy probability of damage states. Besides, the influences of the membership function types on the results of the fuzzy seismic fragility analysis can be ignored.