Abstract: Lunar satellite is the most popular concept in the recent and future lunar exploration missions, such as Clementine (NASA, 1994), Lunar Prospector (NASA, 1998) and Smart-1 (ESA, 2002). The scientific goals of these missions, such as the global mapping, the determination for the high resolution of the lunar gravitational potential and the lunar mass, require the precise orbit determination support. Comparing with that of the Earth satellite, the motion of the lunar satellite could still be regarded as a perturbed two-body problem. However, it is complicated by the fact that the Moon, as the central body of the satellite, is different from the Earth, such as it is rotating slowly, it has no atmosphere, and its harmonic coefficient of the oblate, J₂, isn't much larger than other harmonic coefficients. All of these bring us both the advantage and the disadvantage in building an analytical perturbing theory for the lunar satellite. And the corresponding analytical solution of the motion is not the same as that of the Earth satellite. It is not convenient to use the analytical theory directly during the procedure of orbit determination. A numerical method that combines with some properties of the analytical method seems to be a wise choice, which is investigated in this paper. Since the lunar satellite is usually in the polar orbit for finding the water, we use a set of non singular orbital elements to discuss the precise orbit determination for lunar satellite, which is suitable for the lunar orbiter with 0≤e <1 and 0°≤ i < 180°. Then we could take the most advantage of both the analytical method and the numerical method. An overall strategy and make a simulation with the range data is given in the paper. The result is in agreement with the theoretical analysis and suggests that this strategy have a practical prospect.