Rigid Transformation Estimation in Point Registration Procedures: Study of the Use of Mahalanobis Distances

Keywords: point registration, rigid body transformation, Mahalanobis distance, ICP, ECMPR, polynomial systems

This document explores the problem of the estimation of the rigid body transformation that minimizes the sum of squares between matched points of dimension n where Mahalanobis distances are used. This problem is stated as a constrained minimization problem and the Lagrange multipliers framework is used. The solution is given as a system of polynomial equations. Among the solution of that system is the solution of the estimation problem. The 2-D and the 3-D cases are explicitely studied and simplified equations are given for those cases. In addition results are given concerning the optimal choice of the Lagrange multiplier.

d'un ensemble de couples de points de dimension n lorsque sont utilisées des distances de Mahalanobis. Ce problème est énnoncé comme un problème de minimization sous contraintes et la technique des multiplicateurs de Lagrange est utilisée. La solution de ce problème est donnée sous la forme d'un système d'équations polynomiales. Parmi les solutions de ce système se trouve la solution du problème d'estimation. Les cas 2-D et 3-D sont étudiés explicitement et des systèmes d'équations simplifiés sont présentés pour ces cas précis ainsi que des résultats d'optimalité concernant le choix des multiplicateurs de Lagrange.

Mots-clés : alignement de points, transformation rigide, distance de Mahalanobis, ICP, ECMPR, résolution de systèmes polynomiaux 1 Estimating the optimal rigid body transformation using Mahalanobis distances 1.1 Problem setting

Assume that there exist two corresponded point sets {p i } and {q i } along with a set of Mahalanobis n×n matrices Π i , i = 1 . . . N , such that they are related by :

q i ←-Rp i + t + v i (1)
where R is a standard n×n rotation matrix, t is a n-D translation vector and v i a n-D noise vector.

The objective is to find the optimal transformation [R * , t * ] that maps the set {p i } onto {q i } which minimizes the Mahalanobis least squares error criterion given by :

Σ 2 (R, t) N i=1 v T i Π i v i = N i=1 [q i -(Rp i + t)] T Π i [q i -(Rp i + t)] (2)
For a point Rp i + t, the projection matrix, Π i authorizes, for instance, the point q i to be on a line or a plane since the Mahalanobis projection can assign a small or zero weight to a displacement along such a line or plane. Using Mahalanobis distances is more general, it allows to give n different wheights for the displacements in the n orthogonal directions of any reference system. It is possible to penalize and authorize displacements locally knowing the local geometry of the objects. Therefore the point registration process is more robust since two different points of the same linear subspace, for instance, can be associated. It is also more realistic : the Euclidean distance model assume that it is possible to model the world with points. However this kind of model is true only if an infinite precision and memory is available to encode the world representation. Therefore in the Euclidean approache an incompleteness problem arises du to the finite representation : the points of the second cloud never correspond to the points of the first cloud.

To ensure that R is a rotation matrix it is equivalent that two conditions are fulfilled :

R T R = I (3) det(R) = 1 (4)
The first equation (eq. 3) states that the rotation matrix is orthogonal and the second one (eq. 4) implies that there is a conservation of the angles between figures along the transformation. This is now possible to cast the problem into a constrained minimization framework which Lagrangian is :

L (R, t) = Σ 2 (R, t) + Tr Ä L[R T R -I] ä + λ d (det(R) -1) ( 5
)
where L is a matrix of Lagrange multipliers. Now the objective is to find the zero of the derivatives of L (R, t) with respect to t and R.

1.2 t * as a function of R * As for the Euclidean case, there are no constraints on t. Therefore it is possible to get rid of t by substituting an expression of t * as a function of R * .

∂ L (R, t) ∂ t = ∂ Σ 2 (R, t) ∂ t = -2 j Π j q j -(Rp j + t) Thus, setting ∂ L (R,t)
∂ t to zero leads to the following equation :

j Π j q j -Rp j - j Π j t = 0
and therefore t * can be expressed as a function of R * :

t * = j Π j -1 j Π j q j -R * p j (6)
In the following we note :

K = j Π j -1 (7) l = j KΠ j q j (8)
We remark that as each Π i is symmetric, j Π j is also symmetric and necessarily, K is also symmetric. For the same reasons, K is also semi-definite positive or definite positive wether the Π i 's are semi-definite or definite.

With these notations :

t * = l -K j Π j R * p j . (9) 1.3 Expansion of Σ 2
We drop the star in the notations of t * and R * until the end of the section for clarity. But as we have for a given rotation the optimal translation, we search to minimize a Σ 2 that takes the optimal translation into account.

INRIA

Let us note

x i = q i -(Rp i + t) such that Σ 2 = i x T i Π i x i .
And let us replace t by its expression as a function of R (eq. 9) in x i :

x i = q i -Rp i -l + j KΠ j Rp j . 1.3.1 Expansion of x T i Π i x i Let us note, temporary, m = j KΠ j Rp j ( 10
)
such that :

x i = q i -Rp i -l + m.
Now we can expand x T i Π i x i . We obtain 16 summands that we can group into 10 ones thanks to the fact that x T Π i y = y T Π i x since Π i is symmetric.

The first group of 3 summands does not contain R :

κ i = q T i Π i q i -2q T i Π i l + l T Π i l (11)
Each summand in the second group of 7 contains R and can be rewritten using Kronecker (see 5.1) and vec operators (see 5.2) :

-2q T i Π i Rp i = -2q T i p T i ⊗ Π i vec (R) (12) 2q T i Π i m = 2q T i Π i KM vec (R) (13) (Rp i ) T Π i Rp i = vec (R) T p i p T i ⊗ Π i vec (R) (14) -2(Rp i ) T Π i m = -2 vec (R) T M T K T p T i ⊗ Π i vec (R) (15) 2l T Π i Rp i = 2l T p T i ⊗ Π i vec (R) (16) -2l T Π i m = -2l T Π i KM vec (R) (17) m T Π i m = vec (R) T M T K T Π i KM vec (R) (18)
where

M = j p T j ⊗ Π j . ( 19
)
All these relations are based upon the following facts, using eq. 75 then theorem 7 :

Rp i = vec (Rp i ) = (p T i ⊗ I n×n ) vec (R) . (20)

Yguel

Moreover :

Π i (p T i ⊗ I n×n ) = p T i ⊗ Π i (21) [p i ⊗ I n×n ] p T i ⊗ Π i = p i p T i ⊗ Π i (22)
such that using in addition theorem (5), all the preceding equations are easily derived.

Using theorem (5) and theorem (6), it is possible to rewrite eq. ( 12). Since :

q T i p T i ⊗ Π i = ([p i ⊗ Π i ] q i ) T = vec Π i q i p T i T (23)
eq. ( 12) is :

-2q T i Π i Rp i = -2 vec Π i q i p T i T vec (R) . ( 24
)

Sum over i

To obtain Σ 2 , let sum over i the preceeding expressions :

-summing over i eq. ( 13), i q T i Π i K appears which is [ i KΠ i q i ] T and one recognizes l T (eq. 8) such that :

i 2q T i Π i KM vec (R) = 2l T M vec (R) (25)
-summing over i eq. ( 15), one recognizes M such that :

i -2 vec (R) T M T K T [p i ⊗ Π i ] vec (R) = -2 vec (R) T M T KM vec (R) (26)
-summing over i eq. ( 16), one also recognizes M such that :

i 2l T [p i ⊗ Π i ] vec (R) = 2l T M vec (R) (27)
-summing over i eq. ( 17), one recognizes

K -1 = i Π i such that : i -2l T Π i KM vec (R) = -2l T M vec (R) (28)
-same equality in eq. ( 18) gives :

i vec (R) T M T K T Π i KM vec (R) = vec (R) T M T KM vec (R) (29)
Eq. 27 and eq. 28 make zero ; eq. 29 and eq. 26 makes -vec (R) T M T KM vec (R) such that :

INRIA

Theorem 1 (Σ 2 as a function of vec (R)) It is possible to derive Σ 2 for a given rotation and the associated optimal translation :

Σ 2 = vec (R) T N -M T KM vec (R) + 2 l T M -o T vec (R) + κ where K = ñ i Πi ô -1 l = K i Πiq i M = i p T i ⊗ Πi N = i p i p T i ⊗ Πi o = vec Ç i Πiq i p T i å
and κ is a constant a that does not depend on R. a κ = i κ i as defined in eq. ( 11)

With these notations the translation is :

t * = l -KM vec (R * ) . ( 30
)
Let us note :

A î N -M T KM ó (31) b î o -M T l ó (32)
using these matrices the constrained minimization problem is stated as :

               R * = arg min vec(R) vec (R) T A vec (R) -2b T vec (R) ,
with constraints :

R * T R * = I n×n , det R * = 1 . ( 33
)
One can notice several wothwhile facts :

-A is symmetric, -A is semi-definite positive.

-the gradient of Σ 2 with respect to vec (R) is :

∂ Σ 2 ( vec (R)) ∂ vec (R) = 2 vec (R) T A -2b T (34)

Yguel

Proof :

The first property comes from direct calculus : N is symmetric by definition and as K is symmetric, M T KM is also symmetric and so is A.

The second property comes from the fact that Σ 2 ( vec (R)) ≥ 0. If A were not semi-definite positive, one could find an eigenvector u of A with a negative eigenvalue -µ. Theorem 1 states that

Σ 2 (u) = -µ u 2 -2b T u + κ (35)
which is quadratic with respect to the coefficients of u. Thus it would exist a vector αu colinear with u with sufficiently great length that would made Σ 2 (αu) < 0, since at the limit only the coefficient of u 2 matters. And so we have a contradiction.

The third property comes from differential calculus and the fact that A is symmetric.

✄ ✂ ✁

Theorem 2 (Equivalence with a least-squares formulation) The minimization of Σ 2 as a function of vec (R) has an equivalent least-squares formulation :

           R * = arg min vec(R) F vec (R) -u 2 , with constraints : R * T R * = I n×n , det R * = 1 . with A = F T F b = F T u
Proof :

Let us rewrite :

F vec (R) -u 2 = [F vec (R) -u] T [F vec (R) -u] = vec (R) T F T F vec (R) -2u T F vec (R) + u T u INRIA Therefore minimizing F vec (R) -u 2 with respect to vec (R) is equivalent to minimizing : vec (R) T F T F vec (R) -2u T F vec (R) .
Since A is symmetric semi-definite positive, it can be factorized as F T F, therefore a direct identification with eq. ( 33) gives :

