Dual-parameter regularization method in three-dimensional ionospheric reconstruction

Ionospheric tomography based on the total electron content (TEC) data along the ray path from Global Navigation Satellite Systems (GNSS) satellites to ground receivers is a typical ill-posed inverse problem. The regularization method is an effective method to solve this problem, which incorporates prior constraints to approximate the real ionospheric variations. When two or more prior constraints are used, the corresponding multiple regularization parameters are introduced in the cost functional. Assuming that the ionospheric spatial variations can be separable in the horizontal and vertical directions, different prior constraints are used in each direction, and the dual-parameter regularization algorithm is established to reconstruct the three-dimensional ionospheric electron density in the present paper. To make the reconstruction results comprehensively reflect the observation information and background (prior) information, it is crucial to determine the optimal regularization parameters. The linear model function method is used to choose these regularization parameters. Both an ideal test and a real test show that this regularization algorithm can effectively improve the background model output.


Introduction
The ionosphere is an important part of the earth's environment, significantly influencing the propagation of electromagnetic waves through reflection and absorption.It is generally accepted that radio waves up to 10 GHz can be affected by the ionosphere to some extent when they propagate through the ionosphere.
The ionosphere has extremely complex temporal and spatial variations.Nowadays, the amount of ionospheric measurements steadily increases, and their accuracies continually improve.However, when ionospheric data are collected into certain temporal-spatial bins, some bins have rather sparse measurements, or even have no measurements.
The Global Navigation Satellite Systems (GNSS) ground beacon receiver has the advantages of being low cost, having wide distribution, and having operational simplicity.It can provide a great deal of ionospheric total electron content (TEC) data along the ray path.Moreover, the threedimensional electron density can be reconstructed through the ionospheric tomography technique by using these TEC data (e.g., Austen et al., 1988;Fougere, 1995;Na and Lee, 1991;Raymund et al., 1990), which can greatly enrich the ionospheric data resource.Due to the limitations of the receiver-satellite geometry, the TEC observation in the horizontal direction is limited and the measurement is incomplete, and the ionospheric tomography is a typical ill-posed problem.
The regularization method is an effective method to solve this ill-posed problem by incorporating some prior constraints to approximate the real ionospheric electron density variations.The classical Tikhonov regularization method uses a single constraint to treat ill-posed problems, and naturally has a single regularization parameter.The regularization parameter is applied to balance the weights between background information and real measurements, and different regularization parameters can lead to different reconstruction results.Many methods have been proposed to determine the regularization parameter, such as unbiased predictive risk estimation (UPRE; Mallows, 1973), generalized crossvalidation (GCV; Golub et al., 1979), the L-curve method (Hansen and O'Leary, 1993), and the damped Morozov discrepancy principle (Kunisch, 1993).Chen et al. (2008) have analyzed the superiority of the multiparameter regularization over the single-parameter regularization.When two or more prior constraints are imposed on the cost functional, the reconstruction accuracy may be improved further, and multiple regularization parameters are introduced accordingly.The question of how to optimally determine multiple regularization parameters is an important research avenue in the study of regularization algorithms.The recently proposed model function method is a simple and practical way to choose multiple regularization parameters.
In this paper, the spatial variations of the ionosphere are separated into the horizontal and vertical directions, and a dual-parameter regularization algorithm is established by incorporating different constraints in each direction.The linear model function method is adopted to determine the optimal regularization parameters, and the ideal test and real test are carried out to validate the effectiveness of this algorithm.
2 Data and regularization method

Data
The reconstruction area covers 35-50 • N in latitude, 280-300 • E in longitude, and 100-1000 km in altitude.The spatial interval is 0.5 • in latitude, 1 • in longitude, 20 km from 100 km to 500 km, and 50 km from 500 km to 1000 km in altitude.There are 18 000 grids in total.
The dual-frequency GNSS receiver can provide continuous phase and pseudorange observations with a sample interval of 30 s. Ionospheric TEC can be derived by using phase observations to smooth pseudorange observations.The elevation cutoff angle is 15 • , and the accumulation time for TEC data used at one tomography case is 15 min.The ionospheric residual observations and dual-frequency pseudorange ob- servations (Melbourne-Wübbena combination) are used to detect the cycle slips and gross errors.A modified singlelayer ionospheric model and a spherical harmonic expansion function model are adopted to determine the differential code biases (Jin et al., 2012).Only the satellite-receiver rays that propagate through both the upper and lower boundary in the altitudinal direction of the reconstruction area are used here.Eight ground receivers located in the reconstruction area are chosen, and their geographical positions are displayed in Fig. 1 and are unevenly distributed.The phase and pseudorange data were downloaded from the IGS website.

Regularization method
The different reconstruction results can be obtained when different regularization constraints are incorporated.Here, the ionospheric spatial variations are assumed to be separable in the vertical and horizontal directions.In the vertical direction, the Gaussian correlation constraint is used, and the correlation distance is derived from the statistical results of Yue et al. (2007a).The correlation distance increases exponentially with altitude.It is about 20 km at the ionospheric E layer and F layer, and is about 500 km at the height of 2000 km.In the horizontal direction, two different regularization terms are imposed (denoted as regularization method I and regularization method II) as follows.

Regularization method I
The electron density in the horizontal direction is constrained by a multipoint finite-difference approximation of the twoorder Laplace operator (Hobiger et al., 2008).For a grid on a certain layer, the constraint operator is .The constraint operator for a grid at the northern, southern, western, and eastern boundary of the layer is , and the constraint operator for a grid at the northwest, northeast, southwest, and southeast corner of the layer is .Then the smoothness constraint matrix H can be constructed.To obtain the estimate of the ionospheric electron density x, the following cost functional should be minimized (Huang and Wu, 2011): where A is the observation matrix, its element a ij represents the length of the ith ray propagating through the j th grid, y is the column vector consisting of TEC measurements data, and R is the observation error covariance matrix, assuming that the measurement error is independent (Bust et al., 2004;Yue et al., 2007b) and is taken as 2 TECU, so it is a diagonal matrix.x b is the background electron density, B v is the background error covariance matrix in the vertical direction, and H is the smoothness constraint matrix.O is the corresponding covariance matrix of the smoothness constraint, assuming that it is a diagonal matrix and the value of its diagonal elements are the square of the average background electron density value in each layer.α and β are regularization parameters, and superscript T represents the matrix transpose operation.The dimension of matrix B v is very large, and its inverse operation is time-consuming.To avoid calculating the inverse of this matrix, the Cholesky decomposition method is used; that is, 1) is rewritten as (2)

