An Improved Coordinate Registration for Over-the-Horizon Radar Using Reference Sources

: In the target localization of skywave over-the-horizon radar (OTHR), the error of the ionospheric parameters is one main error source. To reduce the error of ionospheric parameters, a method using both the information of reference sources (e.g., terrain features, ADS-B) in ground coordinates and the corresponding OTHR measurements is proposed to estimate the ionospheric parameters. Describing the ionospheric electron density proﬁle by the quasi-parabolic model, the estimation of the ionospheric parameters is formulated as an inverse problem, and is solved by a Markov chain Monte Carlo method due to the complicated ray path equations. Simulation results show that, comparing with using the a prior value of the ionospheric parameters, using the estimated ionospheric parameters based on four airliners in OTHR coordinate registration process, the ground range RMSE of interested targets is reduced from 2.86 to 1.13 km and the corresponding improvement ratio is up to 60.39%. This illustrates that the proposed method using reference sources is able to signiﬁcantly improve the accuracy of target localization.


Introduction
Operating in the high frequency band, skywave over-the-horizon radar (OTHR) leverages ionospheric reflection to detect and localize moving targets in the areas beyond horizon and plays a vital role in both military missions and civil applications [1,2]. Coordinate registration (CR), resulted from the propagation of radio signals through ionosphere, is an essential part of OTHR target localization. Roughly speaking, CR performs mapping between slant coordinates and ground coordinates. Specifically, the measurement received by an OTHR which includes slant range and azimuth, is in slant coordinates, while the position of a target in ground coordinates, such as latitude and longitude, is needed in general. Hence, the CR process is mandatory for target localization. The map from ground coordinates to slant coordinates is necessary as well in the target tracking procedure [3]. As a result, CR severely affects target localization accuracy of OTHR.
In the common OTHR systems, the parameters used in CR, are provided by the frequency management system (FMS), which is composed of backscatter sounder and mini-radar systems, omni-directional spectrum monitors and directional antennas, and vertical/oblique incidence sounders and channel scattering function equipment [4]. The ionospheric parameters provided by FMS are substituted into the CR equation and then used for the mapping between slant coordinates and ground coordinates of OTHR.
With different modeling approaches of the ionosphere, the existing OTHR CR approaches can be classified into two categories.
In the first category, a large body of OTHR target tracking literature used simplified OTHR measurement geometry models, where the virtual ionospheric heights (VIHs, often denoted as, e.g., h F2, h F1, h E in the ionospheric research community) are the only used parameter in CR. Assuming the layers of the ionosphere as ideal mirrors, the straight ray paths from the transmitter to the target and from the target to the receiver are assumed to be reflected at the virtual height. According to whether considering the curvature of the earth, the simplified models consist of the planar model and the spherical model. The former one was proposed in [5] and has been widely adopted [3,6,7]. The latter one was proposed in [8], which is more accurate and complex comparing with the former one. The CR errors of the two models caused by the inaccurate VIHs were discussed in [9,10]. Meanwhile, according to the assumptions of VIHs, the existing OTHR target tracking approaches can be classified into four groups. The approaches in the first group assumed that the VIH of each layer is a known constant [11][12][13][14]. The second group, assumed that the VIH of each layer is an unknown constant which could be estimated by a joint optimization scheme [10,15]. The approaches in the third group made the assumption that the VIH of each layer is a random variable or a random process in time [16,17]. Different from the previous three groups, the variation of the ionosphere with location and the spatial correlation of VIHs were considered in [7].
Few works belonging to the second category of OTHR CR approaches have considered full ionospheric models, of which parameters mainly refer to the spatial distribution of electron density. Instead of straight lines, the ray path in the ionosphere is a curve which depends on the operating frequency and the transmit elevation angle of OTHR and electron density. The transmitted signal is continuously bent or refracted as it travels through the ionized plasma in ionosphere, which mainly depends on the electron density and can be mathematically described by ray tracing [18]. In [19], the joint multiple target ground track estimation and data association was considered, where CR is determined by ray tracing with given plasma frequency profile described by the vector consisting of the critical frequencies, heights, and thicknesses of ionospheric layers. In [20,21], the multiquasi-parabolic model was used in CR of OTHR Nostradamus target tracking. In [22], a maximum-likelihood probabilistic multi-hypothesis tracker for OTHR was presented, where CR is performed by ray-tracing with a 3-D ionospheric regional model from the International Reference Ionosphere.
All of the above-mentioned work used the ionospheric parameters provided by FMS as the direct input of OTHR target tracking although the uncertainty of the ionospheric parameters were considered in some of the work. In those work, a problem may arise that the ionospheric parameters provided by FMS are not accurate or consistent with the true ionospheric parameters which are supposed to be used in OTHR target tracking. This is because, on the one hand, the deployment of FMS equipment is restricted to available and limited areas. For example, FMS equipment cannot be deployed in the sea or hostile areas. The number of FMS equipment deployed is limited by financial budget. The ionospheric parameters of unavailable areas are often obtained by interpolation with those of available areas based on statistical ionospheric models. On the other hand, FMS equipment may operate independently with OTHR, resulting in inconsistency on the number of propagation modes and inaccuracy on the parameters for each propagation mode. Therefore, to improve the accuracy of CR, alternative methodologies have to be sought.
One way to tackle the problem above is using reference sources. The reference sources may include HF beacons or transponders [23,24], terrain features (e.g., coastlines, islands and mountains) [25][26][27][28][29], microwave radars [30], and the position information of cooperating targets, such as airliners equipped with automatic dependent surveillancebroadcast (ADS-B) and ships equipped with automatic identification system (AIS). While using the reference sources in CR, their locations in the ground coordinates are assumed to be known. For example, Weijers et al. [24] described the CR experiments using HF beacons around the US east coast, which demonstrated that using beacons can significantly improve CR. A 3-year investigation was conducted using the Wide Aperture Research Facility experimental OTHR to determine if terrain features could be used as CR benchmarks [25], the results of which showed that terrain features could be used in the same way with beacons. In these work, correction factors are obtained by comparing the OTHR measurements of the reference sources with their locations in ground coordinates, which are used to correct the location of the targets of interest surrounding to the reference sources. One advantage of using reference sources is that the received signal of the reference sources is propagated via the same channel with the target of interest. However, this heuristic approach dose not utilize the ionospheric model and signal propagation mechanism, which brings limited improvement to CR.
In this paper, considering the full ionospheric model, we provide a novel method for OTHR CR by inferring the ionospheric parameters using the reference sources and their measurements received by OTHR. The estimated ionospheric parameters are used to the CR of the targets of interest. Describing the ionospheric electron density profile by a quasiparabolic model, we formulate the estimation of ionospheric parameters as a Bayesian inverse problem. Due to the complexity of the OTHR ray path equation, the posterior probability distribution function (PDF) of the ionospheric parameters is approximated by a Markov chain Monte Carlo (MCMC) sampling method, of which convergence and effectiveness are verified by simulation. In summary, the contributions of this paper are as follows.
• For the first time, we propose to use reference sources and their corresponding OTHR measurements to infer the ionospheric parameters using the full ionospheric model, potentially providing improvement to CR. • Considering the complexity of the OTHR ray path equation, an MCMC-based approach is implemented to estimate the ionospheric parameters. The estimated ionospheric parameters are then used to the CR of the target of interest, leading to a significant improvement.
The remainder of this paper is organized as follows. Section 2 consists of assumptions, models, and problem formulation of CR for OTHR target localization. The proposed CR method is depicted in Section 3. Simulation results are presented in Section 4 followed by concluding remarks in Section 5.

