Tomography inversion with Gaussian bayesian priori with position dependent standard deviation

Abstract

Methods for determining a physical property of a subterranean target volume, as well as systems and computer-readable media for performing the methods, are provided. According to the method, clues are mined from geographic data, and the method further relates to the improvement of sustainability and environmental development: we create a safe and livable world together. The method includes performing tomography inversion on the obtained empirical propagation time data. Performing the tomographic inversion includes obtaining a posterior probability density function corresponding to each of the first plurality of cells, where obtaining the posterior probability density function includes obtaining a prior velocity model, and defining, for each receiver pair in the subset, a prior probability density function indicative of a velocity of each of the first plurality of cells, a modeled propagation time is determined using initial velocity model values associated with two or more cells traversed by the surface wave path, and a posterior probability density function of the velocity model is determined based on the modeled propagation time and a prior probability density function, wherein the prior probability density function depends on the distance of each cell from a first point on the surface.

Description

Technical Field

[0001] Exploring geographic data, this invention relates to improvements in sustainability and environmental development: together we create a safe and livable world. More specifically, this disclosure relates to methods and systems for analyzing target areas below the Earth's surface. In particular, this disclosure provides methods and systems for presenting models of subsurface target areas with higher precision. Background Technology

[0002] There is a persistent and widespread need for systems and methods for determining underground geotechnical parameters. In particular, there is a need for systems and methods capable of simulating the characteristics of target volumes below the surface to provide useful information for infrastructure planning and foundation design. Determining underground geotechnical properties early in the planning stages of a construction project reduces uncertainty during site selection, design, and construction. This, in turn, reduces delays, cost overruns, and unnecessary consumption of material resources (such as concrete) throughout the construction process and the asset's lifecycle.

[0003] A key parameter for determining the soil properties in a volume of interest is the shear modulus G and the shear velocity Vs. Shear velocity is the speed at which a shear wave passes through a material, controlled by the material's shear modulus. The relationship between shear velocity and shear modulus is determined by… Define, where It is the density of the material. Therefore, the measurement of shear rate provides valuable insights into the material properties of underground soil regions.

[0004] Surface wave spectral analysis (SASW) and surface wave multichannel analysis (MASW) are examples of techniques for collecting surface wave information that can be used to determine material properties in subsurface volumes. In both techniques, surface vibrations are measured from either a passive source (surface vibrations caused by environmental noise) or an active source (e.g., a falling heavy object), and the dispersion of the recorded surface waves is also recorded. ReMi (refractive micro-vibration) is another surface wave technique used to measure surface waves from environmental seismic noise recorded at the surface to infer material properties in subsurface regions.

[0005] Invasive downhole and inter-well techniques can also be used to determine the material properties of target subsurface areas. In both methods, a receiver located in the borehole measures recorded waves from an active source located elsewhere. In downhole techniques, one source or receiver is located subsurface within the borehole, while the other is located at the surface. In cross-hole techniques, the source is located in a first borehole, and the receiver is located in a second borehole. In both techniques, the propagation of the recorded waves is studied to infer the properties of the material through which the waves originating from the source pass.

[0006] Invasive techniques used to measure the material properties of subsurface areas often present logistical challenges, particularly in urban or hard-to-access environments, and the amount of data obtained using these techniques alone is often prohibitively expensive. These techniques allow for one-dimensional or some two-dimensional screening of subsurface properties, requiring significant resources (equipment and personnel) and are associated with safety risks and negative environmental impacts. Conversely, current surface techniques may lack the accuracy and reliability of more invasive analytical techniques. Summary of the Invention

[0007] According to a first aspect of this disclosure, a method for determining the physical properties of a subsurface target volume is provided. The method includes receiving a plurality of signals detected by a plurality of receivers arranged on a surface above the subsurface target volume, wherein each corresponding signal of the plurality of signals is recorded by a corresponding receiver of the plurality of receivers; cross-correlating the signals of the plurality of receivers to obtain empirical propagation time data of surface waves between a plurality of receiver pairs; receiving at least one subsurface signal detected by at least one subsurface receiver at a first location within the subsurface target volume, wherein the first location defines a first point on the surface above the subsurface volume; and determining a model of the physical properties of the subsurface target volume, the model including a first plurality of cells, wherein each cell has a corresponding physical property value, and wherein the plurality of cells at least span a region of the surface. As used herein, “empirical” data or information means data or information obtained through experience, i.e., real signals collected by actual receivers located on the surface above the subsurface target volume. Determining the model includes selecting a subset of the plurality of receiver pairs, wherein the surface wave path between each source and corresponding receiver passes through two or more cells of the first plurality of cells, and performing tomographic inversion on the empirical propagation time data between each receiver pair in the subset to obtain physical property values ​​for each cell.

[0008] The receivers used in this method can be seismic detectors, accelerometers, seismographs, vibration sensors, and / or transducers. The receivers can collect data over a considerable period. For example, ambient seismic noise can be continuously measured over five days. This extended recording time allows for the full extraction of surface wave information from ambient seismic noise recorded on or near the surface of the target area. In turn, as described in more detail below, the (processed) surface wave information can be used for tomographic inversion to obtain a subsurface shear wave velocity model.

[0009] Performing tomographic inversion involves obtaining a posterior probability density function corresponding to each of the first plurality of cells. Obtaining the posterior probability density function includes obtaining an initial velocity model for each of the first plurality of cells. For each receiver pair in the subset, the propagation time for modeling is determined using the initial velocity model values ​​associated with two or more cells traversed by the surface wave path. The posterior probability density function of the velocity model is then determined based on the modeled propagation time and the prior probability density function, where the prior probability density function depends on the distance of each cell from a first point on the surface. By performing tomographic inversion on empirical propagation time data from multiple receiver pairs in this manner, a three-dimensional model of the physical properties of the subsurface target volume can be obtained. This particular approach also considers the positional dependencies between points, resulting in a smoother, more realistic model that (relatively) excludes physically impossible jumps in geophysical properties. This approach also allows for the generation of higher-resolution models without increasing computation time to infeasible levels. By using a prior probability density function with a distance-dependent standard deviation, the tomographic inversion imposes a probabilistic preference on the subsurface characteristics of the intrusive measurement.

[0010] In some examples, methods for determining in-situ physical property values ​​of a subsurface target volume include recording at least one subsurface signal detected by at least one subsurface receiver at a second location within the subsurface target volume, wherein the second location is further defined by a second point on a surface above the subsurface volume, and wherein a prior probability density function depends on the distance of each model cell from the second point on the surface. In some examples, the method further includes receiving multiple subsurface signals corresponding to multiple locations, each location further defined by a point on the surface, wherein a prior probability density function depends on the distance of each cell from the multiple points on the surface. This forces a probabilistic preference for each in-situ physical property value for each measurement location.

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

[0012] In some examples, performing tomographic inversion further includes iteratively perturbing the initial velocity model to form one or more perturbations of the initial velocity model, selecting a perturbation as the updated velocity model, determining the likelihood of the updated velocity model based on the modeled propagation time and obtained empirical propagation time data, determining the posterior probability based on the likelihood of the updated velocity model and the prior probability density function, and calculating the acceptance probability based on the posterior probability of the updated velocity model and the posterior probability of the initial velocity model. In some examples, before determining the likelihood of the updated velocity model, performing tomographic inversion includes, for each receiver pair in the subset, determining the modeled propagation time using the updated velocity model values ​​associated with two or more cells traversed by the surface wave path, or determining the updated velocity model values ​​from the initial velocity model based on the previously calculated surface wave path. For example, the steps of iteratively perturbing the initial velocity, selecting the updated velocity model, determining the modeled propagation time, determining the likelihood of the updated velocity, and calculating the acceptance probability continue until a first termination condition is met. The updated velocity model from each iteration can be used as the initial velocity model for the next iteration. In some examples, determining the acceptance probability involves accepting the updated model if its posterior probability is greater than that of the initial model. In other examples, determining the acceptance probability also involves obtaining the ratio of the updated velocity model's posterior probability to the initial velocity model's posterior probability, generating a random number between 0 and 1, and accepting the updated model if the ratio is greater than the random number, and reverting to the initial model if the ratio is less than the random number.

[0013] In some examples, the steps of perturbing the initial velocity model, selecting an updated velocity model, determining the likelihood of the updated velocity model, and calculating the acceptance probability are performed using a reversible skip Markov chain Monte Carlo algorithm, and the cells are Voronoi cells, initialized with a random number of cells when the initial velocity model is obtained. In some examples, the perturbation can be one or more of the following: velocity changes of cells in a plurality of cells, cell position changes of a plurality of cells, birth of a cell core, death of a cell core, change of a first noise parameter; or change of a second noise parameter.

