NL2033825B1 - Tomographic inversion - Google Patents

Abstract

The disclosure provides a method for determining physical properties of a subsurface target volume, comprising: receiving a plurality of signals detected by receivers arranged on a surface above the subsurface target volume; determining a first model of target volume, wherein the first model comprises a plurality of cells each associated with one or more initial physical property value; selecting a plurality of ray paths between pairs of receivers, wherein each ray path traverses two or more cells; determining an empirical dispersion function for each ray path of the plurality of ray paths using the plurality of signals; determining a modelled dispersion function for each ray path using the initial physical property values; determining an error value indicative of the difference between the modelled and empirical dispersion functions for each ray path; and determining a second model using the error value with new physical property values for the plurality of cells.

Description

TOMOGRAPHIC INVERSION

TECHNICAL FIELD

[0001] The disclosure relates to methods and systems for analysing a target region beneath a surface of the earth. More particularly, the disclosure relates to a method and system for determining one or more material properties of the subsurface target region based on noise measured at or near the surface. 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 sub-surface ground parameters. In particular, there is a need for systems and methods that can be used to model the properties of a target volume beneath the surface of the earth to provide information useful for infrastructure planning. Determination of sub-surface ground properties during the early planning phase of construction projects reduces uncertainty during the location determination, foundation design, and construction phases of a project. This in turn reduces delays, overspend, and unnecessary use of material resources (e.g. concrete) during construction.

[0003] One key parameter for the determination of ground characteristics in a volume of interest is shear-modulus and shear-velocity Vs. The shear velocity Vs is the velocity at which a shear wave moves through the material and is controlled by the shear modulus of the material. The relationship between shear-velocity and shear modulus G is defined by Vs = ¥G/p, where p is the density of the material.

Measurement of Vs therefore provides a valuable insight to the material properties of a sub-surface ground region.

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

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

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

OVERVIEW

[0007] According to a first aspect of the present disclosure, there is provided a method for determining one or more physical properties 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; determining a first model of the subsurface target volume, wherein the first model comprises a plurality of cells each associated with one or more initial physical property value; selecting a plurality of ray paths between pairs of receivers of the plurality of receivers, wherein each ray path traverses two or more cells of the plurality of cells; determining an empirical dispersion function for each ray path of the plurality of ray paths using the plurality of signals; determining a modelled dispersion function for each ray path of the plurality of ray paths using the initial physical property values associated with the two or more cells of the first model traversed by that ray path; determining an error value indicative of the difference between the modelled dispersion function and the empirical dispersion function for each ray path; and determining a second model using the error value, wherein the second model comprises new physical property values for the plurality of cells. By performing the features of determining a modelled dispersion function and comparing with the empirical dispersion function (by way of an error wave) for a section of ray paths, as opposed to individually per cell, the accuracy of determining the one or more physical property of the subsurface target volume is increased. This is because the update to the first model to produce the second model is performed in a single process for a ray path, rather than using an initial tomographic process to find dispersion functions for individual cells using many ray paths which introduces additional approximation or inaccuracy. This 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.

[0008] The determining the first model may comprise one or more of: receiving the first model from a user interface, e.g. as input by a user; retrieving the first model from a memory; calculating the first model using the signals; or generating the first model, such as randomly generating the first model or having a default first model. The cells of the plurality of cells of the first model may create a uniform square grid, rectangular grid, or other tessellating shape grid. The cells of the plurality of cells of the first model may have a first dimension of less than 20m, less than 10m, less than 5, less than 2m, or less than 1m. Typically, the initial physical property value is defined as a function of depth, or a function of frequency (which carries depth information for Rayleigh wave propagation).

[0009] The selecting the plurality of ray paths between pairs of receivers of the plurality of receivers does not exclude there being additional receivers located between any given selected pair of receivers.

Although each ray path of the selection of ray paths traverses two or more cells of the plurality of cells, this does not exclude other unselected ray paths also being used in combination with the ray paths which traverse two or more cells. The determining the empirical dispersion function for each ray path may include cross-correlating the signals from the respective pair of receivers for that ray path.

[0010] The plurality of signals may include data specifying the locations of the receivers, either absolute location (latitude, longitude) or relative location to a common reference point of the receivers. For example, each signal may include data specifying the location of the receiver at which that signal was detected. The method may include determining the location of each receiver of at least the pairs of receivers of selected ray paths, e.g. determining using the plurality of signals, retrieving from a memory, receiving from a user interface, etc.

[0011] The error value may be any single value, set of parameters, list of values, matrix, etc. which is indicative of the difference between the empirical dispersion function and modelled dispersion function.

