Artiﬁcial Neural Network-Based Ionospheric Delay Correction Method for Satellite-Based Augmentation Systems

: Ionospheric delay is a critical error source in Global Navigation Satellite Systems (GNSSs) and a principal aspect of Satellite Based Augmentation System (SBAS) corrections. Grid Ionospheric Vertical Delays (GIVDs) are derived from the delays on Ionosphere Pierce Points (IPPs), which are observed by SBAS reference stations. SBAS master stations calculate ionospheric delay corrections by several methods, such as planar ﬁt or Kriging. However, when there are not enough IPPs around an Ionosphere Grid Point (IGP) or the IPPs are unevenly distributed, the ﬁtting accuracy of planar ﬁt or Kriging is unsatisfactory. Moreover, the integrity bounds of Grid Ionospheric Vertical Errors (GIVEs) are overly conservative. Since Artiﬁcial Neural Networks (ANNs) are widely used in ionospheric research due to their self-adaptation, parallelism, non-linearity, robustness, and learnability, the ANN method for GIVD and GIVE derivation is proposed in this article. Networks are separately trained for IGPs, and ﬁve years of historical data are applied on network training. Principal Component Analysis (PCA) is applied for dimensionality reduction of geomagnetic and solar indices, which is employed as a network input feature. Furthermore, the GIVE algorithm of the ANN method is derived based on the distribution of the residual random variable. Finally, experiments are conducted on 12 IGPs over the East China region. Under normal ionospheric conditions, compared with the planar ﬁt and Kriging methods, the residual reduction of the ANN method is approximately 15%. The ANN method ﬁts the ionospheric delay residual error better. The percentage of GIVE availability under 2.7 m increases at least 25 points in comparison to Kriging. Under disturbed conditions, due to a lack of training samples, the ANN method is incompetent compared with planar ﬁt or Kriging.


Introduction
Global Navigation Satellite System (GNSS)-based positioning services have been extensively used in Safety-of-Life (SoL) aeronautical applications, containing aircraft en route navigation, precision approach, and automatic landing. Nevertheless, the positioning performance of fundamental GNSS services is restricted by various error sources: satellite orbit error, satellite clock error, and ionospheric delay correction error, which make it unable to meet the requirements placed on accuracy, integrity, availability, and continuity of GNSSs for civil aviation users [1,2]. To this end, many countries and regions have established their Satellite-Based Augmentation System (SBAS) to provide GNSS users with additional differential correction information: satellite orbit correction, satellite clock correction, and ionospheric delay correction, and with integrity augmentation information: User Differential Range Error (UDRE) and Grid Ionospheric Vertical Error (GIVE) [3]. The Wide-Area Augmentation System (WAAS) has been fully operational since July 2003, providing SoL services for North and part of South America [4]. The Japanese Multifunctional Satellite Augmentation System (MSAS), European Geostationary Navigation Overlay Service (EGNOS), and Indian GPS-Aided GEO Augmented Navigation (GAGAN) were declared operational and certified for SoL services in 2007 [5], 2011 [6], and 2014 [5], nonlinear adaptive processing systems. Theoretically, a fully connected ANN with at least three layers can approach any continuous function as long as it is given sufficient hidden nodes. However, training small datasets on networks with excessive nodes and layers may cause overfitting problems. Since geomagnetic and solar indices are highly correlated, which will enhance the number of parameters and computational complexity, the PCA technique is applied. PCA is a multivariate statistical analysis method for data dimensionality reduction. It uses an orthogonal transformation to linearly project the observed data of a series of related variables into linearly uncorrelated variables, which are called Principal Components (PCs). One or several PCs that contribute the most to the data variation are sufficient to retain the information while others are abandoned.
The rest of the article is organized as follows. The next section introduces the experiment and data preparation. Then, planar fit, Kriging, and ANN methods for GIVD and GIVE derivation are elaborated in detail. Results of the ANN method under normal and disturbed ionospheric conditions are presented and compared with planar fit and Kriging. After that, the discussion is given. The last section summarizes the entire work and puts forward the prospect for further research.

Experimental Data
The BDSBAS service is an indispensable part of the BeiDou Navigation Satellite System (BDS); it is under construction and will provide single frequency service through a B1C signal for users in China and surrounding areas. Since the reference stations and experimental data are limited, twelve IGPs located in the East China region with dense airlines are selected to test the performance of the proposed method. The ANN method requires reference stations to provide observation data of several years for training, however, there are few stations meeting the requirements in the China region. Therefore, the experimental stations and IGPs in this article are applied to verify the performance of the proposed method. Furthermore, the applicability of the ANN method has no limit in terms of application regions, and further study on the large regional span of China will be replicated when reference stations are available.
Center for Orbit Determination in Europe (CODE) final GIMs are picked as the reference IGP vertical delay, which was 11 days' delay calculated with a Spherical Harmonic (SH) model. It is considered precise and used as a reference for many studies with 1-hour temporal resolution [13,30,32].

