Method for determining a model of a physical property of a subsurface target volume

Abstract

There is provided methods 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 dividing parts of the methods into subtasks, wherein the subtasks are distributed between a network of computational devices of the network and executed in parallel. Unlocking insights from Geo-Data, the present invention further relates to improvements in sustainability and environmental developments: together we create a safe and liveable world.

Description

METHOD FOR DETERMINING A MODEL OF A PHYSICAL PROPERTY OF A SUBSURFACE TARGET VOLUME TECHNICAL 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 volume beneath a surface of the earth and with improved efficiency.  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 shear-modulus 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 / ρ), where ρ 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, surface-level 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 the recorded waves is 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 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 model of a physical property of a subsurface target volume using a network of computational devices is provided. The model comprises a plurality of cells and each cell has a respective physical property value. 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. The method comprises cross-correlating signals between the plurality of receivers to obtain empirical travel time data for surface waves between a plurality of receiver pairs. The method comprises selecting a set of the plurality of receiver pairs, each respective receiver pair having a respective surface wave path between receivers of the receiver pair. The method comprises performing tomographic inversion on the empirical travel time data between receivers of each receiver pair in the set to obtain a respective physical property value for each cell. The method comprises collating the obtained physical property values into the model. One or more of the cross-correlating, the selecting, the performing tomographic inversion, and the collating is divided into subtasks, wherein the subtasks are distributed between the computational devices of the network and executed in parallel.

[0008] The network may be configured to provide cloud-based computing or, in other words, the computational devices may be part of a cloud computing service. The plurality of computational devices may be an entire set of computational devices that provide a cloud-computing environment, or a subset of computational devices that provide the cloud-computing environment. Accordingly, using a network of computational devices means that at least a portion of the method is performed by a plurality of computational devices connected across the network. For example, the tomographic inversion may be performed by a plurality of computational devices in the network.

[0009] The receiving the plurality of signals detected by the plurality of receivers may include receiving the plurality of signals from a user device in communication with the network, e.g., over the internet, or by accessing a database on which the plurality of signals are stored. In some examples, the method may further include detecting, by the plurality of receivers, the plurality of signals to be used in the method and may comprise pre-processing the detected signals for use in the method. The signals detected by the plurality of the receivers may be ambient seismic signals.

[0010] 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 target volume. Accordingly, the physical property of the subsurface target volume may be shear wave velocity.

[0011] The cross-correlating signals between the plurality of receivers may comprise cross-correlating between a portion of the total number of pairs of receivers, i.e., not every possible pair. As used herein, “empirical” data or information refers to such data or information obtained empirically, i.e., by means of real signals collected at physical receivers positioned on the surface above the subsurface target volume.

[0012] The subtasks into which one or more the cross-correlating, the selecting, the performing tomographic inversion, and the collating is divided into are independent, i.e., each subtask can be completed without requiring other subtasks to be completed first and without requiring additional information from other subtasks, such that the subtasks can be executed in parallel. Accordingly, being executed in parallel means that multiple (i.e., at least two, but there could be 100s or more) subtasks are performed concurrently. The subtasks may be performed anywhere in the network of computational devices (subject to the configuration and control of the network). Each subtask may be associated with metadata to label the subtask and / or provide information describing the subtask. For example, the metadata may include an instruction where to save the result. Each subtask may comprise writing the result of the subtask to a designated storage location, or may send the result to another function for saving or actioning. The result of each subtask may be a physical property value, or set of physical property values to be collated into the model.

[0013] By dividing one or more parts of the method into subtasks to be executed in parallel, the efficiency of the method is improved. For example, using ambient seismic signals typically requires a longer measurement period than signals from active testing. This results in a huge amount of data to be processed in order to determine the sought model of the physical property (e.g., shear wave velocity) of the subsurface target volume. Using the network of computational devices as explained above, with particular parts of the above method divided into subtasks for execution in parallel, can significantly reduce the time to determine the model (depending on the available computational resources). Dividing into subtasks also improves efficiency by arranging the computation taking advantage of the independence of certain parts of the computation.

[0014] The tomographic inversion is a computationally intensive part of the method and therefore is particularly suitable for increased computational speed and efficiency. Tomographic inversion may comprise a tomography part and an inversion part, and either or both of these parts may be divided into subtasks. Additionally or alternatively, the improved speed and efficiency can be achieved in any of the cross-correlating, the selecting, the performing tomographic inversion, and the collating, or any combination thereof. The receiving the plurality of signals may also be divided into subtasks in a similar manner. A given subtask may include part of the cross-correlating, the selecting, the performing tomographic inversion, and the collating, e.g., performing all of these for one frequency subrange. In some examples, each of the one or more of the cross-correlating, the selecting, the performing tomographic inversion, and the collating may be individually divided into subtasks, e.g., the subtasks performing only the cross-correlating and new subtask being created for the selecting or tomographic inversion according to a different division of subtasks.

