Munoz Victor
email: victor.munozcuartas@enit.fr

Marianne Perrin
email: marianne.perrin@iut-tarbes.fr

Welemane Hélène
email: helene.welemane@enit.fr

Cantarel Arthur
email: arthur.cantarel@iut-tarbes.fr

Karama Moussa
email: moussa.karama@enit.fr

APPLICATION OF NON DESTRUCTIVE TESTING TO THE DETECTION OF AERONAUTICAL DEFECTS IN COMPOSITE STRUCTURES

Keywords: Non-destructive testing, Ultrasonics, Infrared thermography, Composites, Simulated defects

A study of two Non-destructive Testing methods (NDT) was carried out in specimens with different kinds of simulated defects. Ultrasonics test (US) and Infrared Thermography (IRT) were applied with the aim to evaluate the detectability and the accuracy of each method. These techniques have acquired great importance in the aeronautics industry because they allow to control the aerostructures without intervening in their physical and mechanical integrity. In the second part of the study, a comparison of both techniques was achieved in order to analyse their limits and advantages. It appeared that detectability of defects was much better in a sample with flat-bottomed holes defects in the case of Ultrasonic Test. However it was found that Infrared Thermography is much more limited to the thickness of the specimen than the ultrasonic waves. On the other hand, defects were all revealed with IRT in a sandwich composite including Teflon inserts, which was not the case for US.

INTRODUCTION

The use of composites in the aerospace industry has increased dramatically since the 1970s. The primary benefits that composite components can offer are the reduction of weight and the simplification of assembly. When they are in service, different types of mechanical and thermal loads are applied to these structures. They generate internal stresses. It is in this way that delamination or disbonding may result in the aerostructure. Defects such as random inclusions or undesirable material have also been found during the manufacture process of composite materials [START_REF] Garnier | The detection of aeronautical defects in situ on composite structures using Non Destructive Testing[END_REF]. These flaws may lead to stress concentration with serious consequences. It is therefore important to test the composite structure to ensure their integrity. These defects have to be revealed in order to make the correct maintenance or replacement. The faster the damage or defect is detected, the safer the aerostructure is.

Many techniques have been used in order to detect such defects in the operational configuration. Liquid penetrant inspection [START_REF] Cherfaoui | Essais non destructifs[END_REF], magnetic control inspection [START_REF] Cherfaoui | Essais non destructifs[END_REF], eddy currents control [START_REF] Cherfaoui | Essais non destructifs[END_REF], radiographic testing [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF] and shearography [START_REF] Cherfaoui | Essais non destructifs[END_REF] are some of the methods used to evaluate the material health, without altering its properties. These kinds of controls are called the Non-Destructive Testing (NDT). The most widespread technique is the ultrasonic testing. It is the only technique certified by the aeronautics industry [START_REF] Boro-Djordjevic | Ultrasonic characterization of advanced composite materials[END_REF]. Ultrasonic Test (US) uses high frequency waves to conduct examinations and make measurements. Ultrasonic inspection is used for flaw detection/evaluation, dimensional measurements and material characterization [START_REF] Budnik | Essais non destructifs -Contrôle par ultrasons[END_REF]. Some studies have already shown the detectability and the accuracy of a flaw in size, shape and depth for different kinds of composite structures used in the aeronautic industry [START_REF] Garnier | The detection of aeronautical defects in situ on composite structures using Non Destructive Testing[END_REF][START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF].

A full-field measurement technique based on Infrared Thermography (IRT) is also used in the NDT field. Infrared thermography is a nondestructive, non-intrusive, non-contact mapping of thermal patterns or thermograms, on the surface of objects. IRT is more widely used in recent years for structural investigation [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF]. The main difference about IRT with regard to US is that the former does not give the information about the depth of defects. However, detection is faster than US. In fact, the detection of defects with IRT is immediate; whilst that in US, analysis takes about 10 minutes. Some works aimed at comparing different infrared thermography configurations have detailed the accuracy of defects detection [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF][START_REF] Péronnet | Non destructive investigation of defects in composite structures by three infrared thermographic techniques[END_REF].

The aim of this work is to evaluate two different specimens by ultrasonic test and infrared thermography. Tests on a carbon-epoxy laminated composite and a carbon-epoxy-glass sandwich composite were carried out to check the accuracy and detectability of their simulated defects by using both NDT methods. Specimens have different kinds of defects and were made in different materials used in the aeronautic industry. Defects have various diameters and locations on the material so that the size limit and deep detection can be determined. The second objective is to compare the data provided by each technique and to find out the advantages of each one to show their possible complementary relative to each other. Finally, a conclusion about the intention to capitalize the information obtained by NDT methods is presented.

SPECIMENS AND EXPERIMENTAL SETUP 2.1 Specimens

The specimen used to compare the detectability and accuracy between the IRT and the US is a carbon/epoxy laminated composite with flat bottomed holes defects and variable thickness. It is called specimen A (figure 1). This sample is divided in 2 sections. The first one is 4 mm wide in 50 mm long. The second section is 130 mm long with 8 mm wide. This sample contains 15 holes with diameters ranging from 2 mm to 10 mm. These holes have different depths (from 0.52 mm to 7.5 mm). The aim of these holes is to simulate discontinuity on the specimen that is considered as a continuous medium. The second sample, called specimen B, was only tested with infrared thermography because of its material properties. This specimen is a carbon-glassepoxy skin with foam-core sandwich composite including Teflon inserts. Its skin is 0.66 mm thick. It has been demonstrated that the foam-core is a sound insulation; thereby the ultrasonic waves cannot penetrate this specimen and be detected by the receiving transducer [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF]. This specimen contains 18 defects. These simulated flaws are Teflon inserts ranging from 3 mm to 6 mm. They are located at the skin section. These inserts were placed during the stacking of the tissues so that they simulate as well as possible a delamination within the sample. The manufacturing drawing of this specimen is shown in figure 2.

Ultrasonics

In ultrasonics, the sound energy is introduced and propagates through the material in the shape of waves. When there is a discontinuity (such as inclusion, delamination or disbonding) in the wave path, part of the energy will be reflected back from the flaw surface. The reflected wave signal is turned into an electrical signal by the transducer and is displayed on a screen. The configuration used on this US testing is called reflection mode (figure 3). One transducer is used to send and receive the wave at the same time.

Figure 3. US configuration

This technique allows to get information about the material elasticity by measuring the wave velocity [START_REF] Baste | Evaluation de la matrice d'élasticité des composites orthotropes par propagation ultrasonore en dehors des plans principaux de symétrie[END_REF]. The depth, the size and the shape of a detected defect can also be obtained by this technique.

In this study, an Omniscan 32: 128 PR is used with a 5L64-NW1 multi-elements transducer connected to it. The transducer is both ultrasonic source and receiver at the same time. The frequency of it is 5 MHz, which is the common frequency in aeronautics NDT [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF]. For specimen A, the wave velocity is 2983 m s -1 with a standard deviation of 5.7 m s -1 . This experiment was carried out by contact; therefore it is necessary to use a coating gel between the sample and the transducer. Ultrasonic testing enables twodimensional mapping, accordingly a coding system is used to detect the position of the transducer. Specimen A was placed in a support to hold it while the experiment was being carried out. Figure 4 shows the ultrasonics device.

Infrared Thermography

Infrared thermography (IRT) provides maps and movies of the surface thermal field of specimens; in this way, local changes in surface temperature indicate subsurface defects. Thermal waves are propagated through the material. When they pass through a medium with different thermal conductivity, a thermal gradient is generated in the specimen. Indeed, the two mediums have different thermal conductivities. These differences are captured by an infrared camera enabling the heat emitted by the sample to be converted to temperature. It is measured on the front of the specimen. Thermal two-dimensional mapping is created and heterogeneities can be detected [START_REF] Garnier | The detection of aeronautical defects in situ on composite structures using Non Destructive Testing[END_REF][START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF][START_REF] Péronnet | Non destructive investigation of defects in composite structures by three infrared thermographic techniques[END_REF].

This technique allows to get information about the position, the size and the shape of the flaw. So, delamination, disbonding or crack networks can be detected [START_REF] Garnier | The detection of aeronautical defects in situ on composite structures using Non Destructive Testing[END_REF].

The heat source is a 1000 W halogen lamp. It is used in two configurations: transmission and reflection mode (figure 5). Specimen A and B were tested in both modes (figure 5a and5b). The infrared sensor was an IR camera with thermal resolution of 20 mK.

A 30 s and 50 Hz movie was recorded in which 15 s corresponds to the heating of the sample and the 15 last seconds to its cooling. These times are long enough to heat the sample and make defects to appear in the thermal mapping. Defects are thermally visible while the heating and cooling of the sample. The treatment of the movie was done with the ALTAIR software. The analysis of the movie is done in relative mode with the purpose to remove the ambient temperature and all environmental noise. Lock-in IRT was also used on both specimens. In lock-in IRT temperature modulation induced at the surface of the inspected component from the outside propagates as a harmonic thermal wave. Lock-in IRT allows a better defect inspection than common IRT and it is less sensitive to environmental conditions [START_REF] Choi | Quantitative determination of a subsurface defect of reference specimen by lock-in infrared thermography[END_REF]. The connection of the lock-in IRT is detailed on figure 6. A signal generator is also needed in addition to the classical IRT.

RESULTS

Ultrasonics

As said on section 2, only specimen A was tested with ultrasonics, because foam-core properties in specimen B prevent wave propagate.

