Stephane Redon

Ming C Lin

A Fast Method for Local Penetration Depth Computation

This paper presents a fast method for determining an approximation of the local penetration information for intersecting polyhedral models. As opposed to most techniques, this algorithm requires no specific knowledge of the object's geometry or topology, or any pre-processing computations. In order to achieve real-time performance even for complex, nonconvex models, we decouple the computation of the local penetration directions from the computation of the corresponding local penetration depths: for any pair of intersecting objects, we partition the penetrating zones into coherent regions, and we determine for each of these regions a local penetration direction. Then, for each of these regions, we estimate a local penetration depth in the previously computed penetration direction. This method has been implemented and tested on a 2.0 GHz Pentium PC with a NVIDIA GeForce FX 5900 graphics card with AGP 4× and 768MB of RAM. The examples indicate that a meaningful local penetration information can be computed for complex models and challenging intersection scenarios within a few milliseconds.

Introduction

The penetration depth information, i.e. the minimum translation required to separate two objects, is often useful for computing an appropriate collision response, or for backtracking simulation to the first time of contact. The exact penetration depth, however, is typically very difficult to compute, and several approximate methods have been introduced in the past. Although efficient, these methods often require some specific knowledge of the objects geometries (e.g. convexity) or topologies (e.g. closed manifold), or special pre-processing computations [Cam97, GS87, KR92, Ber01, KLM02, OJSL04]. Moreover, some of these methods might not always be fast enough for real-time interaction with complex models.

In this paper, we present a fast method that computes an approximation of the local penetration information and requires no pre-processing nor specific knowledge of the objects' geometry and topology. We only assume that a collision detection function has processed a pair of objects made of triangles (general triangle soups with no connectivity information), and has output a list of pairs of intersecting triangles as well as the corresponding intersection segments.

Essentially, we simplify the problem by attempting to compute meaningful local penetration information. Moreover, we decouple the computation of the local penetration directions and the corresponding local penetration depths: first, we partition the penetrating zones into coherent regions, and we determine a local penetration direction for each of these regions. Then, for each of these regions, we estimate a local penetration depth along the computed penetration direction.

We do not attempt to compute the exact penetration depth, but to provide estimates of local penetration depth as efficiently as possible. For this reason, we make several simplifying assumptions and design heuristics which work best when these assumptions are applicable, and still provide meaningful estimates in degenerate or complex cases (e.g. non-closed objects). We believe that, although not exact, our algorithm provides a sufficiently good approximation of the local penetration directions and depths, especially when interactive performance is crucial.

The overall pipeline of our algorithm consists of three stages:

1. Intersection segment clustering. We assume that the collision detection algorithm determines a list of intersection segments. The first stage clusters these intersection segments into intersection curves, which define locally coherent penetration regions.

2. Intersection curve fitting and clustering. When the intersection segments have been clustered into intersection curves, we perform a final clustering and estimate an local penetration direction for each intersection curve.

3.

Local penetration depth computation. For each intersection curve and its corresponding local penetration direction, we determine a local penetration depth, i.e. the minimum signed distance required to locally separate the objects in this direction.

Figure 1 shows the penetration information computed by our algorithm on two intersecting Stanford bunnies (48,000 triangles each). In this example, the penetration information is computed in less than three milliseconds.

Intersection segment clustering

This section describes the first stage of our algorithm to compute the local contact information. We assume that the collision detection algorithm has output a list of n s intersection segments, where each segment corresponds to a nondegenerate intersection between two non-coplanar triangles. This assumption corresponds to the usual output of a triangle/triangle intersection test [START_REF] Moller | A Fast Triangle-Triangle Intersection Test[END_REF]. More precisely, each Figure 1: Two polyhedral bunnies (48,000 triangles each) intersect. Our algorithm clusters the intersection segments into a single intersection curve (left, highlighted in the center part), and determines the corresponding penetration direction and depth, which allows to separate the bunnies (right), in less than three milliseconds.

intersection segment is described by its two endpoints p 1 and p 2 .

In the first stage of our algorithm, we attempt to determine locally coherent penetration regions by clustering these intersection segments into intersection curves. Since we do not make any assumptions on the object's topology, there is no guarantee that the intersection segments should be or even can be clustered into coherent intersection curves. This is typically the case for models consisting of non-manifold surfaces and even for "soups of polygons". Thus, we design a technique that works well in practice for most cases.

Even for simple objects, the problem of clustering the intersection segments into intersection curves is difficult, mainly because we do not assume any knowledge of the object geometry or topology. In particular, the intersection curves could be open, if the intersecting objects themselves are not closed, or contain loops. Moreover, due to finite precision computations, two adjacent intersection segments on an intersection curve might not have a common endpoint. Also, some numerical issues in the collision detection stage might lead to some missing intersection segments, and thus result in incomplete sets of intersection segments for the penetration depth computation stage.

The simple yet effective technique we present proceeds as follows. Assume c -1 curves have already been determined, and we want to determine a c-th curve from the remaining intersection segments, i.e. the free intersection segments. One of these free intersection segments is selected and arbitrarily oriented (i.e. the beginning and the end of the segment are arbitrarily set). This segment is the first one in the c-th curve, and is also the current end segment of the curve. We build the curve by repeatedly appending a free intersection segment to the current end segment of the curve. More precisely, we append the free intersection segment for which one of the two endpoints is the closest to the end endpoint of the current end segment. This closest segment is removed from the list of the free segments, and becomes the new current end segment. The construction of the c-th curve continues as long as there exists a free intersection segment that is sufficiently close to the current end segment. When the construction of the c-th curve stops, the algorithm stops, provided that there is no free intersection segment left. Otherwise, a new free intersection segment is selected and arbitrarily oriented, and a new curve, the (c + 1)-th one, is constructed.

Despite its simplicity, a naive implementation of this algorithm would yield a quadratic complexity in the number of intersection segments, essentially due to the proximity query step, which determines the closest segment to another segment. Another reason for the quadratic complexity is in the management of the list of free intersection segments, when proper care is not taken. This is not practical when many intersection segments have been output by the collision detection functions, as we must perform the penetration depth estimation in real time for many interactive applications (e.g. games, VR, etc).

In order to achieve a nearly linear runtime complexity in the number of intersection segments, we use a combination of a heap and a hashtable. The heap is used to manage the list of free intersection segments, and is implemented as a doubly-linked list stored in an array. The hashtable is used to perform the proximity queries in nearly linear time using a spatial partitioning method. Basically, we subdivide the region surrounding the intersection segments into uniform grids and we store the non-empty space cells in a hashtable, to reduce the memory consumption. Although the method is well-known, we carefully design the data structures so as to be able to perform insertion in the hashtable in constant time, and the proximity queries in nearly constant time (cf. online appendix1 ). Using these data structures, the clustering of the intersection segments can be performed in nearly linear time in the number of intersection segments. When a new set of intersection segments has to be clustered, each segment is stored in the heap, while each of its endpoints is stored in the hashtable. When the algorithm tries to determine the closest free endpoint to a specific endpoint e, the hashtable is used to traverse the list of endpoints which are contained in the same grid cell as e and determine the closest one. Assuming each grid cell contains at most a few endpoints, the clustering is performed in linear time.

Intersection curve fitting and clustering

When the first stage of our algorithm completes, the intersection segments have been clustered into intersection curves. We then proceed to the second stage of our algorithm, in which we perform a final clustering and estimate an average penetration direction for each intersection curve.

This second stage is performed using oriented bounding boxes (OBBs). For each of the intersection curves output by the first stage of our algorithm, we determine a fitting OBB using principal component analysis [GLM96]. For each intersection curve, a 3 × 3 covariance matrix of the positions of the curve's endpoints is computed. This matrix is diagonalized using an iterative method [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF], and the eigenvectors of the covariance matrix indicate the OBB axes. The OBB center and dimensions are then computed to obtain the smallest OBB which has these axes and contains the curve's endpoints.

In a first step, the OBBs are used to perform a final clustering of the intersection curves into larger intersection curves. This final clustering helps connecting curve segments that were not connected during the previous stage, either because of missing intersection segments or because two endpoints that should have been connected were stored in two different grid cells. Basically, we repeatedly merge two intersection curves as long as the directional L 2 distance between their fitting OBBs is below some threshold OBB (cf. Section 5). Whenever two curves are merged into a new curve, we determine a new fitting OBB around the new curve. The process stops when no curves can be merged. Although this step is quadratic in the number of intersection curves, the number of curves at the beginning of stage two is typically very small (as opposed to the number of intersection segments), thus merging curves is typically very fast.

When the final set of intersection curves is obtained, we define the average penetration direction for each curve to be the axis of the fitting OBB which corresponds to the smallest dimension of the OBB, i.e. the axis with the smallest associated eigenvalue. The rationale for this choice is that, for simple intersection geometry and a small amount of penetration, the actual penetration direction, i.e. the one corresponding to the direction of the minimal translation required to separate the objects, is well approximated by the normal of the plane fitting the intersection curve (cf. Section 5).

Local penetration depth computation

As stated in the introduction, the efficiency of our approach comes from the fact that we address the problem of estimating the minimal penetration depth by decoupling the estimation of the penetration direction and the determination of the corresponding directional penetration depth. The first two stages of our algorithm are crucial in helping to determine potential local coherent penetration regions and corresponding penetration directions, by making use of the only knowledge available: the intersection segments.

The third and final stage of our algorithm simply deter- mines a directional penetration depth for each intersection curve. When the directional penetration depth has been determined for each curve, the penetration depth per intersection segment endpoint is set to be the directional penetration depth of the belonging curve, offset by the signed distance of the endpoint to the curve's fitting plane.

In order to estimate the directional penetration depths, we use the graphics hardware. Essentially, for a given pair of objects and a given intersection curve, we render the objects in the penetrating region and perform readbacks of the zbuffer to determine the maximal penetration depth at pixel resolution, in the direction previously determined.

Assume the penetration direction is ±n. The direction sign is still undetermined because we do not know which one of the two directions +n or -n would require the smallest translation to separate the objects. We thus estimate the penetration in both directions and select the smallest one.

