Next Article in Journal
Experimental CanSat Platform for Functional Verification of Burn Wire Triggering-Based Holding and Release Mechanisms
Next Article in Special Issue
VHF Omnidirectional Range (VOR) Experimental Positioning for Stratospheric Vehicles
Previous Article in Journal
Challenges of Ablatively Cooled Hybrid Rockets for Satellites or Upper Stages
Previous Article in Special Issue
Tuning of NASA Standard Breakup Model for Fragmentation Events Modelling
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Comparative Analysis of Multi-Epoch Double-Differenced Pseudorange Observation and Other Dual-Satellite Lunar Global Navigation Systems

1
Department of Engineering Technology, College of Technology, University of Houston, 306 Technology 2 Bldg., Calhoun Road, Houston, TX 77004, USA
2
Department of Electrical & Computer Engineering, College of Technology, University of Houston, 306 Technology 2 Bldg., Calhoun Road, Houston, TX 77004, USA
3
Department of Astronautics and Aeronautics, Chubu University, 1200 Matsumoto-cho, Kasugai 487-8501, Japan
4
Department of Aeronautics and Astronautics, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
*
Author to whom correspondence should be addressed.
Submission received: 21 May 2021 / Revised: 8 July 2021 / Accepted: 14 July 2021 / Published: 15 July 2021
(This article belongs to the Special Issue New Space: Advances in Space Science and Engineering)

Abstract

:
In this study, dual-satellite lunar global navigation systems that consist of a constellation of two navigation satellites providing geo-spatial positioning on the lunar surface were compared. In our previous work, we proposed a new dual-satellite relative-positioning navigation method called multi-epoch double-differenced pseudorange observation (MDPO). While the mathematical model of the MDPO and its behavior under specific conditions were studied, we did not compare its performance with other dual-satellite relative-positioning navigation systems. In this paper, we performed a comparative analysis between the MDPO and other two dual-satellite navigation methods. Based on the difference in their mathematical models, as well as numerical simulation results, we developed useful insights on the system design of dual-satellite lunar global navigation systems.

1. Introduction

In recent years, communication and navigation architecture for lunar exploration programs has been of great interest [1]. In particular, the estimation of a rover vehicle’s position on the lunar surface is one of the key technologies for the successful operation of the rover, mapping resources, and making scientific observations on the lunar surface. It is well-known that cold-trapped volatiles, including water-ice, in lunar Permanently Shadowed Regions (PSRs) could be a high priority resource for future space exploration. As PSRs never receive direct sunlight, visual-odometry based navigation methods, such as simultaneous localization and mapping (SLAM), will be considerably constrained. Therefore, some alternative is needed to realize long and efficient exploration of the PSRs. Additionally, we aim to provide navigation information to multiple users on the lunar surface as the locations of various resources are not known precisely, and wide-range exploration by multiple small rovers is considered as promising approach [2]. With these two trends in mind, multiple-user navigation system that can be used in PSRs is of immediate demand.
As one feasible approach to establish the multiple-user navigation system for the lunar surface applications, several groups have studied the use of weak signals, i.e., the spill-over of the beams irradiated from global navigation satellite systems (GNSS) that serve Earth surface and proximity: the weak signals technology was investigated by [3] for the first time, and applications to the lunar navigation were extensively studied by a great deal of research [4,5,6,7,8,9,10]. While they are capable of providing navigation signals at the middle latitude of the lunar surface, they are not available at the far side and polar regions of the Moon due to invisibility. As an alternative method, constellations of global navigation satellites around the Moon have been studied [11,12]. Additionally, a combination of the weak signals and a relatively small constellation of global navigation satellites around the Moon has been studied [13]. While they can provide geo-spatial positioning to the entire Moon and its proximity, the transportation cost to inject satellites into multiple lunar orbits, as well as the ground station cost to operate a large number of lunar satellites, is not affordable at the early stage of the lunar exploration programs.
In an attempt to reduce the cost of global navigation satellite systems, dual-satellite lunar global navigation that consists of a constellation of two navigation satellites have been studied. Originally, dual-satellite global navigation methods have been studied for Earth GNSS in the literature [14,15] and applying them to lunar GNSS domain [16,17]. We also proposed a new dual-satellite relative-positioning navigation method called multi-epoch double-differenced pseudorange observation (MDPO) [18]. While the mathematical model of the MDPO and its behavior under specific conditions were studied, we did not sufficiently compare its performance with other dual-satellite global navigation systems. In this paper, we show a comparative analysis between the MDPO and other selected navigation methods. More specifically, we study and compare the following three navigation methods, (1) MDPO, (2) joint time difference of arrival and frequency difference of arrival (TDOA–FDOA) [15], and (3) two-way ranging [19], and discuss the pros and cons of each method.
These three dual-satellite navigation methods use different types of observations, namely passive ranging, passing ranging, and Doppler, or active ranging (two-way ranging), as shown in Table 1. As most dual-satellite navigation methods can be classified into one of these types of observations, a comparative analysis and evaluation of these methods will provide a benchmark of dual-satellite relative-positioning lunar GNSS. For instance, recent work in joint Doppler and ranging (JDR) [16], which converts a differenced Doppler shift into a pseudo-pseudorange using the Law of Cosines and integrates it with pseudorange observation, can be classified as evolved families of joint TDOA–FDOA: otherwise as evolved families of two-way ranging, if it employs two-way ranging observation instead of pseudorange observation.
Apart from these three types, dual-satellite navigation can be established only with Doppler observation. For example, in [17], it was successfully shown that Doppler Based Autonomous Navigation (DBAN) can operate with as few as one lunar orbiter and a reference station and enable autonomous positioning of crewed missions. However, we will rule out Doppler-only navigation from this comparative analysis, as it requires a long observation period and is not able to provide position estimates as quickly as the other methods do.
In this research, the target of our study comprises a micro-sized satellite and rover systems. In that case, power generation capability is limited by size and, consequently, not compatible with a high-standard clock source, such as the deep space atomic clock (DSAC) [20]. When the clock bias of space/user segment is not ignorable, or when satellite orbit determination error is not ignorable, the existing joint TDOA–FDOA method must be updated to a double-differenced form [21], the so-called double-differenced Time Of Arrival (TOA)–Frequency Of Arrival (FOA), to cope with the clock bias and orbit determination errors. Likewise, the two-way ranging method must be updated to a single-differenced form, the so-called single-differenced two-way ranging, to cope with the orbit determination errors. These updates have been taken into account to give a fair comparison.
The selected three navigation methods have different characteristics in terms of navigation accuracy and system complexity. The comparative study of these three navigation methods is shown in Table 1, which could assist designers to choose an appropriate method for their own purposes. For example, MDPO can provide navigation information to multiple users at a time through passive ranging but requires observations from multiple epochs. Double-differenced TOA–FOA can provide navigation information to multiple users at a time with observation from a single epoch but requires Doppler observation and pseudorange observation. Single-differenced two-way ranging can provide relatively high-accuracy navigation information with observation from a single epoch but only to a single user per one set of radio signals at a time, and also requires an active ranging, i.e., radio signal power emission at the user segment. We also compared the navigation accuracy of these three methods by numerical simulations under the selected conditions.
This paper consists of the following sections. In Section 2, we discuss the assumptions for our study. In Section 3, the mathematical models of the three navigation methods are presented. In Section 4, the achievable user position accuracies of three navigation methods are analyzed by numerical simulation and comparative studies are discussed. In Section 5, we summarize the key insights by analyzing the simulation results, and provide suggestions from a system design point of view.

2. Assumptions

In recent years, NASA and private sectors have extensively studied and developed unmanned micro mobile robots for lunar surface exploration [22,23]. The target of our study comprises a micro-sized satellite and rover system whose power generation capability is limited by size and, consequently, not compatible with the deep space atomic clock (DSAC). In this case, the best current clock technology that is compatible with the micro-sized satellite is the Chip Scale Atomic Clock (CSAC).
As reported in [24], while CSAC can suppress the frequency instability of the clock down to about 1 part per billion (ppb) for 24 h, CSAC incurs several tens to hundreds of meters of error in pseudorange observations after 24 h, which further increases over time. As a result, using CSAC inevitably requires pseudorange-based navigation systems to conduct frequent estimations of the satellite and/or user clock bias using Earth ground stations, which is very challenging in lunar GNSS due to the limitation of the availability and number of earth ground stations that are capable of Earth–Moon distance communication. In summary, the following assumptions were used:
The bias of the satellite clock is not ignorable due to the limited capacity of micro-satellites;
The bias of the rover clock is not ignorable due to the limited capacity of micro-rovers;
The satellite orbit determination error is not ignorable due to the limitation of the availability and number of Earth ground stations.

3. Mathematical Models of the Three Navigation Methods

In this section, we derive the formulas for the three navigation methods listed in Table 1. These navigation methods commonly use pseudorange observation and/or pseudodoppler observation. Therefore, we first introduce the formulas of pseudorange observation and pseudodoppler observation later to derive the formulas of the three different navigation methods.

3.1. Pseudorange Observation

