Localization and Tracking of Discrete Mobile Scatterers in Vehicular Environments Using Delay Estimates ^†

Abstract

This paper describes an approach to detect, localize, and track moving, non-cooperative objects by exploiting multipath propagation. In a network of spatially distributed transmitting and receiving nodes, moving objects appear as discrete mobile scatterers. Therefore, the localization of mobile scatterers is formulated as a nonlinear optimization problem. An iterative nonlinear least squares algorithm following Levenberg and Marquardt is used for solving the optimization problem initially, and an extended Kalman filter is used for estimating the scatterer location recursively over time. The corresponding performance bounds are derived for both the snapshot based position estimation and the nonlinear sequential Bayesian estimation with the classic and the posterior Cramér–Rao lower bound. Thereby, a comparison of simulation results to the posterior Cramér–Rao lower bound confirms the applicability of the extended Kalman filter. The proposed approach is applied to estimate the position of a walking pedestrian sequentially based on wideband measurement data in an outdoor scenario. The evaluation shows that the pedestrian can be localized throughout the scenario with an accuracy of $0.8$ m at 90% confidence.

1. Introduction

With the trends of increasing urbanization and increasing automation in road transportation, the demand for improvements in vehicular safety technologies is steadily growing. In particular the mixed-traffic environment shared by many different users including cars, motorcycles, cyclists, and pedestrians is challenging for any means of automated transport. To safely route through vehicular environments, therefore, timely and reliable information about other road users is required. In this regard, the exchange of user specific information, like position and velocity, enhances safety on roads by supporting mutual awareness [1,2,3]. This cooperative approach requires road users to be equipped with actively probing devices to determine their user specific information and to exchange the data. However, many road users do not carry any devices, i.e., they are non-cooperative. Thus, for a reliable and comprehensive situational awareness, further methods and sensor technologies are needed, accounting for non-cooperative road users. With regard to current automated and autonomous vehicles, the perception of their surrounding environment mainly relies on locally mounted sensors, including radar and lidar sensors, as well as camera based systems [4]. Due to physical properties, however, locally mounted perception sensors exhibit a series of critical limitations, like the limited performance in adverse lighting conditions [3]. Thus, to improve the reliability and to extend the awareness range of the local perception sensors, infrastructure based systems have been suggested, for instance based on dedicated radar sensors [5] and cameras [6]. Apart from remaining technical challenges, the main drawback of these systems is their deployment. For sufficiently sensing only a limited area like an urban intersection, several sensors would have to be mounted, which results in high deployment and maintenance costs.

Therefore, a localization system is introduced in [7], which reuses signals from vehicular communications infrastructure for passive radar application. It is assumed that road users and other objects affect the radio spectrum by inducing delayed and Doppler shifted MPC. This means that the characteristics of the MPC correspond to location and dynamics of the road users. Thus, sensing the wireless propagation characteristics between the links of existing communication networks allows for detecting and localizing road users. Similarly, the authors in [8] propose the usage of existing vehicular radio links for detecting and localizing road users. In addition to static links from communication infrastructure, they propose to include also mobile devices as possible network nodes. The localization accuracy of both systems strongly depends on precise location information of the individual network nodes. Thus, an incorporation of mobile nodes requires to account for location uncertainties, which can degrade the localization performance.

The proposed sensing systems of [7] and [8] follow the idea of PCL, i.e., the usage of arbitrary signals as illuminators of opportunity [9]. With its origin in aeronautics, PCL was mainly used for observing targets far from the sensing network [9,10,11,12]. The localization requirements, for instance, from aeronautic and maritime applications allow high integration times and low bandwidths. Recently, the interest in passive localization approaches for short-range and passive indoor localization is increasing [13,14,15]. The main challenge of short-range PCL are the dynamics of targets moving in the proximity of the sensing network. Thereby, signals from Wi-Fi access points with comparatively large bandwidths are beneficial for localization, since lower integration times are required [7,13]. In this regard, the authors of [16] propose to use even ultra-wideband signals for PCL. Furthermore, they introduce different target tracking algorithms and analyze the algorithms in simulations. Another localization approach was introduced in [17], applying an iterative, nonlinear least-squares approach for position estimation. The applicability of the approach is demonstrated using wideband measurement data.

On the basis of [17], this work extends the localization approach by sequential target tracking, namely an EKF. For the evaluation of the tracking algorithm, the PCRLB is derived. Finally, the proposed localization approach is applied to wideband measurement data. The main contributions of this paper are:

• The signal processing for a localization approach to localize moving, non-cooperative objects recursively, using delay estimates from a network of spatially distributed transmitting and receiving nodes.
• The derivation of performance bounds including the CRLB on position estimation and the PCRLB on nonlinear sequential Bayesian estimation.
• The validation of the applicability of an EKF for the introduced localization problem via Monte Carlo simulations and a comparison to the PCRLB.
• The application of the proposed localization approach to wideband measurement data of an outdoor experiment for localizing a walking pedestrian.

The remainder of this paper is structured as follows. Section 2 introduces the network structure together with the measurement model used for position estimation. In Section 3, the signal processing is described in order to localize mobile scatterers, which comprises a procedure to extract the measurement vector, a snapshot based position estimator, as well as a nonlinear sequential Bayesian estimator. Corresponding performance bounds for positioning and for nonlinear sequential Bayesian estimation are derived in Section 4. Based on an exemplary measurement setup, the localization approach is evaluated in Section 5, first theoretically by using the derived performance bounds and second by applying the approach to channel measurements. Finally, Section 6 concludes the paper by summarizing the findings.

2. Network and Measurement Model

This section introduces the network structure and the measurement model used for the position estimation. Refer to a widely distributed network of K transmitting and L receiving nodes, where both transmitting and receiving nodes are assumed to be static at known locations ${\mathbf{r}}_k^{\mathrm{Tx}}$, $k\in \{1,\cdots ,K\}$, and ${\mathbf{r}}_l^{\mathrm{Rx}}$, $l\in \{1,\cdots ,L\}$, respectively. Receiving nodes can be collocated with transmitting nodes or individually placed. The network link configuration determines the index set $\mathcal{P}$, where each link $(k,l)\in \mathcal{P}$ is composed of the k-th transmitting node and the l-th receiving node. Accordingly, fully meshed networks result in a maximum number of $\left|\mathcal{P}\right|=KL$ links.

Each transmitting node emits known signals $s_k\left(t\right)$ with period $T_{\mathrm{p}}$, which allows for measuring the CIR at the receivers [18]. Thereby, the received signal of each link in the network is modeled as a superposition of a finite number of scaled and delayed replica of the transmit signal. These comprise the LoS and a discrete number of MPCs due to reflection and scattering. Following [19], the MPCs are further differentiated according to the dynamics of the scattering objects, i.e., static and mobile. Eventually, the CIR for a pair of transmitting and receiving nodes ${\mathrm{Tx}}_k$ and ${\mathrm{Rx}}_l$ can be expressed as

$$h_{kl}(t,\tau )=h_{kl}^{\mathrm{LoS}}(t,\tau )+h_{kl}^{\mathrm{S}}(t,\tau )+h_{kl}^{\mathrm{M}}(t,\tau ),$$

(1)

where $h_{kl}^{\mathrm{LoS}}(t,\tau )$, $h_{kl}^{\mathrm{S}}(t,\tau )$, and $h_{kl}^{\mathrm{M}}(t,\tau )$ denote the contribution of the LoS, the sum of $R_{kl}$ discrete static MPC and of $Q_{kl}$ discrete mobile MPC to the CIR, respectively.

