Adaptive Smoothness Constraint Ionospheric Tomography Algorithm.

Ionospheric tomography reconstruction based on global navigation satellite system observations is usually an ill-posed problem. To resolve it, an adaptive smoothness constraint ionospheric tomography algorithm is proposed in this work. The new algorithm performs an adaptive adjustment for the constrained weight coefficients of the tomography system. The computational efficiency and the reconstructed quality of ionospheric imaging are improved by using the new algorithm. A numerical simulation experiment was conducted in order to validate the feasibility and superiority of the algorithm. The statistical results of the reconstructed errors and the comparisons of ionospheric profiles confirmed the superiority of the new algorithm. Finally, the new algorithm was successfully applied to reconstruct three-dimensional ionospheric images under geomagnetic quiet and geomagnetic disturbance conditions over Hunan province. The tomographic results are reasonable and consistent with the general behavior of the ionosphere. The positive and negative phase storm effects are found during geomagnetic storm occurrence.


Introduction
The ionosphere is the upper atmosphere between 60 km and 1000 km above the earth. The temporal and spatial variations of the ionosphere affect the propagation of radio signals [1][2][3], especially during the occurrence of geomagnetic storms. Therefore, it is necessary to detect and understand the characteristics of ionosphere activity under different geomagnetic conditions. Ionospheric electron density (IED) and vertical total electron content (VTEC) are the key physical parameters for the studies of the ionospheric activity [4]. However, VTEC models can only reconstruct the two-dimensional ionospheric structure in the cross-section of latitude and longitude. In late 1980s, Austen et al. [5] introduced the computerized ionospheric tomography technique to reconstruct the IED distributions by using the observations of the navy navigation satellite system. However, the longitudinal variation is ignored because the ground receivers are aligned along the fixed longitudinal chain. To reconstruct the ionospheric structure, three-dimensional ionospheric variations need to be studied. Using simulated high-orbit satellite data, Kunitsyn et al. [6] successfully reconstructed the three-dimensional ionospheric structure. The advent and development of the global navigation satellite system (GNSS) has provided a new avenue for monitoring and studying ionosphere. The construction of more ground-based GNSS networks and the application of more space-based GNSS observations make it possible to reconstruct the three-dimensional IED distributions. Therefore, the GNSS-based ionospheric tomography technique has attracted the interest of many scholars [7][8][9][10][11][12][13][14][15][16]. Theoretical and experimental research has been conducted, and some meaningful results have been achieved. Leitinger et al. [17]

Tomographic Theory
The slant total electron content (STEC) along the ray path can be expressed as: STEC = p N e (s)ds (1) where N e (s) is the IED distribution, and p represents the ray path between a satellite and a receiver.
To reconstruct the three-dimensional IED distribution using the STEC data, the selected region is divided into small cubic voxels. The ionosphere is discretized into i layers, j layers, and k layers in longitudinal, latitudinal, and altitudinal directions, respectively. The number sequences of the voxels are shown in Figure 1. Then, Equation (1) can be simplified as: Equation (2) is generally represented as a matrix notation: where m represents the number of the input STEC data; n is the total voxel numbers, which equals i jk; y is the vector of the m known STEC measurements; A is a coefficient matrix, the element A cd is the length of the cth path in the dth voxel; x is a vector of the unknown IED; and e is a column vector related to the discretization errors and measurement noises [31]. is the length of the th path in the th voxel; is a vector of the unknown IED; and is a column vector related to the discretization errors and measurement noises [31].

Adaptive Smoothness Constraint MART
To solve the ill-posed problem of ionospheric tomographic inversion, it is necessary to select a reasonable algorithm. Compared with the conventional MART algorithm, the SCMART is very attractive, since it overcomes the disadvantages of the MART to some extent. In general, the MART algorithm is iterated cyclically and can be implemented as follows [25]: where ( ) refers to the ℎ member of the vector after the th round iteration; is the th row of the vector ; is the th row of the matrix ; and is the relaxation parameter, which is set to be 0.2 in this work.
In general, the IED variation is smooth and continuous, and the correlation of IED between voxels increases with the distance shortening in the same horizontal layer. The Gauss weighted function is used to create horizontal constraints in the SCMART [30]. For each horizontal layer, we construct a constraint matrix, which is given as

