A reflected waveform inversion method based on random gradient sampling

Abstract

The application relates to a reflection waveform inversion method based on random gradient sampling, comprising the following steps: S1, acquiring original observed seismic data, and preprocessing the original observed seismic data, and acquiring real seismic records based on the original observed seismic data after preprocessing; S2, constructing an initial P-wave velocity model; S3, adopting the reflection waveform inversion method based on random gradient sampling, and updating the model based on the real seismic records to acquire a final P-wave velocity model. Compared with the prior art, the application provides a method for restoring deep velocity structure only by relying on reflected waves instead of direct waves, and the updating amount only contains a reflected wave tomography part without the interference of a migration response, which is more helpful for restoring the deep velocity structure, the inversion precision is higher, and the practical range is wider.

Description

Technical Field

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

[0002] Reflection waveform inversion is a seismic waveform inversion technique that uses reflected seismic data to obtain subsurface velocity parameters. The process involves generating seismic waves in the field using an artificial source and recording the reflected waves back to the surface. An initial subsurface model is then hypothesized, and the reflected waves are simulated based on this model. These simulated waves are then compared with the observed reflected waves. The difference between the two waveforms reflects the difference between the hypothesized initial model and the actual subsurface model. Applying this difference to the hypothesized initial model gradually makes the hypothesized model approximate the actual subsurface model. However, due to limitations in acquisition technology, when the given initial model deviates significantly from the actual model, iterative updates cannot gradually bring the initial model closer to the actual model. Instead, it will fall into an erroneous numerical range. This situation indicates that the inversion has fallen into a local extremum. The specific reason is that the currently acquired seismic data often lacks low-frequency components, which can easily lead to periodic jumps when comparing data. ("Periodic jumps" refers to the phenomenon where, when the phase difference between the simulated seismic wave corresponding to the initial model and the observed seismic wave corresponding to the actual model is greater than 1 / 2 period, we will mistakenly identify data from one reflection layer as data from another reflection layer. Generally, the lower the frequency of the data, the less likely this phenomenon is to occur.)

[0003] To address the periodic jump problem, a multi-scale inversion strategy is proposed, which involves inverting from low-frequency data to high-frequency data. This method can prioritize the extraction of lower-frequency components in the data and avoid premature recovery of high-frequency information, but it does not, in principle, solve the problem of missing low-frequency data.

[0004] Based on this, Chinese patent application CN115657131A provides a method for full waveform inversion of elastic waves based on stochastic gradient sampling. This method uses stochastic gradient sampling to perform full waveform inversion of elastic waves, updating the elastic parameter model based on real seismic records in the frequency domain to obtain the final elastic parameter model, thus achieving full waveform inversion. Although it can reduce the influence of periodic jump phenomena, the direct wave component is the main active component in this method. However, in actual exploration, the penetration depth of direct waves is too shallow, limiting the accuracy of inversion for deep velocity structures. To recover the deep velocity structure, reflected waves from the data must be utilized. When using reflected waves for inversion, the model update generated by this method includes two parts: the reflected wave tomographic component and the reflected wave migration response component. The reflected wave tomographic component mainly recovers the background components of the model, while the migration response component mainly recovers the high-frequency components in the model. The latter has an energy level an order of magnitude higher than the former, which can easily lead to the inversion prematurely falling into outliers.

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

[0006] The purpose of this invention is to overcome the shortcomings of the existing technology and provide a reflection waveform inversion method based on stochastic gradient sampling. The method does not rely on direct waves and achieves accurate recovery of depth, velocity and structure during earthquakes by inverting the reflected waves.

[0007] The objective of this invention can be achieved through the following technical solutions:

[0008] This invention provides a reflection waveform inversion method based on stochastic gradient sampling, characterized by comprising:

[0009] S1. Obtain raw observed seismic data, preprocess the raw observed seismic data, and obtain real seismic records based on the preprocessed raw observed seismic data.

[0010] S2. Construct the initial P-wave velocity model;

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

[0012] As a preferred technical solution, the method for obtaining the final P-wave velocity model by iteratively executing the following steps is as follows:

[0013] S31. Perform forward modeling based on the P-wave velocity model to obtain simulated seismic records. and forward propagation wave field {u n(x,z,t)}, (n=1,2,3…Ns), where g represents the spatial position of the detector, t represents the wave field propagation time, (x,z) represents the spatial 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 k-th iteration, the P-wave velocity model is the P-wave velocity model inverted in the (k-1)-th iteration;