A = F T F , b = F T u . ✄ ✂ ✁
This makes the link with the work in [START_REF] Viklands | Algorithms for Weighted Orthogonal Procrustes Problem and other Least Squares Problems[END_REF] where the general leastsquares problem, presented above, was approached using a numerical method close to that described in the following.

3-D Minimization

For any rotation matrix in 3-D there exists two unit quaterions. Let

ξ =     x y z w     = xi + yj + zk + w a quaternion such that x 2 + y 2 + z 2 + w 2 = 1 then the equivalent rotation matrix is R(ξ) =   1 -2y 2 -2z 2 2xy -2zw 2xz + 2yw 2xy + 2zw 1 -2x 2 -2z 2 2yz -2xw 2xz -2yw 2yz + 2xw 1 -2x 2 -2y 2   (36)
This is now possible to cast the problem into a constrained minimization framework which Lagrangian is :

L (ξ, t) = Σ 2 (R(ξ), t) + λ(ξ T ξ -1) ( 37
)
The problem is far more simpler in terms of the number of equations, since only one Lagrange multiplier is used.

Equivalence with a polynomial system

Knowing the gradient of Σ 2 with respect to vec (R) (eq. 34), we use the chain rule and the Jacobian of the change of variable from ξ to vec (R) :

Yguel D vec(R)(ξ) = 2               0 -2y -2z 0 y x w z z -w x -y y x -w -z -2x 0 -2z 0 w z y x z w x y -w z y -x -2x -2y 0 0               2J , (38)
to find the gradient of Σ 2 with respect to ξ :

∂ Σ 2 (ξ) ∂ ξ = ∂ Σ 2 ( vec (R)) ∂ vec (R) D vec(R)(ξ) = 4 î vec (R(ξ)) T A -b T ó J (39)
Finally, setting the Lagrangian to zero gives the system to solve :

ß 2J T [A vec (R(ξ)) -b] + λξ = 0 1 2 x 2 + y 2 + z 2 + w 2 -1 = 0 (40)
It is a system of polynomial equations of the third degree with five equations and five unknowns including the Lagrange multiplier. This kind of system has no closed form as far as we have investigate. Therefore numerical thechniques must be used to solve it. It is interesting to note that multiplying the first equation of the system (40) by ξ T gives an expression for λ :

λξ T ξ = ξ T 2J T [b -A vec (R(ξ))] λ = [2Jξ] T [b -A vec (R(ξ))] (41)
Substituing eq. ( 41) into the first equation of system (40) gives a polynomial system of the fifth degree with four equations and four unknowns :

J T [A vec (R(ξ)) -b] -(ξ T J T [A vec (R(ξ)) -b]) ξ = 0 . ( 42
) Let φ J T [A vec (R(ξ)) -b] (43)
then eq. ( 42) becomes :

φ -(ξ T φ) ξ = 0 ( 44
)
where ξ T φ has a geometric meaning as the scalar product between ξ and φ.

As ξ has unit norm, a geometric view of the previous equation is that φ and ξ are parallel and have same orientation.

INRIA

Numerical approach

We can rewrite this system as f = 0 and the unknowns as the vector x = x y z w λ

T . Thus the Jacobian of the system, obtained using the chain rule for Hessian matrices, is defined as ∂ f (i) ∂ x(j) i,j is :

ñ 4J T AJ + G + λI 4×4 ξ ξ T 0 ô ( 45
)
where G is a 4×4 matrix whose entries are defined by :

G(i, j) î vec (R(ξ)) T A -b T ó h i,j and
h i,j ∂ vec (R) (ξ) ∂ ξ(i)∂ ξ(j) .
We have h i,j = h j,i such that :

h 1,1 h 1,2 h 1,3 h 1,4 h 2,2 h 2,3 h 2,4 h 3,3 h 3,4 h 4,4               0 0 0 0 -4 0 0 0 -4                             0 2 0 2 0 0 0 0 0                             0 0 2 0 0 0 2 0 0                             0 0 0 0 0 2 0 -2 0                             -4 0 0 0 0 0 0 0 -4                             0 0 0 0 0 2 0 2 0                             0 0 -2 0 0 0 2 0 0                             -4 0 0 0 -4 0 0 0 0                             0 2 0 -2 0 0 0 0 0                             0 0 0 0 0 0 0 0 0              
gives all the h i,j . Given a good guess, the system can be solved using an iterative numerical method such as the Newton algorithm.

Choice of a guess

Ideal case

In the above form, it is clear that the problem has at least one solution. Since the unit quaternions form a sphere which is compact and since Σ 2 is polynomial in ξ and hence continuous, it exists a minimum of the function on the quaternion sphere. It is worthwhile to notice that as Σ 2 is positive, when there is no noise, i.e.

v i = 0, ∀i RR n➦ 6684
in eq. ( 2), the correct rotation matrix achieves a global minimum among all the matrices too. Since this minimum is unique as the solution of the linear equation ∇Σ 2 ( vec (R)) = 0, (34), the correct rotation matrix is also given by vec

(R) = A -1 b (46) R = mat 3×3 A -1 b . (47)

Perturbed case

In presence of noise, these relations do not hold anymore. Nevertheless we can use them to define a guess. It is not possible to use mat 3×3 A -1 b as a guess in general since it is not even a rotation matrix and hence has no quaternion equivalent. However for a given matrix, we know the closest orthogonal matrix for the Frobenius norm [START_REF] Schönemann | A generalized solution of the orthogonal procrustes problem[END_REF]. Much better in [START_REF] Kanatani | Analysis of 3-d rotation fitting[END_REF] is given an explicit algorithm to find the closest rotation matrix that can be adapted to our case : I therefore propose to use as the guess the quaternion

Algorithm 1 -Closest Rotation Matrix A, b // Inputs 2: U, D, V ← the SVD of mat 3×3 A -1 b // mat 3×3 A -1 b = UDV T with D diagonal and U and V orthogonal R ← U   1 0 0 0 1 0 0 0 det(UV T )   V T 4: return R
equivalent to the closest rotation matrix, given by the above algorithm. This choice gives very good results in practice when the noise is reasonable and the matchings are correct.

Remarks

It is worth noticing that in Viklands (2006, paper1, page 41) a smarter algorithm is proposed that have much stronger convergence properties than the standard Newton minimization. It uses a, so called, Newton on manifold procedure that aims to find optimal trajectories on the sphere of rotation matrices. However, that algorithm is also local (it converges also to a local optimum). Moreover, it is slower than the brute force Newton algorithm proposed above for simple cases (when the matchings are almost correct as in a pursuit scenario when initialization is always close to the optimum).

Closed form solution in the 2-D case

Theorem (1) is also valid when the space is 2-D. Also the search of the rotation does not require here to involve quaternions. In 2-D, any rotation matrix R can be written :

INRIA R = ï α -β β α ò ( 48
)
where α is the real part and β is the imaginary part of the complex number :

ξ = α + iβ ( 49
)
such that α 2 + β 2 = 1. This constraint is exactly the determinant constraint stated above. This is now possible to cast the problem into a constrained minimization framework which Lagrangian is :

L E (ξ, t) = Σ 2 E (R(ξ), t) + λ(ξ T ξ -1) (50)
(it has the exact same form than in the 3-D case).

Minimization

The Jacobian matrix of the change of variable from ξ to vec (R) is :

D vec(R)(ξ) =     1 0 0 1 0 -1 1 0     (51)
Therefore using the chain rule and the Jacobian of the change of variable from ξ to vec (R) (eq. 51) in 2-D, the system to solve, here, is :

   D T vec(R)(α,β) [A vec (R(α, β)) -b] + λ ï α β ò = 0 α 2 + β 2 -1 = 0 (52)
By expanding all the matrix calculus, it gives :

   α(x 1 + λ) + βy 1 -z 1 = 0 αy 1 + β(y 2 + λ) -z 2 = 0 α 2 + β 2 = 1 (53)
where

x 1 = A(1, 1) + A(4, 1) + A(1, 4) + A(4, 4) , y 1 = A(2, 1) -A(3, 1) + A(2, 4) -A(3, 4) , z 1 = b(1) + b(4) , y 2 = A(2, 2) -A(3, 2) -A(2, 3) + A(3, 3) , z 2 = b(2) -b(3) .
Now it is possible to obtain α β T as a function of λ :

ï α β ò = 1 (x 1 + λ)(y 2 + λ) -y 2 1 ï -z 2 y 1 + z 1 (y 2 + λ) z 2 (x 1 + λ) -z 1 y 1 ò . (54)

Yguel

With that system of equations, it is clear that z 1 = z 2 = 0 implies α = β = 0 and indeed the problem has no solution. By remplacing α and β into the last equation of the system (53) we obtain :

P (λ) = λ 4 + 2δλ 3 + δ 2 -2∆ 3 -z 2 1 -z 2 2 λ 2 +2(z 1 ∆ 1 + z 2 ∆ 2 -δ∆ 3 )λ + ∆ 2 3 -∆ 2 1 -∆ 2 2 = 0 , (55)
where

δ = x 1 + y 2 , ∆ 1 = z 2 y 1 -z 1 y 2 , ∆ 2 = z 1 y 1 -z 2 x 1 , ∆ 3 = y 2 1 -y 2 x 1 . P (
λ) is a 4 degree polynomial with respect to λ. The problem is solved by finding the roots of P which agrees with the result of Censi ( 2008). λ is necessarily real thus the searched Lagrange multiplier is found by taking the real root that minimizes Σ 2 .

Theorem 3 (Choice of the optimal Lagrange multiplier in the 2-D case) It is not necessary to compute the energy associated with a Lagrange multiplier to obtain the global minimizer because the optimal Lagrange multiplier is always given by the greatest Lagrange multiplier. Let (R1, λ1) and (R2, λ2) be two solutions of the system (52), the link with the associated energy values Σ 2 (R1) and Σ 2 (R2) is given by the following inequality :

if λ1 ≥ λ2 : 1 2 [λ1 -λ2] ≤ Σ 2 (R2) -Σ 2 (R1) (56)
Proof :

The steps of the proof given here are inspired by the proof of theorem (1) in [START_REF] Gander | Least squares with a quadratic constraint[END_REF], unfortunately this theorem does not apply here, even if the problem studied in that work is very close to our case. As D vec(R)(α,β) is constant let us note it D in this proof. The main observation is that :

