Method for cross-correlating receiver signals

Abstract

A computer-implemented method of cross-correlating signals detected by receivers is disclosed. The method comprises: receiving first and second signals respectively detected by first and second receivers arranged on a surface; converting the first signal and the second signal from a time domain to a frequency domain; determining a downsampling factor; downsampling the converted first signal and the converted second signal based on the downsampling factor in the frequency domain; and cross-correlating the downsampled first and second signal to obtain a frequency domain Green's function. Unlocking insights from Geo-Data, the present invention further relates to improvements in sustainability and environmental developments: together we create a safe and liveable world.

Description

METHOD FOR CROSS-CORRELATING RECEIVER SIGNALSFIELD

[0001] The disclosure relates to methods and systems for analysing a target region beneath a surface of the earth. More particularly, the disclosure relates to a method and system for determining one or more ground properties of the sub-surface target region based on ambient noise measured at or near the surface. In particular, the disclosed methods and systems relate to cross-correlation of ambient noise signals measured by receivers. Unlocking insights from Geo-Data, the present invention further relates to improvements in sustainability and environmental developments: together we create a safe and liveable world.BACKGROUND

[0002] There is a general and ongoing need for systems and methods for determining sub-surface ground parameters. In particular, there is a need for systems and methods that can be used to model the properties of a target volume beneath the surface of the earth to provide information useful for infrastructure planning. Determination of sub-surface ground properties during the early planning phase of construction projects reduces uncertainty during the location determination, foundation design, and construction phases of a project. This in turn reduces delays, overspend, and unnecessary use of material resources (e.g. concrete) during construction.

[0003] One key parameter for the determination of ground characteristics in a volume of interest is shear-modulus and shear-velocity Vs. The shear velocity Vs is the velocity at which a shear wave moves through the material and is controlled by the shear modulus of the material. The relationship between shear-velocity and shear modulus G is defined by Vs = ^G / p, where p is the density of the material. Measurement of Vs therefore provides a valuable insight to the material properties of a subsurface ground region.

[0004] Spectral analysis of surface waves (SASW) and multi-channel analysis of surface waves (MASW) are both examples of techniques for gathering surface wave information that can be used in the determination of ground properties in a sub-surface volume. In both of these techniques, surfacelevel vibrations resulting from an active source (e.g. a weight drop) are measured from and the dispersion of the resulting surface waves is studied. ReMi (Refraction Microtremor) is another surfacelevel technique that uses ambient noise and surface waves to infer ground properties of a sub-surface region based on the observation of ambient noise at the surface.

[0005] Down-hole and cross-hole techniques can also be used to determine ground properties of a sub-surface region. In both of these methods, a receiver located in a bore hole measures waves received from an active source located elsewhere. In a down-hole technique, one of the source and the receiver is located at a sub-surface location within the bore hole and the other of the source and the receiver is located at the surface. In a cross-hole technique, a source is located in a first bore hole, with a receiver located in a second bore hole. In both down-hole techniques, the propagation and dispersionof the received waves are studied to infer the properties of the material through which the waves from the receiver have travelled.

[0006] Invasive techniques for measuring material ground properties of a sub-surface region can often present logistical challenges such as long duration of the processes and thus low cost efficiency, e.g. long process set-up, long acquisition times and / or heavy machinery, equipment and processes. Invasive techniques are often particularly undesirable, especially in urban or inaccessible environments, and are often prohibitively expensive. Invasive techniques may also be unfriendly to the environment, e.g. cause disturbance to the local fauna. Conversely, current surface-level techniques may lack the accuracy and reliability of more invasive analysis techniques.

[0007] Additionally, there is a need to reduce the amount of processing required in the determination of sub-surface ground properties as known techniques require significant computing overhead.SUMMARY

[0008] According to a first aspect of the present disclosure, there is provided a computer-implemented method of cross-correlating signals detected by receivers. The receivers, which may also be considered to be sensors, may be geophones, accelerometers, seismometers, vibration sensors and / or transducers. The method comprises: receiving first and second signals respectively detected by first and second receivers arranged on a surface; converting the first signal and the second signal from a time domain to a frequency domain; determining a downsampling factor; downsampling the converted first signal and the converted second signal based on the downsampling factor in the frequency domain; and cross-correlating the downsampled first and second signal to obtain a frequency domain Green’s function.

[0009] This method makes computation of the Green’s functions for receiver pairs more computationally efficient and reduces memory usage. Converting into the frequency domain means that the cross-correlation can be computed much faster. In addition, downsampling the converted first signal and the second signal reduces amount of data to be cross-correlated by the downsampling factor, meaning the cross-correlation of the first signal and the second signal requires less computational overhead and memory usage.

[0010] In some implementations, the step of determining a downsampling factor comprises the steps of: estimating a minimum velocity of surface waves, Vmin, (optionally a minimum velocity of Rayleigh waves) that travel through a medium beneath the surface; determining a minimum frequency,of the received first signal and / or the received second signal; determining a maximum distance, / ), between receivers of an array comprising the first and second receivers; determining a maximum timelag, Tmax, based on the maximum distance, D, the minimum velocity of surface waves, Vmin, and the minimum frequency, obtaining a window time length, T, of the received first signal and / or the received second signal; determining the downsampling factor, M, based on the window time length, T, and the maximum timelag, Tmax.

[0011] The downsampling factor is determined in such a way that a reduction in computational overhead and memory usage is achieved while avoiding aliasing of the signals being processed.

[0012] In some implementations, determining a maximum timelag, Tmaxcomprises determining a maximum timelag, Tmaxaccording to:D kT1max > — - T7 1 ' - f -min Jl whereD is the maximum distance between receivers of an array comprising the first and second receivers;Vmlnis the estimated minimum velocity of surface waves that travel through a medium beneath the surface; k is a coefficient accounting for cross-correlated signal duration; / i is the minimum frequency of the received first signal and / or the received second signal.