Results analyses are made through two parameters: detectability and accuracy. 14 of the 15 holes (93%) were detected by US. The only defect that was not detected is located at 1 mm deep and 2 mm diameter. It is noticed that defects located near the contact surface (less than 1 mm) are hardly distinguished from the signal input peak. Nevertheless defects at 0.52 mm deep were visible in the C-Scan mapping (figure 7). On the other hand, the detected defects were measured on average 1 mm larger or smaller than they really are. Defects located at 7.5 mm deep (3 and 6 mm diameter) were measured 1 mm smaller than their real size, and those located at 5.2 mm deep (2, 3 and 6 mm diameter) were measured at the correct dimension. 10 mm diameter defects were all measured 1 mm larger than their real size. Table 1 shows the measured dimension for each range of defects. Figure 7 shows the C-Scan mapping for this specimen in which measurements were done.

The size of defects was measured with a visual rule.

In this case it is said that there is a defect where the colour of the health part changes and becomes a circular spot with another colour. Other works [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF] have obtained similar detectability with a better accuracy for the same specimen. The gap between the real size and the measured size was 0.5 mm. Those measurements were done by using the -6 dB method. It consists in measuring the defect where the backwall echo decreases of 50% [START_REF] Budnik | Essais non destructifs -Contrôle par ultrasons[END_REF].

Infrared thermography

Specimen A

As in US test, two parameters were taken into account to analyse the results: detectability and accuracy. A single picture must be selected from the recorded film. That selection depends on the quality of the image for each defect. In this study, for specimen A, 2 different images in transmission mode (figure 8) were selected to measure 2 ranges of defects located in two different sections of the specimen. These ranges of defects are located at 0.52 mm and 2.8 mm deep respectively and they are in the 4 mm width section and 8 mm width section of the sample. Figure 8 shows that only 6 from the 15 defects, were detected using classical IRT in transmission mode. In figure 8b it is noticed that heat saturation is generated after a short period of time. The 4 mm width section is rapidly saturated and in the 8 mm width section, no other defect is revealed with the increasing temperature.

The size of the measures defects has the same gap than in US testing. They were measured 1 mm larger than they really are. Table 2 shows the detailed size of the defects. 6 mm diameter defect at 2.8 mm deep was detected but no measurement could be done because the temperature difference between damaged and healthy parts was not enough to be measured. It is also noticed that the number of detected defects decreases with depth. While the defect is deeper, more difficult is the detection.

In lock-in IRT, the number of detected defects in specimen A increased. Best results were obtained in reflection mode and with 0.1 Hz frequency. However, no measurement of their size was done because of the poor temperature contrast. In addition to classical IRT (with transmission mode), 3 mm diameter defect at 2.8 mm deep and 5.2 mm deep defects (6 and 10 mm diameter) were visible as shown in figure 9.

Specimen B

In the case of specimen B, lock-in IRT leads to similar results as classical IRT. All defects were detected and they could be measured in diameter (figure 10). Better results were obtained at 0.1 Hz for both specimens tested in lock-in IRT. It was found that the textured surface of specimen B was confounded with the defects in IRT mapping. Measurements were therefore less accurate than classical IRT in specimen A.

Defects located at the camera side were easily distinguished whereas, defects located on the opposite side of the camera, were difficultly detected. That is way both skins of the sandwich were tested. In addition, the dimension of the measurements was 2 or 3 mm larger than the real size of the defect. Table 3 shows the measured size for each range of defects. Péronnet et al. [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF] tested the same sample and found a gap between real size and measured size ranging from 1.5 mm to 3.5 mm. This stands in accordance with the present study.

