William Schueller
email: william.schueller@inria.fr

Vittorio Loreto

Pierre-Yves Oudeyer

Complexity Reduction in the Negotiation of New Lexical Conventions

Keywords: language emergence, active learning, multi-agent model, control of complexity growth

In the process of collectively inventing new words for new concepts in a population, conflicts can quickly become numerous, in the form of synonymy and homonymy. Remembering all of them could cost too much memory, and remembering too few may slow down the overall process. Is there an efficient behavior that could help balance the two? The Naming Game is a multi-agent computational model for the emergence of language, focusing on the negotiation of new lexical conventions, where a common lexicon self-organizes but going through a phase of high complexity. Previous work has been done on the control of complexity growth in this particular model, by allowing agents to actively choose what they talk about. However, those strategies were relying on ad hoc heuristics highly dependent on fine-tuning of parameters. We define here a new principled measure and a new strategy, based on the beliefs of each agent on the global state of the population. The measure does not rely on heavy computation, and is cognitively plausible. The new strategy yields an efficient control of complexity growth, along with a faster agreement process. Also, we show that short-term memory is enough to build relevant beliefs about the global lexicon.

Motivations

Lexical conventions constitute an important element of social interactions. They can emerge, evolve, or be learnt within a population, without necessarily having a centralized control. In other words, they can be negotiated through local interactions between individuals. In practice, this happens continuously in human societies, being the spread of new words and conventions, the acquisition of those conventions by infants or other learners, or even the emergence of new forms of communication. Despite the high complexity of the processes involved, humans deal with these issues quite efficiently.

Learning of high complexity tasks in individuals can in general be facilitated by an active control of the complexity of learning situations , often driven by intrinsic motivation, like for example maximization of the learning progress [START_REF] Gottlieb | Information-seeking, curiosity, and attention: computational and neural mechanisms[END_REF][START_REF] Baldassarre | Intrinsically motivated learning systems: an overview[END_REF][START_REF] Barto | Intrinsic motivation and reinforcement learning[END_REF][START_REF] Oudeyer | Intrinsic motivation systems for autonomous mental development[END_REF]. This type of mechanism is also argued to be an evolutionary advantage for cognitive abilities [START_REF] Oudeyer | How evolution may work through curiosity-driven developmental process[END_REF], and can also be found in lexicon acquisition at the individual level [START_REF] Partridge | Young Childrens Self-Directed Information Gathering on Touchscreens[END_REF]. But does it have a significant impact on population-wide learning and conventions negotiation dynamics?

The Naming Game [START_REF] Steels | Stochasticity as a source of innovation in language games[END_REF][START_REF] Wellens | Adaptive Strategies in the Emergence of Lexical Systems[END_REF][START_REF] Loreto | Statistical physics of language dynamics[END_REF][START_REF] Ke | Self-organization and selection in the emergence of vocabulary[END_REF] is an adapted framework to test this hypothesis. It is a class of multi-agent models of language emergence and evolution, where pairs of randomly selected individuals try to communicate by referring to some pre-defined meanings using words. At the beginning, they do not share any convention about word-meaning associations. Through repeated decentralized interactions, a common lexicon self-organizes. However, the process can be slow and pass through a high-complexity phase where agents memorize a lot of conflictual information, in the form of synonyms and homonyms.

It has already been shown that active learning mechanisms can increase convergence speed towards a shared lexicon in different language emergence models [START_REF] Cornudella | How intrinsic motivation can speed up language emergence[END_REF][START_REF] Schueller | Active control of complexity growth in naming games: Hearer's choice[END_REF]. The main idea behind those mechanisms is to allow agents to actively choose the topic of their communication, based on information collected during their past interactions and driven by control of complexity growth. However, the algorithms used so far are based on ad hoc heuristics, constrained interaction scenarios and can depend heavily on fine-tuning of parameters themselves depending on population size and number of words and meanings.

In previous work, an approximation of the global state is built by each agent using the information of past interactions, in the form of an average vocabulary of the population [START_REF] Oliphant | Learning and the emergence of coordinated communication[END_REF][START_REF] De | The evolution of conventions in multiagent systems[END_REF]. Is it possible to design a new principled algorithm for an active topic choice based on such a representation? Could decisions be driven by both the compatibility of an agent's own lexicon with this average vocabulary, and a reduction of both their complexities? Such an algorithm should rely on a time scale for the memory of past interactions: Indeed, in the case of uncentralized negotiation of a lexicon, conflictual conventions will necessarily appear and have to be forgotten in order to converge to a functional global vocabulary. Remembering them could slow down the self-organization process.

In this work, we define a principled measure of correlation between an agent's lexicon and a local approximation of the average lexicon of the population. We build a strategy driven by the maximization of this value without being computationally hard, to be cognitively plausible. We study and discuss the impact of this strategy on convergence time and complexity growth, depending on a time scale used for memory.

Methods The Naming Game

We define here precisely the Naming Game model that we used (see fig. 1 for an overview). We need to explicit: -The interaction scenario itself -How agents represent their lexicon -How they update their lexicon at the end of each interaction It is a simple modification of the standard Naming Game scenario [START_REF] Loreto | Statistical physics of language dynamics[END_REF][START_REF] Wellens | Adaptive Strategies in the Emergence of Lexical Systems[END_REF].

Figure 1: Illustration of the Naming Game model. Out of a population of simulated agents, two are picked and try to communicate, using/inventing words to refer to meanings. After repeated such interactions, a common lexicon self-organizes. In this example, there are M=2 possible meanings, W =4 possible words and N=12 agents.

Interaction process We re-use a previously defined interaction process called Speaker's Choice [START_REF] Schueller | Active control of complexity growth in naming games: Hearer's choice[END_REF]. It allows one of the interacting agents, called the speaker, choose actively the topic of the interaction.

Each interaction involves two agents, that are picked randomly from the population. One of them is assigned the role speaker, and the other the role hearer. The speaker chooses a topic and picks up a word for this topic. If it does not have a word associated so far to the meaning used as topic, it just invents a new meaning-word association. It utters this word, which is interpreted by the hearer as a meaning, if it knows this word. If the interpreted meaning is the same as the topic, i.e. the meaning intended by the speaker, the communication is considered successful. Otherwise, it is considered a failure. See fig. 2 for a detailed illustration of the interaction process.

Vocabulary Representation Vocabularies, or lexicons, are a set of associations between meanings and words. In this work, we consider only a finite set of meanings M and a finite set of words W . In this context, vocabularies can be represented as associations matrices, where each row corresponds to a meaning, and each column to a word. This representation has been extensively used in related work [START_REF] Oliphant | Learning and the emergence of coordinated communication[END_REF][START_REF] Steels | Stochasticity as a source of innovation in language games[END_REF][START_REF] Ke | Self-organization and selection in the emergence of vocabulary[END_REF]. Two parts of the lexicon are distinguished, the coding or production part, which maps a meaning to a set of words weighted by probabilities of usage, and a decoding or interpretation part, mapping a word to a set of meanings that can be interpretated from this word, also weighted by probabilities.

We represent the vocabulary of an agent A as a matrix V (A) of size M ×W , with values of 1 for each word-meaning association used by the agent. Each agent starts with an empty vocabulary, a matrix filled with zeros. The coding matrix V c (A) and decoding matrix V d (A) are derived from V (A) by normalizing respectively over rows and columns:

V c (A) mw = V (A) mw ∑ w V (A) mw V d (A) mw = V (A) mw ∑ m V (A) m w (1)
Normalization factors are used only if V (A) mw = 0. In practice, when coding a meaning m, a word w i is sampled using the distribution (V c (A) mw ) w∈W . When decoding a word w, a meaning m j is interpreted, sampled from the distribution V d (A) mw m∈M . In our case, these distributions are uniform either on the set of words associated to m for coding, or on the set of meanings associated to w for decoding. Those 2 sets change over time, during the vocabulary update.

Vocabulary Update Policy At the end of each interaction, each agent takes into account the result of the interaction by modifying its lexicon. There exists various policies that have been described and studied in previous work [START_REF] Wellens | Adaptive Strategies in the Emergence of Lexical Systems[END_REF]. We are using the one called Minimal Naming Game.

In this policy, updates work this way: when the communication fails, both agents add the used word-meaning association (meaning used as a topic by the speaker, and word uttered by the speaker) to their lexicon, and do nothing if they already had it. If the communication is successful, not only do they add this association to their respective lexicons, they also remove any conflicting synonyms and homonyms. See fig. 3 for an illustration of the update policy in both cases.

Typically, among existing policies, Minimal NG and another one called Basic Lateral Inhibition are used: they are more realistic as they allow synonymy/homonymy and yield faster agreement. Moreover, Minimal NG has been shown to yield similar dynamics as Basic Lateral Inhibition, yet being simple and not depending on any parameter, while the latter depends on 3. This is the reason why we are using the Minimal NG as vocabulary update policy.

Figure 3: Vocabulary update (Minimal NG). Failure (when the hearer interpreted the word as another meaning than the topic): the word-meaning association used by the speaker is added to the hearer's vocabulary. The speaker adds it as well if it was just invented. Success (when the hearer interpreted correctly the word as the topic): both hearer and speaker remove synonyms/homonyms in conflict with the word-meaning association used during the interaction. In both examples the topic is the apple, and the word rimi.

Measures

The self-organization process happening while simulating the Naming Game has complex dynamics, and goes through various states before reaching global consensus. We talk about those dynamics as a convergence process, towards the state where all agents share the exact same lexicon, with exactly one word for each meaning without synonymy and homonymy. This state is stable, lexicons will not change anymore whatever are the modalities of the interaction -which agent is the speaker, which is the hearer, and which meanings and words are used. Convergence and stability for different types of Naming Games has been proved analytically [START_REF] De | The evolution of conventions in multiagent systems[END_REF]. In this paper, we do not focus on whether the model converges or not, but on the speed and complexity properties of the dynamics before convergence. Measures for each of those aspects, used to describe the system while in this intermediate state, were defined in previous work [START_REF] Loreto | Statistical physics of language dynamics[END_REF]. We distinguish local measures -accessible to each agent-from global measures, computed on the whole population.

TCS: Theoretical Communicative Success The Theoretical Communicative Success is a measure of distance to the fully converged state. First, for each meaning, we can consider the probability of having a successful communication when using this meaning as a topic, given a state of the population. The TCS is the average of those probabilities, over all possible meanings. In the case of Random Topic Choice, this measure coincides with the general probability of having a successful interaction. By definition, it is a global measure, not accessible to individual agents. To retrieve its value, we can either estimate it using a snapshot of the population and a Monte Carlo method with random topic choice, or compute it. To detail the exact computation formula, we need to first define the probability of success between two given vocabularies of agents A and B. As detailed in the previous section, a vocabulary has 2 components: a coding part, used to find words associated to a meaning, and a decoding part, used to find meanings associated to a word. For vocabulary V (A), we would then have the 2 matrices V c (A) and V d (A). If A is the speaker and B the hearer, A is coding and B decoding, hence the formula of the probability of success in this case, averaged over all possible meanings:

TCS s (A, B) = 1 M ∑ m ∑ w V c (A) m,w •V d (B) m,w (2)
Because before an interaction we do not necessarily know which agent will be the speaker and which will be the hearer, the 2 situations (A speaker and B hearer / B speaker and A hearer) are to be considered, as equiprobable. The final value TCS(A, B) is the mean of TCS s (A, B) and TCS s (B, A).

To scale up to population level, one can compute an average vocabulary for the whole population V (P), and then the probability of success for an interaction between this lexicon and itself. For a large enough population, this value is indeed a good approximation of the probability of success. V (P) is an element-wise average of the lexicon matrices of all agents.

When using random topic choice, this value abruptly goes from 0 to 1 after a certain number of interactions. These dynamics can be seen on fig. 4, where the random topic choice is represented -among active strategies that are explained in a following section. In practice, we use Monte Carlo estimation for the values at population level over time, and the exact computation for the active topic choice strategy (see following section), as it requires more precision and the population vocabulary is already built.

Local Complexity

The starting state of an agent's vocabulary is empty (all-zero matrices), and the end state is identical coding and decoding matrices, with exactly one distinct word per meaning. But between those 2 situations, through which states goes the vocabulary? How much conflictual information (synonymy and homonymy) has to be considered?

For each agent, we can define a local complexity measure, by counting the number of distinct associations present in the vocabulary. In our case, this is exactly the sum of all elements of the matrix V (A). At the beginning of a simulation, while the vocabulary is empty, this measure equals 0. At the end, its value is the number of meanings M. When using random topic choice, there is a fast growth to a maximum, before a slow decrease to the final value M (can be seen in fig. 4). This measure is nearly proportional to the minimal memory needed to represent the lexicon (as a sparse matrix or a list of word-meaning associations), and therefore should remain low in a cognitively plausible situation.

Active Topic Choice Strategy

The main contribution of this work is the definition of an active strategy for the choice of the topic in each interaction, by comparison to the usual choice of picking meanings randomly (with a uniform distribution over the space of meanings). The strategy has to be local, i.e. use only information available to the agent, namely its own vocabulary and results of past interactions it was involved in.

To both converge quickly and control complexity, behavior should be driven by maximization at each interaction of the Theoretical Communicative Success. However, this value is a global measure, therefore not accessible at agent level. Agents only sample information about the global state of the population, or the average vocabulary V (P), through their interactions as hearer or speaker.

The strategies for active topic choice found in previous work are separated in two levels of decision [START_REF] Schueller | Active control of complexity growth in naming games: Hearer's choice[END_REF]. First, a decision between exploring a new meaning (that is associated to no words in the vocabulary so far) and choosing (exploiting) a meaning among those already used before. Then, if exploiting, deciding which known meaning to use depending on past interaction results.

The strategy introduced in this work keeps those two levels, while basing both decisions on a new measure called Local Approximated Probability of Success using a local representation of V (P).

LAPS: Local Approximated Probability of Success

Here, we define an approximation of V (P), V (P), using information sampled by agents during their interactions. We construct independently the coding and decoding parts V c (P) and V d (P). For every meaning m (and every word w), we use a sliding window over the recent past interactions -of maximal length τ, the time scale parameter-and count the number of times it is associated to each word w (or meaning m ). This value divided by τ is the local estimation of the probability of an other agent coding m using w (or decoding w as m ). With this, we retrieve the values of both matrices V c (P) and V d (P).

Let M c (m) be the memory of the past interactions where m was the topic, if there has been T m such interactions. w t denotes the word used during the t th interaction of the agent using the meaning m. We can now build V c (P):

M c (m) = (w t ) 1≤t≤T m V c (P) mw = T m ∑ t=T m -τ+1 δ w,w t τ (3)
Similarly, by defining M d (w) be the memory of the past interactions where w was the topic, with T w such interactions, we can build V d (P):

