Mihai Bogdan Luca
email: lucami@mta.ro

Alexandru Serbanescu

Stéphane Azou
email: stephane.azou@univ-brest.fr

Gilles Burel
email: gilles.burel@univ-brest.fr

A NEW COMPRESSION METHOD USING A CHAOTIC SYMBOLIC APPROACH

Keywords: compression algorithm, symbolic dynamics, probabilistc Bernoulli generator, binary initial condition

In this article, a new compression method is investigated using a particular chaotic modulation type. The symbolic approach presented associates to every informational sequence a trajectory in the state space of the chaotic generator. We are also introducing a new type of chaotic generator adapted to the probability distribution of the informational sequence, and we prove that using the described generator the theoretical compression performances attain the optimal entropy compression. Finally we confirm the theoretical analysis with performance tests for different sequences with different symbolic probabilities.

I. INTRODUCTION

The problem of loss less data compression has received an increased attention both by the soft making companies but also by free-lancer researchers who have developed various types of algorithms, starting from static ones to the more powerful adaptive ones [START_REF]Arithmetic Coding and Decoding Methods and Related Systems[END_REF].

In opposition with the data compression algorithms, the use of chaotic signals is relatively new in the area of telecommunications and few studies have been made with the application of chaotic signals and generators specifically for data compression and coding. A larger attention has received the cipher technology and just lately by the introduction of the symbolic dynamics approach [START_REF] Schweizer | Symbolic Dynamics for Processing Chaotic Signals-I: Communication and Coding[END_REF], [START_REF] Schweizer | Symbolic Dynamics for Processing Chaotic Signals-II: Communication and Coding[END_REF] and the channel coded system considered, have really open the use of Piece Wise Linear Markov maps (PWLM) to the coding methods area.

This paper objective is to provide a method of a loss less compression algorithm and to demonstrate that the proposed algorithm achieves the optimal entropy criteria as level of compression.

The paper is organized as follows. Section II describes briefly the symbolic dynamics method and considers its application to a Bernoulli shift map. In section III the introduction of the probabilistic Bernoulli generator is considered, and finally in section IV the compression algorithm is presented and its optimality is proven.

II. THE SYMBOLIC DYNAMICS

Symbolic dynamics is concerned with the analysis of symbolic sequences that have been generated by a dynamical system. The adaptation of the definitions of symbolic dynamics to chaotic systems allows to partition the infinite number of finite length trajectories into a finite number of trajectory-sets.

This definition given in [START_REF] Schweizer | Symbolic Dynamics for Processing Chaotic Signals-I: Communication and Coding[END_REF] describes almost completely the general symbolic dynamics approach applied to chaotic sequences. In this section we will not consider all the mathematical apparatus employed in demonstrating the general statistic properties but a particular application to the Bernoulli shift map given by:

