A method for generating azimuth reflection angle gathers, computer equipment and computer readable storage medium

Abstract

The present invention discloses a method for generating azimuth reflection angle gathers, a computer device and a computer-readable storage medium. Among them, a method for generating reflection angle gathers includes recording the excitation amplitude information of the forward wave field of the earthquake source in the reverse time migration process based on the excitation amplitude imaging condition, the excitation amplitude information includes the excitation time corresponding to the excitation amplitude and the direction vector corresponding to the excitation amplitude; using high-dimensional wavelet function combined with the excitation amplitude information to approximately reconstruct the forward wave field of the earthquake source, and applying the cross-correlation imaging condition to obtain the imaging result; calculating the direction vector of the reverse wave field of the detection point in the process of reverse time extension; calculating the azimuth and reflection angle of the wave field illumination of the underground imaging point by combining the direction vector corresponding to the recorded excitation amplitude with the direction vector of the reverse wave field of the detection point; assigning the imaging result to the corresponding azimuth and reflection angle to obtain the azimuth reflection angle gather. The azimuth reflection angle gather generated by this method has higher imaging quality.

Description

Technical Field

[0001] The present invention relates to the field of seismic exploration technology, and in particular to a method for generating azimuth reflection angle gathers, a computer device and a computer-readable storage medium. Background Art

[0002] The underground angle domain common imaging point gather is an important output in the seismic migration imaging process. Angle domain common imaging point gather (ADCIG) is a seismic gather corresponding to a common imaging point indexed by the azimuth and reflection angle corresponding to the illumination of the underground seismic wave field, referred to as angle gather. It builds a bridge between seismic migration imaging and seismic attribute inversion, such as migration velocity analysis (Migration velocity analysis, MVA), amplitude variation with angle analysis (Amplitude variation with angle analysis, AVA), and inversion of velocity and impedance disturbances.

[0003] The reverse time migration algorithm is the most advanced depth migration imaging algorithm currently, which can be applied to laterally strongly inhomogeneous media. Generating azimuth reflection angle gathers in the process of three-dimensional reverse time migration is a hot research direction at present, but its industrial application has been limited by computational efficiency. The low computational efficiency mainly comes from the access to the forward wave field of the source in reverse time migration, or the need to reconstruct the forward wave field of the source through the stored wave field boundary during the reverse time extension process. The excitation amplitude imaging condition can effectively solve the problem of low computational efficiency.

[0004] However, in the three-dimensional case, the quality of the azimuth reflection angle gathers generated using the excitation amplitude imaging condition is difficult to guarantee. In particular, high-density data acquisition in three dimensions is very expensive. Due to the loss of waveforms in the excitation amplitude, the azimuth reflection angle gathers present incomplete and discontinuous event axes. This is because the imaging results are expanded and binned into different azimuth and reflection angle gathers, and the effective signals in the azimuth reflection angle gathers are too sparse, while the angle calculation is not completely accurate, resulting in low quality of the azimuth reflection angle gathers generated using the excitation amplitude imaging condition. Summary of the invention

[0005] Based on this, it is necessary to provide a method for generating azimuth reflection angle gathers, a computer device and a computer-readable storage medium to solve the technical problem of low quality of azimuth reflection angle gathers generated using excitation amplitude imaging conditions in the prior art.

[0006] The present invention provides a method for generating an azimuth reflection angle gather, comprising:

[0007] In the reverse time migration process based on the excitation amplitude imaging condition, the excitation amplitude information of the forward wave field of the source is recorded, and the excitation amplitude information includes the excitation time corresponding to the excitation amplitude and the direction vector corresponding to the excitation amplitude;

[0008] Using a high-dimensional wavelet function in combination with the excitation amplitude information to approximately reconstruct the source forward wave field, and applying a cross-correlation imaging condition to obtain an imaging result;

[0009] In the process of reverse time extension, the direction vector of the reverse propagation field of the detection point is calculated;

