## 2 MSFGWO Algorithm Model
### Fitness Function DGDOP  
Due to the adoption of the Doppler positioning method in LEO positioning, the basic observation equation for Doppler positioning is given by Equation (8) [30]:  
\[\dot{\rho}=(v_{r}-v^{s})\frac{x^{s}-x_{r}}{\parallel x^{s}-x_{r}\parallel}+c( \delta t_{r}-\dot{\delta t}^{s})+\varepsilon_{\dot{\rho}}, \tag{8}\]  
where \(\dot{\rho}\) represents the pseudo-range rate, which is the Doppler observation in units of m/s; \(\dot{\delta t}^{s}\) and \(\dot{\delta t}_{r}\) represent the temporal derivatives of the satellite clock bias and receiver clock bias, respectively; \(\varepsilon_{\dot{\rho}}\) represents random noise in the Doppler observations; \(x_{r}=\left[x_{r}\quad y_{r}\quad z_{r}\right]^{\text{T}}\) and represents the position and velocity of the receivers, respectively; \(x^{s}=\left[x^{s}\quad y^{s}\quad z^{s}\right]^{\text{T}}\) and \(v^{s}=\left[v^{s}_{x}\quad v^{s}_{y}\quad v^{s}_{z}\right]\) represent the position and velocity of the satellite, respectively.  
For n satellites, this can be written in matrix form, as follows:  
\[\mathbf{G}\cdot\Delta\mathbf{u}=\mathbf{b}, \tag{9}\]  
\[\begin{bmatrix}\frac{\partial\dot{\rho}^{s1}}{\text{d}x}\quad\frac{\partial \dot{\rho}^{s1}}{\partial y}\quad\frac{\partial\dot{\rho}^{s1}}{\partial x}&1 \\ \frac{\partial\dot{\rho}^{s2}}{\partial x}\quad\frac{\partial\dot{\rho}^{s2}}{ \partial y}\quad\frac{\partial\dot{\rho}^{s2}}{\partial x}&1\\ \vdots&\vdots&\vdots&\vdots\\ \frac{\partial\dot{\rho}^{sn}}{\partial x}\quad\frac{\partial\dot{\rho}^{sn}}{ \partial y}\quad\frac{\partial\dot{\rho}^{sn}}{\partial x}&1\\ \end{bmatrix}_{x=x_{k}}\begin{bmatrix}x_{k+1}-x_{k}\\ y_{k+1}-y_{k}\\ z_{k+1}-z_{k}\\ \delta t_{rk+1}-\delta t_{rk}\end{bmatrix}=\begin{bmatrix}\dot{\rho}^{s_{1}}( \mathbf{u})-\dot{\rho}^{s_{1}}(\mathbf{u}_{k})\\ \dot{\rho}^{s_{2}}(\mathbf{u})-\dot{\rho}^{s_{2}}(\mathbf{u}_{k})\\ \vdots\\ \dot{\rho}^{s_{n}}(\mathbf{u})-\dot{\rho}^{s_{n}}(\mathbf{u}_{k})\end{bmatrix}, \tag{10}\]  
In Equations (9) and (10), \(\mathbf{G}\) represents the coefficient matrix of the positioning system; \(\frac{\partial\dot{\rho}^{sn}}{\partial x}\), \(\frac{\partial\dot{\rho}^{sn}}{\partial y}\), and \(\frac{\partial\dot{\rho}^{sn}}{\partial x}\) represent the Doppler observations \(\dot{\rho}^{s_{n}}\) of the receiver with respect to satellite n in terms of partial derivatives with respect to x, y, and z, respectively, at an initial estimate point \(\mathbf{u}_{k}\); \(\mathbf{b}\) is the Doppler shift change amount; \(\dot{\rho}^{s_{n}}(\mathbf{u})\) and \(\dot{\rho}^{s_{n}}(\mathbf{u}_{k})\) represent the actual measured Doppler value and predicted value for satellite n by the receiver, respectively.  
Figure 1: Position update process for grey wolves in the GWO algorithm.  
The matrix \(\mathbf{Q}\) is defined as an amplification factor for Doppler measurement errors.  
\[\mathbf{Q}=\left(\mathbf{G}^{\mathrm{T}}\mathbf{G}\right)^{-1}=\begin{bmatrix}q_{ 11}&q_{12}&q_{13}&q_{14}\\ q_{21}&q_{22}&q_{23}&q_{24}\\ q_{31}&q_{32}&q_{33}&q_{34}\\ q_{41}&q_{42}&q_{43}&q_{44}\end{bmatrix}. \tag{11}\]  
Therefore, the DGDOP of Doppler positioning can be determined accordingly:  
\[\mathrm{DGDOP}=\sqrt{q_{11}+q_{22}+q_{33}+q_{44}}. \tag{12}\]