A method, device and readable storage medium for generating amplitude-preserving gathers based on least squares reverse time migration

Abstract

The present invention discloses a method, device and readable storage medium for generating amplitude-preserving angle gathers based on least squares reverse time migration. Among them, a method for generating amplitude-preserving angle gathers based on least squares reverse time migration includes extracting angle gathers as initial imaging results in reverse time migration, using the initial imaging results as input, obtaining simulated seismic data using a wave equation forward operator based on Kirchhoff approximation, substituting the residuals between the simulated seismic data and the observed data into the reverse time migration imaging operator to obtain a gradient image, and using the conjugate gradient method to iteratively solve the least squares inversion problem to optimize the angle gathers. This method allows the imaging amplitude to represent the reflection coefficient related to the angle when using LSRTM for quantitative interpretation of complex structures, that is, it is presented in the form of imaging gathers, which can be directly used in seismic interpretation of amplitude-varying-angle (AVA). The use of the LSRTM inversion imaging method can make the extracted angle gathers more amplitude-preserving, and can improve the accuracy and reliability of AVA analysis.

Description

Technical Field

[0001] The present invention relates to the field of seismic exploration technology, and in particular to a method, device and readable storage medium for generating amplitude-preserving gathers based on least squares reverse time migration. Background Art

[0002] Seismic depth migration aims to generate accurate images of subsurface structures from observed seismic reflection data. During this process, the amplitude fidelity of the seismic images is expected to be preserved so that the amplitude information in the images can be used to invert and interpret the physical properties of the subsurface medium. However, in practice, seismic images exhibit uneven illumination, influenced by the observation system and velocity variations in the overlying medium. In particular, illumination variations can be dramatic under complex overlying media, such as in structures beneath salt domes. In these cases, full-wave equation migration methods, such as reverse-time migration, are necessary. To obtain physically meaningful amplitudes, one approach is to model seismic migration as an inverse problem. Assuming a known background velocity model, the goal is to find a reflection coefficient model that minimizes the difference between the simulated and observed seismic data. This imaging method, which solves linear inverse problems, is called least-squares reverse-time migration (LSRTM).

[0003] In recent years, the development of LSRTM has greatly improved the quality of stacked imaging of complex structures. Existing LSRTM methods generally assume that the subsurface reflectivity is independent of the reflection angle. However, the reflection coefficient, which is related to the elastic properties of the subsurface medium, depends on the reflection angle and sometimes also on the azimuth. The amplitude in the stacked image is independent of the angle and is a combination of the reflection coefficient amplitudes at different observation angles. Therefore, the amplitude in the stacked image does not have a clear physical meaning, resulting in the amplitude versus angle analysis (AVA) problem of not preserving the amplitude, thereby reducing the analysis accuracy and reliability of AVA. Summary of the Invention

[0004] Based on this, it is necessary to provide a method, device and readable storage medium for generating amplitude-preserving angle gathers based on least squares reverse time migration to solve the problem that the AVA technology in the existing technology does not preserve amplitude, thereby reducing the analysis accuracy and reliability of AVA.

[0005] The present invention provides a method for generating amplitude-preserving gathers based on least squares reverse time migration, comprising the following steps:

[0006] Extracting angle gathers as initial imaging results in reverse time migration;

[0007] Taking the initial imaging results as input, the simulated seismic data are obtained by using the wave equation forward operator based on Kirchhoff approximation.

[0008] Substituting the residual between the simulated seismic data and the observed data into a reverse time migration imaging operator to obtain a gradient image;

[0009] The conjugate gradient method is used to iteratively solve the least squares inversion problem to optimize the angle gathers.

[0010] Furthermore, the original common shot point seismic records are collected;

[0011] The collected original common shot point seismic records are preprocessed, and the common shot seismic records obtained after processing are called observed seismic data, which are recorded as D(x r ,x s ; t), where x r Indicates the coordinates of the detection point, x s represents the coordinates of the earthquake source, and t represents the time;

