Abstract
This 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. The 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. This 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. The simulation results corroborate the superior localization performance of the proposed method with respect to the conventional two-step approaches and the iterative DPD method.
1. 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. The 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 satellite-challenged 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]. The 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. The 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]. These 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 Taylor-series 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.
The 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–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.
This 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.
2. Signal Model and Problem Formulation
Consider L stationary observers placed at positions , 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. The 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]. This 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 , 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 , containing complex envelopes of the signals collected by the lth observer array at time 0 ≤ t ≤ T, is expressed aswhere sl(t) is the source signal at the input of the lth observer array, and the channel fading has been incorporated in sl(t). is the circularly white Gaussian noise vector mixed through the sensors. The sources and noises are assumed to be uncorrelated and to have a mean of zero. is the array steering vector, where θl and 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 = nTs where Ts is the sampling period, the nth sample of can be expressed aswhere N represents the number of samples. sl(n) and are the nth samples of sl(t) and , 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:
Therefore, we can parameterize by and and thus rewrite (2) by
Given the signal model above, the problem that we address now is to efficiently estimate the position vector comprising the X-coordinate and Y-coordinate in a single step from the samples without explicitly computing the azimuth angles.
3. 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.
3.1. Optimization Model
According to (4), we get the ML-based function by a sum over NL terms, which is shown as follows:where consists of the arrival elevation angles at L observers and contains all the sampled signals received by L observers. Therefore, the joint ML estimate of can be given by
We first estimate by minimizing the cost function and hence obtain
After substituting the estimated to (5), the cost function is reduced to
Here, is the orthogonal projector onto the space of with being the identity matrix.
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.
3.2. 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 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 bywhere contains the estimated parameters. Let us defineas 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 is an importance function with respect to , the integral in (13) can be approximated by the following averaging [21–24]:where is the normalized importance weight with . is the kth realization of distributed according to , and K is the realization number. Usually, we want to choose as a simple function of so that the realization of can be easily generated. Moreover, should be chosen close to 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 . Then, using the conclusion in [28], we can approximate by the first-order Taylor expansion as follows:where , , and signifies the infinitesimal term of .
Applying (15) to (10) leads to
Defining and , we discard the terms that result in contributions of and thus obtain
As we attempt to construct an importance function with respect to the real parameter , we should transform the terms in to be real-valued. Therefore, (17) is represented aswith and . Because it is desired to find that is close to in (12), according to the above approximation, the importance function can be constructed aswhere and are short for and . Assubstituting (20) to (19) leads towhere can be regarded as a constant because its expression has no unknowns. As we attempt to generate distributed according to , is a kind of PDF, whose integral should be equal to one. By comparing (21) with the expression of Gaussian PDF, we make , and thus 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. Therefore, we divide into two parts. The importance function for the transmitter position is constructed aswithand the importance function for the elevation angles is constructed aswith
It should be noted that the choice of the above Gaussian distributions enjoys the following features:(1)They are determined only by means and covariances, and therefore and can be generated easily and separately(2)The joint PDF of and is similar to the function , which reduces the variance of the IS estimates
3.3. DPD Using Importance Sampling
After and are generated based on the derived and , according to , the importance weights for the transmitter position and elevation angles can be computed bywhere and are the kth samples from their Gaussian importance functions. To avoid the overflow since both the numerator and denominator are exponentials, we can replace and with and as follows:
Denote the normalized and by and , respectively. Applying and 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.
4. 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 denotes the noise power. The likelihood function for the N samples received by the L observers is expressed as
Then, according to the conclusion in [17], the CRB for parameter is given by
For the convenience of derivation, we define a new vector , and thus rewrite (33) byin which , and it can be divided intowhere , , and 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 aswhere is the orthogonal projector onto the space of .
5. Results and Discussion
The 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]. The simulated two-step methods first estimate azimuth and elevation angles at each observer independently using the 2-dimensional (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 pseudo-linear weighted least square (PLWLS) method [10], respectively. The initial values for the proposed and TS methods are provided by the PLWLS estimator.
The 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 positions are [−1000, 1000, 0]T (km), [1000,1000,0]T (km), and [1000,−1000,0]T (km), respectively. Each observer is equipped with a uniform circular array (UCA) comprising M = 9 sensors in the X-Y plane. The radius of the UCA is 1.5λ, where λ denotes the wavelength. The HF transmitter is placed at [200, 500, 0]T [km]. The transmitted signals are reflected by the ionosphere and arrive at observers with elevation angles (degree). The 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. The 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. This 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. The 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.

(a)

(b)

(a)

(b)
Then, as the SNR varies from −12 dB to 10 dB, we compare the estimation accuracy of our method with those of the two-step and the iterative DPD methods. The 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. This 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.

(a)

(b)
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 two-step 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. This 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.

(a)

(b)
6. Conclusion
In this paper, we have proposed a single-step geolocation method for an OTH transmitter observed by widely separated arrays. The 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. This 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.
Abbreviation
OTH: | Over-the-horizon |
NLOS: | Non-line-of-sight |
HF: | High-frequency |
MUSIC: | Multiple signal classification |
DOA: | Direction-of-arrival |
TOA: | Time-of-arrival |
BOL: | Bearings-only localization |
DPD: | Direct position determination |
ML: | Maximum likelihood |
SNR: | Signal-to-noise ratio |
IS: | Importance sampling |
TDOA: | Time-difference-of-arrival |
PDF: | Probability density function |
CRB: | Cramér–Rao bound |
TS: | Taylor series |
PLWLS: | Pseudo-linear weighted least square |
RMSE: | Root mean square error. |
: | Transpose |
: | Conjugate transpose |
: | Moore–Penrose inverse |
: | Composition of block diagonal matrix |
: | Euclidean norm |
: | Sets of the real-valued matrix |
: | Sets of the complex-valued matrix |
: | Real part |
: | Imaginary part |
: | Determinant |
: | The orthogonal projector matrix onto the null space of the matrix |
: | The th element of the vector. |
Data Availability
The simulation data implemented on MATLAB used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 61901526, 62071029, 62171469, and 61772548), China Postdoctoral Science Foundation (Grant nos. 2016M592989 and 2019M663997), and the Key Scientific and Technological Research Project in Henan Province (Grant nos. 192102210092 and 192102210117).