A joint inversion method based on surface wave and direct p-wave

Abstract

The application discloses a joint inversion method based on surface waves and direct P waves, relates to the field of seismic tomography, and comprises data acquisition and processing, separate surface wave and P wave model construction, joint inversion algorithm design and inversion result visualization. The surface wave part adopts a frequency-vector wavenumber domain transformation method to extract a dispersion curve, and solves a model correction amount optimization update of a subsurface medium shear wave velocity through Voronoi inversion; the direct P wave part uses a fast marching method in a PYGIMLI program to calculate a ray path and solves a longitudinal wave velocity through a regularization inversion. Finally, through the joint inversion algorithm, the results of the two wave types are comprehensively analyzed and processed, the velocity and geological layering parameters of the subsurface medium in the road collapse area are accurately estimated, the advantages of the surface wave inversion and the direct P wave inversion are combined, the seismic wave data are comprehensively analyzed and processed, the limitations of the single inversion technology are overcome, and the data noise resistance, the inversion accuracy and the resolution of the subsurface 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 in 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:

[0006] S1: Acquire seismic data for the study area, including surface wave background noise data and seismic records directly reaching P waves;

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

[0008] 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.

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

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

[0011] 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.

[0012] Preferably, the preprocessing of the surface background noise data in S2 specifically includes: segmenting the original noise data into one-minute data segments, performing resampling, bandpass filtering, normalization and spectral whitening processing in sequence, and calculating the cross-correlation function between detectors based on the preprocessed data.

[0013] Preferably, in S3, the surface wave inversion 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 amount of the shear wave velocity model, continuously updating the shear wave velocity model to obtain the final shear wave velocity model and dispersion residual.

[0014] Preferably, the formula for surface wave inversion in S3 is expressed as follows:

[0015] ;

[0016] in, This is 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.

[0017] Preferably, the dispersion residual is expressed as:

[0018] ;

[0019] in, For shear wave velocity model; For the observation dispersion of surface waves; For the first Theoretical dispersion of forward iteration.

[0020] Preferably, 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 obtaining the final P-wave velocity model and travel time residual.

[0021] The preferred theoretical travel time for a direct P-wave is expressed as follows:

[0022] ;

[0023] in, The slowness of the initial P-wave velocity model. For the ray path; This is the initial value for the ray path.

[0024] Preferably, the travel time residual is expressed as:

[0025] ;

[0026] 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.

[0027] Preferably, the objective function in S5 is expressed as:

[0028] ;

[0029] in, To obtain the L2 norm; A These are weighting coefficients; For the target model; The prior reference model for the target model; is the damping factor.

[0030] Preferably, in S6, the iteration terminates if one of the following conditions is met:

[0031] (1) The absolute value of the difference between the objective function values ​​of two consecutive iterations is less than a preset threshold;

[0032] (2) The number of iterations reaches the preset maximum number of iterations.

[0033] Therefore, this application employs a joint inversion method based on surface waves and direct P-waves. By combining surface wave and direct P-wave data and comprehensively analyzing and processing seismic wave data, it can more comprehensively estimate the physical parameters of the subsurface medium in road collapse areas, improving data noise resistance, inversion result accuracy, and resolution of shallow subsurface structures. It can provide reliable subsurface structure information under complex geological conditions and is particularly suitable for road collapse monitoring and early warning. Compared to traditional single-waveform inversion methods, this application offers higher inversion accuracy and applicability. Attached Figure Description

[0034] Figure 1 This is a flowchart illustrating a joint inversion method based on surface waves and direct P-waves in this application.

[0035] Figure 2 This is a schematic diagram of the surface wave inversion results in an embodiment of this application;

[0036] Figure 3 This is a schematic diagram of the inversion results of the direct P-wave in the embodiments of this application, where (a) is a P-wave velocity model diagram with ray paths and (b) is a P-wave velocity model diagram.

[0037] Figure 4 This is a schematic diagram of the joint inversion results in an embodiment of this application;

[0038] Figure 5 This is a comparative model diagram in the embodiments of this application. Detailed Implementation

[0039] 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.

[0040] 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.

[0041] The terms "comprising" or "including," as used in this application, mean that the element preceding the term encompasses the element listed after it, 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.

[0042] Example 1:

[0043] A joint inversion method based on surface waves and direct P-waves, such as Figure 1 As shown, it includes the following steps:

[0044] S1: Acquire seismic data for the study area, including surface wave background noise data and seismic records directly reaching P waves;

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

[0046] The preprocessing of surface wave background noise data includes: segmenting the raw noise data acquired by the geophones into 1-minute data segments, and sequentially performing resampling, bandpass filtering, normalization, and spectral whitening; calculating the cross-correlation function between the geophones based on the preprocessed data; and using the open-source Pyrefra program to filter the seismic records of the same survey line and then picking up the first arrival time (semi-automatic extraction, with the picking error set at two sampling periods) to calculate the "shot-receiver distance-travel time" data.

[0047] 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 longitudinal wave velocity model.

[0048] Based on the geological characteristics of the study area and the layout of the linear array, a two-dimensional initial model of the area was constructed. The initial model includes the velocity and density distribution of the subsurface medium, as well as the fixed positions of the seismic source and detectors.

[0049] Surface wave inversion is based on the surface wave dispersion curve. A sensitivity matrix is ​​constructed, and the least squares method is used iteratively at the Voronoi points to solve for the correction factor of the shear wave velocity model. This process continuously updates the shear wave velocity model, yielding the final shear wave velocity model and dispersion residuals. Figure 2 As shown.

[0050] The formula for surface wave inversion is expressed as:

[0051] ;

[0052] 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.

[0053] The dispersion residual is expressed as:

[0054] ;

[0055] in, For shear wave velocity model; For the observation dispersion of surface waves; For the first Theoretical dispersion of forward iteration.

[0056] In practical applications, the model can be continuously updated by solving the incremental model at the Voronoi point 10 times, and the 10th iteration is finally selected as the final shear wave velocity model and dispersion residual.

[0057] The direct P-wave inversion specifically includes: calculating the ray path using the fast travel method, iteratively calculating the residual between the theoretical and actual travel times using the regularized least squares method, and obtaining the final P-wave velocity model and travel time residuals, such as... Figure 3 As shown.

[0058] The theoretical travel time of a direct P-wave is expressed as follows:

[0059] ;

[0060] in, The slowness of the initial P-wave velocity model. For the ray path; This is the initial value for the ray path.

[0061] The travel time residual is expressed as:

[0062] ;

[0063] 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.

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

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

[0066] 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.

[0067] 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:

[0068] ;

[0069] 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.

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

[0071] .

[0072] 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.

[0073] During the inversion process, the objective function parameters are optimized by solving the least squares method, and the model parameters are updated step by step to keep the error within an acceptable range. Information such as the velocity and density of the subsurface medium is then obtained, and the joint inversion results are as follows: Figure 4 As shown.

[0074] Specifically, the criterion for determining the convergence of the objective function is: the difference between the objective function values ​​of two consecutive iterations is less than a preset threshold. ( (or the number of iterations reaches the preset maximum value), ( Stop iterating when the value is ≥30.

[0075] A velocity profile of the underground structure is generated based on the joint inversion results.

[0076] Example 2:

[0077] This embodiment designs a comparative model to verify the effectiveness of the method provided in this application. Specific results are as follows: Figure 5 As shown, by Figure 5 It can be seen that, compared with the traditional single-waveform inversion method, the method provided in this application has higher inversion accuracy and applicability.

[0078] Therefore, this application adopts the above-mentioned joint inversion method based on surface waves and direct P waves. By combining the advantages of surface waves and P waves in their respective characteristics, and through comprehensive analysis and processing of seismic wave data, the noise resistance of the data, the accuracy of the inversion, and the resolution of shallow underground structures are improved.

[0079] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of this application and not to limit them. Although this application has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of this application, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of this 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. The surface wave inversion 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. The direct P-wave inversion 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; 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; The objective function is expressed as: ； 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; The constraint relationship between the shear wave velocity model and the longitudinal wave velocity model is expressed as follows: ； For the shear wave velocity model, For the longitudinal wave velocity model; 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 formula for surface wave inversion in S3 is expressed as follows: ； in, This is 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.

4. The joint inversion method based on surface waves and direct P-waves as described in claim 3, 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.

5. The joint inversion method based on surface waves and direct P-waves as described in claim 4, 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.

6. The joint inversion method based on surface waves and direct P-waves as described in claim 5, 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.

7. The joint inversion method based on surface waves and direct P-waves as described in claim 6, 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