Nowadays there is no specific method for measuring the real size of a defect from an IRT mapping. The criterion clearly depends on the operator. In some works, the size of defects has been calculated using a system of grey image measurement [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite: détection de défauts par mesures de champs[END_REF].

ANALYSIS

A comparative analysis was done between both techniques. Table 4 summarizes the limits of detection for each technique regarding different parameters that were also measured while doing the tests.

The C-Scan mapping in US shows that detectability is more affected by the diameter of defects rather than thickness. This detection becomes more difficult when the defect is close to the detecting a b surface and to the backwall (less than 1 mm). In this case the defect echo is coincident with the echo of the detection surface.

On the other hand, it is clearly noticed that IRT is very limited to large thicknesses. Defects are not visible beyond 2.8 mm in a carbon-epoxy laminated composite because heat diffusion inside the material does not reach the inspected surface in transmission mode (figure 5a). However, by the use of lock-in IRT, the thickness reached comes to 5.2 mm in reflexion mode (figure 5b). At any rate, diameter of defects also affects the detectability in IRT. Figure 8b and 9 show that flaws located in the same range of deepness with different diameter are not all detected. The heat diffusivity inside the material does not lead to enough temperature difference in the surface of the sample. In the case of lock-in thermography, low frequencies are more convenient for detecting the flaws in both specimens. The use of a sinusoidal signal with low frequency (0.1 Hz) allows to optimize the wave spread in the specimen thickness. Other works [START_REF] Péronnet | Non destructive investigation of defects in composite structures by three infrared thermographic techniques[END_REF] with the same sample found that the best frequency for detectability of defects is 0.05 Hz.

For the carbon-glass-epoxy skin with foam-core sandwich composite, the use of infrared thermography is quite difficult since the foam is a thermal insulator. It is recommended to control the 2 skins of the sandwich separately and in reflection mode (figure 5b). Otherwise, the heat flux may not arrive to the detection surface.

It was found that detection time in IRT is immediate; however the analysis time is almost the same for both techniques due to the visual methods used to calculate the defects diameter. In US the detection time is longer because the whole surface has to be scanned by the transducer to get the C-Scan mapping.

CONCLUSION

The subject of this study was to compare the ultrasonics test and the infrared thermography as NDT techniques. By using 2 different samples, it was shown that both methods have the same accuracy regarding the diameter of the flaws, which is not the case for detectability parameter. 93% of defects in a carbon/epoxy laminated composite were detected by the US while only 60% of them were detected by the IRT.

In the case of a carbon-glass-epoxy skin with foamcore sandwich composite, all defects were detected by IRT. Results were similar using both, classical IRT and lock-in IRT for this specimen, but the measured size was less accurate than for the laminated composite. The textured surface was confused with the defects. It is also remarkable that no defect was detected by the use of US in this specimen due to absorption of ultrasonics waves by the foam-core.

Based on the results obtained with this study, it is noted that both techniques could be complementary. Firstly, IRT could give the information about the presence of defects and their position; then an US test would give the information about their size, their kind and also their depth.

In future works, mechanical characterization on carbon-epoxy laminated composites will be studied by calculating the elastic constants from the velocity of the ultrasonic waves. Then NDT techniques will be used for monitoring the damage evolution. Static and dynamic tests (tensile, fatigue) will be carried out and at the same time, NDT techniques (Ultrasonics, Infrared Thermography and Acoustic Emissions) will monitor the damage evolution.

Figure 1 .

 1 Figure 1. Manufacturing drawing of the specimen A

Figure 2 .

 2 Figure 2. Manufacturing drawing of the specimen B 2.2 Experimental setup

Figure 4 .

 4 Figure 4. Ultrasonic measurement

Figure 5 .

 5 Figure 5. Infrared measurement

Figure 6 .

 6 Figure 6. Lock-in thermography connection

Figure 7 .

 7 Figure 7. C-Scan mapping on specimen A

Figure 8 .

 8 Figure 8. IRT mapping of specimen A Defects at 0.52 mm deep were measured using image 8a and those at 2.8 mm deep using image 8b.In figure8bit is noticed that heat saturation is generated after a short period of time. The 4 mm width section is rapidly saturated and in the 8 mm width section, no other defect is revealed with the increasing temperature.

Figure 9 .

 9 Figure 9. Lock-in IRT in specimen A

Figure 10 .

 10 Figure 10. Lock-in IRT in specimen B

Table 1 .

 1 Defect measurement on specimen A using US

			Depth (mm)
		0.52	1	2.8 5.2 7.5
	Real diameter (mm)		Measured diameter (mm)
	10	11	-	11 11 11
	6	7	-	6	6	5
	3	4	-	4	3	2
	2	-	Undetected	-	2	2

Table 2 .

 2 Defect measurement on specimen A using IRT

			Depth (mm)
		0.52 1	2.8	5.2	7.5
	Real
			Measured diameter (mm)
	(mm)
	10	11 -	11	Undetected Undetected
	6	7 -No measurable Undetected Undetected
	3	4 -	Undetected	Undetected Undetected
	2	-2	-	Undetected Undetected

Table 3 .

 3 Defect measurement on specimen B using lock-in IRT

	Real diameter (mm) Measured diameter (mm)
	6	8
	3	5

Table 4 .

 4 Characteristics of US and IRT

	Parameter	US	IRT
	Accuracy in size	1 mm (larger)	1 to 2 mm (larger)
	Minimal diameter	2 mm from 1 mm deep	2 mm from 0 mm deep
	Detection in depth	More than 8 mm in laminated	2.8 mm for a 6 mm diameter defect in
		composites	laminated composites
	Detectability in the laminated	14/15. Non detected a 2 mm defect at	6/15. Non detected defects beyond 5.2
	composite	1 mm deep	mm deep and 3 mm diameter
	Detectability in foam-core sandwich	-	18/18 by analysing both sides of the
	composites		sample
	Time for getting results	About 10 min	Instantly
	General limitation	Detection of defects close to the	Important specimen thickness
		inspection surface and the backwall


Habib Ammari
email: habib.ammari@ens.fr

Elie Bretin
email: bretin@cmap.polytechnique.fr

Josselin Garnier
email: garnier@math.jussieu.fr

Abdul Wahab
email: wahab@cmap.polytechnique.fr

Time Reversal Algorithms in Viscoelastic Media *

Keywords: AMS subject classifications. 35L05, 35R30; Secondary 47A52, 65J20 Key words. Elastic wave propagation, Time reversal algorithms, Attenuation correction

In this paper, we consider the problem of reconstructing sources in a homogeneous viscoelastic medium from wavefield measurements using time-reversal algorithms. Our motivation is the recent advances on hybrid methods in biomedical imaging. We first present a modified time-reversal imaging algorithm based on a weighted Helmholtz decomposition and justify mathematically that it provides a better approximation than by simply time reversing the displacement field. Then, we investigate the source inverse problem in an elastic attenuating medium. We provide a regularized time-reversal imaging which corrects the attenuation effect at the first order.

Introduction

Waves in loss-less media are invariant under time transformation t → -t. This simple observation has provided very promising techniques in a variety of domains including biomedical imaging [START_REF] Fink | Acoustic time-reversal mirrors[END_REF], seismology [START_REF] Larmat | Timereversal imaging of seismic sources and application to the great Sumatra earthquake[END_REF], material analysis [START_REF] Prada | Time reversal techniques in ultrasonic nondestructive testing of scattering media[END_REF], land-mine detection [START_REF] Norville | Time-reversal focusing of elastic surface waves[END_REF], telecommunication [START_REF] Lerosey | Time-reversal of electromagnetic waves and telecommunication[END_REF] and underwater acoustics [START_REF] Fink | Time reversed acoustics[END_REF]. In time-reversal algorithms, an output wave is time reversed and retransmitted into the medium. By time invariance and reciprocity properties, the retransmitted wave retraces its path back through the medium and converges to the location of initial sources. This convergence of the reverted wave holds whatever the complexity of the underlying medium, which can be homogeneous or scattering, dispersive or non-dispersive. See for instance [START_REF] Ammari | An Introduction to Mathematics of Emerging Biomedical Imaging[END_REF][START_REF] De Rosny | Time reversal of electromagnetic waves[END_REF][START_REF] Fink | Time reversed acoustics[END_REF][START_REF] Fink | Acoustic time-reversal mirrors[END_REF][START_REF] Fouque | Wave Propagation and Time Reversal in Randomly Layered Media[END_REF][START_REF] Tanter | Time reversing waves for biomedical Applications[END_REF] and references therein for comprehensive details and discussions on time reversal. See also [START_REF] Ammari | Time reversal in attenuating acoustic media[END_REF][START_REF] Ammari | Transient elasticity imaging and time reversal[END_REF][START_REF] Catheline | Time-reversal of elastic waves in soft solids[END_REF][START_REF] Xu | Time reversal and its application to tomography with diffraction sources[END_REF] for applications of time-reversal techniques in biomedical imaging.

The robustness and simplicity of time-reversal techniques make them an ideal choice to resolve source localization problems. These inverse problems have been of significant interest in recent years and find numerous applications in different fields, particularly in biomedical imaging [START_REF] Ammari | Noise Source Localization in an Attenuating Medium[END_REF][START_REF] Ammari | Mathematical modelling in photoacoustic imaging of small absorbers[END_REF][START_REF] Ammari | Mathematical models and reconstruction methods in magneto-acoustic imaging[END_REF][START_REF] Ammari | Transient anomaly imaging by the acoustic radiation force[END_REF][START_REF] Phung | Time reversal focusing of the initial state for Kirchhoff plate[END_REF][START_REF] Teng | Some theoretical problems on variational data assimilation[END_REF]. If the sources are temporally localized, the source localization problems are actually equivalent to find initial states of a system governed by differential equations from the observations over some finite interval of time.

In this work, we consider the problem of reconstructing sources in a viscoelastic medium from wavefield measurements using time-reversal methods. Our motivation is the recent advances on hybrid methods in biomedical imaging exploiting elastic properties of the soft tissues [START_REF] Greenleaf | Selected methods for imaging elastic properties of biological tissues[END_REF]. Examples of these hybrid methods include magnetic resonance elastography [START_REF] Ammari | A method of biological tissues elasticity reconstruction using magnetic resonance elastography measurements[END_REF][START_REF] Catheline | Time-reversal of elastic waves in soft solids[END_REF], transient elasticity imaging [START_REF] Ammari | Transient elasticity imaging and time reversal[END_REF], shear wave imaging [START_REF] Sarvazyan | Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics[END_REF] and acoustic radiation force imaging [START_REF] Ammari | Transient anomaly imaging by the acoustic radiation force[END_REF][START_REF] Bercoff | The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force[END_REF]. The envisaged problem is quite challenging, indeed, because time reversibility of the wave equation breaks down in lossy media. Further, if not accounted for, these losses produce serious blurring in source reconstruction using classical time-reversal methods. In this paper, we use a thermo-viscous law model for the attenuation losses. We refer, for instance, to [START_REF] Kowar | Photoacoustic imaging taking into account attenuation[END_REF][START_REF] Kowar | Causality analysis of frequency dependent wave attenuation[END_REF][START_REF] Bretin | On the Green function in visco-elastic media obeying a frequency power-law[END_REF] for detailed discussions on the attenuation models in wave propagation and their causality properties.

The main contributions of this paper are twofold. We first provide a modified timereversal imaging algorithm in non-attenuating media based on a weighted Helmholtz decomposition and justify both analytically and numerically that it provides a better approximation than by simply time reversing the displacement field. Then, we give a regularized time-reversal imaging algorithm for source reconstruction in attenuated media and show that it leads to a first-order approximation in terms of the viscosity parameters of the source term. We present a variety of numerical illustrations to compare different time-reversal algorithms and to highlight the potential of our original approach.

Time reversal in homogeneous elastic media without viscosity

Let us consider the homogeneous isotropic elastic wave equation in a d-dimensional open medium:

     ∂ 2 u ∂t 2 (x, t) -L λ,µ u(x, t) = dδ 0 (t) dt F(x), (x, t) ∈ R d × R, u(x, t) = ∂u ∂t (x, t) = 0, x ∈ R d , t < 0, (2.1)
where L λ,µ u = µ∆u + (λ + µ)∇(∇ • u).

(2.2)

Here (λ, µ) are the Lamé coefficients of the medium and the density of the medium is equal to one. The aim in this work is to design efficient algorithms for reconstructing the compactly supported source function F from the recorded data g(y, t) = u(y, t), t ∈ [0, T ], y ∈ ∂Ω ,

where Ω is supposed to strictly contain the support of F. We are interested in the following time-reversal functional:

I(x) = ˆT 0 v s (x, T )ds, x ∈ Ω, (2.4)
where the vector field v s is defined as the solution of

     ∂ 2 v s ∂t 2 (x, t) -L λ,µ v s (x, t) = dδ s (t) dt g(x, T -s)δ ∂Ω (x), (x, t) ∈ R d × R, v s (x, t) = ∂v s ∂t (x, t) = 0, x ∈ R d , t < s.
(2.5)

Here, δ ∂Ω is the surface Dirac mass on ∂Ω and g = u on ∂Ω×R is the measured displacement field.

The time-reversal imaging functional I is usually implemented to reconstruct the source distribution in an elastic medium [START_REF] Catheline | Time-reversal of elastic waves in soft solids[END_REF][START_REF] Norville | Time-reversal focusing of elastic surface waves[END_REF][START_REF] Phung | Time reversal focusing of the initial state for Kirchhoff plate[END_REF]. It is motivated by the time reversibility property of the elastic waves. In a general setting, however, it is not sure that it provides a good reconstruction of the source distribution F. Indeed the problem is that the recorded displacement field at the surface of the domain is a mixture of pressure and shear wave components. By time reversing and backpropagating these signals as in (2.4) a blurred image is obtained due to the fact that the pressure and shear wave speeds are different.

In this work, we first present a modified time-reversal imaging functional I, and justify mathematically that it provides a better approximation than I of the source F. This new functional I can be seen as a correction based on a weighted Helmholtz decomposition to I (which is considered as an initial guess). In fact, we find the compressional and the shear components of I such that

I = ∇ × ψ I + ∇φ I . (2.6)
Then we multiply these components with c p = √ λ + 2µ and c s = √ µ, the pressure and the shear wave speeds respectively. Finally, we define I by

I = c s ∇ × ψ I + c p ∇φ I . (2.7)
We rigorously explain why should this new functional be better than the original one. We substantiate this argument with numerical illustrations.

In the sequel, we define respectively the Helmholtz decomposition operators H (2.8)

Time-reversal imaging analysis

In order to establish some results about time reversal using I and I, we use the following integral formulation based on the elastic Green's tensor.

Integral formulation

Let us introduce the outgoing Green's tensor G ω,0 associated to the elastic wave equation with a point source at 0:

(L λ,µ + ω 2 )G ω,0 (x) = -δ 0 I, x ∈ R d . (2.9)
It can be expressed in the form [1,[START_REF] Ammari | Polarization and Moment Tensors: with Applications to Inverse Problems and Effective Medium Theory[END_REF] G ω,0

(x) = 1 µκ 2 s κ 2 s G s ω,0 (x)I + D G s ω,0 -G p ω,0 (x) , x ∈ R d , (2.10)
where

I = (δ ij ) d i,j=1 , D = ( ∂ 2 ∂xi∂xj ) d i,j=1
, and κ 2 s = ω 2 µ and κ 2 p = ω 2 λ+2µ are the shear and the pressure wavenumbers respectively. Here, G α ω,0 (x) is the fundamental solution of the Helmholtz operator ∆ + κ 2 α in R d with α = p, s. For example when d = 3, we have

G α ω,0 (x) = exp(iκ α |x|) 4π|x| , α = p, s. (2.11)
The functional I(x) defined by (2.4) can be expressed in the form [START_REF] Ammari | Time reversal in attenuating acoustic media[END_REF]

I(x) = ℜe 1 2π ˆRd ˆR ω 2 ˆ∂Ω G ω (x, y)G ω (y, z)dσ(y) dωF(z) dz = 1 4π ˆRd ˆR ω 2 ˆ∂Ω G ω (x, y)G ω (y, z) + G ω (x, y)G ω (y, z) dσ(y) dωF(z) dz, (2.12)
where we have introduced the outgoing Green's tensor for a point source at y:

G ω (x, y) = G ω,0 (x -y). (2.13)
We also introduce the decomposition of G ω,0 into shear and compressional components as

G ω,0 (x) = G p ω,0 (x) + G s ω,0 (x), (2.14) with G p ω,0 = - 1 ω 2 DG p ω,0 and G s ω,0 = 1 ω 2 κ 2 s I + D G s ω,0 . (2.15)
We can extend Helmholtz operator H p and H s to tensors G as follows:

H p [G]p = H p [Gp] and H s [G]p = H s [Gp] for all vectors p.
Note that G p ω,0 and G s ω,0 satisfy, respectively

(L λ,µ + ω 2 )G p ω,0 = H p [-δ 0 I] and (L λ,µ + ω 2 )G s ω,0 = H s [-δ 0 I] . (2.16)
Consequently, the Helmholtz decomposition of I can be derived explicitly

I(x) = H p [I](x) + H s [I](x), (2.17)
with

H p [I](x) = 1 4π ˆRd ˆR ω 2 ˆ∂Ω G p ω (x, y)G ω (y, z) + G p ω (x, y)G ω (y, z) dσ(y) dωF(z) dz, and
H s [I](x) = 1 4π ˆRd ˆR ω 2 ˆ∂Ω G s ω (x, y)G ω (y, z) + G s ω (x, y)G ω (y, z) dσ(y) dωF(z) dz.
Finally, the integral formulation of the modified imaging functional I defined by (2.7) reads

I(x) = ℜe 1 2π ˆRd ˆR ω 2 ˆ∂Ω [c s G s ω (x, y) + c p G p ω (x, y)] G ω (y, z)dσ(y) dωF(z) dz .
(2.18)

Helmholtz-Kirchhoff identity

In order to approximate the integral formulation (2.18) we use a Helmholtz-Kirchhoff identity for elastic media. Some of the results presented in this subsection can be found in [START_REF] Wapenaar | Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlation[END_REF][START_REF] Wapenaar | Green's function representations for seismic interferometry[END_REF] in the context of elastodynamic seismic interferometry. Indeed, the elastodynamic reciprocity theorems (Propositions 2.1 and 2.5) are the key ingredient to understand the relation between the cross correlations of signals emitted by uncorrelated noise sources and the Green's function between the observation points.

Let us introduce the conormal derivative ∂u ∂ν (y), y ∈ ∂Ω, of the displacement field u at the surface ∂Ω in the outward unit normal direction n by

∂u ∂ν := λ(∇ • u)n + µ ∇u T + (∇u T ) T n, (2.19)
where T denotes the transpose. Note also that the conormal derivative tensor ∂G ∂ν means that for all constant vectors p,

∂G ω ∂ν p := ∂ [G ω p] ∂ν .
The following proposition is equivalent to [START_REF] Wapenaar | Green's function representations for seismic interferometry[END_REF]Eq. (73)]. Since our formulation is slightly different and this is the first building block of our theory we give its proof for consistency. Moreover, elements of the proof are used in Proposition 2.2.

Proposition 2.1. For all x, z ∈ Ω, we have

ˆ∂Ω ∂G ω (x, y) ∂ν G ω (y, z) -G ω (x, y) ∂G ω (y, z) ∂ν dσ(y) = 2iℑm{G ω (x, z)}. (2.20)
Proof. By reciprocity we have

G ω (y, x) = [G ω (x, y)] T . (2.21)
Additionally, in the homogeneous case we have G ω (y, x) = G ω (x, y), but we will not use this property here. Our goal is to show that for all constant vectors p and q, we have

ˆ∂Ω q • ∂G ω (x, y) ∂ν G ω (y, z)p -q • G ω (x, y) ∂G ω (y, z) ∂ν p dσ(y) = 2iq • ℑm{G ω (x, z)}p.
Taking scalar product of equations

(L λ,µ + ω 2 )G ω (y, x)q = -δ x q and (L λ,µ + ω 2 )G ω (y, z)p = -δ z p
with G ω (y, z)p and G ω (y, x)q respectively, subtracting the second result from the first, and integrating in y over Ω, we obtain

ˆΩ (G ω (y, z)p) • L λ,µ (G ω (y, x)q) -L λ,µ (G ω (y, z)p) • (G ω (y, x)q) dy = p • (G ω (z, x)q) -q • (G ω (x, z)p) = 2iq • ℑm{G ω (x, z)}p.
Using the form of the operator L λ,µ , this gives

2iq • ℑm{G ω (x, z)}p = λ ˆΩ (G ω (y, z)p) • {∇∇ • (G ω (y, x)q)} -(G ω (y, x)q) • ∇∇ • (G ω (y, z)p) dy +µ ˆΩ (G ω (y, z)p) • {(∆ + ∇∇•)(G ω (y, x)q)} -(G ω (y, x)q) • (∆ + ∇∇•)(G ω (y, z)p) dy.
We recall that, for two functions u, v : R d → R d , we have

(∆u + ∇(∇ • u)) • v = ∇ • (∇u T + (∇u T ) T )v - 1 2 (∇u T + (∇u T ) T ) • (∇v T + (∇v T ) T ), ∇(∇ • u) • v = ∇ • (∇ • u)v -(∇ • u)(∇ • v).
Therefore, we find

2iq • ℑm{G ω (x, z)}p = λ ˆΩ ∇ • [∇ • (G ω (y, x)q)](G ω (y, z)p) -∇ • [∇ • (G ω (y, z)p)](G ω (y, x)q) dy +µ ˆΩ ∇ • ∇(G ω (y, x)q) T + (∇(G ω (y, x)q) T ) T G ω (y, z)p -∇ • ∇(G ω (y, z)p) T + (∇(G ω (y, z)p) T ) T G ω (y, x)q dy.
Now, we use divergence theorem and the definition of the conormal derivative to get

2iq • ℑm{G ω (x, z)}p = λ ˆ∂Ω n • [∇ • (G ω (y, x)q)](G ω (y, z)p) -n • [∇ • (G ω (y, z)p)](G ω (y, x)q) dσ(y) +µ ˆΩ n • ∇(G ω (y, x)q) T + (∇(G ω (y, x)q) T ) T G ω (y, z)p -n • ∇(G ω (y, z)p) T + (∇(G ω (y, z)p) T ) T G ω (y, x)q dσ(y) = λ ˆ∂Ω (G ω (y, z)p) • {∇ • (G ω (y, x)q)n} -(G ω (y, x)q) • ∇ • (G ω (y, z)p)n dσ(y) +µ ˆ∂Ω (G ω (y, z)p) • ∇(G ω (y, x)q) T + (∇(G ω (y, x)q) T ) T n -(G ω (y, x)q) • ∇(G ω (y, z)p) T + (∇(G ω (y, z)p) T ) T n dσ(y) = ˆ∂Ω (G ω (y, z)p) • ∂G ω (y, x)q ∂ν -(G ω (y, x)q) • ∂G ω (y, z)p ∂ν dσ(y) = ˆ∂Ω q • ∂G ω (x, y) ∂ν G ω (y, z)p -q • G ω (x, y) ∂G ω (y, z) ∂ν p dσ(y),
which is the desired result. Note that for establishing the last equality we have used the reciprocity relation (2.21).

The proof of Proposition 2.1 uses only the reciprocity relation and the divergence theorem. Consequently, Proposition 2.1 also holds in a heterogeneous medium, as shown in [START_REF] Wapenaar | Green's function representations for seismic interferometry[END_REF].

The following proposition cannot be found in the literature, probably because its application in the context of seismic interferometry has not been identified. It is an important ingredient in the analysis of our improved imaging functional. Note that the proofs of Propositions 2.2 and 2.3 require the medium to be homogeneous (so that H s and H p commute with L λ,µ ), and we cannot expect these propositions to be true in a heterogeneous medium because of mode conversion between pressure and shear waves. Proposition 2.2. For all x, z ∈ Ω, we have

ˆ∂Ω ∂G s ω (x, y) ∂ν G p ω (y, z) -G s ω (x, y) ∂G p ω (y, z) ∂ν dσ(y) = 0. (2.22)
Proof. First, we recall that G p ω (y, x) and G s ω (y, x) are solutions of (2.16). We proceed as in the proof of the previous proposition to find:

ˆ∂Ω ∂G s ω (x, y) ∂ν G p ω (y, z) -G s ω (x, y) ∂G p ω (y, z) ∂ν dσ(y) = ˆΩ H s [-δ x I](y)G p ω (y, z) -G s ω (x, y)H p [-δ z I](y) dy = [H s [-δ 0 I] * G p ω (•, z)](x) -[G s ω (x, •) * H p [-δ 0 I]](z).
Using the fact that

G p ω = H p [G ω ] and H s H p = H p H s = 0 we get H s [H s [-δ 0 I] * G p ω (•, z)] = 0 and H p [H s [-δ 0 I] * G p ω (•, z)] = 0. Therefore, we conclude [H s [-δ 0 I] * G p ω (•, z)](x) = 0. Similarly, we have [G s ω (x, •) * H p [-δ 0 I]](z) = 0, which gives the desired result.
Finally the following proposition shows that the elastodynamic reciprocity theorem (Proposition 2.1) holds for each wave component in a homogeneous medium. Proposition 2.3. For all x, z ∈ R d and α = p, s,

ˆ∂Ω ∂G α ω (x, y) ∂ν G α ω (y, z) -G α ω (x, y) ∂G α ω (y, z) ∂ν dσ(y) = 2iℑm{G α ω (x, z)}. (2.23)
Proof. As both cases, α = p and α = s, are similar, we only provide a proof for α = p. For α = p, indeed, we have

ˆ∂Ω ∂G p ω (x, y) ∂ν G p ω (y, z) -G p ω (x, y) ∂G p ω (y, z) ∂ν dσ(y) = [H p [-δ 0 I] * G p ω (•, z)](x) -[G p ω (x, •) * H p [-δ 0 I]](z).
Using the fact that

G p ω (y, z) = G p ω,0 (y -z) = G p ω,0 (z -y), we can write [H p [-δ 0 I] * G p ω (•, z)](x) = [H p [-δ 0 I] * G p ω,0 (•)](x -z) and [G p ω (x, •) * H p [-δ 0 I]](z) = [G p ω,0 (•) * H p [-δ 0 I]](z -x) = [H p [-δ 0 I] * G p ω,0 (•)]](x -z).
Therefore,

ˆ∂Ω ∂G p ω (x, y) ∂ν G p ω (y, z) -G p ω (x, y) ∂G p ω (y, z) ∂ν dσ(y) = H p [2iℑm{G p ω (x, z)}] = 2iℑm{G p ω (x, z)},
where the last equality results from the fact that H s H p = 0.

Approximation of the conormal derivative

In this subsection, we derive an approximation of the conormal derivative ∂Gω(x,y) ∂ν , y ∈ ∂Ω, x ∈ Ω. In general this approximation involves the angles between the pressure and shear rays and the normal on ∂Ω. This approximation becomes simple when Ω is a ball with very large radius, since in this case all rays are normal to ∂Ω (Proposition 2.4). It allows us to use a simplified version of the Helmholtz-Kirchhoff identity in order to analyze the imaging functional I when Ω is a ball with large radius (Proposition 2.5).

Proposition 2.4. If n = y -x and |x -y| ≫ 1, then, for α = p, s, ∂G α ω (x, y) ∂ν = iωc α G α ω (x, y) + o 1 |x -y| . (2.24)
Proof. It is enough to show that for all constant vectors q,

∂G p ω (x, y)q ∂ν = iωc p G p ω (x, y)q + o 1 |x -y| and ∂G s ω (x, y)q ∂ν = iωc s G s ω (x, y)q + o 1 |x -y| . Pressure component: Recall that G p ω (x, y) = - 1 ω 2 DG p ω (x, y) = 1 c 2 p G p ω (x, y) y -x ⊗ y -x + o 1 |x -y| , so we have G p ω (x, y)q = 1 c 2 p G p ω (x, y) ( y -x • q) y -x + o 1 |x -y| . Therefore, ∂G p ω (x, y)q ∂ν = λ∇ y • (G p ω (x, y)q) n(y) + µ ∇ y (G p ω (x, y)q) T + (∇ y (G p ω (x, y)q) T ) T n(y) = y -x • q c 3 p iωG p ω (x, y) λ y -x • y -xn + 2µ( y -x ⊗ y -x)n + o 1 |y -x| = y -x • q c 3 p iωG p ω (x, y) λn + 2µ( y -x • n) y -x + o 1 |y -x| = y -x • q c 3 p iωG p ω (x, y) λ n -y -x + 2µ y -x • n -1 y -x +iωc p G p ω (x, y)q + o 1 |y -x| .
In particular, when n = yx, we have

∂G p ω (x, y)q ∂ν = iωc p G p ω (x, y)q + o 1 |y -x| .

Shear components:

As

G s ω (x, y) = 1 ω 2 κ 2 s I + D G s ω (x, y) = 1 c 2 s G s ω (x, y) I -y -x ⊗ y -x + o 1 |x -y| , we have G s ω (x, y)q = 1 c 2 s G s ω (x, y) q -y -x • q y -x + o 1 |x -y| . Therefore, ∂G s ω (x, y)q ∂ν = λ∇ y • (G s ω (x, y)q) n(y) + µ ∇ y (G s ω (x, y)q) T + (∇ y (G s ω (x, y)q) T ) T n(y). Now, remark that λ∇. (G s ω (x, y)q) n = λ iω c 3 s G s ω (x, y) q -y -x • q y -x • y -x n + o 1 |x -y| = o 1 |x -y| , and
µ ∇G s ω (x, y)q + ∇G s ω (x, y)q T n = µ iω c 3 s G s ω (x, y) q ⊗ y -x + y -x ⊗ q -2 y -x • q y -x ⊗ y -x n + o 1 |x -y| = µ iω c 3 s G s ω (x, y) y -x • n q + q • n y -x -2 y -x • q y -x • n y -x + o 1 |x -y| = µ iω c 3 s G s ω (x, y) y -x • n -1 q -y -x • q y -x +µ iω c 3 s G s ω (x, y) q • n -y -x • q y -x • n y -x iωc s G s ω (x, y) + o 1 |x -y| .
In particular, when n = yx, we have

∂G s ω (x, y)q ∂ν = iωc s G s ω (x, y)q + o 1 |y -x| .
The following is a direct consequence of Propositions 2.2, 2.3, and 2.4.

Proposition 2.5. Let Ω ⊂ R d be a ball with radius R. Then, for all x, z ∈ Ω sufficiently far from the boundary ∂Ω, we have

ℜe ˆ∂Ω G α ω (x, y)G α ω (y, z)dσ(y) ≃ 1 ωc α ℑm{G α ω (x, z)}, α = p, s, (2.25) ℜe ˆ∂Ω G s ω (x, y)G p ω (y, z)dσ(y) ≃ 0.
(2.26)

Analysis of the imaging functional I

In this subsection, we assume that Ω is a ball of radius R in R d and that the support, supp {F}, of F is sufficiently localized at the center of Ω so that for all x ∈ supp {F} and for all y ∈ ∂Ω

y -x = n(y) + o 1 |y -x| .
Then, we have the following theorem.

Theorem 2.6. Let x ∈ Ω be sufficiently far from the boundary ∂Ω and I be defined by (2.7). Then,

I(x) ≃ F(x). (2.27)
Proof. From (2.18) we have y). Proposition 2.5 allows us to write

