Detection and Tracking of Human Targets using Ultra-Wideband Radar

(1)

UPTEC F16 058

Examensarbete 30 hp

Oktober 2016

Detection and Tracking of Human

Targets using Ultra-Wideband Radar

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten Besöksadress: Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress: Box 536 751 21 Uppsala Telefon: 018 – 471 30 03 Telefax: 018 – 471 30 00 Hemsida: http://www.teknat.uu.se/student

Abstract

Detection and Tracking of Human Targets using

Ultra-Wideband Radar

Andreas Östman

The purpose of this thesis was to assess the plausibility of using two

Ultra-Wideband radars for detecting and tracking human targets. The detection has been performed by two different types of methods, constant false-alarm rate methods and a type of CLEAN algorithm. For tracking the targets, multiple hypothesis tracking has been studied. Particle filtering has been used for the state prediction, considering a significant amount of uncertainty in a motion model used in this thesis project. The detection and tracking methods have been implemented in MATLAB. Tracking in the cases of a single target and multiple targets has been investigated in simulation and experiment. The simulation results in these cases were compared with accurate ground truth data obtained using a VICON optical tracking system.

The detection methods showed poor performance when using data that had been collected by the two radars and post-processed to enhance target features. For single targets, the detections were accurate enough to continuously track a target moving randomly in a controlled area. In the multiple target cases the tracker was not able to distinguish the multiple moving subjects.

ISSN: 1401-5757, UPTEC F16 058 Examinator: Tomas Nyberg Ämnesgranskare: Ping Wu Handledare: Alen Alempijevic

(3)

Detection and Tracking of Human Targets using

Ultra-Wideband Radar

Andreas Östman October 17, 2016

Abstract

The purpose of this thesis was to assess the plausibility of using two Ultra-Wideband radars for detecting and tracking human targets. The detection has been performed by two different types of methods, constant false-alarm rate methods and a type of

CLEAN algorithm. For tracking the targets, multiple hypothesis tracking has been

studied. Particle filtering has been used for the state prediction, considering a significant amount of uncertainty in a motion model used in this thesis project. The detection and tracking methods have been implemented in MATLAB. Tracking in the cases of a single target and multiple targets has been investigated in simulation and experiment. The simulation results in these cases were compared with accurate ground truth data obtained using a VICON optical tracking system.

The detection methods showed poor performance when using data that had been collected by the two radars and post-processed to enhance target features. For single targets, the detections were accurate enough to continuously track a target moving ran-domly in a controlled area. In the multiple target cases the tracker was not able to distinguish the multiple moving subjects.

(4)

Populärvetenskaplig sammanfattning

Detektion och lokalisering av människor används för mängder av olika ändamål. Vanligen görs detta med hjälp av visuell spårning som baseras på optisk teknik såsom kameror. På senare år har dock alternativa tekniker fått mer uppmärksamhet i lokaliseringssammanhang, varav en av dessa tekniker är att använda en särskild sorts RADAR med väldigt stor bandbredd. Den ökade bandbredden ger bra upplösning i avståndsmätningar, och i samband med att tekniken blir billigare har intresset för den här typen av RADAR ökat markant.

Att uppskatta objekts positioner med RADAR består av två delar: Detektion och Track-ing. Detektion är processen som görs på rå radardata och går ut på att räkna ut uppskattade avstånd till objekt baserat på radarsignalens propageringstid. Detektionsproceduren görs separat för varje RADAR som används. Efter att avstånd har uppskattats används dessa i en Tracker för att över tid följa rörliga objekts positioner över en känd karta.

I det här examensarbetet har ett antal detektions- och trackingalgoritmer implementerats för att appliceras på RADAR-data som samlats in från två enheter. För att kunna samla in radardatan implementerades även ett API som behövdes för att kommunicera mellan de båda enheterna och en dator. Både fixa och adaptiva gränsvärdesmetoder har undersökts för detektion. För tracking har en Multipel-hypotesmetod används, vars idé baseras på att hålla hypoteser kring potentiella mål kvar tills fler mätvärden tas emot som kan förkasta eller bekräfta skapade målhypoteser. Allt har implementerats i MATLAB för att snabbt kunna skapa en prototyp av detektions- och trackingsystemet. Det skapade systemet testades på riktig data där fall med allt från ett till flera rörliga mål undersöktes. Det visade sig att detektionsmetoderna var flaskhalsen i systemet givet den rådatan som hade samlats in, då kvalitéten på denna inte var tillräcklig för att kunna upptäcka mål i fler-målsfallen under nämnvärda perioder under den tid som trackingen undersöktes.

(5)

List of Figures

1.1 Framework overview . . . 5

2.1 UWB radar working principle. . . 6

2.2 PulsON 410 unit . . . 8

3.1 Example of multipaths . . . 10

3.2 Example of a Radargram . . . 11

4.1 CFAR cells . . . 14

4.2 CFAR schematic . . . 15

4.3 Example of CFAR applied to signal . . . 16

4.4 ROC for P_{F A} = 1 × 10−4, m = 1, 2 . . . 19 4.5 ROC for PF A = 1 × 10−4, m = 3, 4 . . . 20 5.1 MHT logic overview . . . 22 5.2 Tree structure of MHT . . . 22 5.3 Validation gating . . . 23 5.4 Track graph . . . 25

5.5 Markov models comparison . . . 27

5.6 Particle filtering overview . . . 29

5.7 System detailed view . . . 30

6.1 Main MATLAB classes . . . 33

6.2 Radar-Host communication . . . 33

6.3 Post-processing enhancements . . . 35

6.4 Measurement-to-Measurement association issue . . . 37

6.5 Example of adjacency matrix . . . 38

6.6 Tracking ambiguites . . . 41

6.7 N-scan pruning . . . 42

7.1 Experiment setup . . . 43

8.1 Detection results for one target. . . 55

8.2 Tracking results for one target . . . 56

8.3 Error comparisons for one target . . . 57

8.4 Detection results for two targets. . . 58

(8)

List of Tables

1 Notations . . . 1

2 Abbreviations and acronyms . . . 3

3 Single-target detection and tracking parameters . . . 36

4 Radar specifications and settings used in experiments . . . 44

(9)

Table 1: Notations ξ Received / Measured signal PF A Probability of false alarm

PD Probability of detection

PF C Probability of false censoring

ttof Time-of-flight

B Signal bandwidth

fh Upper frequency limit of signal bandwidth

fl Lower frequency limit of signal bandwidth

∆r Range resolution σ Standard deviation σ2 Variance

x Target state

Pk|k filter state covariance given observations up to time k.

nq Process noise

n Measurement noise

R Measurement noise covariance matrix Q Process noise covariance matrix N Number of particles (in particle filter). w Particle weights

Nlead Leading reference window length (CFAR)

Nlag Lagging reference window length (CFAR)

M Cell-under-test length (CFAR) α CFAR constant

T CFAR treshold

Nclean Neighbouring zero-padding length (CLEAN)

z Observation / Measurement

Σl_k Innovation covariance of track l and time k C(r) Measurement matrix

H(r) Jacobian of C(r)

Gth Innovation gate treshold

d2 Mahalanobis distance (squared) LR Likelihood ratio

(10)

LLR Log-likelihood ratio

S_kl The lth track LLR score at time k H1 Hypothesis of detecting a true target

H₀ Hypothesis of detecting a false target A Adjacency matrix

