Lapo Boschi

Irene Molinari

Michael Reinwald

A simple method for earthquake location by surface-wave time-reversal

The scalar 2-D Helmholtz' equation (i.e., "membrane waves") can be used to model surfacewave propagation in a laterally smooth, lossless half space. Building on this known result, we develop an algorithm to localize earthquake sources based on surface-wave data, via numerical time-reversal on a membrane, where monochromatic waves propagate with the phase velocity of Rayleigh or Love waves at the same frequency. By conducting monochromatic membrane-wave time-reversal simulations at various frequencies and combining the results, broadband time-reversed surface waves can be modeled. Importantly, membrane-wave modeling is computationally much less expensive than three-dimensional surface-wave modeling.

We first explain rigorously the relationship between surface waves and membrane waves. Our mathematical treatment is slightly different from those found in the literature, in that it does not invoke variational principles. We next implement our time-reversal algorithm via spectral elements as well as simple ray tracing. Both implementations account for the effects of lateral variations in phase velocity. We validate the two resulting tools by means of several numerical experiments. This includes synthetic tests, as well as the localization of a virtual source based on a data set of real ambient-noise cross correlations, and the localization of the epicenter of a real earthquake from real, raw data. In this study, applications are limited to Northern Italy and the Alpine arc, where we have access to recent, high resolution phase velocity maps, ambient-noise cross correlations and data from a recent, relatively large earthquake. The accuracy of epicenter location despite non-uniformity in station coverage encourages further applications of our method, in particular to the task of mapping large-earthquake rupture in space and time.

Introduction

Estimates of seismic slip as a function of position and time for a given earthquake are obtained today in different ways, depending on the magnitude and depth of the earthquake, and on the instrumental coverage. Several different types of seismic and geodetic observations are employed. Dense networks of strong-motion accelerometers are currently deployed in seismic regions worldwide; they are designed to record the high-frequency oscillations generated by a nearby event, but they have little sensitivity to the lower frequencies, and cannot be used to constrain the properties of far earthquakes. At the opposite end of the frequency spectrum, data from GPS networks and satellite geodesy are used as constraints of the final slip associated with an earthquake; they provide good resolution of the surface expression of the rupture, but have little or no sensitivity to fault geometry at depth [e.g. Mai et al., 2016, and references therein]. Wherever the coverage provided by nearby instruments is insufficient, local and/or global broadband seismic networks at teleseismic distances are used to image slip. As a general rule, fault geometry is particularly hard to constrain on the basis of seismic data alone, and is determined based on geodetic data or, wherever possible, field geology observations. Once a data set for a given event has been compiled, seismic oscillations and geodetic offsets are translated to slip on the fault via (1) least-squares inversions, (2) seismic time reversal, or (3) the back-projection method.

