Multilayer Ionospheric Model Constrained by Physical Prior Based on GNSS Stations

The need for accurate modeling of the ionosphere plays an important role in the global navigation satellite system (GNSS) positioning. The traditional multilayer VTEC model without prior has been used for modeling the ionospheric delay error. However, it is assumed that the electron density of the ionosphere is compressed into multiple thin layers at fixed heights in the lack of capturing ionospheric physics. In this article, the data enhancement method by virtual observations is proposed to build the constrained multilayer VTEC model to capture physical features from empirical ionospheric models. The extraction methods of physical knowledge have been developed by prior VTEC based on the principal component analysis and model coefficients based on the EBF. The constrained multilayer modeling has been verified based on simulation and real measurement of GNSS data in Yunnan, China, collected from Ground-based GNSS stations by Qianxun on November 3, 2021. The receiver DCB error estimated by the multilayer model with prior constraint is significantly lower than that of the single-layer model and the traditional multilayer model. The experimental test shows that the constrained multilayer model achieves the accuracy of $0.5 \,\rm {TECU}$ for the independent reference station. The dSTEC of the proposed two multilayer models are significantly lower than those of the single-layer model for low elevation angles, and the RMSE of dSTEC is reduced by 63$\%$ with the cutoff elevation angle of $10^{\circ }$. The spatial distribution of the multilayer VTEC model shows consistency with the tomography model to verify vertical feature-capturing capability. Compared with the undifferenced and uncombined precise point positioning without ionospheric constraint, the multilayer-constrained model based on the test data improves the convergence time approximately by 36.55% and 18.78% in the horizontal (H) and up (U) directions, respectively. These results demonstrate that the proposed multilayer models not only improve ionospheric delay estimation precision but also can obtain the VTEC distribution capturing the physical characteristics of the ionosphere. The proposed multilayer models may be valuable for the ionospheric delay modeling of satellite navigation systems under harsh variable ionospheric conditions.


