Near-surface imaging method and system based on wave equation seismic surface wave full-dispersion spectrum inversion

Abstract

The present invention relates to a near-surface imaging method and system for wave equation seismic surface wave full-frequency dispersion spectrum inversion, comprising: collecting data acquisition information of a target area and inputting the data acquisition information into a grid subdivision model to obtain an observation surface wave ground acquisition system; solving the elastic wave equation to obtain simulated surface wave data in the ground acquisition system; calculating an accompanying earthquake source based on the observed surface wave data and the simulated surface wave data; cross-correlating the forward wavefield of the actual earthquake source with the reverse wavefield of the accompanying earthquake source to calculate the gradient of an updated seismic wave velocity model, and updating a network subdivision model based on the model update gradient; inputting the data to be measured into the updated grid subdivision model to obtain near-surface imaging results from wave equation seismic surface wave full-frequency dispersion spectrum inversion. Because the surface wave dispersion spectrum is used as inversion data, the method avoids manually picking dispersion curves, reduces human errors, and improves the efficiency of large-scale seismic data processing.

Description

Technical Field

[0001] The present invention relates to a near-surface imaging method and system for wave equation seismic surface wave full-frequency dispersion spectrum inversion, belonging to the technical field of land exploration. Background Art

[0002] Seismic surface wave inversion imaging is a key method for understanding the Earth's shallow velocity structure, with widespread applications in engineering geological exploration, petroleum exploration, and crust-mantle scale imaging. Current surface wave imaging methods primarily estimate the corresponding subsurface velocity structure using calculated dispersion curves and their corresponding relationship with velocity models. This method assumes a uniform layered structure in the subsurface. The dispersion curve corresponding to the model is calculated using analytical solutions of seismic surface waves. The optimal one-dimensional velocity model is then derived by comparing it with the dispersion curve extracted from the observational data. Common methods for solving the dispersion equation include the generalized reflection-transmission coefficient method and the fast scalar method. Common methods for extracting surface wave dispersion curves include the dual-station method and multi-channel surface wave analysis. With the widespread adoption of high-density seismic acquisition, multi-channel surface wave analysis has become the mainstream method for surface wave imaging. However, it is worth noting that current surface wave dispersion curve inversion methods assume a laterally uniform subsurface velocity structure and require manual extraction of dispersion curves. Consequently, the inversion yields an average velocity profile beneath the survey line, making them unsuitable for complex structures and large-scale seismic data processing.

[0003] The most advanced full-waveform inversion method for elastic surface waves iteratively updates the velocity model by comparing the predicted seismic waveform with the observed waveform. In theory, this method can achieve high-precision seismic imaging of arbitrarily complex structures. However, in practical applications, due to the complexity of dispersive surface wave waveforms, inversion methods that use waveform subtraction as the objective function are prone to falling into local minima. More advanced objective functions such as adaptive waveform inversion, optimal transport waveform inversion, and travel-time difference waveform inversion emphasize matching the phase information of seismic waves. These methods are more applicable to body waves with simple waveforms, but are difficult to use for surface waves with dispersive characteristics. Due to its disadvantages such as strong nonlinearity and high computational complexity, seismic surface wave full-waveform inversion technology is still in the early stages of research and has not yet been applied to the actual processing of large-scale seismic data. Summary of the Invention

[0004] In response to the above problems, the purpose of the present invention is to provide a near-surface imaging method and system for wave equation seismic surface wave full dispersion spectrum inversion. Since the method uses the surface wave dispersion spectrum as the inversion data, it avoids the manual picking of dispersion curves, reduces human errors, and improves the efficiency of large-scale seismic data processing.

[0005] To achieve the above-mentioned purpose, the present invention proposes the following technical scheme: a wave equation seismic surface wave full-frequency dispersion spectrum inversion near-surface imaging method, comprising the following steps: collecting data acquisition information of a target area, and inputting the data acquisition information into a grid subdivision model to obtain an observation surface wave ground acquisition system; solving the elastic wave wave equation to obtain simulated surface wave data in the ground acquisition system; calculating an accompanying earthquake source based on the observed surface wave data and the simulated surface wave data; cross-correlating the forward wavefield of the actual earthquake source with the reverse wavefield of the accompanying earthquake source to calculate the gradient of the seismic wave velocity model update, and updating the network subdivision model according to the model update gradient; inputting the data to be measured into the updated grid subdivision model to obtain the wave equation seismic surface wave full-frequency dispersion spectrum inversion near-surface imaging results.