[0014] In some examples, the prior probability density function is defined by the following equation: , Among them, the prior covariance matrix It is a diagonal matrix, in which the trajectory is formed by the first... i The standard deviation of each model unit is composed of ,matrix N m Including the i The normalization factor of the probability density function of each unit, a vector V Composed of all model unitsn The earthquake parameter values ​​are composed of vectors. This includes measured seismic parameter values. This addresses the problem of how to improve upon the uniform prior used to infer the velocity distribution describing the subsurface velocity structure, since the Gaussian property of the probability density function imposes a probabilistic bias on specific velocity values. Therefore, this reflects a greater grasp of the specific values ​​for a particular region of the model. In some examples, the standard deviation of the Gaussian prior probability density function of the element increases to a threshold L as the distance from the element to the first point on the surface above the subsurface volume increases. This addresses the problem of how to more accurately generalize the influence of measured points on the model, since the parameter L is defined such that the model captures the fact that the first surface point has no influence after a certain distance.

[0015] In some examples, the standard deviation of the prior probability density function includes a scaling factor, where the scaling factor represents the rate at which the influence of the seismic parameter values ​​measured at a first point on the surface decreases with distance. This addresses the problem of how to interpret changes in the velocity structure of the subsurface volume, because using a scaling factor alters the influence of the measurement at the measurement point—a large scaling factor causes the standard deviation of the prior probability density function to decrease rapidly from the measurement point, limiting the influence of the measured seismic parameter values ​​on the model velocity to locations close to that point. Decreasing the scaling factor increases the influence of the measurement point at greater distances. A very small scaling factor is equivalent to simulating a perfectly horizontal, layered subsurface structure. In some examples, the standard deviation of the prior probability density function is defined as... ,in It describes the location. The standard deviation of the Gaussian probability density function at the parameter values. The first point on the surface The standard deviation at that point, A scaling factor representing a deterministic rate, the influence of seismic parameter values ​​measured at the first point on the surface decreasing with distance at said rate, until... Increase to threshold L This solves the problem of how to reflect the uncertainty of velocity relative to the measurement position, because for a small standard deviation, the first probability distribution behaves like a Gaussian distribution with respect to the in-situ velocity measurement v. loc Centered on a point and moving gradually away from it, the standard deviation increases, resulting in a uniform distribution between two fixed boundaries, where the model reflects v. loc right Vi There are no boundaries to influence this. In some examples, the method includes outputting an updated model of the subsurface target volume to an output device.

[0016] According to another aspect of this disclosure, a method for determining the physical properties of a subsurface target volume is provided. The method includes receiving multiple signals detected by a plurality of receivers arranged on a surface above the subsurface target volume, wherein each corresponding signal of the plurality of signals is detected by a corresponding receiver of the plurality of receivers; cross-correlating the signals among the plurality of receivers to obtain empirical propagation time data of surface waves between a plurality of 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 corresponding physical property value, and wherein the two-dimensional grid at least spans a region of the surface. Determining the model includes selecting a subset of the plurality of receiver pairs, wherein the surface wave path between each source and a corresponding receiver traverses two or more cells of the first plurality of cells, and performing tomographic inversion on the empirical propagation time data between each receiver pair in the subset to obtain physical property values ​​for each cell. By performing tomographic inversion on the empirical propagation time data of the plurality of receiver pairs in this manner, a three-dimensional model of the physical properties of the subsurface target volume can be obtained. This is achieved by obtaining a one-dimensional model for each cell of the grid of cells, which collectively provide information on the three-dimensional physical properties.

[0017] In some examples, performing tomographic inversion includes obtaining an initial velocity model comprising initial values ​​indicating the velocity of each of a first plurality of cells; determining, for each receiver pair in a subset, a modeling propagation time using the initial first velocity model values ​​associated with two or more cells traversed by the surface wave path; determining a first error value representing the difference between the modeling propagation time and the corresponding empirical propagation time for each receiver pair in the subset; and determining an updated velocity model based on the error value, wherein the updated velocity model comprises updated values ​​representing the velocity of each of the first plurality of cells. In some examples, determining the modeling propagation time, determining the first error value, and determining the updated velocity model based on the error value are performed iteratively until the error value satisfies a first termination condition. The updated velocity model from each iteration is used as the initial velocity model for the next iteration. In some examples, iterative nonlinear least squares inversion is used to determine the modeling propagation time, determine the first error value, and determine the updated velocity model based on the error value, which is a computationally efficient inversion technique.

[0018] In some examples, performing tomographic inversion further includes inverting the value representing velocity for each of the first plurality of cells to obtain the physical property value of that cell. The inversion may include obtaining an initial model of the physical properties of the subsurface volume, wherein the initial model includes the initial physical property value of each of the first plurality of cells; for each of the first plurality of cells, determining a modeling value representing velocity based on the initial physical property value of that cell; determining a second error value, which represents the difference between the modeling value indicating velocity and the corresponding value representing 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, wherein the updated model of the physical properties of the subsurface volume includes the updated physical property value of each of the first plurality of cells. The determination of the modeling value representing velocity, the determination of the second error value, and the determination of the updated model of the physical properties of the subsurface volume based on the second error are performed iteratively until the second error satisfies a second termination condition, and wherein the updated model of the physical properties of the subsurface volume in each iteration is used as the initial value of the physical properties of the subsurface volume in the next iteration. Computationally efficient least-squares inversion can be used to determine the modeling values ​​representing velocity, to determine the second error value, and to determine an updated model of the physical properties of the subsurface volume based on the second error value.

[0019] In an alternative example, performing tomographic inversion includes: obtaining an initial model of the physical properties of the subsurface volume, wherein the initial model includes initial physical property values ​​for each of the first plurality of cells; determining an empirical dispersion function for each receiver pair in the subset using empirical propagation time data; determining a modeling dispersion function for each receiver pair in the subset using initial physical property values ​​associated with two or more cells traversed by surface waves traveling from the source to the receiver; determining a third error value indicating the difference between the modeling dispersion function and the empirical dispersion function for each receiver pair; and determining 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 updated physical property values ​​for each of the first plurality of cells.

[0020] In some examples, determining the modeling dispersion function for each receiver pair may involve averaging the physical property values ​​of two or more cells through which the surface wave propagates from the source to the receiver, and using the averaged physical property values ​​to calculate the modeling dispersion function. Alternatively, determining the modeling dispersion function for each receiver pair may include calculating the cell dispersion function for each of the two or more cells through which the surface wave travels from the source to the receiver, to provide multiple cell dispersion functions. Calculating each cell dispersion function involves using the physical property values ​​of the corresponding cell and averaging the multiple cell dispersion functions. In some examples, the path of the surface wave between each source and receiver pair is represented by a straight line or a curve. In these examples, the averaging may include determining a weighted average based on a corresponding weight for each of the two or more cells, where the respective weight for each cell corresponds to the segment length of the ray path within that cell. In other examples, the path of the surface wave between each source and receiver pair is represented by Fresnel zones. In these examples, the averaging includes determining a weighted average based on a corresponding weight for each of the two or more cells, where the respective weight for each cell corresponds to the sensitivity value of the Fresnel zone within that cell.

[0021] In some examples, the model for determining the updated subsurface physical properties involves updating the initial physical property values ​​of each cell based on the corresponding weights of the paths taken by surface waves within that cell. For example, the changes in the physical property values ​​of each cell can be weighted as described above.

[0022] In some examples, determining the second model includes determining the partial derivatives of the third error value with respect to one or more physical properties of the first model; determining a sensitivity matrix including the partial derivatives; solving a system of linear equations defined by the sensitivity matrix to determine the expected changes in one or more initial physical property values ​​of the initial model; and updating the initial model based on the expected changes to produce an updated model of the physical properties of the subsurface volume.

[0023] In some examples, an updated model is iteratively executed, based on an initial model, to determine the dispersion function for modeling, a third error value, and the physical properties of the subsurface volume, until the third error value satisfies a third termination condition. The updated model from each iteration is used as the initial model for the next iteration. The termination condition may include one or more of the following: the third error value of the latest iteration is less than a predetermined threshold error value; and the difference between the third error value of the latest iteration and the third error value of the previous iteration is less than a predetermined difference threshold.

[0024] In some examples, nonlinear iterative least squares inversion is used to determine the dispersion function for modeling, to determine the third error value, and to determine the updated physical model.

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

[0026] In some examples, the empirical dispersion function and the modeled dispersion function include the group velocity dispersion function and / or the phase velocity dispersion function. Alternatively, both phase and group data can be used in the tomographic inversion process described herein.

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

[0028] According to another aspect of this disclosure, a system is provided that includes one or more processors and one or more memories, the memories storing computer-readable instructions configured to cause the one or more processors to perform any of the methods disclosed herein.

