Bayesian formulations of inverse problems are attractive for their ability to incorporate prior knowledge and update probabilistic models as new data become available. Markov chain Monte Carlo (MCMC) methods sample posterior probability density functions (pdfs) but require accurate prior models and many likelihood evaluations. Dimensionality-reduction methods, such as principal component analysis (PCA), can help define the prior and train surrogate models that efficiently approximate costly forward solvers. However, for problems like full waveform inversion, the complex input/output relations often cannot be captured well by surrogate models trained only on prior samples, leading to biased results. Including samples from high-posterior-probability regions can improve accuracy, but these regions are hard to identify in advance. We propose an iterative method that progressively refines the surrogate model. Starting with low-frequency data, we train an initial surrogate and perform an MCMC inversion. The resulting posterior samples are then used to retrain the surrogate, allowing us to expand the frequency bandwidth in the next inversion step. Repeating this process reduces model errors and improves the surrogate's accuracy over the relevant input domain. Ultimately, we obtain a highly accurate surrogate across the full bandwidth, enabling a final MCMC inversion. Numerical results from 2D synthetic crosshole Ground Penetrating Radar (GPR) examples show that our method outperforms ray-based approaches and those relying solely on prior sampling. The overall computational cost is reduced by about two orders of magnitude compared to full finite-difference time-domain modeling.