M d (w) = (m t ) 1≤t≤T w V d (P) mw = T w ∑ t=T w -τ+1 δ m,m t τ (4)
Until τ interactions have been done with a given meaning or word, ∑ w V c (P) mw and ∑ m V d (P) mw do not sum to 1. The re-maining probability weight is assumed to be associated with failure. If we would normalize to 1, with a single interaction an agent would already estimate as 100% sure that the same word-meaning association would be used again with the same topic for example. Without the normalization, this happens only after τ interactions. In other words, this reflects lack of information due to small sample size. We define a Local Approximated Probability of Success, or a local equivalent of the Theoretical Communicative Success for an agent A with vocabulary V (A):

LAPS A = TCS(V (A), V (P) A ) (5)
For some vocabulary update policies called lateral inhibition, similar matrices are computed, but used directly as an agent's own representation of the lexicon. This usage does not prevent the complexity burst [START_REF] Wellens | Adaptive Strategies in the Emergence of Lexical Systems[END_REF].

Exploration vs. Exploitation

The first choice of our new strategy is between exploring new meanings or exploiting already known ones. Exploration should happen when agents are confident enough about their agreement with the rest of the population over their known meanings [START_REF] Schueller | Active control of complexity growth in naming games: Hearer's choice[END_REF]. The LAPS measure in itself is a measure of confidence, and the simplest way to take this into account is to only explore when reaching the maximum value K M where K is the number of known meanings and M the total number of meanings in the world. This value can actually be reached, thanks to the sliding window of parameter τ.

Multi-Armed Bandit

The second decision process concerns the exploitation part, when picking the topic among the known meanings. We designed a behavior driven by the increase of the LAPS measure. In other words, agents seek the meaning that would yield the greatest increase of LAPS. However, computing the expectancy of this value ∆LAPS is hard computationally speaking, and therefore not suitable for a model of a cognitive process. We can only consider the process a black box, where following a decision between a finite set of options, a reward value is obtained. This is exactly the definition of the Multi-Armed Bandit problem, associated to a class of reinforcement learning algorithms that have been extensively studied [START_REF] Bubeck | Regret analysis of stochastic and nonstochastic multi-armed bandit problems[END_REF]. The name comes from an analogy with a person trying to maximize their gain while facing a set of slot-machines (also called one armed bandit), and being able to use only one at a time. The probability distribution of the reward of each machine is unknown, and the player has to both collect information by playing and exploit the highest rewarding machine -with limited knowledge of its reward distribution -hence keep balance between exploration and exploitation. In our problem, we can see known meanings as the possible arms, and the reward ∆LAPS. Our case is quite specific, as: 1) distributions are non stationary, 2) they depend on past choices, 3) and the number of arms grows over time (and starts at 0). This specific situation led us to choose an algorithm, where weights associated to each arm undergo a decay over time, which let them stay at the same order of magnitude of the initial weights of new arms [START_REF] Clement | Multi-Armed Bandits for Intelligent Tutoring Systems[END_REF]. Our algorithm depends on 2 parameters: integrated balance between reward-driven exploitation and random exploration between arms through the parameter γ, and time scale n for the decay of weights. As a reward, we consider the increase of LAPS yielded by the interaction, ∆LAPS, or 0 if the latter is negative in order to avoid negative weights. See algorithm 1.

Results

For all simulations, we set N=M=W =40, compute up to 80,000 interactions and take the mean over 8 trials. The situation M=W is the most constrained and complex to solve, as synonymy and homonymy are more probable. We ran simulations for 1 ≤ τ ≤ 50, and set n=τ. For the exploration rate, if the condition 0 < γ 1 is respected, the actual value of γ does not matter much, as its only function is to avoid rare cases where some weights reach a value of 0 and cannot be selected anymore. We set γ=0.01. However, we also ran simulations with pure random choices at the bandit level, to be able to study the influence of each level of our algorithm. This case identifies as γ=1.

The evolution of the TCS and complexity over time is represented on fig. 4, for several values of the time scale τ. They are compared on the same plots with Random Topic Choice. We can see that convergence is faster for low values of τ, the fastest being for τ=2, which is 4 times faster than Random Topic Choice. As for complexity, for all configurations ex- cepted τ=1 values stay below the final level 40. After reaching a first threshold, they increase linearly with time, the slope being smaller for higher values of τ. For τ=1, the maximum value is only half of the maximum reached by Random Topic Choice. It is understandable that τ=1 is an outlier: in this case LAPS is an autocorrelation with the current interaction, by definition older interactions are not taken into account.

On fig. 5, we can see the dependency of convergence time on the parameter τ, plotted for configurations γ=0.01, γ=1 and the value for Random Topic Choice as a reference. Both have dynamics consistently faster than Random Topic Choice for low values of τ, however γ=0.01 performs better. Excepted for τ=1, convergence time increases linearly with τ for both, with a minimum at τ=2, and a smaller slope for γ=0.01.

Discussion

Results show that the new strategy presented in this paper 1) allows fast convergence, 2) controls efficiently complexity growth, 3) its dynamics are consistent and highly correlated with parameter change, 4) the 2 levels of the algorithm each contribute to the increased performance. With τ=2, each agent on average only speaks 15 times about each meaning before convergence (i.e. less than half the population), and information has already been both conveyed between all agents and disambiguated. The linearity of the evolution of TCS and complexity lets think that this algorithm may as well scale efficiently to other values of N, M and W . Compared to previous work, this topic choice algorithm is more robust, and optimal parameters are easier to find. It generalized well to Minimal Naming Game and can be used for all other Naming Game models.

LAPS is coherent from a cognitive point of view, and corresponds to an actual internal confidence about quality of communication with the rest of the population. As stated in the results section, the case τ=1 is an outlier, being a simple autocorrelation with the current interaction. The optimal value τ=2 is then the lowest possible value taking into account past interactions, i.e. V (P) takes the lowest possible memory, which is therefore credible for humans. Further work will be needed to determine for which values of N, M and W τ=2 stays the optimal value.

Figure 2 :

 2 Figure 2: Interaction process: Beforehand, 2 individuals have been randomly selected among a population, an designated as speaker (S) and hearer (H). 1. S chooses a topic, 2. S checks its vocabulary to find or invent an associated word, 3. S utters the word, 4. H guesses the intended meaning, 5. S indicates the intended meaning.

Algorithm 1

 1 LAPSmax Bandit Multi-Armed Bandit algorithm used as a topic choice strategy maximizing the LAPS measure. New arms are created with weights w a equal to the reward r i obtained at the end of an interaction with a new meaning. Require: γ rate of exploration for bandit Require: n scale for weights decay Require: vocabularies V and V , #meanings M K ← |V.known meanings()| if LAPS A = K M then m ← random (V.unknown meanings()) else for a ∈ Arms do wa = w a ∑ j w j p a = (1γ) • wa + γ K end for Sample m ∈ Arms using distribution (p a ) a∈Arms end if return m {Interact using topic m and compute reward r} if m ∈ Arms then w m ← n n+1 • w m + r else Add m to Arms with w m = r end if

Figure 4 :

 4 Figure 4: Theoretical success rate over time and Local Complexity for Random Topic Choice and Active Topic Choice with several values of the time scale τ used in the LAPS measure. N=M=W =40, γ=0.01, mean over 8 trials.

Figure 5 :

 5 Figure 5: Convergence time depending on time scale τ used in LAPS measure, for both γ=0.01 and γ=1, compared with Random Topic Choice. N=M=W =40, mean over 8 trials.

