Feda Almuhisen

Nicolas Durand

Mohamed Quafafou

Detecting Behavior Types of Moving Object Trajectories

Keywords: Formal concepts, Frequent patterns, Behavior, Trajectories, Visualization

Trajectory mining is a challenging and crucial problem especially in the context of smart cities where many applications depend on human behaviors. In this paper, we characterize such behaviors by patterns, where each pattern type represents a particular behavior, e.g. emerging, latent, lost, etc. From GPS raw data, we introduce algorithms that allow computing a formal concept lattice which encodes optimal correspondences between hidden patterns and trajectories. In order to detect behaviors, we propose an algorithm that analyses the evolution of the discovered formal concepts over time. The method generates tagged city maps to easily visualize the resulting behaviors at different spatio-temporal granularity values. Refined or coarse analysis can thus be performed for a given situation. Experimental results using real-world GPS trajectory data show the relevance of the proposed method and the usefulness of the resulting tagged city maps.

Introduction

The concept of Smart City has become essential due to its concern of improving everyday technology in urban life, developing city network infrastructures, building better intelligent transportation systems to optimize traffic and transportation flows, and providing innovative services for citizens. In this context, a lot of positioning technologies, sensor networks, and online social media are used and deployed in cities. These devices and applications generate, by nature, spatiotemporal data which are collected from smartphones carried by people, GPS tracking systems on cars, airplanes, and location-based services provided by social media. Analyzing such data is very important to understand user behaviors for applications such as human mobility understanding, smart transportation, and urban planning [START_REF] Pan | Trace Analysis and Mining for Smart Cities: Issues, Methods, and Applications[END_REF]. Trajectory analysis [START_REF] Zheng | Computing with spatial trajectories[END_REF] is one of the mining tools that are used for this purpose. The analysis of trajectories is suited for a wide range of applications, for instance, monitoring, tracking, territorial management, and security. The handled spatio-temporal data are complex, rich and huge. They also inherit the geolocalization imprecision. Thus, they need to be analyzed by suitable methods based on clustering, classification, frequent patterns, etc. [START_REF] Li | Spatiotemporal Pattern Mining: Algorithms and Applications[END_REF]. Moreover, these methods generate reports and visualize the results for easily understanding.

In this paper, we propose a new method based on formal concept analysis to detect behaviors over time and visualize them on tagged city maps automatically generated. In order to achieve this goal, we segment and label the GPS trajectories according to the specified spatio-temporal granularity values. Then, each trajectory is transformed and represented in a georeference area as a set of grids. After that, this discrete representation of trajectory is used to extract the hidden patterns. From these data, we construct a formal concept lattice to encode the correspondence between the discovered frequent patterns and trajectories. Thus, a formal concept is a pair (Y, X), where X is a pattern and Y is the set of trajectories that pass through the grids belonging to X. Using this relationship between trajectories and grids, we define different types of behaviors (latent, emerging, decreasing, jumping, and lost) to statistically characterize the dynamics and diversity of mobility in relation to urban spaces. Finally, we use these results to automatically generate maps which are tagged by different colors, whereas each color represents a given behavior. These maps can be generated at different spatio-temporal level to refine the users mobility analysis in a given area of the city. We have evaluated our method using the T-drive data set, which is a GPS real taxi trajectories collected by Microsoft Research Asia [START_REF] Yuan | T-drive: Driving directions based on taxi trajectories[END_REF][START_REF] Yuan | Driving with knowledge from the physical world[END_REF]. Our method is generic regarding the application domain of trajectories and can be applied to smart city problems related to moving objects information, e.g., smartphones, bicycles, animals, etc.

The rest of the paper is organized as follows. Section 2 defines notations and basic notions necessary for understanding the paper. Related works are discussed in Section 3. The proposed method that allows an automatic tagging of city maps using patterns is detailed in Section 4. The experimental results are presented in Section 5. Finally, we conclude and present some future work in Section 6.

Background

A trajectory is a sequence of geographical locations of a moving object on a defined time period. A set of n moving object trajectories T rjs={T rj 1 , T rj 2 , ..., T rj n-1 , T rj n }, where each trajectory T rj i has a different length, and can be defined as a sequence of geographical points returned at a specific time stamp t [START_REF] Zheng | Trajectory Data Mining: An Overview[END_REF]. A trajectory T rj i is equal to {p 1 ,p 2 ,...,p k-1 ,p ki } and each point p k ∈ T rj i is defined in three-dimensional space, i.e., p k =(latitude k , longitude k , t k ), where (latitude k , longitude k ) defines the geographical coordinate location pair at specific time t k .

Frequent Patterns

The discovery of frequent patterns is an important task in data mining which has been discussed in [START_REF] Agrawal | Mining Association Rules between Sets of Items in Large Database[END_REF]. It corresponds into finding the sets of items (i.e. attribute values) which appear simultaneously in at least a certain given number of transactions (i.e. objects) recorded in a database. The latter is represented by the formal context D = (T , I, R), where T consists of a set of transactions, I is a set of items, and R ⊆ T ×I is a binary relation. Each pair (t, i) ∈ R means that transaction t ∈ T has a relation with item i ∈ I. A pattern is a subset of I. We used a string character notation for the patterns, for example AB for {A, B}. They are implicitly sorted according to the lexicographic order. The cover of a pattern X is the set of transactions of D that are related to the items of X (i.e., transactions that contain X):

A B D E D 1 t 2 A B D t 3 A B C D D 2 t 4 B C D t 5 A B C cover(X) = {t ∈ T | ∀i ∈ X, (t, i) ∈ R}.
The support of pattern X, i.e., support count, represents the number of transactions that contain X: support(X) = |cover(X)|.

Let us remark that we can also express the support as a percentage of the total number of transactions in D, i.e., |cover(X)| |T | ×100. Pattern X is frequent if its support exceeds (or is equal to) a minimum threshold denoted minsup. The set of all-frequent patterns is denoted FI and is equal to {X ⊆ I | support(X) ≥ minsup}. The notion of frequent closed pattern has been proposed by [START_REF] Pasquier | Efficient Mining of Association Rules Using Closed Itemset Lattices[END_REF]. The pattern X is closed if none of its supersets have the same support as X: ∀Y ⊃ X, support(Y ) < support(X). The frequent closed pattern is a pattern both closed and frequent. Let us note that the set of frequent closed patterns (denoted FCI) can be used to generate all the frequent patterns and their supports, and FCI ⊆ FI.

Example 1 Table 1 shows an example of a transactional database representing a formal context with 5 transactions and 5 items (denoted as A . . . E). If the minsup value is set to 3 (i.e., 60%) then pattern D is frequent because its support is 4 (t 1 , t 2 , t 3 and t 4 ). AE is not frequent since its support is 1 (t 1 ). We have FI={A, B, C, D, AB, AD, BC, BD, ABD}. C is frequent but not closed. BC is frequent and closed because it corresponds to the maximum set of items shared by t 3 , t 4 and t 5 . ABDE is closed but not frequent. We have FCI = {B, AB, BC, BD, ABD}.

Emerging Patterns

Emerging patterns have been introduced in [START_REF] Dong | Efficient Mining of Emerging Patterns: Discovering Trends and Differences[END_REF] as patterns whose support increased in a significant way from one dataset to another. The data represented in Table 2 contains 5 items (denoted A. . . E) and 5 transactions (denoted t 1 . . . t 5 ) distributed in 2 datasets (D 1 and D 2 ). Let us remark that these datasets can be seen as classes. The support of a pattern in a dataset D, denoted by supp(X, D), is the number of transactions which contain X in D. The quantitative evaluation of the contrast between data sets brought by a pattern is measured by its growth rate. The growth rate of a pattern X from D j to D i is defined by [START_REF] Agrawal | Mining Association Rules between Sets of Items in Large Database[END_REF].

GR j,i (X) = |D j | |D i | × supp(X, D i ) supp(X, D j ) . ( 1
)
The more GR j,i (X) is high, the more X characterizes D i compared to D j . Given a threshold value mingr > 1, a pattern X is an emerging pattern from

D j to D i if GR j,i (X) ≥ mingr. If GR j,i (X) = ∞, X is said to be a jumping emerging pattern (JEP) from D j to D i .
Example 2 Let us consider the example of Table 2 with mingr = 3 and verify if the D, ABD and ABDE are emerging from D 2 to D 1 or not. D is not an emerging pattern (GR 2,1 (D) = 1.5). ABD is an emerging

pattern because GR 2,1 (ABD) = 3. ABDE is a JEP (GR 2,1 (ABDE) = ∞).

Formal Concepts