f (x) = ( f {1} (x) = 2x, x ∈ I {1} = [0, 0.5] f {2} (x) = 2x -1, x ∈ I {2} = (0.5, 1] (1)
In the above equation we have considered the disjoint and connected intervals I {n} , with ∩ N n=1 I {n} =Ø, and ∪ N n=1 I {n} = I where in this case N = 2, n = 1...N and I = [0, 1]. To simply define the symbolic dynamics we consider a sequence s = {s i , i = 1...M|s i ∈ S} where S represents the alphabet of symbols available, S {n} ∈ S, n = 1...N so to associate to every symbol S {n} an interval I {n} . In this way we have defined a bijective function between the division ∪ N n=1 I {n} = I of the chaotic function state space to the alphabet ∪ N n=1 S {n} = S. We have used the terminology of alphabet to make easier the explanation of the compression algorithm in the section IV. Now for every interval I {n} the probability generator associates a linear function f {n} : I {n} → I which on the particular interval I {n} is bijective and in consequence invertible, so we can affirm that exists f -1{n} : I → I {n} for every n = 1...N . As an example, let us consider s M = S {2} . We can deterministically identify the interval where the initial condition for the step M has to be. We have no limits for the step M +1, so the interval from which are mapping initially is concerning all the possible values (I). The function who can map any of the values of the interval I to the interval I {2} is given by f -1{2} and it's bijectivity ensures us that for any value in I it will correspond just one value in I {2} and vice versa. So finally for the M th symbol we have :

I {M} = f -1{2} (I) = I {2} (2)
where we have designated with I {i} the interval where the initial condition must be so that the symbolic sequence generated to be {s i , s i+1 , ..., s M }. Considering now that s M-1 = S {1} we will have this time as interval of possible values I {M} , and consequently as interval of possible initial conditions I {M-1} = f -1{1} (I {M} ). In general we can express the relation (2) for any symbol s i which can take any value S n by the following recurrence form:

I {i} = f -1{n} (I {i+1} ) ( 3 )
The application of back propagation relation given in (3) to the Bernoulli map (1) is exemplified in figure 1 for the sequence S {2} , S {1} , S {2} .

Fig. 1. The interval of possible IC deducted through back propagation

A similar approach has been largely developed by Schweizer and Schimming [START_REF] Schweizer | Symbolic Dynamics for Processing Chaotic Signals-I: Communication and Coding[END_REF] for a bias-free Bernoulli map, with the direct application to the noise reduction for PWLM generated chaotic sequences.

III. PROBABILISTIC BERNOULLI GENERATOR

In this section we introduce the probabilistic Bernoulli generator as non-uniform piecewise linear 1 st order Markov map. The non-uniformity comes from the idea that the intervals I {n} can be unequal. Now we consider the alphabet of symbols S with the probability of apparition of each symbol given by P n , n = 1...N. We have both for I {n} and P n the relations:

∪ N n=1 I {n} = I = [0, 1] P N n=1 P n = 1 ∩ N n=1 I {n} = φ P n -probability of disjoint processes (symbols) (4)
Using the relations (4) we can very easily develope a relation, between the discrete probability distribution and the intervals that define the state space for the symbolic dynamics, of the following form:

I {1} = [0, P 1 ] I {n} = (P n-1 , P n-1 + P n ], n = 2...N (5)
and also give a characteristic function for the probabilistic Bernoulli map:

f p (x) =            f {1} p (x) = 1 P 1 x, x ∈ I {1} f {2} p (x) = 1 P 2 (x -P 1 ) , x ∈ I {2} ..... f {N } p (x) = 1 P N ³ x - P N-1 n=1 P n ´, x ∈ I {N} (6)
A representation of this characteristic function is given in figure 2.

( f -1{1} p (x) = P 1 x f -1{n} p (x) = P n x + P n-1 i=1 P i (7)
The expression given in (7) will be the heart of the compression algorithm presented in section IV, employing the symbolic dynamics method.

IV. CHAOTIC COMPRESSION ALGORITHM

To achieve a certain compression level one has to associate to the compression algorithm some of the statistical properties of the informational signal to be compressed, the goal being to reduce the redundancy and so leaving only the information content. Generally the information of a source code s i is defined bylog 2 (p(s i )), and consequently there is defined the average information content over the source alphabet as the entropy of that alphabet (8):

H = - N X n=1 P n log 2 (P n ) ( 8 )
The optimality of the code is considered in the sense of the minimum redundancy, which is defined by the difference between the average codeword length and the entropy of the associated alphabet, and so a code is considered asymptotically optimal if for a given probability distribution, the ratio between the average codeword length to the entropy approaches 1 as entropy goes to infinity.

The purpose of introducing the symbolic dynamics in section II was to prove that starting with piece-wise particular generator we can "code", using a bijective recursive family of functions, a sequence of symbols to a deterministically obtained interval of initial conditions. We will prove that using this method and the probabilistic Bernoulli generator introduced in section III we can attain the entropy limit for the code obtained.

The compression is achieved in finding the best initial condition, coded in binary format, that using the particular form of the generator defined in (6) can obtain with the symbolic dynamics approach give the initial sequence. Now if we want to find the performances of this method we have to refer to the minimum necessary binary sequence to code a initial condition in a particular interval. It is known that to represent a value in the interval I ⊂ [0, 1) of dimension d we needlog 2 d number of bits, so if we can develop a method for calculating the size of the particular interval that the initial condition has to respect, we can provide a compression performance of the algorithm.

We consider as in section II the sequence s = {s i , i = 1...M|s i ∈ S} where S is the alphabet of symbols available, S {n} ∈ S, n = 1...N, with the symbol apparition probability given by:

P n = 1 M card n s i |s i = S {n} o (9)
To this discrete probability distribution, we associate the probabilistic Bernoulli generator given in (6) with the property of partial function invertibility with the form given in eq. ( 7) and the disjunctivity of the definition intervals I {n} .

If we consider the back-propagation method we can provide the same recurrent form, for the interval I {i} of possible values of initial conditions, as it is given in equation (3); the symbol coded was s i = S {n} :

I {i} = f -1{n} p (I {i+1} ) ( 1 0 )
Using the relation above and the linear form of f -1{.} p we can provide also a recurrent expression for the size of the interval I {i} (we are also working with the hypothesis that the symbol coded was s i = S {n} , and 0 ≤ P n < 1):

size(I {i} ) = P n size(I {i+1} ) ( 1 1 )
and give the size of possible initial conditions for all the sequence s :

size(I {1} ) = Y M i=1 P (s i ∈ S) = Y N n=1 (P n ) P n M (12)
we have used the definition for the discrete probability distribution given in (9). Now we can obtain the number of necessary bits to code the initial condition:

-log 2 ³ size(I {1} ) ´= - X N n=1 P n M log 2 (P n ) = M • H (13)
With the relation (13) the demonstration of the optimality of the compression algorithm is given. A similar classical algorithm was introduced in the 60's by a number of researchers and received the name of arithmetic algorithm. It uses also a [0, 1] initial interval but with forward projecting technique and the direct implementations of this algorithm were prohibited by requiring high computational effort, but now there have been proposed some adjustments like a power of 2 alphabet distribution probability representation [START_REF]Arithmetic Coding and Decoding Methods and Related Systems[END_REF].

The decompression algorithm in our case is very simple, using the determined initial condition together with the probabilistic Bernoulli generator defined by {P n } we can generate a chaotic sequence and just make a state phase separation conformed to the intervals {I {n} }. The problem for this algorithm as for the arithmetic compression method is that it doesn't know when to stop so we have to add to the symbols space, a character of end-sequence, which will affect the optimality of compressing process as a whole but by using long sequences the algorithms achieves almost optimality.

To exemplify the functionality of the algorithm we have chosen two sequences with the same length and symbols apparition probability (table 14) to see how the recurrent determination of the interval for the initial condition changes and also the trajectory of the regenerated chaotic sequence (figures 3 a,b).

Considered

Entropy: H = 1.7108, normal coding for the IC (initial conditions):H•M = 30.7944

Coded IC: [0 0 0 1 1 0 0 1 0 0 0 1 0 1 1 0 1 1 0 1 0 1 1 0 0 0 1 1 0 0 1], length 31

Uncoded IC: x 0 = 0.09800471039489 Considered sequence #2: ''ABAACABCACCBACAAC#'', length: M = 18 Coded IC: [0 0 0 1 1 0 0 1 0 0 0 1 0 1 1 0 1 1 0 1 0 1 1 0 0 0 1 1 0 0 1], length 31

Uncoded IC: x 0 = 0.20770334312692 A qualitative explanation of figures 3 a and b, is that using the property of bijectivity of the functions f {n} p (x) ensures us that the regenerated trajectory determined by the initial condition x 0 will be always contained in the back propagation recurrent determination of the I {i} intervals, and limited by the inferior and the superior limits defining them. We can observe that for both seqences the trajectory generation stops when detecting end-sequence symbol. There can be considered other methods for the decrompression, like a limited length signal encoding permiting only a fix number of symbols to be processed each time, but each method can be particulary good for different informational signals with different discrete probability distributions, a more application directed study needing to be made. V. CONCLUSIONS An optimal entropy criteria loss-less compression algorithm was introduced as a possibly new application of chaotic signals to coding procedure. The use of chaotic dynamics and a special type of a chaotic generator, called probabilistic Bernoulli generator were considered to map a sequence into an initial condition associated with the generator. Only the static compression algorithm is considered in this article, an adaptive form and/or enriched Markov chain order can be envisaged to improve performances. The use of the PWLM map permits a fast reconstruction of the initial sequence making it very interesting for real time powerful implementations.

Fig. 2 .

 2 Fig. 2. Bernoulli probabilistic characteristic function Now the function introduced in (6) has the same properties as the characteristic function of the general Bernoulli map presented in section II, and especially the partial defined function f {n} p : I {n} → I is bijective on the particular interval I {n} and in consequence invertible, so we can affirm that exists f -1{n} p

Fig. 3 .

 3 Fig. 3. Trajectory limits to determine the recurrent interval of IC & regenerated trajectory for: a) sequence #1, b) sequence #2


Emeric Bouin
email: bouin@ceremade.dauphine.fr

Franca Hoffmann

Clément Mouhot
email: c.mouhot@dpmms.cam.ac.uk

Exponential decay to equilibrium for a fibre lay-down process on a moving conveyor belt

Keywords: Hypocoercivity, rate of convergence, fibre lay-down, existence and uniqueness of a stationary state, perturbation, moving belt, AMS Class. No: 35B20, 35B40, 35B45, 35Q84

We show existence and uniqueness of a stationary state for a kinetic Fokker-Planck equation modelling the fibre lay-down process in the production of non-woven textiles. Following a micro-macro decomposition, we use hypocoercivity techniques to show exponential convergence to equilibrium with an explicit rate assuming the conveyor belt moves slow enough. This work is an extension of (Dolbeault et al., 2013), where the authors consider the case of a stationary conveyor belt. Adding the movement of the belt, the global Gibbs state is not known explicitly. We thus derive a more general hypocoercivity estimate from which existence, uniqueness and exponential convergence can be derived. To treat the same class of potentials as in (Dolbeault et al., 2013), we make use of an additional weight function following the Lyapunov functional approach in (Kolb et al., 2013).

Introduction

The mathematical analysis of the fibre lay-down process in the production of non-woven textiles has seen a lot of interest in recent years [14,15,6,10,11,4,12]. Non-woven materials are produced in melt-spinning operations: hundreds of individual endless fibres are obtained by continuous extrusion through nozzles of a melted polymer. The nozzles are densely and equidistantly placed in a row at a spinning beam. The visco-elastic, slender and in-extensible fibres lay down on a moving conveyor belt to form a web, where they solidify due to cooling air streams. Before touching the conveyor belt, the fibres become entangled and form loops due to highly turbulent air flow. In [14] a general mathematical model for the fibre dynamics is presented which enables the full simulation of the process. Due to the huge amount of physical details, these simulations of the fibre spinning and lay-down usually require an extremely large computational effort and high memory storage, see 1 [15]. Thus, a simplified two-dimensional stochastic model for the fibre lay-down process, together with its kinetic limit, is introduced in [6]. Generalisations of the two-dimensional stochastic model [6] to three dimensions have been developed by Klar et al. in [10] and to any dimension d ≥ 2 by Grothaus et al. in [7].

We now describe the model we are interested in, which comes from [6]. We track the position x(t) ∈ R 2 and the angle α(t) ∈ S 1 of the fibre at the lay-down point where it touches the conveyor belt. Interactions of neighbouring fibres are neglected. If x 0 (t) is the lay-down point in the coordinate system following the conveyor belt, then the tangent vector of the fibre is denoted by τ (α(t)) with τ (α) = (cos α, sin α). Since the extrusion of fibres happens at a constant speed, and the fibres are in-extensible, the lay-down process can be assumed to happen at constant normalised speed x ′ 0 (t) = 1. If the conveyor belt moves with constant speed κ in direction e 1 = (1, 0), then dx dt = τ (α) + κe 1 .

Note that the speed of the conveyor belt cannot exceed the lay-down speed: 0 ≤ κ ≤ 1. The fibre dynamics in the deposition region close to the conveyor belt are dominated by the turbulent air flow. Applying this concept, the dynamics of the angle α(t) can be described by a deterministic force moving the lay-down point towards the equilibrium x = 0 and by a Brownian motion modelling the effect of the turbulent air flow. We obtain the following stochastic differential equation for the random variable

X t = (x t , α t ) on R 2 × S 1 ,    dx t = (τ (α t ) + κe 1 ) dt, dα t = -τ ⊥ (α t ) • ∇ x V (x t ) dt + A dW t , (1)
where W t denotes a one-dimensional Wiener process, A > 0 measures its strength relative to the deterministic forcing, τ ⊥ = (-sin α, cos α), and V : R 2 -→ R is an external potential carrying information on the coiling properties of the fibre. More precisely, since a curved fibre tends back to its starting point, the change of the angle α is assumed to be proportional to τ ⊥ (α)•∇ x V (x). It has been shown in [12] that under suitable assumptions on the external potential V , the fibre lay down process (1) has a unique invariant distribution and is even geometrically ergodic (see Remark 1.2).

The stochastic approach yields exponential convergence in total variation norm, however without explicit rate. We will show here that a stronger result can be obtained with a functional analysis approach. Our argument uses crucially the construction of an additional weight functional for the fibre lay-down process in the case of unbounded potential gradients inspired by [12,Proposition 3.7].

The probability density function f (t, x, α) corresponding to the stochastic process (1) is governed by the Fokker-Planck equation

∂ t f + (τ + κe 1 ) • ∇ x f -∂ α τ ⊥ • ∇ x V f = D∂ αα f (2)
with diffusivity D = A 2 /2. We state below assumptions on the external potential V that will be used regularly throughout the paper:

(H1) Regularity and symmetry: V ∈ C 2 (R 2 ) and V is spherically symmetric outside some ball B(0, R V ).

(H2) Normalisation: R 2 e -V (x) dx = 1.

(H3) Spectral gap condition (Poincaré inequality): there exists a positive constant Λ such that for any u ∈ H 1 (e -V dx) with R 2 ue -V dx = 0,

R 2 |∇ x u| 2 e -V dx ≥ Λ R 2
u 2 e -V dx.

(H4) Pointwise regularity condition on the potential: there exists c 1 > 0 such that for any

x ∈ R 2 , the Hessian ∇ 2 x V of V (x) satisfies |∇ 2 x V (x)| ≤ c 1 (1 + |∇ x V (x)|).
(H5) Behaviour at infinity:

lim |x|→∞ |∇ x V (x)| V (x) = 0, lim |x|→∞ |∇ 2 x V (x)| |∇ x V (x)| = 0 .
Remark 1.1. Assumptions (H2-3-4) are as stated in [4]. Assumption (H1) assumes regularity of the potential that is stronger and included in that discussed in [4] since (H1) implies V ∈ W 2,∞ loc (R 2 ). Assumption (H5) is only necessary if the potential gradient |∇ x V | is unbounded. Both bounded and unbounded potential gradients may appear depending on the physical context, and we will treat these two cases separately where necessary. A typical example for an external potential satisfying assumptions (H1-2-3-4-5) is given by

V (x) = K 1 + |x| 2 s/2
(3)

for some constants K > 0 and s ≥ 1 [5,12]. The potential (3) satisfies (H3) since

lim inf |x|→∞ |∇ x V | 2 -2∆ x V > 0 ,
see for instance [17, A.19. Some criteria for Poincaré inequalities, page 135]. The other assumptions are trivially satisfied as can be checked by direct inspection. In this family of potentials, the gradient ∇ x V is bounded for s = 1 and unbounded for s > 1.

Remark 1.2. The proof of ergodicity in [12] assumes that the potential satisfies

lim |x|→∞ |∇ x V (x)| V (x) = 0 , lim |x|→∞ |∇ 2 x V (x)| |∇ x V (x)| = 0 , lim |x|→∞ |∇ x V (x)| = ∞ . ( 4
)
Under these assumptions, there exists an invariant distribution ν to the fibre lay-down process (1), and some constants C(x 0 ) > 0, λ > 0 such that

P x 0 ,α 0 (X t ∈ •) -ν T V ≤ C(x 0 )e -λt ,
where • T V denotes the total variation norm. The stochastic Lyapunov technique applied in [12] however does not give any information on how the constant C(x 0 ) depends on the initial position x 0 , or how the rate of convergence λ depends on the conveyor belt speed κ, the potential V and the noise strength A. This can be achieved using hypocoercivity techniques, proving convergence in a weighted L 2 -norm, which is slightly stronger than the convergence in total variation norm shown in [12]. Conceptually, the conditions (4) ensure that the potential V is driving the process back inside a compact set where the noise can be controlled. Our framework (H1-2-3-4-5) is more general than conditions (4) in some aspects (including bounded potential gradient) and more restrictive in others (assuming a Poincaré inequality). The proof in [12] relies on the strong Feller property which can be translated in some cases into a spectral gap; it also uses hypoellipticity to deduce the existence of a transition density, and concludes via an explicit Lyapunov function argument. With our framework (H1-2-3-4-5), and adapting the Lyapunov function argument presented in [12] to control the effect of κ∂ x 1 , we derive an explicit rate of convergence in terms of κ, D and V .

To set up a functional framework, rewrite (2) as

∂ t f = L κ f = (Q -T) f + P κ f , (5)
where the collision operator Q := D∂ αα acts as a multiplicator in the space variable x, P κ is the perturbation introduced by the moving belt with respect to [4]:

P κ f := -κe 1 • ∇ x f ,
and the transport operator T is given by

Tf := τ • ∇ x f -∂ α τ ⊥ • ∇ x V f .
We consider solutions to (5) in the space

L 2 (dµ κ ) := L 2 (R 2 × S 1 , dµ κ ) with measure dµ κ (x, α) = e V (x) + ζκg(x, α) dx dα 2π .
We denote by •, • κ the corresponding scalar product and by • κ the associated norm. Here, ζ > 0 is a free parameter to be chosen later. The construction of the weight g depends on the boundedness of ∇ x V . When it is bounded, no additional weight is needed to control the perturbation, and so we simply set g ≡ 0 in that case. When the gradient is unbounded, the weight is constructed thanks to the following proposition:

Proposition 1.3. Assume that V satisfies (H1) and (H5) and that

lim |x|→∞ |∇ x V | = +∞.
If κ < 1/3 holds true, then there exists a function g(x, α), a constant c = c(κ, D) > 0 and a finite radius

R = R(k, D, V ) > 0 such that ∀ |x| > R, ∀ α ∈ S 1 , L κ (g)(x, α) ≤ -c |∇ x V (x)|g(x, α) , (6)
where L κ is defined by

L κ (h) := D∂ αα h + (τ + κe 1 ) • ∇ x h -τ ⊥ • ∇ x V ∂ α h -(τ • ∇ x V ) h . (7)
The weight g is of the form

g(x, α) := exp βV (x) + |∇ x V (x)|Γ τ (α) • ∇ x V (x) |∇ x V (x)| ,
where the parameter β > 1 and the function Γ ∈ C 1 ([-1, 1]), Γ > 0 are determined along the proof and only depend on κ.

We show in Section 3 the existence of such a weight function g under appropriate conditions following ideas from [12].

We denote C := C ∞ c R2 × S 1 , and define the orthogonal projection Π on the set of local equilibria Ker Q

Πf := S 1 f dα 2π ,
and the mass M f of a given distribution f ∈ L 2 (dµ κ ),

M f := R 2 ×S 1 f dxdα 2π .
Integrating (2) over R 2 × S 1 shows that the mass of solutions of ( 2) is conserved over time, and standard maximum principle arguments show that it remains non-negative for non-negative initial data. The collision operator Q is symmetric and satisfies

∀ f ∈ C, Qf, f 0 = -D ∂ α f 2 0 ≤ 0 , i.e. Q is dissipative in L 2 (dµ 0 ). Further, we have TΠf = e -V τ • ∇ x u f for f ∈ C, with u f := e V Πf ,
which implies ΠTΠ = 0 on C. Since the transport operator T is skew-symmetric with respect to

• , • 0 , L κ f, f 0 = Qf, f 0 + P κ f, f 0 for any f in C.
In the case κ = 0, if the entropy dissipation -Qf, f 0 was coercive with respect to the norm • 0 , exponential decay to zero would follow as t → ∞. However, such a coercivity property cannot hold since Q vanishes on the set of local equilibria. Instead, Dolbeault et al.

[5] applied a strategy called hypocoercivity (as theorised in [17]) and developed by several groups in the 2000s, see for instance [9,8,13,2,3]. The full hypocoercivity analysis of the long time behaviour of solutions to this kinetic model in the case of a stationary conveyor belt, κ = 0, is completed in [4]. For technical applications in the production process of non-wovens, one is interested in a model including the movement of the conveyor belt, and our aim is to extend the results in [4] to small κ > 0.

We follow the approach of hypocoercivity for linear kinetic equations conserving mass developed in [5], with several new difficulties. Considering the case κ = 0, Q and T are closed operators on L 2 (dµ 0 ) such that Q -T generates the C 0 -semigroup e (Q-T)t on L 2 (dµ 0 ). When κ > 0, we use the additional weight function g > 0 to control the perturbative term P κ in the case of unbounded potential gradients; and show the existence of a C 0 -semigroup for L κ = Q -T + P κ (see Section 4.1). Unless otherwise specified, all computations are performed on the operator core C, and can be extended to L 2 (dµ κ ) by density arguments.

When κ = 0, the hypocoercivity result in [5,4] is based on: microscopic coercivity, which assumes that the restriction of Q to (Ker Q) ⊥ is coercive, and macroscopic coercivity, which is a spectral gap-like inequality for the operator obtained when taking a parabolic drift-diffusion limit, in other words, the restriction of T to Ker Q is coercive. The two properties are satisfied in the case of a stationary conveyor belt:

• The operator Q is symmetric and the Poincaré inequality on S 1 ,

1 2π S 1 |∂ α f | 2 dα ≥ 1 2π S 1 f - 1 2π S 1 f dα implies that -Qf, f 0 ≥ D (1 -Π)f 2 0 .
• The operator T is skew-symmetric and for any h ∈ L 2 (dµ 0 ) such that u h = e V Πh ∈ H 1 (e -V dx) and R 2 ×S 1 h dµ 0 = 0, (H3) implies

TΠh 2 0 = 1 4π R 2 ×S 1 e -V |∇ x u h | 2 dx dα ≥ Λ 4π R 2 ×S 1 e -V u 2 h dx dα = Λ 2 Πh 2 0 .
In the case κ = 0, the unique global normalised equilibrium distribution F 0 = e -V lies in the intersection of the null spaces of T and Q. When κ > 0, F 0 is not in the kernel of P κ and we are not able to find the global Gibbs state of (5) explicitly. However, the hypocoercivity theory is based on a priori estimates [5] that are, as we shall prove, to some extent stable under perturbation. Our main result reads:

Theorem 1.4. Let f in ∈ L 2 (dµ κ ) and let (H1-2-3-4-5) hold. For 0 < κ < 1 small enough (with a quantitative estimate) and ζ > 0 large enough (with a quantitative estimate), there exists a unique non-negative stationary state F κ ∈ L 2 (dµ κ ) with unit mass M Fκ = 1. In addition, for any solution f of (2) in L 2 (dµ κ ) with mass M f and subject to the initial condition f (t = 0) = f in , we have

f (t, •) -M f F κ κ ≤ C f in -M f F κ κ e -λκt , (8)
where the rate of convergence λ κ > 0 depends only on κ, D and V , and the constant C > 0 depends only on D and V .

In the case of a stationary conveyor belt κ = 0 considered in [4], the stationary state is characterised by the eigenpair (Λ 0 , F 0 ) with Λ 0 = 0, F 0 = e -V , and so Ker L 0 = F 0 . This means that there is an isolated eigenvalue Λ 0 = 0 and a spectral gap of size at least [-λ 0 , 0] with the rest of the spectrum Σ(L 0 ) to the left of -λ 0 in the complex plane. Adding the movement of the conveyor belt, Theorem 1.4 shows that Ker L κ = F κ and the exponential decay to equilibrium with rate λ κ corresponds to a spectral gap of size at least [-λ κ , 0]. Further, it allows to recover an explicit expression for the rate of convergence λ 0 for κ = 0 (see Step 5 in Section 2.1). In general, we are not able to compute the stationary state F κ for κ > 0 explicitly, but F κ converges to F 0 = e -V weakly as κ → 0 (see the discussion in Section 5). Let us finally emphasize that a specific contribution of our paper is to introduce two (and not one as in [5,4]) modifications of the entropy: 1) we first modify the space itself with the coercivity weight g, then 2) we change the norm with an auxiliary operator following the hypocoercivity approach.