The static network allows for attributing mobile MPC to moving scattering objects in the observed area. Here, the number of mobile MPC is assumed to be identical for each link in the network, i.e., $Q_{kl}=Q$, and the mobile MPC can be uniquely assigned to the moving objects. Consequently, for time and phase synchronized transmitting and receiving nodes, discrete mobile scatterers are contributing to the CIR by

$$h_{kl}^{\mathrm{M}}(t,\tau )=\sum _{q=1}^Q{\alpha }_{klq}^{\mathrm{M}}\left(t\right)e^{-j2\pi f_{\mathrm{c}}{\tau }_{kl}^{\mathrm{M}}\left({\mathbf{r}}_q\left(t\right)\right)}\delta (\tau -{\tau }_{kl}^{\mathrm{M}}\left({\mathbf{r}}_q\left(t\right)\right)),$$

(2)

with ${\alpha }_{klq}^{\mathrm{M}}\left(t\right)$ and ${\tau }_{kl}^{\mathrm{M}}\left({\mathbf{r}}_q\left(t\right)\right)$ representing the complex amplitude and the propagation delay associated with the qth mobile scatterer at location ${\mathbf{r}}_q\left(t\right)={[x_q\left(t\right),y_q\left(t\right)]}^T$, and $\delta \left(·\right)$ representing the Dirac function. The corresponding delay-induced phase shift for center frequency $f_{\mathrm{c}}$ is expressed by the term $e^{-j2\pi f_{\mathrm{c}}{\tau }_{klq}^{\mathrm{M}}\left({\mathbf{r}}_q\left(t\right)\right)}$. For convenience, the notation for time dependence will be omitted in the remainder of this paper and only applied where explicitly needed.

For any pair of transmitting and receiving node, the propagation delay induced by scatterer q is determined by the physical propagation path from transmitter to scatterer and from scatterer to receiver. Thus, given the distances $d_{kq}^{\mathrm{Tx}}$ and $d_{lq}^{\mathrm{Rx}}$ between scatterer and ${\mathrm{Tx}}_k$ and ${\mathrm{Rx}}_l$, respectively, the propagation delay can be expressed as

$$\begin{array}{cc}\hfill {\tau }_{kl}^{\mathrm{M}}\left({\mathbf{r}}_q\right)& =\frac{1}{c}\left(d_{kq}^{\mathrm{Tx}}+d_{lq}^{\mathrm{Rx}}\right)\hfill \\ \hfill & =\frac{1}{c}\left(\parallel {\mathbf{r}}_q-{\mathbf{r}}_k^{\mathrm{Tx}}\parallel +\parallel {\mathbf{r}}_q-{\mathbf{r}}_l^{\mathrm{Rx}}\parallel \right),\hfill \end{array}$$

(3)

where c denotes the speed of light and the operator $\parallel ·\parallel$ denotes the Euclidean norm. For $\left|\mathcal{P}\right|$ network links, linear independent propagation delays induced by the qth scatterer compose a vector ${\mathit{\tau }}_q$, with each vector element ${\tau }_{klq}^{\mathrm{M}}={\tau }_{kl}^{\mathrm{M}}\left({\mathbf{r}}_q\right)$. Thus, for the mapping $\mathit{h}\left({\mathbf{r}}_q\right)={\tau }_q$ with each element defined according to Equation (3), the measurement model is given by

$${\widehat{\mathit{\tau }}}_q={\mathit{\tau }}_q+{\mathbf{w}}_q=\mathit{h}\left({\mathbf{r}}_q\right)+{\mathbf{w}}_q,$$

(4)

with ${\widehat{\mathit{\tau }}}_q$ as corresponding vector of measured delays and ${\mathbf{w}}_q$ as zero-mean white Gaussian noise with covariance matrix ${\mathbf{R}}_q$ defined as diagonal matrix with elements ${\left\{{\sigma }_{klq}^2\right\}}_{(k,l)\in \mathcal{P}}$.

3. Localization and Tracking

As shown above, propagation delays of time-variant MPC inherently contain location information of mobile scatterers being reflected in the CIR. Besides time-variant MPC, the CIR of each link in the network also comprise a static LoS and static MPC as stated in Equation (1). Thus, for localizing and tracking mobile scatterers, time-variant channel components need to be initially extracted. Therefore, the proposed localization approach is composed of three stages—first, a calibration stage to identify and characterize the static LoS and MPC, second, an estimation stage to extract mobile MPC, and third, a tracking stage to estimate the scatterer position recursively.

Both calibration and estimation stage rely on estimates of channel parameters, including complex amplitude and propagation delays of the LoS and MPC. To estimate and track the channel parameters of the CIR, KEST is used [18].

3.1. Calibration Stage

In order to deduce location information of moving objects from CIR measurements, the propagation effects of the static environment need to be known. Therefore, the channel of every link in the network is initially observed over a calibration period $T_{\mathrm{cal}}$. During this period, the environment is assumed to be devoid of any moving object, i.e., $Q=0$. With $T_{\mathrm{g}}$ as the time interval between two adjacent CIR measurements, a total of $\lfloor T_{\mathrm{cal}}/T_{\mathrm{g}}\rfloor$ consecutive CIR are collected for calibration. Given the set of recorded CIR, first, the channel parameters are estimated using KEST. Second, the parameter estimates are clustered with regard to amplitude and delay [20]. For static environments, the resulting clusters correspond to the LoS and to static MPC. Thereby, the vectors ${\overline{\tau }}_{klr}$ and ${\overline{\sigma }}_{klr}$, $r\in \{0,\cdots ,R_{kl}\}$, contain cluster mean and standard deviation for the LoS ($r=0$) and $R_{kl}$ static MPC. Eventually, the vectors ${\overline{\mathit{\tau }}}_{kl}={[{\overline{\tau }}_{kl0},\cdots ,{\overline{\tau }}_{kl{R}_{kl}}]}^T$ and ${\overline{\mathit{\sigma }}}_{kl}={[{\overline{\sigma }}_{kl0},\cdots ,{\overline{\sigma }}_{kl{R}_{kl}}]}^T$ uniquely characterize the static propagation environment between k-th transmitter ${\mathrm{Tx}}_k$ and l-th receiver ${\mathrm{Rx}}_l$. Note that, even though the static environment is only typically changing very slowly, any modification in the propagation condition requires to recalibrate the system. Particularly, objects with strongly reflecting characteristics, will change the conditions significantly, as for example, the placement of a car in the network environment.

3.2. Estimation Stage

Other than during the calibration period, now, additional objects can move within the observed environment. Thus, the CIR and therefore also the estimated channel parameter comprise the propagation effects of moving objects. Similar to before, KEST is applied for estimating channel parameters from incoming CIR. The resulting vector of propagation delay estimates is composed as

$${\widehat{\mathit{\tau }}}_{kl}^{\mathrm{full}}={[\underset{\mathrm{LoS}\mathrm{and}\mathrm{static}\mathrm{MPC}}{\underbrace{{\widehat{\tau }}_{kl0}^{\mathrm{S}},\cdots ,{\widehat{\tau }}_{kl{R}_{kl}}^{\mathrm{S}}}},{\widehat{\tau }}_{kl1}^{\mathrm{M}},\cdots ,{\widehat{\tau }}_{klQ}^{\mathrm{M}}]}^T.$$