[0012] Based on velocity analysis, the migration velocity v(x) is constructed, where x = (x, z) is the spatial coordinate, and the density is not considered, assuming that the density ρ≡1;

[0013] By observing the seismic data D(x r ,x s ; t) spectrum analysis, constructing a wide-band offset wavelet function w(t);

[0014] The forward wave field extension is performed by using the migration wavelet w(t) as the source term of the constant density acoustic wave equation to obtain the source wave field p F (x; t);

[0015] Taking the observed earthquake data D(x r ,x s t) as the boundary condition of the constant density acoustic wave equation to perform reverse time wave field extension and obtain the received wave field p B (x; t);

[0016] For the source wave field p F (x; t) and the received wave field p B (x; t) respectively use the Poynting vector method to estimate the direction vector P of the wave field propagation direction s (x; t) and P r (x; t), according to P s (x; t) and P r The angle between (x and t) is used to calculate the reflection angle θ, and the P s (x; t) and P r The normal vector of the sum vector P(x; t) of (x; t) estimates the formation dip δ;

[0017] For the source wave field p F (x; t) and the received wave field p B(x; t) applying the cross-correlation imaging condition and performing angle binning to distribute the imaging results to the corresponding reflection angles θ to obtain the angle domain common imaging point gathers R(x; θ), which are the initial imaging results.

[0018] The forward operator based on Kirchhoff approximation is used to As the boundary condition, the forward wave field is extended and the seismic wave field is recorded at the detection point to obtain the simulated seismic data d(x r ,x s ;t);

[0019] According to the adjoint state method, the residual d between the simulated seismic data and the observed seismic data is r (x r ,x s ; t) = d(x r ,x s ; t)-D(x r ,x s t) is substituted into the step of generating angle gathers in migration imaging as the boundary condition to obtain the least squares objective functional corresponding to The gradient image ΔR(x,θ)=L T (LR(x,θ)-D);

[0020] Where L represents the Kirchhoff forward operator, L T represents the migration imaging operator for generating angle gathers, and D represents the observation data;

[0021] Iteratively update the diagonal gather using the conjugate gradient method R n (x,θ)=R n-1 (x,θ)+αΔR g (x,θ), α is the update step size calculated by the conjugate gradient method, ΔR g (x,θ) is the conjugate gradient direction.

[0022] Furthermore, the observation of seismic data D(x r ,x s t), the steps of constructing a broadband offset wavelet function w(t) include:

[0023] Perform one-dimensional Fourier transform on each channel of the common shot seismic record along the time direction and calculate the multi-channel average amplitude spectrum;

[0024] Determine the effective frequency band range of seismic records as [ω1,ω2], and design a filter with a passband of [ω1,ω2];

[0025] It is converted to the time domain through inverse Fourier transform, and the obtained time series is used as the offset wavelet w(t).

[0026] Furthermore, the forward wave field extension is performed using the offset wavelet w(t) as the source term of the constant density acoustic wave equation to obtain the source wave field p F The step of (x; t) further includes:

[0027] The forward extension of the source wavefield can be expressed as

[0028]

[0029] Among them, v(x) is the velocity model, x is the underground position coordinate, p F is the source wave field, w(t) is the source wavelet, and the integral of the source wavelet is used as the boundary condition.

[0030] Furthermore, the observed seismic data D(x r ,x s t) as the boundary condition of the constant density acoustic wave equation to perform reverse time wave field extension and obtain the received wave field p B The step of (x; t) further includes:

[0031] The reverse propagation of the received wavefield is

[0032]

[0033] Furthermore, the source wave field p F (x; t) and the received wave field p B (x; t) respectively use the Poynting vector method to estimate the direction vector P of the wave field propagation direction s (x; t) and P r (x; t), according to P s (x; t) and P r The angle between (x and t) is used to calculate the reflection angle θ, and the P s (x; t) and P r The steps of estimating the formation dip angle δ by using the normal vector corresponding to the vector P(x; t) and the vector P(x; t) include:

[0034] For the source wave field p F (x; t) and the received wave field p B (x; t) respectively use the Poynting vector method to estimate the direction vector P of the wave field propagation direction s (x; t) and P r (x; t);

[0035]

[0036]

[0037] According to P s (x; t) and P r The angle between (x and t) is used to calculate the reflection angle θ:

[0038]

[0039] The Poynting vector is spatially smoothed and normalized to obtain

[0040]

[0041] Where Ω is a region around the imaging point, and ε is a regularization parameter to avoid the zero division problem;

[0042] Using P s (x; t) and P r (x; t) and vector P(x; t) are used to calculate the formation dip δ. Let P(x, t) = (P1, P2), and use the vector perpendicular to P(x; t) to represent the formation dip, that is,

[0043] Furthermore, the source wave field p F (x; t) and the received wave field p B The steps of applying the cross-correlation imaging condition to (x; t) and performing angle binning to distribute the imaging results to the corresponding reflection angles θ to obtain the angle domain common imaging point gathers R(x; θ) include:

[0044] The image of the underground structure is obtained by applying the cross-correlation imaging condition based on wavefield separation and the angle binning operation:

[0045]

[0046] Among them, m(x,x s ) is the stacked imaging section, T max is the total receiving time of earthquake records;

[0047]

[0048]

[0049] is the Fourier transform of the wave field in the depth direction.

[0050] Furthermore, the forward operator based on Kirchhoff approximation is used to As the boundary condition, the forward wave field is extended and the seismic wave field is recorded at the detection point to obtain the simulated seismic data d(x r ,x s ; t) step comprises:

[0051]

[0052] This includes the background wave field p F The forward propagation of The scattered wave field p is obtained as the boundary condition B , record the scattered wave field at the detection point position, and obtain the simulated seismic data d(x r ;t;x s )=p B (x r ;t;x s ).

[0053] The present invention also provides a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and is characterized in that when the processor executes the computer program, it implements the steps of the above-mentioned method for generating angle-preserving gathers based on least squares reverse time migration.

[0054] The present invention also provides a computer-readable storage medium storing a computer program. When the computer program is executed by a processor, the computer program implements the steps of the above-mentioned method for generating angle-preserving gathers based on least squares reverse time migration.

[0055] This paper provides a method for generating amplitude-preserving angle gathers based on least-squares reverse-time migration. This method extracts angle gathers from reverse-time migration as initial imaging results. This allows for quantitative interpretation of complex structures using LSRTM, where the image amplitude is presented as an angle-dependent reflection coefficient. This imaging gather can be directly used in seismic interpretation of amplitude-vary-angle (AVA). LSRTM, an inversion imaging method, allows for better amplitude preservation in the extracted angle gathers, improving the accuracy and reliability of AVA analysis and contributing to enhanced accuracy and precision in oil and gas detection, reservoir characterization, and fluid identification. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the structures shown in these drawings without paying any creative work.

[0057] Figure 1 Flowchart of a method according to an embodiment of the present invention.

[0058] Figure 2 Schematic diagram of wave propagation direction according to an embodiment of the present invention.

[0059] Figure 31 and 2 are the true velocity model (a) and the migration velocity model (b) of the Marmousi model according to an embodiment of the present invention.

[0060] Figure 4 The shot point of the embodiment of the present invention is located at 1.0 km and receives the observed earthquake record (a) and the simulated earthquake record (b).

[0061] Figure 5 Comparison of the angle gathers extracted at CDP 600 according to an embodiment of the present invention. (a) is the angle gather extracted using reverse time migration, and (b) is the angle gather generated after 6 iterations of least squares reverse time migration.

[0062] Figure 6 Comparison of stacking imaging results of embodiments of the present invention: (a) stacking imaging result of angle gathers extracted using reverse time migration; (b) stacking imaging result of angle gathers generated based on least squares reverse time migration.

[0063] Figure 7 4 is a comparison of the AVA curves of the embodiments of the present invention.

[0064] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and in conjunction with the embodiments. DETAILED DESCRIPTION

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

