Reflection waveform inversion method based on stochastic gradient sampling

Abstract

The invention relates to a reflection waveform inversion method based on stochastic gradient sampling, and the method comprises the steps: S1, obtaining original observation seismic data, carrying out the preprocessing of the original observation seismic data, and obtaining a real seismic record based on the preprocessed original observation seismic data; s2, constructing an initial longitudinal wave velocity model; and S3, a reflection waveform inversion method based on stochastic gradient sampling is adopted, model updating is performed based on real seismic records, and a final longitudinal wave velocity model is obtained. Compared with the prior art, the method provided by the invention does not depend on direct waves and only depend on reflected waves to recover a deep-layer speed structure, the update quantity only comprises a chromatography part of the reflected waves, interference of offset response does not exist, recovery of the depth speed structure is more facilitated, the inversion precision is higher, and the application range is wider.

Description

Technical Field

[0001] The present invention relates to the field of velocity modeling in exploration seismology, and in particular to a reflection waveform inversion method based on random gradient sampling. Background Art

[0002] Reflection waveform inversion (RWI) can use reflected seismic data to obtain underground velocity parameters. It is a seismic waveform inversion technology. The technical process is to excite seismic waves in the field through an artificial source and record the reflected waves reflected back to the surface. Then, an initial underground model is guessed. The reflected waves are simulated based on the guessed model and compared with the observed reflected waves. The difference between the two waveforms can reflect the difference between the guessed initial model and the true underground model. By applying this difference to the guessed initial model, the guessed model can gradually be made to approach the true underground model. However, due to the limitations of acquisition technology, when a given initial model deviates far from the true model, iterative updates cannot gradually bring the initial model closer to the true model. Instead, it will fall into the wrong numerical range. This situation indicates that the inversion has fallen into a local extreme value. The specific reason is that the seismic data currently collected often lack low-frequency components, which can easily lead to periodic jumps when comparing data ("period jumps" means that when the phase difference between the simulated seismic waves corresponding to the initial model and the observed seismic waves corresponding to the true model is greater than 1 / 2 cycle, we will mistakenly identify data from a certain reflection layer as data from another reflection layer. Generally, the lower the frequency of the data, the less likely it is to produce this phenomenon).

[0003] In order to solve the problem of cycle jumping, a multi-scale inversion strategy is proposed, that is, inversion from low-frequency data to high-frequency data. This method can preferentially extract the lower-frequency components in the data and avoid premature recovery of high-frequency information, but it does not solve the problem of missing low-frequency data in principle.

[0004] Based on this, Chinese patent application CN115657131A provides an elastic wave full waveform inversion method based on random gradient sampling. This method uses random gradient sampling to perform elastic wave full waveform inversion, updates the elastic parameter model based on real frequency domain seismic records, obtains the final elastic parameter model, and implements full waveform inversion. Although it can reduce the impact of periodic skipping, the direct wave component plays a major role in this method. However, in actual exploration, the direct wave penetration depth is too shallow, limiting the accuracy of inversion of deep velocity structure. To recover the deep velocity structure, it is necessary to use the reflection waves in the data. When using the reflection waves in the data for inversion, the model update generated by this method consists of two parts: the reflection wave tomography component and the reflection wave migration response component. The reflection wave tomography component mainly recovers the model background component, while the migration response component mainly recovers the high-frequency components in the model. The latter has an energy that is an order of magnitude stronger than the former, which can easily cause the inversion to fall into outliers prematurely.

[0005] Therefore, providing a method that can accurately invert deep velocity structure is a technical problem that needs to be solved. Summary of the Invention

[0006] The purpose of the present invention is to overcome the defects of the above-mentioned prior art and to provide a reflection waveform inversion method based on random gradient sampling. The method does not rely on direct waves and can accurately restore the deep velocity structure during earthquakes by inverting the reflected waves.

[0007] The purpose of the present invention can be achieved by the following technical solutions:

[0008] The present invention provides a reflection waveform inversion method based on random gradient sampling, which is characterized by comprising:

[0009] S1. Obtaining original observed seismic data, preprocessing the original observed seismic data, and obtaining real earthquake records based on the preprocessed original observed seismic data;

[0010] S2, constructing the initial longitudinal wave velocity model;

[0011] S3. Using a reflection waveform inversion method based on random gradient sampling, the model is updated based on the real seismic records to obtain a final P-wave velocity model.

[0012] As a preferred technical solution, the method for obtaining the final longitudinal wave velocity model is performed by iteratively performing the following steps:

[0013] S31. Perform forward simulation based on the P-wave velocity model to obtain simulated earthquake records and the forward wave field {u n(x, z, t)}, (n = 1, 2, 3 ... Ns), where g represents the spatial position of the detector, t represents the time of wave field propagation, (x, z) represents the spatial position coordinates, and Ns represents the number of shot points; and in the first iteration, the P-wave velocity model is the initial P-wave velocity model; in the kth iteration, the P-wave velocity model is the P-wave velocity model of the k-1th inversion;

[0014] S32, calculating the residual seismic record of the real seismic record and the simulated seismic record of each shot, and calculating the sum of the residual seismic records. If the residual seismic record is less than a preset value, outputting the P-wave velocity model of the current iteration as the final P-wave velocity model; otherwise, executing S33;

[0015] S33, using each residual seismic record as a companion source, performing back propagation of the detection wave field to obtain a back propagation seismic record

[0016] S34, calculating the updating direction of the reverse propagation seismic record and wavelength using a gradient calculation method based on random spatial movement;

[0017] S35 , calculating the step length based on the update direction, updating the longitudinal wave velocity model based on the step length, and then returning to S31 .

[0018] As a preferred technical solution, the forward modeling method is: using the isotropic acoustic wave equation to perform forward modeling based on the longitudinal wave velocity model, the expression is:

[0019]

[0020] Where v(x) represents the longitudinal wave velocity at the spatial coordinate x; u n (x,t;x s ) represents the location at the earthquake source x s The displacement of the particle at time t; represents the Laplace operator; f(t) represents the earthquake source at time t.

[0021] As a preferred technical solution, the method for calculating the residual seismic record is: in, represents the simulated seismic record of the nth shot at spatial position g; Represents the real earthquake record of the nth shot at spatial position g.

[0022] As a preferred technical solution, the method of detecting wavefield back propagation is: using the fourth-order spatial and second-order temporal finite difference methods, implementing the PML absorbing boundary condition at the boundary, using the residual seismic record of each shot as the accompanying source, performing detecting wavefield back propagation, and satisfying:

[0023]

[0024] Where v(x) represents the longitudinal wave velocity at the spatial coordinate x; u n (x,t;x s ) represents the location at the earthquake source x s The displacement of the particle at time t; Represents the Laplace operator; res represents the associated earthquake source, that is, the residual earthquake record.

[0025] As a preferred technical solution, the S34 includes:

[0026] For each shot, the corresponding forward wave field u is simulated n (x,z,t) and back-propagation seismic records

[0027] The range of random spatial movement is defined as: h∈[-k1λ,k1λ], where k1 represents the proportional coefficient; λ represents the longitudinal wave wavelength and v represents the longitudinal wave velocity, f0 is the main frequency of the wavelet;

[0028] Randomly moving the forward wave field and the reverse seismic record based on the spatial movement range;

[0029] The updated direction of the P-wave velocity of each shot is calculated based on the forward wave field and the reverse seismic record after random spatial shift, and the updated direction matrix is generated.

[0030] As a preferred technical solution, the random movement method is:

[0031] u n (x,z,t)|h=u n (x+h,z,t),

[0032]

[0033] Among them, u n (x,z,t) represents the forward wave field; h represents the random spatial movement range; Represents reverse earthquake record.

[0034] As a preferred technical solution, the method for calculating the update direction is:

[0035]

[0036] in, Indicates the update direction of the nth shot; v indicates the longitudinal wave velocity; u n (x,z,t)|h represents the forward wave field after random movement; represents the back propagation seismic record after random movement, and (x,z) represents the spatial position coordinates.

