Reflective focusing-based depth-velocity modeling method and apparatus

Abstract

A reflective focusing-based depth-velocity modeling method, relating to the field of seismic exploration, comprising: determining an actual seismic record and an initial velocity model (S1); on the basis of a seismic observation system, performing forward simulation on the initial velocity model, and obtaining forward wave field data (S2); on the basis of the actual seismic record and the forward wave field data, calculating a focusing function of a reflected wave residual response, calculating a Green's function of the reflected wave residual response on the basis of the focusing function, and, on the basis of the Green's function and the forward wave field data, calculating a gradient of an objective function with respect to a longitudinal wave velocity (S3); determining an updated velocity model on the basis of the gradient of the objective function with respect to the longitudinal wave velocity and an update step size (S4); and on the basis of the actual seismic record and the updated velocity model, determining whether a termination condition has been met (S5), if so, using the finally obtained updated velocity model as a depth-velocity model, and if not, re-performing forward simulation and the following steps. Further related is a reflective focusing-based depth-velocity modeling apparatus.

Description

A depth velocity modeling method and apparatus based on reflection focusing

[0001] Related applications

[0002] This application claims Chinese Patent Application No. 202411178091.8, filed on August 26, 2024, and incorporates the disclosure of the aforementioned patent application as part of this application. Technical Field

[0003] This application relates to the field of seismic exploration, and in particular to a depth velocity modeling method and apparatus based on reflection focusing. Background Technology

[0004] Underground oil and natural gas reserves are important exploration and development targets for meeting my country's energy needs. Seismic exploration is a crucial method for detecting underground oil and natural gas reservoirs. However, the geophysical field has consistently faced the challenge of obtaining low-frequency information in seismic exploration.

[0005] Traditional seismic imaging methods rely on tomography and ray theory, and their results are limited by high-frequency approximation assumptions. These assumptions require that the seismic wave wavelength be relatively small relative to the model's scale of change, that the stratigraphic interfaces be relatively smooth, and that the lateral variation of the seismic wave propagation velocity be small. Therefore, in complex tectonic regions, the high-frequency approximation assumptions do not hold. When traditional seismic imaging methods are applied to complex tectonic regions, they can lead to inaccurate estimations of small-scale anomalies and lateral stratigraphic variations.

[0006] Full waveform inversion (FWI), while attractive as a method for acquiring high-resolution subsurface imaging, is limited in its application in reflection seismology due to its reliance on low-frequency information and a sufficiently accurate initial model. Furthermore, while mid-to-deep velocity information is often contained within reflection wave data, traditional full waveform inversion techniques have limitations in modeling mid-to-deep velocities. To address this issue, reflection waveform inversion (RWI) techniques have been extensively studied in the international seismic exploration field in recent years.

[0007] Reflection wave waveform inversion is a modeling technique in geophysical exploration. Its main purpose is to reconstruct the background velocity model at depth using reflected wave information. Seismic data recorded by seismic instruments contains complex responses of the subsurface medium to seismic waves, including reflection and scattering. Reflection wave waveform inversion involves adjusting the parameters of the subsurface model to make the reflected wave data generated by forward modeling as consistent as possible with the observed reflected wave data. This is achieved by utilizing the differences between the two and using iterative optimization algorithms to modify the velocity along the seismic wave path, iterating repeatedly to achieve a velocity model that meets the required accuracy.

[0008] Current Reflection-In-Wiring (RWI) typically uses migration / inverse migration algorithms to generate reflected waves, a crucial step in RWI's deep modeling process. The reflected wave residuals are then back-projected into the model space for iterative updates to the velocity model. This requires additional migration / inverse migration imaging processing. This technique necessitates the interaction between the background wavefield and the migration imaging results to simulate the required reflected waves; the accuracy of the current velocity model also affects the accuracy of the inverse migration simulation. Therefore, traditional RWI using migration / inverse migration algorithms suffers from the following technical problems: the requirement for accurate simulation of reflected waves is constrained by the accuracy of the background velocity; it requires high-frequency boundary information in the model or necessitates additional migration imaging processing, increasing computational complexity. Summary of the Invention

[0009] This application addresses the problem that existing methods for inverting reflected wave waveforms use migration / inverse migration algorithms to generate reflected waves, which are limited by the accuracy of background velocity, resulting in low accuracy in mid-to-deep velocity modeling. Furthermore, these methods require high-frequency boundary information or additional migration imaging processing, leading to computational complexity.

[0010] To address the aforementioned technical problems, the first aspect of this application provides a depth velocity modeling method based on reflection focusing, comprising:

[0011] Determine the actual earthquake record and initial velocity model;

[0012] Based on the earthquake observation system, forward modeling of the initial velocity model is performed to obtain forward wavefield data;

[0013] Based on the actual seismic records and the forward wavefield data, the focusing function of the reflected wave residual response is calculated using the Marchenko equation; the Green's function of the reflected wave residual response is calculated based on the focusing function; and the gradient of the objective function with respect to the P-wave velocity is calculated based on the Green's function and the forward wavefield data.