Adaptive Smoothness Constraint MART
To solve the ill-posed problem of ionospheric tomographic inversion, it is necessary to select a reasonable algorithm. Compared with the conventional MART algorithm, the SCMART is very attractive, since it overcomes the disadvantages of the MART to some extent. In general, the MART algorithm is iterated cyclically and can be implemented as follows [25]: where x (p) d refers to the dth member of the vector x after the pth round iteration; y c is the cth row of the vector y; A c is the cth row of the matrix A; and λ is the relaxation parameter, which is set to be 0.2 in this work.
In general, the IED variation is smooth and continuous, and the correlation of IED between voxels increases with the distance shortening in the same horizontal layer. The Gauss weighted function is Sensors 2020, 20, 2404 4 of 16 used to create horizontal constraints in the SCMART [30]. For each horizontal layer, we construct a constraint matrix, which is given as where q is the number of voxels in one horizontal layer; and L refers to the L th layer in the altitudinal direction. In the SMART, the element of H L is calculated using the following equation: where D cd is the distance between the cth voxel and the dth voxel; and σ is the smoothing operator.
Since the distance D cd among voxels is stable, we can see that the h cd is invariant in each round of iteration. It influences the imaging quality of IED and the computational efficiency. Taking into account the change of IED after each round of iteration, the h cd can be adaptively adjusted according to the results of the last round of iteration in the ASMART, and then Equation (7) is changed as In Equation (8), Q is the adaptive adjustment factor, which can be describe as Considering all the horizontal layers, the horizontal constraint can be represented as a matrix notation: According to Equation (5), the horizontal constraint matrix is written as In the altitudinal direction, the variation of IED is usually smooth and continuous, and the adaptive constraint matrix G can then be created using Equation (12). There are i j(k − 1) lines and i jk columns in matrix G.
In the altitudinal direction, the variation of IED is usually smooth and continuous, and the adaptive constraint matrix can then be created using Equation (12). There are ( − 1) lines and columns in matrix .
( ) ( ) Therefore, the altitudinal constraint can be represented as a matrix notation: (12) Therefore, the altitudinal constraint can be represented as a matrix notation: Combining Equation (3) with Equations (10) and (13), we can obtain Equation (14) is expressed as

