Yukiko Kenmochi

Lilian Buzer

Akihiro Sugimoto

Ikuko Shimizu

Planar surface segmentation of a point cloud using local geometric patterns

This paper presents a method for segmenting a 3D point cloud into planar surfaces using recently obtained discretegeometry results. In discrete geometry, a discrete plane is defined as a set of grid points lying between two parallel planes with a small distance, called thickness. Contrarily to the continuous case, there exist a finite number of local geometric patterns (LGPs) appearing on discrete planes. Moreover, such a LGP does not possess the unique normal vector but a set of normal vectors. By using those LGP properties, we first reject non-linear points from a point cloud, and then classify non-rejected points whose LGPs can have common normal vectors into a planar-surfacepoint set. From each planar-surface-point set, we also estimate parameters of a discrete plane by minimizing its thickness.

Introduction

This paper presents a method for segmenting a 3D point cloud into planar surfaces using recently obtained discretegeometry results. Conventional approaches are classified into three categories: region-based approach, edgebased approach and hybrid approach. The first one merges points that have similar region properties calculated from their neighboring points such as normal vectors and curvatures [1]. Because calculated properties are sensitive to noise and quantization errors, it is known that they cause over segmentation. In the second approach, edges are searched for such that they separate regions by using depth discontinuities [START_REF] Zhao | Range-data-based object surface segmentation via edges and critical points[END_REF]. Because edges are not always extracted as connected curves, they cause under segmentation, contrary to the first ones. The third approach is hybrid between the two approaches, so that they import both their merits [START_REF] Stamos | 3d model construction using range and image data[END_REF][START_REF] Yokoya | Range image segmentation based on differential geometry: a hybrid approach[END_REF]. For a special case in the third approach, a planar segmentation method is proposed based on locally planar points, considering points that are not locally planar to be potentially edge points [START_REF] Stamos | 3d model construction using range and image data[END_REF].

In this paper, we present a discrete version of the hybrid method. In discrete geometry, a discrete plane is defined as a set of grid points lying between two parallel planes with a small distance, called a thickness [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF]. Contrarily to the continuous case, there exist a finite number of local geometric patterns (LGPs) appearing on discrete planes, called linear LGP [START_REF] Kenmochi | Digital planar surface detection using local geometric patterns[END_REF]. In fact, points that have linear LGP can be considered to be discrete version of locally planar points [START_REF] Stamos | 3d model construction using range and image data[END_REF].

In addition, such a linear LGP does not possess the unique normal vector but a set of normal vectors [START_REF] Coeurjolly | On digital plane preimage structure[END_REF]. By using those LGP properties, we present a segmentation method following the two steps: first reject non-linear points from a point cloud (edge-based part), and then merge non-rejected points whose LGPs have common normal vectors (regionbased part). It thus uses only precalculated look-up tables with respect to LGP, and does not require any parameter setting. Furthermore, our method is less sensitive to noise as well as quantization errors. Indeed linear LGPs already take into account quantization errors for their generation. We show such advantages by applying our algorithm to 3D point clouds such as range images. In order to evaluate our segmentation results, we also present a method for estimating discrete plane parameters from each segmented planar surface by minimizing its thickness. Because the thickness indicates the segmentation inaccuracy, namely the curvedness of a segmented planar surface, we consider that the thinner the thickness, the better the segmentation result.

Non-linear point rejection using LGP

Discrete planes

Let R be the set of real numbers. A plane P in the 3D Euclidean space R 3 is defined by the following expression:

P = {(p, q, r) ∈ R 3 : αp + βq + γr + δ = 0}
where α, β, γ, δ ∈ R. Let Z 3 be the set of grid points whose coordinates are integers in R 3 . A discrete plane, which is a digitization of P, is then defined such that

D(P) = {(p, q, r) ∈ Z 3 : 0 ≤ αp+βq +γr +δ < ω} (1)
where ω = max (|α|, |β|, |γ|), called the thickness [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF].

Linear LGP on discrete planes

We consider a cubical grid-point set Q(x) whose edge length is 2 around a point x ∈ Z 3 such that

Q(x) = {y ∈ Z 3 : x -y ∞ ≤ 1}.
Let us assume that each point in Z 3 has a binary value such as either 1 or 0. Such a pattern of binary points in Q(x) is called local geometric patterns, abbreviated to LGP. There are 2 26 different LGP for Q(x) providing that the central point x always has the fixed value 1. This indicates that x is considered not to be a background point but to be an object point.

Among those different LGPs, we investigated which LGP can appear on discrete planes [START_REF] Kenmochi | Digital planar surface detection using local geometric patterns[END_REF]. This problem is mathematically written as follows. Let F be a set of points whose binary values are 1 in Q(x). If there is a plane P such that

F = D(P) ∩ Q(x) = {(p, q, r) ∈ Q(x) : 0 ≤ αp + βq + γr + δ < ω}, (2)
we say that F forms a discrete plane in Q(x). Therefore, our problem is solved by looking for all possible F, namely

LGP, satisfying [START_REF] Coeurjolly | On digital plane preimage structure[END_REF]. Such LGP are called linear LGP. Since this problem is considered to be the feasibility of the inequalities of (2) for all (p, q, r) ∈ F, we need to check if there are feasible solutions α, β, γ, δ for each different LGP of Q(x). If they exist, such LGP can appear on discrete planes and become linear LGP. However, [START_REF] Kenmochi | Digital planar surface detection using local geometric patterns[END_REF] shows that we can avoid computing the feasibility test for all 2 26 LGPs of Q(x), by taking an approach based on arithmetic planes [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF], which are related to discrete planes. An algorithm is then proposed to generate all linear LGPs, and it is found that there exist only 34 LGPs that appear on discrete planes, called linear LGPs, up to translations, rotations and symmetries, as shown in Fig. 1. Note that they are generated with the constraints

0 ≤ α ≤ β ≤ 1, γ = 1. (3)
In order to visualize the shapes of linear LGPs in Fig. 1, we add polyhedral meshes generated by applying a discretemarching-cube-like method for the 18-neighborhood system [START_REF] Kenmochi | Combinatorial boundary of a 3d lattice point set[END_REF] to a digitized half space. Interior points of objects are designated as black points in the figures.

Locally linear and non-linear points

Experimentally, those linear LGPs can be seen not only on discrete planes but also on discrete smooth surfaces. Intuitively, this is not difficult to understand, since any local surface patch on a smooth surface can be approximated to a planar surface when the size of the patch becomes small. In the discrete space, even if a point has a linear LGP, we are uncertain whether such a point appears on a planar surface or a non-planar surface. Contrarily, if a point has a non-linear LGP, it is exactly a point that never appears on a planar surface. From this reason, if a point has a linear LGP, it is called a locally linear point, otherwise, simply called a non-linear point.

Non-linear point rejection

By simply checking the LGP linearity, we can therefore reject non-linear points from a grid-point set, since the linear LGPs play an important role in filtering linear points. Note that it is realized by looking up the binary table of LGPs (linear or not). In actual experiments, we see that isolated points that are considered to be noise are automatically rejected as well as points around surface edges. However, it is also observed that some points around an edge are not rejected, especially when two adjacent surfaces joining the edge intersect at an obtuse angle. This fact implies that a simple post-processing, such as the connected component labeling [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF] of a non-rejected point set, does not always give satisfactory results for planar surface segmentation.

Planar surface segmentation of locally linear points

In order to solve the above problem, we propose a method using not only the point connectedness but normal vectors derived from LGPs.

Feasible normal vectors of linear LGPs

A linear LGP is a discrete plane patch of D(P) in a bounded space Q(x), denoted by D Q(x) (P). Given a D Q(x) (P), we can find a set of Euclidean planes P such that the digitization of each of those planes in Q(x) is equal to D Q(x) (P). The set of all such Euclidean planes is called the preimage and it is known that the correspondence between discrete plane patches and Euclidean planes is not one-to-one but one-to-many [START_REF] Coeurjolly | On digital plane preimage structure[END_REF]. Because of the one-tomany correspondence, the preimage of D Q(x) (P) is represented by a set of parameters α, β, γ, δ. More precisely, the preimage is obtained as a feasible solution set of the inequality set of (1) for all points of D Q(x) (P). It means that the preimage is given by a convex polytope in the parameter space [START_REF] Coeurjolly | On digital plane preimage structure[END_REF].

Because all interesting parameters in this paper are translation-invariant, we focus on the three parameters α, β, γ indicating the normal vector of P, distinguished from the intercept δ of P. We thus apply the Fourier-Motzkin elimination to the inequality set of (1) for all (p, q, r) ∈ D Q(x) (P), so that a set of feasible normal vec-tors is calculated from each linear LGP. Remark that all calculations are done by using only integers, i.e., they cause no rounding errors; the details are found in [START_REF] Kenmochi | Digital planar surface detection using local geometric patterns[END_REF]. The results are derived in the space (α, β) from linear LGP with the constraints (3). The feasible region for each linear LGP is obtained as a convex polygon in the triangle region whose vertices are (0, 0), (0, 1) and (1, 1) of the space (α, β) because of (3). Figure 2 illustrates that the inequality set of (1) for all (p, q, r) ∈ D Q(x) (P) divides the triangular region in the space (α, β) into triangular or quadrilateral polygons, called normal cells. The feasible region of each linear LGP is given as a set of normal cells that constitutes a convex polygon in the space (α, β). Table 1 shows such a set of normal cells whose union corresponds to a convex polygon representing a set of feasible normal vectors for each linear LGP depicted in Fig. 1.

Discrete Gaussian sphere

The 26 normal cells in Fig. 2 are generated with the constraints (3). We embed these normal cells into the 3D space (α, β, γ) with γ = 1, as illustrated in Fig. 3 (left). The triangle surrounded by thick lines in Fig. 3 (left) corresponds to the triangular region that is the union of normal cells in Fig. 2. Once the normal cells are embedded into the space (α, β, γ), we make the congruous ones by applying to them 48 transformations of rotations and symmetries of a cube of edge length 2, centered at the origin of the 3D space. We see, in Fig. 3 (left), that there are the 48 triangles on the cube, so that the whole cube contains 1248 normal cells. Such a cube is called the cubical Gaussian sphere.

