Complementing regional ground GNSS-STEC computerized ionospheric tomography (CIT) with ionosonde data assimilation

A near-real-time computerized ionospheric tomography (CIT) technique was developed over the East Asian sector to specify the 3-D electron density field. The technique is based on a plethora of Global Navigation Satellite System observables within the region of interest which is bounded horizontally 110°–160°E and 10°–60°N and extending from 80 to 25,000 km in altitude. Prior to deployment, studies validated the CIT results using ionosonde, middle-upper atmosphere radar and occultation data and found the technique to adequately reconstruct the regional ionosphere vertical structure. However, with room for improvement in estimating the peak height and avoiding physically unrealistic negative densities in the final solution, we present preliminary results from a technique that addresses these issues by incorporating CIT results into a data assimilation (DA) technique. The DA technique adds ionosonde bottomside measurements into CIT results, thereby improving the accuracy of the reconstructed bottomside 3-D structure. More specifically, on average CIT NmF2 and hmF2 improve by more than 60%. Further, during analysis, ionospheric electron densities are assumed to be better described by probability log-normal distribution, which introduces the positivity constraint that is mandatory in ionospheric imaging.


Introduction
Across different industry fields, operation and decision making are increasingly based on the precision/accuracy of trans-ionospheric signals. Subsequently, the spatiotemporal high-resolution monitoring and forecasting of the state of the ionosphere is vital. Hitherto, the abundant availability of global navigation satellite system (GNSS) receivers offers a cost-effective way of monitoring the ionosphere in near real time. More so, tomography techniques based on GNSS observables have gained traction since they provide an insight into the three-dimensional (3-D) structure of the ionosphere. These techniques are predominantly suitable for regional ionosphere imaging, wherein high-resolution spatial grids and computation expense are tractable (Seemala et al. 2014;Saito et al. 2017).
Over the East Asian sector, a near-real-time regional 3-D tomography technique, with a latency of about 6 min, has been developed to cover a region 10° N-60° N in latitude, 110° E-160° E in longitude, and 80-25,000 km in altitude. To reduce the computation cost, the voxel grid is set up in such a way that the horizontal spatial resolution is highest over the area with the most GNSS receivers. That is, 1° × 1° for the area 129°-140° E and 33°-39° N; 2° × 2° for the area enclosing the high-resolution region, extending 125°-150° E and 25°-50° N; otherwise, 5° × 5°. Vertically, the highest resolution is set over the ionospheric E-and F-regions where large density variations are expected, i.e., 20, 50, 100, 3000, and 5000 km, for the altitude ranges of 80-600, 600-900, 900-2000, 2000-5000, and 5000-25,000 km, respectively. For a pictorial illustration of the 3-D grid setup, refer to Saito et al. (2017). Saito et al. (2017Saito et al. ( , 2019 analyzed reconstructions from the near-real-time tomography and found that although the technique on average performs well, there was an under-performance in estimating peak parameters, NmF2 and hmF2, that tends to occur below a height range of 270 km in the months of October and November as shown in Fig. 1. In addition, when solving the linear inverse problem, the current tomography technique lacks a positivity constraint such that the final reconstructed 3-D picture is sometimes prone to meaningless negative electron densities. As such, we report on the steps undertaken to address these issues by complimenting the original tomography technique with a data assimilation (DA) procedure, which facilitates the inclusion of ionosonde measurements into the reconstructed 3-D structure. The results presented here will be a demonstration of the capabilities of the newly modified technique rather than the final operational version. We will neither thoroughly review the building blocks of the tomography algorithm nor will we detail previous results. Instead, we encourage the reader to refer to Seemala et al. (2014) and Saito et al. (2017) for the full genesis of the algorithm.
The next section briefly discusses the core of the current real-time operational tomography technique. We then discuss in detail the steps taken to improve the fidelity and present results analyzing the performance. The final section gives concluding remarks and prospects.

Original tomography
In order to sketch a computerized ionospheric tomography (CIT) problem, assume the regional ionosphere to be partitioned into k finite elements, called voxels, � ⃗ n = (n 1 , n 2 , n 3 , ⋯ , n i , ⋯ , n k ) T , each with a uniform electron density distribution (n i ) . Superscript ("T") represents vector transpose. For a datum, here referring to a single observation,y l , that has a linear relationship with � ⃗ n , the goal is to determine a representation of � ⃗ n that minimizes the error l and satisfy the linear constraint Fig. 1 A comparison between an extract of peak parameters from the East Asia real-time 3-D tomography (Tomo) and observations from the middle-upper atmosphere radar (MUR) located 34.85° N, 136.10° E, Shigaraki, Japan. Top and bottom panels represent NmF2 and hmF2, respectively. On average, tomography results track the observed values relatively well, but with improvement possible during the months of August-October (NmF2) and November-January (hmF2) (Saito et al. 2019) where a i are the coefficients that determine the fractional contribution of voxel i to y l . In cases when a column vector ( � ⃗ Y ) of m observations is available at time t, a system of linear equations (SLE) can be formulated and represented in a matrix compact form Equation (2) is a typical linear inverse problem in which a least-squares estimator needs to be selected to minimize the uncertainty. However, in ionosphere linear inverse problems, in which only GNSS-ground receiver (Rx) links are utilized, the geometry constraints limit the inclusion of any horizontal information. In addition, GNSS satellites orbit at relatively low angular velocities, demands that the time window for the collection of measurements should be relatively small to avoid the averaging of ionospheric dynamics. In return, these limitations exacerbate under-determinedness and ill-posedness of (2) (Yeh and Raymund 1991;Raymund et al. 1994). To adequately solve (2), a prior information from other sources such as models is utilized. The accuracy of the prior information and the extent to which a technique depends on such information highly influences fidelity. In the current tomography technique, the dependence on prior information is limited through the use of a constrained leastsquares fit. That is, a weight matrix, ⃗ n, based a prior information restrains the derived electron density from exceeding a certain value. In addition, ⃗ n asserts a continuity condition in which a voxel i is coupled to its neighboring six voxels (q: up, down, north, south, east, and west). C iq determines the coupling "strength" (Seemala et al. 2014). To stabilize the solution, top and bottom boundary electron densities superscripted by bc in (4), are fixed to background model values.
The error function to be minimized is then written as and the solution (⃗ n) that minimizes (4) is expressed as where λ(> 0) is a regularization parameter that needs to be learned from experience. Note that the last term in (5) has a different sign compared to that presented in Saito et al. (2017), which was a derivation error. Nonetheless, because the boundary densities are very small, we found out that this correction does not significantly change the overall 3-D reconstructions.

Improving the vertical profile structure and dealing with negative densities
The fidelity of (5) is subject to the choice of constraint parameters C iq . These are dynamically derived from the NeQuick model (Nava et al. 2008) together with ad hoc coefficients that are "hard-wired" into the algorithm (Seemala et al. 2014). As mentioned earlier, Saito et al. (2017; found that the best agreement between observations and tomography results was recorded from May to July, but with more discrepancies in August-October and November-January for NmF2 and hmF2, receptively, as illustrated in Fig. 1. Better performance in May-July was not surprising since ad hoc coefficients and a set of regularization parameters used in the current tomography were designed or fine-tuned based on data covering the same period. Therefore, to improve the performance over other time periods, it was quite natural to first attempt at learning and adapting a new set of suitable ad hoc coefficients and λ, which could be time-dependent in terms of different months or seasons of the year. Unfortunately, different attempts at this approach resulted in a final solution that was largely unstable for imperceptible changes in the ad hoc coefficients. The results were not good and are not presented here. In addition, there is no literature we could find relating to the step by step procedure on how the original ad hoc coefficients in Seemala et al. (2014) were or should be computed. For these reasons, we opted for a different approach of using ionosonde data as an observed quantity to improve the vertical structure, focusing on NmF2 and hmF2 estimations.

Inclusion of ionosonde data
The first step at including ionosonde data into the algorithm was to design a geometry matrix ( iono ) that would facilitate the process of concatenating ionosonde data to the already existing GNSS (STEC) data. For example, an exaggerated cartoon sketch in Fig. 2 illustrates two ionosonde data points, which are to be included in the analysis, sampled along a vertical profile (left panel). The points mimic the expected case of either a datum lying inside (orange) or outside (blue) the grid. Solid black lines of different lengths d i , connect the ionosonde data point (reference point) to different grid points ( i ) that lie within a specified radius of influence r. Here, r was assumed to be 10°, a value based on studies that previously utilized mid-latitude horizontal correlation lengths (Bust et al. 2004;Ssessanga et al. 2019).
When a datum is located within the grid, its influence is limited to the two horizontal planes (above and below) bounding that datum point. Conversely, for a datum located outside the grid, its influence is limited to the first horizontal plane intersected by the solid lines. Using the assumption that grid points that are close to the ionosonde point should be affected more than those far away, we defined iono entries ( a i iono ) based on the inverse distance weighting criteria,w i , with constraint that enforces a distribution k ∑ i=1 w i = 1 (Bartier and Keller 1996;Ping et al. 2004); To build the geometry matrix, , that is needed in the inversion process, iono matrix is vertically concatenated to We attempted to directly use the newly formed and ⃗ Y in (5). However, the approach did not produce any substantial improvements from results obtained while using only TEC. We figured that because, generally, the TEC observations are much higher in number than ionosonde observations, and all observation types in ⃗ Y are weighted the same, then the resulting electron densities will generally turn out as an average quantity dominated by TEC estimates. To circumvent the TEC dominance in the final solution, we augmented a data assimilation (DA) procedure to the current tomography algorithm. In addition, through the DA step, we Each datum affects the closest grid points (black asterisks) within a predefined radius of influence r. The magnitude of influence at each grid point is determined by fractional weights w i could implement a positivity constraint that is mandatory in ionospheric imaging but missing in the current tomography algorithm.

Data assimilation (DA) step
A full introduction to the application of DA in ionospheric studies is beyond the scope of this paper, rather for background we recommend Bust et al. (2004), Daley and Barker (2000), and Daley (1991) to reach a more complete understanding of the terminologies used in this section. More specifically, we take a different approach and look at the problem from a probabilistic point of view. The cost function in (4) can be formulated based on a simplistic assumption that ⃗ n follows Gaussian statistics. However, in reality, ionospheric densities are better approximated by a log-normal distribution, whose natural logarithm then assumes Gaussian statistics (Bust and Datta-Barua 2014;Garner et al. 2005). Therefore, it is imperative that we avoid expedient direct application of the Gaussian statistics and rather solve for a more accurate Moreover, assuming a log-normal distribution guarantees that the densities are always positive. The penalty of the modification in (7) is that the relationship between and � ⃗ n is no longer direct but rather nonlinear, and so an additional complex linearization step is required as detailed below.
Consider ⃗ X b to represent background model log densities, to be defined later, and assume both observational and model errors to be uncorrelated. Then, we set the new cost function to be minimized as follows R and B are symmetric and positive definite data and model error covariance matrices, respectively. B acts as a solution stabilizer and also defines how new information should be spread within the specified grid. However, in reality, the true format of B is unknown, and its approximations are computationally difficult. Moreover, even an approximation to a full B generally requires an exhausting process of defining optimal correlation lengths (L). If these L values are poorly overestimated, the algorithm will smooth out spatial variabilities, hence reducing fidelity. In this study, we rather use a diagonal B that is tractable (Ssessanga et al. 2018). That is to say, we annul any interaction or spread of information within the voxel field. This "characteristic" follows from the original tomography (5). To improve accuracy and also reduce the computation cost, solutions to � ⃗ n are only computed for voxels intersected by rays during ray-tracing. Of course, this approach has a downside in that the reconstructed image will always have empty patches in areas where measurements are lacking. Every point in space is not necessarily traversed by the rays that link the GNSS satellites to the receivers.
To minimize the cost function in (8), we take the gradient of J( ⃗ X) and at extreme when where ⃗ X a is the analysis log density at optimal (extreme). Equation (9) still contains a nonlinear term ( ⃗ X a ) which we approximate through the tangent linear model (TLM) and an iterative procedure: where ⃗ X j a and ⃗ X j−1 a are analysis log densities at iterations j and j-1, respectively, and Jacobian Substituting (10) into (9) and rearranging, gives 1 3 93 Page 6 of 15 Equation (11) is the DA step, where we have utilized for B and R being symmetric matrices (Zou et al. 1997, "Appendix" A). The initialization of the iterative procedure in (11) requires the use of background estimates ⃗ X b and the initial best estimate to ⃗ X j−1 a (at j = 1). In most cases, ⃗ X b could have been obtained from a physics-based or empirical model (Thompson et al. 2006;Ssessanga et al. 2019). However, such models are generally empirically biased or with inherent errors that sometimes grow rapidly, thus reducing the fidelity of the final image. To mitigate such errors, we avoid this traditional approach and rather compute both ⃗ X b and ⃗ X j−1 a from outputs of the original tomography. That is to say, under the current tomography (5), eight solutions are computed based on a predefined set of carefully selected eight λ parameters (Saito et al. 2017). We then select solutions (⃗ n 2 , ⃗ n 1 ) corresponding to the two least L-2 norm values. Elements (i) of ⃗ X b and ⃗ X 0 a (j = 1) are then computed as x bi = log(|n 2 i |) and x 0 ai = log (|n 1 i |), respectively. The absolute (| ⋅ |) is taken because log-normal distribution is only defined on the positive axis.
In the DA step (11), it is possible to control the influence of any observation used in the analysis based on the set level of observation uncertainty. That is to say, R assumes the sum of the instrumental and representative error covariances and the values or weightings set in this matrix will determine the influence of a particular data type on the final image. Since ionosonde observations ( ⃗ Y iono ) are considered to be of higher accuracy compared to TEC observations ( ⃗ Y TEC ), the former are weighted 40% (determined from experience) less as shown in (12). The quantities p,l , l TEC , and l iono are the pth column of , index of lth observation in column vector ⃗ Y, number of TEC observations, and ionosonde observations, respectively. R is computed as a diagonal matrix with the assumption that observational errors follow Gaussian statics, are not correlated, unbiased, and are from independent and identically distributed distributions, In the process of determining R, a quality control step is also invoked: A datum y l that satisfies the condition |A l ⃗ n 1 − y l | | > 2 is considered a potential outlier, which is then weighted five times more than the other data. is standard deviation of Δ ⃗ Y = ⃗ n 1 − ⃗ Y (Dee et al. 2011;Ssessanga et al. 2019).
Convergence Eq. (11) is considered to have reached optimality if the mean chi-squared value is less or equal to 0.5. Otherwise, the maximum number of iterations is set to 6. From experience, on average the optimal solution is generally obtained within the first three iterations and the computation time is a fraction of the current latency time of 6 min. Therefore, the added DA complexity does not affect the real-time applications of the algorithm.
Inclusion of IRI model To further improve fidelity, the background model (NeQuick) used in computing weight matrix ⃗ n in (5) was revised to include densities (in the range 100-1500 km) from a continuously upgraded and recommended international reference ionosphere (IRI) model (Bilitza et al. 2011). Moreover, the most recent version used in this study, namely, IRI-2016 (Bilitza et al. 2017), has an improved hmF2 model (Satellite and Digisonde Model of the F2-Layer Height (SDMF2)) developed based on a large amount of radio occultation (RO) and data from digisondes (Shubin et al. 2013;Shubin 2015).

Results and discussion
In this section, we verify and assess the performance of the modified tomography, hereafter referred to as Tomo_DA, in comparison with the original tomography technique hereafter referred to as Tomo_Nq.
In our first analysis, we were interested in how the new modifications would affect the fidelity while considering only GNSS ground observations from 200 receivers, marked as blue asterisks in Fig. 3, that are currently utilized in the real-time tomography as the source of data. This result would act as a benchmark before the inclusion of ionosonde observations and also provide confidence into the linearization and design of the Jacobian in (11). Figure 4 shows an example of spatial analysis densities computed from Tomo_DA (black), alongside those from Tomo_Nq (red). For clarity, the abscissa and ordinate axes have been limited to a magnitude of 15 × 10 11 el/m 3 and an altitude of 800 km, respectively. Epoch 01:00 UT day of year 092, 2016, was selected because negative densities were evident in Tomo_Nq reconstructions. Gaps in some profiles are a result of computing solutions to only those voxels intersected by rays. Eminent is that Tomo_DA has a better resolution of the vertical structure. The profiles are less "noisy" and negative densities are no more. At locations where Tomo_Nq has a good estimation of the profile, Tomo_DA does not significantly change the result. Owing to some unresolved bias and probably the choice of background model ( ⃗ X b ), some profiles still exhibit nonrealistic analysis densities below about 200 km height range. This is well illustrated in the bottom first two subplots of Fig. 4. These nonrealistic densities were mainly found in the E-region, which is difficult to resolve while using absolute STEC: The E-region contributes a very small percentage to ground GNSS-STEC (also see Saito et al. 2017). However, the accuracy of absolute GNSS-STEC, currently ingested, can be obtained to within an accuracy of a few TEC units (1 TEC unit = 10 l6 el/m 3 ), thus, in most cases making a resolve of small densities such as those in E-region uncertain. Consequently, we need additional measurements, here ionosonde data, to specify the electron density vertical structure more uniquely.
In what follows next, we examined the accuracy of the original tomography compared to modified tomography, but this time with including ionosonde data. That is, data from eight ionosonde stations marked as red stars in Fig. 3 were assimilated to outputs from the ground-GNSS-STEC tomography inversions to obtain the final solution. See Table 1 for the ionosonde four-character codes, geographic location and coordinates. The analysis focuses on October and November in years 2015-2016, previously pointed out by Saito et al. (2017; as the period when the tomography algorithm had an under-performance in tracking the observed hmF2 values (refer to Fig. 1). The utilized ionosonde data are accessible at http:// wdc. nict. go. jp/ IONO/ HP2009/ ISDJ/ manual_ txt-E. html and ftp:// ftp. ngdc. noaa. gov/ ionos onde/ data/. The former and latter links provide manually scaled peak parameters and autoscaled full bottomside vertical  profiles, respectively. In the algorithm, we give first priority to manually scaled data due to its higher accuracy. However, the most readily available data are auto scaled. Thus, the results presented here should be analyzed with a consideration of error bounds associated with using ionosonde autoscaled data. Bamford et al. (2008) and Stankov et al. (2012) have reported on these error bounds with a 95% probability; foF2 (− 0.75, + 0.85 MHz), foF1(− 0.25, + 0.35 MH z), foE (− 0.35, + 0.40 MHz), h′F2 (− 68, + 67 km), h′F (− 3 8, + 32 km), and h′E (− 26, + 2 km). Figure 5 shows an example of an extract of vertical profiles from the 3-D reconstructions. Each sub-plot panel represents a specific epoch, and each sub-plot represents vertical profiles at a specific station labeled by the first two station code characters to avoid crowding the picture. Red, blue, green, and yellow represent profiles from original tomography with weight matrix ( ⃗ n) computed using NeQuick model (Tomo_Nq), modified tomography with ⃗ n computed using a combination of IRI-2016, and NeQuick (Tomo_IRI), Tomo_DA and observed ionosonde data (Obz), respectively. Because ground-based ionosondes can only measure electron densities to the NmF2 altitude and the topside profiles are estimated using the technique by Reinisch and Huang (2001), the validation analysis hereafter is limited to the ionosphere bottomside structure.
Noticeable in Tomo_IRI results in that inclusion of IRI in the range 100-1500 km has a profound improvement in tracking observed hmF2, but fails to adequately follow the observed NmF2. On the other hand, Tomo_DA results estimate the NmF2 parameter and the bottomside vertical structure better than both Tomo_Nq and Tomo_IRI. However, during epoch 00:00 UT (Local Time = UT + 9), at Okinawa station (OK), we see an enhanced E-region in the Tomo_DA results that is not reflected in observed data and nonpersistent in the epochs that follow. This could be due to the errors in the linearization procedure. That is, in cases where a trans-ionospheric signal traverses a highly nonlinear region, the assumed first-order approximation TLM model might not fully represent the dynamics, hence the distortion in the final image. Furthermore, under highly nonlinear conditions, the assumption that all inherent errors are Gaussian might not entirely be true, so that the behavior and to what the iterative procedure converges too could be obscure. We should also note that the height of ionosonde observations, specifically hmF2 and hmE, used in assimilation is not a measured quantity. It is Fig. 5 A plot of an extract of vertical profiles from the original tomography (red), modified to include IRI (blue), ionosonde data assimilation (green), and observed ionosonde bottomside data (dashed). To avoid crowding the figure, ionosonde stations are labeled by the first two station code characters. The assimilation step clearly improves the bottomside structure, specifically in tracking the NmF2 parameter both vertically and horizontally rather obtained from inversion schemes, such as POLAN (Titheridge 1985) and NHPC (Huang and Reinish 1996;Reinish et al. 2005) that are subject to systematic errors (Šauli et al. 2007). In addition, currently, the minimum set vertical grid spacing is 20 km, which might be inappropriate to resolve small-scale densities variations like those in the E-region.
Although the focus here is on the bottomside structure, it is interesting to note that for some cases, for example, see profiles in the left-hand middle panel in Fig. 5, the combination of TEC and ionosonde data affects the thickness of the profile and deforms the topside. This is rather an important result indicating that we need to include a thickness constraint in (8) or add another data type, maybe radio occultation observations and in situ satellite densities, that will cover and stabilize the topside.
In Fig. 6, we extract peak parameters NmF2 and hmF2 from Tomo_DA and Tomo_Nq results and compare them to observed values sampled at 1-h intervals and covering months October-November 2015-2016. Left and right panels in each plate represent plots of NmF2 and hmF2, respectively, for a particular station. The significance of assimilating ionosonde is evident at all ionosonde stations. Tomo_DA results are scattered better around the optimal lines (dashed lines), with hmF2 exhibiting the best relative improvement from original tomography estimates. Indeed, Tomo_Nq NmF2 and hmF2 show a positive and negative bias, respectively, but with the latter more pronounced. This result corroborates Saito et al. (2017Saito et al. ( , 2019 findings in which Tomo_Nq was found to overestimate the hmF2 values during the October-November month. As mentioned earlier, this Tomo_Nq under-performance is due to the missing horizontal information vital in reproducing the correct vertical structure, particularly the peak density location. Also, we have already noted from the profiles in Fig. 5 that the choice of using NeQuick as the background in the F-region limits the hmF2 fidelity. Also, as seen in Fig. 6, in some cases the modified technique (Tomo_DA) fails to readjust NmF2 to the correct altitude and rather reverts to the background values (Tomo_Nq). For example, see WK546, YG431, and BP440 tomography hmF2 results in the range 400-600 km. We think that in those particular cases, our choice of the data variance, R, could have been too large when compared to B, such that the analysis densities reverted to background values.
At Guam station, GU513, Tomo_Nq NmF2 distribution has a nonrealistic horizontal stratification. This is probably attributed to Guam's geomagnetic location (5.83° N, 143.51° W) within the equatorial ionization anomaly region. This region is expected to have steep latitudinal densities gradients. However, the currently utilized grid assumes a coarse horizontal spacing (5° × 5°) over this region as shown in Fig. 3, in addition to a pre-setup linear-inverse assumption of uniform electron density within Dashed lines represent the 1:1 reference line. At all stations, original tomography overestimates the hmF2 parameters each voxel. Consequently, the distribution of STEC over this location will not reflect a true representation of the ionospheric state and dynamics. Therefore, we recommend that as the algorithm development expands to include other data types to cover remote areas, the grid should be revised according to prior knowledge of the ionospheric spatial variability. Table 1 lists a summary of the statistics from a comparison of observations to estimates. RMSE_Nq(DA), % impr (= 100 × (RMSE_Nq-RMSE_DA)/RMSE_Nq), and chi_Nq(DA) represent original tomography (data assimilation) root mean square error, percentage improvement and a measure of goodness of fit Chi-squared value, respectively. At some stations like Kokubunji and Wakkanai, NmF2 estimates show an RMSE degradation of 44.22% and 20.47%, respectively. This could be due to RMSE as a measure of accuracy being sensitive to outliers; seen to exist in the Kokubunji and Wakkanai Tomo_DA results in the range 20-40 el/m 3 , Fig. 6. These outliers could be a result of errors in the background model ( ⃗ X b ). Thus, we suggest that in the next algorithm upgrade, the background could be computed from the previous epoch analysis solution so that at least information from earlier times is propagated forward. The choice of the hyperparameter β used in defining R during a particular epoch might also play a major role in filtering out data outliers. Therefore, we suggest a future statistical study to determine a set of β parameters that are suitable for different ionosphere conditions.
Rather than using an outlier sensitive criterion such as RMSE, a better measure of best fit between two "competing" or "nested" models would be to look at changes in chi-squared values (% chi_diff = 100 × (Chi_Nq-Chi_DA)/ Chi_Nq). Here, "nested" refers to Tomo_Nq being a subset of Tomo_DA, which has an extended complexity of assimilation. In Table 1, both under NmF2 and hmF2, the smaller chi_DA values compared to Chi_Nq values suggest that Tomo_DA exhibits a better fit to the observed data. For NmF2, Tomo_DA improvement is most evident over stations Jeju, I-Cheon, Beijing and Gaum that are located in areas not densely covered by the 200 ground GNSS receivers (see Fig. 3). For hmF2, the superiority of using DA to include ionosonde data in the analysis is distinct and well pronounced. At all stations, the relative change in the goodness of fit is above 65% and on average 85%.
Specifically of interest is Wakkanai that has the least improvement in the estimation of NmF2 and hmF2, with the exception of Kokubunji NmF2. This is due to the poor initial guess from tomography at Wakkanai during the DA iterative procedure. That is to say, geometry limitations and minimum boundary conditions limit GNSS observations over Wakkanai, such that the tomography solution is generally acute. Contrary, for Kokubunji, which is quasi-central located within the GNSS receiver network, tomography reconstructed densities have a good accuracy, such that assimilation of ionosonde data does not reflect a large improvement (26.3%) compared to other stations. Nonetheless, on average, there is an improvement at all ionosonde stations that are widely distributed within the specified grid. It would not be premature to conclude that including ionosonde data in the analysis is mandatory to improve the accuracy of the overall 3-D ionosphere picture.

Comparison with Swarm satellite data during storm conditions
To investigate the capability of the modified technique under severe ionospheric conditions, we reconstructed densities during the St. Patrick's day solar storm of March 17-18, 2015 and compared them to observations from Swarm satellites. Previous analyses of ionospheric data covering this storm period showed that the Asian sector recorded nearly a 60% reduction in electron content (Astafyeva et al. 2015;Nava et al. 2016). Thus, these results are part of our future analysis that will explore in detail the stresses of this new technique under an extensive coverage of different ionospheric conditions.
The Swarm constellation consists of three identical satellites, in polar orbits at two different altitudes of ≤ 460 km for Alpha(A) and Charlie(C), and ≤ 530 km for Bravo(B). Each satellite has a Langmuir Probe, which facilitates the measurement of in situ electron densities, considered as the truth in our analysis. The considered data have a 2 Hz resolution and are accessible at ftp:// swarm-diss. eo. esa. int. Since satellites A and C orbit at the same altitude, only data from A and B are presented.
In Fig. 7, in situ time-density profiles (red) are plotted together with densities from our reconstructions (blue) and IRI model (black). Focus is on March 18, 2015, when the negative ionospheric response was most pronounced. The satellite tracks during that period are marked red on the maps in the upper right corner. Because our reconstructions are of a lower time resolution, i.e., 15 min, we assumed the ionosphere to be stationary for a period of 10 min. The vertical dashed lines indicate the start (green) and end (magenta) of grouped traces, in which the time data gaps are less than 1-h. Note that the x-axis is not linear.
The climatological model IRI shows the expected monthly average values of terrestrial plasma densities and frequencies (Bilitza et al. 2017). The reconstructed densities are in good agreement with in situ densities. Also noticeable is that due to differences in orbital altitudes, there is a slight difference in the level of electron density distributions observed by satellites A and B. The previously reported negative ionosphere response of about 60% away from the monthly average values is evident in our reconstructions and in situ observations. On average, the difference between in situ and reconstructed densities is less than 0.2 × 10 12 el/ m 3 . This result offers confidence in the reconstructed 3-D structure, specifically for the topside F-region, which is unattainable when using ground-based instruments such as ionosondes.

Conclusions and future prospects
The near-real-time tomography technique developed to cover the East Asian sector has previously been validated as a reliable source for probing the regional ionosphere 3-D structure. However, that technique still allowed improvements, specifically during the time period of October-November. Thus, this study has discussed the steps undertaken to improve the fidelity of the 3-D structure, especially focusing on estimations of NmF2 and hmF2 peak parameters.
We implemented a data assimilation (DA) procedure that enabled adding ionosonde data into the analysis. Furthermore, we managed to impose the mandatory ionospheric imaging positivity constraint through the DA step by assuming that a log-normal distribution approximates the ionospheric densities better.
Results have shown that the modifications to the tomography technique provide a more reliable 3-D structure, with the bottomside having extended the accuracy due to the inclusion of ionosonde observations. The efficiency and superiority of the modified technique are distinct in the estimation of NmF2 and hmF2, but particularly the latter. New hmF2 estimates fit the observations with improvements above 60% over the original tomography. Overall, results have shown the effectiveness of data assimilation over areas that were scantily covered by the GPS data, but now partially covered by ionosonde observations.

Future prospects
This study has proposed a way to complement the GNSS tomography with another observation type, here ionosonde data. As algorithmic development moves forward and more data sources become abundant, we expect to further improve the fidelity in future research.
In the current implementation, assimilation of ionosonde observations utilizes virtual heights that are subject to systematic errors (Šauli et al. 2007). However, since the current scheme caters for nonlinearity, future modifications should consider a more accurate nonlinear assimilation that involves true heights. In addition, the current technique is still subject to post-analysis. However, there is a possibility of extending the assimilation of ionosonde data to a near-real-time process in the future. . Red marks on maps in the upper right corner indicated the satellite tracks during that period. Trace discontinuities indicate periods when either SWARM or reconstructed densities were absent. The vertical dashed lines show the start (green) and end (magenta) of grouped densities with time data gaps less than 1 h. The x-axis has been readjusted for clarity Since the inclusion of different data types into the algorithm has proven to help in improving the overall accuracy of the 3-D picture, our next step will include an investigation into the addition of radio occultation (RO) data in the analysis. RO data are expected to introduce the missing information that covers the ionosphere in a tangential geometry. RO data will also help cover remote areas, specifically over the low and equatorial latitudes, where the distribution of ionosondes and ground GNSS receivers is significantly lacking.
In the E-region where densities are expected to be quite small, we sometimes observed large enhancements that could be related to the acuteness of the inverse problem. Therefore, there is a possibility that we need to investigate our choice of hyperparameters, linearization, and optimal vertical grid resolution that can resolve small densities in relation to a trade-off with computation cost. Further, we propose that instead of ingesting absolute GNSS-STEC that are determined with large uncertainty, to give preference to using time derivatives of GNSS phase differences that are calculated with high accuracy.
The preliminary results presented in this study are based on a dataset covering only two months in the years 2015 and 2016. Thus, care should be taken when generalizing the result. To make conclusive statistical arguments, we will include in future works a repeat of the same analysis while using a larger database that expands to include more ionospheric dynamics.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.