True amplitude seismic wave imaging method based on inclination angle and reflection angle domain illumination compensation

Abstract

The invention relates to a true amplitude seismic wave imaging method based on inclination angle and reflection angle domain illumination compensation. The method comprises the following steps: acquiring the position of a seismic source point and the position of a detection point; according to the position of the seismic source point and the position of the detection point, acquiring a set of the positions of the seismic source and the detection point; performing simulation by taking each position in the set of the positions of the seismic source point and the detection point as the seismic source point, obtaining a wave field excitation moment and an excitation amplitude, and calculating a wave propagation angle corresponding to the excitation moment; according to the excitation amplitude and the wave propagation angle, obtaining an inclination angle and reflection angle domain illumination intensity; and performing normalization correction processing on the amplitudes of the inclination angle and the reflection angle gathers by using the illumination intensity of the inclination angle and the reflection angle domain, and superposing the corrected inclination angle and the reflection angle gathers to obtain a superposed imaging result of illumination equalization. Illumination in a seismic wave imaging result can be more balanced, and the amplitude fidelity is improved.

Description

Technical Field

[0001] The invention relates to the technical field of geophysical exploration, and in particular to a true amplitude seismic wave imaging method based on illumination compensation in the dip angle and reflection angle domain. 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 parameters 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 lateral inhomogeneity. 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 underground structures, thereby affecting the fidelity of imaging amplitude.

[0003] Least squares migration is used to correct the amplitude in seismic imaging results and improve the resolution. However, least squares migration is often solved iteratively through optimization algorithms such as gradient descent, which is 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 inclination angle and reflection angle domain, which can make the illumination in the seismic wave imaging result 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 angle and reflection angle domains, comprising:

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

[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 earthquake source point and the detection point position as an earthquake source point, obtaining the wave field excitation time and excitation amplitude, and calculating the wave propagation angle corresponding to the excitation time;

[0010] According to the excitation amplitude and the wave propagation angle, obtaining the illumination intensity in the inclination angle and reflection angle domain;

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

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

[0013] Each position in the set of the earthquake source point and the detection point position is regarded as the earthquake 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 earthquake source point;

[0014] The time when the maximum amplitude of the seismic wave field occurs when it passes through a certain position 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 an 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 illumination intensity in the inclination domain and the illumination 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 inclination domain is:

[0020]

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

[0022]

[0023] Among them, I d (x,θ d ) and I r (x,θ r ) represent the illumination 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 at 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 zero division, m d (x,θ d ) and m r (x,θ r ) are the dip angle and reflection angle gathers respectively.

[0027] Optionally, obtaining the illumination equalization stacking imaging result includes:

[0028] or

[0029] Among them, 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 gather and stack imaging results of complex structures in the middle and deep layers. (2) The present invention can correct the amplitude of the reflector in the angular gather and stack imaging results, and improve the problem of amplitude distortion caused by uneven illumination at different positions when the wave propagates 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 drawings required for use in the embodiments will be briefly introduced below. 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 creative work.

[0032] Figure 1 It is a flow chart of a true amplitude seismic wave imaging method based on illumination compensation in the dip angle and reflection angle domain 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 3is the 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 wave propagation direction and angle in an embodiment of the present invention;

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

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

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

[0040] Fig. 9 is the initial dip domain imaging gather of an embodiment of the present invention;

[0041] Fig.10 The dip gather after illumination compensation in an embodiment of the present invention;

[0042] Fig.11 is the initial reflection angle domain imaging gather of an embodiment of the present invention;

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

[0044] Fig.13 is the original superposition imaging result of the embodiment of the present invention;

[0045] Fig.14 This is the superimposed imaging result after illumination compensation according to the embodiment of the present invention. DETAILED DESCRIPTION

[0046] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. 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 creative work 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 only corrects the amplitude without considering the resolution problem, and is a robust and feasible alternative solution. The angle-dependent illumination intensity can be derived from the Hessian of the least squares migration, which is used for illumination compensation amplitude correction of angle gathers and stack imaging results. Angle domain illumination compensation can simultaneously correct the amplitude in gathers and stack imaging results for subsequent seismic interpretation and reservoir inversion work.

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

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

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

[0052] Each position in the set of source point and receiver point positions 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] According to the excitation amplitude and wave propagation angle, the illumination intensity in the inclination angle and reflection angle domain is obtained;