(1) Least-squares inversions are based on the representation theorem [e.g. [START_REF] Aki | Quantitative Seismology[END_REF], i.e. the mathematical expression of the physical law relating the geometry of an arbitrary rupture to the resulting deformation at any point of a given medium. Because the spatiotemporal evolution of seismic ruptures is generally very complex, it is not surprising that their solutions tend to be very non-unique, as shown in detail by [START_REF] Mai | The earthquake-source inversion validation (SIV) project[END_REF].

(2) The physics of acoustic or seismic time reversal can be heuristically summarized as follows: a signal is emitted by a source and recorded by multiple receivers; if receivers are then turned into sources, each emitting its own recorded signal (with the corresponding delay) flipped with respect to time, the resulting wave field will "focus" at the original source location [e.g. [START_REF] Fink | Time-reversed acoustics[END_REF]. This means that by recording real data from an unknown source and then conducting the time-reversal exercise numerically, the location of the source could be determined, provided, of course, that the error associated with modeling of propagation is small, that is to say, that the complexity of the medium of propagation is accounted for within a good approximation. From the standpoint of seismology, this amounts to a kinematic, extended-source inversion, with the additional possibility of monitoring the backward propagation of time-reversed waves before focusing at the source. In seismology, applications of time reversal [e.g. [START_REF] Larmat | Time-reversal imaging of seismic sources-application to the Sumatra earthquake[END_REF][START_REF] Larmat | Time-reversal location of glacial earthquakes[END_REF] are hindered by the high computational costs of accurate wave-propagation modeling, unless only very long periods are considered.

(3) The back-projection method as described, e.g., by [START_REF] Ishii | Extent, duration and speed of the 2004 Sumatra-Andaman earthquake imaged by the Hi-Net array[END_REF] and currently employed by many authors in seismology, is usually thought of as a simplification of wavefield reverse-time migration, a tool for imaging structure in reflection seismology. This is in many ways similar to time reversal, but involves some further, fundamental simplifications. Namely, the term back-projection refers to studies where the effects of time-reversed wave propagation are not modeled, but approximately corrected for by stacking the signals recorded by an array of nearby receivers. One of the practical consequences of this is that the physical nature of the computed, time-reversed wave field that focuses at the source remains undefined. Its interpretation in terms of rupture mechanics is complicated by the fact, e.g., that it is not known whether it more closely approximates a slip or a rate of slip [START_REF] Fukahata | Theoretical relationship between back-projection imaging and classical linear inverse solutions[END_REF].

We provide in this study the building blocks of a new algorithm for constraining extendedsource geometry and time evolution. The algorithm is based on the time-reversal concept, and thus overcomes the limitations of P -wave back-projection, but it is designed so as to reduce significantly the computational costs of full-waveform time reversal. One of the key aspects of our method is that surface waves, instead of P waves, are time-reversed and backwardpropagated. This is preferable for several reasons: (i) Surface waves are dispersive, i.e. they "spread out" along the time axis as they travel across the surface of the earth: time reversal turns this process around, enhancing the focusing of backward-propagating waves onto the source. (ii) The problem of surface-wave propagation modeling, although inherently threedimensional, can be reduced to two dimensions by separating the signal into narrow frequency bands [e.g. [START_REF] Tanimoto | Modelling curved surface wave paths: membrane surface wave synthetics[END_REF][START_REF] Tromp | Variational principles for surface wave propagation on a laterally heterogeneous Earth-III. Potential representation[END_REF][START_REF] Peter | Surface wave tomography: global membrane waves and adjoint methods[END_REF][START_REF] Peter | Tomographic resolution of ray and finite-frequency methods: a membrane-wave investigation[END_REF], to be backprojected separately, and subsequently "stacked" together: this reduces the computational costs drastically. (iii) Our knowledge of the three-dimensional structure of the Earth's deep interior, essential to backward-propagate numerically the time-reversed signal, is limited; but surface-wave propagation is confined to the upper mantle, which is relatively well known; recent, robust global phase-velocity maps of Rayleigh-and Love-wave velocities are available in the frequency band relevant to this project, at the global and, where possible, regional scales [e.g. [START_REF] Ekström | A global model of Love and Rayleigh surface wave dispersion and anisotropy, 25-250 s[END_REF][START_REF] Kaestle | Surface-wave tomography of the alps using ambient-noise and earthquake phasevelocity measurements[END_REF]. In seismology/acoustics jargon, point (iii) is equivalent to saying that very accurate surface-wave "Green's functions" are available and will be used to backward-propagate time-reversed surface-wave data. This further enhances focusing of the time-reversed wave field, and thus the robustness and resolution of mapped seismic slip.

We expect our method to be effective over a broad range of epicentral distances. At distances of 30 • or more from the epicenter, surface waves carry more energy than body waves, and they can be easily identified and isolated on seismograms. At shorter epicentral distances, where they are obscured by the body-wave coda, surface waves can still emerge in a time-reversal exercise as a result of focusing: this is confirmed by our results, discussed in sec. 6.3. Today, broadband "full-waveform" information is not routinely utilized by researchers interested in mapping the seismic source. Tentative implementations of fault imaging via seismic-waveform time reversal such as those by [START_REF] Larmat | Time-reversal imaging of seismic sources-application to the Sumatra earthquake[END_REF] and [START_REF] Larmat | Time-reversal location of glacial earthquakes[END_REF] were successful from a theoretical standpoint, but seem too computationally heavy for systematic practical application. Most seismologists only back-project seismogram peaks associated with the arrival of P waves [e.g. [START_REF] Ishii | Extent, duration and speed of the 2004 Sumatra-Andaman earthquake imaged by the Hi-Net array[END_REF]] so as to avoid costly simulations of broad-band seismic-wave propagation in a heterogeneous, three-dimensional medium (the Earth), whose heterogeneity is only approximately known. The only published experiment in surface-wave back-projection that we are aware of is that of [START_REF] Roten | A Rayleigh wave back-projection method applied to the 2011 Tohoku earthquake[END_REF]. While the basic idea of [START_REF] Roten | A Rayleigh wave back-projection method applied to the 2011 Tohoku earthquake[END_REF] is similar to some of the concepts presented here, their approach is essentially a form of back-projection, with the inherent approximations.

We provide in sec. 2 a description of surface-wave propagation in terms of "potentials" [e.g. [START_REF] Udías | Principles of Seismology[END_REF][START_REF] Aki | Quantitative Seismology[END_REF]. [START_REF] Aki | Quantitative Seismology[END_REF] state (box 7.5) that "potentials are of no direct interest, and are awkward to use. . .". We maintain that, as shown e.g. by [START_REF] Tanimoto | Modelling curved surface wave paths: membrane surface wave synthetics[END_REF] or [START_REF] Peter | Surface wave tomography: global membrane waves and adjoint methods[END_REF], there is interest in using potentials, particularly for surface waves. For instance, if only the phase, and not the amplitude of the data is studied, many useful applications (e.g., imaging, backward-propagation) become possible by using the potentials and the associated simple, two-dimensional (2-D) equations, without having to solve the more cumbersome radial equations, or the general three-dimensional equations. This is strictly true within a high-frequency approximation, but applications to real data have often shown that, in practice, this approximation works remarkably well. Our theoretical formulation in sec. 2 is different from that of [START_REF] Aki | Quantitative Seismology[END_REF] in that we use potentials, and from those of [START_REF] Tanimoto | Modelling curved surface wave paths: membrane surface wave synthetics[END_REF] and [START_REF] Tromp | Variational principles for surface wave propagation on a laterally heterogeneous Earth-III. Potential representation[END_REF] in that we do not invoke variational principles.

The main implication of sec. 2 is that the scalar 2-D Helmholtz' equation can be used to model surface-wave propagation in a laterally smooth, lossless half space, confirming earlier results by [START_REF] Tanimoto | Modelling curved surface wave paths: membrane surface wave synthetics[END_REF] and [START_REF] Tromp | Variational principles for surface wave propagation on a laterally heterogeneous Earth-III. Potential representation[END_REF]. In secs. 3 and 4 we accordingly derive the theory of time reversal in a 2-D "acoustic" medium (i.e., a medium whose deformations are described by the 2-D Helmholtz' equation). Finally (sec. 6), theoretical results are validated by direct application to synthetic and real surface-wave data. The applications presented here are limited to two-dimensions; in future work, we shall explore the resolving power of our method in the vertical direction, combining the results of multiple, Love-and Rayleigh-wave 2-D time-reversal simulations conducted at different frequencies.

Surface waves and the two-dimensional Helmholtz equation

The scalar 2-D Helmholtz' equation can be used to model surface-wave propagation in a laterally smooth, lossless half space. We shall give a simplified proof of this fundamental result by briefly summarizing some parts of earlier studies by [START_REF] Tanimoto | Modelling curved surface wave paths: membrane surface wave synthetics[END_REF] [START_REF] Tromp | Variational principles for surface wave propagation on a laterally heterogeneous Earth-III. Potential representation[END_REF]. Let us start by writing the displacement equation for an elastic, isotropic medium in the frequency (ω) domain [Udías, 1999, eq. (2.60)

], ∂ ∂x j λδ ij ∂u k ∂x k + µ ∂u i ∂x j + ∂u j ∂x i = -ρω 2 u j , (1)
where x 1 , x 2 , x 3 are Cartesian coordinates, with the x 3 -axis perpendicular to Earth's surface (which we assume to be flat) and oriented downward; δ ij is Kronecker's delta, ρ denotes density and λ, µ Lamé's parameters. Repeated indices are implicitly summed over. Following [START_REF] Tanimoto | Modelling curved surface wave paths: membrane surface wave synthetics[END_REF], we assume the Earth to be smooth laterally (horizontal derivatives of ρ, λ, µ, etc. are negligible) but not vertically (x 3 -derivatives of the same parameters are not negligible); eq. ( 1) then takes a slightly different form for i=3 with respect to i=1,2; namely

(λ + µ) ∂ ∂x 1,2 ∂u k ∂x k + µ∇ 2 u 1,2 + ∂µ ∂x 3 ∂u 1,2 ∂x 3 + ∂u 3 ∂x 1,2 = -ρω 2 u 1,2 , (2)
and

(λ + µ) ∂ ∂x 3 ∂u k ∂x k + µ∇ 2 u 3 + 2 ∂µ ∂x 3 ∂u 3 ∂x 3 + ∂λ ∂x 3 ∂u k ∂x k = -ρω 2 u 3 . (3)
The displacement equations (1) or ( 2) and (3) are accompanied by the requirement that no tractions exist on the outer surface of the Earth ("free surface" boundary conditions); for an isotropic elastic medium, the stress tensor

σ ij = λ ∂u k ∂x k δ ij + 2µ ∂u i ∂x j + ∂u j ∂x i , (4)
and the zero-traction requirement at the outer (horizontal) surface is equivalent to requiring that σ 13 =σ 23 =σ 33 =0 when x 3 =0. Displacements and stresses are also usually required to be continuous across all discontinuities.

We next introduce a Rayleigh-wave displacement Ansatz in the frequency domain,

u R = U (x 3 , ω)x 3 φ R (x 1 , x 2 , ω) + V (x 3 , ω)∇ 1 φ R (x 1 , x 2 , ω), (5)
where the unit-vectors x 1 , x 2 , x 3 are parallel to the Cartesian axes, and

∇ 1 = x 1 ∂ ∂x 1 + x 2 ∂ ∂x 2 .
The functions U (x 3 , ω) and V (x 3 , ω) control the dependence of surface-wave amplitude on depth; they do not need to be known explicitly at this stage. The function φ R can be thought of as a "Rayleigh-wave potential". For Love waves,

u L = W (x 3 , ω)(-x 3 × ∇ 1 )φ L (x 1 , x 2 , ω), (6)
with φ L the "Love-wave potential", and W (x 3 , ω) playing the same role as U and V above.

It can be seen by inspection of expressions ( 5) and ( 6) that they indeed describe Rayleighand Love-wave motion, respectively. The functions U , V , W need not be specified at this point, but, if only surface-wave solutions are of interest, it must be required that lim

x 3 -→∞ U (x 3 , ω) = 0 ; lim x 3 -→∞ V (x 3 , ω) = 0 ; lim x 3 -→∞ W (x 3 , ω) = 0. ( 7
)
We next use our surface-wave Ansätze (5) and ( 6), together with the mentioned boundary conditions, to simplify and solve the displacement equations ( 2) and (3).

Love waves

We first substitute u in eqs.

(2), (3) with the expression (6) for u L . It is useful to notice that the x 3 -component of u L is 0, and that u L is divergence-free; as a result, eq. ( 3) is always verifed by u L as given by ( 6), whatever the functions W (x 3 ) and φ L (x 1 , x 2 ). After some algebra, the remaining equations are reduced to

µ ∂ 2 W ∂x 2 3 + ∂µ ∂x 3 ∂W ∂x 3 + ρω 2 W ∂φ L ∂x 2 + µW ∂ 3 φ L ∂x 2 1 ∂x 2 + ∂ 3 φ L ∂x 3 2 = 0, (8) µ ∂ 2 W ∂x 2 3 + ∂µ ∂x 3 ∂W ∂x 3 + ρω 2 W ∂φ L ∂x 1 + µW ∂ 3 φ L ∂x 2 2 ∂x 1 + ∂ 3 φ L ∂x 3 1 = 0. ( 9
)
Remember that φ L is only a function of x 1 , x 2 , while µ = µ(x 3 ), ρ = ρ(x 3 ) and W = W (x 3 ).

If we divide eq. ( 8) by µW ∂φ L ∂x 2 we find

µ ∂ 2 W ∂x 2 3 + ∂µ ∂x 3 ∂W ∂x 3 + ρω 2 W µW = - ∂ 3 φ L ∂x 2 1 ∂x 2 + ∂ 3 φ L ∂x 3 2 ∂φ L ∂x 2 , ( 10
)
which can be solved by separation of variables [e.g. Tromp and Dahlen, 1993, sec. 3] since the right-hand side depends only on x 1 , x 2 , and the left-hand side only on x 3 . This means that we can introduce a constant k L such that

µ ∂ 2 W ∂x 2 3 + ∂µ ∂x 3 ∂W ∂x 3 + ρω 2 W = µk 2 L W (11)
and

∂ 3 φ L ∂x 2 1 ∂x 2 + ∂ 3 φ L ∂x 3 2 = -k 2 L ∂φ L ∂x 2 . ( 12
)
It might be noticed that the "radial" eq. ( 11) is equivalent to equation ( 46) of [START_REF] Takeuchi | Seismic surface waves[END_REF], or eq. ( 7.24) of [START_REF] Aki | Quantitative Seismology[END_REF], even though those treatments are limited to plane waves (which affects φ L but not W ).

Applying the same procedure to eq. ( 9) additionally gives

∂ 3 φ L ∂x 1 ∂x 2 2 + ∂ 3 φ L ∂x 3 1 = -k 2 L ∂φ L ∂x 1 , (13)
and a sufficient condition for φ L to solve both ( 12) and ( 13) is the Helmholtz' equation

∂ 2 φ L ∂x 2 1 + ∂ 2 φ L ∂x 2 2 = -k 2 L φ L . (14)
The boundary conditions can also be simplified when applied to our Love-wave Ansatz:

it follows from ( 4) and ( 6) that, for Love waves, σ 33 =0, and

σ 13 = -σ 23 = µ ∂W ∂x 3 ∂φ L ∂x 1 .
The zero-traction boundary condition at the outer surface thus reduces to ∂W ∂x 3 = 0 at x 3 = 0. (15)

Love-wave radial equation

Several different approaches to the (semi-analytical or numerical) solution of the "radial" equation ( 11) are reviewed in sections 7.1 and 7.2 of [START_REF] Aki | Quantitative Seismology[END_REF], starting with a simple one-layer-over-half-space model and then generalizing to the cases of an arbitrary number of layers, and of continuous velocity and density profiles. We need not repeat here the detailed treatment of [START_REF] Aki | Quantitative Seismology[END_REF], but it is useful to point out some of its essential implications.

Equation ( 11) is a second-order ordinary differential equation, whose general solution thus contains two arbitrary constants. Two boundary conditions must be taken into account: eqs.

(7) and (15). These two equations allow in principle to determine both arbitrary constants.

If the Earth is modeled as a set of one or more uniform, horizontal layers overlying a half space, then within each layer i we have ∂µ ∂x 3 = 0, and eq. ( 11) is simplified to the Helmholtz' equation

∂ 2 W ∂x 2 3 + (ρ i ω 2 -µ i k 2 L )W = 0, ( 16
)
where ρ i and µ i denote the (constant) values of density and rigidity within layer i, respectively.

Each layer adds one second-order equation, and therefore two arbitrary constants to the problem, but also one interface with the associated two continuity conditions (on W and

The parameters ω and k L , however, have not been specified, and, as a consequence, one cannot simply identify a unique solution for W to be substituted into the Ansatz (6).

According to [START_REF] Aki | Quantitative Seismology[END_REF], this problem is solved in general as follows: (i) a numerical value ω 0 is assigned to ω; (ii) a numerical, "trial" value is likewise assigned to k L ;

(iii) the selected numerical values ω 0 and k L are substituted into eq. ( 11) which can then be integrated numerically, or via the "propagator matrix" method [e.g. Aki and Richards, 2002, sec. 7.2.2], starting with W =0 at large depth x 3 ; (iv) it is verified whether condition ( 15) is met at x 3 =0; (v) if this condition is not met, eq. ( 11) is integrated again, with the same ω 0 but a different trial value for k L ; (vi) if, instead, the condition ( 15) is met, the whole process is repeated for a new value ω 0 , until the frequency range of interest is entirely covered.

It is found that a discrete set of one or more (depending on ω 0 ) values of k L for which the free-surface boundary condition is met can be determined [e.g. Aki and Richards, 2002, figures 7.2, 7.3]. These values are dubbed "eigenvalues" in analogy with free-oscillation theory, and each corresponds to a different solution, or "mode," for W . If more than one eigenvalue exist at a given frequency, the mode corresponding to the largest k L eigenvalue is referred to as "fundamental mode," followed by "higher modes" ("overtones").

Helmholtz' equation for the Love-wave potential φ L

The parameters k L and ω in eq. ( 14) must be substituted with one of the eigenvalues of k L , and with the corresponding value ω 0 , respectively, before this equation is solved for φ L .

Substitution of φ L (x 1 , x 2 , ω 0 ) and of the corresponding W (x 3 , ω 0 ) into expression (6) yields a monochromatic Love-wave solution. The process can be iterated at each frequency ω 0 for which the eigenvalues k L and eigenfunctions W have been determined as described in sec. 2.1.1.

Notice that, for a monochromatic wave of frequency ω 0 , eq. ( 14) coincides with the 2-D wave equation with wavespeed

ω 0 k L . The curve k L = k L (ω 0 ) is thus the "dispersion curve"
describing how surface-wave phase velocity depends on frequency.

It is easy to show that a monochromatic plane wave φ L would solve eq. ( 14), and in fact most seismology textbooks replace φ L (and φ R ) with plane-wave formulae in the surface-

wave Ansätze [e.g. [START_REF] Aki | Quantitative Seismology[END_REF]. In view of the applications to be discussed here, however, circular (cylindrical) surface waves are more relevant. This case can be described, in our formulation, starting with the known solution G 2D to the Green's problem associated with equation ( 14), obtained, e.g., in Appendix E of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF];

G 2D (x 1 , x 2 , ω)
is clearly not a monochromatic wave, but the response of the medium to a monochromatic point source can be obtained, according to eq. (E34) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], by timedomain convolution or frequency-domain multiplication of G 2D (x 1 , x 2 , ω) with a sinusoidal signal δ(ω -ω 0 ).

Rayleigh waves

In analogy with section 2.1, we next substitute u in eqs. ( 2), (3) with the expression (5) for u R . This results, after some algebra, in the system of equations

µ ∂ 2 V ∂x 2 3 + ∂µ ∂x 3 U + ∂V ∂x 3 + (λ + µ) ∂U ∂x 3 + ρω 2 V ∂φ R ∂x 1 +(λ+2µ)V ∂ 3 φ R ∂x 3 1 + ∂ 3 φ R ∂x 2 2 ∂x 1 = 0, ( 17
) µ ∂ 2 V ∂x 2 3 + ∂µ ∂x 3 U + ∂V ∂x 3 + (λ + µ) ∂U ∂x 3 + ρω 2 V ∂φ R ∂x 2 +(λ+2µ)V ∂ 3 φ R ∂x 3 2 + ∂ 3 φ R ∂x 2 1 ∂x 2 = 0, ( 18
) (λ + 2µ) ∂ 2 U ∂x 2 3 + 2 ∂µ ∂x 3 ∂U ∂x 3 + ∂λ ∂x 3 ∂U ∂x 3 + ρω 2 U φ R + (λ + µ) ∂V ∂x 3 + ∂λ ∂x 3 V + µU ∂ 2 φ R ∂x 2 1 + ∂ 2 φ R ∂x 2 2 = 0, (19)
which, again, can be solved by the method of separation of variables. After dividing it by (λ + µ) ∂V ∂x 3 + ∂λ ∂x 3 V + µU , eq. ( 19) can be separated into

(λ + 2µ) ∂ 2 U ∂x 2 3 + 2 ∂µ ∂x 3 ∂U ∂x 3 + ∂λ ∂x 3 ∂U ∂x 3 + ρω 2 U = k 2 R (λ + µ) ∂V ∂x 3 + ∂λ ∂x 3 V + µU (20)
and the Helmholtz' equation

∂ 2 φ R ∂x 2 1 + ∂ 2 φ R ∂x 2 2 = -k 2 R φ R , (21)
where k R is, at this point, an arbitrary constant. If one then substitues eq. ( 21) into ( 17) and ( 18), it becomes apparent that a sufficient condition for both of them to be solved is given by

(λ + µ) ∂U ∂x 3 + µ ∂ 2 V ∂x 2 3 + ∂µ ∂x 3 U + ∂V ∂x 3 + ρω 2 V -k 2 R (λ + 2µ)V = 0. ( 22
)
The "radial" eqs. ( 20) and ( 22) form a linear system of second-order ordinary differential equations that can be solved to determine U and V . Since two equations and two unknown functions are now involved, the solution is more cumbersome, but qualitatively similar to the Love-wave case of sec. 2.1.1. Again, as shown by [START_REF] Takeuchi | Seismic surface waves[END_REF] and [START_REF] Aki | Quantitative Seismology[END_REF] for the plane-wave case, a set of Rayleigh-wave "modes" can be found by numerical integration: each mode is defined by a frequency ω 0 and a value of k R for which (20) and ( 22) are solved, and the boundary conditions met. The definitions of fundamental mode and overtone given in sec. 2.1.1 naturally holds also for Rayleigh waves.

The discussion of sec. 2.1.2 on the Love-wave potential φ L also applies to the Rayleighwave potential φ R , which is controlled by the Helmholtz' equation ( 21); in analogy with sec. 2.1.2, k R can be interpreted as the ratio between the frequency and phase velocity of the corresponding Rayleigh-wave mode.

Reciprocity theorem in two dimensions

Consider the non-homogeneous 2-D Helmholtz' equation

∇ 2 1 p(x 1 , x 2 , ω) + ω 2 c 2 p(x 1 , x 2 , ω) = -iωq(x 1 , x 2 , ω), ( 23
)
where p could represent the displacement of a stretched membrane (whose density and tension determine the parameter c), and the forcing term -iωq a pressure exerted on the membrane per unit of surface density [e.g., Kinsler et al., 1999, secs. 4.2 and 4.8]. Here and in the following we denote by f (ω) the Fourier transform of a generic function f (t), and by i the imaginary unit. The following mathematical treatment makes it convenient to denote forcing as -iωq(x 1 , x 2 , ω).

Let us define a vector v = -1 iω ∇ 1 p, such that

∇ 1 p + iωv = 0. ( 24
)
Substituting eq. ( 24) into ( 23), we then find

∇ 1 • v + iω c 2 p -q = 0. ( 25
)
The following treatment follows closely that of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], who summarized earlier results by, e.g., [START_REF] Wapenaar | Green's function representations for seismic interferometry[END_REF] and [START_REF] Snieder | Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves[END_REF], limited to three-dimensional space. Let us consider a surface S bounded by the closed curve ∂S. (∂S is just an arbitrary closed curve within a 2-D medium, and generally does not represent a physical boundary.) Let q A (x 1 , x 2 , ω), p A (x 1 , x 2 , ω) and v A (x 1 , x 2 , ω) denote a possible combination of the fields q, p and v co-existing at (x 1 , x 2 ) in S and ∂S. A different forcing q B would give rise, through eqs. ( 24) and ( 25), to a different "state" B, defined by p B (x 1 , x 2 , ω)

and v B (x 1 , x 2 , ω).
A useful relationship between the states A and B, known as "reciprocity theorem", is obtained by combining eqs. ( 24) and (25) as follows,

S d 2 x [(24) A • v * B + (24) * B • v A + (25) A p * B + (25) * B p A ] = 0, ( 26
)
where x=(x 1 , x 2 ), d 2 x = dx 1 dx 2 , and * denotes complex conjugation. (24) A is short for the expression one obtains after substituting p = p A (x, ω) and v = v A (x, ω) into the left-hand side of eq. ( 24), etc. Namely, (24

) A • v * B = ∇ 1 p A • v * B + iωv A • v * B ( 27
) (24) * B • v A = ∇ 1 p * B • v A -iωv * B • v A (28) (25) A p * B = ∇ 1 • v A p * B + iω c 2 p A p * B -q A p * B ( 29
) (25) * B p A = ∇ 1 • v * B p A - iω c 2 p * B p A -q * B p A . ( 30
)
After substituting expressions ( 27) through (30) into eq. ( 26), the latter simplifies to

S d 2 x (∇ 1 p A • v * B + ∇ 1 p * B • v A + ∇ 1 • v A p * B + ∇ 1 • v * B p A ) = S d 2 x (q A p * B + q * B p A ) . ( 31
)
The integrand at the left-hand side of (31) can be further simplified via the relationship

∇ 1 • (p A v * B ) = ∇ 1 p A • v * B + ∇ 1 • v * B p A (
which naturally holds also if A and B are swapped).

We next apply the 2-D version of the divergence theorem to the resulting expression, and eq.

(31) collapses to

∂S dx (p A v * B + p * B v A ) • n = S d 2 x (q A p * B + q * B p A ) , ( 32
)
where n is a unit vector everywhere perpendicular to ∂S. For instance, "Green's identity" (4.22) of [START_REF] Baker | The Mathematical Theory of Huygens' Principle[END_REF], or the "reciprocity theorem of the correlation type", eq.

(5) of [START_REF] Wapenaar | Green's function representations for seismic interferometry[END_REF], are three-dimensional versions of eq. ( 32).

Application of the reciprocity theorem to impulsive point sources: exact equations

Let us next consider the states A and B resulting from the impulsive forcing terms q A = δ(x -x A ) and q B = δ(x -x B ), respectively, with x A , x B two arbitrary locations on S. It

follows that p A = G 2D (x, x A , ω) and p B = G 2D (x, x B , ω), with G 2D the Green's function
corresponding to a 2-D membrane excited by a nonzero right-hand side in eq. ( 23), and eq.

(

) then implies that v A,B = -1 iω ∇ 1 G 2D (x, x A,B , ω). 24
G 2D is the solution of the non-homogeneous eq. ( 23) with q = δ(x -x A,B ). Based on eq.

(E34) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF],

G 2D (x, x A,B , ω) = R 2 d 2 x G 2D (x, x , ω)(-iω)δ(x -x A,B ) = -iωG 2D (x, x A,B , ω), (33)
where G 2D (x, x , ω) is the Green's function associated with a nonzero initial velocity at x , derived explicitly e.g. by [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF]. To translate the time-domain formula of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF] into frequency domain, it is useful to notice that eq. (E34)

of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF] involves the time-domain convolution of G 2D with the nonhomogeneous term (forcing term) of the wave equation, and to remember that a convolution in the time-domain maps to a product in the frequency domain.

Replacing p A,B and v A,B in eq. ( 32) with their expressions in terms of G 2D , and q A,B with a Dirac delta,

G * 2D (x A , x B , ω) -G 2D (x B , x A , ω) = ∂S dx G * 2D (x , x B , ω)∇ 1 G 2D (x , x A , ω) -G 2D (x , x A , ω)∇ 1 G * 2D (x , x B , ω) • n. (34)
Eq. ( 34) can be thought of as the 2-D version of eq. ( 19) in [START_REF] Wapenaar | Green's function representations for seismic interferometry[END_REF] or eq. ( 96) in [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF].

The above treatment holds if x A and x B are not within S; in that case, q A,B are zero within S. The integral at the right-hand side of eq. ( 32) is therefore zero, and so is, as a result, the left-hand side of (34). It follows that the integral at the right-hand side of ( 34 3.2 Application of the reciprocity theorem to impulsive point sources: farfield/high-frequency approximation Equation ( 34) can be simplified by the "far-field" approximation, which requires that the locations x A and x B be separated from one another and from δS by at least a few wavelengths.

We additionally require that xx A ≈ xx B for any point x on ∂S (i.e., x A and x B are both very far from ∂S). The Green's function G 2D can be replaced by its far-field approximation, which reads

G 2D (x, y, ω) ≈ 1 4iπc 3/2 e -i ω|x-y| c -π 4 ω|x -y| (35)
[e.g., Boschi and Weemstra, 2015, eq. ( E17)]. We next take advantage of this approximation to find a simple expression for ∇ 1 G 2D . Let us consider for example ∇ 1 G 2D (x, x A , ω) and call

r = |x -x A |. Then, ∇ 1 G 2D (x, x A , ω) ≈ 1 4iπc 3/2 x 1 ∂ ∂x 1 + x 2 ∂ ∂x 2 e -i( ωr c -π 4 ) √ ωr = 1 4iπc 3/2 x 1 ∂r ∂x 1 + x 2 ∂r ∂x 2 ∂ ∂r e -i( ωr c -π 4 ) √ ωr = 1 4iπc 3/2 iω c + 1 2r e -i( ωr c -π 4 ) √ ωr ∇ 1 r = G 2D (x, x A , ω) iω c + 1 2r ∇ 1 r. (36)
In the far-field approximation, r is large and r -1 is much larger than r -2 : the second term inside square brackets in eq. ( 36) can be neglected. If one takes the origin, e.g., at x A , the condition xx A ≈ xx B implies that both xx A and xx B can be replaced by x, and

∇ 1 r ≈ x |x| .
We are left with

∇ 1 G 2D (x, x A , ω) ≈ iω c G 2D (x, x A , ω) x |x| , (37)
which we can finally substitute into eq. ( 34), to find

G * 2D (x A , x B , ω) -G 2D (x B , x A , ω) ≈ iω c ∂S dx G 2D (x , x A , ω)G * 2D (x , x B , ω) + G * 2D (x , x B , ω)G 2D (x , x A , ω) x |x | • n. ( 38
)
Remember that the closed curve ∂S does not correspond to a physical boundary. We choose it to be circular (we shall see in the following that this assumption does not affect the relevant physical interpretations of our results), so that x |x| = n on ∂S. Eq. ( 38) collapses to

G * 2D (x A , x B , ω) -G 2D (x B , x A , ω) ≈ 2iω c ∂S dx G 2D (x , x A , ω)G * 2D (x , x B , ω), (39)
which is the 2-D counterpart of eq. ( 102) in [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF]. (It is also consistent with eq. ( 65) of the same study, valid for a 2-D medium, derived via the stationary-phase approximation and setting source density to 1.)

4 Implications for surface waves: diffuse-field interferometry, time reversal

We know from sec. 2 that the Rayleigh-and Love-wave potentials φ R , φ L , just like the "membrane-wave" field p, obey the Helmholtz' equation ( 23). It follows that eqs. ( 34) and

(39) continue to be valid if p is replaced by potentials φ R or φ L , and if c is the Rayleighor Love-wave phase velocity at that frequency. We also know that the vertical displacement associated with a Rayleigh wave is proportional to φ R and thus obeys (23) exactly at the frequency ω [e.g. [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF]; slightly more complicated relations exist between

Love-wave displacement (and the horizontal component of Rayleigh-wave displacement) and the Love-wave (Rayleigh-wave) potential, which are given e.g. by [START_REF] Kaestle | Two-receiver measurements of phase velocity: Cross-validation of ambient-noise and earthquake-based observations[END_REF]. In summary, the results of sec. 3 can be applied to the propagation of seismic surface waves, which will be the focus of the remainder of this study.

Eqs. ( 34) and its approximate version (39) describe the physics underlying both ambientnoise interferometry and acoustic/seismic time reversal. Analogies between these two techniques were first discussed by [START_REF] Derode | How to estimate the Green's function of a heterogeneous medium between two passive sensors? Application to acoustic waves[END_REF].

Analogy with diffuse-field interferometry

In the context of diffuse-field interferometry, the far-field eq. ( 39) is invoked more often than its exact counterpart (34). The points x A , x B in eq. ( 39) are taken to represent the locations of two receivers, while the points on δS are thought of as point sources. The right-hand side of (39) is the cross-correlation of the signal recorded at receiver x A with that recorded at receiver x B , averaged (integrated, "stacked"...) over all sources. It is usually assumed that sources are approximately distributed along a closed curve surrounding the receivers in the far field. (If that is the case, it has also been shown that the cross-correlation of signals generated by different sources that act simultaneously will tend to cancel out; see [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF] for a more detailed discussion.) Eq. ( 39) then implies that the receiver-receiver crosscorrelation at its right-hand side coincides approximately with the imaginary part of the frequency-domain Green's function G 2D at the left-hand side. Since G 2D is real in the timedomain, and nonzero only at positive times, G 2D (x A , x B , t) is determined without ambiguity by the imaginary part of its Fourier transform G 2D (x A , x B , ω) [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF].

It follows that the surface-wave Green's function between two locations can be reconstructed from the cross-correlation of a diffuse surface-wave field recorded at those locations.

Recall now that, at the beginning of sec. 3.1, the assumption has been made that sources be impulsive. In practice, this amounts to selecting q(x, ω) = δ(x -x A,B ), i.e., in the time domain, q(x, t) = δ(x -x A,B )δ(t). Let us next consider the case of arbitrary, unspecified time-dependence h(t) of the source signal, i.e. q(x, t) = δ(x -x A,B )h(t). Eq. ( 33) here was obtained from eq. (E34) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], replacing the generic source signal there with the right-hand side of the non-homogeneous Helmholtz' eq. ( 23). The solution p corresponding to an arbitrary source signal h(t) is therefore obtained by updating the righthand side of eq. ( 23), which in the frequency domain now reads -iωδ(x -x A,B )h(ω). After substituting this into eq. (E34) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], we find

p(x, x A,B , ω) = -iωG 2D (x, x A,B , ω)h(ω), (40)
which replaces our eq. ( 33).

Substituting (40) and q(x, ω) = δ(x -x A,B )h(ω) into eqs. ( 24) and ( 32), we find that introducing the time-dependence h(t) of the source boils down to multiplying both sides of eq. ( 34), and therefore (39), by the squared Fourier spectrum |h(ω)| 2 .

In ambient-noise interferometry, this means that if all noise sources had the same spectrum then the cross-correlation of recorded ambient signal would also exhibit that spectrum (squared): consequently, we would not be reconstructing the Green's function but rather its time-domain convolution with the source-related term |h(ω)| 2 . Indeed, it is well known that the spectrum of seismic ambient-noise cross correlation is dominated by peaks that correspond to the spectrum of oceanic microseisms [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF]

Surface-wave time reversal

If eq. ( 34) is to be used as an illustration of time-reversal acoustics [e.g., [START_REF] Fink | Time-reversed acoustics[END_REF], x B should be thought of as the location of a source; G 2D (x , x B , ω) is the Fourier-transform of the signal generated at x B and recorded by a far-away receiver at x ; its complex-conjugate

G * 2D (x , x B , ω)
is the Fourier transform of the same signal, reversed in time. Imagine that the time-reversed signal be then emitted from x and recorded at another point x A : this amounts to convolving (in the frequency domain, multiplying) the time-reversed signal with the Green's function G 2D (x A , x , ω). Eq. ( 39) then shows that by repeating time reversal and propagation ("backward in time") for all points x on ∂S and summing all the resulting traces at x B , the imaginary part of the Green's function between x B and x A is obtained.

If ∂S is in the near field of x A , x B , the approximate eq. ( 39) should be replaced by (34), which is the 2-D version of eq. (3) in [START_REF] Fink | Time-reversal acoustics in complex environments[END_REF]. In practice, this means that to reconstruct the Green's function between x A and x B one needs to (i) time-reverse (in the frequency domain, take the complex-conjugate of) the signal G 2D (x , x B , ω) emitted by x B and recorded at x ; (ii) take the spatial derivative of the time-reversed signal in the n direction at

x , i.e. n • ∇ 1 G * 2D (x , x B , ω); (iii) convolve the time-reversed signal G * 2D (x , x B , ω) with the dipole response (see appendix) n • ∇ 1 G 2D (x A ,
x , ω) between x A and x ; (iv) convolve its spatial-derivative with the impulse response between the same two points; (v) sum the two signals obtained at (iii) and (iv). In other words, rather than simply backward-propagating the signal recorded at receivers on ∂S, as in the far-field case, we must backward-propagate the sum of a dipole and a monopole source, to which the initial signal itself and its spatial derivative are "fed", respectively.

The backward-propagated signal so obtained coincides, approximately (if eq. ( 39) is implemented) or exactly (eq. ( 34)), with the difference

G * 2D (x A , x B , ω) -G 2D (x B , x A , ω) = -2i [G 2D (x A , x B , ω)] . (41)
To understand the physical meaning of this expression, let us take its inverse Fourier transform (F -1 ). It follows from eqs. (B6) and (B9) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF] that

F -1 {-2i [G 2D (x A , x B , ω)]} = G 2D (x A , x B , t) -G 2D (x A , x B , -t), (42)
similar, e.g., to eq. ( 6) of [START_REF] Fink | Time-reversal acoustics in complex environments[END_REF]. Consider an arbitrary observation point x A within ∂S, and recall that G 2D is nonzero only at positive time. As t grows from -∞ to 0, only the second term at the right-hand side of ( 42) is nonzero, which means that x A records a time-reversed Green's function. In space, a time-reversed impulsive circular wave converging towards the original source location x B is observed. As t-→0, the value of x A for which the field is maximum approaches x B , where the backward-propagating circular wave eventually "focuses." As t grows from 0 to ∞, only the first term at the right-hand side of ( 42) is nonzero, and x A records a regular Green's function with inverted sign. That is to say, another circular wave is emitted from x B after focusing.

In the words of [START_REF] Fink | Time-reversal acoustics in complex environments[END_REF], "if we were able to create a film of the propagation of the acoustic field during" propagation of the signal from the original source at x B to receivers on ∂S, "the final result could be interpreted as a projection of this film in the reverse order, immediately followed by a reprojection in the initial order." [START_REF] Fink | Time-reversal acoustics in complex environments[END_REF] notes that acoustic time reversal, as described here, does not involve the "time reversal of the source," and in "an ideal time-reversed experiment, the initial active source (that injects some energy into the system) must be replaced by a sink (the time reversal of a source)," i.e., "a device that absorbs all arriving energy without reflecting it."

Result ( 41) is limited to impulsive signals. If the signal emitted at x B is an arbitrary function of time, h(t), the signal recorded at each receiver location x is the convolution

G 2D (x , x B , ω)h(ω).
Accordingly, let us replace G * 2D (x , x B , ω) at the right-hand side of eq. ( 34) with the convolution G * 2D (x , x B , ω)h * (ω). Since the function h does not depend on any other variable but t (or ω in the frequency domain), it can be pulled out of the x -integral; it then follows from eq. ( 34) itself that

h * (ω) [G * 2D (x A , x B , ω) -G 2D (x B , x A , ω)] = ∂S dx h(ω)G 2D (x , x B , ω) * ∇ 1 G 2D (x , x A , ω) -G 2D (x , x A , ω)∇ 1 h(ω)G 2D (x , x B , ω) * • n. ( 43
)
If one denotes s(x , x B , ω) = h(ω)G 2D (x , x B , ω), eq. ( 43) takes the more compact form

h * (ω) [G * 2D (x A , x B , ω) -G 2D (x B , x A , ω)] = ∂S dx s * (x , x B , ω)∇ 1 G 2D (x , x A , ω) -G 2D (x , x A , ω)∇ 1 s * (x , x B , ω) • n. ( 44
)
Alternatively, the far-field approximation (37) can be applied to (43), which, in analogy with sec. 3.1, collapses to

h * (ω) [G * 2D (x A , x B , ω) -G 2D (x B , x A , ω)] ≈ 2iω c ∂S dx s * (x , x B , ω)G 2D (x , x A , ω) . ( 45
)
Eqs. ( 44) and ( 45) stipulate that the same results obtained above for impulsive signals also apply to arbitrary signals h(t), except that in this case the backward-propagating Green's function is convoluted with the time-reversed signal, h * (ω) or h(-t), itself. Importantly, the backward-propagating wave field focuses, again, on the source location x B .

The so-called membrane-wave approach is based on the horizontal/radial decoupling of the equation of motion illustrated in sec. 2, where it is shown that the membrane eq. ( 23) holds for both the Love-and Rayleigh-wave potentials φ L and φ R . In sec. 3, some properties of the solution of ( 23), that naturally apply to both φ R and φ L , are derived. Their most important implication in the context of our study is explained in sec. 4: the theory of acoustic time reversal as developed by, e.g., [START_REF] Fink | Time-reversal acoustics in complex environments[END_REF] holds on a flat membrane, and, as a result, the timereversed potentials φ L , φ R can be otbained from eqs. ( 44) or ( 45), i.e., by time-reversing and backward propagating the potentials associated with the recorded waveforms. Surface-wave time reversal then consists of (i) extracting φ L and φ R from the data, for a broad and dense set of surface-wave fundamental and higher modes; (ii) determining radial eigenfunctions (U and V , or W ) for each mode; (iii) backward propagating φ L and φ R for each mode; (iv) combining potentials with radial eigenfunctions at all available frequencies, via eqs. ( 5) and ( 6), to find the time-reversed displacements u R and u L .

An important limitation of this procedure, as discussed in some detail in sec. 4.2, is that the time-reversed wave field necessarily includes an impulse propagating away from the reconstructed source location after focusing. This is not a problem for point sources (or of less-than-wavelength spatial extent, as in this study), but the time-reversed wave field at each point of a finite-extent source will include a non-physical contribution that cannot easily be subtracted from it, and that pollutes images of seismic slip. It should be noted that backprojection methods suffer from the same problem, although this is rarely (if ever) discussed.

This issue will have to be addressed in future work. One possible strategy would be to subdivide the source-imaging process into two steps. First, time reversal could be interrupted before focusing occurs: this way, the surface-wave field in the immediate vicinity of the source could be reconstructed. In a second step, the reconstructed near-field displacement could be treated as data in a classic linear inverse problem, based on the representation theorem [e.g. [START_REF] Ide | Slip inversion[END_REF]: the unknown being slip on the fault. The accuracy of near-field displacement as reconstructed by time reversal would significantly reduce nonuniqueness.

Only monochromatic, fundamental-mode Rayleigh-wave propagation is implemented here.

At each frequency of interest ω, propagation of the corresponding sinusoidal Rayleigh wave is modeled in the time domain. It is apparent from eq. ( 5) that, at frequency ω, φ R is directly proportional to the vertical component of displacement, narrow-band-pass-filtered around ω; i.e., before time reversal, φ R (ω) can be obtained by the vertical component of the displacement by simply multiplying it by 1/U (ω). The implementation of time reversal is exactly the same for Love waves (except that membrane-wave propagation of the Lovewave potential must naturally be modeled in a Love-wave phase velocity map); the Lovewave potential φ L , however, needs to be extracted from the transverse component of cross correlations, which will require some more subtle data processing to be addressed in future work.

Accordingly, we do not yet reconstruct time-reversed displacements from time-reversed potentials. This requires that the eigenfunctions U , V and W be computed for a selected reference model. Because the crust/lithosphere depth range (i.e., the depth range of interest to surface-wave propagation) is characterized by large lateral heterogeneity, it is likely that a 3-D reference model will need to be employed, through the implementation of "local" radial eigenfunctions [e.g. [START_REF] Boschi | New images of the Earth's upper mantle from mea-surements of surface wave phase velocity anomalies[END_REF]. Studying the focusing of the Rayleigh-wave vertical component at various frequencies is, however, sufficient to verify the feasibility of our approach, which is the main goal of the present study. In the following, we model the propagation of time-reversed surface-wave potentials via two different approaches: ray theory and the spectral-element method.

In the ray-theory case, the value of G 2D for any given source and receiver position is determined approximately by tracing the ray between source and receiver, computing the propagation time along such ray path, and shifting by such time the signal prescribed at a source. Rays are traced by means of the algorithm described by [START_REF] Fang | Direct inversion of surface wave dispersion for three-dimensional shallow crustal structure based on ray tracing: methodology and application[END_REF]. Geometrical spreading is accounted for approximately by simply multiplying the signal by the inverse squared root of the source-receiver distance, according, e.g., to eq. (E17) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF].

In the spectral-element case, following [START_REF] Tape | Finite-frequency tomography using adjoint methods. Methodology and examples using membrane surface waves[END_REF], SPECFEM2D [START_REF] Komatitsch | Introduction to the spectral element method for threedimensional seismic wave propagation[END_REF] is used to simulate the propagation, on a stretched, flat membrane, of a displacement perpendicular to the unperturbed membrane surface. Displacement is generated by prescribing a point force/acoustic pressure (rather than an initial displacement as in the ray-theory case), which implies, importantly, that our comparison between ray-theory and SPECFEM2D results is only qualitative. Additionally, to model wave propagation via SPECFEM2D, we need to project our spherical-Earth phase-velocity map onto a flat surface. This is done via a transverse Mercator projection centered at 12 • E, 46 • N. Errors are introduced near the corners of the region of interest, that will (slightly) alter modeled waveforms and might reduce the quality of focusing: the flat-membrane approach is adequate to the feasibility study presented here, but curved membranes will have to be implemented for future applications.

Validation

We test both ray-theory and spectral-element methods on synthetic membrane-wave data, on ambient-noise vertical-component cross correlations (which are theoretically equivalent to recordings of Rayleigh-wave impulse responses) and on vertical-component recordings of a 5.6-magnitude earthquake. To make sure that we can rely on robust, high-resolution Rayleighwave phase-velocity maps and a dense station coverage, we select Northern Italy, including most of the Alpine mountain range, as our study region. This area is characterized by complex tectonics, and at this scale surface waves are difficult to identify as they are hidden in the body-wave coda: if we can validate our theory in such an unfavorable situation, we can then expect that it will hold also at teleseismic scales. Furthermore, by limiting the experiments presented here to a relatively small region, we reduce the associated computational costs.

Earthquake data were downloaded from EIDA (http://www.orfeus-eu.org/data/eida/)

and from all permanent broadband stations that recorded the earthquake within the region via the INGV data center. Time-domain cross correlations were computed as described by [START_REF] Molinari | Imaging 3D crustal structure in the presence of deep sedimentary basins: Po-Plain (Italy) case study[END_REF].

As a general rule, computational costs are much reduced with respect to typical 3-D wavepropagation modeling applications in seismology. A time-reversal simulation, such as the ones shown in the following, involves one single run of SPECFEM2D with multiple sources (one per station), which requires about two hours on a single CPU. Ray-theory based simulations are much cheaper: a time-reversal simulation can be completed in less than two minutes on similar hardware.

Synthetic tests

Theoretical traces associated with a selected point-source location and a realistic station distribution in the region of interest are obtained via ray theory and SPECFEM2D. The source signal h(t) is a Ricker wavelet as implemented in SPECFEM2D, Butterworth-filtered between 6 and 26 s. Membrane waves are propagated through the 16s Rayleigh-wave phasevelocity map of [START_REF] Kaestle | Surface-wave tomography of the alps using ambient-noise and earthquake phasevelocity measurements[END_REF], shown here in Fig. 1b. While only one particular surface-wave mode is implemented for this synthetic test, it is understood that the exact same procedure can be applied in the calculation of other Rayleigh-and Love-wave fundamental modes and overtones. No random noise is added to the synthetic signal.

Because the station distribution is nonuniform, the curve ∂S and, as a consequence, the vector n in eq. ( 34) are not uniquely defined. We avoid this difficulty by replacing eq.

(34) with its far-field approximation (39), which can be implemented without specifying n.

Preliminary experiments show that, despite the small size of the study region, the location of the backward propagating wave field's focus is not visibly affected by this approximation. We plan to find ways to implement (34) exactly in future work, but we believe that the simplified approach employed here is adequate to the scope of this article. We accordingly time-reverse the traces, and propagate them backward in time, essentially implementing the right-hand side of eq. ( 45). Again, waves are propagated through the map of Fig. 1b.

We obtain a pair of animations, one based on ray theory and the other on SPECFEM2D.

Samples of both are shown in Fig. 2. Fig. 3 shows the prescribed and reconstructed signal at the known source location, again for both methods. While the backward propagating wave fields differ because of the mentioned physical differences in the implementation (excitation by initial displacement vs. point force, curved membrane vs. Mercator projection), it appears from Fig. 2c,d that the maximum amplitude with respect to time and position in both simulations corresponds to the known, initial source location. This confirms the validity of surface-wave time reversal as a tool to localize/image a seismic source, despite the severe nonuniformity in receiver distribution. The maximum is less pronounced in the spectral-element simulation, resulting in normalized amplitudes throughout the simulation to be larger than ray-theory-based amplitudes. After focusing, in the absence of an "acoustic sink" (sec. 4.2), a non-physical wave field propagates away from the source. There are no major differences in the quality of focusing achieved by ray-theory vs. SPECFEM2D backward-propagation. In panel a, the difference between forcing term and reconstructed signal is explained by the fact that, in SPECFEM2D, displacement is initiated by prescribing a point force, rather than a displacement as in the ray-theory case.

Ambient-noise cross correlations

Cross-correlations of ambient signal form a perfectly suited data set to validate a sourcelocalization method: each cross-correlation is an approximation for the corresponding receiverreceiver Green's function, and the location of both receivers is naturally well known. We select station LSD.GU (Fig. 5) as our virtual "test" source, and time-reverse and backward propagate ambient-noise based Green's functions associated with it. (Noise cross correlations will be described in a separate study [START_REF] Molinari | Imaging 3D crustal structure in the presence of deep sedimentary basins: Po-Plain (Italy) case study[END_REF].) This amounts to implementing the right-hand side of eq. ( 39) via our two algorithms. We first Butterworth-filter verticalcomponent cross correlations around 16 s (low and high corner frequencies corresponding to periods of 26 and 6 s, respectively), and, as in sec. 6.1, propagate time-reversed signal through the Rayleigh-wave phase-velocity map of Fig. 1b. The results of this exercise are shown in Fig. 4 and Fig. 5. Again, despite the poor azimuthal station coverage in this example, the time-reversed wave field focuses quite precisely on the virtual source, in both the ray-theory and SPECFEM2D implementations. In both cases, the maxima of the timereversed wave field at the known source location is correctly achieved at t=0. Similar to sec. 6.1, non-physical signal naturally emerges after focusing.

If station coverage were uniform and the noise-based Green's functions perfectly reconstructed, the time-reversed signal at LSD.GU (Fig. 5) should closely approximate an impulse, which is not the case. We have seen, however, from the results of sec. 6.1 and in particular Fig. 3, that the source time function can be reconstructed well even when the station coverage is poor. We infer that artifacts in the trace of Fig. 5 result from inaccuracies in the reconstructed Green's function. This is not surprising because, while the phase of Green's functions is reconstructed well by ambient-noise cross correlation, their amplitude probably is not [e.g. [START_REF] Ekström | Determination of surface-wave phase velocities across USArray from noise and Aki's spectral formulation[END_REF].

We next iterate the ray-theory procedure for the 4-to-10 s and 20-to-30 s period bands, and show the results in Figs. 6 and7. Membrane-wave propagation is modeled using phasevelocity maps at 6 and 25 s periods (Fig. 1a,c), respectively, again from [START_REF] Kaestle | Surface-wave tomography of the alps using ambient-noise and earthquake phasevelocity measurements[END_REF].

The quality of focusing is comparable to the intermediate-period case (Fig. 4), and can be considered high, in view of the nonuniformity of station distribution. This result confirms that our algorithm can be applied to a variety of surface-wave modes, and, as far as reliable phase-velocity and Green's function estimates are available, is fairly independent of period, and of the width of the passband.

Recordings of the Emilia earthquake of May 29, 2012

We apply our ray-theory-and SPECFEM2D-based algorithms to vertical-component recordings of the magnitude Mw=5.6 (MI=5.8) Emilia earthquake of May 29, 2012, 7:00:03 AM.

These data, discussed in detail by [START_REF] Molinari | Development and Testing of a 3D Seismic Velocity Model of the Po Plain Sedimentary Basin, Italy[END_REF], are shown here in Fig. 8. Traces are filtered around 16 s, the same way as in sec. 6.2, before time-reversal and backward propagation; propagation is modeled according to the 16s Rayleigh-wave phase velocity map of Fig. 1b. Results are summarized in Figs. 9 and 10. Early time-steps (e.g., Fig. 9a,b) are characterized by the emergence of time-reversed late arrivals, that we believe to be associated with reverberations, e.g. at the sharp boundaries between Po plain and surrounding mountain ranges. This signal does not focus sharply anywhere on our membrane, and can accordingly be neglected in this context. Direct-arrival surface waves, on the other hand, do focus at the known epicenter location in both our implementations (Fig. 9c,d). Similar to secs. 6.1 and 6.2, non-physical signal again emerges after focusing (Fig. 9e,f).

Fig. 10 shows that the maximum amplitude of the reconstructed vertical displacement at the epicenter occurs at t=22s according to spectral-element time reversal; this delay with respect to the reported earthquake origin time is comparable with the considered surface-wave period, and, in order of magnitude, with typical discrepancies between body-and surfacewave-based estimates of rupture times. The ray-theory simulation results in multiple maxima between 0 and 50 s. All this presumably reflects the complexity of surface-wave generation at the source, as well as errors introduced by the mentioned, non-physical propagation of the time-reversed wave field after focusing.

We repeat ray-theory time reversal in the 4-to-10 s and 20-to-30 s passbands, and show the results in Figs. 11 and12. Focusing of the time-reversed wave field is less sharp both in space and time (although, interestingly, in late snapshots of the time-reversal simulation (Fig. 11e,f), wavefronts are nicely centered on the earthquake epicenter). We ascribe the loss in source localization accuracy to the significant reduction in the width of the passbands, with respect to the previously discussed, 6-to-26 s simulation: we had anticipated in sec. 1 that focusing of the time-reversed wave field is enhanced by combining as many time-reversed surface-wave modes as possible. In our future work, we plan to more rigorously take advantage of this effect, multiplication surface-wave potential and their horizontal gradients by the radial Figure 10: Time-reversed signal at the epicenter of the Emilia earthquake as reconstructed by SPECFEM2D (blue curve) and ray theory (red) time reversal. Again, we define t=0 as the earthquake origin time; t should be intepreted as in Fig. 9, i.e. negative t corresponds to time after focusing in a time-reversal simulation.

eigenfunctions U (ω), V (ω), W (ω) according to eqs. ( 5) and ( 6), before integrating over the entire surface-wave frequency range.

Importantly, our analysis of time-reversed earthquake data shows that even at relatively short epicentral distances, where they are obscured by the body-wave coda, surface waves can still emerge in a time-reversal exercise. Focusing of the backward-propagated signal at the source can be thought of as a form of constructive interference. For time-reversed waves emitted at various station locations to interfere constructively, their backward propagation has to be modeled correctly. In our approach, time-reversed seismograms are filtered around one surface-wave frequency, and backward-propagated via the known Green's function (i.e.

phase-velocity map) for that frequency. In other words, only the propagation of time-reversed signal associated with surface waves at that frequency is modeled correctly, and it is only this signal that will contribute to "constructive interference" and to focusing of the timereversed wave field. Accordingly, circular wave fronts that can be associated with body-wave signal, and that do not focus at the epicenter (or elsewhere) are visible in Fig. 9a,b. We infer that surface-wave time reversal can indeed function as a source-imaging method also at relatively small epicentral distances, independently of how clearly surface waves can be identified visually on seismograms.

Conclusions

By taking advantage of the theory of surface-wave "potentials," we have reduced the problem of surface-wave propagation to two dimensions ("membrane waves"). We have shown that 3-D wave fields can then be reconstructed from monochromatic 2-D ones, once radial surface-wave eigenfunctions (sec. 2) are known; in this study, however, we only studied the propagation of surface-wave potentials in 2-D. We implemented a surface-wave time-reversal algorithm that can rely on either spectral-element or ray-theory models of wave propagation. In both cases, the theory is validated by application to real seismometer arrays in Central Europe. First, a synthetic test is implemented by computing approximately monochromatic membrane-wave seismograms at all receiver positions, from an arbitrary selected source location in Northern

Italy. In a second experiment, synthetic traces are replaced by approximate Green's functions, obtained by cross-correlating the real ambient signal recorded at one station of the array with difficulties in disentangling surface-wave signal from the body-wave coda. Importantly, our experiment described in sec. 6.3 suggests that time reversal and backward propagation using the surface-wave Green's function result in focusing of surface waves at the epicenter even at distances less than teleseismic, where surface waves carry less energy than body waves and their coda. These results encourage further applications of our method, in particular to the task of mapping, in space and time, rupture processes associated with large earthquakes.

  ) is zero if x A and x B are not within S [Baker and Copson, 1950, sec. 6.2].

Figure 1 :

 1 Figure 1: Rayleigh-wave phase-velocity maps of Kaestle et al. [2017], at periods of (a) 6 s, (b) 16 s and (c) 25 s.

Figure 2 :

 2 Figure 2: Snapshots of the ray-theory (left) and SPECFEM2D (right) synthetic-data timereversal simulations (sec. 6.1). Station locations are denoted by red triangles, the source location by a yellow circle. We define t=0 as the time when the source experiences the maximum displacement according to the Ricker wavelet in the forward simulations; backward propagation starts at the time corresponding to the last data sample employed in our exercise, and time increments in time-reversal simulations are considered to be negative. For each of the two time-reversal simulations, amplitudes are normalized to the maximum value obtained in the simulation, corresponding to source location at t=0. Snapshots a and b are taken at time t=65 s; c and d at t=0 s, e and f at t=-35 s. As explained in sec. 6, ray-theory and SPECFEM2D wave fields can be compared only qualitatively. Snapshots c and d show that the time-reversed wave field focuses at the "epicenter" location.

Figure 3 :

 3 Figure3: Synthetic test of sec. 6.1: Normalized time-reversed and backward-propagated displacement (dashed red curves) computed at the known location of the source, via (a) SPECFEM2D and (b) ray theory. In both cases, the known source time function is shown (blue curve) for comparison. In panel a, the difference between forcing term and reconstructed signal is explained by the fact that, in SPECFEM2D, displacement is initiated by prescribing a point force, rather than a displacement as in the ray-theory case.

Figure 4 :Figure 5 :

 45 Figure4: Snapshots of the ray-theory (left) and SPECFEM2D (right) time-reversal simulations of real noise cross-correlations described in sec. 6.2, in the 6-to-26 s period band. This is similar to Fig.2, but synthetic traces are replaced by cross-correlations of ambient data recorded at station LSD.GU (yellow circle) and all other stations whose locations are denoted by red triangles. Ambient-noise cross-correlations approximate the Green's function for each station pair, and, in this exercise, station LSD.GU can accordingly be thought of as a "virtual source." Snapshots a and b are taken at time t=82 s; c and d at t=0 s, e and f at t=-38. Snapshots c and d show that the time-reversed wave field indeed focuses at the location of station LSD.GU.

Figure 6 :

 6 Figure 6: Snapshots of time-reversal simulations of real noise cross-correlations, in the 4-to-10 s (left) and 20-to-30 s (right) period bands. As in Fig. 4, cross-correlated data were recorded at station LSD.GU (yellow circle) and all other stations whose locations are denoted by red triangles. Snapshots were selected at the same times as in Fig. 4. Snapshots c and d show that, also in these period bands, the time-reversed wave field focuses at the location of station LSD.GU.

Figure 7 :

 7 Figure 7: Same as Fig. 5, but traces obtained (via ray theory only) in the period bands 4-to-10 s (blue curve), 6-to-26 s (red) and 10-to-20 s (green) are shown.

Figure 8 :

 8 Figure 8: Normalized vertical-component recordings of the Mw=5.6 (MI=5.8) Emilia earthquake of May 29, 2012 [e.g.[START_REF] Molinari | Development and Testing of a 3D Seismic Velocity Model of the Po Plain Sedimentary Basin, Italy[END_REF], that we time-reverse and backwardpropagate as discussed in sec. 6.3. The vertical axis corresponds to epicentral distance, and each trace is plotted about its associated epicentral distance. All traces are Butterworthfiltered around 16 s as described in sec. 6.3.

Figure 9 :

 9 Figure 9: Snapshots of the ray-theory (left) and SPECFEM2D (right) time-reversal simulations of real earthquake data described in sec. 6.3. Again, the locations of stations utilized in the time-reversal simulation are denoted by red triangles, while the earthquake epicenter is marked by a yellow circle. We define t=0 as the earthquake origin time as reported by the Centro Nazionale Terremoti at INGV. Snapshots a and b are taken at time t=76 s; c and d at t=22 s, e and f at t=-50. Snapshots c and d show that the time-reversed wave field focuses at the epicenter of the earthquake.

Figure 11 :

 11 Figure 11: Snapshots of ray-theory time-reversal simulations of real earthquake data, in the 4-to-8 s (left) and 20-to-30 s (right) period bands. Snapshots were selected at the same times as in Fig. 9. All symbols are defined as in Fig. 9.

Figure 12 :

 12 Figure 12: The ray-tracing based trace of Fig. 10 (red curve, 6-to-26s period band) is compared to analogous traces obtained for the 4-to-8s (blue) and 20-to-30s (green) bands. Each trace is normalized to its maximum.

;

  [START_REF] Stehly | A study of the seismic noise from its long-range correlation properties[END_REF]. In many derivations of ambient-noise theory, the signals emitted by different noise sources are simply assumed to be random and uncorrelated, which results in the |h(ω)| 2 factor canceling out [e.g.[START_REF] Campillo | Seismic imaging and monitoring with ambient noise correlations[END_REF].

Acknowledgments

We are indebted to Julien de Rosny, Jean-Paul Montagner, Piero Poli, Kees Wapenaar for their insightful advice. Emanuel Kaestle provided his phase-velocity maps while Piero Poli shared his ambient-noise cross correlations. Several figures were created with the Generic Mapping Tools software [Wessel and Smith, 1991], and we made wide use of the ObsPy software [Krischer et al., 2015]. LB and MR are supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement 641943 (WAVES network). IM is supported by the Swiss National Science Foundation SINERGIA Project CRSII2-154434/1 (Swiss-AlpArray).

Appendix: dipole sources

The term "dipole source" refers here to the superposition of two impulsive point sources of opposite sign, located at two different points separated by a very small distance d. In this study, the concept of dipole emerges from the physical interpretation (sec. 4.2) of equation (34), relating the time-reversed, backward propagating wave field to the signals initally recorded by a receiver array. In the general context of wave physics, dipole sources are used, e.g., to formulate a modern, "corrected" version of Huygens' principle [START_REF] Baker | The Mathematical Theory of Huygens' Principle[END_REF][START_REF] Miller | Huygens's wave propagation principle corrected[END_REF].

The mathematical expression for a dipole source can be obtained by first writing the forcing term q defined in sec. 3 as the sum of two equal source distributions f (x, ω) shifted in space by the vector d (of magnitude d) and switched in sign one with respect to the other, i.e.

A first-order Taylor expansion around the point x S then gives

Substituting expression (47) into (46), we find

Finally, the sought expression is found by replacing f with a Dirac δ(x -x S ); since q accordingly becomes infinitely large at x S and zero elsewhere, the magnitude of d ceases to have meaning and d can be replaced by the corresponding unit vector d, so that

[e.g., [START_REF] Wapenaar | Elastic Wave Field Extrapolation-Redatuming of Single and Multi-Component Seismic Data[END_REF], sec. I.3.1].

Let us next find a simple expression for the response of a medium to dipole forcing. Recall that we have introduced the Green's function G 2D (x, x S , ω) in sec. 3 as the solution of eq. ( 23) with q(x, ω) =δ(x -x S ), i.e.

Applying the operator d • ∇ 1 to both sides of eq. ( 50) yields

We infer from eq. ( 51) that the solution of (23) with q(x, ω

Alternatively, [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF] (eqs. (E1)-(E3)) define the Green's function G 2D in the time domain as the solution of

Applying, again, d • ∇ 1 to both sides of eqs. ( 52)-( 54), we find that d • ∇ 1 G 2D (x, x S , t) is the field resulting from a dipole initial velocity at x S . This, or rather its Fourier transform ω), is what we call "dipole response" throughout this study.


Michael O'byrne
email: obyrnemj@tcd.ie

Franck Schoefs
email: franck.schoefs@univ-nantes.fr

Bidisha Ghosh
email: bghosh@tcd.ie

Vikram Pakrashi
email: v.pakrashi@ucc.ie

A Comparison between Pixel Intensity and Texture Analysis based Damage Detection Techniques

Keywords: Image processing, damage detection, comparative study, texture analysis

Imaging based damage detection techniques are increasingly being utilized alongside traditional visual inspection methods to provide owners/operators of infrastructure with an efficient source of quantitative information for ensuring their continued safe and economic operation. However, choosing a suitable damage detection technique that will provide a high level of quantitative information is often a timely and challenging prospect as there exists a wide array of algorithms currently available. These algorithms may be partitioned into one of two groups; pixel intensity based methods and texture analysis based methods. The algorithms in each group are naturally suited to different applications, depending largely on whether the damaged region under consideration is more separable from the background based on colour or on texture. This paper compares two algorithms, one from each category, which have previously been proposed in the domain of Infrastructure Maintenance Management (IMM). The algorithms are applied to a range of scenes featuring various forms of damage and their performance is investigated for best detection on the basis of performance coordinates in the Receiver Operating Characteristic (ROC) space. The general superiority of the pixel intensity based approach over the texture analysis based approach is demonstrated, in particular when considering High Dynamic Range (HDR) imagery.

Introduction

Regular inspections of structures are vital to ensure that they remain safe and fit for purpose. Currently, many structures are assessed using a regime of visual inspections. The levels of inspection cover a range of detail. It can be a cursory check, a principal inspection involving detailed examination of all surfaces, or special inspections where the employment of invasive, semi-invasive and non-invasive tests may become necessary. Most inspections still depend on visual observations and the quality of data collected largely depends on the ability of the inspectors to observe and objectively record details of damage. The approach is prone to considerations such as operator boredom, lapses in concentration, subjectivity, and fatigue, which contribute to the variability and reduced accuracy of visual inspections [START_REF] Agin | Computer Vision Systems for Industrial Inspection and Assembly[END_REF][START_REF] Komorowski | Role of Enhanced Visual Inspections in the New Strategy for Corrosion Management[END_REF]. With this in mind, it is often desirable to incorporate a second, more objective and quantitative, source of information about the health condition of a structure. Non-Destructive Techniques (NDT) often provide the only method of obtaining such information.

The information obtained from the NDT may be fed into an Infrastructure Management System (IMS), which can help the decision makers to make more effective and informed judgments when allocating resources towards the correction of deficiencies and when choosing an appropriate future course of action. This aspect has attracted a growing interest in recent years as the importance of life cycle optimisation and the related financial benefits continue to be recognised [START_REF] Jr | Cost Optimization of Prestressed Concrete Bridges[END_REF][START_REF] Schoefs | Comparison of Additional Costs for Several Replacement Strategies of Randomly Ageing Reinforced Concrete Pipes[END_REF]. For a well calibrated IMS, it is important that the input information is accurate and comprehensive. This requires selecting the most suitable NDT technique, which for a given application is not always readily apparent as a measure of the onsite performance of an NDT technique remains a pertinent question in the majority of cases [START_REF] O'byrne | Regionally Enhanced Multi-Phase Segmentation Technique[END_REF]. The choice of NDT will largely depend on the damage to be detected and will require an in-depth knowledge of the advantages and limitations associated with each option. Even NDTs which fall under the same category can produce markedly varying results. Unfortunately, many inspections proceed without a proper evaluation of all available options and as a result an inferior level of information is obtained.

In recent times, image processing based damage detection methods have been increasingly considered as a viable NDT option. This is due to a number of reasons. Image processing methods use inexpensive and readily available equipment (i.e. a standard digital camera), and they do not require the inspector to undertake extensive training. Furthermore, advances in camera technology mean that rich detailed imagery of damaged components can be acquired. Additionally, visual inspections almost always capture photographs to include in the inspection report to corroborate the inspector's comments; however, these photographs are rarely exploited to their fullest potential in either a qualitative or a quantitative fashion. The primary limitations image processing methods are the lack of penetration below the surface of the material and the requirement of good visibility and lighting conditions.

Applying damage detection algorithms to the photographs can locate and quantify visible mechanical damage on the surface of infrastructural elements with minimal human supervision. Physical properties of the identified damage, such as the size and shape characteristics, may be easily extracted with knowledge of a real world scale. The quantitative nature of the data obtained from image analysis is important and naturally lends itself to numerous applications. It is helpful for developing new damage models, or strengthening existing ones, which are used to forecast the rate of propagation of damage as the structure continues to operate.

Most image processing based damage detection algorithms consist of segmentation followed by subsequent classification of the segmented regions. Ideally, the segmentation methodology should identify and accurately define all regions of interest in an image whilst minimizing the inclusion of extraneous regions. In reality, perfect segmentation is difficult to achieve given the inherent chromatic and luminous complexities encountered in natural scenes. Segmentation algorithms use either pixel intensity (colour) information or texture information to isolate similar regions in an image. The effectiveness of colour based segmentation algorithms and texture based segmentation algorithms will vary according to the surface and damage type under consideration as certain damages are more separable from the undamaged surface based on either their colour or texture attributes. This paper evaluates and compares the performance of a pixel intensity based method [START_REF] O'byrne | Regionally Enhanced Multi-Phase Segmentation Technique[END_REF] and a texture analysis based method [START_REF] O'byrne | Texture Analysis Based Damage Detection of Ageing Infrastructural Elements[END_REF] which have previously been proposed in the domain of IMM. Both methods are applied to four different scenes featuring various damage forms, lighting conditions, viewing angles, resolutions etc. Additionally, High Dynamic Range (HDR) imagery is adopted as a protocol, as proposed by [START_REF] Ghosh | High Dynamic Range Image Processing for Non-Destructive-Testing[END_REF], in an attempt to optimise the detection accuracy of each method.

The following section provides an overview of each of the detection techniques as well as providing a brief background to HDR. This section also details the performance evaluation process which is based on performance points plotted in the Receiver Operating Characteristic (ROC) space. Section 3 introduces the imagery used while Section 4 presents the results obtained from each method. Section 5 concludes the paper.

Methodology

A comparison study between two previously proposed image processing based damage detection methods has been conducted in this paper; a pixel intensity based segmentation method and a texture analysis based method. The comparison involves two steps. Firstly, each detection method is applied to the four Standard Dynamic Range (SDR) images and the corresponding HDR images. Secondly, the detection accuracy obtained from each method is calculated for each image and a comparison is then drawn. The detection methods are briefly described below.

Pixel Intensity based Segmentation Method

Pixel intensity based segmentation algorithms may be grouped into four major categories: thresholding, edge detection using gradient information, region growing, and hybrid methods [START_REF] Abdel-Qader | Segmentation of Thermal Images for Non-Destructive Evaluation of Bridge Decks[END_REF]. The method used in this paper, known as REMPS (Regionally Enhanced Multi-Phase Segmentation), is a hybrid method which integrates three feature detection stages. The first stage involves the application of the Sobel edge detector on a pre-processed image in order to form closed geometries corresponding to objects in a scene. Statistical properties are then calculated for each closed geometry, which are used as the basis for the clustering based filtering phase. This phase retains closed geometries deemed to represent damaged regions. Finally, Support Vector Machines SVMs are used to identify pixels having intensity values characteristic of damaged zones. These pixels are then applied locally to the filtered closed geometries in order to improve the definition of the detected damaged regions. REMPS attempts to utilise the advantages of these three independent techniques most effectively and extract their mutual benefits. For instance, the robustness and generality of the Sobel edge detector serves as a natural precursor to the closed geometry clustering stage. This clustering stage performs well at classifying the presence of damage, however, it is only after the pixel supplementation stage that the shape and size characteristics of the retained closed geometries are sufficiently realised. A flowchart illustrating the order of the feature detection methods is presented in Figure 1(a).

Texture Analysis based Segmentation Method

Texture is an innate property of surfaces which, for human observers, texture may be qualified by terms such as fine, coarse, smooth, rippled, molled, irregular, or lineated [START_REF] R M Haralick | Textural Features for Image Classification[END_REF]. There are numerous ways to quantify texture; wavelet analysis, Laws' texture energy, First Order Statistics (FOS) and Grey Level Co-occurrence Matrix (GLCM). The texture analysis based segmentation method used in this paper combines both FOS and GLCM statistics. The texture dependent statistics, or texture measures, are calculated at every pixel in an image using a sliding window approach. Four statistics were derived from a GLCM; angular second moment, homogeneity, contrast and correlation, while a further six texture measures were calculated directly from the pixel intensity values from the original image. These were Shannon entropy, mean, variance, range, skewness and kurtosis. The 10 texture measures calculated at each pixel location form a feature vector. Non-linear SVM models are used to classify pixels as either damaged or undamaged based on the feature vector corresponding to that pixel. Two SVM classification models were considered; a Custom-Weighted Iterative (CWI) model and a 4-Dimensional Feature Space (4DFS) model in which the feature vectors were mapped to a four dimensional feature space. The 4DFS model was generally shown to produce better segmentation results so this model was used for the purposes of this comparison. The method is illustrated in the following flowchart Figure 1(b).

High Dynamic Range (HDR)

HDR is a set of techniques that allow a greater dynamic range of luminance values between the brightest and darkest regions of an image than standard digital images. SDR images can typically only accommodate a very limited range bracket of the full tonal spectrum in a real world scene. Therefore, a dynamic range bracket would have to be chosen in the knowledge that all luminance values outside the range would not be represented correctly. The broad principle behind HDR imagery is that multiple SDR images of the same scene, each taken at a different exposure, and thus capturing a different range bracket of the tonal spectrum, may be merged to form one HDR image that has a wider dynamic range [START_REF] Reinhard | Color Imaging: Fundamentals and Applications[END_REF]. Combining SDR images can be done using various merging algorithms [START_REF] Debevec | Recovering High Dynamic Range Radiance Maps From Photographs[END_REF].

The usefulness of HDR imagery as an imaging protocol may be observed in Figure 2 which depicts the pixel intensity values plotted along a profile line for a normally exposed SDR image and the corresponding HDR image (the SDR and HDR images in question are from Figure 4 (1) and ( 5) respectively). It may be observed that HDR image exhibits an increased tonal range compared to the SDR image, thus offering an altogether superior information content.

Performance Evaluation

The performance of each method is evaluated through the use of performance points in the Receiver Operating Characteristic (ROC) space. The ROC space allows for a convenient means for characterising and comparing the performance of NDTs in various conditions [START_REF] Rouhan | Probabilistic Modeling of Inspection Results for Offshore Structures[END_REF] and has been recently expanded to image detection [START_REF] Pakrashi | ROC Dependent Event Isolation Method for Image Processing Based Assessment of Corroded Harbour Structures[END_REF]. For any NDT, the Detection Rate (DR) along with the accompanying Misclassification Rate (MCR), or alternatively known as Probability of Detection (PoD) and Probability of False Alarm (PFA) in the field of probability space and decision theory, are determined by comparing the corroded regions detected with a visually segmented image, which acts as the control. The DR and MCR are represented as a percentage between 0% and 100%. Each (MCR,DR) pair formed a coordinate in the ROC space.

There are a few measures for comparing segmentation performance [START_REF] Hui | Image Segmentation Evaluation: A Survey of Unsupervised Methods[END_REF]. In this paper, a measure of the performance was attained through the use of the α-δ method [START_REF] Baroth | Construction Reliability[END_REF][START_REF] Schoefs | The αδ Method for Modelling Expert Judgment and Combination of NDT Tools in RBI Context: Application to Marine Structures[END_REF]. This method relies on calculating the angle, α, and the Euclidean distance, δ, between the best performance point, defined as an ideal NDT with 100% detection and 0% misclassification rates and represented in the ROC space with coordinates (0,1) and the considered point to give a measure of the performance of the considered point. As this paper is not devoted to risk analysis where the shape to the ROC acts as a key factor, only the delta, δ, parameter is required. A low value for δ is indicative of a strong performing technique.

Data Analysis

The pixel intensity and texture based segmentation methods were applied to four images featuring various forms of damage on the surface of infrastructural elements. In order to provide a more meaningful comparison, the images were chosen to reflect a broad range of surfaces, damage forms, viewing angles, lighting conditions and image resolutions as shown in Figure 3. The sample images in the figure depict, (1) pitting corrosion on metal sheet piling in marine conditions, (2) marine growth on the surface of underwater steel pile wharf, (3) corroded metal sheeting , and (4) an exposed concrete bridge deck through wear of pavement surfacing. The HDR images are also displayed (5 -8). It is readily apparent that the HDR images offer more detail than their SDR counterparts. The next section presents the results obtained from each detection method performed on these images.

Results

Both detection methods required two SVM parameters as input parameters which ultimately had an influence on the detection accuracy. In order to achieve a fair comparison, a parameter search was performed according to the method outlined in each paper to find sufficiently optimized parameters. Following this optimization procedure, the methods were applied to each image in Figure 3 resulting in the detected regions shown in Figure 4. Their performance levels are summarized in Table 1 and the corresponding performance points are plotted in the ROC space in Figure 5.

Discussion

It may be noted from the detected regions in Figure 4 that the pixel intensity based method was quite successful for the majority of cases with the exception of the marine growth image. The poor detection results for this image may be explained by the fact that the damaged regions throughout the image were not characterized by one single colour. Instead they took on numerous contrasting shades which often overlapped with the non-damaged background. Generally however, this method proved effective at locating the presence of damage as well as accurately defining the shape and size of damaged regions.

The texture based method was not as effective in terms of detection accuracy as the pixel intensity based method. The only case where it was on a par with the pixel intensity method was for the image of the marine growth. It also suffered from detecting many small spurious regions unlike the pixel intensity based method which a 'cleaner' and more homogenous detection.

Analysis of the δ values in Table 1 reveals that HDR imagery does indeed improve the accuracy of the pixel intensity based method. Texture based segmentation on the other hand did not appear to benefit from the use of HDR as on two occasions, the best detection results emerged when using SDR images while on the other two occasions, HDR images provided better results. This suggests that adopting a HDR protocol is not especially relevant when using texture analysis as a detection method.

Furthermore, the aspect of computational time may be an important factor for some, especially in cases where a large batch of images are required to be processed. It was found that the colour based segmentation method had a superior computational efficiency over the texture based algorithm, which is due to the fact that texture must be calculated by considering a collection of neighbouring pixels around each pixel while colour based segmentation techniques typically need only consider each pixel intensity value independently.

Conclusion

This paper presents a comparison between two image based damage detection methods. This study was necessitated by the vast array of damage detection methods currently available, which makes choosing an appropriate method a timely and daunting prospect. This paper gauges the effectiveness of two contrasting approaches, thereby providing an insight into the expected performance levels that can be attained by anybody who is considering incorporating an image processing technique in their inspection routine. The two methods were based on detection via pixel intensity (colour) information and detection via texture information. While both of these methods are not necessarily reflective of all pixel intensity and texture analysis based methods currently available, the comparison does serve to underline how each approach may typically respond, especially when the damage type is known beforehand.

The main forms of surface damage encountered on ageing infrastructural elements (corrosion, leaching, etc.) are often characterised to a greater extent by the change in colour from the undamaged surface than by a change in texture. With this in mind, it is not surprising that the pixel intensity based method demonstrated a higher degree of success for the majority of the cases explored in this paper. Texture based segmentation fared worse at isolating damage, so it may be can be classified as suitable for specific applications such as marine growth where the damaged regions do not share a distinct colour but are characterised more so by their rougher texture than the surroundings.

The presented results indicate that improvements can be made to the detection accuracy of the pixel intensity based method by adopting a HDR protocol. The texture analysis based method however did not experience any noticeable gain by using HDR.

Figure 1 .

 1 Figure 1. (a) Flowchart for the pixel intensity based segmentation method, and (b) flowchart for the texture analysis based segmentation method

Figure 2 .

 2 Figure 2. The increased tonal range during damage detection for an HDR image as compared to an SDR image with normal

Figure 3 .

 3 Figure 3. SDR images of various forms of damage (1 -4), and the corresponding HDR images (5 -8)

Figure 3 .

 3 Figure 3. Detected regions for the pixel intensity and texture analysis based segmentation methods

Figure 4 .

 4 Figure 4. Comparison of detection methods through the use of performance points in the ROC space

Table 1 . Detection accuracy of the pixel intensity and texture analysis methods Sample Image Pixel Intensity based Segmentation

 1

	Texture based
	Segmentation

Acknowledgements

The authors wish to thank the Irish Research Council for Science, Engineering and Technology (IRCSET) for providing grant to support this research and CAPACITES/IXEAD society for the practical assistance.


Matthew J Phillips

Frédéric Delsuc

David Penny

Genome-scale phylogeny and the detection of systematic

biases.

Supplementary Information

 [START_REF] Felsenstein | Phylogenies from molecular sequences: inference and reliability[END_REF]). The magnitude tests include (c.)

treeness: the proportion of uncorrected minimum-evolution tree distance attributed to internal (as opposed to external) branches and (d.) relative composition variability (RCV; [START_REF] Phillips | The root of the mammalian tree inferred from whole mitochondrial genomes[END_REF]. The chi-square tests and RCV values were determined with constant sites excluded, as these sites buffer against inferring composition bias among the sites that contribute directly to apparent phylogenetic signal.

16/01/04, page 2 Figure S1. Minimum evolution tree from pairwise AT-GC base-frequency differences.

Following [START_REF] Kuhner | A simulation of phylogeny algorithms under equal and unequal evolutionary rates[END_REF]