Method for generating reflection angle domain point spread function based on LSTM network

Abstract

The invention discloses a method for generating a reflection angle domain point spread function based on an LSTM network, and belongs to the technical field of geophysical exploration, and the method comprises the following steps: S1, selecting a reference point, and constructing a label through ray tracing and wave equation forward modeling; s2, training an LSTM (Long Short Term Memory) network as a ray tracing agent model of the whole space of an imaging domain; s3, predicting a travel time gradient and an amplitude factor through a ray tracing agent model; s4, combining the seismic source and the geophone in pairs, and calculating a reflection angle through the travel time gradient; s5, solving an illumination wave number vector, and solving a wave number domain point spread function in combination with the illumination intensity; s6, performing inverse Fourier transform on the wavenumber domain point spread function to obtain a spatial domain point spread function; and S7, solving a deconvolution problem to optimize the angle gather. According to the method, the calculation efficiency of solving the Green function through seismic wave ray tracing can be improved, the imaging response described by the point spread function is eliminated from the imaging result, and the balance of the resolution of the imaging result and the illumination intensity can be improved.

Description

Technical Field

[0001] The present invention relates to the field of geophysical exploration technology, and in particular to a method for generating a reflection angle domain point spread function based on an LSTM network. Background Art

[0002] In the prior art, migration imaging is a key link in reflection seismic wave exploration, and the results of seismic wave imaging are usually extended to the angle domain for further analysis and interpretation. However, due to the limited observation aperture and complex overlying strata, the imaging results are often not ideal. The uneven illumination of seismic waves leads to unreliable amplitude variation patterns in the angular gathers. The band limitation of seismic data and the stretching of sub-waves in large-angle gathers lead to low resolution problems.

[0003] Using high-dimensional deconvolution of the point spread function to implement least squares migration in the imaging domain is a feasible solution. By extending the concept of least squares migration in the imaging domain to the angle domain, a point spread function that depends on the reflection angle can be constructed and used for the angle domain common imaging point gather, that is, the optimization of the angle gather. The basic element for constructing the point spread function is the Green's function that describes the wave propagation process, or the approximate Green's function. Therefore, further in-depth research is urgently needed. Summary of the invention

[0004] The purpose of the present invention is to provide a method for generating a point spread function in the reflection angle domain based on an LSTM network, which can improve the computational efficiency of solving the Green's function by seismic wave ray tracing, eliminate the imaging response described by the point spread function from the imaging result, and improve the resolution of the imaging result and the balance of the illumination intensity.

[0005] To achieve the above object, the present invention provides a method for generating a reflection angle domain point spread function based on an LSTM network, comprising the following steps:

[0006] S1. Select reference point {x ref}, construct labels through ray tracing and wave equation forward modeling;

[0007] S2. Train an LSTM network F Λ As a ray tracing proxy model for the entire space of the imaging domain;

[0008] S3, predicting travel time gradient and amplitude factors through ray tracing proxy model;

[0009] S4, combining the source and the receiver points in pairs, and calculating the reflection angle by the travel time gradient;

[0010] S5, solving the illumination wave number vector, and solving the point spread function in the wave number domain in combination with the illumination intensity;

[0011] S6, performing inverse Fourier transform on the point spread function in the wavenumber domain to obtain the point spread function in the space domain;

[0012] S7. Solve the deconvolution problem and optimize the processing of angle gathers.

[0013] Preferably, in S1, at the reference point {x ref}, the label of the direction vector on the reference point is obtained by ray tracing, the label of the amplitude factor is obtained by forward simulation of the wave equation, and the label of the direction vector and the amplitude factor corresponding to the starting and ending positions of the seismic ray are constructed {(x ref ,x s ),(k x ,k y ,k z ,A)}, where x s is the earthquake source location, x ref is the end point of the ray, k=(k x ,k y ,k z ) is the wave propagation direction vector, and A is the Green's function amplitude factor.

[0014] Preferably, the Green function amplitude factor A is obtained by picking up the wave field excitation amplitude in the wave equation forward simulation, and the excitation amplitude is the maximum amplitude of the wave propagating through a certain position.

[0015] Preferably, in S2, the forward simulation process of the ray tracing proxy model is expressed as F Λ {(x,x s )}=(k x ,k y ,k z ,A).

[0016] Preferably, in S3, the source position x of the traversal observation system is s and the detection point position x r , the ray tracing proxy model is used to predict the wave propagation direction of the ray between it and any imaging point x, that is, the travel time gradient ▽τ(x,x s ) and ▽τ(x,x r ).

[0017] Preferably, in S4, due to therefore, Here, θ represents the reflection angle.