[0037] As a preferred technical solution, the method for calculating the step length includes:

[0038]

[0039] Among them, λ k represents the step size of the kth iteration; ε represents the proportional coefficient used to update the model; v represents the longitudinal wave velocity; g RSS (x,z) represents the update direction matrix.

[0040] As a preferred technical solution, the method for updating the longitudinal wave velocity model is:

[0041] v k+1 =v k +λ k g RSS (x,z),

[0042] Among them, v k+1 represents the k+1th longitudinal wave velocity model; v k represents the k-th longitudinal wave velocity model; λ k represents the step length; g RSS (x,z) represents the update direction matrix.

[0043] Compared with the prior art, the present invention has the following beneficial effects:

[0044] 1) The method provided by the present invention does not rely on direct waves, but only on reflected waves to restore the deep velocity structure during earthquakes. The update only includes the reflected wave tomography part and does not synchronously update the reflected wave offset response component. Therefore, the problem of the inversion falling into outliers prematurely due to the update of the reflected wave offset response component is avoided, which helps to recover the background components of the model and has strong practicality.

[0045] 2) The present invention utilizes random spatial shift approximation of a single reference model wavefield to realize the calculation of several random sampling model wavefields, without the need to additionally calculate the gradients of multiple models to expand the search space, and without the need to additionally store the gradients of multiple random sampling models.

[0046] 3) The present invention can use a poor initial velocity model or seismic data lacking low-frequency components to perform reflection waveform inversion, thereby reducing the adverse effects of the cycle skipping phenomenon on the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 is a flow chart of the method of the present invention;

[0048] Figure 2This is a real earthquake record of Example 2 of the present invention;

[0049] Figure 3 Schematic diagram of the first generation update direction (negative gradient) of conventional reflection waveform inversion in Example 2 of the present invention;

[0050] Figure 4 This is a schematic diagram of conventional reflection waveform inversion results in Example 2 of the present invention;

[0051] Figure 5 Schematic diagram of the first generation update direction (negative gradient) of the stochastic gradient reflection waveform inversion according to embodiment 2 of the present invention;

[0052] Figure 6 This is a schematic diagram of the random gradient reflection waveform inversion results of Example 2 of the present invention;

[0053] Figure 7 This is a schematic diagram comparing the results of random gradient reflection waveform inversion + conventional reflection waveform inversion in Example 2 of the present invention. DETAILED DESCRIPTION

[0054] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0055] Unless otherwise defined, the technical or scientific terms used in this application should have the ordinary meaning understood by a person of ordinary skill in the technical field to which this application belongs. The words "one", "a", "the" and the like used in this application do not indicate a limit on quantity and may indicate the singular or plural. The terms "include", "comprise", "have" and any variations thereof used in this application are intended to cover non-exclusive inclusions; for example, a process, method, system, product or device that includes a series of steps or modules (units) is not limited to the listed steps or units, but may also include steps or units that are not listed, or may also include other steps or units that are inherent to these processes, methods, products or devices. The words "connect", "connected", "coupled" and the like used in this application are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. The word "multiple" used in this application refers to two or more. "And / or" describes the association relationship of associated objects, indicating that three relationships can exist. For example, "A and / or B" can mean: A exists alone, A and B exist at the same time, and B exists alone. The character " / " generally indicates that the objects before and after are in an "or" relationship. The terms "first", "second", "third", etc. involved in this application are only used to distinguish similar objects and do not represent a specific order for the objects.

[0056] Example 1

[0057] In some underground structures that are complex and have large velocity variations, especially when strong low-speed anomalies occur, complex wave phenomena will occur. In actual exploration areas, on the one hand, due to the absorption and attenuation of low-frequency seismic data, and on the other hand due to technical limitations, the detector cannot receive low-frequency signals. The lack of low-frequency seismic data can easily cause periodic jumps in nonlinear inversion and fall into local extremes, seriously affecting the accuracy of the inversion results. Conventional reflection waveform inversion is difficult to adapt to such problems. Therefore, overcoming the periodic jump phenomenon and avoiding falling into local extremes in reflection waveform inversion is an urgent problem to be solved. In order to solve the problems existing in the prior art, the present invention proposes a reflection waveform inversion method based on random gradient sampling. In the framework of the traditional reflection waveform inversion program, it is only necessary to increase the random spatial movement of the wave field, and use the wave field calculation after random spatial movement to update the direction. Other contents remain basically unchanged, and the overall process is simple. The method flow is as follows. Figure 1 shown.

