A method and system for modeling background velocities using reflected wave travel times from a layer-peeling wave equation

Abstract

This invention belongs to the technical field of exploration seismology and relates to a method and system for modeling background velocity of reflected waves using the layer-stripping wave equation. The method includes: performing reverse-time migration imaging on seismic data based on an initial background velocity model to obtain a depth-domain imaging profile; picking continuous target geological interfaces from the depth-domain imaging profile; performing Rytov approximate reflected wave equation travel-time inversion based on the initial background velocity model and the target geological interfaces to update the background velocity model at the current layer; repeating the above steps until velocity inversion at all target layers is completed. The modeling process in this invention is stable and reliable, overcoming the sensitivity of traditional methods to data noise and travel-time estimation errors, and significantly improving the reliability of the modeling results.

Description

Technical Field

[0001] This invention relates to a method and system for modeling the background velocity of reflected waves in the layer-stripping wave equation, belonging to the technical field of mid-to-deep velocity modeling in exploration seismology. Background Technology

[0002] Seismic velocity modeling is the core and most challenging aspect of seismic imaging. The actual velocity modeling process relies on separating the low-to-medium wavenumber and high wavenumber components of the model. This involves first establishing the macroscopic background velocity using kinematic information, then using imaging or full waveform inversion to establish the high wavenumber components, ultimately creating a velocity model with a full / high wavenumber spectrum. With rapid economic development and significant adjustments to international energy development strategies, the increasing difficulty of onshore exploration and the decreasing availability of resources have led to a shift in exploration focus towards offshore, especially deep-sea oil and gas exploration. As deep-sea oil and gas exploration deepens, wave equation reflected wave travel-time inversion, a crucial technique for mid-to-deep velocity modeling, has gradually become a powerful tool for mid-to-deep background velocity modeling in deep-sea exploration due to its sensitivity to long-wavelength components of the model.

[0003] Traditional data-domain reflected wave travel time inversion implicitly derives the linear relationship between travel time perturbations and velocity perturbations using a link function, i.e., the reflected wave travel time sensitive kernel function, for iterative model updates. However, the link function simultaneously includes travel time and pressure field, causing the traditional reflected wave travel time sensitive kernel function to contain both observed back-reflected wave amplitude and travel time (phase) terms. This means the traditional reflected wave travel time sensitive kernel function is simultaneously affected by variations in field data amplitude and estimated travel time errors. When the data signal-to-noise ratio is low or the estimated reflected travel time difference contains errors, traditional reflected wave travel time inversion may yield erroneous results. Even a slight travel time error can lead to incorrect gradients, ultimately resulting in incorrect velocity inversion results. Summary of the Invention

[0004] To address the aforementioned problems, the present invention aims to provide a method and system for modeling the background velocity of reflected waves in the layer-stripping wave equation. The modeling process is stable and reliable, overcoming the problems of traditional methods being sensitive to data noise and travel time estimation errors, and greatly improving the reliability of the modeling results.

[0005] To achieve the above objectives, the present invention proposes the following technical solution: a method for modeling background velocity of reflected waves using the wave equation of layer stripping, comprising: performing reverse time migration imaging on seismic data based on an initial background velocity model to obtain a depth domain imaging profile; picking continuous target geological interfaces from the depth domain imaging profile; performing Rytov approximate reflected wave equation travel-time inversion based on the initial background velocity model and the target geological interfaces to update the background velocity model at the current layer; repeating the above steps until velocity inversion at all target layers is completed.

[0006] Furthermore, the seismic data is preprocessed, including at least one of the following: denoising, suppression of multiple waves, and removal of direct waves.

