Geological surface wave three-dimensional shear wave velocity model establishment method based on small-scale rectangular coordinate system

Abstract

The invention relates to the technical field of geological survey, in particular to a three-dimensional shear wave velocity model establishment method for geological surface waves based on a small-scale rectangular coordinate system. The method comprises the steps of 1, assuming a three-dimensional transverse wave speed of geological surface waves at the current moment; 2, extracting a one-dimensional velocity model; 3, obtaining theoretical surface wave travel time based on a two-dimensional eikonal equation under a rectangular coordinate system; 4, constructing a residual error between the observation surface wave travel time and the theoretical surface wave travel time; 5, constructing a disturbance model of the two-dimensional group velocity; 6, converting the surface wave travel time residual into a surface wave travel time residual vector; 7, calculating an approximate solution of the disturbance quantity of the shear wave velocity model; 8, obtaining a three-dimensional transverse wave velocity model of the geological surface wave at the current moment; and 9, returning to the step 3 to continue iteration until inversion convergence. According to the method, errors generated in small-scale region application of an original algorithm can be reduced, and the method is more suitable for defects of a surface wave imaging boundary region of a small-scale region.

Description

Technical Field

[0001] The invention relates to the technical fields of geological survey and passive surface wave imaging, and in particular to a method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system. Background Art

[0002] Seismic surface waves propagate along the Earth's surface, and their propagation speed varies with frequency (called surface wave dispersion). Because their different frequency components have a constraining effect on the velocity structure at different depths underground, they are widely used to study shear wave velocity structures at global, regional, and engineering scales.

[0003] Generally, in the study of large-scale structures (tens of kilometers to thousands of kilometers), the method based on surface wave dispersion imaging usually adopts a two-step strategy, that is, firstly, the phase velocity diagram (or group velocity diagram) of the corresponding period under the two-dimensional grid point is constructed by the surface wave travel time of different periods; secondly, the one-dimensional (1-D) shear wave velocity structure under the grid point is inverted according to the phase velocity dispersion curve (or group velocity dispersion curve) obtained at each grid point; finally, the 1-D shear wave velocity structure under all grid points is interpolated into a three-dimensional shear wave velocity model. This strategy is also used to improve the horizontal resolution of shallow surface wave exploration. For the two-step method, it requires more periodic data of observations to obtain higher resolution phase velocity (or group velocity) images to improve the horizontal resolution of the shear wave velocity structure, while the periods with less actual observations are ignored, thus reducing the utilization of data; in addition, since the inversion of a single grid point cannot take into account the correlation of the shear wave velocity structure under adjacent grid points, if the single-point dispersion of some grid points is disturbed, false anomalies will appear.

[0004] In order to solve the above problems brought about by the two-step method, the surface wave direct imaging algorithm developed in recent years can skip the phase velocity (or group velocity) imaging and directly invert the surface wave travel time under different periods based on the ray tracing idea to obtain the underground three-dimensional shear wave velocity structure. This method is called the one-step method or direct surface wave imaging algorithm.

[0005] Surface wave dispersion imaging plays an important role in studying the velocity structure of underground media at different scales. The surface wave travel time direct imaging (DSurfTomo) algorithm developed in recent years has been widely used in the study of three-dimensional shear wave velocity structure at different scales (from several kilometers to thousands of kilometers) because it improves the utilization rate of dispersion data and considers the lateral correlation of underground inhomogeneous medium structure. Usually, for small-scale structural imaging with a detection depth of 3km (such as urban planning and mineral exploration), an independent plane rectangular coordinate system is usually used.

[0006] However, DSurfTomo is based on the fast marching method in spherical coordinates to perform ray tracing and calculate theoretical surface wave travel time, which is more suitable for regional underground structure imaging research. In order to expand DSurfTomo to smaller-scale research areas, this paper proposes a three-dimensional direct imaging algorithm based on rectangular coordinates to reduce the error caused by the surface wave direct imaging algorithm directly using spherical coordinates and improve the accuracy of the inversion results. Numerical experiments show that it can be applied to smaller array spacing (100m). Summary of the invention

[0007] In order to solve the above problems, the present invention provides a method for establishing a geological surface wave three-dimensional shear wave velocity model based on a small-scale rectangular coordinate system, comprising:

[0008] Step 1: Assume that the three-dimensional shear wave velocity of the geological surface wave at the current moment is M current ;