Formal concept analysis (FCA) is a theory of data analysis identifying the conceptual structures within data sets. The closure properties and its capability of discovering inherent hierarchical structures give it an advantage to be used to analyze the different pattern relationships and build better mining algorithms [START_REF] Jay | A data mining approach for grouping and analyzing trajectories of care using claim data: the example of breast cancer[END_REF][START_REF] Poelmans | Formal concept analysis in knowledge processing: A survey on applications[END_REF]. Given a formal context D, there is a unique ordered set which describes the inherent lattice structure defining natural groupings and relationships among the transactions and their related items. This structure is known as a concept lattice or Galois lattice [START_REF] Ganter | Formal concept analysis: mathematical foundations[END_REF]. Each element of the lattice is a couple (T, I) which consists of a set of transactions (i.e., the extent) and a set of items (i.e., the intent). Each couple (called formal concept) must be a complete couple with respect to R, which means that the following mappings (noted f and g) hold. For T ⊆ T and I ⊆ I, we have: [START_REF] Agrawal | Mining Association Rules between Sets of Items in Large Database[END_REF]

f (T ) = {i ∈ I|∀t ∈ T , (t, i ) ∈ R} and (2) g(I) = {t ∈ T |∀i ∈ I , (t, i ) ∈ R}. f (
T ) returns items common to all transactions t ∈ T , while g(I) returns transactions that have at least all items i ∈ I. The idea of maximally extending the sets is formalized by the mathematical notion of closure in ordered sets. The operators h 1 = f •g and h 2 = g •f are the Galois closure operators. Let X be a pattern of items, if h 1 (X)=X, then X is a closed pattern. A formal concept is composed of a closed pattern of items and of the set of transactions containing this closed pattern. The concept lattice, denoted by L, is the set of formal concepts organized by the subsumption relation ≤.

Given two formal concepts (T 1 , I 1 ) and (T 2 , I 2 ), (T 1 , I 1 )

≤ (T 2 , I 2 ) ⇐⇒ T 1 ⊆ T 2 et I 2 ⊆ I 1 . (T 1 , I 1 ) (respectively (T 2 , I 2 )
) is a subconcept (respectively superconcept) of (T 2 , I 2 ) (respectively (T 1 , I 1 )). The frequent concept lattice is formed using the formal concepts that have at least minsup transactions in their extent.

Example 3

The line diagram in Figure 1 is an example of visual representation for the concept lattice computed from Table 1. Each node corresponds to a formal concept. The nodes upwards in this diagram are labeled by the attributes used to name the concept, and the nodes downwards are labeled by the objects used to name the concept. Note that the information reading for this diagram can be as follow: an object o is described by an attribute m if and only if there is an ascending path from the node named o to the node named by m. For example, the node labeled by t 5 corresponds to the formal concept ({t 3 , t 5 }, ABC). The node labeled by D represents ({t 1 , t 2 , t 3 , t 4 }, BD).

In this paper, we propose a new method to detect and visualize user behavior by analyzing trajectories over time. The originality of our method is to leverage formal concepts, frequent patterns and emerging patterns. Before presenting our method, we present the related work in the next section.

Related Work

In this section we first review some of the existing work about trajectory analysis by focusing on pattern mining methods. Then, in order to encode the correspondence between the discovered patterns and trajectories, we study existing work on formal concepts in the context of trajectory analysis. Finally, we review some work related to emerging patterns which capture trends over time, and are considered as particular behaviors.

Trajectory Analysis

Recently, a lot of studies focused on GPS trajectory analysis and providing extensive details about it [START_REF] Zheng | Computing with spatial trajectories[END_REF][START_REF] Zheng | Trajectory Data Mining: An Overview[END_REF]. Other works [START_REF] Hu | A comprehensive survey of recommendation system based on taxi gps trajectory[END_REF][START_REF] Bao | Recommendations in location-based social networks: a survey[END_REF] benefit from trajectory analysis and mining techniques to build better social networks and travel recommendation systems. Authors in [START_REF] Lin | Mining GPS data for mobility patterns: A survey[END_REF] provided a survey that discusses recent results on mining mobility patterns from GPS data. In [START_REF] Chen | Research on patterns mining method for moving objects[END_REF], the authors give a general view for recent publications that discuss trajectory analysis methods and patterns mining issues, by reviewing the methods and algorithms recently used in this field in detail and comparing their different approaches. From literature, trajectory analysis methods can be summarized as clustering method, which concern to find a group of trajectories that share common sub-traces to measure the closeness using distance measures between trajectories, then use density-based clustering methods in order to cluster the trajectories [START_REF] Lee | 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]. In [START_REF] Giannotti | Trajectory pattern mining[END_REF], the authors defined T-patterns (Trajectory patterns) as concise descriptions of frequent behaviors, in terms of both space (i.e., the regions of space visited during movements) and time (i.e., the duration of movements). They also proposed a method for extracting them using a sequential pattern mining algorithm. Alternatively, there is a growing interest in using semantic trajectory analysis techniques and methods which give a meaning to moving objects trajectory in terms of why, how, and for which purpose they move. [START_REF] Parent | Semantic trajectories modeling and analysis[END_REF] provided a survey that focuses on this term. Many tools and frameworks have also been developed to help in analyzing trajectory. An interesting tool, called M oveM ine introduced in [START_REF] Li | Movemine: Mining moving object databases[END_REF], integrates functions that include moving objects pattern mining and trajectory mining.

In the literature there are many work proposed extracting trajectory patterns using spatial information and the aim differs according to the research goal. In [START_REF] Lee | Mining frequent trajectory patterns in spatial-temporal databases[END_REF], the authors used a graph-based mining algorithm to extract the frequent trajectory patterns in a spatialtemporal database by scanning the database to generate a mapping graph and trajectory information list. Then, they use a depth first algorithm to mine frequent patterns. We can note that this method outperforms two classes of algorithms: Apriori-based and PrefixSpan algorithms. In another way, the authors of [START_REF] Benkert | Reporting flock patterns[END_REF] discussed mining flock patterns which are a type of spatio-temporal patterns that represent moving object groups close to each other in a certain time segment. They discussed algorithms related to extract this type of patterns [START_REF] Geng | Trajectory pattern mining in practice -algorithms for mining flock patterns from trajectories[END_REF] and proposed a modified algorithm called F P M . However, the limitation in such research was that it depends on the type of trajectory that we deal with as the users do not know which type of patterns they deal with. One of the works that has to be mentioned and solves this limitation is the work of [START_REF] Lee | A unifying framework of mining trajectory patterns of various temporal tightness[END_REF] where the authors proposed a unifying framework for extracting trajectory patterns from asynchronous movement. It has the potential of navigating the patterns at all different levels of temporal tightness based on location, time, or both and defines what they called a U T -pattern that helps in understanding the interaction between moving objects and group dynamics. In another way, in [START_REF] Cesario | Trajectory Data Analysis Over a Cloud-Based Framework for Smart City Analytics[END_REF], the authors mine the frequent user's behavior from the socio-geographic data by finding frequent regions and change its representation from movements between points into movement between dense regions. Then, they extract trajectory patterns from these dense regions in the form of associative rules. It is important to mention that most of the previous discussed research use clustering and pattern discovery methods to extract frequent behavior, while in our method we extract the proposed behavior types by leveraging formal concept analysis and emerging pattern discovery.

In order to discretize the trajectory data and to establish areas, i.e., regions, we use the grid partitioning concept. Several studies have been conducted to process spatial data using grids. For example, in [START_REF] Li | A trajectory data clustering method based on dynamic grid density[END_REF], the authors proposed a trajectory data clustering method based on dynamic grid division-merger limit. They map the trajectories into equal sized squares determined by experimental analysis. The authors of [START_REF] Xue | Solving the data sparsity problem in destination prediction[END_REF] divided the map of a city into uniform grids and compare their partitioning method within other methods. They conclude that the best prediction result is given by k-d tree based partitioning strategy due to the even distribution of points. Furthermore, the authors of [START_REF] Zhang | Sparse trajectory prediction based on multiple entropy measures[END_REF] proposed a spatial iterative grid partition algorithm. The space is partitioned into uniform grids multiple times recursively until a desired grid granularity is reached, resulting in a more balanced number of points in each cell and better predictions. Unlike the previous work, we use uniform grids because we perform a spatio-temporal analysis by extracting the frequent closed patterns to detect changes and evolution between two time periods.

Finally, the paper [START_REF] Yuan | Discovering Urban Functional Zones Using Latent Activity Trajectories[END_REF] introduces a new concept called latent activity trajectory which captures different urban socio economic activities in different locations and a chronological order. The authors segment the regions according to the major roads, then cluster them into functional zones based on a topic modeling approach. More precisely, they used a kernel density estimation to identify the intensity for each functional zone. Thus, only latent behaviors are discussed. Our work defines new types of behaviors, e.g., emerging, decreasing, lost, etc.