Ipp Vertical Delay Extraction
The original Receiver Independent Exchange (RINEX) format data of seven reference stations are open access from the CDDIS website. To extract vertical delays of IPPs and their latitudes and longitudes, GPS observation files (O files) and navigation files (N files) are indispensable. Meanwhile, Differential Code Biases (DCBs) of satellites and stations are also needed; monthly satellite C1C(C1)-C1W(P1) and C1W(P1)-C2W(P2) DCBs are obtained from CODE; daily receiver C1-P2 DCBs are averaged from GNSS analysis centers' Ionosphere Map Exchange (IONEX) files. Reference stations' O files contain 30 s GPS observations of pseudocode C1, P2, and carrier phase L1C(L1), L2W(L2). For a fixed receiver, the observations are expressed as: (1) where R i is the true distance between the receiver and the i-th satellite antenna phase center; I i and T i are ionosphere and troposphere delays; c is the speed of light in vacuum; τ i and τ are satellite and receiver clock errors; D i j and B i j represent the phase delays and advances caused by instrument biases of the i-th satellite at j-th frequency; D j represents the receiver hardware delay at j-th frequency; λ and N i are wavelength and ambiguity of the carrier phase; E i terms are the residuals in the GPS measurements.
Taking the difference of dual-frequency code observations and carrier phase observations separately and neglecting the residuals, the results are shown as: where DCB i = D i 1 − D i 2 and DCB = D 1 − D 2 represent the satellite and receiver C1-P2 DCBs.
The pseudocode measurements are noisy but not ambiguous, while the carrier phase measurements are accurate but have cycle ambiguities. The Hatch filter is adopted to calculate the smoothed code measurements [33], and filtering process could be expressed as: where k is the epoch index and K = 12 (360 s) is the smoothing window width. Before applying the Hatch filter, cycle-slips must be detected and removed from the data, or the cycle ambiguity terms N i in ∆ i L will not cancel each other out between epochs. For a specific frequency, the ionospheric delay of electromagnetic waves could be expressed as a function of slant TEC: Replacing ∆ i C with the filtered ∆ i C,Hatch , and expressing the ionospheric delay terms as the above equations, we end up with the calculation equation of STEC: RINEX N files help to calculate IPPs' longitudes and latitudes at each epoch; meanwhile, the elevation angles of each satellite are known. Since the thin shell model is inaccurate for low elevation angle TEC estimation, IPP elevation angles less than 15°are abandoned. STECs of qualified IPPs are transferred to VTEC by multiplying an obliquity factor F: where Re = 6378.14 km is the radius of the Earth, H = 350 km is the height of the ionosphere thin shell, and A i is the elevation angle which has a cutoff value of 15°. Thus far, the geographic coordinates and vertical delays of IPPs at times corresponding to reference CODE GIMs are extracted from IGS RINEX files.

Methodology
In this section, three SBAS GIVD methods and their corresponding GIVE methods are interpreted in detail. The eligible IPPs are determined according to the IPP search algorithm for Kriging proposed by Sparks [20]. The target number of IPPs is 30, within a maximum fitting radius of 2100 km.

Planar Fit Method
All IPP locations participating in calculation need to be transferred from the Geodetic Coordinate System to the Earth-Centered, Earth-Fixed (ECEF) coordinate system. The transformation equation is given by: where φ i and λ i are, respectively, the latitude and longitude of IPPs. The planar fit method is the most widely employed SBAS ionospheric delay method as a result of its outstanding performance and low computational cost. The ionosphere TEC dispersion could be described properly by a planar trend under normal conditions [19,20]. The estimated ionospheric delay plane (IGP as origin) contains three unknown coefficients and could be expressed as: where ∆X i = X i − X IGP is the distance separating the IPP from the IGP;ê andn are, respectively, the unit vectors in east and north directions of the local Cartesian coordinate system whose origin is at the IGP, specified as follows: where φ IGP and λ IGP are the latitude and longitude of the IGP. The planar model coefficients a 0 a 1 a 2 are estimated by applying the least square fitting on the vertical delays at IPPs. Before applying the least square fitting, the geometric matrix G and the covariance matrix of the observations W should be determined as follows: where σ i denotes the variance of the ionosphere delay measurements at IPPs; σ i,j denotes the covariance between two IPP measurements caused by sharing the same satellite or receiver; σ deco = 0.35 m is the delay covariance associated with widely separated (uncorrelated) IPPs. I = I 1 I 2 · · · I N T represents the vertical delay vector of IPPs and, finally, the IGP vertical delay is estimated using the least square fitting as:

