Chromatographic inversion

Abstract

A method for determining a subterranean target volumetric physical property, and a system and computer readable medium for performing the method are provided. The method includes receiving a plurality of signals detected by a plurality of receivers disposed on a surface above a subterranean target volume, where each respective signal of the plurality of signals is detected by a respective receiver of the plurality of receivers, cross-correlating the signals between the plurality of receivers, and transmitting the cross-correlated signals to the subterranean target volume. The method includes determining a subsurface target volume to obtain empirical travel time data of surface waves between a plurality of virtual source-receiver pairs, determining a model of a physical property of the subsurface target volume, the model including a first plurality of cells arranged in a first two-dimensional grid, where each cell has a respective physical property value, and where each cell has a respective physical property value. The two-dimensional grid at least crosses an area of the earth surface. Unlocking from the perspective of geographic data, the invention further relates to improvements in sustainability and environmental development: we together create a safe and livable world with their hands.

Description

Technical Field

[0001] The present disclosure relates to methods and systems for analyzing a target area below the surface of the earth. More specifically, the present disclosure relates to methods and systems for determining one or more material properties of a subsurface target area based on noise measured at or near the surface. In particular, the present disclosure provides methods and systems for providing a model of a subsurface target area with improved accuracy. Unlocking insights from geodata, the present invention further relates to improvements in sustainability and environmental development: Together, we create a safe and livable world. Background Art

[0002] There is a general and ongoing need for systems and methods for determining subsurface formation parameters. In particular, there is a need for systems and methods that can be used to model the properties of a target volume below the surface of the earth to provide information useful for infrastructure planning. Determining subsurface formation properties during the early planning stages of a construction project reduces uncertainty during the project's site selection, foundation design, and construction phases. This in turn reduces delays, cost overruns, and unnecessary use of material resources (such as concrete) during construction.

[0003] A key parameter for determining formation properties in a volume of interest is the shear modulus and the shear velocity Vs. The shear velocity Vs is the velocity at which a shear wave travels through a material and is controlled by the shear modulus of the material. The relationship between the shear velocity and the shear modulus G is defined as Vs = √G / ρ, where ρ is the density of the material. Thus, measurement of Vs provides valuable insights into the material properties of a subsurface formation area.

[0004] Spectral analysis of surface waves (SASW) and multichannel analysis of surface waves (MASW) are both examples of techniques for acquiring surface wave information that can be used to determine the material properties of a subsurface volume. In both of these techniques, surface-level vibrations are measured, which originate from passive sources (surface vibrations caused by environmental noise sources) or active sources (such as a heavy object drop), and the resulting dispersion of the surface waves is studied. ReMi (Refraction Microtremor) is another surface-level technique that uses environmental noise and surface waves to infer the material properties of a subsurface area based on observations of environmental noise at the surface.

[0005] Downhole and crosshole techniques can also be used to determine the material properties of a subsurface region. In both methods, a receiver located within a borehole measures waves received from an active source located elsewhere. In downhole techniques, one of the source and the receiver is located at a subsurface location within the borehole, and the other of the source and the receiver is located at the surface. In crosshole techniques, the source is located in a first borehole and the receiver is located in a second borehole. In both downhole techniques, the propagation and dispersion of the received waves are studied to infer the properties of the material through which the waves received by the receiver have traveled.

[0006] Intrusive techniques for measuring the material properties of a subsurface region can generally pose logistical challenges, especially in urban or inaccessible environments, and are generally extremely expensive. In contrast, current surface-level techniques may lack the accuracy and reliability of more intrusive analytical techniques. Summary of the Invention

[0007] According to a first aspect of the present disclosure, a method for determining physical properties of an underground target volume is provided. The method includes: receiving a plurality of signals detected by a plurality of receivers disposed on the ground above the underground target volume, wherein each respective signal of the plurality of signals is detected by a respective receiver of the plurality of receivers; performing cross-correlation on the signals among the plurality of receivers to obtain empirical travel time data of surface waves between a plurality of virtual source-receiver pairs; and determining a model of physical properties of the underground target volume, the model including a first plurality of cells arranged in a first two-dimensional grid, wherein each cell has a respective physical property value, and wherein the two-dimensional grid at least spans a region of the surface of the earth. Although reference is made to virtual source-receiver pairs, the positions of these virtual source-receiver pairs may coincide with the physical arrangement of the receivers. That is, when ambient waves travel through the target volume and are recorded by the plurality of receivers, the wave traveling between two receivers can be considered as traveling from a virtual source (i.e., the first receiver) to a receiver (i.e., the second receiver). Thus, the virtual source-receiver pairs may spatially coincide with the positions of the receivers placed on the surface of the target volume. As used herein, "empirical" data or information refers to such data or information obtained by experience, i.e., obtained from real signals collected at real receivers (such as geophones or accelerometers) on the surface of the earth above the underground target volume. Determining the model includes selecting a subset of the plurality of source-receiver pairs, wherein the surface wave path between each source and the respective receiver passes through two or more of the first plurality of cells, and performing tomographic inversion on the empirical travel time data between the source-receiver pairs in the subset to obtain the physical property values of each cell. By performing tomographic inversion on the empirical travel time data of the plurality of source-receiver pairs in this way, a three-dimensional model of the physical properties of the underground target volume can be obtained. This is achieved by obtaining a one-dimensional model for each cell in the cell grid, and these cells together provide information on the three-dimensional physical properties. This particular method also takes into account the position dependence between points, thereby generating a smoother and more realistic model, i.e., (relatively) avoiding physically impossible mutations in geophysical properties. This method also allows for the generation of a higher resolution model without increasing the computational time to an infeasible level.

[0008] In some examples, performing tomographic inversion includes: obtaining an initial velocity model that includes initial values indicative of the velocity of each of a first plurality of cells; for each source-receiver pair in a subset, determining a modeled travel time using initial first velocity model values associated with two or more cells through which a surface wave path passes; determining a first error value indicative of the difference between the modeled travel time and the corresponding empirical travel time for each source-receiver pair in the subset; and determining an updated velocity model based on the error value, where the updated velocity model includes updated values indicative of the velocity of each of the first plurality of cells. In some examples, determining the modeled travel time, determining the first error value, and determining the updated velocity model based on the error value are iterated until the error value meets a first end condition. The updated velocity model of each iteration is used as the initial velocity model for the next iteration. This provides an iterative process that allows for obtaining a high-accuracy velocity model (the velocity values of each cell in the model). In some examples, determining the modeled travel time, determining the first error value, and determining the updated velocity model based on the error value are performed using iterative non-linear least squares inversion, which is a computationally efficient inversion technique. This is feasible because the determination of the updated velocity is provided for each wave path and has a lower computational intensity compared to common alternative methods (such as the Monte Carlo method).

[0009] In some examples, performing tomographic inversion further includes: for each of the first plurality of cells, performing an inversion of the value indicative of velocity to obtain a value of the physical property of the cell. The inversion may include: obtaining an initial model of the physical properties of the subsurface volume, where the initial model includes initial physical property values of each of the first plurality of cells; for each of the first plurality of cells, determining a modeled value indicative of velocity based on the initial physical property value of the cell; determining a second error value that indicates the difference between the modeled value and the corresponding value, where the modeled value indicates velocity and the corresponding value indicates the velocity obtained from the updated velocity model; and determining an updated model of the physical properties of the subsurface volume based on the second error value, where the updated model of the physical properties of the subsurface volume includes updated physical property values of each of the first plurality of cells. Determining the modeled value indicative of velocity, determining the second error value, and determining the updated model of the physical properties of the subsurface volume based on the second error value are iterated until the second error value meets a second end condition, and where the updated model of the physical properties of the subsurface volume of each iteration is used as the initial for the physical properties of the subsurface volume of the next iteration. Determining the modeled value indicative of velocity, determining the second error value, and determining the updated model of the physical properties of the subsurface volume based on the second error value may be implemented using computationally efficient least squares inversion. This is feasible because the determination of the dispersion function is performed for each wave path and has a lower computational intensity compared to common alternative methods (such as the Monte Carlo method).

[0010] In an alternative example, performing tomographic inversion includes: obtaining an initial model of the physical properties of a subsurface volume, where the initial model includes initial physical property values for each of a first plurality of cells; determining an empirical dispersion function for each source-receiver pair of a subset using empirical travel time data; determining a modeled dispersion function for each source-receiver pair of the subset using the initial physical property values associated with two or more cells through which a surface wave travels from the source to the receiver; determining a third error value indicative of the difference between the modeled dispersion function and the empirical dispersion function for each source-receiver pair; and determining an updated model of the physical properties of the subsurface volume, where the updated model of the physical properties of the subsurface volume includes updated physical property values for each of the first plurality of cells. By performing the determination of the characteristics of the modeled dispersion function and comparing it with the empirical dispersion function (through error waves) of a portion of the source-receiver pairs (rather than comparing each cell individually), the accuracy of determining one or more physical properties of the subsurface target volume is improved. This is because updating the first model to produce the second model is performed as a single process for a given wave path, rather than using multiple wave paths in an initial tomographic process to find the dispersion function of a single cell, which can introduce additional approximations or inaccuracies.

