A surface wave tomography method and system considering terrain effect

Abstract

The present invention relates to a surface wave tomography method and system that considers terrain effects, wherein the method comprises: obtaining seismic surface wave dispersion data to be inverted; using the shortest path method to process the seismic surface wave dispersion data to be inverted to obtain the propagation trajectory of seismic waves underground, and updating the node travel times on all trajectories; using the velocity-density empirical relationship in rock physics experiments to establish a functional relationship between the shear wave velocity of the surface wave direct inversion method and the longitudinal wave velocity and medium density; modeling the node travel times based on the functional relationship and solving to obtain the inversion result. The present invention solves the problem of large deviation in the inversion of undulating terrain structures by existing surface wave imaging methods, and the inversion model makes the inversion result of the geological structure more accurate and reliable compared to the model that does not consider the terrain effect.

Description

Technical Field

[0001] The present invention relates to the technical field of surface wave tomography, and in particular to a surface wave tomography method and system taking terrain effect into consideration. Background Art

[0002] The changes in subsurface elastic parameters are crucial to understanding the dynamic processes in the shallow and deep layers of the Earth. With the development of short-period dense array technology, the increasing amount of surface wave data has made surface wave tomography a powerful tool for obtaining these elastic parameters. The dispersion effect of surface waves is one of the most important properties in the surface wave propagation process, which makes it possible to simulate the surface wave travel time based on a ray tracing problem in two-dimensional space. Commonly used surface wave travel time tomography methods include large-circle surface wave tomography, eikonal equation tomography, and finite frequency tomography. However, in areas with large terrain undulations, most of these traditional methods ignore the impact of terrain on travel time, mainly because scholars have assumed that the scale of tomography can be ignored relative to the wavelength and propagation depth of the surface wave. For example, scholars have shown through spectral element method simulation that on a profile with a length of 200 kilometers, a 2-kilometer terrain undulation only underestimates the phase velocity of a 3-second period by 0.7%.

[0003] However, with the application of short-period dense seismic arrays and distributed fiber optic sensors, the wavelength of short-period surface waves and the scale of imaging are comparable. In areas with drastic topographic changes, the influence of topography cannot usually be ignored, such as detectors deployed on the Qinghai-Tibet Plateau and seafloor seismometers in coastal areas. Therefore, in this case, the impact of topography on forward simulation or data processing must be considered. Some scholars have proposed using an improved quasi-bending method to develop a surface wave ray tracing system that takes into account topographic effects, but this method has certain uncertainties, especially in areas with large topographic changes. Summary of the invention

[0004] To solve the above problems, an object of the embodiments of the present invention is to provide a surface wave tomography method and system that takes terrain effects into consideration.

[0005] A surface wave tomography method considering terrain effect, comprising:

[0006] Obtaining seismic surface wave dispersion data to be inverted;

[0007] The seismic surface wave dispersion data to be inverted are processed by using the shortest path method to obtain the propagation trajectory of the seismic wave in the underground, and the travel time of all nodes on the trajectory is updated;

[0008] By using the empirical velocity-density relationship in rock physics experiments, a functional relationship is established between the shear wave velocity of the surface wave direct inversion method, the longitudinal wave velocity and the medium density.

[0009] The node travel time is modeled based on the functional relationship and solved to obtain an inversion result.

[0010] Preferably, the method of using the velocity-density empirical relationship in rock physics experiments to establish a functional relationship between the shear wave velocity, the longitudinal wave velocity and the medium density in the surface wave direct inversion method includes:

[0011] The Brochet formula is used to determine the functional relationship between the shear wave velocity, the longitudinal wave velocity and the medium density:

[0012] α=F α (β),ρ=F ρ (β)

[0013] Among them, α represents the longitudinal wave velocity, β represents the shear wave velocity, ρ represents the medium density, and F a (β) represents the functional mapping relationship between shear wave velocity and longitudinal wave velocity, F ρ (β) represents the functional mapping relationship between shear wave velocity and medium density.

[0014] Preferably, modeling the node travel time based on the functional relationship and solving to obtain an inversion result includes:

[0015] Modeling the node travel time based on the functional relationship to obtain a three-dimensional Frechet derivative formula;

[0016] The three-dimensional Frechet derivative formula is converted into a matrix form and a regularization term is added, and the inversion result is obtained by solving it through the least squares conjugate gradient method.

[0017] Preferably, the three-dimensional Fréchet derivative formula is:

[0018]