ï α β ò T D T =     α β -β α     = vec (R(α, β)) T . ( 57
)
Let us note now, x i for vec (R(α i , β i )), i = 1, 2. If x 1 , λ 1 and x 2 , λ 2 are solution of the system (52), then :

D T [Ax 1 -b] + λ 1 ï α 1 β 1 ò = 0 , (58)
D T [Ax 2 -b] + λ 2 ï α 2 β 2 ò = 0 . (59)

INRIA

Multiplying eq. ( 58) by

ï α 1 β 1
ò T and eq. ( 59) by

ï α 2 β 2 ò T
to the left, one obtains using the observation (57) :

x T 1 [Ax 1 -b] + λ 1 (α 2 1 + β 2 1 ) = 0 , x T 2 [Ax 2 -b] + λ 2 (α 2 2 + β 2 2 ) = 0 .
Using the norm constraint :

x T 1 Ax 1 -b T x 1 + λ 1 = 0 , (60)
x T 2 Ax 2 -b T x 2 + λ 2 = 0 . (61)
Now substracting eq. ( 61) from eq. ( 60) gives :

λ 1 -λ 2 = (x 2 T Ax 2 -b T x 2 ) -(x 1 T Ax 1 -b T x 1 ) (62)
Multiplying eq. ( 58) by

ï α 2 β 2
ò T and eq. ( 59) by

ï α 1 β 1 ò T
to the left, one obtains using the observation (57) :

x T 2 Ax 1 -b T x 2 + λ 1 1 2 x T 1 x 2 = 0 , (63)
x

T 1 Ax 2 -b T x 1 + λ 2 1 2 x T 1 x 2 = 0 , (64)
noticing that

ï α 1 β 1 ò T ï α 2 β 2 ò = 1 2 x T 1
x 2 Now substracting eq. ( 64) from eq. ( 63) gives :

- 1 2 x T 1 x 2 (λ 1 -λ 2 ) = b T x 1 -b T x 2 (65)
using the equality x T 2 Ax 1 = x T 1 Ax 2 since A is symmetric. Now if we sum eq. ( 62) and eq. ( 65), we obtain :

(λ 1 -λ 2 ) - 1 2 x T 1 x 2 (λ 1 -λ 2 ) = (x T 2 Ax 2 -2b T x 2 ) -(x T 1 Ax 1 -2b T x 1 ) (66) (λ 1 -λ 2 )[1 - 1 2 x T 1 x 2 ] = (x T 2 Ax 2 -2b T x 2 + κ) -(x T 1 Ax 1 -2b T x 1 + κ) (λ 1 -λ 2 )[1 - 1 2 x T 1 x 2 ] = Σ 2 (x 2 ) -Σ 2 (x 1 ) (67)
Using Cauchy-Schwartz inequality we have :

x T 1 x 2 ≤ x 1 x 2 x T 1 x 2 ≤ 1 1 - 1 2 x T 1 x 2 ≥ 1 2 (68)
Therefore eq. ( 67) and eq. ( 68) implies that λ 1λ 2 and Σ 2 (x 2 ) -Σ 2 (x 1 ) have the same sign. Therefore if λ 1 ≥ λ 2 then Σ 2 (x 1 ) ≤ Σ 2 (x 2 ) and the solution Yguel (x 1 , λ 1 ) is a better minimizer. Moreover we have if

λ 1 ≥ λ 2 , then Σ 2 (x 2 ) -Σ 2 (x 1 ) ≥ 1 2 [λ 1 -λ 2 ] . (56)

✄ ✂ ✁

4 Optimality condition for the 3D case

The complete study of the 2-D case led me to investigate more closely the 3-D case and fortunately, I found that a similar result holds.

Theorem 4 (Choice of the optimal Lagrange multiplier in the 3-D case) The 3-D case is similar to the 2-D case : it is not necessary to compute the energy associated with a Lagrange multiplier to obtain the global minimizer because the optimal Lagrange multiplier is always given by the greatest Lagrange multiplier.

Proof :

The main observation is, here, that :

vec (R(ξ)) -vec (I 3×3 ) = Jξ (69)
Let us note now, x i for vec (R(ξ i )), i = 1, 2. If x 1 , λ 1 and x 2 , λ 2 are solution of the system (40) :

2J T [Ax 1 -b] + λ 1 ξ 1 = 0 , ( 70
) 2J T [Ax 2 -b] + λ 2 ξ 2 = 0 . (71)
thus, following closely the steps of the proof of the theorem (3), multiplying each equation ( 70) and ( 71) by ξ T 1 and ξ T 2 in that order then in cross order ; then substracting gives :

λ 1 -λ 2 = 2(x 2 T Ax 2 -b T x 2 -vec (I 3×3 ) T Ax 2 ) -2(x T 1 Ax 1 -b T x 2 -vec (I 3×3 ) T Ax 1 ) (72) -x T 1 x 2 (λ 1 -λ 2 ) = 2(b T x 1 -b T x 2 ) -2( vec (I 3×3 ) T Ax 1 -vec (I 3×3 ) T Ax 2 ) (73)

INRIA

The next step is to sum eq. ( 72) and eq. ( 73) :

(

1 -x T 1 x 2 )(λ 1 -λ 2 ) = 2[(x T 2 Ax 2 -2b T x 2 ) -(x T 1 Ax 1 -2b T x 1 )] (1 -x T 1 x 2 )(λ 1 -λ 2 ) = 2[Σ 2 (x 2 ) -Σ 2 (x 1 )] (74)
Again we have, thanks to the Cauchy-Schwartz inequality, a usefull inequality : 1 -x T 1 x 2 ≥ 0. Now we can conclude that λ 1λ 2 and Σ 2 (x 2 ) -Σ 2 (x 1 ) have the same sign.

Vec operator

Definition 2 (Vec operator) Let A be a m×n matrix and aj its j-th column ; then vec (A) is the mn×1 vector :

vec (A) =     a1 a2 . . . an    
The operator transforms a matrix into a vector by stacking the columns of the matrix one underneath the other.

There is an obvious, but usefull relation : when p is a vector : vec (p) = p. (75)

INRIA

Mat operator

Definition 3 (Mat operator) The mat operator is used to inverse the vec operator and restore a matrix from a vector. However it is necessary to define the dimensions of the output matrix, since different matrices have the same image by the vec operator. Therefore if v is a n dimensional vector such that n = pq mat p×q (v) = M

where M is a p×q matrix such that vec (M) = v.

Relations between the vec operator and the Kronecker product

Theorem 5 (

 5 Kronecker product) Let A be an m×n matrix and B a p×q matrix. The mp×nq matrix defined by  A(1, 1) B . . . A(1, n) B . . . . . . A(m, 1) B . . . A(m, n) B  is called the Kronecker product of A and B and writtenA ⊗ B.The Kronecker product is defined for any pair of matrices A and B whenever their sizes are. For instance if A = a b c and B = I 3×3 then A ⊗ B = aI 3×3 bI 3×3 cI 3×3 = Transpose and Kronecker product) Let A be a m×n matrix and B be a n×q matrix then :[A ⊗ B] T = A T ⊗ B T

Theorem 6 (

 6 Vectorialization of the product of 3 matrices) Let A, B and C three matrices such that the product ABC is defined. Then vec (ABC) = (C T ⊗ A) vec (B) Theorem 7 (Vectorialization of a matrix product) Let A be a m×n matrix and B be a n×q matrix then : vec (AB) = (B T ⊗ I m×m ) vec (A) = (B T ⊗ A) vec I n×n = (I q×q ⊗ A) vec (B) RR n➦ 6684

  Problem setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 t * as a function of R * . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Expansion of Σ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.1 Expansion of x T i Π i x i . . . . . . . . . . . . . . . . . . . . 5 1.3.2 Sum over i . . . . . . . . . . . . . . . . . . . . . . . . . . 6

	Table des matières
	1 Estimating the optimal rigid body transformation using Maha-
	lanobis distances	3
	1.1

RR n➦ 6684


Cédric Desprez
email: cedric.desprez@ifsttar.fr

Panagiotis Kotronis
email: panagiotis.kotronis@ec-nantes.fr

Jacky Mazars
email: jacky.mazars@inpg.fr

Seismic vulnerability assessment of a RC structure before and after FRP retrofitting

Keywords: Structure, Vulnerability, Assessment, Retrofitting, Multifiber, FRP

In earthquake engineering, reliable vulnerability assessment tools suitable to existing reinforced concrete structures are necessary in order to mitigate the seismic risk. Furthermore, amongst the wide range of technical solutions available for structural upgrading, external reinforcement by fiber reinforced polymer (FRP) is often an interesting option. Nevertheless, the use of FRP is limited, one of the reasons being the lack of predictive numerical tools allowing for vulnerability assessment. Based on a case study, this article presents a simplified modeling strategy to assess the seismic vulnerability of an existing reinforced concrete building before and after FRP retrofitting. The structure is simulated using multifiber beam elements and constitutive laws based on damage mechanics and plasticity, while the dynamic characteristics of the numerical model are validated using in-situ ambient vibration records. Nonlinear transient dynamic analysis studies are performed using a synthetic earthquake signal compatible with the Eurocode 8 spectra. Local indicators are adopted to quantify the damage level in the structure before and after FRP retrofitting, in correspondence with the European Macroseismic Scale 98. One of the main conclusions of this study is that the use of local indicators can lead to contradictory assessment results with the evaluations based on global indicators as adopted in HAZUS or Risk-UE.

Introduction

Mitigation through retrofitting of existing reinforced concrete (RC) structures is an important issue in earthquake engineering. Amongst the wide range of available technical solutions, the use of fiber reinforced polymer (FRP) is sometimes adopted. FRP can be seen as an answer to the lack of existing steel bars in RC elements. In order to facilitate the structural design procedure efficient numerical tools are however needed, able to quantify the structural vulnerability before and after retrofitting.

