Filipe Fontanela
email: filipe.fontanela15@imperial.ac.uk

Aurelien Grolet

Loic Salles

Amin Chabchoub

Alan R Champneys

Sophoclis Patsias

Norbert Hoffmann

Dissipative solitons in forced cyclic and symmetric structures

Keywords: Travelling wave excitation, solitons, vibration localisation, cyclic structures

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

Introduction

Localisation of vibrations in cyclic and symmetric structures is an important topic for rotating machines due to high cycle fatigue. In the case of linear vibrations, localisation may occur due to a lack of symmetry, induced e.g. by manufacturing variabilities or wear [START_REF] Hodges | Confinement of vibration by structural irregularity[END_REF][START_REF] Hodges | Vibration isolation from irregularity in a nearly periodic structure: Theory and measurements[END_REF][START_REF] Bendiksen | Mode localization phenomena in large space structures[END_REF]. This problem, usually referred to as a mistuning problem, has attracted much attention in the literature, encompassing areas such as efficient numerical tools for analysis [START_REF] Castanier | Modeling and analysis of mistuned bladed disk vibration: current status and emerging directions[END_REF], experimental investigations [START_REF] Judge | Experimental Investigation of Mode Localization and Forced Response Amplitude Magnification for a Mistuned Bladed Disk[END_REF], and even the use of intentional mistuning to minimise localised vibrations [START_REF] Castanier | Using Intentional Mistuning in the Design of Turbomachinery Rotors[END_REF].

