Sensor Localization Calibration of Ground Sensor Networks with Acoustic Range Measurements

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

Sensor Localization Calibration of Ground Sensor

Networks with Acoustic Range Measurements

Examensarbete utfört i Reglerteknik vid Tekniska högskolan vid Linköpings universitet

av Viktor Deleskog LiTH-ISY-EX--12/4624--SE

Linköping 2012

Department of Electrical Engineering Linköpings tekniska högskola

Linköpings universitet Linköpings universitet

Sensor Localization Calibration of Ground Sensor

Networks with Acoustic Range Measurements

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan vid Linköpings universitet

av

Viktor Deleskog LiTH-ISY-EX--12/4624--SE

Handledare: David Lindgren

foi

Hans Habberstad

foi

Niklas Wahlström

isy, Linköpings universitet

Examinator: Fredrik Gustafsson

isy, Linköpings universitet

Avdelning, Institution Division, Department

Avdelningen för Reglerteknik Department of Electrical Engineering SE-581 83 Linköping Datum Date 2012-09-18 Språk Language Svenska/Swedish Engelska/English   Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport  

URL för elektronisk version

http://www.ep.liu.se

ISBN — ISRN

LiTH-ISY-EX--12/4624--SE Serietitel och serienummer

Title of series, numbering ISSN_—

Titel

Title Kalibrering av Sensorpositioner i Sensornätverk med Akustiska Avståndsmätningar_{Sensor Localization Calibration of Ground Sensor Networks with Acoustic Range} Measure-ments

Författare

Author Viktor Deleskog

Sammanfattning Abstract

Advances in the development of simple and cheap sensors give new possibilities with large sensor network deployments in monitoring and surveillance applications. Commonly, the sensor positions are not known, specifically, when sensors are randomly spread in a big area. Low cost sensors are constructed with as few components as possible to keep price and en-ergy consumption down. This implies that self-positioning and communication capabilities are low. So the question: “How do you localize such sensors with good precision with a feasi-ble approach?” is central. When no information is availafeasi-ble a stafeasi-ble and robust localization algorithm is needed.

In this thesis an acoustic sensor network is considered. With a movable acoustic source a well-defined and audible signal is transmitted at different spots. The sensors measure the time of arrival which corresponds to distance. A two-step sensor localization approach is applied that utilizes the estimated distances.

A novel approach in the first step is presented to incorporate more measurements and gain more position information. Localization and ranging performance is evaluated with simu-lations and data collected at field trials. The results show that the novel approach attains higher accuracy and robustness.

Nyckelord

Keywords sensor network, calibration, sensor localization, multidimensional scaling, ranging, time de-lay estimation

Abstract

Advances in the development of simple and cheap sensors give new possibilities with large sensor network deployments in monitoring and surveillance applica-tions. Commonly, the sensor positions are not known, specifically, when sensors are randomly spread in a big area. Low cost sensors are constructed with as few components as possible to keep price and energy consumption down. This implies that self-positioning and communication capabilities are low. So the ques-tion: “How do you localize such sensors with good precision with a feasible ap-proach?” is central. When no information is available a stable and robust local-ization algorithm is needed.

In this thesis an acoustic sensor network is considered. With a movable acoustic source a well-defined and audible signal is transmitted at different spots. The sensors measure the time of arrival which corresponds to distance. A two-step sensor localization approach is applied that utilizes the estimated distances. A novel approach in the first step is presented to incorporate more measurements and gain more position information. Localization and ranging performance is evaluated with simulations and data collected at field trials. The results show that the novel approach attains higher accuracy and robustness.

Acknowledgments

First I would like to give a big thanks to my supervisors at foi: David Lindgren and Hans Habberstad, for all insightful comments and help during my work with this thesis. With David’s expertise in sensor networks and signal processing and Hans’s knowledge and experience in acoustics made no problems invincible. For my supervisor Niklas Wahlström and examinator Prof. Fredrik Gustafsson at Linköping University I give thanks for the weekly feedback during my work with ideas and further improvements.

I would like to thank David, Niklas, Gustaf Hendeby (foi), and Fredrik for valu-able comments on structure, language, and technical descriptions during the edit-ing of this thesis.

Lastly, I stand in great gratitude to my family and friends who gave me endless joy and support during my studies. Big up!

Linköping, September 2012 Viktor Deleskog

Contents

Notation ix

I Background and Theory

1 Introduction 3 1.1 Problem Formulation . . . 3 1.2 Assumptions . . . 4 1.3 Related Work . . . 4 1.4 Thesis Work . . . 5 1.5 Thesis Outline . . . 5 2 Concepts 7 2.1 Sensor Network . . . 7

2.2 Sensor Localization Calibration . . . 7

2.3 Multidimensional Scaling (mds) . . . 8

2.4 mds and Sensor Localization . . . 8

2.5 Ranging . . . 9

3 Theory 11 3.1 Sound Propagation in Air . . . 11

3.2 Classical Multidimensional Scaling . . . 13

3.3 Least Squres Fitting of Two Point Sets . . . 15

II Models and Estimation

4.2 Time Delay Models . . . 22

4.3 Sensor Models . . . 22

5 Estimation 25

viii CONTENTS

5.1 Least Squares . . . 25

5.2 Fundamental Performance Bounds . . . 26

5.3 Time Delay . . . 27

5.4 Sensor Localization . . . 32

6 Signals 41 6.1 ofdm . . . 41

6.2 Chirp . . . 43

6.3 Birdsong: Black-Throated Loon . . . 44

III Results and Discussion

7.2 Datasets . . . 50

8 Results 55 8.1 Ranging Scenario . . . 55

8.2 Sensor Localization Scenario . . . 55

9 Concluding Remarks 67 9.1 Conclusions . . . 67

9.2 Future Work . . . 68

IV Appendices

A.2 Randomized Configuration . . . 74

A.3 Benchmark Configuration . . . 74

A.4 Outdoors Configuration . . . 75

B Evaluation of Signal parameters 79

Notation

Symbols

Symbol Description

sn _{Position vector of the n:th sensor}

xk Position vector of the k:th source

SN _{Set of N sensors}

XK Set of K sources

y_ki Single toa measurement at sensor i from source k

y_kij Single tdoa measurement between sensors i and j

from source k yi

k Measurement vector with M toa measurements at

sensor i from source k

yij_k Measurement vector with M tdoa measurements be-tween sensors i and j from source k

Rxx(t) Autocorrelation function (acf) Ryx(t) Cross correlation function (ccf)

Φ_xx(ω) Power spectral density (psd) Φyx(ω) Cross spectral density (csd)

E Expectation operator ~ Convolution in time domain

x[n] Discrete signal sampled with frequency fs

DFT{ · } The dft of a discrete signal || · || Euclidean 2-norm

·T _{Matrix transpose}

h(xk; θ) Nonlinear function of state vectorxkparameterized by

the parameter vector θ

x Notation Parameters Parameter Description c Propagation speedh_m s i fs Sampling frequency [Hz] B Bandwidth [Hz] fc Centrum freqency [Hz] Ts Signal duration [s]

br toarange measurement bias [m]

N Number of sensors

K Number of sources

M Number of multiple measurements

Definitions Definition Description E = RT 0 |s(t)|2dt Energy of signal s(t) SNR = E/σ2

e snrof a signal s(t) observed in wgn with variance σe2

SNRdB=