[0019] Where ω represents the surface wave frequency, i, j, and k represent the number of time nodes, shear wave velocity nodes, group velocity or phase velocity nodes, δ represents the disturbance, and d / dβ j is the total partial derivative of the shear wave velocity, c k (ω) represents the group velocity or phase velocity of the corresponding frequency, t i (ω) represents the travel time of the node.

[0020] The present invention also provides a surface wave tomography system taking terrain effect into consideration, comprising:

[0021] A data acquisition module, used to acquire seismic surface wave dispersion data to be inverted;

[0022] A data processing module, used to process the seismic surface wave dispersion data to be inverted using the shortest path method to obtain the propagation trajectory of the seismic wave in the underground, and update the travel time of all nodes on the trajectory;

[0023] Functional relationship building module, used to establish functional relationship between shear wave velocity, longitudinal wave velocity and medium density in surface wave direct inversion method by using velocity-density empirical relationship in rock physics experiment;

[0024] The inversion module is used to model the node travel time based on the functional relationship and solve it to obtain an inversion result.

[0025] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps in the above-mentioned surface wave tomography method considering terrain effects.

[0026] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0027] The present invention relates to a surface wave tomography method and system taking terrain effect into consideration. Compared with the prior art, the present invention solves the problem of large deviation in inversion of undulating terrain structure by the prior surface wave imaging method, and the inversion model makes the inversion result of geological structure more accurate and reliable than the model not taking terrain effect into consideration.

[0028] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, preferred embodiments are given below and described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0030] Figure 1 The flowchart of a surface wave tomography method considering terrain effect according to an embodiment of the present invention.

[0031] Figure 2 1 is a travel time field with a Gaussian mountain one-dimensional model according to an embodiment of the present invention; (a) represents the velocity profile of the one-dimensional model; (b) represents the travel time field, the solid and dashed lines are contour lines with and without terrain, and the black stars and triangles represent the source and station, respectively.

[0032] Figure 33D ray paths of a Gaussian mountain model according to an embodiment of the present invention; (a) ray path without terrain influence, (b) ray path including terrain influence.

[0033] Figure 4 It is a schematic diagram of the three-dimensional shear wave velocity structure inversion principle according to an embodiment of the present invention; (a) represents the terrain and station distribution of the detection plate model, the red triangles are the stations used, and the dotted line marked as A-A' is the profile used to display the inversion results; (b) represents a histogram of the dispersion data of each period used.

[0034] Figure 5 It is a horizontal slice with a depth of 1.2 km according to an embodiment of the present invention. (a) Inversion result without considering the terrain; (b) Inversion result after considering the terrain; (c) True model. The white dotted line indicates the terrain position.

[0035] Figure 6 The inversion results of the AA' profile through the Gauss Mountain according to an embodiment of the present invention are as follows: (a) Inversion result without considering the terrain; (b) Inversion result after considering the terrain; (c) True model. DETAILED DESCRIPTION

[0036] In the description of the present invention, it should be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside", "clockwise", "counterclockwise" and the like indicate orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the referred device or element must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be understood as limiting the present invention.

[0037] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, the meaning of "plurality" is two or more, unless otherwise clearly and specifically defined.

[0038] In the present invention, unless otherwise clearly specified and limited, the terms "installed", "connected", "connected", "fixed" and the like should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or it can be an indirect connection through an intermediate medium, or it can be the internal communication of two components. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0039] See also Figure 1 , a surface wave tomography method considering terrain effect, comprising:

[0040] Step 1: Obtain the seismic surface wave dispersion data to be inverted; In ray-based surface wave tomography, a local method is usually used to determine the dispersion curve (data) of each point on the free surface. This method regards the three-dimensional model as a series of one-dimensional uniform elastic isotropic layers. Then, the Thomson-Haskell matrix method is used to synthesize the dispersion curves of each one-dimensional model to obtain the two-dimensional dispersion pattern on the free surface. The perturbation of the dispersion curve can be converted into elastic parameter perturbation solved by the variational method.

[0041] Step 2: Process the seismic surface wave dispersion data to be inverted using the shortest path method to obtain the propagation trajectory (ray trajectory) of the seismic wave underground, and update the travel time of all nodes on the trajectory;

