Kriging with Unknown Variance Components for Regional Ionospheric Reconstruction

Ionospheric delay effect is a critical issue that limits the accuracy of precise Global Navigation Satellite System (GNSS) positioning and navigation for single-frequency users, especially in mid- and low-latitude regions where variations in the ionosphere are larger. Kriging spatial interpolation techniques have been recently introduced to model the spatial correlation and variability of ionosphere, which intrinsically assume that the ionosphere field is stochastically stationary but does not take the random observational errors into account. In this paper, by treating the spatial statistical information on ionosphere as prior knowledge and based on Total Electron Content (TEC) semivariogram analysis, we use Kriging techniques to spatially interpolate TEC values. By assuming that the stochastic models of both the ionospheric signals and measurement errors are only known up to some unknown factors, we propose a new Kriging spatial interpolation method with unknown variance components for both the signals of ionosphere and TEC measurements. Variance component estimation has been integrated with Kriging to reconstruct regional ionospheric delays. The method has been applied to data from the Crustal Movement Observation Network of China (CMONOC) and compared with the ordinary Kriging and polynomial interpolations with spherical cap harmonic functions, polynomial functions and low-degree spherical harmonic functions. The statistics of results indicate that the daily ionospheric variations during the experimental period characterized by the proposed approach have good agreement with the other methods, ranging from 10 to 80 TEC Unit (TECU, 1 TECU = 1 × 1016 electrons/m2) with an overall mean of 28.2 TECU. The proposed method can produce more appropriate estimations whose general TEC level is as smooth as the ordinary Kriging but with a smaller standard deviation around 3 TECU than others. The residual results show that the interpolation precision of the new proposed method is better than the ordinary Kriging and polynomial interpolation by about 1.2 TECU and 0.7 TECU, respectively. The root mean squared error of the proposed new Kriging with variance components is within 1.5 TECU and is smaller than those from other methods under comparison by about 1 TECU. When compared with ionospheric grid points, the mean squared error of the proposed method is within 6 TECU and smaller than Kriging, indicating that the proposed method can produce more accurate ionospheric delays and better estimation accuracy over China regional area.


Introduction
The ionosphere is the upper part of atmosphere located between from 50 km to 1300 km above the Earth's surface. It has a high density of ions and free electrons that can affect the propagation of electromagnetic radio frequency waves [1]. Ionospheric delay is an important source of errors in Global Navigation Satellite System (GNSS) positioning and navigation, which can be determined from Total Electron Content (TEC) measurements. Accurate ionospheric delay corrections can significantly accelerate the convergence of real-time GNSS ambiguity resolution and thus essentially improve precision and performance of positioning and navigation products for single-frequency users [2].
Ionospheric models can be generally classified into two categories: function-based and grid-based [3]. The former represents regional ionospheric TEC by estimating the coefficients of the employed mathematical functions, such as (generalized) trigonometric series functions [4,5], polynomial functions, low-degree spherical harmonic functions and spherical cap harmonic functions [3]. Mathematical function-based models cannot effectively reflect high-frequency variations in local ionosphere, unless the functions employed possess high frequency components. Grid-based ionosphere models have often been adopted by single-frequency users in wide area augmentation systems such as Wide Area Augmentation System (WAAS) and European Geostationary Navigation Overlay Service (EGNOS). Chao (1997) proposed Inverse Distance Weighted (IDW) functions with the Klobuchar model to compute ionospheric delays in WAAS [6]. Komjathy et al. suggested a planar fitting method to interpolate regional ionospheric grid delays [7]. Then WAAS uses a spatial correlation interpolation scheme of Kriging [8], while EGNOS employs a non-uniform partitioning scheme for ionospheric grid [9]. Moreover, grid-based models can be more effective in detail description of local ionosphere variations.
Ionosphere changes spatially and temporally [8]. Kriging (see e.g., [10]), though originating in the field of mining, is developed to fully account for spatio-temporal information on data [10]. Since then, it has become a powerful tool in geostatistics and spatial statistics to handle spatially and/or temporally correlated and irregularly distributed data and has been widely applied to other fields such as hydrology [11], climatology [12], soil science [13], ecology [14], Geo-Information System (GIS) [15], atmosphere science [16], geophysics [17] and geodesy [18]. Since 2002, it has been shown to be efficient for ionospheric delay estimation as well. Blanch conducted extensive experiments to validate the technique of Kriging for ionospheric estimation by using WAAS ionospheric measurements collected during quiet and disturbed periods [8]. Blanch et al. further developed a hybrid algorithm by combining Kriging and tomography for Satellite Based Augmentation Systems (SBAS) and applied it to post-process ionospheric TEC measurements from the US and Brazil [19]. Wielgosz et al. compared Kriging with multiquadratic models by using GPS observations from five Ohio CORS stations; the results have shown that both methods are suitable for instantaneous regional ionosphere modeling [20]. As one of the Internal GNSS Service (IGS) ionosphere analysis centers, Technical University of Catalonia (UPC) has adopted Kriging to re-process existing UPC Global Ionospheric Map (GIM) products. The results have shown that the Root Mean Squared error (RMS) of the UPC Kriging GIM is about 16% lower than the current UPC GIM and about 2% lower than IGS GIM, where the RMS is the root mean square of the difference among the geometry free linear combination observations and the Slant TEC (STEC) computed by each GIM at the same elevation in a continue arch at two different time steps. Both the standard deviation and the RMS are reduced approximately by 0.3 TEC Unit (TECU) (6%) and 0.1 TECU (3%) over the current UPC GIM products when compared with TOPEX/Poseidon and JASON TEC data, respectively [21]. These studies have shown that the Kriging method is useful and effective for ionospheric TEC estimation.
Sayin et al. compared the performance of the Ordinary Kriging (OK) with the Universal Kriging (UK) using a synthetic data set with different variances and different types of sampling patterns. They observed that for small sampling numbers and with higher variability, OK performs better. However, UK gives better results in case of smaller variances in synthetic surfaces and increasing sample number [22]. Although UK takes the random signals of TEC observations into account, it does not consider measurement errors. Li et al. demonstrated, with the Ground-based Regional Integrity Monitoring System (GRIMS) reference stations in China, that UK can produce a more accurate ionospheric delay correction than the distance-weighted method and a tight confidence bound in the boundary areas [23]. The technique has also been well used to reconstruct the ionosphere critical frequency (f oF2) instantaneous mapping [24][25][26]. Relevant studies can be found in Tierno et al. [27], Deviren et al. [28] and Chen [29].
Trend and signal in delays reflect the variations in the ionosphere i.e., the TEC [8]. The accuracy of measurements and the stochastic model of signals are two basic elements in Kriging. However, TEC measurement noise has often been neglected in the estimation of ionospheric delay, as can be found in the research works mentioned above. Additionally, in the real conditions, TEC measurements can be of different accuracy. Nevertheless, Kriging cannot be able to handle the situation with a number of different unknown variance components of spatial/temporal data. Therefore, as the first motivation of this paper, we will consider the estimation of ionospheric delays with TEC measurements of different (unknown) accuracy. On the other hand, even if measurement noise and the stochastic model of the signal would be fully taken into account in the estimation of ionospheric delays, we have to determine whether the level of measurement noise and the stochastic model of signals are correctly given, since incorrect stochastic models can distort the estimation of trends. The second motivation of this paper is to apply the new Kriging method to real data and to demonstrate how to eliminate the distortion by calibrating stochastic models of measurements and signals. More specifically, we will theoretically extend the basic Kriging principle to the case with measurements of different (unknown) variance components and calibrate factors for signals to balance the stochastic models of both measurements and signals.
This paper is organized as follows: in Section 2, we will briefly outline ordinary Kriging algorithms and semivariogram for the construction of signal stochastics. In Section 3, we will theoretically extend Kriging methods to the case in which the stochastic model of measurements contains a number of unknown variance components. The general formulae to estimate trend parameters, signals and the unknown stochastic models of measurements will be worked out in detail. The proposed new Kriging method will then be adapted for use in ionospheric TEC estimation in Section 4 and its implementation will be outlined as well. A brief overview of the common used ionospheric mapping models are given in Section 5 in order to compare with our introduced method. In Section 6, we will apply the new method developed in Sections 3 and 4 to analyze the data from CMONOC and to evaluate the quality of China Regional Ionospheric Maps (CRIM). Finally, the conclusions derived are summed up in Section 7.

