The multi-stage compound channel, which is a common pattern in natural alluvial rivers and the regulation projects of urban rivers, inevitably freezes in winter when it is situated in cold northern areas with high latitudes. Given that ascertaining the stage–discharge relationship for rivers is the foundation for the development of flood control schemes and water resources management, this study concentrates on proposing an analytical model for predicting the stage–discharge curves of multi-stage ice-covered compound channels. In deducing the analytical model, the cross section of the channel is first segmented into several homogeneous subregions that can be grouped into seven categories according to the geometric characteristics. Through analyzing the momentum transfer between adjacent subregions, the force balance equation for each subregion is then established to get the bulk mean velocity for the corresponding subregion, thereby obtaining the discharge by solving a tridiagonal matrix. Subsequently, measurements from two-stage and three-stage ice-covered compound channel experiments and three sets of experimental data from the literature are used to validate the performance of the proposed model. Good agreement between the predictions and the measured data suggests that the deduced model can accurately estimate the discharge for the multi-stage ice-covered compound channels when the flow depth is given. Finally, sensitivity analysis indicates that Manning's roughness coefficient of the channel bed has a more pronounced impact on the stage–discharge relationship than that of the ice cover. Moreover, when compared to the two-stage ice-covered compound channel, the multi-stage ice-covered compound channel offers greater potential for water resource utilization.