Monsieur Stefan Catheline

Directeur De Recherche

Monsieur Nori Nakata

Monsieur Kees Wapenaar

Julia De

Chantal Van Dinther

Leonard Seydoux

William Lam

Revealing the structure of fault zones provides insights required to assess seismic hazards. Reflection imaging methods, in particular migration methods, are commonly used to image the subsurface. An accurate velocity model is often needed to properly locate the reflectors in depth. Imaging fault zones is challenging due to the complex distribution of velocity. Also, most reflection techniques image fault zones indirectly from the discontinuity of geological layers. Wave diffraction by small-scale heterogeneities, such as cracks and fractures is generally neglected. In this thesis, we propose a matrix imaging approach for geophysical exploration that handles the velocity variation in the crust and allows to resolve heterogeneities of the order of the wavelength throughout the fault zone.

First, I would like to thank my family. Without their love and support over the years none of this would have been possible. I am thankful for everything they have helped me achieve.

Next, I would like to thank my supervisors Michel Campillo and Alexandre Aubry

The reflection matrix approach is inspired from previous studies in ultrasound and optical imaging of complex media. First, ambient noise cross-correlations are used to retrieve the reflection matrix associated with a dense array of geophones. The body wave components of this matrix contains all the information available on the medium at depth. Second, the presented approach does not require a detailed velocity model of the subsurface. A set of matrix operation is actually applied to compensate for the mismatch between the actual wave velocity distribution and its approximate model.

In a first application, we use ambient noise records in the frequency band [10 20] Hz from a dense array in order image the San Jacinto Fault, California, at small scale with an horizontal resolution of 80 m. ZZ cross-correlations are computed and arranged as 2D matrix. By applying time delays, the response matrix is projected to depth, using a homogeneous velocity model. A focused reflection matrix is obtained, that contains the impulse response between virtual sources and receivers at depth. From this matrix, the image of i the medium is built and the resolution of the image can be quantified. Variations between the considered model and the reality result in phase distortions, i.e aberrations, that have detrimental effects on the image of the medium. We develop an aberration correction process that allows to compensate for these distortions by introducing a novel operator, the so-called distortion matrix. 3D images of the first 4 kilometers of the crust are retrieved.

These images reveal the backscattered intensity generated by the heterogeneities in the medium. The location and reflectivity of scatterers are retrieved with a resolution eight times better than the one in free space. Differences in the scattering between the Northwest and the Southeast of the fault were reported with an intense localized damage zone in the Southeast.

In a second application, we image the large scale structure of North Anatolian Fault using [0.1 0.5] Hz horizontal cross-correlations computed between 73 pairs of stations.

A multi-layered velocity model is considered. A local correction of the phase distortions is performed. The scattering structure of the crust and the upper mantle is revealed. Differences in the Moho depth are reported, with a step below the northern branch. Strong scattering is observed in the region lying along the northern strand of the fault, coinciding with the limit of the lithospheric blocks. The scattering in the North extends to 60 km depth, suggesting a shear zone that penetrates in the upper mantle beneath the northern strand. The scattering also reveals the deep structure of the main geological blocks.

In the last part of this thesis, we propose an approach that paves the route towards a 3D passive tomography of the body wave velocity. Finally, all the presented applications confirm the efficiency of the reflection matrix approach in revealing the structure of the subsurface. It provides new insights into the scattering distribution in the Earth. It can be applied to any scale, scattering regime, and frequency bandwith, provided that the spatial sampling of the geophones' array satisfies the Nyquist criterion.

Keywords: wave scattering and diffraction, fault zone structure, passive imaging, reflection matrix, body waves.

Resumé

La structure des zones de faille a une directe implication sur l'évaluation des risques sismiques. Les techniques d'imagerie par réflexion, en particulier les méthodes de migration, sont couramment utilisées pour imager le sous-sol. Un modèle de vitesse précis est souvent nécessaire pour localiser correctement les réflecteurs en profondeur. L'imagerie des zones de faille est un défi en raison de la distribution complexe de la vitesse des ondes de volume en leur sein. De plus, la plupart des techniques de réflexion imagent les zones de faille indirectement à partir de la discontinuité des couches géologiques. Les hétérogénéités à petite échelle, telles que les fissures et les fractures, sont généralement négligées. Dans cette thèse, nous proposons une approche d'imagerie matricielle pour l'exploration géophysique qui compense les variations de vitesse dans la croûte pour imager les hétérogénéités de cette dernière avec une résolution de l'ordre de la longueur d'onde autour des zones de faille. L'imagerie matricielle est inspirée de précédents travaux en imagerie ultrasonore et optique. D'une part, les corrélations du bruit ambiant sont utilisées pour construire la matrice de réflexion associée à un réseau dense de géophones. La composante des ondes de volume dans cette matrice contient toute l'information disponible sur le milieu en profondeur. D'autre part, l'approche présentée ne nécessite pas un modèle de vitesse détaillé du sous-sol. Un ensemble d'opérations matricielles est appliqué pour compenser les différences entre la distribution réelle de la vitesse des ondes et son modèle approximative.

Dans une première application, nous utilisons des enregistrements de bruit ambiant dans la bande de fréquence [10 20] Hz provenant d'un réseau dense de géophones afin iii d'imager la faille de San Jacinto, Californie, à petite échelle avec une résolution transversal de 80 m. Les corrélations ZZ sont calculées et organisées en matrice 2D. En appliquant des lois de retard temporels sur les réponses impulsionnelles entre géophones, la matrice de réponse est projetée en profondeur, en se basant sur un modèle de vitesse homogène.

Une matrice de réflexion focalisée est obtenue et contient les réponses impulsionnelles entre un ensemble de sources et récepteurs virtuels à chaque profondeur. A partir de cette matrice, l'image du milieu est construite et ses aberrations peuvent être quantifiées. La différence entre le modèle de vitesse considéré et la réalité entraînent en effet des distorsions de phase, qui dégradent considérablement la résolution et le contraste de l'image.

Un processus de correction matricielle des aberrations permet ensuite de compenser ces distorsions en introduisant un nouvel opérateur: la matrice distorsion. Une image 3D des quatres premiers kilomètres de la croûte terrestre est obtenue avec une résolution huit fois plus fine que celle attendue en milieu homogène. L'imagerie matricielle se sert du milieu comme d'une lentille diffusante, permettant ainsi d'élargir la taille effective du réseau de géophone. L'image obtenue révèle l'intensité rétrodiffusée par les hétérogénéités du milieu. Des différences significatives de réflectivité sont observées entre le Nord-Ouest et le Sud-Est de la faille avec une zone de dommages intense et localisée dans le Sud-Est qui suggère une génération continue de dommages aux roches au Sud-Est.

Dans un deuxième cas d'étude, nous imageons la structure à grande échelle de la faille Nord Anatolienne avec une résolution transverse de 2 km, en utilisant les corrélations horizontales [0.1 0.5] Hz calculées entre 73 paires de stations. Un modèle de vitesse multicouches est ici considéré. Une correction locale des distorsions de phase est effectuée.

La structure de la croûte et du manteau supérieur est révélée. Des différences dans la profondeur du Moho sont mises en lumière, avec un saut sous la branche nord. Une forte réflectivité du sous-sol est observée dans la région située le long de la branche Nord de la faille, coïncidant avec la limite des blocs lithosphériques. La diffusion dans le nord s'étend jusqu'à 60 km de profondeur, suggérant une zone de cisaillement qui pénètre dans le manteau supérieur sous la branche Nord. L'image obtenue révèle également la structure profonde des principaux blocs géologiques.

Dans la dernière partie de cette thèse, nous proposons une approche qui ouvre la voie à une tomographie passive 3D de la vitesse des ondes de volume. Enfin, toutes les applications présentées confirment l'efficacité de la matrice de réflexion pour révéler la structure du sous-sol. Elle fournit de nouvelles informations sur la distribution des hétérogénéités dans la croûte terrestre. Elle peut être appliquée sur n'importe quelle bande de fréquence, pourvu que l'échantillonnage spatial du réseau de géophones satisfasse au critère de Nyquist.

Mots-clés: diffusion et diffraction des ondes, structure de zone de faille, imagerie passive, matrix de réflexion, ondes de volume. v

Introduction

Natural disasters, in particular earthquakes, have devastating impacts on human lives.

Large earthquakes can result in hundreds of thousands fatalities, in a matter of seconds.

For that reason, fault zones, that are the source of earthquakes, have been the subject of many studies. Characterizing the structure of fault zones and understanding the physical processes inside the Earth is crucial for understanding the earthquakes physics and assessing seismic hazards.

Several geophysical methods were developed in order to image the structure of fault zones. Seismic methods, in particular migration techniques, that study the propagation of seismic waves inside the Earth, are the most common techniques to image the subsurface.

These techniques analyse the reflected and diffracted body waves that return back to the surface after hitting a geological interface or an object having a different velocity and/or density values. To relocate reflectors in depth, migration techniques rely on a correct velocity model. In complex media, and especially in fault zones, the velocity distribution is often not well defined because of the small-scale velocity variations induced by the presence of strong damage (low velocity zones) and the diversity of geological units. As a consequence, locating reflectors in such area becomes much more challenging. Since migration techniques need an accurate seismic velocity model, their application is thus limited in fault zones. Another limitation to these methods is related to the computational cost needed for their implementation.

Faults are generally detected from the discontinuity of interfaces across the slip zone.

Diffractions generated from small-scale heterogeneities, or scatterers, such as fractures and the edge of discontinuous layers, have been disregarded for a long time. Imaging those features is essential for characterizing faults [START_REF] Krey | The significance of diffraction in the investigation of faults[END_REF] and enhancing the resolution [START_REF] Khaidukov | Diffraction imaging by focusing-defocusing: An outlook on seismic superresolution[END_REF]. In this thesis, we aim to image the scattering around fault zones in depth. To do so, inspired by previous works in optics and ultrasound imaging [START_REF] Badon | Approche matricielle de l'imagerie optique des milieux diffusants[END_REF][START_REF] Barolle | Approche matricielle de la tomographie à cohérence optique[END_REF]Lambert et al., 2020a,b) and in the continuity of the work by Blondel [START_REF] Blondel | Matrix approach of seismic imaging[END_REF], we develop a matrix approach suitable for imaging fault zones at different scales and that shows two main advantages. First, this approach only requires an approximate velocity model of the medium. Second, it relies on passive measurements of ambient noise vibrations.

It is now well known that the Green's function (i.e the impulse response) between two points can be obtained by cross-correlating ambient noise recorded at these two points [START_REF] Weaver | Ultrasonics without a source: Thermal fluctuation correlations at mhz frequencies[END_REF][START_REF] Lobkis | On the emergence of the green's function in the correlations of a diffuse field[END_REF]Derode et al., 2003b,a).

This key property has revolutionized the field of seismology since active sources or earthquakes are no longer required for imaging. In this thesis, this approach is leveraged to passively record the response matrix associated with a network of geophones.

Since it contains the set of impulse responses between each geophone, this response matrix holds valuable information on the medium's structure and its mechanical properties. It constitutes the building block of the matrix imaging approach developed in this thesis. In particular, depending on the problem that is considered, it can be projected into different bases. For imaging purpose, a projection in a virtual focused basis is, for instance, particularly relevant. By applying appropriate time delays, i.e focusing operations, to the measured response matrix, a reflection matrix R is synthesized between virtual sources and receivers mapping the subsurface [START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]Lambert et al., 2020a). From this matrix, one can build the confocal image of the medium. One can also quantify the focusing quality and investigate the phase distortions induced by the mismatch between the actual wave velocity distribution and the considered model [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF][START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]Lambert et al., 2020bLambert et al., , 2021b)). This information is crucial as phase distortions can drastically degrade the seismic image and lead to its misinterpretation.

These phase distortions (also referred to as aberrations in the literature), are handled in optics and acoustics through adaptive focusing strategies [START_REF] Freeman | Deformable mirrors for all seasons and reasons[END_REF] image. To improve the focusing process, the focused reflection matrix is projected at its input or output in a dual basis (the Fourier basis) in which wave distortions can be isolated. Such distorted wave-fields form the so-called distortion matrix. Its singular value decomposition enables to decompose the field-of-view into several isoplanatic patches that correspond to the main scatterers contributing to backscattered waves at the considered travel time. The resulting eigenvectors provide the focusing laws that are required to fully correct for aberrations. A 4 km depth 3D image of the fault is obtained, with a resolution that is almost one eighth below the diffraction limit. In the second part of this chapter, a detailed interpretation of the fault zone is presented. The scattering distribution shows differences between the northwest and the southeast parts of the study region. The scattering is localized in the southeast and suggests a geometrical complexity of the fault zone in that area.

The third chapter is dedicated to a large scale application aiming to image the North Anatolian fault zone. Data from the 73-stations DANA (Dense Array for Northern Anatolia) network are used. Horizontal [0.1 0.5] Hz cross-correlations are considered and a layered S-waves velocity model is integrated. The distortion matrix is analysed locally since the medium presents a random and continuous reflectivity. A 60 km depth 3D image is obtained and reveals the scattering structure around this major fault. The scattering is much more localized around the northern strand than around the southern strand. Differences in the thickness of the crust are reported, with a step of the Moho below the northern branch. Sub-Moho structures are also identified, mainly in the North. These findings suggest the existence of a shear zone in the North that reaches the upper mantle.

The deep structure of the main geological blocks is also analysed.

Finally, in chapter IV, inspired by Lambert [START_REF] Lambert | Matrix approach for ultrasound imaging and quantification[END_REF], we demonstrate that axial aberrations can be probed by investigating the focused reflection matrix beyond the ballistic time. Its time dependence can be exploited to estimate the integrated wave velocity inside the medium. From the integrated velocity, the local velocity can potentially be retrieved via an inversion process. Matrix imaging thus paves the way towards a passive tomography of seismic waves. Such large earthquakes result from slip occurring on major fault zones. Major faults, with large accumulated displacement, are structures with complex geometry and heterogeneous composition. The geometry of such systems have an important impact on the fluid flow in the crust [START_REF] Knipe | Faulting processes and fault seal[END_REF] and the future rupture behaviour. Studies also show that the frequency-size and temporal statistics of earthquakes vary with the evolution and the maturity of faults [START_REF] Ben-Zion | Stress, slip, and earthquakes in models of complex single-fault systems incorporating brittle and creep deformations[END_REF][START_REF] Stirling | Fault trace complexity, cumulative slip, and the shape of the magnitude-frequency distribution for strike-slip faults: A global survey[END_REF]. That is why, characterizing the structure of fault zones is crucial to understand the earthquakes physics and provides important information on fault mechanics, seismic and hydraulic properties [START_REF] Caine | Fault zone architecture and fluid flow: Insights from field data and numerical modeling[END_REF][START_REF] Ben-Zion | Characterization of fault zones[END_REF][START_REF] Faulkner | On the structure and mechanical properties of large strike-slip faults[END_REF].

Faults are generally well defined at the surface by a series of traces that correspond to the regions where the movement takes place and are characterized by a structural discontinuity of rocks. However, the subsurface structure of faults is much more complex than predicted by their surface expression. Many disciplines focus on studying the structural complexity of faults in the crust from different perspectives and scales. Fault zones are often described as 3D volumes [START_REF] Riley | Controls on fault damage zone width, structure, and symmetry in the bandelier tuff, new mexico[END_REF] and consist of a fault core and a damage zone [START_REF] Caine | Fault zone architecture and permeability structure[END_REF]. The fault core is a region where most of the displacement is accommodated during the main earthquakes. The damage zone surrounds the fault core and corresponds to a network of fractures and deformation bands [START_REF] Chester | Internal structure and weakening mechanisms of the san andreas fault[END_REF][START_REF] Faulkner | On the internal structure and mechanics of large strike-slip fault zones: field observations of the carboneras fault in southeastern spain[END_REF]. The term damage zone has two slightly different meanings in structural geology [START_REF] Kim | Inferred fluid flow through fault damage zones based on the observation of stalactites in carbonate caves[END_REF][START_REF] Choi | Definition and classification of fault damage zones: A review and a new methodological approach[END_REF]. The first one refer to a highly deformed zone around the fault core [START_REF] Ben-Zion | San andreas fault zone head waves near parkfield, california[END_REF][START_REF] Caine | Fault zone architecture and permeability structure[END_REF], and the other usage is for a wider zone of structures around fault such as fault tips ans step-overs [START_REF] Mcgrath | Damage zone geometry around fault tips[END_REF][START_REF] Flodin | Evolution of a strike-slip fault network, valley of fire state park, southern nevada[END_REF]. In this manuscript, we mainly use the second definition of damage zone.

The damage of the fault zone decays in intensity with distance from the fault [START_REF] Chester | Implications for mechanical properties of brittle faults from observations of the punchbowl fault zone, california[END_REF]. It is pronounced in the first few kilometers of the crust, while at large depth, the high confining pressure suppresses the damage generation and contributes to its healing [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF][START_REF] Ben-Zion | Dynamic rupture on a material interface with spontaneous generation of plastic strain in the bulk[END_REF]. The damage distribution has considerable effects on fault behavior [START_REF] Dor | Geological observations of damage asymmetry in the structure of the san jacinto, san andreas and punchbowl faults in southern cali-fornia: A possible indicator for preferred rupture propagation direction[END_REF] and is the topic of many research studies [START_REF] Vermilye | The process zone: A microstructural view of fault growth[END_REF][START_REF] Joussineau | The evolution of the damage zone with fault growth in sandstone and its multiscale characteristics[END_REF][START_REF] Mitchell | The nature and origin of off-fault damage surrounding strikeslip fault zones with a wide range of displacements: A field study from the atacama fault system, northern chile[END_REF][START_REF] Faulkner | Scaling of fault damage zones with displacement and the implications for fault growth processes[END_REF][START_REF] Savage | Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands[END_REF]. Faults zone are generally associated with low-velocity zones (LVZ) [START_REF] Hough | Fault-zone waves observed at the southern joshua tree earthquake rupture zone[END_REF] that indicate the presence of fractured rocks. It was shown that these zones have a direct influence on rupture dynamics and the amplification of ground motion [START_REF] Harris | Effects of a low-velocity zone on a dynamic rupture[END_REF][START_REF] Li | Low-velocity fault-zone guided waves: numerical investigations of trapping efficiency[END_REF][START_REF] Qiu | Internal structure of the san jacinto fault zone at jackass flat from data recorded by a dense linear array[END_REF].

A number of geophysical methods aim to image and detect the structure of fault zones (see, e.g., review by [START_REF] Tanner | Fault detection[END_REF]. These methods study the properties of rocks such as density, resistivity and seismic velocity. They tend to differ in terms of acquisition and computational costs, and resolution. In this section, we briefly review some of these methods and we discuss their advantages and possible limitations.

I.1.1 Ground-penetration radar (GPR)

This method uses electromagnetic waves emitted at the surface. Echoes are reflected if the electromagnetic wave impedance changes when it hits an interface or an object.

In general, faults are imaged indirectly from offsets of layers across the fault [START_REF] Demanet | The use of geophysical prospecting for imaging active faults in the roer graben, belgium[END_REF]Brandes et al., 2018). GPR enables to characterize fault zones in 3-D with a high spatial resolution in the first few meters of the subsurface. The 3D representation is obtained by combining several 2D GPR lines [START_REF] Grasmueck | 3-d ground-penetrating radar applied to fracture imaging in gneiss[END_REF][START_REF] Christie | Assessing fault displacement and offfault deformation in an extensional tectonic setting using 3-d ground-penetrating radar imaging[END_REF]Brandes et al., 2018). A disadvantage of GPR method is the narrow profile length of the GPR lines that limits the lateral extension and the penetration depth of the imaged area.

Also, applying such method is costly and time consuming.

I.1.2 Electrical resistivity tomography (ERT)

The ERT method allows to obtain cross-sectional (2-D) plot of resistivity enabling to characterize the geological units in the near surface. It can work on a wide range of scale. A large scale application provides a profile of 20 km length and 4 km penetration depth [START_REF] Storz | Electrical resistivity tomography to investigate geological structures of the earth's upper crust[END_REF]. The interpretation of ERT profiles complement other imaging methods, and it is particularly useful to detect the exact location of very shallow faults [START_REF] Demanet | The use of geophysical prospecting for imaging active faults in the roer graben, belgium[END_REF][START_REF] Nguyen | Image processing of 2d resistivity data for imaging faults[END_REF][START_REF] Drahor | Integrated geophysical investigations in a fault zone located on southwestern part of izmir city, western anatolia, turkey[END_REF].

I.1.3 Seismic methods

The characterization of subsurface structures, and especially faults, is done by using the elastic waves that travel into the medium and carry information on the medium's properties. Such waves are generated by seismic active sources or by natural passive sources.

In this section, we briefly discuss active and passive seismic imaging techniques used to image the structural properties of faults zones.

I.1.3.1 Active seismic imaging

Seismic imaging with active sources is the most often used method for imaging faults.

Active seismic techniques are the most reliable techniques to characterize the crust and upper mantle. They are based on the analysis of elastic waves generated by artificially controlled sources such as explosives, vibroseis, airguns, sledgehammer, etc. The waves probe the medium and are then recorded at the surface at known geophones' positions.

The choice of receivers depends on the source used. Their number can vary from tens to thousands depending on the imaging application. Nowadays, the use of large scale 3D arrays of geophones becomes common, targeting the deep crust and the upper mantle [START_REF] Biondi | 3D seismic imaging[END_REF].

The recorded Earth's response should be well sampled at the surface to obtain a detailed structural image that scales from meters to kilometers deep in the crust. These methods are divided into two main categories: refraction and reflection methods. The latter are the most widely used in seismic investigations.

Seismic reflection techniques rely on the analysis of the reflected waves that return to the surface after hitting an interface or object having a different seismic impedance, i.e velocity and density values. These methods use seismic travel time information to build an image of the medium. Note that in a homogeneous medium, the reflector is supposed to be located at depth z = ct/2, with t the two-way travel time and c the wave velocity in the medium. Reflection methods are well suited to image horizontal or gently dipping structures and usually fail to image steeply inclined and vertical structures such as fault zones. Faults are then imaged by the discontinuity of reflectors.

In media with a heterogeneous distribution of velocity and complex wave trajectories, the time depth relation does not hold. When sharp and rapid velocity variations are present, the images of reflectors are distorted leading to potential misinterpretations. A post-processing of seismic data is then required in order to relocate reflectors at their correct position. This processing is know as seismic migration and is particularly useful for imaging dipping reflectors.

I.1.3.1.1 Conventional migration techniques

The basic goal of migration is to produce a structural image of the medium by focusing (migrating) the data inside the medium using the wave equation to recover the reflectivity properties of the Earth. The focusing operation will be further discussed in section I.3.

A large number of migration techniques have been developed through the years, with the improvements of computing power, to address a large variety of imaging problems. These techniques need an accurate velocity estimation of the medium. A review on the migration techniques and their limitations can be found in [START_REF] Gray | Seismic migration problems and solutions[END_REF] and [START_REF] Etgen | An overview of depth imaging in exploration geophysics[END_REF] and the references therein. These methods are divided into two main categories:

integral or ray-based methods [START_REF] French | Computer migration of oblique seismic reflection profiles[END_REF][START_REF] Schneider | Integral formulation for migration in two and three dimensions[END_REF], and wave-equation-based migration methods. Integral methods correspond to Kirchhoff migration and beam migration [START_REF] Hill | Prestack gaussian-beam depth migration[END_REF]. Wave-equation-based migration includes many methods such as reverse time migration (RTM) [START_REF] Baysal | Reverse time migration[END_REF][START_REF] Whitmore | Iterative depth migration by backward time propagation[END_REF][START_REF] Mcmechan | Migration by extrapolation of time-dependent boundary values[END_REF] and full waveform inversion (FWI) [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF][START_REF] Zhang | Passive seismic full waveform inversion using reconstructed body-waves for subsurface velocity construction[END_REF]). An essential step to seismic processing is stacking which corresponds to the process of averaging the seismic traces in order to reduce noise and improve the data quality. Migrations can be done either after stacking (post-stack) [START_REF] Whitmore | Two-dimensional post-stack depth migration: a survey of methods[END_REF] or before stacking (prestack) [START_REF] Biondi | 3-d prestack migration of common-azimuth data[END_REF]. Post-stack migration is faster and has a lower computational cost than pre-stack migration, since in post-stack, migration is done only for one set of stacked data. Pre-stack migrations provide better image quality than post-stack migrations and is suitable in case of strong lateral velocity and large dipping structures.

Seismic imaging started with time migration with a basic goal to position reflectors laterally. With time, migration techniques have evolved with the possibility of positioning reflectors vertically in depth by applying depth migration [START_REF] Huang | Prestack depth migration with acoustic screen propagators[END_REF][START_REF] Muijs | Prestack depth migration of primary and surface-related multiple reflections: Part i-imaging[END_REF]. Depth migration constitutes a more powerful processing technique and gives more structurally correct images than time migration. However, depth migration is more sensitive to the velocity model in order to do a correct time-depth conversion.

Kirchhoff migration [START_REF] Schneider | Integral formulation for migration in two and three dimensions[END_REF][START_REF] Bevc | Imaging complex structures with semirecursive kirchhoff migration[END_REF] is the most popular imaging technique and allows to obtain a good image of the medium. It is mainly designed to account for primary arrivals and does not include multiple reflections. For more complex imaging problems, strong velocity contrasts are handled through RTM method that accounts for the whole wave-field, including multiples, diffractions, and evanescent waves and therefore leads to a much better resolution in the image. The main limitations of RTM is related to its high computational cost and the necessity of using a very accurate velocity model.

Another recent seismic imaging technique is Full Waveform Inversion (FWI) [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF][START_REF] Zhang | Passive seismic full waveform inversion using reconstructed body-waves for subsurface velocity construction[END_REF]. It is a non-linear inverse problem that provides high-resolution velocity models of the subsurface by minimizing the misfit between recorded data and the modeled data. It uses all types of waves, including reflections and multiples, and can resolve several model parameters other than the simple velocity model, such as density and anisotropy [START_REF] Operto | A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice[END_REF]. Inverting different parameters simultaneously is an intensive and challenging task and increases the computational cost of the model reconstruction. It also suffers from cross-talk resulting from the coupling effects between different parameters [START_REF] Pan | A summary of several challenges facing multi-parameter elastic full waveform inversion[END_REF].

I.1.3.1.2 Seismic diffraction

The methods described earlier tend to image the fault indirectly. Recently, an increasing attention has been given to image directly faults as a distribution of small scale objects using the diffracted component of the wave-field. Diffractions arise from features, or scatterers, that are smaller than the wavelength, such as fractures and the edge of discontinuous layers. When interacting with these objects that act as secondary sources, the wave-field is scattered in all directions and passes through different regions of the subsurface. The scattered waves carry valuable information on the features that generate them.

Although diffracted/ scattered waves have been recognized to be relevant in characterizing faults and fractures since the 50's [START_REF] Krey | The significance of diffraction in the investigation of faults[END_REF], they are often neglected in seismic processing and imaging. [START_REF] Khaidukov | Diffraction imaging by focusing-defocusing: An outlook on seismic superresolution[END_REF] even referred to diffractions as "the abandoned stepchildren of traditional seismic processing and imaging". While conventional reflection imaging techniques provide good estimation on the subsurface velocity and reveal the geometry of specular reflectors that are larger than the wavelength, diffrac-tions resolve features that are smaller than the wavelength and are essential for enhancing the resolution [START_REF] Landa | A method for detection of diffracted waves on common-offset sections[END_REF][START_REF] Neidell | Perceptions in seismic imaging part 2: Reflective and diffractive contributions to seismic imaging[END_REF]. They can also lead to a super-resolution imaging [START_REF] Khaidukov | Diffraction imaging by focusing-defocusing: An outlook on seismic superresolution[END_REF][START_REF] Moser | Diffraction imaging in depth[END_REF] with the possibility of imaging features beyond the diffraction limit. Therefore, the detailed structural information provided by diffractions are crucial for seismic interpretations [START_REF] Bansal | Diffraction enhancement in prestack seismic data[END_REF].

One of the main problems in diffraction imaging is the weakness of the diffracted echoes compared with the reflections [START_REF] Klem-Musatov | Theory of seismic diffractions[END_REF]. Diffractions were often considered as noise in seismic migration techniques, attenuated or suppressing during seismic processing [START_REF] Bansal | Diffraction enhancement in prestack seismic data[END_REF][START_REF] Moser | Diffraction imaging in depth[END_REF]. To benefit from diffractions in order to successfully achieve high-resolution imaging, a separation between the diffracted component and the stonger reflected component of the wave-field is needed. Once reflections are suppressed from the data, a focusing operation of diffractions is performed, providing an image of the scatterers distributed inside the medium.

Two different strategies are applied. The first approach consists in separating diffractions and reflections in the data domain before migration [START_REF] Bansal | Diffraction enhancement in prestack seismic data[END_REF][START_REF] Fomel | Poststack velocity analysis by separation and imaging of seismic diffractions[END_REF][START_REF] Schwarz | Accessing the diffracted wavefield by coherent subtraction[END_REF]. The second approach tends to separate them by suppressing specular reflections during the migration scheme [START_REF] Khaidukov | Diffraction imaging by focusing-defocusing: An outlook on seismic superresolution[END_REF][START_REF] Moser | Diffraction imaging in depth[END_REF][START_REF] Koulakov | Distribution of seismic velocities and attenuation in the crust beneath the north anatolian fault (turkey) from local earthquake tomography[END_REF]. Some techniques rely on coherence measurements to separate between reflections and diffractions. These methods were first based on performing a coherence analysis of the diffracted wave-field to extract it from the data [START_REF] Khaidukov | Diffraction imaging by focusing-defocusing: An outlook on seismic superresolution[END_REF][START_REF] Moser | Diffraction imaging in depth[END_REF][START_REF] Koulakov | Distribution of seismic velocities and attenuation in the crust beneath the north anatolian fault (turkey) from local earthquake tomography[END_REF]. Recent methods showed that a better extraction of diffraction is done by assessing the coherence of reflections instead, that are then subtracted from the data [START_REF] Schwarz | Accessing the diffracted wavefield by coherent subtraction[END_REF]Schwarz, 2019b).

Through a comparison between reflection and diffraction images, Schwarz (2019a) illustrated the importance of diffractions in resolving detailed localized features of fault systems and improving the contrast and sharpness of the image. The diffractions complement reflection imaging and provide additional insights into geological interpretation, in particular, for fault imaging [START_REF] Dell | Using seismic diffractions for assessment of tectonic overprint and fault interpretation[END_REF]. Diffraction imaging depends on the quality of focusing that is directly linked to the existence of an accurate migration velocity. [START_REF] Moser | Diffraction imaging in depth[END_REF] confirmed this notion by showing the implications of an incorrect velocity model on diffraction imaging. Also, analysing the diffracting component brings valuable information on the medium's velocity. Many studies used diffractions in order to provide a high-resolution estimation of wave velocity [START_REF] Sava | Wave-equation migration velocity analysis by focusing diffractions and reflections[END_REF][START_REF] Reshef | Post-stack velocity analysis in the dip-angle domain using diffractions[END_REF][START_REF] Decker | Diffraction imaging and time-migration velocity analysis using oriented velocity continuation[END_REF].

Even though diffractions contribute to high-resolution imaging, they remain largely unexplored. More attention should be given to the diffracted/ scattered wave-field that, due to its unique properties, was shown to be essential for the characterization, identification and interpretation of small-scale structures such as faults and fractures.

I.1.3.1.3 Disadvantages of active seismic methods

Despite their success in characterizing subsurface structures, migration-and diffractionbased imaging techniques suffer from two main problems that limit their application.

First, they necessitate information on the medium's properties and a detailed knowledge on the velocity distribution. Yet, in several geological media, the velocity is not perfectly known. As we have seen earlier, fault zones are associated with complex velocity distribution, with low velocity damaged regions. Since migration needs a correct velocity model, conventional techniques fail in imaging the structural complexity of faults.

The second limitation is related to the computational cost of migration algorithm.

Seismic migration techniques require processing a large amount of data increasing the memory demand. Alternative imaging techniques that reduce the computational efforts should be adopted, without sacrificing the resolution of the images.

The described techniques are based on the characterization of the subsurface by using the seismic waves generated by active sources. Seismic exploration surveys based on active sources are not always favorable for environmental reasons and can be very costly. Also, the use of some active sources, such as explosions, is considered potentially dangerous and is restricted near populated cities.

I.1.3.2 Passive seismic imaging

Passive seismic imaging techniques have been adopted as an alternative to active seismic imaging. They are based on the analysis of ambient noise that is recorded at the surface of the Earth. The ambient noise refers to the permanent agitation of the ground surface due to human activities (machinery, cars, trains) or natural forces (wind, ocean waves). The vibrations propagate into the subsurface as seismic waves and can be used for imaging and monitoring.

Noise sources are found everywhere, all day long and all year around, providing a great opportunity to constantly record ambient noise at different locations on the Earth.

With the rapid development of seismological devices and processing techniques, large nodal arrays are deployed to record noise vibrations. Huge amount of noise records are nowadays available, offering the possibility for a wide variety of new applications for geophysical imaging. The main advantage of passive seismology is that active sources are no longer required, which reduces the cost and the environmental impact of geophysical surveys.