Kriging Method
The Kriging-based ionospheric method was introduced in the WAAS Follow-On Release 3 to replace the planar fit [20]. Kriging is a minimum mean square estimator which estimates the IGP delay by a weighted linear combination of IPP delays. The general equation of GIVD calculation applying Kriging is shown as: where w = w 1 w 2 · · · w N T is the weight vector for IPPs. According to the methodology elucidated by Sparks, the weight vector is solved as a constrained optimization problem by the method of Lagrange multipliers. The solution of w could be expressed as: where s = 1 ∆X T ·ê ∆X T ·n is the distance indication vector, s = 1 0 0 for IGP; W K = [M + C] −1 is similar to the W term in the planar fit method, moreover, Kriging describes the ionosphere decorrelation with a structure called the variogram.
The variogram results C and c are, respectively, an N × N matrix and N × 1 vector, depicting the rate at which the values of a random function evaluated at two places become decorrelated as a function of their distance. M is the measurement covariance matrix without considering decorrelation, i.e., M = W −1 − σ 2 deco · I N×N , where I N×N represents the identity matrix.
For ionosphere delay modeling, the observed behavior should be matched as a linear relation on separations close to the origin and a finite limit at infinity. Therefore, an exponential model is used for the variogram. The elements C i,j in C and c i in c are calculated by: where σ decototal = 0.65 m, σ deconomi = 0.3 m, and d deco = 8000 km are, respectively, the delay covariance associated with uncorrelated IPPs, the delay covariance associated with nearly coincident IPPs, and the characteristic decorrelation distance; D i,j = X i − X j denotes the distance between the i-th IPP and j-th IPP, and d i = |X i − X IGP | is the distance from the i-th IPP to the IGP. In addition, when σ decototal and σ deconomi are set to be equal, the Kriging method downgrades to the planar fit.

GIVE Calculation for Planar Fit and Kriging
The Grid Ionospheric Vertical Error (GIVE) is critical information in the SBAS that ensures the integrity and availability of the system. GIVEs provide 99.9% confidence bounds on GIVD residuals [14]. As mentioned above, the planar fit method is a special case of the Kriging method. Thus, GIVE derivations for the planar fit and Kriging methods are similar.
GIVEs protect from the adverse influence of delay estimation, both sampled and undersampled [16]. The GIVE for the Kriging method is given by the following: where K GIVE = 3.291 is the constant that defines the 99.9% confidence interval, σ 2 IGP corresponds to the variance of measurement error, modeling error, and sampling error of a relevant IGP. The composition of σ 2 IGP is described by the following equation: where R 2 irreg is an inflation factor for ionospheric and statistical uncertainty in chi-square statistics associated with the vertical delay estimation. Any calculated value of R 2 irreg less than one is set to one. The inflation factor is computed by [14]: χ 2 Pmd (DoF) corresponds to the chi-square statistic value with degree of freedom DoF in probability Pmd; the degree of freedom DoF equals to the number of IPPs minus the number of parameters of the model; Pmd is the probability of missed detection and is usually set to 10 −3 .
The accuracy of the vertical delay estimation will diminish as the ionospheric disturbance level increases. χ 2 is introduced to inspect the goodness of fit, i.e., whether an estimation satisfies the accuracy requirement [20]. The detection strategy for ionospheric irregularities based on the chi-square is computed as the sum of the squares of the normalized residuals, written as follows: The SBAS also protects against estimation errors due to the undersampled irregularities of ionosphere TEC except for sampling errors. σ undersampled is known as the undersampled threat model. σ undersampled is both a function of the Relative Centroid Metric (RCM) and the fitting radius [15,34]. The RCM represents the ratio of the distance between the IPP centroid and the IGP to the fitting radius. The undersampled threat model applied in this paper is based on the WAAS release 8/9 model for planar fit and Kriging, and the model is shown in Figure 3 [16,34].
By replacing W K with W, and setting σ decototal = σ deconomi = σ deco , we could calculate GIVEs for the planar fit model following the same process.

ANN GIVD Method
Artificial Neural Networks (ANNs), also known as Multi-Layer Perceptrons (MLPs) or Back Propagation Networks (BPNs), are multilayer forward structured networks that consist of fully connected neurons. There is no necessity to determine the mathematical equation of the mapping relationship between inputs and outputs in advance. Only through the training process can ANNs learn the mapping rules and obtain the closest prediction to the expected result when new inputs are given. SBAS ionospheric delay modeling is a complex problem with abundant variations. Historical data assist ANNs in capturing such variations under different times, geographic locations, geomagnetic conditions, and solar conditions. Next, the ANN method will be introduced in sections. Network structure and parameters, as well as its inputs and outputs, will be elaborated first. Then, the Principal Component Analysis (PCA) method for geomagnetic and solar index processing is presented. Finally, the network training strategy is described in detail.

Network Structure and Parameters
The finalized network is formatted to be a three-layer structure. It contains an input layer, a single hidden layer, and an output layer. There are 93 neurons in the input layer, which matches the number of inputs in each training and testing data point. There are 31 hidden neurons in the hidden layer, which is determined by considering the ratio of the training data size and the parameter size. The output layer has only one neuron, which outputs the vertical delay at the relevant IGP.
The ANN implemented in this article is fully connected. That is to say, there are connection weights between each neuron in the previous layer and all neurons in the next layer. Accordingly, the 93 × 31 weight matrix w 1 and 31 × 1 weight vector w 2 are formatted. In addition, after calculating the weighted sum at each neuron, a bias term is added to the result. Biases in the hidden layer are arranged in a 1 × 31 vector b 1 , while the bias in the output layer is represented as b 2 . w 1 , w 2 , b 1 , and b 2 are the network parameters to be learned through the training process.
The 1 × 93 row vector U j denotes the j-th input vector of the network, V j is the corresponding output. An explicit introduction of typical three-layer neural network can be found in [25], and the network transfer function between inputs and outputs could be expressed as follows: Activation functions introduce nonlinear factors to the system and further enhance the expression ability of linear models. The sigmoid function is selected as the hidden layer activation function. As a pure nonlinear function, despite a risk of gradient disappearance and slower converging speed, sigmoid is convenient to take derivatives on GIVE calculations.