Formal Concepts and Trajectories

Formal concept analysis has been used in many applications and disciplines, such as knowledge discovery, information retrieval, etc. In [START_REF] Poelmans | Formal concept analysis in knowledge processing: A survey on applications[END_REF], the authors give an extensive survey of all the papers that apply FCA in different domains such as biology, chemistry, etc. Our concern in these works is trajectory analysis. In [START_REF] Buzmakov | Fca and pattern structures for mining care trajectories[END_REF], the authors provide a framework for analyzing sequential data structure based on FCA and they applied it on a data set of patient trajectories in a health-care system. The work in [START_REF] Jay | A data mining approach for grouping and analyzing trajectories of care using claim data: the example of breast cancer[END_REF] also presented an approach for clustering trajectories of care for breast cancer using FCA. However, this work does not concern moving object trajectory analysis.

To the best of our knowledge, there is no research work that is based on FCA for analyzing moving objects. Our proposition is based on FCA in order to capture the maximal relations between trajectories and regions that they pass through.

Emerging Patterns and Applications

Emerging patterns have different applications. They are highly used for building classifiers. In [START_REF] Li | Instance-based classification by emerging patterns[END_REF], the authors proposed the DeEP classifier, which is an instancebased classifier using the concept of emerging patterns. Researchers try to define new notions related to emerging patterns. For example, in [START_REF] Soulet | Condensed Representation of Emerging Patterns[END_REF], the authors define an exact condensed representation for the emerging patterns in a data set. They proposed their method of extracting emerging patterns that has the best growth rates and they called them strong emerging patterns (SEPs). In [START_REF] Durand | Emerging Overlapping Clusters for Characterizing the Stage of Liver Fibrosis[END_REF], the authors proposed adopting closed emerging patterns (CEPs) for characterizing classes. This method has been applied on hepatitis data to characterize the stage of liver fibrosis. Authors in [START_REF] Asses | A hybrid classification approach based on fca and emerging patterns-an application for the classification of biological inhibitors[END_REF] presented a hybrid classification method in the context of FCA, where each formal context has been augmented according to a class of objects. The class information is used to characterize the concepts whose extents include objects of a single class which allows extracting closed JEP. This approach has been applied in chemistry for classifying biological inhibitors. Our method is more general in terms of using the class of objects to specify the type of patterns. We also define other types of patterns and study the corresponding moving object behaviors.

FCA-based Behavior Discovery

In this section we present a novel method of trajectory analysis that generates moving object behavior tagged maps over time. The proposed method is based on frequent concept extraction and pattern evolution type detection. It is composed of three main steps: (1) Spatiotemporal preprocessing, (2) Pattern and behaviors extraction, and (3) Visualization of geolocalized behaviors (see Figure 2). In Step 1, trajectories are segmented, labeled, and mapped on a raster area G using local georeference according to time and spatial granularity values specified by the user. These discretized data are then used to extract hidden patterns and detect behaviors in the next step. In Step 2, for each time windows ∆t (i.e., a pair of time granularity values (t i-1 , t i )), we compute the frequent concept lattice L according to the minimal support threshold value minsup. Each node of L, i.e., formal concept, encodes the correspondence between a pattern (set of grids) and the set of trajectories that pass through the grids of this pattern. Then, for each formal concept, we analyze the evolution of the pattern by computing an indicative value K. By given K, a minimal threshold value of emergence θ, and an error tolerance value , we detect the type of the pattern that belongs to a predefined set of types: {latent, emerging, decreasing, jumping, lost}. Each pattern type corresponds to a behavior of trajectories. Let us note that Steps 1 and 2 correspond to the off-line process (i.e., precomputed). Step 3 corresponds to the online process, where for user specified time windows ∆t, we use the previous results to automatically generate maps which are tagged using colors, whereas each color represents a given behavior. Table 3 summarizes the important symbols.

Spatio-Temporal Preprocessing

Preprocessing the spatial trajectories is an important step in our methodology. It is performed on raw trajectory data to overcome the size problem. It is also needed to simplify and reduce the number of points with negligible error. These reduced trajectories are used later for improving the efficiency of the visualization process, also for the purpose of mining and analyzing trajectories to extract hidden patterns. In addition to that, the preprocessing step includes detecting and removing incomplete, inaccurate, incorrect, and out of range data points to generate vectors for the moving objects. It is also important to detect a certain time if this object stopped or moved. Therefore, stopping detection can be defined as, the observing of where an object spends a relatively large time without leaving a spatial area,

F C = (Y, X) A formal concept K i (Y, X)
Evolution indicator value of X between time i-1 and time i which is detected by having a sequential redundant spatial position on the generated vector over large period of time. A redundancy of spatial location (i.e., grid index, see Section 4.1.2) means either that the object is not moving or it is moving through the spatial defined coordinates. Thus, in this case we normalize the redundant spatial grid by one grid index in the generated vector. Our main motivation behind the elimination of the redundant grid index is that, in most cases, analysts are interested in the actual movement of an object. In addition, it will affect the pattern extraction later on.

The main steps to process the raw data are (1) Segmentation and labeling trajectories according to time granularity and (2) Mapping trajectories on a raster area according to spatial granularity.

Segment and Label Trajectories according to Time Granularity

Generally, GPS trajectory data are returned every 5-10 seconds and saved in servers. The problem is that the data to be analyzed are huge. So, we perform a trajectory segmentation by predefined time granularity, and we label each segment by its corresponding time class.

The time granularity value is chosen in order to segment each trajectory on set T rjs. It can be for instance minutes, hours, days, year, etc. Each value of granularity is called a class. We note C = {c 1 , c 2 , ...., c k-1 , c k } the set of classes values where c j < c j+1 for 1 ≤ j ≤ k -1. Using the specified time granularity, each trajectory is segmented according to the defined classes. The trajectory of object n after segmentation by C can be defined as follow:

T rj = {T rj c1 n , T rj c2 n , ..., T rj c k n }. If T rj cj n is empty then this trajectory is deleted.
Example 4 Assume that we have n = 100 moving objects in trajectory data T rjs = {T rj 1 , T rj 2 , T rj 3 , . . . , T rj 100 } collected for a duration of a month, for example in October. The time granularity value is the day. The trajectories by class obtained will be T rj cj n where (1 ≤ n ≤ 100) and 1 ≤ j ≤ 31. For instance, the class c 3 corresponds to the 3rd of October, which means that we keep the natural order of time sequence between the days of the month.

Mapping Trajectories on a Raster Area according to Spatial Granularity

The conversion of the geographic coordinates (latitude, longitude) values for the collected trajectories to an image of trajectory area map scale grid points can be achieved by applying two steps: (1) Generate a local georeference system area G for the purpose of mapping (see Algorithm 1), and (2) Map each trajectory to G (see Algorithm 2). Figure 3 illustrates the generated area. LatDp i ← latp i -latitude min ; // Latitude distance in degree from the minimum point in raster area to latp i 6:

LongDp i ← longp i -longitude min ; // Longitude distance in degree from the minimum point in raster area to longp i 7:

g i ← LatDp i /lat step * G.
Columns + LongDp i /long step + 1 ; // Grid index in the raster area for the corresponding geographical point 8:

G trj ← G trj ∪ g i ; 9: end for 10: return G trj ;

Algorithm 1 has the following input values: two points (point min and point max ) that cover the geographical area, we want to map (these points are diagonal and can be extracted from the trajectory real data set), also a spatial granularity width (G width ) and height (G height ) for the grid dimension in kilometers. Lines 2 and 3, calculate the latitude and longitude degrees that cover the mapping area. A conversion of degrees into kilometer map scale is performed in lines 4 and 5 in order to divide the raster area by applying general calculation for distance measure above or below the equator. For example, the general calculation for distance measure in kilometers for country like China which is above east the equator as shown in Figure 3 are as follow: Latk = LatD * 111.320; Longk = 111.320 * LongD * cos((LatD * π)/180). Finally, lines 6 and 7 calculate the number of columns and rows in the generated raster area according to a given grid dimension (spatial granularity) in kilometers, width and height. The algorithm returns the detail of the G raster area (Columns * Rows) cells and its position.

Algorithm 2 shows the way to map a trajectory that contains geographical points into spatial grids in the generated raster area. As input, we need a point min which is the minimum point in the generated raster area, and also a trajectory T rj which we want to map. First, the algorithm calculates the grid step distance long step and lat step value in degree for the covered lon- gitude and latitude in the raster area (lines 1 and 2). Assuming that (lat pi , long pi ) is a point from trajectory T rj, lines 5 and 6 show how to calculate the distances (in degree) that cover the area of mapping, based on the minimum point in the original local georeference discussed above. Then, line 7 calculates the grid index based on the columns and rows number of the whole raster area. The same steps are applied for each point of the trajectory. Finally, the algorithm returns the set of grid indexes G trj .