[0029] According to another aspect of this 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. Attached Figure Description

[0030] Embodiments of this disclosure will now be described by way of example to illustrate various aspects of this disclosure, with reference to the accompanying drawings, wherein: Figure 1 A cross-sectional view of the underground target volume 100 is shown, which shows the surface point 110 that defines the location of the in-situ velocity measurement. Figure 2 Multiple receivers are shown arranged on a surface above the underground target volume; Figure 3 This is a perspective view of the shear wave velocity model; Figure 4A It is a top view of the straight ray path on the surface above the underground target volume, with a shear wave velocity model superimposed on it; Figure 4B It is a top view of the curved ray path on the surface above the underground target volume, with a shear wave velocity model superimposed on it; Figure 4C It is a top view of an elliptical Fresnel kernel, representing the surface wave path above the underground target volume, with a shear wave velocity model superimposed on it; Figure 4D It is a top view of the curved Fresnel kernel, representing the surface wave path above the underground target volume, with a shear wave velocity model superimposed on it; Figure 5This is a schematic diagram illustrating a method for determining one or more physical properties of an underground target volume; Figure 6 This is a schematic diagram of a method for determining the surface wave velocity model; Figure 7 This is a schematic diagram illustrating a method for determining one or more physical properties of an underground target volume; Figures 8A-8C The standard deviation is shown plotted on a surface. The surface includes surface point 110 and additional reference point 112, at which no measurement data is available. Figures 8A to 8C The standard deviations for scaling factors of 160, 50, and 25 are shown respectively.

[0031] Figure 9 Examples of shear wave velocity models obtained by the method described in this paper; and Figure 10 This is a schematic diagram of a computing device suitable for performing the methods described herein. Detailed Implementation

[0032] For detailed explanations, please refer to the following: Figure 1 and Figure 2 A method for measuring the structural properties of soil volume using a receiver is described. Next, refer to... Figure 3-9 A novel method for tomographic inversion using signals detected by a receiver is disclosed. Finally, refer to... Figure 10 A computing device that can be used to perform the disclosed methods is described.

[0033] The following examples are described in the context of a receiver array on the surface of a soil volume and one or more underground receivers to aid understanding. However, it should be understood that the disclosed systems and methods are applicable to a wide variety of receiver types, including but not limited to receivers, accelerometers, velocimeters, seismometers, vibration sensors, and / or transducers. The disclosed methods can be applied to any suitable signal set.

[0034] The methods and systems disclosed herein generally relate to processing signals detected by a receiver on a surface. In a particular example, the receiver is placed on the ground, and the detected signal is environmental seismic noise. Processing these signals provides useful insights into the volumetric structure of surface and subsurface targets on which the receiver is placed, as described in more detail below. The inventors have found that the accuracy and / or resolution of information about subsurface target volumes can be improved by supplementing the signals detected by the receiver on the surface with in-situ measurements that provide additional empirical data. However, it is necessary to correctly represent the impact of the in-situ measurements on cells at different distances from the measurement location. The methods and systems disclosed herein provide a method for performing processing in a computationally efficient and practically feasible manner to produce accurate, high-resolution information about subsurface target volumes.

[0035] Before delving into the details of the disclosed method, some background information related to using receivers to determine surface and subsurface properties will first be provided. Shear modulus is a measure of a material's elastic shear stiffness, representing the deformation of a solid when subjected to a force parallel to one of its surfaces while its opposite surface is subjected to an opposing force. This force and its effect in a target subsurface volume are important parameters studied before and during the design of building and infrastructure projects. To determine the shear modulus of a volume, the shear rate Vs needs to be determined. This, in turn, indicates the stiffness of the subsurface material and its ability to support structures located above and / or through the volume.

[0036] In the context of land studies, two types of waves are typically distinguished: P-waves, in which particles in a volume oscillate along the direction of wave propagation, causing compression and decompression of the land as the wave propagates; and S-waves, which are shear waves, in which particles oscillate in a direction perpendicular to the wave propagation direction.

[0037] P-waves and S-waves are volume waves that propagate in all directions along the bulk of a volume. The interaction of P-waves and S-waves with the Earth's surface produces surface waves that propagate along the Earth's surface. Several types of surface waves can be distinguished. Rayleigh waves are measured and studied in the system and method described herein because it is convenient to measure the vertical component of surface vibrations. However, it is understood that other surface waves (such as Love waves and Scholter waves) can be measured and utilized in the system and method described herein.

[0038] Because surface waves propagate in two dimensions (on the surface), they decay more slowly than volume waves (which propagate in three dimensions). Surface waves typically exist at a depth of one wavelength above the Earth's surface, generally propagate slowly, and have frequencies significantly lower than volume waves.

[0039] The lower attenuation, slower propagation time, and lower frequency of surface waves make their study particularly attractive for determining the shear velocity Vs. Because of the weaker attenuation, the signal strength is better maintained over longer propagation distances. Therefore, the resulting measurements typically have higher signal quality (signal-to-noise ratio) than those obtained from volume wave studies.

[0040] See now Figure 1 The diagram shows a cross-sectional view of a subsurface volume 100. P-waves and S-waves, as volume waves, pass through volume 100. A surface 102, defined by the x and y planes, extends above the subsurface volume. Surface waves propagate along surface 102. In some examples, the subsurface volume 100 may include one or more boreholes, for example, borehole 106 extending from surface 102 through the subsurface volume 100 (i.e., extending in the z-dimensional). Of course, it is understood that in practice, borehole 106 may extend through both the x and y dimensions in addition to the z-dimensional, but for illustrative purposes, boreholes extending only along the z-dimensional are discussed.

[0041] At one example point, a schematic representation of particle oscillations (caused by Rayleigh wave propagation) at the surface above the target subsurface volume is shown. As shown, the oscillations of particle P are partly vertical and partly parallel to the propagation direction. Therefore, the resulting particle motion is essentially elliptical.

[0042] A plurality of receivers 104 are arranged on surface 102 above volume 100. The receivers 104 located on surface 102 can be configured to measure the vertical component of the oscillation schematically shown at this example point.

[0043] Receivers 104 are arranged in a grid array on the surface, extending in two directions. It should be noted that in many cases, the surface above the target region may not be planar. Therefore, the receiver array 104 may not be truly “two-dimensional”, as each receiver may be offset relative to its neighbors in the z-axis. However, for simplicity, this grid arrangement of receivers will be referred to herein as a 2D array.

[0044] It should also be understood that, in the context of underground soil research, the term 2D array can also be used to refer to an array configured to capture 2D information (e.g., a row receiver configured to provide information about a 2D slice of a target area), and 3D array can refer to a receiver grid configured to capture 3D information. However, in the context of the information being described, the term 2D in the receiver context is intended to refer to the two-dimensional configuration of the receiver, rather than the information being collected.

[0045] The surface wave propagating along surface 102 will cause vertical motion at multiple receivers 104 as the wave passes through the surface.

[0046] To determine the shear velocity Vs from observations of surface waves (especially Rayleigh waves), the dispersive behavior of the surface waves can be measured. Surface waves are dispersive, meaning their velocity depends on frequency. Generally, seismic wave velocities increase with depth within the Earth. Therefore, normal surface wave dispersion indicates that surface wave velocities decrease with increasing frequency. By studying the behavior of surface waves on the surface above a subsurface target volume, the material properties of that volume can be determined.

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

[0048] The group velocity of a wave is the speed at which the overall envelope shape of the wave amplitude (called the wave modulation or envelope) propagates in space. Group velocity is equivalent to the speed at which the wave's energy propagates within a volume. It is measured by determining the wave propagation between (synthetic) transducer pairs and is a frequency-dependent point property within a volume, which in turn depends on depth. Group velocity is obtained by measuring the time of flight between a (virtual) source and receiver.

[0049] Phase velocity is the speed at which a specific frequency component of a wave propagates. Therefore, phase velocity is expressed as a function of frequency. To measure phase velocity, at least two measurement nodes are selected to measure the wave propagating through a volume to determine the relative flight time of different frequencies between the receivers. The result is the average phase velocity as a function of frequency over the volume between the two measurement nodes. Phase velocity is obtained as a point in 2D phase space (dispersion spectrum), which is obtained through a 2D transformation (such as tilt superposition, Radon, FK, etc.) of an array of recorded waveforms (time-distance space).

[0050] Refer again Figure 1 Each receiver 104 provides a measurement node for measuring the vertical component of the passing surface wave. Receiver 104 can be configured to measure vibrations caused by environmental seismic noise. That is, the background wave field caused by natural or anthropogenic noise (rather than the pulse points used in active methods, such as explosions or falling hammers).

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