In real engineering applications, the assumption of linear vibrations has to be seen as an approximation, and the dynamics of complex engineering structures is usually a result of multiple linear and nonlinear phenomena. In the case of bladed disks of aircraft engines, nonlinearities may arise e.g. due to friction induced by internal joints [START_REF] Krack | Vibration prediction of bladed disks coupled by friction joints[END_REF][START_REF] Pesaresi | Modelling the nonlinear behaviour of an underplatform damper test rig for turbine applications[END_REF], or due to large deformations induced by strong excitations [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF]. In the nonlinear regime, the processes which lead to localised vibrations may go beyond mistuning. It is well-known, for example, that even perfect cyclic and symmetric structures may localise vibrations due to bifurcations [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Vakakis | Mode Localization in a Class of Multidegreeof-Freedom Nonlinear Systems With Cyclic Symmetry*[END_REF][START_REF] Vakakis | A Multiple-Scales Analysis of Nonlinear, Localized Modes in a Cyclic Periodic System[END_REF][START_REF] King | A very complicated structure of resonances in a nonlinear system with cyclic symmetry: Nonlinear forced localization[END_REF][START_REF] Vakakis | Normal modes and localization in nonlinear systems[END_REF][START_REF] Georgiades | Modal Analysis of a Nonlinear Periodic Structure with Cyclic Symmetry[END_REF]. Moreover, the effects of nonlinearities in nonhomogeneous systems is still a very active topic of research not only in structural dynamics [START_REF] Capiez-Lernout | Mistuning analysis and uncertainty quantification of an industrial bladed disk with geometrical nonlinearity[END_REF], but also in other research fields [START_REF] Ikeda | Intrinsic Localized Modes of Harmonic Oscillations in Nonlinear Oscillator Arrays[END_REF][START_REF] Ivanchenko | Anderson localization or nonlinear waves: A matter of probability[END_REF][19].

This paper focuses on the nonlinear dynamics of cyclic and symmetric structures that are excited by forces that have the form of travelling waves. In the case of bladed disks, such excitations, also known as engine orders, are caused by aerodynamic or structural forces that move along the structure preserving their shape. Analysis is based on a minimal model composed of symmetric Duffing oscillators, cyclically connected to each other, and in the presence of linear velocity-proportional dampers. It has already been shown that, in the conservative regime, the underlying model has localised states composed of envelope solitons [START_REF] Grolet | Travelling and standing envelope solitons in discrete non-linear cyclic structures[END_REF][START_REF] Fontanela | Dark solitons, modulation instability and breathers in a chain of weakly nonlinear oscillators with cyclic symmetry[END_REF]. These results were obtained from analytical solutions of a Nonlinear Schrödinger Equation (NLSE). For the non-conservative case, here a modified NLSE is developed, and analytical solutions describing localised states are not possible anymore. Therefore, general solutions and their stability are assessed numerically. Within this regime, it is demonstrated that localised states, sometimes called dissipative solitons [START_REF] Grelu | Dissipative solitons for mode-locked lasers[END_REF][START_REF] Akhmediev | Dissipative solitons: From optics to biology and medicine[END_REF], branch off from the travelling waves, resulting in complex bifurcation diagrams composed of stable and unstable single and multi-soliton states. Finally, these localised states are continued into the conservative limit and it is demonstrated that they become nonlinear normal modes of the underlying system. These modes are forced to bifurcate from the backbone curve of the travelling wave at finite vibration amplitude.

The paper is organized as follows. In Sec. 2 the minimal model for weakly nonlinear cyclic and symmetric structures is introduced. Moreover, the assumptions necessary to derive a non-conservative NLSE are discussed, and strategies to compute localised solutions and check their stability are presented. General numerical results are shown in Sec. 3. We discuss how localised states branch off from travelling waves in the NLSE framework, and these results are veri-fied through time integration of the full equations. In Sec. 4, localised states obtained from the non-conservative NLSE are continued into the conservative limit. Bifurcation diagrams show that localised states detach from travelling waves at non-zero amplitudes, contrary to the standard analytical predictions obtained from the conservative NLSE solved in an infinite spatial domain. Finally, Sec. 5 discusses the main conclusions and suggests directions for further investigations.

The minimal model and solution methods

The physical system under investigation is displayed in Fig. of N s blocks of mass m, cyclically connected through linear springs k c , and attached to the ground by linear and cubic springs k l and k nl , respectively. The dissipation is modelled by linear viscous dampers c, connected to the ground, while the external force applied to the nth mass is f n . The system in Fig. 1 can be understood as a minimal model of a bladed disk [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF][START_REF] Grolet | Vibration Analysis of a Nonlinear System With Cyclic Symmetry[END_REF] operating in the strong deformation regime, while a positive value of k nl models the induced hardening effect. The displacement of the nth degree of freedom is given by

mü n + c un + k l u n -k c (u n-1 + u n+1 -2u n ) + k nl u 3 n = f n . (1)
This equation can be written in rescaled variables, for convenience, as

ün + γ 2 un + ω 2 0 u n -ω 2 c (u n-1 + u n+1 -2u n ) + ξu 3 n = g n , (2)
where

γ 2 = c/m, ω 2 0 = k l /m, ω 2 c = k c /m, ξ = k nl /m and g n = f n /m.
In this paper, the investigation focuses on responses to travelling wave excitation. The system dynamics induced by these forces is a typical analysis in bladed disks vibrations. Travelling wave excitation is induced in rotating machines e.g. due to aerodynamic forces that move along the structure [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] or due to engine orders [START_REF] Jones | Traveling wave excitation system for bladed disks[END_REF]. Mathematically, a travelling wave excitation can be written as

g n (t) = F 0 exp{i [k(n -1)a -ω f t]} + c.c., (3)
where F 0 is the force amplitude, k is its wave number, a = 2π/N s is the lattice parameter, ω f is the external frequency, i is the imaginary unity, and c.c. represents the complex-conjugate of the preceding expression. It should be noted that the wave number k can only assume integer numbers and, consequently, the external force respects the cyclic symmetry (g Ns+1 = g 1 ).

Linear analysis of the conservative system

In the linear regime (ξ = 0) and in the absence of damping (γ 2 = 0) and external forces (F 0 = 0), the system depicted in Fig. 1 has exact travelling wave solutions of the form

u n (t) = U 0 exp{i[k(n -1)a -ω k t]} + c.c.. (4)
In Eq.( 4) the constant U 0 stands for the travelling wave amplitude, while

ω k = ω 2 0 + 4ω 2 c sin 2 ka 2
is its corresponding natural frequency. Again, the travelling wave solution is only valid when k is an integer and, therefore, the response fulfils the cyclic symmetry (u Ns+1 = u 1 ).

The slowly-varying and weakly nonlinear regime

In the presence of damping (γ 2 = 0) and external forces (F 0 = 0), this paper focuses on the slowly-varying weakly nonlinear regime of Eq. (2). Thus, the system response u n is written as

u n (t) = Ψ(X, T ) exp{i[k(n -1)a -ω k t]} + c.c., (5)
where Ψ(X, T ) is an envelope function which modulates the travelling wave, and is a small parameter. In Eq. ( 5), X and T are continuous variables such as X = x = (n -1)a and T = t. In the undamped case and in the absence of external forces, it is well-known that the evolution of the envelope function Ψ can be approximated by a conservative NLSE [START_REF] Grolet | Travelling and standing envelope solitons in discrete non-linear cyclic structures[END_REF][START_REF] Remoissenet | Waves Called Solitons: Concepts and Experiments[END_REF]. Mathematically, such an equation is obtained after substituting Eq. ( 5) into Eq. ( 2), rejecting higherorder harmonics, and applying a continuum approximation. This assumption works well when the envelope function Ψ varies much more slowly in space and time when compared to the carrier wave period (2π/ω k ) and the lattice parameter (2π/N s ), respectively. For the non-conservative case, a similar approach can be applied if the following additional assumptions are considered: (1) the damping value γ 2 and the newly defined detuning parameter δ ω = ω k -ω f are at order 2 ; and (2) the external force level is at order 3 . As result, a similar equation for the evolution of Ψ is obtained such as

i ∂Ψ ∂τ + P ∂ 2 Ψ ∂η 2 + Q|Ψ| 2 Ψ = -iΓΨ -h exp{iδ ω τ }, (6)
where τ = T , while η = X -c g t is a frame moving with the group velocity c g = dω k dk . In Eq [START_REF] Castanier | Using Intentional Mistuning in the Design of Turbomachinery Rotors[END_REF], the parameter P = 1 2 d 2 ω k dk 2 accounts for the dispersion, Q = -3ξ 2ω k is the nonlinear coefficient, Γ = γ 2 2 models the linear damping effect, and h = F0 2ω k is the external force. The derivation of Eq. ( 6) is presented in the Appendix. For the undamped and unforced case, the right-hand side of Eq. ( 6) is null and the system recovers the conservative NLSE. Within this regime, the NLSE is an integrable system and localised solutions, consisting of envelope solitons, can be obtained analytically (see e.g. Ref. [START_REF] Remoissenet | Waves Called Solitons: Concepts and Experiments[END_REF]). For the forced and damped case Eq. ( 6) loses integrability and localised solutions written analytically are not available in the literature. Therefore, solitons of the damped and driven NLSE are studied numerically in the following.

Solution approach

As usual it is assumed that the travelling response waves are vibrating with the same frequency as the external excitation such as

Ψ(τ, η) = ψ(τ, η) exp{iδ ω τ }. (7)
After substituting Eq. ( 7) into Eq. ( 6), the process leads to the autonomous system

i ∂ψ ∂τ -δ ω ψ + P ∂ 2 ψ ∂η 2 + Q|ψ| 2 ψ = -iΓψ -h. (8)
For the homogeneous ( ∂ 2 ψ ∂η 2 = 0) and standing ( ∂ψ ∂τ = 0) case, solutions of Eq. ( 8) can still be obtained analytically. In the case of localised standing waves, it has been shown, numerically, that they detach from the homogeneous solutions through bifurcations [START_REF] Terrones | Stability and Bifurcation of Spatially Coherent Solutions of the Damped-Driven NLS Equation[END_REF][START_REF] Barashenkov | Existence and stability chart for the ac-driven, damped nonlinear Schrödinger solitons[END_REF][START_REF] Barashenkov | Bifurcation to multisoliton complexes in the ac-driven, damped nonlinear Schrodinger equation[END_REF]. In these previous publications dissipative solitons were computed using the external excitation as a continuation parameter. However, in vibration engineering the standard way to characterise mechanical structures relies on the behaviour of frequency response functions; i.e. keeping h constant while δ ω is varied. Therefore, we solve Eq. ( 8) using δ ω as the continuation parameter.

Another difference from previous studies (see e.g. Refs. [START_REF] Barashenkov | Existence and stability chart for the ac-driven, damped nonlinear Schrödinger solitons[END_REF][START_REF] Barashenkov | Bifurcation to multisoliton complexes in the ac-driven, damped nonlinear Schrodinger equation[END_REF]) is the way we have computed localised solutions as bifurcations from standing waves ( ∂ψ ∂τ = 0) rather than as continuation from the undamped limit. We use the popular continuation code AUTO by first rewriting Eq. (8) as a system of two first-order complex nonlinear differential equations and then solving a two-point boundary-value problem. The theoretical framework implemented in AUTO can be checked in Refs. [START_REF] Doedel | Auto97, Continuation and bifurcation software for ordinary differential equations[END_REF][START_REF] Champneys | A numerical toolbox for homoclinic bifurcation analysis[END_REF] and references therein.

Stability analysis

In this paper, the approach developed in Ref. [START_REF] Barashenkov | Existence and stability chart for the ac-driven, damped nonlinear Schrödinger solitons[END_REF] is used to assess the stability of stationary solutions ψ 0 (η) obtained with AUTO. Thus, the expression

ψ(η, τ ) = ψ 0 (η) + ∆ψ(η, τ ) (9)
is substituted into Eq. ( 8), where ∆ψ is a small perturbation. When the analysis is considered up to order ∆ψ, which results in a linear investigation, the stability analysis yields

J ∂y ∂τ + Γy = Hy. (10)
In Eq. ( 10), the perturbation ∆ψ is split into real and imaginary parts such as y(η, τ ) = (Re{∆ψ} Im{∆ψ}) T , while J and H are the two matrices depending on the real ψ R = Re{ψ 0 } and imaginary ψ I = Im{ψ 0 } parts of ψ 0 written as

J = 0 1 -1 0 , (11)
H = -P ∂ 2 ∂η 2 + δ ω -Q(3ψ 2 R + ψ 2 I ) -2Qψ I ψ R -2Qψ I ψ R -P ∂ 2 ∂η 2 + δ ω -Q(3ψ 2 I + ψ 2 R ) . (12)
After separating the time and space variables via y(η, τ ) = z(η) exp{λτ }, the linear stability problem becomes

µJz(η) = Hz(η), (13)
where µ = λ + Γ. Therefore, if Re{µ} < Γ the solution is stable, otherwise it is unstable. For the special case of homogeneous solutions, stability can be addressed analytically. It has already been proven that flat solutions of Eq. ( 8) can be modulationally unstable [START_REF] Terrones | Stability and Bifurcation of Spatially Coherent Solutions of the Damped-Driven NLS Equation[END_REF][START_REF] Barashenkov | Existence and stability chart for the ac-driven, damped nonlinear Schrödinger solitons[END_REF]. For the case of localised solutions, stability has to be investigated numerically. Therefore, the solution ψ 0 is discretized in space and the second-order derivative in Eq. ( 12) is calculated through finite differences. A detailed explanation about the implementation of this numerical procedure can be checked in Ref. [START_REF] Barashenkov | Existence and stability chart for the ac-driven, damped nonlinear Schrödinger solitons[END_REF].

Numerical results

The cyclic and symmetric system displayed in Fig. 1 is investigated by assuming

N s =48 oscillators, ω 2 0 =1 s -2 , ω 2 c = 1 s -2 , ξ = 0.1 m -2 s -2
, and γ 2 =0.01 s -1 . Figure 2 displays the natural frequency ω k of the linear system, the group velocity c g , the dispersion parameter P , and the nonlinear coefficient Q, as functions of the wave number k. One should note that the system is cyclic and, therefore, the curves in Fig. 2 are symmetric with respect to k=24. Moreover, a transition from positive to negative values of P , at k ≈ 10 m -1 , is highlighted. This boundary separates two different regimes: for P Q < 0 the system is in the defocusing range, while for P Q > 0 the system is in the focusing range. It is well-known (see Refs. [START_REF] Remoissenet | Waves Called Solitons: Concepts and Experiments[END_REF][START_REF] Dauxois | Physics of Solitons[END_REF]) that plane waves in the focusing (defocusing) range and in an infinite domain are unstable (stable) when they are investigated in a conservative NLSE framework. Within the focusing range, solitons are stable solutions of the conservative system. In this section we investigate if soliton solutions also exist in the non-conservative regime. The strategy of solution proposed in Sub. 2.3 is applied to Eq. ( 8). Thus, standing solutions of the NLSE are investigated with AUTO. Due to the cyclic symmetry, any arbitrary solution ψ 0 (η) when translated in space is also a solution of Eq. [START_REF] Pesaresi | Modelling the nonlinear behaviour of an underplatform damper test rig for turbine applications[END_REF]. In order to remove this symmetry in the continuation approach, Eq. ( 8) is solved for η = [0, π] assuming Neumann boundary conditions ( ∂ψ ∂η | η=0 = ∂ψ ∂η | η=π = 0). In this case, any general solution is forced to localise maximums and minimums of energy at the spatial boundaries η = and η = π. Moreover, the solution for the remaining spatial configuration (η = (π, 2π)) is obtained by simply reflecting the solution calculated for η = [0, π]. However, one should note that this implementation restricts the analysis to solutions which are symmetric with respect to η = π.

In the case of a single-hump soliton, the analysis should focus on the maximum amplitude of ψ. This value, when compared to the corresponding homogeneous states, identifies the level of localisation for an specific solution. However, the same measure is not sensitive to the number of humps in a multiple-hump soliton. Therefore, a new measure

E = 2π 0 |ψ| 2 + ∂ψ ∂η 2 + |ψ| 4 dη ( 14
)
is introduced. One should note that E is slightly different from the energy of the envelope function used in Ref. [START_REF] Barashenkov | Bifurcation to multisoliton complexes in the ac-driven, damped nonlinear Schrodinger equation[END_REF].

Figure 3 displays E as function of the external forcing frequency, calculated from the NLSE implementation. The result is obtained by assuming the wave number k=16 m -1 and the external excitation level F 0 =0.007 ms -2 . This configuration leads the parameters of Eq. ( 8) to P =-0.00295 m 2 rad s -1 , Q=-0.00750 m -2 rad -1 s -1 , Γ=0.005 s -1 and h=0.00175 ms -1 . It is possible to verify, from Fig. 3, that the plane waves loose stability at ω f ≈ 2.001 rad/s and localised solutions branch off from the travelling waves. These localised states consist of two hump solitons (see P2 and P3) which, as discussed before, are symmetric with respect to η = π m. The branch detaches at ω f ≈ 2.001 rad/s with stable soliton solutions, but it looses stability after ω f ≈ 2.008 rad/s. The same analysis also shows two other branches of localised solutions: one detaching at ω f ≈ 2.008 rad/s; and another detaching at ω f ≈ 2.010 rad/s.

An important feature, observed in the bifurcation diagram of Fig. 3, is that localised solutions with one hump detaching from travelling waves do not necessarily continue to homogeneous states keeping only one hump as solutions. In fact, they can internally bifurcate from one hump to two humps (see points P1 and P2 in Fig. 3). This behaviour resembles snaking bifurcations [START_REF] Hunt | Cellular buckling in long structures[END_REF][START_REF] Papangelo | Snaking bifurcations in a self-excited oscillator chain with cyclic symmetry[END_REF][START_REF] Papangelo | Multiple spatially localized dynamical states in friction-excited oscillator chains[END_REF], in which a branch of localised solutions fold to generate states with more or less localised elements, or a more or less widely extending zone of localised dynamics. In fact, a similar result has been already observed for solutions of Eq. ( 8), but keeping the detuning parameter δ ω constant and varying the external force level h [START_REF] Barashenkov | Bifurcation to multisoliton complexes in the ac-driven, damped nonlinear Schrodinger equation[END_REF].

In order to verify the predictions, initial states u n (0) and un (0) were computed from the NLSE implementation and the equations of the original physical system (Eq. ( 2)) are calculated with a time-marching approach. Figure 4 illustrates the responses obtained assuming the stable states P1 and P2 of Fig. 3. The initial displacement, in Panel (a), shows a localised hump on a travelling wave background. This hump keeps its shape while it moves around the system with roughly the group velocity, as shown by the displacement u 1 (t) in Panel (b). Panels (c) and (d) show similar results, but starting with a symmetric localised state composed of two humps. In this case, the displacement field u 1 (t) shows a doubly localised state moving around the system with a constant envelope profile.

In the case of higher external forcing amplitude the bifurcation diagram is more complicated. Figure 5 displays results obtained in the same way as for Fig. 3, but calculated assuming a stronger forcing of F 0 =0.010 ms -2 (h=0.0025 ms -1 ). The result in Panel (a) shows that plane waves loose stability at lower values of ω f . Moreover, additional branches of localised solutions detach from the travelling waves. They are mainly composed of three humps solutions, as shown by P3 in Panel (b). However, most solutions are unstable, and the stable branches consist of one and two localised humps states.

Solutions obtained from time integration of states P1 and P2, identified in Fig. 5, are displayed in Fig. 6. The results are similar to the ones presented in Fig. 4, where solutions of one and two localised humps keep their shapes while they move around the system with approximately the group velocity.

In order to check the applicability of the NLSE approach to a different wave 8) for k=16 m -1 and F 0 =0.007 ms -2 . Panel (a) illustrates E as a function of ω f . The black dashed lines denote stable plane wave results, while the red one shows unstable plane wave solutions. The solid lines, which branch off from the travelling wave result, indicate stable (black) and unstable (red) localised solutions. The vertical dotted line shows the linear natural frequency of the system. Panel (b) shows the real (black) and imaginary (blue) parts of Ψ 0 at P1, P2, P3, and P4. the same strategy is applied assuming k=22 m -1 and F 0 =0.0011 ms -2 . This configuration leads the parameters of Eq. ( 8) to P =-0.00378 m 2 rad s -1 , Q=-0.00675 m -2 rad -1 s -1 , Γ=0.005 s -1 and h=0.0025 ms-1. The corresponding results are displayed in Fig. 7.