We now project normal cells tiled on the cubical Gaussian sphere onto a unit sphere centered at the origin, as illustrated in Fig. 3 (right). The unit sphere separated by normal cells is called the discrete Gaussian sphere, because the size of normal cells indicates the resolution of digitized normal vectors calculated from linear LGPs. The triangle surrounded by red lines in Fig. 3 (right) corresponds to the triangle surrounded by thick lines in Fig. 3 (left) that corresponds to the union of normal cells in Fig. 2. In the remainder, we denote G the set of all normal cells on the discrete Gaussian sphere. Remark that we use only integer or rational numbers to calculate all normal cells, which are related to a cubical Gaussian sphere.

Unified discrete Gaussian image

By using the discrete Gaussian sphere, we give a discrete version of extended Gaussian images that are useful for representing surface shapes [START_REF] Horn | Extended gaussian images[END_REF], called unified discrete Gaussian images. Let us first consider a discrete version of the Gaussian image that is the mapping from an object surface point to its normal vector on the Gaussian sphere. Let V be a locally linear point set in Z 3 . For a point x ∈ V, we define a discrete Gaussian image I(x) as the set of normal cells corresponding to the linear LGP of x. Choosing a normal cell c ∈ G, we now consider a point subset of V such that We then obtain the number of points in R(c) for every c ∈ G, called the unified discrete Gaussian image, such that

R(c) = {x ∈ V : c ∈ I(x)}. (4)
u(c) = |R(c)|. (5)
Note that u(c) and R(c) are generated by simply looking up a table such as Table 1.

The concept of unified discrete Gaussian images is similar to that of extended Gaussian images [START_REF] Horn | Extended gaussian images[END_REF]. The differences from extended Gaussian images are the followings: the function ( 5) is defined with respect to a normal cell c on the discrete Gaussian sphere G, instead of a point n on the Gaussian sphere; the value of ( 5) is the number of grid points x such that I(x) includes c, instead of the area of the surface whose normal vector is n. From the definition, we see that our unified discrete Gaussian image represents a distribution of normal cells of a digital object surface. Figure 4 shows an example of the unified discrete Gaussian images for a digitized box. Concerning cell colors on the discrete Gaussian sphere in Fig. 4 (bottom), the darker the blue cell, the larger the value of u(c), and the red cell has the maximum value. The length of the pale blue needle for each cell c also corresponds to the value of u(c). On a digitized box in Fig. 4 (top), red and blue points are locally linear, while green points are non-linear. Note that red points correspond to the red cell in Fig. 4 (bottom). Figure 4 shows that we can extract a set of grid points that belong to a digital plane D(P) by choosing a "right" normal cell, for example, a red one. This is based on the following fact; if (α, β, γ) is a normal vector of D(P), (α, β, γ) is included in the common normal cell(s) of I(x) for all x ∈ D(P).

Algorithm

By using the unified discrete Gaussian image u(c) and the point sets R(c), we present our algorithm for planar surface segmentation from a locally linear point set V. Our problem is formulated as follows; each point x ∈ V is assigned into one of sets S i for i = 1, 2, . . . such that the points in each S i constitutes a connected planar-surface set. From the previous discussions, our method is founded on the following hypothesis: if there is a connected point subset S ⊆ V such that they have a common normal cell for all x ∈ S, S may constitute a discrete plane.

Based on this hypothesis, we present Algorithm 1. we look for the largest connected grid-point set S i , whose points having a common normal cell by using u(c) and R(c). Because each point has several normal cells, our method cannot be processed in parallel with respect to normal cells. It must be a repeated procedure; once we obtain S i , we remove all points of S i from every R(c), modify u(c), and repeat this procedure after the increment of i. Practically, we would like to avoid obtaining a very small surface patch, so that we set a parameter s that is the minimum size for S i .

Algorithm 1 is thus a loop procedure of seeking planar surfaces S i . Each S i is a maximally connected point set, whose points have a common normal cell. Once we find S i , we check the size of S i in Step 11, and if |S i | ≥ s, we remove all points of S i from every R(c) and also modify u(c) in Step 13. After such modification and incrementing i, we seek a new S i . For finding each S i , we look for the maximum connected component C of each R(c), and then set S i to be the maximum among all C. In order to reduce the frequency of calculation of connected components, which is a global operation, we make a priority queue D k of normal cells with u(c) in Step 4. We then repeat dequeue of a normal cell h from D k to obtain the maximum connected component C of R(h) in Step 8. Comparing the size of C with the maximum among those of other normal cells that are already dequeued from D k , we finally obtain the currently maximum point set S l in Step 9. Note that this loop is repeated until the size of R(h) is less than s or more than the size of S l as described in Step 7. For calculating the maximum connected component of R(h), we apply a simple method based on a depth-first strategy by using a queue [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF]. The time complexity is linear with respect to the size of R(h).

Experimental results

We show results of planar surface segmentation from six range images of the same blocks, which are taken from two different viewpoints with three different resolutions. The results are illustrated in Figs. 5 and6. In the cases of Fig. 5, the numbers of valid (measured) points are 12858 for (a), 51740 for (b) and 207448 for (c). Among those valid points, we have 11343 locally linear points for (a), 47034 for (b), and 185566 for (c), respectively. Similarly, in the cases of Fig. 6, the numbers of valid (measured) points are 12142 for (d), 48802 for (e) and 195765 for (f). Among those valid points, we have 10663 locally linear points for (d), 43981 for (e), and 176619 for (f), respectively. Table 2 show the number of locally linear points that are assigned to each segmented planar surface, and their corresponding color in Figs. 5 and6. We see that 10, 9 and 12 planar surfaces are found in Fig. 5 We see in Figs. 5 and6 that non-linear points, colored light green, appear around edges of block faces, and sometimes appear in faces because of noise in the range images. Because we set the minimum surface size s, there are locally linear points that construct no planar surface whose size is not less than s around the points, colored black in the figures. Locally linear points become black, if they are obtained too sparsely in the 3D space to be connected, so that they are seen as some of block faces in Figs. 5 and6. Otherwise, they are considered to have bumps on the whole, even though they are locally linear. From Figs. 5 and6, we can conclude that most of the block surfaces are segmented by our simple algorithm, which require neither complicated parameter setting nor parameter estimation.

Estimation of discrete plane parameters

Formulation

From each segmented planar-surface set S i , we estimate its discrete-plane parameters as follows. In order to simplify our problem, we first consider the case that ω = γ. From (1), we obtain a linear inequality set such that, for all (x, y, z) ∈ S i ,

0 ≤ α ′ x + β ′ y + z + δ ′ ≤ ǫ (6)
where

α ′ = α ω , β ′ = β ω , δ ′ = δ ω , ǫ ≥ 0. A solution set (α ′ , β ′ , δ ′
) is then obtained by minimizing ǫ under the above constraints. In this framework, if ǫ < 1, S i is recognized as a discrete plane patch exactly; otherwise, S i is recognized as a set of grid-points between two parallel planes whose distance is wider than the thickness of a discrete plane. Geometrically, our method looks for two parallel planes such that the distance between them becomes minimum.

For all the other cases such that ω = -γ, β, -β, α, -α, we simply need to modify [START_REF] Horn | Extended gaussian images[END_REF], so that the following inequalities are obtained respectively

0 ≤ -α ′ x -β ′ y -z -δ ′ ≤ ǫ,
0 ≤ α ′ x + y + γ ′ z + δ ′ ≤ ǫ 0 ≤ -α ′ x -y -γ ′ z -δ ′ ≤ ǫ, 0 ≤ x + β ′ y + γ ′ z + δ ′ ≤ ǫ 0 ≤ x -β ′ y -γ ′ z -δ ′ ≤ ǫ,
where γ ′ = γ ω . Practically, we simultaneously use the above 6 types of inequalities to find a parameter set minimizing ǫ.

results

Tables 3 and4 show the estimation results for segmented planar surfaces obtained in the previous section, as illustrated in Figs. 5 and6. We see that the parameter values of α, β and γ that are obtained for the corresponding planar surfaces, segmented from the range images with different resolutions, are very similar. Concerning to the parameter δ, the values in Note that we always set the grid interval to be 1, when we make a grid space from a range image. From Tables 3 and4, we also see that it is rare that ǫ becomes less than 1. In other words, most of our segmented planar surfaces cannot be exactly discrete planes. Moreover, the tables show that the higher the image resolution, the larger the value ǫ. Since each segmented planar surface contains many grid points when the image resolution is high, as seen in Table 2, it can generate a very thicker discrete plane. It is possible that the thickness is related to the spread of grid points in a segmented planar surface as well as the point number. It might be interesting to study how we can reduce the thickness ǫ by changing the image resolution, with the aim of inventing a multiscale method for range image registration by using planar surfaces, for example.

Conclusion

In this paper, we present a discrete version of the hybrid method for planar surface segmentation from a 3D gridpoint set. Our method simply requires two types of look-up tables, such as the binary LGP table (linear or non-linear) and the normal cell list with respect to each linear LGP, and does not require any parameter setting/estimation. The experimental results in Figs. 5 and6 show us that our method is useful for planar surface segmentation from a point cloud, because it takes into account not only quantization errors but also noise. We also present a method for estimating discrete-plane parameters, which is also based on discrete geometry. Our estimation results in Tables 3 and4 show us that exact discrete planes are rarely obtained for practical images, because of their noise. Theoretically, exact discrete planes must be obtained if input is an ideal image, i.e., it does not contain noise, but contains only quantization errors. However, practically, it is no wonder that input image contains noise as well as quantization errors. Therefore, we have to derive a conclusion that we need to eliminate such noise, for example, by reducing image resolutions, before applying our method based on discrete geometry. Because our method is fully discrete and such discreteness helps us to build up a multiscale approach, we will reorient our future work to inventing the multiscale method for range image registration by using discrete planes, for example. We expect that our approach will provide a rough registration result with less computation.

Figure 1 .

 1 Figure 1. The 34 linear LGPs.

Figure 2 .

 2 Figure 2. Normal cells on the αβ-plane with constraint (3).

Figure 3 .

 3 Figure 3. The cubical Gaussian sphere and the discrete Gaussian sphere.

Figure 4 .

 4 Figure 4. A synthetic 3D image of a box (top) and its unified discrete Gaussian image (bottom).

