The identified modal parameters in Operational Modal Analysis (OMA) are always subject to statistical uncertainties from many sources and the uncertainty is also important to assess the quality of the identified results. This study proposes an uncertainty quantification method in OMA for time-varying structures based on the Functional Series Time-dependent AutoRegressive Moving Average (FS-TARMA) model. A multi-stage pseudo-linear optimization scheme is first adopted to obtain the approximate Maximum Likelihood Estimates (MLEs) of FS-TARMA model parameters and the covariance is computed as the inverse Hessian matrix of the Negative Log-Likelihood Function (NLLF) at the MLEs. Finally, the approximate closed-form covariance of modal parameters is derived by the first-order sensitivity method. The proposed method is verified by two numerical examples and an experimental beam. Furthermore, the ability of the uncertainty to extinguish the spurious modes caused by colored noise excitation and the noise modes caused by measurement noise or inappropriate model structure is also demonstrated in these examples.