The error value may be an aggregate or list of error values each specific to a particular ray path. In implementations where the first model is a shear wave velocity model, the first error may be a difference between an initial shear wave velocity value of the first model and a shear wave velocity value resulting from an inversion using the empirical dispersion function, for each cell traversed by the ray path. The determining the second model may then include adding (or subtracting, as appropriate) the difference to the initial shear wave velocity value to result in the new shear wave velocity value which would result in the empirical dispersion function (or at least closer thereto than the initial shear wave velocity value).

[0012] The determining the second model using the error value may include performing an operation onthe initial physical property value based on the error value for each ray path or may include replacing the initial physical property value with a new physical property value which would result in the empirical dispersion function (or at least closer thereto than the initial physical property value).

[0013] In some implementations, the determining the modelled dispersion function may comprise: averaging, for each ray path, the one or more initial physical property values of the two or more cells that ray path traverses to produce an averaged physical property value; and calculating, for each ray path, the modelled dispersion function using the averaged physical property value for that ray path. In some implementations, the determining the modelled dispersion function may comprise: calculating, for each ray path, a cell dispersion function for each of the two or more cells traversed by that ray path to provide a plurality of cell dispersion functions; and averaging, for each ray path, the plurality of cell dispersion functions for that ray path. 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, i.e. the length of the portion of the ray path (straight or curved) which is within the bounds of that cell. By averaging over cells in each ray path, the model is more accurate as it more closely reflects how the ray path between a pair of receivers behaves. This process also provides a feasible way to achieve a high-resolution model, i.e. having smaller cell sizes and therefore more precise resulting model. The weighted average approach also increases the accuracy by accounting for the different affect different cells have on each ray path.

[0014] In some implementations, the determining the second model may comprise updating the initial physical property value of each cell according to a respective segment length of the ray path in that cell.

For example, the change in physical property value of each cell may be weighted as described above.

These techniques mean that the resulting model corresponds more closely to the actual physical properties of the subsurface target volume and may also decrease the time or iterations needed to find the final model to fit the detected signals.

[0015] In some implementations, the determining the second model may comprise: determining partial derivatives of the error value with respect to the one or more physical property 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 first model; updating the first model according to the desired changes to produce the second model. This provides a precise way to implement the changes to the first model to produce the second model, and therefore produces a second model which more closely maps onto the actual physical properties of the subsurface target volume.

[0016] In some implementations, an error value associated with the second model may be less than the error value associated with the first model, i.e. as the second model gets closer to the actual physical properties of the subsurface target volume. In some implementations, the determining the modelled dispersion function based on the first model, determining the error value, and determining the second model may be iterated until the error value meets an end condition, wherein the second model of each iteration is used as the first model of in the next iteration. This further increases the accuracy of the final second model, as each iteration aligns the model more closely to the detected signals (and, in particular, the dispersion functions determined therefrom). The end condition may include one or more of: the error value of a latest iteration is less than a predetermined threshold error value; and a difference between the error value of a latest iteration and the error value of a preceding iteration is less than a predetermined difference threshold. The end conditions control the balance between optimising the second model to the actual physical properties in the subsurface target volume and reducing computational time.

[0017] In some implementations, the determining the modelled dispersion function, the determining the error, and the determining the second model may be performed using a least-squares inversion method.

This is possible because the determination of the dispersion functions for each ray path and provides is less computationally intense compared to common alternatives (such as Monte Carlo methods).

[0018] In some implementations, the empirical dispersion function and the modelled dispersion function may be group velocity dispersion functions and/or phase velocity dispersion functions. That is, the modelled dispersion function and the group velocity dispersion function are both group velocity dispersion functions, they are both phase velocity dispersion functions, or they both include both group velocity and phase velocity dispersion functions. By using both phase velocity and group velocity dispersion functions, a second model is determined with a closer fit to the physical properties of the subsurface target volume.

[0019] In some implementations, the plurality of ray paths may include at least one curved ray path.

This approach acknowledges that the variation in physical properties between pairs of receivers may result in a curved ray path, and therefore performing the methods disclosed herein on curved ray paths improves the accuracy of the resulting model. In general, a combination of straight ray paths and curved 5 ray paths may be used.

[0020] In some implementations, the one or more initial physical property value for each cell may be a shear wave velocity value. This is a useful physical property to model because it has a particularly significant influence on surface wave propagation.

[0021] In some implementations, the method may comprise outputting the second model of the target subsurface volume to an output device. In some implementations, the method may comprise calculating, for one or more cells of the second model, a compressional wave velocity value and/or a density value using the shear wave velocity value. That is, any combination of shear wave velocity, a compressional wave velocity value and a density value may be included in the model.

