Matched Field Processing for complex Earth structure

4 Matched Field Processing (MFP) is a technique to locate the source of a recorded wave 5 ﬁeld. It is the generalization of beamforming, allowing for curved wavefronts. In the standard 6 approach to MFP, simple analytical Green’s functions are used as synthetic wave ﬁelds that 7 the recorded wave ﬁelds are matched against. We introduce an advancement of MFP by 8 utilizing Green’s functions computed numerically for real Earth structure as synthetic wave 9 ﬁelds. This allows in principle to incorporate the full complexity of elastic wave propagation, 10 and through that provide more precise estimates of the recorded wave ﬁeld’s origin. This ap- 11 proach also further emphasizes the deep connection between MFP and the recently introduced 12 interferometry-based source localisation strategy for the ambient seismic ﬁeld. We explore this 13 connection further by demonstrating that both approaches are based on the same idea: both 14 are measuring the (mis-)match of correlation wave ﬁelds. To demonstrate the applicability 15 and potential of our approach, we present two real data examples, one for an earthquake in 16 Southern California, and one for secondary microseism activity in the Northeastern Atlantic 17 and Mediterranean Sea. Tutorial code is provided to make MFP more approachable for the 18 broader seismological community.


22
The ambient seismic field has become an attractive target of seismological studies over the last two 23 decades (Nakata et al., 2019). Interferometry of this complex wave field, combined with increased 24 non-peer reviewed preprint submitted to EarthArXiv station density, has enabled detailed studies of Earth's structure (e.g., Shapiro et al., 2005;Lu et al., 25 2018; Schippkus et al., 2018) and its temporal changes (e.g., Brenguier et al., 2008;Hadziioannou 26 et al., 2011). Such studies rely most commonly on seismic wave fields generated by the interaction 27 between the oceans and the solid Earth, so-called microseisms. Understanding the exact mechanism 28 for this interaction has been a challenge for more than half a century (Longuet-Higgings, 1950; 29 Hasselmann, 1963; Ardhuin et al., 2015) and some open questions remain, e.g., about the emergence inversion, are not applicable to ambient seismic noise due to the complexity of the analysed wave 47 forming the beam (Rost & Thomas, 2002). The quality of the beam is evaluated, giving the so-called 66 beampower. Other formulations of this method exist, e.g., an equivalent cross-correlation approach 67 (Ruigrok et al., 2017). Beamforming has been widely adopted by the seismological community and is wave assumption: sources have to be far from the array, the wave field has to be strictly coherent 75 across stations, and the array geometry limits the resolution capabilities (Rost & Thomas, 2002). 76 For a recent review of beamforming and polarization analysis see Gal & Reading (2019). 77 A new source localisation strategy based on seismic interferometry has been introduced in recent 78 years as an attractive alternative, sometimes referred to as kernel-based source inversion (Ermert 79 et al., 2016). The goal of this approach is not to determine the angle of arrival, but to directly  An inhomogeneous distribution of sources results in asymmetric cross-correlation functions, indicated 84 by the thickness of the wave fronts in Fig. 1c (Paul et al., 2005). In practise, this asymmetry is 85 usually quantified by comparing the causal and acausal part of each correlation function. In the 86 interferometry-based approach, synthetic cross-correlation functions are computed for a given source 87 distribution and compared against cross-correlation functions from real data. The mismatch between 88 non-peer reviewed preprint submitted to EarthArXiv the two is evaluated (e.g., by quantifying amplitude asymmetry), the source model perturbed, and a 89 best-fit source distribution is found via gradient descent in an iterative manner (Ermert et al., 2016). 90 Recent work has focused on improving efficiency (Igel et al., 2021b), the mismatch measure (Sager 91 et al., 2018), or expanding the method to multiple frequencies (Ermert et al., 2021). The advantages 92 of this approach are the stability of results, not as strict requirements on station geometry, and a 93 comprehensive theoretical foundation. Its disadvantages lie in computational cost, treatment of 94 recorded data, through that introduction of assumptions, and loss of temporal resolution.