In the conventional Time Of Arrival (TOA) algorithm, the pseudorange ( ρ ) measurement between one user (user1) and one satellite (satellite1) is presented by the following equation [14];
ρ R S ( t i ) = r R S ( t i ) + c ( d τ R ( t i ) d T S ( t i s ) ) + ω r R S ( t i )
r R S ( t i ) = | X s ( t i s ) X R ( t i ) + d X R s a |
where X s ( t i s ) = ( x S ( t i s ) , y S ( t i s ) , z S ( t i s ) ) is the satellite1 position at the time of signal transmission t i s ; X R ( t i ) = ( x R ( t i ) ,   y R ( t i ) , z R ( t i ) ) is the user1 position at the time of signal reception t i ; c is the speed of light; d τ R is the user clock bias; d T S is the satellite clock bias; d X R s a corresponds to the user1 position transition due to the Sagnac effect; and ω r R S is the range receiver observation error. In this study, we assume that the range receiver observation error ω r R S follows a white Gaussian distribution with the standard deviation of σ ω r .
The coordinate frame of the satellite position and user position is based on a local topocentric frame, i.e., the x-axis points local east, the y-axis points local north, and the z-axis points local up (East-North-Up), hereafter. The equations are formulated using the relative position between the satellite and the user, and both the user positions have a constant rotational offset with respect to the Moon-centered inertial frame. In other words, the user position is changed due to the Moon rotation during the signal traveling time from the satellite to the user, which appears as the Sagnac effect in Equation (2). The Shapiro effect was not considered, as the selected three methods do not aim to have a sub-meter accuracy.

3.2. Pseudodoppler Observation

In the Frequency Of Arrival (FOA) algorithm, the pseudodoppler shift ( Ω R S ) between user1 and satellite1 is given as [14];
Ω R S ( t i ) = f R S ( t i ) + ( d f R ( t i ) d f S ( t i s ) ) + ω d R S ( t i )
f R S ( t i ) = f 0 c { ( V R S ( t i s ) + d V R   s a S ) · ( X s ( t i s ) X R ( t i ) + d X R s a ) T | X s ( t i s ) X R ( t i ) + d X R s a | }
where V R S ( t i s ) = ( V x R S ( t i s ) , V y R S ( t i s ) , V z R S ( t i s ) ) is the satellite velocity relative to the user at a time of t i s ; f 0 is the radio wave frequency; c is the speed of light; d f R is the user frequency bias; d f S is the satellite frequency bias; d V R   s a S corresponds to the satellite relative velocity variation due to the Sagnac effect; and ω d R S is the Doppler receiver observation error. In this study, we assume that the Doppler receiver observation error ω d R S follows a white Gaussian distribution with the standard deviation of σ ω d .

3.3. Multi-Epoch Double-Differenced Pseudorange Observation

Multi-epoch double-differenced pseudorange observation (MDPO) is a multi-user, pseudorange-based navigation algorithm. Using multi-epoch double-differenced observations reduces the number of navigation satellites required from four to two, while dealing with the instability of the satellite clock at the same time, as shown in Figure 1.
In this paper, we only explain important equations which are necessary to clarify the differences between the MDPO and other methods, while the complete formulation derivation of the MDPO can be found in our prior research [18].
MDPO uses double-differenced pseudorange observations to eliminate the clock bias of the space segment and user segment as well as satellite orbit determination error, by subtracting four pseudorange measurements between two users (user1 and user2) and two satellites (satellite1 and satellite2) as shown in Equations (5)–(9):
ρ 1 1 ( t i ) = r 1 1 ( t i ) + c ( d τ 1 ( t i ) d T 1 ( t i 1 ) ) + ω r 1 1 ( t i )
  ρ 1 2 ( t i ) = r 1 2 ( t i ) + c ( d τ 1 ( t i ) d T 2 ( t i 2 ) ) + ω r 1 2 ( t i )  
  ρ 2 1 ( t i ) = r 2 1 ( t i ) + c ( d τ 2 ( t i ) d T 1 ( t i 1 ) ) + ω r 2 1 ( t i )  
  ρ 2 2 ( t i ) = r 2 2 ( t i ) + c ( d τ 2 ( t i ) d T 2 ( t i 2 ) ) + ω r 2 2 ( t i )  
Δ ρ ( t i ) = ρ 1 1 ( t i ) ρ 1 2 ( t i ) ( ρ 2 1 ( t i ) ρ 2 2 ( t i ) ) = r 1 1 ( t i ) r 1 2 ( t i ) ( r 2 1 ( t i ) r 2 2 ( t i ) ) + ω r 1 1 ( t i ) ω r 1 2 ( t i ) ( ω r 2 1 ( t i ) ω r 2 2 ( t i ) ) = Δ r ( t i ) + Δ ω r ( t i )
where Δ ( · ) is the double difference operator. In order to remove the satellite and user clock bias errors effectively from the double-differenced observations, time synchronization between the two receivers is essential. In the common GNSS systems, the time synchronization can be achieved in the position calculation process by estimating the user clock bias at the same time. However, the user clock bias is removed in the double-differenced observation and cannot be estimated. In our proposed method, the time synchronization is assumed to be achieved by the frame synchronization of the navigation message. As the maximum range rate of the pseudorange observation from the low lunar orbiting satellite is about 1.3 km/s, the resulting range error is no larger than 1.3 m if the time synchronization error is maintained under 1 ms. It is acceptable if the targeting navigation accuracy is tens of meters.
In the double difference method, user2 is used as a reference station whose position is fixed and known, and the position of user1 is estimated in relation to the position of user2; i.e., user2′s position is referenced as the origin of navigation (0,0,0). In a lunar navigation system, the lander can be used as a reference station (user2), and its geodetic position is used as the origin of navigation. The geodetic position of the lander must be obtained in advance of the start of the rover navigation by other means, such as identification by satellite image: this is proven technique such as the Lunar Reconnaissance Orbiter Camera (LROC) successfully identified the landing coordinates of China’s Chang’e 5 lander with a reported accuracy of ±20 m [25]. Hereafter, the rover corresponds to user1, and the lander corresponds to user2.
In the MDPO method, multiple double-differenced pseudorange observations, i.e., Δ ρ ( t k ) , , Δ ρ ( t k + N 1 ) , are obtained from multiple epochs, i.e., t i = t k , , t k + N 1 , where N is the number of observation epochs, and k is the epoch number at which the estimation starts. It is important to note that the rover position must be fixed in place during all multi-epoch observations taken, in order to keep the number of estimation parameters lower than the number of observation equations. Otherwise, the rover position cannot be identified deterministically by the MDPO, and the rover position accuracy changes depending on the quality of other navigation information used during multi-epoch observations. Hereafter, X R ( t k ) = ( x R ( t k ) ,   y R ( t k ) , z R ( t k ) ) represents a fixed rover position during multi-epoch observations t k   t k + N 1 .
The rover position can be estimated by solving the following Newton–Raphson equations iteratively. The following equations correspond to ‘2D MDPO’ in [18], which calculates an estimated two-dimensional (X–Y) position by the Newton-Raphson equation and a rover altitude (Z) by using a lunar digital elevation model (DEM) as we discuss later. In 2D MDPO, the number of multi-epoch observations can be reduced to as low as 2 (N = 2). The formulation for three-dimensional position calculation can be found in our previous paper [18]. First, we define a new parameter R for t i = t k , , t k + N 1 :
R ( t i ) = Δ ρ ( t i ) Δ r 0 ( t i ) i = k , , k + N 1
where R is the difference between the measured double-differenced pseudorange value, i.e., Δ ρ , and the calculated double-differenced range on an initial estimated value of the rover position X R 0 ( t k ) , i.e., Δ r 0 . Then, the following equations can be derived:
R = G d X + w
  R = [ R ( t k ) R ( t k + N 1 ) ] T
  d X = [ d x , d y ]
  w = [ Δ ω r ( t k ) Δ ω r ( t k + N 1 ) ] T
  G = [ R x R y ]  
where G in Equation (15) is called an observation matrix, which is equivalent to the Jacobian of R with regard to X .
By solving the least-square problem that minimizes the residual error | R G d X | , an estimated value of d X , defined as d X ^ , is obtained:
d X ^ = ( G T G ) 1 G T R .
Then, a new estimated value X R 1 ( t k ) = ( x R 1 ( t k ) ,   y R 1 ( t k ) ,   z R 1 ( t k ) ) is given by Equation (17), which provides a better fit to the observation.
X R 1 ( t k ) = X R 0 ( t k ) + d X ^ .
This estimation process continues until the number of iterations reaches the designed value n , i.e., X R 1 ,   X R 2 X R n , and then the final estimated value X R n ( t k ) is acquired.
In GNSS terminology, ( G T G ) 1 is known as the dilution of precision (DOP) matrix, which is used to specify error propagation as a mathematical effect of the navigation satellite geometry on positional measurement precision. We define the DOP matrix as
D O P = [ ( σ D O P   11 ) 2 ( σ D O P   1 N ) 2 ( σ D O P   N 1 ) 2 ( σ D O P   N N ) 2 ] = ( G T G ) 1
where σ D O P is the elements of D O P . Using D O P , the achievable rover position error, i.e., d X ^ d X , at a time of t k can be obtained by
U P E ( t k ) = | d X ^ ( t k ) d X ( t k ) | = j = 1 N ( σ D O P   j j ) 2 × σ Δ ω
where σ Δ ω is the standard deviation of double-differenced receiver observation errors and U P E represents user position error, which is the distance between the rover’s true position and an estimated rover position. It is important to highlight that the standard deviation of MDPO’s double-differenced receiver observation errors, i.e., σ Δ ω , is amplified from the standard deviation of the original receiver observation errors, i.e., σ ω r , as a result of the double-differencing process, in particular ω r 1 1 ( t i ) ω r 1 2 ( t i ) ( ω r 2 1 ( t i ) ω r 2 2 ( t i ) ) in Equation (9), and becomes as large as σ Δ ω = σ ω r 2 + σ ω r 2 + σ ω r 2 + σ ω r 2 = 2 σ ω r assuming that the receiver observation errors follow a white Gaussian distribution. Further, by defining G D O P as
G D O P = j = 1 N ( σ D O P   j j ) 2
Equation (19) can be written as
U P E ( t k ) = G D O P × σ Δ ω
As mentioned in the previous section, we assume that the receiver observation errors follow a normal distribution with a zero mean (i.e., Gaussian white noise). As such, U P E also follows a 1D Gaussian distribution, and 95 percent of it lies inside the interval from 2 s to + 2 s , where s is the standard deviation. As a performance index, this research uses 2drms ( 2 s or 95 percent confidence), which is commonly used in two-dimensional position estimation problems:
U P E ( t k ) ( 2 d r m s ) = G D O P × 2 σ Δ ω
Moreover, considering that the U P E value, as well as the G D O P value, changes over time according to the satellite positions relative to the rover, an indicator that represents the overall U P E over the course of the mission time is needed. For this purpose, the T o t a l   U P E is newly defined, along with the T o t a l   G D O P , as below:
T o t a l   U P E ( 2 d r m s ) = 1 m m ( U P E ( t k ) ( 2 d r m s ) ) 2 = T o t a l   G D O P × 2 σ Δ ω
  T o t a l   G D O P = 1 m m ( G D O P ) 2  
