## 1 Introduction  
Low Earth Orbit (LEO) satellite is a valuable platform for Earth observation due to the low altitude, short period, and rapid connection with the ground [1]. The observation data collected by the LEO satellite have been widely used in geosciences with a wide area coverage, such as crustal deformation [2, 3], ocean altimetry [4], Earth gravity field recovery [5, 6], space weather [7, 8] and remote sensing [9, 10]. To ensure successful scientific applications undertaken by the LEO satellite, the precise orbit determination (POD) for the LEO satellite is essential. With the launch of the project Topex/Poseidon, developed jointly by the National Aeronautics and Space Administration (NASA) and the Centre National d'Etudes Spatiales (CNES) in 1992, satellite-borne GPS technology has been utilized to obtain a high-accuracy centimeter-level orbit. The on-board GNSS receivers have been applied to achieve high-accuracy orbital information in the subsequent LEO satellite missions. Kinematic POD (KPOD) of LEO satellites does not rely on any orbital dynamics model and can achieve high-precision results. Many scholars conducted high-precision KPOD. For instance, [PERSON] et al. [11] applied undifferenced KPOD to the CHAMP satellite-borne GPS data, and the root mean squares (RMS) of the results achieved centimeter-level in all Radial, Tangent, and Normal (RTN) directions when compared to dynamic orbit. [PERSON] [12] performed a KPOD solution based on GRACE satellite-borne GPS data of 1 day and realized a centimeter-level orbit accuracy in all RTN directions. [PERSON] and [PERSON] [13] proposed a KPOD solution based on the precise point positioning (PPP) technique for the Swarm satellites, and the accuracy of the solution was consistent with the precise orbit products published by the European Space Agency (ESA). The KPOD of LEO satellites has been widely used in many fields, particularly for global gravity field recovery [14]. With the successful implementation of gravity satellite missions, such as CHAMP (CHAllenging Minisatellite Payload), GRACE (Gravity Recovery And Climate Experiment), and GOCE (Gravity field and steady-state Ocean Circulation Explorer), the recovery in the Earth's gravity field based on precise kinematic orbits of LEO satellites has become a major method. [PERSON] et al. [15] inverted the 60 th-order global time-varying gravity field by using the satellite-borne data of CHAMP. [PERSON] et al. [16] recovered the global gravity field with the 18 months' precision kinematic orbit of Swarm satellites, which could be able to fill the data gap between GRACE and GRACE-FO (GRACE Follow On). [PERSON] et al. [17] calculated the Tongji-Grace 2018 monthly gravity field model based on the GRACE's precise kinematic orbit by using the optimized short-arc method. [PERSON] et al. [18] discuss the possibility of time-variable gravity field detection by using GOCE gravity gradient observations, indicating the potential of the superconducting gravity gradiometer in time-variable gravity detection. In recent years, multi-satellite joint applications have gradually become a main trend due to the advantage of LEO satellites in many aspects. For example, due to the influence of the inclination of the satellite orbit, the re-entry period of the satellite, and the lack of on-board data, the single-satellite-based time-varying gravity field is limited in terms of accuracy, spatial resolution, and temporal resolution [19]. However, this situation can be improved by the multi-satellite observations, which can increase the coverage and density of satellite data. Additionally, when the data from one of the satellites is missing, the data from the other satellites can be well supplemented. However, most of the precise LEO satellite orbits are processed and released by different agencies with different GNSS precise orbit and clock products, which will cause systematic errors in multi-satellite-based applications.