**1. 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–3]. In addition, the high vulnerability of GNSS to interference also seriously affects the security of PNT services [4–6]. Many algorithms have been developed to mitigate the performance degradation of GNSS receivers in dynamic multipath environments [7–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–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

**Citation:** Liu, K.; Yuan, J.; Yan, W.; Yang, C.; Guo, W.; Li, S.; Hua, Y. A Shrink-Branch-Bound Algorithm for eLoran Pseudorange Positioning Initialization. *Remote Sens.* **2022**, *14*, 1781. https://doi.org/10.3390/ rs14081781

Academic Editors: Yuwei Chen, Changhui Jiang, Qian Meng, Bing Xu, Wang Gao, Panlong Wu, Lianwu Guan and Zeyu Li

Received: 4 March 2022 Accepted: 6 April 2022 Published: 7 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–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 pseudorange 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.

#### **2. 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.

#### *2.1. 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].

**Figure 1.** Schematic of eLoran pseudorange measurement principle.

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):

$$TOF = TOA - TOT = T\_P + \Delta ASF(t) + \delta t\_b + t\_a + \eta(t),\tag{1}$$

where *δtb* is the clock deviation between the receiver and the transmitting station, *ta* 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. *TP* is the delay term related to the propagation path as in Equation (2)

$$T\_P = PF + SF + ASF\_\prime \tag{2}$$

where *PF* is the propagation delay of the signal through the atmosphere and is represented by Equation (3)

$$PF = \frac{s \cdot n\_s}{c} \,\prime\tag{3}$$

where *c* is the speed of light in vacuum, *s* is the distance between the signal from the transmitter to the receiver; *ns* 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–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–40]. Now consider the situation after the delay value is calibrated:

$$TOF\_{\varepsilon} = \tau + \delta t\_b + \eta(t),\tag{4}$$

where *TOFc* is the calibrated *TOF*, *τ* is the time delay value of the signal from the transmitting station to the receiver, *δtb* 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 *ta*, *SF* and *ASF* in Equation (1) were calibrated. Multiplying both sides by the speed of light is the following pseudorange observation equation:

$$
\rho = R\_d + \rho\_b + \eta\_\prime \tag{5}
$$

where *ρ* is the pseudorange observation value of the station received by the receiver, *Rd* is the distance between the transmitting station and the receiver, *ρ<sup>b</sup>* is the distance error

caused by the clock deviation between the receiver and the transmitting station, and *η* is the distance error representing all other errors that are difficult to calibrate.

**Figure 2.** The schematic diagram of measured propagation delay. (**a**) time delay in time domain (**b**) the amplitude spectrum of delays in the frequency domain.

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 *Rd* 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* , *ϕ<sup>i</sup>* , and the position of the receiver is (*λ*, *ϕ*). Andoyer–Lambert's great circle distance formula is:

$$R\_d^{(i)} = a\psi^{(i)} + \Delta S^{(i)},\tag{6}$$

$$\begin{cases} \cos\psi^{(i)} = \sin\varrho^{(i)}\sin\varrho + \cos\varrho^{(i)}\cos\varrho\cos(\lambda-\lambda^{(i)})\\ \Delta S = a\frac{f}{4}\Big[\frac{\sin\Psi^{(i)} - \Psi^{(i)}}{1+\cos\varrho^{(i)}}\Big(\sin\varrho+\sin\varrho^{(i)}\Big)^2 - \frac{\sin\Psi^{(i)} + \Psi^{(i)}}{1-\cos\varrho^{(i)}}\Big(\sin\varrho-\sin\varrho^{(i)}\Big)^2\Big] \end{cases} \tag{7}$$

Among them, *λ<sup>i</sup>* , *ϕ<sup>i</sup>* and *λ*, *ϕ* are the longitude and latitude of the transmitting station and the receiver, respectively, and *ψ<sup>i</sup>* 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.

#### *2.2. eLoran Pseudorange Positioning Model and Conventional Positioning Algorithm*

