A three-dimensional reverse-time migration method for petroleum exploration

Abstract

The present application relates to the technical field of petroleum geophysical exploration, and particularly relates to a three-dimensional reverse time migration method for petroleum exploration, comprising the following steps: S1, forward simulation is performed by using a three-dimensional velocity model for forward numerical simulation, and the forward simulation adopts a finite difference method to solve a three-dimensional acoustic wave equation; S2, three-dimensional acoustic wave equation difference decomposition stability condition parameters are determined according to a first-order difference space dispersion relation and stability conditions of a three-dimensional staggered grid format in the finite difference method; S3, the parameter three-dimensional acoustic wave equation difference decomposition stability condition parameters are input into an objective function of optimized difference coefficients to obtain the optimized difference coefficients; and S4, numerical simulation is performed on the three-dimensional acoustic wave equation of the three-dimensional velocity model by using the optimized difference coefficients to obtain a three-dimensional reverse time migration result. When the time sampling step is increased during calculation, the time sampling points can be reduced by several times under the condition that the total extrapolation time is constant, and the calculation time cost and the wave field storage space can be reduced by several times.

Description

Technical Field

[0001] This invention relates to the field of petroleum geophysical exploration technology, and in particular to a three-dimensional reverse time migration method for petroleum exploration. Background Technology

[0002] In recent years, exploration targets have become increasingly complex and concealed, and the requirements for horizontal drilling depth control have become increasingly stringent, as have the requirements for seismic exploration accuracy. Achieving accurate imaging of complex underground geological bodies is of great importance, and the industry is paying more and more attention to the application of the reverse time migration method, which has the advantages of being unrestricted by angle and providing high-precision imaging.

[0003] Since reverse-time migration requires calculating both forward and reverse extended wavefields and storing the forward wavefield, the computational efficiency and memory consumption of wavefield extension are particularly critical. Currently, two-dimensional reverse-time migration has a relatively small computational load, and most computers can meet the memory requirements; however, three-dimensional reverse-time migration will increase the computational load by tens or even hundreds of times compared to two-dimensional migration. A common feasible solution to this problem is to increase the time sampling step size, which can reduce the number of time sampling points by a factor of two for a given total extension time, thereby reducing computational time cost and wavefield storage space by a factor of two. However, due to the limitation of the number of finite difference approximations (CFLs), the stability factor of the difference coefficients solved by conventional methods is not high, and large time sampling intervals cannot be used for wavefield extension. To overcome the limitations of conventional CFLs, it is necessary to improve the stability factor of the difference approximation.

[0004] Currently, methods to improve the stability of difference approximations mainly fall into two categories. The first is to improve the stability of the algorithm by increasing the accuracy of the time-domain difference, such as the Lax-Wendroff (LW) method and the FD method based on time-domain dispersion relations (Liu and Zhang, 2019; Chen, 2007). Unlike traditional methods that obtain spatial FD coefficients through spatial domain dispersion relations, FD methods based on spatiotemporal domain dispersion relations calculate spatial FD coefficients through spatiotemporal domain dispersion relations (Liu, 2009; Liu, 2013), which can improve the accuracy and stability of the algorithm, but the improvement in stability is limited. Taking two-dimensional acoustic wave modeling as an example, the maximum CFL number of the stable model is approximately 0.707. The second method is to remove unstable eigenvalues ​​or perturb them to stabilize the algorithm. Gaffar, by analyzing the root causes of instability in explicit time-domain FD methods for electromagnetic waves, directly derived unstable modes with large eigenvalues. Subsequently, Gao applied the so-called eigenvalue perturbation to acoustic wave FD modeling to remove the stability limit, but this method requires a considerable amount of memory. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of conventional methods in the prior art, such as low stability factors of the differential coefficients and inability to use large time sampling intervals for wavefield extension. This invention provides a three-dimensional reverse time migration method for oil exploration, which achieves a significant reduction in time sampling points, computation time costs, and wavefield storage space when the total extension time is constant.

