Reflector-aided direct locations of multiple signals in the presence of small reflector position biases

Abstract

Direct position determination (DPD) is a single-step method that localizes transmitters from sensor outputs without computing intermediate parameters. It outperforms conventional two-step localization methods, especially under low signal-to-noise ratio conditions. This article proposes a reflector-aided DPD algorithm for multiple signals of known waveforms received by an array observer. In previous studies, reflector-aided localization has always required very precise locations of reflectors. Therefore, the localization performance depends sensitively on accurately knowing each reflector position. This study considers the presence of small biases in reflector locations. To make the problem tractable, we simplify the signal model through an approximation using the first-order Taylor expansion and then directly localize multiple sources in a decoupled manner. Unlike most DPDs that presume noise is spatially uncorrelated, our study imposes no restriction on the correlation structure of noise, allowing this algorithm to be used in more general scenarios. In addition, we derive the Cramér–Rao bound expression and perform an analysis of the direct locations of multiple signals when the reflector positions are assumed accurate but in fact have small biases. Simulation results corroborate the theoretical results and a good localization performance of the proposed algorithm in the presence of small reflector position biases.

Introduction

The problem of passive localization attracts significant interest in various fields such as array signal processing, radar, sonar, seismology, and radio astronomy.^1,2 This study investigates the passive geolocation of multiple transmitters for a specific geolocation architecture called reflector-aided single-platform geolocation (SPG). Introduced by Bar-Shalom and Weiss,^3,4 this architecture has several advantages over the geolocation system comprising multiple receivers. For example, system deployment costs, including the amount of transferred data required, are significantly reduced. ⁴ The SPG system consists of a single receiver equipped with an antenna array and multiple reflectors that are placed at known locations. These reflectors reflect the transmitter’s signal toward the array receiver and therefore produce the signal multipath.

To address the localization problem in the multipath environment, many two-step processing methods have been developed,^5–10 where the measurement parameters (e.g. direction of arrival (DOA) and time of arrival (TOA)) of the multipath signals are first extracted and then the source positions are estimated using statistical or parametric methods. ⁵ Because they involve two separate estimation steps, they cannot guarantee high localization precision for the following two reasons:^11,12 (1) they estimate parameters at each observer separately and independently, as the constraint that the measurements correspond to the same transmitter location is ignored; and (2) data associations are necessary. If the measurements are not related to the correct path, over-modeling occurs, leading to spurious paths ¹³ and hence additional errors. Moreover, two-step processing methods are computationally demanding for locating multipath sources because the computational complexity of the intermediate parameter estimations increases exponentially with the number of paths.

In contrast to the conventional two-step methods, direct position determination (DPD) is a single-step localization method without needing to compute the intermediate parameter values and augments the position estimation with the constraint that all measurements correspond to the same geolocation of the transmitter. In the multipath environment based on a single-bounce channel model, a DPD algorithm was proposed to localize multiple sources when the positions of the main scatterers or reflectors are known. ¹³ It employs a stationary array to directly localize multipath signals by minimizing a multiple signal classification (MUSIC)–type function in terms of source positions without explicitly computing DOAs or TOAs. This method takes advantage of the multipath structure and is in fact reflector-aided, whereas it is suboptimal compared with the maximum likelihood (ML) estimators. Using the reflector-aided SPG system, Bar-Shalom and Weiss,^3,4 developed the ML-based DPD algorithms for a single transmitter, where the transmitted waveform is either known or unknown a priori. Their geolocation system consists of a single stationary or moving receiver, and several static, passive transponders acting as ideal signal reflectors. With accurate positions of the transponders, multipath components are incorporated to implement the localization aided by geometric information. As opposed to the two-step localization methods, the foregoing DPDs directly localize the transmitter from the sensor output by exploiting the location information contained in the multipath components. As a result, they achieved higher localization accuracy, especially under low signal-to-noise ratio (SNR) conditions. Despite this progress, these DPDs encounter two problems. They require very precise knowledge of the reflector locations, and hence, their location performance is very sensitive to the accuracy of each reflector position. However, in reality, the true reflector positions may be different from their nominal values for various reasons, and a slight error in the reflector positions ultimately generates extra errors in source location. Another limitation is that they are vulnerable to spatially colored noise since they presume that the noise is spatially uncorrelated, whereas the noise is frequently correlated along the array with an unknown correlation structure. ¹⁴

In certain applications, the transmitted waveforms are known to the receiver,^15,16 for example, the synchronization and training sequences are transmitted periodically and are known a priori in cellular systems. In this article, we consider the reflector-aided SPG system and concentrate on the localization of multiple signals with known waveforms and unknown complex attenuation under multipath propagation. Similar to previous studies,^3,4,13 the geometric relationship of reflectors, receiver, and transmitters is used to locate multipath signals based on a single-bounce channel model. More practically, we assume that the nominal positions of the reflectors are available, but small biases in the reflector positions exist for the indirect paths of each source.

Our contributions are summarized as follows:

1. We establish a simplified signal model using the first-order Taylor series approximation, and thereby separate the parameters of unknown biases from those of source locations.

2. Based on this simplified model, we propose a single-step location algorithm that applies the weighted least squares (WLS) criterion with the optimal weighting matrix, which is a realization of the ML estimator. This algorithm enjoys the following features:

3. It demonstrates high localization accuracy in the presence of small reflector position biases.

4. It places no restrictions on the model of noise and hence is effective in spatially uncorrelated noise, as well as spatially correlated noise with unknown covariance matrix.

5. It is computationally efficient as it makes use of the known waveforms to decouple the WLS-based localizations of multiple uncorrelated sources.

We evaluate the compact Cramér–Rao bound (CRB) for source position parameters based on the simplified model. The mean square error (MSE) of the locations is also derived when the WLS-based estimator assumes that the reflector positions are accurate but in fact have small biases.

Paper organization

The “Signal model and problem formulation” section establishes the signal model in the presence of small reflector position biases and formulates the problem. The “Decoupled location in the presence of small reflector position biases” section proposes the decoupled solution for the DPD for multiple sources in the reflector-aided SPG system. In the “CRB and MSE” section, we derive the CRB and the MSE without accounting for reflector position biases. In the “Numerical results” section, simulation results are presented forming the basis for a comparison of the proposed algorithm with other location methods. The “Conclusion” section concludes the study with a summary and remarks.

Signal model and problem formulation

In previous studies,^3,4,13 reflector-aided SPG methods require a very precise knowledge of the reflector positions. A slight bias in the reflector positions ultimately generates an extra error in source location. This section shows how to construct the signal model and formulate the DPD problem in the presence of small reflector position biases. We shall first present the multipath signal model and then simplify it assuming that the reflector position biases are small. The simplification is necessary because it separates the biases from the source positions and makes the problem tractable. Finally, we relate the DOAs and TOAs contained in the multipath components to the geometric relationships between reflectors, receiver, and transmitters, and thereby formulate our DPD problem.

Signal model in the presence of reflector position biases

Let us consider a stationary receiver equipped with an antenna array of M isotropic sensors intercepting the transmitted signals. Here, Q transmitters radiate narrowband signals (i.e. the bandwidth of each source signal is small enough compared with the inverse of the propagation time over the array aperture) in the far field of the array. The transmitters and receiver are assumed to be time synchronized. There are N_R reflectors within the area of interest, and the qth source signal arrives at the antenna array via L_q paths. As the single-bounce model is a simplified parametric model, ¹⁷ we consider the single-bounce reflection channel in our work (see Figure 1). Let p _q ( $q=1,2,\dots Q$ ) denotes the position of the qth transmitter, v denotes the position of the receiver, and u _n ( $n=1,2,\dots N_R$ ) denotes the nominal position of the nth reflector. For simplicity and without loss of generality, we consider planar geometry in this article.

Figure 1.

Scenario of the reflector-aided SPG based on a single-bounce channel model.

In this scenario, the observation r(t) containing complex envelopes of the signals collected by the receiver array at time 0 ≤ t ≤ T is expressed as

$$r\left(t\right)=\sum _{q=1}^Q\sum _{l=0}^{L_q-1}{\beta }_{\mathrm{ql}}\mathbf{a}\left({\theta }_{\mathrm{ql}}\right)s_q(t-{\tau }_{\mathrm{ql}})+\mathbf{n}\left(t\right)$$

(1)

where l = 0 refers to the direct path; $s_q(t-{\tau }_{\mathrm{ql}})$ is the complex envelope of the qth source signal delayed by τ_ql, which denotes the propagation time of the lth path of the qth transmitter; β_ql is an unknown complex attenuation coefficient representing the channel effect for the lth path of the qth source; θ_ql signifies the DOA of the lth path of the qth source; a(θ_ql) is the corresponding array response; and n(t) represents the Gaussian noise mixed through this array. The sources and noises are assumed to be uncorrelated and to have a mean of zero.

In reality, the reflector positions are estimated in advance and then used to localize the transmitters.^10,13 Localization inaccuracies of the reflectors produce true reflector positions that are different from the nominal values. Furthermore, even when different signals are reflected by the same reflector, the signals may be received from different positions corresponding to the main scatterers on this reflector. ¹⁰ We therefore assume that there are small biases in the nominal positions of reflectors, and the biases for the indirect paths of any two sources are different. With this assumption, the true reflector position for the lth indirect path of the qth source is expressed as

$${\tilde{\mathbf{u}}}_{\mathrm{ql}}={\mathbf{u}}_{l_q}+\Delta {\mathbf{u}}_{\mathrm{ql}}$$

(2)

where l_q denotes the sequence number of the reflector corresponding to the lth path of the qth source, and Δu _ql represents the unknown bias in the reflector position for this path and is assumed to be sufficiently small. In particular, when the lth path of the qth source and the $l\prime \mathrm{th}$ path of $q\prime \mathrm{th}$ source share the same reflector, we have ${\mathbf{u}}_{l_q}={\mathbf{u}}_{l{\prime }_q}\prime$ and $\Delta {\mathbf{u}}_{\mathrm{ql}}\ne \Delta {\mathbf{u}}_{q\prime l}\prime$ .

With the effect of the reflector position biases, the signal model in equation (1) becomes

$$r\left(t\right)={\displaystyle \sum _{q=1}^Q\left\{{\displaystyle \sum _{l=1}^{L_q-1}{\beta }_{ql}a\left({\tilde{\theta }}_{ql}\right)s_q\left(t-{\tilde{\tau }}_{ql}\right)})+{\beta }_{q0}a\left({\theta }_{q0}\right)s_q\left(t-{\mathrm{\tau }}_{q0}\right)\right\}}+n\left(t\right)$$

(3)

