A method for generating point spread function in reflection angle domain based on LSTM network

Abstract

The present invention discloses a method for generating a reflection angle domain point spread function based on an LSTM network, belonging to the field of geophysical exploration technology, comprising the following steps: S1, selecting reference points and constructing labels through ray tracing and wave equation forward simulation; S2, training an LSTM network as a ray tracing proxy model for the entire imaging domain; S3, predicting traveltime gradients and amplitude factors through the ray tracing proxy model; S4, pairing seismic sources with receivers and calculating reflection angles through traveltime gradients; S5, solving for the illumination wave number vector and solving for the wave number domain point spread function in combination with the illumination intensity; S6, performing an inverse Fourier transform on the wave number domain point spread function to obtain a spatial domain point spread function; and S7, solving a deconvolution problem and optimizing the processing of angle gathers. The present invention can improve the computational efficiency of solving Green's functions 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 illumination intensity.

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 existing technologies, migration imaging is a key step in reflection seismic exploration, and seismic imaging results are typically expanded into the angular domain for further analysis and interpretation. However, limited observation apertures and complex overlying strata often result in suboptimal imaging results. Uneven seismic illumination leads to unreliable amplitude variations in angular gathers. Band-limited seismic data and wavelet stretching in large-angle gathers also lead to low resolution.