[0013] The Nyquist-Shannon sampling theorem indicates that a periodic time signal containing no frequencies higher than fmaxcan be reconstructed using a so called Nyquist-rate sampling dt = — - —2 *fmax without any loss of information. Making use of time-frequency duality, we can write the dual version of the Nyquist-Shannon theorem. It indicates that, in the frequency domain, a periodic signal which inverse Fourier Transform signal has a null amplitude after Tmax(in time), can be reconstructed using a sampling Af = TmaX(in the frequency domain) without any loss of information.

[0014] Practically, as the cross-correlated signal goes from -Tmaxto +Tmax, its length is 2 * Tmax, so we end up with Af = - as the Nyquist-rate sampling. The key then is not to overpass2 Tmax^max this limit in final sampling. As such, the maximum down-sampling factor for the signals in the method requires the determination of Tmax. The down-sampling factor may be applied without any concerns of aliasing, if we stay within those limits for sampling.

[0015] In some implementations, determining the downsampling factor, M, comprises determining the downsampling factor, M, according to:TM < -4 T ‘max whereT is the window time length of the received first signal and / or the received second signal; and Tmaxis the maximum timelag.

[0016] The downsampling factor is determined in such a way that a reduction in computational overhead and memory usage is achieved while avoiding aliasing of the signals being processed. In particular, downsampling the signal reduces the data size of the signal when in the frequency domain without any loss by the determination of downsampling factor M. As will be discussed further herein, downsampling the signal reduces the number of data points for the received signal. As such, any computational processing carried out on the signal, e.g. cross-correlation, may be accomplished with reduced processing and in less time. Moreover, computational processing and memory requirementsassociated with accessing the signal data by a processor from memory, e.g. storing, retrieving or loading, are reduced with smaller data size from a downsampled signal when in the frequency domain. Further, handling any intermediate results which may be calculated after a processing step or all processing steps requires similarly less processing and memory when the signal is downsampled. Therefore, a reduction in hardware requirements and improved efficiency for processing signals, and in particular for the cross correlation of signals, are accomplished.

[0017] In some implementations, a fast Fourier transform is used to convert the first signal and the second signal from a time domain to a frequency domain. By using a fast Fourier transform the complexity of the transformation is reduced from O(n2) to O(n. login)).

[0018] In some implementations, the method may further comprise the step of converting the frequency domain Green’s function from a frequency domain to a time domain.