(5)

Consequently, extracting mobile MPCs from ${\widehat{\mathit{\tau }}}_{kl}^{\mathrm{full}}$ means to sort out static components. The sorting is based on ${\overline{\mathit{\tau }}}_{kl}$ determined in the preceding calibration stage. Particularly, all elements of ${\widehat{\mathit{\tau }}}_{kl}^{\mathrm{full}}$ lying in an interval of ${\overline{\tau }}_{klp}\pm 3{\overline{\sigma }}_{klp}$ are excluded. Here, the 3-$\sigma$ interval ensures that delay estimates assigned as static are not considered as mobile MPC with a probability higher than 99%, assuming that the previously determined amount of static MPCs remains constant. Eventually, the measurement vector of mobile MPCs for the link ${\mathrm{Tx}}_k$ and ${\mathrm{Rx}}_l$ is

$${\widehat{\mathit{\tau }}}_{kl}={[{\widehat{\tau }}_{kl1}^{\mathrm{M}},\cdots ,{\widehat{\tau }}_{klQ}^{\mathrm{M}}]}^T.$$

(6)

Associated elements of ${\widehat{\mathit{\tau }}}_{kl}$ to scatterer q are rearranged over all links in the network, $\mathcal{P}$, which finally determine the measurement vector ${\widehat{\mathit{\tau }}}_q$.

Based on the measurement model in Equation (4), scatterer q can be localized applying maximum likelihood estimation. In particular, a weighted nonlinear least-squares approach [21,22] is used to minimize the cost function

$$\mathbf{L}\left({\mathbf{r}}_q\right)={({\widehat{\mathit{\tau }}}_q-{\mathit{\tau }}_q)}^T{\mathbf{R}}_q^{-1}({\widehat{\mathit{\tau }}}_q-{\mathit{\tau }}_q),$$

(7)

with respect to the unknown position ${\mathbf{r}}_q$, which is expressed as

$${\widehat{\mathbf{r}}}_q=\underset{{\mathbf{r}}_q}{\arg \min } \mathbf{L}\left({\mathbf{r}}_q\right).$$

(8)

The two-dimensional, nonlinear optimization problem of Equation (8) needs to be solved by an iterative approach, since no closed-form solution is existing. Thus, in order to minimize the cost function in Equation (7), the proposed localization approach applies the Levenberg–Marquardt algorithm [22] due to high robustness and fast convergence characteristics. Particularly, the iterative procedure can be written as

$$\begin{array}{c}\hfill {\mathbf{r}}_q^{(i+1)}={\mathbf{r}}_q^{\left(i\right)}+{\left({\mathbf{J}}^T\left({\mathbf{r}}_q^{\left(i\right)}\right){\mathbf{R}}_q^{-1}\mathbf{J}\left({\mathbf{r}}_q^{\left(i\right)}\right)+{\lambda }^{\left(i\right)}\mathbf{I}\right)}^{-1}\\ \hfill {\mathbf{J}}^T\left({\mathbf{r}}_q^{\left(i\right)}\right){\mathbf{R}}_q^{-1}\left({\widehat{\tau }}_q-{\tau }_q^{\left(i\right)}\right),\end{array}$$

(9)

with $\mathbf{J}\left({\mathbf{r}}_q\right)$, $\mathbf{I}$ and ${\lambda }^{\left(i\right)}$ denoting the Jacobian matrix, the identity matrix, and the dampening parameter for iteration step i, respectively. Individual elements of vector ${\mathit{\tau }}_q^{\left(i\right)}$ are calculated according to Equation (3) as ${\tau }_{kl}^{\mathrm{M}}\left({\mathbf{r}}_q^{\left(i\right)}\right)$. The elements of Jacobian matrix $\mathbf{J}\left({\mathbf{r}}_q\right)\in {\mathbb{R}}^{\left|\mathcal{P}\right|\times 2}$, i.e., the partial derivatives of ${\mathit{\tau }}_q$ with respect to ${\mathbf{r}}_q$, are given as

$$\mathbf{J}\left({\mathbf{r}}_q\right)=\left[\begin{array}{cc}\frac{x_q-x_1^{\mathrm{Tx}}}{d_{1q}^{\mathrm{Tx}}}+\frac{x_q-x_1^{\mathrm{Rx}}}{d_{1q}^{\mathrm{Rx}}}& \frac{y_q-y_1^{\mathrm{Tx}}}{d_{1q}^{\mathrm{Tx}}}+\frac{y_q-y_1^{\mathrm{Rx}}}{d_{1q}^{\mathrm{Rx}}}\\ ⋮& ⋮\\ \frac{x_q-x_K^{\mathrm{Tx}}}{d_{Kq}^{\mathrm{Tx}}}+\frac{x_q-x_1^{\mathrm{Rx}}}{d_{1q}^{\mathrm{Rx}}}& \frac{y_q-y_K^{\mathrm{Tx}}}{d_{Kq}^{\mathrm{Tx}}}+\frac{y_q-y_1^{\mathrm{Rx}}}{d_{1q}^{\mathrm{Rx}}}\\ \frac{x_q-x_1^{\mathrm{Tx}}}{d_{1q}^{\mathrm{Tx}}}+\frac{x_q-x_2^{\mathrm{Rx}}}{d_{2q}^{\mathrm{Rx}}}& \frac{y_q-y_1^{\mathrm{Tx}}}{d_{1q}^{\mathrm{Tx}}}+\frac{y_q-y_2^{\mathrm{Rx}}}{d_{2q}^{\mathrm{Rx}}}\\ ⋮& ⋮\\ \frac{x_q-x_K^{\mathrm{Tx}}}{d_{Kq}^{\mathrm{Tx}}}+\frac{x_q-x_2^{\mathrm{Rx}}}{d_{2q}^{\mathrm{Rx}}}& \frac{y_q-y_K^{\mathrm{Tx}}}{d_{Kq}^{\mathrm{Tx}}}+\frac{y_q-y_2^{\mathrm{Rx}}}{d_{2q}^{\mathrm{Rx}}}\\ ⋮& ⋮\\ \frac{x_q-x_1^{\mathrm{Tx}}}{d_{1q}^{\mathrm{Tx}}}+\frac{x_q-x_L^{\mathrm{Rx}}}{d_{Lq}^{\mathrm{Rx}}}& \frac{y_q-y_1^{\mathrm{Tx}}}{d_{1q}^{\mathrm{Tx}}}+\frac{y_q-y_L^{\mathrm{Rx}}}{d_{Lq}^{\mathrm{Rx}}}\\ ⋮& ⋮\\ \frac{x_q-x_K^{\mathrm{Tx}}}{d_{Kq}^{\mathrm{Tx}}}+\frac{x_q-x_L^{\mathrm{Rx}}}{d_{Lq}^{\mathrm{Rx}}}& \frac{y_q-y_K^{\mathrm{Tx}}}{d_{Kq}^{\mathrm{Tx}}}+\frac{y_q-y_L^{\mathrm{Rx}}}{d_{Lq}^{\mathrm{Rx}}}\end{array}\right].$$

(10)

3.3. Tracking Stage

The previous stage provides an approach to localize a moving scatterer based on delay measurements at one specific time instance t. Since the goal of this work is to localize moving scatterers, it is reasonable to additionally take into account the mobility of the object by filtering the state evolution over time. For any mobile scatterer q, the state is defined by