[0014] S32. Calculate the residual seismic record of the real seismic record and the simulated seismic record for 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.

[0015] S33. Using each of the aforementioned residual seismic records as an adjoint source, perform end-field backpropagation to obtain the backpropagation seismic record.

[0016] S34. The update direction of the back-transmitted seismic record and wavelength is calculated using a gradient calculation method based on random spatial movement.

[0017] S35. Calculate the step size based on the update direction, update the P-wave velocity model based on the step size, and then return to S31.

[0018] As a preferred technical solution, the forward modeling method is as follows: Forward modeling is performed based on a longitudinal wave velocity model using the isotropic acoustic wave equation, the expression of which is:

[0019]

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

[0021] As a preferred technical solution, the method for calculating the residual seismic record is as follows: in, This represents the simulated seismic record of the nth shot at spatial location g; This represents the actual seismic record of the nth shot at spatial location g.

[0022] As a preferred technical solution, the method for end-detector wavefield backpropagation is as follows: using the fourth-order spatial and second-order temporal finite difference method, a PML absorbing boundary condition is implemented at the boundary, and the residual seismic record of each shot is used as the adjoint source to perform end-detector wavefield backpropagation, satisfying:

[0023]

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

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

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

[0027] Define the range of random spatial movement as: h∈[-k1λ,k1λ], where k1 represents the scaling factor; λ represents the longitudinal wave wavelength and v represents the longitudinal wave velocity, and f0 is the dominant frequency of the wavelet;

[0028] The forward propagation wavefield and backward propagation seismic records are randomly moved based on the aforementioned spatial movement range;

[0029] The P-wave velocity update direction for each shot is calculated based on the forward propagation wavefield and the reverse propagation seismic records after random spatial movement, and the updated direction matrix is ​​generated by combining them.

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

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

[0032]

[0033] Among them, u n (x,z,t) represents the propagating wave field; h represents the range of random spatial movement. This indicates the reverse transmission of earthquake records.

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

[0035]

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

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

[0038]

[0039] Where, λ k The step size for the k-th iteration is represented by ε; the scaling factor used to update the model is represented by v; and the P-wave velocity is represented by g. RSS (x,z) represents updating the direction matrix.

[0040] As a preferred technical solution, the method for updating the P-wave velocity model is as follows:

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

[0042] Among them, v k+1 This represents the P-wave velocity model for the (k+1)th time; v k Denotes the P-wave velocity model for the k-th time; λ k Indicates step size; g RSS (x,z) represents updating the direction matrix.

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

[0044] 1) The method provided by this invention does not rely on direct waves, but only on reflected waves to recover the deep velocity structure during an earthquake. Moreover, the update only includes the reflected wave tomography part and does not update the reflected wave migration response component synchronously. Therefore, it avoids the problem of the inversion falling into outliers too early due to the update of the reflected wave migration response component, which helps to recover the model background components and has strong practicality.

[0045] 2) This invention uses the random spatial movement of a single reference model wavefield to approximate the calculation of several random sampled model wavefields, without the need to calculate the gradients of multiple models to expand the search space, nor the need to store the gradients of multiple random sampled models.

[0046] 3) This invention can use poor initial velocity models or seismic data lacking low-frequency components to perform reflection waveform inversion, thereby reducing the adverse effects of period jump phenomenon on the inversion results. Attached Figure Description

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

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

[0049] Figure 3 This is a schematic diagram of the first-generation update direction (negative gradient) for conventional reflection waveform inversion in Embodiment 2 of the present invention;

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

[0051] Figure 5 This is a schematic diagram of the first-generation update direction (negative gradient) for stochastic gradient reflection waveform inversion in Embodiment 2 of the present invention;

[0052] Figure 6 This is a schematic diagram of the stochastic gradient reflection waveform inversion result in Embodiment 2 of the present invention;

[0053] Figure 7 This is a schematic diagram comparing the results of stochastic gradient reflection waveform inversion and conventional reflection waveform inversion in Embodiment 2 of the present invention. Detailed Implementation

