Theoretical Concept for a Mobile Underwater Radio-Navigation System Using Pseudolite Buoys

: Reliable and accurate positioning underwater has been the goal of many research activities in the last years. In this paper, we present our design of a self-maintainable underwater navigation system, which provides a simple, low-cost, and ﬂexible solution. In our concept, a set of buoys is equipped with Global Navigation Satellite System (GNSS) receivers as well as underwater pseudolite signal transmitters. Therefore, the buoys are able estimate their position and transmit this information underwater through a pseudolite signal. We particularly designed this signal for underwater ranging for both fresh and salt water. Any suitable underwater receiver is then able to do ranging via these pseudolite signals. In this paper, we theoretically analyze and test our concept using several buoy distributions and diving depths. The optimal arrangement of buoys is determined in terms of accuracy and availability depending on the number of available buoys and targeted water depth for an e ﬃ cient operation. For example, at a targeted depth of 30 m in fresh water, a maximum horizontal position root-mean-square (RMS) error of less than 3 m can be achieved with a set of ﬁve buoys providing a service radius of up to 72 m.


Introduction
A reliable and accurate underwater positioning is mandatory for divers, as well as underwater vehicles [1]. In the past, many different approaches have been proposed and most of them are based on acoustic underwater positioning system [2][3][4][5] or using optical signals [6]. Acoustic waves are mostly used for underwater ranging in salt water as electromagnetic waves are strongly attenuated. Only short ranges are achieved with electrometric waves compared to the ones obtained by acoustic waves, and hence, they are not suitable for underwater ranging methods [7]. A system for the tracking and navigation of underwater vehicles or divers through acoustic distance, or direction measurements, and subsequently by position triangulation can be found in [8]. The measurements are based on a sea-floor baseline transponder network, which limits the area of the system's operation. Depending on the baseline length, there are three broad types: Long Baseline (LBL) System, Ultra Short Baseline (USBL) System, and Short Baseline (SBL) System. Though acoustic waves have clear advantages in low underwater attenuation, they also have some significant drawbacks. They are complex, and many acoustic systems cannot be exported as they are subject to technological restrictions imposed by the respective governments. Acoustic systems are also sensitive to noise, and due to low diffraction suffer from shadow zones created by large objects or seafloor topography. Furthermore, acoustic links can be distancing of buoys in order to provide a maximum service area while maintaining an underwater position standard deviation of less than three meter. This discussion is done for several numbers of available buoys and water depths of operation. Finally, we outline the conclusions of our findings.

Materials and Methods
The concept of our mobile underwater radio-navigation system using pseudolite buoys is a new approach to solve the underwater navigation problem. The system can be split in two parts: the terrestrial localization of each buoy using GNSS, and the underwater localization applying triangulation on pseudolite signals transmitted by each buoy underwater. A scheme of its operation is presented in Figure 1.
The key device is a buoy equipped with a dual-frequency, multi-constellation GNSS antenna and receiver on the top and a pseudolite signal transmitter on the bottom. A set of at least three buoys creates the needed aquatic environment needed for our investigations. The buoy's GNSS receiver determines its position in real time based on dual-frequency, multi-constellation GNSS data. The required antenna and receiver can be COTS products. The goal of our system was to achieve underwater position accuracy within a standard deviation lower than 3 m. The developed navigation system should be portable and low cost. The underwater navigation system should preferably use electromagnetic waves instead of acoustic waves so that they should have a low impact on underwater mammals. At the very least, the rapid attenuation of electromagnetic waves in water prevents any harm to mammals which are outside the service volume of the proposed system. Nevertheless, the principles of the proposed concept can also be utilized for systems with higher accuracy requirements and acoustic signals.

Underwater Signal Design and Ranging
The signal transmitted by each buoy must essentially fulfill two tasks. First, it must allow precise measurement of the distance (range) it has traveled from the buoy to the user's antenna. Second, it should carry the necessary data bits, which provide information about the buoy's coordinates. Ranges and buoys' coordinates are combined by the user to determine its position in an absolute coordinate frame.
Signal propagation conditions of a medium as challenging as water are taken into account when designing a reliable navigation signal system. In the following, we discuss and selection of suitable The key device is a buoy equipped with a dual-frequency, multi-constellation GNSS antenna and receiver on the top and a pseudolite signal transmitter on the bottom. A set of at least three buoys creates the needed aquatic environment needed for our investigations. The buoy's GNSS receiver determines its position in real time based on dual-frequency, multi-constellation GNSS data. The required antenna and receiver can be COTS products.
The goal of our system was to achieve underwater position accuracy within a standard deviation lower than 3 m. The developed navigation system should be portable and low cost. The underwater navigation system should preferably use electromagnetic waves instead of acoustic waves so that they should have a low impact on underwater mammals. At the very least, the rapid attenuation of electromagnetic waves in water prevents any harm to mammals which are outside the service volume of the proposed system. Nevertheless, the principles of the proposed concept can also be utilized for systems with higher accuracy requirements and acoustic signals.

Underwater Signal Design and Ranging
The signal transmitted by each buoy must essentially fulfill two tasks. First, it must allow precise measurement of the distance (range) it has traveled from the buoy to the user's antenna. Second, it should carry the necessary data bits, which provide information about the buoy's coordinates. Ranges and buoys' coordinates are combined by the user to determine its position in an absolute coordinate frame.
Signal propagation conditions of a medium as challenging as water are taken into account when designing a reliable navigation signal system. In the following, we discuss and selection of suitable navigation frequencies in the Low Frequency (LF) / Ultra Low Frequency (ULF) bands. Our design takes into account the propagation of electromagnetic (EM) waves in water at various frequencies and the predominant noise sources in water based upon the comprehensive overview on underwater radio signals by Moore [20]. A brief review of the relevant propagation and noise characteristics is provided in the Appendices A and B of this paper.

