A Computational Efficient Maximum Likelihood Direct Position Determination Approach for Multiple Emitters Using Angle and Doppler Measurements

The unknown parameters of interest are given by

where ${\bf{b}}_{k,l}\buildrel\Delta\over{=}{\left[{{b_{1,k,l}},...,{b_{q,k,l}}}\right]^{T}}$. We focus on the ML estimator that depends on ${{{\bf{\bar{p}}}}}$ and ${{{\bf{\bar{b}}}}}$. Since ${\bf{\bar{w}}}_{l,k}\;\backsim\;\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I}_{NM})$, the ML estimator coincides with the Least Squares (LS) estimator. Thus the log-likelihood function (after dropping the constant terms) can be given by

where ${{\bf{\bar{D}}}_{l,k}}\buildrel\Delta\over{=}{\left[{\bf D}_{l,k}\left({{{\bf{p}}_{1}}}\right){\bf s}_{1,k},...,{\bf D}_{l,k}\left({{{\bf{p}}_{Q}}}\right){\bf s}_{Q,k}\right]}$. The path attenuation vectors that maximize (9) are given by

Substituting (10) back into (9) yields the so-called concentrated likelihood function (CLF) which depends solely on ${{{\bf{\bar{p}}}}}$

Instead of maximizing (11), we can equivalently maximize

Hence, the ML estimate of ${{{\bf{\bar{p}}}}}$ is then obtained from the following multi-dimensional optimization problem

A direct solution of (13) requires $Q\times D$-dimensional search, which is extremely computational challenging. To solve this problem, we resort to the Pincus’ theorem [25]. It provides a mean for performing the nonlinear multi-dimensional optimization and guarantees to produce the global maximum. The Pincus’ theorem states that the vector $\hat{\boldsymbol{\theta}}=[{{\hat{\theta}}_{1}},{{\hat{\theta}}_{2}},\ldots,{{\hat{\theta}}_{\widetilde{Q}}}]$ yields the unique global maximum of a continuous ${\widetilde{Q}}$-dimensional function $f\left({\boldsymbol{\theta}}\right)$, whose ${\tilde{q}}$-th entry is given by

Using a sufficiently large ${\rho}_{0}$ to replace ${\rho}$, the limit in (14) is approximated as

Applying (15) to the optimization problem (13) with ${\boldsymbol{\theta}}={\bf{\bar{p}}}$ and $f\left({\boldsymbol{\theta}}\right)={L_{2}}\left({{\bf{\bar{p}}}}\right)$, and letting $D=2$ for simplifying the exhibition, the expressions for the estimation of the ${{\bf{p}}_{q}}$ are given by

where ${{p}_{q,x}}$ and ${{p}_{q,y}}$ denote x-coordinate and y-coordinate of the $q$-th emitter, respectively. ${\bar{L}_{2}}({\bf{\bar{p}}})$ is the normalized function of ${{e^{{\rho}_{0}{L_{2}}({\bf{\bar{p}}})}}}$ defined as

However, it is difficult to directly deal with the multi-dimensional integral (16a) and (16b). Closely inspecting (17), ${\bar{L}_{2}}({\bf{\bar{p}}})$ satisfies ${\displaystyle\int\cdots\int{{\bar{L}_{2}}({\bf{\bar{p}}})}d{\bf{\bar{p}}}}=1$. Considering ${\bar{L}_{2}}({\bf{\bar{p}}})$ as a pseudo-PDF of ${\bf{\bar{p}}}$, thus the estimation of ${{\bf{p}}_{q}}$ in (16a) and (16b) can be alternatively regarded as statistical expectations

Intuitively, we can approximate the expectations in (18) with the generation of $R$ realizations $\left\{{{{{\bf{\bar{p}}}}^{(r)}}}\right\}_{r=1}^{R}$ distributed according to ${\bar{L}_{2}}({\bf{\bar{p}}})$,

