Costas Panagiotakis

Nikos Pelekis

Ioannis Kopanakis

Emmanuel Ramasso

Yannis Theodoridis

Segmentation and Sampling of Moving Object Trajectories based on Representativeness

Keywords: Trajectory Segmentation, Subtrajectory sampling, Data mining, Moving Object Databases ✦

Moving Object Databases (MOD), although ubiquitous, still call for methods that will be able to understand, search, analyze, and browse their spatiotemporal content. In this paper, we propose a method for trajectory segmentation and sampling based on the representativeness of the (sub-)trajectories in the MOD. In order to find the most representative sub-trajectories, the following methodology is proposed. First, a novel global voting algorithm is performed, based on local density and trajectory similarity information. This method is applied for each segment of the trajectory, forming a local trajectory descriptor that represents line segment representativeness. The sequence of this descriptor over a trajectory gives the voting signal of the trajectory, where high values correspond to the most representative parts. Then, a novel segmentation algorithm is applied on this signal that automatically estimates the number of partitions and the partition borders, identifying homogenous partitions concerning their representativeness. Finally, a sampling method over the resulting segments yields the most representative sub-trajectories in the MOD. Our experimental results in synthetic and real MOD verify the effectiveness of the proposed scheme, also in comparison with other sampling techniques.

INTRODUCTION

N OWADAYS, there is a tremendous increase of Mov- ing Objects Databases (MOD) [START_REF] Guting | Moving Object Databases[END_REF] due to, on the one hand, location-acquisition technologies like GPS and GSM networks and, on the other hand, computer vision based tracking techniques. This explosion of information combines an increasing interest in the area of trajectory data mining and, more generally, knowledge discovery from movement-aware data [START_REF] Giannotti | Mobility, Data Mining and Privacy, Geographic Knowledge Discovery[END_REF]. All these technological achievements require new services, software methods, and tools for understanding, searching, retrieving, and browsing spatiotemporal trajectories content.

In this paper, we tackle a problem combining three different aspects. First of all, we study the problem of alternative representations of trajectories of moving objects (other than the usual sequences of 3-D line segments), according to contextual information that can be automatically derived by the total trajectory population. More specifically, we investigate for an effective way to represent each trajectory by a continuous function that implicitly describes the "representativeness" of each constituent part of it (i.e. a segment) w.r.t. the whole MOD. Given such an intuitive representation, a second interesting arising problem is that of its segmentation in a way that an analyst could gain insight into "representative" (i.e. interesting, dense, frequent) portions (i.e. subtrajectories), but also into "non-representative" parts, which are also of interest in various application scenarios (for example, in detecting movement outliers). On top of the previous issues, and due to the complex nature of the trajectory data and the vast volumes of MOD, a third interesting problem arises; that of "trajectory sampling". This is a very challenging problem where very limited work has been carried out so far. An insightful solution to the problem would be an analyst to be able to supervise the sampling procedure, not only regarding the volume of the sampled dataset, but also the properties of the dataset that reveal the underlying movement patterns of the MOD. In this paper, we argue that this problem can be effectively tackled if interconnected to the previous two discussed problems. In other words, we propose an automatic method for sub-trajectory sampling based on the "representativeness" of the sub-trajectories. In this approach, an analyst may request the top-k representative sub-trajectories that best describe the MOD in an optimised way, where optimization is with respect to the "representativeness". Fig. 1 illustrates an example of a MOD comprised by four trajectories ({T 1 , ..., T 4 }) and the top-2 representative sub-trajectories ({S 1 , S 2 }) that best describe the MOD.

Recently, in the literature there have been proposed several works that try to either efficiently analyze trajectory data or mine movement-aware patterns. In the domain of trajectory segmentation related works deal with the problem locally, by partitioning trajectories in a way as to achieve better database organization [START_REF] Hadjieleftheriou | Efficient indexing of spatiotemporal objects[END_REF] or extract more intuitive local patterns for clustering [START_REF] Han | Trajectory clustering: a partition-and-group framework[END_REF] and classification purposes [START_REF] Lee | Traclass: trajectory classification using hierarchical region-based and trajectory-based clustering[END_REF]. In [START_REF] Anagnostopoulos | Global distance-based segmentation of trajectories[END_REF], a global distance-based approach was proposed for the segmentation of object trajectories that attempts to simplify the trajectories by using minimum bounding rectangles (MBRs) in a way that the original pairwise distances between trajectories are minimized, as such facilitating clustering and classification tasks. Furthermore, there are approaches that attempt to extract global patterns by clustering trajectories [START_REF] Chen | Robust and fast similarity search for moving object trajectories[END_REF]- [START_REF] Pelekis | Similarity search in trajectory databases[END_REF] based on different types of distance functions, while others try to take advantage of the partial (local) grouping of the trajectories [START_REF] Han | Trajectory clustering: a partition-and-group framework[END_REF], [START_REF] Benkert | Reporting flock patterns[END_REF], [START_REF] Li | Clustering moving objects[END_REF]. There are also hybrid approaches that succeed global clustering by using local criteria and focusing in symbolic representations of the trajectories [START_REF] Pelekis | Clustering trajectories of moving objects in an uncertain world[END_REF]. In the domain of sampling, although there are nice solutions for point data [START_REF] Kollios | Efficient biased sampling for approximate clustering and outlier detection in large datasets[END_REF], [START_REF] Nanopoulos | Indexedbased density biased sampling for clustering applications[END_REF], to the best of our knowledge, trajectoryoriented solutions include explorative, user-supervised, clustering-based sampling techniques [START_REF] Andrienko | A visual analytics toolkit for cluster-based classification of mobility data[END_REF], [START_REF] Andrienko | Interactive visual clustering of large collections of trajectories[END_REF] as well as approaches that operate on approximate trajectories [START_REF] Pelekis | Unsupervised trajectory sampling[END_REF], with all of those solutions trying to select trajectories as wholes without being able to identify representative sub-trajectories, a problem completely different, which is exactly the challenge accepted by the current work.

Most of the above mentioned approaches propose different similarity metrics which they are used either for introducing indexing structures for efficient trajectory retrieval [START_REF] Chen | Robust and fast similarity search for moving object trajectories[END_REF], or for clustering purposes [START_REF] Nanni | Time-focused clustering of trajectories of moving objects[END_REF], [START_REF] Pelekis | Similarity search in trajectory databases[END_REF]. Some techniques simplify the given trajectories, focus on spatial criteria and ignore the temporal dimension [START_REF] Han | Trajectory clustering: a partition-and-group framework[END_REF], [START_REF] Lee | Traclass: trajectory classification using hierarchical region-based and trajectory-based clustering[END_REF]. Other proposals either cannot capture the complex nature of the trajectories that follow arbitrary motion patterns [START_REF] Han | Trajectory clustering: a partition-and-group framework[END_REF], or focus on the discovery of very specific definitions of movement patterns [START_REF] Benkert | Reporting flock patterns[END_REF], [START_REF] Li | Clustering moving objects[END_REF], [START_REF] Jeung | Discovery of convoys in trajectory databases[END_REF]. Explorative, usersupervised approaches [START_REF] Andrienko | A visual analytics toolkit for cluster-based classification of mobility data[END_REF], [START_REF] Andrienko | Interactive visual clustering of large collections of trajectories[END_REF] do not provide deterministic solutions as they base on clustering techniques, which in their turn rely on the selection of the clustering algorithm, its parameters and the intrinsic distance function used. On the contrary, in this paper, we do not simplify the MOD by any time-consuming preprocessing steps, we take into account the temporal information, while we do not make any assumption regarding the type of pattern, that is the actual final outcome. Given this noprerequisites setting, we argue that all of the above approaches, as well as those which are dealing with vast volumes of trajectory datasets would benefit if they would be applied in a representative subset (consisting of the representative sub-trajectories) that best describes the whole dataset. In other words, we seek for a methodology to sample those sub-trajectories from a MOD that preserve as much as possible the properties and the mobility patterns hidden in the original MOD. Consider, for example, the domain of visual analytics on movement data [START_REF] Andrienko | Visual analytics tools for analysis of movement data[END_REF] in which it is meaningless to visualize datasets larger than a certain small size, due to the distinguishing properties of the human eye. How would one select a subset so that this would cover the whole space-time and, as such, get the gist of the whole dataset is the challenge set in this paper.

The contributions and merits of our work are summarized below:

1) We propose an index-based global voting method that allows us to represent the representativeness of a trajectory in a MOD as a smooth continuous descriptor. 2) We introduce an algorithm for the automatic segmentation of trajectories into "homogenous" subtrajectories according to their "representativeness" in the MOD. 3) We define the problem of sub-trajectory sampling in a MOD as an optimization problem and we propose a novel solution to tackle the problem. 4) Finally, we conduct a comprehensive set of experiments over synthetic and real trajectory datasets, in order to thoroughly evaluate our approach.

The rest of the paper is organized as follows: Section 2 presents previous work in related domains, while Section 3 sets the scene by presenting the various aspects of the problem and also describing the base of our further developments. Section 4 presents the proposed threestep methodology, i.e. our trajectory voting scheme, the trajectory segmentation technique, and the sub-trajectory sampling method. The experimental results are given in Section 5. Finally, conclusions and discussion are provided in Section 6.

RELATED WORK

In this section, we review existing works in the domains related with the current work. In our setting, representative (sub)-trajectories are a new type of mobility pattern, as such, our discussion includes trajectory pattern mining, segmentation and sampling in MOD.

A MOD consists of spatiotemporal trajectories of moving objects (e.g. humans, vehicles, animals, etc.). In the general case, trajectories are represented as 3-D sequences where each recording encodes the 2-D (two dimensional) geographic location and the 1-D temporal information of mobile objects. During the last decade several approaches have been proposed in the literature so as to enable well-known mining algorithms to operate on trajectories. One such approach is the use of different types of distance functions as the mean to group trajectories into clusters. Some approaches are inspired by the time series analysis domain [START_REF] Vlachos | Rotation invariant distance measures for trajectories[END_REF], [START_REF] Chen | Robust and fast similarity search for moving object trajectories[END_REF], while other exploit on a set of distance operators based on primitive (space and time) as well as derived parameters of trajectories (speed and direction) [START_REF] Pelekis | Similarity search in trajectory databases[END_REF]. An interesting approach also used in our approach is proposed in [START_REF] Frentzos | Index-based most similar trajectory search[END_REF] for the efficient processing of most-similar trajectory (MST) queries. A similar distance function is used in [START_REF] Nanni | Time-focused clustering of trajectories of moving objects[END_REF], where Nanni and Pedreschi adapt the well-known density-based OPTICS [START_REF] Ankerst | Optics: Ordering points to identify the clustering structure[END_REF] clustering algorithm, tailored to work with point data, to a new algorithm for trajectory data, named T-OPTICS (and its varriant TF-OPTICS, which focuses on the discovery of the temporal intervals that lead to best clustering results). The previously mentioned temporal intervals are given by the user, so TF-OPTICS essentially re-executes T-OPTICS on segments of trajectories, obtained by properly clipping the original ones. In comparison with our approach, we automatically segment trajectories to portions based on global criteria (i.e. the representativeness of the trajectory in the MOD). Furthermore, TF-OPTICS mainly clusters whole trajectories and is not tailored to identify patterns of sub-trajectories in an un-supervised way.