[0009] Step 2: Convert the three-dimensional shear wave velocity M obtained in step 1 to current Using trilinear interpolation extraction, a one-dimensional velocity Vs model is obtained;

[0010] Step 3: Substitute the results of the one-dimensional velocity Vs model obtained in step 2 into the CPS software package and calculate the group velocity U of the kth grid point based on the period of the observed dispersion curve. k (f), the group velocity U k (f) Substituting into the two-dimensional eikonal equation in the rectangular coordinate system, we obtain the theoretical surface wave travel time of the i-th surface wave ray path

[0011] Step 4: Based on the observed actual surface wave dispersion curve, obtain the observed surface wave travel time of the i-th surface wave ray path Constructing observed surface wave travel time and theoretical surface wave travel time The residual formula of δt i (f);

[0012] Step 5: Construct a perturbation model of the two-dimensional group velocity;

[0013] Step 6: Substitute the formula in step 5 into step 4 to obtain a new surface wave travel time residual, and convert it into a surface wave travel time residual vector;

[0014] Step 7: Based on the surface wave travel time residual vector in step 6, calculate the approximate solution of the shear wave velocity model disturbance

[0015] Step 8: Take the Add it to the three-dimensional shear wave velocity model assumed in step 1 to obtain the three-dimensional shear wave velocity model of the geological surface wave at the current moment;

[0016] Step 9: Return to step 3 and continue iterating until the inversion converges.

[0017] Preferably, step 2 obtains a one-dimensional velocity Vs model i=2; i = 1 represents the P-wave velocity Vp at the jth grid point in depth; i=3 represents the density of the jth grid point at depth.

[0018] Preferably, the theoretical surface wave travel time of the i-th surface wave ray path in step 3 is

[0019]

[0020] Preferably, the observed surface wave travel time in step 4 is equal to the theoretical surface wave travel time The residual formula of δt i (f) is:

[0021]

[0022] Among them, the right side of the approximate equal sign is unknown, δU k (f) is the group velocity perturbation at the kth grid point, l ik Assign the ray length of the kth grid point on the 2D group velocity diagram to the i-th surface wave ray.

[0023] Preferably, the perturbation model of the two-dimensional group velocity in step 5 is:

[0024]

[0025] Where j is the number of grid points in the depth direction, is the one-dimensional shear wave velocity disturbance.

[0026] Preferably, the formula of step 5 is substituted into step 4 to obtain:

[0027]

[0028] Among them, R α (z j ) and R ρ (z j ) is a scaling factor based on experience, G il It represents the shear wave velocity sensitivity on the ith surface wave ray path and the lth grid point in three-dimensional space. There are K grid points in the horizontal plane and J grid points in the vertical direction. The total number of grid points is K*J, m l represents the disturbance of the shear wave velocity model at the lth grid point, where i=2, for all surface wave travel time residuals, the above equations can be listed to construct the sensitivity matrix G, where Gil It is the element in the i-th row and k-th column of the G matrix;

[0029] δt i The formula of (f) is converted to: Δd = GΔm, where Δd is the residual vector, Δd and G are both known quantities, and the shear wave velocity disturbance vector Δm is obtained.

[0030] Preferably, step 7 specifically includes: calculating an approximate solution of the shear wave velocity model disturbance quantity based on the shear wave velocity disturbance vector Δm in step 6:

[0031]

[0032] In the formula, is the approximate solution of Δm, T represents the transpose, L is the smoothing operator, and λ is the damping factor.

[0033] Preferably, the three-dimensional shear wave velocity model of the geological surface wave at the current moment in step 8 is:

[0034]

[0035] Compared with the prior art, the present invention has the following beneficial effects:

[0036] This paper proposes a method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system to adapt to small-scale surface wave imaging. The effectiveness of the new method is verified by numerical experiments and actual data tests, providing a powerful tool for small-scale shallow fine surface wave exploration and small-scale body wave surface wave joint inversion strategies.

[0037] At present, the only difference between the DSurfTomo algorithm in the rectangular coordinate system proposed in this paper and the DSurfTomo algorithm in the spherical coordinate system is that they use different coordinate systems. The surface wave travel time direct imaging algorithm in the rectangular coordinate system can be directly combined with other inversion algorithms developed based on the rectangular coordinate system (such as magnetotelluric imaging or body wave travel time imaging) to reduce the non-uniqueness brought by the inversion of the physical structure of the underground medium, and avoid the problem of parameter conversion in different coordinate systems during the joint inversion process.