[0052] In some examples, one or more underground receivers are located at one or more subsurface locations below the surface 102. In some examples, an underground receiver may be defined by a single surface point (i.e., only one underground receiver exists at a given x, y coordinate).

[0053] In some examples, where the subsurface volume includes a borehole, one or more subsurface receivers may be arranged along the entire length of the borehole 106 within the volume 100. In this case, multiple receivers 108 are arranged to provide an in-situ velocity distribution by obtaining multiple velocity measurements defined by a point on the surface 110 (i.e., obtaining an in-situ distribution in the z direction at a given x, y coordinate using multiple receivers).

[0054] The underground receiver 108 can be configured to measure underground signals from which field empirical time-propagation data can be obtained. The receiver 108 can be any type of field instrument used to obtain field empirical time-propagation data. For example, field velocity values ​​can be obtained through various techniques, such as geotechnical engineering (seismic CPT) or geophysical (cable logging) invasive methods.

[0055] Figure 2 Multiple pairs of virtual receivers are shown on the surface above the underground target volume. Ray paths 202 between receiver pairs are indicated, specifically the ray paths from a single central receiver near the center of the receiver array and from each other receiver in the array. Background shading and outline rings represent the propagation time field from the central receiver to the other receivers. Corresponding ray paths are also shown between each receiver and all other receivers, i.e., between each pair of receivers; for simplicity, Figure 2 It is not displayed.

[0056] Each receiver pair can provide a signal that is cross-correlated with the signal at the second position from the first position, so as to reconstruct the virtual receiver pair using the principle of interferometry. Specifically, in Figure 1 The cross-correlation of the passive noise measured at the corresponding receiver at the surface can be used to reproduce the response from the subsurface target volume as if it were caused by an impulse point source equivalent to a Green's function.

[0057] In some examples, Figure 1 Receivers 104a and 104b are formed in an array at surface 102 as a first receiver pair. By cross-correlating the signals received at receivers 104a and 104b, receivers 104a and 104b can act as (virtual) source-receiver pairs, where each receiver in the pair records a signal as if the signal originated from the other in the pair.

[0058] In other words, the response received by cross-correlating the records of two receivers can be interpreted as a response measured at one receiver location, as if there were a source at the other receiver. Various methods for determining the Green's function for a virtual receiver pair are known, and the following article provides an overview of each of these methods, in “Seismic Interferometry Tutorial: Part 1 - Basic Principles and Applications”; Geophysics, Vol. 75, No. 5 (September-October 2010; pp. 75A1950-75A209; Wapenaar et al.).

[0059] exist Figure 1 In the array, only one pair of receivers is labeled (104a, 104b). However, it should be understood that for each receiver 104 in the array, every other receiver in the array can act as the other half of the source-receiver pair. In this way, the Green's function for each source-receiver pair can be obtained. The Green's function over multiple dummy source-receiver pairs was investigated to determine the dispersive behavior of the surface waves.

[0060] See Figure 3 The first shear wave velocity model 300 has a cell grid m on the surface above the underground target volume. xy The grid consists of m columns extending along the x-direction. x1 m x2 m x3 Equal to and extending along the y-direction m 1y m 2y m 3y Each cell defines a region of the surface. For example, each cell could define a 5m by 5m square; other example options include a 1m by 1m square or a 10m by 10m square. In other words, the shear wave velocity model comprises multiple cells arranged in a two-dimensional grid 300. This two-dimensional grid covers at least the surface region above the subsurface target volume. Choosing smaller cell regions improves the model's resolution. Each cell also includes a volume extending vertically below the surface region. Model 300 can extend infinitely below the surface region or to a predetermined depth below the surface where the shear wave velocity value significantly affects the propagation of surface waves on the surface. For example, model 300 can be defined to depths up to 50m, 100m, 200m, or 300m.

[0061] In some examples, the first shear wave velocity model 300 is not a mesh of square cells, but rather includes cells with rectangular, rhomboid, or other spliced ​​shapes, including non-uniform shapes or combinations of different shapes. In some examples, the first shear wave velocity model may include spliced ​​Voronoi cells, and a random number of cells are initiated when method 700 is executed to determine the velocity of each cell. This will refer to... Figure 7To provide a more detailed description.

[0062] Each element is associated with a shear wave velocity value, which is an example of a physical property value representing the expected shear wave velocity within the actual subsurface target volume of interest. The shear wave value carries depth information, either because it is constant throughout the entire volume below the element region, or because it indicates how the shear wave value varies with depth (…). Figure 3 The shear wave velocity varies (in the z-direction). For example, the shear wave value defined for each cell can be explicitly a function of depth, either a continuous function or a series of values, each with an associated depth range. In another example, the shear wave velocity value can be a function of frequency, which corresponds to depth information because surface wave propagation is influenced by the physical properties of the subsurface volume up to approximately one wavelength deep. In other words, low-frequency surface waves are influenced by physical properties at deeper depths than high-frequency surface waves.

[0063] In other examples, the model can be a physical property other than shear wave velocity, such as compressive wave velocity, density, elastic modulus, shear modulus, or, if a viscoelastic model is used, optionally a viscosity quality factor Q. s and Q p Generally, a model can define multiple physical property values.

[0064] See Figures 4A to 4D This describes various representations of the wave path between a (virtual) source and receiver. The wave path is the route taken by a surface wave from the source to the receiver. See below for reference. Figure 4A and 4B The wave path can be represented by a straight line or curve connecting the source and receiver. Alternatively, refer to the following... Figure 4C The wave path can be represented by an elliptical Fresnel region (also known as the Fresnel nucleus) between the source and receiver pair. As an alternative, see the following reference... Figure 4D The wave path can be represented by a curved (“banana”) Fresnel zone.

[0065] See Figure 4A The ray path is defined as a straight line between two receivers at surface locations A and B on model 300. The straight-line ray path model assumes a constant lateral velocity (in other words, the velocity varies only with depth). In many cases, the straight-line ray path model provides a sufficiently accurate approximation of the actual path traversed by the wave between the source and receiver. Compared to other more complex techniques, the straight-line path approximation is also mathematically and computationally efficient. According to ray theory, the ray path represents the motion of a surface wave propagating from A to B (or from B to A). Specifically, the ray path is defined as the direction of propagation of the surface wave, either perpendicular to the wavefront in wave theory or perpendicular to the propagation-time profile. Figure 4AThe ray path shown passes through seven elements of model 300, numbered 1 to 7. This ray path consists of seven ray segments, each with a length L. i It depends on the length of the ray path through each unit, i.e. Figure 4A L1 to L7 in the middle.

[0066] See Figure 4B The ray path is defined as the curve between the two receivers at surface locations C and D. Figure 4B The ray path shown passes through six elements of model 300, numbered 1 to 6. This ray path comprises six ray segments, each with a length L. i It depends on the length of the ray path through each unit, i.e. Figure 4B L1 to L6 in model 300. Based on the shear wave velocity values ​​of the elements in model 300, the trajectory of the ray path can be calculated using the minimum time principle (Fermat's principle) between positions C and D.

[0067] See Figure 4C The wave path is represented by an elliptical Fresnel region between two receivers at positions E and F on the surface. The Fresnel region is a geometric ray (e.g., along the gradient direction from receiver (F) to source (E) – the spatial derivative of the propagation time)) Figure 4A The area around the straight line shown.

[0068] More specifically, those skilled in the art will understand that any wave propagating along the path between the (virtual) source at E and the receiver at F will have some off-axis propagation component (not along the straight line connecting E to F). This off-axis propagation component may be deflected by the wave propagation medium, and some of the deflected wave is then directed to the receiver at F. Therefore, the travel time of a wave traveling along the direct (straight) path will differ from that of a wave traveling along the deflected path, meaning that the straight-path wave and the deflected-path wave will arrive at the receiver out of phase. When the phase difference is half the wave period (or one and a half wave periods, or two and a half, or any odd integer multiple of half), the phase difference may cause destructive interference. Conversely, when the phase difference is between 0 and half the wave period (or between 1 and half, or between any integer n-1 and n-½ wave periods), the phase difference will cause constructive interference. In other words, the nth Fresnel zone is defined as a region in which a wave deflected at a point has a phase difference between n-1 and n-½ wavelengths compared to a wave traveling in a straight line between the source and receiver pair. In the examples of this disclosure, only the first-order Fresnel zone can be considered for the wave path between the source and receiver, since the wave components in higher-order Fresnel zones become very small. BD Guenther (2005)'s "Encyclopedia of Modern Optics" describes the general physical principles of the Fresnel zone of any wave propagating through any medium, whereby those skilled in the art will understand that the description of the Fresnel zone applies to any propagation medium. Shibo Xu and Alexey Stovas (2018) "Fresnel Zones in VTI and Orthogonal Media" Fresnel zone in VTI and orthorhombic media The Fresnel zone also describes the waves that pass through the volume of an underground target.

