Estimating rare events in structural engineering demands intensive computational resources, and the complexity of the performance function poses challenges to accurate probability estimation. To efficiently conduct the structural reliability analysis on rare events, we propose a novel variational Bayesian Monte Carlo (VBMC)-based subset simulation (SS) method. The proposed approach can be divided into two aspects: the Monte Carlo simulation (MCS) is first applied to estimate the first-layer unconditional failure probability, then VBMC is employed to estimate subsequent conditional failure probabilities. Finally, the original small failure probability is obtained by the product of the first-layer unconditional failure probability and a series of conditional failure probabilities. The proposed approach inherits the merits of SS as well as VBMC, which converts the tricky rare event estimation into a series of high-frequency event estimations and leverages VBMC instead of the canonical Markov chain to produce the independent and informative next-generation failure-conditional samples. Four case studies, including high-dimensional and discrete state space scenarios, are performed to illustrate the feasibility and generality of the proposed method.