The results in Panel (a) and (b) are very similar to the ones already discussed in Figs. 3 and5. However, no stable solution is observed between ω f ≈ 2.226 rad/s and ω f ≈ 2.231 rad/s. This result suggests that stable configurations within this frequency range might not be symmetric with respect to η = π and, therefore, we can not compute them with the present method.

Solutions from points P1 and P2 are, again, used as initial conditions for time-integration of the physical system. The corresponding solutions are displayed in Fig. 8. It is interesting to notice that, although solutions are symmetric with respect to η = π m, the humps are not necessarily equidistant in space. This feature can be verified from x 1 (t) in Panel (d) of Fig. 8, where the two humps are moving with the same velocity but they are clearly not equally spread over the physical system.

In the case of higher external forces amplitude, the results are very similar to the ones presented in Fig. 5. Figure 9 displays the corresponding results calculated for k=22 m -1 and F 0 =0.016 ms -2 (h=0.0035 ms -1 ).

Again a quite complicated bifurcation diagram is observed. In this case, solutions with one, two, three and four humps detach from travelling waves. However, as already observed in the previous cases, most configurations are unstable. The stable solutions consist, again, of one (P1) and two (P2) hump solitons.

Time-marching results in Fig. 10, obtained from the equations of the physical system, show similar results as in Figs. 4,6, and 8. Localised states which keep their profile move around the system with roughly the group velocity.

Continuation to the conservative limit

The analysis in Sec. 3 demonstrated that localised states branch off from unstable travelling waves. In order to understand how these dissipative solitons are related to the solitons in the conservative limit, a continuation from the driven-dissipative to the conservative case is conducted. To link the regimes, the autonomous NLSE in Eq. ( 8) is rewritten as

i ∂ψ ∂τ -δ ω ψ + P ∂ 2 ψ ∂η 2 + Q|ψ| 2 ψ = -ζ(iΓψ + h), (15)
where ζ is a continuation parameter. One should note that continuing ζ from ζ=1 to ζ=0 allows to follow driven-dissipative solutions of Sec. 3 into the conservative regime. Figure 11 displays the results obtained for the wave number k=16 m -1 . In this case, contrary to analytical solutions of the NLSE in an infinite spatial domain, travelling waves lose stability at non-zero amplitudes. Moreover, at this point, a branch of stable localised solution composed of solitons with one hump detaches from the travelling waves. Initially, the soliton consists 9: Solution obtained from the NLSE implementation in Eq. ( 8) for k=22 m -1 and F 0 =0.016 ms -2 . Panels (a) and (b) illustrate the same quantities as in Fig. 3. a small hump on a travelling wave background. However, when the amplitude of the solution gets larger and nonlinear effects are stronger, this hump tends to localise more. Another important feature observed is that, at higher background amplitudes, other branches of localised states with two, three, four, and five humps detach, successively, from the homogeneous solutions. However, all the computed multi-solitons states are unstable within the dynamic analysed regime.

In order to check the NLSE results, initial conditions resulting from points P1 and P2 in Fig. 11 were computed to solve the equations of the physical system by time integration. Figure 12 displays the corresponding results. It is possible to observe, in Panel (a) of Fig. 12, that the initial conditions consist of one hump modulating the travelling wave. Time-marching results, in Panel (b), show that this hump moves around the system keeping its shape, as expected from the NLSE theory. In Panel (c), where displacements are larger and the system is in a stronger non-linear regime, only a small part of the system vibrates, while the other part stands at rest. Results from time-integration, in Panel (d), show that in this case a zone of localised vibration moves around the system almost preserving its shape. However, it is possible to observe small vibrations in between the successive humps moving around the structure. These small vibrations are not expected from the NLSE results (see P2 in Fig. 11). This discrepancy is mainly due to the two assumptions imposed by the NLSE approximation: (1) the continuous approach, which looses accuracy when solutions are very localised; and (2) due to large displacements, when the single harmonic approximation is not sufficient anymore.

Figure 13 shows analogous results for another wave number k=22 m -1 . The solutions show, again, a one hump stable soliton detaching from unstable travelling waves with a non-zero amplitude. Other localised states with two, three, four, and five humps also detach from the homogeneous solution. However, as in the previous case, these multi-solitons solutions are unstable within the regime analysed.

Figure 14 shows the corresponding results for P1 and P2 identified in Fig. 13. The results show, again, that a localised hump moves around the system preserving its shape.

Conclusions and outlook

This paper focuses on localised oscillations arising from modulation of nonlinear travelling waves in weakly nonlinear cyclic and symmetric structures. The study is based on a minimal model composed of Duffing oscillators in the presence of linear viscous dampers and external travelling wave forces. The system can be understood as a minimal model of a bladed disk operating in the large deformations regime. It is demonstrated that, in the non-conservative regime, localised states composed of single and multi-soliton solutions bifurcate from travelling waves. Moreover, when these localised states are continued to the conservative regime, it is shown that solitons branch off from plane waves at finite amplitudes. future work, more complex physical models need to be studied. Next investigations will focus on more complicated models, e.g. with two degrees of freedom per sector which also include the disk inertia. More realistic models, e.g. based on reduced-order systems obtained from full finite element models, should also be in the focus of future investigations. The phenomenological findings discussed in this paper also need to be compared to small-scale experiments. In addition, localised vibrations arising from other nonlinear models, such as friction and impact, are also an avenue of future research. the following, in order to introduce a continuum approximation, the envelope function at a neighbour degree of freedom Ψ(na ± a) is computed using a Taylor expansion such as

Ψ(na ± a) = Ψ ± a ∂Ψ ∂X + a 2 2 2 ∂ 2 Ψ ∂X 2 + o(> 2 ), (A.5)
where o(> 2 ) stands for terms higher than 2 . Moreover, the following additional assumptions are considered: (1) the damping γ 2 and the detuning parameter δ ω are at order 2 ; and (2) the external force is at order 3 . After substituting the previous expressions into Eq. ( 2), ignoring higher-order harmonics, and considering the dispersion relation ω 2 k = ω 2 0 + 4ω 2 c sin 2 ka 2 , the process leads to

i 2 ∂Ψ ∂T + c g ∂Ψ ∂X + 3 2ω k Λ ∂ 2 Ψ ∂X 2 - ∂ 2 Ψ ∂T 2 -3ξ|Ψ| 2 Ψ + iω k γ 2 Ψ + F 0 exp{i δ ω T } = 0, (A.6)
where c g = 2aω 2

c ω k sin(ka) is the group velocity and Λ = ω 2 c a 2 cos(ka). In order to simplify Eq. (A.6), it is rewritten in a frame η moving with c g such as η = X -c g T and τ = T . Therefore, the following relations are obtained

∂ 2 Ψ ∂η 2 - 3ξ 2ω k |Ψ| 2 Ψ = -i γ 2 2 Ψ - F 0 2ω k exp{iδ ω τ }.
(A.11) Equation (A.11) is an externally forced and linearly damped NLSE, which can be written in the form of Eq. ( 6) by setting P = 1 2

d 2 ω k dk 2 = 1 2ω k (Λ -c 2 g ), Q= -3ξ 2ω k , Γ = γ 2
2 , and h = F0 2ω k .

Figure 1 :

 1 Figure 1: Physical system under investigation. In the left-hand side, an illustration of the full system, while the right-hand side displays any three neighbouring degrees of freedom.

Figure 2 :

 2 Figure 2: Evolution of the system parameters as a function of the wave number k: (a) natural frequency ω k ; (b) group velocity cg; (c) dispersion parameter P ; and (d) nonlinear coefficient Q. The vertical dashed line in Panel (c) illustrates the transition from positive to negative values of P .

Figure 3 :

 3 Figure 3: Solutions obtained from the NLSE implementation in Eq. (8) for k=16 m -1 and F 0 =0.007 ms -2 . Panel (a) illustrates E as a function of ω f . The black dashed lines denote stable plane wave results, while the red one shows unstable plane wave solutions. The solid lines, which branch off from the travelling wave result, indicate stable (black) and unstable (red) localised solutions. The vertical dotted line shows the linear natural frequency of the system. Panel (b) shows the real (black) and imaginary (blue) parts of Ψ 0 at P1, P2, P3, and P4.

Figure 4 :

 4 Figure 4: Results obtained through numerical time integration of the physical system for k=16 m -1 and F 0 =0.007 ms -2 . Panel (a): the spatial displacements at t=0 s calculated at P1 in Fig. 3; Panel (b): time-evolution of the displacement u 1 (t) of the first oscillator. Panels (c) and (d) illustrate the same quantities for P2.

Figure 5 :

 5 Figure 5: Solutions obtained from the NLSE implementation in Eq. (8) for k=16 m -1 and F 0 =0.010 ms -2 . Panels (a) and (b) illustrate the same quantities as in Fig.3.

Figure 6 :

 6 Figure 6: Results obtained through numerical integration of the physical system for k=16 m -1 and F 0 =0.010 ms -2 . Panels (a), (b), (c), and (d) illustrate the same quantities as in Fig. 4.

Figure 7 :

 7 Figure 7: Solution obtained from the NLSE implementation in Eq. (8) for k=22 m -1 and F 0 =0.0011 ms -2 . Panels (a) and (b) illustrate the same quantities as in Fig.3.