$${\mathbf{x}}_q\left(t_n\right)={\left[{\mathbf{r}}_q^T\left(t_n\right),{\mathbf{v}}_q^T\left(t_n\right)\right]}^T,$$

(11)

where ${\mathbf{r}}_q\left(t_n\right)$ and ${\mathbf{v}}_q\left(t_n\right)={\left[v_{q,x}\left(t_n\right),v_{q,y}\left(t_n\right)\right]}^T$ denote the position and the velocity of the scatterer at time instant $t_n$. To describe the state evolution from time instant $t_{n-1}$ to time instant $t_n$, a transition model is applied. Accounting for the mobility of scatterers induced by different road users, such as cars, bikes, and pedestrians, a white noise acceleration model is used [23]. The state equation results in

$${\mathbf{x}}_q\left(t_n\right)={\mathbf{A}}_q{\mathbf{x}}_q\left(t_{n-1}\right)+{\mathbf{n}}_q\left(t_n\right),$$

(12)

with transition matrix ${\mathbf{A}}_q$ and zero-mean white Gaussian process noise ${\mathbf{n}}_q\left(t_n\right)$ with covariance matrix ${\mathbf{Q}}_q$. Being $T_g=t_n-t_{n-1}$, the transition matrix is given by

$${\mathbf{A}}_q=\left[\begin{array}{cccc}1& 0& T_g& 0\\ 0& 1& 0& T_g\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(13)

and the covariance matrix by

$${\mathbf{Q}}_q=\left[\begin{array}{cccc}{\sigma }_q^2\frac{T_g^3}{3}& 0& {\sigma }_q^2\frac{T_g^2}{2}& 0\\ 0& {\sigma }_q^2\frac{T_g^3}{3}& 0& {\sigma }_q^2\frac{T_g^2}{2}\\ {\sigma }_q^2\frac{T_g^2}{2}& 0& {\sigma }_q^2{T}_g& 0\\ 0& {\sigma }_q^2\frac{T_g^2}{2}& 0& {\sigma }_q^2{T}_g\end{array}\right],$$

(14)

with ${\sigma }_q^2$ as process noise intensity of physical dimension $\left[{\mathrm{m}}^2/\mathrm{s}3\right]$, which needs to be set according to application requirements [23].

Given state vector ${\mathbf{x}}_q\left(t_n\right)$, the nonlinear mapping of the measurement model in Equation (4) can be expressed as $\mathit{h}\left({\mathbf{r}}_q\left(t_n\right)\right)=\mathit{h}\left({\mathbf{x}}_q\left(t_n\right)\right)$, equivalently. Accordingly, the measurement model is given by

$${\mathbf{z}}_q\left(t_n\right)=\mathit{h}\left({\mathbf{x}}_q\left(t_n\right)\right)+{\mathbf{w}}_q\left(t_n\right).$$

(15)

Due to the nonlinearity of the measurement model, a recursive Bayesian filter is required, which is able to handle general nonlinear problems. A common implementation of such recursive Bayesian filters is the EKF. Given the nonlinear measurement model in Equation (4) and the linear state model in Equation (12), the EKF results in the following set of equations:

$$\begin{array}{cc}\hfill {\widehat{\mathbf{x}}}_q\left(t_n\right|t_{n-1})& ={\mathbf{A}}_q{\widehat{\mathbf{x}}}_q\left(t_{n-1}\right|t_{n-1}),\hfill \\ \hfill {\mathbf{P}}_q\left(t_n\right|t_{n-1})& ={\mathbf{A}}_q{\mathbf{P}}_q\left(t_{n-1}\right|t_{n-1}){\mathbf{A}}_q^T+{\mathbf{Q}}_q,\hfill \\ \hfill \mathbf{K}\left(t_n\right)& ={\mathbf{P}}_q\left(t_n\right|t_{n-1}){\mathbf{H}}_q^T\left(t_n\right){\left({\mathbf{H}}_q\left(t_n\right){\mathbf{P}}_q\left(t_n\right|t_{n-1}){\mathbf{H}}_q^T\left(t_n\right)+{\mathbf{R}}_q\left(t_n\right)\right)}^{-1},\hfill \\ \hfill {\widehat{\mathbf{x}}}_q\left(t_n\right|t_n)& ={\widehat{\mathbf{x}}}_q\left(t_n\right|t_{n-1})+\mathbf{K}\left(t_n\right)\left({\widehat{\tau }}_q\left(t_n\right)-\mathit{h}\left({\widehat{\mathbf{x}}}_q\left(t_n\right|t_{n-1})\right)\right),\hfill \\ \hfill {\mathbf{P}}_q\left(t_n\right|t_n)& ={\mathbf{P}}_q\left(t_n\right|t_{n-1})-\mathbf{K}\left(t_n\right){\mathbf{H}}_q\left(t_n\right){\mathbf{P}}_q\left(t_n\right|t_{n-1}),\hfill \end{array}$$

(16)

where $\mathbf{K}\left(t_n\right)$ is the Kalman gain, ${\widehat{\mathbf{x}}}_q\left(t_n\right|t_{n-1})$ and ${\mathbf{P}}_q\left(t_n\right|t_{n-1})$ are the predicted state and covariance matrix at time $t_n$, and ${\widehat{\mathbf{x}}}_q\left(t_n\right|t_n)$ and ${\mathbf{P}}_q\left(t_n\right|t_n)$ are the corrected state estimate and covariance matrix after the measurement update. The local linearization of the measurement model around ${\widehat{\mathbf{x}}}_q\left(t_n\right|t_{n-1})$ is denoted in the observation matrix ${\mathbf{H}}_q\left(t_n\right)$. Since the elements of measurement vector ${\mathit{\tau }}_q\left(t_n\right)$ only depend on position and not on velocity, the linearized observation matrix consists of the Jacobian in Equation (10) at the predicted position ${\widehat{\mathbf{r}}}_q\left(t_n\right|t_{n-1})$ and of a $\left|\mathcal{P}\right|\times 2$ zero matrix $\mathbf{0}$, written as

$${\mathbf{H}}_q\left(t_n\right)=\left[\mathbf{J}\left({\widehat{\mathbf{r}}}_q\left(t_n\right|t_{n-1})\right),\mathbf{0}\right].$$

(17)

4. Performance Bounds

For evaluating estimators, typically, theoretical lower bounds on the estimation performance can be used. With regard to the proposed localization approach of this paper, this section first provides the classic CRLB on the error of the position estimate and second the PCRLB as performance bound for unbiased sequential Bayesian estimators. Thereby, the latter bound allows for evaluating recursive Bayesian filters like the proposed EKF.

4.1. Cramér–Rao Lower Bound on Position Estimation

The CRLB is defined as the inverse of the FIM [21]. This means, given the vector parameter ${\mathbf{r}}_q$, the elements of the unbiased estimator ${\widehat{\mathbf{r}}}_q={\left[{\widehat{x}}_q,{\widehat{y}}_q\right]}^T$ satisfy

$$\mathrm{Var}\left({\widehat{x}}_q\right)\ge {\left[\mathbf{F}{\left({\mathbf{r}}_q\right)}^{-1}\right]}_{1,1}={\sigma }_{x_q}^2$$

(18)

and

$$\mathrm{Var}\left({\widehat{y}}_q\right)\ge {\left[\mathbf{F}{\left({\mathbf{r}}_q\right)}^{-1}\right]}_{2,2}={\sigma }_{y_q}^2,$$

(19)

where $\mathrm{Var}\left(·\right)$ denotes the variance of an estimator and the terms ${\left[\mathbf{F}{\left({\mathbf{r}}_q\right)}^{-1}\right]}_{n,n}$, $n=\{1,2\}$, denote the diagonal elements of the inverse FIM $\mathbf{F}\left({\mathbf{r}}_q\right)$.

Inherently, received signals are a function of propagation delays ${\mathit{\tau }}_q$. This holds for every pair of transmitting and receiving nodes in the network. Thus, by applying the chain rule, the FIM $\mathbf{F}\left({\mathbf{r}}_q\right)$ can be alternatively expressed as [21]

$$\mathbf{F}\left({\mathbf{r}}_q\right)=\mathbf{J}{\left({\mathbf{r}}_q\right)}^T\mathbf{F}\left({\mathit{\tau }}_q\right)\mathbf{J}\left({\mathbf{r}}_q\right),$$

(20)

with the Jacobian matrix $\mathbf{J}\left({\mathbf{r}}_q\right)$ as defined in Equation (10) and the FIM $\mathbf{F}\left({\mathit{\tau }}_q\right)\in {\mathbb{R}}^{\left|\mathcal{P}\right|\times \left|\mathcal{P}\right|}$ with respect to delay vector ${\mathit{\tau }}_q$ (Please note also that path loss and phase information can be taken into account for calculating the FIM but is out of the scope of this work). For the linear independent time delays ${\mathit{\tau }}_q$, the Fisher information is well known [21] and the diagonal elements of $\mathbf{F}\left({\mathit{\tau }}_q\right)$ can be written as

$${\left[\mathbf{F}\left({\mathit{\tau }}_q\right)\right]}_{p,p}=\frac{8{\pi }^2{\beta }^2{\mathrm{SNR}}_{pq}}{c^2},$$

(21)

with ${\beta }^2$ as effective bandwidth of the transmit signal. The index $p=1,\cdots ,\left|\mathcal{P}\right|$ provides an enumeration of index set $\mathcal{P}$. Therewith, ${\mathrm{SNR}}_{pq}$ expresses the SNR for the MPC caused by scatterer q in network link p. With ${\mathrm{Tx}}_k$ and ${\mathrm{Rx}}_l$ defining link p, the SNR can be written as

$${\mathrm{SNR}}_{pq}=\left(P_{{\mathrm{Tx}}_k}\frac{G_{{\mathrm{Tx}}_k}{G}_{{\mathrm{Rx}}_l}{\sigma }_q{c}^2}{{\left(4\pi \right)}^3{f}_{\mathrm{c}}^2{\left(d_{kq}^{\mathrm{Tx}}\right)}^2{\left(d_{lq}^{\mathrm{Rx}}\right)}^2}\right)P_{\mathrm{n}}^{-1},$$

(22)

where $P_{{\mathrm{Tx}}_k}$ and $P_{\mathrm{n}}$ denote transmit power and receiver noise power, $G_{{\mathrm{Tx}}_k}$ and $G_{{\mathrm{Rx}}_l}$ express antenna gains at transmitter and receiver, and ${\sigma }_q$ refers to the RCS of scatterer q [7]. Hence, the elements of the FIM in Equation (21) strongly depend on the position of the scatterer and are proportional to

$${\left[\mathbf{F}\left({\mathit{\tau }}_q\right)\right]}_{p,p}\propto {\left({\left(d_{kq}^{\mathrm{Tx}}\right)}^2{\left(d_{lq}^{\mathrm{Rx}}\right)}^2\right)}^{-1}.$$

(23)

Finally, the CRLB of Equations (18) and (19) can be used to evaluate an estimator by comparing it to the RMSE according to inequality

$${\mathrm{RMSE}}_q=\sqrt{E\left[\parallel {\widehat{\mathbf{r}}}_q-{\mathbf{r}}_q{\parallel }^2\right]}\ge \sqrt{{\sigma }_{x_q}^2+{\sigma }_{y_q}^2}.$$

(24)

4.2. Posterior Cramér–Rao Lower Bound for Nonlinear Sequential Bayesian Estimation

The CRLB, as introduced above, allows for evaluating the positioning performance for a specific time instant. Thereby, the system is assumed to be time-invariant. An evaluation of Bayesian filtering and tracking approaches, however, requires a performance bound, which accounts for time-variant systems, for underlying stochastic state space models, and for the incorporation of prior knowledge [24]. Such a performance bound is provided by the PCRLB. The PCRLB is the theoretical performance bound for sequential Bayesian estimators [24,25]. Equivalently to the CRLB, the PCRLB is calculated by the inverse of the posterior FIM ${\mathbf{F}}_{B,q}\left(t_n\right)$ defining the inequality

$$E_{{\mathbf{x}}_q\left(t_n\right)}\left[\left({\widehat{\mathbf{x}}}_q\left(t_n\right)-{\mathbf{x}}_q\left(t_n\right)\right){\left({\widehat{\mathbf{x}}}_q\left(t_n\right)-{\mathbf{x}}_q\left(t_n\right)\right)}^T\right]={\mathbf{M}}_q\left(t_n\right)\ge {\mathbf{F}}_{B,q}{\left(t_n\right)}^{-1},$$

(25)

with $E_{\mathbf{a}}\left[·\right]$ as expectation with respect to the probability density of random variable $\mathbf{a}$. The inequality in Equation (25) implies that the difference between the MSE matrix ${\mathbf{M}}_q\left(t_n\right)$ and the inverse of the posterior FIM is a positive semi-definite matrix. Following [24,25,26,27], the posterior FIM ${\mathbf{F}}_{B,q}\left(t_n\right)$ can be calculated recursively as

$${\mathbf{F}}_{B,q}\left(t_n\right)={\mathbf{D}}_{22,q}\left(t_n\right)-{\mathbf{D}}_{21,q}\left(t_n\right){\left({\mathbf{F}}_{B,q}\left(t_{n-1}\right)+{\mathbf{D}}_{11,q}\left(t_n\right)\right)}^{-1}{\mathbf{D}}_{12,q}\left(t_n\right),$$

(26)

with

$$\begin{array}{cc}\hfill {\mathbf{D}}_{11,q}\left(t_n\right)=& E_{{\mathbf{x}}_q\left(t_{n-1}\right),{\mathbf{x}}_q\left(t_n\right)}\left[-{\Delta }_{{\mathbf{x}}_q\left(t_{n-1}\right)}^{{\mathbf{x}}_q\left(t_{n-1}\right)}lnp\left({\mathbf{x}}_q\left(t_n\right)|{\mathbf{x}}_q\left(t_{n-1}\right)\right)\right],\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\mathbf{D}}_{12,q}\left(t_n\right)=& E_{{\mathbf{x}}_q\left(t_{n-1}\right),{\mathbf{x}}_q\left(t_n\right)}\left[-{\Delta }_{{\mathbf{x}}_q\left(t_{n-1}\right)}^{{\mathbf{x}}_q\left(t_n\right)}lnp\left({\mathbf{x}}_q\left(t_n\right)|{\mathbf{x}}_q\left(t_{n-1}\right)\right)\right]={\mathbf{D}}_{21,q}{\left(t_n\right)}^T,\hfill \\ \hfill {\mathbf{D}}_{22,q}\left(t_n\right)=& E_{{\mathbf{x}}_q\left(t_{n-1}\right),{\mathbf{x}}_q\left(t_n\right)}\left[-{\Delta }_{{\mathbf{x}}_q\left(t_n\right)}^{{\mathbf{x}}_q\left(t_n\right)}lnp\left({\mathbf{x}}_q\left(t_n\right)|{\mathbf{x}}_q\left(t_{n-1}\right)\right)\right]\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & +E_{{\mathbf{x}}_q\left(t_n\right)}\underset{\mathrm{classic}\mathrm{FIM}}{\underbrace{\left[E_{{\mathbf{z}}_q\left(t_n\right)|{\mathbf{x}}_q\left(t_n\right)}\left[-{\Delta }_{{\mathbf{x}}_q\left(t_n\right)}^{{\mathbf{x}}_q\left(t_n\right)}lnp\left({\mathbf{z}}_q\left(t_n\right)|{\mathbf{x}}_q\left(t_n\right)\right)\right]\right]}}.\hfill \end{array}$$