I(x) = 1 4π ˆRd ˆR ω 2 ˆ∂Ω G ω (x, y)G ω (y, z) + G ω (x, y)G ω (y, z) dσ(y) dωF(z) dz, where G ω (x, y) = c s G s ω (x, y) + c p G p ω (x,
I(x) ≃ 1 4π ˆRd ˆR ω 2 ˆ∂Ω G ω (x, y)G ω (y, z) + G ω (x, y) G ω (y, z) dσ(y) dωF(z) dz ≃ 1 4π ˆRd ˆR ω 2 ˆ∂Ω G ω (x, y)G ω (y, z) + G ω (x, y) G ω (y, z) dσ(y) dωF(z) dz.
Proposition 2.4 then gives

I(x) ≃ 1 4π ˆRd ˆR -iω ∂G ω (x, y) ∂ν G ω (y, z) -G ω (x, y) ∂G ω (y, z) ∂ν dωF(z) dz ≃ 1 2π ˆRd ˆR ωℑm [G ω (x, z)] dωF(z) dz ≃ F(x).
The last approximation results from the identity

1 2π ˆR -iωG ω (x, z)dω = δ x (z)I,
which comes from the integration of the time-dependent version of (2.9) between t = 0 - and t = 0 + .

If the unweighted time-reversal imaging function I is used instead of I, then crossed terms remain. Using the same arguments as above we find

I(x) ≃ c s + c p c s c p 1 4π ˆRd ˆR ωℑm [(G p ω + G s ω )(x, z)] dωF(z) dz + c s -c p c s c p 1 4π ˆRd ˆR ωℑm [(G p ω -G s ω )(x, z)] dωF(z) dz ≃ c s + c p 2c s c p F(x) + c s -c p 2c s c p ˆRd B(x, z)F(z)dz, (2.28)
where

B(x, z) = 1 2π ˆR ωℑm [(G p ω -G s ω )(x, z)] dω (2.29)
is the operator that describes the error in the reconstruction of the source F obtained with I when c s = c p . In particular the operator B is not diagonal, which means that the reconstruction mixes the components of F.

Numerical simulations

Here we present numerical illustrations and describe our algorithms for numerical resolution of the source problem to show that I provides a better reconstruction than I.

Description of the algorithm

We describe the algorithm we used for the numerical resolution of the elastic wave equation in 2D:

     ∂ 2 u ∂t 2 (x, t) = [µ∆u + (λ + µ)∇(∇ • u)] (x, t), (x, t) ∈ R 2 × R, u(x, 0) = F(x) and ∂u ∂t (x, 0) = 0. (2.30)
This equation is computed on the box Q = [-L/2, L/2] 2 such that Ω ⊂ Q with periodic boundary conditions. We use a splitting spectral Fourier approach [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF] coupled with a perfectly matched layer (PML) technique [START_REF] Hastings | Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation[END_REF] to simulate a free outgoing interface on ∂Q.

With the notation u = (u 1 , u 2 ) and x = (x 1 , x 2 ), the elastic wave equation can be rewritten as a first order partial differential equation:

∂ t P = AP + BP, where P =     u 1 ∂ t u 1 u 2 ∂ t u 2     , A =     0 1 0 0 (λ + 2µ)∂ 2 x1 + µ∂ 2 x2 0 0 0 0 0 0 1 0 0 (λ + 2µ)∂ 2 x2 + µ∂ 2 x1 0     , and
B =     0 0 0 0 0 0 (λ + µ)∂ x1 ∂ x2 0 0 0 0 0 (λ + µ)∂ x1 ∂ x2 0 0 0     .
This equation is integrated via Strang's splitting method [START_REF] Strang | On the construction and comparison of difference schemes[END_REF]. This splitting approach is known to be of order 2 and reads exp(-

t(A + B)) = exp(-t/2B) exp(-tA) exp(-t/2B) + o(t 2 ).
The first operator A is then computed exactly in the spatial Fourier space. Indeed, the Fourier transform of P A (x, t) exp(-At)P (x) satisfies

               u A,1 (ξ, t) = cos ξ 2 λ,µ,1 t u 1 (ξ) + t sinc ξ 2 λ,µ,1 t ∂ t u 1 (ξ), ∂ t u A,1 (ξ, t) = cos ξ 2 λ,µ,1 t ∂ t u 1 (ξ) -ξ 2 λ,µ,1 sin ξ 2 λ,µ,1 t u 1 (ξ), u A,2 (ξ, t) = cos ξ 2 λ,µ,2 t u 2 (ξ) + t sinc ξ 2 λ,µ,2 t ∂ t u 2 (ξ), ∂ t u A,2 (ξ, t) = cos ξ 2 λ,µ,2 t ∂ t u 2 (ξ) -ξ 2 λ,µ,2 sin ξ 2 λ,µ,2 t u 2 (ξ), with ξ 2 λ,µ,1 = 4π 2 (λ + 2µ)ξ 2 1 + µξ 2 2 , ξ 2 λ,µ,2 = 4π 2 (λ + 2µ)ξ 2 2 + µξ 2 1