[0069] like Figure 4C As shown, this figure represents a first-order Fresnel region approximated by a straight line. This Fresnel region is elliptical, with the (virtual) source location E and the receiver F being the foci of the ellipse. In other examples, such as... Figure 4D As shown, the Fresnel region may be a curved, so-called banana-shaped region, following the path of a curved ray between the virtual source G and the receiver H. Determining this curved Fresnel region follows the same principles as described above. Figure 4A The same considerations apply to the elliptical Fresnel zone, but the differential propagation time and phase difference are relative to the curved ray (see [link to Fresnel zone description]). Figure 4B It is determined by a straight ray between the source and the receiver, rather than by a straight ray between them.

[0070] Figure 4C The Fresnel zone depicted covers eight units of model 300, numbered 1 to 8. This Fresnel zone includes those labeled as... R 1 to R 8The eight subregions are each Fresnel regions within corresponding units 1 to 8. Each Fresnel subregion... R 1 to R 8 A corresponding sensitivity value is calculated uniquely for each Fresnel zone. As those skilled in the art will understand, the sensitivity value is the sensitivity to changes in local shear wave velocity (or any other physical property being modeled). This technique for calculating the sensitivity value for each subregion of the Fresnel zone (i.e., the sensitivity value of the Fresnel zone for each cell covered by the Fresnel zone) is well known to those skilled in the art.

[0071] Figure 4D The Fresnel zone depicted in the diagram also shows the same illustration of a sub-region that passes through units 1 to 8, and is composed of sub-regions. R 1 to R 8 Composition (and) Figure 4C The sub-regions depicted in the text are different. Figure 4D Each sub-region depicted has its own sensitivity value, which those skilled in the art will understand may be related to... Figure 4C The sensitivity values ​​for each sub-region depicted are different.

[0072] For each of the multiple virtual receiver pairs, such as Figure 2 As shown, the corresponding Fresnel regions can be determined. Each Fresnel region consists of multiple sub-regions covered by the Fresnel region for each cell, and each sub-region of each Fresnel region has its own sensitivity value.

[0073] See Figure 5 A method 500 for determining one or more physical properties of an underground target volume includes receiving 502 multiple signals detected by receivers arranged on a surface, such as reference signals. Figure 1 As described above. Typically, one signal is received from each receiver, although in some cases only a subset of signals from all receivers may be received. This signal represents the vertical oscillations of the soil measured at each receiver, caused by environmental seismic noise or surface waves from active sources, as referenced above. Figure 1 The signal may include metadata about the location and time of the recording.

[0074] The signal can be received directly from the receiver or through one or more intermediate devices. For example, a computing device (see reference below). Figure 10The signal (502) can be located locally on the receiver array for transmitting signals via any form of wired or wireless communication. Alternatively, the signal can be received at a location remote from the receiver array to remotely perform method 500, for example, by communicating via the Internet or transmitting a physical computer-readable medium on which a record of the signal is stored. In some examples, additional processing steps can be performed on the signal before or after receiving the signal 502 to optimize the signal for further processing.

[0075] Examples of surface receivers may include speedometers or accelerometers. A specific mechanism of such a receiver includes a ferromagnetic material on a spring that moves within an electric coil in response to surface motion, thereby inducing a current proportional to the measurable ground velocity.

[0076] Examples of underground receivers could include accelerometers or piezoelectric (pressure) sensors.

[0077] Method 500 further includes step 504, cross-correlating the signals between a plurality of receivers arranged on a surface to obtain empirical propagation time data of surface waves between a plurality of virtual receiver pairs. As described above, this step involves cross-correlating passive noise signals measured at a pair of receivers to obtain a correlated Green's function, which represents the wave field as if one of the pair were a virtual source and the other a receiver. Those skilled in the art will well understand that information indicating phase velocity as a function of frequency and information indicating group velocity as a function of frequency can be derived from the cross-correlation process. In step 504, cross-correlation of the signals between the plurality of receivers is performed to obtain group and / or phase propagation time data for a plurality of frequencies. The frequencies may be continuous to obtain group velocity and / or phase velocity dispersion information, or propagation time data for a plurality of finite frequencies may be obtained. Group propagation time refers to the time it takes for a group of waves to propagate from the virtual source to the receiver. Phase propagation time refers to the time it takes for a specific phase component of the wave to propagate from the virtual source to the receiver. Since the distance between the virtual source and the receiver is known, the group and / or phase velocities or slowness (which are the reciprocals of the velocities) for a plurality of frequencies can be derived from the correlated Green's function of the receiver pair. Slowness and propagation time are actually equivalent; the difference is a scaling factor provided by the distance between the source and receiver.

[0078] A model for the physical properties of an underground target volume can be determined based on at least multiple frequencies of group and / or phase propagation time data and in-situ velocity measurements determined by an underground receiver. (See above reference.) Figure 3The model comprises multiple cells arranged in a two-dimensional grid. Each cell in the model has its own physical property values ​​(such as shear wave velocity), which may include a depth distribution. Therefore, in general, the cells of the model provide a three-dimensional model of the physical properties, as each cell provides a one-dimensional model as a function of depth.

[0079] The model for determining the physical properties of the underground target volume involves selecting a subset of over 506 receiver pairs from which empirical propagation time information has been obtained. (See reference above.) Figures 4A-4D The wavepath between each receiver pair traverses two or more cells of the model. In principle, any and all receiver pairs may be suitable for determining the model. However, increasing the number of wavepaths used will increase the required computational power, otherwise leading to very long computation times. Therefore, in practice, it is often beneficial to select a subset of receiver pairs, the number of which to select depending on available computational power or other practical problems. Several principles can be used to limit the number of receiver pairs selected. First, due to the reciprocity of cross-correlation, the wavepath from position A to position B is the same as the wavepath from position B to position A, so only paired receivers need to be selected, regardless of which is considered a virtual source. One way to deselect receiver pairs is to omit those receiver pairs that lack suitable propagation time data, i.e., the cross-correlation of the signals from the pair of receivers does not show an identifiable wave propagating through both receivers that can be used to obtain the propagation time between them. Another way to select wavepaths is to omit pairs with low-quality pickup data, such as those with low signal-to-noise ratios, large uncertainties, or outliers compared to adjacent pickup data. Any pick-up data that would produce non-physical or non-geological results is excluded, as well as any pick-up data that is not reliable for the specific subsurface target volume being investigated. Finally, if the number of remaining pairs still exceeds the number that can be used, a subsample of receiver pairs can be selected, such as selecting one pair out of every n pairs or selecting randomly. In this case, other subsamples of suitable receiver pairs can be used successively to increase the total number of pairs used to exceed memory limits or available computing power. Regardless of the method used to limit the number of receiver pairs (if necessary), the end result is the selection of receiver pairs and corresponding wave paths for which empirical group and / or phase propagation time information has been obtained, such as velocity or dispersion functions at multiple different frequencies (e.g., group velocity dispersion function and / or phase velocity dispersion function).

[0080] As referenced above Figures 4A-4B The wavepath between the (virtual) source and receiver pair can be a ray path extending in a straight line between the source and receiver. In some examples, at least some of the ray paths extend along curves between the respective source and receiver, as referenced above. Figure 4BThe aforementioned method describes a ray path that will pass through at least two cells of the first model, having a ray segment length in each cell it passes through (including the cells where the ray path begins and ends). The curved ray path can be determined based on the initial physical property values ​​of the model's cells. By using curved ray paths determined according to the physical property values ​​of the model's cells, the ray path more closely approximates the actual ray path of the wave propagating between the source and receiver. Therefore, this method improves the accuracy of subsequent processing results because it more closely corresponds to the physical reality of surface waves passing through the subsurface target volume.

[0081] As referenced above Figures 4A-4B In other examples, wave paths are represented by Fresnel zones, which provide the area over which a wave takes many possible paths as it propagates from the source to the receiver. A Fresnel zone provides a region covering two or more cells through which possible waves can propagate, each cell having a corresponding Fresnel zone with an associated sensitivity value defined by the Fresnel zone. The Fresnel zone can be determined based on the initial physical property values ​​of the cells in the model. Because the Fresnel zone defines the area over which surface waves can propagate from the source to the receiver, this approach improves the accuracy of the final model for each cell, as it more closely corresponds to the physical reality of surface waves passing through and scattering within a subsurface target volume.

[0082] The model for determining the physical properties of the subsurface target volume includes step 508 of retrieving one or more subsurface seismic parameter measurements from one or more subsurface receivers 108 defined by a first point on the surface 110. For example, in-situ velocity values ​​may be obtained. In some examples, one or more subsurface receivers 108 may be inserted into the surface to obtain in-situ velocity values. v loc In some examples, one or more underground receivers 108 may be intrusive underground receivers. In some examples, one or more underground receivers 108 may be inserted into borehole 106 to obtain in-situ velocity values. v loc .