Algorithm Validation
The ASCMART is an improvement on the SCMART. However, the SCMART is an improvement of the MART. Instead of the AMART, the MART is introduced to compare with the ASCMART and the SCMART in this work. To confirm the feasibility and superiority of the new algorithm over the MART and the SCMART methods, a numerical simulation scheme was devised. The simulation scenario is as follows: (1) The reconstructed region was Hunan province of southern China (24.
, and the altitude ranged from 100 km to 1000 km. The discretized intervals are 0.5 • in the longitudinal and latitudinal directions. Different from the equal spacing division in the horizontal direction, the unequal spacing division was adopted in the altitudinal direction. The discretized interval was 20 km in the altitudinal range from 200 km to 400 km. However, the discretized interval was 50 km in other altitudinal ranges.
(2) In the numerical simulation experiment, the selected time period was from 00:00 UT to 00:30 UT on 20 June 2015.
(3) The coordinates of the observed satellites and the 124 CORS stations in Hunan Province were used to calculate the crossing distance between the rays and the corresponding voxels, and the coefficient matrix was then constructed. The distributions of the 124 CORS GNSS stations and an ionosonde station were shown in Figure 2.
was 50 km in other altitudinal ranges.
2) In the numerical simulation experiment, the selected time period was from 00:00 UT to 00:30 UT on June 20, 2015.
3) The coordinates of the observed satellites and the 124 CORS stations in Hunan Province were used to calculate the crossing distance between the rays and the corresponding voxels, and the coefficient matrix was then constructed. The distributions of the 124 CORS GNSS stations and an ionosonde station were shown in Figure 2.  To mimic Y simu , the true IED distributions were generated from IRI 2016 model. Therefore, the STEC values without noise can be computed using the following formula: A small amount of random noise E simu , which satisfied the Gaussian distribution, was added to the simulated Y simu in order to obtain more realistic values Y noise .
To facilitate the comparison of the ASCMART method with the SCMART and MART methods, the initial values are needed to operate the three algorithms. In this work, a different model such as NeQuick was selected as the background.
In order to evaluate the accuracy of the tomographic results of the three algorithms, the average density (AVD) error and root mean square (RMS) error were used in this work. They can be written as Using the above discretized scenario, the reconstructed space was divided into 3456 voxels. In the numerical simulation, the IRI 2016 model was applied to simulate the true value of all voxels. At the same time, the initial value of each voxel was obtained from NeQuick model in order to perform the algorithms and distinguish the true value. Figure 3 compares the simulated result of the IRI 2016 model with the tomographic results of the MART, SCMART and ASCMART methods. From Figure 3c, we can see that the tomographic result of the ASCMART algorithm is in good agreement with the true value. This suggests that the ASCMART can be used to reconstruct the IED distribution. Comparing Figure 3c with Figure 3a,b, it can be seen that the result of the ASCMART is more consistent with the true value than those of the SCMART and the MART. the same time, the initial value of each voxel was obtained from NeQuick model in order to perform the algorithms and distinguish the true value. Figure 3 compares the simulated result of the IRI 2016 model with the tomographic results of the MART, SCMART and ASCMART methods. From Figure  3c, we can see that the tomographic result of the ASCMART algorithm is in good agreement with the true value. This suggests that the ASCMART can be used to reconstruct the IED distribution. Comparing Figure 3c with Figure 3a,b, it can be seen that the result of the ASCMART is more consistent with the true value than those of the SCMART and the MART.  Table 1 provides the reconstructed error statistics of the three algorithms based on all voxels. The RMS error of the ASCMART is 7 × 10 9 el/m 3 , which is smaller than those of the MART and SCMART.   Table 1 provides the reconstructed error statistics of the three algorithms based on all voxels. The RMS error of the ASCMART is 7 × 10 9 el/m 3 , which is smaller than those of the MART and SCMART.  Figure 4 illustrates the reconstructed error contour of the three algorithms along different longitudinal chains. From Figure 4, it can be seen that the reconstructed error of the ASCMART is the smallest of the three algorithms. Figure 4c and Table 1 show that the tomographic accuracy of the ASCMART method is the highest of the three methods, which confirms the reliability and superiority of the ASCMART.  Figure 4, it can be seen that the reconstructed error of the ASCMART is the smallest of the three algorithms. Figure 4c and Table 1 show that the tomographic accuracy of the ASCMART method is the highest of the three methods, which confirms the reliability and superiority of the ASCMART. To evaluate the reconstructed efficiency of the three algorithms, the iteration convergence condition was defined as follows in this work: where and represent the true value vector and the final tomographic solution vector of IED distributions, respectively. To evaluate the reconstructed efficiency of the three algorithms, the iteration convergence condition was defined as follows in this work: where x true and x tomo represent the true value vector and the final tomographic solution vector of IED distributions, respectively. In the simulation, the convergence of the iteration was terminated when S < 10 −3 . Since the AMART improved the convergence speed of the MART, the AMART was added to compare with the three algorithms. Figure 5 illustrates the convergence curves of the four algorithms. From Figure 5, we can see that the convergence speed of the ASCMART algorithm is faster than those of the SCMART, AMART, and MART. This means that the computational efficiency of the ASMART is the fastest in the four algorithms.     To test the characteristics of the ASCMART by using actual observations, the observations of 124 CORS stations in Hunan province were introduced in this work, and the sample interval was 30 s. The selected time period was 20 June and 22-23 June 2015. Using the STEC derived from the selected dual-frequency GNSS data, the ASCMART method was applied to reconstruct the IED distributions in Hunan province. The IED temporal-spatial distributions under the two different ionospheric time periods were studied.

Tomographic Reconstruction of IED Distribution Using Actual GNSS Data
The STEC is the input data of the IED tomographic distributions. The reconstructed accuracy of IED depends on the STEC accuracy. Using the differential phases (STEC Φ ) and the differential pseudoranges (STEC p ), an absolute STEC with high accuracy may be obtained by introducing an extra term B L . It can be formulated as follows [32]: Sensors 2020, 20, 2404 9 of 16 If N measurements are obtained during a satellite pass, B L can be modeled using the following equation: Figure 5. Comparison of the convergence speed of ASMART, SCMART, AMART, and MART.  To test the characteristics of the ASCMART by using actual observations, the observations of 124 CORS stations in Hunan province were introduced in this work, and the sample interval was 30 s. The selected time period was June 20 and June 22-23, 2015. Using the STEC derived from the selected dual-frequency GNSS data, the ASCMART method was applied to reconstruct the IED distributions in Hunan province. The IED temporal-spatial distributions under the two different ionospheric time periods were studied.

