An improved method, termed the quasi-point-source method, was proposed to improve the performance of the Green’s function method in simulating the non-stationarity of ground motion time histories. The ratio of the amplitude spectrum with consideration of the fault rupture process to that without consideration was used as a correction coefficient to modify the point-source model. Then the modified point-source model was used to calculate the ground motion amplitude spectrum. The ground motion phase spectrum was expressed as the product of all-pass function and minimum-phase shift function based on the decomposition theory of causal functions. All-pass function and minimum-phase shift function were determined by the phase characteristics of historical records and the causality constraint, respectively. Then the ground motion time histories were obtained from the above amplitude spectrum and phase spectrum through inverse Fourier transform. The method was verified by a hypothetical earthquake and the 2024 Hualien earthquake in Taiwan Province. Results showed that the simulation results could reflect the directivity effect of ground motion, and the non-stationarity of the simulated ground motion time histories accorded with that of historical records. The response spectrum simulation error (Misfit) over the entire period range within 0–10 s was less than 0.3 for 30 stations, and the standard deviation was less than 0.15.