where m is the number of MDPO estimations over the course of the mission time. σ Δ ω is independent of time, and can be excluded from the square root without losing generality.
In the 2D MDPO, the rover altitude z R is pre-estimated using a lunar digital elevation model (DEM). As shown in Equation (25), DEM is a function of longitude and latitude, which are not known at the start. The estimation of sequences proceeds in the following sequence: First, X R 0 ( t k ) is estimated using the rover position before its relocation, i.e., X R 0 ( t k ) = X R ( t k 1 ) = ( x R ( t k 1 ) ,   y R ( t k 1 ) ,   z R ( t k 1 ) ) . Then, a new estimated rover position, i.e., X R 1 ( t k ) , is estimated as X R 1 ( t k ) = ( x R 1 ( t k ) , y R 1 ( t k ) , z R ( t k 1 ) ) by Equation (17). z R is not updated at this moment. After that, the altitude of the rover is updated to z R 1 ( t k ) using x R 1 ( t k ) and y R 1 ( t k ) by Equation (25), i.e., X R 1 ( t k ) = ( x R 1 ( t k ) , y R 1 ( t k ) , z R 1 ( t k ) ). The calculation continues until the number of iterations reaches the designed value, i.e., n .
z R i ( t k ) = z R   D E M ( x R i ( t k ) , y R i ( t k ) )
Here, z R   D E M is a lunar DEM that is a function of latitude and longitude. According to Equation (25), as z R changes along with x R and y R , errors in the X–Y position induce errors in the Z position, which ultimately induce errors in estimated x R and y R , and as a result, the T o t a l   U P E deteriorates stochastically. In our research, we did not apply the case in which the rover altitude changes too rapidly, such as the rover dropping off the cliff or roving on steep slopes. In that case, the T o t a l   U P E would not deteriorate too significantly, which was confirmed by numerical simulations in our prior research [18].

3.4. Double-Differenced TOA–FOA

The conventional TDOA–FDOA approach uses single-differenced TOA observation and FOA observation to cope with the satellite clock bias or the user clock bias, but is not designed to cope with the satellite clock bias and the user clock bias at the same time [15]. Therefore, we updated the conventional TDOA–FDOA into a double-differenced form, the so-called double-differenced TOA–FOA. Double-differenced TOA–FOA can determine the user position using double-differenced pseudorange observations and pseudodoppler observations from a single epoch, as shown in Figure 2.
Double-differenced TOA observations were formulated using Equations (5)–(9). Similarly, double-differenced FOA observations can be formulated:
Ω 1 1 ( t i ) = f 1 1 ( t i ) + ( d f 1 ( t i ) d f 1 ( t i 1 ) ) + ω d 1 1 ( t i )
Ω 1 2 ( t i ) = f 1 2 ( t i ) + ( d f 1 ( t i ) d f 2 ( t i 2 ) ) + ω d 1 2 ( t i )
Ω 2 1 ( t i ) = f 2 1 ( t i ) + ( d f 2 ( t i ) d f 1 ( t i 1 ) ) + ω d 2 1 ( t i )
Ω 2 2 ( t i ) = f 2 2 ( t i ) + ( d f 2 ( t i ) d f 2 ( t i 2 ) ) + ω d 2 2 ( t i )
Δ Ω ( t i ) = Ω 1 1 ( t i ) Ω 1 2 ( t i ) ( Ω 2 1 ( t i ) Ω 2 2 ( t i ) ) = f 1 1 ( t i ) f 1 2 ( t i ) ( f 2 1 ( t i ) f 2 2 ( t i ) ) + ω d 1 1 ( t i ) ω d 1 2 ( t i ) ( ω d 2 1 ( t i ) ω d 2 2 ( t i ) ) = Δ f ( t i ) + Δ ω d ( t i )
In order to remove the satellite and user frequency bias errors effectively from the double-differenced observations, it is assumed that time synchronization is achieved by the frame synchronization of the navigation message as we discussed earlier.
Similar to the MDPO, we can use the Newton–Raphson method to solve the equations. First, we define new parameters R 1 and R 2 at t i = t k :
R 1 ( t k ) = Δ ρ ( t k ) Δ r 0 ( t k )
R 2 ( t k ) = Δ Ω ( t k ) Δ f 0 ( t k )
where R 1 is the difference between the measured double-differenced pseudorange value, i.e.,   Δ ρ , and the calculated double-differenced range on an initial estimated value of the rover position X R 0 ( t k ) , i.e., Δ r 0 ; and R 2 is the difference between the measured double-differenced pseudodoppler value, i.e., Δ ρ d , and the calculated double-differenced Doppler on an initial estimated value of the rover position, i.e., Δ f 0 . Then, the following equations can be derived:
R = [ R 1 ( t k ) R 2 ( t k ) ] T
w = [ Δ ω r ( t k ) Δ ω d ( t k ) ] T
The following process is same as the Equations (11), (13) and (15)–(25) with the exception that the double-differenced receiver observation error, i.e., σ Δ ω , is amplified from the standard deviation of the original receiver observation errors, i.e., σ ω r and σ ω d , during the double-differencing process and becomes as large as σ Δ ω = α × ( σ ω r 2 + σ ω r 2 + σ ω r 2 + σ ω r 2 ) + ( σ ω d 2 + σ ω d 2 + σ ω d 2 + σ ω d 2 ) where α changes depending on the geometrical relationship between the satellites and receivers such as V R S . Under the conditions used in this study, α is small to negligible. It is also important to highlight that the above formulation corresponds to a DEM-aided form that calculates an estimated two-dimensional (X–Y) position using a rover altitude that is given by Equation (25).

3.5. Single-Differenced Two-Way Ranging