[0011] Determining the modeled dispersion function for each source-receiver pair can include: averaging the physical property values of two or more cells through which a surface wave travels from the source to the receiver and using the averaged physical property values to calculate the modeled dispersion function. Alternatively, determining the modeled dispersion function for each source-receiver pair can include: calculating a cell dispersion function for each of two or more cells through which a surface wave travels from the source to the receiver to provide a plurality of cell dispersion functions, where calculating each cell dispersion function includes using the physical property value of the corresponding cell and averaging the plurality of cell dispersion functions. In some examples, the path taken by the surface wave between each source and receiver pair is represented by a straight or curved ray. In these examples, averaging can include: determining a weighted average according to the respective weights of each of the two or more cells, where the respective weights of each cell correspond to a segment length of the ray path in that cell. In other examples, the path taken by the surface wave between each source and receiver pair is represented by a Fresnel zone. In these examples, averaging includes determining a weighted average according to the respective weights of each of the two or more cells, where the respective weights of each cell correspond to the sensitivity value of the Fresnel zone in that cell.

[0012] By averaging over the cells in each wave path, the model is more accurate because it more closely reflects how the wave path between a pair of receivers behaves. This process also provides a viable way to achieve a high-resolution model, which has smaller cell sizes and thus a more accurate resulting model. The weighted averaging method also improves accuracy by considering the different effects of different cells on each wave path.

[0013] In some examples, determining an updated model of the physical properties of the subsurface includes updating the initial physical property values of each cell according to the corresponding weights of the paths taken by the surface waves in that cell. For example, the changes in the physical property values of each cell can be weighted as described above. These techniques mean that the resulting model more closely corresponds to the actual physical properties of the subsurface target volume and can also reduce the time or iterations required to find the final model that fits the detected signal.

[0014] In some examples, determining a second model includes: determining the partial derivatives of a third error value with respect to one or more physical properties of a first model; determining a sensitivity matrix that includes the partial derivatives; solving a system of linear equations defined by the sensitivity matrix to determine the desired changes to one or more initial physical property values of an initial model; and updating the initial model according to the desired changes to produce an updated model of the physical properties of the subsurface volume. This provides an accurate way to implement changes to the initial model to produce an updated model, which results in an updated model that more closely maps to the actual physical properties of the subsurface target volume.

[0015] In some examples, determining the modeled dispersion function, determining the third error value, and determining an updated model of the physical properties of the subsurface volume are iterated until the third error value meets a third end condition, where the updated model of each iteration is used as the initial model for the next iteration. The end condition can include one or more of the following: the third error value of the most recent iteration is less than a predetermined threshold error value; and the difference between the third error value of the most recent iteration and the third error value of the previous iteration is less than a predetermined difference threshold. This further improves the accuracy of the final updated model because each iteration brings the model into closer registration with the detected signal (and in particular the dispersion function determined therefrom). The end condition controls the balance between optimizing the second model to the actual physical properties of the subsurface target volume and reducing the computational time.

[0016] In some examples, determining the modeled dispersion function, determining the third error value, and determining the updated physical model are performed using non-linear iterative least squares inversion. This is possible because the determination of the dispersion function is performed for each wave path and is less computationally intensive compared to common alternative methods such as the Monte Carlo method.

[0017] In some examples, an initial model of the physical properties of a subsurface volume provides initial physical property values in at least two spatial dimensions and optionally in three spatial dimensions. This provides a more accurate and geologically more realistic initial model, which in turn enables a more accurate and geologically more realistic final physical property model to be obtained. A more accurate initial model in two or more spatial dimensions can also achieve a reduction in the number of iterations before an end condition is reached, thereby making the overall process of obtaining a model of the subsurface target volume more efficient. In some examples, the initial model of the physical properties of the subsurface volume is sourced from empirical phase dispersion data obtained from signals detected between a second subset of the plurality of source-receiver pairs. This is useful because a preliminary initial model can be obtained based on the empirical phase velocity information obtained for the plurality of source-receiver pairs. This means that the initial model is empirically based on the subsurface target volume itself to obtain a more accurate initial model for any given subsurface target volume.

[0018] In some examples, the initial physical property values and updated physical properties for each cell are shear wave velocities as a function of depth. This is a very useful physical property for the model because it has a particularly significant effect on surface wave propagation. In some embodiments, the method can include calculating compressional wave velocity values and / or density values for one or more cells in a second model using the shear wave velocity values. That is, any combination of shear wave velocity, compressional wave velocity values, and density values can be included in the model.

[0019] In some examples, the empirical dispersion function and the modeled dispersion function include a group velocity dispersion function and / or a phase velocity dispersion function. This is useful because group or phase dispersion data can be used to obtain a model of the subsurface target volume. Whether group or phase data can be used depends on the availability of the two types of data, the quality of the two types of data, or other factors that may mean that one of the group or phase data is more desirable than the other. Alternatively, both phase and group data can also be used in the tomographic inversion process described herein to produce a more accurate model of the physical properties of the subsurface target volume.

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

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

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

[0023] The disclosed embodiments will now be described by way of example with reference to the accompanying drawings to illustrate aspects of the present disclosure, wherein:

[0024] Figure 1 shows a cross-sectional view of a subsurface target volume;

[0025] Figure 2 shows a plurality of geophones arranged on a surface above the subsurface target volume;

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

[0027] Figure 4A is a top view of ray paths on a surface above the subsurface target volume, with a shear wave velocity model superimposed thereon;

[0028] Figure 4B is a top view of ray paths on a surface above the subsurface target volume, with a shear wave velocity model superimposed thereon;

[0029] Figure 4C is a top view of an elliptical Fresnel kernel representing surface wave paths above the subsurface target volume, with a shear wave velocity model superimposed thereon

[0030] Figure 4D is a top view of a curved Fresnel kernel representing surface wave paths above the subsurface target volume, with a shear wave velocity model superimposed thereon;

[0031] Figure 5 is a schematic diagram of a method for determining one or more physical properties of a subsurface target volume;

[0032] Figure 6 is a schematic diagram of a method for determining a surface wave group velocity model;

[0033] Figure 7 is a schematic diagram of a method for determining one or more physical properties of a subsurface target volume;

[0034] Figure 8 is a schematic diagram of a method for determining one or more physical properties of a subsurface target volume

[0035] Figure 9 is an example of a shear wave velocity model obtained from the method described herein; and

[0036] Figure 10 is a schematic diagram of a computing device suitable for performing the methods described herein. DETAILED DESCRIPTION

[0037] This detailed description refers toFigure 1 and Figure 2 describes a method for measuring the structural properties of a ground volume using a receiver, which in the case of the present exemplary embodiment is a geophone. Next, reference is made to Figures 3 to 9 discloses a new method for tomographic inversion using signals detected by a geophone. Finally, reference is made to Figure 10 describes a computing device that can be used to perform the disclosed method.

[0038] For purposes of aiding understanding, the following example will be described in the context of an array of geophones. However, it will be understood that the disclosed systems and methods are applicable to various types of receivers, including but not limited to geophones, accelerometers, velocimeters, seismographs, vibration sensors, and / or transducers. The disclosed methods can be applied to any suitable set of signals.

[0039] The methods and systems disclosed herein generally relate to processing signals detected by geophones on the ground. In one particular example, the geophones are placed on the ground and the detected signals are ambient noise signals. Processing these signals provides useful insights into the structure of the ground and subsurface target volumes where the geophones are deployed, as described in more detail below. Since the number of signals to be processed can be very large, higher accuracy and / or resolution of information about the subsurface target volume can be achieved. However, using a large number of signals also means that computationally efficient methods are useful for ensuring that these signals can be processed within a practical time length. The methods and systems disclosed herein provide a way to perform the processing in a computationally efficient and practical manner to generate accurate, high-resolution information about the subsurface target volume.

[0040] Before delving into the details of the disclosed method, some background knowledge related to using geophones to determine surface and subsurface properties is first introduced. Shear modulus is a measure of the elastic shear stiffness of a material and represents the deformation that occurs in a solid when it is subjected to a force parallel to one of its surfaces and an opposing force on its opposite surface. Such forces and their effects on subsurface target volumes are important parameters to 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 located above and / or passing through the volume.

[0041] In the case of ground studies, two types of waves are typically distinguished: P-waves, where the particles in the volume vibrate in the direction of wave travel, such that the wave causes compression and decompression of the formation as it propagates through the formation. S-waves are shear waves, where the particles vibrate in a direction perpendicular to the direction of wave propagation.

