# CEN/TR 17447:2020

(Main)## Space - Use of GNSS-based positioning for road Intelligent Transport System (ITS) - Mathematical PVT error model

## Space - Use of GNSS-based positioning for road Intelligent Transport System (ITS) - Mathematical PVT error model

This document is written in the frame of WP1.3 of GP-START project. It discusses several models to provide synthetic data for PVT tracks and the ways to analyse and compare the tracks to ensure these are similar to the reality.

## Mathematisches PVT-Fehlermodell

## Espace - Utilisation du positionnement GNSS pour les systèmes de transport routier intelligents (ITS) - Modèle d’erreur mathématique PVT

## Vesolje - Ugotavljanje položaja z uporabo GNSS za cestne inteligentne transportne sisteme (ITS) - Matematični model za napake PVT

### General Information

### Standards Content (Sample)

SLOVENSKI STANDARD

SIST-TP CEN/TR 17447:2020

01-april-2020

Vesolje - Ugotavljanje položaja z uporabo GNSS za cestne inteligentne transportne

sisteme (ITS) - Matematični model za napake PVT

Space - Use of GNSS-based positioning for road Intelligent Transport System (ITS) -

Mathematical PVT error model

Mathematisches PVT-Fehlermodell

Espace - Utilisation du positionnement GNSS pour les systèmes de transport routier

intelligents (ITS) - Modèle d’erreur mathématique PVT

Ta slovenski standard je istoveten z: CEN/TR 17447:2020

ICS:

03.220.20 Cestni transport Road transport

33.060.30 Radiorelejni in fiksni satelitski Radio relay and fixed satellite

komunikacijski sistemi communications systems

35.240.60 Uporabniške rešitve IT v IT applications in transport

prometu

SIST-TP CEN/TR 17447:2020 en,fr,de

2003-01.Slovenski inštitut za standardizacijo. Razmnoževanje celote ali delov tega standarda ni dovoljeno.

---------------------- Page: 1 ----------------------

SIST-TP CEN/TR 17447:2020

---------------------- Page: 2 ----------------------

SIST-TP CEN/TR 17447:2020

TECHNICAL REPORT

CEN/TR 17447

RAPPORT TECHNIQUE

TECHNISCHER BERICHT

February 2020

ICS 03.220.20; 33.060.30; 35.240.60

English version

Space - Use of GNSS-based positioning for road Intelligent

Transport System (ITS) - Mathematical PVT error model

Espace - Utilisation du positionnement GNSS pour les Mathematisches PVT-Fehlermodell

systèmes de transport routier intelligents (ITS) -

Modèle d'erreur mathématique PVT

This Technical Report was approved by CEN on 8 December 2019. It has been drawn up by the Technical Committee

CEN/CLC/JTC 5.

CEN and CENELEC members are the national standards bodies and national electrotechnical committees of Austria, Belgium,

Bulgaria, Croatia, Cyprus, Czech Republic, Denmark, Estonia, Finland, France, Germany, Greece, Hungary, Iceland, Ireland, Italy,

Latvia, Lithuania, Luxembourg, Malta, Netherlands, Norway, Poland, Portugal, Republic of North Macedonia, Romania, Serbia,

Slovakia, Slovenia, Spain, Sweden, Switzerland, Turkey and United Kingdom.

CEN-CENELEC Management Centre:

Rue de la Science 23, B-1040 Brussels

© 2020 CEN/CENELEC All rights of exploitation in any form and by any means Ref. No. CEN/TR 17447:2020 E

reserved worldwide for CEN national Members and for

CENELEC Members.

---------------------- Page: 3 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Contents Page

European foreword . 3

1 Scope . 4

2 Normative references . 4

3 Terms and definitions . 4

4 List of acronyms . 4

5 Approach . 5

6 Input data . 7

7 Model . 10

7.1 Introduction . 10

7.2 Discussion about models . 11

7.3 AR model . 16

8 Analysis tools . 29

8.1 General . 29

8.2 State of the art on data set similarity assessment . 30

8.3 Validation tool description . 34

9 Conclusions . 48

Annex A (normative) Neural network . 50

A.1 Neuron network concept . 50

Bibliography . 59

2

---------------------- Page: 4 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

European foreword

This document (CEN/TR 17447:2020) has been prepared by Technical Committee CEN-CENELEC/TC 5

“Space”, the secretariat of which is held by DIN.

Attention is drawn to the possibility that some of the elements of this document may be the subject of

patent rights. CEN shall not be held responsible for identifying any or all such patent rights.

3

---------------------- Page: 5 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

1 Scope

This document is written in the frame of WP1.3 of GP-START project. It discusses several models to

provide synthetic data for PVT tracks and the ways to analyse and compare the tracks to ensure these are

similar to the reality.

2 Normative references

The following documents are referred to in the text in such a way that some or all of their content

constitutes requirements of this document. For dated references, only the edition cited applies. For

undated references, the latest edition of the referenced document (including any amendments) applies.

EN 16803-1:2016, Space — Use of GNSS-based positioning for road Intelligent Transport Systems (ITS) —

Part 1: Definitions and system engineering procedures for the establishment and assessment of

performances

3 Terms and definitions

For the purposes of this document, the terms and definitions given in EN 16803-1:2016 apply.

ISO and IEC maintain terminological databases for use in standardization at the following addresses:

• ISO Online browsing platform: available at http://www.iso.org/obp

• IEC Electropedia: available at http://www.electropedia.org/

4 List of acronyms

ANOVA Analysis of variance

AR Autoregressive

ARMA Autoregressive moving average

CDF Cumulated distribution function

CET Central european time

DFT Direct Fourier transform

DOP Dilution of precision

FFT Fast Fourier transform

GNSS Global navigation satellite system

GPS Global positioning system

HDOP Horizontal dilution of precision

