In this work, two dynamical models are formulated to describe the secular dynamics of navigation satellites moving in the medium Earth orbit (MEO) and geosynchronous orbit (GSO) regions. In the dynamical models, the leading terms of the Earth’s oblateness and the luni-solar gravitational perturbations are considered. For convenience, the orbits of the Sun and the Moon are described in the geocentric ecliptic reference frame, where the regression of nodal line and precession of apsidal line of the lunar orbit can be approximated as linear functions of time. The disturbing function acting on navigation satellites is analytically averaged over the mean motions of both the satellite and the third body (the Sun or the Moon). Explicit expressions of the averaged disturbing function are provided in the geocentric ecliptic and equatorial reference frames, corresponding to averaged model 1 and averaged model 2, respectively. It is found that there are seven resonant arguments in averaged model 1, while there are thirty-two resonant arguments in averaged model 2. The associated resonance curves corresponding to these resonant arguments in each model form the dynamical backbone in the phase space, organizing secular behavior of navigation satellites. At last, the averaged models are numerically compared to the associated non-averaged model, and simulation results indicate that (a) the averaged models formulated in the geocentric ecliptic and equatorial reference frames are identical, and (b) both of these two averaged models are applicable in predicting secular behavior of navigation satellites.