[0006] To achieve the above-mentioned objectives, the present invention provides the following technical solution:

[0007] A three-dimensional reverse time migration method for oil exploration, the method comprising:

[0008] S1 uses the three-dimensional velocity model used in forward numerical simulation to perform forward simulation, and the forward simulation uses the finite difference method to solve the three-dimensional acoustic wave equation.

[0009] S2 determines the stability condition parameters of the differential equation of the three-dimensional acoustic wave based on the first-order difference space dispersion relation and stability condition of the three-dimensional staggered mesh scheme in the finite difference method;

[0010] S3 inputs the stability condition parameters of the three-dimensional acoustic wave equation difference into the objective function of the optimization difference coefficients to obtain the optimized difference coefficients;

[0011] S4 uses the optimized difference coefficients to perform a numerical simulation of the three-dimensional acoustic wave equation of the three-dimensional velocity model, and obtains the three-dimensional reverse time migration result.

[0012] Furthermore, the first-order difference space dispersion relation of the three-dimensional staggered mesh includes:

[0013]

[0014]

[0015] Where β is half the product of the total wavenumber and the spatial step size; c m denoted as difference coefficients; m represents the finite difference order; M represents the maximum finite difference order.

[0016] Furthermore, the stability conditions of the three-dimensional staggered mesh include:

[0017] Where C represents the CFL number, and ψ represents the function value of the spatial partial derivative operator in the frequency domain at the Nyqusit wavenumber.

[0018] Specifically, the stability conditions of the three-dimensional staggered mesh are obtained according to the von Neumann stability analysis method, including:

[0019]

[0020]

[0021]

[0022] This represents the frequency domain function obtained after Fourier transforming the spatial partial derivative difference operator with 2M-order precision; β x β represents half the product of the wavenumber in the x-direction and the spatial step size; y β represents half the product of the y-direction wavenumber and the spatial step size; z Δt represents half of the product of the wavenumber in the z-direction and the spatial step size; h represents the spatial step size of the difference operator.

[0023] Specifically, the objective function for optimizing the difference coefficients includes:

[0024]

[0025] Where, β i It is the maximum value that ψ(βi) approximates.

[0026] Specifically, the stability condition parameters of the three-dimensional acoustic wave equation differential decomposition include parameters B, Δβ, M, and η; where Δβ represents the size of the transition interval, η represents the maximum permissible error, B represents the effective frequency band coverage under the control of the maximum permissible error η, and M represents the maximum finite difference order.

[0027] It should be noted that β in the objective function i The range is [0,π].

[0028] Furthermore, β i The value range is [0, π]. The value range is further divided into three intervals: the effective interval [0, B], the transition interval [B, B+Δβ], and the control interval [B+Δβ, π].

[0029] A system for a three-dimensional reverse time migration method for oil exploration, the system storing a computer program that, when executed by a processor, implements the aforementioned three-dimensional reverse time migration method for oil exploration.

[0030] An apparatus based on a three-dimensional reverse time migration method system for oil exploration includes at least one processor and a memory communicatively connected to the at least one processor; the memory stores instructions executable by the at least one processor, which, when executed by the at least one processor, enable the at least one processor to perform the aforementioned three-dimensional reverse time migration method for oil exploration.

[0031] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0032] By utilizing the method proposed in this patent that significantly improves stability, the maximum number of CFLs satisfying the stability condition is greatly increased during forward numerical simulation. Therefore, the time sampling step can be increased during calculation, which can multiply the number of time sampling points when the total extension time is constant, thereby multiplying the computation time cost and wavefield storage space. Attached image description:

[0033] Figure 1 This is a schematic diagram of the stability factor of the conventional difference approximation;

[0034] Figure 2 This is a schematic diagram illustrating the variation of the stability factor with B;

[0035] Figure 3 This is a flowchart illustrating the method provided by the present invention;