[0019] In some implementations, an inverse fast Fourier transform is used to convert the frequency domain Green’s function from a frequency domain to a time domain. By using an inverse fast Fourier transform the complexity of the transformation is reduced from O(n2) to O(n. login .

[0020] In some implementations, the method may further comprise the step of outputting the time domain Green’s function.

[0021] In some implementations, the method may further comprise the step of outputting the frequency domain Green’s function.

[0022] In some implementations, the first and second receivers are part of an array of receivers arranged on the surface.

[0023] In some implementations, the method is repeated for each pair of receivers in the array.

[0024] In some implementations, the receivers comprise geophones.

[0025] In some implementations, each signal in the plurality of signals is indicative of ambient noise measured at or near the surface by a respective receiver.

[0026] In some implementations, the method may further comprise determining a model of ground properties of a sub-surface target region based on the Green’s functions. The ground property may comprise shear wave velocity, Vs, or shear modulus. The model may be a 3D model. In one example, the method comprises determining one or more velocity dispersion profiles as a function of sensed surface wave frequency, based on the determined cross-corelated signals; generating a 3- dimensional phase velocity dispersion profile from the first and second phase velocity dispersion profiles; determining a plurality of group velocity dispersion profiles as a function of sensed surface wave frequency; generating a respective shear wave velocity model as a function of depth from the 3 dimensional phase velocity dispersion profile and each group velocity dispersion profile; and generating a 3 dimensional shear wave velocity model as a function of depth from the respective shear wave velocity models. The shear wave velocity model may be used in a workflow process for determining the location or design requirements of one or more construction projects.

[0027] According to another aspect of the present disclosure, there is provided a system comprising one or more processors and one or more memories having stored thereon computer-readable instructions configured to cause the one or more processors to perform any of the methods disclosed herein.

[0028] According to another aspect of the present disclosure, there is provided a computer-readable medium comprising instructions, that, when executed by one or more data processing apparatus, cause the one or more data processing apparatus to perform any of the methods disclosed herein.

[0029] According to another aspect of the present disclosure, there is provided a computer program comprising instructions which, when the program is executed by a computer, cause the computer to perform any of the methods disclosed herein.BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Disclosed implementations will now be described by way of example to illustrate aspects of the disclosure and with reference to the accompanying drawings, in which:Figure 1 shows a cross-sectional view of a sub-surface target region with a plurality of receivers arranged at the surface;Figure 2 shows a plot identifying the ray paths between a plurality of receivers;Figure 3 shows two example signals respectively detected by two receivers of the grid array shown in Figure 1 ;Figure 4 shows a method for computing the cross-correlations;Figure 5 shows the example signals of Figure 3 converted into the frequency domain;Figure 6 shows downsampled frequency domain signals of Figure 5;Figure 7 shows the downsampled frequency domain signals of Figure 6 and the resultant crosscorrelation of those signals in the form of a frequency domain Green’s function;Figure 8 shows the time domain Green’s function of Figure 7 once it has been converted back into the time domain;Figure 9A shows an example method for determining the downsampling factor M;Figure 9B shows an example determination of Tmaxin accordance with the method shown in Figure 9A.Figure 10 shows a 3D model of seismic shear-wave velocity (Vs) for a sub-surface target region; andFigure 11 shows a block diagram of a computing device which can be used to implement the disclosed methods.DETAILED DESCRIPTION

[0031] This detailed description describes, with reference to Figures 1 and 2, an approach to measuring structural properties of a ground volume using geophones. Next, with reference to Figures 3-9B, novel methods for performing cross-correlation of signals detected by geophones are disclosed. An example output that may be generated from the obtained Green’s functions is shown in Figure 10. Finally, a computing device that may be used to perform the disclosed methods is described with reference to Figure 11 .

[0032] The following examples will be described in the context of a geophone array, to aid understanding. It will, however, be appreciated that the disclosed systems and methods are applicable to a variety of receiver types, including but not limited to geophones, accelerometers, seismometers, vibration sensors and / or transducers. The disclosed methods may be applied to cross-correlation of any suitable set of signals.

[0033] The methods and systems disclosed herein relate generally to cross-correlation of signals detected by geophones on a surface. In one particular example, the geophones are placed on a ground surface and the detected signals are ambient noise signals. Cross-correlation of these signals provides useful insight into the structure of the surface and sub-surface volume on which the geophones are placed, as described in more detail below. Due to the potentially very large number of signals being cross-correlated, however, methods are needed to ensure that cross-correlation of these signals is computationally feasible. The disclosed methods and systems provide such a mechanism for performing the required cross-correlations in a computationally efficient and practically feasible manner.

[0034] Before turning to the details of the disclosed cross-correlation methodology, some background relating to determination of surface and sub-surface properties using geophones will first be provided.

[0035] The present disclosure describes systems and method for determining the ground properties of a sub-surface volume. Although, it will be appreciated that the methods herein may be applied to model physical properties other than shear velocity, such as compressional wave velocity, density, elastic modulus, or, if a viscoelastic model is being used, optionally also viscosity coefficients Qsand QP. In general, the model may define multiple physical property values, the focus of the following detailed description will be the determination of shear velocity, and the related quantity, shear modulus.

[0036] Shear modulus is a measure of the elastic shear stiffness of a material and represents the deformation of a solid when it experiences a force parallel to one of its surfaces while its opposite face experiences an opposing force. Such forces and their effects in sub-surface ground volumes, are an important parameter for study before and during the design of building and infrastructure projects. To determine the shear modulus of a volume, the shear velocity, Vs, is determined. This in turn gives an indication of the stiffness of the sub-surface material, and its ability to support structures extending above and / or through the volume.

[0037] In the context of ground study, two types of waves are generally distinguished: P-waves, in which particles in the volume oscillate in the direction of movement of the wave, cause a compression and de-compression of the ground as the waves propagate through the ground. S-waves are shear waves, in which particles oscillate in a direction perpendicular to the direction of propagation of the waves.

[0038] P-waves and S-waves are body waves and propagate in all directions through the body of the volume. The interaction of P- and S-waves with the earth’s surface generates surface waves, which propagate along that surface. Several types of surface waves can be distinguished. In the systems and methods described herein, Rayleigh waves are measured and studied because it is convenient to measure the vertical component of surface vibrations. However, it will be appreciated that othersurface waves (e.g. Love waves) may be measured and harnessed in the systems and methods described herein.

[0039] Because surface waves propagate in 2D (at the surface), they attenuate less rapidly than body waves (which propagate in 3D). Surface waves are generally present within a depth range of one wavelength from the surface, generally travel more slowly and have a predominantly lower frequency than body waves. This lower attenuation, slower travel time and lower frequency of surface waves makes their study particularly attractive for the purposes of determining shear velocity, Vs. Since the surface waves have lower attenuation, the signal strength is better maintained through over a longer travel distance. The resulting measurement results therefore generally have a higher signal quality (signal-to-noise ratio) than body wave studies.

[0040] The methods and systems of the present disclosure harness the study of surface waves to provide insight into the ground properties of a sub-surface target volume. In general terms, the present disclosure provides methods for analysing one or more ground properties, such as the shearwave velocity, Vs, of a sub-surface target region. The method involves receiving a data set indicative of ambient noise at a surface above a sub-surface target region, analysing the background wavefield to study the dispersive behaviour of the surface waves measured at the surface, and determining the ground properties of the sub-surface target region based on the study of the dispersive behaviour of the observed waves.

[0041] Referring now to Fig. 1 , a cross-sectional view of a sub-surface volume 100 is shown. A surface 102 extends above the sub-surface volume. Surface waves propagate along the surface 102, as shown at point A, where a schematic representation of a particle oscillation (due to Rayleigh wave propagation) at the surface above the target sub-surface volume is shown. As illustrated, the oscillation of the particle P is partly vertical, partly in the direction of propagation and movement is therefore substantially ellipsoid.

[0042] At a surface 102 above the volume 100, a plurality of receivers 104a, 104b is arranged in a 2D array. The receivers (collectively referred to as 104) may be geophones configured to measure the vertical component of the surface waves propagating across the surface 102. However, it will be appreciated that the propagation of surface waves may also be measured with alternative sensing means, for example: accelerometers, seismometers, vibration sensors, and / or transducers.

[0043] The receivers 104 are configured to measure ambient noise. That is, the wavefield present due to background noise (as opposed to noise from an active source such as a hammer drop or explosion). The background noise may be natural (e.g. due to lapping ocean waves, wind, and other naturally occurring vibrations) or cultural (e.g. due to human activity, including traffic, machinery, etc.).

[0044] As shown in Figure 1 , the receivers 104 are arranged in a grid array at the surface, the grid extending in two directions. It should be noted that the surface above the target region may, in many cases, not be planar. The array of receivers 104 may therefore not be truly “2-dimensional” (each receiver may be offset from its neighbours in the grid in the z-direction, as well as the x- and y- directions). However, such a grid arrangement of receivers will be referred to as a 2D array herein. It will also be appreciated that in the context of sub-surface ground study, the term 2D array may also be used to denote an array configured to capture 2D information (e.g. a line of receivers configured toprovide information regarding a 2D slice of a target region) and a 3D array may refer to a grid of receivers configured to capture 3D information. However, in the context of the present information, the term 2D in the context of receivers is intended to refer to the 2-dimensional configuration of the receivers rather than the information gathered.

[0045] Moreover, although the receivers 104 are shown at the surface in Figure 1 , the receivers 104 may also be disposed near the surface (e.g. partially buried, or immediately beneath the surface). In the context of the present application, ‘at the surface’ will be understood to mean at or near the surface, such that the receivers measure surface waves.

[0046] In the schematic shown in Figure 1 , the receivers 104 are arranged in a regular grid, with the inter-receiver distance equal across surface 102. However, in some implementations, the receivers 104 are arranged with varying density across the surface 102.

[0047] The array of receivers 104 shown in Figure 1 can be used to gather a data set indicative of ambient noise at a surface above a sub-surface target region for which a model of shear velocity Vs is desired that is the subject of the present application.

[0048] To determine the shear velocity from the observation of surface waves (in particular Rayleigh waves), the dispersive behaviour of the surface waves is studied. Surface waves are dispersive, i.e. their velocity is dependent on frequency. Since seismic velocities increase with depth in the earth, normal surface wave dispersion shows a decrease of surface-wave velocity with increasing frequency. It is by studying the behaviour of surface waves at a surface above a volume that the ground properties of the volume can be determined.

[0049] There are two ways to measure the velocity of dispersive surface waves and a distinction is made between the determination of group velocity or phase velocity.

[0050] The group velocity of a wave is the velocity with which the overall envelope shape of the wave's amplitudes - known as the modulation or envelope of the wave - propagates through space. The group velocity is equivalent to the speed with which the energy of the wave propagates through the volume and is measured by determining the wave propagation between two points. The group velocity is obtained as a time-of-flight (that is, travel time) measurement between a (virtual) source and a receiver.

[0051] The phase velocity is the velocity at which the phase of any one frequency component of the wave travels. As such, the phase velocity is expressed as a function of frequency. To measure the phase velocity, at least two measurement nodes (e.g. receivers) are chosen to measure the waves propagating through the volume to determine relative time-of-flight between the nodes for different frequencies. The result is the phase velocity as a function of frequency as averaged over the volume between the two measurement nodes. Phase velocity is acquired as a point in 2D phase space (dispersion spectrum) which is itself obtained by a 2D transform (such as slant-stack, Radon, FK, or the like) of an array of recorded waveforms (time-distance space). The wavelength of the surface wave is indicative of its propagation depth. As a result, the phase velocity as a function of frequency is indicative of the shear velocity Vs as a function of depth.

[0052] The receivers 104 at surface 102 provide the measurement nodes for the study of the surface wave behaviour as described above. However, the receivers 104 described above are configured torecord passive noise. Therefore, the measurement nodes of Figure 1 do not represent a real point source and associated receiver. However, receivers 104 can act as a virtual source receiver pair, as explained below.

[0053] Receivers 104a and 104b form a first receiver pair in the array at surface 102. Neither receiver 104a nor 104b represents a point source for noise recorded at the other of receiver 104a, 104b. However, by cross-correlating the received signal at receiver 104a and 104b, receivers 104a and 104b can act as a (virtual) source receiver pair, where each receiver of the pair records a signal as though the signal had originated at the other of the pair. The cross-correlation is performed using the principle of interferometry.

[0054] The cross-correlation of passive noise measured at respective pairs of receivers at the surface shown in Figure 1 can be used to reproduce a response from the sub-surface target volume, as if it were induced by an impulse point source, which is equal to Green’s function.

[0055] In other words, a response that is received by cross-correlating two receiver recordings can be interpreted as a response that would have been measured at one of the receiver locations as if there were a source at the other. Various approaches to determining the Green’s function for a virtual sourcereceiver pair are known, with an overview of the approaches described in “Tutorial on Seismic Interferometry: Part 1 - Basic Principles and Applications”; GEOPHYSICS. Vol. 75, No 5 (Sept-Oct 2010; P.75A195075A209; Wapenaar et al.).

[0056] In Figure 1 , only one pair of receivers is labelled (104a, 104b). However, it will be appreciated that for each receiver 104 in the array, every other receiver in the array may act as the other half of a source receiver pair. In this manner, the Green’s function for each source receiver pair may be obtained. The Green’s functions across the plurality of virtual source receiver pairs is studied to determine the dispersive behaviour of the surface waves.

[0057] Figure 2 shows a plurality of virtual source-receiver pairs across a surface above a sub-surface region of interest. Ray paths 206 between source-receiver pairs are indicated. Note that the ray paths here indicate the propagation of a wave (modelled from) a virtual source at a receiver in a centre of the plot to each of the other receivers in the array. The receivers are not labelled individually but are represented with an inverted triangle in the plot. Each receiver can similarly act as a point source. As can be seen from Figure 2, the array of receivers allows source receiver pairs defining ray paths that extend in different directions, and with varying inter-receiver spacing. The background shading and contour rings indicate the travel-time field from the central receiver to the other receivers 104. There are also corresponding ray paths between each receiver and all other receivers, i.e. between every pair of receivers. These ray paths are not depicted in Figure 2 for simplicity.

[0058] As described above, cross-correlating signals obtained by receivers (e.g. geophones 104 in Figure 1) and obtaining resulting Green’s functions enables useful sub-surface investigation to be performed. In practice, however, there may be up to several thousand receivers 104 provided in each array. To obtain a complete data set, the signal measured at each receiver 104 must be cross-correlated with the signal measured at each other receiver 104 in the array, as well as with itself. Each measured signal is typically recorded over a time period of a day or so, at a sampling frequency that is usually between 100Hz and 500Hz. As a result, the signal for each receiver, which may in some contexts bereferred to as a “data trace”, may comprise on the order of 20 million samples, when taking a typical sampling frequency of 250Hz. It will be appreciated that, as the number of receivers increases, crosscorrelation of the signals, each containing perhaps 20 million samples, rapidly becomes very computationally expensive. In particular, a prohibitive amount of memory is required to cross-correlate the entire array. Merely as an example, taking an array of 1000 receivers would result in 10002cross correlations, with each receiver signal comprising around 2 x 107samples as noted above. This results in a computation involving around 2 x 1013data points. As a result of this very large number of data points, cross-correlation of a large array of receiver signals may be practically infeasible within an acceptable timeframe for all but the most powerful computing systems. The approach disclosed herein addresses this problem by providing methods which make the cross-correlation more efficient and computationally feasible for more commonplace computing hardware.

[0059] The disclosed approach will now be described with reference to Figures 3-9B. Figure 3 depicts two example signals stand s}- respectively detected by two receivers 104 of the grid array shown in Figure 1. Typically, each signal is indicative of the ambient noise measured at or near the surface by the respective receiver which recorded the signal. Hence, signal stis representative of the ambient noise measured by receiver i and signal s}- is representative of the ambient noise measured by receiver j. These are examples of two signals that can be cross-correlated to obtain a Green’s function for enabling useful sub-surface investigation to be performed.