10 log₁₀SNR snrin logarithmic decibel (dB) scale

e∼ N (0, σ2_{) e is Gaussian distributed with zero mean and}

covari-ance σ2 ¯x =

1

N

PN

i=0x[i] Sample mean estimator of data vector

Notation xi Abbreviations

Abbreviation Description

acf Autocorrelation Function ccf Cross Correlation Function

cmds Classical Multidimensional Scaling crlb Cramér-Rao Lower Bound

csd Cross Spectral Density dft Discrete Fourier Transform foi Swedish Defence Research Agency gcc Generalized Cross Correlation gps Global Positioning System

lplse Linear Phase Least Squares Estimator

ls Least Squares

lse Least Squares Estimator

mc Monte Carlo

mds Multidimensional Scaling nls Nonlinear Least Squares

ofdm Orthogonal Frequency Division Multiplexing psd Power Spectral Density

rmse Root Mean Square Error scc Simple Cross Correlator snr Signal to Noise Ratio

svd Singular Value Decomposition tde Time Delay Estimation tdoa Time Difference of Arrival

toa Time of Arrival

tpi Triple Parabolic Interpolation wgn White Gaussian Noise

Part I

1

Introduction

A sensor network is a set of connected and cooperating sensors. A deployed sen-sor network can be used in various applications, often in the area of surveillance and monitoring activity around buildings and objects worthy of protection. A sensor can be any type of device that can measure signals in their environment, e.g. sound waves, ground vibrations, or light. With the development of simple and inexpensive sensors expands the possibilities for large networks with many sparsely spread sensors to achieve better coverage and precision.

If an unknown object, e.g. humans or vehicles, is detected by the sensors mea-suring the signals generated by the object. With these signals the object can be located by combining the measurements and the positions of the sensors. So the position of the sensors is central to locate the target. Different applications have different accuracy requirements on the positions. Regardless accuracy require-ment, the sensor network should provide a process to determine the positions that focuses on simplicity and effectiveness.

This chapter presents the problem formulation, assumptions, related work and the outline of the thesis.

1.1

Problem Formulation

In this thesis the problem of calibrating sensor positions in acoustic ground sen-sor networks is studied. The main goal is to develop a calibration method that from erroneous range measurements estimates the positions of the sensor nodes independent of any prior knowledge of the positions. The problem can therefore be divided into two sub-problems:

4 1 Introduction

(i) Estimate range at far distances from time delay measurements based on acoustic wave propagation.

(ii) Determine the sensor positions, without any prior information of the net-work topology, by using erroneous range estimates.

Both algorithms should be evaluated by simulations and by data collected at field trials.

1.2

Assumptions

To keep the work into the time frame of the thesis some assumptions have been made to fulfil the time constraint:

• The sensor nodes are synchronized.

• The atmosphere is considered as a homogeneous acoustic channel with known propagation speed.

• The acoustic source is considered to be stationary while emitting the acoustic signal.

• An emitted signal from the source can be sensed by all sensor nodes in the network.

1.3

Related Work

The problem of sensor localization have been investigated since early 2000’s. In [Patwari et al., 2005], IEEE Signal Processing Magazine, various algorithms and fundamental issues in the area are summarized and discussed.

In [Moore et al., 2004] the concept of robust quadliterals is introduced to localize sensors with noisy sensor-to-sensor distances using trilateration. Sensors are di-vided into clusters and are able to be localized up to a global rotation and transla-tion which fits in our problem setting. The localizatransla-tion algorithm is divided into multiple phases. This approach is showed to be sound where the first phase finds rough positions and the second refines them to reduce effects of measurement noise using numerical optimization. In [Shang et al., 2004] a algorithm based on Multidimensional scaling (mds) is applied to localize the sensors instead of trilateration. The problem is formulated as least squares problem and solved by Singular value decomposition (svd). Finding the svd is a fast operation com-pared to ordinary iterative least squares solvers.

The mds approach is used in this thesis but a novel extension is presented to lower the influence of measurement noise by incorporating more measurements that are available in our problem setting.

In this thesis the acoustic propagation delay is used to estimate range. The con-cept of time delay estimation is introduced in [Knapp and Carter, 1976] which

1.4 Thesis Work 5 formulates the problem and presents and derives commonly applied time delay estimation algorithms. In [Gustafsson et al., 2010] a Orthogonal frequency divi-sion multiplexing (ofdm) modulated signal is used by a time delay estimation algorithm that exploit the connection between phase and delay to estimate the time delay with least squares theory. Benefits gained with the algorithm pre-sented in [Gustafsson et al., 2010] are sub-sample resolution and quality metric. In [Ahlström and Falk, 2001] Triple parabolic interpolation (tpi) is presented to achieve sub-sample resolution when integer sample estimation algorithms are used. Experiments of outdoors acoustic time delay estimation is presented in [Ash and Moses, 2005] where interesting signal parameters and their impact on estimation performance are evaluated.

Sound propagation is covered in [Jönsson and Nilsson, 2007] and [Blackstock, 2000] where fundamental concepts of sound waves and physical acoustics are presented and derived.

1.4

Thesis Work

This report is a part of a master’s thesis at the department of Electrical Engineer-ing (isy) at LinköpEngineer-ing University under Prof. Fredrik Gustafsson at the division of Automatic Control. The work has been performed at the Swedish Defence Research Agency (foi) under supervision of Dr. David Lindgren, and Hans Hab-berstad.

This thesis work is a part of foi Centre for Advanced Sensors, Multisensors and Sensor Networks (focus) funded by the Swedish Governmental Agency for Inno-vation Systems (vinnova), and the Knowledge Foundation (KK-stiftelsen).

1.5

Thesis Outline

This thesis is divided into different levels where the top levels: • Background and Theory

– Introduction – Concepts – Theory

• Models and Estimation – Models

– Estimation – Signals

• Results and Discussion – Data Collection

6 1 Introduction

– Results – Summary

are the three main cornerstones. This structure is used to give a better overview of the distinct parts of thesis. Background and Theory a brief introduction to the problem, concepts, and the interesting theory behind existing methods and derivations are covered. Later in Models and Estimation the model and sensor definitions, and derived estimators are covered. In Results and Discussion the data collection, results, and summary are covered.

Two appendices are enclosed in this thesis covering sensor network configura-tions and evaluation of signal parameters.

2

Concepts

This chapter introduces some important concepts, e.g sensor networks, sensor localization, mds, and ranging, that are applied and investigated in the thesis.

2.1

Sensor Network

A Sensor Network (sn) consists in general of multiple connected nodes [Gustafs-son, 2010]. The nodes can communicate either by wire, or wireless communica-tion channels dependent on applicacommunica-tion. A node can either be a sensor (receiver) or a target (transmitter). When operational, the sensors detect a target and give some type of observation of it, a measurement.

A measurement can be some sort of range (toa, tdoa), bearing (aoa), or energy (RSS) related to the target. A common application of a sensor network is to esti-mate a position of an observed target by utilizing the measurements, [Lindgren et al., 2010].

2.2

Sensor Localization Calibration

Calibration is the problem of finding or estimating the unknown configuration parameters of the sensor nodes in a sensor network [Gustafsson, 2010]. With sen-sor localization calibration, the unknown position parameters of the sensen-sor nodes are estimated. Calibration is important in a sensor network application because it has to be performed before it can be used for its proposed task. Calibration often involves dedicated experiments with known premisses where parameters are estimated or measured.