[0015] In some examples, the subtasks are frequency subtasks, wherein each frequency subtask corresponds to a respective frequency subrange of a total frequency range of the plurality of signals. The tomography part of tomographic inversion may be divided into frequency subtasks. For surface waves travelling through the subsurface target volume, the frequency carries depth information within the subsurface target volume. This is because surface wave propagation is influenced by the physical properties of the subsurface volume up to approximately one wavelength depth, e.g., surface waves at low frequencies are affected by physical properties at deeper depths than surface waves at high frequencies. Therefore the frequency subtasks correspond to computations for different depths.

[0016] In some examples, the subtasks are model cell subtasks, wherein each model cell subtask corresponds to a respective cell of the model. For example, each subtask is associated with a cell of the model to be determined and the cross-correlating and / or selecting may be divided such that each model cell subtask selects pairs of receivers according to which pairs’ surface wave paths pass through the relevant cell of the model. The inversion part of tomographic inversion may then be performed on each cell of the model individually to obtain the respective physical property value.

[0017] In some examples, the subtasks are path subtasks, wherein each path subtask corresponds to a respective surface wave path. For example, the tomographic inversion may be performed for an entire surface wave path between a pair of receivers, and the result then used to determine or modify the physical property values across the cells that the wave path traverses in the model. By using this approach, fewer tomographic inversion calculations may be required, decreasing the compute time.

[0018] In some examples, the subtasks are travel time subtasks, wherein each travel time subtask corresponds to respective obtained empirical travel time data. For example, the selecting and / or performing tomographic inversion may be divided into subtasks according to the travel times obtained during the cross-correlating.

[0019] In some examples, the subtasks are receiver subtasks, wherein each receiver subtask corresponds to a receiver of the plurality of receivers. For example, the cross-correlating and / or the selecting may be divided into subtasks according to receiver.

[0020] In some examples, the subtasks may be a combination of frequency subtasks, model cell subtasks, path subtasks, and travel time subtasks. The overall method may use several different types of subtasks, e.g. different types of subtasks for different parts of the method, or different types of subtasks in the same part. One subtask may be classified as more than one of a frequency subtask, a model cell subtask, a path subtask, and a travel time subtask, i.e. being for one particular frequency subrange and / or one wave path, etc.

[0021] In some examples, the subtasks are recorded in a work queue, wherein each subtask has a position in the work queue, wherein a result of each completed subtask is recorded in a result queue at a result position corresponding to the position in the work queue for the completed subtask. For example, the positions in the work queue may be numbered sequentially and the corresponding result position is the position number in the result queue equal to the position number of the work queue. In another example, the result position may be determined by an algorithm based on the position in the work queue. By retaining a structure to the results that reflects the input work queue as described above, the collating values into the model can be particularly efficient, as it reduces or avoids time spent searching for the relevant values to record in the model.

[0022] In some examples, the performing tomographic inversion comprises: obtaining an initial model of having initial physical property values for the plurality of cells; for each receiver pair of the set, determining a modelled travel time using the initial physical property values associated with the cells of the model that are traversed by the surface wave path between receivers of the receiver pair; and determining new physical property values for the plurality of cells based on the modelled travel times and empirical travel times. In some examples, the initial model may be chosen based on a coarser model (i.e., having larger cells) determined previously using less precise methods, or based on approximate expected physical property values (e.g. based on the typical geological properties of the location), or a default value. The determining new physical property values may including comparing modelled travel times and empirical travel times and refining the model to reduce the difference between the modelled travel times and the empirical travel times, e.g. by using a least-squares method or Markov-chain Monte Carlo method.

[0023] According to a second aspect of the present disclosure, there is provided a method for determining a model of a physical property of a subsurface target volume. The method may be performed by a user device. The method comprises sending, to a network of computational devices, 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. The method comprises sending, to the network, instructions that cause the network to perform a method as described above in relation to the first aspect of the present disclosure and optionally any of the examples for doing so as described above. The method comprises receiving, from the network, the model of the physical property of the subsurface target volume. The sending may be from the user device to the network and / or the receiving may be at the user device from the network.

[0024] In some examples, the instructions comprise computation parameters including one or more of: dimensions of the subsurface target volume; dimensions of the plurality of cells of the model; a target precision of the physical property value; a target compute time; an amount of computation resource to use; a number of subtasks to create; and a type of subtask to create. The dimensions of the subsurface target volume may be the full height, width, and depth of the volume and may include the shape of the volume if it is not cuboid. The dimensions of the plurality of cells of the model may include how fine the model sought is, which controls the precision of the resulting model and also the amount of computation and time required (where a finer model, having smaller cells, will require more computation). The dimensions of the plurality of the cells may also be defined indirectly by the total number of cells, or the total number of columns and rows (which define the 2D surface grid above the target volume, wherein the value(s) of each cell in a given row and column includes how the value varies throughout the depth below that cell). The target compute time will be used, optionally in combination with the model size and cell size, to determine how much computation resource, e.g. amount of processing power, is required to determine the method in the target compute time. The amount of computation resource to use may directly define the number and power of processing units. The number and type of subtasks sets the parallelisation parameters in order to increase the speed and efficiency of the method. These can be defined directly in the instructions, or the instructions may provide a set of conditions to apply.

