Positioning with medium frequency R-Mode

2 PARAMETRIC ESTIMATION

In this section, we introduce two methods for estimating the phase of a continuous wave (CW). We need a signal model for this. For a single tone, this can be described as:

1

with the three parameters frequency 𝜔_𝑐, amplitude 𝑏₀, and phase 𝜃.

It is a major task in signal processing to deliver good estimation for this parameter. We consider the frequency as known. However the amplitude and especially the phase have to be estimated in the medium-frequency range. As we are interested in a discrete sampled signal, the time 𝑡 can be quantized to 𝑡_n with constant sample rate 𝑡 as:

2

In this case, the bound can be given by:

3

for a complex signal model (Rife & Boorstyn, 1974). Here is the estimated phase parameter, 𝜎² is the variance in the power density function (PDF) of the noise, and 𝑁 is the total number of samples.

To obtain phase estimation, we suggest using a maximum likelihood estimation (Rife & Boorstyn, 1976). The likelihood function is given as:

4

with:

5

where 𝑋 is defined as the input vector of samples and 𝑌 as the corresponding Hilbert transform. Both vectors have the length 𝑁 and are indexed with 𝑛. 𝑘 is the number of tones to estimate, with each tone defined according to our signal model (1), hereby 𝑖 indicates the current frequency to evaluate within the sum. To find the parameters that maximize 𝐿, we need to maximize 𝐴(𝜔). Therefore, we obtain our estimates by:

6

and:

7

We assume that we know the frequency 𝜔_𝑐 of the tone and obtain the estimates for the phase and for the amplitude. To estimate the range, we are interested in the phase estimates for different tones in parallel, so the presented approach does not scale very well and leads to relatively high computational effort.

However, it is widely known that a discrete Fourier transform (DFT) can be used for a parametric estimation (Kay, 1993; Rife & Boorstyn, 1976; Stoica et al., 1989), since we can describe 𝐴(𝜔) from (4) using the DFT base function:

8

with:

9

where 𝐾 describes the point in the DFT spectrum. This yields to the matrix representation:

10

which is one column of the transformation matrix of the DFT. With the fast Fourier transform (FFT), a computational-effective method is known to estimate 𝑁 >bins at once. The center frequencies of these bins are equally distributed at intervals of Δ𝜔 and need to match the continuous wave frequencies that we want to estimate. This otherwise leads to power leakages and reduced accuracy of the results. This is the main difference of the maximum likelihood estimate which can be calculated for any arbitrary frequency but is limited to one frequency per calculation.

Both estimation techniques result in complex numbers where the phase represents the needed phase information at the beginning of the observed time interval with length 𝑇. Therefore the positioning solution would generated with at least T s old data.

3 SYSTEM BIAS AND PROPAGATION DELAY

After the signal has been generated at the transmitter site, it has to propagate to the receiver and pass different system components. Each of them will cause a delay or a distortion of the R-Mode signal.

Figure 2 shows the signal path as a generalized block diagram. The R-Mode signal is generated in a well-synchronized R-Mode modulator which provides the aiding carriers with a zero crossing (rising) at full seconds. From the modulator, the signal is passed to the transmitter followed by the antenna tuning unit (ATU) and the antenna as the last part of the transmitter chain. It was shown that each of these elements cause a phase offset. The characterization of these effects is the subject of ongoing research (Grundhöfer et al., 2019; Grundhöfer & Gewies, 2020).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2

Block diagram of medium-frequency signal propagation

After transmission, the signal propagates as a ground-wave to the receiver where it is received by an antenna (see Figure 2). From there, it is passed to an active band-pass filter with adjustable gain and, finally, to the software-defined radio (SDR). As for the transmitter site, each of the components introduces a certain phase shift.

In addition, the signal propagates along a second path towards the ionosphere where it is attenuated significantly by the E-layer of the ionosphere at daytime. At night, this absorbing layer disappears and the signal is reflected at higher layers. This propagation path is called a sky-wave. At the receiver site, the ground-wave and sky-wave superimpose leading to a degenerate phase estimate at night (Johnson & Swaszek, 2014c). Currently, different mitigation techniques are under discussion to make the service available during the night (Grundhöfer et al., 2020). Solving this problem is beyond the scope of this paper.