[0042] P-waves and S-waves are body waves and travel in all directions within the volume of the body. The interaction of P-waves and S-waves with the Earth's surface generates surface waves that travel along the Earth's surface. Surface waves can be divided into several types. In the systems and methods described herein, Rayleigh waves are measured and studied because they facilitate the measurement of the vertical component of surface vibrations. However, it will be understood that other surface waves (such as Love waves) can be measured and utilized in the systems and methods described herein.

[0043] Since surface waves travel in two dimensions (at the Earth's surface), they attenuate more slowly than body waves (which travel in three dimensions). Surface waves typically exist within a depth range of one wavelength from the Earth's surface, travel generally more slowly than body waves, and have significantly lower frequencies than body waves.

[0044] This lower attenuation, slower travel time, and lower frequency of surface waves make their study particularly attractive for the purpose of determining the shear velocity Vs. Since surface waves have lower attenuation, the signal strength is better maintained over longer travel distances. As a result, the resulting measurements typically have a higher signal quality (signal-to-noise ratio) than body wave studies.

[0045] Now referring Figure 1 , a cross-sectional view of a subsurface volume 100 is shown. P-waves and S-waves travel through the volume 100 as body waves. The Earth's surface 102 extends above the subsurface volume. Surface waves travel along the Earth's surface 102.

[0046] At an example point, a schematic representation of particle vibrations (caused by the propagation of Rayleigh waves) at the Earth's surface above the target subsurface volume is shown. As shown, the vibrations of particle P are partially perpendicular to the propagation direction and partially parallel to the propagation direction. As a result, the resulting particle motion is generally elliptical.

[0047] A plurality of geophones 104 are arranged on the Earth's surface 102 above the volume 100. The geophones 104 located on the Earth's surface 102 can be configured to measure the vertical component of the vibrations schematically shown at the example point.

[0048] The geophones 104 are arranged in a grid array on the Earth's surface, which extends in two directions. It should be noted that the Earth's surface above the target area may not be planar in many cases. Therefore, the geophone 104 array may not be truly "two-dimensional" because each geophone may be offset in the z-direction from its adjacent geophones in the grid. However, for simplicity, this grid arrangement of geophones will be referred to herein as a two-dimensional array.

[0049] It will also be understood that in the context of subsurface formation studies, the term "two-dimensional array" may also be used to denote an array configured to capture two-dimensional information (e.g., a row of receivers configured to provide information about a two-dimensional slice of a target area), while a "three-dimensional array" may refer to a grid of receivers configured to capture three-dimensional information. However, in the context of the current information, the term "two-dimensional" in the context of receivers is intended to refer to the two-dimensional configuration of the receivers, rather than the information collected.

[0050] When the surface wave passes through the surface 102, the surface wave traveling through the surface 102 will cause vertical motion at a plurality of geophones 104.

[0051] In order to determine the shear velocity Vs by observing surface waves (especially Rayleigh waves), the dispersion behavior of the surface waves can be measured. Surface waves are dispersive, which means that their velocity depends on the frequency. Generally, seismic velocity increases with depth in the Earth. Therefore, normal surface wave dispersion will exhibit a decrease in surface wave velocity with increasing frequency. By studying the behavior of surface waves at the surface above the subsurface target volume, the material properties of that volume can be determined.

[0052] There are two methods for measuring the velocity in dispersive surface waves, and a distinction is made between determining the group velocity or the phase velocity.

[0053] The group velocity of a wave is the velocity at which the overall envelope shape of the wave amplitude - called the modulation or envelope of the wave - propagates through space. The group velocity is equivalent to the propagation of the wave energy through a volume of space. The group velocity is measured by determining the wave propagation between (synthetic) transducer pairs, and the group velocity is a frequency-dependent point property in the volume that depends on depth. The group velocity is obtained as a result of the time-of-flight measurement between the (virtual) source and the receiver.

[0054] The phase velocity is the velocity at which a particular frequency component of a wave travels. Therefore, the phase velocity is expressed as a function of frequency. To measure the phase velocity, at least two measurement nodes are selected to measure the wave propagating through the volume to determine the relative time-of-flight between geophones for different frequencies. The result is the phase velocity, which is a function of the frequency averaged over the volume between the two measurement nodes. The phase velocity is obtained as a point in a two-dimensional phase space (dispersion spectrum), which is obtained by performing a two-dimensional transformation (such as slant stack, Radon, FK, etc.) on an array of recorded waveforms (time-distance space).

[0055] Referring again to Figure 1 , each geophone 104 provides a measurement node for measuring the vertical component of the passing surface wave. The geophone 104 can be configured to measure the vibrations caused by ambient noise. That is, the background wave field caused by natural or man-made noise (rather than impulse points such as explosions or hammer blows used in active methods).

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

[0057] Figure 2 Multiple virtual source-receiver pairs on the surface above the subsurface target volume are shown. The ray paths 202 between the source-receiver pairs are marked, in particular the ray paths starting from a single central geophone near the center of the geophone array and each of the other geophones in the array. The background shading and contour rings represent the travel-time field from the central geophone to the other geophones. There are also corresponding ray paths between each individual geophone and all the other geophones, i.e., between each pair of geophones, which are not depicted for simplicity. Figure 2 in the figure.

[0058] Each pair of geophones can provide a signal at a first location for cross-correlation with a signal at a second location, thereby reproducing a virtual source-receiver pair using the principle of interferometry. Specifically, Figure 1 the cross-correlation of the passive noise measured at the corresponding geophone pairs shown at the surface can be used to reproduce the response from the subsurface target volume as if it were caused by a pulsed point source, which is equivalent to the Green's function.

[0059] In other words, cross-correlating the responses received by two receivers can be interpreted as the response that would be measured at one of the receiver locations as if there were a source at the other receiver location. There are various known methods for determining the Green's function of a virtual source-receiver pair, and an overview of these methods is described in the following reference, "Tutorial on Seismic Interferometry: Part 1 – Basic Principles and Applications"; Geophysics, Vol. 75, No. 5 (September - October 2010; pp. 75A195 - 75A209; Wapenaar et al.).

[0060] Reference Figure 3 , the first shear-wave velocity model 300 has cells m on the surface above the subsurface target volume xy grid, including columns m extending in the x-direction x1 , m x2 , m x3 etc., and rows m extending in the y-direction 1y , m 2y , m 3yetc. Each cell defines a region of the earth's surface. For example, each cell can define a square of 5 meters by 5 meters; other example options include a square of 1 meter by 1 meter or a square of 10 meters by 10 meters. In other words, the shear wave velocity model includes a plurality of cells arranged in a two-dimensional grid 300. The two-dimensional grid at least spans the area of the earth's surface above the underground target volume. Selecting a smaller cell area increases the resolution of the model. Each cell also includes a volume that extends vertically below the area of the earth's surface. The model 300 can extend infinitely below the area of the earth's surface or to a predetermined depth below the earth's surface, for which the shear wave velocity value has a significant effect on the propagation of surface waves propagating at the earth's surface. For example, the model 300 can be defined to a depth of up to 50 meters, 100 meters, 200 meters, or 300 meters.

[0061] In some examples, the first shear wave velocity model 300 is not a grid of square cells, but includes cells having other tessellation shapes such as rectangles, rhombuses, or combinations including non-uniform or different shapes.

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

[0063] In other examples, the model can 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 used, optionally also the viscosity quality factors Qs and Qp. Generally, the model can define multiple physical property values.

[0064] Reference Figures 4A to 4D , depicts various representations of the wave path between a (virtual) source and a receiver. The wave path is the path taken by a surface wave from the source to the receiver. As described below with reference to Figure 4A and 4B The wave path can be represented by a straight or curved ray connecting the source and the receiver. Alternatively, as described below with reference to Figure 4CAs described, the wave path can be represented by an elliptical Fresnel zone (also known as the Fresnel kernel) located between a source-receiver pair. As an alternative, as described below with reference to Figure 4D the wave path can be represented by a curved ("banana-shaped") Fresnel zone.

[0065] With reference to Figure 4A , the ray path is defined as the straight line between two geophones at surface locations A and B on model 300. The direct ray path model assumes a laterally constant velocity (in other words, the velocity only varies with depth). In many cases, the direct ray path model provides a sufficiently accurate approximation of the actual path traveled by the wave between the source and the receiver. Compared to other more complex techniques, the direct ray path approximation is also more efficient mathematically and computationally. According to ray theory, the ray path characterizes the motion of the surface wave as it travels from A to B (or from B to A). Specifically, the ray path is defined as the propagation direction of the surface wave, i.e., the direction perpendicular to the wavefront in wave theory or perpendicular to the travel time contour lines. Figure 4A The ray path shown passes through seven cells of model 300, which are numbered 1 to 7. This ray path consists of seven ray segments, and each ray segment has a length L depending on the length of the ray path traveling through each cell i , i.e., Figure 4A L1 to L7 in

[0066] With reference to Figure 4B , the ray path is defined as a curve between two geophones at surface locations C and D. Figure 4B The ray path shown passes through six cells of model 300, numbered 1 to 6. This ray path consists of six ray segments, and each ray segment has a length L depending on the length of the ray path traveling through each cell i , i.e., Figure 4B L1 to L6 in

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

[0068] More specifically, and as will be understood by those skilled in the art, any wave propagating along the path between the (virtual) source at point E and the receiver at point F will have some off-axis propagating wave components (not along the straight line connecting E and F). The off-axis propagating wave components may be deflected by the wave propagation medium, and some of the deflected waves will be directed to the receiver at F. Therefore, the travel time of the wave traveling along the direct (straight) path is different from that of the wave traveling along the deflected path, which means that the direct path wave and the deflected path wave will arrive at the receiver at different phases. When the phase difference is half a wave period (or one and a half wave periods, or two and a half wave periods, or any odd integer number of half periods), the phase difference may cause destructive interference. On the contrary, when the phase difference is between zero and half a wavelength period (or between one and one and a half wavelength periods, or between any integer n - 1 and n - 1 / 2 wavelength periods), the phase difference will cause constructive interference. In other words, the nth Fresnel zone is defined as the region in which the waves deflected at a point will have a phase difference between n - 1 and n - 1 / 2 wavelengths with respect to the wave traveling along the straight line between the source and the receiver pair. In the example of the present disclosure, for the wave path between the source and the receiver, only the first-order Fresnel zone may be considered because the wave components in the higher-order Fresnel zones become very small. B.D. Guenther (2005) described the general physical principle of the Fresnel zone of any wave propagating through any medium in the "Encyclopaedia of Modern Optics", and those skilled in the art will understand that the description of the Fresnel zone therein applies to any propagation medium. Shibo Xu and Alexey Stovas (2018) also described the Fresnel zone of the wave traveling through the underground target volume in "Fresnel zone in VTI and orthorhombic media".

[0069] As Figure 4C shown therein, which represents the first-order Fresnel zone for the direct ray approximation, the Fresnel zone is elliptical, where the (virtual) source position E and the receiver F are the foci of the ellipse. In other examples, such as Figure 4D shown therein, the Fresnel zone may be curved, i.e., a so-called banana-shaped region, along the curved ray path between the virtual source G and the receiver H. Determining such a curved Fresnel zone follows the same considerations as described above for the Figure 4A elliptical Fresnel zone, except that the differential travel time and the phase difference are determined for the curved ray (see Figure 4B ).

[0070] Figure 4CThe Fresnel zone shown in [reference] covers eight cells of the model 300, numbered 1 through 8. The Fresnel zone contains eight sub-regions labeled R1 through R8, each of which is the region of the Fresnel zone within the corresponding cell 1 through 8. Each sub-region R1 through R8 of the Fresnel has a corresponding sensitivity value, which is uniquely calculated for the Fresnel zone. Those skilled in the art will understand that the sensitivity value is the sensitivity to local shear wave velocity changes (or any other physical property being modeled). This technique for calculating the sensitivity values of the individual sub-regions of the Fresnel zone (i.e., the Fresnel zone sensitivity values for the individual cells covered by the Fresnel zone) is well known to those skilled in the art.

[0071] For Figure 4D the curved Fresnel zone depicted in [reference] also shows the same sub-region illustration, which passes through cells 1 through 8 and consists of sub-regions R1 through R8 (different from the sub-regions depicted in Figure 4C [reference]). Figure 4D Each sub-region shown in [reference] has a corresponding sensitivity value, and those skilled in the art will understand that these sensitivity values may be different from the sensitivity values of each sub-region depicted in Figure 4C [reference].

[0072] For each of a plurality of virtual source-receiver pairs such as those shown in Figure 2 [reference], a corresponding Fresnel zone can be determined. Each Fresnel zone consists of a corresponding plurality of sub-regions for each cell covered by that Fresnel zone, and each of those sub-regions in each region has a corresponding sensitivity value.

[0073] Referring to Figure 5 [reference], a method 500 for determining one or more physical properties of a subsurface target volume includes receiving 502 a plurality of signals detected by geophones arranged on the surface of the earth, e.g., as described above with reference to Figure 1 [reference]. Typically, one signal is received from each geophone, but in some cases, signals may only be received from a subset of the total number of geophones. These signals represent the vertical vibrations of the ground measured at the corresponding geophones, as caused by surface waves from ambient noise or an active source, as described above with reference to Figure 1 [reference]. These signals may contain metadata regarding the recording location and time.

[0074] The signals can be received directly from the geophones or through one or more intermediate devices. For example, a computing device (as described below with reference to Figure 10The method 500 may be performed remotely from the geophone array, such as by Internet communication or by transmitting a physical computer-readable medium on which a record of the signal is stored. In some examples, additional processing steps may be performed on the signal before or after receiving 502 the signal to optimize the signal for further processing.

[0075] Example geophones may include velocimeters or accelerometers. Examples of specific mechanisms of such geophones include: a ferromagnetic mass on a spring moves within an electrical coil in response to ground surface motion, thereby inducing a measurable current proportional to the ground surface velocity.

[0076] The method 500 also includes step 504 of cross-correlating the signals between the plurality of geophones to obtain empirical travel time data for surface waves between a plurality of virtual source-receiver pairs. As described above, this step involves cross-correlating the passive noise signals measured at a pair of receivers to obtain a correlated Green's function that represents the wavefield as if one of the pair of receivers were a virtual source and the other of the pair were a receiver. As will be appreciated by those skilled in the art, information indicating phase velocity as a function of frequency and information indicating group velocity as a function of frequency may be obtained by the cross-correlation process. In step 504, the signals between the plurality of geophones are cross-correlated to obtain group and / or travel time data for a plurality of frequencies. These frequencies may be continuous in order to obtain group velocity and / or phase velocity dispersion information, or 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 a virtual source to a receiver. Phase travel time refers to the time taken for a particular phase component of a wave to travel from a virtual source to a receiver. Since the distance between the virtual source and the receiver is known, the group velocity and / or phase velocity or slowness (which is the inverse of the velocity) at multiple frequencies can be derived from the relevant Green's function for the source-receiver pair. Slowness and travel time are actually equivalent, with the scaling factor provided by the distance between the source and the receiver.

[0077] Based at least on the group and / or phase travel time data for the plurality of frequencies, a model of the physical properties of the subsurface target volume may be determined. Figure 3 As described, the model comprises a plurality of cells arranged in a two-dimensional grid. Each cell of the model has a corresponding physical property value (e.g., shear wave velocity), which may include a depth profile. In other words, each cell of the model may comprise a one-dimensional model of the shear wave velocity as a function of depth or frequency (which indicates depth). Since each cell provides a one-dimensional model as a function of depth, collectively, the cells of the model provide a three-dimensional model of the physical property.

[0078] A model for determining the physical properties of an underground target volume involves selecting a subset of more than 506 source-receivers for which empirical travel-time information has been obtained. As referred to above Figures 4A - 4D and as described, the wave path between each source-receiver pair passes through two or more cells of the model. In principle, any and all source-receiver pairs may be applicable to determining the model. However, increasing the number of wave paths used increases the computational power required or alternatively results in a very long computation time. Thus, in practice, it is often beneficial to select a subset of source-receiver pairs, the number to be selected depending on the available computational power or other practical considerations. A variety of principles may be used to limit the number of pairs selected. First, due to the reciprocity of cross-correlation, the wave path from position A to position B is the same as the wave path from position B to position A, so only pairs of geophones need to be selected without regard to which is considered the virtual source. The first method of deselected pairs is to ignore those for which there is no suitable travel-time data, i.e., the cross-correlation of the signals from that geophone pair does not show an identifiable wave passing through the two geophones that can be used to calculate the travel time between them. Another method of selecting wave paths is to ignore pairs with low-quality picks, such as those with a low signal-to-noise ratio, with a large uncertainty, being outliers compared to adjacent picks, etc. Any pick that gives non-physical or non-geological results, as well as any pick that would not be credible for the particular underground target volume being investigated, will also be excluded. Finally, if the remaining number of (source-receiver) pairs still exceeds what is actually available, a subsample of source-receiver pairs may be selected, such as every n pairs, or randomly, etc. In this case, other suitable subsamples of source-receiver pairs may be used in succession to increase the total number of pairs used, but staying within the memory limit or available computational power. Regardless of the method used to limit the number of source-receiver pairs (if necessary), the end result is the selected source-receiver pairs and the corresponding wave paths for which empirical group and / or phase travel-time information has been obtained, such as the velocities for multiple different frequencies or dispersion functions (e.g., group velocity dispersion function and / or phase velocity dispersion function).

[0079] As referred to above Figure 4A and 4B and as described, the wave path between the (virtual) source and the receiver pair may be a ray path that extends in a straight line between the source and the receiver. In some examples, at least some of the ray paths extend in a curve between the respective source and receiver, as referred to above Figure 4BAt least some of the selected wave paths will pass through at least two cells of the first model, having ray segment lengths in each of the cells through which they pass (including the cells in which the ray path starts and ends). The curved ray path can be determined based on the initial physical property values of the cells of the model. By using the curved ray path determined based on the physical property values of the cells of the model, the ray path more closely follows the actual ray path along which the wave will travel between the source and the receiver. Thus, this method improves the accuracy of the results of subsequent method processes because it more closely corresponds to the physical reality of the surface wave passing through the subsurface target volume.

[0080] As described above with reference to Figures 4A - 4B In other examples, the wave path is represented by a Fresnel zone, which provides a region of multiple feasible paths that the wave may take as it travels from the source to the receiver. The Fresnel zone provides a region covering two or more cells, through which possible waves will pass, each cell having a region of the corresponding Fresnel zone with a related sensitivity value defined by the Fresnel zone. The Fresnel zone can be determined based on the initial physical property values of the cells of the model. Since the Fresnel zone defines the region of possible paths taken by the surface wave as it travels from the source to the receiver, this method improves the accuracy of the resulting model for each cell because it more closely corresponds to the physical reality of the surface wave passing through and scattering within the subsurface target volume.

[0081] The model for determining the physical properties of the subsurface target volume also includes performing tomographic inversion on the travel time data for each selected source-receiver pair. Generally, tomographic inversion is used to derive the physical property values of the individual cells in the model from the group or phase travel time data obtained for each of a plurality of source-receiver pairs. As described above, the physical property values of the individual cells can be a single value or provided as a physical property as a function of depth. Details of various processes for performing tomographic inversion will be described below with reference to Figures 6 - 8 The details of performing tomographic inversion. As described throughout this document, tomographic inversion can be performed based on empirical group travel time information or empirical phase dispersion information. In the case of phase dispersion, phase velocity or slowness information can be converted to equivalent travel time information so that tomographic imaging can be performed based on such phase travel time information.

[0082] In some examples, the tomographic inversion step 508 involves Figure 6 the process 600 depicted in Figure 6A method of traveltime tomography of the traveltime information obtained in step 504 for each source-receiver pair selected in step 506 is shown. As described in more detail below, this process maps the traveltime information from the individual source-receiver pairs tomographically into the cells of the physical property model. Thus, the result of this process is an empirical model of the group velocity and / or phase velocity (depending on the traveltime data used), where the individual cells of the model have corresponding group velocity values or phase velocity values. This process 600 is performed for each of a plurality of frequencies to obtain the group velocity values or phase velocity values for each cell for each frequency. Process 600 is the first stage (i.e., the tomography stage) of a two-stage tomographic inversion.

[0083] Process 600 includes step 600 of obtaining an initial velocity model. The initial velocity model includes initial values of the group velocity and / or phase velocity for each of the individual cells for each of a plurality of frequencies. The initial velocity model sets initial (group or phase) velocity values for the individual cells (for each given frequency), which process 600 will refine iteratively using the empirical traveltime data. Thus, the initial velocity model and the corresponding initial velocity values do not have to be of high accuracy or high resolution, although a more accurate initial model can make the tomography process 600 faster or more accurate in mapping the traveltimes to the individual cells. A more accurate initial model can also reduce the risk of finding a local minimum rather than a global minimum in the iterative gradient descent method, although such problems can be addressed using Monte Carlo methods. In some examples, the initial velocity model is determined based on user input, such as according to historical data or map information indicating possible physical property values for the entire subsurface target volume. Alternatively, typical values of the group / phase velocity arbitrarily selected for the individual cells can also be used as a starting point. In these examples, any model can be selected based on an estimate of the physical properties of the subsurface target volume.

[0084] The tomography process 600 also includes step 604: using the initial group or phase velocity model to determine the modeled group or phase traveltime for each selected source-receiver pair. This step is performed by identifying the wave path from the source to the receiver (e.g., a straight ray, a curved ray, or an elliptical / curved Fresnel zone) and identifying the cells through which the wave path passes. Then, the modeled traveltime based on the initial velocity model is determined based on the known distance between the source and the receiver and the velocity values of the individual cells through which the wave path passes.

[0085] The tomography process 600 also includes step 606: determining an error value that indicates the difference between the modeled travel times (determined in step 604) and the empirical travel times (obtained in step 504) for each selected source-receiver pair. For example, the error value can be the simple difference between the modeled travel time and the empirical travel time at each frequency for each source-receiver pair, also known as the residual. In some examples, the error value can be a combination of all differences between the modeled travel times and the empirical travel times for each source-receiver pair. Typically, determining the error value is part of an iterative process, such as a least squares inversion method, where the square of the residual is calculated in order to minimize it through iteration. Other forms of inversion processes use different error values in order to provide feedback to the initial group velocity model.

[0086] Process 600 also includes step 608: using the error value to determine an updated velocity model. The updated model is generally the same in all respects as the initial velocity model, except for the new velocity values associated with at least some cells. In other words, the updated model is an updated version of the initial velocity model that takes into account the determined error value between the empirical and modeled group or phase travel times for each source-receiver pair. This feedback process may involve least squares methods, Markov chain Monte Carlo methods, or other inversion techniques to iteratively update the initial model based on the updated model. Process 600 can be repeated for each of a plurality of finite frequencies.

[0087] Steps 604 through 608 are generally all part of a subroutine of the first stage tomography process 600, which is then iterated according to an inversion method, such as least squares inversion. Figure 6 The dashed arrows in show the iterative nature of process 600, which indicates that the updated model determined at step 608 is used to determine new modeled travel times for each source-receiver pair at step 604. In other words, each time the error value obtained from the initial model and the resulting modeled travel times for each source-receiver pair is used to determine an updated model, the resulting updated model is 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 drops below an absolute value threshold of the difference between the modeled travel time and the empirical travel time, or drops below a threshold of the proportional difference between the modeled travel time and the empirical travel time. Another end condition that can be used alone or in combination with the threshold error value end condition is that the modification made to the initial model in order to generate an updated model for the iteration drops below a threshold amount or proportion. In this way, if the iteration reaches a set minimum error value, the iteration can end because subsequent iterations will not significantly improve the accuracy.

[0088] Once tomography of the first stage process 600 has been performed to obtain a velocity model containing group velocity values or phase velocity values for each cell (for each frequency), the tomographic inversion of step 508 can enter the second stage (inversion) of the two-stage tomographic inversion to obtain a final model of the physical properties of the target subsurface volume. The second stage inversion process 700 is depicted as Figure 7 shown. Similar to the first stage tomography process 600, the second stage inversion process 700 is performed for each of a plurality of frequencies to obtain group velocity values or phase velocity values for each cell at each frequency.

[0089] The inversion process 700 includes a first step 702 of obtaining an initial physical model. The initial physical model is an initial model of the physical properties of the subsurface volume, where the initial model contains initial physical property values for each cell in a first plurality of cells. The initial physical model can be the shear wave velocity model 300 as described above with reference to Figure 3 and / or an initial physical model having any of the features or variations as described above with reference to Figure 3 The initial model sets initial physical property values for each cell in the model, and method 700 will use the group / phase velocity model obtained from the iterative tomography process 600 to refine these initial physical property values through further iterative processes. Thus, the initial physical model and the initial physical property values do not necessarily have to be a highly accurate or high-resolution model of the target subsurface volume, although a more precise first model can improve the expected accuracy of the final result of the method or reduce the computational time to reach the final result.

[0090] In some examples, the initial physical model is determined based on user input, such as historical data or map information indicating possible physical property values of the subsurface target volume. Alternatively, the initial physical model can be determined using received signals, such as by inverting the group velocity or phase velocity dispersion curves between geophones, which can be calculated from the signal cross-correlation as described above, in order to find using a coarser grid or a faster method. As a last resort, for each cell, arbitrarily selected typical values of physical properties can be used as a starting point, which can be based on the estimated physical properties of the subsurface target volume. In some examples, the initial physical model can be determined based on empirical phase dispersion data between source-receiver pairs. The phase velocity information is obtained as a point in a two-dimensional phase space (dispersion spectrum), which can be obtained by performing a two-dimensional transform on a waveform array (time-distance space), which can be a Radon transform, a slope-relaxation transform, FK, or any other suitable method obvious to those skilled in the art. A dispersion curve is generated that shows the phase velocity as a function of frequency, which can provide a preliminary three-dimensional model of the shear velocity in the three-dimensional model. However, since the acquisition of the phase velocity is essentially a measurement of the frequency-dependent average velocity between nodes, the resolution of this model is relatively low. Therefore, a physical model derived only from empirical phase dispersion data can be used as an initial model for inversion, which also takes into account empirical group dispersion data to obtain a more accurate final three-dimensional model of the subsurface target volume.

[0091] Process 700 also includes step 704 of determining the modeled surface wave velocity for each cell. More specifically, a forward modeling approach is employed to derive the phase velocity and / or group velocity for each cell in the initial physical model using the corresponding physical property values for a given cell. For example, the physical property value for each cell can be the shear wave velocity. The value of the shear wave velocity for each cell can be used to calculate the corresponding phase velocity dispersion function (and thus the phase velocity at a given frequency) using the propagation matrix method, which was proposed 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). Those skilled in the art will be well aware of this method and other methods that can be used to forward model the group or phase velocity (or dispersion function) from shear wave velocity values (or functions of depth).

[0092] The inversion process also includes step 706 of determining an error value that indicates the difference between the modeled velocity for each cell of the model (determined in step 704 based on the initial physical model) and the empirical velocity (obtained from process 600). The error value for each cell can be determined in a similar manner as described above for step 606, where the error value between the modeled travel time and the empirical travel time for each source-receiver pair is determined. As with tomographic process 600, generally, determining the error value in step 706 is part of an iterative process, such as in a least squares inversion method, where the square of the residual is calculated and minimized through iteration. 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.

[0093] Process 700 also includes step 708 of determining an updated physical model based on the error value determined at step 706. The updated physical model is generally the same in all respects as the initial physical model, except for the 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, where the determined error value between the empirical surface wave velocity and the modeled surface wave velocity of each cell of the model is taken into account. This feedback process may involve least squares methods, Markov chain Monte Carlo methods, or other inversion techniques to iteratively update the initial physical model based on the updated physical model.

[0094] Steps 704 through 708 are generally all part of a subroutine of the second stage inversion process 700, which is then iterated according to an inversion method, such as least squares inversion. Figure 7 The dashed arrow in indicates the iterative nature of process 700, which indicates that the updated model determined at step 708 is used to determine the new surface wave velocity of each cell at 704. In other words, each time an updated model is determined using the error value obtained from the initial model and the modeled velocities of the individual cells, 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 the end criteria described above with respect to Figure 6 the end criteria described above.

[0095] Now refer to Figure 8 the alternative tomographic inversion process 800 is described. This tomographic inversion process is a single stage process, where tomography and inversion are combined. This means that the iterative inversion stage for finding the physical property values of the model is performed for each individual wave path between the corresponding source-receiver pairs, rather than for each individual cell of the model. This means that in the single stage process 800, there is no need for an initial tomography stage (process 600) to map the surface wave velocity between source-receiver pairs to the cells, because the physical property model inversion process performs the inversion on each individual wave path rather than on each individual cell. In other words, the two-stage tomographic inversion first uses travel time tomography to map the travel time information along the wave path to surface wave velocity information for the cell grid. This two-stage tomographic inversion then performs a cell-by-cell inversion of the surface wave velocity information (surface wave velocity as a function of frequency) to obtain the corresponding shear wave velocity function (as a function of depth) for each cell. In contrast, as described below, the single stage process performs a direct inversion of the surface wave velocity information for each individual wave path (for each source-receiver pair) to obtain the shear wave velocity function for that wave path, which is mapped in a tomographic manner to each individual cell of the model.

[0096] Process 800 first obtains an initial physical model, in a manner similar to Figure 7The same as that described in step 702. The initial model sets initial physical property values for each unit of the model, and process 800 will refine it using the empirical travel time data of each source-receiver pair. In some examples, the initial physical model can be determined based on the empirical phase dispersion data between the source-receiver pairs.

[0097] Process 800 also includes determining 804 the empirical dispersion function for each wave path (i.e., each source-receiver pair), specifically, for each wave path selected in the selection 506 part of the method. The dispersion function can be a group velocity dispersion function or a phase velocity dispersion function (or both can be used). The dispersion function can be determined according to any common method in the art as described above with reference to Figure 1 that.

[0098] Process 800 also includes determining 806 the modeled dispersion function for each wave path using the first model. For example, the value of the shear wave velocity of the first model can be used, and a model of a stack of homogeneous layers of finite thickness (covering a homogeneous half-space) can be used to calculate the phase velocity dispersion function using the propagation matrix method proposed 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).