where ${\tilde{\theta }}_{\mathrm{ql}}$ and ${\tilde{\tau }}_{\mathrm{ql}}$ are the true DOA and TOA deviating from their nominal values, θ_ql and τ_ql. Therefore, we have

$$\mathbf{a}\left({\tilde{\theta }}_{\mathrm{ql}}\right)=\mathbf{a}({\theta }_{\mathrm{ql}}+\Delta {\theta }_{\mathrm{ql}})$$

(4)

and

$${\tilde{\tau }}_{\mathrm{ql}}={\tau }_{\mathrm{ql}}+\Delta {\tau }_{\mathrm{ql}}$$

(5)

in which Δθ_ql and Δτ_ql signify the corresponding offsets contributed by the reflector position biases.

Assuming that small reflector position biases result in small TOA offsets, Δτ_ql, can be included in the phase shift. Hence, equation (3) is rewritten as

$$r\left(t\right)={\displaystyle \sum _{q=1}^Q\left\{{\displaystyle \sum _{l=1}^{L_q-1}{\beta }_{ql}{e}^{-\mathrm{j}2\mathrm{\pi }{f}_c\Delta {\tau }_{ql}}a\left({\tilde{\theta }}_{ql}\right)s_q\left(t-{\tau }_{ql}\right)}+{\beta }_{q0}a\left({\theta }_{q0}\right)s_q\left(t-{\tau }_{q0}\right)\right\}}+n\left(t\right)$$

(6)

where f_c is the carrier frequency of the transmitted signals.

Next, we partition the observed signal into K sections each of length T/K. Sampling the signals in each section at time intervals $t=j{T}_s(j=1,2,\dots ,J)$ , where $T_s$ is the sampling period, the samples in the kth section are given by

$$\begin{array}{l}r\left(j,k\right)={\displaystyle \sum _{q=1}^Q\left\{{\displaystyle \sum _{l=1}^{L_q-1}{\beta }_{ql}{e}^{-\mathrm{j}2\mathrm{\pi }{f}_c\Delta {\tau }_{ql}}a\left({\tilde{\theta }}_{ql}\right)s_q\left(j{T}_s-{\tau }_{ql},k\right)}+{\beta }_{q0}a\left({\theta }_{q0}\right)s_q\left(j{T}_s-{\tau }_{q0},k\right)\right\}}+n\left(j,k\right)\\ k=1,2,\dots ,K,j=1,2,\dots ,J,\end{array}$$

(7)

where $s_q(j{T}_s-{\tau }_{\mathrm{ql}},k)$ and $\mathbf{n}(j,k)$ are the jth samples of s_q(t − τ_ql) and n(t) in the kth section.

Let us define ${\overline{\mathbf{a}}}_0(j,{\mathbf{p}}_q)=\mathbf{a}({\theta }_{q0})e^{-\mathrm{j}2\mathrm{\pi }j{\tau }_{q0}/J{T}_s}$ and ${\tilde{\mathbf{a}}}_l(j,{\mathbf{p}}_q)=\mathbf{a}({\tilde{\theta }}_{\mathrm{ql}})e^{-\mathrm{j}2\mathrm{\pi }(j{\tau }_{\mathrm{ql}}/J{T}_s+f_c\Delta {\tau }_{\mathrm{ql}})}$ . Note that, although $\mathbf{a}({\theta }_{q0})e^{-\mathrm{j}2\mathrm{\pi }j{\tau }_{q0}/J{T}_s}$ and $\mathbf{a}({\tilde{\theta }}_{\mathrm{ql}})e^{-\mathrm{j}2\mathrm{\pi }(j{\tau }_{\mathrm{ql}}/J{T}_s+f_c\Delta {\tau }_{\mathrm{ql}})}$ are not explicitly expressed in terms of ${\mathbf{p}}_q$ , ${\mathbf{p}}_q$ is implicit in the DOAs and TOAs. Therefore, these terms can be parameterized by the position of the qth source, and the relationships between the DOAs and TOAs and the source positions are established later. Assuming that T/K >> max{τ_ql} and the signals are stationary, the discrete Fourier transform (DFT) of r(j, k) can be expressed as

$$\overline{\mathbf{r}}(j,k)=\sum _{q=1}^Q\{\sum _{l=1}^{L_q-1}{\beta }_{\mathrm{ql}}{\tilde{\mathbf{a}}}_l(j,{\mathbf{p}}_q)+{\beta }_{q0}{\overline{\mathbf{a}}}_0(j,{\mathbf{p}}_q)\}{\overline{s}}_q(j,k)+\overline{\mathbf{n}}(j,k)k=1,2,\dots ,K,j=1,2,\dots ,J$$

(8)

where ${\overline{s}}_q(j,k)$ and $\overline{\mathbf{n}}(j,k)$ signify the $j\mathrm{th}$ DFT coefficients of $s_q(j{T}_s,k)$ and $\mathbf{n}(j,k)$ , respectively.

From equation (8), it is clear that the source positions, complex attenuation coefficients, and the reflector position biases have a nonlinear relationship in the observations, making the localization problem complicated. Furthermore, identifying all the unknown parameters is hard as the reflector position biases add to the numerous unknowns. To eliminate this difficulty in the ensuing, we apply an efficient approximation to simplify the signal model.

Simplification on signal model

Assuming small reflector position biases, we now come to simplify our signal model. As $\mathbf{a}({\tilde{\theta }}_{\mathrm{ql}})$ is nonlinearly related to the reflector position biases, we expand $\mathbf{a}({\tilde{\theta }}_{\mathrm{ql}})$ in a first-order Taylor series expansion

$$\mathbf{a}\left({\tilde{\theta }}_{\mathrm{ql}}\right)\approx \mathbf{a}\left({\theta }_{\mathrm{ql}}\right)+{\mathbf{d}}_{\mathrm{ql}}{\mathbf{g}}_{\mathrm{ql}}^{\mathrm{T}}\Delta {\mathbf{u}}_{\mathrm{ql}}$$

(9)

in which ${\mathbf{d}}_{\mathrm{ql}}={\frac{\partial \mathbf{a}(\theta )}{\partial \theta }|}_{\theta ={\theta }_{\mathrm{ql}}}$ denotes the first-order derivative of $\mathbf{a}(\theta )$ with $\theta ={\theta }_{\mathrm{ql}}$ , and ${\mathbf{g}}_{\mathrm{ql}}={\frac{\partial {\theta }_{\mathrm{ql}}}{\partial \mathbf{u}}|}_{\mathbf{u}={\mathbf{u}}_{l_q}}$ represents the first-order partial derivative of ${\theta }_{\mathrm{ql}}$ with $\mathbf{u}={\mathbf{u}}_{l_q}$ . Using equation (9) yields an approximation to ${\tilde{\mathbf{a}}}_l(j,{\mathbf{p}}_q)=\mathbf{a}({\tilde{\theta }}_{\mathrm{ql}})e^{-\mathrm{j}2\mathrm{\pi }(j{\tau }_{\mathrm{ql}}/J{T}_s+f_c\Delta {\tau }_{\mathrm{ql}})}$

$$\begin{array}{c}{\tilde{\mathbf{a}}}_l(j,{\mathbf{p}}_q)\approx \mathbf{a}\left({\theta }_{\mathrm{ql}}\right)e^{-\mathrm{j}2\mathrm{\pi }(j{\tau }_{\mathrm{ql}}/J{T}_s+f_c\Delta {\tau }_{\mathrm{ql}})}\hfill \\ \multicolumn{1}{c}{+{\mathbf{d}}_{\mathrm{ql}}{\mathbf{g}}_{\mathrm{ql}}^{\mathrm{T}}\Delta {\mathbf{u}}_{\mathrm{ql}}{e}^{-\mathrm{j}2\mathrm{\pi }(j{\tau }_{\mathrm{ql}}/J{T}_s+f_c\Delta {\tau }_{\mathrm{ql}})}}\\ \approx ({\overline{\mathbf{a}}}_l(j,{\mathbf{p}}_q)+{\gamma }_{\mathrm{ql}}{\overline{\mathbf{d}}}_l(j,{\mathbf{p}}_q))e^{-\mathrm{j}2\mathrm{\pi }{f}_c\Delta {\tau }_{\mathrm{ql}}}\hfill \end{array}$$

(10)

for $l=1,2,\dots ,L_q-1$ , where ${\gamma }_{\mathrm{ql}}={\mathbf{g}}_{\mathrm{ql}}^{\mathrm{T}}\Delta {\mathbf{u}}_{\mathrm{ql}}$ , ${\overline{\mathbf{a}}}_l(j,{\mathbf{p}}_q)=\mathbf{a}({\theta }_{\mathrm{ql}})e^{-\mathrm{j}2\mathrm{\pi }j{\tau }_{\mathrm{ql}}/J{T}_s}$ , and ${\overline{\mathbf{d}}}_l(j,{\mathbf{p}}_q)={\mathbf{d}}_{\mathrm{ql}}{e}^{-\mathrm{j}2\mathrm{\pi }j{\tau }_{\mathrm{ql}}/J{T}_s}$ . Note that ${\overline{\mathbf{a}}}_l(j,{\mathbf{p}}_q)$ and ${\overline{\mathbf{d}}}_l(j,{\mathbf{p}}_q)$ can be parameterized by ${\mathbf{p}}_q$ because ${\tau }_{\mathrm{ql}}$ in $\mathbf{a}({\theta }_{\mathrm{ql}})e^{-\mathrm{j}2\mathrm{\pi }j{\tau }_{\mathrm{ql}}/J{T}_s}$ and ${\mathbf{d}}_{\mathrm{ql}}{e}^{-\mathrm{j}2\mathrm{\pi }j{\tau }_{\mathrm{ql}}/J{T}_s}$ are related to ${\mathbf{p}}_q$ .

Substituting equation (10) into equation (8), we can approximate the frequency-domain sample as

$$\begin{array}{c}\overline{\mathbf{r}}(j,k)\approx \sum _{q=1}^Q\{\sum _{l=1}^{L_q-1}{\beta }_{\mathrm{ql}}{e}^{-\mathrm{j}2\mathrm{\pi }{f}_c\Delta {\tau }_{\mathrm{ql}}}({\overline{\mathbf{a}}}_l(j,{\mathbf{p}}_q)+{\gamma }_{\mathrm{ql}}{\overline{\mathbf{d}}}_l(j,{\mathbf{p}}_q))+{\beta }_{q0}{\overline{\mathbf{a}}}_0(j,{\mathbf{p}}_q)\}{\overline{s}}_q(j,k)+\overline{\mathbf{n}}(j,k)\approx \sum _{q=1}^Q{\mathbf{A}}_q(j,{\mathbf{p}}_q)\mathit{\beta }{\prime }_q{\overline{s}}_q(j,k)+\overline{\mathbf{n}}(j,k)k=1,2,\dots ,K,j=1,2,\dots ,J\hfill \end{array}$$

