An Improved Propagation Prediction Method of Low-Frequency Skywave Fusing Fine Channel Parameters

: Low-frequency communication constitutes a vital component of essential communication systems, serving a pivotal role in remote radio communication, navigation, timing, and seismic analysis. To enhance the predictive precision of low-frequency skywave propagation and address the demands of engineering applications, we propose a high-precision prediction method based on the ITU-R P.684 wave-hop theory and real-time environmental parameter forecasts. This method features several distinctive attributes. Firstly, it employs real-time ionospheric prediction data instead of relying on long-term ionospheric model predictions. Secondly, it utilizes a detailed map of land–sea surface electrical characteristics, surpassing the simplistic land–sea dichotomy previously employed. Compared with measured data, the findings demonstrate that we attained a reasonable propagation pattern and achieved high-precision field strength predictions. Comparatively, the improved method exhibits an improvement in the time and spatial domains over the ITU-R P.684 standard. Finally, the improved method balances computational efficiency with enhanced prediction accuracy, supporting the advancement of low-frequency communication system design and performance evaluation.


Introduction
Low-frequency (LF) communication operates within a frequency band defined by the International Telecommunication Union (ITU), spanning from 30 kHz to 300 kHz, with wavelengths ranging from 1000 m to 10,000 m [1].Renowned for its extensive coverage; capability of spanning tens of thousands of kilometers; and remarkable penetration through various media like rocks, soil, and seawater, LF communication serves as a vital tool for command and control communication across vast land areas, including deep bunkers and tunnels.Additionally, LF propagation remains resilient against anthropogenic interference and nuclear detonation, making it important for communications under extreme environmental conditions.To establish and ensure the reliability of LF communication systems, in-depth studies on low-frequency propagation characteristics have always been a theme of the past, present, and future.
In 1925, Appleton proposed the magneto-ionic theory [2], laying the theoretical foundation for the propagation of radio waves in the ionosphere.Subsequently, Booker extended this theory [3], and Bracewell discovered the relationship between ionospheric variations and solar angles [4].The Cavendish Laboratory conducted significant research on LF propagation.Finally, Budden systematically summarized these findings in his book Radio Waves in the Ionosphere, contributing to the continuous improvement and maturation of LF propagation theory [5].Based on the work of these pioneers, Wang combined the latest database to analyze in detail four LF skywave propagation prediction methods, comparing them with the measured field strength, and discussed the factors affecting skywave propagation in 1999 [6].In 2000, Warrington and Jones used the summertime ionospheric model to predict the LF skywave field strength and compared it with the measured field strength to determine the validity of the model [7].In 2002, Liu et al. calculated and compared the geometric parameters of LF skywave propagation based on the distribution models of different atmospheric refractive indices [8].In 2004, Wakai used the wave-hop method to predict the field strength of skywaves, ground waves, and synthetic waves within 60 kHz to 500 kHz within 4000 km and then compared it with the measured field strength to verify the effectiveness of the method [9].In 2006, Wakai used the same method to predict the skywave field strength of 40 kHz and 60 kHz in Japan and compared it with the measured field strength [10].In 2008, Wang qualitatively and quantitatively discussed the seasonal variation in LF skywave field strength based on a huge LF skywave database [11].In 2010, Shigeru and Kurihara proposed the development of a receiver system with higher sensitivity and resolution to reduce interference due to the existence of high-order ionospheric reflection mode in LF propagation and multiple noises in long-distance propagation [12].In 2012, Feng and his team proposed a formula for calculating the reflection coefficient of vertically polarized waves under isotropic conduction conditions based on the oblique reflection propagated by low ionospheric LF propagation, which was well verified [13].In 2017, Pal et al. qualitatively explained LF propagation anomalies by studying the effect of sudden stratospheric warming events in 2009 on LF signals in the Earth-Ionosphere Waveguide (EIWG) [14].In the same year, Zhou et al. combined the Finite-Difference Time-Domain (FDTD) method with the Parabolic Equation (PE) method to analyze the propagation characteristics of LF radio waves at short propagation distances with high propagation angles [15].In 2019, Melchinov predicted relatively accurately the seasonal variation in LF radio waves on the radio propagation path of frozen soil [16].In 2020, based on the wave-hop theory and FDTD method, Zhou et al. studied the time-delay characteristics of a one-hop skywave in an EIWG within 200 km under different ionospheric conditions and different calculation methods [17].In 2021, Gasdia and Marshall developed software (released as the LongwaveMode Propagator.jl) similar to the Long Wave Propagation Code (LWPC) but more stable in mode selection than LWPC [18].In 2022, Mu et al. verified the effects of multiple multi-path delay estimation algorithms in different channel environments and obtained the delay characteristics of LF multi-hop skywaves [19].In the same year, Zhao et al. combined the atmospheric model with the International Reference Ionosphere (IRI) model and then used the improved FDTD method to analyze the characteristics of LF skywave signals [20].In 2023, Gao et al. proposed an LF propagation prediction model that considered the geomagnetic effect based on the wave-hop method, which was used to predict the one-hop skywave field strength under different conditions [21].
To address the demands of engineering applications while balancing computational complexity and efficiency and to refine the prediction methods for 40 kHz LF skywave propagation, this study adopts the ITU-R P.684 method (identified as ITU in the following text) as a foundation.Utilizing the Kriging method, we reconstruct real-time data from four ionospheric parameter observation stations and further enhance the electrical characteristic parameters using a global 16-level classification of surface properties.With these inputs, we predict the field strength of 40 kHz signals across a 2000 km range.Subsequently, we conduct spatiotemporal and comprehensive analyses, comparing different methodologies.

