A Reconstruction Method for Ionospheric f o F 2 Spatial Mapping over Australia

: To improve the accuracy of predicting the ionospheric critical frequency of the F2 layer ( f o F 2 ), a reconstruction method for the spatial map of the ionospheric f o F 2 based on modiﬁed geomagnetic dip coordinates is proposed. Based on the strong correlation between the ionospheric f o F 2 and geomagnetic coordinates, the variation function of ionospheric distance is built. In the end, the spatial map of the ionospheric f o F 2 is predicted by solving the Kriging equation. The results show that the regional characteristics of the ionospheric f o F 2 analyzed by the proposed method are consistent with the observations. Compared with the reconstructed value of f o F 2 using traditional geographic coordinates, the root-mean-square error (RMSE) in high solar activity years decreased by 0.43 MHz, and the relative RMSE decreased by 5.48%; The RMSE decreased by 0.35 MHz during low solar activity which is 5.99% lower to relative RMSE. The research results provide support for high-frequency communication frequency selection.


Introduction
The ionosphere is the ionized portion of the upper atmosphere (the thermosphere, part of the mesosphere, and the exosphere). It extends from about 60 to beyond 1000 km and completely encircles the Earth. The primary plasma source for the ionosphere is the photoionization of neutral molecules via solar EUV and soft X-ray radiation, although other production processes may dominate in certain regions. The ions produced undergo chemical reactions with the neutrals, recombine with the electrons, diffuse to higher or lower altitudes, or are transported via neutral wind effects [1]. According to the vertical distribution characteristics of electron density, the ionosphere can be divided into D, E, and F layers, of which the F layer is usually divided into the F1 layer and the F2 layer.
The ionosphere is a time-varying dispersion channel, which varies with seasons, day and night, and regions and is affected by space environments such as solar activity and geomagnetic fields. High-frequency communications are affected by the ionosphere [2,3]. The ionospheric plasma structure on a wide range of space and time scales can have a destructive impact on space application systems (such as navigation and communication systems operated by Global Navigation Satellite Systems (GNSS)) [4], and its spatiotemporal dynamic changes lead to fluctuations in communication channels, which directly affect the reliability and stability of ionospheric communication, this is also the reason why the ionosphere and its communication research has become a research hotspot in the past hundred years [5]. Existing research also shows that the trend scenario of the ionosphere is now more complete than before, but some gaps and differences still need to be further studied [6]. The critical frequency of the F2 layer (f o F 2 ), which represents the minimum frequency of electromagnetic waves that can penetrate the ionosphere vertically, is one of the critical parameters of the ionosphere and can be used to understand ionospheric dynamics and structure [7], as well as the influence of ionospheric weather on communication and navigation [8]. f o F 2 has an irregular and non-uniform spatial distribution, and the unknown data in adjacent regions can be interpolated or extrapolated from the limited known ionospheric observations. This regional ionospheric reconstruction problem is critical in regional ionospheric reporting and prediction [9].
Many organizations and scholars have been continuously studying ionospheric regional reconstruction technology for a long time. At present, the ionospheric spatial features description methods based on different empirical or mathematical methods include inverse distance weight (IDW) interpolation [10], spline interpolation [11], spherical harmonic function [12,13], neural network (NN) [14], deep learning (DL) [15], and Kriging method [16]. Among them, the IDW method lacks the geophysical mechanism and is far from the observation station, so the accuracy is low. The spherical Lagrange function is only orthogonal on the whole sphere and cannot effectively represent the signal in the restricted region on the sphere [17]. Green function is the most robust, efficient, and accurate method [18]. The surface spline interpolation based on the green function has good approximation performance. Due to its high accuracy, simplicity, flexibility, and other advantages, it has become the mainstream method and a potent tool for multivariate approximation [19].
To achieve high-accuracy regional ionospheric foF2 reconstruction, a new f o F 2 reconstruction method for the Australian region is proposed by introducing geomagnetic coordinates and improving the Kriging method. In the proposed method, the variation function of ionospheric distance is determined based on geomagnetic coordinates, and the estimated value is obtained by solving the improved Kriging equation. The proposed method is validated and analyzed by using the data collected from 7 ionospheric stations in Oceania.