I. INTRODUCTION
T HE Earth's ionosphere is the atmosphere within 60 − 2000 km from the ground, and it contains a large number of free electrons. The ionospheric delay is considered one of the main error sources in the process of navigation signal propagation, which can be expressed by the electron density integral on the navigation signal transmission path between satellite and receiver, namely slant total electron content (STEC) [1]. In navigation and positioning applications, dual-frequency and triple-frequency global navigation satellite system (GNSS) users can correct the ionospheric influence using an ionosphere-free combination while single-frequency users usually use an ionospheric model, including empirical model such as Klobuchar and the NeQuick model [2], [3], ionospheric total electron content (TEC) model such as the VTEC modeling [4], ionospheric STEC estimate model such as the location-based linear interpolation model and so on [5], [6].
The traditional ionospheric TEC modeling method is the single-layer vertical TEC (VTEC) model, which assumes the STEC is converted to VTEC by the single-layer mapping function (MF) [4]. The global ionospheric map (GIM) products provided by European Space Agency (ESA), Center for Orbit Determination in Europe (CODE), and Chinese Academy of Sciences (CAS) are the single-layer VTEC model using the spherical harmonic (SH) function [4], [5], [6], [7], [8], [9], [10], [11]. However, the single-layer model does not take into account the characteristics of the ionosphere varying with heights [12]. The simplified MF assumes that the electron density around the ionospheric pierce point (IPP) is symmetrically distributed but the actual distribution of the electron density is not uniform in the ionosphere. The ionospheric multilayer model has been widely studied to capture vertical features of the ionosphere and improve the accuracy of ionospheric delay estimation. The ionospheric multilayer model can be divided into two types according to the reconstruction elements of the model. The first type is the ionospheric tomography method of electron density, and the second one is the multilayer density integration VTEC model. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ The ionospheric tomography method uses observation between satellite and receiver to retrieve three-dimensional electron density by the voxel-based and the basis function methods [13]. The tomographic ionosphere model (TOMION) has been developed at Universitat Politècnica de Catalunya (UPC) since 1995, which is a two-layer voxel-based model [5], [14]. Meanwhile, a global scale three-layer electron density model using the basis function method in 1999 [6] has been constructed by Jet Propulsion Laboratory (JPL), which uses the cubic spline basis function horizontally in longitudes and latitudes and uses the piecewise constant vertically [6]. The data-driven tomography method combined with compressed sensing (CS), called the CS-PCA tomography method, has been proposed to realize quasi-real-time electron density reconstruction and alleviate the ill-posed inversion caused by sparse observation [15].
Distinguishing from tomography, the multilayer VTEC model assumes that electron density in the ionosphere is compressed into multiple thin layers at several fixed heights [16], [17], [18], [19], [20]. The main difference between the multilayer VTEC model and the ionospheric tomography model is the inversion parameters. The inversion parameter of the CS-PCA model is the three-dimensional electron density while the inversion parameter of the multilayer ionospheric model is the two-dimensional VTEC. Compared with the tomography model, the improvement of the multilayer ionospheric model is that the inversion parameters of the multilayer VTEC model are much less than those of the tomography model. The multilayer VTEC model requires less computational cost and relatively low complexity widely adopted in practical GNSS applications. According to implementation strategies for describing the vertical features, the multilayer VTEC modeling consists of two typical scenarios. One is to realize multilayer VTEC model by modifying the MF based on the Chapman layer assumption. The modified MF was proposed by German Aerospace Center (DLR) in 2013, which decrease the modeling error by approximately 50% [16]. Based on the modified MF method, CAS further constructed the global ionospheric multilayer VTEC model [17]. The Barcelona Ionospheric Mapping Function (BIMF) is capable of modifying the MF to reduce the modeling error using the thin-layer assumption [18]. The second type of multilayer VTEC model assumes multiple thin layers vertically rather than a single shell without modifying the MF but simply using the same MF function in each layer. Such a VTEC modeling scenario with multiple thin layers has been operated to generate the GIM products at JPL, and the spline function is selected in the horizontal direction at each layer [19]. Afterward, the ionospheric model based on the two-layer SH function was studied by CAS in 2018 [20]. Compared with the multilayer VTEC model with modified MF, the multilayer VTEC model without modified MF is easier to implement and has lower computing cost. However, there exist challenges in the multilayer VTEC model without modified MF, that is, the VTEC value of each layer tends to be unstable and randomly unrealistic negative due to the lack of vertical information [21], [22].
In view of the negative value or abnormally large values of the VTEC model, progressive solutions have been made to deal with this unrealistic issue [21], [22], [23], [24]. Electronic Navigation Research Institute proposed the residual optimization method to separate the vertical ionospheric delay for each layer [21]. Zhang et al. [22] proposed an inequality-constrained least squares method for global ionospheric single-layer modeling to eliminate the negative value of VTEC in 2012, which shows that the accuracy is consistent with that of the final IGS product without negative VTEC values [22]. Yasyukevich et al. [23] developed a bounded variable least squares fitting algorithm for single-layer modeling to obtain the nonnegative VTEC as well as the differential code biases DCBs. Recently, Maruyama et al. [24] constructed a two-layer VTEC model using a neural network with a nonlinear activation function to realize nonnegative VTEC of each layer. However, most of current methods to eliminate the negative value of VTEC belong to mathematical models without considering actual ionospheric physics.
This article aims to build the multilayer VTEC model to capture vertical physical features and improve the accuracy and stability of the ionospheric modeling process. The virtual observation based on the ionospheric model is incorporated into the observation model by the data augmentation method. We develop two physical prior constrained multilayer models, including VTEC constraints extracted by principal component analysis (PCA) and prior model coefficients extracted by emulated basis function (EBF) method. The simulation experimental test shows that the physical prior constrained multiple layer models improve the stability and accuracy and captures ionospheric features at each layer as the tomography model [15]. The structure of this article is as follows. Section II focuses on describing the constrained multilayer VTEC modeling method. Section III compares the performance of proposed ionospheric multilayer models based on GPS data in Yunnan China from Ground-based GNSS stations by Qianxun. Then, the performance of the proposed method is evaluated quantitatively by the dSTEC analysis and the positioning results of the undifferenced and uncombined precise point positioning (UCPPP) with ionospheric constraints. Finally, Section VI summarizes the conclusion and discussion.

II. METHODOLOGY
This section briefly introduces the method of ionospheric delay extraction. The multilayer ionospheric model with constraints is explained in detail with emphasis on prior construction and regularization methods. The cross-validation methods will be explained.

A. STEC Extraction
In this section, the UCPPP method is used to obtain ionospheric delay data for known fixed station coordinates [25]. The basic GNSS observation models for dual-frequency pseudorange (also known as code-range) P and carrier phase L are defined as follows: where P i and L i pseudorange and phase observations at the ith frequency f i ; ρ true geometric range from the satellite to receiver's antenna phase center in meters; dt k , dt s receiver k and satellite s's clock bias in seconds, respectively; corresponding measurement noises including multipath errors, receiver noise, etc; The "coarse" ionospheric delay I s k including the biases is extracted by UCPPP for fixed known station coordinates: where I s k is the "clean" ionospheric delay without the biases, where the STEC from the satellite s to the receiver k is defined as integral of electron density N e along the path l by STEC s k = s k N e dl.