8 2 Concepts

Sensor localization calibration can be performed in numerous ways where sensor capabilities and resources define the limits. By equipping sensors with a simple Global Positioning System(gps) receiver the absolute position can be determined with accuracy up to 5 - 15 m. To equip all sensors with such receiver, is a both costly and energy consuming way to achieve accurate positions. Other scenarios include reference nodes which are nodes with known positions [Patwari et al., 2005]. If the other sensors can relate their measurements to the available refer-ence nodes they can estimate their positions using the known absolute positions. In [Moore et al., 2004, Sottile and Spirito, 2008] methods are presented that rely on inter-sensor communication where each sensor can measure the distance be-tween neighbouring sensors and with trigonometric relations calculate their po-sition.

The positions that are estimated with a localizing algorithm can either be defined in a world fixed coordinate system or in a locally fixed coordinate system. De-pending on prior knowledge of the sensor network configuration the estimated positions can be correct up to a arbitrary point transformation (translation and rotation/reflection).

2.3

Multidimensional Scaling (

MDS

)

Multidimensional Scaling (mds) is a family of data analysis methods. These methods calculate a low-dimensional representation of high-dimensional data by using a defined model. It started in the 1950s with Torgerson who introduced the general topic of mds applied in psychometrics1[Young, 1987]. His method was the first systematic algorithm that provided a multidimensional map of points from erroneous inter point Euclidean distances, known as metric mds.

A mds model can be defined in different spaces. In general the model is a simple algebraic equation with some geometric interpretation, e.g. Euclidean distance in a Euclidean space. The essence of mds is to visualize complex data in a spatially representation. The data can be represented by a matrix whose elements indicate some relation to each other, e.g. correlations, distances, and similarities.

2.4

MDS

and Sensor Localization

In this thesis the derived mds method Classical Multidimensional Scaling (cmds) [Torgerson, 1965] will be used as basis to estimate sensor positions utlilzing erre-nous range measurements [Shang et al., 2004]. Localization problems are often formulated as a standard optimization problem. To solve such problems itera-tive solvers are generally used where the initial start parameters are crucial to guarantee convergent solutions. With cmds it can instead be solved in closed form.

2.5 Ranging 9

2.5

Ranging

Ranging is the problem of estimating the distance between two points, e.g. be-tween a transmitter and a receiver. Here we will estimate the time delay when an acoustic signal is transmitted from an acoustic source and received by an acous-tic sensor through air. With the estimated time delay the range can be found by multiplying with the propagation speed of sound in air. Various methods on Time Delay Estimation (tde) exist where [Knapp and Carter, 1976] covers the most used variants which is based on estimating the Cross-Correlation Function (ccf) between the transmitted and the received signal. In [Gustafsson et al., 2010, Chan et al., 1978] the estimated Cross Spectral Density (csd) between the trans-mitted and received signals is used to exploit the relationship between phase shift and time delay.

3

Theory

This chapter presents and describes the theory and methods behind the devel-oped algorithms. By studying the theory it will be easier to comprehend the algorithms and to discover potential improvements and constraints.

3.1

Sound Propagation in Air

Sound is longitudinal waves that propagate through a material medium, often a gas. The speed of sound, denoted c, that we are interested in is the speed at which this wave propagates through air. From the acoustic wave equation [Jönsson and Nilsson, 2007, p. 96] ∂2s ∂t2 = 1 κρ ∂2s ∂x2 (3.1)

we can see from the general wave equation that the propagating speed is actually

c2= 1

κρ =⇒ c =

1

√_κρ (3.2)

where κ is the compressibility factor and ρ is the density of the gas. The com-pressibility factor can be defined by assuming that sound waves propagate as an adiabatic process. During an adiabatic process the entropy in the medium remains constant, isentropic. If we assume that the energy losses are negligible in air, which is common in acoustics, the entropy is constant [Blackstock, 2000, p. 32-33]. For an adiabatic process it follows that

P Vγ= k (3.3)

12 3 Theory

Table 3.1:Coefficients for air Coefficient Value

γ 1.40

M 0.029 kg/mol

R R0

M = 287 J/kgK

where k is a constant, P is pressure, V is molar volume and γ is the ratio of specific heats. By writing (3.3) as

P = k · V−γ

we can take the derivative of P with respect to V as

∂P

∂V = −γ · k · V−γ−1

and replace the constant gives

∂P

∂V = −γ · k · V−γ−1= −γ · P V

γ_{· V}−γ−1_{= −γP V}−1 _(3.4)

By definition the compressibility factor is defined as [Jönsson and Nilsson, 2007, p. 97]

κ≡ −_V1 ∂V_∂P = −_V1 _{−γP V}1 ₋₁ !

= 1

γP (3.5)

where the derivative ∂V /∂P is replaced with the one in (3.4). So the sound speed in a gas is according to (3.2) c = _√1 κρ = s γP ρ (3.6)

Through experiments the expression in (3.6) holds so the assumptions above are correct. By applying the perfect gas law

P = RρT (3.7)

where T is the absolute temperature and R is the universal gas constant divided by the molecular weight of air the speed of sound in air can be simplified

c = s γP ρ = P ρ = RT  =pγRT . (3.8)

In this application the gas is air so we know the constant values of γ and R, see Table 3.1. One interesting detail is that sound speed in a gas varies as the square root of absolute temperature.

3.2 Classical Multidimensional Scaling 13 101 102 103 104 105 10−4 10−2 100 102 104 Frequency [Hz] Sound Absorption [dB/100*m]

Figure 3.1: Sound absorption in air as a function of frequency at20◦C per 100 m, [Blackstock, 2000]

and by the absorption of sound in the atmosphere. In [Blackstock, 2000, p. 513-516] the absorption coefficient α [nepers/m] in the atmosphere is described to be dependent on frequency. In Figure 3.1 the absorption coefficient α [dB/100 · m] is showed as a function of frequency. So signals with high frequency content gets attenuated because of sound absorption in the atmosphere. If you are near a strike of a lightning during a thunderstorm you both hear a low-frequency boom and high-frequency sound but at far distances you just hear the boom. So the high-frequency sound is gone because of different sound absorptions.

3.2

Classical Multidimensional Scaling

In the cmds framework we have distance data which we want to visualize in a low-dimensional geometric representation, often 2-D or 3-D. The data is Eu-clidean distancesD ∈ RN×N _{amongst nodesX = (x1, . . . , xN)}T_{. The nodes are}

Cartesian points xi = (x1, . . . , xm)T ∈ Rm in a m-dimensional Euclidean space.

The information from the distances D is sufficient to exactly reconstruct the nodesX up to translation, rotation, and reflection [Torgerson, 1965].

First assume that the nodes inX are known, then the outer product, B ∈ Rm×m_,

ofX is given by

B = XXT _(3.9)

which can be connected to the distances inD as

d_ij2 = (xi− xj)T(xi− xj) =xTi xi− 2xTi xj+xTjxj= bii− 2bij+ bjj (3.10)

where the matrix elements are defined as

dij := ||xi− xj|| bij :=xTi xj.

14 3 Theory

So if the outer product matrixB can be calculated from the distances D, the nodes could be found by factoringB as in (3.9). In [Torgerson, 1965] it is shown that B can be calculated by double-centering D2_{, whereD}2_{is the element-by-element}