The rest of the paper deals with the case κ > 0 and is organised as follows. In Section 2, we prove the main hypocoercivity estimate. This allows us to establish the existence of solutions to (2) using semigroup theory and to deduce the existence and uniqueness of a steady state in Section 4 by a contraction argument. In Section 3, we give a detailed definition of the weight function g that is needed for the hypocoercivity estimate in Section 2.

Hypocoercivity estimate

Following [5] we introduce the auxiliary operator A := (1 + (TΠ) * (TΠ)) -1 (TΠ) * , and a modified entropy, i.e. a hypocoercivity functional G on L 2 (dµ κ ):

G[f ] := 1 2 f 2 κ + ε 1 Af, f 0 , f ∈ L 2 (dµ κ )
for some suitably chosen ε 1 ∈ (0, 1) to be determined later. It follows from

[5] that | Af, f 0 | ≤ f 2 0 . Also, f 2 0 ≤ f 2 κ by construction of µ κ , and hence G[•] is norm-equivalent to • 2 κ : ∀ f ∈ L 2 (dµ κ ), 1 -ε 1 2 f 2 κ ≤ G[f ] ≤ 1 + ε 1 2 f 2 κ , (9)
In this section, we prove the following hypocoercivity estimate:

Proposition 2.1. Assume that hypothesis (H1-2-3-4-5) hold and that 0 < κ < 1 is small enough (with a quantitative estimate). Let f in ∈ L 2 (dµ κ ) and f = f (t, x, α) be a solution of (2) in L 2 (dµ κ ) subject to the initial condition f (t = 0) = f in . Then f satisfies the following Grönwall type estimate:

d dt G[f (t, •)] ≤ -γ 1 G[f (t, •)] + γ 2 M 2 f , (10)
where γ 1 > 0, γ 2 > 0 are explicit constants only depending on κ, D and V .

Note that the estimate ( 10) is stronger than what is required for the uniqueness of a global Gibbs state, and represents an extension of the estimate given in [4]. When applied to the difference of two solutions with the same mass, (10) gives an estimate on the exponential decay rate towards equilibrium.

Proof of Proposition 2.1

Differentiate in time G[f ] to get d dt G[f ] = D 0 [f ] + D 1 [f ] + D 2 [f ] + D 3 [f ] ,
where the entropy dissipation functionals D 0 , D 1 , D 2 and D 3 are given by

D 0 [f ] := Qf, f 0 -ε 1 ATΠf, Πf 0 -ε 1 AT(1 -Π)f, Πf 0 + ε 1 TAf, (1 -Π)f 0 + ε 1 AQf, Πf 0 , D 1 [f ] := ε 1 AP κ f, Πf 0 + ε 1 P * κ Af, Πf 0 , D 2 [f ] := P κ f, f 0 , D 3 [f ] := κζ R 2 ×S 1 L κ (f )f g dx dα 2π .
Note that the term LAf, f 0 vanishes since it has been shown in [5] that A = ΠA and hence Af ∈ Ker Q. Further, Tf, f = 0 since T is skew-symmetric. We estimate the entropy dissipation of the case κ = 0 as in [4]:

# Step 1: Estimation of D 0 [f ].
We will show the boundedness of D 0 , which is in fact the dissipation functional for a stationary conveyor belt. We thus recall without proof in the following lemma some results from [4]. Lemma 2.2 (Dolbeault et al. [4]). The following estimates hold:

Qf, f 0 ≤ -(1 -Π)f 2 0 , AT(1 -Π)f 0 ≤ C V (1 -Π)f 0 , AQf 0 ≤ D 2 (1 -Π)f 0 , TAf 0 ≤ (1 -Π)f 0 .
In order to control the contribution ATΠf, Πf 0 in D 0 , we note that ATΠ = (1 + (TΠ) * TΠ) -1 (TΠ) * TΠ shares its spectral decomposition with (TΠ) * TΠ, and by macroscopic coercivity

(TΠ) * TΠf, f 0 = TΠf 2 0 = TΠ(f -M f e -V ) 2 0 ≥ Λ 2 Π(f -M f e -V ) 2 0 .
Hence,

ATΠf, f 0 ≥ Λ/2 1 + Λ/2 Π(f -M f e -V ) 2 0 .
Now, recalling Lemma 2.2 and using

Π(f -M f e -V ) 2 0 = Πf 2 0 -M 2 f , we estimate D 0 [f ] ≤(ε 1 -D) (1 -Π)f 2 0 + ε 1 λ 2 (1 -Π)f 0 Πf 0 -ε 1 γ 2 Πf 2 0 -M 2 f , with λ 2 := C V + D/2 > 0 and γ 2 := Λ/2 1+Λ/2 > 0. # Step 2: Estimation of D 1 [f ].
We now turn to the entropy dissipation functional D 1 , which we will estimate using elliptic regularity. Instead of bounding AP κ , we apply an elliptic regularity strategy to its adjoint, as for AT(1 -Π) in [4]. Let f ∈ L 2 (dµ 0 ) and define h := (1 + (TΠ) * TΠ) -1 f so that u h = e V Πh satisfies

Πf = e -V u h + ΠT * T e -V u h = e -V u h - 1 2 ∇ x • e -V ∇ x u h .
We have used here the fact that in the space L 2 (dµ 0 ):

T = τ • ∇ x -∂ α τ ⊥ • ∇ x V , T * = -τ • ∇ x + τ ⊥ • ∇ x V ∂ α -(τ • ∇ x V ) . Then A * f = TΠh = e -V τ • ∇ x u h ,
and since the adjoint for •, • 0 of the perturbation operator P κ is given by

P * κ = -P κ -P κ V , it follows that (AP κ ) * f 2 0 = κ τ • ∇ x (e 1 • ∇ x u h )e -V 2 0 = κ 2 2 R 2 ×S 1 e -V |τ • ∇ x (e 1 • ∇ x u h ) | 2 dµ 0 = κ 2 2 R 2 e -V |∇ x (e 1 • ∇ x u h ) | 2 dx ≤ κ 2 2 ∇ 2 x u h 2 L 2 (e -V dx) ≤ κ 2 2 C 2 V Πf 2 0 ,
where in the last inequality we have used an elliptic regularity estimate. This estimate turns out to be a particular case of [4, Proposition 5 and Sections 2-3], where the positive constant C V is the same as in Lemma 2.2 reproduced from [4]. This concludes the boundedness of AP κ ,

AP κ f 0 ≤ κ C V √ 2 Πf 0 ≤ κ C V √ 2 f 0 . (11)
Using a similar approach for the operator P * κ A, we rewrite its adjoint as

A * P κ f = TΠ h ,
where we define h := (1 + (TΠ) * TΠ) -1 P κ f for a given f ∈ L 2 (dµ 0 ), or equivalently

e -V u h - 1 2 ∇ x • e -V ∇ x u h = ΠP κ f = P κ Πf .
Multiplying by u h and integrating over R 2 , we have

u h 2 L 2 (e -V dx) + 1 2 ∇ x u h 2 L 2 (e -V dx) = -κ R 2 e 1 • ∇ x (Πf ) u h dx = κ R 2 (Πf ) e 1 • ∇ x u h dx ≤ κ R 2 ∇ x u he -V /2 Πf e V /2 dx ≤ κ ∇ x u h L 2 (e -V dx) Πf 0 ≤ 1 4 ∇ x u h 2 L 2 (e -V dx) + κ 2 Πf 2 0 .
This inequality is a H 1 (e -V dx) → H -1 (e -V dx) elliptic regularity result. Hence,

A * P κ f 2 0 = TΠh 2 0 = 1 2 ∇ x u h 2 L 2 (e -V dx) ≤ 2κ 2 Πf 2 0 ,
and so we conclude

P * κ Af 0 ≤ √ 2κ (1 -Π)f 0 ≤ √ 2κ f 0 . (12)
Combining ( 11) and ( 12), the entropy dissipation functional D 1 is bounded by

D 1 [f ] ≤ κε 1 C V √ 2 + √ 2 f 2 0 = 2κλ 1 f 2 0 ,
where we defined λ

1 := 1 2 C V √ 2 + √ 2 . # Step 3: Estimation of D 2 [f ].
Using integration by parts, we have

P κ f, f 0 = κ 2 R 2 ×S 1 (e 1 • ∇ x V ) f 2 e V dx dα 2π .
The estimation of this term goes differently depending on the boundedness of

∇ x V . If ∇ x V is bounded, we write D 2 [f ] ≤ | P κ f, f 0 | ≤ κ 2 ∇ x V ∞ f 2 0 = κ 2 ∇ x V ∞ f 2 κ ,
where we have used

f κ = f 0 , since g ≡ 0. Assume now that |∇ x V | → ∞ as |x| → ∞.
Thanks to the choice of g, we have the estimate

D 2 [f ] ≤ | P κ f, f 0 | ≤ κ 2 R 2 ×S 1 |∇ x V |f 2 e V dx dα 2π ≤ κ 2 C 3 R 2 ×S 1 f 2 g dx dα 2π , (13)
with