[0058] S1. Obtain original observed seismic data, preprocess the original observed seismic data, and obtain real earthquake records based on the preprocessed original observed seismic data.

[0059] In detail, the steps of preprocessing the original observed seismic data in the present invention include denoising, filtering, etc.

[0060] S2. Construct an initial longitudinal wave velocity model.

[0061] S3. Using the reflection waveform inversion method based on random gradient sampling, the model is updated based on real seismic records to obtain the final P-wave velocity model.

[0062] S31. Perform forward simulation based on the P-wave velocity model to obtain simulated earthquake records and the forward wave field {u n (x, z, t)}, (n = 1, 2, 3…Ns), where g represents the spatial position of the detector, t represents the time of wave field propagation, (x, z) represents the spatial position coordinates, and Ns represents the number of shot points. In the first iteration, the P-wave velocity model is the initial P-wave velocity model; in the kth iteration, the P-wave velocity model is the P-wave velocity model of the k-1th inversion.

[0063] Specifically, the forward modeling is performed based on the longitudinal wave velocity model using the isotropic acoustic wave equation, which is expressed as follows:

[0064]

[0065] Where v(x) represents the longitudinal wave velocity at the spatial coordinate x; u n (x,t;x s ) represents the location at the earthquake source x s The displacement of the particle at time t; represents the Laplace operator; f(t) represents the earthquake source at time t.

[0066] S32. Calculate the residual seismic record of the real seismic record and the simulated seismic record of each shot, and calculate the sum of the residual seismic records. If the residual seismic record is less than a preset value, output the P-wave velocity model of the current iteration as the final P-wave velocity model; otherwise, execute S33.

[0067] The specific expression for calculating the residual is: in, represents the simulated seismic record of the nth shot at spatial position g; Represents the real earthquake record of the nth shot at spatial position g.

[0068] S33, taking each residual seismic record as an accompanying source, performing the detection end wave field back propagation, and obtaining the back propagation seismic record

[0069] Because this process is non-causal, its wavefield propagates from the future to the past, so the propagation control equations used are consistent with the forward propagation equations. Specifically, using the fourth-order spatial and second-order temporal finite difference methods, PML absorbing boundary conditions are implemented at the boundaries, and the residual seismic records of each shot are used as the accompanying source to perform the detection wavefield back propagation, satisfying:

[0070]

[0071] Where v(x) represents the longitudinal wave velocity at the spatial coordinate x; u n (x,t;x s ) represents the location at the earthquake source x s The displacement of the particle at time t; Represents the Laplace operator; res represents the associated earthquake source, that is, the residual earthquake record.

[0072] S34. Use a gradient calculation method based on random spatial movement to calculate the update direction of the back-propagation seismic record and wavelength.

[0073] S341, simulate the corresponding forward wave field u for each shot n (x,z,t) and back-propagation seismic records

[0074] S342. Define the random spatial movement range as: h∈[-k1λ,k1λ], where k1 represents the proportional coefficient; λ represents the longitudinal wave wavelength and v represents the longitudinal wave velocity, and f0 is the main frequency of the wavelet.

[0075] S343. Randomly move the forward wave field and the reverse seismic record based on the spatial movement range.

[0076] Unlike seismic model inversion based on direct waves, which spatially shifts the forward and reverse wavefields in the same direction when calculating model updates, the present invention relies solely on reflected waves for model updates, requiring the forward and reverse wavefields to be spatially shifted in opposite directions when calculating model updates. This difference arises from the fact that, when modeling based on direct waves, the forward-scattered normal direct wave and the reverse direct wave, both of which are scattered forward in the same direction, must be cross-correlated. However, the forward reflected wave is backscattered, while the reverse wavefield is forward-scattered. Therefore, the update calculation differs from that for direct waves, requiring the forward and reverse wavefields to shift in opposite directions.