Most of applications are based on computing the cross-correlation function between noise records obtained at two stations to retrieve the Green's function between these two stations. The convergence of the cross-correlation function towards the Green's function is ensured for a homogeneous distribution of noise sources [START_REF] Froment | Estimation of the effect of nonisotropically distributed energy on the apparent arrival time in correlations[END_REF]. In section I.2.2, we briefly discuss the principles of Green's function retrieval from noise correlations. Cross-correlation methods have been used in the past few years for numerous applications to image and monitor active faults. Since seismic noise is dominated by surface waves, their contribution to the Green's function is generally easy to extract from noise cross-correlations [START_REF] Shapiro | Emergence of broadband rayleigh waves from correlations of the ambient seismic noise[END_REF]. Studies also reported the extraction of the body wave reflections from ambient noise (Roux et al., 2005a;[START_REF] Draganov | Reflection images from ambient seismic noise[END_REF]Poli et al., 2012a), which is essential for deep structure imaging. The imaging techniques are therefore divided in two main categories: surface-wave tomography and body wave reflection imaging. More details on ambient noise cross-correlations and applications can be found in [START_REF] Nakata | Seismic ambient noise[END_REF]. In the following we discuss the basics of these [START_REF] Shapiro | High-resolution surfacewave tomography from ambient seismic noise[END_REF]Sabra et al., 2005a;[START_REF] Stehly | Tomography of the alpine region from observations of seismic ambient noise[END_REF][START_REF] Yao | Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to se tibet[END_REF][START_REF] Ritzwoller | Ambient noise tomography with a large seismic array[END_REF]. It was also applied at smaller scale to image the complex medium below volcanoes (Brenguier et al., 2007;[START_REF] Stankiewicz | Lake toba volcano magma chamber imaged by ambient seismic noise tomography[END_REF].

Ambient noise surface wave tomography has been also extensively applied to study fault zones. The velocity distribution around faults gives direct information on the damage structure that has relevant influence on the seismic ground motion and crustal fluid flow.

The resolution and the depth of ambient noise tomographic images strongly depend on the distance between stations and the frequency content of the seismic waves [START_REF] Lin | Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band seismic array[END_REF]. Traditional tomographic techniques, that use seismic waves generated by earthquakes, are appropriate to image the crustal structure at depths below 2 km [START_REF] Allam | Seismic velocity structures in the southern California plateboundary environment from double-difference tomography[END_REF], the high frequency content of seismic noise opens the way to resolve near-surface features with a high-resolution. Ambient noise tomography is therefore essential to complement the large scale studies based on earthquake measurements. It was applied at large scale (hundreds of kilometers) to obtain the velocity structure of the crust around major fault systems, such as the San Andreas, San Jacinto and North Anatolian faults [START_REF] Zigone | Seismic tomography of the southern california plate boundary region from noise-based rayleigh and love waves[END_REF][START_REF] Barak | San andreas fault dip, peninsular ranges mafic lower crust and partial melt in the salton trough, southern california, from ambient-noise tomography[END_REF][START_REF] Warren | Crustal velocity structure of central and eastern turkey from ambient noise tomography[END_REF], but it can also be used at smaller scale (few kilometers) [START_REF] Roux | The san andreas fault revisited through seismicnoise and surface-wave tomography[END_REF]. By taking advantage of arrays with high density station, ambient noise tomography allows to resolve high-resolution shear wave velocity in the top 500 m of the crust [START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF][START_REF] Mordret | Shallow three-dimensional structure of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic array[END_REF][START_REF] Zigone | Imaging subsurface structures in the san jacinto fault zone with high-frequency noise recorded by dense linear arrays[END_REF]. Low velocity zones were identified and shown to correspond to the damage zone near the surface traces, and to sedimentary basins [START_REF] Hillers | Focal spot imaging based on zero lag cross-correlation amplitude fields: Application to dense array data at the San Jacinto fault zone[END_REF].

I.1.3.2.2 Body wave imaging

Even though ambient noise is dominated by surface waves, body waves can also be retrieved from noise records (Roux et al., 2005a;[START_REF] Draganov | Retrieval of reflections from seismic background-noise measurements[END_REF]. The extraction of body waves from noise cross-correlations is challenging because of their weak amplitudes compared to surface waves. Since body waves are sensitive to velocity fluctuations, they allow to construct high-resolution images of deep targets (Poli et al., 2012b) and provide a great opportunity to reveal the deep structure of the crust and the mantle around major faults. Although they offer a better vertical resolution than surface waves, body wave applications are still very limited in number.

Body waves were clearly observed in seismic noise at low and high frequencies. Low frequency body waves are mainly generated by the ocean activity (Poli et al., 2012b) while high frequency body waves (> 1 Hz) are mainly related to cultural activity such as vehicle traffic (Brenguier et al., 2019). The extraction of body waves has been reported at different scales. Global scale reflections give insight on the Earth's interior (Boué et al., 2013;[START_REF] Poli | Imaging the d reflector with noise correlations[END_REF], regional scale applications target the upper mantle and deep interfaces such as the transition noise (410 km and 660 km deep) (Poli et al., 2012b), while local scale applications aim to better understand the structure of the crust and image shallow interfaces [START_REF] Draganov | Reflection images from ambient seismic noise[END_REF][START_REF] Draganov | Seismic exploration-scale velocities and structure from ambient seismic noise (¿ 1 hz)[END_REF][START_REF] Nakata | Body wave extraction and tomography at long beach, california, with ambient-noise interferometry[END_REF]. [START_REF] Nakata | Body wave extraction and tomography at long beach, california, with ambient-noise interferometry[END_REF] performed the first high-resolution tomographic images from ambient noise P-waves recorded at the surface. It was also shown that reflectivity images recovered from ambient noise body waves are comparable with those created with active sources (Fig. I.2) [START_REF] Draganov | Reflection images from ambient seismic noise[END_REF][START_REF] Panea | Case study of retrieval of body-wave reflections using ambientnoise seismic interferometry with extra noise recordings[END_REF]. This implies that ambient noise body waves constitute a reliable technique to image the crust, thereby providing an alternative to active seismic imaging methods. We will see in the following chapters that ambient noise body waves can be used to build the reflectivity images around major fault zones. Such images will reveal to retrieve the scattering structure of fault zone with an unprecedented resolution.

I.1.3.2.3 Monitoring

Since ambient noise measurements can be repeated, they can be used to continuously monitor the crust. The Green's function retrieved from noise cross-correlations constitutes a low cost reliable tool to study the temporal and spatial variations of the elastic properties of the subsurface. The principle is to monitor the velocity variations of complex geological structures such as volcanoes [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at merapi volcano, indonesia[END_REF] and active fault zones [START_REF] Chen | Distribution of seismic wave speed changes associated with the 12 may 2008 mw 7.9 wenchuan earthquake[END_REF]. Reservoir monitoring and environmental changes observations [START_REF] Mordret | Monitoring southwest greenland's ice sheet melt with ambient seismic noise[END_REF] are also among the seismic noise-based applications.

Tectonic events and volcanic eruptions are associated with stress and pore-pressure perturbations. These perturbations result in changes in the seismic velocity that appear as time shifts in the measured waveforms [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: An application to the calaveras fault, california[END_REF]. The time shifts are detected from the ambient noise cross-correlation functions [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at merapi volcano, indonesia[END_REF]Brenguier et al., 2008a).

Coda-based approach is the standard technique adopted for monitoring changes in the crust. Seismic coda, that corresponds to the scattered waves that arrive after the ballistic waves, have less sensitivity to noise source changes than ballistic waves [START_REF] Colombi | On the temporal stability of the coda of ambient noise correlations[END_REF] and is useful to track weak small changes of the seismic velocity (Brenguier et al., 2008a). Recently, a new approach is proposed to monitor temporal changes using ballistic waves reconstructed from noise correlations [START_REF] Mordret | Noise-based ballistic wave passive seismic monitoring-part 2: surface waves[END_REF][START_REF] Mordret | Noise-based ballistic wave passive seismic monitoring-part 2: surface waves[END_REF].

An instantaneous drop of the seismic velocity has been observed after large earthquakes [START_REF] Peng | Temporal changes of shallow seismic velocity around the karadere-düzce branch of the north anatolian fault and strong ground motion[END_REF]. Monitoring fault zones is essential to provide important information on the earthquake cycle and the healing process of the damage at depth.

I.2 Seismic ambient noise

Seismic imaging is based on the study of seismic waves that propagate inside the Earth and that hold important information on the medium's properties. Seismic waves are mainly generated by earthquakes or by controlled active sources. However, earthquakes are distributed in time and are limited to few regions. On the other hand, the use of active sources is prohibited in urban areas and can be extremely costly. An alternative is to use ambient noise vibrations to study the Earth's response continuously and at every location.

Seismic noise correlations constitute a reliable approach to estimate the Green's function between two receivers at know positions. This approach, also known as seismic interferometry, was proved successful for developing numerous passive imaging techniques. As we have seen in the previous section, the number of noise based applications has tremendously increased in the last two decades. In the following, we give a brief review on the origin of ambient noise and discuss the method of reconstructing the Green's function from the cross-correlations of noise records.

I.2.1 Noise origin

The Earth's vibrations were observed on the first seismograms and were historically known as 'microseisms'. Later on, authors started to categorize noise based on its origin. 'Microseisms' term is then used to refer to noise generated by natural sources while 'microtremors' to characterize vibrations caused by human activities. More details on the origin of noise sources and the historical overview of the use of ambient noise can be found in [START_REF] Bonnefoy-Claudet | The nature of noise wavefield and its applications for site effects studies: A literature review[END_REF].

The origin of ambient noise differs with the frequency band [START_REF] Gutenberg | Microseisms[END_REF][START_REF] Asten | Geological control on the three-component spectra of rayleigh-wave microseisms[END_REF]. The low frequency noise (below 1 Hz) is generally caused by natural sources.

The energy content of low frequency noise is often divided in three main categories. The 'hum', that is found at long periods (below 30 s), originates from infragravity waves [START_REF] Kobayashi | Continuous excitation of planetary free oscillations by atmospheric disturbances[END_REF][START_REF] Suda | Earth's background free oscillations[END_REF]. The primary microseism peak is observed at period between 10 and 20 s and founds its origin from the interaction between ocean waves and the bottom of the ocean [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF]. The secondary microseism peak is generally around 5s. It corresponds to the strongest peak in the noise spectrum and results from the interference of swells travelling in opposite directions [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF]. The study of low frequency (< 1 Hz) seismic noise shows clear seasonal variability related with the variability of oceanic activity between winter an summer [START_REF] Stehly | A study of the seismic noise from its long-range correlation properties[END_REF]. The high frequency noise (above 1 Hz) is associated with the wide variety of human activities, including airplanes, trains and road traffic. The activities are different for different locations on the Earth. Noise level amplitudes are higher in populated urban cities. It also exhibits a temporal variation between day and night, weekdays and weekends [START_REF] Mcnamara | Ambient noise levels in the continental united states[END_REF]. Recently, and during the COVID-19 pandemic, the lockdown effects were assessed with a large number of surface and boreholes stations around the world [START_REF] Lecocq | Global quieting of highfrequency seismic noise due to covid-19 pandemic lockdown measures[END_REF]. A global quieting of the Earth has been reported in the high-frequency regime (> 1 Hz).

I.2.2 Noise cross-correlations and Green's function

Cross-correlation technique is the most used technique applied to seismic ambient noise.

Under certain conditions, noise cross correlation converges towards the Green's function between receiving stations [START_REF] Shapiro | Emergence of broadband rayleigh waves from correlations of the ambient seismic noise[END_REF]. The Green's function is defined as a solution of the wave-equation and corresponds to the signal that would be recorded by one station if a pulse source was placed at the other. Therefore, it gathers all the information on the surrounding medium and is essential for imaging purposes.

The Green's function retrieval via noise cross-correlation has been investigated in several fields of wave physics. Noise cross-correlation method was first applied to helioseismology to measure travel times of acoustic waves propagating in the Sun [START_REF] Duvall | Time-distance helioseismology[END_REF][START_REF] Giles | A subsurface flow of material from the sun's equator to its poles[END_REF]. Long before, the idea of using the ambient noise measurements to retrieve surface waves properties was proposed by [START_REF] Aki | Space and time spectra of stationary stochastic waves, with special reference to microtremors[END_REF].

With ultrasound laboratory experiments, [START_REF] Weaver | Ultrasonics without a source: Thermal fluctuation correlations at mhz frequencies[END_REF] and [START_REF] Lobkis | On the emergence of the green's function in the correlations of a diffuse field[END_REF] proved that the correlation of diffuse thermal noise recorded at two transducers attached at the surface of an aluminium sample yields the Green's function be-tween these two transducers. This result is valid for a randomization produced either by a distribution of sources or by a multiple scattering inside the medium. These observations were later explained theoretically through many demonstrations, both in closed and open media. Derivations of the Green's function from noise correlation were demonstrated by the fluctuation-dissipation theorem [START_REF] Callen | Irreversibility and generalized noise[END_REF][START_REF] Van Tiggelen | Green function retrieval and time reversal in a disordered world[END_REF][START_REF] Godin | Emergence of the acoustic green's function from thermal noise[END_REF], by the stationary phase theorem for scattered waves [START_REF] Snieder | Extracting the green's function from the correlation of coda waves: A derivation based on stationary phase[END_REF] and by the reciprocity theorem for elastic waves [START_REF] Wapenaar | Retrieving the elastodynamic green's function of an arbitrary inhomogeneous medium by cross correlation[END_REF]. These demonstrations are summarized in several review papers [START_REF] Campillo | Phase and correlation inrandom'seismic fields and the reconstruction of the green function[END_REF][START_REF] Weaver | Diffuse fields in ultrasonics and seismology[END_REF][START_REF] Larose | Correlation of random wavefields: An interdisciplinary review[END_REF][START_REF] Gouédard | Convergence of the two-point correlation function toward the green's function in the context of a seismic-prospecting data set[END_REF]Wapenaar et al., 2010a,b;[START_REF] Campillo | Seismic imaging and monitoring with ambient noise correlations[END_REF].

In the simple case of a homogeneous medium with attenuation and a uniform distribution of noise sources, Roux et al. (2005b) used the stationary phase theorem to demonstrate the relation between the temporal derivative of the correlation function C AB (t) between two points A and B and the causal and anti-causal (time-reversed) Green's function, G AB (t) and G AB (-t), respectively, between these two points, such as:

∂ ∂t C AB (t) ∝ G AB (t) -G AB (-t) (I.1)
The theorem was then generalised to heterogeneous media [START_REF] Weaver | Diffuse fields in open systems and the emergence of the green's function (l)[END_REF][START_REF] Snieder | Extracting the green's function of attenuating heterogeneous acoustic media from uncorrelated waves[END_REF] and to cavities [START_REF] Weaver | Ultrasonics without a source: Thermal fluctuation correlations at mhz frequencies[END_REF]. For more details, the reader can refer to the review of [START_REF] Campillo | Seismic imaging and monitoring with ambient noise correlations[END_REF]. In a more didactic way, Derode et al. (2003b,a) interpreted the reconstruction of the Green's function through an analogy between field-field correlation and a fictitious time reversal experiment [START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF].

The correlation of random noise was successfully applied in seismology [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF][START_REF] Schuster | Interferometric/daylight seismic imaging[END_REF], in acoustics [START_REF] Malcolm | Extracting the green function from diffuse, equipartitioned waves[END_REF][START_REF] Larose | Reconstruction of rayleigh-lamb dispersion spectrum based on noise obtained from an air-jet forcing[END_REF][START_REF] Sabra | Structural health monitoring by extraction of coherent guided waves from diffuse fields[END_REF], in medical imaging [START_REF] Gallot | Passive elastography: shear-wave tomography from physiological-noise correlation in soft tissues[END_REF], in oceanography [START_REF] Roux | Extracting coherent wave fronts from acoustic ambient noise in the ocean[END_REF]Sabra et al., 2005b), in seismic exploration [START_REF] Draganov | Retrieval of reflections from seismic background-noise measurements[END_REF] and even with electromagnetic waves whether it be at microwave [START_REF] Davy | Green's function retrieval and passive imaging from correlations of wideband thermal radiations[END_REF] or optical frequencies [START_REF] Badon | Retrieving time-dependent green's functions in optics with low-coherence interferometry[END_REF]. In seismology, [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF] were the first to demonstrate that the surface wave Green's function can be reconstructed by correlating multiply scattered coda waves recorded in Mexico. Shortly after, [START_REF] Paul | Empirical synthesis of time-asymmetrical green functions from the correlation of coda waves[END_REF] applied the same technique at the regional scale to seismic coda recorded in Alaska. They found that the causal and anticausal parts of the Green's function differ in amplitude. This was explained by the anisotropy of the diffuse field. Such asymmetry in the correlations can be used to study the distribution and the directivity of noise sources [START_REF] Stehly | A study of the seismic noise from its long-range correlation properties[END_REF].

Other than the seismic coda, seismic noise correlations were first computed by [START_REF] Shapiro | Emergence of broadband rayleigh waves from correlations of the ambient seismic noise[END_REF]. They were shown to converge to the Green's function of the surface waves. The retrieval of the Green's function from noise correlations is not always simple and specific processing of the data should be performed [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF][START_REF] Groos | Performance of different processing schemes in seismic noise cross-correlations[END_REF].

The successful reconstruction of the Green's function and the large amount of noise data led to numerous applications in seismology for imaging the crust at different scales.

At first, extensive works have been performed to extract surface wave information from noise correlation, because seismic noise mainly contains surface waves [START_REF] Shapiro | High-resolution surfacewave tomography from ambient seismic noise[END_REF]Sabra et al., 2005a;[START_REF] Zheng | Surface wave tomography of china from ambient seismic noise correlation[END_REF][START_REF] Stehly | Tomography of the alpine region from observations of seismic ambient noise[END_REF][START_REF] Yao | Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to se tibet[END_REF][START_REF] Nishida | Global surface wave tomography using seismic hum[END_REF][START_REF] Ritzwoller | Ambient noise tomography with a large seismic array[END_REF][START_REF] Wang | Ambient noise surface wave tomography of marginal seas in east asia[END_REF][START_REF] Lu | High-resolution surface wave tomography of the european crust and uppermost mantle from ambient seismic noise[END_REF].

Body waves can also be retrieved from noise correlations. The first identification of body waves from noise correlations was done by Roux et al. (2005a). Using a dense array in the Parkfield area, California, they demonstrated the presence of P-waves and Rayleigh waves in the correlations functions. Following these observations, studies showed the possibility of extracting body wave reflections from the deep Earth leading to a high resolution imaging of the crust and upper mantle (Poli et al., 2012a,b;[START_REF] Nishida | Global propagation of body waves revealed by cross-correlation analysis of seismic hum[END_REF]Boué et al., 2013;[START_REF] Lin | Extracting seismic core phases with array interferometry[END_REF][START_REF] Retailleau | Ambient seismic noise imaging of the lowermost mantle beneath the North Atlantic Ocean[END_REF]. Because seismic noise is omnipresent, it is also possible to use ambient noise cross-correlations to continuously monitor fault zones (Brenguier et al., 2008b(Brenguier et al., , 2019) ) and volcanoes [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at merapi volcano, indonesia[END_REF]Brenguier et al., 2008a). More details on the imaging and monitoring applications can be found in sections I.1.3.2.

I.2.3 Analogy with the time reversal experiment

To understand how one can obtain the Green's function between two points by crosscorrelating the wave-fields recorded at these two points, we refer to the simple heuristic approach of Derode et al. (2003b,a). This approach is based on an analogy between a cross-correlation product and a time reversal experiment [START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF].

First we consider N noise sources S i emitting uncorrelated white noises n i (t), such that:

n i (t) ⊗ n j (-t) = σ 2 δ(t)δ ij (I.
u A (t) = i n i (t) ⊗ G S i A (t) u B (t) = i n i (t) ⊗ G S i B (t) (I.3)
where the symbol ⊗ accounts for the correlation product. The cross-correlation C AB (τ )

between the received responses at A and B, can be expressed as:

C AB (τ ) = N i=1 N j=1 G S i A (-t) ⊗ n i (-t) ⊗ G S j B (t) ⊗ n j (t) (I.4)
where τ is the correlation time lag. G S i A describes the propagation from a source at S i to a receiver at A. Because of the spatial reciprocity, G S i A = G AS i . This means that the signal recorded by A when S i emits a pulse is identical to the signal recorded at S i following a pulse emitted at A. Using the spatial and temporal incoherence of noise sources (Eq. I.2), Eq. I.4 can thus be rewritten as:

C AB (τ ) = σ 2 N i=1 G AS i (-t) ⊗ G S i B (t) (I.
u B (t) = N i=1 G AS i (-t) ⊗ G S i B (t) (I.6)
The time-reversed signal is thus proportional to the cross-correlation of a random wave-

field
u B (t) = G AB (t) + G AB (-t) (I.7)
Combining Eq. I.6 and I.7 we obtain:

N i=1 G AS i (-t) ⊗ G S i B (t) ∝ G AB (t) + G AB (-t) (I.8)
This relation implies through a time-reversal analogy that the cross-correlation of a random wave-field measured at A and B is equal to the Green's function between A and B.

This fundamental result holds as long as the random wave-field results from an isotropic

S i A B G S A (t) i G S B (t) i S i A B G AS (t) i ∑G AS (-t) × G S B (t) i i i a b c t= S i A B t= t< G AB (-t) ∑G AS (-t) × G S B (t) = G AB (t)+G AB (-t) i i d S i A B t= t> C AB (τ)=G AS (-t) × G S B (t) i i

I.3 Migration techniques

In this section, we first describe the wave-field extrapolation technique that is an essential step for the migration process before tackling the redatuming operation. In a second step, the main limitations of these different methods are highlighted.

I.3.1 Wave-field extrapolation

Reflection imaging techniques aim to collect the response of the Earth at its surface and then use this response to build an image of the medium. When seismic waves are generated by a seismic source, they propagate into the subsurface where they are reflected by planar interfaces or specular objects and scattered by smaller-scale heterogeneities. The reflected and scattered waves are recorded by receivers, generally located at the surface of the Earth. If the recorded waves are projected backward into the subsurface, they will focus back at the position of the secondary source that reflected or scattered them. This operation is known as wave-field extrapolation [START_REF] Berkhout | Wave field extrapolation techniques in seismic migration, a tutorial[END_REF] and it is the principle of seismic migration techniques, no matter which method is used. In this section, we will briefly describe seismic extrapolation. More details can be found in [START_REF] Berkhout | Wave field extrapolation techniques in seismic migration, a tutorial[END_REF].

Wave field extrapolation is a process that uses the wave equation to transform the data recorded at a surface into a wave-field at virtual receivers position located in another plane, for instance inside the medium. Forward extrapolation is the operation that con-sists in going away from the source while backward extrapolation consists in converging towards the source.

Different methods can be used in wave-field extrapolation. First, and before dis- In this section, we will review the derivation of Kirchhoff and Rayleigh integrals for forward and backward extrapolation. The full derivation can be found in [START_REF] Berkhout | Wave field extrapolation techniques in seismic migration, a tutorial[END_REF].

cussing
We consider a closed surface S and a homogeneous medium with a wave velocity c.

Let us consider the Green's theorem [START_REF] Berkhout | Wave field extrapolation techniques in seismic migration, a tutorial[END_REF]:

ˆV (f ∆g -g∆f )dV = ˛S(f ∇g -g∇f ) • ndS (I.9)
We choose f = Φ(r, ω) the measured wave-field in the frequency domain at position r, denoted with the Cartesian coordinates r(x, y, z), that satisfies the source-free Helmholtz equation:

∇ 2 Φ + k 2 Φ = 0 (I.10)
where k = ω/c is the wavenumber. We define g = G(r, r A , ω) the causal Green's function at point A inside S, that is a solution of the Helmholtz equation:

∇ 2 G + k 2 G = -4πδ (r -r A ) (I.11) with r A (x A , y A , z A ).
For a homogeneous medium, the causal Green's function G is expressed by:

G (r A , r, ω) = e -ik||r A -r|| 4π||r A -r|| (I.12)
If we substitute Eq. I.10 and Eq. I.11 in the Green's theorem (Eq. I.9), the first version of Kirchhoff integral is obtained:

Φ(r A , ω) = - 1 4π ˛ Φ(r, ω) ∂ ∂n G (r A , r, ω) -G (r A , r, ω) ∂ ∂n Φ(r, ω) dS (I.13)
Kirchhoff's theory gives a mathematical formulation to the Huygens-Fresnel principle.

Eq. I.13 states that the knowledge of the wave-field and its normal derivative on a closed that is also a solution of Eq. I.11 (the symbol* stands for phase conjugate). Similarly, a second version of the Kirchhoff integral can be derived [START_REF] Wapenaar | Inverse extrapolation of primary seismic waves[END_REF]:

Φ(r A , ω) = - 1 4π ˛ Φ(r, ω) ∂ ∂n G * (r A , r, ω) -G * (r A , r, ω) ∂ ∂n Φ(r, ω) dS (I.14)
This second version is used for backward or inverse extrapolation and it constitutes the basic equation for image formation. When waves propagate inside the medium, they are reflected to the surface by scattering objects acting as secondary sources. The wave-field at the surface is recorded at time t. To build the reflectivity image of the medium, the wave-field needs to be calculated at depth at earlier times. We thus need to go backwards in time to converge towards the secondary source. In the frequency domain, a backward propagation is modelled by the anti-causal Green's function G * (r, r A , ω) (Fig. I.5a). Equivalent of time reversal in the frequency domain, phase conjugation can be seen a matched filter that compensates for the phase accumulated by the wave during its forward propagation between the secondary source and the measurement surface.

In seismic imaging, the closed surface condition does not hold because measurements cannot be done all around the Earth but they are carried locally along almost planar interfaces. Moreover, in practice, it is difficult to measure Φ(r, ω) and its normal derivative ∂φ/∂n at the same time. Consequently, Kirchhoff integrals are difficult to implement.

By setting boundary conditions for G(r, r A , ω) on planar acquisition surface S, such that G(r, r A , ω) = 0 and ∂ ∂n G (r A , r, ω) = 0, one can derive, respectively, from Eq. I.13 the Rayleigh integrals for forward-propagation (Fig. I.5b) defined as:

Φ(r A , ω) = - 1 2π ˆS Φ (r, ω) ∂ ∂n G (r A , r, ω) dS (I.15)
and

Φ(r A , ω) = 1 2π ˆS G (r A , r, ω) ∂ ∂n Φ(r, ω) dS (I.16)
For backward extrapolation, the Rayleigh integrals are given by:

Φ(r A , ω) = - 1 2π ˆS Φ (r, ω) ∂ ∂n G * (r A , r, ω) dS (I.17)
and

Φ(r A , ω) = 1 2π ˆS G * (r A , r, ω) ∂ ∂n Φ(r, ω) dS (I.18)
Applying Eqs. I.17 and I.18 consists in back-propagating waves towards a scattering object acting as secondary source inside the medium (Fig. I.5c). In other words, the wavefield produced by secondary sources is back-propagated and estimated at the virtual receiver position r A . This corresponds to the principle of numerical focusing at reception that enables the construction of an image of all the secondary sources (i.e scatterers and reflectors) inside the medium. a) Forward (Eq. I.13) and backward (Eq. I.14) extrapolation with Kirchhoff integral. Assuming a source outside the volume V , the wave-field and its normal derivative are recorded at surface S. The blue part of the surface corresponds to when the wave-field is entering into the volume, and the red part represents the surface where the wave-field travels outside the volume. G and G * denote the causal and anti-causal Green's functions, respectively. (b) Forward extrapolation using the Rayleigh integrals (Eqs. I.15 and I.16). The wavefield generated by a source at the surface and scattered by the heterogeneities of the medium is recorded on a planar surface S. The causal Green's function G is used to obtain the wave-field at point A(r A ). (c) Backward extrapolation using the Rayleigh integrals (Eqs. I.17 and I.18). A secondary source generates a wave-field that propagates inside the medium. The wave-field produced by the secondary source is back-propagated and estimated at the virtual receiver A using the anti-causal Green's function G * .

We have seen how a wave-field recorded at the Earth's surface can be back-propagated inside a homogeneous medium. In the next chapters, Eqs. I.17 and I.18 will be used to retrieve the reflectivity of the medium. For heterogeneous media, the Green's function is more complicated and Eq. I.12 will not hold anymore.

I.3.1.2 Extrapolation in the wavenumber-frequency domain

In the previous section, we have seen how the wave-field can be extrapolated in the spacefrequency domain. In this section we will represent the extrapolation in the wavenumberfrequency domain. We will focus only on the backward extrapolation since it enables to build the image of the medium.

We consider again the wave equation in Eq. I.10 such as:

∂ 2 Φ ∂x 2 + ∂ 2 Φ ∂y 2 + ∂ 2 Φ ∂z 2 + k 2 Φ = 0 (I.19)
Its Fourier transform with respect to x and y is expressed by:

∂ 2 Φ ∂z 2 + k 2 -k 2 x -k 2 y Φ = 0 (I.20)
where Φ = Φ(k x , k y , z, ω) corresponds to the wave-field decomposed into a set of plane waves. The Fourier coordinates k x , k y and k z are linked to the frequency ω through the dispersion relation:

k 2 x + k 2 y + k 2 z = ω 2 c 2 (I.21)
By integrating Eq. I.23, the solution of Eq. I.21 can be written as:

Φ (k x , k y , z A ) = Φ (k x , k y , z) e -i ( ω c ) 2 -(k 2 x +k 2 y )∆z , for k 2 x + k 2 y ≤ ω 2 /c 2 (I.22)
and

Φ (k x , k y , z A ) = Φ (k x , k y , z) e -( ω c ) 2 -(k 2 x +k 2 y )∆z , for k 2 x + k 2 y > ω 2 /c 2 (I.23)
where ∆z = z A -z.

Eqs. I.22 and I.23 show that the wave-field extrapolation can be done by considering plane waves that are easily extrapolated in the plane wave basis by a simple multiplication of the wave-field by the phase term exp (ik z ∆z). This method is an adequate tool in a layered medium with a velocity variation in the z direction by applying a phase shift corresponding to the wave velocity in each layer. However, it cannot account for lateral velocity variations.

The aforementioned extrapolation methods were derived for homogeneous and layered media, but small scale-velocity variations were not handled. These variations, if not taken into account in the extrapolation operation, will lead to strong distortions of the wave-field and a misfocusing of waves. That's why, migration always requires an estimation of the velocity, prior to the imaging process, in order to obtain accurate images of the subsurface.

I.3.2 Active seismic imaging: Redatuming

In seismic investigations, as we have seen in section I.1.3.1, data acquisition is done by using an active source and recording the wave-field at several receivers' positions along one single line or forming a 2D array. The survey can be repeated several times by changing the position of the source. The measurements obtained from the set of sources (emitters) and the set of receivers (detectors) can be organized as a 2D data matrix [START_REF] Berkhout | Pushing the limits of seismic imaging, part i: Prestack migration in terms of double dynamic focusing[END_REF], the so-called reflection matrix. The same principle can be applied to passive measurements described in section I.2, where, through ambient noise cross-correlations, the Green's function between a virtual source and a receiver can be retrieved. The crosscorrelation matrix of noise recorded by each station shall converge towards the reflection matrix associated with the network of geophones.

To obtain the image of the medium, migration aims at bringing waves into focus [START_REF] Berkhout | Wave field extrapolation techniques in seismic migration, a tutorial[END_REF]. Inspired by medical imaging and confocal microscopy techniques, migration can be formulated through the concept of focusing. In ultrasound imaging, a set of incident waves is also used to insonify the medium. Each reflected wave-field is recorded by an ). This operation is known as numerical focusing in reception. Therefore, DAS beamforming enables a simultaneous focusing at emission and reception for every point inside the medium, i.e. every pixel of the image. This is the principle of confocal imaging, where the focusing points in emission and detection are identical.

As we have seen in section I.3.1, time shifts need to be applied to the wave-field to compensate for propagation effects and focus inside the medium. This operation is thus analogous to DAS beamforming techniques used in ultrasound imaging and that we just described. In fact, numerical focusing is similar to the stacking operation that involves a move-out correction and a summation [START_REF] Berkhout | Pushing the limits of seismic imaging, part i: Prestack migration in terms of double dynamic focusing[END_REF]. Migration process can therefore be expressed as a double focusing operation done in the frequency domain via simple matrix multiplications [START_REF] Berkhout | Wave field extrapolation techniques in seismic migration, a tutorial[END_REF] ). This process is commonly referred to as "redatuming" [START_REF] Berkhout | A unified approach to acoustical reflection imaging. II: The inverse problem[END_REF][START_REF] Schuster | A theoretical overview of model-based and correlation-based redatuming methods[END_REF]: by applying time shifts to the data, the wave-field appears to be emitted and recorded by virtual sources and receivers inside the medium. The focusing operation at emission and reception is repeated for each point in the subsurface and a set of virtual sources and receivers is obtained (Fig.

I.7c).

The obtained responses can be grouped in a 2D matrix whose lines and columns correspond to the positions of the synthesized sources and receivers. This matrix will be referred to as the focused reflection matrix R F in this thesis. The diagonal elements of R F (co-located focusing points at emission and reception) provide the image of the medium. When the focus point at emission differs from the focus point at reception, the bifocal version of migration is obtained, yielding to the full reflection information of the medium. Since the goal of migration is to provide a structural image of the medium, only the diagonal elements of the reflection matrix are generally calculated. Yet, as we will show in this thesis, the off-diagonal coefficients of R F contain a lot of information on both aberrations and multiple scattering, that constitute the fundamental limits of reflection imaging methods.

I.3.3 Limits of seismic imaging

We have described earlier the principles of image formation under several hypothesis: a homogeneous or layered velocity model, and the single scattering assumption. However, in complex media, these different assumptions are no longer checked. In reality, and in seismic imaging, building a reliable estimator of the medium reflectivity then becomes a difficult task. In the crust and mantle, the velocity of waves shows a gradual increase with depth. In specific geological environments, such as faults and volcanoes, the distribution of velocity can not be illustrated by a simple stratified model. Important lateral velocity variations may be present as well as localized low velocity zones. Imaging such media is thus much more challenging.

I.3.3.1 Multiple scattering

In reflection seismic imaging, the travel times of waves between the surface and each focusing point are obtained under the single-scattering assumption hence assuming a homogeneous weakly scattering medium (Born approximation). The single scattering assumption is valid if the seismic waves are only reflected once by each heterogeneity of the medium (Fig. I.8a). In other words, the path length should be smaller than the scattering mean free path s , the typical distance between two scattering events. Assuming a homogeneous medium, the travel times of the back-scattered waves are then computed by measuring the distance between the object and the surface: t = 2z/c. However, when the travel time becomes of the order of mean free time τ S = l s /c, waves are reflected several times through their travel between the source and the receiver (Fig. I.8b). These events appear in the seismic records and are know as the multiple scattering contribution. This contribution is often seen as a nightmare for imaging. The corresponding arrivals can be misinterpreted and considered as single scattering arrival, resulting in an incoherent background on the image.

In most configurations, the multiple scattering contribution is considered as negligible compared to the single scattering contribution. However, in complex area where short-scale and strong fluctuations of the wave velocity arise, multiple scattering becomes predominant and methods shall be used to remove or attenuate this contribution. Migration stacking, as focused beamforming techniques, leads to enhance the single scattering contribution compared to the multiple scattering background. The singly-scattered echoes actually interfere constructively, while multiply-scattered waves are incoherent with each other. However, the penetration depth of conventional imaging techniques remains limited to a few s . Recently the matrix approach was successfully applied to seismology in a multiple scattering regime and enabled to separate single and multiple scattering contributions extending the imaging depth up to 10 s [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF].

I.3.3.2 Multiple reflections (multiples)

A common problem in seismic imaging is the presence of multiple reflections. This problem occurs in layered media with strong density or velocity contrasts. Instead of being reflected just once off a contrast, waves bounce multiple times between reflectors inside the Earth before being sent to the surface (Fig. I.8c).

In most imaging methods, waves are treated as if they are reflected just once inside the Earth. In reality, multiple reflections, i.e multiples, are also present. One can differentiate between surface-related multiples when the waves bounce between the Earth's surface and an interface between two geological layers, and internal multiples when waves do not reflect off the Earth's surface (Fig. I.8c). Surface and internal multiples are often neglected and not handled properly, leading to erroneous migration images and artefacts.

As a result, reflectors are placed at incorrect positions. Consequently, the true structure of the medium is blurred and the images misinterpreted.

In order to obtain accurate images of the subsurface, multiples have to be properly compensated when focusing operations are applied. Several methods were developed to handle multiples. Among these methods, the Marchenko method [START_REF] Wapenaar | Marchenko imaging[END_REF] and its many extensions [START_REF] Wapenaar | Single-sided marchenko focusing of compressional and shear waves[END_REF][START_REF] Slob | Theory for marchenko imaging of marine seismic data with free surface multiple elimination[END_REF][START_REF] Reinicke | Elastodynamic single-sided homogeneous green's function representation: Theory and numerical examples[END_REF] were developed and showed promising results. The Marchenko method aims to retrieve the upgoing and downgoing Green's functions between a virtual source in the subsurface and a receiver at the surface from the reflection response recorded at the surface. This is done by solving Marchenko equations to obtain the focusing functions that link the recorded wavefield to the Green's functions. The obtained Green's functions contain internal multiples and are used to obtain a redatumed reflection response that is free of internal multiples. To do so, multidimensional deconvolution (MDD) is applied to the retrieved upgoing and downgoing Green's functions [START_REF] Wapenaar | Marchenko imaging[END_REF][START_REF] Broggini | Data-driven wavefield focusing and imaging with multidimensional deconvolution: Numerical examples for reflection data with internal multiples[END_REF]. The MDD method is referred to as conventional Marchenko method.

More recently, redatuming by double focusing was developed and was shown to produce cleaner images than the MDD method when applied to 2D field data [START_REF] Staring | Source-receiver marchenko redatuming on field data using an adaptive double-focusing method[END_REF][START_REF] Wapenaar | Marchenko redatuming, imaging and multiple elimination, and their mutual relations[END_REF].

I.3.3.3 Phase distortions

In the previous sections, we showed that multiple reflections and multiple scattering were fundamental obstacles for deep imaging. However, even at shallower depths, seismic waves can undergo phase distortions due to long-scale fluctuations of the wave velocity.

Such phenomena induce aberrations that lead to loss of resolution and contrast in the beamformed image.

Actually, migration requires a correct wave velocity model. Figure I.9: Aberrations effects on focusing inside the medium. During focusing in emission and reception, waves travel through heterogeneities while propagating between the Earth's surface and the focusing plane, and thus suffer from aberrations. Consequently, the dimensions of the synthesized virtual sources and receivers are enlarged compared to their dimensions in absence of aberrations (diffraction limit).

is not true in reality. The fluctuations of seismic velocity in the Earth is challenging for seismic migration techniques and the modeling of wave propagation. If not taken into account during focusing operations, lateral velocity variations strongly distort the propagating wave-fronts resulting in a degradation of the image resolution and contrast. In other words, with an incorrect velocity model, waves will not focus correctly inside the medium and the coherent summation process is no longer optimal.

The distortions of the wave-fronts is commonly known as aberrations. Both incident and reflected waves undergo such phase distortions during their travel between the Earth's surface and the focusing plane. While, in absence of aberrations, the dimensions of the virtual source or receiver are only limited by diffraction, aberrations enlarge and distort their amplitude distribution (Fig. I.9). As said before, it results in a badly resolved and contrasted image, as we will show in this thesis.

A way to handle this problem in conventional migration techniques is to estimate migration velocity. Migration is usually combined with methods that estimate the accurate velocity of the subsurface, increasing the computational cost of the imaging process. Seismic tomography is the most popular method to build the velocity model. More recently, Full Waveform Inversion was shown to provide high-resolution velocity models of the subsurface [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF].

Several approaches have been developed in optics and acoustics to overcome aberrations effects and correct the focusing process. Astronomy was the first research field to deal with aberrations. While observing stars, the turbulence in the atmosphere affects the trajectory of light and results in the spreading of the stars images [START_REF] Tyson | Principles of adaptive optics[END_REF]. Aberrations correction relies first on the measurements of the wave-front distortions and then on the compensation of these distortions using Deformable Mirrors (DM) whose surfaces are deformed to be adapted with the measured wave-front [START_REF] Freeman | Deformable mirrors for all seasons and reasons[END_REF]. This is similar to apply a phase conjugate of the wave-front in the Fourier domain and constitutes the principle of adaptive optics. Similarly, in ultrasound imaging, time delays can be added to correct for low order aberrations undergone by the waves while propagating inside the medium, in emission and in reception [START_REF] Flax | Phase-aberration correction using signals from point reflectors and diffuse scatterers: Basic principles[END_REF]. ). One important property of an aberrating layer is its coherence length. This parameter is crucial since it indicates the area over which the aberration can be considered as spatially invariant. This is the so-called isoplanatic patch [START_REF] Loos | Turbulence of the upper atmosphere and isoplanatism[END_REF][START_REF] Roddier | Adaptive optics in astronomy[END_REF]. All the focusing points inside that region are subject to the same aberration at emission and at reception, resulting in strong correlations between the phase distortions undergone by each wave-front focused on or coming from this area. To better understand this phenomenon, we consider two focusing points corresponding to the same isoplanatic patch. The same phase law is applied for both focusing points (Fig.

I.10c).

Outside this specific region, the phase law does not compensate for aberrations anymore (Fig. I.10d). This phenomenon is known as isoplanatism [START_REF] Roddier | Adaptive optics in astronomy[END_REF] or memory effect [START_REF] Feng | Correlations and fluctuations of coherent wave transmission through disordered media[END_REF][START_REF] Freund | Memory effects in propagation of optical waves through disordered media[END_REF][START_REF] Liu | Physical picture of the optical memory effect[END_REF]. The size of isoplanatic volumes decreases with the heterogeneity of the medium and the penetration depth. The ratio between the imaged and isoplanatic volumes yields the number of aberration laws that should be determined in order to improve the focusing in every part of the image. This crucial aspect will be investigated in the next chapters.

We have seen that to compensate for low order aberrations in ultrasound imaging, additional time delays should be inserted in the incident and backscattered wave-field to focus inside the medium. In the next section, we will review the techniques adopted in ultrasound to overcome high-order aberrations. These techniques constitute the state of the art of a novel approach that has been recently developed by the team of Alexandre Aubry at Institute Langevin and applied to the fields of optics [START_REF] Badon | Approche matricielle de l'imagerie optique des milieux diffusants[END_REF][START_REF] Barolle | Approche matricielle de la tomographie à cohérence optique[END_REF], seismology [START_REF] Blondel | Matrix approach of seismic imaging[END_REF] and recently in acoustics [START_REF] Lambert | Matrix approach for ultrasound imaging and quantification[END_REF]. It was shown that the approach successfully improves the image quality in medical imaging and paves the way towards a reliable retrieval of the Earth's internal structure at an unprecedented resolution in seismology.

I.4 Reflection matrix approach of wave imaging

In the previous sections, we reviewed how a seismic wave-field recorded at the Earth's surface can be focused numerically inside the medium. Focusing is achieved by applying time delays or phase shifts in the frequency domain and allows to retrieve the reflectivity of the medium. In ultrasound imaging, focusing strategies through complex media have been developed for decades. These techniques are based on a matrix formalism of wave propagation and do not require prior knowledge on the medium. The time reversal experiment, that was introduced by Fink and his team [START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF][START_REF] Fink | Time-reversal waves and super resolution[END_REF], followed DORT (french acronym for "Decomposition of the Time Reversal Operator") method that was introduced by Prada and Fink [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF][START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF], constitute the foundation of matrix imaging in complex media. These works led to the development of a reflection matrix approach dedicated to selective focusing and target detection though aberrating and multiple scattering media.

In the following, we propose a brief review on time reversal and the decomposition of the time-reversal operator (DORT). Finally, we will discuss more recent developments of matrix imaging that is studied in the subsequent chapters to image the structure of two majors fault zones.

I.4.1 The time reversal concept

Focusing inside heterogeneous media is challenging. Inspired by the concept of phaseconjugate mirrors in optics (section I.3.3.3), time reversal concept was developed in acoustics to enable the focusing of waves through complex media. It relies on the timereversal invariance and reciprocity of the wave equation. In order to understand this prop- The emitted wave-field refocus back at the original source position. Figures are adjusted from [START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF].

erty, we consider the acoustic wave equation:

∆Φ - 1 c 2 ∂ 2 Φ ∂t 2 = 0 (I.24)
with Φ the wave-field and c the speed of acoustic waves. This equation contains a secondorder time derivative, consequently, if Φ(r, t) is a solution of Eq. I.24 then Φ(r, -t) is also a solution. This means that for every wave-field Φ(r, t) generated by a source at time t = 0 that can be subsequently reflected, scattered and refracted by the medium, a backward propagating wave-field Φ(r, -t) exists and retraces these paths to refocus at the source at t = 0 [START_REF] Fink | Self focusing in inhomogeneous media with time reversal acoustic mirrors[END_REF][START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF]. This is the principle of focusing waves inside an unknown medium without the need of prior knowledge on the propagating medium.

This concept was illustrated with the time-reversal cavity (TRC). The TRC experiment, illustrated in Figures are adjusted from [START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF].

then it is time-reversed and back-propagated inside the medium. The time-reversed wavefield compensates for the distortions induced by the heterogeneous medium and focuses back at the original source position.

Practically, the TRC is difficult to implement since it is not always possible to surround the medium along a closed surface. TRC is replaced by limited angular aperture time-reversal mirror (TRM) that consists of a linear array of transducers located at one side of the medium of interest [START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF][START_REF] Fink | Time-reversal mirrors[END_REF]. Similarly to the TRC, the focusing on a target through an inhomogeneous medium is done through three steps: (i)a wave-front is transmitted into the medium, (ii) the target acts as a secondary source and generates a backscattered wave-field that is recorded by the array of transducers, and (iii) the transducers time-reverse the recorded wave-field and send it back into the medium. The generated wave-field propagates through the heterogeneities of the medium and focuses on the target. In a medium with several targets, the process can be iterated to converge toward the most reflective target since the wave-field focusing on the weak targets vanishes [START_REF] Prada | The iterative time reversal mirror: A solution to selffocusing in the pulse echo mode[END_REF][START_REF] Prada | Theory of iterative time reversal acoustic mirrors[END_REF][START_REF] Prada | The iterative time reversal process: Analysis of the convergence[END_REF].

The time reversal principle can be understood as a phase conjugation process. Consider Φ(r, ω) the temporal Fourier transform of Φ(r, t), then the temporal Fourier transform Φ(r, -t) is Φ * (r, ω). Consequently, emitting a time-reversed wave-field is equivalent to apply the phase conjugate of the incident wave-field that compensates the aberrations and allows to focus inside the medium. For more details on the time reversal concept and experiments, the reader is referred to the review papers of [START_REF] Fink | Time-reversed acoustics[END_REF] and [START_REF] Fink | Acoustic time-reversal mirrors[END_REF].

The time-reversal focusing constitutes a very robust tool to focus waves in media with high order multiple scattering [START_REF] Derode | Robust acoustic time reversal with high-order multiple scattering[END_REF]. In fact, in multiple scattering regime, the medium acts as a lens and the angular aperture of the array is widened due to scattering, resulting in a better resolution than the theoretical diffraction-limited resolution that is in the order of λz/D, with D the experimental array aperture. These observations open the way for numerous applications using the TRM concept in many fields, such as medical imaging [START_REF] Ikeda | An image reconstruction algorithm using phase conjugation for diffractionlimited imaging in an inhomogeneous medium[END_REF][START_REF] Dorme | Focusing in transmit-receive mode through inhomogeneous media: The time reversal matched filter approach[END_REF], non-destructive testing [START_REF] Chakroun | Time reversal processing in ultrasonic nondestructive testing[END_REF] and ocean acoustics [START_REF] Hodgkiss | A long-range and variable focus phase-conjugation experiment in shallow water[END_REF][START_REF] Kim | Robust time reversal focusing in the ocean[END_REF].

I.4.2 Decomposition of the time reversal operator (DORT)

We have seen that the TRM is adapted to focus on a single target, or the most reflective one, inside a heterogeneous medium. In order to focus on several reflectors in a multitarget medium, the DORT method has been introduced [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF][START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF]. It consists in a matrix approach of the iterative time-reversal process. This method is based on the diagonalization of the time-reversal operator. To understand the concept, we will describe the time-reversal process under a matrix formalism. As already describe in section I.3.2, the reflection matrix contains the response of the medium between a set of sources and receivers. We define R(ω) the temporal Fourier transform of the response matrix. The system is considered linear and time-invariant. The reflected wave-field Φout (ω) is linked to the incident wave-field Φin (ω) through the following relation:

Φ(0) out = R × Φ(0) in (I.25)
In the frequency domain, the time-reversal operation consists in considering the phase conjugate of the output signal Φ * (0) out as the new input wave-field. The resulting wave-field at iteration 1 is:

Φ(1) out = R × Φ * (0) out = [R × R * ] × Φ * (0) in (I.26)
The operator [R × R * ] corresponds to the time-reversal operator. Due to source-receiver reciprocity, the reflection matrix is symmetric, hence R * = R † , where the symbol † stands for transpose conjugate. The operator [R × R † ] is Hermitian and can be diagonalized. DORT method consists then in performing an eignevalue decomposition of RR † , or equivalently a singular value decomposition of R. The singular value decomposition of RR † reads:

RR † = UΣU † (I.27)
where Σ is a diagonal matrix containing the real positive eigenvalues σ i in a decreasing [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF]. U is a unitary matrix whose columns correspond to the eigenvectors of RR † . In the single scattering regime and for point-like target, it has been shown that the number of significant eigenvalues is equal to the number of targets inside the medium. Moreover, there is a one-to-one association between each target and each eigenstate of RR † . The phase conjugate of the each eigenvector should be applied to focus on the corresponding target. The ranking of the eigenvalues scales as the reflectivity of the targets: the first eigenvalue (i.e the largest) corresponds to the most reflective target. [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF].

order σ 1 > σ 2 > • • • > σ N (
elevation [m] elevation [m] x [m] x [m] y [ m ] y [ m ]
The efficiency of DORT in separating several scatterers was demonstrated for homogeneous and heterogeneous media [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF][START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF]. It constitutes a powerful tool for target detection in imaging. It was applied for non-destructive testing to detect flaws in solid [START_REF] Prada | Flaw detection in solid with the dort method[END_REF][START_REF] Prada | Detection and imaging in complex media with the dort method. Imaging of complex media with acoustic and seismic waves[END_REF], in medical imaging and acoustics [START_REF] Robert | Time reversal operator decomposition with focused transmission and robustness to speckle noise: Application to microcalcification detection[END_REF][START_REF] Cochard | Ultrasonic focusing through the ribs using the dort method[END_REF], in underwater acoustics [START_REF] Mordant | Highly resolved detection and selective focusing in a waveguide using the dort method[END_REF][START_REF] Lingevitch | Time reversed reverberation focusing in a waveguide[END_REF], in optics [START_REF] Popoff | Exploiting the time-reversal operator for adaptive optics, selective focusing, and scattering pattern analysis[END_REF] and in electromagnetism [START_REF] Micolau | Dort method as applied to ultrawideband signals for detection of buried objects[END_REF][START_REF] Antoine | Far field modeling of electromagnetic time reversal and application to selective focusing on small scatterers[END_REF].

DORT was also applied in acoustics for target detection in multiple scattering regime combined with a single-scattering filter that extracts the single scattered waves when multiple scattering dominates (Aubry and Derode, 2009a,b). This approach was also extended to optics [START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF], and recently, it has been applied to seismology to image through multiple scattering in the case of the Erebus volcano in Antarctica [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]. The internal structure of the volcano was retrieved using ambient noise body waves (Fig.

I.13).

Those matrix methods are powerful to detect and image targets in a strong multiple scattering regime. Nevertheless the image still suffers from aberrations and is far to be diffraction limited. During the last few years, a lot of efforts have been made to develop post-processing methods for aberration correction under a matrix framework. This approach will be referred to as matrix imaging in the following.

I.4.3 Matrix imaging

After the pioneering works of Fink and colleagues on the iterative time-reversal process [START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF][START_REF] Fink | Time-reversal mirrors[END_REF] and on the decomposition of the time-reversal operator [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF][START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF]Aubry and Derode, 2009a,b), the matrix approach has been extended to imaging purposes in heterogeneous media since the last few years.

The principle of this method is to gather information on the medium from the response matrix recorded between a set of emitters and receivers. It has been shown that this technique is a powerful tool to improve the focusing in all kind of scattering regimes. It allows to overcome the problem of distributed aberrations that pollute imaging process in many fields of wave physics [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2020b).

As described in the previous section, focusing operations can be applied to the response matrix of the medium to obtain a focused reflection matrix that contains the response between virtual sources and receivers inside the medium [START_REF] Aubry | Ultrasonic imaging of highly scattering media from local measurements of the diffusion constant: Separation of coherent and incoherent intensities[END_REF][START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF][START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF]Lambert et al., 2020aLambert et al., , 2021a)). The calculated matrix contains all the information on the medium. A set of matrix operations is thus applied to characterize the propagation medium and retrieve reflectivity images with a close-to-ideal resolution at any pixel.