[0083] One or more underground receivers 108 can obtain in-situ velocity values ​​using any suitable known method, such as geotechnical engineering (seismic CPT) or geophysical (cable logging) invasive methods. v loc .

[0084] Method 500 also includes determining and obtaining the first point of in-situ velocity value. x loc , y loc corresponding standard deviation Step 510.

[0085] In some examples, standard deviation Defined based on the surface quality of a subsurface target volume of 100, to influence the distance from the first point on the surface, in-situ velocity measurements have a probabilistic influence on this. In some cases, the standard deviation... It is based on x loc , y loc The in-situ velocity value is obtained from the point. In a further example, the standard deviation... The values ​​were selected from a set of known experimental values ​​that the inventors found to produce good results. For example, the inventors found that a predetermined standard deviation of 0.01 to 0.1 produced good results. In some examples, the standard deviation was set to 0.04.

[0086] In step 510, the scaling factor is determined. (Compared to the standard deviation) Similarly, scaling factor It is defined based on the surface quality of a volume of 100, to influence the distance with a probabilistic effect on in-situ velocity measurements. Scaling factor The scaling factor can be defined in the range of 1 to 200, where a scaling factor of 1 can represent a perfect surface, i.e., on which the influence of the measurement point will never decrease. Refer to step 708 of method 700. Figures 8A to 8C The scaling factor is described in more detail.

[0087] Steps 508 and 510, and method 500, involve retrieving in-situ parameter measurements corresponding only to a first point on the surface. However, it should be understood that in-situ surface seismic parameter measurements can be obtained from more than one point, in which case the standard deviation and scaling factor corresponding to each point are obtained. For example, a second in-situ parameter measurement can be obtained by receiving at least one subsurface signal detected by at least one subsurface receiver at a second location within the subsurface target volume, wherein the second location defines a second point on the surface above the subsurface volume, and wherein the prior probability density function depends on the distance of each cell to the second point on the surface.

[0088] In some examples, multiple underground signals corresponding to multiple locations can be received, each location defining a point on the surface, where the prior probability density function depends on the distance of each cell to multiple points on the surface.

[0089] The model for determining the physical properties of the subsurface target volume also includes performing a 512-step tomographic inversion on the propagation time data for each of the selected receiver pairs, while constraining the inversion using the in-situ velocity values ​​obtained in step 508. Generally, tomographic inversion is used to derive the physical property values ​​for each cell of the model from the group or phase propagation time data obtained for each of the multiple receiver pairs. As mentioned above, the physical property values ​​for each cell can be a single value, a distribution, or physical properties provided as a function of depth. See below for reference. Figure 6-7 This describes in detail the various procedures involved in tomographic inversion. It is understood that a prior probability density function related to the measurement location is used as a reference below. Figure 7 The starting point of this method can be any suitable tomographic inversion process. As described throughout this document, tomographic inversion can be performed based on empirical group propagation time information or empirical phase dispersion information. In the case of phase dispersion, phase velocity or slowness information can be converted into equivalent propagation time information, thereby enabling tomographic imaging of this phase propagation time information.

[0090] In some examples, tomographic inversion step 512 involves Figure 6 The process 600 shown below, as described in more detail below, maps the propagation time information from each receiver pair to cells of a physical property model via tomography. Therefore, the result of this process is an empirical model of group velocity and / or phase velocity (depending on the propagation time data used), with each cell of the model having its own group velocity value or phase velocity value or distribution. This process 600 is performed for each of a plurality of frequencies to obtain the group velocity value or phase velocity value for each cell at each frequency.

[0091] Process 600 is the first stage (i.e., the tomographic stage) of the two-stage tomographic inversion. Process 600 includes the 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 cell at each of a plurality of frequencies. The initial velocity model sets initial (group or phase) velocity values ​​for each cell (for each given frequency), which process 600 refines using empirical propagation time data during iterations. Therefore, the initial velocity model and the corresponding initial velocity values ​​are not necessarily high-precision or high-resolution, although a more accurate initial model can make the tomographic process 600 faster or more accurate in mapping propagation time to each cell. A more accurate initial model can also reduce the risk of finding local minima rather than global minima in the iterative gradient descent method, although this problem can be addressed using Monte Carlo methods. In some examples, the initial velocity model is determined based on user input, such as historical data or map information indicating possible physical property values ​​on the subsurface target volume. Alternatively, arbitrarily chosen typical values ​​of the group / phase velocities can be used as the starting point for each cell. In these examples, an arbitrary model can be selected based on estimates of the physical properties of the subsurface target volume. The tomographic imaging process 600 also includes a step 604 of determining the modeled group or phase propagation time for each selected receiver pair using an initial group or phase velocity model. This is achieved by identifying the wave path from the source to the receiver (e.g., as a straight line, curve, or elliptical / surface Fresnel zone) and identifying the cells traversed by the wave path. The propagation time based on the initial velocity model is then determined based on the known distance between the source and receiver and the velocity value of each cell traversed by the wave path.

[0092] The tomographic imaging process 600 also includes a step 606, determining an error value that represents the difference between the modeling propagation time (determined in step 604) and the empirical propagation time (obtained in step 504) for each selected receiver pair. For example, the error value could be a simple difference between the modeling and empirical propagation time for each receiver pair at each frequency, also known as a residual. In some examples, the error value could be a combination of all differences between the modeling and empirical propagation time for each receiver pair. Generally, determining the error value is part of an iterative process; for example, in least-squares inversion methods, the square of the residual is calculated to minimize the square of the residual through iteration. Other forms of inversion processes use different error values ​​to provide feedback to the initial group or phase velocity model.

[0093] Process 600 also includes a step 608 of using error values ​​to determine an updated velocity model. Except for the new velocity values ​​associated with at least some cells, the updated model is substantially identical to the initial velocity model in all respects. In other words, the updated model is an updated version of the initial velocity model that takes into account determined error values ​​between the experience of each receiver pair and the modeled group or phase propagation time. This feedback process may involve least squares, Markov chain Monte Carlo, 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.

[0094] Steps 604 to 608 are typically all parts of the subroutine of the first-stage tomographic imaging process 600, which is then iterated according to an inversion method such as least squares inversion. Figure 6 The dashed arrows in the diagram illustrate the iterative nature of process 600, indicating that the updated model determined in step 608 is used to determine the propagation time of the new modeling for each receiver pair in step 604. In other words, each time an updated model is determined using the initial model for each receiver pair and the resulting propagation time of the modeling, the resulting updated model is then used as the initial model for the next iteration. Iteration continues until the error value reaches a termination condition, for example, the error value falls below a threshold absolute value of the difference between the modeling propagation time and the empirical propagation time, or falls below a threshold proportional difference between the modeling propagation time and the empirical propagation time. Another termination condition, which can be used alone or in combination with the threshold error value termination condition, is that the change to the initial model during the iteration to produce the updated model is less than a threshold amount or a threshold proportion. Thus, if the iteration reaches a set minimum error value, the iteration can end because further iterations will not significantly improve accuracy.

[0095] Once the tomographic imaging of the first stage process 600 has been performed to obtain a velocity model including group or phase velocity values ​​for each cell (each frequency), the tomographic inversion of step 512 can proceed to 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 of the inversion process 700 is as follows: Figure 7 As shown. Similar to the first-stage tomographic imaging process 600, a second-stage inversion process 700 is performed for each of the multiple frequencies to obtain the shear velocity value for each cell. The inversion process 700 is iterated until the number of iterations satisfies a first termination condition, and the updated velocity model of each iteration is used as the initial velocity model for the next iteration. For example, the first termination condition could be the number of hundreds to thousands of iterations performed within the region. Of course, any appropriate number of iterations can be determined before the iterative process 700.

[0096] 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, wherein the initial model includes initial physical property values ​​for each element in a first plurality of elements. The initial physical model may be as described in the reference above. Figure 3 The shear wave velocity model 300, and / or the model with reference above. Figure 3 The initial physical model defines any features or variations described. The initial model sets initial physical property values ​​for each cell of the model, and method 700 refines these initial physical property values ​​through further iterations using the group / phase velocity model obtained from the iterative tomography process 600. Therefore, the initial physical model and initial physical property values ​​are not necessarily high-precision or high-resolution models of the subsurface target volume, although a more accurate initial model can improve the expected accuracy of the final result or reduce the computation time to achieve the final result.