[0054] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0055] Unless otherwise defined, the technical or scientific terms used in this application shall have the ordinary meaning understood by one of ordinary skill in the art to which this application pertains. The terms “a,” “an,” “an,” “the,” and similar words used in this application do not indicate quantity limitation and may indicate singular or plural. The terms “comprising,” “including,” “having,” and any variations thereof used in this application are intended to cover non-exclusive inclusion; 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 not listed, or may include other steps or units inherent to these processes, methods, products, or devices. The terms “connected,” “linked,” “coupled,” and similar words used in this application are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. “Multiple” used in this application refers to two or more. “And / or” describes the relationship between related objects, indicating that three relationships may exist; for example, “A and / or B” can represent: A alone, A and B simultaneously, and B alone. The character " / " generally indicates that the preceding and following objects are in an "or" relationship. The terms "first," "second," and "third" used in this application are merely to distinguish similar objects and do not represent a specific ordering of the objects.

[0056] Example 1

[0057] In some complex underground structures with large velocity variations, especially with strong low-velocity anomalies, complex wave phenomena can occur. In actual exploration areas, due to the absorption and attenuation of low-frequency seismic data, and the limitation of technology, detectors cannot receive low-frequency signals. The lack of low-frequency seismic data easily leads to periodic jumps in nonlinear inversion, causing it to fall into local extrema and severely affecting the accuracy of the inversion results. Conventional reflection waveform inversion is difficult to adapt to these problems. Therefore, overcoming periodic jumps and avoiding local extrema in reflection waveform inversion is an urgent problem to be solved. To address the problems of existing technologies, this invention proposes a reflection waveform inversion method based on stochastic gradient sampling. Within the framework of traditional reflection waveform inversion procedures, only a random spatial shift of the wavefield needs to be added. The updated direction is calculated using the wavefield after the random spatial shift, while other aspects remain largely unchanged. The overall process is simple, and the method flow is as follows: Figure 1 As shown.

[0058] S1. Obtain raw observed seismic data and preprocess the raw observed seismic data to obtain real seismic records based on the preprocessed raw observed seismic data.

[0059] In detail, the preprocessing steps for raw observed seismic data in this invention include denoising and filtering.

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

[0061] S3. The reflection waveform inversion method based on stochastic gradient sampling is adopted to update the model based on real seismic records and obtain the final P-wave velocity model.

[0062] S31. Perform forward modeling based on the P-wave velocity model to obtain simulated seismic records. and forward propagation wave field {u n (x,z,t)}, (n=1,2,3…Ns), where g represents the spatial position of the detector, t represents the wave field propagation time, (x,z) represents the spatial 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 k-th iteration, the P-wave velocity model is the P-wave velocity model inverted in the (k-1)-th iteration.

[0063] Specifically, forward modeling is performed using the isotropic acoustic wave equation based on the longitudinal wave velocity model, and its expression is:

[0064]

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

[0066] S32. Calculate the residual seismic record of each shot's real seismic record and simulated seismic record, and calculate the sum of the residual seismic records. If the residual seismic record is less than the 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 as follows: in, This represents the simulated seismic record of the nth shot at spatial location g; This represents the actual seismic record of the nth shot at spatial location g.

[0068] S33. Using each residual seismic record as a companion source, perform end-field backpropagation to obtain the backpropagation seismic record.

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

[0070]

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

[0072] S34. The gradient calculation method based on random spatial movement is used to calculate the update direction of the back-transmitted seismic record and wavelength.

[0073] S341. Simulate the corresponding propagation wave field u for each shot. n (x,z,t) and reverse-propagating seismic records

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

[0075] S343. Randomly move the forward propagation wavefield and reverse propagation seismic records based on the spatial movement range.

[0076] Unlike seismic model inversion calculations that rely on direct waves, where the forward and reverse propagating wavefields are spatially shifted in the same direction, this invention relies solely on reflected waves for model updates. Therefore, calculating the update requires shifting the forward and reverse propagating wavefields in opposite directions. This difference arises because relying on direct waves requires cross-correlating forward-scattered normal direct waves and reverse propagating direct waves in the same direction. However, forward-scattered reflected waves are backscattered, while reverse propagating wavefields are forward-scattered. Therefore, calculating the update differs from calculating the update based on direct waves; it requires shifting the forward and reverse propagating wavefields in opposite directions.

[0077] The specific update method is as follows:

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

[0079]

[0080] Among them, u n (x,z,t) represents the propagating wave field; h represents the range of random spatial movement. This indicates the reverse transmission of earthquake records.

[0081] S344. Calculate the P-wave velocity update direction for each shot based on the forward propagation wavefield and reverse propagation seismic records after random spatial movement, and aggregate them to generate the update direction matrix.

[0082] The method for calculating the updated direction of the nth shot is as follows:

[0083]

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

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