[0022] In some implementations, at least one receiver may comprise; a geophone; a ferromagnetic mass on a spring moving within an electric coil in response to the movement of the surface inducing an electric current proportional to ground velocity; a velocimeter; an accelerometer; a seismometer; a vibration sensor; or a transducer. In general, the plurality of receivers may be all of the same type of receiver or may be a mix of any combination of types of receivers.

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

[0024] 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 DESCRPTION OF THE DRAWINGS

[0025] Disclosed implementations will now be described by way of example to illustrate aspects of the disclosure and with reference to the accompanying drawings, in which:

Figure 1 shows a cross-sectional view of a sub-surface target region with a plurality of receivers arranged at the surface;

Figure 2 shows a plot identifying the ray paths between a plurality of receivers;

Figure 3 is a perspective view of a shear wave velocity model;

Figure 4A is an overhead view of a 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 ray path on a surface 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 an example of a shear wave velocity model resulting from the methods described herein;

Figure 7 is a schematic diagram of a computing device suitable for performing the methods described herein.

DETAILED DESCRIPTION

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

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

[0028] The methods and systems disclosed herein relate generally to processing signals detected by geophones on a surface. In one particular example, the geophones are placed on a ground surface and the detected signals are ambient noise signals. Processing of these signals provides useful insight into the structure of the surface and subsurface target volume on which the geophones are placed, as described in more detail below. Due to the potentially very large number of signals being processed, an increased accuracy and/or resolution of information about the subsurface target volume can be achieved. However, using a large number of signals also means that computationally efficient methods are useful to ensure that processing of the signals is feasible within a practical length of time. 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.

[0029] Before turning to the details of the disclosed methodology, some background relating to determination of surface and subsurface properties using geophones 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.

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

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

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

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

[0034] Referring now to Fig. 1, a cross-sectional view of a sub-surface volume 100 is shown. A surface 102 extends above the sub-surface volume. Surface waves propagate along the surface 102, as shown at point A, where a schematic representation of a particle oscillation (due to Rayleigh wave propagation) at the surface above the target sub-surface volume is shown. As illustrated, the oscillation of the particle

P is partly vertical, partly in the direction of propagation and movement is therefore substantially ellipsoid.

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

[0036] The receivers 104 are configured to measure ambient noise. That is, the wavefield present due to background noise (as opposed to noise from an active source such as a hammer drop or explosion).

The background noise may be natural (e.g. due to lapping ocean waves, wind, and other naturally occurring vibrations) or cultural (e.g. due to human activity, including traffic, machinery, etc.).

[0037] As shown in Figure 1, the receivers 104 are arranged in a grid array at the surface, the grid extending in two directions. It should be noted that the surface above the target region may, in many cases, not be planar. The array of receivers 104 may therefore not be truly “2-dimensional” (each receiver may be offset from its neighbours in the grid in the z-direction, as well as the x- and y-directions).

However, such a grid arrangement of receivers will be referred to as a 2D array herein. It will also be appreciated that in the context of sub-surface ground study, the term 2D array may also be used to denote an array configured to capture 2D information (e.g. a line of receivers configured 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.

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

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

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

[0041] To determine the shear velocity from the observation of surface waves (in particular Rayleigh waves), the dispersive behaviour of the surface waves is studied. Surface waves are dispersive, i.e. their velocity is dependent on frequency. Since seismic velocities increase with depth in the earth, normal surface wave dispersion shows a decrease of surface-wave velocity with increasing frequency.

It is by studying the behaviour of surface waves at a surface above a volume that the ground properties of the volume can be determined.

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

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

[0044] The phase velocity is the velocity at which the phase of any one frequency component of the wave travels. As such, the phase velocity is expressed as a function of frequency. To measure the phase velocity, at least two measurement nodes (e.g. receivers) are chosen to measure the waves propagating through the volume to determine relative time-of-flight between the nodes for different frequencies. The result is the phase velocity as a function of frequency as averaged over the volume between the two measurement nodes. Phase velocity is acquired as a point in 2D phase space

(dispersion spectrum} which is itself obtained by a 2D transform (such as slant-stack, Radon, FK, or the like) of an array of recorded waveforms (time-distance space). The wavelength of the surface wave is indicative of its propagation depth. As a result, the phase velocity as a function of frequency is indicative of the shear velocity Vs as a function of depth.

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

[0046] Receivers 104a and 104b form a first receiver pair in the array at surface 102. Neither receiver 104a nor 104b represents a point source for noise recorded at the other of receiver 104a, 104b.

However, by cross-correlating the received signal at receiver 104a and 104b, receivers 104a and 104b can act as a (virtual) source receiver pair, where each receiver of the pair records a signal as though the signal had originated at the other of the pair. The cross-correlation is performed using the principle of interferometry.

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

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 source- 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).

