Tomographic inversion by means of a gaussian bayesian prior with a location dependent standard deviation

Abstract

There is provided a method for determining a physical property of a subsurface target volume, as well as systems and computer readable media for performing the same. 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. The method comprises performing tomographic inversion on obtained empirical travel time data. Performing tomographic inversion comprises obtaining a posterior probability density function corresponding to each cell of the first plurality of cells, wherein obtaining the posterior probability density function comprises obtaining a prior velocity model and defining a prior probability density function indicative of velocity for each of the cells of the first plurality of cells, for each receiver pair of the subset, determining a modelled travel time using the initial velocity model values associated with the two or more cells that are traversed by the surface wave path and determining a posterior probability density function of the velocity model based on the modelled travel times and the prior probability density function wherein the prior probability density function is dependent on a distance of each of the cells from the first point on the surface.

Description

Tomographic inversion by means of a Gaussian Bayesian prior with a location dependent standard deviationTECHNICAL FIELD

[0001] Unlocking insights from Geo-Data, the present invention relates to improvements in sustainability and environmental developments: together we create a safe and liveable world. More particularly, the disclosure relates to methods and systems for analysing a target region beneath a surface of the earth. In particular, the disclosure provides methods and systems for providing a model of a subsurface target region with improved accuracy.BACKGROUND

[0002] There is a general and ongoing need for systems and methods for determining subsurface 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 the information useful for infrastructure planning and foundation design. Determination of subsurface ground properties during the early planning phase of construction projects reduces uncertainty during the location determination, design, and construction phases of a project. This in turn reduces delays, overspending, and unnecessary use of material resources (e.g., concrete) during construction and an asset’s lifecycle.

[0003] One key parameter for the determination of ground properties in a volume of interest is shearmodulus G and shear velocity Vs. The shear velocity 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 is defined by Vs = ^G / p, where p is the density of the material. Measurement of shear velocity 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 material properties in a subsurface volume. In both of these techniques, surfacelevel vibrations are measured, either from a passive source (vibrations in the surface as a result of ambient sources of noise) or an active source (e.g., a weight drop), and the dispersion of the recorded surface waves is recorded. ReMi (Refraction Microtremor) is another surface-wave technique surface measuring waves from ambient seismic noise recorded at the surface to infer material properties of a subsurface region.

[0005] Invasive down-hole and cross-hole techniques can also be used to determine material properties of a target subsurface region. In both of these methods, a receiver located in a bore hole measures waves recorded from an active source located elsewhere. In a down-hole technique, one of the source or the receiver is located at a subsurface location within the bore hole and the other of the source or 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 techniques, the propagation of therecorded waves are studied to infer the properties of the material through which the waves originating from the source have travelled.

[0006] Invasive techniques for measuring material properties of a subsurface region can often present logistical challenges, especially in urban or inaccessible environments, and are often prohibitively expensive to obtain the amount of data required to rely on these techniques alone. These techniques allow one-dimensional or something two-dimensional screening of the properties of the subsurface, require significant resources (equipment and professional staff), and are associated with safety risks and a negative environmental footprint. Conversely, current surface-level techniques may lack the accuracy and reliability of more invasive analysis techniques.OVERVIEW

[0007] According to a first aspect of the present disclosure, a method for determining a physical property of a subsurface target volume is provided. The method comprises receiving a plurality of signals detected by a plurality of receivers arranged on a surface above the subsurface target volume, wherein each respective signal of the plurality of signals is recorded by a respective receiver of the plurality of receivers, cross-correlating signals between the plurality of receivers to obtain empirical travel time data for surface waves between a plurality of receiver pairs, receiving at least one subsurface signal detected by at least one subsurface receiver at a first location within the subsurface target volume, wherein the first location defines a first point on the surface above the subsurface volume and determining a model of the physical property of the subsurface target volume, the model comprising a first plurality of cells, wherein each cell has a respective physical property value and wherein the plurality of cells spans at least the area of the surface. As used herein, “empirical” data or information refers to such data or information obtained empirically, i.e by means of real signals collected at real receivers positioned on the surface above the subsurface target volume. Determining the model comprises selecting a subset of the plurality of receiver pairs, wherein a surface wave path between each source and respective receiver traverses two or more cells of the first plurality of cells, and performing tomographic inversion on the empirical travel time data between each receiver pair in the subset to obtain the physical property value for each cell.

[0008] The receivers used in the present methods may be geophones, accelerometers, seismometers, vibration sensors and / or transducers. The receivers may collect data over a significant period of time. For example, the ambient seismic noise may be measured consecutively over a period of five days. This longer recording time results in the adequate retrieval of surface wave information from the ambient seismic noise recorded at or near the surface of the target region. In turn, and as described in more detail below, the (processed) surface wave information can be used in a tomographic inversion to obtain a shear wave velocity model of the subsurface.

[0009] Performing the tomographic inversion comprises obtaining a posterior probability density function corresponding to each cell of the first plurality of cells, wherein obtaining the posterior probability density function comprises obtaining an initial velocity model for each of the cells of the first plurality of cells, for each receiver pair of the subset, determining a modelled travel time using the initial velocitymodel values associated with the two or more cells that are traversed by the surface wave path and determining a posterior probability density function of the velocity model based on the modelled travel times and the prior probability density function, wherein the prior probability density function is dependent on a distance of each of the cells from the first point on the surface. By performing a tomographic inversion on the empirical travel time data for a plurality of receiver pairs in this way, a three-dimensional model of the physical property of the subsurface target volume can be obtained. This specific approach also takes into account the locational dependence between points, resulting in a smoother and more realistic model, i.e. being (relatively) free from physically impossible jumps in geophysical properties. This approach also allows for a higher resolution model to be produced without increasing the computational time to unfeasible levels. By using a prior probability density function with a distance dependent standard deviation for the tomographic inversion, the prior enforces a probabilistic preference for the invasively measured subsurface properties.

[0010] In some examples, the method for determining an in-situ physical property value of a subsurface target volume comprises recording at least one subsurface signal detected by at least one subsurface receiver at a second location within the subsurface target volume, wherein the second location is additionally defined by a second point on the surface above the subsurface volume, wherein the prior probability density function is dependent on a distance of each of the model cells from the second point on the surface. In some examples, the method further comprises receiving a plurality of subsurface signals corresponding to a plurality of locations, each location additionally defined by a point on the surface, wherein the prior probability density function is dependent on a distance of each of the cells from the plurality of points on the surface. This enforces a probabilistic preference for each of the in-situ physical property values at each of the measurement locations.

[0011] In some examples, at least one of the one or more subsurface receivers are inserted into the subsurface target volume. In some examples, at least one of the one or more subsurface receivers are inserted into a borehole. This enables in-situ physical property values obtained from any appropriate invasive technique to be incorporated into the prior probability density function.

[0012] In some examples, performing the tomographic inversion further comprises iteratively, perturbating an initial velocity model to form one or more perturbations of the initial velocity model, selecting a perturbation as an updated velocity model, determining a likelihood of the updated velocity model based on the modelled travel times and the obtained empirical travel time data, determining a posterior probability based on the likelihood of the updated velocity model and the prior probability density function, and calculating an acceptance probability based on a posterior probability of the updated velocity model and the posterior probability of the initial velocity model. In some examples, before determining a likelihood of the updated velocity model, performing the tomographic inversion comprises either for each receiver pair of the subset, determining a modelled travel time using the updated velocity model values associated with the two or more cells that are traversed by a surface wave path or determining updated velocity model values based on the previously computed surface wave paths from the initial velocity model. For example, the steps of perturbating the initial velocity, selecting an updated velocity model, determining a modelled travel time, determining a likelihood of the updated velocity and calculating an acceptance probability are iterated until a number of iterations meetsa first end condition. The updated velocity model of each iteration may be used as the initial velocity model of the next iteration. In some examples, determining the acceptance probability comprises accepting the updated model if the posterior probability of the updated model is greater than the posterior probability of the initial model. In some examples, determining the acceptance probability further comprises obtaining a ratio of the posterior probability of the updated velocity model to the posterior probability of the initial velocity model, generating a random number between 0 and 1 , accepting the updated model if the ratio is greater than the random number and reverting to the initial model if the ratio is smaller than the random number.