Assume we begin by estimating the penetration depth in direction d = +n. We proceed in two stages. First, we set the center of the viewing frame to be the center of the OBB which fits the intersection curve, and we set the view direction to be d. We then render the first object and read the z-buffer back into a first two-dimensional array. Then, we reverse the view direction to -d and render the second object, and read the z-buffer back into a second two-dimensional array. In each case, the view limits are set such as to bound the intersection curves and contain the whole object behind the fitting plane (cf. Figure 2). The directional penetration depth is set to the maximum of the pairwise differences of the two penetration depth arrays. The efficiency of this last stage essentially depend on the complexity of the objects to render, and the size of the z-buffer region used to determine the penetration depth. The size of the z-buffer also determines the precision of the penetration depth determination (cf. Section 5).

Implementation and Examples

Parameters

Our algorithm uses one tunable parameter per stage, whose value can be set depending on the application, the objects, and the desired performance. In the following, we indicate how these parameters might be set.

• Number of grid cells. The first parameter occurs in

the first stage of our algorithm. If the number of cells is too small, then the performance of the first stage tends to become quadratic because of the proximity query step. In contrast, if the total number of grid cells is too high, then two segments which belong to the same curve might not be connected (especially if some segments are missing). In order to determine the number of grid cells automatically, we note that a curve is mono-dimensional by definition, and thus there is on average a linear relationship between the number of intersection segments n s and the number of grid cells for each side of the bounding box. Consequently, we express the number of grid cells for each side n g as a linear function of the number of intersection segments: n s = kn g (total number of cells in the grid is n 3 g ). Our preliminary tests have shown that a good value for k is 50 for our benchmarks.

• OBB/OBB distance threshold. This parameter determines when two intersection curves should be merged.

Discussion

If the threshold is too large, there is a risk that two curve segments belonging to the same curve would not be detected. However, if the threshold is too small, two intersection curves might be merged without a need for it. In both cases, the quality of the penetration depth direction is degraded. Depending on the application, this parameter might be used to control the locality of the contact information. Assuming that on average at most one intersection segment is locally missing in a curve, the distance threshold can be set to the average length of an intersection segment, which can be precomputed from the average size of a triangle in the models involved in the application.

• Z-buffer size. Both the precision and the performance of the third stage depend on the size of the region used to render the objects and determine the local penetration depth. If this region is too small, there is a risk that the penetration depth won't take into account the sharp features of the penetration regions. In contrast, if this region is too large, it might take too long to read back. In our tests, we have used a constant buffer size of 64x64 (128x128 for the entire test, which involves four renderings). We are currently investigating different methods to determine the region size adaptively. One of the methods can be simply setting it proportionally to the projected size of the penetration region, in order to have a constant precision when the objects are rasterized. We note that if reading-back the depth-buffer becomes a bottleneck, occlusion queries can be used to perform the depth estimation directly on the graphics hardware [OJSL04].

Benchmarks

This section presents various examples obtained with our algorithm. The benchmarks have been carried out on a 2.0GHz Pentium PC with a NVIDIA GeForce FX 5900 graphics card, AGP4x, and 768MB of RAM. Several polyhedral models have been used for the tests, and we report the most challenging ones performed on a pair of Stanford bunnies (48,000 triangles each, cf. Figure 1)and a pair of Stanford dragons (about 300,000 triangles each, cf. Figure 3.c). These complex models give rise to several interesting intersection scenarios.

Figure 3 describes three such general scenarios. In all these cases, our algorithm determines the appropriate intersection curves and provides meaningful contact information, which allow us to separate the objects. The intersection curves within the OBBs are displayed, as well as the local penetration vectors. We note that, even in case of multiple intersection regions (Fig. 3(a)) or in case of degenerate intersection curves (such as the open curve in Fig. 3(b)), our algorithm is able to return useful penetration information.

In order to determine the performance of our algorithm in practice, we have generated random intersecting configurations for the two bunnies, and measured the time taken by each stage in each case.

Figure 4 reports the time required to cluster the intersection segments into intersection curves (in msec), depending on the number of intersection segments. Figure 4 clearly indicates that the carefully implemented data structures allow us to perform the clustering in linear time, depending on the number of intersection segments. Moreover, the timings in Figure 4 also show the efficiency of the first stage: one thousand intersection segments can be clustered in roughly one millisecond.

The timings of the second stage do not exhibit any simple structure, because they depend on the number of segments per initial curve, on the number of initial curves, and on the number of curves after the final clustering. For the same set of intersecting configurations, we have observed that the average time required by the second stage is 0.33 msec, with a worst case slightly less than 1 msec.

Figure 5 shows the time required by the third (also last) stage of our algorithm in milliseconds (msec), depending on the number of intersection curves to process. The graph clearly reports the expected linear-time performance, when the depth-buffer size and the pair of objects do not change. On average, the time taken by stage three for the two Stanford bunnies is about 2 msec per intersection region. Note that this includes the time to set the buffers, the view-points, render the objects, read the depth-buffer back, and perform the per-pixel depth computation.

Discussion

We have described a fast method to determine meaningful local penetration directions and depths very efficiently. Our algorithm has the following limitations:

• Locality. Our algorithm only determines local penetration information, which could be a potential problem for consistency in collision response. But, it might be pos- sible to use information acquired at previous time steps to alleviate the problem.

• Topology assumptions. We have assumed that the local intersection between two objects is a non-degenerate polygonal curve with no loops. Thus, our algorithm does not handle cases where more than two intersection segments endpoints meet at the same point (more than one curve will be reported).

• Penetration direction approximation. Our algorithm assumes that the local penetration direction can be well approximated by the smallest axis of the OBB which bounds the intersection curve.

• Penetration distance approximation. The directional penetration depth is approximate. In particular, there is no guarantee that the local penetration can be completely removed by only using the local penetration information returned by an approximation algorithm. In the example of Figure 6, some portion of the penetrating cylinder is not fully contained in the bounding box used to estimate the penetration depth in the final stage of our algorithm, and the full penetration is not detected.

• Image-precision errors. Our algorithm makes use of the graphics hardware to determine the directional penetration depth, up to pixel resolution. This can prevent fine, sharp features to be taken into account when determining the penetration depth.

In practice, however, we observed that our algorithm returns meaningful and useful penetration information even for difficult intersection scenarios involving complex, general (nonconvex) objects. Furthermore, the algorithm was able to determine local penetration information within just a few milliseconds.

Figure 6: In this deep penetration example, our method cannot retrieve the exact penetration depth, because some portion of the penetrating cylinder is not fully contained in the bounding box used to estimate the penetration depth in the final stage of our algorithm.

Web Information

Details of our data structures for intersection segment clustering are described in an appendix available at http://www.acm.org/jgt/papers/RedonLin06, and on the project web page at http://gamma.cs.unc.edu/LFPD.

Figure 2 :

 2 Figure 2: The penetration depth is determined by rendering both objects in both potential penetration directions. (a) First object, direction +n. (b) second object, direction -n. These first two renderings enable us to determine the penetration depth in direction d = +n. (c) first object, direction -n. (d) second object, direction +n. These last two renderings allow us to determine the penetration depth in direction d = -n. The actual penetration depth is the smallest of the two rendering passes.

Figure 3 :

 3 Figure 3: Three penetration depth scenarios. (a) Even in case of a complex intersection scenario, our algorithm determines a satisfying partitioning of the intersection segments, and the corresponding local penetration information. (b) Holes in the base of the bunnies (left) lead to open intersection curves appropriately detected by our algorithm (center), still allowing to precisely separate the objects (right). (c) Two dragon models (300,000 triangles each) intersect (left). Our algorithm clusters 1,491 intersection segments in a single curve (center) and determines a penetration information which allows to precisely separate the objects (right), in 9.5 milliseconds.

Figure 4 :

 4 Figure 4: Total time (in msec) spent in clustering the intersection segments depends on the number of intersection segments. Carefully designed data structures allow us to perform the clustering in linear time.

Figure 5 :

 5 Figure 5: Total time (in msec) spent in the third stage of our algorithm and the determination of the local penetration depth depend on the number of intersection curves to process. When the size of the depth-buffer is fixed, the computation time per intersection curve is nearly constant.

The data structure details are in an appendix available online at the web site listed at the end of the paper.Journal of Graphic Tools

Journal of Graphic Tools

Acknowledgements

The authors wish to thank Stanford University for the bunny and dragon models, Miguel A. Otaduy for insightful discussions and his reviewing of a first version of this paper, and the reviewers for their thorough and useful comments on improving this paper. This project is supported in part by the Office of Naval Research VIRTE Program, National Science Foundation, U.S. Army Research Office, and Intel Corporation.


Jacky Cresson

Anna Szafrańska

DIFFUSION EQUATIONS WITH SPATIALLY DEPENDENT COEFFICIENTS AND FRACTAL CAUER-TYPE NETWORKS

published or not. The documents may come

Diffusion equations with spatially dependent coefficients and fractal Cauer-type networks

Diffusion equations with spatially dependent coefficients were recently used to model powerlaw behaviors observed in Hydrogen storage experiments by V. Tartaglione in ( [START_REF] Tartaglione | Analyse de la nature stochastique et modléisation interne des systèmes générant des comportements de type puissance[END_REF],Chap. 4, 4.2.4). The type of diffusion equation used in this work is very particular and is related under a specific spatial discretization to Cauer type networks (see [START_REF] Tartaglione | Analyse de la nature stochastique et modléisation interne des systèmes générant des comportements de type puissance[END_REF],Chap.2,2.6). It was first discussed by J. Sabatier in [START_REF] Sabatier | Beyond the particular case of circuits with geometrically distributed components for approximation of fractional order models: application to a new class of model for power law type long memory behaviour modelling[END_REF]. These models seems to provide a good alternative to more conventional modeling using fractional models, meaning particular form of the transfer function admitting a fractional form (see [START_REF] Tartaglione | Analyse de la nature stochastique et modléisation interne des systèmes générant des comportements de type puissance[END_REF], equation (4.10)).

The aim of the present work is to provide complete proofs and a self-contained presentation for relations and results presented in [START_REF] Sabatier | Beyond the particular case of circuits with geometrically distributed components for approximation of fractional order models: application to a new class of model for power law type long memory behaviour modelling[END_REF][START_REF] Tartaglione | Analyse de la nature stochastique et modléisation interne des systèmes générant des comportements de type puissance[END_REF]. In particular, we are interested in a global theorem which precises the connection existing between Cauer-type networks and in particular fractal Cauer-type networks and diffusion equation with spatially dependent coefficients. Moreover, we also prove the functional relations satisfied by the transfer function of a fractal Cauer-type network allowing us to characterize the power-law type behaviors which can be expected.

