In recent years, with the reduction of the cost of microsatellites, the development of commercial rockets and the multi-satellite launching technology, the construction of large-scale constellations in low-Earth orbit (Mega-Constellations) has become a development trend. Since the motion of LEO satellites is affected by perturbations such as non-spherical gravitational fields and atmospheric drag, as well as the uncertainty of actuators, measurement systems, and dynamic models, it is easy to cause divergence of constellation configurations. The station-keeping control of the satellites is crucial for the stable operation of the mega-constellation. Aiming at this problem, this paper proposes an uncertainty propagation approach based on semi-analytical and Monte Carlo for LEO Mega-Constellations. Under the assumption that initial uncertainty on the osculating trajectory is Gaussian distribution, through hypothesis testing analysis, the uncertainty propagation simulations of a single satellite suggest that the satellite argument of latitude and the relative phase of co-plane satellites can be both considered as Gaussian distributions with zero means. Multi-group Monte Carlo simulations with product-based least-squares surface fitting establish an approximate mapping between initial and terminal errors. The mapping provides an efficient method for deviation prediction and can be used to design the station-keeping control strategy.