[0010] The azimuth angle and reflection angle of the wave field illumination of the underground imaging point are calculated by combining the direction vector corresponding to the recorded excitation amplitude with the direction vector of the reverse wave field of the detection point;

[0011] The imaging results are assigned to the corresponding azimuth angles and reflection angles to obtain azimuth reflection angle gathers.

[0012] Furthermore, in the reverse time migration process based on the excitation amplitude imaging condition, the step of recording the excitation amplitude information of the forward wave field of the source further includes:

[0013] Collect the original common shot point seismic records, pre-process the collected original common shot point seismic records, and the pre-processed original common shot seismic records are offset shot gather records, denoted as D(x r ,x s ; t), where x r Indicates the coordinates of the detection point, x s represents the earthquake source coordinates, and t represents the time;

[0014] 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;

[0015] By recording the offset shot gather D(x r ,x s ; t) spectrum analysis, constructing a wide-band offset wavelet function w(t);

[0016] The offset wavelet function w(t) is used as the source term of the constant density acoustic wave equation to carry out forward wave field extension, and the forward wave field p of the source is obtained. F (x; t), record the excitation amplitude A(x|x s ) and its corresponding excitation time T(x|x s );

[0017] Calculate the direction vector k of the forward wave field of the earthquake source by using the Poynting vector method s (x, t), and record the direction vector k corresponding to the excitation moments (x,T(x|x s )).

[0018] Further, the offset shot gather record D(x r ,x s ; t) spectrum analysis, the steps of constructing a wide-band offset wavelet function w(t) include:

[0019] The offset shot gather record D(x r ,x s t) perform one-dimensional Fourier transform along the time direction for each channel, calculate the multi-channel average amplitude spectrum, and determine the offset shot gather record D(x r ,x s ; The effective frequency band range of t) is [ω1,ω2];

[0020] A filter with a passband of [ω1, ω2] is designed and converted to the time domain through inverse Fourier transform. The obtained time series is used as the offset wavelet function w(t).

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

[0022] The forward propagation of the source forward wave field is described as:

[0023]

[0024] Where v(x) represents the velocity model, x represents the underground position coordinate, and p F is the forward wave field of the earthquake source, w(t) is the earthquake source wavelet, and the integral of the earthquake source wavelet is used as the boundary condition.

[0025] Furthermore, the step of approximately reconstructing the source forward wave field by using a high-dimensional wavelet function in conjunction with the excitation amplitude information, and applying a cross-correlation imaging condition to obtain an imaging result includes:

[0026] In the process of reverse time extension of the back propagation wave field, the azimuth angle α of the propagation direction of the forward propagation wave field of the earthquake source is calculated according to the direction vector corresponding to the excitation time. s and the pitch angle β s ;

[0027] Using the excitation amplitude and the corresponding azimuth angle α s , pitch angle β s Construct a high-dimensional wavelet function as an approximate representation of the Green function to obtain the reconstructed source forward wave field

[0028] Furthermore, the step of calculating the direction vector of the reverse propagation field of the detection point during the reverse time extension process includes:

[0029] Using offset shot gather records D(x r ,x s t) as the boundary condition of the constant density acoustic wave equation to carry out reverse time wave field extension, and obtain the reverse propagation wave field p of the detection point B (x; t).

[0030] Furthermore, the step of calculating the azimuth and reflection angle of the wave field illumination of the underground imaging point by combining the direction vector corresponding to the recorded excitation amplitude with the direction vector of the back-transmitted wave field of the detection point comprises:

[0031] For the back propagation field p B (x; t) The direction vector k of the wave field propagation direction is estimated using the Poynting vector method r (x; t);

[0032] According to the recorded excitation time direction vector k s (x,T(x|x s )), when t=T(x|x s ), through k s (x; t) and k r The angle between (x; t) is used to calculate the reflection angle θ and azimuth

[0033] Furthermore, the step of assigning the imaging result to the corresponding azimuth angle and the reflection angle to obtain an azimuth reflection angle gather includes:

[0034] Reconstructing the forward wave field of the earthquake source and the back propagation field p of the detection point B (x; t) cross-correlation imaging conditions are applied;

[0035] Perform angle binning to assign the imaging results to the azimuth of the underground observation system At the reflection angle θ, we get the azimuth reflection angle gather

[0036] The present invention also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor implements the steps of the above-mentioned method for generating azimuth reflection angle gathers when executing the computer program.

[0037] The present invention also provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps of the above-mentioned method for generating an azimuth reflection angle gather are implemented.

[0038] The present invention provides a method for generating azimuth reflection angle gathers. Since there is no need to reconstruct the forward wave field of the earthquake source in the reverse time delay through an effective boundary storage strategy, the method greatly improves the computational efficiency of generating azimuth reflection angle gathers. Compared with the excitation amplitude imaging condition, since we use high-dimensional wavelet approximation to reconstruct the main component of the forward wave field of the earthquake source, the azimuth reflection angle gathers generated by this method have higher imaging quality and stronger algorithm robustness. The method proposed by the present invention can be used for the efficient generation of azimuth reflection angle gathers, and can provide the necessary high-quality and high-fidelity pre-stack angle domain data volume for AVO / AVA analysis and pre-stack elastic parameter inversion, which helps to improve the accuracy and precision of oil and gas detection, reservoir characterization and fluid identification. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] 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 or the description of the prior art 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 the structures shown in these drawings without paying creative work.

[0040] Figure 1 Schematic diagram of the azimuth and reflection angle of the underground wave field illumination direction in an embodiment of the present invention.

[0041] Figure 2 Schematic diagram of wave propagation direction and propagation angle in an embodiment of the present invention.

[0042] Figure 3 is a diagram showing the process of obtaining an approximate forward wave field of the earthquake source by reconstructing the forward wave field of the earthquake source using a high-dimensional wavelet function in an embodiment of the present invention, wherein Figure 3(a) is the forward wave field of the earthquake source, Figure 3(b) is the excitation amplitude of the forward wave field of the earthquake source at the moment 1s, Figure 3(c) is the direction vector corresponding to the excitation amplitude of the forward wave field of the earthquake source at the moment 1s, and Figure 3(d) is an approximate representation of the forward wave field of the earthquake source reconstructed by a high-dimensional wavelet.

[0043] FIG. 4 shows the true velocity and offset velocity model of the artificial synthesis model in an embodiment of the present invention, wherein FIG. 4( a ) is the true velocity model and FIG. 4( b ) is the offset velocity model.

[0044] FIG5 is a comparison diagram of azimuthal reflection angle traces on synthetic seismic data in an embodiment of the present invention, wherein FIG5(a) is an azimuthal reflection angle trace set generated by excitation amplitude imaging conditions using synthetic model data, and FIG5(b) is an azimuthal reflection angle trace set generated by the method of the present invention using synthetic model data.

[0045] FIG6 is a comparison diagram of azimuth reflection angle gather imaging results on actual data in an embodiment of the present invention, wherein FIG6(a) is an azimuth reflection angle gather generated by excitation amplitude imaging conditions based on actual data, and FIG6(b) is an azimuth reflection angle gather generated by the method of the present invention based on actual data.

[0046] Figure 7 This is a diagram showing the superimposed imaging result of actual data in an embodiment of the present invention.

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

[0048] 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.

[0049] 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 components under a certain specific posture (as shown in the accompanying drawings). If the specific posture changes, the directional indication will also change accordingly.

[0050] In addition, the descriptions of "first", "second", etc. in the present invention are only used 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 it must be based on the ability of ordinary technicians in this field to implement. When the combination of technical solutions is contradictory or cannot be implemented, it should be deemed that such combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0051] In some embodiments, a method for generating an azimuth reflection angle gather comprises the following steps:

[0052] S100. In a reverse time migration process based on an excitation amplitude imaging condition, excitation amplitude information of a forward wave field of a seismic source is recorded, where the excitation amplitude information includes an excitation time corresponding to the excitation amplitude and a direction vector corresponding to the excitation amplitude.

[0053] S200, using high-dimensional wavelet function combined with excitation amplitude information to approximately reconstruct the source forward wave field, and applying cross-correlation imaging conditions to obtain imaging results.

[0054] S300, calculating the direction vector of the reverse propagation wave field of the detection point during the reverse time extension process.

[0055] S400, calculating the azimuth and reflection angle of the wave field illumination of the underground imaging point by combining the direction vector corresponding to the recorded excitation amplitude with the direction vector of the back-transmitted wave field of the detection point.

[0056] S500, assigning imaging results to corresponding azimuth angles and reflection angles to obtain azimuth reflection angle gathers.

[0057] This method is an efficient method for generating azimuthal reflection angle gathers based on the approximate representation of the wave field in high-dimensional phase space. Since there is no need to reconstruct the forward wave field of the earthquake source in the reverse time delay through the effective boundary storage strategy, this method greatly improves the computational efficiency of generating azimuthal reflection angle gathers. Compared with the excitation amplitude imaging condition, since we use the high-dimensional wavelet approximation to reconstruct the main component of the forward wave field of the earthquake source, the gathers generated by this method have higher imaging quality and stronger algorithm robustness. The method proposed in the present invention can be used for the efficient generation of azimuthal reflection angle gathers, and can provide the necessary high-quality and high-fidelity prestack angle domain data body for AVO / AVA analysis and prestack elastic parameter inversion, which is helpful to improve the accuracy and precision of oil and gas detection, reservoir characterization and fluid identification, and promote the industrial application of angle gather imaging technology, which is of great significance for the integration of geology and seismic engineering, efficient and low-cost exploration.

[0058] Specifically, S110, collect original common shot point seismic records, pre-process the collected original common shot point seismic records, and the pre-processed original common shot seismic records are offset shot gather records, denoted as D(x r ,x s ; t), where x r Indicates the coordinates of the detection point, x s represents the source coordinates, and t represents time. Preprocessing includes cutting off direct waves, removing refracted waves, ghost waves or ground roll waves, etc.

[0059] S120. Construct the migration velocity v(x) based on the velocity analysis, where x=(x,z) is the spatial coordinate, and the density is not considered, assuming that the density ρ≡1.

[0060] S130, by offsetting the shot collection record D(x r ,x s ; t) spectrum analysis, construct the wide-band migration wavelet function w(t). r ,x st) Perform one-dimensional Fourier transform along the time direction for each channel, calculate the multi-channel average amplitude spectrum, and determine the offset shot gather record D(x r ,x s ; t) is in the effective frequency band range of [ω1, ω2], a filter with a passband of [ω1, ω2] is designed, and it is converted to the time domain by inverse Fourier transform, and the obtained time series is used as the offset wavelet function w(t).

[0061] S140, using the offset wavelet function w(t) as the source term of the constant density acoustic wave equation to perform forward wave field extension, and obtain the source forward wave field p F (x; t), record the excitation amplitude A(x|x s ) and its corresponding excitation time T(x|x s ). The forward propagation of the source forward wave field is described as:

[0062]

[0063] Where v(x) represents the velocity model, x represents the underground position coordinate, and p F is the forward wave field of the earthquake source, w(t) is the earthquake source wavelet, and the integral of the earthquake source wavelet serves as the boundary condition.

[0064] S150, calculate the direction vector k of the forward wave field of the earthquake source by the Poynting vector method s (x, t), and record the direction vector k corresponding to the excitation moment s (x,T(x|x s )).

[0065] Record the excitation amplitude A(x|x s ), the excitation time T(x|x s ), and the corresponding direction vector k s (x|x s ,t):

[0066]

[0067] A(x|x s )=p F (x|x s ,T(x|x s ))

[0068]