[0025] In some examples, the method comprises receiving, from the network, a request for further instructions due to the network reaching a programmed checkpoint. Correspondingly, the method described above in relation to the first aspect may comprise sending, to the user device, the request for further instructions due to the network reaching the programmed checkpoint. The further instructions may be provided by the user device based on information stored locally at the user device, or may be input by a user via a user interface at the user device. By using programmed checkpoints, the method can reduce wasted time proceeding based on incorrect assumptions, only to provide an incorrect result that then requires the user to reinitiate the whole method.

[0026] According to a third aspect of the present disclosure, a user device is provided. The user device comprises 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 operations comprising a method as described above according to the second aspect, and optionally any of the examples for doing so as described above.

[0027] According to a fourth aspect of the present disclosure, a computer-readable medium is provided. The computer-readable medium comprises instructions, that, when executed by the user device as described above according to the third aspect, cause the one or more processors to perform operations comprising a method according to the second aspect (and optionally any of the examples for doing so as described above). The computer-readable medium may be transitory or non-transitory.

[0028] According to a fifth aspect of the present disclosure, a system is provided, the system comprising a network of computational devices configured to perform a method as described above according to the first aspect, and optionally any of the examples for doing so as described above. In some examples, the system comprises the user device as described above according to the third aspect.

[0029] According to a sixth aspect of the present disclosure, a computer-readable medium is provided. The computer readable medium comprises instructions, that, when executed by the system as described above according to the fifth aspect (and optionally any of the examples for doing so as described above), cause the network to perform operations comprising a method according to the first aspect (and optionally any of the examples for doing so as described above). The computer readable medium may be transitory or non-transitory. BRIEF DESCRIPTION OF THE DRAWINGS 

[0030] Disclosed implementations will now be described by way of example to illustrate aspects of the disclosure and with reference to the accompanying drawings, in which:Figure 1 shows a cross-sectional view of a subsurface target volume;Figure 2 shows a plurality of receivers arranged on a surface above a subsurface target volume;Figure 3 is a schematic diagram of a shear wave velocity model;Figure 4A is an overhead view of a straight wave 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 wave 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 cloud-computing environment;Figure 6 is a schematic diagram of a method for determining a model of a physical property of a subsurface target volume;Figure 7 is a schematic diagram of a process used in the method of Figure 6;Figure 8 is a schematic diagram of a method for determining a model of a physical property of a subsurface target volume; Figure 9 is an example of a shear wave velocity model resulting from the methods described herein;Figure 10 is a schematic diagram of a computing device suitable for performing the methods described herein; andFigure 11 is a schematic diagram of a system suitable for performing the methods described herein. DETAILED DESCRIPTION OF THE DRAWINGS 

[0031] This detailed description describes, with reference to Figures 1 to 4, an approach to measuring structural properties of a ground volume using receivers. Next, with reference to Figures 5 to 9, methods are described for determining a model of a physical property of a subsurface target volume using a network of computational devices. Finally, a user device and a system that may be used to perform the disclosed methods are described with reference to Figures 10 and 11.

[0032] The following examples will be described in the context of a receiver array on a surface above a subsurface volume, 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.

[0033] The methods and systems disclosed herein relate generally to processing signals detected by receivers on a surface. In one example, the receivers are a arranged on a 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. Due to the division of the method into subtasks that are distributed between computational devices of the network and executed in parallel, the speed and efficiency of computation is increased, making processing even a very large number and duration of signals feasible within a practical length of time.

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

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

[0036] 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 and Scholte waves) may be measured and harnessed in the systems and methods described herein.

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

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

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

[0040] At an example point, a schematic representation of a particle oscillation (due to Rayleigh wave propagation) at a surface above a subsurface target 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.

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

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