Network Inputs and Outputs
Three scalars and three vectors are stacked horizontally in the 93-dimensional input vector U j . The formation of U j could be expressed in the following form: where α j is the modified day of year; β j is the modified time of day; γ j denotes the index of geomagnetic and solar conditions processed by PCA. Φ j , Λ j , and I j are 30-dimensional row vectors, which represent the relative latitudes, the relative longitudes, and the ionospheric delays of 30 IPPs.
Considering that the sigmoid function has an apparent nonlinear interval between −5 and 5, all input features are scaled to this interval. Since the location of the subsolar point in the tropic affects the local ionosphere activity, day of year is chosen as an input feature. Similarly, the time of day also contains the subsolar point position information, and thus is also set as an input feature. However, directly using such day and time values is meaningless. They need to be modified to match the subsolar point locations. The modification function is shown in Figure 4.  Components of Φ j , Λ j , and I j are not required to be modified, they are scaled to the −5 to 5 interval by applying a min-max normalization method. The method could be expressed as the following equations: where φ, λ, and I denote the components of Φ j , Λ j , and I j ; the parameters with subscript min and max represent the min and max values in the entire dataset. Since φ and λ represent latitude and longitude differences between IPP and IGP, they could be negative, and absolute values are used in the denominators. The details concerning γ j are introduced in the next subsection. The 1-dimensional network output V j is the predicted ionospheric vertical delay at the relevant IGP, in other words, the predicted GIVD.

Principal Component Analysis
As mentioned before, geomagnetic and solar conditions also affect ionosphere TEC dramatically. Hence, geomagnetic indices Kp, Ap, and Dst, solar activity indices Sunspot Number (SSN), and solar radio flux 10.7 (F10.7) are taken into account when applying the ANN method. However, these five indices are highly correlated; when entirely set as network inputs, they will enhance the number of parameters and computational complexity. Thus, the PCA technique is applied to reduce the dimension of the indices from five to one and without missing major information.
γ org denotes the original geomagnetic and solar index matrix covering the total of 30,130 training and validation data points. The indices are Kp, Ap, Dst, SSN, and F10.7 collected from the OMNI data (spacecraft-interspersed, near-Earth solar wind data). The composition of γ org is shown as follows: where K p, Ap, etc. are column vectors of the indices covering the entire training and validation data used for PCA fitting.
Data normalization is needed before calculating principal components, which means subtracting the column average from each element and dividing by the standard deviation. Taking the Kp column as an example, the normalization process could be expressed as the equations: where Kp nor,j and Kp j are j-th components of normalized vector K p nor and the original K p.
The other columns are operated similarly, and γ nor represents the indices matrix after all columns are normalized.
Since the data have been normalized, the 5 × 5 sample correlation matrix R nor could be easily computed by the following equation: the eigenvalue λ k reflects the contribution rate of the corresponding principal component to the integrity. By solving the characteristic equation of the sample correlation matrix R nor , five eigenvalues are settled. Their corresponding eigenvectors b k are determined by solving each eigenvector equation. The calculation equations are shown in the following form: where || is to calculate the determinant of matrix; I 5×5 denotes the identity matrix. Reforming λ k in descending order and calculating the contribution rate of each eigenvalue as a percentage, we obtain 57.00%, 29.94%, 8.56%, 2.90%, and 1.60%. In this article, the first principal component that keeps 57% of the information is selected as an input feature of the ANN (selecting more is also feasible). λ 1 and b 1 represent the first eigenvalue and its corresponding eigenvector, feature γ j is the j-th component of feature vector γ. We obtain γ from the following equation: where the denominator means calculating the maximum absolute value in vector γ nor b 1 . Testing or real-time data are also transformed by applying the fitted PCA parameters, i.e., mean, variance, and eigenvector. By applying PCA, the input dimensions are reduced substantially, while the primary information is retained.