squared distance matrix ofD. A matrix is double-centered as follows: for each element subtract column and row mean, adding the grand mean, and divide by -2. Each element bij inB is calculated as

bij = −1₂                      d_ij2 −_N1 N X k=1 d_ik2 | {z } Row mean −_N1 N X k=1 d_kj2 | {z } Column mean + 1 N2 N X g=1 N X h=1 d2_gh | {z } Grand mean                      . (3.12)

The nodesX can then be recreated through the (svd) of B. B = U ΣVT _{= U ΣU}T

ˆX = U√Σ (3.13)

where U = V as B is symmetric. U is an orthogonal unitary matrix and Σ is a diagonal matrix containing the singular values Σ = diag(σ1, σ2, . . . , σN). The

singular values are ordered as σ1 ≥ . . . ≥ σN. In (3.13) the solution ˆX have

dimen-sion N because of the high-dimendimen-sional distance data in D. From the beginning we know that the distances are generated from points in a m-dimensional space. The best lowrank approximation, in Least squares (ls) sense, of the original N -dimenional solution is given by using the m first singular values and left singular vectors inB as Bm= UmΣmUmT ˆX = Um p Σm (3.14) where Um contains the m first left singular vectors and Σm = diag(σ1, . . . , σm)

contains the m first singular values. The quality of the low-dimensional approxi-mation can be calculated with

gof= Pm i=1|σi| PN i=1|σi| (3.15) which is the Goodness of fit (gof). The solution ˆX in (3.14) gives the best approx-imation

Bm= ˆX ˆXT (3.16)

ofB in ls-sense [Torgerson, 1965]. The method is summarized in Algorithm 1. Example 3.1: cmds and Swedish cities

One simple example to demonstrate possible applications of cmds is to find the locations of cities by the distance between them. Here we will apply this on the distances between six cities in Sweden: Stockholm, Gothenburg, Malmö, Linköping, Jönköping, and Kiruna. Table 3.2 shows the distances between the cities and with cmds their 2-D/3-D locations can be calculated.

3.3 Least Squres Fitting of Two Point Sets 15 Algorithm 1 Classical Multidimensional Scaling, cmds

1: Form distance matrixD from known or estimated distances between nodes. 2: Calculate the squared double centered matrixB from D.

3: Calculate the svd ofB, B = U ΣUT.

4: Keep the m largest singular values and their singular vectors and define them as Σmand Um.

5: ˆX = U_m√Σ_mare the recreated nodes with dimension m.

Figure 3.2 shows the recreated map of the cities and the goodness of fit which gives the best fit when m ≥ 2 (2-D points), so two dimensions is sufficient to represent the distances as Euclidean points, which is intuitive because the differ-ence in elevation between the cities is much lower than the distances amongst them. Note that a one-dimenional solution would also represent a rather good fit because Sweden is more elongated than wide.

The recreated map has been rotated to a consistent view of the true map due to the nature of cmds.

Table 3.2: Distances between cities in kilometers, calculated from decimal degrees (wgs84)

Stockholm Gothenburg Malmö Linköping Jönköping Kiruna

Stockholm 0 397.9 513.3 174.8 284.7 955.8 Gothenburg • 0 242.2 228.8 131.2 1204.1 Malmö • • 0 349.0 252.0 1414.5 Linköping • • • 0 109.9 1078.2 Jönköping • • • • 0 1163.1 Kiruna • • • • • 0

3.3

Least Squres Fitting of Two Point Sets

Consider two point setsX = (x1, . . . , xN) andY = (y1, . . . , yN) where the points

x and y are two-dimensional Cartesian points x, y ∈ R2×1_{. The relation between}

the points can be described by

yi = Rxi+ T + ei (3.17)

where R ∈ R2×2 _{is a transformation matrix (rotation or reflection), T ∈ R}2×1 _is

translation vector, and ei ∈ R2×1 is a noise vector. The goal is find the Least

squares estimator(lse) of R and T  (R, T ) = arg min R,T N X i=1 ||yi − (Rxi+ T )||2 . (3.18)

16 3 Theory Stockholm Gothenburg Malmö Linköping Jönköping Kiruna

(a)Recreated map

1 2 3 4 5 6 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 1.005 Dimensions GOF (b)Goodness of fit

Figure 3.2:Recreated map of Swedish cities and gof

The lse is often found by iterative algorithms but instead a closed-form solution involving the svd can be used [Arun et al., 1987].

If ˆR and ˆT are the lse to (3.18) then Y and ˆY = (ˆy1, . . . , ˆyN), where ˆyi = ˆRxi, have

the same centroid, i.e

¯y = ¯ˆy (3.19) where ¯y := 1 N N X i=1 yi (3.20) ¯ˆy := 1 N N X i=1 ˆyi = ˆR¯x + ˆT (3.21) ¯x := 1 N N X i=1 xi. (3.22) Denote ˜xi :=xi− ¯x (3.23) ˜yi :=yi − ¯y . (3.24)

3.3 Least Squres Fitting of Two Point Sets 17 Algorithm 2 Finding ˆR using the svd

1: Calculate the "covariance" matrix H ∈ R2×2

H =: N X i=1 ˜xi˜yTi (3.27) 2: Find the svd of H H = U ΣVT (3.28)

3: Calculate the lse ˆR as

ˆR = V UT _(3.29)

The minimum in (3.18) can instead be reformulated as ˆR = arg min R N X i=1 ||˜yi− R˜xi||2 (3.25)

where ˆR is the optimal transformation matrix in ls sense. [Arun et al., 1987]. The original ls problem is therefore divided into two parts:

1) Find the lse ˆR that minimizes (3.25) with Algorithm 2. 2) The translation ˆT is found by

ˆ

T = ¯y− ˆR¯x . (3.26)

A closed-form solution to the problem in (3.18) is presented involving the svd. The found orthogonal matrix ˆR can potentially contain reflections in addition to rotations. If a correct rotation matrix is required some changes in Algorithm 2 can be implemented, see [Arun et al., 1987], but is not considered in this thesis.

Part II

4

Models

With a model, we intend to find a mathematical relationship that describes a process that we want to investigate or predict its future behavior. Sensors receive signals that are emitted from an acoustic source causing a time delay of the signal. From these delayed signals we intend to find the propagation delay to estimate the travelled distance. Therefore it is required that three types of models needs to be defined to simulate the processes:

Signal model Describes the properties of the emitted signal from the acoustic source.

Time delay models Describes how the transmitted signal is delayed when trav-eling from source to sensor.

Sensor models Describes how a sensor interprets time delay measurements into distances in the presence of simulated measurement noise.

4.1

Signal model

A signal, s(t), is emitted from an acoustic source with specific parameters. The frequency content of the signal is defined in the frequency set

Ω∈ [f_c− B/2, f_c+ B/2] (4.1)

and non-zero in the interval [0, Ts] where

• B: Bandwidth, [Hz] • fc: Center frequency, [Hz]

• Ts: Signal duration, [s]

22 4 Models

The energy E of a signal is finite and is calculated as E =

Ts Z

0

|s(t)|2dt (4.2)

where the physical unit is discarded. The sampled version of the signal s(t) is defined as s[n], where it is sampled with the sampling frequency fs ≥ 2(fc+ B/2)

(Nyquist rate).

4.2

Time Delay Models