[0014] The update velocity model is determined based on the gradient of the objective function with respect to the P-wave velocity and the update step size;

[0015] Based on the actual seismic records and the updated velocity model, determine whether the termination condition is met. If yes, use the final updated velocity model as the depth velocity model; otherwise, re-execute the forward modeling and subsequent steps.

[0016] A second aspect of this application provides a depth velocity modeling device based on reflection focusing, comprising:

[0017] Modeling units are used to determine the actual seismic record and the initial velocity model;

[0018] The forward modeling unit is used to perform forward modeling of the initial velocity model based on the seismic observation system to obtain forward wavefield data;

[0019] The gradient calculation unit is used to calculate the focusing function of the reflected wave residual response based on the Marchenko equation according to the actual seismic record and the forward wave modeling data; calculate the Green's function of the reflected wave residual response according to the focusing function; and calculate the gradient of the objective function with respect to the P-wave velocity according to the Green's function and the forward wave modeling data.

[0020] The correction unit is used to determine the updated velocity model based on the gradient of the objective function with respect to the P-wave velocity and the update step size;

[0021] The judgment unit is used to determine whether the termination condition is met based on the actual seismic record and the updated velocity model. If yes, the final updated velocity model is used as the depth velocity model; otherwise, the forward modeling unit, gradient calculation unit, and correction unit are restarted.

[0022] A third aspect of this application provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the method described in any of the foregoing embodiments.

[0023] A fourth aspect of this application provides a computer storage medium having a computer program stored thereon, which, when executed by a processor of a computer device, implements the method described in any of the foregoing embodiments.

[0024] The fifth aspect of this application provides a computer program product comprising a computer program that, when executed by a processor of a computer device, implements the method described in any of the foregoing embodiments.

[0025] The depth velocity modeling method and apparatus based on reflection focusing provided in this application obtains forward wavefield data by performing forward modeling of an initial velocity model based on a seismic observation system. Based on actual seismic records and forward wavefield data, the focusing function of the reflected wave residual response is calculated based on the Marchenko equation, and the Green's function of the reflected wave residual response is calculated based on the focusing function. The gradient of the objective function with respect to the P-wave velocity is calculated based on the Green's function and the forward wavefield data. The Green's function can be used to replace the reflected wave response wavefield obtained by migration / inverse migration algorithms, without being constrained by the accuracy of the background velocity, and without the need for additional computationally intensive migration / inverse migration processing, thus improving the determination accuracy and efficiency of the velocity model.

[0026] To make the above and other objects, features and advantages of this application more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. Attached Figure Description

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

[0028] Figure 1 shows a flowchart of the depth velocity modeling method based on reflection focusing according to an embodiment of this application;

[0029] Figure 2 shows a schematic diagram of the longitudinal wave velocity model according to an embodiment of this application;

[0030] Figure 3 shows a schematic diagram of the initial velocity model of an embodiment of this application;

[0031] Figure 4 shows a flowchart of the forward modeling wavefield data determination process according to an embodiment of this application;

[0032] Figure 5 shows a flowchart of the focus function calculation process according to an embodiment of this application;

[0033] Figure 6 shows a schematic diagram of reconstructing the Green's function of the ascending wave according to an embodiment of this application;

[0034] Figure 7 shows a schematic diagram of reconstructing the down-wave Green's function according to an embodiment of this application;

[0035] Figure 8 shows a flowchart of the calculation process of the gradient of the objective function with respect to the longitudinal wave velocity in an embodiment of this application;

[0036] Figure 9 shows a schematic diagram of the gradient of the objective function with respect to the longitudinal wave velocity at any spatial point in an embodiment of this application;

[0037] Figure 10 shows a flowchart of the update determination process according to an embodiment of this application;

[0038] Figure 11 shows a schematic diagram of the final speed model established in the embodiment of this application;

[0039] Figure 12 shows a schematic diagram of the velocity model established by the conventional waveform inversion method in an embodiment of this application;

[0040] Figure 13 shows a structural diagram of the depth velocity modeling device based on reflection focusing according to an embodiment of this application.

[0041] Figure 14 shows a structural diagram of a computer device according to an embodiment of this application. Detailed Implementation

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

[0043] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this application described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, apparatus, product, or device that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or devices.

[0044] This specification provides the operational steps of the methods described in the embodiments or flowcharts, but based on conventional or non-inventive labor, more or fewer operational steps may be included. The order of steps listed in the embodiments is merely one possible execution order among many and does not represent the only possible execution order. In actual system or device products, the methods shown in the embodiments or drawings can be executed sequentially or in parallel.