Training Strategy
In the previous subsections, each element of input vector U j is elaborated. The entire dataset U with dimension 33,601 ×93 is separated into a 30,130 ×93 training and validation matrix U train/valid , a 3427 × 93 testing matrix under normal conditions U ntest , and a 44 × 93 testing matrix under disturbed conditions U dtest . The splitting specific of each partition is introduced in Section 2.1. Moreover, considering the data volume and splitting strategy adopted by other works [25,32], U train/valid is randomly split into the training data U train and the validation data U valid in the ratio 8 : 1. Y train , Y valid , Y ntest , and Y dtest are the corresponding reference IGP vertical delay vectors.
Training data are the data that contribute and determine the network parameters through the training process; validation data are used to evaluate the training performance and to avoid overfitting problems, while test data are the simulation of real-time data that test the fitting effect of the ANN model. Before training, the loss function needs to be defined. Loss function pertains to the error to be minimized. In this article, the Mean Absolute Error (MAE) is set to be the loss function and calculated as follows: where V j , which is mentioned in Equation (21), represents the j-th output of the network, and Y j is the corresponding reference. batch represents the batch size in each iteration. As long as GPU memory is available, batch size could be set to the dimension of the entire training data, i.e., batch = 26,703. The training strategy incorporates a feed forward propagation process that determines the loss, an error backpropagation process that calculates the gradients and updates the parameters, and iteration of the above. In order to minimize the loss, each component in the network parameter matrices and vectors w 1 , w 2 , b 1 ,and b 2 has to be modified towards its negative gradient. The gradient is calculated by taking the partial derivative of the loss function with respect to the relevant parameter.
The Adam optimizer is applied to update parameters and minimize loss [35]. Adam is derived from the adaptive moment estimation, it is a variant of the gradient descent algorithm. The learning rate of each parameter is dynamically adjusted by using the first and second order moment estimation of the gradient, which keeps the parameters stable. The parameter update rule of the Adam optimizer could be expressed as: where t is the epoch index; u t andû t+1 are, respectively, biased and corrected first moments; v t andv t+1 are biased and corrected second raw moments; g t denotes the gradient at the t th epoch; lr = 0.001 is the learning rate; β 1 = 0.9 and β 2 = 0.999 are exponential decay rates; = 10 −8 is constant to avoid a zero denominator. Validation data U valid and their corresponding reference Y valid are used to evaluate the training performance, they detect overfitting and help to decide the training epoch number. Epoch is the number of iterations of parameter updates before training stops, it is usually set on when the validation loss starts to rise.
When the training process finishes, all epochs and parameter matrices w 1 , w 2 , b 1 , and b 2 (containing 2946 parameters) are determined and recorded. With U ntest , U dtest , or real-time data passing through, the network outputs the predicted IGP vertical delay results, i.e., GIVD.

GIVE Calculation for ANN Method
When the training process of an ANN is completed, the nonlinear model will have a fixed relationship between its inputs and output. We use a generalized function f (U j , P) to represent this relationship and give the following equations: where P is the set that contains all network parameters. If P * represents the true parameter set for the function which models the system, there only exists a Gaussian random error between network prediction and the reference: where each j is an independently and identically distributed Gaussian random variable with zero mean and variance σ 2 . We minimize the loss function L and derive P as a good estimation of P * by applying the error backpropagation algorithm in the training process. Hence, P is very close to the true values P * , and a Taylor expansion to the first order could be used to approximate f (U j , P) in terms of f (U j , P * ): The prediction residual is defined by Y j − V j , and given Equations (33)- (35), the ionospheric delay residuals could be expressed as: To derive the integrity bound, the expected value, as well as the variance of the residual random variable, needs to be calculated by applying the following equations: Since there is statistical independence between random variables j and P, they are decoupled in calculating expectation and variance. Rivals and Personnaz discussed the distribution of P − P * in detail for neural networks. P − P * is also a Gaussian random variable with zero mean and variance σ 2 (∆ T ∆) −1 [36]. ∆ denotes the Jacobian matrix of the training data (and/or validation data) with respect to the parameters, and has the following form: where n and p are data point number of the training (and/or validation) set and parameter number of the network. For the well-trained network with a large amount of data, each parameter in P is within a small neighborhood of the corresponding element in P * . Therefore, we could replace P * with P in partial derivative calculations. Moreover, P is also regarded as an approximation of P * in calculating the unbiased estimation of σ 2 : If the validation data Y valid and U valid are used to estimate σ 2 , the parameter DoF p does not need to be subtracted from n in the denominator.
To sum up, the prediction residual of the ANN is also a Gaussian random variable with zero mean and variance σ 2 + σ 2 δ j T (∆ T ∆) −1 δ j . GIVE for the neural network GIVD method is derived from the following equation: where t α/2 n−p is the t-distribution inflation factor that gives the (1 − α) × 100% confidence bound for the predicted GIVD. t α/2 n−p = t 0.0005 ∞ = 3.291 for a 99.9% integrity level. Since the planar fit and Kriging methods apply empirical models, the GIVE derivation has to independently consider various threat factors for safety-critical purposes. In contrast, the ANN method derives GIVEs based on modeling of the residual random variable and statistical analysis of historical data. Therefore, the error sources and threat factors are integrally considered in the GIVE method. Specifically, the measurement errors of station receivers are contained in the σ 2 term, and the data points that have large decorrelation or undersampled threat will provide large gradients to the network parameter vector and are further reflected by the δ j T (∆ T ∆) −1 δ j term in Equation (40).

