Acoustic full waveform inversion for 2-D ambient noise source imaging

SUMMARY We present a method for estimating seismic ambient noise sources by acoustic full waveform inversion of interstation cross-correlations. The method is valid at local scales for laterally heterogeneous media, and ambient noise sources conﬁned to the Earth’s surface. Synthetic tests performed using an actual ﬁeld array geometry, are used to illustrate three unique aspects of our work. First: the method is able to recover noise sources of arbitrary spatial distribution, both within and outside the receiver array, with high ﬁdelity. This holds true for complex velocity models and does not require a good initial guess for inversion, thereby addressing an outstanding issue in the existing research literature. Second: we analyse the extent of biases in source inversion that arise due to inaccurate velocity models. Our ﬁndings indicate that source inversion using simpliﬁed (e.g. homogeneous) velocity models may work reliably when lateral variations in velocity structure are limited to 5 or 10% in magnitude, but is vitiated by strong variations of 20% or higher, wherein the effect of strong scattering and/or phase distortions become signiﬁcant. Finally, our technique is implemented without the adjoint method, which is


INTRODUCTION
The study of the ambient seismic field, or seismic ambient noise as it is popularly known, is now firmly entrenched in mainstream seismological research.Ambient seismic sources can shed light on such natural phenomena as ocean wave coupling with the seafloor (e.g.Juretzek & Hadziioannou 2016), sediment transport in rivers (Tsai et al. 2012), glacier hydrology and dynamics (e.g.Aso et al. 2017;Labedz et al. 2022), tropical cyclones (e.g.Retailleau & Gualtieri 2019) and underground hydrothermal activity (e.g.Cros et al. 2011).On the other hand, seismologists are widely interested in using ambient noise as a tool for studying Earth structure, typically by applying interferometric techniques to extract meaningful signals from noise recordings (e.g.Shapiro & Campillo 2004).Even in this case, source information is essential because the spatial distribution of noise sources effectively determines whether inter-station Green's functions can be accurately recovered from noise cross-correlations (Roux et al. 2005;Snieder 2004).There is ample evidence for biases in structure estimation, arising from realistic distributions of ambient noise sources on Earth (e.g.Kimman & Trampert 2010;Yao & van der Hilst 2009).If one chooses instead to leverage the power of full waveform inversion (FWI), the assumption of Green's function retrieval is dropped and both sources and structure must be simultaneously estimated (Sager et al. 2018b;Zhou et al. 2022).Thus, regardless of one's particular interest in seismic ambient noise, the ability to determine the strength and locations of the noise sources is of vital importance.
Traditional, computationally cheap methods for locating ambient sources include beamforming (e.g.Gal et al. 2015;Gerstoft & Tanimoto 2007), matched-field processing (e.g.Cros et al. 2011) and backprojection (e.g.Liu et al. 2016) as well as cross-correlation based imaging (e.g.Ermert et al. 2016;Tian & Ritzwoller 2015).The last few years have witnessed the emergence of FWI methods, which seek to match observed noise cross-correlations with theoretically modelled Acoustic FWI of ambient noise sources 3 ones, incorporating their finite frequency sensitivity to spatially distributed sources (Fichtner et al. 2017;Tromp et al. 2010).This paper focusses on FWI to recover noise source distribution with the assumption of a known structure model (Datta et al. 2019;Ermert et al. 2017Ermert et al. , 2021;;Igel et al. 2021;Xu et al. 2019Xu et al. , 2020)).These 'source inversion' techniques are useful in their own right, and essential components in the toolkit for the larger problem of full waveform noise cross-correlation tomography (e.g.Sager et al. 2020).FWI, by definition, entails wave-equation based forward modelling.However, different approximations to the seismic wave equation (e.g.acoustic vs. elastic) or different assumptions about the medium of propagation, lead to a family of methods achieving varying degrees of modelling rigour.From analytical modelling in homogeneous or laterally homogeneous media (Datta et al. 2019;Xu et al. 2019Xu et al. , 2020) ) to numerical simulations in spherically symmetric (Igel et al. 2021) or full 3-D Earth models (Ermert et al. 2017(Ermert et al. , 2021) ) at global scales, a variety of methods have been proposed.In this study, we present an approach using acoustic modelling in 2-D media that incorporates laterally heterogeneous structure information.Xu et al. (2020) used classic surface wave analysis on ambient noise data (Bensen et al. 2007) to estimate a uniform phase velocity required for forward modelling.With our approach, one can go a step further by using the 2-D phase or group velocity maps obtained from ambient noise surface wave tomography (Shapiro 2019).Independent of our choice of modelling scheme, another highlight of this study, achieved by building on the work of Datta et al. (2019), is the ability to localize sources outside the sensor array.Previous non-global studies, deploying a range of strategies for data measurement and misfit definition, have had limited success in this regard.For example, Xu et al. (2019) reported that both traveltime and waveform inversion are only able to estimate rough source directions, rather than actual locations and shapes, when sources are outside the array.Ermert et al. (2020) demonstrated recovery of outside-array sources for a 1-D medium, but for a waveform-difference misfit function, which admittedly suffers from high sensitivity to velocity model inaccuracies.Their 'asymmetry misfit', discussed below, is robust with respect to velocity inaccuracies but offers poor resolution outside the sensor array.Against the backdrop of these challenges, our method recovers external source shapes and sizes fairly accurately, provided a moderate inter-sensor path density is avail-able.It does not require a good initial source model, and can tolerate 2-D velocity model errors of up to 10%.
The two aforementioned features of our method -incorporation of 2-D structure, integrity of resolution outside the station array -lend themselves readily to a scrutiny of the source-structure trade-off in seismic interferometry (Fichtner 2014), which is not very well documented.Xu et al. (2019) reported the failure of source inversion in case of inaccurate structure models, but only for homogeneous halfspace models devoid of lateral variations.Other studies have argued that misfit functions defined using measurements of waveform asymmetry are very weakly sensitive to lateral structural variations (Igel et al. 2021;Sager et al. 2018a).These studies quantify cross-correlation waveform asymmetry via a logarithmic ratio of energies in the the causal and acausal branches.In our approach, we measure waveform energy on a single cross-correlation branch (rather than the ratio of energies on the two branches), and find that structure (velocity) information does impact source inversion.
In the following, we first present the methodology for forward and inverse modelling (Section 2), followed by synthetic tests (Section 3) which form the basis of our focussed Conclusions (Section 4) and a more general Discussion (Section 5).Datta et al. (2019) introduced a method for estimating noise source directionality in a homogeneous medium.The limitation of estimating only directions was imposed by the choice of model parameterization, and the homogeneous medium assumption allowed for cross-correlations to be modelled analytically.In this study, we extend their method on both fronts.A similar but less restrictive model parameterization allows actual source locations and shapes to be estimated, and integration with a numerical solver eliminates the need to assume a homogeneous medium.However the choice of measurement, calculation of kernels and inversion strategy, remain essentially unchanged.

