Shallow velocity model construction method based on joint inversion

Abstract

The invention discloses a shallow velocity model construction method based on joint inversion, and relates to the technical field of geophysical exploration. According to the method, high-frequency Rayleigh wave information and far-earthquake cross convolution waveforms in natural earthquakes are introduced into joint inversion, the influence of parameters on a speed structure model is considered, interference of other types of waves on Rayleigh waves is avoided, noise segments dominant by Rayleigh wave energy are effectively screened out, and the method is suitable for large-scale popularization and application. The difficulty of extracting the ellipticity of the high-frequency Rayleigh wave from a single waveform is effectively solved; the three-dimensional velocity structure model of the shallow imgabs0 # and the three-dimensional velocity structure model of the shallow imgabs1 # of the whole earth crust can be effectively simulated through joint inversion, especially for a sedimentary layer, research of deep structure scientific problems of a research area and exploration and development of oil and gas resources are facilitated, the seismic oscillation simulation precision can be improved, and earthquake resistance and disaster reduction work in an earthquake area can be assisted.

Description

Technical Field

[0001] The invention relates to the technical field of geophysical exploration, and in particular to a method for constructing a shallow velocity model based on joint inversion. Background Art

[0002] Post-earthquake surveys show that basins have a significant amplification effect on ground motion. The physical property difference between low-speed, loose sedimentary layers and underlying bedrock is generally considered to be the direct cause of increased earthquake damage. It significantly prolongs the duration of ground motion in the basin and amplifies the amplitude of ground motion. Further studies have also shown that the thickness of the sedimentary layer and the shallow velocity (P wave and S wave) will significantly affect strong ground motion. Among them, the thickness of the sedimentary layer mainly controls the low-frequency resonance behavior of ground motion, while the shallow velocity mainly affects the ground motion of 1Hz or even higher frequency. In the study of engineering site effects, the average shear wave velocity at a depth of 30 meters near the surface is mainly considered. However, in order to predict the seismic wave propagation effect caused by near-surface structure, especially to estimate the relatively lower frequency (period>1s) ground motion in the basin, it is necessary to accurately determine the three-dimensional P wave and S wave velocity structure of the sedimentary layer above the bedrock. Previous researchers have conducted detection and research on basin structure based on active and passive source seismic detection and gravity methods, and achieved excellent results. However, affected by the distribution of geophysical detection points, the existing regional velocity models cannot meet the needs of earthquake site effect assessment and earthquake disaster prevention and reduction. Especially at the regional scale, the shear wave velocity and the corresponding longitudinal and shear wave velocity ratio in the shallow layer (0-10km) have not been well resolved. Summary of the invention

[0003] In view of the above-mentioned deficiencies in the prior art, the present invention provides a method for constructing a shallow velocity model based on joint inversion, which solves the problem that the prior art cannot well solve the shallow shear wave velocity and the corresponding longitudinal and shear wave velocity ratio on a regional scale.

[0004] In order to achieve the above-mentioned object of the invention, the technical solution adopted by the present invention is: A method for constructing a shallow velocity model based on joint inversion is provided, which comprises the following steps: S1, obtaining seismic data from fixed stations and mobile stations in the area to be measured and preprocessing them to obtain preprocessed seismic data; S2, constructing the background noise cross-correlation function based on the preprocessed seismic data, and using the Eikonal imaging method to obtain the corresponding high-frequency Rayleigh wave phase velocity distribution; S3, using the FDFA algorithm and the HVSR method to process the preprocessed seismic data to obtain the corresponding Rayleigh wave ellipticity; S4, performing data processing and screening on the pre-processed seismic data to obtain a corresponding teleseismic cross-convolution waveform; S5. Based on the high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity and teleseismic cross-convolution waveform, a shallow velocity structure model is constructed using joint inversion.

