Joint inversion method based on surface wave and direct P wave

Abstract

The invention discloses a joint inversion method based on a surface wave and a direct P wave, and relates to the field of seismic tomography, and the method comprises the steps: data collection and processing, independent surface wave and P wave model construction, joint inversion algorithm design, and inversion result visualization. For the surface wave part, a frequency dispersion curve is extracted by adopting a frequency-vector wave number domain transformation method, and the shear wave speed of the underground medium is optimized and updated by solving the model correction through Voronoi inversion; the direct P wave part uses a fast marching method in a PYGIMLI program to calculate a ray path and regularization inversion to solve a longitudinal wave speed. Finally, parameters such as speed and geological stratification of an underground medium in a road collapse area are accurately estimated through a joint inversion algorithm and integration of results of the two wave modes, seismic wave data are comprehensively analyzed and processed by combining the advantages of surface wave inversion and direct P wave inversion, the limitation of a single inversion technology is overcome, and the method is suitable for large-scale popularization and application. And the noise immunity of the data, the inversion precision and the resolution of the underground shallow structure are improved.

Description

Technical Field

[0001] This application relates to the field of seismic tomography technology, and in particular to a joint inversion method based on surface waves and direct P waves. Background Technology

[0002] Seismic tomography is a technique that reconstructs subsurface structures by analyzing the propagation characteristics of seismic waves. The surface wave signal affects the shear wave velocity of the medium (…). Highly sensitive to changes, high-resolution models of shallow geological structures can be constructed using the cross-correlation function of background noise over long time series or signals excited by artificial seismic sources. Model. The frequency-Bessel transform method is an advanced technology developed in recent years that can efficiently extract multi-order surface wave information from seismic background noise, significantly improving the accuracy of dispersive imaging and the optimization effect of the array.