[0013] In some examples the steps of perturbating the initial velocity model, selecting an updated velocity model, determining a likelihood of the updated velocity model and calculating an acceptance probability are performed using the reversible-jump Markov chain Monte-Carlo algorithm and the cells are Voronoi cells and a random number of cells are initialised when obtaining the Initial velocity model. In some examples, the perturbations may be one or more of: a velocity change of a cell of the plurality of cells, a location change of a cell of the plurality of cells, a cell nucleus birth, a cell nucleus death, a change in a first noise parameter; or a change in a second noise parameter.

[0014] In some examples the prior probability density function is defined by:wherein the prior covariance matrix Cpis a diagonal matrix where the trace consists of a standard deviation of an ith model cell, C (i, i) = a2, matrix Nmcontains normalization factors of the probability density function of the ith cell, vector V is comprised of the seismic parameter values of all model cells n and vector V contains the measured seismic parameter values. This solves the problem of how to improve upon a uniform prior for inferring the velocity distribution describing the sub-surface velocity structure since the Gaussian nature of the probability density function imposes probabilistic preference for specific velocity values. As such, it is reflected that we’re more certain of specific values in specific regions of the model. In some examples, the standard deviation of the Gaussian prior probability density function of a cell increases up to a threshold L as the distance from the cell to the first point on the surface above the subsurface volume increases. This solves the problem of how to more accurately encapsulate the measurement point influence on the model as defining parameter L allows the model to capture that there is no influence of the first surface point after a certain distance.In some examples, the standard deviation of the prior probability density function comprises a scaling factor, wherein the scaling factor is indicative of the rate at which the influence of the seismic parameter value measured at the first point on the surface decreases with distance. This solves the problem of how to account for variations in the velocity structure of the subsurface volume as the use of a scaling factor allows the influence of the measurement at a measurement point to be altered - large scaling factors cause the standard deviation of the prior probability density function to rapidly decrease away from the measurement point such that the influence of the measured seismic parameter values on the velocity of the model is only at locations close to the point. Reducing the scaling factor increases the measurement point influence at distances further away from the measurement point. A very small scaling factor is equivalent to modelling a perfectly horizontal layered subsurface. In some examplesthe standard deviation of prior probability density function is defined as<7V= min aloc* exp (a L , where av. is the standard deviation of theGaussian probability density function describing the parameter value at location xt,y , alocthe standard deviation at the first point on the surface xloc,y and a represents a scaling factor which determines the rate with which the influence of the seismic parameter value measured at the first point on the surface decreases with distance as GVincreases until threshold value L. This solves the problem of how to reflect uncertainty of the velocity with respect to the measurement location since for smaller standard deviations the first probability distribution behaves like a Gaussian, centred around an in-situ velocity measurement viocand progressively further away from the location the standard deviation increases such that the distribution behaves uniformly between two fixed bounds wherein the model reflects that vioc has no influence on v. In some examples, the method comprises, outputting the updated model of the subsurface target volume to an output device.

[0015] According to another aspect of the present disclosure, there is provided a method for determining a physical property of a subsurface target volume. The method comprises receiving a plurality of signals detected by a plurality of receivers arranged on a surface above the subsurface target volume, wherein each respective signal of the plurality of signals is detected by a respective receiver of the plurality of receivers, cross-correlating signals between the plurality of receivers to obtain empirical travel time data for surface waves between a plurality of receiver pairs, and determining a model of the physical property of the subsurface target volume, the model comprising a first plurality of cells arranged in a first two-dimensional grid, wherein each cell has a respective physical property value and wherein the two-dimensional grid spans at least the area of the surface. Determining the model comprises selecting a subset of the plurality of receiver pairs, wherein a surface wave path between each source and respective receiver traverses two or more cells of the first plurality of cells, and performing tomographic inversion on the empirical travel time data between each receiver pair in the subset to obtain the physical property value for each cell. By performing a tomographic inversion on the empirical travel time data for a plurality of receiver pairs in this way, a three-dimensional model of the physical property of the subsurface target volume can be obtained. This is by means of obtaining a onedimensional model for each cell of the grid of cells, which collectively provide information on the physical property in three dimensions.

[0016] In some examples, performing the tomographic inversion comprises obtaining an initial velocity model comprising initial values indicative of velocity for each of the cells of the first plurality of cells, for each receiver pair of the subset, determining a modelled travel time using the initial first velocity model values associated with the two or more cells that are traversed by the surface wave path; determining a first error value indicative of the difference between the modelled travel times and the corresponding empirical travel times for each receiver pair in the subset; and determining an updated velocity model based on the error value, wherein the updated velocity model comprises updated values indicative of velocity for each of the first plurality of cells. In some examples, determining the modelled travel time, determining the first error value, and determining the updated velocity model based on the error valueare iterated until the error value meets a first end condition. The updated velocity model of each iteration is used as the initial velocity model of the next iteration. In some examples, determining the modelled travel time, determining the first error value, and determining the updated velocity model based on the error value are performed using iterative non-linear least-squares inversion, which is a computationally efficient inversion technique.

[0017] In some examples, performing the tomographic inversion further comprises, for each cell of the first plurality of cells, performing an inversion of the value indicative of velocity to obtain the physical property value for the cell. The inversion may comprise obtaining an initial model of the physical property of the subsurface volume, wherein the initial model comprises initial physical property values for each of the cells of the first plurality of cells; for each cell of the first plurality of cells, determining modelled values indicative of velocity based on the initial physical property value for the cell; determining a second error value indicative of the difference between the modelled values indicative of velocity and the corresponding values indicative of velocity obtained from the updated velocity model; and determining an updated model of the physical property of the subsurface volume based on the second error value, wherein the updated model of the physical property of the subsurface volume comprises updated physical property values for each of the first plurality of cells. Determining the modelled values indicative of velocity, determining the second error value, and determining the updated model of the physical property of the subsurface volume based on the second error value are iterated until the second error value meets a second end condition, and wherein the updated model of the physical property of the subsurface volume of each iteration is used as the initial of the physical property of the subsurface volume of the next iteration. Determining the modelled values indicative of velocity, determining the second error value, and determining the updated model of the physical property of the subsurface volume based on the second error value may be performed using a computationally efficient leastsquares inversion.

[0018] In alternative examples, performing tomographic inversion comprises: obtaining an initial model of the physical property of the subsurface volume wherein the initial model comprises initial physical property values for each of the cells of the first plurality of cells; determining an empirical dispersion function for each receiver pair of the subset using the empirical travel time data; determining a modelled dispersion function for each receiver pair of the subset using the initial physical property values associated with the two or more cells that are traversed by a surface wave travelling from the source to the receiver; determining a third error value indicative of the difference between the modelled dispersion function and the empirical dispersion function for each receiver pair; and determining an updated model of the physical property of the subsurface volume based on the third error value, wherein the updated model of the physical property of the subsurface volume comprises updated physical property values for each of the first plurality of cells.

[0019] In some examples determining the modelled dispersion function for each receiver pair may comprise averaging, the physical property values of the two or more cells that are traversed by a surface wave travelling from the source to the receiver, and calculating the modelled dispersion function using the averaged physical property value. Alternatively, determining the modelled dispersion function for each receiver pair may comprise calculating a cell dispersion function for each of the two or more cellsthat are traversed by a surface wave travelling from the source to the receiver to provide a plurality of cell dispersion functions, wherein calculating each cell dispersion function comprises using the physical property value of the respective cell, and averaging the plurality of cell dispersion functions. In some examples, the path taken by a surface wave between each source and receiver pair is represented by a straight or curved ray. In these examples, the averaging may comprise determining a weighted average according to a respective weight for each cell of the two or more cells, wherein the respective weight for each cell corresponds to a segment length of the ray path in that cell. In other examples, the path taken by a surface wave between each source and receiver pair is represented by a Fresnel zone. In these examples, the averaging comprises determining a weighted average according to a respective weight for each cell of the two or more cells, wherein the respective weight for each cell corresponds to a sensitivity value of the Fresnel zone in that cell.

[0020] In some examples, determining the updated model of the physical property of the subsurface comprises updating the initial physical property value of each cell according to a respective weighting of the path taken by a surface wave in that cell. For example, the change in physical property value of each cell may be weighted as described above.

[0021] In some examples, determining the second model comprises determining partial derivatives of the third error value with respect to the one or more physical properties of the first model; determining a sensitivity matrix comprising the partial derivatives; solving a linear system of equations defined by the sensitivity matrix to determine desired changes to the one or more initial physical property values of the initial model; and updating the initial model according to the desired changes to produce the updated model of the physical property of the subsurface volume.

[0022] In some examples, determining the modelled dispersion function based on the initial model, determining the third error value, and determining the updated model of the physical property of the subsurface volume are iterated until the third error value meets a third end condition, wherein the updated model of each iteration is used as the initial model of in the next iteration. The end condition may include one or more of; the third error value of a latest iteration is less than a predetermined threshold error value; and a difference between the third error value of a latest iteration and the third error value of a preceding iteration is less than a predetermined difference threshold.

[0023] In some examples, the determining the modelled dispersion function, the determining the third error value, and the determining the updated physical model are performed using non-linear iterative least squares inversion.

[0024] In some examples, the initial model of the physical property of the subsurface volume provides initial physical property values in at least two spatial dimensions, and optionally in three spatial dimensions. In some examples, the initial model of the physical property of the subsurface volume is derived from empirical phase dispersion data obtained from signals detected between a second subset of the plurality of receiver pairs.

[0025] In some examples the empirical dispersion function and the modelled dispersion function comprise group velocity dispersion functions and / or phase velocity dispersion functions. Alternatively, both phase and group data can be used in the tomographic inversion processes described herein.

[0026] In some examples, the method comprises, outputting the updated model of the subsurface target volume to an output device.

[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 program comprising instructions which, when the program is executed by a computer, cause the computer to perform any of the methods disclosed herein.BRIEF DESORPTION OF THE DRAWINGS

[0029] 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 subsurface target volume 100 showing surface points 110 defining in-situ velocity measurement locations;Figure 2 shows a plurality of receivers arranged on a surface above a subsurface target volume;Figure 3 is a perspective view of a shear wave velocity model;Figure 4A is an overhead view of a straight ray path on a surface above a subsurface target volume, with a shear wave velocity model superimposed thereon;Figure 4B is an overhead view of a curved ray path on a surface above a subsurface target volume, with a shear wave velocity model superimposed thereon;Figure 4C is an overhead view of an elliptical Fresnel kernel representing a surface wave path above a subsurface target volume, with a shear wave velocity model superimposed thereonFigure 4D is an overhead view of a curved Fresnel kernel representing a surface wave path above a subsurface target volume, with a shear wave velocity model superimposed thereon;Figure 5 is a schematic diagram of a method for determining one or more physical properties of a subsurface target volume;Figure 6 is a schematic diagram of a method for determining a surface wave velocity model;Figure 7 is a schematic diagram of a method for determining one or more physical properties of a subsurface target volume;Figure 8A - 8C shows the standard deviation, ov. , plotted on a surface comprising a surface point 110 and an additional reference point 112 at which no measurement was obtained. Figures 8A to 8C show standard deviations with scaling factors 160, 50 and 25 respectively.Figure 9 is an example of a shear wave velocity model resulting from the methods described herein; andFigure 10 is a schematic diagram of a computing device suitable for performing the methods described herein.DETAILED DESCRIPTION OF THE DRAWINGS

[0030] This detailed description describes, with reference to Figures 1 and 2, an approach to measuring structural properties of a ground volume using receivers. Next, with reference to Figures 3-9, novel methods for performing tomographic inversion using signals detected by receivers are disclosed. Finally, a computing device that may be used to perform the disclosed methods is described with reference to Figure 10.

[0031] The following examples will be described in the context of a receiver array on a surface of the ground volume and one or more subsurface receivers, 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 receivers, accelerometers, velocimeters, seismometers, vibration sensors and / or transducers. The disclosed methods may be applied to any suitable set of signals.

[0032] The methods and systems disclosed herein relate generally to processing signals detected by receivers on a surface. In one particular example, the receivers are placed on a ground surface and the detected signals are ambient seismic noise signals. Processing of these signals provides useful insight into the structure of the surface and subsurface target volume on which the receivers are placed, as described in more detail below. The inventors have found that by supplementing the signals detected by receivers on the surface with in-situ measurements providing additional empirical data, an increased accuracy and / or resolution of information about the subsurface target volume can be achieved. However, there is a requirement to correctly represent the impact of in-situ measurements on cells at varying distances from where these measurements were taken. The methods and systems disclosed herein provide an approach for performing processes to produce accurate, high-resolution information about the subsurface target volume in a computationally efficient and practically feasible manner.

[0033] Before turning to the details of the disclosed methodology, some background relating to determination of surface and subsurface properties using receivers will first be provided. 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 subsurface target 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 subsurface material, and its ability to support structures positioned above and / or through the volume.

[0034] 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, causing 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.

[0035] 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, whichpropagate 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 other surface waves (e.g. Love waves and Scholte waves) may be measured and harnessed in the systems and methods described herein.

[0036] 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.

[0037] 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 over a longer travel distance. The resulting measurement results therefore generally have a higher signal quality (signal-to-noise ratio) than body wave studies.

[0038] Referring now to Fig. 1 , a cross-sectional view of a subsurface volume 100 is shown. P- and S- waves travel through the volume 100 as body waves. A surface 102 extends above the subsurface volume, defined by an x, y plane. Surface waves propagate along the surface 102. In some examples, the subsurface volume 100 may comprise one or more boreholes, for example, a borehole 106 extends from the surface 102, through the subsurface volume 100 (i.e extends in a z dimension). Of course, it is understood that in practise, the borehole 106 may move through an x and a y dimension in addition to a z dimension, however for illustrative purposes a borehole extending through just the z dimension is discussed.

[0039] At an example point, a schematic representation of a particle oscillation (due to Rayleigh wave propagation) at a surface above a target subsurface volume is shown. As illustrated, the oscillation of the particle P is partly vertical and partly parallel to the direction of propagation. The resulting particle movement is therefore substantially ellipsoidal.

[0040] At a surface 102 above the volume 100, a plurality of receivers 104 are arranged. Receivers 104, located at surface 102 can be configured to measure the vertical component of the oscillation shown schematically at the example point.

[0041] 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” because each receiver may be offset from its neighbours in the grid in the z-dimension. However, such a grid arrangement of receivers will be referred to as a 2D array herein for simplicity.

[0042] 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 to provide 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.

[0043] A surface wave travelling across surface 102 will cause vertical movement at a plurality of receivers 104 as waves travel across the surface.

[0044] To determine the shear velocity, Vs, from the observation of surface waves (in particular Rayleigh waves), the dispersive behaviour of the surface waves can be measured. Surface waves are dispersive, which means that their velocity is dependent on frequency. Usually, seismic velocities increase with depth in the earth. As a consequence, 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 the subsurface target volume that the material properties of the volume can be determined.

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

[0046] 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 at which the energy of the wave propagates through the volume. The group velocity is measured by determining the wave propagation between a (synthetic) transducer pair and is a frequency-dependent point property in the volume, which is dependent on depth. The group velocity is obtained as a time-of-flight measurement between a (virtual) source and a receiver.

[0047] The phase velocity is the speed at which a particular frequency component of a wave travels. As such, the phase velocity is expressed as a function of frequency. To measure the phase velocity, at least two measurement nodes are chosen to measure the waves propagating through the volume to determine relative time of flight between the receivers 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 obtained by a 2D transform (such as slant-stack, Radon, FK, or the like) of an array of recorded waveforms (timedistance space).

[0048] Referring to again to Figure 1 , each receiver 104 provides a measurement node for measuring the vertical component of passing surface waves. The receivers 104 can be configured to measure vibrations due to ambient seismic noise. That is, the background wavefield due to natural or man-made noise (rather than an impulse point such as an explosion or hammer drop used in active methods).

[0049] By cross-correlating the passive noise signals measured at a pair of receivers, the Green’s function for the pair can be obtained, which represents the wavefield as if one of the pair were a virtual source and the other of the pair were a receiver.

[0050] In some examples one or more sub-surface receivers are arranged at one or more sub-surface locations below the surface 102. In some examples one sub-surface receiver may be defined by one surface point (i.e only one sub-surface receiver is found at a given x,y coordinate).

[0051] In some examples, wherein the subsurface volume comprises a borehole, one or more subsurface receivers may be arranged throughout a length of the borehole 106 within the volume 100. In such cases, a plurality of receivers 108 are arranged to give an in-situ velocity profile by obtaining multiple velocity measurements defined by one point on the surface 110 (i.e. a plurality of receivers are used to obtain an in-situ profile in a z direction at a given x,y coordinate).

[0052] The subsurface receivers 108, may be configured to measure a subsurface signal from which in-situ empirical time travel data can be obtained. The receiver 108 may be any type of in-situ instrument for obtaining in-situ empirical time travel data. For example, in-situ velocity values may be obtained through a variety of techniques such as geotechnical (seismic CPT) or geophysical (wireline logging) invasive methods.

[0053] Figure 2 shows a plurality of virtual receiver pairs across a surface above a subsurface target volume. Ray paths 202 between receiver pairs are indicated, in particular, the ray paths from a single central receiver near the centre of the array of receivers and each other receiver in the array. The background shading and contour rings indicate the travel-time field from the central receiver to the other receivers. There are also corresponding ray paths between each receiver and all other receivers, i.e. between every pair of receiver, which are not depicted in Figure 2 for simplicity.

[0054] Each pair of receivers can provide a signal at a first location to be cross-correlated with a signal at a second location to reproduce a virtual receiver pair using the principle of interferometry. In particular, 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 subsurface target volume, as if it were induced by an impulse point source, which is equal to Green’s function.

[0055] In some examples, receivers 104a and 104b in Figure 1 form a first receiver pair in the array at surface 102. 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.

[0056] 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 receiver 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.).

[0057] 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 are studied to determine the dispersive behaviour of the surface waves.

[0058] With reference to Figure 3, a first shear wave velocity model 300 has a grid of cells mxyacross the surface above the subsurface target volume, comprising columns mxi, mX2, mX3, etc. extending in the x-direction and rows miy, rri2y, msy, etc extending in the y-direction. Each cell defines an area of the surface. For example, each cell may define a 5m by 5m square; other example options include a 1 m by 1 m square, or a 10m by 10m square. In other words, the shear wave velocity model comprises a plurality of cells arranged in a two-dimensional grid 300. The two-dimensional grid spans at least the area of the surface above the subsurface target volume. Choosing a smaller cell area increases the resolution of the model. Each cell also includes a volume extending vertically below the area of the surface. The model 300 may extend infinitely below the area of the surface or to a predetermined depth below thesurface for which the shear wave velocity value has a notable influence on the propagation of surface waves propagating on the surface. For example, the model 300 may be defined up to a depth of 50m, 100m, 200m or 300m.

[0059] In some examples, the first shear wave velocity model 300 is not a grid of square cells but instead comprises cells having a rectangular shape, a rhombic shape or otherwise tessellating shapes including non-uniform shapes or a combination of different shapes. In some examples, the first shear wave velocity model may comprise tessellating Voronoi cells and a random number of cells are initiated upon performing a method 700 to determine the velocity of each cell. This is described in more detail with reference to Figure 7.

[0060] Each cell is associated with a shear wave velocity value, an example of a physical property value, which represents the expected value of shear wave velocity in the actual subsurface target volume of interest. The shear wave value carries depth information, either in that the shear wave value is constant throughout the volume below the cell area or in how the shear wave value varies with depth (the z direction in Figure 3). For example, the defined shear wave value for each cell may be explicitly a function of depth, either a continuous function or a series of values with associated ranges of depth for each value. In another example, the shear wave velocity value may be a function of frequency, which corresponds to depth information because surface wave propagation is influenced by the physical properties of the subsurface volume up to approximately one wavelength depth. In other words, surface waves at low frequencies are affected by physical properties at deeper depths than surface waves at high frequencies.

[0061] In other examples, the model may be a physical property other than shear wave velocity, such as compressional wave velocity, density, elastic modulus, shear modulus, or, if a viscoelastic model is being used, optionally also viscosity quality factors Qsand QP. In general, the model may define multiple physical property values.

[0062] With reference to figures 4A to 4D, various representations of a wave path between a (virtual) source and a receiver are depicted. The wave path is the path taken by a surface wave travelling from the source to the receiver. As described below with reference to Figure 4A and 4B, a wave path may be represented by a straight or curved ray connecting the source and receiver. Alternatively, as described below with reference to Figure 4C, the wave path may be represented by an elliptical Fresnel zone, also known as a Fresnel kernel, between the source and receiver pair. As another alternative, and as described below with reference to Figure 4D, the wave path may be represented by a curved (“banana”) Fresnel zone.

[0063] With reference to Figure 4A, a ray path is defined as a straight line between two receivers at surface locations A and B on the model 300. A straight ray path model assumes a laterally constant velocity (in other words, the velocity only changes with depth). The straight ray path model provides a sufficiently accurate approximation for the actual path traversed by a wave between a source and receiver in many situations. The straight ray path approximation is also mathematically and computationally efficient compared to other more complex techniques. The ray path signifies the motion of a surface wave as it travels from A to B (or B to A) according to ray theory. In particular, a ray path is defined as the direction of propagation of a surface wave, i.e. the direction perpendicular to wave frontsin wave theory or perpendicular to travel-time contours. The ray path shown in Figure 4A traverses seven cells of the model 300, numbered 1 to 7. The ray path comprises seven ray segments, each having a length Li according to the length for which the ray path travels through each cell, i.e. Li to L? in Figure 4A.

[0064] With reference to Figure 4B, a ray path is defined as a curved line between two receivers at surface locations C and D. The ray path shown in Figure 4B traverses six cells of the model 300, numbered 1 to 6. The ray path comprises six ray segments, each having a length Li according to the length for which the ray path travels through each cell, i.e. Li to Le in Figure 4B. The trajectory of the ray path may be calculated using the principle of least time (Fermat’s principle) between locations C and D based on the shear wave velocity values of the cells in model 300.

[0065] With reference to Figure 4C, a wave path is represented by an elliptical Fresnel zone between two receivers at surface locations E and F. The Fresnel zone is the region around a geometrical ray (such as the straight ray depicted in Figure 4A) that follows the direction of the gradient (spatial derivative of travel time) from the receiver (F) to the source (E).

[0066] In more detail, and as would be appreciated by the skilled person any wave propagating along a path between a (virtual) source at E and a receiver at F will have some component of the wave that propagates off-axis (not along the straight line connected E to F). This component of the wave that propagates off-axis may be deflected by the wave propagation medium, and some of the deflected wave is then directed towards the receiver at F. The travel time of a wave travelling along a direct (straight) path will therefore differ from a wave travelling along a deflected path, meaning that the direct-path wave and deflected-path wave will arrive at the receiver out of phase. The phase difference may lead to destructive interference when the phase difference is half of a wave period (or one and a half wave periods, or two and a half, or any half of an odd-integer amount). Conversely, the phase difference will lead to constructive interference when the phase difference is between zero and one half of a wave period (or between one and one half, or between any integer n-1 and n - % wave periods). In other words, the n-th Fresnel zone is defined as the region in which a wave that is deflected at one point will arrive between n-1 and n- % wavelengths out of phase with a wave travelling along a straight line between a source and receiver pair. In examples of the present disclosure, only the first order Fresnel zone may be considered for the wave path between the source and the receiver, since components of the wave in higher order Fresnel zones become vanishingly small. The general physical principles of Fresnel zones for any wave propagating through any medium are described by B.D Guenther (2005) “Encyclopaedia of Modern Optics”, in which the description of Fresnel zones would be understood by the skilled person to apply to any propagation medium. Fresnel zones for waves travelling through a subsurface target volume are also described by Shibo Xu and Alexey Stovas (2018) “Fresnel zone in VTI and orthorhombic media”.

[0067] As illustrated in Figure 4C, which represents a first order Fresnel zone for a straight-ray approximation, the Fresnel zone is an elliptical shape, with the (virtual) source location E and receiver F being the focal points of the ellipse. In other examples, such as that depicted in Figure 4D, the Fresnel zone may be a curved, so-called banana shaped region that follows the path of a curved ray between a virtual source G and receiver H. Determining such a curved Fresnel zone follows the sameconsiderations as described above with reference to the elliptical Fresnel zone of Figure 4A, except that the differential travel time and phase difference is determined with respect to a curved ray (c.f. Figure 4B) instead of a straight-line ray between the source and receiver.

[0068] The Fresnel zone depicted in Figures 4C covers eight cells of the model 300, numbered 1 to 8. The Fresnel zone comprises eight sub-regions labelled 7?xto R8, each region being the region of the Fresnel zone in the respective cell 1 to 8. Each sub-region 7?xto R8of the Fresnel has a respective sensitivity value that is calculated uniquely for the Fresnel zone. As would be understood by the skilled person, the sensitivity value is a sensitivity to the change in local shear wave velocity (or any other physical property that is being modelled. Such techniques for calculating the sensitivity value for each sub-region of the Fresnel zone (i.e. the sensitivity value forthe Fresnel zone for each cell that the Fresnel zone covers) are well known to the skilled person.

[0069] The same illustration of sub-regions is also shown for the curved Fresnel zone depicted in Figure 4D, which passes through cells 1 to 8 and consists of sub-regions ft jto R8(different from the subregions depicted in Figure 4C). Each of the sub-regions depicted in Figure 4D has a respective sensitivity value, which the skilled person would understand may be different from the sensitivity values for each sub-region depicted in Figure 4C.

[0070] For each of a plurality of virtual receiver pairs, such as those depicted in Figure 2, a respective Fresnel zone may be determined. Each Fresnel zone consists of a respective plurality of sub-regions for each cell that is covered by the Fresnel zone, and each of those sub-regions for each zone has a respective sensitivity value.

[0071] With reference to Figure 5, a method 500 for determining one or more physical properties of a subsurface target volume includes receiving 502 a plurality of signals detected by receivers arranged on the surface, e.g. as described with reference to Figure 1 . In general, one signal is received from each receiver, although signals from only a subset of the total number of receivers may be received in some circumstances. The signals show the vertical oscillations of the ground measured at the respective receiver as caused by surface waves from ambient seismic noise or an active source, as described above with reference to Figure 1 . The signals may include metadata about the location and time of recording.

[0072] The signals may be received directly from the receivers or via one or more intermediary devices. For example, a computing device (as described below with reference to Figure 10), may be situated locally to the array of receivers for the signals to be transmitted via any form of wired or wireless communication. Alternatively, the signals may be received at a location far from the location of the array of receivers for performing the method 500 remotely from the array of receivers, e.g. by internet communication or transporting a physical computer-readable medium with a recording of the signals saved thereon. In some examples, additional processing steps may be performed on the signals before or after the receiving 502 of the signals to optimise the signals for further processing.

[0073] Example surface receivers can include velocimeters or accelerometers. An example of a particular mechanism for such receivers includes a ferromagnetic mass on a spring moving within an electric coil in response to movement of the surface, thereby inducing an electric current proportional to ground velocity which can be measured.

[0074] Example sub-surface receivers can include accelerometers or piezo-electric (pressure) sensors.

[0075] The method 500 also comprises a step 504 cross-correlating the signals between the plurality of receivers arranged on the surface to obtain empirical travel time data for surface waves between a plurality of virtual receiver pairs. As described above, this step involves cross-correlating the passive noise signals measured at a pair of receivers to obtain a correlated Green’s function, which represents the wavefield as if one of the pair were a virtual source and the other of the pair were a receiver. As would be well understood by the skilled person, information indicative of phase velocity as a function of frequency and information indicative of group velocity as a function of frequency can be derived from the cross-correlation process. In step 504, the cross-correlation of the signals between the plurality of receivers is carried out to obtain group and / or phase traveltime data for a plurality of frequencies. The frequencies may be continuous so as to obtain group velocity and / or phase velocity dispersion information, or the traveltime data may be obtained for a plurality of finite frequencies. Group traveltime refers to the time taken for a group wave to travel from the virtual source to the receiver. Phase traveltime refers to the time taken for a specific phase component of a wave to travel from the virtual source to the receiver. Since the distance between the virtual source and receiver is known, the group and / or phase velocity or slowness for a plurality of frequencies (which is the inverse of velocity) can be derived from the correlated Green’s function for a receiver pair. Slowness and traveltime are effectively equivalent, with a scaling factor provided by the distance between source and receiver.

[0076] Based on at least the group and / or phase travel time data for the plurality of frequencies and the in-situ velocity measurement determined by the sub-surface receivers, a model of the physical property of the subsurface target volume can be determined. As described above with reference to Figure 3, the model comprises a plurality of cells arranged in a two-dimensional grid. Each cell of the model has a respective physical property value (such as a shear wave velocity), which may be include depth profile. Collectively, the cells of the model therefore provide a three-dimensional model of the physical property, since each cell provides a one-dimensional model as a function of depth.

[0077] Determining the model of the physical property of the subsurface target volume comprises selecting 506 a subset of the plurality of receiver pairs for which empirical traveltime information has been obtained. As described above with reference to Figures 4A - 4D, a wave path between each receiver pair traverses two or more cells of the model. In principle, any and all receiver pairs may be suitable for determining the model. However, increasing the number of wave paths used will increase the computational power needed or else cause very long computation times. Therefore, in practice, it is usually beneficial to select a subset of receiver pairs, with the number to be selected depending on the computational power available or other practical concerns. Several principles can be used to limit the number of selected pairs. As a preliminary point, due to the reciprocity of cross-correlation, the wave path from a location A to location B is the same as location B to location A and therefore only the choice of paired receivers is needed without regard to which is treated as the virtual source. A first way to deselect pairs is to omit those which do not have suitable traveltime data, i.e. the cross-correlation of the signals from that pair of receivers does not show an identifiable wave passing through both receivers which can be used to find the travel time therebetween. Another way to select wave paths is to omitpairs with low-quality picks, e.g. having a low signal-to-noise ratio, having large uncertainty, being an outlier compared to adjacent picks etc. Any picks that would give non-physical or non-geological results are also excluded, as well as any picks that would be implausible for the specific subsurface target volume under investigation. Lastly, if the remaining number of pairs is still more than can feasibly be used, a subsample of receiver pairs may be chosen, e.g. every nthpair or a random selection, etc. In such a situation, other subsamples of suitable receiver pairs may be used in succession to increase the total number of pairs used with exceeding memory limitations or available computational power. Regardless of what approaches are used to limit the number of receiver pairs (if necessary), the end result is a selection of receiver pairs and corresponding wave paths for which empirical group and / or phase traveltime information, such as velocities for a plurality of different frequencies or dispersion functions (e.g. a group velocity dispersion function and / or a phase velocity dispersion function) has been obtained.

[0078] As described above with reference to Figures 4A - 4B, a wave path between a (virtual) source and receiver pair may be a ray path that extends in a straight line between the source and receiver. In some examples, at least some of the ray paths extend in a curved line between the respective source and receiver, as described above with reference to Figure 4B. At least some of the selected wave paths will traverse at least two cells of the first model, having a ray segment length in each cell which it traverses (including the cells which the ray path starts and ends in). Curved ray paths may be determined according to initial physical property values for the cells of the model. By using curved ray paths determined according to the physical property values of the cells of the model, the ray paths more closely follow the actual ray path that a wave would travel in between the source and the receiver. Therefore, this approach improves the accuracy of the results of the later method processes because it corresponds more closely to the physical reality of surface waves passing through the subsurface target volume.

[0079] As described above with reference to Figures 4A - 4B, in other examples the wave path is represented by a Fresnel zone which provides a region of many possible paths taken by the wave as it travels from the source to the receiver. The Fresnel zone provides a region that covers two or more cells through which a possible wave traverses, with each cell having a respective region of the Fresnel zone with an associated sensitivity value defined by the Fresnel zone. The Fresnel zone may be determined based on initial physical property values for the cells of the mode. Since the Fresnel zone defines a region of possible paths taken by a surface wave travelling from the source to the receiver, this approach improves the accuracy of the resulting model for each cell because it corresponds more closely to the physical reality of surface waves passing through the subsurface target volume and scattering within the volume.

[0080] Determining the model of the physical property of the subsurface target volume comprises a step 508 of retrieving one or more sub-surface seismic parameter measurements from one or more subsurface receivers 108, defined by a first point on the surface 110. For example, in-situ velocity values may be obtained. In some examples, the one or more subsurface receivers 108 may be inserted into the surface to obtain in-situ velocity values, vloc. In some examples, the one or more subsurface receivers 108 may be invasive subsurface receivers. In some examples, the one or more subsurface receivers 108 may be inserted into a borehole 106 to obtain in-situ velocity values, v,()C,

[0081] The one or more subsurface receivers 108 may obtain in-situ velocity values, vloc, through any appropriate known method, for example through geotechnical (seismic CPT) or geophysical (wireline logging) invasive methods.

[0082] The method 500 further comprises a step 510 of determining a standard deviation, aloc, , corresponding to the first point xloc,ylocat which the in-situ velocity values are obtained.

[0083] In some examples, the standard deviation, aloc, is defined based on a quality of the surface of the subsurface target volume 100 to impact the distance from a first point on the surface to which the in-situ velocity measurements have probabilistic influence. In some examples, the standard deviation, aloc, is obtained based on the in-situ velocity values obtained at point xloc,yloc. In further examples, a standard deviation is ol0Cis picked from a set of known experimental values that the inventors have found yield good results. For example, the inventors have found a pre-determined standard deviation value of 0.01 to 0.1 yields good results. In some examples, the standard deviation is set to 0.04.

[0084] At step 510, a scaling factor is determined. As with the standard deviation, aloc, the scaling factor, a, is defined based on a quality of the surface of the volume 100 to impact the distance which the in-situ velocity measurements have probabilistic influence. The scaling factor, a, may be defined within a range of 1 to 200, where a scaling factor of 1 may be representative of a perfect surface i.e. a surface on which the measurement points influences never decreases. The scaling factor is described in more detail with reference to step 708 of method 700 and Figures 8A to 8C.

[0085] Steps 508 and 510 and hence method 500, pertain to retrieving in-situ parameter measurements corresponding to only a first point on the surface, however, it should be understood that in-situ surface seismic parameter measurements can be obtained from more than one point, in which cases a standard deviation and scaling factors corresponding to each point is obtained. For example, second of in-situ parameter measurements may be obtained by receiving at least one subsurface signal detected by at least one subsurface receiver at a second location within the subsurface target volume, wherein the second location defines a second point on the surface above the subsurface volume, wherein the prior probability density function is dependent on a distance of each of the cells from the second point on the surface.

[0086] In some examples receiving a plurality of subsurface signals corresponding to a plurality of locations may be received, each location defining a point on the surface, wherein the prior probability density function is dependent on a distance of each of the cells from the plurality of points on the surface.

[0087] Determining the model of the physical property of the subsurface target volume further comprises performing 512 tomographic inversion on the travel time data for each of the selected receiver pairs whilst constraining the inversion using the in-situ velocity values obtained in step 508. Generally speaking, the tomographic inversion is used to derive a value of the physical property for each cell of the model, from group or phase travel time data obtained for each of a plurality of receiver pairs. As described above, the physical property value for each cell may be a single value, a distribution or provided as the physical property as a function of depth. Details of various processes for carrying out the tomographic inversion are described below with reference to Figure 6 - 7. It is understood that any appropriate process of carrying out the tomographic inversion could be followed, using the measurement location dependent prior probability density function as a starting point to the method described belowwith reference to Figure 7. As described throughout this text, the tomographic inversion can be carried out based on empirical group traveltime information or empirical phase dispersion information. In the context of phase dispersion, phase velocity or slowness information may be converted to equivalent traveltime information so that the tomography can be performed on such phase traveltime information.

[0088] In some examples, the tomographic inversion step 512 involves process 600 depicted in Figure 6. As described in more detail below, this process tomographically maps the traveltime information from each receiver pair into the cells of the physical property model. The result of the process is therefore an empirical model of group and / or phase velocity (depending on the traveltime data used), with each cell of the model having a respective group or phase velocity value or distribution. This process 600 is carried out for each of a plurality of frequencies to obtain group or phase velocity values for each cell for each frequency.

[0089] Process 600 is the first stage (i.e. the tomography stage) of a two-stage tomographic inversion. Process 600 comprises a step 600 of obtaining an initial velocity model. The initial velocity model comprises initial values of group velocity and / or phase velocity for each cell for each of a plurality of frequencies. The initial velocity model sets initial (group or phase) velocity values for each cell (for each given frequency), which the process 600 will refine using the empirical traveltime data in an iterative process. Accordingly, it is not essential for the initial velocity model and respective initial velocity values to be a highly accurate or high-resolution, although a more accurate initial model may render the tomographic process 600 faster or more accurate at mapping traveltimes to each cell. A more accurate initial model may also reduce the risk of finding a local-minimum rather than a global minimum in a iterative gradient-descent method, although such an issue could be solved using Monte-Carlo methods. In some examples, the initial velocity model is determined based on a user input, e.g. according to historic data or map information indicating possible physical property values across the subsurface target volume. Alternatively, an arbitrarily chosen typical value of the group / phase velocity can be used for each cell as a starting point. In these examples, the arbitrary model may be selected based on an estimation of the physical properties of the subsurface target volume. The tomographic process 600 further comprises the step 604 of determining modelled group or phase travel times for each of the selected receiver pairs using the initial group or phase velocity model. This is carried out by identifying the wave path from the source to the receiver (for example as a straight ray, curved ray or elliptical / curved Fresnel zone) and identifying the cells that are traversed by the wave path. The modelled travel time based on the initial velocity model is then determined based on the known distance between the source and receiver and the velocity values of each cell traversed by the wave path.

[0090] The tomographic process 600 further comprises the step 606 of determining an error value indicative of the difference between the modelled travel times (determined at step 604) and the empirical travel times (obtained at step 504) for each of the selected receiver pairs. For example, the error value could be a simple difference between the modelled and empirical traveltimes, also called a residual, at each frequency for each receiver pair. In some examples, the error value could be a combination of all the differences between modelled and empirical traveltimes for every receiver pair. In general determining the error value is part of an iterative process, for example, in a least-squares inversion method where the squares of the residuals are calculated in order to be minimised through iterations.Other forms of inversion processes use different error values in order to provide feedback to the initial group or phase velocity model.

[0091] The process 600 further comprises the step 608 of determining an updated velocity model using the error value. The updated model is generally in all aspects the same as the initial velocity model except for new velocity values associated with at least some of the cells. In other words, the updated model is an updated version of the initial velocity model taking into account the determined error value between the empirical and modelled group or phase traveltimes for each receiver pair. This feedback process may involve a least-squares method, Markov-chain Monte Carlo method or other inversion technique to iteratively update the initial model based on an updated model. Process 600 may be repeated for each of a plurality of finite frequencies.

[0092] Steps 604 to 608 are typically all part of a subroutine of the first-stage tomographic process 600, which is then iterated according to the inversion method such as least-squares inversion. The iterative nature of the process 600 is illustrated by the dashed arrow in Figure 6, which indicates that the updated model determined at step 608 is used to determine new modelled traveltimes for each receiver pair at step 604. In other words, each time an updated model is determined using the error value resulting from the initial model and resulting modelled traveltimes for each receiver pair, the resulting updated model is then used as the initial model for the next iteration. The iteration continues until the error value reaches an end condition, for example, the error value falling below a threshold absolute value of the difference between the modelled and empirical traveltimes, or falling below a threshold of proportional difference between the modelled and empirical traveltimes. Another end condition, which could be used alone or in combination with the threshold error value end condition, is that the changes to the initial model to produce the updated model for an iteration falls below a threshold amount or proportion. This is so that if the iteration reaches a settled minimum error value, the iteration can end as further iterations will not make significant accuracy improvements.

[0093] Once the tomography of the first-stage process 600 has been carried out to obtain a velocity model comprising a group or phase velocity value for each cell (for each frequency), the tomographic inversion of step 512 can proceed to the second-stage (the inversion) of the two-stage tomographic inversion to obtain the resultant model of the physical property of the target subsurface volume. The second-stage of the inversion process 700 is depicted in Figure 7. As with the first-stage tomographic process 600, the second stage inversion process 700 is carried out for each of a plurality of frequencies to obtain shear velocity values for each cell. Inversion process 700 is iterated until a number of iterations meets a first end condition, and wherein the updated velocity model of each iteration is used as the initial velocity model of the next iteration. For example, the first end condition may be a number of iterations in the region of hundreds to thousands of iterations. Of course, any appropriate number of iterations may be determined before iterating through process 700.

[0094] The inversion process 700 comprises a first step 702 of obtaining an initial physical model. The initial physical model is an initial model of the physical property of the subsurface volume, wherein the initial model comprises initial physical property values for each of the cells of the first plurality of cells. The initial physical model may be a shear wave velocity model 300 as described above with reference to Figure 3, and / or be an initial physical model having any of the features or variations described abovewith reference to Figure 3. The initial model sets initial physical property values for each cell of the model, which the method 700 will refine using the group / phase velocity model obtained from iterative tomographic process 600 using a further iterative process. Accordingly, it is not essential for the initial physical model and initial physical property values to be a highly accurate or high-resolution model of the subsurface target volume, although a more accurate first model may increase the expected accuracy of the end result of the method or decrease the computational time to reach the end result.

[0095] Process 700 further comprises a step 704 of determining modelled surface wave velocities for each cell based on the initial model. For each receiver pair of the subset, determining a modelled travel time using the initial velocity model values associated with the two or more cells that are traversed by the surface wave path. In more detail, a forward modelling approach is used to derive the phase and / or group velocities for each cell of the initial physical model, using the respective physical property value for the given cell. For example, the physical property values of each cell may be a shear wave velocity. The values of shear wave velocity of for each cell can be used to calculate a corresponding phase velocity dispersion function (and thus a phase velocity for a given frequency) using the propagation matrix method introduced by Thomson ((1950). “Transmission of elastic waves through a stratified solid medium”, Journal of applied Physics, 21 (2), 89-93) and Haskell ((1953) “The dispersion of surface waves on multilayered media". Bulletin of the seismological Society of America, 43(1), 17-34). The skilled person would be well aware of this and other approaches that can be used to forward model group or phase velocities (or dispersion functions) from shear wave velocity values (or functions of depth).

[0096] The inversion process further comprises the step 706 of determining an error value indicative of the difference between the modelled velocities (determined at step 704 based on the initial physical model) and the empirical velocities (obtained from process 600). The error value may be determined in a similar manner to that described above in relation to step 606 in which an error value is determined between modelled and empirical traveltimes. As with the tomographic process 600, in general determining the error value at step 706 is part of an iterative process, for example, in a least-squares inversion method where the squares of the residuals are calculated in order to be minimised through iterations.

[0097] Other forms of inversion processes use different error values in order to provide feedback to the initial model of the physical properties of the subsurface target volume, for example Bayesian inference may be used to determine an error value to provide feedback on the initial model of the physical properties of the subsurface target volume provided in step 702. In such examples, an error value indicative of the difference between the modelled velocities (determined at step 704 based on the initial physical model) and the empirical velocities (obtained from process 600) for each cell of the model may be a likelihood. Step 706 may comprise determining a likelihood of the modelled velocities using Bayesian inference on the modelled travel times (determined at step 704) and the obtained empirical travel time data (obtained from process 600). The likelihood is a known function used in Bayesian inference and is well defined in the art.

[0098] The process 700 further comprises the step 708 of determining an updated physical model based on the error value determined at step 706. For example, step 708 may comprise determining aposterior probability based on the likelihood of the updated velocity model and the prior probability density function.

[0099] The prior probability density function is defined as - V) , wherein a prior covariance matrix C„ is adiagonal matrix with a trace comprising of a standard deviation of an ith model cell, Cp(i, i) = o- ., matrix Nmcontains normalization factors of the probability density function of the ith cell, vector V is comprised of the seismic parameter values of all model cells n and vector V comprises the measured seismic parameter values.

[0100] The inventors found that by defining a location dependent standard deviation ov. of the prior probability density function describing the parameter value at a location x^ , the prior probability function behaves like a Gaussian at locations close to the in-situ measurement location (imposing probabilistic preference for specific parameter values (the in-situ physical parameters) and progressively further away approaches a uniform distribution.

[0101] The standard deviation of the prior probability density function is defined as<7V= min L , where alocis the standard deviation at thefirst point on the surface xloc,y and a represents a scaling factor which determines the rate with which the influence of the first point on the surface decreases with distance as <JVincreases until threshold value L. As discussed above with reference to step 510 of Figure 5, the scaling factor is determined based on the confidence of the x,y influence of the in-situ velocity measurements as the distance from the point of measurement increases. Figure 8A - 8C shows the standard deviation, ov. , plotted on a surface comprising a first surface point 110 and an additional reference point 112 at which no measurement was obtained. Figures 8A to 8C show standard deviations ov. with scaling factors 160, 50 and 25 respectively, Each of Figures 8A to 8C have ol0C= 0.04 and threshold L = 12.

[0102] Figure 8A shows the standard deviation, ov. with a scaling factor of 160. The standard deviation increases rapidly with distance away from the first point on the surface 1 10 such that the prior information only imposes constraining velocity values at and near where they were measured at the first surface point 110. The additional reference point 112 is out of reach of the influence of the first surface point 112 and therefore constrained by a uniform probability density function (i.e. not constrained by the in-situ velocity measurements at all).

[0103] Figure 8B shows the standard deviation, ov. with a scaling factor of 50. The standard deviation increases more modestly, extending the influence for the velocity values measured at the first surface point 110 on the modelling domain. The additional reference point 112 is still out of reach of the influence of the first surface point 110 and therefore constrained by a uniform probability density function.

[0104] Figure 8C shows the standard deviation, ov. with a scaling factor of 25. The standard deviation increase slowly and as such, the additional reference point 112 falls within the influence for the velocity values measured at the first surface point 110 on the modelling domain. Therefore, the probability density function at point 112 displays Gaussian like behavior imposing a probability preference for thevelocity value measured at the first surface point 110. As discussed above with reference to Figure 5, the value for the scaling factor, a is determined based on the confidence on the influence of the in-situ velocity values observed at the first surface point 110.

[0105] The updated physical model is generally in all aspects the same as the initial physical model except for new physical property values associated with at least some of the cells. In other words, the updated model is an updated version of the initial physical model taking into account the determined error value between the empirical and modelled surface wave velocities. This feedback process may involve a least-squares method, Markov-chain Monte Carlo method or other inversion technique to iteratively update the initial physical model based on an updated physical model. In some examples, step 708 comprises determining a posterior probability density function of the velocity model based on the likelihood determined in step 706, the prior probability density function and the initial model provided in step 702 using Bayesian inference. Based on the initial model and the posterior probability density function, an updated physical model may then be determined based on an acceptance probability wherein the acceptance probability indicates whether the initial or posterior probability density function is accepted and used as a starting point in the next round of iterations.

[0106] In such examples, calculating the acceptance probability comprises obtaining a ratio of the posterior probability of the updated velocity model to the posterior probability of the initial velocity model, generating a random number between 0 and 1 and accepting the updated model if the ratio is greater than the random number or reverting to the initial model if the ratio is smaller than the random number.

[0107] In some examples calculating the acceptance probability comprises accepting the larger of the two probabilities.

[0108] Steps 704 to 708 are typically all part of a subroutine of the second-stage inversion process 700, which is then iterated according to the inversion method such as least-squares inversion, or Markov Chain Monte Carlo method, such as the reverse jump Markov Chain Monte Carlo Method (rj-McMC method). The iterative nature of the process 700 is illustrated by the dashed arrow in Figure 7, which indicates that the updated model determined at step 708 is used to determine new velocities for each cell at 704. In other words, each time an updated model is determined using the error value resulting from the initial model and resulting modelled velocities for each cell, the resulting updated model is then used as the initial model for the next iteration. The iteration continues until the error value reaches an end condition, such as an end criterion described above in relation to Figure 6.

[0109] For example, when iterating steps 704 to 708 according to the rj-McMC method, the first shear wave velocity model may comprise tessellating Voronoi cells and a random number of cells are initiated upon performing a method 700 to determine the velocity of each cell. In such examples, step 706 may further comprise perturbating the velocity model to form one or more perturbations of the initial velocity model and then selecting a perturbation as an updated model. The perturbations may be one or more of: a velocity change of a cell of the plurality of cells; a location change of a cell of the plurality of cells; a cell nucleus birth; a cell nucleus death;a change in a first noise parameter; or a change in a second noise parameter.

[0110] The likelihood and posterior probability density function may then be determined using the perturbed model.

[0111] With reference to Figure 9, an example output of the methods described herein is a final shear wave velocity model showing a subsurface target volume. The subsurface target volume extends in x and y directions to a depth (z direction) of 100m. The values of shear wave velocity are shown by the shading in the plot and the transition between regions of different shear wave velocity are visible, which indicates the different composition or structure of portions of the subsurface target region. Using the methods described herein, a final shear wave velocity model can be determined with higher resolution and accuracy without resulting in an unfeasible compute time. Such a subsurface model can be used to better understand the suitability of the subsurface target volume for supporting man-made structures on top of or in the subsurface target volume.

[0112] Optionally, additional physical properties of the subsurface target volume can be determined from the final shear wave velocity model, e.g. calculated using one or more of equations 1 and 2. These further physical properties can be recorded in the shear wave velocity model itself, or outputted separated to an output device

[0113] .In additional to or alternatively to shear wave velocity, other physical properties may be used which are related to shear wave velocity. For example, the longitudinal wave (P-wave) velocity Vp and shear wave (or transverse I S-wave) velocity Vs are related to other physical properties of the ground elastic modulus X, shear modulus p, and density p by the following equations from linear elasticity theory:Equation 1 (longitudinal wave velocity):VP =JEquation 2 (shear wave velocity):

[0114] Figure 10 shows a block diagram of one implementation of a computing device 1000 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 ofinstructions (sequential or otherwise) that specify actions to be taken by that machine. 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.

[0115] The example computing device 1000 includes a processor 1002, a main memory 1004 (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 1006 (e.g., flash memory, static random access memory (SRAM), etc.), and a secondary memory (e.g., a data storage device 1018), which communicate with each other via a bus 1030.

[0116] Processor 1002 represents one or more general-purpose processors such as a microprocessor, central processing unit, or the like. More particularly, the processor 1002 may be a complex instruction set computing (CISC) microprocessor, reduced instruction set computing (RISC) microprocessor, very long instruction word (VLIW) microprocessor, processor implementing other instruction sets, or processors implementing a combination of instruction sets. Processor 1002 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 1002 is configured to execute the processing logic (instructions 1022) for performing the operations and steps discussed herein.

[0117] The computing device 1000 may further include a network interface device 1008. The computing device 1000 also may include a video display unit 1010 (e.g., a liquid crystal display (LCD) or a cathode ray tube (CRT)), an alphanumeric input device 1012 (e.g., a keyboard or touchscreen), a cursor control device 1014 (e.g., a mouse or touchscreen), and an audio device 1016 (e.g., a speaker).

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

[0119] The data storage device 1018 may include one or more machine-readable storage media (or more specifically one or more non-transitory computer-readable storage media) 1028 on which is stored one or more sets of instructions 1022 embodying any one or more of the methodologies or functions described herein. The instructions 1022 may also reside, completely or at least partially, within the main memory 1004 and / or within the processor 1002 during execution thereof by the computer system 1000, the main memory 1004 and the processor 1002 also constituting computer-readable storage media.

[0120] 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, forexample, 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.

[0121] 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.

[0122] 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 a special-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.

[0123] 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 or to perform certain operations described herein.

[0124] 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).

[0125] Unless specifically stated otherwise, as apparent from the following discussion, it is appreciated that throughout the description, discussions utilizing terms such as "receiving”, “determining”, “comparing”, “calculating”, “averaging,” “identifying”, “updating”, “solving”, “outputting” 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.

[0126] 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 withreference to the appended claims, along with the full scope of equivalents to which such claims are entitled.

Claims

CLAIMS1 . A method for determining a physical property of a subsurface target volume, the method comprising: receiving a plurality of signals detected by a plurality of receivers arranged on a surface above the subsurface target volume, wherein each respective signal of the plurality of signals is detected by a respective receiver of the plurality of receivers; cross-correlating signals between the plurality of receivers to obtain empirical travel time data for surface waves between a plurality of receiver pairs; receiving at least one subsurface signal detected by at least one subsurface receiver at a first location within the subsurface target volume, wherein the first location defines a first point on the surface above the subsurface volume; and determining a model of the physical property of the subsurface target volume, the model comprising a first plurality of cells, wherein each cell has a respective physical property value and wherein the plurality of cells spans at least the area of the surface, wherein determining the model comprises: selecting a subset of the plurality of receiver pairs, wherein a surface wave path between each source and respective receiver traverses two or more cells of the first plurality of cells; performing tomographic inversion on the empirical travel time data between each receiver pair in the subset to obtain the physical property value for each cell, wherein performing the tomographic inversion comprises obtaining a posterior probability density function corresponding to each cell of the first plurality of cells, wherein obtaining the posterior probability density function comprises: obtaining an initial velocity model; and defining a prior probability density function indicative of velocity for each of the cells of the first plurality of cells; for each receiver pair of the subset, determining a modelled travel time using the initial velocity model values associated with the two or more cells that are traversed by the surface wave path; and determining a posterior probability density function of the velocity model based on the modelled travel times and the prior probability density function, wherein the prior probability density function is dependent on a distance of each of the cells from the first point on the surface.

2. The method of claim 1 further comprising: receiving at least one subsurface signal detected by at least one subsurface receiver at a second location within the subsurface target volume, wherein the second location defines a second point on the surface above the subsurface volume, wherein the prior probability density function is dependent on a distance of each of the cells from the second point on the surface.

3. The method of claim 1 wherein performing the tomographic inversion further comprises iteratively: perturbating an initial velocity model to form one or more perturbations of the initial velocity model; selecting a perturbation as an updated velocity model; determining a likelihood of the updated velocity model based on the modelled travel times and the obtained empirical travel time data; determining a posterior probability based on the likelihood of the updated velocity model and the prior probability density function; calculating an acceptance probability based on a posterior probability of the updated velocity model and the posterior probability of the initial velocity model.

4. The method of claim 3, further comprising, before determining a likelihood of the updated velocity model, either: for each receiver pair of the subset, determining a modelled travel time using the updated velocity model values associated with the two or more cells that are traversed by a surface wave path; or determining updated velocity model values based on the previously computed surface wave paths from the initial velocity model.

5. The method of claims 3-4, wherein the steps of perturbating the initial velocity, selecting an updated velocity model, determining a modelled travel time, determining a likelihood of the updated velocity and calculating an acceptance probability are iterated until a number of iterations meets a first end condition, and wherein the updated velocity model of each iteration is used as the initial velocity model of the next iteration.

6. The method of claims 4-5, wherein calculating the acceptance probability comprises: accepting the updated model if the posterior probability of the updated model is greater than the posterior probability of the initial model.

7. The method of claim 5-6 wherein calculating the acceptance probability comprises: obtaining a ratio of the posterior probability of the updated velocity model to the posterior probability of the initial velocity model; generating a random number between 0 and 1 ;accepting the updated model if the ratio is greater than the random number; reverting to the initial model if the ratio is smaller than the random number.

8. The method of claims 5-7 wherein the first end condition is a predetermined number of iterations.

9. The method of claims 5 - 8, wherein the steps of perturbating the initial velocity, selecting an updated velocity model, determining a likelihood of the updated velocity and calculating an acceptance probability are performed using the reversible-jump Markov chain Monte-Carlo algorithm.

10. The method of any preceding claim, wherein the perturbations may be one or more of: a velocity change of a cell of the plurality of cells; a location change of a cell of the plurality of cells; a cell nucleus birth; a cell nucleus death; a change in a first noise parameter; or a change in a second noise parameter.11 . The method of any preceding claim wherein the prior probability density function is defined by:wherein the prior covariance matrix Cpis a diagonal matrix where the trace consists of a standard deviation of an ith model cell, Cp(i, i) =, matrix Nmcontains the normalization factors of the probability density function of the ith cell, vector V is comprised of the seismic parameter values of all model cells n and vector V contains the measured seismic parameter values.

12. The method of any preceding claim wherein the standard deviation of the prior probability density function is defined as<7V= min aloc* exp (a L , where av. is the standard deviation of theGaussian probability density function describing the parameter value at location xt,y , alocthe standard deviation at the first point on the surface xloc,y and a represents a scaling factor which determines the rate with which the influence of the seismic parameter value measured at the first point on the surface decreases with distance as a,, increases until threshold value L.

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 method of any preceding claim14. 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 method of any of claims 1-12.

Images