Recently, Pelekis et al. [START_REF] Pelekis | Clustering trajectories of moving objects in an uncertain world[END_REF] proposed another approach, called CenTR-I-FCM, taking into advantage of local patterns in time dimension as the base to identify global clusters of whole approximate/symbolic trajectories. In comparison with the current work, this approach also utilizes a global but static and predefined temporal segmentation of trajectories. In addition, in this work trajectories are symbolically represented as intuitionistic fuzzy vectors and not as sequences of 3-D line segments. This approach also aims at clustering trajectories as a whole with special care for handling uncertainty.

Authors in [START_REF] Han | Trajectory clustering: a partition-and-group framework[END_REF] proposed TRACLUS, a partition-andgroup framework for clustering 2-D trajectories which enables the discovery of common sub-trajectories, based on a trajectory partitioning algorithm that uses the minimum description length principle. The core of the framework uses a variant of the DBSCAN algorithm that operates on the partitioned directed line segments. Finally, the notion of the representative trajectory of a cluster is defined. The fundamental differences of TR-ACLUS with our approach is that, a) the temporal information is not considered in [START_REF] Han | Trajectory clustering: a partition-and-group framework[END_REF], b) the segmentation is performed per trajectory and it does not use global criteria, c) the identified clusters of segments conform to straight movement patterns and cannot identify complex (e.g. snake-like) patterns, which are usual in real world applications, and d) the proposed algorithm for identifying the representative trajectory is defined per cluster and it is a synthetic trajectory computed by an averaging technique. Instead, we automatically select the top-k most representative sub-trajectories from the given dataset.

Other related works include techniques for flock patterns [START_REF] Benkert | Reporting flock patterns[END_REF], moving clusters [START_REF] Li | Clustering moving objects[END_REF] and convoys [START_REF] Jeung | Discovery of convoys in trajectory databases[END_REF]. Although these approaches provide lucid definitions of the mined patterns, they are rather rigorous and sensitive to parameters, while they are not in the context of automatic trajectory segmentation according to the representativeness of the trajectories, as well as sub-trajectory sampling.

In [START_REF] Sacharidis | On-line discovery of hot motion paths[END_REF], hot motion paths (frequently traveled trails of numerous moving objects) are detected in a distributed system under the assumption that the moving objects can communicate with a central unit (coordinator). The goal of this work is in the same context with only one of the aims of our research (i.e. the one concerning trajectory representativeness). However, this approach focuses on the online processing and maintenance of hot motion paths by imposing a sliding time window of size W , which excludes from consideration any locations received more than W time units ago.

In [START_REF] Anagnostopoulos | Global distance-based segmentation of trajectories[END_REF], distance-based criteria have been proposed for the segmentation of trajectories. The distance between two trajectories is defined using their MBR representation. The segmentation problem is designed as a maximization problem, that attempts to create MBRs in such a way that the original pairwise distances between trajectories are minimized. We also provide a method for trajectory segmentation having global criteria but this splitting is entailed by the "representativeness" of the trajectory, which is not simply defined in terms of distance but via a novel voting method.

In [START_REF] Panagiotakis | Trajectory voting and classification based on spatiotemporal similarity in moving object databases[END_REF], we have proposed an approach for expressing the "representativeness" in MOD via a voting process that is applied for each segment of a given trajectory. Moreover, a simplistic trajectory ranking method selects the trajectories of highest voting and as such ignores trajectories in low density regions of the MOD. In the current work, we improve the voting scheme originally proposed in [START_REF] Panagiotakis | Trajectory voting and classification based on spatiotemporal similarity in moving object databases[END_REF] with the use of a metric distance function that further allows us to propose an index-based algorithm for the efficient implementation of the voting scheme. More importantly, we improve [START_REF] Panagiotakis | Trajectory voting and classification based on spatiotemporal similarity in moving object databases[END_REF] by handling trajectory segmentation and sub-trajectory sampling aspects.

Trajectory sampling is a new topic in the spatiotemporal database literature with a variety of applications including MOD summarization, visualization, searching and retrieval. Although very interesting as a problem and with great potentials in the MOD domain, to the best of our knowledge there is no related work that tackles the problem in a tailored way to the complex nature of trajectory datasets. In [START_REF] Kollios | Efficient biased sampling for approximate clustering and outlier detection in large datasets[END_REF], [START_REF] Nanopoulos | Indexedbased density biased sampling for clustering applications[END_REF] the authors propose techniques for density biased sampling (DBS) in point sets that obviously can not be applied to trajectory data. In our approach, we extend the idea of DBS in a way that density properties as well as the similarity of trajectories segments is taken into account in a way that allows us to effectively apply sampling in the MOD domain. In [START_REF] Pelekis | Unsupervised trajectory sampling[END_REF] we have proposed an approach for unservised trajectory sampling, however this work operates on approximate, symbolic trajectories (i.e. represented as sequences of cells of a spatio-temporal grid wherefrom the moving object pass), while it focuses on the selection of whole trajectories and not portions of them (i.e. sub-trajectories produced by a segmentation method). These two observations make the problem addressed in the current work completely different.

An interesting approach focusing on the visualization of large trajectories datasets and which uses standard sampling techniques has been proposed in [START_REF] Andrienko | A visual analytics toolkit for cluster-based classification of mobility data[END_REF], [START_REF] Andrienko | Interactive visual clustering of large collections of trajectories[END_REF]. The authors use uniform sampling as a starting step to minimize the volume of the trajectories that can be clustered by the T-OPTICS algorithm [START_REF] Nanni | Time-focused clustering of trajectories of moving objects[END_REF]. At a subsequent step, they either select a certain percent of trajectories from a cluster using distance-based thresholds or they apply a partitioning on each cluster so as to identify spherical sub-clusters, each of which is represented by its medoid. Obviously, the outcome of this approach could be reached by directly applying K-Means or K-Medoid algorithm and then choosing the mean or medoid of each cluster as the representatives that will be sampled. Of course, this choice has the disadvantage that the number of clusters should be known apriori. This approach, which is basically a stratified sampling technique, where strata are produced by the clustering algorithm, has the limitation that it is user-supervised and it depends on the results of the clustering. Notwithstanding our methodology is automated, which is a self-evident advantage, in our experimental study we show the superiority of our approach in comparison to uniform random and stratified sampling techniques used in [START_REF] Andrienko | A visual analytics toolkit for cluster-based classification of mobility data[END_REF], [START_REF] Andrienko | Interactive visual clustering of large collections of trajectories[END_REF], by using standard sampling evaluation metrics.

SETTING THE SCENE

In this section, we set the scene of the various aspects of the problem that this paper addresses, and concurrently we present the stepping-stones where our subsequent developments base on. Let us assume a MOD D = {T 1 , T 2 , • • • , T N } of N trajectories, where T k denotes the k-th trajectory of the dataset, k ∈ {1, 2, ..., N }. We assume that the objects are moving in the xy plane. Let p k (i) = (x k (i), y k (i), t k (i)) be the i-th point, i ∈ {1, 2, ..., L k }, of k-th trajectory, where L k denotes the number of points of k-th trajectory. x k (i), y k (i) and t k (i) denote the 2-D location and the time coordinate of point p k (i), respectively.

Similar to the work of [START_REF] Pfoser | Novel approaches to the indexing of moving object trajectories[END_REF], we consider linear interpolation between successive sampled points p k (i), p k (i+1), so that each trajectory consists of a sequence of 3-D line segments e k (i) = p k (i)p k (i + 1), where each line segment represents the continuous moving of the object during sampled points.

Symbols Definitions

D Given MOD, D = {T 1 , T 2 , • • • , T N } T k k -th trajectory of MOD L k Number of points of T k p k (i) i-point of T k , p k (i) = (x k (i), y k (i), t k (i)) x k (i) x-spatial coordinate of p k (i) y k (i) y-spatial coordinate of p k (i) t k (i) time coordinate of p k (i) e k (i) 3-D line segment, e k (i) = p k (i)p k (i + 1) l k (i) Lifespan of e k (i), l k (i) = t k (i + 1) -t k (i) V k Voting trajectory descriptor of T k LP k Number of T k sub-trajectories P k (i) i -th sub-trajectory of T k V P k (i) Voting descriptor of i -th sub-trajectory of T k in D Nl k (i) Normalized lifespan vector of i -th sub-traj. of T k in D S Sub-trajectory sampling set of D V P k (i) Voting descriptor of i -th sub-traj. of T k in S SR(S)
Representativeness function of S The goal of this work include:

• the automatic segmentation of the given trajectories T k , into "homogenous" sub-trajectories according to their "representativeness" in MOD and • sampling of the most representative sub-trajectories of the MOD. Fig. 2 illustrates a scheme of the proposed system architecture. The following three sections formalize the issues of trajectory representativeness, trajectory segmentation, and sub-trajectory sampling, respectively. Table 1 summarizes the symbols' definitions used in this work.

Trajectory line segment representativeness

In this research, the "representativeness" in MOD is defined by extending the definition of density biased sampling (DBS) in point sets [START_REF] Kollios | Efficient biased sampling for approximate clustering and outlier detection in large datasets[END_REF] for trajectory segments. According to DBS, the local density for each point of the set is approximated by the number of points in a region, divided by the volume of the region. In our case, the "representativeness" of a trajectory segment is defined by the number of the objects that follow this segment along with time, space and direction.

Technically, "representativeness" is calculated by a voting process that is applied for each segment e k (i) of the given trajectory T k , improving a preliminary version of the proposed method, presented in [START_REF] Panagiotakis | Trajectory voting and classification based on spatiotemporal similarity in moving object databases[END_REF]. Thus, e k (i) will be voted by the trajectories of MOD, according to the distance of e k (i) to each trajectory. The sum of these votes is related to the number of trajectories that are close and similar to e k (i), having the number of trajectories in MOD as upper limit. We avoid to give each segment the ability of voting, because, in such a case, long trajectories moving around the same area could vote many times. Moreover, under this definition, voting has the physical meaning of how many objects co-move (i.e. co-location and co-existence) for a period of time. Thus, the voting results will be used to detect the representative paths and sub-trajectories.