, and sinc(t) = sin(t)/t.

The second operator B is also integrated exactly. We have

P B (x, t) = exp(-Bt)P (x) =     u 1 (x) ∂ t u 1 (x) -t (λ + µ) ∂ x1 ∂ x2 u 2 (x) u 2 (x) ∂ t u 2 (x) -t (λ + µ) ∂ x1 ∂ x2 u 1 (x)     .
This global algorithm appears to be stable under a classical condition of the form

δ t ≤ c(λ, µ)δ 2 x ,
where δ t and δ x denote respectively the time and the spatial step of discretization. Here c(λ, µ) is a constant which depends only on Lamé coefficients λ and µ.

The functional I(x) requires also a Helmholtz decomposition algorithm. As the support of the function I(x) is included in Ω ⊂ Q, we apply a Dirichlet boundary condition on ∂Q. This decomposition is numerically obtained with a fast algorithm [START_REF] Wiegmann | Fast Poisson, fast Helmholtz and fast linear elastostatic solvers on rectangular parallelepipeds[END_REF] based on a symmetry principle and a Fourier-Helmholtz decomposition algorithm.

Experiments

In the sequel, for numerical illustrations, Ω is taken as a unit disc. Its boundary is discretized by 1024 sensors. Each solution of elastic wave equation is computed over (x, t) ∈ [-L/2, L/2] 2 × [0, T ] with L = 4 and T = 2. We use a step of discretization given by dt = T /2 13 and dx = L/2 9 .