[0042] (1) In the SPM method, a two-dimensional (curved) free-form surface is discretized using quadrilateral elements. This type of discretization allows the free expression of topographic changes above. The vertex of each element is called the primary node. (2) In order to accurately determine the position of each ray segment, auxiliary nodes are defined on one side of each element. The number of auxiliary nodes is N. x and N y It indicates that the node density can be adjusted according to the specific application, and the speed of the auxiliary node is obtained by bilinear interpolation from the speed of the main node. The combination of main nodes and auxiliary nodes can increase the ray coverage of each unit and reduce the total number of units. It is worth noting that although the present invention uses a regular quadrilateral unit grid, other options such as triangular or irregular quadrilateral grids can also be a combination of main-auxiliary nodes. (3) The nodes in the discretization, including the main nodes and auxiliary nodes, are respectively represented by the global numbering system {x n}, where the subscript n is an integer between 1 and n (indicating the total number of nodes), and the continuity of travel time along adjacent unit grids takes into account the frequency-dependent travel time of each node, represented by {t n(ω)}; for a particular event ray trace, the event location is assigned zero travel time, while all nodes are assigned infinite travel time. In addition, all unit grids are assigned with the label "far". Subsequently, the travel times of all nodes can be calculated on the source unit grid Ω S Calculated in.

[0043] Step 3: Using the velocity-density empirical relationship in rock physics experiments, establish a functional relationship between the shear wave velocity of the surface wave direct inversion method, the longitudinal wave velocity and the medium density;

[0044] Since the dispersion curve is more sensitive to the change of shear wave velocity, the change of shear wave is usually inverted separately, and then the Brochet formula is used to determine the functional relationship between the P-wave velocity and density and the S-wave velocity:

[0045] α=F α (β),ρ=F ρ (β)(1)

[0046] Among them, α represents the longitudinal wave velocity, β represents the shear wave velocity, ρ represents the medium density, ω represents the surface wave frequency, i, j, k represent the number of time nodes, the number of shear wave velocity nodes, the number of group velocity or phase velocity nodes, respectively, and δ represents the disturbance.

[0047] Step 4: Model the node travel time based on the functional relationship and solve it to obtain an inversion result.

[0048] Furthermore, the step 4 comprises:

[0049] The node travel time is modeled based on the functional relationship to obtain a three-dimensional Fréchet derivative formula; wherein the three-dimensional Fréchet derivative formula is:

[0050]

[0051] Where ω represents the surface wave frequency, i, j, and k represent the number of time nodes, shear wave velocity nodes, group velocity or phase velocity nodes, δ represents the disturbance, and d / dβ j is the total partial derivative of the shear wave velocity, c k (ω) represents the group velocity or phase velocity of the corresponding frequency, t i (ω) represents the travel time of the node.

[0052] The three-dimensional Frechet derivative formula is converted into a matrix form and a regularization term is added, and the inversion result is obtained by solving the least squares conjugate gradient method;

[0053] In practical applications, formula (2) can be reorganized into a matrix form:

[0054] Gδβ=δt(3)

[0055] Formula (3) can be further solved by adding regularization terms:

[0056]

[0057] Where G represents the operator that calculates the observed data from the model parameters, t represents the travel time, β represents the shear wave velocity, L and I are the second-order and zero-order Tikhonov regularization terms, γ and θ are the corresponding coefficients, and δ represents the perturbation. This system of linear equations can be effectively solved by the least squares conjugate gradient (LSQR) method. The current model (the model at the nth iteration) is the iterative update formula:

[0058] β (n+1) =β (n) +δβ (n) (5)

[0059] In one embodiment, the goal of the present invention is to simulate a simple one-dimensional travel time field model to observe the terrain effect. The model used in this test is a rectangular block with horizontal dimensions of 10×10 kilometers and a depth of 20 kilometers. The terrain of this area is a Gaussian mountain with a maximum central height of 0.5 kilometers and a standard deviation of 0.5 kilometers (centered in the middle of the area). Figure 2 (b)). The velocity model below this region is a simple one-dimensional isotropic layered model ( Figure 2 (a)). The elastic parameters of the media above and below the mountain are the same. Next, the phase velocity travel-time field with a period of 1 second is calculated and synthesized using SPM.

[0060] Walking time Figure 2 As shown, when the seismic wave does not interact with the Gaussian Mountain, the travel time corresponding to the theoretical travel time (dashed line) is not affected by the terrain (solid line). However, when the seismic wave interacts with the Gaussian Mountain, the travel time is delayed and deviates from the concentric circles. Figure 2 and Figure 3 This is illustrated by the ray paths in Figure 2. The rays orbit around the obstacle and the travel-time field and ray paths show similar characteristics to those in the low-speed region.