B. Multilayer Ionospheric Model With Constraints
The single-layer ionospheric model assumes that all free electrons of the ionosphere are located in an infinitely thin layer at a certain ionospheric height h 0 . The STEC can be converted into the VTEC using the ionospheric mapping function MF as where ϕ and λ are the latitude and longitude at ionospheric piercing point (IPP), and θ is the satellite elevation angle at IPP. The single-layer VTEC model describes the horizontal variation of the ionosphere, which can be expressed by the combination of basis functions where B j is the jth horizontal basis function, c j is the corresponding coefficient, and M is the number of basis functions.
The mapping function MF depends on the satellite elevation θ and the ionospheric shell height h 0 commonly written as [4] (9) where R e denotes the mean Earth radius and α is the model coefficient. The value for α is 1 for the single-layer MF (SLM) model, and 0.9782 for the modified SLM model [4].
The multilayer VTEC model can be constructed by assuming that the vertical structure of the ionosphere is composed of multiple thin shells. The satellite-to-receiver slant path intersects each ionospheric shell with the height of h i at the ionospheric pierce point IPP i , where i = 1, 2, . . ., N means the number of shells. The STEC converted from the multilayer VTEC model can be expressed as where MF(h i , θ i ) and VTEC i denotes the MF and the basis function combination of the ith shell. Here, θ i , ϕ i , and λ i are the elevation angle, the latitude, and longitude at IPP i , respectively. Then, one can obtain where B i,j is the jth horizontal basis function of the ith shell with corresponding coefficients c i,j . Here, M i is the number of horizontal basis function of the ith shell. The horizontal basis function B i,j is selected according to the ionospheric characteristics.
In general, the SH function is adopted for global modeling and the polynomial (POLY) function is for regional modeling. The POLY expresses the regional VTEC as a function of geographic latitudes and sun angles at IPPs of this region over a period of time. It is formulated as [4] with S − S 0 = (λ − λ 0 ) + (t − t 0 ) π 12 , ϕ 0 and S 0 are the geographic latitude and solar angle at the regional center, respectively; λ and λ 0 are the geographic longitudes at IPP and regional center point, respectively; (t − t 0 ) is the difference between the observation epoch and the middle epoch during the modeling period. Here, E o,p is the model coefficient with the order o and p, namely the c i,j in (11), and (ϕ − ϕ 0 ) o (S − S 0 ) p is the basis function term, namely B i,j in (11). The maximum model order O and P of the basis function corresponds to latitude and longitude, respectively. Therefore, (O + 1) × (P + 1) is equal to the number of basis function M for each layer.
The virtual observation is proposed to incorporate the prior from the ionospheric model in order to overcome unreasonable vertical density profile. The virtual observation I s, weight k,IPP i of the ith layer is obtained by the following formula: IPP i is the prior VTEC value of the ith layer, and β IPP i is the prior ratio of STEC of the ith layer to the total STEC with the mapping function MF IPP i of the ith layer. The VTEC prior IPP i can be obtained from the prior model by the extraction methods as shown in the next section. (namely H up i ), the fixed IPP height of the thin shell is selected to be h i . The mapping function MF IPP i at each height h i are calculated based on (9) by replacing h 0 . The physical prior VTEC prior IPP i to calculate the prior ratio β IPP i can be extracted from prior models as will be introduced ahead.
We start to derive the physical prior constraint from the original observation equation where y is the K dimensional vector of the original observation from I s k , x is the D dimensional of the unknown parameters, including the model coefficients c i,j , DCB k and DCB s , and ε is the error associated with observation noise and system noise. Here, the matrix A contains the basis function terms B i,j and the coefficients of the corresponding receiver's and satellite's DCB for the observation (the coefficient of satellite's DCB for current observations is 1, otherwise it is 0). In addition, the constraint of the satellites' DCB zero-mean condition is used.
The observation equation corresponding to the virtual measurement in (13) can be further formulated as where y ω is the W -dimensional vector of the virtual observation from I s, weight k,IPP i , and A ω is the W × D virtual design matrix constrained by physical prior. The virtual design matrix A ω includes the basis function terms B i,j for the VTEC model just at ith layer of the ray, and the coefficients of the corresponding receiver's and satellite's DCB for the observation (the coefficient of satellite's DCB for current observation is β IPP i , otherwise it is 0). Here, ε w is the random error of observation constrained by physical knowledge, and ε ω = β IPP i · ε.
The dataset augmentation method is used to realize the ionospheric multilayer VTEC model with physical constraints. The dataset augmentation is a regularization method, which can dramatically reduce the generalization error of a machine learning model [26]. Specifically, virtual observations are constructed by (13) and then are added to the training set. The regularized objective function corresponding to the observation equation by combining physical constraints as virtual observations can be expressed as where λ c represents the relative contribution of the virtual model to real measurements by the norm penalty term. The weight hyperparameter λ c can be adjusted by the grid search method.
The unequally-weighted least square method is used to estimate the unknown vector x in (16).

C. Prior VTEC Construction Method
We design two methods to extract the prior VTEC VTEC prior IPP i to calculate the prior ratio β IPP i . Here, both PCA and EBF will be introduced ahead.

1) Principal Component Analysis (PCA):
The prior VTEC adopted in (13) is obtained by the PCA [27]. The task of PCA is to achieve the dimension reduction of high-dimensional data by linearly projecting onto a lower-dimensional space to make the reconstruction error minimal. Here is a brief description as follows.
1) Design the prior density matrix D N e . The electron density matrix D N e within a certain time period can be extracted from ionospheric physical or statistical models. The n × 1 dimensional electron density vector N e,i is extracted along the vertical direction at the center of the region with a certain temporal resolution, and the electron density vectors of t time steps are collected and described as the n × t matrix D N e = (N e,1 , N e,2 , . . ., N e,t ), N e,i is the n × 1-dimensional vector. 2) Principal Components Extraction ofN e . The principal componentN e can be extracted from D N e by PCA. The covariance matrix Σ is For the n × n covariance matrix Σ, there exists an orthogonal matrix Ψ of the eigenvectors and a diagonal matrix Λ with the eigenvalues, such that The eigenvalues of Λ are sorted from large to small. The percentage of the variance r accounted for the ith eigenvector is expressed as where z i means the ith eigenvalue of Λ. E p is defined as the threshold value for eigenvalues selection. The largest m eigenvalues (m < n) whose sum percentage is greater than E p are selected, namely m i=1 r i ≥ E p . Thus, the main features ofN e are calculated by linearly combining the m eigenvectors Ψ according to their eigenvalue weightsN where ψ i is the ith eigenvector of Ψ. 3) Reconstruction of the prior VTEC. The VTEC prior VTEC prior can be obtained by integrating the principle components of electron densityN e within the height of ionospheric ith layer where are the bottom height and upper height of the ith layer, respectively.