In the following, it is assumed that all system components at the transmitter and receiver site are stable with respect to the introduced phase offsets. Summing over these offsets will result in a transmitter, receiver site, and aiding carrier-dependent system bias 𝜃_system,𝑖. The index 𝑖 denotes in the following the aiding carrier 𝑖. The propagation causes a time delay 𝑡_p,𝑖 of the carrier. Both effects together cause an accumulated phase delay at the receiver of 𝜃_system,𝑖 +𝜔_𝑖𝑡_p,𝑖 for the carrier with angular frequency 𝜔_𝑖.

Due to the periodic nature of trigonometric functions, we have no access to the accumulated phase at the receiver site for the signal model in (1). The difference between the transmitted and received signal is 𝜃_𝑖 is reduced by a number 𝑚_𝑖 of completed cycles 2𝜋 since transmission. This leads to ambiguities in the estimation of the range, which can be described in terms of the system bias and the propagation delay as:

11

Both, 𝑡_p,𝑖 and 𝑚_𝑖 depend on the carrier and the distance to the transmitter.

The system bias of the transmitter sites in the R-Mode Baltic testbed is unknown. Therefore, this has to be estimated during a calibration step in the field. This step will also provide the number 𝑚_𝑖. With 𝜃_cal,𝑖 as a calibration parameter, we can deduce 𝑡_p,𝑖 as:

12

with:

13

This approach works well as long as the rover does not change position so much that 𝜃_𝑖 reaches 0 or 2𝜋. To prevent unrealistic jumps of 𝑡_p,𝑖, the phase 𝜃_𝑖 has to be tracked after fixing 𝑚_𝑖 so that it can grow above 2𝜋 or decrease below 0.

The pseudorange 𝑃_𝑖 between the rover and the transmitter site can be calculated using the phase velocity 𝑐_ground,𝑖 of the ground-wave with:

14

However, as the R-Mode signal propagates as a ground-wave, the phase velocity depends on the frequency and conductivity of the ground (Guard, 1992; Wait, 1998). For propagation over sea water we assume a phase velocity of 299,616,913 m/s based on calculation done by GRWave (Rotheram, 1981). This value is valid for our test area, as we take the conductivity of the Baltic Sea as a reference. Due to a different ground conductivity of land, the phase velocity differs for a path over land.

If the propagation path is composed of different segments of land and sea paths then (14) is obviously wrong. But if we assume there is little or no change in the propagation speed over land, we can use (14) as an approximation and put the nearly constant error as an additional bias into the calibration system bias 𝜃_system,𝑖. Another bias occurs due to the movement of the vessel. As the estimation averages over time and the distance varies over time, our estimation has a velocity-dependent bias.

4 LEAST SQUARES

As described in Section 3, the signal propagates as a ground-wave. Consequently, the geometrical distance between transmitter and receiver cannot be described as a line-of-sight distance (direct and straight line) which is given by the classical Euclidean distance formula. For the MF R-Mode system, the geometrical distance is given by the so-called geodesic curve or simply geodesic (Jekeli, 2016) which describes the minimum distance between two points on the Earth, assumed to be a perfect ellipsoid. The WGS84 is used as a reference system in order to localize the vessels on the sea. Given the non-linearity of the problem, a least squares (LS) adjustment algorithm that is used to estimate the receiver position and clock bias as described in Groves (2008).

The system of equations contains three unknowns represented by the user geodetic latitude 𝜙_𝑢, the geodetic longitude 𝜆_𝑢, and the receiver clock bias 𝛿𝑡 with respect to the R-Mode system time. The pseudorange 𝑃 that is obtained as an output from the phase tracker, which also considers the calibration bias, can be modeled as follows:

15

where 𝑖 represents the station index, 𝜌 is the geometrical distance (geodesic), and 𝜖 is a noise term. The geodesic 𝜌 can be computed by using Vincenty’s approach which is based on an iterative algorithm (Vincenty, 1975). It also provides estimates for the forward and backward bearing angles that can be used to directly define the geometry matrix.