[0003] A feasible solution is to implement imaging-domain least-squares migration using high-dimensional deconvolution of the point spread function. Extending the concept of imaging-domain least-squares migration to the angular domain allows the construction of a point spread function that depends on the reflection angle and can be used to optimize angular gathers of common imaging points in the angular domain. The fundamental element in constructing the point spread function is the Green's function, or an approximation of the Green's function, which describes the wave propagation process. 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 the reference point {x ref , constructing 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 an inverse Fourier transform on the wavenumber domain point spread function to obtain a spatial domain point spread function;

[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 amplitude factor corresponding to the starting and ending positions of the seismic ray is constructed {(x ref ,x s ),(k x ,k y , k z,A)}, where x s is the earthquake source position, 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 out the wave field excitation amplitude from the wave equation forward simulation, where 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 propagation direction between it and any imaging point x is predicted by the ray tracing proxy model, that is, the travel time gradient and

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

[0018] Preferably, in S5, the illumination wave number vector That is, the sum of the travel time gradients of the source and the receiver, corresponding to the direction vector of the normal to the illumination dip angle.

[0019] Preferably, the effective frequency band is traversed to construct the scattering wave number vector k d =-ωn d, and the illumination intensity on the corresponding scattering 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] Here, U(x) represents the neighborhood centered on x.

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

[0024] (1) The present invention can obtain the point spread function with a relatively low computational cost, and after the LSTM network is trained, it can be used to construct the reflection angle domain point spread function 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 the network prediction is about 96%.

[0025] (2) The present invention can correct the amplitude of the reflector in the imaging results, 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 improving the horizontal and vertical 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 This is a step diagram of a method for generating a reflection angle domain point spread function based on an LSTM network according to 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 This is a schematic diagram of predicting Green's function information of seismic waves through neural networks;

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

[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] Figure 9 is the local reflection coefficient model;

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

[0038] Figure 11 is the initial imaging angle gather;

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

[0040] Figure 13 is the original stacking imaging result;

[0041] Figure 14 is the superimposed imaging result after optimization processing;

[0042] Figure 15 is the comparison of single-channel reflection coefficients;

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

[0044] Figure 17 is the comparison of the transverse average wavenumber spectra;

[0045] Figure 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 with reference to the accompanying drawings and embodiments.

[0047] Example 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 wavenumber 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 the 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 amplitude factor corresponding to the starting and ending positions of the seismic ray is constructed {(x ref , x s ), (k x , k y , k z , A)}, where x s is the earthquake source position, 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 the wave field excitation amplitude from 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 Λ Serves 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 expressed as 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 as follows 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. At these reference points, the travel time gradient is solved by ray tracing 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. These factors are combined to construct a label sample set for training the neural network.

[0053] S3. Predict traveltime gradient and amplitude factors via ray tracing proxy model.

[0054] The source position x of the traversal observation system s and the detection point position x r , the ray propagation direction between it and any imaging point x is predicted by the ray tracing proxy model, that is, the travel time gradient and 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 by the travel time gradient. therefore,

[0056] like Figure 6 As shown, the angle definition of the wave propagation direction, θ s represents the incident wave angle, θ g Represents the reflected wave angle, θ d represents the inclination angle 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 That 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 Vector k in d .

[0059] Traverse the effective frequency band and construct the scattering wave number vector k d =-ωn d , and the illumination intensity on the corresponding scattering 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 in Figure 2, the constructed point spread function in the wavenumber domain is obtained.

[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, which is expressed as:

[0061]

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

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

[0064] Figure 9 is a local reflection coefficient model, Figure 10 is a local seismic image generated by two-dimensional convolution of the point spread function and the localized reflection coefficient. Figure 11 In the middle is the initial imaging angle gather, Figure 12 The angular gathers after optimization using the point spread function are shown in Figure 2. 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 the θ is in the best agreement with the initial imaging result. It can be seen that the quality of the angle gathers has been improved after the optimization process, and the phase axis is more clearly focused.

[0065] Figure 13 is the stacking imaging result of the original angle gather, Figure 14 This is the result of the optimized gather stacking imaging. It can be seen that the resolution of the optimized stacking imaging result has been improved, 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 Figure 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] Figure 16 and 17 The middle curves 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. Figure 18A set of Amplitude versus angle variation (AVA) curves of the gathers are compared, 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 characteristics in the gathers are better consistent with the theoretical values ​​after the point spread function optimization processing, maintaining the characteristics of the reflection wave amplitude varying with the angle, which is beneficial to the subsequent reservoir identification, parameter inversion and other work. Table 1 Prediction results of different neural networks

[0067] 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％

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

[0069] Table 2 Comparison of calculation time

[0070]

[0071] Table 2 shows a comparison of computational time. Using ray tracing for all calculations would require 3.54 hours, while using a neural network approach to establish the label set, train, and predict takes a total of 0.92 hours, improving computational efficiency by approximately 4 times.

[0072] Therefore, the present invention adopts the above-mentioned method of generating a reflection angle domain point spread function based on an LSTM network, utilizes a bidirectional long short-term memory (LSTM) network to reconstruct the Green's function, 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 a label consisting of a travel time gradient and an amplitude, thereby replacing the solution of the Green's function, constructing an angle-dependent PSF, and using the PSF for inversion to optimize the angle set. Numerical experiments conducted on a three-dimensional synthetic model show that this 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 are improved. The response of the amplitude to angle variation is corrected, and the amplitude characteristics are more reliable for further analysis and interpretation.

[0073] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit the same. 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 solutions of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for generating a point spread function in the reflection angle domain based on an LSTM network, characterized in that: The following steps are involved: S1, select the 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 an inverse Fourier transform on the wavenumber domain point spread function to obtain a spatial domain point spread function; S7, solve the deconvolution problem and optimize the processing of angle gathers; 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 amplitude factor corresponding to the starting and ending positions of the seismic ray is constructed {(x ref , x s ), (k x , k y , k z , A)}, where x s is the earthquake source position, x ref is the end point of the ray, k=(k x , k y , k z ) is the wave propagation direction vector, A is the Green function amplitude factor; The Green function amplitude factor A is obtained by picking out the wave field excitation amplitude from the wave equation forward simulation. The excitation amplitude is the maximum amplitude of the wave propagating through a certain position.

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 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).

3. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 2, characterized in that: In S3, the source position x of the traversal observation system s and the detection point position x r , the ray propagation direction between it and any imaging point x is predicted by the ray tracing proxy model, that is, the travel time gradient and 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 S4, due to therefore, Here, θ represents the reflection angle.

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 S5, the illumination wave number vector That is, the sum of the travel time gradients of the source and the receiver, corresponding to the direction vector normal to the illumination dip angle.

6. The method for generating a reflection angle domain point spread function based on an LSTM network according to claim 5, characterized in that: Traverse the effective frequency band and construct the scattering wave number vector k d =-ωn d , and the illumination intensity on the corresponding scattering 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.

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 S6, the wavenumber domain point spread function Perform inverse Fourier transform to obtain the spatial domain point spread function The expression is: Here, U(x) represents the neighborhood centered on x.

Images