[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 the source point and the receiver point position is simulated as a source point to obtain the wave field excitation time and excitation amplitude, including:

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

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

[0058] Specifically, taking each position 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 epicenter 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 using the finite difference format. Do forward modeling, 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 to the source point x G as the boundary condition of the wave equation.

[0059] The excitation time is the time 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 according to the wave propagation angle, wherein 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 according to 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 , the 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 illumination intensity in the inclination domain is:

[0067]

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

[0069]

[0070] Among them, I d (x,θ d ) and I r (x,θ r ) represent the illumination 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 at the detection point.

[0071] Furthermore, the dip and reflection angle domain illumination intensity is used to perform normalization correction on the dip and reflection angle gather amplitudes, including:

[0072]

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

[0074] Furthermore, the stacked imaging results of obtaining illumination balance include:

[0075] or

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

[0077] The present invention will be further described below in conjunction with 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 function are estimated by the excitation amplitude and the wave field propagation direction at the excitation time. Then, the Green 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 a stacked imaging result. Figure 1 As shown, the specific steps include:

[0079] S1, input migration velocity model v(x), and source point location 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 control the kinematic information of seismic wave propagation accurately. The velocity model v(x) used for testing is shown in Figure 2 middle.

[0080] S2, scanning observation system to obtain a set S of source and receiver point locations G ={x G}={x s}∪{x g}.

[0081] S3, taking each position 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 epicenter G (x,x G ,t), where x is the underground location coordinate and t is the time.

[0082] Finite difference wave equation forward modeling refers to the use of finite difference numerical solution format to solve the acoustic wave equation Do forward modeling, 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 to the source point x G as the boundary condition 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 time is the time 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 at the time of excitation 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 moment of excitation 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] According to the geometric relationship between the incident and reflected rays, the inclination angle θ of the geological body is calculated d and the seismic wave reflection angle θ r . Figure 5 The definition of wave propagation direction and angle is shown in s represents the propagation direction of the incident wave from the earthquake source, i.e., the direction vector of the earthquake source point, k g represents 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, i.e., 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 center 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. Fig. 9 The middle is the initial dip gather, Fig.10 The middle one is the dip gather after illumination compensation. Fig.11 The initial reflection angle gathers are shown in Fig.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 Fig.13 The initial stacking imaging results are shown in Fig.14 The stacked imaging results after illumination compensation are shown in .

[0094] observe Figure 6 , Figure 7 and Figure 8 By comparing the illumination intensity in the middle and deep layers, it is found that the illumination intensity in the middle and deep layers is relatively weak, and the illumination intensity under complex structures and strong reflection layers is generally relatively weak. 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 middle and deep positions indicated by the arrows are enhanced. The present invention can enhance the angle gathers and stacked imaging results of complex structures in the middle and deep layers, can 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 only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A true amplitude seismic wave imaging method based on illumination compensation in the dip angle and reflection angle domain, characterized in that: include: Get the location of the earthquake source point and the location of the receiver point; 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 earthquake source point and the detection point position as an earthquake source point, obtaining the wave field excitation time and excitation amplitude, and calculating the wave propagation angle corresponding to the excitation time; According to the excitation amplitude and the wave propagation angle, obtaining the illumination intensity in the inclination angle and reflection angle domain; The amplitude of the dip angle and reflection angle gathers is normalized and corrected by using the illumination intensity in the dip angle and reflection angle domains, and the corrected dip angle and reflection angle gathers are superimposed to obtain a superimposed imaging result with balanced illumination.

2. The true amplitude seismic wave imaging method based on illumination compensation in the dip angle and reflection angle domain according to claim 1, characterized in that: Simulating each position in the set of the source point and the receiver point position as a source point to obtain the wave field excitation time and excitation amplitude includes: Each position in the set of the earthquake source point and the detection point position is regarded as the earthquake 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 earthquake source point; The time when the maximum amplitude of the seismic wave field occurs when it passes through a certain position 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 angle and reflection angle domain 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 angle and reflection angle domain according to claim 3, characterized in that: According to the excitation amplitude and the wave propagation angle, obtaining the illumination intensity in the inclination angle and reflection angle domain comprises: Calculating the geological body inclination and seismic wave reflection angle according to the wave propagation angle; The illumination intensity in the inclination domain and the illumination 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.

5. The true amplitude seismic wave imaging method based on illumination compensation in the dip angle and reflection angle domain according to claim 4, characterized in that: The illumination intensity in the inclination domain is: The illumination intensity in the reflection angle domain is: Among them, I d (x,θ d ) and I r (x,θ r ) represent the illumination 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 at the detection point.

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

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

Images