Jian Han
email: jian.han@ec-nantes.fr

Christophe Dano
email: christophe.dano@ec-nantes.fr

Pierre-Yves Hicher

Zhen-Yu Yin
email: zhenyu.yin@gmail.com.

Strain-Rate Dependency of Shear Strength for A Highly Overconsolidated Clay

The paper summarizes the results from undrained triaxial compression tests on Merville clay, a highly overconsolidated clay at Merville in the north of France. The specimens were prepared in four different overconsolidation ratios (OCR=7, 14, 28, 56), and for each OCR, undrained shear tests were performed using three or four axial strain rates varying from 0.26 %/h to 25.89 %/h. The strain-rate dependency of the undrained stress level q/q f , the pore pressure, the effective stress path and the stress ratio q/p' were investigated. The undrained shear strength was found to increase with the strain rate. The shear band formation during loading was observed for all strain rates. For each OCR test series, a linear relationship between the normalized peak deviatoric stress and the strain rate was obtained in the semilogarithm plane. As OCR increasing, the slope of best-fitting line for the linear relationship increases. Therefore, the strain rate influence on the normalized undrained shear strength is greater for higher OCRs.

INTRODUCTION

The effect of strain rate on mechanical properties of soils has important engineering implications. In particular, parameters inferred from in situ tests (estimated strain rate a  a  =10 3 -10 5 %/h) and standard laboratory measurements (conducted at a  a  =0.5-5 %/h) must be carefully considered during the selection of design parameters fo Díaz-Rodríguez et al. 2009).

Various studies reported in the literature show that different soils, for example, sands, weak rocks and soft clays, have different responses to changes in strain rate [START_REF] Casagrande | Effect of rate of loading on the strength on clays and shales at constant water content[END_REF][START_REF] Lo | Strength anisotropy and time effects of two sensitive clays[END_REF][START_REF] Lefebvre | Rate effects and cyclic loading of sensitive clays[END_REF]. Some comprehensive researches on reconstituted clay have been done by [START_REF] Sheahan | Rate-dependent undrained shear behavior of saturated clay[END_REF] and [START_REF] Zhu | Strain-rate-dependent stress-strain behavior of overconsolidated Hong Kong marine clay[END_REF], with OCR varying from 1 to 8. Progress in quantifying and understanding undrained strain rate effects requires additional comprehensive data on a variety of soil types. There are not many reliable data available on the time-dependent behaviour of stiff clays with different OCRs, particularly natural stiff clays. Some laboratory tests for London clay with different strain rates have been conducted by [START_REF] Sorensen | Influence of structure on the time-dependent behaviour of a stiff sedimentary clay[END_REF]. But the influence of OCR on the strength behaviour was not analyzed.

Thus, this paper focuses on investigating the influence of strain rate and OCR on the strength behaviour of a natural stiff clay.

EXPERIMENTAL PROGRAM

Description of the clay testedSs

Merville clay was deposited in the Flanders Region which are now located at the northern France, Belgium and the southeast of England. It was covered by tertiary sedimentation. The soil surface then stood probably about 200 m above the present surface of the clay. The formations overlying the clay and the top of it were later eroded. This erosion process was sustained to Quaternary era, when Flandrian alluvial soils were deposited on the clay remaining in place which has thus an effective vertical stress much lower than that was applied to the Pliocene. Thus, the clay is known as a stiff, heavily overconsolidated clay [START_REF] Josseaume | Engineering properties of the Flanders clay at Dunkerque and Calais[END_REF]. This study deals with the undrained stress-strain behaviour of Merville clay sampled at a depth greater than 3 meters (thickness of the clayey silt layer above the tested clay). Undisturbed samples at a depth of 9-11 m were selected. The clayey specimens have a uniform appearance in the form of a dark gray material. Table 1 shows typical physical properties of the tested clay.

Test program and procedure

Tests performed in this study include 1D compression test at constant strain rate and isotropically consolidated undrained triaxial compression (denoted CIUC) tests of overconsolidated Merville clay sheared under a constant strain rate.

1D compression test was done using a Wykeham Farrance testing machine which is a strain controlled machine. The range is 0.00001-9.99999 mm/min. CIUC tests were carried out using a computer-controlled GDS triaxial testing system.

The specimen for 1D compression test was 70 mm in diameter and 19 mm in high. The axial strain rate used in this test was 0.63 %/h (0.002 mm/min). During loading, the drainage was permitted from both sides of the specimen. The strain rate was chosen small enough so that no excess pore pressure was generated during loading.

The specimens for undrained triaxial compression tests were 35 mm in diameter and 70 mm in high. Four different axial strain rates were selected: 0.26 %/h, 0.52 %/h, 2.57 %/h and 25.89 %/h. The B value was checked and found to be greater than 96%. Specimens were isotropically consolidated at different effective consolidation pressures (=50 kPa, 100 kPa, 200 kPa and 400 kPa, respectively) for more than 72 h. After the isotropic consolidation, the cell pressure was held constant, and specimens were sheared to failure at a constant rate of axial strain. The pore pressure was measured at the bottom of the specimen.

When carrying out triaxial tests, we used filter paper strips as side drains which were placed on the lateral surface of each specimen. The purpose was to ensure that the pore water redistribution was completed within a shorter period of time, so that even at a strain rate of 25.89 %/h the uniformity of the pore pressure within the specimen could be assured.

TEST RESULTS

One-dimensional compression test

Fig. 1 shows the result of one-dimensional compression test, giving the void ratio e versus the logarithm of the vertical effective stress ' v .

Fig. 1. One-dimensional compression test on Merville clay.

From Fig. 1, the compression index C c was measured equal to 0.711, and the swelling index increases (C s2 > C s1 ) as the material is subjected to high stresses due to structural damage.

As mentioned above, the undisturbed samples are from the depth of 9-11 m, which means that the in situ vertical effective stress varied from 100 kPa to 150 kPa because of the variation of ground water table. The preconsolidation pressure obtained from this one-dimensional compression test, following the Casagrande method, is approximately 2750 kPa, giving an overconsolidation ratio of more than 18.

Overview of CIUC tests