[0099] Specifically, this can be carried out using any available solver that employs the modal approximation method, such as that written by Herrmann in ((2013) "Computer programs in seismology: An evolving tool for instruction and research", Seismological Research Letters, 84(6), 1081–1088). This method calculates the phase velocity Vphase at a specific frequency ω from a stack of homogeneous layers of finite thickness (i.e., a one-dimensional velocity profile). The first step involves determining the boundaries between which the phase velocity will lie, for example, between the minimum and maximum Vs present in the model. For any value of ω, trial values of the wavenumber k are then tested and iteratively varied in a set increment until an eigenvalue solution is found. Subsequently, this eigenvalue can be used to calculate the solution of the eigenfunction, which will be tested for the iteratively varying values of Vphase to find the root of the solution and thus the phase velocity dispersion function.

[0100] Herrmann's modal approximation code takes as input a list of layered elastic parameters and a frequency value of interest and outputs a dispersion function that is the solution to an eigenvalue problem defined by the matrix propagation method proposed by Haskell and Thomson. The Herrmann solver is widely used and can be obtained from the Earthquake Center at Saint Louis University as a software package - the "Computer Programs in Seismology" page at https: / / www.eas.slu.edu / eqc / eqccps.html, along with a large number of manuals and instructions for using the software (in addition to the 2013 paper "Computer programs in seismology: An evolving tool for instruction and research" mentioned above). Alternative solvers can also be obtained via GitHub, such as https: / / github.com / xin2zhang / MCTomo based on 3D Monte Carlo tomography (Zhang, X., Curtis, A., Galetti, E., and de Ridder, S., 2018. "3-D Monte Carlo surface wave tomography", Geophysical Journal International, 215(3), 1644–1658), or Keurfon Luu's https: / / github.com / keurfonluu / disba (also available at zenodo.org under the title keurfonluu / disba: disba v0.5.1, which, according to its documentation page, implements "a subset of codes from Computer Programs in Seismology (CPS) in Python compiled just-in-time with numba").

