Saint- Clair Chabert-Liddell

Nicolas Bez

Pierre Gloaguen

Sophie Donnet

Stéphanie Mahévas

Auto-encoding GPS data to reveal individual and collective behaviour

Keywords: Bhattacharyya, Collective behaviour, Conditional variational autoencoder, Latent representation of trajectories, Network for trajectory data, Fishing vessels clustering

We propose an innovative and generic methodology to analyse individual and collective behaviour through individual trajectory data. The work is motivated by the analysis of GPS trajectories of fishing vessels collected from regulatory tracking data in the context of marine biodiversity conservation and ecosystem-based fisheries management. We build a low-dimensional latent representation of trajectories using convolutional neural networks as non-linear mapping. This is done by training a conditional variational autoencoder taking into account covariates. The posterior distributions of the latent representations can be linked to the characteristics of the actual trajectories. The latent distributions of the trajectories are compared with the Bhattacharyya coefficient, which is well-suited for comparing distributions. Using this coefficient, we analyse the variation of the individual behaviour of each vessel during time. For collective behaviour analysis, we build proximity graphs and use an extension of the stochastic block model for multiple networks. This model results in a clustering of the individuals based on their set of trajectories. The application to French fishing vessels enables us to obtain groups of vessels whose individual and collective behaviours exhibit spatio-temporal patterns over the period 2014-2018.

Introduction

The conservation and sustainable use of marine resources is a world priority that requires an accurate assessment of the status of exploited marine populations. Two sources of information are key to estimating the status of exploited marine populations around the world, namely scientific surveys and commercial fisheries' statistics. Both sources provide a partial picture of the reality. For instance, scientific surveys, when they exist, only occur at some specific time of the year and miss seasonal issues. Commercial fisheries statistics are provided all year long by the fishermen, but their use in the stock estimation process often assumes that fishing data are independent. However, fishermen collaborate or compete to increase their yields, so that the above independence assumption does not hold, with possible, but not well anticipated, detrimental effects on the assessment of the status of the exploited marine populations [START_REF] Poos | An "experiment" on effort allocation of fishing vessels: the role of interference competition and area specialization[END_REF][START_REF] Fulton | Human behaviour: the key source of uncertainty in fisheries management[END_REF]. This assumption has recently been investigated. [START_REF] Joo | Identifying partners at sea from joint movement metrics of pelagic pair trawlers[END_REF] established interaction networks between fishing vessels, studying pairs of their trajectories (called dyads) with proximity and coordination indicators (see Joo et al., 2018, for a review of indicators). Using the graphlet correlation distance, Roux et al. (2023) tested whether their corresponding interaction networks were compatible with the assumption that the trajectories were independent or not. However, the authors acknowledge that a weakness of these approaches lies in the subjective and limited choices of indicators used to describe dyads and measure pairwise distance.

The general motivation of the present work is therefore to develop a method for analysing trajectories that does not require defining indicators for dyads prior to the analysis, but that statistically reveals relevant metrics for comparing trajectories. Trajectories are characterized by a two-dimensional time series (geographical positions at each point in time). We propose to learn from the observed trajectories a low-dimensional space in which some sub-regions can be related to particular behaviours. Once projected in this low-dimensional space, trajectories can then be compared using any relevant metric tools (distance, means, clusters. . . ). These comparisons may subsequently allow us to derive groups of vessels sharing similar behaviours.

To learn a low-dimensional representation of the trajectories [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF], we propose to use a probabilistic hierarchical model. We assume that the trajectories follow a probabilistic distribution whose parameters arise from a low dimensional latent space through a non-linear transformation (mapping), called decoder. The prior distributions in the latent space are assumed to be known, and the resulting posterior distributions can then be recovered by learning an inverse mapping. In practice, we use a variational approach [START_REF] Blei | Variational inference: A review for statisticians[END_REF] to approximate the posterior distribution. The variational distribution is restricted to the Gaussian family, whose parameters are assumed to be non-linear functions of the observations. The mapping from observations to the parameters of the variational distribution is called the encoder. Both the decoder and the encoder are parameterized by neural networks and their parameters are learnt jointly by maximizing a proxy of the log-likelihood.

