# An extensive analysis and calibration of the Modular Aggregation Algorithm across three categories of for GNSS trajectories data sources

Marie-Dominique van Damme, Yann Méneroux, Ana-Maria Olteanu-Raimond

► To cite this version:

Marie-Dominique van Damme, Yann Méneroux, Ana-Maria Olteanu-Raimond. An extensive analysis and calibration of the Modular Aggregation Algorithm across three categories of for GNSS trajectories data sources. Laboratoire sciences et technologies de l'information géographique. 2024. hal-04697576v3

HAL Id: hal-04697576

https://hal.science/hal-04697576v3

Submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Université
Gustave Eiffel

LASTI

<figure><img src="image_3.png" /></figure>

Lab on Geographic Information Science for
sustainable development and smart cities

<figure><img src="image_4.png" /></figure>

INSTITUT NATIONAL
DE L'INFORMATION
GÉOGRAPHIQUE
ET FORESTIÈRE

<figure><img src="image_5.png" /></figure>

ÉCOLE NATIONALE
DES SCIENCES
GÉOGRAPHIQUES

<figure><img src="image_6.png" /></figure>

# An extensive analysis and calibration of the Modular Aggregation Algorithm across three categories of GNSS trajectories data sources

Marie-Dominique Van Damme*, Yann Méneroux, Ana-Maria Olteanu-Raimond

Univ Gustave Eiffel, IGN-ENSG, LASTIG

73 Avenue de Paris

Saint-Mandé

France

29 novembre 2025

<footnote>* All authors contributed equally to this research.
version 2.0</footnote>
# Table des matières

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>1.1</td><td>Data availability</td><td>3</td></tr><tr><td>2</td><td>Data and material</td><td>4</td></tr><tr><td>2.1</td><td>Synthetic GNSS trajectories</td><td>4</td></tr><tr><td>2.2</td><td>Multi-sensors and multi-canopy traces acquisition</td><td>7</td></tr><tr><td>2.3</td><td>Crowdsourced traces</td><td>8</td></tr><tr><td>3</td><td>Formalization of the original algorithm</td><td>9</td></tr><tr><td>3.1</td><td>Algorithm initialization and termination</td><td>9</td></tr><tr><td>3.2</td><td>First step : trajectory matching</td><td>9</td></tr><tr><td>3.3</td><td>Second step : representative selection on trajectory sections</td><td>10</td></tr><tr><td>3.4</td><td>Third step : aggregation of representative points</td><td>10</td></tr><tr><td>4</td><td>MIAA : a Modular and Iterative Aggregation Algorithm</td><td>11</td></tr><tr><td>4.1</td><td>Component 1 - choosing the master trajectory</td><td>11</td></tr><tr><td>4.2</td><td>Component 2 - matching trajectories with the master trajectory</td><td>12</td></tr><tr><td>4.3</td><td>Component 3 - choosing the representative position of each homologous points</td><td>13</td></tr><tr><td>4.4</td><td>Component 4 - aggregating the representative position</td><td>13</td></tr><tr><td>5</td><td>Extensive analysis and calibration</td><td>14</td></tr><tr><td>5.1</td><td>Termination of the algorithm</td><td>14</td></tr><tr><td>5.2</td><td>Algorithm calibration : further details</td><td>15</td></tr><tr><td>5.2.1</td><td>Step 1 : Create reference tracks</td><td>15</td></tr><tr><td>5.2.2</td><td>Step 2 : Create sets of simulated GNSS trajectories</td><td>16</td></tr><tr><td>5.2.3</td><td>Step 3 : Compute aggregated trajectory</td><td>16</td></tr><tr><td>5.2.4</td><td>Step 4 : Evaluate error between the aggregated and ground truth track</td><td>18</td></tr><tr><td>5.3</td><td>Crowd-sourced trajectories dataset</td><td>20</td></tr><tr><td>5.4</td><td>Perspective : performance under extreme conditions</td><td>21</td></tr></table>
Table des figures

<table><tr><td>1</td><td>Example of a non-realistic simulation of a GNSS trajectory with a white noise process, completely missing out the true correlation pattern of GNSS measurements.</td><td>5</td></tr><tr><td>2</td><td>Illustration of the three error components in the GNSS trajectories. From left to right : (1) a long wave-length process describing coordinate system errors, (2) an intermediate wave-length process describing GNSS observation errors (auto-correlated in space and time) and (3) a white noise process (e.g. heat, vibrations, electronic noise).</td><td>5</td></tr><tr><td>3</td><td>Example of 5 synthetic GNSS trajectories (blue) generated on a common ground truth track (dashed line).</td><td>6</td></tr><tr><td>4</td><td>Set of 50 trajectories collected with five sensors in dense forest and ground truth route.</td><td>7</td></tr><tr><td>5</td><td>Trajectories in three different contexts with spatial constraint : (a) ridge (terrain constraint), (b) series of sharp (infrastructure constraint), (c) heterogeneous trail shape (scale variation).</td><td>8</td></tr><tr><td>6</td><td>(a) A matching link between trajectory <i>X</i> and master trajectory <i>R</i>. (b) A connected component of the bipartite graph associated with the matching. (c) representative and aggregation positions for a reference position of the master trajectory <i>R</i>.</td><td>10</td></tr><tr><td>7</td><td>The four components of modular and iterative aggregation algorithm for GNSS trajectories.</td><td>11</td></tr><tr><td>8</td><td>Heuristics for choosing a trajectory master.</td><td>11</td></tr><tr><td>9</td><td>Heuristics for choosing a similarity measure for matching trajectories</td><td>12</td></tr><tr><td>10</td><td>Heuristics for options in [C3] for DTW L<sub>∞</sub></td><td>13</td></tr><tr><td>11</td><td>Heuristics for options in [C4] for DTW L<sub>∞</sub></td><td>13</td></tr><tr><td>12</td><td>Counter-example showing a case of non-convergence</td><td>14</td></tr><tr><td>13</td><td>Aggregated trajectories for two samples of 3 and 20 trajectories and computed with Discrete Fréchet distance (in pink) and DTW-L<sub>2</sub> distance (in blue)</td><td>17</td></tr><tr><td>14</td><td>Position error measurement for each samples of N trajectories generated</td><td>18</td></tr><tr><td>15</td><td>Shape deviation measurement for each samples of N trajectories generated</td><td>19</td></tr><tr><td>16</td><td>Aggregated trajectory result on Crowd-sourced dataset</td><td>20</td></tr></table>
# 1 Introduction