2) Emulated Basis Function (EBF):
The VTEC prior IPP i in (13) can be obtained by the ionospheric emulator, which simulates the prior information of model coefficients for each layer as where  (8), the observation equation corresponding to the prior VTEC model is as follows: where y prior i is the D VTEC,i of the ith layer; x prior i is the prior model coefficients of the ith layer, namely the C i in (22). A prior i is the prior design matrix composed of the model basis function terms B i,j . Then, the least squares estimation is used to estimate the prior model coefficients x i prior basically C i . 3) Construction of prior VTEC: The VTEC prior IPP i can be calculated using (22) with COEF. The first PCA method focuses on obtaining the main features of ionospheric electron density, which can be simply implemented with low computational cost. The second EBF method simulates the prior coefficients of the multilayer model, which represents the spatial (horizontal and vertical) and temporal gradient of the ionosphere.

D. Capacity Hyperparameter Determination Based on Cross Validation
The model's capacity is defined as the ability of the model to fit a wide variety of functions. In this work, neural network parameters such as the number of layers and the order of the polynomial function of the VTEC model are both called the capacity hyperparameters, which both affect the model capacity.
Here, N and M i represent the capacity hyperparameters of the multilayer model. In order to determine the optimal capacity hyperparameters and eliminating statistical uncertainty, the Leave-One-Out Cross-Validation (LOOCV) method and the Bayesian Optimization (BO) method are adopted here.
1) Leave-One-Out Cross-Validation: All stations in the modeling area are divided into three categories of the training set, the validation set, and the test set. On one hand, the stations in the training set are used to minimize the training error and obtain the parameter coefficient of the model c i,j in (11). The stations in the test set, namely independent reference stations, do not participate in the process of modeling, which are only for the performance evaluation of the model. On the other hand, the stations in the validation set are used to determine the capacity hyperparameters N and M i of the model. The LOOCV method is used to repeat training and validation on different datasets divided from the original dataset [29]. First of all, it is assumed that there are N s stations in total, of which R stations are used as the test set R Test . Second, K stations closest to the center of R Test are selected to form the validation set K Validation . Third, the LOOCV method is used to evaluate the performance of the model. Specifically, one station in K Validation is selected as the validation station, and the remaining N s − R − 1 stations are as the training stations. The performance of K validation stations is used to update the capacity hyperparameters to achieve the best performance of the validation set. Finally, the performance of the model is evaluated using the test set R Test .
2) Evaluation Methods of Model Performance: The performance of multilayer models can be evaluated by the differential STEC (dSTEC) analysis method. Also, we introduce criteria for the temporal continuous property of the VTEC modeling.
a) dSTEC analysis: The dSTEC analysis eliminates the biases of the receiver and uses the external products to correct the satellite biases. The dSTEC between any satellite k and the reference satellite r can be defined as follows: where dSTEC s k is the differential STEC of the receiver k and the satellite s. Here, STEC r k is the STEC of receiver k and the reference satellite r, i.e. with the largest elevation angle as Sat r C in Fig. 1.
Actually, the quality assessment is performed by calculating disparity between observation and multilayer VTEC model for the test set ddSTEC = dSTEC vtec − dSTEC obs (25) where dSTEC vtec can be obtained by (10) and dST EC obs is obtained from observations for reference. The root-mean-square error (RMSE) of ddSTEC is adopted to evaluate the accuracy as where d is the number of differential observations for the test site. b) Temporal Gradient: The multilayer VTEC model provides the vertical integral value of electron density in the ionosphere within each layer. Therefore, the VTEC value is in essence nonnegative and is also expected to satisfy temporalspatial continuous characteristics of the ionosphere under quite solar conditions. The time gradient of VTEC as an indicator is defined as where ∇ t VTEC i is the time gradient of VTEC at the ith layer by previous and subsequent continuous epochs. c) Loss function for capacity hyperparameters optimization: BO [28] method is used to determine the optimal capacity hyperparameters of the multilayer VTEC model. The capacity hyperparameter to be optimized includes the model layers N and the number of model basis functions of each layer M 1 , . . ., M i , M N in (11). The loss function for capacity hyperparameter optimization is expressed as where∇ t VTEC is the mean of the temporal gradient ∇ t VTEC i over all epochs. In the loss function, the dSTEC error RMSE (N, M 1 , . . ., M N ) is used to control the accuracy of the model, and the delta VTEC∇ t VTEC i is used to ensure that the optimized model conforms to the ionosphere physics. The BO toolbox is adopted to optimize such hyperparameters [30]. The implementation flowchart of the physics-constrained multilayer model is illustrated in Fig. 2. The procedure of the proposed model includes the following five steps: 1) UCPPP with known fixed station coordinates to extract the ionospheric data; 2) PCA method or EBF method to construct virtual observation; 3) unequally-weighted LS method to calculate the model parameter coefficients; 4) BO to optimize the capacity hyperparameters based on LOOCV; 5) performance evaluation of UCPPP with the ionospheric multilayer VTEC constraints.