An acoustic signal, s(t), is emitted from a remote source and observed by a spa-tially separated sensor in the presence of white Gaussian noise (wgn) e(t) during the time T = Ts+ Tmax where Tmax is the maximum considered time delay. The

delayed signal can then be modelled as

y(t) = as(t− τ) + e(t), 0 ≤ t ≤ T (4.3) where τ corresponds to the time delay and a is a possible attenuation of the signal. The signal s(t) is assumed to be uncorrelated with the noise e(t), E s(t)e(t) ≡ 0. We can also extend this model to cover the case when a signal is received at two spatially separated sensors as

y1(t) = s(t − τ) + e1(t), 0 ≤ t ≤ T (4.4a)

y2(t) = st − (τ + δ) + e2(t), 0 ≤ t ≤ T (4.4b)

where δ is the difference in time arrival between the sensors during the time

T = Ts+ Tmax. The signal s(t) is assumed to be uncorrelated with the two noise

terms, e1(t) and e2(t).

4.3

Sensor Models

A sensor gives some type of observation of a target xk which is interesting to

us. A typical sensor output yk is modelled by a function h( · ) of a target xk and

parameterized by a set of parameters θ in the presence of noise ek which is wgn

with zero mean and covariance σ2

e

yk = h(xk; θ) + ek, e∼ N (0, σe2) (4.5)

where the parameters can consist of necessary properties of the sensor, e.g sensor position and measurement bias. The target xk can generate multiple

measure-ments at the sensor which will be denoted with boldface yk = {yk}M1 ∈ R1×M,

4.3 Sensor Models 23

4.3.1

Time of Arrival (

)

Using the time delay model in (4.3) an acoustic signal from sourcexkis modeled

to reach sensorsi _{at time}

hT OA(xk;si, br) = 1_c||xk− si|| + br (4.6)

where the model is parameterized by sensor positionsi_{and measurement bias b} r

and c is the known propagation speed. A toa sensor then gives timing measure-ments as

τ_ki = hT OA(xk;si, br) + eik, eik ∼ N (0, Rik) (4.7)

where ei

kis independent wgn. This toa sensor can be converted to a range sensor

by multiplying the timing τi

kwith the propagation speed c r_ki = c · τi

k = c · hT OA(xk;si, br) + vki, vik∼ N (0, c2Rik) (4.8)

where vi

k is wgn. Figure 4.1a shows the ranging measurements rki and r j k from

two toa sensors.

4.3.2

Time Difference of Arrival (

)

Using the time delay model in (4.4) the difference in time arrival between sensors si_andsj_{when a signal is transmitted from sourcex}

kwill be on the form hT DOA(xk;si, sj) = 1_c ||xk− si|| − ||xk− sj||



(4.9) where the model is parameterized by the two sensor positionssi _andsj_{. A tdoa}

sensor then gives timing measurements as

τ_kij = hT DOA(xk;si, sj) + eijk, e ij k ∼ N (0, R ij k) (4.10) where e_kij = ei k − e j

k is wgn. This tdoa sensor can also be converted to a range

sensor by multiplying the timing τ_kij with the propagation speed c

r_kij = c · τ_kij = c · hT DOA(xk;si, sj) + vij_k, v_kij ∼ N (0, c2Rij_k) (4.11)

where v_kij is wgn. Figure 4.1b shows the ranging measurement r_kij from a tdoa sensor.

24 4 Models xk si sj r_ki r_kj

(a) Two toa sensors outputs the range measurementsr_ki andr_kj.

xk

si

sj r_kij

(b)A tdoa sensor outputs the range measurementr_kijwhich is the time-difference of arrival be-tween sensorssiandsj.

Figure 4.1:Range measurements from two toa sensors and one tdoa sensor generated by targetxk.

5

Estimation

This chapter presents in detail the developed and applied estimation algorithms; ranging and sensor localization, and a brief introduction to the least squares (ls) estimation framework and fundamental estimation bounds.

5.1

Least Squares

When to find the parameters of a mathematical model the theory of system iden-tification can be applied. System ideniden-tification is to find the model parameters which best describes the observed data [Gustafsson et al., 2010, chapter 6]. A very common model structure is the linear Gaussian model

yk = φTkθ +ek, ek∼ N (0, R), k = 1, 2, . . . , N (5.1)

where ykis the observed data at time instant k, φkis the known regression vector, θ is the unknown model parameters andekis wgn with variance R. So the goal is

to estimate the unknown parameters θ given the observed data. The least squares methodminimizes the cost function

V (θ) = 1 N N X k=1 (yk− φTkθ)2= 1 N N X k=1 k(θ)2 (5.2)

where k(θ) is the residual or model error at time instant k. The lse ˆθ can be

formulated as

ˆθ = arg min

θ

V (θ) (5.3)

where the operator arg min is the minimizing argument. So ˆθ is the estimate of

θ that minimizes the cost function in (5.2). To keep the problem compact it can

26 5 Estimation

instead be formulated with vector notation YN =                y1 y2 .. . yN                ,ΦN=                φT₁ φT₂ .. . φT_N                , EN =                e1 e2 .. . eN                (5.4) so (5.1) can be written as YN =ΦNθ +EN (5.5)

With vector notation the ls cost function takes the form of

Vls (θ) = N X k=1 (yk− φTkθ)T(yk− φTkθ) = (YN− ΦNθ)T(YN− ΦNθ) (5.6)

and the lse

ˆθls = arg min θ Vls (θ) = (ΦT NΦN)−1ΦTNYN (5.7)

which is the minimizing argument of (5.6). ls problems can efficiently be solved using numerical methods discussed in [Gustafsson et al., 2010, p.224-225] where the QR factorization is suggested instead of calculating the matrix inverse in (5.7).

5.2

Fundamental Performance Bounds

In estimation theory the Cramér-Rao Lower Bound (crlb) expresses a lower bound on the variance of any unbiased estimator [Kay, 1993, p. 27]. An unbi-ased estimator ˆθ has the property

E[ ˆθ] = θ0 _θ0_{∈ Θ (5.8)}

where Θ denotes the range of possible values of the true parameters θ0_{. In}

prac-tice the crlb provides a benchmark on the theoretical lower bound of any unbi-ased estimator.

The crlb is given by the Fisher Information Matrix (fim) I(θ) [Gustafsson, 2010, p. 81]. With the assumption of Gaussian measurement errors e ∼ N0, R the fimequals

I (θ) = JT_(θ)R−1_J(θ)

J(θ) =∇θh(θ)

(5.9) where ∇θh(θ) is the Jacobian evaluated at θ. The crlb is given by

Cov( ˆθ) = E(θ0− ˆθ)(θ0− ˆθ)T _{≥ I}−1_(θ0_{) (5.10)}

When investigating localization performance the Root Mean Square Error (rmse) gives a practical measure of the mean position error in meters. In two dimensions

5.3 Time Delay 27 the achievable position rmse is given by

rmse= q Eh(θ0 1− ˆθ1)2+ (θ02− ˆθ2)2i = q tr Cov( ˆθ) ≥qtrI−1_(θ0_{) (5.11)}

which is implied from the crlb.

When multiple independent measurements are available their corresponding in-formation matrices can be added because inin-formation is additive [Gustafsson, 2010, p. 81] I1:M = I1+ I2+ . . . + IM= M X k=1 Ik Cov( ˆθ) ≥ I−1 1:M (5.12) and with more information the better estimation accuracy.