[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 travel time measurement between a (virtual) source and a receiver as a function of frequency, i.e., a dispersion function.

[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 travel time 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 (time-distance 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 ambient 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] Figure 2 shows a plurality of virtual receiver pairs across a surface above a subsurface target volume. Wave paths 202 between receiver pairs are indicated, in particular, the wave 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 wave paths between each receiver and all other receivers, i.e. between every pair of receivers, which are not depicted in Figure 2 for simplicity.

[0051] 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 ambient 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.

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

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

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

[0055] With reference to Figure 3, a first shear wave velocity model 300 has a grid of cells mxy across the surface above the subsurface target volume, comprising columns mx1, mx2, mx3, etc. extending in the x-direction and rows m1y, m2y, m3y, 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. 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 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.

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

[0057] 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 velocity 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.

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

[0059] With reference to figures 4A to 4B, 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 wave path connecting the source and receiver. Alternatively, 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, the wave path may be represented by a curved Fresnel zone.

[0060] With reference to Figure 4A, a wave path is defined as a straight line (sometimes referred to as the ‘ray’ path according to ray theory) between two receivers at surface locations A and B on the model 300. A straight wave path model assumes a laterally constant velocity (in other words, the velocity only changes with depth). The straight wave 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 wave path approximation is also mathematically and computationally efficient compared to other more complex techniques. The wave path signifies the motion of a surface wave as it travels from A to B (or B to A) according to wave theory. In particular, a wave 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 wave path shown in Figure 4A traverses seven cells of the model 300, numbered 1 to 7. The wave path comprises seven wave segments, each having a length Li according to the length for which the wave path travels through each cell, i.e. L1 to L7 in Figure 4A.

[0061] With reference to Figure 4B, a wave path is defined as a curved line between two receivers at surface locations C and D. The wave path shown in Figure 4B traverses six cells of the model 300, numbered 1 to 6. The wave path comprises six wave segments, each having a length Li according to the length for which the wave path travels through each cell, i.e. L1 to L6 in Figure 4B. The trajectory of the wave 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.

[0062] With reference to Figure 5, the methods described herein may be performed in a cloud-computing environment 500. In some examples, the cloud-computing environment 500 has a virtual state machine 510 for coordinating the computing processes in the cloud-computing environment 500 as defined by a user, a virtual machine 520, a virtual storage 530, and an accessible storage 540. The virtual machine 520 emulates a processor and communicates with the virtual state machine 510 to receive instructions what to perform and to send updates on progress and results, e.g. a status report and / or the location of results data. The virtual storage 530 to emulate a memory for which the virtual machine 520 interacts with to perform its computations. The virtual storage 530 may comprise one or more virtual SSD (solid state drive) and / or one or more virtual HDD (hard disk drive) according to the requirements of the processes run by the virtual machine 520. The accessible storage 540 is an object storage where that can be accessed by a user device, e.g. over the internet via an API (application programming interface) gateway. A user of the cloud-computing environment 500 can perform the method described herein using the network of computational devices (such as described below with reference to Figure 11), on which the cloud-computing environment 500 is based, by defining the virtual state machine 510, choosing the number and type of virtual machines 520 and virtual storages 530, inputting the input data (i.e. the plurality of signals detected by the plurality of receivers above the subsurface target volume), and extracting the result data (i.e. the model of the physical property of the subsurface target volume).

[0063] With reference to Figure 6, a method 600 for determining one or more physical properties of a subsurface target volume includes receiving 602 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.

[0064] The signals may be received directly from the receivers or via one or more intermediary devices. In some examples, additional processing steps may be performed on the signals before or after the receiving 602 of the signals to optimise the signals for further processing.

[0065] Example receivers 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.

[0066] The method 600 also comprises cross-correlating 604 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 ambient 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. The cross-correlating 604 of the signals between the plurality of receivers is carried out to obtain group and / or phase travel time 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 travel time data may be obtained for a plurality of finite frequencies. Group travel time refers to the time taken for a group wave to travel from the virtual source to the receiver. Phase travel time 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 travel time are effectively equivalent, with a scaling factor provided by the distance between source and receiver.

[0067] Based on at least the group and / or phase travel time data for the plurality of frequencies, 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.

[0068] Determining the model of the physical property of the subsurface target volume comprises selecting 606 a set of the plurality of receiver pairs for which empirical travel time information has been obtained. As described above with reference to Figures 4A and 4B, a wave path between each receiver pair typically 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 longer 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 intended for use, 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 travel time 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 omit pairs 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 desired to be used, a subsample of receiver pairs may be chosen, e.g. every nth pair 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. 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 travel time 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.

[0069] As described above with reference to Figures 4A and 4B, a wave path between a (virtual) source and receiver pair may be a wave path that extends in a straight line between the source and receiver. In some examples, at least some of the wave 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 path segment length in each cell which it traverses (including the cells which the wave path starts and ends in). Curved wave paths may be determined according to initial physical property values for the cells of the model. By using curved wave paths determined according to the physical property values of the cells of the model, the wave paths more closely follow the actual wave 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.

[0070] As described above with reference to Figures 4A and 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 corresponds more closely to the physical reality of surface waves passing through the subsurface target volume and scattering within the volume.

[0071] Determining the model of the physical property of the subsurface target volume further comprises performing 608 tomographic inversion on the travel time data for each of the selected receiver pairs. The tomographic inversion is used to derive a value of the physical property for each cell of the model, from group and / 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. It is understood that any appropriate process of carrying out the tomographic inversion could be followed. The tomographic inversion can be carried out based on empirical group travel time information or empirical phase dispersion information. In the context of phase dispersion, phase velocity or slowness information may be converted to equivalent travel time information so that the tomography can be performed on such phase travel time information.

[0072] In some examples, performing the tomographic inversion comprises two stages, first a tomography stage and then an inversion stage. In some examples, the tomography stage is divided into frequency subtasks. Alternatively or additionally, the inversion stage is divided into cell subtasks.

[0073] The tomography stage tomograhpically maps the travel time 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 travel time data used), with each cell of the model having a respective group or phase velocity value. This process is carried out for each of a plurality of frequencies to obtain group or phase velocity values for each cell for each frequency. This is the first stage (i.e. the tomography stage) of the two-stage tomographic inversion.

[0074] The tomography process comprises 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 tomography process will refine using the empirical travel time 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 tomography process faster or more accurate at mapping travel times to each cell. A more accurate initial model may also reduce the risk of finding a local-minimum rather than a global minimum in the iterative gradient-descent method, although such an issue can 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.

[0075] The tomography process further comprises determining modelled group or phase travel times for each of the selected source-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.

[0076] The tomography process further comprises determining an error value indicative of the difference between the modelled travel times and the empirical travel times for each of the selected source-receiver pairs. For example, the error value could be a simple difference between the modelled and empirical travel times, also called a residual, at each frequency for each source-receiver pair. In some examples, the error value could be a combination of all the differences between modelled and empirical travel times for every source-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 velocity model.

[0077] The tomography process further comprises 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 travel times for each source-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. The tomography process may be repeated for each of a plurality of finite frequencies.

[0078] The determining the modelled travel time, the determining and error value and determining an updated velocity model are typically all part of a subroutine of the first stage tomography process, which is then iterated according to the inversion method such as least-squares inversion. The iterative nature of the tomography process is such that the updated model determined is used to determine new modelled travel times for each source-receiver pair. In other words, each time an updated model is determined using the error value resulting from the initial model and resulting modelled travel times for each source-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 travel times, or falling below a threshold of proportional difference between the modelled and empirical travel times. 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.

[0079] Once the first-stage tomography process 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 can proceed to the second-stage, i.e., the inversion process of the two-stage tomographic inversion to obtain the resultant model of the physical property of the subsurface target volume. As with the first-stage tomography process, the second-stage inversion process is carried out for each of a plurality of frequencies to obtain group or phase velocity values for each cell for each frequency.

[0080] The inversion process comprises obtaining an initial physical property model. The initial physical property 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 above with reference to Figure 3. The initial model sets initial physical property values for each cell of the model, which the inversion process will refine using the group / phase velocity model obtained from iterative tomography process described above, and using a further iterative process. Accordingly, it is not essential for the initial physical property 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.

[0081] In some examples, the initial physical property 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 initial physical model may be determined using the received signals, e.g. by performing an inversion of group velocity or phase velocity dispersion curves between receivers, 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., which may be based on estimated physical properties of the subsurface target volume. In some examples, the initial physical model may be determined based on empirical phase dispersion data between source-receiver pairs. Phase velocity information is acquired as a point in 2D phase space (a dispersion spectrum) which may be obtained by a 2D transform of an array of waveforms (time-distance space), the transform may be Radon, slant-slack, FK or any other suitable methods apparent to the skilled person. The generation of dispersion curves showing the phase velocity as a function of frequency may provide a preliminary three-dimensional model of the shear velocities in a 3D model. However, since the acquisition of the phase velocity is essentially an average frequency-dependent velocity between measurement positions, the resolution of this model is relatively low. Thus, a physical model derived solely from empirical phase dispersion data may be used as an initial model for an inversion which also takes into account empirical group dispersion data to obtain a more accurate final 3D model of the subsurface target volume.

[0082] The inversion process further comprises determining modelled surface wave velocities for each cell based on the initial model. 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).