This modelling framework is known as a variational auto-encoder (VAE; [START_REF] Kingma | Auto-encoding variational bayes[END_REF][START_REF] Rezende | Stochastic backpropagation and approximate inference in deep generative models[END_REF] and can be seen as non-linear extensions of probabilistic PCA [START_REF] Bishop | Bayesian pca[END_REF][START_REF] Lawrence | Probabilistic non-linear principal component analysis with gaussian process latent variable models[END_REF]. Compared to a classical auto-encoder [START_REF] Kramer | Nonlinear principal component analysis using autoassociative neural networks[END_REF], the agreement between the variational distribution and the prior on the latent space acts as a regularization term, which prevents over-fitting and brings more structured latent representations. Adding covariates to the model is straightforward and leads to a conditional variational auto-encoder (CVAE; [START_REF] Kingma | Auto-encoding variational bayes[END_REF], which ensures that the latent representations do not reflect spurious correlation with these known covariates.

The choice of a neural network architecture for both mappings is guided by the nature of the data (spatio-temporal) and by the question (exploration of the trajectory as a whole, i.e. large-and medium-scale properties). In this work, we specifically use Convolutional Neural Networks (CNN). The interested reader may refer to recent reviews on deep generative models and neural network architectures for spatio-temporal data [START_REF] Gao | Generative Adversarial Networks for Spatio-temporal Data: A Survey[END_REF][START_REF] Graser | Deep learning from trajectory data: a review of deep neural networks and the trajectory data representations to train them[END_REF].

Quantifying similarity between latent representations of trajectories is usually done through the computation of indices using point estimates (e.g. Euclidean distance or cosine similarity on the mean of the variational distributions). We chose here to exploit the full benefits of the probabilistic trajectory modelling, and we quantify the similarity between latent representations with the Bhattacharyya Coefficient (BC; [START_REF] Bhattacharyya | On a measure of divergence between two statistical populations defined by their probability distribution[END_REF], a measure of overlap between probability distributions. Using the full distribution allows capturing the uncertainty associated with the latent representation of each data trajectory. The overlap coefficient can be interpreted as a quantification of similarity between daily individual spatio-temporal behaviours.

Comparison of trajectories using similarity between their representations in the low-dimensional latent space is useful in several ways. First, comparing the trajectories of a single fishing vessel for different fishing days (adjacent or not) helps to quantify the variation of its daily individual behaviour. Alternatively, for a given fishing day, if we compare the trajectories of a group of vessels we can derive a proximity network, i.e. a binary graph where a node is a vessel and an edge between two vessels denotes a similarity above a pre-defined threshold. Using an extension of the Stochastic Block Model (SBM; [START_REF] Holland | Stochastic blockmodels: First steps[END_REF][START_REF] Lee | A review of stochastic block models and extensions for graph clustering[END_REF], for a review) for collections of networks (Chabert-Liddell et al., In press), these daily networks can be analysed to derive groups of vessels with similar behaviours, to understand how the different groups relate to each other and how they change over time periods.

The general approach outlined in this article and specifically, the concept of encoding individual trajectories to later achieve the clustering of individuals, could be of interest in analysing trajectories of different types of objects (or individuals). The method is therefore presented in a general context, and details specific to our case study are set and discussed in the applications and results.

The manuscript is organised as follows. Section 2 describes the trajectory data, introduces the conditional variational auto-encoder (CVAE) model used for analysis, and presents the tools employed for analysing the latent representations. Then the results are divided into several subsections. Section 3.1 covers the selection of the latent space for representation, and Section 3.2 provides intuitive insights about these representations. Section 3.3 deals with the analysis of individual behaviours, while Section 3.4 focuses on collective behaviours. In this latter section, we present the construction of proximity graphs from the latent representations of fishing tra-jectories and the application of network analysis tools for grouping fishing vessels. A discussion on the application of our method to fishery science case study and some perspectives conclude the paper in Section 4.

Material and methods

Data set

The data consists of M = 20512 daily GPS trajectories from B = 33 different fishing vessels over a 5-year period. The number of trajectories per vessel ranges from 140 to 968 (median 701). All the trajectories start and end at the same point, a fishing harbour in French Brittany. The positions of the vessels are obtained from the mandatory Vessel Monitoring System (VMS, [START_REF] Hinz | Confidentiality over fishing effort data threatens science and management progress[END_REF], recording at a roughly hourly rate. To synchronize all the trajectories (i.e., have exactly the same timestamp), we use a linear interpolation of the position recorded by the VMS as recommended by Joo et al. (2018). Hence, each trajectory is a matrix of dimension d Y = H × D. In this study, the number of hours is H = 24 and we consider 2d trajectories (geographical positions expressed in longitude and latitude), so D = 2.

We denote

y m ∈ R d Y the m th trajectory (1 ≤ m ≤ M = 20512
), and y = {y m } m=1,...,M the complete data set. In the following, when necessary, we will use the notation y (b,t) to specifically refer to the trajectory of the vessel b at day t, if it exists, where 1 ≤ b ≤ B, and t spans over the whole period of the data set. In a general setting, these trajectory data can be enriched with additional covariates. For each trajectory 1 ≤ m ≤ M , we denote x m ∈ R dx the set of available covariates. In this paper, the available covariate is the day of the year at which the trajectory occurs. For j ∈ 1; 365 , the time of the year is given by the vector cos 2π×j 365 , sin 2π×j 365 (and consequently d x = 2) to take into account the periodicity.

Methods for learning the latent representations 2.2.1 Probabilistic generative model

We consider a probabilistic framework where the observed trajectories are realizations of random vectors Y = {Y m } m=1,...,M following a probabilistic distribution characterized by a highly non-linear transformation of the realization of low dimension latent random vectors. Formally, in this paper, we suppose that:

Z m ind. ∼ p 0 := N (0, I d Z ) , 1 ≤ m ≤ M Y m |Z m = z m ; x m ind. ∼ N (g(z m , x m ; θ); σ 2 I d Y )
(1) where • N (µ, Σ) denotes the normal distribution with mean µ and matrix of variance Σ;

• I d denotes the identity matrix of dimension d;

• g(•, •; θ) is a function from R d Z × R dx to R d Y , depending on some unknown parameter θ ∈ R d θ ;
• σ 2 > 0 is a variance parameter whose role is discussed in Section 2.2.2.

The key part of this generative model is the function g(•, •; θ), which maps the latent space to the space of trajectories. This function, often called decoder, must be chosen within a class of non-linear mappings, parameterised by θ, typically neural networks. The precise choice for this function is discussed in Section 2.2.3. The main learning objective is then to estimate θ from observed y.

Variational inference

For the probabilistic latent model given by Equation ( 1), the log-likelihood is given by integrating over the whole latent space, namely: log p(y; x, θ) = log p(y|z; x, θ)p 0 (z)dz.

The integral on the right-hand side has no explicit expression, and is usually hard to approximate by numerical approximations or Monte Carlo methods, due to the dimension of z. An alternative inference method is the variational approach [START_REF] Blei | Variational inference: A review for statisticians[END_REF]. This approach exploits a surrogate of the likelihood function, the evidence lower bound (ELBO). For any probability distribution function q(z|y; x, ϕ), parameterised by some parameter ϕ ∈ R d ϕ : log p(y; x, θ) ≥ log p(y; x, θ) -D KL (q(•|y; x, ϕ)∥p(•|y; x, θ))

(2) = E q(•|y;x,ϕ) [log p(y|Z; x, θ)] -D KL q(•|y; x, ϕ)∥p 0 (•)

=:ELBO(θ,ϕ) , (3)
where D KL stands for the Kullback-Leibler divergence. The inequality becomes an equality in the single case where q(•|y; x, ϕ) = p(•|y; x, θ), and maximising the ELBO with respect to ϕ (for a fixed θ) is equivalent to minimising the Kullback-Leibler divergence between q(•|y; x, ϕ) and p(•|y; x, θ) as seen in Equation (2). Therefore, by maximising the ELBO jointly with respect to both θ and ϕ, for a given family of variational distribution q(•|y; x, ϕ), one obtains θ that maximises a proxy of the log-likelihood and ϕ that gives the best member of the family chosen to approximate the distribution p(•|y; x, θ).

A flexible and standard framework for the variational family is to assume normality, i.e., that q(z m |y m ; x m , ϕ) is the p.d.f. of a normal distribution in dimension d Z , with mean µ m ∈ R d Z such that µ m = µ(y m , x m ; ϕ) and a diagonal covariance matrix whose diagonal is noted γ 2 m = γ 2 (y m , x m ; ϕ). Both means and variances, therefore, depend on observed data through a mapping parameterised by ϕ. This mapping h(•,

•; ϕ) : (y m , x m ) → (µ m , log γ 2 m
) is usually called the encoder. The loss function to minimise is -ELBO, taking its form in Equation (3), comprises two terms whose expressions in the Gaussian cases come easily: a recon-struction error term,

-E q(•|y;x,ϕ) [log p(y|Z; x, θ)] = 1 2σ 2 M m=1 E q(•;ϕ) ∥y m -g(Z m , x m ; θ)∥ 2 2 + M d y 2 log σ 2 +O(1),
(4) and a regularization term which takes advantage of the fact that the dimensions in the latent space are independent and that the divergence is the sum of the divergence by dimension,

D KL (q(•|y; x, ϕ)∥p 0 (•)) = 1 2 M m=1 d Z k=1 γ 2 k + µ 2 k -log γ 2 k -1 . ( 5
)
The minimization of the negative ELBO is then done through a stochastic gradient descent algorithm with mini-batches using the ADAM optimizer [START_REF] Kingma | A method for stochastic optimization[END_REF]. During the optimization, for gradient back-propagation purposes we use the reparameterization trick [START_REF] Kingma | Auto-encoding variational bayes[END_REF][START_REF] Rezende | Stochastic backpropagation and approximate inference in deep generative models[END_REF], each entry of the encoder z m is simulated independently as z m = γ m ϵ + µ m where ϵ is the realization of a N (0, I d Z ) (cf. Figure 1).

In accordance with the deep learning literature, the model with covariates is called a conditional variational auto-encoder (CVAE) and the one without covariates, where h and g are not functions of x, a variational auto-encoder (VAE).

Remark about the variance parameter σ 2 In the context of variational inference, the variance parameter σ 2 is difficult to estimate and is usually fixed [START_REF] Skafte | Reliable training and estimation of variance networks[END_REF]. If fixed, for a Gaussian distribution, choosing the variance corresponds to choosing the importance of the regularization term compared to the reconstruction term [START_REF] Higgins | beta-VAE: Learning basic visual concepts with a constrained variational framework[END_REF]. Higher variance leads to higher regularization, as can be seen in Equation (4). When σ 2 is fixed to 1, the reconstruction term of the negative ELBO corresponds to the L2 loss.

Neural network architecture

The practical implementation of the framework described above needs a parametric family of non-linear mapping for the encoder and the decoder. This is done by choosing appropriate neural network architectures.

It is classic to view a trajectory as an image on a two-dimensional grid at the cost of losing the time stamps in the trajectories [START_REF] Olive | Deep trajectory clustering with autoencoders[END_REF][START_REF] Liang | An unsupervised learning method with convolutional auto-encoder for vessel trajectory similarity computation[END_REF]. In the present study, the temporal dimension of the trajectory is important for assessing the synchronization or time lag of collective behaviour. So as not to lose it, we use the time stamps of the trajectory as pixels and the spatial position as channels.

A 2-dimensional trajectory of length H can be viewed as a one-dimensional image of H pixels with 2 channels.Using this representation, a common architecture used in deep generative models is a convolutional neural network (CNN). For generating trajectories, we build from the architecture proposed by [START_REF] Radford | Unsupervised representation learning with deep convolutional generative adversarial networks[END_REF] and successfully used by [START_REF] Roy | Using generative adversarial networks (gan) to simulate central-place foraging trajectories[END_REF] to generate bird foraging trajectories.

Figure 1: Model and its neural network architectures. Each data is transformed into an encoder input, the covariates are duplicated to match the length of the 2-dimensional trajectory. For each box, the vertical axis is ordered (image pixels) while the horizontal axis is unordered (channels). The input of each neural network is in blue, the covariates in orange, the 3 hidden layers in green and the output in red.

This architecture was adapted to a lower dimensional problem and to handle data with covariates. The two covariates x are added to the trajectory image on two additional channels and replicated to match the length of the trajectory. The architecture of the encoder is then chosen to be as symmetric as possible with respect to the decoder. A synthetic view of the model and its architecture is shown in Figure 1. More details about the architecture and pseudocode for PyTorch implementation [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] are provided in Appendix A.

Methods for the analysis of the latent representations

Quantifying similarities in the latent space

Each trajectory y m , possibly knowing its date x m , is mapped through the encoder h(•, •; φ) to a Gaussian distribution with mean μm and variance γ2 m I d Z (therefore assuming independence). In this latent space, the similarity between two encoded trajectories is quantified by a measure of overlap between distributions, the Bhattacharyya coefficient (BC; [START_REF] Bhattacharyya | On a measure of divergence between two statistical populations defined by their probability distribution[END_REF], also referred to as the Hellinger affinity), which integrates the entire variational distributions and has an expression in closed-form for independent multivariate Gaussian distributions. The BC between two fishing vessels b and b ′ fishing on day t and t ′ , is given by:

BC t,t ′ b,b ′ = z q(z|y (b,t) ; x (b ′ ,t ′ ) , φ) • q(z|y (b ′ ,t ′ ) ; x (b ′ ,t ′ ) , φ)dz. ( 6
)
BC tt ′ bb ′ ∈ (0, 1] and is equal to one if and only if the two distributions are the same. The closed-form expression of the BC between two multivariate Gaussian distributions with diagonal covariance matrix is the following:

BC t,t ′ b,b ′ = d Z k=1 2γ (b,t),k γ(b ′ ,t ′ ),k γ2 (b,t),k + γ2 (b ′ ,t ′ ),k exp - 1 4 (μ (b,t),k -μ(b,t),k ) 2 γ2 (b,t),k + γ2 (b ′ ,t ′ ),k
.

Analysis of the similarity between trajectories is then performed at the individual level by looking at the overlap coefficient of a vessel b at two different times t and t ′ , noted BC t,t ′ b := BC t,t ′ b,b , or at the collective level by looking at the overlap of two vessels b and b ′ on the same day t, noted BC t b,b ′ := BC t,t b,b ′ . The overlap coefficient is further binarized to look for high overlapping trajectory representations, by computing BC t,t ′ b,b ′ ≥ s where s ∈ (0, 1] is a users' defined threshold (set to 0.8 in the present study). Finally, the behaviour stability index of a vessel over a set of couples of days S is computed as the mean of the binarized coefficients:

SI S b = 1 |S| (t,t ′ )∈S 1 BC t,t ′ b ≥s . (7)
Similarly, within the fleet, we define a general stability index computed over all possible trajectories irrespective of the vessel to serve as a baseline for comparing the stability of individual behaviours: b , where S d = {(t, t ′ ) : t ′ = t+1}, gives information on the proportion of similar trajectories from one day to the next (daily stability), while DSI = SI S d measures the daily stability within the fleet.

SI S = 2 B(B -1) B b,b ′ =1,b̸ =b ′ 1 |S| (t,t ′ )∈S 1 BC t,t ′ b,b ′ ≥s , (8)

Tools for analysing proximity graphs

A collection of proximity graphs is built to analyse the collective behaviour. During a time period, two vessels are said to be in proximity if the distributions of their trajectories in the latent space for the same day are highly overlapped, i.e. if the qquantile of their BC distribution for the time period is greater than a given threshold s. That is, for each time period S ⊂ 1; T , the adjacency matrix of the proximity graph is defined as:

A S bb ′ =    1 if 1 |S| t∈S 1 BC t b,b ′ ≥s ≥ q, 0 otherwise. (9)
Only vessels with observed trajectories during this time period are represented in the graph. Self-interactions and interaction between vessels with no common fishing days are encoded as NA.

An extension of the Stochastic Block Model (SBM, [START_REF] Holland | Stochastic blockmodels: First steps[END_REF]) is then used to analyse the proximity graphs, called networks hereafter. SBM allows grouping fishing vessels exhibiting similar connectivity patterns with the rest of the fleet. Fishing vessels from the same group are said to be equivalent. Hence, each network is described in a condensed way by its mesoscale structure (the proportion of each group and the probability of connection within and between each group). Extension of the SBM for collection of networks (colSBM, Chabert-Liddell et al., In press), potentially enables finding a common mesoscale structure describing all the networks of the collection, and to determine if this structure is sufficient to describe the whole collection (if not, the collection is divided into sub-collections of networks with similar mesoscale structure). This model does not take into account the vessel ID. Mathematical details of the model are provided in Appendix B.

A workflow of the data analysis, explaining how the different methods introduced in this section interact, is given in Figure 2. The aim is to choose the model with the smallest-dimensional latent space that is sufficient to describe the trajectory data. Figure 3 shows the loss (-ELBO) and the regularization term (D KL ) for the VAE and the CVAE with d Z = 2, . . . , 5. From d Z = 3 for the CVAE and d Z = 4 for the VAE, the loss does not decrease much with additional dimensions, and the D KL term stops increasing meaning that the structure of the latent space does not depart from the prior distribution anymore, and hence adding dimensions to the latent space is not necessary to describe the trajectories. From now on, we consider d Z = 3 for the CVAE model and d Z = 4 for the VAE model.

Ordering the dimensions of the latent space

As the dimensions of the variational distribution are independent (diagonal covariance matrix), the BC between two variational distributions is the product of the BC between each dimension of the two distributions. We can then look at the distribution of the BCs for each dimension and order the dimensions of the latent by level of importance (Figure 4,left). If BCs are close to one, it means that the distributions have a very high overlap on this dimension and that, this dimension has a low impact on the total metric and does not bring much information on differentiating the trajectories in the latent space.

The dimensions where BC are close to 1 correspond also to the dimensions where the variational distributions are close to the prior (D KL massively close to 0; Figure 4, right). These dimensions are thus only useful to describe a few unusual trajectories, which correspond to the one in the tail of the distribution for both BC and D KL . It happens that while the analysis of the outputs of CVAE (resp. VAE) relies on three (resp. four) dimensions, the first two (resp. three) dimensions are sufficient to represent the trajectories in the latent space without significant loss of information.

Linking latent representations to trajectories

With the objective to explore how representations in the latent space are related to trajectories in the geographical space and to choose between VAE and CVAE frameworks, 3 vessels with very contrasting individual behaviours have been selected (vessels number 2, 6 and 19).

Choosing between CVAE and VAE

For each of the three selected vessels, we investigate the time series of their variational distributions on each dimension of the latent space for both the VAE and the CVAE models (Figure 5), their trajectories and the corresponding plots of the mean values in the first three dimensions of the latent space (Figure 6).

For the VAE model, µ (b,•),3 exhibits a very clear sinusoidal seasonality (Figure 5), taking negative values in winter (around -2) and positive values in summer (around 2). This effect is taken into account by the covariates of the CVAE for which times series do not show any seasonal variations anymore. In this regard, a CVAE model with a three-dimensional latent space and time as covariate provides equivalent outputs than a VAE approach with four dimensions. Given that the third dimension of the CVAE is practically useless in terms of graphical representations, this is a first indication for promoting CVAE against VAE.

Both the trajectories and the mode of their latent representations are progressively less compact for vessel 2, 6 and 19 (Figure 5 and grey points in Figure 6) indicating that trajectories with increasing length (i.e. trajectories whose point furthest from the port is getting farther and farther away) are mapped by expanding points in the latent space for both CVAE and VAE. Vessel 2 exhibits only two kinds of behaviours in its trajectories in December, one sailing west and the other one sailing southeast. The fishing behaviours in June are more diverse but one trajectory going west overlaps the trajectory that occurred in December while being slightly longer. These two particular trajectories are overlapped in the CVAE latent space near the origin, while they correspond to two different points in VAE. A similar pattern emerges for vessel 19, the points corresponding to the long trajectories to the south-east direction in June and December overlap for the CVAE (near (-1.5, 0)) but not for the VAE (near (2, 0)).

If the same intrinsic fishing behaviour produces longer trajectories in summer than in winter due to a change in the distribution of the target species, then a CVAE approach that translates these two different trajectories into a single point in latent space provides a better mapping of behaviours than a VAE approach does. Hereafter, we will thus only consider the latent space of the CVAE model.

Deeper analysis of the latent space for CVAE

We show how the dimensions of the latent space are related to the length of the trajectories and to an east-west and north-south gradient. For each trajectory data, we look at the parameters µ and γ 2 of its variational distribution from the CVAE model. In Figure 7 each trajectory of the dataset (before normalization) is coloured according to µ and γ 2 for each dimension. Interpretation for each parameter follows:

1. µ 1 varies on a longitudinal gradient with µ 1 ranging from smallest (negative) values (µ 1 ≈ -2) to largest positive values (µ 1 ≈ 1) from east to west, with µ 1 = 0 corresponding to trajectories going south-west from the harbour.

2. µ 2 varies along a latitudinal gradient, with trajectories from north to west corresponding to the smallest values (µ 2 ≈ -2), while trajectories far to the south have the largest values (µ 2 > 2).

3. µ 3 does not exhibit a clear spatial pattern, except from some trajectories going west with the largest values (µ 3 > 2). 4. γ 2 1 is very low for each trajectory, indicating good confidence in the ability to characterize the direction of the trajectory using the first dimension of the latent space.

γ 2

2 is the smallest for the trajectory going south and the largest for the easternmost trajectories. The second dimension of the latent space will have difficulty distinguishing the latitudinal position of trajectories in the eastern part of the region. It corresponds to an area with a low density of observations. 6. γ 2 3 is the largest (close to one) for short trajectories going south of the harbour.

Those trajectories have also µ 3 ≈ 0 meaning that the variational distribution did not depart from the prior for these trajectories, and this dimension is not informative. The longest trajectories to the west have much lower variance in their variational distribution. Their corresponding means are µ 1 ≈ 1 and µ 3 > 2.

To further quantify the relationship between the latent representations and the actual trajectories, we evaluate the sensitivity of two key trajectory characteristics, i.e. the maximal distance to the harbour (denoted D max ) and the main direction of the trajectory (denoted theta max , the angle corresponding to the maximally distant point in the trajectory), with respect to the two first dimensions (z 1 and z 2 ) of the latent space and to the month of the year (coded as cos(m), sin(m)). To this end, a grid with z 1 spanning from -2 to 1.2, z 2 from -2 to 2 and m corresponding to each month of the year is built (z 3 is fixed to 0). Trajectories are simulated by using the decoder g(•, •; θ) obtained after training the CVAE model. From these trajectories, the response variables D max and θ max are computed. Then, these two quantities are Figure 7: Trajectories coloured by the mode of the parameters of their variational distribution in the latent space (µ above and γ 2 below) on each of the 3 dimensions of the CVAE. The red lozenge represents the harbour (geographical position of vessel departure and arrival), which is located at (-4.28, 47.78).

Figure 8: Sensitivity of trajectory characteristics (D max : distance to the port; θ max : main direction) to the conditioning variable and the first 2 dimensions. predicted using a random forest (with the ranger R package; [START_REF] Wright | ranger: A fast implementation of random forests for high dimensional data in C++ and R[END_REF] having (z 1 , z 2 , cos(m), sin(m)) as input. The estimated sensitivity index is computed using the Sobol MDA index, designed to be a proxy of the total Sobol index for random forests accounting for dependency between variables [START_REF] Bénard | Mean decrease accuracy for random forests: inconsistency, and a practical solution via the sobol-mda[END_REF]. This index ranges from 0 to 1, equalling 0 if there is independence between the predictor and the response, and 1 if the response is a deterministic function of the predictor. Figure 8 displays these indices for the four explanatory variables. As expected, one can see that the maximal distance and main directions are strongly sensitive to the first dimension, and, to a lesser extent, to the second dimension. We also can see that the distance to the port is not sensitive to the time of the year, but the direction is affected by this covariate.

In summary, the latent space allows to reliably discriminate trajectories along a positive east-west gradient on the first dimension, northward trajectories with low values on the second dimension and long trajectories with high values on the third dimension. On the other hand, short and/or southward trajectories are more difficult to identify in latent space.

Exploring individual behaviours

Using the individual stability index and the global stability index at the fleet scale, we investigate the behaviour of the vessels. For all fishing vessels DSI b > GSI b (Figure 9) indicating that fishermen tend to reproduce similar trajectories from one day to the next and don't change drastically fishing grounds every day. The two indices are correlated, but heterogeneity can be found in the variety of individual behaviours between vessels. For instance, vessels 2, 3 and 5 have high global stability, but vessel 2 has by far the most stable trajectories between following days (DSI 2 > 0.65). Vessels have very different daily stabilities, ranging from 0.27 to 0.67, with a fleet average of about 0.45. In comparison, the global stability of each vessel is very low (ranging from 0.025 to 0.175) indicating that vessels have many kinds of individual behaviour, i.e. they globally often change of trajectories.

As expected, amongst the three selected vessels mentioned in section 3.2, the vessel with the most variable time series (vessel 19 in Figure 5) is also the one with the smallest temporal stabilities (DSI 19 and GSI 19 ). It explores widely the fishing area and its trajectories change frequently. Conversely, the one with the most stable time series and compact latent space (vessel 2 in Figure 6) is the one with the largest global and daily stability indices (see Figure 9). The behaviour of vessel 6 lies in between and is close to the average individual stability (green triangle in Figure 9).

At the fleet level, the vessels have a very low daily stability (about 0.16, black diamond at the bottom left of Figure 9). Each of the vessels has on average more day-to-day stable trajectories than the fleet (DSI b > DSI). This is not always true for distant time trajectories, as about a third of the fleet has GSI b < GSI. This indicates that vessels have many different behaviours that are explored below using network analyses.

Exploring collective behaviour

Choosing adequate time spans to analyse collective behaviour is of primary importance. Relying on the CVAE outputs, a daily time step is likely to result in highly unstable behaviour. Conversely, considering the entire period (5 years) masks seasonal behaviour or yearly trends. The quarter is a common time period used to define fishing behaviour (ICES, 2022), and therefore we build 20 (one per quarter over 5 years) proximity graphs from the latent representations of the vessels' trajectories using the Bhattacharyya coefficients (see section 2.3.2). The adjacency matrices of the proximity graphs are shown in Appendix C, Figure 12.

Clustering the vessels into groups

We show that the structure of the proximity graphs is stable over time and that the fleet contains two main groups of vessels exhibiting exclusive spatio-temporal behaviours.

To analyse the networks, we use the colSBM model introduced in Section 2.3.2. The model selection criterion does not favour the partition of the collection of 20 proximity networks into sub-collections. It means that the structures of the networks are similar between quarters and that the common mesoscale structure shown in Figure 10 is sufficient to describe the 20 quarterly proximity graph. The fleet is divided into 5 groups of fishing vessels labelled A to E (Figure 10, right) and ordered in the matrix summary of the mesoscale structure from top to bottom, left to right (Fig. 10, left). Groups B and D are two exclusive communities of highly internally connected fishing vessels, representing about a third of the fleet each. A and E are smaller groups that are on the periphery of respectively groups B and D, with group E regrouping vessels nearly without mutual connection. It is worth noticing that there are almost no connections between groups A, B and groups D, E.Finally, C is a small group of fishing vessels with heterogeneous behaviours and connecting with both sides of the fleet. The memberships of fishing vessels to group D are very stable across quarters, with however ephemeral switch to the peripheral group E (Figure 10). On the contrary, vessels from the rest of the fleet switch more often between groups C, B and A.

The analysis of the fishing behaviour in terms of daily and global stability can be further highlighted by the group memberships (Figure 9). Groups B and D representing two-thirds of the fleet have on average the same daily stabilities of their trajectories (DSI B = 0.46 and DSI D = 0.47), with strong variabilities around. However, fishing vessels from group B have nearly systematically more diverse trajectories, and lower global stability, than the ones from group D (GSI B = 0.06 < GSI D = 0.11). Fishing vessels in the peripheral groups (A and E) have much lower daily stabilities (DSI A = 0.38 and DSI E = 0.37) and behave alike group B for global stability.

Analysis of the trajectories and their latent representations by groups

Yet, we can investigate the specificities of trajectories within groups. The distributions of the mean of the first two dimensions of the latent representations of the trajectories for the five groups are shown in the lower part of Figure 11a, while the trajectories are shown in the upper part. Both are very different for groups D and E from the ones of groups A, B and C. We can briefly summarize the trajectory behaviour for each group as follows:

A, B Latent representations of trajectories of these two groups look very similar.

Both have a wide range of values in the latent space, with high concentrations around the same place. There are some differences between the two groups. Group B takes more extreme values in the latent space (for instance on the top left, which corresponds to the long trajectory to the south). But, as seen in Figure 10, both groups are composed of the same fishing vessels. In particular, they make long trajectories.

C Looks very similar to group A on the trajectory space, but the concentration is more uniform in the latent space.

D Values for the first dimension are taken on a short interval from -1 to 0.5, with a high concentration around (0, 0). Notice how the distribution is almost complementary from the one of group B on the latent space, and included on the trajectory space. The fishing vessels of this group make shorter trajectories.

E Takes lower values in the second dimensions and higher in the first dimensions than the other groups. Its fishing vessels have the northernmost trajectories in the trajectory space. They are the only ones to consistently explore this fishing ground.

In the initial dataset, we had catches in values (euros) for each daily trajectory (fishing trip). We investigated whether these groups differed in terms of gross revenue. The main source of gross value for the fleet consists of fishing crustaceans (Nephrops norvegicus), but there are side significant incomes which come from fishes (mainly monkfish (Lophius), hake (merluccius merluccius), sole (solea solea), megrim (Lepidorhombus), and john dory (Zeus faber)). The gross values per fishing trip are variable between the years and also between vessels within the group. They are significantly different between groups A, B, C and groups D, E (1400€ on average for group D and 2400€ on average for groups B and A; Figure 11b). The vessels from group E have median incomes close to the one of group D, but in some quarters have much higher incomes in the upper quantiles.

As the trajectories of groups D and E are less distant from the port and therefore generate lower fuel costs, we may wonder whether these vessels are not looking for a cost-revenue trade-off. It would be interesting to explore net revenues, which unfortunately were not available for this study.

Discussion and perspectives

We proposed the utilization of deep learning methods, especially variational autoencoders, to learn a low-dimensional representation (a probability distribution) of trajectory data. This approach allowed for the transformation of complex trajectory data into a more interpretable representation facilitating the analysis of the individual and collective spatio-temporal behaviours. Temporal variability in the latent representations was attenuated by using the day of the year of a trajectory as a covariate in the variational auto-encoder. This statistical conditioning also allowed to decrease the number of dimensions of the latent space from 4 to 3. In our study, we used an overlap index, the Batthacharyya Coefficient, between the distributions of the representations to quantify the similarities in individual behaviour and build proximity graphs. These proximity graphs facilitated the clustering of fishing vessels, resulting in groups of vessels based on the stability of individual behaviours. The exploration of these groups has highlighted differences in spatial and economic patterns. It was observed that groups exhibiting shorter trajectories have lower gross values per trip, most likely revealing a cost/revenue trade-off strategy. Additionally, these groups exhibited lower overall effort (cost) per trip, suggesting a more efficient utilization of resources.

In fisheries, one of the aims of building groups of vessels with similar behaviour is to identify reference groups (reference fleets) in order to estimate indices of abundance from their catches, which reveal inter-annual variations in fish abundance. The quality of these abundance indices depends on good sampling by fishermen of the species' distribution area. One possibility is to use these groups to construct a stratified estimator that takes into account the spatial structure of exploration by vessels to remove bias in current estimators and yet improve stock assessment and management decisions [START_REF] Poos | An "experiment" on effort allocation of fishing vessels: the role of interference competition and area specialization[END_REF][START_REF] Alglave | Combining scientific survey and commercial catch data to map fish distribution[END_REF].

Whether these indices can be improved by this method or not relies on the assumption of trajectory independence [START_REF] Bez | The catchability in a spatial context: a simulation exercise[END_REF]. We could use the generative aspect of the CVAE model to generate trajectories to test a null hypothesis of trajectory independence. This aspect of the model is not present in this paper, except to relate the characteristics of the trajectories to the dimensions of the latent space.

We now discuss the choice of the neural network architectures, the indices to calculate proximity between latent representations and the hypotheses that were made when comparing behaviours and building proximity networks.

In this work, we used Convolutional Neural Networks (CNN) as the neural network architecture for both the encoder and the decoder. An alternative would have been to use Recurrent Neural Networks, in particular a Long-Short Term Memory (LSTM) architecture, as they are classic for temporal data. However, [START_REF] Roy | Using generative adversarial networks (gan) to simulate central-place foraging trajectories[END_REF] point to the inability of LSTM in a generative model to consistently generate trajectories which start and end at the same position (cycle), which is a common feature in all the trajectories of our data set. In addition, LTSM are particularly appropriate to i) capture short-scale behaviour while we are more interested in medium and large-scale characteristics of the trajectory and ii) provide an operational prediction of a vessel's next position based on its past trajectory [START_REF] Capobianco | Deep learning methods for vessel trajectory prediction based on recurrent neural networks[END_REF], for instance, to avoid collisions [START_REF] Jia | Motion primitives learning of ship-ship interaction patterns in encounter situations[END_REF].

Because of the heavy parametrization of the nonlinear mapping learnt by the VAE, the embedding usually lacks interpretability. In order to circumvent this drawback, some authors make a tradeoff between reconstruction error and interpretability by still using a non-linear encoder but a linear decoder (Svensson et al., 2020, e.g.), while others aim at disentangling the dimension of the latent space to represent features (Higgins et al., 2017, e.g.), as it is difficult to relate a dimension to a feature of the data. The trajectory of a fishing vessel being smooth and regular, starting and ending at the same point in space, we were able to obtain very low-dimensional representations. This allowed for post hoc analysis to interpret the latent representations and relate them to the trajectories. Low-dimensional space was also beneficial for graphical representation as by properly ordering the dimension of the latent space, we did not need to rely on further dimensional reduction tools (Van der Maaten and Hinton, 2008, e.g. t-sne) to visualize the latent space.

We decided to use a notion of overlap between distributions (BC) to compare trajectories in the latent space, mainly for interpretation purposes. If willing to have a metric to compare distributions, many exist in closed-form for independent multivariate Gaussian. For instance, BC is related to the Hellinger distance by the following relationship: H(p, q) = (1 -BC(p, q)) 1/2 . We also set a threshold s on the BCs to define a high overlap in the latent space. In terms of comparing individual behaviours, we are just interested in proximity, as having two different or very different individual behaviour do not change the qualitative analysis (for instance if a trajectory is going far to the west, it does not matter if the other trajectory is going to the south or to the east, both are the fruits of different behaviours than the first trajectory). Hence, a need to binarize the indices. If willing to use network analysis tools to quantify the collective behaviour and group the fishing vessels without using binary networks, would be to use a unimodal continuous distribution on (0, 1) (e.g. a Beta distribution with constrained parameters or a continuous Bernoulli distribution Loaiza-Ganem and Cunningham (2019)) as the emission distribution for the edges.

Finally, we set a quantile q when building similarity networks over a defined time period. A solution would be to study networks constructed with the same temporal scale as the data, i.e. daily proximity graphs. Acknowledging the necessity to distinguish between noise or temporary local effects and collaborative effects, we plan on extending existing approaches to reconstruct networks from unreliable observations [START_REF] Le | Estimating a network from multiple noisy realizations[END_REF][START_REF] Young | Bayesian inference of network structure from unreliable data[END_REF], by accounting for the assortative community structure bias induced by the use of proximity graphs. It will help us to better account for collaborative effects in similarity networks, which can further enhance the reliability and accuracy of the analysis. The analysis would still need to be conducted at an appropriate time period to define a fishing fleet.

Aknowledgements

This work was financially supported by France Filière Pêche (FFP) as part of the MACCOLOC project. independent and distributed as:

A (q,y) bb ′ |W (q,y) br W (q,y) b ′ r ′ = 1 ∼ B(α rr ′ ),
where B is the Bernoulli distribution. As the networks are undirected, we set A (q,y) b ′ b = A (q,y) bb ′ . The above defines the SBM model, SBM B (q,y) (R, α, ρ). Then the collection of networks is distributed as: A = (q,y)

A (q,y) , where A (q,y) ∼ SBM B (q,y) (R, α, ρ).

For R fixed, the inference of the model parameters and the latent group are done through the maximization of a variational approximation of the log-likelihood. A penalized model selection criterion is used to choose the number of groups R (Chabert-Liddell et al., In press).

C Quarterly proximity graphs

  GSI b = SI Sg b , where S g = {(t, t ′ ) : t ̸ = t ′ } indicates the proportion of similar trajectories per vessel globally over time (global stability). Respectively, GSI = SI Sg quantifies the global stability of the fleet. DSI b = SI S d

Figure 2 :

 2 Figure 2: Workflow. 1) Learning latent representations of trajectory data using (conditional) variational auto-encoders. 2) Analysis of the relation between the latent representations and the trajectory data. 3) Analysis of daily and global stability of individual behaviours using the Bhattacharyya Coefficients computed from the latent representations. 4) Analysis of the collective behaviour with the Stochastic Block Model using proximity graphs between individuals built from the Bhattacharyya Coefficients computed from the latent representations.

