The structure of L625 Rayleigh lidar system, principles of atmospheric density measurement which is based on the Rayleigh scattering theory and iterative method to retrieve middle atmosphere density profiles were described respectively. In order to prove the effectiveness of iterative algorithm, simulation analysis was made. Error analysis indicated that contributions to the uncertainty of the retrieved molecular number density mainly derived from signal to noise ratio of the return signals and molecular number density at the reference point. Meanwhile, main error sources of the return signals and method to determine the reference point were given. Finally, according to 16 years’ routine observations of L625 Rayleigh lidar, monthly distribution and annual average of molecular number density from 1997 to 2004 and 2006 to 2012 in the middle atmosphere were presented. Results show that molecule number density distribution is sparse in winter and it becomes thick in summer. The yearly distribution of molecular number density seems to be relatively stable. Comparison with the model density which is obtained from 1976 USA standard atmospheric model indicates that density ratio profile of L625 Rayleigh lidar to the model density fluctuates between 1.05 and 1.13.