[0066] It should be noted that all directional indications in the embodiments of the present invention (such as up, down, left, right, front, back, etc.) are only used to explain the relative position relationship, movement status, etc. between the various components under a certain specific posture (as shown in the accompanying drawings). If the specific posture changes, the directional indication will also change accordingly.

[0067] In addition, the descriptions of "first", "second", etc. in the present invention are only for descriptive purposes and cannot be understood as indicating or implying their relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include at least one of the features. In addition, "and / or" in the full text includes three solutions. Taking A and / or B as an example, it includes technical solution A, technical solution B, and technical solution that satisfies both A and B. In addition, the technical solutions between the various embodiments can be combined with each other, but they must be based on the ability of ordinary technicians in this field to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be deemed that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0068] In some embodiments, a method for generating amplitude-preserving angle gathers based on least squares reverse time migration includes: applying the Poynting vector method to extract angle gathers as initial imaging results during the reverse time migration process, using an angle-dependent Kirchhoff forward operator with the angle gathers as input to obtain simulated seismic data, substituting the residuals between the simulated seismic data and the observed data into the reverse time migration imaging operator to obtain a gradient image, and using the conjugate gradient method to iteratively solve the least squares inversion problem to optimize the angle gathers.

[0069] Specifically, the Poynting vector method, in which the Poynting vector refers to the energy flux density vector in the electromagnetic field. It is used here as a directional vector to indicate the propagation direction of the seismic wave field. The "Kirchhoff" forward operator refers to the wave equation forward operator derived via the Kirchhoff integral method. The forward operator is the process of obtaining simulated seismic data through finite-difference simulation of the wave equation based on model parameters (here, the subsurface structure reflectivity model, i.e., the seismic image).

[0070] The conjugate gradient method (CGM) lies between the steepest descent method and the Newton method. It utilizes only first-order derivative information, overcoming the slow convergence of the steepest descent method while avoiding the Newton method's drawbacks of requiring the storage, calculation, and inversion of the Hesse matrix. The CGM is not only one of the most useful methods for solving large linear systems of equations, but also one of the most efficient algorithms for solving large-scale nonlinear optimization problems. Among various optimization algorithms, the CGM is a very important one. Its advantages include low memory requirements, step-wise convergence, high stability, and the absence of any external parameters. The CGM was first proposed by Hestenes and Stiefle. Building on this, Fletcher and Reeves further developed the CGM for solving nonlinear optimization problems. Because the CGM does not require matrix storage, offers advantages such as rapid convergence and quadratic termination, it has now become widely used in practical applications. The conjugate gradient method is a typical conjugate direction method. Each of its search directions is conjugate to each other, and these search directions d are simply a combination of the negative gradient direction and the search direction of the previous iteration. Therefore, it requires less storage and is easy to calculate.