With the growing interest in outdoor recreation, there has been a surge in sport associations proposing new group outdoor activities and offering new opportunities for people to explore and engage with nature in exciting ways. This expansion contributes to the economic development of local communities and, more generally, to tourism. However, the expansion of outdoor activities could have negative consequences on the preservation of biodiversity. Indeed, scientific ecologists mentioned "the landscape of fear" to describe the negative effects on wildlife (e.g., chamois) disturbance [GPML22] or loss of plant species diversity.

In this context, it becomes crucial to better estimate the human pressure in time and space to guide informed management and conservation actions as well as to compute SDGs indicators. Such data measuring the presence of people, animals, and their interactions or protected plants are collected by local and regional stakeholders. However, a previous report showed a lack of official data and disparities between countries at global scale. In parallel, recent research has demonstrated the potential of crowdsourcing and citizen science data to enhance the representation of human pressure on ecosystems such as GNSS trajectories collected and shared openly by practitioners on different platforms (e.g. CamptoCamp, IGNRando, OutdoorVision, OpenStreetMap, Geo-life) as well as data provided by citizens through Citizen Science (CS) initiatives.

In this context, our aim is to develop innovative methods and tools to assess human pressures on ecosystems and provide new solutions to help stakeholders mitigate these impacts while promoting sustainable tourism and outdoor activities. More concretely, the goal is to build a framework allowing integration of spatiotemporal data representing a network of outdoor activities derived from road network and GNSS trajectories. We aim to enhance the semantic of this network with human frequency and activities types, biodiversity information (e.g. rest areas for animals) as well as human-wildlife encounters. One of the challenges is to deal with data heterogeneity originating from multiple sources across different spatial and temporal scales and from various types of data producers, including national, regional, and local authorities, research communities, and citizens (through participatory science or volunteered geographic information). Both input and output data, as well as the proposed methods, will be described in a knowledge graph to facilitate data retrieval, dissemination, and reproducibility.

This work is part of the research project IntForOut (Multisource spatial data INTEGRation FOR the Monitoring of Ecosystems under the pressure of OUTdoor recreation) funded by the French National Research Agency (n° ANR-23-CE55-0003).

This technical report aims to complement the conference paper [VDMOR24] by providing additional experiments or further details that could not be included in the paper.

## 1.1 Data availability

- This extensive analysis is based on the work described in [VDMOR24];

- The open-source implementation of the MIAA algorithm in the tracklib library [MvD24];

- Multi-sensors and multi-canopy traces acquisition with the ground truth [RMC+25]
## 2 Data and material

Since our goal is to reconstruct the geometry of the path from a given subset of GPS trajectories, disregarding the time dimension, each dataset used in this study contains trajectories following the same route, travelling in the same direction, and having approximately the same start and end points. In this study, the data used consist of a collection of datasets that respect these conditions, aligning with those outlined in [ED14].

The scenarios of greatest interest for studying the aggregation method are paths under forest cover, due to the noise effect, and those containing a series of sharp handling challenges in aggregation. Moreover, our analysis will be carried out by using both synthetic and real trajectories. Thus three categories of data are considered:

1. Realistic Synthetic GNSS Trajectories: These simulated trajectories provide a controlled environment to test the algorithm's performance under various predefined conditions, allowing us to assess its accuracy and robustness;

2. Multi-Sensor Trajectories: data collected with different sensors (professional, watched and smartphone applications) to evaluate how well the algorithm handles the heterogeneity of sensor data. Multi-sensors trajectories are acquired by following a repeatable data collection protocol, we defined;

3. Real world data (i.e. crowdsourced datasets), more specifically three crowdsourced datasets categories (terrain-constrained itinerary, infrastructure-constrained and multi-scale); this choice is not strictly «metrological» but rather qualitative : to investigate graphically how the algorithm performs on real typical cases.

### 2.1 Synthetic GNSS trajectories

In order to assess the metrological performances of the algorithm (i.e. its ability to reconstruct accurately the common path followed by all the individual sample trajectories) it is required to test it on a large array of configurations, with a substantial amount of GPS trajectories. This raises two operational problems :

~~ 1. First, to compare the estimated trajectory with the real route actually followed by the individual samples, it is required to perform a costly and time-consuming on-the-field survey to measure the ground truth. This is especially problematic since moderately to densely covered forest areas are of interest, where the direct use of precision carrier-phase GNSS with Real-Time Kinematic fast and convenient procedure is often impossible. Therefore, ground truth has to be measured with classical surveying traverse, tied to an absolute reference point, possibly located hundreds of meters away; spending a full day of topometric survey for each case study is not realistic. ~~

~~ 2. Secondly, for each case study, individual GNSS trajectories must be collected by walking several times along the ground truth route. Because GNSS error is known to be auto-correlated in time [Rob93], for a realistic acquisition of data, trajectories must ideally be sampled at different times of day, which increases again drastically the effort needed to collect real data. ~~

To overcome these limitations, it was first decided to proceed to extensive experimentation of the algorithms on simulated GNSS trajectories. This enables to simulate as many case studies as needed, with potentially unlimited number of GNSS trajectories, and a readily available ground truth track to compare the results with.

This methodology however, requires an accurate modeling of auto-correlation error of GNSS trajectories, to avoid non-realistic simulations, as depicted for example in Fig. 1, which can also result in topological errors [Bon02, Vau97].

GNSS errors were then modeled through their covariance function
<figure><img src="image_7.png" /><figcaption>FIGURE 1 – Example of a non-realistic simulation of a GNSS trajectory with a white noise process, completely missing out the true correlation pattern of GNSS measurements.</figcaption></figure>

$$ \gamma(s_1, s_2) = \mathrm{Cov}(X(s_1), X(s_2)) $$

describing the statistical covariance between positioning errors $X(s_1)$ and $X(s_2)$ at two locations $s_1, s_2 \in \mathbb{R}_+$ (described trough their curvilinear abscissa along the ground truth trajectory).

Further, the error $X$ is supposed to be a second-order stationary process (hence described only by the difference $s_2 - s_1$), and is modelled to take into account different error components in the GNSS trajectory measurement process, as illustrated on Fig 2.