Regularization method II
The electron density constraint in the horizontal direction is taken as Gaussian correlation, and the correlation distance is also derived from the statistical results of Yue et al. (2007a); that is, the horizontal correlation distance is about 16 • .The cost functional of the regularization method II is The symbols in this equation are consistent with the formula mentioned above.B h is the background error covariance matrix in the horizontal direction.To avoid calculating the inverse of the large dimension matrix, the Cholesky decomposition method is adopted; that is, The above equation can be rewritten as

Regularization parameter selection
The linear model function method in the framework of the damped Morozov discrepancy principle is used to determine the regularization parameters (Wang, 2012), and its basic idea is constructing a linear function to locally approximate the original function at each iteration step, greatly reducing the calculation time.In the following, we use this method to determine the optimal regularization parameters in the cost functional (2) as an example.
According to the damped Morozov discrepancy principle, the cost functional (2) can be rewritten as where γ > 1 and κ > 1 are damping coefficients, c ≥ 1 is a constant, and δ is the error level.For fixed regularization parameters α and β, the minimum value of the cost functional ( 2) is denoted as Then its deviation equation is In the kth iteration step for solving this equation, the linear model function , where T k , C k , and D k are constants needed to solve in each iteration step, and then The flow chart of determining the regularization parameters is shown in Fig. 2. For certain α k and β k , the corresponding electron density  0) , and α k+1 can be obtained by α k+1 = ωα k , where ω ∈ (0, 1) is a fixed constant.When the iteration stop criterion, where ε 1 and ε 2 are constants, is satisfied, the derived α k+1 and β k+1 are the optimal regularization parameters, and the corresponding x is the reconstructed electron density; otherwise, the α k+1 and β k+1 are set as the initial values to repeat the above steps.In the following test, we take γ = 6.5, κ = 3.5, ε 1 = 10 −5 , ε 2 = 10 −8 , and ω = 0.618.