CROSS-CORRELATION MODELLING AND INVERSION
Acoustic FWI of ambient noise sources 5

Forward modelling
The foundation of the method is a forward model for computing the cross-correlation C(x α , x β ) between any pair of receivers α and β: where ω is angular frequency, P is the source power spectrum, G is the (scalar) Green function for the medium, σ is the source strength, and the position coordinate x is limited to a horizontal plane (approximating an area on the Earth's surface).The asterisk denotes complex conjugation.In case of elastic modelling, appropriate components of the elastodynamic Green tensor need to be used instead of the scalar G. Equation 1 follows from the formulation of Tromp et al. (2010), subject to the following assumptions (e.g.Malkoti et al. 2021;Xu et al. 2019): (i) all ambient seismic sources are confined to the Earth's surface, so that the integral in ( 1) is a surface integral (ii) all sources have the same spectral shape P (ω), so that the power spectral density (S) of noise sources can be separated into its positional and frequency contributions, i.e.S(x, ω) = P (ω)σ(x).
The frequency-domain expression (1) can be evaluated semi-analytically when closed form solutions for G are available -such as in a homogeneous (Datta et al. 2019) or 1-D (Malkoti et al.

2021) medium.
In anticipation of numerically solving the wave equation for heterogeneous media, we recast the expression for crosscorrelation in equation ( 1) in terms of the wavefields u: with u(x, x α ; ω) = P (ω) 1/2 G(x, x α ; ω).Note that source-receiver reciprocity has also been invoked for computational efficiency, turning every receiver location into a source (e.g.Hanasoge 2014; Xu et al. 2019).The monochromatic wavefields in equation ( 2) can be computed in the frequency domain (e.g.Kumar et al. 2022).However, we use a time-domain finite difference solver implemented using "Devito" (Louboutin et al. 2019;Luporini et al. 2020), a Python package that can generate optimized codes based on symbolic input of governing equations using SymPy (Meurer et al. 2017).We use a uniform finite-difference stencil in Devito that is second-order in time and fourth-order in space.Equation (2) in the time domain is We obtain u(x, x α ; τ ) as the solution to the acoustic wave equation with a source wavelet equal to the inverse Fourier transform of P (ω) 1/2 .Provided that the computed wavefields u(x, x α ) can be stored in memory, N numerical simulations are required to obtain all possible cross-correlations for an N -receiver array.In the inverse problem of estimating σ(x), these N simulations, once performed at the start, are not required to be repeated every iteration, because the velocity model is held fixed.

Model parameterization
As in Datta et al. (2019), the model space is parameterized using a set of 2-D Gaussian basis functions, B j (x).In this study, we introduce a non-negative parameterization to ensure that the iterative optimization scheme (Section 2.3) does not lead to unphysical, negative values for σ(x): m j being the basis function coefficients that are directly inverted for.We note that alternate means of enforcing positivity are available, such as in Xu et al. (2019).A second feature of our parameterization is that the basis functions are present uniformly throughout the domain, rather than merely in a ring surrounding the receiver array (Datta et al. 2019).The size and spacing of the Gaussian basis are user-controlled parameters, held fixed during inversion.The default specifications used in this study (Section 3 and Table 3.1), lead to a total of 625 basis functions, as shown in Figure 1.

Inversion Strategy
We follow the inversion strategy detailed in Datta et al. (2019, Appendix A), subject to minor changes necessitated by the model parameterization, equation (4).To summarize, we use a misfit Acoustic FWI of ambient noise sources 7 function defined by Hanasoge ( 2013): where the index i runs over receiver pair combinations, and E, the measurement, is the waveform energy in a time window of interest w(t) on either the positive or negative cross-correlation branch: When working with real data, one can define a measurement window aimed at encompassing the dominant signal arrivals on each correlation branch (Datta et al. 2019).However in this synthetic study, windowing is not necessary and we let w(t) span the entire positive or negative correlation branch.
Local optimization techniques require the gradient of χ, or a misfit kernel K satisfying δχ = − K(x)δσ(x)d 2 x.This misfit kernel is the sum of individual source kernels for each receiver pair (K i ), weighted by the corresponding misfit (Hanasoge 2013): Using (4), we get the following expression for the gradient (g) of χ: K, and therefore g, can be obtained by either evaluating (7), which requires all the individual source kernels K i , or by adjoint techniques, which build K from adjoint wavefields and 'event kernels' (Tromp et al. 2010), without access to the individual K i .We use the former, non-adjoint approach.Using an analogy from earthquake traveltime tomography, this is akin to the difference between calculating individual banana doughnut kernels (e.g.Dahlen et al. 2000) or not (Tromp et al. 2005).
Our approach remains computationally feasible because numerical simulations are not required to obtain the K i in every iteration, as explained in Section 2.1.Furthermore, we parallelise the implementation over the individual receivers.The advantage of this inversion strategy is that in addition to the gradient, it gives us access to the Jacobian through the individual kernels.
Combining ( 7) and ( 8), the elements of our Jacobian matrix (J) are: The Jacobian can in turn be used to build an approximation to the Hessian operator, paving the way for optimization methods such as the Gauss-Newton method, which we use.For details of how this is done, and complete derivations of equations ( 8) and ( 9), the reader is referred to Appendix A of Datta et al. (2019).Here we wish to highlight two subtleties: (i) Our least-squares inversion is regularized or damped in the usual way, by the L 2 norm of model misfit, weighted by a damping parameter determined by L-curve analysis (Hansen & O'Leary 1993).This corresponds to a diagonal model covariance matrix.Additional smoothing, which we do not use in this synthetic study, can be imposed through a banded model covariance matrix.
(ii) The entire workflow described in this section, from measurement to optimization, is executed separately for positive and negative cross-correlation branches.At the end of each iteration, the model update from the two sets of optimizations, is averaged to obtain the updated model.

SYNTHETIC TESTS
We conduct a series of synthetic tests using the field array geometry of Datta et al. (2019), which corresponded to an exploration seismic deployment by Shell.This choice is motivated by the Acoustic FWI of ambient noise sources 9 relatively large number of receivers (289) in the array, which allows us to explore the effects of different array sizes and path densities, by randomly picking different subsets of receivers.We present results in this paper for two sub-array geometries, one with 20 and the other with 50 receivers (Figure 2, parts b and c, respectively).Synthetic tests with more than 50 receivers did not produce significantly different results.
Our default modelling domain size is 50 km × 50 km, and simulations are performed with P (ω) set to a Gaussian centred at 0.2 Hz, with a standard deviation of 0.05 Hz.Assuming a wavespeed of 2 km/s to define a homogeneous medium, we have a dominant wavelength (λ 0 ) of 10 km.The shortest wavelength is ≈ 5 km when velocity perturbations are introduced (Section 3.1), so we choose a conservative grid spacing of 0.5 km, which ensures a minimum of 10 grid points per wavelength in all cases.The size of the model basis functions (5 km across) is comparable to the shortest wavelength.Table 3.1 summarizes the various parameters used in our simulations.
We present results for two test source models, TSM-1 (Fig. 2d) and TSM-2 (Fig. 2g).In actual field scenarios, TSM-1 may represent, for example, roads situated close to or passing through, the receiver array.TSM-2 is more of a toy model, containing sources of varying strengths located entirely outside the array.The test models are built independently, without using the parameterization (4), which is reserved for the inverse problem.
For both test models, "observed data" are synthetically generated by forward modelling and adding noise to the resulting cross-correlation waveforms.We add Gaussian noise with a standard deviation equal to 10% of the average RMS value of all the computed crosscorrelations (e.g.Ermert et al. 2020).The frequency band of the noise (0.05-1 Hz) is kept wider than that of the power spectrum P (ω) used in the simulations (∼ 0.05 − 0.35 Hz).
For inversion, a quasi-uniform source distribution (Figure 1 and 2a) is used as a prior as well as the initial model.This initial model has very low amplitude compared to the test models, but it cannot be zero because of the logarithmic misfit function (5) used in the inversion.Finally, we note that inversions may be performed with either analytical or numerical modelling, because a homogeneous medium is assumed.For consistency with Section 3.1, we have presented results A. Datta, B. Shekar, P.L. Kumar obtained with numerical modelling.The results with analytical Green functions were similar to those in Figure 2 and have not been shown.
Figure 2 shows that the inversion recovers true source locations and shapes fairly well, for both test models.As expected, results with 50 receivers are better than those with 20, but the improvement is modest, especially for sources lying outside the receiver array (TSM-2).The pixellated appearance of the inversion results as compared to the test models, is due to the parameterization (4) and the size of basis functions used (Table 3.1).In this way, our model parameterization dictates the resolution of the inverse problem, for a given frequency band and array geometry.To put this limitation in perspective, it is important to remember that our method, in its current implementation, is unsuitable for distances over which Earth's sphericity becomes important.
In this paper we choose to focus on the question of how source inversions are impacted by unmodelled lateral heterogeneities in velocity structure.For the sake of uniformity, we adhere to the simulation parameters of Table 3.1, so that inversion results can be benchmarked against those of Figure 2.  to the test models in order to illustrate its low amplitude.

Impact of heterogeneous velocity structure
To investigate the effects of heterogeneous structure, we use two types of velocity models: (i) a single low velocity anomaly (Figure 3) -this is a smooth model containing 'long wave-length' velocity variations (length scale λ 0 ), but with an average velocity lower than 2 km/s.It represents a baseline shift relative to the homogeneous velocity model used in Fig. 2).
(ii) a checkerboard (Figure 4) -here the average velocity is equal to that of the homogeneous model (no baseline shift), but there are 'short-wavelength' velocity variations (length scale ∼ λ 0 ).
With both velocity models, we consider perturbation amplitudes of 5, 10 and 20% relative to the reference homogeneous velocity of 2 km/s.Test data are generated using these heterogeneous velocity models.Inversions are then performed both with and without knowledge of the accurate velocity models.In this section, we present tests performed with the 50-receiver array geometry only.
When the velocity information provided in inversion is accurate (Fig. 3(d-f)), we obtain results that are all nearly identical to the homogeneous velocity case.The inverse modelling can tolerate errors in velocity information up to 10%, but artefacts appear for the case of 20% velocity heterogeneity (Fig. 3(g-i)).Also, sources outside the array are not retrieved accurately in this case.
With the checkerboard velocity model, inversion artefacts, at 20% heterogeneity, are more pronounced (Figure 4), reflecting the complications of wave propagation through short-wavelength structural variations.In Figure 5, we compare a few waveform fits from the long wavelength anomaly model (Fig. 3i) with those from the checkerboard model (Fig. 4i) at 20% velocity perturbation.The long-wavelength anomaly primarily affects the phase of the propagating waves.For receiver pairs separated by large distances (e.g.Fig. 5c), the observed and predicted waveforms are separated by more than half-a-cycle (∼ 2.5 s), leading to poor inversion results: this is the well known phenomenon of cycle skipping (Virieux & Operto 2009).We note that the energy-based misfit function is not immune to cycle skipping, but analysis of misfit functions in this context is beyond the scope of this study.The short-wavelength checkerboard model, on the other hand, produces significant multiple scattering.Scattering is visible in the cross-correlation waveforms (e.g.inversion is performed using a homogeneous velocity model.In all cases, the starting (source) model for inversion is the same as in Figure 2a.
regions are reasonably well recovered, for velocity inaccuracies of up to 10%.For the case of unmodelled 20% velocity heterogeneity, the low velocity anomaly produces some distortion of sources, while strong aretefacts appear with the checkerboard velocity model.We also performed tests with a long wavelength checkerboard model and verified that the phase perturbations due to positive and negative anomalies compensate each other, leading to improved inversion results even at 20% velocity perturbation.
(a) It should be emphasised that the σ(x) starting model is the same (see Fig. 2a) for all tests in this section.On the other hand, the optimal value of the damping parameter used in inversion, is determined separately in each case by independent L-curve analysis.In this paper, we have considered two contentious issues in the ambient seismic source inversion literature -source resolution outside the sensor array (empirical constraints only), and sensitivity of source inversion to velocity structure heterogeneity.Based on our synthetic tests, we conclude that: (i) it is possible to accurately recover source distributions up to a certain distance from the sensor array, as long as the forward modelling theory (Section 2.1) is valid.
(ii) accurate modelling of velocity structure, within the acoustic approximation, has limited impact on ambient noise source inversion.Detrimental effects of assuming a simplistic velocity model are seen only when the assumed model is a very poor approximation to the true structure.
For media with weak lateral heterogeneity, even the simplest approximation of a homogeneous medium (with a reference wavespeed accurate to within 10%) may be acceptable.It becomes unacceptable when true wavespeed variations exceed 20% of the assumed reference value, due to unmodelled scattering and/or phase distortions.This has a bearing on realistic scenarios where the velocity structure is not known a priori.With real data, one can expect unmodelled lateral heterogeneity to have a greater impact on source inversion, because real data would be beset with 'modelling error' (which was absent in the synthetic data of Section 3).
These conclusions are specific to the particular technique used, as elaborated below.

DISCUSSION
We have introduced an acoustic full waveform inversion technique for ambient noise source distribution.Although more sophisticated techniques, which account for Earth's three-dimensional, elastic (and anelastic) structure already exist (e.g.Ermert et al. 2017Ermert et al. , 2021)), our method occupies a niche within the field.Compared to 3-D elastic FWI, it incorporates a significant amount of modelling rigour, at a fraction of the computational cost.Our implementation does not require writing to disk during simulations, and most of our simulations were performed on a desktop computer with 64 GB RAM.A second convenience of our technique is that it can leverage results from
For our chosen simulation parameters, we have verified with 'spike tests' that it is possible, in parts of the model domain, to accurately recover a source represented by a single basis function.A formal resolution analysis (e.g.Xu & Mikesell 2022), which would yield resolution as a function of spatial position, is beyond the scope of this paper.Some examples of waveform fits involved in the inversions of Figure 2, are shown in AppendixA, where we address the question of how much noise our inverse method can tolerate.Reasonable inversion results, with only a marginal drop in quality, have been obtained with up to 50% added noise.However, for the sake of clarity in the next section, we use 10% noise amplitude as the default setting throughout the paper.Another question which arises at this point is with regard to the effect of distance -up to what distance from the array can sources be faithfully recovered by inversion?The tests of Figure2contain sources at a maximum distance of 25 km from the array centre.In Appendix B, we enlarge the domain to incorporate sources at larger distances, and find that far-away sources can be significantly distorted, as observed byErmert et al. (2020). oise

Figure 2 .
Figure 2. Source inversion results with homogeneous velocity models.(a) Starting model for all inversions -same as Figure 1 but with different normalization as described below.(b)-(c) Tomographic setup used for inversion: map with 20 and 50 receivers selected respectively (red triangles), along with inter-receiver ray paths (grey lines).(d) TSM-1, shown with the 50 receivers selected in (c).(e)-(f) Inversion results for TSM-1, using 20 and 50 receivers respectively.(g) TSM-2, shown with 50 receivers and with a few receivers labelled for association with subsequent figures.(h)-(i) Inversion results for TSM-2, using 20 and 50 receivers respectively.All source distributions are shown normalized by their maximum value (to focus on relative source strength), with the exception of the starting model (a), which is normalized with respect

Figure 3 .
Figure 3. Source inversion results with 2-D velocity models.Top row: Velocity models with a perturbation amplitude of (a) 5% (b) 10% (c) 20%.Middle row: Inversion results for TSM 1, obtained using an accurate velocity model, i.e. velocity models (a)-(c) are used to generate synthetic test data, as well as in source inversion leading to the corresponding results (d)-(f).Bottom row: Results of inversion performed without an accurate velocity model, i.e. velocity models (a)-(c) are used to generate synthetic test data, but source

Figure 4 .
Figure 4. Similar to Figure 3, but using checkerboard-style velocity models.

Figure 5 .Figure 6 .
Figure 5. Fit to crosscorrelograms for selected receiver pairs (highlighted in Figure 2g) corresponding to the single low velocity anomaly (a -c) and checkerboard velocity models (e -g), both with a velocity perturbation magnitude of 20%.Cycle skipping is readily apparent in (c).Part of the waveforms are highlighted in (g) to show the effect of scattering.(d), (h) Inversion statistics: after inversion (orange outline bars), no. of measurements in the central histogram bins (fit accuracy at least 72%) is approximately 60%for the single anomaly (d) and 55% for the checkerboard model (h).This refers to all measurements, i.e. positive and negative branches combined.The blue filled bars represent the measurements at the start of the inversion, i.e. for the initial model.For reference, 'good inversions', such as plots (d-f) in Figures3 and 4, put > 99% of measurements in the central bin of similarly constructed histograms.

Table 1 .
Simulation parameters used.Disambiguation: here σ denotes the standard deviation of 2-D Gaussian functions used as model basis.