In general, for an over-determined system with 𝑚 observations, the LS solution is (Groves, 2008; Subirana et al., 2013):

16

where 𝑙 is the algorithm recursion index. 𝐱_𝑙 represents the update estimate (𝜙_𝑢𝑙, 𝜆_𝑢𝑙, 𝛿𝑡_𝑙)^T, whereas 𝐱_𝑙−1 is the previous estimate (𝜙_𝑢𝑙−1,𝜆_𝑢𝑙−1,𝛿𝑡_𝑙−1)^T of the unknowns. 𝐇 is the geometry matrix and 𝐖 is a weighting matrix whereas 𝐌 ∈ ℝ^3×3 represents a transformation matrix that can be obtained as follows (Groves, 2008):

17

with 𝑅_𝑀(𝜙) meridional or north-south curvature radius, 𝑅_𝐸(𝜙) transverse or east-west curvature radius, and 𝜙′ angle subtended by an east-west great circle that subtends a unit longitude change. These values depend on the latitude of the user and are calculated as follows (Groves, 2008):

18

19

20

where 𝑎 and 𝑒 are the semi-major axis and eccentricity of the WGS84 ellipsoid.

In our approach, the geometry matrix 𝐇 ∈ ℝ^𝑚×3 containing the unitary line-of-sight vector components in the navigation frame can be obtained exploiting the bearing angle 𝛼, provided by the Vincenty formula, as (RTCM Special Commitee 127, 2016):

21

In general, the weighting matrix 𝐖 ∈ ℝ^𝑚×𝑚 takes into account the difference in the data quality (Subirana et al., 2013). Since the characterization of the measurements error has not been performed for the MF R-Mode transmitters, the weighting matrix is assumed to be a diagonal matrix where all the stations are treated equally (i.e., all values on the diagonal are set to 1). Last but not least, we define 𝐲_𝑙−1 as pre-fit residuals obtained as follows:

22

where represents the geodesic for the previous estimated point to 𝑖-th station.

The algorithm presented here was implemented in Python and integrated into the processing which is described in more detail in Section 5.2.

5 TEST SETUP

5.1 Hardware

The current setup for the R-Mode positioning test is comprised of two main components: an R-Mode receiver developed in the German Aerospace Center (DLR; Grundhöfer &Gewies, 2018), which is able to track all MF R-Mode signals of transmitters in the service area where it is currently located, and a Saab R5 GNSS receiver, which makes use of the Atlas correction service, Hemisphere, to achieve a position accuracy of 8 cm using a real-time kinematic (RTK) positioning approach. The R5 receiver worked during the trial in the position fixed mode such that it provided accurate positioning for the calibration of the R-Mode receiver and for the reference trajectory for later data analysis.

The current hardware evolution of the DLR R-Mode research receiver utilizes off-the-shelf components to create an easy-to-replicate receiver design. Figure 3 shows the wiring of the components as a block diagram. In the first stage, an active bandpass filter is used as an anti-aliasing filter, which also amplifies the signal by 30 dB. In our setup, we used a Krohn-Hite FBM3004 with a passband of between 285 kHz to 325 kHz.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3

Block diagram of MF R-Mode receiver

The heart of the setup is built around the Ettus X310 SDR (Ettus Research, 2019) with a LFRX daughterboard. It is important to notice that this SDR is designed for differential-ended or quadrature-modulated signals.

The front end described here is single-ended with coax cable as interconnections. Therefore, each daughterboard supports up to two separated front ends, which can be used to compare different antennas. The sampled data of SDR will be processed at the computer as described in Section 5.2.

For performance tests of the R-Mode system, we utilized an external timing source, an LL-3760 GPS timing receiver from Lange-Electronic, that provides a pulse-per-second (PPS) and 10-MHz signal as a time reference to the SDR.

As an antenna, we use an E-field antenna called Boni-Whip, which is a short dipole with active matching to the desired frequency band.

5.2 Software

The software for processing the obtained sample is organized into four main tasks. Figure 4 shows an overview of how each task interacts with the following one. Therefore, each routine has a name, an input interface, and an output interface, utilizing TCP sockets and ZeroMQ.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

