A Single-Step Method for Over-the-Horizon Geolocation Using Importance Sampling

0is paper investigates the geolocation for an over-the-horizon (OTH) transmitter observed by widely separated arrays. We propose a maximum likelihood (ML) based direct position determination (DPD) method to directly locate the transmitter in a single step by exploiting the position information embedded in azimuth angles. 0e Monte Carlo importance sampling (IS) technique is employed to find an approximate global solution to this DPD problem, where the importance function analogous to Gaussian distribution is derived.0is enables the transmitter to be precisely located with low complexity in a noniterative manner. Additionally, we derive the Cramér–Rao bound (CRB) expression for the investigated problem. 0e simulation results corroborate the superior localization performance of the proposed method with respect to the conventional two-step approaches and the iterative DPD method.


Introduction
Passive transmitter localization with antenna arrays is a significant task in various fields such as signal processing, wireless communication, navigation, and radio astronomy [1]. Over-the-horizon (OTH) geolocation has many applications, such as navigation and the long-range reconnaissance and surveillance [2]. Here, we consider the scenario of OTH geolocation of a high-frequency (HF: 2 MHz to 30 MHz) transmitter observed by widely separated arrays.
e geolocation system with widely separated arrays has been developed and extensively applied for civil and military purpose during some considerable time. It can be used as the reserve navigation system for ships and airplanes [3], or as the reconnaissance and surveillance backup in the satellitechallenged environments [4].
Traditional OTH geolocation methods by widely separated arrays involve two-step processing [5]. In the first step, each observer array enables the estimation of the azimuth and elevation angles of the non-line-of-sight (NLOS) HF signal paths. Hereby, the well-known multiple signal classification (MUSIC) method [6] is usually used. In the second step, one kind of methods embeds the estimated azimuth and elevation angle measurements into the mathematical models of the ionosphere layers to obtain the location of the transmitter [7,8]. e authors in [7] addressed the joint estimation of the OTH target location and height by employing the diversity of multipath signal and structure of a 2D array. However, the atmospheric refraction effects and inaccuracies in the ionosphere model [9] that are frequently encountered will generate the larger location errors. To deal with the atmospheric refraction effects, the vertical atmospheric refractivity profile was modeled with a quadratic polynomial of the altitude [8] and then used to locate the target. Note that the refractivity parameters of this model are determined using known position information of auxiliary targets. Compared with these aforementioned methods, another kind of method is more practical and less complex in structure, which uses azimuth angles to locate the transmitter via a suitable bearings-only localization (BOL) algorithm [10] without the precise knowledge of the current ionosphere model. To further improve the performance, cooperative positioning has been studied with user collaboration [11,12], where the factor graph is often used to perform the location based on direction-of-arrival (DOA) and time-of-arrival (TOA) measurements. In cooperative positioning, the sensor nodes are synchronous and can communicate with each other through wireless links. However, in the OTH geolocation for long-range reconnaissance and surveillance, the sensors or observers are widely located, yielding the difficulty of communication between sensors and time synchronization. Note that the aforementioned two-step location methods cannot guarantee high localization precision as they extract measurements at each observer separately and independently, thus ignoring the constraint that the measurements correspond to the same transmitter location.
In contrast to the conventional two-step methods, direct position determination (DPD) exploits this constraint to improve performance. It uses the sampling data without estimating the intermediate parameters and ends with a cost function that relies on the locations of transmitters. e subspace data fusion (SDF) based DPD estimators were proposed in [13,14], in which the positions of multiple sources can be directly estimated from a cost function fusing the noise subspaces of all observer arrays. Another alternative of SDF-based DPD is the minimum-variance distortionless response (MVDR) based DPD method [15,16]. ese two kinds of DPD algorithms can localize multiple sources via a low-dimensional search, but their performance will be greatly degraded at low SNRs. As maximum likelihood (ML) estimators can approach the associated Cramér-Rao bounds (CRBs), many ML-based DPDs can be found in [5,17,18], where the location accuracy has shown to be significantly higher than that of the two-step methods, especially under low signal-to-noise ratio (SNR) conditions. However, when solving ML estimators without known waveforms, a large number of parameters yield a significant computation effort. Hence, an iterative decoupled algorithm for ML-based DPD was proposed for multiple sources to reduce the complexity [19], but it still requires a 2D-dimensional or 3D-dimensional search of the position. In order to make the DPD method more computationally efficient, an iterative scheme for DPD problem using Taylorseries method was developed in [20] to solve the position rather than the commonly used grid search. It is noteworthy that iterative solutions cannot guarantee to converge without an efficient initialization or under large noises. e importance sampling (IS) technique, a powerful Monte Carlo tool, has attracted great interest recently due to its superior performance on parameter estimation, such as time-difference-of-arrival (TDOA) and joint DOA-Doppler estimation and, more recently, the two-step location [21][22][23][24][25]. It can maximize the ML function in a computationally efficient manner by applying the global maximization theorem of Pincus [26]. In light of this, the IS technique can offer an efficient way to solve the ML-based DPD problems.
is paper investigates a DPD method using the IS technique for OTH geolocation of a transmitter observed by widely separated arrays. To the best of our knowledge, existing passive DPD methods do not consider the OTH geolocation by widely separated arrays. Our study employs the DPD technique to enhance the OTH geolocation of a HF transmitter and thus has considerable practical significance. Under the general assumption that the transmitter is placed on the ground and the signal waveforms are unknown, we formulate an ML-based function in terms of the transmitter position and elevation angles at all observers. As most existing DPDs are accomplished via an exhaustive grid search with high computation complexity or via iteration suffering from convergence problem, to make our method computationally attractive, we employ the IS technique to achieve an approximate global solution to the prescribed ML function, and parameters of interest are efficiently estimated in a noniterative manner. As our DPD method uses the IS technique to directly solve the transmitter position and elevation angles from the sampling data received by several observer arrays, which is more complicated compared with previous study, we derive an importance function analogous to Gaussian distribution, and therefore the required IS realizations can be easily generated according to this Gaussian distributed probability density function (PDF). Finally, we derive the CRB expression for the OTH DPD problem and test the localization performance of the proposed method through computer simulations.