The paper is organized as follows: In Section 2, we define the class of diffusion equations under study as well as the equations obtained under a spatial discretization. Section 3 deals with Cauer type ladder networks and the connection with discretizations of specific diffusion equations. Section 4 gives an explicit expression for the transfer function of a Cauer type ladder network and section 5 specify all the previous results for Fractal Cauer type ladder networks. Section 6 is devoted to a characterization of the fractional behavior of fractal Cauer type ladder network and Section 7 gives some perspectives of this work.

Diffusion equation with variable coefficient

Let u : R + × R + → R, u = u(t, z)
, is a solution of the following diffusion equation with spatially dependant coefficients:

(1)

∂ t u = γ(z)∂ z β(z)∂ z u
where γ : R + → R, β : R + → R and let u(0, z) = u 0 (z).

Taking the Laplace transform of (1)

(2)

U (s, z) = L (u)(s),
where L is the Laplace transform, we obtain

(3) sU = γ(z)∂ z β(z)∂ z U .
In [START_REF] Cresson | From fractal R-L ladder networks to the diffusion equation[END_REF], following a previous work of A. Oustaloup [START_REF] Oustaloup | Diversity and Non-integer Differentiation for System Dynamics[END_REF], we have proved that discretization over a geometric space-scale of the inverse Fourier transform of the diffusion equation with constant coefficient can be interpreted as a fractal RL ladder network.

In the following, we give the properties satisfied by a space-discretization of the diffusion equation (3) over a uniform space-scale.

2.1. Space discretization of diffusion equation. Let h > 0, we denote by T z the uniform space-scale defined by the set of points z k = kh, k ∈ N. Let φ(s, z) = β(z)∂ z U , then we consider the following space-discretizations φ h , U h of functions φ and U (4)

φ h (s, kh) = β(kh) 1 h U h (s, (k + 1)h) -U h (s, kh) sU h (s, kh) = γ(kh) 1 h φ h (s, kh) -φ h (s, (k -1)h) .
Therefore, (φ h , U h ) satisfy for all k ∈ N, k ≥ 1, the following relations

(5) φ h (s, kh) = β(kh)∆ + [U h ](s, kh) sU h (s, kh) = γ(kh)∆ -[φ h ](s, kh),
where ∆ + , ∆ -are the forward and backward discrete operators defined for all g ∈ C(T z , R d ) by ( 6)

∆ + [f ](z k ) = f (z k+1 ) -f (z k ) h and ∆ -[f ](z k ) = f (z k ) -f (z k-1 ) h ,
respectively.

The space-discretization of equation ( 3) over the space scale T z is then given by ( 7)

sU h (s, kh) = γ(kh)∆ -β(kh)∆ + [U h ] (s, kh), k ≥ 1.
The next section gives formal expression for the recursive relation between voltage and current in a Cauer-type ladder network. We will see that the discrete form of the diffusion equation ( 7) over the space scale T z can be interpreted as a specific Cauer-type ladder network. We introduce two kinds of the basic relations of the impedance with the current I and the voltage U .

Cauer type ladder networks

• series type relation:

S 1 S 2 U = U 1 + U 2 = ZI, Z = Z 1 [S 1 ] + Z 2 [S 2 ] I U 1 U 2 U
• parallel type relation:

3.2. Cauer type RC networks. Let us introduce the recursive synthesis of the Cauer type RC network.

According to Figure 1 we can observe that the RC circuit can be described by the recurrent pattern presented on Figure 2.

We use the classical current and voltage relationships in the components of RC electronic circuit with series-parallel topology (Figures 3 and4).

S 1 S 2 U = Z(I 1 + I 2 ) = ZI, 1 Z = 1 Z 1 [S 1 ] + 1 Z 2 [S 2 ] I I 1 I 2 U R 1 R 2 C 1 S 3 C 2 I 0 I 1 I 2 S 1 S 2 U 0 U 1 U 2 Figure 1.
Recursive synthesis of the RC circuit with the series-parallel topology

R k S k+1 C k I k-1 I k U k-1 Figure 2. Recurrent pattern of RC circuit Let U ∈ C(R + × N; R), I ∈ C(R + × N; R), R ∈ C(N; R), C ∈ C(N; R).
We consider a Cauer type RC network (Figure 2) with

I k = I(s, k), U k = U (s, k), R k = R(k), C k = C(k). C k Z[C k ](I k-1 -I k ) = U k , k ≥ 1 I k-1 -I k U k Figure 3. Current-voltage relationship in capacitor R k U k-1 -U k = Z[R k ]I k-1 , k ≥ 1 I k-1 U k-1 -U k
Z[R k ](s) = R k , Z[C k ](s) = 1
sC k for all k ≥ 0. Then the classical relations for elements behaviour of the Cauer type RC network can be written in the form:

(8) - 1 R σ (k) ∆ + [U ](s, k) = I(s, k), k ≥ 0, - 1 C(k) ∆ -[I](s, k) = sU (s, k), k ≥ 1.
where R σ (k) = R(k + 1). From (8) we derive for all k ≥ 1

(9) sU (s, k) = - 1 C(k) ∆ -- 1 R σ (k) ∆ + [U ](s, k) (s, k).
This equation is very close to the discrete form of a diffusion equation with spatially dependents coefficients (7) discuss in the previous section. We precise this connection in the next Section.

3.3.

Interpretation of a diffusion equation as a Cauer type network. Comparing (7) and (9), we see that (7) can be obtained as a Cauer type network by choosing appropriately the functions R and C in C(N, R).

Precisely, let T z = {z k = kh} k≥0 . We introduce the mapping π h : N → T z defined for all k ∈ N by π h (k) = kh. We have the following diagram

T z N R f ∈ C(T; R) π h f • π h ∈ C(N; R)
allowing to associated to any function f ∈ C(T z , R) a function in C(N, R).

We then choose R and C such that the following relations are satisfied

(10) - 1 C = β • π h , β ∈ C(T, R), (11)
- 1 R σ = γ • π h , γ ∈ C(T, R).
We want to select specific Cauer-type ladder networks which are able to exhibits fractional behaviors. In order to do that, we first derive an explicit form of the transfer function associated to a Cauer type ladder networks.

H(s, k) = I(s, k) U (s, k) , k ≥ 0.
Using (8) we obtain for all k ≥ 0

(12) H(s, k) = 1 R(k + 1) 1 + 1 sR(k + 1)C(k + 1) 1 + 1 sC(k + 1) H(s, k + 1)
.

Denoting by Z C (s) = Z[C](s) = 1 sC , (12
) can be rewritten as

(13) H(s, k) = 1 R(k + 1) 1 + Z C(k+1) (sR(k + 1)) 1 + Z C(k+1) (s)H(s, k + 1)
.

Continued fraction.

In order to write clearly the recursive formula, we introduce the following form of the continued fractions:

(14) [a 1 , . . . , a k ] = a 1 1 + a 2 1 + a 3 1 + • • • 1 + a k-1 1 + a k
The continued fraction (14) can be rewritten by construction as:

(15) [a 1 ,
(17) H(s, 0) = 1 R(1) , Z C(1) (R(1)s), Z C(1) (R(2)s), . . . , Z C(k) (R(k)s), Z C(k) (R(k + 1)s), . . .
Proof. Using ( 13) and ( 14) we can write

H(s, 0) = 1 R(1) , Z C(1) (R(1)s), Z C(1) (s)H(s, 1)
.

and for k ≥ 1 we have

(18) H(s, k) = 1 R(k + 1) , Z C(k+1) (R(k + 1)s), Z C(k+1) (s)H(s, k + 1)
.

Replacing H(s, 1) in H(s, 0) by expression (18) for k = 1, we have:

H(s, 0) = 1 R(1) , Z C(1) (R(1)s), Z C(1) (s) 1 R(2) , Z C(2) (R(2)s), Z C(2) (s)H(s, 2) .
Form (15) and from

1 α Z C(1) (s) = Z C(1) (αs) we obtain H(s, 0) = 1 R(1) , Z C(1) (R(1)s), Z C(1) (R(2)s), Z C(2) (R(2)s), Z C(2) (s)H(s, 2)
.

By recurrence we obtain the form (17). □

Fractal Cauer type ladder networks

Fractal Cauer type ladder networks are specific type of Cauer ladder networks which can be defined using few parameters and exhibiting a rich dynamics. We give an explicit form for the transfer function of a Fractal Cauer type ladder network. A Cauer type ladder network is said fractal if there exists two constants σ and ρ, such that we have the scaling relation

(19) R(k + 1) = σR(k), C(k + 1) = ρC(k).

Fractal Cauer transfer function.

Then we can simplify the expression of the transfer function H(s, 0) given by ( 17), as follows

(20) H(s, 0) = 1 R 1 , Z(s), 1 σ Z(s), 1 ρσ Z(s), . . . , 1 ρ n σ n Z(s), 1 ρ n σ n+1 Z(s), . . . with Z(s) = Z C(1) (R(1)s).
(21) r((k + 1)h) = σ(h) r(kh), c((k + 1)h) = ρ(h) c(kh), ∀ k≥0 .
These two conditions are of the form: for a given ς : R → R find f : R → R such that for all h > 0 we have

(22) f ((k + 1)h) = ς(h)f (kh), ∀ k≥0 .
From ( 22) for k = 0 we have f (h) = ς(h)f (0), and we obtain a definition of function ς in ( 22):

(23) ς(h) = f (h) f (0) Posing g(ζ) = f (ζ) f (0)
, the equation ( 22) with ( 23) reads as follows (24) ∀ h>0 g((k + 1)h) = g(h)g(kh), ∀ k≥0 .

This functional relation has very particular solution:

Theorem 5.1. A continuous real valued function g is a solution of (24) if anf only if it there exists a constant λ ∈ R such that

(25) g(ζ) = e ζλ .
Proof. First we proof that the functional relation (24) is satisfied. This proof consists with three steps:

1) ∀ k,k ′ ∈N g((k + k ′ )ζ) = g(kζ)g(k ′ ζ); 2) ∀ a,b∈Q g(a + b) = g(a)g(b); 3) ∀ u,v∈R g(u + v) = g(kζ)g(k ′ ζ);
The first equality is proved by induction. Indeed, we have

(26) g((k + k ′ )z) = g((k + k ′ -1)z + z) = g((k + k ′ -1)z)g(z).
Replacing g((k

+ k ′ -1)z) by g((k + k ′ -2)z)g(z) we obtain (27) g((k + k ′ )z) = g((k + k ′ -2)z)g(z)g(z). But g(z)g(z) = g(2z) so (28) g((k + k ′ )z) = g((k + k ′ -2)z)g(2z)
.

