High-Accuracy Positioning Based on Pseudo-Ranges: Integrated Difference and Performance Analysis of the Loran System

The Long Range Navigation (Loran) system as a backup of the Global Navigation Satellite System (GNSS) is a good choice. The dominant deterioration factors of position accuracy are the pseudo-range measurement errors and the geometric dilution of precision (GDOP). This paper focuses on the algorithm integrated difference with pseudo-ranges to improve the position accuracy. Firstly, the theoretical prediction of propagation delay and raw measurement are compared. The results show that the measured pseudo-range consists of a constant term and a temporal term, which reflect the propagation situation along the true path. Secondly, a position solution algorithm based on a pseudo-range and difference is presented, exceeding the limit of a single chain. Finally, some simulation tests are implemented utilizing the new proposed position algorithm to verify the differential performance. This method can reduce the GDOP conveniently through increasing the number of transmitters. In view of the amplitude and characteristics of errors in measurement, systematic error and random noise are distinguished and discussed. The absolute accuracy responds to the pseudo-range bias that is different from geometric distance and repeatable accuracy is mainly influenced by random noise. The difference method can improve the absolute accuracy via the correction degree without changing the geometry of the transmitters.


Introduction
Accurate position is a significant topic that should be intensively studied. Generally, it is easy to get the required accuracy using the Global Navigation Satellite System (GNSS). However, the vulnerability of GNSS to unintentional and intentional interference signals can be found frequently nowadays [1][2][3]. For national security and economic effectiveness, a reliable and complementary navigation system is needed desperately [4]. The suitability of the Loran for a backup navigation system has been evaluated and reported [5]. Since then, as a leading role, the United States Coast Guard (USCG) has tried to improve the performance of the Loran system dramatically by the modernization of equipment [6], which mainly included the improvement of transmitters and reference clocks. The improved effect and application are also very obvious [7]. Besides the US, the United Kingdom (UK), South Korea and other countries have also introduced much research on alternatives to the GNSS backup system [8][9][10][11]. They all identified Loran as the terrestrial radio navigation system that has the potential to fulfil those requirements. China has also carried out relevant research [12][13][14]. Meanwhile, China's major national infrastructure project for developing terrestrial Loran began recently. Three new transmitters using the Loran-C mechanism will be built in the west region of China, which will construct a full coverage timing and positioning network. In order to meet the stated performance [15] as a backup, Generally, it is important to convert the propagation delay to the ground truth range in positioning [16]. If the truth range is known, the propagation delay is derived from the distance divided by the speed of light in a vacuum. However, the Loran signal propagates out radially from the transmitter by groundwave, traveling parallel to the surface of the Earth. So, it does not travel at the free space speed of light, rather, it is slowed by the atmosphere and the surface of the Earth. The groundwave accumulates a delay compared to the expected propagation time over the same distance in free space. This delay is the accumulation of three theoretical components or factors [16][17][18][19][23][24][25]. One part is the PF, which is the time delay of the long wave signal propagation through the atmosphere. The second part is the seawater SF, which is the additional delay due to the signal traveling over the sea. The last one is the ASF, the additional delay due to the signal traveling over the land. The specific formula can be seen in Equation (1).
where PF = R d * n s C , the light speed C in the vacuum is 299,792,458 m/s and R d is the distance of signal propagation between the transmitter and the receiver. The refractive index of the atmosphere n s , Sensors 2020, 20, 4436 3 of 15 assumed as a constant, means that the speed of the signal is a fraction lower than the speed of light in a vacuum. At present, it is difficult to distinguish the SF and the ASF, so SF is often incorporated with ASF as a function of distance R d , carrier frequency f , relative dielectric ε r and equivalent conductivity σ e , as shown in Equation (2): where ω = 2π f corresponds to the Loran carrier frequency of 100 kHz and argW is the phase of signal attenuation function W about the Loran. It is known that the value of SF + ASF is a smaller modified term comparing to PF [30]. According to the parameters in [30], SF + ASF via distance for each type of ground is given in Figure 1. When the propagation distance increases the SF + ASF also increases. At a constant distance, SF + ASF is the least for average sea water.
Sensors 2020, 20, x FOR PEER REVIEW 3 of 15 , assumed as a constant, means that the speed of the signal is a fraction lower than the speed of light in a vacuum. At present, it is difficult to distinguish the and the ，so is often incorporated with as a function of distance , carrier frequency , relative dielectric and equivalent conductivity , as shown in Equation (2): where = 2 corresponds to the Loran carrier frequency of 100 kHz and is the phase of signal attenuation function about the Loran. It is known that the value of + is a smaller modified term comparing to [30]. According to the parameters in [30], + via distance for each type of ground is given in Figure 1. When the propagation distance increases the + also increases. At a constant distance, + is the least for average sea water. In the above model, the total propagation delay is calculated using some simplification, similarly to the result of Brunavs' formula. The result of the calculation is constant if the ground truth distance is known and the dielectric constant with equivalent conductivity is determined. At this level, the effect of climate and weather is not taken into account, or of ground elevation on the index of refraction of the atmosphere on the surface of the ground. Meanwhile, it is not enough to consider the dielectric and equivalent conductivity changes due to freezing rain and snowy weather. These variations are often lumped together for the user. Real-time factors, such as climate and weather, will result in a distinct difference between predicted value and true propagation delay. That is to say, the total propagation delay is not accurate, although it is possible to predict it from Equations (1) and (2).