[0097] Process 700 further includes a step 704 of determining the modeled surface wave velocity for each cell based on the initial model. For each receiver pair in the subset, the modeled propagation time is determined using the initial velocity model values ​​associated with the two or more cells through which the surface wave path passes. More specifically, using a forward modeling approach, the phase velocity and / or group velocity of each cell in the initial physical model are derived using the corresponding physical property values ​​of a given cell. For example, the physical property value for each cell could be the shear wave velocity. The shear wave velocity value 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 developed by Thomson (1950) "Propagation of Elastic Waves in Layered Solid Media" (… Transmission of elastic waves through a stratified solid medium "", Journal of Applied Physics, 21(2), 89-93) and Haskell (1953), "Dispersion of surface waves in multilayer media ( The dispersion of surface waves on multilayered media This is introduced in the American Seismological Society Bulletin, 43(1), 17-34. This will be clear to those skilled in the art as well as other methods that can be used to derive model group or phase velocity (or dispersion function) from shear wave velocity values ​​(or depth functions).

[0098] The inversion process also includes a step 706 of determining an error value, which represents the difference between the modeling velocity (determined in step 704 based on the initial physical model) and the empirical velocity (obtained from process 600). The error value can be determined in a similar manner to that described above in conjunction with step 606, where an error value between the modeling propagation time and the empirical propagation time is determined. Similar to tomographic process 600, determining the error value in step 706 is generally part of an iterative process, for example, in a least-squares inversion method, calculating the square of the residuals in order to minimize the residuals through iteration.

[0099] Other forms of inversion processes use different error values ​​to provide feedback to the initial model of the physical properties of the subsurface target volume. For example, Bayesian inference can be used to determine the error values ​​to provide feedback on the initial model of the physical properties of the subsurface target volume provided in step 702. In these examples, the error value representing the difference between the modeling rate of each element of the model (determined in step 704 based on the initial physical model) and the empirical rate (obtained from process 600) might be a likelihood. Step 706 may include determining the likelihood of the modeling rate using Bayesian inference on the propagation time of modeling (determined in step 704) and the obtained empirical propagation time data (obtained from process 600). Likelihood is a known function used in Bayesian inference and is well defined in the art.

[0100] Process 700 further includes a step 708 of determining an updated physical model based on the error value determined in step 706. For example, step 708 may include determining the posterior probability based on the likelihood of the updated velocity model and the prior probability density function.

[0101] The prior probability density function is defined as follows: Among them, the prior covariance matrix It is a diagonal matrix whose trace is given by the first... i The standard deviation of each model unit is composed of ,matrix N m Including the i The normalization factor of the probability density function of each unit, a vector V Composed of all model units n The earthquake parameter values ​​are composed of vectors. This includes measured seismic parameter values.

[0102] The inventors discovered that by defining the description of location x i , y i Location-related standard deviation of the prior probability density function for the parameter value The prior probability function behaves like a Gaussian function near the in-situ measurement location (applying a probability preference to specific parameter values ​​(in-situ physical parameters) and gradually moving away from them, approaching a uniform distribution).

[0103] The standard deviation of the prior probability density function is defined as follows: in It is the first point on the surface. The standard deviation, This represents the scaling factor, which determines the scaling factor as the speed of data increases. The rate at which the influence of the first point on the surface decreases with increasing distance increases until a threshold L is reached. (See above reference.) Figure 5 As described in step 510, the scaling factor is determined based on the confidence level of the influence of the x and y values ​​of the in-situ velocity measurements as the distance from the measurement point increases. Figures 8A-8C Showing standard deviation It is plotted on a surface containing a first surface point 110 and an additional reference point 112 that has not been measured. Figures 8A to 8C The standard deviations with scaling factors of 160, 50, and 25 are shown respectively. . Figures 8A to 8C Each of them has And threshold L=12.

[0104] Figure 8A The standard deviation with a scaling factor of 160 is displayed. The standard deviation increases rapidly with distance from the first point on surface 110, so that prior information only imposes constraint on velocity values ​​at and around the location measured at the first surface point 110. The additional reference point 112 is not affected by the first surface point 112 and is therefore constrained by a uniform probability density function (i.e., not constrained by in-situ velocity measurements at all).

[0105] Figure 8B The standard deviation with a scaling factor of 50 is displayed. The increase in standard deviation is relatively mild, amplifying the influence of the velocity value measured at the first surface point 110 on the modeling domain. The additional reference point 112 remains unaffected by the first surface point 110 and is therefore constrained by the uniform probability density function.

[0106] Figure 8C The standard deviation with a scaling factor of 25 is displayed. The standard deviation increases slowly, therefore, the additional reference point 112 falls within the influence range of the velocity value measured at the first surface point 110 on the modeling domain. Therefore, the probability density function at point 112 exhibits Gaussian-like behavior, imposing a probabilistic preference on the velocity value measured at the first surface point 110. (See above reference...) Figure 5 The scaling factor The value is determined based on the confidence level of the influence of the in-situ velocity value observed at the first surface point 110.

[0107] Aside from the new physical property values ​​associated with at least some units, the updated physical model is substantially identical to the initial physical model in all respects. In other words, taking into account the determined error between the empirical surface wave velocity and the modeled surface wave velocity, the updated model is an updated version of the initial physical model. This feedback process may involve least squares, Markov chain Monte Carlo, or other inversion techniques to iteratively update the initial physical model based on the updated physical model. In some examples, step 708 includes using Bayesian inference to determine the posterior probability density function of the velocity model based on the likelihood determined in step 706, the prior probability density function, and the initial model provided in step 702. Based on the initial model and the posterior probability density function, the updated physical model can then be determined based on the acceptance probability, where the acceptance probability indicates whether the initial or posterior probability density function is accepted and used as the starting point for the next round of iterations.

[0108] In such an example, the probability of acceptance is calculated by the ratio of the posterior probability of the updated velocity model to the posterior probability of the initial velocity model. A random number between 0 and 1 is generated. If the ratio is greater than the random number, the updated model is accepted; if the ratio is less than the random number, the model is returned to the initial model.

[0109] In some examples, the calculation of the acceptance probability includes the larger of the two probabilities of acceptance.

[0110] Steps 704 to 708 are typically all parts of the subroutines of the second-stage inversion process 700, and are then iterated over according to an inversion method (such as least squares inversion) or a Markov chain Monte Carlo method (such as the reverse jump Markov chain Monte Carlo method (rj-McMC method)). Figure 7 The dashed arrows in the diagram illustrate the iterative nature of process 700, indicating that the updated model determined in step 708 is used to determine the new velocity for each unit in 704. In other words, each time an updated model is determined using the error value generated by the initial model and the obtained modeling velocity for each unit, the resulting updated model is then used as the initial model for the next iteration. Iteration continues until the error value reaches the termination condition, as described above. Figure 6 The aforementioned termination criteria.

[0111] For example, when iterating steps 704 to 708 according to the rj-McMC method, the first shear wave velocity model may include splicing Voronoi elements and initiating a random number of elements while executing method 700 to determine the velocity of each element. In these examples, step 706 may also include perturbing the velocity model to form one or more perturbations of the initial velocity model, and then selecting a perturbation as the updated model. The perturbation may be one or more of the following: The velocity change of the units in the plurality of units; The position of the unit among the plurality of units changes; The birth of the unit core; The death of the unit core; The change in the first noise parameter; or The second noise parameter changes.

[0112] Then, the perturbation model can be used to determine the likelihood and posterior probability density function.

[0113] See Figure 9 The example result of the method described in this paper shows the final shear wave velocity model of the subsurface target volume. The subsurface target volume extends to a depth of 100 m in the x and y directions (z direction). The shear wave velocity values ​​are represented by shading in the figure, and the transition between regions with different shear wave velocities is visible, indicating the different composition or structure of the subsurface target region. Using the method described in this paper, the final shear wave velocity model can be determined with higher resolution and accuracy without resulting in unreasonable computation time. This subsurface model can be used to better understand the suitability of the subsurface target volume for man-made structures supporting the top or interior of the subsurface target volume.

[0114] Alternatively, additional physical properties of the subsurface target volume can be determined from the final shear wave velocity model, for example, by calculation 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.

[0115] Besides shear wave velocity, or as an alternative, other physical properties related to shear wave velocity can be used. For example, the longitudinal wave (P-wave) velocity Vp and the shear wave (or transverse wave / S-wave) velocity Vs are related to the soil elastic modulus. shear modulus and density Other physical properties are related to the following equations in linear elasticity theory: Equation 1 (P-wave velocity): Equation 2 (Shear wave velocity): Figure 10A block diagram of one implementation of a computing device 1000 is shown, in which a set of instructions can be executed to cause the computing device to perform any or more methods discussed herein. In alternative implementations, the computing device may be connected (e.g., networked) to other machines in a local area network (LAN), intranet, extranet, or the Internet. The computing device may operate as a server or client machine in a client-server network environment, or as a peer machine in a peer-to-peer (or distributed) network environment. The computing device may be a personal computer (PC), 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) specifying the action to be taken by the machine. Furthermore, although only a single computing device is shown, the term "computing device" should also be considered as a collection of any machines (e.g., computers) that individually or jointly execute a set (or more) of instructions to perform any or more methods discussed herein.

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

