A Joint Tomographic Inversion Method for P-wave and S-wave Kernel Clustering of Shallow Surface Seismic Waves

Abstract

This invention discloses a joint tomographic inversion method for shallow surface seismic P-wave and S-wave nuclear clustering, belonging to the field of geophysical exploration and seismic imaging technology. The invention includes: comprehensively utilizing P-wave and S-wave travel time data to establish a joint inversion model, and introducing nuclear clustering constraints to enhance the spatial consistency of the two wave velocity structures; wherein the P-wave observation system adopts a combination of well-excitation and well-received data and ground-excitation and ground-received data, and the S-wave observation system adopts ground-excitation and ground-received data, realizing multi-scale joint imaging. This invention maps samples to a high-dimensional space through kernel functions, which can strengthen inter-class differences, significantly improve nonlinear clustering capabilities, and enhance the accuracy of the joint P-wave and S-wave inversion model. It features reasonable structural constraints, stable imaging results, and high interface resolution, and is suitable for high-precision velocity structure inversion in shallow surfaces.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration and seismic imaging technology, and relates to a shallow surface seismic P-wave and S-wave joint tomographic inversion method, specifically a shallow surface seismic P-wave and S-wave kernel clustering joint tomographic inversion method. Background Technology

[0002] Shallow surface seismic exploration is an important tool for studying the structure and physical properties of near-surface media, and is widely used in geological disaster prevention, urban foundation evaluation, and the detection of complex geological bodies such as karst and faults. Current technologies typically employ seismic tomography to invert seismic wave travel times to obtain subsurface velocity structures. Based on wave type, it can be divided into two categories: P-wave (longitudinal wave) tomography and S-wave (short-span wave) tomography.

[0003] However, inversion using P-waves or S-waves alone is prone to non-uniqueness and instability: P-wave velocity is sensitive to lithology but struggles to characterize shear properties; S-wave velocity reflects the mechanical properties of the medium but is highly sensitive to noise. To improve model reliability, joint P-wave and S-wave inversion is typically employed, enhancing the stability of the results through complementary constraints between the two types of data. However, existing joint inversion methods often employ constraints such as cross-gradients, resulting in limited structural correlation and difficulty in effectively distinguishing geological interfaces from homogeneous regions. Joint inversion methods based on fuzzy C-means (FCM) clustering are unsatisfactory for parametric clustering of shallow geophysical data because FCM is only applicable to near-circular (spherical) clusters and has poor noise resistance. Mapping samples to a high-dimensional space using kernel functions can strengthen inter-cluster differences, significantly improve nonlinear clustering capabilities, and enhance the accuracy of seismic P-wave and S-wave joint inversion models. Therefore, a method for inverting shallow geophysical data based on kernel clustering is urgently needed. Summary of the Invention

[0004] To address the shortcomings of existing technologies, the present invention aims to provide a method for shallow surface seismic P-wave and S-wave nuclear clustering combined tomographic inversion.

[0005] The objective of this invention is achieved through the following technical solution: a shallow surface seismic P-wave and S-wave nuclear clustering joint tomographic inversion method, comprising the following steps:

[0006] (1) Obtain travel time observation data of seismic P-waves (P-waves) and S-waves (S-waves) in shallow surface areas;

[0007] (2) The observation system based on the earthquake P-wave and S-wave travel time observation data calculates the theoretical travel time field of P-wave and S-wave, and establishes the ray path sensitivity matrix based on the theoretical travel time field to form the travel time equation of P-wave and S-wave.

[0008] (3) Based on the travel time equations of the P-wave and S-wave, a joint inversion objective function of P-wave and S-wave containing a kernel clustering constraint term is constructed, wherein the kernel clustering constraint term is used to constrain the consistency between the P-wave and S-wave velocity structure and the geological unit division.

[0009] (4) Based on the objective function of the joint inversion of P-wave and S-wave, an augmented equation set is established, and the augmented equation set is solved iteratively. The P-wave and S-wave slowness models are updated iteratively until convergence, and the P-wave velocity model and S-wave velocity model are obtained.

[0010] (5) Output P-wave velocity model and S-wave velocity model for shallow surface velocity structure inversion and geological structure interpretation.

[0011] Furthermore, the seismic P-wave and S-wave travel time observation data are obtained using different observation systems. Specifically, the seismic P-wave travel time observation data are obtained by a combination of well-excitation and well-received data and ground-excitation and ground-received data, while the seismic S-wave travel time observation data are obtained by ground-excitation and ground-received data.

[0012] Furthermore, during the ground excitation and ground reception process, in order to ensure that the ground rays can reach the bottom of the probe and obtain good ray coverage, the ratio of the ground survey line length to the probe depth is 3:1-5:1, thereby ensuring the effective inversion of the shallow surface velocity structure.

[0013] Furthermore, the theoretical travel time fields of the longitudinal and transverse waves are calculated using the Fast Marching Method (FMM).

[0014] Furthermore, the construction of the travel time equations is specifically as follows: based on the ray path sensitivity matrices of the P-wave and S-wave, the slowness model increments of the P-wave and S-wave are multiplied by the corresponding sensitivity matrices, so that the results correspond to the travel time residual vectors of the P-wave and S-wave respectively, thereby establishing the linearized travel time equations of the P-wave and S-wave.

[0015] Furthermore, the objective function for the joint inversion of P-waves and S-waves is: based primarily on the fitting of P-wave and S-wave travel time data, and simultaneously includes model regularization constraints and kernel clustering constraints. The kernel clustering constraints are used to maintain the similarity of P-wave and S-wave velocity structures within the same geological unit and to enhance the separability at the interfaces of different geological units.

[0016] Furthermore, the regularization constraint term adopts second-order Laplacian regularization to smooth the velocity field and reduce local oscillations.

[0017] Furthermore, the augmented equations are specifically constructed based on minimizing the objective function of the joint inversion of P-waves and S-waves. The augmented equations include travel time data fitting terms, model regularization terms, and kernel clustering constraint terms.

[0018] Furthermore, the augmented equations are solved using the Least Squares QR iterative algorithm based on the least squares criterion to obtain the update quantities of the P-wave and S-wave slowness models.

[0019] The present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method described thereon.

[0020] The beneficial effects of this invention are as follows: by introducing a kernel clustering constraint term into the P-wave and S-wave joint tomographic inversion model, the P-wave and S-wave velocity structures maintain similarity within the same geological unit and separability at the interface, thereby enhancing the structural consistency and geological interpretability of the joint inversion results. Furthermore, this method can simultaneously process P-wave data from cross-aperture excitation-cross-aperture reception and surface excitation-surface reception, and, combined with surface S-wave observations, achieve multi-scale, multi-system joint imaging, making it suitable for high-precision inversion of shallow surface velocity structures and identification of geological interfaces. Attached Figure Description

[0021] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0022] Figure 1 This is a flowchart of the method of the present invention;

[0023] Figure 2 This is a schematic diagram of the Vp real model provided in an embodiment of the present invention;

[0024] Figure 3 This is a schematic diagram of the Vs real model provided in an embodiment of the present invention;

[0025] Figure 4 This is a schematic diagram of RMS error convergence for P-wave data provided in an embodiment of the present invention.

[0026] Figure 5 This is a schematic diagram of RMS error convergence for shear wave data provided in an embodiment of the present invention.

[0027] Figure 6 This is a schematic diagram of the Vp inversion results of kernel clustering constraints provided in an embodiment of the present invention;

[0028] Figure 7 This is a schematic diagram of the Vs inversion results of kernel clustering constraints provided in an embodiment of the present invention;

[0029] Figure 8This is a schematic diagram of the Vp inversion results without kernel clustering constraints provided in an embodiment of the present invention;

[0030] Figure 9 This is a schematic diagram of the Vs inversion results without kernel clustering constraints provided in an embodiment of the present invention;

[0031] Figure 10 This is a schematic diagram of the Vp inversion result of FCM clustering constraints provided in an embodiment of the present invention;

[0032] Figure 11 This is a schematic diagram of the Vs inversion results of FCM kernel clustering constraints provided in an embodiment of the present invention. Detailed Implementation

[0033] To enable those skilled in the art to better understand the technical solutions in the embodiments of this application, the technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art should fall within the protection scope of the embodiments of this application.

[0034] This invention provides a method for shallow surface seismic P-wave and S-wave nuclear clustering joint tomography inversion, such as... Figure 1 As shown, it includes the following steps:

[0035] Step 1: Obtain travel time observation data of seismic P-waves (P-waves) and S-waves (S-waves) in shallow surface areas;

[0036] Step 2: The observation system calculates the theoretical travel time fields of the P-wave and S-wave based on the earthquake P-wave and S-wave travel time observation data, and establishes the ray path sensitivity matrix based on the theoretical travel time fields to form the travel time equations of the P-wave and S-wave.

[0037] Step 3: Based on the travel time equations of the P-wave and S-wave, construct a joint inversion objective function for P-wave and S-wave that includes a kernel clustering constraint term, wherein the kernel clustering constraint term is used to constrain the consistency between the P-wave and S-wave velocity structure and the geological unit division.

[0038] Step 4: Based on the joint inversion objective function of P-wave and S-wave, establish an augmented equation set and iteratively solve the augmented equation set, iteratively updating the P-wave and S-wave slowness models until convergence, to obtain the P-wave velocity model and the S-wave velocity model.

[0039] Step 5: Output the P-wave velocity model and S-wave velocity model for shallow surface velocity structure inversion and geological structure interpretation.

[0040] In a preferred embodiment, step one specifically involves: establishing as follows: Figure 2 and Figure 3 The Vp and Vs models shown were tested, among which... Figure 2 This represents a realistic Vp model with three geological layers and velocities of 400m / s, 700m / s, and 1000m / s, including undulating bedrock faults and other structural features. Figure 3 Represents the Vs real model, and Figure 2 The Vp real models have the same underground structure and stratum distribution, with velocities of 200m / s, 388m / s, and 600m / s, respectively.

[0041] Travel time observation data of seismic P-waves (P-waves) and S-waves (S-waves) in shallow surface areas were acquired, and 1% Gaussian noise was added to simulate the real environment. The P-wave observation system employed a combination of well-excitation and well-received data, and the S-wave observation system used a ground-excitation and ground-received method. During the ground-excitation and ground-received process, to ensure that the ground rays reached the bottom of the probe and achieved good ray coverage, the survey line length for both models was 120m, and the probe depth was 24m, i.e., the ratio of the ground survey line length to the probe depth was 5:1, thus ensuring the effective inversion of the shallow surface velocity structure.

[0042] As a preferred embodiment, step two specifically involves: using a linearly increasing velocity model as the initial inversion model, with the velocity range of the Vp initial inversion model being 400-1000 m / s and the velocity range of the Vs initial inversion model being 200-600 m / s; calculating the theoretical travel time fields of the P-wave and S-wave using the fast travel method; establishing the ray path sensitivity matrix; and forming the travel time equations of the P-wave and S-wave, thus providing a forward modeling basis for subsequent joint inversion.

[0043] For each waveform, the travel time residual can be expressed as:

[0044] ;

[0045] in, and These are the ray path sensitivity matrices for longitudinal and transverse waves, respectively; , These are the increments for the slowness (reciprocal of velocity) models of P-waves and S-waves, respectively; This is the corresponding travel time residual vector.

[0046] To obtain a high-precision and stable travel time field, this invention employs the Fast Marching Method (FMM) to calculate the theoretical travel time of each grid point. The FMM is based on the Eikonal equations:

[0047] ;

[0048] in, To make a appearance, It is a slow-degree field.

[0049] By employing a one-way propagation update strategy, the fast travel method can obtain a stable travel time distribution under a given velocity model, which can be used to construct the ray path sensitivity matrix. and The forward modeling system established in this way provides the basis for data fitting terms for the joint inversion objective function, enabling the synchronous inversion of P-wave and S-wave travel time data.

[0050] As a preferred embodiment, the objective function for the joint P-wave and S-wave inversion defined in step three is specifically as follows:

[0051] ;

[0052] Among them, the first item Item 2 These represent the travel time fitting errors for the longitudinal and transverse waves, respectively; the third term... This is the regularization term; the fourth term. This is a kernel clustering constraint term used to maintain structural consistency between P-wave and S-wave models. and These are the weight coefficients for regularization constraints and clustering constraints, respectively.

[0053] Preferably, the regularization constraint term adopts second-order Laplacian regularization, which can further smooth the velocity field and reduce local oscillations, as shown in the following formula:

[0054] ;

[0055] The kernel clustering constraint term achieves adaptive clustering of the model by calculating the similarity between the P-wave and S-wave models in a high-dimensional kernel space. It is defined as follows:

[0056] ;

[0057] Where N is the number of model units; C is the number of cluster categories; is the membership degree of the i-th unit belonging to the j-th class; m is the fuzzy index (usually taken as 2); The Gaussian kernel is used as the kernel mapping function to measure the nonlinear similarity between P-wave and S-wave models. The Gaussian kernel is expressed as follows:

[0058] ;

[0059] in, This represents two data points (which can correspond to the feature vectors of two model units in the inversion). ; This represents the Euclidean distance between the two points in the original feature space; The kernel width parameter controls the rate at which similarity decays. The value ranges from 0 to 1, and the closer the value is to 1, the more similar the two points are.

[0060] Therefore, when and The closer they get, The smaller the value, the closer the exponent term is to 1, indicating that the two points belong to similar geological structures; conversely, the larger the distance, the closer the exponent term is to 1. The larger the value, the closer the kernel value is to 0, indicating a significant difference between the two. This constraint term applies clustering coupling to the P-wave and S-wave velocity fields during the inversion process, making the P-wave and S-wave velocity inversion models converge within the same geological unit, while maintaining significant differences at geological interfaces.

[0061] As a preferred embodiment, step four specifically involves: using the LSQR (Least Squares QR) algorithm to solve the augmented equations and iteratively updating the P-wave and S-wave slowness models until convergence;

[0062] The augmented equation that needs to be solved can be expressed as:

[0063] ;

[0064] in, For regularized difference operators, This is a coupling operator derived from the clustering membership matrix.

[0065] In this embodiment, the maximum number of iterations is set to 30. The iteration process is as follows: Figure 4 and Figure 5 As shown, Figure 4 The variation of P-wave data error with the number of iterations, Figure 5 The RMS error of the shear wave data changes with the number of iterations. As the iterations proceed, the RMS error of both types of data decreases and eventually stabilizes.

[0066] In a preferred embodiment, step five specifically involves: outputting the longitudinal wave velocity model and the transverse wave velocity model as shown below. Figure 6 and Figure 7 As shown, Figure 6 The Vp inversion results are based on kernel clustering constraints. Figure 7 The Vs inversion results with kernel clustering constraints show that the subsurface structure and stratigraphic distribution of both are highly consistent with the real model, and can be used for shallow surface velocity structure inversion and geological structure interpretation.

[0067] As a preferred embodiment, to verify the effectiveness of the method of the present invention, the joint tomographic inversion results of P-waves and S-waves with kernel clustering constraints are compared with the joint inversion results without kernel clustering constraints. Figure 8 and Figure 9As shown, Figure 8 , Figure 9 The images show the Vp and Vs inversion results without kernel clustering constraints, respectively. The velocity distribution within the geological unit is relatively discrete, and the interface location and morphology are not clear enough. In contrast, the inversion results with kernel clustering constraints show more consistent P-wave and S-wave velocity distributions within the same geological unit, clearer structural boundaries, and interface thickness variations that better reflect geological realities.

[0068] As a further preferred embodiment, the method of the present invention is compared with the traditional fuzzy C-means (FCM) clustering joint inversion of P- and S-waves. Figure 10 and Figure 11 As shown, in the Vp and Vs inversion results obtained under FCM constraints, due to the limited adaptability of FCM to irregular structures and noise, significant velocity differences still exist within some geological units. In contrast, kernel clustering can better maintain the consistency within homogeneous bodies and enhance the separability between different tectonic units. The above comparisons demonstrate that the P-wave and S-wave kernel clustering combined tomographic inversion method proposed in this invention significantly improves the accuracy and reliability of shallow surface velocity structure inversion.