Figure 8 :

 8 Figure 8: Results obtained through numerical integration of the physical system for k=22 m -1 and F 0 =0.0011 ms -2 . Panels (a), (b), (c), and (d) illustrate the same quantities as in Fig. 4.

Figure 10 :

 10 Figure 10: Results obtained through numerical integration of the physical system for k=22 m -1 and F 0 =0.016 ms -2 . Panels (a), (b), (c), and (d) illustrate the same quantities as in Fig. 4.

Figure 11 :

 11 Figure 11: Localised states calculated for the conservative physical system for k=16 m -1 . Panel (a) displays the bifurcation diagram, where the black dashed line represents stable homogeneous solutions, while the red one shows the unstable regime. Full lines, detaching from travelling waves, represent stable (black) and unstable (red) localised solutions. Panel (b) shows the real (black) and imaginary (blue) parts of the ψ 0 solutions at P1, P2, P3, and P4 of Panel (a).

Figure 12 :

 12 Figure 12: Results obtained through numerical integration of the conservative physical system k=16 m -1 . Panels (a), (b), (c), and (d) illustrate the same quantities as in Fig. 4.

Figure 13 :

 13 Figure 13: Results obtained for the conservative system assuming k=22 m -1 . Panels (a) and (b) illustrate the same quantities as in Fig. 11.

Figure 14 :

 14 Figure 14: Results obtained through numerical integration of the conservative physical system for k=22 m -1 . Panels (a), (b), (c), and (d) illustrate the same quantities as in Fig. 4.

  1. It consists

	𝑢 𝑛-1		𝑢 𝑛		𝑢 𝑛+1
	𝑘 𝑐	𝑘 𝑐	𝑘 𝑐		𝑘 𝑐
	𝑚	𝑘 𝑙	𝑚	𝑘 𝑙	𝑚	𝑘 𝑙
		𝑘 𝑛𝑙		𝑘 𝑛𝑙		𝑘 𝑛𝑙
		𝑐		𝑐		𝑐
	𝑓 𝑛-1		𝑓 𝑛		𝑓 𝑛+1

Acknowledgements

The authors thank Dr. Igor Barashenkov for discussions and advice regarding the non-conservative NLSE. F. Fontanela is funded by the Brazilian National Council for the Development of Science and Technology (CNPq) under the grant 01339/2015-3. The authors thank Rolls-Royce plc for allowing this publication.

Appendix A. Derivation of the non-conservative NLSE

The weakly nonlinear regime of the system described by Eq. ( 2) is considered assuming a travelling wave excitation such as

where the external force frequency ω f is decomposed using the linear natural frequency ω k and a detuning parameter δ ω such as ω F = ω k -δ ω . In order to derive a NLSE in the presence of linear viscous damping and external travelling wave excitation, it is assumed that the system is in the weakly nonlinear regime and experiencing travelling waves modulated by an envelope function such as

where X = x and T = t are two slowly-varying functions. Therefore, the corresponding velocity un and acceleration ün of the nth mass are written as .4)


Cédric Pradalier
email: cedric.pradalier@georgiatech-metz.fr

Pascal Pomarede

Jan Steckel
email: jan.steckel@uantwerpen.be

On-plate localization and mapping for an inspection robot using ultrasonic guided waves: a proof of concept

This paper presents a proof-of-concept for a localization and mapping system for magnetic crawlers performing inspection tasks on structures made of large metal plates. By relying on ultrasonic guided waves reflected from the plate edges, we demonstrate that it is possible to recover the plate geometry and robot trajectory to a precision comparable to the signal wavelength. The approach is tested using real acoustic signals acquired on test metal plates using lawn-mower paths and random-walks. To the contrary of related works, this paper focuses on the practical details of the localization and mapping algorithm.

I. INTRODUCTION

This paper aims at demonstrating the interest of using ultrasonic guided waves to support inspection robots operating on structures made of metal plates. Such structures include in particular ship outer hulls and large storage tanks, as depicted in 1. On such a structure, localization with respect to individual plates can be beneficial to precisely triangulate defects such as corrosion patches [START_REF] Druet | Passive SHM system for corrosion detection by guided wave tomography[END_REF], [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF] or even attempt acoustic tomography as in [START_REF] Cailly | Three dimensional ultrasonic imaging of mechanical components by inversion[END_REF]. In combination with odometry and an external localization system (laser theodolite, Ultra-Wide Band (UWB) beacons...), this would also lead to precise absolute localization of the inspection results. On metal plates, Ultrasonic Guided Waves (UGWs) can mostly be generated using piezo-electric transducers in contact with the plate. When the relation between the frequency, the plate thickness and the wave velocity is right, these waves propagate radially around the emitter through the plate material, like ripples around a stone thrown into a pool. For metal plates relevant to inspection tasks, frequencies in the range of hundreds of kHz are typically appropriate, with wave length of the order of a few centimeters.

In this paper, we are considering short bursts of a few wave cycles emitted from a transducer mounted on a mobile platform. When encountering the plate edges, these waves are reflected normally and these reflections are measured by a transducer co-located with the emitter to be converted into electrical signals. From these signals, it is possible to identify the reflected waves and from them recover the distance to the edge. In a practical case, while working on a rectangular plate, the acquired signal may include reflections from multiple edges as well as spurious detections. These sequences of edge detections are essentially range-only measurements to a line. Our intent is to use these ranges to both reconstruct the geometry of the plate and the localization of the emitter-receiver device over time. In the robotic community, this is obviously a SLAM problem with a nontrivial data-association challenge to identify which edge has been observed at a given time. Although a number of works have considered the very similar problem of room shape reconstruction from acoustic echoes as a SLAM problem, most of them only consider simulated measurements and do not describe exhaustively their map management strategies.

In summary, the contributions of this paper are the following:

1) a demonstration of the applicability of a sparsity-based reconstruction of the arrival times of the ultrasonic reflections in the plate using an L 1 -regularized least squares approach; 2) a demonstration of the applicability of room reconstruction methodologies to on-plate mapping with a detailed description of the map management strategy, from initialization to landmark addition and outlier removal; 3) a proof-of-concept for the application of FastSLAM [START_REF] Montemerlo | Fastslam: A factored solution to the simultaneous localization and mapping problem[END_REF] to the context of on-plate localization for inspection tasks.

II. RELATED WORKS

a) Non-Destructive Evaluation: Ultrasound-based inspection technique is the most common tool for Non-Destructive Evaluation (NDE) on metal plates. In the most common set-up, the transducers are put in transmission or in echo mode. In the first one, two transducers are aligned along the same axis, one transducer is used as an emitter and the second receives the signal after propagation through the sample. In the second mode, one transducer acts both as emitter and receiver and the signal is recorded after reflecting on the sample. This is typically used with high frequencies (i.e. 5MHz) for thickness measurements. Alternatively, this paper consider the use of Lamb waves, a specific type of UGWs. Those waves propagate along the metal plate, potentially for large distances. In industry, they have been used to detect defects in pipelines, rails, or large structures with success [START_REF] Su | Guided lamb waves for identification of damage in composite structures?: A review[END_REF], [START_REF] Cawley | The use of lamb waves for the long range inspection of large structures[END_REF], but never deployed on autonomous mobile systems. Another technique is to use signals obtained during the inspection from different positions localize defects using acoustic tomography [START_REF] Druet | Passive SHM system for corrosion detection by guided wave tomography[END_REF], [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF]. This technique is usually used in the Structural Health Monitoring (SHM) field with sensors embedded in the structures [START_REF] Larose | Passive correlation imaging of a buried scatterer[END_REF], [START_REF] Chehami | Detection and localization of a defect in a reverberant plate using acoustic field correlation[END_REF]. However, for such a technique to work, it is critical to know the location of the emitters and receivers with a good precision, which is a challenge for magnetic crawlers operating on an inspected structure. In general, the use of UGWs for defect detection is an active area of research in the NDE community, in particular for acoustic tomography [START_REF] Cailly | Three dimensional ultrasonic imaging of mechanical components by inversion[END_REF] or imaging discontinuities [START_REF] Walaszek | New testing method for the rapid an flexible ultrasonic imaging of large metal structures of composites by combining guided wave and matrix phased array technologies[END_REF]. However, most of this research assumes a known localization of the emitter and receiver with respect to the inspected plate, and to the authors' knowledge, there are no work considering that the UGWs themselves could both provide the defect detection capabilities and help inferring the plate geometry and sensor localization.

b) Echo detection: Time delay estimation is often performed using a matched filter approach, which allows high-resolution time-delay estimation under the condition that the signal bandwidth is sufficiently high, and that the auto-correlation function of the used signals is sufficiently narrow [START_REF] Turin | An introduction to matched filters[END_REF]. Generalizations of the matched filter to overcome spectral imperfections are often applied in situations where the signal spectrum is not flat, with the most notable example of this approach the Generalized Cross Correlation -Phase Transform. The GCC-PHAT approach is the de-facto standard in time-delay estimation, given that the signal bandwidth is sufficiently high [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF]. More recently, sparsity-based methods using an L 1 -regularized least squares approach have surfaced as super-resolution time-delay estimation, which have been extensively applied to ultrasonic time-delay estimation [START_REF] Steckel | Sparse decomposition of in-air sonar images for object localization[END_REF], [START_REF] Fontaine | Determining biosonar images using sparse representations[END_REF]. This approach promises the detection of highly-overlapping reflections, even under low Signal-to-Noise Ratio (SNR) conditions. In our approach, we will use this sparsity-based approach to solve the timeof-arrival estimation problem.