[0006] Furthermore, the grid subdivision model is a finite element discrete grid model of comprehensive observation system information and seismic wave velocity established based on the arrangement of earthquake sources and detectors and the surface undulation in the target area.

[0007] Furthermore, the elastic wave equation is:

[0008]

[0009] Among them, Ψ = (v1, v2, v3, σ1, σ2, σ3, σ4, σ5, σ6) is a vector containing three particle velocities and six stresses, T is the transposed matrix, C is the stiffness matrix, E is the spatial difference matrix, f is the source vector, ρ is the density, and I3 is the unit matrix.

[0010] Furthermore, the method for calculating the accompanying earthquake source based on the observed surface wave data and the simulated surface wave data is: through linear Radon transformation, the observed surface wave data and the simulated surface wave data are transformed into the frequency-phase slowness domain, and the optimal transmission distance of the frequency-phase slowness spectrum of the observed surface wave data and the simulated surface wave is calculated, and the accompanying earthquake source is determined according to the optimal transmission distance.

[0011] Furthermore, the observed surface roll data and the simulated surface roll data are transformed into the frequency-phase slowness domain by performing a time Fourier transform on the common shot gathers through a linear Radon transform; and integrating each frequency component in the time Fourier transformed common shot gathers along the spatial two-dimensional phase slowness direction, wherein the integral formula is:

[0012]

[0013] Where D(f,x,y) is the time domain Fourier transform of the original common shot gather, C(f,s x ,s y ) is the calculated surface wave frequency-phase slowness spectrum, f is the frequency, s x is the source position in the x direction, sy is the source position in the y direction, x is the x direction in space, y is the y direction in space, x min is the minimum offset in the x direction, x max is the maximum offset in the x direction, y min is the minimum offset in the y direction, y max is the maximum offset in the y direction.

[0014] Furthermore, the calculation formula of the accompanying earthquake source is:

[0015]

[0016] Where φ(m) is the dispersion spectrum error based on cross-correlation, m is the shear wave velocity model, and d p (m) is the forward seismic wave field, is the real part, iFFT is the inverse Fourier transform; Adj is the adjoint Radon transform; C p (m) is the predicted surface wave dispersion spectrum, C o is the observed surface wave dispersion spectrum.

[0017] Furthermore, the calculation formula of the cross-correlation-based dispersion spectrum error is:

[0018] φ(m)=-<|C P (m,f,s x )|,|C o (f,s x )|>

[0019] Among them, |C P (m,f,s x )| is the amplitude spectrum after normalization of the predicted surface wave frequency-phase slowness spectrum; |C o (f,s x )| is the amplitude spectrum after normalization of the observed surface wave frequency-phase slowness spectrum.

[0020] Furthermore, the method for calculating the gradient of the model iteration and updating the network subdivision model according to the model iteration gradient is as follows: taking the companion source as input, using the spectral element method to back-propagate the companion wavefield; at time zero, cross-correlating the back-propagated wavefield of the companion wavefield with the forward propagated wavefield of the actual source to obtain the gradient of the model iteration; iteratively updating the grid subdivision model by the LBFGS method, and the model update formula is: m = m0-λH -1 g, where λ is the update step size, H is the approximate Hessian matrix, m is the updated seismic wave velocity model; m0 is the initial seismic wave velocity model; g is the update gradient of the seismic wave velocity model; the grid subdivision model is iteratively updated until the convergence condition is met, and the final network subdivision model is output.

[0021] The present invention also discloses a wave equation seismic surface wave full-frequency dispersion spectrum inversion near-surface imaging system, comprising: an actual surface wave observation module, used for collecting data acquisition information of a target area, and inputting the data acquisition information into a grid subdivision model to obtain a seismic surface wave ground observation system; a simulation surface wave module, used for solving the elastic wave wave equation, and obtaining simulated surface wave data in the ground acquisition system; an accompanying source generation module, used for calculating an accompanying source based on the observed surface wave data and the simulated surface wave data; a network subdivision model updating module, used for cross-correlating the forward wave field of the actual source with the reverse wave field of the accompanying source to calculate the gradient of the model iteration, and updating the network subdivision model according to the model iteration gradient; and a result output module, used for inputting the data to be measured into the updated grid subdivision model, and obtaining the wave equation seismic surface wave full-frequency dispersion spectrum inversion near-surface imaging results.