[0060] The present inventors have identified a method for cross-correlating signals from two receivers 104 in a faster, more computationally efficient manner, and in which the amount of memory required is reduced. This method is depicted in Figure 4, and will now be explained.

[0061] A method 400 for computing the cross-correlations will now be described with reference to Figure 4. The method 400 of Figure 4 may be performed by one or more computing devices, as will be explained in more detail below.

[0062] The method 400 of Figure 4 begins, at step 401 , by receiving first and second time domain signals stand s}- respectively detected by first and second receivers arranged on a surface, for example two of the receivers 104 of Figure 1. In this example, both signals are ambient noise signals.

[0063] At step 403, the first signal stand the second signal s}- are converted from a time domain to a frequency domain. In an example implementation, a fast Fourier transform may be used to convert the first signal stand the second signal s}- from a time domain to a frequency domain to enable faster computation of the cross-correlation. In other words, St= FFTfSf) and Sj = FFT(sj), where Stis the first signal in the frequency domain and Sj second signal in the frequency domain.

[0064] By using a fast Fourier transform algorithm, the complexity of the transformation is reduced from O( / V2) to O(N. log(N ) (as compared to the straightforward Fourier Transform). Further, by processing the signal in the frequency domain, the cross-correlation may be completed in N operations, as opposed to N2in the time domain (for N samples signals). This appears when using Equations 1 and 2 below: in the time domain, Equation 1 needs N operations for each timelag t value, whereas, in the frequency domain, only N operations are needed to cross-correlate the two signals (i.e. one operation per frequency value). A cross-correlation in the time domain is as follows:where cl}is the cross-correlation of a first signal stobtained by a first receiver with a second signal s7obtained by a second receiver, t is the timelag for which the cross-correlation function is computed, and T is the integration variable (time).