Results
Experiments are conducted on an Intel i5-10400F CPU and Nvidia 1050Ti GPU. The software environment is based on Python 3.7.4 and TensorFlow 2.4.1. The training process is operated on 26,703 pairs of training data U train , Y train and validated on 3427 pairs of validation data U valid , Y valid , which are approximately in a ratio of 8:1. The training process lasts 72 s. Different from the planar fit and Kriging methods, neural network models are trained for each IGP shown in Figure 1. Figure 5 shows the decay of validation loss along with training epochs. It is noticed from Figure 5 that loss functions concerning epochs are not smooth, and light oscillations occur between epochs. After 10,000 epochs, the final validation losses (possibly not the best in history) are between 15 and 20 cm.

Under Normal Ionospheric Conditions
Test data under normal ionospheric conditions are selected from the year 2021, and contain a total of 3427 data points from January to May. Details are introduced in Section 2.1. Before presenting the derivation results, we firstly define the residual evaluating indicators MAE, RMSE, and R by giving the following equations: where V j represents the predicted GIVD for ANNs or fitted GIVD for other methods; Y j is the corresponding reference which is acquired from CODE GIM final products.  To intuitively reveal the fitting accuracy and residual distribution of each existing and proposed method, the following Figure 6 shows the predicted GIVDs by applying various methods as a function of their corresponding references. The residual (reference minus the predicted) distribution histograms are generated. The evaluating indicator R for each IGP and each method is also calculated and presented in Figure 7.  It is known that the fitting accuracy is affected by the accuracy of IPP vertical delay extraction. The Hatch filter-smoothed code measurements are adopted in this connection.
Moreover, the distribution of IPPs relative to the IGP is also a dominant factor that affects the predicted accuracy. Once locations of the reference stations and constellations used are determined, the typical distributions of the IPPs are hardly changed. Figure 8 presents the typical distributions of IPPs around each IGP. GIVE provides the confidence bounds for GIVD residuals and is an essential index of SBAS integrity. To a certain extent, the reliability of GIVE estimation is more significant than the accuracy of GIVD estimation in the SBAS. Figures 9-11 present the absolute values of the GIVD residuals with corresponding GIVEs applying the planar fit, Kriging, and ANN methods for contrastive purposes. For each IGP, ANN data are firstly printed, then followed by the planar fit, and then by Kriging data. Table 3 records the average and maximum values of the bounding error ratios for different methods. The bounding error ratio is defined as the absolute residual divided by its GIVE [37]. 1  10  20  30  40  50  60  70  80  90  100  110  120  130  140  149  0   3   6   9   12   ANN 55 (m)   1  10  20  30  40  50  60  70  80  90  100  110  120  130  140    The Grid Ionospheric Vertical Error Indicator (GIVEI) is the discrete form of GIVE representation for message convenience [11]. Figure 12 presents the GIVE histograms of the planar fit, Kriging, and ANN GIVE methods. The GIVE partitions on the abscissa axis are divided according to GIVEI 3 to 11. Figure 13 shows the GIVEI availability distributions under normal ionospheric conditions. GIVE values below 2.7 m and 6.0 m in percentages are listed in the following Table 4. The ANN method has at least 25 percentage points more on most IGPs for 2.7 m GIVE availability.

Under Disturbed Ionospheric Conditions
Test data under disturbed ionospheric conditions are selected from 2012 and 2013, containing a total of 44 data points. Data for disturbed conditions are sparse for several reasons: total ionospheric disturbed days are rare in the past decade; days before 2016 are selected due to not overlapping in time with training and validation data; all observation files of seven stations need to be intact on selected days.
Due to the small amount of test data under disturbed conditions, reference prediction figures and residual histograms are not convincing. Therefore, we just put forward the figure of predicted GIVDs by various methods as well as their references, data point by data point. Residual indicators MAE and RMSE are presented in tables. Figure 14 shows the predicted GIVDs by multiple methods and the CODE GIM references of each data point. Table 5 records the resulting MAE and RMSE of the fitting residuals under disturbed conditions by applying various methods.   Figure 15 presents the absolute vertical delay residuals and corresponding GIVEs derived by three methods under disturbed ionospheric conditions. Meanwhile, the averaged and maximum bounding error ratios of GIVE under abnormal conditions are shown in Table 6.

Case of Insufficient IPPs
When the number of IPPs around an IGP for GIVD or GIVE derivation is insufficient, GIVDs and GIVEs of IGPs could be estimated if the number of IPPs is larger than N min = 10 by applying the planar fit and Kriging methods [16,20].
In this part, the resolvent of GIVDs and GIVEs by applying the ANN method for insufficient IPPs is given. IGPs 54 (25°N, 110°E), 68 (30°N, 105°E), and 83 (35°N, 105°E) that are distributed on the edge of the selected station network region are tested. Based on historical data, the minimum number of available IPPs around IGP 54 is 25, by applying the search algorithm proposed by Sparks [20]. The minimum numbers of available IPPs around IGP 68 and IGP 83 in history are 20 and 15, respectively.
Since the ANN has a fixed structure for a certain IGP, the input number of the network is also invariant. Therefore, we trim the network inputs for IGPs with insufficient IPPs. Specifically, this is carried out by trimming the dimensions of the input vectors Φ j , Λ j , and I j to the minimum available IPP numbers, and adjusting the number of hidden layer neurons accordingly. As a result, the structures of ANNs for IGP 54, IGP 68, and IGP 83 are 78-26-1, 63-21-1, and 48-16-1.
The experimental results on IGP 54, IGP 68, and IGP 83 by applying the planar fit, Kriging, and ANN methods are given below. When IPP numbers are larger than the minimum in history, the planar fit and Kriging methods adopt all the eligible IPPs. However, the ANN only uses the nearest minimum number of IPPs and abandons the others. Figures 16 and 17 show the absolute values of the vertical delay residuals and corresponding GIVEs under normal and disturbed conditions, respectively. Table 7 presents the MAE and RMSE of the fitting residuals under different ionospheric conditions, and Table 8 records the average and maximum bounding error ratios of GIVEs.  1  10  20  30  40  50  60  70  80  90  100  110  120  130  140