III. RESULTS AND ANALYSIS
In order to evaluate the performance of the proposed multilayer VTEC model, the simulation and real observation are adopted. The ionosphere at low latitudes is relatively active; thus, we select Yunnan province in China as the modeling region. The network of 22 GNSS reference stations from the Ground-based GNSS System is depicted in Fig. 3. The 22 stations include test site (red triangle), validation sites (purple diamonds), and training sites (black dots). Among the training sites, the green circle covers five validation sites (purple diamonds) and one test site (independent reference stations) for validation. The disturbance storm time (Dst) index and the Kp index during 2021 are analyzed. The maximum value of the Kp index is greater than 5 and the minimum value of Dst is less than −55 nT on November 3, 2021. The ionospheric model used for constructing the virtual observation is NeQuick2 model, which is more suitable for GNSS practical application [3]. The NeQuick2 is a quick-run ionospheric electron density model particularly designed for trans-ionospheric propagation applications. Table I  the bottom height and upper height of the i layer based on previous research and measured experience [5], [6], [20]. Taking the three-layer model as an example, the ionosphere is empirically divided into three layers vertically, i.e., 80 km − 500 km, 500 km − 800 km, and 800 km − 2000 km. The electrons in three layers are assumed to be concentrated at the fixed altitude, e.g., 300 km, 750 km, and 1000 km [5], [6].
The polynomial function is applied to the regional modeling in Yunnan. It is assumed that the orders of latitude and longitude are the same, so the number of basis functions M i is determined by (1 + O i ) 2 according to (12). Therefore, the number of basis functions M i can be indirectly obtained by optimizing the order of basis function O i in the BO process. The weight hyperparameter λ c in (16) is 0.01 after optimization by the grid search method.
In order to verify the performance of the proposed methods, the single-layer VTEC model and traditional multilayer VTEC model without physical prior are adopted for comparison. The detailed processing strategies of different methods are summarized in Table II. Here, the single-layer model and traditional multilayer model use the least square method to calculate the model parameters while the multilayer model physical constraints use the unequally-weighted least square method. The adopted data are a whole day of GPS data on November 3, 2021 with an epoch interval of 30 s and the cutoff elevation angle of 30 • . In addition, the updating frequency of model coefficients is one group per hour, and the smoothing window is used between each time period.