<figure><img src="image_8.png" /><figcaption>FIGURE 2 – Illustration of the three error components in the GNSS trajectories. From left to right : (1) a long wave-length process describing coordinate system errors, (2) an intermediate wave-length process describing GNSS observation errors (auto-correlated in space and time) and (3) a white noise process (e.g. heat, vibrations, electronic noise).</figcaption></figure>

Generation of GNSS trajectories was done independently on each of the two planimetric components, with a methodology described in [Rip09] and also employed in [MMALGOR23] : with a random generator, we sampled n i.i.d. unit-variance and zero-mean gaussian values, compiled in a vector $x$. It can easily be shown that, for any positive-definite matrix $\Sigma \in \mathbb{R}^{n \times n}$, the random vector $y = Ax$ where $A$ is a Cholesky factor of $\Sigma$, is a realization of a correlated random vector $Y$ having covariance matrix $\Sigma$. The covariance matrix $\Sigma$ is formed with $\Sigma_{ij} = \gamma(s_j - s_i)$ (Cf. example in Figure 3).

**Application** : In [VDMOR24], the used methodology is the following : for each case study, a reference track
<figure><img src="image_9.png" /><figcaption>FIGURE 3 – Example of 5 synthetic GNSS trajectories (blue) generated on a common ground truth track (dashed line).</figcaption></figure>

is simulated (or extracted from an existing topographic database) and is considered as the ground truth track from which all GNSS trajectories are simulated. The error between the estimated and ground truth track is then evaluated, which in turn, enables to assess the sensitivity of the algorithm to all its parameters. Similar methodologies have been used for example by [BHLS15] and [ZY15].

In our experimentation, trajectories have been generated with a noise generated with a 100 m range Gaussian Process and with a 50 cm-amplitude for reference error, completed by a noise generated with a 5 m-amplitude exponential covariance process (GPS error) ([GBTY05]) with a 50 m correlation scope and completed by 1 m white noise process (vibrations, electronic noise, etc.).

Which translates to the following code using the `Tracklib` library ([MvD24]):

```python
# Generate a track
tkl.seed(123)
base_lacets = tkl(generate(0.4, dt=10)
chemin = tkl.noise(base_lacets, 20, tkl.SincKernel(20),
                    direction=tkl.MODE_DIRECTION_ORTHO)[:3]
chemin = chemin[80:250]
chemin.scale(10)

# Generate a collection of five identical tracks
N = 5
tracks3 = tkl.core.TrackCollection([chemin]*N)

# Generate noise to get realistic GNSS trajectory
tracks3.noise(0.5, tkl.GaussianKernel(100))
tracks3.noise(5, tkl.ExponentialKernel(50))
tracks3.noise(1, tkl.DiracKernel())

# Visualisation
plt.plot(chemin.getX(), chemin.getY(), color="black", linewidth=-1, linewidth=1)
for track in tracks3:
    plt.plot(track.getX(), track.getY(), color="royalblue", linewidth=-1, linewidth=1)
```
## 2.2 Multi-sensors and multi-canopy traces acquisition

To assess the impact of canopy and the sensors on precision accuracy, we defined and implemented a data collection protocol.

First, based on the literature, we identified three types of canopy (i.e. open area, moderate coverage, and heavy coverage). Additionally, we delineated five types of sensors (i.e. mobile phone equipped with Visio-Rando application, Polar GPS device, Garmin GPS device, Keymaze device and professional Ublox GPS sensor chip).

Second, for each type of canopy, the following data collection protocol was defined :

- Identification of areas without spatial constraints (e.g. bridge, stream, unobstructed)

- Identification within a rectangle of an Origin/Destination route with moderate winding and approximate length of 300 m.

- Collect five round trips following the route exactly.

Third, the field work was done by two of the authors of this paper. The placement of the sensors is also relevant according to the literature [BKT11]. To limit this effect, the GPS watches were worn on the wrist, the professional GPS devices were carried in a bag with the antenna positioned externally, and the mobile phones were held in the hand. Data collection has been carried out during the summer season. In total, for the three types of canopy and five sensors, 150 trajectories are collected.

**Application** : In [VDMOR24], our experiments focus only on dense forest 50 trajectories collected with five sensors are shown in Figure 4.

<figure><img src="image_10.png" /><figcaption>FIGURE 4 – Set of 50 trajectories collected with five sensors in dense forest and ground truth route.</figcaption></figure>

The ground truth route is obtained through a topometric survey conducted by a group of students as part of a topometric project, under the supervision of the authors of this paper and their teachers [CLP24]. Absolute positioning is conducted through GNSS differential static positioning of a set of reference points located in open-sky conditions, between 300 and 500 meters from the surveyed route. Topometric determination of the route geometry is performed with surveying traverse. The output ground truth is sampled with about 42 points (i.e. about 1 point every 7 meters) with an absolute positioning accuracy of 5 mm in each 3D axis (1σ) (Cf. Figure 4).
## 2.3 Crowdsourced traces

To validate the results obtained with synthetic GNSS tracks and data acquired according to the proposed protocol, we have chosen to test the algorithms on crowdsourced trajectories, as these will be used to derive pressure and route frequency indicators for the end-users.

To this end, our third experiment was carried out considering tracks downloaded from Visorando<fnref n="1" /> and from Wikiloc<fnref n="2" />, websites offering donwloading tracks published online by contributors. The GPS sensors used are therefore unknown, as well as data changes, context conditions.

In this experiment, we focused on hiking activity in the mountain area of the Pralognan Valley in the French Alps. Specifically, from all the trajectories, we selected some of them having specific and challenging configurations and spatial constraints: terrain constraints (e.g. ridge, river), infrastructure support (e.g. switchback trail) and variations of the shape of the path (e.g. a straight line followed by a series of twists and turns with varying distances between them). Note that, in a dataset, all the selected trajectories have the same origin/-destination (Cf. Figure 5). Note that, traces may have been previously filtered via a simplification algorithm like Douglas-Peucker. Indeed, all geometries offer the same characteristics like vertex of bends, turns in route. If this was not the case, traces would be more dissimilar in shape.

<figure><img src="image_11.png" /><figcaption>FIGURE 5 – Trajectories in three different contexts with spatial constraint: (a) ridge (terrain constraint), (b) series of sharp (infrastructure constraint), (c) heterogeneous trail shape (scale variation)</figcaption></figure>