HPE Horizontal position error

IGS International GNSS service

ITS Intelligent transport systems

KS Kolmogorov–Smirnov

MFNN Multilayer feedforward neural networks

NED Northeast down

NN Neural network

4

---------------------- Page: 6 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

OBU On board unit

PVT Position velocity time

SBAS Satellite based augmentation system

SVM Support vector machine

UTC Universal time coordiante

VDOP Vertical dilution of precision

WGS84 World geodetic system 1984

WP Work package

5 Approach

The objective of the WP1.3 model the PVT tracks similar enough to the reality to be used for application

validation. The approach used is to twofold: to find a model and to build the tools to compare the

simulated data to the real data.

The following figure shows the overall logic of the work to be undertaken in WP 1.3:

Figure 1 — WP1.3 logic

The blue boxes represent the raw data that can be stored in files. The red boxes are the computations that

need to be carried out to get new data. Finally, the green boxes correspond to the outputs of this work

package that need to be documented.

5

---------------------- Page: 7 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

In some more detail:

— Reference data set: These are the reference trajectories computed with extra sensors and is expected

to have errors at centimetric level or better.

— Real data sets: These are the trajectories computed by the receiver under test.

— Error calculation: The error calculation computes the difference between the position computed by

the receiver and the one from the reference trajectory. Note that this step will be used to define the

parameters to be modelled. For instance in the Ifsttar paper [RD2], the model was based on the

absolute error and its angle. But one could use the position in XYZ (WGS84 coordinates). Since the

behaviour in the vertical axis is usually different, one could use a NED reference system or change

the reference to an axis in the direction of the movement another across the direction of the

movement and the third would be the vertical.

— Real error files: These files should contain the errors for each receiver trajectory.

— Error models: The error models will be described in the documentation and should depict the error

of the position depending on the previous errors, random variables (which can follow any

distribution) and parameters of the distributions.

— it algorithm: The fit algorithm should calculate the parameters that define the models.

— Model parameters: Storage of the models parameters.

— Synthetic error generation: Once the parameters are found, by starting the random variables with

different seed, one will get synthetic errors that should be equivalent to the real ones.

— Error analysis: As can be on the figure, the error analysis box appears twice: for the real errors and

for the synthetic ones. These analysis will be compared to validate the models. So far, the analysis

used have been the CDF (cumulative distribution function) and autocorrelation. As far as the CDF will

be representative for the error distribution in the position domain, in the time domain,

autocorrelation may not be enough to observe the behaviour of the error. It is suggested to compute

the DFT (discrete Fourier transform) to analyse the error in the frequency domain and observe if

there are errors at particular frequencies. Autocorrelation analysis will be kept for comparison with

previous models.

— Comparison: the Error analysis of real and synthetic data are compared to validate the models.

Synthetic data generator: this step adds the synthetic errors with the reference trajectory to create a new

data set that can be used for the sensitivity analysis.

6

---------------------- Page: 8 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

6 Input data

The input data to observe and models has been provided by Q-Free. It contains a header and then the

data with each column separated by a space. The header contains the list of fields

# SAVE Position Export

#

# Equipment: 1, seagull_datalogger_1

# Export Date: Tue Dec 13 13:58:20 CET 2016

#

# Space separated field list

#

# Equipment id (6 is reference)

# Drive id

# Leg id

# Date yyyy-mm-dd (UTC)

# Time hh:mm:ss.fff (UTC)

# Gps week

# Gps time of week, seconds

# OBU Timestamp, seconds

# GNSS Good flag (1 is good, 0 is bad)

# GNSS Fix Type

# Latitude

# Longitude

# Ellipsode height (meters)

# Geoid height, height over sea level (meters)

# Speed over ground (m/s)

# Course over ground (deg)

# Yaw (deg/second)

# Pitch (deg/second)

# Roll (deg/second)

# Horizontal pos error (meter)

# Vertical pos error (meter)

# sN, st.dev of position, in north direction (meters)

# sE, st.dev of position, in east direction (meters)

# sD, st.dev of position, in down direction (meters)

# HDOP

# PDOP

# NED n speed (m/s)

# NED e speed (m/s)

# NED d speed (m/s)

# Number of satellites used in solution

# Average signal quality

# SBAS Active

# rN integrity, north (m)

# rE integrity, east (m)

# rD integrity, down (m)

#

7

---------------------- Page: 9 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

To do the analyses and modelling of the PVT error the outputs computed by the receiver in a specific

proprietary way are not used. For instance, the horizontal and vertical position errors are calculated by

the receiver but one does not know how. Therefore, all the positions in latitude, longitude and height are

compared to the receiver #6 which is the reference. The other fields used are:

— Equipment Id and leg to separate the data sets;

— GNSS good flag and type fix to use only the positions with 4 (four) or more satellites that do not

use additional sensors (i.e.: type of fix 4 or 6).

For more details, the GNSS fix types are:

— NO_FIX = 0, → red

— ONE_SV = 1, → purple

— TWO_SV = 2, → blue

— THREE_SV = 3, → orange

— FOUR_PLUS_SV = 4, → green

— LS_2D = 5, → cyan

— LS_3D = 6, → dark green

— DEAD_RECKONING = 7, → yellow

— PSR_DIFF = 8, → brown

— SENSOR_FUSION = 9 → dark purple

The following figure shows the behaviour one can expect of a receiver (Reference is black) with a tunnel

in the middle of the path:

8

---------------------- Page: 10 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 2 — Oslo input data plotted

Some other issues have been observed in the input data. The following figure shows the histograms of

the error a data set from Frankfurt. The 6 (six) plots show the same data with different histogram bin

widths:

9

---------------------- Page: 11 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 3 — Frankfurt errors histograms

As it will be seen in the following clauses, the chosen model tends to a Gaussian distribution for data sets