The previous discussion implies that we have to compute the distance (d(e k (i), e m (j))) between two 3-D line segments. More specifically, we need to be able to identify the trajectories that are the closest (i.e. Nearest Neighbors -NNs) to a 3-D line segment, and this process should be done in an efficient and preferably in an incremental way, in the sense that the k-th NN can be obtained with very little additional work once the (k-1)-th NN has been found. These requirements imply that data should be indexed and, as such, the distance function should be obey the metric properties. To meet the above prerequisites, in this framework we use and extend the approach for incremental continuous NN search algorithms proposed by Frentzos et al. in [START_REF] Frentzos | Algorithms for nearest neighbor search on moving object trajectories[END_REF]. In this approach, which assumes that data are indexed by one of the two following members of the R-tree family for trajectory data (namely, the TB-tree [START_REF] Pfoser | Novel approaches to the indexing of moving object trajectories[END_REF] or the 3D-R-tree [START_REF] Theodoridis | Spatio-temporal indexing for large multimedia applications[END_REF]), the distance between 3-D line segments is defined as their Minimum Euclidean Horizontal distance (M EH, for short) during their common lifespan (see Fig. 3(a)) 1 . As proved in [START_REF] Frentzos | Algorithms for nearest neighbor search on moving object trajectories[END_REF], the Euclidean horizontal distance function between two 3-D line segments follows the quadratic form D(t) = √ at 2 + bt + c, where a, b, c are the factors of this trinomial (real numbers, a ≥ 0), and whose minimum is a closed formula that can be computed in O [START_REF] Guting | Moving Object Databases[END_REF]. Note that, by simply defining the distance between two 3-D line segments using M EH distance, will make the distance function dependable to a single value and not to the whole lifespan of the segments and, as such, very sensitive to small variations. As a counter-example consider the case where two objects move on the same road segment following different directions. At some timepoint they will meet, so their M EH will be (very close to) zero, but obviously this is not representative of their co-motion. As in our approach we aim at transforming trajectories into 1-D continuous descriptors that will smoothly describe the representativeness of the trajectories in the database, a more consolidated approach for measuring the distance between 3-D segments is needed.

Intuitively, an ultimate distance function would make use of the fact that for a single 3-D segment, different portions of it may have different NNs of other segments, and the sizes of these portions are taken into account. To this end, from the palette of algorithms proposed in [START_REF] Frentzos | Algorithms for nearest neighbor search on moving object trajectories[END_REF] we base on the best-first historical continuous trajectory NN IncT rajectoryN N Search algorithm (HCN N Q T , for short). The HCN N Q T (D, Q T , Q per ) query over the query moving object with trajectory Q T during a time period Q per returns a list of triplets for each of the k-NN, i.e.

LoT kN N =< • • • , (R j , T j , [t j-start , t j-end )), • • • > consisting of
the time-varying real value R j along with a moving object T j (belonging in database D) and the corresponding time period [t j-start , t j-end ) for which the nearest distance between Q T and T j stands. These timevarying real values R j are, in any time instance of their lifetime, smaller or equal to the distance between any 1. If there does not exist common lifespan, the distance is given by the diameter of the dataset. moving object in D and the query trajectory Q T . The time periods [t j-start , t j-end ) are mutually disjoint and their union forms Q per . With e j we denote the portion of the segment that corresponds to each triplet in LoT kN N .

Recall that, in our case, we want to apply such a NN query for each segment e k (i) of a trajectory, so following the notation of the HCN N Q T query, the Q T is the e k (i) and Q per becomes the respective temporal period of l k (i). To exemplify the proposed k-NN, in Fig. 3(b) we depict the 2-NN of the query segment e 1 . The 1-NN list includes e 2 for the interval [t 1 , t 3 ) and e 3 for the interval [t 3 , t 6 ]; the 2-NN list includes e 3 for the interval [t 1 , t 3 ) and e 1 for the interval [t 3 , t 6 ].

Given a segment e k (i) of a trajectory we calculate its distance from the NN segments e j provided by a LoT kN N , by computing the definite integral of the timevarying distance R j between the two segments during the same period, following the approach proposed by Frentzos et al. in [START_REF] Frentzos | Index-based most similar trajectory search[END_REF]:

d(e k (i), e j ) = t j-end tj-start R j (t)dt (1)
As R j follows the trinomial previously discussed, and as proved in [START_REF] Frentzos | Index-based most similar trajectory search[END_REF] its integral can be computed in O [START_REF] Guting | Moving Object Databases[END_REF] by distinguishing between two cases for the value of the non-negative factor a (a = 0 and a > 0). As also proved in [START_REF] Frentzos | Index-based most similar trajectory search[END_REF], d(e k (i), e j ) can be efficiently computed by approximating the integral with the Trapezoid Rule (

(Rj (tj-start)+Rj (t j-end ))•(t j-end -tj-start) 2
), providing also bounds for the approximation error.

Given the previous setting, let V k of L k -1 components be considered the voting trajectory descriptor ("representativeness" value) along the T k line segments, computed by using the previous distance function. Each component of this signal V k (i) corresponds to the number of votes of e k (i), i ∈ {1, 2, ..., L k -1} of T k . Section 4.1 describes the proposed global voting method that estimates V k (i). (c) Segmentation of T 1 using the proposed method (three sub-trajectories). (d) The two most representative subtrajectories according to the proposed method.

Trajectory Segmentation

In this section, the trajectory segmentation is presented. The goal of trajectory segmentation is a trajectory partitioning into sub-trajectories of homogenous representativeness irrespectively of their shape complexity.

Let V P k of LP k partitions be considered the voting sub-trajectory representativeness of T k trajectory. Let P k (i), i ∈ {1, ..., LP k } be the i-th sub-trajectory of T k , where LP k denotes the number of partitions of T k . Then according to previous definition, V P k (i) is the vector (descriptor) of votes along the line segments of P k (i) partition, showing how many objects follow each line segment of sub-trajectory. This value will turn out to be useful for sampling purposes, as it will be presented in the section that follows.

Fig. 4 illustrates an example of a MOD comprised by four trajectories ({T 1 , ..., T 4 }). In this figure, the time dimension has been ignored for visualization reasons. The segmentation of T 1 according to the TRACLUS method [START_REF] Han | Trajectory clustering: a partition-and-group framework[END_REF] and the proposed method are illustrated in Figs. 4(b) and 4(c), respectively. Indeed, most of the proposed techniques in literature like TRACLUS first simplify the trajectories into large line segments and then apply a grouping of similar of line segments, see six subtrajectories on Fig. 4(b). On the contrary, according to our approach, the segmentation is applied on the original trajectory by taking into account its neighborhood in the rest of the MOD. In addition, there is not any constraint on sub-trajectory shape complexity, yielding a segmentation that is related only on line segment representativeness (voting results), see three sub-trajectories on Fig. 4(c). In Fig. 4(c), the line thickness is related with the line segments representativeness (voting results). It holds that the sub-trajectory representativeness values will be almost the same along the line segments of each sub-trajectory (see Fig. 4(c)). According to the proposed methodology, the top-2 representative sub-trajectories ({S 1 , S 2 }) that best describe the MOD are illustrated in Fig. 4(d).

Subtrajectory Sampling

The trajectory segmentation provides homogenous subtrajectories concerning their representativeness. The goal of sub-trajectory sampling is to select the most representative subset of these sub-trajectories. By selecting the higher voting segments as in [START_REF] Panagiotakis | Trajectory voting and classification based on spatiotemporal similarity in moving object databases[END_REF], which sounds to be an obvious decision, the high density regions of the MOD will be oversampled, resulting in a non-representative sample when the aim is to cover the whole space-time of the MOD as much as possible.

On the contrary, in this paper, we propose the sampling to be done by minimizing a formula (see Equation 3), taking into account the votes (i.e. representativeness) V P k (i) in the original MOD as well as the votes V P k (i) in the sampling MOD and the vector of the lifespans of the line segments divided by the lifespan of the trajectory (N l k (i)). V P k (i) is computed in sampling MOD using the voting method similar to V P k (i) computation in the original MOD. Thus, it takes into account the fact that the sampling set should contain dissimilar representative sub-trajectories. So, our goal is that the sampling set should contain high voting trajectories of the MOD and, at the same time, should cover the whole 3-D space as much as possible. This means that the resulting set should avoid to contain similar sub-trajectories.

Let S denote the sampling set, so that S k (i) is one, if P k (i) sub-trajectory belongs to the sampling set, and zero otherwise. According to the previous properties, the number of moving objects of the original MOD that are represented in sampling set SR(S), should be maximized over the time. This is formalized in Equations 2 and 3.