[0069] Specifically, the Poynting vector method, the Poynting vector refers to the energy flux density vector in the electromagnetic field. Different from the excitation amplitude imaging condition in the related art, this step records the direction vector corresponding to the excitation moment in addition to the excitation amplitude and the corresponding excitation moment.

[0070] S210, in the process of reverse time extension of the back propagation wave field, the azimuth angle α of the propagation direction of the forward propagation wave field of the earthquake source is calculated according to the direction vector corresponding to the excitation time. s and the pitch angle β s (A schematic diagram of the wave field propagation direction is shown in Figure 2 middle):

[0071]

[0072]

[0073] α(x|x s )=α s (x|x s ,T(x|x s ))

[0074] β(x|x s )=β s (x|x s ,T(x|x s )).

[0075] S220, using the excitation amplitude and the corresponding azimuth angle α s , pitch angle β s Construct a high-dimensional wavelet function as an approximate representation of the Green function, and use the excitation amplitude and its propagation angle to reconstruct the forward wave field of the earthquake source.

[0076]

[0077] Among them, the high-dimensional wavelet function is expressed as:

[0078]

[0079] k c is the modulation wave number. Generally, in order to satisfy the wavelet admissibility condition, ||k is selected c ||2>=5.6. The translation, scaling and rotation of the wavelet function are expressed as b is the translation factor, a is the scale factor, is the rotation operator. The process of approximate reconstruction of the forward wave field of the earthquake source is shown in Figure 3.

[0080] In the related art, the maximum amplitude (excitation amplitude) of the forward wave field of the earthquake source at each imaging point and its corresponding occurrence time (excitation time) are recorded. In the process of reverse time extension of the reverse wave field at the detection point, the excitation amplitude is used to represent the forward wave field of the earthquake source, and the deconvolution imaging condition is applied to generate the seismic image. Compared with the cross-correlation imaging condition in the related art, the excitation amplitude imaging condition can improve the storage efficiency by more than three orders of magnitude and the calculation efficiency by about one third.

[0081] S310, record D(x r ,x s t) as the boundary condition of the constant density acoustic wave equation to carry out reverse time wave field extension, and obtain the reverse propagation wave field p of the detection point B (x; t). The reverse propagation of the detector wave field can be described as:

[0082]

[0083] The Poynting vector method is used to estimate the direction vector k of the wave field propagation direction for the detection wave back propagation field. r (x; t):

[0084]

[0085] According to k s (x; t) and k r The angle between (x; t) is used to calculate the reflection angle θ and azimuth

[0086]

[0087]

[0088] Where k = k s +k r , n x =(1,0,0) is the unit vector in the x direction. The azimuth and reflection angle of the underground wave field illumination direction are shown in the figure below: Figure 1 Different from the reverse time migration algorithm in the related art, which uses the process of reverse time extension of the wave field to calculate the direction vector of the forward wave field of the earthquake source to solve the wave field illumination angle, the present invention uses the direction vector of the forward wave field of the earthquake source corresponding to the recorded excitation time and the direction vector of the reverse wave field to solve the wave field illumination angle.

[0089] S410, for reconstructing the forward wave field of the earthquake source and the back propagation field p of the detection point B (x; t) Cross-correlation imaging conditions are applied.

[0090] S510, perform angle binning operation to distribute the imaging results to the azimuth of the underground observation system At the reflection angle θ, we get the azimuth reflection angle gather

[0091] Specifically, it can be described as applying the cross-correlation imaging condition based on wave field separation and the angle binning operation to obtain the image of the underground structure:

[0092]

[0093] Among them, T max is the total reception time of the earthquake record.

[0094] and is the Fourier transform of the wave field in the depth direction. max / 2 is used as the total imaging time to improve the computational efficiency without affecting the imaging quality. The cross-correlation imaging condition based on wavefield separation is applied to remove the low-wavenumber noise generated during the imaging process. At the same time, the source energy is used for normalization to perform illumination compensation as a preconditioning operator.

[0095] Numerical simulation results