In this way, the complicated integrations can be approximated by the sample mean estimates. Clearly, according to the law of large numbers [36], the sample mean estimate will converge to the corresponding expectation as $R\to\infty$. Thus the increase of $R$ results in the decrease of the variance of ${{\hat{p}}_{q,x}}$ and ${{\hat{p}}_{q,y}}$. Unfortunately, ${\bar{L}_{2}}({\bf{\bar{p}}})$ is still highly non-linear and cannot be practically used to generate $\left\{{{{{\bf{\bar{p}}}}^{(r)}}}\right\}_{r=1}^{R}$. We shall resort to the IS concept and approximately find a simple PDF for realization generation.

Equations (16a) and (16b) are equivalent to the following forms respectively

where ${{{{\bar{L}}_{im}}({\bf{\bar{p}}})}}$ is the other pseudo-PDF called importance function [27] [37]. ${{{{\bar{L}}_{im}}({\bf{\bar{p}}})}}$ is generally designed as a simple function of ${\bf{\bar{p}}}$ for easily generating realizations $\left\{{{{{\bf{\bar{p}}}}^{(r)}}}\right\}_{r=1}^{R}$, which is also desired as similar as possible to ${{{{\bar{L}}_{2}}({\bf{\bar{p}}})}}$. Thus (20a) and (20b) are regarded as computing expectation of transformed random variables

where

Now we turn into the approximate choice for $\bar{L}_{im}(\bar{\mathbf{p}})$. In order to generate $\left\{{{{{\bf{\bar{p}}}}^{(r)}}}\right\}_{r=1}^{R}$ easily, intuitively ${{{{\bar{L}}_{im}}({\bf{\bar{p}}})}}$ shall be designed separable in terms of the $Q$ emitter locations

where ${{{\bar{\mathcal{L}}}_{q}}\left({{{\bf{p}}_{q}}}\right)}$ is a pseudo-PDF of ${{{\bf{p}}_{q}}}$. ${{{{\bar{L}}_{im}}({\bf{\bar{p}}})}}$ in form (23) can be regarded as a joint pseudo-PDF corresponding to multiple mutually independent random variables, which can be easily generated using $\left\{{{{\bar{\mathcal{L}}}_{q}}\left({{{\bf{p}}_{q}}}\right)}\right\}_{q=1}^{Q}$. Hence, considering the design of ${{{{\bar{L}}_{im}}({\bf{\bar{p}}})}}$ shall be appropriate approximation of ${{{{\bar{L}}_{2}}\left({{\bf{\bar{p}}}}\right)}}$ satisfies (23).

Revisiting (12), ${L_{2}}\left({\bf{\bar{p}}}\right)$ can be approximated as a separable function under the condition that ${{\bf{\bar{D}}}_{l,k}^{H}{{{\bf{\bar{D}}}}_{l,k}}}$ is replaced by an invertible diagonal matrix. It can be observed that the entry in the $i$-th row and the $g$-th column of matrix ${{\bf{\bar{D}}}_{l,k}^{H}{{{\bf{\bar{D}}}}_{l,k}}}$ can be expressed as

The diagonal elements are

We verify statistically that the off-diagonal entries of ${{\bf{\bar{D}}}_{l,k}^{H}{{{\bf{\bar{D}}}}_{l,k}}}$ are much smaller compared to $\left[{{\bf{\bar{D}}}_{l,k}^{H}{{{\bf{\bar{D}}}}_{l,k}}}\right]_{i,i}$ with a high probability for almost all possible emitter loactions. To this end, let us define

as the ratio of the off-diagonal entries over diagonal entries of ${{\bf{\bar{D}}}_{l,k}^{H}{{{\bf{\bar{D}}}}_{l,k}}}$. Then, a large number of random variable vectors ${{{\bf{p}}_{i}}}$ and ${{{\bf{p}}_{g}}}$ are generated, which are independent identically distributed (IID) in $U\left[{-100,100}\right]^{2}$ Km, to compute the complementary cumulative distribution function (CCDF) of ${\delta_{i,g}}$.