c) Localization and Mapping: Outside of the field of NDE, the closest problem to on-plate localization and mapping is acoustic room reconstruction. In this sub-field of acoustics and signal processing, there are a number of works trying to reconstruct the shape of a room based on acoustic echoes. Earlier works (e.g. [START_REF] Salvati | Sound source and microphone localization from acoustic impulse responses[END_REF]) assumed that multiple static sources where located in a polygonal room and that a moving receiver was collecting data. [START_REF] Peng | Room shape reconstruction with a single mobile acoustic sensor[END_REF] addressed the problem of a moving emitter-receiver device to reconstruct a room shape using a geometric solution. However, it used simulated data and did not assume any measurement noise. In a groundbreaking work, [START_REF] Kreković | Echoslam: Simultaneous localization and mapping with acoustic echoes[END_REF] showed that SLAM techniques could be use to estimate the shape of a convex room. This work however was only tested in simulation, and as a result assumed that all echoes were matched to room edges and that no spurious echoes were present. A series of work expanded on [START_REF] Kreković | Echoslam: Simultaneous localization and mapping with acoustic echoes[END_REF]. [START_REF] Krekovic | Look, no beacons! optimal all-in-one echoslam[END_REF] expresses the problem as an optimization but is again tested only in simulation and does not address explicitly the ambiguous matching of measurements to room walls. In all the works above the map management strategy is not described exhaustively: in particular they do not explicit how the geometry of the room is initialized and how map edges are included when first detected or removed when proven to be outliers. On a more practical standpoint, [START_REF] Pradhan | Smartphone-based acoustic indoor space mapping[END_REF] and [START_REF] Song | Indoor pedestrian selfpositioning based on image acoustic source impulse using a sensor-rich smartphone[END_REF] build on the article above to build a practical room reconstruction system running on a smart-phone. The latter would be used to emit sound pulses that would reflect from the room walls. In [START_REF] Song | Indoor pedestrian selfpositioning based on image acoustic source impulse using a sensor-rich smartphone[END_REF], a rectangular environment is assumed, with known sizes, and the focus is on localization. In [START_REF] Pradhan | Smartphone-based acoustic indoor space mapping[END_REF] on the other hand, both the localization and the wall estimation are considered with only the first echo but very little information is given on the map management: when to add a new wall, how to handle outliers, etc.

In the context of SLAM, data association has always been a critical steps, as described in reference works such as [START_REF] Neira | Data association in stochastic mapping using the joint compatibility test[END_REF] or [START_REF] Bailey | Data association for mobile robot navigation: A graph theoretic approach[END_REF]. FastSLAM [START_REF] Montemerlo | Fastslam: A factored solution to the simultaneous localization and mapping problem[END_REF] is an alternative formulation of the SLAM problem where the state is estimated by sampling and where each sample contains an estimate of the trajectory and the map collected along this trajectory. FastSLAM also allows to handle uncertain data association by including multiple samples for the various association hypotheses. This is particularly important for the problem at hand since the validity of wall hypotheses can sometimes only be estimated several steps after their initialization. To our knowledge, FastSLAM has not been used previously in the context of the acoustic room reconstruction and localization problem.

d) Summary: This work will use UGWs to build on the room reconstruction techniques within a FastSLAM framework. While moving an emitter-receiver pair on a metal plate, this will allow to recover the trajectory of the device and the geometry of the plate. To support fast and robust echo detection, we make use of an L 1 -regularized least squares approach to signal-delay estimation using a sparse signal model. This approach allows very accurate time-delay estimation of multiple overlapping reflections to be performed, even in low SNR conditions. The next section will describe how these components combine together to demonstrate the feasibility of using UGW for on-plate localization and plate geometry inference.

III. METHODOLOGY

A. Notations and assumptions

In this paper we are considering a mobile unit transporting an acoustic emitter-receiver pair on a metal plate. At every time step k, the emitter sends an acoustic pulse s k b (t). Under mild assumptions of the acoustic properties of the ultrasonic probing system, we can model the received transducer signal s k r (t) as a linear system through convolution:

s r (t) = x(t) * s b (t) (1)
with x(t) the environments impulse response and s r (t) hopefully containing information about reflections on the plate edges. The superscript k is omitted when unambiguous. We assume that the plate is an homogeneous material (steel, aluminium,...) of a constant thickness. We assume the plate to be a convex polygon, but we do not assume it to be rectangular, even though this is to be expected in any industrial case. Because the edges are linear, we will also assume that only the orthogonal reflection of the signal is picked up by the receiver. Any other reflection may lead to secondary echoes after bouncing from several edges but these are neglected in our work.

Out of the reflected signal s k r (t), we will denote as {r i (k), i = 1 . . . n(k)} the set of detected ranges to the plate edges, assuming that not all the edges are detected all the time and that spurious detection may be included in the list.

Additionally, because we assume the transducers are carried by a mobile device on the plate, we assume that some level of odometry is available to estimate the displacement between measurements. Because the inspection crawlers we are considering are moving on vertical structures (boat hulls, storage tanks), they can easily embed an accelerometer from which the crawler heading with respect to gravity can be observed with a good enough precision. Also, by construction, to avoid falling from the inspected structures, these crawlers have to have an extremely good adherence and incur very little slippage. We can hence assume the odometry locally very precise. These two assumptions (known heading and very precise odometry) are an important distinction with respect to the room reconstruction works.

In the next section we will describe first how these echoes can be robustly extracted from the signal and then how they can be integrated with the proprioceptive measurements into a localization and mapping framework.

B. Echo detection

One of the key enabling technologies for the on-plate SLAM approach is the robust and accurate detection of ultrasonic echoes in the metallic plate. Due to the physics of the sound propagation in the metallic plate, only signals with limited bandwidth can be used for probing the underlying material. This limited bandwidth constraint complicates accurate time-of-arrival estimation for the reflections in the medium, removing classical signal processing techniques such as matched filtering from the set of applicable techniques. We apply a super-resolution technique based on L 1 -regularized least-squares and the concept of sparsity to solve this estimation problem. The underlying hypothesis of this approach is that most of the metal plate does not reflect the ultrasonic signals, and that only discontinuities in the medium cause reflections to occur. In that case, the environment has only a limited number of reflections, and can be considered sparse.

Under the sparsity assumption of the environment, we can consider x(t) from eq. 1 to be sparse. After discretization of the problem, this allows us to cast the following L 1 -regularized least squares problem:

min x |D • x -s r | 2 2 + λ • |x| 1 (2)
with x the discretized vector-representation of the impulse response x(t) of size [n t × 1], s r the discretized received acoustic signal of size [n t × 1], where n t is the number of time-samples used in the discretization process. The vector norms |a| p denote the L p -norm of vector-object a, with the well-known fact that the L 1 -norm promotes sparsity of the least-squares solution. The matrix D is called a dictionary matrix of size [n t ×n t ], containing time-shifted copies of the emitted echo signal s b (t):

D =       s b (t) 0 0 ... 0 0 s b (t) 0 ... 0 0 0 s b (t) ... 0 ... ... ... ... ... 0 0 0 ... s b (t)       (3)
This minimization problem is a convex problem which can easily be solved using open-source toolboxes such as the Matlab CVX toolbox [START_REF] Grant | Cvx: Matlab software for disciplined convex programming[END_REF].

C. From echoes to edge hypotheses

A single echo r i (k) provides information about the distance to an object. As will be described later, this information is useful to refine the parameters of the edge the echo reflected from. However, this is too ambiguous to initialize a new edge: it could be any tangent to the circle centered on the current pose with radius r i (k).

In [START_REF] Kreković | Echoslam: Simultaneous localization and mapping with acoustic echoes[END_REF] on the other hand, the authors show that 2 echoes taken at different known position are sufficient to define two line hypotheses: there are only two lines tangent to the two circles defined by these radius and centered on the known poses. A third echo can be used to disambiguate between the two lines if the trajectory is not parallel to the line. Beside requiring three perception steps and the corresponding delay, using three echoes requires considering all the O(n 3 ) 3tuples.

Instead, we decided to use only the last two echoes and deal with the multiple hypotheses in a later step. There are several arguments for this choice: first the disambiguation with 3 echoes is only useful when the path is not parallel to the edge, which is actually a very common case in a practical deployment where the robot use its accelerometer to realize vertical or horizontal transects. Second, even on a non-parallel path, a clear disambiguation require a large translation between the second and third echo to compensate for the perception noise. Last, when moving on a direction normal to the edge, two measurements are sufficient for a unique solution. Both cases are illustrated in fig. 2.

Because of our combinatorial generation of edge hypotheses, a method is required to identify inconsistent edges and remove them. To this end, we take advantage of the hypothesis that only orthogonal reflections are reaching the receiver. Hence, with the assumption of a convex plate, if we consider two edges L 1 and L 2 observed from pose P and P i the projection of P on L i , then L 2 cannot intersect the If one of this condition is true, then both edges cannot exist simultaneously on a convex plate because observing one edge would require the ultrasonic waves to cross the other one, which is not physically possible. In the following, this test will be denoted as the consistency check between two edge hypotheses.

D. FastSLAM integration

FastSLAM [START_REF] Montemerlo | Fastslam: A factored solution to the simultaneous localization and mapping problem[END_REF] is a solution to the simultaneous localization and mapping problem which is particularly relevant in the context of uncertain data association. In essence, FastSLAM relies on a particle filter in the localization space and every particle holds an hypothesis of the map that can be inferred from the trajectory of this particle. In order to integrate FastSLAM into our localization and plate geometry inference problem, this section will define the particle state, the chosen map representation, the particle initialization strategy, how the internal maps are updated and the evaluation of a particle based on the measured echoes.

1) FastSLAM state: Because of our assumption of a known orientation, we only consider the sensor coordinates (x, y) on the plate as system state. Orientation could be added to the state without significant changes to the framework.

As a map representation, we decided to use a list of infinite lines represented by an angle and offset. Several options were considered and ultimately rejected. Inspired from traditional occupancy grids [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF], we first considered storing a grid map where every cell would represent either the likelihood of being on an edge or the likelihood of being inside the plate. Both representations work well in a mapping only framework but in a SLAM context they lack robustness. Using the edge likelihood presents a challenge for removing outliers and artifacts. Using the likelihood of being inside the plate works well when only using the first echo (i.e. the closest edge) but does not scale to using echoes from multiple edges. Using a parameterization as an oriented rectangle was also evaluated but rejected because first, it only works for rectangles and second, the estimator tended to forget about one side of the rectangle when only the other one was observed for too long. In the end, estimating independent lines offers more flexibility and provides sufficient information about the plate geometry. In practice, a line is represented by a pair (θ, b) that define the line equation:

cos(θ) • x + sin(θ) • y + b = 0 (4)
A particle P i can then be described as a trajectory hypothesis associated with an estimated map defined as a set of lines:

P i = [(x i , y i ) i=1...k , {(θ j , b j ), j = 1 . . . n i }] (5)
2) FastSLAM initialization: Most of the earlier papers [START_REF] Kreković | Echoslam: Simultaneous localization and mapping with acoustic echoes[END_REF], [START_REF] Krekovic | Look, no beacons! optimal all-in-one echoslam[END_REF] on room reconstruction do not discuss the initialization of the Bayesian filter, which is an important element of a practical implementation. In our case, we rely on the edge hypotheses generated from section III-C. After waiting for the second set of range measurements, we estimate a number of edge hypotheses from which we extract maximally consistent sets. These sets are built using a dynamic programming approach inspired from JCBB/JCDA [START_REF] Neira | Data association in stochastic mapping using the joint compatibility test[END_REF] which will be omitted here for the sake of page limits. As a result, every particle is initialized by sampling around the zero position and randomly selecting a consistent set of edges from the maximally consistent sets. Including all the edge hypotheses in all the particles would be a viable alternative given the outlier removal decision described below. It would however add more ambiguity than necessary in the estimation.