Acknowledgments

The IdEx program (Univ. de Bordeaux) allowed W. Schueller to visit V. Loreto. We thank Miguel Ibañez Berganza and Benjamin Clément for the fruitful discussions.

Source code

The Python code used for the simulations of this paper is available as open source software: https://github.com/wschuell/notebooks cogsci2018


Farouk Ghallabi

Fawzi Nashashibi
email: fawzi.nashashibi@inria.fr

Ghayath El-Haj-Shhade

Marie-Anne Mittet
email: marie-anne.n.mittet@renault.com

LIDAR-Based Lane Marking Detection For Vehicle Positioning in an HD Map

. Moreover, map-based localization using HD maps became an interesting source of information for intelligent vehicles. In this paper, a Map-based localization using a multilayer LIDAR is proposed. Our method mainly relies on road lane markings and an HD map to achieve lane-level accuracy. At first, road points are segmented by analysing the geometric structure of each returned layer points. Secondly, thanks to LIDAR reflectivity data, road marking points are projected onto a 2D image and then detected using Hough Transform. Detected lane markings are then matched to our HD map using Particle Filter (PF) framework. Experiments are conducted on a Highway-like test track using GPS/INS with RTK correction as ground truth. Our method is capable of providing a lane-level localization with a 22 cm cross-track accuracy.

I. INTRODUCTION

Precise vehicle localization is an important task for autonomous driving and advanced driver assistance systems (ADAS). Indeed, precisely estimating the position of the ego-vehicle with respect to the driving lane is crucial for many ADAS functionalities such as Lane Keeping Assist (LKA). This localization level has been introduced in [START_REF] Edmap | Enhanced Digital Mapping Project Final Report[END_REF] as «WHEREINLANE» application. The report [START_REF] Edmap | Enhanced Digital Mapping Project Final Report[END_REF] shows that such an application requires a localization accuracy of 50cm which is not affordable by the majority of current localization systems. Existing localization approaches can be divided into two categories: global and relative localization [START_REF] Thrun | Robust Monte Carlo Localization for Mobile Robots[END_REF]. Relative localization techniques are mainly focusing on pure odometry such as visual odometry [START_REF] Nistér | Visual odometry[END_REF], laser odometry [START_REF] Zhang | LOAM : Lidar Odometry and Mapping in Real-time[END_REF] or SLAM-Based approaches [START_REF] Durrant-Whyte | Simultaneous localization and mapping: Part I[END_REF], [START_REF] Bailey | Simultaneous localization and mapping (SLAM): Part II[END_REF] such as Visualbased SLAM [START_REF] Fuentes-Pacheco | Visual simultaneous localization and mapping: a survey[END_REF] and Laser-based SLAM [START_REF] Moosmann | Velodyne SLAM[END_REF], [START_REF] Nuchter | 6D SLAM-3D Mapping Outdoor Environments[END_REF]. Global localization is mostly determined by global navigation satellite systems (GNSS). However, the position accuracy of automotive GNSS solutions is no better than 2-3 m in opensky environments [START_REF] Zogg | GPS: Essentials of Satellite Navigation: Compendium: Theorie and Principles of Satellite Navigation[END_REF], which is considered to be insufficient for autonomous driving. Alternatively, localization on a prior map, usually known as Map-Based Localization or Mapmatching, has gained interest since maps can be pre-built very accurately. Current Map-based localization can be divided according to the type of the sensor used: passive or active. Vision-based (passive sensor) techniques [START_REF] Buczko | Efficient Global Localization Using Vision and Digital Offline Map[END_REF], [START_REF] Zhu | Real-time global localization with a pre-built visual landmark database[END_REF], [START_REF] Chausse | Vehicle localization on a digital map using particles filtering[END_REF], [START_REF] Košecká | Global localization and relative positioning based on scale-invariant keypoints[END_REF], [START_REF] Wu | Vehicle localization using road markings[END_REF] match visual features or landmarks to map attributes. Visual landmarks can be of different types: lane sides [START_REF] Chausse | Vehicle localization on a digital map using particles filtering[END_REF], traffic signs painted on the road such as arrows, pedestrian crossing and speed limits [START_REF] Wu | Vehicle localization using road markings[END_REF], traffic signs (vertical poles) [START_REF] Li | Localization for intelligent vehicle by fusing mono-camera, low-cost GPS and map data[END_REF], or feature points such as SIFT [START_REF] Košecká | Global localization and relative positioning based on scale-invariant keypoints[END_REF]. Despite promising localization results, the main restrictions of vision techniques are the dependency on external light conditions and the sensitivity to shadows and illumination noise. On the other hand, LIDAR-based approaches (active sensor) have proven to be more accurate than vision techniques. Indeed, LIDAR sensors are much less sensitive to light conditions and provide a 3D representation of the environment with centimetric accuracy. However, processing the whole LIDAR points is time consuming and is not suitable for real time applications. Therefore, 3D point cloud is usually projected into 2D representations such as 2D orthographic reflectivity map [START_REF] Levinson | Robust vehicle localization in urban environments using probabilistic maps[END_REF], precise height map [START_REF] Fu | An efficient scan-to-map matching approach for autonomous driving[END_REF], [START_REF] Pollard | Step and curb detection for autonomous vehicles with an algebraic derivative-based approach applied on laser rangefinder data[END_REF], or a combination with other information such as colors, curvatures and normals [START_REF] Chong | Mapping with synthetic 2D LIDAR in 3D urban environment[END_REF]. Similarly to vision, recent research has been conducted on lidar-based road marking detection [START_REF] Hata | Road marking detection using LIDAR reflective intensity data and its application to vehicle localization[END_REF], [START_REF] Kammel | Lidar-based lane marker detection and mapping[END_REF], [START_REF] He | Lane marking detection based on convolution neural network from point clouds[END_REF]. Road markings are reflective objects and can be detected using intensity laser data. Since road markings occur on the road surface, it is first necessary to perform road segmentation. Road segmentation methods often rely on detecting curbs such as road-edges. Then, marking points are extracted from the road plane according to their high reflectivity values. In [START_REF] Zhang | LIDAR-based road and road-edge detection[END_REF], the road and road-edge points are extracted using the elevation information of the LIDAR, and then validated on the ground plane for robustness. In [START_REF] Kammel | Lidar-based lane marker detection and mapping[END_REF], an occupancy grid map of elevation difference is computed. Grid cells whose elevation difference is smaller than a threshold are clustered as ground cells. The work in [START_REF] Hata | Robust Curb Detection and Vehicle Localization in Urban Environments[END_REF] is an extension to [START_REF] Montemerlo | Junior: The stanford entry in the urban challenge[END_REF] where the authors developed a curb detection method based on ring compression analysis. At first, LIDAR point cloud are structured in a polar grid. Then, distances between rings are computed and compared to the expected distance when rings intercept a flat plane. If the computed distance is below a threshold value then the grid cell is classified as curb.

In this paper, we propose a lane markings detection method for localization within an HD map using a Multilayer LIDAR. At first, road segmentation is performed by analysing the geometry of each returned layer points. Then, high reflective road points are projected onto a 2D grid which is always centred around the vehicle position and with a cell size of 0.15m ×0.15m. A line detection using Hough Transform [START_REF] Duda | Use of the Hough transformation to detect lines and curves in pictures[END_REF] is implemented on the resulted image. To validate our approach within a localization system, a map-matching algorithm using particle filter and lane markings as features is implemented. Experiments were conducted on a highwaylike test track at high vehicle speeds ≈ 80 kph. The laser range finder used in our experiments is a Velodyne VLP-32C. Moreover, a highly accurate GPS with RTK correction signals (ixblue: ATLANS-C) was used as ground truth for evaluation purposes. Finally, a comparative analysis between our localization system and an automotive u-blox GNSS is provided based on a cross-track metric [START_REF] Héry | Map-based Curvilinear Coordinates for Autonomous Vehicles[END_REF].