Ideal test
In the ideal test, the real position between GNSS satellites and ground receivers on 8 March 2013 is used to establish the observation matrix A. The background electron density and simulated true electron density are both provided by the IRI2012 model (Bilitza et al., 2014).In this model, the IG (ionospheric global) index influences the electron density peak height, and the Rz (sunspot number) index influences the electron density peak density.To make the background electron density have a large difference from the simulated true electron density, the input parameters for the background model are set to IG−20, Rz−20 (denoting background I) and IG+20, Rz+20 (denoting background II).The simulated TEC measurement is derived by the length of the satellitereceiver ray propagating through each grid multiplied by the corresponding electron density (i.e., forward problem).To avoid the "inverse crime" and yielding unrealistically optimistic results, the grid sizes in the forward and inverse problem are not the same (Kaipio and Somerrsalo, 2005).The spatial steps in the forward problem are 0.5 • in latitude, 0.5 • in longitude, and 20 km in altitude.Due to the fact that the real measurements are inevitably subject to observational errors, noises are artificially imposed on the simulated TEC observations.
Figure 3 shows the reconstruction results using regularization method I under the background I condition, and the background electron density is lower than the simulated true electron density.From top to bottom are the longitudelatitude slices at different altitudes (270 km, 290 km, 310 km, and 330 km), the altitude-latitude slices at different longi-  The average relative error is 12.89 %, and the standard deviation of the relative error is 8.61 %.Even though there are many satellite-receiver rays propagating through the certain grid, its reconstructed relative error is still large under some situations.The TEC measurement error caused by the different grid sizes in the forward and inverse problem are not the same along different satellite-receiver ray paths, leading to the large relative error in this region.
Figure 4 shows the average relative error of the reconstructed electron density (red dashed line) and standard deviation of the reconstructed relative error (blue dashed line) on 8 March, and the average relative error of the background electron density (red solid line) and standard deviation of the  background relative error (blue solid line) are also superimposed for reference.After using the regularization method I, the average relative reconstruction error is significantly reduced compared with that of the background model, but the standard deviation of the relative reconstruction error increases during some periods.
Figure 5 is the same as Fig. 3, but for the reconstruction results under background II conditions.The background electron density is larger than the simulated true electron density.The regularization parameters are 0.0034 and 0.1502, and the maximum absolute error of the reconstructed electron density is about 1.4×10 11 electrons m −3 and mainly occurs near the southeast boundary of the reconstruction area.The maximum relative error is 46.79 %, and the large relative error mainly occurs at the height of about 170 km in the center and northeast of the reconstruction area.The average relative reconstruction error is 14.83 %, and the standard deviation of the reconstructed relative error is 9.20 %. Figure 6 is the same as Fig. 4, but for the average relative error and standard deviation under background II conditions.The average relative error of the reconstructed electron density is less than that of the background electron density at these periods, but the standard deviation of the reconstructed relative error is sometimes greater than that of the background relative error.
The reconstruction results using regularization method II under the background I situation are shown in Fig. 7, and the  Figure 8.The comparisons between the average relative error and standard deviation of the reconstructed relative error using regularization method II and the average relative error and standard deviation of the relative error using the single-parameter regularization method under background I conditions.background electron density is lower than the simulated true electron density.The regularization parameters are 0.0012 and 5 × 10 −5 , and the maximum absolute error of the reconstructed electron density is about 5 × 10 10 electrons m −3 and mainly occurs near the center of the reconstruction area.The maximum relative error is 115.07 %, the large relative error mainly occurs at the height of about 170 km in the center and north of the reconstruction area, the average relative reconstruction error is 11.04 %, and the standard deviation of the reconstructed relative error is 10.33 %.
Figure 8 shows the average relative error of the reconstructed electron density (red dashed line) and standard deviation of the reconstructed relative error (blue dashed line) on 8 March.For comparison, the average relative error (red solid line) and standard deviation of the relative error (blue solid line) using the single-parameter regularization method are also shown.The cost functional of the single-parameter regularization method should be minimized to derive the approximate solution of the ionospheric electron density: where B is the background error covariance matrix, constructed in the same way as in Bust et al. (2004).The Gaussian correlation and the Gaussian elliptic correlation are used www.ann-geophys.net/36/1255/2018/Ann.Geophys., 36, 1255-1266, 2018  Figure 10.The comparisons between the average relative error and standard deviation of the relative error using regularization method II and the average relative error and standard deviation of the relative error using the single-parameter regularization method under background II conditions. in the vertical and horizontal direction, respectively, and the correlation distances are taken from Yue et al. (2007a).The average relative error and standard deviation of the reconstructed electron density using the dual-parameter regularization method are generally smaller than those of using singleparameter regularization method.
Figure 9 is the same as Fig. 7, but for the reconstruction results under the background II conditions.The regularization parameters are 0.0012 and 5.14 × 10 −4 , and the maximum absolute error of the reconstructed electron density is about 6×10 10 electrons m −3 and mainly occurs near the center of the reconstruction area.The maximum relative error is 55.32 %, the large relative error occurs mainly at the height of about 170 km in the center and northwest of the reconstruction area, the average relative error is 12.21 %, and the standard deviation of relative error is 8.93 %. Figure 10 is the same as Fig. 8, but for the reconstruction results under the background II condition, the reconstructed electron density using the dual-parameter regularization method show further improvement compared with that using the single-parameter regularization method.
Regularization method II can transfer the observation information in one place to a nearby place by the action of the background error covariance.The regularization method I transfers the observation information under the action of the Laplace operator, and its influence area is relatively limited.Overall, when the observation data are sparse, the regularization method II has a better performance.In the following, only the reconstruction results using regularization method II are shown in the real measurements test.