In this study, four series of undrained triaxial compression tests with OCR=7, 14, 28, and 56 were performed. Table 2 presents the results of tests on overconsolidated specimens. It gives the pre-shear consolidation stress condition (effective consolidation pressure ' 0 , actual OCR and back pressure u b ), and stress-strain properties at the peak deviatoric stress (failure) and at maximum excess pore pressure during shearing. Values of interest for comparative study in this paper are presented in Table 2, where u m is the maximum excess pore pressure,  au is the axial strain corresponding to u m , q f is the peak deviatoric stress, and  aq and u f are the axial strain and excess pore pressure corresponding to q f , respectively. For each test, the axial strain  au is less than the  aq , which means the maximum excess pore pressure u m is reached earlier than the peak shear stress q f .

Tests at OCR=7

Fig. 2a presents the relationship between stress level and axial strain for the test series with OCR=7. Strain localization and shear band formation were observed within these overconsolidated specimens even at the lowest strain rate. The axial strain corresponding to the peak stress increases with an increase in the strain rate. Fig. 2b presents the relationship between excess pore pressure level and axial strain. The axial strain corresponding to the maximum pore pressure increases with an increase in the strain rate.

The relationship between stress ratio and axial strain is shown in Fig. 2c, where p' = (' 1 +2' 3 )/3 is the mean effective stress. These three curves almost coincide during the first stages of the loading and then diverge before the peak of the stress-strain curve. Fig. 2d shows that higher strain rates correspond to larger peak stresses. The best-fitting line for the relationship between the normalized peak stress and the axial strain rate can be represented by a straight line in a semi-logarithmic plane.

Tests at OCR=14, 28 and 56

Figs. 3 to 5 present the relationship between stress level and axial strain (a), the relationship between excess pore pressure level and axial strain (b), the relationship between stress ratio and axial strain (c), and the relationship between normalized peak stress and axial strain rate (d) for the test series with OCR=14, 28 and 56, respectively. Strain localization and shear band formation were observed in these overconsolidated specimens for the tests at the three strain rates. The values of the axial strain corresponding to the peak stress and the maximum pore pressure are almost equal. Figs. 3a to 5a show that the stress-strain curves are almost identical before reaching the peak stress, which means that the evolution of stress level is almost independent of the strain rate before reaching the peak stress. The same phenomenon is observed in Figs. 3b to 5b for the evolution of the pore pressure.

Figs. 3d to 5d show that higher strain rate leads to higher peak stress. The bestfitting line for the relationship between the normalized peak stress and the axial strain rate can also be represented by a straight line in a semi-log diagram.

Analysis of the test results

Fig. 6a presents the effective stress paths at different confining pressures (50 kPa, 100 kPa, 200 kPa and 400 kPa) and at different strain rates (0.26 %/h, 2.57 %/h and 25.89 %/h). We obtained three parallel failure lines for the three different strain rates (Fig. 6a). Thus, the apparent friction angle, ', appears independent of the strain rate.

Its mean value is 24°.

Fig. 6b shows the relationships between the normalized undrained shear strength (q f /' 0 ) and the axial strain rate ( a  a  ) summarized from Fig. 2d, Fig. 3d, Fig. 4d, and Fig. 5d. For each OCR, the set of test points is fitted with a straight line in a semi-log diagram, and the equations of these lines are shown in the figure. The slope of each line reflects the rate of the undrained strength increase with the strain rate. To characterize the influence of the strain rate on the shear strength, [START_REF] Zhu | Strain-rate-dependent stress-strain behavior of overconsolidated Hong Kong marine clay[END_REF] suggested that the undrained shear strength varies almost linearly with the logarithm of the strain rate, expressed as follows:

    f 0 / ' log q a q           a        (1)
where q  is the slope of each straight line shown in Fig. 2d For OCR=7, 14, 28 and 56, the values of the slopes are 0.172, 0.263, 0.418 and 0.717, respectively. Fig. 6c shows the relationship between q  and OCR, in which q  appears to increase with an increase in OCR. It can be concluded that for Merville clay the strain rate influence on the normalized undrained shear strength is greater for higher OCRs.

CONCLUSIONS

A total of thirteen undrained triaxial tests on isotropically consolidated Merville clay with OCR up to 56 were carried out at different axial strain rates. The influence of OCR and strain rate on the undrained shear strength behaviour was studied.

The following conclusions can be drawn from this study:

(1) For undrained triaxial tests, higher strain rate results in higher undrained shear strength for OCRs ranging from 7 to 56 for Merville clay within the strain-rate range tested. Strain localization and shear band formation develop systematically during shearing at any strain rate.

(2) There are three parallel failure lines for the three different strain rates. Thus, the apparent friction angle of Merville clay appears independent of the strain rate. Its mean value is equal to 24°.

(3) For each OCR test series, the relationship between the normalized peak strength and the axial strain rate can be fitted by a straight line in a semi-logarithmic plot.

(4) For each OCR test series, the value of the strain-rate coefficient, q  , calculated from Eq.( 1) was determined by the slope of the fitting line shown in Fig. 6b. The value of q  increases with OCR increasing. Hence, it can be concluded that for Merville clay the strain rate influence on the normalized undrained shear strength is more pronounced for higher OCRs.

Further experiments will be carried out on reconstituted samples, in order to obtain a better understanding of the strain-rate effect and the influence of OCR in intact and reconstituted samples of Merville clay.  and OCR.

Fig. 2 .

 2 Fig. 2. Tests at OCR=7: (a) Stress level versus axial strain. (b) Excess pore pressure level versus axial strain. (c) Stress ratio versus axial strain. (d) Normalized peak shear stress versus axial strain rate.

Fig. 4 .

 4 Fig. 3. Tests at OCR=14: (a) Stress level versus axial strain. (b) Excess pore pressure level versus axial strain. (c) Stress ratio versus axial strain. (d) Normalized peak shear stress versus axial strain rate.

Fig. 5 .

 5 Fig. 5. Tests at OCR=56: (a) Stress level versus axial strain. (b) Excess pore pressure level versus axial strain. (c) Stress ratio versus axial strain. (d) peak shear stress versus axial strain rate.

  Fig. 6. (a) Undrained effective stress paths at different strain rates. (b) Relationship between normalized peak shear stress and axial strain rate. (c) Relationship between q

Table 1 . Physical Properties of Undisturbed Merville Clay.

 1

	Liquid limit	97.21%
	Plastic limit	38.61%
	Plasticity index (PI)	58.62%
	Natural water content	32.52.2%
	Degree of saturation	100%
	Specific gravity of solids	2.650.2%
	Percent finer than 2 m	26%
	m)	2.250.08

Table 2 . Results of CIUC Tests on Stiff Merville Clay.

 2

					At maximum excess	At peak deviatoric
	OCR	' 0 (kPa)		u b (kPa)	pore pressure  au u m	 aq	stress q f	u f
					(%)	(kPa)	(%)	(kPa)	(kPa)
	(a) OCR=7
	6.86	401	0.26	500	2.31	168	2.61	631	164
	6.88	400	2.57	500	2.54	168	3.14	690	160
	6.88	400	25.89	500	3.06	194	3.64	769	185
	(b) OCR=14
	13.82	199	0.26	501	2.99	117	3.49	427	107
	13.75	200	2.57	500	3.09	117	3.33	484	116
	13.75	200	25.89	501	3.30	116	3.70	532	114
	(c) OCR=28
	27.50	100	0.26	700	2.37	82	3.70	371	65
	27.23	101	0.52	699	2.38	82	3.20	387	78
	27.23	101	2.57	700	2.47	89	2.96	415	81
	27.78	99	25.89	701	2.61	78	3.10	455	74
	(d) OCR=56
	53.92	51	0.26	499	1.60	41	2.92	312	35
	55.00	50	2.57	500	1.55	37	2.50	337	34
	56.12	49	25.89	501	2.13	37	2.69	383	31

a  a  (%/h)

ACKNOWLEDGMENTS

This research was financially supported by China Scholarship Council (CSC), the ANR project SOLCYP. These supports are greatly appreciated.


Alexandre Munnier
email: alexandre.munnier@univ-lorraine.fr

Karim Ramdani
email: karim.ramdani@inria.fr

Calderón cavities inverse problem as a shape-from-moments problem

Keywords: Inverse conductivity problem, shape reconstruction, Laplace equation, integral equations, moments problem

In this paper, we address a particular case of Calderón's (or conductivity) inverse problem in dimension two, namely the case of a homogeneous background containing a finite number of cavities (i.e. heterogeneities of infinitely high conductivities). We aim to recover the location and the shape of the cavities from the knowledge of the Dirichlet-to-Neumann (DtN) map of the problem. The proposed reconstruction method is non iterative and uses two main ingredients. First, we show how to compute so-called generalized Pólia-Szegö tensors (GPST) of the cavities from the DtN of the cavities. Secondly, we show that the obtained shape from GPST inverse problem can be transformed into a shape from moments problem, for some particular configurations. However, numerical results suggest that the reconstruction method is efficient for arbitrary geometries.

Introduction

Let be given a simply connected open bounded set Ω in R 2 with Lipschitz boundary Γ. Let σ be a positive function in L ∞ (Ω) and consider the elliptic boundary value problem:

-∇ • (σ∇u f ) = 0 in Ω (1.1a)

u f = f on Γ. (1.1b)
Calderón's inverse conductivity problem [START_REF] Calderón | On an inverse boundary value problem[END_REF] is to recover the conductivity σ knowing the Dirichletto-Neumann (DtN) map f -→ ∂ n u f |Γ of the problem. In a recent work [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF], the authors investigated this problem in the particular case of piecewise conductivity with infinitely high contrast (see for instance Friedman and Vogelius [START_REF] Friedman | Identification of small inhomogeneities of extreme conductivity by boundary measurements: a theorem on continuous dependence[END_REF] who considered this problem in the case of small inclusions). Combining an integral formulation of the problem with tools from complex analysis, they proposed an explicit reconstruction formula for the geometry of the unknown cavity. However, due to the crucial use of the Riemann mapping theorem, the proposed approach was limited to the case of a single cavity. The aim of this paper is to investigate the case of a multiply connected cavity. More precisely, we suppose that Ω contains a multiply connected domain

Ω Γ ω 1 ω 2 ω 3 γ 1 γ 2 γ 3
ω = ∪ N k=1 ω k
, where the open sets ω k , for k = 1, . . . , N are non intersecting simply connected domains with C 1,1 boundaries γ k and ω ⊂ Ω (see Figure 1). We denote by γ = ∪ N k=1 γ k and by n the unit normal to Γ ∪ γ directed towards the exterior of Ω \ ω.

For every f in H 1 2 (Γ), let (u f , c f ) ∈ H 1 (Ω \ ω) × R N , with c f := (c f 1 , . . . , c f N ) T , be the solution of the Dirichlet problem:

-∆u f = 0 in Ω \ ω (1.2a) u f = f on Γ (1.2b) u f = c f k on γ k , k = 1, . . . , N, (1.2c)
with the additional circulation free conditions:

γ k
∂ n u f dσ = 0, k = 1, . . . , N.

(1.2d)

By following the proof given in the Appendix of [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF] for the case of a single cavity (N = 1), it can be easily shown that this elliptic transmission problem is well-posed and that its solution can be seen as the limit solution obtained by considering problem (1.1) for a piecewise constant conductivity and letting the constant conductivity inside the cavities tend to infinity (at the same speed).

The inverse problem investigated in this paper can be formally stated as follows (the exact functional framework will be made precise later on): knowing the Dirichlet-to-Neumann (DtN) map Λ γ : f -→ ∂ n u f , how to reconstruct the multiply connected cavity ω ? Roughly speaking, one can distinguish in the literature two classes of approaches for shape identification: iterative and non iterative methods (see for instance the survey paper by Potthast [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF]). In the first class of methods, one computes a sequence of approximating shapes, generally by solving at each step the direct problem and using minimal data. Among these approaches, we can mention those based on optimization [START_REF] Borcea | Electrical impedance tomography[END_REF][START_REF] Capdeboscq | Imaging by modification: numerical reconstruction of local conductivities from corresponding power density measurements[END_REF], on the reciprocity gap principle [START_REF] Kress | Nonlinear integral equations and the iterative solution for an inverse boundary value problem[END_REF][START_REF] Ivanyshyn | Nonlinear integral equations for solving inverse boundary value problems for inclusions and cracks[END_REF][START_REF] Cakoni | Integral equation methods for the inverse obstacle problem with generalized impedance boundary condition[END_REF], on the quasi-reversibility [START_REF] Bourgeois | A quasi-reversibility approach to solve the inverse obstacle problem[END_REF][START_REF] Bourgeois | The "exterior approach" to solve the inverse obstacle problem for the Stokes system[END_REF] or on conformal mapping [START_REF] Akduman | Electrostatic imaging via conformal mapping[END_REF][START_REF] Kress | Inverse Dirichlet problem and conformal mapping[END_REF][START_REF] Haddar | Conformal mappings and inverse boundary value problems[END_REF][START_REF] Haddar | Conformal mapping and an inverse impedance boundary value problem[END_REF][START_REF] Haddar | Conformal mapping and impedance tomography[END_REF][START_REF] Kress | Inverse problems and conformal mapping[END_REF][START_REF]A conformal mapping method in inverse obstacle scattering[END_REF]. The second class of methods covers non iterative methods, generally based on the construction of an indicator function of the inclusion(s). These sampling/probe methods do not need to solve the forward problem, but require the knowledge of the full DtN map. With no claim as to completeness, let us mention the enclosure and probe method of Ikehata [START_REF] Ikehata | Reconstruction of the shape of the inclusion by boundary measurements[END_REF][START_REF] Ikehata | Reconstruction of the support function for inclusion from boundary measurements[END_REF][START_REF]On reconstruction in the inverse conductivity problem with one measurement[END_REF][START_REF] Ikehata | Numerical method for finding the convex hull of an inclusion in conductivity from boundary measurements[END_REF][START_REF] Erhard | A numerical study of the probe method[END_REF], the linear sampling method [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF][START_REF] Colton | A simple method using Morozov's discrepancy principle for solving inverse scattering problems[END_REF][START_REF] Cakoni | A qualitative approach to inverse scattering theory[END_REF]], Kirsch's Factorization method [START_REF] Brühl | Numerical implementation of two noniterative methods for locating inclusions by impedance tomography[END_REF][START_REF] Hanke | Recent progress in electrical impedance tomography[END_REF][START_REF] Kirsch | The factorization method for a class of inverse elliptic problems[END_REF] and Generalized Polya-Szegö Tensors [START_REF] Ammari | Reconstruction of small inhomogeneities from boundary measurements[END_REF][START_REF]Generalized polarization tensors, inverse conductivity problems, and dilute composite materials: a review[END_REF][START_REF]Polarization and moment tensors[END_REF][START_REF] Ammari | Generalized polarization tensors for shape description[END_REF][START_REF] Kang | Construction of conformal mappings by generalized polarization tensors[END_REF]. The reconstruction method proposed in this paper is non iterative and can be decomposed into two main steps. First, we show that the knowledge of DtN map gives access to the so-called Generalized Pólya-Szegö Tensors (GPST) of the cavity. This is done (see §. 3.1) by adapting to the multiply connected case the boundary integral approach proposed in [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF] for a simply connected cavity. The second step is to transform this shape from GPST problem into a (non standard) shape from moments problems (see §. 3.2). Reconstructing the geometry of the cavities amounts then to reconstructing the support of a density from the knowledge of its harmonic moments. Our reconstruction algorithm is then obtained by seeking a finite atomic representation of the unknown measure. Let us emphasize that we have been able to justify the connection between the shape from GPST problem and the shape from moments problem only in some particular cases (for a single cavity, for two disks and for small cavities). However, the reconstruction method turns out to be numerically efficient for arbitrary cavities. The paper is organized as follows. We collect some technical material from potential theory in Section 2. The reconstruction method is described in 3. Section 4 is devoted to the proof of Theorem 3.6. Finally, examples of numerical reconstructions are given in Section 5.

Background on potential theory

This section aims to revisit the results from potential theory given in [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF]Section 2.1.] in the context of a multiply connected cavity. For the proofs, we refer the reader to [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF] and to the books of McLean [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF], Steinbach [START_REF] Steinbach | Numerical approximation methods for elliptic boundary value problems[END_REF] or Hsiao and Wendland [START_REF] Hsiao | Boundary integral equations[END_REF] for more classical material.

G(x) = - 1 2π log |x|
the fundamental solution of the operator -∆ in R 2 . We pay careful attention to state the results in a form that includes multiply connected boundaries.

Single layer potential

We define the function spaces