(11)

in which

$$\begin{array}{c}{\mathbf{A}}_q(j,{\mathbf{p}}_q)=[{\overline{\mathbf{a}}}_0(j,{\mathbf{p}}_q),\dots ,{\overline{\mathbf{a}}}_{L_q-1}(j,{\mathbf{p}}_q),{\overline{\mathbf{d}}}_1(j,{\mathbf{p}}_q),\dots ,{\overline{\mathbf{d}}}_{L_q-1}(j,{\mathbf{p}}_q)]\in {\mathbb{C}}^{M\times (2L_q-1)}\hfill \\ \mathit{\beta }{\prime }_q={[{\beta }_{q0},\beta {\prime }_{q1},\dots ,\beta {\prime }_{q(L_q-1)},{\gamma }_{q1}\beta {\prime }_{q1},\dots ,{\gamma }_{q(L_q-1)}\beta {\prime }_{q(L_q-1)}]}^{\mathrm{T}}\in {\mathbb{C}}^{(2L_q-1)\times 1}\hfill \end{array}$$

(12)

with $\beta {\prime }_{\mathrm{ql}}={\beta }_{\mathrm{ql}}{e}^{-\mathrm{j}2\mathrm{\pi }{f}_c\Delta {\tau }_{\mathrm{ql}}}(l=1,2,\dots ,L_q-1)$

The expression in equation (11) can be written in a more compact form

$$\overline{\mathbf{r}}(j,k)=\mathbf{\Phi }\left(j\right)\mathbf{\Gamma }\overline{\mathbf{s}}(j,k)+\overline{\mathbf{n}}(j,k)k=1,2,\dots ,K,j=1,2,\dots ,J$$

(13)

where

$$\begin{array}{c}\overline{\mathbf{s}}(j,k)={[{\overline{s}}_1(j,k),{\overline{s}}_2(j,k),\dots ,{\overline{s}}_Q(j,k)]}^{\mathrm{T}}\in {\mathbb{C}}^{Q\times 1}\hfill \\ \mathbf{\Phi }\left(j\right)=[{\mathbf{A}}_1(j,{\mathbf{p}}_1),{\mathbf{A}}_2(j,{\mathbf{p}}_2),\dots ,{\mathbf{A}}_Q(j,{\mathbf{p}}_Q)]\in {\mathbb{C}}^{M\times L}\hfill \\ \mathbf{\Gamma }=\mathrm{blkdiag}\{\mathit{\beta }{\prime }_1,\mathit{\beta }{\prime }_2,\dots ,\mathit{\beta }{\prime }_Q\}\in {\mathbb{C}}^{L\times Q}\hfill \end{array}$$

(14)

with $L=2\sum _{q=1}^Q{L}_q-Q$ .

Remark 1

Note that the unknown reflector position biases are incorporated into the complex attenuation coefficients for the indirect paths (i.e. $\beta {\prime }_{\mathrm{ql}}$ and ${\gamma }_{\mathrm{ql}}\beta {\prime }_{\mathrm{ql}}$ for $l=1,2,\dots ,L_q-1$ ). Based on the above approximation, we separate the unknown biases from the source positions, thereby making the problem tractable. As ${\mathbf{A}}_q(j,{\mathbf{p}}_q)$ is related to only the nominal positions of the reflectors, we then exploit the information contained in the multipath components without the effect of the unknown biases.

Problem formulation

As our work focuses on directly locating the multiple signals in the reflector-aided SPG system, we now express the nominal DOAs and TOAs in terms of the locations of the transmitters and receiver and the nominal locations of the reflectors

$$\{\begin{array}{c}{\theta }_{q0}=\arctan \frac{{\mathbf{p}}_q\left(1\right)-\mathbf{v}\left(1\right)}{{\mathbf{p}}_q\left(2\right)-\mathbf{v}\left(2\right)}\hfill \\ \multicolumn{1}{c}{{\theta }_{\mathrm{ql}}=\arctan \frac{{\mathbf{u}}_{l_q}\left(1\right)-\mathbf{v}\left(1\right)}{{\mathbf{u}}_{l_q}\left(2\right)-\mathbf{v}\left(2\right)}l=1,2,\dots ,L_q-1\forall q}\\ \multicolumn{1}{c}{{\tau }_{q0}=\frac{1}{c}\Vert {\mathbf{p}}_q-\mathbf{v}\Vert ,}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{ql}}=\frac{1}{c}(\Vert {\mathbf{p}}_q-{\mathbf{u}}_{l_q}\Vert +\Vert {\mathbf{u}}_{l_q}-\mathbf{v}\Vert )l=1,2,\dots ,L_q-1}\end{array}$$

(15)

Furthermore, the following assumptions are adopted:

• A1. The nominal positions and total number of the reflectors, N_R, are known, whereas the number of paths for each source, L_q, is unknown.

• A2. The signal waveforms are known and uncorrelated.

• A3. The noise vector $\overline{\mathbf{n}}(j,k)$ is a circularly complex Gaussian process with zero mean. It is temporally uncorrelated with its second-order moments given by

$$\{\begin{array}{c}\mathbb{E}\left[\overline{\mathbf{n}}(j,k){\overline{\mathbf{n}}}^{\mathrm{H}}(i,m)\right]={\mathbf{R}}_{\mathrm{n}}{\delta }_{j,i}{\delta }_{k,m}\hfill \\ \multicolumn{1}{c}{\mathbb{E}\left[\overline{\mathbf{n}}(j,k){\overline{\mathbf{n}}}^{\mathrm{T}}(i,m)\right]=\mathbf{O}}\end{array}$$

(16)

where ${\mathbf{R}}_n$ is assumed to be positive definite but is otherwise unknown. ${\delta }_{j,i}$ represents the Kronecker Delta function, namely, if $i=j$ , ${\delta }_{j,i}=1$ , and otherwise ${\delta }_{j,i}=0$ .

Given the above assumptions, the problem addressed is to estimate the source positions given the frequency-domain samples of the K sections (see equation (13)). Different from the conventional two-step methods, we attempt to locate multiple transmitters in a single step using the geometric relation (see equation (15)) without explicitly computing TOAs and DOAs.

Decoupled location in the presence of small reflector position biases

We now have established the approximated signal model that accounts for small reflector position biases. In this section, as the transmitted waveforms are known to the location system in certain applications,^15,16,18 we follow the idea of Weiss and Amar ¹⁵ to decouple the locations of multiple sources assisted by the known waveforms. Note that the unknown noise structure is considered in our study.

WLS-based function

For notational convenience, we define $\tilde{\mathbf{p}}=[{\mathbf{p}}_1^{\mathrm{T}},{\mathbf{p}}_2^{\mathrm{T}},\dots ,{\mathbf{p}}_Q^{\mathrm{T}}{]}^{\mathrm{T}}$ and $\tilde{\mathit{\beta }}=[{\mathit{\beta }\prime }_1^{\mathrm{T}},{\mathit{\beta }\prime }_2^{\mathrm{T}},\dots ,{\mathit{\beta }\prime }_Q^{\mathrm{T}}{]}^{\mathrm{T}}$ . Then, the unknown parameters are estimated by minimizing the WLS-based function

$$V(\tilde{\mathbf{p}},\tilde{\mathit{\beta }})=\sum _{j=1}^J\sum _{k=1}^K{(\overline{\mathbf{r}}(j,k)-\mathbf{B}\left(j\right)\overline{\mathbf{s}}(j,k))}^{\mathrm{H}}\mathbf{W}(\overline{\mathbf{r}}(j,k)-\mathbf{B}\left(j\right)\overline{\mathbf{s}}(j,k))$$

(17)

where

$$\mathbf{B}\left(j\right)=\mathbf{\Phi }\left(j\right)\mathbf{\Gamma }$$

(18)

depends on $\tilde{\mathbf{p}}$ and $\tilde{\mathit{\beta }}$ . $\mathbf{W}$ is a positive definite weighting matrix and is set to

$${\mathbf{W}}^{\mathrm{opt}}={\mathbf{R}}_n^{-1}$$

(19)

as the optimal weighting matrix. As R _n is assumed to be unknown, we use its consistent estimation computed from the sample data. This does not influence the asymptotic performance. ¹⁹ It is straightforward to estimate R _n from

$${\hat{\mathbf{R}}}_n=\frac{1}{\mathrm{KJ}}\sum _{j=1}^J\sum _{k=1}^K(\overline{\mathbf{r}}(j,k)-\hat{\mathbf{B}}\left(j\right)\overline{\mathbf{s}}(j,k)){(\overline{\mathbf{r}}(j,k)-\hat{\mathbf{B}}\left(j\right)\overline{\mathbf{s}}(j,k))}^{\mathrm{H}}$$

(20)

where $\hat{\mathbf{B}}(j)$ is the unstructured estimate of $\mathbf{B}(j)$ . An extension to the derivation in Li et al. ¹⁸ yields

$$\hat{\mathbf{B}}\left(j\right)={\hat{\mathbf{R}}}_{\mathrm{sr}}^{\mathrm{H}}\left(j\right){\hat{\mathbf{R}}}_{\mathrm{ss}}^{-1}\left(j\right)$$

(21)

Here, ${\hat{\mathbf{R}}}_{\mathrm{ss}}(j)=\frac{1}{K}\sum _{k=1}^K\overline{\mathbf{s}}(j,k){\overline{\mathbf{s}}}^{\mathrm{H}}(j,k)$ and ${\hat{\mathbf{R}}}_{\mathrm{sr}}(j)=\frac{1}{K}\sum _{k=1}^K\overline{\mathbf{s}}(j,k){\overline{\mathbf{r}}}^{\mathrm{H}}(j,k)$ . Substituting equations (21) into (20), we have

$${\hat{\mathbf{R}}}_{\mathrm{n}}=\frac{1}{J}\sum _{j=1}^J\{{\hat{\mathbf{R}}}_{\mathrm{rr}}\left(j\right)-{\hat{\mathbf{R}}}_{\mathrm{sr}}^{\mathrm{H}}\left(j\right){\hat{\mathbf{R}}}_{\mathrm{ss}}^{-1}\left(j\right){\hat{\mathbf{R}}}_{\mathrm{sr}}\left(j\right)\}$$

(22)

in which ${\hat{\mathbf{R}}}_{\mathrm{rr}}(j)=\frac{1}{K}\sum _{k=1}^K\overline{\mathbf{r}}(j,k){\overline{\mathbf{r}}}^{\mathrm{H}}(j,k)$ . We now replace the weighting matrix W with ${\hat{\mathbf{R}}}_n^{-1}$ in equation (17) and obtain the WLS-based function with the optimal weighting matrix

$$V(\tilde{\mathbf{p}},\tilde{\mathit{\beta }})=\sum _{j=1}^J\sum _{k=1}^K{(\overline{\mathbf{r}}(j,k)-\mathbf{B}\left(j\right)\overline{\mathbf{s}}(j,k))}^{\mathrm{H}}{\hat{\mathbf{R}}}_n^{-1}(\overline{\mathbf{r}}(j,k)-\mathbf{B}\left(j\right)\overline{\mathbf{s}}(j,k))$$

(23)

Decoupled solution for multiple sources

Solving equation (23) straightforwardly requires a high-dimensional search that is often impractical, and therefore, we exploit the information of known waveforms to ease the computation. Following the algebraic steps similar to those in Weiss and Amar, ¹⁵ the objective function in equation (23) is proportional to

$$V\prime (\tilde{\mathbf{p}},\tilde{\mathit{\beta }})=\mathrm{tr}\left\{\sum _{j=1}^J{(\mathbf{B}\left(j\right)-\hat{\mathbf{B}}\left(j\right))}^{\mathrm{H}}{\hat{\mathbf{R}}}_{\mathrm{n}}^{-1}(\mathbf{B}\left(j\right)-\hat{\mathbf{B}}\left(j\right)){\hat{\mathbf{R}}}_{\mathrm{ss}}\left(j\right)\right\}$$

(24)

Given that ${\hat{\mathbf{R}}}_{\mathrm{ss}}(j)$ is asymptotically diagonal as the signals are uncorrelated, we rewrite equation (24) as

$$\begin{array}{l}{V}^{\prime }\left(\tilde{p},\tilde{\mathit{\beta }}\right)=\mathrm{tr}\left\{{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{q=1}^Q{\left[{\widehat{R}}_{ss}(j)\right]}_{qq}{\left(b_q\left(j\right)-{\widehat{b}}_q\left(j\right)\right)}^{\mathrm{H}}{\widehat{R}}_{\mathrm{n}}^{-1}\left(b_q\left(j\right)-{\widehat{b}}_q\left(j\right)\right)}}\right\}\\ ={\displaystyle \sum _{q=1}^Q\left\{{\displaystyle \sum _{j=1}^J{\left(b_q\left(j\right)-{\widehat{b}}_q\left(j\right)\right)}^{\mathrm{H}}{\left[{\widehat{R}}_{ss}\left(j\right)\right]}_{qq}{\widehat{R}}_{\mathrm{n}}^{-1}\left(b_q\left(j\right)-{\widehat{b}}_q\left(j\right)\right)}\right\}}\end{array}$$