Block diagram of signal processing chain

In the first step we have the SDR interface. This routine controls the SDR and its functionality. All the main settings for wiring the time source and antennas need to be defined here. Moreover, it is necessary to choose the desired channel and sample rate.

The interface hands over the samples to the right phase estimator. To simplify the processing, this routine creates chunks of samples for a period of 1 s. Each chunk is timestamped on the first sample. The size of the chunk is chosen according to the minimum observation time required in the receiver.

The chunks are handed via a TCP port into the phase estimation routine. Here, the phase is determined with a parametric estimation. To deliver a phase estimation near to real time, the execution time of any estimate is crucial. Therefore, we decided to implement the FFT as phase estimator on our current receiver design. According to Section 2, this corresponds to a maximum likelihood estimation with a bandwidth of 1 Hz for 1 s of observation time.

The cyclic nature of the R-Mode signal with its wave-length of about 1 km leads to jumps in the estimated phase after a few minutes of sailing. These jumps are then seen in the form of ambiguities in distance estimates. This complicates data processing when we want to reuse information from earlier epochs. Therefore, we implemented a phase tracking algorithm that compensates any phase jumps over the upper or lower phase limit of our estimator.

It adds multiples of 2𝜋 to the phase estimate so that the resulting phase never changes more then 𝜋 from one epoch to the next. So the tracked phase starts between 0 and 2𝜋 but can reach very large values after hours of continuous signal tracking.

The tracked phase needs to be corrected for any bias described in Section 3. This bias is obtained during the calibration step. Here, the difference of the expected phase, which is calculated based on the known own position (e.g., from GNSS), and the measured phase is averaged over a defined time and taken as a bias. It is different for every station and every continuous wave component. For the R-Mode measurement on land, a similar approach was proposed (Johnson et al., 2020), but they assumed a fixed position. The obtained biases are used to correct all future estimations.

After the calibration, GNSS data are no longer required, as we expect that the bias is constant during our measurements. The speed of light for propagation over sea water is used to calculate a range, which is handed over in the next processing routine. This is done by using a publisher subscriber pattern from zeroMQ (Sústrik, 2020), as the execution time is crucial to prevent memory overflow and the ranges drastically reduce memory consumption. The pattern works as a buffer and enables distribution of the ranges to different positioning algorithms in the future.

The positioning routine implements the least squares algorithm presented in Section 4. The obtained position and clock error is saved with a timestamp in a comma-separated values (CSV) file and called a first visualization. This file can be read later on by the visualization routine to examine the data more closely and present it in an attractive way.

6 RESULTS

In August 2020, we conducted a measurement campaign in the R-Mode Baltic testbed with the Fyrbyggaren, a buoy-tender vessel provided by the SMA. As described in Section 5, the DLR R-Mode receiver and the Saab R5 GNSS RTK receiver were set up on the ship. The DLR receiver was continuously tracking the signals of the MF R-Mode stations with the best signal quality while the Saab R5 receiver continuously provided a highly accurate ship position which was used for real-time calibration of the R-Mode receiver and as a reference for the performance estimation of R-Mode. Furthermore, the DLR receiver was able to perform R-Mode real-time positioning at sea. The results of two data sets, where the ship was sailing through the center of the testbed, are presented in this paper.

The track of the Fyrbyggaren for both data sets are shown in Figure 5 together with the R-Mode transmitter sites of the three strongest signals (Holmsjö, Hoburg, and Rozewie) in that region. The first data set is surrounded by a white rectangle and the second by a yellow rectangle. An inset shows in detail the deviation of the R-Mode-based track in blue from the RTK track in red. There is a gap between both data sets because we performed experiments with a different sampling rate over some hours.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

Track comparison between R-Mode and RTK (The map data was provided by OpenStreetMap©contributors [Online])

During the measurements, the ship was sailing for several hours between the three transmitters, not far away from the center between them. The horizontal dilution of precision (HDOP) was better than 1.7 for all times which is a good precondition to obtain accurate position estimations.