Due to their simplicity, structural engineers often adopt nonlinear static pushover analysis schemes. Although recent developments improved the accuracy and the applicability of the method, see for example the multimodal analysis considering the influence of several modes in the structural response [START_REF] Chopra | A modal pushover analysis procedure to estimate seismic demands for buildings: theory and preliminary evaluation[END_REF][START_REF] Chopra | A modal pushover analysis procedure to estimate seismic demands for buildings: theory and preliminary evaluation[END_REF], the adaptive pushover that deals with stiffness degradation [START_REF] Antoniou | Advantages and limitations of adaptive and non-adaptive force-based pushover procedures[END_REF][START_REF] Kalkan | Adaptive modal combination procedure for nonlinear static analysis of building structures[END_REF]o rt h e work of [START_REF] Penelis | 3D Pushover analysis: the issue of torsion[END_REF]Kappos (2002)andFajfar et al. (2005) introducing torsional effects, a numerical strategy based on a pushover analysis cannot be used as a general tool. It is not adequate for example for the structural assessment of 3D irregular structures submitted to severe earthquake loadings or when higher order modes dominate the nonlinear structural behavior. For these cases, a nonlinear dynamic analysis is more appropriate.

The main limitations of a refined 3D nonlinear dynamic analysis are the computational cost, the complexity of the finite element mesh and the huge amount of the obtained numerical results. This paper presents a simplified nonlinear dynamic analysis strategy to assess the vulnerability of an existing RC structure before and after FRP retrofitting. The numerical model is done using tri-dimensional multifiber Timoshenko beam elements [START_REF] Guedes | A fibre Timoshenko beam element in CASTEM[END_REF]) and constitutive laws based on damage mechanics (La Borderie 1991; La [START_REF] Borderie | Damage mechanics model for reinforced concrete structures under cyclic loading[END_REF]) and plasticity [START_REF] Menegotto | Method of analysis of cyclically loaded reinforced concrete plane frames including changes In: geometry and non-elastic behaviour of elements under combined normal force and bending[END_REF]. The use of beam elements leads to a reduced number of degrees of freedom and thus to acceptable computational times. The combination of multifiber beams for the structural elements (beams, columns, walls) and advanced nonlinear constitutive models for concrete and steel has been proven very efficient in recent research studies [START_REF] Mazars | Numerical modelling for earthquake engineering: the case of lightly RC structural walls[END_REF]Kotronis et al. 2005;[START_REF] Grange | Numerical modelling of the seismic behaviour of a 7-story building: NEES benchmark[END_REF][START_REF] Grange | The effects of soil-structure interaction on a reinforced concrete viaduct[END_REF]. Comparison of experimental with the numerical results showed that this simplified modeling strategy is able to reproduce global (forces/displacements) but also local results (strains), especially when bending behavior is predominant (as for the case of the studied building).

The Grenoble City Hall in France (GCH) is chosen as case study. In order to validate the dynamic characteristics of the numerical model, results coming from the permanent monitoring of the building are used [START_REF] Michel | Vulnérabilité sismique de l'échelle du bâtiment à celle de la ville[END_REF][START_REF] Michel | Full-scale dynamic response of an RC building under weak seismic motions using earthquake recordings, ambient vibrations and modeling[END_REF]. The structural assessment is performed with and without FRP retrofitting following the Eurocode 8 (EC8) prescriptions. Local quantities (damage variables and plastic strains) are adopted to quantify the damage level in the structure in accordance with the European Macroseismic Scale (EMS 98). The HAZUS (AT C -4 0 1996) and the Risk-UE method [START_REF] Risk-Ue | An advance approache to earthquake risk scenarios with applications to different European towns[END_REF] are finally used and the results are compared with the ones coming from the local analysis.

Simplified modeling strategy

The Grenoble City Hall (GCH) (Fig. 1) is a RC building constructed in 1966. It is composed of two distinct parts, a tower and a peripheral building. Only the tower is studied hereafter. Grenoble's area is classified as a moderate seismic zone, subject to local site effects. The seismic level is medium with a nominal acceleration between 1.6 and 3 m/s 2 according to the French seismic mapping [START_REF] Planseisme | French seismic intensity areas[END_REF].

Geometrical and material characteristics

The tower is 52 m high, 43 m long and 13 m wide (Fig. 1). The upper part is an eleven-story frame supported by a large transfer deck with pre-stressed concrete beams (Fig. 2a). The lower part of the tower is composed of a 12 m height hall and an underground level. The main structural frame of the tower is made of four RC piers (cores) containing the lift shafts and stairwells (Fig. 2b,c). From the basement to the deck the piers are of "Square type" and of "U type" above the deck till the top. Except from the top level (11th) which contains RC walls, all the other levels (1st to 10th) follow a beam-column frame design. The foundations of the tower are RC circular deep piles.

The material characteristics of the GCH being unknown, the adopted numerical values (Table 1) are based on the usual values corresponding to the 1960-1970 period.

Numerical meshing

In order to decrease the number of the necessary degrees of freedom (DOF), 3D Timoshenko multifiber beam elements are used for the spatial discretization (Fig. 3) [START_REF] Guedes | A fibre Timoshenko beam element in CASTEM[END_REF]; Figure 4a shows the finite element mesh of the tower. Nonlinear multifiber beam elements are used for the piers, beams and columns. The slabs and the walls at the top level (11th) are discretized using linear elastic Kirchhoff plate elements. Linear beam elements are adopted for the pre-stressed transfer deck. The total number of DOF is 58848. Details about the discretization of the different parts are provided hereafter.

Piers, columns and beams

In a typical level, all piers and columns are modeled using 5 Timoshenko multifiber beams. In the underground and the first level above the deck where a higher demand level is expected, 15 elements are used to improve accuracy. In order to decrease the necessary number of fibers, the section of the piers is simplified (Fig. 4b); openings are not precisely discretized but the total area and flexural inertia are equal to the exact ones. The sections of the beams and columns are identical to the real ones. The multifiber sections are discretized using 5-20 cm edge triangular elements representing the concrete fibers and square elements representing the steel fibers (Fig. 4c). Each steel fiber represents a 5 steel bars group with an area of approximately 5 cm 2 .

Transfer deck

Although the GCH is relatively well documented, details about the structural design of the transfer deck are unknown. The numerical model of the transfer deck is therefore based on architectural plans and consists of four longitudinal beams (2 m × 0.6 m) and three transverse beams (2 m × 1.15 m) (Fig. 4b). Due to the pre-stressed effect an elastic linear behavior is assumed.

Top story

Considering the specific design of this story that contains RC walls and presents a higher stiffness, linear elastic Kirchhoff plate elements are introduced in the numerical model (Fig. 4d).

Boundary conditions

Due to the existence of pile foundations, the numerical model is assumed fixed at the base. Ambient vibrations measurements provided a verification of this assumption for small dynamic loadings [START_REF] Michel | Full-scale dynamic response of an RC building under weak seismic motions using earthquake recordings, ambient vibrations and modeling[END_REF]).

Constitutive laws

A constitutive model suitable for cyclic loadings is associated at each fiber of the multifiber elements. The La Borderie damage mechanics law is used for concrete (La Borderie 1991; La [START_REF] Borderie | Damage mechanics model for reinforced concrete structures under cyclic loading[END_REF]) and the Menegotto & Pinto plasticity model for steel [START_REF] Menegotto | Method of analysis of cyclically loaded reinforced concrete plane frames including changes In: geometry and non-elastic behaviour of elements under combined normal force and bending[END_REF]. The concrete model takes into account the opening and closing of cracks while the steel model is able to consider the buckling of the reinforced bars. Shear is considered linear. A short description of the La Borderie and Menegotto-Pinto models is presented hereafter.

La Borderie's constitutive model

This concrete model is based on damage mechanics and takes into account the crack opening and closing (Fig. 5) during cyclic loading. The model's general formulation is tridimensional (3D), but only the uniaxial (1D) version is used herein, according to the multifiber modeling needs. An overview of the adopted parameters is presented in Table 2. The total strain (ε) is defined as the sum of an elastic (ε e ) and an anelastic part (ε a )(Eqs.1, 2,and3).

ε = ε e + ε a (1) ε e = σ + E(1 -D 1 ) + σ - E(1 -D 2 )
(2)

ε a = β 1 .D 1 E (1 -D 1 ) F ′ (σ ) + β 2 .D 2 E (1 -D 2 )
(3) σ + and σ -are, respectively, the tensile and compressive stresses (Eqs. 4 and 5); E is the Young's modulus; and β 1 and β 2 are material constant parameters controlling the anelastic strains in tension and compression, respectively;

if σ>0thenσ + = σ and σ -= 0( 4 )

if σ<0thenσ + = 0andσ -= σ (5)
F ′ (σ ) is a function that controls crack opening and closing (Eq. 6); -σ f is the stress crack closure value (Fig. 5).

F ′ (σ ) = dF dσ (6)
if σ ≥ 0thenF(σ ) = σ ⇒ F ′ (σ ) = 1 if -σ f ≤ σ<0thenF(σ ) = σ 1 + σ 2σ f ⇒ F ′ (σ ) = 1 + σ σ f if σ ≤-σ f then F(σ ) = -σ f 2 ⇒ F ′ (σ ) = 0
D 1 and D 2 (Eq. 7) are the damage variables due to tension and compression respectively, varying from 0 (no damage) to 1 (completely damaged material). D i is driven by Y i (Eqs. 8 and 9). Y 0i is the initial damage threshold and A i and B i are constants directly identified from uniaxial tension and compression tests.

D i = 1 - 1 1 + [A i (Y i -Y oi )] B i with i = 1 (tension) or i = 2 (compression) (7) Y 1 = (σ + ) 2 2E(1 -D 1 ) 2 + β 1 .F(σ ) E(1 -D 1 ) 2 (8) Y 2 = (σ -) 2 2E(1 -D 2 ) 2 + β 2 .σ E(1 -D 2 ) 2 (9)

Menegotto-Pinto constitutive model

The cyclic behavior of the steel bars is simulated using the classical Menegotto-Pinto model with cinematic hardening (Fig. 6) and possibility to take into account buckling by introducing a negative modulus slope in compression depending on the transverse steel reinforcement spacing. The monotonic behavior is defined through the initial Young's modulus (E s ), the plastic threshold (ε sy ,σ sy ), the ultimate strength and strain (ε u ,σ u ) and the yielding slope (E h ). The unloading and reloading process, is guided by analytic relations (Eqs. 10-13) corresponding to a set of curves ranging between the elastic and the yielding asymptotes. An overview of the adopted parameters is presented in Table 3.

A fac = 0.006 A6 fac = 620 ρ fac = 20
Stiffness ratios between cyclic yielding and elastic regime B fac = 0.00357

