The terrain slope is one of the most important surface characteristics for quantifying the Earth surface processes. Space-borne LiDAR sensors have produced high-accuracy and large-area terrain measurement within the footprint. However, rigorous procedures are required to accurately estimate the terrain slope especially within the large footprint since the estimated slope is likely affected by footprint size, shape, orientation, and terrain aspect. Therefore, based on multiple available datasets, we explored the performance of a proposed terrain slope estimation model over several study sites and various footprint shapes. The terrain slopes were derived from the ICESAT/GLAS waveform data by the proposed method and five other methods in this study. Compared with five other methods, the proposed method considered the influence of footprint shape, orientation, and terrain aspect on the terrain slope estimation. Validation against the airborne LiDAR measurements showed that the proposed method performed better than five other methods (R2 = 0.829, increased by ~0.07, RMSE = 3.596°, reduced by ~0.6°, n = 858). In addition, more statistics indicated that the proposed method significantly improved the terrain slope estimation accuracy in high-relief region (RMSE = 5.180°, reduced by ~1.8°, n = 218) or in the footprint with a great eccentricity (RMSE = 3.421°, reduced by ~1.1°, n = 313). Therefore, from these experiments, we concluded that this terrain slope estimation approach was beneficial for different terrains and various footprint shapes in practice and the improvement of estimated accuracy was distinctly related with the terrain slope and footprint eccentricity.