A Distortion Matrix Framework for High-Resolution Passive Seismic 3-D Imaging: Application to the San Jacinto Fault Zone, California
Rita Touma, Thibaud Blondel, Arnaud Derode, Michel Campillo, Alexandre Aubry
AA Distortion Matrix Framework for High-Resolution Passive Seismic 3D Imaging: Application tothe San Jacinto Fault Zone, California
Rita Touma,
1, 2
Thibaud Blondel, Arnaud Derode, Michel Campillo, and Alexandre Aubry ∗ ISTerre, Universit´e Grenoble Alpes, Maison des G´eosciences, BP 53, F-38041 Grenoble, France Institut Langevin, ESPCI Paris, PSL University, CNRS, Univ Paris Diderot,Sorbonne Paris Cit´e, 1 rue Jussieu, F-75005 Paris, France (Dated: August 5, 2020)Reflection seismic imaging usually suffers from a loss of resolution and contrast because of the fluctuationsof the wave velocities in the Earth’s crust. In the literature, phase distortion issues are generally circumventedby means of a background wave velocity model. However, it requires a prior tomography of the wave velocitydistribution in the medium, which is often not possible, especially in depth. In this paper, a matrix approach ofseismic imaging is developed to retrieve a three-dimensional image of the subsoil, despite a rough knowledge ofthe background wave velocity. To do so, passive noise cross-correlations between geophones of a seismic arrayare investigated under a matrix formalism. They form a reflection matrix that contains all the information avail-able on the medium. A set of matrix operations can then be applied in order to extract the relevant information asa function of the problem considered. On the one hand, the background seismic wave velocity can be estimatedand its fluctuations quantified by projecting the reflection matrix in a focused basis. It consists in investigatingthe response between virtual sources and detectors synthesized at any point in the medium. The minimizationof their cross-talk can then be used as a guide star for approaching the actual wave velocity distribution. Onthe other hand, the detrimental effect of wave velocity fluctuations on imaging is overcome by introducing anovel mathematical object: The distortion matrix. This operator essentially connects any virtual source insidethe medium with the distortion that a wavefront, emitted from that point, experiences due to heterogeneities. Atime 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 forany point of the underground. Applied to seismic data recorded along the Clark branch of the San Jacinto faultzone, the present method is shown to provide an image of the fault until a depth of 4 km over the frequencyrange 10-20 Hz with a transverse resolution of 80 m. Strikingly, this resolution is almost one eighth below thediffraction limit imposed by the geophone array aperture. The heterogeneities of the subsoil play the role ofa scattering lens and of a transverse wave guide which increase drastically the array aperture. The contrast isalso optimized since most of the incoherent noise is eliminated by the iterative time reversal process. Beyondthe specific case of the San Jacinto Fault Zone, the reported approach can be applied to any scales and areas forwhich a reflection matrix is available at a spatial sampling satisfying the Nyquist criterion.
Waves constitute a powerful means to non destructivelyprobe an unknown medium. Indeed, wave propagation is fullydetermined by the wave equation and boundary conditions.As stated by diffraction theory in the acoustic approximation,knowing the incident wave-field and the internal propertiesin terms of density ρ and celerity c theoretically allows tocompute the wave field everywhere and at any time inside themedium. This is the so-called ”forward” or ”modeling” prob-lem. Conversely, when the medium is unknown, the lack ofcelerity and density knowledge makes it impossible to com-pute the spatio-temporal evolution of the wave field. Yet thisevolution can be known, at least at the boundary, through ex-perimental measurement of the wave-field scattered by themedium. The inverse problem then consists in deducing themedium internal properties from the recording of the wavefield at its surface. A first way to do so is to assume a ve-locity and density background, solve the forward problem tocompute the time-dependent signal that would be backscat-tered if the background model was true, and iteratively updatethis model to minimize the difference with the actual record-ings. Another way is to directly back-propagate the scattered ∗ [email protected] echoes to reflectors inside the medium. This also amountsto updating a background model since a reflector is nothingelse than a variation in acoustic impedance ρc . In both strate-gies, referred to as ”inversion” and ”migration”, respectively,a celerity macro model is required and the purpose is to com-pute variations from this model under the assumption that theyare small (Born approximation). If they are not, the reflectedwave-field may be subject to aberrations and multiple scatter-ing that the macro model fails at modeling. These issues leadto distorted images, lack of resolution and unphysical features,which are very detrimental to the imaging process.In seismic exploration, these issues are important becausein most cases the celerity is non constant in space and its dis-tribution is unknown. Most geological settings actually con-sist of several layers of rocks and sediments with distinct me-chanical properties as well as location-dependent thickness,faulting and strata organization. These may be difficult toestimate without previous geologic expertise of the subsur-face, especially in areas with high lateral mechanical stressthat bend and break the layers and make them superimpose.When trying to retrieve details at the physical limit, the pre-vious knowledge required to build reliable images may be al-ready fairly demanding. On the one hand, the inversion prob-lem cannot be solved if the initial celerity model is too farfrom the reality because the regularization procedure can end a r X i v : . [ phy s i c s . g e o - ph ] A ug up stuck in a local minimum. On the other hand, migrationtechniques would lead to loss of resolution due to phase dis-tortions and a blurred image due to multiple scattering. Thatbeing said, the question that naturally arises and which thepresent work aims at addressing is: how to retrieve an accu-rate image when little to no previous knowledge on the spatialvariations of the wave speed is available?To cope with this issue, our strategy is to develop a ma-trix approach of seismic wave imaging. In a linear scatter-ing medium, a reflection matrix relates the input and outputon a single side of the medium. It contains all the relevantinformation on the medium as it fully describes wave propa-gation inside the scattering medium. In the last decade, theadvent of multi-element arrays with controllable emitters andreceivers has opened up a route towards the ability of mea-suring a reflection matrix in the case where the input andoutputs points are located on the same side of the scatter-ing medium. In particular, the reflection matrix has beenshown to be of great interest for detection and imaging pur-poses in scattering media, whether it be in acoustics [1–3] oroptics [4, 5]. The reflection matrix contains the set of inter-element impulse responses recorded between each array el-ement. It has already been shown to be a powerful tool forfocusing in multi-target media [6, 7], as well as for separat-ing single and multiple scattering [5, 8] in strongly scatteringmedia. These matrix methods have been successfully appliedto geophysics in the extremely challenging case of the Ere-bus volcano in Antarctica [9]. However, although the maininternal structures are successfully revealed by removing themultiple scattering background, the matrix image still suffersfrom the phase distortions induced by the long-scale fluctua-tions of the seismic velocity. To retrieve a diffraction-limitedresolution, the concept of distortion matrix has been recentlyintroduced by two seminal works in ultrasound imaging [10]and in optical microscopy [11]. Inspired by the pioneeringwork of Robert and Fink [1], the distortion matrix D is definedbetween a set of incident plane waves [12] and a set of pointsinside the medium [1]. It contains the deviations from an idealreflected wavefront which would be obtained in the absence ofinhomogeneities. As shown by recent studies [10, 11], a timereversal analysis of the D -matrix allows to synthesize virtualreflectors in depth. This process can then be leveraged forunscrambling the phase distortions undergone by the incidentand reflected wavefronts.Overcoming such phase distortions would be especiallyvaluable for geophysical applications given the stratifiedstructures of the environments of interest. Migration tech-niques in Fourier domain have actually been very popular forimaging in layered media [13, 14], however they only holdfor 1D celerity models with no lateral variations. Subsequentworks have focused on adapting these techniques to take intoaccount increasing lateral velocity variations, at the cost ofmore numerical and computational complexity [15–17]. Con-trary to these well-established methods, the matrix approachdoes not require any assumption on the structures and on thevelocity distribution inside the medium, while being fairlylight on the computational aspect. The present paper aims atstudying the relevance of the matrix approach for geophysical imaging.Coherent sources (vibrating trucks, explosives, etc. ) canbe used in shallow subsurface ( < et al. [32] showedthe possibility of mapping the upper mantle discontinuities (at410 and 660 km of depth) by extracting body waves reflec-tion from ambient noise, while Retailleau et al. [33] mappeda region of the core-mantle boundary at about 2900 km depth.In this paper, inspired by the reflection matrix approachdeveloped by Blondel et al. [9] and based on noise cross-correlations, the distortion matrix approach is extended to sat-isfy seismic imaging purposes. The method is applied to SanJacinto fault zone (SJFZ) site. Fault zones are indeed amongthe most challenging media for seismic imaging given theirhighly localized and abrupt variations of mechanical proper-ties, extensive fractures and damage zones. In that respect,the SJFZ is the most seismically active fault zone in South-ern California [34]. It accounts for a large portion of the platemotion in the region [35, 36]. A highly complex fault-zonestructure with prominent lateral and vertical heterogeneities atvarious scales have already been highlighted in previous stud-ies [37–39]. In particular, maps of the P and S wave velocities, V P and V S , have been inverted from earthquake arrival timesfor a depth range of 2-20 km [37, 40]. Surface wave tomo-graphic images built from noise correlations revealed the ve-locity structure in the top 7 km of the complex plate boundaryregion at a resolution of about ten kilometers [38]. To com-plement these regional studies and provide structural featuresin the first few kilometers with an improved resolution, ambi-ent noise at higher frequency up to 10 Hz was analysed fromdata recorded by a dense rectangular array deployed aroundthe Clark branch of the SJFZ [39, 41, 42]. In particular,Zigone et al. [43] used ambient noise cross-correlations in the2-35 Hz frequency range to derive a velocity model in the top100 m with a resolution of 50 m.Imaging deeper the fault area at such resolution is challeng-ing because of the damage and the complex distribution ofsmall-scale heterogeneities. Yet, a much larger penetrationdepth can be expected by taking advantage of the reflectedbulk waves. To do so, the matrix approach of seismic imag-ing is particularly useful since it only requires a rough ideaof the mean wave velocity. Besides, it shall provide a three-dimensional image of the subsoil acoustic impedance insteadof just the wave velocity. To implement this matrix approach,we take advantage of a spatially dense array of geophonesdeployed over the damage zone of SJFZ [41]. Noise cross-correlations are used to retrieve the impulse responses be-tween the geophones. The associated passive reflection ma-trix is then investigated to image the first few kilometers ofthe crust by virtue of body waves emerging from noise corre-lations. As a whole, the process we present in this paper canbe analyzed as a combination of six building blocks:• (B1) A Fourier transform of the recorded signals yieldsa set of response matrices K ( f ) associated with thedense array of geophones.• (B2) Based on a rough estimate of velocity c , a doublefocusing operation is performed both at emission andreception by means of simple matrix operations. A setof focused reflection matrices R ( f, z ) are obtained atany arbitrary depth z below the surface.• (B3) A coherent sum of these matrices over the fre-quency bandwidth yields a broadband reflection matrix R ( z ) at any depth z .• (B4) By projecting the input or output entries of thismatrix in the far-field, the distorted component D ofthe reflected wave-field can be extracted.• (B5) A virtual iterative time reversal process is appliedto the matrix D to extract the phase distortions un-dergone by the incident or reflected wave-fields duringtheir travel from the Earth surface to the focal plane.• (B6) The whole process converge towards the focusinglaws that shall be applied at input or output of the re-flection matrix in order to compensate for aberrations.As a result of these six steps, an in-depth confocal image ofthe SJFZ is built. While conventional migration methods leadto a badly resolved image of the SJFZ subsoil, the matrix ap-proach clearly reveals sedimentary layers close to the surface( z < m) and several geological layers at larger depth( m < z < m). The layers structure is shown to bedifferent on each side of the fault. Large dip angles are alsohighlighted in the vicinity of the fault. A structural interpreta-tion of the obtained images can be finally built on the existingliterature about SJFZ. I. REFLECTION MATRIXA. Response matrix between geophones
The data used in this study has been measured from May 7,2014 to June 13, 2014 by a spatially dense Nodal array con-sisting of 1108 vertical geophones straddling the Clark Branchof SJFZ, southeast of Anza [41]. Figure 1a shows the locationof the 1108 vertical geophones organized as a 600 m ×
700 mgrid with inter-station distances δu x ∼
10m and δu y ∼ − , from which cross-correlation has been per-formed after whitening in the 10-20 Hz range with time lags ranging from − s to +5 s. This provides an estimate of theimpulse response between every pair of geophones. Each geo-phone is denoted by an index i and its position s i . The impulseresponse between stations i and j is noted k ij ( t ) , with t thetime lag. The set of impulse responses form a time-dependentresponse matrix K ( t ) .Given the high density of the network, neighbouring geo-phones belong to the same coherence area of seismic noise.The characteristic dimension of this area is indeed of λ/ ∼
50 m which is larger than the interstation distance δu . This isresponsible for a strong auto-correlation signal around t = 0 for geophones located in the same coherence area. This peakis proportional to the seismic noise power and does not ac-count 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 theweight of the corresponding impulse responses k ij ( t ) whoseassociated geophones i and j are contained in the same coher-ence area (see Section S1 in supporting information).The impulse responses exhibit several direct arrivals thathave already been investigated by Ben-Zion et al. [41] andRoux et al. [39]. Ballistic waves, likely direct inter-station S-wave and P-wave, arrive before the Rayleigh wave at apparentvelocities larger than 1000 m/s. Roux et al. [39] used iterativedouble beamforming to map the phase and group velocitiesof Rayleigh waves across the fault in the 1-5 Hz frequencybandwidth. Subsequently, Mordret et al. [42] inverted thesedispersion curves to build a 3-D shear wave velocity modelaround the Clark fault down to 500 m depth. Assuming the V p /V s ratio to be a linear function of depth, the following av-eraged value was found for the P-wave velocity over the top800 m: V p ∼ V p inthe SJFZ region is not available beyond this shallow layer. Asa consequence, in the present study, we will use an approx-imated homogeneous P-wave velocity model of c = 1500 m/s. This choice will be validated and discussed a poste-riori by a minimization of the aberration effects in the 3Dimage (see Supplementary Fig. S1). We are not interestedin the ballistic component of the wave-field but rather in itsscattered contribution due to reflections by the in-depth struc-ture along the fault. These vertical echoes are mainly asso-ciated with P-waves since only the vertical component of theGreen’s functions is considered in this study. Unlike our pre-vious study on the Erebus volcano [9], the scattered wave-field consists of a single scattering contribution which is apriori largely predominant compared to the multiple scat-tering background. This will be confirmed a posteriori bythe reflection matrix features (see Sec. I B). Singly-scatteredechoes can then be taken advantage of to build a 3D imageof the subsoil reflectivity. This local information can be re-trieved from K ( t ) by applying appropriate time delays to per-form focusing in post-processing, both in emission and recep- a b r in r out La t i t ude -116.596 -116.594 -116.592 -116.59 -116.588 -116.586Longitude m m FIG. 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 SanJacinto Fault (Southern California). Each row along the x-direction is composed of ∼
55 sensors with a pitch of 10 m, and the nominalseparation 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 t/ ) yields the impulse response betweenvirtual geophones placed at these two points. The same operation is repeated for any couple of points in the focal plane and yields the focusedreflection matrix R . tion. While focusing in emission consists in applying propertime delays in the recorded seismic data so that they con-structively interfere at an arbitrary position at depth, focus-ing in reception consists in applying proper time delays in therecorded seismic data so that the information coming from anarbitrary position at depth constructively interfere. Based onthe Kirchoff-Helmholtz integral, such a focusing operation isstandard in exploration seismology and referred to as reda-tuming [47–49]. However, in the present case, the stronglyheterogeneous distribution of the seismic wave velocities in-duces strong phase distortions that degrade this imaging pro-cess. A prior quantification and correction of these phase dis-tortions is thus required to reach a diffraction-limited resolu-tion and an optimized contrast for the image. As we will see,a matrix formalism is a well-matched tool to locally capturesuch information. B. Focused reflection matrix
The reflection matrix can be defined in general as an en-semble of responses, each response linking one vector to an-other vector. The type of vector coordinates will be referredto as bases . They can be spatial coordinates (hence the vec-tor refers to an actual point within or at the surface of themedium, see Fig. 1b) or wave vector coordinates. Variousbases are involved in this work : ( i ) the recording basis ( u ),whose elements are the positions of the geophones., ( ii ) thefocused basis ( r ) which corresponds to the positions of vir-tual geophones at which focusing at emission or reception isintended; and ( iii ) the Fourier basis ( k ). Because of linearityand time-invariance, seismic data can be projected from therecording basis to the focused basis by a simple matrix prod- uct. In the frequency domain, simple matrix products allowseismic data to be easily projected from the recording basisto the focused basis where local information on the mediumproperties can be extracted [5, 9, 50].Consequently, we first apply a temporal Fourier transformto the response matrix to obtain a set of monochromatic ma-trices K ( f ) . To project K ( f ) into the focused basis, we thendefine a free-space Green’s matrix, G ( f ) , which describesthe propagation of waves between the geophones and focusedbasis. Its elements correspond to the 3D Green’s functionwhich connects the geophone’s transverse position u to anyfocal point defined by its transverse position r and depth z ina supposed homogeneous medium: G ( r , u , z, f ) = e − j πf √ (cid:107) r − u (cid:107) + z /c π (cid:113) (cid:107) r − u (cid:107) + z (1) K ( f ) can now be projected both in emission and receptionto the focused basis via the following matrix product at eachdepth z [9, 50]: R ( z, f ) = G ∗ ( z, f ) × K ( f ) × G † ( z, f ) , (2)where the symbols ∗ , † and × stands for phase conjugate,transpose conjugate and matrix product, respectively. Equa-tion (2) simulates focused beamforming in post-processing inboth emission and reception. Each coefficient of the focusedreflection matrix R ( z, f ) involves pairs of virtual geophones, r in = ( x in , y in ) and r out = ( x out , y out ) , which are locatedat the same depth z (Fig. 1b). For broadband signals, bal-listic time-gating can be performed to select only the echoesarriving at the ballistic time t B in the focused basis [50]: t B = ( (cid:107) r in − u in (cid:107) + (cid:107) r out − u out (cid:107) ) /c . Under a matrix for-malism, this time-gating can be performed by means of a co-herent sum of R ( f ) over the frequency bandwidth ∆ f . Ityields the broadband focused reflection matrix R ( z ) = (cid:90) f + f − df R ( z, f ) , (3)where f ± = f ± ∆ f / and f is the central frequency.Each element of R ( z ) contains the complex amplitude ofthe wave that would be detected by a virtual detector lo-cated 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 − at the central frequency f . Importantly, the broadband fo-cused reflection matrix synthesizes the responses between vir-tual geophones which have a greatly reduced axial dimension δz = cδt compared to their stretching δz = 2 λ/ sin θ inthe monochromatic regime [51]. θ = arctan( (cid:68) / z ) is themaximum angle under which the geophone array is seen fromthe common mid-point and (cid:68) ∼ m, the characteristicsize of the geophone array. As a consequence, considering abroadband reflection matrix R ( z ) will significantly improvethe accuracy and axial resolution of the subsequent analysis.For the sake of a lighter notation, we will omit, in the follow-ing, the dependence in z but keep in mind that the focusedreflection matrix differs at each depth.Figure 2a displays one example of the broadband focusedreflection matrix R at depth z = z = 3600 m (Fig. 2a); this is very differentfrom the Erebus volcano for which the reflection matrix dis-played a fully random feature [9]. Here a significant part ofthe backscattered energy is still concentrated in the vicinity ofthe diagonal of the focused reflection matrix. This indicatesthat single scattering dominates at this depth: The beam is fo-cused, scattered just once, and focused in reception. On thecontrary, a broadening of the back-scattered energy outsidethe diagonal would mean that the beam undergoes aberrationand/or multiple scattering. In fact, the diagonal elements of R ( r in = r out ) correspond to what would be obtained fromconfocal imaging: transmit and receive focusing are simulta-neously performed on each point in the medium. A confocalimage can thus be obtained from the diagonal elements of R ,computed at each depth: (cid:73) ( r , z ) ≡ | R ( r , r , z ) | . (4)Figure 2c displays the 2D confocal image built from the di-agonal of the reflection matrix in Fig. 2a at time t = 4 . s,hence at an effective depth z = c t/ m. Some scat-tering structures seem to arise at different locations along thefault but confocal imaging is extremely sensitive to aberrationissues. One thus has to be very careful about the interpreta-tion of a raw confocal image. This observation is confirmedby Figure 2d that displays a slice of the confocal B-scan ofthe SJFZ underground. Each speckle grain in this image oc-cupies a major part of the field-of-view. Hence, aberrations seem to be pretty intense at large depths (beyond 1500 m) andthe inner structure of the SJFZ cannot be deduced from a basicconfocal image.Fortunately, the matrix R contains much more informationthan a single confocal image. In particular, focusing qualitycan be assessed by means of the off-diagonal elements of R .To understand why, R can be expressed theoretically as fol-lows [50]: R = H (cid:62) out × Γ × H in , (5)where the symbol (cid:62) stands for transpose. The matrix Γ de-scribes the scattering process inside the medium. In the sin-gle scattering regime, Γ ( z ) is diagonal and its coefficientsmap the local reflectivity γ ( r ) of the subsoil. H out and H in arethe output and input focusing matrices, respectively. Theircolumns correspond to the transmit or receive point spreadfunctions (PSFs), i.e. the spatial amplitude distribution ofthe focal spots around the focusing point r in or r out . Forspatially-invariant aberration, H out ( r , r out ) = H out ( r − r out ) and H in ( r , r in ) = H in ( r − r in ) . The previous equation canthen be rewritten in terms of matrix coefficients as follows: R ( r out , r in ) = (cid:90) d r H out ( r − r out ) γ ( r ) H in ( r − r in ) . (6)This last equation confirms that the diagonal coefficients of R , i.e. an horizontal slice of the confocal image, result from aconvolution between the medium reflectivity γ and the prod-uct of the input and output PSF, H out × H in .Interestingly, the off-diagonals terms in the reflection ma-trix can be exploited to estimate the imaging PSF, and therebyassess the quality of focusing. To that aim, the relevant ob-servable is the intensity distribution along each antidiagonalof R , I ( r m , ∆ r ) = R ( r m − ∆ r , r m + ∆ r )= (cid:90) d r (cid:48) H out ( r (cid:48) − ∆ r ) γ ( r (cid:48) + r m ) H in ( r (cid:48) + ∆ r ) . (7)All couple of points on a given antidiagonal have the samemidpoint r m = ( r out + r in ) / , but different spacings ∆ r =( r out − r in ) / . Whatever the nature of the scattering medium,the common midpoint intensity profile is a direct indicatorof the local PSF. However, its theoretical expression differsslightly depending on the characteristic length scale l γ of thereflectivity γ ( r ) at the ballistic depth and the typical width δ (0) in/out of the PSFs [50]. In the specular scattering regime( l γ >> δ (0) in/out ), the common-midpoint amplitude is directlyproportional to the convolution between the coherent outputand input PSFs, [ H out ⊗ H in ] (2∆ r ) (the symbol ⊗ standsfor convolution). In the speckle regime ( l γ << δ (0) in/out ),the common midpoint intensity I ( r m , ∆ r ) is directly propor-tional to the convolution between the incoherent output andinput PSFs, (cid:2) | H out | ⊗ | H in | (cid:3) (2∆ r ) . In the present case,the subsoil of SJFZ can be assumed as a sparse scatteringmedium. It means that only a few bright and coherent re-flectors emerge at each depth. This hypothesis will be ver-ified a posteriori with the three-dimensional image we will a ∆ y [ m ] ∆ x [m] c ddb r in r out b y x m m b - - - -60-40-200 T i m e [ s ] D e p t h [ m ] . . y [ m ] - x [ m ] d FIG. 2. (a) Focused reflection matrix R in the focal plane at depth z = R in (a)whose common mid-point exhibits the maximum confocal signal. The white circle accounts for the theoretical resolution cell imposed by thegeophone array aperture. (c) Confocal image extracted from the diagonal of R in (a). (d) Vertical slice of the 3D confocal image obtainedby 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 islinear and applies to panels (a,c). The color scale on the right is in dB and applies to panels (b,d). obtain. For an isolated scatterer, the common mid-point in-tensity at its position scales as the product between the twoPSFs, H out ( r m − ∆ r ) × H in ( r m − ∆ r ) . Therefore, the en-ergy spreading in the vicinity of each scatterer position shallenable one to probe the spatial extension of the PSF. As thescatterer position is a priori unknown, the imaging PSF willbe, in practice, probed by considering the antidiagonal whosecommon mid-point exhibits the maximum confocal signal.Figure 2b shows the corresponding common midpoint in-tensity profile for the matrix R displayed in Fig.2a. It shows asignificant spreading of energy over off-diagonal coefficientsof R . This effect is a direct manifestation of the aberrationssketched in Fig. 1b. Indeed, in absence of aberration, all theback-scattered energy would be contained in a diffraction-limited confocal focal spot, H (∆ r ) = sinc ( π ∆ r/δ ) , with δ = λ/ (2 sin θ ) . The ideal -6dB main lobe width (or fullwidth at half maximum) is roughly equal to δ ∼ m. Thisdiffraction-limited resolution is depicted by a white circle inFig. 2b. Here the characteristic size of the main central lobeis δ (0) in/out ∼ z = 3600 m. Hence, the back-scatteredenergy spreads well beyond the diffraction limit. Besides thiscentral lobe, few secondary lobes also emerge in Fig. 2b dueto the gap between the velocity model and the actual seismicwave velocity distribution in SJFZ. As shown in Supplemen-tary Section S2, these secondary lobes are strongly affectedby our choice of c . Hence this observable can be used foroptimizing our wave propagation model. As shown by Sup-plementary Fig. S2, the value c = 1500 m/s is the seismicwave velocity that clearly minimizes the level of these sec-ondary lobes. Despite this optimization, the focusing quality remains farfrom being ideal because of the heterogeneous distribution of c in the subsoil. In the following, we will show how this fun-damental issue can become a strength since it can enlarge vir-tually the aperture angle under which the geophone array isseen, thereby leading to an enhanced spatial resolution. II. DISTORTION MATRIX
To that aim, a new operator is introduced: The so-calleddistortion matrix D [10, 11]. This operator essentially con-nects each virtual geophone with the distortion exhibited bythe associated wave front in the far-field. The D -matrix isthus equivalent to a reflection matrix but in a moving frame, i.e centered around each input focusing beam. This change offrame will allow us to unscramble the contribution of phaseaberrations from the medium reflectivity. Last but not least,it will be shown to be particularly efficient for spatially dis-tributed aberrations. While conventional adaptive focusingtechniques 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 ad-dress them independently. A. Reflection matrix in a dual basis
The reflection matrix R is first projected into the Fourierbasis in reception. To that aim, we define a free-space trans-mission matrix T which corresponds to the Fourier trans- ^ Δk out/in r r r b c d e wavefront alignment iterative time reversal processing coherent virtual scatterer point-li k e virtual scatterer r r r Rinput & output focusing output/inputadaptive focusing r r r a r r r de-scan of input-focal spots SVD of δ C in/out distortion matrix D out / D in r r r input/output exchange f R out / R in SVD of D out / D in ^ focused matrix R (p) output/intput far- f ield projection ^ δ in/out (p) U out/in(p) U out/in(p) δ FIG. 3.
Time reversal analysis of the distortion matrix . (a) In the matrix imaging scheme, each point in the focal plane is probed by meansof a focused beam at input and output. (b) A far-field projection of the focused reflection matrix [Eq. (9)] yields the dual-basis matrix R . (c)By removing the geometrical tilt of each reflected wave-front [Eq. (12)], a distortion matrix D is obtained. D is equivalent to a reflectionmatrix but with a static input PSF H ( r ) scanned by moving scatterers (Eq. 15). (d) The SVD of the distortion matrix enables the synthesis ofa virtual coherent reflector of scattering distribution | H ( r ) | (Eq. 18) and an estimation of the aberration phase transmittance ˜ H ( p ) for eachisoplanatic patch p in the field-of-view [Eq. (21)]. (e) This estimation can be refined by considering, in a second iteration of the process, theSVD of the normalized time reversal operator δ ˆC ( p ) [Eq. (32)]. 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 anovel focused reflection matrix R ( p ) . form operator. Its elements link any transverse component k = ( k x , k y ) of the wave vector in the Fourier space to thetransverse coordinate r = ( x, y ) of any point in an ideal ho-mogeneous medium: T ( k , r ) = exp ( i k . r ) . (8)Each matrix R can now be projected in the far field at its out-put via the matrix product R out = T × R , (9)Each column of the resulting matrix R out = [ R ( k out , r in )] contains the reflected wave-field in the far-field for each inputfocusing point r in . Injecting Eq. 5 into the last equation yieldsthe following expression for R out : R out = T out × Γ × H in , (10)where T out = T × H (cid:62) out is the output transmission matrixthat describe wave propagation between the focused basis andthe Earth surface in the Fourier basis. In terms of matrix co-efficients, the last equation can be rewritten as follows: R out ( k out , r in ) = (cid:90) d r T out ( k out , r ) γ ( r ) H in ( r , r in ) . (11)In a multi-target medium, the reflection matrix can be lever-aged to focus iteratively by time reversal processing on eachscatterer. This is the principle of the DORT method [Frenchacronym for Decomposition of the Time Reversal Operator;Prada et al. [6, 52]]. Derived from the theoretical study of theiterative time reversal process, it relies on the singular valuedecomposition of R out . In the single scattering regime, a one-to-one association exists between each eigenstate of R out andeach scatterer. Each singular value is directly equal to the re-flectivity γ i of the corresponding scatterer. Each output eigen-vector yields the wave-front, T ∗ out ( k out , r i ) , that should be ap- plied from the far-field in order to selectively focus on eachscatterers’ position r i . The DORT method is thus particularlyuseful for selective focusing in presence of aberrations [6, 52]or target detection in a multiple scattering regime [3, 9]. How-ever, this decomposition is of poor interest for diffraction-limited imaging since each input eigenvector yields an imageof the scatterer, H in ( r i , r in ) , that is still hampered by aberra-tions. B. Definition and physical interpretation of the distortionmatrix
To cope with the fundamental issue of imaging, a novel op-erator, the so-called distortion matrix D , has been recentlyintroduced [10, 11]. Inspired by previous works in ultrasoundimaging [1, 53], it relies on the decomposition of the reflectedwavefront into two contributions (see Fig. 3): ( i ) a geometriccomponent which would be obtained for a point-like target at r in in a perfectly homogeneous medium (represented by theblack dashed line in Fig. 3b) and which can be directly ex-tracted from the reference matrix T , and ( ii ) a distorted com-ponent due to the mismatch between the propagation modeland reality (Fig. 3c). The principle of our approach is to iso-late the latter contribution by subtracting, from the reflectedwave-front, its ideal counterpart. Mathematically, this oper-ation can be expressed as a Hadamard product between thenormalized reflection matrix R out and T ∗ , D out = R out ◦ T ∗ , (12)which, in terms of matrix coefficients, yields D out ( k out , r in ) = R out ( k out , r in ) T ∗ ( k out , r in ) . (13)The matrix D out = [ D ( k out , r in )] connects any input focalpoint r in to the distorted component of the reflected wave-fieldin the far-field. Removing the geometrical component of thereflected wave-field in the Fourier plane as done in Eq. (12)amounts to a change of reference frame. While the originalreflection matrix is recorded in the Earth’s frame (static un-derground scanned by the input focusing beam, see Fig. 3b),the matrix D out can be seen a reflection matrix in the frame ofthe input focusing beam (moving subsoil insonified by a staticfocusing beam, see Fig. 3c). C. Decomposition of the distortion matrix over isoplanaticpatches
To confirm this intuition, D out is now expressed mathemat-ically. To do so, we will decompose the matrix D out into aset of sub-matrices D ( p ) out , such that D out = (cid:80) p D ( p ) out . Eachmatrix D ( p ) out is associated with a distinct isoplanatic patch p in the field-of-view. An isoplanatic patch is, by definition,an area over which the aberrations are spatially-invariant. Ineach of them, phase distortions can be modelled by: ( i ) aspatially-invariant input PSF, H ( p ) in ( r , r (cid:48) ) = H ( p ) in ( r − r (cid:48) ) ,in the focused basis; ( ii ) a far-field phase screen of trans-mittance ˜H ( p ) out = [ ˜ H ( p ) out ( k out )] in the Fourier basis, where ˜ H ( p ) out ( k out ) = (cid:82) d r H ( p ) out ( r ) exp( − i k out r ) is the 2D Fouriertransform of the corresponding output PSF H ( p ) out ( r ) . Theoutput transmission matrix T ( p ) out can then be expressed asan Hadamard product between ˜H ( p ) out and T , the free-spacetransmission matrix, T ( p ) out = ˜H ( p ) out ◦ T . (14)In each isoplanatic patch, the injection of Eq. (10) intoEq. (12) gives the following expression for the coefficients of D ( p ) out : D ( p ) out ( k out , r in ) = (cid:90) d r (cid:48) T ( p ) out ( k out , r (cid:48) ) H ( p ) in ( r (cid:48) , ) γ ( r in + r (cid:48) ) (15)where r (cid:48) = r − r in is the position vector in a new coordi-nate system centered around the input focusing point r in . Inany isoplanatic patch p , the distortion matrix [Eq. (15)] canbe interpreted as the result of a fictitious process: Imagine abeam with γ as PSF and impinging onto a (fictitious) scat-terer at the origin, with reflectivity distribution H in ( r ) ; theresulting scattered wavefield in the direction k out would begiven by Eq. (15) [see the analogy with Eq. (11)]. The prob-lem is to isolate the propagation term T out from the scattering( γ ) and input beam ( H in ) terms, and estimate the aberration.In that context, iterative time reversal brings a solution : itis well known that it eventually converges the wavefront thatovercomes aberrations and optimally focuses on the brightestpoint of the (fictitious) scatterer [6, 52]. In the next section,we show in details how to implement this idea, to estimateaberration for each isoplanatic patch, and finally build a high-resolution image of the subsoil. III. TIME REVERSAL ANALYSIS OF THE DISTORTIONMATRIX
At each depth, a time reversal analysis of the distortionmatrix is performed. This iterative method consists of dif-ferent steps that we will describe below. At each iteration,a virtual scatterer is synthesized at input or output throughthe distortion matrix concept. In the first setp, a singularvalue decomposition of D ( p ) out decomposes the field-of-view ( i.e. , the transverse size of the focal plane) into a set of iso-planatic patches. The corresponding eigenvectors yield an es-timation of the aberration transmittance ˜ H ( p ) out over each iso-planatic patch. Their phase conjugate provide the focusinglaws that enable a (partial) compensation for the phase dis-tortions undergone by the incident and reflected waves duringtheir travel between the geophone array and the focal plane.By alternatively applying the same aberration correction pro-cess at input, the size of the virtual scatterer can be graduallyreduced. It converges to optimal focusing laws that will, ul-timately, provide a high-resolution mapping of the SJFZ sub-soil. A. Output distortion matrix and isoplanatic patches
Mathematically, the iterative time reversal process is equiv-alent to the singular value decomposition of D out that can bewritten as follows: D out = U out × Σ × V † in (16)or, in terms of matrix coefficients, D ( k out , r in ) = N (cid:88) p =1 σ ( p ) out U ( p ) out ( k out ) V ( p ) ∗ in ( r in ) . (17) Σ is a diagonal matrix containing the real positive singularvalues σ p in a decreasing order σ > σ > .. > σ N . U out and V in are unitary matrices whose columns, U ( p ) out = [ U ( p ) out ( k out )] and V ( p ) in = [ V ( p ) in ( r in )] , correspond to the output and inputsingular vectors, respectively. In the present case, we ex-pect a one-to-one association between each sub-matrix D ( p ) out and each eigenstate of D out . Each isoplanatic patch beingmapped onto a single virtual scatterer, the associated matrix D ( p ) out should be, in first approximation, of unity rank.To check this intuition, the correlation matrix (or time-reversal operator), C ( p ) out = D ( p ) out × D ( p ) † out , can be investi-gated. Its eigenvalue decomposition is actually equivalent tothe SVD of D ( p ) out [10, 11, 54]. Using Eq. (15), the coefficientsof each matrix C ( p ) out can be expressed as follows: C p out ( k , k (cid:48) ) = ρ p ˜ H p out ( k ) ˜ H ( p ) ∗ out ( k (cid:48) ) (cid:104) ˜ H ( p ) in (cid:126) ˜ H ( p ) in (cid:105) ( k − k (cid:48) ) , (18)where ρ p is the overall patch reflectivity and the symbol (cid:126) stands for a correlation product. The correlation term inEq. (18) results from the Fourier transform of the scattering
15 200.080.090.10.110.120.130.14 a rank n o r m a l i z e d s i n g u l a r v a l u e s b c d e V in | F o c a l p l a n e F o c a l p l a n e F o u r i e r d o m a i n F o u r i e r d o m a i n (1) | V in | (2) | U out | (1) | U out | (2) FIG. 4. Singular value decomposition of the distortion matrix D out at time t = 4 . s and depth z = ˆ σ p . (b,c) Modulus of input eigenvectors V (1) in and V (2) in . (d,e) Modulus of output eigenvectors U (1) out and U (2) out . distribution | H ( p ) in ( r ) | of the virtual scatterer in Eq. (15). Itssupport is the correlation width ∆ k out of the aberration phasetransmittance that scales as the inverse of the spatial extension δ (0) in of the input PSF intensity | H ( p ) in | : ∆ k out = λz/δ (0) in . Ifthe virtual scatterer was point-like, this correlation term wouldbe constant and the matrix C ( p ) out would be of rank 1. Its maineigenvalue σ p is then equal to ρ p and the associated eigenvec-tor U ( p ) out directly yields the aberration transmittance ˜H p out . Inpractice, the virtual scatterer is of finite size and the correla-tion term in Eq. (18) is not negligible. The rank of C ( p ) out and D ( p ) out then scales as the number of resolution cells mappingthe virtual scatterer [55, 56]: Q = ( δ (0) in /δ ) . Q is typi-cally equal to for the local PSF | H ( p ) in | displayed in Fig.2b.Among these Q eigenmodes of D ( p ) out , only the first eigenstate U ( p ) out is of interest since it maximizes the backscattered energyby focusing at the center of the virtual scatterer.In the following, we thus expect an overall distortion matrix D out of rank P × Q with the following singular value distri-bution: ( i ) P largest singular values associated with the maineigenstates U ( p ) out of each patch p ; ( ii ) a set of ( P − × Q ofsmaller singular values associated with secondary eigenstatesthat focus on the edges of the virtual scatterers. Figure 4aconfirms this prediction by displaying the normalized singu-lar values, ˆ σ p = σ p / (cid:113)(cid:80) Ni =1 σ i , of D out at depth z = 3600 m( t = 4 . s). As expected, a few singular values potentially as-sociated with the signal subspace seem to predominate over acontinuum of eigenvalues characteristic of the noise subspace.However, it is difficult to determine the effective rank of thesignal subspace. This issue can be circumvented by comput-ing the Shannon entropy (cid:72) of the singular values [10, 11],such that (cid:72) = − N (cid:88) i =1 ˆ σ i log (cid:0) ˆ σ i (cid:1) . (19)The Shannon entropy can be used as an indicator of the rankof the signal subspace. In the case of Fig.4, the singular valuesof D out display an entropy (cid:72) (cid:39) . . As the effective rank of D out scales as P × Q and that Q ∼ , this means that thenumber P of isoplanatic patches is roughly equal to .Figure 4 shows the corresponding eigenstates of D out . The two first input singular vectors, V ( p ) in , exhibit a rough imageof the corresponding isoplanatic patch p (see Figs.4b and c).This observation can be explained by deriving their theoreticalexpression [11]: V ( p ) in ( r in ) = (cid:90) d r γ ( r ) H ( p ) in ( r , r in ) . (20)The input eigenvectors provide an underground image stillpolluted by input aberrations. They are thus of limited inter-est. On the contrary, the output eigenvectors U ( p ) out display thecorresponding phase distortion in the Fourier plane [11]: U ( p ) out ( k out ) ∝ ˜ H ( p ) out ( k out ) (cid:104) ˜ H ( p ) in (cid:126) ˜ H ( p ) in (cid:105) ( k out ) . (21)The eigenvectors U ( p ) out provide the aberration trans-mittance ˜ H ( p ) out modulated by its autocorrelation function (cid:104) ˜ H ( p ) in (cid:126) ˜ H ( p ) in (cid:105) (see Fig 3d). This last term tends to limitthe support of the eigenvector U ( p ) out to ∆ k out . Figures 4dand e confirm this theoretical prediction by showing the mod-ulus of U (1) out and U (2) out , respectively. Both singular vectorscover a restricted and different angular domain in the Fourierspace ( ∆ θ out = ∆ k out /k ∼ o ). To circumvent this issue,one trick is to use only the phase of these eigenvectors U ( p ) out (Figs. 5a and b) by considering the normalized vector ˆU ( p ) out ,such that ˆ U ( p ) out ( k out ) = U ( p ) out ( k out ) / | U ( p ) out ( k out ) | . (22)Indeed, if we make the realistic hypothesis of a real and pos-itive autocorrelation function (cid:104) ˜ H ( p ) in (cid:126) ˜ H ( p ) in (cid:105) in Eq. (21), thenormalized vector ˆU ( p ) out should, in principle, provide the aber-ration transmittance ˜H ( p ) out . In practice, noise degrades the es-timation of ˜H ( p ) out outside the coherence area ∆ k out . In the fol-lowing, we will note δ ˜H ( p ) out = ˆU ( p ) ∗ out ◦ ˜H ( p ) out , the residual phasemismatch between ˆU ( p ) out and ˜H ( p ) out .Despite this phase mismatch, the phase conjugate of ˆU ( p ) out can be used as a focusing law to compensate (at least partially)0 P S F P S F a be fc d m m m m ° ° +π−π U out(2) * U out(1) * ^ ^ V in(1) ^ V in(2) ^ FIG. 5. Output phase distortion correction at time t = 4 . s anddepth z = ˆU (1) out and ˆU (2) out , re-spectively. (c,d) Imaging PSF extracted from the main antidiago-nals of the focused reflection matrices R (1) and R (2) , respectively(Eq. 25). The color scale is in dB. (e,f) Confocal images built frominput vectors ˆV (1) in and ˆV (2) in , respectively [Eq. (24)]. for phase distortions at the output. The resulting input vector, ˆV ( p ) in = ˆU ( p ) † out D out , (23)yields an image of the subsoil reflectivity. ˆV (1) in and ˆV (2) in aredisplayed in Figs. 5e and f, respectively, at depth z = 3600 m.The comparison with the initial confocal image (Fig. 2c) illus-trates the benefit of our matrix approach. While the originalimage displays a random speckle feature across the field-of-view, ˆV (1) in reveals a complex structure in the vicinity of thesurface traces of the Clark Fault. On the contrary, the image ˆV (2) in spans along an oblique direction compared to the direc-tion of the fault. While it shows a more diffuse image of theClark Fault at the center of the field-of-view, ˆV (2) in reveals astrong scattering structure on the west of the fault. The differ-ence between these images and the input eigenvectors V ( p ) in (Figs. 4b and c) can be understood by expressing the coeffi-cients of ˆV ( p ) in using Eqs. 14, 15 and 21: ˆ V ( p ) in ( r in ) = (cid:90) d r δH ( p ) out ( r , r in ) γ ( r ) H ( p ) in ( r , r in ) . (24)where δH ( p ) out is the novel output PSF related to the resid-ual phase transmittance δ ˜ H ( p ) out through a Fourier transform.While the input eigenvectors V ( p ) in display an incoherent im-age of the subsoil strongly hampered by input aberrations (Eq. 20), the vectors ˆV ( p ) in yield a confocal image with an opti-mized focusing at the output. A gain in resolution and contrastis thus obtained and confirms our previous analysis.To quantify the gain in image quality, the local imaging PSFshould be investigated. To that aim, one can notice that theinput vector ˆV ( p ) in is actually equivalent to the confocal imagethat would be extracted from a focused reflection matrix R ( p ) after output aberration compensation (Fig. 3f): R ( p ) = ( ˆU ( p ) out ◦ T ) † × R out (25)The focused matrix is useful since its antidiagonals canbe used to probe the local imaging PSF across the field-of-view [Eq. (7)]. The resulting PSFs are displayed in Figs. 5cand d. It should be compared with the initial imaging PSF(Fig. 2b). While the original PSF exhibits a distorted centrallobe that spans over almost four diffraction-limited resolutioncells (white circle in Fig. 2b), the corrected PSF shows a cen-tral lobe thinner than the diffraction limit. This striking resultwill be discussed further. Note, however, the occurence of astrong secondary lobe and an incoherent background that canbe accounted for by the subsistence of input aberrations. Thelatter ones are tackled in the next section. B. Input distortion matrix and super-resolution
Aberrations undergone by the incident waves are now com-pensated by building input distortion matrices D ( p ) in (Fig. 3c).To that aim, the reflection matrices R ( p ) should be projectedin the Fourier plane at emission: D ( p ) in = T † ◦ (cid:16) R ( p ) × T (cid:62) (cid:17) . (26)The SVD of each matrix D ( p ) in is performed. Among all theoutput eigenvectors of these matrices, one is of special in-terest: their p th eigenvector U ( p ) in . Its modulus is shown for p = 1 and in Fig. 6a and b, respectively. The input eigenvec-tors U ( p ) in exhibit a much larger angular aperture ( ∆ θ in ∼ o )than the output eigenvectors U ( p ) out ( ∆ θ out ∼ o ). To under-stand this effect, the eigenvectors U ( p ) in can be expressed, infirst approximation, as follows [11]: U ( p ) in ( k in ) ∝ ˜ H ( p ) in ( k in ) ˜ H ( p ) in ( k in ) (cid:104) δ ˜ H ( p ) out (cid:126) δ ˜ H ( p ) out (cid:105) ( k in ) . (27)As mentioned before, the correlation term, δ ˆ H ( p ) out (cid:126) δ ˆ H ( p ) out ,can be seen as the Fourier transform of the virtual scatterersynthesized from the output focal spots in the distortion matrix D ( p ) in (see Fig. 3c). The width δ ( p ) out of the corrected output PSF δH ( p ) out (Figs. 5c and d) being much thinner than the originalone δ (0) in at the input (Fig. 2b), the correlation width ∆ k ( p ) in ∼ λz/δ ( p ) out of the incident wave-field is larger than the outputcorrelation width ∆ k ( p ) out ∼ λz/δ (0) in of the reflected wave-fieldat the previous step.1 P S F P S F ca bg he f m m m ° ° ° ° d m +π−π U in | (1) | U in | (2) U in(1) * ^ U in(2) * ^ V out(1) V out(2) ^ ^ FIG. 6. Input phase distortion correction at time t = 4 . s and depth z = U (1) in and U (2) in .(c,d) Output focusing laws derived from the phase conjugation ofnormalized eigenvectors ˆU (1) in and ˆU (2) in , respectively. (e,f) ImagingPSF extracted from the main antidiagonals of the focused reflectionmatrices R (1) and R (2) , respectively [Eq. (30)]. The color scale is indB. (g,h) Confocal images built from input vectors ˆV (1) out and ˆV (2) out ,respectively [Eq. (29)]. The phase of the first input singular vectors U ( p ) in (Figs. 6cand d) are thus better estimators of ˜H ( p ) in than the original in-put eigenvectors U ( p ) out . The phase mismatch between themwill be noted δ ˜H ( p ) in in the following. Despite this residualphase error, the normalized eigenvectors, ˆU ( p ) in , can be usedas focusing laws to compensate for input aberrations. The re-sulting output vectors, ˆV ( p ) out = D ( p ) in × ˆU ( p ) ∗ in , (28)yield images of the subsoil reflectivity. The result is displayedin Figs. 6g and h at depth z = 3600 m. Their comparison withthe previous images (Figs. 5e and f, respectively) illustratesthe benefit of a simultaneous aberration correction at inputand output. While the first vector ˆV (1) out yields a refined im-age of the complex structure lying in the first isoplanatic patch(Fig. 6g), the second eigenvector ˆV (2) out now clearly highlights a coherent reflector belonging to a second isoplanatic patchon the west of the fault (Fig.6h). The gain in image qualitycan be understood by expressing the coefficients of ˆV ( p ) out us-ing Eqs. (14), (15) and (21): ˆ V ( p ) out ( r in ) = (cid:90) d r δH ( p ) out ( r , r in ) γ ( r ) δH ( p ) in ( r , r in ) . (29)where δH ( p ) in is the novel input PSF related to the residualphase transmittance δ ˜ H ( p ) in . While the previous images ˆV ( p ) in still suffer from input aberrations [Eq. (24)], ˆV ( p ) out yields anearly perfect confocal image.The resolution and contrast enhancements are now quanti-fied. As before, it can be done by extracting the local imagingPSF from the antidiagonals of the updated reflection matrices R ( p ) after input and output aberration compensation: R ( p ) = (cid:104) U ( p ) ∗ in ◦ D ( p ) in (cid:105) × T ∗ (30)An imaging PSF can be deduced for each isoplanatic patch p by considering the antidiagonal of R ( p ) whose commonmid-point corresponds to the maximum of each image ˆV ( p ) out [Eq. (7)]. The result is displayed in Figs. 6g and h. The strongsecondary lobes exhibited by the imaging PSF at the previousstep (Figs. 5c and d) have been fully suppressed thanks to thecompensation of aberrations at the input. Strikingly, the spa-tial extension δ ( p ) in of the imaging PSF at -6dB is of the orderof 150 m. This value is much thinner than the expected res-olution cell δ ∼ m, depicted as a white dashed line inFigs. 6e and f. The observed resolution is one-fourth of thetheoretical limit based on the geophone array aperture. Tworeasons can be invoked to account for that surprising result.The first possibility is an overestimation of the wave speed c atshallow depth. This value actually dictates the maximum spa-tial frequency exhibited by body waves induced at the Earthsurface. However, the chosen value c = 1500 m/s seems inagreement with recent measurements of the P-wave velocityin SJFZ [44]. The other explanation is that the heterogeneoussubsoil acts as a lens between the geophones and the focalplane either by scattering or wave guiding. As in a time rever-sal experiment, scattering by small-scale heterogeneities canwiden the angular aperture of the focused beam [57]. Alter-natively, the local gradient of the background seismic velocityin the vicinity of the fault can constitute a wave guide throughwhich seismic waves can be channeled [58]. Similarly to scat-tering, reflections on wave guide boundaries can also enlargethe effective array aperture: This is the so-called kaleido-scopic effect [59]. In both cases, if the induced aberrationsare properly compensated, the spatial extension δ ( p ) in/out of theimaging PSF can be improved compared to the diffraction-limited resolution δ predicted for a homogeneous medium.2 P S F P S F a bie f j ° ° c dg h o u t p u t o u t p u t P S F i n p u t P S F i n p u t +π−π m m m m δ U out(1) * ^ δ U in(1) * δ U out(2) * ^ δ U in(2) * δ V out(1) ^ δ V out(2) ^ FIG. 7. Correction of residual phase distortions by making the virtualscatterer point-like. (a,b,c,d) Output and input focusing laws derivedfrom the phase conjugation of the normalized eigenvectors δ ˆU (1) out , δ ˆU (2) out , δ ˆU (1) in and δ ˆU (2) in , respectively. (e,f,g,h) Imaging PSFs ex-tracted 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 confocalimages built from the diagonal of the final reflection matrices R ( p ) ,respectively [Eq. (36)]. C. Normalized correlation matrix and compensation forhigh-order aberrations
Despite these remarkable properties, the imaging PSFs inFigs. 6e and f still exhibit high-order aberrations that resultin an incoherent background of -20 dB beyond the resolutioncell. To compensate for these residual aberrations and againimprove the image quality, a last step consists in consideringa normalized correlation matrix in order to make the virtualscatterer point-like [10] (see Fig. 3d).It first consists in projecting the updated reflection matrix R ( p ) (Eq. 30) in the Fourier basis at the output: R ( p ) out = T × R ( p ) . Then, a novel output distortion matrix δ D ( p ) out from thefocused reflection matrix R ( p ) resulting from Eq. 30: δ D ( p ) out = R ( p ) out ◦ T ∗ , (31)The corresponding correlation matrix, δ C ( p ) out = δ D ( p ) out × δ D ( p ) † out , should be investigated. By analogy with Eq. (18), its coefficients can be expressed as follows: δC ( p ) out ( k out , k (cid:48) out ) ∝ ρ p δ ˜ H ( p ) out ( k out ) δ ˜ H ( p ) ∗ out ( k (cid:48) out ) × (cid:104) δ ˜ H ( p ) in (cid:126) δ ˜ H ( p ) in (cid:105) ( k out − k (cid:48) out ) . (32)As previously highlighted in Sec. III A [Eq. (18)], the correla-tion term in Eq. (32) prevents a proper estimation of the aber-ration phase transmittance over the whole angular spectrum.To circumvent that issue, the correlation matrix coefficientscan be normalized as follows [10] δ ˆ C ( p ) out ( k out , k (cid:48) out ) = δC ( p ) out ( k out , k (cid:48) out ) / | δC ( p ) out ( k out , k (cid:48) out ) | . (33)As illustrated by Fig. 3d, this operation makes the virtualreflector point-like [10]. Indeed, if we make the realis-tic hypothesis of a real and positive autocorrelation function (cid:104) δ ˜ H ( p ) in (cid:126) δ ˜ H ( p ) in (cid:105) in Eq. 32, the coefficients of δ C ( p ) out can ac-tually be expressed as follows: δC ( p ) out ( k out , k (cid:48) out ) ∝ δ ˜ H ( p ) out ( k out ) δ ˜ H ( p ) ∗ out ( k (cid:48) out ) . (34)Such a matrix is equivalent to the time reversal operator asso-ciated with a point-like reflector at the origin [6, 52]. In thatcase, the matrix δ ˆC ( 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 ) (35)The phase of the eigenvectors δ U ( p ) out are displayed in Fig. 7aand b, respectively. Compared to the Fresnel zone fringesexhibited by the first-order aberrations in U ( p ) out (Fig. 5a andb), δ U ( p ) out displays higher-order aberrations, thereby leadingto more complex focusing law. It can be directly applied tothe entries of R ( p ) out to compensate for residual aberrations. Anupdated focused reflection matrix R ( p ) is then obtained R ( p ) = T † × (cid:104) δ U ( p ) ∗ out ◦ δ D ( p ) out (cid:105) (36)The resulting focused reflection matrices enables an estima-tion of the imaging PSFs displayed in Fig. 7e and f over eachisoplanatic patch. This second correction step drastically re-duces the spatial extension of the PSF width compared to theprevious step (Figs.6e and f). The incoherent background isbelow -30 dB beyond the theoretical resolution cell ( δ = 600m at depth z = 3600 m). Moreover, the FWHM of the PSF isnow of the order of the wavelength: δ ( p ) out ∼ λ ∼ m. Strik-ingly, the observed resolution is one-sixth of the theoreticallimit based on the geophone array aperture.As before, the process can be iterated by exchanging the fo-cused and Fourier basis at input and output. A novel input dis-tortion matrix δ D ( p ) in is built for each isoplanatic patch p . Ad-ditional input focusing laws δ U ( p ) ∗ in are extracted through theeigenvalue decomposition of the normalized correlation ma-trix δ ˆC ( p ) in . Their phases are displayed in Figs. 7c and d. Theyare flatter than their output counterparts (Figs. 6c and d), indi-3cation that our approach gradually converges towards a finiteaberration phase law. The normalized eigenvectors, δ ˆU ( p ) in ,can be used as input focusing laws that should be applied to D ( p ) in to compensate for residual input aberrations. The result-ing images of the subsoil reflectivity are displayed in Figs. 7i and j at z = 3600 m over the two main isoplanatic patches.Their comparison with the previous images (Figs. 6g and h,respectively) illustrates the benefit of iterating the aberrationcorrection process. The first isoplanatic patch highlights thepresence of three distinct reflectors that arise in the vicinityof the fault’s surface traces (Fig.7i). The second isoplanaticpatch yields a highly-resolved image of a coherent reflector onthe west of the Clark fault (Fig.7i). Compared to the previousimage in Fig. 7j, the contrast is drastically improved. This ob-servation can be understood by looking at the correspondingPSFs (Figs.7g and h). The incoherent background is below-35 dB beyond the theoretical resolution cell ( δ = 600 m atdepth z = 3600 m). Last but not least, the FWHM of the PSFnearly reaches the diffraction limit δ ( p ) in ∼ m. Remarkably,the observed resolution is nearly one-eighth of the theoreti-cal limit based on the geophone array aperture. As mentionedin the previous section, the scattering events induced by thesmall-scale heterogeneities or the wave guiding induced bythe gradient of the seismic velocity around the fault can bothaccount for this unexpected resolution.The comparison of these final images with the original one(Fig. 2c) illustrates the relevance of a matrix approach for deepseismic imaging. In the next section, the three-dimensionalstructure of the SJFZ is revealed by combining the imagesderived at each depth. At last, based on the derived high-resolution 3D images, a structural interpretation of the SJFZis provided. IV. THREE-DIMENSIONAL IMAGING OF THE SANJACINTO FAULT ZONE
Having demonstrated in Sec. III how to correct phase dis-tortions at each depth, a 3D image of the subsurface can nowbe uncovered. To that aim, each isoplanatic patch should berecombined at each depth. The resulting image (cid:73) M is a co-herent sum between the diagonal coefficients of the correctedreflection matrices R ( p ) : (cid:73) M ( r , z ) ≡ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) P (cid:88) p =1 R ( p ) ( r , r , z ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (37)Figure 8 shows a slice of this 3D image with the same ori-entation as the original b-scan displayed in Fig. 2c. Whilethe raw image (cid:73) is completely blurred, the gain in resolutionprovided by the matrix aberration correction process revealsan image of the SJFZ subsoil with a refined level of details.Based on the surface traces displayed in Fig. 1a, the fault isexpected to spread orthogonally to the x -axis. The choice ofthe profiles’ orientation in Figs. 2c and 8 is dictated by ourwillingness to highlight the fault blocks on the right and leftside of the slip interface. -70-50-30-101030 T i m e [ s ] D e p t h [ m ] . . y [ m ] - x [ m ] FIG. 8. Vertical slice of the 3D confocal image obtained from stackedcorrected images derived at each depth. The slice orientation is cho-sen to be normal to the fault plane. The color scale is in dB.
The SJFZ is a major continental strike-slip fault system thataccumulated through history a significant slip of tens of kilo-meters. Such large fault systems induce zones of stronglydamaged materials [60, and references therein] and the dam-age is expected to be pronounced in the shallow crust. Thanksto the distortion matrix correction, Fig. 8 shows the two mainsignatures of the fault in diffraction imaging. First, it revealsa damage zone associated with a dense distribution of scatter-ers and highlighted by a shaded area in Fig. 8. The damagevolume extends from the surface to 2000 m below with a sec-tion decreasing in depth. Second, beyond a depth of 2000m, the matrix image displays strong echoes associated withthe presence of sub-horizontal reflectors. The correspondingstrata layers are located at different depths on both sides of thefault. This discontinuity seems to indicate the fault location indepth, although the damage area is no longer discernible be-yond 2000 m. This example thus shows both the efficiency ofour approach to cope with phase distortions and the potentialof passive seismic imaging to study the fine structure of activefaults in depth.In this section, only a slice of the 3D volume has beenshown to illustrate the drastic improvements granted by a ma-trix approach of seismic imaging (Fig. 8). Its comparison withthe raw confocal image highlights the benefit of a drastic gainin resolution (Fig. 2c). A more detailed interpretation of thesubsurface properties will be the focus of a future study. In-deed, a detailed interpretation of the obtained 3D image and a4confrontation with previous studies may provide informationon the structural properties of the fault zone. More gener-ally, the high resolution and penetration depth of our matriximaging method will play an important role in detecting ac-tive faults, evaluating their long-term behaviour and, conse-quently, following the deformation of the Earth’s crust. It willalso be an essential tool for understanding earthquake physicsand evaluating seismic hazard.
V. CONCLUSION
Inspired by pioneering works in optical microsopy [5, 11]and ultrasound imaging [1, 10], a novel matrix approach toseismic imaging is proposed in this paper. Taking advantageof the reflection of bulk seismic waves by heterogeneities indepth, it can be applied to both active or passive seismic imag-ing. The strength of this approach lies in the fact that it workseven when the velocity distribution of the subsoil is unknown.By projecting the seismic data either in the focused basis or inthe Fourier plane, it takes full advantage of all the informationcontained in the collected data. For aberration correction, pro-jection of the reflection matrix into a dual basis allows the iso-lation of the distorted component of the reflected wave-field.Seen from the focused basis, building the distortion matrix D consists in virtually shifting all the input or output focalspots at the onto the same virtual location. An iterative timereversal analysis then allows to unscramble the phase distor-tions undergone by the wave-front during its travel betweenthe Earth surface and the focal plane. A one-to-one associa-tion is actually found between each eigenstate of D and eachisoplanatic patch in the field-of-view. More precisely, eachsingular vector in the Fourier space yields the far-field focus-ing law required to focus onto any point of the correspondingisoplanatic patch. A confocal image of the subsoil reflectiv-ity can then be retrieved as if the underground had been madeperfectly homogeneous. However, it should be noted that theavailable velocity model is often a rough approximation of re-ality and the estimated depth is likely to be inaccurate. Yeta correct model would dilate the subsboil 3D image up ordown but would not change significantly its transverse evo-lution. Moreover, the matrix approach can be easily incorpo-rated in any iterative procedure that aims at mapping the ve-locity distribution of bulk seismic waves in the underground.This would, in turn, correct the image from its axial aberra-tions.In this paper, as a proof-of-concept, the case of the San Jac-into Fault zone is considered. While a raw confocal image suf-fers from an extremely bad resolution due to the strong lateralvariations of the seismic velocities in the vicinity of the fault,our matrix approach provides a high-resolution image. Strik-ingly, its resolution is almost one eighth below the diffractionlimit imposed by the geophone array aperture. The hetero-geneities of the subsoil play the role of a scattering and/orchanneling lens which increases drastically the effective arrayaperture. The matrix image reveals a damage volume partic-ularly pronounced in the shallow crust (¡2000 m). At largerdepth, the 3D image of the fault exhibits the fault blocks on the right and left side of the slip interface. A more detailed in-terpretation of the obtained image will be the focus of a futurestudy.Besides seismic fault zones, a matrix approach of passiveimaging is particularly suited to the study of volcanoes. Whilethe multiple scattering problem in volcanoes has been recentlytackled under the reflection matrix approach [9], the distortionmatrix would be particularly powerful to restore a diffraction-limited image of the in-depth structure of a volcano. Addi-tionally, our matrix method can advantageously be applied toreal-time geophysical imaging, similarly to adaptive optics inastronomy. This may provide a valuable tool for monitor-ing oil or gas reservoirs when drilling, extracting hydrocar-bon or injecting CO . Finally, this method can be taken ad-vantage of for marine exploration where the variations of sealocal properties can distort the acoustic wave-front depend-ing on the gradients of temperature and salinity, not to men-tion streams and swirls. The matrix approach presented herecan compensate for the phase distortions undergone by the re-ceived echoes and improve the interpretation of the data withno more required knowledge than the speed of sound in wa-ter. Although not leveraged in the present paper, Lambert etal. [10] have also demonstrated the relevance of the matrixapproach at removing multiple reverberations between stratalayer interfaces. This also may be useful for any marine orseismic exploration purpose. ACKNOWLEDGMENTS
The authors are grateful for funding provided by LABEXWIFI (Laboratory of Excellence within the French Pro-gram Investments for the Future, ANR10LABX24 andANR10IDEX000102 PS and by TOTAL R&D. We acknowl-edge the support from the European Research Council (ERC)under the European Union Horizon 2020 research and inno-vation program (grant agreement No 742335, F-IMAGE andgrant agreement No. 819261, REMINISCENCE).
DATA AVAILABILITY
The seismic data used in this study can be obtained fromthe Data Management Center of the Incorporated ResearchInstitutions for Seismology.5
Supplementary Information
This supplementary material includes details about: ( i ) theprior filtering applied to the reflection matrix in the geophonebasis; ( ii ) the choice of the wave velocity model based on thefocused reflection matrix. S1. PRIOR FILTERING OF THE REFLECTION MATRIX.
Figure S1a displays the response matrix K ( t ) at timelag t = 0 . s. Surprisingly, this matrix is dominatedby a predominant signal along its diagonal. To under-stand the origin of this effect, let us first recall that eachcoefficient K ( u out , u in , t ) of the response matrix is com-puted from the cross-correlation of the seismic noise wave-field ψ ( u , τ ) recorded by geophones located at u = u out ) and u = u in ) . Under appropriate wave-field conditions,coda cross-correlation converges towards the Greens function G ( u out , u in , t ) between receiving stations as if one of them( u in ) had become a source: K ( u out , u in , t ) = lim T →∞ T (cid:90) T 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 sig-nal K ( u in , u in , τ ) gives rise to an intense peak at lag time t = 0 , that physically corresponds to the seismic noise powerspectral density measured at point u in . Only the non zerolag time contribution carries information on the reflectivityat depth and is thus of interest for imaging purposes. How-ever, this zero time peak gathers most of the energy content inthe retrieved signal, and can have a detrimental impact on ouranalysis. In the present situation, the limited frequency band-width makes this initial pulse duration far from being negligi-ble ( δt ∼ / ∆ f = λ/ . Given the high densityof the geophones network, neighbour stations thus belong tothe same coherence area. Therefore the corresponding cross-correlation signals, that lie close to the diagonal of K , are alsodominated by this autocorrelation peak.To get rid of this central pulse, a gaussian mask is appliedto each element of the raw reflection matrix K (cid:48) ( u out , u in , t ) = K ( u out , u in , t ) × (cid:20) − exp (cid:18) − (cid:107) u out − u in (cid:107) λ (cid:19)(cid:21) where λ is the wave-length of body waves at the central fre-quency (here λ = 100 m). This filter gradually penalizes theimpulse responses computed for interstation distances smallerthan λ/ , the coherence length of the body waves in the 10Hz–20 Hz working frequency range.Figure S1b displays the filtered version K (cid:48) of the measuredreflection matrix K shown in Figure S1a. The filtered wave-field is no longer dominated by the seismic noise autocorre- lation. Note the strong gap in terms of order of magnitudebetween the diagonal coefficients of K (Fig. S1a) and the ac-tual wave-field patterns in K (cid:48) (Fig. S1b). This pre-filteringoperation is thus critical for an accurate interpretation of theeffective Earth’s response.Each column of K (cid:48) corresponds to the wavefield that wouldbe recorded by the set of geophones if a pulse were emittedfrom one geophone at u in . Figure S1c displays the wavefieldgenerated from the central geophone position and interpolatedat the surface between all receiving geophones. This wave-field clearly exhibits the contribution of direct Rayleigh wavesthat emerges along the white dashed circle. Interestingly, weobserve fast travelling waves (see for instance the top left-hand corner of the array) that correspond to the signature ofthe body waves propagating below the Earth’s surface at largerspeed than the surface waves tracked by the geophones. Theyovertake the surface waves and eventually get first to the edgeof the array, after reflection on some shallow structures or re-fraction at overcritical angles. S2. CHOICE OF THE WAVE VELOCITY MODEL.
In the accompanying paper, the reflected body waves areused to build a high-resolution in-depth image of the SJFZ.To that aim, the choice of the initial wave velocity model iscrucial. As mentioned in the accompanying paper, this choicecan be enlightened by the properties of the focused reflectionmatrix R . In particular, a common mid-point intensity profilecan be extracted from the main antidiagonal of R that exhibitsthe maximum confocal (diagonal) signal. Figure S2 showsthe corresponding intensity profile at several times of flight( t = 1 . , 1.5 and 2.65 s, from top to bottom) and for differentseismic wave velocity ( c = 1000 , 1500 and 2000 m/s, fromleft to right) in our propagation model. To make the com-parison quantitative, each intensity profile is displayed as afunction of spatial coordinates normalized by δ , the expectedresolution for each velocity model at the corresponding depth.In each panel, the focal spot consists of a main lobe centeredaround the input focusing point and a random distribution ofsecondary lobes. The dimension of the main lobe is roughlytwice the diffraction limit prediction δ , depicted by a whitecircle in each panel Fig. S2. This loss of resolution is a man-ifestation of the aberrations induced by the gap between thevelocity model and the actual seismic wave velocity distribu-tion in the SJFZ underground. However, the dimension ofthis main lobe is not significantly affected by our choice of c .Hence this observable cannot be used for optimizing our wavepropagation model. On the contrary, the amplitude and spatialextension of the secondary lobes strongly depend on the wavevelocity model. The value c = 1500 m/s is the seismic wavevelocity that clearly minimizes the level of these secondarylobes in Fig. S2.The choice is validated a posteriori by the quality of theconfocal image obtained under our matrix approach. FigureS3 shows a slice of this image for different wave velocity mod-els ( c = 1000 , 1500 and 2000 m/s). The slice orientation isthe same as the b-scans displayed in the accompanying pa-6 -300300 E a r t h s u r f a c e - y o u t [ m ] -3 Earth surface - x out [m] a b c
Geophones - Emission s in G e o ph o n e s - R e c e p � o n s o u t Geophones - Emission s in G e o ph o n e s - R e c e p � o n s o u t FIG. S1. Reflection matrix in the geophones basis. (a) Raw reflection matrix K at time t = 0 . s. (b) Filtered reflection matrix K (cid:48) at time t = 0 . s. (c) Reflected wavefield deduced from the spatial interpolation of the central column of K (cid:48) ( u in = ), depicted by a dashed orangeline in panel b. The white dashed circle accounts for the expected position of the direct surface wave-front for a Rayleigh wave speed of 350m/s [39]. ∆ x/ δ ∆ y / δ ∆ x/ δ ∆ x/ δ ∆ x/ δ ∆ x/ δ ∆ x/ δ ∆ y / δ ∆ y / δ ∆ y / δ ∆ y / δ ∆ y / δ ∆ y / δ ∆ y / δ ∆ x/ δ ∆ x/ δ ∆ x/ δ c = 1000 m/s c = 1500 m/s c = 2000 m/s t = 1.05 s t = 1.5 s t = 2.65 s ∆ y / δ -4 44 -4 44 -4 44-4 44 -4 44 -4 44-3.5 3.53.5 3.5 -3.5 3.53.5 -3.5 3.53.5 a b c FIG. S2. Imaging PSF deduced from the antidiagonal of R whose common mid-point exhibits the maximum confocal signal. This imaagingPSF is shown at different times of flight ( t = 1 . , 1.5 and 2.65 s, from top to bottom) for different wave velocity models: c = 1000 m/s (a), c = 1500 m/s (b), c = 2000 m/s (c). Each panel displays the modulus of the reflected wave-field normalized by its maximum. The whitecircle accounts for the theoretical resolution cell (disk of radius δ ) imposed by the geophone array aperture. per. While the raw image (cid:73) is completely blurred (Fig. 2c ofthe accompanying paper), a gain in resolution is provided bythe matrix aberration correction process. However, the imagequality strongly depends on the choice of the wave velocitymodel. This sensitivity can be explained by the fact that ax-ial 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.S3a and b confirms thata model with c = 1000 m/s clearly underestimates the actualbody wave velocity in the depth range considered in this work.While Fig. S3b ( c = 1500 m/s) clearly highlights strata lay-ers at different depths on both sides of the fault, Fig. S3a( c = 1000 m/s) shows a blurred image of the subsoil. For7 c = 2000 m/s (Fig.S3c), the strata structure of the SJFZ un-derground is partially revealed but with a worse resolution andlower contrast than for c = 1500 m/s (Fig.S3b) until an echotime t =4 s. Nevertheless, note that, beyond that time, the twolast interfaces seem to be better resolved for c = 2000 m/s.This indicates that, not surprisingly, the seismic wave veloc-ity increases with depth and that c = 1500 m/s is probably not the optimal wave velocity for t > s. A mapping of thewave velocity could be indeed possible through the focusedreflection matrix approach as already demonstrated in ultra-sound imaging [50]. However, the extension of this methodto seismology is beyond the scope of this paper. The matrixmapping of the bulk seismic wave velocity will be tackled infuture works. [1] J.-L. Robert and M. Fink, J. Acoust. Soc. Am. , 866 (2008).[2] A. Aubry and A. Derode, Phys. Rev. Lett. , 084301 (2009).[3] S. Shahjahan, A. Aubry, F. Rupin, B. Chassignole, andA. Derode, Appl. Phys. Lett. , 234105 (2014).[4] S. Kang, S. Jeong, H. Choi, W. Ko, T. D. Yang, J. H. Joo, J.-S.Lee, Y.-S. Lim, Q.-H. Park, and W. Choi, Nat. Photonics , 1(2015).[5] A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, andA. Aubry, Sci. Adv. , e1600370 (2016).[6] C. Prada and M. Fink, Wave Motion , 151 (1994).[7] S. M. Popoff, A. Aubry, G. Lerosey, M. Fink, A. C. Boccara,and S. Gigan, Phys. Rev. Lett. , 263901 (2011).[8] A. Aubry and A. Derode, J. Appl. Phys. , 044903 (2009).[9] T. Blondel, J. Chaput, A. Derode, M. Campillo, and A. Aubry,J. Geophys. Res.: Solid Earth , 10,936 (2018).[10] W. Lambert, L. A. Cobus, T. Frappart, M. Fink, and A. Aubry,Proc. Nat. Sci. Acad. , 14645 (2020).[11] A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, andA. Aubry, Sci. Adv. , eaay7170 (2020).[12] G. Montaldo, M. Tanter, J. Bercoff, N. Benech, and M. Fink,IEEE Trans. Ultrason., Ferroelectr., Freq. Control , 489(2009).[13] R. H. Stolt, Geophysics , 23 (1978).[14] J. Gazdag, Geophysics , 1342 (1978).[15] J. Gazdag and P. Sguazzero, Geophysics , 124 (1984).[16] J. F. Claerbout, Imaging of of the Earth’s Interior, 2nd ed.(Blackwell Scientific Publications, 1996).[17] B. Biondi, 3D Seismic Imaging (Society of Exploration Geo-physicists, 2006).[18] R. L. Weaver and O. I. Lobkis, Phys. Rev. Lett. , 134301(2001).[19] M. Campillo and A. Paul, Science , 547 (2003).[20] A. Derode, E. Larose, M. Tanter, J. de Rosny, A. Tourin,M. Campillo, and M. Fink, J. Acoust. Soc. Am. , 2973(2003).[21] K. Wapenaar, Phys. Rev. Lett. , 254301 (2004).[22] R. Snieder, Phys. Rev. E , 046610 (2004).[23] E. Larose, A. Margerin, L.and Derode, B. van Tiggelen,M. Campillo, N. Shapiro, A. Paul, L. Stehly, and M. Tanter,Geophysics , SI11 (2006).[24] N. M. Shapiro, M. Campillo, L. Stehly, and M. H. Ritzwoller,Science , 1615 (2005).[25] K. G. Sabra, P. Gerstoft, P. Roux, W. A. Kuperman, and M. C.Fehler, Geophys. Res. Lett. , n/a (2005).[26] Y. Yang, M. H. Ritzwoller, A. L. Levshin, and N. M. Shapiro,Geophys. J. Int. , 259 (2007).[27] P. Roux, K. G. Sabra, P. Gerstoft, W. Kuperman, and M. C.Fehler, Geophys. Res. Lett. (2005).[28] D. Draganov, K. Wapenaar, W. Mulder, J. Singer, andA. Verdel, Geophys. Res. Lett. (2007).[29] D. Draganov, X. Campman, J. Thorbecke, A. Verdel, andK. Wapenaar, Geophysics , A63 (2009). [30] P. Poli, H. Pedersen, and M. Campillo, Geophys. J. Int. ,549 (2012).[31] E. Ruigrok, X. Campman, D. Draganov, and K. Wapenaar,Geophys. J. Int. , 339 (2010).[32] P. Poli, M. Campillo, H. Pedersen, L. W. Group, et al., Science , 1063 (2012).[33] L. Retailleau, P. Bou´e, L. Li, and M. Campillo, Geophys. J. Int. , 1339 (2020).[34] E. Hauksson, W. Yang, and P. Shearer, Bull. Seism. Soc. Am. , 2239 (2012).[35] H. O. Johnson, D. C. Agnew, and F. K. Wyatt, J. Geophys.Res.: Solid Earth , 23951 (1994).[36] E. O. Lindsey and Y. Fialko, J. Geophys. Res.: Solid Earth ,689 (2013).[37] A. A. Allam and Y. Ben-Zion, Geophys. J. Int. , 1181(2012).[38] D. Zigone, Y. Ben-Zion, M. Campillo, and P. Roux, Pure Appl.Geophys. , 1007 (2014).[39] P. Roux, L. Moreau, A. Lecointre, G. Hillers, M. Campillo,Y. Ben-Zion, D. Zigone, and F. Vernon, Geophys. J. Int. ,980 (2016).[40] A. Allam, Y. Ben-Zion, I. Kurzon, and F. Vernon, Geophys. J.Int. , 978 (2014).[41] Y. Ben-Zion, F. L. Vernon, Y. Ozakin, D. Zigone, Z. E. Ross,H. Meng, M. White, J. Reyes, D. Hollis, and M. Barklage,Geophys. J. Int. , 370 (2015).[42] A. Mordret, P. Roux, P. Bou´e, and Y. Ben-Zion, Geophys. J.Int. , 896 (2019).[43] D. Zigone, Y. Ben-Zion, M. Lehujeur, M. Campillo, G. Hillers,and F. L. Vernon, Geophys. J. Int. , 879 (2019).[44] P.-E. Share, P. T´aboˇr´ık, P. ˇStˇepanˇc´ıkov´a, J. Stemberk, T. K.Rockwell, A. Wade, J. R. Arrowsmith, A. Donnellan, F. L. Ver-non, and Y. Ben-Zion, Geophys. J. Int. , 781 (2020).[45] G. Hillers, P. Roux, M. Campillo, and Y. Ben-Zion, J. Geophys.Res.: Solid Earth , 8048 (2016).[46] L. Qin, Y. Ben-Zion, H. Qiu, P. Share, Z. Ross, and F. Vernon,Geophys. J. Int. , 98 (2018).[47] A. J. Berkhout, Imaging of Acoustic Energy by Wave Field Extrapolation,3rd ed. (Elsevier, 1984).[48] J. R. Berryhill, Geophysics , 2064 (1984).[49] A. J. Berkhout and C. P. A. Wapenaar, J Acoust. Soc. Am. ,2017 (1993).[50] W. Lambert, L. A. Cobus, M. Couade, M. Fink, and A. Aubry,Phys. Rev. X , 021048 (2020), arXiv:eprint number here.[51] M. Born and E. Wolf, Principles of optics (Seventh edition)(Cambridge University Press, Cambridge, 2003).[52] C. Prada, S. Manneville, D. Spoliansky, and M. Fink, J. Acoust.Soc. Am. , 2067 (1996).[53] T. Varslot, H. Krogstad, E. Mo, and B. A. Angelsen, J. Acoust.Soc. Am. , 3068 (2004).[54] C. Prada and J.-L. Thomas, J. Acoust. Soc. Am. , 235(2003). T i m e [ s ] D e p t h [ m ] . . y [ m ] - x [ m ] a b c FIG. S3. Vertical slice of the 3D matrix images computed under the matrix approach for different seismic wave velocity models: c = 1000 m/s (a), c = 1500 m/s (b), c = 2000 m/s (c). The slice orientation is chosen to be normal to the fault plane. The color scale for each imageis in dB.[55] J.-L. Robert and M. Fink, J. Acoust. Soc. Am. , 218 (2009).[56] A. Aubry, J. de Rosny, J.-G. Minonzio, C. Prada, and M. Fink,J. Acoust. Soc. Am. , 2746 (2006).[57] A. Derode, P. Roux, and M. Fink, Phys. Rev. Lett. , 4206(1995).[58] Y.-G. Li and P. C. Leary, Bulletin of the Seis- mological Society of America , 1245 (1990),https://pubs.geoscienceworld.org/bssa/article-pdf/80/5/1245/2706830/BSSA0800051245.pdf.[59] P. Roux and M. Fink, J. Acoust. Soc. Am. , 2631 (2001).[60] Y. Ben-Zion and C. G. Sammis, Pure Appl. Geophys.160