A true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains

Abstract

The present invention relates to a true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains, comprising: obtaining a source point position and a receiver point position; obtaining a set of source and receiver point positions based on the source point position and the receiver point position; simulating each position in the set of source and receiver point positions as a source point, obtaining a wave field excitation moment and an excitation amplitude, and calculating a wave propagation angle corresponding to the excitation moment; obtaining an illumination intensity in the dip and reflection angle domains based on the excitation amplitude and the wave propagation angle; performing normalization correction processing on the amplitude of the dip and reflection angle gathers using the dip and reflection angle domain illumination intensity, superimposing the corrected dip and reflection angle gathers, and obtaining a superimposed imaging result with equalized illumination. The present invention can make the illumination in the seismic wave imaging result more balanced and improve the amplitude fidelity.

Description

Technical Field

[0001] The present invention relates to the field of geophysical exploration technology, in particular to a true amplitude seismic wave imaging method based on inclination and reflection angle domain illumination compensation. Background Art

[0002] In the prior art, seismic data interpretation and reservoir identification require that seismic data have reliable amplitude information so that the physical properties of the underground medium can be inverted based on the amplitude information. Therefore, the fidelity of the amplitude in the seismic data migration imaging results is crucial. Reverse time migration uses the full-wave wave equation as the wave field extension operator, which can simulate various complex wave phenomena such as rotation waves and prismatic waves, and has no inclination restrictions, so it can effectively image complex media with strong laterally inhomogeneous properties. However, even if reverse time migration uses a true amplitude wave field extension operator, it is still not enough to achieve true amplitude seismic wave imaging. This is because the migration is only an adjoint operator of the forward simulation process, not a true inverse operator. Complex overlying strata and limited migration aperture will lead to uneven illumination of the underground structure, thereby affecting the fidelity of the imaging amplitude.

[0003] Least-squares migration is used to correct amplitudes in seismic imaging results and improve resolution. However, it is often solved iteratively using gradient descent-type optimization algorithms, which are computationally expensive and sometimes unstable. Summary of the Invention

[0004] The purpose of the present invention is to provide a true amplitude seismic wave imaging method based on illumination compensation in the dip angle and reflection angle domains, which can make the illumination in the seismic wave imaging results more balanced and improve the amplitude fidelity.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains, comprising:

[0007] Get the location of the earthquake source and the location of the receiver;

[0008] According to the earthquake source point position and the detection point position, obtaining a set of earthquake source point and detection point position;

[0009] Simulating each position in the set of the source point and the receiver point position as a source point, obtaining the wave field excitation time and excitation amplitude, and calculating the wave propagation angle corresponding to the excitation time;

[0010] Obtaining illumination intensity in the inclination angle and reflection angle domains according to the excitation amplitude and the wave propagation angle;

[0011] The dip angle and reflection angle domain illumination intensity is used to perform normalization correction processing on the amplitude of the dip angle and reflection angle gathers, and the corrected dip angle and reflection angle gathers are superimposed to obtain a superimposed imaging result with equalized illumination.

[0012] Optionally, simulating each position in the set of the source point and the receiver point positions as a source point to obtain the wavefield excitation time and excitation amplitude includes:

[0013] Each position in the set of the source point and the receiver point position is regarded as the source point, and the acoustic wave equation is forward simulated by a finite difference numerical solution format to obtain the seismic wave field of the source point;

[0014] The time when the seismic wave field passes through a certain position and the maximum amplitude occurs is obtained as the excitation time, and the maximum amplitude is obtained as the excitation amplitude.

[0015] Optionally, calculating the wave propagation angle corresponding to the excitation moment includes: using a Poynting vector method to calculate the excitation amplitude direction vector corresponding to the excitation moment, and calculating the wave propagation angle according to the direction vector.

[0016] Optionally, obtaining the illumination intensity in the inclination angle and reflection angle domain according to the excitation amplitude and the wave propagation angle includes:

[0017] Calculating the geological body inclination and seismic wave reflection angle according to the wave propagation angle;

[0018] The lighting intensity in the inclination domain and the lighting intensity in the reflection angle domain are calculated according to the inclination angle of the geological body, the reflection angle of the seismic wave and the excitation amplitude.