Method of ITU
Skywave propagation refers to the propagation of radio waves back to the ground receiving point following reflection or refraction through the ionosphere.The amplitude and phase of LF skywaves undergo significant variations during day and night and across seasons.Moreover, they are susceptible to irregular and abnormal changes due to the influence of solar and geophysical events.ITU recommends utilizing the wave-hop method for calculating the field strength of LF skywaves.This method directly addresses the field at the receiving point without delving into the intricacies of the propagation process.In the wave-hop method, the total field strength of the received signal is a synthesis of ground wave field strength and skywave field strength.Particularly in short-range propagation, ground wave propagation exhibits distinct advantages.As the propagation distance increases, the skywave gradually assumes a dominant position.In this paper, the 40 kHz field strength is calculated based on the wave-hop method to validate and improve the LF skywave propagation prediction method.Within 2,000 km, there are four possible modes of skywave propagation.Figure 1a,b both depict one-hop skywave propagation.In (a), radio waves emitted by the transmitter are directly received by the receiver following ionospheric reflection.In (b), radio waves emitted by the transmitter initially diffuse along the ground to point A, then incidentally reach the ionosphere tangentially, reflecting at point P and subsequently reaching point B on the ground tangentially.The waves then continue to diffuse along the ground to the receiver.Figure 1c,d illustrate two-hop skywave propagation.In (c), radio waves emitted by the transmitter are received by the receiver after two ionospheric reflections and one ground reflection.In (d), radio waves emitted by the transmitter initially diffuse along the ground to point A, then incidentally reach the ionosphere tangentially.After two ionospheric reflections and one ground reflection, they reach point B on the ground tangentially and proceed to diffuse along the ground to the receiver.The bold parts in Figure 1b,d show the path where the diffraction occurs.
Skywave propagation refers to the propagation of radio waves back to the ground receiving point following reflection or refraction through the ionosphere.The amplitud and phase of LF skywaves undergo significant variations during day and night and acros seasons.Moreover, they are susceptible to irregular and abnormal changes due to the in fluence of solar and geophysical events.ITU recommends utilizing the wave-hop method for calculating the field strength of LF skywaves.This method directly addresses the field at the receiving point without delving into the intricacies of the propagation process.I the wave-hop method, the total field strength of the received signal is a synthesis o ground wave field strength and skywave field strength.Particularly in short-range prop agation, ground wave propagation exhibits distinct advantages.As the propagation dis tance increases, the skywave gradually assumes a dominant position.In this paper, the 4 kHz field strength is calculated based on the wave-hop method to validate and improv the LF skywave propagation prediction method.Within 2,000 km, there are four possibl modes of skywave propagation.Figure 1a,b both depict one-hop skywave propagation In (a), radio waves emitted by the transmitter are directly received by the receiver follow ing ionospheric reflection.In (b), radio waves emitted by the transmitter initially diffus along the ground to point A, then incidentally reach the ionosphere tangentially, reflectin at point P and subsequently reaching point B on the ground tangentially.The waves the continue to diffuse along the ground to the receiver.Figure 1c,d illustrate two-hop sky wave propagation.In (c), radio waves emitted by the transmitter are received by the re ceiver after two ionospheric reflections and one ground reflection.In (d), radio wave emitted by the transmitter initially diffuse along the ground to point A, then incidentall reach the ionosphere tangentially.After two ionospheric reflections and one ground re flection, they reach point B on the ground tangentially and proceed to diffuse along th ground to the receiver.The bold parts in Figure 1b,d show the path where the diffractio occurs.The field strength E at the receiving point is the synthesis of the field generated by the ground diffraction wave and the field generated by the skywave, that is, where E g is the ground wave field strength [22] and E s m is the effective field strength of the m-hop skywave.For a short vertical dipole antenna, the effective field strength of the m-hop skywave is given in the following: where P t is the radiated power of the transmitter (kW) and L is the skywave path length, which is expressed as follows: where d is the great circle distance (km) between the transmitter and the receiver, R e is the effective earth radius (km), and ψ is the elevation angle, which is expressed explicitly as follows: where h r is the ionospheric reflection height, which is expressed as follows: where f is the frequency (kHz), f b is the fundamental frequency (set at 10 kHz), f O E is the critical frequency of the E layer [23], f min (the minimum of f O E at a point), f max (the maximum of f O E at a point), f 0 (f O E when the cosine value of the solar zenith angle χ is 0), h max is the maximum peak height of layer E (set at 100 km), y min is the minimum half-thickness of layer E (set at 10 km), and y max is the maximum half-thickness of layer E (set at 30 km).
In Formula (2), D is the ionospheric focusing factor [24], which is expressed as follows: where R(ρ) is the correction factor, k is the number of wavelengths (=2π/λ), and θ is the angular distance between the launching points and receiving points; when the wave is tangent to the ground, θ is θ c , which can be expressed as follows: (1) where λ is the wavelength (km) and H (1) where R m,n is the ionospheric reflection coefficient of the n-th reflection point of the m-hop skywave.It is a function of f cos(i) and cos(χ) and is calculated via interpolation under given conditions [25].i is the ionospheric incidence angle, which is where F t and F r are the antenna factors of the transmitter and receiver, respectively.These factors are derived based on different conductivities σ and dielectric constant ε and then interpolated according to the elevation angle ψ and the frequency f.R g is the ground reflection coefficient, which is expressed as follows: where E D is the diffracted field strength, E F is the free-space field strength (which can be calculated with the "LFMF-Smooth Earth (v1.1)" software described in ITU-R P.368), and n is the complex dielectric constant, which is expressed as follows: The ionosphere serves as the upper boundary of LF skywave propagation.The precision of predicting its reflection height directly influences the accuracy of LF skywave propagation forecasts.In view of the real-time performance of ionospheric detection data, we used the ionospheric reconstruction method to achieve the real-time prediction of ionospheric parameters so as to improve the accuracy of LF skywave propagation prediction.According to the method of ITU, the ionospheric reflection height is related to f O E. By calculating surrounding detection data, the f O E can be reconstructed, enabling the determination of the ionospheric reflection height for LF propagation reflection.Here, we choose the Kriging method for real-time prediction.
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 data.It uses known sample values and a variance function describing the spatial correlation of the observed sample (calculated by the difference between two pairs of observations at a given distance) to determine unknown values at different spatial locations [26].
When using the Kriging interpolation method to reconstruct, the procedure involves establishing the variance function to obtain spatial autocorrelation of the variables.Subsequently, statistical correlation between observed values is estimated, and the target position is reconstructed.In ionospheric parameter reconstruction, the Kriging interpolation method is specifically applied using "ionospheric distance" to construct the variance function in the Kriging equations.The reconstructed result of the improved Kriging method is defined as the linear weighting of the known ionospheric f O E values, which can be expressed as follows: where is the observed value of f O E at the i-th station, and λ i and φ i are the latitude and longitude of the i-th station, respectively.The weight coefficient is calculated using the Lagrange multiplier algorithm.To ensure the mathematical expectation of the error of the estimated value zero and the sum of squares between the reconstructed value and the observed value minimum, the optimal weight coefficient w i can be obtained from the linear Kriging equation: where µ is the Lagrange factor, r ij is the "ionospheric distance" between station i and j, r j0 is the "ionospheric distance" between station j and the reconstructed location, and the sum of all weighting coefficients w i is 1 to ensure unbiased results.
The described process outlines the reconstruction process of the ionospheric parameter f O E. The core is to use the variance function established through "ionospheric distance" to identify spatial correlation, determine the weight coefficient for each station, and finally obtain the reconstructed value for the target position through the linear combination of the known observed value and the obtained weight coefficients.The variation function of ionospheric distance directly influences the reliability of f O E reconstruction.This paper adopts the modified Euclidean distance with a scale factor under standard geographic coordinates [27], expressed as follows: where s is the scale factor.In summary, the specific implementation process of the f O E reconstruction method is as follows: based on the determined ionospheric distance calculation method, the semivariance function of the ionospheric distance between observation stations is estimated.
Simultaneously, the ionospheric distance vector between the unknown location (λ, φ) and the known observation station (λ 0 , φ 0 ) is calculated.This calculation includes the estimation of the semi-variance function R 0 , namely: Remote Sens. 2024, 16, 2241 7 of 20 where T represents the transpose of the matrix.Formula ( 17) is used to calculate the weight W, which is expressed through vector representation: The f O E for any location is calculated using Formula ( 16).Based on this, the ionospheric reflection height of the corresponding region can be obtained using Formulas ( 6)-( 8).In East Asia, for example, observational data from four stations, Kokubunji, Okinawa, Wakkanai, and Yamagawa (in Figure 2), are used.The variation pattern of ionospheric reflection height h r at 12 Japan Standard Time (JST) and the corresponding f O E of different methods can be compared.As shown in Figure 3, the region colors from dark to light represent h r and f O E from low to high.It can be seen that the distribution of f O E and h r calculated based on the ITU method shows a trend of north-south change (in (a) and (b)).The f O E and h r distribution were calculated using the measured data of the four stations, and the Kriging method shows a trend of change from northeast to southwest (in (c) and (e)).Accordingly, the different distributions of the two show a changing trend from northwest to southeast.At least in East Asia, the distribution presented by the method based on measured data in this paper is closer to the real situation than that offered using the ITU method [28].
[ ] where T represents the transpose of the matrix.Formula ( 17) is used to calcu weight W, which is expressed through vector representation: The fOE for any location is calculated using Formula ( 16).
Based on this, the ionospheric reflection height of the corresponding region obtained using Formulas ( 6)-(8).In East Asia, for example, observational data fr stations, Kokubunji, Okinawa, Wakkanai, and Yamagawa (in Figure 2), are used.T iation pattern of ionospheric reflection height hr at 12 Japan Standard Time (JST) corresponding fOE of different methods can be compared.As shown in Figure 3, th colors from dark to light represent hr and fOE from low to high.It can be seen distribution of fOE and hr calculated based on the ITU method shows a trend o south change (in (a) and (b)).The fOE and hr distribution were calculated using th ured data of the four stations, and the Kriging method shows a trend of chan northeast to southwest (in (c) and (e)).Accordingly, the different distributions of show a changing trend from northwest to southeast.At least in East Asia, the dist presented by the method based on measured data in this paper is closer to the re tion than that offered using the ITU method [28].