[0007] Furthermore, the method for Rytov approximation of the travel time inversion of the reflected wave equation is as follows: Based on the initial background velocity model and the target geological interface, Born forward modeling is performed to obtain the background wave field, scattered wave field, and reflected wave field of the back-migration at the source end; the cross-correlation function of the observed reflected wave and the reflected wave field of the back-migration is calculated, and the time difference is estimated by the maximum value of the cross-correlation function; based on the reflected wave field of the back-migration and the time difference, the accompanying source is calculated by Rytov approximation; based on the accompanying source, the background wave field and scattered wave field of the back propagation at the detector end are calculated; the zero-delay cross-correlation of the forward propagation wave field at the source end and the background wave field and scattered wave field of the back propagation at the detector end is calculated to obtain the gradient of the reflected wave travel time inversion; based on the gradient of the reflected wave travel time inversion, the step size for each update is obtained, and the background velocity model is updated.

[0008] Furthermore, the method for obtaining the background wave field, scattered wave field, and reverse-migrated reflected wave field at the source end is as follows: the background wave field at the source end is obtained by numerically solving the time-domain acoustic wave equation using the finite difference method; the scattered wave field at the source end is obtained by Born forward modeling; the scattered wave field at the detector location is recorded, and the reverse-migrated reflected wave field is obtained.

[0009] Furthermore, the formula for calculating the cross-correlation function is:

[0010] C t (x s ,x r ;τ)=∫d obs (x s ,x r ;t+τ)∫d cal (x s ,x r ;t)dt

[0011] Where, d obs (x s ,x r ;t+τ) represents the observed reflected wave; d cal (x s ,x r ;t) represents the reflected wave field with reverse offset, x s x represents the coordinates of the earthquake source; r Let C be the coordinates of the detector, τ be the time shift for calculating the cross-correlation, and C be the time shift for calculating the cross-correlation. t (x s ,x r ;τ) is the cross-correlation function.

[0012] Furthermore, the time difference is obtained using reverse or forward tracing in the dynamic image warping algorithm. The formula for calculating the time difference is:

[0013]

[0014] Wherein, Δt(x) s ,x r () is the time difference.

[0015] Furthermore, the method for calculating the adjoint source using the Rytov approximation is as follows: take the first and second derivatives of the reverse-shifted reflected wave field with respect to time; obtain the regularization factor of the adjoint source based on the reverse-shifted reflected wave field and the second derivative; and calculate the adjoint source based on the first derivative, the regularization factor, and the time difference.

[0016] Furthermore, the formula for calculating the regularization factor is:

[0017] E=∫ü cal (x s ,x r ;t)u cal (x s ,x r ;t)dt

[0018] The formula for calculating the accompanying source is:

[0019]

[0020] Where E is the regularization factor, u cal (x s ,x r ;t) is the reflected wave field with reverse offset, It is the first derivative of the reflected wave field with reverse offset, ü cal (x s ,x r ;t) is the second derivative of the reflected wave field with reverse offset, Δt(x) s ,x r ) is the time difference, and f(t) is the associated source.

[0021] Furthermore, the formula for calculating the gradient of the travel time inversion of the reflected wave is as follows:

[0022]

[0023] Where v(x) is the velocity value at any spatial position x, ü0(x,t;x s ) is the second derivative of the forward propagating wave field at the source end, üü(x,t;x s ) is the second derivative of the scattered wave field at the source end; λ(x,t; x) s ) is the background wave field propagating backward at the detector end, δλ(x,t;xs ) is the scattered wave field propagating in the reverse direction at the detector end.

[0024] This invention also discloses a layer-stripping wave equation reflection wave travel time background velocity modeling system, comprising: a depth domain imaging profile acquisition module, used to perform reverse time migration imaging on seismic data according to an initial background velocity model to obtain a depth domain imaging profile; a target geological interface picking module, used to pick continuous target geological interfaces from the depth domain imaging profile; a velocity model updating module, used to perform Rytov approximate reflection wave equation travel time inversion according to the initial background velocity model and target geological interfaces to update the background velocity model at the current layer; and a velocity inversion module, used to complete the velocity inversion at all target layers.

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