[0083] The inversion process further comprises determining an error value indicative of the difference between the modelled velocities and the empirical velocities for each cell of the model. The error value may be determined for each cell in a similar manner to that described above in relation to the tomography process in which an error value is determined between modelled and empirical travel times for each source-receiver pair. As with the tomography process, 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 model of the physical properties of the subsurface target volume.

[0084] The inversion process further comprises determining an updated physical model based on the determined error value. The updated physical property 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 for each cell of the model. 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.

[0085] The determining modelled surface wave velocities, determining the error value, and determining the updated physical property model are typically all part of a subroutine of the second-stage inversion process, which is then iterated according to the inversion method such as least-squares inversion. The iterative nature of the inversion process is such that the updated model determined is used to determine new surface wave velocities for each cell. 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 the tomography process.

[0086] In an alternative to the two-stage tomographic inversion, the performing 608 tomographic inversion may be a single-stage process in which the tomography and inversion are combined. This means that the iterative inversion stage to find the physical property values for the model is carried out for each wave path between a respective source-receiver pair, rather than for each cell of the model. This means that the first-stage tomography process for mapping surface wave velocities between source-receiver pairs to cells is not required in the single-stage process because the physical property model inversion process carries out the inversion on each wave path, rather than each cell. In other words, the two-stage tomographic inversion first uses travel time tomography to map travel time information along wave paths into surface wave velocity information for a grid of cells. The two-stage tomographic inversion then carries out a cell-by-cell inversion of the surface wave velocity information (surface wave velocity as a function of frequency) to obtain a corresponding shear wave velocity function (as a function of depth) for each cell. In contrast, the one-stage process carries out a direct inversion of the surface wave velocity information for each wave path (for each source-receiver pair) to obtain a shear wave velocity function for the wave path, which is tomographically mapped to each cell of the model.