To do so, a novel operator, the so-called, distortion matrix D has been introduced.

It is inspired by the FDORT method developed by Robert and Fink in ultrasound imaging [START_REF] Robert | Time reversal operator decomposition with focused transmission and robustness to speckle noise: Application to microcalcification detection[END_REF]. FDORT is a generalization of DORT method, where the singular value decomposition is applied to the response matrix between virtual sources inside the medium and sensors at the surface. This method is particularly suitable for aberration correction in medical imaging where the echographic images are distorted due to variations of the speed-of-sound [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF].

D-matrix essentially connects any focal point inside the medium with the distortion that a wave front, emitted from that point, experiences due to heterogeneities. The phase laws that should be applied to compensate aberrations and optimize the focusing process are estimated from a singular value decomposition of the D-matrix. This novel concept has been tested on well-controlled systems in ultrasound (Lambert et al., 2020b) and I.5. CONCLUSION optics [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF], for a random distribution of unresolved scatterers (speckle) and specular reflectors, respectively. It was also applied to in vivo imaging of human soft tissues (Lambert et al., 2021b). Compared to conventional adaptive focusing techniques, its strength lies in its ability to extract an adapted focusing for each isoplanatic patch in the field-of-view. The transmission matrix between each sensor used to probe the medium and any voxel inside it is actually retrieved. This matrix constitutes the Holy Grail for imaging as it enables a close-to-ideal resolution and optimized contrast for any pixel of the image.

In this thesis, I will show, in the next chapters, how the distortion matrix concept can be a powerful tool in passive seismology by exploiting body waves reflected by the heterogeneities of the Earth's crust. Despite a rough wave velocity model, matrix imaging will succeed in revealing the structure of fault zones by revealing both main specular reflectors and distribution of more localized inhomogeneities in fault zones, at different scales and over different frequency bandwidths.

I.5 Conclusion

Since faults are the sources of destructive earthquakes, earthquake seismology has know an exponential growth, especially with the emerging of passive seismology. This first chapter was dedicated to highlight the complexity of fault zones and to review a number of notions that are crucial to image and understand the structure of faults.

Numerous conventional techniques are used to image the Earth's crust. But these techniques needs information and accurate velocity model on the medium, which limits their applications. Adaptive focusing strategies and the time-reversal concept constitute a way to overcome these challenges and pave the way to develop new imaging techniques for geophysical applications. Another promising approach is the extraction of the Green's function from the cross-correlation of passive seismic measurements.

The aforementioned notions, combined together, establish a new passive seismic imaging technique in reflection that can be applied to any heterogeneous geological environment without prior knowledge on the medium's properties. The main objective of this thesis is to exploit the recent developments made in ultrasound imaging and optical microscopy in seismology. The principle of matrix imaging will be detailed in the following chapters.

Two applications are discussed aiming to image the structure of two major fault zones: This chapter is devoted to the application of the reflection matrix approach at a small scale, for the case of San Jacinto Fault. The two sections present two published papers:

San
the first paper introduces the imaging technique and the second paper discussed the interpretation of the structural images around the fault. The subsurface here exhibits a sparse distribution of scattering.

II.1 A Distortion Matrix Framework for High-Resolution

Passive Seismic 3D Imaging: Application to the San Jacinto Fault Zone, California R. Touma at al.]Rita Touma

II.1.1 Abstract

Reflection seismic imaging usually suffers from a loss of resolution and contrast because of the fluctuations of the wave velocities in the Earth's crust. In the literature, phase distortion issues are generally circumvented by means of a background wave velocity model.

However, it requires a prior tomography of the wave velocity distribution in the medium, which is often not possible, especially in depth. In this paper, a matrix approach of seismic imaging is developed to retrieve a three-dimensional image of the subsoil, despite a rough knowledge of the background wave velocity. To do so, passive noise cross-correlations between geophones of a seismic array are investigated under a matrix formalism. They form a reflection matrix that contains all the information available on the medium. A set of matrix operations can then be applied in order to extract the relevant information as a function of the problem considered. On the one hand, the background seismic wave velocity can be estimated and its fluctuations quantified by projecting the reflection matrix in a focused basis. It consists in investigating the response between virtual sources and detectors synthesized at any point in the medium. The minimization of their cross-talk can then be used as a guide star for approaching the actual wave velocity distribution. On the other hand, the detrimental effect of wave velocity fluctuations on imaging is overcome by introducing a novel mathematical object: The distortion matrix. This operator essentially connects any virtual source inside the medium with the distortion that a wave-front, emitted from that point, experiences due to heterogeneities. A time reversal analysis of the distortion matrix enables the estimation of the transmission matrix that links each real geophone at the surface and each virtual geophone in depth. Phase distortions can then be compensated for any point of the underground. Applied to passive seismic data recorded along the Clark branch of the San Jacinto fault zone, the present method is shown to provide an image of the fault until a depth of 4 km over the frequency range 10-20 Hz with an horizontal resolution of 80 m. Strikingly, this resolution is almost one eighth below the diffraction limit imposed by the geophone array aperture. The heterogeneities of the subsoil play the role of a scattering lens and of a transverse wave guide which increase drastically the array aperture. The contrast is also optimized since most of the incoherent noise is eliminated by the iterative time reversal process. Beyond the specific case of the San Jacinto Fault Zone, the reported approach can be applied to any scales and areas for which a reflection matrix is available at a spatial sampling satisfying the Nyquist criterion.

II.1.2 Introduction

Waves constitute a powerful means to non destructively probe an unknown medium. Indeed, wave propagation is fully determined by the wave equation and boundary conditions. As stated by diffraction theory in the acoustic approximation, knowing the incident wave-field and the internal properties in terms of density ρ and celerity c theoretically allows to compute the wave field everywhere and at any time inside the medium. This is the so-called "forward" or "modeling" problem. Conversely, when the medium is unknown, the lack of celerity and density knowledge makes it impossible to compute the spatiotemporal evolution of the wave field. Yet this evolution can be known, at least at the boundary, through experimental measurement of the wave-field scattered by the medium.

The inverse problem then consists in deducing the medium internal properties from the recording of the wave field at its surface. A first way to do so is to assume a velocity and density background, solve the forward problem to compute the time-dependent signal that would be backscattered if the background model was true, and iteratively update this model to minimize the difference with the actual recordings. Another way is to directly back-propagate the scattered echoes to reflectors inside the medium. This also amounts to updating a background model since a reflector is nothing else than a variation in acoustic impedance ρc. In both strategies, referred to as "inversion" and "migration", respectively, a celerity macro model is required and the purpose is to compute variations from this model under the assumption that they are small (Born approximation). If they are not, the reflected wave-field may be subject to aberrations and multiple scattering that the macro model fails at modeling. These issues lead to distorted images, lack of resolution and nonphysical features, which are very detrimental to the imaging process.

In seismic exploration, these issues are important because in most cases the celerity is non constant in space and its distribution is unknown. Most geological settings actually consist of several layers of rocks and sediments with distinct mechanical properties as well as location-dependent thickness, faulting and strata organization. These may be difficult to estimate without previous geologic expertise of the subsurface, especially in areas with high lateral mechanical stress that bend and break the layers and make them superimpose.

When trying to retrieve details at a diffraction-limited resolution, the previous knowledge required to build reliable images may be already fairly demanding. On the one hand, the inversion problem cannot be solved if the initial celerity model is too far from the reality because the regularization procedure can end up stuck in a local minimum. On the other hand, migration techniques would lead to loss of resolution due to phase distortions and a blurred image due to multiple scattering. That being said, the question that naturally arises and which the present work aims at addressing is: how to retrieve an accurate image when little to no previous knowledge on the spatial variations of the wave speed is available?

To cope with this issue, our strategy is to develop a matrix approach of seismic wave imaging. In a linear scattering medium, a reflection matrix relates the input and output on a single side of the medium. It contains all the relevant information on the medium as it fully describes wave propagation inside the scattering medium. In the last decade, the advent of multi-element arrays with controllable emitters and receivers has opened up a route towards the ability of measuring a reflection matrix in the case where the input and outputs points are located on the same side of the scattering medium. In particular, the reflection matrix has been shown to be of great interest for detection and imaging purposes in scattering media, whether it be in acoustics [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF]Aubry and Derode, 2009a;[START_REF] Shahjahan | A random matrix approach to detect defects in a strongly scattering polycrystal: How the memory effect can help overcome multiple scattering[END_REF] or optics [START_REF] Kang | Imaging deep within a scattering medium using collective accumulation of single-scattered waves[END_REF][START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF]. The reflection matrix contains the set of inter-element impulse responses recorded between each array element. It has already been shown to be a powerful tool for focusing in multi-target media [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF][START_REF] Popoff | Exploiting the time-reversal operator for adaptive optics, selective focusing, and scattering pattern analysis[END_REF], as well as for separating single and multiple scattering (Aubry and Derode, 2009b;[START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF] in strongly scattering media. These matrix methods have been successfully applied to geophysics in the extremely challenging case of the Erebus volcano in Antarctica [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF].

Seismic data redatuming leads to the synthesis of a focused reflection matrix R containing the impulse responses between a set of virtual geophones mapping the underground to be imaged. This matrix is of particular interest for imaging: its confocal component i.e., diagonal (or close to diagonal) elements result from single scattering, whereas multiple scattering is responsible for a spreading outside the diagonal. By applying an adaptive confocal filter and iterative time reversal to the redatumed data, most of the multiple scattering background is removed, thereby revealing main internal structures of the volcano.

Yet, the resulting "confocal" image still suffers from the phase distortions induced by the long-scale fluctuations of the seismic velocity. In the present paper, we show that the off-diagonal coefficients of R can be taken advantage of instead of being tossed away.

Phase distortions are overcome, which improves the confocal image of the subsoil with a diffraction-limited transverse resolution.

To that aim, we will rely on the distortion matrix concept that has been recently introduced by two seminal works in ultrasound imaging (Lambert et al., 2020b) and in optical microscopy [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]. Inspired by the pioneering work of [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF], the distortion matrix D is defined between a set of incident plane waves [START_REF] Montaldo | Coherent plane wave compounding for very high frame rate ultrasonography and transient elastography[END_REF] and the set of virtual geophones inside the medium [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF]. It contains the deviations from an ideal reflected wave-front which would be obtained in the absence of inhomogeneities. As shown by recent studies [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2020b), a time reversal analysis of the D-matrix allows to synthesize virtual reflectors in depth. This process can then be leveraged for unscrambling the phase distortions undergone by the incident and reflected wave-fronts.This matrix imaging approach has been shown to be particularly robust since it applies to all kind of scattering regimes: point-like targets (Lambert et al., 2020b), specular reflectors [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF] or randomly distributed scatterers (Lambert et al., 2020b) etc. The range of wave velocity distributions on which this matrix method can be applied is vast. It goes from the simple case of multi-layered media (Lambert et al., 2020b) to strongly heterogeneous media (Lambert et al., 2020b;[START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF] displaying both lateral and axial variations of the wave velocity. These proof-of-concept studies made on both synthetic samples, exvivo or in-vivo tissues have thus shown the power of the matrix imaging approach to solve very different imaging problems whatever the nature of waves and length scales. It only requires a spatial sampling of the recorded wave-field that meets the Nyquist criterion.

The present paper aims at demonstrating the relevance of this matrix imaging approach for geophysical imaging.

Indeed, overcoming phase distortions induced by wave velocity variations would be especially valuable for geophysical applications given the stratified structures of the environments of interest. Migration techniques in Fourier domain have actually been very popular for imaging in layered media [START_REF] Stolt | Migration by Fourier transform[END_REF][START_REF] Gazdag | Wave equation migration with the phase-shift method[END_REF], however they only hold for 1D celerity models with no lateral variations. Subsequent works have focused on adapting these techniques to take into account increasing lateral velocity variations, at the cost of more numerical and computational complexity [START_REF] Gazdag | Migration of seismic data by phase shift plus interpolation[END_REF][START_REF] Claerbout | Imaging of of the Earth's Interior[END_REF][START_REF] Biondi | 3D seismic imaging[END_REF]. Contrary to these well-established methods, the matrix approach only requires an estimation on the velocity distribution inside the medium, while being fairly light on the computational aspect. The present paper aims at studying the relevance of the matrix approach for geophysical imaging.

Coherent sources (vibrating trucks, explosives, etc.) can be used in shallow subsurface (<1 km) imaging. Incoherent signals (seismic noise) can also be taken advantage of for imaging purposes. It was shown, twenty years ago, how a coherent information can be extracted from this incoherent seismic noise. Under appropriate wave-field condi-tions, the cross-correlation of seismic noise recorded by two stations was actually shown to yield the impulse response between them [START_REF] Weaver | Ultrasonics without a source: Thermal fluctuation correlations at mhz frequencies[END_REF][START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF]Derode et al., 2003b;[START_REF] Wapenaar | Retrieving the elastodynamic green's function of an arbitrary inhomogeneous medium by cross correlation[END_REF][START_REF] Snieder | Extracting the green's function from the correlation of coda waves: A derivation based on stationary phase[END_REF][START_REF] Larose | Correlation of random wavefields: An interdisciplinary review[END_REF], providing new opportunities to develop imaging techniques without using active sources.

As surface waves dominate ambient noise, most papers on the topic aimed at extracting surface wave properties from ambient noise correlations [START_REF] Shapiro | High-resolution surfacewave tomography from ambient seismic noise[END_REF]Sabra et al., 2005a;[START_REF] Yang | Ambient noise rayleigh wave tomography across europe[END_REF]. However, a few studies also reported the retrieval of body wave reflection from noise correlations (Roux et al., 2005a;[START_REF] Draganov | Retrieval of reflections from seismic background-noise measurements[END_REF][START_REF] Draganov | Reflection images from ambient seismic noise[END_REF]Poli et al., 2012a). Reflected body waves contain information about the subsurface and allow the imaging of deep structures with an improved resolution [START_REF] Ruigrok | High-resolution lithospheric imaging with seismic interferometry[END_REF].

Strikingly, Poli et al. (2012b) showed the possibility of mapping the upper mantle discontinuities (at 410 and 660 km of depth) by extracting body waves reflection from ambient noise, while Retailleau et al. ( 2020) mapped a region of the core-mantle boundary at about 2900 km depth.

In this paper, inspired by the reflection matrix approach developed by Blondel et al.

(2018) and based on noise cross-correlations, the distortion matrix approach is extended to satisfy seismic imaging purposes. Compared to previous works [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2020a), the matrix imaging method is here refined to take the best advantage of the nature of scatterers in geophysics (sparse scattering). The method thus developed is applied to San Jacinto fault zone (SJFZ) site. Fault zones are indeed among the most challenging media for seismic imaging given their highly localized and abrupt variations of mechanical properties, extensive fractures and damage zones. In that respect, the SJFZ is the most seismically active fault zone in Southern California [START_REF] Hauksson | Waveform relocated earthquake catalog for southern california (1981 to june 2011)[END_REF]. It accounts for a large portion of the plate motion in the region [START_REF] Johnson | Present-day crustal deformation in southern California[END_REF][START_REF] Lindsey | Geodetic slip rates in the southern San Andreas Fault system: Effects of elastic heterogeneity and fault geometry[END_REF].

A highly complex fault-zone structure with prominent lateral and vertical heterogeneities at various scales have already been highlighted in previous studies [START_REF] Allam | Seismic velocity structures in the southern California plateboundary environment from double-difference tomography[END_REF][START_REF] Zigone | Seismic tomography of the southern california plate boundary region from noise-based rayleigh and love waves[END_REF][START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF]. In particular, maps of the P and S wave velocities, V P and V S , have been inverted from earthquake arrival times for a depth range of 2-20 km [START_REF] Allam | Seismic velocity structures in the southern California plateboundary environment from double-difference tomography[END_REF][START_REF] Allam | Seismic velocity structure in the hot springs and trifurcation areas of the San Jacinto fault zone, California, from doubledifference tomography[END_REF]. Surface wave tomo-graphic images built from noise correlations revealed the velocity structure in the top 7 km of the complex plate boundary region at a resolution of about ten kilometers [START_REF] Zigone | Seismic tomography of the southern california plate boundary region from noise-based rayleigh and love waves[END_REF]. To complement these regional studies and provide structural features in the first few kilometers with an improved resolution, ambient noise at higher frequency up to 10 Hz was analysed from data recorded by a dense rectangular array deployed around the Clark branch of the SJFZ [START_REF] Ben-Zion | Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone[END_REF][START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF][START_REF] Mordret | Shallow three-dimensional structure of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic array[END_REF].

In Imaging deeper the fault area at such resolution is challenging because of the damage and the complex distribution of small-scale heterogeneities. Yet, a much larger penetration depth can be expected by taking advantage of the reflected bulk waves. To do so, the matrix approach of seismic imaging is particularly useful since it only requires a rough idea of the mean wave velocity. Besides, it shall provide a three-dimensional image of the subsoil acoustic impedance instead of just the wave velocity. To implement this matrix approach, we take advantage of a spatially dense array of geophones deployed over the damage zone of SJFZ [START_REF] Ben-Zion | Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone[END_REF]. Noise cross-correlations are used to retrieve the impulse responses between the geophones [START_REF] Ben-Zion | Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone[END_REF][START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF]. The associated passive reflection matrix is then investigated to image the first few kilometers of the crust by virtue of body waves emerging from noise correlations. As a whole, the process we present in this paper can be analyzed as a combination of six building blocks:

• (B1) A Fourier transform of the recorded signals yields a set of response matrices K(f ) associated with the dense array of geophones.

• (B2) Based on a rough estimate of velocity c 0 , a double focusing operation is performed both at emission and reception by means of simple matrix operations. A set of focused reflection matrices R(f, z) are obtained at any arbitrary depth z below the surface.

• (B3) A coherent sum of these matrices over the frequency bandwidth yields a broadband reflection matrix R(z) at any depth z.

• (B4) By projecting the input or output entries of this matrix in the far-field, the distorted component D of the reflected wave-field can be extracted.

• (B5) A virtual iterative time reversal process is applied to the matrix D to extract the phase distortions undergone by the incident or reflected wave-fields during their travel from the Earth surface to the focal plane.

• (B6) The whole process converge towards the focusing laws that shall be applied at input or output of the reflection matrix in order to compensate for aberrations.

As a result of these six steps, an in-depth confocal image of the SJFZ is built. While Each row along the x-direction is composed of ∼55 sensors with a pitch of 10 m, and the nominal separation between the rows in the y-direction is 30 m. Seismic ambient noise was recorded over more than one month, in May-June 2014. (b) Adaptive focusing at emission and reception on two points r in and r out of the focal plane (z = c 0 t/2) yields the impulse response between virtual geophones placed at these two points. The same operation is repeated for any couple of points in the focal plane and yields the focused reflection matrix R.

Given the high density of the network, neighbouring geophones belong to the same coherence area of seismic noise. The characteristic dimension of this area is indeed of λ/2 ∼ 50 m which is larger than the inter-station distance δu. This is responsible for a strong auto-correlation signal around t = 0 for geophones located in the same coherence area. This peak is proportional to the seismic noise power and does not account for the impulse response between neighbour geophones. To prevent this artifact from spoiling the subsequent analysis, a prior filter has been applied to the data in order to reduce the weight of the corresponding impulse responses k ij (t) whose associated geophones i and j are contained in the same coherence area (see Supplementary Section S1).

The impulse responses exhibit several direct arrivals that have already been investigated by [START_REF] Ben-Zion | Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone[END_REF] and [START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF]. Ballistic waves, likely direct inter-station S-wave and P-wave, arrive before the Rayleigh wave at apparent velocities larger than 1000 m/s. [START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF] used iterative double beamforming to map the phase and group velocities of Rayleigh waves across the fault in the 1-5 Hz frequency bandwidth. Subsequently, [START_REF] Mordret | Shallow three-dimensional structure of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic array[END_REF] inverted these dispersion curves to build a 3-D shear wave velocity model around the Clark fault down to 500 m depth. Assuming the V p /V s ratio to be a linear function of depth, the following averaged value was found for the P-wave velocity over the top 800 m: V p ∼1500 m/s. More recently, the P-wave velocity distribution in the 100 m-thick shallow layer has also been inverted using travel time data associated with active shots [START_REF] Share | Characterizing the uppermost 100 m structure of the San Jacinto fault zone southeast of Anza, California, through joint analysis of geological, topographic, seismic and resistivity data[END_REF]. Low-velocity structures were detected, associated with a shallow sedimentary basin [START_REF] Hillers | Focal spot imaging based on zero lag cross-correlation amplitude fields: Application to dense array data at the San Jacinto fault zone[END_REF][START_REF] Mordret | Shallow three-dimensional structure of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic array[END_REF][START_REF] Share | Characterizing the uppermost 100 m structure of the San Jacinto fault zone southeast of Anza, California, through joint analysis of geological, topographic, seismic and resistivity data[END_REF] and a fault zone trapping structure [START_REF] Ben-Zion | Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone[END_REF][START_REF] Qin | Internal structure of the San Jacinto fault zone in the trifurcation area southeast of Anza, California, from data of dense seismic arrays[END_REF].

To the best of our knowledge, an accurate model of V p in the SJFZ region is not available beyond this shallow layer. As a consequence, in the present study, we will use an approximated homogeneous P-wave velocity model of c 0 = 1500 m/s. This choice will be validated and discussed a posteriori by a minimization of the aberration effects in the 3D image (see Supplementary Fig. S3). We are not interested in the ballistic component of the wave-field but rather in its scattered contribution due to reflections by the in-depth structure along the fault. The beamformed echoes used in our matrix imaging process are mainly associated with P-waves since only the vertical component of the impulse responses between geophones is considered in this study. Unlike our previous study on the Erebus volcano [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF], the scattered wave-field consists of a single scattering contribution which is a priori largely predominant compared to the multiple scattering background. This will be confirmed a posteriori by the reflection matrix features (see Sec. II. 1.3.2). Singly-scattered echoes can then be taken advantage of to build a 3D image of the subsoil reflectivity. This local information can be retrieved from K(t) by applying appropriate time delays to perform focusing in post-processing, both in emission and reception. While focusing in emission consists in applying proper time delays in the recorded seismic data so that they constructively interfere at an arbitrary position at depth, focusing in reception consists in applying proper time delays in the recorded seismic data so that the information coming from an arbitrary position at depth constructively interfere. Based on the Kirchoff-Helmholtz integral, such a focusing operation is standard in exploration seismology and referred to as redatuming [START_REF] Berkhout | Imaging of Acoustic Energy by Wave Field Extrapolation[END_REF][START_REF] Berryhill | Wave equation datuming before stack[END_REF][START_REF] Berkhout | A unified approach to acoustical reflection imaging. II: The inverse problem[END_REF]. However, in the present case, the strongly heterogeneous distribution of the seismic wave velocities induces strong phase distortions that degrade this imaging process. A prior quantification and correction of these phase distortions is thus required to reach a diffraction-limited lateral resolution and an optimized contrast for the image. As we will see, a matrix formalism is a well-matched tool to locally capture such information.

II.1.3.2 Focused reflection matrix

The reflection matrix can be defined in general as an ensemble of responses, each response linking one vector to another vector. The type of vector coordinates will be referred to as bases. They can be spatial coordinates (hence the vector refers to an actual point within or at the surface of the medium, see Fig. II.1b) or wave vector coordinates. Various bases are involved in this work :(i) the recording basis (u), whose elements are the positions of the geophones; (ii) the focused basis (r) which corresponds to the positions of virtual geophones at which focusing at emission or reception is intended; and (iii) the Fourier basis (k). Because of linearity and time-invariance, seismic data can be projected from the recording basis to the focused basis by a simple matrix product. In the frequency domain, simple matrix products allow seismic data to be easily projected from the recording basis to the focused basis where local information on the medium properties can be extracted [START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]Lambert et al., 2020a).

Consequently, we first apply a temporal Fourier transform to the response matrix to obtain a set of monochromatic matrices K(f ). To project K(f ) into the focused basis, we then define a free-space Green's matrix, G 0 (f ), which describes the propagation of waves between the geophones and focused basis. Its elements correspond to the causal 3D Green's functions which connect the geophone's transverse position u to any focal point defined by its transverse position r and depth z in a supposed homogeneous medium:

G 0 (r, u, z, f ) = e -j2πf √ r-u 2 +z 2 /c 0 4π r -u 2 + z 2 (II.1)
K(f ) can now be projected both in emission and reception to the focused basis via the following matrix product at each depth z [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]Lambert et al., 2020a):

R(z, f ) = G * 0 (z, f ) × K(f ) × G † 0 (z, f ) , (II.2)
where the symbols * , † and × stand for phase conjugate, transpose conjugate and matrix product, respectively. Equation (II.2) simulates focused beamforming in post-processing in both emission and reception. The choice of causal Green's function in the propagation matrix G 0 (Eq. II.1) implies that beamformed echoes are associated with downgoing waves at the input and up-going waves at the output (see Fig. It yields the broadband focused reflection matrix

R(z) = ˆf+ f - df R(z, f ), (II.3)
where f ± = f 0 ± ∆f /2 and f 0 = 15 Hz is the central frequency. Each element of R(z)

contains the complex amplitude of the wave that would be detected by a virtual detector located at r out = (x out , y out ) just after a virtual emitter at r in = (x in , y in ) emits a brief pulse of length δt = ∆f -1 at the central frequency f 0 . Importantly, the broadband focused reflection matrix synthesizes the responses between virtual geophones which have a greatly reduced axial dimension δz = cδt compared to their stretching δz 0 = 2λ/ sin 2 θ in the monochromatic regime (Born and Wolf, 2003). θ = arctan(D/2z) is the maximum angle under which the geophone array is seen from the common mid-point and D ∼ 700 m, the characteristic size of the geophone array. As a consequence, considering a broadband reflection matrix R(z) will significantly improve the vertical resolution of the subsequent analysis. For the sake of a lighter notation, we will omit, in the following, the dependence in z but keep in mind that the focused reflection matrix differs at each depth. [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]. This indicates that single scattering dominates at this depth: The beam is focused, scattered just once, and focused in reception. On the contrary, a broadening of the back-scattered energy outside the diagonal would mean that the beam undergoes aberration and/or multiple scattering.

In fact, the diagonal elements of R (r in = r out ) correspond to what would be obtained from confocal imaging: transmit and receive focusing are simultaneously performed on each point in the medium. A confocal image can thus be obtained from the diagonal elements of R, computed at each depth:

I (r, z) ≡ |R (r, r, z)| 2 .
(II.4) Fortunately, the matrix R contains much more information than a single confocal image. In particular, focusing quality can be assessed by means of the off-diagonal elements of R. To understand why, R can be expressed theoretically as follows (Lambert et al., 2020a(Lambert et al., , 2021a)):

0 0 0 3 -0 0 0 2 -0 0 0 1 - Depth[m] 0 -1700 x [m]
R = H out × Γ × H in , (II.5)
where the symbol stands for transpose. The matrix Γ describes the scattering process inside the medium. In the single scattering regime, Γ(z) is diagonal and its coefficients map the local reflectivity γ(r) of the subsoil.

H out = [H out (r, r out )] and H in = [H in (r, r in )]
are the output and input focusing matrices, respectively. Their columns correspond to the transmit or receive point spread functions (PSFs), i.e. the spatial amplitude distribution of the focal spots around the focusing point r in or r out . For spatially-invariant aberration,

we have

H out/in (r, r out/in ) = H out (r -r out/in ). (II.6)
In that case, the previous equation can then be rewritten in terms of matrix coefficients as follows:

R(r out , r in ) = ˆdrH out (r -r out )γ(r)H in (r -r in ). (II.7)
This last equation confirms that the diagonal coefficients of R, i.e. an horizontal slice of the confocal image, result from a convolution between the medium reflectivity γ and the product of the input and output PSF,

H out × H in .
Interestingly, the off-diagonal terms in the reflection matrix can be exploited to estimate the imaging PSF, and thereby assess the quality of focusing. To that aim, the relevant observable is the intensity distribution along each antidiagonal of R,

I(r m , ∆r) = R(r m -∆r, r m + ∆r) = ˆdr H out (r -∆r)γ(r + r m )H in (r + ∆r). (II.8)
All couple of points on a given antidiagonal have the same midpoint r m = (r out + r in )/2

, but different spacings ∆r = (r out -r in )/2. Whatever the nature of the scattering medium, the common midpoint intensity profile is a direct indicator of the local PSF.

However, its theoretical expression differs slightly depending on the characteristic length scale l γ of the reflectivity γ(r) at the ballistic depth and the typical width

