## 1 Introduction  
Low Earth orbiters (LEOs) have been widely applied in Earth observation systems, such as remote sensing, ocean altimetry, atmosphere exploration, and Earth gravity field determination. The applications in these fields require high accuracy, reliability, and real-time performance of satellite orbits. Therefore, the precise orbit determination (POD) of LEOs is a popular topic. Using on-board GPS receivers for the POD of LEOs is the most effective method thanks to its low cost and high precision, and this method has been employed in many LEO missions, such as TOPEX/Poseidon [1], Jason-1/2 [2, 3, 4], CHAllenging Ministatellite Payload (CHAMP) [5, 6], Gravity field Recovery and Climate Experiment (GRACE) [7, 8], TerraSAR-X [9], MetOp-A [10], gravity field and steady-state ocean circulation explorer (GOCE) [11, 12], SWARM [13, 14, 15], and Sentinel-3A [16]. There are four classical orbit determination methods: the kinematic method, the reduced-kinematic method, the dynamic method, and the reduced-dynamic (RD) method. The kinematic method estimates the satellite position of each epoch only from GPS code and carrier-phase observations, and involves no orbit dynamic models [7, 11, 17, 18, 19, 20, 11]. Thus, the estimated precise kinematic orbit is often used to recover the Earth's gravity field [21, 22]. However, based on the kinematic method, the estimated orbits may be of poor quality or even impossible to determine because of unfavorable conditions such as data interruption, poor geometric conditions, and erroneous measurements [23]. In order to reduce or remove the \"spikes\" and \"jumps\" in the kinematic solutions, the reduced-kinematic method applies the constraints on the kinematic epoch to epoch position differences with respect to corresponding differences in the a priori dynamic orbit [17; 24]. The dynamic method makes use of [PERSON]'s second law of motion to establish the relationship between GPS observations and a satellite's orbital parameters (the initial conditions and dynamic parameters) with the known dynamic models. The known dynamic models include the forces that act on the satellite, such as Earth's gravity, tidal forces, luni-solar gravitation, atmospheric drag, and solar radiation pressure, which are used to predict the satellite's positions and enable the averaging of measurements from different epochs and the propagation of orbits across data gaps [23]. Thus, the orbits determined via the dynamic method are smoother than those determined via the kinematic approach [6; 23]. However, it is difficult to precisely obtain all various dynamic models (e.g., atmospheric drag) for LEOs because of the complex space environment, which limits the accuracy of the dynamic orbits to some extent. Moreover, an arc-length of 24 or 30 h is usually used for orbit integration in the dynamic method [11]; the unmodeled signals cannot be ignored for POD of LEOs. [PERSON]. [25]. and [PERSON] et al. [26] introduced the RD approach for POD of GPS satellites based on estimating empirical accelerations or pseudo-stochastic parameterization with velocity pulses. [PERSON] et al. [27] and [PERSON] et al. [28] proposed the RD approach for the POD of LEOs, in which pseudo-stochastic orbit parameters were introduced into the equation of motion to absorb the errors caused by the unmodeled dynamic signals. [PERSON] et al. [29], [PERSON] and [PERSON] [24], and [PERSON] et al. [23] made some improvements and employed the RD technique with CHAMP or GRACE data. [PERSON] et al. [6] presented and developed several methods for pseudo-stochastic orbit parameterization, such as instantaneous velocity changes (pulses), piecewise-constant accelerations, and continuous piecewise-linear accelerations. The RD approach, which combines the advantages of the kinematic positioning approach and the dynamic method, has been considered a powerful and efficient method for the POD of LEOs [30].