A Shrink-Branch-Bound Algorithm for eLoran Pseudorange Positioning Initialization

: Currently, eLoran is the ideal backup and supplement for global navigation satellite systems. The time synchronization accuracy between stations in the eLoran system has improved, providing conditions for eLoran pseudorange positioning. The pseudorange positioning of eLoran is a nonlinear least-squares problem and the location of the eLoran transmitting stations may cause the above problem to be non-convex. This makes the conventional pseudorange positioning algorithm strongly depend on the initial value when solving the eLoran pseudorange positioning. We propose a shrink-branch-bound (SBB) algorithm to solve the eLoran pseudorange positioning initialization problem. The algorithm ﬁrst uses a shrink method to reduce the search space of the position estimator. Then, optimization is performed using a branch and bound algorithm within the shrunk region, where a trust region reﬂective algorithm is used for the lower bound process. The algorithm can help the receiver to complete the initial positioning without any initial value information. Simulation experiments verify that the algorithm has a success rate of more than 99.5% in solving the initialization problem of eLoran pseudorange positioning, and can be used as an initialization algorithm for pseudorange positioning problems for eLoran or other long-range terrestrial-based radio navigation system.


Introduction
Global navigation satellite system (GNSS) provides all-weather, all-day positioning, navigation, and timing (PNT) services in most outdoor environments. However, in cities or canyons, GNSS performance can degrade due to multipath or poor visibility [1][2][3]. In addition, the high vulnerability of GNSS to interference also seriously affects the security of PNT services [4][5][6]. Many algorithms have been developed to mitigate the performance degradation of GNSS receivers in dynamic multipath environments [7][8][9][10]. However, these algorithms can only improve receiver performance under certain conditions, and it is still difficult for GNSS receivers to work properly in scenarios with fewer visible satellites, such as cities or canyons. Geomagnetic, Wifi, Doppler, and pseudolite-based positioning technologies have been developed for GNSS denial scenarios [11][12][13][14], but these technologies can only provide positioning services in small areas, which cannot meet the positioning requirements of large cities or canyon scenes. In recent years, the eLoran system has regained attention due to its unique system performance, which is expected to solve the existing problems of GNSS [15,16]. The eLoran system is a terrestrial-based radio navigation system that transmits navigation information through a pulse signal with a carrier frequency of 100 kHz. The signal frequency band transmitted by the eLoran system is low and the transmission power is high. Therefore, the eLoran system has the advantages of wide coverage and good anti-interference performance, making it a good backup for GNSS [17][18][19].The traditional Loran navigation system uses a hyperbolic positioning method based on the time difference of arrival (TDOA) [20]. The receiver can only use the stations in a single chain for positioning. Therefore, it has the disadvantage of poor geometric dilution of precision (GDOP), limiting its positioning accuracy. In addition, the TDOA observations include delay errors along the two propagation paths, which makes it difficult to measure and remove abnormal propagation delays. This positioning method cannot directly solve the clock deviation between the receiver and the transmitting station. The eLoran positioning method is based on pseudorange measurement and uses a circular positioning method based on the time of arrival (TOA). This method has the following advantages. First, the receiver uses the signals of multiple chains and multiple stations for positioning, which significantly improves the GDOP factor. Second, the receiver can directly complete the clock error calculation. Third, it can be easily integrated with the wireless positioning system, which helps build an integrated world-ground PNT system [21,22]. Due to limited conditions, the eLoran positioning failed to attract attention in the past. With the transformation and upgrading of eLoran stations, the time between stations in different chains has been synchronized to Universal Time Coordinated (UTC) through technologies such as optical fiber, and the time synchronization accuracy reaches the nanosecond level, providing the basis for the use of eLoran positioning technology. In addition, the application of digital technology in eLoran receivers has improved their sensitivity, which allows them to receive signals from multiple chains and stations simultaneously. Owing to this technical background, the Loran positioning method has regained attention in recent years.
Groves briefly introduced the Loran pseudorange positioning method and pointed out that it was processed by analogy with GNSS-related methods [23]. Yan analyzed the feasibility of Loran pseudorange positioning and the influence of additional secondary factor (ASF) errors on various errors in pseudorange positioning [24]. Kim used the eLoran pseudorange measurements from multiple chains for positioning and performed real-world testing [25]. Peterson and Fang studied the integrated positioning of eLoran and GNSS and pointed out that eLoran pseudorange positioning is a necessary condition for integrated positioning [22,26]. In the above-mentioned literature, eLoran pseudorange positioning is regarded as a nonlinear least squares problem, and local optimization algorithms such as Newton-Raphson algorithm (NR) are used to solve it. However, the eLoran system is not specifically designed for pseudorange positioning, and the location of the eLoran transmitting station may make the problem non-convex. In addition, the nonlinear term in the eLoran pseudorange function is a complex nonlinear function with trigonometric functions, which may cause an ill-condition problem when using the first-order or secondorder derivation information to optimize the objective function. Therefore, for many existing nonlinear least squares algorithms, when the selected initial values are inaccurate, convergence problems to local solutions or erroneous convergence results occur. This initial value dependence affects the ability of the receiver to locate autonomously and causes the receiver to experience localization errors under cold start. At present, there is no literature on the problem of eLoran pseudorange positioning under insufficient initial value information.
This study proposes a shrink-brand-bound (SBB) algorithm to solve the eLoran pseudorange positioning problem. The algorithm first obtains the shrunk region of the estimator through the shrink algorithm. The positioning problem is then solved within this compressed feasible region using a branch-and-bound algorithm, where a trust region reflective algorithm is used for each bound process [27,28]. The SBB algorithm has a global optimization capability and can achieve accurate positioning solutions without initial value information. The algorithm avoids the problem faced by the traditional nonlinear least-squares method by relying on the initial value when solving the Loran pseudor-ange positioning, which further improves the Loran positioning technology based on pseudorange measurement.
The rest of the paper is organized as follows. In Section 2, first, we describe the eLoran pseudorange measurement method and the error in the pseudorange. Then, we build a mathematical model of the eLoran pseudorange positioning and analyze the shortcomings of the NR algorithm in solving it. The principle of the SBB algorithm and the details of each part of the algorithm are introduced. In Section 3, we evaluate the performance of the SBB algorithm and other nonlinear least squares algorithms in solving the eLoran pseudorange positioning problem without initial value information through simulation experiments. Finally, we present the main conclusions of this paper.