[0087] One or more of the cross-correlating, the selecting, the performing tomographic inversion, and the collating, as described above with reference to Figure 6, involves a parallel processing routine in a cloud-computing environment 500, such as described above with reference to Figure 5. In some examples, the parallel processing routine divides part of the method into subtasks according to frequency, as described below with reference to Figure 7.

[0088] With reference to Figure 7, a parallel processing routine 700 based on frequency comprises determining 702 a plurality of frequency subranges f1, f2, f3, … fN within the received (602 in Figure 6, as described above) plurality of signals. The subranges may be provided by a user device at the same time as receiving 602 the plurality of signals, e.g., in metadata accompanying the plurality of signals and chosen by a user sending the plurality of signals to the cloud-computing environment 500. In some examples, the virtual state machine 502 has been programmed with a protocol to determine the frequency subranges from the plurality of signals, e.g. dividing the total frequency range of the plurality of signals into N (e.g. 100 or some other predetermined number) subranges of equal bandwidth, or subranges on a logarithmic scale. In some examples, the number and bandwidth of the subranges is determined based on a desired compute time and / or a desired amount computational resource to be used.

[0089] The frequency parallel processing routine 700 divides part of the method into subtasks T1, T2, T3, … TN, where each subtask is associated with one of the plurality of frequency subranges f1, f2, f3, … fN. The subtasks T are recorded in a work queue 704, each subtasks position in the work queue 704 is based on its subtask number. The work queue 704 may be 1-dimensional queue or an array with two or more dimensions.

[0090] Each subtask includes any information required to independently carry out, for its designated frequency subrange, the part of the method being divided into subtasks. In some examples, the part of the method being divided into subtasks includes all of the cross-correlating, the selecting, the performing tomographic inversion, and the collating. In this instance, these parts of the method can be performed for a particular frequency subrange in the same manner as described above with reference to Figure 6. In some examples, the subtask includes information regarding where the result of the subtask should be recorded in the resulting physical property model 300 during the collating 610. Since, as explained above, the surface waves are influenced by physical properties of the medium within a depth range of approximately one wavelength from the surface, the result of each frequency subtask provides the physical property values for different depths in the model (e.g. comparing the results of two frequency subranges indicates the effect of the layer in the subsurface target volume that only significantly influences the lower frequency subrange).

[0091] The virtual machine 520 of the cloud-computing environment 500 can perform multiple processes at the same time (using virtual storage 530 for any short term storage during the processing), since it is created by a distributed network of computational devices. This can be understood as a worker pool 706 having a number of workers W1, W2, W3, … WP, with P representing the number of independent processes that the virtual machine 520 can perform concurrently. The number of workers P may be determined by the virtual state machine 510 at the start of the method or may be determined by a user input according to how much computational resources the user wishes to use. Each worker W in the worker pool 706 picks up a subtask T and executes the actions associated with that subtask and records the result R into a result queue 708. Accordingly, the compute time is reduced compared to conventional techniques, by leveraging the large processing resources available in the cloud-computing environment across a network of computational devices. Once all the subtasks T are complete, the result queue 708 will be complete with results R1, R2, … RN. The workers W may select subtasks T to perform by each free worker selecting the next available subtask in the work queue 704 (the subtask with the lowest subtask number), or according to some other protocol, or even at random. When a worker has completed its subtask, it will select a new subtask from the work queue until all the subtasks have been completed. The results R are recorded into the result queue 708 at a position corresponding to the subtask number, so that the order and structure of the subtasks in the work queue 704 is preserved in the result queue 708. This means that it will be quicker and simpler for the results to be extracted from the result queue 708. In some examples, each result R may comprise multiple values that will be used to determine the model of the physical property.

[0092] In some examples, a worker W may split its subtask further into smaller divisions, for example, if each worker corresponds to a virtual CPU with four cores, four processes can be performed simultaneously.