As the reviewers suggested in [VDMOR24], considering a broader approach based on real world data would be interesting. This could be done by generalizing the first context: a route with terrain constraints is a corridor, covered by an acceptable number of traces, for which we can guarantee they have actually followed the path. For example, in addition to ridge paths, constrained route segments can include trails along the mountainside, trail sections bordered by fences or rock walls. This can be addressed by defining a protocol to extract such trajectories, involving running spatial and topological queries to identify such areas, and then cross-referencing them with traces dataset. To ensure the representativeness of such constraint contexts, we use crowdsourced data issued from collaborative platforms. Eventually, the geometry of the road section stemming from an accurate topographic database is used as the ground truth. Alternatively, ground truth dataset may be collected directly on the field with topometric surveying. The algorithm will be run on all extracted constrained route segments, providing a quasi-metrological validation on real data. Although this approach is relevant, we could not further develop it due to the paper length constraints, so it will be explored in future works.

<footnote n="1">1. www.visorando.com</footnote>

<footnote n="2">2. www.wikiloc.com</footnote>
## 3 Formalization of the original algorithm

This section formalizes the algorithm proposed [ED14] lacking in a mathematical description to better describe its properties.

Considering a set of trajectories, where each trajectory is defined as an *n* ordered points :

$$X = (x_i)_{i=1..n}, \text{ with } x_i \in \mathbb{R}^2 \text{ and GPS records in a 2D space.}$$

Let *d* be a distance (Euclidian, Manhattan, etc.) between those points :

$$d: \mathbb{R}^2 \times \mathbb{R}^2 \rightarrow \mathbb{R}_+.$$

Note that heights and timestamps for example, could also be considered with **x**<sub>*i*</sub>, **y**<sub>*i*</sub> ∈ ℝ<sup>3</sup> or ℝ<sup>4</sup>, and defining *d* accordingly to measure distance between three-dimensional timestamped records.

An **aggregated trajectory**, noted *AX*, is defined as the best geometric representation of a set of trajectories *X* following exactly the same route defined from an origin to a destination :

$$\underline{\mathcal{AX}} = (\overline{xj}), j=1..m,$$

where *xj* represents the aggregated points of matched points for the *j*-th point in the master trajectory.

We defined an **accurate aggregated trajectory**, noted *AAX*, an aggregated trajectory that optimizes a quality criterion *Q*, with respect to the ground truth *G* (unknown). More formally, *AAX* optimizes *E[Q]*, where *E* is the expected value of the quantity *Q*. The main challenge is to define *Q*, and then to compute it in an approximate and satisfactory way, which is our focus.

One might say, a trace **T**<sub>1</sub> is a good partial representation of **T**<sub>2</sub> if **Q**<sub>**T**<sub>2</sub>→**T**<sub>1</sub></sub> is minimal with *Q* being the square root of the mean of the squared distances between each point in **T**<sub>1</sub> and its closest neighbor on **T**<sub>2</sub>.

Thus, the quality of the trace *AAX* can then be evaluated from the average of the partial qualities :

$$\mathcal{AX} = (\mathcal{Q}_{\mathcal{AX} \rightarrow G} + \mathcal{Q}_{G \rightarrow \mathcal{AX}}) \div 2, \text{ with } \mathcal{G} \text{ the ground truth.}$$

### 3.1 Algorithm initialization and termination

As mentioned before, this algorithm uses an iterative refinement approach that improve the existing solution at each step, continuously performing multiple matches. And each iteration is composed of three steps : trajectory matching, representative selection on trajectory sections and aggregation of representative points.

But first algorithm initialization consists of choosing a first trajectory, will be called the master trajectory, *R* = (r<sub>*j*</sub>)<sub>*j*=1..*m*</sub>. This process involves defining a set of reference positions ((r<sub>*j*</sub>)<sub>*j*=1..*m*</sub>), which are all the vertices of the master polyline. Selecting a master trajectory, to start the algorithm, avoids having to match all the trajectories in pairs.

And finally, the iterative process stops when the difference between the aggregated trace and the master trace is below a minimum threshold.

### 3.2 First step : trajectory matching

The points of each trajectory will be matched with the master trajectory. To be more precise, a matching between the trajectory *X* and the master trajectory *R* is a bipartite graph (*X*, *R*, *μ*) with no isolated vertex and where *μ* are the set of matching links. Note that the connectivity constraint imposes that each point in *X* is linked to at least one points in *R*, and reciprocally. We denote by *M(n, m)* the set of all possible matching between two trajectories containing *n* and *m* points respectively.

An **ordered matching** between trajectories *X* and *R* is a matching with the additional constraint that there should be no pair of crossing links, i.e. that we cannot find two links (x<sub>*i*</sub>, r<sub>*j*</sub>) and (x<sub>*k*</sub>, r<sub>*l*</sub>) with *i* > *j* and *k* < *l*. We denote by *M<sup>+</sup>(n, m)* the set of all possible ordered matchings between two trajectories containing *n* and
m points respectively. Note that the definition of an *ordered matching* imposes that $(x_1, r_1)$ and $(x_n, r_m)$, i.e. start and end points of trajectories $\mathcal{X}$ and $\mathcal{R}$ should be matched together respectively.

Matching algorithms based on dynamic programming approach (*Dynamic Time Warping matching* (DTW), *Discrete Fréchet distance*) are used to calculate ordered matching $\mu \in M^+(n, m)$. The $L_p$-norm optimal dynamic time warping matching between the trajectory $\mathcal{X}$ and the master trajectory $\mathcal{R}$ is given by :

$$DTW_p(\mathcal{X}, \mathcal{R}) \in \underset{\mu \in M^+(|\mathcal{X}|, |\mathcal{R}|)}{\operatorname{argmin}} \sum_{i=1}^{|\mu|} d(\mathbf{x}_{\mu_i(1)}, \mathbf{r}_{\mu_i(2)})^p$$

where, for any index $i \in [1, |\mu|]$, the point of index $\mu_i(1)$ of trajectory $\mathcal{X}$ is linked to the point of index $\mu_i(2)$ of $\mathcal{R}$.