C fac = 0.5
Ratio between transverse steel reinforcement (TSR) spacing and steel bar diameter

f ald = 5
Volumetric mass (Kg/m 3 ) ρ s = 7, 800

σ * = bε * + 1 -b (1 + (ε * ) R ) 1/R ε * (10) σ * = σ s -σ r σ o + σ r (11) ε * = ε s -ε r ε o + ε r (12) R = R o - A i ξ A j + ξ Tension : i = 1and j = 2; Compression : i = 3and j = 4 (13) (σ s ,ε s ) is the studied point; (σ o ,ε o )
is the crossing point of the elastic and yielding slopes; (σ r ,ε r ) are the coordinates of the previous point of load reversion; b is the E h /E s ratio; R is a shape parameter; ζ is the ratio between the maximum reached strain during loading; ε o , R o , A i and A j are material constants that can be obtained from experimental results. A detailed description of the model is presented in [START_REF] Menegotto | Method of analysis of cyclically loaded reinforced concrete plane frames including changes In: geometry and non-elastic behaviour of elements under combined normal force and bending[END_REF].

Numerical model validation

The structure's ambient vibrations due to wind, road traffic or small seismic ground motions, are monitored using a permanent accelerometers network. Six tri-directional accelerometric stations are used inside the building and one outside. The response of the structure in terms of accelerations, velocities and displacements is known and thus its main natural frequencies and stiffness evolution [START_REF] Michel | Vulnérabilité sismique de l'échelle du bâtiment à celle de la ville[END_REF][START_REF] Michel | Full-scale dynamic response of an RC building under weak seismic motions using earthquake recordings, ambient vibrations and modeling[END_REF].

In order to validate the finite element model, a modal analysis is performed and the numerical natural frequencies are compared to the ones coming from the ambient measurements. It is found that up to the fifth frequency the numerical model reproduces correctly the insitu measurements (Fig. 7) [START_REF] Michel | Full-scale dynamic response of an RC building under weak seismic motions using earthquake recordings, ambient vibrations and modeling[END_REF]. A low motion earthquake record (Fig. 8)is also used to compare the response of the numerical model with the in-situ ambient vibration records. The signal is tri-directional with a maximal peak ground acceleration of 0.023 m/s 2 in the horizontal direction. Comparisons of the maximum inter-story drifts (Fig. 9a) and the Fourier transforms of the top displacement signals (Fig. 9b) show again a good correlation [START_REF] Michel | Full-scale dynamic response of an RC building under weak seismic motions using earthquake recordings, ambient vibrations and modeling[END_REF]. 4), a synthetic accelerogram is created (Fig. 10a) using empiric green functions [START_REF] Causse | Calibrating median and uncertainty estimates for a practical use of empirical Green's functions technique[END_REF]). This signal is tri-directional (uncoupled) and fits the EC8 design spectra (Fig. 10b). The maximum peak ground acceleration is approximately equal to 3 m/s 2 in the two horizontal directions and 2m/s 2 in the vertical direction.

Nonlinear dynamic analysis

The response of the structure under the synthetic tri-directional EC8 earthquake is studied hereafter in terms of global (displacement, drift) and local (damage, strain) quantities.

Flexural behavior

The maximum horizontal top displacements are found equal to 132.3×10 -3 and 127×10 -3 m respectively in the longitudinal (X) and transverse (Y) directions (Fig. 11). The vertical displacement in the center of the top slab doesn't exceed 1 cm, a relative small value for a building of 28.7 m long, leading to small strain values. This last result comforts the choice of a linear elastic constitutive law for the slabs.

Stiffness is visibly higher in the lower part of the tower (rectangular piers) (Fig. 12a), something that is certainly due to the change of the piers' section (Fig. 2). The flexural behavior of the structure is obvious looking at the maximum inter-story drift diagram (Fig. 12b). Drifts The maximum shear stress is found equal to 1.36 MPa at the base of the "U-type" piers, a value lower than the Eurocode 2 prescriptions limits. The use of Timoshenko multifiber beam elements with a linear shear assumption seems thus adequate.

The structure bends mainly according to the first and second flexural modes, (60 and 30 % of the effective modal mass respectively). The inflexion point is situated just above the transfer deck.

Torsion effects

Despite the 3D loading and the small dissymmetry of the building (differences between the stairwell and the lift shaft cores) no significant torsional effect is observed. Moreover, as shown in Fig. 13 the torsion moment at the base of the piers is very small in comparison with the two bending moments.

Damage distribution in concrete

The two damage variables D 1 and D 2 of the concrete constitutive model are directly used as indicators of damage due to traction and compression respectively. The two variables vary between 0 (no damage) and 1 (totally damaged material).

Although damage due to traction is not representative of the building capacity, as it appears for low levels of loading and evolves fast, it gives an idea of the areas of the structure where the seismic energy is dissipated. The distribution of the D 1 variable at the final step of the calculation is represented in Fig. 14a. Tensile macro-cracks (D 1 > 0.9) are localized at the base of the piers and at the first four levels above the deck. Its history evolution indicates that energy is dissipated mainly during the first 10 s of the signal.

High values of the compression damage variable D 2 (D 2 > 0.9) are representative of concrete crushing, spalling and possible reinforcement bucking. The distribution of the D 2 variable at the final step is found small, with a maximum value limited at D 2 =0.25andis located at the base of the "U-type" pier. Compression damage is thus not significant for the chosen earthquake signal.

Strains in steel bars

Figure 14b shows the strains envelop in the steel bars of the pier P4. Two yielding zones are apparent: the first is located at the base of the rectangular pier. The reached value is small, but is situated in a critical area for the stability of the building. The maximum strain value (4.3 × 10 -3 ) is reached in the first level above the deck. Even if these yielding values are not excessive, it is important to take into account that the building, constructed at the 60 s, does not comply with the capacity design philosophy.

3.3 Damage level assessment: initial RC structure HAZUS (ATC-40 1996)andRisk-UELM2(RISK 2003) compute damage level from global indicators, as the inter-story drift or the top displacement, depending of the building typology (ie: materials, geometry…). Although these methods are widely used, they are not well adapted for the specific building where the dynamic behavior is modified (with the FRPstrengthening) without significant changes in the building geometry. In this study, the chosen modeling strategy allows obtaining global but also local indicators as damage in concrete and yielding of the steel bars. For this reason, the adopted assessment methodology in the following is based on equivalence between the European Macro-seismic Scale classification (EMS 1998) and the local variables according to the work of Lang [START_REF] Lang | Seismic vulnérability of existing building[END_REF]), (Table 5; Fig. 15). The results of the HAZUS (ATC-40 1996) and the Risk-UE methods [START_REF] Risk-Ue | An advance approache to earthquake risk scenarios with applications to different European towns[END_REF] are also presented in the last section of the article in order to highlight the contradictions that can be found using global and/or local variables.

According to the results of the nonlinear dynamic analysis, concrete cracking appears at the underground level and at the four first levels above the deck. Reinforcement bars yield 5, the damage level can be therefore identified as Moderate (level 2/5).

Vulnerability assessment after FRP retrofitting

A FRP retrofitting solution is studied hereafter with the objective, because of the particular importance of the GCH, to decrease the damage expectance level from moderate (2/5) to negligible (1/5) (Table 5;Fig. 15).

The nonlinear dynamic analysis presented in Sect. 3 highlights two critical areas in the structure: the base of the tower subject to maximal bending moments, and the first four floors above the deck where the geometry of the piers changes from "Rectangular-type" to "U-type". In these areas, concrete cracks (tension damage) and reinforcement bars yield. No important values of compression damage are observed and shear stresses are acceptable. The bending modes are mainly excited and torsional effects are negligible.

Following the previous remarks, the adopted design solution is to retrofit (numerically) the four main piers (P1 to P4, Fig. 2b) in order to increase their flexural capacity and to avoid steel bar yielding. This can be achieved by bonding FRP on both sides of the piers as FRP has a function similar to an external reinforcement. Other retrofitting solutions, such as columns confinement, are not studied hereafter as damage due to compression or shear is not significant.

FRP design solution no. 1

The chosen retrofitted areas are the four main piers at the underground level and at the four levels above the deck (Fig. 16). The FRP location takes into account the opening access in the piers. 5 mm thick sheets are introduced numerically on the four levels above the deck, and 4 mm thick sheets at the basement level. The thickness of the FRP sheets is calculated using the numerical interaction diagram (normal forces-bending moments) at the base of the rectangular and "U" piers (Fig. 17) as follows: The results coming from the EC8 synthetic signal are compared to the "service design limit" corresponding to the steel yield threshold (0.2 %) and 60 % of the concrete strength in compression (18 MPa). Adding FRP in the pier section allows defining a new "service design limit with FRP". The thickness is chosen such as all the numerical values are situated inside the new interaction diagram. Table 6 shows the retained FRP properties (typical values are used).

Modeling strategy for FRP retrofitting

In a finite element model based on multifiber beams, FRP can be simply introduced using additional fibers in the beam section (Fig. 18). Perfect bonding is assumed (no delamination) and a linear elastic brittle constitutive law is adopted for the additional fibers [START_REF] Dufour | Modélisation du comportement dynamique d'une structure à murs porteurs en béton armé renforcé à l'aide de tissu à fibres de carbone[END_REF]LESSLOSS 2006;[START_REF] Desprez | Analyse et réduction de la vulnérabilité sismique des structures existantes: renforcement par collage de tissus de fibres de carbone (TFC)[END_REF]. The comparison of the experimental results with the numerical ones for a monotonic loading proves the performance of the approach. It is important to notice however that this FRP modeling strategy is well adapted for flexural reinforcement but not for FRP column confinement. For this last case, a novel simplified modeling strategy and a new global constitutive law can be found in [START_REF] Desprez | Damage model for FRP-confined concrete columns under cyclic loading[END_REF].

Nonlinear dynamic analysis of the retrofitted structure

A modal analysis of the elastic retrofitted structure showed that the main natural frequencies remained unchanged. This is due to the small amount of the FRP sheets used. Nevertheless, as it is shown hereafter, the nonlinear response of the building is found extensively modified. More specifically, Fig. 19a, b shows that the floor displacements and inter-story drifts in the longitudinal direction (X) are higher for the retrofitted than for the non-retrofitted configuration. The transverse direction (Y) displacements and drifts remain almost unchanged. In both directions, the maximum values of curvatures are now shifted upper on the structure (Fig. 19c).