Algorithm 1 :

 1 Planar surface segmentation input : a unified discrete Gaussian image u(c), point sets R(c), and a minimum surface size s output: planar-surface point sets S i for i = 1, 2, 3, . . . begin 1 initialize a label such that l = 0; D k of normal cells with 4 priorities of values u(c); increment l and initialize S l = ∅; 5 set h to be the highest priority cell in D k and 6 remove it from D k ; while |R(h)| > max(s -1, |S l |) do 7 set C to be the maximum connected 8 component of R(h); if |C| > |S l | then set S l = C; 9 reset h to be the highest priority normal 10 cell in D k and remove it from D k ; if |S l | ≥ s then 11 forall c such that u(c) = 0 and 12 R(c) ∩ S l = ∅ do reset R(c) = R(c) \ S l and 13 u(c) = |R(c)|; until |S l | < s; 14 return S i for i = 1, 2, . . . , l -1;

Figure 5 .Figure 6 .

 56 Figure 5. Planar surface segmentation results from range images of blocks, which are taken from the same viewpoint, with different resolutions: the image sizes are 160 × 120 (a), 320 × 240 (b), 640 × 480 (c). The minimum surfaces sizes s are set to be 100 (a), 500 (b), and 1000 (c), respectively.

Fig. 5 (

 5 b) and Fig.6(e) (resp. Fig.5 (c) and Fig.6 (f)) is twice (resp. four times) as large as that of