big enough. However, the errors observed in Frankfurt data set, are not likely to follow a normal

distribution. When one observes the errors in the time domain, one can see that the error distribution

changes during the test. This is, the test starts in small error environment (an almost open-sky) and ends

in a great error environment (urban canyon).

Figure 4 — Frankfurt errors vs. time

7 Model

7.1 Introduction

Clause 7 presents the model chosen to create synthetic data representative of real PVT environments.

The first subclause discusses the previous model made by Ifsttar and Q-Free and the ARMA models

proposed in this document.

The second subclause presents the AR model chosen to fit the real data and the way its parameters are

calculated.

10

---------------------- Page: 12 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

7.2 Discussion about models

7.2.1 Laplace-steps and Cauchy-random walk

A previous model discussed is the one by Ifsttar and Q-Free that have already worked on the data sets to

be studied (See RD17). The following analysis intends to analyse which are the error characteristics to be

simulated and the limits of the model.

The model is twofold: on the one hand the absolute value of the horizontal error is calculated and on the

other the angle between the reference heading and the error direction. The radius (absolute error) is

calculated as follows:

R n with probability 1−p

( )

Rn 1 { + C n

r

L with probability p

r

p is the probability of a Bernoulli trial that in case of failure keeps the radius and in case of success will

draw a new radius from a Laplace distribution (positive values only) plus a Cauchy noise that accumulates

every time a step takes place.

NOTE On the radius calculation:

— Depending on the probability of the Bernoulli trial one can observe arcs of circle which is not a GNSS

positioning known behaviour.

— The probability of the Bernoulli trial allows mastering the correlation of the process since one can

manage the proportion of strongly correlated samples and those which are not.

— The Laplace distribution allows managing the mean and variance of the error.

— When one adds a cumulative error, the variances add up: This means that the cumulative error will

tend to infinity with time. On top of that the Cauchy distribution has undefined variance. So in the model a

minimum and maximum values need to be identified to avoid divergence.

Finally the error is filtered by an autoregressive filter.

The angle is calculated with a random walk:

∅+n 1 = ∅ nS+ n

∅

The noise chosen is a Cauchy distribution truncated with a min and a max values similarly to the Cauchy

distribution for the radius.

7.2.2 ARMA models

7.2.2.1 Introduction

The previous subclause has already discussed about an autoregressive filter. The idea of this subclause is

to present the characteristics of these models without the complex inputs (Laplace, Cauchy, random

walk…) of the previous subclause.

To do so, the input noise of the processes shown below will be a normal noise and the differences in the

outputs will entirely be due to the ARMA filter.

11

+=

---------------------- Page: 13 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

7.2.2.2 Definition

ARMA stands for autoregressive moving average. The autoregressive part means that it depends on the

previous values of the process and a weighting of the previous inputs.

An ARMA process can be defined as follows:

X a X +… + b N + bN + bN +…

nn1 −1 0 n 11n− 22n−

where

X can be the position in 1 (one) axis or the radius of the error;

N is a noise;

a and b are real numbers;

n is the index of the series.

One can use the Z transform to obtain:

Bz

( )

X= N

Az

( )

where

A and B are the polynomials with a and b coefficients;

The zeroes of B are the points at which there is no frequency response;

The zeroes of A are the points at which there are peaks in the response.

As a brief reminder, the roots of A need to be in the unit circle for the system to be stable. The closest to

the unit circle the bigger the impact of the root on the output.

7.2.2.3 Examples

7.2.2.3.1 General

Subclause 7.2.2.3 shows a few examples of the outputs that can be obtained with an ARMA process.

7.2.2.3.2 Gaussian noise

A Gaussian noise is hardly an ARMA process since it does not depend on the past information, but it is

important to see the response of a system when there is no filtering yet.

Its formula is the following:

X = bN

nn0

The outputs shown are the direct Fourier transform (DFT) and the time series. These two plots should

illustrate the outputs of the model.

12

=

---------------------- Page: 14 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 5 — DFT of a white noise

It is important to see that there is no specific resonance point: all the frequencies react just the same. Note

that the mean is not centred at zero. The frequencies observed go up to 0,5 Hz since the GNSS positions

are updated every second.

In the time domain, one observes the normal distribution behaviour.

Figure 6 — White noise against time

13

---------------------- Page: 15 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

7.2.2.3.3 AR process

For an AR process, there is no mean of the previous results. The solution depends only on the current

input and some previous results. This first example is the first order AR, so the output depends only on

the previous term.

X a X+ bN

nn1 −1 0 n1

Figure 7 — DFT of an AR process

One can observe that there is a peak close to the zero. This means that there is a root of the A polynomial

close to the unit circle. Indeed the root is (z-a) with 0 < a < 1. When a is close to 1, X is close to X

n n−1

with an added noise. This means that the correlation is strong between close samples. In the example

provided a is 0,8 and b is 1.

1 0

In the time domain, it can be observed that the samples change in some correlated way:

14

=

---------------------- Page: 16 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 8 — AR process against time

7.2.2.3.4 Peak example

The peak example is meant to show the behaviour of the DFT when the A polynomial contains complex

roots. To do so, one needs at least a second order AR process.

X= aX++aX b N

nn11−−22n 0 n

Figure 9 — DFT of AR process with 2 (two) poles

15

---------------------- Page: 17 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 10 — AR process with 2 (two) poles against time

In the example provided, a is −1,29, a is 0,64 and b is 1.

1 2 0

7.3 AR model

7.3.1 Model choice

The model chosen is an autoregressive one because it reduces the number of parameters needed for the

model and it can fit the data set errors calculated from a DFT point of view. The following figures show

the DFT transforms observed in different environments:

16

---------------------- Page: 18 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 11 — Frankfurt across error DFT

Figure 12 — Oslo across error DFT

17

---------------------- Page: 19 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 13 — Trondheim across error DFT

Figure 14 — Troms across error DFT

The previous figures show a behaviour similar to the AR model with a peak close to the zero that goes

down asymptotically to a certain value.

18

---------------------- Page: 20 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

7.3.2 Model parameters

7.3.2.1 Regression parameters

The AR equation is the following:

X a X+ bN

nn1 −1 0 n

The previous expression shows a model based on 2 (two) parts:

— aX : which represents the correlation (it makes one error depend on the previous one);

11n−

— bN : a random noise.

0 n

In order to figure the shape of the DFT one needs to find a . The equation for all N is:

1

XX N

2 1 1

X a X+ bN

3 12 0 2

XN

X

N NN−−11

One can remove the error term and solve the over-determined system with the mean-squares estimator:

−1

XX

12

a=XX……X X XX X X

( ) ( )

1 1 2 NN−−1 2 12 1 3

X X

N−1 N

N−1

xx

ii+1

∑

i=1

a =

1

N−1

2

x

i

∑

i=1

Check that a < 1 .

1

7.3.2.2 Noise parameter

7.3.2.2.1 Gaussian constant Noise parameter

Once the a has been figured out, one can calculate b . To do so, one needs to apply some basic

1 0

knowledge about GNSS positioning: the error will not go to infinity and it will not converge to a certain

number. From a statistics point of view, this means that the variance is constant.

Starting with the AR equation:

X a X+ bN

nn1 −1 0 n

22

var X var a X+ b N var a X+ var b N a var X+ b var N

nn 11−−0 nn 11 0 n 1 n−1 0 n

The variance of a normal standard distribution is 1. If we want the variance to remain constant:

var X var X var X

nn−1

19

==

= = =

=

=

=

---------------------- Page: 21 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

being the variance of the process regardless of the time.

var X

22

var X a var X+ b

10

2

b

0

var X =

2

1− a

1

In our case the calculation is made the other way round: we have the process variance that is calculated

from the real data, and one needs to calculate the model. Therefore:

2

b var X 1− a

( )

01

The noise parameter has been computed in 3 (three) cases:

— all the data;

— 99 % of the data (excluding the 1 % of biggest errors);

— 95 % of the data (excluding the 5 % of biggest errors).

The results provided by validation tool (see Clause 8 Analysis tools) are the following:

— all the data: 0 synthetic data sets are considered not similar to the real ones;

— 99 % of the data: 56,9 % not similar;

— 95 % of the data: 42,3 % not similar.

Tests have been performed with Oslo data sets which are representative of the urban environment.

7.3.2.2.2 Gaussian sliding noise parameter

From the previous subclause, it is understood that the variance remains constant through all the test.

However, it is observed that this is not the case in real data. Therefore, the model has been modified to

take into account the changing variance of the GNSS error.

n−1

2

1

ˆ

var X X− X

( )

n ∑ i

50

in−50

Indeed, the variance at each moment is calculated taking into account the last 50 values of the error.

Of course, b needs to be calculated each second of the data set and the formula remains almost the same:

0

2

b var X 1− a

( )

01nn

One of the main issues regarding the conception of a model is that one has to define the reality to model.

It is usually characterized by a type of receiver configuration and an environment. As far as the receiver

is easily identified, the environment is not. One of the parameters that describes this different

environments is the variance observed in each of those. Having a moving variance allows accounting for

the changes of environment end be representative of all the environments of the data set. Note that

although the scenarios have categories such as “urban”, it is impossible to guarantee that the buildings

around will always have the same height or be made out of the same material.

20

=

=

=

=

=

---------------------- Page: 22 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

With this sliding noise, the results provided by the analysis tool (see Clause 8 Analysis tools) are that only

7,3 % of the data sets generated are not similar to real ones they should emulate.

The problem of this model is that one needs data sets to build the moving variance instead of having a

defined distribution of a set of parameters. This makes the model very deterministic, since one will

always have the peaks at the same times. As this model is not kept, no further tests have been done with

variances taking into account more or less than 50 error values.

7.3.2.2.3 Experimental distribution

As mentioned previously, the model equation is:

X a X+ bN

nn1 −1 0 n

Then

N−1

xx

∑ ii+1

i=1

a =

1

N−1

2

x

i

∑

i=1

which is not time dependant. So, if we model properly,

We then model bN X− a X bN

0 nn 11n− 0 n

we should have a distribution that resembles the reality at the end. The first approximation is to get the

experimental distribution for the ahead and across axes.

The tool validation results says that only 4,9 % of the synthetic data are not similar to the real data that

it has to emulate.

Tests have been performed with Oslo data sets which are representative of the urban environment.

This is a good result, but it has one issue, when one plots the ahead errors against the across errors.

Figure 15 — Horizontal error

As can be seen on the figure, the errors are as big as expected, but they are always on the axes. This is

because the big errors are very rare and it is very improbable to have big error on both axes at the

same time.

21

=

=

---------------------- Page: 23 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

7.3.2.2.4 Experimental distribution with polar reference

7.3.2.2.4.1 Correlation factor for horizontal error

The first approximation is to compute the correlation of the horizontal error. Then one computes the

vector of the non-correlated error. This is, the difference between the horizontal error and the part of the

error due to the previous one. Finally, the cumulated distribution of the norm of non-correlated errors is

computed.

Figure 16 — Horizontal error correlation

To build the synthetic data sets, the opposite operation needs to be computed.

The norm of the error is retrieved randomly from the cumulated distribution and an angle at random

value between zero and 2 π. The vector with that norm and angle is added to the previous error times the

correlation factor (a1 of the AR model) to get the new error vector.

22

---------------------- Page: 24 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 17 — Error generation with horizontal error correlation

This model provides 65,04 % of non-similar synthetic data sets when compared to the real data sets that

should be similar.

Tests have been performed with Oslo data sets which are representative of the urban environment.

7.3.2.2.4.2 Correlation accounted differently for both axes

The reasoning is similar to the previous one, but in this case the errors are built considering that the

correlation is different for each axis:

23

---------------------- Page: 25 ----------------------

SIST-TP CEN/TR 17447:2020

CEN/TR 17447:2020 (E)

Figure 18 — Axis based correlation

Then the error is built accounting for the different correlations as well:

Figure 19 — Error generation with axis correlation

The number of non-similar synthetic data sets is 2,4 % compared to the real ones.

Tests have been performed with Oslo data sets which are representative of the urban environment.

7.3.2.2.4.3 Correlation between radius

**...**

SLOVENSKI STANDARD

kSIST-TP FprCEN/TR 17447:2019

01-november-2019

Vesolje - Ugotavljanje položaja z uporabo GNSS za cestne inteligentne transportne

sisteme (ITS) - Matematični model za napake PVT

Space - Use of GNSS-based positioning for road Intelligent Transport System (ITS) -

Mathematical PVT error model

Mathematisches PVT-Fehlermodell

Espace - Utilisation du positionnement GNSS pour les systèmes de transport routier

intelligents (ITS) - Modèle d’erreur mathématique PVT

Ta slovenski standard je istoveten z: FprCEN/TR 17447

ICS:

03.220.20 Cestni transport Road transport

33.060.30 Radiorelejni in fiksni satelitski Radio relay and fixed satellite

komunikacijski sistemi communications systems

35.240.60 Uporabniške rešitve IT v IT applications in transport

prometu

kSIST-TP FprCEN/TR 17447:2019 en,fr,de

2003-01.Slovenski inštitut za standardizacijo. Razmnoževanje celote ali delov tega standarda ni dovoljeno.

---------------------- Page: 1 ----------------------

kSIST-TP FprCEN/TR 17447:2019

---------------------- Page: 2 ----------------------

kSIST-TP FprCEN/TR 17447:2019

TECHNICAL REPORT

FINAL DRAFT

FprCEN/TR 17447

RAPPORT TECHNIQUE

TECHNISCHER BERICHT

September 2019

ICS 03.220.20; 33.060.30; 35.240.60

English version

Space - Use of GNSS-based positioning for road Intelligent

Transport System (ITS) - Mathematical PVT error model

Espace - Utilisation du positionnement GNSS pour les Mathematisches PVT-Fehlermodell

systèmes de transport routier intelligents (ITS) -

Modèle d'erreur mathématique PVT

This draft Technical Report is submitted to CEN members for Vote. It has been drawn up by the Technical Committee

CEN/CLC/JTC 5.

CEN and CENELEC members are the national standards bodies and national electrotechnical committees of Austria, Belgium,

Bulgaria, Croatia, Cyprus, Czech Republic, Denmark, Estonia, Finland, France, Germany, Greece, Hungary, Iceland, Ireland, Italy,

Latvia, Lithuania, Luxembourg, Malta, Netherlands, Norway, Poland, Portugal, Republic of North Macedonia, Romania, Serbia,

Slovakia, Slovenia, Spain, Sweden, Switzerland, Turkey and United Kingdom.

Recipients of this draft are invited to submit, with their comments, notification of any relevant patent rights of which they are

aware and to provide supporting documentation.

Warning : This document is not a Technical Report. It is distributed for review and comments. It is subject to change without

notice and shall not be referred to as a Technical Report.

CEN-CENELEC Management Centre:

Rue de la Science 23, B-1040 Brussels

© 2019 CEN/CENELEC All rights of exploitation in any form and by any means Ref. No. FprCEN/TR 17447:2019 E

reserved worldwide for CEN national Members and for

CENELEC Members.

---------------------- Page: 3 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

Contents Page

European foreword . 3

1 Scope . 4

2 Normative references . 4

3 Terms and definitions . 4

4 List of acronyms . 4

5 Approach . 5

6 Input data . 7

7 Model . 10

7.1 Introduction . 10

7.2 Discussion about models . 11

7.2.1 Laplace-steps and Cauchy-random walk . 11

7.2.2 ARMA models . 11

7.3 AR model . 16

7.3.1 Model choice . 16

7.3.2 Model parameters . 19

8 Analysis tools . 29

8.1 General . 29

8.2 State of the art on data set similarity assessment . 30

8.3 Validation tool description . 34

8.3.1 GP-START project context . 34

8.3.2 Validation tool overview . 34

8.3.3 Used data set similarity metrics . 35

8.3.4 Data set similarity decision process . 42

8.3.5 Data set similarity results . 44

9 Conclusions . 47

Annex A (normative) Neural network . 49

A.1 Neuron network concept . 49

A.1.1 General . 49

A.1.2 Neural network weight estimation . 51

Bibliography . 57

2

---------------------- Page: 4 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

European foreword

This document (FprCEN/TR 17447:2019) has been prepared by Technical Committee CEN-

CENELEC/TC 5 “Space”, the secretariat of which is held by DIN.

This document is currently submitted to the Vote on TR.

3

---------------------- Page: 5 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

1 Scope

This document is written in the frame of WP1.3 of GP-START project. It discusses several models to

provide synthetic data for PVT tracks and the ways to analyse and compare the tracks to ensure these are

similar to the reality.

2 Normative references

The following documents are referred to in the text in such a way that some or all of their content

constitutes requirements of this document. For dated references, only the edition cited applies. For

undated references, the latest edition of the referenced document (including any amendments) applies.

EN 16803-1:2016, Space — Use of GNSS-based positioning for road Intelligent Transport Systems (ITS) —

Part 1: Definitions and system engineering procedures for the establishment and assessment of

performances

3 Terms and definitions

For the purposes of this document, the terms and definitions given in EN 16803-1:2016 apply.

ISO and IEC maintain terminological databases for use in standardization at the following addresses:

• ISO Online browsing platform: available at http://www.iso.org/obp

• IEC Electropedia: available at http://www.electropedia.org/

4 List of acronyms

ANOVA Analysis of variance

AR Autoregressive

ARMA Autoregressive moving average

CDF Cumulated distribution function

CET Central european time

DFT Direct Fourier transform

DOP Dilution of precision

FFT Fast Fourier transform

GNSS Global navigation satellite system

GPS Global positioning system

HDOP Horizontal dilution of precision

HPE Horizontal position error

IGS International GNSS service

ITS Intelligent transport systems

KS Kolmogorov–Smirnov

MFNN Multilayer feedforward neural networks

NED Northeast down

NN Neural network

4

---------------------- Page: 6 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

OBU On board unit

PVT Position velocity time

SBAS Satellite based augmentation system

SVM Support vector machine

UTC Universal time coordiante

VDOP Vertical dilution of precision

WGS84 World geodetic system 1984

WP Work package

5 Approach

The objective of the WP1.3 model the PVT tracks similar enough to the reality to be used for application

validation. The approach used is to twofold: to find a model and to build the tools to compare the

simulated data to the real data.

The following figure shows the overall logic of the work to be undertaken in WP 1.3:

Figure 1 — WP1.3 logic

The blue boxes represent the raw data that can be stored in files. The red boxes are the computations that

need to be carried out to get new data. Finally, the green boxes correspond to the outputs of this work

package that need to be documented.

5

---------------------- Page: 7 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

In some more detail:

— Reference data set: These are the reference trajectories computed with extra sensors and is expected

to have errors at centimetric level or better.

— Real data sets: These are the trajectories computed by the receiver under test.

— Error calculation: The error calculation computes the difference between the position computed by

the receiver and the one from the reference trajectory. Note that this step will be used to define the

parameters to be modelled. For instance in the Ifsttar paper [RD2], the model was based on the

absolute error and its angle. But one could use the position in XYZ (WGS84 coordinates). Since the

behaviour in the vertical axis is usually different, one could use a NED reference system or change

the reference to an axis in the direction of the movement another across the direction of the

movement and the third would be the vertical.

— Real error files: These files should contain the errors for each receiver trajectory.

— Error models: The error models will be described in the documentation and should depict the error

of the position depending on the previous errors, random variables (which can follow any

distribution) and parameters of the distributions.

— it algorithm: The fit algorithm should calculate the parameters that define the models.

— Model parameters: Storage of the models parameters.

— Synthetic error generation: Once the parameters are found, by starting the random variables with

different seed, one will get synthetic errors that should be equivalent to the real ones.

— Error analysis: As can be on the figure, the error analysis box appears twice: for the real errors and

for the synthetic ones. These analysis will be compared to validate the models. So far, the analysis

used have been the CDF (cumulative distribution function) and autocorrelation. As far as the CDF will

be representative for the error distribution in the position domain, in the time domain,

autocorrelation may not be enough to observe the behaviour of the error. It is suggested to compute

the DFT (discrete Fourier transform) to analyse the error in the frequency domain and observe if

there are errors at particular frequencies. Autocorrelation analysis will be kept for comparison with

previous models.

— Comparison: the Error analysis of real and synthetic data are compared to validate the models.

Synthetic data generator: this step adds the synthetic errors with the reference trajectory to create a new

data set that can be used for the sensitivity analysis.

6

---------------------- Page: 8 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

6 Input data

The input data to observe and models has been provided by Q-Free. It contains a header and then the

data with each column separated by a space. The header contains the list of fields

# SAVE Position Export

#

# Equipment: 1, seagull_datalogger_1

# Export Date: Tue Dec 13 13:58:20 CET 2016

#

# Space separated field list

#

# Equipment id (6 is reference)

# Drive id

# Leg id

# Date yyyy-mm-dd (UTC)

# Time hh:mm:ss.fff (UTC)

# Gps week

# Gps time of week, seconds

# OBU Timestamp, seconds

# GNSS Good flag (1 is good, 0 is bad)

# GNSS Fix Type

# Latitude

# Longitude

# Ellipsode height (meters)

# Geoid height, height over sea level (meters)

# Speed over ground (m/s)

# Course over ground (deg)

# Yaw (deg/second)

# Pitch (deg/second)

# Roll (deg/second)

# Horizontal pos error (meter)

# Vertical pos error (meter)

# sN, st.dev of position, in north direction (meters)

# sE, st.dev of position, in east direction (meters)

# sD, st.dev of position, in down direction (meters)

# HDOP

# PDOP

# NED n speed (m/s)

# NED e speed (m/s)

# NED d speed (m/s)

# Number of satellites used in solution

# Average signal quality

# SBAS Active

# rN integrity, north (m)

# rE integrity, east (m)

# rD integrity, down (m)

#

7

---------------------- Page: 9 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

To do the analyses and modelling of the PVT error the outputs computed by the receiver in a specific

proprietary way are not used. For instance, the horizontal and vertical position errors are calculated by

the receiver but one does not know how. Therefore, all the positions in latitude, longitude and height are

compared to the receiver #6 which is the reference. The other fields used are:

— Equipment Id and leg to separate the data sets;

— GNSS good flag and type fix to use only the positions with 4 (four) or more satellites that do not

use additional sensors (i.e.: type of fix 4 or 6).

For more details, the GNSS fix types are:

— NO_FIX = 0, → red

— ONE_SV = 1, → purple

— TWO_SV = 2, → blue

— THREE_SV = 3, → orange

— FOUR_PLUS_SV = 4, → green

— LS_2D = 5, → cyan

— LS_3D = 6, → dark green

— DEAD_RECKONING = 7, → yellow

— PSR_DIFF = 8, → brown

— SENSOR_FUSION = 9 → dark purple

The following figure shows the behaviour one can expect of a receiver (Reference is black) with a tunnel

in the middle of the path:

8

---------------------- Page: 10 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

Figure 2 — Oslo input data plotted

Some other issues have been observed in the input data. The following figure shows the histograms of

the error a data set from Frankfurt. The 6 (six) plots show the same data with different histogram bin

widths:

9

---------------------- Page: 11 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

Figure 3 — Frankfurt errors histograms

As it will be seen in the following clauses, the chosen model tends to a Gaussian distribution for data sets

big enough. However, the errors observed in Frankfurt data set, are not likely to follow a normal

distribution. When one observes the errors in the time domain, one can see that the error distribution

changes during the test. This is, the test starts in small error environment (an almost open-sky) and ends

in a great error environment (urban canyon).

cou

Figure 4 — Frankfurt errors vs. time

7 Model

7.1 Introduction

Clause 7 presents the model chosen to create synthetic data representative of real PVT environments.

The first subclause discusses the previous model made by Ifsttar and Q-Free and the ARMA models

proposed in this document.

The second subclause presents the AR model chosen to fit the real data and the way its parameters are

calculated.

10

---------------------- Page: 12 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

7.2 Discussion about models

7.2.1 Laplace-steps and Cauchy-random walk

A previous model discussed is the one by Ifsttar and Q-Free that have already worked on the data sets to

be studied (See RD17). The following analysis intends to analyse which are the error characteristics to be

simulated and the limits of the model.

The model is twofold: on the one hand the absolute value of the horizontal error is calculated and on the

other the angle between the reference heading and the error direction. The radius (absolute error) is

calculated as follows:

R n with probability 1−p

( )

Rn 1 { + C n

r

L with probability p

r

p is the probability of a Bernoulli trial that in case of failure keeps the radius and in case of success will

draw a new radius from a Laplace distribution (positive values only) plus a Cauchy noise that accumulates

every time a step takes place.

NOTE On the radius calculation:

— Depending on the probability of the Bernoulli trial one can observe arcs of circle which is not a GNSS

positioning known behaviour.

— The probability of the Bernoulli trial allows mastering the correlation of the process since one can

manage the proportion of strongly correlated samples and those which are not.

— The Laplace distribution allows managing the mean and variance of the error.

— When one adds a cumulative error, the variances add up: This means that the cumulative error will

tend to infinity with time. On top of that the Cauchy distribution has undefined variance. So in the model a

minimum and maximum values need to be identified to avoid divergence.

Finally the error is filtered by an autoregressive filter.

The angle is calculated with a random walk:

∅+n 1 = ∅ nS+ n

∅

The noise chosen is a Cauchy distribution truncated with a min and a max values similarly to the Cauchy

distribution for the radius.

7.2.2 ARMA models

7.2.2.1 Introduction

The previous subclause has already discussed about an autoregressive filter. The idea of this subclause is

to present the characteristics of these models without the complex inputs (Laplace, Cauchy, random

walk…) of the previous subclause.

To do so, the input noise of the processes shown below will be a normal noise and the differences in the

outputs will entirely be due to the ARMA filter.

11

+=

---------------------- Page: 13 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

7.2.2.2 Definition

ARMA stands for autoregressive moving average. The autoregressive part means that it depends on the

previous values of the process and a weighting of the previous inputs.

An ARMA process can be defined as follows:

X a X +… + b N + bN + bN +…

nn1 −1 0 n 11n− 22n−

where

X can be the position in 1 (one) axis or the radius of the error;

N is a noise;

a and b are real numbers;

n is the index of the series.

One can use the Z transform to obtain:

Bz

( )

X= N

Az

( )

where

A and B are the polynomials with a and b coefficients;

The zeroes of B are the points at which there is no frequency response;

The zeroes of A are the points at which there are peaks in the response.

As a brief reminder, the roots of A need to be in the unit circle for the system to be stable. The closest to

the unit circle the bigger the impact of the root on the output.

7.2.2.3 Examples

7.2.2.3.1 General

Subclause 7.2.2.3 shows a few examples of the outputs that can be obtained with an ARMA process.

7.2.2.3.2 Gaussian noise

A Gaussian noise is hardly an ARMA process since it does not depend on the past information, but it is

important to see the response of a system when there is no filtering yet.

Its formula is the following:

X = bN

nn0

The outputs shown are the direct Fourier transform (DFT) and the time series. These two plots should

illustrate the outputs of the model.

12

=

---------------------- Page: 14 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

Figure 5 — DFT of a white noise

It is important to see that there is no specific resonance point: all the frequencies react just the same. Note

that the mean is not centred at zero. The frequencies observed go up to 0,5 Hz since the GNSS positions

are updated every second.

In the time domain, one observes the normal distribution behaviour.

Figure 6 — White noise against time

13

---------------------- Page: 15 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

7.2.2.3.3 AR process

For an AR process, there is no mean of the previous results. The solution depends only on the current

input and some previous results. This first example is the first order AR, so the output depends only on

the previous term.

X a X+ bN

nn1 −1 0 n1

Figure 7 — DFT of an AR process

One can observe that there is a peak close to the zero. This means that there is a root of the A polynomial

close to the unit circle. Indeed the root is (z-a) with 0 < a < 1. When a is close to 1, X is close to X

n n−1

with an added noise. This means that the correlation is strong between close samples. In the example

provided a is 0,8 and b is 1.

1 0

In the time domain, it can be observed that the samples change in some correlated way:

14

=

---------------------- Page: 16 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

Figure 8 — AR process against time

7.2.2.3.4 Peak example

The peak example is meant to show the behaviour of the DFT when the A polynomial contains complex

roots. To do so, one needs at least a second order AR process.

X= aX++aX b N

nn11−−22n 0 n

Figure 9 — DFT of AR process with 2 (two) poles

15

---------------------- Page: 17 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

Figure 10 — AR process with 2 (two) poles against time

In the example provided, a is −1,29, a is 0,64 and b is 1.

1 2 0

7.3 AR model

7.3.1 Model choice

The model chosen is an autoregressive one because it reduces the number of parameters needed for the

model and it can fit the data set errors calculated from a DFT point of view. The following figures show

the DFT transforms observed in different environments:

16

---------------------- Page: 18 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

Figure 11 — Frankfurt across error DFT

Figure 12 — Oslo across error DFT

17

---------------------- Page: 19 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

Figure 13 — Trondheim across error DFT

Figure 14 — Troms across error DFT

The previous figures show a behaviour similar to the AR model with a peak close to the zero that goes

down asymptotically to a certain value.

18

---------------------- Page: 20 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

7.3.2 Model parameters

7.3.2.1 Regression parameters

The AR equation is the following:

X a X+ bN

nn1 −10 n

The previous expression shows a model based on 2 (two) parts:

— aX : which represents the correlation (it makes one error depend on the previous one);

11n−

— bN : a random noise.

0 n

In order to figure the shape of the DFT one needs to find a . The equation for all N is:

1

XX N

2 1 1

X a X+ bN

3 12 0 2

X XN

N NN−−11

One can remove the error term and solve the over-determined system with the mean-squares estimator:

−1

XX

12

a=XX……X X XX X X

( ) ( )

1 1 2 NN−−1 2 12 1 3

X X

N−1 N

N−1

xx

ii+1

∑

i=1

a =

1

N−1

2

x

i

∑

i

=1

Check that a < 1 .

1

7.3.2.2 Noise parameter

7.3.2.2.1 Gaussian constant Noise parameter

Once the a has been figured out, one can calculate b . To do so, one needs to apply some basic knowledge

1 0

about GNSS positioning: the error will not go to infinity and it will not converge to a certain number. From

a statistics point of view, this means that the variance is constant.

Starting with the AR equation:

X a X+ bN

nn1 −10 n

22

var X var a X+ b N var a X+ var b N a var X+ b var N

nn 11−−0 nn 11 0 n 1 n−1 0 n

The variance of a normal standard distribution is 1. If we want the variance to remain constant:

var X var X var X

nn −1

19

==

= = =

=

=

=

---------------------- Page: 21 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

var X being the variance of the process regardless of the time.

22

var X a var X+ b

10

2

b

0

varX =

2

1− a

1

In our case the calculation is made the other way round: we have the process variance that is calculated

from the real data, and one needs to calculate the model. Therefore:

2

b var X 1− a

( )

01

The noise parameter has been computed in 3 (three) cases:

— all the data;

— 99 % of the data (excluding the 1 % of biggest errors);

— 95 % of the data (excluding the 5 % of biggest errors).

The results provided by validation tool (see Clause 8 Analysis tools) are the following:

— all the data: 0 synthetic data sets are considered not similar to the real ones;

— 99 % of the data: 56,9 % not similar;

— 95 % of the data: 42,3 % not similar.

Tests have been performed with Oslo data sets which are representative of the urban environment.

7.3.2.2.2 Gaussian sliding noise parameter

From the previous subclause, it is understood that the variance remains constant through all the test.

However, it is observed that this is not the case in real data. Therefore, the model has been modified to

take into account the changing variance of the GNSS error.

n−1

2

1

ˆ

var X X− X

( )

n ∑ i

50

in−50

Indeed, the variance at each moment is calculated taking into account the last 50 values of the error.

Of course, b needs to be calculated each second of the data set and the formula remains almost the same:

0

2

b var X 1− a

( )

01nn

One of the main issues regarding the conception of a model is that one has to define the reality to model.

It is usually characterized by a type of receiver configuration and an environment. As far as the receiver

is easily identified, the environment is not. One of the parameters that describes this different

environments is the variance observed in each of those. Having a moving variance allows accounting for

the changes of environment end be representative of all the environments of the data set. Note that

although the scenarios have categories such as “urban”, it is impossible to guarantee that the buildings

around will always have the same height or be made out of the same material.

20

=

=

=

=

=

---------------------- Page: 22 ----------------------

kSIST-TP FprCEN/TR 17447:2019

FprCEN/TR 17447:2019 (E)

With this sliding noise, the results provided by the analysis tool (see Clause 8 Analysis tools) are that only

7,3 % of the data sets generated are not similar to real ones they should emulate.

The problem of this model is that one needs data sets to build the moving variance instead of having a

defined distribution of a set of parameters. This makes the model very deterministic, since one will

always have the peaks at the same times. As this model is not kept, no further tests have been done with

variances taking into account more or less than 50 error values.

7.3.2.2.3 Experimental distribution

As mentioned previously, the model equation is:

X a X+ bN

nn1 −10 n

Then

N−1

xx

∑ ii+1

i=1

a =

1

N−1

2

x

i

∑

i=1

We then model bN X− a X which is not time dependant. So, if we model bN properly,

0 nn 11n− 0 n

we should have a distribution that resembles the reality at the end. The first approximation is to get the

experimental distribution for the ahead and across axes.

The tool validation results says that only 4,9 % of the synthetic data are not similar to the real data that

it has to emulate.

Tests have been performed with Oslo data sets which are representative of the urban environment.

This is a good result, but it has one issue, when one plots the ahead errors against the across errors.

Figure 15 — Horizontal error

As can be seen on the figure, the errors are as big as expected, but they are always on the axes. This is

because the big errors are very rare

**...**

## Questions, Comments and Discussion

## Ask us and Technical Secretary will try to provide an answer. You can facilitate discussion about the standard in here.