By iteration, we obtain for all i = 0, . . . , k ′ , that

(29) g((k + k ′ -i)z + iz) = g((k + k ′ -i)z)g(iz).
In particular, for i = k ′ , we obtain

(30) g((k + k ′ )z) = g(kz)g(k ′ z).
As z can be negative, the result extend to k, k ′ ∈ Z.

The second equality is a simple computation using equality 1). We consider two rational numbers of the form m/n and p/q with n, q ∈ N * and m, n ∈ N. We have

(31) g m n + p q = g mq + np nq = g (mq + np) 1 nq
Using equality 1), we then obtain

(32) g m n + p q = g mq 1 nq g np 1 nq = g m n g p q .
The last equality is proved by density using the continuity of g. Let u, v ∈ R + and u n , v n , n ∈ N be two sequences of rational numbers converging to u and v respectively. Then, for all n ∈ N, we have by equality 2) that ( 33)

g(u n + v n ) = g(u n )g(v n ).
Taking the limit when n goes to +∞ and by continuity of g, we obtain for all u, v ∈ R that (34)

g(u + v) = g(u)g(v).
As g is a positive function, the previous equality means that g is a continuous morphism of the additive group (R, +) to the multiplicative group (R + , ×). It is known in this case that g must be of the form g(z) = e λz for a given λ ∈ R. □

Using the previous Theorem, we obtain the following structure result:

Theorem 5.2. Continuous positive real valued functions R and C satisfy equality (10) and (11), respectively, if and only if there exists two real constants λ R and

λ C such that R(z) = R 0 e λ R z and C(z) = C 0 e λ C z
As a consequence, Fractal Cauer type networks are associated to diffusion equations characterized by β(z) = -β 0 e -λ R z and γ(z) = -γ 0 e -λ C z .

We then obtain the following connection between Fractal Cauer type networks and generalized diffusion equations: Theorem 5.3. Fractal Cauer type networks with characteristics λ R and λ C are in one to one correspondence for all h > 0 with the discretization in space along the space-scale {kh}, k ∈ N of the diffusion equation with non-constant diffusion coefficient

(35) ∂ t u = -γ 0 e -λCz ∂ z -β 0 e -λ R z ∂ z u .
corresponding to an approximation of (39).

In the next Section, we characterize the fractional behaviors compatible with the functional relation (45). 6.2. Fractional behaviors. We look for solutions of (45) of the form (46) g(Z(s), σ, ρ; ϵ) = K(σ, ρ; ϵ)Z(s) n(σ,ρ) .

Our goal is to determine the possible values for the exponent n(σ, ρ).

Replacing g by its expression (46) in (45), we obtain:

(47) K(σ, ρ; ϵ)K(ρ, σ; ϵ)Z(s) n(σ,ρ)+n(ρ,σ) σ -n(ρ,σ) = Z(s) 1 + ϵ .
As a consequence, we obtain by identification the following system:

(48) K(σ, ρ; ϵ)K(ρ, σ; ϵ) = σ n(ρ,σ) 1 + ϵ , n(σ, ρ + n(ρ, σ) = 1.
The first equation of (48) must be valid by exchanging the role of σ and ρ due to the symmetry of the left hand side. As a consequence, we must have:

(49) σ n(ρ,σ) 1 + ϵ = ρ n(σ,ρ) 1 + ϵ .
As n(ρ, σ) = 1 -n(σ, ρ) by the second equation of (48), we deduce by replacing in (49) that (50) σ = (σρ) n(σ,ρ) .

We then deduce the following theorem: Theorem 6.2. Fractional behaviors of the form (46) satisfying the simplified functional relation (45) possess a fractional exponent given by (51) n(σ, ρ) = ln(σ) ln(σρ) .

The previous result appear in ([6],Appendix 2, (A2.28)) as a consequences of relations satisfied by a formal expansion of the function g (see [START_REF] Sabatier | Beyond the particular case of circuits with geometrically distributed components for approximation of fractional order models: application to a new class of model for power law type long memory behaviour modelling[END_REF],Theorem 1).

Conclusion and perspectives

The previous results can be used in a variety of situations as already stressed by J. Sabatier in [START_REF] Sabatier | Beyond the particular case of circuits with geometrically distributed components for approximation of fractional order models: application to a new class of model for power law type long memory behaviour modelling[END_REF] and V. Tartaglione et al. in [START_REF] Tartaglione | Analyse de la nature stochastique et modléisation interne des systèmes générant des comportements de type puissance[END_REF][START_REF] Sabatier | Fractional Behaviours Modelling. Analysis and Application of Several Unusual Tools[END_REF] covering application to the modeling of diffusion in materials used for the storage of Hydrogen. The functional relations satisfied by the transfer function of a Fractal Cauer type ladder networks is very specific and the previous discussion is only a first step in its study. In particular, two mathematical directions deserve more attention:

• J. Sabatier in ( [START_REF] Sabatier | Beyond the particular case of circuits with geometrically distributed components for approximation of fractional order models: application to a new class of model for power law type long memory behaviour modelling[END_REF],Theorem 1) and V. Tartagione et al. in [START_REF] Sabatier | Fractional Behaviours Modelling. Analysis and Application of Several Unusual Tools[END_REF] discuss a result concerning the asymptotic behavior of the transfer function using an explicit expression (see [START_REF] Sabatier | Fractional Behaviours Modelling. Analysis and Application of Several Unusual Tools[END_REF], Theorem 7.1. and details in Appendix B). The proof is based on an assumed form for the transfer function which is not proved in the paper. A complete and detailed discussion of this result is then needed. • J. Sabatier [START_REF] Sabatier | Beyond the particular case of circuits with geometrically distributed components for approximation of fractional order models: application to a new class of model for power law type long memory behaviour modelling[END_REF] discuss distributions, i.e. relations between the parameters of the Cauer type networks or equivalently different space scale discretization of a diffusion equation with spatially dependent coefficient, which are going beyond geometric distribution (see [START_REF] Sabatier | Beyond the particular case of circuits with geometrically distributed components for approximation of fractional order models: application to a new class of model for power law type long memory behaviour modelling[END_REF],p.250) and leading to power law behaviour. It will be interesting to give a full characterization of the possible distributions as well as the form of the associated coefficients corresponding to the class discuss in [START_REF] Sabatier | Beyond the particular case of circuits with geometrically distributed components for approximation of fractional order models: application to a new class of model for power law type long memory behaviour modelling[END_REF]. From the application point of view, it would be nice to give full characterization of the diffusion properties of some specific material like COF (Covalent Organic Framework) or MOF (Metal-Organic Framework) (see [START_REF] Kumar Gangu | A review on contemporary Metal-Organic Framework materials[END_REF] for a review) using the previous approach.

3. 1 .

 1 Classical relations on electronic circuits. Let us denote by Z and Y the impedance and admittance, respectively. We have Z = 1 Y . In the following, by S we denote an equivalent to one of any component of electronic circuit: resistance (R), inductance (L), capacitance (C).

Figure 4 .

 4 Figure 4. Current-voltage relationship in resistance

4 .

 4 The transfer function of a Cauer tpe ladder network 4.1. Recursive expansion of the transfer function. Let us denote by H : C(R + × N; R) the transfer function

5. 1 .

 1 Fractal Cauer ladder networks. We consider a Cauer type ladder network with coefficients R k and L k for k ≥ 1. As usual, we identify R k and L k as the values of two functions R : N → R and L : N → R at point k ∈ N.

5. 3 .

 3 Diffusion equation associated to a fractal Cauer type network. Using (10) and (11), we want to explicit the functions β and γ such that R and C satisfy the fractality condition (19). By definition, R and C depend on h. As a consequence, the fractal constant σ and ρ are functions of h. Instead of searching β and γ, we look for two continuous functions r(ζ) and c(ζ) satisfying r(kh) = R(k), c(kh) = C(k). Then the scaling relations (19) read as follows for all h > 0

  a 2 , . . . , a k ] = [a 1 , [a 2 , . . . , a k ]] Continued fraction and the transfer function. Based on relations (15) and (16) we are able to introduce the following lemma: Lemma 4.1. The closed formula for transfer function H ∈ C(R + × N) has the form

	and we have
	(16)	α[a

1 , a 2 , . . . , a k ] = [αa 1 , a 2 , . . . , a k ].

4.3.

Acknowledgements References Acknowledgements A. Szafrańska thanks the National Science Center for the financial support in the form of the research project No. 2021/05/X/ST1/00332 and J. Cresson thanks the GDR CNRS no. 2043, Géométrie différentielle et Mécanique and the fédération MARGAUx (FR 2045) for supports and discussion with J. Sabatier, V. Tartaglione and C. Farges at the IMS Bordeaux.

Fractional behavior of fractal Cauer networks

A fractional behavior is associated to a transfer function of the form (36) H(s) = cs ν , where c and ν are two constants. The main problem is to determine if a fractal Cauer type systems exhibits a fractional behavior and if yes, to characterize the exponent.

We first derive structural relations satisfied by the transfer function of fractal Causer networks and we then determine the possible fractional behaviors.

6.1. Structural relations. Let us denote by g(z, σ, ρ) the function

Then, using equation (20), we have

We have the following functional relation satisfied by the function g: Theorem 6.1. The function g satisfies

It is important to notice the permutation between σ and ρ in the formula (39). Formula (6.1) was first stated by J. Sabatier in ([6], Property 1 p.249) without a formal proof. As we have Z(σs) = Z(s)σ -1 , we obtain (41) [Z(s)σ -1 , . . . , Z(s)(σρ) -n , Z(s)σ -1 (σρ) -n , . . . ] = [Z(σs), Z(σs)ρ -1 , . . . , Z(σs)(ρσ) -n , Z(σs)ρ -1 (ρσ) -n , . . . ], = g(Z(σs), ρ, σ).

This concludes the proof. □

The previous functional relation is difficult to study directly. We can look for a simplified version of it following an idea already used by A. Oustaloup in [[2],p.208]:

Let us assume that s is such that g satisfies

Let 1 >> ϵ > 0 be such that (44) | g(Z(s), ρ, σ) | -1 < ϵ.

We consider the simplified functional relation (45) g(Z(s), σ, ρ)g(Z(σs), ρ, σ) = Z(s) 1 + ϵ ,


Impact of contaminants on formability of mechanically recycled PET:

Application to bottle blowing process Many researchers have addressed the subject of PET recycling. They found that recycled plastic suffers degradation during the recycling process which impacts both its physical and mechanical properties (decrease of molecular mass, increasing of crystallinity and rigidity, reduction of elongation at break, loss of clarity, yellowness) [3][4][5]. Very few scientists show also the bad impact of contaminants on physical and mechanical characteristics of RPET (increasing of crystallinity, reduction of elongation at break, loss of clarity, yellowness) [6][7][8]. And just one publication [9] has explored the impact of recycling and of contaminants on the implementation of new bottles.

The main objective of this present study is to understand the impact of contaminants on the stretching behavior of RPET and to adapt blowing machine parameters. Two mechanically RPET from two different countries are characterized to determine the nature of contaminants, compared to a virgin PET resin. Thanks to optical microscopy, dark particles, referred as "black specks", are visible and are more numerous in recycled material. The recycling process and the presence of contaminants enhance crystallization ability in static conditions, which has the consequence of reducing the forming range of the material. For the three studied PET, the α-transition are very close. Nevertheless, regarding ISBM process, RPET are more easily deformed than virgin because they can absorb more energy. These first results let thinking that recycled PET present their own thermo-mechanical behavior, that will result in an adaptation of blowing machine parameters.


Clément Laroche

Madalina Olteanu

Fabrice Rossi

Pesticide concentration monitoring: investigating spatio-temporal patterns in left censored data

Keywords: Pesticide concentration monitoring: investigating spatio-temporal patterns in left censored data pesticide concentration monitoring, left censored data, change-point detection, anomaly detection, Pareto front, water pollution, prosulfocarb

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Monitoring the environmental pollution is of great interest for public authorities, important adverse health-effects being well documented nowadays [START_REF] Khopkar | Environmental pollution monitoring and control[END_REF][START_REF] Marchant | Monitoring urban environmental pollution by bivariate control charts: New methodology and case study in santiago, chile[END_REF][START_REF] Nougadère | How dietary risk assessment can guide risk management and food monitoring programmes: The approach and results of the french observatory on pesticide residues (anses/orp)[END_REF]. National health agencies are thus much concerned with monitoring ambient levels and quantifying the concentration of various pollutants in given environmental areas.

At the same time, modelling environmental pollution data is a complex issue, due to several reasons, some intrinsic to the types of data under study, some specific to the data collection process implemented in different countries. Firstly, pollutant concentration levels are measured by sensors which have generally detection and quantification limits: the corresponding data are then left-censored. Secondly, the data is usually skewed to the right, with long tails hinting high concentrations. Thirdly, in numerous situations the data is irregularly sampled because of measurement practices, and is often multivariate, since various pollutant levels are monitored. Fourthly, pollution is monitored in various locations, each location possibly using different sensors, yielding a significant spatial heterogeneity. Notice that the last two problems are specific to some applications such as the one considered in the present paper. Indeed many countries have implemented very strict data collection protocols that ensure regular measurement rates on standardized sensors for a selection of pollutants1 .

How to handle left-censored and right-skewed data is therefore one of the first aspects to consider when modelling environmental data. A rich literature has been developed on this topic during the last thirty years, and may be roughly divided into three categories of approaches: substitution methods (censored data is imputed using some values chosen a priori or via a generative model), parametric methods (maximum likelihood estimates are computed under the hypothesis that the data comes from some log-normal, Weibull, Gamma, exponential, or other log-logistic distribution), and non-parametric methods (Kaplan-Meier or hazard function estimates). Detailed reviews of the various approaches are available for instance in [START_REF] Antweiler | Evaluation of statistical treatments of left-censored environmental data using coincident uncensored data sets: I. summary statistics[END_REF]; Authority (2010); [START_REF] Canales | Methods for handling left-censored data in quantitative microbial risk assessment[END_REF]; [START_REF] Gillespie | Estimating population distributions when some data are below a limit of detection by using a reverse kaplan-meier estimator[END_REF]; [START_REF] Hewett | A comparison of several methods for analyzing censored data[END_REF]; [START_REF] Mitra | Analysis of left censored data from the generalized exponential distribution[END_REF]; [START_REF] Shoari | Toward improved analysis of concentration data: embracing nondetects[END_REF].