(29)

With ${\nabla }_{\mathbf{a}}$ denoting the partial derivatives with respect to $\mathbf{a}$, the operator ${\Delta }_{\mathbf{b}}^{\mathbf{a}}={\nabla }_{\mathbf{b}}{\nabla }_{\mathbf{a}}^T$ gives the corresponding second-order partial derivatives. Furthermore, $p\left(\mathbf{a}\left(t_n\right)|\mathbf{b}\left(t_n\right)\right)$ is the conditional probability distribution of random variable $\mathbf{a}\left(t_n\right)$ given $\mathbf{b}\left(t_n\right)$ at time instant $t_n$. Note that the term inside the second expectation of ${\mathbf{D}}_{22,q}\left(t_n\right)$, as highlighted by a brace, refers to the definition of the classic FIM [25]. For state ${\mathbf{x}}_q\left(t_n\right)$ given in Equation (11), the $4\times 4$ FIM is

$$\tilde{\mathbf{F}}\left({\mathbf{x}}_q\left(t_n\right)\right)=\left[\begin{array}{cc}\mathbf{F}\left({\mathbf{r}}_q\left(t_n\right)\right)& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right],$$

(30)

where $\mathbf{F}({\mathbf{r}}_q\left(t_n\right)$ denotes the FIM of Equation (20). Thus, referring to the linear transition model introduced in Section 3.3 with transition matrix ${\mathbf{A}}_q$ given in Equation (13) and white Gaussian process noise with covariance matrix ${\mathbf{Q}}_q$ given in Equation (14), Equations (27)–(29) result in

$$\begin{array}{cc}\hfill {\mathbf{D}}_{11,q}\left(t_n\right)=& {\mathbf{A}}_q^T{\mathbf{Q}}_q^{-1}{\mathbf{A}}_q,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\mathbf{D}}_{12,q}\left(t_n\right)=& -{\mathbf{A}}_q^T{\mathbf{Q}}_q^{-1}={\mathbf{D}}_{21,q}{\left(t_n\right)}^T,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\mathbf{D}}_{22,q}\left(t_n\right)=& {\mathbf{Q}}_q^{-1}+E_{{\mathbf{x}}_q\left(t_n\right)}\left[\tilde{\mathbf{F}}\left({\mathbf{x}}_q\left(t_n\right)\right)\right].\hfill \end{array}$$