Refined Land-Sea Surface Electrical Characteristics
As depicted in Figure 4a,b, the blue squares denote sea areas, while the other colored squares represent land areas.Various colors signify different classifications of electrical characteristics, detailed in Table 1.For example, blue squares indicate a permittivity of 80 F/m and a conductivity of 5 S/m for corresponding sea areas.Figure 4a shows that, according to the ITU method, only blue and dark green areas are distinguished in East Asia due to their simplistic and rather coarse land-sea classification.To enhance prediction accuracy, we adopt a global distribution of 16-level land-sea surface electrical characteristics.This provides a more detailed division of land types, such as dry land and wetland, whose electrical characteristics vary and do not align with the ITU's assumptions.As depicted in Figure 4b, the improved method reveals "Light Blue, Green, Grey, and Purple" areas in East Asia, each representing distinct electrical characteristics.This results in a more detailed classification compared to the ITU method.For example, the electrical characteristics of Japan's main islands undergo changes following the improvement.

Refined Land-Sea Surface Electrical Characteristics
As depicted in Figure 4a,b, the blue squares denote sea areas, while the other colored squares represent land areas.Various colors signify different classifications of electrical characteristics, detailed in Table 1.For example, blue squares indicate a permittivity of 80 F/m and a conductivity of 5 S/m for corresponding sea areas.Figure 4a shows that, according to the ITU method, only blue and dark green areas are distinguished in East Asia due to their simplistic and rather coarse land-sea classification.To enhance prediction accuracy, we adopt a global distribution of 16-level land-sea surface electrical characteristics.This provides a more detailed division of land types, such as dry land and wetland, whose electrical characteristics vary and do not align with the ITU's assumptions.As depicted in Figure 4b, the improved method reveals "Light Blue, Green, Grey, and Purple" areas in East Asia, each representing distinct electrical characteristics.This results in a more detailed classification compared to the ITU method.For example, the electrical characteristics of Japan's main islands undergo changes following the improvement.