[0048] In Figure 1, only one pair of receivers is labelled (104a, 104b). However, it will be appreciated that for each receiver 104 in the array, every other receiver in the array may act as the other half of a source receiver pair. In this manner, the Green’s function for each source receiver pair may be obtained.

The Green's functions across the plurality of virtual source receiver pairs is studied to determine the dispersive behaviour of the surface waves.

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

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

[0051] The starting model for the solution of the inverse problem may be determined in multiple ways since it sets initial physical property values, for example for Vs, within the sub-surface target volume, to be refined by solution of the inverse problem. The starting model can comprise historical data based on known local or regional data, predicted values based on selective invasive test methods, coarse surface level studies or a theoretical model based on local or regional historical knowledge. Although the starting model is to be improved upon and need not be a particularly accurate representation of the sub-surface target region, the more accurately the starting model reflects the physical properties of the sub-surface target region, the more accurate the final model of the ground properties of the target region will be.

Therefore, improvements in the starting model can provide an improved output for the methods and systems described herein.

[0052] Once a starting model has been selected, the inverse problem can be solved via inversion, as will be explained in more detail below.

[0053] With reference to Figure 3, a first shear wave velocity model 300 has a grid of cells mx across the surface above the subsurface target volume, comprising columns mx, Mx2, Ms, etc. extending in the x-direction and rows muy, may, May, 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 1m by 1m square, or a 10m by 10m square. 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 the surface 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.

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

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

[0056] In other examples, the model may 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 Qs and Qs. In general, the model may define multiple physical property values.

[0057] With reference to Figure 4A, a ray path is defined as a straight line between two geophones at surface locations A and B on the model 300. 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 fronts in 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 L: according to the length for which the ray path travels through each cell, i.e. Li to Ly in Figure 4A.

[0058] With reference to Figure 4B, a ray path is defined as a curved line between two geophones 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. L1 to Ls 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.

[0059] 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 geophones arranged onthe surface, e.g. as described with reference to Figure 1. In general, one signal is received from each geophone, although signals from only a subset of the total number of geophones may be received in some circumstances. The signals show the vertical oscillations of the ground measured at the respective geophone as caused by surface waves from ambient 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.

[0060] The signals may be received directly from the geophones or via one or more intermediary devices. For example, a computing device (as described below with reference to Figure 7), may be situated locally to the array of geophones 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 geophones for performing the method 500 remotely from the array of geophones, 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.

[0061] Example geophones can include velocimeters or accelerometers. An example of a particular mechanism for such geophones 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.

[0062] The method 500 also comprises determining 504 a first model of the subsurface target volume.

The first model may be a shear wave velocity model 300 as described above with reference to Figure 3, and/or be a first model having any of the features or variations described above with reference to

Figure 3. The first model sets initial physical property values for each cell of the first model, which the method 500 will refine using the received plurality of signals. Accordingly, it is not essential for the first 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.

[0063] In some examples, the first 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, the first model may be determined using the received signals, e.g. by performing an inversion of group velocity or phase velocity dispersion curves between geophones, which can be calculated from cross-correlation of the signals as described above, to find a using a coarser grid or quicker method. As a last resort, an arbitrarily chosen typical value of the physical property can be used for each cell as a starting point.

[0064] The method 500 also comprises selecting a plurality of ray paths out of the total number of possible ray paths between pairs of geophones. In principle, any ray path between two geophones could be suitable for including the method 500. However, increasing the number of ray 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 ray paths between geophones, with the number of rays paths to be selected depending on the computational power available or other practical concerns. Several principles can be used to limit the number of selected ray paths. As a preliminary point, due to the reciprocity of cross-correlation, the ray path from a location A to location B is the same as location B to location A and therefore only the choice of paired geophones is needed without regard to which is treated as the virtual source. A first way to deselect ray paths is to omit ray paths which do not have suitable group dispersion picks, i.e. the cross-correlation of the signals from that pair of geophones does not show an identifiable wave passing through both geophones which can be used to find the travel time therebetween. Another way to select ray paths is to omit ray paths with low-quality picks, e.g. having a low signal-to-noise ratio, having large uncertainty, being an outlier compared to adjacent picks etc.

Lastly, if the remaining number of ray paths is still more than can feasibly be used, a subsample of ray paths may be chosen, e.g. every nt" ray path or a random selection, etc. In such a situation, other subsamples of suitable ray paths may be used in succession to increase the total number of ray paths used with exceeding memory limitations or available computational power. Regardless of what approaches are used to limit the number of ray paths (if necessary), the end result is a selection of ray paths for which dispersion functions (e.g. a group velocity dispersion function and/or a phase velocity dispersion function) can be determined using the signals from the geophones at each end of the respective ray path.