Note that the discrete Fréchet distance (the method used by [ED14]) is similar, but not identical to DTW, since Discrete Fréchet distance seeks to minimize the maximum of distances between the two curves, whereas DTW aims at minimizing the sum of these distances. In fact, the discrete DTW distance given in the formula above degenerates to Fréchet distance when $p$ grows to the infinity.

<figure><img src="image_12.png" /><figcaption>FIGURE 6 – (a) A matching link between trajectory $\mathcal{X}$ and master trajectory $\mathcal{R}$. (b) A connected component of the bipartite graph associated with the matching. (c) representative and aggregation positions for a reference position of the master trajectory $\mathcal{R}$.</figcaption></figure>

A connected component of the bipartite graph $(\mathcal{X}, \mathcal{R}, \mu)$ corresponds exactly to one or more matching links with $ni$ positions of $\mathcal{X}$ and $mj$ positions of $\mathcal{R}$. For example, in Figure 6-a and 6-b, the $j$-th connected component is composed of vertex $(rj_1)$ from $\mathcal{R}$, vertices $(xi_1$ et $xi_2)$ and the two associated vertex.

### 3.3 Second step : representative selection on trajectory sections

At the end of matching process, positions of trajectory $\mathcal{X}$ can be linked to many points of $\mathcal{R}$ and reciprocally. Thus, the second step of algorithm iteration is to choose for each connected component previously established, a representative position of each group of vertices from the trajectory $\mathcal{X}: \overline{xj}$. This ensures that the final aggregation is not too much influenced. If the segment of the trajectory is composed of several vertices, the representative is calculated from all the positions of the segment, otherwise the single position is chosen. [ED14], take the center of gravity of the segment. In Figure 6-c, the red circle $\overline{xj}$ represents the center of gravity of the vertices $xi_1$ and $xi_2$).

### 3.4 Third step : aggregation of representative points

At this point, all positions of the master trajectory $(\mathbf{r}_j)_{j=1..m}$ have a unique homologous point $\overline{xj}$ from each trajectory $\mathcal{X}$. The positions of the new merged trajectory are calculated from the median of matching points. The median aggregation operator is better suited to handle outliers. [ED14] add a constraint in this step of the algorithm : each position aggregated must also be part of existing trajectory positions, so as not to be located in an unlikely place.

At last, we iterate using the aggregated trajectory as the new master trajectory. The algorithm stops when the distance between two subsequent estimation of the aggregated trajectory is below a predefined threshold. The distance used is the pointwise $L_2$ distance.
# 4 MIAA: a Modular and Iterative Aggregation Algorithm

Knowing the diversity of similarity and aggregation method for trajectories, we adapt the existing algorithm proposed in [ED14] by transforming it into a modular one. Modularity has substantial advantages such as flexibility (i.e. a measure can be easily replaced by another), scalability (i.e. new measures can be easily added), and more relevant to the goal of this study, experiment testing (i.e. study the behavior and the influence of measures and parameters in different contexts and with different data).

<figure><img src="image_13.png" /><figcaption>FIGURE 7 – The four components of modular and iterative aggregation algorithm for GNSS trajectories.</figcaption></figure>

Figure 7 illustrates the proposed Modular and Iterative Aggregation Algorithm (MIAA) for GNSS trajectories which is composed by four components with new options. Thus, this proposal will create many variants of the algorithm. Note that, each component corresponds to a step of the algorithm.

The open-source implementation of the MIAA algorithm have been integrated into the Tracklib library ([MvD24]).

## 4.1 Component 1 - choosing the master trajectory

This component is representing the first step of the workflow. We propose three options (Cf. Figure 8) to get the master trajectory. The first, already proposed in [ED14], selects as master, the trajectory whose length is closest to the median of the lengths of all trajectories to be aggregated. We add, to this initial option, two new options : (1) the master trajectory is one that minimizes the sum of distances to other trajectories and (2) the master trajectory is randomly selected from the set of trajectories to be aggregated. These features have the advantages to study the influence of the choice of the master trajectory on the final result.

<figure><img src="image_14.png" /><figcaption>FIGURE 8 – Heuristics for choosing a trajectory *master*</figcaption></figure>
## 4.2 Component 2 - matching trajectories with the master trajectory

This component is linked to the second step of the algorithms and contains different measures to compute the distance between two trajectories. Among the many similarity measures, we considered four measures : in addition to the matching based on the discrete Frechet distance, we add as new variants, two methods parameterized with $L_p$-norm ($p \in 1, 2$) and the nearest neighbour matching.

<figure><img src="image_15.png" /><figcaption>FIGURE 9 – Heuristics for choosing a similarity measure for matching trajectories</figcaption></figure>

Clearly, the Nearest Neighbour distance differs from the others because it does not preserve the order of position in matching between the two trajectories, as described in the Figure 9 (bottom right image).

As mentioned in the section « Formalization of the original algorithm », Discrete Fréchet distance seeks to minimize the maximum of distances between the two curves. Once the maximum distance link is fixed, the other matching links are placed in the remaining spots. In Figure 9, the matching links are either grouped for distant points (Circled area 4) or create a large offset (Circled area 3).

DTW, both $L_1$ and $L_2$, aims at minimizing the sum of these distances, the matching links are therefore spaced more evenly, as shown in the top 2 images in the Figure 9. The difference between the $L_1$ and $L_2$ distances in matching is minimal. The $L_2$ distance is the straight-line distance, while the $L_1$ distance follows the x-axis and y-axis. In Figure 9, the difference between these two distances becomes significant when the links between the points are close to the diagonals of the two axis (Circled area 1 and circle area 2).
## 4.3 Component 3 - choosing the representative position of each homologous points

As mentioned in the second step of the iteration of the aggregation algorithm, it is possible that multiple points on a trajectory are matched to the same point in the reference trajectory. Then it is required to reduce these points to a single representative position. This can be done with different methods : center of gravity, position with median time, position furthest from the master trajectory.

<figure><img src="image_16.png" /><figcaption>FIGURE 10 – Heuristics for options in [C3] for DTW $L_\infty$</figcaption></figure>

With the Discrete Fréchet distance ($L_\infty$) in C2, Figure 10 describe the three options. Unsurprisingly, the higher the link cardinality, the more this option will influence the position of the representative. The choices of the center of gravity and the position furthest from the master trajectory (Circled area 2) correspond to a spatial strategy, unlike the option of the position with median time. In Figure 10, since the points of the track are spaced at fixed time intervals, the difference between option 1 (Circled area 1) and option 3 (Circled area 3) is the same.