[0018] Preferably, in S5, the illumination wave number vector n d =▽τ(x,x s )+▽τ(x,x r ), which is the sum of the travel time gradients of the source and the receiver, and corresponds to the direction vector normal to the illumination dip angle.

[0019] Preferably, the effective frequency band range is traversed to construct the scattering wave number vector k d =-ωn d , and the illumination intensity on the corresponding scattered wave number vector is obtained based on the Green function amplitude information Where s(ω) represents the frequency domain wavelet, A s and A r They represent the amplitude factor of the source wave field and the amplitude factor of the receiver wave field, respectively. They are obtained by forward modeling of the wave equation and are superimposed on the wave number vector corresponding to the wave number spectrum of the point spread function. in, is the point spread function in the wavenumber domain.

[0020] Preferably, in S6, the wavenumber domain point spread function Perform inverse Fourier transform to obtain the spatial domain point spread function The expression is:

[0021]

[0022] Among them, U(x) represents the neighborhood centered on x.

[0023] Therefore, the beneficial effects of the present invention using the above-mentioned method for generating a reflection angle domain point spread function based on an LSTM network are as follows:

[0024] (1) The present invention can obtain a point spread function with a relatively low computational cost, and after the LSTM network is trained, it can be used to construct a point spread function in the reflection angle domain at any position in the imaging space. Compared with the conventional ray tracing method, the computational efficiency of the present invention can be increased by about 4 times, and the accuracy of network prediction is about 96%.

[0025] (2) The present invention can correct the amplitude of the reflector in the imaging result, improve the problem of amplitude imbalance caused by uneven illumination at different positions when the wave propagates underground, and improve the fidelity of the amplitude response with angle changes.

[0026] (3) The present invention can improve the blurring effect in the imaging results due to the limited frequency band of the source wavelet, thereby achieving the effect of improving the lateral and longitudinal resolution of the imaging results.

[0027] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 It is a step diagram of a method for generating a reflection angle domain point spread function based on an LSTM network of the present invention;

[0029] Figure 2 is the real speed model used for testing in the present invention;

[0030] Figure 3 is the reflection coefficient model used for testing in the present invention;

[0031] Figure 4 It is a schematic diagram of predicting the Green’s function information of seismic waves through neural networks;

[0032] Figure 5 It is a schematic diagram of the seismic observation system and the direction of the rays;

[0033] Figure 6 It is the definition of the direction and angle of wave propagation;

[0034] Figure 7 is the point spread function in the wavenumber domain;

[0035] Figure 8 is the point spread function in the spatial domain;

[0036] Fig. 9 is the local reflection coefficient model;

[0037] Fig.10 It is a local seismic image synthesized by point spread function;

[0038] Fig.11 is the initial imaging angle gather;

[0039] Fig.12 is the angle gather after point spread function optimization;

[0040] Fig.13 is the original stacking imaging result;

[0041] Fig.14 It is the stacked imaging result after optimization processing;

[0042] Fig.15 It is the comparison of single-channel reflection coefficients;

[0043] Fig.16 is the comparison of the longitudinal average wavenumber spectra;

[0044] Fig.17 is the comparison of the lateral average wavenumber spectra;

[0045] Fig.18 It is a comparison of the response of amplitude to angle change. DETAILED DESCRIPTION

[0046] The technical solution of the present invention is further described below through the accompanying drawings and embodiments.

[0047] Embodiment 1

[0048] The present invention provides a method for generating a point spread function in the reflection angle domain based on an LSTM network. The point spread function is constructed in the wave number domain. First, the amplitude and propagation angle information of the Green's function are estimated at sparsely distributed reference points. Then, the LSTM network is trained to predict the Green's function in the entire space. Finally, the Green's function is used to construct a point spread function that depends on the reflection angle. Figure 1 As shown, the specific steps include:

[0049] S1. Select reference point {x ref}, construct labels through ray tracing and wave equation forward modeling.

[0050] At the reference point {x ref}, the label of the direction vector on the reference point is obtained by ray tracing, the label of the amplitude factor is obtained by forward simulation of the wave equation, and the label of the direction vector and the amplitude factor corresponding to the starting and ending positions of the seismic ray are constructed {(x ref ,x s ),(k x ,k y ,k z ,A)}, where x s is the earthquake source location, x ref is the end point of the ray, k=(k x ,k y ,k z ) is the wave propagation direction vector, and A is the Green's function amplitude factor. The Green's function amplitude factor A is obtained by picking up the wave field excitation amplitude in the wave equation forward simulation. The excitation amplitude is the maximum amplitude of the wave propagating through a certain position.

[0051] S2. Train an LSTM network F Λ As a ray tracing proxy model for the entire space of the imaging domain.