Measured Propagation Delay and Assessment
In order to describe the true value, a modified formula is presented: where ∆ ( ) accounts for all of the time varying aspects, reflecting the impact of various real-time factors. In other words, ∆ ( ) depends on the weather and ground conductivity along the propagation path over time. Any changes in climate and the conductivity would change ∆ ( ). If the receiver position is not accurate or the propagation path is irregular, the fluctuations would also be contained in ∆ ( ). Similarly, we also divided the propagation delay into two parts: a spatial component , which is the constant part of ( ), and a temporal component ∆ ( ). A better way to get propagation delay is to measure it directly. For this purpose, a Loran receiver is used. The measured propagation delay is the difference between the Time of Arrival ( ) and In the above model, the total propagation delay is calculated using some simplification, similarly to the result of Brunavs' formula [22]. The result of the calculation is constant if the ground truth distance is known and the dielectric constant with equivalent conductivity is determined. At this level, the effect of climate and weather is not taken into account, or of ground elevation on the index of refraction of the atmosphere on the surface of the ground. Meanwhile, it is not enough to consider the dielectric and equivalent conductivity changes due to freezing rain and snowy weather. These variations are often lumped together for the user. Real-time factors, such as climate and weather, will result in a distinct difference between predicted value and true propagation delay. That is to say, the total propagation delay is not accurate, although it is possible to predict it from Equations (1) and (2).