[0077] The specific update method is:

[0078] u n (x,z,t)|h=u n (x+h,z,t),

[0079]

[0080] Among them, u n (x,z,t) represents the forward wave field; h represents the random spatial movement range; Represents reverse earthquake record.

[0081] S344. Calculate the updated direction of the longitudinal wave velocity of each shot based on the forward wave field and the reverse seismic record after random spatial movement, and collectively generate an updated direction matrix.

[0082] Among them, the method for calculating the update direction of the nth shot is:

[0083]

[0084] in, Indicates the update direction of the nth shot; v indicates the longitudinal wave velocity; u n (x,z,t)|h represents the forward wave field after random movement; represents the back propagation seismic record after random movement, and (x,z) represents the spatial position coordinates.

[0085] S35. Calculate the step size based on the updated direction, update the longitudinal wave velocity model based on the step size, and then return to S31.

[0086] S351, calculation step size:

[0087]

[0088] Among them, λ k represents the step size of the kth iteration; ε represents the proportional coefficient used to update the model; v represents the longitudinal wave velocity; g RSS (x,z) represents the update direction matrix.

[0089] S352. Update the longitudinal wave velocity model:

[0090] v k+1 =v k +λ k g RSS (x,z),

[0091] Among them, v k+1 represents the k+1th longitudinal wave velocity model; v k represents the k-th longitudinal wave velocity model; λ k represents the step length; g RSS (x,z) represents the update direction matrix.

[0092] Example 2

[0093] In this embodiment, the model of the two-dimensional Gaussian ball low-speed anomaly is used as the real model, such as Figure 2As shown in the figure, the density in the model is 2000 kg / m3 and remains constant. It has a total of 401×201 grids with a grid spacing of 10 m×10 m. A two-dimensional isotropic medium acoustic wave forward modeling was performed on this model. A total of 40 shots were simulated, with the shots evenly distributed on the surface and a shot spacing of 100 m. Receivers were evenly distributed on the surface at 10-m intervals. The excitation source was an 18 Hz Ricker wavelet. The seismic record was received for a 2.0-second period with a 1-millisecond interval. A uniform velocity model was used as the initial model for the inversion.

[0094] In order to verify the accuracy and feasibility of the method provided by the present invention, based on the initial model provided above, this embodiment uses the random gradient-based reflection waveform inversion (RSS-RWI) method provided by the present invention and the conventional reflection waveform inversion (CRWI) method to conduct comparative experiments.

[0095] In detail, CRWI and RSS-RWI are used to calculate the update direction (negative gradient) of the first generation model, where the results of CRWI are as follows: Figure 3 and Figure 4 As shown, from Figure 1 The real model provided shows that it should be updated in the negative direction of velocity near the center of velocity anomaly. Figure 3 It can be seen that the update direction of CRWI at the corresponding position is positive, which is not in line with expectations, indicating that the inversion is moving in the wrong direction and a cycle jump phenomenon occurs. Figure 4 It can be seen that the CRWI inversion results have periodic jumps and fall into local extreme values. The inversion results have strong velocity discontinuities, which are far from the true model, making subsequent seismic interpretation difficult.

[0096] Figure 5 and Figure 6 To use the results of RSS-RWI, Figure 5 It can be seen that the update direction of RSS-RWI at the corresponding position is negative, which is in line with expectations, making the inverted model closer to the true model; Figure 6 It can be seen that the results are more consistent with the background characteristics of the real model and can provide a better initial model for reflection waveform inversion.

[0097] In summary, it can be seen that the method provided by the present invention can invert the reflected waveform more accurately than the conventional inversion method.

[0098] In order to further verify that the method provided by the present invention still has advantages in the face of poor initial model and low-frequency data missing, the reflection wave inversion result after performing RSS-RWI is used as the new initial model to perform CRWI for accuracy verification. Figure 7 As shown, it can be seen that Figure 1The model in is very close, and there is no cycle jump. Instead, a high-precision and high-resolution inversion model is obtained, that is, the method provided by the present invention has superiority.