5.3

Time Delay

From the background on ranging in Section 2.5 two commonly and easy to im-plement estimation methods are chosen to estimate the time delay parameter in the models defined in Section 4.2. The first method is based on the Cross Corre-lation Function(ccf) which is the Simple Cross Correlator (scc) and the second is based on the Cross Spectral Density (csd) which is the Linear Phase Least Squares Estimator(lplse).

5.3.1

Cross Correlation Method - SCC

With the time delay model (4.3) we can compute the ccf

Rys(t) = E y(l)s(l − t) = Eh [as(l − τ) + e(l)] s(l − t)i = (5.13)

aRss(t − τ) + E independent z }| { e(l)s(l− t) | {z } =0 = aRss(t − τ) = aRss(t) ~ δ(t − τ) (5.14)

where δ(t) is the Dirac delta function and ~ denotes convolution in time domain. So it turns out that the ccf is a time shifted attenuated copy of the Autocorrelac-tion funcAutocorrelac-tion(acf) Rss(t). At time t = τ the ccf should take its maximum value.

Due to sampled signals and finite observation time we must estimate the ccf [Gustafsson et al., 2010, p. 96] ˆRys[k] = 1 N N−k−1 X m=0 y[m + k]s[m], k = 0, 1, . . . , N− 1 (5.15)

28 5 Estimation −5000−1 0 5000 −0.5 0 0.5 1 Lag [samples] Normalized correlation 0 50 100 150 200 −1 −0.5 0 0.5 1 Lag [samples] Normalized correlation

Figure 5.1:Estimated ccf of the transmitted and received signal. Left: Dur-ing the whole duration of received signal. Right: Zoomed in around the maximum value atτ = 100.

where we observe N samples sampled with frequency fs. We can then form the

time delay estimate ˆτ as

ˆ∆ = arg max

k

ˆRys[k]

ˆτ = fsˆ∆

(5.16) where ˆ∆ is the estimated delay in samples. (5.16) is the Simple Cross Correlator (scc) estimate. In theory the wgn term disappears due to independence w.r.t. the transmitted signal but in a real application that is not the case. A band pass filter applied on the received signal y would increase the snr and just keep the signal information in the frequency set Ω.

Example 5.1: scc Estimator

A signal s(t) with bandwidth B = 500 Hz and center frequency fc = 1 kHz is

transmitted to a receiver. The receiver observes the signal with snr 10 dB as the BP-filtered signal y(t). The true delay between the transmitter and the receiver is 100 samples. With scc the delay is estimated to ˆτ = 100 samples which is the true delay. Figure 5.1 shows the ccf ˆRys[k].

Sub-sample Resolution

In some cases the maximum value of the ccf is not attained at the true delay since the resolution of the delay estimate ˆτ in (5.16) is limited to multiples of 1

f_s.

In [Ahlström and Falk, 2001] the method Triple Parabolic Interpolation (tpi) is presented to interpolate the ccf around the maximum peak. The method fits a quadratic polynomial with three points, the peak ( ˆτ) and the surrounding two points. The interpolated ccf around the found peak at ˆ∆ is modelled as

5.3 Time Delay 29 98 98.5 99 99.5 100 100.5 101 101.5 102 0.92 0.93 0.94 0.95 0.96 0.97 0.98 Lag [samples] Normalized correlation

Figure 5.2:Interpolated ccf with tpi. Solid: ccf. Dotted: ccf +tpi

where the unknown parameters in θ are estimated using ls. The sub-sample estimate ˆ∆tpican is estimated by setting the derivative of ˆRys,tpiw.r.t. the sample

lag k to zero d dk ˆRys,tpi[k]     _{k= ˆ∆tpi} = ˆθ2+ 2 ˆθ3ˆ∆tpi= 0 =⇒ ˆ∆tpi= − ˆθ2 2 ˆθ3 ˆτtpi= fsˆ∆tpi (5.18) where ˆτtpiis the scc + tpi time delay estimator.

Example 5.2: Sub-sample Resolution

Consider the same signal parameters as in Example 5.1 but with 1001₃ samples delay. The standard scc estimates the delay to ˆ∆ = 100 samples but with scc + tpi we get the correct delay ˆ∆tpi= 1001₃samples. Figure 5.2 illustrates the ccf

(solid) where circles represent integer samples and the tpi interpolated version (dotted) where the black circle is the estimated sub-sample delay ˆ∆tpi.

5.3.2

Linear Phase Method - LPLSE

Still consider the time delay model (4.3). By taking the Fourier transform of the ccfin (5.13) we get the theoretical Cross Spectral Density (csd)

Φ_ys(ω) = ae−jωτΦ_ss(ω) (5.19)

between the signals y(t) and s(t). With the csd, an estimate of the time delay can be provided from relations in the frequency domain. The method is based on the description made in [Gustafsson et al., 2010, p. 111,232-235]. The phase of the csd Φys(ω) should be a straight line with constant slope τ. From the signal

model described in Chapter 6 which are defined in a bounded frequency set Ω with center frequency fcand bandwidth B the phase of the csd can be modeled

30 5 Estimation as ϕ(ω) = arg(Φys(ω)) =        φ0− ωτ + eφ(ω) if ω ∈ Ω eφ(ω) otherwise (5.20) where φ0is the random initial phase and eφis additive phase noise which comes

from the noise term in (4.3). The model describes a linear phase curve in the fre-quency set Ω with an initial offset and random outside. The signals are sampled so instead we use a discrete phase model

ϕk= arg  Φ_ys[ω_k] =        φ0− ωk∆+ eφ[ωk] if ωk ∈ Ω eφ[ωk] otherwise (5.21) where ∆ is the time delay in samples and

Φ_ys[ω_k] = Y [ω_k]S∗[ω_k]

Y [ωk] = DFTny[n]o S[ωk] = DFTns[n]o

ωk = 2πf_Nsk , k = 0, 1, . . . , N− 1

(5.22)

from N signal observations. With the linear model the lse of the parameters can be found using the linear estimation framework in Section 5.1

YN=ΦNθ YN =               ϕ1 ϕ1 .. . ϕN               ,ΦN=               1 −ω1 1 −ω2 .. . ... 1 −ωN               , θ ="φ0 ∆ # (5.23)

where ω1, . . . , ωN ∈ Ω. The Linear Phase Least Squares Estimator (lplse) can

then be calculated as

ˆθ = (ΦT

NΦN)−1ΦTNYN

ˆτ = fsˆ∆ .

(5.24) The variance of the time delay estimate ˆτ can be calculated from ls theory [Gustafs-son et al., 2010, p. 233-234] Var( ˆτ) ≈ Bσe2/2 aB4/12 = 6σ2 e aB3 (5.25)

where a is possible signal attenuation or amplitude. Example 5.3: lplse Estimator

A signal s(t) with bandwidth B = 500 Hz and center frequency fc = 1 kHz is

transmitted to a receiver. The receiver observes the signal with snr 10 dB as the signal y(t). The true delay between the transmitter and the receiver is 100 sam-ples. With lplse the delay is estimated to ˆ∆ = 99.77 ± 0.41 samsam-ples. Figure 5.3

5.3 Time Delay 31 0 1000 2000 3000 4000 5000 −80 −60 −40 −20 0 20 40 Frequency [Hz] Phase [rad] 800 900 1000 1100 1200 −4 −2 0 2 4 Frequency [Hz] Phase [rad]

Figure 5.3: Phase of the csd. Left: Φys[ωk] phase. Right: Measured (thin)

and estimated linear phase (thick) in frequency set Ω.

shows the measured and the estimated linear phase of the csd.

5.3.3

Estimation Performance

To study the efficiency of the three estimators: scc, scc + tpi, and lplse they will be compared against the derived crlb in (5.26) when estimating a unknown parameter of a deterministic signal observed in wgn [Kay, 1993, p. 53-56]

Var( ˆτ) ≥ σe2 EF2 = 1 SNR · F2 (5.26) where SNR := E σe2 (5.27) is the signal to noise ratio (snr) defined as the quotient of signal energy and noise variance. F2_{is the Mean square bandwidth (msb)}

F2:= R∞ −∞(2πf )2|S(f )|2df R∞ −∞|S(f )|2df (5.28) of the signal s(t). S(f ) is the Fourier transform of s(t) and f is continuous-time frequency. The expression in (5.26) can be rewritten as [Gustafsson, 2010, p. 400]

Var( ˆτ) ≥ 1 8π2_BT

sfc2·SNR

(5.29) which is easier to interpret by the signal parameters. From (5.29) we see that the lower bound is inversely proportional to snr and the signal parameters of s(t). To compare the theoretical variance of the lplse in (5.25) with crlb in terms of signal snr one can recall that the energy of a spectrum with constant amplitude

32 5 Estimation 0 10 20 30 40 50 10−20 10−15 10−10 10−5 100 SNR [dB] Variance [s 2 ] CRLB LPLSE

Figure 5.4: crlband calculated lplse variance with respect to snr. Signal: ofdmwith parameters:f_c= 1000 Hz, B = 500 Hz, and T_s = 0.3 s.

a and bandwidth B isE = aB. The signal snr can then be calculated as SNR = E

σe2

= aB

σe2

, (5.30)

so it turns out that (5.25) can be approximated as

Var( ˆτ) ≥_B2_·6_SNR (5.31)

where the amplitude of signal s(t) is considered as constant within its defined frequency set Ω.

Figure 5.4 illustrates the crlb and the approximative variance of the lplse in (5.25) with respect to snr. We see that the lplse does not attain the crlb despite high snr.

5.4

Sensor Localization

Consider a sensor network with N sensors SN _{= {s}n_}N

1, wheresn= (x1, x2)T ∈ R2,

spread over an area A. We also have K acoustic sources XK = {xk}1K, wherexk =

(x1, x2)T ∈ R2, around the area A. Neither the positions of sensors nor sources

are known.

The sensors measure range to a source, toa timing, and range difference between two sensors, tdoa timing, as described in Section 4.3.

The parameters that are of highest interest to estimate are the sensor positionsSN

while the source positionsXk and the toa measurement bias br are considered

5.4 Sensor Localization 33 r_ki r_kij toa tdoa ˆSN ˆSN ˆX K ˆXK ˆb_r initial step refinement step

Figure 5.5: Localization algorithm workflow from measurements to esti-mated parameters.

Due to the nature of the problem structure with many parameters and nonlin-ear relationships the algorithm needs to be divided into two steps. Commonly applied localizing algorithms relies on nonlinear minimization which is often solved by numerical iterative solvers, e.g. Gauss-Newton and Newton-Raphson. Iterative solvers requires a initial guess of the unknown parameters. If no param-eters are known the solver tends to diverge or find a local minima. But if the initial guess is close enough to the true solution a iterative solver would in the best case find the global minima, which corresponds to the true solution. The sensor localization algorithm will therefore be divided into two steps:

Initial step Calculates an estimate that is close to true solution.

Refinement step Refines the initial estimate with nonlinear optimization. The initial step will estimate positions that are in a sufficient neighbourhood, initial estimate, of the true solution by applying a extended approach using the cmdsalgorithm. The refinement step uses the initial estimate as prior informa-tion and with numerical optimizainforma-tion find the global soluinforma-tion, refined estimate. Because of the nature of the problem the final solution is defined in a locally co-ordinate system which is correct up to arbitrary transformation and translation. Instead of applying the original approach of cmds, described in [Shang et al., 2004], a novel approach will be presented, called the novel cmds approach, which is a own contribution. The novel approach result in more accurate estimates of the sensors positions and provides more information to integrate by refinement step estimation algorithms.

In Figure 5.5 the whole sensor localization workflow can be studied which inte-grates both toa and tdoa measurements divided into two localization steps.

34 5 Estimation

5.4.1

Initial step

With cmds based localization methods, see Section 3.2, the sensors are localized given that the Euclidean distances

dij = ||si− sj|| (5.32)

are known [Shang et al., 2004]. With the novel approach we also want to incor-porate the knowledge that we measure the distances between each source and sensors, toa measurements. This will probably give more accurate sensors posi-tion estimates and sources posiposi-tion estimates as a bonus. The idea is to augment the distance matrixD with source to sensor distances

Dk aug=                 d_k1

D

which r K unique distance matrices, {Dk

aug}K1. By applying cmds, Algorithm 1,

on each extended distance matrixDk

augwill therefore give estimates of the sensor

positions ˆSN

k and source position ˆxk. So when all unique estimates are calculated

we get a total of K estimates of the sensor positions which are stacked in the set, {ˆSN

k }K1, and one source positions estimate, ˆXk.

In the nature of cmds these estimates are centered in arbitrary origins. By calcu-lating a sequential mean value of the sensor positions set {ˆSN

k}1K0 for each source

xk, will result in a mean value of the sensor positions defined as ¯ˆSN. For each k:th

sensor positions estimate transform it to fit the sequential mean, see Section 3.3 with a transformationT( · ) that includes both rotation/reflection and translation. Using a sequential mean separate origins of the K sensor positions estimates are avoided and possible influences of measurement noise is lowered.

Unfortunately the available sensors do not give inter-sensor distance measure-ments dij instead they must be approximated. One approach is to estimate the

distances utilizing the tdoa range measurements in (4.11) and by applying the triangle inequality, see Figure 5.6, we get an underestimated estimate ˆdij of the

distance dij. First we define the noise-free and zero bias range variants ˜rki and ˜r ij k

of the range measurements.

According to the triangle inequality we can state the following 0 ≤ ˜ri k− ˜r j k = ˜r ij k ≤ ||si − sj|| . (5.34)

Assume that infinitely many sources are spread evenly around the network lim

K→∞maxk∈K ˜r

ij

5.4 Sensor Localization 35 xk si sj r_ki _rj k dij

Figure 5.6: Illustration of the triangle inequality, ||x

k − si|| − ||xk − sj||

   ≤ ||si_{− s}j_||

where the limit approaches the true distance dij. So

ˆ

dij = max k∈K ˜r

ij

k (5.36)

is a plausible approximation of the true distance dij between sensorssiandsj.

In Section 4.3 it is stated the each target xk generates M range measurements:

rij_k andri

k. Instead of using one single range measurement the sample mean ¯r ij k

and ¯ri

k of the measurement vectors would lower the influence of the noise when

estimating the distances

ˆ dij = max k∈K ¯r ij k ˆ d_ki = ¯ri k (5.37) With the approximated distances the augmented distance matrices { ˆDk

aug}K1 can

be formed and the novel localization method can be applied, as described in Al-gorithm 3. Compared with the original cmds approach the novel approach have some advantages:

1) Lower the influence of measurement noise. 2) Higher precision of sensor position parameters. 3) Ability to estimate source position parameters.