[0019] Optionally, the lighting intensity in the tilt angle domain is:

[0020]

[0021] The illumination intensity in the reflection angle region is:

[0022]

[0023] Among them, I d (x,θ d ) and I r (x,θ r ) represent the lighting intensity in the inclination angle and reflection angle domain, x is the underground position coordinate, θ d is the inclination of the geological body, x s is the location of the earthquake source, x g is the detection point position, A is the amplitude factor, δ is the Dirac function, θ r is the seismic wave reflection angle, θ s and θ gis the abbreviated information of the excitation angles of the source and receiver points θ(x,x s ),θ(x,x g ), A(x,x s ) is the excitation amplitude at the earthquake source, A(x,x g ) is the excitation amplitude of the detection point.

[0024] Optionally, performing normalization correction processing on the amplitude of the dip angle and reflection angle gather using the illumination intensity in the dip angle and reflection angle domain includes:

[0025]

[0026] in, and The distribution is the amplitude-corrected dip and reflection angle gathers, ε is a regularization parameter to avoid division by zero, m d (x,θ d ) and m r (x,θ r ) are the dip angle and reflection angle gathers respectively.

[0027] Optionally, obtaining the illumination-equalized stacked imaging result includes:

[0028] or

[0029] Where m(x) is the superposition imaging result of illumination equalization.

[0030] The beneficial effects of the present invention are as follows: (1) The present invention can enhance the angular gathers and stacked imaging results of complex structures at medium and deep depths. (2) The present invention can correct the amplitude of reflectors in the angular gathers and stacked imaging results, thereby improving the amplitude distortion problem caused by uneven illumination at different locations when waves propagate underground. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] 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. 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 these drawings without paying any creative work.

[0032] Figure 1 This is a flow chart of a true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains according to an embodiment of the present invention;

[0033] Figure 2 is a migration velocity model used for testing according to an embodiment of the present invention;

[0034] Figure 3The excitation amplitude data of a shot in an embodiment of the present invention;

[0035] Figure 4 The wave field propagation angle data of a shot at the excitation moment of an embodiment of the present invention;

[0036] Figure 5 Definition of the wave propagation direction and angle according to the embodiment of the present invention;

[0037] Figure 6 is the illumination intensity in the tilt domain according to the embodiment of the present invention;

[0038] Figure 7 is the illumination intensity in the reflected angle domain according to the embodiment of the present invention;

[0039] Figure 8 is the superimposed illumination intensity obtained after angle superposition according to an embodiment of the present invention;

[0040] Figure 9 is the initial dip domain imaging gather according to an embodiment of the present invention;

[0041] Figure 10 This is the dip gather after illumination compensation according to an embodiment of the present invention;

[0042] Figure 11 is the initial reflection angle domain imaging gather according to an embodiment of the present invention;

[0043] Figure 12 This is the reflection angle gather after illumination compensation according to an embodiment of the present invention;

[0044] Figure 13 is the original superimposed imaging result of an embodiment of the present invention;

[0045] Figure 14 This is the superimposed imaging result after illumination compensation in an embodiment of the present invention. DETAILED DESCRIPTION

[0046] 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. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0047] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0048] Illumination compensation, which corrects only for amplitude without considering resolution, is a robust alternative solution. The angle-dependent illumination intensity can be derived from the Hessian of the least-squares migration and used for amplitude correction in angle gathers and stacked imaging. Angle-domain illumination compensation simultaneously corrects amplitude in both gathers and stacked imaging, enabling subsequent seismic interpretation and reservoir inversion.

[0049] Therefore, this embodiment provides a true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains, including:

[0050] Get the location of the earthquake source and the location of the receiver;

[0051] According to the location of the earthquake source point and the location of the receiver point, a set of the location of the earthquake source and the location of the receiver point is obtained;

[0052] Each location in the set of source and receiver locations is simulated as a source point, the wave field excitation time and excitation amplitude are obtained, and the wave propagation angle corresponding to the excitation time is calculated;

[0053] Obtain the illumination intensity in the inclination and reflection angle domains based on the excitation amplitude and wave propagation angle;