Two-way ranging method determines Time of Arrival (TOA) of radio signal in round trip and then calculates the distance between the nodes by multiplying the round-trip time by the speed of light.
When the two-way ranging radio signal is initiated at the satellite side and the user is fixed on the lunar surface, as shown in Figure 3, the pseudorange observation is presented as
ρ R S ( t e m i t S ) = r R S ( t e m i t S ) + ω r R ( t r e c e i v e R ) + ω r S ( t e m i t S + t r o u n d t r i p S R + t p r o c e s s i n g S R )
r R S ( t e m i t S ) = | X s ( t e m i t S ) X R + d X R s a | + | X s ( t e m i t S + t r o u n d t r i p S R + t p r o c e s s i n g S R ) X R + d X R s a |
where t e m i t S is the time of signal emission at the satellite; t r e c e i v e R is the time of signal reception at the rover; t r o u n d t r i p S R is a signal round-trip time that is the sum of the onward and return signal traveling time; and t p r o c e s s i n g R is the time to process the signal and rebroadcast at the rover, which are also visually shown in Figure 3. ω r R and ω r S are the range receiver observation error at the user and satellite, respectively. In this study, we assume that the range receiver observation error ω r R and ω r S follow a white Gaussian distribution.
In Equation (36), r R S is written as a function of t e m i t S based on the fact that t r o u n d t r i p can be presented as a function of t e m i t S as long as the satellite orbit dynamics, orbital parameters, the signal processing time t p r o c e s s i n g S R , as well as an estimate of the user position, X R are provided.
Two pseudorange observations between two users (user1, user2) and a satellite can be written as
ρ 1 S ( t i ) = r 1 S ( t i ) + ω r 1 S ( t i )
ρ 2 S ( t i ) = r 2 S ( t i ) + ω r 2 S ( t i )
where t i corresponds to t e m i t S of Equation (35) and ω r R S ( t i ) corresponds to ω r R ( t r e c e i v e R ) + ω r S ( t e m i t S + t r o u n d t r i p + t p r o c e s s i n g S R ) of Equation (35) assuming ω r R and ω r S follow a white Gaussian distribution and are independent of time. The standard deviation of ω r R S ( t i ) is defined as σ ω r .
A method called single difference is used to effectively remove the satellite orbit determination error, i.e., d X s a t   O D S as we discuss in Section 3.6, by subtracting two pseudorange observations between two users (user1, user2) and a satellite:
Δ ρ S ( t i ) = ρ 1 S ( t i ) ρ 2 S ( t i ) = r 1 S ( t i ) r 2 S ( t i ) + ω r 1 S ( t i ) ω r 2 S ( t i ) = Δ r S ( t i ) + Δ ω r S ( t i )
where Δ ( · ) is the single difference operator.
To effectively remove the satellite orbit determination error at the moment of signal emission, i.e., t e m i t S , the clock t i of two pseudoranges, i.e., ρ 1 1 ( t i ) and ρ 2 1 ( t i ) , must be synchronized. This can be easily achieved, for instance, when pseudorange signals are initiated from satellite side at the request of the user. In that case, pseudorange observations are to be obtained at the satellite side, and then transferred to the user by telemetry: the so-called telemetry ranging [26].
Furthermore, the timing of returned-signal reception at the satellite side, i.e., t e m i t S + t r o u n d t r i p S R + t p r o c e s s i n g S R , slightly differs among four pseudorange observations due to the difference in user1 and user2 position, as well as their different signal processing delay times. In order for single difference to effectively remove the satellite orbit determination error at the moment of returned-signal reception, the difference among d X s a t   O D S ( t e m i t 1 + t r o u n d t r i p 1 1 + t p r o c e s s i n g 1 1 ) , d X s a t   O D S ( t e m i t 1 + t r o u n d t r i p 1 2 + t p r o c e s s i n g 1 2 ) , d X s a t   O D S ( t e m i t 2 + t r o u n d t r i p 2 1 + t p r o c e s s i n g 2 1 ) , and d X s a t   O D S ( t e m i t 2 + t r o u n d t r i p 2 2 + t p r o c e s s i n g 2 2 ) must be negligible. This assumption practically holds, unless the distance between user1 and user2 becomes largely apart, or their signal processing delay times are largely different, which holds at least under the simulated conditions of this research.
Single-differenced two-way ranging can determine the user position using single-differenced pseudorange observations from two satellites, i.e., S = 1 ,   2 , at a single epoch as shown in Figure 4. Similar to the other two methods, we can use the Newton–Raphson method to solve the equation. First, we define a new parameter R S at t i = t k :
R S ( t k ) = Δ ρ S ( t k ) Δ r S   0 ( t k ) S = 1 ,   2
where R S is the difference between the measured single-differenced pseudorange value, i.e., Δ ρ S , and the calculated single-differenced range of s-th satellite on an initial estimated value of the rover position X R 0 ( t k ) , i.e., Δ r S   0 . Then, R and w are derived as follows:
R = [ R 1 ( t k ) R 2 ( t k ) ] T  
w = [ Δ ω r 1 ( t k ) Δ ω r 2 ( t k ) ] T  
The following process is the same as the Equations (11), (13) and (15)–(25), with the exception that the single-differenced receiver observation error, i.e., σ Δ ω , is used instead of the double-differenced receiver observation error, i.e., σ Δ ω . The standard deviation of single-differenced receiver observation errors, i.e., σ Δ ω , is amplified from the standard deviation of the original receiver observation errors, i.e., σ ω r , during the single-differencing process and becomes as large as σ Δ ω = σ ω r 2 + σ ω r 2 = 2 σ ω r assuming that the receiver observation errors follow a white Gaussian distribution. It is also important to highlight that the above formulation corresponds to a DEM-aided form that calculates an estimated two-dimensional (X–Y) position using a rover altitude that is given by Equation (25).

3.6. Systematic Errors

In an actual situation, with the presence of other systematic errors as shown in this section, the discussed achievable T o t a l   U P E in Equation (23) will increase. In this section, the theoretical background of systematic errors, as well as their impacts on the T o t a l   U P E , is discussed. As the impact of such errors on the T o t a l   U P E cannot be predicted analytically, we used a numerical simulation, reported in the following section, to quantitatively determine the resulting T o t a l   U P E .

3.6.1. Satellite Orbit Determination Error

In the derived formulas, the pseudorange ρ is calculated on the basis of pre-estimated satellite positions X s = ( x S , y S , z S ) . In an actual situation, satellite orbit determination is not perfect, and pre-estimation of the satellite position entails some error relative to the true positions ( d X s a t   O D S ). According to a general satellite orbit determination process, the error is decomposed along with the satellite velocity direction (Along), satellite zenith direction (Radial), and cross-track direction (Cross). In this simulation, the orbit determination error is defined along with the Along, Radial, and Cross directions and then converted into a user frame:
d X s a t   O D S ( t i ) = T × ( d A l o n g ( t i ) , d R a d i a l ( t i ) , d C r o s s ( t i ) )
where T is a coordinate transformation matrix from the Along, Radial, and Cross directions to a topocentric frame. The definition of the topocentric frame is explained in the previous chapter. In multilateration theory, only satellite orbit determination error in the line-of-sight direction (rover to satellite) matters, and other directions have almost no impact on the rover position error. In the three navigation methods, the line-of-sight direction error is effectively eliminated by either the double-difference or the single-difference equation, along with the satellite, the rover, and the lander clock biases.

3.6.2. Time Tag Error

In the estimation process of the satellite position at a given time, the time tag of the receiver is used to propagate the estimated satellite positions. As we described earlier, the time tag of the receiver clock is initialized by the frame synchronization of the navigation message. In this case, the time tag has some error due to the signal propagation delay between the satellite and user receiver, as well as the processing delay of the navigation message. As a result, an estimated satellite position X s = ( x S , y S , z S ) is deteriorated by the receiver clock bias d τ R ( t i ) , and has some error relative to the true positions ( d X t i m e   t a g S ), such as
d X t i m e   t a g S ( t i ) = ( V x R S ( t i ) , V y R S ( t i ) , V z R S ( t i ) ) × d τ R ( t i )
where ( V x R S , V y R S , V z R S ) is a pre-estimated satellite relative velocity in a topocentric frame. Essentially, the time tag error is mostly eliminated from the estimation by either the double-difference or single-difference equation, except for the ‘difference’ of two user time tags. In the numerical simulation, we only modeled the difference of two user time tags without losing generality.

3.6.3. Signal Processing Delay Time Uncertainty

For the single-differenced two-way ranging, the time to process the signal and rebroadcast at the rover, i.e., t p r o c e s s i n g R , may not be precisely known and have some uncertainty. In this study, uncertainty of the signal processing delay time is modeled as d t p r o c e s s i n g R .

3.6.4. DEM Information Error

As reported in [27,28], current lunar DEM information is developed from remote-sensing data and, as a result, is not perfect. Therefore, the DEM error d z R   D E M , which is the difference between the true rover vertical position z R   t r u e and a pre-given rover vertical position z R   D E M , has a fixed, unknown, non-random bias. The DEM error leads to a position estimation error in the X–Y plane ( x R , y R ) . The impact of the DEM model error on the X–Y position estimation accuracy stochastically changes depending on the satellite position and velocity in relation to the rover and lander position.
d z R   D E M = z R   t r u e z R   D E M

3.6.5. Other Systematic Errors

In the general context of navigation satellite systems, other systematic errors must be considered, such as ionospheric delay, tropospheric delay, antenna phase characteristics, and multipath. However, such errors are negligible, or not detrimental to the rover position estimation in lunar surface navigation systems. Ionospheric delay and tropospheric delay are deemed negligible. Antenna phase characteristics appear in the same way and are almost negligible. We assume that multipath can be suppressed by antenna design as there are few high objects in the surroundings of the rover and lander on the lunar surface. Therefore, these errors can be ignored, and were not considered in this research.

3.7. Design Parameters

3.7.1. DOP and Availability

The spatial position of two satellites is one of the most important design parameters that directly impacts the rover position accuracy. In order to acquire an accurate user position, a small DOP value is required.
In this analysis, we assume the satellite formation that has two satellites placed in the same orbital planes with a phase difference (i.e., the difference in argument of latitude of two satellites): this formation is the most desirable arrangement to keep the relative position of two satellites. In that case, to reduce the DOP value, a large phase difference is preferable. In comparison, to keep both satellites in the rover’s view for a long time, a small phase difference is desirable. As a result, these two requirements conflict with each other, and both impacts must be carefully considered to find the best compromise point in the satellite trajectory selection. In our comparative analysis, we used two performance index parameters, the T o t a l   G D O P and a v a i l a b i l i t y , where a v a i l a b i l i t y is the percentage of time at which both satellites are in the rover’s view to the total mission time.
The relationship between the T o t a l   G D O P and a v a i l a b i l i t y change depending on the navigation methods, as well as the satellite orbital parameters. Table 2 and Table 3 show the T o t a l   G D O P and a v a i l a b i l i t y comparisons among the three navigation methods under different orbital conditions: in this study, without losing generality, the rover/lander position were fixed to the south-pole ( 90 deg latitude) and the satellite orbit inclination was fixed to 110 deg while the orbital attitudes of two satellites were changed (300 and 2100 km), and the phase differences between the two satellites were also changed (5, 10, 15, and 25 deg).
We found that single-differenced two-way ranging provided the smallest DOP. This is as single-differenced observation has a better quality by nature than double-differenced observation, due to fewer differencing processes. On the other hand, double-differenced TOA–FOA tends to result in the largest DOP due to the low quality of double-differenced pseudodoppler observations. Double-differenced TOA–FOA and single-difference two-way ranging have a larger a v a i l a b i l i t y than the MDPO, as they only require single-epoch observation, while the MDPO requires multiple-epoch observations. The T o t a l   G D O P and a v a i l a b i l i t y changes depending on the attitudes of the satellites, as well as the phase difference between the two satellites. In general, higher orbits provide a larger T o t a l   G D O P and a v a i l a b i l i t y , and larger phase differences provide a smaller T o t a l   G D O P and a v a i l a b i l i t y .