[0061] This embodiment can qualitatively analyze the deviation caused by the terrain that may appear in the inversion. It can be observed that ignoring the influence of the terrain in the initial model may lead to a time delay between the observed data and the synthetic data. In areas with changing terrain, when inverting the phase / group velocity pattern, the time delay will cause low-velocity anomalies in two-dimensional space. Therefore, this type of anomaly caused by considering the terrain is a pseudo low-velocity anomaly and is superimposed on the real anomaly of each stage. Therefore, if imaging is performed in a medium with undulating terrain, the inverted velocity model that does not consider the terrain will be biased and difficult to interpret.

[0062] The above method is applied to a numerical model example to perform a 3D shear wave velocity structure inversion. The study area has a horizontal dimension of 10 × 10 km, including a 1 km high Gaussian mountain with a maximum width of 2 km. Figure 4 (a) shows the distribution of stations and topography in the region. The Rayleigh wave phase / group velocity data used for inversion is derived from the cross-correlation of 87 stations, with a time span from 0.2 seconds to 7.4 seconds. The histogram of the dispersion data is shown in Figure 4 (b) as shown.

[0063] The check plate model is constructed by adding high and low anomaly disturbances in a 2×2 km grid to the previous 1D model, to a depth of 15 km. The velocity in the Gaussian mountain area (above 0 km) remains the same as at 0 km. To illustrate the influence of the topography, the mountain is located directly above a high velocity anomaly. The SPM method is then used to obtain observational data in the form of synthetic surface wave travel times. Finally, the check plate model is obtained from the initial 1D model using a 3D direct inversion method.