The second aspect to consider is spatio-temporal heterogeneity. Air pollution data has received, for instance, a great deal of attention, and several modelling approaches have been proposed in the literature. Some are based on temporal regression models combined with kriging [START_REF] Lindström | A flexible spatio-temporal model for air pollution with spatial and spatio-temporal covariates[END_REF][START_REF] Sampson | Pragmatic estimation of a spatio-temporal air quality model with irregular monitoring data[END_REF], while others use latent variables and co-clustering approaches [START_REF] Bouveyron | Co-clustering of multivariate functional data for the analysis of air pollution in the south of france[END_REF]. Nevertheless, these approaches do not include the fact that monitoring data is not normally distributed, and is usually left-censored. In the specific field of pesticide concentration monitoring, several recent papers address the spatio-temporal issue from an exploratory point of view (see for instance [START_REF] Aznar | Spatio-temporal distribution of pyrethroids in soil in mediterranean paddy fields[END_REF][START_REF] Ccanccapa | Spatio-temporal patterns of pesticide residues in the turia and júcar rivers (spain)[END_REF][START_REF] Figueiredo | Spatio-temporal variation of outdoor and indoor pesticide air concentrations in homes near agricultural fields[END_REF].

If one focuses specifically on temporal heterogeneity, a common approach to deal with it is to use a change-point based segmentation. Assuming the data is strictly stationary, conditionally to a (possibly unknown) number of change-points and their associated locations, change-point analysis aims at identifying the number of changepoints (also known asbreaks), their locations, and the characteristics of the probability distribution within each temporal segment. Widely used in a variety of applications [START_REF] Basseville | Detection of abrupt change: Theory and application[END_REF][START_REF] Chen | Parametric statistical change point analysis: With applications to genetics[END_REF][START_REF] Lévy-Leduc | Detection and localization of change-points in high-dimensional network traffic data[END_REF][START_REF] Liu | Change-point detection method for clinical decision support system rule monitoring[END_REF][START_REF] Reeves | A review and comparison of changepoint detection techniques for climate data[END_REF], and, in particular, for environmental pollution monitoring [START_REF] Costa | Change-point detection in environmental time series based on the informational approach[END_REF], change-point detection is a reference technique for time series segmentation. The present manuscript relates to the offline framework, by supposing the full data has been recorded, and the segmentation is done posterior wise. For recent and detailed reviews of the offline changepoint detection, the reader may refer to [START_REF] Bardet | Change-point detection, segmentation, and related topics[END_REF][START_REF] Truong | Selective review of offline change point detection methods[END_REF]. While the literature on change-point detection is abundant, applications to spatial data are somewhat limited. An early example of such method can be found in [START_REF] Majumdar | Spatio-temporal changepoint modeling[END_REF] while recent advances in a setting close to ours are presented in [START_REF] Chen | S3t: A score statistic for spatiotemporal change point detection[END_REF]. As far as we know, none of the existing change-point detection method for spatial data applies to irregularly sampled and sparse data (on the temporal axis).

The present manuscript tackles the issue of pesticide concentration monitoring, and introduces a new methodology which integrates both the specific left-censored distribution of the data, and the spatio-temporal context. The main goal is to identify contextual anomalies, both from a temporal and a spatial point of view. The proposed method builds on a parametric model for left-censored and right-skewed distributions, and combines it with a change-point detection step and a clustering step.

Change-point detection is used for modelling temporal heterogeneity, by assuming a piece-wise stationary distribution on the series of maximum values, for a given time resolution. It produces temporal segments in which the pesticide concentrations are assumed to follow a stationary distribution.

Clustering is then used for modelling the expected spatial homogeneity while integrating geographical constraints such as river networks, wind directions, etc. Indeed, as geological, terrain and climatic characteristics of an area can influence the dispersion of a chemical substance and on its potential use in the case of e.g. a pesticide, concentrations are expected to be somewhat correlated in small scale regions that are homogeneous in terms of influencing characteristics. Especially in the application presented here, which relates to the investigation of pollutants in surface waters, it is interesting to take into account the hydrographic structure of the region (as in e.g. [START_REF] Chen | S3t: A score statistic for spatiotemporal change point detection[END_REF]. Indeed, if a high concentration of a substance is detected at a certain point in time, traces of this substance should be found later downstream. This hypothesis is accounted for by building clusters of measuring stations according to their proximities measured via the hydrographic network.

Conditionally to the temporal segment detected by the change-point procedure, and to the spatial cluster detected by the clustering procedure, one may analyse the data and identify contextual anomalies.

The rest of the manuscript is organised as follows: in Section 2, the generative model assumed for environmental pesticide monitoring data is described; the proposed method for estimating and handling this model from observed data is detailed in Section 3; a detailed example on data collected by French authorities in Val de Loire region is fully illustrated in Sections 4 and 5.

Data collection procedure and associated generative model

We study specifically in this paper a non homogeneous data collection process for pesticide use monitoring. It is represented by a generative model with two levels. The first level, a.k.a. the fine-grain level, consists of a network made of monitoring stations, where each station is associated to an irregularly sampled time-series. The second level, a.k.a. the coarse-grain level, summarises the maximum recorded values throughout the network, for a specified temporal resolution, and assumes a piece-wise stationary distribution.

Monitoring stations network

We consider a network of monitoring stations used to collect concentration measurements at irregularly sampled instants. The stations are represented by an undirected graph G = (V, E), which vertices V = (v i ) 1≤i≤N are the monitoring stations and which weighted edges E are links between stations that are directly comparable. The aim of the graph is to represent expert knowledge about expected measurement homogeneity. When two stations are connected in G, their measurements can be compared directly: a small edge weight assumes simultaneous measurements to be close, while a large one allows for significant differences. Shortest paths in the graph can be used to compare stations that are not directly connected, using the total weight of the paths to measure non homogeneity. This approach is inspired by methods developed for signal processing on graphs [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]), but we use a dissimilarity based weighting rather than the classical similarity based one. This graph based representation is very flexible and can be used to model different types of spatial homogeneity. For instance, the focus of the present paper is the monitoring of water concentration of pesticides and thus dissimilarities between stations will be computed based on the network of rivers on which they are situated (see Section 3.2). Other modelling approaches may use a different graph considering for instance dominant wind directions relevant for air diffusion of pollutants.

Data collection

Each station v i is supposed to be associated to a time series (y ij , t ij ) 1≤j≤p i , where p i is the number of sampled data points at v i , and y ij is the concentration level of some pollutant at time t ij . All measurements y ij are left-censored by some threshold q ij , representing the quantification limit. Quantification limits depend on the machines used at each station and at each time instant, hence depend both on the station v i and on the collection instant t ij . Furthermore, quantification limits are supposed to be known, fixed quantities.

Summarising the above notations and hypotheses, a data set sampled from the stations network is given by a collection of measurements and associated quantification limits, and denoted D = (y ij , t ij , q ij ) 1≤j≤p i 1≤i≤N .

Notice that in practical applications, we expect to have a rather small number of measurements for each station, i.e. to have small values for the p i . In addition, we do not expect the measurement instants to be shared among the stations. See Section 4.1 for examples. From the complete representation of the data D, one may derive an aggregated, coarser representation. First, an adapted temporal resolution for the phenomenon at study is selected. For instance, in the case of the present paper, a daily resolution is considered. Second, the selected resolution is used to build a time series of increasing instants (τ k ) 1≤k≤K , at which at least one observation is available in the data collection. We denote t ij ∈ τ k the fact that the observation time t ij is compatible with τ k at the specified resolution, e.g. that the observation y ij was made during the day τ k .

Third, once (τ k ) 1≤k≤K has been computed, one may introduce a coarse-grain, global series, summarising the maximum values recorded within the temporal resolution with

y k = max {y ij | t ij ∈ τ k } .
(1)

For instance, for a daily aggregation level, y k is the largest value among all the measurements that took place during day τ k . Notice that (y k ) 1≤k≤K is left-censored as the consequence of the censoring of the underlying values. The quantification limit for y k is denoted q k , with

q k = max {q ij | t ij ∈ τ k } . (2)
The coarse representation of D is then

D = (y k , τ k , q k ) 1≤k≤K .
(3)

A piece-wise stationary model for the coarse-grain time series

In order to model the global use of the substance under monitoring, a piece-wise stationary generative model is introduced for the coarse data set D. The model is based on the following assumptions:

• there are L * > 0 change-points producing L * + 1 stationary intervals defined by

0 = η * 0 < η * 1 < . . . < η * L * < η * L * +1 = K; • the observations (y k ) 1≤k≤K are realisations of K independent random variables (Y k ) 1≤k≤K ; • when k ∈ [η * l-1 + 1, η * l ], Y k is distributed according to a left-censored parametric distribution Q with interval dependent parameters θ *
l and a left-censoring threshold q k , which is a known constant.

Notice that the model only accounts for the concentrations y k but not for the instants and the quantification limits which are supposed deterministic quantities.

Methods

We are interested in finding anomalies in data collected according to this spatiotemporal model. The proposed methodology combines two different homogeneity models. The temporal aspect is based on the piece-wise stationary model proposed in Section 2.3, while the spatial aspect is based on the graphical representation introduced in Section 2.1. In a first step, we estimate the parameters of the temporal model. In a second step, the homogeneity assumptions represented by the graph of stations is used to detect stations with anomalous measurements with respect to close stations in a given stationary temporal segment.

Piece-wise stationary model estimation via change-point detection

The coarse-grained data D is segmented using a change-point detection approach applied to the model introduced in Section 2.3. Since both the number and the location of the change-points are unknown, we shall optimise a penalised cost function, and seek to estimate the number of change-points L * , the change-point locations η = (η * l ) 1≤l≤L * , and the parameters θ = (θ * l ) 1≤l≤L * +1 . The estimates write as

( L, η, θ) = arg min L,η,θ C(L, η, θ; D). (4)
The penalised cost is given by

C(L, η, θ; D) = L+1 l=1 -ln L Q (θ l ; y η l-1 +1 , . . . , y η l ) + β K (L + 1)D, (5)
where L Q (θ l ; y η l-1 +1 , . . . , y η l ) is the likelihood of the l-th segment for the distribution Q, β K is the penalty to apply at the addition of new segment, and D is the dimension of the parameter vectors θ l . Notice that the penalty β K depends on the size of the data K.

For fixed values of L and of η, C(L, η, θ; D) is maximized by setting θ to the maximum likelihood estimate (MLE), θ M LE (L, θ). Thus, the optimisation problem may be further written as

( L, η) = arg min L,η C(L, η, θ M LE (L, η); D). (6)
The number of change-points L and the associated locations η are obtained by applying the PELT procedure [START_REF] Killick | Optimal detection of changepoints with a linear computational cost[END_REF], which improves the optimal partitioning approach through a lower, linear complexity. The choice of the penalty term β K is driven by the CROPS algorithm [START_REF] Haynes | Computationally efficient changepoint detection for a range of penalties[END_REF], which computes all optimal segmentations as the penalty varies over some interval. Eventually, the final penalty value is selected using an elbow rule heuristic as proposed in (Lung-Yut-Fong, Lévy-Leduc, & Cappé, 2015): segmentation scores are plotted against their corresponding number of change-points L. One looks for the number of breaks L that minimizes the sums of squares of two linear models respectively fitted on the L ≥ L and the L ≤ L. The penalty value associated to L is selected.

Spatial clustering

In any of stationary intervals identified in the previous step, the measurements are assumed to be consistent with the homogeneity assumptions represented by the graph G = (V, E). A natural way of assessing the actual regularity of the measurements would be to use graph signal processing techniques (see e.g. [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. However the irregular, unaligned, sparse and censored nature of the measurements at each station, prevents the use of such methods.

The measurements are also incompatible with techniques designed to detect anomalous clusters in a graph (see for instance Arias-Castro, Candès, & Durand, 2011).

To circumvent this problem, we propose to leverage the graphical representation to build spatial aggregates and to assess homogeneity at this aggregated level. This corresponds to clustering the stations using the graph structure. We proceed as follows.

Nodes of each connected component of the graph G = (V, E) are clustered using a Ward hierarchical clustering method implemented on the shortest path distance computed from the edge weights. Those component specific hierarchies are combined in a global one in a greedy way. The initial global clustering of V is obtained by assigning all vertices in a connected component to the same cluster. Subsequent levels of the global hierarchy are obtained by replacing the clusters of a connected component by the next refined level of the local hierarchy. At each step of the refinement, we select the component that reduce the most the inertia of the clustering. We use the standard definition of inertia given for the clustering P = (C 1 , . . . , C M ) by

W (P) = M m=1 1 |C m | v i ,v j ∈C k d 2 ij , (7)
where d 2 ij is the square of the shortest path distance in G between vertices v i and v j , and |A| denotes the cardinality of set A. Clustering with a small inertia contain clusters that group close monitoring stations according to the graph G.

To select the final clustering in the hierarchy, we use the same decision rule as 3.1. This time, the inertia of the clustering is plotted against the corresponding number of clusters M . We look for the number of breaks M * that minimizes the sums of squares of two linear models respectively fitted on the M ≥ M * and the M ≤ M * .

Notice that we rely on a simple graph clustering approach for two main reasons. Firstly, we do not expect graphs of monitoring stations to exhibit the specific characteristics of complex networks (such as very high degree vertices, small diameter, etc. see e.g. [START_REF] Newman | The structure and function of complex networks[END_REF]) that justify the use of techniques such as maximal modularity clustering (see e.g. Fortunato, 2010, for a survey). On the contrary, simpler approaches that interpret shortest paths weights as dissimilarities should be sufficient (see e.g. [START_REF] Schaeffer | Graph clustering[END_REF]. Secondly, we work on relatively small graphs with even smaller connected components and we do not face computational issues associated to hierarchical clustering.

Anomaly detection

Two types of anomalous clusters are targeted: either clusters with anomalous stations, or wholly anomalous clusters. Clusters containing anomalous stations are detected by studying the homogeneity of the measurements provided by the stations in a given spatial cluster. Anomalous clusters of stations are detected by simply pooling all measurements of each cluster to estimate the local use of the substance and detect large rates. We derive in this section two anomaly scores covering those cases.

For the first case, we need to assess the homogeneity of the measurements of the stations in a spatial cluster for a stationary time interval. As pointed out previously, the number of measurements provided by a single station is usually quite small, especially when we consider a single stationary interval. As a consequence classical distances between empirical distributions are not appropriate, mainly because the measurements of two stations do not have any value in common. Then the Kolmogorov-Smirnov statistics will be essentially driven by the number of observed values rather than the actual values, while other quantities, such as the Jensen-Shannon divergence, cannot be properly estimated (see appendix B). For this reason, we propose to use the Wasserstein w 1 distance [START_REF] Villani | Optimal transport: old and new[END_REF]. For two discrete distributions on R, it is expressed as the L 1 -distance between their cumulative distribution functions and is therefore simple to compute.

The measurement homogeneity of the clusters obtained in Section 3.2 is therefore defined as the mean within cluster empirical Wasserstein average distance of a station measurements to the others. Denoting C m the m-th cluster and |C k | the number of stations present in C m , w 1 (y i , y j ) the empirical 1-Wasserstein distance between the data of stations v i and v j , this quantity is expressed as

Wk = 1 |C m |(|C m | -1) 1≤j≤|Cm| 1≤i≤|Cm|,i =j w 1 (y i , y j ). ( 8
)
The second type of potentially anomalous clusters are simply associated to the presence of quantified measurements and high values of concentration. Thus we estimate for each spatial cluster C m the parameters of distribution Q (see Section 2.3) on the pooled measurements obtained from all the stations of the cluster during the chosen stationary interval. From those parameters, we compute a statistics, denoted Īm , used as a proxy for the intensity of the measurements (see Section 5.3 for an example). Hence we consider a low concentration to be the normal case, but we do not define a threshold between normal clusters and abnormal ones.

Each cluster C m is therefore characterised by two values ( Wm , Īm ). To select potentially anomalous clusters, we use a multi-objective optimisation approach, considering that both characteristics are equally interesting. Following [START_REF] Kießling | Chapter 28 -foundations of preferences in database systems[END_REF], we say that X k = ( Wm , Īm ) is Pareto dominated by X l = ( Wl , Īl ), and we write X m ≺ X l if and only if ( Wm < Wl ) and ( Īm ≤ Īl ) or ( Wm ≤ Wl ) and ( Īm < Īl ) .

The level 1 Pareto optimal front is the set of maximal points for ≺. Level b with b > 1 is defined recursively as the optimal Pareto front computed for the set of points that do not belong to the optimal Pareto front of levels 1, . . . , b -1. Therefore clusters in the level 1 Pareto front are remarkable in the sense that there is no other cluster with higher heterogeneity and more extreme measurements. We define these clusters as anomalous. Pareto front and levels are evaluated using the Skyline algorithm [START_REF] Borzsony | The skyline operator[END_REF][START_REF] Endres | Scalagon: an efficient skyline algorithm for all seasons[END_REF].

Data presentation

The methodology introduced in the above sections will be illustrated next using a case study on the prosulfocarb concentration (for Biotechnology Information, 2022) in Centre-Val de Loire. This chemical compound is mainly used as a herbicide in field crops, with a typical period of active use in autumn. The monitoring of its concentrations in surface waters has been subject to increasing attention due to its aquatic ecotoxicology (Agriculture & Environment Research Unit (AERU) at the University of Hertfordshire, 2021; ANSES (Agence nationale de sécurité sanitaire de l'alimentation, de l'environnement et du travail), 2018).

Time period and geographical area selection

Prosulfocarb usage was banned in France before 2007. A market re-authorisation was issued by the French Observatory on Pesticide Residues (now part of the ANSES2 ) in 2009. Since then, two modifications of the authorisation for use have been put in place, in November 2018 and in November 2019 respectively. Both changes consist in restrictions of use, one imposing specific equipment for application, the other restricting the application schedule in the presence of non-target crops next to the treated area. Motivated by these changes in regulation, the time period chosen for our study spans from January 1, 2007, to September 8, 2020.

Moreover, our study focuses on the geographical area of French Centre-Val de Loire region. Indeed, between 2009 and today, the annual mass of prosulfocarb sold in this region exploded, making it rise from the 17th most sold substance in 2009 to the 4th in 2017 (see Figure 10 in Appendix C.2). This region is also characterised by high concentrations of prosulfocarb target crops (such as the Beauce plains) (see Figure 9 in Appendix C.1). Many target crops (cereal crops) are also concentrated in the region. These two elements combined guarantee a significant use of the product in this area. Thus, we expect significant variations in concentration values in this area during this period.

Data about surface water quality in France is available from the French Biodiversity Agency (français de la biodiversité, n.d.). We collected from the site the data selected above 3 . These choices led to a data set D comprising 337 monitoring stations that performed 11,002 measurements. Each measurement is described by the monitoring station ID, the sampling date, the quantification limit (LOQ), and the concentration measurement value, if the concentration exceeds the LOQ. Indeed, as this is usually the case when measuring the concentration of a chemical substance [START_REF] Bernal | Limit of detection and limit of quantification determination in gas chromatography[END_REF][START_REF] Currie | Nomenclature in evaluation of analytical methods including detection and quantification capabilities (IUPAC recommendations 1995)[END_REF], the measurements are submitted to

• a limit of detection (LOD) which is the smallest concentration of the substance in a test sample that can be reliably distinguished from zero;

• a limit of quantification (LOQ) which is the smallest concentration of the substance in a test sample which can be measured reliably.

The LOD is always lower than the LOQ and both quantities act as censoring values. In the data used in this work, the LOD is unknown: the left censoring phenomenon corresponds therefore to the LOQ of the measuring stations. When both limits are known, one can adapt the model proposed in Section 2.2 to take both of them into account: this would translate into a slightly more complex likelihood as the one derived in Appendix A as we need to consider three cases (when the concentration is between 0 and the LOD, when the concentration is between the LOD and the LOQ, and finally when the concentration is observed and larger than the LOQ). Among the 11,002 recorded measurements during the period of interest, only 12.37% were above the quantification limit. Figure 1 shows the distribution of the number of measurements per station: the mean and median number of samples collected by each monitoring station are respectively 32 and 19. This illustrates that sampling rates are different across stations, most of them making few measures, and the monitoring process is heterogeneous. One may note here that despite the aggregation process, the coarse series remains irregularly sampled, and that for about two thirds of the days included in the studied time span, no measurements were made.

Graphical representation of the station network

The stations network G = (V, E) introduced in Section 2.1 is built using the hydrographic map of the Centre-Val de Loire region. Indeed, once the monitoring stations are geolocalized through their GPS coordinates, one still has to compute the edges between them, as well as the associated weights.

For the data at hand, edges are determined using the river network. A database provided by the French National Institute of Geographic and Forest Information (IGN) (de l'Information Géographique et Forestière, 2021) contains a fine-grained description of rivers, encoded as sequences of hydrographic sections (or river sections). River sections are segments with constant geographic and hydrographic attributes.

The procedure used for computing the edges in the stations network based on the river network may be summarised as follows:

1. One starts by building a river network R = (S, H), where the vertices S are made of the connecting points between the river sections, and the edges H contain all sections. Each edge is thus naturally weighted by the length (in meters) of the corresponding river section.

2. Each monitoring station v i in V is assigned to the closest node si in the river network R, by minimizing the geographical distance between the station v i and all connecting points si = min s∈S d(v i , s).

3. Given two stations v i , v j ∈ V and their associated connecting points si , sj ∈ S, an edge will be generated between v i and v j if there exists at least one path between si and sj . Furthermore, the weight associated to an edge (v i , v j ) is equal to the length of the shortest path between si and sj .

One may notice at this point that the above procedure may result into an unconnected graph, with several connected components.

For illustration, Figure 3 displays a subgraph involving 129 stations of interest that will be used later in Section 5.2. The resulting subgraph is not fully connected and exhibits 5 distinct connected components.

Results

Temporal segmentation

First, the coarse-grained time series D in Figure 2 is segmented using the change-point detection procedure described in Section 3.1. Since the data distribution is right-skewed and shows heavy tails, a Weibull distribution is selected for the parametric distribution Q (see Appendix A for technical details).

Under the assumption above, one should fit Weibull distributions within each stationary temporal segment. In order to limit the number of parameters to estimate, and also to avoid numerical issues rapidly induced by the large number of censored data, the shape parameter σ in the Weibull distributions is supposed not to vary with the change-points. σ is thus constant throughout the series, and is estimated globally under a stationary hypothesis. The only parameter supposed to be varying at each change-point is therefore the rate of the Weibull distribution, say λ. From the application point of view, this simplifying assumption corresponds to the hypothesis that the differences in usage and diffusion of the prosulfocarb among the different users is captured by the shape parameter, and should not vary much over time. On the contrary, the overall average usage of prosulfocarb varies, and this dependency is captured by changes in the rate parameter.

Hence, after computing the MLE of σ, σMLE , over the whole time series, change-points and rate parameters over each temporal segment are estimated by minimizing the cost function in Equation 10. Let us remark here that the estimated value of the shape parameter is σMLE = 0.3. This confirms the data has a heavier tail than an exponential distribution (σ=1), and that the assumption of using Weibull distributions for our data is appropriate.

In the change-point detection procedure, the penalty value for the PELT algorithm was calibrated using a large range of values explored according to the CROPS algorithm.

The range, inspired by the BIC criterion, was set to [ log (K) 5 , 5 log(K)] where K is the number of daily maximum concentrations available, here K = 1, 808. Note that when using the BIC penalty in change point detection, the penalty term written in section 3.1 becomes : β K (L + 1)D = D 2 log(K)(L + 1) = 1 2 log(K)(L + 1). The range chosen allows to screen an interval of penalties containing the BIC penalty.

The penalty calibration procedure resulted in 14 different segmentations, with a number of change-points ranging from 1 to 26. The best segmentation is selected using the elbow method, as illustrated in Figure 11 in Appendix C.3. This amounts to a temporal segmentation with L = 10 change-points, illustrated in Figure 4.

According to Figure 4, the usage of prosulfocarb in Centre-Val de Loire shows different patterns throughout time. Before 2016, most of the values are not quantified, and there are almost no change-points detected. Starting with 2016, two regimes of pesticide usage appear to emerge, and correspond respectively to the periods of intensive usage of prosulfocarb and to the off-peak periods. Indeed, the starting dates of the peak periods coincide with the season where the substance is spread, which is Autumn. The emergence of this two-regime pattern, alternating high concentration values during the peak periods and low concentration values during the off-peaks, is correlated with an important increase in the prosulfocarb sales as shown in Figure 10 in Appendix C.2.

Spatial segmentation

The second step of the analysis consists in the spatial segmentation using the graph-based clustering on the monitoring stations network. Because not all monitoring stations are active during a temporal segment, the clustering procedure described in Section 3.2 is slightly modified: spatial clustering is applied to the subgraph induced by the active stations only.

For illustration, let us focus on a specific temporal segment. An off-peak period, spanning between February 7, 2017 and September 14, 2017 was selected. This period was identified as a homogeneous temporal segment by the change-point detection procedure. This period is highlighted by the black rectangle in Figure 4.

During the selected period, 129 monitoring stations only produced at least one measure. The spatial clustering algorithm was applied with a number of potential clusters varying between 5 and 35 (the minimum number of clusters is equal to the number of connected components in the subgraph induced by the active stations). The optimal number of clusters was selected using the elbow method applied to the inertia curve. According to this criterion, illustrated in Figure 12 in Appendix C.3, the best solution is made of a 10-clusters configuration. The spatial segmentation for the selected off-peak period is illustrated in Figure 5.

To check the relevance of the homogeneity assumption formulated in section 3.2, we computed the within cluster average empirical pairwise Wasserstein distance and observe that for 8 clusters out of 10, this indicator is less than 0.0015, whereas the global average pairwise Wasserstein distance for the 129 stations is 0.003. This suggests that the distance chosen for our station graph is indeed a good proxy of the homogeneity in the concentration space. Additional comments can be made when we look at the geography of the region. Some clusters are overlapping with hydro-ecoregions. Hydro-ecoregions are geographic entities in which hydrographic ecosystems share common characteristics. The criteria defining them combine properties of geology, terrain and climate [START_REF] Wasson | Approche régionale de la typologie des eaux courantes et éléments pour la définition des peuplements de référence d'invertébrés[END_REF]. The borders of those regions are drawn in grey in Figure 5. This ensures that the substances will have homogeneous dispersion properties on these clusters (see clusters 2, 3, 4, 6 and 10 for instance). As expected the biggest component in Figure 3 is the most segmented. Some among them are easy to identify, for instance clusters 9 corresponds to the Indre river. Cluster 1 is identified as the most western part of the Loire and its tributaries mainly the Vienne and the Creuse rivers. Clusters 5,7 and 8 are a little bit harder to identify. If one look closely at the map of the region, there is a high presence of small channels all across this part of the region.

Anomalous cluster identification

Following the methodology proposed in 3.3, the scaling parameter λ k of the aggregated data of each spatial cluster found in Section 5.2 was estimated. The statistics Īk was set to 1/ λk . The Pareto front involving the two descriptors Wk and Īk was computed. It led to the cluster ranking displayed in Figure 6 using the rPref package [START_REF] Roocks | Computing Pareto Frontiers and Database Preferences with the rPref Package[END_REF]. We recall that the selected time segment corresponds to a period of non-use of prosulfocarb. From this it can be deduced that finding quantified measurements of the substance during this period is an anomaly. Three clusters stood out with a Pareto front level of 1. Among them we can find from left to right on Figure 6 : • a cluster which has recorded a very high concentration value but is composed of one station only, it corresponds to cluster 10 in Figure 5. Four measurements were made at that station among which one was quantified. It leads to a higher quantification rate than in the other clusters and high value of 1/ λk .

• a cluster where one can find different profiles of stations. Some of them recorded high concentration values while the others did not provide any quantified measurement. It corresponds to cluster 4 in Figure 5. Nine stations compose this cluster. These stations performed 68 measurements among which five were quantified. The quantified measurements can be found in three different stations.

• a cluster composed of 24 stations that recorded 99 samples. This cluster recorded the maximum level of concentration during the selected time segment. Only six measurements were quantified and they are distributed over three stations. It corresponds to cluster 5 in Figure 5. The low rate of quantification leads to a low value of 1/ λk but the heterogeneity between concentration profiles is higher given that the station that recorded the maximum is very abnormal.

Figure 7 displays the Pareto front levels on the station map.

It is interesting to note that the Pareto front level is not uniformly distributed in the region. The three anomalous clusters are located in the east of the region. It could be related to the agricultural practices and land use. For the sake of the argument, we present in Appendix C.1 the map of barley and wheat crops in Centre-Val de Loire. In future works, we shall investigate the spatial correlation between anomalous clusters and areas with high concentration of these crops.

Conclusion

A new methodology for investigating abnormal signals in pesticide concentration data was introduced in the present manuscript. It takes into account both the left-censored data distribution due to quantification limits, and the spatio-temporal heterogeneity due to measurements made at different stations with irregular frequency. The case study illustrated above yields promising results. Indeed, change-point detection shows clear temporal patterns and identifies periods of intensive prosulfocarb use. Focusing on specific temporal segments and combining those with the spatial clusters allows to characterise concentrations homogeneity, and peak values present in each cluster. Spatio-temporal patterns may be thus highlighted using Pareto front levels, with contiguous clusters, for instance, emerging in the Eastern part of Centre Val de Loire. However, this paper also shows that the non-regular and non-standardized collection of concentration data requires the development of a complex methodology. Therefore, comparison with existing work on spatiotemporal signal analysis is difficult [START_REF] Hamdi | Spatiotemporal data mining: a survey on challenges and open problems[END_REF]. The exploratory aspect of this approach also makes comparison with existing methods difficult. This issue should encourage the various analytical laboratories performing the measurements to standardize station instrumentation and time their sampling rhythm to obtain regular time steps and homogeneous data precision. These first results are encouraging and open up new perspectives for improvements and future developments. This application may allow some operational advances in chemical analysis by a government agency. In this particular case, we indicated that Anses changed the prosulfocarb's approval in October 2017. The time segment we study in this paper precedes this period, covering the period from February to October 2017. In subsequent years, we observe a seasonal signal. For example, we could look at the other off-peak time segments in subsequent years to see if there is still as much prosulfocarb in the periods of non-use. This would help the agency measure the effectiveness of the market approval change.

In Section 5, we observe that the number of samples of each cluster is an important factor that affects the detection of anomalies. In fact, the detection of cluster 10 is not due to the observation of an anomalous or important quantity, but to a strong quantification of the samples. This quantification results from the fact that there are only four measured values among which one only is quantified. Since this result is not particularly satisfactory, there are several paths and ideas to tackle this problem. First, one could imagine changing the nature of the abnormality criteria. For example, one could simply look at the average of the Wasserstein distances of the immediate neighbors in the station graph. This would allow one to include more local information in the analysis of a cluster. Second, one could also change the structure of the station graph. Based on additional information, such as that found in studies of the distribution of pesticides in the air, one could edit the edges of the station graph. For example, one could connect stations that are linked by parcels of the same crop. This would have the effect of grouping components that are not connected in the current graph. The weight of these new edges could be determined by the numerous studies of pesticide dispersal in other environmental compartments (land use, air, groundwater ...etc) during spraying [START_REF] Payraudeau | Modelling pesticides transfer to surface water at the catchment scale: a multi-criteria analysis[END_REF][START_REF] Rozemeijer | The groundwater contribution to surface water contamination in a region with intensive agricultural land use (noord-brabant, the netherlands)[END_REF][START_REF] Tong | Modeling the relationship between land use and surface water quality[END_REF]. Another option to perform the whole procedure is to couple change point detection with clustering and perform joint spatiotemporal segmentation instead of treating temporal and spatial heterogeneity separately. It may also be mentioned that it would be interesting to look at a multivariate approach. Co-occurrence effects between different substances are well studied [START_REF] Baas | Comparison and evaluation of pesticide monitoring programs using a process-based mixture model[END_REF][START_REF] Schreiner | Pesticide mixtures in streams of several european countries and the USA[END_REF]. A new modeling on this subject could bring a better understanding of the field results.

A Change-point detection in left-censored Weibull distributions

We instantiate in this section the general model and methodology presented in the Section 2.3 and Section 3.1 in the particular case of the Weibull distribution.

If Y is a random variable distributed according to a Weibull distribution with a rate parameter λ, a shape parameter σ, the p.d.f. of Y is given by

f Y (y; λ, σ) = I y≥0 σλ (λy) σ-1 e -(λy) σ ,
where I is the indicator function. The c.d.f. of Y is given by

F Y (y; λ, σ) = I y≥0 1 -e -(λy) σ .
Then if Y is subjected to left censoring with a deterministic threshold q, its likelihood function is given by L(λ, σ; y) = σλ (λy) σ-1 e -(λy) σ Iy>q 1 -e -(λq) σ Iy=q .

This likelihood function does not provide closed form estimators for λ and σ but any classical continuous optimisation algorithm can be used to obtain reasonable approximation of those estimators. (as in the fitdistrplus R package (Delignette-Muller & Dutang, 2015) used in the present work).

For change point detection, we use a fixed and global shape parameter and a segment specific rate parameter. More precisely, for a segment [η l-1 + 1, η l ] we use the rate parameter λ l (and the global shape parameter σ). Then under the assumptions of Section 2.3, the negative log likelihood of y η l-1 +1 , . . . , y η l is given by

-L Q (λ l , σ; y η l-1 +1 , . . . , y η l ) = - η l k=η l-1 +1 I y k =q k log 1 -e -(λ l q k ) σ + η l k=η l-1 +1 I y k >q k ((λ l y k ) σ -log(σλ l ) -(σ -1) log(λ l y k )) , (9)
where the q k are the known censoring thresholds.

We denote σ the shape parameter estimated globally. Then the cost function is given by

C σ (L, η, λ; D) = L+1 l=1 -ln L Q (λ l , σ; y η l-1 +1 , . . . , y η l ) + β K (L + 1), (10)
where we use σ as a subscript to emphasize the fact that the shape parameter is not estimated in each segment. This cost function can be analyzed as proposed in Section 3.1 by plugging in the MLE estimator of λ obtained by maximizing L Q (λ l , σ; y η l-1 +1 , . . . , y η l ) for l = 1, . . . , L + 1. Notice that despite the fact σ is fixed during those calculations, we still need to rely on an optimisation algorithm as no closed form formula can be derived to compute the rate parameters under censoring.

B Assessing the homogeneity of the measurements of the stations

The Wasserstein distance was chosen over the Kolmogorov-Smirnov or the Jensen-Shannon metric. It has the advantage of integrating in the distance calculation both the differences between the probabilities of observing different values but also the distances between those values. This is a critical point which is illustrated on a simple simulated example provided by Figure 8. We show here three monitoring stations that have quite different behaviors. Those different behaviors are obvious both on in the temporal representation and in the histograms. However, the Kolmogorov-Smirnov distance between stations 1 and 3 is equal to the Kolmogorov-Smirnov distance between stations 1 and 2. This distance cannot capture the fact that station recorded higher concentration values than station 3. On the contrary, the Wasserstein distance between stations 1 and 3 is smaller than the Wasserstein distance between stations 1 and 2.

Computing information theoretic distances/dissimilarities such as the Jensen-Shannon divergence requires estimating densities for the distributions observed at the stations. As noted earlier, few concentration records (and even fewer quantified ones) are available at the level of a station and within a time period. Therefore, density estimations based on such a small number of observations are unreliable. October 13, 2016;February 7, 2017;September 14, 2017;January 19, 2018;October 5, 2018;January 28, 2019;October 11, 2019;March 25, 2020. The black rectangle corresponds to the selected temporal segment in section 5.2

Figure 1 :

 1 Figure 1: Distribution of the number of measurements per station.

Figure 2 :

 2 Figure 2: Plot of daily maximum concentrations

Figure 4 :

 4 Figure 3: Connected components of 129 stations subgraph

Figure 7 :

 7 Figure 7: Pareto front displayed on the stations map

Figure 8 :Figure 12 :

 812 Figure 8: Example of three stations data. The data were simulated.

See for example the air quality data reported by UK-AIR https://uk-air.defra.gov.uk/.

ANSES stands for Agence nationale de sécurité sanitaire de l'alimentation, de l'environnement et du travail, i.e., French Agency for Food, Environmental and Occupational Health & Safety.

Data exported in September 2020 using http://www.naiades.eaufrance.fr/acces-donnees#/ physicochimie/resultats?debut=09-01-2007&fin=08-09-2020&regions=24&parametres= 1092&fractions=23&supports=3&qualifications=1

C Supplementary figures

C.1 Regional map of crops

The regional map of crops provided in Figure 9 have been produced using data from the registre parcellaire graphique produced by the IGN (de l'Information Géographique et Forestière, 2020).

C.3 Elbow methods

Year

Annual sold weight in the region (kg)

Number of breaks Cost of segmentation

Figure 11: Elbow method illustration for the temporal segmentation. The optimal penalty value according to this approach is 8.36.