## 2 Materials and Methods
### Orbit Estimation with an Unscented Batch Filtering Method
#### 2.2.2 Orbit Estimation Based on an Unscented Batch Filtering Method
At the chosen epoch \(t_{k}\), the propagated state, covariance, and sigma sets are directly given by their corresponding estimated values  
\[\overline{\mathbf{x}}_{k}=\hat{\mathbf{x}}_{k},\,\,\overline{\mathbf{P}}_{k} =\hat{\mathbf{P}}_{k},\,\,\overline{\mathbf{x}}_{i,k}=\hat{\mathbf{x}}_{i,k}. \tag{11}\]  
At each observation epoch \(t_{j}\), the theoretical measurements \(\boldsymbol{\gamma}_{i,j}\) corresponding to each propagated sigma point are first computed through the nonlinear function \(\mathbf{h}\), as follows:  
\[\boldsymbol{\gamma}_{i,j}=\mathbf{h}\Big{(}\overline{\mathbf{x}}_{i,j},t_{j} \Big{)},\,\,i=0,\ldots,n;\,\,j=1,\ldots,\,N\,\,\text{and}\,\,j\
eq k. \tag{12}\]  
The predicted measurement vector \(\overline{\mathbf{z}}_{j}\) at each observation epoch \(t_{j}\) is then obtained as a weighted sum of the values of \(\boldsymbol{\gamma}_{i,j}\) for all propagated sigma points at the corresponding epoch:  
\[\overline{\mathbf{z}}_{j}=\sum_{i=0}^{n}\omega_{i}\boldsymbol{\gamma}_{i,j}, \,\,j=1,\ldots,\,N\,\,\text{and}\,\,j\
eq k. \tag{13}\]  
Given that both \(\boldsymbol{\gamma}_{i,j}\) and \(\overline{\mathbf{z}}_{j}\) originate from the state at epoch \(t_{k}\), the sets consisting of \(\boldsymbol{\gamma}_{i,j}\) and \(\overline{\mathbf{z}}_{j}\) for all measurement moments (with the exception of the \(t_{k}\) moment) are denoted by \(\widetilde{\boldsymbol{\gamma}}_{i,k}\) and \(\overline{\mathbf{Z}}_{k}\), respectively,  
\[\widetilde{\boldsymbol{\gamma}}_{i,k}=\left[ \begin{array}{c}\boldsymbol{ \gamma}_{i,1}\\ \boldsymbol{\gamma}_{i,2}\\ \vdots\\ \boldsymbol{\gamma}_{i,N}\end{array} \right],\,\,\overline{\mathbf{Z}}_{k}= \left[ \begin{array}{c}\overline{\mathbf{z}}_{1}\\ \overline{\mathbf{z}}_{2}\\ \vdots\\ \overline{\mathbf{z}}_{N}\end{array} \right]. \tag{14}\]Then, the covariances for the propagated measurements \(\overline{\mathbf{P}}_{k}^{Z}\) can be derived as follows:  
\[\overline{\mathbf{P}}_{k}^{Z}=\sum_{i=0}^{n}\omega_{i}\big{(}\widetilde{\mathbf{ \gamma}}_{i,k}-\overline{\mathbf{Z}}_{k}\big{)}\big{(}\widetilde{\mathbf{\gamma }}_{i,k}-\overline{\mathbf{Z}}_{k}\big{)}^{T}+\widetilde{\mathbf{R}}_{k} \tag{15}\]  
where \(\widetilde{\mathbf{R}}_{k}\) is the measurement noise matrix and the cross-correlation matrix \(\overline{\mathbf{P}}_{k}^{zZ}\) is given as follows:  
\[\overline{\mathbf{P}}_{k}^{zZ}=\sum_{i=0}^{n}\omega_{i}\big{(}\widetilde{ \mathbf{x}}_{i,k}-\overline{\mathbf{x}}_{k}\big{)}\big{(}\widetilde{\mathbf{ \gamma}}_{i,k}-\overline{\mathbf{Z}}_{k}\big{)}^{T}. \tag{16}\]  
Furthermore, the filter gain \(\mathbf{K}_{k}\) is calculated according to the following expression:  
\[\mathbf{K}_{k}=\overline{\mathbf{P}}_{k}^{xZ}\big{(}\overline{\mathbf{P}}_{k }^{Z}\big{)}^{-1}. \tag{17}\]  
The measurement residual is determined as the difference between the actual measurement (\(\widetilde{\mathbf{Z}}\)) and the predicted measurement (\(\overline{\mathbf{Z}}_{k}\)) at all measurement moments  
\[\Delta\widetilde{\mathbf{z}}_{k}=\widetilde{\mathbf{Z}}-\overline{\mathbf{Z}} _{k}. \tag{18}\]  
The satellite augmented state vector at the selected epoch \(t_{\text{k}}\) is updated as indicated below  
\[\mathbf{\hat{x}}_{k}=\overline{\mathbf{x}}_{k}+\mathbf{K}_{k}\Delta\widetilde {\mathbf{z}}_{k}. \tag{19}\]  
An iterative procedure is applied to estimate the satellite's augmented state vector \(\mathbf{\hat{x}}_{k}\).