Materials and Methods
In this section, we first introduce the pseudorange measurement technology and the error in the pseudorange. Secondly, we construct the mathematical model of eLoran pseudorange positioning and analyze the advantages and disadvantages of the traditional NR algorithm. Finally, we give the principle of the SBB algorithm and the details of each part of the algorithm.

Principle of eLoran's Pseudorange Measurement and Error Analysis
The eLoran positioning technology based on pseudo-range measurement includes two parts: pseudorange measurement technology and positioning algorithm. This section briefly describes the basic principles of pseudorange measurement technology and the error analysis in pseudorange measurement.
The basic principle of eLoran pseudorange measurement is shown in Figure 1. The receiver obtains the signal propagation delay or time of flight (TOF) by measuring the difference between the signal time of arrival (TOA) and the signal time of transmission (TOT). Usually, a certain characteristic point on the eLoran signal is selected as the TOT, such as the initial point or the zero-crossing point in the third circle. The TOA is obtained through the process of a group repetition period, carrier synchronization, and cycle identification. More details can be found in the references [24,29,30].
squares method by relying on the initial value when solving the Loran pseudorange p sitioning, which further improves the Loran positioning technology based on pseudo ange measurement.
The rest of the paper is organized as follows. In Section 2, first, we describe the eLora pseudorange measurement method and the error in the pseudorange. Then, we build mathematical model of the eLoran pseudorange positioning and analyze the shortcom ings of the NR algorithm in solving it. The principle of the SBB algorithm and the deta of each part of the algorithm are introduced. In Section 3, we evaluate the performance the SBB algorithm and other nonlinear least squares algorithms in solving the eLor pseudorange positioning problem without initial value information through simulatio experiments. Finally, we present the main conclusions of this paper.

Materials and Methods
In this section, we first introduce the pseudorange measurement technology and t error in the pseudorange. Secondly, we construct the mathematical model of eLoran pse dorange positioning and analyze the advantages and disadvantages of the traditional N algorithm. Finally, we give the principle of the SBB algorithm and the details of each pa of the algorithm.