H 1 2 (γ) := H 1 2 (γ 1 ) × H 1 2 (γ 2 ) × • • • × H 1 2 (γ N ), H -1 2 (γ) := H -1 2 (γ 1 ) × H -1 2 (γ 2 ) × • • • × H -1 2 (γ N ),
which are Hilbert spaces when respectively endowed with the norms

q 1 2 ,γ = q 1 2 1 2 ,γ 1 + • • • + q N 2 1 2 ,γ N 1 2 , ∀ q = (q 1 , . . . , q N ) ∈ H 1 2 (γ), q -1 2 ,γ = q1 2 -1 2 ,γ 1 + • • • + qN 2 -1 2 ,γ N 1 2 , ∀ q = (q 1 , . . . , qN ) ∈ H -1 2 (γ).
Definition 2.1. For every q = (q 1 , . . . , qN ) ∈ H -1 2 (γ), we denote by S γ q the single layer potential associated with the density q.

Given q = (q 1 , . . . , qN ) ∈ H -1 2 (γ), it is well-known that the single layer potential S γ q defines a harmonic function in R 2 \ γ. Furthermore, if qk ∈ L 2 (γ k ) for all k ∈ {1, . . . , N }, we can write:

S γ q(x) = γ G(x -y)q(y) dσ y = N k=1 S γ k qk (x), ∀ x ∈ R 2 \ γ.
The single layer potential defines a bounded linear operator from H -1 2 (γ) into H 1 oc (R 2 ), and the asymptotic behavior of S γ q reads as follows (see for instance [43, p. 261])

S γ q(x) = - 1 2π q, 1 -1 2 , 1 2 ,γ log |x| + O(|x| -1 ), (2.1)
where we have set for every q ∈ H -1 2 (γ) and p ∈ H

1 2 (γ): q, p -1 2 , 1 2 ,γ := N k=1 qk , p k -1 2 , 1 2 ,γ k , in which •, • -1 2 , 1
2 ,γ k stands for the duality brackets between H -1 2 (γ k ) and H 1 2 (γ k ). This shows in particular that S γ q has finite Dirichlet energy (i.e. ∇(S γ q) ∈ (L 2 (R 2 ) 2 ) for all q = (q 1 , . . . , qN ) in the function space

H(γ) := H(γ 1 ) × H(γ 2 ) × . . . × H(γ N ),
where for every k ∈ {1, . . . , N }:

H(γ k ) := {q k ∈ H -1 2 (γ k ) : qk , 1 -1 2 , 1 2 ,γ k = 0}. Remark 2.2. It is worth noticing that condition q ∈ H(γ) is only a sufficient condition ensuring ∇(S γ q) ∈ (L 2 (R 2 )) 2
, since a necessary and sufficient condition is: q, 1 -1 2 , 1 2 ,γ = 0. However, considering conditions (1.2d), H(γ) is clearly the appropriate function space to tackle the cavity problem (1.2).

We also recall that the single layer potential satisfies the following classical jump conditions

[S γ q ] γ k = 0, [∂ n (S γ q)] γ k = -q k , ∀ q = (q 1 , . . . , qN ) ∈ H -1 2 (γ). (2.2)
In the above relations, we have used the notation [u]

γ k = u i | γ k -u e | γ k and [∂ n u] γ k = (∂ n u i )| γ k - (∂ n u e )| γ k
, where u i and u e denote respectively the restrictions of a given function u to the interior and exterior of γ k . Let us emphasize that these classical formulae (and other trace formulae detailed below) may be impacted by different conventions concerning the definition of the fundamental solution, the unit normal or the jumps. Let us focus now on the trace of the single layer potential.

Definition 2.3. For every q = (q 1 , . . . , qN ) ∈ H -1 2 (γ), we denote by q = S γ q ∈ H 1 2 (γ) the trace on γ of the single layer S γ q: q = (q 1 , . . . , q N ),

q i = N j=1 S ij qj , where S ij qj = Tr γ i (S γ j qj ), ∀ i, j ∈ {1, . . . , N }. Note that S γ = (S ij ) 1 i,j N : H -1 2 (γ) → H 1 2 (
γ) defines a bounded linear operator with weakly singular kernel. Hence, if q = (q 1 , . . . , qN )

∈ H -1 2 (γ) is such that qk ∈ L ∞ (γ) for all k ∈ {1, . . . , N }, then S γ q(x) = γ G(x -y)q(y) dσ y , ∀ x ∈ γ.
Using Green's formula and the asymptotics (2.1), we can easily prove the identity

q, S γ q -1 2 , 1 2 ,γ = R 2 |∇(S γ q)| 2 dx < +∞, ∀ q ∈ H(γ). (2.3)
According to [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF]Theorem 8.12], S γ :

H -1 2 (γ) → H 1 2 (
γ) defines a strictly positive-definite operator on H (since H ⊂ {q ∈ H -1 2 (γ) : q, 1 -1 2 , 1 2 ,γ = 0}). It is also known (see [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF]Theorem 8.16]) that S γ is boundedly invertible if and only if the logarithmic capacity of γ (see [43, p. 264] for the definition) satisfies Cap(γ) = 1. From now on, and without loss of generality, let us assume that the diameter of Ω is less than 1 (otherwise, it suffices to rescale the problem), which implies in particular that Cap(Γ) < 1 and Cap(γ) < 1 (see [50, p. 143] and references therein).

In order to characterize the image of Ĥ(γ) by S γ , we need to introduce the following densities.

Definition 2.4. For every k ∈ {1, . . . , N }, we define the unique density:

êk := (ê k 1 , . . . , êk N ) ∈ H -1 2 (γ)
such that S γ êk is constant on γ for every ∈ {1, . . . , N } and satisfying the circulation conditions:

êk , 1 -1 2 , 1 2 ,γ = δ k, , = 1, . . . , N
where δ k, denotes Kronecker's symbol.

The existence and uniqueness of such functions êk is ensured by Lemma A.1 of the Appendix (simply take f = 0 and b as the k-th element of the canonical basis of R N ). Furthermore, the family {ê k , k = 1, . . . , N } is obviously linearly independent in H -1 2 (γ) and thus, so is the family

{e k := Sê k , k = 1, . . . , N } in H 1 2 (γ).
Proposition 2.5. The operator S γ defines an isomorphism from Ĥ(γ) onto

H(γ) := {q ∈ H 1 2 (γ) : êk , q -1 2 , 1 2 ,γ = 0, k = 1, . . . , N }.
Proof. We only need to prove that H(γ) = S γ ( Ĥ(γ)). Let q ∈ H -1 2 (γ) and set q := S γ q. We note that for all i ∈ {1, . . . , N }:

êi , q -1 2 , 1 2 ,γ = q, e i -1 2 , 1 2 ,γ = N j=1 qj , e i j -1 2 , 1 2 ,γ j = N j=1 e i j qj , 1 -1 2 , 1 2 ,γ j ,
where the last equality follows from the fact that e i j is constant on each boundary γ j . The matrix (e i j ) 1 i,j N being invertible, we have q ∈ Ĥ(γ) if and only if q ∈ H(γ).

The above result allows us to use the linear operator:

S γ : q ∈ H(γ) -→ q ∈ H(γ),
to identify any density q ∈ H(γ) with the trace q := S γ q ∈ H(γ).

Throughout the paper, we will systematically use this identification, using the notation with (respectively without) a hat on single layer densities of H(γ) (respectively traces of single layer densities).

Definition 2.6. For all p, q ∈ H(γ), we set:

p, q -1 2 ,γ = p, q 1 2 ,γ = p, q -1 2 , 1
2 ,γ . Obviously, using these inner products, the isomorphism S γ turns out to be an isometry between the spaces H(γ) and H(γ):

q 2 -1 2 ,γ = q 2 1 2 ,γ = R 2 |∇(S γ q)| 2 dx, ∀ q ∈ H(γ).
The following orthogonal projections will be needed in the next sequel.

Definition 2.7. Let Π γ and Π γ denote respectively the orthogonal projections from

H 1 2 (γ) into H(γ) and from H -1 2 (γ) into H(γ).
It can be easily checked that:

∀ q ∈ H -1 2 (γ) : Π γ q = q - N k=1 qk , 1 -1 2 , 1 2 ,γ k êk ∀ q ∈ H 1 2 (γ) : (Π γ q) k = q k -êk , q -1 2 , 1 2 ,γ 1 |γ k , k = 1, . . . , N.
Definition 2.8. We denote by Tr γ the classical trace operator (valued into H 1 2 (γ)), and by Tr 0 γ when it is left-composed with the orthogonal projection onto H(γ): Tr 0 γ := Π γ Tr γ . Let us recall a useful characterization of the norm chosen on H(γ) (the proofs of the assertions stated below are given in [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF]Section 2.1.] for the case of a simply connected cavity, but they can be easily extended to the multiply connected case). We define the quotient weighted Sobolev space:

W 1 0 (R 2 ) = {u ∈ D (R 2 ) : ρu ∈ L 2 (R 2 ), ∇u ∈ (L 2 (R 2 )) 2 }/R
, where the weight is given by

ρ(x) := 1 + |x| 2 log(2 + |x| 2 ) -1 , x ∈ R 2 ,
and where the quotient means that functions of W 1 0 (R 2 ) are defined up to an additive constant. This space is a Hilbert space once equipped with the inner product:

u, v W 1 0 (R 2 ) := R 2 ∇u • ∇v dx.
For q ∈ H(γ), and according to (2.3), we have

S γ q ∈ W 1 0 (R 2 ) and q -1 2 ,γ = q 1 2 ,γ = S γ q W 1 0 (R 2 ) , ∀ q ∈ H(γ).
To conclude this subsection, let us recall that the normal derivative of the single layer potential is not continuous across γ, as (with obvious matrix notation, and where the signs + andrefer respectively to the trace taken from the exterior and the interior of γ):

(∂ n (S γ q)) |γ ± = ± 1 2 + L * γ q, where L * γ = L * ij 1 i,j N : H -1 2 (γ) → H -1 2 (γ) (2.4)
where for smooth densities

(L * ij qj )(x) = γ j ∂ nx G(x -y)q j (y) dσ y , x ∈ γ i . (2.5)

Double layer potential

For more details about the results recalled in this section, we refer the interested reader to the monographs by Hackbusch [START_REF] Hackbusch | Integral equations[END_REF]Chapter 7], Kress [START_REF] Kress | Linear integral equations[END_REF]Chapter 7] or Wen [START_REF] Wen | Conformal Mappings and Boundary Value Problems[END_REF]Chapter 4].

Definition 2.9. For every q = (q 1 , . . . , q N ) ∈ H 1 2 (γ), we denote by D γ q the double layer potential associated with the trace q.

Given q = (q 1 , . . . , q N ) ∈ H 1 2 (γ), it is well-known that the double layer potential

D γ q(x) = γ ∂ ny G(x -y)q(y) dσ y = N k=1 D γ k q k (x), ∀ x ∈ R 2 \ γ
defines a harmonic function in R 2 \ γ whose normal derivative across γ is continuous, but whose trace is not continuous. More precisely, we have

Tr γ ± (D γ q) = ∓ 1 2 + L γ q, where L γ = (L ij ) 1 i,j N : H 1 2 (γ) → H 1 2 (γ)
, the adjoint of the operator defined in (2.4)-(2.5), is given by

(L ij q j )(x) = γ j ∂ ny G(x -y)q j (y) dσ y , x ∈ γ i .
One has in particular

[D γ q ] γ k = q. ( 2.6)
Later (see the proof of Theorem 4.6), we will need to compute the inner product of two densities p, q ∈ H using double layer potentials. This is provided by the two next lemmas, which deal respectively with the cases of real and complex valued densities.