Vspace Measurement space

G Graph corresponding to entries of A V Set of all vertices of G

v Single vertex in V

E Set of all edges (connections between vertices in G) xk x-position of target at time k

yk y-position of target at time k

˙

xk x-component of velocity of target at time k

˙

yk y-component of velocity of target at time k

P (.) Probability of (.)

(11)

Table 2: Abbreviations and acronyms CFAR Constant False-Alarm rate

API Application Programming Interface SVD Singular-Value Decomposition PCA Principal-Component Analysis LADAR LAser Detection And Ranging

GCMLD Generalized Censored Mean Level Detector MHT Multiple Hypothesis Tracking

TO-MHT Track-oriented Multiple Hypothesis Tracking LLR Log-likelihood ratio

MOTA Multiple Object Tracking Accuracy MOTP Multiple Object Tracking Precision MWIS Maximum Weighted Independent Set EKF Extended Kalman Filter

(12)

1 Introduction

1.1 Background

Detection and tracking of human targets are useful for numerous applications. One of the commonly-used techniques for the detection and tracking is LADAR (Laser Detection And Ranging) that makes use of visual tracking. These techniques, however, are of very limited use if visible conditions are poor. An alternative to conventional approaches to tracking is to use e.g. infra-red imaging or radar for the same purposes. In particular a certain kind of radar called ultra wideband (UWB) radar has gained increased interest for tracking purposes. UWB radar is a type of radar that operates over very wide frequency bands. The wide frequency characteristics make it promising as it provides very good resolution in range [1]. Due to the small size and low power transmissions of UWB radars, multiple sensors can cost efficiently be positioned to cover limited areas. As the transmitted pulses of UWB radar is below the noise floor of most common commercial radio-frequency products, it proves to be suitable for indoor survailance systems without causing interference while also being harmless to people. UWB has already proven to be useful for tracking using range-only measurements in both controlled and more complex environments [2, 3, 4]. A complicated issue on multi-target detection and tracking is the ambiguities that arise when there are many successive false alarms or missed detections. Multi-sensor-multi-target tracking also gives rise to complex target association problems of determining measurement-to-measurement correspondence, which, using range-only measurements, makes any feature extraction for this association non-trivial. Different approaches to these issues have been suggested. In [5] a combination of Kalman filtering and expectation maximization (EM) shows promising results, however, the number of targets needs to be known a priori. One of the more novel techniques for tracking is to use probability hypothesis density (PHD) filters [6]. In [3] the author uses UWB radar with a constant false alarm rate method for detection and Gaussian mixture-PHD for tracking. The technique shows good accuracy for both one and two target scenarios that are presented in the present work. In [7] the association problem is addressed by attempting to extract slow-time features of the targets by using video time density functions (VTDF). This does, however, require knowledge of the Doppler shifts of the targets.

1.2 Purpose and goals

The goal of this thesis is to assist an existing framework (used for visual tracking) using depth cameras in situations where external factors might prevent them from being used. Furthermore, the detection and tracking was thought to be used as an estimator of crowdness in different environments.

1.3 Tasks and methodology

For acheiving the goal of the thesis, a well-adopted technique in the tracking community called Multiple Hypothesis Tracking [8] is considered with the goal of naturally resolving some of the data association issues that arise. Two pulse UWB radars are used for data aquisition and several detection techniques are reviewed and implemented. A flowchart of the different parts of the proposed detection and tracking framework can be seen in Figure 1.1. In the system overview, the data is collected by two radar modules, indicated by the arrow to the node “Signal post-processing”. In this node the collected scans are modified to enhance scan information regarding potential targets. The processed scan data is then forwarded to the “Detection”-node where different detection algorithms are applied and estimated range to targets are computed. The estimated range data from the two radars are then input to the tracker, represented by the node “Tracking” where attempts to predict the trajectory of the supposed targets are performed. The output of the tracking is a set of proposed trajectories

(13)

over time (in 2D) for the targets walking in the surveilled area. The scope of this thesis covers all sub-parts shown in the figure.

Radar measurements Signal post-processing (Target enhance-ment) Detection Tracking User output (2D target trajectories) Figure 1.1: High level overview of the proposed framework.

1.4 Outline of the thesis

The thesis starts with basic knowledge about radar where details are explained for UWB system types. In Section 4 the different detection algorithms that have been implemented are explained in detail. Section 5 covers the target tracking framework and related algo-rithms. Details regarding the experiments performed and results are found in Section 7 and 8 respectively. Lastly, conclusions are drawn and discussion of presented results is found in Section 9.

2 Radar

Radar (RAdio Detection and Ranging) is an object-detection and ranging system that uses radio waves to detect objects and determine their range, angle, or velocity. It usually consists of a transceiver (transmitter/receiver) and one or two antennas to generate electromagnetic radio waves and capture returns from targets. The distance r to any target can be computed using the round-trip time-of-flight of the radio signal as follows

r = cttof

2 (2.1)

where c is the speed of electromagnetic signal, t_tof is the time of flight during which the wave “flies” from transmitter to the target and then back the receiver. The factor 1₂ is due to the round-trip. A very simplistic visualization of the basic principle of a radar is shown in Figure 2.1.

(14)

(a) Bistatic radar. A pulse is transmitted by one antenna on one of the units. The pulse is reflected or pass through targets and reaches a spatially separated receiving antenna of a second unit.

(b) Monostatic radar. A pulse is transmitted by one antenna on the unit, reflected by the target and received by another antenna on the same unit that transmitted the pulse.

Figure 2.1: Simplified working principle of the UWB pulse radar. One unit transmits pulses which are reflected by objects in the surveilled area and received by the same (monostatic) or other (bistatic/multistatic) radar(s).

2.1 Ultra Wideband Radar

2.1.1 Characteristics

A system can be of different types in terms of bandwidth, e.g., narrow band, wideband, and ultra wideband (UWB). The type of a system can be defined based on the fractional bandwidth of the system [1], given by

Bf = 2

fh− fl

fh+ fl

(2.2) where f_l and f_h are the lower and upper cutoff frequencies of the transmitted signal, respec-tively. As of standards set by IEEE [9], a system is said to be of UWB type if Bf ≥ 0.2.

Common operating ranges for ultra wideband radar typically lies within the range of 100 MHz to 10 GHz [10].

For pulse radars transmitting sinusoidal pulses of duration τ (which is the case of the radars considered in this thesis), the -3dB bandwidth can be approximated as B = 1_τ. The range resolution ∆r depends on the pulse duration according to ∆r = cτ₂ [1]. Using the approximation of the bandwidth, the range resolution can be expressed as

∆r = c

2B (2.3)

which indicates that resolution improves with increased bandwidth. For a bandwidth of several GHz the range resolution is in the order of centimeters.

As UWB devices operate over wide bands electromagnetic energy is radiated over po-tentially occupied spectral bands. To avoid interference with existing systems the UWB radiation for commercial product is limited to a power spectral density of -41.3 dBm/MHz.

(15)

2.2 Monostatic and Bistatic radar

Radar systems available can be of different types. When the receiver and transmitter are separated the radar is said to be operating in bistatic mode. If a radar unit is transmitting and receiving, it is said to operate in monostatic mode. It is possible to use one transmitter together with more than one receiver. The system is then said to be operating in multistatic mode.

