Accurate prediction of impact forces arising from cavitation bubble-cloud collapse is crucial for biomedical technologies, particularly non-invasive ultrasound-based therapies including histotripsy. When collapse occurs near a solid boundary, complex multiscale phenomena—including rapid shock-wave emission and re-entrant jet formation—emerge and necessitate highly resolved computational fluid dynamics (CFD) simulations. The associated computational cost, however, renders parameter studies across an unbounded configuration space impractical. To address this challenge, we propose a hierarchical multi-fidelity strategy where fidelity levels are defined by grid resolution or, more broadly, by distinct data sources. The long-term vision is to leverage single-bubble simulations as low-fidelity surrogates for costly bubble-cloud simulations. We introduce a Multi-Fidelity Variational Heteroscedastic Gaussian Process (MFVHGP) framework to predict peak wall pressures, which explicitly accommodates heteroscedasticity, a critical requirement given that uncertainty in cavitation dynamics exhibits extreme variance during the collapse stage. Inter-fidelity correlations are captured via a product kernel and the Linear Model of Coregionalization (LMC), with fidelity defined by mesh resolution. We present an initial validation of this framework using multi-resolution single-bubble data. Training data are generated using the high-order compressible solver JAX-Fluids, while varying the ambient pressure and the non-dimensional stand-off distance. Computational efficiency is achieved by coupling abundant low-fidelity simulations on coarse grids (20 cells/$R_0$) with a limited set of high-fidelity simulations on fine grids (50 cells/$R_0$). The results demonstrate that the MFVHGP effectively learns inter-mesh correlations and reconstructs the high-fidelity response using fewer high-fidelity samples than single-fidelity Gaussian processes. By augmenting scarce high-fidelity data with inexpensive low-fidelity simulations, the framework improves predictive accuracy and captures key physical transitions, including nonlinear variations in impact intensity as the bubble approaches the wall. Overall, the proposed multi-fidelity approach provides a robust, scalable surrogate model for impact-load prediction and establishes the foundation for extending data-driven surrogates to more complex, multi-bubble collapse regimes.