Measured Propagation Delay and Assessment
In order to describe the true value, a modified formula is presented: where ∆ASF(t) accounts for all of the time varying aspects, reflecting the impact of various real-time factors. In other words, ∆ASF(t) depends on the weather and ground conductivity along the propagation path over time. Any changes in climate and the conductivity would change ∆ASF(t).
If the receiver position is not accurate or the propagation path is irregular, the fluctuations would also be contained in ∆ASF(t). Similarly, we also divide the propagation delay into two parts: a spatial component T p , which is the constant part of T p (t), and a temporal component ∆ASF(t).
A better way to get propagation delay is to measure it directly. For this purpose, a Loran receiver is used. The measured propagation delay N is the difference between the Time of Arrival (TOA) Time of Transmission ( ) of Loran pulses coherently modulated on the carrier ground wave. The measurement principle is shown in Figure 2.
where and are the clock bias and the internal processing delays of the receiver, respectively, and ( ) represents the random residual error due to the unavoidable noise in artificial measurement. In the Loran frequency band, the dominant source of noise is atmospheric noise, which is caused by lighting discharges and man-made noise and local interference from, for example, switch-mode power supplies. Other sources of noise may include transmitter pulse timing jitter or receiver-related noise [31]. It is likely the accuracy of itself is a function of , consisting of a constant part ′ (i.e., ′ = + + ) different from and a temporal part ∆ ( ) + ( ). Figure 3 presents a raw measurement of in a fixed position (109.7503°E, 38.2503°N) for four days, receiving the signal of a BPL transmitter (109.5431°E, 34.9486°N) in Pu Cheng (PC, a city in Shaanxi Province, China). Here, abnormal measurements are removed from the results in advance. The propagation environment is rather homogeneous. Obviously, varys via time, as seen from the blue line in Figure 3a. A mean value is also plotted by the red line, giving the constant part of . The average value of the observation is 1253.145 µ s, while the predicted constant is 1223.860 µ s based on the parameters of average land [30] using Equations (1) and (2). Depending upon the navigation application, the difference originating from the effect of + is negligible, because each transmitter is synchronized and the user receiver clock bias is common. Due to the influence of random noise on raw measurements, polynomial fitting is used to derive the trend of observation, seen as the green curve in Figure 3a.  In order to compare the amplitude of each component, Figure 3b gives the decomposition of ′ , ∆ ( ) and ( ). It is clear that ′ , the mean value over four days, is the determinant of position and ∆ ( ) is much smaller than the former. However, its variation as time advanced is no more τ in TOT(τ) represents the transmitter's frame and r in TOA(r) is the receiver time frame. The propagation delay N between receiver and transmitter is defined by Equation (4) below: where δT r and P r are the clock bias and the internal processing delays of the receiver, respectively, and ξ(t) represents the random residual error due to the unavoidable noise in artificial measurement.
In the Loran frequency band, the dominant source of noise is atmospheric noise, which is caused by lighting discharges and man-made noise and local interference from, for example, switch-mode power supplies. Other sources of noise may include transmitter pulse timing jitter or receiver-related noise [31]. It is likely the accuracy of N itself is a function of t, consisting of a constant part N (i.e., N = T p + δT r + P r ) different from T p and a temporal part ∆ASF(t) + ξ(t). The propagation environment is rather homogeneous. Obviously, N varies via time, as seen from the blue line in Figure 3a. A mean value is also plotted by the red line, giving the constant part of N. The average value of the observation is 1253.145 µs, while the predicted constant T p is 1223.860 µs based on the parameters of average land [30] using Equations (1) and (2). Depending upon the navigation application, the difference originating from the effect of δT r + P r is negligible, because each transmitter is synchronized and the user receiver clock bias is common. Due to the influence of random noise on raw measurements, polynomial fitting is used to derive the trend of observation, seen as the green curve in Figure 3a. in ( ) represents the transmitter's frame and in ( ) is the receiver time frame. The propagation delay between receiver and transmitter is defined by Equation (4) below: where and are the clock bias and the internal processing delays of the receiver, respectively, and ( ) represents the random residual error due to the unavoidable noise in artificial measurement. In the Loran frequency band, the dominant source of noise is atmospheric noise, which is caused by lighting discharges and man-made noise and local interference from, for example, switch-mode power supplies. Other sources of noise may include transmitter pulse timing jitter or receiver-related noise [31]. It is likely the accuracy of itself is a function of , consisting of a constant part ′ (i.e., ′ = + + ) different from and a temporal part ∆ ( ) + ( ). Figure 3 presents a raw measurement of in a fixed position (109.7503°E, 38.2503°N) for four days, receiving the signal of a BPL transmitter (109.5431°E, 34.9486°N) in Pu Cheng (PC, a city in Shaanxi Province, China). Here, abnormal measurements are removed from the results in advance. The propagation environment is rather homogeneous. Obviously, varys via time, as seen from the blue line in Figure 3a. A mean value is also plotted by the red line, giving the constant part of . The average value of the observation is 1253.145 µ s, while the predicted constant is 1223.860 µ s based on the parameters of average land [30] using Equations (1) and (2). Depending upon the navigation application, the difference originating from the effect of + is negligible, because each transmitter is synchronized and the user receiver clock bias is common. Due to the influence of random noise on raw measurements, polynomial fitting is used to derive the trend of observation, seen as the green curve in Figure 3a. In order to compare the amplitude of each component, Figure 3b gives the decomposition of ′ , ∆ ( ) and ( ). It is clear that ′ , the mean value over four days, is the determinant of position and ∆ ( ) is much smaller than the former. However, its variation as time advanced is no more In order to compare the amplitude of each component, Figure 3b gives the decomposition of N , ∆ASF(t) and ξ(t). It is clear that N , the mean value over four days, is the determinant of position and ∆ASF(t) is much smaller than the former. However, its variation as time advanced is no more than 30 ns no matter whether it increased or decreased. The amplitude of ξ(t) is similar to ∆ASF(t), whose histogram of each day is plotted in Figure 4. The standard deviations (STDs) of calculation for ξ(t) are 10.343 ns, 13.521 ns, 11.463 ns and 9.647 ns, respectively. Besides, normal distribution fitting curves with the same STDs are also presented, which are used to explain the characteristic of ξ(t).
Sensors 2020, 20, x FOR PEER REVIEW 5 of 15 than 30 ns no matter whether it increased or decreased. The amplitude of ( ) is similar to ∆ ( ), whose histogram of each day is plotted in Figure 4. The standard deviations (STDs) of calculation for ( ) are 10.343 ns, 13.521 ns, 11.463 ns and 9.647 ns, respectively. Besides, normal distribution fitting curves with the same STDs are also presented, which are used to explain the characteristic of ( ).