[0099] In addition, the present invention also provides a reflection waveform inversion device based on random gradient sampling, including a central processing unit (CPU), which can perform various appropriate actions and processes according to computer program instructions stored in a read-only memory (ROM) or loaded from a storage unit into a random access memory (RAM). The RAM can also store various programs and data required for device operation. The CPU, ROM, and RAM are connected to each other via a bus. An input / output (I / O) interface is also connected to the bus.

[0100] Many components in a device are connected to the I / O interface, including: input units, such as a keyboard and mouse; output units, such as various types of displays and speakers; storage units, such as magnetic disks and optical disks; and communication units, such as network cards, modems, and wireless communication transceivers. The communication unit allows the device to exchange information / data with other devices via computer networks such as the Internet and / or various telecommunication networks.

[0101] The processing unit performs the various methods and processes described above, such as methods S1 to S3. For example, in some embodiments, methods S1 to S3 may be implemented as a computer software program, which is tangibly contained in a machine-readable medium, such as a storage unit. In some embodiments, part or all of the computer program may be loaded and / or installed on the device via a ROM and / or a communication unit. When the computer program is loaded into the RAM and executed by the CPU, one or more steps of methods S1 to S3 described above may be performed. Alternatively, in other embodiments, the CPU may be configured to execute methods S1 to S3 by any other appropriate means (e.g., by means of firmware).

[0102] The functions described above herein may be performed, at least in part, by one or more hardware logic components. For example, and without limitation, exemplary types of hardware logic components that may be used include: field programmable gate arrays (FPGAs), application specific integrated circuits (ASICs), application specific standard products (ASSPs), systems on chip (SOCs), complex programmable logic devices (CPLDs), and the like.

[0103] The program code for implementing the method of the present invention can be written in any combination of one or more programming languages. Such program code can be provided to a processor or controller of a general-purpose computer, a special-purpose computer, or other programmable data processing device so that when the program code is executed by the processor or controller, the functions / operations specified in the flow chart and / or block diagram are implemented. The program code can be executed entirely on the machine, partially on the machine, as a stand-alone software package, partially on the machine and partially on a remote machine, or entirely on a remote machine or server.