As shown in Fig.1, $\left[{{\bf{\bar{D}}}_{l,k}^{H}{{{\bf{\bar{D}}}}_{l,k}}}\right]_{i,i}$ is indeed dominant compared to corresponding off-diagonal elements, since ${\delta_{i,g}}$ almost has a zero probability to exceed 0.15. Therefore, a valid approximation for (24) can be obtained

Substituting (27) into (12), ${L_{2}}$ is approximately in proportion to the superposition of $Q$ separated terms

where

Then, the importance function (23) is transformed to

where ${\rho_{1}}$ is the new design parameter involving factor $\frac{1}{M}$. It should be noted that the choices of ${\rho_{0}}$ and ${\rho_{1}}$ affect the performance of the proposed estimator. ${{\bar{L}}_{im}}\left({{\bf{\bar{p}}}}\right)$ can indeed be factorized in terms of the $Q$ emitter locations as (23), where ${{{\bar{\mathcal{L}}}_{q}}\left({{{\bf{p}}}}\right)}$ is expressed as

Thus the generation of ${{{{\bf{\bar{p}}}}^{(r)}}}$ can be approximated as generating $\left\{{{{{\bf{p}}}^{(r)}_{q}}}\right\}_{q=1}^{Q}$ with bivariate distribution (31), which can be implemented separately and run in parallel with a much faster and less complex execution.

One specific way for generating required vector realizations is based on the inverse probability transformation applied in [36]. Let ${F_{\bf{X}}}\left({\bf{x}}\right)$ be the cumulative distribution function (CDF) of random variable vector ${\bf{X}}$ and generating $R$ IID uniform random numbers, i.e., $\left\{{{{{u}}^{(r)}}}\right\}_{r=1}^{R}\in U\left[{0,1}\right]$. Thus the variable vector realization ${{{{{\bf{x}}}}^{(r)}}}$ can be generated as

However, this way still requires a two-dimensional search process for generating ${{{{\bf{p}}}^{(r)}_{q}}}$. Recalling (31) that depicts a joint distribution, thus ${\bar{\cal L}_{q}}\left({\bf{p}}\right)$ can be factorized as the product of a marginal PDF and a conditional PDF as follows

where ${\bar{\cal L}_{p_{q,x}}}\left({p_{x}}\right)$ is the marginal PDF of x-coordinate $p_{q,x}$ and ${\bar{\cal L}_{{p_{q,y}}\left|{{p_{q,x}}}\right.}}\left({{p_{y}}\left|{{p_{x}}}\right.}\right)$ is the conditional PDF of y-coordinate $p_{q,y}$ given $p_{q,x}$. The definition of ${\bar{\cal L}_{p_{q,y}}}\left({p_{y}}\right)$ and ${\bar{\cal L}_{{p_{q,x}}\left|{{p_{q,y}}}\right.}}\left({{p_{x}}\left|{{p_{y}}}\right.}\right)$ is in the same manner. ${\bar{\cal L}_{p_{q,x}}}\left({p_{x}}\right)$ can be given by

which is used to generate the realization ${{{{{p}}}_{q,x}^{(r)}}}$. Then the conditional PDF of ${{p_{q,y}}}$ given ${{{{{p}}}_{q,x}^{(r)}}}$ can be expressed as following

