Conventional land cover classification methods based on airborne multispectral light detection and ranging (MS-LiDAR) data have insufficient capability for the cooperative utilization of spatial-spectral information or too high dimensions of features in the joint utilization of various features. This study proposed a spatial-spectral joint segmentation algorithm for airborne MS-LiDAR point clouds based on the multivariate Gaussian mixture model (GMM). First, radiometric correction, anomaly removal, and data fusion were performed for the original multi-band independent point clouds, forming multispectral point clouds that presented spatial locations and their multi-band spectral information. Then, spatial-spectral feature vectors were constructed using the extracted multispectral and elevation features of laser points. Meanwhile, the unit and scale differences among different types of features were eliminated through feature normalization and discretization. Subsequently, a GMM was built to fit the multimodal distribution of objects in the spatial-spectral feature space. Accordingly, the response levels of laser points to various objects were obtained, and the classification of various objects was determined according to the principle of maximum responsiveness. Finally, a 3D majority voting method was designed to optimize the segmentation results. The effectiveness and feasibility of the proposed algorithm were verified through experiments based on surveyed Optech Titan MS-LiDAR data. The experimental results show that the multivariate GMM combined with multi-band intensity features and elevation features yielded an overall segmentation accuracy of 93.57% and a Kappa coefficient of 0.912. The results also indicate that the high-accuracy segmentation of MS-LiDAR point clouds can be achieved by only combining four-dimensional features. This study provides a new approach for comprehensively utilizing the multispectral and spatial information in MS-LiDAR data.