Figure 1 presents a first experiment with Lamé parameters (λ, µ) = (1, 1). The first (top) line corresponds to the 2 components of the initial source F. The second line corresponds to the data g(y, t) = u(y, t) recorded over (y, t) ∈ ∂Ω × [0, T ]. Note that the shear and pressure waves are mixed in the recorded signal and it seems difficult to separate them. The third line corresponds to the imaging functional I(x). This example clearly shows that the reconstruction of the source F is not so accurate with classical time-reversal imaging. Moreover, in the last row, the plots represent the modified imaging functional I(x) where the reconstruction is quite better and nearly optimal.

Figure 2 shows another example with different Lamé parameters (λ, µ) = (10, 1). The same conclusion holds.

In Figure 3, we use a "less localized" (large) source distribution. We observe some artifacts in the reconstruction of the imaging functional I. We can also observe, from the recorded data, that the pressure and shear waves are very much "mixed" with each other. We expect that the artifacts in the reconstruction are the consequence of such mixing. In this situation, we do not have a real orthogonality between the two waves on ∂Ω.

To conclude, I provides a better reconstruction of the sources than I. However, in certain cases, the reconstructions by I are not optimal, and need further improvements.

Time-reversal algorithm for viscoelastic media

In this section, we investigate the source inverse problem in an elastic attenuating medium. We provide an efficient regularized time-reversal imaging algorithm which corrects for the leading-order terms of the attenuation effect. Consider the viscoelastic wave equation in a open medium Ω ⊂ R d with d = 2, 3, i.e.,

     ∂ 2 u a ∂t 2 (x, t) -L λ,µ u a (x, t) - ∂ ∂t L η λ ,ηµ u a (x, t) = dδ 0 (t) dt F(x), (x, t) ∈ R d × R, u a (x, t) = ∂u a ∂t (x, t) = 0, t < 0, (3.1)
where the viscosity parameters η µ and η λ are positive constants and account for losses in the medium. As in the acoustic case [START_REF] Ammari | Time reversal in attenuating acoustic media[END_REF], the strategy of time reversal is to consider the functional

I a (x) = ˆT 0 v s,a (x, T )ds, (3.2)
where v s,a should be the solution of the adjoint (time-reversed) viscoelastic wave equation, i.e.,

           ∂ 2 v s,a ∂t 2 (x, t) -L λ,µ v s,a (x, t) + ∂ ∂t L η λ ,ηµ v s,a (x, t) = dδ s (t) dt g a (x, T -s)δ ∂Ω (x), (x, t) ∈ R d × R, v s (x, t) = ∂v s ∂t (x, t) = 0 x ∈ R d , t < s. (3.3)
Further, the idea is to enhance image resolution using I a , as in purely elastic media, where

I a = c s H s [I a ] + c p H p [I a ] .
Unfortunately, the adjoint viscoelastic problem is severely ill-posed. Indeed, the frequency component is exponentially increasing due to the presence of the anti-damping term +∂ t L η λ ,ηµ v s,a , which induces instability. Therefore, we need to regularize the adjoint problem by truncating high-frequency components either in time or in space.

Let us introduce the outgoing Green's tensor G a,ω associated to the viscoelastic wave equation

(L λ,µ -iωL η λ ,ηµ + ω 2 )G a,ω (x, y) = -δ y (x)I, x, y ∈ R d , (3.4)
and let G -a,ω be the adjoint viscoelastic Green's tensor, that is, the solution to

(L λ,µ + iωL η λ ,ηµ + ω 2 )G -a,ω (x, y) = -δ y (x)I, x, y ∈ R d . (3.5)
We introduce an approximation v s,a,ρ of the adjoint wave v s,a by

v s,a,ρ (x, t) = - 1 2π ˆ|ω|≤ρ ˆ∂Ω iωG -a,ω (x, y)g a (y, T -s)dσ(y) exp (-iω(t -s)) dω, (3.6
) where ρ ∈ R + is the cut-off parameter. The regularized time-reversal imaging functional defined by

I a,ρ (x) = ˆT 0 v s,a,ρ (x, T )ds, (3.7)
can be written as

I a,ρ (x) = ˆ∂Ω ˆT 0 ∂ ∂t G -a,ρ (x, y, T -s)g a (y, T -s)ds dσ(y), (3.8)
where