Lemma 2.10. Given p, q ∈ H(γ), let u := S γ p and v := S γ q denote the single layers respectively associated with p, q ∈ H(γ). Let u be the function defined in R 2 \ γ by

u =    u -in ω u + in ω c := R 2 \ ω.
where u -(respectively u + ) denotes a harmonic conjugate of u in ω (respectively in ω c ). Similarly, we define the function v associated to v. Then, the functions u and v admit double layer representation formulae associated to two densities p, q ∈ H(γ):

u := D γ p, v := D γ q.
Moreover, we have

p, q 1 2 ,γ = -∂ n u, q -1 2 , 1 2 ,γ .
(2.7)

Proof. The harmonic conjugate u -of u in ω (i.e. the harmonic function such that u + i u -is holomorphic in ω) exists and is uniquely defined up to a constant (and this does not affect (2.7)). The existence of the harmonic conjugate u + of u in ω c is ensured by the fact that u = S γ p is circulation free on γ (since p ∈ H(γ)). Moreover, thanks to Cauchy-Riemann's equations we have

∇u = (∇ u) ⊥ in ω ∪ ω c . (2.8)
Since u admits by assumption a single layer representation, it is continuous across γ, and hence, the same holds for its tangential derivative. According to (2.8), this yields the continuity of ∂ n u across γ (note that trace of u is not continuous across γ). From classical integral representation formula for harmonic functions, this implies that u admits a double layer representation: u := D γ p, with p = [ u] γ . Defining similarly v := D γ q, we obtain using (2.8) that

p, q 1 2 ,γ = R 2 ∇u • ∇v = ω ∇ u • ∇ v + ω c ∇ u • ∇ v = -∂ n u, [ v] -1 2 , 1 2 ,γ ,
where the last equality follows from Green's formula. Equation (2.7) can then be deduced from the fact that

[ v] = [D γ q] = q.
In the rest of the paper, we still denote by •, • 1 2 ,γ the hermitian inner product on H 1 2 (γ) seen as a complex Hilbert space. For instance, for p = p 1 + ip 2 and q = q 1 + iq 2 we will have

p, q 1 2 ,γ = p 1 , q 1 1 2 ,γ + p 2 , q 2 1 2 ,γ + i p 2 , q 1 1 2 ,γ -p 1 , q 2 1 2 ,γ .
(2.9) Lemma 2.10 admits then the following counterpart for complex-valued densities.

Lemma 2.11. Given p = p 1 + ip 2 and q = q 1 + iq 2 with p j , q j ∈ H (j = 1, 2), let us set

u := S γ p = u 1 + iu 2 , v := S γ q = v 1 + iv 2
where u j = S γ pj and v j := S γ qj . Define on R 2 \ γ the complex-valued functions u = u 1 + i u 2 and v = v 1 + i v 2 where u j and v j (j = 1, 2) are respectively the harmonic conjugates of u j and v j (as defined in Lemma 2.10). Then, the functions u and v admit double layer representation formulae associated to two complex-valued densities p, q ∈ H(γ):

u := D γ p, v := D γ q.
Moreover, we have

p, q 1 2 ,γ = -∂ n u, q -1 2 , 1 2 ,γ , (2.10)
where

∂ n u = ∂ n u 1 + i∂ n u 2 .
Proof. From (2.7) and (2.9), we immediately get that

p, q 1 2 ,γ = -∂ n u 1 , q1 -1 2 , 1 2 ,γ -∂ n u 2 , q2 -1 2 , 1 2 ,γ -i ∂ n u 2 , q1 -1 2 , 1 2 ,γ -∂ n u 1 , q2 -1 2 , 1 2 ,γ
which is the claimed result.

Let us conclude by recalling the relation (in dimension two) between the double layer potential and the (complex) Cauchy transform. In the sequel, we identify a point x = (x 1 , x 2 ) of the plane with the complex number z = x 1 + ix 2 . The Cauchy transform of a density q defined on γ is given by:

C γ q(z) := 1 2iπ γ q(ζ) ζ -z dζ, z ∈ C \ γ.
Remark 2.12. It is worth noticing that the Cauchy transform can be easily computed in the following particular cases using Cauchy integral formula and Cauchy integral theorem.

1. If q is the trace of a holomorphic function F on ω, then

C γ q(z) =    F (z) z ∈ ω, 0 z ∈ ω c := C \ ω. 2. If q(z) = (z -a) -1 , with a / ∈ ω, then C γ q(z) =    (z -a) -1 z ∈ ω, 0 z ∈ ω c := C \ ω. 3. If q(z) = (z -a) -1 , with a ∈ ω, then C γ q(z) =    0 z ∈ ω, -(z -a) -1 z ∈ ω c := C \ ω.
It turns out (see for instance [25, p. 254] or [39, p. 100]) that for real-valued traces q, the double layer potential coincides with the real part of the Cauchy integral. Therefore, for complex-valued densities q, the Cauchy transform and the double layer potential are related via following formulae:

D γ (Re q) = 1 2 Re (C γ (q + q)) (2.11a) D γ (Im q) = 1 2
Im (C γ (qq)) .

(2.11b)

The above relations can be summarized in the identity:

D γ q = 1 2 C γ q + 1 2 C γ q, ( 2.11c)
or equivalently:

D γ q = 1 2 C γ q + 1 2 C γ q. (2.11d)
3 The reconstruction method

From DtN measurements to GPST

The first step of the proposed reconstruction method is to to recover the GPST of the cavities from the DtN measurements.

Definition 3.1. Let K γ Γ and K Γ γ be the operators:

K γ Γ : q ∈ H(Γ) -→ Tr 0 γ (S Γ q) ∈ H(γ), K Γ γ : p ∈ H(γ) -→ Tr 0 Γ (S γ p) ∈ H(Γ),
where Tr 0 γ and Tr 0 Γ are given in Definition 2.8. We define the boundary interaction operator between Γ and γ as the operator

K := K Γ γ K γ Γ .
One of the main results of [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF] (see Theorem 3.1) is to provide a relation between the measurements and the boundary interaction operator K := K Γ γ K γ Γ in the case of a single cavity. The preliminary results of Section 2 immediately lead to a generalization of this relation to the case of multiple cavities (see Theorem 3.3). For the proofs, we refer the interested reader to [46, Section 2.2.].

Proposition 3.2. The boundary interaction operators K γ

Γ and K Γ γ enjoy the following properties:

1. If p ∈ H(γ), then q := Tr Γ (S γ p) belongs to H(Γ).

Operators K γ

Γ and K Γ γ are compact, one-to-one and dense-range operators. Moreover, for every functions q ∈ H(Γ) and p ∈ H(γ), we have:

K γ Γ q, p 1 2 ,γ = q, K Γ γ p 1 2 ,Γ . (3.1)
3. The norms of the operators K γ Γ and K Γ γ are strictly less that 1.

Note that the first assertion in Proposition 3.2 shows that Tr 0 Γ can be replaced by Tr Γ in the definition of K Γ γ . Going back to the DtN operator Λ γ of problem (1.2), and according to (1.2d), we have by Green's formula

∂ n u f , 1 -1 2 , 1 2 ,Γ = -∂ n u f , 1 -1 2 , 1
2 ,γ = 0, which shows that Λ γ is valued in H(Γ). Considering data f ∈ H(Γ), we can thus define the DtN operator Λ γ as follows:

Λ γ : f ∈ H(Γ) -→ ∂ n u f ∈ H(Γ). (3.2)
Let us denote by Λ 0 the DtN map Λ γ in the case where ω = ∅ (the cavity free problem).

Theorem 3.3. The two following bounded linear operators in H(Γ):

R := S Γ (Λ γ -Λ 0 ) and K := K Γ γ K γ Γ ,
satisfy the following equivalent identities:

R = (Id -K) -1 K, K = (Id + R) -1 R. (3.3)
The above result shows that the knowledge of the DtN maps Λ γ and Λ 0 respectively corresponding to the cases with and without the cavities, entirely determines the boundary interaction operator K. Using (3.1), it is worth reformulating the second identity in (3.3) in a variational form:

K γ Γ f, K γ Γ g 1 2 ,γ = (Id + R) -1 Rf, g 1 2 ,Γ , ∀ f, g ∈ H(Γ). (3.4)
This identity can be used to compute the entries of the so-called polarization tensors of the multiply connected cavity. To make this statement precise, let us introduce the following definition.

Definition 3.4. Identifying x = (x 1 , x 2 ) in R 2 with the complex number z = x 1 + ix 2 , we define for every m 1, the harmonic polynomials of degree m:

P m 1 (x) = Re (z m ) and P m 2 (x) = Im (z m ) .
We define as well

Q m 1,Γ := Tr 0 Γ (P m 1 ) and Q m 2,Γ := Tr 0 Γ (P m 2 ) (3.5)
where the projected trace operator Tr 0 Γ is defined in Definition 2.8 (recall that this projector simply amounts to adding a suitable constant to the considered function). Finally, we set

Q m Γ := Q m 1,Γ + iQ m 2,Γ . (3.6)
The crucial point about these polynomials Q m ,Γ , = 1, 2, lies in the fact that since they are traces of harmonic functions, we have

K γ Γ (Q m ,Γ ) = Q m ,γ
, and hence, applying formula (3.4) with Q m ,Γ , we obtain that for all m, m 1 and all , = 1, 2:

K γ Γ Q m ,Γ , K γ Γ Q m ,Γ 1 2 ,γ = Q m ,γ , Q m ,γ 1
2 ,γ . These quantities are strongly connected with the so-called Generalized Pólya-Szegö Tensors (GPST) used in [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF] to reconstruct a single cavity. Unfortunately, the reconstruction method proposed there rests on the Riemann mapping theorem, which does not apply in the multiply connected case considered in this paper. We propose in the next section a new reconstruction method to recover the geometry of the cavities from the available data, namely the quantities

Q m ,γ , Q m ,γ 1 2 ,γ . Hence we can access to the complex quantities Q m γ , Q m γ 1
2 ,γ (see (2.9)). In fact, as already pointed out in [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF], the terms

Q m γ , Q 1 γ 1
2 ,γ contain enough information to recover the unknown geometry.

Towards a harmonic moments problem

