Enhancement of MF R-Mode ranging accuracy by exploiting measurement-based error mitigation techniques

Nowadays, the maritime sector strongly relies on global navigation satellite systems (GNSS) for the provision of position, navigation, and timing (PNT) information. The standard functionality of several ships’ bridge instruments depends on such information, which becomes critical to have a safe and reliable navigation. Nevertheless, the possibility of GNSS outages combined with unintentional and intentional interference to GNSS signals, which have been increasing over the last years, can severely threaten the nominal activity of the crew onboard of vessels. For these reasons, alternative position, navigation, and timing (APNT) systems become fundamental in order to provide operational continuity. R(anging)-Mode is a terrestrial alternative system to GNSS for the maritime domain. In one of its possible implementation, it exploits synchronized medium frequency (MF) signals transmitted by maritime radio beacons. In the MF band, the radio waves are attenuated and distorted by the change of electrical properties of the ground along the propagation path. Moreover, these signals are affected by terrain elevation changes and large metallic infrastructures which induce an overall loss in the accuracy of the system. In this paper, we propose a mitigation technique based on the estimation of these additional error sources by using in-field range measurements. The theoretical details are explained in depth and the technique is validated with real data gathered on a measurement campaign. A clear and significant performance improvement is achieved on the range accuracy when the proposed approach is applied.

vacuum propagation, based on the knowledge of the atmospheric parameters and the ground conductivity maps.
Most often, only rough ground conductivity maps exist (ITU 2015). When using these for the generation of AGDF maps, this helps to reduce the ranging error in an R-Mode receiver but significant errors remain. Additionally, a delay may also be introduced by the variation in the terrain elevation along the path, resulting in a longer travelled distance. Lastly, large metallic infrastructures, such as bridges or power line plants, can produce a further signal distortion which is complex to be modelled (U.S. Coast Guard 1992).
In order to compensate for the aforementioned additional biases and modelling mismatch, one solution is to exploit direct measurements in the field. This second approach enhances the prediction quality, as it will be shown. In this paper, we show the impact of unmodelled delay sources and propose a novel methodology to generate a correction function dependent on the azimuth which improves the ranging accuracy. Additionally, a simplified technique with a reduced computational complexity is proposed along with the comparison with the main approach.
The paper is divided as follows: Section 2 briefly describes the theory on the AGDF. In Section 3, the new methodology is illustrated in detail, with the description of the main assumptions and validity of the technique. Here, the simplified approach is also described. In Section 4, we apply the proposed approaches on real data from a measurement campaign, presenting the achieved improvement with respect to the application of AGDF correction alone. Section 5 contains a relevant discussion about the operational implications related to the research discoveries. Finally, Section 6 concludes the paper.