[0005] Furthermore, the preprocessed seismic data in step S1 include Z-component continuous waveform data, three-component continuous waveform data and teleseismic event data; wherein, the preprocessing of the Z-component continuous waveform data and the three-component continuous waveform data adopts resampling, removal of seismic events and removal of mean, and the preprocessing of the teleseismic event data adopts resampling and removal of mean; the teleseismic event data are waveform data of seismic events with a magnitude of not less than 5.5 and a distance from the epicenter of 30° to 90°.

[0006] Further, step S2 includes the following steps: S2-1. Cut the 10 Hz continuous noise data into 1-hour segments every day, perform cross-correlation calculation and superposition on the Z component continuous waveform data of each station pair every hour, and obtain the ZZ component background noise cross-correlation function between each station pair on the same day, i.e., NFC (Noise Cross-Correlation Function); S2-2, perform superposition processing on each NFC to obtain the ZZ component NCF waveform data between each station pair during the entire observation time; S2-3, picking up the Rayleigh wave dispersion information between the station pairs corresponding to the ZZ component NCF waveform data; S2-4. Use the Eikonal imaging method to process the Rayleigh wave dispersion information and obtain the high-frequency Rayleigh wave phase velocity distribution.

[0007] Further, step S3 includes the following steps: S3-1, performing time window segmentation processing on the three-component continuous waveform data to obtain segmented three-component continuous waveform data; S3-2, using the FDFA algorithm to calculate the three-component continuous waveform data after segmentation in each time window, and obtain the phase difference between the horizontal component and the vertical component in each time window; S3-3, based on the phase difference between the horizontal component and the vertical component in each time window, the time windows are screened, and the time windows with phase differences in the range of 75° to 115° are retained to obtain screened time windows; and the screened time windows are used as the band dominated by Rayleigh waves; S3-4. Use the HVSR method to calculate the band dominated by Rayleigh waves and obtain the corresponding Rayleigh wave ellipticity.

[0008] Furthermore, the calculation formula of the Rayleigh wave ellipticity in step S3-4 is:

[0009] in, represents the frequency corresponding to the band dominated by Rayleigh waves, represents the Rayleigh wave ellipticity, , Represents two orthogonal horizontal components at frequencies The spectrum of Indicates the vertical component at frequency spectrum.

[0010] Further, step S4 includes the following steps: S4-1, performing time window processing on the teleseismic event data to obtain processed teleseismic event data; S4-2, screening and superimposing the processed teleseismic event data to obtain superimposed teleseismic event data; S4-3. Perform a cross-convolution waveform operation on the superimposed teleseismic event data to obtain a corresponding teleseismic cross-convolution waveform, and construct a corresponding cross-convolution error function.

[0011] Furthermore, the specific process of step S4-2 is: The ENZ three-component data of the processed teleseismic event data are rotated using the rotation matrix to obtain the RTZ component, and the P-wave R component data and Z component data are extracted; The P-wave R component data and Z component data with waveform similarity not less than the threshold and signal-to-noise ratio greater than 5 are superimposed to obtain the body wave stable waveform of each station, that is, the superimposed teleseismic event data.

[0012] Furthermore, the calculation formula of the teleseismic cross-convolution waveform in step S4-3 is:

[0013]

[0014]

[0015] in, represents the teleseismic cross-convolution waveform, , represent the vertical and radial components of the teleseismic cross-convolution waveform, respectively. Indicates the time, , Respectively indicate time No. The theoretical vertical waveform and theoretical radial waveform of the earth model, represents the total number of earthquake events, represents the summation function, , denote the vertical and radial Green's functions of the Earth, respectively, , Respectively represent The instrument response and source time function of each earthquake event, , Respectively represent Vertical and radial component waveforms of a seismic event.

[0016] Further, step S5 includes the following steps: S5-1. Determine the weight parameters of high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity and teleseismic cross-convolution waveform according to the research objectives; S5-2, constructing a joint inversion error function based on the weight parameters of step S5-1; S5-3. Obtain the shear wave velocity of the sedimentary layer and upper crust in the area to be measured and longitudinal wave velocity ; S5-3, based on the joint inversion error function, the thickness of the sedimentary layer, and The high-frequency Rayleigh wave phase velocity, high-frequency Rayleigh wave ellipticity and teleseismic cross-convolution waveform are jointly inverted using the Bayesian Monte Carlo method as inversion parameters to obtain the one-dimensional and Speed ​​structure; S5-4, constructing a joint inversion error function based on the cross convolution error function; S5-5. One-dimensional interpolation of each station by grid interpolation within the region and The velocity structure is processed to obtain the shallow three-dimensional and Velocity structure model.

[0017] The beneficial effects of the present invention are as follows: by introducing teleseismic cross-convolution waveforms in the joint inversion and taking into account the influence of parameters on the velocity structure model, the interference of other types of waves on Rayleigh waves is avoided, and the noise fragments dominated by Rayleigh wave energy are effectively screened out, which effectively solves the difficulty of extracting Rayleigh wave ellipticity from a single waveform; the shallow three-dimensional of the entire crust can be effectively simulated through time inversion. and The velocity structure model, especially for sedimentary layers, is beneficial to the research of deep structural scientific issues and the exploration and development of oil and gas resources in the study area. It can also improve the accuracy of seismic motion simulation and assist in earthquake-resistant disaster reduction work in earthquake-prone areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 is a method flow chart; Figure 2The data waveform diagram and curve diagram of the present invention; wherein, (a) is a 24-hour three-component continuous waveform diagram; (b) is a phase difference curve diagram of the horizontal and vertical components; (c) is a comparison curve diagram of the Rayleigh wave ellipticity values ​​before and after screening. DETAILED DESCRIPTION

[0019] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0020] like Figure 1 As shown in FIG. 1 , a method for constructing a shallow velocity model based on joint inversion includes the following steps: S1. Obtain the seismic data from the fixed and mobile stations in the area to be tested and preprocess them to obtain the preprocessed seismic data. The area to be tested is North China, where strong earthquakes are very active inland. This study uses seismic data recorded by 156 fixed stations of the seismic network and 466 ChinArrayIII broadband mobile stations in a certain area. The mobile seismic station records for about 2 years, which can provide sufficient continuous waveforms and seismic data.

[0021] We plan to collect and organize drilling data and artificial seismic exploration results in the study area, and select data related to sediment thickness, shear wave velocity, and the ratio of longitudinal and transverse wave speeds The data are compared with the one-dimensional models of nearby stations to verify the reliability of the model.

[0022] The preprocessed seismic data in step S1 include Z-component continuous waveform data, three-component continuous waveform data and teleseismic event data; wherein, the preprocessing of the Z-component continuous waveform data and the three-component continuous waveform data adopts resampling, removal of seismic events and removal of mean, and the preprocessing of the teleseismic event data adopts resampling and removal of mean; the teleseismic event data are waveform data of seismic events with a magnitude of not less than 5.5 and a distance from the epicenter of 30° to 90°.

[0023] S2. Construct a background noise cross-correlation function based on the preprocessed seismic data, and use the Eikonal imaging method to obtain the corresponding high-frequency Rayleigh wave phase velocity distribution; the high-frequency Rayleigh wave phase velocity distribution includes the Rayleigh wave high-frequency phase velocity and its azimuthal anisotropy distribution characteristics.

[0024] The ZZ component NCF waveforms between each station are calculated to obtain the short-period Rayleigh surface wave dispersion information. In order to ensure that the 2-10s Rayleigh surface wave signal can be extracted, the 10Hz continuous noise data is used to perform cross-correlation calculations every hour, and then superimposed to obtain the NCF of one day. Finally, the high-frequency noise cross-correlation function (NCF) of the ZZ component of each day is superimposed to improve the signal-to-noise ratio, such as the NCF superposition processing algorithm proposed by Li et al in 2018. Under normal circumstances, 200-250 days of superimposed NCF can obtain a relatively stable empirical Green's function. Since the station recording time is more than 2 years, it is sufficient to obtain high-quality NCF.

[0025] Step S2 includes the following steps: S2-1. Cut the 10 Hz continuous noise data into 1-hour segments every day, perform cross-correlation calculation and superposition on the Z component continuous waveform data of each station pair every hour, and obtain the ZZ component background noise cross-correlation function between each station pair on the same day, i.e., NFC (Noise Cross-Correlation Function); S2-2, perform superposition processing on each NFC to obtain the ZZ component NCF waveform data between each station pair during the entire observation time; S2-3, picking up the Rayleigh wave dispersion information between the station pairs corresponding to the ZZ component NCF waveform data; S2-4. The Eikonal imaging method is used to process the Rayleigh wave dispersion information and obtain the high-frequency Rayleigh wave phase velocity distribution. The Eikonal imaging method can not only obtain the two-dimensional distribution of surface wave phase velocity and anisotropy, but also take into account the possible ray path bending phenomenon in the path direction compared to the traditional method, thereby obtaining more accurate measurement results.

[0026] S3. The preprocessed seismic data are processed using the FDFA algorithm and the HVSR method to obtain the corresponding Rayleigh wave ellipticity. The FDFA algorithm can filter out the noise fragments dominated by the Rayleigh wave energy, effectively solving the difficulty of extracting the Rayleigh wave ellipticity from a single waveform and avoiding interference from other types of waves (such as body waves, Love waves, etc.).

[0027] Step S3 includes the following steps: S3-1, perform time window segmentation processing on the three-component continuous waveform data to obtain segmented three-component continuous waveform data; wherein the time window segmentation processing is the same as step S3-2, using the FDFA algorithm to calculate the three-component continuous waveform data after segmentation in each time window, and obtain the phase difference between the horizontal component and the vertical component in each time window; S3-3, based on the phase difference between the horizontal component and the vertical component in each time window, the time windows are screened, and the time windows with phase differences in the range of 75° to 115° are retained to obtain screened time windows; and the screened time windows are used as the band dominated by Rayleigh waves; like Figure 2 As shown in (a), the 24-hour three-component continuous waveform includes two earthquake events, one at 9 hours and the other at 17 hours. The continuous waveform is divided into 5-minute time windows, and then the phase difference between the horizontal component and the vertical component in each time window is calculated by the FDPA method, as shown in Figure 2 As shown in (b). It can be seen that the phase difference in different time windows is quite different, especially the phase difference in the time window where the two earthquake events are located varies greatly, which is consistent with the characteristics of body waves. The time window with a phase difference in the range of 75° to 115° is selected as the band dominated by Rayleigh waves, and then the Rayleigh wave ellipticity is calculated by the HVSR method; among them, Figure 2 The area between the two horizontal dashed lines in (b) is the selected dominant band of Rayleigh waves. Figure 2 As shown in (c), there is a significant difference between the results of using the FDPA method to screen the band dominated by Rayleigh waves and those without screening, indicating that the FDPA method can effectively screen out Rayleigh waves and avoid interference from other types of waves (such as body waves, Love waves, etc.). Figure 2 The HSVR corresponding to the black curve in (c) refers to the ellipticity at each period obtained by directly calculating the ellipticity using the HVSR method, while Figure 2 The FDPA-HSVR corresponding to the red curve in (c) refers to the ellipticity at each period obtained by calculating the ellipticity using this method.

[0028] S3-4. Use the HVSR method to calculate the band dominated by Rayleigh waves and obtain the corresponding Rayleigh wave ellipticity.

[0029] The calculation formula of the Rayleigh wave ellipticity in step S3-4 is:

[0030] in, represents the frequency corresponding to the band dominated by Rayleigh waves, represents the Rayleigh wave ellipticity, , Represents two orthogonal horizontal components at frequencies The spectrum of Indicates the vertical component at frequency spectrum.

[0031] S4, performing data processing and screening on the pre-processed seismic data to obtain a corresponding teleseismic cross-convolution waveform; In order to obtain a higher resolution and more reliable shallow velocity structure model in North China, body wave data can be introduced to analyze the shallow structure and The body wave data include teleseismic cross-convolution waveform, body wave amplitude ratio and receiver function. In areas with thick sedimentary layers, the receiver function is strongly disturbed by multiple waves of the sedimentary layer, which will increase the instability of the receiver function deconvolution calculation. In addition, the teleseismic cross-convolution waveform avoids the deconvolution of the deconvolution calculation and is more stable, so the teleseismic cross-convolution waveform is selected.

[0032] Step S4 includes the following steps: S4-1, performing time window processing on the teleseismic event data to obtain processed teleseismic event data; S4-2, screening and superimposing the processed teleseismic event data to obtain superimposed teleseismic event data; The specific process of step S4-2 is: The ENZ three-component data of the processed teleseismic event data are rotated using the rotation matrix to obtain the RTZ component, and the P-wave R component data and Z component data are extracted; The P-wave R component data and Z component data with waveform similarity not less than the threshold and signal-to-noise ratio greater than 5 are superimposed to obtain the body wave stable waveform of each station, that is, the superimposed teleseismic event data.

[0033] S4-3. Perform a cross-convolution waveform operation on the superimposed teleseismic event data to obtain a corresponding teleseismic cross-convolution waveform, and construct a corresponding cross-convolution error function.

[0034] The teleseismic cross-convolution waveform is obtained by cross-convolution of the teleseismic event waveform of the vertical component and the radial component with the theoretical waveform. That is, the calculation formula of the teleseismic cross-convolution waveform in step S4-3 is:

[0035]

[0036]

[0037] in, represents the teleseismic cross-convolution waveform, , represent the vertical and radial components of the teleseismic cross-convolution waveform, respectively. Indicates the time, , Respectively indicate time No. The theoretical vertical waveform and theoretical radial waveform of the earth model, represents the total number of earthquake events, represents the summation function, , denote the vertical and radial Green's functions of the Earth, respectively, , Respectively represent The instrument response and source time function of each earthquake event, , Respectively represent The vertical and radial component waveforms of the earthquake event. When the model approaches the real earth model, the two cross convolution functions in the above two calculation formulas tend to be consistent.

[0038] Since the Green's function of teleseismic events (30° to 90°) remains almost unchanged in the sedimentary layer, all teleseismic event waveforms recorded by the station are superimposed to obtain a stable waveform, and then two cross-convolution functions are calculated to improve the efficiency of joint inversion. At the same time, the cross-convolution error function is defined The calculation formula is:

[0039] in, , Respectively represent the vertical and radial component waveforms after superposition, , They represent the time window ranges of direct P waves and multiple reflection waves from the sediment layer, , Respectively indicate time No. The theoretical radial waveform and theoretical vertical waveform after superposition of the models, represents the absolute value, Represents the integral function.

[0040] S5. Based on the high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity and teleseismic cross-convolution waveform, a shallow velocity structure model is constructed using joint inversion. Joint inversion of multiple data sets can achieve sensitivity complementarity, thereby obtaining a higher-resolution velocity model.

[0041] Step S5 includes the following steps: S5-1. Determine the weight parameters of high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity and teleseismic cross-convolution waveform according to the research objectives; S5-2, constructing a joint inversion error function based on the weight parameters of step S5-1; Joint inversion error function The calculation formula is:

[0042] in, , denote the measured cycle numbers of Rayleigh wave ellipticity and high-frequency Rayleigh wave phase velocity, respectively. , , represent the uncertainties corresponding to the high-frequency Rayleigh wave phase velocity, Rayleigh wave ellipticity and teleseismic cross-convolution waveform, respectively. , denote the observed high-frequency Rayleigh wave phase velocity and Rayleigh wave ellipticity, respectively. , represent the high-frequency Rayleigh wave phase velocity and Rayleigh wave ellipticity obtained in step S2 and step S3 respectively, Represents a sum function.

[0043] S5-3. Obtain the shear wave velocity of the sedimentary layer and upper crust in the area to be measured and longitudinal wave velocity ; S5-3, based on the joint inversion error function, the thickness of the sedimentary layer, and As inversion parameters, the Bayesian Monte Carlo method is used to jointly invert the high-frequency Rayleigh wave phase velocity, high-frequency Rayleigh wave ellipticity and teleseismic cross-convolution waveform. 400,000 forward models are randomly generated within the weight parameter range, and the error function value of each forward model is calculated. Finally, only the 1,000 models with the smallest error value are retained. Then, the average value of these 1,000 models is calculated to obtain the one-dimensional and Speed ​​structure; S5-4. One-dimensional interpolation of each station by grid interpolation within the region and The velocity structure is processed to obtain the shallow three-dimensional and Velocity structure model.

[0044] In the actual inversion process, the sedimentary layer and the upper crust are divided into several small layers according to actual needs, and the shear wave velocity in each small layer is It can be determined by four parameters: top shear wave velocity , bottom shear wave velocity , longitudinal and transverse wave velocity ratio and thickness. Assuming the shear wave velocity of the sedimentary layer increases linearly with depth, while the shear wave velocity in the crust It is obtained by interpolation using the B-spline algorithm.

[0045] In summary, the present invention introduces teleseismic cross-convolution waveforms in the joint inversion, and takes into account the influence of parameters on the velocity structure model, avoids the interference of other types of waves on Rayleigh waves, effectively screens out the noise fragments dominated by Rayleigh wave energy, and effectively solves the difficulty of extracting Rayleigh wave ellipticity from a single waveform; through time inversion, the shallow three-dimensional of the entire crust can be effectively simulated. and The velocity structure model, especially for sedimentary layers, is beneficial to the research of deep structural scientific issues and the exploration and development of oil and gas resources in the study area. It can also improve the accuracy of seismic motion simulation and assist in earthquake-resistant disaster reduction work in earthquake-prone areas.

Claims

1. A method for constructing a shallow velocity model based on joint inversion, characterized in that: The following steps are involved: S1, obtaining seismic data from fixed stations and mobile stations in the area to be measured and preprocessing them to obtain preprocessed seismic data; S2, constructing the background noise cross-correlation function based on the preprocessed seismic data, and using the Eikonal imaging method to obtain the corresponding high-frequency Rayleigh wave phase velocity distribution; S3, using the FDFA algorithm and the HVSR method to process the preprocessed seismic data to obtain the corresponding Rayleigh wave ellipticity; S4, performing data processing and screening on the pre-processed seismic data to obtain a corresponding teleseismic cross-convolution waveform; S5. Based on the high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity and teleseismic cross-convolution waveform, a shallow velocity structure model is constructed using joint inversion.

2. The method for constructing a shallow velocity model based on joint inversion of high-frequency surface waves and teleseismic body waves according to claim 1, characterized in that: The preprocessed seismic data in step S1 include Z-component continuous waveform data, three-component continuous waveform data and teleseismic event data; wherein, the preprocessing of the Z-component continuous waveform data and the three-component continuous waveform data adopts resampling, removal of seismic events and removal of mean, and the preprocessing of the teleseismic event data adopts resampling and removal of mean; the teleseismic event data are waveform data of seismic events with a magnitude of not less than 5.5 and a distance from the epicenter of 30° to 90°.

3. The method for constructing a shallow velocity model based on joint inversion of high-frequency surface waves and teleseismic body waves according to claim 2, characterized in that: The step S2 comprises the following steps: S2-1, cut the 10 Hz continuous noise data into 1 hour segments every day, perform cross-correlation calculation and superposition on the Z component continuous waveform data of each station pair every hour, and obtain the ZZ component background noise cross-correlation function between each station pair on the same day, i.e., NFC; S2-2, perform superposition processing on each NFC to obtain the ZZ component NCF waveform data between each station pair during the entire observation time; S2-3, picking up the Rayleigh wave dispersion information between the station pairs corresponding to the ZZ component NCF waveform data; S2-4. Use the Eikonal imaging method to process the Rayleigh wave dispersion information and obtain the high-frequency Rayleigh wave phase velocity distribution.