The second step of the reconstruction is to recast the shape-from-GPST problem (namely reconstructing γ from the quantities

Q m γ , Q 1 γ 1
2 ,γ ) as a moments problem. Let us recall that the classical moments problem consists in recovering an unknown measure with support in K ⊂ C from its moments. The literature on this problem is very rich and covers a wide range of questions (solvability, uniqueness and reconstruction) and settings (dimension one or higher dimensions, arbitrary measures or measures absolutely continuous with respect to the Lebesgue measure, full or partial (harmonic or truncated) set of moments,...). Proposing a complete review is thus clearly beyond the scope of this paper. Let us simply make a few comments and quote some references. As far as we know, there is no framework to tackle this problem in full generality and most contributions address particular issues. The solvability of the moment problem in the cases K = [0, +∞[, K = {|z| = 1}, K = R and K = [0, 1] has been answered by the classical theorems of Stieltjes, Toeplitz, Hamburger and Hausdorff (see for instance Curto and Fialkow [START_REF] Curto | Recursiveness, positivity, and truncated moment problems[END_REF] and references therein). In dimension two, Davis studied the reconstruction of a triangle from four moments [START_REF] Davis | Plane regions determined by complex moments[END_REF], Milanfar et al. [START_REF] Milanfar | Reconstructing polygons from moments with connections to array processing[END_REF] investigated the reconstruction of arbitrary polygons and Putinar [START_REF] Putinar | A two-dimensional moment problem[END_REF] provided solvability conditions for the moments problem in the complex plane and pointed out some nice connections with quadrature domains (two-dimensional domains which are uniquely determined by finitely many of their moments, see [START_REF]Quadrature domains and their applications[END_REF]). Finally, two-dimensional shape reconstruction algorithms have been also proposed [START_REF] Golub | A stable numerical method for inverting shape from moments[END_REF][START_REF] Gustafsson | Reconstructing planar domains from their moments[END_REF][START_REF] Milanfar | Shape reconstruction from moments: theory, algorithms, and applications[END_REF][START_REF] Gustafsson | Bergman polynomials on an archipelago: estimates, zeros and shape reconstruction[END_REF]. For the problem studied in this paper, the connection between the shape-from-GPST problem and the shape-from-moments problem is based on the existence of a Borel measure ν supported in the cavity ω such that (z denotes the variable in the complex plane)

1 m Q m γ , Q 1 γ 1 2 ,γ = ω z m-1 dν, ∀ m 1. (3.7)
We have been able to prove this formula only in two special cases: for a simply connected cavity and in the case of two disks. Still, the numerical results (see Section 5) obtained for arbitrary cavities using the reconstruction method based on this formula (see subsection 3.3) are conclusive. This is why we state the following conjecture.

Conjecture 3.5. Let ω = ∪ N k=1 ω k be a multiply connected set, where the sets ω k , for k = 1, . . . , N are non intersecting simply connected open domains with C 1,1 boundaries γ k and ω ⊂ Ω. We set γ = ∪ N k=1 γ k . Then, there exists a Borel measure ν supported in ω such that:

f, Q 1 γ 1 2 ,γ = ω F dν, for every trace f ∈ H 1 2 (γ) of a holomorphic function F in ω.
As mentioned above, we have obtained the following result, whose proof is given in Section 4. In particular, the expression of the measure ν is given therein. Theorem 3.6. Conjecture 3.5 is true for N = 1 (single cavity) and when ω is constituted of two non intersecting disks.

The main interest of identity (3.7) lies in the fact that it bridges two classical inverse problems, namely the Calderòn's conductivity inverse problem and the historical moments problem. However, our moment problem has two features that makes it difficult to solve. First, we only have at our disposal the harmonic moments (i.e. the sequence ω z m dν, and not the doubly indexed sequence ω z m zn dν). Second, the involved measure is generally not of the form dν = 1 ω dx which is the most studied in the literature.

In order to have some insight on what this measure might represent, let us consider the particular case where the cavity is constituted of a collection of small inclusions. There is a wide literature dealing with this case and we refer the interested reader to the book by Ammari and Kang [START_REF] Ammari | Reconstruction of small inhomogeneities from boundary measurements[END_REF]. For small disks of centers z i and radii ερ i , i = 1, . . . , N (ε > 0 being a small parameter), and for the boundary conditions considered in this work, the DtN admits the following asymptotic expansion (apply for instance [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF]Theorem 2.1], in the case of non moving disks):

Λ ε = Λ 0 + ε 2 Λ 2 + O(ε 3 ), where Λ 2 f Γ , g Γ -1 2 , 1 2 ,Γ = 2π N i=1 ρ 2 i F (z i )G (z i ). (3.8)
In the above formula, f Γ , g Γ ∈ H (3.8), the above formula reads

that R ε = S Γ (Λ ε -Λ 0 ) = ε 2 S Γ Λ 2 + O(ε 3 ) and thus K ε = (Id + R ε ) -1 R ε = ε 2 S Γ Λ 2 + O(ε 3 ). Hence, we have f γε , Q 1 γε 1 2 ,γε = K ε f Γ , Q 1 Γ 1 2 ,Γ = ε 2 Λ 2 f Γ , Q 1 Γ -1 2 , 1 2 ,Γ + O(ε 3 ) Thanks to
f γε , Q 1 γε 1 2 ,γε = 2πε 2 N i=1 ρ 2 i F (z i ) + O(ε 3 ) = ω F dν + O(ε 3 ), provided we set ν = n i=1 c i δ z i , with c i = 2π(ερ i ) 2 .
This formula suggests that the measure ν can be approximated by an atomic measure, involving the centers and the radii of the unknown disks.

The reconstruction algorithm

Going back to the general case of an arbitrary cavity, it is thus natural to seek an approximation of the measure ν appearing in Conjecture 3.5 in the form ν * = n i=1 c i δ z i for some integer n, where the weights c i are positive and the points z i are distinct. To do so, we equalize the first 2n complex moments:

τ m = n i=1 c i z m i , ∀ m = 0, • • • , 2n -1 (3.9)
where

τ m := ω z m dν, ∀ m = 0, • • • , 2n -1.
The non linear system (3.9) with 2n unknowns is usually referred to as Prony's system. To solve it, we follow the method proposed by Golub et al. [START_REF] Golub | A stable numerical method for inverting shape from moments[END_REF]. For all integer 0 m 2n -1, let:

V =        1 1 • • • 1 z 1 z 2 • • • z n . . . . . . . . . . . . z n-1 1 z n-1 2 • • • z n-1 n        , H 0 =        τ 0 τ 1 • • • τ n-1 τ 1 τ 2 • • • τ n . . . . . . . . . . . . τ n-1 τ n • • • τ 2n-2        , H 1 =        τ 1 τ 2 • • • τ n τ 2 τ 3 • • • τ n+1 . . . . . . . . . . . . τ n τ n+1 • • • τ 2n-1       
.

Setting c = (c i ) 1 i n and z = (z i ) 1 i n , we easily see that if the pair (z, c) solves System (3.9), then the points z i , 1 i n, are the eigenvalues of the generalized eigenvalue problem

H 1 x = λ H 0 x. (3.10)
Once the (z i ) 1 i n have been determined, the weights (c i ) 1 i n can be easily computed (by solving the linear system (3.9) or computing the diagonal elements of the matrix

V -1 H 0 V -T ).
Considering the particular case of small disks described in subsection 3.2, the cavity will be reconstructed by drawing the disks of centers (z i ) 1 i n and radii (ρ i ) 1 i n , with ρ i = |c i |/2π (see Section 5 for some examples).

Proof of Theorem 3.6 4.1 Connected obstacle of arbitrary shape

If the cavity has only one connected component, then according to the Riemann mapping theorem, its boundary γ can be parametrized as follows t ∈] -π, π] → φ(e it ), where

φ : z → a 1 z + a 0 + m -1 a m z m ,
maps the exterior of the unit disk D onto the exterior of ω (see the book of Pommerenke [47, p. 5]). Let us emphasize that |a 1 | is the logarithmic capacity of γ and can be chosen such that a 1 > 0, while a 0 is the conformal center of ω. According to [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF]Lemma 3.5] (see formula (3.11a)), we have

Q m γ , Q 1 γ 1 2 ,γ = 2a 1 π -π e -it φ m (e it ) dt, m 1.
This formula can be easily generalized for every trace f ∈ H 1 2 (γ) of a holomorphic function F in ω as follows:

f, Q 1 γ 1 2 ,γ = 2a 1 π -π e -it F (φ(e it )) dt. ( 4.1)
Theorem 4.1. Assume that ω is a simply connected domain with C 1,1 boundary γ. Then, there exists a C ∞ diffeomorphism η that maps ω onto the unit disk D such that:

f, Q 1 γ 1 2 ,γ = 4a 2 1 ω F (ξ)|Dη| |dξ|.
In particular, Conjecture 3.5 is true for N = 1 with dν = 4a 2 1 |Dη| |dξ|.

Remark 4.2. Note that the diffeomorphism η is not a conformal mapping (not even a harmonic function) and satisfies η(ξ) = φ -1 (ξ) on γ = ∂ω.

Proof of Theorem 4.1. Integrating by parts in (4.1), we get that:

f, Q 1 γ 1 2 ,γ = 2a 1 π -π φ (e it )F (φ(e it )) dt = 2a 1 π -π φ (e it ) |φ (e it )| F (φ(e it ))
|φ (e it )|dt.

Denoting ψ = φ -1 that conformally maps the exterior of ω onto the exterior of the unit disk D, we obtain:

f, Q 1 γ 1 2 ,γ = 2a 1 γ |ψ (ξ)| ψ (ξ) F (ξ)d|ξ| = 2a 1 γ ψ(ξ)F (ξ)ψ(ξ) |ψ (ξ)| ψ (ξ) d|ξ|,
where we have used the fact that |ψ(ξ)| = 1. In the right hand side, we recognize:

ψ(ξ) |ψ (ξ)| ψ (ξ) = n 1 (ξ) + in 2 (ξ),
where n = (n 1 , n 2 ) is the unit normal to γ. We deduce that:

Re F (ξ)ψ(ξ) |ψ (ξ)| ψ (ξ) = ∂ n Re (F (ξ)) Im F (ξ)ψ(ξ) |ψ (ξ)| ψ (ξ) = ∂ n Im (F (ξ)). Hence 1 2a 1 f, Q 1 γ 1 2 ,γ = γ Re (ψ(ξ))∂ n Re (F (ξ)) + Im (ψ(ξ))∂ n Im (F (ξ))d|ξ| + i γ Re (ψ(ξ))∂ n Im (F (ξ)) -Im (ψ(ξ))∂ n Re (F (ξ))d|ξ|.
At this point, we would like to integrate by parts but ψ is not defined inside ω. We introduce the harmonic (but non holomorphic) diffeomorphism:

δ(z, z) = a 1 z + a 0 + a -1 z + a -2 z2 + . . . defined from D onto ω and such that δ(z, z) = φ(z) for z ∈ ∂D.
We define as well its inverse η = δ -1 (which in not even harmonic but only C ∞ ) and we have:

η(ξ, ξ) = ψ(ξ) for ξ ∈ γ.
Tolerating a slight abuse of notation, we shall write η(x, y) = η 1 (x, y)+iη 2 (x, y) = η(ξ, ξ), for ξ = x+iy.

We can now integrate by parts to get:

1 2a 1 f, Q 1 γ 1 2 ,γ = ω ∇η 1 (x, y) • ∇Re (F (ξ)) + ∇η 2 (x, y) • ∇Im (F (ξ))|dξ| + i ω ∇η 1 (x, y) • ∇Im (F (ξ)) -∇η 2 (x, y) • ∇Re (F (ξ))|dξ|.
Recall that ∇Re (F (ξ)) = (Re (F (ξ)), -Im (F (ξ))) and ∇Im (F (ξ)) = (Im (F (ξ)), Re (F (ξ))). We deduce that:

1 2a 1 f, Q 1 γ 1 2 ,γ = ω Re (F (ξ)) ∂η 1 ∂x + ∂η 2 ∂y + Im (F (ξ)) - ∂η 1 ∂y + ∂η 2 ∂x |dξ| + i ω Re (F (ξ)) ∂η 1 ∂y - ∂η 2 ∂x + Im (F (ξ)) ∂η 1 ∂x + ∂η 2 ∂y |dξ|,
which can be simply rewritten as:

1 2a 1 f, Q 1 γ 1 2 ,γ = ω F (ξ) ∂η 1 ∂x + ∂η 2 ∂y + i ∂η 1 ∂y - ∂η 2 ∂x |dξ|.
On the one hand, we have:

∂η ∂x = ∂η 1 ∂x + i ∂η 2 ∂x ∂η ∂y = ∂η 1 ∂y + i ∂η 2 ∂y .
Hence, we deduce that (see [18, p. 3])

∂η 1 ∂x + ∂η 2 ∂y = Re ∂ ∂x -i ∂ ∂y η = 2Re ∂η ∂ξ ∂η 1 ∂y - ∂η 2 ∂x = Re ∂ ∂y + i ∂ ∂x η = -Im ∂ ∂x -i ∂ ∂y η = -2Im ∂η ∂ξ ,
and therefore:

1 4a 1 f, Q 1 γ 1 2 ,γ = ω F (ξ) ∂η ∂ξ |dξ|. (4.2)
On the other hand, using the chain rule formulae (see [18, p. 5]), we have:

∂ ∂z (η • δ) = ∂η ∂ξ ∂δ ∂z + ∂η ∂ ξ ∂δ ∂ z = 1 ∂ ∂ z (η • δ) = ∂η ∂ξ ∂δ ∂ z + ∂η ∂ ξ ∂δ ∂z = 0.
Eliminating ∂η ∂ ξ from the above relations, we get that

∂δ ∂z 2 - ∂δ ∂ z 2 ∂η ∂ξ = ∂δ ∂z . ( 4.3)
But, on the one hand, a direct calculation shows that ∂δ ∂z = a 1 , while on the other hand, we know (see [18, p. 5]) that for a harmonic diffeomorphism, we have:

∂δ ∂z 2 - ∂δ ∂ z 2 = |Dδ| = |Dη| -1 .
Substituting these relations in (4.3), we obtain that

∂η ∂ξ = a 1 |Dη|.
The conclusion follows then from (4.2).

The case of two disks

We assume here that ω = ω 1 ∪ ω 2 , where ω 1 and ω 2 are two disks centered at z 1 and z 2 respectively and with radii ρ 1 and ρ 2 . We set ρ := |z 1 -z 2 | and we denote by ω c j = C \ ω j (j = 1, 2). Proposition 4.3. There exists two sequences of points (z 1,n ) n 0 and (z 2,n ) n 0 and four sequences of complex numbers (κ 1,n ) n 1 , (κ 2,n ) n 1 , (α 1,n ) n 0 and (α 2,n ) n 0 such that the following properties hold true:

1. z 1,n ∈ ω 1 and z 2,n ∈ ω 2 for all n 0.

The densities defined by

q j,0 (z) = 2(z -z j ),
and for all n 1:

q 1,n (z) = κ 1,n z -z 2,n-1 + α 1,n , q 2,n (z) = κ 2,n z -z 1,n-1 + α 2,n ,
are such that for j = 1, 2 the series n 0 q j,2n + q j,2n+1 converges in H 1 2 (γ j ) to some limit q j .

3. The limits q 1 ∈ H 1 2 (γ 1 ) and q 2 ∈ H 1 2 (γ 2 ) are such that

D γ 1 q 1 (ξ) + D γ 2 q 2 (ξ) =    ξ -z 1 in ω 1 ξ -z 2 in ω 2 .
The proof of this result is constructive: the sequences (z j,n ) n 0 , (κ j,n ) n 1 and (α j,n ) n 0 (j = 1, 2) are respectively given by formulae (4.4), (4.5) and (4.9). In order to prove the above Proposition, let us introduce some further notation and prove two lemmas. Setting:

Λ j = ρ j ρ 2 j = 1, 2
let us define the two real sequences (λ 1,n ) n 0 and (λ 2,n ) n 0 by:

λ 1,0 = λ 2,0 = 1,
and for n 1

λ 1,n = 1 - Λ 2 λ 2,n-1 λ 2,n = 1 - Λ 1 λ 1,n-1 ,
It is worth noticing that

λ 1,n = 1 - Λ 2 1 - Λ 1 λ 1,n-2 .
Lemma 4.4. The real sequences (λ 1,n ) n 0 and (λ 2,n ) n 0 are convergent. More precisely, we have:

lim n→+∞ λ 1,n = 1 2 1 + Λ 1 -Λ 2 + √ Λ and lim n→+∞ λ 2,n = 1 2 1 + Λ 2 -Λ 1 + √ Λ ,
where

Λ := |Λ 1 -Λ 2 | 2 + 1 -2(Λ 1 + Λ 2 ) > 0.
Moreover, for every n 0:

1 -Λ 2 < 1 2 1 + Λ 1 -Λ 2 + √ Λ < λ 1,n 1, 1 -Λ 1 < 1 2 1 + Λ 2 -Λ 1 + √ Λ < λ 2,n 1.
Proof. The convergence of the two sequences (λ j,n ) n 0 (j = 1, 2) depends on the sign of Λ. Noticing that Λ 2 /Λ 1 = (ρ 2 /ρ 1 ) 2 and rewriting Λ as:

Λ = Λ 2 2 ρ 2 1 ρ 2 2 -1 2 -2Λ 2 ρ 2 1 ρ 2 2 + 1 + 1,
we deduce that:

Λ = 1 - (ρ 1 -ρ 2 ) 2 ρ 2 1 - (ρ 1 + ρ 2 ) 2 ρ 2 ,
and therefore that Λ > 0.

Since (1 -

Λ 1 + Λ 2 ) > 0 and Λ = (1 -Λ 1 + Λ 2 ) 2 -4Λ 2 , we have: √ Λ -(1 -Λ 1 + Λ 2 ) < 0.
It follows that:

(1 -Λ 1 + Λ 2 ) 2 -4Λ 2 + Λ -2(1 -Λ 1 + Λ 2 ) √ Λ = 2 √ Λ( √ Λ -(1 -Λ 1 + Λ 2 )) < 0.
We deduce first that:

4Λ 2 > (1 -Λ 1 + Λ 2 ) 2 + Λ -2(1 -Λ 1 + Λ 2 ) √ Λ = (1 -Λ 1 + Λ 2 - √ Λ) 2 ,
and hence that

1 -Λ 2 < 1 2 1 + Λ 1 -Λ 2 + √ Λ .
Let us define the two complex sequences:

z 1,n = (1 -λ 2,n )z 2 + λ 2,n z 1 , z 2,n = (1 -λ 1,n )z 1 + λ 1,n z 2 . (4.4)
We also introduce the complex sequences (κ 1,n ) n 1 defined by:

κ 1,1 = 2ρ 2 2 κ 1,2 = - 2(ρ 1 ρ 2 ) 2 (z 2 -z 1 ) 2 , κ 1,n+2 = Λ 1 Λ 2 (λ 2,n λ 1,n-1 ) 2 κ 1,n , n 1. (4.5)
Note that we have the following relations for all n 1:

κ 1,n+2 = (1 -λ 1,n+1 )(1 -λ 2,n ) λ 1,n-1 λ 2,n κ 1,n = Λ 1 Λ 2 (λ 1,n-1 -Λ 1 ) 2 κ 1,n . (4.6)
We define similarly a sequence (κ 2,n ) n by exchanging the role of the boundary indices. The sequences (κ 1,n ) n and (κ 2,n ) n can be deduced from each other according to the following identities:

κ 1,n+1 = - ρ 2 2 (z 2 -z 1 ) 2 λ 2 2,n-1 κ 2,n = - (z 2 -z 1 ) 2 |z 2 -z 1 | 2 (1 -λ 1,n ) λ 2,n-1 κ 2,n (4.7
)

κ 2,n+1 = - ρ 2 1 (z 1 -z 2 ) 2 λ 2 1,n-1 κ 1,n = - (z 1 -z 2 ) 2 |z 2 -z 1 | 2 (1 -λ 2,n ) λ 1,n-1 κ 1,n . (4.8)
Lemma 4.5. There exists a constant C > 0 and a constant 0 < δ < 1 such that:

|κ 1,n | + |κ 2,n | Cδ n ∀n 1.
Proof. According to Lemma 4.4 and to relation (4.6), it suffices to prove that:

0 < 4Λ 1 Λ 2 1 -Λ 1 -Λ 2 + √ Λ 2 < 1.
Direct computations lead to:

(1 -Λ 1 -Λ 2 + √ Λ) 2 -4Λ 1 Λ 2 = 2 √ Λ √ Λ + (1 -Λ 1 -Λ 2 ) ,
and the quantity in the right hand side is clearly positive. The conclusion follows.

We are now in position to prove the claims detailed in Proposition 4.3.

Proof of Proposition 4.3.

1. According to the estimate of Lemma 4.4, we have:

|z 2,n -z 2 | = ρ(1 -λ 1,n ) < ρ Λ 2 = ρ 2 ,
and therefore, the point z 2,n ∈ ω 2 for every n 0. Of course, similar arguments show that z 1,n ∈ ω 1 . 2. Let the densities q 1,n be defined as in Proposition 4.3, with

α 1,n = - 1 2 κ 1,n λ 1,n-1 (z 1 -z 2 )
. (4.9)

According to lemmas 4.4 and 4.5, the series of functions:

n 0 κ 1,2n z -z 2,2n-1 + α 1,2n + κ 1,2n+1 z -z 2,2n + α 1,2n+1 ,
converges in H 1 (ω 1 ) and hence its trace on γ 1 converges in H 1 2 (γ 1 ) to some limit q 1 . The same reasoning holds for j = 2. 3. According to Remark 2.12, the densities q j,0 and qj,0 (j = 1, 2) admit the following Cauchy transforms:

C γ j q j,0 (ξ) =    2(ξ -z j ) in ω j 0 in ω c j and C γ j qj,0 (ξ) =        0 in ω j - 2ρ 2 j ξ -z j in ω c j ,
where, for the last equation, we have used the relation:

ξ -z j = ρ 2 j ξ -z j on γ j .
Based on formulae (2.11), we deduce that:

D γ j q j,0 (ξ) =        ξ -z j in ω j - ρ 2 j ξ -z j in ω c j .
In particular, we have

D γ 1 q 1,0 (ξ) + D γ 2 q 2,0 (ξ) =          (ξ -z 1 ) - ρ 2 2 ξ -z 2 in ω 1 (ξ -z 2 ) - ρ 2 1 ξ -z 1 in ω 2 .
Similarly, and using once again Remark 2.12, we have for every n 1:

C γ 1 q 1,n (ξ) =    κ 1,n ξ -z 2,n-1 + α 1,n in ω 1 0 in ω c 1 .
(4.10)

In order to compute C γ 1 q1,n (ξ), we note that for every ξ ∈ γ 1 , we have:

1 ξ -z 2,n-1 = 1 ρ 2 1 ξ-z 1 + z 1 -z 2,n-1 (4.11) = 1 z 1 -z 2,n-1   1 - ρ 2 1 z 1 -z 2,n-1 1 ξ -z 1 - ρ 2 1 z 1 -z 2,n-1    .
(4.12)

But notice now that, on the one hand:

z 1 -z 2,n-1 = λ 1,n-1 (z 1 -z 2 ),
and on the other hand:

z 1 - ρ 2 1 z 1 -z 2,n-1 = z 1 - ρ 2 1 λ 1,n-1 (z 1 -z 2 ) = z 1 - ρ 2 1 ρ 2 λ 1,n-1 (z 1 -z 2 ) = 1 - Λ 1 λ 1,n-1 z 1 + Λ 1 λ 1,n-1 z 2 = λ 2,n z 1 + (1 -λ 2,n )z 2 = z 1,n .
Plugging the last two relations in (4.12) and using the identity