3.7.2. Satellite Orbital Parameters and Systematic Errors Related to DEM

As indicated by Equations (25) and (45), using DEM information in the estimation process induces errors in the X–Y position estimates. There is some correlation between the error and satellite positions: a larger elevation angle from the rover plane to the satellite position tends to lead to a larger X–Y position error: This can be explained by looking at Equation (2). Equation (2) can be reformatted as r R S ( t i ) = | X s ( t i s ) X R ( t i ) + d X R s a | = ( x S ( t i s ) x R ( t i ) + d x R   s a ) 2 + ( y S ( t i s ) y R ( t i ) + d y R   s a ) 2 + ( z S ( t i s ) z R ( t i ) + d z R   s a ) 2 and when the elevation angle is large, z S ( t i s ) z R ( t i ) becomes larger in relation to x S ( t i s ) x R ( t i ) and y S ( t i s ) y R ( t i ) and, consequently, the projection of the Z direction error on the X–Y plane is greater.
The mathematical process R = G d X + w , and its converted form d X = ( G T G ) 1 G T ( R w ) , indicates that the error in the estimates does not appear linearly with respect to the elevation angle, as the process takes the double-differenced or single-differenced form of pseudo-observations but multiplies them by the pseudo-inversed observation matrix, which is also a function of the satellite and user positions. However, we can exploit some useful findings with respect to the satellite orbital parameters: for instance, the higher the satellite altitude at the same orbit inclination, the larger the error. Again, the error cannot be analytically calculated, and we need to use a numerical simulation on a case-by-case basis.
Table 4 shows the user position error with different satellite altitudes, which was calculated using the same numerical simulation as in Section 4.2. The satellite altitude was changed among 300, 600, 900, and 2100 km, while the other parameters were set the same as in Section 4. In this section, satellite orbital positions were created without considering the precise cis-lunar dynamics, and the orbital parameters did not change due to perturbations from the initial set. The result indicates that the user position accuracy deteriorated immediately along with the satellite orbit altitude and that low lunar orbits are suitable for these DEM-aided radio-triangulation-based relative-positioning systems.

4. Comparative Analysis of Three Navigation Methods

As discussed in Section 3.6, the user position accuracy is subject to systematic errors stochastically. Therefore, we require numerical simulation to compare the accuracy of the three navigation methods.

4.1. Overview of Numerical Simulation

The simulation code was initially developed in our prior research [18] and was updated to incorporate the double-differenced TOA–FOA and single-differenced two-way ranging. The simulation code is open to the public and can be accessed at [29].
Figure 5 provides an overview of the simulation system. First, a rover trajectory in the X–Y direction, i.e., a time-series dataset of x R and y R , was created, and then a rover position in the Z direction, i.e., z R , was also created using the lunar DEM data z R   D E M . Then, by adding the DEM error ( d z R   D E M ) to a created rover trajectory, the true rover position X R   t r u e was developed. For lunar DEM data, we used [30], which is 5 m resolution DEM data for latitude from 87.5 deg to 90 deg. The DEM error dataset, i.e., d z R   D E M , was prepared at a 1 m grid interval. In other words, the DEM data changed every 5 m grid, while the DEM error data changed every 1 m grid. The true rover altitude, i.e., the z-component of X R   t r u e , was estimated using the DEM value and DEM error value of the closest grid point from its horizontal location, respectively. For example, if the rover is horizontally located at ( x R , y R ) = (11.3 m, 3.5 m), it refers to the DEM data of the point ( x R , y R ) = (10.0 m, 5.0 m) and the DEM error data of the point ( x R , y R ) = (11.0 m, 3.0 m) to calculate the true rover altitude.
Next, the true satellite trajectory X t r u e S was prepared separately. A precise cis-lunar dynamics model takes into account the gravity model of the Moon of degree 40, as well as the gravity from the Earth and the Sun, which was used to generate satellite trajectory data in the topocentric frame, whose origin is at the lander position.
The true velocity, true range and true doppler were calculated using the true satellite, rover, and lander positions while taking into account the Moon rotation during the signal traveling time between the satellites and rover/lander, i.e., d X R s a and d V R   s a S . Signal processing delay time, i.e., t p r o c e s s i n g S R , was set to zero without losing generality. Then, by adding the receiver observation errors, and the signal processing delay time uncertainty d t p r o c e s s i n g S R for the single-differenced two-way ranging method, to a true range and true doppler, the pseudorange observation ρ R S ( t i ) and pseudodoppler observation ρ d R S ( t i ) were prepared.
By adding satellite orbit determination error d X s a t   O D S and time tag error d X t i m e   t a g S to the satellite’s true position X t r u e S , the observed satellite position X o b S was prepared. Then, the three methods respectively calculate an estimated rover position X R   e s t using the pseudorange observation ρ R S ( t i ) , pseudodoppler observation ρ d R S ( t i ) , observed satellite position X o b S , and lunar DEM data z R   D E M over the course of the simulation period. Finally, the true rover position X R   t r u e and the estimated rover position X R   e s t were compared to evaluate the estimation accuracy.

4.2. Other User-Set Conditions

Table 5 summarizes the general parameters used in the simulation. The total simulation period was set to 15,000 min, assuming a two-week-long mission. The range measurement resolution at the user pseudorange receiver was set to 0.4 m, and the Doppler measurement resolution at the user pseudodoppler receiver was set to 0.2 Hz, assuming a typical space GNSS receiver specification with a conservative safety margin. The initial rover position and lander position were set to ( 90 deg latitude), assuming a south-pole mission. The rover trajectory was created dynamically by changing the rover position after each observation epoch according to the defined traveling distance and the random heading direction specified in Table 5.
MDPO requires pseudorange observations from two epochs, while other navigation methods require pseudorange and/or pseudodoppler observation from a single epoch, and the interval of observations was set to 0.5 min. Hence, it took 1.0 min for the MDPO method to estimate the rover position, and 0.5 min for the other navigation methods to estimate the rover position. The rover position was fixed during the observation epoch(s), then the rover position was changed in the following 0.5 min and then stopped for another observation epoch(s), which continued over the course of the simulation period. In addition, the rover moved only when both orbiters were in view.
Table 6, Table 7, Table 8, Table 9 and Table 10 show the systematic and random error statistics used in the simulation. The realism of these values is discussed as follows: Table 6 summarizes the values of the satellite orbit determination error, i.e., Δ A l o n g ,   Δ R a d i a l , and Δ C r o s s , which are defined in Equation (43). The satellite orbit determination error consists of white noise and systematic error modeled as a sinusoidal function with the period of the satellite orbit. The values were chosen by adding a sufficient margin to the reference data from the Lunar Reconnaissance Orbiter (LRO) project [31]. Table 7 summarizes the value of the time tag error, which is defined as the difference of two user time tags. Time tag error consists of offset and random walk error: the offset component represents a residual time tag error after the frame synchronization. The random walk component was reset to zero and increased until the next frame synchronization. We assumed the frame synchronization takes place in every orbital period. Table 8 summarizes the value of DEM model error, i.e., d z R   D E M . The value was determined based on the actual lunar DEM data by adding a sufficient margin: the accuracy of the best existing DEM data in a vertical direction is about 3 m within a ±60–deg latitude and about 10 m near polar regions [27,28]. The DEM error is derived from calibration errors between multiple sensors and practically consists of white noise and offset according to Figure 4 of [27]. Table 9 summarizes the magnitude of receiver observation errors used in the simulation, i.e., ω r and ω d . Table 10 summarizes the value of signal processing delay time uncertainty, i.e., d t p r o c e s s i n g S R . The signal processing delay time uncertainty consists of white noise and systematic error modeled as a sinusoidal function with the period of the lunar rotation, assuming that the systematic noise is derived from thermal variation of the user radio, which coincides with the lunar thermal environment variation due to the Sun elevation transition of the landing point. Table 11 shows the satellite orbital parameters used in the simulation: two satellites were placed in the 110 deg–300 km (inclination–altitude) circular orbits with 15 deg phase difference. It is important to note that argument of latitude was defined instead of the argument of periapsis and the true anomaly as they were circular orbits. The same parameters were used in the following simulations unless otherwise mentioned.
To secure the statistical accuracy, a Monte Carlo simulation was conducted 100 times, and averaged data are presented for each specific scenario. The rover trajectory and model errors were renewed and created with every simulation.

4.3. Numerical Simulation Result