Data Acquisition
In this research, we acquired the 40 kHz measurement data published by the ITU and referred to the radio wave transmitter specifications in Table 2.The transmitter position was set at (37.4°N, 140.8°E), and field measurements were conducted in November 2013 using a ship equipped with a receiver.The receiver fitted with a short vertical dipole antenna performed LF propagation data acquisition every hour.Finally, a total of 78 point-

Data Acquisition
In this research, we acquired the 40 kHz measurement data published by the ITU and referred to the radio wave transmitter specifications in Table 2.The transmitter position was set at (37.4 • N, 140.8 • E), and field measurements were conducted in November 2013 using a ship equipped with a receiver.The receiver fitted with a short vertical dipole antenna performed LF propagation data acquisition every hour.Finally, a total of 78 pointto-point link measurements from 55 receiver position points within a range of 2000 km were obtained.

Analysis of Improved Method
According to the improved method proposed in this paper, the variation in propagation parameters is compared before and after applying the improved prediction method for the mentioned stations and their surrounding areas.
In Figure 5, the first line shows the variation in ψ, the second line shows the variation in i, the third line shows the variation in D, the fourth line shows the variation in F r , and the fifth line shows the variation in R. In each line, the first column shows the parameter variation of the ITU method, the second column shows the parameter variation of the improved method (identified as the PRO method), and the third column shows the difference between ITU and PRO.The region's color from light to dark shows the parameter values from high to low.
As observed in Figure 5, the introduction of the fine-channel parameter correction method exhibits noticeable differences compared to the traditional ITU method: (1) For most links in the eastern part of the region, the elevation angle calculated via the PRO method is increased.(2) For most links in the eastern part of the region, the ionospheric incidence angle calculated via the PRO method is reduced.(3) The ionospheric focusing factor and antenna factor are different between the two methods.(4) The ionospheric reflection coefficient slightly decreases around the transmitting point after correction.
Whether the prediction accuracy of the field strength brought by the PRO method is improved will be analyzed in the following chapters.In order to facilitate the analysis, Absolute Error (AE) and Root-Mean-Square Error (RMSE) are selected as the evaluation criteria.The specific calculation methods are as follows: where E OBS is the measured field strength (identified as OBS) and n is the total number of samples.
PRO method is increased.(2) For most links in the eastern part of the region, the ionospheric incidence angle calculated via the PRO method is reduced.(3) The ionospheric focusing factor and antenna factor are different between the two methods.(4) The ionospheric reflection coefficient slightly decreases around the transmitting point after correction.
Whether the prediction accuracy of the field strength brought by the PRO method is improved will be analyzed in the following chapters.In order to facilitate the analysis, Absolute Error (AE) and Root-Mean-Square Error (RMSE) are selected as the evaluation criteria.The specific calculation methods are as follows: where EOBS is the measured field strength (identified as OBS) and n is the total number of samples.

Analysis of Time Domain Characteristics
As shown in Figure 6, we gathered the 40 kHz skywave propagation field strength received by the receiver at location R as a fixed reception point.The data were collected on a stationary ship between 5 and 6 November 2013 and included changes over a 24 h period.
the PRO method.(g) D of the ITU method.(h) D of the PRO method.(i) The difference in D b the ITU method and the PRO method.(j) Fr of the ITU method.(k) Fr of the PRO method difference in Fr between the ITU method and the PRO method.(m) R of the ITU method.the PRO method.(o) The difference in R between the ITU method and the PRO method.