Real observation
In the real measurements cases, the effectiveness of the regularization method II is tested under the quiet and active geomagnetic activity conditions.The ionosonde measurement from MHJ station is used as the independent validation data.These data were obtained from GIRO and were manually scaled by the SAO software, and they can provide accurate electron peak density and peak height.Due to the topside ionospheric electron density profiles derived by the extrapolation method, the electron density data below the height of 500 km are used to validate the reconstruction result.In the real reconstruction, it is assumed that the electron density in each grid does not vary within a 15 min interval.
The geomagnetic activity on 8 March was relatively quiet.The reconstruction results are shown in Fig. 11, in which the red line is the reconstructed electron density profile, the blue line is the background electron density profile, and the green line is the observed electron density profile from ground ionosonde.There are no ionosonde data at 01:00 and 12:00 UT.Overall, the reconstructed peak electron den-sity using the regularization method II is much closer to the ionosonde measurements compared with that from the background model, and the reconstructed electron density peak height is basically similar to the background value.This is mainly because the satellite-receiver geometry has limitations, and the TEC data contain much more ionospheric structure information in the horizontal direction than in the vertical direction.Only GNSS TEC data have limited influences on the changes in altitudinal resolution in the background model.
On 17 March the geomagnetic activity was relatively active, and the spatial correlation distance used here is the same as that under the quiet geomagnetic activity condition.The reconstruction results are shown in Fig. 12, and there is no validation data from the MHJ ionosonde at 06:00 UT.Compared with the background model, the regularization method can distinctly improve the accuracy of the electron peak density, but these measurements are still quite different from the ionosonde measurements, which may be related to the accuracy of background error covariance.Due to the limitations of spatial and temporal resolution of ionospheric observation data and the complexity of the ionospheric variations, the background error covariance is not totally known and needs to be studied further, especially under the ionospheric disturbed condition.Moreover, the observation error is artificially assumed in the present paper, and so the real spatial and temporal variations of the measurement error also need some more researches.