C 3 := sup x∈R 2 |∇ x V |e V g -1 ,
which is finite by (H5).

# Step 4: Estimation of D 3 [f ].
We start by recalling that this estimate is only relevant when ∇ x V is unbounded. Indeed, in the opposite case, D 3 [f ] = 0 since g ≡ 0 by definition. By the identity

R 2 ×S 1 L κ (f )f g dx dα = 1 2 R 2 ×S 1 L κ (g)f 2 dx dα -D R 2 ×S 1 |∂ α f | 2 g dx dα
with L κ as defined in (7), we have

D 3 [f ] ≤ κζ 1 2 R 2 ×S 1 L κ (g)f 2 dx dα 2π . (14)
Proposition 1.3 allows us to control the g-weighted L 2 -norm outside some fixed ball. More precisely, take R > 0 in (6) large enough s.t.

|∇ x V | ≥ 1 for all |x| > R, then R 2 ×S 1 L κ (g)f 2 dx dα 2π ≤ S 1 |x|<R L κ (g)f 2 dx dα 2π -c S 1 |x|>R |∇ x V |f 2 g dx dα 2π ≤ S 1 |x|<R (L κ (g) + cg) e -V f 2 e V dx dα 2π -c R 2 ×S 1 f 2 g dx dα 2π ≤ C 4 (R) f 2 0 -c R 2 ×S 1 f 2 g dx dα 2π , (15)
where C 4 (R) := sup |x|≤R |L κ (g) + cg|e -V .

Remark 2.3. Observe here that one could take advantage of the growth of ∇ x V by playing with the cut-off parameter R and keeping track of min |x|≥R |∇ x V | in the negative term. It could lead to more optimal constants but we chose instead to vary the parameter ζ in front of the coercivity weight g in the measure µ κ for simplicity.

#

Step 5: Putting the four previous steps together.

Combine the previous steps into

D 0 [f ] + D 1 [f ] ≤(ε 1 -D) (1 -Π)f 2 0 + ε 1 λ 2 (1 -Π)f 0 Πf 0 -ε 1 γ 2 Πf 2 0 -M 2 f + 2κλ 1 f 2 0 = -(D -ε 1 -2κλ 1 ) (1 -Π)f 2 0 + ε 1 λ 2 (1 -Π)f 0 Πf 0 -(ε 1 γ 2 -2κλ 1 ) Πf 2 0 + ε 1 γ 2 M 2 f ≤ -D -ε 1 -2κλ 1 - ε 1 λ 2 b 2 (1 -Π)f 2 0 -ε 1 γ 2 -2κλ 1 - ε 1 λ 2 2b Πf 2 0 + ε 1 γ 2 M 2 f ≤ -2ξ(κ) f 2 0 + ε 1 γ 2 M 2 f ,
by Young's inequality with the choice b = λ 2 /γ 2 , and where we used the fact that (1

-Π)f 2 0 + Πf 2 0 = f 2 0 .
Here, ξ(κ) is explicit, and given by

ξ(κ) := 1 2 min D -ε 1 1 + λ 2 2 2γ 2 , ε 1 γ 2 2 -κλ 1 = Dγ 2 2 2 γ 2 2 + 2γ 2 + λ 2 2 -κλ 1 ,
since the minimum in the first term is realised when the two arguments are equal, fixing ε 1 = 2Dγ 2 / γ 2 2 + 2γ 2 + λ 2 2 . Note that this choice of ε 1 satisfies ε 1 < D and ε 1 < 1. Choosing κ small enough ensures ξ(κ) > 0. From this analysis we conclude

D 0 [f ] + D 1 [f ] ≤ -2ξ(κ) f 2 0 + ε 1 γ 2 M 2 f . (16)
Let us now add the control of

D 2 + D 3 . If ∇ x V is bounded, g ≡ 0 and D 3 = 0: d dt G[f ] = D 0 [f ] + D 1 [f ] + D 2 [f ] ≤ -(4ξ(κ) -κ ∇ x V ∞ ) 1 2 f 2 κ + ε 1 γ 2 M 2 f ≤ -γ 1 G[f ] + ε 1 γ 2 M 2 f
by the norm equivalence (9). Here, we defined

γ 1 := 4ξ(κ) -κ ∇ x V ∞ 1 + ε 1 > 0 .
When ∇ x V is unbounded, ( 13)-( 14)-( 15)-( 16) imply

d dt G[f ] =D 0 [f ] + D 1 [f ] + D 2 [f ] + D 3 [f ] ≤ -2ξ(κ) f 2 0 + ε 1 γ 2 M 2 f + κ 2 C 3 R 2 ×S 1 f 2 g dx dα 2π + κζ 2 C 4 (R) f 2 0 -c R 2 ×S 1 f 2 g dx dα 2π = - 1 2 (4ξ(κ) -κζC 4 (R)) f 2 0 - κζ 2 c - C 3 ζ R 2 ×S 1 f 2 g dx dα 2π + ε 1 γ 2 M 2 f ≤ - 1 2 min 4ξ(κ) -κζC 4 (R), c - C 3 ζ f 2 κ + ε 1 γ 2 M 2 f ≤ -γ 1 G[f ] + ε 1 γ 2 M 2 f
again by norm equivalence ( 9), and where we defined

γ 1 := 1 1 + ε 1 min 4ξ(κ) -κζC 4 (R), c - C 3 ζ > 0 .
This requires ζ > 0 to be large enough, and the upper bound for κ should be chosen accordingly:

ζ > C 3 c , 4ξ(κ) -κζC 4 (R) > 0 .
In order to maximise the rate of convergence to equilibrium given κ, D and V , one can optimise γ 1 over ζ whilst respecting the above constraints.

Remark 2.4. The condition γ 1 > 0 translates into an explicit upper bound on κ. More precisely, we require ξ(κ) > κu/4 where u := ∇ x V ∞ in the case of a bounded potential gradient, and u := ζC 4 (R) otherwise. This condition is satisfied for small enough κ:

0 ≤ κ < ε 1 γ 2 (4λ 1 + u) = 2Dγ 2 2 (4λ 1 + u)(γ 2 2 + 2γ 2 + λ 2 2 )
which also implies ξ(κ) > 0. Recall that Proposition 1.3 requires κ < 1/3 in the case of unbounded potential gradients. These conditions provide a range of κ for which Proposition 2.1 holds.

The coercivity weight g

In this section, we define the function g in such a way that it allows us to control the loss of weight in the perturbation operator P κ . When ∇ x V is bounded, we do not need any extra weight since then we may control the perturbation thanks to the stationary weight e V , and so we set g ≡ 0 in that case. When it is not, Proposition 1.3 provides a suitable weight function g by constructive methods.

Proof of Proposition 1.3

The proof is strongly inspired from [12], however our weight is different since we work in an L 2framework rather than in an L 1 one. Assuming ∇ x V is unbounded, we seek a weight g of the form

g(x, α) = exp βV (x) + |∇ x V (x)|Γ τ (α) • ∇ x V (x) |∇ x V (x)| ,
where the parameter β > 1 and the function Γ ∈ C 1 ([-1, 1]), Γ > 0 are to be determined. We define

Y (x, α) := τ (α) • ∇ x V (x) |∇ x V (x)| , Y ⊥ (x, α) := τ ⊥ (α) • ∇ x V (x) |∇ x V (x)| ,
and split the proof into four steps: 1) we rewrite statement (6) using the explicit expression of the weight g, 2) we simplify the obtained expression using assumption (H5), 3) we prove the equivalent statement obtained in Step 2 by defining a suitable choice of Γ(•) and β, and 4) we demonstrate that it is indeed possible to choose suitable parameters for the calculations in Step 3 to hold, fixing explicit expressions where possible.