Signal Model and Problem Formulation
Consider L stationary observers placed at positions u l ∈ R 3 (l � 1, 2, . . . , L), each of which is equipped with an antenna array composed of M isotropic sensors. A transmitter on the ground is assumed to radiate narrowband HF signals in the far field of the arrays. e signals are reflected by the ionosphere and received by the arrays. For simplicity, we only consider one dominant NLOS path component in our work since others may be heavily attenuated and can be ignored in some cases [27].
is scenario is described schematically in Figure 1. Let the ground plane be the Z plane, and thus the position of the transmitter can be denoted by a 3 × 1 vector of coordinates p � [p T , 0] T , whose Z-coordinate equals to zero as the transmitter is placed on the ground.
Assuming that the time duration T is sufficiently short during which the signal path remains unchanged, the observation r l (t), containing complex envelopes of the signals collected by the lth observer array at time 0 ≤ t ≤ T, is expressed as where s l (t) is the source signal at the input of the lth observer array, and the channel fading has been incorporated in s l (t). n l (t) ∈ C M×1 is the circularly white Gaussian noise vector mixed through the sensors. e sources and noises are assumed to be uncorrelated and to have a mean of zero. a l (θ l , φ l ) ∈ C M×1 is the array steering vector, where θ l and φ l represent the arrival azimuth and elevation angles of the signal at the lth observer array.
After being downconverted to the base band and sampled at t � nT s where T s is the sampling period, the nth sample of r l (t) can be expressed as where N represents the number of samples. s l (n) and n l (n) are the nth samples of s l (t) and n l (t), respectively. Without apriori information of the ionosphere model, the elevation angle has no evident relationship with the position, whereas the azimuth angle is still a function of the observer and transmitter positions: erefore, we can parameterize a l (θ l , φ l ) by p and φ l and thus rewrite (2) by r l (n) � a l p, φ l s l (n) + n l (n), n � 1, 2, . . . , N.
(4) Given the signal model above, the problem that we address now is to efficiently estimate the position vector p comprising the X-coordinate and Y-coordinate in a single step from the samples r l (n) without explicitly computing the azimuth angles.