[0093] In some examples, the frequency parallel processing routine 700 includes combining 710 the results R1, R2, RN from the work queue into a single result. This may be part of the collating 610 of the overall method 600.

[0094] Although the above example of a parallel processing routine is a frequency parallel processing routine 700, the same principles apply to creating subtasks according to model cells, wave paths, travel times, or other parameters of the method 600 for determining the physical model.

[0095] For example, when creating subtasks based on model cells, the number of subtasks will correspond to the number of cells in the model to be determined. For example, this could on the order of thousands, e.g. 5000. By dividing the process into subtasks that are executed in parallel by the workers W, the time it takes to determine the model of physical property values is significantly reduced. For example, in an example having approximately 500 workers, the computation time will be on the order of 500 times faster. In some examples, each model cell subtasks include selecting pairs of receivers, and the corresponding cross-correlated signals, which have surface wave paths traversing the designated model cell for that subtask. The travel time information may then be combined to produce the group and / or phase velocity information for that cell and an inversional calculation performed, as described above with reference to Figure 6, to obtain the respective physical property value for that cell. The result queue 708 after all model cell subtasks have been completed will then contain the final information for the model of the physical value of the subsurface target value.

[0096] In some examples, when dividing into subtasks based on surface wave paths, each path subtasks may be designated to one pair of receivers. This may be before the cross-correlating 604, so that each cross-correlation is done by a separate path subtask, or after the selecting 606 only for the selected pairs of receivers (and their corresponding surface wave paths). Path subtasks are well suited for the one-stage tomographic inversion process described above, as the inversion calculation is done for an entire surface wave path, and the calculated value then applied to each of the model cells that the surface wave path traverses (optionally with weightings according the proportion of the surface wave path in each model cell, Li as shown in Figures 4A and 4B).

[0097] In general, the subtasks into which one or more parts of the method are divided may be any set of instructions that, when executed, contribute to the performance of the method for determining a model of a subsurface target volume. The subtasks may include the relevant data (e.g. signal, part of a signal, or intermediate result of the method) or a location from where to retrieve the data to be processed in that subtask. The subtask may also include a label to identify the subtask and how it should be processed.

[0098] With reference to Figure 8, a method 800 for determining a model of a physical property of a subsurface target volume can be performed by a user device interacting with a cloud-computing environment 500 created by a network of computational devices. The method comprises sending a plurality of signals to the network of computational devices. In some examples, this is through an API gateway or saving the plurality of signals in an accessible storage 540 in the cloud-computing environment 500.

[0099] The method comprises sending instructions to the network instructions to perform a method as described above with reference to Figure 6. This may be done by an orchestrator program that manages HTTP requests to and from the cloud-computing environment 500. The instructions may comprise a file address where to locate the plurality signals for performing the method, if already saved on the accessible storage 540. The instructions may also comprise the specifics of the model to be determined (dimensions and cell size) and / or an initial model from where to start the tomographic inversion. The instructions may also include user input defining the level of accuracy desired (e.g. how small a discrepancy between the modelled and empirical travel times needed to end the inversion process) the amount and type of computational resource to use (e.g. the number of virtual processors for the virtual machine 520 to run, the numbers of subtasks T to use, and / or the number of workers W to use) or a desired compute time from which the virtual state machine 510 can calculate the necessary amount of computational resources to deploy in order to achieve that compute time.

[00100] The method further comprises receiving the determined model from the network, in response to the computation devices of the network performing the method as described above with reference to Figure 6 in the cloud-computing environment 500. Accordingly, the method on the user device controls the parallelisation parameters in the cloud-computing environment 500 to achieve a fast and efficient determination of the model.

[00101] In some examples, the user device works through a virtual private cloud for providing secure traffic to and from the network of computational devices.

[00102] In some examples, the virtual state machine 510 in the cloud-computing environment 500 has pre-programmed checkpoints for user input. This can be useful to test the parallelisation methods described herein, because a single part of the overall process can be performed and reviewed without the rest. In some examples, the checkpoint is caused by identification of an unexpected event in the procedure, at which point the virtual state machine 510 can halt the process and request further instruction from the user device how to proceed.

[00103] In some examples, the method comprises pre-processing the plurality of signals into a format to be used by the network, or the instructions include how to pre-process the plurality of signals.

[00104] 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 high resolution and accuracy and a shorter 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.

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

[00106] 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 / S-wave) velocity Vs are related to other physical properties of the ground elastic modulus λ, shear modulus μ, and density ρ by the following equations from linear elasticity theory:Equation 1 (longitudinal wave velocity):  Equation 2 (shear wave velocity):  

[00107] 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 a method as described above with reference to Figure 8 as the user device, 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.

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

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

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

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

[00112] The data storage device 1018 may include one or more computer-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.