A. Simulation Test
In order to assess the performance of the proposed method, numerical simulation is adopted. Here, the test site is at the center. First, the STEC for simulated is calculated based on the NeQuick2 model while the STEC for the experimental test is the real observations obtained from the UCPPP algorithm. Second, the DCB of receiver and satellite extracted from the actual observations by the VTEC modeling are added to the simulated STEC to construct the "coarse" ionospheric delay. Third, VTEC models are obtained by the proposed methods and the comparison methods, and the performance of the model is compared by using the dSTEC error, the receiver DCB error, and temporal VTEC evolution.
where DCB real i means the truth value of the ith receiver DCB, and DCB vtec i means the estimated value by the VTEC modeling. The DCB estimation error of each modeling station needs to be subtracted from the corresponding error of the reference station (i.e., i = 1) to obtain the interstation DCB estimation error. The RMSE of DCB errors is 0.3085 TECU, 0.2557 TECU, 0.2265 TECU, and 0.1706 TECU for four different models. The bottom panel shows the temporal distribution of VTEC value over different epochs. The time gradient percentage of VTEC, namely VTEC max grad , is adopted to evaluate the temporal continuity of the model based on (27) The time gradient percentage of VTEC is nearly 0.951%, 54.949%, 0.658%, and 0.666% for four different models in Fig. 4.  , and VTEC max grad for different models based on simulated GNSS data on November 3, 2021. The left, middle, and right panels depict the dSTEC error, the receiver DCB error, and the time gradient of VTEC. Fig. 4 shows that the proposed models have good comprehensive performance. The single-layer model has the largest dSTEC error and DCB estimation error. The traditional twolayer model has lower dSTEC error and DCB estimation error than the single-layer model. However, the VTEC of each layer in the traditional two-layer model has negative and discontinuous values. For the two-layer models constrained by physical prior, the DCB estimation error is the smallest and the VTEC of each layer satisfies the physical property of ionospheric space-time continuity. Although the dSTEC errors of the proposed models are slightly higher than that of the traditional two-layer model, the accuracy can be improved by increasing the model capacity. Fig. 5 shows the simulated results of four models at the independent reference station located in the regional center to verify the performance of the proposed models with increased model capacity hyperparameters . Here, (a) Fig. 6 shows the histogram of three indexes for different models based on simulated GNSS data on November 3, 2021: RMSE ddSTEC , DCB error i , and VTEC max grad . It is noted that the multilayer VTEC models constrained by physical prior can obtain the optimal performance in each index by appropriately increasing the model capacity. First, the three-layer models constrained by physical prior with O 1 = 5, O 2 = 4, and O 3 = 3 have the lowest dSTEC error. Second, the three-layer model constrained by EBF with O 1 = 5, O 2 = 4, and O 3 = 3 has the lowest DCB estimation error. Third, the VTEC gradient value of the traditional multilayer model is abnormally large, which indicates that the model does not satisfy the physical characteristics. However, the multilayer models constrained by physical prior always maintain the physical characteristics of the ionosphere.  Fig. 7 shows the mean of receiver DCB error for different models based on simulated GNSS data on November 3, 2021. Here, the left and right panels depict the mean over different stations and the mean over different hours. It is noted that the single-layer model has the largest error in each station and each time period, and the error of the two-layer model without prior is slightly lower than that of the single-layer model at most casts. The DCB error of multilayer model constrained by physical prior is significantly lower than that of the single-layer and the twolayer model without prior. The three-layer model constrained by EBF with O 1 = 5, O 2 = 4, and O 3 = 3 has the lowest DCB estimation error.

