## 1 Introduction  
As GPS carrier phase observations are much more accurate than pseudorange observations, carrier phase measurements are usually used as observations for the precise GPS-derived orbit determination of low-Earth orbit (LEO) satellites [1; 2; 3; 4; 5]. Compared to ground GPS data, due to the high-speed motion and complicated observation environment of LEO satellites, cycle slips occur frequently in LEO satellite GPS carrier phase data, which are undesirable and must be detected and repaired correctly to maintain the continuity of GPS carrier phase observations. The difficulties of cycle slip detection are as follows: (i) the geometric distance between receiver and satellite varies greatly; (ii) the ionospheric delay varies quickly. How to eliminate the influence of geometric distance and ionospheric delay, and construct an effective quantity that is sensitive to the cycle slips is the key point of cycle slip detection. In the past, a variety of cycle slip detection methods for dual-frequency GPS signals were proposed. According to the strategy to eliminate the geometric distance, these methods can be divided into two categories, described below. The first category encompasses linear combination methods, where the linear combinations of observations at different frequencies are formed to eliminate the geometric distance, including the Hatch-Melbourne-Wubbena (HMW) combination [6; 7; 8] and the geometry-free phase combination [9]. The HMW combination is widely used to detect cycle slips for dual-frequency GPS data of both static and kinematic receivers. However, this method is not suitable for data with the same size of cycle slips on both \(L_{1}\) and \(L_{2}\) frequencies. In order to detect all potential cycle slips, a complementary method must be adopted, such as the traditional TurboEdit method [10]. However, TurboEdit may fail to detect small cycle slips of 1-2 cycles due to the large pseudorange noise [11]. The geometry-free phase combination is also usually used to detect cycle slips along with the HMW combination. This combination, however, contains ionospheric residuals and has a wavelength of only 5.4 cm for cycle slips, which is only effective for the cases when the ionosphere is stable. Some modifications of the phase combination with different techniques were also proposed, such as polynomial fitting [12; 13], Kalman filter [14; 15; 16], higher-order time difference [17; 18], cycle slip detection by integrating an inertial navigation system and GPS [19; 20], double- or triple-epoch difference residuals [21], and Bayesian theory [22]. In addition, [PERSON] [23] applied the ionospheric total electron content rate (TEGR) to detect cycle slips. [PERSON] et al. [11] proposed a method based on a forward- and backward-moving window averaging algorithm integrated with a second-order time-difference phase ionospheric residual (STPIR) algorithm. As the backward-moving technique is used, it is utilized for post-processing. Having smooth phase observations is an underlying assumption of many of these cycle slip detection methods. This assumption can be violated under high ionospheric activities [24]. The second category encompasses estimation methods, where the first-order or second-order time-difference observations of carrier phase and pseudorange observations are formed; then, the position variations are estimated along with the receiver clock variation, the ionospheric delay variations, and the cycle slips [25; 26]. Due to the numerous unknown parameters, the introduction of reasonable constraints of the ionospheric delay variations is necessary. Incorrect cycle slip determination may be caused when there are irregular rapid ionospheric variations. [PERSON] et al. [27] developed a geometry-based ionosphere-weighted (GBIW) model dedicated to efficiently estimating cycle slips in scenarios with active ionospheric conditions.