The remainder of this paper is structured as follows: Section II describes the proposed road segmentation method. Section III illustrates the developed map-matching algorithm with a brief description of the used HD map. An experimental evaluation is presented in section IV and we conclude with perspectives and future work in section V.

II. ROAD SEGMENTATION

In order to detect lane markings, it is necessary at first to segment road points from the whole point cloud since road markings are located on the road surface. Once the road is segmented, an intensity-based thresholding is applied to extract road marking points. Our segmentation method is based on a geometric analysis of the layer scan. A layer scan is a set of captured layer points from a full physical rotation of the LIDAR. Under the assumptions of flatness and smoothness of highway road surfaces, the geometric pattern of a layer scan intercepting the road can be locally modelled by a circular arc. However, due to non smooth and discontinuous surfaces outside road boundaries such as grass, bushes and gravel, the geometry of some regions of the scan is scattered (figure 1.a) and cannot be modelled. Differently to curb-based approaches, our method is capable of detecting change of surface type (from smooth to nonsmooth and vice versa) even in the absence of distinctive geometric curbs. For illustration purposes and with no loss of generality, we consider the case of one layer scan, also called ring. In addition, we consider the polar representation of points (r, θ, ϕ), were r is the range of the returned point, θ is the layer vertical angle and ϕ is the azimuth (horizontal) angle.

A. Ring geometry analysis

Let L = {p i = (r i , θ, ϕ i )} be the set of ring points parametrized in polar coordinates. So far, L is discretized to a set of contiguous slices s j (figure 1.b) that are defined by the following system:

     s j = {p k = (r k , θ, ϕ k ), ϕ j min ≤ ϕ k ≤ ϕ j max } ϕ j max = j × ϕ res , ϕ j min = (j -1) × ϕ res (1)
where, ϕ res is the slice angular resolution, j ∈ {1, 2, .., M } and 2π = ϕ res × M . As previously described, if a ring slice is intercepting the road then its geometry can be approximated by a circular arc. However, a slice outside the road is generally intercepting discontinuous, non flat and non smooth surfaces such as vegetations and therefore is not geometrically consistent with a road slice model. Consequently, the proposed segmentation method is based on fitting a circular arc to each slice s j , then checking the residual fitting score to determine whether s j is a circular arc or not. The computation of the residual fitting score goes through three steps:

• Calculation of the mean altitude µ z of slice points. This represents the horizontal z-plan to which the expected circular arc sj should belong.

• Computation of the center and the radius (Cj, Rj) of sj as follows: Cj = (0, 0, µz), Rj = µz × cot(θ) • Compute the normalized sum of slice point projections to sj which is also called the fitting score value. The projection function for a given point

p k = (r k , θ, ϕ k ) to sj = (Cj, Rj) is: proj(p k , sj) = (Rj -r k ) 2 + (r k × sin(θ) -µz) 2

B. Slice-edge detection

The proposed road segmentation method makes use of the previously computed fitting score values. At first, the 2D (x, y) horizontal plan is partitioned into four quadrants Q 1,2,3,4 as depicted in figure (2.a). This partition allows us to simplify the segmentation process and to parallelize the computation. Indeed, by definition, the LIDAR ring crosses road boundaries at most once within each quadrant. Therefore, the problem can be formulated as executing a quadrant search for the slice that defines the road boundary, which we call in this paper a slice-edge.

To do that, an iterative process is implemented for each quadrant. Starting from an initial position on the ring section (blue dots in figure 2.a), the process iterates over the slices until crossing the first slice (red dots in figure 2.a) whose score value is greater than a threshold T 1 . The azimuth angle that corresponds to the detected slice-edge is called: cut angle ψ. As a result, slices coming after ψ are clustered as non road slices. However, in practical situations, some slices can have high scores even when belonging to the road surface. This is due in general to a local point cloud distortion caused by vehicle dynamics or road surface irregularities. A Gaussian filter is applied to each quadrant in order to smooth the data. For example, in figure (2.b) and with a fixed threshold of 0.05m, the cut angle (red) in the original data is approximately 42.5°whereas in the smoothed data (black) it is approximately 50°.

D. Detection of lane markings

Once road points have been segmented, the intensity LIDAR data is used to extract marking points. The extraction is based on the application of an intensity threshold T 2 as we claim that most of high reflectivity returns from road surface are stemming from road markings. A 2D grid, centred around the ego-vehicle, is constructed with a cell size equals to 0.15m. The cell size is set taking into account road marking widths on highways. The maximum intensity value of points that fall within a grid cell is stored. Grid reflectivity data is then transformed to a 2D intensity image to which a line detection using a Hough Transform [START_REF] Duda | Use of the Hough transformation to detect lines and curves in pictures[END_REF] is applied for lane markings detection. We chose the polar representation (r, θ) of lines as it is relevant for localization. Indeed, r and θ are interpreted as the vehicle lateral distance and heading to the lane marking. For robustness, a set of constraints are used to filter out the false detections. These constraints are summarized in the following:

1) Lines on the driving direction: the detection process only searches for lines that are approximately parallel to the driving direction of the ego-vehicle. This can be achieved by defining constraints on the θ search of the Hough Transform.

2) Parallel lines: Consider M the number of detected lines and N k the number of parallel lines to a given line l k . The filtering algorithm picks the line given by:

l max = arg max k∈{1,..,M } N k
and discards all the lines that are not parallel to it.

3) Line fusion: : For two detected lane markings l 1 = (r 1 , θ 1 ) and

l 2 = (r 2 , θ 2 ), if r 1 -r 2 < then a new line is created l new = a × l 1 + b × l 2 where a + b = 1. For our case a = b = 1 2 .

E. Tracking of lane markings

A tracking algorithm using Kalman filter [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF] and inertial measurements has been implemented to track lane markings. A standard Kalman filter has proven to be sufficient for our case as we used a simplified prediction model given by the following:

r k = r k-1 -v × dt × sin(w × dt) + ν r,k θ k = θ k-1 -w × dt + ν θ,k (2)
where ν r,k = N (0, σr) and ν θ,k = N (0, σ θ ) are gaussian distributions with zero mean and variances (σ r , σ θ ) and v, w are the vehicle speed and yaw rate, respectively. The model validity is guaranteed under the assumption of low θ values. This assumption holds for highway use cases where vehicles encounter small changes of heading while driving. Figures

III. MAP-BASED LANE-LEVEL LOCALIZATION

The detected lane markings are used as input features to the map matching algorithm. To achieve precise vehicle localization, a High Definition (HD) Map, provided by a mapmaker, with a lane-level information and absolute accuracy of 2 cm has been used. The prototype map is composed of many links that represent lane markings. Each link is a sequence of segments expressed in WGS 84 geodetic system.

A. Map-matching

In order to localize the vehicle on the previously described map, a map matching algorithm has been implemented under the Particle Filter (PF) framework [START_REF] Dellaert | Monte carlo localization for mobile robots[END_REF]. Particle filtering is able to represent multi-modal non-Gaussian distributions which is convenient for global localization. A set of particles (Xj, Wj) that represent vehicle pose hypotheses and their relative weights are defined. The vehicle pose is considered in our study as in 2D space: (x j , y j , γ j ). The particle filter algorithm is implemented in three main steps: prediction, update and resmapling steps.