δ (0)
in/out of the PSFs (Lambert et al., 2020a). In layered media (l γ >> δ (0) in/out ), the common-midpoint amplitude is directly proportional to the convolution between the coherent output and input PSFs, [H out ⊗ H in ] (2∆r) (the symbol ⊗ stands for convolution). In the speckle regime

( l γ << δ (0)
in/out ), the common midpoint intensity I(r m , ∆r) is directly proportional to the convolution between the incoherent output and input PSFs, bert et al., 2021a). In the present case, the subsoil of SJFZ can be assumed as a sparse scattering medium. It means that only a few bright and coherent reflectors emerge at each depth. This hypothesis will be verified a posteriori with the three-dimensional image we will obtain. For an isolated scatterer, the common mid-point intensity at its position scales as the product between the two PSFs, H out (r m -∆r) × H in (r m -∆r). Therefore, the energy spreading in the vicinity of each scatterer position shall enable one to probe the spatial extension of the PSF. As the scatterer position is a priori unknown, the imaging PSF will be, in practice, probed by considering the antidiagonal whose common mid-point exhibits the maximum confocal signal. to the gap between the velocity model and the actual seismic wave velocity distribution in SJFZ. These side lobes are analogous to cycle skipping effects that occur in fullwaveform inversion [START_REF] Yao | Tackling cycle skipping in fullwaveform inversion with intermediate data[END_REF]. As shown in Supplementary Section S3, they are strongly affected by our choice of c 0 . Hence this observable can be used for optimizing our wave propagation model. As shown by Supplementary Fig. S3, the value c 0 = 1500 m/s is the seismic wave velocity that clearly minimizes the level of these secondary lobes.

[|H out | 2 ⊗ |H in | 2 ] (2∆r) (Lam-
Despite this optimization, the focusing quality remains far from being ideal because of the heterogeneous distribution of c in the subsoil. In the following, we will show how this fundamental issue can become a strength since it can enlarge virtually the aperture angle under which the geophone array is seen, thereby leading to an enhanced horizontal resolution.

II.1.4 Distortion matrix

To that aim, a new operator is introduced: The so-called distortion matrix D [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2020b). This operator essentially connects each virtual geophone with the distortion exhibited by the associated wave front in the far-field. The D-matrix is thus equivalent to a reflection matrix but in a moving frame, i.e centered around each input focusing beam. This change of frame will allow us to unscramble the contribution of phase aberrations from the medium reflectivity. Last but not least, it will be shown to be particularly efficient for spatially distributed aberrations. While conventional adaptive focusing techniques are only effective over a single isoplanatic patch, the typical area over which aberrations are spatially-invariant, the D-matrix is an adequate tool to discriminate them and address them independently.

II.1.4.1 Reflection matrix in a dual basis

The reflection matrix R is first projected into the Fourier basis in reception. To that aim, we define a free-space transmission matrix T 0 which corresponds to the Fourier transform operator. Its elements link any transverse component k = (k x , k y ) of the wave vector in the Fourier space to the transverse coordinate r = (x, y) of any point in an ideal homogeneous medium:

T 0 (k, r) = exp (ik • r), (II.9)
where the symbol • stands for the scalar product between the vectors k and r. Each matrix R can now be projected in the far field at its output via the matrix product

R out = T 0 × R, (II.10)
Each column of the resulting matrix R out = [R(k out , r in )] contains the reflected wave-field in the far-field for each input focusing point r in . Injecting Eq. II.5 into the last equation yields the following expression for R out :

R out = T out × Γ × H in , (II.11)
where T out = T 0 × H out is the output transmission matrix that describe wave propagation between the focused basis and the Earth surface in the Fourier basis. In terms of matrix coefficients, the last equation can be rewritten as follows:

R out (k out , r in ) = ˆdrT out (k out , r)γ(r)H in (r, r in ).
(II.12)

In a multi-target medium made of a few bright scatterers, the reflection matrix can be leveraged to focus selectively on each scatterer. This is the principle of the DORT method [French acronym for Decomposition of the Time Reversal Operator; [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF]; [START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF]]. Mathematically, the DORT method relies on a singular value decomposition (SVD) of the reflection matrix. Physically, the singular vectors of R out are indeed shown to be the time reversal invariants of the system, i.e the wave-fronts on which would converge an iterative time reversal process, i.e a succession of time reversal operations on the reflected wave-field recorded by the array. In the single scattering regime, a one-to-one association actually exists between each eigenstate of R out and each scatterer. Each singular value is directly equal to the reflectivity γ i of the corresponding scatterer. Each output eigenvector yields the wave-front, T * out (k out , r i ), that should be applied from the far-field in order to selectively focus on each scatterers' position r i .

The DORT method is thus particularly useful for selective focusing in presence of aberrations [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF][START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF] or target detection in a multiple scattering regime [START_REF] Shahjahan | A random matrix approach to detect defects in a strongly scattering polycrystal: How the memory effect can help overcome multiple scattering[END_REF][START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]. However, this decomposition is of poor interest for diffraction-limited imaging since each input eigenvector yields an image of the scatterer, H in (r i , r in ), that is still hampered by aberrations.

II.1.4.2 Definition and physical interpretation of the distortion matrix

To cope with the fundamental issue of imaging, a novel operator, the so-called distortion matrix D, has been recently introduced [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2020b). Inspired by previous works in ultrasound imaging [START_REF] Varslot | Eigenfunction analysis of stochastic backscatter for characterization of acoustic aberration in medical ultrasound imaging[END_REF][START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF], it relies on the decomposition of the reflected wave-front into two contributions ). The principle of our approach is to isolate the latter contribution by subtracting, from the reflected wave-front, its ideal counterpart. Mathematically, this operation can be expressed as a Hadamard (element-wise) product between the normalized reflection matrix R out and T * 0 ,

D out = R out • T * 0 , (II.13)
which, in terms of matrix coefficients, yields For spatially-invariant aberrations [Eq. II.6], the D-matrix coefficients can be derived from Eqs. (II.12) and II.14:

D out (k out , r in ) = R out (k out , r in )T * 0 (k out ,
D(k out , r in ) = Hout (k out ) ˆdrγ(r + r in )H in (r)e ikout.r , (II.15)
where Hout is the aberration transmittance, that is to say the 2D Fourier transform of the output PSF H out : Hout (k out ) = ´drH out (r)e -ikout.r . Hout is the key for optimal focusing since its phase conjugate directly provides the focusing law that needs to be used to overcome the aberrations induced by the medium heterogeneities.

To extract Hout , the D-matrix is the right tool. Equation (II.15) can actually be interpreted as the result of the following fictitious process: Imagine a beam with γ as PSF and impinging onto a (fictitious) scatterer at the origin, with reflectivity distribution H in (r)

[see Fig. II.3(c)]; the resulting scattered wave-field in the direction k out would be given by Eq. (II.15) [see the analogy with Eq. (II.12)]. The problem is to isolate the aberration transmittance Hout from the scattering (γ) and input beam (H in ) terms. In that context, iterative time reversal brings a solution : it is well known that it eventually converges the wave-front that overcomes aberrations and optimally focuses on the brightest point of the (fictitious) scatterer [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF]. In the next section, we show how to implement this idea, to estimate aberration, and finally build a high-resolution image of the subsoil.

II.1.5 Time reversal analysis of the distortion matrix

in (b) and V

(2) in (c), and output eigenvectors, U

out (d) and U

(2)

out (e), are shown with a linear color scale.

The subsequent time reversal analysis consists of different steps that we will describe below. At each iteration, a virtual scatterer is synthesized at input or output through the distortion matrix concept. In the first step, a SVD of D out decomposes the field-ofview (i.e., the transverse size of the focal plane) into a set of isoplanatic patches. The corresponding eigenvectors yield an estimation of the aberration transmittance over each isoplanatic patch. Their phase conjugate provide the focusing laws that enable a (partial) compensation for the phase distortions undergone by the reflected waves during their travel between the focal plane and the geophone array. By alternatively applying the same aberration correction process at input, the size of the virtual scatterer can be gradually reduced. It converges to optimal focusing laws that will, ultimately, provide a high-resolution mapping of the SJFZ subsoil.

II.1.5.1 Output distortion matrix and isoplanatic patches

At each depth, a time reversal analysis of the distortion matrix is performed. The first step consists in a SVD of the output distortion matrix D out :

D out = U out × Σ × V † in (II.16)
or, in terms of matrix coefficients,

D(k out , r in ) = N i=1 σ i U (i) out (k out )V (i) * in (r in ).
(II.17)

Σ is a diagonal matrix containing the real positive singular values σ i in a decreasing or-

der σ 1 > σ 2 > • • • > σ N .
U out and V in are unitary matrices whose columns, U

(i) out = [U (i) out (k out )] and V (i) in = [V (i) in (r in )]
, correspond to the output and input singular vectors, respectively. For spatially-invariant aberrations, the physical meaning of this SVD can be intuitively understood by considering the asymptotic case of a point-like input focusing beam [H in (x) = δ(x)]. In this ideal case, Eq. (II.15) becomes D(k out , r in ) = Hout (u out )γ(r in ). Comparison with Eq.(II.17) shows that D out is then of rank 1 -the first output singular vector U [START_REF] Yao | Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to se tibet[END_REF] out yields the aberration transmittance Hout while the first input eigenvector V [START_REF] Yao | Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to se tibet[END_REF] in directly provides the medium reflectivity. In reality, the input PSF H in is of course far from being point-like. Moreover, aberration is not laterally invariant across the focal plane. The matrix D out is thus not singular and its spectrum displays a continuum of singular values. As shown in Supplementary Section S3, only the two first eigenstates are of interest for imaging in this specific case. In the following, we will thus restrict our study to the corresponding signal subspace.

In a medium displaying a complex wave velocity distribution, the SVD of D out can provide a decomposition of the field-of-view into several isoplanatic patches. On the one hand, each input singular vector, V (2) in focus onto two disjoint areas at depth z = 3600 m. On the other hand, each output singular vector U (p) out is linked to the corresponding aberration transmittance H(p) out [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]:

U (p) out (k out ) ∝ H(p) out (k out ) H(p) in H(p) in (k out ), (II.18)
where the symbol stands for the correlation product. However, U

out is also mod-

Û (p) out (k out ) = U (p) out (k out )/|U (p) out (k out )|. (II.19)
Indeed, if we make the realistic hypothesis of a real and positive autocorrelation function Despite this phase mismatch, the phase conjugate of Û(p) out can be used as a focusing law to compensate (at least partially) for phase distortions at the output (Fig. II.3f).

H(p) in H(p)
Updated focused reflection matrices can indeed be obtained as follows: in can indeed be expressed as follows: Two reasons can be invoked to account for that surprising result. The first possibility is an overestimation of the wave speed c at shallow depth. This value actually dictates the maximum spatial frequency exhibited by body waves induced at the Earth surface. However, the chosen value c 0 = 1500 m/s seems in agreement with P-wave velocity measured in the shallow (<100 m) fault zone at the site under study [START_REF] Share | Characterizing the uppermost 100 m structure of the San Jacinto fault zone southeast of Anza, California, through joint analysis of geological, topographic, seismic and resistivity data[END_REF]. The other explanation is that the heterogeneous subsoil acts as a lens between the geophones and the focal plane either by scattering or wave guiding. As in a time reversal experiment, scattering by small-scale heterogeneities can widen the angular aperture of the focused beam [START_REF] Derode | Robust acoustic time reversal with high-order multiple scattering[END_REF]. [START_REF] Khaidukov | Diffraction imaging by focusing-defocusing: An outlook on seismic superresolution[END_REF] also discussed that the scattering component of the wave-field holds information on small-scale subsurface heterogeneities and therefore contributes to the high resolution or even super-resolution of the migrated images. Alternatively, the local gradient of the background seismic velocity in the vicinity of the fault can constitute a wave guide through which seismic waves can be channeled [START_REF] Li | Fault zone trapped seismic waves[END_REF]. Similarly to scattering, reflections on wave guide boundaries can also enlarge the effective array aperture: This is the so-called kaleidoscopic effect [START_REF] Roux | Time-reversal in a waveguide[END_REF]. In both cases, if the induced aberrations are properly compensated, the spatial extension δ (p) in/out of the imaging PSF can be improved compared to the diffractionlimited lateral resolution δ 0 predicted for a homogeneous medium. Several reasons can be invoked to account for this striking result. The scattering events induced by the smallscale heterogeneities in the damage zone or the wave guiding induced by the gradient of

R p = ( Û(p) out • T 0 ) † × R out (II.20) New confocal images, I p (r, z) = |R p (r, r, z)| 2 ,
U (p) in (k in ) ∝ H(p) in (k in ) H(p) in (k in ) δ H(p) out δ H(p) out (k in ). (
δD (p) out = [T 0 × R p ] • T * 0 , (II.24)
The corresponding correlation matrix, δC

(p) out = δD (p) out × δD (p) †
out , is then computed and its coefficients are normalized, such that: out ∼ λ ∼ 100 m. Strikingly, the observed lateral resolution is one-sixth of the theoretical limit based on the geophone array aperture.

δ Ĉ(p) out (k out , k out ) = δC (p) out (k out , k out )/|δC
As before, the process can be iterated by exchanging the focused and Fourier basis at input and output. A novel input distortion matrix δD

II.1.6 Three-dimensional imaging of the San Jacinto Fault zone

Having demonstrated in Sec. II.1.5 how to correct phase distortions at each depth, a 3D image of the subsurface can now be uncovered. To that aim, each isoplanatic patch should be recombined at each depth. The resulting image I M is a coherent sum between the diagonal coefficients of the corrected reflection matrices R p :

I M (r, z) ≡ P p=1 R p (r, r, z) 2 .
(II.28)

In the present case, the number P of isoplanatic patches is found to be equal to 2 over the whole depth range. The fault structure, showing a different wave velocity distribution on each side of the fault, probably explains this peculiar behavior. For a more complex distribution of seismic wave velocity (such as an anticline), aberrations will also distort the coherence surface and the associated isochronous volume.

The distortion matrix method will then correct aberrations along this coherence surface but the vertical (axial) aberrations and the distortion of the coherence surface will remain unchanged. To cope with this issue, a first option is to consider the impulse responses be-tween virtual geophones located at different depths. Such a 3D reflection matrix approach will be addressed in future works. A second strategy is to include the matrix approach in any inversion scheme that aims at mapping the velocity distribution of bulk seismic waves in the underground. This would, in turn, correct the image from its axial aberrations. Indeed, a wave velocity model closer to reality and a more accurate approximation of the propagating Green's function (Eq. II.1) will always improve the final image, especially its vertical resolution.

Second, our approach relies on a scalar model of wave propagation. In the present case, this approximation is justified by the fact that we are only considering the vertical components of the computed impulse responses between geophones. As the bulk waves are mostly propagating in the vertical direction, the detected echoes are mainly associated with P waves. Nevertheless, possible shear wave conversion during wave propagation can occur. Moreover, a more contrasted image could be obtained if we were able to consider the other components of the Green's functions and also take into account the presence of shear waves in our propagation model. The method could thus be refined in the near future by taking into account both P-and S-waves as well as potential wave conversion between them induced by scattering. However, these aspects are out-of-scope for a first demonstration in the context of the seismic imaging of a fault zone.

Despite the limits of the propagation model used in our matrix approach, the 3D image of the SJFZ subsoil exhibits striking properties that are subject to physical interpretation.

In a layered medium, constant horizontal slowness implies that the structure outside of the lateral extension of the receiver area cannot be imaged. On the contrary, in This value actually dictates the maximum spatial frequency exhibited by body waves induced at the Earth surface. However, the chosen value c 0 = 1500 m/s seems in agreement with P-wave velocity measured in the shallow (<100 m) fault zone at the site under study [START_REF] Share | Characterizing the uppermost 100 m structure of the San Jacinto fault zone southeast of Anza, California, through joint analysis of geological, topographic, seismic and resistivity data[END_REF]. Alternatively, the local gradient of the background seismic velocity in the vicinity of the fault can constitute a wave-guide through which seismic waves can be channeled [START_REF] Li | Fault zone trapped seismic waves[END_REF]. Reflections on wave guide boundaries can also enlarge the effective array aperture: This is the so-called kaleidoscopic effect [START_REF] Roux | Time-reversal in a waveguide[END_REF]). As mentioned above, a last hypothesis is that the heterogeneous subsoil acts as a lens between the geophones and the focal plane either by scattering or wave guiding. As in a time reversal experiment, scattering by small-scale heterogeneities can widen the angular aperture of the focused beam [START_REF] Derode | Robust acoustic time reversal with high-order multiple scattering[END_REF]. Here the damage area is particularly heterogeneous and can play the role of scattering lens. Moreover, its location near the surface and its finite thickness (∼ 1000 m) implies the existence of an angular memory effect even for multiple scattering speckle [START_REF] Freund | Memory effects in propagation of optical waves through disordered media[END_REF][START_REF] Feng | Correlations and fluctuations of coherent wave transmission through disordered media[END_REF][START_REF] Katz | Looking around corners and through thin turbid layers in real time with scattered incoherent light[END_REF][START_REF] Katz | Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations[END_REF]. In that configuration, multiple scattering manifests itself as high-order aberrations associated with relatively small isoplanatic patches. Such highorder aberrations can be corrected by our matrix method in the plane wave basis [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]. This justifies a posteriori the choice of this basis for the correction of aberrations. Moreover, plane wave beamforming is particularly adequate in a multi-layered medium since aberrations are laterally-invariant in that frame. Beyond the specific case of SJFZ, note that the choice of basis for the aberration correction is flexible (Lambert et al., 2021b). The reflection matrix can be ideally projected onto any aberrating layer in the subsoil. This choice shall be dictated by the local topography and any prior knowledge on the local distribution of seismic wave velocities in the zone under study.

Multiple scattering and/or wave-guiding effects could also explain the optimal velocity (c 0 = 1500 m/s) found for the wave propagation model. Initially, this choice was justified by the focusing quality and image resolution reached at the end of the matrix imaging process (see Supplementary Figs. S3 andS4). Nevertheless, a physical interpretation can now be provided to account for it. Our hypothesis is that a weak wave velocity model can enable the time gating of multiply-scattered or guided waves that are, ultimately, transmitted until the focal plane where they can be efficiently scattered by subsoil heterogeneities. Indeed, such distorted paths are associated with larger echo times than the ballistic waves usually considered by reflection imaging methods. As mentioned above, they can enlarge the angular aperture of our imaging system and account for a transverse resolution much better than what would be expected if we had used direct ballistic waves. In other words, the use of a relatively weak wave velocity allows us to adapt a rough homogeneous model to a medium that displays a strongly heterogeneous wave velocity distribution. This remains of course an hypothesis and a further analysis will be required in order to prove rigorously the origin for the transverse super-resolution in the context of fault seismic imaging.

II.1.8 Conclusion

Inspired by pioneering works in optical microsopy [START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF][START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF] and ultrasound imaging [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF]Lambert et al., 2020b), a novel matrix approach to seismic imaging is proposed in this paper. Taking advantage of the reflection of bulk seismic waves by heterogeneities in depth, it can be applied to both active or passive seismic imaging. The strength of this approach lies in the fact that it works even when the velocity distribution of the subsoil is unknown. By projecting the seismic data either in the focused basis or in the Fourier plane, it takes full advantage of all the information contained in the collected data. For aberration correction, projection of the reflection matrix into a dual basis allows the isolation of the distorted component of the reflected wavefield. Seen from the focused basis, building the distortion matrix D consists in virtually shifting all the input or output focal spots at the onto the same virtual location. An iterative time reversal analysis then allows to unscramble the phase distortions undergone by the wave-front during its travel between the Earth surface and the focal plane. A one-to-one association is actually found between each eigenstate of D and each isoplanatic patch in the field-of-view. More precisely, each singular vector in the Fourier space yields the farfield focusing law required to focus onto any point of the corresponding isoplanatic patch.

A confocal image of the subsoil reflectivity can then be retrieved as if the underground had been made homogeneous.

In this paper, as a proof-of-concept, the case of the San Jacinto Fault zone is considered. While a raw confocal image suffers from an extremely bad horizontal resolution due to the strong lateral variations of the seismic velocities in the vicinity of the fault, our matrix approach provides a high-resolution image. Strikingly, its lateral resolution is almost one eighth below the diffraction limit imposed by the geophone array aperture. This surprising property may be accounted for by the heterogeneities of the subsoil that can play the role of a scattering and/or channeling lens which increases drastically the effective array aperture. The matrix image reveals a damage volume particularly pronounced in the shallow crust (<1000 m). At larger depth, the 3D image of the fault exhibits the fault blocks on the right and left side of the slip interface. A more detailed interpretation of the obtained image will be the focus of a future study.

Besides seismic fault zones, a matrix approach of passive imaging is particularly suited to the study of volcanoes. While the multiple scattering problem in volcanoes has been recently tackled under the reflection matrix approach [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF], the distortion matrix would be particularly powerful to restore a diffraction-limited image of the in-depth structure of a volcano. Additionally, our matrix method could advantageously be applied to real-time geophysical imaging, similarly to adaptive optics in astronomy.

This may provide a valuable tool for monitoring oil or gas reservoirs when drilling, extracting hydrocarbon or injecting CO 2 . However, more research is needed to see whether the current method would succeed in these complex areas that exhibit such strong variations of the background wave velocity. Finally, matrix imaging can be taken advantage of for marine exploration where the variations of sea local properties can distort the acoustic wave-front depending on the gradients of temperature and salinity, not to mention streams and swirls. The matrix approach presented here can compensate for the phase distortions undergone by the received echoes and improve the interpretation of the data with no more required knowledge than the speed of sound in water. Although not leveraged in the present paper, Lambert et al. (2020b) have also demonstrated the relevance of the matrix approach at removing multiple reverberations between strata layer interfaces. This also may be useful for any marine or seismic exploration purpose. by geophones located at u = u out and u = u in . Under appropriate wave-field conditions, coda cross-correlation converges towards the Green's function G(u out , u in , t) between receiving stations as if one of them (u in ) had become a source:

II.1.9 Supplementary Material

K(u out , u in , t) = lim T →∞ 1 T ˆT 0 dτ ψ(u in , τ )ψ * (u out , t + τ ) = G(u out , u in , t) -G(u out , u in , -t) (S1)
For the special case of u in ≡ u out , the autocorrelation signal K(u in , u in , τ ) gives rise to an intense peak at lag time t = 0, that physically corresponds to the seismic noise power spectral density measured at point u in . Only the non zero lag time contribution carries information on the reflectivity at depth and is thus of interest for imaging purposes. However, this zero time peak gathers most of the energy content in the retrieved signal, and can have a detrimental impact on our analysis. In the present situation, the limited frequency bandwidth makes this initial pulse duration far from being negligible (δt ∼ 1/∆f = 0.1 s). It also gives rise to strong side lobes that pollute the relevant signal in the full time range (see Fig. S1a). Additionally, the ambient seismic wave exhibits a coherence length close to λ/2. Given the high density of the geophones network, neighbour stations thus belong to the same coherence area. Therefore the corresponding cross-correlation signals,

that lie close to the diagonal of K, are also dominated by this autocorrelation peak.

To get rid of this central pulse, a gaussian mask is applied to each element of the raw reflection matrix

K (u out , u in , t) = K(u out , u in , t) × 1 -exp - u out -u in 2 2λ 2
where λ is the wave-length of body waves at the central frequency (here λ = 100 m).

This filter gradually penalizes the impulse responses computed for interstation distances smaller than λ/2, the coherence length of the body waves in the 10 Hz-20 Hz working frequency range.

Figure S1b displays the filtered version K of the measured reflection matrix K shown in Fig. S1a. The filtered wave-field is no longer dominated by the seismic noise autocorrelation. Note the strong gap in terms of order of magnitude between the diagonal coefficients of K (Fig. S1a) and the actual wave-field patterns in K (Fig. S1b). This pre-filtering operation is thus critical for an accurate interpretation of the effective Earth's response.

Each column of K corresponds to the wave-field that would be recorded by the set of geophones if a pulse were emitted from one geophone at u in . Figure S1c displays the wave-field generated from the central geophone position and interpolated at the surface between all receiving geophones. This wave-field clearly exhibits the contribution of direct Rayleigh waves that emerges along the white dashed circle. Interestingly, we observe fast travelling waves (see for instance the top left-hand corner of the array) that correspond to the signature of the body waves propagating below the Earth's surface at larger speed than the surface waves tracked by the geophones. They overtake the surface waves and eventually get first to the edge of the array, after reflection on some shallow structures or refraction at overcritical angles.

S2 Singular value decomposition of the distortion matrix

To decompose the field-of-view into several isoplanatic patches, a SVD of D out is performed. Mathematically, this operation is equivalent to the eigenvalue decomposition of the correlation matrix C = D out × D † out :

C = U out × Σ 2 × U † out (S2)
where the eigenvalues of C are the square norm of the singular values σ i of D out and its eigenvectors are the output singular vectors

U (i)
out . In presence of multiple isoplanatic patches, C can be decomposed into a set of sub-matrices C (p) , such that C = p C (p) .

Each matrix C (p) is associated with a distinct isoplanatic patch p in the field-of-view. In each of them, phase distortions can be modelled by: (i) a spatially-invariant input PSF,

H (p) in (r, r ) = H (p)
in (r -r ), in the focused basis; (ii) a far-field phase screen of transmittance

H(p) out = [ H(p) out (k out )].
Using Eq. II.18, the coefficients of C (p) can be expressed as follows: out would be of rank 1. The corresponding eigenvalue σ 2 1 would then be equal to ρ (p) and the associated eigenvector U (p,1) out would directly yield the aberration transmittance H(p) out . However, in practice, the virtual scatterer is of finite size and the correlation term in Eq. ( S3) is not negligible. The rank of C (p) out then scales as the number of transverse resolution cells mapping the virtual scatterer [START_REF] Robert | The prolate spheroidal wave functions as invariants of the time reversal operator for an extended scatterer in the Fraunhofer approximation[END_REF][START_REF] Aubry | Gaussian beams and legendre polynomials as invariants of the time reversal operator for a large rigid cylinder[END_REF]: out is of interest since it maximizes the backscattered energy by focusing at the center of the virtual scatterer (see Supplementary Section S3).