Principle of eLoran's Pseudorange Measurement and Error Analysis
The eLoran positioning technology based on pseudo-range measurement includ two parts: pseudorange measurement technology and positioning algorithm. This sectio briefly describes the basic principles of pseudorange measurement technology and t error analysis in pseudorange measurement.
The basic principle of eLoran pseudorange measurement is shown in Figure 1. T receiver obtains the signal propagation delay or time of flight (TOF) by measuring t difference between the signal time of arrival (TOA) and the signal time of transmissio (TOT). Usually, a certain characteristic point on the eLoran signal is selected as the TO such as the initial point or the zero-crossing point in the third circle. The TOA is obtaine through the process of a group repetition period, carrier synchronization, and cycle ide tification. More details can be found in the references [24,29,30]. The eLoran signal is mainly propagated by ground waves and its propagation pr cess is affected by terrain, weather, and other conditions. Interference and noise also affe the measurement during the receiver measurement process, so the , which co tains various additional time delay items, is not the true distance [31,32], as shown Equation (1): , ( The eLoran signal is mainly propagated by ground waves and its propagation process is affected by terrain, weather, and other conditions. Interference and noise also affect the TOA measurement during the receiver measurement process, so the TOF, which contains various additional time delay items, is not the true distance [31,32], as shown in Equation (1): where δt b is the clock deviation between the receiver and the transmitting station, t a is the receiver delay, η(t) is the delay deviation caused by interference and noise in the TOF measurement process, and ∆ASF(t) is the time-related delay due to the ground wave propagation process time-varying factors such as weather. T P is the delay term related to the propagation path as in Equation (2) T P = PF + SF + ASF, where PF is the propagation delay of the signal through the atmosphere and is represented by Equation (3) where c is the speed of light in vacuum, s is the distance between the signal from the transmitter to the receiver; n s is the refractive index of the atmosphere, which represents the ratio of the signal propagation speed in the atmosphere lower than the propagation speed in a vacuum. SF is the propagation delay of the signal through the entire seawater path, which is mainly related to the conductivity of the propagation path. ASF represents the propagation delay of eLoran signal caused by passing through a heterogeneous path of non-full seawater, which is mainly affected by parameters such as distance, surface impedance of the propagation path, and topography. ASF is an important factor affecting the positioning accuracy of eLoran, and it is often calibrated by eLoran differential station or ASF map [33][34][35][36].
In Equation (1), η(t) and ∆ASF(t) are time-related delay items, which are difficult to calibrate. Figure 2 shows the statistical graph of the raw TOF value obtained by the receiver over time. The signal in the picture was transmitted from the Pucheng transmitting station (109.5438 • E, 34.95043 • N) and received in Lintong (109.2221 • E, 34.3686 • N). The fluctuation of the blue line in Figure 1 represents the TOF, which is affected by noise interference and its standard deviation is approximately 9 ns. The red line is the fitted curve of the data shown in blue, representing the fluctuation value with a standard deviation of approximately 10 ns. In order to present these time delays more clearly, we use the Fourier transform to analyze the spectrum of Figure 2a [37], and the obtained spectrum amplitude is shown in Figure 2b. In Figure 2b, we omit the spectrum after 0.001 Hz because its amplitude is too small. Among them, the amplitude at the lowest frequency is about 12 ns, which represents the deviation of the fitted curve in Figure 2a, that is, the delay introduced by ∆ASF(t). Other amplitudes due to measurement noise or interference are around 6 ns. As regards the delay error caused by measurement noise and interference η(t), it is difficult to correct, so we uniformly regard it as noise. The error caused by ∆ASF(t) is often as high as more than 10 ns, so in high-precision eLoran positioning applications, the ASF prediction model is often used for calibration.
The propagation delay error calibration technology is essential for achieving highprecision positioning. There has been considerable research on this aspect [38][39][40]. Now consider the situation after the delay value is calibrated: where TOF c is the calibrated TOF, τ is the time delay value of the signal from the transmitting station to the receiver, δt b is the clock deviation between the receiver and the transmitting station, η(t) is the observation error introduced by the receiver due to timevarying factors such as interference, noise and ∆ASF(t). The t a , SF and ASF in Equation (1) were calibrated. Multiplying both sides by the speed of light is the following pseudorange observation equation: where ρ is the pseudorange observation value of the station received by the receiver, R d is the distance between the transmitting station and the receiver, ρ b is the distance error is shown in Figure 2b. In Figure 2b, we omit the spectrum after 0.001 Hz because its amplitude is too small. Among them, the amplitude at the lowest frequency is about 12 ns, which represents the deviation of the fitted curve in Figure 2a, that is, the delay introduced by . Other amplitudes due to measurement noise or interference are around 6 ns. As regards the delay error caused by measurement noise and interference , it is difficult to correct, so we uniformly regard it as noise. The error caused by is often as high as more than 10 ns, so in high-precision eLoran positioning applications, the ASF prediction model is often used for calibration.  It is worth noting that the eLoran signals mainly propagate through ground waves, and the transmitter and receiver are usually not within the line-of-sight range, so R d cannot be calculated directly using the Euclidean distance formula but needs to be calculated using the great circle distance. The great circle refers to the shortest distance between two points on the surface of a sphere or ellipsoid. The Andoyer-Lambert formula is commonly used in the navigation field to calculate the distance between two points on the earth [41,42]. Suppose the position of the i-th station of eLoran is λ i , ϕ i , and the position of the receiver is (λ, ϕ). Andoyer-Lambert's great circle distance formula is: cos Among them, λ i , ϕ i and λ, ϕ are the longitude and latitude of the transmitting station and the receiver, respectively, and ψ i is the geocentric angle between the i-th eLoran station and the receiver. f and a are the basic geodetic parameters based on WGS-84; the former is the flattening of the ellipsoid, and the latter is the major axis radius of the reference ellipsoid.

eLoran Pseudorange Positioning Model and Conventional Positioning Algorithm
The eLoran pseudorange positioning is solving the estimator x = ϕ λ δt T .
Since the eLoran positioning is a plane positioning system, we only estimate the longitude λ and latitude ϕ. The principle of eLoran pseudorange positioning is shown in Figure 3. Each circle takes the transmitting station as the center and the calibrated pseudorange observation between point A and each transmitting station as the radius. The circles represent all possible solutions to the pseudorange observation of Equation (5). Since x contains three unknowns, the pseudorange observation equations of at least three stations are required to determine x.  When we have no less than three pseudorange observation equations, we obtai by solving the following equation set: The superscript of Equation (7) represents the eLoran station number. Owing to existence of noise in the pseudorange observations, Equation (7) is often transformed i the following least-squares problem: , .
Equation (9) is the basic mathematical model of eLoran pseudorange positioning. NR algorithm is widely used to solve the above problems. The algorithm linearizes Equation (9) through Taylor's formula and transforms it into a linear least-squares pr lem. The basic process is as follows: First, we perform Taylor's first-order expansion of Equation (9) at , and obta When we have no less than three pseudorange observation equations, we obtain x by solving the following equation set: The superscript of Equation (7) represents the eLoran station number. Owing to the existence of noise in the pseudorange observations, Equation (7) is often transformed into the following least-squares problem: Equation (9) is the basic mathematical model of eLoran pseudorange positioning. The NR algorithm is widely used to solve the above problems. The algorithm linearizes Equation (9) through Taylor's formula and transforms it into a linear least-squares problem. The basic process is as follows: First, we perform Taylor's first-order expansion of Equation (9) at x k−1 , and obtain: where Remote Sens. 2022, 14, 1781 Then, using the linear least-squares algorithm, the result is: Finally, the state estimator is: The advantage of this method is that it is simple, and if a suitable initial value x 0 is selected, the convergence speed is fast and the solution is accurate. However, F(x) is affected by the geometry of eLoran stations and may have local minima. Consider a special case, as shown in Figure 4, in which Tr represents the transmitting station, A is the test point, and the four stations are in linear distribution; a common feature as stations are often built along the coastline. It can be seen from the contour line of the function F(λ, ϕ) that there is a local minimum value W in F(x). This means that when using local optimization algorithms such as the NR algorithm [43] or the Levenberg-Marquardt (LM) algorithm [44] to solve the above problem, an inappropriate initial point will cause the algorithm to converge to a local minimum. We will confirm this with a simulation in Section 3. In addition, since the great-circle distance function contained in the eLoran pseudorange equation is a nonlinear term with trigonometric functions, which means that the optimization using the first-order and second-order derivation information of the objective function may face the problem of ill-condition, thereby converging to an erroraneous result. In view of this, it is necessary to design a global optimization algorithm to satisfy the positioning solution in the case of eLoran receiver cold-start.
Then, using the linear least-squares algorithm, the result is: .
Finally, the state estimator is: The advantage of this method is that it is simple, and if a suitable initial value is selected, the convergence speed is fast and the solution is accurate. However, is affected by the geometry of eLoran stations and may have local minima. Consider a special case, as shown in Figure 4, in which Tr represents the transmitting station, A is the test point, and the four stations are in linear distribution; a common feature as stations are often built along the coastline. It can be seen from the contour line of the function , that there is a local minimum value W in . This means that when using local optimization algorithms such as the NR algorithm [43] or the Levenberg-Marquardt (LM) algorithm [44] to solve the above problem, an inappropriate initial point will cause the algorithm to converge to a local minimum. We will confirm this with a simulation in Section 3. In addition, since the great-circle distance function contained in the eLoran pseudorange equation is a nonlinear term with trigonometric functions, which means that the optimization using the first-order and second-order derivation information of the objective function may face the problem of ill-condition, thereby converging to an erroraneous result. In view of this, it is necessary to design a global optimization algorithm to satisfy the positioning solution in the case of eLoran receiver cold-start.

The Shrink-Branch-Bound Algorithm
We define the eLoran positioning solution as the following optimization problem:

The Shrink-Branch-Bound Algorithm
We define the eLoran positioning solution as the following optimization problem: where F : D → R is the objective function, and F is defined in Equation (9). D is the feasible region of x, or search space. λ and ϕ in x have the following constraints The above boundary constraints represent the range of latitude and longitude coordinates of the earth. Since ρ b (δt) and δt have a linear relationship, the selection of the initial value of δt has no effect on the optimization process, so there is no need to consider the range of δt. From now on, we will refer D only to the feasible regions of λ and ϕ.
The SBB algorithm is a modification of the BB algorithm for the eLoran positioning problem. Before introducing the SBB algorithm, the BB algorithm needs to be described first. To solve the problem P, the BB algorithm first obtains a feasible solution as the optimal solution x ∈ D through a certain algorithm, and then iteratively divides the search space D into smaller subsets D s1 , D s2 , . . . , D sn . In each iteration process, when a solution x 1 with a better objective function value can be found in a subset D si , the current solution is updated to x = x 1 , and the subset is divided into smaller subsets; the above process is repeated. If no solution in the subset is better than x , the subset is pruned. When no subset can be pruned, x is the optimal value of P, and the iteration stops. The pseudocode for the generic BB algorithm is given in Algorithm 1 [28,45]. if a solution x 1 ∈ x ∈ D s F(x) < F x can be found, then x = x 1

5.
if D s cannot be pruned: 6.
Remove D s from L 9. Return x The proposed SBB algorithm adds the process of shrinking the feasible region based on the BB algorithm and designs the corresponding branching strategy, bounding method, and pruned strategy according to the eLoran positioning problem. The basic flow chart of the SBB algorithm is shown in Figure 5. We introduce the SBB algorithm from the shrink method and the BB algorithm.

The Shrink Method
From the basic principle of the BB algorithm, the search space D affects the amount of computation of the algorithm. If D can be shrunk, the subsequent BB algorithm can be significantly simplified. The range of D given by Equation (9) is derived from the range of latitude and longitude of the earth. Due to the limited coverage of the eLoran station, we can reduce D according to this feature.
The transmitting power of the eLoran transmitting station is usually fixed, and the eLoran receiver can receive signals from 800 km to 2500 km away from the transmitting station owing to the difference in the propagation path. When the receiver receives signals from multiple stations, it must be within the intersection of the coverage areas of these transmitters. Setting the range of this intersection as , Figure 6 shows the basic schematic for determining . The observable stations are TR1, TR2, TR3, and TR4. The prime vertical arc length between TR2 and TR4 is , which can be estimated by Equation (18); the meridian arc length between TR1 and TR3 is . It is estimated by Equation (19). , Assuming that the maximum working distance between the receiver and the transmitting station is , the range in the blue box of Figure 6 is . Figure 6a,b show under different conditions. The value of can be calculated by Equation (20).

The Shrink Method
From the basic principle of the BB algorithm, the search space D affects the amount of computation of the algorithm. If D can be shrunk, the subsequent BB algorithm can be significantly simplified. The range of D given by Equation (9) is derived from the range of latitude and longitude of the earth. Due to the limited coverage of the eLoran station, we can reduce D according to this feature.
The transmitting power of the eLoran transmitting station is usually fixed, and the eLoran receiver can receive signals from 800 km to 2500 km away from the transmitting station owing to the difference in the propagation path. When the receiver receives signals from multiple stations, it must be within the intersection of the coverage areas of these transmitters. Setting the range of this intersection as D s , Figure 6 shows the basic schematic for determining D s . The observable stations are TR1, TR2, TR3, and TR4. The prime vertical arc length between TR2 and TR4 is W, which can be estimated by Equation (18); the meridian arc length between TR1 and TR3 is L. It is estimated by Equation (19).
where is the equivalent radius of the earth under the WGS-84 model. ̄ is the average latitude of the four stations. In practical applications, the setting of Q does not need to be precise but can be set as the maximum propagation distance according to the receiver performance and actual environment. In addition, Equation (20) is a general equation not limited to the two cases shown in Figure 6a,b. Therefore, once the receiver has identified the station information, Equation (20) can be used to calculate .

The Branch and Bound Method in SBB Algorithm
The proposed branch-and-bound algorithm is as follows: First, a feasible solution ̑ of F on Ds is obtained through a shrink algorithm, and ̑ is assumed to be the global optimal solution. Then, we divide into and and calculate the lower bounds and of function F on feasible domains and . We compare and and retain the subset that has a lower bound , where i=1,2. Thereafter, we compare the order of and . If the order of is smaller than , we update the solution ̑ , and divide again and repeat the above steps. If and are of the same order, or the order of is less than , then ̑ is the global optimal solution. The pseudocode of the SBB algorithm is shown in Algorithm 2. , and , where ∈ 1,2 , 6. If , Then * , the iteration ends; Assuming that the maximum working distance between the receiver and the transmitting station is Q, the range in the blue box of Figure 6 is D s . Figure 6a,b show D s under different conditions. The value of D s can be calculated by Equation (20).
where R e is the equivalent radius of the earth under the WGS-84 model. ϕ is the average latitude of the four stations. In practical applications, the setting of Q does not need to be precise but can be set as the maximum propagation distance according to the receiver performance and actual environment. In addition, Equation (20) is a general equation not limited to the two cases shown in Figure 6a,b. Therefore, once the receiver has identified the station information, Equation (20) can be used to calculate D s .

The Branch and Bound Method in SBB Algorithm
The proposed branch-and-bound algorithm is as follows: First, a feasible solution x of F on D s is obtained through a shrink algorithm, and x is assumed to be the global optimal solution. Then, we divide D s into D s1 and D s2 and calculate the lower bounds F 1 (x 1 ) and F 2 (x 2 ) of function F on feasible domains D s1 and D s1 . We compare F 1 and F 2 and retain the subset D si that has a lower bound F i , where i = 1,2. Thereafter, we compare the order of F s and F i . If the order of F i is smaller than F s , we update the solution x = x 1 , and divide D si again and repeat the above steps. If F s and F i are of the same order, or the order of F s is less than F i , then x is the global optimal solution. The pseudocode of the SBB algorithm is shown in Algorithm 2.
Line 3 of the pseudocode is the branch strategy and we adopt the binary branch scheme as shown in Figure 7. The basic division principle is to make a vertical line at the midpoint of the broadest side of D s to bisect D s . Because the number of local minima on the F function is small, there is no need to divide D s too much, and this binary branch strategy can effectively reduce the amount of calculation without losing the accuracy of the algorithm. Algorithm 2 SBB Algorithm 1. Shrinking D to D s , using Equation (20) 2. Take the initial value x 0 ∈ D, use TRR algorithm to calculate F s ( x ) = {min(F)|x ∈ D s , x 0 } 3. Branch D s into D s1 and D s2 . 4. Calculate F 1 = {min(F)|x ∈ D s1 , x 0 } and F 2 = {min(F)|x ∈ D s2 , x 0 } and their corresponding solutions x 1 and x 2 . 5. F k = min{F 1 , F 2 } and F s , where k ∈ {1, 2}, 6. If F s F k < µ, Then x * = x , the iteration ends; 7. If F s F k > µ, Then F s = F k , D s = D sk , x = x k , and repeat steps 3-5.

Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 20
Line 3 of the pseudocode is the branch strategy and we adopt the binary branch scheme as shown in Figure 7. The basic division principle is to make a vertical line at the midpoint of the broadest side of to bisect . Because the number of local minima on the F function is small, there is no need to divide too much, and this binary branch strategy can effectively reduce the amount of calculation without losing the accuracy of the algorithm. In lines 2 and 4 of the pseudocode, it is necessary to calculate the lower bound of the objective function in the specified feasible region, that is, to solve the following mathematical equation: Ds is determined by Equation (20). Equation (21) is a nonlinear least-squares problem with box constraints, which can be solved by the trust region reflective (TRR) algorithm. Based on the trust region algorithm, the trust region reflective method transforms the boundary-constrained optimization problem into an unconstrained optimization problem through reflection transformation so that each iteration result satisfies the boundary constraints [27]. The TRR algorithm uses the function to fully approximate the behavior of the function in the neighborhood N of , and find the tentative step s in this neighborhood. The pseudocode of the TRR algorithm is shown in Algorithm 3. In lines 3 and 4 of the pseudocode of Algorithm 3, the trust region model to be solved is as follows: where is the gradient of the current , is the Hessian matrix or the approximation of the Hessian matrix of , is the trust region, and ‖ ‖ is the 2-norm. For the solution of Equation (22), please refer to the literature [46,47]. Details of the reflection transformation method in line 5 can be found in the literature [27]. The approximation factor of to in line 6 can be given by Equation (23): , When is greater than the set value , it means that the current approximation effect of to is good and the update step is . Otherwise, the trust region radius needs to be adjusted, the trust region sub-problem solved again, and the above process repeated. In lines 2 and 4 of the pseudocode, it is necessary to calculate the lower bound of the objective function F in the specified feasible region, that is, to solve the following mathematical equation: D s is determined by Equation (20). Equation (21) is a nonlinear least-squares problem with box constraints, which can be solved by the trust region reflective (TRR) algorithm. Based on the trust region algorithm, the trust region reflective method transforms the boundary-constrained optimization problem into an unconstrained optimization problem through reflection transformation so that each iteration result satisfies the boundary constraints [27]. The TRR algorithm uses the function q(s) to fully approximate the behavior of the function F(x) in the neighborhood N of x k , and find the tentative step s in this neighborhood. The pseudocode of the TRR algorithm is shown in Algorithm 3. In lines 3 and 4 of the pseudocode of Algorithm 3, the trust region model to be solved is as follows: where g is the gradient of the current F(x k ), H is the Hessian matrix or the approximation of the Hessian matrix of F(x k ), N is the trust region, and is the 2-norm. For the solution of Equation (22), please refer to the literature [46,47]. Details of the reflection transformation method in line 5 can be found in the literature [27]. The approximation factor ρ k of q(s k ) to F(s k ) in line 6 can be given by Equation (23): When ρ k is greater than the set value µ, it means that the current approximation effect of q(s k ) to F(x k ) is good and the update step is x k+1 = x k + N k . Otherwise, the trust region radius N k needs to be adjusted, the trust region sub-problem solved again, and the above process repeated. Algorithm 3 TRR Algorithm 1. Initial x 0 , N 0 and µ 2. While g(x k ) > µ 3. Build a trust region model q(s) 4. Solve the trust region subproblem, and get s k 5. If s k / ∈ D s 6.
Perform a reflection transform on s k 7. Calculate the approximation ρ k of q(s k ) to F(s k ) and update s k or x k 8. Return x The TRR algorithm can make full use of the feature that the BB algorithm divides the feasible region. When the feasible region is divided, the box constraints will continue to shrink, and the probability of the trust region algorithm converging to the global optimal value will continue to increase. Using the TRR algorithm to obtain the lower bound of F under different feasible regions, the following inequalities must be satisfied. where Lines 5 and 6 of Algorithm 2 are the verification phase. We use µ to verify the convergence process, and µ can be a constant less than 5. When F s /F k < µ, it means that F s and F k are of the same order, and the current iteration value is close to converging to the global optimal value, and the iteration ends. Otherwise, the above branch and bound process needs to be repeated.

Complexity Analysis
The main computational complexity of the proposed SBB algorithm is related to the number of branch iterations N and the convergence accuracy ε. In each iteration, the main computational complexity is related to the update of the bounding process of F(x). More specifically, when we set the norm of the gradient of the solution to be ∇F ≤ ε, the upper bounds of the complexity required to solve steps (2) and (4) are O ε −2 and O 2ε −2 , respectively [48]. Considering the number of branch iterations N, the upper bound of the complexity of the SBB algorithm is O (2N + 1)ε −2 . The upper bounds of the complexity of the following algorithms are shown in the Table 1.

NR [49]
O kn 2 m LM [50] O ε −2 Dogleg/TTR [48] O The above table shows the upper bound of the computational complexity of different algorithms. Among them, the LM algorithm, the Dogleg algorithm, and the TTR algorithm are all Cauchy-related algorithms or Newton-like algorithms, and the upper bound of their complexity is O ε −2 . The complexity of the NR algorithm is related to the number of iterations and the matrix calculation, where k is the number of iterations required, and m and n represent the dimensions of the estimator and the number of equations, respectively. It can be found that the complexity of the SBB algorithm compared with other algorithms mainly lies in N. Since we have shrunk D to D s , this makes the number of branches N usually small, and we will confirm this in simulation experiments.

Results
The SBB algorithm is used to solve the initialization problem of eLoran pseudorange positioning. Therefore, the evaluation of the algorithm is mainly from two aspects. First, the algorithm should still be able to solve the position correctly when no initial value is available, which means that given a random initial value, the algorithm should be able to solve the position accurately. Secondly, the computational complexity of the algorithm should be at a reasonable level so that it can be implemented in the receiver. Based on the above evaluation criteria, this section is organized as follows: we first set the simulation parameters according to the actual station distribution. Then, the performance of various algorithms in solving the eLoran pseudorange positioning problem is compared. Finally, the reliability of the SBB algorithm was verified through simulation.

Simulation Parameter Settings
Assuming that the receiver at point A receives the signals from the four eLoran transmitting stations shown in Table 2, the calibrated pseudorange observations and geodesic distance values between point A and eLoran stations are shown in Table 3, and the atmospheric refractive index n s is 1.000315. Where the calibrated pseudorange observations ρ are as described by Equation (5), they only include the clock deviation δt and the observation error η caused by time-varying delay factor. We set the clock error δt to be 5 µs and η follows a normal distribution, that is, η ∼ N(0, 50).  Figure 8 shows the location of the transmitter station and receiver on the map. To clearly show the influence of (ϕ, λ) on F(x), the contour of F(ϕ, λ) is shown in Figure 9, where δt is set to a known value. The four black contours in Figure 9 are, respectively, surrounded by the solution sets of the four observation equations. The contour shape shows the non-convexity of F(ϕ, λ), which is mainly related to the topology of the transmitting station. Take A as the test point and select the four positions shown in Table 3 as the initial value points. Since δt 0 has no effect on the optimization process, it will always be set to 0 in subsequent simulations.  Since the initialization problem of eLoran has not been studied in the is a lack of competing algorithms for performance comparison. To this end commonly used nonlinear least squares methods, namely, the NR algorit berg-Marquardt (LM) algorithm, and the trust region Dogleg algorithm to the SBB algorithm. The NR algorithm is a commonly used algorithm in po widely used in various pseudorange positioning scenarios. Its advantage culation is simple, and if the initial is suitable, it will converge quickly. C  Since the initialization problem of eLoran has not been studied in the is a lack of competing algorithms for performance comparison. To this end commonly used nonlinear least squares methods, namely, the NR algorit berg-Marquardt (LM) algorithm, and the trust region Dogleg algorithm t the SBB algorithm. The NR algorithm is a commonly used algorithm in po widely used in various pseudorange positioning scenarios. Its advantage culation is simple, and if the initial is suitable, it will converge quickly. Since the initialization problem of eLoran has not been studied in the literature, there is a lack of competing algorithms for performance comparison. To this end, we select four commonly used nonlinear least squares methods, namely, the NR algorithm, the Levenberg-Marquardt (LM) algorithm, and the trust region Dogleg algorithm to compare with the SBB algorithm. The NR algorithm is a commonly used algorithm in positioning and is widely used in various pseudorange positioning scenarios. Its advantage is that the calculation is simple, and if the initial is suitable, it will converge quickly. Currently, only this algorithm is mentioned in the existing papers to solve the eLoran localization problem. The LM algorithm is an algorithm that combines the steepest descent method and Newton's method, and is currently widely used in nonlinear least squares. It is characterized by considering the stability of the steepest descent method and the fast convergence characteristics of Newton's method. This algorithm is a benchmark algorithm for solving nonlinear least squares problems based on the derivation algorithm, and is widely used in various scenarios. The LM algorithm can represent a series of scenarios based on the derivation algorithm to demonstrate the problem of solving the eLoran localization problem based on the first-order derivation and the second-order derivation algorithm.

Analysis and Comparison of Simulation Results
The trust region Dogleg algorithm is representative of another large class of algorithms for solving nonlinear least squares algorithms. It is different from the line search algorithm; the algorithm first sets the step size, and then determines the search direction. The advantage of this algorithm is that it does not require a line search process when solving complex nonlinear least squares problems. Furthermore, even if the condition number of the objective function is poor, it is easy to introduce second-order information of the function. The above three algorithms represent the three most commonly used ideas for solving nonlinear least squares problems. The results are shown in Table 4.  The data in red are the incorrect results, and the data in black are the correct results.
In Table 4, the data in red are the incorrect results, and the data in black are the correct results. The results of all algorithms may be incorrect due to the selection of initial values, except for the SBB algorithm. Among them, both the LM and Dogleg algorithms converge to (31.2167, 103.7164), which is the local minimum L shown in Figure 9. In addition, when the initial value point is close to Point A, both the LM and Dogleg algorithms converge correctly; when the initial value point is close to the local minimum point L, all the results of the above two are incorrect. The erroneous results of the NR algorithm may go beyond the feasible region D, mainly because the convergence of the NR algorithm may be out of control due to the lack of line search. Results from the above table verify that we need a global optimization algorithm to solve the eLoran pseudorange positioning problem when the initial value is not available.
We analyze how the SBB algorithm can always converge to the correct result, regardless of the change in the initial value.
Consider the shrink method of the SBB algorithm. Without loss of generality, we set Q in Equation (20) to 3000, and the reduced feasible region D s is shown as the red box in Figure 10. It can be seen that D s has been significantly reduced compared to D, which reduces the subsequent computation of the SBB algorithm.
To observe the global optimization performance of the SBB algorithm more clearly, Tables 5 and 6 show the iterative process of branch and bound under some initial value points. less of the change in the initial value.
Consider the shrink method of the SBB algorithm. Without loss of generali Q in Equation (20) to 3000, and the reduced feasible region is shown as the r Figure 10. It can be seen that has been significantly reduced compared to reduces the subsequent computation of the SBB algorithm.   It can be seen from Tables 5 and 6 that, as the feasible region is continuously shrunk and divided, the SBB algorithm gradually converges to close to the global minimum.
To further verify the performance of the SBB algorithm, we designed the following simulation experiments: we randomly selected 1000 locations within D s as test points and used the above mentioned algorithms to solve for these locations. Note that these locations were chosen to keep the GDOP as consistent as possible to avoid the impact of GDOP on location accuracy. For each algorithm, the initial value was randomly selected in D and D s . When the positioning error was lower than the set threshold, the solution was successful. The statistical results of the success rate of these algorithms in solving these 1000 positions are shown in Figure 11. simulation experiments: we randomly selected 1000 locations within as test points and used the above mentioned algorithms to solve for these locations. Note that these locations were chosen to keep the GDOP as consistent as possible to avoid the impact of GDOP on location accuracy. For each algorithm, the initial value was randomly selected in and . When the positioning error was lower than the set threshold, the solution was successful. The statistical results of the success rate of these algorithms in solving these 1000 positions are shown in Figure 11.  As shown in Figure 11, the LM and Dogleg algorithms have a success rate of 55% in Figure 11b, while in Figure 11a, the success rates of the two are only 25% and 30%, respectively. This shows that the two algorithms depend strongly on the selection of initial values. The NL algorithm has the lowest success rate, and its solution probabilities are 5% and 35%, respectively, under the two initial value selection schemes. The main reason for the poor performance of the NL algorithm is that it lacks a line search process compared to the LM and Dogleg. The solution success rates of the SBB algorithm under the two initial value selection schemes are 99.9% and 99.5%, respectively, showing good global optimization performance. The possible reason for the failure of the SBB algorithm is that the algorithm will converge to the local minimum value when x 0 is selected very close to the local minimum value. Thus, when x 0 is selected in D, there is a smaller probability of selecting points close to the local minimum. Therefore, the success rate of the algorithm will be improved under x 0 ∈ D compared to under x 0 ∈ D s . To avoid choosing a point near the local minimum as the initial value when using the SBB algorithm, we can choose a point far away from all possible solutions as the initial value point, such as (0, 0).
Computational complexity affects the performance of an algorithm. The previous analysis of the complexity of the SBB algorithm showed that the number of branches, N, has an important impact on the complexity of the SBB algorithm. The figure shows the statistical graph of the number of branch iterations, N, required by the SBB algorithm to complete the positioning solution each time in 1000 positioning simulation experiments. Figure 12 shows that the SBB algorithm needs at most two branch iterations to complete the solution, and even only one branch is required in most cases. Comparing Figure 12a,b, it can be found that the probability that the latter requires two branches to solve is 46%, which is much higher than the 24% of the former. This is because when the initial value is randomly selected in D s , there will be a higher probability of selecting the point close to the local minimum, which makes it converge to the global optimal value after two branches. the solution, and even only one branch is required in most cases. Comparing Figure 12a,b, it can be found that the probability that the latter requires two branches to solve is 46%, which is much higher than the 24% of the former. This is because when the initial value is randomly selected in Ds, there will be a higher probability of selecting the point close to the local minimum, which makes it converge to the global optimal value after two branches.

Discussion
Using pseudorange measurements for positioning in the eLoran system can make full use of the available eLoran stations, thereby expanding the coverage of the eLoran system and improving the positioning accuracy of the system. An important problem with eLoran pseudorange positioning, however, is that the geometric distribution of available eLoran transmitting stations may cause the positioning problem to be non-convex. This makes the existing pseudorange positioning algorithms such as NR algorithms extremely dependent on the selection of initial value. In practical positioning applications, it is difficult for the receiver to obtain reliable initial values in many cases. Therefore, conventional

Discussion
Using pseudorange measurements for positioning in the eLoran system can make full use of the available eLoran stations, thereby expanding the coverage of the eLoran system and improving the positioning accuracy of the system. An important problem with eLoran pseudorange positioning, however, is that the geometric distribution of available eLoran transmitting stations may cause the positioning problem to be non-convex. This makes the existing pseudorange positioning algorithms such as NR algorithms extremely dependent on the selection of initial value. In practical positioning applications, it is difficult for the receiver to obtain reliable initial values in many cases. Therefore, conventional positioning algorithms may converge to wrong solutions due to lack of reliable initial values. At present, there is no literature to study the eLoran pseudorange localization initialization problem.
We transformed the eLoran pseudorange positioning into a nonlinear least squares problem with box constraints and proposed the shrink-branch-bound algorithm (SBB), a global optimization algorithm that can achieve accurate positioning without any initial value. The SBB algorithm first obtains the shrunk region of the estimator through the shrink method. The positioning problem is then solved within this shrunk feasible region using a branch-and-bound algorithm, where a trust region reflective algorithm is used for each bound process. We verified the performance of this method through simulation experiments. The results show that the success rate of the SBB algorithm to solve the position is more than 99.5%, when no initial value is available. However, the success rate of other conventional nonlinear least squares algorithms (such as LM algorithm, Dogleg algorithm) in this case is only around 50%. These results confirm that our proposed SBB algorithm can help the receiver to obtain correct positioning results when no initial value is available.
For the eLoran receiver, both the accuracy of the positioning algorithm and the computational complexity need to be considered. The computational complexity of the SBB algorithm is comparable to traditional Newton-based methods or Cauchy-related methods, which means that it can be implemented in the receiver.

Conclusions
eLoran is the ideal backup and supplement to GNSS systems. The improved accuracy of time synchronization between eLoran stations provides conditions for eLoran pseudorange positioning. We proposed a shrink-branch-bound (SBB) algorithm to solve the eLoran pseudorange positioning problem when the receiver has no initial value available. We verified the performance of the SBB algorithm through simulation experiments. The results show that the success rate of SBB algorithm in converging to the correct result without initial value is over 99.5%, which is more than 40% higher than that of conventional nonlinear least squares algorithms such as LM algorithm and Dogleg algorithm.
The proposed SBB algorithm is expected to make up for the defect that the existing eLoran pseudorange localization algorithm may converge to wrong results when no initial value is available, so it can be used as a cold-start algorithm for eLoran receivers. Therefore, the focus of follow-up research is to combine the SBB algorithm with the existing highprecision positioning algorithms, which is expected to further improve the positioning accuracy and reliability of the eLoran system under high dynamic conditions

Data Availability Statement:
Restrictions apply to the availability of these data. The ownership of data belongs to the National Time Service Center (NTSC). These data can be available from the corresponding author with the permission of NTSC.