Models and Problem Formulation
In this section, we describe the models we introduce for ionospheric electron density profile, ray path equations, reference sources, and OTHR measurement. Based on these models, we present a problem formulation of ionospheric parameters inference as a Bayesian inverse problem. Then the CR of the targets of interest is formulated as a forward model evaluation given the estimated ionospheric parameters.

The Ionospheric Electron Density Profile Model
In reality, as we mentioned in Section 1, the OTHR signal is continuously bent or refracted as it travels throughout the ionized plasma in ionosphere, which is mainly influenced by the electron density and could be mathematically described by ray tracing technique [18]. Regarding description of the ionospheric environment through which the signal is propagating, there are two categories of models, ionospheric profile models and 3D synthetic models. Ionospheric profile models use explicit mathematical expression to describe variation of the ionospheric electron density along the altitude, including liner model, exponential model, parabolic model and quasi-parabolic model, etc. There are a few parameters in those models which make them easy to be used in OTHR, while with limited accuracy. 3D synthetic models, most of which are global empirical models, are largely based on data, and use appropriate mathematical functions to represent the characteristic variation patterns seen in long data records [31]. The representatives of these models are Bent model, NeQuick model, International Reference Ionosphere model, etc.
The performance of these models depends the accumulation of the ionospheric data; they may not work well in regions with limited equipment.
For the convenience of derivation and expression, we adopt the widely used quasiparabolic model [32] that allows the exact derivation of ray path equation for the ionospheric electron density profile. Although the quasi-parabolic model was first proposed in 1968 [32], it is still widely adopted in the literature [33][34][35]. In the quasi-parabolic model, given a location (latitude and longitude), the electron density in a layer of the ionosphere is defined as: where N e (r) is the electron density at a distance r from the center of the earth, r b is the value of r at the bottom of the layer, y m is the semi-thickness of the layer and r m = r b + y m , N m is the maximum electron density, i.e., the electron density at r m . Correspondingly, if the thermal motion of the electrons is ignored, the critical frequency of the layer is calculated as f c 2 = N m e 2 /(4π 2 0 m * ) [36], where e, 0 and m * are electron charge, permittivity of free space and electron mass, respectively. By the fact that the electron charge e = 1.602 × 10 −19 C, the electron mass m * = 9.1094 × 10 −31 kg and the permittivity of free space 0 = 8.8542 × 10 −12 A 2 s 4 m −3 kg −1 , one can get f 2 c ≈ 80.6N m . From Equation (1), it is seen that, given the quasi-parabolic model and the ionospheric parameter vector ζ = [r b , f c , y m ], N e (r) is uniquely determined.
In this paper, we consider using the OTHR measurements of reference sources via single-mode propagation, i.e., the propagation from the transmitter to a reference source and the propagation from the backscattered signal to the receiver are both supported by a single ionospheric layer (e.g., E-E mode, F-F mode). By using the information of reference sources in ground coordinates and the corresponding OTHR measurements, we aim to estimate the ionospheric parameter vector ζ = [r b , f c , y m ] of each layer. Note that as we mentioned above, more advanced and precise models on ionospheric electron density profile are available in the literature. The method in this paper maybe extended to those models that may have more parameters. For example, if the electron density model consisting of a Chapman layer with tilts, rippels and gradients is applied, the parameters to be estimated may include critical frequency at the equator, height of the maximum electron density at the equator, amplitude of periodic variation of squared critical frequency with latitude, period of variation of squared critical frequency with latitude, coefficient of linear variation of squared critical frequency with latitude, and tilt of the layer, etc, [18]. If the CCIR-67 model is adopted, the to-be-estimated parameters are the coefficients of the model [37]. Aiming to provide real-time ionospheric characteristics, data assimilation using observations from vertical incidence ionosondes, oblique ionosondes [38], etc, can be exploited before applying our method. We leave the extensions to the utilization of the OTHR measurements of reference sources via mixed-mode propagation (e.g., E-F mode, F-E mode) and more complicated models of electron density profile as our future work.