[0026] 1. The modeling accuracy of the present invention is high: its inversion method is based on the forward modeling engine of the wave equation, which accurately considers the finite frequency band effect in the process of seismic wave propagation, and the theoretical resolution can reach the scale of Fresnel bodies.

[0027] 2. The scheme of the present invention has high computational efficiency and is easy to parallelize: it uses the adjoint state method to calculate the gradient of the travel time inversion of the reflected wave, and the number of wave equation simulations required is only related to the number of sources. There is no need to explicitly calculate the gradient of all shot-receiver pairs, which greatly reduces the computation time.

[0028] 3. The modeling process of the present invention is stable and reliable, overcoming the problem of traditional methods being sensitive to data noise and time estimation errors, and greatly improving the reliability of the modeling results. Attached Figure Description

[0029] Figure 1 This is a flowchart of a method for modeling the background velocity of reflected waves in the layer-peeling wave equation according to an embodiment of the present invention;

[0030] Figure 2 This is a diagram of a real velocity model in one embodiment of the present invention;

[0031] Figure 3 This is an initial velocity model diagram in one embodiment of the present invention;

[0032] Figure 4 This is a surface co-offset atlas image of the initial model in one embodiment of the present invention;

[0033] Figure 5 This is a background velocity modeling result diagram output after three iterations in one embodiment of the present invention;

[0034] Figure 6This is a result of pre-stack depth migration of the background velocity model output after three iterations in one embodiment of the present invention;

[0035] Figure 7 This is a velocity profile diagram of the background velocity output after three iterations in one embodiment of the present invention, at x = 2.25km, 3.75km, 5.25km, and 6.75km respectively;

[0036] Figure 8 This is the surface common offset gather of the background velocity model image output after three iterations in one embodiment of the present invention;

[0037] Figure 9 This is an initial velocity model diagram in another embodiment of the present invention;

[0038] Figure 10 This is a background velocity model diagram output in another embodiment of the present invention;

[0039] Figure 11 This is a speed update graph in another embodiment of the present invention;

[0040] Figure 12 This is a pre-stack depth imaging profile of the initial model in another embodiment of the present invention;

[0041] Figure 13 This is a pre-stack depth imaging profile based on the velocity model obtained using the method of this invention in another embodiment of the present invention. Detailed Implementation

[0042] To enable those skilled in the art to better understand the technical solutions of the present invention, the present invention is 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 terminology used is for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0043] To address the impact of data noise and estimated travel time errors in existing technologies, this invention proposes a method and system for modeling the travel time background velocity of reflected waves using a layer-stripping wave equation. This method is based on the Rytov approximation for layer-stripping reflected wave travel time inversion. Since the Rytov approximation describes the linear relationship between complex phase and model parameter perturbations, it can explicitly derive the linear relationship between travel time perturbations and velocity perturbations. The travel time-sensitive kernel function of the reflected wave based on the Rotov approximation is determined solely by wave theory, eliminating or reducing the impact of data noise and estimated travel time errors on model updates. The Rytov approximation reflected wave travel time inversion proposed in this invention eliminates the influence of data noise and is insensitive to travel time errors, thus providing a stable and reliable method for reflected wave travel time inversion. The following detailed description, in conjunction with the accompanying drawings and embodiments, illustrates the invention in detail.

[0044] Example 1

[0045] This embodiment discloses a method for modeling the background velocity of reflected waves in the layer-stripping wave equation, such as... Figure 1 As shown, it includes:

[0046] S1 uses the initial background velocity model to perform reverse time migration imaging on seismic data to obtain depth domain imaging profiles.

[0047] Seismic data undergoes preprocessing before imaging. The preprocessing process includes at least one of the following: denoising, suppression of multiple waves, and removal of direct waves.

[0048] S2 automatically or manually picks a continuous target geological interface from the depth domain imaging profile.

[0049] S3 performs Rytov approximate reflection wave equation travel-time inversion based on the initial background velocity model and the target geological interface to update the background velocity model at the current stratum;