Proposed DPD
In this section, we first formulate an ML-based function for the OTH geolocation and then apply an efficient IS algorithm to solve the prescribed ML function.

Optimization Model.
According to (4), we get the MLbased function by a sum over NL terms, which is shown as follows: where φ � [φ 1 , φ 2 , . . . , φ L ] T consists of the arrival elevation angles at L observers and s contains all the sampled signals received by L observers. erefore, the joint ML estimate of p, φ, and s can be given by We first estimate s by minimizing the cost function and hence obtain After substituting the estimated s l (n) to (5), the cost function is reduced to Here We then concatenate the observed vectors and projectors at all observer arrays and thus form Consequently, the ML-based function is expressed in a more compact form: We notice that the solution to this function is not simple because of the stray parameter φ and nonlinearity of unknowns in this function. Considering that exhaustive search is impractical, we resort to the IS technique to solve our DPD problem.

Choice of Importance Function for DPD.
To avoid direct maximization, we come to introduce the global optimization based on Pincus theorem [26], which provides a means of Mathematical Problems in Engineering performing the nonlinear optimization and can guarantee to find the global maximum. According to the Pincus theorem [26], the global minimizer of the ML estimator in (10) can be implemented by contains the estimated parameters. Let us define as a pseudo-PDF, and thus (11) can be rewritten as In practice, it is difficult to directly evaluate the (2 + L)-dimensional integrals. To make it tractable, we trade this problem by using the IS technique. Assuming that q(η) is an importance function with respect to η, the integral in (13) can be approximated by the following averaging [21][22][23][24]: (14) where w(η) � w(η)/ K k�1 w(η k ) is the normalized importance weight with w(η) � p(η)/q(η). η k is the kth realization of η distributed according to q(η), and K is the realization number. Usually, we want to choose q(η) as a simple function of η so that the realization of η can be easily generated. Moreover, q(η) should be chosen close to p(η) to improve the efficiency of IS sampling [23]. Inspired by this, to find an appropriate importance function for our DPD problem, we will derive an importance function analogous to Gaussian distribution, thus making IS realizations easily generated according to the Gaussian distributed PDF.
To begin with the derivation, initial estimates of the transmitter location and elevation angles are given, which are denoted by η (0) . en, using the conclusion in [28], we can approximate Π ⊥ (η)r(n) by the first-order Taylor expansion as follows: where G n (η) � (z Π ⊥ (η)r(n) /zη T ), h n (η) � G n (η)η− Π ⊥ (η)r(n), and o(‖η − η (0) ‖ 2 ) signifies the infinitesimal term of ‖η − η (0) ‖ 2 . Applying (15) to (10) leads to Defining we discard the terms that result in contributions of o(‖η − η (0) ‖ 2 ) and thus obtain As we attempt to construct an importance function q(η)with respect to the real parameter η, we should transform the terms in f ML (η) to be real-valued. erefore, (17) is represented as with G(η (0) ) � [Re T G(η (0) ) , Im T G(η (0) ) ] T and h(η (0) ) � [Re T h(η (0) ) , Im T h(η (0) ) ] T . Because it is desired to find q(η) that is close to p(η) in (12), according to the above approximation, the importance function can be constructed as where G and h are short for G(η (0) ) and h(η (0) ). As substituting (20) to (19) leads to can be regarded as a constant because its expression has no unknowns. As we attempt to generate η distributed according to q(η), q(η) is a kind of PDF, whose integral should be equal to one. By comparing (21) with the expression of Gaussian PDF, we make ����������������� det ((G T G) − 1 /2λ 1 ) ), and thus q(η) is equivalent to a Gaussian PDF: whose mean and covariance matrix are given by We observe that λ 1 is a parameter to adjust the covariance of the Gaussian distribution and then impacts the dispersion degree of the IS samples. Given the conclusion in [24] which reveals that the choice of λ 1 appears to be sufficiently robust to the SNR when the performance approaches the optimum, we just choose λ 1 � 0.5.
Considering that the magnitudes and units of the position and angles are quite different, it is reasonable to decouple them because their importance weights have different effect on the final results. erefore, we divide q(η) into two parts. e importance function for the transmitter position is constructed as and the importance function for the elevation angles is constructed as It should be noted that the choice of the above Gaussian distributions enjoys the following features: (1) ey are determined only by means and covariances, and therefore p and φ can be generated easily and separately (2) e joint PDF of p and φ is similar to the function p(η), which reduces the variance of the IS estimates 3.3. DPD Using Importance Sampling. After p and φ are generated based on the derived q p (p) and q φ (φ), according to w(η) � p(η)/q(η), the importance weights for the transmitter position and elevation angles can be computed by where p k and φ k are the kth samples from their Gaussian importance functions. To avoid the overflow since both the numerator and denominator are exponentials, we can replace w p (η k ) and w φ (η k ) with w p ′ (η k ) and w φ ′ (η k ) as follows:

Mathematical Problems in Engineering
Denote the normalized w p ′ (η k ) and w φ ′ (η k ) by w p ′ (η k ) and w φ ′ (η k ), respectively. Applying w p ′ (η k ) and w φ ′ (η k ) to (14), the position and elevation angles can be estimated by Remark 1. Based on the Pincus theorem, the global optimum can be obtained when λ ⟶ ∞, but this is impossible to achieve in practice. Fortunately, we can approximate λ to a sufficiently large value. Simulations will show that the location performance is not sensitive to the value of λ as long as it is sufficiently large and does not lead to numerical problem.
Remark 2. Note that we have performed the first-order linearization using Taylor expansion to find an importance function which approximates the pseudo-PDF related with ML-based function, and this will not have negative effect on the estimation performance since the objective function (i.e., ML-based function) is not approximated, which still plays a significant role as used in the importance weight (see (29)). Furthermore, our method differs greatly from the iterative methods like Taylor-series method [29], where the Taylor expansion is used to directly linearize the objective function and thus estimate the final result in an iterative manner.

Deterministic CRB
Relying on the considered model as shown in (4), we introduce the parameter vector ρ comprising all the unknown real parameters for this DPD problem: where σ 2 n denotes the noise power. e likelihood function for the N samples received by the L observers is expressed as en, according to the conclusion in [17], the CRB for parameter ω � [η T , Re T s { }, Im T s { }] T is given by Re za l p, φ l s l (n) zω T H za l p, φ l s l (n) For the convenience of derivation, we define a new vector x(ω) � [a T 1 (p, φ 1 )s 1 (n), . . . , a T L (p, φ L )s L (n), . . . , a T 1 (p, φ 1 )s 1 (N), . . . , a T L (p, φ L )s L (N)] T , and thus rewrite (33) by in which Σ ω � (zx(ω)/zω T ), and it can be divided into where Σ η , Σ Re s { } , and Σ Im s { } are defined similar to Σ ω . As the unknown waveforms are stray parameters in our study, to reduce the CRB of stray parameters and thus to obtain a more compact CRB matrix, we apply the scheme of the CRB derivation in [30]. Although the scenario considered in [30] is different, where the known signals are assumed to propagate in multipath environments and are observed by a single array, we can follow the steps similar to those in [30] and get the CRB for the transmitter position and elevation angles as where Π ⊥ Σ Re s { } is the orthogonal projector onto the space of Σ Re s { } .