Tomographic Reconstruction of IED Distribution Using Actual GNSS Data
The STEC is the input data of the IED tomographic distributions. The reconstructed accuracy of IED depends on the STEC accuracy. Using the differential phases ( ) and the differential To obtain high-accuracy STEC, the preprocessing of the GNSS data is very important. The method presented by Blewitt [33] was employed to clean outliers and repair cycle slips in GNSS observations. The differential phases was formed to eliminate the clock bias and the tropospheric error. In general, the instrumental biases of satellites and receivers were usually fixed in a single day. In this work, the instrumental bias was fitted using the least square technique.
Using the STEC data derived from the selected GNSS observations, the ASCMART method was applied to reconstruct the three-dimensional IED distributions on geomagnetic quiet and geomagnetic storm days. The geomagnetic activity was very quiet on 20 June 2015. Figure 7 illustrates the three-dimensional IED distributions on the day. From Figure 7, it can be seen that the peak height of IED rises from 280 km to 360 km between 00:00 UT and 04:00 UT, and the peak height then falls to 320 km from 08:00 UT to 16:00 UT. At 20:00 UT, the peak height of IED falls to 280 km. The reconstructed images reflect the variations of the ionospheric vertical structure. Meanwhile, Figure 7 shows that the IED values first increase from 00:00 UT to 08:00 UT. As time goes on, the IED values gradually fall. At 20:00 UT, the IED values are the minimum. This is consistent with the normal change laws in daytime and nighttime over Hunan province, as well as the fact that the IED variation depends mainly on the intensity of solar radiation. In addition, we can see that the values of IED in the north of Hunan province are smaller than those in the south as a whole. This suggests that the IED distributions has strong correlation with latitude.
To further verify the reliability of the tomographic results of the ASCMART, the ionospheric peak density hmF2 and the peak height NmF2 obtained from the ionosonde were compared with those of the ASCMART algorithm. In this work, we selected the ionosonde data from Shaoyang station, whose geographic location is shown in Figure 2. The comparison is shown in Figure 8. The comparisons confirm that the hmF2 and NmF2 extracted from the ASCMART algorithm have a good agreement with those obtained from the ionosonde as a whole. This suggests that the reconstructed results of the ASCMART are reliable.  To evaluate the superiority of the ASCMART, the IED reconstructed error can be obtained from the difference between ionosonde data and the results of the three algorithms at the geographic location of Shaoyang station. Table 2 shows the error statistics of the ASCMART, SCMART, and MART. From Table 2, it can be seen that the AVD error and the RMS of the ASCMART are the smallest among the three algorithms. This validates the ASCMART's superiority compared to the other two algorithms. To evaluate the superiority of the ASCMART, the IED reconstructed error can be obtained from the difference between ionosonde data and the results of the three algorithms at the geographic location of Shaoyang station. Table 2 shows the error statistics of the ASCMART, SCMART, and MART. From Table 2, it can be seen that the AVD error and the RMS of the ASCMART are the smallest among the three algorithms. This validates the ASCMART's superiority compared to the other two algorithms. As mentioned above, a strong magnetic storm occurred between 15:00 UT on June 22 and 12:00 UT on June 23, 2015. We reconstructed the three-dimensional images of the IED distribution over Hunan province by using the ASCMART method and the GNSS data in this period. Figure 9 shows the time-series variation of the ionospheric structure on the selected geomagnetic storm days.  As mentioned above, a strong magnetic storm occurred between 15:00 UT on 22 June and 12:00 UT on 23 June 2015. We reconstructed the three-dimensional images of the IED distribution over Hunan province by using the ASCMART method and the GNSS data in this period. Figure 9 shows the time-series variation of the ionospheric structure on the selected geomagnetic storm days. Figure 10 illustrates the difference between the storm-time IED and the reference values, which are equal to the 15-day average before the storm occurrence. From Figure 10, it can be seen that the positive phase storm appeared between 260 km and 360 km at 16:00 UT on 22 June 2015. After that, the IED was almost unchanged between 20:00 UT on 22 June and 00:00 UT on 23 June. As time went on, the positive and the negative phase storm varied with time at different altitude. At 04:00 UT on 23 June, the strong positive phase storm appeared at altitudes of 200 km, 240 km, 280 km, and 320 km in the northwest of the reconstructed region, where the IED values increased obviously. At 08:00 UT on 23 June, the strong negative phase storm occurred at the altitude of 280 km in the north of Hunan province. The IED values increased slightly at the altitude of 240 km and 280 km but decreased at the altitude of 320 km and 360 km in the southwest of the reconstructed region. Afterwards, the IED showed almost no change except for a small-area increase at the altitude of 280 km and 320 km. This fact suggests that the reconstructed results of the ASCMART can reflect the storm phase variations with the altitude.
To further characterize the temporal variations of the total electron content (TEC) maps during the storm, the TEC values were obtained by using the GNSS observations in Hunan province during the storm occurrence, and we calculated the differences between the above TEC values and the 15-day average TEC values before the storm occurrence. Figure 11 presents the time-series variations of the differential TEC values in the reconstructed region on 22 and 23 June 2015. Figure 11 indicates that the differential TEC values showed a small increase between 16:00 UT and 18:00 UT on 22 June. As time went on, the differential TEC values increased significantly from 02:00 UT to 06:00 UT on 23 June. However, the differential TEC values decreased greatly from 08:00 UT to 10:00 UT on 23 June and subsequently recovered to the normal state. The variation of the differential TEC values is generally in agreement with that of the IED change.  Figure 10 illustrates the difference between the storm-time IED and the reference values, which are equal to the 15-day average before the storm occurrence. From Figure 10, it can be seen that the positive phase storm appeared between 260 km and 360 km at 16:00 UT on June 22, 2015. After that, the IED was almost unchanged between 20:00 UT on June 22 and 00:00 UT on June 23. As time went UT on June 23, the strong negative phase storm occurred at the altitude of 280 km in the north of Hunan province. The IED values increased slightly at the altitude of 240 km and 280 km but decreased at the altitude of 320 km and 360 km in the southwest of the reconstructed region. Afterwards, the IED showed almost no change except for a small-area increase at the altitude of 280 km and 320 km. This fact suggests that the reconstructed results of the ASCMART can reflect the storm phase variations with the altitude. To validate the reliability of the ASCMART during geomagnetic storms occurrence, the comparison of hmF2 and NmF2 was made between the ionosonde and the ASCMART. Figure 12 shows the comparison results; it can be seen that the hmF2 and NmF2 of the ASCMART have good agreement with those obtained from the ionosonde. This validates the feasibility and the reliability of the ASCMART to reconstruct three-dimensional images of the ionosphere under geomagnetic disturbance. day average TEC values before the storm occurrence. Figure 11 presents the time-series variations of the differential TEC values in the reconstructed region on June 22 and 23, 2015. Figure 11 indicates that the differential TEC values showed a small increase between 16:00 UT and 18:00 UT on June 22. As time went on, the differential TEC values increased significantly from 02:00 UT to 06:00 UT on June 23. However, the differential TEC values decreased greatly from 08:00 UT to 10:00 UT on June 23 and subsequently recovered to the normal state. The variation of the differential TEC values is generally in agreement with that of the IED change.  To validate the reliability of the ASCMART during geomagnetic storms occurrence, the comparison of hmF2 and NmF2 was made between the ionosonde and the ASCMART. Figure 12 shows the comparison results; it can be seen that the hmF2 and NmF2 of the ASCMART have good agreement with those obtained from the ionosonde. This validates the feasibility and the reliability of the ASCMART to reconstruct three-dimensional images of the ionosphere under geomagnetic disturbance.  Table 3 shows the error statistics of the four algorithms. Accordingly, we can see that the reconstructed accuracy of the ASCMART is higher than that of the SCMART and MART algorithms. This confirms that the ASCMART has superiority over the SCMART and MART.  Table 3 shows the error statistics of the four algorithms. Accordingly, we can see that the reconstructed accuracy of the ASCMART is higher than that of the SCMART and MART algorithms. This confirms that the ASCMART has superiority over the SCMART and MART.

Conclusions
In order to resolve the ill-posed problem of GNSS-based ionospheric tomography and the limitation of the conventional MART and the SCMART, a new ionospheric algorithm is presented in this work. The algorithm surmounts the dependence on the initial values for the voxels that have rays traversing them by imposing a horizontal and vertical smoothness constraint on the tomographic system, which improves the reconstructed accuracy of the IED distributions. Furthermore, the computational efficiency is improved by adaptively adjusting the constraint weight coefficients. The reliability of the GNSS-based tomographic algorithm was further tested by using GNSS data under geomagnetic quiet and geomagnetic disturbance conditions. The time-series variations of the three-dimensional ionospheric structure were studied, and the results were reasonable and consistent with the general behavior of the ionosphere. The positive and negative phase storm effects were found during geomagnetic storm occurrence. Comparison with profiles obtained from ionosonde data showed good agreement, and the validity of the tomography results were confirmed.