- Email: [email protected]
Source point discovery through high frequency asymptotic time reversal
Journal of Computational Physics 231 (2012) 4643–4661
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Source point discovery through high frequency asymptotic time reversal Jean-David Benamou a,⇑, Francis Collino b, Simon Marmorat a a b
INRIA, B.P. 105, 78153 Le Chesnay Cédex, France CERFACS, 42 Avenue Gaspard Coriolis, 31100 Toulouse, France
a r t i c l e
i n f o
Article history: Received 18 January 2011 Received in revised form 2 March 2012 Accepted 3 March 2012 Available online 4 April 2012 Keywords: Linear waves Geometrical optics Time reversal
a b s t r a c t Given local frequency domain wave data, the Numerical Micro-Local Analysis (NMLA) method (Benamou et al., 2004) [5] and its recent improved version (Benamou et al., 2011) [4] gives a pointwise numerical approximation of the number of rays, their slowness vectors and corresponding wavefront curvatures. With time domain wave data and assuming the source wavelet is given, the method also estimates the travel-time. The paper provides a non technical presentation of the improved NMLA algorithm and presents a numerical application which can be interpreted as a high frequency asymptotic version of the classical time reversal method (Borcea et al., 2003) [7]. A detailed technical presentation of the algorithm is available in Benamou et al. (2011) [4] and more numerical experiments can be found in Collino and Marmorat (2011) [15]. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction Source point or scatterer location identification is a classical inverse problem in computational wave propagation. It is of great practical interest in, for example, non destructive testing of alloys defect, kidney stones location discovery [28], underwater communications [23]. It is a basic step in inverse synthetic aperture radar [8]. See also [24] for the same problem of source discovery. The popular time reversal (TR) method (see [7] for a mathematical study and references) solves this inverse problem. It is rapidly summarized in the following steps: (i) record a wave signal on some array of ‘‘transducers’’ for ‘‘some’’ time, (ii) time reverse this signal and re-emit it in the medium. Then, physics observes and maths proves that the time reversed signal will exhibit refocusing properties towards the scatterers or the source points depending on the set up of the experiment. We propose a time reversal method based on the High Frequency (HF) asymptotic model of wave propagation, also known as Geometric Optics (GO). The HF model only makes sense with smooth deterministic media when the characteristic wavelength of the signal is smaller than the variation scale of the acoustic speed. We will therefore NOT be concerned with the widely studied ‘‘super-resolution’’ effect of the focusing observed in random media [7]. The goal of the paper is twofold: Provide a non technical description of the Numerical Micro-Local Analysis algorithm (NMLA) [5] and its recent improved second order version [4]. We apply it to analyze the underlying GO ray directions of the signal recorded in step (i). The ⇑ Corresponding author. E-mail address: [email protected] (J.-D. Benamou). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.03.012
4644
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
NMLA algorithm works in the frequency domain. We therefore use classic pre and post processing techniques to extend it to time domain and also to treat a signal recorded on linear arrays of transducers. We show that the use of time domain data simplify multiple rays discrimination and that frequency redundant information provides a travel-time estimate. Illustrate numerically that the output ray direction and travel-time are accurate enough, at least with synthetic noisy data, to shoot the ray back to its unknown source location. The notion of ‘‘shooting back’’ implies that the travel-time, the actual propagation time, is reversed in the same way as in step (ii) above. Finally, we compare the ray directions to those found with more classical signal processing methods like Plane Wave Destructor [21] or the Local Slant Stack methods [8,19]. Potential applications of the NMLA method are of course not limited to source point discovery. Here is a non exhaustive list: Micro-local finite element methods try to discretize the geometric optic quantities [1,17]. If known, an analytic solution for the phase / can be use to remove oscillations and the size of the discretization [12,25,22]. GO is simple and cheap when the underlying wave field is not too complicated: no diffraction, reflections, caustics, etc. In these more complicated regime, GO has to manage several complex interacting ansatz. As simple propagating waves can be treated locally by GO, coupling GO and full wave fields across physical interfaces or in a multi-scale framework make sense [26,2]. A reliable process to invert the GO ansatz is needed. In Geophysics ‘‘Cleaning’’ data by filtering out low frequencies either to use faster HF models or simply because only high frequencies are really relevant (migration) or nicer for inversion (see [30,29] for example). The Direction of Arrival problem (DOA) is a classic signal processing problem [27] where one tries to find the direction of the source of a signal impinging on a, generally linear, array of antennas. Source point reduced models for far field representations [9]. In order to make the paper self contained, it is organized as follows: Section 2 recalls the HF model derivation and formulates the HF asymptotic time reversal source discovery inverse problem. Section 3 reviews the refined NMLA algorithm [4] and its extension to time domain data recorded on a linear array of ‘‘transducers’’. Section 4 presents numerical source discovery from synthetic data. 2. Forward and inverse problem 2.1. High frequency acoustic wave equation We start from the classical acoustic wave equation in the time domain. The wave field u(x, t) is generated by a finite time wavelet emitting from a source point:
8 < c21ðxÞ @ 2tt uðx; tÞ Duðx; tÞ ¼ dðx xs ÞW xs ðtÞ; 8ðx; tÞ 2 R2 ½0; TÞ; :
uðx; t ¼ 0Þ ¼ @ t uðx; t ¼ 0Þ ¼ 0;
ð1Þ
8x 2 R2 ;
where
[0, T) is the simulation time (one may consider T = +1). x ´ c(x) is the acoustic speed law function on R2 . As we deal with HF–GO models, we will consider smooth speeds. xs is the source point location. t # W xs ðtÞ is the source time wavelet. The radian frequency x represents the dual variable of time in the Fourier transform denoted^: (we assume u 0 for t < 0):
^ ðx; xÞ ¼ u
Z
1
uðx; tÞeixt dt:
ð2Þ
0
The source wavelet, denoted W xs ðtÞ, determines the frequency content of the solution. It depends explicitely on the ‘‘characteristic pulsation (or frequency)’’ xs. We define it as
W xs ðtÞ ¼ Wðxs tÞ; where W(t) is a shape function with support in [0, 1]. Assuming further that W is absolutely integrable and has zero mean, the following properties are satisfied:
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
1 ; SupportðW xs Þ ¼ 0; xs c x : c x ð xÞ ¼ 1 W W s For
ð3Þ ð4Þ
xs
xs
> 0, there exist two positive constants cþ Z
xRI
c x ðxÞj2 dx 6 jW s
Z
4645
c x ðxÞj2 dx; jW s
and c such that
I ¼ ½c xs ; cþ xs :
R
ð5Þ
So, when xs tends to infinity, (3) shows that the source time support shrinks and (4 and 5) that the frequency content grows and but is shifted towards infinity. 2.2. Frequency domain and ansatz The HF model uses a formal asymptotic representation of the solution in which the unknowns are independent of the radian frequency x. It is easier to derive and interpret it in the frequency domain: Assume T = +1 and Fourier transform equation (1):
x2 ^ c x ðxÞ: ^ ðx; xÞ ¼ dðx xs Þ W uðx; xÞ Du s 2
c ðxÞ
ð6Þ
We also assume that c(x) = v is constant outside a compact in X and add the Sommerfeld radiation boundary condition at infinity
^ ðx; xÞ lim @ jxj u
jxj!1
ix
v
^ ðx; xÞ ¼ 0: u
ð7Þ
Problems (6) and (7) are well posed [16]. In the homogeneous case, c(x) = v constant everywhere, the analytic solution of (6) and (7) is given by
^ ðx; xÞ ¼ u
i ð1Þ jx xs j c W xs ðxÞ; H0 x 4 v
ð8Þ
ð1Þ
where H0 ðuÞ is the Hankel function of first type and order 0. Plugging in (8) the following asymptotic estimate (cf. [31])
sffiffiffiffiffiffiffi 2 iðqp4Þ ; qÞ e
ð1Þ H0 ð
q!1
pq
we obtain an asymptotic approximation of the solution in the form:
c ð xÞ W s ^ ðx; xÞ ’ Aðx; xs Þeix/ðx;xs Þ pxffiffiffiffiffiffiffiffiffi ffi ; u ix
ð9Þ
where
Aðx; yÞ ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
v
8pjx yj
/ðx; yÞ ¼
;
jx yj
v
:
ð10Þ
sj The large parameter is x jxx v . It represents the number of wavelengths between the source and the observation point. In the time domain (9) corresponds to the ‘‘ray’’ solution
1 Aðx; xs Þ uray ðx; tÞ ’ pffiffiffiffiffiffi W 2 ðxs ð/ðx; xs Þ tÞÞ;
xs
1
where W 2 ðtÞ is
ð11Þ
RW b ðxÞ pffiffiffiffiffiffi eixt dx, the half derivative of the wavelet. ix
The ‘‘ansatz’’ idea, apparently going back to a 1911 paper of Debye on ‘‘the article of Sommerfeld’’ (see [20]), consists in generalizing the a priori form (9) to the general heterogeneous case and search for a formal asymptotic solution of the Helmholtz equation in the form
^ ray ðx; xÞ ’ eixuðxÞ u
1 A ðxÞ P j
xÞ j
j¼0 ði
:
ð12Þ
pffiffiffiffiffiffiffiffiffiffi We assume the source wavelet and the half primitive (symbol ix) have been removed, see Section 3.1. In general, only the first A0 term of the amplitude series is used. Oscillations are formally enforced in the phase and travel at travel-time u. The radian frequency x is the large parameter.
4646
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
2.3. Geometric optics–Ray Tracing We give a quick overview to make the paper self contained, for more details and references see [20]. Assuming a solution taking the ansatz form (12) locally away from the source. One can plug it into
^ ðx; xÞ þ Du
x2 ^ uðx; xÞ ¼ 0
c2 ðxÞ
ð13Þ
and expand in inverse power of x. Canceling the coefficients of the first and second order term of the series gives respectively the Eikonal equation and amplitude equations
(
1 jruðxÞj ¼ cðxÞ ;
2ruðxÞ rA0 ðxÞ þ A0 ðxÞDuðxÞ ¼ 0:
ð14Þ
These equations are local but global solutions are constructed using a Lagrangian method called Ray Tracing. The travel-time u is computed along the integral curves of the gradient called rays and denoted s ? y(s, y0). Each ray may be parameterized by its initial position y0 and, assuming the travel-time and its gradient are known at this location, one solves the following system of ODEs:
8d > < ds yðs; y0 Þ ¼ ruðyðs; y0 ÞÞ :¼ pðs; y0 Þ; yð0; y0 Þ ¼ y0 ; d pðs; y0 Þ ¼ 12 rn2 ðyðs; y0 ÞÞ; pð0; y0 Þ ¼ ruðy0 Þ given; ds > :d 2 u ðyðs; y ÞÞ ¼ n ðyðs; y ÞÞ; uðy0 Þ given; 0 0 ds
ð15Þ
where y(s, y0) is the ray shot from y0 and s is a parameterization of the curve. p = r/(y) is called the ‘‘slowness’’ vector and gives the ray direction. n ¼ 1c is called the slowness of the medium. When the velocity is homogeneous (c(x) = v) around a source point, the analytic asymptotic form (10) can be used to define the initial conditions: For s small enough the level curve W C ¼ fy0 ; s:t:/ðy0 Þ ¼ jy0 vxs j ¼ sg is a circle around the source sÞ . Letting s go to 0 we consider that the HF asymptotic form of a source point corresponds point and pð0; y0 Þ ¼ ruðy0 Þ ¼ vðyky0 x xs j 0
to rays leaving the source point xs isotropically with the same amplitude. We then switch the parameterization of the rays to the initial take off angle y0 ? h0. The global HF solution is formed of all rays; this is the solution, for all h0 2 [0, 2p[, of
8d > < ds yðs; h0 Þ ¼ ruðyðs; h0 ÞÞ :¼ pðs; h0 Þ; yð0; h0 Þ ¼ xs ; ~ d pðs; h0 Þ ¼ 12 rn2 ðyðs; h0 ÞÞ; pð0; h0 Þ ¼ sv0 ; ds > :d uðyðs; h0 ÞÞ ¼ n2 ðyðs; h0 ÞÞ; uðxs Þ ¼ 0; ds
ð16Þ
where ~s0 ¼ ðcos h0 ; sin h0 Þ. Ray Tracing (16) can be extended to compute solutions of the amplitude equation (see [3] for instance). In smooth heterogeneous media, the Lagrangian travel-time function remains smooth but it is quite common for rays to cross. Then, solutions to the Eulerian equation (14) become multi-valued. In the simplest case, the HF ansatz remains a superposition/sum of elementary (12) ansatz. See [18,3] for more details and references on multi-valued solutions. The NMLA algorithm assume that its data, the Helmholtz solution, can be locally represented as a general HF ansatz of the form
^ ðx; xÞ ’ u
N P
A0;n ðxÞeixun ðxÞ ;
ð17Þ
n¼1
where N is the number of rays going through x and n the subscript labeling the different branches of the solution. We kept only the zeroth order term of the amplitude series of (12). 2.4. Source point discovery with HF time reversal: formulation and method of resolution In order to simplify the presentation of the source point discovery and resolution method, we will restrict to the case N = 1. We give NMLA computations for N = 2 in Section 4.2.2 and discuss ray separation using time and direction windows in Appendix C. The source point discovery problem can be formally formulated as follows: Let u be the solution of (1). The source point location xs is the unknown we are looking for. The wavelet W is assumed to be given (see Appendix B for a discussion on the method when it is also unknown). We fix an observation point x0 ¼ xo1 ; xo2 such that c(x) = v is constant for all x ¼ ðx1 ; x2 Þ; x2 > xo2 . We assume further that there is a time t0 such that u can be represented, in a neighborhood of (t0, x0) as an asymptotic ray solution in the sense of (11):
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661 1 AðxÞ uðx; tÞ ’ pffiffiffiffiffiffi W 2 ðxs ðuðxÞ tÞÞ þ duðx; tÞ:
xs
4647
ð18Þ
In (18), we have defined two new auxiliary unknowns A and u which satisfy the GO equations (14) and the Ray Tracing equations (16). We assume the remainder du and all its derivatives to be negligible with respect to u as xs ? +1. Our method to locate xs first uses the local data (18) to compute the ray unknowns u(x0) and pdata ¼ rudata ðx0 Þ. Then, we 0 use the global reversible rays equations to back propagate the ray to the source point from these initial conditions. This is accomplished by changing ds in ds (16) and use the initial conditions obtained in (i) for this particular ray:
8 d > < ds yðsÞ ¼ ruðyðsÞÞ :¼ pðsÞ; yð0Þ ¼ x0 ; d pðsÞ ¼ 12 rn2 ðyðsÞÞ; pð0Þ ¼ pdata 0 ; ds > : d 2 uðyðsÞÞ ¼ n ðyðs; h0 ÞÞ; uðxs Þ ¼ udata ðx0 Þ: ds
ð19Þ
The travel-time u(y(s)) will then decrease and we stop integrating (19) when it gets to 0. We have ‘‘time reverted’’ the ray and y(s) has reached the source xs. We emphasize that the medium slowness n ¼ 1c must be known. The ray unknowns u(x0) and pdata ¼ rudata ðx0 Þ are computed using the NMLA algorithm. This is our main contribution 0 ^ ðx; xÞ on a local observation circle around the choand it is described in Section 3. The method uses frequency domain data u sen observation point x0. Common numerical discretizations of wave fields directly do not provide data on a circular array around an ‘‘observation’’ point. Regarding real data, circular arrays do not seem to be a common acquisition device either. Interpolation is a solution but it needs data on a volume covering the acquisition circle which typically is several wavelengths large. Considering the most common acquisition format is a linear array, like for instance seismograms recorded on an ideal flat surface in geophysics, we propose here to use one way wave extrapolation to reconstruct the circular array from a seismogram from this type of data. It implies that rays of interest must be transverse to the line of acquisition but this is not very restrictive. We summarize below the data pre-processing steps to generate NMLA compliant data: The acoustic signal (18) is recorded on a line or a linear array which contains a chosen observation point x0 ¼ xo1 ; xo2 . In o 2D, it gives time data u x1 ; x2 ; t for t in some time interval [0, T] and x1 on the line (we denote x = (x1, x2)). As we deal with time data it is possible to preselect a given time window centered around a chosen time of interest t0 by ~ x1 ; xo2 ; t ¼ u x1 ; xo2 ; t gðt t0 Þ were g is a smooth cut-off window with support larger than the source waveconsidering u let’s W xs . ^ ~ x1 ; xo2 ; x ¼ F t!x u ~ x1 ; xo2 ; t to go to frequency domain and choose a frequency x Next we Fourier transform the data: u such that the frequency content of the source is significant, typically close to xs. ^ ~ x1 ; xo2 ; x to reconstruct the field and its gradient in a neighborhood of x0, (see We use a paraxial extrapolation in x2 of u Appendix A). Remark 1. All the procedure only makes sense if the HF representation (18) is relevant. If not it simply means that the HF model is not suitable to approximate the solution, i.e. there are no rays in the solution. Remark 2. Classical antenna localization algorithms (GPS, cell phone) use estimated travel-times from different sources with known locations (satellites) to draw spheres which intersections ‘‘triangulate’’ the antenna position. The wave signal directions between the antenna and the sources are not known. Our approach is based on finding both the direction and the travel-time from a single source and then back propagate ‘‘exactly’’. Repeating the resolution method at different observations points and using different rays is of course also possible. 3. Numerical micro-local analysis NMLA does not directly invert (17) but a system of linear equations obtained through the successive steps described in the following subsections. 3.1. Source deconvolution When the wavelet W xs is known, one can perform the following source deconvolution on the data
^ ðx; xÞ u
pffiffiffiffiffiffiffiffiffiffi ix ^ ðx; xÞ: u c x ðxÞ W
ð20Þ
s
When the source wavelet is unknown, the method still works but the travel-time cannot be estimated, see Appendix B. 3.2. Plane wave approximation In the GO regime, we expect the phase to be real and with variations on a scale larger than x. Within a few wavelength of x0 the constant amplitude approximation A0(x) ’ A0(x0) will remain of order x1 .
4648
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Using a Taylor expansion on the travel-time and in the vicinity of x0 we have
1 2
uðxÞ ’ uðx0 Þ þ ðx x0 Þ ruðx0 Þ þ ðx x0 ÞT Huðx0 Þðx x0 Þ þ
ð21Þ
The usual first order ‘‘plane wave’’ approximation is, for N = 1 in (17),
^ ðx; xÞ ’ Bðx0 Þeixðxx0 Þruðx0 Þ ; u
ð22Þ
where
Bðx0 Þ ¼ A0 ðx0 Þeixuðx0 Þ : ^ can be identified as the N-rays ansatz In the more general case where N rays cross at x0, we assume u
^ ðx; xÞ ’ u
N P
Bn ðx0 Þeixðxx0 Þrun ðx0 Þ :
ð23Þ
n¼1
3.3. Relaxation Our goal is to identify (Bn, run) and N the number of rays which appear in Eq. (23). In order to simplify this non-linear inverse problem we instead consider the linear problem of finding the amplitude distribution.{Bm, m = 1, . . . , M} over a dis~m ¼ m ¼ 1; . . . ; Mg with cretization of the set of 2-D unit vectors, fd
~m ¼ ðcos hm ; sin hm Þ hm ¼ m 2p ; d M
m ¼ 1; . . . ; M
ð24Þ
The ~: notation will always indicate unit vectors. As the norm of the gradient of the travel-times is given by the Eikonal equation,
jrun ðx0 Þj ¼
1 cðx0 Þ
for all n;
the relaxed version of (23) is
^ ðx; xÞ ’ u
M P
~m i x ðxx0 Þd
Bm ðx0 Þe cðx0 Þ
ð25Þ
:
m¼1
The size of the discretization M is discussed below. 3.4. The observable and the linear system ^ ðx; xÞ on a scaled circle around x0. The discretization (24) was used to construct a In our first paper [5] we evaluated u finite dimensional linear system (see the geometry os the support on Fig. 1). It was subsequently inverted to obtain the amplitude distribution Bm. The data for the linear system was given by
~ ; xÞ r ¼ acðx0 Þ ^ ðx0 þ r d U a;m ¼ u m
x
j~sj ¼ 1:
ð26Þ
This choice sometimes led to unstable inversions and a L2 regularization was used in [5]; a different L1 regularization was proposed in [24]. The new version proposed in [4] still uses the circle geometry but the sampled data is in impedance form: n o cðx0 Þ @ u ^ ^ . The resulting NMLA filter is stable independently of its parameters. The new observable therefore is: þu ix @r
Fig. 1. Geometry of the observable support, a sample point in the direction ~s. r ¼ a cðxx0 Þ.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Um ¼
^ cðx0 Þ @ u ~ m ; xÞ þ u ~m ; xÞ; ^ ðx0 þ r d ðx0 þ rd ix @r
r¼
acðx0 Þ : x
4649
ð27Þ
Note that the parameter a is the ratio between the radius of the circle r and 2pk, k the wavelength. Modifying (25) according to this new observable and setting the vectors
U a ¼ fU m g Ba ¼ fBm g m ¼ 1; . . . ; M and the matrix
~m þ 1Þ eiad~m d~l g l; m ¼ 1; . . . ; M: K a ¼ fð~s d
ð28Þ
The linear system to be inverted is
U a ¼ K a Ba :
ð29Þ
3.5. Ka pseudo-inverse and stability result The spectrum of the linear operator Ka and its continuous version have been carefully studied in [15]. In particular it can n o be diagonalized on the Fourier basis el ð~sÞ ¼ p1ffiffiffiffi eilh~s : 2p l2Z
K a el ¼ Dl ðaÞel ;
‘
D‘ ðaÞ ¼ 2pi ðJ ‘ ðaÞ
0 iJ‘ ð
aÞÞ;
where J‘(a) is the Bessel function of order ‘ and argument a and J 0‘ ðaÞ its derivative. Bessel functions decrease more than exponentially for ‘ above a threshold which depends on a (see Fig. 2 for an example). The operator Ka is compact and not invertible in any classical functional space, but a stable and normalized pseudo-inverse can be computed as follows:
fBm g :¼
1 ^‘ gÞ; F 1 ðfb 2LðaÞ þ 1
^‘ ¼ H‘ F ðfU m gÞ ; b ‘
ð30Þ
where F is the M discrete Fourier transform on our discrete set of directions (24). H‘ ¼ D1 if j‘j < L(a) and 0 else. ‘ L(a) corresponds to a truncation threshold. It is shown in [4] that if
LðaÞ ¼ minfa; a þ a1=3 2:5g
ð31Þ
then the norm of the normalized inverse operator defined by (30) is bounded by 1, independently of a:
supjBm j 6 supjU m j m
m
The proof is technical and based on bounding below the eigenvalues D‘(a) for ‘ < L(a). The inversion is stable under this condition. The sampling M must be larger than 2L(a) + 2. 3.6. Exact solutions for plane wave (linear) data When the data is a pure plane wave, i.e. the linearization performed in Section 3.2 is exact, an analytical form of the out~ the direction of the plane wave (h~ the corresponding angle) and A the amplitude. The put of NMLA is available. We denote d d asymptotic HF content consists in one ray going through the observation point x0 with the same direction. The observable data (26) is given by
0
Fig. 2. The stars show jJ‘ ðaÞ iJ ‘ ðaÞj (y-axis) for integer values of ‘ (x-axis) and a = 30.
4650
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 3. A constant ~~
U plwa ¼ Beiaddm ; m
1 d
~ at the observation point. curvature wavefront going through the observable geometry with ray direction d
m ¼ 1; . . . ; M:
The outputs of the NMLA procedure (30) are
8 ^plwa :¼ H‘ F fU plwa ¼ Bei‘hd~ ; >
sin ððLðaÞþ1ÞhÞ > : Bplwa ¼ BSa ðhd hm Þ m ¼ 1; . . . ; M; Sa ðhÞ ¼ ð2LðaÞþ1Þ sin2 h : m ð2Þ
ð32Þ
We see that the inverse (30) is normalized such that the Bm of maximum modulus matches the amplitude of the plane wave. The Sinc like function Sa approaches a dirac mass when L(a) increases. In practice, the discretization M must be large enough to catch the maximum of h ? Sa(hd h). As we will see in the next section this is not the case when the asymptotic wavefront has curvature or just adding noise. Then, the maximum of the peaked function is moved and a confidence interval can be estimated thanks to the stability result [4]. 3.7. Second order curvature correction Fig. 8 illustrates, in the case of a source point solution, the focusing around the ray direction produced by increasing a (and thus L(a), see (31)). However, the focusing property fails when linearization error is too large as on the last plot.1 It grows with the observation circle radius (also depending on a, see (27)). Increasing the number of samples M does not help. The linearization error main contribution comes from the second order term in the Taylor expansion (21). We did not manage to include this term in the NMLA filter nor to use the transport Eq. (14). Instead, we first assume there is only one ray (N = 1). This is not a very restrictive assumption as both time and direction windows can be used to select a particular ray (see Appendix C). We then make the simplest non linear approximation on the HF asymptotic solution:
we look for a ð1Þ sj locally constant curvature as if the wave-field solution was the fundamental solution Gðx; xs Þ ¼ 4i H0 x jxx . This virtual cðx0 Þ ~ where 1 and d ~ are respectively the sought after curvature of the front and direction source point xs is such that x0 xs ¼ dd d of the ray locally near the observation point, see Fig. 3. ~ plus x x a vector The argument x xs is decomposed into the sum of the main propagation direction ðx x0 ¼ ddÞ 0 s turning around x0: x0 xs ¼ r~s. Fig. 4 presents this decomposition. Thus, we use
U a ð~sÞ ’
sffiffiffiffiffiffiffiffiffiffiffi 8pd i ð1Þ x ~ þ r~sj : A0 ðx0 Þeixðuðx0 ÞdÞ H0 jdd cðx0 Þ 4 cðx0 Þ
ð33Þ
The high frequency amplitude and phase above are corrected as the global wave-field is not necessarily globally a point source solution in an homogeneous medium. We use the known HF asymptotic of the Hankel function, see (10). 1
See also Section 3.5.2 [5].
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4651
~ þ r~s. Fig. 4. Observable geometry data in the direction ~s on the circle is written as dd
Inspired by techniques from the Fast Multipole Method applied to electromagnetics [10,11], one can then consider r as a perturbation in the argument of the Hankel function and expand it. Using technical asymptotic expansions, we show in [4] ^ Fourier modes of the NMLA filter (30) is given by that a second order accurate approximation of the b
^‘ ’ Aeixuðx0 Þ ei b
ð‘2 14Þ
‘hd~ þ
2c
ð34Þ
;
xd is the large parameter and h is the angle associated to d. ~ This is to be compared to the exact solution for one where c ¼ cðx ~ d 0Þ plane wave linear data (32). We see that we get a c (curvature dependent) second order correction. From there we simply remark that
b^‘ 1 log i b^0
!
’ ‘hd~ þ
‘2 14 2c
ð35Þ
is a parabolic function of the mode number ‘ and the coefficients hd~ and 21c can be fitted by a least square approximation, see Fig. 5. We summarize the 2nd order correction of the NMLA algorithm: 1. Get an estimate hmI of hd~ over the set of all discretized directions (24) using (30) output:
mI
such that jBmI j ¼ maxjBm j: m
ð36Þ
Fig. 5. Stars are the (35) values and the solid line the fitted parabola as a function of ‘ (x-axis). The data correspond to the test case presented in Section 4.3.
4652
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 6. Travel-time for different windows, x-axis is the center of the picked time window, y-axis gives the travel-time computed using (40).
2. A least square parabolic fitting of
^‘ b 1 ‘ ! log ei‘hmI ^ i b0
! ’ a0 þ a1 ‘ þ a2 ‘2
ð37Þ
yields (a0, a1, a2) the constant, linear and quadratic coefficients of the parabola. Based on (35) the linear term coefficient is the angle correction: a1 ¼ dh ¼ hmI hd (notice the hmI phase shift in (37)) and the quadratic term coefficient depends on the mean curvature of the front: a2 ¼ 21c. 3. Correct NMLA amplitudes
0
1 ð‘2 14Þ 0 1 @ ^ i ‘dh 2c A: b‘ e Bm :¼ F
ð38Þ
Fig. 9 illustrates the improved focusing and increase amplitude at the maximum.
3.8. Travel-time estimate We mentioned in Section 2.4 the possibility of preselecting a time window around a hand picked time t0. The travel-time of the ray corresponding to the wave recorded in this time window is necessarily close to t0. The wavelet is however not a perfect wavefront and may spread over a given period of time (around 1 s in our numerical illustration Fig. 17).
Fig. 7. Seismogram with time (y-axis) and position (x-axis) picking.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4653
NMLA combined with the use of time data can provide an accurate travel-time estimation. The idea is to use the implicit redundancy of the HF travel-time over the frequency range. More precisely we have, again assuming perfect asymptotic one ray ansatz data,
Fig. 8. One source point: jBmj (y-axis) versus hm (x-axis). Varying a from top to bottom: 10, 20, 40. Vertical line is the exact ray angle. Obtained from 20% noisy data.
4654
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
BmI ðxÞ ’ A0 expix/ðx0 Þ :
ð39Þ
The x dependence is explicit in the notation since we are going to use the variation with respect to frequency,
/ðx0 Þ ’
@ x BmI : iBmI
ð40Þ
The x derivative is obtained by applying both procedures, NMLA and second order curvature correction. to the x derivative of the data. The derivative data can be computed using finite difference or using the Fourier transform of t u(t, .). Using the case study of Section 4.3, we can show that this method is more accurate than ‘‘hand picking’’ t0: we recall that we use time data only in a time window centered around this ‘‘hand picked’’ t0. The width of the window must be larger than the width of the source wavelet. In Fig. 6, we plot the travel-time estimated by (40) against t0. As long as the time window contains significant part of the signal, the travel-time is constant and correct and therefore does not depend on the picking time. 4. Numerical illustration: source point discovery 4.1. Generation of synthetic data We generate synthetic data using the following procedure: A standard FDTD scheme, second order in time and space, is used to simulate the waves generated by one or more source points. Discretization in space is chosen to maintain a minimum of 20 points per wavelength for the cutting frequency Fc = 3 xs ⁄(space step is h ¼ mincðxÞ ), in our test case Fc = 90 Hz. Time discretization satisfy the usual CFL condition (time step Fc is Dt ¼ pffiffi h ). 2maxcðxÞ
A ‘‘Ricker’’ (Gaussian derivative) is the source wavelet. The source pulsation xs (see Section 2.3) we use is 60p rad/s. We use a perfect reflecting boundary condition @@m u ¼ 0 at the bottom of the domain and Higher Order Absorbing Conditions [13] elsewhere. We record over time a seismogram (just u) at all discretization points on a line. Once the observation time t0 and surface position x0 have been picked (see Fig. 17), a smooth window in space and time is used to single out one event. We Fourier transform all the time lines, then choose a frequency x such that the frequency content of the source is significant. Typically close to xs. The NMLA frequency data is then computed on the observation circle using the one way method described in Appendix A. 4.2. Homogeneous medium 4.2.1. one source point Fig. 7 depicts the seismogram recorded over time (vertical axis). Each vertical line is the signal t ? u(t, x) recorded at a grid point x on the surface (horizontal axis). The horizontal deviation of these lines is set proportional to u. It generates oscilla-
Fig. 9. One source point: jBmj (y-axis) versus hm (x-axis). Solid line: NMLA - Solid line with + markers: NMLA + curvature correction - Vertical line: exact ray direction.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4655
tions and a shade at times when waves (also called events) reach the surface. In this particular simulation there is only one source point but two successive events due to the reflecting boundary condition at the bottom. We manually select an observation point on the surface x0 (vertical line). We also manually pick a time t0 (horizontal line) which approximately gives the travel-time of the event. Using a time window centered at t0, we can eliminate other events and retain only one significant wavefront and ray in the solution. Travel-time estimation is also discussed in Section 3.8. Fig. 8 shows the impact of a, the radius of the observation circle counted in 2p wavelength, on the basic NMLA filter (30). If a is too small, we do not have enough Fourier modes to focus, note the angle is still correctly estimated for one ray, however for more than one ray (as in the next section) a = 10 would not allow to observe cleanly separated maxima when the ray directions are close. If a is too large, then the curvature error induced by the plane wave approximation pollutes the result. Fig. 9 shows the output of the NMLA procedure (30) and also the second order correction of Section 3.7. After correction the amplitudes are better focused and the maximum more accurate. Table 1 quantifies the error on the angle and the source point discovery. Local Radon transform (LSS) completely fails while NMLA does a better job than plane wave destructor (PWD) (references on these two methods in Section 4.3). Table 1 HF ray direction analysis and source point discovery. Method
Ray angle error
NMLA basic NMLA 2nd order LSS PWD
0.5988° 0.0184° failed 0.1917°
Error (distance to source point) 0.0011 0.0071
Fig. 10. One source point (0.1, 0.2): Angle error (solid line) and Travel-time error (dashed line) in Log10 scale versus the position of the observation point x0 on the ‘‘surface’’ xo2 ¼ 1.
Fig. 11. Seismogram with time (y-axis) and position (x-axis) picking.
4656
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 12. Two source points: jBmj (y-axis) versus hm (x-axis). Dashed line is NMLA output for 2 source points. Solid line is NMLA + Gaussian window to enable separation a = 2.5 (more details in Appendix C).
Fig. 13. jBmj (y-axis) versus hm (x-axis). White noise: 20% solid line - 40% dashed line.
Fig. 14. jBmj (y-axis) versus hm (x-axis). Correlated noise: 20% solid line - 40% dashed line.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 15. Heterogeneous medium: velocity distribution.
Fig. 16. Synthetic data numerical simulation, snapshots at 4 successive times.
4657
4658
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 17. Seismogram with time (y-axis) and position (x-axis) picking.
Fig. 10 shows the error on angle and travel-time when applying NMLA second order to several position on the surface y = 1. 4.2.2. Two source points Basic NMLA is able to recover well separated multiple rays. We treat here the case of two but it also works fine for more, see the original paper [5,15]. Fig. 11 shows the recorded seismogram for two source points. The picked observation point is at equal travel-time from the source points. Thus, the time window will contain two wavefronts and two rays. Fig. 12 shows the NMLA basic filter result on the left and the application of the Gaussian window convolution (47) (see Appendix C). The filter removes the possible interferences caused by the trailing oscillations the ray contributions to B can be considered separately. One can go back ^ of each ray contribution and perform the curvature correction and traveltime estimation (Sections 3.7 to the fourier modes b and 3.8). Table 2 HF ray direction analysis and source point discovery. Method
Ray angle
Error (distance to source point)
NMLA basic NMLA 2nd order LSS PWD
90.2970° 90.5503° 91.002° 90.7418°
0.0069 0.0878 0.0201
Fig. 18. Ray backward (downward) propagation with initial conditions specified in Table 2. The circle is centered at the true source position. Source discovered with, from left to right, NMLA/PWD/LSS ray angle at observation.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4659
4.2.3. Homogeneous medium noisy data - stability In this section we present numerical results (Figs. 13 and 14) which illustrate the stability claim on the algorithm (Section 3.5) in the two source points case. The NMLA filter (30) easily filters out white noise (see Fig. 13) as it is precisely designed to single out coherent signals. In order to test correlated noise, we perturb the data with the solution collected on a circle around the observation point but with a different radius (see Fig. 14). 4.3. Heterogeneous medium We now consider the smooth speed velocity distribution c(x) of Fig. 15. The + sign is the source point location and our observation point is the cross at the surface. Fig. 16 shows several snapshots of the FDTD simulation used to generate the synthetic seismogram in Fig. 17. Table 2 displays the ray direction found with several methods. NMLA basic and the second order correction. The Plane Wave Destructor Method PWD [21] and the Local Slant Stack LSS (Radon transform) [8,19] are popular methods in Geophysics and radar processing but are sensitive to the necessary windowing. Note that PWD only identified plane waves by construction; the Radon transform can be generalized to identify parabolas [6] but we did not tested it. Note that both PWD and LSS can be considered as image processing methods and do not use the wave models. Fig. 18 illustrates the source discovery through the proposed HF TR method (Section 2.4). NMLA second order yields a u(x0) = 0.4602 s. travel-time. We use this information together with the three different ray directions found in Table 2 to shoot the rays backwards. A circle centered around the true source position highlights the better performance of NMLA. 5. Conclusion The second order version of NMLA method is a stable local HF wave analysis tool. It is based on the true wave HF models and robust even to correlated noise. The companion paper [4] also tries to mathematically establish how close in the near field this method is going to work. It is completely automatic and has no other tuning parameters than a, the radius of the observation circle measured in 2p wavelength. A mathematical study of the error depending on a and the discretization parameter M is done in [4]. It also proposes an optimal law for a depending on the front curvature. The source deconvolution is only needed for the travel-time estimation. The curvature correction requires that only one ray is present in the data. We briefly explain in Appendix C how a gaussian filter can be use to separate multiple ray components with sufficiently distinct directions. This part should be more thoroughly investigated in a forthcoming paper. A numerical example close to a caustic was presented in [5]. Even though rays directions were recovered, this topic certainly deserves a more serious investigation. We think that NMLA can be extended to 3D and other linear models including elastic waves (see [15]). Appendix A. Seismogram data/ FT/ one-way extrapolation of linear array data We assume for simplicity that:
We are in 2D and x = (x1, x2). ^ 0 ðx1 ; xÞ is given on the surface x2 ¼ xo2 . The data u Our observation point is on the data surface x0 ¼ xo1 ; xo2 . The medium speed is constant c(x) = v in a layer containing the observation circle of radius r (see (26)).
We recall here the one way extrapolation procedure. More details, references and also the more complicated case of one way extrapolation in heterogeneous medium, can be found in [14]. First, we Fourier transform the data and Eq. (13) with respect to (x1, x2) and denote the Fourier space variable (k1, k2). Fourier transform rules that the partial derivatives @ xj corresponds in Fourier space to a multiplication by, respectively, i kj, that is
x2 ^ 2 2 ^ 0 ðk1 ; xÞ given on x2 ¼ xo2 : ðk1 ; k2 ; xÞ ¼ 0 u k1 þ k2 2 u
v
ð41Þ
^ notation and drop the x dependence for simplicity. If we restrict ourselves to real ‘‘up-going’’ k2 > 0 Fourier We keep the u mode, this is the one way assumption, and after factorization, Eq. (41) can be simplified into
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!
k2
x2 c2
2 ^ ðk1 ; k2 Þ ¼ 0: k1 u
ð42Þ
The significant real Fourier modes of the solutions have to satisfy the (so called) ‘‘dispersion’’ relation characterizing the paraxial operator:
4660
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x2
k2 ¼
c2
2
k1 :
ð43Þ
The solution can be computed as a function of (k1, x2). Indeed after an inverse Fourier transform in x2, we get the ODE (we ^ notation again) abusively keep the u
^ k1 ; xo2 ¼ u ^ ðk1 ; x2 Þ ¼ k2 u ^ ðk1 ; x2 Þ u ^ 0 ðk1 Þ; i@ x2 u
ð44Þ
where k2 is given by (43), with trivial solution The exact solution of (44) at depth x2 is
^ 0 ðx; k1 Þeik2 x2 : ^ ðk1 ; x2 Þ ¼ u u The solution at any point (x1, x2) is recovered by inverse Fourier transform from k1 to x1. This procedure generates the data on a circular array for up-going rays in constant medium at the price of two 1D Fourier transforms. It could also be used to filter out down-going rays. As we restrict to observation point near the x2 ¼ xo2 depth of the linear array, it is possible to use a x1 windows which size depends on the finite speed v and the resulting dependence cone. Appendix B. When the wavelet W xs is unknown In that case, the source deconvolution (20) is not possible so we just remove the half primitive factor
^ ðx; xÞ u
pffiffiffiffiffiffiffiffiffiffi ^ðx; xÞ: ixu
ð45Þ
This implies that the NMLA algorithm expected output is multiplied by the frequency domain wavelet
B
c x ðxÞ: BW s
ð46Þ
The curvature correction method of Section 3.8 is unchanged since the amplitudes are removed in (37). The travel-time estimation method (40) does not work since the amplitude BmI contains now the wavelet with an unknown x dependent phase. Appendix C. Time and ray direction windows As shown in Fig. 17 we assume that, after picking the observation point x0 ¼ xo1 ; xo2 on the surface data (the vertical red2 line), we select an observation time in the seismogram (the horizontal red line). This can be done by the user or an ad hoc automatic signal processing procedure. The advantage of time picking, is that rays crossing at the observation point with different travel-time may be separated using a time window. In Fig. 17 we see that a second later wavefront passes the observation point. If we just Fourier transform in time and then perform NMLA, we expect the (Bm)’s to characterize a two ray ansatz. Using a smooth time window, we are able to isolate for instance only the first arrival. We then only have one ray and can perform the curvature correction as explained in Section 3.7. However multiple rays with same or close travel-times may cross at the observation point. Then, the time window is unable reduce the data to only one ray (N = 1). The NMLA filter output is the sum of cardinal sine Sa functions centered at the different ray directions. Even if the rays directions are well separated interferences make it impossible to locate ray directions as the maximum of the two components. A way to fix this problem is to remove the trailing oscillations of Sa by a convolution of B with a smooth Gaussian weight function for instance. The convolution does not affect the Sa components maximum position. The convolution is simply done using an additional multiplicative filter in Fourier space before the last Fourier transform in (30):
^‘ b
^‘ w ^‘ b
^ ‘g ¼ F fw
n o
a 2 exp2hm :
ð47Þ
The parameter a can be used to tune the Gaussian filter and determines the resolution of the filter or the minimum width between ray direction in order to separate them. A large a yields a small window but a must be small enough to remove the high frequency oscillations. We heuristically found that a = 2.5 seems optimal for our test cases. Once ray components have been separated one can go back to the Fourier modes and perform the curvature correction (Section 3.7). Fig. 12 shows the filter effect comparison for two source points in homogeneous medium. References [1] T. Abboud, J.-C. Nédélec, B. Zhou, Méthode des équations intégrales pour les hautes fréquences, C.R. Acad. Sci. Paris Sér. I Math. 318 (2) (1994) 165–170. [2] M. Balabane, V. Tirel, Décomposition de domaine pour un calcul hybride de l’équation de Helmholtz, C.R. Acad. Sci. Paris Sér. I Math. 324 (3) (1997) 281–286. 2
For interpretation of colour in Fig. 17, the reader is referred to the web version of this article.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4661
[3] Jean-David Benamou, An introduction to Eulerian geometrical optics (1992–2002), J. Sci. Comput. 19 (1–3) (2003) 63–93 (Special issue in honor of the sixtieth birthday of Stanley Osher). [4] Jean-David Benamou, Francis Collino, and Simon Marmorat. Second order stable numerical microlocal analysis of harmonic plane wave and source point wavefields, submitted for publication.
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Source point discovery through high frequency asymptotic time reversal Jean-David Benamou a,⇑, Francis Collino b, Simon Marmorat a a b
INRIA, B.P. 105, 78153 Le Chesnay Cédex, France CERFACS, 42 Avenue Gaspard Coriolis, 31100 Toulouse, France
a r t i c l e
i n f o
Article history: Received 18 January 2011 Received in revised form 2 March 2012 Accepted 3 March 2012 Available online 4 April 2012 Keywords: Linear waves Geometrical optics Time reversal
a b s t r a c t Given local frequency domain wave data, the Numerical Micro-Local Analysis (NMLA) method (Benamou et al., 2004) [5] and its recent improved version (Benamou et al., 2011) [4] gives a pointwise numerical approximation of the number of rays, their slowness vectors and corresponding wavefront curvatures. With time domain wave data and assuming the source wavelet is given, the method also estimates the travel-time. The paper provides a non technical presentation of the improved NMLA algorithm and presents a numerical application which can be interpreted as a high frequency asymptotic version of the classical time reversal method (Borcea et al., 2003) [7]. A detailed technical presentation of the algorithm is available in Benamou et al. (2011) [4] and more numerical experiments can be found in Collino and Marmorat (2011) [15]. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction Source point or scatterer location identification is a classical inverse problem in computational wave propagation. It is of great practical interest in, for example, non destructive testing of alloys defect, kidney stones location discovery [28], underwater communications [23]. It is a basic step in inverse synthetic aperture radar [8]. See also [24] for the same problem of source discovery. The popular time reversal (TR) method (see [7] for a mathematical study and references) solves this inverse problem. It is rapidly summarized in the following steps: (i) record a wave signal on some array of ‘‘transducers’’ for ‘‘some’’ time, (ii) time reverse this signal and re-emit it in the medium. Then, physics observes and maths proves that the time reversed signal will exhibit refocusing properties towards the scatterers or the source points depending on the set up of the experiment. We propose a time reversal method based on the High Frequency (HF) asymptotic model of wave propagation, also known as Geometric Optics (GO). The HF model only makes sense with smooth deterministic media when the characteristic wavelength of the signal is smaller than the variation scale of the acoustic speed. We will therefore NOT be concerned with the widely studied ‘‘super-resolution’’ effect of the focusing observed in random media [7]. The goal of the paper is twofold: Provide a non technical description of the Numerical Micro-Local Analysis algorithm (NMLA) [5] and its recent improved second order version [4]. We apply it to analyze the underlying GO ray directions of the signal recorded in step (i). The ⇑ Corresponding author. E-mail address: [email protected] (J.-D. Benamou). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.03.012
4644
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
NMLA algorithm works in the frequency domain. We therefore use classic pre and post processing techniques to extend it to time domain and also to treat a signal recorded on linear arrays of transducers. We show that the use of time domain data simplify multiple rays discrimination and that frequency redundant information provides a travel-time estimate. Illustrate numerically that the output ray direction and travel-time are accurate enough, at least with synthetic noisy data, to shoot the ray back to its unknown source location. The notion of ‘‘shooting back’’ implies that the travel-time, the actual propagation time, is reversed in the same way as in step (ii) above. Finally, we compare the ray directions to those found with more classical signal processing methods like Plane Wave Destructor [21] or the Local Slant Stack methods [8,19]. Potential applications of the NMLA method are of course not limited to source point discovery. Here is a non exhaustive list: Micro-local finite element methods try to discretize the geometric optic quantities [1,17]. If known, an analytic solution for the phase / can be use to remove oscillations and the size of the discretization [12,25,22]. GO is simple and cheap when the underlying wave field is not too complicated: no diffraction, reflections, caustics, etc. In these more complicated regime, GO has to manage several complex interacting ansatz. As simple propagating waves can be treated locally by GO, coupling GO and full wave fields across physical interfaces or in a multi-scale framework make sense [26,2]. A reliable process to invert the GO ansatz is needed. In Geophysics ‘‘Cleaning’’ data by filtering out low frequencies either to use faster HF models or simply because only high frequencies are really relevant (migration) or nicer for inversion (see [30,29] for example). The Direction of Arrival problem (DOA) is a classic signal processing problem [27] where one tries to find the direction of the source of a signal impinging on a, generally linear, array of antennas. Source point reduced models for far field representations [9]. In order to make the paper self contained, it is organized as follows: Section 2 recalls the HF model derivation and formulates the HF asymptotic time reversal source discovery inverse problem. Section 3 reviews the refined NMLA algorithm [4] and its extension to time domain data recorded on a linear array of ‘‘transducers’’. Section 4 presents numerical source discovery from synthetic data. 2. Forward and inverse problem 2.1. High frequency acoustic wave equation We start from the classical acoustic wave equation in the time domain. The wave field u(x, t) is generated by a finite time wavelet emitting from a source point:
8 < c21ðxÞ @ 2tt uðx; tÞ Duðx; tÞ ¼ dðx xs ÞW xs ðtÞ; 8ðx; tÞ 2 R2 ½0; TÞ; :
uðx; t ¼ 0Þ ¼ @ t uðx; t ¼ 0Þ ¼ 0;
ð1Þ
8x 2 R2 ;
where
[0, T) is the simulation time (one may consider T = +1). x ´ c(x) is the acoustic speed law function on R2 . As we deal with HF–GO models, we will consider smooth speeds. xs is the source point location. t # W xs ðtÞ is the source time wavelet. The radian frequency x represents the dual variable of time in the Fourier transform denoted^: (we assume u 0 for t < 0):
^ ðx; xÞ ¼ u
Z
1
uðx; tÞeixt dt:
ð2Þ
0
The source wavelet, denoted W xs ðtÞ, determines the frequency content of the solution. It depends explicitely on the ‘‘characteristic pulsation (or frequency)’’ xs. We define it as
W xs ðtÞ ¼ Wðxs tÞ; where W(t) is a shape function with support in [0, 1]. Assuming further that W is absolutely integrable and has zero mean, the following properties are satisfied:
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
1 ; SupportðW xs Þ ¼ 0; xs c x : c x ð xÞ ¼ 1 W W s For
ð3Þ ð4Þ
xs
xs
> 0, there exist two positive constants cþ Z
xRI
c x ðxÞj2 dx 6 jW s
Z
4645
c x ðxÞj2 dx; jW s
and c such that
I ¼ ½c xs ; cþ xs :
R
ð5Þ
So, when xs tends to infinity, (3) shows that the source time support shrinks and (4 and 5) that the frequency content grows and but is shifted towards infinity. 2.2. Frequency domain and ansatz The HF model uses a formal asymptotic representation of the solution in which the unknowns are independent of the radian frequency x. It is easier to derive and interpret it in the frequency domain: Assume T = +1 and Fourier transform equation (1):
x2 ^ c x ðxÞ: ^ ðx; xÞ ¼ dðx xs Þ W uðx; xÞ Du s 2
c ðxÞ
ð6Þ
We also assume that c(x) = v is constant outside a compact in X and add the Sommerfeld radiation boundary condition at infinity
^ ðx; xÞ lim @ jxj u
jxj!1
ix
v
^ ðx; xÞ ¼ 0: u
ð7Þ
Problems (6) and (7) are well posed [16]. In the homogeneous case, c(x) = v constant everywhere, the analytic solution of (6) and (7) is given by
^ ðx; xÞ ¼ u
i ð1Þ jx xs j c W xs ðxÞ; H0 x 4 v
ð8Þ
ð1Þ
where H0 ðuÞ is the Hankel function of first type and order 0. Plugging in (8) the following asymptotic estimate (cf. [31])
sffiffiffiffiffiffiffi 2 iðqp4Þ ; qÞ e
ð1Þ H0 ð
q!1
pq
we obtain an asymptotic approximation of the solution in the form:
c ð xÞ W s ^ ðx; xÞ ’ Aðx; xs Þeix/ðx;xs Þ pxffiffiffiffiffiffiffiffiffi ffi ; u ix
ð9Þ
where
Aðx; yÞ ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
v
8pjx yj
/ðx; yÞ ¼
;
jx yj
v
:
ð10Þ
sj The large parameter is x jxx v . It represents the number of wavelengths between the source and the observation point. In the time domain (9) corresponds to the ‘‘ray’’ solution
1 Aðx; xs Þ uray ðx; tÞ ’ pffiffiffiffiffiffi W 2 ðxs ð/ðx; xs Þ tÞÞ;
xs
1
where W 2 ðtÞ is
ð11Þ
RW b ðxÞ pffiffiffiffiffiffi eixt dx, the half derivative of the wavelet. ix
The ‘‘ansatz’’ idea, apparently going back to a 1911 paper of Debye on ‘‘the article of Sommerfeld’’ (see [20]), consists in generalizing the a priori form (9) to the general heterogeneous case and search for a formal asymptotic solution of the Helmholtz equation in the form
^ ray ðx; xÞ ’ eixuðxÞ u
1 A ðxÞ P j
xÞ j
j¼0 ði
:
ð12Þ
pffiffiffiffiffiffiffiffiffiffi We assume the source wavelet and the half primitive (symbol ix) have been removed, see Section 3.1. In general, only the first A0 term of the amplitude series is used. Oscillations are formally enforced in the phase and travel at travel-time u. The radian frequency x is the large parameter.
4646
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
2.3. Geometric optics–Ray Tracing We give a quick overview to make the paper self contained, for more details and references see [20]. Assuming a solution taking the ansatz form (12) locally away from the source. One can plug it into
^ ðx; xÞ þ Du
x2 ^ uðx; xÞ ¼ 0
c2 ðxÞ
ð13Þ
and expand in inverse power of x. Canceling the coefficients of the first and second order term of the series gives respectively the Eikonal equation and amplitude equations
(
1 jruðxÞj ¼ cðxÞ ;
2ruðxÞ rA0 ðxÞ þ A0 ðxÞDuðxÞ ¼ 0:
ð14Þ
These equations are local but global solutions are constructed using a Lagrangian method called Ray Tracing. The travel-time u is computed along the integral curves of the gradient called rays and denoted s ? y(s, y0). Each ray may be parameterized by its initial position y0 and, assuming the travel-time and its gradient are known at this location, one solves the following system of ODEs:
8d > < ds yðs; y0 Þ ¼ ruðyðs; y0 ÞÞ :¼ pðs; y0 Þ; yð0; y0 Þ ¼ y0 ; d pðs; y0 Þ ¼ 12 rn2 ðyðs; y0 ÞÞ; pð0; y0 Þ ¼ ruðy0 Þ given; ds > :d 2 u ðyðs; y ÞÞ ¼ n ðyðs; y ÞÞ; uðy0 Þ given; 0 0 ds
ð15Þ
where y(s, y0) is the ray shot from y0 and s is a parameterization of the curve. p = r/(y) is called the ‘‘slowness’’ vector and gives the ray direction. n ¼ 1c is called the slowness of the medium. When the velocity is homogeneous (c(x) = v) around a source point, the analytic asymptotic form (10) can be used to define the initial conditions: For s small enough the level curve W C ¼ fy0 ; s:t:/ðy0 Þ ¼ jy0 vxs j ¼ sg is a circle around the source sÞ . Letting s go to 0 we consider that the HF asymptotic form of a source point corresponds point and pð0; y0 Þ ¼ ruðy0 Þ ¼ vðyky0 x xs j 0
to rays leaving the source point xs isotropically with the same amplitude. We then switch the parameterization of the rays to the initial take off angle y0 ? h0. The global HF solution is formed of all rays; this is the solution, for all h0 2 [0, 2p[, of
8d > < ds yðs; h0 Þ ¼ ruðyðs; h0 ÞÞ :¼ pðs; h0 Þ; yð0; h0 Þ ¼ xs ; ~ d pðs; h0 Þ ¼ 12 rn2 ðyðs; h0 ÞÞ; pð0; h0 Þ ¼ sv0 ; ds > :d uðyðs; h0 ÞÞ ¼ n2 ðyðs; h0 ÞÞ; uðxs Þ ¼ 0; ds
ð16Þ
where ~s0 ¼ ðcos h0 ; sin h0 Þ. Ray Tracing (16) can be extended to compute solutions of the amplitude equation (see [3] for instance). In smooth heterogeneous media, the Lagrangian travel-time function remains smooth but it is quite common for rays to cross. Then, solutions to the Eulerian equation (14) become multi-valued. In the simplest case, the HF ansatz remains a superposition/sum of elementary (12) ansatz. See [18,3] for more details and references on multi-valued solutions. The NMLA algorithm assume that its data, the Helmholtz solution, can be locally represented as a general HF ansatz of the form
^ ðx; xÞ ’ u
N P
A0;n ðxÞeixun ðxÞ ;
ð17Þ
n¼1
where N is the number of rays going through x and n the subscript labeling the different branches of the solution. We kept only the zeroth order term of the amplitude series of (12). 2.4. Source point discovery with HF time reversal: formulation and method of resolution In order to simplify the presentation of the source point discovery and resolution method, we will restrict to the case N = 1. We give NMLA computations for N = 2 in Section 4.2.2 and discuss ray separation using time and direction windows in Appendix C. The source point discovery problem can be formally formulated as follows: Let u be the solution of (1). The source point location xs is the unknown we are looking for. The wavelet W is assumed to be given (see Appendix B for a discussion on the method when it is also unknown). We fix an observation point x0 ¼ xo1 ; xo2 such that c(x) = v is constant for all x ¼ ðx1 ; x2 Þ; x2 > xo2 . We assume further that there is a time t0 such that u can be represented, in a neighborhood of (t0, x0) as an asymptotic ray solution in the sense of (11):
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661 1 AðxÞ uðx; tÞ ’ pffiffiffiffiffiffi W 2 ðxs ðuðxÞ tÞÞ þ duðx; tÞ:
xs
4647
ð18Þ
In (18), we have defined two new auxiliary unknowns A and u which satisfy the GO equations (14) and the Ray Tracing equations (16). We assume the remainder du and all its derivatives to be negligible with respect to u as xs ? +1. Our method to locate xs first uses the local data (18) to compute the ray unknowns u(x0) and pdata ¼ rudata ðx0 Þ. Then, we 0 use the global reversible rays equations to back propagate the ray to the source point from these initial conditions. This is accomplished by changing ds in ds (16) and use the initial conditions obtained in (i) for this particular ray:
8 d > < ds yðsÞ ¼ ruðyðsÞÞ :¼ pðsÞ; yð0Þ ¼ x0 ; d pðsÞ ¼ 12 rn2 ðyðsÞÞ; pð0Þ ¼ pdata 0 ; ds > : d 2 uðyðsÞÞ ¼ n ðyðs; h0 ÞÞ; uðxs Þ ¼ udata ðx0 Þ: ds
ð19Þ
The travel-time u(y(s)) will then decrease and we stop integrating (19) when it gets to 0. We have ‘‘time reverted’’ the ray and y(s) has reached the source xs. We emphasize that the medium slowness n ¼ 1c must be known. The ray unknowns u(x0) and pdata ¼ rudata ðx0 Þ are computed using the NMLA algorithm. This is our main contribution 0 ^ ðx; xÞ on a local observation circle around the choand it is described in Section 3. The method uses frequency domain data u sen observation point x0. Common numerical discretizations of wave fields directly do not provide data on a circular array around an ‘‘observation’’ point. Regarding real data, circular arrays do not seem to be a common acquisition device either. Interpolation is a solution but it needs data on a volume covering the acquisition circle which typically is several wavelengths large. Considering the most common acquisition format is a linear array, like for instance seismograms recorded on an ideal flat surface in geophysics, we propose here to use one way wave extrapolation to reconstruct the circular array from a seismogram from this type of data. It implies that rays of interest must be transverse to the line of acquisition but this is not very restrictive. We summarize below the data pre-processing steps to generate NMLA compliant data: The acoustic signal (18) is recorded on a line or a linear array which contains a chosen observation point x0 ¼ xo1 ; xo2 . In o 2D, it gives time data u x1 ; x2 ; t for t in some time interval [0, T] and x1 on the line (we denote x = (x1, x2)). As we deal with time data it is possible to preselect a given time window centered around a chosen time of interest t0 by ~ x1 ; xo2 ; t ¼ u x1 ; xo2 ; t gðt t0 Þ were g is a smooth cut-off window with support larger than the source waveconsidering u let’s W xs . ^ ~ x1 ; xo2 ; x ¼ F t!x u ~ x1 ; xo2 ; t to go to frequency domain and choose a frequency x Next we Fourier transform the data: u such that the frequency content of the source is significant, typically close to xs. ^ ~ x1 ; xo2 ; x to reconstruct the field and its gradient in a neighborhood of x0, (see We use a paraxial extrapolation in x2 of u Appendix A). Remark 1. All the procedure only makes sense if the HF representation (18) is relevant. If not it simply means that the HF model is not suitable to approximate the solution, i.e. there are no rays in the solution. Remark 2. Classical antenna localization algorithms (GPS, cell phone) use estimated travel-times from different sources with known locations (satellites) to draw spheres which intersections ‘‘triangulate’’ the antenna position. The wave signal directions between the antenna and the sources are not known. Our approach is based on finding both the direction and the travel-time from a single source and then back propagate ‘‘exactly’’. Repeating the resolution method at different observations points and using different rays is of course also possible. 3. Numerical micro-local analysis NMLA does not directly invert (17) but a system of linear equations obtained through the successive steps described in the following subsections. 3.1. Source deconvolution When the wavelet W xs is known, one can perform the following source deconvolution on the data
^ ðx; xÞ u
pffiffiffiffiffiffiffiffiffiffi ix ^ ðx; xÞ: u c x ðxÞ W
ð20Þ
s
When the source wavelet is unknown, the method still works but the travel-time cannot be estimated, see Appendix B. 3.2. Plane wave approximation In the GO regime, we expect the phase to be real and with variations on a scale larger than x. Within a few wavelength of x0 the constant amplitude approximation A0(x) ’ A0(x0) will remain of order x1 .
4648
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Using a Taylor expansion on the travel-time and in the vicinity of x0 we have
1 2
uðxÞ ’ uðx0 Þ þ ðx x0 Þ ruðx0 Þ þ ðx x0 ÞT Huðx0 Þðx x0 Þ þ
ð21Þ
The usual first order ‘‘plane wave’’ approximation is, for N = 1 in (17),
^ ðx; xÞ ’ Bðx0 Þeixðxx0 Þruðx0 Þ ; u
ð22Þ
where
Bðx0 Þ ¼ A0 ðx0 Þeixuðx0 Þ : ^ can be identified as the N-rays ansatz In the more general case where N rays cross at x0, we assume u
^ ðx; xÞ ’ u
N P
Bn ðx0 Þeixðxx0 Þrun ðx0 Þ :
ð23Þ
n¼1
3.3. Relaxation Our goal is to identify (Bn, run) and N the number of rays which appear in Eq. (23). In order to simplify this non-linear inverse problem we instead consider the linear problem of finding the amplitude distribution.{Bm, m = 1, . . . , M} over a dis~m ¼ m ¼ 1; . . . ; Mg with cretization of the set of 2-D unit vectors, fd
~m ¼ ðcos hm ; sin hm Þ hm ¼ m 2p ; d M
m ¼ 1; . . . ; M
ð24Þ
The ~: notation will always indicate unit vectors. As the norm of the gradient of the travel-times is given by the Eikonal equation,
jrun ðx0 Þj ¼
1 cðx0 Þ
for all n;
the relaxed version of (23) is
^ ðx; xÞ ’ u
M P
~m i x ðxx0 Þd
Bm ðx0 Þe cðx0 Þ
ð25Þ
:
m¼1
The size of the discretization M is discussed below. 3.4. The observable and the linear system ^ ðx; xÞ on a scaled circle around x0. The discretization (24) was used to construct a In our first paper [5] we evaluated u finite dimensional linear system (see the geometry os the support on Fig. 1). It was subsequently inverted to obtain the amplitude distribution Bm. The data for the linear system was given by
~ ; xÞ r ¼ acðx0 Þ ^ ðx0 þ r d U a;m ¼ u m
x
j~sj ¼ 1:
ð26Þ
This choice sometimes led to unstable inversions and a L2 regularization was used in [5]; a different L1 regularization was proposed in [24]. The new version proposed in [4] still uses the circle geometry but the sampled data is in impedance form: n o cðx0 Þ @ u ^ ^ . The resulting NMLA filter is stable independently of its parameters. The new observable therefore is: þu ix @r
Fig. 1. Geometry of the observable support, a sample point in the direction ~s. r ¼ a cðxx0 Þ.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Um ¼
^ cðx0 Þ @ u ~ m ; xÞ þ u ~m ; xÞ; ^ ðx0 þ r d ðx0 þ rd ix @r
r¼
acðx0 Þ : x
4649
ð27Þ
Note that the parameter a is the ratio between the radius of the circle r and 2pk, k the wavelength. Modifying (25) according to this new observable and setting the vectors
U a ¼ fU m g Ba ¼ fBm g m ¼ 1; . . . ; M and the matrix
~m þ 1Þ eiad~m d~l g l; m ¼ 1; . . . ; M: K a ¼ fð~s d
ð28Þ
The linear system to be inverted is
U a ¼ K a Ba :
ð29Þ
3.5. Ka pseudo-inverse and stability result The spectrum of the linear operator Ka and its continuous version have been carefully studied in [15]. In particular it can n o be diagonalized on the Fourier basis el ð~sÞ ¼ p1ffiffiffiffi eilh~s : 2p l2Z
K a el ¼ Dl ðaÞel ;
‘
D‘ ðaÞ ¼ 2pi ðJ ‘ ðaÞ
0 iJ‘ ð
aÞÞ;
where J‘(a) is the Bessel function of order ‘ and argument a and J 0‘ ðaÞ its derivative. Bessel functions decrease more than exponentially for ‘ above a threshold which depends on a (see Fig. 2 for an example). The operator Ka is compact and not invertible in any classical functional space, but a stable and normalized pseudo-inverse can be computed as follows:
fBm g :¼
1 ^‘ gÞ; F 1 ðfb 2LðaÞ þ 1
^‘ ¼ H‘ F ðfU m gÞ ; b ‘
ð30Þ
where F is the M discrete Fourier transform on our discrete set of directions (24). H‘ ¼ D1 if j‘j < L(a) and 0 else. ‘ L(a) corresponds to a truncation threshold. It is shown in [4] that if
LðaÞ ¼ minfa; a þ a1=3 2:5g
ð31Þ
then the norm of the normalized inverse operator defined by (30) is bounded by 1, independently of a:
supjBm j 6 supjU m j m
m
The proof is technical and based on bounding below the eigenvalues D‘(a) for ‘ < L(a). The inversion is stable under this condition. The sampling M must be larger than 2L(a) + 2. 3.6. Exact solutions for plane wave (linear) data When the data is a pure plane wave, i.e. the linearization performed in Section 3.2 is exact, an analytical form of the out~ the direction of the plane wave (h~ the corresponding angle) and A the amplitude. The put of NMLA is available. We denote d d asymptotic HF content consists in one ray going through the observation point x0 with the same direction. The observable data (26) is given by
0
Fig. 2. The stars show jJ‘ ðaÞ iJ ‘ ðaÞj (y-axis) for integer values of ‘ (x-axis) and a = 30.
4650
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 3. A constant ~~
U plwa ¼ Beiaddm ; m
1 d
~ at the observation point. curvature wavefront going through the observable geometry with ray direction d
m ¼ 1; . . . ; M:
The outputs of the NMLA procedure (30) are
8 ^plwa :¼ H‘ F fU plwa ¼ Bei‘hd~ ; >
sin ððLðaÞþ1ÞhÞ > : Bplwa ¼ BSa ðhd hm Þ m ¼ 1; . . . ; M; Sa ðhÞ ¼ ð2LðaÞþ1Þ sin2 h : m ð2Þ
ð32Þ
We see that the inverse (30) is normalized such that the Bm of maximum modulus matches the amplitude of the plane wave. The Sinc like function Sa approaches a dirac mass when L(a) increases. In practice, the discretization M must be large enough to catch the maximum of h ? Sa(hd h). As we will see in the next section this is not the case when the asymptotic wavefront has curvature or just adding noise. Then, the maximum of the peaked function is moved and a confidence interval can be estimated thanks to the stability result [4]. 3.7. Second order curvature correction Fig. 8 illustrates, in the case of a source point solution, the focusing around the ray direction produced by increasing a (and thus L(a), see (31)). However, the focusing property fails when linearization error is too large as on the last plot.1 It grows with the observation circle radius (also depending on a, see (27)). Increasing the number of samples M does not help. The linearization error main contribution comes from the second order term in the Taylor expansion (21). We did not manage to include this term in the NMLA filter nor to use the transport Eq. (14). Instead, we first assume there is only one ray (N = 1). This is not a very restrictive assumption as both time and direction windows can be used to select a particular ray (see Appendix C). We then make the simplest non linear approximation on the HF asymptotic solution:
we look for a ð1Þ sj locally constant curvature as if the wave-field solution was the fundamental solution Gðx; xs Þ ¼ 4i H0 x jxx . This virtual cðx0 Þ ~ where 1 and d ~ are respectively the sought after curvature of the front and direction source point xs is such that x0 xs ¼ dd d of the ray locally near the observation point, see Fig. 3. ~ plus x x a vector The argument x xs is decomposed into the sum of the main propagation direction ðx x0 ¼ ddÞ 0 s turning around x0: x0 xs ¼ r~s. Fig. 4 presents this decomposition. Thus, we use
U a ð~sÞ ’
sffiffiffiffiffiffiffiffiffiffiffi 8pd i ð1Þ x ~ þ r~sj : A0 ðx0 Þeixðuðx0 ÞdÞ H0 jdd cðx0 Þ 4 cðx0 Þ
ð33Þ
The high frequency amplitude and phase above are corrected as the global wave-field is not necessarily globally a point source solution in an homogeneous medium. We use the known HF asymptotic of the Hankel function, see (10). 1
See also Section 3.5.2 [5].
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4651
~ þ r~s. Fig. 4. Observable geometry data in the direction ~s on the circle is written as dd
Inspired by techniques from the Fast Multipole Method applied to electromagnetics [10,11], one can then consider r as a perturbation in the argument of the Hankel function and expand it. Using technical asymptotic expansions, we show in [4] ^ Fourier modes of the NMLA filter (30) is given by that a second order accurate approximation of the b
^‘ ’ Aeixuðx0 Þ ei b
ð‘2 14Þ
‘hd~ þ
2c
ð34Þ
;
xd is the large parameter and h is the angle associated to d. ~ This is to be compared to the exact solution for one where c ¼ cðx ~ d 0Þ plane wave linear data (32). We see that we get a c (curvature dependent) second order correction. From there we simply remark that
b^‘ 1 log i b^0
!
’ ‘hd~ þ
‘2 14 2c
ð35Þ
is a parabolic function of the mode number ‘ and the coefficients hd~ and 21c can be fitted by a least square approximation, see Fig. 5. We summarize the 2nd order correction of the NMLA algorithm: 1. Get an estimate hmI of hd~ over the set of all discretized directions (24) using (30) output:
mI
such that jBmI j ¼ maxjBm j: m
ð36Þ
Fig. 5. Stars are the (35) values and the solid line the fitted parabola as a function of ‘ (x-axis). The data correspond to the test case presented in Section 4.3.
4652
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 6. Travel-time for different windows, x-axis is the center of the picked time window, y-axis gives the travel-time computed using (40).
2. A least square parabolic fitting of
^‘ b 1 ‘ ! log ei‘hmI ^ i b0
! ’ a0 þ a1 ‘ þ a2 ‘2
ð37Þ
yields (a0, a1, a2) the constant, linear and quadratic coefficients of the parabola. Based on (35) the linear term coefficient is the angle correction: a1 ¼ dh ¼ hmI hd (notice the hmI phase shift in (37)) and the quadratic term coefficient depends on the mean curvature of the front: a2 ¼ 21c. 3. Correct NMLA amplitudes
0
1 ð‘2 14Þ 0 1 @ ^ i ‘dh 2c A: b‘ e Bm :¼ F
ð38Þ
Fig. 9 illustrates the improved focusing and increase amplitude at the maximum.
3.8. Travel-time estimate We mentioned in Section 2.4 the possibility of preselecting a time window around a hand picked time t0. The travel-time of the ray corresponding to the wave recorded in this time window is necessarily close to t0. The wavelet is however not a perfect wavefront and may spread over a given period of time (around 1 s in our numerical illustration Fig. 17).
Fig. 7. Seismogram with time (y-axis) and position (x-axis) picking.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4653
NMLA combined with the use of time data can provide an accurate travel-time estimation. The idea is to use the implicit redundancy of the HF travel-time over the frequency range. More precisely we have, again assuming perfect asymptotic one ray ansatz data,
Fig. 8. One source point: jBmj (y-axis) versus hm (x-axis). Varying a from top to bottom: 10, 20, 40. Vertical line is the exact ray angle. Obtained from 20% noisy data.
4654
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
BmI ðxÞ ’ A0 expix/ðx0 Þ :
ð39Þ
The x dependence is explicit in the notation since we are going to use the variation with respect to frequency,
/ðx0 Þ ’
@ x BmI : iBmI
ð40Þ
The x derivative is obtained by applying both procedures, NMLA and second order curvature correction. to the x derivative of the data. The derivative data can be computed using finite difference or using the Fourier transform of t u(t, .). Using the case study of Section 4.3, we can show that this method is more accurate than ‘‘hand picking’’ t0: we recall that we use time data only in a time window centered around this ‘‘hand picked’’ t0. The width of the window must be larger than the width of the source wavelet. In Fig. 6, we plot the travel-time estimated by (40) against t0. As long as the time window contains significant part of the signal, the travel-time is constant and correct and therefore does not depend on the picking time. 4. Numerical illustration: source point discovery 4.1. Generation of synthetic data We generate synthetic data using the following procedure: A standard FDTD scheme, second order in time and space, is used to simulate the waves generated by one or more source points. Discretization in space is chosen to maintain a minimum of 20 points per wavelength for the cutting frequency Fc = 3 xs ⁄(space step is h ¼ mincðxÞ ), in our test case Fc = 90 Hz. Time discretization satisfy the usual CFL condition (time step Fc is Dt ¼ pffiffi h ). 2maxcðxÞ
A ‘‘Ricker’’ (Gaussian derivative) is the source wavelet. The source pulsation xs (see Section 2.3) we use is 60p rad/s. We use a perfect reflecting boundary condition @@m u ¼ 0 at the bottom of the domain and Higher Order Absorbing Conditions [13] elsewhere. We record over time a seismogram (just u) at all discretization points on a line. Once the observation time t0 and surface position x0 have been picked (see Fig. 17), a smooth window in space and time is used to single out one event. We Fourier transform all the time lines, then choose a frequency x such that the frequency content of the source is significant. Typically close to xs. The NMLA frequency data is then computed on the observation circle using the one way method described in Appendix A. 4.2. Homogeneous medium 4.2.1. one source point Fig. 7 depicts the seismogram recorded over time (vertical axis). Each vertical line is the signal t ? u(t, x) recorded at a grid point x on the surface (horizontal axis). The horizontal deviation of these lines is set proportional to u. It generates oscilla-
Fig. 9. One source point: jBmj (y-axis) versus hm (x-axis). Solid line: NMLA - Solid line with + markers: NMLA + curvature correction - Vertical line: exact ray direction.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4655
tions and a shade at times when waves (also called events) reach the surface. In this particular simulation there is only one source point but two successive events due to the reflecting boundary condition at the bottom. We manually select an observation point on the surface x0 (vertical line). We also manually pick a time t0 (horizontal line) which approximately gives the travel-time of the event. Using a time window centered at t0, we can eliminate other events and retain only one significant wavefront and ray in the solution. Travel-time estimation is also discussed in Section 3.8. Fig. 8 shows the impact of a, the radius of the observation circle counted in 2p wavelength, on the basic NMLA filter (30). If a is too small, we do not have enough Fourier modes to focus, note the angle is still correctly estimated for one ray, however for more than one ray (as in the next section) a = 10 would not allow to observe cleanly separated maxima when the ray directions are close. If a is too large, then the curvature error induced by the plane wave approximation pollutes the result. Fig. 9 shows the output of the NMLA procedure (30) and also the second order correction of Section 3.7. After correction the amplitudes are better focused and the maximum more accurate. Table 1 quantifies the error on the angle and the source point discovery. Local Radon transform (LSS) completely fails while NMLA does a better job than plane wave destructor (PWD) (references on these two methods in Section 4.3). Table 1 HF ray direction analysis and source point discovery. Method
Ray angle error
NMLA basic NMLA 2nd order LSS PWD
0.5988° 0.0184° failed 0.1917°
Error (distance to source point) 0.0011 0.0071
Fig. 10. One source point (0.1, 0.2): Angle error (solid line) and Travel-time error (dashed line) in Log10 scale versus the position of the observation point x0 on the ‘‘surface’’ xo2 ¼ 1.
Fig. 11. Seismogram with time (y-axis) and position (x-axis) picking.
4656
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 12. Two source points: jBmj (y-axis) versus hm (x-axis). Dashed line is NMLA output for 2 source points. Solid line is NMLA + Gaussian window to enable separation a = 2.5 (more details in Appendix C).
Fig. 13. jBmj (y-axis) versus hm (x-axis). White noise: 20% solid line - 40% dashed line.
Fig. 14. jBmj (y-axis) versus hm (x-axis). Correlated noise: 20% solid line - 40% dashed line.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 15. Heterogeneous medium: velocity distribution.
Fig. 16. Synthetic data numerical simulation, snapshots at 4 successive times.
4657
4658
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
Fig. 17. Seismogram with time (y-axis) and position (x-axis) picking.
Fig. 10 shows the error on angle and travel-time when applying NMLA second order to several position on the surface y = 1. 4.2.2. Two source points Basic NMLA is able to recover well separated multiple rays. We treat here the case of two but it also works fine for more, see the original paper [5,15]. Fig. 11 shows the recorded seismogram for two source points. The picked observation point is at equal travel-time from the source points. Thus, the time window will contain two wavefronts and two rays. Fig. 12 shows the NMLA basic filter result on the left and the application of the Gaussian window convolution (47) (see Appendix C). The filter removes the possible interferences caused by the trailing oscillations the ray contributions to B can be considered separately. One can go back ^ of each ray contribution and perform the curvature correction and traveltime estimation (Sections 3.7 to the fourier modes b and 3.8). Table 2 HF ray direction analysis and source point discovery. Method
Ray angle
Error (distance to source point)
NMLA basic NMLA 2nd order LSS PWD
90.2970° 90.5503° 91.002° 90.7418°
0.0069 0.0878 0.0201
Fig. 18. Ray backward (downward) propagation with initial conditions specified in Table 2. The circle is centered at the true source position. Source discovered with, from left to right, NMLA/PWD/LSS ray angle at observation.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4659
4.2.3. Homogeneous medium noisy data - stability In this section we present numerical results (Figs. 13 and 14) which illustrate the stability claim on the algorithm (Section 3.5) in the two source points case. The NMLA filter (30) easily filters out white noise (see Fig. 13) as it is precisely designed to single out coherent signals. In order to test correlated noise, we perturb the data with the solution collected on a circle around the observation point but with a different radius (see Fig. 14). 4.3. Heterogeneous medium We now consider the smooth speed velocity distribution c(x) of Fig. 15. The + sign is the source point location and our observation point is the cross at the surface. Fig. 16 shows several snapshots of the FDTD simulation used to generate the synthetic seismogram in Fig. 17. Table 2 displays the ray direction found with several methods. NMLA basic and the second order correction. The Plane Wave Destructor Method PWD [21] and the Local Slant Stack LSS (Radon transform) [8,19] are popular methods in Geophysics and radar processing but are sensitive to the necessary windowing. Note that PWD only identified plane waves by construction; the Radon transform can be generalized to identify parabolas [6] but we did not tested it. Note that both PWD and LSS can be considered as image processing methods and do not use the wave models. Fig. 18 illustrates the source discovery through the proposed HF TR method (Section 2.4). NMLA second order yields a u(x0) = 0.4602 s. travel-time. We use this information together with the three different ray directions found in Table 2 to shoot the rays backwards. A circle centered around the true source position highlights the better performance of NMLA. 5. Conclusion The second order version of NMLA method is a stable local HF wave analysis tool. It is based on the true wave HF models and robust even to correlated noise. The companion paper [4] also tries to mathematically establish how close in the near field this method is going to work. It is completely automatic and has no other tuning parameters than a, the radius of the observation circle measured in 2p wavelength. A mathematical study of the error depending on a and the discretization parameter M is done in [4]. It also proposes an optimal law for a depending on the front curvature. The source deconvolution is only needed for the travel-time estimation. The curvature correction requires that only one ray is present in the data. We briefly explain in Appendix C how a gaussian filter can be use to separate multiple ray components with sufficiently distinct directions. This part should be more thoroughly investigated in a forthcoming paper. A numerical example close to a caustic was presented in [5]. Even though rays directions were recovered, this topic certainly deserves a more serious investigation. We think that NMLA can be extended to 3D and other linear models including elastic waves (see [15]). Appendix A. Seismogram data/ FT/ one-way extrapolation of linear array data We assume for simplicity that:
We are in 2D and x = (x1, x2). ^ 0 ðx1 ; xÞ is given on the surface x2 ¼ xo2 . The data u Our observation point is on the data surface x0 ¼ xo1 ; xo2 . The medium speed is constant c(x) = v in a layer containing the observation circle of radius r (see (26)).
We recall here the one way extrapolation procedure. More details, references and also the more complicated case of one way extrapolation in heterogeneous medium, can be found in [14]. First, we Fourier transform the data and Eq. (13) with respect to (x1, x2) and denote the Fourier space variable (k1, k2). Fourier transform rules that the partial derivatives @ xj corresponds in Fourier space to a multiplication by, respectively, i kj, that is
x2 ^ 2 2 ^ 0 ðk1 ; xÞ given on x2 ¼ xo2 : ðk1 ; k2 ; xÞ ¼ 0 u k1 þ k2 2 u
v
ð41Þ
^ notation and drop the x dependence for simplicity. If we restrict ourselves to real ‘‘up-going’’ k2 > 0 Fourier We keep the u mode, this is the one way assumption, and after factorization, Eq. (41) can be simplified into
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!
k2
x2 c2
2 ^ ðk1 ; k2 Þ ¼ 0: k1 u
ð42Þ
The significant real Fourier modes of the solutions have to satisfy the (so called) ‘‘dispersion’’ relation characterizing the paraxial operator:
4660
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x2
k2 ¼
c2
2
k1 :
ð43Þ
The solution can be computed as a function of (k1, x2). Indeed after an inverse Fourier transform in x2, we get the ODE (we ^ notation again) abusively keep the u
^ k1 ; xo2 ¼ u ^ ðk1 ; x2 Þ ¼ k2 u ^ ðk1 ; x2 Þ u ^ 0 ðk1 Þ; i@ x2 u
ð44Þ
where k2 is given by (43), with trivial solution The exact solution of (44) at depth x2 is
^ 0 ðx; k1 Þeik2 x2 : ^ ðk1 ; x2 Þ ¼ u u The solution at any point (x1, x2) is recovered by inverse Fourier transform from k1 to x1. This procedure generates the data on a circular array for up-going rays in constant medium at the price of two 1D Fourier transforms. It could also be used to filter out down-going rays. As we restrict to observation point near the x2 ¼ xo2 depth of the linear array, it is possible to use a x1 windows which size depends on the finite speed v and the resulting dependence cone. Appendix B. When the wavelet W xs is unknown In that case, the source deconvolution (20) is not possible so we just remove the half primitive factor
^ ðx; xÞ u
pffiffiffiffiffiffiffiffiffiffi ^ðx; xÞ: ixu
ð45Þ
This implies that the NMLA algorithm expected output is multiplied by the frequency domain wavelet
B
c x ðxÞ: BW s
ð46Þ
The curvature correction method of Section 3.8 is unchanged since the amplitudes are removed in (37). The travel-time estimation method (40) does not work since the amplitude BmI contains now the wavelet with an unknown x dependent phase. Appendix C. Time and ray direction windows As shown in Fig. 17 we assume that, after picking the observation point x0 ¼ xo1 ; xo2 on the surface data (the vertical red2 line), we select an observation time in the seismogram (the horizontal red line). This can be done by the user or an ad hoc automatic signal processing procedure. The advantage of time picking, is that rays crossing at the observation point with different travel-time may be separated using a time window. In Fig. 17 we see that a second later wavefront passes the observation point. If we just Fourier transform in time and then perform NMLA, we expect the (Bm)’s to characterize a two ray ansatz. Using a smooth time window, we are able to isolate for instance only the first arrival. We then only have one ray and can perform the curvature correction as explained in Section 3.7. However multiple rays with same or close travel-times may cross at the observation point. Then, the time window is unable reduce the data to only one ray (N = 1). The NMLA filter output is the sum of cardinal sine Sa functions centered at the different ray directions. Even if the rays directions are well separated interferences make it impossible to locate ray directions as the maximum of the two components. A way to fix this problem is to remove the trailing oscillations of Sa by a convolution of B with a smooth Gaussian weight function for instance. The convolution does not affect the Sa components maximum position. The convolution is simply done using an additional multiplicative filter in Fourier space before the last Fourier transform in (30):
^‘ b
^‘ w ^‘ b
^ ‘g ¼ F fw
n o
a 2 exp2hm :
ð47Þ
The parameter a can be used to tune the Gaussian filter and determines the resolution of the filter or the minimum width between ray direction in order to separate them. A large a yields a small window but a must be small enough to remove the high frequency oscillations. We heuristically found that a = 2.5 seems optimal for our test cases. Once ray components have been separated one can go back to the Fourier modes and perform the curvature correction (Section 3.7). Fig. 12 shows the filter effect comparison for two source points in homogeneous medium. References [1] T. Abboud, J.-C. Nédélec, B. Zhou, Méthode des équations intégrales pour les hautes fréquences, C.R. Acad. Sci. Paris Sér. I Math. 318 (2) (1994) 165–170. [2] M. Balabane, V. Tirel, Décomposition de domaine pour un calcul hybride de l’équation de Helmholtz, C.R. Acad. Sci. Paris Sér. I Math. 324 (3) (1997) 281–286. 2
For interpretation of colour in Fig. 17, the reader is referred to the web version of this article.
J.-D. Benamou et al. / Journal of Computational Physics 231 (2012) 4643–4661
4661
[3] Jean-David Benamou, An introduction to Eulerian geometrical optics (1992–2002), J. Sci. Comput. 19 (1–3) (2003) 63–93 (Special issue in honor of the sixtieth birthday of Stanley Osher). [4] Jean-David Benamou, Francis Collino, and Simon Marmorat. Second order stable numerical microlocal analysis of harmonic plane wave and source point wavefields, submitted for publication.