The Principle of Kriging Spatial Interpolation
Kriging originates from the field of mining and was developed to deal with spatio-temporally correlated data [30]. Since then, it has become a standard and powerful method in geostatistics and spatial statistics [10] and has found a variety of applications in all the subjects of study where (regularly and/or irregularly) spatially and/or temporally correlated data are routinely encountered. It is based on the variability and spatial correlation of regionalized variables to determine the weights of sampling points distributed around the point to be estimated, according to the principle of unbiased and optimal estimation. Finally, the value of the estimated point is obtained by using the linear combination of data samples. The spatial variability and correlation of data is described by using spatial covariance function or semivariogram. Very often, we do not know the semivariogram in advance. Instead, it must be practically estimated from the original data set. In the remainder of this section, we will closely follow Cressie to briefly outline the basic principle of the ordinary Kriging. In geostatistics, we usually assume the stationarity for spatially distributed data. To be specific, given a spatial random function Z(x), we assume (i) that its expectation is constant and does not depend on the location x and (ii) that the correlation function between any two points depends solely on their distance. These two assumptions can be mathematically described equivalently as follows: (see e.g., Cressie [23]), where the function γ(h) is called semivariogram, which characterizes the spatial correlation of random function Z(x), and h = h is the distance between the two spatial points. If the assumption of isotropy is removed, then h will not be a scalar of distance but should be replaced by the vector between the two points. If we further remove the assumption of homogeneity, in this case, the semivariogram should be written in its most general form as γ(x + h, x). Based on the above assumptions and given a semivariogram function, the basic Kriging method is to find the best linear unbiased estimator or interpolator with minimum variance. Given a set of measurements Z(x 1 ), Z(x 2 ), ...,Z(x N ), we can interpolate the valueẐ(x 0 ) of a given (non-measured) point x 0 by constructing a linear combination of the measurements as follows: where λ i are the unknown coefficients to be determined. In general, we also require that these coefficients be non-negative. The expectation of the difference between the linear interpolation Equation (3) and the signal at the point x 0 can be written as follows: Since we require the linear interpolator Equation (3) be unbiased, and bearing the condition Equation (1) in mind, we must have: Under condition Equation (5), we can further compute the error variance forẐ(x 0 ) − Z(x 0 ), which is given below: where γ x i , x j is the semivariogram between points x i and x j . To construct the optimal interpolator Equation (3), we require minimum error variance Equation (6) under the constraint Equation (5) of unbiasedness. Following Cressie [23], we construct the augmented objective function: where µ is the Lagrange multiplier. By computing the partial derivatives of the objective function Equation (7) with respect to λ i and setting them to zero, we have the normal equations: from which, together with the equality condition Equation (5) of unbiasedness, we can readily obtain the following linear system of equations: By solving the Equation (9), we have the solution of λ i and µ, and as a result, can use them to construct the OK interpolator Equation (3) and to complete the estimation ofẐ(x 0 ). Here γ ij (i, j = 1, 2 ,..., N) stands for γ x i , x j in Equation (8), and γ i0 stands for the semivariogram between the i-th measured point and the point to be interpolated, for the conciseness of notations. The estimated error variance of the optimal OK interpolator can be derived by applying the error propagation law to Equation (3) and is simply given as follows: Equations (9) and (10) can also be written alternatively in matrix form as follows: where: Therefore, the weights λ i and the Lagrange multiplier µ can be estimated by solving Equation (11).