Ray Path Equations
Given the above quasi-parabolic model, ray path equations can be derived. Figure 1 illustrates the geometry of a ray trajectory in the quasi-parabolic layer, where the notations are defined as follows.
• β: ray elevation angle, i.e., the angle between ray path and the line that is parallel to the horizon; • β 0 : initial ray elevation angle, i.e., the value of β when r = r 0 ; • r t : r at the top of the ray path; Center of Earth Figure 1. The geometry of the ray trajectory [32].
For ease of exposition, define the ratio of the OTHR operating frequency f w to the critical frequency of the layer f c as f = f w / f c . Accordingly, the refractive index of the layer at r is defined as [32]: By Snell's law [32], for the ray in Figure 1, the product rµ cos β is constant at any position in the layer, i.e., rµ cos β = r 0 cos β 0 . Therefore, considering the geometry of the ray segments, the reflection by the ionosphere and the symmetry of the ray, the ground distance R D and the group path distance R P that the ray traverses are calculated as [32]: where R D is the distance measured along the earth's surface, and R P is the distance that calculated as signal transit time multiplied by the light speed c, and s t is real ray path length at r t . Given the definitions of all notations in the above equations, and based on standard integration table [39], the integrations in Equations (3) and (4) are calculated as: where γ = cos −1 (r 0 /r b cos β 0 ) by the law of sines. The variables A, B and C are written as: According to the above geometry, at the ray apogee, one has: The maximum of the ray height r t,max = −B/(2A). For ray tracing using other ionospheric electron density profile models, the reader is referred to [18,40,41].