[0101] In addition to, or as an alternative to, shear wave velocity, other physical properties related to shear wave velocity can be used. For example, the compressional wave (P-wave) velocity Vp and the shear wave (or transverse wave / 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 elastic theory:

[0102] Equation 1 (compressional wave velocity):

[0103]

[0104] Equation 2 (Shear Wave Velocity):

[0105]

[0106] In summary, there are established methods for determining the phase velocity dispersion function from a stack of homogeneous layers of finite thickness (constant in the x and y directions) and the associated shear wave velocity values (or associated physical properties) for each layer. Thus, one way to use the methods developed by Thomson, Haskell, and Herrmann is to use multiple wave paths between geophones that intersect within a particular cell of the model and combine the dispersion functions of the wave paths to determine the dispersion relation for that particular cell. The dispersion relation for that cell can then be used to "invert" using the homogeneous layer method with the shear wave velocity value for that cell, i.e., to improve the model values in that cell to better generate the dispersion relation present for that cell (and perform a similar operation for all cells individually). However, the present inventors have developed a different method in which, instead of inverting each cell individually, the entire wave path is inverted using multiple shear wave velocity values in the cells through which the wave path passes. This method has many advantages, including increased computational efficiency, allowing for more geophones in the geophone array and / or a higher resolution model, and increased accuracy in model determination. For example, since the shear wave velocity model is improved in a single-step manner from the dispersion function present between the geophones, this eliminates the additional step of performing travel-time tomography on the wave paths for each cell, which is an additional source of approximation and thus inaccuracy since the inversion is independent for each cell (and thus the interdependence of the cross wave paths is lost in the inversion step). In addition to eliminating the additional step in the process and thus reducing the amount of computational power required, performing the inversion in a single step also means that the least squares method can be used instead of the more computationally intensive Monte Carlo inversion method. This in turn means a reduction in the inversion calculation time, or for the same amount of calculation time, more geophone signals can be processed, or a higher resolution model can be generated.

[0107] In some examples, determining the modeled dispersion function 806 for each wave path is done by averaging the shear wave velocity values of the cells through which the wave path passes. Specifically, for the shear wave velocity, the values of 1 / Vs are averaged because the travel time through each cell is inversely proportional to the velocity. For other physical properties of the cells, such as density, it may be appropriate to simply arithmetically average the values in each cell. In either case, the values are still a function of depth. The average shear wave velocity (or other average physical property) is then used in the Hermann mode approximation to determine the modeled dispersion function for the entire wave path. This averaging process can be implemented as an additional function in the solver software before running the inversion program.

[0108] In some examples, the averaging of the shear wave velocity values is a weighted average, which is weighted according to the ray segment lengths in each cell, as shown in Equation 3, which applies to both straight and curved wave paths. For wave paths represented by the Fresnel zone, the weights will also include an additional term that corresponds to the sensitivity of the wave path to local velocity changes at different parts of the Fresnel zone (i.e., there is higher sensitivity at the center of the zone than at the periphery). Using a weighted average improves the accuracy of the calculation because it more closely corresponds to the actual values that the wave path will experience between geophones.

[0109] Equation 3:

[0110]

[0111] Where:

[0112] Vs A->B is the average shear wave velocity (as a function of depth) of the wave path from A to B;

[0113] N is the number of cells of the first model through which the wave path from A to B passes;

[0114] L i is the ray segment length of the i-th cell along the wave path;

[0115] L A->B is the total path length from A to B, i.e., the sum of all L i ;

[0116] Vs i is the shear wave velocity value (as a function of depth) of the i-th cell along the wave path.

[0117] In some examples, instead of first averaging the shear wave velocity values along the wave path length, a modal approximation method is performed to generate a dispersion function for each element, and then the element dispersion functions of the individual elements through which the wave path passes are averaged to generate a modeled dispersion function for the wave path. This can be done by averaging the phase velocity dispersion function of each element along the wave path with respect to 1 / Vphase (or conversely, averaging the group velocity dispersion function 1 / Vgroup), optionally including weights similar to the above averaging. This averaging process can be implemented in the solver software as an additional function after running the inversion program.

[0118] The method further includes determining an error value at 808, which indicates the difference between the modeled dispersion function (determined from the model at 806) and the empirical dispersion function (determined from the detected signal at 804) for each wave path (i.e., each source-receiver pair). The error value can be determined for each wave path in a manner similar to that described above with respect to steps 606 and 706. Similar to the two-stage tomography processes 600 and 700, typically, determining the error value at 808 is part of an iterative process, such as calculating the square of the residuals in a least squares inversion method in order to minimize it through iteration. 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.

[0119] Process 800 further includes step 810 of determining an updated physical model based on the error value determined at 808. The updated physical model is generally the same in all respects as the initial physical model, except for the new physical property values associated with at least some of the elements. In other words, the updated model is an updated version of the initial physical model, where the determined error value between the empirical surface wave velocity and the modeled surface wave velocity of each element in the model is taken into account. This feedback process may involve least squares methods, Markov chain Monte Carlo methods, or other inversion techniques to iteratively update the initial physical model based on the updated physical model.

[0120] In some examples, the determined error value indicates that the empirical dispersion function along the wave path shows phase velocity (or group velocity) values that are greater than the modeled dispersion function over the entire function or in one or more frequency bands. Thus, this indicates that the value of the shear wave velocity along this wave path is less than the true physical properties of the subsurface target volume. Therefore, using the error value to determine the second model includes recording a new shear wave velocity value for the elements through which the wave path passes, which is greater than the corresponding value in the first model. Conversely, if the empirical dispersion function has a smaller value than the modeled dispersion function, the value of the shear wave velocity along this wave path is greater than the true value in the subsurface volume, and the model should be updated accordingly.

[0121] One way to update the first model to determine the second model is to find the empirical average shear wave velocity of the wave path through the inversion of the dispersion function, and scale each shear wave velocity value from the first model according to the difference between the empirical average shear wave velocity and the average shear wave velocity of the first model for that wave path. More specifically, the update factor for the first model for each wave path can be defined as the following ratio: the difference between the empirical average shear wave velocity and the average shear wave velocity of the first model for that wave path; and the average shear wave velocity of the first model for that wave path. In practice, this requires an arithmetic average of 1 / Vs (also known as "slowness") because cells with smaller shear wave velocities will have a greater impact on the average shear wave velocity over the entire wave path. Optionally, updating the first model to produce the second model includes weighting the variation of the shear wave velocity values of each cell according to the ray segment length of the wave path in each cell, as explained in the method for modeling the dispersion function part of the determination 510 method 500 above.

[0122] In some examples, the shear wave values of the first model are updated using the sensitivity matrix method to generate the second model. In this method, the partial derivatives of the error value (or residual) with respect to the shear wave velocity and / or other parameters of the first model are determined, and the other parameters are, for example, the boundary depth between shear wave velocity values in a given cell, the compressional wave velocity, the density, etc. The partial derivatives are combined into a sensitivity matrix, and a system of linear equations is solved to determine how to change the parameters of the first model so as to reduce the overall error value of the second model with the updated values.

[0123] Once the updates for the values of the first cell along the first wave path (e.g., ray or Fresnel zone) are calculated, these updates can be directly applied to generate the second model. Alternatively, the updates for the values of all cells affected by each wave path calculation and inversion can be combined and applied to the values of the first model all at once to generate the second model.

[0124] Steps 806 to 810 are generally part of the subroutine of a single-stage tomographic inversion process 800, which then iterates according to the inversion method (e.g., least squares inversion). Figure 8 The dashed arrow in indicates the iterative nature of process 800, and this arrow indicates that the updated model determined at step 810 is used to determine the new surface wave velocity of each cell at 806. In other words, each time the error value obtained from the initial model and the modeled velocities of each cell is used to determine the updated model, the resulting updated model will then be used as the initial model for the next iteration. The iteration continues until the error value reaches an end condition, such as the end criteria described above regarding Figure 6 and Figure 7 described.

[0125] Refer to Figure 8, although parts of process 800 have been sequenced and numbered, the method is not limited to the specific sequence of parts presented herein. For example, there is no requirement that either the empirical dispersion function or the first modeled dispersion function be performed before the other.

[0126] In some examples, once the final second model is determined, e.g., once the iteration satisfies an end condition, the final second model can be sent to an output device. For example, the final second model can be displayed on a display or other user interface, or sent via wired or wireless communication to another computing device for further processing or display.

[0127] Reference Figure 9 , an example output of the method described herein is the final shear wave velocity model, which shows the subsurface target volume. The subsurface target volume extends to a depth of 100 meters (z-direction) in the x and y directions. The values of the shear wave velocity are shown by shading in the figure, and the transitions between regions of different shear wave velocities can be seen, which indicate different compositions or structures of parts of the subsurface target area. Using the method described herein, the final shear wave velocity model can be determined with higher resolution and accuracy without resulting in an infeasible computation time. Such subsurface models can be used to better understand whether the subsurface target volume is suitable to support an artificial structure on top of or within the subsurface target volume.

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

[0129] Figure 10A block diagram of one embodiment of a computing device 1000 is shown, in which a set of instructions can be executed to cause the computing device to perform any one or more of the methods discussed herein. In alternative embodiments, the computing device may be connected (e.g., networked) to other machines on a local area network (LAN), intranet, extranet, or the Internet. The computing device may operate in a client-server network environment with the capabilities of a server or a client machine, or as a peer machine in a peer-to-peer (or distributed) network environment. The computing device may be a personal computer (PC), tablet computer, set-top box (STB), personal digital assistant (PDA), cellular phone, network device, server, network router, switch, or bridge, or any machine capable of executing a set of instructions (sequentially or otherwise) that specify actions to be taken by that machine. Further, although only a single computing device is shown, the term "computing device" should also be understood 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 methods discussed herein.

[0130] 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 an auxiliary memory (e.g., data storage device 1018), which communicate with each other via a bus 1030.

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

[0132] Computing device 1000 may also include a network interface device 1008. Computing device 1000 may also 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 a touch screen), a cursor control device 1014 (e.g., a mouse or a touch screen), and an audio device 1016 (e.g., a speaker).