[0036] Figure 4 This is a schematic diagram comparing wavefield snapshots of three-dimensional uniform models using different difference coefficients in an embodiment of the present invention;

[0037] Figure 5 This is a schematic diagram illustrating the test results of the present invention applied to some three-dimensional reverse engineering models. Detailed Implementation

[0038] The present invention will be further described in detail below with reference to experimental examples and specific embodiments. However, this should not be construed as limiting the scope of the above-mentioned subject matter of the present invention to the following embodiments; all technologies implemented based on the content of the present invention fall within the scope of the present invention.

[0039] This invention provides a three-dimensional reverse time migration method for oil exploration that breaks through the limitation of the number of CFLs in the staggered grid differential approximation. It can reduce the number of time sampling points by a factor of two when the total extension time is constant, thereby reducing the computation time cost and wavefield storage space by a factor of two.

[0040] Example 1

[0041] This embodiment provides a three-dimensional reverse time migration method for oil exploration, the method comprising:

[0042] S1 uses the three-dimensional velocity model used in forward numerical simulation to perform forward simulation, and the forward simulation uses the finite difference method to solve the three-dimensional acoustic wave equation.

[0043] S2 determines the stability condition parameters of the differential equation of the three-dimensional acoustic wave based on the first-order difference space dispersion relation and stability condition of the three-dimensional staggered mesh scheme in the finite difference method;

[0044] S3 inputs the stability condition parameters of the three-dimensional acoustic wave equation difference into the objective function of the optimization difference coefficients to obtain the optimized difference coefficients;

[0045] S4 uses the optimized difference coefficients to perform a numerical simulation of the three-dimensional acoustic wave equation of the three-dimensional velocity model, and obtains the three-dimensional reverse time migration result.

[0046] It should be noted that the stability condition parameters of the three-dimensional acoustic wave equation differential decomposition include parameters B, Δβ, M and η; where Δβ represents the size of the transition interval, η represents the maximum permissible error, B represents the effective frequency band coverage under the control of the maximum permissible error η, and M represents the maximum finite difference order.

[0047] The stability conditions of the three-dimensional staggered mesh are obtained based on the von Neumann stability analysis method. By analyzing the shortcomings of the traditional method that only constrains the stability conditions in the small wavenumber interval, a method for constraining stability in a larger wavenumber interval is proposed.

[0048] The stability condition for the three-dimensional staggered mesh obtained by the Neumann stability analysis method is:

[0049]

[0050] Based on the dispersion relation of the first-order difference space of the staggered grid, we have:

[0051] β is related to the wave number and the spatial frequency of the wave. k represents the wave number, h represents the spatial step size of the difference operator; c m The coefficients are represented as difference coefficients; m represents the finite difference order; M represents the maximum finite difference order.

[0052] C is the CFL number (Courant-Friedrichs-Lewy Condition). C is directly proportional to the wave propagation speed in the medium and the time step of the difference operator, and inversely proportional to the spatial step of the difference operator. The expression for C is:

[0053] There is a function ψ related to β.

[0054] Substituting the function ψ into the stability conditions of the three-dimensional staggered mesh obtained by the von Neumann stability analysis method, we get the following conditions:

[0055]

[0056] The conventional optimization finite difference coefficient relationship is as follows:

[0057]

[0058] Where B is the effective frequency band coverage under the control of the maximum permissible error ε.

[0059] By using the objective function in the above formula to solve for the difference coefficients, we can obtain β. i The difference coefficients within the range [0, B] have a maximum permissible error of ε = 1e-4. For example... Figure 1 As shown, LS represents the conventional optimized finite difference coefficients, and TEA represents the traditional Taylor (TE) expansion method. As the finite difference order M increases, the stability factor S of LS is always less than that of TEA. It can be seen that the stability factor S of the above objective function is lower than that of the traditional Taylor (TE) expansion method.