[0003] The longitudinal wave velocity of the P-wave directly reaching the subsurface medium ( The surface wave inversion technique is sensitive to structural and stratigraphic interfaces, and velocity structures can be inferred by analyzing travel time data from the source to the receiver. However, surface wave inversion alone has limitations in detail characterization due to its limited vertical resolution and insufficient sensitivity to stratigraphic interfaces. While direct P-wave inversion is highly sensitive to interface location and morphology, its ability to reveal shallow structures and shear wave velocity structures is relatively weak. More importantly, neither method can simultaneously acquire the required data when used alone. and These two key parameters. Summary of the Invention

[0004] To address the aforementioned challenges, this application provides a joint inversion method based on surface waves and direct P-waves. By simultaneously considering the propagation characteristics of surface waves and direct P-waves, and combining the dispersion characteristics of background noise from surface waves with the travel time data of P-waves, a joint inversion is performed to estimate parameters such as the velocity of the subsurface medium and geological stratification in the study area.

[0005] To achieve the above objectives, this application provides a joint inversion method based on surface waves and direct P-waves, comprising the following steps: S1: Acquire seismic data for the study area, including surface wave background noise data and seismic records directly reaching P waves; S2: Preprocess the surface wave background noise data, calculate the cross-correlation function between detectors, and extract the surface wave dispersion curve using the frequency-Bessel transform method; pick up the first arrival time of the direct P-wave seismic record to obtain the travel time data of the direct P-wave. S3: Based on the surface wave dispersion curve and the direct P-wave travel time data, surface wave inversion and direct P-wave inversion are performed respectively to obtain the shear wave velocity model and the P-wave velocity model. S4: taking the grid of the S-wave velocity model as a common grid, mapping the P-wave velocity model to the common grid by interpolation, and constructing a target model parameter space; S5: constructing a joint inversion target function based on the surface wave dispersion data and the direct P-wave travel time data; solving the target model parameters by minimizing the target function, and performing joint inversion of the surface wave and the direct P-wave; S6: determining whether the target function converges, if not, updating the parameters of the target function and re-performing joint inversion until convergence, and ending the inversion.

[0006] Preferably, the specific steps of the surface wave background noise data preprocessing in S2 include: segmenting the original noise data into one-minute data segments, sequentially performing resampling, band-pass filtering, normalization and spectral whitening processing, and calculating the cross-correlation function between the detectors based on the preprocessed data.

[0007] Preferably, the surface wave inversion in S3 is based on the surface wave dispersion curve, a sensitivity matrix is constructed, the least square method is used to iteratively solve the S-wave velocity model correction at the Voronoi point, the S-wave velocity model is constantly updated, and the final S-wave velocity model and dispersion residual are obtained.

[0008] Preferably, the calculation formula of the surface wave inversion in S3 is: ; Wherein, is the sensitivity matrix, is the initial S-wave velocity model of the surface wave inversion, is the theoretical dispersion data obtained by forward calculation based on the initial S-wave velocity model.

[0009] Preferably, the dispersion residual is represented as: ; Wherein, is the S-wave velocity model; is the observed dispersion of the surface wave; is the theoretical dispersion of the n-th forward iteration.

[0010] Preferably, the direct P-wave inversion in S3 specifically includes: using the fast marching method to calculate the ray path, using the regularized least square method to iteratively calculate the residual between the theoretical travel time and the actual travel time, and obtaining the final P-wave velocity model and travel time residual.

[0011] Preferably, the theoretical travel time of the direct P-wave is represented as: ; Wherein, is the slowness of the initial P-wave velocity model, is the ray path;​ Initial value of the ray path.

[0012] Preferably, the travel time residual is represented as: ; wherein, is the travel time residual of the direct P-wave; is the theoretical travel time of the direct P-wave in the n th iteration; is the observed travel time of the direct P-wave; , , is the sensitivity matrix of the direct P-wave.

[0013] Preferably, the objective function in S5 is represented as: ; wherein, is the L2 norm; A is the weighting coefficient; is the target model; is the prior reference model of the target model; is the damping factor.

[0014] Preferably, the judgment criterion of convergence in S6 satisfies one of the following conditions to terminate the iteration: (1) the absolute value of the difference between the objective function values of two consecutive iterations is less than a preset threshold value; (2) the number of iterations reaches a preset maximum number of iterations.

[0015] Therefore, the application adopts the above-mentioned joint inversion method based on surface waves and direct P-waves, which can more comprehensively estimate the physical parameters of the underground medium in the road collapse area by combining the data of surface waves and direct P-waves, comprehensively analyzing and processing the seismic wave data, improving the anti-noise performance of the data, the accuracy of the inversion result, and the resolution of the shallow underground structure. It can provide reliable underground structure information under complex geological conditions, and is particularly suitable for road collapse monitoring and early warning. Compared with the traditional single wave type inversion method, the application has higher inversion accuracy and applicability. BRIEF DESCRIPTION OF DRAWINGS

[0016] Figure 1 is a flowchart of a joint inversion method based on surface waves and direct P-waves in the application; Figure 2 is a schematic diagram of the result of surface wave inversion in the embodiment of the application; Figure 3 is a schematic diagram of the result of inversion of direct P-waves in the embodiment of the application, wherein (a) is a longitudinal wave velocity model diagram with ray paths, and (b) is a longitudinal wave velocity model diagram; Figure 4This is a schematic diagram of the joint inversion results in the embodiments of this application; Figure 5 This is a comparative model diagram in the embodiments of this application. Detailed Implementation

[0017] The following detailed description of the embodiments of this application provided in the accompanying drawings is not intended to limit the scope of the claimed application, but merely to illustrate selected embodiments of the application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without inventive effort are within the scope of protection of this application.

[0018] Unless otherwise defined, the technical or scientific terms used in this application shall have the ordinary meaning as understood by a person of ordinary skill in the art to which this application pertains.

[0019] The terms "comprising" or "including," as used in this application, mean that the element preceding the term encompasses the element listed after the term, and do not exclude the possibility of encompassing other elements as well. The terms "inner," "outer," "upper," and "lower," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, and are only for the convenience of describing this application and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this application. When the absolute position of the described object changes, the relative positional relationship may also change accordingly. In this application, unless otherwise expressly specified and limited, the term "attached," etc., should be interpreted broadly. For example, it can refer to a fixed connection, a detachable connection, or an integral part; it can refer to a direct connection or an indirect connection through an intermediate medium; it can refer to the internal communication of two elements or the interaction relationship between two elements. Those skilled in the art can understand the specific meaning of the above terms in this application according to the specific circumstances.

[0020] Example 1: A joint inversion method based on surface waves and direct P-waves, such as Figure 1 As shown, it includes the following steps: S1: Acquire seismic data for the study area, including surface wave background noise data and seismic records directly reaching P waves; S2: Preprocess the surface wave background noise data, calculate the cross-correlation function between detectors, and extract the surface wave dispersion curve using the frequency-vector wavenumber domain transform method; perform first arrival time picking on the direct P-wave seismic record to obtain the travel time data of the direct P-wave; The specific pre-processing of the surface wave background noise data includes: segmenting the original noise data collected by the detector into 1-minute data segments, sequentially performing resampling, band-pass filtering, normalization and spectral whitening processing, and calculating the cross-correlation function between the detectors based on the pre-processed data. The first arrival time of the seismic record of the same line is picked up after filtering using the open-source Pyrefra program (semi-automatic extraction, picking error selection two sampling periods), and the "offset-time" data is calculated.

[0021] S3: based on the surface wave dispersion curve and the direct P-wave travel time data, respectively, surface wave inversion and direct P-wave inversion are performed to obtain a shear wave velocity model and a P-wave velocity model; Based on the geological characteristics of the research area and the arrangement of the linear array, a two-dimensional initial model of the region is constructed. The initial model includes the velocity, density distribution of the underground medium, and the fixed position of the source and the detector.

[0022] The surface wave inversion is based on the surface wave dispersion curve to construct a sensitivity matrix, and the least squares method is used to iteratively solve the shear wave velocity model correction at the Voronoi point, and the shear wave velocity model is constantly updated to obtain the final shear wave velocity model and dispersion residual error, as shown in Figure 2

[0023] The calculation formula of the surface wave inversion is: ; Wherein, is the sensitivity matrix, is the initial shear wave velocity model of the surface wave inversion, is the theoretical dispersion data obtained based on the initial shear wave velocity model.

[0024] The dispersion residual error is represented as: ; Wherein, is the shear wave velocity model; is the observed dispersion of the surface wave; is the theoretical dispersion of the th iteration.

[0025] In specific applications, the model increment can be solved at the Voronoi point by iterating 10 times to constantly update the model, and finally the 10th iteration is selected as the final shear wave velocity model and dispersion residual error.

[0026] The direct P-wave inversion specifically includes: using the fast marching method to calculate the ray path, using the regularized least squares method to iteratively calculate the residual error between the theoretical travel time and the actual travel time, and obtaining the final P-wave velocity model and travel time residual error, as shown in Figure 3

[0027] ​​The theoretical travel time of a direct P-wave is expressed as follows: ; in, The slowness of the initial P-wave velocity model. For the ray path; This is the initial value for the ray path.

[0028] The travel time residual is expressed as: ; in, For the travel time residual directly reaching the P wave; To reach the P wave The theoretical timekeeping of the next iteration; For observation travel time directly reaching the P-wave; , This is the sensitivity matrix for direct P-wave transmission.

[0029] S4: Using the shear wave velocity model grid as a common grid, the longitudinal wave velocity model interpolation is mapped onto the common grid to construct the target model parameter space; S5: Construct a joint inversion objective function based on surface wave dispersion data and direct P-wave travel time data; solve the target model parameters by minimizing the objective function to perform joint inversion of surface waves and direct P-waves; During the joint inversion process, surface wave data is mainly used to invert shallow structures, while P-wave data is used to detect deep strata interface structures. The velocity, density, and other parameters of the subsurface medium are calculated through iterative optimization.

[0030] Given an objective function The inversion results of surface waves and direct P-waves, along with the regularization term of the target model, are expressed in the following form: ; in, To obtain the L2 norm; A These are weighting coefficients; For the target model; The prior reference model for the target model; As the damping factor, through Solving using the curve method, selecting The inflection point of the curve corresponds to The value is used as the optimal parameter.

[0031] The relationship between the shear wave velocity model and the longitudinal wave velocity model is expressed as follows: .

[0032] S6: judging whether the objective function converges, if not, updating the parameters of the objective function and re-performing joint inversion until convergence, and ending the inversion.

[0033] In the inversion process, the optimization parameters of the objective function are solved by the least square method, and the model parameters are updated gradually to make the error in an acceptable range to solve the velocity, density and other information of the underground medium. The joint inversion result is as shown in Figure 4 .

[0034] Specifically, the determination criterion of the convergence of the objective function is that the difference value of the objective function values of two continuous iterations is less than a preset threshold ( ), or the iteration number reaches a preset maximum value ( ≥30) to stop iteration.

[0035] Based on the joint inversion result, the velocity profile of the underground structure is generated.

[0036] Embodiment 2 This embodiment designs a comparison model to verify the effectiveness of the method provided in the application. The specific results are as shown in Figure 5 . As can be seen from Figure 5 , compared with the traditional single wave type inversion method, the method provided in the application has higher inversion accuracy and applicability.

[0037] Therefore, the application adopts the above-mentioned joint inversion method based on surface wave and direct P wave, which combines the advantages of surface wave and P wave in their respective characteristics, and improves the anti-noise performance of the data, the inversion accuracy and the resolution of the underground shallow structure by comprehensively analyzing and processing the seismic wave data.

[0038] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the application but not to limit it. Although the application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the application can still be modified or replaced by equivalents, and these modifications or replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the application.

Claims

1. A joint inversion method based on surface waves and direct P-waves, characterized in that, Includes the following steps: S1: Acquire seismic data for the study area, including surface wave background noise data and seismic records directly reaching P waves; S2: Preprocess the surface wave background noise data, calculate the cross-correlation function between detectors, and extract the surface wave dispersion curve using the frequency-Bessel transform method; pick up the first arrival time of the direct P-wave seismic record to obtain the travel time data of the direct P-wave. S3: Based on the surface wave dispersion curve and the direct P-wave travel time data, surface wave inversion and direct P-wave inversion are performed respectively to obtain the shear wave velocity model and the P-wave velocity model. S4: Using the shear wave velocity model grid as a common grid, the longitudinal wave velocity model interpolation is mapped onto the common grid to construct the target model parameter space; S5: Construct a joint inversion objective function based on surface wave dispersion data and direct P-wave travel time data; solve the target model parameters by minimizing the objective function to perform joint inversion of surface waves and direct P-waves; S6: Determine whether the objective function has converged. If it has not converged, update the parameters of the objective function and perform joint inversion again until convergence, at which point the inversion ends.

2. The joint inversion method based on surface waves and direct P-waves as described in claim 1, characterized in that, The specific preprocessing of the surface background noise data in S2 includes: segmenting the original noise data into one-minute data segments, performing resampling, bandpass filtering, normalization and spectral whitening in sequence, and calculating the cross-correlation function between detectors based on the preprocessed data.

3. The joint inversion method based on surface waves and direct P-waves as described in claim 2, characterized in that, The surface wave inversion described in S3 is based on the surface wave dispersion curve. A sensitivity matrix is ​​constructed, and the least squares method is used at the Voronoi point to iteratively solve the correction of the shear wave velocity model. The shear wave velocity model is continuously updated iteratively to obtain the final shear wave velocity model and dispersion residual.

4. The joint inversion method based on surface waves and direct P-waves as described in claim 3, characterized in that, The formula for surface wave inversion in S3 is expressed as follows: ； in, For the sensitivity matrix, This is the initial shear wave velocity model for surface wave inversion. The data represents the theoretical dispersion data obtained from forward modeling of the initial shear wave velocity model.

5. The joint inversion method based on surface waves and direct P-waves as described in claim 4, characterized in that, The dispersion residual is expressed as: ； in, For shear wave velocity model; For the observation dispersion of surface waves; For the first Theoretical dispersion of forward iteration.

6. The joint inversion method based on surface waves and direct P-waves as described in claim 1, characterized in that, The direct P-wave inversion in S3 specifically includes: calculating the ray path using the fast travel method, iteratively calculating the residual between the theoretical travel time and the actual travel time using the regularized least squares method, and iteratively updating to obtain the final P-wave velocity model and travel time residual.

7. The joint inversion method based on surface waves and direct P-waves as described in claim 6, characterized in that, The theoretical travel time of a direct P-wave is expressed as follows: ； in, The slowness of the initial P-wave velocity model. For the ray path; This is the initial value for the ray path.

8. The joint inversion method based on surface waves and direct P-waves as described in claim 7, characterized in that, The travel time residual is expressed as: ； in, For the travel time residual directly reaching the P wave; To reach the P wave The theoretical timekeeping of the next iteration; For observation travel time directly reaching the P-wave; , This is the sensitivity matrix for direct P-wave transmission.

9. The joint inversion method based on surface waves and direct P-waves as described in claim 1, characterized in that, The objective function described in S5 is expressed as follows: ； in, To obtain the L2 norm; A These are weighting coefficients; For the target model; The prior reference model for the target model; It is the damping factor.

10. The joint inversion method based on surface waves and direct P-waves as described in claim 9, characterized in that, In S6, the convergence criterion is to terminate the iteration if one of the following conditions is met: (1) The absolute value of the difference between the objective function values ​​of two consecutive iterations is less than a preset threshold; (2) The number of iterations reaches the preset maximum number of iterations.

Images