[0117] Processor 1002 represents one or more general-purpose processors, such as microprocessors, central processing units, etc. More specifically, 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 combinations of instruction sets. Processor 1002 may also be one or more special-purpose processors, such as application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), digital signal processors (DSPs), network processors, etc. Processor 1002 is configured to execute processing logic (instructions 1022) to perform the operations and steps discussed herein.

[0118] The computing device 1000 may also include a network interface device 1008. The 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).

[0119] It is obvious that Figure 10 Some functions of the computer device 1000 shown may be absent. For example, one or more computing devices 1000 may not require a display device 1010 (or any associated adapter). This may be the case, for example, for a particular server-side computer device 1000 that is only used for its processing power and does 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.

[0120] Data storage device 1018 may include one or more machine-readable storage media (or more specifically, one or more non-transient computer-readable storage media) 1028 on which one or more sets of instructions 1022 are stored, embodying any one or more methods or functions described herein. During execution of instructions 1022 by computer system 1000, instructions 1022 may also reside wholly or at least partially within main memory 1004 and / or processor 1002, which also constitute computer-readable storage media.

[0121] The various methods described above can be implemented by a computer program. The computer program may include computer code arranged to instruct a computer to perform one or more functions of the various methods described above. The computer program and / or code for performing these 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 temporary or non-temporary. 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 via 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 disk, random access memory (RAM), read-only memory (ROM), rigid disk, and optical disk, such as CD-ROM, CD-R / W, or DVD.

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

[0123] A "hardware component" is a tangible (e.g., non-transient) physical component (e.g., a group or one or more processors) capable of performing certain operations, which can be configured or arranged in some physical manner. A hardware component may include dedicated circuitry or logic permanently configured to perform certain operations. A hardware component may be or include dedicated processors, such as field-programmable gate arrays (FPGAs) or ASICs. A hardware component may also include programmable logic or circuitry that is temporarily configured by software to perform certain operations.

[0124] Therefore, the phrase “hardware component” should be understood to include tangible entities that can be physically constructed, permanently configured (e.g., hardwired) or temporarily configured (e.g., programmed) to operate or perform certain operations described herein.

[0125] Furthermore, modules and components can be implemented as firmware or functional circuitry within a hardware device. Additionally, modules and components can be implemented in any combination of hardware devices and software components, or solely in software (e.g., code stored or otherwise embodied in a machine-readable or transportable medium).

[0126] Unless otherwise explicitly stated, it will be apparent from the following discussion that throughout the description, the use of terms such as “receive,” “determine,” “compare,” “calculate,” “average,” “identify,” “update,” “solve,” and “output” refers to the actions and processes of a computer system or similar electronic computing device that manipulate and convert data represented as physical (electronic) quantities in computer system registers and memories into other data represented as physical quantities in computer system memory or registers or other such information storage, transmission, or display devices.

[0127] It should be understood that the above description is intended to be illustrative and not limiting. Many other implementations will become apparent to those skilled in the art upon reading and understanding the above description. While this disclosure has been described with reference to specific example implementations, it should be recognized that this disclosure is not limited to the described implementations, but can be modified and varied within the spirit and scope of the appended claims. Therefore, the specification and drawings should be considered illustrative rather than restrictive. Consequently, the scope of this disclosure should be determined by reference to the appended claims and the full scope of their equivalents.

Claims

1. A method for determining the physical properties of an underground target volume, the method comprising: Receive multiple signals detected by multiple receivers arranged on a surface above the underground target volume, wherein each of the multiple signals is detected by a corresponding receiver among the multiple receivers; Cross-correlation of the signals between the multiple receivers is performed to obtain empirical propagation time data of surface waves between multiple receiver pairs; Receive at least one underground signal detected by at least one underground receiver at a first location within an underground target volume, wherein the first location defines a first point on a surface above the underground volume; and A model for determining the physical properties of the underground target volume, the model comprising a plurality of elements, each element having a corresponding physical property value, and wherein the plurality of elements at least span a region of the surface, wherein determining the model includes: Select a subset of the plurality of receiver pairs, wherein the surface wave path between each source and the corresponding receiver passes through two or more of the first plurality of cells; Tomographic inversion is performed on the empirical propagation time data between each receiver pair in the subset to obtain the physical property values ​​for each cell. The process of performing the tomographic inversion includes obtaining a posterior probability density function corresponding to each of the first plurality of units, wherein obtaining the posterior probability density function includes: Obtain the initial velocity model; and Define a prior probability density function that indicates the velocity of each of the first plurality of units; For each receiver pair in the subset, the modeled propagation time is determined using initial velocity model values ​​associated with two or more cells traversed by the surface wave path; and The posterior probability density function of the velocity model is determined based on the propagation time of the model and the prior probability density function. The prior probability density function depends on the distance of each cell from the first point on the surface.

2. The method according to claim 1, further comprising: Receive at least one underground signal detected by at least one underground receiver at a second location within the underground target volume, wherein the second location defines a second point on the surface above the underground volume, and wherein the prior probability density function depends on the distance of each cell from the second point on the surface.

3. The method according to claim 1, wherein, Performing the tomographic inversion also includes iteratively: The initial velocity model is perturbed by one or more perturbations to form the initial velocity model; Choose the perturbation as the updated velocity model; The likelihood of the updated velocity model is determined based on the propagation time of the model and the obtained empirical propagation time data; The posterior probability is determined based on the likelihood of the updated velocity model and the prior probability density function. The acceptance probability is calculated based on the posterior probability of the updated velocity model and the posterior probability of the initial velocity model.

4. The method of claim 3, further comprising any of the following steps performed before determining the likelihood of the updated velocity model: For each receiver pair in the subset, the modeled propagation time is determined using updated velocity model values ​​associated with the two or more cells traversed by the surface wave path; or Based on previously calculated surface wave paths, updated velocity model values ​​are determined from the initial velocity model.

5. The method according to claims 3-4, wherein, The steps of perturbing the initial velocity, selecting the updated velocity model, determining the propagation time for modeling, determining the likelihood of the updated velocity, and calculating the acceptance probability are iterated until the number of iterations satisfies the first termination condition, wherein the updated velocity model of each iteration is used as the initial velocity model for the next iteration.

6. The method according to claims 4-5, wherein, Calculating the acceptance probability includes: If the posterior probability of the updated model is greater than the posterior probability of the initial model, then the updated model is accepted.

7. The method according to claims 5-6, wherein, Calculating the acceptance probability includes: The ratio of the posterior probability of the updated velocity model to the posterior probability of the initial velocity model; Generate a random number between 0 and 1; If the ratio is greater than the random number, then the updated model is accepted; If the ratio is less than the random number, then return to the initial model.

8. The method according to claims 5-7, wherein, The first termination condition is a predetermined number of iterations.

9. The method according to claims 5-8, wherein, The Monte Carlo algorithm for reversible jump Markov chains is used to perform the steps of perturbing the initial velocity, selecting the updated velocity model, determining the likelihood of the updated velocity, and calculating the acceptance probability.

10. The method according to any of the preceding claims, wherein, The disturbance can be one or more of the following: The velocity change of the units in the plurality of units; The positional changes of the units among the plurality of units; The birth of the unit core; The death of the unit core; Changes in the first noise parameter; or The second noise parameter changes.

11. The method according to any of the preceding claims, wherein, The prior probability density function is defined by the following equation: ， Among them, the prior covariance matrix It is a diagonal matrix, where the trace is formed by the first... i The standard deviation of each model unit is composed of ,matrix N m Including the i The normalization factor of the probability density function of each unit, a vector V Composed of all model units n The earthquake parameter values ​​are composed of vectors. It contains measured seismic parameter values.

12. The method according to any of the preceding claims, wherein, The standard deviation of the prior probability density function is defined as: ， in It describes the location. The standard deviation of the Gaussian probability density function at the parameter values. The first point on the surface The standard deviation at that point, This represents the scaling factor, which determines the scaling factor as the speed of data increases. Increase until the threshold L The rate at which the influence of the seismic parameter value measured at the first point on the surface decreases with distance.

13. A system comprising: - One or more processors; - One or more memories having computer-readable instructions stored thereon, the instructions being configured to cause the one or more processors to perform operations including the method described in any of the preceding claims.

14. 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 operations including the method of any one of claims 1-12.

Images