Choice of Navigation Frequency and Bandwidth
The signal propagation characteristics given in Appendix A suggest that underwater EM waves are not an attractive option for use as they generally require large antenna systems, high power, and very low frequencies [16,17]. The need for usage of very low frequencies limit the selection of possible bandwidths and data transfer rates: since higher frequencies are attenuated significantly more than lower frequencies and the transmitted signal that does not have a very narrow bandwidth will arrive in a distorted form at the receiver. Clearly, a signal such as the GPS L1 C/A code containing all frequencies from 0-1 MHz could not be transmitted or retransmitted by a buoy as some frequencies would be attenuated significantly more than others, and the received signal would be an extremely distorted one. Instead, a very simple navigation signal with a narrow bandwidth should be targeted, which transmits only a few bits per second.
When choosing the preferable transmit frequency for a system such as ours, certain difficulties are faced: on the one hand, signals with short wavelengths are undesirable because of their strong attenuation. In a conducting medium, transmitting the signals over the distance of several wavelengths is not at all possible [16]. On the other hand, using a large wavelength may be a better option due to weaker attenuation, but larger antennas will be required. As an intermediate solution, we decided to use the moderate wavelength λ = 50.0 m (equivalently, a skin depth of δ = 7.96 m). Depending on the type of water, this corresponds to a frequency of f U = 1.0 kHz (sea water/ULF) and f V = 80.0 kHz (fresh water/LF). For these frequencies, attenuation is still acceptable (1.09 dB/m) and far-field assumptions such as shown in Equation (A1) can serve as a first-order approximation.

Pilot Signal and Ranging Accuracy
Ranging in GPS is classically achieved by estimating the navigation signal's time-delay, which it experienced when travelling through space. Time-delay estimation is done by detecting rising or falling signal edges of a known code sequence (pilot) and comparing them to a local reference. The accuracy of such code-delay-based ranging depends on the "sharpness" of the rising and falling signal edges, hence on the signal bandwidth. The approximate order of magnitude of achievable ranging accuracy is given by c/B, where c is the speed of light and B is the signal bandwidth. For the 1 MHz L1 C/A signal, this conforms to 300 m. Obviously, an underwater EM signal with bandwidth of the order of Hz or even kHz does not result in sufficient accuracy.
We propose ranging based on the signal amplitude instead of the signal delay as a solution. We exploit a simple principle: the lower the received signal power, the further away the signal source. Unlike in free space, the amplitude attenuation of a received waveform is a strong indicator of the range r it has travelled. Amplitude decreases only with r −1 in air or vacuum, but with r −1 e −r/δ in a conducting medium such as water, where δ is the medium's skin depth.
Consider the signal block diagram in Figure 2. The buoy's pilot signal is represented by Remote Sens. 2020, 12, 3636 5 of 23 with the effective transmit power P and ω = 2πf, with the technical frequency f as selected in Section 2.1.1. We identify the range-dependent attenuation and phase in Equation (A1): Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 23 The buoy's pilot signal is represented by with the effective transmit power and = 2 , with the technical frequency as selected in Section 2.1.1. We identify the range-dependent attenuation and phase in Equation (A1): Notably, the phase also carries information about the unknown range, which we shall exploit. Now, the signal received by an antenna with receiving cross-section (or effective aperture) 0 (in units of m 2 ) located at a distance from the transmitter can be expressed as: The unknown term ( ) is any kind of interference, which is mainly atmospheric noise. After passing through a low-noise amplifier (LNA), the signal is low-pass filtered with bandwidth and is converted from its carrier frequency to an in-phase/quadrature (I/Q) baseband signal. We assume that the complex-valued baseband signal is sampled with a high-resolution logarithmic quantization analog-to-digital (A/D) converter during an observation interval of length 0 ( , ) = ( ) ( ) + ( ), 0 < ≤ where ( ) is now a complex-valued zero-mean noise process. Typical noise levels in the relevant LF/ULF band are briefly discussed in Appendix B. For simplicity, we assume proper Gaussian atmospheric noise with power spectral density 0 We use the EXtended Invariance Principle (EXIP) of maximum likelihood estimation to estimate from ( , ) in a two-step approach. Notably, the phase also carries information about the unknown range, which we shall exploit. Now, the signal received by an antenna with receiving cross-section (or effective aperture) A 0 (in units of m 2 ) located at a distance r from the transmitter can be expressed as: The unknown term w(t) is any kind of interference, which is mainly atmospheric noise. After passing through a low-noise amplifier (LNA), the signal is low-pass filtered with bandwidth B and is converted from its carrier frequency to an in-phase/quadrature (I/Q) baseband signal. We assume that the complex-valued baseband signal is sampled with a high-resolution logarithmic quantization analog-to-digital (A/D) converter during an observation interval of length T 0 z(r, t) = PA 0 a(r)e jφ(r) + w(t), 0 < t ≤ T 0 (5) where w(t) is now a complex-valued zero-mean noise process. Typical noise levels in the relevant LF/ULF band are briefly discussed in Appendix B. For simplicity, we assume proper Gaussian atmospheric noise with power spectral density N 0 We use the EXtended Invariance Principle (EXIP) of maximum likelihood estimation to estimate r from (r, t) in a two-step approach.