[0038] In addition, compared with the surface wave travel time direct imaging algorithm in the spherical coordinate system, the method proposed in the present invention can be applied to the local independent coordinate system used in engineering exploration, which is convenient for joint interpretation with exploration data such as drilling. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0040] In the attached picture:

[0041] Figure 1 It is the flow chart of the three-dimensional direct imaging algorithm of surface wave travel time in rectangular coordinate system;

[0042] Figure 2 (a) shows the array with an aperture of 196 km in different coordinate systems; Figure 2 (b) is the distribution of lateral grid points of the velocity model with an aperture of 196 km in different coordinate systems;

[0043] Figure 3 The average relative error of theoretical surface wave travel time in different coordinate systems varies with the size of the array aperture. (The top horizontal axis represents the minimum station spacing, the bottom horizontal axis represents the array aperture, and the vertical axis represents the average relative error of theoretical travel time. The black line with squares and the black line with dots represent the surface wave travel time errors in the rectangular coordinate system and the spherical coordinate system, respectively).

[0044] Figure 4 This is the distribution map of the Yangquan array in different coordinate systems. Figure 4 The left picture is the 2000 National Geodetic Coordinate System; Figure 4 The right figure is an independent plane rectangular coordinate system. (The black triangle represents the seismograph, and the four survey lines are the AA', BB', CC', and DD' profile positions respectively).

[0045] Figure 5 The distribution diagrams of station distance error and surface wave travel time error in different coordinate systems. (a) The distribution histogram of station distance absolute error; (b) The distribution histogram of station distance relative error; (c) and (d) are the distribution histograms of theoretical travel time absolute error in spherical coordinate system and rectangular coordinate system respectively; (e) and (f) are the distribution histograms of theoretical travel time relative error in spherical coordinate system and rectangular coordinate system respectively.

[0046] Figure 6 The distribution diagram of the relative error of surface wave travel time in different coordinate systems as the distance between stations changes. (a), (b), and (c) are all spherical coordinate systems; (d), (e), and (f) are all rectangular coordinate systems. (The horizontal axis represents the distance between stations, the vertical axis represents the theoretical surface wave travel time, and the colored dots represent the relative error of surface wave travel time).

[0047] Figure 7 The distribution of lateral inversion grid points in different coordinate systems. (a) Spherical coordinate system; (b) Cartesian coordinate system; (the triangles are stations and the points are grid points).

[0048] Figure 8Comparison of the detection board experiment under different coordinate systems. (a) The real Vs model in the spherical coordinate system; (b) The real Vs model in the rectangular coordinate system; (c) The Vs model restored when there is no distance error in the spherical coordinate system; (d) The Vs model restored in the rectangular coordinate system; (e) The Vs model restored when there is a distance error in the spherical coordinate system. (The horizontal slice depth is 600m, and the triangle represents the station).

[0049] Fig. 9 Comparison of surface wave direct imaging results in spherical coordinate system and rectangular coordinate system. (a) Vs slices at different depths; (b) Vs profiles along different survey lines. (AA', BB', CC', DD' survey lines are located at Figure 4 ).

[0050] Fig.10 The percentage difference distribution of the Vs model inversion results in different coordinate systems. (a) Vs percentage difference slices inverted at different depths; (b) Vs percentage difference along different survey lines. (AA', BB', CC', DD' survey lines are located at Figure 4 ). DETAILED DESCRIPTION

[0051] The following combination Figure 1-Figure 10 The preferred embodiments of the present invention are described. It should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention.

[0052] A method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system comprises:

[0053] Step 1: Assume that the three-dimensional shear wave velocity of the geological surface wave at the current moment is M current ;

[0054] Step 2: Convert the three-dimensional shear wave velocity M obtained in step 1 to current Using trilinear interpolation extraction, we get a one-dimensional velocity Vs model i=2; i = 1 represents the P-wave velocity Vp at the jth grid point in depth; i = 3 represents the density of the jth grid point at depth;

[0055] Step 3: Substitute the results of the one-dimensional velocity Vs model obtained in step 2 into the CPS software package and calculate the group velocity U of the kth grid point based on the period of the observed dispersion curve. k (f), the group velocity U k (f) Substituting into the two-dimensional eikonal equation in the rectangular coordinate system, we obtain the theoretical surface wave travel time of the i-th surface wave ray path

[0056]