Example 5.4: Sensor localization and cmds

Consider a sensor network with 8 sensors and 10 sources distributed in an area of 10 m2_{. The measurement noise is set to σ}

e= 0.1 m. With the novel approach the

sensors localization rmse is lower than with the original approach, see Figure 5.7. Note that sources are estimated with localization rmse 0.289 m.

36 5 Estimation

Algorithm 3 Novel approach to estimate 2-D positions of sensors and sources with cmds utilizing tdoa and toa measurements

for k = 1 → K do

Form the augmented distance matrixDk

aug

Calculate estimates SN

k , xkby applying cmds in Algorithm 1 withDkaug

if k=1 then . Initiate sequential mean with estimates from first source

¯ˆSN _{:= ˆS}N

1

Store source positions in vector ˆX1← ˆx1

else

Update sequential mean ¯ˆSN _{← ¯ˆS}N₊ 1

k+1T(ˆSNk ) − ¯ˆSN



Store transformed source positions in vector ˆXk ← T(ˆxk)

end if end for

5.4.2

Refinement step

The second step in the sensor localization algorithm is to minimize the difference between model output and measurements to achieve possibly better accuracy of the sensor positions. With the given range measurements and the sensor mod-els in Section 4.3 we can formulate an optimization problem that will minimize the residuals by finding the lse of the unknown parameters. With the novel approach on cmds localization above both estimates of the sensors and sources position are available. Therefore the toa sensor model will be considered, de-scribed in Section 4.3.1. To keep the structure compact for each target xk the