[0050] The method for performing the travel time inversion of the Rytov approximate reflected wave equation is as follows:

[0051] S3.1 performs Born forward modeling based on the initial background velocity model and the picked target geological interface to obtain the background wave field, scattered wave field and reverse migration reflection wave field at the source end, and stores the background wave field, scattered wave field and reverse migration reflection wave field data at the source end.

[0052] The method for obtaining the background wave field, scattered wave field, and reverse-migrated reflected wave field at the source end is as follows:

[0053] The background wave field at the source end is obtained by numerically solving the time-domain acoustic wave equation using the finite difference method.

[0054] While obtaining the background wave field at the source end through time-axis extension, the scattered wave field at the source end is simultaneously obtained through Born forward modeling.

[0055] Record the scattered wave field at the detector location and obtain the reverse-shifted reflected wave field.

[0056] S3.2 Calculate the cross-correlation function of the observed reflected wave and the reflected wave field of the reverse migration, and automatically estimate the time difference by using the maximum value of the cross-correlation function.

[0057] Within the shot-receiver cycle, the cross-correlation function for each shot-receiver pair is calculated. The formula for calculating this cross-correlation function is as follows:

[0058] C t (x s ,x r ;τ)=∫d obs (x s ,x r ;t+τ)∫d cal (x s ,x r ;t)dt

[0059] Where, d obs (x s ,x r ;t+τ) represents the observed reflected wave; d cal (x s ,x r ;t) represents the reflected wave field with reverse offset, x s x represents the coordinates of the earthquake source; r Let τ be the coordinates of the detector, and τ be the time shift for calculating the cross-correlation.