Position Solution Method Integrated Difference
The difference method is usually used to correct common errors between the user and referent station. A scheme of the difference correction system in this paper is shown in Figure 5. Based on the above illustration, the procedure using difference is:  Consistency test of reference station receiver and user receiver, to determine the processing delay deviation between two receivers.  Identify the accuracy position of a reference station and get the reference station's pseudorange, then derive the difference correction based on the accurate location.  Obtain pseudo-range observations from at least three transmitters.  Calculate the correction of each pseudo-range, then use the least squares solution to find the location. The details are listed in the following part. The measured propagation delays between the transmitting stations and receiver are converted into pseudo-ranges by multiplying the speed. A pseudo-range's equation for referent station is given, based on the delay relationship, which is used to calculate the difference correction: where corresponds to the th transmitter and = 1,2,3. .. .. Other parameters have the same meaning as above. Here, it is necessary that the local clock of the reference station is consistent with

Position Solution Method Integrated Difference
The difference method is usually used to correct common errors between the user and referent station. A scheme of the difference correction system in this paper is shown in Figure 5.
Sensors 2020, 20, x FOR PEER REVIEW 5 of 15 than 30 ns no matter whether it increased or decreased. The amplitude of ( ) is similar to ∆ ( ), whose histogram of each day is plotted in Figure 4. The standard deviations (STDs) of calculation for ( ) are 10.343 ns, 13.521 ns, 11.463 ns and 9.647 ns, respectively. Besides, normal distribution fitting curves with the same STDs are also presented, which are used to explain the characteristic of ( ).

Position Solution Method Integrated Difference
The difference method is usually used to correct common errors between the user and referent station. A scheme of the difference correction system in this paper is shown in Figure 5. Based on the above illustration, the procedure using difference is:  Consistency test of reference station receiver and user receiver, to determine the processing delay deviation between two receivers.  Identify the accuracy position of a reference station and get the reference station's pseudorange, then derive the difference correction based on the accurate location.  Obtain pseudo-range observations from at least three transmitters.  Calculate the correction of each pseudo-range, then use the least squares solution to find the location. The details are listed in the following part. The measured propagation delays between the transmitting stations and receiver are converted into pseudo-ranges by multiplying the speed. A pseudo-range's equation for referent station is given, based on the delay relationship, which is used to calculate the difference correction: where corresponds to the th transmitter and = 1,2,3. .. .. Other parameters have the same meaning as above. Here, it is necessary that the local clock of the reference station is consistent with Based on the above illustration, the procedure using difference is: • Consistency test of reference station receiver and user receiver, to determine the processing delay deviation between two receivers. • Identify the accuracy position of a reference station and get the reference station's pseudo-range, then derive the difference correction based on the accurate location.