Reconstruction Method
The Kriging method is a widely used interpolation algorithm in ionospheric parameter reconstruction, which can provide the best linear unbiased estimation and accurately describe the spatial structure of ionospheric parameters.
When using the Kriging interpolation method for reconstruction, the variation function is first determined to obtain the spatial autocorrelation of variables, the statistical correlation between observations is estimated, and then the target location [20]. The reconstructed result of the Kriging model based on MGD coordinates is defined as the linear weighting of the value of the known ionospheric sampling observations (such as f o F 2 ), and the formula is as follows: where,f o F 2 is the f o F 2 reconstruction value of any location, N is the maximum number of observation stations, w n is the weight coefficients of the nth observation station, f o F 2 (λ n , ϕ n ) is the observation value of the nth station, λ n and ϕ n are the latitude and longitude of the nth observation station respectively. The Lagrange multiplier algorithm is used to calculate the weight coefficients. To make the mathematical expectation of the error of the estimated value of f o F 2 (λ, ϕ) zero and to minimize the sum of squares between the reconstructed value and the observed value, the optimal weight coefficients w n can be solved from the linear Kriging equation [21]: where, N is the maximum number of observation stations used for spatial interpolation, µ is the Lagrange factor, r ij is the ionospheric distance between stations i and j, r i0 is the ionospheric distance between station i and the reconstructed station, and the sum of all weight coefficients is 1 to ensure the unbiased reconstruction results: The above process is the reconstruction of the ionospheric parameter. Its core is to use the variation function established by ionospheric distance to find the spatial correlation, determine the weight coefficients of each ionospheric observation station, and finally solve the reconstructed value of the target location by a linear combination of the known observed value and the calculated weight coefficients [22].
In the current common reconstruction methods, the ionospheric distance is mostly calculated using the modified Euclidean distance with scale factor (SF) in the standard geographic coordinate system, which is the ionospheric distance determined by Stanislawska when he introduced the Kriging interpolation method into the ionospheric field [23]. Considering that the ionosphere is jointly controlled by the direction of the Earth's rotation axis and the geomagnetic field structure, geomagnetic parameters play an important role in the main variations of the ionosphere, especially in the F2 layer, and the variation trend of the F2 layer parameters is strongly correlated with geomagnetic coordinates [24]. Therefore, the model no longer adopts simple traditional geographical coordinates (TGG), but geomagnetic coordinates can be adopted. Modified geomagnetic dip coordinates (MGD) have been proven to be very effective in modeling ionospheric parameters, especially suitable for the ionospheric F2 layer plasma characteristics [25][26][27]. Therefore, the ionospheric distance is determined based on the MGD as follows: where, (ϑ i , φ i ) and (ϑ j , φ j ) correspond to the geomagnetic longitude and latitude of the ith position (λ i , ϕ i ) and the jth position (λ j , ϕ j ) respectively; s is the scale factor, reflecting the different variations of ionospheric parameters in the specific region at the magnetic longitude and latitude. In summary, the specific implementation process of the above reconstruction method is as follows: Based on the determined ionospheric distance calculation method, the semivariation function of the ionospheric distance between observation stations is estimated: Atmosphere 2023, 14, 1399 4 of 12 At the same time, the ionospheric distance vector between the unknown location (λ i , ϕ i ) and the known observation station (λ 0 , ϕ 0 ) can be calculated, including the semivariation function (R 0 ): R 0 = r 10 r 20 · · · r N0 T (6) Formulas (2) and (3) are used to calculate the weight coefficients W, which is expressed by vector representation: Finally, the ionospheric parameters at any position are calculated using Equation (1).