[0057] Step 4: Based on the observed actual surface wave dispersion curve, obtain the observed surface wave travel time of the i-th surface wave ray path Constructing observed surface wave travel time and theoretical surface wave travel time The residual formula of δt i (f);

[0058]

[0059] Among them, the right side of the approximate equal sign is unknown, δU k (f) is the group velocity perturbation at the kth grid point, l ik Assign the ray length of the kth grid point on the two-dimensional group velocity diagram to the i-th surface wave ray;

[0060] Step 5: Construct a perturbation model of the two-dimensional group velocity:

[0061]

[0062] Where j is the number of grid points in the depth direction, is the one-dimensional shear wave velocity disturbance;

[0063] Step 6: Substitute the formula from step 5 into step 4 to obtain:

[0064]

[0065] Among them, R α (z j ) and R ρ (z j ) is a scaling factor based on experience, G il It represents the shear wave velocity sensitivity on the ith surface wave ray path and the lth grid point in three-dimensional space. There are K grid points in the horizontal plane and J grid points in the vertical direction. The total number of grid points is K*J, m l represents the disturbance of the shear wave velocity model at the lth grid point, where i=2, for all surface wave travel time residuals, the above equations can be listed to construct the sensitivity matrix G, where G il It is the element in the i-th row and l-th column of the G matrix;

[0066] δt i The formula of (f) is converted to: Δd = GΔm, where Δd is the residual vector, Δd and G are both known quantities, and the shear wave velocity disturbance vector Δm is obtained;

[0067] Step 7: Based on the shear wave velocity disturbance vector Δm in step 6, calculate the approximate solution of the shear wave velocity model disturbance

[0068]

[0069] In the formula, is the approximate solution of Δm, T represents the transpose, L is the smoothing operator, and λ is the damping factor.

[0070] Step 8: Take the Add it to the three-dimensional shear wave velocity model assumed in step 1 to obtain the three-dimensional shear wave velocity model of the geological surface wave at the current moment:

[0071]

[0072] Step 9: Return to step 3 and continue iterating until the inversion converges.

[0073] Simulation experiment:

[0074] 1. Analysis of the accuracy of numerical calculation of surface wave theoretical travel time

[0075] By constructing arrays with different aperture sizes, the accuracy of wave travel time in different coordinate systems is analyzed to verify the reliability of numerical calculation of wave travel time in the rectangular coordinate system.

[0076] like Figure 2 As shown in the figure, the array consists of 25 stations (triangles) with an aperture of about 196 km. It is a part of the Himalayan Phase II mobile array and is distributed in Baoji City and surrounding areas. https: / / pyproj4.github.io / pyproj / stable / ) The WGS-84 (EPSG:4326) coordinate system is projected onto the UTM (Universal Transverse Mercator Grid System) coordinate system (EPSG:32648, 48 represents the belt number) to obtain an independent plane rectangular coordinate system with the center of the array as the origin. Subsequently, the one-dimensional Vs velocity model was taken out from the USTClitho2.0 model (Han et al., 2022) according to the central longitude and latitude coordinates of the array, and its depth grid was interpolated to 0, 2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 30, 35, 40, 50, and 60 km, corresponding to the velocities of 2.85, 3.03, 3.20, 3.34, 3.43, 3.52, 3.55, 3.58, 3.61, 3.63, 3.69, 3.79, 3.90, 4.01, 4.28, and 4.55 km / s.

[0077] In order to ensure the consistency of model parameterization in different coordinate systems as much as possible, this paper sets the number of horizontal grid points of the model in different coordinate systems to 24×24, which corresponds to a grid size of about 7.6km×7.7km, and the number of vertical grid points is 16. This paper scales the array into arrays of different aperture sizes, and similarly follows the above steps to parameterize the model, construct arrays of different scales, and analyze the accuracy of wave travel time calculation in different coordinate systems.