[0060] Through the cross-correlation function C t (x s ,x r The maximum value of τ is used to automatically estimate the time difference. A key to estimating a stable and reliable time difference is to limit the lateral variation of the time difference to be too large; the lateral variation of the time difference should be less than a pre-set threshold. In this embodiment, the reverse or forward tracing technique in the dynamic image warping (DIW) algorithm is used to obtain the time difference during the automatic estimation process. The specific formula for calculating the time difference is:

[0061]

[0062] In the above formula, when x s =x r When, Δt(x) s ,x r = 0. This relationship is determined by the zero offset travel time invariant, which theoretically guarantees the stability of the reflection time difference estimation method proposed in this embodiment. The time difference obtained by the travel time estimation method in this embodiment is the globally optimal path, which can avoid the cumbersome manual picking of reflected wave travel time.

[0063] S3.3 Based on the reflected wave field and time difference of the reverse migration, the accompanying source is calculated by the Rytov approximation;

[0064] The method for calculating the adjoint source using the Rytov approximation is as follows:

[0065] The reflected wave field u is offset in the opposite direction cal (x s ,x r The first derivative of t with respect to time and second derivative ü cal (x s ,x r ;t);

[0066] The regularization factor of the adjoint source is obtained from the reflected wave field of the inverse offset and the second derivative.

[0067] The formula for calculating the regularization factor is:

[0068] E=∫ü cal (x s ,x r ;t)u cal (x s ,x r ;t)dt

[0069] The adjoint source is calculated based on the first derivative, regularization factor, and time difference. The formula for calculating the adjoint source is as follows:

[0070]

[0071] Where E is the regularization factor, u cal (x s ,x r ;t) is the reflected wave field with reverse offset, It is the first derivative of the reflected wave field with reverse offset, ü cal (x s ,x r ;t) is the second derivative of the reflected wave field with reverse offset, Δt(x) s ,x r ) is the time difference, and f(t) is the associated source.

[0072] To illustrate the advantages of the method for calculating the accompanying source in this embodiment, the following describes a prior art method for calculating the accompanying source, with the formula for calculating the accompanying source as follows:

[0073]

[0074] Comparing the accompanying source f(t) obtained in this embodiment with the accompanying source f in the prior art ′(t) It can be observed that the existing method includes the amplitude term of the reflection data, i.e., u obs (x s ,x r Since f′(t) is affected by noise in the data, noise in the data can negatively impact the inversion of traditional methods. For data with low signal-to-noise ratio, traditional reflected wave travel time inversion methods may yield incorrect inversion results. On the other hand, f′(t) contains a time delay (phase delay) term, Δt. The error in the reflected wave travel time estimation will result in a time-shifted (phase-shifted) signal, which may lead to an incorrect reflected wave travel time gradient, affecting the final inversion result. In practical applications, noise in the data can often be largely eliminated through preprocessing, and the impact of travel time estimation error on gradient calculation is the main factor. Due to the time-varying characteristics of reflected wave travel time, estimating these time-varying characteristics is a highly challenging problem, as the estimated reflected wave travel time difference often contains significant errors.

[0075] Compared with existing technologies, the accompanying source in this embodiment is not affected by noise in the data. More importantly, the new accompanying source does not contain time delay or phase delay terms, and therefore is not sensitive to errors in travel time estimation. The scheme in this embodiment provides a stable and reliable method for modeling the background velocity of reflected waves in the wave equation.

[0076] S3.4 Calculate the back-propagating background wave field and scattered wave field at the detector end based on the accompanying source;

[0077] S3.5 Calculate the zero-delay cross-correlation between the forward propagating wave field at the source end and the reverse propagating background wave field and scattered wave field at the detector end to obtain the gradient of the reflected wave travel time inversion.

[0078] The formula for calculating the gradient of the travel time inversion of the reflected wave is:

[0079]

[0080] Where v(x) is the velocity value at any spatial position x, ü0(x,t;x s ) is the second derivative of the forward propagating wave field at the source end, ü(x,t;x s ) is the second derivative of the scattered wave field at the source end; λ(x,t; x) s ) is the background wave field propagating backward at the detector end, δλ(x,t;x s ) is the scattered wave field propagating in the reverse direction at the detector end.

[0081] S3.6 Calculate the update step size for each step based on the gradient obtained from the travel time inversion of the reflected wave, and update the background velocity model.

[0082] S4 Repeat steps S1-S4 until velocity inversion is completed at all target levels.

[0083] As exploration and development work progresses, seismic exploration is gradually expanding from land to the ocean, especially the deep sea. Towed cable observation systems are the most economical and efficient acquisition method, and the pressure data collected by these systems constitutes the main body of marine seismic data. However, due to the limited length of the towed cables, marine seismic data often has a small observation aperture, making it difficult to obtain the structure of mid-to-deep target areas using conventional first-arrival modeling techniques. Therefore, reflection wave modeling methods have become the main technical tool for mid-to-deep velocity modeling. Since traditional wave equation modeling techniques are sensitive to data noise and travel time estimation errors, this embodiment provides a more stable and reliable wave equation reflection wave travel time mid-to-deep modeling method, meeting the common needs of both field exploration and laboratory data processing.

[0084] Example 2

[0085] Based on the same inventive concept, such as Figure 2 In this embodiment, a two-dimensional complex theoretical model is used as the real model. This model has 750×240 grids with a grid spacing of 15m×15m, and the maximum and minimum velocities are 2000m / s and 4200m / s, respectively. A forward modeling simulation of acoustic waves is performed on this model, with a total of 125 shots. Shot points are evenly distributed on the ground surface with a shot spacing of 90m, and the first shot is at a horizontal position of 0m. Detector points are evenly distributed on the ground surface at 15m intervals. Shots are fired from both sides, with a maximum offset of 3600m and a minimum offset of 0m per shot. Three target layers were selected in the modeling process. Specifically:

[0086] Seismic data is acquired and preprocessed to avoid damaging the kinematic characteristics of the seismic data, such as filtering and suppressing multiple waves, and direct waves are removed.

[0087] An initial background velocity model is generated based on prior subsurface information. The preprocessed seismic data is then input into the processor, which performs reverse time migration imaging based on the initial background velocity model to obtain a depth domain imaging profile.

[0088] Automatically or manually pick a continuous target geological interface from a depth domain imaging profile;

[0089] Born forward modeling was performed based on the initial background velocity model and the picked target geological interface to obtain the background wave field, scattered wave field and reverse migration reflected wave field at the source end.

[0090] Based on the reflected wave field and time difference of the reverse offset, the accompanying source is calculated by the Rytov approximation, and the background wave field and scattered wave field of the back propagation at the detector end are calculated.

[0091] The zero-delay cross-correlation between the forward propagating wave field at the source end and the reverse propagating background and scattered wave fields at the detector end is calculated to obtain the gradient for the travel time inversion of the reflected wave.

[0092] The step size for each update is obtained based on the gradient obtained from the travel time inversion of the reflected wave, and the background velocity model is updated accordingly.

[0093] Repeat the above steps until velocity inversion is completed at all target levels. The velocity inversion results are displayed through output devices, such as monitors.

[0094] Figure 3 This is the initial velocity model diagram in this embodiment. The following is obtained through this initial velocity model diagram: Figure 4 The surface common offset atlas shown is shown. Figure 3 The fact that the phase axes in the middle generally tilt upwards indicates that the initial model velocity is too low, stemming from a coarse, constant gradient initial model (i.e., ...). Figure 3 Starting from the model, the inversion results of three rounds of iterative Rytov approximation reflection wave travel time inversion in this embodiment are shown in the figure below. Figure 5 As shown, the depth migration imaging results obtained using the inverted velocity model are shown in the figure. Figure 6 , refer to Figure 2 The real model shows that the velocity modeling results are highly similar to the real model, and most of the velocity layers can be reconstructed.

[0095] Figure 7 These are velocity profiles of the background velocity output after three iterations in this embodiment, at x = 2.25km, 3.75km, 5.25km, and 6.75km respectively; from Figure 7 The inversion results, specifically the velocity profiles at x = 2.25 km, 3.75 km, 5.25 km, and 6.75 km, show that the background velocity obtained after three iterations using the method in this embodiment is very close to the true velocity, as indicated by the solid line in the figure. The inversion results from three iterations (i.e.,...) Figure 5 The obtained surface co-offset atlas is shown in the following figure. Figure 8 As shown in the figure, the in-phase axes are generally flattened, indicating that the inversion results obtained in this embodiment are highly reliable.

[0096] It should be noted that this embodiment successfully reconstructs the background velocity information in the mid-to-deep regions using only reflected waves with a small offset. The process in this embodiment does not involve manually picking up the travel time of the reflected waves, nor does it involve manually picking up the residual curvature of the imaging gather. It is highly automated, stable, and reliable, demonstrating the great potential of the method in practical applications.

[0097] Example 3

[0098] Based on the same inventive concept, this embodiment uses towed cable data collected in a certain work area in the East my country Sea as an example to verify the method for modeling the background velocity of reflected waves in the mid-layer stripping wave equation of the present invention. Figure 9 This is the initial velocity model diagram in this embodiment. This initial velocity model is a constant gradient model with a grid spacing of 12.5 meters × 12.5 meters. The towed cable data consists of 360 shots with a shot spacing of 225 meters. Each shot is received by 315 detectors with a channel spacing of 12.5 meters. The maximum offset of the observation system is 4000 meters, and the minimum offset is 0 meters. The specific process is similar to that in Embodiment 2, the only difference being that only one target layer is selected in this embodiment. The background velocity model diagram output in this embodiment is as follows: Figure 10 As shown, the speed update amount in this embodiment is as follows: Figure 11 As shown.

[0099] from Figure 10 and Figure 11 As can be seen from the data, the velocity has been significantly updated using the velocity modeling method in this embodiment, with a maximum value of approximately 500 m / s. The velocity update also indicates that the initial model's velocity was significantly underestimated. Figure 12 This is a pre-stack depth imaging profile of the initial model in this embodiment; from Figure 12 As can be seen, the pre-stack depth imaging profile of the initial model has poor continuity, making it difficult to identify geologically significant stratigraphic layers, and contains obvious artifacts of offset arcs caused by low velocities. Further interpretation is impossible based on this profile. This invention uses a mid-layer stripping wave equation reflected wave travel time background velocity modeling method to generate a velocity model. The pre-stack depth imaging profile generated based on this velocity model is shown below. Figure 13 As shown. From Figure 13 As can be seen from the figure, the continuity of the in-phase axis is greatly improved, and the in-phase axis of the imaging profile shows strong energy. Similarly, the method in this embodiment does not involve manually picking up the travel time of the reflected wave, nor does it involve manually picking up the residual curvature of the imaging gather. It has a high degree of automation, is stable and reliable, and proves the effectiveness of the invention in practical applications.

[0100] Example 4

[0101] Based on the same inventive concept, this embodiment discloses a system for modeling the travel time and background velocity of reflected waves in a layer-stripping wave equation, comprising:

[0102] The depth domain imaging profile acquisition module is used to perform reverse time migration imaging on seismic data based on the initial background velocity model to obtain a depth domain imaging profile.

[0103] The target geological interface picking module is used to pick up continuous target geological interfaces from depth domain imaging profiles.

[0104] The velocity model update module is used to perform Rytov approximate reflection wave equation travel-time inversion based on the initial background velocity model and the target geological interface, and update the background velocity model at the current layer.

[0105] The velocity inversion module is used to perform velocity inversion at all target levels.

[0106] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0107] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0108] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0109] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0110] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific embodiments of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the protection scope of the claims of the present invention. The above content is only a specific embodiment of this application, but the protection scope of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in this application should be covered within the protection scope of this application. Therefore, the protection scope of this application should be determined by the protection scope of the claims.