(25)

where ${\mathbf{b}}_q(j)$ and ${\hat{\mathbf{b}}}_q(j)$ are the $q\mathrm{th}$ columns of $\mathbf{B}(j)$ and $\hat{\mathbf{B}}(j)$ , respectively. From equation (18), ${\mathbf{b}}_q(j)$ is expressed as

$${\mathbf{b}}_q\left(j\right)={\mathbf{A}}_q(j,{\mathbf{p}}_q)\mathit{\beta }{\prime }_q$$

(26)

With equation (25), we define another weighting matrix ${\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}(j)=[{\hat{\mathbf{R}}}_{\mathrm{ss}}(j){]}_{\mathrm{qq}}{\hat{\mathbf{R}}}_{\mathrm{n}}^{-1}$ for the qth signal, which is proportional to the SNR with respect to the jth frequency component of the qth signal. Therefore, upon substituting equations (26) into (25), the locations of the Q transmitters can be decoupled into Q lower dimensional minimization problems

$$\{{\hat{\mathbf{p}}}_q,{\widehat{\mathit{\beta }\prime }}_q\}=\underset{\mathbf{p},\mathit{\beta }\mathbf{\prime }}{\mathrm{argmin}}{V}_q(\mathbf{p},\mathit{\beta }\prime )q=1,2,\dots ,Q$$

(27)

with

$$V_q(\mathbf{p},\mathit{\beta }\prime )=\sum _{j=1}^J{({\mathbf{A}}_q(j,\mathbf{p})\mathit{\beta }\prime -{\hat{\mathbf{b}}}_q\left(j\right))}^{\mathrm{H}}{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}\left(j\right)({\mathbf{A}}_q(j,\mathbf{p})\mathit{\beta }\prime -{\hat{\mathbf{b}}}_q\left(j\right))$$

(28)

Moreover, we concatenate all the frequency components and rewrite equation (28) in the more compact form

$$V_q(\mathbf{p},\mathit{\beta }\prime )={({\mathbf{A}}_q\left(\mathbf{p}\right)\mathit{\beta }\prime -{\hat{\mathbf{b}}}_q)}^{\mathrm{H}}{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}({\mathbf{A}}_q\left(\mathbf{p}\right)\mathit{\beta }\prime -{\hat{\mathbf{b}}}_q)$$

(29)

where

$$\begin{array}{c}{\mathbf{A}}_q\left(\mathbf{p}\right)={[{\mathbf{A}}_q^{\mathrm{T}}(1,\mathbf{p}),{\mathbf{A}}_q^{\mathrm{T}}(2,\mathbf{p}),\dots ,{\mathbf{A}}_q^{\mathrm{T}}(J,\mathbf{p})]}^{\mathrm{T}}\hfill \\ {\hat{\mathbf{b}}}_q={[{\hat{\mathbf{b}}}_q^{\mathrm{T}}\left(1\right),{\hat{\mathbf{b}}}_q^{\mathrm{T}}\left(2\right),\dots ,{\hat{\mathbf{b}}}_q^{\mathrm{T}}\left(J\right)]}^{\mathrm{T}}\hfill \\ {\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}=\mathrm{blkdiag}\{{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}\left(1\right),{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}\left(2\right),\dots ,{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}\left(J\right)\}\hfill \end{array}$$

(30)

Thus, the solution of $\mathit{\beta }{\prime }_q$ can be first derived by minimizing equation (29)

$$\hat{\mathit{\beta }}{\prime }_q={\left\{{\mathbf{A}}_q^{\mathrm{H}}\left(\mathbf{p}\right){\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}{\mathbf{A}}_q\left(\mathbf{p}\right)\right\}}^{-1}{\mathbf{A}}_q^{\mathrm{H}}\left(\mathbf{p}\right){\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}{\hat{\mathbf{b}}}_q$$

(31)

Substituting equation (31) back into equation (29) yields the following function in terms of only the source position

$$V_q\left(\mathbf{p}\right)={\hat{\mathbf{b}}}_q^{\mathrm{H}}{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\mathbf{\Pi }}_{{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\mathbf{A}}_q(\mathbf{p})}^{\perp }{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\hat{\mathbf{b}}}_q$$

(32)

Consequently, we establish the optimization problem to locate the $q\mathrm{th}$ source

$${\hat{\mathbf{p}}}_q=\underset{\mathbf{p}}{\mathrm{argmin}}{V}_q\left(\mathbf{p}\right)q=1,2,\dots ,Q$$

(33)

In this way, each source position can be determined via a two-dimensional search within the area of interest. Note that it is impossible to construct ${\mathbf{A}}_q(j,\mathbf{p})\in {\mathbb{C}}^{M\times (2L_q-1)}$ as L_q is unknown. To overcome this difficulty, one straightforward solution is to replace ${\mathbf{A}}_q(j,\mathbf{p})$ with

$${\mathbf{A}}_{\mathrm{R}}(j,\mathbf{p})=[{\overline{\mathbf{a}}}_0(j,\mathbf{p}),\dots ,{\overline{\mathbf{a}}}_{N_{\mathrm{R}}}(j,\mathbf{p}),{\overline{\mathbf{d}}}_1(j,\mathbf{p}),\dots ,{\overline{\mathbf{d}}}_{N_{\mathrm{R}}}(j,\mathbf{p})]\in {\mathbb{C}}^{M\times (2N_{\mathrm{R}}+1)}$$

(34)

After ${\hat{\mathbf{p}}}_q$ is estimated, we obtain $\hat{\mathit{\beta }}{\prime }_q$ by substituting $\mathbf{p}={\hat{\mathbf{p}}}_q$ into equation (31)

$$\hat{\mathit{\beta }}{\prime }_q={\left\{{\mathbf{A}}_{\mathrm{R}}^{\mathrm{H}}\left({\hat{\mathbf{p}}}_q\right){\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}{\mathbf{A}}_{\mathrm{R}}\left({\hat{\mathbf{p}}}_q\right)\right\}}^{-1}{\mathbf{A}}_{\mathrm{R}}^{\mathrm{H}}\left({\hat{\mathbf{p}}}_q\right){\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}{\hat{\mathbf{b}}}_q$$

(35)

where

$${\mathbf{A}}_{\mathrm{R}}\left(\mathbf{p}\right)={[{\mathbf{A}}_{\mathrm{R}}^{\mathrm{T}}(1,\mathbf{p}),{\mathbf{A}}_{\mathrm{R}}^{\mathrm{T}}(2,\mathbf{p}),\dots ,{\mathbf{A}}_{\mathrm{R}}^{\mathrm{T}}(J,\mathbf{p})]}^{\mathrm{T}}$$

(36)

The elements of $\hat{\mathit{\beta }}{\prime }_q$ corresponding to the spurious paths approach zeros, and therefore, the over-modeling problem is solved intuitively without prior knowledge of the path number.

To summarize this decoupled DPD, we provide the procedure displayed in Algorithm 1.

Algorithm 1.

Decoupled DPD Method.

1: Choose a suitable $K\in {\mathbb{Z}}^{+}$ , and $J\in {\mathbb{Z}}^{+}$ ( ${\mathbb{Z}}^{+}$ denotes the set of positive integers)

2: forj← 1, Jdo

3:  Compute ${\hat{\mathbf{R}}}_{\mathrm{ss}}(j)$ , ${\hat{\mathbf{R}}}_{\mathrm{sr}}(j)$ , and ${\hat{\mathbf{R}}}_{\mathrm{rr}}(j)$

4: end for

5: Estimate $\hat{\mathbf{B}}(j)$ and ${\hat{\mathbf{R}}}_{\mathrm{n}}$ according to equations (21) and (22)

6: for q← 1, Qdo

7:  Define the area of interest and a suitable grid of locations ${\mathbf{p}}_q^1,{\mathbf{p}}_q^2,\dots ,{\mathbf{p}}_q^{N_{\mathrm{G}}}$ , where $N_{\mathrm{G}}$ is the number of grids

8:  forn← 1, NGdo

9:  Evaluate ${\mathbf{A}}_{\mathrm{R}}({\mathbf{p}}_q^n)$ according to equation (36)

