A method for constructing a shallow velocity model based on joint inversion

Abstract

The present invention discloses a method for constructing a shallow velocity model based on joint inversion, which relates to the technical field of geophysical exploration. In the joint inversion of the present invention, high-frequency Rayleigh wave information and teleseismic cross-convolution waveforms in natural earthquakes are introduced, and the influence of parameters on the velocity structure model is considered. At the same time, the interference of other types of waves on Rayleigh waves is avoided, and noise segments dominated by Rayleigh wave energy are effectively screened out, effectively solving the difficulty of extracting the high-frequency Rayleigh wave ellipticity from single-station waveforms; through joint inversion, the shallow three-dimensional and velocity structure models of the entire crust can be effectively simulated, especially for sedimentary layers, which is not only conducive to the research of scientific problems of deep structures in the study area and the exploration and development of oil and gas resources, but also can improve the accuracy of ground motion simulation and assist in earthquake-resistant and disaster-reduction work in earthquake areas.

Description

Technical Field

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

[0002] Seismic post-disaster investigations have shown that basins have a significant amplifying effect on ground motion. The physical property difference between the low-velocity, loose sedimentary layer and the underlying bedrock is generally considered to be the direct cause of the increased seismic damage. It significantly lengthens the duration of ground motion within the basin and amplifies the amplitude of ground motion. Further research has also shown that both the sedimentary layer thickness and the shallow velocity (P-wave and S-wave) significantly affect strong ground shaking. Among them, the sedimentary layer thickness mainly controls the low-frequency resonance behavior of ground motion, while the shallow velocity mainly affects ground motion at 1 Hz and even higher frequencies. 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, to predict the seismic wave propagation effect caused by the near-surface structure, especially to estimate the relatively lower-frequency (period > 1 s) ground motion in the basin, it is necessary to accurately determine the three-dimensional P-wave and S-wave velocity structures of the sedimentary layer above the bedrock. Previous studies have detected and studied the basin structure based on methods such as active and passive source seismic exploration and gravity, and remarkable results have been achieved. However, affected by the distribution of geophysical exploration points, the existing regional velocity models still cannot meet the needs of seismic ground motion site effect assessment and earthquake prevention and disaster reduction. Especially at the regional scale, the shear wave velocity in the shallow layer (0 - 10 km) and the corresponding P-wave to S-wave velocity ratio have not been well resolved. Summary of the Invention

[0003] Aiming at the above deficiencies in the prior art, the method for constructing a shallow velocity model based on joint inversion provided by the present invention solves the problem that the existing technology cannot well resolve the shear wave velocity in the shallow layer and the corresponding P-wave to S-wave velocity ratio at the regional scale.

[0004] To achieve the above invention purpose, the technical solution adopted by the present invention is as follows:

[0005] A method for constructing a shallow velocity model based on joint inversion is provided, which includes the following steps:

[0006] S1. Obtain the seismic data of fixed stations and mobile stations in the area to be measured and perform preprocessing to obtain the preprocessed seismic data;

[0007] 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;

[0008] S3. Use the FDFA algorithm and the HVSR method to process the preprocessed seismic data to obtain the corresponding Rayleigh wave ellipticity;

[0009] S4. Process and screen the preprocessed seismic data to obtain the corresponding teleseismic cross-correlation waveforms.

[0010] S5. Based on the high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity, and teleseismic cross-correlation waveforms, use joint inversion to construct a shallow velocity structure model.

[0011] Furthermore, the preprocessed seismic data in step S1 includes Z-component continuous waveform data, three-component continuous waveform data, and teleseismic event data; among them, for the preprocessing of Z-component continuous waveform data and three-component continuous waveform data, resampling, removing seismic events, and removing the mean are adopted, and for the preprocessing of teleseismic event data, resampling and removing the mean are adopted; the teleseismic event data are seismic event waveform data with a magnitude of not less than 5.5 and an epicentral distance between 30° and 90°.

[0012] Furthermore, step S2 includes the following steps:

[0013] S2-1. Cut the 10Hz continuous noise data per day into one-hour segments, perform cross-correlation calculation and superposition on the Z-component continuous waveform data of each hour for each pair of stations, and obtain the ZZ-component background noise cross-correlation function between each pair of stations on the same day, that is, NFC (Noise Cross-Correlation Function).

[0014] S2-2. Perform superposition processing on each NFC to obtain the ZZ-component NCF waveform data between each pair of stations during the entire observation time.

[0015] S2-3. Pick up the Rayleigh wave dispersion information between the pairs of stations corresponding to the ZZ-component NCF waveform data.

[0016] S2-4. Use the Eikonal imaging method to process the Rayleigh wave dispersion information to obtain the high-frequency Rayleigh wave phase velocity distribution.

[0017] Furthermore, step S3 includes the following steps:

[0018] S3-1. Perform time window segmentation processing on the three-component continuous waveform data to obtain the segmented three-component continuous waveform data.

[0019] S3-2. Use the FDFA algorithm to calculate the segmented three-component continuous waveform data in each time window respectively to obtain the phase difference between the horizontal component and the vertical component in each time window.

[0020] S3-3. Screen the time windows based on the phase difference between the horizontal component and the vertical component in each time window, retain the time windows with the phase difference in the range of 75° to 115°, and obtain the screened time windows; and use the screened time windows as the frequency bands dominated by Rayleigh waves.

[0021] S3-4. Calculate the corresponding Rayleigh wave ellipticity using the HVSR method for the frequency bands dominated by Rayleigh waves.

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

[0023]

[0024] where, represents the frequency corresponding to the frequency band dominated by Rayleigh waves, represents the Rayleigh wave ellipticity, 、 represent the spectra of two orthogonal horizontal components at frequency and represents the spectrum of the vertical component at frequency .

[0025] Furthermore, step S4 includes the following steps:

[0026] S4-1. Perform time window slicing on the teleseismic event data to obtain the processed teleseismic event data.

[0027] S4-2. Screen and stack the processed teleseismic event data to obtain the stacked teleseismic event data.

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

[0029] Furthermore, the specific process of step S4-2 is:

[0030] Use the rotation matrix to rotate the ENZ three-component data of the processed teleseismic event data to obtain the RTZ components, and extract the P-wave R-component data and Z-component data.

[0031] Stack the P-wave R-component data and Z-component data with the waveform similarity not less than the threshold and the signal-to-noise ratio greater than 5 to obtain the body wave stable waveform of each station, that is, the stacked teleseismic event data.

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

[0033]

[0034]

[0035]

[0036] Among them, represents the teleseismic cross-correlation waveform, , respectively represent the vertical and radial components of the teleseismic cross-correlation waveform, represents the time, , respectively represent the theoretical vertical waveform and theoretical radial waveform of the th Earth model at time represents the total number of seismic events, represents the summation function, , respectively represent the vertical and radial Green's functions of the Earth, , respectively represent the instrumental response and source time function of the th seismic event, , respectively represent the vertical and radial component waveforms of the th seismic event.

[0037] Furthermore, step S5 includes the following steps:

[0038] S5-1. Determine the high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity, and the weight parameters of the teleseismic cross-correlation waveform according to the research objective;

[0039] S5-2. Construct a joint inversion error function based on the weight parameters in step S5-1;

[0040] S5-3. Obtain the shear wave velocity and the P-wave velocity of the sedimentary layer and the upper crust in the area to be measured; S5-3. Based on the joint inversion error function, take the thickness of the sedimentary layer, and as inversion parameters, and use the Bayesian Monte Carlo method to jointly invert the high-frequency Rayleigh wave phase velocity, high-frequency Rayleigh wave ellipticity, and the teleseismic cross-correlation waveform to obtain the one-dimensional and velocity structure of each station;

[0041] S5-4. Construct a joint inversion error function based on the cross-correlation error function;

[0042] S5-5. Through grid interpolation within the region, for the one-dimensional and Process the velocity structure to obtain the shallow three-dimensional of the area to be measured and velocity structure model.

[0043] The beneficial effects of the present invention are as follows: In the joint inversion, teleseismic cross-convolution waveforms are introduced, and the influence of parameters on the velocity structure model is considered. The interference of other types of waves to Rayleigh waves is avoided, and the noise segments dominated by Rayleigh wave energy are effectively screened out, effectively solving the difficulty of extracting the Rayleigh wave ellipticity from single-station waveforms; Through time reversal, the shallow three-dimensional and velocity structure model can be effectively simulated. Especially for the sedimentary layer, it is beneficial to the research of deep structural scientific problems in the study area and the exploration and development of oil and gas resources, and can also improve the accuracy of ground motion simulation and assist in earthquake-resistant and disaster-reduction work in earthquake areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 is a flowchart of the method;

[0045] Figure 2 are the data waveform diagram and curve diagram of the present invention; among them, (a) is a 24-hour three-component continuous waveform diagram; (b) is a phase difference curve diagram of horizontal and vertical components; (c) is a comparison curve diagram of Rayleigh wave ellipticity values before and after screening. DETAILED DESCRIPTION OF THE INVENTION

[0046] The following describes the specific embodiments of the present invention to facilitate those skilled in the art of the present technology to understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art of the present technology, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are within the scope of protection.

[0047] As Figure 1 shown, a method for constructing a shallow velocity model based on joint inversion includes the following steps:

[0048] S1. Obtain the seismic data of fixed stations and mobile stations in the area to be measured and perform preprocessing to obtain the preprocessed seismic data; The area to be measured is the North China region where intracontinental strong earthquakes are very active. In this study, seismic data recorded by 156 fixed stations of a seismic network and 466 broadband mobile stations of the ChinArray III phase in a certain region are used. The recording duration of the mobile seismic stations is about 2 years, which can provide sufficient continuous waveforms and seismic data.

[0049] It is planned to collect and collate drilling data and artificial seismic exploration results in the study area, and select those related to sedimentary layer thickness, shear wave velocity and P-S wave velocity ratio data, and compare this data with the one-dimensional models of nearby stations to verify the reliability of the models.

[0050] The preprocessed seismic data in step S1 includes Z-component continuous waveform data, three-component continuous waveform data, and teleseismic event data; among them, for the preprocessing of the Z-component continuous waveform data and the three-component continuous waveform data, resampling, removing seismic events, and removing the mean are adopted, and for the preprocessing of the teleseismic event data, resampling and removing the mean are adopted; the teleseismic event data is the waveform data of seismic events with a magnitude of not less than 5.5 and an epicentral distance between 30° and 90°.

[0051] S2. Construct a cross-correlation function of ambient noise 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 high-frequency phase velocity of Rayleigh waves and its azimuthal anisotropy distribution characteristics.

[0052] Calculate the ZZ-component NCF waveforms between each pair of stations to obtain the dispersion information of short-period Rayleigh surface waves. To ensure that Rayleigh surface wave signals in the range of 2 - 10 s can be extracted, continuous noise data at 10 Hz is used to perform cross-correlation calculations every hour respectively, and then stacked to obtain the NCF for one day. Finally, the high-frequency noise cross-correlation function (NCF) of the ZZ-component for each day is subjected to corresponding stacking processing to improve the signal-to-noise ratio, such as the stacking processing algorithm of NCF proposed by Li et al in 2018. Usually, a relatively stable empirical Green's function can be obtained from the stacked NCF for 200 - 250 days. Since the recording duration of each station is more than 2 years, it is sufficient to obtain high-quality NCF.

[0053] Step S2 includes the following steps:

[0054] S2-1. Cut the 10-Hz continuous noise data for each day into segments of 1 hour each, perform cross-correlation calculations and stacking on the Z-component continuous waveform data for each hour of each pair of stations, and obtain the ZZ-component ambient noise cross-correlation function between each pair of stations on the same day, that is, NFC (Noise Cross-Correlation Function);

[0055] S2-2. Perform stacking processing on each NFC to obtain the ZZ-component NCF waveform data between each pair of stations during the entire observation time;

[0056] S2-3. Pick up the Rayleigh wave dispersion information between the pairs of stations corresponding to the ZZ-component NCF waveform data;

[0057] S2-4. Use the Eikonal imaging method to process the Rayleigh wave dispersion information to 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 takes into account the possible ray path bending phenomenon in the path direction compared with the traditional method, so as to obtain more accurate measurement results.

[0058] S3. Use the FDFA algorithm and the HVSR method to process the preprocessed seismic data to obtain the corresponding Rayleigh wave ellipticity; the FDFA algorithm can screen out the noise segments dominated by Rayleigh wave energy, effectively solve the difficulty of extracting Rayleigh wave ellipticity from single-station waveforms, and avoid the interference of other types of waves (such as body waves, Love waves, etc.).

[0059] Step S3 includes the following steps:

[0060] S3-1. Perform time window segmentation processing on the three-component continuous waveform data to obtain the segmented three-component continuous waveform data; among them, the time window segmentation processing is related to the steps

[0061] S3-2. Use the FDFA algorithm to calculate the segmented three-component continuous waveform data in each time window respectively to obtain the phase difference between the horizontal component and the vertical component in each time window;

[0062] S3-3. Based on the phase difference between the horizontal component and the vertical component in each time window, screen the time windows, retain the time windows with the phase difference in the range of 75° to 115°, and obtain the screened time windows; and use the screened time windows as the bands dominated by Rayleigh waves;

[0063] As Figure 2 shown in (a) below, the 24-hour three-component continuous waveform, which includes two seismic events, near 9 hours and 17 hours respectively. The continuous waveform is segmented by 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 respectively, as Figure 2 shown in (b) below. It can be seen that the phase differences in different time windows are quite different, especially the phase differences in the time windows where the two seismic events are located change greatly, which conforms to the characteristics of body waves. Screen the time windows with the phase difference in the range of 75° to 115° as the bands dominated by Rayleigh waves, and then calculate the Rayleigh wave ellipticity by the HVSR method; among them, Figure 2 the two horizontal dotted lines in (b) below are the screened bands dominated by Rayleigh waves. As can be seen from Figure 2 shown in (c) below, there are obvious differences between the results of screening the bands dominated by Rayleigh waves by the FDPA method and not screening, indicating that the FDPA method can effectively screen out Rayleigh waves and avoid the interference of other types of waves (such as body waves, Love waves, etc.). Among them, Figure 2The 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 through this method.

[0064] S3-4. Use the HVSR method to calculate the Rayleigh wave-dominated band to obtain the corresponding Rayleigh wave ellipticity.

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

[0066]

[0067] where represents the frequency corresponding to the Rayleigh wave-dominated band, represents the Rayleigh wave ellipticity, , represent the spectra of two orthogonal horizontal components at frequency , represents the spectrum of the vertical component at frequency .

[0068] S4. Perform data processing and screening on the preprocessed seismic data to obtain the corresponding teleseismic cross-correlation waveform;

[0069] To obtain a higher-resolution and more reliable shallow velocity structure model in North China, body wave data can be introduced to provide important constraints on the shallow structure in North China and The body wave data includes teleseismic cross-correlation waveforms, body wave amplitude ratios, and receiver functions. Since in areas with thick sedimentary layers, the receiver functions are strongly interfered by sedimentary multiples, which will increase the instability of receiver function deconvolution calculations; and since teleseismic cross-correlation waveforms avoid the deconvolution of deconvolution calculations and are more stable, teleseismic cross-correlation waveforms are thus selected for introduction.

[0070] Step S4 includes the following steps:

[0071] S4-1. Perform windowing processing on the teleseismic event data to obtain the processed teleseismic event data;

[0072] S4-2. Perform data screening and stacking on the processed teleseismic event data to obtain the stacked teleseismic event data;

[0073] The specific process of step S4-2 is:

[0074] Use the rotation matrix to rotate the ENZ three-component data of the processed teleseismic event data to obtain the RTZ components, and extract the P-wave R-component data and Z-component data;

[0075] Superimpose 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, and obtain the body wave stable waveform of each station, that is, the teleseismic event data after superposition.

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

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

[0078]

[0079]

[0080]

[0081] Where represents the teleseismic cross-convolution waveform, , represent the vertical component and radial component of the teleseismic cross-convolution waveform respectively, represents the time, , represent the theoretical vertical waveform and theoretical radial waveform of the th earth model at time represents the total number of seismic events, represents the summation function, , represent the vertical and radial Green's functions of the earth respectively, , represent the instrument response and source time function of the th seismic event respectively, , represent the vertical and radial component waveforms of the th seismic event respectively. When the th model approaches the true earth model, the two cross-convolution functions in the above two calculation formulas tend to be consistent.

[0082] Since the Green's function of teleseismic events (30° to 90°) remains almost unchanged in the sedimentary layer, the waveforms of all teleseismic events recorded at 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, define the calculation formula for the cross-convolution error function as:

[0083]

[0084] Among them, and respectively represent the vertical and radial component waveforms after superposition, and respectively represent the time window ranges of the direct P-wave and the multiple reflected waves in the sedimentary layer, and respectively represent the th theoretical radial waveform and the theoretical vertical waveform after the superposition of the model, represents the absolute value, represents the integral function.

[0085] S5. Based on the high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity, and teleseismic cross-correlation waveforms, use joint inversion to construct a shallow velocity structure model. Joint inversion of multiple datasets can achieve complementary sensitivities, thereby obtaining a higher-resolution velocity model.

[0086] Step S5 includes the following steps:

[0087] S5-1. According to the research objective, determine the weight parameters of the high-frequency Rayleigh wave phase velocity distribution, Rayleigh wave ellipticity, and teleseismic cross-correlation waveforms;

[0088] S5-2. Based on the weight parameters in step S5-1, construct a joint inversion error function;

[0089] The joint inversion error function is calculated by the formula:

[0090]

[0091] Among them, and respectively represent the number of measurement periods of the Rayleigh wave ellipticity and the high-frequency Rayleigh wave phase velocity, and and respectively represent the uncertainties corresponding to the high-frequency Rayleigh wave phase velocity, Rayleigh wave ellipticity, and teleseismic cross-correlation waveforms, and respectively represent the observed high-frequency Rayleigh wave phase velocity and Rayleigh wave ellipticity, and respectively represent the high-frequency Rayleigh wave phase velocity and Rayleigh wave ellipticity obtained in steps S2 and S3 respectively, represents the summation function.

[0092] S5-3. Obtain the shear wave velocity and the P-wave velocity ; S5-3. Based on the joint inversion error function, taking the thickness of the sedimentary layer, and as inversion parameters, using the Bayesian Monte Carlo method to jointly invert the phase velocity of high-frequency Rayleigh waves, the ellipticity of high-frequency Rayleigh waves and the teleseismic cross-convolution waveform. Randomly generate 400,000 forward models within the weight parameter range, calculate the error function values of each forward model, and finally only retain the 1,000 models with the smallest error values, and then average these 1,000 models to obtain the one-dimensional and velocity structure of each station;

[0093] S5-4. Process the one-dimensional and velocity structure of each station through grid interpolation in the region to obtain the shallow three-dimensional and velocity structure model of the area to be measured.

[0094] In the actual inversion process, the sedimentary layer and the upper crust are divided into several small layers according to actual needs. The shear wave velocity in each small layer can be determined by four parameters: the top shear wave velocity , the bottom shear wave velocity , the shear wave to P-wave velocity ratio and the thickness. Assume that the shear wave velocity of the sedimentary layer increases linearly with depth, while the shear wave velocity in the crust is interpolated according to the B-spline algorithm.

[0095] To sum up, the present invention introduces the teleseismic cross-convolution waveform in the joint inversion, 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 segments dominated by Rayleigh wave energy, and effectively solves the difficulty of extracting the ellipticity of Rayleigh waves from single-station waveforms; through time reversal, it can effectively simulate the shallow three-dimensional and velocity structure model of the entire crust. Especially for the sedimentary layer, it is beneficial to the research of deep structural scientific problems in the study area and the exploration and development of oil and gas resources, and can also improve the accuracy of ground motion simulation and assist in earthquake-resistant and disaster-reduction work in earthquake 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; 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; 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, performing a cross-convolution waveform operation on the superimposed teleseismic event data to obtain a corresponding teleseismic cross-convolution waveform, and constructing a corresponding cross-convolution error function; 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.

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 1, 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 frequencies The spectrum of Indicates the vertical component at frequency spectrum.

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 1, 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.

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 1, 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.

Images