[0086] S351, Calculate the step size:

[0087]

[0088] Where, λ k The step size for the k-th iteration is represented by ε; the scaling factor used to update the model is represented by v; and the P-wave velocity is represented by g. RSS (x,z) represents updating the direction matrix.

[0089] S352, Update the P-wave velocity model:

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

[0091] Among them, v k+1 This represents the P-wave velocity model for the (k+1)th time; v k Denotes the P-wave velocity model for the k-th time; λ k Indicates step size; g RSS (x,z) represents updating the direction matrix.

[0092] Example 2

[0093] In this embodiment, a two-dimensional Gaussian sphere low-velocity anomaly model is used as the real model, such as... Figure 2As shown, the density in the model is 2000 kg / m³ and remains constant. It has 401 × 201 grids with a grid spacing of 10 m × 10 m. A two-dimensional isotropic medium acoustic forward modeling simulation is performed on this model, with a total of 40 shots. The shot points are uniformly distributed on the surface with a shot spacing of 100 m, and the receivers are uniformly distributed on the surface at 10 m intervals. The excitation source is a Ricker wavelet with a dominant frequency of 18 Hz, the seismic record reception length is 2.0 seconds, and the interval is 1 millisecond. A uniform velocity model is set as the initial model for the inversion.

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

[0095] In detail, the update direction (negative gradient) of the first-generation model was calculated using CRWI and RSS-RWI respectively. The results using CRWI are as follows: Figure 3 and Figure 4 As shown, from Figure 1 The provided real model shows that, near the center of the velocity anomaly, it should be updated in the negative velocity direction; however, from... Figure 3 It is evident that the CRWI update direction at the corresponding position is positive, which is contrary to expectations, indicating that the inversion is heading in the wrong direction and a periodic jump phenomenon has occurred; from Figure 4 As can be seen, the CRWI inversion results exhibit periodic jumps, fall into local extrema, and show strong velocity discontinuities, which are far from the actual model, making subsequent seismic interpretation difficult.

[0096] Figure 5 and Figure 6 To use the results of RSS-RWI, from Figure 5 As can be seen, the update direction of RSS-RWI at the corresponding position is negative, which is as expected, making the inverted model closer to the true model; from Figure 6 As can be seen, 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, the method provided by this invention can more accurately invert the reflection waveform compared to conventional inversion methods.

[0098] To further verify the superiority of the method provided in this invention when faced with poor initial models and missing low-frequency data, the reflected wave inversion results after performing RSS-RWI were used as a new initial model to perform CRWI for accuracy verification. The results obtained are shown in the figure. Figure 7 As shown, it can be seen that it is similar to Figure 1The model is very close to the model in the original paper, and there is no periodic jump. Instead, a high-precision and high-resolution inversion model is obtained, which shows that the method provided by this invention has advantages.

[0099] Furthermore, this invention also provides a reflection waveform inversion device based on stochastic gradient sampling, including a central processing unit (CPU), which can execute various appropriate actions and processes according to computer program instructions stored in read-only memory (ROM) or loaded from a storage unit into random access memory (RAM). The RAM can also store various programs and data required for device operation. The CPU, ROM, and RAM are interconnected via a bus. Input / output (I / O) interfaces are also connected to the bus.

[0100] Multiple components in the device are connected to the I / O interface, including: input units such as keyboards and mice; output units such as various types of displays and speakers; storage units such as disks and optical discs; and communication units such as network interface cards (NICs), modems, and wireless transceivers. The communication unit allows the device to exchange information / data with other devices through computer networks such as the Internet and / or various telecommunications networks.

[0101] The processing unit executes 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 computer software programs 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 ROM and / or a communication unit. When the computer program is loaded into 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 suitable means (e.g., by means of firmware).

[0102] The functions described above in this document can be performed, at least in part, by one or more hardware logic components. For example, exemplary types of hardware logic components that can be used, without limitation, include: Field Programmable Gate Arrays (FPGAs), Application-Specific Integrated Circuits (ASICs), Application Standard Products (ASSPs), System-on-Chip (SoCs), Complex Programmable Logic Devices (CPLDs), and so on.