For the prediction step, each particle hypothesis at time k is estimated from the previous estimation at time k -1 and from odometry data (velocity v and yaw rate w). Different vehicle dynamic models can be found in the literature [START_REF] Schubert | Empirical evaluation of vehicular models for ego motion estimation[END_REF]. In our case, we considered a constant velocity and yaw rate model given by the following system (3):

     γ k = γ k-1 + w × dt + νw x k = x k-1 + v × dt × cos(γ k ) + νx y k = y k-1 + v × dt × sin(γ k ) + νy (3)
where νw = N (0, σw), νx = N (0, σx), νy = N (0, σy) are gaussian noise distributions. Moreover, GNSS signals are used to initialize the filter.

For the update step, particle weights are recursively updated from the detected lane markings and the prototype map M . For each particle Xj, the nearest segment for each map link (lane marking) on the road section is searched. Then, its relative distance and heading with respect to Xj are computed. We will note these parameters as m j l = (r j,l , θ j,l ), where l is the nearest segment for a given map link. The update step is composed of two parts:

1) Weight calculation from the Map: If the projection of a particle to the map is not within the driving road section, then a zero weight is assigned.

2) Weight calculation from lane markings: for a given observation data z k = (r k , θ k ) and a matched link segment m j l = (r j,l , θ j,l ) a likelihood function similar to [START_REF] Chausse | Vehicle localization on a digital map using particles filtering[END_REF] is computed as follows:

p(m j l |z k ) = e (r k -r j,l ) 2 2σ 2 r + e (θ k -θ j,l ) 2 2σ 2 θ
Where, σ r and σ θ define the measurement noise variances. Taking into account all the observations and all the matched links, the calculated weight is given by:

W j = α k l p(m j l |z k )
Where α is a weight normalization factor.

For the resampling step a systematic resampling strategy with a threshold of 50% effective particles has been implemented to carry on particle degeneracy,. Indeed, according to [START_REF] Douc | Comparison of resampling schemes for particle filtering[END_REF], systematic resampling is found to provide comparable results to other state of the art methods. However, it is often preferred thanks to the simplicity of implementation. II) to evaluate the robustness of our approach. The map matching algorithm has been developed under the ROS framework.

IV. EXPERIMENTAL RESULTS

B. Road marking detection

For the road segmentation part, a histogram-based analysis was conducted to set the threshold T 1 (figure 6). The histogram illustrates the score values for Q 1 -slices. One can notice that the histogram peak corresponds to slices for which the score values are less than 0.1m, which is expected as most of the slices lie on the road surface. Consequently, the segmentation threshold value T 1 is set to 0.1m.

C. Cross-track error results

In this paper, we adopted a cross-track metric presented in [START_REF] Héry | Map-based Curvilinear Coordinates for Autonomous Vehicles[END_REF] to evaluate our localization system. To the best of our knowledge, most of state of the art localization systems are evaluated in terms of absolute error. In our approach, we propose an evaluation based on a cross-track metric since it is more relevant for localization within an HD map.

As previously stated, the map is composed of many links (lane markings) in polylines form. For each position: p gt (ground truth), p est (our estimation) and p ubx (u-blox) , the shortest signed distances to map links of the driving road section are computed, let them be (l i 1 , l i 2 , ..l i N ), where the superscript i denotes one of the three positions. The evaluation of the cross-track error for our estimate and the u-blox position is given by the following: Finally, a velocity dependent error evaluation is given in Table II. This is particularly interesting to evaluate the effect of point cloud distortion of spinning LIDARs in motion.