The advantages of using a multi-static setup are twofold: More reflections from targets could potentially reach other transceivers and the use of separated radars makes the system less vulnerable in military applications. The caveat of using a bistatic setup is that the synchronization is more complex and the transmitted power needed would increase if the transmitter and receivers are separated far away [11]. There are two radars available for use in this thesis project and both were operated in monostatic mode instead of bi-static, as the latter choice would result in measurements from only one device in this specific case, which would make any kind of accurate tracking more or less impossible.

2.3 PulsON 410

The UWB radars used in this thesis are the PulsON 410 (Figure. 2.2), manufactured by the U.S company Time domain [12]. As is illustrated in the figure, it consits of two antenna ports that allows it to transmit and receive over either the same or different antennas. The power is supplied at 12 V via the power jack marked in Figure 2.2b. The radar module consists of an extremely accurate RF front-end that operates over the frequency band 3.1 GHz - 5.3 GHz with a center frequency of 4.3 GHz, which according to Equation (2.3) gives a range resolution ∆r ≈ 6.7 cm. It should be emphazied, however, that as detection needs to be performed on the reflected signals from the radar, the resolution is significantly less in practice. The device can operate in dedicated range mode (not as a radar), which, however, requires multiple devices as ranging is done directly between the units. It is in this mode that resolution becomes extremely good. The transmit power of the unit can be specified within a certain range according to application needs.

The PulsON 410 is connects to the host via USB - micro USB. The host is a computer which communicates with the radar module over API:s specified by manufacturer. Depending on what mode the user needs to operate the device in (Radar, Ranging and Communications, Channel analysis) separate API:s are used for issuing requests. The only API that has been used in this thesis is the Monostatic Radar Module API (MRM), which allows the device to operate as a monostatic radar. Each device comes with two separate antenna connectors and pair of broadspec omnidirectional antennas ship with each device. If necessary it is possible to change the antennas to other types.

(16)

(a) PulsON 410 viewed from side. The antennas are the two large green components on the device, connected to the SMA connectors.

(b) PulsON 410 viewed from above.

Figure 2.2: The PulsON 410 Radar, Ranging And Communication module.

3 Data acquisition and post-processing

3.1 Acquiring data

The radar modules function by accepting commands from a host (PC or laptop) and returns responses (or data, if issued to scan) to these commands. Each scan issued by the host to the radar is divided into “bins”, where each bin is 1.907 picoseconds (ps). The host specifies a start and a stop time of the scan and the embedded software of the radar converts this

(17)

into the correct number of bins. The radar scans are gathered by the device using multiple samplers in paralell acting as a rake receiver. The scan data is quantized into partitions of 96 readings that takes 5859.36 ps (also called “quanta”, the rake sampler size). This means that the start and stop times that the host specifies will be rounded if they are not given in even multiples of quanta. What determines the duration of a radar scan is besides the quanta is the pulse integration index (PII). As the radars are transmitting coherent pulses (constant in phase difference), it is possible to use pulse integration for improving the Signal-to-noise ratio (SNR) of the return signal. Based on the assumption that the noise of the system is independent and identically distributed (IID), the idea of pulse integration is to accumulate samples for each transmitted pulse with the assumption that the signal returns are significantly larger than the noise. This results in the accumulated return pulses having a higher SNR than without using pulse integration.

The final scan time (given in µs) is a function of the quanta and the PII, given by Equation (3.1) [12]. The scan time defines the maximum update rate that can be achieved for a given configuration.

Tscan= #quanta × 0.792µs × 2P II (3.1)

When the radar has sampled the reflected pulses they are converted to relative ampli-tudes. The output to the user are these relative amplitudes over the specified scan duration. The reflected pulses will be distorted and noisy versions of the transmitted ones, with reflec-tions not only originating from the direct target but also delayed reflecreflec-tions of the targets (multipaths). Antenna coupling and static objects that cause reflections of transmitted sig-nal in the surveilled area (commonly referred to as clutter) could also make the detection a challenging task and needs to be accounted for.

A signal model that has been used in this thesis to describe a received signal (or mea-surement) at time t is given by Equation (3.2)

ξ(t) = w(t) + n(t) (3.2) where ξ(t) is the received signal, w(t) is the returned signals from the objects in the surveilled area and can be expressed as w(t) = h(t)∗s(t) where ∗ denotes convolution, h(t) is the channel impulse response and n(t) is noise. The expression can similarly be written in the state space form shown in Equation (3.3), with ξ being the received signal, C(r) representing the measurement equation where r is the distance to the target, and n representing measurement noise. When considering movement in two dimensions, the measurement equation can be considered to consist of contributions from an x and y component of the target location. Details regarding the state space model and how it is extended in 2D will be explained in Section 5

ξ = C(r) + n (3.3) As radar scans are digitized by analog-to-digital converter (ADC), the signal model can readily be extended and expressed as a digital signal with m denoting the scan number (i.e M scans are performed where each scan consists of K samples) given by

ξk,m= wk,m+ nk,m (3.4)

where k = {0, 1, . . . , K − 1, K} and m = {0, 1, . . . , M − 1, M }.

The relatively simple signal model that has been presented serves as basis when presenting the detection and tracking methods in subsequent sections. In practice, multiple signal paths (multipaths) caused by multiple reflections can influence the received signal significantly. This is shown in a typical collected radar scan illustrated in Figure 3.1. In the collected scan, there are multiple potential targets and “shadowing” that appears to be multipath delays.

(18)

The actual number of targets for this specific scan is one. Without further processing it would be difficult to state something valuable about the number of targets.

0 200 400 600 800 1000 1200 Sample -4 -3 -2 -1 0 1 2 3 4 Relative Amplitude ×104 Potential target Probable multipaths Antenna Coupling

Figure 3.1: A typical radar return signal which may contain direct reflections and multiple reflections (multipaths) from targets.

3.2 Visualizing data

To get an idea of what the aquired data look like, so called Radargrams are commonly used. These plots show how the amplitude of the samples aquired in each individual scan change over time. The dimension showing the samples of the scan is usually called Fast time (order of nano-seconds) and the dimension showing successive scans is called Slow time (in the order of seconds). The radargrams serve as a good means to get an intuitive feel of where potential targets might be located in the surveiled scene and can be useful to get a picture of potential clutter present in the scene. A typical radargram for a set of 1200 successive scans with a single target moving in circles in a controlled area is shown in Figure 3.2. The x-axis is the scan samples (which here have been converted to distance by the two-way time of flight) and the y-axis is successive scans. The periodic intensities of the returned amplitudes is apparent because of the circular motion the target was doing.

(19)

Figure 3.2: Radargram from one sensor of a single target.

3.3 Background subtraction

To perform automatic detection and tracking, the raw data collected by the radars first needs to be processed to enhance targets that might be present in the surveilled area. This is done in two steps: Background subtraction and moving target enhancement. As the transmission power of UWB is very weak and noise-like, methods of signal processing needs to be applied to increase the SNR. There exist various background subtraction techniques for processing raw radar data [13]. Three techniques that have been tested in this project are

• Averaging and subtracting static background (similar to basic image foreground detec-tion).

• Singular Value Decomposition (SVD) • PCA (Principal Component Analysis)

Basic background subtraction by averaging radar scans of the empty scene and subtracting from the non-empty scene, together with moving target enhancement and the pulse integra-tion performed by the radar proved to provide sufficient results for performing detecintegra-tions. As of this, the SVD and PCA were not further considered. It should be mentioned that using long pulse integrations is not very favorable if the targets are moving fast, as subsequent pulse returns could be very different and the result would not benefit much from the integration.