SR gain (k, i) = |P k (i)| j=1 V P k (i)(j) • N l k (i)(j) • (1 -V P k (i)(j)) (2) SR(S) = N k=1 LP k i=1 S k (i) • SR gain (k, i) (3)
SR gain (k, i) expresses the gain of SR(S) if we add in sampling set the i th sub-trajectory of the k th trajectory of the MOD. |P k (i)| denotes the number of line segments of the ith sub-trajectory of T k . The values V P k (i)(j), N l k (i)(j) and V P k (i)(j) denote for the votes in D, the normalized lifespan and the votes in S of the j th line segment of i th sub-trajectory of k th trajectory of the MOD, respectively. In other words, the number of objects of the MOD that follows the sub-trajectories of S is maximized over the time. Fig. 4(d) illustrates the top-2 representative sub-trajectories S 1 and S 2 according to the proposed method (S 1 is the top-1st since it represents all four objects, while S 2 is the top-2nd since it represents two objects. Let us suppose that the normalized lifespan of S 1 is 1 2 and the normalized lifespan of S 2 is 1 3 . It holds that S 1 will be selected first, since the SR gain of S 1 is about 3 2 (see Equation 2). S 2 will be selected next, since the SR gain of S 2 is about 2 3 (see Equation 2). In both of the cases, the corresponding V P k (i)(j) is zero, since S 1 and S 2 are defined into different time periods. Next, the SR gain of other sub-trajectories will be zero or a negative value. Therefore, the proposed sampling method will select these two sub-trajectories, covering the time and space of the MOD.

METHODOLOGY

In this section, the proposed methodology is presented, consisting of the trajectory voting, the trajectory segmentation and the sub-trajectory sampling methods.

Global Voting Method

This section describes the Global Voting Algorithm (GVA). The input of the algorithm is a MOD D = {T 1 , T 2 , • • • , T N } indexed by a R-tree-like structure such as the TB-tree [START_REF] Pfoser | Novel approaches to the indexing of moving object trajectories[END_REF] or the 3D-R-tree [START_REF] Theodoridis | Spatio-temporal indexing for large multimedia applications[END_REF], as described in [START_REF] Frentzos | Algorithms for nearest neighbor search on moving object trajectories[END_REF], a trajectory T k ∈ D and an intrinsic parameter σ > 0 of the method. The output of the method is the vector V k of L k -1 components that can be considered as a trajectory descriptor along the line segments e k (i), i ∈ {1, 2, ..., L k -1} of trajectory T k (recall from Section 3, that L k denotes the number of points of T k trajectory). As such, each component of the vector V k (i) corresponds to the number of votes (representativeness) for each e k (i) of T k .

According to the problem formulation presented in Section 3, for each line segment e k (i) of T k , the proposed GVA algorithm incrementally identifies the NN segments of other trajectories T j ∈ D, j = k. For each set of NN segments represented by the list of segments/triplets in LoT kN N , the distances d(e k (i), e j ) from the corresponding segments are computed. These distances are used to define the voting function V (e k (i), LoT kN N ), which quantifies the representativeness of the line segment to a LoT kN N . In the literature, a lot of voting functions have been proposed, like step functions or continuous functions [START_REF] Patterson | Artificial Neural Networks[END_REF]. In this work, we have selected to use the continuous function of a gaussian kernel, which is widely used in a variety of applications of pattern recognition [START_REF] Yuan | Adaptive spherical gaussian kernel in sparse bayesian learning framework for nonlinear regression[END_REF]. Formally,

V (e k (i), LoT kN N ) = ∀ej ∈LoT kN N ,T k =Tj e -d 2 (e k (i),e j ) 2•σ 2 (4)
Note that for data collected from GPS devices where the segments are very small (due to the high sampling rate), in practice the previous function degenerates to the computation of a single gaussian kernel, as the HCNN of a segment results in a LoT kN N containing only one NN. This is actually verified in our experiments where we used real GPS datasets. However, in cases where the datasets are highly compressed (in applications where storage cost is important), e.g. by an approach like the one proposed in [START_REF] Meratnia | Spatiotemporal compression techniques for moving point objects[END_REF], this may have an influence in the smoothness of the V k (i) descriptor. We leave such a study as future work. The control parameter σ > 0 shows how fast the function ("voting influence") decreases with distance. Given the previous assumption, and according to Equation 4, it holds that 0 ≤ V (e k (i), LoT kN N ) ≤ 1. If d(e k (i), e j ) is close to zero, the voting function gets its maximum value, i.e. 1. This means, that there exists a line segment of T j that is being (in time, space and direction) very close to e k (i). Otherwise, if d(e k (i), e j ) is high, e.g. greater than 5 • σ, the voting function results in almost 0, meaning that T j is very far away from e k (i).

The use of a continuous voting function, like the gaussian kernel, gives smooth results for small changes on parameters (σ in our case), and the possibility to get decimal values as results of voting process increasing the robustness of the method. However, σ depends on space units the object movements of MOD and it is difficult to tune it. We have solved this problem by estimating σ as the percentage (e.g. 0.1%) of dataset diameter (maximum space distance). This percentage can be kept almost constant for every dataset. Finally, V k (i) is computed by getting the sum of votes for all of the nearest neighbor segments of trajectories T j ∈ D, j = k, according to GVA (see Algorithm 1). Note that for each e k (i) (line 1) the algorithm starts a loop (loop 3) wherein each time it incrementally retrieves the next LoT kN N (line 5, function Get N ext HCN N ), for which its representativeness is computed. This measure is added to the V k (i) if the current voting is significant (i.e. larger than a small value (ε is a very small number comparing to one, e.g. ε = 10 -3 ), lines 6,7), otherwise the loop is terminated (line 9), which means that we stop retrieving subsequent NN. The loop is also terminated if there are no other NN given the spatiotemporal criteria of e k (i) (line 11). Given the above discussion, a nice property that holds for the proposed local trajectory descriptor V k is that it smoothly changes over the trajectory segments.

input : An (indexed) database D = {T 1 , T 2 , • • • , T N }, a trajectory T k in D and the parameters σ, ε. output: Voting vector V k of T k for i = 1 to L k -1 do 1 V k (i) = 0 2 repeat 3 LoT kN N = 4 Get N ext HCN N Q T (D, e k (i), l k (i)) V (e k (i), LoT kN N ) = 5 ∀ej ∈LoT kN N ,T k =Tj e -d 2 (e k (i),e j ) 2•σ 2 if V (e k (i), LoT kN N ) > ε then 6 V k (i) = V k (i) + V (e k (i), LoT kN N )

Trajectory Segmentation Method

Having presented the voting procedure in the previous section, the next step is to provide a solution to the trajectory segmentation problem defined in section 3.2. For this purpose we propose the Trajectory Segmentation Algorithm (TSA) (see Algorithm 2). The input of the algorithm is the normalized trajectory voting signal V k , and two intrinsic parameters w, τ of the method. The normalization is done by dividing V k by the maximum over all V k , thus bounding V k ≤ 1. The output of the method is the segmentation P k of T k into LP k partitions, where LP k is automatically estimated by the proposed scheme.

The method uses two sequential sliding signal windows W 1 and W 2 of w samples estimating the sample when the "difference" between the two windows is maximized.

This methodology has been successfully applied on sound signal segmentation [START_REF] Panagiotakis | A speech/music discriminator based on rms and zero-crossings[END_REF] and P Phase Picking of seismic signals [START_REF] Panagiotakis | Automatic P-Phase Picking Based on Local-Maxima Distribution[END_REF]. To facilitate the discussion, Fig. 5 illustrates the two sequential sliding windows W 1 (light gray horizontal lines) and W 2 (heavy gray vertical lines) locating at sample n on the given voting signal V k . First, as the two windows slide, the two means m 1 , m 2 and two variances σ 2 1 , σ 2 2 of two sequential signal windows and W 1 , W 2 locating at sample n are estimated, respectively (see Fig. 5) (see lines 1-7 of Algorithm 2). The next equations define m 1 and σ 2 1 ; m 2 and σ 2 2 are similarly defined in window W 2 .

m 1 = 1 w i∈W1 V k (i) (5)
σ 2 1 = 1 w i∈W1 (V k (i) -m 1 ) 2 (6)
The symmetric Mahalanobis Distance [START_REF] Panagiotakis | Automatic P-Phase Picking Based on Local-Maxima Distribution[END_REF], [START_REF] Mahalanobis | On the generalized distance in statistics[END_REF], presented in Equation 7, is used to measure the distance between the two windows locating at sample n.

d(n) = (m 1 -m 2 ) 2 • ( 1 σ 2 1 + 1 σ 2 2 ) ( 7
)
The symmetric Mahalanobis Distance has been selected since it outperforms other frequently used distances like symmetric Kullback-Leibler (KL2) or Bhattacharyya [START_REF] Zhou | From sample similarity to ensemble similarity: Probabilistic distance measures in reproducing kernel hilbert space[END_REF] in segmentation problem, which has been shown in [START_REF] Panagiotakis | Automatic P-Phase Picking Based on Local-Maxima Distribution[END_REF].

A sample n is characterized as border point (start of a new partition), The first statement ensures that the difference of two windows W 1 , W 2 will be high enough while the second selects the sample n, where d(n) is locally maximized, meaning that the difference of two windows W 1 , W 2 is locally the highest. Both statements are related with the number of partitions LP k , while the second ensures that the minimum partition size is w samples (w 3-D line segments) [START_REF] Panagiotakis | Automatic P-Phase Picking Based on Local-Maxima Distribution[END_REF].

Parameter w sets the minimum size of a partition, so w depends on the given dataset and the user preferences. In other words, w expresses the minimum number of line segments that can define a sub-trajectory. In addition, w is analogous to sampling rate of the trajectories (e.g. if a MOD has double sampling rate, then w should be multiplied by two). TSA needs to estimate mean and variance measures, thus we need at least two samples (w ≥ 2). Low values on w can affect the robustness of mean and variance estimation yielding false alarms (over-segmentation). High values on w gives more robust results, and it will affect the results of the method, only if there is a sub-trajectory with length less than w that will not be detected. According to our experiments, when w ∈ [START_REF] Lee | Traclass: trajectory classification using hierarchical region-based and trajectory-based clustering[END_REF][START_REF] Andrienko | A visual analytics toolkit for cluster-based classification of mobility data[END_REF] most of sub-trajectories were robustly detected without important false alarms.

Regarding τ , it should be a positive number close to zero, in order to be sure that TSA will detect all the sub-trajectories. According to our experiments, when τ ∈ [0.001, 0.1] most of sub-trajectories were detected without important false alarms, since this parameter is related with the segmentation sensitivity of our method. As τ increases, the number of sub-trajectories reduces. It holds that τ can be set as a positive number close to zero (e.g. 0.01), due to the fact that in the first step, we perform normalization by dividing V k by the maximum over all V k , thus bounding V k ≤ 1.

Subtrajectory Sampling

In the previous sections, we have presented our methodology for segmenting the trajectories of a MOD into subtrajectories using the votes gathered for the MOD. In this section, we exploit on this knowledge in order to select the top representative sub-trajectories to be the result of a sampling process. In particular, we propose the Subtrajectory Sampling Algorithm (SSA). SSA solves the sampling problem defined in section 3.3 (see Algorithm 3). The input of the algorithm is the set of sub-trajectories of the MOD P k as estimated by TSA, the voting V P k (i) input : A normalized trajectory voting vector V k and two parameters w and τ . output: The segmentation P k of LP k partitions. and the normalized lifespan N l k (i) vectors of the trajectory segments. The output of the method is the subtrajectory sampling set S consisting of M samples. M can be given as input to the method or (more interestingly) it can be automatically estimated by the proposed scheme.

for n = w + 1 to L k -w -1 do 1 m 1 = 1 w n-1 i=n-w V k (i) 2 m 2 = 1 w n-w-1 i=n V k (i) 3 σ 2 1 = 1 w n-1 i=n-w (V k (i) -m 1 ) 2 4 σ 2 2 = 1 w n-w-1 i=n (V k (i) -m 2 ) 2 5 d(n) = (m 1 -m 2 ) 2 • ( 1 σ 2 1 + 1 σ 2 2 ) 6 end 7 LP k = 1, a = 1 8 P k (1) = [a : L k -1] 9 for n = w + 1 to L k -w -1 do 10 m = max i∈{n-w,...,n+w} d(i) 11 if d(n) > τ ∧ d(n) ≥ m then 12 P k (LP k ) = [a : n -1] 13 P k (LP k + 1) = [n : L k -1]
The goal of SSA is the maximization of the number of sub-trajectories SR(S) of the original MOD that are represented in the sampling set (see Equation 3). The complexity of an exhaustive algorithm that would search for all the possible solutions in order to maximize Equation 3 is O( N M ). On the other hand, our proposed algorithm suboptimally solves the problem in O(N • M ) iterations by applying iterative optimization.

SSA starts with an empty sampling set (S k (i) = 0), where S k (i) is defined in section 3.3. In each iteration step, SSA adds in sampling set an unselected subtrajectory of MOD that maximizes Equation 3. This is equivalent with the maximization of SR(S) gain SR gain (k, i) (see Equation 2). Recall that SR gain (k, i) expresses the gain of SR(S) if we add in sampling set the i th sub-trajectory of k th trajectory of the MOD. According to the proposed algorithm, it holds that SR(S) gain is a monotonically decreasing function as sampling size increases. Since V P k (i)(j) ≥ 0, and recalling Equation 2, it holds that:

SR gain (k, i) ≤ |P k (i)| j=1 V P k (i)(j) • N l k (i)(j) (8)
Let SR gain zero (k, i) denotes the sum of Equation 8, that corresponds to the SR gain (k, i) when S = ∅.

SR gain zero (k, i) = |P k (i)| j=1 V P k (i)(j) • N l k (i)(j) (9)
Therefore, an efficient way to reduce the computation cost of the maximization of SR gain , is to keep the sub-trajectories' indexes in a vector (vec) sorted in descending order by SR gain zero (line 6 of Algorithm 3). Instead of computing SR gain for each sub-trajectory of the MOD in order to find the maximum, we get the j th , j ∈ {1, 2, ..., N k=1 LP k } (line 9 of Algorithm 3) sub-trajectory from the vector and we compare its voting with current highest SR gain (SR max gain ) (lines 11-15 of Algorithm 3). First, SR max gain is set to minus infinity (line 8 of Algorithm 3). Let the j th sub-trajectory of vec be the i th sub-trajectory of k th trajectory of the MOD (line 10 of Algorithm 3). This loop terminates, if SR max gain will be higher than the V P k (i), because it holds that the SR gain of each sub-trajectory from the rest sub-trajectories from the list would be lower than its SR gain zero (see Equation 8) and lower than SR max gain , since we sort them in descending order by SR gain zero . Then, the v th sub-trajectory of u th trajectory of the MOD that corresponds to SR max gain is added on sampling set (line 18).

SSA terminates when the size of the sampling set reaches M (which is an input parameter). Alternatively, the algorithm terminates if SR gain (k, i) is lower than a given threshold ( , a positive number close to zero, see line 19 of Algorithm 3). If this threshold is set to zero, it means that SR(S) has been maximized, as the latter is increased in each step by SR max gain . An advantage of the proposed method is that it provides a deterministic solution in contrary with other probabilistic techniques [START_REF] Kollios | Efficient biased sampling for approximate clustering and outlier detection in large datasets[END_REF], [START_REF] Nanopoulos | Indexedbased density biased sampling for clustering applications[END_REF] that provide a randomly constructed sampling set trying to fit it to a desired distribution or usersupervised, explorative approaches [START_REF] Andrienko | A visual analytics toolkit for cluster-based classification of mobility data[END_REF], [START_REF] Andrienko | Interactive visual clustering of large collections of trajectories[END_REF] that base on clustering, which in turn depends on the algorithm itself, its parameters and the distance function adopted.

Computational Complexity Issues

Concerning the complexity of GVA, and given the use of the R-tree-like structures, the computational cost for each line segment e k (i), of T k is O(log( L•N )), where L denotes the mean number of trajectory points. Executing GVA for each trajectory of the database, the total computation cost is

O( L • N • log( L • N )).
Concerning the complexity of TSA, the computational cost for the segmentation of a trajectory T k is O(L k ), since the mean and the variance of a sliding window can be estimated recursively in O( 1) by the mean and the variance of the sliding window of the previous step. For example, let m 1 be the mean of the window W 1 = V k (n-W : n-1) and ḿ1 be the mean of the next window

Ẃ1 = V k (n -W + 1 : n). Then, it holds that ḿ1 = m 1 • w -V k (n -W ) + V k (n) w ( 10
)
input : The voting vector V P k (i), the normalized lifespan vector N l k (i) and the segmentation P k (i), ∀k ∈ {1, ..., N }, i ∈ {1, ..., LP k }, . The size of sampling set M . output: The sampling set S k (i). If we perform TSA for each trajectory of the MOD, then the total computation cost is O( L • N ).

Concerning the complexity of SSA, the sorting of

V P costs O(N • L LP • log(N • L LP ))
, where LP denotes the mean number of trajectory segments and L LP denotes the mean number of sub-trajectory points. Next, in the worst case, the method computes SR gain (k, i)

(M • N k=1 LP k = M •LP •N ) times. SR gain (k, i) computa- tion requires the computation of V k (S). The cost of V k (S) is O( L LP •log( L•M ))
, since the maximum size of sampling set is M and the mean number of a sub-trajectory points is L LP . Therefore, in the worst case, the overall cost of the proposed segmentation-and-sampling method is

O(N • L LP •log(N • L LP )+ L•M •N •log( L•M )).
In the best case, the break of line 16 of the algorithm can stop the intrinsic loop in two (O(1)) instead of N • L steps. In this case, the cost turns out to be Ω(N

• L LP •log(N • L LP )+M •log( L•M )).
Conclusively, the most computationally intensive part of the proposed method is the GVA with O( L • N • log( L • N )) complexity. In turn, the most time consuming step in GVA is the search of the nearest neighbors of a trajectory in a given time period. In order to make the application of our approach feasible to large datasets, we have adopted efficient Continuous Nearest Neighbor query processing techniques [START_REF] Frentzos | Algorithms for nearest neighbor search on moving object trajectories[END_REF], where trajectories are indexed by R-tree-like structures. Scalability experiments under various cases for such queries have been presented in [START_REF] Frentzos | Algorithms for nearest neighbor search on moving object trajectories[END_REF], where it has been shown their applicability in large datasets, with an almost linear behavior with the size of dataset. Actually, this conclusion is in accordance with the above theoretical analysis of the computational complexity of the proposed method.

On the effect of the MOD extension

As already mentioned, the proposed method is deterministic, which implies that different invocations for a given MOD will have the same result. This is a crucial and distinct characteristic of our approach w.r.t. other sampling approaches. In this section, we discuss the behavior of the proposed methods when the MOD is extended either in space or time dimensions. More specifically, we study the following scenario: In a given dataset S, we add trajectories which come from a different spatial or temporal space (extension of the MOD). The question is whether the GVA, TSA and SSA are affected by such an extension?

According to this scenario, the results of the voting procedure (GVA) will not change concerning the given dataset S. The new trajectories do not affect the trajectory descriptors of the trajectories of S, since they exist in different spatial or temporal space (see Equation 4). Therefore, the results of TSA will be exactly the same concerning the given dataset S. Similarly, the new sampling set will contain the same sub-trajectories as the sampling set of S (when the algorithm terminates if SR gain is lower than a given threshold), with some additional samples selected from the new trajectories, since the input of the SSA algorithm concerning S remains the same. In other words, the sampling set of the union of two distinct (in space and/or in time dimension) MOD, is the same with the union of the sampling sets, when the sampling process is performed to each MOD independently. Therefore, the effectiveness of the proposed method is not affected by the extension of the MOD.

EXPERIMENTAL RESULTS

In this section, the experimental results of our performance study are presented. In order to evaluate the effectiveness and robustness of the proposed scheme, we have performed two different experiments, with synthetic and real datasets. The implementation of the proposed algorithms was done in Matlab. A linear scan instead of an R-tree-based index scan was performed for HCNN search. However, this does not affect the effectiveness of our approach but only scalability. In other words, the results of GVA produce the same results either by using the R-tree-based index-based implementation, or by adopting a linear scan for the HCNN search.

In all of presented experiments, we have used τ = 0.01 and w = 8. Similar results were obtained when we set w equal to a percentage (e.g. 20% -25%) of the mean number of line segments per trajectory. In addition, we have run several experiments that show that the quality of sampling results not sensitive to those parameters. These parameters affect only the results of segmentation (e.g. when τ and w decrease, TSA will gives an oversegmentation of the MOD).

Experimenting with Synthetic Data

In order to measure the performance of segmentation and sampling for movements that can be created by vehicles, we have created a synthetic MOD (SMOD). Note that in order to evaluate our approach, well known synthetic trajectory generators [START_REF] Brinkhoff | A framework for generating network-based moving objects[END_REF] cannot be utilized as we would not be aware of the ground truth of the generated trajectories. As such, assume the following points A(0, 0), B(1, 0), C(2, -1), D(2, 1) and E(4, 0) be the destination nodes of a simple graph upon which objects move. We assume that the half of the objects are moving with normal speed (5 units per second) and the rest of them are moving with high speed (10 units per second). The objects are moving under the following scenario (rules), for a lifetime of one second.

• There are four one-directional roads (B → D, D → E, E → C, C → B) and two bidirectional roads (A B, B E). • At t = 0 sec, all objects of MOD start from (in a real application, very close to) point A. Thus, the first destination of all objects is point B. • When an object arrives at a destination point, it ends its trajectory with a probability of 10%. Otherwise, it continues with the same speed to the next point according to the road network rules. If there exist more than one possible next points, it decides randomly the next destination. As an example, if an object moves from A to B, then the next destination could be A, D or E with the same probability 30% (the rest 10% is the case where the object ends its trajectory at B). Fig. 6 illustrates the 2-D map of SMOD with the four one-directional and the two bidirectional roads.

According to the rules that the generated trajectories should obey, there exist 78 distinct sub-trajectories in the MOD (ground truth). Table 2 illustrates the spatial (first column) and temporal (second column) projection

Path Time Support A → B [0, 0.1] 50% A → B [0, 0.2] 50% B → A [0.2, 0.4] 15% B → A [0.1, 0.2] 15% B → D [0.1, 0.26] 15% B → D [0.2, 0.52] 15% B → E [0.1, 0.4] 15% B → E [0.2, 0.8] 15%
TABLE 2: Some representative sub-trajectories and their support in the synthetic MOD (ground truth).

of the top-8 representative sub-trajectories, along with their support (third column), i.e. the average percentage of objects in the MOD that "vote" for the corresponding sub-trajectory. These percentage values have been computed under the aforementioned assumptions of the previous paragraph. The percentage of objects of MOD that follows a path could change between experiments due to random decisions of the prementioned scenario.

According to this synthetic MOD, the possible endpoints of each trajectory partition are points A, B, C, D, and E (ground truth). The above dataset is ideal for the purposes of experimentation: It consists of a high number of predefined sub-trajectories with each one being associated with a known representativeness (as ground truth).

In order to measure the stability of our method to noise effects, we have added Gaussian white noise of different Signal to Noise Ratio (SNR) levels, measured in db, to spatial coordinates of synthetic MOD. The synthetic MOD with additive noise of SNR = 50 db and SNR = 30 db projected in 2-D spatial and 3-D spatiotemporal space is illustrated in Fig. 7.

According to this synthetic MOD, there exist eight different paths that are defined by four one-directional roads and two bidirectional roads. We have used 600 moving objects with 1934 sub-trajectories. The mean value and the standard deviation for number of line segments per trajectory (L k ) are 82.03 and 30.95, respectively. The proposed GVA/TSA gives between 1800 and 2200 sub-trajectories under any case of additive noise. The average accuracy of trajectory segmentation method is measured using the average distance error (in spatial coordinates) between estimated segment borders and the ground truth borders. In order to get a normalized error, we divide the average distance error by the half of each corresponding sub-trajectory length getting the normalized segmentation error in the range [0, 1]. Fig. 8 plots this segmentation error for SMOD under additive noise of various SNR. Concerning the second trajectory, in Fig. 9(d), the samples 1, 10, 30, 45, 80, 90 and 100 (horizontal axis) correspond to ground truth ends of borders. In this case, we get errors on borders' approximation less than 3 samples due to additive white noise.

Next, we present the results of the proposed sampling method, in other words the results of the proposed GVA/TSA/SSA three-step methodology. Fig. 10 illustrates results of sampling method for synthetic MOD with additive noise of SN R = 50db. It holds that SR(S) gain (see Fig. 10(b)) is almost zero when the sampling set size is greater than 60 meaning that the sampling set of the top-60 representative sub-trajectories is an appropriate minimum number of sampling set size. We get 60 sub-trajectories instead of 78 (ground truth), since some sub-trajectories are not perfectly segmented, e.g. see the sub-trajectory A → B → A of Fig. 9(a). Figs. 10(c), 10(d) illustrate sampling/visualization of the dataset using the top-60 representative sub-trajectories under the proposed method using a colormap according to SR(S) gain and line segments representativeness, respectively. The most representative sub-trajectories belong to two bidirectional roads (A B (629 movements), B E (523 movements)) which always appear to have high traffic, while the less representative sub-trajectories be-long to a one-directional road (C → B) which has the lowest traffic (66 movements) mainly appeared during two time periods (see Fig. 7(b)). The sub-trajectories of sampling set are projected in spatiotemporal 3-D space using a colormap according to SR(S) gain (Fig. 11(c)) or line segments representativeness (10(d)). We have used dark blue-to-dark red colormap according to the corresponding to segments representativeness (dark red color for high values, dark blue color for low values). Moreover, the proposed method is not affected by the trajectories' shape, yielding high performance results for both the non-straight and straight movements. Fig. 11 illustrates results of sampling method for synthetic MOD with additive noise of SN R = 30db. The selected subtrajectories represent well the dataset, getting similar results with the case of SN R = 50db.

Experimenting with Real Data

Moreover, we have evaluated the proposed scheme using two real MODs: temporal domain, called "Milano-I" and "Milano-II", respectively. The "Milano-I" MOD contains 5910 trajectories all existing inside a temporal period of 100 minutes duration. For this MOD the GVA/TSA produces 9835 sub-trajectories. The mean value and the standard deviation for number of line segments per trajectory (L k ) are 31.04 and 29.41, respectively. The "Milano-II" consists of 11500 randomly selected trajectories from the original dataset, spanning inside a temporal period of 5 days duration. The mean value and the standard deviation of the number of line segments per trajectory (L k ) are 30.77 and 31.34, respectively. For this MOD the GVA/TSA produces 12895 sub-trajectories. As such, the "Milano-II" is a very sparse in time dataset and we have used it in order to show the limitations of our method. Fig. 15 illustrates the trajectories of "Athens trucks" MOD projected in 2-D spatial space ignoring time dimension (Fig. 15(a)) and in 3-D spatiotemporal space (Fig. 15(b)). Fig. 16 illustrates the trajectories of "Milano-II" MOD projected in 2-D spatial space ignoring time dimension (Fig. 16(a)) and in 3-D spatiotemporal space (Fig. 16(b)) (the corresponding figures for "Milano-I" MOD are quite similar). It is clear that the full visual-ization cannot be effective, due to the large number of projected trajectories in almost the same time and space. In order to solve this problem, we have used the results of SSA sampling the most representative sub-trajectories of the MOD. Fig. 12 illustrates the functions SR(S) and SR(S) gain, as sampling set size increases for the "Athens trucks" MOD. It holds that SR(S) gain is low when the sampling set size is greater than 150, meaning that the most appropriate number of sampling set size is at least 150. SR(S) gain always decreases as sampling size increases. Similarly, Fig. 13 illustrates the functions SR(S) and SR(S) gain, as sampling set size increases for the "Milano-I" MOD. It holds that SR(S) gain is very low, when the sampling set size is greater than 120, meaning that the most appropriate number of sampling set size is about 120. Fig. 14 illustrates the functions SR(S) and SR(S) gain, as sampling set size increases for the "Milano-II" MOD. This is a very sparse in time dataset, since the maximum value of SR(S) gain is about 10, while the dataset contains more than 11.000 trajectories. It holds that SR(S) gain is less than one, when the sampling set size is greater than 150.

Fig. 15(c) illustrates a sampling/visualization of the dataset using the top-150 representative sub-trajectories under the proposed SSA method for "Athens trucks" MOD. We have used dark blue-to-dark red colormap according to the corresponding to SR(S) gain (dark red color for high values, dark blue color for low values). Similarly, Fig. 15(d) uses a dark blue-to-dark red colormap according to the corresponding to subtrajectories representativeness V P k (i) (dark red color for high values, dark blue color for low values) for "Athens trucks" MOD. The low values of sub-trajectories representativeness is due to the fact that "Athens trucks" MOD is sparse on its temporal dimension. Fig. 16(c) and 16(d) illustrate a sampling/visualization of the dataset using the top-150 representative sub-trajectories under the proposed SSA method for "Milano-II" MOD using colormaps according to SR(S) gain and V P k (i), respectively. The sub-trajectory sampling set includes the perimeter of trajectories set as well as some intrinsic trajectories.

Comparison with other sampling techniques

The results of the proposed sub-trajectory sampling method have been compared with standard sampling techniques (random and stratified sampling) as these have been proposed in the literature [START_REF] Thompson | Sampling[END_REF]. The sampling techniques have been evaluated on the same dataset of sub-trajectories that is provided by the proposed TSA. According to random sampling, each sub-trajectory is chosen randomly and entirely by chance, such that each sub-trajectory has the same probability of being chosen at any stage during the sampling process. According to stratified sampling, the sub-trajectories need first to be divided into strata (groups). The only relevant approach that has been proposed in the literature for MOD is [START_REF] Andrienko | A visual analytics toolkit for cluster-based classification of mobility data[END_REF], which uses the state-of-the-art T-OPTICS algorithm [START_REF] Nanni | Time-focused clustering of trajectories of moving objects[END_REF]. As such, in a way, we compare our method with the state-of-the-art stratified sampling technique proposed so far in the literature [START_REF] Andrienko | A visual analytics toolkit for cluster-based classification of mobility data[END_REF]. Following the stratification, the number of sub-trajectories that are selected from each of the strata is analogous to the size of the stratum (i.e. the number of sub-trajectories in this group). Finally, the sub-trajectories are randomly selected from each stratum. In order to measure the performance of sampling, we have used the root mean square error (RMSE) metric, that is defined by the square root of the mean of minimum squared distances between the line segments of dataset (D) and the line segments of the

e u (v) is the line segment from sampling set that minimizes d(e k (i), e u (v)). Intuitively, the smaller the RMSE metric the better the sampling, as this implies better coverage of the space-time of the MOD, or equivalently, it means that less "useless" sub-trajectories have been sampled, where by "useless" we mean trajectories that are already represented in the sampling set by some other sub-trajectory. According to the experimental results, the proposed method outperforms the other sampling techniques in all the cases of synthetics and "Milano-I" MODs. Concerning, the "Athens trucks" MOD, the proposed method outperforms the other sampling techniques in more than 85% of the cases. This is due to the fact that this MOD is temporally sparse. When a MOD is temporally sparse, most of the sub-trajectories have low representativeness, since the number of common time periods between different trajectories is low (see Equation 1). Therefore, in such cases, a sampling method that takes into account the time, as the proposed method does, could not always give the best results, since the selected samples do not represent other samples due to low representativeness values of sub-trajectories. In other words, the sparser in time a MOD is, the closer the behavior of our sampling method is with the random sampling. This is rational, as in highly sparse MOD (where an example of a highly sparse MOD could be the one in which each trajectory exists in different time periods) it is like omitting one of the two dimensions under which the problem definition is formulated, which obviously is a different kind of problem. As shown in Fig. 17(e), the above intuition is verified in this experiment for the "Milano-II" MOD, which is a very sparse in time. Nevertheless, the effectiveness of the proposed method is still high, since it outperforms or gives almost the same results with the other sampling techniques in more than 75% of the cases. In general, by using the standard evaluation metric of RMSE, which is independent to any distribution that a dataset may follow, the results clearly show that our methodology covers better that space-time of a MOD.

CONCLUSIONS

In this paper, we have discussed the problem of finding representative sub-trajectories in a MOD. Especially, we addressed this issue by segmentation and sub-trajectory sampling based on global spatiotemporal similarity of trajectories. In particular, we have proposed three algorithms: GVA, TSA and SSA for trajectory voting, segmentation and sub-trajectory sampling, respectively. GVA extends the density biased sampling (DBS) from point sets [START_REF] Kollios | Efficient biased sampling for approximate clustering and outlier detection in large datasets[END_REF] to trajectory segments providing a local trajectory descriptor per line segment that is related to line segment representativeness. Next, TSA automatically and effectively estimates the number of sub-trajectories and their borders, separating each trajectory of MOD into homogenous partitions concerning their representativeness. Finally, SSA is applied over the resulting partitions providing the most representative sub-trajectories of the MOD, also taking into account that high density regions of the MOD should not be oversampled. SSA can be automatically terminated by thresholding the number of moving objects of the original MOD that are represented in sampling set SR(S). Moreover, the index-based voting algorithm, which is the computationally most expensive step in our framework, and the polynomial computational cost of the proposed algorithms makes the scheme applicable to large databases. In our approach, contrary to related work, the temporal dimension of the MOD is taken into consideration, while there is not any inherent constraint on sub-trajectory complexity and shape, yielding trajectory segmentation and sub-trajectory sampling that are related only to representativeness. We have evaluated the proposed method under real and synthetic databases, and the experimental results show the effectiveness and robustness of the proposed scheme.

As future work, we plan to investigate the applicability of the proposed method for (sub-)trajectory clustering. The idea is that MOD clustering can be provided concurrently with MOD sampling. It holds that each sub-trajectory of the sampling set has been voted by different sub-trajectories of the MOD (cluster), under the minimization of the objective function proposed in the current work. Therefore, each sub-trajectory of the sampling set can be considered as a cluster representative (i.e. a seed around which a cluster is formatted). This is a different tactic as the one followed in [START_REF] Han | Trajectory clustering: a partition-and-group framework[END_REF]. In the same context, outliers [START_REF] Lee | Trajectory outlier detection: A partition-and-detect framework[END_REF] can be discriminated from low values in voting sub-trajectory descriptor (V k ).

Fig. 1 :

 1 Fig. 1: (a) A MOD of four trajectories. (b) The two most representative sub-trajectories.

Fig. 2 :

 2 Fig. 2: Scheme of the proposed system architecture.

Fig. 3 :

 3 Fig. 3: (a) The distance between 3-D line segments is defined as their Minimum Euclidean Horizontal distance during their common lifespan. (b) The 2-NN of the query segment e 1 .

Fig. 4 :

 4 Fig. 4: (a) A MOD of four trajectories. (b) Segmentation of T 1 using the TRACLUS method (six sub-trajectories). (c) Segmentation of T 1 using the proposed method (three sub-trajectories). (d) The two most representative subtrajectories according to the proposed method.

12 Algorithm 1 :

 121 Global Voting Algorithm (GVA).

Fig. 5 :

 5 Fig. 5: The two sequential sliding windows W 1 (light gray horizontal lines) and W 2 (heavy gray vertical lines) locating at sample n on the given voting signal V k .

1 )

 1 if d(n) > τ and 2) d(n) ≥ d(i), i ∈ {nw, ..., n + w} (see lines 12-17 of Algorithm 2).