The eLoran pseudorange positioning is solving the estimator **x** = *ϕλδt <sup>T</sup>* . 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**.

**Figure 3.** eLoran pseudorange positioning principle.

When we have no less than three pseudorange observation equations, we obtain **x** by solving the following equation set:

$$\begin{cases} \rho^{(1)} - R\_d^{(1)}(\boldsymbol{\varrho}, \lambda) - \rho\_b(\delta t) = 0 \\ \rho^{(2)} - R\_d^{(2)}(\boldsymbol{\varrho}, \lambda) - \rho\_b(\delta t) = 0 \\ \cdots \\ \rho^{(n)} - R\_d^{(n)}(\boldsymbol{\varrho}, \lambda) - \rho\_b(\delta t) = 0 \end{cases} \tag{8}$$

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:

$$\min \{ F(\mathbf{x}) \} = \min \left\{ \sum\_{i=1}^{N} \left[ \rho^{(i)} - R\_d^{(i)}(\boldsymbol{\varphi}, \boldsymbol{\lambda}) - \rho\_b(\delta t) \right]^2 \right\}. \tag{9}$$

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:

$$A \cdot \Delta \mathbf{x} = B\_{\prime} \tag{10}$$

where

$$\mathbf{H} = \begin{bmatrix} \frac{\partial R^1\_{d,k-1}}{\partial \phi} & \frac{\partial R^1\_{d,k-1}}{\partial \lambda} & 1\\ \frac{\partial R^2\_{d,k-1}}{\partial \phi} & \frac{\partial R^2\_{d,k-1}}{\partial \lambda} & 1\\ \vdots & \vdots & \vdots\\ \frac{\partial R^n\_{d,k-1}}{\partial \phi} & \frac{\partial R^n\_{d,k-1}}{\partial \lambda} & 1 \end{bmatrix} \tag{11}$$

$$\Delta \mathbf{x} = \begin{bmatrix} \varrho\_k - \varrho\_{k-1} \\ \lambda\_k - \lambda\_{k-1} \\ \delta t\_k - \delta t\_{k-1} \end{bmatrix} \tag{12}$$

$$\mathbf{B} = \begin{bmatrix} \rho^1 - \begin{pmatrix} R\_{d,k-1}^1 + \rho\_{b,k-1} \\ \rho^2 - \begin{pmatrix} R\_{d,k-1}^2 + \rho\_{b,k-1} \end{pmatrix} \\ \vdots \\ \rho^n - \begin{pmatrix} R\_{d,k-1}^n + \rho\_{b,k-1} \end{pmatrix} \end{bmatrix} . \tag{13}$$

Then, using the linear least-squares algorithm, the result is:

$$\Delta \mathbf{x} = \left(\mathbf{H}^T \mathbf{H}\right)^{-1} \mathbf{H}^T \mathbf{B}. \tag{14}$$

Finally, the state estimator is:

$$
\widehat{\mathbf{x}}\_{\mathbf{k}} = \mathbf{x}\_{\mathbf{k}-1} + \Delta \mathbf{x}.\tag{15}
$$

The advantage of this method is that it is simple, and if a suitable initial value **x**<sup>0</sup> 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.

**Figure 4.** Contour map of *F*(*ϕ*, *λ*) when the transmitting stations are linearly distributed.

#### *2.3. The Shrink-Branch-Bound Algorithm*

We define the eLoran positioning solution as the following optimization problem:

$$\{\mathbf{x}^\*|F(\mathbf{x}^\*) = F\_{\min}(\mathbf{x}), \mathbf{x} \in D\},\tag{16}$$

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

$$\mathbf{x} \in D = \left\{ \begin{array}{c} -\pi \le \lambda \le \pi, \\ -\pi/2 \le \varphi \le \pi/2 \end{array} \right. \tag{17}$$

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 *Ds*1, *Ds*2, ... , *Dsn*. In each iteration process, when a solution **x**<sup>1</sup> with a better objective function value can be found in a subset *Dsi*, 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].


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.

**Figure 5.** A diagram of the shrink-branch-bound algorithm.