[0060] Stability is optimized only within the wavenumber interval β∈[0, B] while maintaining algorithm accuracy. However, no stability constraints are applied to the larger wavenumber interval β∈[B, π], where ψ will vary with β. The increase is almost proportional to the increase in C, thus leading to max Decreasing the value reduces stability.

[0061] To obtain a higher C, this patent attempts to adjust β in the objective function. i The range is extended to [0,π], and the maximum permissible error is appropriately reduced accordingly, resulting in the following new objective function:

[0062]

[0063]

[0064] In the formula f(β) i ) is about β i The function, because ψ(β) i ) is a smooth function, while Φ(β) is a smooth function. i Since β is not a smooth function, it will produce a large approximation error, especially near β = B.

[0065] Determine the transition function: In order to reduce Φ(β) i The error caused by non-smoothness needs to be addressed by establishing an interval [B, B+Δβ] such that Φ(β) i f(β) smoothly transforms from [0,B] to [B+Δβ,π] i Combining the two equations above, the new objective function can be rewritten as:

[0066]

[0067] Because of β i ∈[0,B],β i It is ψ(β) i To maximize the stability of the difference approximation, the maximum value of ψ(β) is approximated. i ) throughout β iThe maximum approximation value within the interval [0,π] is B, then f(β) i )=β i The final form of the new objective function is:

[0068]

[0069] Where β i The value range of is [0, π]. This method can further divide the wavenumber interval into three intervals: the effective interval [0, B], the transition interval [B, B+Δβ], and the control interval [B+Δβ, π].

[0070] Compared to traditional methods that optimize difference coefficients and control stability only within the range [0, B], the method proposed in this embodiment not only inherits the advantages of traditional methods in optimizing difference coefficients and improving stability within the effective range [0, B], but also, by applying the aforementioned constraint control to the function ψ in a larger wavenumber range and the control range [B+Δβ, π], it enables stability optimization across the entire wavenumber range [0, π]. This ensures that the value of the function ψ does not increase with the increase of β in the large wavenumber interval [B+Δβ,π], thereby improving the stability of the traditional control method.

[0071] like Figure 2 As shown, under the conditions of M=8 and ε=1e-3, the stability factor S obtained by solving the above objective function using the least squares method varies with B. It can be seen that LS-Old represents the traditional optimization objective function and LS-New represents the new optimization objective function. As B decreases, the stability factor S of the new optimization objective function gradually becomes greater than that of the traditional optimization objective function. Therefore, under the condition of acceptable error, the stability of the new optimization objective function is greatly improved at medium and low wavenumbers.

[0072] like Figure 4 The image shows wave field snapshots obtained by three-dimensional numerical simulation using the staggered grid finite difference method with time steps of 1ms, 2ms, and 4ms, respectively, using the traditional method and the new method proposed in this patent. The three-dimensional uniform model is set to a size of 4km×4km×4km, the P-wave velocity is 4000m / s, and the spatial step size is h=10m. The wave field snapshots obtained are obtained by the traditional objective function and the new objective function with optimized difference coefficients using the finite difference method.

[0073] It can be seen that the new method is stable in wave field simulation when the time sampling interval Δt = 4ms, while the old method is already unstable when Δt = 2ms. This proves that the method proposed in this patent can effectively improve the stability of solving the three-dimensional acoustic wave equation by the finite difference method.

[0074] like Figure 5As shown, a model was constructed with a size of 3km×3km×2km, a maximum P-wave velocity of 4500m / s, a spatial step size of h=10m, a time sampling interval of Δt=4ms, and optimized parameters of M=8, ε=1E-3, and B=0.3.