14 a = n 15 LP k = LP k + 1 18 Algorithm 2 :

 151182 Trajectory Segmentation Algorithm (TSA).

for k = 1 to n do 1 for i = 1 to LP k do 2 S 10 if 20 Algorithm 3 :

 11210203 k (i) = 0 //Empty sampling set S k (i) = 0 ∧ SR gain (k, i) > SR max gain then 11 SR max gain = SR gain (k, i) max gain > V P k (i) then break Subtrajectory Sampling Algorithm (SSA).

Fig. 6 :

 6 Fig. 6: The 2-D map of SMOD with the four onedirectional and two bidirectional roads.

Fig. 9 Fig. 7 :

 97 Fig. 7: The trajectories of our synthetic MOD (SMOD) with additive noise of SNR = 50 db projected in (a) 2-D spatial space ignoring time dimension and (b) spatiotemporal 3-D space. The trajectories of our synthetic MOD with additive noise of SNR = 30 db projected in (c) 2-D spatial space ignoring time dimension and (d) spatiotemporal 3-D space.

Fig. 8 :

 8 Fig. 8: The segmentation error for SMOD under additive noise of various SNR.

Fig. 9 :Fig. 10 :

 910 Fig. 9: Two trajectories of the synthetic MOD used for the segmentation method evaluation projected in spatiotemporal 3-D space (a) with additive noise of SN R = 50db, (b) with additive noise of SN R = 30db. (c), (d) The dotted lines shows the segmentation results projected on d(n) and on normalized voting signal V k for the trajectory (a) and (b), respectively.

