## 2 Data and Methods
### CIT Methods  
CIT is a classical inverse problem where the ionospheric electron density distribution in the tomographic region can be reconstructed based on the discrete inversion theory using the STEC data on the ray path from the GNSS receivers to the GNSS satellites, in which, the STEC can be expressed by the integral of the electron density on the ray path, as follows:  
\[STEC=\int_{l}Ne(h)dh \tag{1}\]  
Figure 1: Global distribution of ground-based global navigation satellite systems (GNSS) stations. \begin{table}
\begin{tabular}{c c} \hline \hline
**Items** & **Design Parameters** \\ \hline Orbital altitude & 300 km \\ Number of satellites & 15/30/60/96/192 \\ Number of orbital planes & 3/6/12/12/12 \\ Satellites per plane & Varies \\ Orbital inclination & 86.5\({}^{\circ}\) \\ Right ascension of ascending node (RAAN) & Evenly spaced \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Example design parameters of the low earth orbit (LEO) satellite constellation (300 km altitude). where \(STEC\) is the total electron content on the ray path, \(l\) is the ray path from the GNSS receivers to the GNSS satellites, and \(Ne(h)\) is the electron density at the altitude \(h\). In GNSS tomography, the voxel-based model is a commonly used method, which generally divides the tomography region into several discrete voxels, then composes the observation equation based on the distance of each ray through the voxel and the STEC, and finally solves the value of the electron density in each grid using the mathematical method as below:  
\[ \begin{bmatrix}STEC_{1}\\ STEC_{2}\\ \cdots\\ STEC_{m}\end{bmatrix} = \begin{bmatrix}A_{11}&A_{12}&\cdots&A_{1n}\\ A_{21}&A_{22}&\cdots&A_{2n}\\ \cdots&\cdots&\cdots&\cdots\\ A_{m1}&A_{m2}&\cdots&A_{mm}\end{bmatrix} \bullet \begin{bmatrix}x_{1}\\ x_{2}\\ \cdots\\ x_{n}\end{bmatrix} + \begin{bmatrix}e_{1}\\ e_{2}\\ \cdots\\ e_{m}\end{bmatrix}  \tag{2}\]  
where \(m\) denotes \(m\) observed rays, \(n\) denotes the number of voxels in the tomography region, \(\mathrm{Amn}\) denotes the intercept value of the \(m\)th observed ray in the \(n\)th voxel, and \(x_{n}\) and \(e_{n}\) denote the value of the electron density and the error in the \(n\)th voxel, respectively. The most commonly used algorithm in the process of tomography reconstruction of the global ionospheric IED distribution is the algebraic reconstruction method. This algorithm uses an empirical model to obtain the electron density as the initial background field, and iteratively solves sequentially for the electron density in the voxel \(x\). Since the algebraic reconstruction algorithm relies heavily on the initial background field, the ionospheric tomography results are strongly influenced by the initial value, resulting in a loss of accuracy. Therefore, in this paper we use the regularized least squares (RLS) method to perform the laminar inversion reconstruction. The RLS algorithm adds the laminar observational equations to the constraint matrix to solve the ill-posed problem of inversion, and the commonly used regularization method is Tikhonov regularization, as follows:  
\[\mathrm{min}_{x}\|Ax-y\|^{2}+\lambda\|Rx\|^{2} \tag{3}\]  
where \(A\) is the intercept matrix in the ionospheric tomography observation equation, \(x\) is the electron density value to be solved, \(y\) is the STEC on the ray path, \(\lambda\) is the regularization parameter to control the weight of the regularization term, the L-curve method is used to obtain the optimal value of \(\lambda\), and \(R\) is the constraint matrix. In this paper, the constraint matrix consists of three parts, which are the first-order difference matrix (for smoothing the solution), the second-order difference matrix (for smoothing the small curvature solution), and the physical background constraint matrix (for constraining the altitude-direction electron density profile).