Analysis of Time Domain Characteristics
As shown in Figure 6, we gathered the 40 kHz skywave propagation field st received by the receiver at location R as a fixed reception point.The data were co on a stationary ship between 5 and 6 November 2013 and included changes over period.Figure 7 compares the variation in the predicted field strength and the measure strength within 24 h between the ITU and PRO methods.As shown in the figure, t dicted field strength exhibits smooth changes at night, while significant fluctuation during the day.At the same time, we can see that the time domain characteristics night before and after improvement demonstrate good consistency.This is beca value of fOE during the nighttime is minimal, making it difficult for observation s to detect fOE and leading to a loss of nighttime data.This loss presents significan lenges for predictions based on measured data.To address this issue, the paper the steady-state value of fOE during the nighttime, which closely aligns with t method.This method allows for the prediction of a steady-state value for field st However, because the measured field strength is only transient and easy to fluctu fluctuation of the predicted field strength at night is obviously weaker than that measured field strength.In contrast, the predicted field strength fluctuates signif during the day before and after the improvement.It can be seen that the predicte strength fluctuates greatly between 5:00-8:00 and 15:00-19:00.This is because fOE w duce large fluctuations in these periods (near sunrise and sunset), which will have impact on the ionospheric reflection height.As a result, the predicted field streng tuates more significantly, which is well reflected in both the ITU and PRO metho pecially between 7:00 and 8:00, the predicted field strength difference based method in this paper reaches 3.55 dB.However, fOE will basically remain in a stead without major fluctuations during 9:00-14:00, and the corresponding iono Figure 7 compares the variation in the predicted field strength and the measured field strength within 24 h between the ITU and PRO methods.As shown in the figure, the predicted field strength exhibits smooth changes at night, while significant fluctuations occur during the day.At the same time, we can see that the time domain characteristics of the night before and after improvement demonstrate good consistency.This is because the value of f O E during the nighttime is minimal, making it difficult for observation stations to detect f O E and leading to a loss of nighttime data.This loss presents significant challenges for predictions based on measured data.To address this issue, the paper adopts the steady-state value of f O E during the nighttime, which closely aligns with the ITU method.This method allows for the prediction of a steady-state value for field strength.However, because the measured field strength is only transient and easy to fluctuate, the fluctuation of the predicted field strength at night is obviously weaker than that of the measured field strength.In contrast, the predicted field strength fluctuates significantly during the day before and after the improvement.It can be seen that the predicted field strength fluctuates greatly between 5:00-8:00 and 15:00-19:00.This is because f O E will produce large fluctuations in these periods (near sunrise and sunset), which will have a great impact on the ionospheric reflection height.As a result, the predicted field strength fluctuates more significantly, which is well reflected in both the ITU and PRO methods.Especially between 7:00 and 8:00, the predicted field strength difference based on the method in this paper reaches 3.55 dB.However, f O E will basically remain in a steady state without major fluctuations during 9:00-14:00, and the corresponding ionospheric reflection height will not change significantly.Therefore, it is reasonable for the predicted field strength to remain at a steady-state value.reflection height will not change significantly.Therefore, it is reasonable for the predicted field strength to remain at a steady-state value.Figure 8 is a radar map of the variation of the error over a 24 h period, where th height of the green column is the magnitude of the error at the corresponding time.It can be seen that the column height during the day is generally reduced after improvement However, due to the consistency of the time domain characteristics at night, the column height at the corresponding time remains the same.Considering the difference between the daytime and nighttime models, we present the RMSE of the measured field strength and the predicted field strength during the daytime and night based on the two method in Table 3.It is evident that compared with the ITU method, the PRO method enhance prediction accuracy during the daytime.
Through the analysis, the RMSE between the measured field strength and the pre dicted field strength of the PRO method within 24 h is about 1.78 dB (where the RMSE between the measured field strength and ITU is about 1.91 dB), about 6.81% higher than the ITU method.Figure 8 is a radar map of the variation of the error over a 24 h period, where the height of the green column is the magnitude of the error at the corresponding time.It can be seen that the column height during the day is generally reduced after improvement.However, due to the consistency of the time domain characteristics at night, the column height at the corresponding time remains the same.Considering the difference between the daytime and nighttime models, we present the RMSE of the measured field strength and the predicted field strength during the daytime and night based on the two methods in Table 3.It is evident that compared with the ITU method, the PRO method enhances prediction accuracy during the daytime.reflection height will not change significantly.Therefore, it is reasonable for the predicted field strength to remain at a steady-state value.Figure 8 is a radar map of the variation of the error over a 24 h period, where the height of the green column is the magnitude of the error at the corresponding time.It can be seen that the column height during the day is generally reduced after improvement.However, due to the consistency of the time domain characteristics at night, the column height at the corresponding time remains the same.Considering the difference between the daytime and nighttime models, we present the RMSE of the measured field strength and the predicted field strength during the daytime and night based on the two methods in Table 3.It is evident that compared with the ITU method, the PRO method enhances prediction accuracy during the daytime.
Through the analysis, the RMSE between the measured field strength and the predicted field strength of the PRO method within 24 h is about 1.78 dB (where the RMSE between the measured field strength and ITU is about 1.91 dB), about 6.81% higher than the ITU method.Through the analysis, the RMSE between the measured field strength and the predicted field strength of the PRO method within 24 h is about 1.78 dB (where the RMSE between the measured field strength and ITU is about 1.91 dB), about 6.81% higher than the ITU method.

Analysis of Spatial Domain Characteristics
As shown in Figure 9, the receiver, following the survey ship along the Great Circle course from 8 to 10 November 2013, was used to receive 40 kHz skywaves from the JJY 40 kHz transmitter.That is, the position of the receiver changes with time.Different receiver position points are marked in the figure, which can also be seen as different receivers in different locations.