•
Obtain pseudo-range observations from at least three transmitters.

•
Calculate the correction of each pseudo-range, then use the least squares solution to find the location.
The details are listed in the following part. The measured propagation delays between the transmitting stations and receiver are converted into pseudo-ranges by multiplying the speed. A pseudo-range's equation for referent station is given, based on the delay relationship, which is used to calculate the difference correction: N i (t)C = n s R di + C(SF i + ASF i + ∆ASF i (t)) + C(δT r + P r ) + Cξ i (t). (5) where i corresponds to the ith transmitter and i = 1, 2, 3 . . . Other parameters have the same meaning as above. Here, it is necessary that the local clock of the reference station is consistent with the time system of the transmitter or that the time deviation is known and all the transmitters in one chain are synchronized or the inherent difference is known in advance. The relative difference correction, removing ξ i (t), is expressed as: In the previous work, the database of differential correction is necessary, which cost a large number of measurements [19]. In this paper, according to the analysis of SF i + ASF i in Section 2.1.1 and the amplitude of ∆ASF i (t) in Section 2.1.2, differential correction for user position is derived depended on distance relation. A simple data map around the referent station is calculated. The pseudo-range integrated difference equation for the user is given as: where C(SF i + ASF i + ∆ASF i (t)) on the left side of Equation (7) is the correction corresponding to different transmitters. It is noted that the environmental features of the path from the transmitter to the reference station are similar to those of the path to the user's receiver. If the paths are very different, the difference method is not suitable.
In the geodetic coordinate system, the Loran system transmitter's position (ϕ i , λ i ) is known. The receiver position (ϕ, λ), R di and δT r + P r are three unknown parameters to be sought. It is necessary to measure three pseudo-ranges used to at least find a position. The algorithm is similar to the GNSS pseudo-range solution based on all available measurements using least squares. Differing from the former algorithm [22], another processing method for the distance between the transmitter and the user is presented in this paper, and details can be seen in Appendix A. With the help of matrix inversion, the main linearized observation equation is thus given by: where A is the design matrix, X is the unknown quantity and B is the observation vector eliminating the difference correction. The least squares solution is derived as: The assumed position P 0 (ϕ 0 , λ 0 , T u0 ) is updated by the increment (∆ϕ, ∆λ, ∆T u ) and the process is repeated. If the new update is less than 1mm, then no further iterations are considered necessary and we arrive at a "best" Loran fix.

Results and Discussion
Raw measurements of propagation delay or pseudo-range contain a number of biases and errors. In order to illustrate the practicability of the pseudo-range method, some detailed simulation tests are discussed in the following parts. The simulation environment is a seawater path, which is partly homogeneous. Another propagation type is considered in the Discussion.

Geometric Dilution of Precision
In the Loran positioning system, the GDOP between the receiver and the transmitters is a dominant deterioration factor of positioning accuracy. The impact of geometry on the accuracy performance of a ranging system is well understood [32]. An effective method for improving position accuracy is to optimize the geometrical placement of transmitters, thereby decreasing the GDOP. On the other hand, another way to improve performance is by increasing number of stations involved in the position solution when the transmitter station's location has been confirmed, However, it is difficult to diminish the GDOP using traditional hyperbolic navigation position fixing due to chain restrictions. The pseudo-range method proposed in this paper is convenient to progress two chains or more, not limiting the number of transmitters. As an example, two chains on the east coast of China are used to analyze the GDOP. The details are listed in Table 1. The GDOP is calculated for three and four sites, respectively, using the geodesic and its azimuth, shown in Figure 6. An area with GDOP ≤ 20 is presented. The left one is for three sites (M, X, Y) and the right one is for four sites (M, X, Y, Z). It is obvious that the coverage area of the three sites is smaller than for the four sites. That is to say, the value of the GDOP for the four sites is smaller than for the three sites on fixed points. Table 2 lists three points within coverage, arbitrarily including their coordinates and corresponding results, which are for later analysis.  The GDOP is calculated for three and four sites, respectively, using the geodesic and its azimuth, shown in Figure 6. An area with GDOP ≤ 20 is presented. The left one is for three sites (M, X, Y) and the right one is for four sites (M, X, Y, Z). It is obvious that the coverage area of the three sites is smaller than for the four sites. That is to say, the value of the GDOP for the four sites is smaller than for the three sites on fixed points. Table 2 lists three points within coverage, arbitrarily including their coordinates and corresponding results, which are for later analysis.