[0133] It will be apparent that Figure 10 some features of the computer device 1000 shown in may be missing. For example, one or more computing devices 1000 may not require a display device 1010 (or any associated adapter). For example, this may be the case for certain server-side computer devices 1000 that are only used for their processing functions and do not need to display information to a user. Similarly, a user input device 1012 may not be required. In its simplest form, the computer device 1000 includes a processor 1002 and a memory 1004.

[0134] The data storage device 1018 may include one or more machine-readable storage media (or more specifically, one or more non-transitory computer-readable storage media) 1028, on which one or more instruction sets 1022 are stored, and these instruction sets 1022 embody any one or more of the methods or functions described herein. The instructions 1022 may also reside, in whole or at least in part, in the main memory 1004 and / or the processor 1002 during their execution by the computer system 1000, and the main memory 1004 and the processor 1002 also constitute computer-readable storage media.

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

[0136] In one embodiment, the modules, components, and other features described herein may be implemented as discrete components or integrated into the functions of hardware components (such as ASICS, FPGAs, DSPs, or similar devices).

[0137] "Hardware component" refers to a tangible (e.g., non-transitory) physical component (e.g., a set of one or more processors) that can perform a specific operation and can be configured or arranged in a specific physical manner. A hardware component can include dedicated circuitry or logic that is permanently configured to perform a specific operation. A hardware component can be or include a dedicated processor, such as a field programmable gate array (FPGA) or an ASIC. A hardware component can also include programmable logic or circuitry that is temporarily configured by software to perform a specific operation.

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

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