The simulation results for the three navigation methods are shown in Table 12. Due to the systematic errors discussed in Section 3.6, the T o t a l   U P E ( 2 d r m s ) became larger than the product of the T o t a l   G D O P and twice the standard deviation of differenced receiver observation error, i.e., T o t a l   G D O P × 2 σ Δ ω for the MDPO and double-differenced TOA–FOA with σ Δ ω = 0.4 m and σ Δ ω = 0.2 m, respectively, and T o t a l   G D O P × 2 σ Δ ω for the single-differenced two-way ranging with σ Δ ω = 0.28 m. The deterioration by systematic errors in T o t a l   U P E in relative proportion is greater in the single-differenced two-way ranging than in other two navigation methods under the selected condition.
The three navigation methods provided different T o t a l   G D O P and a v a i l a b i l i t y and, as a result, the T o t a l   U P E and total traveling distance were also different. In general, single-differenced two-way ranging outperformed the other two methods with respect to the user position accuracy due to a smaller T o t a l   G D O P . The total traveling distance of the MDPO was shorter than that of the other two methods due to a longer observation period.
Figure 6 shows examples of the estimated rover trajectory overlaying the true rover trajectory of the three different navigation methods, as well as the distribution of the user position error between the true rover positions and the estimated rover positions. According to Figure 6, under the condition of the selected orbital parameters shown in Table 10, the error distribution does not have a large anisotropy, but may become more anisotropic for other cases, depending on the satellite orbital parameters, initial rover/lander position, and DEM error.

4.4. Numerical Simulation Result with Increased Reveiver Observation Noise

It is important to highlight once again that the user position accuracy is dependent on the magnitude of the receiver observation error. The receiver observation error is caused by several factors including thermal noise, symbol timing offsets in signal processing, and environmental factors. Additionally, the signal processing delay time uncertainty is subject to thermal noise and environmental factors as well. The simulation results with increased receiver observation errors and signal processing delay time uncertainty by a factor of 10, i.e., σ ω r = 2.0 m, σ ω d = 1.0 Hz, the magnitude of ω p r o c e s s is 200.0 ns, and the magnitude of c p r o c e s s is 200.0 ns, as shown in Table 13.
T o t a l   U P E of the three methods were increased by a factor smaller than 10, and the factors were slightly different among the three methods. This is essentially, as mentioned in Section 4.3, due to the systematic errors that affect T o t a l   U P E differently depending on navigation methods.

5. Discussion

In general, the three navigation methods have different characteristics in terms of navigation accuracy and system complexity. Therefore, the system designer must understand the difference and representative performance of these three navigation methods to choose an appropriate method based on the desired specifications.
Through the numerical simulation, we quantitatively confirmed an achievable position accuracy for the three navigation methods under the selected orbital condition. From a user position accuracy point of view, single-differenced two-way ranging outperformed the other two navigation methods. The drawback of the single-differenced two-way ranging is power efficiency, which requires transmitting power at the rover side to reply each receiving radio signal from two satellites. Additionally, the method requires transmitting power at the satellite side to send radio signals to multiple users respectively, i.e., radio signals to multiple rovers and at least one lander.
Based on the simulation results in Section 4.3, if the mission requires a navigation accuracy as high as 30 m, but only for a single rover, and allows a radio signal emission at a rover, then single-differenced two-way ranging is the best choice. On the other hand, if the mission requires the provision of navigation information to multiple users, the MDPO or double-differenced TOA–FOA could be a more efficient option depending on the required accuracy, total traveling distance, and the desired system complexity, i.e., either requiring range sensors or range and Doppler sensors. For instance, the MDPO is the best choice from power efficiency point of view when the mission requires a multi-user navigation system, and the required user position accuracy is about 50 m.
A combination of two navigation methods could be considered to compensate their weaknesses. For example, having the MDPO and single-difference two-way ranging on the same satellite can change the configuration between providing several tens of meters of navigation accuracy to multiple users, or providing higher than thirty meters of navigation accuracy to a single user, depending on the mission needs, without launching another set of satellites.
Single-differenced two-way ranging can also be achieved with a single satellite using multi-epoch observation at the expense of availability. We will study that in our future work.

6. Conclusions

In this paper, we studied and compared three dual-satellite lunar navigation systems that consist of a constellation of two navigation satellites. Dual-satellite navigation systems play key roles in establishing a low-cost navigation platform around the Moon. While several dual-satellite navigation methods have been studied, we focused on the comparison of three navigation methods, MDPO, double-differenced TOA–FOA, and single-differenced two-way ranging, as these three methods represent three different types in terms of observation data, i.e., passive ranging, passive ranging and doppler, and active ranging, into which most dual-satellite navigation methods can be classified.
First, we derived the mathematical models of these three methods step by step, to clarify the differences among the three navigation methods. Next, we confirmed the achievable user position accuracy of the three navigation methods by numerical simulation under the selected orbital conditions. Based on the numerical simulation results, we discussed the advantages and disadvantages of the three navigation methods and provided a guideline to select one or a combination of these three navigation methods depending on the mission requirements: From a user position accuracy point of view, single-differenced two-way ranging outperformed the other two navigation methods, while single-differenced two-way ranging requires larger power consumption at the rover side as well as at the satellite side. If the mission requires the provision of navigation information to multiple users, the MDPO or double-differenced TOA–FOA could be a more efficient option, depending on the desired specifications such as required accuracy, total traveling distance, and the desired system complexity, i.e., either requiring range sensors or range and Doppler sensors.
Furthermore, a combination of two navigation methods could be considered to compensate their weaknesses. For instance, having the MDPO and single-difference two-way ranging on the same satellite enables the users to choose between two configurations, providing several tens of meters of navigation accuracy to multiple users, or providing higher than thirty meters of navigation accuracy to a single user, without launching another set of satellites.

Author Contributions

Conceptualization, T.T., T.E. and S.N.; methodology, T.T.; software, T.T.; validation, T.T., T.E. and S.N.; formal analysis, T.T.; investigation, T.T., T.E. and S.N.; resources, T.T.; data curation, T.T.; writing—original draft preparation, T.T.; writing—review and editing, T.T., T.E., S.N. and H.M.; visualization, T.T.; supervision, T.T., T.E., S.N. and H.M.; project administration, T.T., T.E., S.N. and H.M.; and funding acquisition, N/A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