10:  Compute $V_q({\mathbf{p}}_q^n)$ according to equation (32) with ${\mathbf{A}}_q(\mathbf{p})$ replaced by ${\mathbf{A}}_{\mathrm{R}}({\mathbf{p}}_q^n)$

11: end for

12: Find the grid point for which $V_q$ is the minimum, and this point is the source location ${\hat{\mathbf{p}}}_q$

13: end for

CRB and MSE

To provide a benchmark for the highest localization accuracy that can be achieved, this section first evaluates the CRB regarding only the source position parameters based on the simplified signal model. Previous studies required very precise knowledge of the reflector positions. It is therefore of interest to examine the sensitivity of the localization precision to the biases of the reflector positions. To this end, we derive the MSE of the WLS-based estimator assuming that the reflector positions are accurate, but in fact the reflector positions have small biases.

CRB

To derive a compact CRB formula for only the source position parameters, we first state a proposition that has been proved in Pekka et al. ²⁰ This proposition is crucial in evaluating the position-related block of CRB.

Proposition 1

Assume that $\mathbf{CRB}(\mathbf{\omega })$ is the CRB matrix for the real vector $\mathbf{\omega }$ and define a real vector as $\mathbf{\omega }\prime =\mathbf{H}\mathbf{\omega }$ , where H is an invertible matrix. Then, the CRB matrix for vector $\mathbf{\omega }\prime$ is given by $\mathbf{CRB}(\mathbf{\omega }\prime )=\mathbf{H}\mathbf{CRB}(\mathbf{\omega }){\mathbf{H}}^{\mathrm{T}}$ .

Recall the multipath signal model in the presence of small reflector position biases based on the frequency domain presented earlier

$$\overline{\mathbf{r}}(j,k)=\sum _{q=1}^Q{\mathbf{A}}_q(j,{\mathbf{p}}_q)\mathit{\beta }{\prime }_q{\overline{s}}_q(j,k)+\overline{\mathbf{n}}(j,k)k=1,2,\dots ,K,j=1,2,\dots ,J$$

(37)

Relying on this model, an augmented parameter vector $\mathbf{\rho }$ that includes the full parameter vector is defined as

$$\mathbf{\rho }={[{\mathbf{\sigma }}^{\mathrm{T}},{\tilde{\mathbf{\omega }}}^{\mathrm{T}}]}^{\mathrm{T}}={[{\mathbf{\sigma }}^{\mathrm{T}},{\tilde{\mathbf{p}}}^{\mathrm{T}},\mathrm{R}{\mathrm{e}}^{\mathrm{T}}\{\tilde{\mathit{\beta }}\},{\mathrm{Im}}^{\mathrm{T}}\left\{\tilde{\mathit{\beta }}\right\}]}^{\mathrm{T}}$$

(38)

where $\mathbf{\sigma }$ is parameterized by the real parameters in ${\mathbf{R}}_n$ , and

$$\tilde{\mathbf{\omega }}={[{\tilde{\mathbf{p}}}^{\mathrm{T}},\mathrm{R}{\mathrm{e}}^{\mathrm{T}}\left\{\tilde{\mathit{\beta }}\right\},{\mathrm{Im}}^{\mathrm{T}}\left\{\tilde{\mathit{\beta }}\right\}]}^{\mathrm{T}}$$

(39)

According to Li et al., ¹² the CRB matrix for parameter $\tilde{\mathbf{\omega }}$ is computed from

$$\mathbf{CRB}\left(\tilde{\mathbf{\omega }}\right)=\frac{1}{2}{\mathrm{Re}}^{-1}\left\{\sum _{j=1}^J\sum _{k=1}^K\left\{{\left(\frac{\partial \sum _{q=1}^Q{\mathbf{A}}_q(j,{\mathbf{p}}_q)\mathit{\beta }{\prime }_q{\overline{s}}_q(j,k)}{\partial {\tilde{\mathbf{\omega }}}^{\mathrm{T}}}\right)}^{\mathrm{H}}{\mathbf{R}}_n^{-1}\left(\frac{\partial \sum _{q=1}^Q{\mathbf{A}}_q(j,{\mathbf{p}}_q)\mathit{\beta }{\prime }_q{\overline{s}}_q(j,k)}{\partial {\tilde{\mathbf{\omega }}}^{\mathrm{T}}}\right)\right\}\right\}$$

(40)

We now divide $\tilde{\mathbf{\omega }}$ into

$$\tilde{\mathbf{\omega }}={[{\mathbf{\omega }}_1^{\mathrm{T}},{\mathbf{\omega }}_2^{\mathrm{T}},\dots ,{\mathbf{\omega }}_Q^{\mathrm{T}}]}^{\mathrm{T}}$$

(41)

with

$${\mathbf{\omega }}_q={[{\mathbf{p}}_q^{\mathrm{T}},\mathrm{R}{\mathrm{e}}^{\mathrm{T}}\left\{\mathit{\beta }{\prime }_q\right\},{\mathrm{Im}}^{\mathrm{T}}\left\{\mathit{\beta }{\prime }_q\right\}]}^{\mathrm{T}}$$

(42)

From equation (40), as the sources are assumed to be uncorrelated, the CRB for the parameters of the qth source is

$$\mathbf{CRB}\left({\mathbf{\omega }}_q\right)=\frac{1}{2}{\mathrm{Re}}^{-1}\left\{\sum _{j=1}^J\sum _{k=1}^K\left\{{\left(\frac{\partial \left({\mathbf{A}}_q(j,{\mathbf{p}}_q)\mathit{\beta }{\prime }_q{\overline{s}}_q(j,k)\right)}{\partial {\mathbf{\omega }}_q^{\mathrm{T}}}\right)}^{\mathrm{H}}{\mathbf{R}}_n^{-1}\left(\frac{\partial \left({\mathbf{A}}_q(j,{\mathbf{p}}_q)\mathit{\beta }{\prime }_q{\overline{s}}_q(j,k)\right)}{\partial {\mathbf{\omega }}_q^{\mathrm{T}}}\right)\right\}\right\}$$

(43)

For notational convenience, let us define

$$\begin{array}{l}{\Lambda }_q\left(p_q\right)=\mathrm{blkdiag}\left\{A_q\left(1,p_q\right),A_q\left(2,p_q\right),\dots ,A_q\left(J,p_q\right)\right\}\\ {\overline{s}}_q\left(k\right)={\left[{\overline{s}}_q\left(1,k\right),{\overline{s}}_q\left(2,k\right),\dots ,{\overline{s}}_q\left(J,k\right)\right]}^{\mathrm{T}}\\ {\overline{s}}_q={\left[{\overline{s}}_q^{\mathrm{T}}\left(1\right),{\overline{s}}_q^{\mathrm{T}}\left(2\right),\dots ,{\overline{s}}_q^{\mathrm{T}}\left(K\right)\right]}^{\mathrm{T}}\end{array}$$

(44)

Hence, equation (43) is expressible as

$$\begin{array}{c}\mathbf{CRB}\left({\mathbf{\omega }}_q\right)\hfill \\ =\frac{1}{2}{\mathrm{Re}}^{-1}\left\{\sum _{k=1}^K\left\{{\left(\frac{\partial \left({\mathbf{\Lambda }}_q\left({\mathbf{p}}_q\right)({\mathbf{I}}_J\otimes \mathit{\beta }{\prime }_q){\overline{\mathbf{s}}}_q\left(k\right)\right)}{\partial {\mathbf{\omega }}_q^{\mathrm{T}}}\right)}^{\mathrm{H}}({\mathbf{I}}_J\otimes {\mathbf{R}}_n^{-1})\left(\frac{\partial \left({\mathbf{\Lambda }}_q\left({\mathbf{p}}_q\right)({\mathbf{I}}_J\otimes \mathit{\beta }{\prime }_q){\overline{\mathbf{s}}}_q\left(k\right)\right)}{\partial {\mathbf{\omega }}_q^{\mathrm{T}}}\right)\right\}\right\}\hfill \\ =\frac{1}{2}{\mathrm{Re}}^{-1}\left\{{\left(\frac{\partial \left(({\mathbf{I}}_K\otimes \left({\mathbf{\Lambda }}_q\left({\mathbf{p}}_q\right)({\mathbf{I}}_J\otimes \mathit{\beta }{\prime }_q)\right)){\overline{\mathbf{s}}}_q\right)}{\partial {\mathbf{\omega }}_q^{\mathrm{T}}}\right)}^{\mathrm{H}}{\tilde{\mathbf{R}}}_n^{-1}\left(\frac{\partial \left(({\mathbf{I}}_K\otimes \left({\mathbf{\Lambda }}_q\left({\mathbf{p}}_q\right)({\mathbf{I}}_J\otimes \mathit{\beta }{\prime }_q)\right)){\overline{\mathbf{s}}}_q\right)}{\partial {\mathbf{\omega }}_q^{\mathrm{T}}}\right)\right\}\hfill \end{array}$$

(45)

where ${\tilde{\mathbf{R}}}_n={\mathbf{I}}_{\mathrm{JK}}\otimes {\mathbf{R}}_n$ . Denoting $\mathbf{\Lambda }({\mathbf{\omega }}_q)=({\mathbf{I}}_K\otimes ({\mathbf{\Lambda }}_q({\mathbf{p}}_q)({\mathbf{I}}_J\otimes \mathit{\beta }{\prime }_q))){\overline{\mathbf{s}}}_q$ , $\mathbf{CRB}({\mathbf{\omega }}_q)$ takes the form

$$\begin{array}{c}\mathbf{CRB}\left({\mathbf{\omega }}_q\right)=\frac{1}{2}{\mathrm{Re}}^{-1}\left\{{\left(\frac{\partial \mathbf{\Lambda }\left({\mathbf{\omega }}_q\right)}{\partial {\mathbf{\omega }}_q^{\mathrm{T}}}\right)}^{\mathrm{H}}{\tilde{\mathbf{R}}}_{\mathrm{n}}^{-1}\left(\frac{\partial \mathbf{\Lambda }\left({\mathbf{\omega }}_q\right)}{\partial {\mathbf{\omega }}_q^{\mathrm{T}}}\right)\right\}\hfill \\ =\frac{1}{2}{\mathrm{Re}}^{-1}\left\{{\left({\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{\omega }}_q}\right)}^{\mathrm{H}}{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{\omega }}_q}\right\}\hfill \end{array}$$

(46)

in which ${\mathbf{\Lambda }}_{{\mathbf{\omega }}_q}$ can be divided into

$${\Sigma }_{{\mathrm{\omega }}_q}=\left[{\Sigma }_{p_q},{\Sigma }_{\mathrm{Re\{}{\mathit{\beta }\prime }_q\mathrm{\}}},{\Sigma }_{\mathrm{Im}\left\{{\mathit{\beta }\prime }_q\right\}}\right]\mathrm{=}\left[{\Sigma }_{p_q},{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}},\mathrm{j}{\Sigma }_{\mathrm{Re}\left\{{\mathit{\beta }\prime }_q\right\}}\right]$$