Reference Sources and Targets Modeling
As we mentioned in Section 1, the reference sources may include HF beacons or transponders, terrain features, microwave radars, and airliners information from ADS-B and ships information from AIS.
Assume that there are L reference sources in the OTHR surveillance area. The ground range of reference source l (l = 1, . . . , L) is denoted as ρ l . The measurement of the ground range of the reference source l can be written as: where u l is the Gaussian measurement noise with variance U l . The variance is different for different type of reference sources. For example, the measurement noise of ADS-B is small comparing with the measurement noise of OTHR. Therefore, it can be ignored. So are the HF beacons and terrain features. The measurement set of all reference sources is written as: Denote the ground range of the targets of interest in the OTHR surveillance area bȳ ρ n , where n = 1, . . . , N, and N is the number of targets. Note that the CR of each target is carried out separately. Some reference sources, such as airliners and ships could be targets of interest as well.

OTHR Measurement Modeling
Based on the above models, the OTHR measurement models of targets and reference sources are given here. We focus on the range accuracy of targets in CR; hence the related measurement is the group path distance in Equation (6). For each reference source or target, the measurement equation is written as, where v is Gaussian white noise with variance V.
The measurement sets of all reference sources and interested targets are written as y = {y 1 , . . . , y L } andȳ = {ȳ 1 , . . . ,ȳ N }, respectively. Meanwhile, the initial elevation angle sets of the rays that reach each reference source and interested target are written as It is worth noting that, only few OTHR, such as French Nostradamus [21] can measure the elevation angle. Generally, the initial elevation angle β 0 is unknown as we assume in this paper. Meanwhile, although the relationship between OTHR measurement y and reference source or target ground range ρ in Equation (11) is not explicitly shown, the initial elevation angle β 0 determines whether the signal can be reflected back to the receiver via the reference source or the target, and OTHR measurement y is implicitly related to the reference source or the target ground range ρ by β 0 . In other words, the ground distance R D of the ray that reaches the reference source or the target is equal to the ground range of the reference source or the target, i.e., where φ is defined as in Equation (5).

Problem Formulation
Our aim is to improve the accuracy of CR by inferring the ionospheric profile model parameters vector ζ from reference source data z and their OTHR measurements y and then performing the CR of the interested targets based on the inferred ionospheric parameters ζ and the measurements of the targetsȳ. Note that an important preliminary procedure to apply our method is to associate the reference sources with their corresponding OTHR measurements. This procedure is out of the scope of this paper. The reader is referred to [24][25][26]. In other words, we assume that the pair (y l , z l ) of a reference source and its corresponding OTHR measurement is known.
Based on the above-mentioned models, the block diagram of our problem, which consists of a typical inverse problem [42] and a forward model evaluation problem, is shown in Figure 2. By Bayes' rule, given the OTHR measurements y and reference source data z, the joint posterior PDF of ionospheric parameters vector ζ and the initial elevation angle vector β 0 is written as:

OTHR transmitter
where π 0 (ζ, β 0 ) is the joint prior distribution of ζ and β 0 . The analytical solution of Equation (13) can not directly obtained because of the implicit OTHR measurement function.
As will be shown in the next section, we resort to a sampling-based method.
Once the estimation ofζ is derived, the CR of the target of interest is performed, which includes two steps. The first step is to find the elevation angleβ 0 with OTHR measurement y andζ through Equation (6) based on a dichotomy method. The second step is to calculate the ground distance with elevation angleβ 0 andζ through Equation (5).

Solutions
In this section, the estimation of the ionospheric parameters and the CR process are introduced in detail. Firstly, the joint posterior PDF of ionospheric parameters ζ and the initial elevation angle β 0 is derived. Then, MCMC-based estimation of ζ is proposed, and the ESAI sampling method and its convergence evaluation are followed.