(20)

Because of this, the integration factor of the radars were set to a value of P II = 8 (mini-mum being 6, maxi(mini-mum 15), equalling 28 _{= 256 pulses being per sample for each slow-time}

scan. Each time the PII is doubled, the SNR improves with 3 dB, with the catch that each slow-time scan takes longer to finish.

3.4 Moving target enhancement

After background subtraction had been performed, the SNR of the signal returns still proved to not be good enough for performing target detection. As the aim is to perform tracking of moving targets, a simple procedure of subtracting successive scans to enhance changes over time was applied. For all samples ξk,m, the resulting signal ˜ξk,m becomes

˜

ξk,m= ξk,m− ξk−1,m (3.5)

4 Target detection

The detection of targets in a scene basically becomes the task of finding prominent peaks in the reflected signals while at the same time discarding noise and the peaks that arise from non-targets. This is a non-trivial task and there exist several techniques for approaching this problem [14, 15, 16]. Two techniques have been considered in this thesis, namely

• CLEAN Algorithm

• Constant False Alarm Rate (CFAR)

In the following sections each detection method will be described in more detail.

4.1 CLEAN algorithm

The CLEAN algorithm was first presented in [17] as a way of performing deconvolution on images aquired in radio astronomy.

The idea of using CLEAN is to estimate the channel impulse response h(t). As the transmitted signal is known (or can be measured), knowing h(t) would imply that it is possible to extract a noisy target signal. Given the signal model presented in Equation (3.2), the received signal can be described as the convolution of the channel impulse response with the transmitted signal

ξ(t) = h(t) ∗ s(t) + n(t) (4.1) where ∗ denotes the convolution operator. Finding h(t) is the same as deconvolving the channel impulse response (CIR) with the transmitted signal. This is generally an ill-posed problem that is often unstable. Assuming white noise and recalling that taking the Fourier transform of Equation (4.1) yields ξ(f ) = H(f )S(f ), the Fourier transform of the channel im-pulse response can be computed as H(f ) = _{S(f )}ξ(f ) which is unstable as S(f ) → 0. The CLEAN algorithm applies an iterative approach to avoid this ill-posedness to yield an approximation of the CIR [18].

In the standard CLEAN deconvolution procedure, the sampled channel impulse response hk is found by first computing the cross-correlation in Equation (4.2) of the transmitted and

received signal in a loop. A “dirty” map is initialized as equivalent to the received signal. This map is used to subtract samples from as highly correlated samples are found. The iterations continue until the cross-correlation goes below a pre-set treshold. During each iteration, the maximum value of the cross-correlation of a complete scan (that is, all K samples) is computed and the samples is added to a vector containing the estimated channel impulse response. Whenever a peak is found to exceed the treshold, it is removed from the signal and N neighbouring samples of the peak are padded with zeros. This process continues until

(21)

there are no more peaks found that exceed the treshold. The result of the algorithm is a vector with samples at values where ξ_k and s_k are most correlated.

Rξs = ξk∗ sk (4.2)

A version of the algorithm has been suggested and used for estimating UWB indoor chan-nel impulse responses in mild multipath environments [19]. As with the standard CLEAN pro-cedure, the algorithm used in [19] assumes the transmitted signal characteristics are known. As this could be measured in many scenarios, a modified version of the CLEAN algorithm [15] that does not require any knowledge of the transmitted signal was however instead chosen and implemented. Whenever the CLEAN algorithm is mentioned from here on, it will imply this modified version of the algorithm. With some loss of performance, the modified version computes the power of the the received signal instead of computing the cross-correlation between the transmitted and received signals.

The estimated time-of-arrival is the corresponding time at which each peaks was detected in the CLEAN iterations. A large peak indicates high likelihood of a target being present. Knowing the scan duration, the time corresponding to the found sample indices can readily be converted to range measurements. The performance of the detection results using this algorithm depends heavily on a “good” choice of the treshold, together with a suitable choice of N . To avoid targets being filtered out in multi-target cases, N needs to be chosen with specific care in these scenarios. The specific implementation of the CLEAN algorithm steps is given in Algorithm 1.

Procedure 2 The modified version of the CLEAN algorithm. The modifications makes it possoble to apply the method without any knowledge about the transmitted signal.

Input: r,N,T_hd . received signal, padding size, Peak treshold Output: TOA . Time of arrival

1: _{procedure clean}

2: T OA ← [] . Initialize empty list

3: p ← r2 . Power of received signal

4: while 1 do

5: τn← argmaxn(p) . Index of largest peak

6: αn← p(τn) . Amplitude of largest power peak

7: if αn> Thd then

8: T OA ← [T OA τn] . store current peak index

9: p(τn− N : τn+ N ) ← 0 . Pad neighbouring samples

10: else

11: break

12: return T OA

4.2 Constant False Alarm Rate

Constant false alarm rate (CFAR) detection is a common form of adaptive algorithm used in radar technology to detect the returns from targets against a background of noise, clutter and interference [16]. The CFAR algorithm determines the threshold of power above which any return can be considered to probably come from a target. It can be implemented in different methods, e.g., cell-averaging CFAR (CA-CFAR), generalized censored mean level detector CFAR (GCMLD-CFAR), greatest-of CFAR and least-of CFAR. In this thesis, CA-CFAR and GCMLD-CA-CFAR algorithms are going to be addressed in detail. CA-CFAR methods takes a sliding window-approach on estimating peaks in the received signal. This is done under condition that the false-alarm rate P_{F A} cannot exceed some user-defined value. This

(22)

property of CFAR methods makes them useful in alarm-critical applications, such as missile defence radar systems.

The CFAR method expect a square-law input signal, which means that the power of the collected radar scans must be computed. Some important design parameters of CFAR are the center window cell under test (CUT), guard windows (G) surrounding the CUT and leading and lagging reference windows (Nlead,lag) surrounding the guard cells. The concept of samples

that represent reference cells, guard cells and CUT is shown in Figure 4.1. The number of samples that these windows consist of does not need to be the same and “good” choices of window lengths is commonly determined empirically. As CFAR methods per design attempts to keep the false-alarm rate constant, this also becomes an important design parameter when applying the method. As the noise plus interference in the reference cells vary, so does the treshold that the CUT power is compared with. This results in the constant false-alarm rate being ensured.

Reference cells Guard cells CUT Guard cells Reference cells

Figure 4.1: The different cells of the sliding window in the CFAR method. The reference cells consists of N_lead and N_lag number of samples, while the guard cells consists of N_G samples each.

The method estimates the noise power in each of the reference windows N_lead,lag and compares the CUT power Z to a treshold T . The general expression of the treshold is given by Equation (4.3)

T = ˆg(Nlead, Nlag)α (4.3)

where ˆg(Nlead, Nlag) is the estimated noise power and α is the CFAR constant. This constant

is analytically derived and depends on how the CFAR is implemented [20]. The detection procedure can be summarized as binary hypothesis test

Z ≷1

0

T (4.4)

where 1 corresponds to a detection and 0 corresponds to no detection.

A typical implementation of the CFAR algorithm is illustrated in Figure 4.2. The samples of the square-law signal are denoted by X, where X_k=1:K,m = ˜ξ_k=1:K,m2 . The reader is reminded that the reference cells are denoted N_lead,lag and the guard cells given by G. The complete CFAR window consists of the reference cells, guard cells surrounding the CUT together with the CUT itself, as shown in Figure 4.2. This window is sliding over all samples of the squared signal. The algorithm is applied individually for each scan and the time-of-flight at sample indices exceeding the treshold is computed and converted to range measurements. The obtained ranges are considered target detections. To get an idea of what it looks like when applying CFAR to real radar data, a CFAR iteration on a radar signal is shown in Figure 4.3. In this figure it is clearly seen that the signal exceeds the treshold and detections will be registered around samples ≈ 490 − 520. The general CFAR algorithm is for convenience presented in Algorithm 3.

(23)

Input signal (X_k= ˜ξ_k2) X1 XN₂ G CUT G X2+1N XN

Reference samples Reference Samples

Noise power estimate

.

CFAR constant α Comparator Detection / No detection ˆ

σ2_i

αˆσ_i2

Figure 4.2: Schematic of CFAR method.

Procedure 4 Detailed steps in pseudo-code for the general CFAR algorithm (i.e no specific CFAR constant computation is shown in the algorithm).

Input: x,i,N,M,G . Scan, Sample index, Reference cell length, sample cell length, guard cell length

Output: H . Binary hypothesis that equals 1 if detection, 0 if no detection.

1: _{procedure CFAR} . Ensure M is even or 1 and N is even.

2: r ← x2 . Needs square-law input to CFAR algorithm

3: if M == 1 then

4: CU T ← i

5: else

6: CU T ← r(i − M/2 : i + M/2) . If at edge, form as many samples as possible.

7: Wlag ← r(i − M/2 − G − N/2 : i − M/2 − G) . Form lagging window of size N/2. If

at edge, form as many samples as possible.

8: Wlead ← r(i + M/2 + G : i + M/2 + G + N/2) . Form leading window of size N/2. If

at edge, form as many samples as possible.

9: . Depending on the CFAR method used, the noise power estimation and CFAR constant is different. For sake of brevity, CA-CFAR will be shown.

10:

11: α ← N (P

−1 N

F A − 1) . CA-CFAR constant

12: σ_i2← mean(W_lag+ Wlead) . Noise power estimate

13: Thd← σ2iα . Treshold to compare CUT power with

14: if mean(CU T ) ≥ T_hd then

15: H ← 1

16: else

17: H ← 0

(24)

400 500 600 700 800 900 1000 1100 Fast-time sample -4 -2 0 2 4 6 8 10 12 14 16 Power ×105 Return signal Adaptive Treshold

Figure 4.3: CA-CFAR illustrated on a real radar return signal for N_lead= 80, G = 30, Nlag=

30, PF A = 1 × 10−5 (one single slow-time scan). The treshold (in red) is shown overlaid on

the return signal (in blue). 4.2.1 Cell-Averaging CFAR

CA-CFAR is one of the most widely used CFAR algorithms, and the algorithm is implemented in the MATLAB standard library. In CA-CFAR the noise power is estimated as the average power in the reference cells. Assumptions are made on the noise statistics when deriving α for CA-CFAR. This degrades the practical performance of the method slightly, as all of the assumptions are rarely fulfilled. To derive an expression for the CA-CFAR constant, it is assumed the received signal power is exponentially distributed according to Equation (4.5). Note that X here represents the square-law output (i.e. the power).

pX(X) = 1 σ2 i e −X σ2_i _, _{X ≥ 0} _(4.5)

where σ_i2 is the average noise power that is estimated from the reference cells. By integrating the probability distribution given in Equation (4.5) one can get an expression for the proba-bility of detection P_D. From this, a probability of false alarm P_{F A} can be obtained. Details regarding these calculations can be found in Appendix A.1. The resulting treshold is [20]

TCA = αCAσˆi2 (4.6)

in which the CA-CFAR constant is given by αCA = N [P

−1 N

F A − 1] (4.7)

where N is the complete length of the reference windows N_lead,lag and σi is the estimated

interference power in the reference cells given by σ2_i = _N1 PN

n=1rn, where R = {r1, r2, . . . , rN}

(25)

CA-CFAR performs well in homogenous environments and single-target situations. If the surroundings are considered heterogenous or if there are multiple targets closely spaced, returns that appear in the reference windows of the CFAR method will increase the treshold, potentially making a target detection in the CUT fail. There is also a chance of sudden noise-floor level changes which could result in an increased number of false alarms or targets being missed in the transition regions of the noise-floor levels [20].

4.2.2 Generalized Censored Mean Level Detector CFAR

To counter the problems that the CA-CFAR experiences in multi-target scenarios and chang-ing noise-floor levels, several different methods for CFAR have been proposed, e.g., called Robust CFAR methods [20, 21, 22]. One of these that do not require any prior knowledge on the number of targets in the scene is the GCMLD-CFAR [22]. The idea of GCMLD-CFAR is to “censor” out strong interferers in the reference cells to get a more accurate measure of the actual background noise level. The “level” of censoring that is performed depends on a selected parameter called probability of false censoring, PF C. Typically this value is chosen

to be the same as P_{F A}[22]. The procedure of censoring in GCMLD is summarized as follows: The samples in the reference cell are ranked according to magnitude to yield ordered samples q1 ≤ q2 ≤ · · · ≤ qK where qi are reference cell samples and K is the number of

samples in the cell being censored. This ranking is done seperately for the leading and lagging reference cells. The method then proceeds by assuming that the lowest ordered sampled q₁ in the window under consideration is an estimate of the background noise. The next sample q2 is compared to q1T1, where T1 is a scaling constant that depends on PF C. If q2 ≥ T1q1,

the sample q₂ is declared an interfering target and will be censored out (removed from the reference cell). If q2 ≤ T1q1, the sample q2 is decided to correspond to only noise and the

censoring continues by forming the sum of the samples considered so far s₂= q1+ q2. Next,

q3 is compared to T2s2. The procedure then continues in the same manner as for previous

samples until the first sample determined to correspond to an interfering targets is found, or until all samples have been checked. The number of interfering targets are kept track of in each window being censored as m₁ and m₂ respectively, with the total number of estimated interfering targets being m = m1+ m2.

As CFAR methods normally operates on single-pulse returns, it was slightly modified as pulse integration with the radars used in this thesis was required. When performing the cen-soring only the current reference cell sample that exceeded the treshold was removed, instead of all remaining samples. This proved to provide significant performance gains compared to applying the regular GCMLD censoring scheme.

The censoring procedure that is applied to each reference window individually is presented in pseudo-code in Table 5.

(26)

Procedure 5 Modified censoring procedure for GCMLD-CFAR

Input: q . Samples of reference window sorted by magnitude. Output: ˜q, mi . Censored samples, estimated interferring targets.

1: _{procedure censoring} 2: M ← length(q)

3: s ← q1 . First sample considered

4: m ← 0 . Initialize interferring targets as 0

5: k ← 2 . Index of next sample being considered

6: Tc← αF C,M −k . Compute all tresholds

7: while true do

8: T ← Tc,k . Get current treshold

9: if qk ≥ T s then

10: qk← 0 . Censor only current sample instead of all remaining

11: mi← mi+ 1 . Increment target count . In the non-modified version

the censoring would return here. Instead it now continues checking if remaining samples exceed the treshold.

12: else 13: s =Pmin(k+1,M ) i=k qi 14: k ← k + 1 15: if k ≥ M then break 16: return ˜q, m

As was the case with the CA-CFAR method, the P_{F A} is derived from P_D, which is obtained by integrating the probability distribution of the assumed signal power. Some details regarding derivation of the treshold can be found in Appendix A.1. The interested reader should refer to [22] for a thorough derivation. The resulting constant and corresponding treshold is given by Equation (4.8).

TGCMLD = ˆσ_i,censored2 (4.8)

in which the GCMLD-CFAR constant is given by

αGCMLD= (N − m)P

1 m−N

F A − 1 (4.9)

where N is the length of the reference windows and ˆσ_i,censored2 is the estimated noise power (averaged) in the reference cell after censoring has been performed.

4.2.3 CFAR performance

Depending on the CFAR method used the performance can vary a lot for different Signal-To-Interference Noise Ratios (SINR). A useful tool for visualizing exactly how much the P_D differs with regard to SINR and parameter choices for different CFAR methods are Receiver Operating Characteristics (ROC). The ROC curves show the probability of detection given SINR, N, M, G and P_{F A}. The ROC curves for some different of parameter selections are shown Figure 4.4a - 4.5b to illustrate the degradation in performance of CA-CFAR when multiple interferring targets are present.

(27)

0 5 10 15 20 25 30 35 SINR [dB] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Probability of detection GCMLD-CFAR CA-CFAR Neyman-Parson Detector

(a) Detection curves for PF A = 1 × 10−4, N = 80 and m = 1. As the

number of interferring targets m = 1, the performance of the CA-CFAR (blue) and GCMLD-CFAR (red) are the same (curves are overlaid in plot). 0 5 10 15 20 25 30 35 SINR [dB] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Probability of detection GCMLD-CFAR CA-CFAR Neyman-Parson Detector

(b) Detection curves for PF A= 1 × 10−4, N = 80 and m = 2.

Figure 4.4: ROC curves for CA-CFAR (blue), GCMLD-CFAR (red) and the optimal Neyman-Pearson detector (yellow) (for reference).

(28)

(a) Detection curves for PF A= 1 × 10−4, N = 80 and m = 3.

(b) Detection curves for PF A= 1 × 10−4, N = 80 and m = 4.

Figure 4.5: ROC curves for CA-CFAR (blue), GCMLD-CFAR (red) and the optimal Neyman-Pearson (yellow) detector (for reference).

5 Target tracking

As the goal of the thesis was to track multiple targets simultaneously a systematic approach to solve the multiple-target and data association problem called Multiple Hypothesis Tracking (MHT) was chosen. MHT is well adopted in the tracking community and has been around since the 1980’s [8]. The tracking is performed in a centralized manner, meaning each radar sends its measurements to a central processing unit where the detection, data fusion and tracking are performed. In the next sections the core concepts of a version of the MHT framework will be explained.

(29)

5.1 Measurement and Motion Model

A radar transmitter emits pulses. Then the pulses propagate and get reflected by objects (including targets) in the surveilled area. The reflections (returns) are measured by one or more radar receivers. The 1D measurements obtained from the radar can be represented by the measurement equation given in Equation (3.3), repeated here for convenience.

ξ = C(r) + n

where ξ is the measurement, n is the measurement noise, r is the distance to a target, and C(r) denotes the measurement matrix, given by (in discrete time)

C(rk) =

q

x2_k+ y2_k (5.1) where k denotes the sample index. The measurement equation C(r_k) is non-linear and the states of the targets that are to be tracked is hidden (it is in 2D, observations in 1D) which makes standard techniques such as linear Kalman-filtering non applicable. To obtain a measurement model that can be used in practice, the non-linear model C(r_k) can be linearized around (xk, yk) by means of the first order Taylor expansion under the assumption

that the targets are moving with linear motion (i.e. constant velocity model). This results in a linearized measurement equation (the Jacobian of C(r_k)) given by

H(rk) = xk √ x2 k+y2k yk √ x2 k+yk2 T (5.2) As an accurate model of the target motion is unkown, a constant-velocity Brownian motion model is assumed. This model can, in disrete time, be described in 1D by Equation (5.3)

rk+1= rk+ σvk (5.3)

where r represents the range and σ_v is the process noise. Both the measurement noise variance σ_r2 and process noise variance σ_v2 had to be empirically chosen as sufficient data related to target velocity and detection errors was not available. The measurement variance was selected based on how accurate the detections appeared, while the process noise variance was tested for many different values in the range 0.01 ∼ 1.0 and fine-tuned to yield good performance. Details regarding to the state prediction is further explained in Section 5.2.4.

5.2 Track-Oriented Multiple Hypothesis Tracking

Track-oriented multiple hypothesis tracking (TO-MHT) is a favored algorithm for tracking in difficult environments with high computational resources. In MHT, targets are said to move along “tracks”. A track is a sequence of hypotheses formed after receiving multiple observations. The core concept of MHT is that, as observations are received at a certain time step, hypotheses are formed and kept until later time instants [8]. The idea is that additional received observations eventually should resolve any data association problems that might be appearent at first. The TO-MHT has a tree structure where hypotheses make up the nodes of the tree. As new observations are received and associated, the tree adds or extends branches by adding or extending hypotheses. Tracks are “paths” of hypotheses in the tree from the root hypothesis (first associated observation) to the most recent ones. As the beliefs are propagated at each time instant there is a potential exponential growth in the number of hypotheses. Clearly this is not a very tractable property of the MHT and there are multiple ways to deal with this issue that will be further discussed in Section 5.2.3. An overview of the MHT logic can be seen in Figure 5.1 and the tree structure can be viewed in Figure 5.2.

(30)

Note that the term “track hypotheses” will be used, meaning a sequence of hypotheses that constitute a track. Observations Gating Track Creation Score & N-scan pruning Track prediction

Figure 5.1: General logic of MHT.

1 2 1 2 4 3 2 0 scan k − 2 scan k − 1 scan k

(a) The tree-structure of the MHT. The numbers inside the nodes correspond to different ob-servations where 0 denotes the “dummy” observation (no mea-surements associated). 1 2 1 2 4 3 2 0 scan k − 2 scan k − 1 scan k

(b) One of the tracks highlighted in orange in the example hypoth-esis tree.

Figure 5.2: MHT tree structure. 5.2.1 Gating

Gating is the procedure of associating measurements with tracks. For existing tracks the association is based on statistics computed from the last prediction of the state. The mea-surement is considered to be inside the gate if the so called squared Mahalanobis distance is less than a chosen treshold [23],

d2 = (ˆxl_k− y_i_k)T(Σ_kl)−1(ˆxl_k− y_i_k) ≤ Gth (5.4)

where ˆxl

k is the state prediction, yik is the i

th _{observation, G}

th is the gate size and Σlk is the

innovation covariance of the lth track at time k. The innovation covariance is computed as Σl_k= HP_k−1l HT + R (5.5) where R is the measurement covariance matrix and P_k−1 is the state covariance of the track estimate up to the last received measurement. The innovation gate size Gth is χ2 distributed

with n degrees of freedom, where n is the dimensionality of the states. The gate size can thus be chosen according to desired probability of observations falling within the gate. In the case of gating in 2D, Equation (5.4) corresponds to an ellipse centered in the tracks predicted positions.

When a measurements falls within the validation gate of a track, the track initializes a hypothesis related to the measurement. The same measurements could potentially be gated

(31)

to several tracks, generating multiple hypotheses. This ambiguity is however resolved at a later stage in the MHT framework. If a measurement does not fall within the gate of any track or if there are no existing tracks, a new track tree is initialized with a hypothesis associated with the corresponding measurement. As a convention, a dummy observation is added to tracks that are not associated with any measurements at a given time. The gating procedure is illustrated in Figure 5.3. Note that if multiple measurements are received at a given time step, this means Mahalanobis distance will be computed for each observation-prediction pair and compared to the gating size.

1 2

Figure 5.3: Gating procedure of associating observations with existing tracks. The track is shown with the line and cross, while the observations are shown as the circles. Here observation “2” is associated with the track and the other observation that did not fall within the gate would initiate a new track.

5.2.2 Track Scoring

Each node in the tree of hypotheses have an associated score. This score is based on the Likelihood Ratio (LR) which is a probabilistic measure of how likely a track is to be valid. The general definition of a LR is the ratio of the following probabilities [24]

LR = p(D|H1)P0(H1) p(D|H0)P0(H0)

(5.6) where H1,0 is the hypotheses of a true target and false alarm respectively, D is the received

data, p(D|H_i) is the pdf (probability density function) of the received data assuming Hi

is correct and P0(Hi) is the a priori probability of Hi. Per definition this is equal to the

probability of detecting a true target P_T over the probability of detecting a false alarm P_F. In practice it is convenient to work with the logarithm of likelihood ratio (LLR), which based on eq. (5.6) is defined as

LLR , ln(PT PF

) (5.7)

Let i denote the ith observation where i = 0 corresponds to a missed observation. Using an approach of track scoring presented in [23], the LLR score S_kl at time k can be computed recursively as

S_kl = S_k−1l + ∆S_kl (5.8) where the difference in score is given by

∆S_kl = ( ≈ ln(1 − P_D), i = 0 ln(Vs 2π) − 1 2ln(|Σ l k|) − d 2 2, i 6= 0 (5.9) where V_s is the measurement space which in the case of 2D tracking is the surveilled area. Upon creating a new track, the score was initialized with the same score as for a missed observation. As new hypotheses are added to previous ones, the scores are accumulated. This implies that at each bottom node of the hypotheses tree, the sum of the LLR score of

(32)

the corresponding branches (tracks) will be available. The LLR can also readily be converted to probability of a track being valid using [24]

Ptrack=

eLLR

1 + eLLR (5.10)

The probability P_track is used when evaluating if a track should be deleted or not. A typical choice of tresholds used for track deletion is ≈ 1 × 10−3 [24].

5.2.3 Track Pruning

As previously mentioned the number of generated hypotheses grows expontentially. To avoid this unmanagable amount of hypotheses, deletion of tracks has to be performed in some systematic manner. This is done by Pruning (deleting) unlikely tracks. Pruning can be done with regard to different properties [24], two that have been considered in this thesis are Score pruning and Global pruning.

Score pruning At each iteration, the tracks are scored according to eq. (5.8) and the probability of a track being valid is calculated using eq. (5.10). Tracks that have Ptrack

less than a certain treshold will be deleted.

Global pruning Global pruning is performed to find a set of tracks with maximum LLR scores that are compatible. For a track to be compatible with another, they cannot share observations with any other track at a given time. Tracks that diverge from the global tracks (set of “best” sequence of hypotheses) N steps up the tree are pruned from the hypotheses trees. This is also commonly called “N-scan pruning” [24], where N is an integer value with typical choices being between of 3-5.

The procedure of finding a compatible set of tracks with maximum score can be formulated as a Maximum Weighted Independent Set (MWIS) problem, which in general is a NP-hard optimization problem [23]. The tree structure of the track hypotheses can be represented as a undirected graph G. Each track Tlin the graph is represented by a vertex vl∈ V where V

is the set of all vertices (tracks). If tracks T_l and Tj share observations, they are connected

by edges (l, j) ∈ E in the graph. Each vertex also has a weight w_l, which is the LLR score of the track. An independent set of maximum weight is a set of track hypotheses that together yield the largest sum of weights, while having no shared observations at a given time. This optimization problem can thus be formulated as

max

v

X

l

wlvl, vl+ vj ≤ 1, ∀(l, j) ∈ E, vl∈ {0, 1} (5.11)

where the inequality of the vertices simply states that a given track (and hence an observation) can only be used once. An example graph of tracks with some given weights and sequence of observations can be seen in Figure 5.4. The MWIS is marked in orange.

(33)

T1(6) 5 2 4 T2(12) 3 5 4 T3(3) 5 1 2 T4(21) 5 3 1 T5(3) 1 2 3

Figure 5.4: Track graph constructed from tree representation. The MWIS is shown in orange, where the associated observations at successive time instants are shown above each vertex. The vertex weights are shown in parentheses.

There are multiple ways to solve the MWIS optimization problem with varying algorithm performance. In [24] the author suggests finding the N best tracks by use of Lagrangian Relaxation or using a modified version of the Viterbi algorithm. It was however instead of using the aforementioned alternatives chosen to solve the MWIS using two different ap-proaches presented in the field of graph theory: One approximative method using a Greedy Algorithm called GWMIN and one method solving for the exact solution using the Bron-Kerbosch algorithm. Greedy algorithms are sub-optimal in the sense that they will find a solution satisfying the independent set criteria, however the sum of the track scores may not be the global maximum. It was shown in [25] that a lower bound on the sum of the weight of the independent set can be guaranteed, given by

Wtot ≥

X

v∈V (G)

W (v)

d(v) + 1 (5.12) where W_tot is the total weight of the independent set, V (G) is the set of all vertices in the graph G, W (v) denotes the weight of a specific vertex and d(v) denotes the degree of a vertex (number of connecting edges).

The Bron Kerbosh algorithm is an algorithm presented in 1973 for finding all cliques in undirected graphs [26]. A clique is a set of vertices that share edges in a graph. The problem of finding independent sets of vertices is the complement to the problem of finding cliques, which makes the extension to the Bron Kerbosh algorithm intuitive (simply solving the original clique problem on the complement of the original graph provided to the algorithm). An extension to the original algorithm for handling graphs with weighted vertices was presented in [27], which is the equivalent of the MWIS problem that arises in TO-MHT.

5.2.4 Track Prediction

The part of the MHT framework that states where the targets are assumed to be located at a given time is the track predictions. The predictions are necessary to propagate a target position in time and is also used during the next iteration of the MHT to compute the valida-tion gate for associating new measurements. To perform tracking certain filtering techniques needs to be applied. In this thesis the targets are being tracked in 2D Cartesian Coordinates. The state of a targets being tracked is represented by

(34)

A technique that is commonly used in MHT implementations is Extended Kalman Filtering (EKF) [28]. In this thesis, however, it was decided to perform the state prediction using particle filters.

The use of using particle filtering instead of EKF was motivated by the large uncertainty in the assumption of a motion model. Whilst an EKF would perform badly given a non-accurate motion model, a particle filter would not suffer from the same issues as the same strict assumption of motion is not implied. As discussed in Section 5.1, the represented states are two-dimensional with the measurements being 1D and measurement equation non-linear. This fact makes it possible to represent the states by a Hidden Markov Model, which can be solved by use of particle filtering.

The concept of particle filtering will be explained in Section 5.4. Each of the initiated tracks in the MHT will predict its state with a separate particle filter. The gating procedure (computing Mahalanobis distance) is also affected on the accuracy of the state predictions, with poor predictions resulting in an increased gate size.

5.2.5 Track Merging and splitting

The situation under which two targets combine together to form a larger target is called merging, while the situation under which targets are allowed to break into smaller target groups is called splitting. Normally in multiple-target tracking systems the events of tracking same targets need to be handled. This is typically done by merging tracks that share similar state estimates or associated observations over an extended period of time [24]. As the radars only measure distance, any splitting of tracks would be very difficult as no specific features of the targets would be available for resolving merge-split cases later in time. Due to this issue, merging and splitting were not considered in greater detail in this thesis.

5.2.6 Track Clustering

A common way of making a MHT implementation more efficient is to form clusters of hy-potheses. The most straightforward way of forming clusters is to group all tracks that are incompatible (share observations at a given time or originate from the same root observa-tion). This feature naturally arises when creating and extending the hypotheses trees as all branches in the same tree is incompatible per definition, as they share the same root node. Further formation of clusters can account for tracks of different trees that share observations. With clusters formed, the gating and scoring can be done independently from other clusters thus making it suitable for parallell processing implementations.

5.3 Hidden Markov Models

Before introducing particle filters it is useful to know why they are needed. To breifly motivate this, a means of modeling stochastic systems called Markov Models should to be mentioned. Markov Models is a way of modeling changes in states of a stochastic system. The model assumes that a system obeys the Markov Property, stating that future states are conditionally independent of all previous states except the current. Let a sequence of discrete states be denoted by {x₁, x2, . . . }. Let X = {X1, X2, . . . } be the possible sequence of states the

random system has taken on up to the current time. The Markov Property can then be expressed as

P (Xk+1= xk+1|X0 = x0, X1 = x1, . . . , Xk−1= xk−1, Xk = xk)

= P (Xk+1= xk+1|Xk= xk)

(5.14) The sequence of states Xi:k the random system takes on is called a Markov Chain. For a

(35)

Markov Model, however, does not provide this privilege, instead only observations related to the current state of the system can be observed (thus the term hidden). The difference between a regular Markov Model and Hidden Markov Model is illustrated in Figure (5.5).

X1 X2 X3

(a) Stationary Markov Model.

X1 X2 X3

z1 z2 z3

(b) Hidden Markov Model.

Figure 5.5: Stationary Markov Model and Hidden Markov Model. In the SMM the states are directly observable whereas for a HMM, only observations related to the current state of the system can be obtained.

In the context of this thesis, the 2D positions of the targets being tracked can be considered the hidden states. The observations obtained containing information about these hidden states are the distances measured by the radars. If the observations available in a HMM up to time k are denoted by z_1:k, then finding the target states is equivalent to solving the filtering problem of the HMM. Ultimately the goal is to find the posterior pdf p(xk|z1:k).

By using basic inference rules, likelihood weighting, normalization and a resampling step, particle filtering provides a systematic approach to solving the problem. An introduction to particle filtering can be found in [29], and a more thorough explanation is presented in [30]. The basics will be presented in the following section.

5.4 Particle Filtering

Particle filtering is a type of Monte-Carlo method for approximating the pdf of the state of a system described by a HMM model. The generality of HMMs is not restricted to linear systems so that approximations of solutions to the nonlinear filtering problems can likewise be obtained.

The state (x, y, ˙x, ˙y) is represented by a set of N samples {xi_k−1}N

i=1 of different weights

called particles, where i denotes the ith particle. The type of particle filter considered in this thesis is called Sequential Importance Resampling (SIR). The idea is to approximate the target distribution p(x) by using samples drawn from a proposed distribution q(x). As q(x) is just an assumption of the desired distribution, this discrepancy needs to be accounted for by weighting each sample according to wi ∝ π(x

i₎

q(xi₎, where π(x) is a function proportional to

p(x), that can be evaluated but from which it might be difficult to draw samples. Following the steps in [31] it can be shown that the weights can be recursively computed as

wi_k∝ w_k−1i p(zk|x i k)p(xik|xik−1) q(xi k|xik−1, zk) (5.15)

(36)

and the sought after posterior filtering pdf can be approximated as p(xk|z1:k) ≈ N X i=1 wi_kδ(xk− xik) (5.16)

where δ is the Dirac delta function.

The particles are resampled according to the updated weights to avoid that all but a few particles gets neglegible weights. The steps of the particle filtering algorithm can be summa-rized as follows [29, 30, 31], assuming at the start the N particles {xi_k−1}N

i=1∼ p(xk−1|z1:k−1),

where z_k is the measurement at time k.

Prediction - Propagate the particles using the prior pdf (assumed motion model) as proposal distribution q(xi_k|xi

k−1, zk) = p(xk|xik−1). This will yield a set of propagated

particles {xi_k}N i=1.

Weight Update - Compute weights ˜wi

k for each particle. By inserting the choice of

q(.) in eq. (5.15) the weights are computed as the likelihood of a measurement evaluated for a certain particle, given by ˜w_ki = p(zk|xik). This distribution is commonly called

“sensor model” and reflects values that the state is likely to take on. The weights are normalized, with a normalized weight being denoted by w_ki.

Resample - A new set of particles {xi∗_k}N

i=1 are drawn according to the probabilities

given by the normalized weights, repeated N times. These particles will more likely be distributed according to {xi∗_k}N

i=1∼ p(xk|z1:k), which is the sought after distribution.

A major design parameter of the particle filters is the number of particles used. It is almost certain that weak convergence holds [30], meaning limN →∞p(xˆ 1:k|z1:k) = p(x1:k|z1:k), where

ˆ

p(x1:k|z1:k) denotes the estimated state distribution. An illustrative step-by-step overview of

the algorithm is shown in Figure 5.6. Particles with equal weights are propagated (row one to two) using the transition model and evaluated with the latest measurement zk using the

likelihood distribution (sensor model) p(z_k|x_k) to compute new weights (row two to three). The weighted particles are resampled changing the shape of the posterior distribution of the state, better accounting for the latest received information (row four and five). More practical aspects of implementing the particle filtering is explained in Section 6.4.4.

(37)

(xi_k,_N1) p(zk|xk) (xi k, wki) (xi∗_k,_N1) (xi_k+1,_N1)

Figure 5.6: One full iteration of particle filtering.

5.5 System overview

As the main parts of from pre-processing to detection have been presented, a complete system overview including the individual steps is presented in Figure 5.7.

(38)

UWB radar measurements Background subtraction Motion filtering Detection (CFAR/CLEAN) Grouping detections Gating Track creation Score pruning N-Scan pruning Prediction Output best tracks Post-processing Detection Tracking

Figure 5.7: Detailed view of steps included in the detection and tracking framework. The two monostatic radars send their data to the host computer on which all the dashed blocks are performed.

Detection and Tracking of Human Targets using Ultra-Wideband Radar

Examensarbete 30 hp

Oktober 2016

Detection and Tracking of Human

Targets using Ultra-Wideband Radar

Abstract

Detection and Tracking of Human Targets using

Ultra-Wideband Radar

Andreas Östman

Detection and Tracking of Human Targets using

Ultra-Wideband Radar

Populärvetenskaplig sammanfattning

Contents

List of Figures

List of Tables

1

Introduction

2

Radar

3

Data acquisition and post-processing

4

Target detection

5

Target tracking