The simulation of spatially varying non-stationary seismic waves is a crucial consideration for the seismic analysis of engineering structures. The traditional approach often incurs substantial computational expenses and suffers from low efficiency due to the extensive implementation of large-scale matrix factorizations at numerous sampling points in both time and frequency domains. In this study, an enhanced approach is developed for the efficient simulation of non-stationary seismic waves. The core of this approach lies in the utilization of random frequency sampling coupled with time-domain interpolation to reduce matrix factorizations, and non-uniform fast Fourier transform for expediting the summation of trigonometric functions. A parametric analysis is conducted to investigate the effect of vital parameters, including the decoupling strategy of time-frequency spectrum, the number of random frequencies, and the number of interpolation knots over time, on both the efficiency and accuracy of the enhanced approach. Consequently, the preferred values for these parameters are determined. By considering the simulation of seismic waves acting on a group of bridges, the enhanced approach is validated in the view of correlation functions. The results demonstrate a high level of fidelity of the simulated seismic waves and proves the effectiveness of the enhanced approach in the efficient simulation of spatially varying non-stationary earthquake ground motions. The developed approach can be used for the simulation of non-stationary seismic waves for a wide range of engineering structures.