[0078] For arrays of different aperture sizes, the first station is taken as the source ( Figure 2 The phase velocity with a period of 0.1 s is forward calculated by looping each grid point and is used as the velocity field for surface wave travel time calculation. Then, the surface wave travel time from the source to each station is calculated according to the rapid marching method.

[0079] In this paper, the average relative error of surface wave travel time is defined as Where N represents the number of all stations, t syn (i) represents the theoretical surface wave travel time from the source to the i-th station, which is obtained by the ratio of the distance from the source to the i-th station to the phase velocity (the two-dimensional phase velocity diagram obtained for a three-dimensional isotropic layered medium is uniform); t num (i) is the surface wave travel time from the source to the i-th station in the horizontal direction calculated by numerical calculation, which is calculated by the rapid march method.

[0080] Figure 3 The average relative error of surface wave travel time calculated by the fast marching method for arrays of different scales in different coordinate systems is shown, ranging from 0.04% to 0.12%, indicating that the theoretical travel time of surface waves calculated by the fast marching method is highly accurate under the given grids in the two coordinate systems. The experiment verifies the correctness of the travel time calculated by the fast marching method in the rectangular coordinate system. Whether in the rectangular coordinate system or the spherical coordinate system, the numerical error caused by the fast marching algorithm to solve the eikonal equation is inevitable, which is mainly related to the size of the grid division and the singularity of the source area. In the inversion process, in order to balance the calculation efficiency and calculation accuracy, a suitable grid needs to be selected.

[0081] 2. The impact of different distance calculation formulas on the calculation of surface wave theoretical travel time

[0082] In actual dispersion data extraction, the distance between two stations is a key parameter for extracting dispersion curves using the dual-station method. For small-scale arrays (<10km), the distance represents the straight-line distance between two points; for large-scale arrays, the earth is no longer a plane, and the distance between two points can be approximated as the great circle path length. The DSurfTomo package assumes that the earth is a perfect sphere and uses the Haversine formula to calculate the great circle path distance between two points. For small-scale areas, how big is the difference between the distance between two points calculated in this way and the distance calculated using the local plane coordinate system? How much impact does it have on theoretical surface wave forward modeling? To explore the above issues, this paper selects the existing technology, namely the data received by the Yangquan array deployed in the Sijiazhuang mine field area of ​​Yangquan Coal Mine (Huang Yuqi et al., 2021) for further analysis.

[0083] The existing technology uses the one-month continuous background noise data recorded by the Yangquan array, and performs seismic surface wave background noise tomography based on the DSurfTomo program package, and obtains a fine three-dimensional shear wave velocity structure at a depth of 0 to 800 meters. Finally, the relationship between seismic wave velocity and gas content is established based on rock physics experiments, and the corresponding gas content prediction is made for the work area. 2360 group velocity and 2091 phase velocity dispersion curves were obtained from the background noise cross-correlation data, with a period range of 0.2 to 1.4s and a velocity range of 1.0 to 1.5km / s. The Yangquan array has an aperture of about 6km, an average station spacing of about 0.5km, and contains 96 short-period seismometers. The actual position of the array's substations is measured by the Zhonghaida D8-RTK surveyor with an accuracy of up to millimeters. The array position information measured by RTK (Real-Time Kinematic) includes geographic coordinate information and plane coordinate information of the substations in an independent coordinate system.

[0084] Figure 4 The distribution of Yangquan station arrays in the 2000 National Geodetic Coordinate System and the independent plane rectangular coordinate system is shown.

[0085] In this paper, Haversine formula and Euler distance formula are used to calculate the great circle distance D of all station pairs. sph and the straight-line distance D car , and calculate the absolute distance error value |D of the Haversine formula sph -D car | and relative distance error |D sph -D car | / D car Distribution map( Figure 5 a and Figure 5b). The maximum absolute distance error of the platform distance calculated based on the Haversine formula is about 12m, the average is about 3.6m, and the maximum relative error of the absolute distance is 0.25%. At the same time, based on the phase velocity and group velocity data, the initial model, and the given period, this paper uses DSurfTomo in the spherical coordinate system and the rectangular coordinate system to respectively forward calculate the surface wave travel time T sph and T car , and the travel time T true for comparison. Figure 5 c and Figure 5 d shows the distribution of absolute travel time errors of surface waves in spherical coordinates and rectangular coordinates, respectively. The average values ​​of the absolute travel time errors of surface waves in spherical coordinates and rectangular coordinates are approximately 4.7 ms and 2.6 ms, respectively, and the maximum values ​​are 18.7 ms and 11.5 ms, respectively. Figure 5 e and Figure 5 f shows the distribution of the absolute relative travel time error of the wave in the spherical coordinate system and the rectangular coordinate system, where the average relative travel time error given by the rectangular coordinate system is 0.11%, and the average relative travel time error given by the spherical coordinate system is 0.19%. Furthermore, the relative travel time error of the wave in different coordinate systems varies with the distance between stations (Euler distance) as shown in Figure 6 The relative travel time errors within 0.5% in the two coordinate systems are distributed in the range of 0.5km to 6.5km between station pairs ( Figure 6 a and Figure 6 d) In the spherical coordinate system, the larger relative errors of surface wave travel time (>0.5%) are mainly concentrated in the interval between stations 0.5km to 4.2km ( Figure 6 b and Figure 6 c) In the rectangular coordinate system, the larger relative errors of surface wave travel time (>0.5%) are mainly concentrated in the range of 0.5km to 1.8km between station pairs. Figure 6 e and Figure 6 f).