The first data set was logged for roughly 90 minutes in the morning, about 1.5 hours after sunrise, on August 31. Clearly visible in the inset of Figure 5, the R-Mode position follows the RTK position with some systematic scatter southwards as the ship was sailing at different speeds and the sea waves were small.

The second data set was logged on the same day for slightly less than five hours after noon. In contrast to the first data set, the ship sailed faster and had stronger waves from the direction we were sailing. Figure 5 shows that the propagation path of the signals from the Hoburg and Rozewie stations stayed most of the time over sea, whereas for the Holmsjö signal, some variations of the signal delay was expected because of the changing length of the land path component.

In both cases, the three transmitters of Holmsjö, Hoburg, and Rozewie were used for positioning with the LS approach. Groß Mohrdorf, a fourth station could also be received (see Figure 6), but due to the smaller signal-to-noise ratio (SNR) and the fact that we do not weight the stations in our current approach, the position became worse when including Groß Mohrdorf.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6

Overview of measured frequency spectrum at 5:50 UTC

At the beginning of both data sets, an initialization over 30 s took place. For data set one, this was a quasi-static calibration and for data set two, a dynamic calibration (speed of about 10 kn superimposed with movement of the ship due to sea waves) to estimate the biases in the system of transmitter, receiver, and propagation path for each CW of all R-Mode transmitters in view. To reduce the impact of noisy phase estimations of one CW on the position, we always considered both CW per transmitter site in the position estimation. We performed positioning with six tracked signals of the three nearest R-Mode transmitters.

To evaluate the positioning performance of the R-Mode receiver, the horizontal distance between the R-Mode and the reference position was calculated and considered as horizontal error. Of particular relevance were the maximum (𝜖_𝑀), mean and 95 % bound (𝜖₉₅) of the horizontal error which were used for performance assessment.

It is important to mention that the antennas of the R-Mode and GNSS receiver were mounted on a horizontal traverse at a distance of 25 cm on top of the ship’s mast. This distance has been neglected in the initialization step and in the performance evaluation. Therefore, a horizontal error of up to 50 cm can be expected when the ship changes the heading after initialization.

6.1 First data set

The first data set combines quasi-static and dynamic measurements. In particular, at the beginning of the data set, the ship was supposed to be in a static condition. However, some movements were caused by the presence of low wind speed and small waves. For this reason, this scenario is referred to as quasi-static. The second part of the measurements was conducted with 5 kn speed, a scenario referred to as low speed whereas the last part was conducted with 10 kn speed and referred to as a medium-speed condition. The difference between the three aforementioned cases can be seen in Figure 7(b) that presents the speed of the ship obtained from the GNSS receiver.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 7

(a) Horizontal position error and (b) RTK-based speed of the ship over the time for the first data set

Figure 7(a) shows the horizontal error of the R-Mode receiver over the entire data set. As expected, it coincides with the R-Mode and RTK track of Figure 5 where both are very close together. It appears that most of the epochs are characterized by an error significantly smaller than 10 m but the presence of repetitive spikes that overtake this value can be noted. A few outliers with more than 20 m are visible.

Figure 8 illustrates the pseudorange error (the difference between the R-Mode pseudorange and RTK-based calculated range) for the three used stations and for both CWs, designated as CW1 (low-frequency carrier) and CW2 (high-frequency carrier). The presence of fast variations which represent the noise on the pseudorange residual can be noted for the three stations. This is significantly smaller than 10 m most of the time.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 8

Pseudorange residuals of the first data set for (a) Holmsjö, (b) Hoburg, and (c) Rozewie

On top of the noise, a small time-varying bias can be observed for all stations, which is most likely caused by the changes in the signal propagation path over land, which has so far not been modeled in the receiver (see Section 3). This assumption is justified because the bias on time scales of several minutes is correlated for the CWs of one station. This is expected because changes in the propagation path affect both signals very similarly due to their nearly identical frequencies. On the other hand, the biases are different among different stations due to the fact that the signals travel along different paths and therefore are affected by distinct effects.

The pseudorange error for Holmsjö in Figure 8(a) provides the explanation for the spikes in the position error in Figure 7(a). It shows the same pattern. Every minute, the pseudorange error of both CWs increases for 7 s by 10 m to 15 m. Unfortunately, an explanation for this issue has not been found.