G -a,ρ (x, y, t) = 1 2π ˆ|ω|≤ρ G -a,ω (x, y) exp (-iωt) dω. (3.9)
Remark 3.1. Let S ′ be the space of tempered distributions, i.e., the dual of the Schwartz space S of rapidly decreasing functions [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators. I. Distribution Theory and Fourier Analysis[END_REF]. The function v s,a,ρ (x, t) can be identified as the solution of the following wave equation:

∂ 2 v s,a,ρ ∂t 2 (x, t) -L λ,µ v s,a,ρ (x, t) + ∂ ∂t L η λ ,ηµ v s,a,ρ (x, t) = S ρ dδ s (t) dt g a (x, T -s)δ ∂Ω (x), (3.10)
where the operator S ρ is defined on the space S ′ by

S ρ [ψ] (t) = 1 2π ˆ|ω|≤ρ exp (-iωt) ψ(ω)dω, (3.11)
with ψ(ω) = ˆR ψ(t) exp(iωt) dt.

(3.12)

Green's tensor in viscoelastic media

As in Section 2, we decompose G ±a,ω in the form

G ±a,ω = G s ±a,ω + G p ±a,ω , (3.13)
where G s ±a,ω and G p ±a,ω are respectively the fundamental solutions of

(L λ,µ ∓ iωL η λ ,ηµ + ω 2 )G α ±a,ω (x, y) = H α [-δ y I] , α = p, s. (3.14)
Let us also introduce the decomposition of the operator L λ,µ into shear and pressure components as

L λ,µ = L s λ,µ + L p λ,µ , and L η λ ,ηµ = L s η λ ,ηµ + L p η λ ,ηµ , (3.15)
where

L s λ,µ u = c 2 s [∆u -∇(∇ • u)] and L p λ,µ u = c 2 p ∇(∇ • u), (3.16)
and

L s η λ ,ηµ u = ν 2 s [∆u -∇(∇ • u)] and L p η λ ,ηµ u = ν 2 p ∇(∇ • u). (3.17)
Here, ν 2 s = η µ and ν 2 p = η λ + 2η µ . Therefore, the tensors G s ±a,ω and G p ±a,ω can also be seen as the solutions of

(L α λ,µ ∓ iωL α η λ ,ηµ + ω 2 )G α ±a,ω (x, y) = H α -δ y I (x), α = p, s. (3.18)
The corrected regularized time-reversal imaging functional defined by

I a,ρ = c s H s [I a,ρ ] + c p H p [I a,ρ ] , (3.19)
is then given by

I a,ρ (x) = ˆ∂Ω ˆT 0 ∂ ∂t c p G p -a,ρ (x, y, T -s) + c s G s -a,ρ (x, y, T -s) g a (y, T -s)ds dσ(y), (3.20) where G α -a,ρ (x, y, t) = 1 2π ˆ|ω|≤ρ G α -a,ω (x, y) exp (-iωt) dω, α = p, s. (3.21)
In the next subsection we express the relationship between the data g a and the ideal measurements g obtained in the non-attenuated case. Doing so, we prove with the help of a new version of the Helmholtz-Kirchhoff identity that a regularized image of the source F can be obtained.

Attenuation operator and its asymptotics

Recall that u and u a are respectively the solutions of the wave equations

∂ 2 u ∂t 2 (x, t) -L λ,µ u(x, t) = dδ 0 (t) dt F(x), (3.22)
and

∂ 2 u a ∂t 2 (x, t) -L λ,µ u a (x, t) - ∂ ∂t L η λ ,ηµ u a (x, t) = dδ 0 (t) dt F(x), (3.23)
with the initial conditions

u(x, t) = u a (x, t) = ∂u ∂t (x, t) = ∂u a ∂t (x, t) = 0, t < 0. (3.24)
We decompose u and u a as

u = u s + u p = H s [u] + H p [u] and u a = u s a + u p a = H s [u a ] + H p [u a ]. (3.25)
The Fourier transforms u α ω and u α a,ω of the vector functions u α and u α a are respectively solutions of

ω 2 + L α λ,µ u α ω = iωH α [F] and κ α (ω) 2 + L α λ,µ u α ω = i κ α (ω) 2 ω H α [F] , α = p, s, (3.26) where κ
α (ω) = ω 1 -iων 2 α /c 2 α . (3.27)
In particular, it implies that

u s a = A ν 2 s /c 2 s [u s ] and u p a = A ν 2 p /c 2 p [u p ] , (3.28)
where A a , for a > 0, is the attenuation operator

A a [φ] (t) = 1 2π ˆR κ a (ω) ω ˆR φ(s) exp{iκ a (ω)s}ds exp{-iωt}dω, (3.29)
with κ a (ω) = ω √ 1 -iωa .

We also define the operator A -a,ρ by

A -a,ρ [φ] (t) = 1 2π ˆR+ φ(s) ˆ|ω|≤ρ κ -a (ω) ω exp{iκ -a (ω)s} exp{-iωt}dω ds, (3.30)
which is associated with

κ -a (ω) = ω √ 1 + iωa . Moreover, its adjoint operator A * -a,ρ reads A * -a,ρ [φ] (t) = 1 2π ˆ|ω|≤ρ κ -a (ω) ω exp{iκ -a (ω)t} ˆR+ φ(s) exp{-iωs}ds dω. (3.31)
Then, we have following results from [START_REF] Ammari | Time reversal in attenuating acoustic media[END_REF][START_REF] Ammari | Photoacoustic imaging for attenuating acoustic media[END_REF].

Proposition 3.2.

• Let φ(t) ∈ S ([0, ∞)), S being the Schwartz space. Then

A a [φ] (t) = φ(t) + a 2 (tφ ′ ) ′ (t) + o(a) as a → 0. (3.32) • Let φ(t) ∈ D ([0, ∞)), where D ([0, ∞)) is the space of C ∞ -functions of compact support on [0, ∞).
Then, for all ρ,

A * -a,ρ [φ] (t) = S ρ [φ] (t) - a 2 S ρ (tφ ′ ) ′ (t) + o(a) as a → 0. (3.33) • Let φ(t) ∈ D ([0, ∞)).
Then, for all ρ,

A * -a,ρ A a [φ] (t) = S ρ [φ] (t) + o(a) as a → 0. (3.34)
We extend the operators A a , A -a,ρ and A * -a,ρ to tensors G, i.e., for all vectors

p ∈ R d , A a [Gp] = A a [G]p, A -a,ρ [Gp] = A -a,ρ [G]p, and A * -a,ρ [Gp] = A * -a,ρ [G]p.
By the definition of the operators A a and A -a,ρ , we have for α = p, s:

∂G α a ∂t (x, y, t) = A ν 2 α /c 2 α ∂G α ∂t (x, y, •) (t), (3.35) ∂G α -a,ρ ∂t (x, y, t)
= A -ν 2 α /c 2 α ,ρ ∂G α ∂t (x, y, •) (t). (3.36)
Recall that g = u and g a = u a on ∂Ω × R. It then follows from (3.28) that, for α = p, s,

A * -ν 2 α /c 2 α ,ρ A ν 2 α /c 2 α [g α a ] = S ρ [g α ] + o(a), (3.37)
where

g α a = H α [g a ], g α = H α [g]. (3.38)
Identity (3.37) proves that A * -ν 2 α /c 2 α ,ρ is an approximate inverse of A ν 2 α /c 2
α . Moreover, it plays a key role in showing that the regularized time reversal algorithm provides a firstorder correction of the attenuation effect.

Helmholtz-Kirchhoff identity in an attenuating media

In this subsection we derive a new Helmholtz-Kirchhoff identity in elastic attenuating media. For doing so, let us introduce the conormal derivatives ∂u ∂νa and ∂u ∂ν-a as follows:

∂u ∂ν ±a := λ(∇ • u)n + µ(∇u T + (∇u T ) T )n ∓ iω η λ (∇ • u)n + η µ (∇u T + (∇u T ) T )n .
(3.39) Note also that for a tensor G the conormal derivative ∂G ∂ν±a means that for all constant vectors p, ∂G ∂ν ±a

p := ∂ [Gp] ∂ν ±a . (3.40)
The following properties hold.

Proposition 3.3. For all x, z ∈ Ω, we have

ˆ∂Ω ∂G s -a,ω (x, y) ∂ν -a G p a,ω (y, z) -G s -a,ω (x, y) ∂G p a,ω (y, z) ∂ν a dσ(y) = 0. (3.41)
Proof. Note that

J := ˆ∂Ω ∂G s -a,ω (x, y) ∂ν -a G p a,ω (y, z) -G s -a,ω (x, y) ∂G p a,ω (y, z) ∂ν a dσ(y) = ˆ∂Ω ∂G s -a,ω (x, y) ∂ν -a G p a,ω (y, z) -G s -a,ω (x, y) ∂G p a,ω (y, z) ∂ν -a dσ(y) = ˆΩ L λ,µ G s -a,ω (x, y) + iωL η λ ,ηµ G s -a,ω (x, y) G p a,ω (y, z) dy -ˆΩ G s -a,ω (x, y) L λ,µ G p a,ω (y, z) + iωL η λ ,ηµ G p a,ω (y, z) dy, = ˆΩ L λ,µ G s -a,ω (x, y) + iωL η λ ,ηµ G s -a,ω (x, y) G p a,ω (y, z) dy -ˆΩ G s -a,ω (x, y) L λ,µ G p a,ω (y, z) -iωL η λ ,ηµ G p a,ω (y, z) dy.
Since G p a,ω (x, y) and G s -a,ω (x, y) are solutions of equations (3.14) with α = p, s, respectively, it follows that

J = [H s [-δ 0 I] * G p a,ω (•, z)](x) -[G s -a,ω (x, •) * H p [-δ 0 I]](z).
As in the proof of Proposition 2. We now give an approximation of the attenuating co-normal derivative.

Proposition 3.4. If n = yx, then for α = p, s, we have

∂G α ±a,ω (x, y) ∂ν ±a ≃ ic α ω 2 κ α ∓ (ω) G α ±a,ω (x, y), (3.42)
where

κ α ∓ (ω) = ω 1 ∓ iων 2 α /c 2 α . (3.43)
Proof. Indeed, notice that

G α ±a,ω (x, y) = κ α ∓ (ω) ω 2 G α κ α ∓ (ω) (x, y), α = p, s. (3.44)
Then, from Proposition 2.4, we obtain

∂G α ±a,ω (x, y) ∂ν ±a ≃ κ α ∓ (ω) ω 2 ic 2 α κ α ∓ (ω) c α G α κ α ∓ (ω) (x, y) ∓ iω iν 2 α κ α ∓ (ω) c α G α κ α ∓ (ω) (x, y) ≃ ic α κ α ∓ (ω) 1 ∓ iων 2 α /c 2 α G α ±a,ω (x, y) ≃ ic α ω 2 κ α ∓ (ω)
G α ±a,ω (x, y).

In particular, the following estimate holds as a direct consequence of Propositions 3.3 and 3.4.

Analysis of the regularized time-reversal imaging algorithm

The aim of this subsection is to justify that the regularized time-reversal imaging functional I a,ρ provides a correction of the attenuation effect.

Theorem 3.6. The regularized time-reversal imaging functional I a,ρ satisfies

I a,ρ (x) = I ρ (x) + o ν 2 s /c 2 s + ν 2 p /c 2 p , (3.46)
I ρ (x) := ˆ∂Ω ˆT 0 c s ∂ ∂t G s (x, y, s) + c p ∂ ∂t G p (x, y, s) S ρ [g(y,
•)] (s)ds dσ(y), (3.47) where S ρ is defined by (3.11).

Proof. We can decompose the functional I a,ρ as follows: where

I a,ρ (x) = ˆ∂Ω ˆT 0 ∂ ∂t c p G p -a,ρ (x, y, T -s) + c s G s -a,ρ (x, y, T -s) g a (y, T -s)ds dσ(y) = ˆ∂Ω ˆT 0 ∂ ∂t c p G p -a,ρ (x, y, s) + c s G s -a,ρ (x,
I αβ a,ρ (x) = ˆ∂Ω ˆT 0 ∂ ∂t c α G α -a,ρ (x, y, s) g β a (y, s)ds dσ(y), α, β ∈ {p, s}.
Similarly we can decompose the functional I ρ as

I ρ (x) = I ss ρ (x) + I sp ρ (x) + I ps ρ (x) + I pp ρ (x),
with

I αβ ρ (x) = ˆ∂Ω ˆT 0 ∂ ∂t [c α G α (x, y, s)] S ρ [g β (y,
•)](s)ds dσ(y), α, β ∈ {p, s}.

The first term I ss a,ρ (x) satisfies

I ss a,ρ (x) = ˆ∂Ω ˆT 0 A -ν 2 s /c 2 s ,ρ ∂ ∂t [c s G s (x, y, s)] A ν 2 s /c 2 s [g s (y, •)] (s)ds dσ(y) = ˆ∂Ω ˆT 0 ∂ ∂t [c s G s (x, y, s)] A * -ν 2 s /c 2 s ,ρ A ν 2 s /c 2 s [g s (y, •)] (s)ds dσ(y) = ˆ∂Ω ˆT 0 ∂ ∂t [c s G s (x, y, s)] S ρ [g s (y, •)] (s)ds dσ(y) + o ν 2 s /c 2 s = I ss ρ (x) + o ν 2 s /c 2 s ,
by using Proposition 3.2. Similarly, we get

I pp a,ρ (x) = I pp ρ (x) + o ν 2 p /c 2 p .
Moreover, the coupling terms I sp a,ρ and I ps a,ρ vanish. Indeed, thanks to Proposition 3.5, we have

I sp a,ρ (x) = 1 2π ˆRd ˆ|ω|<ρ ω 2 ˆ∂Ω c s G s -a,ω (x, y) G p a,ω (y, z) dσ(y) dωF(z) dz ≃ 0,
and

I ps a,ρ (x) = 1 2π ˆRd ˆ|ω|<ρ ω 2 ˆ∂Ω c p G p -a,ω (x, y) G s a,ω (y, z) dσ(y) dωF(z) dz ≃ 0.
Proposition 2.5 shows that we also have

I sp ρ (x) = 1 2π ˆRd ˆ|ω|<ρ ω 2

ˆ∂Ω

[c s G s ω (x, y)] G p ω (y, z) dσ(y) dωF(z) dz ≃ 0,

I ps ρ (x) = 1 2π ˆRd ˆ|ω|<ρ ω 2 ˆ∂Ω [c p G p ω (x, y)] G s ω (y, z) dσ(y) dωF(z) dz ≃ 0,
which concludes the proof.

It is straightforward to check that

I ρ (x) ρ→∞ -→ I(x) ≃ F(x), (3.48)
by Theorem 2.6. Therefore, I a,ρ provides a first-order correction in terms of ν 2 s /c 2 s + ν 2 p /c 2 p of the attenuation effect. Moreover, the imaging functional I a,ρ can be seen as the timereversal functional I defined by (2.7) applied to A * -ν 2 α /c 2 α ,ρ g α a , α = p, s. As shown in (3.34), the regularized operators A * -a,ρ give a first-order approximation of the inverse of A a for a = ν 2 α /c 2 α . It would be very interesting to construct higher-order reconstructions in terms of the attenuation effect using higher-order approximations of the inverse of the operator A a . The problem is more challenging than the one in the scalar case [START_REF] Ammari | Time reversal in attenuating acoustic media[END_REF] because of the coupling between the shear and pressure components. Note finally that, if one applies the time-reversal functional I to the data g a directly, then one finds

I(x) = ˆ∂Ω ˆT 0 ∂ ∂t [c p G p (x, y, s) + c s G s (x, y, s)] [g a (y, •)] (s) ds dσ(y) = 1 2π ˆRd ˆR ω 2 ˆ∂Ω [c p G p ω (x, y) + c s G s ω (x, y)] G a,ω (y, z)dσ(y) dωF(z)dz, (3.49)
which gives an error of the order of ν 2 s /c 2 s + ν 2 p /c 2 p as can be seen from the expansion (3.32).

Numerical simulations

In this section we present numerical illustrations and describe our algorithms for numerical resolution of the source problem to show that I a,ρ provides a better reconstruction than Ĩ, where the attenuation effect is not taken into account.

Description of the algorithm

In the expression of I a,ρ , the solution v s,a,ρ (x, t) is very difficult to obtain numerically. Therefore, we prefer to regularize the problem by truncating high-frequency components in space instead of time, in contrast with our theoretical analysis. This can be seen as an approximation ṽs,a,ρ (x, t) of v s,a,ρ (x, t) defined as the solution of

∂ 2 ṽs,a,ρ ∂t 2 (x, t) -L λ,µ ṽs,a,ρ (x, t) + ∂ ∂t L η λ ,ηµ ṽs,a,ρ (x, t) = ∂δ s (t) ∂t X ρ [g a (•, T -s)δ ∂Ω ] (x),
where the operator X ρ is defined by

X ρ [g a (•, T -s)δ ∂Ω ] (x) =
ˆ|k|≤ρ ˆ∂Ω g a (y, T -s)e -2iπk•y dσ(y) e 2iπk•x dk δ ∂Ω (x).

The operator X ρ , as the operator S ρ , truncates high frequencies but in the space variable.

To compute the solution of the viscoelastic wave equation in two dimensions

∂ 2 u a ∂t 2 (x, t) -L λ,µ u a (x, t) ± L η λ ,ηµ u a (x, t) = 0,
we use the same algorithm presented in the unattenuated case, i.e., we use a larger box Ω ⊂ Q = [-L/2, L/2] 2 with periodic boundary condition and again a splitting spectral Fourier approach coupled with a PML technique to simulate a free outgoing interface on ∂Q.

Experiments

In the sequel, for numerical illustrations, Ω is taken to be a unit disc. Its boundary is discretized by 2 11 sensors. Each solution of elastic wave equation is computed over (x, t) ∈ [-L/2, L/2] 2 × [0, T ] with L = 4 and T = 2. We use a step of discretization given by dt = T /2 13 and dx = L/2 9 .

Figure 4 presents a first experiment with Lamé parameters (λ, µ) = (1, 1) and attenuation coefficients (η λ , η µ ) = (0.0002, 0.0002). The first line corresponds to the two components of the initial source F. The second line corresponds to the reconstruction of F without taking into account the attenuation effect. The imaging functional I(x) appears to be blurred due to coupling effects. The three last lines correspond to a reconstruction of F using the imaging functional I a,ρ for different values of ρ. We clearly observe a better reconstruction of the source F than using the functional I provided that the regularization parameter ρ is chosen appropriately large in order to insure the good resolution of the reconstruction.

  p and H s by H p [I] := ∇φ I and H s [I] := ∇ × ψ I .

Figure 1 : 14 - 1 -

 1141 Figure 1: Test 1: Source reconstruction; Comparison between the imaging functionals I and I with Lamé constants (λ, µ) = (1, 1). From top to bottom: First line: the source F; Second line: recorded data g(y, t); Third line: imaging functional I; Last line: imaging functional I. 14

Figure 2 :

 2 Figure 2: Test 2: Source reconstruction; Comparison between the imaging functionals I and I with Lamé constants (λ, µ) = (10, 1); First line: the source F; Second line: recorded data g(y, t); Third line: imaging functional I; Last line: imaging functional I. 15

Figure 3 :

 3 Figure 3: Test 3: Source reconstruction; Comparison between the imaging functionals I and I with Lamé constants (λ, µ) = (1, 1) and less localized source than in Figure 1; First line: the source F; Second line: recorded data g(y, t); Third line: imaging functional I; Last line: imaging functional I. 16

  2 one can show that [H s [-δ 0 I] * G p a,ω (•, z)](x) = 0 and [G s -a,ω (x, •) * H p [-δ 0 I]](z) = 0, which completes the proof of the proposition.

Proposition 3 . 5 .

 35 Let Ω ⊂ R d be a ball with radius R. Then for all x, z ∈ Ω sufficiently far from boundary ∂Ω, we haveℜe ˆ∂Ω G s -a,ω (x, y)G p a,ω (y, z)dσ(y) ≃ ℜe ˆ∂Ω G p -a,ω (x, y)G s a,ω(y, z)dσ(y) ≃ 0. (3.45)

Figures 5 and 6

 6 Figures 5 and 6 present two other examples of reconstruction using I a,ρ . The same observation holds.

Figure 4 :

 4 Figure 4: Test 4: Comparison between I and I a,ρ in a viscoelastic medium; The parameters are (λ, µ) = (1, 1) and (ν 2 s /c 2 s , ν 2 p /c 2 p ) = (0.0002, 0.0002); First (top) line: initial condition; Second line: without correction of attenuation; Last lines: with I a,ρ and ρ = 15, 20, 25.

Figure 5 :

 5 Figure 5: Test 5: Comparison between I and I a,ρ in a viscoelastic medium; The parameters are (λ, µ) = (1, 1) and (ν 2 s /c 2 s , ν 2 p /c 2 p ) = (0.00005, 0.00005); First line: initial condition; Second line: without correction of attenuation; Last lines: with I a,ρ and ρ = 15, 20, 25.

Figure 6 :

 6 Figure 6: Test 6: Comparison between I and I a,ρ in a viscoelastic medium; The parameters are (λ, µ) = (1, 1) and (η λ , η µ ) = (0.00005, 0.00005); First line: initial source F; Second line: reconstruction of F using I; Three last lines: reconstruction of F using I a,ρ with respectively ρ = 25, 30, 35.

* This work was supported by the ERC Advanced Grant Project MULTIMOD-267184.