[0045] In one embodiment of this application, a depth velocity modeling method based on reflection focusing is provided to address the problems of low accuracy in mid-to-deep velocity modeling caused by existing reflection wave waveform inversion methods that use migration / inverse migration algorithms to generate reflection waves, which are limited by the accuracy of background velocity. Furthermore, these methods require high-frequency boundary information or additional migration imaging processing, resulting in computational complexity. Specifically, as shown in Figure 1, the depth velocity modeling method based on reflection focusing includes:

[0046] S1, determine the actual earthquake record and initial velocity model.

[0047] In detail, actual seismic records can be obtained through reflection wave receiving devices installed in deep underground reservoirs, including at least reflection wave data. Specifically, the reflection wave receiving device can be, for example, a reflection wave receiving sensor. This sensor, installed in the deep underground reservoir, receives reflection wave data and shot records, and then sends this data to a server. Alternatively, it can obtain actual seismic records based on the reflection wave data and then send these records to the server. The server is communicatively connected to the reflection wave receiving sensor, receiving the shot records and reflection wave data or actual seismic records sent by the sensor. Upon receiving reflection wave data, the server determines the actual seismic data based on this data. The initial velocity model is obtained by analyzing the shot records received by the reflection wave receiving device through server analysis. The specific analysis process includes: processing the shot records using methods including but not limited to denoising, static correction, dynamic correction, deconvolution, etc.; performing stacking velocity analysis on the processed data using NMO (Normal Moveout Correction) to obtain the P-wave velocity model; and converting the P-wave velocity model to the depth domain to obtain the initial velocity model. In this step, the P-wave model is transformed into the depth domain using the Dix equation.

[0048] In some implementations, the longitudinal wave velocity model is shown in Figure 2, and the initial velocity model is shown in Figure 3. In Figures 2 and 3, the horizontal axis represents distance, the vertical axis represents depth, and the legend on the right represents velocity.

[0049] S2, based on the seismic observation system, forward modeling of the initial velocity model is performed to obtain forward wavefield data.

[0050] In detail, an earthquake observation system describes the arrangement of the relative spatial positions of the earthquake source and receiver in earthquake exploration.