Analysis of Spatial Domain Characteristics
As shown in Figure 9, the receiver, following the survey ship along the Great Circle course from 8 to 10 November 2013, was used to receive 40 kHz skywaves from the JJY 40 kHz transmitter.That is, the position of the receiver changes with time.Different receiver position points are marked in the figure, which can also be seen as different receivers in different locations.Figure 10 shows the variation in the predicted field strength with distance d between the two methods from JST = 0 to JST = 23 [29].The horizontal axis represents the distance and the vertical axis represents the magnitude of the field strength.The spatial domain characteristics indicate that the predicted field strength exhibits significant ups and downs when the propagation distance is relatively close.The predicted field strength gradually demonstrates a gentle fluctuation trend as the distance increases.During the day, the pre dicted field strength shows a slight decline near 500 km, while during the night, a sub stantial decline occurs near 700 km.Through comparison, it is evident that the spatia domain characteristics before and after improvement remain consistent during the night while significant changes are observed during the day.
Figure 11 is a comparison of errors, where the column height is the magnitude o errors, and the horizontal axis is the corresponding receiving point (from near to far).I can be seen that the height of the column at medium and short distances is higher than that at long distances.According to statistics, 45 receiving points have an error difference of less than 1 dB, with 21 receiving points having the same error before and after improve ment, particularly at night.The error of the 13th receiving point (about 700 km) is the largest, about 20 dB.
It is apparent that in the medium and short distance, where the receiving point is located, significant fluctuations and large errors are observed, possibly due to interference Figure 10 shows the variation in the predicted field strength with distance d between the two methods from JST = 0 to JST = 23 [29].The horizontal axis represents the distance, and the vertical axis represents the magnitude of the field strength.The spatial domain characteristics indicate that the predicted field strength exhibits significant ups and downs when the propagation distance is relatively close.The predicted field strength gradually demonstrates a gentle fluctuation trend as the distance increases.During the day, the predicted field strength shows a slight decline near 500 km, while during the night, a substantial decline occurs near 700 km.Through comparison, it is evident that the spatial domain characteristics before and after improvement remain consistent during the night, while significant changes are observed during the day.
Figure 11 is a comparison of errors, where the column height is the magnitude of errors, and the horizontal axis is the corresponding receiving point (from near to far).It can be seen that the height of the column at medium and short distances is higher than that at long distances.According to statistics, 45 receiving points have an error difference of less than 1 dB, with 21 receiving points having the same error before and after improvement, particularly at night.The error of the 13th receiving point (about 700 km) is the largest, about 20 dB.two methods.It can be seen that the error is larger at near and medium distances.Ho ever, at medium distances, the performance of the PRO method is significantly better th that of the ITU method.At far distances, the RMSE for both methods is greatly reduc and shows consistency, but the performance of the PRO method needs to be improv compared with the ITU method.Through the analysis, the accuracy of the PRO method improved by about 3.69% compared to the ITU method.ITU method.Red solid line: predicted field strength of the PRO method.Horizontal axis: distance (km).Vertical axis: field strength (dB (μV/m)).

Comprehensive Analysis
Figure 12 shows the overall comparison between the measured field strength and the predicted field strength based on the two methods and AEs between them.Through analysis, the RMSE between the predicted field strength based on the PRO method and the measured field strength of 78 improved point-to-point links is about 5.94 dB (where the RMSE between the two based on the ITU method is 6.17 dB), about 3.73% higher than that of the ITU method.It is apparent that in the medium and short distance, where the receiving point is located, significant fluctuations and large errors are observed, possibly due to interference caused by the strong ground wave in this distance range.However, as the distance increases, the ground wave is significantly weakened, leading to a reduction in field strength fluctuations and smaller errors.Table 4 shows the RMSE between the measured field strength and the predicted field strength at near, medium, and far distances for the two methods.It can be seen that the error is larger at near and medium distances.However, at medium distances, the performance of the PRO method is significantly better than that of the ITU method.At far distances, the RMSE for both methods is greatly reduced and shows consistency, but the performance of the PRO method needs to be improved compared with the ITU method.Through the analysis, the accuracy of the PRO method is improved by about 3.69% compared to the ITU method.