[0065] A cross-correlation in the frequency domain can be found using the equation:where C,7is the Fourier transform of ctj, Stis the Fourier transform of shand Sj is the Fourier transform of s . The operator * represents the complex conjugate. is the angular frequency.

[0066] Figure 5 depicts the example signals Stand Sj of Figure 3 once they have been converted into the frequency domain, at step 403.

[0067] Returning to Figure 4, at step 405, a downsampling factor M is determined. Determination of the downsampling factor M is described in greater detail with reference to Figure 9A below. Importantly, the downsampling factor M is determined in such a way that loss of information is avoided.

[0068] At step 407, the first converted frequency domain signal Stand the second converted frequency domain second signal Sj are downsampled in the frequency domain based on the downsampling factorM. This means that the amount of data to be processed is globally reduced by a factorM, reducing the memory and processing requirements of the cross-correlation process.

[0069] In an example implementation, the first converted frequency domain signal stand the second converted frequency domain signal sj can be downsampled as follows:St= St(l; M + 1; 2M + 1; 3M + 1; kM + 1; ...; ) , S}= 5,(1; M + 1; 2M + 1; 3M + 1; kM + 1; ... ; N), where k e Wand max{kM + 1} < N (3) where Stis the downsampled first converted frequency domain signal and S}downsampled second converted frequency domain signal. N is the total number of datapoints in the converted frequency domain first signal Stand the converted frequency domain second signal Sj .

[0070] In other words, for each of the frequency domain signals Stand Sj, one data point for every M data points of the signal is sampled. It will be apparent that once the frequency domain signals Stand Sj have been downsampled, storage and processing requirements for operations using those signals, such as cross-correlating the signals, are reduced.

[0071] Figure 6 depicts the downsampled frequency domain signals Stand S}of Figure 5 obtained through step 407 of Figure 4, once they have been downsampled based on the downsampling factorM.

[0072] At step 409, the downsampled first and second signals Stand S}are cross-correlated to obtain a frequency domain Green’s function as follows:where Gl}is the Green’s function obtained from cross-correlation of the downsampled first and second signals Stand S}. The operator * represents the complex conjugate. is the angular frequency.

[0073] As discussed, performing cross-correlation in the frequency domain reduces the number of operations from O( / V2) in the time-domain to O(N). Using the downsampling factor M, the complexity is further reduced by a factor of M.

[0074] Figure 7 depicts the downsampled frequency domain signals Stand S}of Figure 6 and the resultant cross-correlation of those signals in the form of a frequency domain Green’s function

[0075] After step 409, the optional step of converting the frequency domain Green’s function G^ to the time domain can be carried out to obtain a time domain Green’s function G *.7. . An inverse fast Fourier transform may be used for this purpose. In other words, GBy using an inverse fast Fourier transform the complexity of the transformation is reduced from 0(nN2) to 0(nN. log(nN)).