[0052] The forward simulation process of the ray tracing proxy model is represented by F Λ {(x,x s )}=(k x ,k y ,k z ,A). The velocity model used for numerical testing is as follows Figure 2 As shown, the reflection coefficient model is Figure 3 shown. Figure 4 The figure shows a schematic diagram of predicting the Green's function through a neural network, where the red ones are reference reflection points. Ray tracing is used to solve the travel time gradient at these reference points as the label of the wave propagation direction. The Green's function amplitude factor is obtained by finite difference forward simulation of the wave equation and picking up the excitation amplitude. The label sample set is constructed together for training the neural network.

[0053] S3. Predict travel time gradient and amplitude factor through ray tracing proxy model.

[0054] The source position x of the traversal observation system s and the detection point position x r , the ray tracing proxy model is used to predict the wave propagation direction of the ray between it and any imaging point x, that is, the travel time gradient ▽τ(x,x s ) and ▽τ(x,x r ).like Figure 5 As shown, k s represents the propagation direction of the incident wave from the earthquake source, k g Indicates the propagation direction of the reflected wave at the detection point.

[0055] S4, combine the source and the receiver points in pairs, and calculate the reflection angle through the travel time gradient, because therefore,

[0056] like Figure 6 The figure shows the angle definition of the wave propagation direction, θ s is the incident wave angle, θ g Represents the reflected wave angle, θ d represents the inclination of the geological body, θ represents the reflection angle,

[0057] S5. Solve the illumination wave number vector and solve the point spread function in the wave number domain in combination with the illumination intensity.

[0058] Illumination wave number vector n d =▽τ(x,x s )+▽τ(x,x r ), which is the sum of the travel time gradients of the source and the receiver, corresponds to the direction vector of the normal to the illumination dip angle, that is, Figure 5 and Figure 6 The vector k in d .

[0059] Traverse the effective frequency band range and construct the scattering wave number vector k d =-ωn d , and the illumination intensity on the corresponding scattered wave number vector is obtained based on the Green function amplitude information Where s(ω) represents the frequency domain wavelet, A s and A r They represent the amplitude factor of the source wave field and the amplitude factor of the receiver wave field, respectively. They are obtained by forward modeling of the wave equation and are superimposed on the wave number vector corresponding to the wave number spectrum of the point spread function. in, is the point spread function in the wave number domain. Figure 7 As shown, the constructed point spread function in the wavenumber domain is shown.

[0060] S6, wave number domain point spread function Perform inverse Fourier transform to obtain the spatial domain point spread function like Figure 8 As shown, this is the generated spatial domain point spread function, expressed as:

[0061]

[0062] Among them, U(x) represents the neighborhood centered on x.

[0063] S7. Solve the deconvolution problem and optimize the processing of angle gathers.

[0064] Fig. 9 is a local reflection coefficient model, Fig.10 It is a local seismic image generated by two-dimensional convolution of the point spread function and the localized reflection coefficient. Fig.11 In the middle is the initial imaging angle gather, Fig.12 The angular gathers after optimization using the point spread function are shown in Figure 1. The optimization problem can be expressed as solving the objective function That is, to solve a reflection angle domain imaging result m θ , so that it is consistent with the point spread function The image generated by the high-dimensional convolution of is in the best agreement with the initial imaging result. It can be seen that after the optimization process, the quality of the angle gather is improved and the event axis is more clearly focused.

[0065] Fig.13 is the stacking imaging result of the original angle gather, Fig.14 This is the stacked imaging result of the optimized gather. It can be seen that the resolution of the stacked imaging result has been improved after the optimization, the illumination of the deep layer has been compensated, and the illumination of the shallow and deep layers has become more balanced. Extract one of the data for comparison, such as Fig.15 As shown, the black curve is the theoretical value of the reflection coefficient, the blue curve is the initial imaging result, and the red curve is the result after point spread function deconvolution optimization processing. It can be seen that the amplitude fidelity of the imaging result is better after the optimization processing, which verifies the effectiveness of the illumination compensation.

[0066] Fig.16 and 17 The middle ones are comparisons of the longitudinal and transverse average wavenumber spectra of the stacked imaging results, where the blue curve is the wavenumber spectrum of the initial imaging result, and the red curve is the wavenumber spectrum of the imaging result after inversion optimization. It can be seen that both the transverse and longitudinal wavenumber spectra are broadened, confirming the effectiveness of the resolution improvement. Fig.18A set of Amplitude versus angle variation (AVA) curves of the gathers are compared in the figure, where the red curve is the theoretical value, the blue curve is the initial imaging result, and the red curve is the result after high-dimensional deconvolution optimization processing and correction. It can be seen that the AVA features in the gathers are better consistent with the theoretical values ​​after being optimized using the point spread function, and the characteristics of the reflection wave amplitude varying with the angle are maintained, which is beneficial to subsequent reservoir identification, parameter inversion and other work.