Claims

1. A method for modeling the background velocity of reflected waves in a layer-stripping wave equation, characterized in that, include: Based on the initial background velocity model, the seismic data is subjected to reverse time migration imaging to obtain a depth domain imaging profile. Pick up continuous target geological interfaces from the depth domain imaging profile; Based on the initial background velocity model and the target geological interface, the travel-time inversion of the Rytov approximate reflected wave equation is performed to update the background velocity model at the current stratum. Repeat the above steps until velocity inversion is completed at all target levels; The method for performing the travel-time inversion of the Rytov approximate reflected wave equation is as follows: Based on the initial background velocity model and the target geological interface, Born forward modeling was performed to obtain the background wave field, scattered wave field and reverse-migrated reflected wave field at the source end; Calculate the cross-correlation function between the observed reflected wave and the reverse-shifted reflected wave field, and estimate the time difference using the maximum value of the cross-correlation function; Based on the reflected wave field and time difference of the reverse offset, the accompanying source is calculated by the Rytov approximation. The back-propagating background wave field and scattered wave field at the detector end are calculated based on the accompanying source. The zero-delay cross-correlation between the forward propagation wave field at the source end and the reverse propagation background wave field and scattered wave field at the detector end is calculated to obtain the gradient of the travel time inversion of the reflected wave. The step size for each update is obtained based on the gradient of the reflected wave travel time inversion, and the background velocity model is updated accordingly. The method for obtaining the background wave field, scattered wave field, and reverse-migrated reflected wave field at the source end is as follows: The background wave field at the source end is obtained by numerically solving the time-domain acoustic wave equation using the finite difference method. The source-side scattered wave field was obtained by Born forward modeling. Record the scattered wave field at the detector location and obtain the reverse-shifted reflected wave field.