Figure 3 :

 3 Figure 3: Loss (-ELBO) and regularization term (DKL) for the first 1000 epochs for the VAE and the CVAE with 2 to 5 latent dimensions.

Figure 4 :

 4 Figure 4: Empirical distribution of BC between couples of variational distribution and Kullback-Leibler divergence between variational distribution and normal Gaussian for each dimension of the latent space.

Figure 5 :

 5 Figure 5: Time series of the latent variational distributions represented by the curves t →, µ (b,t),k and t → µ (b,t),k ± 1.96γ (b,t),k for vessels b = 2, 6, 19 and k = 1, . . . , d Z .

Figure 6 :

 6 Figure 6: Trajectories and latent space for fishing vessels 2, 6 and 19. The highlighted trajectories and points correspond to the one of June 2015 in orange and December 2015 in blue. The x-axis, y-axis and size of points in the latent space correspond respectively to µ (•,•),1 , µ (•,•),2 and µ (•,•),3 .

Figure 9 :

 9 Figure 9: Proportions of trajectories with high overlap between their latent representations (BC ≥ 0.8) for the 33 fishing vessels. On the x-axis, the proportion for adjacent days trajectories (DSI b ) and on the y-axis, the proportions for all couples of trajectories (GSI b ). Each vessel (each point) is identified by its number and the proportion of group memberships found by the colSBM (refer to the right of Figure 10) given by the pie-chart colours. The three contrasted vessels (2, 6, 19) used to illustrate the exploration of individual behaviour are written and framed in red. The black diamond shows the proportions computed for the whole fleet (DSI and GSI), while the black triangle represents the average of individual stabilities. Averages of individual stabilities for each group are shown in black-boxed squares.

