Network-Based Positioning Using Last Visited Cells Report

(1)

Master of Science Thesis in Communication Systems

Department of Electrical Engineering, Linköping University, 2016

Network-based positioning

using Last Visited Cells

report

(2)

Tor Olofsson LiTH-ISY-EX–17/5040–SE Supervisor: Trinh van Chien

isy, Linköping University Yuxin Zhao

Ericsson Research, Ericsson AB

Examiner: Mikael Olofsson

isy, Linköping University

Communication Systems Department of Electrical Engineering

(3)

Abstract

The positioning performance with the lvc (Last Visited Cells) report is evaluated and compared with extended reports with signal strength data. The lvc report contains cell identities and time spent in the last cells listened to. This is an off-line data source and the purpose of the positioning is to extract information about users’ whereabouts, which for example can be used to optimize the cellular network or vehicular traffic. The positioning evaluation is done in Matlab with a log-distance model, a fingerprinting algorithm, and a new lvc specific algorithm. A particle filter and a particle smoother is used to process simulated lvc reports and extended reports with different amount of information. The results are com-pared and evaluated with regard to the positioning accuracy and the information density of the reports.

(4)

(5)

List of Figures

2.1 Tracking area . . . 6

2.2 Three directional antenna gains . . . 14

2.3 Trilateration . . . 15

2.4 Fingerprinting. . . 16

2.5 Closest fingerprint estimation. . . 19

2.6 Antenna gain difference . . . 23

2.7 Angle probability . . . 24

2.8 Voronoi diagram . . . 25

2.9 Voronoi algorithm vectors . . . 26

2.10 Voronoi probability . . . 28

2.11 Particle filter measurement update . . . 32

2.12 Particle filter time update . . . 33

2.13 Particle filter resampling . . . 34

2.14 Particle smoother transition density . . . 38

3.1 Process overview . . . 40

3.2 Map with trajectories . . . 42

3.3 Log-distance least squares fit . . . 44

4.1 Error cdf for the lvc report. . . 49

4.2 Comparison of the lvc report log-distance error for PF and PS . . 49

4.3 Error cdf for the rsrp report. . . 52

4.4 Comparison of the rsrp report fingerprinting error for PF and PS 52 4.5 Error cdf for the proximity report. . . 53

4.6 Comparison of cdf for the reports. . . 55

(8)

2.1 lvc report example . . . 7

3.1 Constant simulation parameters . . . 41

4.1 Position rmse and cdf for the lvc report . . . 48

4.2 Error statistics for the lvc report . . . 49

4.3 Position rmse and cdf for the rsrp report using the log-distance model . . . 50

4.4 Position rmse and cdf for the rsrp report using fingerprinting . . 51

4.5 Error statistics for the rsrp report . . . 51

4.6 Position rmse and cdf for the proximity report . . . 54

4.7 Error statistics for the proximity report . . . 55

(9)

List of Algorithms

2.1 Selecting fingerprint for position . . . 18

2.2 Fingerprint evaluation . . . 21 2.3 Voronoi values dP E, dSE . . . 27 2.4 Bayesian Filter . . . 30 2.5 Particle Filter . . . 35 2.6 Particle Smoother . . . 36 ix

(10)

(11)

Notation

Abbreviations

Notation Meaning

bs Base Station

cdf Cumulative Density Function

dod Direction of Departure

fbr Front to Back Ratio

ffbsi Forward Filter Backward Simulator

fp Fingerprinting

gps Global Positioning System

hpbw Half Power Beam Width

ldm Log-distance model

los Line of Sight

lte Long Term Evolution

lut Look-up Table

lvc Last Visited Cells

mu Measurement Update

nlos Non-Line of Sight

nan Not a Number

pdf Probability Density Function

pf Particle Filter

ps Particle Smoother

rmse Root Mean Square Error

rsrp Reference Signal Received Power

rss Received Signal Strength

ta Tracking Area

tdoa Time Difference of Arrival

tu Time Update

ue User Equipment

wcdma Wideband Code Division Multiple Access

(12)

Mathematical notation

Notation Meaning

Rnx Set of real numbers in n_xdimensions

p(y|x) Posterior probability for y given x

f (x) Dynamic model function

h(x) Measurement function

yk Measurement at time k in particle filter

eh Measurement error for h(x)

yi Measurement for cell i

ci Cell identity for cell i

xk System state at time k

xi_k State of particle i at time k

wi_k Weight of particle i at time k ˆxk Estimated state at time k

Gp Path loss

GA Antenna gain

φ Horizontal angle between site and a position

σ_RSRP2 Shadowing variance

dP E Distance from point P to edge E

dSE Distance from seed S to edge E ˆhi _{Fingerprint vector with index i} ˆhi

j Fingerprint value of ˆhifor cell j

dth Fingerprinting distance threshold

Pth Threshold rsrp for serving cell transition

yth Threshold rsrp for a threshold indicator

ymin Smallest observable rsrp

N Number of measurement cells

M Number of reference cells

Ns Number of reference sites

Nc Number of cells in the network

Nh Number of fingerprints

Ts Sample interval

v uevelocity

φN(x|σ2) Normal pdf with zero mean ΦN(x|σ2) Normal cdf with zero mean

kxk2 Euclidean distance or 2-norm of x

δ(x) Dirac delta function

O(x) Computational complexity x

(13)

1

Introduction

Information about the movement of users has potential to be used for many pur-poses. The obvious use for movement data gathered in cellular networks is opti-mizing the cellular network and enhancing the user experience. The 3GPP pro-motes research contributing to the cause of minimization of drive tests. A drive test is an expensive and limited way of analyzing the network. If the results from a drive test can be achieved with other means, considerable costs for network operators can be avoided.

Positioning by the network can be done in real time or in off-line mode. Real time positioning is important for locating people in need during emergency calls, and in USA it is an announced goal from the Federal Communications Commis-sion to increase accuracy and capabilities of handsets [1]. Off-line positioning is of course not helpful in emergencies, but can instead be valuable for large data analyses.

Analyses on data from positioning can for example be used for traffic plan-ning to reduce the time spent waiting in queues and to avoid traffic jams. Less time in traffic means less pollutants which in the end would prevent many cases of disease, apart from making the saved time available for productive work. It is no surprise that companies seek to use accessible information to produce addi-tional value.

For online cell phones, there are many sophisticated methods available for de-termining position. One example is the uplink time difference of arrival (TDOA). The case with off-line cell phones and information flowing in a single direction has more limited options.

The Last Visited Cells report is a way for cell phones to send stored off-line positioning data to the network. When a cell phone is not used for calls or inter-net connectivity, the inter-network knows only roughly in what area the phone is [2]. While the phone idly moves, it saves a list of the network cells it listens to and

(14)

how long it listened to each cell. This information is saved to the Last Visited Cells list which then can be requested from the phone when it actively connects to the network. The Last Visited Cells report and its use in off-line positioning will be the focus of this thesis.

1.1 Positioning methods

Positioning in wireless networks is an active research field. Signal strength based positioning with a path loss model is one of the options investigated [3]. The posi-tion accuracy indoors using a path loss model has been reported to be mostly the same when using measured signal strength values and when applying a thresh-old to them to achieve a binary proximity report indicating which network nodes are strong [4].