Fig. 11 :

 11 Fig. 11: Results of sampling method for synthetic MOD with additive noise of SN R = 30db. (a) The function SR(S) under different sizes of sampling set. (b) The SR(S) gain as sampling set increases for synthetic MOD. (c), (d) The 60 most representative sub-trajectories are projected in spatiotemporal 3-D space using a colormap according to (c) SR(S) gain (d) line segments representativeness.

Fig. 12 :Fig. 13 :

 1213 Fig. 12: (a) The function SR(S) under different sizes of sampling set for "Athens trucks" MOD. (b) The SR(S) gain as sampling set increases for "Athens trucks" MOD.

Fig. 14 :

 14 Fig. 14: (a) The function SR(S) under different sizes of sampling set for "Milano-II" MOD. (b) The SR(S) gain as sampling set increases for "Milano-II" MOD.

Fig. 15 :

 15 Fig. 15: The trajectories of "Athens trucks" dataset (1.100 traj.) projected in (a) 2-D spatial space ignoring time dimension and (b) spatiotemporal 3-D space. (c), (d) The top 150 representative sub-trajectories of the dataset projected in 2-D spatial space using a colormap according to (c) SR(S) gain (d) sub-trajectories representativeness V P k (i).

Fig. 16 :

 16 Fig. 16: The trajectories of "Milano-II" dataset (11500 traj.) projected in (a) 2-D spatial space ignoring time dimension and (b) spatiotemporal 3-D space. (c), (d) The top 150 representative sub-trajectories of the dataset projected in 2-D spatial space using a colormap according to (c) SR(S) gain (d) sub-trajectories representativeness V P k (i).