[0054] The amplitude of the dip and reflection angle gathers is normalized and corrected using the illumination intensity in the dip and reflection angle domains. The corrected dip and reflection angle gathers are superimposed to obtain a stacked imaging result with balanced illumination.

[0055] Furthermore, each position in the set of source and receiver positions is simulated as a source point to obtain the wave field excitation time and excitation amplitude, including:

[0056] Each position in the set of source and receiver points is regarded as a source point, and the acoustic wave equation is forward simulated by a finite difference numerical solution format to obtain the seismic wave field of the source point;

[0057] The moment when the seismic wave field passes through a certain position and the maximum amplitude occurs is obtained as the excitation moment, and the maximum amplitude is obtained as the excitation amplitude.

[0058] Specifically, taking each location in the set as the earthquake source point, the wave equation finite difference forward simulation is performed to obtain the x G is the seismic wave field P at the earthquake source G (x,x G ,t), where x is the underground position coordinate and t is the time; the finite difference forward modeling of the wave equation refers to the numerical solution of the acoustic wave equation by the finite difference format. Do forward simulation, where P represents the seismic wave field P(x,x G ,t), represents the spatial Laplace operator, which is f s (t) represents the source wavelet, applied at the source point x G as the boundary conditions of the wave equation.

[0059] The excitation moment is the moment when the wave field passes through a certain position and the maximum amplitude occurs. The excitation amplitude is the maximum amplitude A(x,x G )=P G (x,x G ,T e ).

[0060] Furthermore, calculating the wave propagation angle corresponding to the excitation moment includes: using a Poynting vector method to calculate the excitation amplitude direction vector corresponding to the excitation moment, and calculating the wave propagation angle according to the direction vector.

[0061] Specifically, the corresponding excitation amplitude direction vector is calculated by the Poynting vector method as The corresponding propagation angle is where n z is the unit vector in the z direction.

[0062] Furthermore, according to the excitation amplitude and the wave propagation angle, the illumination intensity in the inclination angle and reflection angle domain is obtained, including:

[0063] Calculate the geological body inclination and seismic wave reflection angle based on the wave propagation angle, where the wave propagation angle includes the source point excitation angle information and the receiver point excitation angle information;

[0064] The illumination intensity in the inclination domain and the illumination intensity in the reflection angle domain are calculated based on the inclination angle of the geological body, the reflection angle of the seismic wave and the excitation amplitude, wherein the excitation amplitude includes the excitation amplitude at the source point and the excitation amplitude at the receiver point.

[0065] Specifically, the inclination angle θ of the geological body is calculated based on the geometric relationship between the incident and reflected rays. d and the seismic wave reflection angle θ r , read the source point and the receiver point excitation angle information θ(x,x s ),θ(x,x g ) is abbreviated as θ s and θ g , then the inclination of the geological body seismic wave reflection angle

[0066] Furthermore, the lighting intensity in the tilt angle domain is:

[0067]

[0068] The illumination intensity in the reflected angle area is:

[0069]

[0070] Among them, I d (x,θ d ) and I r (x,θ r ) represent the lighting intensity in the inclination angle and reflection angle domain, x is the underground position coordinate, θ d is the inclination of the geological body, x s is the location of the earthquake source, x g is the detection point position, A is the amplitude factor, δ is the Dirac function, θ r is the seismic wave reflection angle, θ s and θ g is the abbreviated information of the excitation angles of the source and receiver points θ(x,x s ),θ(x,x g ), A(x,x s ) is the excitation amplitude at the earthquake source, A(x,x g ) is the excitation amplitude of the detection point.

[0071] Furthermore, the normalization correction processing of the dip and reflection angle gather amplitudes using the illumination intensity in the dip and reflection angle domains includes:

[0072]

[0073] in, and The distribution is the amplitude-corrected dip and reflection angle gathers, ε is a regularization parameter to avoid division by zero, m d (x,θ d ) and m r (x,θ r ) are the dip angle and reflection angle gathers respectively.

[0074] Furthermore, obtaining the illumination-balanced stacked imaging results includes:

[0075] or

[0076] Where m(x) is the superposition imaging result of illumination equalization.

[0077] The present invention will be further described below with reference to the accompanying drawings:

[0078] This embodiment provides a true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains. First, the amplitude and propagation angle information of the Green's function are estimated by the excitation amplitude and the wave field propagation direction at the excitation moment. Then, the Green's function information is used to construct the illumination intensity in the dip and reflection angle domains. Finally, the illumination intensity is used to correct the amplitude in the dip and reflection angle gathers to obtain the stacked imaging result. Figure 1 As shown, the specific steps include:

[0079] S1, input migration velocity model v(x), and source point position x s and the detection point position x g Among them, the migration velocity model is generally obtained through seismic wave tomography inversion. As the background velocity model of migration imaging, it can accurately control the kinematic information of seismic wave propagation. The velocity model v(x) used for testing is shown in Figure 2 middle.

[0080] S2, the scanning observation system obtains a set S of source and receiver position G ={x G}={x s}∪{x g}.

[0081] S3, taking each position in the set as the earthquake source point, perform finite difference forward simulation of the wave equation to obtain the x G is the seismic wave field P at the earthquake source G (x,x G ,t), where x is the underground location coordinate and t is the time.

[0082] Finite difference forward modeling of wave equations refers to the numerical solution of acoustic wave equations using finite difference format. Do forward simulation, where P represents the seismic wave field P(x,x G ,t), represents the spatial Laplace operator, which is f s (t) represents the source wavelet, applied at the source point x G as the boundary conditions of the wave equation.

[0083] S4. Pick up the wave field excitation time T during the wave field forward simulation process e and the excitation amplitude A(x,x G ), calculate the excitation time T e The corresponding direction vector k G (x,x G ) and the wave propagation angle θ(x,x G ).

[0084] The excitation moment is the moment when the wave field passes through a certain position and the maximum amplitude occurs. The excitation amplitude is the maximum amplitude A(x,x G )=P G (x,x G ,T e ), Figure 3The excitation amplitude data of a shot are shown in Figure 2. The corresponding excitation direction vector is calculated by the Poynting vector method and is The corresponding propagation angle is where n z is the unit vector in the z direction, Figure 4 The wave propagation angle data of a shot at the excitation moment is shown in FIG.

[0085] S5, with the source position x G As the index, the excitation amplitude A(x,x G ) and angle information θ(x,x G ) is stored on the hard disk.

[0086] S6. Rescan the observation system and calculate the angular domain illumination intensity.

[0087] When taking x G is the earthquake source point x s Or detection point x g When the source point excitation amplitude A(x,x s ), the detection point excitation amplitude A(x,x g ) and the focal point angle information θ(x,x s ), detection point angle information θ(x,x g ).

[0088] Calculate the inclination angle θ of the geological body based on the geometric relationship between the incident and reflected rays d and the seismic wave reflection angle θ r . Figure 5 The definition of wave propagation direction and angle is shown in s It represents the propagation direction of the incident wave from the earthquake source, that is, the direction vector of the earthquake source point, k g Indicates the propagation direction of the reflected wave at the detection point, that is, the detection point direction vector, k d Represents the normal direction vector of the reflector, θ s Indicates the incident wave propagation angle, that is, the source point excitation angle information, θ g Indicates the propagation angle of the reflected wave, that is, the excitation angle information of the detection point, the reflection angle: inclination: Solve the lighting intensity and do angle binning and lighting intensity accumulation to get the lighting intensity I d (x,θ d ), I r (x,θ r ), the process can be expressed as:

[0089]

[0090] Among them, I d (x,θ d ) and Ir (x,θ r ) represent the illumination intensity in the inclination angle and reflection angle domains, respectively. Figure 6 and Figure 7 The illumination intensity in the inclination angle and reflection angle domains are shown in Fig. Figure 8 The middle is the superimposed lighting intensity obtained after angle superposition.

[0091] S7. Correct the imaging amplitude using the illumination intensity.

[0092] Use lighting intensity I d (x,θ d ), I r (x,θ r ) for the dip gather m d (x,θ d ) and reflection angle gather m r (x,θ r ) is normalized and corrected, which can be expressed as and Here, ε is a regularization parameter to avoid division by zero. and That is, the dip and reflection angle gathers after amplitude correction. Figure 9 The middle is the initial dip gather, Figure 10 The middle one is the dip gather after illumination compensation. Figure 11 The initial reflection angle gathers are shown in Figure 12 The reflection angle gather after illumination compensation is shown in .