Along with the increasing computer power available, the particle filter has become a popular solution to non-linear filtering problems. The particle filter achieves good results in many tracking applications, and enables map matching using large databases [5]. The database approach is commonly called fingerprint-ing, and has been reported useful for positioning handovers to estimate travel time along with crude position [6]. Accuracy of positioning using signal strengths can be significantly improved with fingerprinting compared to using a simple path loss model [7].

Earlier research at Ericsson has investigated the performance when combin-ing a Direction of Departure algorithm with the Time Advance command to get a position estimate in cellular networks [8]. While the Direction of Departure approach is possible to utilize in the models in this thesis, the Time Advance measurement needs to be substituted with another measure of distance that is not dependent on bidirectional communication with the network.

1.2 Aim

The aim of this thesis is to evaluate the positioning ability with the Last Visited Cells report and investigate what positioning improvements can be made by us-ing reports with additional signal strength information. This will be done by answering the following questions.

• What positioning performance can be achieved with the lvc report? • How is the positioning improved with additional signal strength

informa-tion?

1.3 Scope

The Last Visited Cells report is defined in a technical specification of LTE release 13, which will be the reference network standard in this thesis. Evaluations of

(15)

1.4 Outline of the report 3 algorithms will be done in Matlab with generic computer hardware. Models will be based on a simple log-distance model and a fingerprinting algorithm.

1.4 Outline of the report

The outline for the remainder of the report is as follows.

• Chapter 2 starts with necessary definitions for cellular networks and contin-ues with relevant theory. The measurement reports which are investigated are defined with a probabilistic framework. Algorithms for deriving posi-tion probabilities from the reports are then presented and explained. • Chapter 3 describes the method in which the theory and the algorithms

de-fined in Chapter 2 are used to simulate measurement reports and evaluate the positioning performance of the reports.

• Chapter 4 presents the results of the evaluation achieved from following the method in Chapter 3. The results are displayed in figures and tables. • Chapter 5 discusses the results and the methods used.

(16)

(17)

2

Theory

The methods used in this thesis have a theoretical background and make certain assumptions and simplifications of the environment. This section serves as a basis for the theory behind the methods.

Models in this section present probability functions p(y|x) which are used in a particle filter. The particle filter uses the probabilities generated by the models to estimate the positions x of a simulated cell phone.

2.1 Cellular Networks

Following sections in this master thesis requires some knowledge about cellular networks. This section defines some of the crucial concepts.

For the average person, what may come first to mind when thinking of wire-less networks is the User Equipment (ue), which usually is a cell phone. The ue acts as a client and is served by the network. While the network is static, the ue is dynamic and moves around and is only spuriously connected to the network.

The network connection of a ue is granted by a Base Station (bs), which is located on a site. Put simply, each bs has one or more antennas, and each antenna provides connectivity to a cell. A cell comprises the area where ues are served by the same antenna. The ue estimates the signal strength in a cell with the Reference Signal Received Power (rsrp).

In LTE, multiple cells are grouped together into a larger Tracking Area (ta) [9]. When a ue is idly loitering, the network is only interested in knowing its ta and does not need to know the specific cell the ue is listening to. In this case, the uestays silent while moving between cells belonging to the same ta, and only reports to the network when moving to a new ta. This way, the ta is used to reduce the need for the ue to report its serving cell to the network.

When the network needs to contact the ue, it sends a paging message from 5

(18)

Figure 2.1: Tracking area. Two tracking areas are shown with a trajectory running through them. The trajectory crosses the boundaries of the hexago-nal cells multiple times, but only crosses a tracking area boundary once. The tracking area crossing is indicated by an arrow.

multiple cells in the ta to be sure that the ue will receive the message. The reduction of upstream data and thus lower power consumption at the ue is paid for by extra downstream capacity use at the bs. An illustration of cells belonging to different tas is shown in Figure 2.1.

2.2 Last Visited Cells

When the ue moves to a new ta, it may report a list of its Last Visited Cells (lvc) to the network [10]. The lvc report contains entries of cell identities and the time spent in cells while moving through the ta. There is a maximum of 16 entries that can be stored in the list. If a new entry is added when the list is full, the oldest entry is removed from the list. An example of a lvc report is shown in Table 2.1 on the next page.

The lvc report can be interpreted as a list of transitions between cells. For example, if the cell visited before has ID 5 and the next cell has ID 4, it can be interpreted as a transition between cells 5 → 4. Each entry in the lvc report thus represents a transition between two cells in two consecutive entries. The transi-tions corresponding to a report is shown in the rightmost column in Table 2.1 on the facing page. This cell transition interpretation is used in this thesis to estimate positions for the lvc report.

To be able to receive calls, a ue always listens to at least one cell while waiting to be paged by the network. This serving cell being listened to is assumed to be the strongest cell with a signal to the ue, or at least not much weaker than the strongest cell.

(19)

2.3 Reports 7 Table 2.1: lvc report example. An example of what an lvc report could look like is displayed. Each line contains the time spent and what cell was visited during this time. The entries can also be interpreted as transitions between two cells. The same cell can be visited multiple times and therefore occur more than once. The list also does not need to be exactly 16 entries long.

Index Time spent Cell ID Transition

1 109 5 — 2 244 4 5 → 4 3 111 3 4 → 3 4 267 5 3 → 5 5 250 7 5 → 7 6 94 8 7 → 8 7 268 9 8 → 9 8 239 7 9 → 7 9 210 3 7 → 3

cell at all times, the ue selects the strongest cell as its new serving cell when it is perceived as stronger than the current serving cell. When this event with a new serving cell occurs, the ue adds a cell transition to its lvc list.

To reduce the amount of cell transitions, a threshold for minimum signal strength difference to the current serving cell is used [11]. This threshold acts as a hysteresis and avoids filling the list with numerous transitions between two neighbouring cells with approximately the same signal strength. To further re-duce the amount of transitions, there is also a time-to-trigger which defines the minimum time for which the threshold condition must hold before a transition actually is done.

2.3 Reports

In this thesis, three different types of reports for serving cell transitions are used. Each time entry in a report is described by a measurement vector y which is used as input to a filter. The Last Visited Cells report is the first report. Two other reports using rsrp values are also defined. These two reports are designed so they might contain more and better information than what is available in the limited lvc report.

The measurement vector for the lvc report is defined as

y= [cnew cold c1 · · · cM], (2.1)

where cnewand coldidentifies the old and the new serving cell for a transition, and

ci identifies reference cells which should be weaker than cnewfor i = 1, . . . , M. The second type of report is the rsrp report. This report describes the N strongest cells with cell identities and an rsrp value for each cell. This variant

(20)

with explicit rsrp values for cells is used in [7, 12, 4]. The measurement vector is defined as

y= [c₁ y1 c2 y2 · · · cN yN cN +1 · · · cN +M], (2.2) where ciis the cell identity and yiis the rsrp for cell i in the report. Additionally,

M cells which should be weaker than yN are included.