Discussion
Scrutinizing Table 2 and Figure 5, the ANN prediction residuals in the MAE of test data under normal conditions are similar to the validation MAE losses and even better. This phenomenon shows that the ANN captures the fitting regulations by learning numerous data and does not overfit them. The planar fit and Kriging methods perform similarly, which verifies the statement made by Sparks in [20]. The MAE and RMSE residual reductions in the ANN method are approximately 15% relative to planar fit and Kriging under normal ionospheric conditions among tested IGPs.
From Figure 8, which describes the distribution of IPPs, it is noticed that low latitude IGPs 55, 56, and 57 have dispersed surrounding IPPs and deviated IPP focuses. Accordingly, the residual distributions of the planar fit and Kriging methods emerge in a non-Gaussian form in the histograms shown in Figure 7. In contrast, the residual distributions for the other six IGPs at higher latitudes with concentrative and centric surrounding IPPs are more likely to meet the Gaussian distribution. However, the residuals of GIVDs predicted by the ANN method almost perfectly meet the Gaussian distribution regardless of the IPP scattering or eccentricity. Therefore, the ANN method maintains the fitting performance of IGPs at the edges of the SBAS service area. Furthermore, Figure 6 shows that the ANN method achieves the most centralized reference prediction plots and maximal R values, which symbolizes the best correlations between predicted and reference vertical delays.
It is found from Figures 9-11 that, in general, GIVE fluctuates daily, and most residuals also follow the trend. GIVD residuals and GIVE crests often emerge when the solar altitude angle reaches maximum. Since low latitude IGPs 55, 56, and 57 have more active ionospheric states, the top values of GIVEs during the tested period (149 days) are between 9 m and 12 m. In contrast, the maximum GIVEs of other IGPs at higher latitudes are less than 8 m.
From Figures 9-11, it is also noticed that although the ANN method is completely different from the planar fit and Kriging methods in the principle of GIVE derivation, they give similar GIVE estimations at corresponding time points. The consistency of GIVE results makes the methods mutually confirm their rationality. It is also noticed that ANN GIVEs for active ionospheric status (GIVE peaks) are more conservative than other methods. On the contrary, the flat parts of ANN GIVEs are more radical than that of the planar fit and Kriging methods, i.e., the ANN GIVE method provides lower confidence bounds on static ionospheric status. Nevertheless, the average bounding error ratio presented in Table 3 proves that the ANN GIVE method is the most conservative under normal ionospheric conditions among the three. The maximum ratios in Table 3 imply that no residual exceeds the GIVE bounds.
Combining Figures 12 and 13, it could be found that on most IGPs, the ANN GIVE method provides smaller error bounds than other methods. In addition, IGP 71 has perfect IPP distributions in most cases. Therefore, the GIVEs derived by the planar fit and Kriging methods are better than the ANN method.
In the GIVEI availability in Figure 13, the ANN method ascends more rapidly than the planar fit and Kriging methods at the early stage (small GIVEI), while at the later stage (large GIVEI), planar and Kriging methods rise more rapidly and reach 100% earlier. This property accords with Figures 9-11 that show that the ANN GIVE method provides radical bounds for small errors and conservative bounds for large errors.
Under disturbed ionospheric conditions, the prediction accuracy of the ANN method falls behind the other two methods which is corroborated by Figure 14 and Table 5. The Kriging method attains the best accuracy. The fitting residuals of Kriging are 10% lower than that of the planar fit method on most IGPs. It is found from Figure 14 that the ANNpredicted GIVDs are basically consistent with the references. Nevertheless, the prediction errors at several ionospheric delay peaks approach 3 m, which tremendously drag down the total accuracy. Under disturbed conditions, ANNs are unable to make good predictions when large distinctions exist between testing data and training data. The losses of minimal samples are sacrificed so that the network maintains a minimum total loss during the training process.
In Figure 15, the varying GIVE trends of different methods are still consistent, which is reflected by the locations of peaks and troughs. Consequently, the ANN GIVE method is effective under extremely active ionospheric conditions. However, under disturbed conditions, the conservative GIVE of the ANN method is obviously manifested, which will restrict the system's availability. The 2.7 m GIVE availabilities averaged from all IGPs of the planar fit, Kriging, and ANN methods are 14.90%, 18.94%, and 5.05%, and the 6.0 m availabilities are 88.64%, 90.15%, and 72.22%, respectively, under disturbed conditions. In Table 6, the average bounding error ratios are increased compared to normal conditions. The rising ratios in planar fit and Kriging methods are conventional due to the design of the GIVE algorithm. However, in the ANN GIVE method, the ascending average bounding error ratios are deteriorated by anomalously enlarged residuals. Although the maximum bounding error ratios of the ANN method are still much less than 100%, there exist risks in the ANN GIVE method when lacking training data. In addition, it is noted that in order to deal with the influence of ionospheric random disturbance in the SBAS, some schemes have been proposed to mitigate or even overcome, as much as possible, the adverse effects due to mathematical methods (modeling, interpolation, etc.) [38,39].
Based on the results provided by Figures 16 and 17 and Tables 7 and 8, in the case of insufficient IPPs, the fitting accuracy of planar fit and Kriging methods obviously declines, which can be evaluated by MAE and RMSE. However, the accuracy of the ANN method declines slightly due to the assistance of time, solar, and geomagnetic information. In addition, GIVEs derived by multiple methods are even more conservative, and the 2.7 m and 6.0 m GIVE availabilities of the ANN method are better than the other two methods on most IGPs.