The concrete damage distribution D 1 (Fig. 20a) shows a more extended area of damage due to tension in comparison with the non-retrofitted structure (Fig. 14). Tensile cracks are significant at the base of the tower and above the deck and they reach the fifth level. Damage distribution due to compression is similar with the non-retrofitted case, with a small increase on the maximal value (D 2 = 0.3).

The addition of FRP limits the strains and avoids the creation of plastic hinges inside the reinforced areas (Fig. 20b). However, new plastic zones appear just above the retrofitted elements. The discontinuity of the FRP sheets leads to new areas of strain location in the structure and thus to new sources of potential vulnerability. This clearly indicates the necessity of performing nonlinear dynamic analysis to choose the optimum retrofitting design solution.

A Fast Fourier Transform (FFT) of the top displacement time histories for the retrofitted and the non-retrofitted structures highlights a shift in frequencies that is less important when FRP is introduced (Fig. 21). The retrofitted structure is subject to more significant accelerations (Fig. 22) and to higher internal forces that result to increased maximum displacements values.

Damage level assessment: FRP design solution no. 1

At the final step of the calculation, the structure presents cracks at the base and on the first five levels above the deck. A moderate yielding of the steel bars (ε max = 2.617 × 10 -3 ) appears above the retrofitting areas (Fig. 20b) while no concrete compression failure or steel bar buckling is observed. It is also important to notice that strains in FRP reinforcement are small and thus the risk of debonding, delamination or failure is limited. According to the classification presented in Table 5, the damage level is however In order to achieve this goal, two new optimized design solutions are studied in the following sections.

Retrofitting design optimization

Damage level assessment: FRP design solution no. 2

Figure 23 shows the nonlinear dynamic response of the structure in the X direction for the case of the second design solution, where FRP are placed till the fifth level above the deck. The maximal displacements and inter-stories remain similar to the results of the first retrofitting design solution while the strain localization zones are simply shifted at the end of the retrofitted area (Fig. 24;Table7). The damage assessment level is still "Moderate" (level 2/5) using local indicators and the EMS98 scale.

Damage level assessment: FRP design solution no. 3

From the previous results, it is obvious that the position of the plastic hinges is driven by the discontinuities between the retrofitted and the non-retrofitted areas. For the design solution no. 3, the FRP's thicknesses are adjusted considering the internal forces at each level. This leads to a progressive decrease of the FRP thickness of 0.5 mm at each level. Furthermore, 7). Nevertheless, the progressive development of the curvature in the structure (Fig. 25c) limits the maximum strain values and cracks in concrete, allowing a more uniform energy dissipation (Fig. 26a); a smooth pattern of the strain localization zones is achieved (Fig. 26b). D 1 indicator is now distributed from the base of the building to the sixth level, no yielding of reinforced bars is observed and no cracks in concrete due to compression. The small D 2 increase comes from a better use of the compressed concrete part. According to the methodology of Table 5 the damage level in the EMS 98 scale is now considered "Negligible" (Level 1/5). Building type : concrete shear wall (C2)-high (H)

Local versus global indicators

Table 7 presents an overview of the main characteristics of the structural response for all the different FRP design solutions. Top displacements and inter story-drifts (global variables) are found significantly increased in X direction whereas the damage level, based on the local state of materials, is found decreased.

In the following sections, the damage levels found using local indicators (Table 7)a r e compared to the damage levels from global indicators (HAZUS and Risk-UE).

Damage assessment using HAZUS

According to the HAZUS methodology (ATC-40 1996), the damage level of the building is defined from the maximum inter-story drift. Table 8 present the HAZUS damage levels for the GCH, which is classified hereafter as a high level concrete building with shear walls.

Considering the initial GCH structure, without FRP, the maximum inter-story drifts in X and Y directions (Table 7) lead to a Moderate damage level (Table 8). After the first FRP With uy = 0.9 u -0.7 y ; U=ultime; y=yield (1) Trend evolution of damage indicators from the previous tested case (2) Damage indicator location considering for example the HAZUS method and the first FRP design case, the maximal drift is 5.937 × 10 -3 (Table 7), leading to a Moderate damage level (Tables 8). The lower and upper thresholds are respectively 3.0 × 10 -3 and 8.0 × 10 -3 . Considering the lower threshold as 0 % and the upper one as 100 % of the possible progression in the considered damage level, the 5.937 × 10 -3 value is thus situated at 58 % strengthening, inter-story drifts are increased in the X direction, and are almost unchanged in the Y direction. According to the HAZUS assessment, the damage level remains Moderate. The same conclusions are drown for the second and third retrofitting configurations.

Considering the local approach however, results are different for the 3rd retrofitting scheme, where the damage level is classified as Slight (Table 10). Furthermore using the HAZUS approach, the retrofitted configurations lead to an increase of the maximal inter-story drift and thus we get closer to the upper damage threshold. This trend is opposite to the results obtained with the local approach, where damage decreases with FRP retrofitting (Table 10). This contradiction on the damage assessment based on the HAZUS and the local methods can lead to significant errors.

Damage assessment using Risk-UE

To apply the Risk-UE LM2 (RISK 2003) methodology, pushover analyses are needed. The damage levels thresholds depend then on the top displacement at "elastic" and "ultimate" limits of the building (Table 9). The static target displacement is then commonly used to define the reached damage level. However, having a deeper knowledge of the structure motion, the damage level is here defined from the maximal top displacement of the structure during dynamic analyses. In this study, the damage thresholds have been defined using pushover analysis on the GCH detailed in [START_REF] Desprez | Analyse et réduction de la vulnérabilité sismique des structures existantes: renforcement par collage de tissus de fibres de carbone (TFC)[END_REF], available for the initial structure and the FRP designs 1 and 3.

Considering the initial GCH structure, without FRP, the maximum top displacement in X and Y directions (Table 7) lead to an Important damage level (Table 9). After the first FRP strengthening, top displacements are increased in the X direction, and are slightly decreased in the Y direction. Due to the higher capacity of the building however, the Risk-UE damage level thresholds are higher in the retrofitted configuration and the damage level becomes Moderate. The same conclusions are drown for the second and third retrofitting configurations. The Risk-UE damage indicator trend for the studied configurations is in accordance with the results of the local analysis: with better retrofitting design, the damage index decreases (Tables 7 and10). However, considering the resulting classifications for the non-retrofitted and the 3rd retrofitting scheme results are not identical (Important vs. moderate, moderate vs. slight respectively, see Table 10). The use of global or local quantities can therefore lead to contradictory conclusions depending on the adopted assessment methodology.

A new indicator for the damage state assessment

Considering the local indicators as the accurate way for structure assessment, it is however still necessary to provide efficient global indicators. Indeed, in a post-seismic situation, global indicators can be used in concordance with in-situ monitoring; allowing to quickly evaluate the damage state of a structure. Global indicators are also the most practicable way for vulnerability assessment in a design context, as for the engineering field.

In this article, it is shown that the use of FRP leads to changes in the dynamic behavior of the structure without significantly impacting its geometry (the location of the plastic hinges and the values of the internal forces can change). The use of global indicators for the damage assessment, based on defined building typologies, can however lead to inaccurate diagnostics. It seems that a new building typology should be created regarding to the "FRP retrofitted RC building"; nevertheless, this would be very difficult as there is no unique retrofitting configuration. Another solution could be to define a new global indicator bridging the gap between local and global fields. The curvature could play a part of this role as it is directly linked with the local strains (Figs. 19,20,25).

Conclusions and perspectives

This paper presents a simplified method to assess the vulnerability of a RC structure before and after FRP retrofitting:

• The structure is simulated using multifiber beam elements and uniaxial constitutive laws for concrete and steel taking into account the cyclic behavior of the materials. • The comparison between ambient vibrations records and the numerical response shows the good performance of the model.

• The structure is submitted to a synthetic signal in accordance with the Eurocode 8 requirements and the results of the nonlinear dynamic analysis are presented. • Local indicators as concrete damage variables and steel strains are used to assess the damage level of the building following a methodology based on the EMS 98 and Lang [START_REF] Lang | Seismic vulnérability of existing building[END_REF] research work. • FRP retrofitting is added to the numerical model and optimized solutions are studied in order to decrease the vulnerability level of the structure.

The study highlights the advantages of a multifiber beam modeling strategy to perform 3D dynamic nonlinear numerical calculations with low computational demands.

Fig. 1

 1 Fig. 1 CGH: view of the tower and the peripheral building

Fig. 2

 2 Fig. 2 GCH. a Transverse view of the structure, b structural frame of a typical level: 1st to 10th and c structural frame of the hall

Fig. 3

 3 Fig. 3 Modeling strategy-multifiber beam concept

Fig. 5

 5 Fig. 5 La Borderie model-Uniaxial stress-strain relation

Fig. 6

 6 Fig. 6 Menegotto & Pinto model-Uniaxial stress-strain relation

Fig. 7

 7 Fig. 7 GCH modeling-modal analysis; in-situ versus model

Fig. 9

 9 Fig. 9 CGH response to the low motion signal. a Maximum value of the inter-story drift; in-situ versus model and b Fourier transform of the top displacement signal; in-situ versus model

Fig.

  Fig. Synthetic EC8 earthquake. a Tri-directional signal and b signal spectra versus EC8 design spectra

Fig.

  Fig. CGH response to the synthetic EC8 earthquake-Displacements history at the top of the structure

Fig. 13

 13 Fig. 13 CGH response to the synthetic EC8 earthquake-time history of the bending (Mx, My) and torsion (Mz) moments at the base of the structure-Pier P1

Fig. 14

 14 Fig. 14 CGH response to the synthetic EC8 earthquake. a Distribution of D 1 (tension damage), final step and b Envelop of strains in the steel bars of the pier P4; yielding threshold at 0.2 %

Fine

  Fig.15EMS 98 damage level(EMS 1998)

Fig. 16

 16 Fig. 16 GCH (FRP design no. 1)-Location of FRP sheets: a in the structure, b in a "U-type" pier and c in a "square-type" pier

Fig. 17

 17 Fig. 17 GCH (FRP design no. 1)-Interaction diagram: normal force (Nz) versus flexural moment (Mx) at the base of a "U-type" pier (P4)

Fig. 18

 18 Fig. 18 Introducing additional FRP fibers in a multifiber beam section