[0096] A numerical experiment was conducted using an artificial synthetic model and a set of actual offshore seismic data. The azimuth reflection angle gathers were generated using the excitation amplitude imaging condition and the method of the present invention, and the two were compared. Figure 4 shows the true velocity and the offset velocity model of the artificial synthetic model. Using the true velocity model ( Figure 4a ), the observed seismic data is obtained by forward modeling a constant density acoustic wave equation. The source uses a Ricker wavelet with a main frequency of 20Hz. The observed seismic data contains a total of 260 shots, the shot spacing in the vertical and horizontal survey lines is 100m, and the first shot position is 200m. The detectors are distributed on the model surface at a spacing of 10m. The total recording time is 3s. The reflection angle range of the azimuth reflection angle gather is 0~40°, the angle sampling interval is 2°, and it is divided into 6 azimuths with an interval of 30°.

[0097] By extracting the gathers at the position of 1.3 km horizontally and 0.5 km vertically for comparison, the results are shown in Figure 5. Figure 5a The middle one is the angular gather obtained by excitation amplitude imaging condition. Figure 5b The angle gathers are generated by the method proposed in this invention. Figure 5a Compared with the results in 5b, we can see that due to the loss of waveform information in the excitation amplitude, Figure 5aThe angular gathers in the image are adversely affected, and the imaging quality is low. At the same time, each imaging point corresponds to a different direction vector. However, due to the complex overlapping wave fields, the calculation of the direction vector may not be completely correct. Since different parts of a waveform may calculate different direction vectors, the waveform fragmentation is assigned to different angles in the gather. By convolving a high-dimensional wavelet function with the excitation amplitude, the forward wave field of the source can be approximately reconstructed. Therefore, Figure 5b The gathers in the 3D image show complete waveforms, continuous phase axes, high resolution and high signal-to-noise ratio, and the imaging quality is greatly improved. Figure 6 compares the imaging results of the angle gathers of actual offshore data. It can be seen that Figure 6a The information in the angle gathers generated by the excitation amplitude imaging condition is too sparse and is interfered by noise, making it difficult to pick up information such as amplitude changes and curvature in the event axis, which is not conducive to subsequent seismic prestack attribute inversion. Figure 6b The angular gathers obtained by the method of the present invention have complete and continuous event axes, high signal-to-noise ratio, and the stacking imaging results are as follows: Figure 7 As shown in . In the numerical experiment, 8 RTX3090 GPU graphics cards were used for seismic imaging, and a parallel calculation strategy in the gun domain was adopted. When the excitation amplitude imaging condition was used, the total calculation time was 46.6 hours. When the method proposed in the present invention was used, 58.3 hours of calculation time was required, and the runtime was only increased by 25% compared with the excitation amplitude imaging condition, which is far less than the calculation complexity of the traditional reverse time migration algorithm.

[0098] 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 changes made by using the contents of the present invention specification and drawings under the inventive concept of the present invention, or directly / indirectly applied in other related technical fields are included in the patent protection scope of the present invention.

Claims