[0064] Figure 5 and Figure 6 The horizontal slice model (depth 1.2 km) and vertical section (AA') model are obtained by 3D direct inversion. It is worth noting that by adding the influence of terrain (such as Figure 5 (b) and Figure 6 (b) shows that the velocity structure in both horizontal and vertical sections can be successfully recovered. In addition, it can be seen from the results that if the influence of terrain is not considered, such as Figure 5 As shown in (a), the high-speed anomaly located directly below Mount Gauss is inverted to a low-speed anomaly. Figure 6 Similar features are also observed in (a). This phenomenon is consistent with the previous discussion. If the delayed travel time due to topography is not considered in the inversion imaging, it will produce low-velocity pseudo anomalies in the velocity model.

[0065] The method of the present invention is based on a simple formula and is easy to implement. In addition, the use of a heap sort algorithm in the update step allows the efficient identification of unit grids with the minimum travel time, with a time complexity of O(subsequence). Secondly, SPM is flexible in the choice of grid structure and does not rely on a specific coordinate system. Therefore, the range of grid configurations (triangular or quadrilateral grids) can be used to adapt to different modeling scenarios. Thirdly, the accuracy of travel time calculation can be improved by increasing the density of auxiliary nodes along the sides of each unit grid. In addition, once the propagation time is updated, the ray path between the source and any given node can be automatically determined. Finally, SPM allows "delayed calculation", in which case, once all travel times of the required stations are updated, it can be stopped without calculating the entire travel time. The present invention solves the problem of large deviation in the inversion of undulating terrain structures by the existing surface wave imaging method through the above-mentioned scheme, and the inverted model makes the inversion result of the geological structure more accurate and reliable relative to the model that does not consider the terrain effect.

[0066] The present invention also provides a surface wave tomography system taking terrain effect into consideration, comprising:

[0067] A data acquisition module, used to acquire seismic surface wave dispersion data to be inverted;

[0068] A data processing module, used to process the seismic surface wave dispersion data to be inverted using the shortest path method to obtain the propagation trajectory of the seismic wave in the underground, and update the travel time of all nodes on the trajectory;

[0069] Functional relationship building module, used to establish functional relationship between shear wave velocity, longitudinal wave velocity and medium density in surface wave direct inversion method by using velocity-density empirical relationship in rock physics experiment;

[0070] The inversion module is used to model the node travel time based on the functional relationship and solve it to obtain an inversion result.

[0071] An embodiment of the present invention further provides a computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a processor, the various processes of the above-mentioned surface wave tomography method embodiment considering terrain effects are implemented, and the same technical effect can be achieved. To avoid repetition, it will not be described here.

[0072] The above is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any technical solution that can be easily thought of by a person skilled in the art within the technical scope disclosed by the present invention should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the protection scope of the claims.

Claims

1. A surface wave tomography method considering terrain effect, characterized in that: include: Obtaining seismic surface wave dispersion data to be inverted; The seismic surface wave dispersion data to be inverted are processed by using the shortest path method to obtain the propagation trajectory of the seismic wave in the underground, and the travel time of all nodes on the trajectory is updated; In the SPM method, the two-dimensional free-form surface is discretized using quadrilateral elements. This type of discretization allows the free expression of the terrain changes above. The vertex of each element is called the main node. In order to accurately determine the position of each ray segment, auxiliary nodes are defined on one side of each element. The number of auxiliary nodes is N. x and N y Indicates that the speed of the auxiliary node is obtained by bilinear interpolation from the speed of the main node; the combination of the main node and the auxiliary node can increase the ray coverage of each unit and reduce the total number of units; the nodes in the discretization, including the main node and the auxiliary node, are represented by the global numbering system {x n }, where the subscript n is an integer between 1 and n, indicating the total number of nodes. The continuity of travel time along adjacent unit grids takes into account the frequency-dependent travel time of each node, and is represented by {t n (ω)} represents; for the ray tracing of a particular event, the location of the event is assigned zero travel time, while all nodes are assigned infinite travel time; In addition, all unit grids are specified with the tag "far". Then, the travel times of all nodes can be calculated on the source unit grid Ω. S Medium calculation; By using the empirical velocity-density relationship in rock physics experiments, a functional relationship is established between the shear wave velocity of the surface wave direct inversion method, the longitudinal wave velocity and the medium density. The node travel time is modeled based on the functional relationship and solved to obtain an inversion result.

2. A surface wave tomography method considering terrain effect according to claim 1, characterized in that: The method uses the velocity-density empirical relationship in rock physics experiments to establish a functional relationship between the shear wave velocity of the surface wave direct inversion method, the longitudinal wave velocity and the medium density, including: The Brochet formula is used to determine the functional relationship between the shear wave velocity, the longitudinal wave velocity and the medium density: α=F α (b),ρ=F ρ (b) Among them, α represents the longitudinal wave velocity, β represents the shear wave velocity, ρ represents the medium density, and F α (β) represents the functional mapping relationship between shear wave velocity and longitudinal wave velocity, F ρ (β) represents the functional mapping relationship between shear wave velocity and medium density.

3. A surface wave tomography method considering terrain effect according to claim 2, characterized in that: The step of modeling the node travel time based on the functional relationship and solving the model to obtain an inversion result includes: Modeling the node travel time based on the functional relationship to obtain a three-dimensional Frechet derivative formula; The three-dimensional Fréchet derivative formula is: Where ω represents the surface wave frequency, i, j, and k represent the number of time nodes, shear wave velocity nodes, group velocity or phase velocity nodes, δ represents the disturbance, and d / dβ j is the total partial derivative of the shear wave velocity, c k (ω) represents the group velocity or phase velocity of the corresponding frequency, t i (ω) represents the travel time of the node; The three-dimensional Frechet derivative formula is converted into a matrix form and a regularization term is added, and the inversion result is obtained by solving it through the least squares conjugate gradient method.

4. A surface wave tomography system considering terrain effects, characterized in that: include: A data acquisition module, used to acquire seismic surface wave dispersion data to be inverted; A data processing module, used to process the seismic surface wave dispersion data to be inverted using the shortest path method to obtain the propagation trajectory of the seismic wave in the underground, and update the travel time of all nodes on the trajectory; In the SPM method, the two-dimensional free-form surface is discretized using quadrilateral elements. This type of discretization allows the free expression of the terrain changes above. The vertex of each element is called the main node. In order to accurately determine the position of each ray segment, auxiliary nodes are defined on one side of each element. The number of auxiliary nodes is N. x and N y Indicates that the speed of the auxiliary node is obtained by bilinear interpolation from the speed of the main node; the combination of the main node and the auxiliary node can increase the ray coverage of each unit and reduce the total number of units; the nodes in the discretization, including the main node and the auxiliary node, are represented by the global numbering system {x n }, where the subscript n is an integer between 1 and n, indicating the total number of nodes. The continuity of travel time along adjacent unit grids takes into account the frequency-dependent travel time of each node, and is represented by {t n (ω)} represents; for the ray tracing of a particular event, the location of the event is assigned zero travel time, while all nodes are assigned infinite travel time; In addition, all unit grids are specified with the tag "far". Then, the travel times of all nodes can be calculated on the source unit grid Ω. S Medium calculation; Functional relationship building module, used to establish functional relationship between shear wave velocity, longitudinal wave velocity and medium density in surface wave direct inversion method by using velocity-density empirical relationship in rock physics experiment; The inversion module is used to model the node travel time based on the functional relationship and solve it to obtain an inversion result.

5. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps in a surface wave tomography method considering terrain effects as claimed in any one of claims 1 to 3 are implemented.

Images