Fig. 19

 19 Fig. 19 GCH (FRP design no. 1)-retrofitted versus non-retrofitted configuration. Response of the pier P4 in the X direction: a displacements, b inter-story drifts, c curvatures

Fig. 20

 20 Fig. 20 GCH (FRP design no. 1)-a Distribution of D 1 (tension damage), final step and b envelop of strains in the steel bars of the pier P4 (yielding threshold at 0.2 %)

Fig. 21

 21 Fig. 21 GCH (FRP design no. 1)-FFT of the top displacement time histories-with and without FRP retrofitting Vs. elastic model: a longitudinal direction X and b transverse direction Y

Fig. 23

 23 Fig. 23 GCH (FRP design no. 2)-FRP strategy no. 1 versus FRP strategy no. 2. Response of the pier P4 in the X direction: a displacements, b inter-story drifts, c curvatures

Fig. 24

 24 Fig. 24 GCH (FRP design no. 2)-a Distribution of D 1 (tension damage), final step and b envelop of strains in the steel bars of the pier P4; FRP strategy no. 1 versus FRP strategy no. 2

  Fig. 25 GCH (FRP design no. 3)-optimized retrofitting versus non-retrofitted configuration. Response of the pier P4 in the X direction: a displacements, b inter-story drifts, c steel strains d curvatures

Fig. 26

 26 Fig. 26 GCH (FRP design no. 3). a Distribution of D 1 (tension damage), final step and b envelop of strains in the steel bars of the pier P4; Optimized retrofitting versus non-retrofitted configuration

Table 1

 1 Concrete and steel characteristics

	Concrete	Compression strength	Tension strength		Young modulus
		30 MPa	3 MPa		32 GPa
	Steel	Yield limit	Young modulus	Failure	Strain at failure
		400 MPa	200 GPa	461 MPa	9 %

Table 3

 3 Steel model parameters

	Parameters	Tension	Compression
	Young's modulus (GPa)	E a = 200
	Poisson ratio	ν s = 0.3
	Elastic limits	σ sy = 400 MPa, ε sy = 0.002
	Yielding threshold	ε sh = 0.03
	Ultimate limits	ε su = 0.09	σ su = 461 MPa
	Material coefficients	A1 fa = 18.5	A2 fa = 0.15

Table 4

 4 GCH-parameters for the generation of the synthetic EC8 earthquake

	Area seismic level	Medium
	Building class	IV
	Soil class	C

Table 5

 5 Equivalence between EMS98 and local variables[START_REF] Lang | Seismic vulnérability of existing building[END_REF]

	EMS 98 damage level	Identification
	1-Negligible to slight damage (no structural damage,	Point of onset of cracking
	slight non-structural damage)

Table 6 GCH

 6

	properties	retrofitting-FRP	Thickness of the layers t f	Ultimate strength f fu	Ultimate strain ε fu	Young modulus
			mm	MPa		E f GPa
			1	1,400	0.0133	105

Table 7

 7 GCH-maximum responses in the pier P4 for the different FRP design solutions Indicator (Max) Non-retrofitted structure FRP design no. 1 FRP design no. 2 FRP design no. 3

	Top disp X (m)	132.3 × 10 -3	161.6 × 10 -3	162.6 × 10 -3	163.5 × 10 -3
	Top disp Y (m)	127.0 × 10 -3	116.3 × 10 -3	116.0 × 10 -3	115.6 × 10 -3

Table 9

 9 Risk-Ue method-damage level for the Grenoble City Hall Damage level Maximum top head displacement : (m)

		RISK-UE	Non-retrofitted	FRP	FRP	FRP
		(LM2)-	structure	design	design	design
		Threshold		no. 1	no. 2	no. 3
	Slight	>0.7 y	> 26 × 10 -3	> 26 × 10 -3 -	> 26 × 10 -3
	Moderate	>0.7 y + 0.05 uy > 40 × 10 -3	> 61 × 10 -3 -	> 90 × 10 -3
	Important	>0.7 y + 0.2 uy > 81 × 10 -3	> 165 × 10 -3 -	> 282 × 10 -3
	Ve r y	>0.7 y + 0.5 uy > 163 × 10 -3	> 373 × 10 -3 -	> 666 × 10 -3
	Important
	Destruction	>0.7 y + uy	> 300 × 10 -3	> 719 × 10 -3 -	> 1, 305 × 10 -3