Derivation and Framework
Given ζ and β 0 of all rays that reach the reference sources, and assuming that the OTHR measurements y are independent over different reference sources, the conditional PDF (measurement likelihood) of the OTHR measurements y is written as: Similarly, assuming reference sources data z are independent, then the conditional PDF (measurement likelihood) of z is written as: As the measurement noises v l in Equation (11) and u l in Equation (10) are both zero mean and Gaussian, the measurement likelihoods are written as: Apparently, ζ and β 0 are independent, i.e., π 0 (ζ, β 0 ) = π 0 (ζ)π 0 (β 0 ), where π 0 (ζ) is assumed to be Gaussian with meanm 0 and covariance matrixΣ 0 . Initial elevation angles are independent as well. We assume that each initial elevation angle follows a continuous uniform distribution over a certain interval. That is, π 0 (β 0 ) = ∏ l=1,...,L π 0 (β l 0 ), and π 0 (β l 0 ) = U [λ 1 , λ 2 ]. As the denominator in Equation (13) is a constant, Equation (13) is simplified as: which is then expressed as: where the operator D(x, Σ) denotes the operation D(x, Σ) = −x T Σ −1 x. Solving Equation (19) directly is intractable. This is because, firstly, from Equations (5), (6), (11) and (12), no explicit expression exists between the function ψ(·) and the function φ(·), (i.e., the reference source ground distance ρ and the OTHR measurement y). Secondly, β 0 is a hidden variable. The fact that β 0 and ζ are coupled in the OTHR measurement process makes their relationship be expressed implicitly by equation ρ l = φ(ζ, β l 0 ). As the elevation angle vector β 0 is also unknown in the ionospheric parameters estimation, there are L + 3 variables to be estimated. If the measurement noise u l of reference sources (e.g., airliners with ADS-B, ships with AIS, HF beacons, and terrain features) in the ground coordinates is small, we let z l ≈ ρ l to reduce the computational cost. In this case, given ζ, β l 0 can be found based on z l through a dichotomy method.
Accordingly, the dimension of the space of the estimated parameters is reduced from L + 3 to 3. Equation (19) is then simplified as: If the measurement noise of the reference sources in the ground coordinates is too big to be ignored, which means the approximation z l ≈ ρ l cannot be applied, then simultaneous estimation of both the ionospheric parameters and the initial elevation angles is needed. We leave it as our future work.
Because of the complexity of the ray path equations mentioned above, a samplingbased method is adopted to solve Equation (21), which generates a set of samples to approximate the posterior PDF (21), and the estimation of ζ can be obtained. More specifically, Markov Chain Monte Carlo (MCMC) [43], which is a powerful algorithm framework for handling complicated inference problems, is used in this paper to estimate ζ.
MCMC constructs a Markov chain (an ensemble of chains) that has (have) the desired distribution (e.g., the posterior PDF (21)) as its (their) equilibrium distribution. Then the samples of the desired distribution can be obtained by recording states from the chain. The simplest and most popular MCMC method is Metropolis-Hastings (MH) algorithm, of which the basic idea is as follows. Given the distribution h(ζ), a proposal distribution q(ζ |ζ) and an initial value of ζ, an MH step involves sampling a candidate value ζ according to q(ζ |ζ) [44]. The Markov chain then moves towards ζ with the acceptance probability a(ζ |ζ) = min{1, [h(ζ)q(ζ |ζ)] −1 h(ζ )q(ζ|ζ )}; otherwise, it remains at ζ. This process is repeated until enough samples are generated.
Denote the obtained samples by {s i } N s i=1 . The minimum variance estimation of ζ is: To use the samples of the Markov chain when it reaches its equilibrium distribution, which brings more accurate estimation, burn-in strategy is often used when exploiting MCMC-based method. That is, the samples generated at the beginning of an MCMC run, e.g., the first half, are throw away [45]. Additionally, the thinning technique [46,47] is used in order to reduce correlation between samples. That is, instead of keeping all samples, only every ξ th (e.g., ξ = 10) iteration is kept. At last, perform convergence test of the chain. If the chain converges, the generated samples are used to estimate the ionospheric parameters by Equation (22); otherwise, find alternative samplers. We summarize the complete procedure for the estimation of ionospheric parameters ζ as below. (1) Initialization. Using the prior distribution of the ionospheric parameters to initialize the sampler.
Sampling. Draw samples from the posterior PDF in Equation (21). The detailed sampling method will be descried in the following Section 3.2.
(3) Burn-in. The beginning part, e.g., the first half, of the samples are throw away.
(4) Thinning. Only every ξ th (e.g., ξ = 10) sample is kept. (5) Convergence evaluation. Determine if the Markov chain converges to the distribution (21) based on the method will be introduced in the following Section 3.3. If the Markov chain converges, go to next step; otherwise, find alternative samplers and go to Step 2.
Note that no MCMC method exists which guarantees that the Markov chain converges to the desired distribution [48]. In Step 5 of the above procedure, if the number of samples exceeds the given maximum number of samples and the convergence of the Markov chain is not reached, one may need to adjust the parameters of the sampling method or resort to other sampling methods.
Note that the posterior PDF (21) is implicitly defined. We found that the standard MH method might need a long chain to sample from Equation (21), which is very timeconsuming. Alternatively, an improved MCMC algorithm, namely ensemble samplers with affine invariance (ESAI) [49], is adopted to estimate ionospheric parameters and will be detailed in the next section.