3) Particle evaluation: Given {r i (k), i = 1 . . . n(k)} the set of echoes measured at time k, we evaluate a particle based on its ability to explain the measurements. For a line L i,j = (θ i,j , b i,j ) in particle P i , the expected measured range is

d i,j = | cos(θ i,j ) • x i + sin(θ i,j ) • y i + b i,j |
From this range, the likelihood of measurement r i (k) given L i,j is expected to follow a Gaussian distribution centered on d i,j , with a standard deviation consistent with the echo detection uncertainty. If the highest likelihood over all the lines is lower than a threshold, then the measurement is considered "unexplained" and allocated a low probability P 0 .

For a complete set of echoes, the evaluation of a particle will then be the product of the likelihood of the n(k) likelihoods of the independent measurements:

L(P i ) = n(k) i=1 max P 0 , max j (P (r i (k) | P i L i,j ))
In this formulation, the product has n(k) terms which make the evaluation of different particles comparable. Performing the product over the set of lines of each particle would result in a varying number of terms, which is incompatible with importance sampling.

Note that for the purpose of the map update, if a measurement is considered explained, the index of the line leading to highest likelihood is recorded as j k,i . The list of unexplained measurement is also stored to create additional line hypotheses.

4) Map update:

The map update stage has three purposes, first for "explained" measurements, the associated line needs to be updated to account for the new piece of information. Second, the "unexplained" measurements are used to create additional line hypotheses. Third, consistency checks are used to eliminate line hypotheses which are no longer supported by the observations. a) Line update: Knowing a measurement r i (k) and its associated line in particle P i , L i,j k,i , one needs to update the corresponding (θ j , b j ). The definition of d i,j suggest the use of an Extended Kalman Filter. Even though such a filter is feasible, the individual measurements can often be explained equivalently by changing θ j or b j , which prevents a precise convergence of the filter. To sidestep this issue, an alternative is to keep a record of the sensor poses and ranges associated with this line and use a non-linear minimization of the following cost function:

C j (θ, b) = (6) k s=1 (cos(θ) • x i (s) + sin(θ) • y i (s) + b) 2 -r(s) 2
Although much more expensive, this approach converges without bias even when the measurements are taken on a monotonic walk along the edge.

b) New lines: For unexplained measurements, we need to combine them with the measurements from the previous time step to create edge hypotheses (sec. III-C). Every pair made of one previous measurement and a new unexplained measurement is used to generate one or two edge hypotheses. All the previous measurements are used for additional robustness against incorrectly associated measurements at the previous time step. All these hypotheses are added to the line set of the current particles. c) Map clean-up: In a final stage of the map update, all the lines in a particle line set are checked for pairwise consistency. If two lines are deemed inconsistent and one of them has been observed (i.e. associated with a measurement) significantly more than the other since its creation, then the least observed line is marked for deletion. In our implementation, this criterion is defined as being observed two additional times. In this stage we also mark for deletion lines that were not re-observed enough since their creation. For instance, we delete a line that is more than 10 stepold but has been observed less than 3 times. The simple delayed outlier removal described in this paragraph is what makes possible the generation of many edge hypotheses from the set of detected echoes. This is critical since many of these hypotheses cannot be ruled out until the robot moves significantly and changes direction. In comparison, in [START_REF] Kreković | Echoslam: Simultaneous localization and mapping with acoustic echoes[END_REF] the authors made the hypotheses of a random walk of the agent, which is a much more informative path, but also much less realistic for a robotic crawler.

IV. RESULTS

A. Experimental setup

At the time of this writing, the experimental setup is not yet integrated on a robotic crawler. Instead, in order to test our framework in a setup as close as possible from reality, we used a pair of emitter-receiver piezo-electric transducer on two different aluminium plates (plate 1: 470x470x5mm and plate 2: 600x450x6mm) and moved them by hand on the vertices of a regular grid with 30 to 40mm spacing. At every position, the response to 10 ultrasonic scans were averaged and recorded with their acquisition position. A scan in this context is the emission of two periods of a 100kHz sine wave and the recording of 400µs of analog signal at 1Msample/second. To simulate a sweep of a plate by a robotic crawler, the corresponding sequence of measurements is selected from the database and presented to the SLAM framework, with the theoretic displacement between grid cells used as odometry. We believe that, apart from the size of the plates which is smaller than real ones, this setup is a correct representation of a real system in terms of signal quality and overall input precision.

B. Echo detection

Figure 4 illustrates the echo detection process. Panel a) shows the emitted signal, which is a 2-cycle burst at 100kHz. The reflected signal from the plate structure can be seen in panel b). It shows the pickup of the emission in the beginning of the signal, and a reflection due to coupling mismatches. Then, a series of echoes is apparent. Using the L 1 -based approach, the impulse response can be reconstructed (see panel b, orange trace). The method reconstructs the major echo components of the impulse response, while ensuring that the solution is sparse. Panel c) shows the reconstruction error. The fact that this residual is not zero can be explained by the fact that the direction-dependent filtering of the transducer is not taken into account during the reconstruction phase. However, as the remaining experiments show, this crude approach to the echo detection process is sufficient for the plate-estimation algorithm. C. Edge hypotheses Fig. 5 gives an example of computed edges from the combination of the two first measurements on plate 2. As can be observed, the number of edge hypotheses is relatively large but a few of them are definitely good estimates of the real top and bottom edges. At this position, lateral edges are not observed correctly. The level of uncertainty in the putative edges is what calls for the FastSLAM framework and the delayed decision the validity of the edges.

D. Mapping performance

Figure 6 shows one anecdotal example of a mapping run on plate 2 by displaying the highest-ranked particle (out of 16 in this case), its line set and its trajectory. The reconstructed path is clearly visible and the stability of the estimate can be observed when performing a second sweep of the path (iterations 225 and 239). The convergence of the estimated lines to the true plate outline can also be observed, with a final error in the order of 2cm on the line offsets b (see eq. 4). Even though this run has a lower precision than the average case, 2cm is still a very acceptable precision given that the sensors are positioned by hand and that the wave length in this plate is approximately 3cm. As a final evaluation, we evaluated the precision and repeatability of our approach through 100 repetitions in three scenarios. Scenario 1 and 2 consisted of a lawn-mover path through the plate 1 and 2. The length of the sweeping path being 81 steps and 108 steps respectively. In a third scenario, we used plate 1 and simulated localization and plate geometry estimation over a random walk. In the three cases, the real measurements described above were used. However, a single database of measurements was used for every repetition. The table below presents a summary of the resulting precision at different steps in the process. It is clear that the initial estimates after observing only half of the plate are still very uncertain for all the scenarios, but at the end of the simulation all lines from plate 1 are estimated with a precision consistent with the signal wavelength. The poorer performance on plate 2 are due to the somewhat larger size of the plate and to the presence of artificial defects that are acting as reflectors and creating detection artifacts.

V. CONCLUSIONS

This paper presented a proof-of-concept for a localization and plate geometry inference framework for a magnetic crawler performing inspection of structures assembled out of metal plates. The results show that there is enough information in the reflected signals to achieve a good localization and mapping precision as long as a sufficient coverage of the plate is performed. The next steps will be to embed this framework on a robotic platform as presented in fig. 1, test on larger plate, improve the overall system robustness and consider an active sensing strategy to recover the plate geometry faster and even more reliably.

Fig. 1 .

 1 Fig. 1. Typical inspection conditions for magnetic crawlers operating on a structure assembled out of welded steel plates [Source: RoboPlanet]

Fig. 3 .

 3 Fig. 3. Experimental setup for metal plate 1. The emitter and receiver transducer can be seen on the bottom right of the figure.

Fig. 4 .

 4 Fig. 4. Illustration of the echo detection process. Panel a) shows the emitted signal. Panel b) shows the echo signal (blue trace) and the reconstructed impulse response (orange trace). Panel c) shows the reconstruction residual. The remaining energy in that residual can be explained by the mismatch in the used signal model D and the real transducer.

Fig. 5 .

 5 Fig. 5. Edge hypotheses generated on plate 2 from detected echoes r(1) (blue dot) and r(2) (purple dot), without consistency check. The purple rectangle is the outline of the true plate.

Fig. 6 .

 6 Fig. 6. Evolution of the map representation for the highest-ranked particle at step 1, 12, 36, 124, 225 and 239. The purple frame represent the outline of the true plate. The red dot is the current estimated sensor pose and the red line is the history of its estimated trajectory. The darkness of the line is proportional to the number of times they have been observed.

Scenario

  Angle error [rad] Offset error [m] Scenario 1, step 50 0.017 ± 0.057 0.004 ± 0.032 Scenario 2, step 50 -0.029 ± 0.154 -0.061 ± 0.199 Scenario 3, step 50 -0.008 ± 0.084 -0.031 ± 0.133 Scenario 1, step 239 0.005 ± 0.020 0.002 ± 0.010 Scenario 2, step 239 -0.007 ± 0.091 -0.015 ± 0.126 Scenario 3, step 239 -0.001 ± 0.040 0.002 ± 0.045

This work was supported by the EC under grant agreement 123456


J L Le Bras

M Le Goff

B Deschamps

F Le Pennec

A Péden

D Bourreau
email: daniel.bourreau@enst-bretagne.fr

S Toutain

HIGH EFFICIENCY PLANAR FREQUENCY MULTIPLIER AT MILLIMETRE WA VES AND QUASI-OPTICAL INVESTIGATIONS FOR (SUB) MILLIMETRE APPLICATIONS

This paper presents the design and measurement of a hybrid varactor frequency multiplier in the millimetre wave range and the first investigations on a quasi-optical frequency doubler. For the planar circuit, a multilayer technology was used for the synthesis of the source and load impedances of the varactor. The conversion tosses are about 12 dB at 62 GHz with an output power of 10.8 dBm.

The measurement test bench and results on passive circuits in quasi-optical technology are presented. A quasi-optical doubler has also been investigated and designed and the first results are presented.

INTRODUCTION