ρ 2 1 (λ 1,n-1 ) 2 (z 1 -z 2 ) 2 κ 1,n = -κ 2,n+1 , we obtain that for ξ ∈ γ 1 : κ 1,n ξ -z 2,n-1 = -2α 1,n + κ 2,n+1 ξ -z 1,n . Hence q1,n (ξ) = κ 2,n+1 ξ -z 1,n -α 1,n .
We deduce first that:

C γ 1 q 1,n (ξ) =        -α 1,n in ω 1 - κ 2,n+1 ξ -z 1,n in ω c 1 ,
and next, using (4.10), that:

D γ 1 q 1,n (ξ) = 1 2 C γ 1 q 1,n (ξ) + 1 2 C γ 1 q 1,n (ξ) =          1 2 κ 1,n ξ -z 2,n-1 in ω 1 - 1 2 κ 2,n+1 ξ -z 1,n in ω c 1 .
After similar computations for D γ 1 q 2,n , we finally obtain that

D γ 1 q 1,n (ξ) + D γ 2 q 2,n (ξ) =          1 2 κ 1,n ξ -z 2,n-1 - 1 2 κ 1,n+1 ξ -z 2,n in ω 1 1 2 κ 2,n ξ -z 1,n-1 - 1 2 κ 2,n+1 ξ -z 1,n in ω 2 .
It is then easy to verify that:

n k=0 j=1,2 D γ j q j,2k (ξ) + D γ j q j,2k+1 (ξ) =          (ξ -z 1 ) - 1 2 κ 1,2n+2 ξ -z 2,2n+1 in ω 1 (ξ -z 2 ) - 1 2 κ 2,2n+2 ξ -z 1,2n+1 in ω 2 .
Letting n go to infinity and invoking again Lemma 4.5, we obtain the claimed result.

We can now prove the main results of this section.

Theorem 4.6. There exists two sequences of complex weights (c 1,n ) n 0 and (c 2,n ) n 0 such that for every

f = (f 1 , f 2 ) ∈ H 1 2 (γ 1 ∪ γ 2
) with f j the trace of a holomorphic function F j in ω j , there holds

f, z 1 2 ,γ = j=1,2 n 0 c j,n F j (z j,n ).
In other words, we have

f, z 1 2 ,γ = ω F dν,
where the measure ν is defined by: ν = j=1,2 n 0 c j,n δ z j,n .

Proof. Let us first prove the following identity:

f, z 1 2 ,γ = j=1,2 -i γ j q j (z)F j (z) dz. ( 4

.13)

We use here the notation of Lemma 2.11. If

F j = F 1 j + iF 2 j , let F j = F 1 j + i F 2 j where F j ( = 1, 2) is the harmonic conjugate of F j . Obviously, since F j is harmonic in ω j , we have F 1 j = F 2 j and F 2 j = -F 1 j
in ω j , and hence F j = -iF j in ω j . Similarly, we have z -z j = -i(z -z j ) = -iD γ q = D γ (-iq) in ω j . Therefore, we have q = -iq. Formula (2.10) shows that

f, z 1 2 ,γ = -∂ n F , q -1 2 , 1 2 ,γ = -∂ n F, q -1 2 , 1
2 ,γ = -j=1,2 γ j q j (z)∂ n F j (z) ds z .

Since ∂ n F j (z) = ∂ n F 1 j (z) + i∂ n F 2 j (z) = F j (z)n(z) and dz = in(z)ds z , the above formula immediately yields (4.13). As q j = n 0 q j,2n + q j,2n+1 , equation (4.13) implies that f, z 1 2 ,γ = j=1,2 n 0 -i γ j q j,2n (z)F j (z)dz -i γ j q j,2n+1 (z)F j (z)dz.

The density q j,2n+1 being the trace of a holomorphic function in ω j , the last terms in the right hand side vanishes. Let us focus on the term:

-i

γ 1 q 1,2n (z)F 1 (z)dz = -i γ 1 κ 1,2n 1 z -z 2,2n-1 F 1 (z)dz.
But for z ∈ γ 1 , we have:

1 z -z 2,2n-1 = 1 λ 1,n-1 (z 1 -z 2 ) 1 - ρ 2 1 λ 1,2n-1 (z 1 -z 2 ) 1 z -z 1,2n ,
and therefore:

-i

γ 1 q 1,2n (z)F 1 (z)dz = iρ 2 1 κ 1,2n λ 2 1,2n-1 (z 1 -z 2 ) 2 γ 1 1 z -z 1,2n F 1 (z)dz = - 2πρ 2 1 κ 1,2n λ 2 1,2n-1 (z 1 -z 2 ) 2 F 1 (z 1,2n ) = 2πκ 2,2n+1 F 1 (z 1,2n ),
where the last equality follows from (4.8). Adding the contribution corresponding to γ 2 , we finally obtain that f, z 1 2 ,γ = 2π j=1,2 n 0 κ 2,2n+1 F 1 (z 1,2n ) + κ 1,2n+1 F 2 (z 2,2n ) , and the proof is complete.

Numerical tests

We collect in this section some numerical experiments illustrating the feasibility of the proposed reconstruction method. For the sake of clarity, we first sum up the different steps of the simple reconstruction algorithm:

1. Compute a numerical approximation of the operator R = S Γ (Λ γ -Λ 0 ).

2. Fix an integer n 1 and compute for 0 m 2n -1:

τ m := ω z m dν = 1 m + 1 Q m+1 γ , Q 1 γ 1 2 ,γ = 1 m + 1 Q m+1 Γ , (Id + R) -1 RQ 1 Γ 1 2 ,Γ .
3. Following the method described in §. 3.3, solve Prony's system:

n i=1 c i z m i = τ m , ∀ m = 0, • • • , 2n -1.
to determine the positive weights (c i ) 1 i n and (distinct) points (z i ) 1 i n .

4. Plot the disks of centers (z i ) 1 i n and radii (ρ i ) 1 i n , with ρ i = |c i |/2π.

We refer the interested reader to [START_REF] Munnier | Conformal mapping for cavity inverse problem: an explicit reconstruction formula[END_REF]Section 4] for technical details about the implementation of steps 1 and 2 of the algorithm. For practical reconstructions, a natural question is how to determine the number n of atoms (disks) to be used. From our numerical experiments, there is no clear answer to this issue. However, increasing n generally yields reconstructions of better quality. This fact is illustrated in Figures 2,3 and 4 which show respectively examples of reconstructions for a rectangular cavity, a clover shaped cavity and a multiply connected cavity (with three connected components).

The following remarks are worth being mentioned:

• Increasing the number of atoms from n to n + 1 does not result in just adding an additional disk. Indeed, this leads to a new Prony's system and hence, to a completely new distribution of disks.

• For a given value of n, we can obtain atoms with zero radii which seem to be randomly distributed outside the cavity. See for instance Figure 2(d), where among the 22 atoms, only 13 have non zero radii.

• Sometimes, spurious atoms can be observed (see Figure 5

Figure 1 :

 1 Figure 1: The multiply connected cavity.

1 2 (

 2 Γ) denote respectively the traces on Γ of holomorphic functions F, G defined in the cavity ω ε . With the notation of Theorem 3.3, this expansion shows

  (a)), but they disappear when n increases (see Figure5(b)).

(a) n = 1 (b) n = 5 (

 15 c) n = 10 (d) n = 22

Figure 2 :

 2 Figure 2: Example of reconstruction of a rectangle cavity for different values of n, the number of atoms (the red dots represent the centers of the disks). Note that some disks are degenerate (i.e. have zero radii).

(a) n = 1 (b) n = 3 (

 13 c) n = 8 (d) n = 17

Figure 3 :

 3 Figure 3: Example of reconstruction of a clover shaped cavity for different values of n.

(a) n = 1 (b) n = 3 (

 13 c) n = 9 (d) n =

Figure 4 :

 4 Figure 4: Example of reconstruction of a multiply connected cavity for different values of n.

A Appendix

The next lemma generalizes to the case of a multiply connected boundary the result given in [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF]Lemma 8.14] for a simply connected boundary. Our proof is slightly different from the one given there, although it also uses the Fredholm alternative. The notation are those of Section 2 and we recall that the assumption Cap(γ) = 1 is supposed to hold true.

Then, there exists a unique density q = (q 1 , Ê . . . , qN ) ∈ H -1 2 (γ) and a unique c = (c 1 , . . . , c N ) T ∈ R N satisfying the system of equations system of equations

Proof. Let us introduce the operator A defined on

in which

.

Using this notation, system (A.1) simply reads

Clearly, A defines a bounded operator from

it is clear that A is a Fredholm operator (C is a finite rank operator and S γ is boundedly invertible). Hence, (A.1) is uniquely solvable if and only if the corresponding homogeneous problem admits only the null solution. Let then (q, c) ∈ H -1 2 (γ) × R N satisfying A q c = 0 and denote by u := S γ q the single layer potential corresponding corresponding to q. Then, according to the first equation in (A.1) (with f = 0), u is constant on γ and thus constant in ω. On the other hand, the second equation in (A.1) implies that q ∈ H(γ) and hence u vanishes outside ω. Consequently, qk = -[∂ n u] γ k = 0 and c = -S γ q = 0.


Ping Han

Haoliang Gan

Weikun He

Aircraft Attitude Estimation Based on Central Difference Kalman Filter

Keywords: Attitude estimation, Quaternion, CDKF, UKF, EKF, SPKF

When the extended Kalman filter (EKF) is applied in the aircraft attitude estimation, two defects exist: one is compu tational complexity; the other is large linearization error. Aiming at these problems, central difference Kalman filter (CDKF) based on Stirling interpolation formulation is applied to the low-cost aircraft attitude estimation system which is of less accurate and high noisy sensors. First, the nonlinear mathematical model of aircraft attitude based on quaternion is established, then CDKF is applied to attitude estimation. Experimental results with real flying data show that CDKF is superior to the commonly used EKF method and unscented Kalman filter (UKF). The algorithm not only improves the attitude estimation precision and stability effectively, but also avoids the computing burden of Jacobian matrices. In addition, it is more simple and easy to implement, because it has only one adjustable parameter instead of three in the UKF circumstances.

I. INTRODUCTION

Aircraft attitude can be defined as the orientation relation ship between the body frame and the earth frame, usually described by a set of Euler angles such as yaw, pitch and roll. It provides not only the pilot with the aircraft navigation information, but also a three-dimensional attitude reference for the autopilot, fire control system, the radar antenna and aerial camera and other airborne equipments. Therefore, it is a very important parameter for safety flying and in the above applications.

Aircraft attitude is usually measured by the strapdown inertial navigation system [l] . First, the nonlinear math ematical model of aircraft attitude based on quaternion is established, then the appropriate filter is utilized to estimate the attitude. For the nonlinear aircraft attitude estimation problem, extended Kalman filter (EKF) is probably the most widely used at present. However, it is only reliable for systems which are almost linear and the Jacobian matrices are hard to obtain. In fact, these difficulties arise from its use of linearization [START_REF] Schaub | Analytical Mechanics of'Space Systems[END_REF] . To address this issue, unscented Kalman filter (UKF) was proposed by Julier and Uhlmann [START_REF] Julier | Unscented filtering and nonlinear estimation[END_REF] . UKF can be applied directly to nonlinear systems, and it uses a num ber of deterministic sampling points to capture the posterior mean and covariance of the pertinent Gaussian approximate densities, which can be accurate up to the second order of any nonlinearity [START_REF] Yang | Gaussian Sum Particle Filter for Spacecraft Attitude Estimation[END_REF] . However, three scalar scaling parameters Daniel Alazard, Franc,;ois Defay, Yves Briere Institut Superieur de l' Aeronautique et de I'Espace, France (a, /3,�) are employed by the UKF. These parameters are always selected based on the nonlinearities of the system models [START_REF] Merwe | Sigma-point Kalman .filters /iJr probabilistic inference in dynamic state[END_REF], and the choice of the parameters' values will affect the filter's estimate precision.