Example 5 Let us assume that we want to generate a local georeference for the following inputs: point min = (36.2355, 106.65184), point max = (36.2380153, 106.653277), and a spatial granularity G height = G width = 0.04km. Using Algorithm 1, we obtain 4 columns and 7 rows for the raster area as shown in Figure 4. Now, let us consider 4 points from T rj 1 as shown in Figure 4: p 1 = (36.23639831, 106.6520197), p 2 = (36.2364881, 106.6524688), p 3 = (36.2366678, 106.6526485), and p 4 = (36.237117, 106.652738). In order to return the corresponding grid indexes for these points we follow the steps in Algorithm 2. Let us take p 1 and p 2 . We know that G.Columns = 4, G.Rows = 7, G.LongD = 0.001437 and G.LatD = 0.0025153, which are previously calcu- lated using Algorithm 1. From these values, we calculate lat step , long step which are the values of distance in degree according to the spatial granularity given in our georeference system. Here, lat step = 0.000359324 and long step = 0.000359329. Then, we calculate the latitude and longitude degrees that cover the area corresponding to the given points p 1 and p 2 : LatD p1 =0.000898311, LongD p1 =0.000179664, g 1 = 9 ; LatD p2 =0.000898311, LongD p2 =0.000538993, g 2 = 10. By applying the same procedure to the rest of the points we will get the following values: g 3 = 15 and g 4 = 19. The trajectory will be defined according to the generated grid indexes as follow T rj = {9, 10, 15, 19}. Let us note that in case of having redundancy occurrence of the grid index in the same trajectory sequentially, this means that either the object is stopped or its moving in the same spatial granularity coordinates. In this case, we keep only one of the redundant grid index. Figure 5 and Table 4 show an example of 5 trajectories after mapping using the local georeference system, within its corresponding class time in two specified labels as {Sunday, M onday}.

Patterns and Behaviors Extraction

The behaviors detection associated to the trajectories over the time windows is performed in two main steps:

(1) Extraction of the frequent concept lattice and (2) Detection of behavior type for each frequent formal concept. Algorithm 3 summarizes these steps. The next two sections discuss in detail the algorithms.

Frequent Concept Lattice Computation

The computation of the frequent concept lattice (see section 2.3) allows to obtain both the grid closed patterns and their corresponding sets of trajectories. The first step to compute this lattice is to create the formal context D = (T , G, R) for a specified time windows ∆t (see Algorithm 3, line 1). T is the set of trajectories for the studied time window. I a set of grid index

G = m i=1 G i .
For each trajectory t and each grid index g ∈ G t , (t, g) is added to R to indicate that the trajectory t has a relation with grid index g. Then, the frequent concept lattice is computed according to the specified minimal threshold minsup (see Algorithm 3, line 2). Let us remark that CHARM L algorithm [START_REF] Zaki | Efficient algorithms for mining closed itemsets and their lattice structure[END_REF] coming from the pattern community is especially adapted for this computation. 5 presents the formal concepts of the frequent concept lattice obtained with minsup = 2 (i.e., 40%). For example, the formal concept (trj 1 , trj 2 , trj 3 , trj 5 , ; 9 10) means that the trajectories trj 1 , trj 2 , trj 3 , and trj 5 all pass through the grid indexes 9 and 10 (9 10 is a frequent closed pattern).

Detection of Behavior Types

After the computation of the frequent concept lattice for the specified time window, we detect the type of the closed pattern (i.e., the intent) of each frequent formal concept by studying the time classes of the corresponding trajectories (i.e., the extent). The type of a closed pattern reflects the behavior of users who have generated the corresponding trajectories. We define five types: emerging, decreasing, latent, jumping and lost. During the studied time window, the emerging type means that the presence of the pattern increased in the trajectories. The decreasing type means that the presence of the pattern decreased in the trajectories. The latent type means that the presence of the pattern is quite similar. The jumping type means the pattern which is absent, appeared. The lost type means that the pattern disappeared.

Let F C = (Y, X) be a formal concept, θ be the minimal threshold value of emergence, and be the error tolerance. For the time window [t i-1 , t i ], the type of the grid index pattern X (i.e., the intent of F C) is detected by computing K i (see Equations 2 and 3).

K i (Y, X) = count(Y, t i ) count(Y, t i-1 ) (2)
Where count(Y, t j ) is the number of trajectories of Y (i.e., the cardinality if the extent of F C) labeled by the class corresponding to t j .

Type =                                      If K i (Y, X) = 0 ± then the type of X is lost. (If ((K i (Y, X) > θ) ∧ (θ = 1)) ∨ ((K i (Y, X) ≥ θ) ∧ (θ > 1))) then the type of X is emerging. If K i (Y, X) < θ then the type of X is decreasing. If K i (Y, X) = 1 ± then the type of X is latent. If K i (Y, X) = +∞
then the type of X is jumping.

(3)

Example 7 Table 5 presents the result of computing K i and of the behavior type detection for the frequent formal concepts obtained from the example of Table 4 with minsup = 2 (i.e., 40%), θ = 1 and = 0. Here, t i corresponds to M onday and t i-1 to Sunday. Let us consider the first frequent formal concept (trj 1 ,trj 2 ,trj 3 , trj 5 ; 9 10). Two trajectories are in the Sunday class (trj 1 and trj 2 ) and the other two are in the M onday class (see Table 4). Thus, K i is equal to 1 and the type of the pattern {9 10} is latent. There are two latent patterns ({9 10} and {10 15}), three jumping patterns ({10 11}, {9 10 11} and {10 11 15}), one emerging pattern ({10}), and one decreasing pattern ({9 10 15}).

Visualization of Geolocalized Behaviors

This section explains the method that we proposed in order to visualize the behaviors identified for a specified time window (see Section 4.2). The visualization method is based on the spatial granularity value in kilometers, the number of columns and rows in the raster area generated in Algorithm 1. Algorithm 4 shows the steps that we follow in order to extract the geographical dimensions for certain grid patterns of the formal concept (see Section 4.2), which we will use later for the visualization method. Let f c ty be a frequent formal concept extended by its type, (intent f c , extent f c , T ype). In order to start the visualization method we need to extract the real dimension for each grid in the formal concept. Each grid g in the intent f c is defined by min, max dimensions based on the raster area. Then each detected type is assigned to a specific color. Hence, when reading the grid index and the type, the color is given to the mapped grid dimension area in real map. It is important to note that the type of each grid item is chosen according to the majority vote value, which is calculated by considering all the formal concepts that contains this item in its intent. In the case of equality, several strategies can be used. For instance, we can randomly choose from the obtained types or we can give a special type color. Let us remark that we can use the neighbors of the grid item to pick the type according to the majority vote, but we may also face the equality problem within the item neighbors. We chose to use a special color to indicate ambiguous type.

Example 8