(33)

Inserting Equations (31)–(33) into Equation (26) and applying the matrix inversion lemma results in the recursive expression

$${\mathbf{F}}_{B,q}\left(t_n\right)={\left({\mathbf{Q}}_q+{\mathbf{A}}_q{\mathbf{F}}_{B,q}{\left(t_{n-1}\right)}^{-1}{\mathbf{A}}_q^T\right)}^{-1}+E_{{\mathbf{x}}_q\left(t_n\right)}\left[\tilde{\mathbf{F}}\left({\mathbf{x}}_q\left(t_n\right)\right)\right].$$

(34)

Given prior information according to the probability density function $p\left({\mathbf{x}}_q\left(t_0\right)\right)$, the initial FIM ${\mathbf{F}}_{B,q}\left(t_0\right)$ can be calculated as

$${\mathbf{F}}_{B,q}\left(t_0\right)=E_{{\mathbf{x}}_q\left(t_0\right)}\left[-{\Delta }_{{\mathbf{x}}_q\left(t_0\right)}^{{\mathbf{x}}_q\left(t_0\right)}lnp\left({\mathbf{x}}_q\left(t_0\right)\right)\right],$$

(35)

which is used to initialize the recursion in Equation (34).

5. Case Study

In this section, the proposed localization approach is analyzed for the example of a walking pedestrian. Together with the network structure, the measurement setup is presented first. Based on the introduced static network, the localization performance is evaluated theoretically using the performance bounds of Section 4. Finally, the localization approach is applied to channel measurements.

5.1. Network and Measurement Setup

The analyzed static network consists of $K=1$ transmitting and $L=4$ receiving nodes. As shown in Figure 1a, the network nodes are individually placed with an inter-node distance between 10 $\mathrm{m}$ and 15 $\mathrm{m}$. The center of the network forms a transmitting node with a small directional antenna. Thereby, the main beam of the transmit antenna is oriented towards a designated observation area between the receiving nodes, as shown in Figure 1b, illustrating the gain of the transmit antenna in the network area [28]. With the transmit antenna $\mathrm{Tx}$ chosen as the origin of the coordinate system, the overall considered observation area spans from $-7.5$ m to 30 m in the x-direction and from $-10$ m to 26 m in the x-direction. The selected scenario considers a walking pedestrian as single mobile scatterer ($Q=1$), walking a wide circle within the illuminated area of the static network setup. As the experiment was conducted on an apron close to surrounding buildings, an elevated tachymeter was used as ground truth system. Using the tachymeter together with a high precision reflector prism, the individual static antenna locations of the network were determined prior to the experiment. For recording the ground truth during the experiment, the pedestrian wore a helmet, on which the reflector prism was mounted.

As a measurement system, the Medav RUSK-DLR wideband channel sounder was used [29]. Therefore, the measured data are CIR between transmit and receive antennas, i.e., network nodes, respectively. A full summary of the corresponding measurement parameter settings provides

Table 1. Measurement parameters.

Parameter	Value
Center frequency $f_{\mathrm{c}}$	$5.2$ GHz
Bandwidth B	120 MHz
Signal period $T_{\mathrm{p}}$	$3.2$ μs
Measurement rate $T_{\mathrm{g}}$	$2.048$ ms
Transmit power $P_{\mathrm{Tx}}$	37 dBm
Antenna gain $G_{\mathrm{Tx}}$	9 dBi (small directional [28])
Antenna gain $G_{{\mathrm{Rx}}_{1-4}}$	8 dBi (toroidal, omni-directional)

.

5.2. Theoretical Performance Evaluation