2.
It is realized in digital signal processing by using an integrate-and-dump (I&D) filter. Now the range r fromγ is estimated. The maximum likelihood estimator is given by: A graphic interpretation of this approach is shown in Figure 3 in the form of a Nyquist plot (location curve) in logarithmic polar coordinates. The two outputs from the digital signal processing are shown in Figure 2: Re(γ) and Im(γ) can be mapped to a polar plot with phase and amplitude. The receiver will estimate the range depending on nature of the two outputs. Each traveled wavelength corresponds to one complete circle of 360 • , and the phase ambiguity can be resolved via the amplitude. As a noisy measurement will typically never synchronize with the shown curve, a criterion such as Equation (7) is necessary to obtain the range estimate. When many noisy measurements are available, they can be averaged to increase accuracy.
A graphic interpretation of this approach is shown in Figure 3 in the form of a Nyquist plot (location curve) in logarithmic polar coordinates. The two outputs from the digital signal processing are shown in Figure 2: Re( ) and ( ) can be mapped to a polar plot with phase and amplitude. The receiver will estimate the range depending on nature of the two outputs. Each traveled wavelength corresponds to one complete circle of 360°, and the phase ambiguity can be resolved via the amplitude. As a noisy measurement will typically never synchronize with the shown curve, a criterion such as Equation (7) is necessary to obtain the range estimate. When many noisy measurements are available, they can be averaged to increase accuracy. The EXIP states that for sufficiently long observation times 0 (long averaging), the above twostep approach has the same performance as a direct (but computationally costly) estimation of the range from the received signal. The variance of the range estimate is given by: This means that the ranging accuracy is dependent on the range itself: it is more difficult to get range estimation from a faraway buoy. The standard deviation ("accuracy") of ranging is given by ̂= √Var .
In Table 1 we suggest a set of meaningful/feasible system parameters and the accuracy is plotted as a function of distance in Figure 7.  The EXIP states that for sufficiently long observation times T 0 (long averaging), the above two-step approach has the same performance as a direct (but computationally costly) estimation of the range from the received signal. The variance of the range estimate is given by: This means that the ranging accuracy is dependent on the range itself: it is more difficult to get range estimation from a faraway buoy. The standard deviation ("accuracy") of ranging is given by σr = √ Varr. In Table 1 we suggest a set of meaningful/feasible system parameters and the accuracy is plotted as a function of distance in Figure 7.
It can be seen that good ranging accuracy (on the order of a few meters) is possible for ranges up to 90 m. For smaller ranges r < λ, the figure suggests that millimeter accuracy is possible. However, in this case, far-field conditions are clearly violated so that these results are overoptimistic and should be interpreted with care. However, we can state at this point that ranging with our system is theoretically feasible.
Finally, note from Equation (8) that a trade-off between transmit power and observation time can be made. The power P can be reduced at the cost of a longer time T; hence, a slower dynamic response of the system. To transmit the necessary bits containing the buoy's absolute position in the GPS coordinate frame, we transmit a second signal given by: Pilot and data signals do not interfere due to the orthogonality of sine and cosine. The phase of the sine wave alternates by ± π 2 with a rate of 1/T b , where T b is the bit duration, depending on whether a binary 0 or 1 is transmitted. The data signal and the necessary receiver's signal processing are illustrated in Figure 4. As the phase ∠γ = ϕ(r) is already known from the pilot module, the data can be coherently demodulated, without using a quadrature baseband channel. The I&D filter output β is used by a hard decision device to decide which datum was transmitted. In the case β < 0, we decide for "1", or otherwise we decide for "0". It can be seen that good ranging accuracy (on the order of a few meters) is possible for ranges up to 90 m. For smaller ranges < , the figure suggests that millimeter accuracy is possible. However, in this case, far-field conditions are clearly violated so that these results are overoptimistic and should be interpreted with care. However, we can state at this point that ranging with our system is theoretically feasible.
Finally, note from Equation (8) that a trade-off between transmit power and observation time can be made. The power can be reduced at the cost of a longer time ; hence, a slower dynamic response of the system.

