The full waveform decomposition technique is significant for LiDAR ranging. It is challenging to extract the parameters from non-Gaussian shaped waveforms accurately. Many parametric models (e.g. the Gaussian distribution, the lognormal distribution, the generalized normal distribution, the Burr distribution, and the skew-normal distribution) were proposed to fit sharply-peaked, heavy-tailed, and negative-tailed waveforms. However, these models can constrain the shape of the waveform components. In this article, the Gaussian convolution model is established. Firstly, a set of Gaussian functions is calculated to characterize the system waveform so that asymmetric and non-Gaussian system waveforms can be included. The convolution result of the system waveform and the target response is used as the model for fitting the overlapped echo. Then a combination method of the Richardson–Lucy deconvolution, layered iterative, and Gaussian convolution is introduced to estimate the initial parameters. The Levenberg–Marquardt algorithm is used for the optimization fitting. Through experiments on synthetic data and practical recorded coding LiDAR data, we compare the proposed method with two decomposition approaches (Gaussian decomposition and skew-normal decomposition). The experiment results revealed that the proposed method could precisely decompose the overlapped non-Gaussian heavy-tailed waveforms and provide the best ranging accuracy, component fitting accuracy, and anti-noise performance. However, the traditional Gaussian and skew-normal decomposition methods can not fit the components well, resulting in inaccurate range estimates.