Method Validation
The f o F 2 used in this paper can be obtained from the ftp://ftp-out.sws.bom.gov.au/ wdc/iondata/medians (accessed on 1 May 2023), which contains ionospheric monthly medians obtained from manually scaled hourly ionospheric data. The seven stations selected in Australia are Brisbane, Canberra, Learmonth, Perth, Darwin, Townsville, and Hobart. Figure 1 shows the geographical and geomagnetic location distributions of the stations. Table 1 shows the station's corresponding geographical and geomagnetic coordinates, where the expressions of geomagnetic latitude and modified geomagnetic dip longitude are shown in Equation (8) and Equation (9), respectively: (a) geomagnetic latitude [28]: where λ m is geomagnetic latitude, λ and ϕ are the geographic latitude and longitude, respectively, of the apex, and λ 0 and ϕ 0 are the geographic latitude and longitude, respectively, of the north pole of the earth-centered dipole.
(b) Modified Geomagnetic Dip [27]: where I is the magnetic dip angle, defined by the downward and horizontal components of the field, and λ is geographical latitude. (c) Modified Geomagnetic Dip Longitude [27]: where µ is the modified geomagnetic dip angle, λ m is the geomagnetic latitude, ϕ is the geographic longitude of the apex, and ϕ 0 is the geographic longitude of the north pole of the earth-centered dipole. It can be seen from the Figure 1 as follows: (1) The spatial distribution of the MGD coordinates is significantly different from that of the TGG coordinates; (2) The geomagnetic longitude marked by the blue line has a specific deviation from the geographical longitude line, and the greater the deviation is toward the east; (3) The modified geomagnetic dip latitude marked by the red line significantly differs from the geographical latitude, especially in the low and high-latitude regions.  To estimate the influence of different coordinates and different distance SFs on the reconstruction results, root-mean-square error (RMSE) and relative root-mean-square error (RRMSE) are selected as the evaluation criteria of the proposed method, and the formula is as follows: (a) RMSE: where, f' is the reconstructed value of the ionospheric foF2; f is the observed value of the ionospheric foF2, and H is the statistical number.
For the stations shown in Figure 1, based on the MGD and the observed foF2 values of the seven ionospheric stations from 2013 to 2017, different distance SFs s were analyzed, and the reconstructed RMSEs were shown in Figure 2a below. It can be seen that the RMSE  To estimate the influence of different coordinates and different distance SFs on the reconstruction results, root-mean-square error (RMSE) and relative root-mean-square error (RRMSE) are selected as the evaluation criteria of the proposed method, and the formula is as follows: (a) RMSE:  Figure 2a below. It can be seen that the RMSE decreases gradually and then increases slightly as the SF s increases. When the s is 9.5, the RMSE is the smallest, so the optimal SF is determined. In this case, the corresponding ionospheric distance is shown in Figure 2b, which fundamentally differs from the traditional distance.
Atmosphere 2023, 14, x FOR PEER REVIEW 6 of decreases gradually and then increases slightly as the SF s increases. When the s is 9.5, th RMSE is the smallest, so the optimal SF is determined. In this case, the correspondin ionospheric distance is shown in Figure 2b, which fundamentally differs from the trad tional distance. According to the optimal SF determined above, the foF2 in the Oceania region is r constructed. The typical reconstruction diagram is shown in the following Figure 3.  According to the optimal SF determined above, the f o F 2 in the Oceania region is reconstructed. The typical reconstruction diagram is shown in the following Figure 3.
Atmosphere 2023, 14, x FOR PEER REVIEW 6 of 12 decreases gradually and then increases slightly as the SF s increases. When the s is 9.5, the RMSE is the smallest, so the optimal SF is determined. In this case, the corresponding ionospheric distance is shown in Figure 2b, which fundamentally differs from the traditional distance. According to the optimal SF determined above, the foF2 in the Oceania region is reconstructed. The typical reconstruction diagram is shown in the following Figure 3.  As can be seen from Figure 3, the maximum value appears near the equator, and the f o F 2 gradually decreases with the increase of latitude. With the change of time, the maximum value region has a drift property near the equator, and the drift property is the same in different years. Generally, the f o F 2 of noon is higher than the value of midnight, and the f o F 2 of high solar activity is higher than that of low solar activity. These characteristics are consistent with the trend of ionospheric characteristics in the southern hemisphere, proving this method's effectiveness and applicability.
Further, according to the observation data of the seven stations, the cross-validation method was adopted to reconstruct the 24-h f o F 2 of the seven stations one by one and analyze the error. The following Figure 4 shows the comparison between the observed values and the reconstructed values at noon in low solar activity years (a~g) and high solar activity years (h~n).
Atmosphere 2023, 14, x FOR PEER REVIEW 7 of 12 As can be seen from Figure 3, the maximum value appears near the equator, and the foF2 gradually decreases with the increase of latitude. With the change of time, the maximum value region has a drift property near the equator, and the drift property is the same in different years. Generally, the foF2 of noon is higher than the value of midnight, and the foF2 of high solar activity is higher than that of low solar activity. These characteristics are consistent with the trend of ionospheric characteristics in the southern hemisphere, proving this method's effectiveness and applicability.
Further, according to the observation data of the seven stations, the cross-validation method was adopted to reconstruct the 24-h foF2 of the seven stations one by one and analyze the error. The following Figure 4 shows the comparison between the observed values and the reconstructed values at noon in low solar activity years (a~g) and high solar activity years (h~n).  In Figure 4, the left column is the comparison of the high solar activity year, and the right column is the comparison of the low solar activity year. It can be seen from Figure 4 that the reconstruction accuracy of Canberra, Darwin, and Learmonth stations in low solar activity years is the highest. The observed results are consistent with the reconstruction results. In contrast, the reconstruction accuracy of other stations at night will show slight deviations, and the reconstruction results in other periods are consistent with the observed values. In high solar activity years, the reconstruction results of Brisbane, Canberra, Hobart, and Townsville stations have high accuracy, while the other stations have certain deviations around noon, and the reconstruction results of other times are relatively accurate. Combining the above two groups of curves, the period with high reconstruction accuracy accounted for 81.9%, and the period with slight deviation accounted for 18.1%.
To further illustrate the advantages of using the Kriging method based on TGG coordinates (Kriging-TGG) for ionospheric f o F 2 modeling, we compare the results with those based on MGD coordinates (Kriging-MGD). Figure 5a shows the ionospheric distance calculated using TGG coordinates. Figure 5b shows the absolute value of the difference between the ionospheric distance calculated using MGD coordinates and the ionospheric distance calculated using TGG coordinates, which shows a particular gap in the ionospheric distance calculated under the two different coordinates.
Atmosphere 2023, 14, x FOR PEER REVIEW 8 of 12 In Figure 4, the left column is the comparison of the high solar activity year, and the right column is the comparison of the low solar activity year. It can be seen from Figure 4 that the reconstruction accuracy of Canberra, Darwin, and Learmonth stations in low solar activity years is the highest. The observed results are consistent with the reconstruction results. In contrast, the reconstruction accuracy of other stations at night will show slight deviations, and the reconstruction results in other periods are consistent with the observed values. In high solar activity years, the reconstruction results of Brisbane, Canberra, Hobart, and Townsville stations have high accuracy, while the other stations have certain deviations around noon, and the reconstruction results of other times are relatively accurate. Combining the above two groups of curves, the period with high reconstruction accuracy accounted for 81.9%, and the period with slight deviation accounted for 18.1%.
To further illustrate the advantages of using the Kriging method based on TGG coordinates (Kriging-TGG) for ionospheric foF2 modeling, we compare the results with those based on MGD coordinates (Kriging-MGD). Figure 5a shows the ionospheric distance calculated using TGG coordinates. Figure 5b shows the absolute value of the difference between the ionospheric distance calculated using MGD coordinates and the ionospheric distance calculated using TGG coordinates, which shows a particular gap in the ionospheric distance calculated under the two different coordinates. The absolute value of the difference between the ionospheric distance was calculated using MGD coordinates, and the ionospheric distance was calculated using TGG coordinates.
To measure the performance of the two Kriging methods, more specifically, the RMSE and RRMSE of the two methods are calculated respectively according to Equations (11) and (12). Figures 6 and 7 show the RMSE in summer when solar activity is high and winter when solar activity is low, respectively. Among them, the hollow circle is the RMSE obtained by the Kriging-TGG, and the solid circle is the RMSE obtained by the Kriging-MGD. The size or color of the circle can represent the size of the RMSE. For any station, the RMSE obtained using the Kriging-MGD is smaller than that obtained using the Kriging-TGG. The absolute value of the difference between the ionospheric distance was calculated using MGD coordinates, and the ionospheric distance was calculated using TGG coordinates.
To measure the performance of the two Kriging methods, more specifically, the RMSE and RRMSE of the two methods are calculated respectively according to Equations (11) and (12). Figures 6 and 7 show the RMSE in summer when solar activity is high and winter when solar activity is low, respectively. Among them, the hollow circle is the RMSE obtained by the Kriging-TGG, and the solid circle is the RMSE obtained by the Kriging-MGD. The size or color of the circle can represent the size of the RMSE. For any station, the RMSE obtained using the Kriging-MGD is smaller than that obtained using the Kriging-TGG. Figure 8 shows the difference in RMSE statistically obtained using the two Kriging models, with the horizontal axis representing the size of RMSE and the vertical axis representing the station name. The length of the green column represents the difference in RMSE when the solar activity is high in summer, which shows that the difference in RMSE between the two Kriging models is 0.43 MHz at the maximum and 0.16 MHz at the minimum. The length of the red column represents the difference in RMSE in winter when the solar activity is low, which shows that the difference in RMSE between the two models at Darwin station is the largest, 0.71 MHz. The difference in RMSE between the two models at Townsville station is the smallest, 0.03 MHz. Using the Kriging-TGG, the average RMSE of the seven stations is 0.81 MHz when the solar activity is high in summer and 0.90 MHz when the solar activity is low in winter. Using the Kriging-MGD, the average RMSE of the seven stations is 0.50 MHz when the solar activity is high in summer, and the average RMSE in winter when solar activity is low is 0.55 MHz, which decreases by 0.31 MHz and 0.71 MHz, respectively, compared with the Kriging-TGG coordinates.   Figure 8 shows the difference in RMSE statistically obtained using the two Kriging models, with the horizontal axis representing the size of RMSE and the vertical axis representing the station name. The length of the green column represents the difference in RMSE when the solar activity is high in summer, which shows that the difference in RMSE between the two Kriging models is 0.43 MHz at the maximum and 0.16 MHz at the minimum. The length of the red column represents the difference in RMSE in winter when the solar activity is low, which shows that the difference in RMSE between the two models at Darwin station is the largest, 0.71 MHz. The difference in RMSE between the two models at Townsville station is the smallest, 0.03 MHz. Using the Kriging-TGG, the average RMSE of the seven stations is 0.81 MHz when the solar activity is high in summer and 0.90 MHz when the solar activity is low in winter. Using the Kriging-MGD, the average RMSE of the seven stations is 0.50 MHz when the solar activity is high in summer, and the average RMSE in winter when solar activity is low is 0.55 MHz, which decreases by 0.31 MHz and 0.71 MHz, respectively, compared with the Kriging-TGG coordinates.   Figure 8 shows the difference in RMSE statistically obtained using the two Kriging models, with the horizontal axis representing the size of RMSE and the vertical axis representing the station name. The length of the green column represents the difference in RMSE when the solar activity is high in summer, which shows that the difference in RMSE between the two Kriging models is 0.43 MHz at the maximum and 0.16 MHz at the minimum. The length of the red column represents the difference in RMSE in winter when the solar activity is low, which shows that the difference in RMSE between the two models at Darwin station is the largest, 0.71 MHz. The difference in RMSE between the two models at Townsville station is the smallest, 0.03 MHz. Using the Kriging-TGG, the average RMSE of the seven stations is 0.81 MHz when the solar activity is high in summer and 0.90 MHz when the solar activity is low in winter. Using the Kriging-MGD, the average RMSE of the seven stations is 0.50 MHz when the solar activity is high in summer, and the average RMSE in winter when solar activity is low is 0.55 MHz, which decreases by 0.31 MHz and 0.71 MHz, respectively, compared with the Kriging-TGG coordinates.  Figure 9 shows the RRMSE obtained statistically using two Kriging models. The cylinder's length represents the RRMSE's size, and the blue circle corresponds to the right vertical axis, representing the difference between the RRMSE obtained statistically by the Kriging-TGG and the RRMSE obtained statistically by the Kriging-MGD. Figure 9a shows the statistical diagram of RRMSE in high summer solar activity years, which shows that for any station, the RRMSE obtained by the Kriging-MGD is smaller than that obtained by the Kriging-TGG. The difference in RRMSE of the Hobart station is the smallest, while that of the Brisbane station is the largest. Figure 9b shows the statistical diagram of RRMSE in winter low solar activity years. Except for the Townsville station, the RRMSE obtained by the Kriging-MGD is smaller than that obtained by the Kriging-TGG. The RRMSE obtained by the Kriging-MGD is 1.33% larger than that obtained by the Kriging-TGG. Atmosphere 2023, 14, x FOR PEER REVIEW 10 of 12  Figure 9 shows the RRMSE obtained statistically using two Kriging models. The cylinder's length represents the RRMSE's size, and the blue circle corresponds to the right vertical axis, representing the difference between the RRMSE obtained statistically by the Kriging-TGG and the RRMSE obtained statistically by the Kriging-MGD. Figure 9a shows the statistical diagram of RRMSE in high summer solar activity years, which shows that for any station, the RRMSE obtained by the Kriging-MGD is smaller than that obtained by the Kriging-TGG. The difference in RRMSE of the Hobart station is the smallest, while that of the Brisbane station is the largest. Figure 9b shows the statistical diagram of RRMSE in winter low solar activity years. Except for the Townsville station, the RRMSE obtained by the Kriging-MGD is smaller than that obtained by the Kriging-TGG. The RRMSE obtained by the Kriging-MGD is 1.33% larger than that obtained by the Kriging-TGG. Generally, the average RRMSE of the seven stations is 19.98% in summer when the solar activity is high, 18.81% in winter when the solar activity is low, and 14.50% in summer when the Kriging-MGD is applied. The average RRMSE in winter when solar activity   Figure 9 shows the RRMSE obtained statistically using two Kriging models. The cylinder's length represents the RRMSE's size, and the blue circle corresponds to the right vertical axis, representing the difference between the RRMSE obtained statistically by the Kriging-TGG and the RRMSE obtained statistically by the Kriging-MGD. Figure 9a shows the statistical diagram of RRMSE in high summer solar activity years, which shows that for any station, the RRMSE obtained by the Kriging-MGD is smaller than that obtained by the Kriging-TGG. The difference in RRMSE of the Hobart station is the smallest, while that of the Brisbane station is the largest. Figure 9b shows the statistical diagram of RRMSE in winter low solar activity years. Except for the Townsville station, the RRMSE obtained by the Kriging-MGD is smaller than that obtained by the Kriging-TGG. The RRMSE obtained by the Kriging-MGD is 1.33% larger than that obtained by the Kriging-TGG. Generally, the average RRMSE of the seven stations is 19.98% in summer when the solar activity is high, 18.81% in winter when the solar activity is low, and 14.50% in summer when the Kriging-MGD is applied. The average RRMSE in winter when solar activity Generally, the average RRMSE of the seven stations is 19.98% in summer when the solar activity is high, 18.81% in winter when the solar activity is low, and 14.50% in summer when the Kriging-MGD is applied. The average RRMSE in winter when solar activity is low is 12.82%; compared with the Kriging-TGG, RRMSE decreased by 5.48% and 5.99%, respectively.
According to the above analysis, the Kriging-MGD has advantages over the Kriging-TGG. In other words, using geomagnetic coordinates instead of geographic coordinates to calculate ionospheric distance can realize ionospheric f o F 2 reconstruction with higher accuracy.

Conclusions
In this paper, by introducing modified geomagnetic dip coordinates, we proposed a high-accuracy f o F 2 regional reconstruction model over the Australian region. The regional characteristics are consistent with the objective characteristics of the ionosphere, and the results can reflect the physical characteristics of ionospheric f o F 2 at different locations. In addition, we compared the Kriging-MGD with the Kriging-TGG. The experimental results show that the mean RMSE and RRMSE of Kriging-MGD are 0.50 MHz and 14.50%, respectively, in summer when solar activity is high. The RMSE and RRMSE of Kriging-TGG are 0.81 MHz and 19.98% respectively. When solar activity is low in winter, the average RMSE and RRMSE of Kriging-MGD and Kriging-TGG are 0.55 MHz and 12.82%, respectively, and 0.90 MHz and 18.81%, respectively. The results show that compared with Kriging-TGG, Kriging-MGD has more minor errors and higher accuracy. This technology can be extended to the global region and provide essential support for optimizing longdistance communication systems' frequency.