C (p) out (k, k ) = ρ (p) H(p) out (k) H(p) * out (k ) H(p) in H(p) in (k -k ), (
Q = δ (0) in /δ 0 2 . Q is typically equal to 4 for the local PSF H (p) in 2 displayed in
In the following, we thus expect an overall distortion matrix D out of rank P × Q with the following singular value distribution: (i) P largest singular values associated with the main eigenstates U (p) out of each patch p; (ii) a set of (P -1) × Q of smaller singular values associated with secondary eigenstates that focus on the edges of the virtual scatterers. Figure II.4a confirms this prediction by displaying the normalized singular values, σi = σ i / N j=1 σ 2 j , of D out at depth z = 3600 m (t B = 4.8 s). As expected, a few singular values potentially associated with the signal subspace seem to predominate over a continuum of eigenvalues characteristic of the noise subspace. However, it is difficult to determine the effective rank of the signal subspace. This issue can be circumvented by computing the Shannon entropy H of the singular values [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2020b), such that

H = - N i=1 σ2 i log 2 σ2 i . (S4)
The Shannon entropy can be used as an indicator of the rank of the signal subspace.

In the case of Fig II .4, the singular values of D out display an entropy H 9.2. As the effective rank of D out scales as P × Q and that Q ∼ 4, this means that the number P of isoplanatic patches is roughly equal to 2. Hence, only the two first eigenstates of D out should be considered in the present case. As shown below, the higher-order eigenstates do not bring any supplementary information on aberrations and are shown to be useless for any imaging purpose.

S3 High-order eigenstates of the D-matrix.

In this section, the higher-order eigenstates of D out are investigated theoretically and shown to be higher-order modes of the two first ones.

In

in (Fig. S2a), seems associated with the first isoplanatic patch (see the comparison with V

(1) in in Fig. II.4b), the fourth one, V (4) in (Fig. S2b), is also linked to the second isoplanatic patch (see the comparison with V In this Supplementary Section, we explain how these higher-order eigenstates of D out are linked to the first-order ones and why they are not useful for an imaging purpose.

To that aim, we will first study theoretically the SVD of D , mapping the virtual reflector [START_REF] Robert | The prolate spheroidal wave functions as invariants of the time reversal operator for an extended scatterer in the Fraunhofer approximation[END_REF].

The corresponding eigenvectors U (i)

out can be expressed as the Hadamard product between the far-field aberration phase law H (p) out and eigenmodes

W (i) in of the correlation kernel, H(p) in H(p) in = [ H(p) in H(p) in (k -k )], such that U (i) out = H(p) out • W (i) in (S5)
The shape of the eigenmodes

W (i)
in depends on the exact form of the correlation function

H(p) in H(p) in .
For instance, a sinc correlation function imply 3D prolate spheroidal eigenfunctions [START_REF] Robert | The prolate spheroidal wave functions as invariants of the time reversal operator for an extended scatterer in the Fraunhofer approximation[END_REF]; a Gaussian covariance function leads to Hermite-Gaussian eigenmodes [START_REF] Aubry | Gaussian beams and legendre polynomials as invariants of the time reversal operator for a large rigid cylinder[END_REF]. The first eigenmode is given, in first approximation, by [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]:

W (1) in (k) = H(p) in H(p) in (k) (S6)
As the correlation function

H(p) in H(p) in (k)
is, in first approximation, real and positive, a general trend is that the first eigenmode W [START_REF] Yao | Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to se tibet[END_REF] in (k) displays a nearly constant phase. This is a very important property since it means that the phase of the first eigenvector U [START_REF] Yao | Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to se tibet[END_REF] out is a direct estimator of H(p) out . The higher rank eigenvectors W (i)

in are more complex and exhibit a number of lobes that scales with their rank i. The corresponding eigenvectors U (i) out do not bring a priori any useful information compared to the fundamental one U

(1) out as their far-field phase law is modulated by eigenfunctions W (i) in (k) that display a number (i -1) of sign reversals. As a consequence, they cannot be used for any imaging purpose since the associated PSF will be made of i lobes [START_REF] Robert | The prolate spheroidal wave functions as invariants of the time reversal operator for an extended scatterer in the Fraunhofer approximation[END_REF][START_REF] Aubry | Gaussian beams and legendre polynomials as invariants of the time reversal operator for a large rigid cylinder[END_REF].

On the contrary the PSF associated with U [START_REF] Yao | Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to se tibet[END_REF] out shows a single central lobe. In the general case (several isoplanatic patches in the field-of-view), the rank of D out and C out scales as P × Q, with P the number of isoplanatic patches. By estimating this rank from the Shannon entropy of the singular values (Eq. S4) and the number Q from the imaging PSF (Fig. II.2b), we can estimate the number of significant isoplanatic patches (here roughly equal to 2). We thus expect that the two first eigenstates of D out and C out are associated with distinct isoplanatic patches and that higher-rank eigenstates are just higher-order eigenmodes derived from the two fundamental eigenstates.

Nevertheless, from the singular value spectrum displayed in Fig. II.4a, it is difficult to state that the third and fourth eigenstates should not be considered. Their singular values clearly emerge from the continuum of lowest singular values that can be seen as a "noise" background. A first check consists in looking at the input eigenvectors V in , are shown to focus on the same areas but with a more complex pattern (see Fig. S3a andb, respectively). To show the link between these high-order eigenstates and the first-order ones, one can investigate the associated output eigenvectors U (i) out (see Fig. S3c andd). Indeed, if they are associated with the same aberration phase law H(p) out , this product should be real and exhibits a sign reversal characteristic of the second order mode W This detailed analysis confirms that the third and fourth eigenstates correspond to second-order modes of the first ones. Hence they do not bring any supplementary information on subsoil-induced aberrations, hence they should not be considered for matrix imaging. Beyond this example at this specific time-of-flight, it should be noted that the same behavior is observed over the whole depth range. Two main isoplanatic patches dominate the singular value spectrum. The fault structure probably explains this peculiar behavior: The two different wave velocity distribution on each side of the fault implies the presence of two main isoplanatic patches in depth.

S4 Normalized correlation matrix

In the last part of the matrix imaging process (Section 4.3), an output residual distortion matrix δD (p) out is considered. By analogy with Eq. S3, the corresponding correlation matrix, δC

(p) out = δD (p) out × δD (p) †
out , can be expressed as follows:

δC (p) out (k out , k out ) ∝ ρ p δ H(p) out (k out )δ H(p) * out (k out ) δ H(p) in δ H(p) in (k out -k out ). (S7)
As previously highlighted in Sec. II.1.5.1 [Eq. ( S3)], the correlation term in Eq. (S7) prevents a proper estimation of the aberration phase transmittance over the whole angular spectrum. To circumvent that issue, the correlation matrix coefficients can be normalized [Eq. ( 25)]. As illustrated by Fig. 3e, this operation makes the virtual reflector pointlike (Lambert et al., 2020b). Indeed, if we make the realistic hypothesis of a real and positive autocorrelation function δ

H(p) in δ H(p)
in in Eq. (S7), the coefficients of δC (p) out can actually be expressed as follows:

δC (p) out (k out , k out ) ∝ δ H(p) out (k out )δ H(p) * out (k out ). (S8)
Such a matrix is equivalent to the time reversal operator associated with a point-like reflector at the origin [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF][START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF]. In that case, the matrix δ Ĉ(p) out is of rank 1 and its eigenvector δU

(p)
out yields the residual aberration phase transmittance:

δU (p) out (k out ) = δ H(p) out (k out ). (S9)
The phase conjugate of δU

(p)
out yields the additional focusing law that is applied to the reflection matrix in order to overcome the residual output aberrations [Eq. ( 27)]. S5 Choice of the wave velocity model.

∆ / y δ 0 ∆ / y δ 0 ∆ / x δ 0 ∆ / x δ 0 ∆ / x δ 0 ∆ / x δ 0 ∆ / x δ 0 ∆ / x δ 0 ∆ / x δ 0 ∆ / x δ 0 ∆ / x δ 0 ∆ / y δ 0 ∆ / y δ 0 ∆ / y δ 0 ∆ / y δ 0 ∆ / y δ 0 ∆ / y δ 0
In the accompanying paper, the reflected body waves are used to build a high-resolution in-depth image of the SJFZ. To that aim, the choice of the initial wave velocity model is crucial. As mentioned in the accompanying paper, this choice can be enlightened by the properties of the focused reflection matrix R. In particular, a common mid-point intensity profile can be extracted from the main antidiagonal of R that exhibits the maximum confocal (diagonal) signal. Figure S3 shows the corresponding intensity profile at several times of flight (t B = 1.05, 1.5 and 2.65 s, from top to bottom) and for different seismic wave velocity (c 0 = 1000, 1500 and 2000 m/s, from left to right) in our propagation model. To make the comparison quantitative, each intensity profile is displayed as a function of spatial coordinates normalized by δ 0 , the expected resolution for each velocity model at the corresponding depth. In each panel, the focal spot consists of a main lobe centered around the input focusing point and a random distribution of secondary lobes.

The dimension of the main lobe is roughly twice the diffraction limit prediction δ 0 , depicted by a white circle in each panel Fig. S3. This loss of resolution is a manifestation of the aberrations induced by the gap between the velocity model and the actual seismic wave velocity distribution in the SJFZ underground. However, the dimension of this main lobe is not significantly affected by our choice of c 0 . Hence this observable cannot be used for optimizing our wave propagation model. On the contrary, the amplitude and spatial extension of the secondary lobes strongly depend on the wave velocity model. The value c 0 = 1500 m/s is the seismic wave velocity that clearly minimizes the level of these secondary lobes in Fig. S3.

The choice is validated a posteriori by the quality of the confocal image obtained under our matrix approach. For c 0 = 2000 m/s (Fig. S4c), the strata structure of the SJFZ underground is partially revealed but with a worse resolution and lower contrast than for c 0 = 1500 m/s (Fig. S4b) until an echo time t=4 s. Nevertheless, note that, beyond that time, the two last interfaces seem to be better resolved for c 0 = 2000 m/s. This indicates that, not surpris- ingly, the seismic wave velocity increases with depth and that c 0 = 1500 m/s is probably not the optimal wave velocity for t > 4 s. A mapping of the wave velocity could be indeed possible through the focused reflection matrix approach as already demonstrated in ultrasound imaging (Lambert et al., 2020a). However, the extension of this method to seismology is beyond the scope of this paper. The matrix mapping of the bulk seismic wave velocity will be tackled in future works. and reflectivity of scatterers for seismic waves starting with a simple homogeneous background velocity model of the medium. The derived three-dimensional image of the fault zone resolves lateral variations of scattering properties in the region within and around the surface fault traces, as well as differences between the Northwest (NW) and the Southeast (SE) parts of the study area. A localized intense damage zone at depth is observed in the SE section, suggesting that a geometrical complexity of the fault zone at depth induces ongoing generation of rock damage.

II.2.2 Introduction

Earthquakes are among the most destructive natural disasters. Although earthquakes are generally unpredictable, some aspects of their behavior such as the likelihood of being arrested and statistically-preferred propagation direction can be estimated from structural properties of fault zones (see e.g. [START_REF] Chester | Internal structure and weakening mechanisms of the san andreas fault[END_REF]; [START_REF] Wesnousky | Seismological and structural evolution of strike-slip faults[END_REF]; Ben-Zion ( 2008)). Fault zones that are the structural manifestation of earthquakes evolve during deformation and have generally complex properties (e.g. [START_REF] Mitchell | The nature and origin of off-fault damage surrounding strikeslip fault zones with a wide range of displacements: A field study from the atacama fault system, northern chile[END_REF]; [START_REF] Ben-Zion | Characterization of fault zones[END_REF]). Characterizing the geometrical and seismic properties of fault zones can provide important information for assessing likely past and future rupture properties. Fault zones have also strong impact on fluid flow in the lithosphere [START_REF] Knipe | Faulting, fault sealing and fluid flow in hydrocarbon reservoirs: an introduction[END_REF].

Fault zones are manifested at the surface by several main fault traces that accommodate the bulk of the long term slip. They are characterized by lineaments, topography and various geometrical complexities. The identification of fault traces is done by field observations (major line of fracturing, offset in geological units), along with remote sensing techniques that analyze ground deformation after major earthquakes obtained from satellites and aircrafts such as Interferometric Synthetic Aperture Radar (InSAR, [START_REF] Massonnet | The displacement field of the landers earthquake mapped by radar interferometry[END_REF]) and subpixel correlation of optical images (SPOT, [START_REF] Binet | Horizontal coseismic deformation of the 2003 bam (iran) earthquake measured from spot-5 thr satellite imagery[END_REF]). Fault zone properties below the surface are obtained by seismic and other geophysical imaging techniques.

Fault zones have hierarchical damage structures that evolve during the fault zone activ-ity and have several general elements (e.g. [START_REF] Chester | Internal structure and weakening mechanisms of the san andreas fault[END_REF]; [START_REF] Rockwell | High localization of primary slip zones in large earthquakes from paleoseismic trenches: Observations and implications for earthquake physics[END_REF]). The principle slip zone is a highly localized thin layer (0.01 -0.1 m thick) that accommodates most of the fault slip and is characterized by ultra cataclasite rock particles. The principle slip zone is bounded by a core damage zone (inner damage zone) that is typically about 100 m wide and asymmetrically located on one side of the slip zone of large faults [START_REF] Lewis | Shallow seismic trapping structure in the san jacinto fault zone near anza, california[END_REF][START_REF] Dor | Geological observations of damage asymmetry in the structure of the san jacinto, san andreas and punchbowl faults in southern cali-fornia: A possible indicator for preferred rupture propagation direction[END_REF]. The core damage is surrounded by a broader zone of reduced damage intensity (referred to as outer damage zone) that may extend for several km on each side of the fault. Properties of the fault zone damage provide information on statistical tendencies of local earthquake ruptures, operating dynamic stress field, energy dissipation and more (e.g. [START_REF] Manighetti | Evidence for self-similar, triangular slip distributions on earthquakes: Implications for earthquake and fault mechanics[END_REF]; [START_REF] Mitchell | The nature and origin of off-fault damage surrounding strikeslip fault zones with a wide range of displacements: A field study from the atacama fault system, northern chile[END_REF]; [START_REF] Xu | Properties of inelastic yielding zones generated by in-plane dynamic ruptures-ii. detailed parameter-space study[END_REF]). For that reason, a number of seismic and other methods have been developed to provide detailed information on fault zone structures.

Among the seismic imaging techniques, reflection seismology is generally pertinent to image planar horizontal layers with a very high resolution, and provides indirect imaging of faults by the offset of sedimentary layers. It relies on the analysis of seismic waves that are sent back towards the surface after being reflected or scattered by subsurface structures with strong impedance contrasts. The recorded wave-field is composed of reflected waves that result from the interaction of seismic waves with planar reflectors such as layer boundaries, and diffracted waves from small-scale geological objects such as cracks and fractures. To obtain structural information on the subsurface, migration techniques are applied aiming mainly to relocate reflectors and scatterers in depth or in time (see, e.g., review by [START_REF] Etgen | An overview of depth imaging in exploration geophysics[END_REF]. Migration of seismic wave-fields requires an accurate velocity model of the Earth. Errors and biases in the velocity model produce artefacts and defocusing due to phase distortions in the migrated images [START_REF] Moser | Diffraction imaging in depth[END_REF].

Reflection information is often gathered from seismic surveys where seismic energy propagating into the medium is generated by man-made sources such as vibrators, explosives, etc. In the last decade or so, passive methods based on the ambient seismic noise have been developed to substitute active imaging techniques (e.g. [START_REF] Campillo | Seismic imaging and monitoring with ambient noise correlations[END_REF]. Cross-correlation of passive traces recorded at two receivers allows retrieving the Green's function between these two receivers. In other words, the resulting correlation is comparable to the seismogram that would be obtained at the first receiver if there is a source located at the second receiver's location. This approach is referred to as seismic interferometry (Wapenaar et al., 2010a). Noise cross-correlations are often used to retrieve the surface waves component of the Green's function (Shapiro and Campillo, 2004, and later works). Extracting body waves contributions to the Green's function is more difficult (Poli et al., 2012a). However, it has been shown that ambient noise cross-correlations can be used to image deep targets inside the Earth with body wave reflections [START_REF] Draganov | Retrieval of reflections from seismic background-noise measurements[END_REF]Poli et al., 2012b). Inspired by works done in ultrasound imaging (Aubry and Derode, 2009b) and optical microscopy [START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF], a reflection matrix approach was introduced to geophysics and used body wave reflection from coda-wave cross-correlations to image the complex medium below Erebus volcano -Antarctica [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]. The reflection matrix approach allows to image both specular reflectors, corresponding to the boundaries between layers of different propagation velocities, and non-specular reflections generated by a distribution of localized inhomogeneities. The matrix imaging approach has then been developed to overcome phase distortions for multi-layered media (Lambert et al., 2020a) and strongly heterogeneous media [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2021a). Using the same matrix formalism, [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF] analysed ambient noise recorded at a dense array (1108 vertical geophones in 600 m × 700 m configuration) to image subsurface properties of the San Jacinto Fault zone (SJFZ) southeast of Anza, California. While [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF] dealt with imaging problems in the multiple scattering regime in the case of volcanoes, the main challenge in fault zones imaging is the presence of phase distortions, also referred to as aberrations, due to the strong structural heterogeneities within and around the fault zones.

The present paper follows the approach of [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF] to derive more detailed 3D images of the SJFZ at depth and to resolve the shallow structure with unprecedented resolution. [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF] computed noise cross-correlations in the (10-20) Hz frequency range. Whereas the cross-correlations provide a response matrix between sources and receivers located at the surface, the reflection matrix contains the response between sources and receivers that are virtually moved inside the medium by performing focus-ing operations. This process is generally known as redatuming [START_REF] Berkhout | A unified approach to acoustical reflection imaging. II: The inverse problem[END_REF] and it allows local information on the medium's reflectivity to be retrieved. To project the data in the virtual focused basis, an homogeneous transmission matrix is used with a constant velocity model of 1500 m/s. This velocity was chosen to optimize the focusing as discussed in [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. The main advantage of this method is that it only requires an approximate estimation of the medium's velocity. Although an incorrect velocity model produces phase distortions (aberrations) in the propagated data, the reflection matrix approach allows us to account and correct these aberrations through the distortion matrix concept [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2020bLambert et al., , 2021b)). From that matrix, the distorted component is extracted and is used to virtually focus back waves inside the medium (see Appendix A). As a result of the correction process, the resolution of the final subsurface images is drastically improved and a three-dimensional image of the subsurface reflectivity is revealed. The derived images represent reflectivity maps of the medium beneath the Clark branch of the SJFZ at the Sage Brush Flat site [START_REF] Ben-Zion | Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone[END_REF][START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF]. The site under study is located in the complex trifurcation area southeast of Anza, California (Fig. II.9a). The locations of the surface fault traces are derived from recent detailed studies of the surface geological mapping and shallow geophysical imaging [START_REF] Wade | Geologic and structural characterization of shallow seismic properties along the san jacinto fault at sage brush flat, southern california[END_REF][START_REF] Share | Characterizing the uppermost 100 m structure of the San Jacinto fault zone southeast of Anza, California, through joint analysis of geological, topographic, seismic and resistivity data[END_REF]. The basic goal of this paper is to interpret the obtained scattering images that exhibit features with higher lateral resolution than conventional seismic investigations. The variability and attenuation of scattered intensity within and around the major fault zone are also discussed.

II.2.3 The San Jacinto Fault Zone

The 230 km-long San Jacinto Fault Zone is the most seismically active fault zone in southern California [START_REF] Hauksson | Waveform relocated earthquake catalog for southern california (1981 to june 2011)[END_REF] and is one of several major right-lateral strike-slip fault zones over which the North American-Pacific plate boundary is distributed in southern California. The SJFZ branches from the San Andreas fault at Cajon Pass and was formed 1 -2 million years ago, presumably in response to the geometrical complexities on the San Andreas Fault in the transverse ranges (e.g. [START_REF] Matti | Extension and contraction within an evolving divergent strikeslip fault complex: The san andreas and san jacinto fault zones at their convergence in southern california. The San Andreas fault system: Displacement, palinspastic reconstruction, and geologic evolution[END_REF]). The SJFZ represents a less mature evolutionary stage in the life of a large continental strikeslip structure than the San Andreas fault. Approximately 24 km of slip has been accommodated by the SJFZ since the latest Pliocene to early Pleistocene [START_REF] Dorsey | Quaternary landscape evolution in the san jacinto fault zone, peninsular ranges of southern california: Transient response to strike-slip fault initiation[END_REF], with estimated slip rates that vary along strike between 8 -20 mm/yr [START_REF] Rockwell | 3,000 years of groundrupturing earthquakes in the anza seismic gap, san jacinto fault, southern california: Time to shake it up[END_REF][START_REF] Fialko | Interseismic strain accumulation and the earthquake potential on the southern san andreas fault system[END_REF]. The SJFZ has varying surface complexity and seismicity along its strike. The Anza section to the northwest of the SGB site consists of a single strand, the Clark fault, with relatively regular geometry and low current background microseismicity [START_REF] Sanders | A seismotectonic analysis of the anza seismic gap, san jacinto fault zone, southern california[END_REF]. The trifurcation area of the SJFZ where the SGB site is located (Fig

II.2.4 3D scattering volume

Migration techniques are known to be powerful tools for imaging strong reflecting boundaries. These boundaries are identified by discontinuities of acoustic impedance in the subsurface and are characterized by specular returns in seismic records. Less interest has been accorded to the non-specular component that arises from small-scale geological objects [START_REF] Khaidukov | Diffraction imaging by focusing-defocusing: An outlook on seismic superresolution[END_REF]. The energy generated by such small objects is commonly referred to as diffractions. Non-specular energy holds valuable information on the local heterogeneities in the medium (Schwarz, 2019a, and references therein). While many conventional migration techniques, such as Gaussian beam migration and reverse-time migration, honor the diffracted component, such contributions are generally difficult to analyze: they are often suppressed due to the processing done in conventional seismic methods or masked by specular reflections whose amplitudes are much larger than the scattered components [START_REF] Kozlov | Imaging scattering objects masked by specular reflections[END_REF]. They can also be considered as noise in several migration techniques. Keeping the non-specular component in the analysis allows retrieving signatures of localized scatterers such as cracks or inclusions that lack lateral continuity. Imaging such features whose size is of the order and even smaller than the seismic wavelength contributes significantly to seismic interpretation [START_REF] Schwarz | Coherent diffraction imaging for enhanced fault and fracture network characterization[END_REF].

Fault surfaces and zones with increased fractures density are non-specular objects for surface sensors [START_REF] Kanasewich | Imaging discontinuities on seismic sections[END_REF]. Imaging such features is a challenge in most conventional seismic exploration surveys. Several studies have discussed the necessity of distinguishing between diffractions and reflections, and provided techniques to separate them [START_REF] Khaidukov | Diffraction imaging by focusing-defocusing: An outlook on seismic superresolution[END_REF][START_REF] Moser | Diffraction imaging in depth[END_REF][START_REF] Bakhtiari Rad | Common-reflection-surfacebased prestack diffraction separation and imaging[END_REF]. Diffraction imaging is performed usually by suppressing the specular reflections so that the migrated image contains the diffracting component that have been isolated. The reflection matrix technique allows us to image, without performing any prior filtering, not only specular reflectors but mainly non-specular backscattered energy that directly gives insight into rock properties at the subsurface. The contribution of specular and non-specular features is distinguishable in the derived images.

The constant velocity model of 1500 m/s chosen for the redatuming operation in the initial study of [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF] is highly approximate especially for the deep structure.

It was chosen to optimize the focusing at depth with an "apparent" velocity that increases the effective aperture of the geophone array to exploit the actual contribution of multiscattered paths as discussed in [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. We recall that a higher velocity will only stretch the detected features vertically. The images are presented as a function of depth for the chosen background velocity and two-way travel time.

In the subsequent subsections, we provide a detailed description of cross-sections will lead to enlightening interpretations, they should be also taken with caution. First, these images have been corrected from transverse aberrations but axial aberrations can subsist. To cope with this issue, a more accurate 2D velocity model can be implemented in the future. Secondly, the aberration correction process applied to these images only identifies two isoplanatic patches at each depth. This is probably a perfectible point in such complex heterogeneous media with strong lateral velocity variations. Other less reflective features belonging to other isoplanatic patches are thus possibly not revealed by our method in Fig. II.10 and II.11.

A compensation method for attenuation, described in C, is applied to the reflectivity maps. The results are plotted on a Cartesian grid, with the origin located at the center of the array and the x-axis orthogonal to the fault traces. The colors represent the backscattered intensity plotted in dB. The local distribution of heterogeneities and discontinuities of material boundaries are revealed with maximal focused intensity. We first present results associated with shallow materials and then discuss deeper structures.

We compared our results with previous studies investigating the same region. Most of the studies conducted at the SGB site utilize surface waves to resolve shallow features and the velocity distribution below the dense array [START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF][START_REF] Hillers | Focal spot imaging based on zero lag cross-correlation amplitude fields: Application to dense array data at the San Jacinto fault zone[END_REF][START_REF] Mordret | Shallow three-dimensional structure of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic array[END_REF] or fault zone head waves to image lateral variations of lithology and damage across the fault [START_REF] Ben-Zion | Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone[END_REF][START_REF] Qin | Internal structure of the San Jacinto fault zone in the trifurcation area southeast of Anza, California, from data of dense seismic arrays[END_REF]. To our knowledge, the deep scattering structure has not been examined and analyzed yet.

II.2.4.1 Images of the shallow fault zone

The high frequency cross-correlations of noise recorded by the dense array allow resolving features in the top few hundred meters of the crust. 2015) and [START_REF] Qin | Internal structure of the San Jacinto fault zone in the trifurcation area southeast of Anza, California, from data of dense seismic arrays[END_REF]. This zone is characterized by significant low velocities and an intense localized damage producing reflections [START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF][START_REF] Hillers | Focal spot imaging based on zero lag cross-correlation amplitude fields: Application to dense array data at the San Jacinto fault zone[END_REF][START_REF] Mordret | Shallow three-dimensional structure of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic array[END_REF]. The reflective layer SW of the fault observed in Fig II .10c coincides with the local sedimentary basin reported by [START_REF] Ben-Zion | Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone[END_REF], [START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF] and [START_REF] Hillers | Focal spot imaging based on zero lag cross-correlation amplitude fields: Application to dense array data at the San Jacinto fault zone[END_REF].

Many studies of the San Jacinto Fault Zone observed an asymmetric rock damage across the fault [START_REF] Lewis | Shallow seismic trapping structure in the san jacinto fault zone near anza, california[END_REF][START_REF] Dor | Geological observations of damage asymmetry in the structure of the san jacinto, san andreas and punchbowl faults in southern cali-fornia: A possible indicator for preferred rupture propagation direction[END_REF][START_REF] Qin | Internal structure of the San Jacinto fault zone in the trifurcation area southeast of Anza, California, from data of dense seismic arrays[END_REF][START_REF] Wade | Geologic and structural characterization of shallow seismic properties along the san jacinto fault at sage brush flat, southern california[END_REF]. The These observations are also in agreement with (Sharp, 1967, Fig.3) where the SW of the SJFZ at Table Mountain exhibits more complexity than the NE region. The scattering zones in both sections are dipping toward the SE where the fault zone is more localized.

In the following, we will show that the reflection matrix contains much more information than the medium's reflectivity. It also provides a direct insight into the variations of the scattering properties and the attenuation of scattered waves in the fault region.

II.2.5 Lateral variations of intensity

Seismic waves propagating inside the Earth give direct insight on the nature and properties of rocks. While travelling through complex heterogeneous media, waves suffer from seismic attenuation. In this section we briefly discuss seismic wave attenuation principles and how attenuation is accounted for in our matrix formalism, and therefore in the images obtained.

Seismic attenuation describes the decay of energy experienced by seismic waves while they propagate. Amplitudes are easily altered by several factors such as geometrical spreading, scattering and absorption (intrinsic or anelastic attenuation) [START_REF] Shapiro | Seismic attenuation by scattering: theory and numerical results[END_REF][START_REF] Knipe | Faulting, fault sealing and fluid flow in hydrocarbon reservoirs: an introduction[END_REF]. Evaluation of the attenuation due to scattering and intrinsic absorption has been the subject of considerable studies in seismology [START_REF] Aki | Analysis of the seismic coda of local earthquakes as scattered waves[END_REF][START_REF] Sato | Seismic wave propagation and scattering in the heterogeneous earth[END_REF]. Estimating attenuation properties provides complementary information to seismic velocity distribution, and can be particularly useful in fault zone studies to obtain a better

II.2.6 Discussion and conclusions

The presented results provide detailed images of seismic properties in the 3D volume around the San Jacinto fault at the SGB site. We used one month of ambient seismic noise recorded by a dense array deployed at SGB site around the Clark branch of SJFZ.

The high frequency seismic data provided by the spatially dense array allows us to resolve features near the surface with high resolution. The reflectivity maps representing slices of the 3D scattering volume are obtained through the reflection matrix procedure developed in [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. Body wave reflections from ambient noise correlations are used to image the fault zone structure up to 4 km below the surface. These images reveal the backscattered intensity generated by the distribution of heterogeneities in the medium.

Fault zones are very complex regions with extensive fracturing and damage that can reach the bottom of the seismogenic zone in some places as seen in tomographic studies around large faults (e.g. [START_REF] Allam | Seismic velocity structure in the hot springs and trifurcation areas of the San Jacinto fault zone, California, from doubledifference tomography[END_REF]). Tomographic and other imaging studies provide average properties of rock volumes, but do not resolve the presence, location and intensity of scatterers that are imaged with the reflection matrix method. The strong variations of velocities and significant attenuation in fault zone regions present challenges for conventional imaging techniques. However, our reflection matrix approach allows us to derive the distribution of scatterers inside the medium with an approximate velocity model of the medium. To that aim, a powerful aberration correction process is performed and provides high resolution images of the subsurface. However, associating a reflector with a specific depth remains dependent on the reference velocity model, so the reflectivity maps are displayed as a function of an effective depth and the observed two-way travel time.

A significant advantage of the matrix approach is that focusing inside the medium enables the imaging of not only specular reflectors but also of scattering objects such as cracks and fractures. While many methods consider the diffracting and scattering components as noise in the seismic data, and tend to remove these components to image discontinuous layers, the current approach takes advantage of the scattering in the complex fault zone to resolve features of the order of the wavelength. This constitutes one of the main strength of the method. The images in Figs. II.10 and II.11 show both discontinuities of some layers that are signature of a large fault, along with lateral and axial variations in the backscattered intensity induced by cracks and other small-scale heterogeneities. The axial variations of the backscattered intensity in Fig. II.12 also reveal systematic differences in scattering properties in the region within and around the surface fault traces relative to the outside volume, as well as differences between the NW and SE portions of the study area. The results are consistent with more localized intense damage zone at depth in the SE section where the SJFZ enters the Horse canyon, and more diffuse rock damage to the NW where the SJFZ is in a linear valley. The higher damage at depth in the SE section also suggests a geometrical complexity at depth leading to an ongoing generation of rock damage that is overprinted on older healed damage [START_REF] Sharp | San jacinto fault zone in the peninsular ranges of southern california[END_REF]. The response matrix is projected onto a focused basis at depth z both at emission (in) and reception (out). The focusing points synthesize virtual sources and receivers scanning every point in the virtual plane. In the presence of heterogeneities, the waves are distorted and the focusing is not optimal. (b) Confocal imaging principle: The image of the medium is obtained when input and output focusing are performed at the same location in the focused basis. (c) Far-field projection of the focused reflection matrix: For each virtual source, the reflected wave-front is investigated in the far-field. The phase distortions are identical for nearby virtual sources belonging to the same isoplanatic patch. (d) The phase conjugate of the aberration phase law enables a fine compensation of aberrations over the corresponding isoplanatic patch.

sponds to the location of virtual geophones inside the medium. A constant velocity model of 1500 m/s is chosen for the transmission matrix G 0 (f ) that describes the propagation of waves between the surface and the focused basis. At each frequency, the response matrix is projected into the focused basis both in emission and reception (Fig. Aa). The broadband focused reflection matrix R rr (z) is then computed by calculating the coherent sum of the focused reflection matrices over the considered frequency range. This last operation amounts to a ballistic time gating of singly-scattered waves at time t ∼ 2z/c 0 . The coefficients R (r out , r in ) of this matrix correspond to the responses between a set of virtual sources r in = (x in , y in , z) and receivers r out = (x out , y out , z) at each depth z. Among all these coefficients, the diagonal elements (r in = r out ), are of particular interest since they provide a confocal image of the underground (Fig. Ab). The other coefficients of R rr are also useful since the focusing quality can be directly quantified by the spreading of the backscattered energy over its off-diagonal elements. In the presence of strong phase distortions (aberrations), resulting from a very complex seismic velocity distribution, each input and output focal spot can actually spread well beyond the diffraction limit, giving rise to a loss of resolution and contrast of the confocal image. Interestingly, the reflection matrix can also be used to retrieve in post-processing the shape of wave-fronts that would allow a perfect compensation for these phase distortions. To do so, one can exploit the spatial correlations that exist between the phase distortions urdergone by the reflected wave-fronts induced by nearby virtual geophones. This coherence area is referred to the isoplanatic patch (Fig. Ac). The idea is to project the focused reflection matrix at input or output in a basis that maximizes the size of such isoplanatic patches in the focal plane. For a multi-layered basis, the most adequate basis is the Fourier basis (k) that amounts to projecting the reflected wave-field in the far-field. The aforementioned spatial correlations are then leveraged by extracting the distorted component of each reflected wave-field. It results in a distortion matrix whose singular value decomposition provides a decomposition of the field-of-view into a set of isoplanatic patches with the associated aberration phase transmittances. The phase conjugate of each aberration phase law provides the focusing law that should be used to compensate for aberrations on each isoplanatic patch (Fig. Ad). Such a focusing law is applied here to the reflection matrix from the dual basis to obtain an updated focused reflection matrix. A corrected confocal image is obtained by considering the diagonal coefficients of this new focused reflection matrix.

The whole process can be iterated at input and output in order to refine our estimation of the focusing law over each isoplanatic patch. To that aim, the focused reflection matrix should be projected in the far-field (Fourier basis k) by switching, at each iteration, between input or output.

B Specular vs. diffuse reflection

The finite size of the geophone array has an impact on the field-of-view that can be imaged by the reflection matrix approach. To quantify this effect, two synthetic tests have been performed: (i) considering planar reflectors at a given depth z = 3600m (specular scattering, see Fig. Ba1); (ii) considering randomly distributed scatterers at the same depth (diffuse scattering, see Fig. Bb1). In each case, an homogeneous effective wave velocity c 0 = 1500m/s is considered. Fig. Ba2 and Bb2 show the results of the numerical simulations.

Fig. Ba2 shows a synthetic confocal image obtained for the specular reflector. The field-of-view is limited by the size of the array. The fraction of the incident wave-front that spreads beyond the transverse size of the array gives rise to reflected waves that cannot be recorded by the array of geophones (red arrows in Fig. Ba1). Fig. Bb2 displays the field-of-view of our matrix imaging method in the diffuse scattering regime. This intensity distribution is obtained by averaging the synthetic confocal image over a number N = 50 realizations of disorder. Unlike for specular reflectors, incident wave-fronts are reflected in different directions (arrows in Fig. Bb1). Consequently, the field-of-view spreads well beyond the transverse size of the geophone's array. The latter observation explains the wide field-of-view obtained in the 3D-scattering images of the accompanying paper (Fig. II.11). The reflection matrix method mainly images a distribution of small heterogeneities that are located in the damage area at shallow depths (z < 1500m) and that reside between each geological layer at larger depths.

C Time-gain compensation

The energy losses increase with the heterogeneity of the medium. The propagation matrix used to project the raw data to the focused basis accounts for geometrical spreading whereas the effect of absorption and scattering are ignored. In particular, in strong scattering regime, these losses can strongly degrade the contrast of the images at larger depth.

To overcome these problems in the imaging process and in order to visualize the entire field-of-view, it is mandatory to compensate for the amplitude drop in the 3D-images (Fig. II.10 and II.11), especially in the shallow crust. This time gain compensation is done manually by multiplying the intensity profiles by an increasing function with depth, that is in this case the reciprocal of the mean intensity calculated at each depth.

An estimation of the seismic wave attenuation at the SGB site can be directly measured through the backscattered intensity at each depth. The energy is obtained by calculating the mean square of the intensity of the image pixels, i.e. the intensity of the diagonal of the reflection matrix. The amplitude decay of the energy is expected to follow an exponential decaying as a function of depth (Fig. C).

Figs. Ca and b displays the natural logarithm applied to the mean intensity as a func-tion of depth calculated from the raw reflection matrix and the corrected reflection matrix respectively (aberration correction process disclosed in [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]). In the logarithmic scale, the mean confocal energy is expressed as a decreasing linear line. The characteristic length of attenuation, that corresponds to the inverse of the slope of the log-energy decay, is also shown.

A change of slope is observed in both plots corresponding to depth of 750 m (t = 1 s). While a pronounced decrease in intensity occurs until 750 m, a more gentle slope is noticed below 750 m. The presence of stratigraphic boundaries, absorbing sediments and intense fracturing is a possible cause of the high attenuation observed at shallow layers.

Another thing to notice is the increase by around twice the characteristic length values after correction. We recall that the energy decay describes the losses that the seismic waves undergo while propagating inside the medium. In case of strong inhomogeneities and velocity variation, the wave-field is heavily distorted and the focusing inside the medium This chapter presents a large scale application of the reflection matrix approach that consists in imaging the deep structure of the NAFZ. Compared to the previous chapter that dealt with a sparse scattering medium, the subsurface here exhibits a continuous reflectivity made of specular reflectors and randomly distributed heterogeneities. Matrix methods are refined to optimize the focusing process in such a scattering regime. A 3D diffraction-limited image of NAFZ is finally obtained and a structural interpretation is provided.

III.1 Abstract

Imaging the structure of major fault zones is essential for our understanding of crustal deformations and their implications on seismic hazards. Investigating such complex regions presents several limitations, including the variation of seismic velocity due to the diversity of geological units and the cumulative damage caused by earthquakes. Conventional imaging techniques are in general sensitive to the velocity model. Here we apply a passive matrix approach to image the North Anatolian Fault using the Dense Array for Northern Anatolia (DANA) consisting of 73 stations with a nominal inter-stations spacing of ∼ 7 km that was deployed for 18 months in the region of the 1999 Izmit earthquake. In this region, the fault splits into two strands separating three major crustal blocks. We used the body waves from ambient noise cross-correlations in the [0.1 0.5] Hz frequency band to build the reflection matrix gathering all the information on the medium. A set of matrix operations are performed to project this matrix in a focused basis where the confocal image of the medium is built. Due to lateral velocity variations, the image suffers from phase distortions and is not interpretive. To overcome this issue, a compensation of these aberrations is performed through the distortion matrix concept. This matrix, that connects each point of the medium to the distortion of the wave-front, is investigated locally to estimate the transmission matrix between the surface and any point inside the medium.

The phase conjugate of this transmission matrix yields a set of optimized focusing laws to retrieve a 3D diffraction-limited image of the subsurface. The final image reveals the scattering structure of the crust and upper mantle around the NAFZ reaching 60 km depth with a transverse resolution of 2 km. The results show that most of the scattering is associated with the northern branch, that cuts throughout the crust and penetrates into the upper mantle. Our image displays first evidence that the Armutlu-Almacik (AA) block is narrower in depth than it appears at the surface in the East of DANA array. Sub-Moho structures are also identified below AA and Istanbul zone (IZ) in the North. The passive imaging method provides new insights into the structure of the NAFZ. Our results show strong evidence that the northern branch of the NAFZ cuts the crust and extends to the upper mantle. The structure of the major geological units is also revealed.

III.2 Introduction

The North Anatolian Fault zone (NAFZ) is one of the major continental right-lateral strike slip faults, and forms a border between the Eurasian continent and the Anatolian block.

With an extremely well developed surface expression, it is one of the most active faults in the Eastern Mediterranean region [START_REF] Sengör | The north anatolian transform fault: its age, offset and tectonic significance[END_REF][START_REF] Barka | The north anatolian fault zone[END_REF]. It is over 1600 km long and extends from eastern Turkey in the east to Greece in the west and historically has III.2. INTRODUCTION been subjected to many destructive earthquakes [START_REF] Ambraseys | The seismicity of Turkey and adjacent areas: a historical review[END_REF][START_REF] Stein | Progressive failure on the north anatolian fault since 1939 by earthquake stress triggering[END_REF]. The seismic activity of such large faults constitutes a continuous hazard/threat to the surrounding regions and big cities, especially Istanbul city located to the West of the fault.

Faults are well defined at the surface by the localized deformation and displacement delineating the fault traces, but the deep structure and the heterogeneities' distribution remain poorly understood [START_REF] Vauchez | Faults (shear zones) in the earth's mantle[END_REF]. The understanding of such major fault systems and seismic hazard requires a characterization of the geometrical and seismic properties of the crust and upper mantle. A large number of geological and geophysical studies discussed the complexity of fault zones. They are not only confined in the mid crust, they penetrate deep into the crust and may extend to the upper Mantle [START_REF] Handy | Continental fault structure and rheology from the frictional-to-viscous transition downward[END_REF]. If so, faults develop into shear zones, representing a zone of localized deformation in which the displacement occurs.

Seismic imaging techniques, especially seismic reflection, refraction, and tomographic methods, constitute a very powerful tool to characterize fault zones and report the variation of the properties of the crust and the upper mantle. They rely on the study of wave propagation inside the Earth that is governed by the density and elastic properties of the rocks. To properly probe the medium, waves should be generated by a dense distribution of seismic sources. Conventional seismic exploration techniques use either the earthquakes as seismic source, or explosions and vibrators to generate seismic waves in regions with weak seismicity. Because of the limitations and high cost of these methods, there was a need for alternative imaging approaches that would not rely on any coherent source. In the 2000's, the extraction of deterministic information about the Earth structure from ambient seismic noise revolutionized the field of seismology (e.g. [START_REF] Campillo | Seismic imaging and monitoring with ambient noise correlations[END_REF].

It was shown that the cross-correlation of ambient noise recorded at two stations provides the Green's function between these two stations [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF]. The reflection response of the medium is then retrieved and can be applied to build tomographic or structural images of the Earth. Because ambient noise is predominated by surface waves, their

Green's function component can be easily extracted [START_REF] Shapiro | Emergence of broadband rayleigh waves from correlations of the ambient seismic noise[END_REF]. It has been proved that body wave reflections can also be retrieved from ambient noise cross-correlations [START_REF] Wapenaar | Retrieving the elastodynamic green's function of an arbitrary inhomogeneous medium by cross correlation[END_REF][START_REF] Draganov | Retrieval of reflections from seismic background-noise measurements[END_REF]Poli et al., 2012a,b). Body waves contain valuable information on the structure of the medium and can be investigated to obtain high-resolution images of the crust and the mantle [START_REF] Retailleau | Ambient seismic noise imaging of the lowermost mantle beneath the North Atlantic Ocean[END_REF].

Faults are imaged indirectly through strong velocity contrasts in tomographic profiles [START_REF] Mordret | Shallow three-dimensional structure of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic array[END_REF], or through the offset of geological layers observed in reflectivity images [START_REF] Qian | Imaging the active faults with ambient noise passive seismics and its application to characterize the huangzhuang-gaoliying fault in beijing area, northern china[END_REF]. Here we aim to characterize the crustal structure of the NAFZ by imaging the major reflective features and the distribution of heterogeneities in the fault region. We use the body waves from ambient noise records to obtain 3D images of the fault by applying a reflection matrix approach originally developed in ultrasound imaging [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF]Lambert et al., 2020a) and optics [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF].

It was shown to provide high-resolution images of the Earth in complex areas, such as volcanoes [START_REF] Giraudat | Passive seismic matrix imaging of la soufrière of guadeloupe volcano[END_REF] and fault zones [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. While most conventional reflection techniques are strongly sensitive to the available velocity model, the reflection matrix approach is robust with respect to the inaccuracy of this model. An approximate velocity distribution is actually sufficient since a time-reversal analysis of seismic data enables a local and adapted auto-focus on each part of the subsurface image.

Our approach is based on the passive measurement of the reflection matrix K associated with a network of geophones. It contains the set of impulse responses between each geophone extracted from cross-correlations of seismic noise. Based on the available velocity model, a focused reflection (FR) matrix is built by redatuming applied to K [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]Lambert et al., 2020a). It contains the impulse response between virtual sources and receivers synthesized inside the medium. This FR matrix is powerful as it provides: (i) a confocal image of the subsurface along its diagonal; (ii) a local quantification of aberrations in the vicinity of each virtual source from its off-diagonal elements. These aberrations originate from the mismatch between the velocity model and the actual velocity distribution. The fluctuations of wave velocities actually induce phase distortions on the focused wave-fronts that may result in blurry images and artifacts. To overcome these detrimental effects, a novel operator, the distortion matrix, has been introduced. It contains the phase distortions of the reflected wave-fronts with respect to the propagation model. A time-reversal analysis of this matrix allows to retrieve the wave-

III.2. INTRODUCTION

fronts that can be used to perfectly focus on any virtual detector in the field-of-view of the geophone's network. The distortion matrix concept has been applied to restore the focusing in all kind of scattering regimes: in ultrasound, for randomly distributed scatterers (Lambert et al., 2020b(Lambert et al., , 2021b)); in optics for specular reflectors [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF];

and recently in seismology to image the San Jacinto Fault Zone scattering structure that exhibited a sparse distribution of scatterers [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF][START_REF] Touma | Distribution of seismic scatterers in the san jacinto fault zone, southeast of anza, california, based on passive matrix imaging[END_REF]. In this latter work, despite the strong variation of velocity in the complex fault zone, 3D images of the first 4 km of the crust were retrieved. Specular reflectors and scattering objects such as cracks and fractures were imaged with an unexpected resolution. This paper aims to demonstrate the applicability of the matrix approach at crustal scale in the NAFZ region. In contrast with [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF], where a homogeneous velocity model is used, in this work we integrate a layered velocity model to reduce the aberration level and obtain a better timeto-depth conversion. Moreover, as we will see, the NAFZ subsurface does not exhibit a distribution of sparse scatterers but rather a continuous reflectivity distribution made of specular reflectors and randomly distributed heterogeneities. We will take advantage of this continuity to adopt a local analysis of the distortion matrix in order to retrieve to the transmission matrix between the surface and any point inside the medium. This transmission matrix is of great interest since its phase conjugate provides the optimized focusing laws that need to be applied to the reflection matrix in order to retrieve a diffractionlimited image of the subsurface. As we will show, the obtained image will extend to 60 km depth below the surface with a transverse resolution of 2 km.

To image the crustal structure of NAFZ , we use data from the Dense Array for Northern Anatolia (DANA, 2012) that was deployed over the western segment of the fault, in the latest rupture region during the 1999 Izmit (M = 7.6) and Düzce (M = 7.2) earthquakes [START_REF] Barka | The surface rupture and slip distribution of the 17 august 1999 izmit earthquake (m 7.4), north anatolian fault[END_REF][START_REF] Akyuz | Surface rupture and slip distribution of the 12 november 1999 duzce earthquake (m 7.1), north anatolian fault, bolu, turkey[END_REF]. where most of the slip is accommodated [START_REF] Barka | The north anatolian fault zone[END_REF][START_REF] Reilinger | Gps constraints on continental deformation in the africa-arabia-eurasia continental collision zone and implications for the dynamics of plate interactions[END_REF], has been subjected to a series of major earthquakes in the last century. It developed more seimic activity than the southern branch and lately ruptured during Izmit Earthquake. On the contrary, the latest rupture of the southern branch dates back to the fifteenth century (Am- in the center, between the two fault strands. The Istanbul zone consists of transgressive sediments overlaying a granitic and metamorphic crystalline basement of a late-Cambrian age [START_REF] Chen | Geochronology of the karadere basement (nw turkey) and implications for the geological evolution of the istanbul zone[END_REF]. The Armutlu-Almacik block is formed of mixed pre-Jurassic units and metamorphosed sediments and marbles of unknown age and provenance [START_REF] Yılmaz | Geological evolution of the late mesozoic continental margin of northwestern anatolia[END_REF][START_REF] Okay | Tethyan sutures of northern turkey[END_REF]. The Sakarya Terrane in the South consists of an accretion complex of Jurassic sedimentary and volcanic rocks overlying a metamorphic basement [START_REF] Yılmaz | Geological evolution of the late mesozoic continental margin of northwestern anatolia[END_REF].

This study reveals the 3D scattering structure of the medium below this major fault with a transverse and axial resolution of 2 km and 800 m, respectively. It does not only image planar interfaces, but provides a direct insight on the heterogeneities that mainly sit in the vicinity of the strands. The observed results complement previous studies conducted in the region. A step in the Moho is detected below the Northern branch, and several sub-Moho structures are observed in the North confirming that the northern branch penetrates in the upper mantle. The southern strand does not have a strong signature in the scattering profiles. The evolution of the AA block in the crust is described. The block seems to be thinner at depth in the East of DANA array than it appears at the surface.

III.3 The reflection matrix approach III.3.1 The reflection matrix in the geophones basis

To apply our imaging technique, we used the ambient seismic noise recorded at DANA (see Fig. dows. Each window is band-pass filtered between 0.1 and 0.5 Hz after applying a spectral whitening between 0.01 and 1 Hz [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF]. The cross-correlation between every pair of stations is computed over the 1h windows and finally stacked to obtain the mean cross-correlation function. Although, considering noise in higher frequency range would allow, in principle, to improve the resolution of the images, matrix imaging requires the Nyquist criterion to be fulfilled: the inter-station distance (7 km) shall be of the order of a half-wavelength. Considering a bulk wave velocity c 0 = 1700 m/s near the surface, this criterion led us to choose the 0.1 -0.5 Hz frequency range (λ = 5.7 km at the central frequency). The ambient noise energy in the frequency band considered in this study comes from the secondary microseisms (5-10 s period band) originated by the ocean [START_REF] Stehly | A study of the seismic noise from its long-range correlation properties[END_REF] and constitutes one of the most energetic parts of the seismic noise.

The symmetric cross-correlations can be stacked in a time-dependant response matrix K(t), where 'in' and 'out' denote downgoing and upgoing directions, respectively. One element k ij (t) of this matrix corresponds to the impulse response between geophones i and j. In other words, k ij (t) contains the reflection response recorded at receiver i if there were a virtual source at j. This reflection response includes not only direct arrivals, but also all the seismic wave-field scattered by the medium. The latter component can be used to retrieve local information on the medium's reflectivity. The scattered wave-field itself can be decomposed into two components: (i) a single-scattering contribution that is useful for imaging the scattering structure of the medium, and (ii) a multiple scattering contribution that is very detrimental for imaging purposes and results in spurious arrivals and ghosts in the reflectivity images. [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF] showed how multiple scattering can be overcome through a reflection matrix approach. It allowed them to reveal the internal structure of a highly heterogeneous volcano (Erebus, Antarctica) in which multiple scattering was originally largely predominant. As we will see in the next section, in the present case, we are not in the same situation. The single scattering component is dominant compared to the multiple scattering contribution.

An image of the medium reflectivity can be obtained by applying a double focusing operation to K(t) [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF][START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. It consists in backpropagating the response measured at the surface into wave-fields below the surface as if there were sources and receivers inside the medium. This is similar to the "wave-field extrapolation" concept that forms the basis of the migration process [START_REF] Berkhout | Wave field extrapolation techniques in seismic migration, a tutorial[END_REF]. It requires performing beamforming operations both at emission and reception. On the one hand, focusing in emission consists in applying appropriate time delays to the emitted sources so that waves constructively interfere and focus on one point inside the medium.

Physically, this operation amounts to synthesize a virtual source inside the medium. On the other hand, focusing in reception is carried out by applying proper time delays to the received signals so that can constructively interfere. As in emission, this focusing operation can be seen as the synthesis of a virtual receiver inside the medium. This operation is known as "redatuming" in seismology [START_REF] Berkhout | A unified approach to acoustical reflection imaging. II: The inverse problem[END_REF] and consists of virtually moving sources and receivers from the surface to the medium below (Fig. III.2a). Generally, an image of the sub-surface is built by considering the response of virtual source and receiver placed at the same location (Fig. III.2b). The principle of matrix imaging is to decouple both locations (Lambert et al., 2020a).

Mathematically, the response between virtual sources and receivers is obtained from the response matrix at the surface through the Green's function, that describes the propagation between each geophone and each point inside the medium. A model for the wave velocity distribution inside the medium is required. In this study, and since the horizontal EE cross-correlation functions are considered, only an estimation of the S-wave velocity is required. Unlike [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF] that considered a homogeneous P-wave velocity model, here a layered 1-D S-wave velocity model is used for the focusing process (see Table III.1). To that aim, a combination of two models derived by [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF], for the first 5 km, and by [START_REF] Karahan | Crustal structure at the western end of the north anatolian fault zone from deep seismic sounding[END_REF], for deeper layers. This layered model will allow a better time-depth conversion and will limit the aberration level in the subsurface image.

However, because of the complexity of the fault zone in the studied region and the presence of three major tectonic blocks, a lateral variation in the crustal structure is expected. Therefore, the considered velocity model does not account for the lateral changes of seismic velocity. The mismatch between our model and the actual velocity distribution results in a biased focusing process inside the medium and produce phase distortions (aberrations). The reflectivity image is then deteriorated and can lead to an incorrect interpretation of the crustal structure. We show in the following that the reflection matrix approach allows to overcome such limitations and optimize the focusing process for each pixel of the image.

Layer Depth (km) V S (m s -1 )

1 0 -1 1700 2 1 -3 5 2500 3
3 5 -14 3200 4 14 -26 3500 5

26 -40 3600 6 40 -60 4300

Table III.1: 1-D S-wave velocity model. Depth and V S following [START_REF] Karahan | Crustal structure at the western end of the north anatolian fault zone from deep seismic sounding[END_REF] and [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF].

III.3.2 Propagator from the geophones to the focused basis

The focusing operations described in the previous section allow to obtain a reflection matrix containing the impulse responses between a set of virtual sources and receivers mapping the medium. We now show how such reflection matrix can be obtained through simple matrix operations. This matrix can be expressed in three different bases: (i) the geophones basis where the matrix represents the cross-correlations between all pairs of stations located at s(x s , y s , 0), (ii) the focused basis corresponding to the location r(x, y, z)

of virtual sources and receivers synthesized by the focusing operations and in which the image of the medium is built, and (iii) the spatial Fourier basis k(k x , k y ) that will be first used for wave-field extrapolation and then for aberration correction.

Switching between bases can be easily achieved by simple matrix products in the frequency domain. Consequently, a temporal Fourier transform is first applied to the measured response matrix K(t). For each frequency in the bandwidth of interest, a monochromatic matrix K(f ) is obtained. The forward and backward extrapolation of the response matrix is performed through decomposing the wave-field into plane waves [START_REF] Berkhout | Wave field extrapolation techniques in seismic migration, a tutorial[END_REF]. Indeed, plane waves can be easily extrapolated by applying a simple phase shift.

The decomposition of the wave-field recorded at z = 0 into plane waves can be done by considering the spatial Fourier transform operator T 0 = [T 0 (s, k)] that connects the geophone's position s to the transverse component k = (k x , k y ) of the wave vector:

T 0 (s, k) = exp (i(k x • s x + k y • s y )) , (III.1)
The symbol • denotes a scalar product. We then define the spatial transfer function exp (-ik

(i)
z ∆z i ), that models the ballistic propagation of the seismic wave through the i th layer with k

(i) z = 2πf c i 2 -k 2
x -k 2 y , the longitudinal component of the wave vector k in the i th layer, c i the wave velocity in the i th layer of our model and ∆z i its thickness.

To propagate the plane waves from the surface to depth z, we define the wave-field extrapolator, F(z, f ) = [F (k, z, f )], as the product of the spatial transfer function of the N layers above the considered depth as follows:

F (k, z, f ) = exp(-ik (N ) z (z -z N )) * N -1 i=1 exp (-ik (i) z ∆z i ), (III.2)
where z i is the depth at which starts the i th layer. The far field propagator, P(z, f ) = [P (s, r, z, f )], that models wave propagation from the Earth surface to depth z can finally be expressed as follows:

P(z, f ) = [T 0 • F(z, f )] × T † 0 (z). (III.3)
where the symbols † and • stands for transpose conjugate and Hadamard product, respectively. T 0 (z) = [T 0 (r, k, z)] is the free-space transmission matrix that basically connects the plane wave basis with the focused basis at depth z. T 0 (z) corresponds to the Fourier transform operator. It connects the transverse wave vector k = (k x , k y ) of each plane wave to any focusing point r defined by its transverse coordinates (x, y) at depth z:

T 0 (r, k) = exp (ik • r). (III.4)
To avoid aliasing during the change of basis between the plane wave and focused bases, a Shannon criterion should be respected. The spatial frequencies (k x , k y ) are supposed to fulfill the following conditions:

k x < 1 2δ 0
, and k y < 1 2δ 0 (III.5) with δ 0 the transverse resolution in the focal plane that corresponds to the distance between the focusing points r. δ 0 is chosen ∼ 2 km ∼ λ min /2 to circumvent spatial aliasing in the focused basis. Theoretically, δ 0 is supposed to be ∼ λ min /4 [START_REF] Balondrade | Approche matricielle de la microscopie optique: quantification et correction des aberrations[END_REF], but due to the small aperture of the array, and to reduce the computational cost, λ min /2 was found to be sufficient to avoid aliasing.

III.3.3 The focused reflection matrix

As we will see now, the propagator P defined by Eq. III.3 enables a direct projection of the response matrix from the geophones' basis to the focused basis based on the wave velocity model given in Table III.1. Each monochromatic response matrix K(f ) can be projected in the focused basis both at input and output by applying appropriate phase shifts (Eq. III.2) associated with downgoing waves at input and upgoing waves at output is a four-dimension matrix concatenated in 2D as a set of blocks [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF].

As in the SJFZ study [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF], it appears that the backscattered energy is concentrated in the vicnity of the diagonal of the focused reflection matrix. This indicates that single scattering dominates at this depth. However, the spreading of the backscattered energy outside the diagonal indicates the presence of phase distortions (aberrations) due to lateral variations of shear velocity that was not accounted for in the focusing operations.

Nevertheless, one can build an image of the medium reflectivity at effective depth z by considering the diagonal elements of the FR matrix that correspond to identical focusing points at input and output (Fig. III.2b). This is the principal of confocal imaging: The matrix Γ describes the scattering process in the focused basis. This matrix is diagonal in the single scattering regime, and its coefficients then correspond to the true reflectivity γ(r, z) at depth z. T(z) = [T (k, r, z)] is the true transmission matrix that describes wave propagation between the focused and the Fourier bases.

I (z) = |R (r,
As shown by Lambert et al. (2020a), a theoretical expression of R kk can be derived in the single scattering regime. Interestingly, the norm-square of its coefficients R(k out , k in , z) is shown to be independent of aberrations. It directly yields the spatial frequency spectrum of the scattering medium at depth z: focused basis for a local quantification of the focusing quality. The FR matrix actually allows to quantify the level of aberrations. To understand more quantitatively the manifestation of aberrations in R rr , it can be decomposed as follows (Lambert et al., 2020a(Lambert et al., , 2021a;;[START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]:

|R (k out , k in , z)| 2 = |γ (k out + k in , z)| 2 . ( III
R rr (z) = H (z) × Γ(z) × H(z), (III.12)
or, in terms of matrix coefficients, R(r out , r in ) = ˆdrH(r, r out )γ(r)H(r, r in ).

(III.13) H = T 0 × T is the focusing matrix that contains the transmit/receive point spread function H(r, r in/out ) around each focusing point r in/out in the field-of-view. The amplitude distribution of H(r, r in/out ) accounts for the lateral extent of each vitual source/detector at r in/out . Eq III.13 confirms that the diagonal elements of R rr form an estimator of the subsurface reflectivity: it results from the convolution between the true reflectivity and H 2 , the product between the input and output PSFs.

To assess the aberration level, i.e. the quality of focusing, we now investigate the off-diagonal elements of the FR matrix that provide an estimation of the imaging PSF for any point of the image. Each coefficient along an anti-diagonal of R rr share the same midpoint r m = (r out + r in )/2 and its norm-square yields the common mid-point (CMP) intensity profile (Lambert et al., 2020a;[START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF].

I(r m , ∆r) = |R(r m -∆r/2, r m + ∆r/2)| 2 . (III.14)
In a sparse scattering regime [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF], the CMP intensity profile has to be considered in the vicinity of the brightest structures in the field-of-view. Here, on the contrary, the reflectivity of the medium displays a speckle-like feature. A statistical approach can be adopted by performing a spatial average of I(r m , ∆r) over the whole fieldof-view. In the speckle regime, the averaged CMP is proportional to the convolution between the incoherent output and input PSF, independently from the medium's reflectiv-ity : I(r m , ∆r) ∝ |H| 2 * |H| 2 (Lambert et al., 2020a). Consequently, by looking at the mean intensity of the anti-diagonal elements, we can probe the spatial extension of the PSF and estimate the mean transverse resolution of the image. In the following we will show that it is possible to restore an optimal resolution and a high contrast for the image through the distortion matrix concept.

III.3.6 The distortion matrix

The distortion matrix was already introduced in ultrasound (Lambert et al., 2020b(Lambert et al., , 2021b)), optics [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF] and seismology [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. Several applications proved the efficiency of this matrix in overcoming aberrations and improve the image quality. The distortion matrix concept is here reviewed and will be applied to retrieve the structure of the NAFZ.

In order to improve the focusing, the aberrations should be investigated and corrected.

To isolate the aberrations, the first step is to chose a basis in which the distortion of the wave-front is the most spatially-invariant. In a horizontally multi-layered medium as NAFZ, the far field-basis is the most adequate since plane waves are the propagation invariants in this geometry. A far-field projection implies performing a spatial Fourier transform to the FR matrix. We first study this projection at the output of R rr : The reflected wave-field can be understood as a sum of two components: (i) a geometric component described by the reference matrix T 0 , containing the ideal wave-front generated from a source r in in a perfectly homogeneous medium (dashed black curves in In terms of matrix coefficients, the distortion matrix can be written as:

D(k out , r in , z) = R(k out , r in , z)T * 0 (k out , r in ).
(III.17)

The matrix D kr connects any input virtual source r in to the distorted component of the reflected wave-field measured in the far-field (Fig. III.4d). By removing the geometrical component of the reflected wave-field, spatial correlations are highlighted between distorted wave-fields induced by neighbour virtual sources r in [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]. Such correlations are a manifestation of a spatial invariance of aberrations over areas generally referred to as isoplanatic patches [START_REF] Roddier | Adaptive optics in astronomy[END_REF]. As shown in previous papers [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]Lambert et al., 2020bLambert et al., , 2021b)), the idea is to exploit these correlations in order to extract the aberration phase law associated with each isoplanatic patch in the field-of-view. This will be the focus of the next section.

III.4 Matrix approach for local aberration correction

Recent works in seismology [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF] and optics [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF] have

shown that for certain scattering regimes (specular reflectors or sparse scattering), there was a one-to-one association between the eigenstates of D and the isoplanatic patches present in the field-of-view. However, in a more general case of a random reflectivity (speckle regime), which is the case in the NAFZ, this property does not hold and local distortion matrices should be considered over smaller areas in which the isoplanatic hypothesis is ideally fulfilled (Lambert et al., 2021b).

In this section, a correction process suitable for the case of the NAFZ is introduced.

The process is outlined by five steps: (i) projection of the broadband FR matrix at input or output into the plane wave basis, (ii) the realignment of the transmitted or reflected wavefronts to form the distortion matrix D, (iii) the truncation of D into overlapping patches, (iv) the singular value decomposition of D or of its normalized correlation matrix to extract an aberration phase law for each patch, (v) the application of the phase conjugate of the wave-front distortion as a focusing law to correct for output and input aberrations.

All of these steps are repeated by exchanging input and output bases.

III.4.1 Isoplanicity and local distortion matrices

A local analysis of the distortion matrix is here presented. As we will see, this local analysis will provide an estimator T of the transmission matrix T between the far-field and any focusing point mapping the subsurface. where W ∆r (r) is a spatial window function such that W ∆r (r) = 1 for |x| < ∆x and |y| < ∆y, and zero elsewhere. Ideally, wave-front distortions should be invariant over each region, meaning that the virtual sources r in associated with each region belong to the same isoplanatic patch. However, in practice, this hypothesis is not fully verified.

The isoplanatic length actually scales as the typical transverse dimension over which the wave velocity fluctuates. On the one hand, the dimensions ∆r of the window function should therefore be reduced to cover the smallest isoplanatic region as possible in order to provide a local and sharp measurement of aberrations. On the other hand, it should also be large enough to include a sufficient number of realizations of disorder in order to unscramble the effect of aberrations from the medium's reflectivity (Lambert et al., 2020b;[START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF]. To reach a good estimation of the aberration phase law, the The isoplanatic condition is now assumed to be fulfilled over each region of spatial extension ∆r. Aberrations can thus be considered as invariant : H(r, r in/out ) = H(rr in/out , r p ) (Lambert et al., 2021a). Under this assumption, the coefficients of D kr (r p ) matrix can be expressed as follows (Lambert et al., 2021b):

D(k out , r in , r p ) = H(k out , r p ) ˆdrγ(r + r in , z)H(r, r p )e ikout.r , (III.19)
where H(k out , r p ) is the aberration transmittance and is the Fourier transform of the output PSF H(r, r p ): H(k out , r p ) = ´drH(r, r p )e -ikout.r . Eq. III.19 can be seen as a product between two terms: the output aberration transmittance and a virtual source term modulated by the medium's reflectivity γ(r). In the following, we will show how these two terms can be discriminated in order to get a proper estimation of the aberration trans-mittance H(k out ; r p ) at each point r p . This quantity is crucial since it will provide an estimator of the transmission matrix, T = T 0 • H, whose phase conjugate will directly provide the focusing laws to be applied for an optimized focus on each patch. To retrieve H(k out , r p ), we will take advantage of the correlations that exist between each distorted wave-fied in each patch (Fig. III.4d). In practice, this can be done through an iterative time reversal analysis of each matrix D kr (r p ).

III.4.2 Iterative time reversal analysis of the distortion matrix

Mathematically, it consists in a SVD of each local distortion matrix D kr (r p ):

D kr (r p ) = U out (r p ) × Σ(r p ) × V in (r p ) † (III.20)
where Σ is a diagonal matrix containing the real positive singular values σ i in a decreasing

order σ 1 > σ 2 > • • • > σ N .
U out and V in are unitary matrices whose columns, U

(i) out = [U (i) out (k out )] and V (i) in = [V (i)
in (r in )], correspond to the output and input singular vectors, respectively. Under the isoplanicity condition, D kr (r p ) is, in first approximation, of rank 1 and the normalized output singular vector Ûout (r p ) = [U [START_REF] Yao | Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to se tibet[END_REF] out (k out , r p )/|U (1) (k out , r p )|] yields an estimation of the aberration transmittance H(r p ) (Lambert et al., 2020b(Lambert et al., , 2021b)). shows a diffraction limited central lobe but also a diffuse component that spreads well beyond the diffraction limit. The latter contribution can be explained by the subsistence of aberrations at the input of the reflection matrix.

The input aberrations can be corrected through the same correction process described above, but after exchanging the focused and the Fourier basis between the input and output [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. To that aim, the dual reflection matrix

R rk (z) = [R(r out , k in , z)]
is built from the corrected FR matrix: R rk = R (c) rr × T t 0 . The input distortion matrix D rk (z) = [D(r out , k in , z)] (Eq. III.16) is then deduced. Instead of using as before a SVD of the distortion matrix to retrieve the aberration transmittance, it is more efficient to consider the normalized correlation matrix of the input distortion matrix [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF].

To that aim, the correlation matrix, C kk (r p ) = D rk (r p ) † × D rk (r p ), is computed and its Fourier bases between input and output [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. The eigenvectors of the normalized correlation matrix provide the phase laws that should be used to compensate for the residual high order aberrations, first at the output and then at the input. After The gains in contrast (dB) are reported in Table III The 3D image of the medium around the NAFZ is obtained by combining the images derived at each depth. Four iterations of the correction process at each depth are performed as described previously. The scattering generated by the heterogeneities of the medium yields to energy losses with depth. Consequently, drop of amplitude is observed in the 3D images. In order to compensate for the energy losses, the intensity in the cross-sections is normalized by the mean intensity calculated at each depth. In the next section, the 3D scattering structure of the NAFZ is revealed at several location below DANA network, and a structural interpretation of the crust is provided with comparison with previous studies.

III.5 3D structure of the NAFZ

Our study area, the North Anatolian Fault (NAF) Zone in the broad sense, corresponds to the contact between two continental blocks, Istanbul block in the North (Rhodope-Pontide fragment [START_REF] Yılmaz | Geological evolution of the late mesozoic continental margin of northwestern anatolia[END_REF]) and Sakarya block in the South, whose seismic studies indicate that they correspond to lithospheres of different thicknesses [START_REF] Kounoudis | Seismic tomographic imaging of the eastern mediterranean mantle: implications for terminalstage subduction, the uplift of anatolia, and the development of the north anatolian fault[END_REF]. This contact zone includes the Armutlu block which is a remnant of the intrapontide suture, possibly associated with an ancient ocean (S ¸engör and [START_REF] Yilmaz | Tethyan evolution of turkey: a plate tectonic approach[END_REF] and its margins [START_REF] Akbayram | Evidence for a minimum 52±1 km of total offset along the northern branch of the north anatolian fault in northwest turkey[END_REF]. The boundaries of these units appear to be well defined in the South with the southern strand of the NAF while the northern boundary of the Armutlu block is hidden under recent sedimentary cover. The role of NAF in the structuring is of course very important. Since about 4 million years ago, a strike slip motion between the continental blocks reaches at least 88 km [START_REF] Akbayram | Evidence for a minimum 52±1 km of total offset along the northern branch of the north anatolian fault in northwest turkey[END_REF].

Most of which (67 km) occurred along the broad shear zone at the boundary of the Istanbul block that includes the surface expressions of the recent 1999 earthquakes [START_REF] Rockwell | Lateral offsets on surveyed cultural features resulting from the 1999 izmit and duzce earthquakes, turkey[END_REF].

Nevertheless, it is difficult to assert that the surface traces of recent earthquakes actually correspond to a contact at depth, nor to the position of a lithospheric fault, or a mantle shear zone. We therefore propose to produce scattering images that allow us to characterize the materials in terms of their heterogeneity, and in particular the zones of very strong deformation around the fault at depth. Note that due to its original nature as a subduction complex, the Armutlu block would be a heterogeneous feature as evidenced at the surface by intensely deformed rocks of deep origin, such as the Sapanca complex and its serpentinites [START_REF] Akbayram | Early cretaceous closure of the intra-pontide ocean in western pontides (northwestern turkey)[END_REF].

An interpretation of the geological features can be provided from the obtained scat- The 3D scattering volume details the structure of the medium below DANA array and shows clear differences between the major crustal blocks. In the following we highlight these differences and we describe the deep structures of the blocks revealed by our images.

Based on the scattering distribution, we determine the signature of the fault at depth. We also assess the lateral variations in the lower crust below the fault strands that give direct insight on the crust-mantle boundary location (the Moho). We detect the Moho depth in several regions and also heterogeneities in the upper mantle.

III.5.1 Previous geophysical studies on the NAFZ

Because of its considerable seismic activity in the last 100 years, a large number of studies have been conducted on the NAFZ in order to image its structure and understand its mechanical properties. Tectonic complexity in the area and differences in crustal composition, structure and properties between the blocks had been reported. In particular, strong velocity contrasts were found across the fault strands by several tomographic [START_REF] Salah | Seismic velocity structure along the western segment of the north anatolian fault zone imaged by seismic tomography[END_REF][START_REF] Koulakov | Distribution of seismic velocities and attenuation in the crust beneath the north anatolian fault (turkey) from local earthquake tomography[END_REF][START_REF] Papaleo | Seismic tomography of the north anatolian fault: New insights into structural heterogeneity along a continental strike-slip fault[END_REF][START_REF] Papaleo | Constraints on north anatolian fault zone width in the crust and upper mantle from s wave teleseismic tomography[END_REF] and full waveform inversion studies [START_REF] Fichtner | The deep structure of the north anatolian fault zone[END_REF][START_REF] ¸ubuk-Sabuncu | 3-d crustal velocity structure of western turkey: Constraints from full-waveform tomography[END_REF]. Local earthquake tomography [START_REF] Koulakov | Distribution of seismic velocities and attenuation in the crust beneath the north anatolian fault (turkey) from local earthquake tomography[END_REF] as well as S-and P-wave teleseismic tomography show that low velocity zones and high-attenuation are found below the surface traces of the SNAFZ and NNAFZ [START_REF] Papaleo | Seismic tomography of the north anatolian fault: New insights into structural heterogeneity along a continental strike-slip fault[END_REF][START_REF] Papaleo | Constraints on north anatolian fault zone width in the crust and upper mantle from s wave teleseismic tomography[END_REF]. These zones are correlated with the presence of fluids, fractured zones and fault damage [START_REF] Allam | Seismic velocity structures in the southern California plateboundary environment from double-difference tomography[END_REF]. The crust of Istanbul and Armutlu Blocks is characterized by high velocities while SZ shows relatively low velocities [START_REF] Koulakov | Distribution of seismic velocities and attenuation in the crust beneath the north anatolian fault (turkey) from local earthquake tomography[END_REF][START_REF] Papaleo | Seismic tomography of the north anatolian fault: New insights into structural heterogeneity along a continental strike-slip fault[END_REF][START_REF] Papaleo | Constraints on north anatolian fault zone width in the crust and upper mantle from s wave teleseismic tomography[END_REF][START_REF] Taylor | Near-surface structure of the north anatolian fault zone from rayleigh and love wave tomography using ambient seismic noise[END_REF].

Magnetotelluric (MT) measurements [START_REF] Tank | Magnetotelluric imaging of the fault rupture area of the 1999 izmit (turkey) earthquake[END_REF] along two profiles across the NAFZ show a variation in crustal conductivity between the South and the North. SZ in the South is associated with a lower resistivity (500 Ωm) than the IZ (≥ 1000 Ωm). On the contrary, the AA block in the center exhibits a conductive zone (30 -50 Ωm).

Several studies determined the crustal thickness at the east of the Sea of Marmara.

P-wave receiver function (RF) regional studies [START_REF] Zor | The crustal structure of the east anatolian plateau (turkey) from receiver functions[END_REF][START_REF] Vanacore | Moho structure of the anatolian plate from receiver function analysis[END_REF] detected a difference in the average Moho depth between the west (29-32 km) and the east (34-35 km). The westward crustal thinning was also reported in RF studies using data from DANA network [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF] with a deepening of the Moho beneath the AA block. A larger Moho depth in the IZ (∼ 40 km) was observed in the DANA array region based on analysis of P-wave transfer functions [START_REF] Frederiksen | Crustal thickness variations and isostatic disequilibrium across the north anatolian fault, western turkey[END_REF],

RF study [START_REF] Jenkins | Crustal thickness variation across the sea of marmara region, nw turkey: A reflection of modern and ancient tectonic processes[END_REF] and teleseismic scattering tomography [START_REF] Rost | Structure of the northwestern north anatolian fault zone imaged via teleseismic scattering tomography[END_REF].

Studies find strong evidence on the extension of the NNAFZ from the crust to reach the upper mantle [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF][START_REF] Papaleo | Constraints on north anatolian fault zone width in the crust and upper mantle from s wave teleseismic tomography[END_REF][START_REF] Jenkins | Crustal thickness variation across the sea of marmara region, nw turkey: A reflection of modern and ancient tectonic processes[END_REF][START_REF] Rost | Structure of the northwestern north anatolian fault zone imaged via teleseismic scattering tomography[END_REF], and also in the east of our study region [START_REF] Fichtner | The deep structure of the north anatolian fault zone[END_REF]. This observation is explained by a localized shear zone that extends into the mantle. On the contrary, no similar observations have been reported for the SNAFZ. ), that also appears as negative amplitudes in the receiver function study [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF], in agreement with ambient noise surface waves analysis [START_REF] Taylor | Near-surface structure of the north anatolian fault zone from rayleigh and love wave tomography using ambient seismic noise[END_REF].

The AA block appears to be less heterogeneous in the east compared to the west. The It appears to be the same structure identified by [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF]. It seems to be continuous and emerges in the western profile also.

The IZ profile (Fig. III.7a) displays complex structures appearing at 25 and 30 km (indicated by 2) consistent with the scattering tomography of [START_REF] Rost | Structure of the northwestern north anatolian fault zone imaged via teleseismic scattering tomography[END_REF]. In contrast with the IZ profile, the crust of the AA block (Fig. III.7b) appears to be relatively less heterogeneous. A near-surface reflective layer is observed and is consistent with the low velocity layer identified at 10 km (marked as 1) in the AA block [START_REF] Rost | Structure of the northwestern north anatolian fault zone imaged via teleseismic scattering tomography[END_REF]). Another reflective feature is visible around 20 km (marked as 2), also detected by [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF].

The SZ profile in Fig. III.7c shows relatively less reflectivity in the upper crust compared to the IZ and AA profiles, indicating that the crust of SZ is less complex than the two other blocks. This observation is also supported by the receiver function profiles from [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF]. This may be associated with the small slip rate of the southern strand [START_REF] Meade | Estimates of seismic potential in the marmara sea region from block models of secular deformation constrained by global positioning system measurements[END_REF], and therefore the SZ block incorporates less damage than the other blocks that are around the northern strand.

III.5.4 Lower crust and Moho

In the West (Fig. III.6a), a scattering lower crust is identified around 32 km. Its limit corresponds with the presence of the Moho (indicated by the red dashed line). The Moho is clearly detected below the southern branch at 32 km in the SZ, and deepens to 35 km beneath the AA block. The Moho reflectivity is less clear in the center and east of the profile, coinciding with the location of the northern branch. This observation is consistent with the decrease in amplitude in the receiver function study [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array[END_REF] beneath the IZ. The weak Moho reflection below the northern branch was also reported by the auto-correlation study [START_REF] Taylor | Crustal imaging across the north anatolian fault zone from the autocorrelation of ambient seismic noise[END_REF] and the absence of Moho signature was noted in the S-wave teleseismic tomography [START_REF] Papaleo | Constraints on north anatolian fault zone width in the crust and upper mantle from s wave teleseismic tomography[END_REF]. These findings support our result and suggest that the NNAFZ extends to large depth and is rooted in the mantle.

In comparison with the western profile, the eastern profile ( We note the overall clear concentration of scattering strength at depth in the region surrounding the northern strand. This suggest that the fault is rooted in a region beneath or slightly north of the northern branch, that is also corresponding to the limit of the Istanbul continental unit (Rhodope-Pondite fragment).

III.5.6 Fault dipping and Shear zone

In the scattering images, the signature of the fault strands is highlighted by the termination of crustal and sub-crustal structures. In the western profile, the structures terminate near both the strands, with a weak reflectivity around the southern strand. No reflective structures appear just below the southern strand in the Eastern profile, where the middle crust of IZ and AA block is nearly homogeneous.

As previously discussed, most of the reflectivity is concentrated around the northern branch, illustrating the fact that the northern branch is more active than the southern branch. The scattering below the northern branch reaches at least 60 km, confirming that the NNAF extends through the entire crust and penetrates into the upper mantle.

III.6 Discussion

The scattering volume shows distinct differences between the West and the East and between the North and the South. The observed variations are related to the variety of the geological structures of different tectonic origins forming the fault system. The cumulative slip of the fault is estimated to exceed 80 km [START_REF] Armijo | Westward propagation of the north anatolian fault into the northern aegean: Timing and kinematics[END_REF][START_REF] Bohnhoff | Maximum earthquake magnitudes along different sections of the north anatolian fault zone[END_REF], which results in the accumulation of damage in the region (Ben-Zion and Sammis, III.6. DISCUSSION The variation in the dipping of the fault between the West and the East can be explained as a consequence of the two fault strands converging to completely connect around 30.7 • E.

III.7 Conclusion

The matrix approach presented in this paper provides unprecedented view of the NAFZ.

We used seismic noise data from a dense deployment over the rupture region of the 1999 Izmit earthquake. By applying focusing operations to the ambient noise crosscorrelations, both at emission and reception, a focused reflection matrix is built in a virtual focused basis inside the medium. Virtual sources and receivers are synthesized, enabling to scan the medium and allow to take full advantage of all the information contained in the collected data, beyond the simple image of the subsurface. The strength of this approach is that it does not require an accurate velocity model. An approximate model is sufficient. Even though the differences between the considered model and the reality result in strongly aberrated image of the subsurface, our approach enables to correct for distortions in order to improve the focusing. To that aim, a distortion matrix is built between the focused and far-field bases. It contains the distortion of the wave-field with respect to our propagation model for each focusing point in the medium. A time-reversal analysis of this distortion matrix then provides the transmission matrix between the Earth's surface and the voxels mapping the medium. The latter matrix is a key quantity for imaging in heterogeneous media since its phase conjugate provides the phase laws to optimally focus everywhere inside the medium.

Previous works on passive seismic imaging have shown that the distortion matrix concept can be applied to sparse media made of a few isolated scatterers [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF][START_REF] Touma | Distribution of seismic scatterers in the san jacinto fault zone, southeast of anza, california, based on passive matrix imaging[END_REF]. Other works in ultrasound and optics proved that it is also efficient in the case of specular objects [START_REF] Badon | Distortion matrix concept for deep imaging in optical coherence microscopy[END_REF] and random single scattering media (Lambert et al., 2020b). The goal of the present article is to show that the method can be successfully applied for aberration correction in such scattering regimes also in seismology. Inspired by previous works in ultrasound (Lambert et al., 2020b(Lambert et al., , 2021b)), the method applied here differs with our previous work [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF] in several aspects. First, a layered velocity model is integrated during the focusing operations, enabling a better time-depth conversion. Secondly, the medium of the NAFZ displays both specular reflectors and a random distribution of scatterers. While in a sparse scattering regime, each eigenvector of the distortion matrix provides a focusing law associated with each scatterer of the medium, such analysis is not suitable in the case of a random reflectivity. A local aberration correction should thus be applied for the subsurface of the NAFZ. The idea is to take advantage of the local isoplanicity of the aberration phase law in each region of the medium. This approach was shown to be an appropriate way to handle lateral velocity variation and to compensate locally for phase distortions. It is done by retrieving the output and input phase law corresponding to each region and use it to estimate the true transmission matrix of the medium. Surprisingly, the phase laws were found to be different between the output and input. This inconsistency requires further analysis but might be solved by gradually reducing the size of the targeted patches at each step of of the time reversal scheme (Lambert et al., 2021b). Nevertheless, albeit perfectible, passive matrix imaging succeed in improving the resolution by a factor of ∼ 6 compared with the original image. Only high-order aberrations subsist at the end of the process. They are probably associated with isopalanatic area smaller that the size of the spatial windows used to truncate the field-of-view.

The scattering structure of the crust and upper mantle of the NAFZ continental strike slip fault was revealed. The 60 km depth profiles, having a transverse resolution of 2 km, show terminations of crustal discontinuities mainly below the northern branch. The localized scattering around the NNAF is consistent with the fact that it is the most seismically active fault and that it ruptured during the last 7.6 Izmit earthquake. We observe anomalies in the lower crust below the northern branch, with a step of the Moho coinciding with the surface location of this branch in the East of DANA network. Moreover, the scattering extends to the upper mantle in the North. All these observations are consistent with previous studies and suggest that the shear zone of the NNAF is localized in the crust and extends to the upper mantle. The fault dips to the North in the West and appears to be sub-vertical in the East.

The profiles show variation in the scattering structure between the West and the East.

The lack of reflectivity below the southern strand, together with the homogeneity of the crust in the central zone located between the two fault strands at the surface, suggest that the AA block is narrower in depth to the East than it appears at the surface. This block nearly disappears around 30.37 • E. Such observation needs of course further investigations and should be confirmed with other imaging techniques.

Finally, we show that continuous recordings of the seismic wave-field paves the way to develop new imaging techniques that provide deeper insight into the lithospheric structure around fault zones. We also prove that the presented technique constitutes a powerful passive imaging technique that can be applied to any scale when little knowledge on the distribution of the seismic velocity is available. It can also be applied in active imaging.

The reflection matrix contains information beyond the simple image of the medium. This approach also paves the way towards mapping the velocity distribution inside the medium.

This will be the focus of a future study. the overall distortion matrix, defined in Eq. III.16, is analysed as a whole. In other words, the SVD of the global distortion matrix is performed at each depth. The region over which the correction process is done includes the whole field-of-view.

Of course, the global distortion matrix is not of rank 1, but here we only consider the first eigenmode in the correction process. In fact, higher-order eigenmodes were here found to be linked with first eingenmode and do not bring supplementary information [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. At a first iteration, the transmission matrix T is estimated from the normalized singular vector Ûout of the global distortion matrix D as defined by Eq. III.21. The corrected FR matrix is then deduced from Eq. III.22. The aberration correction process can now be iterated by exchanging the focused and the Fourier basis between input and output as described previously. Starting from the second iteration, the eigenvectors of the normalized correlation matrix of the global distortion matrix are used to compensate for residual aberrations. lobe of energy outside the diffraction limited resolution is not negligible. This result shows that the global aberration correction is not optimal, even after four iteration steps.

A global correction is now adopted at each depth. A 3D scattering volume is obtained.

In In the previous chapters, we have seen the efficiency of the reflection matrix approach in retrieving the scattering images around fault systems. To do so, the reflection matrix was studied at the ballistic time. In the geophone's basis, the ballistic time is the roundtrip delay time that a singly-scattered wave will take to go from the geophones to the focal plane where it is scattered once before coming back towards the geophone's array.

Of course, the expected ballistic time is valid if and only if the propagation model is correct. If not, a time shift exists between the actual ballistic time and its prediction. In this chapter, inspired by the PhD work of William Lambert [START_REF] Lambert | Matrix approach for ultrasound imaging and quantification[END_REF], we will investi-IV.1. THE TIME FOCUSED REFLECTION MATRIX gate the reflection matrix beyond the ballistic time to investigate the focusing process in the spatio-temporal domain and exploit it to measure locally the velocity integrated over the path from the geophone's basis to a given focusing point in the subsurface.

This integrated velocity can lead to the estimation of a local velocity by solving an inversion problem. The access to this time shift allows to quantify the axial aberrations which allows to convert the travel time into a depth at which the scattering event takes place. In this preliminary study, the case of the North Anatolian fault is considered.

IV.1 The time focused reflection matrix

Analysing the reflection matrix beyond the focusing time first consists in focusing at emission at a point r in to create a virtual source. The medium is then scanned by performing focusing at reception at every location to create virtual receivers r out but, contrary to the previous chapters, an additional time delay δt is added (or subtracted) between the input and output focusing times. Mathematically, this can be done via an inverse Fourier transform of the focused monochromatic reflection matrices R rr (f ) defined in Eq. III.6, such as:

R rr (δt) = ˆf2 f 1 df R rr (f ) exp (i2πf δt), (IV.1)

with P(f ) the far field propagator over a volume of focusing points in depth. δt is a time vector including negative and positive time delays to follow the converging and diverging wave around each focusing point r in . R rr (δt) is a 7D matrix, with x, y and z coordinates for focusing points at emission and reception, and the time dimension. It gathers much more information on the medium than the FR matrix considered at the ballistic time. To simplify the analysis in the following, the matrix is calculated only in a single (y, z) plane (x in = x out will be fixed). In particular, it is investigated in order to retrieve information on the medium's local velocity. IV.2 Propagation movies IV.2.1 Correct wave velocity model

In the following, the time focused reflection matrix is studied in the virtual source frame:

the virtual source location is fixed. A change of basis is applied so that the transmitted response is described around the focusing point at ∆r = r out (x in , y out , z out )-r in (x in , y in , z in ).

We first consider the case of a homogeneous wave velocity model c 0 = 3500 m/s. The de-scanned reflection matrix R ∆r,r (δt) is calculated for a virtual source located on an echogenic structure for a two way travel time t = 27.8 s. that due to reciprocity, the analysis could be performed in the virtual receiver's frame, where the time focused reflection matrix is calculated for a virtual receiver r out measuring the wave-field generated by multiple virtual sources r in around the virtual receiver.

For δt < 0, we observe a coherent wave propagating downward and converging to- the wave needs more time to focus, i.e. to reach the virtual receiver that is at the same location than r in . Consequently, the focusing appears below the virtual source.

During the PhD of William Lambert [START_REF] Lambert | Matrix approach for ultrasound imaging and quantification[END_REF], it was shown that ∆z in is a direct indicator of the travel time shift and therefore of the error in the velocity model.

The integrated wave velocity c can thus be estimated such as: The time focused reflection matrix can thus be used to probe axial aberrations and, more precisely, to determine the integrated wave velocity at a focusing point where an echogenic structure is present. In the following, we show how such a measurement can be done more routinely without having to compute the full time-dependent FR matrix.

IV.3 Focusing map

As explained earlier, the location of the coherent wave is perfectly known, based on the wave velocity model. It fulfills the relation of Eq. IV.2 for a homogeneous model. To follow the propagation of the coherent wave, and instead of looking at the propagation movie, we select the signals of the time focused reflection matrix that satisfy Eq. IV.2.

We obtain a focusing map defined as:

F (∆r out , r in ) = R (∆r out , r in , δt w (∆r out )) (IV.5)

We obtain an image of the focal spot for each input focusing point. Such a focusing map corresponds to the coherent wave in the coherent wave-front frame. displays the normalized amplitude along ∆y = 0. In this plot, the shift ∆z is easily measured at the maximum of the curves: ∆z in = 4.8 km and ∆z in = -6.4 km, for c 0 = 3300 m/s and c 0 = 3700 m/s, respectively. From Eq. IV.4, the integrated velocity is estimated to be ∼ 3490 m/s and 3420 m/s, for c 0 = 3300 and c 0 = 3700 m/s, respectively.

In order to find the exact correct integrated velocity, the same analysis is done for several homogeneous velocity models. The results are displayed in Until now, we have studied the time focused reflection matrix assuming that the medium is homogeneous. We have seen that the focusing matrix is a relevant observable to measure the integrated wave velocity at a given point. However, seismologists are rather interested in the local wave velocity in order to map the local mechanical properties of the subsurface. A first step towards this goal is to consider a layered velocity model instead of a homogeneous one.

IV.4 Extension to a multi-layered model

The goal of this section is to check the accuracy of the wave velocity model considered in chapter III. Similarly to the process described for the case of a homogeneous model, the time focused reflection matrix is computed in the virtual source frame by considering the layered velocity model displayed in table III. With this preliminary study, we tested the reliability of the velocity model at a specific travel time. Nevertheless, we are still far from being able to provide a map of the local wave velocity. To do so, one should first be able to measure an integrated wave velocity at any point in the subsurface. In principle, this should be possible even in random speckle regions but, in that case, the propagation movies should be extracted through a SVD of the de-scanned reflection matrix in order to extract a coherent wave over several realizations of disorder (i.e several focusing points r in ) [START_REF] Lambert | Matrix approach for ultrasound imaging and quantification[END_REF]. Once an estimation of the integrated wave velocity is obtained everywhere in the subsurface, an inverse problem should be solved in order to deduce a map of the local wave velocity [START_REF] Burdick | Upper mantle heterogeneity beneath north america from travel time tomography with global and usarray transportable array data[END_REF][START_REF] Jakovljevic | Local speed of sound estimation in tissue using pulse-echo ultrasound: Model-based approach[END_REF][START_REF] Ali | Distributed phase aberration correction techniques based on local sound speed estimates[END_REF].

IV.5 Conclusion

In the previous chapters, we have seen the impact of a wrong velocity model on the images. The resulting aberrations strongly degrade the image quality. In addition to that, a wrong velocity model prevent the positioning of the scatterers in the right location in depth, and the axial axis of the image is often a temporal axis. Based on a hypothesis on the velocity model, we presented in chapters II and III a technique that allows to overcome transverse aberrations and drastically improve the estimation of the medium reflectivity.

Here we have shown how a study of the focused reflection matrix in the time domain paves a way towards a passive tomography of the wave velocity across the medium.

While the focused reflection matrix was used in the previous chapters to build the image of the medium at the ballistic time, we introduce here a time focused reflection matrix that enables to study the focusing process. In particular, the focusing map F is studied and provides a direct estimation on the error introduced by the considered velocity model. From the shift of the focal spot with respect to the focal plane, one can conclude if the integrated velocity is overestimated, correct or underestimated.

In this preliminary study, where we tested the accuracy of the layered velocity model for the case of the NAFZ that we used in the chapter III, for a specific travel time. This easily implemented approach appears to be very promising in the field of seismology.

From a mapping of the integrated velocity across the whole field-of-view, the local velocity, that is essential to obtain accurate image of the subsurface, can be retrieved via an inversion problem.

The integrated wave velocity model can also be used as a first step prior to the imaging scheme to reduce axial aberrations and to relocate reflectors correctly along the axial dimension.

This work is of course not complete, but promising perspectives are here presented.

Further analysis and theoretical work are still required to implement this approach and build a robust passive tomography of the local velocity exhibited by body waves.

the recorded seismic data into the subsurface is done by simple matrix operations, using an estimated velocity model. A focused reflection matrix is obtained in a virtual plane inside the medium at the focusing time, i.e.time at which the incident wave reaches the focusing point. This matrix measures the impulse responses between virtual sources and receivers located within the medium. From this matrix, one can retrieve the scattering image at the considered depth. This configuration is referred to as confocal image, where the virtual source and receiver are located at the same location. But due to the mismatch between the considered velocity model and the reality, the waves suffer from phase distortions, i.e aberrations during their propagation, resulting in distorted images. To handle this problem, we introduced a novel operator, the so-called distortion matrix, that is similar to the focused reflection matrix studied in a dual basis. This matrix gives access to the distorted component of the waves in order to correct these images. This is the principle of adaptive optics and time-reversal experiments. We show that the distorted component allows to compensate for aberrations and restore the focusing inside the medium. Moreover, from the reflections matrix, the focusing quality can be assessed and the accuracy of the obtained image quantified from the common mid-point (CMP) intensity profile obtained from the anti-diagonal coefficients of the focused reflection matrix. The technique allows also to go beyond reflectivity imaging and to quantify the seismic velocity inside the subsurface. In this configuration, the focused reflection matrix is studied beyond the ballistic time.

The thesis is divided in four chapters. In the first chapter, we discussed the structure of faults and its complexity. We also reviewed conventional reflection imaging techniques adopted to image fault zones and their limitations, as well as the main notions that establish the basis of reflection imaging. We reviewed the concept of passive imaging and the reconstruction of the Green's function from noise cross-correlations. We discussed also the principle applications of these techniques that are used to image and monitor the Earth. Finally, we proposed a brief review on the first milestones of matrix imaging.

In the second chapter, we applied the matrix approach to image the San Jacinto fault using a small scale array with more than 1000 geophones. ZZ cross-correlation in the [10 20] Hz frequency band are used and a homogeneous P-waves velocity is considered. The medium displays sparse scattering properties, a few reflectors contributing to the reflected wave-field at each travel time. We successfully developed and applied the distortion matrix in this context for aberration correction to overcome small-scale velocity variations.

The singular value decomposition of the distortion matrix decomposes the field-of-view into a set of isoplanatic patches, each one being mainly associated with a scatterer inside the medium. The eigenvectors of D give access to a phase law that should be applied to compensate for the phase distortions, thereby leading to an optimization of the focusing process. The small-scale heterogeneities (in the order of the wavelength ∼ 100 m) in the first 4 km of the crust are resolved with a transverse super resolution. Such resolution is provided by the scattering of the medium that acts as a lens, thereby increasing the array effective aperture. Also the strong heterogeneities and the damage around the fault zone increase the dimensions of the imaged field-of-view, which extends beyond the size of the geophones' array. Differences between the Northwest and the Southeast of the fault were reported with an intense localized damage zone in the Southeast. Also, lateral variations of scattering properties within and around the surface fault traces are observed.

In the third chapter, the approach is successfully applied in the case of the North Ana- andB2). Further analysis will be carried out in this direction to better understand the relation between the scattering, the earthquakes' distribution and the damage intensity.

A continuation of the current work is to further study the time focused reflection matrix in order to build a 3D velocity model for the NAFZ region. In the future, we may also expect to improve the construction of the body waves Green's function from noise correlations and see its implication on the final images. Also, it could be interesting to combine the results obtained from both vertical and horizontal cross-correlations. Another aspect is to consider P and S-waves conversions in the imaging scheme.

Another interesting future opportunity is to test the method with other array configurations. A possible application is to apply the reflection matrix approach to image the medium between two small seismic networks. Another aspect is a monitoring application in which we study the temporal changes of the scattering inside a medium with the presented approach.

Another common issue in imaging that was not addressed in this thesis is the multiple scattering of waves. These events could appear as a coherent noise in the reflection matrix and may result in spurious arrivals and ghosts in the reflectivity images. But it was shown that, by virtue of the matrix approach, it is possible to separate single and multiple scattering contribution (Aubry and Derode, 2009b;[START_REF] Badon | Smart optical coherence tomography for ultra-deep imaging through highly scattering media[END_REF]. This is done by applying a confocal filter to the reflection matrix prior to any aberration correction [START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]. Alternatively, if precisely compensated, the multiple scattering contribution can be exploited for deep imaging and super-resolution.

Last but not least, multiple reflections were also not studied in the present thesis.

Specific multiples' attenuation techniques, such as Marchenko method [START_REF] Wapenaar | Marchenko imaging[END_REF][START_REF] Staring | Source-receiver marchenko redatuming on field data using an adaptive double-focusing method[END_REF][START_REF] Wapenaar | Marchenko redatuming, imaging and multiple elimination, and their mutual relations[END_REF], could be applied to the measured reflection response prior to aberration correction in media where multiples can not be neglected.

Together, multiples' attenuation techniques and the distortion matrix concept can lead to improved images of the subsurface.

  zones structure, imaging and detection . . . . . . . . . . . . . I.1.1 Ground-penetration radar (GPR) . . . . . . . . . . . . . . . . . I.1.2 Electrical resistivity tomography (ERT) . . . . . . . . . . . . . I.1.3 Seismic methods . . . . . . . . . . . . . . . . . . . . . . . . . I.1.3.1 Active seismic imaging . . . . . . . . . . . . . . . . I.1.3.2 Passive seismic imaging . . . . . . . . . . . . . . . I.2 Seismic ambient noise . . . . . . . . . . . . . . . . . . . . . . . . . . I.2.1 Noise origin . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.2.2 Noise cross-correlations and Green's function . . . . . . . . . . I.2.3 Analogy with the time reversal experiment . . . . . . . . . . . I.3 Migration techniques . . . . . . . . . . . . . . . . . . . . . . . . . . I.3.1 Wave-field extrapolation . . . . . . . . . . . . . . . . . . . . . I.3.1.1 Kirchhoff and Rayleigh integrals for wave-field extrapolation (summation approach) . . . . . . . . . . . I.3.1.2 Extrapolation in the wavenumber-frequency domain . I.3.2 Active seismic imaging: Redatuming . . . . . . . . . . . . . . I.3.3 Limits of seismic imaging . . . . . . . . . . . . . . . . . . . . 38 I.3.3.1 Multiple scattering . . . . . . . . . . . . . . . . . . . 38 I.3.3.2 Multiple reflections (multiples) . . . . . . . . . . . . 39 I.3.3.3 Phase distortions . . . . . . . . . . . . . . . . . . . . 40 I.4 Reflection matrix approach of wave imaging . . . . . . . . . . . . . 44 I.4.1 The time reversal concept . . . . . . . . . . . . . . . . . . . . 44 I.4.2 Decomposition of the time reversal operator (DORT) . . . . . . 47 I.4.3 Matrix imaging . . . . . . . . . . . . . . . . . . . . . . . . . . 50 I.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 I.1 Fault zones structure, imaging and detection Earthquakes constitute one of the main natural disasters in the world. With a magnitude that can exceed M>8, earthquakes can cause a tremendous number of fatalities and property losses. They can also trigger tsunamis, landslides and fire. The 21 st century has witnessed many devastating earthquakes such as the Indian Ocean Earthquake and Tsunami, 2004 (M9.1) and Haiti earthquake, 2010 (M7) with more than 200, 000 fatalities each.

  methods and show few examples of applications aiming to image fault zones structures. I.1.3.2.1 Surface wave tomography Surface-wave tomography relies on the extraction of Rayleigh and Love waves from the cross-correlation of noise records. The group and phase velocity dispersion curves are then inverted in order to obtain shear wave velocity images of the subsurface. This technique has been successfully applied in several regions, taking advantage of the increasing deployment of dense arrays (Fig. I.1)

Figure I. 1 :

 1 Figure I.1: Ambient noise tomography for California showing 7.5-s-period Rayleigh waves group velocity map. The figure is taken from Shapiro et al. (2005).

Figure I. 2 :

 2 Figure I.2: Comparison between migrated reflection image obtained from reflections retrieved from noise correlations (a) and from active-source gather (b). Coinciding structures are highlighted in red. The figure is taken from Draganov et al. (2009).

Figure I. 3 :

 3 Figure I.3: Surface waves velocity variation obtained from ambient noise crosscorrelations observed along the San Andreas fault after San Simeon (2003) and Parkfield (2004) earthquakes. The figure is taken from Brenguier et al. (2008b).

  Using continuous measurements of ambient seismic noise, Brenguier et al. (2008b) reported a decrease by 0.04% in seismic velocities along the San Andreas fault after the San Simeon (2003) and Parkfield (2004) earthquakes (Fig. I.3).

  2) with σ 2 the power spectrum density of the noise and δ(t) the Dirac distribution and δ ij the Kronecker symbol. The noise sources surround two receivers at A and B recording the displacement (Fig. I.4a),such as:

  5) Now let us consider the following time reversal experiment. A point source at A is surrounded by a time-reversal cavity (TRC) of N devices (Fig. I.4b). Each time-reversal device S i records the wave-field emitted by A, time-reverses it (G AS i (t) → G AS i (-t)) and sends it back into the medium (Fig. I.4c). The recorded wave-field at B is equal to the convolution between the emitted signals and the response of the medium such as:

  between A and B (Fig. I.4c). The time-reversed wave-field can also be expressed by describing the wave-field at B at each step of the time reversal process. During the forward step (Fig. I.4b), at time t = 0, A sends a pulse that propagates everywhere in the medium including in B where the field received is G AB (t), and is eventually recorded on every point S i . After time reversal (Fig. I.4c), the wave should exactly go backwards: it should hit B first and refocus on A at time t = 0, which implies that the field received in B at times t < 0 is exactly G AB (-t), the anti-causal part of the Green's function. Once the pulse has refocused on A, it does not stop but diverges again from A and gives rise, at times t > 0, to G AB (t) in B (Fig. I.4d). The response at B thus consists of the sum of the causal and the acausal part of the Green's function, according to:

Figure I. 4 :

 4 Figure I.4: Analogy between the cross-correlation between A and B and the time reversal experiment. (a) The correlation between wave-fields emitted by S i and recorded at A and B. (b) Response recorded at S i if there is a source at A. (c) The time-reversed responses recorded at S i are emitted back into the medium and focus back at A at t = 0 after hitting B at t < 0. (d) The wave-field diverges again from A and hits again B at t > 0. The wave-field recorded at B thus corresponds to the Green's function between A and B and is equal to the cross-correlation of the random wave-field between these two points.

  extrapolation techniques, one should differentiate between two approaches considered in seismic imaging: modeling and migration. These two techniques are closely related and are based on the wave equation. Modeling, which involves forward extrapolation of downgoing and upgoing waves, is a simulation of the recorded wave-field, i.e. simulation of wave propagation inside the medium of interest, whereas migration is essentially based on the inverse (backward) propagation of the recorded wave-field. It aims at back-propagating the wave-field and provides an estimation of the medium's reflectivity. In both cases, wave-field extrapolation can be expressed through Kirchhoff integral or in the wavenumber-frequency domain. I.3.1.1 Kirchhoff and Rayleigh integrals for wave-field extrapolation (summation approach)

  surface S, allow to compute the wave-field at any point A inside the volume V enveloped by S. This version of the Kirchhoff integral is used for forward extrapolation in which the wave-field is calculated at later times than when it was measured on the surface S (Fig. I.5a). That's why this is the causal Green's function that appears in Eq. I.13.(Fig. I.5a). Now instead of G(r, r A , ω), we consider the anti-causal Green function G * (r, r A , ω)

  Figure I.5: Kirchhoff and Rayleigh integrals. (a) Forward (Eq. I.13) and backward (Eq. I.14) extrapolation with Kirchhoff integral. Assuming a source outside the volume V , the wave-field and its normal derivative are recorded at surface S. The blue part of the surface corresponds to when the wave-field is entering into the volume, and the red part represents the surface where the wave-field travels outside the volume. G and G * denote the causal and anti-causal Green's functions, respectively. (b) Forward extrapolation using the Rayleigh integrals (Eqs. I.15 and I.16). The wavefield generated by a source at the surface and scattered by the heterogeneities of the medium is recorded on a planar surface S. The causal Green's function G is used to obtain the wave-field at point A(r A ). (c) Backward extrapolation using the Rayleigh integrals (Eqs. I.17 and I.18). A secondary source generates a wave-field that propagates inside the medium. The wave-field produced by the secondary source is back-propagated and estimated at the virtual receiver A using the anti-causal Green's function G * .

  Figure I.5: Kirchhoff and Rayleigh integrals. (a) Forward (Eq. I.13) and backward (Eq. I.14) extrapolation with Kirchhoff integral. Assuming a source outside the volume V , the wave-field and its normal derivative are recorded at surface S. The blue part of the surface corresponds to when the wave-field is entering into the volume, and the red part represents the surface where the wave-field travels outside the volume. G and G * denote the causal and anti-causal Green's functions, respectively. (b) Forward extrapolation using the Rayleigh integrals (Eqs. I.15 and I.16). The wavefield generated by a source at the surface and scattered by the heterogeneities of the medium is recorded on a planar surface S. The causal Green's function G is used to obtain the wave-field at point A(r A ). (c) Backward extrapolation using the Rayleigh integrals (Eqs. I.17 and I.18). A secondary source generates a wave-field that propagates inside the medium. The wave-field produced by the secondary source is back-propagated and estimated at the virtual receiver A using the anti-causal Green's function G * .

Figure I. 6 :

 6 Figure I.6: Double focusing operation in ultrasound imaging can be interpreted as a wave extrapolation in seismic imaging. (a) Focusing at emission. (b) Focusing at reception.

  Figure I.8: Schematic description of the single scattering regime, multiple scattering regime and multiple reflections (multiples). (a) Sketch of the single scattering regime where the seismic waves are only reflected once by each heterogeneity. (b) Sketch of the multiple scattering regime in which waves are reflected several times through their travel between the source and the receiver. The scattering mean free path s is represented. (c) Cartoon representing surface (red arrows) and internal (green arrows) multiples observed in a marine setting.

  Adaptive focusing principle in seismic imaging is illustrated in Fig. I.10. The wavefront distortions are measured by the array of geophones (Fig. I.10a). The phase conjugate of the aberration is applied in order to focus on a specific point inside the medium (Fig. I.10b

  Figure I.10: Adaptive focusing principle in seismic imaging for aberration correction. (a) The wave-front distortions are measured by the surface array. (b) Phase distortions compensation consists of applying the phase conjugate of the aberration in order to focus in a specific isoplanatic patch inside the medium. (c) The same phase law is applied to any focusing point inside the same isoplanatic region. (d) Outside this specific region, the phase law does not hold anymore.

  Figure I.11: Time-reversal cavity. (a) The distorted wave-field generated by a source inside a heterogeneous medium is recorded by an array of transducers over a closed surface. (b) The recorded signals are time-reversed and back-propagated inside the medium. The emitted wave-field refocus back at the original source position. Figures are adjusted from[START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF].

  Figure I.12: Principle of a time reversal mirror. (a) A wave-front is transmitted from the array of transducers through a heterogeneous medium behind which a single target is hidden. (b) The target generates a back-scattered wave-field that is distorted while propagating through the medium and recorded by the array of transducers. (c) The array sends back the time-reversed wave-field into the medium, thereby compensates for the distortions induced by the heterogeneous medium and finally focus ideally on the target. Figures are adjusted from[START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF].

  Figure I.13: DORT applied to retrieve 3D images of Erebus volcano, Antarctica. (a) Random speckle confocale image. (b) Coherent combination of the response matrix eigenvectors. Figures are adjusted from[START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF].

  Figure I.13: DORT applied to retrieve 3D images of Erebus volcano, Antarctica. (a) Random speckle confocale image. (b) Coherent combination of the response matrix eigenvectors. Figures are adjusted from[START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF].

  Jacinto fault [Chap. II] and the North Anatolian Fault [Chap. III]. Imaging the scattering structure of the crust, and even the upper mantle, is made possible through body waves reflections that are retrieved from ambient noise cross-correlations. The reflection matrix also provides a new promising tool to measure the local seismic velocity [Chap. IV]. Contents II.1 A Distortion Matrix Framework for High-Resolution Passive Seismic 3D Imaging: Application to the San Jacinto Fault Zone, California . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.1.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.1.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.1.3 Reflection matrix . . . . . . . . . . . . . . . . . . . . . . . . . II.1.3.1 Response matrix between geophones . . . . . . . . . II.1.3.2 Focused reflection matrix . . . . . . . . . . . . . . . II.1.4 Distortion matrix . . . . . . . . . . . . . . . . . . . . . . . . . II.1.4.1 Reflection matrix in a dual basis . . . . . . . . . . . II.1.4.2 Definition and physical interpretation of the distortion matrix . . . . . . . . . . . . . . . . . . . . . . . II.1.5 Time reversal analysis of the distortion matrix . . . . . . . . . . II.1.5.1 Output distortion matrix and isoplanatic patches . . . II.1.5.2 Input distortion matrix and super-resolution . . . . . II.1.5.3 Normalized correlation matrix and compensation for high-order aberrations . . . . . . . . . . . . . . . . . II.1.6 Three-dimensional imaging of the San Jacinto Fault zone . . . . II.1.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.1.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.1.9 Supplementary Material . . . . . . . . . . . . . . . . . . . . . S1 Prior filtering of the reflection matrix . . . . . . . . . S2 Singular value decomposition of the distortion matrix S3 High-order eigenstates of the D-matrix. . . . . . . . S4 Normalized correlation matrix . . . . . . . . . . . . S5 Choice of the wave velocity model. . . .

  particular, Zigone et al. (2019) used ambient noise cross-correlations in the 2-35 Hz frequency range to derive a velocity model in the top 100 m with a resolution of 50 m.

  Figure II.1: (a) Map of the 1108 (10 Hz) geophones installed in a 600 m × 700 m configuration above the Clark branch (red lines) of the San Jacinto Fault (Southern California).Each row along the x-direction is composed of ∼55 sensors with a pitch of 10 m, and the nominal separation between the rows in the y-direction is 30 m. Seismic ambient noise was recorded over more than one month, in May-June 2014. (b) Adaptive focusing at emission and reception on two points r in and r out of the focal plane (z = c 0 t/2) yields the impulse response between virtual geophones placed at these two points. The same operation is repeated for any couple of points in the focal plane and yields the focused reflection matrix R.

  II.1b). Each coefficient of the focused reflection matrix R(z, f ) involves pairs of virtual geophones, r in = (x in , y in ) and r out = (x out , y out ), which are located at the same depth z (Fig.II.1b).For broadband signals, ballistic time-gating can be performed to select only the echoes arriving at the ballistic time t B in the focused basis(Lambert et al., 2020a): t B = ( r in -u in + r out -u out )/c 0 . Under a matrix formalism, this time-gating can be performed by means of a coherent sum of R(f ) over the frequency bandwidth ∆f = 10 Hz.

Figure

  Figure II.2a displays one example of the broadband focused reflection matrix R at depth z =3600 m. In the case of SJFZ, it appears that a part of the backscattered energy is still concentrated in the vicinity of the diagonal of the focused reflection matrix at z = 3600 m (Fig. II.2a); this is very different from the Erebus volcano for which the reflection matrix displayed a fully random feature[START_REF] Blondel | Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica[END_REF]. This indicates that

Figure

  Figure II.2: (a) Focused reflection matrix R in the focal plane at depth z =3600 m. (b) Imaging PSF deduced from the antidiagonal of R in (a) whose common mid-point exhibits the maximum confocal signal. The white circle accounts for the theoretical transverse resolution cell imposed by the geophone array aperture. (c) Confocal image extracted from the diagonal of R in (a). (d) Vertical slice of the 3D confocal image obtained by combining the diagonal of R at each depth. The slice orientation is chosen to be normal to the fault. The color scale on the bottom left is linear and applies to panels (a,c). The color scale in panels (b,d) is in dB.

Figure

  Figure II.2c displays the 2D confocal image built from the diagonal of the reflection matrix in Fig. II.2a at ballistic time t B = 4.6 s, hence at an effective depth z = c 0 t B /2 = 3600 m. Some scattering structures seem to arise at different locations along the fault but confocal imaging is extremely sensitive to aberration issues. One thus has to be very careful about the interpretation of a raw confocal image. This observation is confirmed by Fig. II.2d that displays a cross-sectional view of the SJFZ underground. Each speckle grain in this image occupies a major part of the field-of-view. Hence, aberrations seem to be pretty intense at large depths (beyond 1500 m) and the inner structure of the SJFZ cannot be deduced from a basic confocal image.

Figure

  Figure II.2b shows the corresponding common midpoint intensity profile for the matrix R displayed in Fig. II.2a. It shows a significant spreading of energy over off-diagonal coefficients of R. This effect is a direct manifestation of the aberrations sketched in Fig. II.1b. Indeed, in absence of aberration, all the back-scattered energy would be

δ

  Figure II.3: Time reversal analysis of the distortion matrix. (a) In the matrix imaging scheme, each point in the focal plane is probed by means of a focused beam at input and output. (b) A far-field projection of the focused reflection matrix [Eq. (II.10)] yields the dual-basis matrix R. (c) By subtracting from each reflected wave-front the geometrical phase law which would be obtained for a perfectly homogeneous medium of wave velocity c 0 [Eq. (II.13)], a distortion matrix D is obtained. D is equivalent to a reflection matrix but with a static input PSF H(r) scanned by moving scatterers (Eq. II.15). (d) The SVD of the distortion matrix enables the synthesis of a virtual coherent reflector of scattering distribution |H(r)| 2 (see Supplementary Section S2) and an estimation of the aberration phase transmittance H(p) for each isoplanatic patch p in the field-of-view [Eq. (II.18)]. (e) This estimation can be refined by considering, in the second part of the process, the SVD of the normalized time reversal operator δ Ĉ(p) (see Supplementary Section S4). This operation makes the virtual reflector point-like and the estimation of H(p) more precise. (f) The phase conjugate of H(p) yields the focusing law to scan the corresponding isoplanatic patch p and synthesize a novel focused reflection matrix R p [Eq. II.27].

(

  Fig. II.3b) and which can be directly extracted from the reference matrix T 0 , and (ii) a distorted component due to the mismatch between the propagation model and reality (Fig.II.3c). The principle of our approach is to isolate the latter contribution by subtract-

  r in ). (II.14) The matrix D out = [D(k out , r in )] connects any input focal point r in to the distorted component of the reflected wave-field in the far-field. Removing the ideal phase law predicted by our propagation model from the reflected wave-field in the Fourier plane as done in Eq. (II.13) amounts to a change of reference frame. While the original reflection matrix is recorded in the Earth's frame (static underground scanned by the input focusing beam, see Fig. II.3b), the matrix D out can be seen a reflection matrix in the frame of the input focusing beam (moving subsoil insonified by a static focusing beam, see Fig. II.3c).

Figure

  Figure II.4: Singular value decomposition of the distortion matrix D out at time t B = 4.8 s and depth z =3600 m. (a) Plot of the normalized singular values σp . The two first eigenstates form the relevant signal subspace for imaging. The moduli of the corresponding input eigenvectors, V

  Figure II.4a confirms this prediction by displaying the normalized singular values, σi = σ i / N j=1 σ 2 j , of D out at depth z = 3600 m (t B = 4.8 s).

  onto the corresponding isoplanatic patch p. Figures II.4b and c confirm this assertion by showing that V (1) in and V

  Output phase distortion correction at time t B = 4.8 s and depth z =3600 m. (a,b) Output focusing laws derived from the phase conjugation of normalized eigenvectors Û(1) out and Û(2) out , respectively. (c,d) Imaging PSF extracted from the main antidiagonals of the focused reflection matrices R 1 and R 2 , respectively (Eq. II.20). The color scale is in dB. (e,f) Corresponding confocal images I 1 and I 2 , respectively. ulated by the autocorrelation function H(p) in H(p) in (see Fig II.3d). This last term is a manifestation of the finite size δ (p) in of the virtual reflector (Fig. II.3d) that tends to limit the support of the eigenvector U (p) out to ∆k (p) out ∼ λz/δ (p) in . Figures II.4d and e confirm this theoretical prediction by showing the modulus of U (1) out and U (2) out , respectively. Both singular vectors cover a restricted and different angular domain in the Fourier space (∆θ (p) out = ∆k (p) out /k 0 ∼ 10 o ). To circumvent this issue, one trick is to use only the phase of these eigenvectors U (p) out (Figs. II.5a and b) by considering the normalized vector Û(p) out , such that

Figure

  Figure II.6: Input phase distortion correction at time t B = 4.8 s and depth z =3600 m. (a,b) Modulus of input eigenvectors U (1) in and U (2) in . (c,d) Output focusing laws derived from the phase conjugation of normalized eigenvectors Û(1) in and Û(2) in , respectively. (e,f) Imaging PSF extracted from the main antidiagonals of the focused reflection matrices R 1 and R 2 , respectively [Eq. (II.23)]. The color scale is in dB. (g,h) Corresponding confocal images I 1 and I 2 , respectively.

  Figs. II.6g and h at depth z = 3600 m. Their comparison with the previous images (Figs. II.5e and f, respectively) illustrates the benefit of a simultaneous aberration correction at input and output. While I 1 yields a refined image of the complex structure lying in the first isoplanatic patch (Fig. II.6g), the second eigenvector I 2 now clearly highlights a coherent reflector belonging to a second isoplanatic patch on the west of the fault (Fig. II.6h).

Figure

  Figure II.7: Correction of residual phase distortions by making the virtual scatterer pointlike (t B = 4.8 s, z = 3600 m). (a,b,c,d) Output and input focusing laws derived from the phase conjugation of the normalized eigenvectors δ Û(1) out , δ Û(2) out , δ Û(1) in and δ Û(2) in , respectively. (e,f,g,h) Imaging PSFs extracted from the main antidiagonal of the focused reflection matrices R p after application of the additional focusing law displayed in (a,b,c,d), respectively. The color scale is in dB. (i,j) Final confocal images built from the diagonal of the final reflection matrices R p , respectively [Eq. (II.27)].

  out , k out )|. (II.25) As illustrated by Fig. II.3e, this operation makes the virtual reflector point-like. Such a matrix is actually equivalent to the time reversal operator associated with a point-like reflector at the origin [Lambert et al. (2020b), see Supplementary Section S4]., In that case, the matrix δ Ĉ(p) out is ideally of rank 1 and its eigenvector δU in Fig. II.7a and b, respectively. Compared to the Fresnel zone fringes exhibited by the first-order aberrations in U (p) out (Fig. II.5a and b), δU (p) out displays higher-order aberrations, thereby leading to more complex focusing law. Updated focused reflection matrices R p are then deduced: R p = T † 0 × T 0 • δU enables an estimation of the imaging PSF over each isoplanatic patch (see Figs. II.7e and f). This second correction step drastically reduces the spatial extension of the PSF width compared to the previous step (Figs. II.6e and f). The incoherent background is below -30 dB beyond the theoretical transverse resolution cell (δ 0 = 600 m at depth z = 3600 m). Moreover, the FWHM of the PSF is now of the order of the wavelength: δ (p)

  for each isoplanatic patch p. Additional input focusing laws δU (p) * in are extracted through the SVD of the normalized correlation matrix δ Ĉ(p) in . The corresponding wave-fronts are displayed in Figs. II.7c and d. They exhibit phase fluctuations that are reduced compared to their output counterpart (Figs. II.7a and b), an indication that our approach gradually converges towards a finite aberration phase law. The normalized eigenvectors, δ Û(p) in , can be used as input focusing laws that should be applied to D (p) in to compensate for residual input aberrations. The resulting images of the subsoil reflectivity are displayed in Figs. II.7 i and j at z = 3600 m over the two main isoplanatic patches. Their comparison with the previous images (Figs. II.6g and h, respectively) illustrates the benefit of iterating the aberration correction process. The first isoplanatic patch highlights the presence of three distinct reflectors that arise in the vicinity of the fault's surface traces (Fig. II.7i). The second isoplanatic patch yields a highly-resolved image of a coherent reflector on the west of the Clark fault (Fig. II.7i). Compared to the previous image in Fig. II.7j, the contrast is drastically improved. This observation can be understood by looking at the corresponding PSFs (Figs. II.7g and h). The incoherent background is below -35 dB beyond the theoretical transverse resolution cell (δ 0 = 600 m at depth z = 3600 m). Last but not least, the FWHM of the PSF nearly reaches the diffraction limit δ (p) in ∼ 80 m. Remarkably, the observed lateral resolution is nearly one-eighth of the theoretical limit based on the geophone array aperture. As mentioned in the previous section, the reasons for this spectacular resolution will be discussed in Sec. II.1.7. The comparison of these final images with the original one (Fig. II.2c) illustrates the relevance of a matrix approach for deep seismic imaging. In the next section, the threedimensional structure of the SJFZ is revealed by combining the images derived at each depth. At last, based on the derived high-resolution 3D images, a structural interpretation of the SJFZ is provided.

Figure

  Figure II.8: Vertical slice of the 3D confocal image obtained from stacked corrected images derived at each depth. The slice orientation is chosen to be normal to the fault plane. From this image, the fault location can be circumscribed and is represented by the shaded area. The color scale is in dB.

Figure

  Figure II.8 shows a slice of the final 3D image I M (r, z) with the same orientation as the original cross-sectional view displayed in Fig. II.2d. While the raw image I is completely blurred, the gain in horizontal resolution provided by the matrix aberration correction process reveals an image of the SJFZ subsoil with a refined level of details. Based on the surface traces displayed in Fig. II.1a, the fault is expected to spread orthogonally to the x-axis. The choice of the profiles' orientation in Figs. II.2d and II.8 is dictated by our willingness to highlight the fault blocks on the right and left side of the slip interface. The SJFZ is a major continental strike-slip fault system that accumulated through history a significant slip of tens of kilometers. Such large fault systems induce zones of strongly damaged materials (Ben-Zion and Sammis, 2003, and references therein) and the damage is expected to be pronounced in the shallow crust. Thanks to the distortion matrix correction, Fig. II.8 shows the two main signatures of the fault in diffraction imaging. First, it reveals a damage zone associated with a dense distribution of scatterers. The damage volume extends from the surface to 1000 m below with a section decreasing in depth. Second, beyond a depth of 1000 m, the matrix image displays strong echoes associated with the presence of sub-horizontal reflectors. The corresponding strata layers are located at different depths on both sides of the fault. This discontinuity seems to indicate the fault location in depth (see the shaded area in Fig. II.8), although the damage area isno longer discernible beyond 1000 m. This example thus shows both the efficiency of our approach to cope with phase distortions and the potential of passive seismic imaging to study the fine structure of active faults in depth.In this section, only a slice of the 3D volume has been shown to illustrate the drastic improvements granted by a matrix approach of seismic imaging (Fig.II.8). Its comparison with the raw confocal image highlights the benefit of a drastic gain in horizontal resolution (Fig.II.2c). A more detailed interpretation of the subsurface properties will be the focus of a future study. Indeed, a detailed interpretation of the obtained 3D image and a confrontation with previous studies may provide information on the structural properties of the fault zone. More generally, the high transverse resolution and penetration depth of our matrix imaging method will play an important role in detecting active faults,

  Fig. II.8, the horizontal strata layers are imaged on each side of the fault over a field-of-view much larger than the geophone array dimension. Two reasons can account for this surprising result. First, the scattering between each layer may not be only specular but also induced by a distribution of localized inhomogeneities at each layer interface. The transverse images shown at depth z=3600 m in Fig. II.7 confirm this by highlighting the presence of four localized scattering structures in the first isoplanatic patch [Fig. II.7(i)]. Second, the large extension of the image in Fig. II.8 can also be due to the scattering induced by the strongly heterogeneous damage area. The effective geophone array aperture can thus be increased by scattering and so is the imaged area. Last but not least, the matrix image of the SJFZ fault shows a transverse superresolution highlighted by the imaging PSF displayed in Fig. II.7(e)-(h). A first reason for this striking result could be an overestimation of the wave speed c at shallow depth.

  FigureS1adisplays the response matrix K(t) at time lag t = 0.50 s. Surprisingly, this matrix is dominated by a predominant signal along its diagonal. To understand the origin of this effect, let us first recall that each coefficient K(u out , u in , t) of the response matrix is computed from the cross-correlation of the seismic noise wave-field ψ(u, τ ) recorded

  S3)where ρ(p) is the overall patch reflectivity and the symbol stands for a correlation product. The correlation function in Eq. (S3) results from the Fourier transform of the scattervirtual scatterer shown in Fig. II.3d. Its support is the correlation width ∆k out of the aberration phase transmittance that scales as the inverse of the spatial extension δ scatterer was point-like, this correlation function would be constant and the matrices C (p)

Fig

  Fig. II.2b. Among these Q eigenmodes of C (p) out , only the first eigenvector U (p,1)

  Fig. II.4a of the accompanying paper, the SVD of D out shows four main singular values on top of a "noise" continuum. However, only the two first eigenstates (Fig. II.5)are then used to compensate for aberrations over two main isoplanatic patches. Fig.S2displays the two following eigenstates of rank 3 and 4. While the third input singular vector, V

  . II.4c). Similarly, the ouptput singular vectors U(3) out (Fig.S2c) and U (4) out (Fig. S2d) map onto the same angular domain as U Singular value decomposition of the distortion matrix D out at time t B = 4.8 s and depth z =3600 m. (a,b) Modulus of input eigenvectors V in,3 and V (4) in,4 . (d,e) Modulus of output eigenvectors U out,3 and U out,4 . (e) Hadamard product between U out,3 and U * out,1 . (f) Hadamard product between U out,4 and U * out,2 . (Fig. II.4d) and U (2) out (Fig. II.4e), respectively.

  single isoplanatic patch p before extending those predictions to the general case of transverselyvarying aberrations. Over each isoplanatic patch p, the corresponding correlation ma-

  to map onto the corresponding isoplanatic patches. Fig. II.4b and c indeed show that the two first eigenstates map onto two disjoint areas. This observation proves that the two first eigenstates are associated with different isoplanatic patches. On the contrary, the third and fourth eigenstates, V

  The link between the first and third eigenstates is confirmed by Fig.S3ethat displays the real part of the Hadamard product U

Figure S3 :

 S3 Figure S3: Imaging PSF deduced from the antidiagonal of R whose common mid-point exhibits the maximum confocal signal. This imaaging PSF is shown at different times of flight (t B = 1.05, 1.5 and 2.65 s, from top to bottom) for different wave velocity models: c 0 = 1000 m/s (a), c 0 = 1500 m/s (b), c 0 = 2000 m/s (c). Each panel displays the modulus of the reflected wave-field normalized by its maximum. The white circle accounts for the theoretical resolution cell (disk of radius δ 0 ) imposed by the geophone array aperture.

  FigureS4shows a slice of this image for different wave velocity models (c 0 = 1000, 1500 and 2000 m/s). The slice orientation is the same as the b-scans displayed in the accompanying paper. While the raw image I 0 is completely blurred (Fig.2cof the accompanying paper), a gain in resolution is provided by the matrix aberration correction process. However, the image quality strongly depends on the choice of the wave velocity model. This sensitivity can be explained by the fact that axial aberrations are not tackled by our matrix approach so far. Hence, an optimized wave velocity model allows to properly capture all the echoes back-scattered at the focal depth z in R(z). The comparison between Fig.S4a and bconfirms that a model with c 0 = 1000 m/s clearly underestimates the actual body wave velocity in the depth range considered in this work. While Fig. S4b (c 0 = 1500 m/s) clearly highlights strata layers at different depths on both sides of the fault, Fig. S4a (c 0 = 1000 m/s) shows a blurred image of the subsoil.

  Figure S4: Vertical slice of the 3D matrix images computed under the matrix approach for different seismic wave velocity models: c 0 = 1000 m/s (a), c 0 = 1500 m/s (b), c 0 = 2000 m/s (c). The slice orientation is chosen to be normal to the fault plane. The color scale for each image is in dB.

Figure

  Figure II.9: (a) Topographic map of Trifurcation area of San Jacinto fault zone. The red square marks the studied area. (b) Map of the geophones of the dense array at SGB site (black dots). The Clark fault traces are represented by the red lines. The blue lines indicate the locations of the cross-sections represented in Fig. II.10 and II.11.

  . II.9) is associated with branching of the Clark fault in the Anza section into three major faults: a continuation of the Clark fault and the Buck Ridge and Coyote Creek faults. The Trifurcation area has a broad zone of high seismicity rates that include five earthquakes with magnitudes around 5 since 2001. The geometrical properties of the seismicity in the trifurcation area are very complex and consist of a diffuse pattern in the top 5 km that changes to more localized structures dipping to the NE below 6 km, along with zones of seismicity that are orthogonal to the main strike of the SJFZ[START_REF] Ross | Abundant off-fault seismicity and orthogonal structures in the san jacinto fault zone[END_REF].

  taken from the 3D scattering volume. Those images are obtained after correcting the aberrations induced by the mismatch between the velocity model used to perform the focusing operation and the data. Although the images shown in Fig. II.10 and Fig. II.11

  Fig. II.10 shows vertical slices of the 3D volume with a close-up view of the first 400 m: two slices perpendicular to the fault (AA', BB') and two slices parallel to the fault traces (CC', DD'). For the sake of simplicity, the cross-sections are labeled North, South, West and East, respectively. The slices are plotted in logarithmic dB scale and reveal the backscattered intensity rising from highly reflective features detected through the aberration correction process. The location of the each cross-section is indicated by blue lines in Fig. II.9b. In the first two vertical slices, the white lines refer to the location of the three main sub-parallel strands represented in red in Fig. II.9b.

Figure

  Figure II.10: Shallow cross-sections of the 3D scattering volume. Vertical slices oriented perpendicular to the fault traces. North/South denote Northwest and Southeast. The main fault strands are represented by the bold white lines. (a) Profile A-A', (b) Profile B-B'. Vertical slices oriented parallel to the fault traces. West/East denote Southwest and Northeast. (c) Profile C-C', (d) Profile D-D'. The color scale is in dB. The white dashed lines correspond to the spatial extension of the array.

1 (Figure

 1 Figure II.11: Deep cross-sections of the 3D scattering volume. Vertical slices oriented perpendicular to the fault traces. (a) Profile E-E' Northwest (NW) of the array, (b) Profile F-F' Southeast (SE) of the array. Vertical slices oriented parallel to the fault traces. (c) Profile G-G' Southwest (SW), (d) Profile H-H' Northeast (NE). The color scale is in dB. The white dashed lines correspond to the spatial extension of the array.

Figure

  Fig.II.12a2. The linear regression is also plotted (dashed lines) to visualize the slope of the energy decay. We only show the energy decay corresponding to the images after the correction process. The first noticeable feature is that the fault zone is associated with higher intensity values resulting from non-specular energy transmitted from the localized damage. This is consistent with the overall amplification of seismic waves in low velocity fault zone layers (e.g.[START_REF] Kurzon | Ground motion prediction equations in the san jacinto fault zone: Significant effects of rupture directivity and fault zone amplification[END_REF]). The energy in the fault zone with a steep slope distinguishes itself from the surrounding western and eastern regions where the decay slopes are more gradual. The fault zone is defined by major fracturing and crushed rocks, and consequently it is characterized by a rapid energy decay. In the neighboring regions where the damage is more distributed, less attenuation is observed in comparison with the fault zone. These results are consistent with significantly lower values of attenuation coefficients generally found within fault zones by waveform fitting of trapped waves (e.g.[START_REF] Lewis | Shallow seismic trapping structure in the san jacinto fault zone near anza, california[END_REF];[START_REF] Qin | Internal structure of the San Jacinto fault zone in the trifurcation area southeast of Anza, California, from data of dense seismic arrays[END_REF]).We also compare the intensity decay between the SE and NW sections 4 km below the array. Fig. II.12b2 shows plots of the backscattered intensity averaged across the two regions delineated in Fig. II.12b1. We choose to plot the measured intensity values without the linear regression. The discrepancy observed in the backscattered intensity distribution primarily reflects the difference in subsurface structure between the two regions that was

FigFigure

  Fig. II.13 summarized schematically the obtained imaging results for both the inner and outer damage zones in the area. The fault traces in Fig. II.9a suggest a broader and less intense fault zone in the NW than in the SE. Surface observations are consistent with the fact that the main fault and the surrounding core damage are more localized in the SE where the principle slip zone is delineated by intense damage and fracturing extending down to 3 to 4 kilometers. In the NW, the scattering in the FZ is only observed in the first kilometer indicating a shallow less intense and diffuse damage. The results in Figs. II.11c and d show that the damage distribution is more complex to the SW side of the fault exhibiting a more pronounced outer damage in the SW than in the NE side of the fault.The lateral variations of the fault structure between the NW and SE are consistent with the transition from the Anza section of the SJFZ associated with a single major fault trace, to a complex fault zone in the trifurcation area with several traces at the surface. The high

Figure B :

 B Figure B: Impact of specular and diffuse reflection on the field-of-view of the confocal image. (a) Smooth mirror-like surface. (a.1) Sketch showing the reflected waves traveling with the same angle as the incident waves. Some of the reflected rays are not captured by the array. (a.2) Synthetic confocal image of a specular reflector at z = 3600 m. The white dashed square corresponds to the spatial extension of the array. (b) Random scattering medium. (b.1) Sketch displaying distributed heterogeneities at depth. An incident ray is scattered in all directions enabling a confocal image that extends well beyond the lateral dimensions of the geophone array. (b.2) Mean synthetic confocal image at z = 3600 m estimated from several iterations obtained in the case of a random medium. The white dashed square corresponds to the spatial extension of the array. The imaged field-of view extends well beyond the transverse size of the array.

  fails. The increase in the slope provided by the aberration correction process indicates an improvement in the focusing operation and consequently an enhancement of the detec-time reversal analysis of the distortion matrix . . . . . 156 III.5 3D structure of the NAFZ . . . . . . . . . . . . . . . . . . . . . . . . 161 III.5.1 Previous geophysical studies on the NAFZ . . . . . . . . . . . 163 III.5.2 Scattering profiles . . . . . . . . . . . . . . . . . . . . . . . . 164 III.5.3 Crustal structure . . . . . . . . . . . . . . . . . . . . . . . . . 165 III.5.4 Lower crust and Moho . . . . . . . . . . . . . . . . . . . . . . 167 III.5.5 Sub-Moho structure . . . . . . . . . . . . . . . . . . . . . . . 169 III.5.6 Fault dipping and Shear zone . . . . . . . . . . . . . . . . . . . 170 III.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 III.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 III.8 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 A Global aberration correction . . . . . . . . . . . . . . . . . . . 175

  The dense array was installed temporarily between May 2012 and October 2013. It consists of 73 3-component broadband seismometers, 66 stations arranged along 11 east-west lines and 6 North-South lines forming a rectangular grid and covering an area of 35 km by 70 km with a nominal inter-stations spacing of ∼7km (Fig. III.1a). Seven additional stations were deployed east of the rectangular array in a semi-circle shape. In this region, the fault splits into two major strands: the northern (NNAFZ) and southern (SNAFZ) strands (Fig. III.1b). The northern strand,

  braseys, 2002). The fault delineates three tectonic blocks (Fig. III.1b): (i) the Istanbul Zone (IZ) situated North of the northern branch, (ii) the Sakarya zone (SZ) situated to the South of the southern branch and (iii) the Armutlu-Almacik crustal block (AA) located

  III.1) to compute the cross-correlation functions of horizontal EE component over the 18 months of recording period. First, the data were down-sampled at 25 Hz and corrected from instrument response. Then, the data were split into one-hour win-

Figure III. 1 :

 1 Figure III.1: Study region and location of DANA array (DANA, 2012). (a) Map of the study region and location of DANA array geophones (black triangle. The surface traces of the NAFZ are represented by the red lines (Emre et al., 2018). The blue lines indicate the locations of the cross sections represented in Figs. III.6, III.7 and III.8: Western profile (W), central profile (C), eastern profile (E), Istanbul zone profile (IZ), Armutlu-Almacik profile (AA) and Sakarya zone profile (SZ). (b) Geological map of the region (adjusted from Taylor et al. (2019) and Akbayram et al. (2016). ). The major geological blocks are represented: Istanbul zone (IZ) in the North, Armutlu-Almacik (AA) block in the center and Sakarya zone (SZ) in the South. The Adapazari and Pamukova basin location are indicated by AB and PB, respectively.

( 1 df

 1 Fig. III.2a). Under a matrix formalism, this can be done by the following matrix product:R rr (z, f ) = P † (z, f ) × K(f ) × P * (z, f ),(III.6)A broadband FR matrix R rr (z) is obtained by integrating R rr (z, f ) over the frequency bandwidth:R rr (z) = ˆf2 f P † (z, f ) × K(f ) × P * (z, f ),(III.7) with f 1 = 0.1 Hz and f 2 = 0.5 Hz. The symbol * stands for phase conjugate. This operation amounts to do an inverse Fourier transform at time t = 0 over the frequency band [0.1 0.5] Hz. Physically, it corresponds to a ballistic time gating that tends to select singly-scattered waves associated with a scattering event in the focal plane. Each element of R rr = [R(r out , r in , z)] contains the wave-field that would be recorded by a virtual geophone located at r out = (x out , y out , z) if a virtual source at r in = (x in , y in , z) emits a pulse of length δt = ∆f -1 at the central frequency f 0 , with ∆f = f 2 -f 1 = 0.4 Hz. An example of the FR matrix R rr is shown at depth z = 25 km in Fig. III.2c. R rr

Figure

  Fig. III.2d shows the resulting 2-D image I at z = 25 km retrieved from the diagonal of the FR matrix in Fig. III.2c. Such an image is an estimator of the reflectivity at depth z with a transverse resolution that would be ideally only limited by diffraction: δx = λ/(2 sin θ) ∼ 8.2 km (Born and Wolf, 2003) where θ = arctan(D/2z) is determined by the size of the array D = 50 km, and corresponds to the maximum angle under which a focusing point sees the geophones' array. By stacking the confocal image computed at

Figure

  Fig. III.2d.

Fig

  Fig. III.2e displays the averaged CMP intensity over the whole field-of-view. Even though the layered velocity model was optimized in order to obtain a good focusing, the lateral variation of velocity results in a spreading of energy over the anti-diagonal elements. In an ideal homogeneous medium, the CMP intensity profile would be depicted by a diffraction-limited PSF: I(r m , ∆r = H 2 0 * H 2 0 (∆r), with H 0 (∆r) = sinc (π∆r/δx), with δx the diffraction-limited focal spot. The averaged CMP intensity profile (Fig. III.2e) should be compared with the averaged CMP intensity profile obtained in an ideal layered medium. The latter is displayed in Fig. III.2f. The spatial extension of the ideal PSF is found to be of the order of δx. Hence, the resolution of the image will be compared with respect to the diffraction limited resolution at the considered depth δx that is represented by the white circle in Fig. III.2e and f. The CMP intensity profile shows an important spreading of energy outside the diffraction limited resolution which is a manifestation of aberrations. This observation originates from the fluctuations of the wave velocity distribution compared to our model in Table III.1. It results in a loss of resolution and contrast in the confocal image displayed in Fig. III.2d and f.

R

  kr (z) = T 0 × R rr (z), (III.15) where T 0 is the Fourier transform operator defined by Eq. III.4. R kr (z) = [R(k out , r in , z)] connects each input focusing point r in to the reflected wave-field in the far-field (Fig. III.4a).

Fig

  Fig. III.4a); (ii) a distorted component due to the mismatch between the velocity modeland the reality. The latter component refers to the phase distortions that should be isolated from the reflected wave-field in order to be properly compensated. This can be done by subtracting the ideal wave-front that would be obtained in absence of aberrations (i.e the geometric component) from each reflected wave-front induced by each input focusing wave at r in . Such operation can be expressed mathematically via a Hadamard product between R kr (z) and T 0 . It yields the distortion matrix D kr :

Figure

  Figure III.4: Local aberration correction: correction over two isoplanatic patches. (a)Far-field projection of the FR matrix in a heterogeneous medium. For nearby sources, associated with a same isoplanatic patch, the phase distortions of the reflected wave-field are identical. (b) For each focusing point, the reflected wave-field is stored along a line of R kr . (c) A set of local distortion matrices, each matrix being associated with an isoplanatic region (1 and 2), are obtained at each depth after removing the geometrical component. (d) The local distortion matrices are equivalent to a reflection matrix in which the far-field wave-fronts have been realigned with each other. (e) The singular value decomposition of D or of its normalized correlation matrix Ĉ yields an aberration phase law for each isoplanatic patch by synthesizing a virtual coherent reflector. (f) The aberration phase law obtained over each patch enables to estimate the transmission matrix T. The phase conjugate of T provides the focusing laws to compensate for the output and input phase distortions over each isoplanatic patch.

  number of input focusing points in each region should be at least four times larger than the number of resolution cells mapping the aberrated initial focal spot (Fig. III.2d) (Lambert et al., 2020b). The area covered by the initial focal spot being of 22 × 14 km 2 (Fig. III.2e), the extent of the window is chosen to be 55 × 55 km 2 .

Figure

  Figure III.5: Aberration correction process. Output (1) and input (2) phase distortion correction at z = 25 km. (a) Confocal image obtained after applying the aberration output ( Û * out ) (c1) and input ( Û * in ) (c2) phase laws corresponding to each region. The correlation coefficients between the different phase laws and the central one are displayed below each phase mask. (b) The averaged common midpoint intensity profile obtained after each aberration correction step. (3) Final step of the correction process: final confocal image (a3) and averaged common midpoint intensity profile (b3). The red lines represent the NAFZ fault traces at the surface.

  each iteration, the FR matrix is updated. The confocal image obtained after the additional two iteration steps and the averaged CMP intensity profile are displayed in Fig. III.4a3 and b3, respectively. The averaged CMP intensity profile (Fig. III.4b3) shows a nearly diffraction-limited focal spot with almost all the energy concentrated in the vicinity of the focusing point proving that most of the initial aberrations have been compensated (see the comparison with Fig. III.4b2). The residual background is probably associated with highorder aberrations whose coherence length (isoplanatic area) is smaller than the size of the window functions W (r). The gain in contrast at each step of the correction process can be quantified from the intensity ratio between the corrected and the initial image (Fig. III.2e).

  The efficiency of the local correction is illustrated in Appendix A by comparing it with a global correction, when no division of the field-of-view is done. The averagedCMP intensity profile inFig. A obtained with a global correction, displays an important residual background that covers an area of 35 × 20 km 2 while with the local correction, residual background in the averaged CMP intensity profile the focal spot spreads over an area of 20 × 20 km 2 (Fig. III.4b3).

Fig

  Fig. III.6a shows the same cross-section as in Fig. III.2g. The differences between the corrected image and the raw image are pronounced. While in the raw image, no clear structures and layers are visible, the corrected image reveals features delineating three blocks which were hidden in the original image.

Figure

  Figure III.6: Vertical North-South profiles of the 3D scattering volume oriented perpendicular to the fault traces (their locations are shown in Fig. III.1a). (a) Western profile at 30.2 • E, (b) eastern profile at 30.37 • E. The color scale refers to the scattering intensity. It is normalized by the maximum value of the scattering energy inside the volume. The locations of the southern (SNAF) and northern (NNAF), and the major crustal blocks (SZ: Sakarya zone, AA: Armutlu-Almacik and IZ: Istanbul zone) are labeled. The interpreted location of the fault at depth are indicated by red arrows. The interpreted crustal and subcrustal structure are numbered and indicated by black dashed lines. The position of the Moho is indicated by red dashed lines.

  tering distribution. From the surface traces displayed in Fig. III.1a, the northern and southern strands connect to the East of DANA array. Armutlu-Almacik (AA) block is then narrower in the East than in the West.

  Profiles have been extracted from the 3D volume and are plotted in Figs. III.6 and III.7, displaying two North-South and three East-West profiles, respectively. The locations of the cross sections are indicated in Fig III.1a. The North-South profiles cross the two fault strands, while the East-West profiles span the three geological blocks.

Figure

  Figure III.7: Vertical East-West profiles of the 3D scattering volume. The slices span the three crustal terranes. (a) Istanbul zone (IZ) profile at 40.83 • N, (b) Armutlu-Almacik (AA) profile at 40.58 • N and (c) Sakarya zone (SZ) profile at 40.37 • N. The locations of the profiles are shown in Fig. III.1a. The color scale refers to the scattering intensity. It is normalized by the maximum value of the scattering energy inside the volume. The interpreted crustal and sub-crustal structure are numbered and indicated by black dashed lines. The position of the Moho is shown by red dashed lines.

  Fig. III.7, it can be noticed that the crust of the Istanbul zone shows a strong heterogeneity.The crust beneath the AA block appears to be relatively less heterogeneous with more scattering to the West.In all the profiles, we observe a high scattering zone between 25 and 40 km depth, which we identify as the heterogeneous lower crust. The reflectivity of the lower crust varies between the North and the South, mainly below the two fault strands. This reflectivity is observed to be deeper in some regions compared to others, for example in the North (seeFig. III.6b). At greater depths (> 40 km), structures in the mantle are observed, mainly beneath IZ and AA block.

  same observation was made by[START_REF] Rost | Structure of the northwestern north anatolian fault zone imaged via teleseismic scattering tomography[END_REF]. A shallow structure is identified at 13 km around 40.6 • (marked by 2 in Fig. III.6b). Deeper in the crust, the reflectivity is localized near the northern branch with no clear indication of the presence of a southern strand in the crust. The most noticeable feature in the eastern profile is the highly reflective structure at depth between 18 and 25 km in the IZ (indicated by 3 in Fig. III.6b), dipping northward.

  Fig. III.6a. Similar to the western profile, sub-Moho structures are identified in the eastern profile beneath AA block (Fig. III.6b). A structure is detected between 45 and 50 km and another one between 55 and 60 km (indicated by 4 and 5 respectively). The latter coincides with the transition between AA and IZ blocks and can be interpreted as backscattered energy from a heterogeneous shear zone in the upper mantle or from the truncation of interfaces across the fault. Here also, there is a strong evidence that the NNAF extends to greater depths and is rooted in the upper mantle. Whereas the lack of significant scattered energy in SZ indicates that the SSAF does not reach the upper mantle. In the E-W profiles, sub-Moho structures are also detected. Below the IZ (Fig. III.7a), a reflective structure rises between 50 and 60 km (marked by 3) highlighting the complexity of the upper Mantle. Below AA block (Fig. III.7b), a reflective layer 45 and 50 km (marked as 4) indicates the presence of a sub-Moho structure appearing in the western, central and eastern parts of this profile. Even below the SZ profile (Fig. III.7c), a relatively weak sub-Moho reflection is observed between 50 and 54 km (marked as 1).

  It appears to dip northward (Fig. III.6a). Sub-Moho scattering structures in the proximity of the northern branch originate from the discontinuities of layers and geological features and from the damage associated with a localized shear in the mantle of the NNAF.The Moho continuity below the southern strand in the western and eastern profiles is supported by[START_REF] Jenkins | Crustal thickness variation across the sea of marmara region, nw turkey: A reflection of modern and ancient tectonic processes[END_REF]. S-wave tomography[START_REF] Papaleo | Constraints on north anatolian fault zone width in the crust and upper mantle from s wave teleseismic tomography[END_REF] shows a termination of the low velocity zone at Moho depth beneath the southern strand. Together these findings, with the evident signature of the Moho in our reflectivity profile in the South, indicate that the SNAF is confined in the crust and does no extend to the upper mantle.

Figure

  Figure III.8: Interpreted vertical North-South profiles of the 3D scattering volume oriented perpendicular to the fault traces (the locations of which are shown in Fig. III.1a. (a) Western profile at 30.2 • E, (b) Central profile at 30.3 • E and (c) Eastern profile at 30.37 • E. The color scale refers to the scattering intensity. It is normalized by the maximum value of the scattering energy inside the volume. The locations of the southern (SNAF) and northern (NNAF), and the major crustal blocks (SZ: Sakarya zone, AA: Armutlu-Almacik and IZ: Istanbul zone) and the upper (UC), middle (MC) and lower crust (LC) with the upper mantle (UM) are labeled. The interpreted location of the fault at depth are indicated by black thick lines. The shaded area delineated by black dashed lines correspond to the possible location of the shear zone beneath the northern strand.

  Figure A: Global aberration correction at z = 25 km. (a) The obtained confocal image and (b) the corresponding averaged common midpoint profile intensity. The red lines represent the NAFZ fault traces at the surface.

  Figure B: Vertical North-South profiles of the 3D scattering volume obtained after performing four iterations of a global aberration correction. The profiles are oriented perpendicular to the fault traces (the locations of which are shown in Fig. III.1a). (a) Western profile at 30.2 • E, (b) eastern profile at 30.37 • E. The color scale refers to the scattering intensity. It is normalized by the maximum value of the scattering energy inside the volume. The locations of the southern (SNAF) and northern (NNAF), and the major crustal blocks (SZ: Sakarya zone, AA: Armutlu-Almacik and IZ: Istanbul zone) are labeled.

1

 1 Fig. B two scattering profiles are plotted. These profiles are chosen at the same location as the profiles in Fig. III.6. The differences between the two corrections are noticeable. In Fig. B, only few structures are observed and the discontinuities below the fault strands are less apparent in comparison with Fig. III.6. An explanation of this result is that the eigenvector of the global distortion matrix tends to only grasp the low order aberrations that are invariant across the field-of-view and cannot capture the higher-order aberrations associated with smaller isoplanatic area. We conclude that, in the present case (random and continuous distribution of reflectivity), a local aberration correction leads to a better estimation of the transmission matrix, in particular its high order aberration components that rapidly fluctuate across the fieldof-view. Focusing laws are therefore optimized for each point in the medium, leading a far better resolution and contrast for the final image of NAFZ. The time focused reflection matrix . . . . . . . . . . . . . . . . . . . 179 IV.2 Propagation movies . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 IV.2.1 Correct wave velocity model . . . . . . . . . . . . . . . . . . . 181 IV.2.2 Impact of a wrong velocity model on the propagation movie . . 183 IV.3 Focusing map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 IV.4 Extension to a multi-layered model . . . . . . . . . . . . . . . . . . . 188 IV.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189

Figure IV. 1 :

 1 Figure IV.1: Propagation movie obtained from the time focused reflection matrix. (a) Vertical profile of the original 3D image below DANA network at 30.36 • E with the location of the input focusing point r in (white cross). (b)-(f) Propagation movie from R ∆r,r (δt) matrix at different times. (1,2) The amplitude and the real part of the matrix, respectively. The red curve indicated the location of the theoretical coherent wave.

  Figs. IV.1(b)-(f) shows the amplitude and the real part of this propagation movie. Note

Figure

  Figure IV.3: (a,b,c) Propagation movie at δt = 0 associated with an echogenic reflector for a two travel time t = 27.8 s, calculated for three different homogeneous velocity model c 0 = 3300, 3500 and 3700 m/s, respectively. The dashed black curve illustrates the convergence and divergence of the wave. (d,f) Propagation movie at a specific focusing time δt for c 0 = 3300 and 3700 m/s, respectively (1,2) Amplitude and real part of the propagation movie, respectively. (e) Vertical profile of the original 3D image below DANA network at 30.36 • E with the location of the input focusing point r in .

  shift ∆z in of the focal spot in Fig. IV.3d, the integrated velocity c is estimated to be ∼ 3500 m/s, which confirms that the homogeneous wave velocity c 0 = 3500 m/s originally considered in Fig. IV.1 was a good estimation of the integrated velocity for the considered travel time. We also observe that for an overestimated velocity, the focusing occurs at smaller depth than r in and at negative times (Fig. IV.3f). The wave is already diverging at δt = 0 (black dashed curve in Fig. IV.3c).

FigureFigure

  Figure IV.4: Focusing maps. (a,b,c) Amplitude of the focusing maps (Eq. IV.5) associated with a reflective structure represented in Fig. IV.3e, calculated for three different homogeneous velocity models c 0 = 3300, 3500 and 3700 m/s, respectively. (b1) Real part of the focusing map F computed for c 0 = 3500 m/s. (d) Normalized amplitude of the three focusing maps F along the axial axis ∆y = 0.

Fig

  Fig. IV.4 displays the amplitude of the focal spots for the three velocity models considered earlier. Its real part is also shown for c 0 = 3500 m/s. For the case of a correct velocity model (Fig. IV.4b), the amplitude reveals a focal spot that is centered around ∆r = 0, while the real part shows an axial phase shift at ∆r = 0. The phase shift originates from the Gouy phase shift that undergoes a focused wave-field when it goes through the focusing point (Feng and Winful, 2001). In 3D, this Gouy phase shift of π is actually what we observe in Fig. IV.4. For an incorrect velocity, the spot is shifted downward (Fig. IV.4a) or upward (Fig. IV.4c), whether the velocity is underestimated or over estimated, respectively. The maximum of the spot indicates where the waves focus and allows to estimate the shift ∆z. Fig. IV.4d

Figure

  Figure IV.6: Focusing map obtained for a multi-layered velocity model. (a,b) Amplitude of the focusing maps R w associated with a reflective structure represented in Fig. IV.3e, calculated for layered velocity model considered to build images in chapter III. (b) Normalized amplitude of the focusing maps along the z-axis axis for ∆x = ∆y = 0.

1 .

 1 The focusing point is chosen for the same two way travel time t = 27.8 s at 30.36 • E, as in Sec. IV.2.1 (Fig. IV.4a). The focusing map associated with the virtual source is displayed in Fig. IV.6a and the normalized amplitude plot along z for ∆x = ∆y = 0 is shown in Fig. IV.6b. The focal spot appears to be slightly shifted downward. This indicates that the integrated velocity predicted by the considered model is slightly smaller than the correct one. To confirm this observation, the integrated velocity of the layered model can be calculated from Eq. IV.3: c (z in ) = 3275 m/s which is indeed smaller than 3475 m/s. This observation shows that the layered velocity model is not perfect. Lateral variations around the fault zone between the true velocity and the layered velocity model are expected and have already been revealed by the transverse aberrations highlighted in chapter III.

tolian

  Fault. We used ambient noise measurements recorded at 73 stations and covering an area of 35 km by 70 km. EE cross-correlations are computed in [0.1 0.5] Hz frequency range. A layered S-waves velocity model is considered. Contrary to SJFZ, the medium here displays both specular reflectors and a random distribution of scatterers. Matrix imaging is adapted to this scattering context and a local analysis of the distortion matrix is performed. It leads to an estimation of the transmission matrix between the geophones and the set of focusing points that map the underground. The knowledge of this matrix is a significant breakthrough in seismic imaging since it provides optimized focusing laws for each point of the medium. A diffraction limited image of NAFZ is obtained until a depth of 60 km at a transverse resolution of 2 km. It reveals the structure of the crust and the upper mantle below NAFZ. The crust thickness is reported at several locations and a step of the Moho is found in the North. We also find strong evidence of a localized shear zone that penetrates the crust and reaches the upper mantle below the northern branch.The structure of the major geological blocks in the region is also discussed. localized in Sapanca lake region. The highest fractal dimension values are also observed in this region. Looking at the scattering distribution, we notice that the seismicity is distributed above or on the boarders of the regions with highest reflectivity (Fig. V.1B1

  Figure V.1: A) Map view of the earthquakes' template with color coded fractal dimension. The red box corresponds to the extension of the scattering volume shown in chapter III. B) Depth distribution of the events along parallel and perpendicular cross-section (boxes in A) with color coded fractal dimension. The events are projected onto the boxes' central axis. The scattering profiles at the center of the boxes are represented. Courtesy of Eric Beauce.

  : (i) focusing in emission to synthesize a focused virtual source at r in (x in , y in , z in ) inside the medium (Fig. I.7a), and (ii) focusing in reception to synthesize the response of a virtual receiver at r out (x out , y out , z out ) inside the medium (Fig. I.7b

  . . . . . . . 103 II.2 Distribution of seismic scatterers in the San Jacinto Fault Zone, southeast of Anza, California, based on passive matrix imaging . . . 106 II.2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 II.2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 II.2.3 The San Jacinto Fault Zone . . . . . . . . . . . . . . . . . . . . 110

II.2.4 3D scattering volume . . . . . . . . . . . . . . . . . . . . . . . 112 II.2.4.1 Images of the shallow fault zone . . . . . . . . . . . 114 II.2.4.2 Deeper and larger scale structure . . . . . . . . . . . 116 II.2.5 Lateral variations of intensity . . . . . . . . . . . . . . . . . . 119 II.2.6 Discussion and conclusions . . . . . . . . . . . . . . . . . . . 122 II.2.7 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 A Imaging and aberration correction . . . . . . . . . . . 125 B Specular vs. diffuse reflection . . . . . . . . . . . . . 128 C Time-gain compensation . . . . . . . . . . . . . . . 130 D Supplementary Data . . . . . . . . . . . . . . . . . . 131

  1,2 , Thibaud Blondel 2 , Arnaud Derode 2 , Michel Campillo 1 ,

	Alexandre Aubry 2
	1 ISTerre, Université Grenoble Alpes, Maison des Géosciences, BP 53, F-38041 Greno-
	ble, France.

2 Institut Langevin, ESPCI Paris, PSL University, CNRS, Univ Paris Diderot, Sorbonne Paris Cité, 1 rue Jussieu, F-75005 Paris, France. Published in Geophysical Journal International 2021 (doi: 10.1093/gji/ggab133)

  .2.

	Correction steps Contrast (dB)
	1	5
	2	7.71
	3	7.75
	4	8.2
	Table III.2: Result of the correction process.

Acknowledgements

Acknowledgments The authors are grateful for funding provided by LABEX WIFI (Laboratory of Excellence within the French Program Investments for the Future, ANR-10-LABX-24 and ANR-10-IDEX-0001-02 and by TOTAL R&D. We acknowledge the support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No 742335, F-IMAGE and grant agreement No. 819261, REMINISCENCE). Acknowledgments We acknowledge support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No 742335, F-IMAGE and grant agreement No. 819261, REMINISCENCE). YBZ acknowledges support from the Department of Energy (Award DE-SC0016520). The paper benefited from constructive comments by two anonymous referees and Editor Rebecca Bendick. Acknowledgments This project has received funding from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No

DATA AVAILABILITY

The seismic data used in this study can be obtained from the Data Management Center of the Incorporated Research Institutions for Seismology.

II.2.1 Abstract

Fault zones are associated with multi-scale heterogeneities of rock properties. Large scale variations may be imaged with conventional seismic reflection methods that detect offsets in geological units, and tomographic techniques that provide average seismic velocities in resolved volumes. However, characterizing elementary localized inhomogeneities of fault zones, such as cracks and fractures, constitutes a challenge for conventional techniques. Resolving these small-scale heterogeneities can provide detailed information for structural and mechanical models of fault zones. Recently, the reflection matrix approach utilizing body wave reflections in ambient noise cross-correlations was extended with the introduction of aberration corrections to handle the actual lateral velocity variations in the fault zone [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF]. Here this method is applied further to analyze the distribution of scatterers in the first few kilometers of the crust in the San Jacinto Fault Zone at the Sage Brush Flat (SGB) site, southeast of Anza, California. The matrix approach allows us to image not only specular reflectors but also to resolve the presence, location

DATA AVAILABILITY

The seismic data used in this study can be obtained from the [START_REF] Vernon | Sage brush flats nodal experiment. international federation of digital seismograph networks[END_REF]

II.2.7 Appendices

A Imaging and aberration correction

In this section, we review the reflection matrix technique [START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-d imaging: application to the san jacinto fault zone, california[END_REF] applied to obtain the 3D scattering images of the accompanying paper. First we describe the virtual focusing method inside the medium and then we recall the aberration correction process. Ambient noise cross-correlations between the stations are computed in the 10-20 Hz range. These correlations provide the impulse response between virtual sources (emission) and receivers (reception) at the surface. In other words, the cross-correlations form a response matrix K(t) of the underground associated with the geophones' array. This matrix contains considerable of information about the medium, especially the distribution of its local reflectivity. This information is retrieved from the response matrix by performing focusing operations, both in emission and reception. Focusing consists in applying appropriate time delays, at emission and reception, so that the seismic waves resulting from a scattering event at the focal point interfere constructively. The goal of this operation is to virtually move the sources and receivers inside the medium onto a virtual plane. It is commonly referred to as redatuming [START_REF] Berkhout | A unified approach to acoustical reflection imaging. II: The inverse problem[END_REF] and can be easily implemented in the frequency domain via simple matrix products. The first step is to apply a temporal Fourier transform to the correlation matrix to obtain K(f ).

To perform focusing, we define at each depth a basis of focal points r, which corre-

D Supplementary Data

Supplementary data associated with this article can be found in the online version. the existence of specular reflectors that prevent from a clear discrimination between input and output aberrations [START_REF] Barolle | Approche matricielle de la tomographie à cohérence optique[END_REF]; (ii) the limited size of spatial windows that does not allow a full convergence of the SVD towards the aberration phase law. The latter issue could be solved in principle by starting the aberration correction process with larger spatial windows and gradually reducing their size at each iteration of the time reversal process (Lambert et al., 2021b). Further works will be done to clarify these imperfections of the iterative time reversal process.

Nevertheless, the input phase laws can be used to estimate the input transmission matrix Tin = Ûin • T 0 . The FR matrix can then be updated:

A corrected confocal image is deduced from the diagonal elements of R

III.8 Appendix A Global aberration correction

In this appendix, we validate our choice of choosing a local aberration correction for the case of the NAFZ. To do so, we show the result of a global aberration correction in which

General conclusion and perspectives

General conclusion

In this thesis, I presented a passive matrix approach for imaging in complex media, such as fault zones. This method addresses the principle issue of imaging in heterogeneous media:

the seismic velocity variations inside the Earth. It uses ambient noise cross-correlations to reveal the scattering structure around two major fault zones, at different scales: the San Jacinto fault and the North Anatolian fault.

In the last two decades, the use of ambient noise vibrations have seen a major increase in geophysical applications. The cross-correlations of noise records were shown to converge towards the Green's function of the medium. Because the ambient noise is omnipresent, cross-correlations method constitutes an alternative low-cost method to using active sources. Body waves have been retrieved from the cross-correlations at several scales, leading to numerous applications for deep Earth's imaging. Inspired by previous works in acoustics, optics and seismology [START_REF] Prada | Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media[END_REF][START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF][START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF]Aubry and Derode, 2009a,b;[START_REF] Badon | Approche matricielle de l'imagerie optique des milieux diffusants[END_REF][START_REF] Barolle | Approche matricielle de la tomographie à cohérence optique[END_REF][START_REF] Blondel | Matrix approach of seismic imaging[END_REF]Lambert et al., 2020a,b), we developed a passive imaging technique that provides unprecedented insight on the scatterers distribution in the subsurface.

The idea is to retrieve from noise cross-correlations all the available information on the medium's structure. Conventional migration techniques require an accurate velocity model of the subsurface, limiting their application to regions where the velocity is perfectly known. The method presented in this thesis can be applied in any medium with little knowledge on the velocity distribution. First, a back-propagation of cross-correlations of Finally in the last chapter (Chap. IV), we proposed a route towards a mapping of the 3D wave velocity distribution, thereby going beyond a simple reflectivity image. To that aim, a temporal analysis of the reflection matrix is done in the virtual source frame where we study the time response recorded at virtual receivers when a virtual source is synthesized inside the medium. The propagation of the coherent wave is studied via the focusing map. This analysis provides information of the integrated velocity inside the medium and paves the way towards an estimation of the local body wave velocity.

To summarize, the results presented in this thesis provide unprecedented view of two major fault zones contributing to our understanding of the structure and dynamics of fault zones. Beyond the specific case of fault zones, this matrix approach could be used to image volcanoes. It can also be applied in marine explorations. Due to its flexibility and low computational cost, the presented matrix approaches can be revolutionary in geophysical investigations.

Perspectives

The results presented in this thesis are very promising. Of course, these results should be compared with future studies conducted on the two fault zones. A short term future perspective of this work is to combine the scattering images around the NAFZ with a new earthquake catalog built by Beauce and colleagues (in prep.).