B. Experimental Test
In order to further verify the accuracy and robustness of the proposed constrained multilayer methods in practical applications, the real observation GPS measured data in Yunnan province has been examined for two scenarios: one test site in the center (station number 15) and three test sites at the edge of the region (station number 2, 6, and 19) in Fig. 3 The RMSE of dSTEC are 0.7845 TECU, 0.4907 TECU, and 0.4923 TECU for three different models when the test station is located in the regional center while the RMSE of dSTEC are 1.1897 TECU, 0.9917 TECU, and 0.9949 TECU for these models when the test stations are located at the regional edges. It is noted that the RMSE of dSTEC for the test site at edges of the region are about two times compared with these for the test site in the center of the region. The total VTEC distributions of the single-layer model and the proposed models are similar. The VTEC distributions of each layer of the proposed models conform to the ionospheric physics and have basically the same trend. The dSTEC value fluctuates greatly from UTC time 4:00 to UTC time 12:00, which is positively correlated with the VTEC value. Fig. 9 shows dSTEC over different elevation angles at the independent ground reference station for different models based on actual GNSS data on November 3, 2021. The left panel depicts the RMSE of dSTEC for test site binned by elevation angle, and the right panel depicts the dSTEC distribution over different elevation angles. Here, the cutoff elevation angle of 10 • is adopted in order to analyze the improvement of the proposed models at different elevations. The advantage of the three-layer models constrained by physical prior can be clearly observed for low elevation angle. It is noted that the statistical errors of the proposed two models are the same over different elevation angles, which are significantly lower than those of the single-layer model. The RMSE of dSTEC over different elevation angles of the single-layer model is 1.5406 TECU while that of the three-layer model constrained by PCA is 0.5722 TECU, which is reduced by 63 %. The error reduction of the proposed models at low elevations is better than that at high elevations for the test sites in the center and at the edge of the modeling region. It is due to the fact that the nonuniform variation of electron density along the ray path becomes obvious at low elevation angles with a large obliquity factor. The single-layer model only has the fixed obliquity factor, which cannot express the nonuniform variation, whereas the three-layer models constrained by prior not only have three obliquity factors to present the nonuniform variation but also use the prior physical information to supplement the vertical structure.
2) Comparison With Tomography Results: In this section, the CS-PCA tomography method is used to build the threedimensional electron density model, which is used as a comparison with the multilayer method proposed in this article. The main differences between the multilayer VTEC model and the ionospheric tomography model are the inversion parameters. The similarities between the proposed model and the CS-PCA model include the GNSS observation data collected from Ground-based GNSS stations by Qianxun and ionospheric  features extracted from NeQuick2 model using the similar concept of the data-driven method in two models. and third panels are the VTEC distribution at first, second, and third layers, respectively. The fourth panel is the total VTEC distribution. Here, the first column is the VTEC distribution obtained by integrating the electron density of tomography reconstructed by the CS-PCA method, which is used as the reference value [15]. The VTEC distribution of the tomography model is continuously distributed with longitudes and latitudes at each layer. The total VTEC distribution of the multilayer models constrained by physical prior is consistent with that  of the tomography model. The three-layer models constrained by PCA and the three-layer models constrained by EBF have similar layered VTEC distributions, which shows similar dSTEC errors of the two models. There may have two reasons for the consistency of the two models. First, the two methods are based on the NeQuick2 model to obtain the ionospheric vertical structure for the construction of virtual observations. Second, the weight of virtual observations as the regularization term is low, and the model is still mainly determined by the common actual observations. The VTEC value at the third layer of the tomography model is lower than that of the two models, which may be caused by the difference between the upper boundary height of the tomographic model at 1400 km and the VTEC model at 2000 km.  except that the first layer VTEC and the third layer VTEC of the tomography model fluctuate greatly from UTC time 04:00 to 12:00. The disparity still depends on the VTEC modeling and previous tomography method, which will be further investigated in the future.
3) UCPPP With Ionospheric Constraint: To evaluate the accuracy of the ionospheric model, the positioning accuracy of UCPPP algorithm constrained by the external ionospheric information is tested. The ionospheric delay of the independent reference station is converted from the multilayer VTEC model based on the MF. The constraint on UCPPP with external ionospheric information can be implemented similar as [31] where z + measurement matrix with ionospheric enhancement; z measurement matrix of UCPPP without constraint; H design matrix; y state matrix; v measurement errors; R origin covariance matrix; z 0 ionospheric enhancement information obtained from tomography model; N constraint ionospheric matrix; v 0 error of ionospheric enhancement; R + corresponding augmented covariance matrix; P 0 covariance matrix of ionospheric delay; The external ionospheric information is added to the UCPPP to analyze its influence on the positioning accuracy and convergence time. The parameters of UCPPP are shown in Table III.  Fig. 8. The error of single-layer model-constrained UCPPP represented by blue lines fluctuates more than that of the unconstrained UCPPP in some periods, indicating that the ionospheric prior information with large error cannot improve the positioning accuracy, and lead to the decrease in positioning accuracy. According to the positioning results in Fig. 12 for the 10th reconvergence cycle, the proposed three-layer models with physical constraint require only 6 epochs to reach convergence, the single-layer model requires 30 epochs, and the unconstrained UCPPP takes nearly 50 epochs to reach convergence. In order to accurately evaluate the influence of different ionospheric constraints on positioning, the positioning accuracy and convergence time in the convergence of 21 times per day are counted and averaged.
where error E and error N are the positioning error on the direction of East and North, respectively. The criteria for calculating the convergence time is that the positioning error converges to 10 cm for the first time and lasts for 10 epochs. According to Table IV, the accuracy of positioning parameters of UCPPP constrained by models is better than that of UCPPP without constraint, except that the error in the U direction of UCPPP constrained by the single-layer model is greater than that of UCPPP without constraint. The positioning parameters of the three-layer model with physical constraints are significantly better than those of the single-layer model, and the positioning performance of the three-layer model constrained by PCA is slightly better than that of the multilayer model constrained by EBF. Compared with UCPPP without ionospheric constraints, the multilayer model constrained by PCA improves the convergence time by 36.55% and 18.78% in the H and U directions, respectively.