[0051] Forward modeling wavefield data includes: the first simulated wavefield u(x,t) generated at any source. i |x S The second simulated wave field u excited at the detector d (x,t i |x R ) and simulated seismic records R at the detector cal (x R ,t i |x S The direct wave u determined according to the second simulated wave field d (x R ,t i |x), the simulated reflected wave field δu(x) obtained from the second simulated wave field j ,t i |x sSeismic wave field refers to the spatial and temporal distribution of seismic waves as they propagate through the subsurface medium.

[0052] Forward modeling is an important numerical simulation technique in geophysics, simulating the propagation of seismic waves in subsurface media. In this step, the finite-difference forward modeling method based on the scalar wave equation can be used to simulate the forward wave field during forward propagation in time. Forward propagation in time refers to the change in time from an earlier moment to a later moment.

[0053] The finite difference forward modeling method for scalar wave equations is a numerical simulation method used to simulate the propagation of sound waves in underground media. The finite difference method transforms the partial differential equations in the scalar wave equations into finite difference calculations, dividing the continuous spatiotemporal domain into discrete grid points, and simulating sound wave propagation by solving the difference equations.

[0054] S3. Based on actual seismic records and forward wavefield data, calculate the focusing function of the reflected wave residual response using the Marchenko equation; calculate the Green's function of the reflected wave residual response based on the focusing function; and calculate the gradient of the objective function with respect to the P-wave velocity based on the Green's function and forward wavefield data.

[0055] In detail, the focusing function includes the upward focusing function f + (x R ,t i |x) and downward focusing function f - (x R ,t i |x). The Green's function includes the ascending wave Green's function G of the reflected wave residual response. + (x,t i ;x R The Green's function G of the down-flowing wave in the residual response of the reflected wave. - (x,t nt -t i ;x R Where x represents any point in space, t i Indicates the sampling time, nt represents the number of seismic data samples corresponding to time t, and x represents the sampling time. S The x-coordinate represents the spatial location coordinates of the earthquake source. R This indicates the spatial coordinates of the detector.

[0056] In geophysical exploration, seismology, and rock physics, the Green's function is a mathematical function that represents the seismic wave response at one point in response to an earthquake or other wave excitation at another point. The reflected wave from the ground and its corresponding Green's function can be correlated using a focusing function, which focuses the seismic wave as it propagates from the ground into the subsurface medium. The relationship between the Green's function, the focusing function, and the reflected wave is expressed as follows:

[0057] In the above formula, t is a certain sampling time, t' is any time when the sampling process is integrated, G is the Green's function, f is the focusing function, R is the reflected response received by the ground surface, and x s ,x R The coordinates of the source and the detector are represented, respectively, and x represents the spatial coordinates of any point in the velocity model. The superscript "+" represents a downflow wave, and the superscript "-" represents an upflow wave. Upflow and downflow waves are two basic wave forms in seismic wave propagation. An upflow wave is a seismic wave that propagates upwards along the subsurface structure. A downflow wave is a seismic wave that propagates from the source towards the subsurface structure. The Green's function G is defined as a solution to the scalar wave equation; for example, G(rec,t|sou) represents the impulse point source response to the source sou observed at the observation location rec. The symbol after "|" indicates the source location for this Green's function.

[0058] The focusing function f and the reflected wave R(x) R ,t|x S The relationship between these factors is represented by the Marchenko equation, which forms the basis for determining the focusing function. Once the focusing function is obtained by solving the Marchenko equation, the decomposed Green's function G can be retrieved by substituting the upward and downward focusing functions into Equations 1 and 2. + and G - .

[0059] S4. Determine the update velocity model based on the gradient of the objective function with respect to the P-wave velocity and the update step size.

[0060] In this step, the model update amount is defined by the product of the gradient of the objective function with respect to the P-wave velocity and the update step size. The updated velocity model is obtained by subtracting the model update amount from the current model. Specifically, the velocity model is updated using the following formula: v(x) n+1 =v(x) n -α n ·g(x) n ;

[0061] Where n represents the number of iterations, v(x) n+1 Let v(x) represent the model for the (n+1)th iteration. n Let α represent the velocity model for the nth iteration. n Let g(x) represent the step size. n This represents the gradient of the objective function with respect to the P-wave velocity at the nth iteration.

[0062] In some implementations, the step size α n It can be calculated using parabolic interpolation.

[0063] S5: Determine whether the termination condition is met based on the actual seismic records and the updated velocity model. If yes, use the final updated velocity model as the depth velocity model. If not, repeat steps S2 to S5.

[0064] Specifically, the termination condition could be, for example, that the value of the objective function corresponding to the actual seismic record and update velocity model is less than a threshold value. If the value is less than this threshold, the actual seismic record and update velocity model is deemed to meet the termination condition. It should be noted that other conditions can also be used, and this specification does not limit this. Further, after determining the depth velocity model, seismic imaging is performed using this model to obtain subsurface images, which are then displayed to the user via interaction with a display device. Alternatively, after determining the depth velocity model, seismic imaging can be performed using this model to obtain subsurface images to determine subsurface oil and gas reservoir information. This subsurface oil and gas reservoir information could include, for example, whether an oil and gas reservoir is included, and / or, the oil and gas reservoir reserves.

[0065] This application can use the Green's function to replace the reflected wave response wave field obtained by the migration / inverse migration algorithm, without being constrained by the accuracy of the background velocity, and without the need for additional migration / inverse migration processing with huge computational cost, thereby improving the determination accuracy and efficiency of the velocity model.

[0066] In one embodiment of this application, as shown in Figure 4, the above-mentioned S2 performs forward modeling of the initial velocity model based on the seismic observation system to obtain forward wavefield data, including:

[0067] S21, Based on the earthquake observation system, perform forward modeling on the initial velocity model to obtain the velocity data for any earthquake source x. S The first simulated wave field u(x,t) excited at point i |x S The second simulated wave field u excited at detector xR d (x,t i |x R And extract the simulated seismic record R at the detector. cal (x R ,t i |x S ).

[0068] Where, x S The x-coordinate represents the spatial location coordinates of the earthquake source. R The coordinates of the detector's spatial position are represented by x, and the coordinates of any point in the velocity model are represented by t. i This indicates the i-th time sampling point, and the symbol after | indicates the source location corresponding to the forward modeling.

[0069] S22, based on the reciprocity principle, the second simulated wave field ud (x,t i |x R Let u be the direct wave from the underground space point to the detector. d (x R ,t i |x); for the first simulated wavefield u(x,t) i |x S Select The part is the simulated reflected wave field δu(x) j ,t i |x s ).

[0070] Where v(x) is the velocity value at x, and Δt is the seismic sampling rate.

[0071] In detail, in the field of seismic exploration, the reciprocity principle is a fundamental concept that describes a symmetrical property of seismic waves propagating in a medium. Simply put, the reciprocity principle states that in a linear, isotropic medium, if seismic waves are generated and received between two points, the recorded seismic information (such as waveform and arrival time) will be numerically identical regardless of how the positions of the generating and receiving points are exchanged.

[0072] Specifically, suppose in a seismic exploration experiment, a seismic wave is emitted from point A to point B and its response is recorded; according to the principle of reciprocity, if we reverse this and emit a seismic wave from point B to point A and record the response at point A, then the results obtained from the two experiments will be comparable or the same.

[0073] The reciprocity principle can help simplify seismic data acquisition and processing. For example, it can be used to reduce computational load when performing forward modeling of seismic waves. Furthermore, in seismic imaging technology, the reciprocity principle is one of the foundations for constructing inversion algorithms.

[0074] In one embodiment of this application, as shown in Figure 5, the above-mentioned S3 calculates the focusing function of the reflected wave residual response based on the Marchenko equation according to the actual seismic record and forward wavefield data, including:

[0075] S31, based on actual earthquake records R obs (x R ,t i |x S ) and simulated seismic records R in forward wavefield data cal (x R ,t i |x S The residual response of the reflected wave is calculated.

[0076] This step is implemented using actual earthquake records R. obs (x R ,t i |x S Subtract simulated earthquake record R cal (x R ,t i |x S The residual response of the reflected wave is calculated and denoted as R(x). R ,t i |x S ) = R obs (x R ,t i |x S )-R cal (x R ,t i |x S ).

[0077] S32, the reflected wave residual response is used as the reflection response of the Marchenko equation system, and the direct wave u in the forward modeling wavefield data is used as the reflection response. d (x R ,t i |x) is used as the input to the Marchenko equation system, and the upward focusing function f is obtained by iteratively solving the Marchenko equation. + (x R ,t i |x) and downward focusing function f - (x R ,t i |x).

[0078] In detail, a direct wave refers to a seismic wave that is generated from a designated location and arrives directly at the receiver without being reflected.

[0079] The process of solving the Marchenko equation in this step can be referred to existing techniques and will not be described in detail here.

[0080] After obtaining the above upward focusing function f + (x R ,t i |x) and downward focusing function f - (x R ,t i After |x), the upward focusing function f will be applied. + (x R ,t i |x), downward focusing function f - (x R ,t iSubstituting |x) and the reflected wave residual response into Equations 1 and 2, we can obtain the decomposed downwave Green's function G. + (x,t i |x R and the ascending wave Green's function G - (x,t nt -t i |x R ). t i nt and nt are the time corresponding to the i-th time sampling point and the number of seismic data samples, respectively.

[0081] In some implementations, the ascending wave Green's function G + (x,t i |x R As shown in Figure 6, the down-wave Green's function G - (x,t nt -t i |x R As shown in Figure 7. In Figures 6 and 7, the vertical axis represents time, and the horizontal axis represents the spatial offset distance.

[0082] In one embodiment of this application, as shown in FIG8, S3 calculates the gradient of the objective function with respect to the P-wave velocity based on the Green's function and forward wave field data, including:

[0083] S33, for the first simulated wavefield u(x,t) i |x S ) and simulated reflected wave field δu(x j ,t i |x s The central difference scheme is used to calculate the second-order partial derivative with respect to time at each time sampling point.

[0084] In detail, the first simulated wave field u(x,t) i |x S ) and simulated reflected wave field δu(x j ,t i |x s The second-order partial derivative of can be calculated using the following formula:

[0085] Among them, t i nt and nt are the time corresponding to the i-th time sampling point and the number of seismic data samples, respectively, and Δt is the seismic sampling rate.

[0086] S34, using the second-order partial derivatives and Green's function as input, calculates the gradient of the objective function with respect to the P-wave velocity at any spatial point. In this step, G... + (x,t|x R ), G - (x,-t|x R), ü(x,t i |x S ) and δü(x,t i |x S ) is the input.

[0087] In some implementations, the gradient of the objective function with respect to the P-wave velocity at any spatial point is calculated using the second-order partial derivatives and the Green's function as inputs, including:

[0088] The gradient of the objective function with respect to the P-wave velocity at any spatial point is calculated using the following formula (as shown in Figure 9):

[0089] Where x represents any spatial point, g(x) represents the gradient of the objective function with respect to the P-wave velocity at any spatial point, and ü(x,t) i ;x s ) denotes the second derivative of the first simulated wave field, δü(x,t) i ;x s G represents the second derivative of the simulated reflected wave field. + (x,t i ;x R G represents the ascending-wave Green's function of the reflected wave residual response. - (x,t nt -t i ;x R Let denot be the Green's function of the reflected wave residual response, v(x) be the velocity value at x, and t be the velocity value at x. i The sampling time is represented by nt, which represents the number of seismic data samples at time t.

[0090] In one embodiment of this application, as shown in FIG10, step S5, determining whether the termination condition is met based on the actual earthquake record and the update rate model, includes:

[0091] S51, perform forward modeling on the updated velocity model again to obtain the corrected simulated seismic record.

[0092] In this step, forward modeling can be performed using the finite difference method of the scalar wave equation. The corrected simulated seismic record can be expressed as R. cal (x R ,t i |x S )′.

[0093] S52, based on the corrected simulated earthquake record R cal (x R ,t i |x S )′ and actual earthquake record R obs (x R ,t i|x S ), calculate the time correlation function between the simulated seismic source and the detector.

[0094] In some implementations, the time correlation function between the simulated seismic source and the detector is calculated using the following formula:

[0095] Among them, C t (τ,h=|x S -x R |) represents the time correlation function, Δt represents the time sampling interval of the seismic record, and x S Indicates the spatial location of the simulated earthquake source, x R The τ represents the spatial location of the simulated geophone, the τ range is the time range of the seismic record, and nt represents the number of seismic data samples corresponding to time t.

[0096] In detail, the time correlation function reflects the relationship between simulated and actual seismic records at different velocities. If the velocity model is accurate, the correlation function will reach a maximum at τ = 0. If the velocity model is inaccurate, the energy will not focus at τ = 0, and the correlation function will then focus at Δτ(h = |x S -x R |) is the maximum. The realization of reflected wave waveform inversion utilizes this focusing characteristic of the data domain to calculate the velocity model update amount through correlation functions.

[0097] S53, Calculate the objective function based on the time-dependent function.

[0098] In detail, this step includes: calculating the time τ corresponding to the maximum value of the time-dependent function, and assigning τ to Δτ (h = |x S -x R |); According to Δτ(h=|x S -x R |), thus obtaining the following objective function:

[0099] Where h represents the offset between the shot point and the receiver point, and Δτ represents the focusing error, i.e., the correlation function C. t (τ,h=|x S -x R When |) reaches its maximum value, the time shift of the reflected focusing energy from τ=0 is represented by n, where n represents the iteration number, and J n This represents the objective function.

[0100] In detail, the maximum value of the time correlation function indicates that the simulated earthquake record is most similar to the actual earthquake record at that point. The corresponding time value represents the time required to shift the simulated earthquake record by τ to achieve the greatest similarity to the actual earthquake record. Therefore, τ measures the difference in travel time between the simulated and actual earthquake records. τ is assigned to Δτ(h=|x S -x R |), thereby calculating the travel time difference between all simulated and actual earthquake records.

[0101] The objective function is based on the mean square error that minimizes the difference in travel time between simulated and actual seismic records, reflecting the magnitude of the difference in travel time between simulated and actual seismic records caused by velocity errors.

[0102] S54. Determine whether the objective function meets the preset precision. If yes, the termination condition is met; otherwise, the termination condition is not met.

[0103] When implementing this step, the preset precision is, for example, the current objective function value J. n The value is less than a preset multiple of the initial iteration objective function J1, for example, J1*0.001. The objective function is then checked to determine if it meets the accuracy requirement, i.e., whether it is less than J1*0.001. If yes, the termination condition is met; otherwise, the preset condition is not met, and the updated velocity model is used as the new initial velocity model. Steps S2 to S5 are repeated until the accuracy requirement is met.

[0104] In one embodiment, the actual velocity model is shown in Figure 2, the velocity model established using the method of this application is shown in Figure 11, and the velocity model established using the conventional waveform inversion method is shown in Figure 12. By comparing Figures 2, 11, and 12, it can be seen that this application can more accurately estimate the reflected wave data by using the focusing function, and the velocity model established using the method of this application is closer to the actual velocity model.

[0105] Based on the same inventive concept, this application also provides a depth velocity modeling device based on reflection focusing, as described in the following embodiments. Since the principle and method of solving the problem by the depth velocity modeling device based on reflection focusing are similar, the implementation of the depth velocity modeling device based on reflection focusing can refer to the depth velocity modeling method based on reflection focusing, and repeated details will not be elaborated further. Specifically, as shown in Figure 13, the depth velocity modeling device based on reflection focusing includes:

[0106] Modeling unit 1301 is used to determine the actual earthquake record and the initial velocity model;

[0107] Forward modeling unit 1302 is used to perform forward modeling of the initial velocity model based on the seismic observation system to obtain forward wavefield data;

[0108] The gradient calculation unit 1303 is used to calculate the focusing function of the reflected wave residual response based on the Marchenko equation, according to the actual seismic record and forward wavefield data; calculate the Green's function of the reflected wave residual response based on the focusing function; and calculate the gradient of the objective function with respect to the P-wave velocity based on the Green's function and forward wavefield data.

[0109] The correction unit 1304 is used to determine the updated velocity model based on the gradient of the objective function with respect to the P-wave velocity and the update step size.

[0110] Judgment unit 1305 is used to determine whether the termination condition is met based on the actual seismic record and the updated velocity model. If yes, the final updated velocity model is used as the depth velocity model; otherwise, the forward modeling unit, gradient calculation unit, and correction unit are restarted.

[0111] This application can use the Green's function to replace the reflected wave response wave field obtained by the migration / inverse migration algorithm, without being constrained by the accuracy of the background velocity, and without the need for additional migration / inverse migration processing with huge computational cost, thereby improving the determination accuracy and efficiency of the velocity model.

[0112] In one embodiment of this application, as shown in FIG14, computer device 1402 may include one or more processors 1404, such as one or more central processing units (CPUs), each of which may implement one or more hardware threads. Computer device 1402 may also include any memory 1406 for storing information of any kind, such as code, settings, data, etc. Non-limitingly, for example, memory 1406 may include any type of RAM, any type of ROM, flash memory device, hard disk, optical disk, etc. More generally, any memory can use any technology to store information. Further, any memory may provide volatile or non-volatile retention of information. Further, any memory may represent a fixed or removable component of computer device 1402. In one case, when processor 1404 executes associated instructions stored in any memory or combination of memories, computer device 1402 may perform any operation of the associated instructions. Computer device 1402 also includes one or more drive mechanisms 1408 for interacting with any memory, such as hard disk drive mechanisms, optical disk drive mechanisms, etc.

[0113] Computer device 1402 may also include an input / output module 1410 (I / O) for receiving various inputs (via input device 1412) and providing various outputs (via output device 1414). A specific output mechanism may include a presentation device 1416 and an associated graphical user interface 1418 (GUI). In other embodiments, the input / output module 1410 (I / O), input device 1412, and output device 1414 may be omitted, and the device may function solely as a computer device within a network. Computer device 1402 may also include one or more network interfaces 1420 for exchanging data with other devices via one or more communication links 1422. One or more communication buses 1424 couple the components described above together.

[0114] Communication link 1422 can be implemented in any way, such as via a local area network, a wide area network (e.g., the Internet), a point-to-point connection, or any combination thereof. Communication link 1422 may include any combination of hardwired links, wireless links, routers, gateway functions, name servers, etc., governed by any protocol or combination of protocols.

[0115] This application also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, performs the steps of the above-described method.

[0116] This application also provides computer-readable instructions, wherein when a processor executes the instructions, the program therein causes the processor to perform the method of any of the foregoing embodiments.

[0117] It should be understood that in the various embodiments of this application, the order of the above-mentioned processes does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0118] It should also be understood that, in the embodiments of this application, the term "and / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, and B existing alone. Additionally, the character " / " in this application generally indicates that the preceding and following related objects have an "or" relationship.

[0119] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed in this application can be implemented in electronic hardware, computer software, or a combination of both. To clearly illustrate the interchangeability of hardware and software, the composition and steps of each example have been generally described in terms of functionality in the foregoing description. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0120] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0121] In the embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the couplings or direct couplings or communication connections shown or discussed may be indirect couplings or communication connections through some interfaces, apparatuses, or units, or they may be electrical, mechanical, or other forms of connection.

[0122] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of the embodiments of this application, depending on actual needs.

[0123] Furthermore, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit.

[0124] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0125] This application uses specific embodiments to illustrate the principles and implementation methods of this application. The description of the above embodiments is only for the purpose of helping to understand the method and core ideas of this application. At the same time, for those skilled in the art, there will be changes in the specific implementation methods and application scope based on the ideas of this invention. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A depth velocity modeling method based on reflection focusing, characterized in that, include: Determine the actual earthquake record and initial velocity model; Based on the earthquake observation system, forward modeling of the initial velocity model is performed to obtain forward wavefield data; Based on the actual seismic records and the forward wavefield data, the focusing function of the reflected wave residual response is calculated using the Marchenko equation; the Green's function of the reflected wave residual response is calculated based on the focusing function; and the gradient of the objective function with respect to the P-wave velocity is calculated based on the Green's function and the forward wavefield data. The update velocity model is determined based on the gradient of the objective function with respect to the P-wave velocity and the update step size; Based on the actual seismic records and the updated velocity model, determine whether the termination condition is met. If yes, use the final updated velocity model as the depth velocity model; otherwise, re-execute the forward modeling and subsequent steps.

2. The method as described in claim 1, characterized in that, Based on the seismic observation system, forward modeling of the initial velocity model is performed to obtain forward wavefield data, including: Based on the earthquake observation system, a forward modeling simulation was performed on the initial velocity model to obtain the first simulated wave field u(x,t) generated at any earthquake source. i |x S The second simulated wave field u excited at the detector d (x,t i |x R And extract the simulated seismic record R at the detector. cal (x R ,t i |x S ); Based on the reciprocity principle, the second simulated wave field u d (x,t i |x R Let u be the direct wave from the underground space point to the detector. d (x R ,t i |x); for the first simulated wavefield u(x,t) i |x S Select The part is the simulated reflected wave field δu(x) j ,t i |x s ); Where, x S The x-coordinate represents the spatial location coordinates of the earthquake source. R The coordinates of the detector's spatial position are represented by x, and the coordinates of any point in the velocity model are represented by t. i Let | represent the i-th time sampling point, the symbol after | is the source location corresponding to the forward modeling, v(x) is the velocity value at x, and Δt is the seismic sampling rate.

3. The method as described in claim 2, characterized in that, Based on the actual seismic records and the forward wavefield data, the focusing function of the reflected wave residual response is calculated using the Marchenko equation, including: Based on the actual earthquake records and the simulated earthquake records R in the forward wavefield data cal (x R ,t i |x S The residual response of the reflected wave was calculated. The reflected wave residual response is used as the reflection response of the Marchenko equation system, and the direct wave u in the forward modeling wavefield data is used as the reflection response. d (x R ,t i |x) is used as the input to the Marchenko equation system, and the upward focusing function f is obtained by iteratively solving the Marchenko equation. + (x R ,t i |x) and downward focusing function f - (x R ,t i |x).

4. The method as described in claim 2, characterized in that, Calculating the gradient of the objective function with respect to the P-wave velocity based on the Green's function and the forward wavefield data includes: For the first simulated wave field u(x,t) i |x S ) and simulated reflected wave field δu(x j ,t i |x s Using the central difference scheme, the second-order partial derivative with respect to time at each time sampling point is calculated; Using the second-order partial derivatives and Green's function as inputs, the gradient of the objective function with respect to the longitudinal wave velocity at any spatial point is calculated.

5. The method as described in claim 4, characterized in that, Using the second-order partial derivatives and Green's function as input, the gradient of the objective function with respect to the P-wave velocity at any spatial point is calculated, including: The gradient of the objective function with respect to the P-wave velocity at any spatial point is calculated using the following formula: Where x represents any spatial point, g(x) represents the gradient of the objective function with respect to the P-wave velocity at any spatial point, and ü(x,t) i ;x s ) denotes the second derivative of the first simulated wave field, δü(x,t) i ;x s G represents the second derivative of the simulated reflected wave field. + (x,t i ;x R G represents the ascending-wave Green's function of the reflected wave residual response. - (x,t nt -t i ;x R Let denot be the Green's function of the reflected wave residual response, v(x) be the velocity value at x, and t be the velocity value at x. i The sampling time is represented by nt, which represents the number of seismic data samples at time t.

6. The method as described in claim 1, characterized in that, The updated velocity model is determined based on the gradient of the objective function with respect to the P-wave velocity and the update step size, including: The velocity model is updated using the following formula: v(x) n+1 =v(x) n -α n ·g(x) n ; Where n represents the number of iterations, v(x) n+1 Let v(x) represent the model for the (n+1)th iteration. n Let α represent the velocity model for the nth iteration. n Let g(x) represent the step size. n This represents the gradient of the objective function with respect to the P-wave velocity at the nth iteration.

7. The method as described in claim 2, characterized in that, Determine whether the termination conditions are met based on the actual earthquake records and update rate model, including: The updated velocity model was re-performed in forward modeling to obtain the corrected simulated seismic record; Based on the corrected simulated earthquake record R cal (x R ,t i |x S )′ and the actual earthquake record R obs (x R ,t i |x S ), calculate the time correlation function between the simulated seismic source and the detector; Calculate the objective function based on the time-related function; Determine whether the objective function meets the preset precision. If it does, the termination condition is met; otherwise, the termination condition is not met.

8. The method as described in claim 7, characterized in that, Based on the corrected simulated earthquake record R cal (x R ,t i |x S )′ and the actual earthquake record R obs (x R ,t i |x S ), calculate the time correlation function between the simulated seismic source and the detector, including: The time correlation function between the simulated seismic source and the detector is calculated using the following formula: Among them, C t (τ,h=|x S -x R |) represents the time correlation function, Δt represents the seismic record time sampling interval, and xS represents the spatial location of the simulated seismic source. R The τ represents the spatial location of the simulated geophone, the τ range is the time range of the seismic record, and nt represents the number of seismic data samples corresponding to time t.

9. The method as described in claim 8, characterized in that, Based on the time-dependent function, calculate the objective function, including: Find the time τ corresponding to the maximum value of the time correlation function, and assign τ to Δτ(h=|x S -x R |); According to Δτ(h=|x S -x R |), thus obtaining the following objective function: Where h represents the offset between the shot point and the receiver point, and Δτ represents the focusing error, i.e., the correlation function C. t (τ,h=|x S -x R When |) reaches its maximum value, the time shift of the reflected focusing energy from τ=0 is represented by n, where n represents the iteration number, and J n This represents the objective function.

10. A depth velocity modeling device based on reflection focusing, characterized in that, include: Modeling units are used to determine the actual seismic record and the initial velocity model; The forward modeling unit is used to perform forward modeling of the initial velocity model based on the seismic observation system to obtain forward wavefield data; The gradient calculation unit is used to calculate the focusing function of the reflected wave residual response based on the Marchenko equation according to the actual seismic record and the forward wave modeling data; calculate the Green's function of the reflected wave residual response according to the focusing function; and calculate the gradient of the objective function with respect to the P-wave velocity according to the Green's function and the forward wave modeling data. The correction unit is used to determine the updated velocity model based on the gradient of the objective function with respect to the P-wave velocity and the update step size; The judgment unit is used to determine whether the termination condition is met based on the actual seismic record and the updated velocity model. If yes, the final updated velocity model is used as the depth velocity model; otherwise, the forward modeling unit, gradient calculation unit, and correction unit are restarted.

11. A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the method according to any one of claims 1 to 9.

12. A computer storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor of the computer device, it implements the method of any one of claims 1 to 9.

13. A computer program product, characterized in that, The computer program product includes a computer program that, when executed by a processor of a computer device, implements the method of any one of claims 1 to 9.

Images