4. The method for constructing a shallow velocity model based on joint inversion of high-frequency surface waves and teleseismic body waves according to claim 2, characterized in that: The step S3 comprises the following steps: S3-1, performing time window segmentation processing on the three-component continuous waveform data to obtain segmented three-component continuous waveform data; S3-2, using the FDFA algorithm to calculate the three-component continuous waveform data after segmentation in each time window, and obtain the phase difference between the horizontal component and the vertical component in each time window; S3-3, based on the phase difference between the horizontal component and the vertical component in each time window, the time windows are screened, and the time windows with phase differences in the range of 75° to 115° are retained to obtain screened time windows; and the screened time windows are used as the band dominated by Rayleigh waves; S3-4. Use the HVSR method to calculate the band dominated by Rayleigh waves and obtain the corresponding Rayleigh wave ellipticity.

5. The method for constructing a shallow velocity model based on joint inversion of high-frequency surface waves and teleseismic body waves according to claim 4, characterized in that: The calculation formula of the Rayleigh wave ellipticity in step S3-4 is: in, represents the frequency corresponding to the band dominated by Rayleigh waves, represents the Rayleigh wave ellipticity, , Represents two orthogonal horizontal components at frequency The spectrum of Indicates the vertical component at frequency spectrum.

6. The method for constructing a shallow velocity model based on joint inversion of high-frequency surface waves and teleseismic body waves according to claim 2, characterized in that: The step S4 comprises the following steps: S4-1, performing time window processing on the teleseismic event data to obtain processed teleseismic event data; S4-2, screening and superimposing the processed teleseismic event data to obtain superimposed teleseismic event data; S4-3. Perform a cross-convolution waveform operation on the superimposed teleseismic event data to obtain a corresponding teleseismic cross-convolution waveform, and construct a corresponding cross-convolution error function.

7. The method for constructing a shallow velocity model based on joint inversion of high-frequency surface waves and teleseismic body waves according to claim 6, characterized in that: The specific process of step S4-2 is as follows: The ENZ three-component data of the processed teleseismic event data are rotated using the rotation matrix to obtain the RTZ component, and the P-wave R component data and Z component data are extracted; The P-wave R component data and Z component data with waveform similarity not less than the threshold and signal-to-noise ratio greater than 5 are superimposed to obtain the body wave stable waveform of each station, that is, the superimposed teleseismic event data.

8. The method for constructing a shallow velocity model based on joint inversion of high-frequency surface waves and teleseismic body waves according to claim 6, characterized in that: The calculation formula of the teleseismic cross-convolution waveform in step S4-3 is: in, represents the teleseismic cross-convolution waveform, , represent the vertical and radial components of the teleseismic cross-convolution waveform, respectively. Indicates the time, , Respectively indicate time No. The theoretical vertical waveform and theoretical radial waveform of the earth model, represents the total number of earthquake events, represents the summation function, , denote the vertical and radial Green's functions of the Earth, respectively, , Respectively represent The instrument response and source time function of each earthquake event, , Respectively represent Vertical and radial component waveforms of a seismic event.

9. The method for constructing a shallow velocity model based on joint inversion of high-frequency surface waves and teleseismic body waves according to claim 2, characterized in that: The step S5 comprises the following steps: S5-1. Determine the weight parameters of high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity and teleseismic cross-convolution waveform according to the research objectives; S5-2, constructing a joint inversion error function based on the weight parameters of step S5-1; S5-3. Obtain the shear wave velocity of the sedimentary layer and upper crust in the area to be measured and longitudinal wave velocity ; S5-3, based on the joint inversion error function, the thickness of the sedimentary layer, and The high-frequency Rayleigh wave phase velocity, high-frequency Rayleigh wave ellipticity and teleseismic cross-convolution waveform are jointly inverted using the Bayesian Monte Carlo method as inversion parameters to obtain the one-dimensional and Speed ​​structure; S5-4. One-dimensional interpolation of each station by grid interpolation within the region and The velocity structure is processed to obtain the shallow three-dimensional and Velocity structure model.

Images