ESAI Sampling
ESAI is particularly useful for sampling badly scaled distribution and can produce independent samples with a much shorter autocorrelation time than the standard MH method [49]. The details of ESAI sampling for the estimation of the ionospheric parameters are as follows. Let S = {s i (t)} i=1,...,I;t=1,...,T e represent I chains which are also called stochastic processes of "walkers" that move around randomly in the space of ζ, where T e is the length of each walker. For each walker i, its state is initialized as s i (1) =m 0 , which is the mean of the prior distribution of ζ. For the choice of the number and the length of the walkers, the read is referred to [50].
Those walkers sample separately to generate samples from the distribution (21), which are denoted as the states of those walkers. The update of each walkers is according to the following rules [49].
The new state of i-th walker depends on the current states of the other I − 1 walkers, i.e., the complementary set S [i] = {s j (t), ∀j = i}. To update the state of the i-th walker s i (t), a walker s j (t) is randomly and uniformly chosen from the complementary set S [i] , and then the new candidate state of s i is computed as: where d is a random variable sampled from the distribution g(d). A feasible distribution of g(d) is [49]: where the parameter τ is adjustable to improve performance. If the new state of the walker s * i is accepted with probability: where M = 3 is the dimension of the parameter space and h(·) is the distribution in Equation (21), set s i (t + 1) = s * i ; otherwise, set s i (t + 1) = s i (t). The detailed single update step at t of all walkers is summarized in Algorithm 1. By repeating the update step in Algorithm 1, we can obtain the samples drawn from Equation (21). (21) 1: for i = 1, . . . , I do 2: Randomly select a walker s j (t) with equal probability from the complementary set S [i] ; 3: Sample d from the distribution g(d) in Equation (24) (23); 5: Compute β * 0 = [β 1, * 0 , . . . , β L, * 0 ] by setting ζ = s * i and Equation (20); 6: Compute

Algorithm 1 A single update step of all walkers to sample from Equation
Sample r from U [0, 1]; 8: if r ≤ a then 9: 10: else 11: 12: end if 13: end for

Convergence Evaluation
It is necessary to determine the convergence and the length of the Markov chain when using MCMC methods, for which various methods are available. For a recent review, the reader is referred to [51]. In this paper, the Gelman-Rubin (GR) diagnostic [52] and the graphical method [51] are used to heuristically evaluate the convergence of the sampler.
For ESAI, the GR diagnostic relies on multiple ensembles {S 1 , S 2 , . . . , S m }, each of which starts at different initial points and contains n samples. It is worth noting that, each ensemble S i is independently sampled using ESAI which includes I chains. One should not compute GR statistic using multiple chains in the same ensemble because the chains are not independent. Using those ensembles, GR calculates two estimators of the variance of S, namely within-ensemble variance estimate [52]: where S ij is the jth samples in ensemble S i , and the pooled variance estimate [52] where is the between-ensemble variance estimate, andS i andS are the ith ensemble mean and the overall mean, respectively. Then the multivariate potential scale reduction factor (MPSRF) is used for convergence diagnostic, which is denoted byR and given as [52]: where λ 1 is the largest eigenvalue of the matrix (W −1 B)/n. Convergence is reached when R sufficiently approaches to one. The trace plot, which is a time series plot that shows the realization of the Markov chain at each iteration against the iteration numbers, is used to visualize how the Markov chain is moving around the parameter space. That is, how well the Markov chain is mixed which means how well it is close to its equilibrium distribution. The chain is likely to be converged only if the random-walk process crosses the domain of ζ many times during the MCMC run [48].
In addition, the autocorrelation time is also adopted to evaluate the convergence of the Markov chain. The lag-t autocorrelation is defined to be the correlation between the samples t steps apart. The autocorrelation plot shows the values of the lag-t autocorrelation function (ACF) C e (t) against increasing t values. For fast-mixing Markov chains, lag-t autocorrelation values drop down to zero quickly as t increases. The autocorrelation function C e (t) is defined as [50]: In practice, one can estimate C e (t) for a Markov chain of N s samples as: whereζ is the mean of all samples and calculated in Equation (22).