[0076] Figure 8 depicts the time domain Green’s function G of Figure 7 once it has been converted back into the time domain from the frequency domain.

[0077] Turning now to Figure 9A, as mentioned above this depicts an example method 900 for determining the downsampling factor M (step 405 of the method 400 shown in Figure 4). The method 900 of Figure 9A may be performed by one or more computing devices, for example the same one or more computing devices carrying out the method 400 of Figure 4, as will be explained in more detail below. In order to determine the downsampling factor M, we first need to determine the time lag Tmax.

[0078] The Nyquist-Shannon sampling theorem indicates that a periodic time signal containing no frequencies higher than fmaxcan be reconstructed using a so called Nyquist-rate sampling dt = — - — 2 *fmax without any loss of information. Making use of time-frequency duality, we can write the dual version of the Nyquist-Shannon theorem. It indicates that, in the frequency domain, a periodic signal which inverse Fourier Transform signal has a null amplitude after Tmax(in time), can be reconstructed using a sampling Af = (in the frequency domain) without any loss of information. Tmax

[0079] Practically, as the cross-correlated signal goes from -Tmaxto +Tmax, its length is 2 * Tmax, so we end up with Af = - as the Nyquist-rate sampling. The key then is not to overpass2 Tmax ^max this limit in final sampling. As such, the maximum down-sampling factor for the signals in the method requires the determination of Tmax. The down-sampling factor may be applied without any concerns of aliasing, if we stay within those limits for sampling.

[0080] The downsampling factor M may be determined by first determining the maximum time lag Tmaxwhich is the maximum time for wave signal propagation and obtaining T which is the window time length of the received signal. Using the dual version of the Shannon-Nyquist sampling theorem, we can ensure that no aliasing of the signal will occur if the signal is sampled in the frequency domain at a frequency rate down to Af = - . That is, as the original signal is sampled at a 8 Jf = - rp frequency J ^‘ max ^‘ max1rate, the frequency can be downsampled by a factor M which is less than or equal to °f TSo the down samplingafactor can be chosen as an integaer such that M < — Sf = U -max , which sets the upper limit for which down sampling factor may be applied to ensure no aliasing in the data will occur.

[0081] The method 900 of Figure 9A begins, at step 901 , by estimating a minimum velocity Vmlnof surface waves, such as Rayleigh waves, travelling through a medium beneath the surface on which the first and second receivers (i.e. the receivers providing the first and second signals stand sy) are arranged. This estimate can be based on known or assumed parameters of the medium beneath the surface. This estimate can be checked at a later point by ensuring that the effective signal obtained after cross-correlation is not truncated at its Tmaxvalue, which would mean the frequency downsampling has created aliasing.

[0082] At step 903, a minimum frequencyof the received first signal and / or the received second signal is determined. This can be obtained by analysing the first and second signals stand sy .

[0083] At step 905, the maximum distance / ) between receivers in the array to which the first and second receivers, providing the first and second signals stand sy , belong is determined.

[0084] At step 907, a maximum timelag Tmaxis determined based on the distance / ) the minimum velocity of Rayleigh wavelengths Vmlnand the minimum frequency f .

[0085] In an example implementation, determining the maximum timelag Tmaxcomprises determining a maximum timelag Tmaxaccording to:where D is the aforementioned maximum distance, Vminis the aforementioned estimated minimum velocity, f is the aforementioned minimum frequency and fc is a coefficient accounting for crosscorrelated signal duration, fc can be 1 or, up to 3 or more if a greater safety buffer is desired.

[0086] Determining maximum timelag Tmaxin line with Equation 5 ensures that no aliasing in time will occur.

[0087] At step 909, a window time length T of the received first signal and / or the received second signal stand sy is obtained. The window time length is the length of the signal being processed, which may be up to several days as noted above. In some examples, it may be a 1 hour window or any suitable window length that allow for stacking of one or more windows together after processing.

[0088] At step 911 , the downsampling factor M is determined based on the window time length T and the maximum timelag Tmax.

[0089] In an example implementation, determining the downsampling factor M comprises determining the downsampling factor M according to:where T is the aforementioned window time length and Tmaxis the aforementioned maximum timelag.

[0090] Figure 9B depicts an example determination of Tmaxin accordance with the method 900 shown in Figure 9A. Figure 9B depicts a number of signals from a number of receivers 104 stacked on top ofeach other and labelled sequentially on the y-axis and propagating through time on the x-axis, in accordance with examples described above. The minimum velocity of surface waves, Vmin, in this example is 150 m / s. The value for the coefficient accounting for cross-correlated signal duration k is 2. The maximum distance, D, between receivers of an array is 1000 m. The minimum frequency,of the received first signal and / or the received second signal is 1 Hz. Applying Equation 5 to these values, a Tmax°f 8.67 s is determined. As can be seen in Figure 9B, choosing Tmaxin this manner ensures no signal data is lost when downsampling as all signals fall within the Tmaxwindow. It can be seen that choosing a larger coefficient accounting for cross-correlated signal duration k increases the safety buffer provided.

[0091] The Green’s functions that may be obtained through the above method may be used in a variety of ways to investigate the structure and properties of the surface and sub-surface volume on which the receivers are placed. The details of these implementations are beyond the scope of the present disclosure and will be apparent to a skilled person familiar with the use of Green’s functions in the context of seismology.

[0092] In at least one implementation, following acquisition of data indicative of ambient noise, and the analysis of said ambient noise data to model a plurality of response signals at a plurality of virtual source-receiver pairs, the received data may now be used to determine the ground properties of the sub-surface target volume by solution of an inverse problem, in which the response at the receivers (e.g. receivers 104) is known but the ground properties of the sub-surface target volume is not (yet) known.