[0104] In the context of the present invention, machine-readable medium can be a tangible medium that can contain or store a program for use with an instruction execution system, device or equipment or used in combination with an instruction execution system, device or equipment. Machine-readable medium can be a machine-readable signal medium or a machine-readable storage medium. Machine-readable medium can include, but is not limited to, electronic, magnetic, optical, electromagnetic, infrared or semiconductor systems, devices or equipment, or any suitable combination of the foregoing. More specific examples of machine-readable storage media can include electrical connections based on one or more lines, portable computer disks, hard disks, random access memories (RAM), read-only memories (ROM), erasable programmable read-only memories (EPROM or flash memory), optical fibers, portable compact disk read-only memories (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination of the foregoing.

[0105] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in the present invention, and such modifications or substitutions are intended to be within the scope of protection of the present invention. Therefore, the scope of protection of the present invention shall be subject to the scope of protection of the claims.

Claims

1. A reflection waveform inversion method based on random gradient sampling, characterized in that: include: S1. Obtaining original observed seismic data, preprocessing the original observed seismic data, and obtaining real earthquake records based on the preprocessed original observed seismic data; S2, constructing the initial longitudinal wave velocity model; S3. Using a reflection waveform inversion method based on random gradient sampling, the model is updated based on the real seismic records to obtain a final P-wave velocity model.

2. The reflection waveform inversion method based on random gradient sampling according to claim 1 is characterized in that: The final P-wave velocity model is obtained by iteratively performing the following steps: S31. Perform forward simulation based on the P-wave velocity model to obtain simulated earthquake records and the forward wave field {u n (x, z, t)}, (n = 1, 2, 3 ... Ns), where g represents the spatial position of the detector, t represents the time of wave field propagation, (x, z) represents the spatial position coordinates, and Ns represents the number of shot points; and in the first iteration, the P-wave velocity model is the initial P-wave velocity model; in the kth iteration, the P-wave velocity model is the P-wave velocity model of the k-1th inversion; S32, calculating the residual seismic record of the real seismic record and the simulated seismic record of each shot, and calculating the sum of the residual seismic records. If the residual seismic record is less than a preset value, outputting the P-wave velocity model of the current iteration as the final P-wave velocity model; otherwise, executing S33; S33, using each residual seismic record as a companion source, performing back propagation of the detection wave field to obtain a back propagation seismic record (n=1,2,3…Ns); S34, calculating the updating direction of the reverse propagation seismic record and wavelength using a gradient calculation method based on random spatial movement; S35 , calculating the step length based on the update direction, updating the longitudinal wave velocity model based on the step length, and then returning to S31 .

3. The reflection waveform inversion method based on random gradient sampling according to claim 2 is characterized in that: The forward modeling method is to use the isotropic acoustic wave equation to perform forward modeling based on the longitudinal wave velocity model, and the expression is: Where v(x) represents the longitudinal wave velocity at the spatial coordinate x; u n (x,t;x s ) represents the location at the earthquake source x s The displacement of the particle at time t; represents the Laplace operator; f(t) represents the earthquake source at time t.

4. The reflection waveform inversion method based on random gradient sampling according to claim 2, characterized in that: The method for calculating the residual seismic record is: in, represents the simulated seismic record of the nth shot at spatial position g; Represents the real earthquake record of the nth shot at spatial position g.

5. The reflection waveform inversion method based on random gradient sampling according to claim 2, characterized in that: The method of detecting wavefield back propagation is to use the fourth-order spatial and second-order temporal finite difference methods, implement the PML absorbing boundary condition at the boundary, use the residual seismic record of each shot as the accompanying source, and perform detecting wavefield back propagation to meet the following requirements: Where v(x) represents the longitudinal wave velocity at the spatial coordinate x; u n (x,t;x s ) represents the location at the earthquake source x s The displacement of the particle at time t; Represents the Laplace operator; res represents the associated earthquake source, that is, the residual earthquake record.

6. The reflection waveform inversion method based on random gradient sampling according to claim 2, characterized in that: The S34 includes: For each shot, the corresponding forward wave field u is simulated n (x,z,t) and back-propagation seismic records The range of random spatial movement is defined as: h∈[-k1λ,k1λ], where k1 represents the proportional coefficient; λ represents the longitudinal wave wavelength and v represents the longitudinal wave velocity, f0 is the main frequency of the wavelet; Randomly moving the forward wave field and the reverse seismic record based on the spatial movement range; The updated direction of the P-wave velocity of each shot is calculated based on the forward wave field and the reverse seismic record after random spatial shift, and the updated direction matrix is generated.

7. The reflection waveform inversion method based on random gradient sampling according to claim 6, characterized in that: The random movement method is: u n (x,z,t)|h=u n (x+h,z,t), Among them, u n (x,z,t) represents the forward wave field; h represents the random spatial movement range; Represents reverse earthquake record.

8. The reflection waveform inversion method based on random gradient sampling according to claim 6, characterized in that: The method for calculating the update direction is: in, Indicates the update direction of the nth shot; v indicates the longitudinal wave velocity; u n (x,z,t)|h represents the forward wave field after random movement; represents the back propagation seismic record after random movement, and (x,z) represents the spatial position coordinates.

9. The method for reflection waveform inversion based on random gradient sampling according to claim 2, characterized in that: The method for calculating the step size includes: Among them, λ k represents the step size of the kth iteration; ε represents the proportional coefficient used to update the model; v represents the longitudinal wave velocity; g RSS (x,z) represents the update direction matrix.

10. The reflection waveform inversion method based on random gradient sampling according to claim 2, characterized in that: The method for updating the P-wave velocity model is: v k+1 =v k +λ k g RSS (x,z), Among them, v k+1 represents the k+1th longitudinal wave velocity model; v k represents the k-th longitudinal wave velocity model; λ k represents the step length; g RSS (x,z) represents the update direction matrix.

Images