[0071] The method for generating amplitude-preserving angle gathers based on least squares reverse time migration provided by the present invention has better amplitude preservation of the extracted angle gathers, and the illumination is more balanced (more balanced illumination means Figure 5 The amplitude variation between small and medium angle gathers and large angle gathers is more continuous and smooth. Figure 7 It can also be observed from the curve comparison in Figure 6 The amplitude on the left and right sides of the middle image is improved and balanced with the illumination energy in the middle position. The improvement in resolution means Figure 5 and Figure 6 The sidelobes of the mid-waveform are suppressed, and blurring effects are partially eliminated. This improves resolution and reduces the amplitude loss problem associated with uneven illumination and imperfect migration imaging operators when seismic waves propagate through complex media, significantly enhancing the accuracy and reliability of AVA analysis. Therefore, ADCIGs provide the high-quality, high-fidelity prestack angular domain data necessary for AVO / AVA analysis and prestack elastic parameter inversion, helping to improve the accuracy and precision of oil and gas detection, reservoir characterization, and fluid identification.

[0072] like Figure 1 FIG. 4 shows the steps of a method for generating amplitude-preserving gathers based on least squares reverse time migration.

[0073] Specifically, S1, collect the original common shot point seismic records, and preprocess the collected original common shot points. The preprocessing includes removing the direct wave, the refracted wave, the ghost wave or the ground roll wave, etc. The common shot seismic records obtained after processing are called the observed seismic data, which is recorded as D(x r ,xs ; t), where x r Indicates the coordinates of the detection point, x s represents the coordinates of the earthquake source, and t represents the time.

[0074] S2. Construct the migration velocity v(x) through velocity analysis, where x = (x, z) is the spatial coordinate, ignoring the density, and assuming that the density ρ≡1.

[0075] S3, by observing the earthquake data D(x r ,x s ; t) spectrum analysis, constructing a wide-band offset wavelet function w(t).

[0076] The specific steps include: performing one-dimensional Fourier transform on the common shot seismic records along the time direction, calculating the multi-channel average amplitude spectrum, determining the effective frequency band range of the seismic records as [ω1, ω2], designing a filter with a passband of [ω1, ω2], converting it to the time domain through inverse Fourier transform, and using the obtained time series as the offset wavelet w(t).

[0077] S4. Use the offset wavelet w(t) as the source term of the constant density acoustic wave equation to perform forward wave field extension and obtain the source wave field p F (x; t).

[0078] The forward extension of the source wave field is described as:

[0079]

[0080] Among them, v(x) is the velocity model, x is the underground position coordinate, p F is the source wave field, w(t) is the source wavelet, and the integral of the source wavelet is used as the boundary condition.

[0081] S5, using the observed earthquake data D(x r ,x s t) as the boundary condition of the constant density acoustic wave equation to perform reverse time wave field extension and obtain the received wave field p B (x; t).

[0082] The reverse propagation of the received wavefield is

[0083]

[0084] S6, source wave field p F (x; t) and the received wave field p B (x; t) respectively use the Poynting vector method to estimate the direction vector P of the wave field propagation direction s (x; t) and P r (x; t), according to P s (x; t) and Pr The angle between (x and t) is used to calculate the reflection angle θ, and the P s (x; t) and P r The formation dip δ is estimated by the normal vector of the sum vector P(x; t).

[0085] The specific steps of S6 include:

[0086] For the source wave field p F (x; t) and the received wave field p B (x; t) respectively use the Poynting vector method to estimate the direction vector P of the wave field propagation direction s (x; t) and P r (x; t);

[0087]

[0088]

[0089] According to P s (x; t) and P r The angle between (x and t) is used to calculate the reflection angle θ:

[0090]

[0091] The Poynting vector is spatially smoothed and normalized to obtain

[0092]

[0093] Where Ω is a region around the imaging point, and ε is a regularization parameter to avoid the zero division problem;

[0094] Using P s (x; t) and P r (x; t) and vector P(x; t) are used to calculate the formation dip δ. Let P(x, t) = (P1, P2), and use the vector perpendicular to P(x; t) to represent the formation dip, that is, The wave field propagation direction is as follows Figure 2 shown.

[0095] S7, source wave field p F (x; t) and the received wave field p B (x; t) applies the cross-correlation imaging condition and performs angle binning to distribute the imaging results to the corresponding reflection angle θ to obtain the angle domain common imaging point gather R(x; θ), which is the initial imaging result.

[0096] The specific steps of S7 include:

[0097] The image of the underground structure is obtained by applying the cross-correlation imaging condition based on wavefield separation and the angle binning operation:

[0098]

[0099] Among them, m(x,x s ) is the stacked imaging section, T max is the total receiving time of earthquake records;

[0100]

[0101]

[0102] is the Fourier transform of the wave field in the depth direction.

[0103] S8, using the forward operator based on Kirchhoff approximation to As the boundary condition, the forward wave field is extended and the seismic wave field is recorded at the detection point to obtain the simulated seismic data d(x r ,x s ; t); Since the receiving wave field is unknown, the receiving wave field is obtained by reverse time propagation through the observed data or data residuals as the boundary conditions of the wave equation. In this scenario, the forward operator based on the Kirchhoff approximation has an input of the underground reflectivity model (i.e., seismic image) and an output of the scattered wave field (i.e., the receiving wave field obtained by simulation). The simulated seismic data is obtained by recording the receiving wave field at the surface receiving point position. In order to verify the reliability of the generated seismic image, the purpose is to obtain simulated seismic data through forward simulation based on the reflectivity model, and to verify its similarity with the observed data in order to evaluate the reflectivity model. This is a general paradigm for solving inverse problems. The forward process obtains the observed data according to the model parameters, where the observed data is unknown. The inversion process obtains the model parameters according to the observed data, where the model parameters are the content to be solved. Therefore, the method for calculating the reflection angle is changed to use the source wave field p F The direction vector P of (x; t) s It is expressed as the angle between (x; t) and the normal vector P(x; t) of the formation dip.

[0104] The specific steps of S8 include:

[0105]

[0106] This includes the background wave field p F Forward propagation of The scattered wave field p is obtained as the boundary condition B , record the scattered wave field at the detection point position, and obtain the simulated seismic data d(x r ;t;x s )=pB (x r ;t;x s ).

[0107] S9. According to the adjoint state method (the adjoint state method is a method for calculating the gradient of the objective function when solving the optimization problem under the constraints of partial differential equations), the residual d between the simulated seismic data and the observed seismic data is converted to r (x r ,x s ; t) = d(x r ,x s ; t)-D(x r ,x s t) is substituted into S4-S7 of the migration imaging generated angle gather as the boundary condition to obtain the least squares objective functional corresponding to The gradient image ΔR(x,θ) is the gradient image ΔR(x,θ)=L T (LR(x,θ)-D);

[0108] Where L represents the Kirchhoff forward operator, L T represents the migration imaging operator for generating angle gathers, and D represents the observation data.

[0109] Since the scattered wave field is unknown, the method for calculating the reflection angle in step S6 is changed to use the source wave field p F The direction vector P of (x; t) s The angle between (x; t) and the normal vector P(x; t) of the formation dip angle δ is expressed. P(x; t) can be expressed by δ as P(x, t) = (-sinδ, cosδ), then the reflection angle can be expressed as:

[0110]

[0111] S10, iteratively update the diagonal gather using the conjugate gradient method R n (x,θ)=R n-1 (x,θ)+αΔR g (x,θ), α is the update step size calculated by the conjugate gradient method, αR g (x,θ) is the conjugate gradient direction.

[0112] S11. Repeat steps S8-S10 until the data residual converges to an acceptable level or reaches a certain number of iterations.

[0113] In this embodiment, the Marmousi model is used to conduct numerical experiments, and RTM and LSRTM are used to extract angle gathers and compare them. Figure 3As shown, the true velocity and offset velocity models are given. The lateral width and depth of the model are 7.5 km and 3.75 km respectively, with a grid spacing of 10 m and a grid size of 750*375. The source starts at 50.0 m on the left side of the model, with a shot spacing of 75.0 m and a total of 100 shots. A bilateral receiving observation system is used, with 750 geophones arranged in each grid in the horizontal direction. The depth of both the source and the geophone is 10.0 m. The source uses a Ricker wavelet with a main frequency of 20 Hz. The angular range of ADCIGs is 0-60°, with an angular sampling interval of 2°.

[0114] Figure 4 (a) shows the common shot seismic record received at a shot point located 1.0 km away.

[0115] First, the Poynting vector method is applied in reverse time migration to extract the angle gathers as the initial imaging results, such as Figure 5 As shown in (a); then, unlike the conventional LSRTM which uses a forward operator based on the Born approximation to take the stacked image as input to obtain simulated seismic data, this embodiment uses a forward operator based on the Kirchhoff approximation to take the angle gather as input to obtain simulated seismic data.

[0116] like Figure 4 As shown in (b), the least squares inversion problem is solved by the conjugate gradient method to optimize the angle gathers, and the results are as follows Figure 5 (b) Compare Figure 5 From the results in (a) and 5(b), we can see that at the circled position, LSRTM improves the continuity of the angular gather amplitude and the balance of illumination. At the position indicated by the arrow, the weak signal in the original RTM angular gather is enhanced in the LSRTM angular gather, and the resolution is improved. The angular gathers are superimposed to obtain the stacked imaging result as shown in the figure. Figure 6 As shown in Figure 6 (a) and 6(b) are the superimposed images of RTM and LSRTM respectively. By comparison, we can see that, as shown by the arrows, LSRTM effectively improves the imaging resolution, partially alleviates the side lobe effect of the in-phase axis sub-wave, makes the illumination more balanced, suppresses the noise, and makes the imaging amplitude more fidelity. Figure 7 Comparing the curves of reflection coefficient changing with angle, we can see that the AVA curve obtained by LSRTM is in better agreement with the theoretical curve.

[0117] The above description is only a preferred embodiment of the present invention and does not limit the patent scope of the present invention. All equivalent structural transformations made by using the contents of the present invention description and drawings under the inventive concept of the present invention, or direct / indirect application in other related technical fields are included in the patent protection scope of the present invention.

Claims

1. A method for generating amplitude-preserving gathers based on least squares reverse time migration, characterized in that: The steps include: Extracting angle gathers as initial imaging results in reverse time migration; Taking the initial imaging results as input, the wave equation forward operator based on Kirchhoff approximation is used to obtain simulated seismic data; Substituting the residual between the simulated seismic data and the observed data into a reverse time migration imaging operator to obtain a gradient image; The conjugate gradient method is used to iteratively solve the least squares inversion problem to optimize the angle gathers; The steps include: Collect original common shot point seismic records; The collected original common shot point seismic records are preprocessed, and the common shot seismic records obtained after processing are called observed seismic data, which are recorded as D(x r ,x s ; t), where x r Indicates the coordinates of the detection point, x s represents the coordinates of the earthquake source, and t represents the time; The depth domain migration velocity v(x) is constructed through velocity analysis, where x = (x, z) is the spatial coordinate, and the density is not considered, assuming that the density ρ≡1; By observing the seismic data D(x r ,x s ; t) spectrum analysis, constructing a wide-band offset wavelet function w(t); The forward wave field extension is performed by using the migration wavelet w(t) as the source term of the constant density acoustic wave equation to obtain the source wave field p F (x; t); Taking the observed earthquake data D(x r ,x s t) as the boundary condition of the constant density acoustic wave equation to perform reverse time wave field extension and obtain the received wave field p B (x; t); For the source wave field p F (x; t) and the received wave field p B (x; t) respectively use the Poynting vector method to estimate the direction vector P of the wave field propagation direction s (x; t) and P r (x; t), according to P s (x; t) and P r The angle between (x and t) is used to calculate the reflection angle θ, and the P s (x; t) and P r The normal vector of the sum vector P(x; t) of (x; t) estimates the formation dip δ; For the source wave field p F (x; t) and the received wave field p B (x; t) Apply the cross-correlation imaging condition and perform angle binning to distribute the imaging results to the corresponding reflection angle θ to obtain the angle domain common imaging point gather R(x; θ), and use the angle domain common imaging point gather R(x; θ) as the initial imaging result; The wave equation forward operator based on Kirchhoff approximation is used to As the boundary condition, the forward wave field is extended and the seismic wave field is recorded at the detection point to obtain the simulated seismic data d(x r ,x s ;t); According to the adjoint state method, the residual d between the simulated seismic data and the observed seismic data is r (x r ,x s ; t) = d(x r ,x s ; t)-D(x r ,x s t) is substituted into the step of generating angle gathers in migration imaging as the boundary condition to obtain the least squares objective functional corresponding to The gradient image ΔR(x,θ)=L T (LR(x,θ)-D); Where L represents the Kirchhoff forward operator, L T represents the migration imaging operator for generating angle gathers, and D represents the observation data; Iteratively update the diagonal gather using the conjugate gradient method R n (x,θ)=R n-1 (x,θ)+αΔR g (x,θ), α is the update step size calculated by the conjugate gradient method, ΔR g (x,θ) is the conjugate gradient direction.

2. The method according to claim 1, characterized in that The observed seismic data D(x r ,x s ; t), and the steps of constructing a wide-band offset wavelet function w(t) include: Perform one-dimensional Fourier transform on each channel of the common shot seismic record along the time direction and calculate the multi-channel average amplitude spectrum; Determine the effective frequency band range of seismic records as [ω1,ω2], and design a filter with a passband of [ω1,ω2]; It is converted to the time domain through inverse Fourier transform, and the obtained time series is used as the offset wavelet w(t).

3. The method according to claim 2, characterized in that The forward wave field extension is performed using the offset wavelet w(t) as the source term of the constant density acoustic wave equation to obtain the source wave field p F (x; The step of t) also includes: the forward continuation of the source wave field is described as Among them, v(x) is the velocity model, x is the underground position coordinate, p F is the source wave field, w(t) is the source wavelet, and the integral of the source wavelet is used as the boundary condition.

4. The method according to claim 3, characterized in that The observed seismic data D(x r ,x s t) as the boundary condition of the constant density acoustic wave equation to perform reverse time wave field extension and obtain the received wave field p B (x; t) step also includes: The back propagation of the received wavefield is described as 5. The method according to claim 4, characterized in that The source wave field p F (x; t) and the received wave field p B (x; t) respectively use the Poynting vector method to estimate the direction vector P of the wave field propagation direction s (x; t) and P r (x; t), according to P s (x; t) and P r The angle between (x and t) is used to calculate the reflection angle θ, and the P s (x; t) and P r The steps of estimating the formation dip angle δ from the normal vector of the sum vector P(x; t) include: For the source wave field p F (x; t) and the received wave field p B (x; t) respectively use the Poynting vector method to estimate the direction vector P of the wave field propagation direction s (x; t) and P r (x; t); According to P s (x; t) and P r The angle between (x and t) is used to calculate the reflection angle θ: The Poynting vector is spatially smoothed and normalized to obtain Where Ω is a region around the imaging point, and ε is a regularization parameter to avoid the zero division problem; Using P s (x; t) and P r (x; t) and vector P(x; t) are used to calculate the formation dip δ. Let P(x, t) = (P1, P2), and use the vector perpendicular to P(x; t) to represent the formation dip, that is, 6. The method according to claim 5, characterized in that The source wave field p F (x; t) and the received wave field p B The steps of applying the cross-correlation imaging condition to (x; t) and performing angle binning to distribute the imaging results to the corresponding reflection angles θ to obtain the angle domain common imaging point gathers R(x; θ) include: The image of the underground structure is obtained by applying the cross-correlation imaging condition based on wavefield separation and the angle binning operation: Among them, m(x,x s ) is the stacked imaging section, T max is the total receiving time of earthquake records; is the Fourier transform of the wave field in the depth direction.

7. The method according to claim 6, characterized in that The forward operator based on Kirchhoff approximation is used to As the boundary condition, the forward wave field is extended and the seismic wave field is recorded at the detection point to obtain the simulated seismic data d(x r ,x s ; The steps of t) include: This includes the background wave field p F The forward propagation of The scattered wave field p is obtained as the boundary condition B , record the scattered wave field at the detection point position, and obtain the simulated seismic data d(x r ;t;x s )=p B (x r ;t;x s ).

8. A computer device comprising a memory and a processor, wherein the memory stores a computer program, When the processor executes the computer program, the steps of the method for generating angle-preserving gathers based on least squares reverse time migration according to any one of claims 1 to 7 are implemented.

9. A computer-readable storage medium storing a computer program, wherein when the computer program is executed by a processor, the computer program implements the steps of the method for generating angle-preserving gathers based on least squares reverse time migration according to any one of claims 1 to 7.

Images