(47)

Appendix 2 provides the derivation of the expressions for the sub-matrices in equation (47). Note that $\mathbf{CRB}({\mathbf{\omega }}_q)$ does not exhibit a block-diagonal structure, and therefore complicates the computation of the compact position-related block of the CRB, which is of most interest in our study. To facilitate the subsequent derivation, we use Proposition 1 to reparametrize the unknown parameters so that the new CRB matrix becomes block-diagonal.

We now introduce a new parameter vector

$${\mathit{\omega }\prime }_q={\left[p_q^{\mathrm{T}},{\left(\mathrm{Re}\left\{{\mathit{\beta }\prime }_q\right\}+\mathrm{Re}\left\{U_q\right\}{p}_q\right)}^{\mathrm{T}},{\left(\mathrm{Im}\left\{{\mathit{\beta }\prime }_q\right\}+\mathrm{Im}\left\{U_q\right\}{p}_q\right)}^{\mathrm{T}}\right]}^{\mathrm{T}}$$

(48)

where

$${\mathbf{U}}_q={\left({\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{\mathrm{Re}\{\mathit{\beta }{\prime }_q\}}\right)}^{+}{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{p}}_q}$$

(49)

From equation (48), the relationship between vectors $\mathbf{\omega }{\prime }_q$ and ${\mathbf{\omega }}_q$ is given by

$$\mathbf{\omega }{\prime }_q={\mathbf{H}}_q{\mathbf{\omega }}_q$$

(50)

where

$${\mathbf{H}}_q=\left[\begin{array}{ccc}\mathbf{I}\hfill & \mathbf{O}\hfill & \mathbf{O}\hfill \\ \multicolumn{1}{c}{\mathrm{Re}\left\{{\mathbf{U}}_q\right\}}& \mathbf{I}\hfill & \mathbf{O}\hfill \\ \multicolumn{1}{c}{\mathrm{Im}\left\{{\mathbf{U}}_q\right\}}& \mathbf{O}\hfill & \mathbf{I}\hfill \end{array}\right]$$

(51)

In light of the proposition, the CRB matrix for vector $\mathbf{\omega }{\prime }_q$ is expressed as

$$\begin{array}{c}\mathbf{CRB}\left(\mathbf{\omega }{\prime }_q\right)={\mathbf{H}}_q\mathbf{CRB}\left({\mathbf{\omega }}_q\right){\mathbf{H}}_q^{\mathrm{T}}\hfill \\ =\frac{1}{2}{\mathbf{H}}_q{\mathrm{Re}}^{-1}\left\{{\left({\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{\omega }}_q}\right)}^{\mathrm{H}}{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{\omega }}_q}\right\}{\mathbf{H}}_q^{\mathrm{T}}\hfill \\ =\frac{1}{2}{\mathrm{Re}}^{-1}\left\{{\left({\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{\omega }}_q}{\mathbf{H}}_q^{-1}\right)}^{\mathrm{H}}\left({\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{\omega }}_q}{\mathbf{H}}_q^{-1}\right)\right\}\hfill \end{array}$$

(52)

where the sub-matrix form of ${\mathbf{H}}_q^{-1}$ can be deduced from equation (51)

$${\mathbf{H}}_q^{-1}=\left[\begin{array}{ccc}\mathbf{I}\hfill & \mathbf{O}\hfill & \mathbf{O}\hfill \\ \multicolumn{1}{c}{-\mathrm{Re}\left\{{\mathbf{U}}_q\right\}}& \mathbf{I}\hfill & \mathbf{O}\hfill \\ \multicolumn{1}{c}{-\mathrm{Im}\left\{{\mathbf{U}}_q\right\}}& \mathbf{O}\hfill & \mathbf{I}\hfill \end{array}\right]$$

(53)

Combining equations (47) and (53) yields, after some algebraic manipulations

$${\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{\omega }}_q}{\mathbf{H}}_q^{-1}=[{\mathbf{\Pi }}_{{\tilde{\mathbf{R}}}_{\mathrm{n}}^{-1/2}{\mathbf{\Lambda }}_{\mathrm{Re}\{\mathit{\beta }{\prime }_q\}}}^{\perp }{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{p}}_q},{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{\mathrm{Re}\{\mathit{\beta }{\prime }_q\}},j{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{\mathrm{Re}\{\mathit{\beta }{\prime }_q\}}]$$

(54)

Then substituting equation (54) back into equation (52), we have

$$\mathbf{CRB}\left(\mathbf{\omega }{\prime }_q\right)=\frac{1}{2}{\left[\begin{array}{cc}{\mathbf{G}}_{q1}\hfill & \mathbf{O}\hfill \\ \multicolumn{1}{c}{\mathbf{O}}& {\mathbf{G}}_{q2}\hfill \end{array}\right]}^{-1}$$

(55)

where

$${\mathbf{G}}_{q1}=\mathrm{Re}\left\{{\mathbf{\Lambda }}_{{\mathbf{p}}_q}^{\mathrm{H}}{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Pi }}_{{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{\mathrm{Re}\{\mathit{\beta }{\prime }_q\}}}^{\perp }{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{p}}_q}\right\}$$

(56)

and

$$G_{q2}=\left[\begin{array}{cc}\mathrm{Re}\left\{{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}}^{\mathrm{H}}{\tilde{R}}_n^{-1}{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}}\right\}& -\mathrm{Im}\left\{{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}}^{\mathrm{H}}{\tilde{R}}_n^{-1}{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}}\right\}\\ \mathrm{Im}\left\{{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}}^{\mathrm{H}}{\tilde{R}}_n^{-1}{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}}\right\}& \mathrm{Re}\left\{{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}}^{\mathrm{H}}{\tilde{R}}_n^{-1}{\Sigma }_{\mathrm{Re}\{{\mathit{\beta }\prime }_q\mathrm{\}}}\right\}\end{array}\right]$$

(57)

Now, the CRB for $\mathbf{\omega }{\prime }_q$ is block-diagonal, and the position-related block of CRB is obtained as

$$\mathbf{CRB}\left({\mathbf{p}}_q\right)=\frac{1}{2}{\mathbf{G}}_{q1}^{-1}=\frac{1}{2}\mathrm{R}{\mathrm{e}}^{-1}\left\{{\mathbf{\Lambda }}_{{\mathbf{p}}_q}^{\mathrm{H}}{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Pi }}_{{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{\mathrm{Re}\{\mathit{\beta }{\prime }_q\}}}^{\perp }{\tilde{\mathbf{R}}}_n^{-1/2}{\mathbf{\Lambda }}_{{\mathbf{p}}_q}\right\}$$

(58)

The trace of $\mathbf{CRB}({\mathbf{p}}_q)$ provides the lower bound of the location MSE for any unbiased estimators. Note that equation (58) is also the MSE of the proposed algorithm, as the WLS-based estimator with the optimal weighting matrix is asymptotically equivalent to the ML estimator for Gaussian noise.

MSE without accounting for reflector position biases

Next, the MSE for the locations is examined when using an optimum algorithm that assumes no reflector position biases, but in fact the reflector positions have small biases. We use a WLS-based estimator to analyze the corresponding MSE. This WLS-based estimator is computed following a similar procedure to that provided in a previous section, except that the reflector position biases are not considered here. We notice that if we replace ${\mathbf{A}}_{\mathrm{R}}(j,\mathbf{p})$ (see equation (34)) with

$${\overline{\mathbf{A}}}_{\mathrm{R}}(j,\mathbf{p})=[{\overline{\mathbf{a}}}_0(j,\mathbf{p}),{\overline{\mathbf{a}}}_1(j,\mathbf{p}),\dots ,{\overline{\mathbf{a}}}_{N_{\mathrm{R}}}(j,\mathbf{p})]\in {\mathbb{C}}^{M\times (N_{\mathrm{R}}+1)}$$

(59)

the obtained estimator becomes the one without accounting for reflector position biases. As the purpose of this subsection is to study the sensitivity to the reflector position biases, we assume that there is no spurious path and hence avoid any effects encountered under the condition of over-modeling. Therefore, the objective function for the qth source becomes

$${\overline{V}}_q(\mathbf{p},{\hat{\mathbf{b}}}_q)={\hat{\mathbf{b}}}_q^{\mathrm{H}}{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\mathbf{\Pi }}_{{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\overline{\mathbf{A}}}_q(\mathbf{p})}^{\perp }{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\hat{\mathbf{b}}}_q$$

(60)

where ${\overline{\mathbf{A}}}_q(\mathbf{p})=[{\overline{\mathbf{A}}}_q^{\mathrm{T}}(1,\mathbf{p}),{\overline{\mathbf{A}}}_q^{\mathrm{T}}(2,\mathbf{p}),\dots ,{\overline{\mathbf{A}}}_q^{\mathrm{T}}(J,\mathbf{p}){]}^{\mathrm{T}}$ with ${\overline{\mathbf{A}}}_q(j,\mathbf{p})=[{\overline{\mathbf{a}}}_0(j,\mathbf{p}),{\overline{\mathbf{a}}}_1(j,\mathbf{p}),\dots ,{\overline{\mathbf{a}}}_{L_q-1}(j,\mathbf{p})]\in {\mathbb{C}}^{M\times L_q}$ . Now, we represent ${\hat{\mathbf{b}}}_q(j)$ as a sum of three vectors

$${\hat{\mathbf{b}}}_q\left(j\right)={\overline{\mathbf{b}}}_q\left(j\right)+{\mathbf{e}}_q\left(j\right)+{\mathbf{\epsilon }}_q\left(j\right)$$

(61)

in which

$$\begin{array}{c}{\overline{\mathbf{b}}}_q\left(j\right)={\overline{\mathbf{A}}}_q(j,{\mathbf{p}}_q){\mathit{\beta }″}_q\hfill \\ {\mathit{\beta }″}_q={[{\beta }_{q0},{\beta \prime }_{q1},\dots ,\beta {\prime }_{q(L_q-1)}]}^{\mathrm{T}}\hfill \end{array}$$

(62)

Here, ${\mathbf{e}}_q(j)$ is a vector associated with the reflector position biases, and ${\mathbf{\epsilon }}_q(j)$ is a stochastic vector associated with measurement noise. Denote ${\overline{\mathbf{b}}}_q=[{\overline{\mathbf{b}}}_q^{\mathrm{T}}(1),{\overline{\mathbf{b}}}_q^{\mathrm{T}}(2),\dots ,{\overline{\mathbf{b}}}_q^{\mathrm{T}}(J){]}^{\mathrm{T}}$ , and therefore we have ${\overline{\mathbf{b}}}_q={\overline{\mathbf{A}}}_q({\mathbf{p}}_q){\mathit{\beta }″}_q$ . Using this result, it can be proved that