[0065] In some examples, each ray path extends in a straight line between two geophones, as described above with reference to Figure 4A. In some examples, at least some of the ray paths extend in a curved line between two geophones, as described above with reference to Figure 4B. At least some of the selected ray 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). By using curved ray paths determined according to the physical property values of the cells of the first model, the ray paths more closely follow the actual ray path that a wave would travel in between the pair of geophones.

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.

[0066] In some examples, instead of a single straight or curved ray path, a Fresnel zone is determined between the pair of geophones, which covers a range of possible ray paths which a surface wave may travel along between geophones. For the purposes of the later method processes, the Fresnel zone will also be described as a ‘ray path’.

[0067] The method 500 further comprises determining 508 an empirical dispersion function for each ray path, in particular, each ray path selected in the selecting 506 part of the method. The dispersion function may be a group velocity dispersion function or a phase velocity dispersion function (or both may be used). The group velocity dispersion function may be determined as described above with reference to Figure 1, according to any of the usual methods in the art. As a high-level summary, this generally includes cross-correlation between the two signals from the pair of geophones, identifying from the signals a wave which passes both locations of the geophones, finding the time delay at which the wave appears in the cross-correlated signal, and dividing the distance between the geophones by the time delay. This can be done at many different frequencies and plotted as a dispersion curve showing the variation of group velocity or phase velocity as a function of frequency for that ray path. This process is done for each of the ray paths selected in the selecting 506 part of the method 500.

[0068] The method 500 further comprises determining 510 a modelled dispersion function for each ray path using the first model. For example, the values of shear wave velocity of the first model can be used to calculate a phase velocity dispersion function using the propagation matrix method introduced by

Thomson ((1950). “Transmission of elastic waves through a stratified solid medium”, Journal of

[0069] 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), using a model of a stack of homogenous layers with finite thickness (overlying a homogenous half-space).

[0070] In particular, this can be performed using any available solver using a modal approximation method, such as written by Herrmann ((2013) “Computer programs in seismology: An evolving tool for instruction and research” Seismological Research Letters, 84(6), 1081-1088). This approach computes the phase velocity Vphase for a particular frequency w from a stack of homogeneous layers with finite thickness, i.e. a 1-dimensional velocity profile. The first step consists of establishing the boundaries of where the phase velocity will lie between, e.g. the minimum and maximum Vs present in the model. For an arbitrary value of w, a trial value of wavenumber k is then tested and iteratively changed with a set increment until the eigenvalue solution is found. This eigenvalue can subsequently be used to calculate the solution to the eigenfunction, which will be tested against an iteratively changing value of Vphase to find the root of the solution and therefore a phase velocity dispersion function.

[0071] Hermanns modal approximation code takes as inputs a list of layered elastic parameters and the frequency values of interest and outputs a dispersion function as the solution of the eigenvalue problem defined by the matrix propagator method introduced by Haskell and Thomson. Herrmann’s solver is widely used and available as a software package from Saint Louis University Earthquake

Center, “Computer Programs in Seismology” page hitps.//www.eas.stu.edu/eqc/eqeeps.htmi, along with extensive manuals and instructions for using this software (in addition to the 2013 paper “Computer programs in seismology: An evolving tool for instruction and research” references above). Alternative solvers are also available via GitHub, such as htips://github.com/xin2zhang/MCTomo based on 3D

Monte Carlo tomography (Zhang, X., Curtis, A., Galetti, E., & de Ridder, S., 2018. “3-D Monte Carlo surface wave tomography’, Geophysical Journal International, 215(3), 1644-1658), or https://github.com/keurfoniuu/disba by Keurfon Luu (also available at zenodo.org titled keurfonluu/disba: disba v0.5.1 and which, according to its documentation page, implements “a subset of codes from

Computer Programs in Seismology (CPS) in Python compiled just-in-time with numba”.

[0072] In addition to shear wave velocity, other physical properties may be used for the model which are related to shear wave velocity. For example, the longitudinal wave (P-wave) velocity Vp and shear wave (or transverse / S-wave) velocity Vs are related to other physical properties of the ground elastic modulus A, shear modulus p, and density p by the following equations from linear elasticity theory:

Equation 1 (longitudinal wave velocity): yo [rm p

Equation 2 (shear wave velocity): -§ p

[0073] In summary, there are established methods for determining a phase velocity dispersion function (and/or group velocity dispersion function) from a stack of homogenous layers (constant is x and y directions) with finite thickness and associated shear wave velocity values (or related physical properties) for each layer. Accordingly, one way to use the methods developed by Thomson and Haskell and Herrmann is to use multiple ray paths between geophones that intersect within a particular cell of the model and combine the dispersion functions of the ray paths to determine a dispersion function for that particular cell. The dispersion function for that cell can then be used for ‘inversion’ using the homogenous layers approach using the shear wave velocity values for that cell, i.e. refining the model values in that cell to better produce the dispersion function found for that cell (and doing likewise for all cells individually). However, the inventors have developed a different approach wherein, instead of performing the inversion on each cell, the inversion is done for a whole ray path using the multiple shear wave velocity values in the cells that the ray path traverses. This approach has numerous advantages, including improved computational efficiency, allowing for geophone arrays with more geophones and/or higher resolution models, and increased accuracy of the model determination. For example, because the shear wave velocity model is refined in a single step from the dispersion functions found between geophones, this removes the extra step of travel-time tomography using the ray paths for each cell, which is an extra source of approximation and therefore inaccuracies because the inversion is independent for each cell (therefore losing the interdependence of crossing rays in the inversion step)

In addition to removing an extra step in the process, thereby reducing the amount of computational power needed, performing the inversion in a single step also means that a least-squares method can be used rather than more computationally intensive Monte-Carlo inversion methods. This in turn means that the inversion compute time is reduced or, for the same amount of compute time, more geophone signals can be processed or a higher resolution model can be produced.

[0074] In some examples, the determining 510 the modelled dispersion function for each ray path is performed by averaging the shear wave velocity values of the cells which that ray path traverses. In particular, for shear wave velocity, the 1/ Vs values are averaged because the travel-time through each cell is inversely proportional to velocity. For other physical properties of the cells, such as density, a simple arithmetic mean of the values in each cell may be suitable. In either case, the values are still a function of depth. The averaged shear wave velocity (or other averaged physical property) is then used in the Herrmann modal approximation method to determine the modelled dispersion function for the whole ray path. This averaging process can be implemented in the solver software as an additional function prior to running the inversion routine.

[0075] In some examples, the averaging of the shear wave velocity values is a weighted average, weighted according to the ray segment length in each cell, as set out in equation 6, which applies to straight ray paths or curved ray paths. For Fresnel zone ray paths, the weights include a term corresponding to sensitivity of the ray path to variations in the local velocity at different parts of the

Fresnel zone (i.e. there is higher sensitivity at the centre of the zone than the peripheries). Using a weighted average increases the accuracy of the calculation as it more closely corresponds to the actual values that the ray path will experience between the pair of geophones.

Equation 3:

Lay

Vsqup Lap el Vs where:

Vsas is the average shear wave velocity for the ray path A to B (as a function of depth};

N is the number of cells of the first model that the ray path A to B traverses;

Li is the ray segment length of the in cell along the ray path;

Lag is the total ray path length from A to B, i.e. the sum of all L;

Vsi is the shear wave velocity value for the itt cell along the ray path (as a function of depth).

[0076] In some examples, rather than first averaging the shear wave velocity values along the ray path length, the modal approximation method is performed to produce a dispersion function for each cell and then the cell dispersion functions for each cell which the ray path traverses are averaged to produce the modelled dispersion function for the ray path. This can be done by averaging 1 / Vphase for phase velocity dispersion function (or conversely 1 / Vgroup for group velocity dispersion function) for each cell along the ray path, optionally including weightings analogously to the averaging described above. This averaging process can be implemented in the solver software as an additional function after running the inversion routine

[0077] The method further comprises determining 512 an error value indicative of the difference between the modelled dispersion function (determined 510 from the model} and the empirical dispersion function (determined 508 from the detected signals) for each ray path. For example, the error value could be a simple difference between the modelled and empirical dispersion functions, also called a residual, at each frequency for each ray path. In some examples, the error value could be a combination of all the differences between modelled and empirical dispersion functions for every ray path. In general, the determining 512 the error value is part of the inversion process, for example, in a least-squares 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 first model.

[0078] The method further comprises determining 514 a second model using the error value. The second model is generally in all aspects the same as the first model except for new shear wave velocity values {or whichever physical property values are used) are associated with at least some of the cells.

In other words, the second model is an updated version of the first model taking into account the determined error value between the empirical and modelled dispersion functions. This feedback process is typically built into the inversion method, e.g. a least-squares method, Markov-chain Monte Carlo method etc. and performed as part of the inversion routine.

[0079] In some examples, the determined error value indicates that the empirical dispersion function along a ray path shows larger values for phase velocity (or group velocity) than the modelled dispersion function, either across the function or in one or more frequency bands. Accordingly, this indicates that the values of shear wave velocity values along that ray path are less than the true physical properties of the subsurface target volume. As such, determining the second model using the error value comprises recording new shear wave velocity values for the cells which that ray path traverses which are greater than the corresponding values in the first model. Conversely, if the empirical dispersion function has smaller values than the modelled dispersion function, then the values of the shear wave velocity values along that ray are greater than the true values in the subsurface volume and the model should be updated accordingly.

[0080] One way to update the first model to determine the second model is to find an empirical shear wave velocity average for the wave path through inversion of the dispersion function, and to scale each of the shear wave velocity values from the first model according to a difference between the empirical shear wave velocity average and the average shear wave velocity of the first model for that wave path.

In more detail, the update factor for the first model for each wave path maybe defined as the ratio of: the difference between the empirical shear wave velocity average and the average shear wave velocity of the first model for that wave path; and the average shear wave velocity of the first model for that wave path. In practice, this entails doing an arithmetic mean of 1 / Vs, also referred to as the ‘slowness’, because cells with smaller shear wave velocity will have a bigger impact on the average shear wave velocity over the whole wave path. Optionally, the updating of the first model to produce the second model includes weighting the change to the shear wave velocity values for each cell according to the ray segment length of the wave path in that cell, according to the approach explained above with reference to determining 510 the modelled dispersion function part of the method 500.

[0081] In some examples, the shear wave values of the first model are updated to produce the second model using a sensitivity matrix method. In this approach, partial derivatives of the error value (or residual) are determined with respect to the shear wave velocity and/or other parameters of the first model, such as depth of boundaries between shear wave velocity values in a given cell, longitudinal wave velocity, density, etc. The partial derivatives are assembled into a sensitivity matrix and a linear system of equations is solved to determine how to change the parameters of the first model such that the overall error value is reduced for the second model with updated values.

[0082] Once the updates to the values of the first cell along a first ray path are calculated, these updates can be applied to produce the second model directly. Alternatively, the updates to values of all cells as influenced by each ray path calculation and inversion can be combined and applied to the values of the first model at once to produce the second model.

[0083] The parts of the method of determining 510 the modelled dispersion function, determining 512 the error value, and determining the second model using the error value are typically all part of a subroutine of the method 500, which is then iterated according to the inversion method such as least- squares gradient-descent inversion. In other words, each time a second model is determined using the error value resulting from the first model and resulting modelled dispersion functions for each ray path, the resulting second model is then used as the first 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 dispersion functions, or falling below a threshold of proportional difference between the modelled and empirical dispersion functions.

Another end condition, which could be used alone or in combination with the threshold error value end condition, is that the changes to the first model to produce the second 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.

[0084] With reference to Figure 5, although the parts of the method 500 have been arranged and numbered in a sequence, the method is not limited to the particular order in which the parts are presented herein. As an example, there is no requirement for either the empirical dispersion function or the first modelled dispersion function to be performed before the other.

[0085] In some examples, once a final second model has been determined, e.g. once the iterations satisfy the end condition, the final second model may be sent to an output device. For example, the final second model may be displayed on a monitor or other user interface, or sent via wired or wireless communication to another computing device for further processing or display.

[0086] With reference to Figure 6, 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.

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

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

[0089] The example computing device 700 includes a processor 702, a main memory 704 (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 706 (e.g., flash memory, static random access memory (SRAM), etc.}, and a secondary memory {(e.g., a data storage device 718), which communicate with each other via a bus 730.

[0090] Processor 702 represents one or more general-purpose processors such as a microprocessor, central processing unit, or the like. More particularly, the processor 702 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 702 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 702 is configured to execute the processing logic (instructions 722) for performing the operations and steps discussed herein.

[0091] The computing device 700 may further include a network interface device 708. The computing device 700 also may include a video display unit 710 (e.g., a liquid crystal display (LCD) or a cathode ray tube (CRT)), an alphanumeric input device 712 {e.g., a keyboard or touchscreen), a cursor control device 714 (e.g., a mouse or touchscreen), and an audio device 716 (e.g., a speaker).

[0092] It will be apparent that some features of computer device 700 shown in Figure 7 may be absent.

For example, one or more computing devices 700 may have no need for display device 710 {or any associated adapters). This may be the case, for example, for particular server-side computer apparatuses 700 which are used only for their processing capabilities and do not need to display information to users. Similarly, user input device 712 may not be required. In its simplest form, computer device 700 comprises processor 702 and memory 704.

[0093] The data storage device 718 may include one or more machine-readable storage media (or more specifically one or more non-transitory computer-readable storage media) 728 on which is stored one or more sets of instructions 722 embodying any one or more of the methodologies or functions described herein. The instructions 722 may also reside, completely or at least partially, within the main memory 704 and/or within the processor 702 during execution thereof by the computer system 700, the main memory 704 and the processor 702 also constituting computer-readable storage media.

[0094] The various methods described above may be implemented by a computer program. The computer program may include computer code arranged to instruct a computer to perform the functions of one or more of the various methods described above. The computer program and/or the code for performing such methods may be provided to an apparatus, such as a computer, on one or more computer readable media or, more generally, a computer program product. The computer readable media may be transitory or non-transitory. The one or more computer readable media could be, for example, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, or a propagation medium for data transmission, for example for downloading the code over the Internet.

Alternatively, the one or more computer readable media could take the form of one or more physical computer readable media such as semiconductor or solid state memory, magnetic tape, a removable computer diskette, a random access memory (RAM), a read-only memory (ROM), a rigid magnetic disc, and an optical disk, such as a CD-ROM, CD-R/W or DVD.

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

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

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

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

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

[00100] It is to be understood that the above description is intended to be illustrative, and not restrictive.

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

Claims (15)

Conclusions 1. A method for determining one or more physical properties of a subsurface target volume, the method comprising: receiving a plurality of signals detected by a plurality of receivers disposed at a surface above the subsurface target volume, each respective one of the plurality of signals being detected by a respective one of the plurality of receivers, determining a first model of the subsurface target volume, the first model comprising a plurality of cells associated with one or more initial physical property values; selecting a plurality of ray paths between pairs of receivers of the plurality of receivers, each ray path passing through two or more cells of the plurality of cells; determining an empirical dispersion function for each ray path of the plurality of ray paths using the plurality of signals; determining a modeled dispersion function for each ray path of the plurality of ray paths using the initial physical property values associated with the two or more cells of the first model through which that ray path passes; determining an error value indicative of the difference between the modeled dispersion function and the empirical dispersion function for each ray path; and determining a second model using the error value, the second model including new physical property values for the plurality of cells. 2. The method of claim 1, wherein determining the modeled dispersion function comprises: averaging, for each ray path, the one or more initial physical property values of the two or more cells through which that ray path passes to produce an average physical property value; and calculating, for each ray path, the modeled dispersion function using the average physical property value for that ray path. The method of claim 1, wherein determining the modeled dispersion function comprises: calculating, for each ray path, a cell dispersion function for each of the two or more cells through which that ray path passes to provide a plurality of cell dispersion functions, and averaging, for each ray path, the plurality of cell dispersion functions for that ray path. 4. A method according to claim 2 or 3, wherein the averaging comprises determining a weighted average according to a respective weight for each cell of the two or more cells, the respective weight for each cell corresponding to a segment length of the ray path in that cell. 5. A method according to any preceding claim, wherein determining the second model comprises: updating the initial physical property value of each cell according to a respective segment length of the ray path in that cell. 6. A method according to any preceding claim, wherein determining the second model comprises: determining partial derivatives of the 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 in the one or more initial physical property values of the first model; updating the first model according to the desired changes to produce the second model. 7. A method according to any preceding claim, wherein an error value associated with the second model is less than the error value associated with the first model. 8. A method according to any preceding claim, wherein determining the modeled dispersion function based on the first model, determining the error value and determining the second model are repeated until the error value satisfies a terminal condition, the second model of each iteration being used as the first model of the next iteration, the terminal condition optionally comprising one or more of the following: the error value of a last iteration is less than a predetermined error threshold; and a difference between the error value of a last iteration and the error value of a previous iteration is less than a predetermined difference threshold. 9. A method according to any preceding claim, wherein determining the modelled dispersion function, determining the error and determining the second model are performed using a least squares inversion method. 10. A method according to any preceding claim, wherein the empirical dispersion function and the modelled dispersion function are group velocity dispersion functions and/or phase velocity dispersion functions. 11. A method according to any preceding claim, wherein the plurality of beam paths comprises at least one curved beam path. 12. A method according to any preceding claim, wherein the method comprises: outputting the second model of the subsurface target volume to an output device. 13. A method according to any preceding claim, wherein the one or more initial physical property values for each cell is a shear wave velocity value, and wherein the method optionally comprises: calculating, for one or more cells of the second model, a compressional wave velocity value and/or a density value using the shear wave velocity value. 14. A method according to any preceding claim, wherein at least one receiver of the plurality of receivers comprises: a geophone; a ferromagnetic mass on a spring that moves within an electrical coil in response to surface motion that induces an electrical current proportional to ground velocity; a speedometer; an accelerometer; a seismometer; a vibration sensor; or a transducer. 15. A system comprising: - one or more processors; one or more memories storing computer-readable instructions configured to cause the one or more processors to perform operations comprising the method of any one of claims 1 to 16; or a computer-readable medium containing instructions which, when executed by one or more data processing devices, cause the one or more data processing devices to perform operations comprising the method of any one of claims 1 to 16.