Conclusions
Ionospheric tomography based on GNSS TEC is a typical ill-posed inverse problem, and the regularization method is used to solve this problem by incorporating prior constraints to approximate the real electron density distribution.Because of the complexity of the spatial and temporal variations of the ionosphere, a single regularization term may not obtain the high-accuracy reconstruction results.Multiple constraints sometimes need to be incorporated, and multiple regularization parameters are introduced accordingly.Here, the ionospheric variations are separated in the horizontal and vertical direction.The cost functional of the dual-parameter regularization method is established, and the regularization parameters are determined by the linear model function method.This regularization algorithm is tested by the ideal cases and real cases, and the results show that it can significantly improve the background model outputs.Moreover, to improve the reconstruction accuracy further, the measurement error and the background error covariance should be studied in detail.

Figure 2 .
Figure 2. The flow chart of determining the regularization parameters by using the linear model function method.

Figure 3 .
Figure 3.The reconstruction results using regularization method I under background I conditions.From top to bottom are the longitudelatitude slices at different altitudes, altitude-latitude slices at different longitudes, and longitude-altitude slices at different latitudes.From left to right are the background electron density, true electron density, reconstructed electron density, reconstructed absolute error, reconstructed relative error, and whether the satellite-receiver ray propagates through the corresponding grid or not.

Figure 4 .
Figure 4.The average relative reconstruction error and standard deviation using regularization method I under the background I situation.The red solid line is the average relative error of the background electron density.The blue solid line is the standard deviation of the background relative error.The red dashed line is the average relative error of the reconstructed electron density.The blue dashed line is the standard deviation of reconstructed relative error.

Figure 5 .
Figure 5.The reconstruction results using regularization method I under background II conditions.From top to bottom are the longitudelatitude slices at different altitudes, altitude-latitude slices at different longitudes, and longitude-altitude slices at different latitudes.From left to right are the background electron density, true electron density, reconstructed electron density, reconstructed absolute error, reconstructed relative error, and whether the satellite-receiver ray propagates through the corresponding grid or not.

Figure 6 .
Figure6.The average relative reconstruction error and standard deviation using regularization method I under the background II situation.The red solid line is the average relative error of the background electron density.The blue solid line is the standard deviation of the background relative error.The red dashed line is the average relative error of the reconstructed electron density.The blue dashed line is the standard deviation of the reconstructed relative error.

Figure 7 .
Figure 7.The reconstruction results using regularization method II under background I conditions.From top to bottom are the latitudelongitude slices at different altitudes, altitude-latitude slices at different longitudes, and longitude-altitude slices at different longitudes.From left to right are the background electron density, true electron density, reconstructed electron density, reconstructed absolute error, reconstructed relative error, and whether the satellite-receiver ray propagates through the corresponding grid or not.

Figure 9 .
Figure9.The reconstruction results using regularization method II under background II conditions.From top to bottom are the longitudelatitude slices at different altitudes, altitude-latitude slices at different longitudes, and longitude-altitude slices at different latitudes.From left to right are the background electron density, true electron density, reconstructed electron density, reconstructed absolute error, reconstructed relative error, and whether the satellite-receiver ray propagates through the corresponding grid or not.

Figure 11 .
Figure 11.The comparisons between the reconstructed electron density profiles near the MHJ station and the measured electron density profiles from ground ionosonde under the relatively quiet geomagnetic activity condition; LT = UT−4.6h.

Figure 12 .
Figure12.The comparisons between the reconstructed electron density profiles near the MHJ station and the measured electron density profiles from ground ionosonde under the relatively active geomagnetic activity condition; LT = UT−4.6h.