[0086] 3. Inversion and comparison of actual surface wave dispersion data

[0087] In order to compare the differences in the inversion results of DSurfTomo for actual dispersion data in spherical coordinates and rectangular coordinates, this paper further uses the existing technology to extract the fundamental Rayleigh wave phase velocity and group velocity data from the Yangquan array to perform surface wave direct imaging.

[0088] In order to avoid the influence of different inversion grid divisions in different coordinate systems, this paper adopts the following strategy for horizontal grid division while ensuring that the vertical grid division is the same as the existing technology: the absolute difference of the maximum latitude (or longitude) of the station is divided equally in the spherical coordinate system; similarly, the absolute difference of the maximum X-axis (or Y-axis) coordinate of the station is also divided equally in the rectangular coordinate system. The number of equally divided grids in the north-south and east-west directions is set to 9 grids and 13 grids respectively. At the same time, in order to ensure that the ray tracing is within the calculation domain, the calculation area is expanded by 4 grids in the north-south and east-west directions ( Figure 7 ). Finally, the grid spacing in the north-south and east-west directions in the spherical coordinate system is approximately 0.005 degrees and 0.0044 degrees, respectively, and the grid spacing in the north-south and east-west directions in the rectangular coordinate system is approximately 558.2 meters and 388.3 meters, respectively.

[0089] Before the actual data inversion comparison, this paper conducts the comparison of the detection plate recovery experiment under different coordinate systems according to the actual observation system distribution. The depth range of the model is set to 0-3km. The chessboard model used in the synthetic test of the two coordinate systems is the same. They are all based on the one-dimensional average shear wave velocity model (Vs = 2km / s), and the abnormal space size is set to 4×4×12 grids in the horizontal and vertical directions, with an amplitude range of 1.7km / s to 2.3km / s. Figure 8 a and Figure 8 As shown in Fig. 2b, the three-dimensional slices of the real chessboard at a depth of 600m in the spherical coordinate system and the rectangular coordinate system are shown respectively. In this paper, the real chessboard is used to calculate the Rayleigh wave travel time consistent with the actual observation system, and 2% random noise is added to the synthetic data. The average shear wave velocity model is used as the initial model for inversion, and the recovery detection plate models in the spherical coordinate system and the rectangular coordinate system are obtained as follows: Figure 8 c and Figure 8 As shown in d, in the absence of distance error, the chessboard model results restored by the two coordinate systems are roughly equivalent. In addition, this paper replaces the synthetic travel time in the spherical coordinate system with the synthetic travel time in the rectangular coordinate system (introducing the relative error of station distance <0.25%), and performs inversion in the spherical coordinate system. The restored chessboard results are as follows Figure 8 As shown in Fig. 5, the overall result is comparable to the recovered detection board result without distance error, indicating that the smaller relative error of station distance (<0.25%) in the actual example of this paper has little effect on the inversion result.

[0090] The comparison results of the above synthetic experiments demonstrate the reliability of the wave travel time direct imaging algorithm in the rectangular coordinate system. Furthermore, this paper applies it to the actual dispersion data of Yangquan and compares it with the results given by the spherical coordinate system. In the actual inversion, the two use the same inversion grid ( Figure 7), the same initial velocity model and maximum number of iterations, where the initial velocity model comes from the existing technology and the maximum number of iterations is set to 10. Finally, the final root mean square residuals of the inversion results of the spherical coordinate system and the rectangular coordinate system are roughly the same, 0.212s and 0.211s, respectively. According to the group velocity and phase velocity depth sensitivity curves and ray distribution diagrams of Huang Yuqi et al. (2021), the results have a higher resolution at a depth of 400-800m. Fig. 9 a shows the Vs distribution comparison of the depth range obtained by the shear wave direct imaging algorithm in different coordinate systems. It can be seen that the high and low velocity anomaly distributions obtained by the two methods show a high degree of consistency. Similarly, the shear wave velocity distribution of the four survey lines in different coordinate systems ( Fig. 9 b) also showed a high degree of consistency. In addition, Fig.10 a gives Fig. 9 a The corresponding Vs slice percentage difference and Fig.10 b gives Fig. 9 b The difference in Vs percentage corresponding to different measurement lines is less than 0.2% overall.