Table 10

 10 GCH-damage level comparison-local approach, HAZUS approach and Risk-UE approach

	Damage assessment	Non-retrofitted	FRP design no. 1	FRP design no. 2	FRP design no. 3
		structure
	Local method
	Damage Level	Moderate	Moderate	Moderate	Slight
	Trend (1)		Decreasing	Decreasing	Decreasing
	HAZUS method
	Damage level	Moderate	Moderate	Moderate	Moderate
	Damage indicator	30 %	58 %	57 %	57 %
	location (2)
	Trend (1)		Increasing	Almost unchanged	Unchanged
	Risk-UE method
	Damage level	Important	Moderate	-	Moderate
	Damage indicator	63%	97%	-	38%
	location (2)
	Trend (1)	-	Decreasing	-	Decreasing

  Comparison between several retrofitting design solutions shows that they can lead to different damage localization zones. The use of different methodologies for vulnerability assessment based on global (HAZUS, Risk-UE) or local(EMS 1998;[START_REF] Lang | Seismic vulnérability of existing building[END_REF] indicators can provide contradictory results. It is therefore necessary to develop concordant indicators. Future work should be focused on studying several earthquake levels and various structures in order to establish a good correlation between local and global damage indexes. Although local variables are most representative of the effective state of a structure, global indicators (top displacements and inter story-drifts) are easier to use. A simplified way to deal with this could be to adopt curvatures as complementary indicators. They can be determined from the first derivative of the inter-story drifts, or the second derivative of the displacements. For flexion dominated structural elements strains are directly linked to curvatures, something that allows bridging the gap between local and global variables.

Acknowledgments

The authors are grateful for the financial support of the French ANR (National Agency for Research) program ARVISE (Analyse et Réduction de la vulnérabilité sismique du bâti existent), Projet ANR-2006-PGCU-007-01. The in-situ records and the synthetic signal used in the dynamics analysis were provided by Phillipe Gueguen and Clotaire Michel (Institut des Sciences de la Terre (ISTERRE), http://isterre. fr).


Hang Le

Florentin Barbier

Ha Nguyen

Natalia Tomashenko

Salima Mdhaffar

Souhir Gahbiche

Bougares Fethi

Benjamin Lecouteux

Didier Schwab

Yannick Estève

Fethi Bougares

ON-TRAC' systems for the IWSLT 2021 low-resource speech translation and multilingual speech translation shared tasks

Introduction

In the two last editions of the IWSLT evaluation campaigns, the ON-TRAC consortium focused on end-to-end offline speech translation and simultaneous speech translation [START_REF] Nguyen | On-trac consortium end-to-end speech translation systems for the iwslt 2019 shared task[END_REF][START_REF] Elbayad | On-trac consortium for end-to-end and simultaneous speech translation challenge tasks at iwslt[END_REF]. In 2021, we chose to focus on low-resource speech translation and multilingual speech translation by using two different kinds of approaches: a cascaded speech-to-text translation (combining source language automatic speech recognition (ASR) and source-to-target text translation) to process the low resource speech translation tasks, and a neural end-to-end model for the multilingual speech translation task. For the low resource task, we investigated the use of speech features extracted by a neural model pretrained by self supervision the wav2vec XLSR-53 model [START_REF] Conneau | Unsupervised cross-lingual representation learning for speech recognition[END_REF] in order to process Swahili lan-guages by a classical hybrid Markovian/neuronal ASR system. The ASR outputs were processed by neural machine text-to-text translation systems dedicated to the two targeted language pairs. For the multilingual speech translation task, we investigated the use of a dual-decoder Transformer that jointly transcribes and translates an input speech. This model was trained in order to translate from multiple source languages to multiple target ones.

The ON-TRAC Consortium is composed of researchers from three French academic laboratories and an industrial partner: LIA (Avignon Université), LIG (Université Grenoble Alpes), LIUM (Le Mans Université), and researchers from Airbus.

Low resource speech translation

The task of the low resource speech translation track was to build the speech transcription/translation system for transcribing and/or translating between the two language pairs:

• Coastal Swahili (swa) to English (eng)

• Congolese Swahili (swc) to French (fra)

ASR system

The same ASR models were used for both test datasets: Coastal Swahili (swa) and Congolese Swahili (swc).

Data

The training corpus for the ASR acoustic model (AM) comprises of several datasets:

• 5k instances for Congolese Swahili speech provided by the IWSLT-2021 organizers1 ;

• a training subset of the ALFFA corpus [START_REF] Gelas | Developments of swahili resources for an automatic speech recognition system[END_REF] (read speech and broadcast news);

• a subset of the IARPA Babel Swahili Language Pack 2 (conversational and scripted telephone speech that spoken in the Nairobi dialect region of Kenya).

The total size of the training corpus is about 74 hours. In our preliminary experiments, we also tried to include a swa dataset (5k instances of Coastal Swahili), provided by the IWSLT-2021 organizers, into the training corpus, but this does not improve the ASR performance. Hence, for the submitted system and for the results reported in the paper, this corpus was not used.

Architecture

In this work, we investigated the impact of using self-supervised learning [START_REF] Baevski | wav2vec 2.0: A framework for self-supervised learning of speech representations[END_REF] on the hybrid ASR HMM/DNN acoustic models, as well as on the pipeline ASR+MT system performance. Self-supervised learning (SSL) has shown to be effective for various speech-related tasks including ASR and MT [START_REF] Schneider | wav2vec: Unsupervised pre-training for speech recognition[END_REF][START_REF] Baevski | wav2vec 2.0: A framework for self-supervised learning of speech representations[END_REF][START_REF] Evain | Lebenchmark: A reproducible framework for assessing self-supervised representation learning from speech[END_REF][START_REF] Nguyen | Investigating Self-Supervised Pre-Training for End-to-End Speech Translation[END_REF] and could be especially beneficial for a low-resource scenario.

We trained several acoustic models (AM) with two different types of input features for comparison:

(1) 40-dimensional high-resolution (hires) MFCC features; and (2) wav2vec 2.0 features [START_REF] Baevski | wav2vec 2.0: A framework for self-supervised learning of speech representations[END_REF] extracted by the multilingual large model pretrained pretrained in 53 languages, XLSR-53-large [START_REF] Conneau | Unsupervised cross-lingual representation learning for speech recognition[END_REF].

The phoneme set and transcriptions were the same as in the work [START_REF] Gelas | Developments of swahili resources for an automatic speech recognition system[END_REF].

The AMs are state-of-the-art factorized time delay neural networks (TDNN-F) [START_REF] Povey | Semi-orthogonal low-rank matrix factorization for deep neural networks[END_REF][START_REF] Peddinti | A time delay neural network architecture for efficient modeling of long temporal contexts[END_REF] and were trained using the Kaldi toolkit [START_REF] Povey | The kaldi speech recognition toolkit[END_REF]. The models have similar topology (except for the input features): 12 TDNN-F layers (1,024-dimensional, with projection dimension of 128) and a 2232-dimensional output layer. The AMs were trained using lattice-free maximum mutual information (LF-MMI) [START_REF] Povey | Purely sequence-trained neural networks for asr based on lattice-free mmi[END_REF] and cross-entropy criteria. Speed and volume perturbations have been applied for data augmentation. 100-dimensional speaker i-vectors were appended to the input features.

We used a 3-gram LM with a 466K vocabulary provided in the ALLFA recipe [START_REF] Gelas | Developments of swahili resources for an automatic speech recognition system[END_REF] 3 .

2 IARPA-babel202b-v1.0d, https://catalog.ldc. upenn.edu/LDC2017S05

3 https://github.com/getalp/ALFFA_ PUBLIC

Neural machine translation system

In order to translate the ASR outputs from source languages to target languages, two neural machine translation systems were built.

Data

For the swa-eng sentence pairs, training dataset for machine translation system includes:

• OPUS4 english-swahili parallel data :

CCAligned and MultiCCAligned (El-Kishky et al., 2020), WikiMatrix, Wikimedia, XLEnt and ParaCrawl.

• 5k parallel swa-eng dataset provided by IWSLT-2021.

The total size of the training dataset for swa-eng is about 3.2M sentence pairs. We applied language identification filtering LangID [START_REF] Lui | langid.py: An off-the-shelf language identification tool[END_REF] keeping only swa-eng sentence pairs with correct English. Sentence pairs where the English side is detected as noisy are removed from the swa-eng training dataset. In total, we filter out about 30% of the original training set and obtains a dataset of 2.2M sentence pairs. As for swc-fra NMT system, training data includes parallel corpora made available by the organizers in addition to the available corpora for this language pair on OPUS website. Overall we used a training set of 1.1 M sentence pairs.

Architecture

We propose an NMT model using long shortterm memory neural networks (LSTMs) (Hochreiter and Schmidhuber, 1997). NMT systems for swa-eng and for swc-fra were trained using the lstm luong wmt en de model template, a standard LSTM Encoder-Decoder architecture with Luongstyle attention [START_REF] Luong | Effective approaches to attention-based neural machine translation[END_REF]. Swa-eng system was built at the subword level using a joint BPE vocabulary of 32768 BPE unit, trained using source and target language. Swc-fra NMT model, on its side, was trained at the word level.

Results

The ASR results in terms of word error rate (WER) are reported in Table 1 on the development datasets for different types of acoustic features. We can see that using wav2vec features significantly decreases the WER and provides about 8% of relative WER reduction for both datasets. The MT results in terms of BLEU [START_REF] Papineni | Bleu: a method for automatic evaluation of machine translation[END_REF] score are reported in Table 3. Notice that while the WER of the outputs of the ASR fed by wav2vec features is lower than the one fed by MFCC features, for the swc-fra language pair, the BLEU score of the translation from the MFCCbased ASR system is higher than the one got on the wav2vec-based ASR. By lack of time, we did not yet investigate the reason of this, but we will do as soon as possible.

Features

Dev swa-eng Dev swc-fra

MFCC

Multilingual speech translation

Speech-to-text translation (ST) consists in translating a speech utterance in a source language to a text in another target language (e.g., English audio to French text). In this section, we describe a multilingual ST system that can translate from multiple source languages to multiple target ones.

Data

The data provided for the multilingual ST task is a subset of the Multilingual TEDx corpus (Salesky

Features

Test swa-eng Test swc-fra wav2vec 2.0 12.9 9.1 et al., 2021), in which there are four source languages (Spanish (es), French (fr), Portuguese (pt), and Italian (it)) and five target languages (the aforementioned source languages plus English (en)).

The sizes of the ASR talks range from 107 hours (it) to 189 hours (es). Translation data is part of the ASR talks for a given source language. Our experiments were performed in the constrained setting where only the provided data for the task is used.

Model architecture

Our system is based on the Dual-decoder Transformer [START_REF] Le | Dualdecoder transformer for joint automatic speech recognition and multilingual speech translation[END_REF] which consists of an encoder and two decoders. This architecture jointly transcribes and translates an input speech. Each of the decoders is responsible for one task (ASR or ST) while interacting with each other. We refer the reader to the paper for further details. We initially followed [START_REF] Le | Dualdecoder transformer for joint automatic speech recognition and multilingual speech translation[END_REF] and used 12 encoder layers, 6 decoder layers, and a hidden dimension of d = 256. However, this model produced poor results. We hypothesize that with this configuration, the model capacity is too large for the dataset described in the previous section. In the end, we ended up using only 6 encoder layers and 3 decoder layers (with the same d = 256).

In addition, we also trained a Transformer model having the same encoder of 6 layers but with only one decoder as the baseline (hereafter called singledecoder model).

Implementation details

For text pre-processing, we normalize the punctuation and build the vocabulary on the concatenation of the transcript and translation text using Senten-cePiece [START_REF] Kudo | Sentencepiece: A simple and language independent subword tokenizer and detokenizer for neural text processing[END_REF] without pretokenization. We used 10k unigram vocabulary as it performed slightly better than a vocabulary of 8k tokens in our preliminary experiments. The speech features are 80-dimensional log Mel filterbank. Utterances having more than 3000 frames are removed for GPU efficiency. We used SpecAugment [START_REF] Park | Specaugment: A simple data augmentation method for automatic speech recognition[END_REF] with Librispeech double (LD) policy for data augmentation.

es-en es-fr es-pt es-it fr-en fr-es fr-pt pt-en pt-es it-en it-es For the target-forcing mechanism, we prepended a language-specific token to the target sequence [START_REF] Inaguma | Multilingual end-toend speech translation[END_REF][START_REF] Le | Dualdecoder transformer for joint automatic speech recognition and multilingual speech translation[END_REF]. In order to provide good initialization for our multilingual ST system, we separately trained a multilingual ASR system and a multilingual MT one on the allowed data. We then used the weights from the pre-trained ASR encoder, ASR decoder and MT decoder to initialize our ST encoder, ASR decoder, and ST decoder, respectively. We also used the obtained multilingual MT model to augment the training data by translating the transcripts to the target languages as well as translating the translations back to the source languages.

Our model was trained for 150 epochs using the Adam optimizer (Kingma and Ba, 2015) with the inverse square root scheduler. We averaged the last 10 checkpoints and used beam search with a beam size of 5 for decoding. The results reported are detokenized case-sensitive BLEU [START_REF] Papineni | Bleu: a method for automatic evaluation of machine translation[END_REF]. Our implementation is based on the FAIRSEQ S2T toolkit [START_REF] Wang | Fairseq S2T: fast speech-to-text modeling with fairseq[END_REF].

Results

Table 5 displays the results on the dev and hidden test sets. One can observe that the Dual-decoder Transformer outperforms the baselines of single decoder on all language pairs except for the pt-es direction where it is surpassed by the single-decoder models. The use of transcripts as additional languages [START_REF] Antonino | One-to-many multilingual end-toend speech translation[END_REF] in the single-decoder model improves the results for 4 out of 11 language pairs. Since we aim to obtain a single end-to-end multilingual ST system that can perform many-tomany translation, we selected the Dual-decoder Transformer for our final submission.

Conclusion

This paper described the ON-TRAC consortium submissions to the low-resource translation task and to the multilingual speech translation task. Our unique ASR system for both Swahili and Congolese Swahili languages uses XLSR-53 wav2vec features as speech representation input. It got the best results on both Swahili languages (respectively 31.25% and 36.75% of WER). The NMT systems used to translated these transcription into respectively to English and to French got BLEU scores of 12.9 (swa→eng) and 9.1 (swc→fr). The Dualdecoder Transformer we used in the multilingual speech translation got promising results. We did not try a specific strategy to handle language pairs without training data. The low results we got on such language pairs confirm that a specific treatment must be applied in these conditions.

Table 2 :

 2 Table2shows the official ASR results on the test datatests. For our submissions, we used wav2vec features only. These ASR results are the best ones among the ASR results submitted by the participants to this task. Official ASR performance (WER,%) for the test datasets of the low-resource task.

	Features	Dev swa Dev swc
	MFCC hires	19.90	29.57
	wav2vec 2.0	18.34	27.29
	Table 1: ASR performance (WER,%) for the develop-
	ment datasets of the low-resource task.
	Features	Test swa Test swc
	wav2vec 2.0	31.25	36.75

Table 3 :

 3 MT performance (BLEU) for the development datasets of the low-resource task.

	hires	13.39	9.60
	wav2vec 2.0	14.19	9.45
	reference text	18.36	14.07

Table 4 :

 4 Official MT performance (BLEU) on the test datasets of the low-resource task for the submitted system.

Table 5 :

 5 BLEU on the dev and hidden test sets. denotes the use of the transcripts in training.

	Training data (hours)	69	11	42	11	50	38	25	59	-	-	-
	Results on the dev sets
	Single decoder	11.56 4.95 19.14 17.98 13.27 12.99 12.21 13.54 8.31 3.88 4.54
	Single decoder	12.75 5.32 21.95 16.82 11.27 12.15 11.01 12.37 10.25 2.24 2.50
	Dual-decoder	18.59 8.02 25.38 19.22 17.81 17.79 15.20 17.63 3.00 4.09 4.81
	Official results on the hidden test sets
	Dual-decoder	20.20 8.20 25.60 11.10 14.40 15.00 14.90 13.20 3.00 4.20 4.60

https://iwslt.org/2021/low-resource

https://opus.nlpl.eu

Acknowledgments

This work was funded by the French Research Agency (ANR) through the ON-TRAC project under contract number ANR-18-CE23-0021.