In this paper, Stirling interpolation in the state estimator is employed for aircraft attitude estimation problems. Stirling interpolation replaces the unscented transform in the UKF algorithm architecture. As proved in [START_REF] Merwe | Sigma-point Kalman .filters /iJr probabilistic inference in dynamic state[END_REF], [START_REF] Ito | Gaussian Filters for Nonlinear Filtering Prob lems[END_REF], Stirling inter polation based central difference Kalman filter (CDKF) has the same or superior performance as the unscented transform based UKF, with one advantage over UKF: The CDKF uses only a single scalar scaling parameter, the central difference interval size h , as opposed to the three (a, /3, �) that the UKF uses [START_REF] Julier | Unscented filtering and nonlinear estimation[END_REF]. As shown in our experiment, the CDKF algorithm can get better perfonnance on the estimate accuracy and error robustness.

The paper is organized as follows. The attitude estimation algorithm based on the CDKF is given in Section II. Section III gives the relationship between quatemion and kinematics, the gyroscope model and establishes the nonlinear system model of aircraft attitude. Section IV presents the experimental results of the EKF, UKF and CDKF, and their estimate performances are compared. Finally, Section V summarizes this paper.

II. CENTRAL DIFFERENCE KALMAN FILTER

The CDKF is the core member of the sigma-point Kalman filters (SPKF) family of algorithms based on the sigma-point approach [START_REF] Merwe | Sigma-point Kalman .filters /iJr probabilistic inference in dynamic state[END_REF]. It uses a symmetric set of sigma points derived by Stirling interpolation to approximate nonlinear functions, whereas the EKF uses the Taylor series. The advantage of the Stirling interpolation is that the calculation of Jacobian matrices is not required, and can be accurate up to the second order of any nonlinearity. In case of Gaussian distributions of the system variables, the mean and covariance can be represented by these sigma points. More details about Stirling polynomial interpolation method for approximating nonlinear models can be found in [START_REF] Norgaard | RavnNew Developments in State EstimationfiJr Nonlinear Systems[END_REF].

Consider the nonlinear dynamic system composed of state equation and observation equation as follows.

(1) where X k ER n is the state vector, Y k E R7n is the observation vector, W k E R n is the state noise vector, and V k E R7n is the measurement noise vector. It is assumed that W k and V k are zero-mean Gaussian noise processes with covariances given by Q k and R k o respectively. The flow of the central difference Kalman filter is formulated as follows.

1) Initialization:

X o = E[xoJ, P o = E[ ( xo -xo) ( xo -xo fl (2)
2) Calculate sigmap oints for time-u p date:

Xk-l = [ Xk-l Xk-l + hVPk-l Xk-l -hVPk-ll (3)
where h is the interval length, if the random variables obey a Gaussian distribution, the optimal value of h is V3.

3) Time-u p date equations:

These sigma points are further passed through the nonlinear function F(.) , such that the predicted sigma points for the discrete time k are derived [START_REF] Yang | Gaussian Sum Particle Filter for Spacecraft Attitude Estimation[END_REF] Finally, the first two moments of the predicted state vector are obtained by linear regression of the transformed sigma points n 2n ""

" rn * x k/k-l = L Wi Xi,k/k-l i=O Pk/k-1 = 2)W�1(x7,k/k-l -X�+i,k/k-l ) (x7,k/k-l - i=l (5) * ) T + C2 ( * + * 2 * ) Xn+i,k/k-l w i Xi,k/k-l Xn+i,k/k-l -XO,k/k-l
P XkYk = J Wf' Pk/k-1 (Y�.n,k/k-l -Y�+12n,k/k-lf (11)
Hence, the update of the estimated means and estimated error covariance at time k are given by: ( 12)

Pk = Pk/k-1 -KkPy,K[ ( 13
)
where Kk is the Kalman gain and can be defined as:

(14) ) Cl where W o -�,w i -W Z -1, ... ,2n ,w o 1 C') h 2-1' 4 h ' Wi-= 4 h 4 ( z = l, ... ,n). The rate sensor, accelerometer and magnetometer measure angular rates, the gravity vector and local magnetic field vector respectively. The accelerometer and magnetometer are used for the rate gyro drift correction in order to improve the precision and stability of attitude estimation.

III. AIRCRAFT ATTITUDE NONLINEAR MODEL X (X7,k/k-l + X�+i,k/k-l -2 X � ,k/k-lf l + Q k-l (
( 8 ) A. Quaternion and Kinematics Furthermore, the mean, covariance and cross-covariance are derived by n 2n

" rn * Yk/k-l = L Wi Yi,k/k-l i=O PYk = L [ W�l (y7,k/k-l -Y�+i,k/k-l ) (Y7,k/k-l - i=l (9) * ) T + C') ( * + * 2 * ) Yn+i,k/k-l Wi-Yi,k/k-l Yn+i,k/k-l -YO,k/k-l x (ytk/k-l + Y�+i,k/k-l -2 y � ,k/k_l ) T l + R k (10)
Quaternion [START_REF] Shuster | A Survey of Attitude Representation[END_REF] is widely used to describe the attitude of an aircraft, because it provides nonsingular attitude descriptions and expresses arbitrary and large rotations of the aircraft [START_REF] John | Survey of Nonlinear Attitude Estimation Methods[END_REF]. The quaternion is a four-dimensional vector defined as follows qu = [a b c dV. Note that the quatemion has the constraint The quaternion kinematic differential equation [START_REF] John | Unscented Filtering for Spacecraft Attitude Estimation[END_REF] can be derived as follows.

(16) where -q P o w = [p q r] T ,and p,q and r are aircraft angular rate around x axis, y axis and z axis respectively in body frame.

To solve the equations of kinematic, we can get quaternion, then make use of the relationship between quatemion and attitude angle, we can obtain the yaw, pitch and roll angles as follows

{ I.(J = atan2(2( ab + cd), 1 -2(b2 + c2)) e = arcsin(2( ac -bd)) 1j; = atan2(2(ad + bc), 1 -2( c2 + d2))

C. State Equations

We choose the quaternion, angular rates and gyro biases as state vector, that is The relationship among state vector is shown in Fig. 2.

We can get the nonlinear state equation from Fig. 2 We can get nine measurements which can be chosen as a measurement vector directly from the rate sensor, accelerome ter and magnetometer. Therefore the measurement vector can be defined as [Gx Gy Gz l and [Mx My Mz l are the gravity vector and the magnetic field vector respectively whose values are known, b x 3 is a 3 order unit matrix , n is the measurement noise vector which is a zero-mean Gaussian white-noise process.

IV. EXPERIMENT RESULTS

In the experiment, the EKF algorithm which is internally integrated in the CHR-6 dm sensor and the SBG AHRS sensor combines data from onboard accelerometer, rate gyros, andmagnetic sensors to produce yaw, pitch, and roll angle estimates, the untreated angular rates, gravity vector and magnetic field vector data. The performance, reliability and ease of use of the SBG AHRS sensor have been kept to the highest levels, but it is more expensive compared with CHR system. The price of CHR system is about 89 dollars [START_REF]CHR-6dm AHRS[END_REF], while SBG system is more than 2400 dollars. In order to verify CDKF algorithm's estimated effect, we utilize the CDKF method to treat the angular rates, gravity vector and magnetic field vector data obtained from the the CHR-6 dm sensor to get attitude esti mation compared with the output of the SBG AHRS sensor. The very well with the output of the SBG AHRS sensor, Since the statistics of nonlinearly transformed Gaussian approximate random variables in CDKF are calculated by the deterministic sampling based sigma-point approach without using the Jacobian matrices and the linearization of the nonlinear model, the uncertainties are propagated well and the accuracy of the state estimation has been improved over the EKF approach. The same advantages also exits in UKF, however, it has only one adjustable parameter in contrast to three in the UKF circumstances, which affect the accuracy of the state estimation. Hence, the performance of CDKF is better than that of UKE

V CONCLUSION

Aiming at the problem that low quality of the data collected from less accurate and high noisy sensors lead to the poor attitude estimation based on EKF, this paper presents an aircraft attitude estimation algorithm by introducing the CDKF method based on Stirling interpolation formulation, Experimental results show that the CDKF algorithm can get better performance on the estimate accuracy and error robustness than that of EKF and UKF. Moreover, it it does not employ the calculation of the Jacobian matrices and the linear approximations to the nonlinear models. So it provides a new way to improve the precision of attitude estimation.

Fig. l .

 l Fig. l. Experimental aircraft equipped with the SBG and CRR sensor

  6) The CHR-6 dm AHRS sensor produced by the CH Robotics company in USA and the SBG AHRS sensor made by the rn _ h 2-n rn _ 1 ( . _

4 )

 4 Calculate sigmap oints for measurement-u p date: 5) Measurement-u p date equations: The predicted measurement points are obtained by transform ing the sigma points through H ( • ) French SBG System company are installed in the model plane which is applied for experimental flight, as shown in Fig, L Each of them contains a three-axis magnetometer, a three-axis angular rate sensor, and a three-axis accelerometer.

Fig. 2 .

 2 Fig. 2. Process Model for Angular Rates and Quatemion

  where I.(J , e and 1j; are yaw,pitch and roll respectively.

  this study, a gyroscope is used as the attitude sensor of the aircraft. The gyroscope is a sensor that measures the angular rate of the aircraft. The gyroscope system can be expressed mathematically by modeling the measured angular velocity as the true angular velocity with an additive bias[11]. The bias dynamics are considered to be driven by a Gaussian white noise process. The gyroscope model can be represented as (1 8 ) where W 9 = [p g qg r g] T is the measured angular velocity, W = [p q r] T is the true angular velocity, b is the drift and b = [bp bq br] T , nb = [nbp nbq nbr] T is a zero-mean Gaussian white-noise process .

  as follows p = P g -bp q = qg -bq r = rg -br bp = nbp bq = nbq br = nbr . 1 qu = 20(w)qu D. Measurement Equation (20)

  y = [Accx Accy Accz Magx Magy Magz P g qg rg] (21)where Accx,Accy and Accz are the gravity vector in the body frame, and ]\;1 agx,lvi agy and ]\;1 agz are the local magnetic field vector in the body frame.According to the basic principle of strapdown inertial navi gation, the relationship between what is measured in the body frame and the known values in the earth frame is:

  ____ � ______ �� __ �� ____ � ____ �� ____ � ""'"

Fig, 3 ,Fig. 4 .Fig, 5 ,Fig. 6 .

 3456 Fig, 3, Yaw and Yaw Estimation Error

Pitch

Fig. 7 .Fig. 8 .

 78 Fig. 7. Partial enlarged drawing of pitch

ACKNOWLEDGMENT This work is supported by the National Natural Science Foundation of China and Civil Aviation Administration of China (Item No. 60979002) and in part by the Fund of Civil Aviation University of China under grant ZXH2012DOOI and 20l2KYE03.