95
Another approach that has gained some popularity in seismology in recent years is Matched 96 Field Processing (MFP). MFP is the generalisation of beamforming to allow arbitrary wavefronts 97 (Baggeroer et al., 1993). This approach has been developed in ocean acoustics, where coherency  In this paper, we introduce an advancement of MFP to incorporate real Earth structure and 112 account for the complexity of seismic wave propagation. In the following, we introduce the standard 113 MFP approach, demonstrate its shortcomings, and present our solution by incorporating realistic 114 Green's functions. We discuss implications of our approach, strategies to cope with some of them, 115 how different disciplines and localisation approaches intersect, and finally demonstrate the applica-  The MFP algorithm is straight-forward: For a given potential source location, a synthetic wave field 119 is computed and matched against the recorded wave field, i.e., the seismograms, taking coherency 120 of the wave fields across stations into account. This match is evaluated and compared against 121 other potential source locations. The potential source location with the highest score or beampower 122 (representing the best-matching synthetic wave field) is the resolved source location.

123
More precisely, spectra d(ω, x j ) are computed from the recorded seismograms at each receiver 124 position x j . The cross-spectral density matrix is computed as with * denoting the complex conjugate. K jk (ω) holds all information about the recorded wave 126 field and encodes its coherency across stations; it contains the cross correlations of the seismograms 127 from all station pairs. Following Bucker (1976), auto correlations are excluded, i.e., only components 128 k = j are computed and later utilized. This is particularly useful in the context of ambient seismic 129 noise, because noise wave fields are likely to be only weakly coherent across stations.

130
The synthetic wave field, i.e., the seismograms expected at each station from the candidate