Position Accuracy
When referring to the accuracy of a positioning system, it is necessary to distinguish its absolute accuracy and repeatable accuracy. According the Loran-C User Handbook [33], the absolute accuracy is defined as the accuracy of a position with respect to the geographic or geodetic coordinates of the Earth. The repeatable accuracy, then, is the accuracy with which a user can return to a position whose coordinates have been measured at a previous time with the same navigational system [31,33]. The true test of the position method is the position accuracy between solutions in a receiver and true fix, Figure 6. The GDOP of user based on azimuth, the left one is for three sites and the right one is for four sites.

Position Accuracy
When referring to the accuracy of a positioning system, it is necessary to distinguish its absolute accuracy and repeatable accuracy. According the Loran-C User Handbook [33], the absolute accuracy is defined as the accuracy of a position with respect to the geographic or geodetic coordinates of the Earth. The repeatable accuracy, then, is the accuracy with which a user can return to a position whose coordinates have been measured at a previous time with the same navigational system [31,33]. The true test of the position method is the position accuracy between solutions in a receiver and true fix, where the absolute accuracy and repeatable accuracy are contained.
Before discussion, a specification of some simulation parameters is illustrated in Table 3. Three preferential user sites within coverage range are presented, with coordinates the same as in Table 2. The distances between the test site and each transmitter using Equations (A2) and (A3) in the Appendix A, as well as PF, are also given. Similarly, the predicted value of SF and ASF, assuming all-seawater paths according to Equation (2), are calculated and are presented in Table 3. Based on the former two items, the pseudo-ranges not containing time variations and clock biases are simulated, as seen in Table 3. Additionally, four stations are considered regardless of the number of participants. Here, the different correction progress is ignored and the performance of difference is stressed.

Position without Difference Correction
The method without difference is tested. Based on the observation characteristics of the propagation delay referred to in Section 2.1.2, the pseudo-range consists of the distance between the user and transmitter, SF and ASF. Here, the temporal ∆ASF is ignored due to a small amplitude. Meanwhile, random noise with an STD of 50ns (i.e., about 15m, which is smaller than the typical value) is added to sampling point 1,000 of each pseudo-range. Actually, the range of noise determined by the signal-to-noise ratio (SNR) of the received signal is temporarily not considered here. All terms in the pseudo-range are considered as measurement biases besides distance related to position. Figure 7 gives the scatter position error using the algorithm proposed in Section 2.2. The longitude and latitude error are represented by the horizontal and vertical axis, respectively. The upper three subgraphs, with respect to the user's position A, B, C, are the results of the three sites composed of M, X, Y and the bottom three subgraphs are the results of the four sites composed of M, X, Y, Z. The red star is the statistical average corresponding to the case without random noise. A significant absolute offset of position emerges due to the impact of SF and ASF. It implys that the pseudo-range suffers from large measurement biases resulting in absolute accuracy in the order of hundreds of meters or more, seen in the 95% error radius. The error's directional components, along longitude and latitude in the form of a statistical average, are listed in Table 4, where the STDs of each component represents the repeatable accuracy. Comparing the upper three subgraphs and the bottom one, the repeatable accuracy is improved because the GDOP for the four sites is smaller than for the three sites, that is, a smaller GDOP suppresses the influence of random noise. However, absolute accuracy does not always improve because it also depends on observation error.