Fig. 17 illustrates

 17 the function RM SE(S, D) under different sizes of sampling set for the synthetic MOD with additive noise of 50 db (see Fig.

17

 17

  (a)), the synthetic MOD with additive noise of 30 db (see Fig.

17

 17

  (b)), for "Athens trucks" MOD (see Fig.

17

 17

  (c)), for "Milano-I" MOD (see Fig.

17

 17

  (d)) and for "Milano-II" MOD (see Fig.

17

 17

  (e)). The blue, red dashed and black dashdot lines have been used to plot SSA, random and stratified sampling RMSE curves.

Fig. 17 :

 17 Fig. 17: The function RM SE(S, D) under different sizes of sampling set for (a) the synthetic MOD with additive noise of 50 db. (b) the synthetic MOD with additive noise of 30 db. (c) for "Athens trucks" MOD. (d) for "Milano-I" MOD. (e) for "Milano-II MOD.

  • C. Panagiotakis is with the Dept. of Commerce and Marketing, Tech.

	Educational Inst. of Crete and with the Dept. Of Computer Science,
	University of Crete, Greece. E-mail: cpanag@csd.uoc.gr
	• N. Pelekis is with the Dept. of Statistics and Insurance Science, University
	of Piraeus, Greece. E-mail: npelekis@unipi.gr

• I. Kopanakis is with the Dept. of Commerce and Marketing, Tech. Educational Inst. of Crete, Greece, Greece. E-mail: i.kopanakis@emark.teicrete.gr • E. Ramasso is with the Dept. of Automatic Control and Micro-Mechatronic Systems, FEMTO-ST Research Inst., France. E-mail: e ramasso@yahoo.fr • Y. Theodoridis is with the Dept. of Informatics, University of Piraeus, Greece. E-mail: ytheod@unipi.gr

TABLE 1 :

 1 Symbol table.

ACKNOWLEDGEMENTS

We thank anonymous reviewers for their very useful comments and suggestions. This research was partially supported by the FP7 ICT/FET Project MODAP (Mobility, Data Mining, and Privacy) funded by the European Union. URL: www.modap.org.


Probabilistic Methods

NMPC -Computational Approach

1: propagate model with a given set of inputs from k•T s .. k•T s +T p 2: J min,k  minimum cost function 3: if J min,k < J min,k-1 , then 4:

apply input with minimum cost 5: else 6:

apply previous input 7: end