Yosuke Kawabata (Intelligent Space Systems Laboratory, The University of Tokyo) and Keidai Iiyama (same) are kindly acknowledged for their technical supports and advice.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Flanegan, M.; Gal-Edd, J.; Anderson, L.; Warner, J.; Ely, T.A.; Lee, C.; Shah, B.; Vaisnys, A.; Schier, J. NASA’s Lunar Communication and Navigation Architecture. In Proceedings of the SpaceOps 2008 Conference, Heidelberg, Germany, 12–16 May 2008. [Google Scholar] [CrossRef]
  2. Laîné, M.; Yoshida, K. Multi-Rover Exploration Strategies: Coverage Path Planning with Myopic Sensing. Ph.D. Thesis, Tohoku University, Sendai, Japan, 2019. [Google Scholar] [CrossRef]
  3. Moreau, M.C.; Davis, E.P.; Carpenter, J.R.; Kelbel, D.; Davis, G.W.; Axelrad, P. Results from the GPS Flight Experiment on the High Earth Orbit AMSAT OSCAR-40 Spacecraft. In Proceedings of the 15th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GPS 2002), Portland, OR, USA, 24–27 September 2002; pp. 122–133. [Google Scholar]
  4. Palmerini, G.B.; Sabatini, M.; Perrotta, G. En route to the Moon using GNSS signals. Acta Astronaut. 2009, 64, 467–483. [Google Scholar] [CrossRef]
  5. Witternigg, N.; Obertaxer, G.; Schönhuber, M.; Palmerini, G.B.; Rodriguez, F.; Capponi, L.; Soualle, F.; Floch, J. Weak GNSS Signal Navigation for Lunar Exploration Missions. In Proceedings of the 28th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS + 2015), Tampa, FL, USA, 14–18 September 2015; pp. 3928–3944. [Google Scholar]
  6. Ashman, B.W.; Parker, J.K.; Bauer, F.H.; Esswein, M. Exploring the Limits of High Altitude GPS for Future Lunar Missions. In Proceedings of the 41st Annual AAS Guidance Control Conference (AAS/GNC), Breckenridge, CO, USA, 1–7 February 2018. [Google Scholar]
  7. Delépaut, A.; Giordano, P.; Ventura-Traveset, J.; Blonski, D.; Schönfeldt, M.; Schoonejans, P.; Aziz, S.; Walker, R. Use of GNSS for Lunar Missions and Plans for Lunar In-Orbit Development. Adv. Space Res. 2020, 66, 2739–2756. [Google Scholar] [CrossRef]
  8. Capuano, V.; Botteron, C.; Leclère, J.; Tian, J.; Wang, Y.; Farine, P. Feasibility study of GNSS as navigation system to reach the Moon. Acta Astronaut. 2015, 116, 186–201. [Google Scholar] [CrossRef]
  9. Capuano, V.; Blunt, P.; Botteron, C.; Tian, J.; Leclère, J.; Wang, Y.; Basile, F.; Farine, P.-A. Standalone GPS L1 C/A Receiver for Lunar Missions. Sensors 2016, 16, 347. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Manzano-Jurado, M.; Alegre-Rubio, J.; Pellacani, A.; Seco-Granados, G.; López-Salcedo, J.A.; Guerrero, E.; García-Rodríguez, A. Use of weak GNSS signals in a mission to the moon. In Proceedings of the 2014 7th ESA Workshop on Satellite Navigation Technologies and European Workshop on GNSS Signals and Signal Processing (NAVITEC), Noordwijk, The Netherlands, 3–5 December 2014; pp. 1–8. [Google Scholar] [CrossRef]
  11. Chen, H.; Liu, J.; Long, L.; Xu, Z.; Meng, Y.; Zhang, H. Lunar far side positioning enabled by a CubeSat system deployed in an Earth-Moon halo orbit. Adv. Space Res. 2019, 64, 28–41. [Google Scholar] [CrossRef]
  12. Iiyama, K. Optimization of Navigation Satellite Constellation and Lunar Monitoring Station Arrangement for Lunar Global Navigation Satellite System (LGNSS). In Proceedings of the 32nd International Symposium on Space Technology and Science (ISTS), Fukui, Japan, 15–21 June 2019. [Google Scholar]
  13. Schonfeldt, M.; Grenier, A.; Delépaut, A.; Giordano, P.; Swinden, R.; Ventura-Traveset, J.; Blonski, D.; Hahn, J. A System Study about a Lunar Navigation Satellite Transmitter System. In Proceedings of the 2020 European Navigation Conference (ENC), Dresden, Germany, 23–24 November 2020; pp. 1–10. [Google Scholar] [CrossRef]
  14. Guo, F.; Fan, Y.; Zhou, Y.; Xhou, C.; Li, Q. Space Electronic Reconnaissance: Localization Theories and Methods; WILEY: New York, NY, USA, 2014. [Google Scholar]
  15. Shilong, W.; Jingqing, L.; Liangliang, G. Joint FDOA and TDOA location algorithm and performance analysis of dual-satellite formations. In Proceedings of the 2010 2nd International Conference on Signal Processing Systems (ICSPS), Dalian, China, 5–7 July 2010; pp. V2-339–V2-342. [Google Scholar] [CrossRef]
  16. Jun, W.W.; Cheung, K.; Lightsey, E.G.; Lee, C. Localizing in Urban Canyons using Joint Doppler and Ranging and the Law of Cosines Method. In Proceedings of the 32nd International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2019), Miami, FL, USA, 16–20 September 2019; pp. 140–153. [Google Scholar] [CrossRef]
  17. Jun, W.; Cheung, K.-M.; Milton, J.; Lee, C.; Lightsey, G. Autonomous Navigation for Crewed Lunar Missions with DBAN. In Proceedings of the 2020 IEEE Aerospace Conference, Big Sky, MT, USA, 7–14 March 2020; pp. 1–13. [Google Scholar] [CrossRef]
  18. Tanaka, T.; Ebinuma, T.; Nakasuka, S. Dual-Satellite Lunar Global Navigation System Using Multi-Epoch Double-Differenced Pseudorange Observations. Aerospace 2020, 7, 122. [Google Scholar] [CrossRef]
  19. Music, M.; Newberry, L.; Peters, D.; Quetin, G.; Stickle, A.; Trescott, R.; Ackerman, K.; Bruun, E.; Cleveland, T.; Culver, T. Luna POLARIS A Lunar Positioning and Communications System. Available online: https://www.researchgate.net/publication/268005025_Luna_POLARIS_A_Lunar_Positioning_and_Communications_System (accessed on 25 June 2021).
  20. Ely, T.A.; Burt, E.A.; Prestage, J.D.; Seubert, J.M.; Tjoelker, R.L. Using the Deep Space Atomic Clock for Navigation and Science. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 2018, 65, 950–961. [Google Scholar] [CrossRef] [PubMed]
  21. Cheung, K.; Lee, C. Differencing Methods for 3D Positioning of Spacecraft. In Proceedings of the 2019 Integrated Communications, Navigation and Surveillance Conference (ICNS), Herndon, VA, USA, 9–11 April 2019; pp. 1–10. [Google Scholar] [CrossRef]
  22. NASA Article: Commercial CubeRover Test Shows How NASA Investments Mature Space Tech. Available online: https://www.nasa.gov/feature/commercial-cuberover-test-shows-how-nasa-investments-mature-space-tech (accessed on 25 June 2021).
  23. Walker, J. Flight System Architecture of the Sorato Lunar Rover. In Proceedings of the International Symposium on Artificial Intelligence, Robotics and Automation in Space (i-SAIRAS 2018), Madrid, Spain, 4–6 June 2018. [Google Scholar]
  24. Rybak, M.M.; Axelrad, P.; Seubert, J.; Ely, T. Chip Scale Atomic Clock–Driven One-Way Radiometric Tracking for Low-Earth-Orbit CubeSat Navigation. J. Spacecr. Rocket. 2021, 58, 200–209. [Google Scholar] [CrossRef]
  25. The LROC Team Computed the Coordinates of China’s Chang’e 5 Lander. Available online: http://lroc.sese.asu.edu/posts/1172?fbclid=IwAR3fNoQuFQEmT6vSGh24yGPyRJAmpIfOln64IJCJg0bxC_A_dV6DNaCLJz0 (accessed on 25 June 2020).
  26. Hennawy, J.; Adams, N.; Sanchez, E.; Srinivasan, D.; Hamkins, J.; Vilnrotter, V.; Xie, H.; Kinman, P. Telemetry ranging using software-defined radios. In Proceedings of the 2015 IEEE Aerospace Conference, Big Sky, MT, USA, 7–14 March 2015; pp. 1–14. [Google Scholar] [CrossRef]
  27. Barker, M.K.; Mazarico, E.; Neumann, G.A.; Zuber, M.T.; Haruyama, J.; Smith, D.E. A new lunar digital elevation model from the Lunar Orbiter Laser Altimeter and SELENE Terrain Camera. Icarus 2016, 273, 346–355. [Google Scholar] [CrossRef] [Green Version]
  28. Smith, D.E.; Zuber, M.T.; Neumann, G.A.; Mazarico, E.; Lemoine, F.G.; Head, J.W., III; Lucey, P.G.; Aharonson, O.; Robinson, M.S.; Sun, X.; et al. Summary of the results from the lunar orbiter laser altimeter after seven years in lunar orbit. Icarus 2017, 283, 70–91. [Google Scholar] [CrossRef] [Green Version]
  29. GitHub Repository. Available online: https://github.com/tskTNK/DualSatLunarNav (accessed on 25 June 2020).
  30. LRO Data Products. Available online: https://lunar.gsfc.nasa.gov/dataproducts.html (accessed on 25 June 2020).
  31. Mazarico, E.; Neumann, G.A.; Barker, M.K.; Goossens, S.; Smith, D.E.; Zuber, M.T. Orbit determination of the Lunar Reconnaissance Orbiter. J. Geod. 2012, 86, 193–207. [Google Scholar] [CrossRef]
