# Autonomous Orbit Determination of LLO Satellite Using DRO-LLO Links and Lunar Laser Ranging
## 3 Autonomous Orbit Determination Algorithm
### Extended Kalman Filter Algorithm  
Considering the nonlinear situation of practical engineering, this article chooses the Extended Kalman Filter (EKF) to process the measurement data of each epoch. Given the initial epoch state parameter \(\mathbf{X}_{0}\) and state covariance matrix \(P_{0}\) of the satellite, the EKF is used to process the inter-satellite measurement data and lunar laser ranging data, achieving the optimal estimation of the estimated state parameters. The definition of the estimated state is as follows:  
\[\mathbf{X}=\left(\mathbf{Y}^{\text{T}},C_{\text{r}}\right)^{\text{T}}=\left(\mathbf{r}^{ \text{T}},\mathbf{v}^{\text{T}},C_{\text{r}}\right)^{\text{T}} \tag{10}\]  
where \(C_{\text{r}}\) is the solar radiation pressure coefficient. Among the non-conservative forces acting on the satellite, solar radiation pressure is the largest in magnitude. However, the physical mechanism of solar radiation pressure indicates that the radiation pressure models of different spacecrafts have their specificities, and the modeling work is relatively complex [27]. Therefore, the solar radiation pressure coefficient is taken as the estimated parameter. The EKF process includes two steps: time update and measurement update. Time update stage: this uses the state parameter \(\mathbf{X}_{i}\) and state covariance matrix \(\mathbf{P}_{i}\) of epoch \(t_{i}\) to calculate the prior state parameter \(\mathbf{X}_{i+1|i}\) and prior state covariance matrix \(\mathbf{P}_{i+1|i}\). The calculation formula is as follows [7; 22; 28]:  
\[ \begin{split}\mathbf{X}_{i+1|i}&=\mathbf{X}(t_{i+1},\mathbf{X} (t_{i})=\mathbf{X}_{i})+\mathbf{\Gamma}_{i+1|i}\mathbf{u}_{i}\\ \mathbf{P}_{i+1|i}&=\mathbf{\Phi}_{i+1|i}\mathbf{P}_{i}\mathbf{ \Phi}_{i+1|i}^{\text{T}}+\mathbf{\Gamma}_{i+1|i}\mathbf{Q}\mathbf{\Gamma}_{i+1|i}^{\text{T }}\end{split}  \tag{11}\]  
where \(\mathbf{Q}\) is the covariance matrix of process noise, with diagonal elements representing process noise compensation values. The compensation values are set based on the magnitude of unmodeled dynamic errors, and adjusting \(\mathbf{Q}\) can prevent the state covariance matrix from converging too quickly. In the simulation test, the magnitude of the diagonal elements in \(\mathbf{Q}\) are determined by the trial-and-error method [20]. The principle relied on is that the posterior O-C residual conforms to the characteristics of Gaussian noise, and the standard deviation is close to the magnitude of the ranging thermal noise. The state transition matrix \(\mathbf{\Phi}_{i+1|i}\) and the process noise transition matrix \(\mathbf{\Gamma}_{i+1|i}\) are represented as follows:  
\[\mathbf{\Phi}_{i+1|i}= \begin{bmatrix}\mathbf{\Phi}_{\text{Y}}(t_{i+1},t_{i})& \frac{\partial\mathbf{Y}(t_{i+1})}{\partial\mathbf{C}_{\text{r}}(t_{i})}\\ 0_{1\times 6}&1\end{bmatrix}  \tag{12}\]  
\[\mathbf{\Gamma}_{i+1|i}= \begin{bmatrix}\mathbf{\Gamma}_{\text{Y}}(t_{i+1},t_{i})\\ 0_{1\times 3}\end{bmatrix}  \tag{13}\]Measurement update stage: this uses the measurement data \(z_{i+1}\) of the current epoch to calculate the state parameter \(\mathbf{X}_{i+1}\) and state covariance matrix \(\mathbf{P}_{i+1}\), with the following calculation formula:  
\[ \begin{array}{rl}\mathbf{K}_{i+1}=&\mathbf{P}_{i+1|i}\mathbf{H}_{i+1}^{\mathrm{T}}\big{(} \mathbf{R}_{i+1}+\mathbf{H}_{i+1}\mathbf{P}_{i+1|i}\mathbf{H}_{i+1}^{\mathrm{T}}\big{)}^{-1}\\ \mathbf{X}_{i+1}=&\mathbf{X}_{i+1|i}+\mathbf{K}_{i+1}\big{(}z_{i+1}-h\big{(}\mathbf{X}_{i+1|i} \big{)}\big{)}\\ \mathbf{P}_{i+1}=&\Big{(}\mathbf{I}-\mathbf{K}_{i+1}\mathbf{H}_{i+1}\Big{)}\mathbf{P}_{i+1|i}\Big{(} \mathbf{I}-\mathbf{K}_{i+1}\mathbf{H}_{i+1}\Big{)}^{\mathrm{T}}+\mathbf{K}_{i+1}\mathbf{R}_{i+1}\bm {K}_{i+1}^{\mathrm{T}}\end{array}