For illustrating the localization performance of the measurement setup, the RMSE is calculated for the considered observation area according to the CRLB as defined in Equation (24). To determine the RMSE, the parameters provided in Table 1 are used. Assuming a transmit signal with rectangular power spectral density, the effective bandwidth is defined as ${\beta }^2=B^2/12$ [21]. The value of the RCS accounts for an object’s reflectivity characteristics, which is influenced by its size, shape, and material. Thus, for a pedestrian being considered as moving object, a typical RCS is 1 $\mathrm{m}$2 [7]. Eventually, the localization RMSE for the measurement setup is shown in Figure 2a. Overall, the shape of the RMSE reveals good localization performance for scatterers located in areas of high transmit antenna gain; see also Figure 1b. Hence, scatterers located close to the transmitting node can be localized very precisely. Qualitatively, this holds also for scatterers located close to receiving nodes. For these areas, the high localization performance can be explained by the received signal strength of the backscattered signal, since the RMSE strongly depends on the SNR. Additionally, the decreasing localization performance for far positions from the network nodes confirms this observation. Besides SNR, the localization performance depends on the system geometry. Effects due to system geometry can be explained by interpreting the scatterer location as the intersection point of multiple ellipses, which correspond to the propagation delays of respective transmitting and receiving nodes. To achieve high localization performance in both x- and y-directions, ideally, ellipses should intersect perpendicular to each other. For scatterers located far from the localization network, the shapes of the individual ellipses differ only marginally from each other. Thus, the intersection angle is very low, which results in an increased localization uncertainty in the direction parallel to the tangent at the intersection point. With a widely distributed network covering a large observation area, the effects of low intersection angles can be avoided. A further geometrical effect can be observed for scatterers located in the proximity of the baseline between transmitting and receiving nodes. Regarding a single link, the location information in the direction parallel to the baseline is very low, and, thus, the localization performance also decreases. Accordingly, in the localization network, the localization performance decreases in the proximity of the individual baselines; see Figure 2a. This performance degradation is independent of the network topology.

Please note that the derived CRLB of Section 4 only depends on waveform and SNR. This means that the influence of any superposition of LoS and MPC on the parameter estimation is not considered [27], and is out of the scope of this paper. Since the superposition of LoS and MPCs strongly impacts the estimation capabilities of its influence on parameter estimation will be included in future research.

In order to evaluate the performance of the EKF proposed in Section 3.3, the PCRLB is calculated for different scenarios within the measurement setup and compared to simulation results of the tracking filter. Three different scenarios are considered, each for a single scatterer moving with constant velocity. The initial absolute velocity is per-mode=symbol $1.41$ $\mathrm{m}$/$\mathrm{s}$ for every scenario. However, each scenario possesses an individual starting position and movement direction. Figure 2b provides an overview over the considered scenarios. Referring to Figure 2a, Scenario $\mathrm{I}$ is characterized by high localization capabilities throughout the whole trajectory. Scenario $\mathrm{II}$, in contrast, crosses an area of poor localization capabilities between $\mathrm{Tx}$ and ${\mathrm{Rx}}_3$. In Scenario $\mathrm{III}$, the trajectory starts in an area of poor localization capabilities and moves towards an area of very high localization capabilities in the main beam the transmitting node. For the calculation of the PCRLB as well as the EKF simulations, similar system and signal parameters as for calculating the static positioning CRLB are used. As process noise intensity for the covariance in Equation (14), a value of ${\sigma }_q^2=$ per-mode=symbol $0.01$ $\mathrm{m}$${}^2$/$\mathrm{s}$3 is defined [23]. The filter state is initialized randomly according to the initial state covariance $\mathbf{P}\left(t_0\right)$ around the initial state $\mathbf{x}\left(t_0\right)$ given by the respective scenario. Thereby, $\mathbf{P}\left(t_0\right)$ is defined as $4\times 4$ diagonal matrix, with an initial variance of per-mode=symbol $0.1$ $\mathrm{m}$${}^2$ on position and per-mode=symbol $0.01$ $\mathrm{m}$${}^2$/$\mathrm{s}$${}^2$ on velocity, in both x- and y-directions. Determining the PCRLB and the EKF performance results requires performing multiple Monte Carlo runs. In each run, the system equations of Section 3.3 are simulated with different samples of the process noise for 10 $\mathrm{s}$, which equals an average walking distance of approximately 14 $\mathrm{m}$. Thus, the multitude of simulation runs allows for approximating the expectation in Equation (33). In order to estimate the MSE matrices in Equation (25), the measurement noise also needs to be sampled for each run. In this study, for every scenario, 5000 Monte Carlo runs are performed with 50 realizations of measurement noise. Hence, the EKF is evaluated $2.5\times {10}^5$ times in each scenario. Thereby, the number of Monte Carlo runs was chosen to achieve statistically stable results for the PCRLB and to achieve results for the EKF, which fluctuate only marginally compared to the absolute RMSE values.

Figure 3 shows the simulation results for the three scenarios. The results for both PCRLB and EKF are given in terms of RMSE, i.e., as the square root of the MSE in Equation (25). The RMSE of PCRLB and EKF are denoted by ${ϵ}_{PCRLB}$ and ${ϵ}_{EKF}$. Given $\mathbf{P}\left(t_0\right)$, the initial RMSE of the positioning error is approximately $0.45$ $\mathrm{m}$ for all scenarios. Due to sequential filtering characteristics, the localization RMSE strongly decreases shortly after the initialization and follows the static performance capabilities along the respective scenario trajectory. Overall, it can be observed that the filter RMSE values are very close to those of the theoretical bound. This holds for the simulation results of each scenario. Particularly, for Scenario $\mathrm{I}$, the results of the EKF converge very fast and reach the bound in less than 1 $\mathrm{s}$; see Figure 3a. In addition, for Scenario $\mathrm{II}$, the EKF converges to the bound after about 3 $\mathrm{s}$ as shown in Figure 3b. With the trajectory crossing an area of poor localization capabilities between $\mathrm{Tx}$ and ${\mathrm{Rx}}_3$, both the PCRLB and the filter RMSE show a temporary increase. For scatterer positions close to the baseline between transmitting and receiving nodes, the elements of the Jacobian in Equation (10) and therewith of the observation matrix in Equation (17) are close to zero. Due to singularities caused by the inversion of the observation matrix, the region close to the baselines between $\mathrm{Tx}$ and ${\mathrm{Rx}}_3$, as well as $\mathrm{Tx}$ and ${\mathrm{Rx}}_4$, is particularly challenging for the filter. For calculating the bound in Equation (34), the expectation of the FIM averages these singularities. The results in Figure 3b confirm this effect, with the filter RMSE diverging from the bound between 1 $\mathrm{s}$ and 2 $\mathrm{s}$. Finally, the filter results for Scenario $\mathrm{III}$ are provided in Figure 3c. The filter RMSE shows a fast convergence to the bound after approximately 2 $\mathrm{s}$. As described above, however, the effects of singularities on the EKF in regions close to $\mathrm{Tx}$-$\mathrm{Rx}$ baselines also influences the filter results of the third scenario. Here, the mean trajectory crosses $\mathrm{Tx}$-${\mathrm{Rx}}_4$ and travels along $\mathrm{Tx}$-${\mathrm{Rx}}_3$. With increasing state covariance, many simulation trajectories lie on or are close to the baselines. Thus, the filter RMSE does not fully converge to the bound in the period after 2 $\mathrm{s}$ of simulation time.

5.3. Measurement Based Evaluation