Position with Difference
Based on the results of Section 3.2.1, the difference correction is used to mitigate the influence of SF and ASF to obtain the actual distance between transmitter and receiver. For example, the amplitude of SF and ASF is reduced to 70%, that is to say, 30% of the bias in the pseudo-range is corrected. The results are presented in Table 4. Similarly, 70% of biases are corrected and are listed in Table 4, too.
It shows that the latitude and longitude direction's positioning average error decreased proportionally and gradually with the increase in correction. Compared with no difference, the latitude average error is reduced to 873.035 m and 370.236 m, respectively, corresponding to 30% and 70% correction for the A site. However, the STD is almost unchanged because the random noise is not corrected. Ideally, difference would correct all the system errors, implying that there is no SF or ASF in the observations. The absolute error is almost nearly zero for both three and four transmitters, as seen in Table 4. Figure 8 presents the scatter plot of position error by the influence of random noise only, using the algorithm proposed in Section 2.2 with difference. The longitude and latitude errors are represented by the horizontal and vertical axis, respectively. The upper three subgraphs are the results of three sites, composed of M, X, Y. The bottom three subgraphs are the results of four sites, composed of M, X, Y, Z. There is almost no position bias from the statistical mean results, which means a good absolute accuracy performance. The 95% error radius for the B site is smaller than for the A and C sites due to its smaller GDOP. Comparing the upper and bottom subgraphs, the 95% error radius decreases for the bottom case, which also responds to the decreased GDOP for the four transmitters. In a practical application, the difference correction is not always or necessarily proportional, so the correction effect based on arbitrary correction is not improved at all times.

Position with Difference
Based on the results of Section 3.2.1, the difference correction is used to mitigate the influence of and to obtain the actual distance between transmitter and receiver. For example, the amplitude of and is reduced to 70%, that is to say, 30% of the bias in the pseudo-range is corrected. The results are presented in Table 4. Similarly, 70% of biases are corrected and are listed in Table 4, too. It shows that the latitude and longitude direction's positioning average error decreased proportionally and gradually with the increase in correction. Compared with no difference, the latitude average error is reduced to 873.035 m and 370.236 m, respectively, corresponding to 30% and 70% correction for the A site. However, the STD is almost unchanged because the random noise is not corrected. Ideally, difference would correct all the system errors, implying that there is no or in the observations. The absolute error is almost nearly zero for both three and four transmitters, as seen in Table 4. Figure 8 presents the scatter plot of position error by the influence of random noise only, using the algorithm proposed in Section 2.2 with difference. The longitude and latitude errors are represented by the horizontal and vertical axis, respectively. The upper three subgraphs are the results of three sites, composed of M, X, Y. The bottom three subgraphs are the results of four sites, composed of M, X, Y, Z. There is almost no position bias from the statistical mean results, which means a good absolute accuracy performance. The 95% error radius for the B site is smaller than for the A and C sites due to its smaller GDOP. Comparing the upper and bottom subgraphs, the 95% error radius decreases for the bottom case, which also responds to the decreased GDOP for the four transmitters. In a practical application, the difference correction is not always or necessarily proportional, so the correction effect based on arbitrary correction is not improved at all times.

Discussion
In Sections 3.2.1 and 3.2.2 of the paper, the main influence of + and random noise on position accuracy is discussed. The effect of ∆ with a smaller amplitude is omitted because its influence is similar to + , which is system bias corresponding to position. It is shown that difference correction can calibrate the system error only, improving the absolute accuracy. However, if the correction scale of the system error in the observed pseudo-range is inconsistent, i.e., all the

Discussion
In Sections 3.2.1 and 3.2.2 of the paper, the main influence of SF + ASF and random noise on position accuracy is discussed. The effect of ∆ASF with a smaller amplitude is omitted because its influence is similar to SF + ASF, which is system bias corresponding to position. It is shown that difference correction can calibrate the system error only, improving the absolute accuracy. However, if the correction scale of the system error in the observed pseudo-range is inconsistent, i.e., all the observation is non-weighted, then a process method considering unequal precision is very necessary. The least squares solution equation with weight could be expressed as follows: Naturally, it is important to determine the weighting matrix. One method is based on the SNR of each signal, seen reference [19] for details. As for random noise, it is eliminated significantly by filtering or smoothing the pseudo-range. The repeatable accuracy is improved using a pseudo-range after filtering, such as the Kalman filtering progress. When the proposed pseudo-range method is used, the user receiver needed to obtain three transmitter signals at least. The position of the new transmitter and geometric placement must be considered before using this new method in order to ensure signal reception.
In the above discussion, the simulation path is seawater, where SF + ASF is the smallest one of all the ground types shown in Figure 1. Clearly, in another propagation type, such as average land, the absolute position bias will be much larger than that for seawater. The difference method is necessary in order to get a higher absolute accuracy. Compared with the actual difference, the degree of correction depends on the correlation of the two paths. One is the path between the transmitter and the reference station. Another one is the path between the same transmitter and the user receiver. A more complicated situation is that the degree of correction is different for all transmitters. A more detailed analysis should be considered in future work.