Fig. 5 (

 5 Fig.5(a) and Fig.6 (d), because of their image resolutions. Note that we always set the grid interval to be 1, when we make a grid space from a range image.From Tables3 and 4, we also see that it is rare that ǫ becomes less than 1. In other words, most of our segmented planar surfaces cannot be exactly discrete planes. Moreover, the tables show that the higher the image resolution, the larger the value ǫ. Since each segmented planar surface contains many grid points when the image resolution is high, as seen in Table2, it can generate a very thicker discrete plane. It is possible that the thickness is related to the spread of grid points in a segmented planar surface as well as the point number. It might be interesting to study how we can reduce the thickness ǫ by changing the image resolution, with the aim of inventing a multiscale method for range image registration by using planar surfaces, for example.

Table 1 . Linear LGPs and their normal cells.

 1

	linear LGP	normal cells
	1	0 25
	2	1 9 11 12
	3	4 5 7 10 23
	4,5	0 1 16 17 18 24
	6,17	2 3 4 5 7 8
	7	2 3 5 8
	8,9	6 9 10 11 14 15 21 23
	10,12	8 19 20 25
	11	8 17 18 19 20
	13,28	2 3 4 5 6 7 9 10 11 12 13 14 15 21 22 23
	14	2 3 6 13 14 15 16 21 22 24
	15	2 3 6 11 12 13 14 22
	16	4 5 7 10 23
	18,19	0 18 19 25
	20,23	0 1 3 8 12 13 16 17 18 19 20 22 24 25
	21,22	3 8 16 17 20 22
	24,25	1 9 11 12 13 14 15 24
	26,34	2 4 5 6 7 10 21 23
	27	2 5 6 7 21 23
	29,30	0 17 18 19 20 25
	31,32	1 12 13 16 22 24
	33	6 9 11 14 15 21

Table 2 . Colors and point numbers of seg- mented planar surfaces in Figs. 5 and 6.

 2

	color	(a)	(b)	(c)	(d)	(e)	(f)
	1 blue	1750	6529	23393	1914	6304	15823
	2 yellow	1557	6125	19822	1280	4254	13809
	3 pink	1536	5573	16645	1171	4006	12879
	4 pale blue	1180	4408	14377	822	3505	12198
	5 orange	595	2627	9375	800	2929	9098
	6 green	579	1901	3022	638	2437	8867
	7 brown	418	1620	2376	612	2321	5009
	8 turquoise	361	1424	1609	220	844	2266
	9 olive	197	799	1419			2169
	10 purple	112		1418			1857
	11 violet			1160			1212
	12 moss green			1036

Table 3 . Parameter estimation results of seg- mented planar surfaces in Fig. 5

 3

				(a)
		ǫ	α	β	γ	δ
	1	2.30846	0.0895522	0.512438	1	-480.756
	2	2.51369	0.444688	1	0.779847	345.085
	3	3.47581	-0.475806	-1.16935	1	-480.306
	4	2.01022	-0.434227	-0.157088	1	-477.374
	5	1.26786	-0.482143	1	-0.946429	-508.393
	6	1.47327	1	0.0178218	-0.558416	-335.198
	7	1.125	0.351563	1	0.664062	274.57
	8	0.605714	1	-0.0114286	-0.554286	-334.337
	9	0.5	1	0	-0.5	-305.5
	10	1.51563	-0.78125	-1.32812	1	-468.875
				(b)
		ǫ	α	β	γ	δ
	1	3.90025	0.0910617	0.51279	1	-962.483
	2	4.89735	0.4446018	1	0.780531	691.458
	3	6.31767	-0.474464	-1.16547	1	-963.488
	4	3.24647	-0.433281	-0.153846	1	-954.964
	5	1.54492	1	0.0212766	-0.554374	-668.843
	6	1.72454	-0.478588	1	-0.956019	-1027.81
	7	2.23761	0.356076	1	0.671419	555.339
	8	0.791798	1	-0.0126183	-0.536278	-651.688
	9	0.717703	1	-0.00956938	-0.564593	-680.029
				(c)
		ǫ	α	β	γ	δ
	1	7.38188	0.086763	0.507398	1	-1926.6
	2	9.66667	0.444444	1	0.777778	1378.33
	3	6.23908	-0.422621	-0.149001	1	-1910.67
	4	10.8503	-0.487362	-1.15165	1	-1924.55
	5	2.1997	1	0.0222728	-0.558241	-1347.4
	6	2.44488	-0.479921	1	-0.955906	-2056.17
	7	2.26689	0.375	1	0.726351	1217.36
	8	0.886656	1	-0.0103611	-0.562951	-1358.17
	9	0.875	1	-0.00961538	-0.572115	-1376.38
	10	0.808511	1	-0.0141844	-0.529255	-1290.44
	11	4.24536	1.81477	0.0910698	1	-1930.37
	12	2.17989	0.39418	1	0.626984	1009.64

Table 3

 3 (b) and Table 4 (e) (resp.Table 3 (c) and Table 4 (f)) are almost twice (resp. four times) as large as those in Table 3 (a) and Table 4 (c). The reason is that the grid space of

Table 4 . Parameter estimation results of seg- mented planar surfaces in Fig. 6

 4

				(d)
		ǫ	α	β	γ	δ
	1	2.53866	-0.481622	-0.539924	1	-505.439
	2	2.78801	-0.235546	1	0.940043	446.441
	3	2.02985	-0.41791	0.238806	1	-489.358
	4	1.02676	1	0.0126761	-0.550235	-331.395
	5	1.92059	-0.505882	-0.638235	1	-464.826
	6	0.760252	1	-0.0115668	-0.539432	-326.766
	7	1.74775	-0.0720721	-1.63964	1	-500.802
	8	0.660584	-0.0720721	1	-0.448905	-245.204
				(e)
		ǫ	α	β	γ	δ
	1	4.0725	-0.486857	-0.541143	1	-1010.1
	2	5.03725	-0.232787	1	0.936131	890.082
	3	3.22654	-0.417907	0.232147	1	-978.429
	4	1.75286	1	0.0111605	-0.549891	-663.497
	5	3.50924	-0.50308	-0.652977	1	-927.704
	6	3.09351	-0.073294	-1.64027	1	-1000.16
	7	1.07945	1	-0.00797011	-0.546077	-660.958
	8	1.20741	-0.237037	1	-0.451852	-493.904
				(f)
		ǫ	α	β	γ	δ
	1	4.88192	-0.462209	-0.519922	1	-2020.35
	2	6.38296	-0.413744	0.231031	1	-1956.85
	3	8.11213	-0.227372	1	0.948015	1804.6
	4	2.76525	1	0.0120939	-0.553442	-1335.69
	5	6.59621	-0.51084	-0.651762	1	-1854.09
	6	5.85422	-0.0745445	-1.64605	1	-1998.81
	7	4.33574	-0.520009	-0.586203	1	-2016.94
	8	0.981103	1	-0.0106891	-0.559076	-1350.38
	9	0.888889	1	-0.0121382	-0.531279	-1294.01
	10	1.06774	1	-0.00811321	-0.55717	-1346.02
	11	0.921053	1	-0.0115132	-0.541118	-1314.11


Ali Shakil

Mohammad Ali Khalighi

Pierre Pudlo

Cyril Leclerc

Dominique Laplace

Alexandre Boudonne

Outlier Detection in Non-stationary Time Series Applied to Sewer Network Monitoring

Keywords: Sewer network monitoring, Smart city, Data cleaning, Outlier detection, Internet-of-Things

We consider the case of data processing for a sewer infrastructure where water drains are equipped with wastelevel sensors, which frequently send the related data to a data processing unit. In order to understand the dynamics of waste accumulation within the whole drain network, the collected data should first be pre-processed by removing the unreliable (or, in other words, noisy) measurements. As we show, the evolution of the waste inside a drain can be modeled by a non-stationary discontinuous time series model. Due to the chaotic aspect of the waste and the hostile conditions under which the sensor should operate, the observed time series can include outliers in the form of peaks, which should be removed from the raw data prior to any data processing. This paper proposes an efficient data cleaning algorithm that makes a good compromise between computational complexity and performance. This latter is evaluated in terms of the probabilities of peak detection (i.e., detecting actual outliers) and false detection (i.e., incorrectly denoting measurements as outliers). A trade-off between these two criteria should be made by setting appropriately the detection threshold (which, in the proposed method, does not depend on the mean or variance of the data). For instance, for a threshold of 2.5, the algorithm provides a correct outlier detection probability of 0.85 and a false detection probability of 2.5×10 -2 . The efficiency of the proposed algorithm is demonstrated by applying it to real measurement data.

I. INTRODUCTION

In a more than ever connected world, the Internet-of-Things (IoT) paradigm offers promising capabilities for efficient resource management and environment protection at large scale. In the smart-city context, IoT allows an efficient use of energy resources such as electricity and water through the deployment of dedicated sensors and connectivity solutions. Combining IoT and artificial intelligence (AI) opens doors for the deployment of efficient smart systems, capable of adapting to the actual environmental/field conditions, ultimately enabling optimal management of the resources [START_REF] Yao | Deep learning for the Internet of Things[END_REF], [START_REF] Helal | Challenges and opportunities for data science and machine learning in IoT systems a timely debate: Part 2[END_REF]. Typical examples are smart grids, smart water resource management and quality monitoring, smart sewerage network, etc. [START_REF] Hossain | Application of big data and machine learning in smart grid, and associated security concerns: A review[END_REF], [START_REF] Junior | Internet of Water Things: A remote raw water monitoring and control system[END_REF], [START_REF] Hassan | Internet of Things A to Z: Technologies and Applications[END_REF]. In particular, smart management of potable water distribution and sewerage systems can optimize resource consumption and reduce the impact on the environment [START_REF] Nasser | A two-layer water demand prediction system in urban areas based on micro-services and LSTM neural networks[END_REF], [START_REF] Mustafa | Applications of IoT and artificial intelligence in water quality monitoring and prediction: A review[END_REF]. Indeed, the ecological impact of an inefficient sewerage system can be devastating, appealing for emergency deployment of adequate measures by the suppliers.

This work focuses on the case of a sewerage system for a large city and considers the deployment of an IoT-enabled smart network for managing the water drains as part of the sewer infrastructure. These drains collect rain and water but are also subject to urban wastes (they are not systematically equipped with bars to prevent waste accumulation inside them). The waste can enter the sewerage system, potentially causing important damages, e.g., obstructing the drain evacuation, resulting in local flooding, or entering the pipes of the sewerage system, damaging then the downstream equipment. Also, too much waste accumulation in the drains results in overflowing in the street. In particular, in some locations in Marseille where the drains of the stormwater network are directly connected to the sea, the rain can directly flush the waste into to the sea throughout the pipes, causing hence disastrous damage to the environment. As such, sewer networks are critical and complex systems that need significant maintenance and surveillance, thus the rapidly increasing deployment of IoT-based solutions over the past few years to ensure quality of life and protection of the environment. Examples include automatic or semi-automatic condition assessment to identify damages in pipes [START_REF] Rayhana | Automated vision systems for condition assessment of sewer and water pipelines[END_REF], water quality monitoring [START_REF] Olatinwo | Enabling communication networks for water quality monitoring applications: A survey[END_REF], detection and prediction of sewer overflow [START_REF] Zhang | Sewer system control using artificial intelligence, hydraulic model, and Internet of Things[END_REF], or measurement of hazardous gases in sewer networks [START_REF] Kasat | Smart city solutions on drainage, unused well and garbage alerting system for human safety[END_REF].

In this paper, we focus on the waste level measurement for the storm-water network, and more specifically, on the processing of the associated collected data for the whole sewerage system. In fact, analyzing these data allows to understand the dynamics of waste accumulation inside the drains, and the appropriate actions that should be taken to ensure proper operation of the whole network. Although this application seems to be rather simple at a first glance, the related data processing can be challenging in practice. Indeed, due to the chaotic aspect of the waste and the hostile conditions in which the sensor is installed, the observed time series can present erroneous peaks. The detection and removal of outliers from the raw data collected from a wireless sensor network is an important preliminary step prior to data analysis [START_REF] Al-Qurabat | An overview of periodic wireless sensor networks to the internet of things[END_REF], [START_REF] Finjan | Important extrema points extraction-based data aggregation approach for elongating the WSN lifetime[END_REF], e.g., in view of developing a reliable mathematical model or training an AI-based algorithm [START_REF] Cook | Anomaly detection for IoT timeseries data: A survey[END_REF].

For this purpose, we propose in this paper an efficient outlier detection method with a low computational complexity. The proposed method allows quick cleaning of the data without any need to labeling them, as well as without any modification of the other data points. This latter is critically important in order not to bias the data, which will otherwise result in information loss, impacting hence the mathematical modelling and the event prediction of the network. The originality of this work is related to the nature of collected data in the form of time series, with the properties of being non-stationary and discontinuous, to which most of the previously-proposed solutions in the literature cannot be applied, as explained in detail later in Section 3.

The remainder of the paper is organized as follows. Section II explains in more detail the considered sewerage system, the sensors, and the nature of the collected data. Next, Section III describes the related work on outlier detection in time series. Then, Sections IV and V present our proposed peak detection algorithm and its performance evaluation, respectively. Concluding remarks and discussions are provided in Section VI.

II. CONSIDERED SEWERAGE SYSTEM

A. Context

Towards the deployment of IoT solutions for smart-city applications [START_REF] Mohanty | Everything you wanted to know about smart cities: The Internet of Things is the backbone[END_REF], [START_REF] Pradhan | IoT for smart city: Improvising smart environment[END_REF], one big project of the French company SUEZ has been the digitalization of the water cycle in the city of Marseille, the second largest city of France with over 1.8 million inhabitants in the metropolis. This includes the instrumentation of the city's stormwater network and, more precisely, the corresponding water drains by means of a wireless sensor network (WSN) to monitor the volume of accumulated waste. The objective is to get a deeper understanding of the dynamics of the network and to intervene in a smart manner (i.e., when and where necessary) for waste removal before it causes environmental, health, olfactory, or visual pollution. In other words, the collected data from the sensors are used to establish mathematical models to predict the waste accumulation throughout the Marseille metropolitan area.

B. The Stormwater Network

Given the size of the metropolitan area and its particularities, the associated water and sewer networks are massive and quite complex. The sewerage system can consist of three types of sub-networks: the wastewater and the rainwater networks, and the combined wastewater/rainwater sewer network (the old infrastructure), see Fig. 1. Here, we focus on the abovementioned second and third networks, which form together the so-called stormwater network.

The role of a stormwater network is to absorb the rainwater and to avoid flooding. To do so, the network includes over 18,000 drains that need to be maintained meticulously, which includes cleaning the waste that accumulates therein over time. The possibility of large quantities of generated waste, the elaborate network, and the hilly topography of the city are factors among others that make this task particularly complex. Today, the maintenance system is inefficient and nonoptimized; the whole process necessitates human intervention and supervision, which includes visual assessment of the waste level inside the drains. The maintenance is currently ensured by the cleaning staff who visit each drain regularly and collect the waste if needed. The installation of the above-mentioned WSN allows improving the efficiency of network maintenance and avoiding unnecessary human intervention. So far, more than 3,500 sensors have been installed and the goal is to reach around 5,000 connected sensors by the end of 2022.

The current network covers about half of the city, i.e., an area of around 128 km 2 .

C. Sensor Network

In the considered network, the waste level inside each drain is monitored using ultrasonic sensors connected through a low-power wide-area network (LPWAN), more specifically a LoRa network [START_REF] Raza | Low power wide area networks: An overview[END_REF], [START_REF] Centenaro | Long-range communications in unlicensed bands: the rising stars in the IoT and smart city scenarios[END_REF]. A few measurements are done per 24h, and the data is sent to a server through LoRa gateways, where the collected data from the ensemble of the sensors are processed. To provide a more practical view of a typical network, Figures 2(a) and 2(b) show an illustration of the actual network deployed around Marseille downtown (the Vieux Port district), and a schematic of the network architecture, respectively. Details on the employed sensors and the other parts of the network are provided in the appendix.

Figure 3 shows a simple illustration of a drain in the form of a box-shape container with some forms of waste inside. An ultrasonic sensor is installed on the top of the drain, pointing toward the bottom. By measuring the delay between an emitted pulse and its echo, the distance to the reflective surface is calculated; this latter can correspond to waste, water or the bottom of the drain (if it is empty). This measured distance, denoted by D, is in fact the data that is then transmitted through the WSN, and which is considered hereafter. Note that a small distance signifies a high waste level W , where W = H d -D with H d denoting the drain depth. Fig. 3: Schematic of an ultrasonic sensor in a typical drain, with some waste (in green). The waste level is measured based on the delay between the emitted pulse and the received echo.

D. Characteristics of the Collected Data

Each sensor measures the waste level within a drain at regular (adjustable) time intervals. The resulting univariate measurement data describe the evolution of the waste level corresponding to the specific drain, which can be modeled as a non-stationary, discontinuous time series. Note, the nonstationarity of the collected data is due to the complex context and the role of numerous exogenous and environmental factors such as rain, wind, topography, etc., which are specific to each drain (and hence, to the corresponding time-series), and in addition, change over time depending on the unknown contextual factors. A few typical examples of such behaviors are:

• Drains with small variations of the waste level, i.e., accumulating gradually with waste, illustrated in the example of Figure 4(a); • Drains with sudden variations of the measured distance, e.g., when exposed to large-size wastes, or when it is cleaned, illustrated in the example of Figure 4(b); • Drains subject to rain where the waste level can either decrease, e.g., flush type behavior where the waste is evacuated in the network or in the streets, or increase, i.e., the rain brings more waste to the drain (these cases are not illustrated for the sake of brevity).

In practice, the measured distances by the sensors are corrupted by noise and undesired outliers. In our case, due to working in a hostile environment subject to dust, humidity, etc., the sensors' measurements have a precision of about ±2 cm at most, which is considered as the measurement noise. Moreover, the measured data may include some peaks (i.e., the outliers), which are either because of the chaotic nature of the waste inside the drain (see the example in Figure 5(a)) or the form of the drain itself, resulting in multiple echos of the We have shown in Figure 6 the measured data of a typical drain sensor during an interval of about one month with two measurements per 24h, for a drain of depth H d ≈ 55 cm. As we can observe, the measured data contains a few high peaks. Note that the occurrence of successive peaks can be due to the form and volume of the waste inside the drain. This specific case study has been selected among the collected data in the following since it represents a relatively difficult case, i.e., containing two waste falls (around days 9 and 22, as indicated in Figure 6), as well as several peaks (including two successive peaks), in a relatively short measurement interval (i.e., less than 15 days from the time the drain is subject to waste accumulation).

Obviously, to study the dynamics of waste accumulation within the network, the collected data should first be preprocessed by removing the above-described outliers, which is the aim of this work. Note that, this can represent a real challenge when considering the substantial amount of collected data in the whole network composed of thousands of sensors.

III. RELATED WORK

As explained above, in this work, we particularly focus on the problem of peak detection and their removal from the raw data. As mentioned in the Introduction, there are numerous previously-proposed works on similar topics, which have considered outlier detection in time series [START_REF] Braei | Anomaly detection in univariate time-series: A survey on the state-of-the-art[END_REF], [START_REF] Blázquez-García | A review on outlier/anomaly detection in time series data[END_REF]. However, most of these solutions are not suitable for use in our specific context, where we are concerned with a large quantity of non-stationary and unlabeled univariate data with a complex evolution over time, depending on a number of different (and even unknown) contextual parameters, e.g., rain, wind, slope of the street, proximity to markets, etc. To elucidate the requirement of developing a new solution for the considered application, we present in the following a brief description of Fig. 6: Example of the collected data from a drain sensor with the measured waste level. some classical but also some recently-proposed methods and explain their inappropriateness for use in our case.

A simple way of removing peaks from data would be to do some sort of low pass filtering to remove the "high frequency" noise. However, this cannot be used in our case due to the nature of the measurements as signal discontinuities, inherent in the data, will be discarded. Even filters that preserve signal discontinuities such as a median filter [START_REF] Justusson | Median filtering: Statistical properties[END_REF] will result in a loss of useful information due to data smoothing since it cannot distinguish between temporarily local events in the data (useful information here) and the noisy peaks. Moreover, data smoothing such as exponential smoothing [START_REF] Hyndman | Forecasting with exponential smoothing. The state space approach[END_REF] will bias any post-classification or prediction of the data, which is, of course, undesirable. For solutions such as adaptive piecewise constant approximation of a time series [START_REF] Dani | Adaptive Threshold for Anomaly Detection Using Time Series Segmentation[END_REF], where outliers are removed within each data segment independently, the performance highly depends on the number of segments initially set, which is difficult to optimize in our case.

Concerning statistical anomaly detection methods proposed in the literature, which are mostly prediction based, the classical approaches based on auto-regressive (AR), moving average (MA), or ARMA models require second order stationarity (i.e., in terms of mean and variance), which is not valid in our case. It is the same for AR integrated MA (ARIMA) method, which requires stationarity in the sense of signal variance [START_REF] Wei | Time Series Analysis: Univariate and Multivariate Methods[END_REF]. Also, the iterative outlier removing approach based on extreme studentized deviate test requires the number of outliers to be known [START_REF] Hochenbaum | Automatic anomaly detection in the cloud via statistical learning[END_REF].

Some other investigated approaches for outlier detection are based on machine learning. There, most of the proposed algorithms such as one-class support vector machines (OC-SVM) [START_REF] Ma | Time-series novelty detection using one-class support vector machines[END_REF] are either supervised or semi-supervised, which require labeled outliers or at least a set of clean data (without outliers), and hence, cannot be applied to the raw data in our case. A few unsupervised outlier detection algorithms have also been proposed for time series, such as peer-group analysis [START_REF] Ferdousi | Unsupervised outlier detection in time series data[END_REF], which consist in characterizing an expected pattern of behavior between similar objects. However, in our case, almost every drain has its own characteristics, and could be considered as independent from the others. Another classical unsupervised technique is that of sub-sequence time series clustering (STSC) [START_REF] Peker | Subsequence time series (sts) clustering techniques for meaningful pattern discovery[END_REF], which consist in applying the Kmeans algorithm to time series using a sliding window. Again, each drain having its specific characteristics, the parameter K (the number of clusters in the K-means algorithm) needs to be set in a personalized way, thus resulting in prohibitive data processing complexity in our case. Another solution can be to use auto-encoder based methods, which is a special artificial neural network used to learn an "encoding" (i.e., an optimal representation of the signal) in an unsupervised manner, e.g., using long-short term memory (LSTM) based auto-encoder [START_REF] Provotar | Unsupervised anomaly detection in time series using lstm-based autoencoders[END_REF]. However, such methods usually necessitate a significant volume of data to optimize the parameters of the underlying neural network, e.g., number of layers, number of cells in each layer, window size, smoothing window size, etc.

IV. OUTLIER CLEANING FROM COLLECTED DATA

In general sense, outlier detection can be considered as a classification problem, where a common approach consists in transforming the time series into a scatter point representation that facilitates separation between normal and abnormal points. Our proposed solution is based on this idea, which is in fact used in some of the methods described in the previous section, where this transformation has been done using sliding windows. Our proposed method, however, realizes this in a much simpler way.

In fact, in order to transform the time series data into a scattered point representation, the idea here is to give a score to each point in order to characterize its "abnormality." Then, different metrics can be used to detect the noisy peaks, as described in the following.

A. Z-Score Method

Z-score is one of the most classical metrics used for anomaly detection [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. Given a data-set of size N , X = [X 1 , • • • , X N ] with mean µ and standard deviation σ, the associated Z-score to each measurement X i , denoted here by Z(X i ), is given by Z(X i ) = (X i -µ)/σ. Note, it implicitly assumes a Gaussian distribution for the data. The higher the Z-score, the higher is the probability that a given measurement corresponds to an outlier. Typically, outlier detection is done using a threshold: if a measurement's Z-score exceeds the threshold, it is considered as an outlier. This will obviously eliminate a certain percentage of the data, for example, a threshold of 2 and 3 results in discarding ∼ 4.6% and ∼ 0.3% of the data, respectively. In our case, considering a rough estimate of noisy peaks of 1% of the collected data, the corresponding threshold is about 2.5, based on the assumption of Gaussian distribution for X.

Since in our case we are concerned with non-stationary data (i.e., changing µ and σ over time), the Z-score method should be applied on a sliding window of appropriate length, depending on the dynamics of waste accumulation, in order to adjust the outlier detection threshold. As an example, we have illustrated in Figure 7 the measured data of Figure 6, where the Z-score method is applied either to the whole data set (see Figure 7(a)) or to a sliding window of length 10 corresponding to 5-day measurements (see Figure 7(b)), considering a threshold of 2.5. Obviously, decreasing the threshold will improve peak detection but it will also result in false peak detection, and vice versa. It is worth mentioning that, as the statistics of waste accumulation are different from one drain to another, the considered threshold should be set in a "personalized" way for every drain and further adapted to environmental and seasonal conditions, for instance, due to the non-stationarity of the waste accumulation process. However, performing such a personalized and adaptive thresholding is challenging and would require, for instance, to go through change point detection in order to delimit the changes in the distribution, resulting in considerable data processing complexity.

B. Opposite Variation Detection

Looking for a better outlier detection, we firstly developed a simple intuitive method, which consist in identifying outliers as corresponding to an "unusual" change in the trend of the measured distance, i.e., a sudden large (negative or positive) variation of D. We will refer to this method as opposite variation detection (OVD). Let us consider the time derivation of X at instance i as follows:

δ i = X i+1 -X i ∆ t , (1)
with ∆ t representing the time interval between two successive measurements. To identify a change in the trend of X, we look for a change in the sign of δ i by defining κ i as follows:

κ i = sign(δ i+1 ) -sign(δ i ) 2 , (2)
where sign(.) is the sign function. Here, κ i takes the values ±1 if the sign of δ i changes, and zero otherwise; X i is considered as a potential peak if κ i = 0. Based on this, we define the outlier detection metric θ i as follows:

θ i = κ i min(|δ i+1 |, |δ i |), (3)
where |.| denotes the absolute value. The idea behind ( 3) is to associate with each point an abnormality score, which is obtained from the absolute value of the local slopes of the D plot (with respect to the previous and the next data points). This latter is the minimum of the two calculated slopes in order to account for signal discontinuities due to normal variations of the waste level in the drain or the measurement noise. This way, "potential peaks" would correspond to relatively large values of θ i .

To explain this more clearly, two examples are shown in Figure 8, where in Figure 8(a) P 1 is a "normal" point whereas in Figure 8(b) P 2 is a peak to detect. Here, both points can be considered as potential peaks since from (2) the corresponding κ are non-zero. The metric θ is in fact calculated so as to distinguish between these two cases. Let us denote the corresponding metric by θ P1 and θ P2 . The slopes of the D plot at P 1 , i.e., δ P1 are denoted by SP 1,1 and SP 1,2 ; likewise the slopes at P 2 are denoted by SP 2,1 and SP 2,2 . Here, for point P 1 , the smallest slope is SP 1,2 , and hence, θ P1 = |SP 1,2 |. Also, for point P 2 , the smallest slope is SP 2,1 , and hence, θ P2 = |SP 2,1 |. We notice that θ P1 < θ P2 , in other words, P 2 is more likely to be a peak than P 1 , which is indeed the case here.

We have further shown in Figure 9 an illustration of outlier detection by calculating the metric θ and applying a threshold to it. Comparing this figure with Figure 6, one can see that for a threshold of 2.5, the OVD method has detected two peaks, but has missed the two others. In fact, the OVD metric is close to zero for the case of successive peaks, which is, of course, unsuitable.

C. Proposed Solution

As discussed in the two previous subsections, both sliding window Z-score-based and OVD methods show limited performance in practice in detecting outliers in the collected data. In fact, these methods are based on metrics which are centered at 0 when applied to "clean" data. Here, we propose a more efficient feature-based solution, that we will call peakpattern-based Z-score (PPZ), by combining the two ideas and applying them to a two-dimensional (2D) representation of data, as described in the following. The idea behind this 2D feature-based approach is to complete the information provided by Z-score method with that obtained from OVD.

We propose to use the two previously-presented metrics to obtain a scatter plot of data, centered at (0, 0). The points corresponding to outliers will be placed distant from the center, which can be distinguished from the clean data by applying a 2D threshold, i.e., in the form of an ellipse in the 2D representation. In other words, for each point X i of the time series, Z-score and θ metrics are calculated, resulting in two vectors Z and θ, respectively. Threshold setting for outlier detection is based on eigen-decomposition of the (2 × 2) covariance matrix K Z,θ , called standard deviational ellipse (SDE) [START_REF] Wang | Confidence analysis of standard deviational ellipse and its extension into higher dimensional euclidean space[END_REF]; the resulting eigenvectors determine the angle of the threshold ellipse, while the axes length of this latter is set by multiplying the square root of the corresponding eigenvalues by a constant ξ.

In fact the proposed approach relies on the simplifying assumption that the time series X, modeled as samples of a random variable, follows a Gaussian distribution whose mean and standard variation are almost constant within the time window of data collection. For instance, the considered period in the provided examples (see Figure 6) is 120 hours, i.e., five days, which corresponds to 10 points with sensing rate of 2 measurements per day. Under this assumption, the calculated Z-score will also follow a Gaussian distribution. As for θ metrics, it will follow a pseudo-Gaussian distribution, i.e., with a higher occurrences of zeros due to κ.

Figure 10 illustrates outlier detection using the proposed PPZ method, applied to the set of collected data from Figure 6. Fig. 9: OVD algorithm applied to the data in Figure 6, considering a threshold of 2.5 (in absolute value) on the calculated metric θ.

Here, the constant ξ is fixed to 2.5, as previously considered for Z-score and OVD methods. As can be seen, each point in the 2D plot corresponds to the calculated Z-score and θ metrics for every measurement (as already shown in Figures 7(b) and 9). The method allows the detection of all peaks, albeit incorrectly identifying the point at Day 4 as a peak.

V. PERFORMANCE EVALUATION OF THE PROPOSED ALGORITHM

The common approach for comparing the performances of different outlier detection methods is to contrast the ROC (Receiver Operating Characteristic) curves. These are obtained upon calculation of the confusion matrix, with elements consisting of the number of "true positives" N TP (i.e., correctly detected outliers), "true negatives" N TN (i.e., correctly notdetected normal points), "false positives" N FP (i.e., normal points identified incorrectly as outliers) and "false negatives" N FN (i.e., undetected outliers). ROC curves are obtained by plotting the so-called "true positive rate" TPR, defined as N TP /(N TP +N FN ), versus the "false positive rate" FPR, defined as N FP /(N FP + N TN ). 1 Performance evaluation can be done based on either simulated or real labeled data. Obviously, the former approach is pertinent only when the simulated data are representative of the real data. This appears to be challenging in our case since each drain has a distinct behavior, related to a number of environmental parameters. Nevertheless, as a preliminary performance study of the proposed algorithm, we applied it to a set of simulated data, generated based on a piecewise constant or a piecewise linear time series (based on some arbitrary parameters and probabilistic laws) to which a Gaussian noise was added to represent the sensor measurement noise as well as a few peak-type points as outliers (e.g., using an exponential 1 Note, other metrics could also be used to evaluate the performance of the algorithms including the so-called "precision" P r defined as N TP /(N TP + N FP ), or the F 1-score, defined as the harmonic mean between TPR and precision, i.e., 2 P r TPR/(P r + TPR).

distribution with an offset). 2 Applying the proposed algorithm to such generated data was quite promising (results are not shown for the sake of brevity). However, this cannot be considered as representative of the performance on the real data, given the complexity of the collected data in practice. Therefore, we have decided to evaluate the algorithm performance by applying it to real (measured) data. This, however, necessitates manual labeling of a data set in order to identify with certitude the peaks present in the data (and to see whether or not they are actually detected through outlier detection).

For this purpose, a data set of two months from around 300 sensors was labeled, representing a total of about 33, 600 points, among which over 460 points (i.e., 1.3% of the data) were labeled as peaks. To label the data manually, for each sensor, the time series was visualized as a scatter plot, where each point was labeled or not as a peak (this was a tedious task, of course). We have contrasted the corresponding ROC curves of the proposed PPZ method with those of Z-score and OVD in Figure 11, where we can notice the superiority of the former. Note, for the PPZ algorithm, the threshold parameter ξ can be set so as to ensure a required minimum TPR or a maximum FPR. A good compromise can be made by setting ξ = 2.5, which results in a TPR of 0.85 and an FPR of 2.5 × 10 2 , as shown in the figure. Note, the performances of the algorithms can further be quantified by calculating the air under each ROC curve. In this case, the airs are 0.92, 0.95 and 0.98 for the Z-score, OVD, and PPZ methods, respectively.

Lastly, concerning the computational complexity of our proposed method, for a sensor with N measurements, the calculation of each metric has a complexity of O(N ). Then, the algorithm calculates the threshold ellipse, which is based on eigen-decomposition of the covariance matrix of dimension (2 × 2), as described in Subsection IV-C, with relatively low complexity [START_REF] Arora | Computational Complexity: A Modern Approach[END_REF]. Meanwhile, the estimation of each of the 2 More specifically, two approaches were used for data generation in these simulations. The first approach consisted of generating data according to a piecewise constant signal level of length N divided into P equal segments. The value of each segment was set to an integer, generated randomly according to a uniform distribution U (0, 100). By the second approach, data were generated according to a piecewise linear signal level variation, where the slope of each segment was generated randomly following a uniform distribution U (-1, 1). In both cases, a normalized Gaussian noise was added to each point to represent the measurement noise; then, P points, selected randomly, were modified by adding an exponential noise with an offset to simulate the outlier peaks. four entries of the covariance matrix entails a complexity on the order of O(N ). So, overall, the proposed method has a computational complexity on the order of O(N ).

VI. CONCLUSIONS AND DISCUSSIONS

We proposed in this paper a new algorithm for peak-type outlier detection applied to sewer network collected data. The proposed PPZ method provides superior performance compared to the classical Z-score method with a relatively low additional compositional complexity. The main idea behind this algorithm has been to enhance the performance of the Zscore metric by augmenting it with a pattern-detection metric designed here to detect peak-type outliers. In other words, the Z-score metric provides a large spectrum outlier detection but necessitates the signal to follow a Gaussian distribution as well as the use of an optimal sliding window length. However, the arbitrary distribution of the signal and its non-stationarity make theses two requirements difficult to achieve in general. On the other hand, by the OVD method, the θ metric was specifically designed to detect peak-type outliers, independently from the window size. As such, combining the complementary metrics of Z-score and OVD through 2D thresholding for outlier detection results in considerable performance improvement, while entailing a relatively low computational complexity.

Despite the advantages of the proposed method, highlighted in the paper, it could have limited performance in specific circumstances in practice. The first consideration concerns the calculation of the metric θ for the time series, which requires that the measurement rate of the sensor should be the same for the considered time interval (remember that the sensing rate of each individual sensor can be modified, if needed, to have more or less precise evolution of the waste accumulation within the corresponding drain). Similarly, the proposed approach assumes that the parameters of the collected data are almost unchanged during the considered interval. Any abrupt change in the properties of the data (e.g., flash function after waste accumulation due to heavy rain, or drain cleaning) will affect the performance of the algorithm. Moreover, the selection of the time window width for the Z-score metric (taken as 120 hours in the presented results) should be adapted for each specific sensor and the dynamics of waste accumulation in the corresponding drain.

To improve the outlier detection performance, future work will consider the use of more elaborate clustering methods (in contrast to the considered 2D SDE), such as one-class-supportvector-machine (OCSVM) [START_REF] Ma | Time-series novelty detection using one-class support vector machines[END_REF], when a sufficient amount of clean (outlier free) data is available, required for algorithm training. This has additionally the advantage of developing a boundary function adapted to each sensor. Also, in the case of recurrent successive peaks, a customized metric can be developed, adapted to such specific patterns.

APPENDIX

The network uses an LPWAN network for data transmission from the sensors, that can be according to LoRa (used here) or Sigfox technologies. The LoRa gateways are collocated with network servers. The distance between a sensor and the corresponding BS in our networks varies between a few meters up to a few hundreds of meters. The collected data are transferred to a cloud via the network server for storage in a database and processing. The users can access the data via a web application connected to the cloud. More details on the network and sensors' specifications are provided in the following.

• Sensor: Ultrasonic level measurement sensor ("Hummbox Level" third generation, made by GreenCityZen Co., based on Microchip's SAM microcontroller with its development tools and programming environment);

Fig. 1 :

 1 Fig. 1: Illustrating the three existing types of sewer networks in the city of Marseille.

Fig. 2 :

 2 Fig. 2: Illustration of the LoRa Network: (a) Part of the sewer monitoring network deployed in the Vieux Port district at downtown of Marseille with green circles showing the location of the connected drains and towers representing LoRa gateways through which the data is sent onto a cloud. (b) Schematic of the network architecture.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Simple illustration of the evolution of the waste level W and the measured distance D by a sensor; H d = 50 cm; (a) slow and steady variations, (b) sudden variations of W .

Fig. 7 :

 7 Fig. 7: Example of outlier detection using Z-score method on the data showed in Figure 6. (a) Z-score of the data calculated on the whole window (no peaks are detected); (b) Z-score calculated on a sliding window of length 10 points (5 days), in this case only one peak is correctly detected. The outlier detection threshold is considered at ±2.5.

Fig. 8 :

 8 Fig. 8: Illustration of the idea behind the θ metric: (a) Normal point, P 1 , (b) a peak to detect, P 2 . We have SP 1,2 < SP 1,1 and SP 2,1 < SP 2,2 . Also, |θ(P 1 )| < |θ(P 2 )|.

Fig. 10 :

 10 Fig. 10: Illustration of outlier detection using the proposed outlier detection algorithm applied to the data of Figure 6 considering ξ = 2.5: (a) scatter plot of Z-score and θ metrics; (b) the resulting detected peaks in the original time series.

Fig. 11 :

 11 Fig. 11: Contrasting ROC curves of the Z-score, OVD, and the proposed PPZ method applied to labeled data.

•

  Cloud: The Things Network (TTN); • Wireless connectivity: LoRa; -Carrier frequency: 868 MHz, -Transmission power: 14 dBm, -Effective Isotropic Radiated Power (EIRP): 16 dBm, -Receiver sensitivity: -148 dBm, -Number of channels: 3 to 8.


Sélim Chefchaouni Moussaoui

Hugo Pousseur

Alessandro Corrêa

Marie-Hélène Abel

Dynamic Context Awareness in Autonomous Navigation

Keywords:

Many studies faced the problem of vehicle autonomous navigation in different fields, but nowadays just a few of them uses all the implicit information coming from the context in which such navigation is occurring. This results in a huge potential information loss that prevents us from adapting the vehicle's behavior to each different situation it may be in. In a previous work, we defined a method to model the static context of navigation using ontologies and take it into account in the command law when performing a local navigation task. In this paper, we extend our model of the context of navigation, and define a software architecture able to update the context dynamically, by using sensor information. The method is tested with real-time experiments on driving simulator. They show that the Context of Navigation can be effectively updated during the navigation and leads to a smarter vehicle's behavior on the road.

I. INTRODUCTION

The problem of autonomous navigation has been studied for years now, and nowadays many functioning solutions for it have been proposed and approved. These solutions generally allow vehicles to perform both local tasks (e.g. going through an intersection or avoiding an obstacle) using sensor information, and global tasks (e.g. following a desired itinerary) using other information sources like maps. By combining local and global tasks, efficient autonomous navigation can be performed nowadays, such that autonomous navigation solutions are now included by car manufacturers.

However, autonomous navigation systems do not take into account the contextual information coming from the vehicle and its environment. It follows that the vehicle has the same behavior regardless of the state of the passengers (e.g. pregnant or injured passenger) or the vehicle components (e.g. battery level). This contextual information which is not directly coming from the sensors constitute what we call the context of navigation. Being able to take it into account in the control law would lead to a smarter navigation that adapts better to the vehicle's environment.

In previous works [START_REF] Faruffini | Vehicle Autonomous Navigation with Context Awareness[END_REF] [2], we have shown that ontologies are an efficient tool to model such contextual information, as they can be used with rules to provide driving suggestions from the modeled context. After that, we developed a control method based on visual servoing which was able to take these suggestions into account in the control law, when performing a local navigation task.

However, the applications of these works were limited, because the context was defined in a static way and could not be updated during the driving. This static definition was sufficient for some elements that are not regularly updated, like the passenger's state. But it appears as a strong limitation when we want to consider contextual elements that may change in real-time, like the presence of an obstacle. The purpose of this paper is to provide an solution to this problem, by proposing a wider context model, and a framework able to update dynamically the knowledge base obtained from the context of navigation ontology and to take it into account in real-time in the control law. It appears as the natural extension of the capabilities of context aware navigation, towards its consideration and its application in autonomous vehicles.

The paper will proceed as follows. We are first going through the state-of-the-art work on context aware navigation in Sec. II. We then propose a new ontology to model the context of navigation that considers the global information in Sec. III. We then remind the handling of local context in Sec. IV. After this, we focus on the applicative aspects, by explaining the dynamic consideration of the context in Sec. V. Finally, Sec. VI is dedicated to the tests and results obtained in simulation and is followed by a conclusion in Sec. VII.

II. RELATED WORK

Some ontology-based solutions have been proposed in the past as the premises of context aware navigation solutions. Many of them were focused on Advanced Driving Assistance Systems (ADAS) solutions, i.e. electronic systems built to help the driver with parking and driving functions. Zhao et al. proposed an ontology-based system to prevent the driver from exceeding the speed limit [START_REF] Zhao | An ontology-based intelligent speed adaptation system for autonomous cars[END_REF] or to help him in intersections [START_REF] Zhao | Ontology-based decision making on uncontrolled intersections and narrow roads[END_REF]. Their model was based on three different ontologies taking into account information on different levels: the Map Ontology, the Control Ontology and the Car Ontology.

More recently, Armand et al. in [START_REF] Armand | Ontology-based context awareness for driving assistance systems[END_REF] modelled the inference rules for an ADAS system to allow the use of reasoners to reason over the ontology representing the contextual information, in the case of ADAS. Thus, the reasoner was able to propose suggestions for the driving assistance system on the decisions it should take. From these previous works, Faruffini et al. proposed a solution to take into account these reasoner suggestions directly in the control law [START_REF] Faruffini | Vehicle Autonomous Navigation with Context Awareness[END_REF]. The approach is a modification of an image-based (or visual servoing) controller that modulates its output to consider more or less the contextual suggestion. This solution has been demonstrated on simulated vehicles, and has proved to be an efficient solution to consider the context in a local navigation law. This approach is reminded in Sec. IV.

From this state of the art, we will propose a generic framework for context-aware autonomous driving, able to infer in either global and local navigation law, and a related ontology for the context model. In this paper, the developed framework will focus on local navigation by taking into account the dynamic evolution of the navigation context.

III. CONTEXT MODELING

In this section, we detail the ontology-based model developed for the context of navigation, which uses elements of the Navigation Context Ontology (NCO) created in [START_REF] Faruffini | Context Modelling applied to the Intelligent Vehicle Navigation[END_REF].

A. Ontologies

Before going through the model, we need to remind some definitions and properties of ontologies and reasoners, that are the tools used to model the context of the vehicle autonomous navigation.

An ontology can be defined as a semantic structure to store real world data with a formal representation. This kind of representation has been originally developed for the semantic web, which aimed to make web pages machine-readable. More recently, the use of ontologies for robotic applications have been considered as stated in the previous section. Given this historical consideration, we model ontologies using the tools originally developed for the semantic web. These tools are the Ontology Web Language (OWL) that is the language describing the ontology, the Semantic Web Rule Language (SWRL) that describes the rules applied to our ontology, and the Protege software [START_REF] Musen | The protégé project: a look back and a look forward[END_REF] as a front-end ontology editor.

To model a system using ontologies, we need to define classes corresponding to the semantic concepts, and properties that are the links between the classes (object properties) or between classes and numerical data (data properties). These abstract definitions are together forming the Terminological Box (TBox). From these definitions, we can define individuals that are equivalent to the class instances. These individuals constitute what we call the Assertional Box (ABox). Together, TBox and ABox form a Knowledge Base (KB).

The main benefit of using an ontology structure is the possibility to exploit reasoners. Reasoners are software that can operate logical reasoning on the KB to make deductions that provide new information on it. That is the reason why an ontology-based system is adapted to our problem and the ones mentioned in the previous section. More generally, ontologybased systems are useful tools for proposing recommendations or suggestions to their users. That is the way ontology has been used in previously mentioned works and is being used in this work.

B. Context model

The Context of Navigation has two components: the Static Context and the Dynamic Context. The ontology we modeled in [START_REF] Faruffini | Context Modelling applied to the Intelligent Vehicle Navigation[END_REF] was focused on the modeling of the Static Context, i.e. elements that are not updated in real-time. These were mainly information on the passengers' states (to adapt the driving to a pregnant or injured passenger for example) or on the vehicle state (its loads for example). Here, we consider again these static elements, and we add information that vary during the trip with respect to the navigation: this is the Dynamic Context, which models for example the vehicle speed or the presence of an obstacle in front of the vehicle. These elements require to be able to update the context dynamically (when the vehicle enters a new type of road), at varying moments during the

Attribute Meaning CP DESCRIPTION NodeID

The ID of the CP, numbered from the beginning to the end of the itinerary Lat

The latitude of the CP Lon

The longitude of the CP Vmax

The maximum speed limit of the road R lanes

The number of available lanes in the road Rway

The road classification (one way road or not) Ra

The CP classification (roundabout or not) N roads

The number of roads at the CP CP GEOMETRY Θ 1 to Θ N roads

The angles for each road at the CP, ordered counterclockwise

D 1 to D N roads
The directions for each road at the CP, ordered counterclockwise travel. These dynamic elements are used on our test scenario (see Sec. VI).

In addition, we extended the definition of our ontology by including classes that model global navigation elements that can be used by a reasoner to provide suggestions for the global navigation. The global navigation elements modeled in our ontology are places and itineraries. A place is simply modeled by its coordinates, from which we can compute the distance to other places or its presence in an itinerary. For itineraries, we use the same model as [START_REF] Kenmogne | A Global Navigation Management Architecture Applied to Autonomous Robots in Urban Environments[END_REF]. In this model, an itinerary is represented as a suite of intersections (called Critical Points) and corridors. The description of intersection attributes is detailed in Tab. I .

We transcribed the information we want to model in OWL using the Protege software. The overview of the ontology classes, properties and individuals is given in Fig. 1. The same implementation is also described in Fig. 2 in the form of a UML class diagram, to understand better the links between classes and properties. This model presents several local and global elements that will be updated at a certain frequency (e.g. the road state) or when a specific event occurs (e.g. the detection of an obstacle). The OntologyHandler module presented in Sec. V is in charge of separating these elements.

In addition to the ontological model, we created a set of SWRL rules which allow the reasoner to suggest a maximum vehicle speed, that is taken into account by our controller. An example rule is given in Sec. IV.

IV. LOCAL CONTEXT AWARENESS

In this section, we remind briefly how we embedded context awareness in a local driving problem.

A. Problem formulation

We want to control our vehicle to make it drive autonomously in the middle of its lane, while avoiding possible obstacles. The control of the vehicle is given as u = [v ω] T where v represents the longitudinal speed of the car and ω the angular velocity. We use the kinematic model of a front-wheel car with Ackerman's approximation for the steering angle ϕ, which is written as follows:

    ẋ ẏ θ φ    =     cos (θ) • cos (ϕ) sin (θ) • cos (ϕ) sin (θ)/l 0     v +     x0 0 0 1     ω (1)
where x and y represent the vehicle's coordinates in a 2D plane, while θ and φ represent its orientation and yaw rate.

In a previous work, D. Fox et al. developed an approach called Dynamic Window Approach (DWA) [START_REF] Fox | The dynamic window approach to collision avoidance[END_REF]. It consists in finding the couple (v, ω) that maximizes the following objective function:

DWA(v, ω) = α • heading(v, ω) + β • dist(v, ω) + γ • velocity(v, ω) (2)
where α, β and γ are parameters to be tuned. The heading function is based on the final orientation of the robot to lead it to the final goal position. The dist function prioritizes movement in areas free of obstacles, while the velocity function prioritizes longitudinal speeds close to the desired constant set point.

This method was then adapted in D. Lima [START_REF] Lima | A hybrid controller for vision-based navigation of autonomous vehicles in urban environments[END_REF] to be able to compute the heading function based on visual information. That is the Image-based Dynamic Window Approach (IDWA). This autonomous navigation controller is adapted in this paper, in order to take into account the context of the navigation. This is done by adding in this controller framework the information provided by a reasoner, as described in the next section.

B. Reasoner output

As said in Sec. III, one of the main advantages of the ontologies is the possibility to use rules that allow reasoning on the KB. Thus, rules can be used to provide suggestions given a KB. In our case, the reasoner results in a suggested speed, that is used by the controller as explained after.

We define a set of rules to provide speed suggestion using the SWRL language. For example, the following rule aims to provide suggested speed w.r.t. the passenger preferences:

EgoVehicle(?v) ∧ SpeedSuggestion(?sg) ∧ hasPassenger(?v, ?p) ∧ hasDrivingStylePreference(?p, ?dsp) ∧ hasRoadType(?rp, "CityRoad") ∧ hasCityRoadSpeedPreference(?sp) ∧ isOnRoadPart(?v, ?rp) → hasSuggestedSpeed(?sg, ?sp)

The suggested speed is then used by the ICDWA robotic control loop, as described in the next section.

C. Image and Context-based Dynamic Window Approach

To take the suggested speed into account in the control law, F. Faruffini et al. defined a new approach called Image and Context-based Dynamic Window Approach (ICDWA) [START_REF] Faruffini | Vehicle Autonomous Navigation with Context Awareness[END_REF].

To allow the ICDWA to take the suggested speed v CN into account, we have modified the IDWA by adding a fourth function to the maximization problem. It results:

ICDWA(v, ω) = α • heading(v, ω) + β • dist(v, ω) + γ • velocity(v, ω) + δ • context(v, ω) (3)
The parameter γ represents at which point the suggested speed is considered by the controller: for γ = 0 the suggestion is not considered at all, for γ = 1 it is the only speed to be considered.

The context(v, ω) function is decoupled as:

context(v, ω) = ε • (context v (v) + context ω (ω)) (4)
where ε ∈ [0, 1] and context v (v) is expressed as follow (the principle is the same for context ω (ω)):

context v (v) =    σ(x) if x ̸ = 0 1 2 + 1 4
x otherwise [START_REF] Armand | Ontology-based context awareness for driving assistance systems[END_REF] where σ(x) represents the sigmoid function and x the difference between the current speed and the speed suggested by the reasoner:

x = v -v CN .

V. DYNAMIC CONTEXT UPDATE

This section focuses on the dynamic update of the context of navigation and our software implementation of this update.

To benefit from the advantages of reasoning, the context used by the reasoner should be dynamically updated. The principle of context update in the case of the ICDWA is shown in Fig. 3. We are now able to update the individuals in the ontology used by the reasoner.

To be able to perform such an update, we developed a Python module named OntologyHandler. It uses the owlready2 API by J-B. Lamy [START_REF] Lamy | Ontology-oriented programming in python with automatic classification and high level constructs for biomedical ontologies[END_REF]. It has the capabilities of initializing a KB from a configuration file, updating the KB in real-time, and get the KB at any time, including the reasoner output. The advantage of owlready2 compared to the methods tested in [START_REF] Faruffini | Context Modelling applied to the Intelligent Vehicle Navigation[END_REF] is that it accesses directly the whole OWL file. For the reasoning part, we use a modified version of the Pellet reasoner [START_REF] Sirin | Pellet: A practical owl-dl reasoner[END_REF] used included in owlready2, as it is the only known reasoner able to infer new data from SWRL rules.

The OntologyHandler module is included into a wider framework for context aware navigation. The framework is based on ROS [START_REF]Robotic operating system[END_REF] and implements local and global contextaware navigation capabilities. The details of the framework architecture are presented Fig. 4.

If the framework is able to handle both local and global navigation, the requirements for these navigation approaches are not the same when it comes to the consideration of contextual information. For global navigation the context must be updated regularly, but updates can be easily separated by several seconds. However, the handling of local navigation tasks requires to update the context with a high frequency, as it has a direct impact on the vehicle's speed. Then, real-time considerations must be tackled.

In the framework architecture presented in Fig. 4, the different modules -including the reasoning module OntologyHandler -are executed in separated threads. This feature allows us to choose the execution frequency for the reasoner, in function of the real-time requirements.

VI. TESTS AND RESULTS

This section aims to show the tests done on our developments about local dynamic context awareness. These tests have been realized under the professional driving simulator SCANeR Studio, ver. 2022.1, which provides a ROS bridge that can be used with our architecture, like we tested on a real vehicle.

A. Test scenario

For the tests, we define a scenario that requires dynamic updates, in order to test the principle of dynamic context awareness. In this scenario, we have an intelligent car, driving in a city (thus the speed limitation is 50km/h), on a straight road, with an obstacle at 150m inf front of it in its lane that should be avoided. The principle of the test scenario is shown in Fig. 5.

The vehicle is simulated in SCANeR Studio and communicates via ROS to our contextual navigation framework. To be able to avoid effectively the obstacle, the vehicle should detect the lanes. For this, we interface the deep-learning-based lane detector by A. Shour et al. [START_REF] Shour | Shared Decision-Making Forward an Autonomous Navigation for Intelligent Vehicles[END_REF] to our simulation environment.

Concerning the context, the KB is initialized from a configuration file, which states there is a single passenger, with the following speed preferences:

• 110 km/h on highways • 70 km/h on countryside roads • 40 km/h on cities • 5 km/h when crossing an obstacle As the scenario alternates between city road an the crossing of an obstacle, it is sufficient to test the principle of our context update.

B. Dynamic context update

We want to test if the software architecture detailed in the previous section allow to update the context efficiently. The software has been implemented in Python for ROS and shows to be functional. As stated, we modulated the execution frequency of the OntologyHandler module. For the tests, we used a frequency of 20Hz, which appears to be sufficient for our scenario.

C. Local context awareness

The ICDWA method was already tested in [START_REF] Faruffini | Vehicle Autonomous Navigation with Context Awareness[END_REF]. Here, we test it with our new scenario and software architecture.

When driving in the city and avoiding the obstacle, the vehicle's speed is adapted from the passenger's preferences. It results to the velocity profiles of Fig. 6. We used varying values for the coefficient γ of the ICDWA function. As the other parameters α, β and δ, the parameter γ is defined statically in the controller.

The profiles of Fig. 6 show that the velocity is reduced by taking into account the passenger preference. Besides, we can see the impact of the coefficient γ on the context awareness. When γ = 0.5, the suggestion is taken partially, especially when it comes to avoiding the obstacle: the speed is reduced to 30km/h only, and the situation becomes dangerous. When we increase the coefficient to have γ = 0.75, the velocity profile matches much more the speed suggested by the reasoner. However, the absolute value of the acceleration is higher for this profile, which is something we want to avoid too. Thus, a compromise has to be found for the choice of the γ coefficient. This study is being performed in our going works.

D. Discussion

From the results presented in Fig. 6, we show that the speed is reduced and does not exceed the driver's speed preferences. The γ parameter seems to be well-tuned between 0.3 and 0.6. The optimal tuning of this parameter could lead to further works on the context consideration.

Due to hardware limitations and a lack of optimization in the bridge between the simulator and the framework, the simulation was performed at a very slow rate (between 4Hz and 10Hz) that caused a few instabilities for the ICDWA controller and high instabilities for the lane detector. To face this limitations, the speeds (both passenger's preferred speed and speed limitation of the road) have been reduced by 5 for the simulation. This issue should be targeted for future implementations on real vehicles.

Our test shows how contextual information can be updated in the ontology of the context of navigation from the vehicle's sensors, integrating dynamic information to the KB, and then inferred back on the vehicle control law through the ICDWA controller. The example displayed in this paper are simple but functional, and it is possible to extend them towards a much more complex and complete context aware navigation. The software framework and the ontology have been designed to be generic enough to handle new kind of information, including information on global navigation, that could be included and tested in the future.

VII. CONCLUSION AND PERSPECTIVES

In this paper, we presented a method to model contextual information and add it dynamically to the navigation law of an autonomous vehicle, in order to provide a smarter vehicle behavior on the road. We also defined a software framework capable of handling the context by interacting with the ontology of the context of navigation in real-time during the navigation of the vehicle. This framework was tested on simulation with a new scenario, and proved that the context can be effectively updated to infer dynamically to the vehicle behavior.

The promising results obtained proved that this approach could improve the autonomous vehicle future developments. In further works, a much wider ontology of the context of navigation could be developed, including information on the global navigation, in addition to new methods to handle this kind of information. Besides, the ICDWA approach used in this work could be improved too, by choosing wisely at which point the contextual suggestions should be considered, with respect to the safety norms and considerations and to the navigation integrity. Finally, the choice of an adequate moment for updating each element of the context could be defined to avoid useless executions of the reasoner. These points are essential for future research and need to be targeted in the future developments, towards the consideration of the semantic modeling of the context of navigation.

Fig. 1 .

 1 Fig. 1. The OWL ontology of the navigation context shown in Protege. The OWL file created in Protege handles the classes (a), data properties (b) and object properties (c). Individuals are created by the OntologyHandler module when initializing or updating the ontology to constitute the knowledge base. Individuals shown in (d) are the one created in the initialization.

Fig. 2 .

 2 Fig. 2. The context of navigation implemented in the OWL ontology, presented under the form of a UML class diagram. This representation is useful to understand in a single diagram the links between the different classes and their properties. In the bottom-left corner are displayed the enumerations and the classes that represent only a data type.

Fig. 3 .

 3 Fig.3. ICDWA controller architecture that is able to update dynamically the context of navigation, using sensor data from the vehicle.

Fig. 4 .

 4 Fig. 4. Details of the software architecture developed for the contextaware navigation framework. The framework is capable of running separated instances of the OWL ontology, and is able to update in real-time the ontology individuals, including those who represent the reasoner output. The different submodules displayed in this figure run as separated threads.

Fig. 5 .

 5 Fig. 5. The principle of our test scenario: the vehicle should avoid the obstacle on its lane.

Fig. 6 .

 6 Fig. 6. Velocity profile obtained for the vehicle with our test scenario, with varying coefficient γ. The velocity is reduced by taking into account the passenger preference.

TABLE I THE

 I DESCRIPTION OF CRITICAL POINTS (CP) ATTRIBUTES. EACH CP

REPRESENTS A ROAD INTERSECTION AND IS PART OF AN ITINERARY.

This research was funded in part by the Labex MS2T ("Laboratory of Excellence" on the Control of Technological Systems-of-Systems) a joint research entity between the Université de Technologie de Compiègne (UTC), the French National Centre for Scientific Research (CNRS) and Sorbonne Universités. Research was carried out in Heudiasyc Laboratory, UMR CNRS 7253, in Compiègne, France.