Simulation and Analysis
In this section, numerical simulations are performed to verify the performance of proposed OTHR CR method based on MCMC. The method is implemented in Matlab R2018b and ran on an Inter (R) Core (TM) i5-4590 CPU@3.3GHz computer.

Scenario Settings
The detailed parameter settings are given as follows. To show the performance of our proposed method, we generate 30-scan data, including the OTHR measurements and position of reference sources and interested targets. The sampling period T s = 20 s. It is assumed that the earth is a regular sphere, i.e., the radius of the earth r 0 = 6371 km is the same everywhere. The measurement noise standard deviation of the slant range The initial states of the airliners and the targets are found in Table 1. The measurements y of 20 airliners received by OTHR at the first scan are shown in Table 2. Table 1. Initial states of 20 airliners and 20 interested targets. The first two lines are the ground range ρ l (km) and ground range rateρ l (km/s) of each airliners. The last two lines are the ground range ρ l (km) and ground range rateρ l (km/s) of each interested targets. In order to visualize the posterior PDF of the ionospheric parameters in Equation (21), we show it in a selected cuboid parameter space with the boundary m 0 ± [8 km, 0.6 MHz, 12 km]. The selected parameter space is discretized into 40 3 grids in the selected domain. The values of Equation (21) at these grid points calculated using the OTHR measurements of the first four airliners in Table 2 are shown in Figure 3. In Figure 3, the first three subgraphs are the planes where y m = 20.9 km, f c = 3.69 MHz, and r b − r 0 = 79.8 km, respectively; the last subgraph is a corner of the selected parameter space in a 3D view.

Convergence Evaluation of ESAI and Ionospheric Parameters Estimation
Before giving the convergence analysis of ESAI, an intuitive illustration of samples generated by ESAI is given. For ESAI, we use 100 walkers. 80,000 samples are sampled and only the latter half are taken due to the burn-in. At the same time, the thinning technique is used in order to reduce the correlation between samples. Finally, 4000 samples are obtained. The samples of the posterior PDF of the ionospheric parameters using the OTHR measurements of the first four airliners shown in Table 2 are shown in Figure 4. The samples are more dense in Figure 4 in the sub parameter space that the posterior PDF is large in Figure 3. Figures 3d and 4d show the samples that correspond to the posterior PDF, which has intuitively verified the effectiveness of ESAI.  respectively. The GR MPSRF is calculated base on these ensembles. The plots of iterativê R at the increment of every 100 iterations are given in Figure 5. It is seen that MPSRF approaches to one rapidly, which means that the used ESAI converges. The convergence analysis of the samples in terms of the graphical method and the autocorrelation time is shown in Figure 6. Figure 6a shows the trace plots of the first 1000 samples and Figure 6b shows the autocorrelation plots of the Markov chain samples. From the track plots, it is seen that the chains are mixed well, i.e., the chains move around the parameter space many times and are not stuck in any region. At the same time, the autocorrelation of the samples is almost negligible after lag 15, which means the chains mix fast. In conclusion, the used ESAI sampler is effective and the parameter setting of the sampler is appropriate.  In order to evaluate the performance of the ESAI-based approach on ionospheric parameters estimation, different numbers of airliners are used. The ionospheric parameters estimation results are shown in Figure 7. Table 3 shows the mean values of the ionospheric parameter estimation in each dimension in 30 scans under different airliner numbers in Figure 7, and the improvement ratio compared with the a prior values. By observing Figure 7 and Table 3, it is seen that, the RMSEs of the ionospheric parameters r b , f c and y m are reduced to 1.73 km, 183.3 kHz and 3.61 km, respectively, when using the ADS-B data of 20 airlines, and the corresponding improvement ratio compared with the a prior values are 65.89%, 8.8%, and 10.74%, respectively. Meanwhile, as we expected, as the number of airliners increases, the estimation error of the ionospheric parameters decreases. Overall, the accuracy of the ionospheric parameters are effectively improved when utilizing the proposed method with reference sources. Table 3. Statistics on estimated ionospheric parameters in Figure 7. IR is the abbreviation of improvement ratio.