## 4.4 Component 4 - aggregating the representative position

So, for each reference position $(\mathbf{r}_j)_{j=1..m}$ on the master trajectory $\mathcal{R}$, we have a representative $\bar{x}_j$ on each trajectory $\mathcal{R}$. Using an aggregation operator, a unique aggregating position can represent each reference position $\mathbf{r}_j$. We propose four aggregate function to calculate aggregating position : the marginal median ([ED14]), the geometric median, the $L_2$ mean and the $L_\infty$ which is the center of the minimum covering circle. In this component, we can add a constraint with the four approaches : whether or not to anchor the aggregated position to an existing position in the dataset's trajectories as explained in step 3 of 3.

<figure><img src="image_17.png" /><figcaption>FIGURE 11 – Heuristics for options in [C4] for DTW $L_\infty$</figcaption></figure>

This step is the one that determines the final position of the points in the track. In Figure 11, we can observe that the marginal median suggests a point among the positions of the representatives of C3, whereas the geometric mean proposes a new position for the trajectory, which can therefore be influenced by outliers.
# 5 Extensive analysis and calibration

## 5.1 Termination of the algorithm

It should be mentioned that, to our knowledge, the termination of the algorithm has not yet been studied in [ED14]. MIAA algorithm is an iterative algorithm, and there is no guarantee that it will converge from an IT point of view. For example, at the end of an iteration, you can return to a previous step (the merged trajectory corresponds to a previous master trajectory), therefore iterations are entering a cycle. We demonstrate this potential drawback with a counter-example described below.

For example, consider these two trajectories shown in Figure 12-a :

$$X = < (68.0, 20.0), (69.0, 22.0), (69.0, 24.0), (67.0, 25.0) >$$

$$Y = < (71.0, 14.0), (71.0, 16.0), (72.0, 19.0), (70.0, 22.0) >$$

<figure><img src="image_18.png" /><figcaption>FIGURE 12 – Counter-example showing a case of non-convergence</figcaption></figure>

The selected configuration is : for the master trajectory, is one that minimizes the sum of distances to other trajectories, the DTW-L1 has been chosen for matching trajectories, the center of gravity is taken for the representative position of each homologous and to aggregate the representative position, median aggregator is used to ultimately be map matched over a position of the trajectories.

In this configuration, the first iteration of the MIAA algorithm returns the aggregated trajectory represented by the red line in Figure 12-b. Now let us take a closer look at the second iteration. The matching of the trajectories with the master trajectory produces a single matching link (one vertex of $X_2$ with one vertex of $R$) except for the first reference position, we have a matching link 1-3 : ($r_1$, $x2_1$, $x2_2$, $x2_3$) (Figure 12-c) (corresponds to these three points of the second trajectory $x_1$ (71.0,14.0), $x2_2$ (71.0,16.0) and $x2_3$ (72.0,19.0)). Consequently the center of gravity of these three points is $\overline{x2_1}$ (214/3, 49/3) that is to say, the representative position associated to $r_1$ (Figure 12-d).

With anchor constraint, we need to match the center of gravity $\overline{x2_1}$ with one of the four vertex : $r_1$, $x2_1$, $x2_2$, $x2_3$, respecting this order (IT implementation). There are two candidates positions with the same distance $r_1$ and $x2_3$ (6.1388...8), so the winner is $r_1$ (first in the list).

At the end of the second iteration, the aggregated trajectory corresponds to the master trajectory of the first iteration. We are locked in a two-iteration cycle.
## 5.2 Algorithm calibration : further details

This section reproduces the online notebook on the Tracklib documentation website "Aggregated Trajectory : position errors and shape deviation estimation". It is inspired by the work presented in [VDMOR24] and also provides some pieces of Python code from Tracklib library ([MvD24]) documentation<fnref n="3" />.

We will examine the aggregation of trajectories using two matching distances : the Fréchet distance and the DTW-L2 distance, and thus see its ability to reconstruct accurately the common path followed by all the individual sample trajectories by comparing the position errors and the shape deviation.

Our experiment will be conducted in four steps :

1. Step 1 : creating a synthetically reference track considered as the ground truth track,

2. Step 2 : creating a set of simulated tracks from the reference track,

3. Step 3 : computing the aggregation track from the set previously constructed (step 2),

4. Step 4 : the error between the estimated and ground truth track is then evaluated.

### 5.2.1 Step 1 : Create reference tracks

For this experimentation, we will examine the trajectory aggregation process for different types of paths to be reconstructed. In the context of mountain hiking, we identified three characteristic path shapes : nearly straight segments, moderately sinuosity road segments, and a zigzagging path composed of a series of switch-backs. To maintain shape consistency throughout the entire road segment, the generated trajectories have a length of approximately 300 meters.

This three commonly mountain path shapes correspond to the reference tracks.

```python
tkl.seed(123)
# ----------------------------------------------------------------------
# Generate the path 'Almost straight'
sentier1 = tkl(generate(0.5, dt=10)[::3]
sentier1.scale(5)

# ----------------------------------------------------------------------
# Generate the path 'Moderate sinuosity'
sentier2 = tkl(generate(0.1, dt=10)[::3]

# ----------------------------------------------------------------------
# Generate the path 'Switchbacks'
base_lacets = tkl(generate(0.4, dt=10)
sentier3 = tkl.noise(base_lacets, 20, tkl.SincKernel(20),
direction=tkl.MODE_DIRECTION_ORTHO)[::3]
sentier3.scale(4)

# ----------------------------------------------------------------------
SHAPES = ['Almost straight', 'Moderate sinuosity', 'Switchbacks']
sentiers = [sentier1, sentier2, sentier3]
```

<figure><img src="image_19.png" /></figure>

<footnote n="3">3. https://tracklib.readthedocs.io/en/latest/usecase/AggregatedTrajectory.html</footnote>
### 5.2.2 Step 2: Create sets of simulated GNSS trajectories

Now that we have our three reference trajectories, we want to create three sets of N simulated GNSS trajectories from them. To generate realistic noise, we used an approach described in [Rip09] and also employed in [MMALGOR23]:

with a random generator, we sampled N i.i.d. unit-variance and zero-mean gaussian values, compiled in a vector x. It can easily be shown that, for any positive-definite matrix $Σ \in R^{n \times n}$, the random vector $y = Ax$ where A is a Cholesky factor of Σ, is a realization of a correlated random vector Y having covariance matrix Σ. The covariance matrix Σ is formed with a (stationary) covariance kernel with three parameters:

- The **type of kernel**: exponential, gaussian, and triangular models are used.

- The **amplitude of noise**: is between 0 and 5 meters, as it is quite uncommon to find building databases with more than 5 m error amplitude. If necessary, the output tables could be extended to handle large errors.

- The correlation **scope** of the noise which roughly speaking describes how far apart two errors would remain correlated (in both amplitude and direction): between 1 m (white noise) and 1000 m (global translation).

In tracklib, you have to create a Kernel and use the noise method on a track.

Trajectories have been generated with a 5 m-amplitude exponential covariance process, completed by a 1 m white noise process, a 50 cm range Gaussian Process for reference error.

```python
N = 20

# Generate 'Almost straight'
tracks1 = tkl.core.TrackCollection([sentier1]*N)
tracks1.noise(5, tkl.ExponentialKernel(1))
tracks1.noise(1, tkl.GaussianKernel(0.5))

# Generate 'Moderate sinuosity'
tracks2 = tkl.core.TrackCollection([sentier2]*N)
tracks2.noise(5, tkl.ExponentialKernel(1))
tracks2.noise(1, tkl.GaussianKernel(0.5))

# Generate 'Switchbacks'
tracks3 = tkl.core.TrackCollection([sentier3]*N)
tracks3.noise(5, tkl.ExponentialKernel(1))
tracks3.noise(1, tkl.GaussianKernel(0.5))
```

<figure><img src="image_20.png" /></figure>

### 5.2.3 Step 3: Compute aggregated trajectory

We have two evaluations to do: compute the aggregated trajectory with Frechet distance and with $L_2$ distance.

Each one is conducted as follows: we generate $N'$ random noisy traces, with $N' \in [1, E]$ and $E \le N$. For each generation, we are going to compare the aggregated trajectory against the reference track (*i.e.* the ground truth).
This code compute an aggregation track for a set of 10 trajectories :

FRECHET = {'Almost straight': tkl.TrackCollection(),
           'Moderate sinuosity': tkl.TrackCollection(), 'Switchbacks': tkl.TrackCollection()}

DTW = {'Almost straight': tkl.TrackCollection(),
      'Moderate sinuosity': tkl.TrackCollection(), 'Switchbacks': tkl.TrackCollection()}

represent_method = tkl.MODE_REP_BARYCENTRE

aggmeth = tkl.MODE_AGG_MEDIAN

cstt = False

master = tkl.MODE_MASTER_MEDIAN_LENGTH

itermax = 25

# create set of 10 trajectories

TAB = set()

while len(TAB) <= 10:
    n = random(0, N-1)
    TAB.add(n)

sets1 = tkl.TrackCollection()

for idx in TAB:
    sets1.addTrack(tracks1[idx])

# Frechet

p = float('inf')

mode = tkl.MODE_MATCHING_FRECHET

central1 = tkl.fusion(sets1, master=master, dim=2, mode=mode, p=p,
                          represent_method=represent_method, agg_method=aggmeth, constraint=cstt)

FRECHET['Almost straight'].addTrack(central1)

...

# DTW

p = 2

mode = tkl.MODE_MATCHING_DTW

central1 = tkl.fusion(sets1, master=master, dim=2, mode=mode, p=p,
                          represent_method=represent_method, agg_method=aggmeth, constraint=cstt)

DTW['Almost straight'].addTrack(central1)

...

<figure><img src="image_21.png" /><figcaption>FIGURE 13 – Aggregated trajectories for two samples of 3 and 20 trajectories and computed with Discrete Fréchet distance (in pink) and DTW-L₂ distance (in blue)</figcaption></figure>
## Shape deviation measurement

Let's start by aligning, with a geometric affine transformation, the aggregated track with the reference track to abstract away from position issues. Then, we finely resample both the aggregated and ground truth trajectories from 120 to 1,000 points, ensuring one point every 50 cm over 500 meters. Finally, we compute the distance *NearestNeighbour* for all positions to estimate shape deviation.

```python
def shapeDeviationMeasure(central, sentier):
    # Align the aggregated track with the reference track
tkl mapping.mapOn(central, sentier, verbose=False)

    # resample to 1000 points
central.resample(npts=1000, mode=1)
sentier.resample(npts=1000, mode=1)

    # compute the distance NearestNeighbour
return tkl compare(central, sentier, tkl.MODE_COMPARISON_NN, p=2)

shapeDeviationMeasure(FRECHET['Almost straight'][s], sentier1)
...
shapeDeviationMeasure(DTW['Switchbacks'][s], sentier3)
```

Which results in the following graphically :

<figure><img src="image_22.png" /><figcaption>FIGURE 15 – Shape deviation measurement for each samples of *N* trajectories generated</figcaption></figure>

The Figure 15 still shows a close similarity between the two curves : shape deviation measurement are quite close together and their profiles are quite similar.
5.3 Crowd-sourced trajectories dataset

Here we present the results obtained on real data. As already mentioned, the analyse is not strictly "metro-
logical" but rather qualitative, a graphic investigation. We define a first context : a ridge which characterises
terrain constraints (Figure 16a), a switchback trail which characterises an infrastructure support (Figure 16b)
and a straight line followed by a series of twists and turns with varying distances between them which cha-
racterises variations of the shape of the path (Figure 16c).

Note that, traces may have been previously filtered via a simplification algorithm like Douglas-Peucker.
Indeed, all geometries offer the same characteristics like vertex of bends, turns in route. If this was not the
case, traces would be more dissimilar in shape.

<figure><img src="image_23.png" /><figcaption>(a) Aggregated trajectory closely follows the ridge path</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) Aggregated trajectory in red follow the most popular route</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c) median operation effect in a series of bends with a
Crowdsourced dataset</figcaption></figure>

<caption>FIGURE 16 – Aggregated trajectory result on Crowd-sourced dataset</caption>