[0140] Unless otherwise explicitly stated, as will be apparent from the following discussion, throughout the specification, discussions using terms such as "receive", "determine", "compare", "calculate", "average", "identify", "update", "resolve / solve", "output", 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 registers and memory into other data similarly represented as physical quantities within the computer system memory or registers or other such information storage, transmission, or display devices.

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

Claims

1. A computer-implemented method for determining physical properties of a subsurface target volume, the method comprising: Receiving a plurality of signals detected by a plurality of receivers disposed on a surface above the subsurface target volume, wherein each respective signal of the plurality of signals is detected by a respective one of the plurality of receivers; Performing cross-correlation on the signals between the plurality of receivers to obtain empirical travel time data of surface waves between a plurality of virtual source-receiver pairs; and Determining a model of the physical properties of the subsurface target volume, the model comprising a first plurality of cells arranged in a first two-dimensional grid, wherein each cell has a respective physical property value, and wherein the two-dimensional grid at least spans a region of the surface, wherein determining the model comprises: Selecting a subset of the plurality of virtual source-receiver pairs, wherein the surface wave path between each virtual source and the respective receiver passes through two or more of the first plurality of cells; Performing tomographic inversion on the empirical travel time data between each virtual source-receiver pair in the subset to obtain the physical property value of each cell.