As stated in Section 3, the channel parameters for all links in the network are estimated using KEST. The estimation results of KEST for the measured CIR over time are shown in Figure 4a–d for each link individually. The figures show the consecutive vectors of delay estimates in Equation (5) over the full measurement time. Static delay estimates, including LoS and static MPCs, are shown in gray. These static delays together with corresponding standard deviations are determined during a preceding calibration phase according to the procedures described in Section 3.1. This characterization of the static propagation environment allows for extracting the mobile MPC, highlighted in color according to the estimated amplitude level. For each link, the delay estimates fit the ground truth data as indicated by black dashed lines. Hence, the results qualitatively confirm the point scattering assumption for pedestrians [29]. Due to the geometrical arrangement of the receiving antennas and limits in the dynamic range of the measurement system, the detection, estimation, and tracking capabilities of KEST differ. With the small directional antenna gain and the orientation of the transmitting antenna, as shown in Figure 1b, the LoS signal power is reduced for all links. The reduced LoS signal power avoids an elevation of the noise floor due to limits of the dynamic range. Noticeably, link $\mathrm{Tx}$-${\mathrm{Rx}}_3$ exhibits particularly good parameter estimation throughout the measurement. On the one hand, this can be explained by the advantageous placement of the receive antenna with respect to the transmit antenna gain, i.e., very low $\mathrm{Tx}$ gain towards ${\mathrm{Rx}}_3$. On the other hand, ${\mathrm{Rx}}_3$ and also ${\mathrm{Rx}}_4$ are not impacted by strong multipath fading like ${\mathrm{Rx}}_1$ and ${\mathrm{Rx}}_2$. Since ${\mathrm{Rx}}_1$ and ${\mathrm{Rx}}_2$ are located closer to the hangar wall, signal reflections off the wall are received with higher power. Examples for such multipath fading provide the estimation results of ${\mathrm{Rx}}_1$ and ${\mathrm{Rx}}_2$, where double reflections from the pedestrian and the hangar wall are clearly visible (see Figure 4a,b). However, apart from the perturbing effect regarding parameter estimation, assignable double reflections from mobile scatterers would contain additional location information. Even though in this work double reflections are not considered as source of information, an exploitation of reflected mobile MPC will be included in future research. Besides dynamic range and multipath fading, the estimation results of Figure 4 clearly show how the presence of LoS and static MPCs impact the composition of measurement vector ${\widehat{\mathit{\tau }}}_q$ and therewith the localization of the mobile scatterer. An unambiguous solution of Equation (8) requires at least three independent measurements. The required static environment mitigation, i.e., the displacement of any delay estimates close to static components (see Section 3.2), however, reduces the amount of delay estimates assigned to a specific mobile scatterer. Thus, a rich static MPC environment reduces the overall localization capabilities with the current approach. This holds particularly for very sparse networks, such as the four-link network considered in this paper, since an outage in two links impedes localization. Besides static MPCs, the LoS also impacts the localization capability due to the so-called blind zone problem [30]. This means that MPC caused by scatterers located close to the baseline between a transmitting and receiving node is hardly detectable. Exemplary, the parameter estimation results for link $\mathrm{Tx}$-${\mathrm{Rx}}_4$ confirm these blind zones, since it is not possible to extract mobile MPCs when the pedestrian is located in the proximity to the link baseline. Thus, for ${\mathrm{Rx}}_4$, no parameter estimates are available from 33 $\mathrm{s}$ to 44 $\mathrm{s}$.

Finally, the localization approach proposed in Section 3 is used to estimate and track the location of the walking pedestrian using the CIR measurements described above. For initialization, the iterative localization procedure of Equation (9) is applied to determine an initial location estimate based on available parameter estimates of the four links at $t_0=0$ s. The initial velocity components in x- and y-directions are randomly chosen and taken from a uniform distribution $\mathcal{U}(0 \mathrm{m}/\mathrm{s},1 \mathrm{m}/\mathrm{s})$, respectively. Given the initial state, the EKF described in Section 3.3 is used for tracking the location of the pedestrian. According to [23], the covariance matrix of the process noise in Equation (14) is determined by the process noise intensity of ${\sigma }_q^2=0.01 {\mathrm{m}}^2/{\mathrm{s}}^3$. The resulting positioning error for the experiment over time is given in Figure 5 and the corresponding CDF in Figure 6. Overall, except for a few outliers, the positioning error remains below 1 $\mathrm{m}$ throughout the experiment. In particular, according to the CDF, more than 97% of positioning errors are below 1 $\mathrm{m}$. The first outlier appears between 20 $\mathrm{s}$ and 25 $\mathrm{s}$. This time period coincides with a sharp left curve, as shown in Figure 1a. Thus, a delayed adaption of the velocity vector by the EKF can be a possible explanation for the exceeding positioning error. A very strong localization performance can be observed between 32 $\mathrm{s}$ and 37 $\mathrm{s}$, when the pedestrian passes the main beam of the transmit antenna. Due to strong reflection and therewith high SNR, the parameters of each link and thus the position can be estimated very well. For the period from 37 $\mathrm{s}$ to 44 $\mathrm{s}$, the error steadily increases. In this period, parameter estimates are almost solely available for ${\mathrm{Rx}}_3$. Hence, the few simultaneous and even erroneous delay estimates can explain the increasing positioning error and the outlier at about 44 $\mathrm{s}$. After leaving the blind zone of $\mathrm{Tx}$-${\mathrm{Rx}}_4$ at about 44 $\mathrm{s}$, the additional link parameter estimates support the positioning performance thereafter. The increasing error after 53 $\mathrm{s}$ is again due to few simultaneous parameter estimates from different network links. Similarly to before, the position estimation is mainly driven by parameter estimates for ${\mathrm{Rx}}_3$.

6. Conclusions

This paper presented a localization approach to detect, locate, and track moving, non-cooperative objects by means of a network of spatially distributed transmitting and receiving nodes. With moving objects affecting the radio spectrum as time-variant MPC, i.e., discrete mobile scatterers, multipath propagation could be exploited for passive localization. Therefore, the localization of mobile scatterers was formulated as a nonlinear optimization problem. For initialization, the proposed approach uses an iterative nonlinear least squares algorithm following Levenberg and Marquardt to solve the optimization problem. Subsequently, an EKF is applied to estimate the scatterer location recursively over time. For both the snapshot based localization and the nonlinear sequential Bayesian estimation problem performance bounds, the classic CRLB on position estimation and the PCRLB were derived.

The proposed approach is evaluated based on a case study, which considers a pedestrian as a single mobile scatterer walking within a network of one transmitting and four receiving nodes. Thereby, a simulation study has shown that the EKF achieves a localization performance very close to the PCRLB. Even though the EKF performance results diverged from the theoretical bound in areas close to $\mathrm{Tx}$-$\mathrm{Rx}$ baselines due to singularities, the study overall confirmed the applicability of the EKF for the localization problem. Particularly, for scatterers remote from any $\mathrm{Tx}$-$\mathrm{Rx}$ baseline, the localization performance of the filter was shown to converge to the PCRLB. Moreover, the approach was applied to wideband measurement data corresponding to the case study of the theoretical analysis. With a resulting localization accuracy of $0.8$ m at 90% confidence, the case study proved that mobile scatterers can be localized using CIR, i.e., exploiting multipath propagation.

The performance of the proposed approach strongly depends on the underlying parameter estimation algorithm and the mitigation of the LoS and static MPC. Since both the parameter estimation and the mitigation of the static components are prone to errors, future work will focus on a direct multilateration approach. The direct usage of CIR for localization and tracking could help to overcome the dependence on individual parameter estimation and would link an object’s mobility information implicitly to the recursive position estimation.