Sensor networks to measure environmental noise at gravitational wave detector sites
Koley, S.
2020
document version
Publisher's PDF, also known as Version of record
Link to publication in VU Research Portal
citation for published version (APA)
Koley, S. (2020). Sensor networks to measure environmental noise at gravitational wave detector sites.
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal ?
Take down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
E-mail address:
Passive seismic data analysis
methods
1
2.1
Passive seismic studies - A historical context
Studies of seismic waves can be categorized broadly into two types. The first one deals with studying the propagation of seismic waves generated by either artificial sources like explo-sions or by introducing vibrations into the subsurface with a shaker that is coupled to the Earth’s surface. The advantage of using such sources is that the phase and the amplitude of the seismic source is known. Seismic waves generated from these sources gets reflected and transmitted at the subsurface layers, and are then received on the surface by sensors which are capable of measuring tiny vibrations on the surface of the earth. These seismic studies are referred to as active studies as they involve using an external seismic source to activate the Earth’s subsurface. It is the body waves generated by these sources that are analyzed and used to infer the propagation characteristics of the medium, like the P-wave and the S-wave velocity. While these controlled-source seismic methods flourish due to the clarity and accuracy of the resultant subsurface image, there are some caveats to the method. Due to the high-frequency content of the source signal (greater than 1.5 Hz), these methods can im-age the subsurface upto a few kilometers and are not suitable for deep crustal level imaging. Moreover, exciting the Earth’s surface with an external seismic source might not be always possible due to the intrusive nature of the method, especially in urban environments. Hence, passive seismic techniques for imaging the subsurface are much more desirable in such situ-ations.
Passive seismic methods use the seismic noise produced by natural and anthropogenic sources to image the subsurface. Seismic waves generated by ocean waves hitting against a shore or those generated by earthquakes are some examples of naturally occurring seismic noise sources. Artificial sources of seismic noise are due to human activities. Cars passing by on a road, ground vibrations created by farming activities or vibrations induced into the ground
1The contents of this chapter have been partially published as ‘Innovations in seismic sensors driven by the search
by road bridges are some of the examples of anthropogenic seismic noise. Seismic waves generated from these sources propagate through the Earth’s surface both as body and surface waves. In active seismic studies, the surface waves generated from the seismic sources are treated as noise, and hence filtered out. In passive seismic studies, we measure the ground motion due to surface waves propagating along the surface of the Earth and use it to infer information about shallow subsurface geological structures. The dominant contribution of surface waves to the seismic wavefield can be due to both the source mechanism and the sta-tionary phase arguments about the noise source. Firstly, for observation points on the Earth’s surface, mostly the surface wave part of the Green’s function is retrieved from ambient-noise sources located at or near the surface, whereas body wave response retrieval requires a distri-bution of noise sources in the depths below the receivers [50]. Secondly, the lack of stationary points for body wave generation also yield weaker body wave retrievals, because it is the sta-tionary part of the wavefield that contributes the most to the retrieved Earth response between receivers [51]. Since the surface waves propagating along the surface of the Earth can be treated as plane waves, array techniques for computing the velocity and azimuth of propaga-tion of these surface waves are well suited [52].
Array studies of microseisms have gained much popularity in recent years, especially after the works of Bromirski et al. in 1999 [53] who showed the usefulness of acquiring micro-seismic data as a substitute to buoy data for monitoring oceanic activities. However, the most important utility of measuring seismic noise on the Earth’s surface with an array of seismic sensors was only known after the work of Campillo & Paul, 2003 [54], Shapiro & Campillo, 2004 [55], and Gouedard et al. in 2008 [56]. They demonstrated the method of extracting the Green’s function of the Earth along the path joining two seismic sensors by cross-correlating simultaneously recorded long duration seismic noise signals measured at the two sensor loca-tions. Work by Wapenaar, 2004 [57] and Wapenaar & Fokkema, 2006 [58] were instrumental in formulating the theoretical framework behind the feasibility of extracting the complete elastodynamic Green’s function of the medium by cross-correlating signals recorded simul-taneously by two seismic sensors at different locations. A case study of the Green’s function extraction and subsequent surface wave tomography is presented in Chapter 5 of this thesis. An important assumption in Green’s function retrieval is that of an isotropic distribution of the seismic noise sources around the pair of sensors. In reality this is never the case. As a result, an understanding of the direction of propagation of the surface waves is necessary in order to quantify the accuracy of the extracted Green’s function of the propagating medium. Array studies of ambient seismic noise refer to the analysis of simultaneously recorded seis-mic noise at different receiver locations in order to obtain the direction and velocity of the surface waves propagating through the array of sensors. Due to the dispersive nature of sur-face waves, as was discussed in Chapter 1, the estimated velocity of propagation varies with frequency. Array studies of microseismic noise [53] dominant in the frequency band 0.05 to 0.3Hz have been previously reported in the work of Yao et al. in 2008 [59], and Kimman et al. in 2012 [60]. Since the wavelength of the surface waves originating from microseisms can range in the order of tens of kilometers, large arrays with sensor separation of at least three times the maximum wavelength of interest are desired for satisfying the far-field criterion of the incident surface waves. In the case of high-frequency surface waves (greater than 1 Hz), seismic sensor array sizes are much smaller (order of kilometers) as the wavelength of high-frequency surface waves are relatively smaller. Some of the array studies of high high-frequency
seismic noise have been reported by Hanneman et al. in 2014 [61], Wang et al. in 2016 [62], Koley et al. in 2017 [63] and Koley et al. in 2018 [64]. Since high frequency surface waves attenuate faster, array analysis of such surface waves is much more difficult and the surface waves recorded between distant sensors often lack coherence. To summarize, the objective of array studies of seismic noise is to first obtain an understanding of the frequency depen-dent propagation characteristics since sources in different frequency band may originate at a different azimuth, and also estimate the dispersion surface waves. The dispersion of surface waves can be further used to estimate a 1D S-wave subsurface velocity model of the region. This chapter starts off with a discussion of wireless seismic sensors used for measuring the seismic ground motion, followed by a description of the array analysis techniques for surface waves like beamforming and extended spatial autocorrelation (ESAC). These methods are used in subsequent chapters of the thesis for analyzing ambient seismic noise data. We then briefly discuss the method used for estimating a 1D S-wave velocity model from the surface wave dispersion curve. Each theoretical section of this chapter is complemented by citing results from a passive seismic survey that was conducted in Groningen, Netherlands. The example passive seismic case study presented in this chapter concerns only Rayleigh waves, as only the vertical component of the seismic noise was measured. Throughout the chapter, whenever we refer to surface waves, we refer to Rayleigh waves, unless otherwise mentioned.
2.2
Wireless seismic sensors
Cabled sensor networks pose several restrictions on flexible seismic array design. Especially in the context of passive seismic surveys that require dense sensor sampling at some places and sparse at others, cabled sensor networks are not desired. In this section we discuss the features of the state of the art wireless seismic sensors that were designed in collaboration with Shell2and Innoseis3. All passive and active seismic data acquired and later used in this
thesis were measured by these wireless seismic sensors.
2.2.1
Sensor overview
The required performance of existing seismic sensors exceeded what was cost effectively available in the market to sample the seismic wavefield with sufficient density. This set in motion a development of novel sensor technology as well as a means to quickly and cost ef-fectively measure seismic vibrations with high density flexible arrays. The sensors were then commercialized by a Nikhef spin-out named Innoseis. Innoseis leveraged the technology and knowledge from fundamental gravitational-wave research to develop an ultralow-power, lightweight, and long-endurance autonomous recording system called Tremornet. Fig. 2.1(a) shows the different sensor components; an assembled version is shown in Fig. 2.1(b). The major challenges in developing autonomous sensors that can acquire data continuously for months are: GPS synchronized time stamping while ensuring low-power usage and also no
2Shell Global Solutions International B.V. 3Innoseis B.V.
GPS antenna
Geophone
Sensor casing
Electronics (FPGA, MCU, ADC, GPS)
Batteries
Spike
(a) (b)
Figure 2.1:(a) A overview of the different components of the Tremornet seismic sensor. (b) Tremornet seismic sensor with all components assembled.
compromise on the sensitivity of the sensors. Moreover, noise of the sensor read-out electron-ics must be low, especially at frequencies below the resonance frequency of the geophone. The sensor’s front end comprises a geophone and an ADC (analog to digital converter) along with an FPGA (field-programmable gate array) implementing CIC filters for the purpose of decimation and sequential combination of operations. The sensor is equipped with a GPS module for acquiring GPS time-stamps at desired time intervals. The GPS location accu-racy is within a meter. At present 4.5 Hz geophones are used, but the choice of geophone completely depends on the purpose and the type of measurement desired. For example in mi-croseismic studies where we want to measure the ground motion below 0.1 Hz, the resonance frequency of the geophones need to be low (about 1 Hz) in order to have sufficient sensitivity below 0.1 Hz. On the contrary, in case of data acquisition applications for active seismic imaging a resonance frequency of 4.55 Hz would be sufficient as the frequency content in such data is in the range 5 to 90 Hz.
Since the sensors are wireless, they can be easily deployed. Data acquisition commences as soon as the sensor is planted vertically on the ground and it has obtained a GPS fix4; typically
accomplished within a minute. Initially, the data is acquired by the ADC at a high sampling rate of upto 1, 024 kS/sec. It is then downsampled to a desired sampling rate using a com-bination of anti-aliasing and CIC filters, which are implemented in the FPGA. The FPGA is suited to perform fixed tasks in a parallel manner at high frequency, and is less suited for complicated instruction sets that might be necessary for other operations. For this purpose, an additional 32-bit RISC MCU (micro-controller unit) is used. The FPGA runs on a 32 MHz crystal which has maximum deviations of 1 ppm every second. Hence, the time stamps given by the FPGA need to be synchronized to GPS time-stamps in order to avoid erroneous time stamping of the acquired seismic data. However, for low-power operation to achieve
long battery lifetime, the GPS needs to be switched-on for only about 0.1% of the total run time. Hence for the time stamping to be synchronized within 5 µs of the GPS time-stamp, methods were developed that can correct for this deviation in time stamping. The timing correction algorithm was tested with various 32 MHz crystals, and for different temperatures in the range −20◦C, to 70◦C. We achieved a timing accuracy of better than 5 µs and a GPS
awake time of less than 0.5% of the total run time. In the process we also tested several of the GPS modules for attaining the desired location accuracy of one meter, and u−blox5modules
were finalized for use.
2.2.2
Vault noise test
After extensive testing in the lab, verification of the Tremornet sensors was done in cooper-ation with the Royal Netherlands Meteorological Institute (KNMI). Seismic measurements were performed at the KNMI Heimansgroeve seismic observatory near Epen in the Nether-lands from 11 to 25 November 2015. The objective of the experiment was to verify the performance of Tremornet sensors by comparing the recordings with the KNMI Streckeisen STS-1 seismometer. The sensor settings like the sampling rate, gain and the full-scale input voltage was determined in order to have a low input noise while ensuring a high dynamic range. Fig. 2.2(a) and (b) show the dynamic range of the Innoseis sensors for five differ-ent gain settings of 0, 6, 12, 18 and 24 dB corresponding to sampling rates of 2 and 4 ms. Between gains of 0 and 18 dB the noise magnitude halves as the full-scale input voltage is also reduced by a factor two. However, the decrease in the noise-level is not significant while increasing the gain from 18 to 24 dB. Hence a gain of 18 dB was used which also ensures a higher dynamic range as compared to a gain of 24 dB.
The relatively low seismic background noise and stable measuring conditions allow the ap-plication of the three-channel correlation technique (Sleeman et al. in 2006 [65]) in order to compute the frequency dependent self-noise of the seismic sensors. This is particularly useful when verifying the self-noise performance of seismic sensors that have an integrated architecture. Fig. 2.3(a) gives a summary of the results. The spectra are calculated from 20 s stretches of data that were windowed with a Blackman-Harris window and then averaged. The transfer function of the respective sensors was used to convert the sensor output to equiv-alent ground motion input. A comparison can be made up to 12.5 Hz, as this is roughly the bandwidth of KNMIs STS-1 data.
In Fig. 2.3(a), we see excellent agreement among the various sensors. It can also be observed that the three-channel correlation technique allows us to accurately extract the sensor self-noise levels. These results show that the sensor self self-noise is less than 1 ng/√Hz at 1 Hz and as low as 0.1 ng/√Hz at the sensor resonance frequency of 4.5 Hz. A plot of the endurance of the sensor’s remaining battery charge and recorded node temperature versus time is shown in Fig. 2.3(b). With a setting of 250 ms sampling rate and a highest gain of 24 dB, the sensors can record continuously for a period of approximately seven weeks.
0 6 12 18 24 Gain (dB) 10-8 10-6 10-4 10-2 100 102
Input voltage range (V)
(a) Sampling rate = 2 ms
0 6 12 18 24 Gain (dB) 10-8 10-6 10-4 10-2 100 102 (b) Sampling rate = 4 ms
Figure 2.2:Dynamic range of the Tremornet recording system based on the full-scale, gain and equivalent noise specifications for sampling rate of (a) 2 ms and (b) 4 ms.
(a) (b)
Figure 2.3: (a) Self-noise measurement of three Tremornet sensors compared against the KNMI STS-1 seismometer. The Peterson’s low and high-noise models are indicated for com-parison. (b) The temperature and the battery status of a Tremornet sensor during a 21 day survey.
2.3
Data Analysis Methods
From a geotechnical standpoint, determination of the S-wave velocity structure of the sub-surface is important for accurately estimating the strong ground motion from earthquakes or explosions. In gravitational wave (GW) community, estimation of the S-wave velocity model is necessary for simulating seismic ground motion, both on the surface of the Earth and in
the subsurface. Since seismometers cannot be positioned at all points near a GW detector, simulations are the only way for estimating the ground motion near a detector. There are two objectives of simulating the ground motion. Firstly, we can compute the ground motion in the subsurface even without having to install a seismometer in a borehole, which is a rather expensive and cumbersome exercise. This information is extremely valuable especially in the context of third-generation GW detectors like Einstein Telescope [18] that are to be built un-derground. Secondly, by using the simulated ground motion we can compute the contribution of the Newtonian noise [66] on the suspended elements of the detector. Newtonian noise is one of the noise sources that limits the low-frequency sensitivity of the current ground based GW detectors and is also foreseen to limit the sensitivity of the third-generation underground GW detectors [21].
In order to simulate the ground motion, we need information about the geology of the region that includes the P-wave velocity, S-wave velocity, density and the quality factor models of the subsurface. Solutions for the subsurface displacement-stress can then be accomplished following Eq. (1.59) stated in Chapter 1. One way of estimating the subsurface parame-ters is by computing the Rayliegh wave dispersion from surface measurements of seismic noise. A detailed formulation of the dependence of subsurface parameters on the Rayliegh wave dispersion was discussed in section 1.2.5 of Chapter 1. The frequency-wavenumber method (f-k) and the ESAC method are the two data processing techniques used to obtain the surface wave phase velocities from ambient seismic noise measurements. The f-k method is also sometimes referred to as beamforming [67]. Beamforming has been widely used in studies by various researchers, for example, Asten & Henstridge, 1984 [68]; Horike, 1985 [69]; Okada et al. in 1987 [70]. On the contrary, the ESAC method, which can estimate surface wave phase velocities more accurately than beamforming [71], was seldom used by geophysicists due to restrictions on the sensor layout and because most of the work has been reported in Japanese. It was not until the work of Ohori et al. in 2002 [71] who reported the work previously done by Okada et al. in 1987 and Chouet et al. in 1998 [72] that the method gained popularity. The ESAC method is an extension to the SAC method originally proposed by Aki in 1957 [73]. However, due to the limitation that the SAC method could only be used for circular arrays, it has been used only by a few researchers. Ohori et al. , 2002 [71] extended the SAC method such that it could be used for arbitrary array geometry and named it ESAC. The ESAC method assumes that over a long recording time of seismic noise, the sources of noise illuminate the array isotropically. The ESAC method does not provide any information about the direction of propagation of the surface waves, unlike the f-k method which is slightly advantageous in this regard. In the next subsections we discuss the two passive seismic data analysis methods and demonstrate their applicability by using the seismic noise data from the Groningen array.
2.3.1
Beamforming
Consider a 2D array of seismic sensors as shown in Fig. 2.4(a). Let a body wave traveling along an azimuth φ measured clockwise from the vertical z direction be incident on the array as shown in Fig. 2.4(b). The component of the body wave on the horizontal x-y plane makes an angle θ, measured anticlockwise from the y-axis.
0 200 400 600 800 1000 1200 1400 1600 0 400 800 1200 1600 x-axis (m) y-axis (m ) (a) (xm,ym) u uy ux ua O y z x (b) sensor location
Figure 2.4:(a) An example of a 2D non-uniform surface array with the cartesian coordinate of the mthsensor represented as (x
m, ym). (b) Decomposition of the body wave slowness
vector ~u.
Let the velocity of the body wave be ~v and the slowness be ~u such that |~u| = 1 |~v|, ~ u |~u| = ~v |~v|. (2.1)
The component of the slowness vector ~u on the horizontal x-y plane is denoted ~ua. Vector ~ua
can be further decomposed along the x and y-axis with the respective components denoted by ~uxand ~uy. Then we can write
ua = | ~ua| =
1 |~v|sin φ, ux= | ~ux| = uasin θ,
uy = | ~uy| = uacos θ. (2.2)
In practice, the vertical angle of incidence φ in case of body waves is unknown and we only estimate the apparent velocity | ~ua|. For surface waves, the problem is much simple, as
φ = 90◦and the apparent velocity equals the true velocity of propagation. Data model vector
Let a seismic noise signal represented as s(x, t) be received by an array of M sensors, which are laid out on a 2D surface, mimicking the surface of the Earth. Under the assumption that the seismic signal can be can be described as a superposition of plane waves, the time delay
τmfor the signal to reach the mthsensor with position vector ~rmwith respect to some chosen
origin can be written as
τm= ~u. ~rm. (2.3)
Then the signal as received by the mthsensor can be expressed as
s( ~rm, t) = s(0, t − ~u. ~rm) (2.4)
where we have just applied a time shift to the source signal s(0, t). If we perform a Fourier decomposition of the signal s( ~rm, t), then we get
F (s( ~rm, t)) = 1 2π Z ∞ −∞ s( ~rm, t)ejω(t−~u. ~rm)dt, = 1 2π Z ∞ −∞ s( ~rm, t)ejω(t−τm)dt, = 1 2π Z ∞ −∞ s( ~rm, t)ejωt(e−jωτm)dt, (2.5)
where j is the imaginary number. From Eq. (2.5) we can infer that in the Fourier domain, the signal received by each sensor can be simply obtained by applying a phase shift to the source signal. The phase shift e−jωτm is again a function of the receiver coordinates (x
m, ym).
Given this benefit of processing the signal in the Fourier domain, from hereon we perform all the data processing in the Fourier domain.
In practice, for any given time window of seismic ground-motion measurement, there may be several sources that contribute to the measured seismic noise. Hence the seismic signal measured by each sensor when expressed in the time domain is a superposition of the ground motion from each of the sources. At first, it might appear that the scheme for processing such signals is complicated. However, under the assumption that every source is uncorre-lated, meaning (s1(t) ? s2(t))(τ ) = 0where ? denotes the correlation operator and τ is the
correlation lag variable, the processing of such signals is much simpler. Hence, we proceed by writing out the theory of beamforming for one source first and then we extend it for the case of multiple uncorrelated sources.
Let a source signal represented as s(0, t) be received by an array of M sensors such that in a desired time window ∆t, the signal is observed by all the sensors. The width of the time window ∆t should not matter from a theoretical standpoint because we split a long timeseries into several of such windows and the consecutive windows have an overlap such that no data are left unused. Now, for processing the data measured in one such time window we first de-fine the array response vector AR(f). It is a column vector of M elements, and each element corresponds to the frequency domain phase shift associated with a particular sensor location. Hence, we can write the array response vector AR(f) as
AR(f) = e−2πjωτ0 e−2πjωτ1 . . . e−2πjωτM −1 (2.6)
where τm−1 = xmu sin θ + ymu cos θis the delay of the signal to reach the mth sensor
location from the origin, and (u, θ) correspond to the slowness of the propagating signal and the azimuth measured clockwise from the y-axis (Fig. 2.4(b)) respectively. If the Fourier transform of the source signal s(0, t) is S(f) then the theoretical data model vector Xth(f )
in the frequency domain can be written as
Xth(f ) = e−2πjf τ0 e−2πjf τ1 . . . e−2πjf τM −1 S(f ). (2.7)
Each element xi(f )of the column vector Xth(f )is the Fourier transform of the signal
re-ceived at the ithsensor and at frequency f. Assuming some noise is present in the signals
within the measurement time window ∆t, we can finally write the observed data model vector Xobs(f )in the frequency domain as
Xobs(f ) = Xth(f ) + N(f ), (2.8) where N(f) is the Fourier component of the observed noise at frequency f during the time of measurement. Similar to Xobs(f ), N(f) has a total of M elements which is the same as
the number of sensors.
Extending the above formulation for N sources acting on an array of seismic sensors, we can write the data model Xmo(f )vector in frequency domain as
Xmo(f ) =AR1(f )S1(f ) +AR2(f )S2(f ) + ... +ARP(f )SP(f ) + N (f ), (2.9)
where ARi(f )corresponds to the array response vector for the ithsource as stated in Eqs.
(2.6) and (2.7). Eq. (2.9) can also be written as a matrix equation of the form
Xmo(f ) =AM(f)S(f) + N(f), (2.10)
where AM(f) = AR1(f ) AR2(f ) ... ARN(f )
, S(f) = S1(f ) S2(f ) . . . SN(f ) .
The matrix AM(f) is also referred to as the Array Response matrix. Theoretical Array Response
Before acquiring ambient seismic noise in the field, designing an appropriate geometry for the seismic sensor array is vital. In this subsection we discuss the concept of Theoretical Array Response that is used throughout the thesis for designing a seismic sensor array that is capable of measuring the direction and slowness of the propagating surface waves in a given
frequency band. Due to the dispersive nature of surface waves, the slowness or the velocity of the surface waves change with frequency. As was explained in the previous chapter, we know that this dispersive nature of the surface waves is a function of the subsurface layer depths, velocities (P and S-wave) and density. Moreover, the amplitude of surface waves decay exponentially as a function of their depth of propagation. As a rule of thumb the sen-sitivity of surface waves deteriorates below a third of the dominant wavelength [74]. Hence we can draw the conclusion that surface waves with larger wavelengths will penetrate deeper in the subsurface and vice versa for smaller wavelengths. Now, given that the P-wave and the S-wave velocity of the subsurface layers increase with depth, we can infer that the velocity of higher wavelength surface waves is larger than that of short wavelength surface waves which are confined to shallower layers of the subsurface.
The idea behind using an array of seismic sensors in the field is to be able to spatially sample the seismic wavefield propagating through the survey area. Similar to sampling signals in the time domain where an ADC is used to digitize the signal at a fixed sampling rate, array processing can be visualized as sampling the seismic wavefield spatially with a network of sensors distributed on the surface of the Earth. Unlike sampling time-domain signals where a signal is sampled only along the time axis, array processing of surface waves involve sam-pling the seismic wavefield along two orthogonal directions that span the horizontal plane of the Earth’s surface. For time-domain sampling, the Nyquist criterion forms the basis for selecting the interval between two successive time samples. Similarly, for spatial sampling of surface waves it is the smallest distance between sensors and their azimuth of orientation on the Earth’s surface that play a key role in deciding whether an array layout is suitable for spa-tially sampling the seismic wavefield. An obvious solution to this is to use an uniform square grid of seismic sensors with a minimum inter-sensor spacing of less than half the smallest wavelength of the surface waves in a given frequency band. However, the number of seismic sensors that would be necessary for executing such a survey is huge and not realistic. Hence, we need to be able to design seismic arrays more efficiently while ensuring that there is no spatial aliasing at high frequencies and a good resolution for low-frequency surface waves. Selection of the frequency band of interest is also important at this stage, as Rayleigh wave phase velocities increase rapidly at low frequencies and would need larger inter-sensor dis-tances.
Several researchers like Asten & Henstridge in 1984 [68] posed restrictions on the minimum and the maximum inter-sensor distances in an array, but there was no definite rule for design-ing an optimal array layout. We use the concept of theoretical array response as proposed by Woods & Lintz in 1973 [75] to determine if an array is suitable for estimating the direction and velocity of propagation of surface waves in a desired frequency band. For a fixed sensor layout and at a fixed frequency f, the array response vector defined in Eq. (2.7) is a function of the slowness u and the direction of wave propagation θ. The theoretical array response ThAR(uK, θK)for a given array layout and corresponding to a fixed frequency surface wave
propagating with slowness uK and along an azimuth θKis then computed as
ThAR(uK, θK) = 1 M[AR(uK, θK)] ∗[AR(u 1, θ1),AR(u2, θ2), ...,AR(uP, θP)], (2.11) = [AR(uK, θK)]∗AMtrial, (2.12)
where ‘*’ denotes the complex conjugate transpose operator and M is the number of sensors in the array. Each column of matrix AMtrialcorresponds to an array response vector for a
trial value of slowness and azimuth. What we perform in Eq. (2.11) is the inner product of the array response vector ARuk,θK with several other array response vectors corresponding
to different values of (u, θ). This is also the reason why we term the matrix AMtrial as a
trial array response matrix, because it comprises the array response vectors for all possible combinations (u, θ) that we consider. For example, if we choose to discretize the azimuth domain beyween 0◦and 360◦)at an interval of 0.5◦and the slowness domain between 0.001
and 0.01 s/m at an interval of 0.0001 s/m, then the total number of trial array response vectors we get is 720 × 91.
The resulting theoretical array response vector ThAR(uK, θK)is a row vector with P
ele-ments, the same as the number of the trial array response vectors used. It is clear that the inner product of AR∗(u
K, θK)with the appropriate trial array response vector AR(uL, θL)
will result in unity only if uK = uL and θK = θL. All the other inner products will be
less than unity. Thus, what we have now is a set of values of the array response for every trial combination (u, θ). Our objective is to design a seismic array, such that for surface waves with given slowness uK and direction of propagation θK, a unique peak is observed
in the (u, θ) domain and only at (uK, θK). There may be two observations in the computed
ThAR(uK, θK). Firstly, the peak observed in the (u, θ) graph may be too broad and hence
lack resolution. It implies that the uncertainty associated with detecting a peak at (uK, θK)
is greater. This is typically the case for low-frequency surface waves if the maximum aper-ture of the array is not large enough as compared to the largest wavelength of interest. The maximum aperture must be at least the same as the largest wavelength of the surface waves in the given frequency band. Secondly, multiple peaks might be observed in the (u, θ) graph in addition to the one that we expect at (uK, θK). This is the case for spatial aliasing and is
ob-served for high-frequency surface waves when the smallest inter-sensor separation is greater than half the smallest wavelength of interest. It is also worth noting that the behavior of the theoretical array response ThAR(uK, θK)might not be the same as ThAR(uK, θM). Although
the slowness of the approaching surface wave is the same in both cases, it is the azimuth of the propagating surface wave that makes the difference. This case is often observed if the designed array is azimuthally asymmetric. Hence, designing an array that has azimuthal symmetry is desired.
We now put forth an example of the computed theoretical array response for the seismic array comprising 96 wireless seismic sensors that was deployed in Groningen for a period of three weeks between January 5 and 23, 2016. This passive seismic study was conducted in collab-oration with Shell6 and Innoseis7. The objective of the survey was to test the performance
of the TremorNet wireless seismic sensors developed by Innoseis7alongside standard cabled
short period seismometers. The array recorded the vertical component of ambient seismic noise continuously for a period of three weeks and hence was suitable for carrying out pas-sive seismic data analysis.
Fig. 2.5(a) shows the array layout on a cartesian plane. Seismic sensors were deployed on a regular grid with a spacing of 100 m between consecutive sensors. The maximum inter-sensor distance was about 1200 m and the minimum was 100 m. Fig. 2.5(b) shows the distribution
6Shell Global Solutions International B.V., Netherlands 7Innoseis B.V., Netherlands
0 500 1000 1500 Inter-sensor separation (m) 0 100 200 300 400 500
Sensor pair count
(c)
0 500 1000 1500
Distance along longitude (m) 0
500 1000 1500
Distance along latitude (m)
(a) 0 500 1000 1500 Inter-sensor separation (m) 0 50 100 150 200 One-way azimuth ( 0) (b) East
Inter-sensor separation bin size = 50 m
Figure 2.5:(a) Seismic array layout at Witteweirum, Groningen shown in cartesian coordi-nates, the one way azimuth θ is measured anticlockwise from the East. (b) Distribution of one-way azimuths and inter-sensor separation corresponding to all sensor pairs. (c) His-togram of the number of sensor pairs corresponding to all inter-sensor separations.
of the inter-sensor distance and the one-way azimuth for all sensor pairs. A well spread distri-bution of inter-sensor distances is observed for all azimuths upto a maximum distance of 800 m. For inter-sensor distances greater than 800 m, the azimuth distribution is good only along the azimuths of 45◦and 135◦which corresponds to the diagonals of the square grid. Sensor
pair counts in excess of 200 are observed for the inter-sensor distance interval of 100 to 800 m, and the pair count gradually tapers off to values of less than 50 for greater inter-sensor distances (Fig. 2.5(c)).
Given the array layout, the next objective is to compute the theoretical array response values following Eq. (2.11). However, before computing the theoretical array response we need to decide on the values of few parameters. They are the slowness and the azimuth of the incident test plane wave and also the slowness-azimuth ranges for constructing the trial array response matrix AMtrial. We consider a plane wave propagating at three different frequencies
of 0.4, 1.0, and 1.7 Hz and impinging the array from an azimuth of 45◦. The slowness at
each of the three frequencies are selected on the basis of previously observed Rayleigh wave phase velocities in the region [60]. The set of trial array response vectors as in Eq. (2.11) are computed in the slowness range 0.00025 to 0.005 s/m at an interval of 0.0001 s/m and in the azimuth range 0 to 359.5◦at an interval of 0.5◦. Hence, the theoretical array response value
is computed corresponding to a total of 48 × 721 trial array response vectors. Fig. 2.6 shows the computed theoretical array response values on a polar graph with slowness u increasing radially outwards and azimuth θ measured clockwise from the North (y-axis). Consequently, the x and the y-axis are labeled as u cos θ and u sin θ respectively. As stated previously, key observations from the theoretical array response plots in Fig. 2.6 are lack of resolution of the observed peak at 0.4 Hz and spatial aliasing at a frequency of 1.7 Hz. The lack of resolution of the peak in the u − θ domain corresponding to the plane wave of frequency 0.4 Hz is due to insufficient array aperture (maximum of 1200 m) as compared to the wavelength of the plane wave at 0.4 Hz which is about 3000 m. The smallest sensor separation in the array is 100m, hence it can only sample plane waves with wavelengths of 300 m or more. As a result for our test case of a plane wave of frequency 1.7 Hz and propagating at a speed of 330 m/s and wavelength of about 194 m, we observe spatial aliasing. Fig. 2.6(c) shows the two
incor-rect peaks observed due to spatial aliasing in the u − θ graph alongside the true peak which is observed at an azimuth of 45◦and a slowness of 0.003 s/m. Thus we can conclude that
with this particular array configuration, we can only estimate the speed of the Rayleigh wave propagation upto a maximum frequency of 1.7 Hz. In practice it turned out that the upper limit in frequency was even lower than 1.7 Hz, since the phase velocity of Rayleigh waves was smaller than 330 m/s. In the frequency regime below 0.4 Hz, the associated uncertainty for obtaining the correct Rayleigh wave speed with a 95% confidence is of the order of a few hundred m/s.
Figure 2.6: Theoretical array response for a plane wave impinging the array of seismic sensors shown in Fig. 2.5(a) corresponding to (a) velocity of 1200 m/s, propagation azimuth of 45◦ at a frequency of 0.4 Hz, (b) velocity of 570 m/s, propagation azimuth of 45◦ at a
frequency of 1.0 Hz, and (c) velocity of 330 m/s, propagation azimuth of 45◦at a frequency
of 1.7 Hz.
Estimating beampower
Recall from Eq. (2.11) that the data model vector Xmo(f )can be written as Xmo(f ) =
AM(f)S(f) + N(f). The vector S(f) comprises the Fourier component of the multiple seismic sources at a particular frequency f. Now from Eq. (2.10) if we want to recover the source signal SL(f )by beamforming, we multiply the weight vector wL = M1ARL(f )and
the data vector Xmo(f ). It must be noted that although we represent Xmo(f )theoretically
in Eq. (2.11), in reality during beamforming the only quantity that is known to the data processor is Xmo(f ). The multiplication of the complex conjugate transpose of the weight
vector wLwith the data model vector Xmo(f )yields
w∗Xmo= 1 M AR∗ LAR1S1+AR∗LAR2S2+ ... +AR∗LARLSL+ ... +AR∗LARPSP+AR∗LN = SL+ 1 ME, (2.13)
where ‘*’ represents the complex conjugate transpose operator and E is the interference plus noise term given by
E =AR∗LAR1S1+ ... +AR∗LARL−1SL−1+ ... +AR∗LARL+1SL+1+ ...
+AR∗LARPSP+AR∗LN. (2.14)
The product in Eq. (2.13) is not a real number but rather a complex quantity. Hence, we multiply the product in Eq. (2.13) with its complex conjugate transpose to obtain a real quantity which is also referred to as the beampower corresponding to a particular weight vector wL. The beampower BP(uL, θL, f )corresponding to a particular weight vector wLis
then given by BP(uL, θL, f ) = 1 L(wLXmo)(w ∗ LXmo) ∗= 1 Lw ∗ LRxxwL (2.15)
where Rxx= XmoX∗mois the data covariance matrix and L is the number of windows into
which the data is split. Up to this point we have always treated the data model vector Xmo(f )
as a column vector where each element corresponds to the Fourier component of the data recorded by a particular sensor and at a particular frequency f. But when we compute the data covariance matrix Rxx, we construct our data model as a matrix where each column of
the matrix represents the data model vector for a given window. So if there are M sensors in the array and the data are split into L overlapping windows, then the data model matrix has the dimension M × L. It must be noted that the size of Rxxmatrix will always be M × M,
irrespective of the number of windows into which the data are split.
In Eq. (2.13) the source term SLis recovered correctly only when multiplied with the
appro-priate weight vector wL. Hence, the multiplication of the data model vector Xmo(f )needs
to be carried out with several weight vectors where each of them corresponds to a partic-ular value of slowness and direction of propagation. This is similar to what we did when computing the theoretical array response. The weight vector that maximizes the beampower BP(u, θ, f) gives us an estimate of the slowness and direction of the surface wave propagating through the array.
Ambient noise measurements at Groningen
Similar to Section 2.3.1 we use the passive seismic data measured in Groningen to demon-strate the applicability of beamforming to real ambient seismic noise recordings. Before proceeding with the beamforming results we briefly discuss the attributes of the ambient seismic noise observed in the region. Fig. 2.7(b) shows the power spectral density (PSD) of the seismic ground velocity measured by sensor 173 shown by the red dot in Fig. 2.7(a). Since the resonance frequency of the sensors is 4.5 Hz, the sensors are sensitive only down to frequencies of 0.2 Hz. Below 0.2 Hz, the noise of the electronics dominate the noise spec-trum. In the frequency band 0.2 to 0.7 Hz, a peak in the daily averaged PSD is observed, which corresponds to the secondary microseism due to oceanic disturbances in the North Sea. The magnitude of the peak varies within an order of magnitude during the 18 days of measurement, with maximum microseismic noise being registered on January 8, 2016. Fig. 2.8(b) shows a spectrogram of the seismic data recorded between January 5 and 11, 2016 by
sensor 173. The magnitude of microseism in the band 0.2 to 0.5 Hz peaks on January 8, 2016 and is correlated to the observed mean significant wave height (MSW) [76] in the North Sea during that day as shown in Fig. 2.8(c). A trough in the ground motion spectrum is observed between 0.7 and 1 Hz which agrees very well with global observations [77].
The frequency band 1.5 to 5 Hz is mostly dominated by ground motion generated by human activities and regional earthquakes of shallow origin. Two peaks, one centered at 3 Hz and the other at 5 Hz are observed during all days of measurement which corresponds to noise from a road running South-West of the seismic array. The noise level in this band varies within half an order of magnitude over all days of measurement. For clarity of the temporal variation of high frequency noise, we plot the spectrogram of the seismic data measured by sensor 173, where the temporal resolution of the spectrogram is one hour and the PSD for each hour of data is computed with a window length of 10 minutes and a 50% overlap between successive windows. Fig. 2.8(a) shows the daily variation of the PSD as a function of frequency and time. Since the high-frequency noise is dominantly man-made, PSDs estimated during the day are higher than that measured at night by an order of magnitude. At frequencies greater than 5 Hz, the noise sources are mostly local and transient. This is also the reason that the diurnal variation of the measured PSD level spans over two orders of magnitude, an order of magnitude higher than that observed for frequencies below 5 Hz.
The ambient noise data recorded by the array of 96 sensors were used for beamforming.
-600 -400 -200 0 200 400 600 800 Distance along longitude (m) 0 200 400 600 800 1000 1200 1400
Distance along latitude (m)
(a) sensor #173 10-1 100 101 Frequency (Hz) 10-16 10-14 10-12 10-10 PSD (m 2/s 4/Hz) (b) Jan 5 Jan 6 Jan 7 Jan 8 Jan 9 Jan 10 Jan 11 Jan 12 Jan 13 Jan 14 Jan 15 Jan 16 Jan 17 Jan 18 Jan 19 Jan 20 Peterson's LNM/HNM
Figure 2.7: (a) Sensor 173 marked with a red dot alongside the sensor layout in cartesian coordinates. (b) Daily average of the PSD of the seismic ground velocity as measured by sensor 173 in the period January 5 to 20, 2016.
Beamforming was performed at first on every hour of data. Hence for a set of 96 sensors and with a window length of 10 minutes and an overlap of 50% between consecutive windows, the size of the data matrix Xmo(f )following Eq. (2.9) was 96 × 11. The beampower for
every hour of data was computed following Eq. (2.17). Trial array response vectors were used in the slowness range 2.5 × 10−4to 5 × 10−3s/m at an interval of 8.33 × 10−6s/m, and
in the azimuth range 0 to 359.5◦at an interval of 0.5◦. Hence the total number of trial array
response vectors wLwas 571×721. Beampowers computed hourly as a function of slowness
and azimuth are then averaged for all hours of the day. As was the case for Fig. 2.6 we plot the beampower as a function of slowness (u) and azimuth (θ) on a polar graph, where the
(a) Hourly Averaged PSD for node 173 (0.1-20.0 Hz)
6 8 10 12 14 16 18 20 22
Days in month of January, 2016 2 4 6 8 10 12 14 16 18 20 Frequency (Hz) -180 -160 -140 -120 -100 -80
(b) Hourly Averaged PSD for node 173 (0.1-2.0 Hz)
6 7 8 9 10 11 12 0.2 0.4 0.6 0.8 1 Diurnal variation in seismic activity 20log 10(m/s/ Hz) Peak microseismic activity
Days in month of January, 2016
5 6 7 8 9 10 11 12 3 6 9 0 12 Frequ ency (Hz) Peri od (s) (c) North sea M SW (m) 3 6 9 0 12 Local earthquake
Figure 2.8: (a) Spectrogram of the seismic ground velocity measured by Sensor 173 for the entire measurement period between January 5 to 22, 2016 with a temporal resolution of one hour and in the frequency band 0.1 to 20 Hz. Typical diurnal variations due to man-made seismic noise is observed. An earthquake of moment magnitude 1.7 epicentered at Zuidbroeke measured on January 7, 2016 is also shown. (b) Zoomed spectrogram in the frequency band 0.1to 1 Hz showing the temporal variation of the secondary microseism energy. (c) Observed MSW of the North Sea (blue curve) between January 5 and 11, 2016 and the dominant wave period (red curve) measured during the same time.
slowness increases radially outwards and the azimuth is measured clockwise from the North. Fig. 2.9 shows the beampower plotted as a function of slowness-azimuth at frequencies of 0.3, 0.5, 0.7, 1.2, 1.4, and 1.7 Hz. In the frequency band 0.2 to 1.2 Hz, beamforming shows the dominant noise direction to be North-West which explains their oceanic origin. However, it is only the secondary microseismic peak that we observe due to the limited sensitivity of the 4.5 Hz geophones. At frequencies above 1.2 Hz, noise originating South of the array ap-pears as semi-circular arc. This corresponds to the noise due to a nearby road that borders the Southern edge of the array. At 1.6 Hz, an alias of the noise source is observed and is marked in Fig. 2.9. The true azimuth and velocity of propagation of the surface waves at 1.6 Hz are 227◦(measured clockwise from North) and 384 m/s. The alias is observed at an azimuth
of 32◦ and a velocity of 277 m/s. An easy way of removing such an aliased surface wave
arrival is to interpret the observations using the theoretical array response estimates (for this case, if we limit the minimum velocity to 300 m/s instead of 200 m/s which is the value used for computing the beampowers as shown in Fig. 2.9, we can remove the alias from the new u − θdomain). However, at high frequencies with the presence of many such aliases in u − θ domain, it is practically impossible to select a specific slowness range without any prior idea of the expected dispersion of the surface waves. Hence, extraction of the dispersion curve for this array was limited to a maximum of 1.7 Hz.
2.3.2
Extended Spectral Autocorrelation
The ESAC method was first introduced by Ling & Okada in 1993 [78], following the work of Aki in 1957 [73]. Since it is based on the Spectral Auto-correlation (SAC) method, we
Figure 2.9:Beamforming output for January 8, 2016 with slowness 2.5 × 10−4to 5 × 10−3 s/m on the radial axis and azimuth in the range 0 to 359.5◦ at frequencies of 0.3, 0.5, 0.7,
1.2, 1.4, and 1.7 Hz.
first discuss the SAC method briefly before describing the ESAC method. The normalized cross-spectra of the seismic data between the mthand nthsensors can be expressed as
Smn(f ) = 1 M PM i=1 iSmn(f )) q 1 M2( PM i=1 iSmm(f ))(P M i=1 iSnn(f )) , (2.16)
where f represents the frequency,iSmn(f )the cross-spectrum for the ithdata segment, and
M the total number of data segments. We consider an array layout with N seismic sensors, such that we have one sensor at the center which coincides with some chosen origin and the remaining sensors are distributed equally spaced in azimuth along the circumference of a circle. In practice, sensors are distributed along the circumference of several such circles of increasing radii. Then the cross-spectrum between the mthsensor with location coordinate
(xm, ym)and the nthsensor with location coordinate (xn, yn)can be rewritten in the form
Smn(f, rmn), where rmn = {(xm− xn)2+ (ym− yn)2}1/2. Considering m = 0, which
corresponds to the sensor at the center of the array and r0n = r, we can then write the
azimuthally averaged cross-correlation function as S(f, r) = 1 N − 1 N −1 X n=1 < S0n(f, r). (2.17)
Following the work of Aki in 1957 [73] the phase velocity c(f) of the propagating surface wave can be obtained from
S(f, r) = J0(
2πf r
c(f )), (2.18)
where J0 is the zero-th order Bessel function. In order to estimate the phase velocity from
velocity is a function of frequency, keeping the frequency constant and then trying to satisfy Eq. (2.18) is deemed suitable. The azimuthal averaging of the cross-spectrum data performed in Eq. (2.17) needs the sensors to be spaced equally in azimuth along the circumference of the circle of radius r, and it is one of the limitations of the SAC method. However, the ESAC method remedies this constraint in the array design by using a slightly modified version of Eq. (2.18). The phase velocity c(f) is found to satisfy the following equation in a least square sense, S0n(f, r0n) = J0 2πf r0n c(f ) , (2.19)
where n = {1, 2, 3, ..N − 1} are the number of sensors surrounding one central sensor positioned at the origin of the array layout. The phase velocity c(f) is estimated by fixing the frequency at a particular value and then fitting the zero-th order Bessel function, which is now only a function of the radial distance r0n. The best fitting Bessel function is obtained by
minimizing the error function E(f ) = N −1 X n=1 h S0n(f, r0n) − J0 2πf r0n c(f ) i2 . (2.20)
The fundamental difference between the SAC method and the ESAC method is in the way the cross-spectrum of the seismic data is used. While the SAC method uses an azimuthal averaging of the computed cross-spectrum between one central sensor and the sensors placed on the circumference of a circle, the ESAC method relies on simply using the temporally averaged cross-spectrum between a central sensor and the rest of the sensors, irrespective of their position in the array layout. In reality, the temporal averaging of the cross-spectrum of the seismic data over long time periods is a substitute to positioning the sensors on a cir-cumference of a circle. Over long periods of ambient noise recording, the assumption that noise sources originate from all directions and can be treated as a stationary random field is the key behind the working of the ESAC method. In other words, the azimuthal averaging performed in the SAC method by positioning sensors at equally spaced azimuth and at a fixed radial distance from the central sensor is equivalent to using just one sensor pair and instead positioning the noise sources equally spaced in azimuth on a circle around the sensor pair. As a result, the SAC method always works even if the noise source distribution around the array layout is not azimuthally isotropic. It is also the reason why the ESAC method might not work if the temporal averaging of the cross-spectrum is over short time spans and the noise source distribution around the sensor pair is not isotropic.
We present here an example of the SAC method by using the passive seismic data recorded from the Groningen array. Although the positioning of the sensors was not along the cir-cumference of circles with increasing radii, but rather on a regular square grid, it was still possible to implement the SAC method. The first step is the selection of a central sensor. The array of 96 sensors deployed in Groningen extends upto a maximum of 1000 m along two orthogonal cartesian directions. Hence, we select sensor 10 approximately located in the center of the array and as shown in Fig. 2.10(a) as the central sensor. This is followed by binning the distance of each sensor in the array from the selected central sensor. A histogram of these inter-sensor distances is computed for a bin size of 10 m, and then the distance bins that comprise greater than 4 sensors are selected for the analysis. A total of 11 such distance
bins are identified and their respective distances from the central sensor are shown with ar-rows marked in the histogram plot in Fig. 2.10(b). For clarity, the circles corresponding to each of these selected distance bins, are also shown in Fig. 2.10(a) with the sensors and the respective circle passing through the group of sensors plotted with a distinct color. The radius of each circle is computed as the mean of the distance of the sensors in that group from the central sensor. Now, since the sensors are partitioned into groups with a bin size of 10 m, all the sensors classified to a particular circle might not exactly lie on the circumference. We also note that for inter-sensor distances above 500 m, the azimuthal coverage of the sensor positions in the last three rings is skewed.
After the selection of the central sensor and classifying the sensor locations into different
cir-0 200 400 600 800 1000 1200
Distance along longitude (m) -800 -600 -400 -200 0 200 400 600
Distance along latitude (m)
100 200 300 400 500 600 Inter-sensor distance (m) 0 2 4 6 8 10 12
Sensor pair count
10
(a) (b)
Figure 2.10:(a) Layout of the sensors used for the SAC method with sensor locations that lie on a circle of fixed radius (error: ±10 m) shown with a distinct color and the central sensor marked as 10. (b) Histogram of the distances of the all the sensors from the central sensor with a distance bin size of 10 m and the arrows with different colors point to each collection of sensors that lie on a circle of fixed radius (shown in (a)).
cles, we now select the data window length and the percentage of overlap for computing the cross-spectrum following Eq. (2.16). We select a window length of 10 minutes and an over-lap of 50% between consecutive windows. First we compute the cross-spectrum between the central sensor and all the sensors that belong to one particular ring as shown in Fig. 2.10(a). The temporal averaging of the cross-spectrum between the central sensor and another sensor on a given circle is carried out for a day of seismic data. Hence, for a window length of 10 minutes and with an overlap of 50% we have a total of 314 cross-spectrum segments for a day of data. After the temporal averaging of the cross-spectrum has been done, we perform a spa-tial averaging by using the temporally averaged cross-spectrum from all the sensors that were classified to a particular ring. The spatially averaged cross-spectra which are just a function of frequency are then assigned to the circle in consideration. This process is then repeated for all the rings, such that at each frequency we have cross-spectrum magnitudes as function of radius of the rings. Fig. 2.11 shows the cross-spectrum magnitudes obtained using Eq (2.17) as a function of radius at six different frequencies of 0.4, 0.6, 0.8, 1.0, 1.2, and 1.6 Hz. The black curve which corresponds to the best fitting zero-th order Bessel function as expressed
in Eq. (2.18) is fitted to the observed cross-correlations by minimizing the error function in Eq. (2.20) by applying the Marquadt least-square method [79]. Repeating this minimization scheme at the desired set of frequencies gives us the phase velocity c(f) as a function of frequency. In Fig. 2.11 we also plot the Bessel functions corresponding to a 10% error in the estimated phase velocity. Most of the observed cross-correlation magnitudes shown as the green dots in Fig. 2.11 lie within the blue and the red dotted curves corresponding to +10% and −10% error in the estimated phase velocity respectively.
0 100 200 300 400 500 600 -1 0 1 Normalized cross-correlation -10 100 200 300 400 500 600 0 1 Fitted J0 J0, 90% vest J0, 110% vest Observed 0 100 200 300 400 500 600 -1 0 1 Normalized cross-correlation 0 100 200 300 400 500 600 -1 0 1 0 100 200 300 400 500 600 Radial distance (m) -1 0 1 Normalized cross-correlation 0 100 200 300 400 500 600 Radial distance (m) -1 0 1 0.4 Hz 0.6 Hz 0.8 Hz 1.2 Hz 1.0 Hz 1.6 Hz
Figure 2.11:Spatially averaged real value of the normalized cross-correlation expressed as function of the radius of the rings shown with the green dots and the best fitting zero-th order Bessel function shown by the black curve at frequencies of 0.4, 0.6, 0.8, 1.0, 1.2, and 1.6 Hz. The dotted red and blue curves show the Bessel functions corresponding to a ±10% error in the estimated phase velocity.
2.4
Rayleigh wave phase velocity estimation
Up to this point we have discussed the two methodologies that can be used for estimating the phase velocity and the direction of propagation of Rayleigh waves. While the ESAC method directly gives us an estimate of the phase velocity c(f) upon solving the minimization prob-lem stated in Eq. (2.20), extracting the Rayleigh wave phase velocity from beamforming is not obvious. Following Eq. (2.15) we know that the beampower is a function of slow-ness magnitude u, azimuth of propagation θ and frequency f. However, in order to obtain the phase velocity as a function of frequency, we need to express the beampower only as a function of u and f. In order to accomplish this in the u − θ plane, we sum the beampower across all values of azimuth for a particular slowness and frequency. Hence, at one particu-lar frequency if the beampower is expressed as a matrix with slowness increasing down the column and azimuth increasing along the row, we add up the matrix along each row to obtain a single vector that comprises the beampower only as a function of slowness. If this process is repeated for all frequencies of interest we again get a matrix, but now instead of azimuth varying along a row, the frequency varies. Following this, we pick the slowness value along
each column of the matrix corresponding to the maximum beampower. Repeating this pro-cess for all columns of the matrix returns the slowness magnitudes for every frequency, which is equivalent to the observed phase velocity of the Rayleigh waves.
Fig. 2.12(a) shows a surface plot of the beampower as a function of phase velocity (1/u) and
Figure 2.12: (a) Azimuthally averaged beampower computed for a day of ambient noise recording on January 12, 2016 plotted as a function of velocity and frequency, with the black curve showing the observed Rayleigh wave dispersion. (b) Observed (red dots) and mean Rayleigh wave phase velocity (solid blue curve) with the respective standard deviations (blue vertical lines) obtained from beamforming analysis of ambient noise recorded between Jan-uary 5 and 22, 2016 at an interval of 0.1 Hz in the frequency band 0.4 to 1.7 Hz. (c) Es-timated Rayleigh wave phase velocity (red dots) plotted as a function of frequency obtained from ESAC analysis of every hour of ambient noise data recorded on January 12, 2016. (d) Observed (red dots) and mean Rayleigh wave phase velocity (solid blue curve) with the re-spective standard deviations (blue vertical lines) obtained from ESAC analysis of ambient noise recorded between January 5 and 22, 2016 at an interval of 0.1 Hz in the frequency band 0.4 to 1.7 Hz.
frequency. The plot shown in Fig. 2.12(a) corresponds to the average beampower for a day of seismic noise recording. The black curve in Fig. 2.12(a) joins the slowness values corre-sponding to the maximum beampower at each frequency. Since we estimate the beampower as a function of u and θ at an interval of 0.1 Hz, we get a total of 14 values of velocity in the frequency band 0.4 to 1.7 Hz. Fig. 2.12(b) shows a scatter plot of the estimated dispersion curves obtained from beamforming analysis of all 18 days of ambient noise recordings. The mean dispersion curve is shown with the blue solid curve and the associated standard devia-tion of the observed phase velocity at each frequency shown with the error bars. Due to the limited maximum aperture of the array, we observe that the error associated with the phase velocity estimation at frequencies below 0.6 Hz is considerably higher as compared to the high frequencies. Fig. 2.12(c) shows a scatter plot of the phase velocities as a function of fre-quency obtained from ESAC analysis of a day of ambient noise recordings. Phase velocities are again estimated at 14 discrete frequencies in the band 0.4 to 1.7 Hz at an interval of 0.1 Hz. Fig. 2.12(d) shows the scatter of the observed phase velocity dispersion obtained from the ESAC analysis of all 18 days of ambient noise recording. The solid blue line shows the mean
phase velocity from all the observations and the standard deviation is shown using the error bars. Comparing the mean phase velocities obtained from beamforming and ESAC analysis, a good agreement is observed for frequencies greater than 0.8 Hz. In the frequency band 0.4 to 0.8 Hz the phase velocities obtained from ESAC show a few jumps especially between 0.6 and 0.7 Hz. This artefact observed in the ESAC results can be ascribed to two factors. Firstly, the noise source in the frequency band is microseismic in origin and has a fixed azimuth of 330◦measured clockwise from the North, and the very assumption that long period ambient noise recordings tend to depict an isotropic noise source distribution does not hold. Secondly, due to inadequate azimuthal coverage of the sensors in the rings with higher inter-sensor sep-aration (greater than 500 m), the spatial averaging of cross-correlation magnitudes inherent with the SAC method does not overcome the directional bias in the observations.
2.5
Rayleigh wave phase velocity inversion strategies
After the surface wave phase velocities have been estimated by using ambient seismic noise recordings, the next step is to compute a subsurface velocity model that accurately mimics the observed surface wave dispersion. In geophysical terms, this type of an analysis is referred to as an inverse problem. The parameters to be estimated in such a problem are the P-wave, S-wave velocities and the density model of the subsurface. These parameters are generally estimated for each subsurface layer, and hence the number of layers in the subsurface model is also unknown. Inversion of dispersion curves to estimate deep crustal S-wave velocity models was first tried by Dorman & Ewing in 1962 [80]. Song et al. in 1989 [81] presented two examples of surface wave phase velocity measurements and studied the sensitivity of the observed dispersion with respect to changes in the model parameters. It was observed that the sensitivity of the P-wave velocities and density to the measured dispersion is significantly lower than the sensitivity to S-wave velocity model. Hence a reliable estimation of only the S-wave velocity model is possible. Turner in 1990 [82] presented the feasibility of obtaining reliable near surface P-wave velocity models from dispersion measurements of both Love and Rayleigh waves. In the process of inverting for the S-wave velocity model there are two important aspects to be considered: the accuracy and the resolution of the model. The maximum wavelength upto which dispersion measurements are available limits the accuracy of the model and the distribution of the dispersion data over the wavelength band of interest controls the resolution of the estimated S-wave velocity model.
For the inversion, the first assumption is that the earth can be represented as a layered sub-surface model such that the observed Rayleigh wave phase velocity is a non-linear implicit function of the subsurface parameters and can be expressed as
F (fj, cj, Vs, Vp, ρ, h) = 0, (2.21)
where fjis the frequency in Hz, cjis the observed Rayleigh wave phase velocity at frequency
fj; Vs= (Vs1, Vs2, ..., Vsn)T is the S-wave velocity vector, with Vsjthe shear wave velocity
of the jthlayer; n is the number of layers; V
p = (Vp1, Vp2, ..., Vpn)T is the P-wave velocity
vector, with Vpj the P-wave velocity of the jth layer; ρ = (ρ1, ρ2, ..., ρn)is the density
vector, with ρj the density of the jthlayer; h = (h1, h2, ..., hn)is the layer thickness vector
the function F with respect to the model parameters Vs, Vp, ρ, hin order to understand the
accuracy of the inversion scheme. Hermann, 1987 [83] implemented a damped least-square method that uses the analytical formulation of the Jacobian matrix and executes the inversion starting with an initial model. At each iteration, the model updates on the basis of a simpler linear problem and a misfit is computed. The updates to the subsurface model are performed based on the misfit that is computed at each iteration. However, when the uncertainties in the observed dispersion curve are higher, derivative-based methods fail. It always gives a solution that is close to the starting model and the solution is stuck at one of the several local minima of the misfit function. As a result for these types of problems in geophysics, inversion schemes that make use of pseudo-random sampling of the parameter space (Monte-Carlo type) are suitable. Some of the popular schemes are simulated annealing proposed by Sen & Stoffa in 1991 [84] and genetic algorithms proposed by Lomax & Snieder in 1994 [85]. One of the more recent ones which we use in this thesis is the Neighborhood algorithm proposed by Sambridge in 1999 [86]. Broadly, the objective of these algorithms is to investigate the entire parameter space, looking for the set of parameters that fit the observed dispersion curve best.
2.5.1
Theoretical dispersion curve computation
The first step in the process of estimating an 1D S-wave velocity model is to be able to compute a Rayleigh wave dispersion curve theoretically, given a 1D subsurface model. This is also commonly referred to as the forward problem. Since the neighborhood algorithm searches through tens of thousands of 1D subsurface models before it can estimate the best fitting ground model, the computation of the theoretical dispersion curve needs to be fast and accurate. The method used for solving the forward problem is the same as described in Section 1.2.5 of Chapter 1. In short, we use the Dunkin’s formulation [45] for propagating the motion-stress tensor from the bottom of the layered half-space to the surface, followed by a fast root search method based on the bisection scheme and Lagrange’s polynomial inter-polation [87]. The results are the f-k pairs that satisfy the stress free boundary condition on the surface and the radiation condition at the bottom of the layered half-space. The forward computation of the dispersion curve is exclusively for Rayleigh waves, as we are interested in processing only the vertical component of the recorded ground motion. Moreover, the forward computation of the dispersion curve is also able to identify the dispersion curves per-taining to the different propagation modes of the Rayleigh waves. In the surface wave case studies discussed in the thesis, we are interested only in the fundamental mode and the first overtone of the Rayleigh waves. As an example, Fig. 2.13(c) shows the theoretical phase ve-locities corresponding to the fundamental mode and the first overtone of the Rayleigh waves for the five layer subsurface structure shown in Fig. 2.13(a) and (b).
2.5.2
Model parametrization
In the neighborhood algorithm discussed previously, the parametrization of the model space is an important step. The important factors to be decided by the user are the number of layers to use for the inversion and also the bounds on the values of Vp, Vsand ρ for each layer.
The number of parameters for the inversion problem should not be too large such that the problem becomes ill-posed. On the contrary, the bounds on the parameter space should not
0 1000 2000 3000 4000 Velocity (m/s) 0 50 100 150 200 250 300 350 400 Depth (m) (a) P-wave S-wave 1.5 2 2.5 3 Density (gm/cc) 0 50 100 150 200 250 300 350 400 Depth (m) (b) 0 5 10 15 20 Frequency (Hz) 0 500 1000 1500 (c) Fundamental mode First overtone
Figure 2.13: (a) An example P and S-wave velocity-depth model. (b) An example density model of the subsurface. (c) Theoretical Rayleigh wave phase velocities for the fundamental mode and the first overtone corresponding to the subsurface model shown in (a) and (b). be too strict such that the solution is stuck at a local minima due to insufficient model-search space. The best strategy is to start with a simple model and then progressively add new layers to the starting model and redo the inversion until a satisfactory fit is observed [88].
The first factor that the user must decide on is the maximum depth to which a reliable sub-surface model can be estimated. This decision is dependent on the frequency band of the observed dispersion curve. Typically, the maximum depth upto which a 1D ground model can be reliably estimated equals about one-third the maximum wavelength of the surface wave [89]. The next factor that the user must decide is the number of layer interfaces that needs to be modeled and their depth ranges. If no a-priori information is available about the geology of the region, the best way to proceed is by trial and error and executing the inver-sion several times in order to get an idea of the depth ranges within which the theoretically computed dispersion curve is the most sensitive. In such a case starting with a layer over half-space and then gradually increasing the number of layers in the model is the best way to proceed. From a theoretical standpoint the dispersion curve is the least sensitive to the subsurface density model [90]. Hence for most applications, the values of density are kept fixed for each layer based on the information available from previous geotechnical studies. The values of Vpand Vsare linked to each other using a bound on the Poisson ratio which
should lie between −0.5 and 0.5. For selecting the search ranges on the values on Vpand Vs
for each layer, we generally start with bigger search space and then gradually diminish the search range based on the observation from several inversion runs. It is also possible to fix the ratio of Vp/Vsfor each layer if the material properties are following the work of Castagna
in 1985 [91].
For the studies discussed in this thesis we have a-priori geological information about the region based on previous seismic studies. In the context of the dispersion measurements from the Groningen array, we parameterize the problem into a four layer model including the half-space. As we already mentioned, the dispersive properties of Rayleigh wave are mostly sensitive to the S-wave velocity structure. As a consequence, we fixed the values of the P-wave velocities and density for each layer. The only unknowns for the problem are the