1. A method for generating azimuth reflection angle gathers, characterized in that: The following steps are involved: In the reverse time migration process based on the excitation amplitude imaging condition, the excitation amplitude information of the forward wave field of the source is recorded, and the excitation amplitude information includes the excitation time corresponding to the excitation amplitude and the direction vector corresponding to the excitation amplitude; Using a high-dimensional wavelet function in combination with the excitation amplitude information to approximately reconstruct the source forward wave field, and applying a cross-correlation imaging condition to obtain an imaging result; In the process of reverse time extension, the direction vector of the reverse propagation field of the detection point is calculated; The azimuth angle and reflection angle of the wave field illumination of the underground imaging point are calculated by combining the direction vector corresponding to the recorded excitation amplitude with the direction vector of the reverse wave field of the detection point; Allocating the imaging results to the corresponding azimuth angles and reflection angles to obtain azimuth reflection angle gathers; In the reverse time migration process based on the excitation amplitude imaging condition, the step of recording the excitation amplitude information of the forward wave field of the source further comprises: The original common shot point seismic records are collected, and the collected original common shot point seismic records are preprocessed. The original common shot seismic records after preprocessing are offset shot gather records, which are recorded as ,in represents the coordinates of the detection point, represents the earthquake source coordinates, Indicates time; Constructing migration velocities based on velocity analysis ,in is the spatial coordinate, ignoring the density, assuming the density ; By recording the offset shot collection Spectral analysis to construct a wideband offset wavelet function ; With the offset wavelet function As the source term of the constant density acoustic wave equation, the forward wave field is extended to obtain the forward wave field of the source , record the excitation amplitude at each position in the calculation area and its corresponding triggering moment ; Calculate the direction vector of the forward wave field of the earthquake source by using the Poynting vector method , and record the direction vector corresponding to the excitation moment ; The offset shot gathering record Spectral analysis to construct a wideband offset wavelet function The steps include: The offset shot gather record Perform one-dimensional Fourier transform along the time direction for each channel, calculate the multi-channel average amplitude spectrum, and determine the offset shot gather record The effective frequency band range is ; Design passband for The filter is converted to the time domain by inverse Fourier transform, and the obtained time series is used as the offset wavelet function .

2. The method for generating azimuth reflection angle gathers according to claim 1, characterized in that: The offset wavelet As the source term of the constant density acoustic wave equation, the forward wave field is extended to obtain the forward wave field of the source The steps also include: The forward propagation of the source forward wave field is described as: ； in, represents the velocity model, Indicates the underground location coordinates, is the positive wave field of the earthquake source, is the source wavelet, and the integral of the source wavelet serves as the boundary condition.

3. The method for generating azimuth reflection angle gathers according to claim 2, characterized in that: The step of approximately reconstructing the source forward wave field by using a high-dimensional wavelet function in conjunction with the excitation amplitude information, and applying a cross-correlation imaging condition to obtain an imaging result comprises: In the process of reverse time extension of the back propagation wave field, the azimuth of the propagation direction of the forward propagation wave field of the earthquake source is calculated according to the direction vector corresponding to the excitation time. and pitch angle ; Use the excitation amplitude and the corresponding azimuth , Pitch angle Construct a high-dimensional wavelet function as an approximate representation of the Green function to obtain the reconstructed source forward wave field .

4. The method for generating azimuth reflection angle gathers according to claim 3, characterized in that: The step of calculating the direction vector of the reverse propagation field of the detection point during the reverse time extension process includes: Recording with offset shot gathers As the boundary condition of the constant density acoustic wave equation, the reverse time wave field is extended to obtain the reverse propagation wave field of the detection point. .

5. The method for generating azimuth reflection angle gathers according to claim 4, characterized in that: The step of calculating the azimuth angle and reflection angle of the wave field illumination of the underground imaging point by combining the direction vector corresponding to the recorded excitation amplitude with the direction vector of the back-transmitted wave field of the detection point comprises: Counter-propagation field Using Poynting vector method to estimate the direction vector of wave field propagation direction ; According to the recorded excitation time direction vector ,when When, through and Calculate the reflection angle and azimuth .

6. The method for generating azimuth reflection angle gathers according to claim 5, characterized in that: The step of assigning the imaging result to the corresponding azimuth angle and the reflection angle to obtain an azimuth reflection angle gather comprises: Reconstructing the forward wave field of the earthquake source and the back propagation field of the detection point Applying cross-correlation imaging conditions; Perform angle binning to assign the imaging results to the azimuth of the underground observation system , Reflection Angle On the azimuth reflection angle gather, we get .

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the steps of a method for generating azimuth reflection angle gathers as described in any one of claims 1 to 6 are implemented.

8. A computer-readable storage medium storing a computer program, wherein the computer program, when executed by a processor, implements the steps of a method for generating an azimuthal reflection angle gather as described in any one of claims 1 to 6.

Images