ep = 1 N N i=1 (l gt i -l p i ) p = { pest, p ubx }
Results of Table II demonstrate that we have similar performances for low and high speed scenarios. This justifies that our road segmentation and lane marking detection algorithms are robust to point cloud distortion by cause of the vehicle movement. More precisely, since the main movement of the vehicle is in the x-direction (small displacement in y- In this paper, a lane-level localization within an HD map is presented. A LIDAR sensor has been used to detect lane markings and to match them with a prototype map. The lane marking detection is implemented in two different steps: road segmentation and Hough transform on an intensity image. Road segmentation is based on Ring geometric analysis which measures the discrepancy of real local geometric features with expected features on the road surface. This is particularly useful to remove environment regions where surfaces are not smooth and are discontinuous (e.g vegetations, gravel). Furthermore, a Hough line transform using an a priori information on the environment is applied to detect lane markings. Finally, a map-matching algorithm has been implemented to validate the detection phase. Current results are promising and sufficient for Highway use-cases.

As a perspective, we believe that the road segmentation method can be advantageous to many other applications that are not covered by this study. For example, a 3D marking map can be directly built from the 3D marking points.

Future work will deal with the longitudinal localization (the along-track error) by using a monocular camera in addition to the LIDAR to detect more Highway landmarks such as traffic signs and poles. Indeed, in contrast to urban environments, highways lack reliable landmarks for longitudinal localization. Thus, we expect that minimizing the alongtrack error to be more challenging.

Fig. 1 .

 1 Fig. 1. Slice representation of a layer scan

Fig. 2 .

 2 Fig. 2. Detection of slice edges

Fig. 3 .

 3 Fig. 3. Segmentation & registration of LIDAR data

  (4.a) and (4.b) depict the projection of the 3D road marking points to a 2D image, detected lane markings are shown in figure (4.c).

Fig. 4 .

 4 Fig. 4. Projection of 3D points to 2D image and lane markings detection

Fig. 5 .

 5 Fig. 5. Renault test track (google earth map)

Fig. 6 .

 6 Fig. 6. Q1-Slice variances histogram Once 3D road points are segmented, LIDAR reflectivity byte values (from 0 to 255) are used to segment road marking points. The reflectivity threshold used in our approach is fixed with respect to the sensor datasheet which assigns a value greater than T 2 = 100 for retroreflector objects such as road markings. Road marking points are projected onto a 2D grid, centered around the vehicule position, of size [0 m, 20 m] × [-15 m, 15 m] and a cell resolution of 0.15 m × 0.15 m. Hough transform is only applied for small θ values (from π 10 to π 10 ) to insure detection of nearly parallel lines. The threshold distance for line fusion is set to 0.5m.

Figures 7 .

 7 Figures 7.a and 7.b depict histograms and time variations of the cross-track error for the proposed method and the ublox system for a vehicle speed = 80 kph. Means and standard deviations are also summarized in TableI. One can notice that, with the use of the laser-based lane markings and inertial data, the mean error is reduced from -1.42 m to 0.04 m and the standard deviation from 1.0 m to 0.22 m.

Fig. 7 .

 7 Fig. 7. Histograms (a) and timeplot (b) of the cross-track errors

TABLE I CROSS

 I

-TRACK ERROR: MEAN AND STD


Xun Wang

Mohamed S Ghidaoui

Pedro J Lee

Regularization for pipeline impulse response extraction with least square deconvolution

Impulse response function (IRF) is frequently used for pipeline defect (e.g., leakage and blockage) detection. IRF can be extracted via a least square (LS) deconvolution method if the input signal is known. However, the LS deconvolution is usually an ill-posed inverse problem, thus the regularization approach is needed to solve this problem. The determination of the regularization parameter is essential for regularization problems and is strongly problem-dependent (no universal method that always produces robust and good results). In the present paper, prevailing regularization methods Generalized Cross Validation (GCV) and L-curve are used and tested to decide the optimal regularization parameter. Pseudo random binary sequence (PRBS) signal, which switches between two constant levels with noise-level amplitude, is used as input signal. Numerical and experimental examples show the importance of the determination of regularization parameter and justify that the proposed methods are efficient for extracting the IRF with PRBS signal.

1

INTRODUCTION

Pipeline defect detection is an important issue in urban water supply systems and have attracted a lot of attention from both theoretical and application aspects [1][2][3][START_REF] Lee | Leak location using the pattern of the frequency response diagram in pipelines: a numerical study[END_REF][START_REF] Covas | Application of hydraulic transients for leak detection in water supply systems[END_REF][START_REF] Meniconi | Anomaly pre-localization in distribution-transmission mains by pump trip: preliminary field tests in the Milan pipe system[END_REF]. The timedomain impulse response function (IRF), which is the output of the system when presented with an impulse input signal, clearly shows the reflections of impulse input signal at defects and is thus frequently used for pipeline defect (e.g., leakage and blockage) detection [START_REF] Nguyen | Least squares deconvolution for leak detection with a pseudo random binary sequence excitation[END_REF][START_REF] Lee | Leak location in pipelines using the impulse response function[END_REF][START_REF] Vitkovsky | Leak and blockage detection in pipelines via an impulse response method. Pumps, electromechanical devices and systems applied to urban water management[END_REF]. IRF also depicts the relationship between the input signal and the resultant output signal and illustrates how the system modifies the input signal as it propagates through the system.

IRF can be extracted by imposing a pseudo random binary sequence (PRBS) signal [START_REF] Newland | An introduction to random vibrations, spectral & wavelet analysis[END_REF][START_REF] Liou | Pipeline leak detection by impulse response extraction[END_REF][START_REF] Lee | Valve design for extracting response functions from hydraulic systems using pseudorandom binary signals[END_REF] that switches between two constant levels. The spectrum of a PRBS signal is approximately flat in a wideband of frequency provided that the sequence length is long enough, meaning that a PRBS signal corresponds to a broadband noise and is thus suitable for low frequency system identification [START_REF] Newland | An introduction to random vibrations, spectral & wavelet analysis[END_REF]. The amplitude of PRBS can be noise-level such that the system operation disturbance and structural fatigue by other wave generation methods (e.g., step or impulse signals by valve operation) can be avoided. For this reason, recent pipeline leakage detection methods proposed in [START_REF] Wang | Spectral-based methods for pipeline leakage detection[END_REF], which is efficient but requires multiple measurements, become easy to implement.

The relationship between designed input signal and measured output signal can be represented by IRF via a convolution, whose discrete form follows a linear model [START_REF] Nguyen | Least squares deconvolution for leak detection with a pseudo random binary sequence excitation[END_REF]. As a matter of fact, linear relationship has been found and used widely in pipeline detection problem [START_REF] Wang | Pipeline leak detection using the matchedfield processing method[END_REF][START_REF] Wang | Identification of multiple leaks in pipeline: Linearized model, maximum likelihood, and super-resolution localization[END_REF][START_REF] Wang | Identification of multiple leaks in pipeline II: Iterative beamforming and leak number estimation[END_REF][START_REF] Zouari | Single extended blockage identification using a model-based matched-field processing approach[END_REF], particularly the discrete form of the linear models have been proposed [START_REF] Keramat | Inverse transient analysis for pipeline leak detection in a noisy environment[END_REF][START_REF] Zhou | Compressive Sensing-Based Multiple-Leak Identification for Smart Water Supply Systems[END_REF]. For example, the pressure head variation due to leaks equals to contribution from different leaks and each contribution equals to multiplication of leak size and a function of leak location [START_REF] Wang | Pipeline leak detection using the matchedfield processing method[END_REF][START_REF] Wang | Identification of multiple leaks in pipeline: Linearized model, maximum likelihood, and super-resolution localization[END_REF][START_REF] Lee | Frequency domain analysis for detecting pipeline leaks[END_REF]. However, in the impulse response extraction problem, the least square (LS) solution of the linear model is ill-posed, thus the regularization problem needs to be considered [START_REF] Nguyen | Least squares deconvolution for leak detection with a pseudo random binary sequence excitation[END_REF].

The determination of the regularization parameter is critical for regularization problems and is strongly problem-dependent. Up to now, there is no universal method for determining the regularization parameter that always produces robust and good results [START_REF] Bauer | Comparing parameter choice methods for regularization of ill-posed problems[END_REF]. In a recent paper [START_REF] Nguyen | Least squares deconvolution for leak detection with a pseudo random binary sequence excitation[END_REF], the Tikhonov regularization for pipeline impulse response extraction with least square deconvolution has been proposed. However, the determination of the regularization parameter in this problem has not been discussed. In the present paper, prevailing regularization methods Generalized Cross Validation (GCV) [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF] and L-curve [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the Lcurve[END_REF][START_REF] Hansen | The use of the L-curve in the regularization of discrete ill-posed problems[END_REF] are used and tested to decide the optimal regularization parameter.

The organization of this paper is as follows. In Section 2, the linear model for IRF extraction is introduced. Section 3 briefly introduces GCV and L-curve methods. Numerical and experimental results are presented in Section 4. Finally, conclusions are drawn in Section 5.

MODEL

A reservoir-pipe-valve system is considered in this paper. A PRBS input signal, denoted as () xt , is used to excite the system. Let () yt denote the system output signal, the input and output signals have the relationship:

       ( ) ( ) ( ) , y t g x t d (1) in which () gt is the IRF of the system. Let  [ ], [ ], [ ],
1, , ,

x n y n g n n N be the discretized ( ), ( ), ( ) x t y t g t , a linear model describes the convolution Eq. (1): ,

y = Xg (2) in which ( [1], , [ ]) , ( [1], , [ ]) TT y y N g g N y = g = (3) and           [1] 0 0 0 [2] [1] 0 0 [3] [2] [1] 0 . [ ] [ 1] [ 2] [1] x xx x x x x N x N x N x X (4)
Here, N is the maximum concerned time-step of IRF. For leak detection problem, in general only the first period of IRF is concerned but both the input signal of PRBS and the measured signal can be much longer.

Given [ ], [ ]

x n y n , the discrete IRF can be obtained by solving the least square (LS) problem with Tikhonov regularization [START_REF] Nguyen | Least squares deconvolution for leak detection with a pseudo random binary sequence excitation[END_REF]:

        22 1 22 ˆarg min ( ) ( ) TT g g y Xg g X X I X y (5)
in which  is the regularization parameter and 2 stands for the L2-norm. Note that according to the Gauss-Markov theorem, the LS solution without regularization (i.e.,   0 ) of Eq. ( 2) is optimal in the sense of minimum variance unbiased estimator if the measurement errors of y are uncorrelated, zero-mean, and have a constant variance. It can be seen from Eq. ( 4) that X is a lower triangular matrix which theoretically has full rank. In practice, however, the matrix X is usually numerically rank deficient and ill-conditioned (in particular for large size of X ), which can be verified via the singular value decomposition (SVD). As a result, the LS solution is ill-posed. Therefore, regularization (i.e., a non-zero  in Eq. ( 5)) has to be applied; in this case, a new problem is solved that does not inherit the properties of the LS solution.

The determination of the regularization parameter  is critical for Tikhonov regularization and is strongly problem-dependent that no universal method always produces robust and good results [START_REF] Bauer | Comparing parameter choice methods for regularization of ill-posed problems[END_REF]. This issue is illustrated by Figure 1. Here, experimental results with data that have been used in Ref. [START_REF] Lee | Valve design for extracting response functions from hydraulic systems using pseudorandom binary signals[END_REF] are shown, where a steel intact pipe with length being 37.5 m is considered and PRBS signal are used to excite the system. The experimental setup, as well as the input PRBS signal and output measurements, can be found in Ref. [START_REF] Lee | Valve design for extracting response functions from hydraulic systems using pseudorandom binary signals[END_REF]. Figure 1 shows the IRF estimation obtained from Eq. ( 5) with various    45 [10 ,10 ] . It is clear from this figure that the IRF extraction is sensitive to  and its determination is important for obtaining an IRF with least noise. To the best of our knowledge, this issue has not been reported in the literature, which is the main goal of the present paper.

GCV AND L-CURVE FOR DETERMINING THE REGULARIZATION PARAMETER

In this section, the GCV [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF] and L-curve [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the Lcurve[END_REF][START_REF] Hansen | The use of the L-curve in the regularization of discrete ill-posed problems[END_REF] methods are briefly introduced to decide the regularization parameter  .

GCV originates from the leave-one-out (LOO) cross-validation method. The basic idea is to minimize the average estimation error of a data point [] yn computed from a (N-1)- dimensional data/model where the n-th line in Eq. (2) left out [START_REF] Allen | The relationship between variable selection and data agumentation and a method for prediction[END_REF]. GCV is a rotationinvariant version of the LOO cross-validation which decides the regularization parameter  by minimizing: where

      2 2 ( ( ) () ( (
   1 ( ) ( ) TT A X X X I X . (7)
L-curve is a graphical tool with a log-log plot of

   2 () g versus    2 ()
y Xg for all the possible  . Note that     (log ( ),log ( )) usually has a distinct L-shape, and the vertical and horizontal parts of the L-curve correspond to overregularization and under-regularization, the corner point of the L-curve is seen as a good estimate of the regularization parameter. This corner point is automatically computed by choosing the maximum point of the curvature  () of the plot     (log ( ),log ( )) .

NUMERICAL AND EXPERIMENTAL RESULTS

In this section, a numerical example is introduced to show the performance of the proposed methods. The pipe length is  1 l km and the pipe diameter is  0.2 D m. The wave speed is  1000 a m/s. Let the coordinates of the upstream and downstream nodes be  0

x and  1000 x m. A leak is assumed to locate at  300 x m and the leak size is

  4 3 10 dL CA m 2 .
A transmitter and a pressure sensor are both located at the downstream node  1000 x m. A PRBS signal is emitted by the transmitter, it propagates along the pipe and is measured by the sensor. Both the signal emission and the measurement last  80 T s. The wave propagation is simulated via the Method of Characteristics (MoC) [START_REF] Chaudhry | Applied Hydraulic Transients[END_REF], where the grid lengths in space and in time are  10

x m and  0.01 t s. Figure 2 SVD is done for the matrix X where  8000 N ; the singular values in descending order normalized by the maximum are shown in Figure 3, which tends to 0 at around 7000-th singular value, meaning that the LS deconvolution is an ill-posed inverse problem. Here, the GCV and L-curve methods are used to regularize this problem. In fact, the IRF can also be obtained from the cross-correlation method [START_REF] Liou | Pipeline leak detection by impulse response extraction[END_REF]:

    1 1 ( ) [ ] [ ].
Nt oi i

R t x i y i t N (8)

Note that in order for Eq. ( 8) to be used as a valid estimate of the IRF, the signal is needed to be broadband, i.e, the spectrum of the input function is flat within the considered frequency range. In this section, numerical results of IRF estimation using four methods are compared, i.e., the cross-correlation method (Eq. ( 8)), LS without regularization (   0 in Eq. ( 5)), GCV, and L-curve. Figure 4 shows the normalized reconstructed IRF. Here, independent and identically distributed Gaussian random noise is added to the measurement and the signal-to-noise ratio (SNR) is 10 dB. Note that the reference signal in the definition of SNR is the average of the output signal, instead of the signal drop due to leak in Refs. [3,[START_REF] Wang | Pipeline leak detection using the matchedfield processing method[END_REF][START_REF] Wang | Identification of multiple leaks in pipeline: Linearized model, maximum likelihood, and super-resolution localization[END_REF]. It can be seen from Figure 4 that the LS without regularization returns a very noisy IRF estimate that the leak cannot be localized. By contrast, from the results of the other three methods, the jump in the IRF due to leak can be clearly observed and GCV returns a less noisy IRF estimate. In order to detect a leak, the error of IRF estimate must be much lower than the jump in IRF due to wave reflection at the leak. For example, as can be seen in Figure 4 where the leak size is m 2 the height of jump is approximately 0.8 (we numerically test this case but the result is not shown in the present paper), such that the leak detection tolerates a higher IRF estimation error and a higher noise level. However, in order to detect a small leak which is the main concern of the leak detection problem, an IRF extraction method with error as low as possible is desirable.

Finally, let us revisit the results shown in Figure 1 that uses the experimental data in Ref. [START_REF] Lee | Valve design for extracting response functions from hydraulic systems using pseudorandom binary signals[END_REF]. Here, the cross-correlation, GCV and L-curve methods are used to extract the IRF; the corresponding results are shown in Figure 5. In this case, the cross-correlation method returns a very noisy result. The Tikhonov regularization parameters  decided by GCV and L-curve are respectively 20 and 210. The fluctuations in the results of GCV and L-curve are relatively small comparing with the main reflection at the boundary. In particular, L-curve returns almost the optimal results in the sense of least fluctuation in IRF estimation, which can be observed by comparing Figure 5(c) with Figure 1. Note that different from the random noise added in the previous simulation example, model uncertainties appear in real experiments, for example the corner sharpness of the PRBS signal. Tikhonov regularization is able to deal with uncertainties while the crosscorrelation method seems to only cancel random noises.

CONCLUSION

This paper proposes the generalized cross validation (GCV) and L-curve methods to decide the optimal Tikhonov regularization parameter in the least square (LS) deconvolution method for impulse response function (IRF) extraction of water pipeline. The methods proposed in this paper are also compared with the cross-correlation method. The performance of these methods is evaluated by the estimation of IRF with pseudo random binary sequence (PRBS) input signal. The numerical and experimental results show that the choice of IRF extraction method is important. The availability of the proposed methods is also justified by experimental data.

ACKKNOWLEDGEMENT

Figure 1 .

 1 Figure 1. Impulse response function obtained from the Eq. (5) with various Tikhonov regularization parameter  .

  (a) shows the input signal of the first two seconds, while Figure 2 (b) displays the output signal for 0-20 s.

Figure 2 .

 2 Figure 2. Normalized (a) input signal (PRBS) and (b) output signal at the downstream end of the pipe.

Figure 3 .

 3 Figure 3. Singular values  n in descending order of X , normalized by their maximum value  1 .

  the jump corresponding to the reflection at the leak is 0.2. For a larger leak, the jump in IRF increases, for example if

  This work has been supported by research grants from the Research Grant Council of the Hong Kong SAR, China (Project No. T21-602/15R).

Figure 4 .

 4 Figure 4. Estimation of IRF normalized by its maximum value. SNR=10 dB.

Figure 5 .

 5 Figure 5. Impulse function obtained from the Eq. (5) using the crosscorrelation, GCV and L-curve methods.