Figure 1. Overview of the multi-epoch double-differenced pseudorange observation (MDPO) method.
Figure 1. Overview of the multi-epoch double-differenced pseudorange observation (MDPO) method.
Aerospace 08 00191 g001
Figure 2. Overview of the double-differenced TOA–FOA method.
Figure 2. Overview of the double-differenced TOA–FOA method.
Aerospace 08 00191 g002
Figure 3. Two-way ranging.
Figure 3. Two-way ranging.
Aerospace 08 00191 g003
Figure 4. Overview of the single-differenced two-way ranging method.
Figure 4. Overview of the single-differenced two-way ranging method.
Aerospace 08 00191 g004
Figure 5. Simulation overview.
Figure 5. Simulation overview.
Aerospace 08 00191 g005
Figure 6. The simulated rover trajectories and user position errors. (a,b) correspond to the MDPO, (c,d) correspond to the double-differenced TOA–FOA, (e,f) correspond to the single-differenced two-way ranging.
Figure 6. The simulated rover trajectories and user position errors. (a,b) correspond to the MDPO, (c,d) correspond to the double-differenced TOA–FOA, (e,f) correspond to the single-differenced two-way ranging.
Aerospace 08 00191 g006
Table 1. Benchmark of dual-satellite lunar navigation systems.
Table 1. Benchmark of dual-satellite lunar navigation systems.
MethodObservation TypeNumber of Supported UsersObservation TimeNavigation Accuracy
Multi-epoch double-differenced pseudorange observationPassive ranging Multi-userUsing observations from at least two epochs 250 m under the condition used in Section 4.3.
Double-differenced time of arrival (TOA)-frequency of arrival (FOA)Passive ranging and Doppler Multi-userUsing observations from a single epoch 3100 m under the condition used in Section 4.3.
Single-differenced two-way ranging Active rangingSingle-user at a time 1Using observations from a single epoch 330 m under the condition used in Section 4.3.
1 The second user needs another set of radio signals separately. 2 For example, 1 min (0.5 min × 2 epochs) such as set in this research. 3 For example, 0.5 min (0.5 min × 1 epoch) such as set in this research.
Table 2. The Total GDOP and availability comparison among three navigation methods under different orbital conditions: two satellites are placed in 300 km circular low lunar orbit with different phase differences (5, 10, 15, and 25 deg). The rover/lander position were fixed to the south pole ( 90 deg latitude), and the satellite orbit inclination was fixed to 110 deg.
Table 2. The Total GDOP and availability comparison among three navigation methods under different orbital conditions: two satellites are placed in 300 km circular low lunar orbit with different phase differences (5, 10, 15, and 25 deg). The rover/lander position were fixed to the south pole ( 90 deg latitude), and the satellite orbit inclination was fixed to 110 deg.
Navigation MethodsPhase Difference [deg]Total GDOPAvailability [%]
MDPO5217.84.1
1063.13.2
1538.52.3
2534.31.4
Double-differenced TOA–FOA5891.96.2
10253.54.7
15159.53.3
25144.62.1
Single-differenced Two-way Ranging58.16.2
102.54.7
151.63.3
251.42.1
Table 3. The Total GDOP and availability comparison among three navigation methods under different orbital conditions: two satellites were placed in 2100 km circular low lunar orbit with different phase differences (5, 10, 15, and 25 deg). The rover/lander positions were fixed to the south pole ( 90 deg latitude), and the satellite orbit inclination was fixed to 110 deg.
Table 3. The Total GDOP and availability comparison among three navigation methods under different orbital conditions: two satellites were placed in 2100 km circular low lunar orbit with different phase differences (5, 10, 15, and 25 deg). The rover/lander positions were fixed to the south pole ( 90 deg latitude), and the satellite orbit inclination was fixed to 110 deg.
Navigation MethodsPhase Difference [deg]Total GDOP Availability [%]
MDPO51336.610.9
10447.810.0
15274.29.0
25203.68.1
Double-differenced TOA–FOA55669.016.3
101899.015.0
151162.913.4
25863.512.1
Single-differenced Two-way Ranging521.616.3
106.715.0
153.813.4
252.712.1
Table 4. One example of correlation between the satellite altitude and Total UPE: two satellites were placed in different circular LLO, 300, 600, 900, and 2100 km. The simulation takes an average of 100 simulation cases for each.
Table 4. One example of correlation between the satellite altitude and Total UPE: two satellites were placed in different circular LLO, 300, 600, 900, and 2100 km. The simulation takes an average of 100 simulation cases for each.
Navigation MethodsSatellite Altitude (km)Total UPE (2drms) (m)
MDPO30057.9
600121.4
900197.3
21001038.3
Double-differenced TOA–FOA300119.8
600245.7
900454.3
21002049.6
Single-differenced Two-way Ranging30027.0
60040.1
90056.3
2100217.8
Table 5. The simulation parameters.
Table 5. The simulation parameters.
ItemsValueUnitRemarks
Simulation Period15,000minApproximately two weeks in Earth time.
Range measurement resolution of the user pseudorange receivers0.4mMinimum observable resolution by the rover and lander receivers.
Doppler measurement resolution of the user pseudodoppler receivers0.2HzMinimum observable resolution by the rover and lander receivers.
Latitude of initial rover/lander position 90deg
Interval of pseudorange/doppler observations0.5min
Rover traveling distance between observations3.75mThe rover travels at 7.5 m/min for 0.5 min between position estimations.
Rover traveling directionRandomdegThe heading direction is selected from three values (+ π 3 , π 3 ,   0 ) randomly.
Table 6. Overview of the satellite orbit determination error used in the simulation.
Table 6. Overview of the satellite orbit determination error used in the simulation.
ItemsTypeValueUnitRemarks
Satellite Orbit Determination Error in the Along Direction d A l o n g ( t i ) = ω O D A l o n g   ( t i ) + c O D A l o n g
White Gaussian random error ω O D A l o n g   100.0m ω O D   t = V a l u e × a random scalar drawn from the standard normal distribution each time.
Systematic error c O D A l o n g   200.0mSystematic error c O D is an output of the sinusoidal function A × s i n ( 2 π x / T ) : the argument x is epoch time, the period T was set equal to the satellite orbital period, and the amplitude A is randomly selected between V a l u e and V a l u e at the beginning of each simulation.
Satellite Orbit Determination Error in the Radial Direction d R a d i a l ( t i ) = ω O D R a d i a l ( t i ) + c O D R a d i a l
White Gaussian random error ω O D R a d i a l   10.0mSame as above.
Systematic error c O D R a d i a l   20.0m
Satellite Orbit Determination Error in the Cross Direction d C r o s s ( t i ) = ω O D C r o s s ( t i ) + c O D C r o s s
White Gaussian random error ω O D C r o s s   100.0mSame as above.
Systematic error c O D C r o s s 200.0m
Table 7. Overview of the time tag error used in the simulation.
Table 7. Overview of the time tag error used in the simulation.
ItemTypeValueUnitRemarks
Time Tag Error   d τ R ( t i ) = c t i m e   t a g + x t i m e   t a g  
Offset error c t i m e   t a g   1.0msOffset error c t i m e   t a g is randomly selected between V a l u e and V a l u e after the time synchronization and fixed until the next time synchronization. The time synchronization takes place in every orbital period.
Random walk x   t i m e   t a g   1.0 × 10 8 ms/minA random walk is a time series model x t i m e   t a g   ( t ) such that x t i m e   t a g   ( t ) = x t i m e   t a g   ( t 1 ) + ω t where ω t is a discrete white noise series. Random walk noise is reset to zero after the time synchronization and increases until the next time synchronization. The time synchronization takes place in every orbital period.
Table 8. Overview of the DEM error used in the simulation.
Table 8. Overview of the DEM error used in the simulation.
ItemTypeValueUnitRemarks
DEM Error d z R   D E M = ω D E M + c D E M
White Gaussian random error ω D E M   10.0m ω D E M   t = V a l u e × a random scalar drawn from the standard normal distribution each time.
Offset error c D E M 5.0mOffset error c D E M is randomly selected between V a l u e and V a l u e at the beginning of each simulation and fixed during the simulation.
Table 9. Overview of the receiver observation error used in the simulation.
Table 9. Overview of the receiver observation error used in the simulation.
ItemTypeValueUnitRemarks
Receiver Observation ErrorRange white Gaussian random error ω r 0.2m ω r = V a l u e × a random scalar drawn from the standard normal distribution each time, i.e., σ ω r   = 0.2 m.
Doppler white Gaussian random error ω d 0.1Hz ω d = V a l u e × a random scalar drawn from the standard normal distribution each time, i.e., σ ω d   = 0.1 Hz.
Table 10. Overview of the signal processing delay time uncertainty used in the simulation.
Table 10. Overview of the signal processing delay time uncertainty used in the simulation.
ItemTypeValueUnitRemarks
Signal Processing Delay Time Uncertainty d t p r o c e s s i n g S R = ω p r o c e s s + c p r o c e s s
White Gaussian random error ω p r o c e s s   20.0ns ω D E M   t = V a l u e × a random scalar drawn from the standard normal distribution each time.
Systematic error c p r o c e s s   20.0nsSystematic error c p r o c e s s is an output of the sinusoidal function A × s i n ( 2 π x / T ) : the argument x is epoch time, the period T was set equal to the lunar rotation period, and the amplitude A is randomly selected between V a l u e and V a l u e at the beginning of each simulation.
Table 11. The satellite orbital parameters used in the simulation.
Table 11. The satellite orbital parameters used in the simulation.
ItemsValueUnit
Satellite 1 perilune altitude300km
Satellite 1 apolune altitude300km
Satellite 1 inclination110deg
Satellite 1 right ascension of the ascending node0deg
Satellite 1 argument of latitude0deg
Satellite 2 perilune altitude300km
Satellite 2 apolune altitude300km
Satellite 2 inclination110deg
Satellite 2 right ascension of the ascending node0deg
Satellite 2 argument of latitude−15deg
Table 12. The numerical simulation results taking an average of 100 simulation cases.
Table 12. The numerical simulation results taking an average of 100 simulation cases.
Navigation MethodsTotal GDOP Total UPE (2drms) (m)Availability (%)Total Traveling Distance (m)
MDPO46.755.33.33753.75
Double-differenced TOA–FOA193.8109.95.05625
Single-differenced Two-way Ranging1.226.35.05625
Table 13. The numerical simulation results with increased receiver observation errors and signal processing delay time uncertainty: σ ω r = 2.0 m, σ ω d = 1.0 Hz, the magnitude of ω p r o c e s s is 200.0 ns, and the magnitude of c p r o c e s s is 200.0 ns.
Table 13. The numerical simulation results with increased receiver observation errors and signal processing delay time uncertainty: σ ω r = 2.0 m, σ ω d = 1.0 Hz, the magnitude of ω p r o c e s s is 200.0 ns, and the magnitude of c p r o c e s s is 200.0 ns.
Navigation MethodsTotal GDOP Total UPE (2drms) (m)Availability (%)Total Traveling Distance (m)
MDPO46.7437.93.33753.75
Double-differenced TOA–FOA193.8922.15.05625
Single-differenced Two-way Ranging1.2249.95.05625
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tanaka, T.; Ebinuma, T.; Nakasuka, S.; Malki, H. A Comparative Analysis of Multi-Epoch Double-Differenced Pseudorange Observation and Other Dual-Satellite Lunar Global Navigation Systems. Aerospace 2021, 8, 191. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace8070191

AMA Style

Tanaka T, Ebinuma T, Nakasuka S, Malki H. A Comparative Analysis of Multi-Epoch Double-Differenced Pseudorange Observation and Other Dual-Satellite Lunar Global Navigation Systems. Aerospace. 2021; 8(7):191. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace8070191

Chicago/Turabian Style

Tanaka, Toshiki, Takuji Ebinuma, Shinichi Nakasuka, and Heidar Malki. 2021. "A Comparative Analysis of Multi-Epoch Double-Differenced Pseudorange Observation and Other Dual-Satellite Lunar Global Navigation Systems" Aerospace 8, no. 7: 191. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace8070191

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop