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

Authors: Jiexin Yin; Ding Wang; Bin Yang; Xin Yang
Journal: Mathematical Problems in Engineering
 (2021)

---

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


---

## Body

## 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

ConsiderL stationary observers placed at positions ul∈ℝ3l=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. 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 p˜=pT,0T, whose Z-coordinate equals to zero as the transmitter is placed on the ground.Figure 1
Scenario of the OTH transmitter observed by multiple arrays.Assuming that the time durationT is sufficiently short during which the signal path remains unchanged, the observation rlt, containing complex envelopes of the signals collected by the lth observer array at time 0 ≤ t ≤ T, is expressed as(1)rlt=alθl,φlslt+nlt,l=1,2,…,L,where sl(t) is the source signal at the input of the lth observer array, and the channel fading has been incorporated in sl(t). nlt∈ℂM×1 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. alθl,φl∈ℂ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 att = nTs where Ts is the sampling period, the nth sample of rlt can be expressed as(2)rln=alθl,φlsln+nln,n=1,2,…,N,where N represents the number of samples. sl(n) and nln are the nth samples of sl(t) and nlt, 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:(3)θl=arctanp2−ul2p1−ul1,l=1,2,…,L.Therefore, we can parameterizealθl,φl by p and φl and thus rewrite (2) by(4)rln=alp,φlsln+nln,n=1,2,…,N.Given the signal model above, the problem that we address now is to efficiently estimate the position vectorp comprising the X-coordinate and Y-coordinate in a single step from the samples rln 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:(5)fMLp,φ,s=∑l=1L∑n=1Nrln−alp,φlsln2,where φ=φ1,φ2,…,φLT consists of the arrival elevation angles at L observers and s contains all the sampled signals received by L observers. Therefore, the joint ML estimate of p,φ,ands can be given by(6)p^,φ^,s^=argminp,φ,sfMLp,φ,s.We first estimates by minimizing the cost function and hence obtain(7)s^ln=al†p,φlrln.After substituting the estimateds^ln to (5), the cost function is reduced to(8)fMLp,φ=∑l=1L∑n=1NΠal⊥p,φlrln2.Here,Πal⊥p,φl=IM−al†p,φlalp,φl is the orthogonal projector onto the space of alp,φl with IM being the M×M identity matrix.We then concatenate the observed vectors and projectors at all observer arrays and thus form(9)rn=r1Tn,r2Tn,…,rLTnT,Π⊥p,φ=blkdiagΠa1⊥p,φ1,Πa2⊥p,φ2,…,ΠaL⊥p,φL.Consequently, the ML-based function is expressed in a more compact form:(10)fMLp,φ=∑n=1NΠ⊥p,φrn2.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 by(11)η^i=limλ⟶∞∫⋯∫ηiexp−λfMLηdη∫⋯∫exp−λfMLηdη,where η=pT,φTT∈ℝL+2×1 contains the estimated parameters. Let us define(12)pη=exp−λfMLη∫⋯∫exp−λfMLηdηas a pseudo-PDF, and thus (11) can be rewritten as(13)η^i=limλ⟶∞∫⋯∫ηipηdη.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–24]:(14)η^i=limλ⟶∞∫⋯∫ηipηqηqηdη≈∑k=1Kηkiw˜ηk,where w˜η=wη/∑k=1Kwη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. Then, using the conclusion in [28], we can approximate Π⊥ηrn by the first-order Taylor expansion as follows:(15)Π⊥ηrn≈Π⊥η^0rn+Gnη^0η−η^0+oη−η^02≈Gnη^0η−hnη^0+oη−η^02,where Gnη=∂Π⊥ηrn/∂ηT, hnη=Gnηη−Π⊥ηrn, and oη−η^02 signifies the infinitesimal term of η−η^02.Applying (15) to (10) leads to(16)fMLη=∑n=1NGnη^0η−hnη^02+oη−η^02.DefiningGη=G1Tη,G2Tη,…,GNTηT and hη=h1Tη,h2Tη,…,hNTηT, we discard the terms that result in contributions of oη−η^02 and thus obtain(17)fMLη≈Gη^0η−hη^02.As we attempt to construct an importance functionqηwith respect to the real parameter η, we should transform the terms in fMLη to be real-valued. Therefore, (17) is represented as(18)fMLη≈ReGη^0η−Rehη^02+ImGη^0η−Imhη^02≈G˜η^0η−h˜η^02,with G˜η^0=ReTGη^0,ImTGη^0T and h˜η^0=ReThη^0,ImThη^0T. 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(19)qη=exp−λ1G˜η−h˜TG˜η−h˜∫⋯∫exp−λ1G˜η−h˜TG˜η−h˜dη,where G˜ and h˜ are short for G˜η^0 and h˜η^0. As(20)G˜η−h˜TG˜η−h˜=η−G˜TG˜−1G˜Th˜TG˜TG˜η−G˜TG˜−1G˜Th˜+h˜Th˜−h˜TG˜G˜TG˜−1G˜Th˜,substituting (20) to (19) leads to(21)qη=αexp−λ1η−G˜TG˜−1G˜Th˜TG˜TG˜η−G˜TG˜−1G˜Th˜,where α=exp−λ1h˜Th˜−h˜TG˜G˜TG˜−1G˜Th˜/∫⋯∫exp−λ1G˜η−h˜TG˜η−h˜dη 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 α=1/2πL/2+1detG˜TG˜−1/2λ1, and thus qη is equivalent to a Gaussian PDF:(22)qη=12πL/2+1detR1/2exp−12η−η¯TR−1η−η¯,whose mean and covariance matrix are given by(23)η¯=G˜TG˜−1G˜Th˜,R=G˜TG˜−12λ1.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 divideqη into two parts. The importance function for the transmitter position is constructed as(24)qpp=12πdetRp1/2exp−12p−p¯TRp−1p−p¯,with(25)p¯=η¯1:2,Rp=R1:2,1:2,and the importance function for the elevation angles is constructed as(26)qφφ=12πL/2detRφ1/2exp−12φ−φ¯TRφ−1φ−φ¯,with(27)φ¯=η¯3:2+L,Rφ=R3:2+L,3:2+L.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 thereforep and φ can be generated easily and separately(2)
The joint PDF ofp and φ is similar to the function pη, which reduces the variance of the IS estimates


### 3.3. DPD Using Importance Sampling

Afterp and φ are generated based on the derived qpp and qφφ, according to wη=pη/qη, the importance weights for the transmitter position and elevation angles can be computed by(28)wpηk∝exp−λfMLηk+12pk−p¯TRp−1pk−p¯,wφηk∝exp−λfMLηk+12φk−φ¯TRφ−1φk−φ¯,where pk 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 wpηk and wφηk with wp′ηk and wφ′ηk as follows:(29)wp′ηk=exp−λfMLηk+12pk−p¯TRp−1pk−p¯−maxl=1,2,…,K−λfMLηl+12pl−p¯TRp−1pl−p¯,wφ′ηk=exp−λfMLηk+12φk−φ¯TRφ−1φk−φ¯−maxl=1,2,…,K−λfMLηl+12φl−φ¯TRφ−1φl−φ¯.Denote the normalizedwp′η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(30)p^=∑k=1Kpkw˜p′ηk,φ^=∑k=1Kφkw˜φ′ηk.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.


## 3.1. Optimization Model

According to (4), we get the ML-based function by a sum over NL terms, which is shown as follows:(5)fMLp,φ,s=∑l=1L∑n=1Nrln−alp,φlsln2,where φ=φ1,φ2,…,φLT consists of the arrival elevation angles at L observers and s contains all the sampled signals received by L observers. Therefore, the joint ML estimate of p,φ,ands can be given by(6)p^,φ^,s^=argminp,φ,sfMLp,φ,s.We first estimates by minimizing the cost function and hence obtain(7)s^ln=al†p,φlrln.After substituting the estimateds^ln to (5), the cost function is reduced to(8)fMLp,φ=∑l=1L∑n=1NΠal⊥p,φlrln2.Here,Πal⊥p,φl=IM−al†p,φlalp,φl is the orthogonal projector onto the space of alp,φl with IM being the M×M identity matrix.We then concatenate the observed vectors and projectors at all observer arrays and thus form(9)rn=r1Tn,r2Tn,…,rLTnT,Π⊥p,φ=blkdiagΠa1⊥p,φ1,Πa2⊥p,φ2,…,ΠaL⊥p,φL.Consequently, the ML-based function is expressed in a more compact form:(10)fMLp,φ=∑n=1NΠ⊥p,φrn2.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 by(11)η^i=limλ⟶∞∫⋯∫ηiexp−λfMLηdη∫⋯∫exp−λfMLηdη,where η=pT,φTT∈ℝL+2×1 contains the estimated parameters. Let us define(12)pη=exp−λfMLη∫⋯∫exp−λfMLηdηas a pseudo-PDF, and thus (11) can be rewritten as(13)η^i=limλ⟶∞∫⋯∫ηipηdη.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–24]:(14)η^i=limλ⟶∞∫⋯∫ηipηqηqηdη≈∑k=1Kηkiw˜ηk,where w˜η=wη/∑k=1Kwη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. Then, using the conclusion in [28], we can approximate Π⊥ηrn by the first-order Taylor expansion as follows:(15)Π⊥ηrn≈Π⊥η^0rn+Gnη^0η−η^0+oη−η^02≈Gnη^0η−hnη^0+oη−η^02,where Gnη=∂Π⊥ηrn/∂ηT, hnη=Gnηη−Π⊥ηrn, and oη−η^02 signifies the infinitesimal term of η−η^02.Applying (15) to (10) leads to(16)fMLη=∑n=1NGnη^0η−hnη^02+oη−η^02.DefiningGη=G1Tη,G2Tη,…,GNTηT and hη=h1Tη,h2Tη,…,hNTηT, we discard the terms that result in contributions of oη−η^02 and thus obtain(17)fMLη≈Gη^0η−hη^02.As we attempt to construct an importance functionqηwith respect to the real parameter η, we should transform the terms in fMLη to be real-valued. Therefore, (17) is represented as(18)fMLη≈ReGη^0η−Rehη^02+ImGη^0η−Imhη^02≈G˜η^0η−h˜η^02,with G˜η^0=ReTGη^0,ImTGη^0T and h˜η^0=ReThη^0,ImThη^0T. 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(19)qη=exp−λ1G˜η−h˜TG˜η−h˜∫⋯∫exp−λ1G˜η−h˜TG˜η−h˜dη,where G˜ and h˜ are short for G˜η^0 and h˜η^0. As(20)G˜η−h˜TG˜η−h˜=η−G˜TG˜−1G˜Th˜TG˜TG˜η−G˜TG˜−1G˜Th˜+h˜Th˜−h˜TG˜G˜TG˜−1G˜Th˜,substituting (20) to (19) leads to(21)qη=αexp−λ1η−G˜TG˜−1G˜Th˜TG˜TG˜η−G˜TG˜−1G˜Th˜,where α=exp−λ1h˜Th˜−h˜TG˜G˜TG˜−1G˜Th˜/∫⋯∫exp−λ1G˜η−h˜TG˜η−h˜dη 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 α=1/2πL/2+1detG˜TG˜−1/2λ1, and thus qη is equivalent to a Gaussian PDF:(22)qη=12πL/2+1detR1/2exp−12η−η¯TR−1η−η¯,whose mean and covariance matrix are given by(23)η¯=G˜TG˜−1G˜Th˜,R=G˜TG˜−12λ1.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 divideqη into two parts. The importance function for the transmitter position is constructed as(24)qpp=12πdetRp1/2exp−12p−p¯TRp−1p−p¯,with(25)p¯=η¯1:2,Rp=R1:2,1:2,and the importance function for the elevation angles is constructed as(26)qφφ=12πL/2detRφ1/2exp−12φ−φ¯TRφ−1φ−φ¯,with(27)φ¯=η¯3:2+L,Rφ=R3:2+L,3:2+L.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 thereforep and φ can be generated easily and separately(2)
The joint PDF ofp and φ is similar to the function pη, which reduces the variance of the IS estimates


## 3.3. DPD Using Importance Sampling

Afterp and φ are generated based on the derived qpp and qφφ, according to wη=pη/qη, the importance weights for the transmitter position and elevation angles can be computed by(28)wpηk∝exp−λfMLηk+12pk−p¯TRp−1pk−p¯,wφηk∝exp−λfMLηk+12φk−φ¯TRφ−1φk−φ¯,where pk 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 wpηk and wφηk with wp′ηk and wφ′ηk as follows:(29)wp′ηk=exp−λfMLηk+12pk−p¯TRp−1pk−p¯−maxl=1,2,…,K−λfMLηl+12pl−p¯TRp−1pl−p¯,wφ′ηk=exp−λfMLηk+12φk−φ¯TRφ−1φk−φ¯−maxl=1,2,…,K−λfMLηl+12φl−φ¯TRφ−1φl−φ¯.Denote the normalizedwp′η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(30)p^=∑k=1Kpkw˜p′ηk,φ^=∑k=1Kφkw˜φ′ηk.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:(31)ρ=ηT,ReTs,ImTs,σn2T,where σn2 denotes the noise power. The likelihood function for the N samples received by the L observers is expressed as(32)fρ=1πσn2MLNexp−1σn2∑l=1L∑n=1Nrln−alp,φlsln2.Then, according to the conclusion in [17], the CRB for parameter ω=ηT,ReTs,ImTsTis given by(33)CRBω−1=2σn2∑l=1L∑n=1NRe∂alp,φlsln∂ωTH∂alp,φlsln∂ωT.For the convenience of derivation, we define a new vectorxω=a1Tp,φ1s1n,…,aLTp,φLsLn,…,a1Tp,φ1s1N,…,aLTp,φLsLNT, and thus rewrite (33) by(34)CRBω=σn22Re−1∂xω∂ωTH∂xω∂ωT=σn22Re−1ΣωHΣω,in which Σω=∂xω/∂ωT, and it can be divided into(35)Σω=Ση,ΣRes,ΣIms=Ση,ΣRes,jΣRes,where Ση, ΣRes, and ΣIms 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(36)CRBη=σn22Re−1ΣηHΠΣRes⊥Ση,where ΠΣRes⊥ is the orthogonal projector onto the space of ΣRes.


## 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 φ=50,55,70T (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.Figure 2
RMSEs of the proposed method with different values ofλ with λ1 = 0.5. (a) Location RMSE. (b) Average RMSE of elevation angles.
(a)(b)Figure 3
RMSEs of the proposed method with different values ofλ1 with λ = 1000. (a) Location RMSE. (b) Average RMSE of elevation angles.
(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 Figure4, 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.Figure 4
Comparison of RMSEs and CRB versus SNR. (a) Location RMSE. (b) Average RMSE of elevation angles.
(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 Table1, 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.Table 1
Average running time in seconds for different SNRs.
Proposed DPDExh. searchIterative DPDMUSIC + TSMUSIC + PLWLSSNR = −10 dB0.451291.04011.14984.66034.3123SNR = 0 dB0.427690.73980.32684.62084.3289SNR = 10 dB0.418490.56120.13634.53984.3326Finally, 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 Figure5 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.Figure 5
Comparison of RMSEs and CRB for different elevation angles. (a) Location RMSE. (b) Average RMSE of elevation angles.
(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.



---

*Source: 5574110-2021-11-26.xml*