[0067] Table 1 Prediction results of different neural networks

[0068] RNN Bi-RNN LSTM Bi-LSTM Mean square error 0.00134 0.00102 0.00094 0.00085 Accuracy 85.74％ 89.53％ 92.15％ 95.96％

[0069] Table 1 shows the comparison of prediction results of different neural networks. It can be seen that the prediction accuracy of the bidirectional LSTM network is the highest.

[0070] Table 2 Comparison of calculation time

[0071]

[0072] Table 2 shows a comparison of computing time. If all calculations are performed using ray tracing, it takes 3.54 hours. However, using the neural network method, the time for comprehensive label set establishment, training, and prediction is only 0.92 hours, which improves the computing efficiency by about 4 times.

[0073] Therefore, the present invention adopts the above-mentioned method of generating a reflection angle domain point spread function based on an LSTM network, reconstructs the Green's function using a bidirectional long short-term memory (LSTM) network, uses a ray tracing method to obtain the propagation direction of the wave (travel time gradient), and accurately calculates the wavefront amplitude using a wave equation forward simulation. The LSTM network is trained by labels consisting of travel time gradients and amplitudes, thereby replacing the solution of the Green's function, constructing an angle-dependent PSF, using the PSF for inversion, and optimizing the angle set. Numerical experiments conducted on a three-dimensional synthetic model show that the method can improve the imaging quality of pre-stack angle gathers and post-stack seismic images, compensate for illumination, and improve the resolution of angle gathers. Both vertical and lateral resolutions have been improved. The response of amplitude changes with angle has been corrected, and the amplitude characteristics are more reliable for further analysis and interpretation.

[0074] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A method for generating a reflection angle domain point spread function based on an LSTM network, characterized in that: The following steps are involved: S1. Select reference point {x ref }, construct labels through ray tracing and wave equation forward modeling; S2. Train an LSTM network F Λ As a ray tracing proxy model for the entire space of the imaging domain; S3, predicting travel time gradient and amplitude factors through ray tracing proxy model; S4, combining the source and the receiver points in pairs, and calculating the reflection angle by the travel time gradient; S5, solving the illumination wave number vector, and solving the point spread function in the wave number domain in combination with the illumination intensity; S6, performing inverse Fourier transform on the point spread function in the wavenumber domain to obtain the point spread function in the space domain; S7. Solve the deconvolution problem and optimize the processing of angle gathers.

2. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 1, characterized in that: In S1, at the reference point {x ref }, the label of the direction vector on the reference point is obtained by ray tracing, the label of the amplitude factor is obtained by forward simulation of the wave equation, and the label of the direction vector and the amplitude factor corresponding to the starting and ending positions of the seismic ray are constructed {(x ref ,x s ),(k x ,k y ,k z ,A)}, where x s is the earthquake source location, x ref is the end point of the ray, k=(k x ,k y ,k z ) is the wave propagation direction vector, and A is the Green's function amplitude factor.

3. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 2, characterized in that: The Green function amplitude factor A is obtained by picking out the wave field excitation amplitude in the wave equation forward simulation. The excitation amplitude is the maximum amplitude of the wave propagating through a certain position.

4. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 3, characterized in that: In S2, the forward simulation process of the ray tracing proxy model is represented by F Λ {(x,x s )}=(k x ,k y ,k z ,A).

5. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 4, characterized in that: In S3, the source position x of the traversal observation system s and the detection point position x r , the ray tracing proxy model is used to predict the wave propagation direction of the ray between it and any imaging point x, that is, the travel time gradient ▽τ(x,x s ) and ▽τ(x,x r ).

6. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 5, characterized in that: In S4, due to therefore, Here, θ represents the reflection angle.

7. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 6, characterized in that: In S5, the illumination wave number vector n d =▽τ(x,x s )+▽τ(x,x r ), which is the sum of the travel time gradients of the source and the receiver, and corresponds to the direction vector normal to the illumination dip angle.

8. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 7, characterized in that: Traverse the effective frequency band range and construct the scattering wave number vector k d =-ωn d , and the illumination intensity on the corresponding scattered wave number vector is obtained based on the Green function amplitude information Where s(ω) represents the frequency domain wavelet, A s and A r They represent the amplitude factor of the source wave field and the amplitude factor of the receiver wave field, respectively. They are obtained by forward modeling of the wave equation and are superimposed on the wave number vector corresponding to the wave number spectrum of the point spread function. in, is the point spread function in the wavenumber domain.

9. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 8, characterized in that: In S6, the point spread function in the wavenumber domain Perform inverse Fourier transform to obtain the spatial domain point spread function The expression is: Among them, U(x) represents the neighborhood centered on x.

Images