Atmospheric and ground wave delay factor
At medium frequencies, ground waves propagate along the earth's surface following its curvature. The attenuation of the wave is a complex valued factor that causes a damping of the waves field strength and a phase delay. It depends on the ground conductivity and permittivity of the surface medium (Wait 1998). The effect can be described by an attenuation function that accounts for ground wave-related effects with respect to free-space propagation. For short distances below 20 km, the extended flat earth solution of the wave equation is used to express the field while for larger distances, the residue series solution is evaluated. By using the ground conductivity, permittivity and distance travelled as input values, the argument of the attenuation function yields the ground wave propagation delay of the wave with respect to vacuum-free space propagation. To correct phase estimates which are used to determine pseudoranges in R-Mode using the vacuum speed of light, the propagation delay can simply be subtracted from the measurement. To obtain the phase delay of a non-homogeneous, segmented propagation path with sections of different electrical properties, the propagation delay of each segment can be concatenated by using the Millington/Pressey (Pressey et al. 1953) method. The correction term derived with the aforementioned method serves as the so-called AGDF for R-Mode.
The composition of a propagation path has to be determined to obtain the AGDF. Since the ITU World Atlas of Ground Conductivities (ITU 2015) provides a database of electrical properties for a large portion of the earth, it can be used to determine the composition by extracting the ground conductivity and distance from segments of equal conductivity along the propagation path from the transmitter to the receiver. For each segment, the attenuation function is calculated using the well-known method introduced by Rotheram (1981aRotheram ( , 1981b and Wait (1998) before concatenating the functions of all segments using the Millington/Presseys method. The solution provided by Rotheram's approach accounts not only for the finite conductivity of a curved earth surface but also for refraction in the exponential atmosphere.
In comparison to Loran-C (U.S. Coast Guard 1992), a terrestrial low frequency navigation system which today is only available in a few parts of the world, the predicted AGDF represents the sum of primary, secondary, and additional secondary factors. In the case of the AGDF, the effects that are modelled with the primary and secondary factors are not computed separately or approximated through polynomials. The separation of these different contributions is not useful since the effects heavily depend on the frequency of the electromagnetic wave. For MF R-Mode, each transmitter uses a separate frequency band, which would require a dedicated polynomial function for the primary and secondary factor contribution. Furthermore, the salinity of the Baltic Sea varies over space and time, which makes a constant factor accounting for pure sea water propagation less useful. Figure 1 provides an example of the AGDF calculated with the ITU-R Ground Conductivity maps for a transmitter located next to Groß Mohrdorf in the north of Germany. The introduced phase delay compared to propagation in vacuum is given here in radians. By observing the gradient in the color map, it is clear that the right side area of the transmitter is affected by a different delay compared to the left side. In particular, a higher delay is expected on the right part due to the presence of additional land (the island of Rugen), which causes a reduction in the propagation speed with respect to the sea.
The accuracy of the correction obtained with the described approach depends strongly on the quality of the ground conductivity maps utilized. If the conductivity value in the maps closely represent reality, the estimated signal delay will be accurate. Unfortunately, the ground conductivity data is quite old and, in general, with low spatial resolution. The conductivity and permittivity of the ground depend of the soil texture, volumetric water content, temperature, and for sea water on the relative salinity. These factors can vary significantly in some areas, which makes the static ground conductivity map less accurate and not parameterizable.
Additionally, this approach does not include a model to compensate for the delay introduced by the variation in the terrain elevation. Also, distortion generated by metallic infrastructures, which may be encountered by the signal along the propagation path, are not incorporated. Therefore, to overcome these limitations and increase the ranging accuracy, the approach described in the next section, based on direct-range measurements, has been established and evaluated.

Generating a direct measurement-based correction function
In this section, we describe how to exploit R-Mode pseudorange measurements in combination with GNSS-based ranges to produce a correction function which increases the R-Mode ranging accuracy. MF R-Mode pseudorange measurements are computed by using phase observation derived from transmitted pilot sinusoidal signals. The concept is similar to the GNSS carrier phase observation and the interested reader can refer to Grundhöfer et al. (2021) for a detailed explanation of the MF R-Mode ranging and positioning principles. The phase variation between transmitter and receiver is equivalent to the travelling time of the signal and therefore the pseudorange ρ R in meters can be expressed as follows: with c the propagation speed in vacuum and t rx and t tx are reception and transmission time respectively. To be precise, phase observation are ambiguous due to the cyclic nature of the sinusoidal signal and a single cycle has an equivalent distance of approximately 1 km. Therefore, an important task of the receiver is the resolution of such an ambiguity. At the time of writing, the DLR software receiver solves the ambiguity by performing an initial calibration process which uses GNSS data, as explained in Grundhöfer et al. (2021). As discussed previously, different impairments affect the signal along its path; hence, the pseudorange model can be expressed as follows with d the geodesic distance between receiver and transmitter, modelled with the Vincenty's formula (Vincenty 1975), b CK the combined receiver and transmitter clock offset, b AGD the delay induced by the atmosphere and ground wave propagation as explained in Section 2. The terrain elevation error is represented by b TE whereas the unmodelled distortions, due to large ferromagnetic obstacles, are described by b O . Last but not least, n is a noise term which depends on the receiver itself and the received signal-to-noise ratio (SNR). With the help of a GNSS receiver, which provides accurate information about the position of the vessel, the transmitter-receiver distance can also be easily calculated and we define the GNSS range as follows where n represents the noise in the range which depends on the working mode of the GNSS receiver and its internal algorithms, going from a few meters for single point positioning (SPP) to a few centimeters for real-time kinematic (RTK) or precise point positioning (PPP). We define the difference between R-Mode and GNSS range as and by substituting (2) and (3) in (4), we obtain with n = n − n . It is clear that ρ represents a noisy version of the delay budget. By using a receiver which is synchronized with the transmitters, the clock bias can be neglected; hence, b CK ≈ 0. Additionally, even in a non-synchronized scenario, this could also be estimated and removed by using transmitter clock offset data correction, which will be broadcasted by the service provider and from the positioning, velocity, and timing (PVT) algorithm of the R-Mode receiver.
We can now apply a filter to remove the noise term n . Different filtering strategies can be applied but for this paper a simple average window forward-backward filtering has been considered. The filtered delay budget, ρ f , can then be represented as follows: as a clean version of the delay components. In principle, if the ground conductivity maps used for the prediction coincide with the real ground conductivity, the bias b AGD can be removed completely. Nevertheless, due to the mismatch between real and assumed values, a residual bias remains, which we define as follows: where b AGDF represents the predicted correction obtained as explained in Section 2. We then define the overall residual error, which we call model error (ME), since it accounts for model mismatch and unknown or unmodelled source of delay (b AGD r , b TE , and b O ), as follows: Finally, a function f ME (α) based on the azimuth angle α (angle formed by the vector connecting the transmitter and the vessel and the transmitter vector pointing to the north) is generated. Such a function can be easily obtained by using interpolation or fitting techniques. In this work, a cubic 1D interpolator has been used.
To clarify the principle of this methodology, we use Fig. 2 to provide a simple explanation to the reader. Suppose the transmitter is located on an island depicted by the gray circle area, the main assumption is that the majority of the delay introduced in the signal is due to the land component of the path. Indeed, on this section, we have higher variability of the ground conductivity, terrain elevation, and potential additional source of distortion. The assumption is that a measurement has been taken on the point P 1 (yellow cross), which has an arbitrary azimuth α. By applying the described process from (4) to (8), one can obtain the ME for that azimuth. Imagine now performing a second measurement on P 2 (red cross) which is located at a longer distance from the transmitter but with the same azimuth angle α of P 1 . Assuming that the ground conductivity maps are accurate on the sea path section, the variation due to the land component should in principle be the same as for P 1 . Therefore, the ME computed for P 1 is also applicable to P 2 .
This description should clarify the overall procedure to compute and apply the proposed correction scheme. With this approach, we try to capture high-order biases introduced in the signal mainly by the land path section between the transmitter and the coast. Therefore, its main limitation is the assumption that no additional land path is encountered between the coast and the vessel. Nevertheless, even in this case, the correction still holds, but it cannot include compensation for the additional ground area which must be then compensated separately. The main advantage of this approach is that, in principle, by only surveying the entire azimuth on the sea once we can apply the correction over a large and extensive maritime area. This clearly reduces not only the effort of planning and conducting measurements on a large area, but also drastically decreases the cost of the survey campaigns.
In this paper, we include a second simplified approach which follows the general idea of generating the ME correction function. In the simplified approach, we assume the R-Mode user has insufficient knowledge of the ground parameters and is not able to create the AGDF correction function for the overland propagation path. In such a case, the map for the delay function is generated based on the assumption of a sea only propagation path from transmitter to the user with constant electrical parameters of the sea water around the transmitter which should be measured in the vicinity of the transmitter. The ground wave propagation delay across sea water is described by a polynomial function, inspired by the secondary factor approach used in Loran (U.S. Coast Guard 1992). The electrical characteristics of sea water are a function of frequency, temperature, and salinity. A detailed explanation is given in ITU (2021), which includes formulas for the computation of conductivity and permittivity based on such parameters. The attenuation function for a given set of frequencies and ground conductivity values is calculated using a modified version of GRWAVE (ITU 2022). Afterwards, a fitting function is derived. The parameters of the polynomial function f Sea of the form are estimated through least squares curve fitting, where x is the distance in kilometers. The resulting polynomial is used to determine the phase delay across sea water as a function of distance from the transmitter. As an example, Table 1 depicts the fitted values of the polynomial parameters for the two continuous wave signal frequencies of Groß Mohrdorf for a salinity of 12.5 PSU and a temperature of 5 • C.
We define this correction term as sea model SM, given that it only models the bias introduced by the speed over sea water. By following the procedure to generate the model error function (from Eqs. 7 to 8), a new correction function, indicated as ME SM , can be created. This function will not only contains the unmodelled error terms but also the error component introduced by the change of ground and water properties along the azimuth.
Similar performance are expected to be observed for these two techniques, at least in areas where the electrical properties of the sea are constant. Moreover, the simplified approach can be particularly useful in areas which present strong deviations between the true ground parameters and the one assumed in the ITU maps.

Results and discussion
In this section, we present and discuss the results obtained, in post-processing, from the application of the suggested technique on real data. The measurement campaign was conducted in February 2021, with the Deneb ship, provided by the German Federal Maritime and Hydrographic Agency (BSH). The ship was equipped with the medium frequency R-Mode receiver designed by the DLR and capable of computing ranges and position of the user (Grundhöfer et al. 2021).The measurement setup was composed by such a receiver connected to a GPS-stabilized rubidium clock, The values are calculated for both the transmitted sinusoidal tones based on a temperature of 5 • C and a salinity of the water of 12.5 PSU used as a reliable and accurate time source. This last component is not needed in the receiver implementation but it is of particular importance for research purposes. Due to the usage of this time source, the unknown time parameter is not anymore relevant assuming an ideal synchronization between the receiver and the transmitter, which is typically equipped with a similar high accuracy clock. For this reason, the term b CK is neglected in the following analysis. Additionally, a GNSS receiver with RTK service was used to record a reference trajectory which we used for assessing the performance of the R-Mode receiver and to perform the initial calibration process, as explained before. The map presented in Fig. 3 contains relevant information about the location of the test area and the transmitter. For this paper, the Groß Mohrdorf transmitter, visible in the map as yellow triangle, has been considered for the analysis. The blue line, indicated as forward path, describes the ship's movement from the shore to the sea on 20th February which was used for generating the correction function whereas the red line, indicated as backward path, represents the trajectory of the ship from the sea to the shore on the same day and was used for validation. The green line is the trajectory of the vessel for the 24th February which was also entirely used as validation datatset.
As explained in Section 3, the first step to start building the correction function is to compute the error between the R-Mode ranges and GNSS-based ranges. For the forward path, this is represented in Fig. 4 as blue line and indicated as ρ. We can see that the initial bias is 0 due to the calibration process. It tends to grow, in the absolute value, presenting several variations over time. We remind that during Fig. 3 Map showing the test scenario. The yellow triangle represents the Groß Mohrdorf transmitter location, the blue line is the forward path (shore to sea) whereas the red line shows the backward path (sea to shore) for the day 20th February. The green line represents the path of the vessel on 24th February  Fig. 4 Comparison between true range error ρ (blue line), b AGDF (red dotted), and b SM (orange dashed) prediction this time, the vessel was moving from the shore to the sea and we assume that the measured bias is due to the change in signal propagation speed caused by the landsea path variation. On top of that, the effect of noise is also clearly visible. On the same figure, the predicted bias (b AGDF ), based on the AGDF maps, is also plotted as red dotted line. It is easy to see that the predicted bias follows the general trend of the true one. Therefore, the AGDF maps are suitable for correcting large errors. Nevertheless, it appears clear that not all the variations are accurately modelled due to the mismatch between the maps and the measurement, plus the possible additional error sources, as discussed in Section 2. Thus, the difference between the two curves represents the component of the bias which is not compensated and that we want to model with the new approach for future usage. In the same figure, the correction term obtained based on SM is visible as orange dashed line. As expected, such a correction does not provide accurate results. It can be noted that the curve is almost linear and has a small decreasing trend which depends on the fact that the receiver-transmitter distance was decreasing over that time.
By applying a 30 -s average window filter in forward-backward mode, the noise is cancelled. Afterwards, by using Eq. 8, the difference between the predicted bias, based on the SM and AGDF, and the new filtered measured error is used to obtain the correction terms b ME SM and b ME which are represented in Fig. 5 as a function of the azimuth angle. The angle is defined in the range (−180, 180] northwards. Finally, an interpolation process is used to build the final functions f ME (α) and f ME SM (α) which can be then used to predict the residual error terms. Several interpolation algorithms might be considered but, in this work, a cubic 1D interpolator was used. We can now test the two functions on the backward path to validate the proposed methodology. Figure 6 includes five curves representing the R-Mode range error without any correction in blue, the predicted error based on the AGDF maps as red dotted line, the predicted error based on AGDF maps and the estimated azimuthal correction f ME (α) depicted as black dash-dot line, the predicted error based on SM in orange dashed, and finally the predicted error based on SM and the estimated azimuthal correction f ME SM (α) in green. The range error and the prediction curves start from zero because they were all initially calibrated. In contrast to Fig. 4, the range error increases over time due to the fact the ship was moving in the opposite direction, as expected. As for the forward path, despite the AGDF-based prediction follows the increasing trend of the range error reducing the bias of the measurement, the faster variation are not properly modelled, resulting in an overall residual error. Also, the predicted bias based on SM does not improve the accuracy of the range, as it is almost constant.
It can be observed that with the use of the ME functions, the prediction quality substantially increases. The fluctuations are followed accurately, especially between 16:00 and 16:30 where there is very good match between the blue, black, and green curves. A small increasing bias between the true error and two ME curves starting from 16:30 is visible. Additionally, the fluctuation between 15:30 and 16:00 is not properly represented. There might be multiple reasons to explain these deviations. One of them is the presence of clock instabilities on the receiver or transmitter side or in both of them. In case the GPS-stabilized clock loses the lock on the satellites, its reference frequency will start drifting away inducing a bias. Nevertheless, given Fig. 6 Range error (blue), AGDF (red dotted), AGDF+ME prediction (black dash-dot), SM (orange dashed), and SM + ME SM (green dash-dot) predictions for the backward path the high grade and quality of the used clock, this is unlikely to be the reason. Further explanations include the possible changes of weather conditions along the path. We know that temperature, pressure, and humidity changes can produce a change in the atmosphere but also in the ground conductivity, especially in case of rains. Also, possible biases might be introduced by receiver or transmitter chain hardware components. Determining the exact cause is a challenging task and a clear answer cannot be provided to the reader due to the lack of information.
To assess the performance of the different techniques, we evaluate the absolute value of the range error and in particular its 95% ( 95 ) and maximum ( M ) value. If the range is used without any correction strategy, M is 72.4 m, as it can be seen in Fig. 6 whereas 95 is 58.9 m. With the SM, the improvements are not significant whereas with the AGDF, a maximum error of 18.4 m and a 95% error of 13.9 m are obtained. Finally, the best results are obtained with the use of the correction functions. When f ME SM (α) is applied, 95 is 7.4 m and M is 10.4 m, whereas by applying f ME (α), we obtain 7.8 m for the 95% error and 10.9 m for the maximum error.
To further validate the methodology, we decided to test the ME functions, generated with the forward path on 20th February, on a second dataset gathered on 24th February. This dataset is significantly larger than the previous one, lasting for approximately 7 h and covering the entire working maneuvers of the ship, as it is visible in Fig. 3. In a similar way to the analysis conducted for the backward path scenario, Fig. 7 presents the R-Mode range error and all the prediction curves with the same line style used in Fig. 6. All the curves start in zero due to the performed calibration, as previously explained. It can be noted that the plot can be divided in three main sections based on the curves' trend. The first section spans from 7:00 to 8:50, a time in which the vessel was moving from the shore to the open sea to conduct its activities. By comparing this first section with Fig. 4, it can be observed that the they all have a decreasing trend. This is indeed expected as the ship was moving in a similar direction. The only exception is the SM prediction, which increases. This can be explained by the fact that the distance is increasing over time; therefore, the predicted bias, which is only based on the constant velocity assumption over sea, increases as well. It is interesting to see that with the application of the ME functions, the variations of the range error are described very well. On the contrary, with the prediction based on the AGDF maps, only the first and the last parts are well fitted whereas a substantial bias is observed in the central part between 7:10 and 7:50.
The second section, from 8:50 to 12:40 is, in general, characterized by a nearly constant range error. In this time slot, the vessel was conducting its activity in a limited area. Therefore, the azimuth angle and the distance were almost constant. Excluding the SM curve, a good overlap of the lines is visible between 8:50 and 11:10, and afterwards, a deviation starts to appear reaching a maximum absolute value difference of about 8 m at around 12:00 and then recovering back at 12:40. As on February 20, the reason for this deviation is unknown. Additionally, we can identify a clear jump event in the rage error at 12:26 lasting for roughly 3 min. Other jumps can be identified in the data which are sometimes difficult to be seen depending on the magnitude of the jump itself. In general, the jumps can be of several meters and the duration can be of several tens of minutes in the wort cases. This is a know issue which affects some transmitters in a non-deterministic way.
Last but not least, the third section, from 12:40 until the end, represents the results obtained with the ship moving southwards. Therefore, an opposite trend is expected to be observed with respect to the first section. This is clearly visible in Fig. 7. As for the first section, the ME functions significantly improve the accuracy of the bias prediction when compared with the AGDF only case. Also, for this section, the SM follows an opposite decreasing trend. In fact, the distance between receiver and transmitter tends to decrease as the ship moved from the sea to the shore.
The analysis showed the following ranging performance of the different error reduction approaches for the second dataset. In the case of range measurements without correction schemes, a maximum error of 43.9 m and a 95% error of 34.5 m are obtained. If the AGDF maps are used, 95 is 16.6 m and M is 26.7 m. The accuracy is further increased when f ME (α) is applied, reducing 95 to 8.0 m and M to 16.7 m. On the contrary, if the SM is used, the accuracy slightly decreases. The 95% error increases to 37.2 m whereas the maximum error rises to 46.7 m. The performance improves if the SM prediction is combined with the ME SM function. In such a case, 95 drops down to 8.2 m while M becomes 17.3 m.
Finally, all the obtained results are summarized in the Table 2.

Operational implications
R-Mode is designed as a maritime navigational backup system which uses existing radio communication systems to broadcast ranging signals in the service areas of the legacy service. In the MF frequency band, the maritime radio beacons can be adapted to the R-Mode concept that they can provide R-Mode signals in a typical range of about 250 km. This is suitable to provide alternative position information especially in congested areas next to the coastlines. The requirements for a backup system are given in the IALA Recommendation R-129 (IALA 2012). They are structured according to different manoeuvres and risk of accidents. For R-Mode of special interest are the coastal waters and port approaches. Here R-129 defines a requirement of 100 m (coastal waters) and 10 m (port approaches) for the horizontal positioning accuracy with a 95% probability. A challenge for R-Mode is that the distribution of the maritime radio beacon was optimized to communicate with the vessels. This implies at least single coverage. For R-Mode, we are using the overlapping service areas of the stations. Depending on the region, the geometry of the transmitter sites is not optimal which could cause higher variances for the position. We expect an R-Mode positioning accuracy of about 1 to 2 times the ranging accuracy for the Baltic Sea. Therefore, accurate ranging is a precondition to fulfil IALA requirements. The paper presented two approaches to perform improved ranging. For the first, a simple model for the only sea signal propagation is used. With the knowledge of the average electrical parameters of the sea water, the delay of the signal over the sea path can be described. In addition, with a measurement campaign, the azimuthal error caused by the land, before the signal hits the sea water, can be estimated by sailing in a circle around the transmitter. This approach should be easy to implement by each service provider for new R-Mode transmitters. The service provider may provide the average sea parameter and azimuthal correction function to the user as an additional service.
The second approach is using AGDF maps, which are based on a model of the electrical properties of the ground at land and sea. The advantage of this approach is that it can already improve ranging without additional measurements as shown in Fig. 6. This approach is suitable for the implementation into an R-Mode receiver if no other information is available. Due to the fact that the land path is sometimes very inhomogeneous, the generation of accurate AGDF maps is not possible. In that case, the measurement of azimuthal error can be used to compensate AGDF map inadequacies. The measurement of that error can be done by the R-Mode service provider for a well-defined AGDF map or by the vessel itself if GNSS is usable.
Shortcomings of the prediction such as the spatial divergence of the true ground conductivity and permittivity values from the static map can be overcome with direct measurements. Also, seasonal variation can easily be captured by AGDF maps and it was shown that, on a time scale of a few days, the measurement-based corrections are stable over space and time.
The results in Table 2 clearly show that the simple sea path model will cause large errors. It works only in combination with the azimuthal correction function. Here a ranging of better than 9 m (95%) could be achieved for both days. With the AGDF maps, the ranging performance was better than 17 m (95%). It could be further reduced to better than 9 m (95%) when using the azimuthal correction function. Based on the analyzed data, we conclude the simple sea path model is not suitable to fulfil the requirements for an electronic backup navigation support for the coastal waters because position accuracies of worse than 100 m can be expected. The AGDF approach seems to meet the 100-m position accuracy requirement for coastal waters. When using azimuthal corrections, we are getting near to the IALA 10-m requirement for port approaches. With a good station geometry, the support of port approaches seems to be feasible.

Conclusion
In this paper, we proposed an innovative approach to increase the accuracy of the MF R-Mode ranging exploiting direct in-field measurements and GNSS positioning information in addition to a simple over sea path and a detailed over sea and land path error modelling approach. The latter is referred to as AGDF. After briefly describing the AGDF and how the signal propagation is affected by the electrical properties of the ground, a model of the range error was presented along with the new proposed methodology to estimate such a model error and generate a correction function dependent on the azimuth angle. Two correction functions were generated accordingly to two different model assumptions. In the first case, only the sea model was considered whereas in the second one, AGDF maps were used. The two approaches have been applied on real data in order to show that they can provide improvements with respect to the sole use of AGDF, which are based on ground conductivity maps. In the presented scenarios, the absolute range error was evaluated and its 95% and maximum values were taken as main performance indicator. In general, with the use of the azimuthal correction functions, an accuracy better then 9 m was achieved whereas the maximum error bound resulted to be 17.3 m. This level of improvement opens the door to the usability of the system where the 10-m horizontal positioning accuracy needs to be met.
The main advantage of azimuthal R-Mode correction functions remains in the fact that only small-scale measurement campaigns are needed to characterize the signal delay, generating a correction function dependent only on the azimuth angle. In principle, only one-way path survey covering all the azimuth angles from the transmitter would be sufficient; therefore, the costs and the time of the survey campaign can be remarkably reduced. Nevertheless, the assumption of no additional land between the coast line and the point of application must be fulfilled in order to obtain accurate ranging. Otherwise, the correction can still be applied but with an expected reduced R-Mode ranging accuracy.
Improving the quality of the AGDF maps might be seen as an alternative way to improve the ranging performance. This different approach would be optimal, allowing to use well-known prediction formulas but has two main challenges. First of all, the generation of high accuracy maps requires a deep and consuming surveying campaign. Secondly, some differences between the estimated value and the reality might remain due to local effects which are of extreme complexity to be modelled.
The results shown in this paper are of particular interest for the national maritime service providers since we demonstrated that the MF signal delay can be directly measured and therefore, a correction service can be established, which provides additional correction information to the user to improve the overall R-Mode service performance.