## 2 Materials and Methods
### Orbit Estimation with an Unscented Batch Filtering Method
#### 2.2.2 Orbit Estimation Based on an Unscented Batch Filtering Method  
The nonlinear system considered in this study consists of the nonlinear dynamics function \(\mathbf{f}(\cdot)\) and the measurement function \(\mathbf{h}(\cdot)\), formulated as follows:  
\[\left\{ \begin{array}{l}\mathbf{x}_{j+1}=\mathbf{f}\big{(}\mathbf{x}_{j},t_{ j}\big{)}\\ \mathbf{z}_{j}=\mathbf{h}\big{(}\mathbf{x}_{j},t_{j}\big{)}+\mathbf{v}_{j}\end{array} \right. \tag{8}\]  
where \(\mathbf{x}_{j}\) denotes the state vector at time \(t_{j}\) with a covariance \(\mathbf{P}_{j}\), and \(\mathbf{z}_{j}\) represents the measurement vector. The additional measurement noise vector, \(\mathbf{v}_{j}\), is modeled as a zero-mean Gaussian random vector with a covariance of \(\mathbf{R}_{j}\). Notably, no process noise is added to the system model, as the iterative batch filtering scheme employed in this work ensures convergence through successive refinements. To account for the effects of continuous low-thrust maneuvers, the state vector is augmented beyond the traditional position and velocity vectors in the ECI frame to include the unknown thrust acceleration vector expressed in the RTN frame. Specifically, only the along-track component, denoted as \(a_{\text{T}}\), is estimated, in accordance with the model assumptions described in Equation (4). The measurementvector consists of the slant range \(\rho\), radar azimuth \(Az\) and elevation \(El\), \(\mathbf{Z}=[\rho\,Az\,El]\), all of which are functions of the satellite's position relative to a ground tracking station. The theoretical calculations of the observations \(\rho\), \(Az\), and \(El\) at specific epochs, based on the satellite's ECI position vector and the geographic coordinates of the ground station, are integrated into the function \(\mathbf{h}=(\cdot)\). For detailed derivations and formulations of these observation equations, readers are referred to the relevant literature, such as ref. [3]. Given a sequence of observations \(\widetilde{\mathbf{Z}}\) over the time series \(t_{j}\) for \(j=1,\,\ldots,\,N\), where \(N\) denotes the total number of observation epochs, the objective is to estimate the augmented state vector at a selected time \(t_{k}\), \(k\in[1,\ldots,N]\), using both the measurement data and the orbital dynamics model. The rest of this section provides a detailed formulation and implementation of the proposed orbit estimation method based on UT. The initial estimates of the satellite's position and velocity, denoted as \(\mathbf{r}_{\text{initial}}\) and \(\mathbf{v}_{\text{initial}}\), are obtained using the Lambert method under the assumption of two-body Keplerian dynamics, with no thrust maneuvers considered. In the absence of prior information regarding the thrust acceleration, its initial value is set to 0 as a general assumption. Accordingly, the initial augmented state estimate \(\hat{\mathbf{x}}_{k}\) and its covariance estimate \(\hat{\mathbf{P}}_{k}\) at the chosen epoch, \(t_{k}\) is set by  
\[\hat{\mathbf{x}}_{k}=[r_{k,\,\text{initial}}\,\,\mathbf{v}_{k,\,\text{initial}} \,\,0],\,\,\hat{\mathbf{P}}_{k}=\mathbf{P}_{k,\,\text{initial}}, \tag{9}\]  
where \(\mathbf{P}_{k,\,\text{initial}}\) is prespecified. The sigma points of \(\hat{\mathbf{x}}_{k}\) are constructed in the way presented in Section 2.2.1, and are denoted by \(\hat{\mathbf{x}}_{k,k}\), \(i=0,\,\ldots,\,\text{n}\). The sigma set \(\hat{\mathbf{x}}_{i,k}\) is then propagated to all measurement moments \(t_{j}\) (\(j=1,\,\ldots,\,N\), \(j\
eq k\)) using the nonlinear function \(\mathbf{f}\) from Equation (4)  
\[\overline{\mathbf{x}}_{i,j}=\mathbf{f}\big{(}\hat{\mathbf{x}}_{i,k},\,\,\,\,t _{j}-t_{k}\big{)},\,j=1,\,\ldots,\,N\,\,\text{and}\,\,j\
eq k \tag{10}\]  
where \(\overline{\mathbf{x}}_{i,j}\) represents the propagated sigma point.