[0093] The starting model for the solution of the inverse problem may be determined in multiple ways since it sets initial physical property values, for example for Vs, within the sub-surface target volume, to be refined by solution of the inverse problem. The starting model can comprise historical data based on known local or regional data, predicted values based on selective invasive test methods, coarse surface level studies or a theoretical model based on local or regional historical knowledge. Although the starting model is to be improved upon and need not be a particularly accurate representation of the sub-surface target region, the more accurately the starting model reflects the physical properties of the sub-surface target region, the more accurate the final model of the ground properties of the target region will be. Therefore, improvements in the starting model can provide an improved output for the methods and systems described herein. Once a starting model has been selected, the inverse problem can be solved via inversion, as explained below in connection with the example model shown in Figure 10.

[0094] Figure 10 shows a 3D shear-velocity Vs model for a sub-surface target region of interest. The region of interest was “imaged” to a depth of around 100m through acquisition of ambient noise data using receivers. X, Y and Z co-ordinates in meters are shown. The recorded noise data signals were cross-correlated in the manner described above to efficiently generate a set of Green’s functions. Tomographic inversion was then performed based on the Green’s functions to obtain the 3D model shown. In particular, the inter-receiver travel-times of the surface wave modes were obtained to derive group velocity dispersion curves as a function of frequency and location. From this, a 3D model of seismic shear-wave velocity (Vs) was obtained, which is shown in Figure 10. Different velocities are represented in different colours (darker areas generally representative of higher velocities) andtransitions between regions of different shear wave velocity are clearly visible. These factors indicate the varying composition and structure of portions of the subsurface target region, providing valuable insight into the constitution of the target region, for example in terms of ground stiffness distribution and soil classification, and thus its suitability as a site for possible construction projects.

[0095] The above description has provided a variety of examples to illustrate the disclosed methods. However, the described arrangements and methods are merely exemplary, and it will be appreciated by a person skilled in the art that various modifications can be made without departing from the scope of the appended claims. In particular all example values for variables are merely exemplary. The use of geophones in the examples and Figures is merely illustrative and other types of receiver can be used to implement the disclosed techniques.

[0096] While various specific combinations of components and method steps have been described, these are merely examples. Components and method steps may be combined or ordered in any suitable arrangement or combination. Components and method steps may also be omitted to leave any suitable combination of components or method steps.

[0097] Figure 11 shows a block diagram of one implementation of a computing device 1100 within which a set of instructions, for causing the computing device to perform any one or more of the methodologies discussed herein, may be executed. In alternative implementations, the computing device may be connected (e.g., networked) to other machines in a Local Area Network (LAN), an intranet, an extranet, or the Internet. The computing device may operate in the capacity of a server or a client machine in a client-server network environment, or as a peer machine in a peer-to-peer (or distributed) network environment. The computing device may be a personal computer (PC), a tablet computer, a set-top box (STB), a Personal Digital Assistant (PDA), a cellular telephone, a web appliance, a server, a network router, switch or bridge, or any machine capable of executing a set of instructions (sequential or otherwise) that specify actions to be taken by that machine.

[0098] Further, while only a single computing device is illustrated, the term “computing device” shall also be taken to include any collection of machines (e.g., computers) that individually or jointly execute a set (or multiple sets) of instructions to perform any one or more of the methodologies discussed herein. More particularly, a number of computing devices can be used to compute cross-correlations of signal data subsets independently and in parallel, as described above. Each computing device may have the structure shown in Figure 11 . Alternatively, a plurality of processors within a single computing device, such as computing device 1100, can perform the independent computations.

[0099] The example computing device 1100 includes a processor 1102, a main memory 1104 (e.g., read-only memory (ROM), flash memory, dynamic random access memory (DRAM) such as synchronous DRAM (SDRAM) or Rambus DRAM (RDRAM), etc.), a static memory 1106 (e.g., flash memory, static random access memory (SRAM), etc.), and a secondary memory (e.g., a data storage device 11 18), which communicate with each other via a bus 1 130.

[0100] Processor 1102 represents one or more general-purpose processors such as a microprocessor, central processing unit, or the like. More particularly, the processor 1102 may be a complex instruction set computing (CISC) microprocessor, reduced instruction set computing (RISC) microprocessor, very long instruction word (VLIW) microprocessor, processor implementing otherinstruction sets, or processors implementing a combination of instruction sets. Processor 1102 may also be one or more special-purpose processors such as an application specific integrated circuit (ASIC), a field programmable gate array (FPGA), a digital signal processor (DSP), network processor, or the like. Processor 1102 is configured to execute the processing logic (instructions 1122) for performing the operations and steps discussed herein.

[0101] The computing device 1100 may further include a network interface device 1108. The computing device 1 100 also may include a video display unit 1110 (e.g., a liquid crystal display (LCD) or a cathode ray tube (CRT)), an alphanumeric input device 11 12 (e.g., a keyboard or touchscreen), a cursor control device 1114 (e.g., a mouse or touchscreen), and an audio device 1 116 (e.g., a speaker).

[0102] It will be apparent that some features of computer device 1100 shown in Figure 11 may be absent. For example, one or more computing devices 1100 may have no need for display device 1 110 (or any associated adapters). This may be the case, for example, for particular server-side computer apparatuses 1100 which are used only for their processing capabilities and do not need to display information to users. Similarly, user input device 1112 may not be required. In its simplest form, computing device 1100 comprises processor 1102 and memory 1104.

[0103] The data storage device 1118 may include one or more machine-readable storage media (or more specifically one or more non-transitory computer-readable storage media) 1128 on which is stored one or more sets of instructions 1122 embodying any one or more of the methodologies or functions described herein. The instructions 1122 may also reside, completely or at least partially, within the main memory 1104 and / or within the processor 1102 during execution thereof by the computer system 1100, the main memory 1104 and the processor 1102 also constituting computer-readable storage media.