[0069] This invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method described thereon.

[0070] The computer-readable storage medium can be an internal storage unit of any data processing device described in any of the foregoing embodiments, such as a hard disk or memory. The computer-readable storage medium can also be any data processing device, such as a plug-in hard disk, smart media card (SMC), SD card, flash card, etc., equipped on the device. Furthermore, the computer-readable storage medium can include both internal storage units of any data processing device and external storage devices. The computer-readable storage medium is used to store the computer program and other programs and data required by the data processing device, and can also be used to temporarily store data that has been output or will be output.

[0071] The above embodiments are only used to illustrate the design concept and features of the present invention, and their purpose is to enable those skilled in the art to understand the content of the present invention and implement it accordingly. The protection scope of the present invention is not limited to the above embodiments. Therefore, all equivalent changes or modifications made based on the principles and design ideas disclosed in the present invention are within the protection scope of the present invention.

Claims

1. A method for shallow surface seismic P-wave and S-wave nuclear clustering combined tomographic inversion, characterized in that, Includes the following steps: (1) Obtain travel time data of seismic P-waves and S-waves in shallow surface areas; (2) The observation system based on the earthquake P-wave and S-wave travel time observation data calculates the theoretical travel time field of P-wave and S-wave, and establishes the ray path sensitivity matrix based on the theoretical travel time field to form the travel time equation of P-wave and S-wave. (3) Based on the travel time equations of the P-wave and S-wave, a joint inversion objective function of P-wave and S-wave containing a kernel clustering constraint term is constructed, wherein the kernel clustering constraint term is used to constrain the consistency between the P-wave and S-wave velocity structure and the geological unit division. (4) Based on the objective function of the joint inversion of P-wave and S-wave, an augmented equation set is established, and the augmented equation set is solved iteratively. The P-wave and S-wave slowness models are updated iteratively until convergence, and the P-wave velocity model and S-wave velocity model are obtained. (5) Output P-wave velocity model and S-wave velocity model for shallow surface velocity structure inversion and geological structure interpretation.

2. The method according to claim 1, characterized in that, The earthquake P-wave and S-wave travel time observation data were obtained using different observation systems. Specifically, the earthquake P-wave travel time data were obtained by a combination of well-excitation and well-received data and ground-excitation and ground-received data, while the earthquake S-wave travel time data were obtained by ground-excitation and ground-received data.

3. The method according to claim 2, characterized in that, The ratio of ground survey line length to detection depth during the ground excitation and ground reception process is 3:1 to 5:

1.

4. The method according to claim 1, characterized in that, The theoretical travel time fields of the longitudinal and transverse waves were calculated using the fast travel method.

5. The method according to claim 1, characterized in that, The specific construction of the travel time equation is as follows: based on the ray path sensitivity matrices of the P-wave and S-wave, the slowness model increments of the P-wave and S-wave are multiplied by the corresponding sensitivity matrices, so that the results correspond to the travel time residual vectors of the P-wave and S-wave respectively, thereby establishing the linearized travel time equations of the P-wave and S-wave.

6. The method according to claim 1, characterized in that, The objective function for the joint inversion of P-waves and S-waves is: based on the fitting of P-wave and S-wave travel time data, and simultaneously includes model regularization constraints and kernel clustering constraints. The kernel clustering constraints are used to maintain the similarity of P-wave and S-wave velocity structures within the same geological unit and to enhance the separability at the interfaces of different geological units.

7. The method according to claim 6, characterized in that, The regularization constraint term adopts second-order Laplacian regularization to smooth the velocity field and reduce local oscillations.

8. The method according to claim 1, characterized in that, The augmented equations are specifically constructed based on minimizing the objective function of the joint inversion of P-waves and S-waves. The augmented equations include travel time data fitting terms, model regularization terms, and kernel clustering constraint terms.

9. The method according to claim 8, characterized in that, The augmented equations are solved using the LSQR iterative algorithm based on the least squares criterion to obtain the update values ​​for the P-wave and S-wave slowness models.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it implements the method as described in any one of claims 1-9.

Images