A frozen orbit is beneficial for observation owing to its stationary apsidal line. The traditional gravitational field model of frozen orbits only considers the main zonal harmonic terms $J_2$ and limited high-order terms, which cannot meet the stringent demands of all missions. In this study, the gravitational field is expanded to $J_{15}$ terms and the Hamiltonian canonical form described by the Delaunay variables is used. The zonal harmonic coefficients of the Earth are chosen as the sample. Short-periodic terms are eliminated based on the Hori–Lie transformation. An algorithm is developed to solve all equilibrium points of the Hamiltonian function. A stable frozen orbit with an argument of perigee that equals neither 90 ${}^{\circ }$ nor 270 ${}^{\circ }$ is first reported in this paper. The local stability and topology of the equilibrium points are obtained from their eigenvalues. The bifurcations of the equilibrium points are presented by drawing their global long-term evolution of frozen orbits and their orbital periods. The relationship between the terms of the gravitational field and number of frozen points is addressed to explain why only limited frozen orbits are found in the low-order term case. The analytical results can be applied to other Earth-like planets and asteroids.