8: k  k+1 9 www.ifr.uni-stuttgart.de             p T t t i i u i x F u J 0 0 , cost function Design Parameters  Control system  plant model, required underlying control loops  control signals -number -quantization -discretization  cost function  Implementation, i.e.
z f y f x f x E y E z E r q p ( E x, E y, E z) (f, q, y geodetic body-fixed ( f u, f v, f w) (p, q, r) www.ifr.uni-stuttgart.de NMPC LQR aircraft   ref z y x y q f , , , , ,   w v u r q p , , , , ,   y q f , , , ,
                                                                      r q p w v u c c c s s c c s c c s s c c c s s s s c c s s s s c c s c c s c s s c s s s c c c z y x f f f E E E q f q f f f q q f q q f q f q f q y q f y f y q f y f y q y q f y f y q f y f y q y q f 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0       d/
q p                                   ) cos( ) sin( ) cos( ) cos( ) sin( ) sin( 0 0 1 q f f q q f y q f    2 2 1 1 u u      a a x      2 1 2 1 , , , a a a a C              ) cos( ) cos( ) cos( ) sin( 0 ) sin( ) cos( 0 ) cos( ) sin( ) cos( ) cos( ) sin( ) sin( 1 q f q f f f q q f q q f C new direction =
 ) ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ( 1 i i i t i pq i i i i m i T i t o t k t c t c t c t c t x t Q t x u J         q f   t pq j k c j j , , , else 0 violated constraint if 0 q f       input constant with intervals time of #  m m i N i t t t i j j i     D     1 ,

  www.ifr.uni-stuttgart.de total flight time 24 h, 17 flights path length 260 km area 8 km 2 ground resolution < 2.6cm Motion Planning -Towards Autonomy  3 dimensional, considering 6 degrees of freedom, 12 states  From a given initial state ("in flight") to a final state ("landing")  Under constraints of -position (obstacles) -velocitiy (flight dynamics) -attitude (flight dynamics) -rates (flight dynamics) Motion Planning Research -Overview  Goerzen/Kong/Mettler, J Intell Robot Syst (2010) 57 A Survey of Motion Planning Algorithms from the Perspective of Autonomous UAV Guidance uncertainties and robustness has not been studied much..  ..majority of methods do not include much discussion about practical application.  When they do, they often only provide simulation results..  ..most of the practical implementations of UAV guidance have been of the hierarchical decoupled control type..

  .uni-stuttgart.de

  .uni-stuttgart.de  Prediction horizon T p  d quantization values of controls, simplest case: +1, 0, -1  m intervals with constant controls, simplest case: Dt 1 , Dt 2  n control signals, simplest case: p and q (roll/pitch rate)





  Prediction over 8 sec  Sample rate 10 Hz  d=81=3 2•2 control options  2 controls (p and q)  3 constant values each (+1/0/-1)  2 time intervals (Dt 1 +Dt 2 =T p )  within constraints, obstacle avoided www.ifr.uni-stuttgart.de  NMPC1  d=3, # of quantized control values  n = m = 2  d n•m = 81 control options  NMPC2  d=9, # of quantized control value  n = m = 2  d n•m = 6561 control options

de NMPC Method(s)

	 "Embedded" into control system
	 Kang&Hedrick, Linear Tracking of a Fixed-Wing UAV Using Nonlinear
	MPC, IEEE Trans. Contr. Syst. Technology, Sept 2009
	 2D (horizontal) motion
	 fixed wing
	 PC104
	 real-time testbench results
	 Rest of the presentation: NMPC with limited computational effort
	www.ifr.uni-stuttgart.de	8

 Waypoint RRT and bounding boxes  straight line and steady-state turns  RRTs with motion primitives  becomes complex with non-steady motion primitives  RRTs with closed loop dynamics included ("towards" NMPC) 7 www.ifr.uni-stuttgart.

www.ifr.uni-stuttgart.de

www.ifr.uni-stuttgart.de

More:

F. Weimer, M. Trittler, A. Joos, M. Gros, A. Posch, W.


Equations of Motion

The equations of motion are:

ẍ = 2 ẏ + x -(1 -µ) x -µ r 3 ps -µ x + 1 -µ r 3 pe + β 1 -µ r 2 ps r s , n 2 n x , ÿ = -2 ẋ + y - 1 -µ r 3 ps + µ r 3 pe y + β 1 -µ r 2 ps r s , n 2 n y , z = - 1 -µ r 3 ps + µ r 3 pe z + β 1 -µ r 2 ps r s , n 2 n z ,
where n = (n x , n y , n z ) is the normal to the surface of the sail with n x = cos(φ(x, y ) + α) cos(ψ(x, y , z) + δ), n y = sin(φ(x, y , z) + α) cos(ψ(x, y , z) + δ), n z = sin(ψ(x, y , z) + δ),

and r s = (x -µ, y , z)/r ps is the Sun -sail direction.

8 / 45

Equilibrium Points

The RTBP has 5 equilibrium points (L i , i = 1, . . . , 5). For small β, these 5 points are replaced by 5 continuous families of equilibria, parametrised by α and δ.

For a fixed small value of β, we have 5 disconnected family of equilibria near the classical L i .

For a fixed and larger β, these families merge into each other. We end up having two disconnected surfaces, S 1 and S 2 , where S 1 is like a sphere and S 2 is like a torus around the Sun.

All these families can be computed numerically by means of a continuation method.

Interesting Missions Applications

Observations of the Sun provide information of the geomagnetic storms, as in the Geostorm Warning Mission.

Sun

Interesting Missions Applications

To ensure reliable radio communication between Mars and Earth even when the planets are lined up at opposite sides of the Sun.

Periodic Motion Around Equilibria

We must add a constrain on the sail orientation to find bounded motion. One can see that when α = 0 and δ ∈ [-π/2, π/2] (i.e. only move the sail vertically w.r.t. the Sun -sail line):

The system is time reversible ∀δ by R : (x, y , z, ẋ, ẏ , ż, t) → (x, -y , z, -ẋ, ẏ , -ż, -t) and Hamiltonian only for δ = 0, ±π/2.

There are 5 disconnected families of equilibrium points parametrised by δ, we call them FL 1,...,5 (each one related to one of the Lagrangian points L 1,...,5 ).

Three of these families (FL 1,2,3 ) lie on the Y = 0 plane, and the linear behaviour around them is of the type saddle×centre×centre.

The other two families (FL 4,5 ) are close to L 4,5 , and the linear behaviour around them is of the type sink×sink×source or sink×source×source.

We focus on ...

We focus on the motion around the equilibrium on the FL 1 family close to SL 1 (they correspond to α = 0 and δ ≈ 0).

We fix β = 0.051689.

We consider the sail orientation to be fixed along time. Let us see the periodic motion around these points for a fixed sail orientation and show how it varies when we change, slightly, the sail orientation. Design station keeping strategy to maintain the trajectory of a solar sail close to an unstable equilibrium point.

Instead of using Control Theory Algorithms, we want to use Dynamical Systems Tools to find a station keeping algorithm for a Solar Sail.

We focus on the previous missions, where the equilibrium points are unstable with two real eigenvalues, λ 1 > 0, λ 2 < 0, and two pair of complex eigenvalues,

ν 1,2 ± i ω 1,2 , with |ν 1,2 | << |λ 1,2 |.
The linear dynamics at the equilibrium point is of the type saddle × centre × centre.

We describe the trajectory of the sail in three reference planes defined by the eigendirections.

(x1, y1) (x2, y2) (x3, y3)
For small variations of the sail orientation, the equilibrium point, eigenvalues and eigendirections have a small variation. We will describe the effects of the changes on the sail orientation on each of these three reference planes.

In the saddle projection of the trajectory:

When we are close to the equilibrium point, p 0 , the trajectory escapes along the unstable direction.

If we change the sail orientation the equilibrium point is shifted. Now the trajectory will escape along the new unstable direction.

We want to find a new sail orientation (α, δ) so that the trajectory will come close to the stable direction of p 0 .

In the centre projection of the trajectory:

A sequence of changes on the sail orientation implies a sequence of rotations around different equilibrium points on the centre projection. As we have seen a sequence of rotations around different equilibrium points can result unbounded. How can we choose the position of the new equilibrium point on the centre projection to keep this projection bounded ? The constants ε min , ε max and d will depend on the mission requirements and the dynamics around the equilibrium point.

Results

We have applied this station keeping strategy to two different mission applications, the Geostorm Warning Mission and the Polar Observer.

For each mission:

We have done a Monte Carlo simulation taking a 1000 random initial conditions.

For each simulation we have applied the station keeping strategy for 30 years.

We have tested the robustness of our strategy including random errors on the position and velocity determination, as well as on the orientation of the sail at each manoeuvre.

Results for the Geostorm (No Errors in Manoeuvres)

XY and XZ and XYZ Projections

Results

We manage to maintain the trajectory close to the equilibrium point for 30 years.

The most significant errors are the ones due to the sail orientation.

This station keeping strategy does not require previous planning as the decisions taken by the sail only depend on its position in the phase space.

The same ideas can be used to design strategies to move along the family of equilibrium points.

Surfing along the family of equilibria

Scheme on the idea to surf along the family of equilibria. The left hand side of equation (1),

L(x) = Dv (x)Ax -Bv (x),
is a linear operator w.r.t v (x) that diagonalizes if A and B are diagonal.

The right hand side of equation (1),

h(x) = g (x, v (x)) -Dv (x)f (x, v (x)), can be expressed as, h(x) = |k|≥2 h k x k
, where h k depend on v j in a known way.

It can be seen that for a fixed degree |k| = n, the h k depend only on the v j such that |j| < n.

Approximating the Centre Manifold (III)

We can solve equation ( 1) in an iterative way, equalising the left and the right hand side degree by degree.

Notice:

It is important to have a fast way to find the h k to compute high degrees.

We do not recommend to expand f (x, y ) y g (x, y ), and then compose with y = v (x). One should find other alternative ways, faster in terms of computational time.

For instance, compute the terms in a recurrent way and use the recurrences to derive an efficient algorithm.

The matrices A and B do not have to be diagonal, but then one must solve a larger linear system at each degree.

On the Centre Manifold

We have computed the centre manifold around different equilibrium points of the FL 1 family up to degree 16.

After this reduction we are in a four dimensional phase space (x 1 , x 2 , x 3 , x 4 ) that is difficult to visualise.

We need to perform suitable Poincaré sections to reduce the phase space dimension and help us visualise the phase space.

For δ = 0 the system is Hamiltonian and we have a first integral. We will use it to reduce the phase space dimension.

For δ = 0 the system is no longer Hamiltonian. We will take a function that varies little along the trajectories as an "approximate first integral", and use it to visualise the phase space.

Dynamics for δ = 0 (x 3 = 0 section)

  Earth's poles, as in the Polar Observer.

(

  Schematic representation of the equilibrium points on Y = 0)

  the Station Keeping Strategy (VI)To control the saddle and centre projection we want the new equilibrium point to satisfy:

  0001 -8e-05 -6e-05 -4e-05 -2e-0001 -8e-05 -6e-05 -4e-05 -2e-

/ 45

/ 45

v 5

Station Keeping around equilibria

Goal:

Design station keeping strategy to maintain the trajectory of a solar sail close to an unstable equilibrium point.

Instead of using Control Theory Algorithms, we want to use Dynamical Systems Tools to find a station keeping algorithm for a Solar Sail.

Idea:

We focus on the linear dynamics around an equilibrium point and study how this one varies when the sail orientation changes.

We want to change the sail orientation (i.e. the phase space) to make the system act in our favour: keep the trajectory close to a given equilibrium point.

Schematic Idea of the Station Keeping Strategy (I)

In the saddle projection of the trajectory:

When we are close to the equilibrium point, p 0 , the trajectory escapes along the unstable direction.

Schematic Idea of the Station Keeping Strategy (I)

In the saddle projection of the trajectory:

When we are close to the equilibrium point, p 0 , the trajectory escapes along the unstable direction.

If we change the sail orientation the equilibrium point is shifted. Now the trajectory will escape along the new unstable direction.

Schematic Idea of the Station Keeping Strategy (II)

In the saddle projection of the trajectory:

With these ideas we can control the instability due to the saddle. We need to take into account the centre projection of the trajectory, as it might grow.

Schematic Idea of the Station Keeping Strategy (III)

In the centre projection of the trajectory:

When we are close to the equilibrium point the trajectory is a rotation.

Schematic Idea of the Station Keeping Strategy (III)

In the centre projection of the trajectory:

When we are close to the equilibrium point the trajectory is a rotation.

If we change the sail orientation the equilibrium point is shifted. Now the trajectory will rotate around the new equilibrium point.

Schematic Idea of the Station Keeping Strategy (III)

In the centre projection of the trajectory:

When we are close to the equilibrium point the trajectory is a rotation.

If we change the sail orientation the equilibrium point is shifted. Now the trajectory will rotate around the new equilibrium point.

A sequence of changes on the sail orientation results in a sequence of rotations around the different equilibrium points.

Schematic Idea of the Station Keeping Strategy (V)

In the centre projection of the trajectory:

A sequence of changes on the sail orientation implies a sequence of rotations around different equilibrium points on the centre projection. As we have seen a sequence of rotations around different equilibrium points can result unbounded. How can we choose the position of the new equilibrium point on the centre projection to keep this projection bounded ?

Results for the Geostorm

We take β = 0.051689 (i.e. a satellite of 130kg mass with a 67m × 67m square sail).

Success

Reduction to the Centre Manifold

Using an appropriate linear transformation, the equations around the fixed point can be written as,

where A is an elliptic matrix and B an hyperbolic one, and f (0, 0) = g (0, 0) = 0 and Df (0, 0) = Dg (0, 0) = 0.

We want to obtain y = v (x), with v (0) = 0, Dv (0) = 0, the local expression of the centre manifold.

The flow restricted to the invariant manifold is

Approximating the Centre Manifold (I)

To find y = v (x) we substitute this expression on the differential equations. Hence, v (x) must satisfy,

We take,

its expansion as power series. Then we solve equation ( 1) to find the coefficients v k up to high degree (|k| = N).

v k x k is a high order approximation of the centre manifold.

gives a high order approximation of the motion on the centre manifold.

Dynamics for δ = 0

Here the first integral is:

We fix a Poincaré section x 3 = 0 to reduce the system to a three dimensional phase space. (Taking x 3 = 0 is like taking Z = 0).

We fix the energy level to determine x 4 and reduce the system to a two dimensional phase space that is easy to visualise. (Taking x 4 (J c , x) is like taking Ż (J c , x)).

We have taken several initial conditions and computed their successive images on the Poincaré section.

Dynamics for δ = 0

Here we take an "approximated first integral" :

We fix a Poincaré section x 3 = 0 to reduce the system to a three dimensional phase space. (Taking x 3 = 0 is similar to taking Z = Z * ).

We fix J c to determine x 4 and reduce the system to a two dimensional phase space that is easy to visualise. (Taking x 4 (J c , x) is like taking Ż (J c , x)).

We have taken several initial conditions and computed their successive images on the Poincaré section.

Station Keeping around Periodic Orbits

We want to extend the ideas on the station keeping around an equilibrium point for the station keeping around a periodic orbit.

Again, the idea is to change the phase space, by changing the sail orientation, so that the system acts as we wish, maintaining the sail's trajectory close to the nominal orbit.

We will focus on the motion around a Halo-type orbit. Although the same ideas should apply to other families of periodic orbits with the same characteristics.

First approach

The first idea is to extend the methods used for equilibrium points to periodic orbits.

First approach

The first idea is to extend the methods used for equilibrium points to periodic orbits.

There is an extra difficulty: while the points are persistent when the orientation of the sail is changed, periodic orbits may not. Hence, a careful study of the dynamics is needed.

Second approach

A second idea is the following: assume that the periodic (for instance, Halo) orbit we want to follow corresponds to α = α 0 and δ = δ 0 .

The computation of variational equations (of suitable order) w.r.t. α and δ gives explicit expressions for the effect of different orientations (close to the reference values α = α 0 , δ = δ 0 ) trajectory.

Hence, we can impose conditions on the "future" of the orbit and find orientations that fulfill them (or show that the condition is unattainable).