Abstract: The epistemic uncertainty of source exists objectively and has a significant impact on the characteristics of near-surface ground motions. In this paper, the multiplicative dimensional reduction method (M-DRM) is applied to solve the ground motion variability of complex sites near faults considering the uncertainty of the source. The uncertainty analysis problem is converted into a finite deterministic analysis to obtain statistical moments of ground motion parameters consistent with Monte Carlo simulation. The deterministic analysis uses the boundary element method to simulate the entire physical process. Based on this method, the mountain-canyon site near a strike-slip fault was discussed as an example. The spatial distribution variability of the peak acceleration (PGA) and peak velocity (PGV) under the coupling of near-fault effect, site effect and source epistemic uncertainty was analyzed, as well as statistical values of spectral acceleration (SA) and velocity large pulse at some surface points. The results show that the M-DRM is suitable for solving random ground motions in complex sites with near faults with uncertainty of source parameters. The scattering effect of the mountain-valley topography on seismic waves and the near-fault effect caused the PGA and PGV means of the hanging wall of the fault increase significantly, and the increase rate could reach 2 times. The uncertainty is transmitted in the complex site, and the transmission characteristics are related to the frequency of seismic waves, resulting in an amplified PGA variability and a decreasing PGV variability. The influence on the peak value of the velocity pulse can reach more than 25%. The influence on the structural response cannot be ignored, and the large-span structures should also consider the spatial distribution difference of ground motion variability.