[00113] The methods described above with reference to Figure 8 may be implemented by a computer program. The computer program may include computer code arranged to instruct a computer to perform the functions of the methods. 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.

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

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

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

[00117] 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 computer-readable medium or in a transmission medium).

[00118] Figure 11 shows a block diagram, of one implementation of a system 1100 within which a set of instructions, for causing a network of computational devices to perform a method as described above as the user device, may be executed. For example, the network comprises a plurality of computational devices 1110 to provide the cloud-computing environment 500 as described above with reference to Figure 5, according to any standard cloud-computing method. The computational devices 1110 may have specific computer-readable storage media 1128 having any of the features of the computer-readable storage media 1028 as described above with reference to Figure 10.

[00119] The computational devices 1110 forming the network are connected by a communication connection 1130, which may be a local connection (e.g. via a Local Area Network (LAN)), an intranet, an extranet, a virtual private cloud, or the Internet, according to any typical cloud-computing arrangement.

[00120] The system may also comprise a physical storage device 1120 for storing data and accessible to the network of computational devices via the communication connection 1130.

[00121] The system 1100 may also comprise a computational device 1000 as described above with reference Figure 10, to perform as a user device according to a method as described above with reference to Figure 8, and computer-readable storage media 1028 as described above with reference to Figure 8.

[00122] In some examples, the connection between the network of computational devices 1110 is separate to the connections with the user device (computational device 1000) and / or the physical storage device 1120.

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

[00124] 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 entitle.

Claims

1. A method for determining a model of a physical property of a subsurface target volume using a network of computational devices, wherein the model comprises a plurality of cells and each cell has a respective physical property value, 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; selecting a set of the plurality of receiver pairs, each respective receiver pair having a respective surface wave path between receivers of the receiver pair; performing tomographic inversion on the empirical travel time data between receivers of each receiver pair in the set to obtain a respective physical property value for each cell; andcollating the obtained physical property values into the model;wherein one or more of the cross-correlating, the selecting, the performing tomographic inversion, and the collating is divided into subtasks, wherein the subtasks are distributed between the computational devices of the network and executed in parallel.

2. A method according to claim 1, wherein the subtasks are frequency subtasks, wherein each frequency subtask corresponds to a respective frequency subrange of a total frequency range of the plurality of signals.

3. A method according to any preceding claim, wherein the subtasks are model cell subtasks, wherein each model cell subtask corresponds to a respective cell of the model.

4. A method according to any preceding claim, wherein the subtasks are path subtasks, wherein each path subtask corresponds to a respective surface wave path.

5. A method according to any preceding claim, wherein the subtasks are travel time subtasks, wherein each travel time subtask corresponds to respective obtained empirical travel time data.

6. A method according to any preceding claim, wherein the subtasks are recorded in a work queue, wherein each subtask has a position in the work queue, wherein a result of each completed subtask is recorded in a result queue at a result position corresponding to the position in the work queue for the completed subtask.

7. A method according to any preceding claim, further comprising recording subtask status and work performed in a central database via an applications programming interface.

8. A method according to any preceding claim, further comprising using a central service operating a daemon observing subtask and work status to take appropriate action to perform new compute tasks when conditions are met.

9. A method according to any preceding claim, wherein the performing tomographic inversion comprises: obtaining an initial model of having initial physical property values for the plurality of cells; for each receiver pair of the set, determining a modelled travel time using the initial physical property values associated with the cells of the model that are traversed by the surface wave path between receivers of the receiver pair; anddetermining new physical property values for the plurality of cells based on the modelled travel times and empirical travel times.

10. A method for determining a model of a physical property of a subsurface target volume, the method comprising:sending, to a network of computational devices, 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; sending, to the network, instructions that cause the network to perform a method according to any preceding claim; andreceiving, from the network, the model of the physical property of the subsurface target volume.

11. A method according to claim 10, wherein the instructions comprise computation parameters including one or more of: dimensions of the subsurface target volume; dimensions of the plurality of cells of the model; a target precision of the physical property value; a target compute time; an amount of computation resource to use; a number of subtasks to create; a type of subtask to create.

12. A method according to claim 10 or 11, wherein the method comprises:receiving, from the network, a request for further instructions due to the network reaching a programmed checkpoint.

13. A user device comprising: one or more processors; andone or more memories having stored thereon computer readable instructions configured to cause the one or more processors to perform operations comprising a method according to any of claims 10 to 12.

14. A computer readable medium comprising instructions, that, when executed by the user device according to claim 14, cause the one or more processors to perform operations comprising a method according to any of claims 10 to 12.

15. A system comprising:a network of computational devices configured to perform a method of any of claims 1 to 9.

16. A system according to claim 15, comprising:a user device according to claim 13.

17. A computer readable medium comprising instructions, that, when executed by a system according to claim 15 or 16, cause the network to perform operations comprising a method according to any of claims 1 to 9.   

Images