The third report is the proximity report. Instead of using explicit values for the measured rsrp of cells, the report contains the N cells which are stronger than or equal to a threshold yth, and M cells which are weaker than the thresh-old. This report is similar to a report in [4], which used a static threshthresh-old. A disadvantage with a static threshold when using the lvc sample times is that occasionally no cells at all are above the threshold, which provides very little accuracy for position prediction. The measurement vector is defined as

y= [c1 c2 · · · cN yth cN +1 · · · cN +M], (2.3) where ci as before is the cell identity. The vector is defined so that the cell cN defines the threshold yth, and therefore is equal to the threshold while N − 1 cells in the list are stronger than the threshold.

While the three reports are described by different vectors, they are all as-sumed to be sampled at the time of a cell transition as defined by the lvc report. The two rsrp reports are therefore regarded as extensions to the lvc report. They both contain cnewas their strongest cell, but cold is not necessarily included.

2.4 Measurement models

The reports are used to estimate positions by passing their measurement vectors y through a filter. In the filter, a probability function p(y|x) is used to extract information from the vectors. This section presents what this function looks like for the different reports.

The probabilities for explicit rsrp measurements will be explored first, and then an adaptation for the cell transitions in the lvc report will be made. A general definition of the probability p(y|x) will be made for the reports, and then this will be used to define more precise specifications for a log-distance method and a fingerprinting method which are described later.

Assume that a measurement y from a ue can be stated as

y = h(x) + eh, (2.4)

where h(x) is a deterministic function giving a predicted rsrp value for a position

x, and the term eh is the error of the measurement y because of noise and an imperfect prediction function. The definition in Equation (2.4) can be rewritten as

eh= y − h(x), (2.5)

which is an equation to calculate the error ehfor a measurement y given a position

x. Using a known distribution for the error eh, the equality Equation (2.5) is used to infer probabilities for measurements.

(21)

2.4 Measurement models 9

2.4.1 RSRP and proximity reports

For the measurements in the rsrp report, y describes a numeric rsrp value. The probability p(y|x) for a measurement y given a position x is then defined using Equation (2.5) on the preceding page as

p(y|x) = p(eh= y − h(x)), (2.6)

which is evaluated with the distribution for eh.

Both the rsrp report and the proximity report expects the rsrp of some cells to be on a certain side of a threshold. For the cells ci in the vector y which are re-lated to a threshold and do not have a numeric rsrp value, an implicit threshold indicator yi is used for the probability calculation. Using the notation in Equa-tion (2.4) on the facing page, an rsrp measurement from a cell is converted to a threshold indicator as y =        0, h(x) + eh< yth 1, h(x) + eh≥ yth , (2.7)

where ehis the error of the actual rsrp measurement, and ythdenotes the rsrp value used as the threshold. A threshold indicator y = 0 infers the inequality

h(x) + eh< yth, (2.8) which is rewritten as

eh< yth− h(x) (2.9) to get an inequality which is used to get the probability p(y|x) as

p(y|x) =        p(eh< yth− h(x)), y = 0 1 − p(eh< yth− h(x)), y = 1 . (2.10)

Now there is a definition for the probability p(y|x) for a single measurement y given a position x in both the case when y is a numeric rsrp value and a threshold indicator. The probability for a measurement vector y in the rsrp report or the proximity report can then be stated as

p(y|x) = N +M

Y

i=1

p(yi|x), (2.11)

where the proper definition of p(yi|x) is used depending on if yiis a numeric rsrp value or a threshold indicator.

If the distribution for ehis assumed Gaussian, the evaluation of the probabil-ity functions in Equation (2.6) and Equation (2.10) can be stated more explicitly. A Gaussian distribution is described by the normal pdf and cdf. For a distribu-tion with a mean of zero, the normal pdf is defined by

φN(x|σ2) = 1

σ√2πe

−x2

(22)

and the normal cdf is defined by Φ_N(x|σ2) = 1 σ√2π x Z −∞ e−2σ2x2 (2.13)

where x is the deviation from the mean, and σ2_{is the variance of the distribution.}

For the measurements y, the error ehacts as the deviation from the mean. Using the normal pdf and cdf, the probability function p(y|x) for numeric rsrp values is written as

p(y|x) = φN(y − h(x)|σe2h), (2.14)

and the probability for a threshold indicator is written as

p(y|x) =        Φ_N(y_th− h(x)|σ_e2 h), y = 0 1 − ΦN(yth− h(x)|σe2h), y = 1 , (2.15) where σ2

ehis the variance of the error eh. These two definitions in Equation (2.14)

and Equation (2.15) are used to calculate the probabilities p(y|x) for measure-ments in the rsrp report and the proximity report, given a prediction function

h(x). The actual definitions of h(x) which are used to evaluate the report

per-formance are specified later in the log-distance method and the fingerprinting method.

With the probability definitions specified for both numeric rsrp values and threshold indicators, the total probabilities for the measurement vectors in Equa-tion (2.11) on the previous page can be stated in expanded form for the rsrp report and the proximity report. Recall that the rsrp report is represented by a vector containing rsrp values and cell identities for N cells, and only cell iden-tities for M cells which are weaker than the N:th cell. The proximity report is represented by a vector containing M cell identities and a threshold yth which separates all but one of the cells into two groups on either side of the thresh-old. The probability p(y|x) for the rsrp vector can be stated using the equality in Equation (2.14) and the inequality in Equation (2.15) as

p(y|x) = N Y i=1 φN(yi = hi(x)|σe2_i) N +M Y i=N +1 Φ_N(y_N− h_i(x)|σ_e2 i), (2.16)

and the probability for the proximity report is stated as

p(y|x) =        N_Y−1 i=1 1 − ΦN(yth− hi(x)|σe2_i)        φN(yth= hN(x)|σe2_N) N +M Y i=N +1 Φ_N(y_th− h_i(x)|σ_e2 i). (2.17)

2.4.2 LVC report

While the rsrp report and the proximity report uses a single measurement of rsrp to create y, the lvc report is assumed to use two measurements for the

(23)

2.4 Measurement models 11 equality

ynew− yold = Pth, (2.18)

which describes the cell transition condition where Pthis the rsrp threshold for a transition. Because the lvc report only contains two cell identities cnew and

cold for the new and the old serving cell, the measurements ynewand yoldcan not directly be used. By substituting ynew and yold using Equation (2.4) on page 8, the equation becomes

hnew(x) + enew− hold(x) − eold = Pth, (2.19) which is rewritten with the error terms on the left-hand side as

enew− eold= Pth− hnew(x) + hold(x), (2.20) where the error terms are used for the probability function

p(cnew, cold|x) = p(enew− eold= 0), (2.21) which determines the probability that the equality in Equation (2.18) is valid.

The lvc report is also used to assure that the cell cnewis the strongest cell for position x. This can for another cell yibe stated as

yi < ynew, (2.22)

which by substitution using Equation (2.4) on page 8 becomes

ei− enew< hnew(x) − hi(x) (2.23) and results in the probability function

