## 2 Materials and Methods
### Second-Order Time Difference of the Satellite's Position  
Based on a reformulation of [PERSON]'s equation of motion, the second-order time difference of the satellite's position \(\Delta\Delta\overset{\rightarrow}{r}(t)\) can be given as by [PERSON] and [PERSON] [33].  
\[\begin{array}{ll}\Delta\Delta\overset{\rightarrow}{r}(t)&=\overset{ \rightarrow}{r}(t+\Delta t)-2\overset{\rightarrow}{r}(t)+\overset{ \rightarrow}{r}(t-\Delta t)\\ &=\int\limits_{t-\Delta t}^{t+\Delta t}(\Delta t-|\tau-t|)\overset{\rightarrow }{r}(\tau)d\tau\\ &=\Delta t^{2}\int\limits_{-1}^{1}(1-|\tau^{\prime}|)\overset{\rightarrow}{r}(t +\tau^{\prime}\Delta t)d\tau^{\prime}\\ &=\Delta t^{2}\int\limits_{-1}^{1}\overset{\rightarrow}{r}(t+\tau^{\prime} \Delta t)d\tau^{\prime}\end{array}, \tag{8}\]where  
\[\kappa=\left\{\begin{array}{ll}1-\tau^{\prime}&0\leq\tau^{\prime}\leq 1\\ 1+\tau^{\prime}&-1\leq\tau^{\prime}<0\end{array}\right.,0\leq\kappa\leq 1\]  
\[\tau^{\prime}=\frac{\tau-t}{\Delta t},\;\;\;t-\Delta t\leq\tau\leq t+\Delta t\;\;.\]  
\(\Delta\Delta(\star)=(\star)(t+\Delta t)-2(\star)(t)+(\star)(t-\Delta t)\) is the second-order time-difference (ST) operator, in which \((t+\Delta t),(t)\) and \((t-\Delta t)\) refer to three consecutive epochs, and \(\overrightarrow{r}(t)\) and \(\vec{\dot{r}}(t)\) are the satellite's position and acceleration at time epoch \(t\).  
Equation (8) is considered as the basic constraint in our work, which shows that the STP value \(\Delta\overrightarrow{r}(t)\) is only related to the acceleration \(\vec{\dot{r}}\) and sampling interval \(\Delta t\). Since \(\int\limits_{-1}^{1}\kappa d\tau^{\prime}=1\), the accuracy requirements of acceleration for different sampling intervals with an \(\Delta\Delta\overrightarrow{r}(t)\) error at the 1-cm level can be given as \(\sigma_{\frac{-}{r}}=\frac{0.01}{\Delta t^{2}}\), as shown in Table 1. Thus, the errors or magnitudes of the accelerations less than the requirements can be practically ignored. In fact, most of the accelerations of LEOs are below \(10^{-4}\) ms\({}^{-2}\)[34].  
The precise dynamic orbit determination of the GPS and LEO satellites shows that it is usually possible to obtain the satellites' acceleration from known orbit dynamic models, with an accuracy better than \(10^{-5}\) ms\({}^{-2}\)[34; 30; 31]. Therefore, the STP values \(\Delta\Delta\overrightarrow{r}^{G}\) and \(\Delta\Delta\overrightarrow{r}_{r}\) can be computed via Equation (8), and used to remove the geometric variations in GPS data without a high-precision a priori orbit.  
What is more, in the real-time application, given a priori orbit error \(\delta r=3m\) (broadcast ephemeris for GPS satellites [35], a real-time GPS code solution for LEO satellites), the max acceleration error \(\delta a\) comes from the Earth's center of gravity.  
\[\delta a\approx\frac{GM}{r^{3}}\delta r\approx\frac{g}{r}\delta r\approx 5 \times 10^{-6}\;\text{ms}^{-2}<10^{-5}\;\text{ms}^{-2},\]  
where \(GM\) is the Earth's gravitational constant, \(r\) is the distance from the satellite to the Earth's center, and \(g\approx 9.8\) ms\({}^{-2}\).  
As Table 1 suggests, \(\delta a\) can be practically ignored, even for 30-s intervals. Thus, the term \(\Delta\Delta(\overrightarrow{e}\cdot d\overrightarrow{r})\) in Equation (4) can be derived precisely, even in real-time applications.