# Step 1: Rewriting the weight estimate (6).

Applying the operator L κ defined in (7) to g, we can compute explicitly

L κ (g) g =D |∇ x V |∂ αα Γ(Y ) + |∇ x V | 2 |∂ α Γ(Y )| 2 + (τ (α) + κe 1 ) • (β∇ x V + ∇ x (|∇ x V |Γ(Y ))) -|∇ x V | 2 Y ⊥ ∂ α Γ(Y ) -|∇ x V |Y . Since ∂ α Γ = Y ⊥ Γ ′ (Y ) and ∂ αα Γ = ∂ α Y ⊥ Γ ′ (Y ) = -Y Γ ′ (Y ) + |Y ⊥ | 2 Γ ′′ (Y ) ,
we get

L κ (g) g =D |∇ x V | -Y Γ ′ (Y ) + |Y ⊥ | 2 Γ ′′ (Y ) + |∇ x V | 2 |Y ⊥ | 2 Γ ′ (Y ) 2 + (τ (α) + κe 1 ) • (β∇ x V + ∇ x (|∇ x V |Γ(Y ))) -|∇ x V | 2 |Y ⊥ | 2 Γ ′ (Y ) -|∇ x V |Y =(β -1 -DΓ ′ (Y ))|∇ x V |Y + κβe 1 • ∇ x V + (τ (α) + κe 1 ) • ∇ x (|∇ x V |Γ(Y )) + |Y ⊥ | 2 D|∇ x V |Γ ′′ (Y ) + |∇ x V | 2 D Γ ′ (Y ) 2 -Γ ′ (Y ) .
+ |Y ⊥ | 2 DΓ ′′ (Y ) + |∇ x V | D Γ ′ (Y ) 2 -Γ ′ (Y ) ≤ -c .

  In order to see which Γ to choose, let us divide by |∇ x V | and denote the diffusion and transport part bydiff(x, α) := (τ (α) + κe 1 ) • ∇ x (|∇ x V |Γ(Y )) |∇ x V | , tran(x) := e 1 • ∇ x V |∇ x V | .Now, we can rewrite the statement of Proposition 1.3: we seek a positive constant c > 0 and a radius R > 0 such that for any α ∈ S 1 and |x| > R, (β -1 -DΓ ′ (Y ))Y + κβtran(x) + diff(x, α)

dα,


P Chevret

M E Zakharia Lyon

ON SPATIAL UNCERTAINTY IN A SURFACE LONG BASELINE POSITIONING SYSTEM

This paper describes a general method for estimating the probability distribution of an acoustic source resulting from noisy measurements in a surface long baseline (SLBL) positioning system. Compared with Monte Carlo simulations this method calculates the mean and covariance at a reduced computational cost, due to a deterministic sampling strategy. Hence, it can easily be checked if a SLBL configuration is well suited to a particular positioning scenario.

INTRODUCTION

High precision positioning is of particular interest in various underwater application areas (marine salvage, offshore survey, mineral exploration, seafloor mapping, ROV tracking...). Numerous sophisticated systems are used for this purpose today but no global solution allowing performances comparable to GPS [1] is available yet. With traditional long baseline (LBL) acoustic positioning systems, a centimeter-level accuracy is sometimes attainable but the time required to deploy and calibrate the transponders array on the seabed become a severe drawback for some applications. Short baseline (SBL) and ultra short baseline (USBL) systems only require an acoustic unit composed of multiple hydrophones, mounted on a survey ship, but suffer from loss of accuracy in shallow water.

The idea to use GPS technology underwater has been a major advance over the past ten years. The first solution, due to Youngberg from the US-Air Force, is a direct transposition of the GPS technique to underwater world: satellites are replaced by surface buoys and the mobile, equipped with a receiver, processes location by measuring time of arrival (TOA) of acoustic pulses transmitted by the buoys. A reverse configuration has been proposed later by Thomas [1]: the principle is based on measuring on a set of buoys, the TOA of an acoustic pulse sent by the vehicle. At regular intervals, each buoy transmits to a survey ship its GPS position and the TOAs of the acoustic pulses. Assuming the sound velocity in seawater is known, the distances lead to the vehicle's location. Throughout this paper we will refer to this method as surface long baseline (SLBL) system. SLBL constitutes a very attractive positioning solution: it is easily deployable, and it can cover a wide operation area with meter-level accuracy.

Obviously, assumptions such as known sensor locations, good quality measurements or wellknown vehicle dynamic are not realistic. Due to complex and non-stationary characteristics of the acoustic environment and variable atmospheric conditions, it is necessary to reason on the basis of inaccurate information. As pointed out by Sutherland et al. [3] in a different context, the combination of approximate measures can lead to large errors in localization, even if robust tracking algorithms are applied. These localization errors can be greatly reduced by a good choice of the location of sensors [5].

In this paper, we describe a general method for estimating the spatial uncertainty in vehicle's localization. Our approach is based on a transformation proposed by Julier et al. [4] for calculating the statistics of a random variable (r.v.) which undergoes a non-linear transformation. Thanks to this method, the mean and covariance of the probability density function (PDF) of the vehicle's location are estimated in an efficient manner. The explicit form of the non-linear functions is not required and the first two moments of the PDF are precisely captured at a reduced computational cost in comparison with Monte Carlo simulations. This permits to decide in advance whether a particular SLBL configuration is able to provide good localization regardless of measurement errors.

The paper is organized as follows. In section 1 we briefly describe the principle of a SLBL. Section 2 examines various causes of error in such a system and how they affect errors in localization. Section 3 presents the Julier procedure to calculate the statistics of a r.v. obtained through a non linear function. Then, in section 4, we apply this method to model the uncertainty resulting from inexact information when using a SLBL. An illustrative example is proposed in section 5 with a towed sonar tracking scenario. Concluding remarks are given in section 6.

SURFACE LONG BASELINE POSITIONING

The SLBL is a recent concept [2] which uses the strong potential of the GPS while taking into account the principle of the LBL. The vehicle to be localized is equipped with a precision clock, synchronous to the GPS time. In order to avoid interrogation-response cycles usually utilized in such systems, the target transmits a signal at preprogrammed instants. Accordingly the sources of errors are decreased and the sampling rate is increased. The buoys can be anchor on the operation area or can be let drifting according to working conditions. Each buoy has to: detect the TOA of the acoustic signal transmitted by the target, date precisely this TOA using the clock reference of the GPS, update its own position using the GPS and transmit the collected information towards a central computer (by VHF in general).

The central computer determines the target coordinates at the transmission time using the whole collected data. Usually, this determination can be realized by hyperbolic or pseudo-circular algorithms. In this paper the pseudo-circular mode has been considered.

To improve its localization, the mobile can transmit its depth. This additional information improves the precision of the system and permits to reduce the number of unknown parameters to two (longitude and latitude). In order to keep simplicity and clearness figures are presented in 2D.

The main advantage of the SLBL, which can be assimilated to a reverse LBL, lies in its high precision without calibration even in long range operation. Moreover, the sensors are located far from any source of noise (except the surface of the sea) and then allow low level transmitted signals. These systems conceived by ACSA are manufactured and marketed in France by ORCA instrumentation. They allow a high rate of sampling in trajectory calculation.

THE EFFECTS OF INACCURATE INFORMATION