135
The match of the two wave fields represented through K jk (ω) and s(ω, x j , x s ) is then estimated 136 through a so-called beamformer or processor. The most straight-forward beamformer is the con-137 ventional beamformer, which in its most compact form in vector notation is often written as (e.g., non-peer reviewed preprint submitted to EarthArXiv Other estimators of beampower exist, and their development is an active field of research (e.g., 144 Capon, 1969; Schmidt, 1986;Cox et al., 1987;Cox, 2000;Gal et al., 2014;Zhu et al., 2020).

145
Beamformers are often classified into conventional (eq. 3), adaptive (e.g., Capon, 1969;Cox et al., 146 1987; Cox, 2000) and sub-space beamformers (e.g., Schmidt, 1986). Adaptive beamformers aim to 147 increase resolution of the beampower distribution by increasing sensitivity to signal, but inherently 148 rely on high signal-to-noise ratio (SNR). The increased resolution is also accompanied by increased 149 computational cost, e.g., the Capon beamformer involves computing the inverse of K jk (ω) (Capon,150 1969). Sub-space detectors such as MUSIC (Schmidt, 1986) involve computation of the eigenvectors were able to resolve multiple sources this way. One of the expressed goals of the approach we 154 introduce in this paper is to be able to locate sources of ambient seismic noise, and as such SNR of this in the context of our approach we introduce here is beyond the scope of this paper. entirely, as is usually done, makes standard MFP equivalent to delay-and-sum beamforming without 171 the plane-wave assumption (Bucker, 1976). More on this in section 2.4.

172
In the simplest possible study target, i.e., a single stationary source in an isotropic, homogeneous 173 non-peer reviewed preprint submitted to EarthArXiv medium with constant velocity v = const and only straight-ray propagation of a single phase, the We propose to use Green's functions computed numerically for Earth structure directly as the syn-202 thetic wave field s(ω, x r , x s ) ("our approach") instead of the analytical form described above (eq. 203 4, "standard approach"). Effects such as dispersion and multiple wave types are then inherently ac-204 non-peer reviewed preprint submitted to EarthArXiv counted for, even for simple 1D media. If the Green's function are computed for a 3D Earth, further 205 effects such as focusing and defocusing, wave-type conversion, and coupling can all be accounted 206 for, increasing the precision of this approach further. 207 We demonstrate our method with synthetic examples for a broadband and a narrowband explosion 208 source (Fig. 2). The setup consists of two small arrays of three stations each that record the wave 209 field emitted by a seismic source located at the surface between them. The medium is a 3D

239
MFP originated as a narrowband localisation technique (Bucker, 1976) and has been adopted for further effects would need to be considered, as described above. One approach to this is empirical  Green's functions for complex Earth structure is expensive, which is why we rely in our analysis on  Therefore, a strategy is required to correct for amplitude terms. Numerically computed Green's 294 functions for real Earth structure inevitably contain amplitude terms. Several approaches may appear 295 reasonable to correct for them: correcting for amplitude decay (Fig. 3b), time-domain normalisation 296 (Fig. 3c), and spectral whitening (Fig. 3d). Without any treatment of amplitudes, the beampower 297 distribution is heavily biased by distance to stations (Fig. 3a). This effect is more pronounced

317
With this approach, we resolve the beampower peak close to the true source (Fig. 3c), but introduce This approach successfully retrieves the correct source location and does not appear to introduce 336 any unwanted biases (Fig. 3d).

337
From our tests, whitening the spectra of the synthetic wave fields (Fig. 3d) appears to be the 338 most advantageous approach, and follows the original formulation of MFP (Baggeroer et al., 1993).

339
It introduces no alteration of the recorded data, eliminates attenuation and spreading effects, removes

355
In ocean acoustics, other terminology is more common for some of these concepts. The dis-356 tribution of beampowers may be called ambiguity surface (Bucker, 1976), intended to express the 357 emergence of sidelobes for narrowband sources (Fig. 2c,d) e −iωt( x j , xs) are used, and only the manner in which t( x j , x s ) is estimated is adapated to use relative 366 distances perpendicular to the plane wavefront (Fig. 1b) instead of distances to potential source 367 locations (Fig. 1d). In that case, MFP is exactly delay-and-sum beamforming, regardless of whether 368 plane waves or curved wavefronts are used (Bucker, 1976 Here, d * k d j is the correlation of the recorded wave fields and s * j s k correlation of the synthetic wave fundamentally the same may not be intuitive at first, especially considering the strikingly different 394 sketches to illustrate them (Fig. 1c,d), and the different language both communities use.

395
To retrieve "reliable" cross-correlation functions of the recorded data in ambient noise seismology, Above, we have already explored the advantages and limitations of using numerically computed 419 synthetic wave fields (Fig. 2) and amplitudes (Fig. 3) in MFP, as well as the emergence of striped 420 interference patterns for narrowband sources (Fig. 2). MFP shows further undesired behaviour 421 under certain conditions that we encounter in real-world applications. Some of these are more 422 non-peer reviewed preprint submitted to EarthArXiv straight-forward to understand in the conceptual framework of convolution introduced above. For beamforming, the impact of an array's geometry on its resolution capability is well studied, and 425 quantified by the array-response function (Rost & Thomas, 2002 beamforming is an unsolved problem, because the parameter space is much larger.

437
Closely related to these considerations is that high waveform coherency is required across stations, Station density has direct impact on the retrieved beampower distribution that is worth pointing out 451 explicitly (Fig. 4a). In a synthetic test, we place four times the stations, with the same coordinates, 452 non-peer reviewed preprint submitted to EarthArXiv on the right side than we do on the left (Fig. 4b). The beampower distribution shows a visual bias does not seem to be desirable here. Single sources can cause prominent interference patterns for narrowband sources (Fig. 2c,d), which 470 depend on station geometry and frequency band. This leads to even more complex, secondary 471 interference when multiple sources are active at the same time. In a synthetic test, we place two 472 narrowband sources that excite identical wave fields simultaneously (Fig. 4b). The second source 473 is placed such that it lies at the edge of a sidelobe of the first source (Fig. 2d). From the retrieved 474 distribution of beampowers it is not at all obvious that two and only two sources are active here, 475 and instead this may be misidentified as a single source close to the left array (Fig. 4b). The  Hills earthquake of 2008-07-29 (Fig. 5). When applying the standard MFP approach, with an 513 assumed velocity v = 3.2 km/s (the best fit in the synthetic test in Fig. 2), we find a relatively good 514 location of the earthquake with 7.7 km distance to the location in the CI catalog (Fig. 5a, SCEDC, explosive source mechanism, we at first find a decrease in location accuracy (Fig. 5b). The retrieved 519 location is 18.3 km away from the CI location. When we incorporate the moment tensor solution 520 from the CI catalog (SCEDC, 2013), trivial to do with our approach, we find an improvement in 521 location accuracy with a distance of only 1.9 km to the CI location (Fig. 5c). This demonstrates one

536
On first order, we find a good match between the standard approach (v = 3.2 km/s, left), our 537 approach (middle), and the distribution of significant wave height (right). A good match between 538 seismic wave excitation and ocean wave activity is expected for the secondary microseism. The com-539 mon explanation is that ocean gravity waves at the water surface, propagating in roughly opposite 540 direction, interact and cause a standing wave that generates a vertically-propagating pressure wave 541 non-peer reviewed preprint submitted to EarthArXiv field in the water column. This pressure wave field then interacts with the ocean bottom, generating 542 seismic waves in the solid Earth (Hasselmann, 1963;Ardhuin et al., 2015).

543
The similarity between the standard approach and our approach is high (Fig. 6 left and middle). Hills earthquake (Fig. 5), improving the accuracy of the synthetic wave field has larger impact. for which we believe our approach should perform a lot better than the standard approach. Still, 554 we do find that beampower distributions retrieved with our approach are more focused on specific 555 regions compared to the standard approach. We do not yet feel comfortable in judging whether 556 these differences are certain to be an improvement in source estimation.

557
Our synthetic tests (Fig. 2) and the Chino Hills earthquake example (Fig. 5) suggests that 558 our approach can be more precise in locating the sources. However, we have to be careful with 559 interpreting these patterns, as we have also demonstrated in synthetic tests (Fig. 4). If our approach 560 will prove to be more precise also for microseisms, we may find that seismic waves are excited in 561 specific regions in the oceans and not distributed homogeneously beneath storm systems. It is 562 important to note here that for now we use an explosion source mechanism for the synthetic wave 563 fields to locate the microseism, which we have already shown to be inadequate for an earthquake 564 (Fig. 5). In the future, we require a strategy to describe and incorporate a source mechanisms 565 appropriate for microseisms. Some insight in how that could be approached has been given by the need for the plane-wave assumption. It is one of the current approaches to locating sources of 571 non-peer reviewed preprint submitted to EarthArXiv ambient seismic noise (Fig. 1). In this study, we advance MFP to better incorporate elastic wave 572 propagation in the Earth by using Green's functions numerically computed for real Earth structure 573 directly as the synthetic wave field that the data is matched against.

574
When amplitudes are considered in MFP, results are biased by amplitude effects such as geo-575 metrical spreading and anelastic attenuation. In the standard approach, this is usually neglected 576 through spectral whitening of the synthetic wave field. We find that this strategy performs best for 577 us as well, and that trying to correct for spreading and attenuation via an amplitude term, as has 578 been suggested before, may not be advisable (Fig. 3). This is especially the case for our approach, 579 where multiple wave types can be considered simultaneously.

580
Two examples on real data showcase the potential of our approach (Figs. 5, 6). In principle, we 581 can use it to search for the source mechanism of a seismic source, as suggested by the improved source 582 location after incoporating the earthquake's moment tensor (Fig. 5). This could be particularly 583 useful in the context of tremor activity, where source mechanism determination is challenging with 584 classical approaches. In a second example, we locate sources of the secondary microseism in the 585 Northern Atlantic and Mediterrenean Sea (Fig. 6). Results from our approach match the standard 586 approach's results closely, likely due to source geometry, narrow frequency band, and our reliance on 587 Green's functions computed for an axisymmetric Earth. Our approach retains the advantage that is 588 not biased by author choice of a constant velocity, and potentially provides higher resolution. 589 We clarify conceptual approaches to MFP and its close connection to the interferometry-based 590 localisation. The striking similarity between them suggests that it may be a worthwhile endeavour to 591 unify them in the future, or at least provide a framework to let the different communities benefit from 592 each others' work. MFP in particular would benefit tremendously from an approach for quantifying 593 its resolution. The lack of such a measure is currently its major disadvantage.  non-peer reviewed preprint submitted to EarthArXiv Standard approach. Emergence of sidelobes due to interference. d) Our approach in the same narrow frequency band. Source is precisely located. non-peer reviewed preprint submitted to EarthArXiv  km distance to the CI location. b) Beampower distribution using numerical Green's functions for an explosive source mechanism. 18.3 km distance to the CI location. c) Beampower distribution using numerical Green's functions for the moment tensor solution in the CI catalog (SCEDC, 2013). 1.9 km distance to the CI location. non-peer reviewed preprint submitted to EarthArXiv