By (35) the realization ${{{{{p}}}_{q,y}^{(r)}}}$ is generated, thus the two-dimensional search process for generating ${{{{\bf{p}}}^{(r)}_{q}}}=[{{{{{p}}}_{q,x}^{(r)}}},{{{{{p}}}_{q,y}^{(r)}}}]^{T}$ is converted to two line search. We have just showed the process for generating ${{{{{p}}}_{q,x}^{(r)}}}$ firstly and then ${{{{{p}}}_{q,y}^{(r)}}}$. The generation of ${{{{{p}}}_{q,y}^{(r)}}}$ using ${\bar{\cal L}_{p_{q,y}}}\left({p_{y}}\right)$ and then ${{{{{p}}}_{q,x}^{(r)}}}$ with ${\bar{\cal L}_{{p_{q,x}}\left|{{p_{q,y}}}\right.}}\left({{p_{x}}\left|p_{y}={{{{{p}}}_{q,y}^{(r)}}}\right.}\right)$ is straightforward.

Recalling the expression of ${{\hat{p}}_{q,x}}$ and ${{\hat{p}}_{q,y}}$ in (21) and using the IS function defined in (30), the position estimator of the $q$-th emitter can be equivalent to estimate the linear mean [35] of generated realizations

where

However, the linear mean just averages all realizations without considering the existence of outlier seeds, which may result in the inevitable estimation bias. An improvement is using the circular mean [38] instead. As mentioned in [32], under the condition that ${\rho_{0}}$ is as high as desired, the circular mean estimation corresponds to the realization that minimizes the Euclidean distance to the true parameter. For a given random variable $Z\in\left[{-\pi,\pi}\right]$, let ${h}\left({Z}\right)$ be a transformation of $Z$, then the circular mean of ${h}\left({Z}\right)$ can be expressed as

where $\left\{{{{Z}^{(r)}}}\right\}_{r=1}^{R}$ are the realizations of random variable $Z$.

Assuming that the x-coordinates and y-coordinates of all emitters are in the range $\left[{p_{x,min},p_{x,max}}\right]$ and $\left[{p_{y,min},p_{y,max}}\right]$, respectively. Then the alternative formulation of the estimators shown in (21) using circular mean are expressed as

where

Note the normalization integral term in $\eta\left({{{{\bf{\bar{p}}}}^{(r)}}}\right)$ will not affect the result of $\angle\;\left\{\cdot\right\}$ in (39a) and (39b), which can be omitted to reduce the computational load. Actually, defining the weighting coefficient as

and replacing $\eta\left({{{{\bf{\bar{p}}}}^{(r)}}}\right)$ with $\eta^{\prime}\left({{{{\bf{\bar{p}}}}^{(r)}}}\right)$ in (39a) and (39b) can further reduce the computational consumption and guarantee the accuracy [28].

For the sake of clarity, we give all the necessary details for the location estimation of all $Q$ emitters by importance sampling as follows:

1. STEP 1:

   Derive separable term ${I_{q}}\left({{{\bf{p}}}}\right)$ by (29) for every $q=1,2,\ldots,Q$ with the complex signal vectors $\left\{{{{\bf{\bar{r}}}}_{l,k}}\right\}_{l=1,k=1}^{L,K}$ observed from all interception intervals.

2. STEP 2:

   Choose ${N_{x}}\times{N_{y}}$ grid points in coordinate range $\left[{p_{x,min},p_{x,max}}\right]\times\left[{p_{y,min},p_{y,max}}\right]$. All grid points $\left({{p_{{n_{x}}}},{p_{{n_{y}}}}}\right)$ have discretization steps $\Delta_{p_{x}}$ and $\Delta_{p_{y}}$ ($n_{x}\in\left\{{1,2,\ldots,N_{x}}\right\},n_{y}\in\left\{{1,2,\ldots,N_{y}}\right\}$). Then, by approximating integrals with discrete sums, the pseudo-PDF ${\bar{\cal L}_{q}}\left({\bf{p}}\right)$ for every $q$ can be approximated as

   $${\bar{\cal L}_{q}}\left({{p_{{n_{x}}}},{p_{{n_{y}}}}}\right)\approx\frac{{{e^{{\rho_{1}}{I_{q}}\left({{p_{{n_{x}}}},{p_{{n_{y}}}}}\right)}}}}{{\sum\nolimits_{{n_{x}}=1}^{{N_{x}}}{\sum\nolimits_{{n_{y}}=1}^{{N_{y}}}{{e^{{\rho_{1}}{I_{q}}\left({{p_{{n_{x}}}},{p_{{n_{y}}}}}\right)}}\Delta_{{p_{x}}}\Delta_{{p_{y}}}}}}}.$$

   (41)