[0104] The various methods described above may be implemented by a computer program. The computer program may include computer code arranged to instruct a computer to perform the functions of one or more of the various methods described above. The computer program and / or the code for performing such methods may be provided to an apparatus, such as a computer, on one or more computer readable media or, more generally, a computer program product. The computer readable media may be transitory or non-transitory. The one or more computer readable media could be, for example, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, or a propagation medium for data transmission, for example for downloading the code over the Internet. Alternatively, the one or more computer readable media could take the form of one or more physical computer readable media such as semiconductor or solid state memory, magnetic tape, a removable computer diskette, a random access memory (RAM), a read-only memory (ROM), a rigid magnetic disc, and an optical disk, such as a CD-ROM, CD-R / W or DVD.

[0105] In an implementation, the modules, components and other features described herein can be implemented as discrete components or integrated in the functionality of hardware components such as ASICS, FPGAs, DSPs or similar devices.

[0106] A “hardware component” is a tangible (e.g., non-transitory) physical component (e.g., a set of one or more processors) capable of performing certain operations and may be configured or arranged in a certain physical manner. A hardware component may include dedicated circuitry or logic that is permanently configured to perform certain operations. A hardware component may be or include aspecial-purpose processor, such as a field programmable gate array (FPGA) or an ASIC. A hardware component may also include programmable logic or circuitry that is temporarily configured by software to perform certain operations.

[0107] Accordingly, the phrase “hardware component” should be understood to encompass a tangible entity that may be physically constructed, permanently configured (e.g., hardwired), or temporarily configured (e.g., programmed) to operate in a certain manner orto perform certain operations described herein.

[0108] In addition, the modules and components can be implemented as firmware or functional circuitry within hardware devices. Further, the modules and components can be implemented in any combination of hardware devices and software components, or only in software (e.g., code stored or otherwise embodied in a machine-readable medium or in a transmission medium).

[0109] Unless specifically stated otherwise, as apparent from the following discussion, it is appreciated that throughout the description, discussions utilizing terms such as "receiving”, “determining”, “identifying”, “estimating”, “obtaining”, “converting”, “downsampling”, “cross-correlating” or the like, refer to the actions and processes of a computer system, or similar electronic computing device, that manipulates and transforms data represented as physical (electronic) quantities within the computer system's registers and memories into other data similarly represented as physical quantities within the computer system memories or registers or other such information storage, transmission or display devices.

[0110] It is to be understood that the above description is intended to be illustrative, and not restrictive. Many other implementations will be apparent to those of skill in the art upon reading and understanding the above description. Although the present disclosure has been described with reference to specific example implementations, it will be recognized that the disclosure is not limited to the implementations described, but can be practiced with modification and alteration within the spirit and scope of the appended claims. Accordingly, the specification and drawings are to be regarded in an illustrative sense rather than a restrictive sense. The scope of the disclosure should, therefore, be determined with reference to the appended claims, along with the full scope of equivalents to which such claims are entitled.

Claims

CLAIMS1 . A computer-implemented method of cross-correlating signals detected by receivers, the method comprising: receiving first and second signals respectively detected by first and second receivers arranged on a surface; converting the first signal and the second signal from a time domain to a frequency domain; determining a downsampling factor; downsampling the converted first signal and the converted second signal based on the downsampling factor in the frequency domain; and cross-correlating the downsampled first and second signal to obtain a frequency domain Green’s function.

2. The computer-implemented method of claim 1 , wherein the step of determining a downsampling factor comprises the steps of: estimating a minimum velocity of surface waves, Vmin, that travel through a medium beneath the surface; determining a minimum frequency, f , of the received first signal and / or the received second signal; determining a maximum distance, D , between receivers of an array comprising the first and second receivers; determining a maximum timelag, Tmax, based on the maximum distance, D, the minimum velocity of surface waves, Vmin, and the minimum frequency, f obtaining a window time length, T, of the received first signal and / or the received second signal; determining the downsampling factor, M, based on the window time length, T, and the maximum timelag, Tmax.

3. The computer-implemented method of claim 2, wherein determining a maximum timelag, Tmaxcomprises determining a maximum timelag, Tmaxaccording to:D kT1max > - 11— f''min JI where4. The computer-implemented method of either of claim 2 or claim 3, wherein determining the downsampling factor, M, comprises determining the downsampling factor, M, according to:TM < — 4 T -1max whereT is the window time length of the received first signal and / or the received second signal; andTmaxis the maximum timelag.

5. The computer-implemented method of any preceding claim, wherein a fast Fourier transform is used to convert the first signal and the second signal from a time domain to a frequency domain.

6. The computer-implemented method of any preceding claim, further comprising the step of converting the frequency domain Green’s function from a frequency domain to a time domain.

7. The computer-implemented method of claim 6, wherein an inverse fast Fourier transform is used to convert the frequency domain Green’s function from a frequency domain to a time domain.

8. The computer-implemented method of either of claim 6 and claim 7, further comprising the step of outputting the time domain Green’s function.

9. The computer-implemented method of any preceding claim, further comprising the step of outputting the frequency domain Green’s function.

10. The computer-implemented method of any preceding claim, wherein the first and second receivers are part of an array of receivers arranged on the surface.11 . The computer-implemented method of claim 10, wherein the method of any preceding claim is repeated for each pair of receivers in the array.

12. The method of any preceding claim, wherein the receivers comprise geophones, optionally configured to measure a vertical component of vibrations at the surface.

13. A system comprising: one or more processors; one or more memories having stored thereon computer readable instructions configured to cause the one or more processors to perform operations comprising the steps of any of claims 1-12.

14. A computer readable medium comprising instructions, that, when executed by one or more data processing apparatus, cause the one or more data processing apparatus to perform operations comprising the steps of any of claims 1-12.

15. A computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the method of any of claims 1 -12.

Images