For sub-millimetre wave applications, like FIRST and SOFIA programs [I], local oscillators at very high frequencies (up to 1 THz) are required for heterodyne receivers. A solution for providing power with low phase noise at sub-millimetre waves associates a low frequency oscillator with one -or more -frequency multipliers. This high frequency oscillator must have sufficient power for efficient driving of the down converter. Waveguide techniques are still widely used for millimetre and sub-millimetre radiometers, but planar circuits (hybrid or MMIC) could allow the drawbacks of waveguide whisker contacted diode multipliers (losses, fragility• ... ) to be overcome. However, no tuning can be done on planar multipliers and accurate modelling of the circuit is required. Quasi-optical solutions could be an interesting choice in the near future.

Il. PLANAR CIRCUIT DESIGN

The source and load impedances and the working conditions of the multiplier are initialised with the Penfield and Rafuse equations [START_REF] Penfield | Varactor Applications[END_REF]. Non-linear simulations are then used to optimise the results and calculate the output power at fondamental and harmonie frequencies (figure 1 ).

Realisations at lower frequencies and sensitivity simulations give important information to control each step in the design procedure such as impedance terminations at fondamental and second harmonie frequencies (for a doubler). For the 28/56 GHz doubler, radial stubs are used to short-circuit the second harmonie at the input and the fondamental frequency at the output in a large frequency band. These stubs are positioned very close to the diode : this is very important to improve the conversion efficiency.

Classical planar technologies such as microstrip or coplanar technologies lead to very low and very high impedance transmission line sections for the diode matching, inducing high losses and then reducing the efficiency. A technology [START_REF] Rius | lntegration of Various Types of Compensated Dielectric Bridges for Millimeter Coplanar Applications[END_REF] with microstrip and multilayer coplanar lines (figure 2) was adopted to combine the advantages of each one. To obtain very low characteristic impedances, a small substrate thickness is used. But at the same time, high values are required. The multilayer technology detailed in figures 2 and 3 allows us to approximately retain the same central conductor width throughout the circuit, cutting the ground plane under the inductive transmission line sections. Classical coplanar lines are used at the accesses for coaxial probe measurements. The doubler circuit is fabricated on a 635 µm Alumina substrate with two metal layers and one dielectric layer. The surface is 6.7 mm x 9.7 mm.

Zg 2jo = 0

Zl 2 1 0 Zg Jjo Zl 3.fo z... z ...

III. RESUL TS FOR THE PLANAR DOUBLER

Experimental results (figure 4) show that the optimal input frequency is 30.8 GHz for maximum output power (10.8 dBm) and efficiency (10.3 dB). The fondamental rejection is 25 dB. But at the 28.25 GHz input frequency, the fundamental rejection is optimal (-48 dB) as expected in simulations.

The electromagnetic simulation (2D and 3D) of the matching circuits was compared with HP-MDS library models used for the multiplier design. Differences in the computed S parameters, that is, the source and Joad impedances, can explain the frequency shift of the optimal output power. The dielectric permittivity of the thin layer should also be extracted for better accuracy.

No instabilities are observed at any frequency or power sweep. To have a good stability forecast many simulations and tests were performed.

IV. QUASI-OPTICAL TECHNIQUE

Recent developments give a solid state response for power capabilities at millimetre and sub-millimetre waves [START_REF] Steup | A Quasi-Optical Doubler-Array", !nt[END_REF]. Fig. 5 illustrates an example of a quasi-optical circuit. Passive components were studied for accurate modelling [START_REF] Le Bras | Quasi-Optical Circuit Measurement Method in the W Band[END_REF]. The grid array can be an amplifier (output at fO) or a multiplier (output at n.fO).

Active planar arrays allow high level output power by recombination in free space with increased reliability. The matching of the active device (diode, transistor) can be done on the active grid, or can be carried out by extemal components such as a polarizer which works like a tuning short circuit and a dielectric slab like an adaptator. For a better input and output isolation, the electric fields at the two ports are cross-polarized. lt is also important for a multiplier to shunt the second harmonie at the input and the fundamental at the output with two separate circuits.

To design a millimetre doubler, the hybrid active grid of figures 6 and 7 with 69 cells, was chosen because of its simplicity ( [ 4 J). The unit cell consists of two quarter wave dipole antennas in quadrature polarisation at the fundamental ( 38 GHz) and at the second harmonie (76 GHz). These antennas are connected via a V aractor diode which is the non-linear component. The active grid (i.e. the multiplier ) is extemally matched with tuning slabs and polarizers which can be shifted during the measurement for optimal performances ( output power, conversion efficiency, ... ). The multiplier grid is printed on 210 µm substrate (i:;,=3.5).

Tuning Slab

Input Bearn ( fo)

A. Quasi-Optical test bench and modelling

For quasi-optical element characterisation and modelling, two test benches were carried out with Gaussian Optic Lens Antennas (GOLA) in the Ka band and in the W band. These benches were validated for S parameter measurements of passive quasi-optical elements such as polarizer, filter and dielectric slab. Simulated and measured S parameters show very good agreement. The complex permittivity of dielectric slabs is also extracted with high precision [START_REF] Legoff | Ka Band Quasi-Optical Test bench using Focusing Homs[END_REF].

Using models or measured S parameters, simulations of simple quasi-optical fonctions (cascaded polarizer and slab for example) in classical CAD software (such as HP MDS) were carried out [START_REF] Le Bras | Quasi-Optical Circuit Measurement Method in the W Band[END_REF]. Non-linear simulations of the whole quasi-optical frequency multiplier were investigated and are still under study. The precise modelling of the component access through the antenna (dipole, ... ) is very important and rather difficult because it has to take into account ail the modes in the structure, the cross-polarisation and the multiple accesses (quasi optical and active components ).

B. Quasi-Optical multiplier measurement

With a Ka band GOLA, a W band GOLA and a Ka band power amplifier (lOW), the active array with input and output polarizers was tested with a 38 GHz input frequency.

In a first experiment, only 9 diodes were mounted on the array. But with so few diodes, the performances are not very attractive. The conversion efficiency is iess than -40d8 at the 76 GHz output frequency. To explain these poor results, it must be pointed out that the 1. 7 cm active area diameter with 9 diodes on the grid is too small in comparison with the 8 cm waist diameter of the input Gaussian beam. So, too little power is intercepted by the active area that is, by the 9 dipole antennas connected to the diodes. Most of the incident power is not transferred to the active device and the efficiency is very low.

The second harmonie output power was measured for different positions of the input polarizer in relation to the active grid. Fig. 8 gives the normalised output power evolution and shows some periodicity. This result is very important because the positioning of the tuning elements must be very precise : there is only 700 µm between the maximal output power position and the minimal output power position. However, this behaviour was observed in preliminary non-linear simulations of the whole multiplier with very simple modelling (an inductance) of the dipole antenna.

VI. CONCLUSION

The design and measurement of a planar hybrid frequency doubler have been presented. The diode matching circuits make use of a coplanar and microstrip multilayer technology for high and low impedance synthesis without discontinuity of the central strip. With a 9.3 % (-10.3 dB) maximum conversion efficiency, it is a candidate for a state of the art of hybrid diode multipliers at this frequency. However, frequency shift was observed.

Quasi-optical technology is very attractive for combining solid state components, tuning circuits and power. Unlike passive circuits, the simulation of an active grid is quite difficult to achieve. Developments are in progress and a W band quasi-optical doubler will be tested in the near future.

Fig. 1 .

 1 Fig. 1. lmpedance optimization at each port.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Multilayer technology process on an Alumina substrate with two metal layers and one dielectric layer.

4 .

 4 Pout(@fO) • -••••• Pout(@2f0) P(dBm)-35 *"--'-� -'--ll--,..�--l��'--l -��-1-��-Output power at fondamental and second harmonie for P in = 20 dBm and V=-9 V.

Fig. 5 .

 5 Fig. 5. Quasi-optical frequency multiplier.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Photograph of the active grid (dimension in cm).

Fig. 8 .

 8 Fig. 8. Normalized output power measurement of the quasi-optical multiplier for ditlèrent positions of the input polarizer in relation to the active grid.

ACKNOWLEDGEMENTS

The authors would like to thank Eric Rius (LEST Brest France) for his planar circuit design expertise.


Alexandre Devaux

Mathieu Brédif

Realtime Projective Multi-Texturing of Pointclouds and Meshes for a Realistic Street-View Web Navigation

Keywords: Image Based Rendering, Projective Texturing, Streetview, WebGL, Point Based Rendering, GIS Concepts:, Information systems → Web applications, •Humancentered computing → Geographic visualization, •Computing methodologies → Image-based rendering, Point-based models

Figure 1: Rendering using mesh and pointCloud. a) Projective texturing on simplified planar geometry (point of view one meter distant from original photography position) b) Original Point cloud c) Projective texturing on Point Cloud. d) Hybrid rendering.

Introduction

City models visualization attained the last decade a high level of quality with realistic fly-throughs even on web platforms. However, displaying realistic virtual environements of real-world places at a pedestrian level is still a challenging task. On the one hand, accurate 3D geometry and calibrated images may be acquired and used, within an expensive process, to model a high quality environment with a sufficiently high level of details. On the other hand, with a more approximate geometry, different image-based rendering (IBR) techniques have coped to create satisfying renderings without requiring a detailed modelisation.

In the early days, Light Field Rendering or Unstructured Lumigraph Rendering [START_REF] Buehler | Unstructured lumigraph rendering[END_REF] allowed the real-time visualization of unstructured sets of cameras with or without relying on an approximate model (so-called geometry proxy). More recently, some * E-Mail for all authors:firstname.lastname@ign.fr impressive results [START_REF] Chaurasia | Depth synthesis and local warps for plausible image-based navigation[END_REF]] have been achieved to provide plausible free-viewpoint navigation using local warps with few input calibrated images, which can cope with unreliable geometry. Google Street View TM uses a sparse set of panoramic images. Transitions between captured viewpoints employ cross-fading and geometry-aware warping to approximate the expected optical flow. Microsoft Photosynth TM displays an unstructured collection of photographs in the reconstructed spatial layout and applies image-space transformations and blending transitions. Despite their relative simplicity, these systems create reasonably compelling 3D experiences and they are quite suited for web platforms. The biggest remaining drawback is that such systems typically restrict viewers to be very close to one of the capture points. They also tend to remove small depth details such as poles or other light urban furnitures in the targeted street-view application.