Construction of Semivariogram
A semivariogram describes the spatial correlation of a random field and plays a key role in spatial Kriging interpolation. In practice, a semivariogram is generally unknown but has to be estimated from spatial measurements. To simplify the numerical representation of spatial correlations among random points under the assumption of stationarity and isotropy, the semivariogram in this case is often represented by using a few parameters [30].
By definition, a semivariogram is half of the variance of the difference between Z(x) and Z(x + h), which is related to the covariance function through the following relationship [30]: Given N(h) pairs of measurements separated at a distance h and under the assumption of stationarity and isotropy, the semivariogram can then be numerically estimated as follows: Practically, all the measurement points may not be regularly distributed in space. In this case, we will have to allow a certain tolerance in order to have a sufficient number of data pairs to estimate semivariogram. Nevertheless, a large number of semivariogram values, as given by Equation (16), may not be convenient to use. Thus, one often selects some appropriate and simplified models to fit and represent semivariogram. The most common used semivariogram models include spherical functions, exponential functions and Gaussian functions [8,30].

Variance Component Estimation (VCE) Based on Collocation
Variance component estimation (VCE) has been one of the most important topics in geodesy [31][32][33][34][35]. To correctly determine the weights of different types of measurements and/or measurements of different precisions, we have to simultaneously estimate both the model parameters and variance components. A number of methods have been proposed both in geodesy and statistics to estimate the unknown parameters along with the unknown stochastic model of measurements. The most widely used methods include Helmert quadratic estimation [32], maximum likelihood [33], the Best Invariant Quadratic Unbiased Estimation (BIQUE) [34], the Restricted Maximum Likelihood Method (REML) [36], and the MInimum Norm Quadratic Unbiased Estimation (MINQUE) [37,38]. Extensions to ill-posed problems can be found in Xu et al. [35], Xu [39] and Eshagh [40].
In this section, we focus on the following collocation model: where y is an n × 1 observation vector. In the literature of ionospheric modeling, the notation I(x) is often used instead of y, to denote observations of ionospheric delays. ε is the corresponding measurement error vector; A is the (n × t) design matrix of the unknown parameters, β is a t × 1 vector that contains (deterministic) unknown trend parameters to be estimated; s is an m × 1 vector of random signals with an associated (n × m) design matrix B.
In general, the least-squares collocation method usually assumes that the variance-covariance matrices of the measurement errors and the signals are given and then uses the measurements y to estimate the parameters β and the random signals s. Yang and Xu extended the collocation model to the case that allows one unknown variance component for the measurements y and one unknown variance component for the random signals s [41]. Yang et al. introduced an adaptive factor into a new adaptive collocation procedure and used the maximum likelihood technique to determine the weights of the signals and measurements [42].
In this paper, we will further extend the collocation model to a very general case with a number of unknown variance components for the stochastic models of both the measurements ε and the random signals s. For simplicity, we will focus on the following stochastic model: where the variance-covariance matrix Σ contains Σ ε and Σ s that are the variance-covariance matrices of the measurements and the random signals, respectively. U si and U εi are the given positive (semi-) definite matrices, and σ 2 si (i = 1,2, . . . ,k s ) and σ 2 εi (i = 1,2, . . . ,k ε ) are the unknown variance components. If the measurements and the random signals can be divided into a number of stochastically independent sub-groups, and if the measurements and the random signals are stochastically independent as well, then Σ s and Σ ε become block-diagonal, which can be rewritten as follows: where P si and P εi are the weight matrices, corresponding, respectively, to the i-th sub-group of the random signals s and the measurement errors ε. In this case, the matrices U si and U εi become: and: If the signal vector s at unmeasured points is also included into the collocation model (16), we then have: If the variance components in the stochastic models is known or given with some initial values, we apply the least squares collocation principle (see e.g., Huang 1992) [43] to model (22) and obtain the following solution:β = A T P y A −1 A T P y ŷ s = Σ s B T P y y − Aβ ŝ = Σ s s B T P y y − Aβ (23) where P y = BΣ s B T + Σ ε −1 ,β is the estimated trend parameter vector,ŝ is the estimated signal vector at the known points,ŝ is the estimated signal vector at the unmeasured points.
Since the variance components are unknown, we now apply variance component estimation to the collocation model (22) to estimate the variance components of the measurements and the random signals. Specifically, we will use the MINQUE method in this section.
To estimate the variance components for the measurements and the random signals, and keeping in mind that the variance component estimation has nothing to do with the unobserved signals, we rewrite the stochastic signals as pseudo observations. Thus, the collocation model (17) can be alternatively presented as follows: where y s and ε s denotes the prior values and errors of the signals, respectively.
Given some initial values σ 2 s0 = σ 2 s10 σ 2 s20 · · · σ 2 sk s 0 T and σ 2 ε0 = σ 2 ε10 σ 2 ε20 · · · σ 2 εk ε 0 T for the unknown variance components, the MINQUE estimate of the variance components can be computed by using the following equations (see e.g., [35,37,43]): where: (28) and the elements of the coefficient matrix G and vector q are given by: tr(·) denotes the trace of a square matrix, and: Σ 0 stands for Σ with initial sets of σ 2 s0 and σ 2 ε0 . After obtaining the estimates ofσ 2 s0 andσ 2 ε0 , we can re-estimateβ,ŝ andŝ in Equation (23) and further use these Equations (22)- (35) to iteratively estimate the variance components. Then the estimationsŷ at unmeasured points can be get by the following Equation (36):ŷ = A β + B ŝ (36) and the corresponding variance-covariance matrix Σŷ is: where A and B are the design matrices, and the formulas of variance matrices in Equation (37) can refer to Huang [43] in detail.