3. STEP 3:

   Compute the marginal PDF ${\bar{\cal L}_{p_{q,x}}}\left({p_{x}}\right)$ from (41)

   $${\bar{\cal L}_{p_{q,x}}}\left({p_{{n_{x}}}}\right)={{\sum\nolimits_{{n_{y}}=1}^{{N_{y}}}{{\bar{\cal L}_{q}}\left({{p_{{n_{x}}}},{p_{{n_{y}}}}}\right)\Delta_{p_{y}}}}},$$

   (42)

   then compute the corresponding CDF

   $${\bar{\cal G}_{p_{q,x}}}\left({p_{{n_{x}}}}\right)={{\sum\nolimits_{{n^{\prime}}=1}^{{n_{x}}}{{\bar{\cal L}_{p_{q,x}}}\left({p_{{n^{\prime}}}}\right)\Delta_{p_{x}}}}}.$$

   (43)

4. STEP 4:

   Generate $R$ uniform random numbers $\left\{{{{u}^{(r)}_{q,x}}}\right\}_{r=1}^{R}\sim U\left[{0,1}\right]^{1}$ and perform linear interpolation to ${\bar{\cal G}_{p_{q,x}}}\left({p_{{n_{x}}}}\right)$ for the generation of $R$ position x-coordinate realizations

   $${p_{q,x}^{(r)}}=\mathop{\arg\min}\limits_{{p_{{x}}}}\left|{{\bar{\cal G}_{p_{q,x}}}\left({p_{{x}}}\right)-{{{u}^{(r)}_{q,x}}}}\right|.$$

   (44)

5. STEP 5:

   Compute the the conditional PDF of ${{p_{q,y}}}$ given ${{{{{p}}}_{q,x}^{(r)}}}$ with linear interpolation

   $${\bar{\cal L}_{{p_{q,y}}\left|{{p_{q,x}}}\right.}}\left({{p_{{n_{y}}}}\left|p_{x}={{{{{p}}}_{q,x}^{(r)}}}\right.}\right)=\frac{{\bar{\cal L}_{q}}\left({{{{{{p}}}_{q,x}^{(r)}}},{p_{{n_{y}}}}}\right)}{{\bar{\cal L}_{p_{q,x}}}\left({{{{{p}}}_{q,x}^{(r)}}}\right)}.$$

   (45)

6. STEP 6:

   Evaluate the CDF ${\bar{\cal G}_{{p_{q,y}}\left|{{p_{q,x}}}\right.}}\left({{p_{{n_{y}}}}\left|{{{{{p}}}_{q,x}^{(r)}}}\right.}\right)$ as similarly as (43) and generate uniform random numbers $\left\{{{{u}^{(r)}_{q,y}}}\right\}_{r=1}^{R}\sim U\left[{0,1}\right]^{1}$.

7. STEP 7:

   Similarly as (44), obtain position y-coordinate realizations $\left\{{p_{q,y}^{(r)}}\right\}_{r=1}^{R}$ with ${\bar{\cal G}_{{p_{q,y}}\left|{{p_{q,x}}}\right.}}\left({{p_{{y}}}\left|{{{{{p}}}_{q,x}^{(r)}}}\right.}\right)$, using linear interpolation technique as well.

8. STEP 8:

   Repeat STEP 2-STEP 7 until finishing the generation of $\left\{{{{{\bf{\bar{p}}}}^{(r)}}}\right\}_{r=1}^{R}$ for all $Q$ emitters.

9. STEP 9:

   Finally estimate the location for the $q$-th emitter