The 95% bound for the positioning error (𝜖₉₅) was mainly driven by the spikes of the Holmsjö pseudorange error. This caused a significant error in about 12% of the epochs and was independent of the speed of the vessel. 𝜖₉₅ is very similar for the three phases of quasi-static, low speed, and high speed with 10.3 m, 7.5 m, and 7.6 m, respectively.

The decrease of 𝜖₉₅ for the dynamic phases can be explained with a decrease of the Holmsjö pseudorange error (see Figure 8) caused by the superposition of the spikes and a negative bias for the time after 5:39 UTC (see Figure 8), which was caused by changes in the propagation path over land. For the entire data set, 𝜖₉₅ is 8.4 m, which is better than the requirement for navigation support of port approaches. Fixing the Holmsjö issue should considerably improve the performance.

The mean error is a better measure for the R-Mode performance in this case, because it is not only affected by the Holmsjö issue. It reflects more than the noise of the pseudorange estimations. With 3.9 m, 3.4 m, and 3.5 m for the quasi-static, low-speed, and high-speed phases, and a value of 3.7 m for the entire data set, it is considerably lower than 𝜖₉₅.

The maximum error 𝜖_𝑀 reflects the outliers in the R-Mode positioning. It is 18.1 m for the quasi-static, 44.7 m for the low speed, and 27.1 m for the medium-speed phase of the experiment. The outliers are caused by jumps in the phase of the CW for a single station [see Figure 7(a) and Figure 8]. Such errors might be suppressed when pseudoranges with unrealistically large jumps are excluded from positioning.

6.2 Second Data set

The second data set refers to a purely dynamic test in which the vessel was traveling with a speed of 10 kn and was experiencing adverse weather conditions with strong waves. This is visible in Figure 9(b), which shows the speed of the vessel obtained with the GNSS receiver. A higher variance of measurement, induced by the heavy sea-waves, can be observed when compared with Figure 7(b) where the same scale is adopted.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 9

(a) Horizontal position error and (b) RTK-based speed of the ship over the time for the second data set

Figure 9(a) illustrates the R-Mode horizontal position error. An error smaller than 20 m characterizes most of the epochs and the presence of the repetitive spikes can also be observed in this case. Compared to 7(a), a time-varying bias can be noted which causes a larger error from 10:20 to 11:15 and from 13:00 up until the end of the experiment. This strong bias is mainly induced by Holmsjö, as we see later.

Similarly to the first data set, Figure 10 presents the pseudorange residuals for the CWs of the transmitters. Hoburg [see 10(b)] was very stable during the measurements with a final small negative bias of approximately 5m. Rozewie [depicted in Figure 10(c)] was also very stable but the noise level was slightly higher compared to Hoburg. As explained previously, the stability of these two stations comes from the fact that the signal traveled mostly over the sea.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 10

Pseudorange residuals of the second data set for (a) Holmsjö, (b) Hoburg, and (c) Rozewie

On the contrary, for Holmsjö [Figure 10(a)], the bias was strongly influenced by the land path component which was fluctuating severely during the movement of the vessel and induced an offset increase from 9:00 to 11:00 UTC and a decrease from 11:00 to 13:30 UTC. The repetitive spikes are clearly visible.

The 95% horizontal error was driven primarily by the spikes, as explained in Section 6.1. In this experiment 𝜖₉₅ was 12.0 m, which is slightly higher than the first data set value, but it is mainly due to the effect of the bias caused by Holmsjö. For the same reason, the mean horizontal error was–at 5.2 m–also larger. These results were expected due to the larger variation in the ship’s position during the 5 h time period and the uncompensated effect of the ground conductivity changes.

The maximum horizontal error 𝜖_𝑀, which can represent outliers, was 20.2 m. Compared to the first data set, jumps were not visible for all transmitters (see Figure 10) and therefore no considerable outliers could be identified.

HOW TO CITE THIS ARTICLE

Grundhöfer L, Rizzi FG, Gewies S, et al. Positioning with medium frequency R-Mode. NAVIGATION.2021;68:829–841. https://doi.org/10.1002/navi.450