Aggregated trajectory produced by the MIAA algorithm (red line in Figure 16b and white line in Figure 16c)
detects the road shape very well, all the bends have been detected. In the context of mountain hiking, we can
see the aggregated trajectory deviates from the road, which may make it possible to quantify the off-tracking.
In Figure 16a, we can observe that some turns are not detected, as the geometry of the aggregated trajectory
tends to generalize and lose certain inflection points.
5.4 Perspective : performance under extreme conditions

Concerning high noise amplitude, we deliberately chose a low one (5 meters) in our experiments to adapt to advancements in GNSS technological (new constellations, new signals, improvements in receiver electronics). However, we have not studied the case of extreme noise conditions, and we could relaunch computations with extreme values, on the same synthetic GNSS trajectories datasets, with an amplitude of the exponential covariance process from 5 meters (Cf. [LBMC23]) to 30 meters (as mentioned in the literature for errors in heavy coverage canopy areas in [PG05]).

Additionally, for outlier traces, particularly off-roads, we conducted empirical tests on the multi-sensors/multi-canopy traces dataset. With few outliers, our tests confirmed there was no impact on the aggregated trajectory, demonstrating its robustness as mentioned in [ED14]. The algorithm computes, for each matched points in the master trajectory, the robust median positions between all trajectories. Thus, we decided not to further investigate robustness, in order to focus on its modularity. We could study the performance of the algorithm with many off-roads in the sample traces to be merged. For example, the robustness analysis would consist in including a proportion p of off-roads in a given number of trajectories, and determine how the quality of the aggregated trajectory decreases as p increases. In our opinion, with the effect of the median, results should not be dramatically modified, at least up to p=25%. If needed, we can run the test and add comments with a figure depicting the aggregated trajectory quality as a function of p.
Références

[BHLS15] Filip Biljecki, Gerard BM Heuvelink, Hugo Ledoux, and Jantien Stoter. Propagation of positional error in 3d gis : estimation of the solar irradiation of building roofs. *International Journal of Geographical Information Science*, 29(12) :2269–2294, 2015.

[BKT11] Henrik Blunck, Mikkel Baun Kjærgaard, and Thomas Skjødeberg Toftegaard. Sensing and classifying impairments of gps reception on mobile devices. In Kent Lyons, Jeffrey Hightower, and Elaine M. Huang, editors, *Pervasive Computing*, pages 350–367, Berlin, Heidelberg, 2011. Springer Berlin Heidelberg.

[Bon02] Olivier Bonin. Modèle d'erreurs dans une base de données géographiques et grandes déviations pour des sommes pondérées : application à l'estimation d'erreurs sur un temps de parcours. PhD thesis, Paris 6, 2002.

[CLP24] P. Calloch, B. Labbé, and L. Peseux. Établissement d'une méthode low-cost de vérité terrain pour l'étude de la précision de traces gnss de randonneurs en forêt. Technical report, École nationale des sciences géographiques (ENSG-Géomatique), 2024.

[ED14] Laurent Etienne and Thomas Devogele. Trajectoires médianes. In *14 ème conférence Extraction et Gestion des Connaissances, Ateliers fouille de données spatiales et temporelles & construction, enrichissement et exploitation de ressources géographiques pour l'analyse de données*, Rennes, France, January 2014. URL: https://hal.science/hal-02110086.

[GBTY05] Dorota Grejner-Brzezinska, Charles Toth, and Yudan Yi. On improving navigation accuracy of gps/ins systems. *Photogrammetric engineering & remote sensing*, 71(4) :377–389, 2005.

[GPML22] Léna Gruas, Clémence Perrin-Malterre, and Anne Loison. From the Crowded Valleys to the Preserved Summits : Mountain Sports Participants' Attitudes Toward Protected Areas in the Sprawling Urban Areas of the Northern French Alps. *Mountain Research and Development*, 42(3) :R50, 2022. doi:10.1659/MRD-JOURNAL-D-21-00001.1.

[LBMC23] Taeyoon Lee, Pete Bettinger, Krista Merry, and Chris Cieszewski. The effects of nearby trees on the positional accuracy of gnss receivers in a forest environment. *PLOS ONE*, 18(3) :1–20, 03 2023. doi:10.1371/journal.pone.0283090.

[MMALGOR23] Y. Méneroux, I. Maidaneh Abdi, A. Le Guilcher, and A.-M. Olteanu-Raimond. Is the radial distance really a distance? an analysis of its properties and interest for the matching of polygon features. *International Journal of Geographical Information Science*, 37(2) :438–475, 2023. doi:10.1080/13658816.2022.2123487.

[MvD24] Yann Méneroux and Marie-Dominique van Damme. Tracklib : a python library with a variety of tools, operators and functions to manipulate GPS trajectories, December 2024. URL: https://hal.science/hal-04356178.

[PG05] Christian Piedallu and Jean-Claude Gégout. Effects of forest environment and survey protocol on gps accuracy. *Photogrammetric Engineering & Remote Sensing*, 71(9) :1071–1078, 2005. URL: https://www.ingentaconnect.com/content/asprs/pers/2005/00000071/00000009/art00006, doi:10.14358/PERS.71.9.1071.

[Rip09] Brian D Ripley. Stochastic simulation. John Wiley & Sons, 2009.

[RMC+25] Ana-Maria Raimond, Yann Meneroux, Paul Calloch, Brice Labbé, Lorine Peseux, Franck Tertre, and Emmanuel Clédat. DATASET : Multi-sensors and Ground Truth Tracks, 2025. doi:10.57745/MOUSLS.

[Rob93] William David Summerfield Roberts. GPS time correlation and its implication for precise navigation. PhD thesis, Newcastle University, 1993.
[Vau97] François Vauglin. *Modèles statistiques des imprécisions géométriques des objets géographiques linéaires*. PhD thesis, Université de Marne-la-Vallée (1991-2019), 1997.

[VDMOR24] Marie-Dominique Van Damme, Yann Méneroux, and Ana-Maria Olteanu-Raimond. A metrological analysis of a modular and iterative aggregation algorithm of gnss trajectories. In *Proceedings of the 32nd ACM International Conference on Advances in Geographic Information Systems*, SIGSPATIAL '24, page 633–636, New York, NY, USA, 2024. Association for Computing Machinery. doi:10.1145/3678717.3691325.

[ZY15] Yongli Zhang and Yuhong Yang. Cross-validation for selecting a model selection procedure. *Journal of Econometrics*, 187(1):95–112, 2015.