## Appendix B Detailed Expressions for the EKF, Matrices, and Coefficient Used in This Paper  
### The Detailed Expressions for the EKF  
The measurement model of EKF is  
\[z(t)=h(x(t),t)+w_{m}(t) \tag{10}\]  
where \(h\) is a nonlinear function of state vector and \(w_{m}\) is the white noise source.  
The discrete form of the measurement model is shown as  
\[z_{k}=H_{k}x_{k}+w_{mk}. \tag{11}\]  
With the measurements, the new optimal state estimate is a linear combination of measurements and previous estimates. Therefore,  
\[\hat{x}_{k}^{+}=L_{k}\hat{x}_{k}^{-}+K_{k}z_{k} \tag{12}\]  
where \(L_{k}\) and \(K_{k}\) are the weight functions to be determined. Substitute Equation (11) into Equation (12):  
\[\hat{x}_{k}^{+}=L_{k}\hat{x}_{k}^{-}+K_{k}H_{k}x_{k}+K_{k}w_{mk}. \tag{13}\]  
Since KF is an unbiased estimation algorithm, the expected value of Equation (13) can be calculated as  
\[L_{k}=I_{k}-K_{k}H_{k}. \tag{14}\]  
Bring Equation (14) into Equation (12):  
\[\hat{x}_{k}^{+}=\hat{x}_{k}^{-}+K_{k}\big{[}z_{k}-H_{k}\hat{x}_{k}^{-}\big{]} \tag{15}\]  
and the error covariance matrix of the measurement can be obtained:  
\[P_{k}^{+}=\left(I-K_{k}H_{k}\right)P_{k}^{-}\left(I-K_{k}H_{k}\right)^{T}+K_{ k}R_{k}K_{k}^{T} \tag{16}\]  
where \(R_{k}\) is the covariance matrix of the measurement noise.  
The best choice of weight function \(K_{k}\) is to minimize the estimation error of \(\hat{x}_{k}^{+}\). In addition, the variance of error estimation is determined by the diagonal element of the error covariance matrix. Therefore, it is necessary to minimize the trace of \(P_{k}^{+}\) matrix with respect to \(K_{k}\), as shown in Equation (10):  
\[\frac{\partial}{\partial K_{k}}\big{[}Tr\big{(}P_{k}^{+}\big{)}\big{]}=0. \tag{17}\]  
Substituting Equation (16) into Equation (17),  
\[-2(I-K_{k}H_{k})P_{k}^{-}H_{k}^{T}+2K_{k}R_{k}=0 \tag{18}\]  
the Kalman gain matrix can be obtained:  
\[K_{k}=P_{k}^{-}H_{k}^{T}\big{(}H_{k}P_{k}^{-}H_{k}^{T}+R_{k}\big{)}^{-1}. \tag{19}\]Update the state vector with the actual measurement vector:  
\[\begin{split}&\hat{x}_{k}^{+}=\hat{x}_{k}^{-}+K_{k}\Big{[}z_{k}-h \Big{(}\hat{x}_{k}^{-},t\Big{)}\Big{]}\\ &=\hat{x}_{k}^{-}+K_{k}r_{k}\end{split} \tag{13}\]  
where \(r_{k}\) is the measurement innovation vector:  
\[r_{k}=z_{k}-h\Big{(}\hat{x}_{k}^{-},t\Big{)}. \tag{14}\]  
Once the state estimate converges to its corresponding true value, the measurement innovation can be modeled as a linear function of the state vector, but direct measurement cannot be modeled as a linear function. Therefore,  
\[r_{k}\approx H_{k}\delta x_{k}^{-}+w_{mk} \tag{15}\]  
where  
\[H_{k}=\left.\frac{\partial h(x,t_{k})}{\partial x}\right|_{x=\hat{x}_{k}^{-}}. \tag{16}\]  
For the tightly coupled integrated navigation in the ECEF coordinate system, \(H_{k}\) can be expressed as  
\[H_{k}=\left(\begin{array}{cccccc}0_{1,3}&0_{1,3}&{h_{\rho}^{1}}^{T}&0_{1,3} &0_{1,3}&1&0\\ 0_{1,3}&0_{1,3}&{h_{\rho}^{2}}^{T}&0_{1,3}&0_{1,3}&1&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\\ 0_{1,3}&0_{1,3}&{h_{\rho}^{m}}^{T}&0_{1,3}&0_{1,3}&1&0\\ 0_{1,3}&{u_{a1}^{e}}^{T}&0_{1,3}&0_{1,3}&0_{1,3}&0&1\\ 0_{1,3}&{u_{a2}^{e}}^{T}&0_{1,3}&0_{1,3}&0_{1,3}&0&1\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\\ 0_{1,3}&{u_{am}^{e}}^{T}&0_{1,3}&0_{1,3}&0_{1,3}&0&1\end{array}\right) \tag{17}\]  
where \({u_{aj}^{e}}^{T}\) is the line of sight unit vector, \({h_{\rho}^{i}}^{T}\) is shown as Equation (18):  
\[{h_{\rho}^{i}}^{T}=\left(\begin{array}{c}\Big{(}R_{N}\Big{(}L_{b}\Big{)}+ \hat{h}_{b}\Big{)}{u_{aj,N}^{e}}\\ \Big{(}R_{E}\Big{(}L_{b}\Big{)}+\hat{h}_{b}\Big{)}\cos\hat{L}_{b}{u_{aj,E}^{e }}\\ -{u_{aj,D}^{e}}^{e}\end{array}\right) \tag{18}\]  
where \(R_{N}\) is the meridian radius of curvature, \(R_{E}\) is the transverse radius of curvature, \(\hat{L}_{b}\) is the geodetic latitude of users, and \(\hat{h}_{b}\) is the geodetic height of users.