Assuming a two dimensional space, a straightforward analysis of the error-free SLBL localization problem yields the acoustic source location: knowledge of the sound speed and various TOA gives the distance of the source with respect to each buoy, and then a set of circles centered on the buoys intersect at a single point. Unfortunately, we must deal with noisy measurements and then the true location lies somewhere in an area of uncertainty, whose size and shape depend upon error amount and buoys positions. As noticed by Sutherland [3] for a robot navigating in an unstructured outdoor environment, it is crucial to precisely estimate this uncertainty.

The following errors will be encountered in a SLBL positioning system: the sound velocity used for calculation is an approximate value; as the signals of GPS satellites are corrupted by various errors (data noise, multipath errors, clock errors, atmospheric delays,…) the buoys locations are uncertain; the hydrophones are not precisely located with respect to the buoys because of sea state fluctuations; estimation of the TOA is subject to environmental noises (reflections, multipath,…).

In order to characterize the spatial uncertainty resulting from these inexact information, a least mean square (LMS) criteria will be employed in the following : given the set of circles coming from the measured distances between the acoustic source and the buoys, the searched position will correspond to the point which minimizes the sum of the squared perpendicular distances to the circles.

Then, the statistics of the source location can be explored by means of Monte Carlo simulations. This approach is illustrated in Fig. 1, which corresponds to 10000 trials, assuming Gaussian error models for the three buoys positions, the sound speed and the TOA. Supposing that the resulting probability distribution is Gaussian, the 95 % confidence ellipse is drawn to reveal the spatial uncertainty of the source location. Now we can ask if a better localization could be achieved with more than three buoys. Numerous simulations have been carried out to answer this question : we can note that the spatial uncertainty does not necessarily decrease by adding buoys, and in most cases a good three buoysconfiguration would perform better than a configuration having more buoys. This characteristic is clearly revealed by the second simulation (Fig. 2), performed with a different configuration of four buoys (same error models): We can observe that the first configuration (3 buoys) is less sensitive to measurement errors than the last one.

The purpose of this paper is to propose an efficient method to estimate the spatial uncertainty in order to find a configuration well suited to the desired positioning scenario.

In the following, the acoustic source location will be considered as a r.v. resulting from the non linear transformation of the r.v. associated to uncertain environmental parameters and noisy measurements (sound celerity, buoys location, TOA,…). The mean and covariance of the source probability density function (pdf) will be calculated through a general method recently developed by Julier et al. [4] for non linear estimation problems. This method avoids the large computational cost usually encountered with Monte Carlo simulation technique, and is more accurate than the standard linearisation process using a truncation of the Taylor expansion of the non linear function.

APPROXIMATING NONLINEAR TRANSFORMATIONS OF PROBABILITY DISTRIBUTION THROUGH JULIER'S METHOD

The general problem of approximating nonlinear transformations of probability distributions can be stated as follows : given a r.v. . To solve this problem, Julier had an intuition that with a fixed number of parameters it should be easier to approximate a Gaussian distribution than it is to approximate an arbitrary non linear transformation.

The procedure can be summarized as follows : 2. Each point is then transformed as

i i χ f Y Y Y Y ; 3.
The mean y is given by the weighted average of the transformed points :

i n i i W Y Y Y Y 2 0 y ; 4.
The predicted covariance yy P is computed as :

T i i n i i yy W y y P Y Y Y Y Y Y Y Y 2 0 .
In comparison with Monte Carlo techniques, the fundamental difference is that the samples are not drawn at random but according to a deterministic algorithm. A very small number of samples is then required to approximate the probability distribution. This method has many others advantages : owing to the parameter the scaling of the fourth and higher order moments can be influenced, the function f, which may be implicit, is not approximated through a truncation of its series expansion, and the implementation is very easy since no evaluation of Jacobians is needed.

ESTIMATION OF SPATIAL UNCERTAINTY IN A SLBL SYSTEM

The Julier's method can be easily applied to estimate the spatial uncertainty of localization using a SLBL. The r.v. x is composed of all the noisy data required to estimate the position of the target. For an SLBL system with N buoys the components of r.v. x are:

the Cartesian coordinates of the N buoys ( i x = longitude, i y = latitude), N 1,..., i the propagation times i t , N 1,..., i the sound velocity c . The r.v.

T T T y x , y

denotes the position of the target where T x and T y represent the longitude and the latitude respectively.

. f is the nonlinear transformation that permits to obtain the position of the target :

2 1 2 2 N i T i T i i y y y x x t c min Arg . f
The method is then applied to estimate the mean and the covariance of y . This result is depicted in a 95% confidence ellipse. To have a measure of the spatial uncertainty, we use the Frobenius norm of the covariance yy P .

ILLUSTRATIVE EXAMPLE

In this section we consider the experimental setup depicted in Fig. 5. An underwater mobile moves on a straight-line trajectory. The position of the target is to be tracked by a SLBL with three buoys (B 1 , B 2 , and B 3 ). The goal of this experiment is to estimate the uncertainty on the current position during the mobile movement by using the second order statistics.

In order to show the efficiency of this method, sinusoidal model is used for the position errors of the buoys, as a model of the sea swell:

T i i i i T T ) sin( ), sin( , , Bi Bi i i y x y x ,
where Bi Bi y x , is the true position of the buoy B i and i and i are independent and uniformly distributed. Although this model is very far from a Gaussian, as shown by Fig. 3, it is seen on Fig. 4 that the method performs well. For a three buoys-SLBL all the uncertain parameters have been represented by the r.v.

, , , Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y },
corresponding to the errors on second TOA 2 t and the celerity. Hence, a great advantage of this method is the capacity to identify the influence of each parameter separately.

CONCLUSIONS

The effects of inexact information have been studied in a SLBL acoustic positioning system. It has been emphasized that the spatial uncertainty is very sensitive to the layout of the buoys network and that this uncertainty does not necessarily reduces by adding buoys. The nonlinear transformations involved in the localization process have been handled owing to the Julier's algorithm to approximate the resulting probability distribution. This method avoids the large computational cost usually encountered with Monte Carlo simulations and provides very accurate mean and covariance. Also, the 'influence' of each noisy parameter upon the spatial uncertainty can be clearly identified. These results should be useful to check in advance the performance of a SLBL configuration, given a particular positioning scenario.

REFERENCE

[1] J. M. Coudeville, H. Thomas, "A primer: Using GPS underwater", Sea Technology, 1998. IEEE Trans. Robotics and Automation, vol. 10, no. 6, pp. 740-754, 1994. [4] J. S. Julier, J. K. Uhlmann, H. F. Durrant-Whyte, «A new approach for filtering nonlinear systems», In Proc. of the American Control Conf., pp. 1628-1632, 1995. [5] L. Collin, "Pré-étude d'un système de localisation de cibles en milieu sous-marin", Rapport de stage, DEA d'électronique, Univ. de Bretagne Occidentale, Brest, France, 1999.

Fig. 1 :

 1 Fig. 1: MonteCarlo simulation / 3 buoys Fig. 2: Monte Carlo simulation / 4 buoys

  x and covariance xx P , we would like to predict the mean y and covariance yy P of a r.v.

Fig. 6

 6 shows the variations of the estimated value of the norm F yy P during the linear displacement of the mobile. It is clearly seen that the localization is more accurate near the middle of the SLBL [5].

Fig. 3 :

 3 Fig.3: Sinusoidal and Gaussian distributions Fig 4: 95% Confidence ellipses. Fig. 5: Simulation mobile trajectory Fig. 6: Measure of localization uncertainty.

  [2] C. Viellard, S. Husy, "Surface long base line in deep offshore applications", In OCEANOLOGY international Conference, Brighton, pp. 489-496, 2000. [3] K. T. Sutherland, W. B. Thompson, «Localizing in unstructured environments: Dealing with the errors»,