[0022] The present invention also discloses a computer-readable storage medium having a computer program stored thereon. The computer program is executed by a processor to implement any of the above-mentioned wave equation seismic surface wave full-dispersion spectrum inversion near-surface imaging methods.

[0023] The present invention has the following advantages due to the adoption of the above technical solution:

[0024] 1. The present invention uses the surface wave dispersion spectrum as inversion data, thereby avoiding manual picking of dispersion curves, reducing human errors, and improving the efficiency of large-scale seismic data processing.

[0025] 2. Since the present invention uses the spectral element method to simulate the seismic wave field, it is suitable for imaging of any complex medium and avoids the lateral averaging effect on the velocity model of the conventional one-dimensional inversion method.

[0026] 3. The present invention proposes an optimal transport dispersion spectrum objective function, which reduces the dependence of the inversion process on the initial velocity model and improves the stability of the inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 is a flow chart of a near-surface imaging method using full dispersion spectrum inversion of seismic surface waves using a wave equation in one embodiment of the present invention;

[0028] Figure 2 Figure 1 is a representative initial shear wave velocity model and inversion result diagram of an embodiment of the present invention. Figure 2 (a) is a schematic diagram of the initial velocity model. Figure 2 (b) Schematic diagram of the inverted shear wave velocity model;

[0029] Figure 3 is a representative surface wave data waveform and corresponding dispersion spectrum in one embodiment of the present invention. Figure 3 (a) is a surface wave waveform with high complexity; Figure 3 (b) is a relatively simple surface wave dispersion spectrum;

[0030] Figure 4 is a concavity and convexity analysis diagram of different types of objective functions proposed in one embodiment of the present invention;

[0031] Figure 5 This is a diagram showing the sensitivity of surface wave leakage energy and intrinsic frequency dispersion information of different orders to changes in underground longitudinal and shear wave velocities in one embodiment of the present invention; DETAILED DESCRIPTION

[0032] In order to enable those skilled in the art to better understand the technical solutions of the present invention, the present invention will be described in detail through specific embodiments. However, it should be understood that the specific embodiments are provided only for a better understanding of the present invention and should not be construed as limiting the present invention. In the description of the present invention, it should be understood that the terms used are for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0033] In order to solve the problems existing in the prior art, such as the difficulty in processing surface waves with dispersion characteristics, and the strong nonlinearity and poor practicality of seismic surface wave full waveform inversion technology, the present invention discloses a wave equation seismic surface wave full dispersion spectrum inversion near-surface imaging method and system. It combines the advantages of traditional surface wave dispersion curve inversion and the latest elastic wave waveform inversion to propose a wave equation surface wave full dispersion spectrum inversion method. The elastic wave wave equation is solved by numerical methods to obtain predicted surface wave data, and the surface wave full dispersion spectrum calculated from the seismic waveform is compared. The shear wave velocity model is iteratively updated with the help of a unified full waveform inversion framework. Since the surface wave dispersion spectrum is used as the inversion data, the manual picking of the dispersion curve is avoided, human error is reduced, and the efficiency of large-scale seismic data processing is improved. The scheme of the present invention is described in detail below through examples with reference to the accompanying drawings.

[0034] Example 1

[0035] This embodiment discloses a near-surface imaging method using the full dispersion spectrum inversion of seismic surface waves using a wave equation. Figure 1 As shown, the following steps are included:

[0036] S1 collects data information of the target area and inputs the data information into the grid subdivision model to obtain the observation surface wave model.

[0037] The grid model is a finite element discrete grid model that integrates observation system information and seismic wave velocity based on the target area's source and detector layout and surface undulation. Figure 2 As shown, Figure 2(a) shows that the mesh subdivision model in this embodiment is a 1D linear model in which the velocity value increases with depth.

[0038] S2 solves the elastic wave equation by collecting data information and obtains the simulated surface wave model. The simulated surface wave model is as follows: Figure 3 (a)

[0039] The elastic wave equation is:

[0040]

[0041] Among them, Ψ = (v1, v2, v3, σ1, σ2, σ3, σ4, σ5, σ6) is a vector containing three particle velocities and six stresses, T is the transposed matrix, C is the stiffness matrix, E is the spatial difference matrix, f is the source vector, ρ is the density, and I3 is the unit matrix.

[0042] S3 calculates the associated earthquake source based on the observed surface wave data and the simulated surface wave data;

[0043] The method for calculating the associated earthquake source based on the observed surface wave data and the simulated surface wave data is as follows: the observed surface wave model and the simulated surface wave model are transformed into the frequency-phase slowness domain by linear Radon transformation. The surface wave dispersion spectrum after transformation is shown as follows: Figure 3 As shown in (b), the optimal transmission distance of the frequency-phase slowness spectrum of the observed surface wave and the simulated surface wave is calculated, and the accompanying earthquake source is determined based on the optimal transmission distance. Figure 3 As shown, Figure 3 (a) The surface wave waveform, that is, the surface wave waveform without transformation is more complex, while the transformed surface wave waveform is more complex. Figure 3 The surface wave waveform in (b) is relatively simple, but still retains the key dispersion information for velocity inversion.

[0044] The method of transforming the observed surface wave model and the simulated surface wave model into the frequency-phase slowness domain is as follows: the common shot gather is subjected to time Fourier transform through linear Radon transform; each frequency component in the common shot gather after time Fourier transform is integrated along the spatial two-dimensional phase slowness direction, and the integral formula is:

[0045]

[0046] Where D(f,x,y) is the time domain Fourier transform of the original common shot gather, C(f,s x ,s y ) is the calculated surface wave frequency-phase slowness spectrum, f is the frequency, s x is the source position in the x direction, s y is the source position in the y direction, x is the x direction in space, y is the y direction in space, x min is the minimum offset in the x direction, xmax is the maximum offset in the x direction, y min is the minimum offset in the y direction, y max is the maximum offset in the y direction. The surface wave dispersion spectrum obtained by transforming to the frequency-phase slowness domain is as follows: Figure 3 (b) shown.

[0047] The calculation formula of the accompanying earthquake source is:

[0048]

[0049] Where φ(m) is the dispersion spectrum error based on cross-correlation, m is the shear wave velocity model, and d p (m) is the forward seismic wave field, R is the real part, iFFT is the inverse Fourier transform; Adj is the adjoint Radon transform; C p (m) is the predicted surface wave dispersion spectrum, C o is the observed surface wave dispersion spectrum.

[0050] The calculation formula of the dispersion spectrum error based on cross-correlation is:

[0051] φ(m)=-<|C P (m,f,s x )|,|C o (f,sx ) |>

[0052] Among them, |C P (m,f,s x )| is the amplitude spectrum after normalization of the predicted surface wave frequency-phase slowness spectrum; |C o (f,s x )| is the amplitude spectrum after the normalization of the observed surface wave frequency-phase slowness spectrum. Figure 4 As shown, Figure 4 This is a graph analyzing the convexity of different types of objective functions proposed in this embodiment. The dispersion spectrum objective function used in this embodiment has better convexity than the objective function used as comparison data, i.e., the objective function based on waveform comparison, and has lower requirements for the accuracy of the initial velocity model.

[0053] S4 performs cross-correlation on the forward wave field of the actual earthquake source and the reverse wave field of the accompanying earthquake source to calculate the gradient of the model iteration, and updates the network segmentation model according to the model iteration gradient;

[0054] The method for calculating the gradient of model iteration and updating the network subdivision model according to the model iteration gradient is as follows: taking the accompanying source as input, using the spectral element method to back-propagate the accompanying wave field; at time zero, the back-propagated wave field of the accompanying source is cross-correlated with the forward propagated wave field of the actual source to obtain the gradient of model iteration; iteratively updating the grid subdivision model using the LBFGS method, and the model update formula is: m = m0-λH-1 g, where λ is the update step size, H is the approximate Hessian matrix, m is the updated seismic wave velocity model; m0 is the initial seismic wave velocity model; g is the update gradient of the seismic wave velocity model; the grid partitioning model is iteratively updated until the convergence condition is met, and the final network partitioning model is output. The final network partitioning model is as follows Figure 2 (b) Compared with Figure 2 The initial velocity model in (a) is a 1D linear model. Figure 2 The inverted shear wave velocity model in (b) has obvious lateral velocity variations.

[0055] S5 inputs the data to be measured into the updated grid subdivision model to obtain the near-surface imaging results of the wave equation seismic surface wave full-frequency dispersion spectrum inversion.

[0056] like Figure 5 As shown, Figure 5 This chart shows the sensitivity of surface wave leakage energy and eigendispersion of different orders to variations in subsurface P- and S-wave velocities. The fundamental-order eigendispersion is more sensitive to shallow S-wave velocity variations and can be used to constrain S-wave velocity inversion. The fundamental, first-order, and second-order leakage energy dispersions are less sensitive to shallow S-wave velocity and cannot effectively constrain S-wave velocity inversion. The fundamental-order eigendispersion is insensitive to shallow P-wave velocity variations and is generally not used for P-wave velocity inversion. The fundamental, first-order, and second-order leakage energy dispersions are more sensitive to shallow P-wave velocity and can be used to constrain P-wave velocity inversion.

[0057] Example 2

[0058] Based on the same inventive concept, the present invention discloses a near-surface imaging system for wave equation seismic surface wave full-dispersion spectrum inversion, comprising:

[0059] The actual surface wave observation module is used to collect data collection information of the target area and input the data collection information into the grid subdivision model to obtain the surface wave ground observation system;

[0060] The surface wave simulation module is used to solve the elastic wave equation through data acquisition information and obtain the simulated surface wave model;

[0061] An accompanying earthquake source generation module is used to calculate the accompanying earthquake source based on the observed surface wave data and the simulated surface wave data;

[0062] The network segmentation model update module is used to perform cross-correlation on the forward wave field of the actual earthquake source and the reverse wave field of the accompanying earthquake source to calculate the gradient of the model iteration and update the network segmentation model according to the model iteration gradient;

[0063] The result output module is used to input the data to be measured into the updated grid subdivision model to obtain the near-surface imaging results of the wave equation seismic surface wave full-frequency dispersion spectrum inversion.

[0064] Example 3

[0065] Based on the same inventive concept, the present invention discloses a computer-readable storage medium having a computer program stored thereon, which is executed by a processor to implement any of the above-mentioned wave equation seismic surface wave full-frequency dispersion spectrum inversion near-surface imaging methods.

[0066] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0067] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the steps in the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0068] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0069] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 The steps for the function specified in one or more boxes.

[0070] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that the specific embodiments of the present invention can still be modified or replaced by equivalents, and any modifications or equivalent replacements that do not depart from the spirit and scope of the present invention should be included within the scope of protection of the claims of the present invention. The above content is only a specific embodiment of the present application, but the scope of protection of the present application is not limited thereto. Any person skilled in the art who is familiar with the technical field can easily think of changes or replacements within the technical scope disclosed in the present application, which should be included within the scope of protection of the present application. Therefore, the scope of protection of the present application should be based on the scope of protection of the claims.

Claims

1. A near-surface imaging method based on wave equation seismic surface wave full dispersion spectrum inversion, characterized in that: The following steps are involved: Collecting data collection information of a target area, and inputting the data collection information into a grid subdivision model to obtain a ground collection system for observing surface waves; Solving elastic wave equations to obtain simulated surface wave data in the ground acquisition system; calculating an accompanying earthquake source based on the observed surface wave data and the simulated surface wave data; Cross-correlating the forward wavefield of the actual earthquake source with the reverse wavefield of the companion earthquake source to calculate the gradient of the seismic wave velocity model update, and updating the grid subdivision model according to the model update gradient; The data to be measured are input into the updated grid model to obtain the near-surface imaging results of the wave equation seismic surface wave full-frequency dispersion spectrum inversion.

2. The near-surface imaging method of wave equation seismic surface wave full dispersion spectrum inversion according to claim 1, characterized in that: The grid subdivision model is a finite element discrete grid model of comprehensive observation system information and seismic wave velocity established based on the arrangement of earthquake sources and detectors in the target area and the surface undulation.

3. The near-surface imaging method of wave equation seismic surface wave full dispersion spectrum inversion according to claim 1, characterized in that: The elastic wave equation is: Among them, Ψ = (v1, v2, v3, σ1, σ2, σ3, σ4, σ5, σ6) is a vector containing three particle velocities and six stresses, T is the transposed matrix, C is the stiffness matrix, E is the spatial difference matrix, f is the source vector, ρ is the density, and I3 is the unit matrix.

4. The near-surface imaging method of wave equation seismic surface wave full dispersion spectrum inversion according to claim 1, characterized in that: The method for calculating the accompanying earthquake source based on the observed surface wave data and the simulated surface wave data is: The observed surface wave data and the simulated surface wave data are transformed into the frequency-phase slowness domain by linear Radon transformation, and the optimal transmission distance of the frequency-phase slowness spectrum of the observed surface wave data and the simulated surface wave is calculated, and the accompanying earthquake source is determined according to the optimal transmission distance.

5. The near-surface imaging method of wave equation seismic surface wave full dispersion spectrum inversion according to claim 4, characterized in that: The method for transforming the observed surface roll data and the simulated surface roll data into the frequency-phase slowness domain is: Perform time Fourier transform on the common shot gathers through linear Radon transform; Each frequency component in the common shot gather after time Fourier transform is integrated along the two-dimensional phase slowness direction. The integral formula is: Where D(f,x,y) is the time domain Fourier transform of the original common shot gather, C(f,s x ,s y ) is the calculated surface wave frequency-phase slowness spectrum, f is the frequency, s x is the source position in the x direction, s y is the source position in the y direction, x is the x direction in space, y is the y direction in space, x min is the minimum offset in the x direction, x max is the maximum offset in the x direction, y min is the minimum offset in the y direction, y max is the maximum offset in the y direction.

6. The near-surface imaging method of wave equation seismic surface wave full dispersion spectrum inversion according to claim 4, characterized in that: The calculation formula of the accompanying earthquake source is: Where φ(m) is the dispersion spectrum error based on cross-correlation, m is the shear wave velocity model, and d p (m) is the forward seismic wave field, is the real part, iFFT is the inverse Fourier transform; Adj is the adjoint Radon transform; C p (m) is the predicted surface wave dispersion spectrum, C o is the observed surface wave dispersion spectrum.

7. The near-surface imaging method of wave equation seismic surface wave full dispersion spectrum inversion according to claim 6, characterized in that: The calculation formula of the cross-correlation-based dispersion spectrum error is: φ(m)=-<|C P (m,f,s x )|,|C o (f,s x )|> Among them, |C P (m,f,s x )| is the amplitude spectrum after normalization of the predicted surface wave frequency-phase slowness spectrum; |C o (f,s x )| is the amplitude spectrum after normalization of the observed surface wave frequency-phase slowness spectrum.

8. The near-surface imaging method of wave equation seismic surface wave full dispersion spectrum inversion according to claim 1, characterized in that: The method of calculating the gradient of the model iteration and updating the meshing model according to the model iteration gradient is: Taking the companion source as input, backpropagating the companion wavefield using the spectral element method; At time zero, the gradient of the model iteration is obtained by cross-correlating the reverse wavefield of the companion wavefield with the forward wavefield of the actual earthquake source; The mesh model is iteratively updated using the LBFGS method. The model update formula is: m = m0 - λH -1 g, where λ is the update step size, H is the approximate Hessian matrix, m is the updated seismic wave velocity model; m0 is the initial seismic wave velocity model; g is the seismic wave velocity model update gradient; Iteratively update the mesh model until the convergence conditions are met and output the final mesh model.

9. A near-surface imaging system based on wave equation seismic surface wave full dispersion spectrum inversion, characterized in that: include: The actual surface wave observation module is used to collect data acquisition information of the target area and input the data acquisition information into the grid subdivision model to obtain the seismic surface wave ground observation system; A surface wave simulation module, used for solving elastic wave equations and obtaining simulated surface wave data in the ground acquisition system; An accompanying earthquake source generation module is used to calculate an accompanying earthquake source based on the observed surface wave data and the simulated surface wave data; a gridding model updating module, configured to perform cross-correlation on the forward wavefield of the actual earthquake source and the reverse wavefield of the companion earthquake source to calculate a gradient of model iteration, and to update the gridding model according to the gradient of model iteration; The result output module is used to input the data to be measured into the updated grid subdivision model to obtain the near-surface imaging results of the wave equation seismic surface wave full-frequency dispersion spectrum inversion.

10. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a computer program, which is executed by a processor to implement the wave equation seismic surface wave full dispersion spectrum inversion near-surface imaging method according to any one of claims 1 to 8.

Images