[0091] The above shows and describes the basic principles, main features and advantages of the present invention. It should be understood by those skilled in the art that the present invention is not limited to the above embodiments, and the above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention to be protected. The scope of protection of the present invention is defined by the attached claims and their equivalents.

Claims

1. A method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system, characterized in that: include: Step 1: Assume that the three-dimensional shear wave velocity of the geological surface wave at the current moment is M current ; Step 2: Convert the three-dimensional shear wave velocity M obtained in step 1 to current Using trilinear interpolation extraction, a one-dimensional velocity Vs model is obtained; Step 3: Substitute the results of the one-dimensional velocity Vs model obtained in step 2 into the CPS software package and calculate the group velocity U of the kth grid point based on the period of the observed dispersion curve. k (f), the group velocity U k (f) Substituting into the two-dimensional eikonal equation in the rectangular coordinate system, we obtain the theoretical surface wave travel time of the i-th surface wave ray path Step 4: Based on the observed actual surface wave dispersion curve, obtain the observed surface wave travel time of the i-th surface wave ray path Constructing observed surface wave travel time and theoretical surface wave travel time The residual formula of δt i (f); Step 5: Construct a perturbation model of the two-dimensional group velocity; Step 6: Substitute the formula in step 5 into step 4 to obtain a new surface wave travel time residual, and convert it into a surface wave travel time residual vector; Step 7: Based on the surface wave travel time residual vector in step 6, calculate the approximate solution of the shear wave velocity model disturbance Step 8: Take the Add it to the three-dimensional shear wave velocity model assumed in step 1 to obtain the three-dimensional shear wave velocity model of the geological surface wave at the current moment; Step 9: Return to step 3 and continue iterating until the inversion converges.

2. The method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system according to claim 1, characterized in that: Step 2: Get the one-dimensional velocity Vs model represents the longitudinal wave velocity Vp at the jth grid point in depth; Represents the density of the jth grid point at depth.

3. The method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system according to claim 2, characterized in that: Theoretical surface wave travel time of the i-th surface wave ray path in step 3 4. The method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system according to claim 3, characterized in that: Observed surface wave travel time and theoretical surface wave travel time in step 4 The residual formula of δt i (f) is: Among them, the right side of the approximate equal sign is unknown, δU k (f) is the group velocity perturbation of the kth grid point, l ik Assign the ray length of the kth grid point on the 2D group velocity diagram to the i-th surface wave ray.

5. The method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system according to claim 4, characterized in that: The perturbation model of the two-dimensional group velocity in step 5 is: Where j is the number of grid points in the depth direction, is the one-dimensional shear wave velocity disturbance.

6. The method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system according to claim 5, characterized in that: Substituting the formula from step 5 into step 4 yields: Among them, R α (z j ) and R ρ (z j ) is a scaling factor based on experience, G il It represents the shear wave velocity sensitivity on the ith surface wave ray path and the lth grid point in three-dimensional space. There are K grid points in the horizontal plane and J grid points in the vertical direction. The total number of grid points is K*J, m l represents the disturbance of the shear wave velocity model at the lth grid point, where For all surface wave travel time residuals, the above equations can be listed to construct the sensitivity matrix G, where G il It is the element in the i-th row and l-th column of the G matrix; δt i The formula of (f) is converted to: Δd = GΔm, where Δd is the residual vector, Δd and G are both known quantities, and the shear wave velocity disturbance vector Δm is obtained.

7. The method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system according to claim 6, characterized in that: Step 7 is specifically as follows: based on the shear wave velocity disturbance vector Δm in step 6, an approximate solution of the shear wave velocity model disturbance is calculated: In the formula, is the approximate solution of Δm, T represents the transpose, L is the smoothing operator, and λ is the damping factor.

8. The method for establishing a three-dimensional shear wave velocity model of geological surface waves based on a small-scale rectangular coordinate system according to claim 7, characterized in that: The three-dimensional shear wave velocity model of the geological surface wave at the current moment in step 8 is:

Images