2. The method for modeling the background velocity of reflected waves in the layer-peeling wave equation as described in claim 1, characterized in that, The seismic data is preprocessed, and the preprocessing process includes at least one of the following: denoising, suppression of multiple waves, and removal of direct waves.

3. The method for modeling the travel time and background velocity of reflected waves according to the layer-peeling wave equation as described in claim 1, characterized in that, The formula for calculating the cross-correlation function is as follows: in, To observe the reflected wave; For the reverse-shifted reflected wave field, The coordinates of the earthquake source; The coordinates of the detector are... It calculates the time shift of cross-correlation. It is a cross-correlation function.

4. The method for modeling the travel time and background velocity of reflected waves according to the layer-peeling wave equation as described in claim 1, characterized in that, The time difference is obtained using the reverse or forward tracing method in the dynamic image warping algorithm, and the formula for calculating the time difference is: in, It's the time difference.

5. The method for modeling the background velocity of reflected waves in the layer-peeling wave equation as described in claim 1, characterized in that, The method for calculating the adjoint source using the Rytov approximation is as follows: The first and second derivatives of the reverse-shifted reflected wave field with respect to time; The regularization factor of the adjoint source is obtained based on the reflected wave field of the reverse offset and the second derivative. The associated source is calculated based on the first derivative, the regularization factor, and the time difference.