p(cnew, ci|x) = p(ei− enew< hnew(x) − hi(x)). (2.24) The probabilities for Equation (2.21) and Equation (2.24) can be stated with the normal pdf and cdf if assuming a Gaussian distribution for the errors eh. Both these equations contains two separate error terms with one distribution each. As the sum of two Gaussian variables is also a Gaussian variable, with a combined variance equal to the sum of the two variances, the probabilities can be interpreted as having a single error using a variance with two parts.

The probability for Equation (2.21) for the transition condition can then be stated as

p(cnew, cold|x) = φN(enew− eold|σnew2 + σold2 ), (2.25) which by using Equation (2.20) for enew− eoldbecomes

(24)

Likewise, the probability for Equation (2.24) on the previous page which com-pares weaker cells to the new serving cell is stated as

p(cnew, ci|x) = ΦN(hnew(x) − hi(x)|σnew2 + σi2). (2.27) The total probability p(y|x) for the measurement vector representing the lvc report can now be stated using these two equations. The vector contains the cell identities cnew and cold for the new and old serving cell, as well as M reference cells which should be weaker than cnew. The probability for the lvc report is stated as

p(y|x) = φN(Pth− hnew(x) + hold(x)|σnew2 + σold2 ) N Y

i=1

Φ_N(h_new(x) − h_i(x)|σ_new2 + σ_i2), (2.28) where Pthas before is the transition threshold.

2.5 Channel propagation models

Line of Sight (los) and Non-Line of Sight (nlos) are used to describe propaga-tion situapropaga-tions in radio communicapropaga-tion. los indicates that there are no obstacles in the direct path between the transmitting end and the receiving end, while nlosindicates that there are obstacles blocking the signal and that the strongest signal path is not necessarily the shortest.

Even if there is los, the signal quality can still be heavily compromised by reflections on surfaces and diffraction effects. Measurements relying on the signal strength and signal propagation time are generally good in los scenarios, while they can be utterly unreliable in a nlos case. To increase robustness of a model against nlos, two separate error distributions for los and nlos can be combined. An example of a two-mode Gaussian mixture is

pE = αN (0, σLOS2 ) + (1 − α)N (0, σN LOS2 ), (2.29) which can be configured for different scenarios by varying the set of parameters [3]. The parameters σ_LOS2 , σ_{N LOS}2 and α in Equation (2.29) describe the variances for the los and nlos distributions and their weights when combined.

2.6 Received Signal Strength

A ue continuously measures the Received Signal Strength (rss) of nearby cells to be able to request a hand-over if the signal of current serving cell becomes too weak. The rss is estimated by calculating the power of specific reference signals, which yield Reference Signal Received Power (rsrp) which will be used as rss in this thesis. The rsrp is dependent on the transmit power Pbs of the base station,

which is affected by losses in the feeder cable, Lc, and antenna gain, GA. During the propagation of the signal from the bs to the ue, a channel and path dependent

(25)

2.6 Received Signal Strength 13 gain will be present, here notated as Gp(d). Using these notations and expressing the gain in dB, the received signal strength Pr at the ue can be written as a sum of many separate terms

Pr = Pbs+ Lc+ GA+ Gp dB. (2.30)

The path loss gain Gpin this case comprises attenuation due to signal spread, multipath interference and shadow fading.

2.6.1 Antenna Gain

The typical site has three directional antennas which each cover 120◦_{. The gain}

of a directional antenna can be expressed as a function of polar coordinates,

GA(φ, θ), where φ is the horizontal angle and θ is the elevation angle relative to the horizontal plane [8]. To simplify, the antenna gain is often considered to be separable into a horizontal and vertical gain component

GA(φ, θ) = gAh(φ) + gAv(θ) dB. (2.31)

where gAh(φ) and gAv(θ) are the horizontal and vertical gain respectively. Information about the antenna gain is typically available as parameters in the antenna data sheet. With a few parameters describing the antenna characteristics, a standardized model can be used to calculate the approximate gain of multiple antennas. Another possibility is to use a look-up table (lut) of the antenna gain values for φ and θ.

In this thesis, the elevation angle θ will be neglected. Only the horizontal gain gAh(φ) in Equation (2.31) will be used, and use the notation GA(φ) for this partial antenna gain. Further, the horizontal angle φ will be referred to as the Direction of Departure (dod).

An example of the horizontal antenna gains of a site with three antennas is displayed in Figure 2.2 on the following page. As can be seen in the figure, each antenna dominates its own interval of angles.

2.6.2 Path loss

To be able to use rsrp measurements for estimating a distance, a path loss model is needed. The path loss can potentially be a very complicated function of the sig-nal propagation environment. In the case that there is little information available about the propagation environment, a simple path loss model is needed.

A common model for propagation losses of radio signals is the log-distance model

Gp(d) = 10B log(d

d0) + C dB. (2.32)

The model in Equation (2.32) uses a constant C which corresponds to the path loss at a reference distance d0. The path loss is then determined by the

(26)

0 50 100 150 200 250 300 350 −45 −40 −35 −30 −25 −20 −15 −10 −5 0 Angle [degree] An tenna gain [dB] Antenna 1 Antenna 2 Antenna 3

Figure 2.2: Three directional antenna gains. The antenna gains for three antennas sharing the same site are shown. Their maximum gains are located on three different angles separated by 120◦_{to optimize coverage.}

distance d and the path loss exponent B, which is dependent on the propagation environment and usually resides in the interval [2, 4].

Path loss prediction is a method used in trilateration when a set of rsrp mea-surements are used to estimate a position. An example of trilateration is dis-played in Figure 2.3 on the facing page.

2.7 Log-distance model

The path loss equation in Equation (2.32) on the previous page is able to make a prediction of the path loss component of the rsrp in Equation (2.30) on the preceding page. Only the path loss Gpand the antenna gain GAare dependent on the distance or the angle, which together describes a position relative to a bs. The other terms of Pr in Equation (2.30) on the previous page are in that sense constant and can be merged into a single constant term. In the following model, the constant is defined as the expected rsrp at distance d0and is denoted A0. By

adding terms for the expected rsrp and the antenna gain GA(φ) for the direction

φ, a complete model for calculating the expected rsrp for an arbitrary position

is acquired as