[0093] The amplitude-corrected gathers are stacked according to the angle to obtain the stacked imaging result with equalized illumination, which can be expressed as: or Figure 13 The initial stacking imaging results are shown in Figure 14 The stacked imaging results after illumination compensation are shown in .

[0094] observe Figure 6 、 Figure 7 and Figure 8 The illumination intensity in the deep layer is found to be weak, and the illumination intensity is generally weaker under complex structures and strong reflective layers. After illumination equalization processing, by comparing the angle gathers and stacked imaging results before and after illumination compensation, it can be found that the imaging results of the deep layer indicated by the arrow are enhanced. The present invention can enhance the angle gathers and stacked imaging results of complex structures in the deep layer, correct the amplitude of the reflectors in the angle gathers and stacked imaging results, and improve the problem of amplitude distortion caused by uneven illumination at different positions when the wave propagates underground.

[0095] The embodiments described above are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by persons skilled in the art should fall within the scope of protection defined by the claims of the present invention.

Claims

1. A true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains, characterized in that: include: Get the location of the earthquake source and the location of the receiver; According to the earthquake source point position and the detection point position, obtaining a set of earthquake source point and detection point position; Simulating each position in the set of the source point and the receiver point position as a source point, obtaining the wave field excitation time and excitation amplitude, and calculating the wave propagation angle corresponding to the excitation time; Obtaining the illumination intensity in the inclination angle and reflection angle domains according to the excitation amplitude and the wave propagation angle, including: Calculating the geological body inclination and seismic wave reflection angle according to the wave propagation angle; Calculating the illumination intensity in the inclination domain and the illumination intensity in the reflection angle domain according to the inclination angle of the geological body, the reflection angle of the seismic wave, and the excitation amplitude; The illumination intensity in the tilt angle domain is: The illumination intensity in the reflection angle region is: Among them, I d (x,θ d ) and I r (x,θ r ) represent the lighting intensity in the inclination angle and reflection angle domain, x is the underground position coordinate, θ d is the inclination of the geological body, x s is the location of the earthquake source, x g is the detection point position, A is the amplitude factor, δ is the Dirac function, θ r is the seismic wave reflection angle, θ s and θ g is the abbreviated information of the excitation angles of the source and receiver points θ(x,x s ),θ(x,x g ), A(x,x s ) is the excitation amplitude at the earthquake source, A(x,x g ) is the excitation amplitude of the detection point; The dip angle and reflection angle domain illumination intensity is used to perform normalization correction processing on the amplitude of the dip angle and reflection angle gathers, and the corrected dip angle and reflection angle gathers are superimposed to obtain a superimposed imaging result with equalized illumination.

2. The true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains according to claim 1, characterized in that: Simulating each position in the set of the source point and the receiver point positions as a source point to obtain the wave field excitation time and excitation amplitude includes: Each position in the set of the source point and the receiver point position is regarded as the source point, and the acoustic wave equation is forward simulated by a finite difference numerical solution format to obtain the seismic wave field of the source point; The time when the seismic wave field passes through a certain position and the maximum amplitude occurs is obtained as the excitation time, and the maximum amplitude is obtained as the excitation amplitude.

3. The true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains according to claim 2, characterized in that: Calculating the wave propagation angle corresponding to the excitation moment includes: using a Poynting vector method to calculate the excitation amplitude direction vector corresponding to the excitation moment, and calculating the wave propagation angle according to the direction vector.

4. The true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains according to claim 1, characterized in that: Performing normalization correction processing on the amplitude of the dip angle and reflection angle gather using the illumination intensity in the dip angle and reflection angle domain includes: in, and The distribution is the amplitude-corrected dip and reflection angle gathers, ε is a regularization parameter to avoid division by zero, m d (x,θ d ) and m r (x,θ r ) are the dip angle and reflection angle gathers, respectively.

5. The true amplitude seismic wave imaging method based on illumination compensation in the dip and reflection angle domains according to claim 4, characterized in that: Obtaining the stacked imaging results of illumination equalization includes: or Where m(x) is the superposition imaging result of illumination equalization.

Images