Data Signal and Bit-Error-Probability
To transmit the necessary bits containing the buoy's absolute position in the GPS coordinate frame, we transmit a second signal given by: Pilot and data signals do not interfere due to the orthogonality of sine and cosine. The phase of the sine wave alternates by ± with a rate of 1/Tb, where Tb is the bit duration, depending on whether a binary 0 or 1 is transmitted. The data signal and the necessary receiver's signal processing are illustrated in Figure 4. Note that if the phase ∠γ = φ(r) is already known from the pilot module, and the data can be coherently demodulated; hence the receiver can save the quadrature baseband channel. The I&D filter output β is used by a hard decision device to decide which datum was transmitted. In the case β < 0, we decide for "1", or otherwise we decide for "0". From the energy per bit = we can obtain the following well-known expression for bit error probability with binary phase-shift keying (BPSK) modulation: Suppose that the buoy transmits two coordinates (in the horizontal plane), each consisting of 12 bits. Therefore, one frame (packet) of coordinates consists of bits. The frame error probability is just given by = 1 − (1 − ) . Both and are shown in Figure 5, where the data rate was chosen as 3 bits/s (i.e., the 2D coordinates are transmitted and updated every 8 s, which should be From the energy per bit E b = PT b we can obtain the following well-known expression for bit error probability with binary phase-shift keying (BPSK) modulation: Remote Sens. 2020, 12, 3636 8 of 23 Suppose that the buoy transmits two coordinates (in the horizontal plane), each consisting of 12 bits. Therefore, one frame (packet) of coordinates consists of N = 24 bits. The frame error probability is just given by P F = 1 − (1 − P e ) N . Both P e and P F are shown in Figure 5, where the data rate was chosen as 3 bits/s (i.e., the 2D coordinates are transmitted and updated every 8 s, which should be enough to cover the slow buoy dynamics). All remaining parameters are as specified in Table 1. Reliable data transmission is possible up to 85 m of range with such a data rate. Again, error probabilities for ranges on the order of the wavelength may be overoptimistic.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 25 Suppose that the buoy transmits two coordinates (in the horizontal plane), each consisting of 12 bits. Therefore, one frame (packet) of coordinates consists of bits. The frame error probability is just given by = 1 − (1 − ) . Both and are shown in Figure 5, where the data rate was chosen as 3 bits/s (i.e., the 2D coordinates are transmitted and updated every 8 s, which should be enough to cover the slow buoy dynamics). All remaining parameters are as specified in Table 1. Reliable data transmission is possible up to 85 m of range with such a data rate. Again, error probabilities for ranges on the order of the wavelength may be overoptimistic.

Multiple Access Strategy and Total System Bandwidth
So far, we have only considered one buoy. However, since we need multiple buoys to allow positioning in all three dimensions of space, we have to ensure that buoys do not interfere with each other. One possible strategy to avoid interference is frequency-division multiple access (FDMA) as illustrated in Figure 6.

Multiple Access Strategy and Total System Bandwidth
So far, we have only considered one buoy. However, since we need multiple buoys to allow positioning in all three dimensions of space, we have to ensure that buoys do not interfere with each other. One possible strategy to avoid interference is frequency-division multiple access (FDMA) as illustrated in Figure 6. Here, the buoys would not all use the same carrier frequency given in Equations (A4) and (A5), but rather transmit on separate frequencies 1, 2, 3, … adjacent to . Each buoy needs a bandwidth equal to the data rate 1/ plus a guard interval ± . Then, the total range of frequencies occupied by an ensemble of buoys is the frequency interval + , − , with the total system bandwidth: For instance, with = = 5Hz and four buoys, we would have a total system bandwidth of = 80 Hz.

Antennas for Underwater Communication and Frequency Allocation
The antenna length should be a fraction of the operating wavelength for an effective radiation. For the most of antennas it is half-wavelength (λ/2) or quarter wavelength (λ/4). As water has high dielectric constant and conductivity, the resulting wavelength will be much shorter than in the air for a given frequency (e.g., for an electromagnetic wave with a frequency of 10 kHz the wavelength is 30 km in space but only 15.8 m in seawater). However, as mentioned before, water is a conductive medium and to minimize very high attenuation loss of signals in water, very low frequencies (i.e., long wavelengths) must be used. A wavelength of 50 m was selected for the proposed navigation system, which results in the operating frequency of 1 kHz in seawater and 80 kHz in fresh water. For such a wavelength, even antennas with a length of a quarter of operating wavelength are too bulky for a portable system. Therefore, electrically small antennas must be used [21]. Such antennas are much shorter than the wavelength of the operating signal and are generally less efficient. They are also more challenging to design than larger full-sized antennas. Impedance matching is needed for electrically short antennas to be able to operate at 50 ohms and it can be achieved with the help of lump components. Many types of antennas for underwater applications are discussed in the literatures [22][23][24][25][26][27]. Most of them are wideband antennas for underwater communications in the MHz region; however, the basic principles of design and performance of these antennas apply also for lower frequencies in the kHz region [28]. These antennas are based mostly on loops, long wires, and Here, the buoys would not all use the same carrier frequency f given in Equations (A4) and (A5), but rather transmit on separate frequencies f 1 , f 2 , f 3 , . . . adjacent to f. Each buoy needs a bandwidth equal to the data rate 1/T b plus a guard interval ±F. Then, the total range of frequencies occupied by an ensemble of K buoys is the frequency interval f + B s 2 , f − B s 2 , with the total system bandwidth: For instance, with 1 T b = F = 5Hz and four buoys, we would have a total system bandwidth of B s = 80 Hz.

Antennas for Underwater Communication and Frequency Allocation
The antenna length should be a fraction of the operating wavelength for an effective radiation. For the most of antennas it is half-wavelength (λ/2) or quarter wavelength (λ/4). As water has high dielectric constant and conductivity, the resulting wavelength will be much shorter than in the air for a given frequency (e.g., for an electromagnetic wave with a frequency of 10 kHz the wavelength is 30 km in space but only 15.8 m in seawater). However, as mentioned before, water is a conductive medium and to minimize very high attenuation loss of signals in water, very low frequencies (i.e., long wavelengths) must be used. A wavelength of 50 m was selected for the proposed navigation system, which results in the operating frequency of 1 kHz in seawater and 80 kHz in fresh water. For such a wavelength, even antennas with a length of a quarter of operating wavelength are too bulky for a portable system. Therefore, electrically small antennas must be used [21]. Such antennas are much shorter than the wavelength of the operating signal and are generally less efficient. They are also more challenging to design than larger full-sized antennas. Impedance matching is needed for electrically short antennas to be able to operate at 50 ohms and it can be achieved with the help of lump components. Many types of antennas for underwater applications are discussed in the literatures [22][23][24][25][26][27]. Most of them are wideband antennas for underwater communications in the MHz region; however, the basic principles of design and performance of these antennas apply also for lower frequencies in the kHz region [28]. These antennas are based mostly on loops, long wires, and dipoles. They offer bi-directional radiation pattern similar to the number "8". However, the radiation characteristics of these antennas significantly change due to conductivity of the medium, for example from freshwater to seawater and in space. Therefore, to allow for accurate range estimate, 3D radiation patterns of the chosen antenna must be simulated or measured for all relevant water environments and their influence considered in the range computations. As stated in [29], antennas which are based on magnetic field outperform electric antennas in mediums with conductivity. Therefore, loop antennas seem to be particularly attractive for underwater applications. In [9], a low loss electrically small loop antenna sealed in a waterproof box was proposed for underwater applications. Loop antennas can be built from a long bendable copper wire and shaped to a N-turn rectangular or circular loop. This shaping can be done already shortly before use allowing for efficient stowage. As the loop antennas can be made out of standard copper wire, a low-cost solution seems to be possible. More accurate estimations of the required volume and costs require extensive full-wave electromagnetic simulations and measurements and are beyond the scope of this theoretical study.
The frequency range from 8.3 kHz to 275 GHz is allocated and therefore transmission only in certain bands with specified requirements is allowed [30]. There are no such restrictions for the proposed operation frequency of 1 kHz in seawater. The second chosen frequency of 80 kHz for the operation in fresh water is allocated to maritime radio navigation and radiolocation. The amount of allowed power is not specified for this band but it is only noted, that the transmission is allowed only on condition that it does not cause harmful interference to other services to which these bands are allocated. Our system should transmit signals only underwater. However, due to radiation pattern characteristics of the used antennas, underwater propagation direction, and partially transparent water/air boundary, a small portion of the navigation signal transmitted underwater would be radiated in air. If this radiated energy is above the locally allowed power level, transmit power can be reduced without accuracy loss by extending the observation time. The proposed underwater navigation system can of course operate at frequencies lower than 80 kHz in fresh water and when some transmission restrictions cannot be met, not allocated frequencies can be used.

Underwater Positioning
In this section, the impact of the expected GNSS position accuracy of each buoy on the underwater navigation is discussed. We also derive the horizontal underwater position error and analyze possible circular setups of buoys in terms of expected horizontal position accuracy and system availability area both depending on the diving depth of the user.

Determination of Underwater Position Accuracy
Each buoy is equipped with a multi-constellation, dual frequency GNSS receiver, and hence, able to compute its own position. The GNSS positioning is an iterative process described detailed and repeatedly in the literature, for instance in [19]. After the self-localization, the buoy sends its position information underwater. The underwater vehicle or diver receives this transmitted signal and can extract two key points from it: (a) the estimated range to the buoy and (b) the buoys transmitted position. To compute a full 3D position independently, at least three buoys have to be within the reception range of the user. Please note that in comparison to the GNSS concept, we do not need to resolve any clock offset since the range to each buoy is estimated in a one-way manner by examine the signal amplitude attenuation at the receiver.
To solve the resulting non-linear underwater user position equations, an iteratively Weighted Least Square (WLS) algorithm is applied: where H B = [u 1 , u 2 , . . . , u N ] T is the geometry matrix containing in each row the transpose of all 3D Line of Sight (LoS) vectors u i pointing from the ith buoy to the underwater user. The underwater vector's range error is denoted by ∆r B and is composed of two main sources: one resulting from the uncertainty of the ith buoy's own position ∆r i x and one caused by the underwater propagation of signal from the ith ∆r i p . Both error sources are discussed below in detail.

Impact of GNSS Position Error
The first component, the GNSS buoy position error, which maps onto the underwater range error, can be expressed as: Thus, ∆r i x one dimensional range error that results after mapping 3D GNSS position error of the kth buoy ∆x i,GNSS (1 : 3) = ∆x i onto the underwater range. In case N buoys are within the maximum range of reception, the range error vector due to the individual GNSS position errors can be expressed as: To protect the GNSS performance, two assumption are made: First, the selection of the GNSS antenna and mounting of the antenna are done such that GNSS signal reflection from the water surface are suppressed and they do not affect GNSS ranging. Second, we assume the usage of GNSS fault detection and exclusion methods. Therefore, faulty GNSS signal can be detected and are not used for the buoy's localization. Therefore, we can expect that the remaining GNSS errors can be modeled as a Gaussian distributed random variable with zero mean and variance matrix Σ r x . It is assumed that all buoys are equipped with the same GNSS receiver and antenna, and they see the same set of satellites with same satellite geometry (i.e., they are placed relatively close). In addition, we can say that all buoys experience the same GNSS position uncertainty, and so, their GNSS position covariance matrix can be assumed to be identical: Σ x = Σ x i , ∀i = 1, . . . , N as a consequence, we can formulate the range covariance matrix based on the GNSS position error as: Please note that each buoy transmits its position and associated position error covariance matrix within the underwater signal range. However, the GNSS clock error will not affect the underwater range error.

Impact of Underwater Signal Propagation Range Error
The second error component is caused by the underwater signal propagation. The error related to the underwater signal propagation ∆r k p depends mainly on the travelled distance of the signal.
For simplicity, we will assume a perfectly homogeneous medium, although the permittivity may vary slightly with temperature and pressure, and although water currents may change salinity locally. Thus, the underwater ranging error can be assumed to follow a Gaussian distribution with zero mean and variance σ 2 r i . Recall Equation (8), we can express this term as where r i is the traveled distance underwater, K r i is a constant factor and δ is skin depth in meter. Additional information on the ranging error can be found in Appendix A. We also assume that the underwater signal transmitter is the same for all buoys, i.e., K r i = K r , ∀i = 1, . . . , N. The covariance matrix of the underwater propagation error is given by: Figure 7 illustrates the standard deviation of this error over the traveled distance. As it can be seen, for a distance larger than 30 m the standard deviation increases exponentially. For ranges above 80 m, the deviation is bigger than 1 m. This means that the position performance decreases significantly with ranges larger than 80 m. To limit the propagation range error, we have selected a maximum range of 100 m for our investigations. Therefore, we bound the standard deviation of the propagation range error to a maximum of 20 m. where is the traveled distance underwater, is a constant factor and δ is skin depth in meter. Additional information on the ranging error can be found in Appendix A. We also assume that the underwater signal transmitter is the same for all buoys, i.e., = , ∀ = 1, … , . The covariance matrix of the underwater propagation error is given by: Figure 7 illustrates the standard deviation of this error over the traveled distance. As it can be seen, for a distance larger than 30 m the standard deviation increases exponentially. For ranges above 80 m, the deviation is bigger than 1 m. This means that the position performance decreases significantly with ranges larger than 80 m. To limit the propagation range error, we have selected a maximum range of 100 m for our investigations. Therefore, we bound the standard deviation of the propagation range error to a maximum of 20 m. Please note that a precise clock might affect this type of ranging error as it will lead to a more precise and stable complex phasor estimation.
Putting everything together in Equation (12), we can reformulate the underwater user position error equation as: The overall underwater position error covariance matrix is then obtained by: The underwater weighting matrix is given by: Please note that a precise clock might affect this type of ranging error as it will lead to a more precise and stable complex phasor estimation.
Putting everything together in Equation (12), we can reformulate the underwater user position error equation as: Remote Sens. 2020, 12, 3636 12 of 23 The overall underwater position error covariance matrix is then obtained by: buoy ∆ ,GNSS 1: 3 = ∆ onto the underwater range. In case N buoys are within the maximum range of reception, the range error vector due to the individual GNSS position errors can be expressed as:

Wrong symbol in Equation 20
(page 12 in PDF) The overall underwater position error covariance matrix is then obtained by: ℇ buoy ∆ ,GNSS 1: 3 = ∆ onto the underwater range. In case N buoys are within the maximum ra reception, the range error vector due to the individual GNSS position errors can be expressed as:

Wrong symbol in Equation 20
(page 12 in PDF) The overall underwater position error covariance matrix is then obtained by: The underwater weighting matrix is given by: To evaluate the overall underwater position error accuracy, we can state the underwater position covariance matrix by: Hence, the underwater position accuracy depends mainly on the distribution of the buoys. As there will be a strong attenuation of the signal power in water, this range helps to ascertain the signal and ranging quality.
Then for evaluating the underwater position, we select the expected horizontal position standard deviation as metric of interest and it is given by: Please note that the underwater depth of the receiver is also estimated and not constraint by any additional sensor or method. If this would be the case, the horizontal position accuracy would be even enhanced. Due to the reduction of estimation parameters, the measurement redundancy increases and the horizontal standard deviation reduces.
In order to assure a certain quality of the underwater position performance, we used the Horizontal underwater Dilution of Precision (HuDOP), which is given by This value provides an idea of how strong the propagation of buoy's range errors onto the underwater position solution is. The higher this value is, the higher the position error will be.

Results
In this section, the expected performance of the above described system is evaluated in terms of horizontal position uncertainty and expected service volume. The values for the underwater signal and water quality are assumed to be the following: a constant magnitude of K r = 5 × 10 −17 and a skin depth of δ = 7.96 m. The full set of fixed underwater signal parameters is given by the Table 2.
Although we chose here a fresh water scenario, all provided results in this paper are also valid for a salt water scenario of an underwater signal carrier frequency of 1.0 kHz. Evidently, the location or placement of the buoys will have a major impact on the underwater performance. We investigated circular distribution of N buoys, where N-1 buoys distributed evenly on a circle of radius r buoy and one is placed in the center of the circle. An example of nine buoys is represented in Figure 8.
The center of operation is located on the Lake Wessling in Germany and the georeferenced location is latitude 48.072898 • N, longitude 11.250898 • E, and at an altitude of 596 m. It is also illustrated in Figure 9. The maximum range of each buoy is assumed as 100 m since otherwise the underwater signals will be too attenuated to reach the underwater user. All investigations are based on GPS and Galileo constellation seen on February 14th, 2020 at 10 am UTC. This day was chosen arbitrary and serves only as an example. Appendix C will give more details about the GNSS reception. In this section, the expected performance of the above described system is evaluated in terms of horizontal position uncertainty and expected service volume. The values for the underwater signal and water quality are assumed to be the following: a constant magnitude of = 5 × 10 and a skin depth of = 7.96 m. The full set of fixed underwater signal parameters is given by the Table 2. Although we chose here a fresh water scenario, all provided results in this paper are also valid for a salt water scenario of an underwater signal carrier frequency of 1.0 kHz. Evidently, the location or placement of the buoys will have a major impact on the underwater performance. We investigated circular distribution of N buoys, where N-1 buoys distributed evenly on a circle of radius buoy and one is placed in the center of the circle. An example of nine buoys is represented in Figure 8. The center of operation is located on the Lake Wessling in Germany and the georeferenced location is latitude 48.072898° N, longitude 11.250898° E, and at an altitude of 596 m. It is also illustrated in Figure  9. The maximum range of each buoy is assumed as 100 m since otherwise the underwater signals will The center of operation is located on the Lake Wessling in Germany and the georeferenced location is latitude 48.072898° N, longitude 11.250898° E, and at an altitude of 596 m. It is also illustrated in Figure  9. The maximum range of each buoy is assumed as 100 m since otherwise the underwater signals will be too attenuated to reach the underwater user. All investigations are based on GPS and Galileo constellation seen on February 14th, 2020 at 10 am UTC. This day was chosen arbitrary and serves only as an example. Appendix C will give more details about the GNSS reception. We consider the underwater positioning system as available only if the signals of at least three buoys are received at the same time by the underwater receiver. The underwater position accuracy is simulated in terms of expected horizontal position error standard deviation as described in the previous section. Different scenarios with different number of buoys and radius of buoys were analyzed. The parameters of different scenarios are presented in Table 3.

Scenario Parameter
Value Range Figure 9. Illustration of the center of operation located in the middle of Lake Wessling, Germany.
We consider the underwater positioning system as available only if the signals of at least three buoys are received at the same time by the underwater receiver. The underwater position accuracy is simulated in terms of expected horizontal position error standard deviation as described in the previous section. Different scenarios with different number of buoys and radius of buoys were analyzed. The parameters of different scenarios are presented in Table 3. For the following results, we only consider underwater position solutions with a HuDOP smaller or equal to 10. The expected horizontal position performance for a water depth of 30 m and buoy distribution radius of 30 m for the different numbers of buoys is illustrated in Figure 10. The locations of each buoy are indicated by a black circle. If less than three buoys are received by the user, the system is not available and the corresponding area stays white. The expected horizontal underwater position standard deviation can go up to 95 m. However, this normally only occurs on the border of the availability area. Additional results for different sets of parameters and a water depth of 30 m can be found in Appendix D. Please note that the horizontal position accuracy is limited by the accuracy of the GNSS receivers. That is, even with no underwater ranging errors, any transmitted position error of the buoys, which results from their GNSS position estimation, would propagate to an underwater position error. It is observed from the figures that the horizontal underwater position accuracy in most of the available area is below 5 m and it drops down very quickly to nearly above 70 m at the edges. The upper value of 95 m is mainly driven by the maximum allowable HuDOP value of 10. In the white areas, the system is unavailable for the underwater positioning system. It is observed from the figures that the horizontal underwater position accuracy in most of the available area is below 5 m and it drops down very quickly to nearly above 70 m at the edges. The upper value of 95 m is mainly driven by the maximum allowable HuDOP value of 10. In the white areas, the system is unavailable for the underwater positioning system.
We have analyzed the performance for water depths from 1 to 100 m. The resulting averaged horizontal position accuracy for the availability area is shown Figure 11. (c) (d) It is observed from the figures that the horizontal underwater position accuracy in most of the available area is below 5 m and it drops down very quickly to nearly above 70 m at the edges. The upper value of 95 m is mainly driven by the maximum allowable HuDOP value of 10. In the white areas, the system is unavailable for the underwater positioning system.
We have analyzed the performance for water depths from 1 to 100 m. The resulting averaged horizontal position accuracy for the availability area is shown Figure 11. The different lines represent the results for different buoy circle radius indicated by the legend. It is observed that for all setups and water depths smaller than 40 m the performance is very similar and shows a similar trend. That is up to the depth of 50 m the average standard deviation is below 10 m. Below this depth the error increase will be exponential, and this highly correlates the underwater propagation error.
As one of our design goals was to achieve a maximum horizontal position error standard deviation less or equal to 3 m, we investigated in particular for this performance. In Figure 12, we see the availability area versus the water depth for all investigated setups. As it is seen, we can only provide a service with the required underwater position performance of less or equal 3 m standard deviation for water depths less than 80 m. The different lines represent the results for different buoy circle radius indicated by the legend. It is observed that for all setups and water depths smaller than 40 m the performance is very similar and shows a similar trend. That is up to the depth of 50 m the average standard deviation is below 10 m. Below this depth the error increase will be exponential, and this highly correlates the underwater propagation error.
As one of our design goals was to achieve a maximum horizontal position error standard deviation less or equal to 3 m, we investigated in particular for this performance. In Figure 12, we see the availability area versus the water depth for all investigated setups. As it is seen, we can only provide a service with the required underwater position performance of less or equal 3 m standard deviation for water depths less than 80 m. underwater propagation error.
As one of our design goals was to achieve a maximum horizontal position error standard deviation less or equal to 3 m, we investigated in particular for this performance. In Figure 12, we see the availability area versus the water depth for all investigated setups. As it is seen, we can only provide a service with the required underwater position performance of less or equal 3 m standard deviation for water depths less than 80 m.

Discussion
Our proposed system is very flexible and has low installation and maintenance costs. It is completely self-calibrating and independent. Moreover, it is world-wide deployable and operational. Based on the above provided simulation results, we are able to obtain optimal radius of the buoy circle depending on the number of available buoys and horizontal underwater position accuracy.
For example, let us have a closer look at the case of five buoys, where four are placed on a circle with radius buoy and one at the center of the circle: From Figure 12a, we can extract the expected availability area depending on the water depth and buoy circle radius for maximum accepted horizontal position accuracy requirement in terms of standard deviation below 3 m. It is observed that the largest availability area of = 16365 m² is available at 30 m water depth with a radius buoy = 40 m. As this is difficult to picture, we introduce here the term equivalent service radius service = ( ). This is the radius leading to the same availability area. In our case, it would be service = 72.2 m. Table 4 summarize the optimal buoy radii to obtain the maximum service area in a water depth of 30 m for different number of buoys. Additionally, the equivalent service radii are also provided. As it is stated, we can achieve a horizontal service area up to 163 m² and 341 m² for 5 buoys with 40 m radius and 11 buoys with 80 m radius are exploited, respectively. The corresponding equivalent service radii would be 72.2 m and 104.2 m.

Discussion
Our proposed system is very flexible and has low installation and maintenance costs. It is completely self-calibrating and independent. Moreover, it is world-wide deployable and operational. Based on the above provided simulation results, we are able to obtain optimal radius of the buoy circle depending on the number of available buoys and horizontal underwater position accuracy.
For example, let us have a closer look at the case of five buoys, where four are placed on a circle with radius r buoy and one at the center of the circle: From Figure 12a, we can extract the expected availability area depending on the water depth and buoy circle radius for maximum accepted horizontal position accuracy requirement in terms of standard deviation below 3 m. It is observed that the largest availability area of A s = 16365 m 2 is available at 30 m water depth with a radius r buoy = 40 m. As this is difficult to picture, we introduce here the term equivalent service radius r service = A s π . This is the radius leading to the same availability area. In our case, it would be r service = 72.2 m. Table 4 summarize the optimal buoy radii to obtain the maximum service area in a water depth of 30 m for different number of buoys. Additionally, the equivalent service radii are also provided. As it is stated, we can achieve a horizontal service area up to 163 m 2 and 341 m 2 for 5 buoys with 40 m radius and 11 buoys with 80 m radius are exploited, respectively. The corresponding equivalent service radii would be 72.2 m and 104.2 m. We also investigated the performances for different water depths. Figure 13 summarizes the equivalent service radius results for different water depths and number of buoys used. The corresponding optimal buoy circle radius is indicated by the color scale. As it can be seen depending on the number of available buoys and water depth of interest, changing the circle radius would result in the highest availability area of the system. We can extract that the deeper the operational area of interest, the smaller the optimal buoy circle radius. The results in this paper are based on theoretical studies only. They provide a first assessment on the potential of the concept. In order to further test and validate this concept, field tests would be mandatory. They would also provide a more resilient assessment of the expected system performance.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 23 studies only. They provide a first assessment on the potential of the concept. In order to further test and validate this concept, field tests would be mandatory. They would also provide a more resilient assessment of the expected system performance. Figure 13. Equivalent service radius versus water depth for different numbers of available buoys; the color scale indicates the optimal radius of the buoy's outer circle.

Conclusions
There are many developed concepts and ready-made underwater navigation systems available. Each of them has some advantages and disadvantages based upon its construction and functioning. Our proposed solution operates on the basis of GNSS systems and underwater ranging with a cleverly designed electromagnetic signal. Our investigations show that theoretically it is possible to develop an underwater navigation system, which is flexible, low cost, and world-wide deployable and operational.
The theoretical background of our system is presented in the paper. The scientific considerations and simulations are explained in detail and characterized. A simplistic LF/ULF navigation signal is designed, which allows estimation of the buoy-terminal range with reasonable accuracy in the highly attenuating transmission medium fresh/sea water. The simulations are performed for different buoys distributions and setups. Our analyses lead to the conclusion that horizontal position accuracy in terms of standard deviation below 3 m for water depths up to 80 m is achievable. For depths higher than 80 m, we can still provide a solution, however with significant loss in accuracy. Furthermore, optimal buoy distribution for the given system requirements and the corresponding service area underwater are determined. For instance, the optimal radius of a circular distribution of 7 buoys (6 evenly distributed on the circle with this radius and one in the center) would be 80 m leading to a service area of 200 m at a water depth of 30 m.

Conclusions
There are many developed concepts and ready-made underwater navigation systems available. Each of them has some advantages and disadvantages based upon its construction and functioning. Our proposed solution operates on the basis of GNSS systems and underwater ranging with a cleverly designed electromagnetic signal. Our investigations show that theoretically it is possible to develop an underwater navigation system, which is flexible, low cost, and world-wide deployable and operational.
The theoretical background of our system is presented in the paper. The scientific considerations and simulations are explained in detail and characterized. A simplistic LF/ULF navigation signal is designed, which allows estimation of the buoy-terminal range with reasonable accuracy in the highly attenuating transmission medium fresh/sea water. The simulations are performed for different buoys distributions and setups. Our analyses lead to the conclusion that horizontal position accuracy in terms of standard deviation below 3 m for water depths up to 80 m is achievable. For depths higher than 80 m, we can still provide a solution, however with significant loss in accuracy. Furthermore, optimal buoy distribution for the given system requirements and the corresponding service area underwater are determined. For instance, the optimal radius of a circular distribution of 7 buoys (6 evenly distributed on the circle with this radius and one in the center) would be 80 m leading to a service area of 200 m at a water depth of 30 m.

Conflicts of Interest:
The authors declare there is no conflict of interest.

Appendix A. Underwater Propagation of Electromagnetic (EM) Waves
The electromagnetic properties of a transmission medium can be characterized by three material constants: • Permittivity ε = ε 0 × ε r (measured in Farad per meter); • Permeability µ = µ 0 × µ r (measured in Henry per meter); • Conductivity σ (measured in Siemens per meter).
Here, the quantities ε 0 and µ 0 are natural constants, while the dimensionless quantities εr and µr give the relative permittivity/permeability compared with free space (in vacuum, ε r = µ r = 1). In water, we generally have µ r ≈ 1, while ε r ≈ 80 for 20 • C. Conductivity is known to depend heavily on the water's salinity: in sea water, we have σ ≈ 4 S/m, while in fresh water σ ≈ 0.05 S/m and in distilled water σ ≈ 0. For simplicity, we will assume a perfectly homogeneous medium, although the permittivity may vary with temperature and pressure, and although water currents may change salinity locally. In [31], conductivity is reported to increase by less than 3% per 1K temperature increase in fresh water, and less in salt water. This model mismatch leads to a slight bias which can be calculated using Equations (2), (3), and (7). For instance, with a vertical temperature gradient 0.5 K/meter in fresh water, the homogeneous model leads to a ranging bias of 0.4 m in fresh water at 10 m depth. However, note that this effect could easily be calibrated, if the receiver is equipped with a temperature sensor.
Consider now an EM wave radiated from a point source located at r = 0, which is oscillating with angular frequency ω = 2πf. The electric field generated by this source at distance r and at time t has the following form (∝ denotes proportionality): Thus, propagation is characterized by a spatial and temporal oscillation, by an exponential attenuation term and a radial attenuation term. The dependency on the specific transmission medium, taking into account the above material constants, can be elegantly and compactly described by a single quantity, the complex-valued wave number: where j = √ −1 is the imaginary unit and k is measured in 1/meter. Interestingly, the real part of this number denoted by (k) gives the number of wavelengths per unit distance, while the imaginary part (k) gives the attenuation per unit distance. From the above equation, we can see that in fresh water or sea water, we have σ ω for all technical frequencies up to the MHz range, hence clearly for all frequencies that are relevant for this work (kHz range). Therefore, we may approximate k as follows: where we introduced the skin depth δ = 2 µωσ . The skin depth is measured in meters and denotes the distance after which the medium has attenuated the wave amplitude by 1 Neper (or the wave power by 8.7 dB). In addition, the wavelength defined as λ = 2π Re(k) is just 2π skin depths. Skin depth, attenuation, and wavelength as functions of frequency for fresh water (σ = 0.05 S/m) and sea water (σ = 4 S/m) are illustrated in Figure A1.  = ( ) (A5) HDOP = , + , (A6) In our case, HDOP is 0.64, which is a very good value and indicates a very good satellite reception scenario.

Appendix D: Accuracy Performance at Water Depth of 30 m
A complete set of investigated buoy's placement and availability areas including their horizontal position accuracy information as color scheme is shown in Tables A1 and A3. The radius between the center buoy and the outer ones varies from 20 m to 80 m and the total number of buoys from 5 to 11. In the white areas, no underwater positioning was possible. The colors from dark blue to yellow indicate the achievable underwater position accuracy in terms of horizontal position error standard deviation given in meters. These results are obtained at a water depth of 30 m.

Appendix D. Accuracy Performance at Water Depth of 30 m
A complete set of investigated buoy's placement and availability areas including their horizontal position accuracy information as color scheme is shown in Tables A1 and A2. The radius between the center buoy and the outer ones varies from 20 m to 80 m and the total number of buoys from 5 to 11. In the white areas, no underwater positioning was possible. The colors from dark blue to yellow indicate the achievable underwater position accuracy in terms of horizontal position error standard deviation given in meters. These results are obtained at a water depth of 30 m. = ( ) (A5) HDOP = , + , (A6) In our case, HDOP is 0.64, which is a very good value and indicates a very good satellite reception scenario.

Appendix D: Accuracy Performance at Water Depth of 30 m
A complete set of investigated buoy's placement and availability areas including their horizontal position accuracy information as color scheme is shown in Tables A1 and A3. The radius between the center buoy and the outer ones varies from 20 m to 80 m and the total number of buoys from 5 to 11. In the white areas, no underwater positioning was possible. The colors from dark blue to yellow indicate the achievable underwater position accuracy in terms of horizontal position error standard deviation given in meters. These results are obtained at a water depth of 30 m. = ( ) (A5) HDOP = , + , (A6) In our case, HDOP is 0.64, which is a very good value and indicates a very good satellite reception scenario.

Appendix D: Accuracy Performance at Water Depth of 30 m
A complete set of investigated buoy's placement and availability areas including their horizontal position accuracy information as color scheme is shown in Tables A1 and A3. The radius between the center buoy and the outer ones varies from 20 m to 80 m and the total number of buoys from 5 to 11. In the white areas, no underwater positioning was possible. The colors from dark blue to yellow indicate the achievable underwater position accuracy in terms of horizontal position error standard deviation given in meters. These results are obtained at a water depth of 30 m. = ( ) (A5) HDOP = , + , (A6) In our case, HDOP is 0.64, which is a very good value and indicates a very good satellite reception scenario.

Appendix D: Accuracy Performance at Water Depth of 30 m
A complete set of investigated buoy's placement and availability areas including their horizontal position accuracy information as color scheme is shown in Tables A1 and A3. The radius between the center buoy and the outer ones varies from 20 m to 80 m and the total number of buoys from 5 to 11. In the white areas, no underwater positioning was possible. The colors from dark blue to yellow indicate the achievable underwater position accuracy in terms of horizontal position error standard deviation given in meters. These results are obtained at a water depth of 30 m. = ( ) (A5) HDOP = , + , (A6) In our case, HDOP is 0.64, which is a very good value and indicates a very good satellite reception scenario.

Appendix D: Accuracy Performance at Water Depth of 30 m
A complete set of investigated buoy's placement and availability areas including their horizontal position accuracy information as color scheme is shown in Tables A1 and A3. The radius between the center buoy and the outer ones varies from 20 m to 80 m and the total number of buoys from 5 to 11. In the white areas, no underwater positioning was possible. The colors from dark blue to yellow indicate the achievable underwater position accuracy in terms of horizontal position error standard deviation given in meters. These results are obtained at a water depth of 30 m.