$${\overline{V}}_q({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)={\overline{\mathbf{b}}}_q^{\mathrm{H}}{\mathbf{W}}_{{\mathrm{s}}_{q}n}^{1/2}{\mathbf{\Pi }}_{{\mathbf{W}}_{{\mathrm{s}}_{q}n}^{1/2}{\overline{\mathbf{A}}}_q({\mathbf{p}}_q)}^{\perp }{\mathbf{W}}_{{\mathrm{s}}_{q}n}^{1/2}{\overline{\mathbf{b}}}_q=0$$

(63)

Considering the nonlinear relationship between the objective function and p _q , we apply a small error analysis to derive the position error and then take the first-order Taylor expansion of the gradient ${\mathbf{v}}_{\mathbf{p}}({\hat{\mathbf{p}}}_q,{\hat{\mathbf{b}}}_q)=\frac{\partial {\overline{V}}_q(\mathbf{p},\mathbf{b})}{\partial \mathbf{p}}|(\mathbf{p},\mathbf{b})=({\hat{\mathbf{p}}}_q,{\hat{\mathbf{b}}}_q)$ around the point $({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)$

$${\mathbf{v}}_{\mathbf{p}}({\hat{\mathbf{p}}}_q,{\hat{\mathbf{b}}}_q)={\mathbf{v}}_{\mathbf{p}}({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)+{\mathbf{V}}_{\mathbf{pp}}({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)({\hat{\mathbf{p}}}_q-{\mathbf{p}}_q)+2\mathrm{Re}\left\{{\mathbf{V}}_{\mathbf{pb}}({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)({\hat{\mathbf{b}}}_q-{\overline{\mathbf{b}}}_q)\right\}$$

(64)

where ${\mathbf{V}}_{\mathbf{pp}}(\mathbf{p},\mathbf{b})={\partial }^2{\overline{V}}_q(\mathbf{p},\mathbf{b})/\partial \mathbf{p}\partial {\mathbf{p}}^{\mathrm{T}}$ and ${\mathbf{V}}_{\mathbf{pb}}(\mathbf{p},\mathbf{b})={\partial }^2{\overline{V}}_q(\mathbf{p},\mathbf{b})/\partial \mathbf{p}\partial {\mathbf{b}}^{\mathrm{T}}$ . Since the estimator ${\hat{\mathbf{p}}}_q$ is the position minimizing the objective function ${\overline{V}}_q(\mathbf{p},{\hat{\mathbf{b}}}_q)$ , we have ${\mathbf{v}}_{\mathbf{p}}({\hat{\mathbf{p}}}_q,{\hat{\mathbf{b}}}_q)=0$ . Applying this result and equations (63) to (64) leads to

$$V_{pp}\left(p_q,{\overline{b}}_q\right)\left({\widehat{p}}_q-p_q\right)+2\mathrm{Re}\left\{V_{pb}\left(p_q,{\overline{b}}_q\right)\left({\widehat{b}}_q-{\overline{b}}_q\right)\right\}=0$$

(65)

Then, substituting equation (61) into equation (65), the localization error of the qth source, $\Delta {\mathbf{p}}_q={\hat{\mathbf{p}}}_q-{\mathbf{p}}_q$ , can be computed by

$$\Delta {\mathbf{p}}_q=-2{\mathbf{V}}_{\mathbf{pp}}^{-1}({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)\mathrm{Re}\left\{{\mathbf{V}}_{\mathbf{pb}}({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)({\mathbf{e}}_q+{\mathbf{\epsilon }}_q)\right\}$$

(66)

where ${\mathbf{e}}_q=[{\mathbf{e}}_q^{\mathrm{T}}(1),{\mathbf{e}}_q^{\mathrm{T}}(2),\dots ,{\mathbf{e}}_q^{\mathrm{T}}(J){]}^{\mathrm{T}}$ and ${\mathbf{\epsilon }}_q=[{\mathbf{\epsilon }}_q^{\mathrm{T}}(1),{\mathbf{\epsilon }}_q^{\mathrm{T}}(2),\dots ,{\mathbf{\epsilon }}_q^{\mathrm{T}}(J){]}^{\mathrm{T}}$ . Following the similar derivation to that in Amar and Weiss, ¹⁶ we have

$${\mathbf{V}}_{\mathbf{pp}}({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)=2\mathrm{Re}\left\{{\overline{\mathbf{B}}}_{{\mathbf{p}}_q}^H{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\mathbf{\Pi }}_{{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\overline{\mathbf{A}}}_q}^{\perp }{\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\overline{\mathbf{B}}}_{{\mathbf{p}}_q}\right\}$$

(67)

$${\mathbf{V}}_{\mathbf{pb}}({\mathbf{p}}_q,{\overline{\mathbf{b}}}_q)=-{\overline{\mathbf{B}}}_{{\mathbf{p}}_q}^H{\mathbf{W}}_{{\mathrm{s}}_{q}n}^{1/2}{\mathbf{\Pi }}_{{\mathbf{W}}_{{\mathrm{s}}_{q}n}^{1/2}{\overline{\mathbf{A}}}_q}^{\perp }{\mathbf{W}}_{{\mathrm{s}}_{q}n}^{1/2}$$

(68)

where ${\overline{\mathbf{B}}}_{{\mathbf{p}}_q}=\frac{\partial {\overline{\mathbf{A}}}_q(\mathbf{p}){\mathit{\beta }″}_q}{\partial {\mathbf{p}}^{\mathrm{T}}}|\mathbf{p}={\mathbf{p}}_q$ and ${\overline{\mathbf{A}}}_q$ is short for ${\overline{\mathbf{A}}}_q({\mathbf{p}}_q)$ . Inserting equations (67) and (68) into equation (66), we derive the covariance matrix of $\Delta {\mathbf{p}}_q$ , ${\mathbf{C}}_{{\mathbf{p}}_q}=\mathbb{E}[\Delta {\mathbf{p}}_q\Delta {\mathbf{p}}_q^{\mathrm{T}}]$ ; see Appendix 3. The expression for ${\mathbf{C}}_{{\mathbf{p}}_q}$ is given by

$${\mathbf{C}}_{{\mathbf{p}}_q}={\mathbf{C}}_{{\mathbf{p}}_q}^{\mathrm{n}}+{\mathbf{C}}_{{\mathbf{p}}_q}^{\mathrm{e}}$$

(69)

where

$${\mathbf{C}}_{{\mathbf{p}}_q}^n=\frac{1}{2K}\mathrm{R}{\mathrm{e}}^{-1}\left\{{\mathbf{F}}_q^H{\mathbf{F}}_q\right\}$$

(70)

is related to the measurement noise from a finite number of samples ( $K\ne \infty$ ), and

$${\mathbf{C}}_{{\mathbf{p}}_q}^{\mathrm{e}}=2\mathrm{Re}\{{({\mathbf{F}}_q^{\mathrm{H}}{\mathbf{F}}_q+{\mathbf{F}}_q^{\mathrm{T}}{\mathbf{F}}_q^{*})}^{-1}{\mathbf{F}}_q^{\mathrm{H}}{\overline{\mathbf{W}}}_{{\mathrm{s}}_{q}n}^{1/2}\mathbf{R}{\prime }_e{\overline{\mathbf{W}}}_{{\mathrm{s}}_{q}n}^{\mathrm{T}/2}{\mathbf{F}}_q^{*}{({\mathbf{F}}_q^{\mathrm{H}}{\mathbf{F}}_q+{\mathbf{F}}_q^{\mathrm{T}}{\mathbf{F}}_q^{*})}^{-1}+{({\mathbf{F}}_q^{\mathrm{H}}{\mathbf{F}}_q+{\mathbf{F}}_q^{\mathrm{T}}{\mathbf{F}}_q^{*})}^{-1}{\mathbf{F}}_q^{\mathrm{H}}{\overline{\mathbf{W}}}_{{\mathrm{s}}_{q}n}^{1/2}{\mathbf{R}}_e{\overline{\mathbf{W}}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\mathbf{F}}_q{({\mathbf{F}}_q^{\mathrm{H}}{\mathbf{F}}_q+{\mathbf{F}}_q^{\mathrm{T}}{\mathbf{F}}_q^{*})}^{-1}\}$$

(71)

represents the additional error resulted from the reflector position biases. Here, $\mathbf{R}{\prime }_e={\mathbf{e}}_q{\mathbf{e}}_q^{\mathrm{T}}$ , ${\mathbf{R}}_e={\mathbf{e}}_q{\mathbf{e}}_q^{\mathrm{H}}$ , and ${\mathbf{F}}_q={\mathbf{\Pi }}_{{\overline{\mathbf{W}}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\overline{\mathbf{A}}}_q}^{\perp }{\overline{\mathbf{W}}}_{{\mathrm{s}}_q\mathrm{n}}^{1/2}{\overline{\mathbf{B}}}_{{\mathbf{p}}_q}$ , with ${\overline{\mathbf{W}}}_{{\mathrm{s}}_q\mathrm{n}}$ being the true value of ${\mathbf{W}}_{{\mathrm{s}}_q\mathrm{n}}$ .

Remark 2

Note that ${\mathbf{C}}_{{\mathbf{p}}_q}^{\mathrm{n}}$ is inversely proportional to the SNR embedded in ${\overline{\mathbf{W}}}_{{\mathrm{s}}_q\mathrm{n}}$ , whereas ${\mathbf{C}}_{{\mathbf{p}}_q}^{\mathrm{e}}$ is proportional to $\mathbf{R}{\prime }_e$ and ${\mathbf{R}}_e$ . As ${\mathbf{e}}_q$ is assumed to be efficiently small in our study, ${\mathbf{C}}_{{\mathbf{p}}_q}^{\mathrm{n}}$ dominates the location MSE of this estimator when the SNR is rather small. However, with an increase in SNR, the effect of the reflector position biases becomes critical. This result is illustrated using the simulations described in the next section.

Numerical results

Monte Carlo simulations were used to examine the localization performance of the newly proposed DPD algorithm and to validate the CRB and MSE analysis. Here, we consider that the receiver is equipped with a uniform linear array consisting of M = 10 sensors located at $\mathbf{v}=(0,0)(\mathrm{km})$ . The adjacent sensors are spaced at 0.5λ intervals, where λ denotes the wavelength. We generate the simulated signal waveforms with a bandwidth of 17 kHz and identical power ${\sigma }_{\mathrm{s}}^2$ , and the noise is modeled as a circularly complex Gaussian process with covariance matrix R _n . In simulations, free-space propagation is assumed, where the path-loss attenuation magnitude is inversely proportional to the squared distance. For different scenarios, each simulation performs 500 Monte Carlo trials, and the location is estimated based on 240 samples (K = 15 and J = 16). We evaluate the performance using the root mean square error (RMSE) of the locations, that is

$$\mathrm{RMS}{\mathrm{E}}_q=\sqrt{\frac{1}{500}\sum _{i=1}^{500}{\Vert {\mathbf{p}}_q-{\hat{\mathbf{p}}}_q^i\Vert }^2}$$

(72)

where q denotes the qth source, and ${\hat{\mathbf{p}}}_q^i$ indicates the estimation of ${\mathbf{p}}_q$ in the ith Monte Carlo trial.

Sensitivity to spurious paths

The first experiment examines the impact of spurious paths on the estimation accuracy of the proposed DPD algorithm. The noise is spatially uncorrelated with covariance matrix ${\mathbf{R}}_n={\sigma }_{\mathrm{n}}^2{\mathbf{I}}_M$ , where ${\sigma }_{\mathrm{n}}^2$ denotes the power. We consider a scenario of one transmitter placed at p₁ = [1,6]^T (km) and one true reflector placed at the nominal position u₁ = [−2,3]^T (km) with an unknown bias of Δu₁₁ = [0.05,0.1]^T (km) (sufficiently small). For the spurious paths within the area of interest, three cases are simulated: no spurious path, two spurious indirect paths, and four spurious indirect paths. Table 1 lists the nominal positions of the spurious reflectors for each case, and Figure 2 summarizes the corresponding geometries, where the solid and dash-dot lines indicate the true and spurious paths, respectively. In all three cases, we evaluate the RMSE of the proposed algorithm and the associated CRB with respect to SNR. As shown in Figure 3, the performance of our algorithm is not sensitive to the spurious paths because there is no significant difference between the RMSE curves for different numbers of spurious paths. In particular, the curve corresponding to the two spurious paths overlaps with the curve corresponding to no spurious path.

Table 1.

The nominal positions of spurious reflectors in three cases.

Case	X (km)	Y (km)
Case 1	None	None
Case 2	–2	1
	–2	5
Case 3	–2	1
	–2	5
	3	2
	3	4

Figure 2.

Problem geometry of a single source (pentagram), a true reflector (triangle), and different spurious reflectors (squares) for three cases: (a) case 1, (b) case 2, and (c) case 3.

Figure 3.

Square root of the CRB and the estimated RMSEs versus SNR in spatially uncorrelated noise for different numbers of spurious paths (L_S denotes the number of spurious paths in the legend).

Location accuracy in the presence of small reflector position biases

In this set of experiments, we consider two correlation structures for noise:

• The noise is spatially uncorrelated with covariance matrix ${\mathbf{R}}_n={\sigma }_{\mathrm{n}}^2{\mathbf{I}}_M$ .

• The noise is spatially correlated with exponentially decaying correlations among the elements, ¹⁴ and the element of ${\mathbf{R}}_n$ is given by

$${\left[{\mathbf{R}}_n\right]}_{\mathrm{lk}}={\sigma }_{\mathrm{n}}^{2}0.9^{|l-k|}\exp \left(\mathrm{j}\frac{\mathrm{\pi }}{2}(l-k)\right)$$

(73)

To make our results as comprehensive as possible, we employ the following estimators for a comparison with the proposed DPD algorithm:

1. The two-step processing estimator—The DOA and TOA estimation using the expectation–maximization algorithm ²¹ with known waveforms, and the least squares (LS) location with the DOA and TOA estimates used as data. The corresponding plot of the results is labeled by “Two-step with known sig.”

2. The DPD algorithm with unknown signal waveforms, ¹³ for which the plot of the results is labeled by “DPD with unknown sig.”

3. The DPD algorithm for a single transmitter with known signal waveform, ³ for which the plot is labeled by “DPD with known sig.”

For the two-step processing estimator, the ideal data association is assumed as this task is not easy in practice. ²² Without the a priori true paths, the LS location in the second step performs exceedingly poorly, whereas the DPD algorithm avoids this problem. In the remaining part, we shall assume that there is no spurious path. The curves of the CRB and the theoretical RMSE without accounting for the reflector position biases are also displayed. The latter provides a benchmark for the localization accuracy when the reflector positions are assumed accurate but in fact have small biases.

First, the case of a single transmitter is considered and the geometry for the layout is the same as that of case 1 in the preceding experiment. In the spatially uncorrelated and correlated noise environments, we compare the RMSE of our algorithm with those of the three prescribed positioning algorithms as SNR increases. As illustrated in Figure 4, our algorithm holds the lowest RMSE, followed by the DPD with known signal, ³ whereas the DPD with unknown signal ¹³ exhibits the worst performance. Notice that the performance advantage of the proposed algorithm becomes increasingly prominent as SNR increases. As reflector position biases are present, other algorithms exhibit irreducible location errors at high SNRs. When SNR equals 30 dB, the improvement gain of the proposed estimator is more than 20 dB compared with the DPD with known signal in uncorrelated noise, whereas this improvement gain reaches almost 30 dB in correlated noise. Moreover, the theoretical RMSE without accounting for reflector position biases and the simulated RMSE of the DPD with known signal match very well in regard to spatially uncorrelated noise, thereby verifying the theoretical result (see equation (69)). For spatially correlated noise, the DPD with known signal has deteriorated severely, but the performance of the proposed algorithm approaches the CRB with insignificant differences. This demonstrates the robustness of the proposed algorithm to unknown colored noise.

Figure 4.

Square root of the CRB, the estimated RMSEs, and the theoretical RMSE without accounting for reflector position biases versus SNR (Q = 1): (a) in spatially uncorrelated noise and (b) in spatially correlated noise.

Next, we examine the performance of the proposed algorithm for multiple transmitters. Without loss of generality, two transmitters placed at p₁ = [1,6]^T (km) and p₂ = [3,2]^T (km) are considered. Figure 5 displays the layout of this experiment. Each signal is received via two paths: the direct path and the indirect path that is created by the same reflector at the nominal position u₁ = [−2,3]^T (km). The unknown reflector position biases are Δu₁₁ = [0.05,0.1]^T (km) and Δu₂₁ = [0.06,0.105]^T (km) for the indirect paths of the two transmitters. We plot the RMSE curves for spatially uncorrelated and correlated noises in Figures 6 and 7. As the DPD with known signal was proposed to locate a single source, we do not present its performance here.

Figure 5.

Problem geometry of two sources and a single true reflector.

Figure 6.

Square root of the CRB, the estimated RMSEs, and the theoretical RMSE without accounting for reflector position biases versus SNR in spatially uncorrelated noise (Q = 2): (a) source 1 and (b) source 2.

Figure 7.

Square root of the CRB, the estimated RMSEs, and the theoretical RMSE without accounting for reflector position biases versus SNR in spatially correlated noise (Q = 2): (a) source 1 and (b) source 2.

From Figures 4, 6, and 7, the theoretical RMSE without accounting for reflector position biases overlaps with the simulated RMSE of the proposed algorithm at low SNRs. However, the gap between them widens as SNR increases, where the effect of the reflector position biases becomes more critical. This is consistent with the analysis in the previous section. Compared with the CRB (or the simulated RMSE of the proposed algorithm), the increases in the theoretical RMSE without accounting for reflector position biases are more than 20 dB for each source location at SNR = 30 dB. Furthermore, it can be seen that our algorithm is considerably superior to the DPD with unknown signals throughout the whole range of SNR. In explanation, the extra information contained in the waveforms is helpful in enhancing the position estimation.

Location accuracy in the presence of larger reflector position biases

The preceding results demonstrate the superiority of the proposed DPD in the presence of small reflector position biases. We now come to study its performance when this small-bias assumption is not satisfied. In this set of experiments, the problem geometry is the same as that presented in Figure 5, and the reflector position biases are expressed as

$$\begin{array}{c}\Delta {\mathbf{u}}_{11}=(1,2)\Delta \sigma \left(\mathrm{km}\right)\hfill \\ \Delta {\mathbf{u}}_{21}=(1.2,2.1)\Delta \sigma \left(\mathrm{km}\right)\hfill \end{array}$$

(74)

where $\Delta \sigma$ is regarded as the bias level. The noise is spatially uncorrelated having a scaled identity matrix covariance, and the SNR is fixed at 5 dB. With $\Delta \sigma$ varying from 0.02 to 0.2 km, we assess the location accuracy of the proposed estimator, the two-step estimator with known signals and known reflector position biases, and the two-step estimator with known signals but unknown reflector position biases. The RMSE curves are plotted in Figure 8. Unsurprisingly, the bias level has no significant effect on the performance of the two-step estimator with known reflector position biases. On the contrary, the estimation accuracy of the two-step estimator with unknown reflector position biases diminishes as $\Delta \sigma$ increases. Up to $\Delta \sigma$  = 0.14 km, our algorithm performs well and its RMSEs remain approximately unchanged. However, the location performance of the proposed algorithm deteriorates for sufficiently large $\Delta \sigma$ , where a model mismatch is encountered, that is, the approximation of the signal model fails. This result confirms the need for the small-bias assumption. Moreover, as our algorithm introduces more unknown parameters to account for small reflector position biases, its performance deteriorates faster than that of the two-step method when $\Delta \sigma$  > 0.16 km.

Figure 8.

Square root of the CRB and the estimated RMSEs versus the bias level in spatially uncorrelated noise (Q = 2): (a) source 1 and (b) source 2.

Conclusion

This study proposed a reflector-aided direct geolocation algorithm of multiple sources received by an observer array. Practical scenarios were taken into consideration, where the nominal positions of reflectors are available but the small reflector position biases exist, and the correlation structure of noise is unknown. The algorithm is computationally efficient in that (1) it avoids the estimation of intermediate parameters, the dimension of which increases with multipath number, and (2) it decouples the locations of the multiple sources into several lower dimensional optimization problems with the information of uncorrelated waveforms. In addition to this algorithm, we derived a compact CRB formula for source positions and the MSE of the WLS-based estimator without accounting for the reflector position biases. The MSE result indicates that inaccurate knowledge of the reflector positions leads to a decrease in source location accuracy. With extensive numerical examples, we showed that the proposed algorithm is not sensitive to over-modeling as it removes spurious paths intuitively. The performance superiority of our algorithm was demonstrated in spatially uncorrelated and correlated noise environments when the reflector position biases are small. However, we found that large biases in the reflector positions generate large errors in the proposed estimator. To deal with this problem, one possible direction for future work is to develop a Bayesian approach given the prior knowledge of the statistical properties (e.g. probability distribution) of the biases.