Compared with conventional optical remote sensing data, which contain only multispectral information and can only realize two-dimensional land cover classification, the advantage of airborne multispectral light detection and ranging (MS-LiDAR) is that it contains both multispectral and spatial information and can realize three-dimensional land cover classification. However, the existing land cover classification methods for airborne MS-LIDAR data require too high feature dimension to distinguish all kinds of objects simultaneously and have high algorithm complexity. So, a stepwise classification algorithm combining spatial and normalized difference ratio index (NDRI) features is proposed. Firstly, the multi-band independent point clouds of airborne multispectral LiDAR are merged to obtain the merged point cloud data with spatial location and their multispectral information. Secondly, based on the elevation consistency of urban ground spatial adjacent points, a ground filtering algorithm under spatial neighborhood growth is used to separate ground and non-ground points. Thirdly, based on the difference of laser reflectance characteristics of different objects, the NDRI index is designed to separate the grass (tree) from the ground (non-ground), and the adaptive optimal NDRI index under the principle of maximum inter-class variance is used to achieve the fine classification of ground and non-ground points. Finally, 3D majority voting is used to alleviate the noise in the previous classification result in order to further optimize classification result. The proposed algorithm makes comprehensive use of the spatial and multispectral features contained in multispectral LiDAR data, and the step-by-step processing method is more convenient to design simple and effective rules according to the characteristics of specific separation objects. The effectiveness and feasibility of the proposed algorithm are tested by using Optech Titan airborne multi-spectral LiDAR data of different scenes. The experimental results show that: (1) The average overall accuracy and Kappa coefficient of the proposed algorithm can reach 90.17% and 0.861, this demonstrates that the proposed algorithm can realize the accurate three-dimensional land cover classification of the multi-spectral LiDAR data in urban areas. (2) The step-by-step processing method adopted by the proposed algorithm is more convenient for designing simple and effective rules according to the characteristics of specific separation objects. The algorithm design is simple and the complexity is low. (3) The availability of NDRI features designed to distinguish between grass (trees) and road (buildings) can provide theoretical support for the design and selection of salient features of other machine learning algorithms.