measurements and model outputs will stacked by sensors in vectors as Yk=h(xk;SN, br) +Vk, Cov(Vk) = Rk Yk =                r1 k r2 k .. . rN k                , h(xk;SN, br) = c ·                htoa(xk;s1, br) htoa(xk;s2, br) .. . htoa(xk;sN, br)                , Vk =                v1 k v2 k .. . vN k                (5.38)

where the measurement noise is wgn with known covariance, Rk. The covariance

matrix is a diagonal matrix with the estimated variances of the measurement vectorsri k in the diagonal Rk =            Var(r1_k) 0 0 0 . .. 0 0 0 Var(rN k )            . (5.39)

5.4 Sensor Localization 37 −10 −5 0 5 10 −5 0 5 x 1 [m] x 2 [m]

(a)Original cmds approach, Sensors rmse_{0.057 m} −10 −5 0 5 10 −5 0 5 x 1 [m] x 2 [m]

(b) Novel cmds approach, Sensors rmse_{0.045 m, Sources rmse 0.289 m}

Figure 5.7: Original and novel cmds approach. Ground truth are denoted with cross, sensor position estimates with circles and source position esti-mates with squares.

To integrate measurements from all sourcesXK we introduce the stacked

struc-ture by sources to drop the index k on measurement and noise vectors Y = H(XK;SN, br) +V , Cov(V) = R Y =                Y1 Y2 .. . YK                , H(XK;SN, br) =                h(x1;SN, br) h(x2;SN, br) .. . h(xK;SN, br)                , V =                v1 v2 .. . vK                (5.40)

where the covarianceR is a diagonal block matrix R =           R1 0 0 0 . .. 0 0 0 RK           . (5.41)

From now on the stacked structure by sources will be used to formulate the opti-mization problem.

From the structure in (5.40) we want to formulate optimization problem whose goal is to find the estimates of sensor positions, source positions and measure-ment bias that minimize the cost function

V (SN, Xk, br) = ₂1Y − H(XK;SN, br)

T

Y − H(XK;SN, br) = 1₂εTε (5.42)

in ls sense where ε = Y − H(XK;SN, br) is the residual vector. The lse of the

unknown parameters are then given by ‚ X_k, SN_{, b} r = arg min Xk,SN,br V (SN, Xk, br) = arg min Xk,SN,br 1 2εTε. (5.43)

38 5 Estimation

Algorithm 4 General Gauss-Newton

1: Given initial value ˆx(0), the functionh(x) and its gradient J(x) = −∂h_∂xT(x). Set

i = 0.

2: Set α(i)= 1. 3: Compute

ˆx(i+1) _{= ˆx}(i)_{+ α}(i)_J(x)JT_(x)−1

J(x)y − h(x)

4: If the cost V ( ˆx(i+1)) > V ( ˆx(i)), set α(i)= α(i)/2 and repeat from step 3.

5: Terminate if the change in cost, the change in estimate, or the size of the

gradient is small enough, or if the number of iterations has reached an upper limit.

6: Otherwise, set i := i + 1 and repeat from step 2.

The linear approach in Section 5.1 can not be applied because of the nonlinear relationships in the cost function. To find the lse we must introduce the Non-linear Least Squares(nls) framework. nls problems are usually solved by using iterative numerical optimization methods. Generally these methods updates the estimates iteratively as

ˆx(i+1)= ˆx(i)+ α(i)f(i) (5.44) where α(i) _{is the step length and f}(i) _{is the search direction in the i:th}

itera-tion. Such methods are characterized by their calculation of the search direcitera-tion. Here the basic Gauss-Newton algorithm will be used [Gustafsson, 2010, p. 49-50] which works as in Algorithm 4.

The Jacobian J(x) is central in Gauss-Newton and is constructed by stacking the gradients of all residuals ε on each other and defined as

J(Xk, SN, br) = ∂ε T ∂(Xk, SN, br)

, J(Xk, SN, br) ∈ R[2(K+N )+1]×MN K (5.45)

which can be calculated symbolically or numerically by approximating the deriva-tives using Finite-difference methods (fdm). Here it will be approximated by applying central differences as

J(i, :)≈H(θ + ei) − H(θ − ei).2 (5.46)

where the parameters are stored in the vector θ = [XKSN br]T for simplification, ei denotes the i:th column of the identity matrix, and  is the step size. The step

size should be small enough to find the linear region ofH(Xk, SN, br), but still

sufficiently large to avoid numerical ill-conditioning.

5.4.3

Localization Performance

The theoretical performance of the refinement step will be investigated by study-ing the sensor localization rmse usstudy-ing crlb theory. In Section 5.2 the FIM for an

5.4 Sensor Localization 39 unbiased estimator with Gaussian measurement errors is calculated as

I (θ) = JT_(θ)R−1_{J(θ) (5.47)}

for the sensor positions θ =SN_{. The crlb implies the following bound on the}

rmse rmse≥ q tr I−1 1:M  θ0 (5.48)

where θ0are the true parameters considering M independent measurements. Fig-ure 5.8 shows the crlb of the localization rmse of the refinement step evaluated with the benchmark configuration, see Section A.3. The localization rmse in-creases log-log linear with measurement noise variance. The rmse dein-creases as more independent measurements are generated from each source, which is intu-itive because more information are available.

40 5 Estimation 10−4 10−3 10−2 10−1 100 101 10−6 10−4 10−2 100 102 σ_e [m] RMSE [m]

(a)Localization rmse with respect to simu-lated measurement noise,M = 10

1 5 10 15 20 10−2 10−1 100 Measurements RMSE [m]

(b)Localization rmse with respect to gen-erated measurements from each source with simulated measurement noise,σe= 0.1 m.