Comprehensive Analysis
Figure 12 shows the overall comparison between the measured field strength and the predicted field strength based on the two methods and AEs between them.Through analysis, the RMSE between the predicted field strength based on the PRO method and the measured field strength of 78 improved point-to-point links is about 5.94 dB (where the RMSE between the two based on the ITU method is 6.17 dB), about 3.73% higher than that of the ITU method.
As shown in Figure 13, we conducted statistics on the overall error between the two methods (for ease of analysis, we removed the absolute value of AE, the measured value minus the predicted value).It can be seen that the error distribution roughly conforms to a normal distribution.Before the improvement, the measured field strength of 48 points is less than the predicted field strength, and the measured field strength of 30 points is greater than the predicted field strength.There are 19 points with an error ranging from −1 to 1 dB.After the improvement, the measured field strength of 46 points is less than the predicted field strength, and the measured field strength of 32 points is greater than the predicted field strength.There are also 19 points with an error of −1 to 1 dB, and 13 points with an error of 1 to 3 dB.
Figure 12 shows the overall comparison between the measured field strength and the predicted field strength based on the two methods and AEs between them.Through analysis, the RMSE between the predicted field strength based on the PRO method and the measured field strength of 78 improved point-to-point links is about 5.94 dB (where the RMSE between the two based on the ITU method is 6.17 dB), about 3.73% higher than that of the ITU method.As shown in Figure 13, we conducted statistics on the overall error between the two methods (for ease of analysis, we removed the absolute value of AE, the measured value minus the predicted value).It can be seen that the error distribution roughly conforms to a normal distribution.Before the improvement, the measured field strength of 48 points is less than the predicted field strength, and the measured field strength of 30 points is greater than the predicted field strength.There are 19 points with an error ranging from −1 to 1 dB.After the improvement, the measured field strength of 46 points is less than the predicted field strength, and the measured field strength of 32 points is greater than the predicted field strength.There are also 19 points with an error of −1 to 1 dB, and 13 points with an error of 1 to 3 dB.Section 3.2 illustrates that variation in ionospheric reflection height impact relevant parameters, ultimately influencing field strength.Considering the method outlined in Section 2, each parameter is interrelated, so improvements in one parameter alone cannot be isolated when assessing field strength changes.Notably, all parameters are directly or indirectly related to ψ, meaning even minor adjustments to ψ can cascade changes across parameters and subsequently impact field strength.Therefore, the alterations resulting from "Real-time reconstruction of ionospheric parameters" primarily stem from ψ.A statistical analysis reveals that the modification in ψ ultimately enhances field strength prediction accuracy by approximately 3.55% in this study.
By analyzing the time and spatial domain characteristics, Figure 14a,b compare the traditional and improved methods based on ITU.Using the transmitter as the center, we forecasted the field strength at JST = 12 in November 2013 over a range of 2000 km, including the Great Circle track.The change in the regional color from dark to light is the magnitude of field strength from low to high.It can be observed that the field strength fluctuation roughly follows the distribution form of concentric rings.When the propagation distance is relatively short, the regional color changes frequently, and the field strength fluctuates significantly.With the increase in distance, the regional color changes less, and the field strength fluctuates less, aligning with the spatial domain characteristics.Upon comparison, as opposed to the traditional method, the predicted field strength obtained with the PRO method exhibits more fluctuations when the propagation distance is Section 3.2 illustrates that variation in ionospheric reflection height impact relevant parameters, ultimately influencing field strength.Considering the method outlined in Section 2, each parameter is interrelated, so improvements in one parameter alone cannot be isolated when assessing field strength changes.Notably, all parameters are directly or indirectly related to ψ, meaning even minor adjustments to ψ can cascade changes across parameters and subsequently impact field strength.Therefore, the alterations resulting from "Real-time reconstruction of ionospheric parameters" primarily stem from ψ.A statistical analysis reveals that the modification in ψ ultimately enhances field strength prediction accuracy by approximately 3.55% in this study.
By analyzing the time and spatial domain characteristics, Figure 14a,b compare the traditional and improved methods based on ITU.Using the transmitter as the center, we forecasted the field strength at JST = 12 in November 2013 over a range of 2000 km, including the Great Circle track.The change in the regional color from dark to light is the magnitude of field strength from low to high.It can be observed that the field strength fluctuation roughly follows the distribution form of concentric rings.When the propagation distance is relatively short, the regional color changes frequently, and the field strength fluctuates significantly.With the increase in distance, the regional color changes less, and the field strength fluctuates less, aligning with the spatial domain characteristics.Upon comparison, as opposed to the traditional method, the predicted field strength obtained with the PRO method exhibits more fluctuations when the propagation distance is relatively short, and the changes induced by the PRO method gradually decrease with the distance increase.Figure 14c shows the AE between (a) and (b).It can be seen that when closer to the transmitter, there are noticeable differences in field strength.However, as the distance from the transmitter increases, the changes in field strength become less significant.

Discussion
In this section, we will follow the method and process outlined, replacing the Kriging method with two other commonly used interpolation methods (Nearest and Linear) to seek model optimization.Tables 5 and 6 show the improvement effects of these three

Discussion
In this section, we will follow the method and process outlined, replacing the Kriging method with two other commonly used interpolation methods (Nearest and Linear) to seek model optimization.Tables 5 and 6 show the improvement effects of these three interpolation methods compared to the ITU method in the time domain, spatial domain, and overall.While the Nearest interpolation method has shown remarkable improvement in the time domain, it has hardly shown any enhancement in the spatial domain; instead, it has increased computational complexity.Linear interpolation, although exhibiting some improvement in the time domain, performs even worse in the spatial domain compared to the ITU method, and it also increases computational complexity.In contrast, the Kriging method utilized in this paper integrates improvements in both time and spatial domains, meeting the requirements of engineering applications and achieving multidimensional enhancement in LF propagation prediction.In conclusion, the Kriging method chosen in this paper demonstrates significant advantages.
However, the PRO method still has several limitations: (1) As mentioned earlier, f O E is difficult to observe at night.Even during the daytime, the observation equipment may be susceptible to various factors such as anthropogenic interference, leading to some reading errors.Therefore, data missingness and misalignment can greatly impact predictions based on the measured data.(2) The effectiveness of the Kriging method is easily influenced by the distances between ionospheric observation stations and between these stations and the target location.