6. The method for modeling the travel time and background velocity of reflected waves according to the layer-peeling wave equation as described in claim 5, characterized in that, The formula for calculating the regularization factor is as follows: The formula for calculating the accompanying source is: in, It is a regularization factor. It is a reverse-shifted reflected wave field. It is the first derivative of the reflected wave field with the reverse offset. It is the second derivative of the reflected wave field with the reverse offset. It's the time difference. It is an accompanying source.

7. The method for modeling the travel time and background velocity of reflected waves according to the layer-peeling wave equation as described in claim 6, characterized in that, The formula for calculating the gradient of the travel time inversion of the reflected wave is as follows: in, It is any spatial location The velocity value at that location, It is the second derivative of the forward-propagating wave field at the source end. It is the second derivative of the scattered wave field at the source end; It is the background wave field propagating in the back direction at the detector end. It is the scattered wave field propagating in the reverse direction at the detector end.

8. A system for modeling the travel time and background velocity of reflected waves in a layer-stripping wave equation, characterized in that, include: The depth domain imaging profile acquisition module is used to perform reverse time migration imaging on seismic data based on the initial background velocity model to obtain a depth domain imaging profile. The target geological interface picking module is used to pick up continuous target geological interfaces from the depth domain imaging profile. The velocity model update module is used to perform Rytov approximate reflection wave equation travel-time inversion based on the initial background velocity model and the target geological interface, and update the background velocity model at the current layer. The velocity inversion module is used to perform velocity inversion at all target levels; The method for performing the travel-time inversion of the Rytov approximate reflected wave equation is as follows: Based on the initial background velocity model and the target geological interface, Born forward modeling was performed to obtain the background wave field, scattered wave field and reverse-migrated reflected wave field at the source end; Calculate the cross-correlation function between the observed reflected wave and the reverse-shifted reflected wave field, and estimate the time difference using the maximum value of the cross-correlation function; Based on the reflected wave field and time difference of the reverse offset, the accompanying source is calculated by the Rytov approximation. The back-propagating background wave field and scattered wave field at the detector end are calculated based on the accompanying source. The zero-delay cross-correlation between the forward propagation wave field at the source end and the reverse propagation background wave field and scattered wave field at the detector end is calculated to obtain the gradient of the travel time inversion of the reflected wave. The step size for each update is obtained based on the gradient of the reflected wave travel time inversion, and the background velocity model is updated accordingly. The method for obtaining the background wave field, scattered wave field, and reverse-migrated reflected wave field at the source end is as follows: The background wave field at the source end is obtained by numerically solving the time-domain acoustic wave equation using the finite difference method. The source-side scattered wave field was obtained by Born forward modeling. Record the scattered wave field at the detector location and obtain the reverse-shifted reflected wave field.

Images