[0075] from Figure 5 Figure a shows the three-dimensional velocity model used in the forward numerical simulation, and Figure b shows the wavefield snapshot obtained from the forward numerical simulation. In Figures a and b, the X and Y axes represent the distance (km) along the X and Y survey lines, respectively, and the Z axis represents the depth (km). Figure c shows the common shot gather, where the X and Y axes represent the shot-receiver distance (km) along the X and Y survey lines, respectively, and the Z axis represents the depth (km). Figure d shows the migration profile, where the X and Y axes represent the point numbers along the X and Y survey lines, respectively, with a step size of 10m, and the Z axis represents the depth (km). From the wavefield snapshot and shot gather, the wavefield snapshot is clear, and the direct and reflected wave characteristics are obvious in the shot gather record. The solid-difference algorithm is stable, and the numerical simulation results are correct. Using a partial three-dimensional inverse model for finite-difference forward numerical simulation and reverse-time migration, the resulting migration profile has a high degree of stratigraphic matching with the initial velocity model, further verifying the reliability of the new method proposed in this patent.

[0076] Large time sampling step size and large CFL number improve the stability of the algorithm. With a fixed total extension time, the number of time sampling points is reduced by a factor of two, resulting in a significant reduction in computation time cost and wavefield storage space. This convenience can save a lot of cost for subsequent reverse time migration.

[0077] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A three-dimensional reverse-time migration method for oil exploration, characterized in that, The method includes: S1 uses the three-dimensional velocity model used in forward numerical simulation to perform forward simulation, and the forward simulation uses the finite difference method to solve the three-dimensional acoustic wave equation. S2 determines the stability condition parameters of the differential equation of the three-dimensional acoustic wave based on the first-order difference space dispersion relation and stability condition of the three-dimensional staggered mesh scheme in the finite difference method; The first-order difference space dispersion relation of the three-dimensional staggered mesh format includes: in, It is half the product of the total wavenumber and the spatial step size. To and The relevant spatial dispersion relation function, For the order of the maximum finite difference, For finite difference order, These are the difference coefficients. Total wavenumber; Let S be the spatial step size of the difference operator; The stability condition of the three-dimensional staggered mesh is obtained according to the von Neumann stability analysis method, which includes: in, express The frequency domain function obtained by Fourier transforming the spatial partial derivative difference operator with first-order precision; express Half of the product of directional wave number and spatial step size; express Half of the product of directional wave number and spatial step size; express Half of the product of directional wave number and spatial step size; Indicates the longitudinal wave velocity of the medium; Indicates the sampling time interval; This represents the spatial step size of the difference operator. Let be the function value of the spatial partial derivative operator in the frequency domain at the Nyquist wavenumber; The stability conditions for the three-dimensional staggered mesh include: in, For CFL numbers, The maximum value of the CFL number; The stability condition parameters of the difference decomposition of the three-dimensional acoustic equation include parameters. , , and ; in, Indicates the size of the transition interval. Indicates the maximum permissible error. Indicates the maximum permissible error Effective frequency band coverage under controlled conditions Represents the order of the maximum finite difference; S3 inputs the stability condition parameters of the differential equation of the three-dimensional acoustic wave into the objective function of the differential coefficients to obtain the optimized differential coefficients; The objective function for optimizing the difference coefficients includes: in, for The approximate maximum value, To optimize the objective function of the difference coefficients, To and The relevant spatial dispersion relation function, For sampling sequence number, This represents the number of sampling points; S4 uses the optimized difference coefficients to perform a numerical simulation of the three-dimensional acoustic wave equation of the three-dimensional velocity model, and obtains the three-dimensional reverse time migration result.

2. The three-dimensional reverse-time migration method for oil exploration according to claim 1, characterized in that, The range of values ​​is The value range is further divided into three intervals, namely the valid interval. transition zone and control range .

3. A system for a three-dimensional reverse-time migration method for oil exploration, characterized in that, The system stores a computer program that, when executed by a processor, implements the method as described in any one of claims 1 to 2.

4. An apparatus based on a three-dimensional reverse-time migration method for oil exploration, characterized in that, It includes at least one processor and a memory communicatively connected to the at least one processor; the memory stores instructions executable by the at least one processor, which, when executed by the at least one processor, enable the at least one processor to perform the three-dimensional reverse time migration method for oil exploration as described in any one of claims 1 to 2.

Images