Conclusions
Based on the wave-hop theory as outlined in the ITU method, this study incorporated real-time ionospheric forecast data in place of the long-term predictions typically used by the ionospheric model.Additionally, a detailed land-sea surface electrical characteristics map was utilized instead of a simple land-sea differentiation.These refinements improved the channel parameters, leading to greater accuracy in predicting LF skywave propagation field strength.By comparing and analyzing measurement data from 55 stations within 2000 km, we established a reasonable skywave propagation pattern and achieved accuracy improvements of 6.81% in the time domain and 3.69% in the spatial domain.Overall, the method showed an approximate improvement of 3.73% based on comprehensive analysis.Combined with Section 4, although the PRO method still has certain limitations, it achieves multi-dimensional accuracy improvements while maintaining computational efficiency, making it highly significant for engineering applications.Future research endeavors may focus on further optimizing prediction results by selecting more suitable observation stations, increasing their number, and refining ground electrical characteristics.Moreover, deeper exploration and analysis of long-distance multi-hop skywave propagation can be pursued by expanding the sample size.

20 Figure 3 .
Figure 3.Comparison of variation pattern of ionospheric reflection height hr and fOE when JST = 12 in East Asia.(a) fOE of the ITU method.(b) fOE of the improved method.(c) The difference in fOE between the ITU method and the improved method.(d) hr of the ITU method.(e) hr of the improved method.(f) The difference in hr between the ITU method and the improved method.

Figure 3 .
Figure 3.Comparison of variation pattern of ionospheric reflection height h r and f O E when JST = 12 in East Asia.(a) f O E of the ITU method.(b) f O E of the improved method.(c) The difference in f O E between the ITU method and the improved method.(d) h r of the ITU method.(e) h r of the improved method.(f) The difference in h r between the ITU method and the improved method.Remote Sens. 2024, 16, x FOR PEER REVIEW 9 of 20

Figure 4 .
Figure 4. Changes in land-sea surface electrical characteristics in East Asia.(a) The land-sea distinction in the ITU method.(b) Refined land electrical characteristics.

Figure 4 .
Figure 4. Changes in land-sea surface electrical characteristics in East Asia.(a) The land-sea distinction in the ITU method.(b) Refined land electrical characteristics.

Figure 5 .
Figure 5.Comparison of the variation in relevant parameters in East Asia.(a) ψ of the ITU method.(b) ψ of the PRO method.(c) The difference in ψ between the ITU method and the PRO method.(d) i of the ITU method.(e) i of the PRO method.(f) The difference in i between the ITU method and the PRO method.(g) D of the ITU method.(h) D of the PRO method.(i) The difference in D between the ITU method and the PRO method.(j) F r of the ITU method.(k) F r of the PRO method.(l) The difference in F r between the ITU method and the PRO method.(m) R of the ITU method.(n) R of the PRO method.(o) The difference in R between the ITU method and the PRO method.

Figure 6 .
Figure 6.Location of the transmitter and fixed receiver.Green •: transmitter location.Blue ▼ receiver location.

Figure 6 .
Figure 6.Location of the transmitter and fixed receiver.Green •: transmitter location.Blue ▼: fixed receiver location.

Figure 7 .
Figure 7.The variation pattern of field strength with time.Blue •: measured field strength at th corresponding moment.Red square line: predicted field strength of the ITU method.Green triangl line: predicted field strength of the PRO method.

Figure 7 .
Figure 7.The variation pattern of field strength with time.Blue •: measured field strength at the corresponding moment.Red square line: predicted field strength of the ITU method.Green triangle line: predicted field strength of the PRO method.

Figure 7 .
Figure 7.The variation pattern of field strength with time.Blue •: measured field strength at the corresponding moment.Red square line: predicted field strength of the ITU method.Green triangle line: predicted field strength of the PRO method.

Figure 9 .
Figure 9. Location of the transmitter and mobile receivers.Green •: transmitter location.Blue ▼ mobile receiver locations.

Figure 9 .
Figure 9. Location of the transmitter and mobile receivers.Green •: transmitter location.Blue ▼: mobile receiver locations.

Figure 10 .
Figure 10.The variation pattern of field strength with distance at different hours.Blue •: measur field strength at the corresponding distance.Blue double-dash line: predicted field strength of t

Figure 10 .
Figure 10.The variation pattern of field strength with distance at different hours.Blue •: measured field strength at the corresponding distance.Blue double-dash line: predicted field strength of the ITU method.Red solid line: predicted field strength of the PRO method.Horizontal axis: distance (km).Vertical axis: field strength (dB (µV/m)).

Figure 13 .
Figure 13.Overall error statistics of the two methods.(a) ITU method.(b) PRO method.Blue column: histogram of frequency distribution.Red solid line: normal distribution curve.

Figure 13 .
Figure 13.Overall error statistics of the two methods.(a) ITU method.(b) PRO method.Blue column: histogram of frequency distribution.Red solid line: normal distribution curve.
Remote Sens. 2024, 16, x FOR PEER REVIEW 17 of 20 distance from the transmitter increases, the changes in field strength become less significant.

Table 1 .
The electric characteristic parameters of the land-sea surface correspond to different color squares.

Table 1 .
The electric characteristic parameters of the land-sea surface correspond to different color squares.

Table 2 .
Specifications for 40 kHz radio wave transmitters.

Table 3 .
RMSE during day and night based on the two methods.

Table 4 .
RMSE at different propagation distance stages based on the two methods.

Table 4 .
RMSE at different propagation distance stages based on the two methods.

Table 5 .
Time improvement effects of three interpolation methods compared to the ITU method.

Table 6 .
Spatial improvement effects of three interpolation methods compared to the ITU method.