Implementation of VCE Integrating with Kriging in Ionospheric Delay Estimation
The observable of STEC can be generated by forming the geometry-free linear combination of pseudorange and phase data. The Vertical TEC (VTEC) can then be derived from STEC at the corresponding ionospheric Intersecting Pierce Point (IPP) by multiplying with the mapping function [44]. In this paper, the Modified Single-Layer Model (MSLM) mapping function is employed [44]. The variance of the TEC observations at different elevations is calculated using the following formula: where e is the elevation angle of the satellite (the cut-off elevation angle is set to 15 • , Σ εi is the corresponding variance of the ionospheric delay observation error; σ 2 0 is the prior variance of the the TEC estimation noise. For ionospheric delay estimation, we set the trend in the Equation (17) to be constant. Moreover, the matrix B is an identity matrix. The observation equation for the vertical ionospheric delay at the IPP located at x with x ∈ φ, λ, denoted by I(x), becomes: For the sake of simplicity in the ordinary Kriging interpolation, we directly use the VTEC observed at IPPs within a certain limit around the Ionospheric Grid Points (IGPs) to interpolate the VTEC according to the algorithm in Section 2.1.
In order to apply Kriging interpolation and VCE for ionospheric estimation, we need to choose some appropriate semivariogram or variance function which can give exact veracious description for ionosphere random behavior. Figure 1 presents the semivariograms of the ionosphere on day-of-year (DOY) 305 in 2014 (2014.11.01) from UTC 00:00 to 22:00 with a 2-h interval which are computed by using Equation (16) based upon CMONOC observational data sets. The dots and curves shown in Figure 1 are the experimental semivariograms and fitting results respectively. According to the shape and behavior of the experimental semivariogram, the Gaussian semivariogram function of the signals is selected: (see e.g., [30]), where C 0 , C and a are the unknown parameters to be estimated. In geostatistics, C 0 is called "nugget" which is linked to the continuity and to the spatial regularity of the regionalized variables. The physical explanations in ionosphere study, the discontinuity at the origin are expounded due to the obliquity error, remaining bias in the measurements and the fast-changing ionosphere [8]. The semivariogram reaches a limiting value, C 0 ± C, termed as "sill" on behalf of the structural variance of spatial variation, and reflect the biggest variation extent of the variables. The two observations are uncorrelated if the distance is beyond a limit. This distance is called the "range". This value is equal to √ 3a for Gaussian function model [30]. These parameters are estimated by using the data in the corresponding interval of time to characterize the spatial correlation and variability. Therefore, for each epoch map, a set of new semivariogram parameters (nugget, sill, range) is obtained which give descriptions for the behavior of the ionospheric TEC. Figure 2 indicates how the semivariogram changes during the campaign days.
The prior variance-covariance of the signals can be computed according to Equation (15) after fitting the experimental semivariogram to the theoretical semivariogram model (40) by using the least squares estimation. Once the semivariogram and variance are computed, with the known and typically diagonal measurement noise matrix, the Kriging and VCE methods can be used to estimate the results and the estimation variance. Using the developed algorithms, the maps for regional TEC in China can be obtained in real-time with the fitted semivariogram function automatically at any desired epoch.
is obtained which give descriptions for the behavior of the ionospheric TEC. Figure 2 indicates how the semivariogram changes during the campaign days.
The prior variance-covariance of the signals can be computed according to Equation (15) after fitting the experimental semivariogram to the theoretical semivariogram model (40) by using the least squares estimation. Once the semivariogram and variance are computed, with the known and typically diagonal measurement noise matrix, the Kriging and VCE methods can be used to estimate the results and the estimation variance. Using the developed algorithms, the maps for regional TEC in China can be obtained in real-time with the fitted semivariogram function automatically at any desired epoch.  In order to estimate the VTEC at a given IGP at a given epoch, we must first select a set of VTEC measurements whose IPPs are distributed within the restricted scope centered upon this IGP. We set the maximum radius and minimum radius as max R = 2000 km and min R = 500 km respectively. It is important to note that the appropriate radius R for searching IPPs at a given IGP must be within a half of the given searching limits and half of the range, as well as the max R . The minimum and maximum number of IPP measurements is set as min N = 5 and max N = 25 based on our experiments. In order to estimate the VTEC at a given IGP at a given epoch, we must first select a set of VTEC measurements whose IPPs are distributed within the restricted scope centered upon this IGP. We set the maximum radius and minimum radius as R max = 2000 km and R min = 500 km respectively. It is important to note that the appropriate radius R for searching IPPs at a given IGP must be within a half of the given searching limits and half of the range, as well as the R max . The minimum and maximum number of IPP measurements is set as N min = 5 and N max = 25 based on our experiments. If there are IPPs fewer than N min within R, no corresponding estimate is computed. If we find the N max IPPs around the estimated IGP, then we stop searching. On the other hand, due to the variation level of ionospheric activity which calls for a change for the degree of the ionosphere measurements correlation, the experimental semivariograms have to be fitted in real-time with the IPPs located within appropriate searching limits, with the distance interval to compute the semivariogram from data and its tolerance set to 100 km and 50 km, respectively. The searching limits here refer to the maximum distance between any two IPPs with correlation. From Figure 1, we can see that in the disturbed ionospheric condition, the semivariograms have smaller range and larger sill than those in the quite ionospheric condition. In the measured data tests, this value of searching limits is set as 3500 km between 14:00 and 22:00 in China local time (=UTC + 8 h), and 4500 km at other moments, as shown in blue vertical dashed lines in Figure 1.

Regional Ionospheric Models
For the purpose of numerical comparisons in this paper, we will briefly outline some ionospheric interpolation methods used for ionospheric mapping, namely, polynomial interpolations and interpolations with low-degree spherical harmonic functions and spherical cap harmonic analysis.
The polynomial ionospheric interpolation method uses a polynomial function of the latitude differentials and the hour angles differentials of the Sun to interpolate ionospheric delays. The mathematical model can be written as follows: (see e.g., [44,45] (41) are set to n = m = 8 for regional modeling and n = m = 1 for local interpolation, respectively. The functional method with low-degree spherical harmonic functions is to represent an ionospheric delay measurement with the latitude-dependent associated Legendre functions and the sum of the longitude-dependent sine and cosine terms. The mathematical expression of a VTEC measurement I v at an IPP with low-degree spherical harmonic functions is given by [3,44]: where ϕ is the geographic latitude of the IPP, s = λ − λ 0 is the sun-fixed longitude of the IPP, λ is the longitude of the IPP, λ 0 is the longitude of the Sun, n max is the maximum degree of the spherical harmonic expansion, P nm is the normalized associated Legendre function of degree n and order m, and C nm and S nm are the spherical harmonic coefficients to be estimated. The degree and order define the resolution of the model, and should be coordinated with the scope of region and the number of measurements [44]. In this paper, the maximum degree is set to 5 experimentally and the total number of the unknown coefficients is (n max + 1) 2 = 36.
Sensors 2017, 17, 468 12 of 23 The regional model on the basis on spherical cap harmonic analysis has also been routinely used in ionospheric modeling [3]. The model with a spherical cap in interval [0, θ 0 ] for regional mapping of VTEC is expressed as: where I v (ϕ c , λ c ) is the VTEC measurement at IPP (ϕ c , λ c ) inside the spherical cap, K max , M and n k(m) are the maximum degree, order of the series and the non-integral degrees, respectively. k is the index of the non-integral degrees n k(m) (0 ≤ k ≤ K max ) computed by an iterative bisection solution. P(cos θ c ) is the normalized associated Legendre function, C km and S km are the unknown normalized spherical cap harmonic coefficients. If the geographic coordinate of the pole of the spherical cap is (θ N , λ N ), a point with the geographic coordinates (ϕ, λ) can be transformed into the spherical cap coordinates as follows (see e.g., [3]): where the co-latitudes θ = 90 • − ϕ, and θ c = 90 • − ϕ c . The modeling area ranges from 5 • to 60 • in latitude and from 70 • to 140 • in longitude, so the pole of the spherical cap (θ N , λ N ) and the half angle θ 0 are determined as (35 • , 105 • ) and 30 • , respectively. In the next section of applications, K max and M are respectively set to 8 and 6 according to Equation (45) to reach the tradeoff between the measurement resolution ω min and computational load, and the total number of the unknown coefficients is More mathematical information on spherical cap harmonic function and the n k(m) solution can be found in Liu [3] and the references therein:

Applications to CMONOC Data and Result Analysis
We have applied Kriging with unknown variance components to reconstruct the ionospheric maps over China, which will be abbreviated as KVCE. We will also compare different methods for ionospheric reconstruction. The local POLYnomial Interpolation (named as IPOLY) approach is also employed. The regional function-based models used include the low-degree SPHerical harmonic function model (SPH), POLYnomial function model (POLY) and Sphere Cap Harmonic Analysis model (SCHA).
The Crustal Movement Observation Network of China (CMONOC) consists of about 260 GPS stations. To demonstrate the construction of a regional real-time ionospheric model over China, 80 stations with a reasonable level of a uniform distribution from CMONOC are selected with data in November 2014 (from DOY 305 to 334). The processing interval is 60s with a cutoff angle of 15 • due to high level of errors in low elevation angles. The regional range is (5 • -60 • N, 70 • -140 • E) with resolution of 2.5 • in latitude and 5 • in longitude with a total of 345 grids. Figure 3 shows the distributions of the selected reference stations (red pentagram points) and IPPs (green points). In the processing, due to the lack of P1 observations at most stations, the CA code pseudoranges are used to extract the ionospheric delays by carrier-to-code leveling process. The satellites differential code biases (DCBs) are calibrated using the products processed by Wuhan University [45], and the receivers DCBs are estimated using the GIM products and the known satellites DCBs. distributions of the selected reference stations (red pentagram points) and IPPs (green points). In the processing, due to the lack of P1 observations at most stations, the CA code pseudoranges are used to extract the ionospheric delays by carrier-to-code leveling process. The satellites differential code biases (DCBs) are calibrated using the products processed by Wuhan University [45], and the receivers DCBs are estimated using the GIM products and the known satellites DCBs.

TEC Maps
The ionospheric delay effect can be reflected by the values of TEC. Therefore, consecutive VTEC maps may allow detecting the local temporal-spatial variations in ionosphere. Based on the observations from CMONOC, China regional VTEC products are generated by applying the ionospheric TEC models mentioned above. The ionospheric activities during the campaign days have similar characteristics that the diurnal variations of the regional mean VTEC distribute in wave shapes with different maximum and minimum. The mean VTEC ranges from 10 to 80 TECU during campaign days. It reaches the maximum values of 55 to 80 TECU during the session from LT 12:00 to

TEC maps
The ionospheric delay effect can be reflected by the values of TEC. Therefore, consecutive VTEC maps may allow detecting the local temporal-spatial variations in ionosphere. Based on the observations from CMONOC, China regional VTEC products are generated by applying the ionospheric TEC models mentioned above. The ionospheric activities during the campaign days have similar characteristics that the diurnal variations of the regional mean VTEC distribute in wave shapes with different maximum and minimum. The mean VTEC ranges from 10 to 80 TECU during campaign days. It reaches the maximum values of 55 to 80 TECU during the session from LT 12:00 to LT 14:00, and approaches the minimum values around 10 to 15 TECU in the local evening and before dawn. Figure 4 presents the examples of China regional instantaneous ionosphere VTEC maps by using the adopted models at the representative epoch with a high solar activity, LT 14:00, when the maximum ionospheric TEC appears commonly. The TEC value can reach nearly 100 TECU in low latitudes. In general, the phenomena that intense ionosphere variations and high values of ionospheric TEC occur frequently in low latitudes are attributed to the geomagnetic disturbances. While in middle latitudes, the ionospheric TEC varies smoothly, ranging from 30 to 50 TECU. Due to the lack of enough IPPs in the low latitude boundary areas, there are no interpolated estimations for KVCE, OK and IPOLY models. That is why there are blank areas in the first three maps of Figure 4. This situation restricts the correction effect of interpolation methods to some extent. Comparatively, as shown in the latter three maps of Figure 4, the function-based models, namely, POLY, SPH and SCHA, can provide the VTEC estimations for these areas, which are computed with their coefficients. But the map obtained from KVCE is similar to that of OK and seems to be a bit smoother compared with the other four approaches.
KVCE, OK and IPOLY models. That is why there are blank areas in the first three maps of Figure 4. This situation restricts the correction effect of interpolation methods to some extent. Comparatively, as shown in the latter three maps of Figure 4, the function-based models, namely, POLY, SPH and SCHA, can provide the VTEC estimations for these areas, which are computed with their coefficients. But the map obtained from KVCE is similar to that of OK and seems to be a bit smoother compared with the other four approaches.  Figure 5 shows the time series of China regional mean VTEC values obtained from different models. In order to avoid overlapping of resultant curves, they are plotted into two groups to conveniently distinguish the differences among them. The distributions show similar ionospheric behaviors in spite of somewhat differences. The maximum VTEC, which can reach nearly 80 TECU, occurs during LT 12:00 to LT 14:00 due to intense activities of ionosphere. Then the VTEC gradually decreases to the minimum, ranging from 10 to 15 TECU. The overall mean and standard deviations of VTEC values for KVCE, OK, IPOLY, POLY, SPH and SCHA during the selected campaign days are equal to (28.2, 14.2) TECU, ( decreases to the minimum, ranging from 10 to 15 TECU. The overall mean and standard deviations of VTEC values for KVCE, OK, IPOLY, POLY, SPH and SCHA during the selected campaign days are equal to (28.2, 14.2) TECU, (27.8, 14.1) TECU, (33.0, 17.0) TECU, (33.7, 17.6) TECU, (34.1, 18.3) TECU and (33.3, 17.7) TECU, respectively. Since the mean VTEC values are different within 6 TECU, and the standard deviations are different by 4 TECU, we may conclude that all the methods under study achieve good agreements of results. Moreover, the results of the lowest mean value of OK may indicate that OK underestimates the ionospheric TEC a little more than KVCE, as compared to IPOLY, POLY, SPH and SCHA. Furthermore, the nearly same standard deviations of VTEC mean values of 14.2 TECU for KVCE and 14.1 TECU for OK indicate almost equivalent stabilities for both methods. The standard deviations of mean VTEC for other methods range from 17.0 to 18.3 TECU and are larger than those of KVCE and OK by about 3 TECU. Owing to the fact that KVCE neither underestimates nor overestimates the VTEC with high stability, this approach seems to be very promising for specifying the TEC over the investigated area.

Local Interpolation Residuals Analysis
The strategy of cross-validation is used to evaluate interpolation performances of different methods of KVCE, OK and IPOLY. For each IPP, its interpolated value of TEC is calculated in the same way as the procedure of IGP TEC estimation. This means that while polynomial is used for the interpolation, only those IPPs around the target IPP are used to calculate the polynomial coefficients. KVCE and OK are performed by following the same procedure. Taking

Local Interpolation Residuals Analysis
The strategy of cross-validation is used to evaluate interpolation performances of different methods of KVCE, OK and IPOLY. For each IPP, its interpolated value of TEC is calculated in the same way as the procedure of IGP TEC estimation. This means that while polynomial is used for the interpolation, only those IPPs around the target IPP are used to calculate the polynomial coefficients. KVCE and OK are performed by following the same procedure. TakingÎ interp to denote the STEC by interpolated VTEC conversion using the MSLM mapping function and I to denote the observed STEC at the same IPP, we can then compute Interpolation Root Mean Squared error (IRMS): Figures 6 and 7 present the IRMS results from KVCE, the ordinary Kriging and polynomial interpolation for DOY 305 to 334. The IRMS from the KVCE ranges from 0.5 to 2.5 TECU, while the IRMS values from the OK and IPOLY are from 1.5 to 3.5 TECU and 1.5 to 3.0 TECU, respectively. As shown in Figure 6, the standard deviations of the IMRS for KVCE, OK and IPOLY are 0.29 TECU, 0.40 TECU and 0.22 TECU, respectively. The comprehensive comparisons about the mean and standard deviation of IRMS indicate that KVCE is superior to OK and IPOLY. In addition, the overall statistical IRMS information can be seen obviously in Figure 7. The IRMS daily mean values of the three approaches are 1.37 TECU, 2.59 TECU and 2.06 TECU, respectively. The mean IRMS of KVCE is smaller than those of OK and IPOLY by about 1.2 TECU and 0.7 TECU, respectively. Besides, the standard deviation of the daily mean IRMS for KVCE is 0.08 TECU, achieving a better stability and/or precision than OK (0.18 TECU) and IPOLY (0.11 TECU). The polynomial interpolation fits the trend with the weights of the measurements, which is effective for estimating the deterministic part (trend) of TEC distribution in local areas; nevertheless, it does not use the random signals in spatial variations of TEC. Since OK uses a constant model as the trend, and since a trend model plays a dominating role in TEC modeling, this explains why the results of local IPOLY with a more flexible trend model are better than those of OK. On the other hand, KVCE not only gives full consideration to the balance between the weights of measurement noises and random signals but also the spatial relationships of the scattered IPP measurements. This probably explains why KVCE performance is superior over OK and IPOLY in local interpolation.
precision than OK (0.18 TECU) and IPOLY (0.11 TECU). The polynomial interpolation fits the trend with the weights of the measurements, which is effective for estimating the deterministic part (trend) of TEC distribution in local areas; nevertheless, it does not use the random signals in spatial variations of TEC. Since OK uses a constant model as the trend, and since a trend model plays a dominating role in TEC modeling, this explains why the results of local IPOLY with a more flexible trend model are better than those of OK. On the other hand, KVCE not only gives full consideration to the balance between the weights of measurement noises and random signals but also the spatial relationships of the scattered IPP measurements. This probably explains why KVCE performance is superior over OK and IPOLY in local interpolation.   precision than OK (0.18 TECU) and IPOLY (0.11 TECU). The polynomial interpolation fits the trend with the weights of the measurements, which is effective for estimating the deterministic part (trend) of TEC distribution in local areas; nevertheless, it does not use the random signals in spatial variations of TEC. Since OK uses a constant model as the trend, and since a trend model plays a dominating role in TEC modeling, this explains why the results of local IPOLY with a more flexible trend model are better than those of OK. On the other hand, KVCE not only gives full consideration to the balance between the weights of measurement noises and random signals but also the spatial relationships of the scattered IPP measurements. This probably explains why KVCE performance is superior over OK and IPOLY in local interpolation.

Regional Modeling Residuals Analysis
POLY, SPH and SCHA can be used to make regional TEC maps. The IPPs' STEC can be computed with the fitted coefficients. Thus, we can also compute the residuals. To compare the performances of POLY, SPH and SCHA, we will use the following measure: where theÎ model is the STEC computed from the models above. For KVCE, only IGP VTEC is estimated locally and VTEC grid models will be constructed here. With the estimates of the IGP VTEC, the IPPs with STEC surrounding IGP within a certain radius can be obtained sequentially. Then the KVCE residuals are computed. For the VTEC grid models from KVCE, the differences between the estimated IPP STEC and the observed IPP STEC are used to compute the modeling precision as follows: where theÎ KVCE is the j-th estimated IPP STEC, M i is the number of IPPs surrounding the i-th IGP, N IGP is the number of IGPs. Figure 8 plots the mean Modeling RMS (MRMS) with the results from KVCE, POLY, SCHA and SPH. The MRMS from KVCE ranges from 0.5 to 3 TECU, with the mean value of 1.49 TECU, while the MRMS values from POLY, SPH and SCHA are about 2.57 TECU, 2.75 TECU and 2.47 TECU, respectively. Obviously, the modeling fitting of KVCE is better than those of SCHA, POLY and SPH methods. This is reasonable because KVCE grid models take local fitting around IGPs while POLY, SCHA and SPH make regional fitting over China. Even if POLY is used both for local interpolation for IGPs with local IPPs around it and for regional fitting over China area, the fitting residuals are at a similar level of around from 2 to 4 TECU (compare Figures 6 and 8). SCHA and SPH have a similar residual level as POLY in the regional fitting over China.

 
where the model Iˆ is the STEC computed from the models above.
For KVCE, only IGP VTEC is estimated locally and VTEC grid models will be constructed here. With the estimates of the IGP VTEC, the IPPs with STEC surrounding IGP within a certain radius can be obtained sequentially. Then the KVCE residuals are computed. For the VTEC grid models from KVCE, the differences between the estimated IPP STEC and the observed IPP STEC are used to compute the modeling precision as follows:  Obviously, the modeling fitting of KVCE is better than those of SCHA, POLY and SPH methods. This is reasonable because KVCE grid models take local fitting around IGPs while POLY, SCHA and SPH make regional fitting over China. Even if POLY is used both for local interpolation for IGPs with local IPPs around it and for regional fitting over China area, the fitting residuals are at a similar level of around from 2 to 4 TECU (compare Figures 8 and 6). SCHA and SPH have a similar residual level as POLY in the regional fitting over China.

TEC RMS Maps
After estimating the ionospheric delays at IGPs, we can further compute the estimation accuracy, namely, RMS or grid sigma σ. The estimation accuracy of KVCE, POLY, SCHA and SPH can be computed by applying the variance-covariance propagation law expressed by Equation (37). The estimation accuracy of OK can be obtained by Equation (12). With the data from this CMONOC network, the daily mean RMS values of the POLY, SCHA and SPH are equal to 10.5 TECU, 9.05 TECU and 9.44 TECU, respectively. These numbers are significantly larger than those of KVCE and OK. Therefore, we focus on the comparison of results between KVCE and OK in this section. Figure 9 presents the nephograms of estimation accuracy over China at four different local time epochs, morning (LT 08:00), noon (LT 14:00), afternoon (LT 18:00) and evening (LT 22:00) on DOY 305, with the subplots in the left and right columns of Figure 9 showing the IGPs RMS results of OK and KVCE, respectively. Due to the lack of IPPs at the boundary, the IGPs in low latitudes have no interpolation values in the blank area shown in Figure 9, which will be called as invalid grids in the following. Figure 9 presents the nephograms of estimation accuracy over China at four different local time epochs, morning (LT 08:00), noon (LT 14:00), afternoon (LT 18:00) and evening (LT 22:00) on DOY 305, with the subplots in the left and right columns of Figure 9 showing the IGPs RMS results of OK and KVCE, respectively. Due to the lack of IPPs at the boundary, the IGPs in low latitudes have no interpolation values in the blank area shown in Figure 9, which will be called as invalid grids in the following. The accuracy of the estimation results depends on the quantity of the measurements and the degree of ionospheric activity. In general, the accuracy in the central region is higher than that along the boundary. For the subplots in the left column of Figure 9, the OK RMS at LT 08:00 and LT 22:00 The accuracy of the estimation results depends on the quantity of the measurements and the degree of ionospheric activity. In general, the accuracy in the central region is higher than that along the boundary. For the subplots in the left column of Figure 9, the OK RMS at LT 08:00 and LT 22:00 is from 0 to 8 TECU, and the characteristics of IGPs RMS distribution are inconspicuous for the quiet ionospheric activities, when compared with the IPPs distribution. At LT 14:00 and LT 18:00, the RMS is within 2 to 20 TECU, and the magnitude of IGPs RMS depends on the density of the IPPs. The pattern of these maps is apparently the same as that of the IPPs distributions. This result indicates that the RMS at noon and in the afternoon is larger than that in the morning and evening, and has evident variation along with the distribution of IPPs. The major reason for this phenomenon is the disturbed ionosphere condition in the daytime, especially at noon and in the afternoon. The KVCE RMS results, as shown in the right column of Figure 9, are within 6 TECU, which are better than those of OK interpolation. The abnormal phenomena appear at LT 14:00 in the KVCE ionosphere map, which shows that the central area has a larger RMS than the boundary areas. A possible reason is likely due to the non-convergence of the VCE estimates. The accuracy obtained in the morning and at night is better than that at noon and in the afternoon. In general, KVCE is better than OK, no matter whether the points are in the boundary area. Moreover, occurrence of the boundary errors have a big impact on the results which are presented in the Figure 9, as well as the following Figures 10 and 11. Certainly, if the points in the boundary area are not taken into account, the differences between Kriging and KVCE will be smaller. Nevertheless, we believe that to fairly compare both methods, it is more reasonable to count on all the points. In other words, the advantage of KVCE becomes even clearer in the boundary area.   Figures 10 and 11 plot the distribution of time series for the epoch mean RMS and daily mean RMS of both methods. Although the patterns of the mean RMS of KVCE and OK are similar, KVCE is clearly significantly better than OK, as can also be further confirmed from the standard deviation statistics plotted in Figure 12. The daily mean RMS of KVCE is less than 1 TECU and that of OK is   Figures 10 and 11 plot the distribution of time series for the epoch mean RMS and daily mean RMS of both methods. Although the patterns of the mean RMS of KVCE and OK are similar, KVCE is clearly significantly better than OK, as can also be further confirmed from the standard deviation statistics plotted in Figure 12. The daily mean RMS of KVCE is less than 1 TECU and that of OK is  Figures 10 and 11 plot the distribution of time series for the epoch mean RMS and daily mean RMS of both methods. Although the patterns of the mean RMS of KVCE and OK are similar, KVCE is clearly significantly better than OK, as can also be further confirmed from the standard deviation statistics plotted in Figure 12. The daily mean RMS of KVCE is less than 1 TECU and that of OK is about 4 TECU. The former is about four times better than the latter in terms of RMS. The abnormal phenomena appear at LT 14:00 in the KVCE ionosphere map, which shows that the central area has a larger RMS than the boundary areas. The possible reason is likely due to the non-convergence of the VCE estimates. Furthermore, the accuracy estimated in the morning and at night is better than that at noon and in the afternoon. Figure 13 also presents the Kp index and SunSpot Number (SSN) showing that the daily sum Kp index and SSN are lower than 30 and 180, respectively. Generally, the ionospheric activity is relatively stable for the campaign days. It is nice confirmation of the KVCE method which gives much more reliable estimations as shown in Figure 12.

Conclusions
We have proposed a Kriging method with unknown variance components to interpolate ionospheric delays for use in real-time, which is an extension of the works by Xu et al. [35] and Xu [39]. As a result, we are able to correctly determine the weighting factors of measurement noises and random signals. The proposed new method has been applied to GPS data collected in 2014 from the Crustal Movement Observation Network of China (CMONOC) and compared with other methods such as ordinary Kriging, polynomial interpolation, low-degree spherical harmonic function models, polynomial function models and spherical cap harmonic analysis models in terms of TEC maps comparisons, local interpolation, regional modeling and VTEC accuracy at ionospheric grid points.
The results show that the ionospheric vertical structures obtained by applying these methods reach a good agreement, with the daily mean VTEC variations ranging from 10 to 80 TECU and the overall mean values from 27.8 to 34.1 TECU during the campaign days. The mean VTEC value from

Conclusions
We have proposed a Kriging method with unknown variance components to interpolate ionospheric delays for use in real-time, which is an extension of the works by Xu et al. [35] and Xu [39]. As a result, we are able to correctly determine the weighting factors of measurement noises and random signals. The proposed new method has been applied to GPS data collected in 2014 from the Crustal Movement Observation Network of China (CMONOC) and compared with other methods such as ordinary Kriging, polynomial interpolation, low-degree spherical harmonic function models, polynomial function models and spherical cap harmonic analysis models in terms of TEC maps comparisons, local interpolation, regional modeling and VTEC accuracy at ionospheric grid points.
The results show that the ionospheric vertical structures obtained by applying these methods reach a good agreement, with the daily mean VTEC variations ranging from 10 to 80 TECU and the overall mean values from 27.8 to 34.1 TECU during the campaign days. The mean VTEC value from KVCE is close to that from OK, which is lower by about 6 TECU than IPOLY, POLY, SPH and SCHA

Conclusions
We have proposed a Kriging method with unknown variance components to interpolate ionospheric delays for use in real-time, which is an extension of the works by Xu et al. [35] and Xu [39]. As a result, we are able to correctly determine the weighting factors of measurement noises and random signals. The proposed new method has been applied to GPS data collected in 2014 from the Crustal Movement Observation Network of China (CMONOC) and compared with other methods such as ordinary Kriging, polynomial interpolation, low-degree spherical harmonic function models, polynomial function models and spherical cap harmonic analysis models in terms of TEC maps comparisons, local interpolation, regional modeling and VTEC accuracy at ionospheric grid points.
The results show that the ionospheric vertical structures obtained by applying these methods reach a good agreement, with the daily mean VTEC variations ranging from 10 to 80 TECU and the overall mean values from 27.8 to 34.1 TECU during the campaign days. The mean VTEC value from KVCE is close to that from OK, which is lower by about 6 TECU than IPOLY, POLY, SPH and SCHA but with a better stability of around 3 TECU. The VTEC maps derived from KVCE and OK are smoother than those of other four methods, but OK tends to underestimate VTEC a little more than KVCE. This fact indicates that KVCE seems to be more promising for reconstructing the TEC over China regional area. In addition, the local interpolation precision of our proposed method is 0.5-2.5 TECU with a mean value of 1.37 TECU, which is smaller than that of the ordinary Kriging and polynomial interpolation by about 1.2 TECU and 0.7 TECU, respectively. The regional modeling precision from KVCE ranges from 0.5 to 3 TECU with an overall mean value of 1.5 TECU, which is smaller than those from the function-based models by about 1 TECU. The estimation accuracy at ionospheric grid points from our new method remains within 6 TECU, with a daily mean of 0.74 TECU and a standard deviation of 0.4 TECU. The results are better than those from the ordinary Kriging, which can even reach 20 TECU under intense solar activities. In addition, the advantage of KVCE becomes even clearer than OK in the boundary area when all the points are counted. The comprehensive analysis results in terms of TEC maps, interpolation, modelling and estimation accuracy with ionospheric grid points have clearly shown that the proposed Kriging method with variance components has the best performance and can produce more rational, and accurate ionospheric TEC than all the other methods used for comparison in this paper.