Combining intrusive geotechnical site investigations with non-intrusive geophysical surveys is a cost-effective approach to producing data with varying levels of accuracy, uncertainty, and different spatial scales to better characterize the site's liquefaction properties. Moreover, the demand for three-dimensional (3D) subsurface models in geotechnical engineering is increasing, but the models contain uncertainties and spatial variability associated with the use of relevant stochastic and geostatistical methods, and it remains a challenging task to obtain reliable liquefaction assessment results and the corresponding damage capacity. This study proposes a data-driven and non-parametric form of 3D multi-source fusion Bayesian compressive sampling (3D MSF-BCS) method for assessing 3D soil liquefaction-induced damage capacity. It consists of three main components: (i) 3D MSF-BCS fusing sparse geotechnical data (e.g., cone penetration test (CPT)) and geophysical data (e.g., multichannel analysis of surface waves (MASW)) for 3D site modeling, (ii) quantifying the accuracy and uncertainty of 3D MSF-BCS, and (iii) Soil liquefaction-induced damage capacity analysis in 3D space. The method was applied to numerical examples and a real case study at the Cresselly Place site, and the results showed that the proposed method performs well.