Conclusions
In this paper, a theoretical predicted propagation delay of the Loran signal and influential factors are discussed. Moreover, the direct measurement of propagation delay between user and transmitter are illustrated. Both theoretical analysis and measured results showed that propagation delay could be divided into a spatial constant term and temporal term. In particular, we presented a new position solution algorithm combining difference using the Loran system conveniently based on pseudo-range measurement. Some conclusions about this method are: The new algorithm, not limited by the chain, is better than the traditional hyperbolic method. Besides, it is capable of decreasing the GDOP because of greater participation of the transmitters in the position solution. The repeatable accuracy could be improved without changing the amplitude of random noise.
In addition, the absolute accuracy corresponded to the measurement's bias in the pseudo-range, such as SF and ASF. So, in the new method using difference, it is convenient to eliminate the measurement bias of the pseudo-range to improve the absolute position accuracy due to the clear error sources. In detail, absolute positioning offset decreased proportionally with the increase in correction degree about the pseudo-range.
The work in this paper could be applied to the Loran receiver with a filter restraining random noise, thereby improving the repeatable accuracy of the Loran system further. It is able to promote the safety of HEA and the development of the Loran receiver. In the future, the Loran system may create greater value within the navigation field. Different from the former algorithms in [22], another processing method is presented in this paper. The pseudo-range equation is expressed as: The algorithm is based on a GNSS-like pseudo-range solution of all available measurements using least squares. In order to describe it clearly, C(SF i + ASF i + ∆ASF i (t)) is neglected due to its smaller amplitude, however, its influence is reflected in the measurement.
In the geodetic coordinate system, the Loran system transmitter's position (ϕ i , λ i ) is known, and the receiver position (ϕ, λ), R di and δT r + P r are three unknown parameters to be sought. The distance between the transmitter and receiver on an ellipsoid could be calculated by many formulas. A universal formula used in the field of navigation is the Andoyer-Lambert formula [35]: where: S 0i is the arc distance of a sphere, a is the radius of the Earth, δ 0i is the geocentric corner between the transmitter and receiver, ∆S i is a correction due to the sphere to ellipsoidal surface, f = (a − b)/a is the oblateness of the ellipsoid and a and b refer to semi-major and semi-minor axes, respectively, of the WGS-84 ellipsoid.
It is necessary to measure three pseudo-ranges at least used to find a position. So, the process is started by an assumed position P 0 (ϕ 0 , λ 0 , T u0 ) and Equation (A1) had to be linearized using a Taylor expansion: where ρ i = N i C and R di represents the distance between the transmitter and the assumed position. With the help of matrix inversion, the linearized observation equation is thus given by: The most important process is to get the partial derivatives in matrix A. However, ∆δ i is very complex, which is a function of δ 0i , ϕ and λ. The partial differential of ϕ is given by: The right part is expressed: Similarly, the partial differential of ∆δ i about λ is given by: Now, all components in matrix A are listed in Equations (A5) and (A6). Based on the expression of the matrix, the least squares solution is derived as: The assumed position P 0 (ϕ 0 , λ 0 , T u0 ) is updated by the increment (∆ϕ, ∆λ, ∆T u ) and the process is repeated. If the new update is less than 1mm, then no further iterations are considered necessary and arrive at a "best" Loran fix. The whole process of position solution based on pseudo-ranges is clear. The accuracy of the derived position solution can be calculated by comparing the ground truth fixes.