In recent years, Full-Waveform Inversion (FWI) has been extensively used to derive high-resolution subsurface velocity models from seismic data. However, due to the nonlinearity and ill-posed nature of the problem, FWI requires a good starting model to avoid producing non-physical solutions. Moreover, conventional optimization methods fail to quantify the uncertainty associated with the recovered solution, which is critical for decision-making processes. Bayesian inference offers an alternative approach as it directly or indirectly evaluates the posterior probability density function. For example, Markov Chain Monte Carlo (MCMC) methods generate multiple sample chains to characterize the solution's uncertainty. Despite their ability to theoretically handle any form of distribution, MCMC methods require many sampling steps; this limits their usage in high-dimensional problems with computationally intensive forward modeling, as is the FWI case. Variational Inference (VI), on the other hand, provides an approximate solution to the posterior distribution in the form of a parametric or non-parametric proposal distribution. Among the various algorithms used in VI, Stein Variational Gradient Descent (SVGD) is recognized for its ability to iteratively refine a set of samples to approximate the target distribution. However, mode and variance-collapse issues affect SVGD in high-dimensional inverse problems. This study aims to improve the performance of SVGD within the context of FWI by utilizing, for the first time, an annealed variant of SVGD and combining it with a multi-scale strategy. Additionally, we demonstrate that Principal Component Analysis (PCA) can be used to evaluate the performance of the optimization process. Clustering techniques are also employed to provide more rigorous and meaningful statistical analysis of the particles in the presence of multi-modal distributions.