Results and Discussion
e purpose of this section is to examine the performance of the proposed DPD method by comparing with the two conventional two-step localization methods and the DPD approach using the iterative scheme proposed in [20]. e simulated two-step methods first estimate azimuth and elevation angles at each observer independently using the 2dimensional (2D) MUSIC algorithm and then determine the source positions based on the estimated azimuth angles at all observers using Taylor-series (TS) method [29] and pseudolinear weighted least square (PLWLS) method [10], respectively.
e initial values for the proposed and TS methods are provided by the PLWLS estimator. e simulations are conducted upon the application of three stationary observers. Without loss of generality, we assume that three observers are located on the ground, whose e baseband signal waveforms are generated as narrowband signals, and channel fading is inversely proportional to the squared distance. We collect N � 200 samples of signals to implement the location, and each simulation performed 500 Monte Carlo runs.
As λ and λ 1 are significant parameters in the IS method, we first come to examine the effect of their values on the performance of the proposed method, and thus the choices of λ and λ 1 can be determined. e realization number K � 1000 is employed in the IS method. At SNR � 0 dB and λ 1 � 0.5, we vary the value of λ from 0 to 2000 and depict the root mean square error (RMSE) of the transmitter position and the average RMSE of three arrival elevation angles in Figure 2. It can be seen that the RMSEs diminish dramatically for λ < 200, whereas the performance remains almost unchanged as λ becomes sufficiently large. is confirms the analysis in Remark 1. Inspired by this result, we choose λ � 1000 in the following simulations. Under the same SNR, we fix λ at 1000 and change λ 1 from 0.1 to 3.3. e corresponding results are illustrated in Figure 3, which shows that the performance of the proposed algorithm is not evidently sensitive to different values of λ 1 . On this basis and without loss of generality, the choice of λ 1 � 0.5 is used in our study.
en, as the SNR varies from −12 dB to 10 dB, we compare the estimation accuracy of our method with those of the twostep and the iterative DPD methods. e corresponding results are shown in Figure 4, where the derived CRBs are presented as the benchmark. We notice that the proposed DPD not only performs better than the TS and PLWLS methods in terms of location RMSE but also outperforms the 2D MUSIC algorithm in terms of the average RMSE of elevation angles at low SNRs. With an increase in SNR, our method can reach the associated CRB. Moreover, it can be seen from Figure 4 that our method using the IS technique has greater robustness performance than the iterative DPD method in the low SNR region. is may be due to the fact that our method solves the ML problem in a noniterative manner, whereas the iterative DPD method cannot guarantee to converge when the noise is sufficiently large. For different SNRs, we evaluate the average computer running time using a laptop with a 2.5 GHz Intel Core CPU. Each prescribed location method and the exhaustive search implementation of the investigated problem is examined. As shown in Table 1, the complexity of the proposed DPD is much lower than those of the exhaustive search and the twostep methods. It significantly decreases the computational cost compared with the exhaustive search due to the   application of the IS technique and the derived Gaussian importance function. Although the proposed DPD costs more time than the iterative DPD when the SNR is larger than 0 dB, the running time of our method using the IS technique is not affected by the SNRs, whereas the iterative DPD requires more iterations and thus much longer time at low SNRs. Finally, we test the localization performance for different arrival elevation angles. When the SNR is −5 dB, the RMSE curves as well as the CRB curve are displayed in Figure 5 with each elevation angle changing from 10 to 80 degrees. We see that both the location RMSEs and the average RMSEs of elevation angles are altered with different arrival elevation angles. Specifically, the location RMSEs get smaller when elevation is low but larger when elevation is high, whereas the average RMSEs of elevation angles are opposite. is can be explained by the fact that all the simulated methods explicitly or implicitly exploit the information in the azimuth angles and the estimation of azimuth angles relies on the values of elevation angles. Notice that the RMSE of elevation angle estimated by iterative DPD deteriorates dramatically when the elevation angles are lower than 20 degree, but our algorithm still exhibits superior location performance over other methods at each elevation angle.

Conclusion
In this paper, we have proposed a single-step geolocation method for an OTH transmitter observed by widely separated arrays. e proposed method directly estimates the transmitter position and elevation angles without the need for prior information of the ionosphere model. IS technique is employed in its solution instead of the commonly used grid search and iteration scheme, which makes our method practically attractive. Simulation results demonstrate that our method outperforms the conventional two-step approaches as well as the iterative DPD, and it asymptotically reaches the corresponding CRB as the SNR increases.
For the OTH geolocation, we can obtain better results if an ellipsoid model is used to describe the shape of the Earth. is can be achieved by adapting the proposed method to the geocentric coordinate system, and we defer it to the extended version of this paper. In this paper, the boldface italic upper case letter denotes the matrix and the boldface italic lower case letter signifies the vector. For convenience, we list the notations used in this paper.