In order to compute the grid dimensions for the concept (trj 1 , trj 2 , trj 3 , trj 5 ; 9 10) in Table 4 in the raster area of Figure 5, we do: Grid 9 : row g = 3, col g = 1. The (lat g min , long g min ) , (lat g max , long

Experiments

In this section, we summarize the experimental results obtained on a real world data set. We present information about this data and about the followed protocol. The results and discussion are then exposed.

Data

We have used a real data set (T-drive) which contains a large amount of taxi GPS trajectories collected by Microsoft Research Asia. These data have been used in [START_REF] Yuan | Driving with knowledge from the physical world[END_REF][START_REF] Yuan | T-drive: Driving directions based on taxi trajectories[END_REF]. The original data set contains the GPS trajectories for 10,357 taxis during the period of one week from Feb2 to Feb8 in 2008. Each file of the data set contains the trajectories of one taxi. Each line in a file has the following fields separated by comma: TaxiID, Date, Time, Longitude and Latitude. The total number of points in this data set is about 15 million and the total distance of the trajectories reaches 9 million kilometers. The (longitude,latitude) was returned every 5 to 10 seconds and the coordinate system used for the returned longitude latitude was the Decimal Degree coordinate system (DD). In order to perform our experiments, we choose 222 taxis randomly.

Experimental Protocol

For the experiments, the methodology has been applied based on three factors, time granularity, minsup value, and spatial granularity. The protocol in general can be summarized as follows, the trajectories have been segmented according to time granularity, all the inaccurate, redundant data have been deleted, then the trajectories were mapped in raster area G by applying Algorithm 1 using the spatial granularity value. Each trajectory has been mapped on the raster area for finally having a vector which contains information about the grid numbers that each trajectory passes through its time window journey (see Algorithm 2). For each pair of time ∆t, and for certain minsup value, we extracted the frequent formal lattice using CHARM L algorithm [START_REF] Zaki | Efficient algorithms for mining closed itemsets and their lattice structure[END_REF]. Then, the type of the pattern contained in each formal concept has been identified (see Algorithm 3). Let us recall that the type of a pattern corresponds to a certain behavior. Finally, we visualized the moving object behavior on a real map by applying Algorithm 4.

Several experiments were applied by changing time granularity (24 hours, 12 hours, 2 hours), spatial granularity (20, 40, 60 meters), and minsup [START_REF] Durand | Frequent Itemset Border Approximation by Dualization[END_REF][START_REF] Lee | Trajectory clustering: A partition-and-group framework[END_REF][START_REF] Pasquier | Efficient Mining of Association Rules Using Closed Itemset Lattices[END_REF][START_REF] Zheng | Computing with spatial trajectories[END_REF]. For the type detection we set to 0 and θ to 1 in Algorithm 3.

In order to show the relevance and the interest of our method, we compare our results with those of a baseline method which uses only statistics instead of patterns. The baseline method computes, for each grid, the total number of the trajectories passing through it. The type detection and the visualization are performed as our method.

All the algorithms were implemented in Java, and all the experiments were performed on a Intel Core i7 2.50 GHz with 16 GB of memory.

Experimental Results

As mentioned before, a study according to different spatial granularities and zones within different pairs of days or half a day and different minsup values has been done to evaluate the proposed method. We present the quantitative and qualitative results of our method, and the comparison with the baseline method. Let us note that the corresponding days in the data are as follow: Feb2 = Saturday, Feb3 = Sunday, Feb4 = Monday, Feb5 = Tuesday, Feb6 = Wednesday, Feb7= Thursday, and Feb8 = Friday.

Quantitative Results

Tables 6,7, and 8 and Figures 7, 6 and8 show the statistical and graphical representation for our experimental results using different spatial granularities 20, 40, 60 meters respectively, for minsup = 10 (i.e, 2.3%). We observe that when we have high spatial granularity dimensions, we have a high number of patterns. This is clearly shown in Figure 9 where the total number of patterns increases whenever we have a high granularity value. Thus, if the spatial granularity value increases, the total number of grids decreases and the data density increases. Consequently, there are more possible common grids between the trajectories, so the probability to find patterns increases. We can conclude that there is a relation between spatial granularity and the total number of generated patterns. Moreover, we observe that the number of lost and jumping patterns is low compared to those of the other types. This is due to the time granularity which is coarse and does not allow to capture fast evolution of behaviors. A pattern is lost or jumping if it is not detected during both the two compared time intervals.

In order to study the effect of minsup values and types, we did another experiment where we change the minsup value for same spatial granularity for one week as shown in Tables 9, 10, and 11. We note from these statistics that when ever the minsup value increases the number of generated pattern decrease and some types may disappear. For example, if we compare the results in Table 8 with those in Table 9, we can notice that whithin minsup=20 (i.e., 4.5%) we have no jumping or lost patterns.

We have also applied the same protocol by changing the time segmentation. Instead of using 24 hours, we have set the time granularity value to 12 hours and 2 hours, while keeping spatial granularity value as 60 meters and minsup value as 2.3%. This percentage corresponds to 10 trajectories in the case of 12 hours (same as 24 hours), and 6 for 2 hours (because the total number of trajectories is lower). Tables 12 and13 show the statistical results for 12 hours and 6 hours, respectively. In Table 13, Feb3(0:00-2:00) means from 0:00:00 to 1:59:59. If we compare these statistics within the ones in Table 8, we can conclude that when ever we decrease the time granularity, we get less total number of patterns. For example, in time window (Feb3,Feb4) Table 8, the total number of patterns is equal to 88,719 while for the same period of time on Table 12 it is equal to 18,444. The decreasing of the number of patterns occurs because the data density is lower. Let us remark that decreasing the minsup leads to more patterns.

Qualitative Results

In order to analyze the results, we captured a different interesting areas shown in Figure 10 in a real map, and visualized the detected behaviors. We give explanation for the results. For the purpose of clarity, we give each pattern type a key letter in the figures where each letter related to a color and type of pattern as follows: (R:Red, decreasing, (G:Green, emerging), (B:Blue, latent), (Y:Yellow, lost), (GR:Gray, jumping and (P:Pink, ambiguous). It is important to say that the percentage of ambiguous patterns (see Figure 11) which are the grids that hold two types of behavior at the same time was very low. For example, the percentage of ambiguous type is 2.24% for (Feb2(pm),Feb3(am)). Some emerging patterns have been captured in the airport area (see Figure 12) and in the Administration area (see Figure 13) during the period of (Feb2,Feb3), and (Feb3,Feb4) respectively. This means that these regions face a higher taxi movement in the second day compared to the first day. The result in airport area can be clarified by the fact that Chinese New Year was in this week, where most of Chinese go back to celebrate the holiday. Taxi movement in administrative area was also higher in Feb4 compared to Feb3 and this can be explained by the fact that Feb4 is a working day while Feb3 is weekend. The majority of patterns shown in the temple of heaven area (see Figure 16) during (Feb6,Feb7) were also emerging, which means that in Feb7 more taxis passed this area comparing to the one in Feb6 where there was a worshiping ceremony at this temple for the Chinese New Year celebration.

In another way, some decreasing patterns have been identified in the Beijing railway station area during (Feb4,Feb5) (see Figure 14). We got the same observation for the Chaoyang golf club area during (Feb5,Feb6) (see Figure 15). It means that these areas have been facing a higher taxi movement in the first day compared to the one in the second day. For example, this behavior can be clarified in Beijing railway area by the fact that Feb4 is the first day of the week so the transportation used will be higher than any other day through the week.

Finally, a latent pattern has been discovered in the region of Temple of heaven during (Feb7,Feb8) (see Figure 17). It means that the taxis movement is approximately similar in these two days. This can be explained by the fact that in these two days temple of heaven (Tiantan park) which is the largest building for religious worshiping in China also a tourist attraction where having the chinese new year ceremonies.

To study the effect of time granularity, we compare Figure 18 which shows the behavior in the region of Beijing Railway during the period of (Feb3,Feb4) where time granularity equals day, within Figures 19 and20 where the time granularity value is 12 hours. We do not visualize the results obtained for 2 hours because the total number of discovered patterns is too low with no important conclusion.

It is clear that when ever we have a smaller time granularity value, we will get more homogeneous behaviors. For example, in Figure 18, the whole behavior of the area is not really clear compared to Figure 19 where the majority of the patterns are emerging. This means that the area face higher taxi movements afternoon compared to the morning period. This can also be verified by Figure 20. We can observe a decreasing behavior which confirms that the traffic on Feb3 pm was higher compared to the morning period of Feb3 and Feb4. We can explain that by, Feb3 is Sunday pm and all the citizens use transportation at the end of weekend or holiday.

Comparative Results

Tables 14, 15 and 16 present the distribution of the grids according to behavior types computed by the baseline method. We observe that the number of jumping and lost grids is very high compared to the number of the other types. These two types are dominant and do not allow to identify the others. Moreover, the total number of grids to visualize is high, which does not help to analyze the results. As we can see in Figure 21, a lot of grids are tagged as "lost" and the other behavior types are hidden. If the number of labeled grids is too high, the map can be impractical to compute and exploit.

In fact, our method avoids these problems as the discovery of frequent closed patterns, i.e., contained in formal concepts, allows us to tune the importance of the considered trajectories by using the minsup threshold value. The more the minsup value is high, the more the discovered patterns are important and in the same time their number is smaller. In conclusion, our method thus detect only the important behavior types and not all types.

Conclusion

Trajectory analysis is a challenging problem especially in the urban context as the world's population will in- crease significantly in the next few decades. The automatic generation of on demand reports, that aggregate results of trajectory analysis and the visualization of geolocalized information, is an important issue that will support urban management, adapted policy definition, public security, etc.

In this paper, we have proposed a new method based on frequent formal concept extraction and the analysis of their evolution over time in order to detect behaviors. We have proposed five types of behavior: emergent, latent, decreasing, jumping, and lost. Our method generates tagged city maps at different spatio-temporal gran-ularity values. Therefore, refined or coarse analysis can be performed for a given situation. Experimental results using real-world GPS trajectory data have shown the relevance of the proposed method and the usefulness of the resulted tagged maps.

In future work, we will extend our system by predicting user behavior in a given area during the next time window. We will also use our recent works on approximation of frequent itemset borders [START_REF] Durand | Frequent Itemset Border Approximation by Dualization[END_REF][START_REF] Durand | Approximation of Frequent Itemset Border by Computing Approximate Minimal Hypergraph Transversals[END_REF] and on probabilistic topic models [START_REF] Naim | Probabilistic Approach for Diversifying Web Services Discovery and Composition[END_REF][START_REF] Naim | Semantic Pattern Mining Based Web Service Recommendation[END_REF] for enrichment and trajectory recommendation purposes.

Fig. 1 :

 1 Fig. 1: Example of concept lattice.

Fig. 2 :

 2 Fig. 2: Overview of the proposed approach.

Fig. 3 :Algorithm 1 1 :

 311 Fig. 3: Georeference system area according to spatial granularity.

Algorithm 2

 2 Mapping a trajectory on a raster area according to spatial granularity Require: T rj : {p 1 , p 2 , ..., p k-1 , p k } // each point contains (latp k , longp k ) point min : (latitude min , longitude min ) // Minimum geographical coordinate in the raster area. G : Raster area Ensure: G trj : Set of grid indexes corresponding to T rj 1: long step ← G.LongD/G.Columns ; // The longitude distance between each column in raster area in degree 2: lat step ← G.LatD/G.Rows ; // The latitude distance between each row in raster area in degree 3: G trj ← {}; 4: for i = 1 to k do 5:

Fig. 4 :

 4 Fig. 4: Geo-mapped trajectory example.

Fig. 5 :

 5 Fig. 5: Example of mapping.

Algorithm 3 1 :

 31 Extraction of frequent concept lattice and detection of patterns type Require: ∆t: Time window G : l j=1 G j which is consists of grid indexes, where l = 1. . . n * m (m is the number of classes, and n is the number of initial trajectories) minsup: Minimum support threshold value θ: Minimal emergence threshold value : Error tolerance value Ensure: L: Frequent concept lattice St: Set of detected types for the formal concepts of L Creation of formal context D =(T , G, R) according to the time window ∆t; 2: Computation of frequent concept lattice L for D according to minsup; 3: Detection of the type of each pattern contained in the formal concepts of L according to θ and ; 4: return L and St ;

Example 6

 6 Let us consider the example of Table 4 as the binary relation R of a formal context D = (T , G, R) where T = {trj 1 , trj 2 , trj 3 , trj 4 , trj 5 }, G = {9, 10, 11, 15, 19}. The time granularity value is one day (Sunday, October 11, 2016 and Monday, October 12, 2016), and each trajectory is labeled by a time class C ∈ {Sunday, M onday}. Table

Fig. 6 :

 6 Fig. 6: Number of patterns for one week and spatial granularity = 20 meters.

Fig. 7 :

 7 Fig. 7: Number of patterns for one week and spatial granularity = 40 meters.

Fig. 8 :

 8 Fig. 8: Number of patterns for one week and spatial granularity = 60 meters.

Fig. 9 :

 9 Fig. 9: Total number of patterns for one week and different spatial granularities.

Fig. 10 :Fig. 12 :

 1012 Fig. 10: Beijing areas of interest. Fig. 11: Ambiguous pattern example.

Fig. 14 :

 14 Fig. 14: Tagged map for the Beijing Railway station area (Feb4,Feb5).

Fig. 15 :

 15 Fig. 15: Tagged map for the Chaoyang Golf club area (Feb5,Feb6).

Fig. 16 :

 16 Fig. 16: Tagged map for the Temple of heaven area (Feb6,Feb7).

Fig. 17 :Fig. 18 :

 1718 Fig. 17: Tagged map for the Temple of Heaven and Danbi Bridge area (Feb7,Feb8).

Fig. 19 :

 19 Fig. 19: Tagged map for Beijing Railway station area in (Feb3(am),Feb3(pm)).

Fig. 20 :

 20 Fig. 20: Tagged map for Beijing Railway station area in (Feb3(pm),Feb4(am)).

Table 1 :

 1 Example

			of transactional database
	tid			items
	t 1	A	B		D	E
	t 2	A	B		D
	t 3	A	B	C	D
	t 4		B	C	D
	t 5	A	B	C
	Table 2: Example of transactional data in different datasets
	tid		items		tsid
	t 1

Table 3 :

 3 Important symbols

	Symbol	Explanation
	T rjs	Set of trajectories
	trj	A trajectory
	p	Geographical point
	long ∆t	Time window (a pair of time granularity values)
	G trj	Set of grids corresponding to trj
	D	Formal context
	minsup	Minimum support threshold value
	θ	Minimal emergence threshold value
		Error tolerance value
	L	Frequent concept lattice
	X	A pattern of grid cells

p , lat p Longitude and latitude of p LongD Longitude and latitude distance LatD covered in degrees Longk Longitude and latitude distance Latk covered in kilometers C Set of time granularity values G Raster area G height Dimensions of one grid G width (according to spatial granularity)

Table 4 :

 4 Example of trajectory data presented within grid indexes and their corresponding classes (the time granularity is the day)

	Trajectory		Grid indexes (G)		Class (C)
	trj 1	9	10		15	19	Sunday
	trj 2	9	10		15		Sunday
	trj 3	9	10	11	15		M onday
	trj 4		10	11	15		M onday
	trj 5	9	10	11			M onday

Table 5 :

 5 Frequent formal concepts and type obtained for the example of Table4(minsup = 2 (40%), θ = 1, = 0)

		Formal Concept	#Sunday #Monday	K	Type
		( trj 1 trj 2 trj 3 trj 5 ; 9 10 )	2	2	1	LATENT
		( trj 1 trj 2 trj 3 trj 4 ; 10 15 )	2	2	1	LATENT
		( trj 3 trj 4 trj 5 ; 10 11 )	0	3	+∞	JUMPING
		( trj 3 trj 5 ; 9 10 11 )	0	2	+∞	JUMPING
		( trj 3 trj 4 ; 10 11 15 )	0	2	+∞	JUMPING
		( trj 1 trj 2 trj 3 trj 4 trj 5 ; 10)	2	3	1.5	EMERGING
		( trj 1 trj 2 trj 3 ; 9 10 15 )	2	1	0.5	DECREASING
	Algorithm 4 Grid geographical dimension extraction
	in raster area
	Require:
		f c ty : Formal concept containing three parts (intent f c ,
		extent f c ,T ype)
		G: Raster area contain (columns*rows) cells, also its lat-
		itude, longitude distance covered in degree (LatD and
		LongD)
	Ensure:
		LG the
		grid column position
	7:	lat g min ← G.lat min +(lat step * (row g -1)); //calculate
		the grid g latitude min value
	8:	lat g max ← lat g min + lat step ; //calculate the grid g lat-
		itude max value
	9:	long g min ← G.long min + (long step * (col g -1)); //cal-
		culate the grid g longitude min value
	10:	long g max ← long g min + long step ; //calculate the grid
		g longitude max value
	11:	add(LG dim , {(lat g min ,long g min ), (lat g max ,long g max )});
		// add the pair grid min max longitude, latitude values
		to the end of the list LG dim
	12: end for
	13: return LG dim ;

dim : Ordered list containing the grid dimension min, max for formal concept intent values 1: LG dim ← {}; 2: long step ← G.LongD/G.Columns; //longitude distance between each column in raster area in degree 3: lat step ← G.LatD/G.Rows; //latitude distance between each row in raster area in degree 4: for all g ∈ intent f c do 5:

row g ← (g -1)/G.Columns + 1; //calculate the grid row position 6:

col g ← ((g -1) mod G.Columns) + 1; //calculate

Table 6 :

 6 Number of patterns according to the type, one week, granularity 20 meters and minsup = 10 (2.3%)

	Time window Emerging Decreasing Latent Lost Jumping Total
	(Feb2,Feb3)	1,699	363	252	0	15	2,329
	(Feb3,Feb4)	2,810	2,011	846	1	3	5,671
	(Feb4,Feb5)	1,674	2,755	819	2	0	5,250
	(Feb5,Feb6)	323	1,594	238	4	0	2,159
	(Feb6,Feb7)	105	289	74	1	0	469
	(Feb7,Feb8)	7	95	18	0	0	120

Table 7 :

 7 Number of patterns according to the type, one week, granularity 40 meters and minsup = 10 (2.3%)

	Time window Emerging Decreasing Latent Lost Jumping	Total
	(Feb2,Feb3)	8,184	1,195	1,074	0	41	10,494
	(Feb3,Feb4)	11,439	9,072	3,642	4	14	24,171
	(Feb4,Feb5)	7,664	11,221	3,094	19	5	22,003
	(Feb5,Feb6)	1,422	7,605	1,107	21	0	10,155
	(Feb6,Feb7)	555	2,043	278	8	1	2,885
	(Feb7,Feb8)	109	685	88	0	0	882

Table 8 :

 8 Number of patterns according to the type, one week, granularity 60 meters and minsup = 10 (2.3%)

	Time window Emerging Decreasing Latent Lost Jumping	Total
	(Feb2,Feb3)	20,729	2,220	2,041	0	312	25,302
	(Feb3,Feb4)	40,442	34,797	13,431	9	40	88,719
	(Feb4,Feb5)	20,767	40,327	10,118	55	3	71,270
	(Feb5,Feb6)	2,658	15,668	1,933	90	0	20,349
	(Feb6,Feb7)	1,290	3,105	524	1	1	4,921
	(Feb7,Feb8)	237	1,421	148	1	0	1,807

Table 9 :

 9 Number of patterns according to the type, one week, granularity 60 meters and minsup = 20 (4.5%)

	Time window Emerging Decreasing Latent Lost Jumping Total
	(Feb2,Feb3)	1,678	99	57	0	0	1,834
	(Feb3,Feb4)	2,304	1,872	447	0	0	4,623
	(Feb4,Feb5)	1,328	2,516	444	0	0	4,288
	(Feb5,Feb6)	148	1,544	96	0	0	1,788
	(Feb6,Feb7)	115	300	30	0	0	445
	(Feb7,Feb8)	12	124	10	0	0	146

Table 10 :

 10 Number of patterns according to the type, one week, granularity 60 meters and minsup = 30 (6.8%)

	Time window Emerging Decreasing Latent Lost Jumping Total
	(Feb2,Feb3)	389	9	5	0	0	403
	(Feb3,Feb4)	587	422	75	0	0	1,084
	(Feb4,Feb5)	261	618	76	0	0	955
	(Feb5,Feb6)	26	303	10	0	0	339
	(Feb6,Feb7)	23	36	3	0	0	62
	(Feb7,Feb8)	0	13	0	0	0	13

Table 11 :

 11 Number of patterns according to the type, one week, granularity 60 meters and minsup = 40 (9%)

	Time window Emerging Decreasing Latent Lost Jumping Total
	(Feb2,Feb3)	103	1	0	0	0	104
	(Feb3,Feb4)	197	140	22	0	0	359
	(Feb4,Feb5)	69	207	18	0	0	294
	(Feb5,Feb6)	4	66	1	0	0	71
	(Feb6,Feb7)	2	4	0	0	0	6
	(Feb7,Feb8)	0	0	0	0	0	0

Table 12 :

 12 Number of patterns according to the type, one week, spatial-granularity 60 meters, time granularity 12 hours and minsup = 10 (2.3%)

	Time window	Emerging Decreasing Latent Lost Jumping Total
	Feb2(am),Feb2(pm)	0	0	0	0	797	797
	Feb2(pm),Feb3(am)	394	2,135	206	33	0	2,768
	Feb3(am),Feb3(pm)	4,276	267	171	0	128	4,842
	Feb3(pm),Feb4(am)	833	6,087	653	69	1	7,643
	Feb4(am),Feb4(pm)	4,557	796	515	2	89	5,959
	Feb4(pm),Feb5(am)	135	3,767	174	125	1	4,202
	Feb5(am),Feb5(pm)	4,339	144	161	2	48	4,694
	Feb5(pm),Feb6(am)	24	3,358	35	153	0	3,570
	Feb6(am),Feb6(pm)	78	889	66	5	0	1,038
	Feb6(pm),Feb7(am)	11	695	17	21	0	744
	Feb7(am),Feb7(pm)	356	24	15	0	12	407
	Feb7(pm),Feb8(am)	6	275	7	24	0	312
	Feb8(am),Feb8(pm)	63	6	8	0	1	78

Table 13 :

 13 Number of patterns according to the type, Feb3, spatial-granularity 60 meters, time granularity 2 hours and minsup = 6 (2.3%) Time window Emerging Decreasing Latent Lost Jumping Total Feb3

Acknowledgments

The authors would like to thank the anonymous reviewers for comments and suggestions on earlier drafts of this paper.


Denise Karamessini

Salomé Poyer

Laurence Charles
email: laurence.charles@univ-amu.fr

Jean-François Lutz
email: jflutz@unistra.fr

Conditions for Self-Archiving 2D sequence-coded oligourethane barcodes for plastic materials labeling

Mixtures of uniform sequence-defined oligourethanes were evaluated as 2D molecular barcodes for labeling three different commodity polymers, namely polystyrene, polyvinylchloride and polyethylene terephthalate. Six different oligourethanes were synthesized by solid-phase iterative synthesis and were coded using a binary monomer alphabet.High-resolution mass spectrometrystudies indicated that all oligomers are uniform and sequence-defined. However, instead of using them as individual codedchains, oligomers with different chain-length mass and sequence were mixed into intentionally polydisperse libraries. In particular, a three-component and a four-component library were created to encode a 2-bytes model binary sequence. These 2D-coded libraries were incorporated in all commodity plastics via a simple solvent casting procedure. Furthermore, in all cases, the oligomer mixtures could be extracted from the host polymer films and deciphered by mass spectrometry, thus opening interesting avenues for anti-counterfeiting and traceability applications.

Introduction

Among the potential applications of sequence-controlled polymers, [1] it has been suggested that anti-counterfeiting technologies could be an interesting area of use. [2] Indeed, a macromolecule containing a perfectly defined sequence of monomers can potentially be used as a label or "barcode", [3] in which molecular information certifying authenticity can be included. For instance, DNA, which contains a defined sequence of nucleotides, is sold by several companies as an anti-counterfeiting solution for the labeling of high-value materials.

Since DNA strands can be selectively hybridized and amplified by polymerase chain reaction, DNA tags can be used in trace amounts for materials labeling. However, DNA is not stable in all conditions and may degrade during materials processing, e.g. in hot-melt extrusion, thermoset curing, photo-curing, 3D printing or sol-gel process. Thus, for many applications, DNA barcodes have to be encapsulated, for example in processable nanoparticles as developped by Grass and coworkers. [4][5] Alternatively, non-natural sequence-defined polymers may also be used as anti-counterfeiting barcodes. [2] Over the last years, we [6][7][8][9][10][11] and others [12][13][14][15] have demonstrated that a wide variety of uniform sequence-defined synthetic macromolecules can be synthesized using multi-step growth approaches. [16] In particular, our group has evidenced that such polymers can be used to store information. For instance, binary information can be easily written in synthetic sequence-defined macromolecules using two monomers that are defined as 0 and 1-bits. [17] In addition, the digital sequences can be easily read using a sequencing methodology [18] such as tandem mass spectrometry (MS/MS). [19][20][21] Among all the different types of digital polymers that have been reported by our group, sequence-coded oligourethanes (or oligocarbamates) appear as a very interesting class of polymers for anti-counterfeiting technologies. Uniform sequence-defined oligocarbamates can be synthesized using a variety of solid-phase iterative approaches. [14,[22][23][24] In a publication reported last year, [25] our group has shown that digital oligourethanes can be synthesized using a simple protecting-group-free iterative approach. [26] Furthermore, it was found in this study that the COOH-terminated oligomers prepared via this orthogonal approach are particularly easy to decode by negative-mode MS/MS sequencing. [25] In other words, sequence-defined oligourethanes can easily be both coded and decoded. Furthermore, classical polyurethanes constitute a well-known class of plastics with widely-described physico-chemical properties, including for example chemical resistance and thermal stability. Therefore, sequence-coded oligourethanes appear as ideal candidates for materials labeling. For instance, we have already demonstrated that oligourethane barcodes can be included in small amounts (typically less than a weight percent) in casted polystyrene films, [25] 3D-printed photo-cured methacrylate resins [25] and acrylate-based intraocular implants. [27] In all cases, it was also demonstrated that the barcodes can be selectively extracted from the host polymer materials and sequenced by MS/MS. Yet, in all these examples, the information sequence was included in a single oligomer chain (i.e. 1D-coding). In another recent study, our group has shown that information can also be written in a set of oligomers of different mass (i.e. 2D-coding). [28] In such a strategy, digital information is stored in a library of chains. The electrospray mass spectrometry (ESI-MS) analysis of these mixtures leads to an intentionally-polydisperse spectrum (i.e. first dimension fingerprint) and the MS/MS sequencing of eachindividual peakof this mass distributiongives access to the sequence-information included in each oligomer (i.e. second dimension fingerprint). This concept is powerful because (i) it does not require the synthesis of long coded oligomers and (ii) it gives access to an increased security level for anti-counterfeiting applications, as compared to standard 1D barcodes. However, this interesting 2D approach was not tested so far for materials labeling, mostly because the reported proof-of-concept was obtained with sequence-coded oligo(alkoxyamine amide)s, [28] which are thermolabile and thus very fragile oligomers. [10] Hence, in the present work, the 2Dcoding strategy was applied to sequence-defined oligourethanes and tested for the labelling of different commodity plastics such as polystyrene (PS), polyvinylchloride (PVC) and polyethylene terephthalate (PET).

Results and Discussion

Two libraries of sequence-coded oligourethanes were used as 2D barcodes for the labelling of thin plastics films (Figure 1). Each oligomer was prepared via a solid-phase orthogonal iterative process, as described in an earlier publication. [25] In this approach, two monomers are used as a binary alphabet. In the formed oligourethanes, a 1-bit corresponds to a methylated synthon, whereas a 0-bit corresponds to its non-methylated analogue (Figure 1).Combining this 14 Da mass difference with the specific dissociation pattern of sequencedefined oligourethanes in the negative ion mode, the nature of these synthons as well as their respective locations along these oligomeric chains can readily be distinguished in MS/MS sequencing, as reported earlier. [25] In the present work, the model digital sequence 01000001 01000010 was encoded in oligurethanes. In order to show the versatility of the 2D barcoding concept, this sequence was written in two different libraries, containing either three of four components, as shown in Figure 1.The four components library contained the oligomers P1-P4 (Figure 1 and Table 1), with sequence α-01, α-000, α-0010 and α-1000010, to be combined from the smallest to the largest one in order to compose the model digital sequence. [28] Alternatively, but using the same writing rules, [28] this sequence was generated with oligomers from the three components library, namely P1, P5 and P6 (Figure 1 and Table 1), with sequence α-01, α-000001 and α-01000010, respectively. In all oligomers, the prefix α in these sequences denotes a short spacer obtained after resin cleavage, as defined in previous works. [25,27] Table 1.ESI-HRMS characterization of the sequence-coded oligourethanes. a All the oligomers composing these libraries were first individually analyzed by electrospray ionizationhigh resolution mass spectrometry (ESI-HRMS). In all cases, the main signal could be assigned to the desired oligomer based on elemental composition derived from accurate measurements performed in the positive ion mode (Table 1).Moreover, all samples were highly monodisperse, although some impurities were sometimes detected as shown by mass spectra recorded in the negative ion mode (Figures S1-S6). Furthermore, the monomer sequence of each digital oligourethane was confirmed by negative mode ESI-MS/MS (Figures S7-S12).

Sequence

As mentioned in the introduction, three model commodity polymers were tested as host-matrices for the 2D-coded oligourethane libraries, namely polystyrene, poly(vinyl chloride) and polyethylene terephthalate. In order to show the versatility of the 2D-barcoding concept, both the three-component and the four-component oligourethane libraries were tested for labelling these three different polymers (i.e. six different samples were studied).For optimal anti-counterfeiting and traceability applications, the amount of barcodes that is included inside a product shall be as low as possible in order to not alter the physico-chemical properties of the host materials and to not be too easily spotted by potential infringers. Still, the embedded amount of the molecular barcode shall be sufficiently high to allow its extraction from the host matrix and its analytical detection. Following these prerequisites, the global oligourethane content was kept below 1 wt% in all studied films (i.e. 0.6-0.7 w% for the three-component library and 0.7-0.9 w% for the four component library).The sequencecoded libraries were incorporated in the polymer matrices using a simple solvent casting procedure. Typically, the host matrix and the oligourethanes constituting the chosen coded library were dissolved in a common good solvent(i.e. tetrahydrofuran for PS and PVC, phenol for PET) and mixed together. In most cases, it was necessary to heat the solutions in order to obtain a good dissolution of the oligourethanes, and eventually of the host matrix (e.g. for PET in phenol). The solutionswere then poured into glass plates and the solventswereevaporated to form homogeneous thin films.For PS and PVC, evaporation was performed at room temperature, whereas heating was necessary to form the PET films. It shall be noted that these mild thermal treatments (i.e. in the 50-90°C range) are not leading to noticeable oligourethane degradation. Indeed, as demonstrated in an earlier work, [25] the thermal decomposition of sequence-coded oligourethanes starts to occur around 150°C. In all cases, dried polymer films -transparent for PS and PVC and white for PET -were easily obtained via this solvent casting procedure. The formed films were then analyzed in order to verify the presence and the readability of the 2D barcodes within the host matrices. First, it was important to verify that the oligourethanes were homogeneously dispersed within the films. Taking into account the Coded library 2:

(0) (α) (1) (0) (α) (1) (0) (0) (0) (0) (0) (α) (1) (1)
very low weight fraction and the very short chain-length of the oligourethanes that wereembeddedin the host matrices, phase segregation shall not occur within the films.

Nevertheless, homogeneity was verified using a simple qualitative 1 H NMR analysis. For this, small pieces of the films were cut at different locations and placed in a deuterated solvent. It shall be noted that the chosen deuterated solvent shall not necessarily fully solubilize the host films but shall solubilize the oligourethanes. The PET and PVC membranes were studied in DMSO-d 6 , whereas the PS membrane was studied in THF-d 6 .For all types of host matrices, the NMR spectra indicated the presence of the oligourethanes in all studied regions of the films, thus ruling out the possibility of an inhomogeneous distribution(Figures S13-S15).More importantly, the possibility to extract the 2D-oligourethane mixtures from the films was studied. To do so, all films were incubated in a methanolic solution of ammonium acetate (3 mM) for 10 min. Methanol is a non-solvent for PS, PVC and PET but it partially solubilizes oligourethanes at room temperature. Composition of this incubation solution allows optimal ionization of oligourethanes in the negative ion mode and was hence specifically selected to enable direct injection of extracts in the ESI source with no prior dilution. Figure 2 shows the MS spectrum that was obtained for the PET membrane labeled with the three-component oligourethanes library. It appears clearly that the three oligomers P1, P5 and P6 of the library can be clearly identified as deprotonated species [M-H] -at m/z values of 375.2, 835.5 and 1079.6, respectively. Due to their intentional mass differences, these three oligomers lead to a detectable first-dimension fingerprint in MS. [28] Furthermore, each peak of this mass distribution can be individually sequenced by MS/MS. Figure 2a-c shows the MS/MS spectra that were obtained for the three oligomers. Their digital sequence can be very easily deciphered, showing that neither solvent casting nor matrix-storage affected their structural integrity. In the present case, the model protection sequence 01000001 01000010 was unequivocally found, as evidenced in Figure 2.Comparable results were obtained for all types of films and all coded libraries (Figures S16-S20). It shall be mentioned that MS spectra obtained for PVC films showed oligourethanes signals of with much lower intensity compared to the other two membranes and exhibited noticeable impurities, which were initially present in the commercial PVC samples as supported by MS data recorded for the extract of a non-labeled PVC film (Figure S21).

Conclusions

In summary, uniform sequence-coded oligourethaneswere used as anti-counterfeit tags in commodity polymer matrices, i.e. PS, PVS and PET. However, instead of using one type of coded oligomer as reported before, [27] multicomponent libraries containing oligomers with different chain-length, sequence and mass were used. These coded libraries could be dispersed in all types of polymer matrices using facile solvent casting procedures.

Furthermore, the oligomer mixtures could be extracted in all cases from the films and analyzed by mass spectrometry. In particular, a combination of MS and MS/MS analysis allows to reconstruct a complex digital sequence included in the whole library. These results are important because they indicate that oligourethane barcoding is not restricted to singlecoded-chains but can be extended to more complex oligomer mixtures. Furthermore, the results of this communication evidence that the concept is relatively universal and can be applied to very different types of host plastics. Indeed, the selective solvent-extraction protocols described in this work are practicaland allow direct MS analysis of the sequencecoded barcodes (i.e. without requiring barcode amplification or other pre-analysis treatments).These results open up interesting new horizons for the identification and traceability of plastic-based or plastic-containing products.

Figure 1 .

 1 Figure 1. Sequence-coded oligourethane libraries used in this work for the labelling of plastic

a

  All accurate mass measurements were performed in positive mode ESI-MS, for either b [M+H] + or c [M+NH 4 ] + .

Figure 2 .

 2 Figure 2.Mass spectrometry analysis of a PET membrane labeled with the three-

Acknowledgements J.F.L. thanks the H2020 program of the European Union (project Euro-Sequences, H2020-MSCA-ITN-2014, grant agreement n°642083) and the CNRS for financial support. The PhD position of D. K. is supported by the ITN Euro-Sequences. L.C. acknowledges support from Spectropole, the Analytical Facility of Aix-Marseille University, by allowing a special access to the instruments purchased with European Funding (FEDER OBJ2142-3341).The SEC results shown in the publication were obtained with the help of the polymer characterization service of the Institut Charles Sadron.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author marc201700426-sup-0001-S1.pdf