The point-based representation is getting more popular since two decades, thanks to the progress of Structure From Motion, lidar technology, and light RGB-D system such as the kinect. The litterature on point-based rendering is therefore pretty dense. The most trivial way to render point clouds is to use one color per point and project it to the screen. It usually creates holes and low color resolution even for dense point clouds. The other classical approach is to mesh the point cloud using [START_REF] Kazhdan | Poisson surface reconstruction[END_REF]] or 3D Delaunay reconstruction offering a hole-free model. Then all images are generally merged into a texture atlas to get a full-resolution unified texture. The usual drawback is to add artefacts to the original geometry and it lacks small details. Potree, a popular WebGL pointcloud renderer, uses an octree to spatialize the data for fast streaming and add compelling rendering algorithms such as Eye-Dome-Lighting (EDL) to improve depth perception without the need of normals on points. Combining point cloud with mesh is less studied. [Pagés et al. 2015] explored the mixing of meshes with splats creating hybrid renderings but the methodology requires heavy pre-processing such as Poisson reconstruction and geometry shape segmentation. Some experiments on the perception of visual artifacts in imagebased rendering have been conducted by [START_REF] Vangorp | Perception of visual artifacts in image-based rendering of fac ¸ades[END_REF]] in urban areas and show for example that when cross-fading between panoramic images, shorter transition durations are preferred. Another side effect we discovered is that the eye is more distracted by deformed textured meshes even though geometrically closer to reality than by very approximate planar geometries inducing more parallax effects (fig. 3). Creating a simplified mesh from the pointcloud using the previously discussed methods requires very refined algorithms working on high quality point clouds to get acceptable results inducing also heavier data to be transmitted and rendered. Another artefact often encountered while creating meshes is the loss of small details. In our case, we found it very hard to create meshes that could keep geometry such as vegetation, poles, balconies, and, more generally, urban furnitures. This is problematic as the user moves through the streets, as it means that the texture will be projected on an incorrect geometry, typically flattened on the ground or fac ¸ades, inducing an unpleasant rendering (fig. 1.a). For this reason, a view-dependant projection of oriented images on a simplified mesh (composed by few quads for the facade and the road such as LOD1 building model) is a very efficient solution as a base rendering. Its simplicity makes it extremely fast to stream and very light to render not to mention its increasing disponibility as OpenData. Another important feature is that it allows easy daily update as new data is acquired it is instantly visible without pre-processing. As in [Brédif 2013], we compute the mesh from 2D building footprints and a low resolution height-map (or Digital Terrain Model, DTM). We extend this method to the web environment by streaming the DTM and the 3D gutter linestrings only and reconstructing the simplified mesh online (in the client) using a constrained Delaunay triangulation.

Our main contributions are (i) to add geometric details and refine occlusion handling, by creating detailed depth maps using the point cloud (ii) to use realtime projective texturing on both the mesh and the point cloud to create a compelling view-dependent hybrid rendering that does not rely on any pointcloud or texture atlas preprocessing. The proposed method thus contributes on reducing parallax distortions, blurring, ghosting and popping errors that usually distort the scene appearance through IBR in a web-friendly manner as the lightweight base mesh and a low-res projective texture are streamed almost instantly, producing a reasonnable view (fig. 1

Overview of the method

Working with mobile-mapping acquired data, the application deals with oriented images taken through the vehicle trajectory, hence a small number of images usually see the same spot. This specificity makes projective texturing an adequate choice. The basic application pipeline is as follow (fig. 2): 1. Get the building footprints around the viewing position v and triangulate to create basic geometry proxy (mesh)

2. Get the projective textures {t k } (images and calibration) and lidar point cloud {c l } close to v.

3. Initialize the depth buffer of each texture t k with a depth-pass on the base mesh and update it as new chunks c l are received.

4. Hybrid (mesh + point cloud) projective texturing with texture distortion and occlusion handling support.

Point Cloud Visualization Optimization

Raw lidar data comes with a lot of information, such as 3D position, reflectance, number of echoes, etc but does not bring any color information. Offline colorization [START_REF] Waechter | Let there be color! Large-scale texturing of 3D reconstructions[END_REF]] can provide good results using oriented images with method classically applied to meshes. Though it usually creates one colour per point or generate huge texture atlases.

Our method does not need any pre-processing to get very pleasing point clouds renderings. We rely on the same online projective texturing we use for the mesh hence mutualizing textures for fast web transmission. Another point to justify the whole methodology of projective texturing here is that the rendering is perfectly equal to the original images when viewed from their viewpoint location, and degrades gracefully in-between viewpoint locations. Pre-processed textured meshes or point clouds are optimized to provide the best global representation of acquired data from any point of view which necessarily distort reality.

No matter the density of the point cloud, setting a unique color per point will always create a subsampled look compared to the original image. We experimented the creation of splats directly in the GPU in a single pass using barycentric coordinates. When receiving points from the server, we create for each point a triangle with colocated vertices but with differing barycentric coordinates so that splats may be created in the vertex shader (as geometry shaders are unavailable in WebGL) and change their projected geometry (as ellipses) in the fragment shader.

Alternatively to using splats, we propose a perspective-correct projective texture mapping of point sprites. Similarly, we approximate the point rendering by a set of planar patches, but instead of rasterizing one or two coplanar triangles per point, a point sprite (i.e. a screenspace square) is used, minimizing the primitive setup and rasterization workloads. The 3D planar support is left unchanged as the plane going through the 3D point P and orthogonal to its normal N . The main idea is to compute the homography H induced by this plane from the synthesized view to the input view [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Using projective geometry notations, we denote the respective input and synthesized view projections as p = M (P -C) and p = M (P -C ) where C, C are the viewpoints, M, M are the view-projection 3x3 matrices and p, p are the image coordinates in screen-pixel and texture coordinates respectively:

H = M I + (C -C )N T N T (P -C) M -1
Implementation-wise, this homography may be computed in the vertex shader, passed as a mat3 varying and evaluated as p tex = Hpscreen in the fragment shader. The vertex shader computation is performed as H = M + E N T N T (P -C) M -1 with:

• View uniforms M -1 and C

• Per-texture uniforms M and E = M (C -C ) • Per-point attributes P and N (N = P -C by default)

Figure 4 shows the large difference of rendering quality when computing one color per point versus texturing the whole point surface. Raw lidar point clouds do not generally come with normals. If unavailable, we could compute a screen-space normal map using a depth-prepass and a normal-map lookup, which would yield some side-benefits but would slow down the rendering on light GPU configurations with an extra pass. We propose a much simpler solution leaving the point sprites at their original behaviour: always facing the camera with an eye-vector normal approximation. Applicationwise, this introduces a screen-space dilation of the point cloud support that fills holes and nicely convey the cylindricity on poles and other thin objects typically acquired by a single lidar scan line (where a normal has no real meaning). Note in figure 4 how high resolution projective textures handle well sparse point clouds.

Occlusions

A classic artefact when navigating through Street View applications are bad occlusions handling. Namely, texturing part of the geometry with an irrelevant texture. Examples are (i) disocclusions (e.g. at crossroads, fig. 6), where a disoccluded geometry is textured by a viewpoint from which it is not visible and (ii) unmodeled geometry (e.g. poles projected on the ground, fig. 1.a).

Per-texture Depth Map

We adapted the shadow mapping-like technique from mesh-only [Brédif 2013] to hybrid rendering (mesh+pointcloud). This is performed by computing an off-line rendering of the scene from each texture position, creating a depth texture that is pixel-wise aligned with the (undistorted) input color image (fig. 5). While rendering, a simple depth (shadow) test is performed between a texturing depthmap lookup and the depth relative to the sampled texture viewpoint. This technique is relatively cheap as it only costs a onetime depth-pass during the initialization phase and a comparison against a texture sample (i.e. a shadow test) in the rendering phase.

Hybrid View-Dependent Multi-Texture Mapping

Compared to other View-Dependent Multi-Texture Mapping approaches, we propose two extensions. First, to address the web context, the projective texture calibrations (position, rotation and projection matrices) are stored in a PostgreSQL/PostGIS database with a Node.js frontend, enabling fast queries for nearest textures from the current view position and pre-fetching based on query regions. Second, the hybrid rendering of an approximate mesh with a more accurate but incomplete pointcloud renders well the added pointcloud protrusion details (e.g. balconies, poles...) but not the depression details where the approximate mesh was overestimated (windows, doors, insets...). To cope with this issue, the mesh is rendered with an offset (e.g. 1m) that lets points belonging to small depressions pass the depth test.

Results and Limitations

The use of both a simplified mesh planar geometry and a point cloud offers a compelling visualization. As we said earlier they are combining in such a way to enhance realism. However, as we use simplified planar geometry computed from external databases (building footprints), its geometry can have an important shift with mobile mapping geometry. One way to deal with that is to bundle adjust the mesh geometry using trajectory information if available. As we load oriented images, they have absolute position informations and usually comes with height above ground information. Hence we can adjust the local altitude of the mesh with this information creating a local Digital Terrain Model (DTM), again, without preprocessing. Another way is to analyze the point cloud estimating the minimal height of the points gives already enough information to adjust the height of the mesh. Those are classical problematics when dealing with different sources of geospatial data.

Conclusion

The solution proposed in this paper contributes to web immersive navigation system enhancing its realism while keeping precise data representation. It can be improved in different ways. Even if there is a need to keep access to raw lidar point clouds for professional usage, pre-processing it for visualization aspect will lighten it for streaming as splats precomputing can save a lot of bandwith in streaming and in GPU memory. We are currently working on improving the streaming of the data, investigating with Shape Resource Container (SRC) [START_REF] Limper | Src -a streamable format for generalized web-based 3d data transmission[END_REF]] which offers great optimization for textured mesh level of detailed transmission and could also work for splats.

Figure 2 :

 2 Figure 2: Overall application pipeline

Figure 3 :

 3 Figure 3: Comparaison of support for projective texturing. Poisson reconstructed mesh from Lidar (left) and simple planar geometry (right)

  .a), while local pointcloud chunks and higher resolution textures are streamed, updating continuously the visualization (fig. 1.d).

Figure 4 :

 4 Figure 4: Comparaison of point cloud rendering. (left) Basic point sprites, 1 color per point. (right) Projective textured point sprites

Figure 5 :

 5 Figure 5: Depth map for one projective texture (right) mixed with original imagery from other projective textures (left).

Figure 6 :

 6 Figure 6: View-dependent texturing with (bottom) and without (top) occlusion handling.

Acknowledgements

The implementation of this work and a sample dataset are available with an open source license at https://github.com/iTowns