2. The computer-implemented method according to claim 1, wherein, Performing the tomographic inversion comprises: Obtaining an initial velocity model that includes initial values indicating the velocity of each of the first plurality of cells, For each virtual source-receiver pair in the subset, determining a modeled travel time using the initial first velocity model values associated with the two or more cells through which the surface wave path passes; Determining a first error value that indicates the difference between the modeled travel time of each virtual source-receiver pair in the subset and the corresponding empirical travel time; and Determining an updated velocity model based on the error value, wherein the updated velocity model includes updated values indicating the velocity of each of the first plurality of cells.

3. The computer-implemented method according to claim 2, wherein, Iterating the determining of the modeled travel time, the determining of the first error value, and the determining of the updated velocity model based on the error value until the error value satisfies a first end condition, and wherein the updated velocity model of each iteration is used as the initial velocity model of the next iteration.

4. The computer-implemented method according to claim 2 or 3, wherein, The determining of the modeled travel time, the determining of the first error value, and the determining of the updated velocity model based on the error value are performed using iterative non-linear least squares inversion.

5. The computer-implemented method according to any one of claims 2-4, wherein, Performing the tomographic inversion further comprises: For each of the first plurality of cells, performing an inversion of the value indicating velocity to obtain the physical property value of the cell.

6. The computer-implemented method according to claim 5, wherein, Performing the inversion comprises: Obtaining an initial model of the physical properties of the subsurface volume, wherein the initial model includes initial physical property values of each of the first plurality of cells; For each of the first plurality of cells, determining a modeled value indicating velocity based on the initial physical property value of the cell; Determining a second error value that indicates the difference between the modeled value representing velocity and the corresponding value representing velocity, the corresponding value representing velocity being obtained from the updated velocity model; and Determine an updated model of the physical properties of the subsurface volume based on the second error value, wherein the updated model of the physical properties of the subsurface volume includes an updated physical property value for each of the first plurality of cells.

7. The computer-implemented method according to claim 6, wherein, Iterate between determining the modeled value indicative of velocity, determining the second error value, and determining an updated model of the physical properties of the subsurface volume based on the second error value until the second error value meets a second end condition, and wherein the updated model of the physical properties of the subsurface volume for each iteration is used as the initial model of the physical properties of the subsurface volume for the next iteration.

8. The computer-implemented method according to claim 6 or 7, wherein, Determining the modeled value indicative of velocity, determining the second error value, and determining an updated model of the physical properties of the subsurface volume based on the second error value are performed using iterative non-linear least squares inversion.

9. The computer-implemented method according to claim 1, wherein, Performing the tomographic inversion includes: Obtain an initial model of the physical properties of the subsurface volume, wherein the initial model includes an initial physical property value for each cell of the first plurality of cells; Use the empirical travel time data to determine an empirical dispersion function for each virtual source-receiver pair in the subset; Use the initial physical property values associated with the two or more cells through which the surface wave travels from the virtual source to the receiver to determine a modeled dispersion function for each virtual source-receiver pair of the subset; Determine a third error value that indicates the difference between the modeled dispersion function and the empirical dispersion function for each virtual source-receiver pair; and Determine an updated model of the physical properties of the subsurface volume based on the third error value, wherein the updated model of the physical properties of the subsurface volume includes an updated physical property value for each of the first plurality of cells.

10. The computer-implemented method according to claim 9, wherein, Determining the modeled dispersion function for each virtual source-receiver pair includes: Averaging the physical property values of the two or more cells through which the surface wave travels from the virtual source to the receiver; Calculating the modeled dispersion function using the averaged physical property values; Or: Calculate a cell dispersion function for each of the two or more cells through which the surface wave travels from the virtual source to the receiver to provide a plurality of cell dispersion functions, wherein calculating each cell dispersion function includes using the physical property value of the corresponding cell; and Average the plurality of cell dispersion functions.

11. The computer-implemented method according to any one of the preceding claims, wherein, The path taken by the surface wave between each virtual source and receiver pair is represented by a straight ray or a curved ray.

12. The computer-implemented method according to claim 11 when dependent on claim 10, wherein, The averaging includes determining a weighted average according to the respective weights of each of the two or more cells, wherein the respective weight of each cell corresponds to the segment length of the ray path in that cell.

13. The computer-implemented method according to any one of claims 1-10, wherein, The path taken by the surface wave between each virtual source and receiver pair is represented by a Fresnel zone.

14. The computer-implemented method according to claim 12 when dependent on claim 10, wherein, The averaging includes determining a weighted average according to the respective weights of each of the two or more cells, wherein the respective weight of each cell corresponds to the sensitivity value of the Fresnel zone in that cell.

15. The computer-implemented method according to any one of claims 9 to 14, wherein, The updated model for determining the physical properties of the subsurface includes updating the initial physical property values of each cell according to the respective weights of the paths taken by the surface waves in each cell.

16. The computer-implemented method according to any one of claims 9 to 15, wherein, Determining the updated model includes: determining the partial derivatives of the third error value with respect to one or more physical properties of the initial model; determining a sensitivity matrix including the partial derivatives; solving a system of linear equations defined by the sensitivity matrix to determine the required changes in the one or more initial physical property values of the initial model; and updating the initial model according to the required changes to produce an updated model of the physical properties of the subsurface volume.

17. The computer-implemented method according to any one of claims 9 to 16, wherein, Iterating the determination of the modeled dispersion function, the determination of the third error value, and the determination of the updated model of the physical properties of the subsurface volume based on the initial model until the third error value satisfies a third termination condition, wherein the updated model of each iteration is used as the initial model of the next iteration.

18. The computer-implemented method according to claim 17, wherein, The termination condition includes one or more of the following: the third error value of the most recent iteration is less than a predetermined threshold error value; and the difference between the third error value of the most recent iteration and the third error value of the previous iteration is less than a predetermined difference threshold.

19. The computer-implemented method according to any one of claims 9 to 18, wherein, Determining the modeled dispersion function, determining the third error value, and determining the updated physical model are performed using non-linear iterative least squares inversion.

20. The computer-implemented method according to any one of claims 6 to 19, wherein, The initial model of the physical properties of the subsurface volume provides initial physical property values in at least two spatial dimensions.

21. The computer-implemented method according to claim 20, wherein, The initial model of the physical properties of the subsurface volume provides initial physical property values in three spatial dimensions.

22. The computer-implemented method according to any one of claims 6 to 21, wherein, The initial physical property value and the updated physical property of each cell are the shear wave velocity as a function of depth.

23. The computer-implemented method according to claim 22, wherein, The method includes: for one or more cells of the updated model, calculating a compressional wave velocity value or depth profile, and / or a density value or depth profile using a shear wave velocity function.

24. The computer-implemented method according to any one of claims 9 to 23, wherein, The empirical dispersion function and the modeled dispersion function include a group velocity dispersion function and / or a phase velocity dispersion function.

25. The computer-implemented method according to any one of claims 6-10, 12 or 14, wherein, The method includes: outputting the updated model of the subsurface target volume to an output device.

26. The computer-implemented method according to any one of claims 6 to 25, wherein, The initial model of the physical properties of the subsurface volume is derived from empirical phase dispersion data obtained from signals detected between a second subset of the plurality of virtual source-receiver pairs.

27. A system, comprising: - one or more processors; - one or more memories storing computer-readable instructions configured to cause the one or more processors to perform the operations of a computer-implemented method including any of the preceding claims.

28. A computer-readable medium including instructions that, when executed by one or more data processing devices, cause the one or more data processing devices to perform the operations of a computer-implemented method including any one of claims 1-26.

Images