[0103] The program code used to implement the methods of the present invention can be written in any combination of one or more programming languages. This program code can be provided to a processor or controller of a general-purpose computer, special-purpose computer, or other programmable data processing device, such that when executed by the processor or controller, the program code causes the functions / operations specified in the flowcharts and / or block diagrams to be implemented. The program code can be executed entirely on the machine, partially on the machine, as a standalone 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 this invention, a machine-readable medium can be a tangible medium that may contain or store a program for use by or in conjunction with an instruction execution system, apparatus, or device. A machine-readable medium can be a machine-readable signal medium or a machine-readable storage medium. Machine-readable media can include, but are not limited to, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatus, or devices, or any suitable combination of the foregoing. More specific examples of machine-readable storage media include electrical connections based on one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fibers, portable compact disk read-only memory (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 these modifications or substitutions should all be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A reflection waveform inversion method based on stochastic gradient sampling, characterized in that, include: S1. Obtain raw observed seismic data, preprocess the raw observed seismic data, and obtain real seismic records based on the preprocessed raw observed seismic data. S2. Construct the initial P-wave velocity model; S3. Using the reflection waveform inversion method based on stochastic gradient sampling, the model is updated based on the real seismic record to obtain the final P-wave velocity model; The method for obtaining the final P-wave velocity model by iteratively executing the following steps: S31. Perform forward modeling based on the P-wave velocity model to obtain simulated seismic records. And the main Tron ,in, Indicates the spatial position of the detector. Indicates the time of wave field propagation. Represents spatial location coordinates, This indicates the number of shot points; and in the first iteration, the P-wave velocity model is the initial P-wave velocity model; in the... In the nth iteration, the P-wave velocity model is the nth k -1 inversion P-wave velocity model; S32. Calculate the residual seismic record of the real seismic record and the simulated seismic record for 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. S33. Using each of the aforementioned residual seismic records as an adjoint source, perform end-field backpropagation to obtain the backpropagation seismic record. ; S34. The update direction of the back-transmitted seismic record and wavelength is calculated using a gradient calculation method based on random spatial movement. S35. Calculate the step size based on the update direction, update the P-wave velocity model based on the step size, and then return to S31. S34 includes: For each shot, the corresponding forward propagation wave field is simulated. and reverse transmission earthquake records ; Define the range of random spatial movement as: ,in Indicates the proportionality coefficient; Indicates the wavelength of the longitudinal wave and , Indicates the longitudinal wave velocity. The dominant frequency of the wavelet; The forward propagation wavefield and backward propagation seismic records are randomly moved based on the aforementioned spatial movement range; The P-wave velocity update direction for each shot is calculated based on the forward propagation wavefield and the reverse propagation seismic records after random spatial movement, and the updated direction matrix is ​​generated by combining them. The method of random movement is as follows: ， in, This represents the forward propagation wave field after random movement; This represents the back-transmitted seismic record after random movement.

2. The reflection waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, The forward modeling method is as follows: Forward modeling is performed based on the longitudinal wave velocity model using the isotropic acoustic wave equations, the expression of which is: in, This represents the longitudinal wave velocity at spatial coordinate x; Indicates the location of the epicenter. The displacement of the particle at time t; Represents the Laplace operator; This indicates the earthquake source at time t.

3. The reflection waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, The method for calculating the residual seismic record is as follows: ,in, Indicates the spatial position of the nth cannon. g Simulated earthquake records; Indicates the spatial position of the nth cannon. g The actual earthquake records.

4. The reflection waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, The method for end-detector wavefield backpropagation is as follows: using the fourth-order spatial and second-order temporal finite difference method, PML absorbing boundary conditions are implemented at the boundary, and the residual seismic record of each shot is used as the adjoint source to perform end-detector wavefield backpropagation, satisfying the following: in, This represents the longitudinal wave velocity at spatial coordinate x; Indicates the location of the epicenter. The displacement of the particle at time t; Represents the Laplace operator; This indicates the accompanying epicenter, i.e., the residual earthquake record.

5. The reflection waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, The method for calculating the update direction is as follows: ， in, Indicates the update direction of the nth shot; Indicates the longitudinal wave velocity; This represents the forward propagation wave field after random movement; This represents the back-transmitted seismic record after random movement. Represents spatial location coordinates.

6. The reflection waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, The method for calculating the step size includes: ， in, This represents the step size of the k-th iteration; This represents the scaling factor used to update the model; Indicates the longitudinal wave velocity; This indicates updating the direction matrix.

7. The reflection waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, The method for updating the P-wave velocity model is as follows: ， in, This represents the longitudinal wave velocity model for the (k+1)th time. This represents the longitudinal wave velocity model for the kth time. This represents the step size of the k-th iteration; This indicates updating the direction matrix.

Images