Pr(d, φ) = A0+ 10B log(d

(27)

2.8 Fingerprinting 15

Figure 2.3: Trilateration.The rsrp from a cell can with a path loss function be translated to a distance. In this figure, three circles are created from rsrp values. A position for the measurements can then be estimated by finding the intersection of the three circles, here indicated by an arrow.

where d and φ describes a relative position.

Equation (2.33) on the facing page will be referred to as the log-distance model. It will be used as a prediction function h(x) introduced in Section 2.3 on page 7 to use for calculating the probabilities p(y|x) for measurement vectors y given a position x. This model will be used to make predictions for the rsrp from cells, using one distinct function for each cell to be able to represent different site positions and values for the parameters A0and B.

2.8 Fingerprinting

An alternative to using an explicit model to predict the rsrp of cells to compare to the reported measurements is using a method called location fingerprinting [7, 13]. Fingerprinting is an algorithm which uses pattern recognition to identify how well new measurements match a database of so called fingerprints from pre-viously collected measurements. While the log-distance model in Equation (2.33) on the facing page does not perform well for significant shadowing variances and nlos conditions, fingerprinting has been shown to perform better for these cases [14]. Using fingerprints, an inverse function can be emulated to estimate a position from measurements of rsrp values. An example of a function being represented by fingerprints is displayed in Figure 2.4 on the next page.

This simple form of machine learning requires a first off-line phase of col-lecting a training data set to create a fingerprint database. In this off-line training phase, an area is selected and fingerprint vectors for positions in the area are sam-pled. No computations are made with the training data set until the fingerprints are compared to measurements in an on-line positioning phase.

(28)

50 100 150 200 250 300 350 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 x [ ◦ ] sin(x)

(a)Continuous function y = sin x.

50 100 150 200 250 300 350 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 x [ ◦ ] sin(x)

(b)Map of values from x to y.

Figure 2.4: Fingerprinting. A continuous function is shown together with a fingerprinting representation of the same function. Fingerprinting creates a map which links samples from the setx to samples in the set y. Note that

the fingerprint map only represents a finite set of mapped value pairs and its precision is limited by the density of samples.

While the positions of sites and transmitting power of cells are known to the network operator and are used in the log-distance model, the fingerprinting method does not need such information. Because of this independence of infor-mation otherwise needed in an explicit model, a fingerprinting database can be created by virtually anyone and not only the network operator.

One example of a crowd sourced database used for location fingerprinting is the Mozilla Location Service [15]. This open service relies on individuals collect-ing data as well as organizations sharcollect-ing information about infrastructure, and allows users to get a position from observations of various networks like cellular networks and Wifi networks.

2.8.1 Training

The fingerprints describe features of the training data, which in this thesis will be the rsrp of cells in the network. Each fingerprint vector describing these features is also associated with a position. Both rsrp and position are continuous variables, and the goal of the training phase is to create a database which links positions to vectors of rsrp values.

There are different approaches to creating this database. One alternative is that the provider of the service collects accurate estimates for all possible posi-tions and under different weather condiposi-tions. This might be expensive and time consuming. Another alternative is to let users with GPS receivers contribute with samples to expand the database, as exemplified with the Mozilla Location Ser-vice. However, with users having different measuring equipment, the data might be less accurate. A third option is to let a simulation framework predict samples

(29)

2.8 Fingerprinting 17 using a model.

The last option is used in this thesis. The same simulation model will be used to generate the fingerprint database and the measurements in the model performance evaluation. This will in a sense be similar to the case where the real world is used to create both the database and evaluation measurements.

A square grid of equidistant points will be assumed for the reference positions of the fingerprints which describe the rsrp features of nearby cells. Each grid point is then associated with a fingerprint vector of rsrp values. The distance between the grid points will be one of the limits for how well a position can be estimated.

As each fingerprint in practice will represent a square area and not just a specific position, the fingerprint should be created from an average of multiple measurements in the area. Additionally, the fingerprints should also be averaged from multiple measurements over time to avoid large influence of temporary fad-ing effects.

Assume like before in Equation (2.4) on page 8 that the observed rsrp y from a cell depend on the position x as

y = h(x) + eh (2.34)

where h(x) represents a deterministic path loss and ehrepresents noise from shad-owing.

An average measurement ˆh(x) is then defined as ˆh =        E(y), y ≥ y_min NaN, y < ymin , (2.35)

where E(y) is the expectancy of y, which here is interpreted as the average of y in the local area around the position x over a period of time. To model the case when a signal is too weak to be measured by a non-ideal receiver, a threshold ymin is used to set a lower bound for the fingerprint measurements. The fingerprint ˆh is then set to NaN, not a number, as an indicator if its rsrp is below the threshold

ymin.

The fingerprint vector ˆh(x) is defined as the collection of averaged measure-ments

ˆh(x) = [ˆh1(x) ˆh2(x) · · · ˆhNc(x)], (2.36)

where ˆhj(x) is the averaged measurement for cell j and NCis the number of cells in the network.

When using a grid with a finite set of points with fingerprints, each finger-print position can instead be indexed with a number. Let pi _{be the position for} the fingerprint at the grid point indexed by i. The fingerprint vector in Equa-tion (2.36) can then be stated as ˆhi _{for grid point i with the elements}

ˆhi _{= [ˆh}i

1 ˆhi2 · · · ˆhiNc], (2.37)

where ˆhi

(30)

2.8.2 Positioning

When the training phase is completed, new measurements of rsrp values can be compared to the fingerprint vectors. The fingerprints will like the log-distance model be used to provide a prediction function h(x) as defined in Section 2.3 on page 7. Unlike the log-distance model, the fingerprints only covers a discrete set of positions x, and also sometimes does not always have a numeric value for each cell. This section describes how these issues are handled.

Because the location fingerprinting method works with a finite set of reference positions and reference vectors of rsrp values, some kind of conversion has to take place to translate a continuous position x to a fingerprint. This allows the use of the probability p(y|x) for a vector y as defined in Section 2.4 on page 8. In the solution chosen here, any position x will be mapped to the closest reference position pi _{and use its fingerprint vector ˆh}i_.

To restrict the possible positions to the training area, a position will be as-signed zero probability if it its distance to the closest reference position is more than dth. This mapping to the closest reference point can be compared to using a set of step functions around the reference positions to generate a partially con-tinuous function. It also has the effect that the integration of the probability over

x∈ R2does not equal unity, which luckily is ignored in a particle filter.

An example of using the closest fingerprint like this is illustrated in Figure 2.5 on the facing page. The algorithm to select a fingerprint vector based on a posi-tion x is described in Algorithm 2.1.

Algorithm 2.1:Selecting fingerprint for position

Position estimation: Using fingerprints ˆhi on positions pi_{, the posterior} probability p(y|x) for a measurement vector y given a position x is determined as following.

1. Calculate the distance of the position x to each of the grid points as

di = x− p i ₂ (2.38)

2. Find the grid point closest to x as ˆi = arg min

i

di (2.39)

3. Calculate the probability p(y|x) as

p(y|x) =       

p(y| ˆhî), d_î≤ dth 0, d_î> dth

(2.40) where dthis a threshold which limits the positioning to the training area. Having defined the translation from continuous positions to the finite set of reference positions in the fingerprinting database, only reference positions will

(31)

2.8 Fingerprinting 19 50 100 150 200 250 300 350 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 x [ ◦ ] sin(x)

(a)Fingerprinted function.

50 100 150 200 250 300 350 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 x [ ◦ ] sin(x)

(b)Partially continuous function.

Figure 2.5: Closest fingerprint estimation.A partially continuous function is created from the fingerprint map by assigning each value of the continuous variablex to the fingerprint of the closest sampled x. The resulting function

is constructed by a sequence of step functions.

be considered in the algorithm for comparing fingerprints. The probability p(y|x) will be calculated as in Section 2.4 on page 8, with the fingerprint hi

jas the value of the prediction function h(x), where i is the index closest to the position x and

j is the cell that y represents. This is correct for when hi_jhas a numeric value. For the case when hi

j is NaN and only known to be below the threshold ymin, a new definition for p(y|x) is needed.

When the fingerprint value is below the threshold ymin, so that ˆhij = NaN, the probability for a measurement y is given by

p(y|ˆhi

j) = p(y|ˆhij = NaN), (2.41)

which by using Bayes’ rule can be rewritten as

p(y|ˆhi j) =

p(ˆhi

j = NaN|y)p(y)

p(ˆhi_j = NaN) , (2.42)

which depends on the prior probability p(y) and the probability p(ˆhi

j = NaN), which both are unknown and hard to find. In [7], another approach is suggested to produce an appropriate likelihood for this case which is possible to calculate using only known variables and has provided satisfactory results.

The value of p(y|ˆhi

j) when ˆhij= NaN is defined so that it satisfies the inequality

p(y|ˆhi

j) ≤ p(y|ˆh i2

j ), (2.43)

where the fingerprint value on the right-hand side ˆhi2

j , NaN is above the thresh-old, while ˆhi

(32)

locations where the difference between y and ˆhi

j could be large, while preferring locations where the difference is known to be limited, which in this case is ˆhi2

j . The proposed method to achieve this for y with a numeric rsrp value is

p(y|ˆhi

j) = min_m_∈M

j

p(y|ˆhm

j ), (2.44)

where the set Mjis defined as

Mj =ni|ˆhij , NaNo . (2.45)

This method assures that the probability in Equation (2.44) is valid for the in-equality in Equation (2.43).

For a threshold indicator y meaning h(x) + eh < yth, the value for h(x) can simply be set to yminso that

p(y|ˆhi

j) = p(y|h(x) = ymin). (2.46)

While a numeric yj does not produce a different probability depending on the sign of its difference to ˆhi

j, the threshold indicator does. This means that it is not unreasonable to use yminfor the threshold indicator, while it is for the numeric measurement y which may yield the same probability for a fingerprint ˆhi

jwith a very high value and a very low value.

The algorithm described above is summarized in Algorithm 2.2 on the next page.

2.9 DOD Voronoi model

In this section, a new method for estimating positions from the lvc report is presented. The method uses two separate models, a Direction of Departure (dod) model to restrict positions with regard to angles and a nearest neighbour (Voro-noi) model to limit distances to sites. Only the information in y = [ynew yold] about the new and the old serving cell is used for signal strength assumptions.

The combination of the dod model and the Voronoi model is combined to a new measurement model used for the lvc report in addition to the log-distance model and fingerprinting. This combined dod and Voronoi model is simply defined by the product of their two individual probabilities, here denoted by

pDOD(y|φ) and pV(y|x). Using the position x to calculate the angle φ between x and the site, the dod and Voronoi model probability function is stated as

p(y|x) = pDOD(y|φ)pV(y|x). (2.51)

2.9.1 Angle probability

With measurements from two cells on the same site, an approximation of the Direction of Departure (dod) can be made. The following method uses the differ-ence in the antenna directions for the cells and consequently the antenna gains to derive an angle from rsrp measurements.

(33)

2.9 DOD Voronoi model 21

Algorithm 2.2:Fingerprint evaluation

Training: Collect a training data set of fingerprint vectors ˆhifor i = 1, . . . , Nhdefined by their position pi and their rsrp values ˆhij for

j = 1, . . . , Nc, where Ncis the number of cells in the reference network and

Nhis the number of fingerprint vectors. Let NaN be the value of ˆhijwhen the observed rsrp is below the threshold ymin.

Select fingerprint: Select a fingerprint ˆhi for position x using Algorithm 2.1

on page 18.

Probability evaluation: The probabilityp(y|x) for a measurement vector y given a fingerprint ˆhi _{is determined as in Section 2.4 on page 8 with the} probability p(y|x) defined as following, with ˆhi

jas the fingerprint for the cell j that y represents.

• If ˆhi

jis a numeric value above the threshold ymin, use

p(y|x) = p(y|ˆhi

j). (2.47)

• If ˆhi

jis NaN, below ymin, use

p(y|x) = p(y|h(x) = ymin), (2.48)

when y is a threshold indicator, and

p(y|ˆhi

j) = min_m_∈M

j

p(y|ˆhm

j), (2.49)

when y is a numeric rsrp and the set Mj is defined as

(34)

The two measurements from two cells i and j can each be written in the form of Equation (2.30) on page 13, giving

Pr,i = Pbs,i+ Lc,i+ GA,i+ Gp,i dB, (2.52)

Pr,j= Pbs,j+ Lc,j+ GA,j+ Gp,j dB. (2.53)

With two antennas sharing virtually the same position, the complicated and typically inaccessible path loss Gpis assumed to be the same for both rsrp mea-surements at the ue. This assumption simplifies the problem by disregarding scenarios where a different nlos paths may be stronger than the los path. In these cases, the actual rsrp measurements by the ue could represent different signal paths with different departure angles, and invalidate the assumption of equal Gp. Still, it is the best assumption possible to make without knowledge of the environment in which the signal propagates.

By assuming equal path loss Gp from the two cell antennas, the difference between Equation (2.52) and Equation (2.53) can be put as

Pr,i− Pr,j = GA,i− GA,j+ Pbs,i− Pbs,j+ Lc,i− Lc,j dB. (2.54)

The transmit power Pbs and the feeder loss Lcare easily available for a bs, and

their differences between cells i and j might from now on be assumed to cancel out to zero. Remaining unknown in the equation, is the difference of the antenna gains

Hij(φ) = GA,i(φ) − GA,j(φ) dB. (2.55)

The probability for an angle φ can be derived from the difference of antenna gain. With information about the direction and the directional gain of each an-tenna, a relative antenna gain function can be constructed as the difference of antenna gain for each φ. A diagram of the relative antenna gain Hij for two of the antennas in Figure 2.2 on page 14 is illustrated in Figure 2.6 on the facing page.

The diagram in Figure 2.6 on the next page shows that the dod can be ex-tracted from where the relative gain intersects the measured rsrp difference. As the relative gain function is continuous, there will be an even number of intersec-tions for any rsrp difference. Two intersecintersec-tions will be found in the normal case. However, it is apparent that one angle is unlikely as it is in the territory of the third antenna, which was not strong enough to be one of the two measured. By comparing the antenna gains of the two measured antennas with the third, it is easy to select the correct angle.

The vector y = [ynew yold] contains information about the cell transition. This information indicates whether the transition was between cells on the same site or if the cells belonged to different sites.

Let GA,i be the expected perceived antenna gain for cell i, and N be the number of cells on the site, where the cells are ordered by decreasing signal strength so that the strongest cell is first. The strongest cell is then the new serving cell ynew which is indexed by 1. For a transition between cells on the

(35)

2.9 DOD Voronoi model 23 0 50 100 150 200 250 300 350 −40 −30 −20 −10 0 10 20 30 40 Angle [degree] An tenna gain [dB] Antenna 1-2

Figure 2.6: Antenna gain difference. The difference between two antenna gains gives a coupling between the received power difference and the depar-ture angle of the signal.

same site, the next strongest cell is yold with index 2. Also assume that the rsrp measurements from different cells are conditionally independent to each other given an angle φ. The joint probability p(GA,1, GA,2|φ) is then simply the product

p(GA,1|φ)p(GA,2|φ) of the individual probabilities. The probability for y given an angle φ is then given as

p(y|φ) = p(GA1(φ) = GA2(φ) + Pth) N Y

i=3

p(GA1(φ) > GAi(φ)) (2.56)

in the case of transition to a cell on the same site, and

p(y|φ) = N Y

i=2

p(GA₁(φ) > GA_i(φ)) (2.57) when transitioning from a cell on another site. Note that a threshold Pthis used to account for transition hysteresis.

Equation (2.56) and Equation (2.57) look like Equation (2.16) on page 10 used in the rsrp measurement model. The difference here is that the path loss Gpwill cancel out and leave only the antenna gains left as the equations only deal with cells on the same site, which are assumed to share the same path loss.

(36)

0 50 100 150 200 250 300 350 0 2 · 10−2 4 · 10−2 6 · 10−2 8 · 10−2 0.1 0.12 0.14 Angle [degree] Probability Antenna 1-2

Figure 2.7: Angle probability. A transition between cells belonging to dif-ferent sites is presented. The angle probability is largest at aproximately 300◦_{, and looks in the proximity of this angle very much like a Gaussian}

distribution, which the rsrp error in the model used.

σ_DOD2 , the probabilities used in the transitions can be written as

p(GA,1(φ) = GA,2(φ) + Pth) = φN(GA,2+ Pth− GA,1|σ2= 2σDOD2 ) (2.58) and

p(GA,1(φ) > GA,i(φ)) = ΦN(GA,1− GA,i|2σDOD2 ), (2.59) where φN(x) and ΦN(x) denote the pdf in Equation (2.12) on page 9 and cdf in Equation (2.13) on page 10 of the normal distribution.

An example of angle probabilities using the above model for a cell transition is shown in Figure 2.7.

2.9.2 Voronoi Diagram

With an estimate of the dod, the set of positions which satisfy the angle is not bounded. To get a more restrictive measurement model, some measure of dis-tance would effectively limit the amount of possible positions.

Assume that each site has one omnidirectional antenna and that all sites trans-mit with the same average power. For any given position, the closest site would not only have the shortest distance but also have the strongest received signal

(37)

2.9 DOD Voronoi model 25 5 5.01 5.02 5.03 5.04 5.05 5.06 5.07 5.08 x ×105 1.56 1.57 1.58 1.59 1.6 1.61 1.62 1.63 1.64 y ×105

Figure 2.8: Voronoi diagram. Starting with a set of points called seeds, a corresponding set of regions of influence is created. All points inside a specific region are then closest to the seed of that region.

strength at the position if assuming free space signal propagation. With this rela-tion between the closest site and the received signal strength, a second bound to the position can be approximated.

Consider a diagram where all points in the two dimensional space belongs to a region of influence associated to a seed. In the case that each region contains the points which are closer to the seed of that region than any other seed, the diagram is called a Voronoi diagram. The Voronoi diagram can then be used as a model for the omnidirectional antenna assumption above. An example of a Voronoi diagram is shown in Figure 2.8.

The position of a cell transition can be said to belong to one of two cases when considering the Voronoi diagram. The position can be either on an edge E of a region if the new cell belongs to a different site, or inside a region if it belongs to the same site.

For these two cases, a simple algorithm is proposed to acquire a pdf for the positions of a cell transition. The algorithm uses the Voronoi diagram to produce a measure of probability for a cell transition y given a position x.

Using a Voronoi diagram based on the position of sites, the following defini-tions are made. The distance from the closest edge E to the point P is defined as

dP E, where P represents the position x. The standard deviation σV is defined as a multiple of dSE, which is the distance from the edge E to the seed S of the region. The magnitude of σV is dependent on how large position errors are expected for a cell transition.

(38)

P

S

E

Figure 2.9: Voronoi algorithm vectors. Given a point P and a region of

influence corresponding to a seedS, the edge point E closest to P is found.

The mean of the normal probability distribution is the distance|EP |, and the standard deviation isc|SE|.

pdfin Equation (2.12) on page 9 as

p(y|x) = φN(dP E|σ2= σV2) (2.60) in the case that the cell transition is between different sites. Using a small vari-ance in Equation (2.60), the position is correctly restricted to the edges of the region, while a larger variance loosens the requirements so that positions in the proximity of edges are allowed.

In the other case when the transition is between cells on the same site, the probability for x is given by the normal cdf in Equation (2.13) on page 10 as

p(y|x) = ΦN(dP E|σ2= σV2) (2.61) when the position x is inside the region, and as

p(y|x) = ΦN(−dP E|σ2= σV2) (2.62) when x is outside of the region. The different signs of dP Emakes the probability larger for points inside the region and smaller for points outside the region. The cdf represents the possibility of being on a specific side of an edge of a region given a variance. Using a small variance, the probability given by Equation (2.61) and Equation (2.62) is 1 inside the region while it is 0 outside the region, correctly implementing the description given earlier with the Voronoi diagram, while a larger variance allows for uncertainty of the assumed edges.

The values of the variables dP E and σV2 are more precisely defined in Algo-rithm 2.3 on the facing page which shows how to calculate the values. An illustra-tion of the algorithm vector distances is displayed in Figure 2.9. The probability for a position x for two example cell transitions is illustrated in Figure 2.10 on page 28.

(39)

2.9 DOD Voronoi model 27

Algorithm 2.3:Voronoi values dP E, dSE

For a specific region of influence, let Vifor i = 1, . . . , M denote the

coordinates of the vertices of the region. The edges of the region are then described by the origins Vi and the edge vectors

Ei = Vi+1− Vi, i = 1, . . . , M. (2.63) For a point P , the values dP Eand dSEare calculated with the following

algorithm.

1. Calculate the vectors V Pi from the vertices Vi for i = 1, . . . , M to the point P

V Pi = P − Vi (2.64)

2. For each i, create the projection ˜pE,iof the vector P Vion the edge vector

Ei

˜pE,i = V P_kEi· Ei

ik2 (2.65)

3. Limit the projections ˜pE,ito the bounds of the edges

pE,i= max{0, min{ ˜pE,i,kEik2}} (2.66)

4. The vector V Ei from the vertex Vito the point on the edge Eiclosest to P is then

V Ei = pE,iEi (2.67)

5. The shortest vector from P to any edge Eiis defined by

EPi = V Pi − V Ei (2.68)

6. Find the index of the edge closest to P

ic= arg min i kEPik2

(2.69) 7. Calculate the distances dP Eand dSE for the closest edge

dP E= E Pic ₂ (2.70) dSE = V Ei_c− (S − Vi_c) ₂ (2.71)

(40)

0 1 1.64 2 5.08 1.62 ×10-3 Probability 3 5.06 y ×105 4 1.6 x ×105 5.04 5 1.58 5.02 5 1.56

(a)Transition between sites.

0 0.2 1.64 0.4 5.08 1.62 Probability 0.6 5.06 y ×105 0.8 1.6 x ×105 5.04 1 1.58 5.02 5 1.56

(b)Transition on same site.

Figure 2.10: Voronoi probability.The position probability distributions for two scenarios are shown. The pdf for transitions between sites is focused along the edges, while the pdf for transitions on the same site is focused around the seed.

2.10 Filtering

In filtering of a dynamic system, a dynamic model is used to describe the be-haviour of the observed system. Together with a measurement model, the dy-namic model is used to combine the information of measurements over time and continuously estimate the current state of the system.

2.10.1 Dynamic Model

In general, the state xk of a system at time k is said to be affected by a process noise vkand can be described with

xk+1= f (xk, vk), vk∼ pvk, x0∼ px0, (2.72)

yk= h(xk) + ek, ek ∼ pek, (2.73)

where yk is a measurement dependent on the state and has a measurement error

ek. The functions f and h describes the dynamic model and the measurement model respectively.

The process noise vk and measurement noise ek will add uncertainty to the filter prediction based on the models. This uncertainty is described by the proba-bilities

xk+1∼ p(xk+1|xk), (2.74)

yk ∼ p(yk|xk). (2.75)

In many linear filter applications, only the measurement model is non-linear while the dynamic model is non-linear. The non-linear dynamic model can generally

(41)

2.10 Filtering 29 be expressed as

xk+1= Axk+ Bvk. (2.76)

In this thesis, the state of interest is the position. A model which tracks posi-tion and velocity with Gaussian acceleraposi-tion noise is used for the dynamic model. Using the linear notation in Equation (2.76) and assuming that the sample time is Ts, the linear model representing position and velocity can be expressed with the matrices

A =             1 0 Ts 0 0 1 0 Ts 0 0 1 0 0 0 0 1             , B =             Ts2/2 0 0 Ts2/2 Ts 0 0 Ts             . (2.77)

2.10.2 Bayesian Filtering

In Bayesian filtering, the goal is to estimate the marginal posterior probability

p(xk|y1:k), also known as the filtering distribution [16]. This is the probability of

the state xkconditioned on all collected observations y up to the time point k. A set of three recursive equations provides a solution to the filtering problem.

First, the filtering distribution is rewritten by expressing y_1:k as the combina-tion ykand y1:k−1and then using Bayes’ rule

p(xk|yk, y1:k−1) = p(yk|xk, y_p(y1:k−1)p(xk|y1:k−1)

k|y1:k−1) . (2.78)

Noting in Equation (2.78) that the observation history y_1:k−1 is redundant when xk is given, the expression for the filtering distribution is reduced to

p(xk|y1:k) = p(yk_p(y|xk)p(xk|y1:k−1)

k|y1:k−1) , (2.79)

which is the first of the three recursive equations and is called the measurement update (mu) equation because it uses the current measurement to update the filtering density.

The second equation calculates the denominator in Equation (2.79), which acts as a normalization constant to make the total probability integrate to unity. By using the law of total probability and removing the redundant observation history y_1:k−1given the state xk, the normalization constant is calculated as

p(yk, xk|y1:k−1) = p(yk|xk, y1:k−1)p(xk|y1:k−1), (2.80)

p(yk|y1:k−1) =

Z

Rnx

p(yk|xk)p(xk|y1:k−1)dxk, (2.81)

which is the integration of the first equation over the entire state space Rn_x_where

(42)

Finally, the recursion is completed by an equation known as the time update (tu) predicting the state xk+1with

p(xk+1, xk|y1:k) = p(xk+1|xk, y1:k)p(xk|y1:k), (2.82)

p(xk+1|y1:k) =

Z

Rnx

p(xk+1|xk)p(xk|y1:k)dxk. (2.83)

The equations describing the mu, its normalizing constant, and the tu are summarized in Algorithm 2.4. The algorithm is stated with the mu first and the tulast. This is easy to follow when the measurements arrive at regular intervals and the tu can be predicted with the correct time difference.

If the measurements passed to the filter algorithm do not have equal interval, the algorithm has to wait in the loop until it can do the prediction for the next measurement. A simple rearrangement of the tu and the mu can in the case of random time intervals avoid this by waiting between loop executions until the time interval is actually known. Such a rearrangement is made in the algorithm for the particle filter presented next.

Algorithm 2.4:Bayesian Filter

Having a set of T measurements yk, the state xkwith nxdimensions is approximated with the following recursion.

For k = 1, . . . , T

1. Measurement update:

p(xk|y1:k) = p(yk_p(y|xk)p(xk|y1:k−1)

k|y1:k−1) , (2.84)

where the normalization constant ckis

2.10.3 Particle Filter

The generic equations provided in the recursive Bayesian filtering algorithm is the solution to the filtering problem. In many applications, there are no easy analytical solutions. Instead, Monte Carlo methods provide a numerical way to solve for models which are neither linear nor Gaussian.

(43)

2.10 Filtering 31 The primary output from a non-linear filter is the posterior distribution

p(xk|y1:k), which can be used to estimate the state xk with the minimum mean squared error (MMSE) estimator

ˆxk|k=

Z

Rnx

xkp(xk|y1:k)dxk. (2.87)

The integral in Equation (2.87) can be approximated with Monte Carlo meth-ods using a finite set of samples to represent the posterior distribution. The sub-stitution problem is then a sum of the finite set of samples, which is simple to calculate.

In a particle filter, there are a set of particles xi_{, i = 1...N representing the} state space, where each particle is associated with a specific state xi _{and a weight}

wi [17]. These particles can be illustrated as a dynamic grid that adapts to the changing probability densities so that it favours state regions of high probabil-ity. At every sample time instance k, the weight associated with each particle is calculated to a relative probability

w_ki_|k−1∝ p(xk = xi|y1:k−1) (2.88)

so that the sum of all weights is unity N X

i=1

wi_k_|k−1= 1. (2.89)

Using the above notations for the particle filter, the filtering density p(xk|y1:k)

from the mu in Equation (2.79) on page 29 can with the finite set of particles xi and their weights wi_{be approximated by}

ˆp(xk|y1:k) = N X i=1 p(yk|xik)wik|k−1δ(xk− xki) ck (2.90)

where δ(x) denotes the Dirac delta function and the normalization constant ckis given by ck = N X i=1 p(yk|xki)wki|k−1. (2.91)

An example of the mu for a particle filter is displayed in Figure 2.11 on the next page.

Just like the notation for wi

k|k−1, the weight wki|k represents the relative

prob-ability for state xi

k to be the current state given the whole measurement history

y_1:k, and can from Equation (2.90) be seen to be

w_ki_|k = p(yk|x i k)wik|k−1

ck

(44)

4.91 4.92 4.93 4.94 4.95 4.96 4.97 4.98 4.99 5 · 105 0 1 2 3 4 5 6 7 8 9 · 10 −2 x Probability p(x_k| y1:k-1) p(y_k| xk) p(x_k| y1:k)

Figure 2.11: Particle filter measurement update.The measurement update combines the prediction distribution with the probability of a new measure-ment. The graph shows approximations of the distributions represented by the particles and weights. The prediction is simply multiplied with the mea-surement probability to achieve the new filtering distribution.

Network-Based Positioning Using Last Visited Cells Report

Master of Science Thesis in Communication Systems

Department of Electrical Engineering, Linköping University, 2016