IV. CONCLUSION AND DISCUSSION
In this article, the constrained multilayer VTEC model is achieved by the data enhancement method to capture vertical physical features of ionosphere and improve the accuracy and stability of the ionospheric model. Two multilayer VTEC models with ionospheric prior are proposed, including the multilayer VTEC model constrained by PCA and multilayer VTEC model constrained by EBF. The implementation process of the models includes the following six steps.
1) UCPPP with known fixed station coordinates to obtain real observation of the ionospheric delay value; 2) PCA method or EBF method to construct virtual observation; 3) the regularized objective function to combine real observations and virtual observations; 4) the unequally-weighted least square method to obtain the model parameter; 5) Cross-validation and Bayesian optimizer to optimize the capacity hyperparameters; 6) performance evaluation by UCPPP with the ionospheric multilayer VTEC constraints. The proposed methods have three advantages over traditional methods. First, the proposed multilayer VTEC models not only use multiple obliquity factors to present the nonuniform variation of electron density along the ray path but also use the prior physical information to supplement the vertical structure. Second, the proposed models combine the prior data from statistical models and real-time observations to improve the accuracy while maintaining the ionospheric physical characteristics. Third, the proposed models can effectively improve the positioning performance of UCPPP when it is added into UCPPP as an ionospheric constraint.
This article compares the performance of the single-layer model, traditional multilayer model without prior, multilayer model constrained by PCA, and multilayer model constrained by EBF based on simulation and real observation GPS data collected in Yunnan region from Ground-based GNSS stations by Qianxun on November 3, 2021. The receiver DCB error estimated by the proposed multilayer models is significantly lower than that of the single-layer model and the traditional multilayer model. The experimental test shows that the constrained multilayer model achieves the accuracy of 0.5 TECU for the independent reference station. The dSTEC of the proposed two multilayer models constrained by prior are significantly lower than those of the single-layer model over different elevation angles with the cutoff elevation angle of 10 • , and the RMSE of dSTEC is reduced by 63%. The spatial distribution of the multilayer VTEC model shows consistency with the tomography model to verify vertical feature-capturing capability. Compared with UCPPP without the ionospheric constraint, the multilayer model constrained by PCA improves the convergence time by 36.55% and 18.78% in the H and U directions, respectively. The multilayer model constrained by PCA and the multilayer model constrained by EBF with the same capacity hyperparameters have similar results. These results show that the proposed multilayer models not only improve ionospheric delay estimation precision but also can obtain the VTEC distribution that captures the physical characteristics of the ionosphere.
The proposed multilayer VTEC models constrained by ionospheric prior can be used in high-precision GNSS technology and are especially important for the detection of the space weather under harsh variable ionospheric conditions. In the future, there are several research directions worth considering, which are as follows: 1) multisystem GNSS data can be used to further improve the temporal and spatial resolution; 2) ionospheric prior information can be constructed from historical data besides the ionospheric model; 3) the form of basis function and the height of each layer can be seen as capacity hyperparameters to realize model optimization based on BO.