Analysis of OTHR CR Results
Given the effectiveness of ESAI in the last section, this section conducts statistical analysis of the results of the OTHR CR scheme proposed in Section 3. With OTHR measurements, the following three cases are considered for comparison.

•
Case 1: The ground ranges of targets are calculated using the a prior of the ionospheric parameters; • Case 2: The ground ranges of targets are calculated using the true ionospheric parameters; • Case 3: The ground ranges of targets are calculated using the estimation of ionospheric parameters with ADS-B data of different number of airliners (i.e., 4, 8, 12, 16, 20 airliners).
The results of the three cases are compared with the true ground ranges of targets to analyze the improvement of the proposed OTHR CR method. Figure 8 shows the RMSE comparison of the target ground ranges during OTHR CR under the three cases. Table 4 shows the RMSE mean values(m) of each curves over 30 scans in Figure 8, and the improvement ratios of Case 3 with respect to Case 1. From Figure 8 and Table 4, it is seen that, comparing with using the a prior ionospheric parameters, the RMSE of ground range is reduced from 2.86 to 1.04 km when using ionospheric parameters which are estimated with ADS-B data of 20 airliners and the corresponding improvement ratio is 63.66%. Overall, the more accurate the ionospheric parameters, the smaller the RMSE of the ground range during OTHR CR. The results show that the error of OTHR CR is reduced when using reference sources. The effectiveness of the proposed OTHR CR method is verified.
We emphasize that this improvement ratio may be relevant to the model of ionospheric electron density profile and is related to how close are the reference sources to the interested targets. In this paper, we use the quasi-parabolic model. As we mentioned in Section 2.1, the method can be applied to other models. Any model is a simplified and an approximate representation of the real and complicated ionosphere using mathematical expressions with parameters. The difference when applying different models of ionospheric electron density profile to our problem is that the to-be-estimated parameters are different. In general, the more complex the model is, the more accurate the representation, the more parameters we need to estimate, and the more computationally expensive the estimation of the parameters. We leave the evaluation of our proposed method when exploiting different models of ionospheric electron density profile as our future work. In our simulation, we have assumed that the reference sources and the interested targets are distributed in a relatively small area such that their corresponding signals are reflected via a neighbouring small ionospheric area. Intuitively, the closer the reference sources to the interested targets, the more accurate of the parameter estimation is. In practice, reference sources may be distributed in a large area or even in the entire OTHR coverage area. In this case, all the information of the reference sources and the corresponding OTHR measurements are collectively used to estimate the parameters together with the model of the ionospheric electron density profile, which potentially provides improvement on the estimation of all the parameters globally.

Conclusions
We proposed a novel method for OTHR CR through using MCMC to estimate the ionospheric parameters with reference sources. The quasi-parabolic model was used to describe ionospheric electron density profile. The estimation of the ionospheric parameters was formulated as a Bayesian inverse problem, which was then solved by an MCMC method, namely ESAI, due to the complicated ray path equation. The entire implementation procedure of ESAI for the estimation problem was described, followed by a description of convergence evaluation of ESAI using two different methods. Then the OTHR CR was performed based on the estimated ionospheric parameters. The simulation results have shown that, the estimated ionospheric parameters are more accurate comparing with the a prior value of ionospheric parameters, which in turn greatly reduced the OTHR CR error of interested targets. For future work, we plan to apply the proposed method to more models of ionospheric electron density profile and evaluate the corresponding improvement on OTHR CR. We will also consider using the OTHR measurements of reference sources via mixed-mode propagation to further reduce CR error.