Conclusions
Ionospheric delay is a critical error source in Global Navigation Satellite Systems (GNSSs) and a principal aspect of the SBAS corrections. Ionospheric vertical delays on IGPs that broadcast to users are derived from the IPPs observed by SBAS reference stations. Consequently, the ionospheric correction method is extremely important for the accuracy, integrity, and availability of navigation systems. However, traditional ionospheric methods just take the spatial specificities into consideration, ignoring the temporal correlations. In order to improve the system's accuracy and integrity, the ANN-based ionospheric method is proposed in this article. Both temporal features and abundant historical data are applied for network training. Furthermore, based on the distribution of the residual random variable, the ANN-based GIVE algorithm is also derived.
In this study, when the number of IPPs around an IGP is above 30, the ANN model has a 93-31-1 three-layer structure. The 93 input features are 2-dimensional time indication, 1-dimensional geomagnetic and solar activity index derived by applying the Principal Component Analysis (PCA) technique, and 30 by 3-dimensional IPP relative latitudes, relative longitudes, and vertical delays. Other characters of the networks include selecting the sigmoid function as the hidden layer activation function and applying the Adam optimizer as the gradient updating algorithm. For each IGP, the network is trained and validated by a total of 30,130 data points selected from 2016 to 2020. Nine IGPs over the East China region are tested under normal and disturbed ionospheric conditions. GIVDs and accompanying GIVEs are calculated, for comparison with the planar fit and Kriging methods.
The derivation accuracy under normal ionospheric conditions is tested on 3427 data points collected from January to May 2021, while the accuracy under disturbed conditions is tested on 44 data points selected from 2012 and 2013. Results show that the MAE and RMSE of the ANN method are approximately 15% less than the planar fit and Kriging methods. Under disturbed conditions, on the contrary, the ANN is incompetent and shows the worst prediction accuracy, because of the minimal samples in the training set. Kriging retains a good performance and the fitting residuals of Kriging are 10% lower than that of planar fit on most IGPs.
On most tested IGPs under normal conditions, the ANN method provides smaller GIVEs and the residuals meet Gaussian distribution better than the other two methods. Specifically, the percentage of GIVE values below 2.7 m by applying the ANN method has at least 25 points more than Kriging. Under disturbed ionospheric conditions, the ANN GIVE method is still effective and conservative. However, the GIVE availability decline is evident in all GIVE Indicator (GIVEI) levels. Even though the residuals are far from bounding limits, the ANN GIVE method involves risks because of the minimal samples of the training dataset.
When the number of IPPs around an IGP is insufficient, the structures of ANNs are changed for three tested IGPs based on their available IPP number. The accuracy of the ANN method declines slightly, and 2.7 m and 6.0 m GIVE availabilities of the ANN method are better than the other two methods on most IGPs. Nevertheless, the regional and boundary performance of the ANN method in China could not be generalized until a large amount of the service area and all the edge IGPs are tested.
As a result, an ionospheric delay threshold should be set which indicates that the ANN model is not reliable. Training strategy optimization for the ANN method aiming at disturbed ionospheric conditions will be the research direction. GIVD and GIVE on IGPs over China by way of the ANN method will be discussed in the future.  Data Availability Statement: GNSS observation files and navigation message files applied in this article are accessible from the NASA's CDDIS website: https://cddis.nasa.gov/archive/gnss/data/ daily/ (accessed on 26 January 2022). IONEX files are also accessible from the CDDIS website: https://cddis.nasa.gov/archive/gnss/products/ionosphere/ (accessed on 26 January 2022). OMNI data are accessible from NASA's SPDF website: https://spdf.gsfc.nasa.gov/pub/data/omni/low_ res_omni/ (accessed on 26 January 2022).