Figure 10 :

 10 Figure 10: Left: Common mesoscale structure of the collection. The dark lines separate the 5 groups, the size of the blocks stands for the block proportions, while the shade of reds (alpha) gives the average number of connections between and within groups. Right: Group membership proportions (Prop) found by the colSBM model for each fishing vessel over the collection of 20 quarterly proximity graphs. The three selected vessels are highlighted (vessels 2, 6 and 19).

Figure 11 :

 11 Figure 11: (a) Top: trajectories ordered by the groups found with the colSBM model. Bottom: distribution of the mean of the latent representations (µ (•,•),1 , µ (•,•),2 ) for the same sets of trajectories. (b) Quantile distribution of the income per fishing day by quarter and by the group found with the colSBM model. The black line is the median, while the shades depict the α and 1 -α quantiles for α ∈ {.01, .1, .25}.

Figure 12 :

 12 Figure 12: Adjacency matrices of the 20 quarterly proximity graphs between fishing vessels. Matrices are ordered by quarter (row) and year (column). Proximity between vessels is printed in black and NA value in red.

Encoder:

Convolutional layer: Conv1d(2+2, 8, kernel_size=4, stride=2) Batchnorm1d( 8

B Details of the probabilistic model for proximity graphs

Let B (q,y) be the number of vessels having done at least a fishing trip on quarter (q, y). The vessels are divided into R groups, which drive their connections. Let (W (q,y) b ) = r if vessel b belongs to group r ∈ {1, . . . , R} in quarter (q, y). The group memberships are independently distributed with probability:

where (ρ r ) r∈ 1;R ) is a probability vector. Given their respective group memberships, the connections between each couple of vessels (b, b ′ ), b < b ′ are assumed to be