A reverse time migration imaging method based on solving adaptive eikonal equation

Abstract

The invention discloses a reverse time migration imaging method based on solving an adaptive eikonal equation, which relates to the field of exploration geophysical technology, including: S1, inputting a model and source parameters; S2, obtaining a first arrival traveltime reference solution; S3, inputting a forward extension of the source, and saving a forward wave field corresponding to the first arrival traveltime; S4, inputting a reverse extension of a seismic record, and saving a reverse wave field corresponding to the first arrival traveltime; S5, applying imaging conditions for imaging. The method alternately scans in each direction based on the Fermat principle, adaptively selects local plane waves, spherical waves, and refraction wave operators, accurately simulates the first arrival traveltime of near, medium, and far migration, and obtains a first arrival traveltime reference solution; then the traveltime reference solution is used as the excitation moment, the corresponding forward wavefield and reverse wavefield are saved, and then the imaging conditions are applied to obtain a high-fidelity migration image, while also maintaining the advantage of efficient calculation.

Description

Technical Field

[0001] The invention relates to the field of exploration geophysical technology, and in particular to a reverse time migration imaging method based on solving an adaptive eikonal equation. Background Art

[0002] Imaging condition methods are key components in reverse time migration, such as excitation amplitude imaging and cross-correlation imaging conditions. Different imaging conditions show different capabilities in acquiring image amplitude, physical validity, and resolution. The cross-correlation imaging condition utilizes the wavefield information at all time points, can handle multipath problems, and achieve high imaging accuracy for complex structures such as salt domes. However, this method requires the storage of all forward wavefields, requiring a large amount of storage space, thereby reducing the computational efficiency of reverse time migration. In contrast, the excitation amplitude imaging condition only relies on the wavefield information at the imaging time, eliminating the necessity for large-scale storage, thereby improving the computational efficiency of reverse time migration. However, the traditional method determines the imaging time by identifying the maximum amplitude of the forward wavefield at each grid point. In cases involving large offsets or multipaths, this method may lead to inaccuracies in the determination of the imaging time, resulting in problems such as low imaging resolution and distortion.

[0003] Therefore, it is urgent to study a method to improve the quality of reverse time migration imaging. Summary of the invention

[0004] To solve the above technical problems, the present invention discloses a reverse time migration imaging method based on solving an adaptive eikonal equation. The method is based on the reverse time migration technology of the excitation amplitude imaging condition of the precise first arrival traveltime. Specifically, the traveltime calculated by the first arrival traveltime calculation method by accurately solving the eikonal equation is used as the excitation moment to improve the traditional excitation amplitude imaging condition, thereby improving the imaging quality of reverse time migration.

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

[0006] A reverse time migration imaging method based on solving an adaptive eikonal equation comprises the following steps:

[0007] S1, input model and source parameters;

[0008] S2, obtaining a reference solution for the initial trip;

[0009] S3, input the forward extension of the earthquake source and save the forward wave field corresponding to the first arrival travel time;

[0010] S4, reverse continuation of input seismic records, and saving the reverse propagation wave field corresponding to the first arrival travel time;

[0011] S5. Apply imaging conditions to perform imaging.

[0012] Optionally, in step S2, in the step of obtaining the first arrival travel time reference solution, the first arrival travel time in each direction is calculated using a plane wave operator, a spherical wave operator and a refraction wave operator respectively and iterated until the result is stable and converted into the first arrival travel time reference solution;

[0013] The calculation formula of the plane wave operator is:

[0014]

[0015] Where, T i,j represents the first arrival time at (i, j), T h and T v Variables set for convenience of expression;

[0016]

[0017] Where, d x and d z are the spatial sampling intervals of the model grid in the horizontal and vertical directions, S i-1,j-1 represents the slowness at (i-1, j-1), a, b and c are variables set for the convenience of expression, the same below;

[0018] The calculation formula of the spherical wave operator is:

[0019]

[0020] In the formula, T0 and S0 are known solutions, T 0,x and T 0,z are the derivatives of T0 with respect to x and z respectively;

[0021] The calculation formula of the refraction wave operator is:

[0022]

[0023] In the formula, min() means finding the minimum value from different variables, T rh and T rv They represent the travel time from point (i-1, j) and (i, j-1) to point (i, j) respectively;

[0024] The travel time of each point is calculated using the plane wave operator, spherical wave operator and refraction wave operator, and the minimum value is selected for iterative update until the initial arrival travel time is stable, and T is obtained. AFDM ;

[0025] First arrival travel time T calculated based on adaptive finite difference method AFDM In order not to ignore the influence of the first arrival wavefront tail part and the observation system layout, the first arrival travel time reference solution T is obtained by the following formula refor:

[0026]

[0027] In the formula, f m is the center frequency of the source wavelet. The factor a is set to reduce the error caused by the arrangement of the observation system, etc. The value range of a is [1,1.25].

[0028] Optionally, in step S3, in the step of inputting the forward continuation of the earthquake source and saving the forward wave field corresponding to the first arrival travel time, the earthquake source is input into the second-order velocity-stress acoustic wave equation for forward simulation, and the formula is as follows:

[0029]

[0030] Where P s is the forward wave field data at any time, v is the velocity vector, x and z are the x and z directions in the 2D plane respectively, t is the time, f(t) is the earthquake source, x s and z s Respectively represent the x-direction and z-direction positions of the earthquake source, δ(xx s ,zz s )f(t) means setting the source at (x s ,z s )

[0031] In the process of forward wave field continuation, the first arrival travel time reference solution obtained in step S2 is used as the excitation time, and the forward wave field corresponding to the excitation time is saved as the excitation amplitude, and the formula is as follows:

[0032] S max (x,z;T)=max{S(x,z;T)},T∈[T re -dt,T re +dt];

[0033] Where S(x,z;T) is the corresponding forward wave field within the time range of T, max{} is the function of taking the maximum value, S max (x,z;T) is the excitation amplitude.

[0034] Optionally, in step S4, input the reverse continuation of the seismic record, and in the step of saving the reverse propagation wave field corresponding to the first arrival travel time, in the second-order velocity-stress acoustic wave equation, input the reverse propagation wave field continuation of the seismic record simulation, and the formula is as follows:

[0035]

[0036] Where P ris the back propagation wave field data at any time, u is the seismic record, which is used as the boundary condition, and the back propagation wave field corresponding to the excitation time in step S3 is saved as R(x, z; T).

[0037] Optionally, in step S5, in the step of applying imaging conditions to image, the forward propagation wave field S at the excitation time saved in step S3 and step S4 is max (x, z; T) and the reverse wave field R(x, z; T) apply the excitation amplitude imaging condition, the formula is as follows:

[0038] I(x,z;s j )=R(x,z;T;s j )·S max (x,z;T;s j );

[0039] In the formula, s j represents the jth earthquake source, I(x,z;s j ) is the offset image of each shot, S max (x,z;T;s j ) is the forward wave field corresponding to the excitation time in each shot, R(x,z; T; s j ) is the back propagation wave field corresponding to the excitation moment in each shot.

[0040] The beneficial effect of the present invention is that the present invention improves the traditional excitation amplitude imaging conditions, adopts an adaptive finite difference method to solve the eikonal equation to calculate the first arrival travel time, the method is based on the Fermat principle to alternately scan in each direction, adaptively selects local plane waves, spherical waves and refraction wave operators, accurately simulates the first arrival travel time of near, medium and far offsets, and obtains the first arrival travel time reference solution; then the travel time reference solution is used as the excitation moment, the corresponding forward wave field and reverse wave field are saved, and then the imaging conditions are applied to obtain a high-fidelity offset image; finally, the Marmousi model is used to perform simulation experiments to prove the accuracy of the method of the present invention. The present invention adopts the above technical solution, while saving a lot of hard disk space and read and write tasks, the saved excitation moment and excitation amplitude are also more continuous and stable, solving the problem of artifacts and distortion of traditional excitation amplitude imaging conditions in complex media. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 A flowchart of a reverse time migration imaging method based on solving an adaptive eikonal equation according to the present invention;

[0042] Figure 2 It is the Marmousi speed model;

[0043] Figure 3A is the reverse time migration result under the cross-correlation imaging condition;

[0044] Figure 3B This is the reverse time migration result of the traditional excitation amplitude imaging condition;

[0045] Figure 3C The reverse time migration result of the improved excitation amplitude imaging condition of the present invention;

[0046] Figure 4A for Figure 3A The enlarged part of the middle area I;

[0047] Figure 4B for Figure 3B The enlarged part of the middle area I;

[0048] Figure 4C for Figure 3C The enlarged part of the middle area I;

[0049] Figure 5A for Figure 3A Enlarged part of middle area II;

[0050] Figure 5B for Figure 3B Enlarged part of middle area II;

[0051] Figure 5C for Figure 3C Enlarged part of middle area II;

[0052] Fig. 6A Comparison of reflection coefficients at x=50m under three imaging conditions.

[0053] Figure 6B Comparison of reflection coefficients under three imaging conditions at x=900m.

[0054] Figure 7 is the velocity model of actual seismic data;

[0055] Fig. 8A is the reverse time migration result under the cross-correlation imaging condition;

[0056] Figure 8B This is the reverse time migration result of the traditional excitation amplitude imaging condition;

[0057] Figure 8C The reverse time migration result of the improved excitation amplitude imaging condition of the present invention;

[0058] Fig.9A for Fig. 8A The enlarged part of the middle area I;

[0059] Fig. 9B for Figure 8B The enlarged part of the middle area I;

[0060] Fig. 9C for Figure 8CThe enlarged part of the middle area I;

[0061] Fig.10 FIG. 4 is a diagram showing the speedup ratio of the present invention relative to the cross-correlation imaging condition. DETAILED DESCRIPTION

[0062] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the invention claimed for protection, but merely represents the selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0063] Example

[0064] A reverse time migration imaging method based on solving the adaptive eikonal equation, such as Figure 1 As shown, the following steps are included:

[0065] S1. Input model and source parameters

[0066] This embodiment sets the Marmousi model, such as Figure 2 As shown, the model has a grid size of 234×663, a spatial sampling interval of 6.25m, a total of 100 shots, each shot interval of 37.5m, the source is a Ricker wavelet with a main frequency of 25HZ, the total simulation time is 3s, and the time sampling interval is 0.5ms.

[0067] S2. Obtain reference solutions for the first trip

[0068] Specifically, it includes: using plane wave operator, spherical wave operator and refraction wave operator to calculate the first arrival travel time in each direction and iterate until the result is stable, and then convert it into the first arrival travel time reference solution;

[0069] The calculation formula of the plane wave operator is:

[0070]

[0071] Where, T i,j represents the first arrival time at (i, j), T h and T v Variables set to facilitate the following expressions;

[0072]

[0073] Where, d x and d z are the spatial sampling intervals of the model grid in the horizontal and vertical directions, S i-1,j-1 represents the slowness at (i-1, j-1), a, b and c are variables set for the convenience of expression, the same below;

[0074] The calculation formula of the spherical wave operator is:

[0075]

[0076] In the formula, T0 and S0 are known solutions, T 0,x and T 0,z are the derivatives of T0 with respect to x and z respectively;

[0077] The calculation formula of the refraction wave operator is:

[0078]

[0079] In the formula, min() means finding the minimum value from different variables; T rh and T rv They represent the travel time from point (i-1, j) and (i, j-1) to point (i, j) respectively;

[0080] The travel time of each point is calculated using the plane wave operator, spherical wave operator and refraction wave operator, and the minimum value is selected for iterative update until the initial arrival travel time is stable, and T is obtained. AFDM .

[0081] First arrival travel time T calculated based on adaptive finite difference method AFDM In order not to ignore the influence of the first arrival wavefront tail part and the observation system layout, the first arrival travel time reference solution T is obtained by the following formula re for:

[0082]

[0083] In the formula, f m is the center frequency of the source wavelet. The factor a is set to reduce the error caused by the arrangement of the observation system, etc. The value range of a is [1,1.25].

[0084] S3, input the forward extension of the source, and save the forward wave field corresponding to the first arrival travel time

[0085] The seismic source is input into the second-order velocity-stress acoustic wave equation for forward modeling. The formula is as follows:

[0086]

[0087] Where P s is the forward wave field data at any time, v is the velocity vector, x and z are the x and z directions in the 2D plane respectively, t is the time, f(t) is the earthquake source, x s and z s Respectively represent the x-direction and z-direction positions of the earthquake source, δ(xx s ,zz s )f(t) means setting the source at (x s ,z s )

[0088] In the process of forward wave field continuation, the first arrival travel time reference solution obtained in step S2 is used as the excitation time, and the forward wave field corresponding to the excitation time is saved as the excitation amplitude, and the formula is as follows:

[0089] S max (x,z;T)=max{S(x,z;T)},T∈[T re -dt,T re +dt];

[0090] Where S(x,z;T) is the forward wave field corresponding to the excitation time, max{} is the function that takes the maximum value, and S max (x,z;T) is the excitation amplitude.

[0091] S4. Input the earthquake record and extend it backwards to save the reverse propagation wave field corresponding to the initial travel time.

[0092] In the second-order velocity-stress acoustic wave equation, the seismic record is input to simulate the back propagation wave field extension, and the formula is as follows:

[0093]

[0094] Where P r is the back propagation wave field data at any time, u is the seismic record, which is used as the boundary condition, and the back propagation wave field corresponding to the excitation time in step S3 is saved as R(x, z; T).

[0095] S5. Imaging using imaging conditions

[0096] Optionally, in step S5, in the step of applying imaging conditions to image, the forward propagation wave field S at the excitation time saved in step S3 and step S4 is max (x, z; T) and the reverse wave field R(x, z; T) apply the excitation amplitude imaging condition, the formula is as follows:

[0097] I(x,z;s j )=R(x,z;T;s j )·Smax (x,z;T;s j );

[0098] In the formula, s j represents the jth earthquake source, I(x,z;s j ) is the offset image of each shot, S max (x,z;T;s j ) is the forward wave field corresponding to the excitation time in each shot, R(x,z; T; s j ) is the back propagation wave field corresponding to the excitation moment in each shot.

[0099] Figure 3A , Figure 3B and Figure 3C They are respectively the offset images of the cross-correlation imaging condition, the traditional excitation amplitude imaging condition and the improved excitation amplitude imaging condition of the present invention. Figure 4B The arrows in the figure indicate the distortion areas in the application of the traditional excitation amplitude imaging conditions. Figure 4C The improved excitation amplitude imaging condition of the present invention greatly improves the signal-to-noise ratio and solves the problem stably and effectively. Figure 4A The imaging effect is comparable. Figure 5A , Figure 5B and Figure 5C The imaging effect of the present invention on the deep structure is better than that of the other two methods as indicated by the dotted line.

[0100] Fig. 6A and Figure 6B In the figure, the black curve is the cross-correlation imaging condition, the blue curve is the excitation amplitude imaging condition, and the red curve is the improved excitation amplitude imaging condition of the present invention. It can be clearly seen that the red curve is closer to the black curve, and the dotted line frame shows the part of the blue curve where there is obvious noise.

[0101] Application Examples

[0102] The method of the present invention is applied to actual seismic data, and the specific implementation steps are as follows:

[0103] S1. Input model and source parameters

[0104] like Figure 7 As shown, the profile is divided into a grid size of 1561×801, the spatial sampling interval is 10m, a total of 88 shots are unevenly distributed on the surface, the recording time of each shot is 6s, the time sampling interval is 0.5ms, and the Ricker wavelet with a main frequency of 25HZ is set as the earthquake source.

[0105] S2-S5 are the same as in the embodiment.

[0106] Fig. 8A , Figure 8B and Figure 8CThey are images under the cross-correlation imaging condition, the traditional excitation amplitude imaging condition and the improved excitation amplitude imaging condition of the present invention. Fig. 8A and Figure 8C The phase axis is obviously higher than Figure 8B More continuous, clearer, and with a wider identifiable range, but due to the high amplitude of the cross-correlation imaging condition, Fig. 8A The event energy in Figure 8C .

[0107] To provide a clearer explanation, this application example will Fig. 8A , Figure 8B and Figure 8C The area I in the figure is enlarged, and we get Fig.9A , Fig. 9B and Fig. 9C .contrast Fig.9A and Fig. 9C It can be found that the imaging results of the method of the present invention are comparable to those of the cross-correlation imaging method. Fig. 9B In comparison, the method of the present invention produces more continuous event axes and higher resolution in deeper layers, thereby proving the effectiveness of the method of the present invention.

[0108] Fig.10 The calculation speed improvement ratio of the present invention relative to the cross-correlation imaging condition is demonstrated. It can be observed that when the model size exceeds 400×400, the speed improvement increases significantly, and the maximum speed improvement can be as high as 129.924%. The calculation efficiency of the present invention is much higher than that of the cross-correlation imaging condition.

[0109] For prestack reverse time migration, the cross-correlation imaging condition requires the storage of all forward and reverse wavefields, which places high demands on hardware, especially when processing large-scale 3D data. For this reason, the excitation amplitude imaging condition is used as an efficient alternative. Although the excitation amplitude imaging condition is very efficient in computation, the extracted maximum amplitude may become chaotic in areas with significant lateral velocity variations and large incident angles, resulting in artifacts and distortion in the migration image. During the stacking process, this chaos causes the interference of the phase axis to destruct, which ultimately affects the quality of the migration image.

[0110] The present invention improves the traditional excitation amplitude imaging conditions, adopts an adaptive finite difference method to solve the eikonal equation to calculate the first arrival travel time, and accurately simulates the first arrival travel time of near, medium and far migration. The travel time calculated in this way is continuous and stable, and the travel time is used as the excitation time, and the corresponding forward wave field and reverse wave field are also saved continuously and stably. Then the imaging conditions are applied to obtain high-fidelity migration images, while also maintaining the advantage of efficient calculation.

[0111] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A reverse time migration imaging method based on solving an adaptive eikonal equation, characterized in that: The steps include: S1, input model and source parameters; S2, obtaining a reference solution for the initial trip; S3, input the forward extension of the earthquake source and save the forward wave field at the time of the first arrival travel; S4, reverse continuation of input seismic records, and saving the reverse propagation wave field corresponding to the first arrival travel time; S5, applying imaging conditions to form images; In step S2, in the step of obtaining the first arrival travel time reference solution, the first arrival travel time in each direction is calculated using the plane wave operator, the spherical wave operator and the refraction wave operator respectively and iterated until the result is stable and converted into the first arrival travel time reference solution; The calculation formula of the plane wave operator is: Where, T i,j represents the first arrival time at (i, j), T h and T v Variables set for convenience of expression; Where, d x and d z are the spatial sampling intervals of the model grid in the horizontal and vertical directions, S i-1,j-1 represents the slowness at (i-1, j-1), a, b and c are variables set for the convenience of expression; The calculation formula of the spherical wave operator is: In the formula, T0 and S0 are known solutions, T 0,x and T 0,z are the derivatives of T0 with respect to x and z respectively; The calculation formula of the refraction wave operator is: In the formula, min() means finding the minimum value from different variables, T rh and T rv They represent the travel time from point (i-1, j) and (i, j-1) to point (i, j) respectively; The travel time of each point is calculated using the plane wave operator, spherical wave operator and refraction wave operator, and the minimum value is selected for iterative update until the initial arrival travel time is stable, and T is obtained. AFDM ; T AFDM The reference solution T for the initial travel time is obtained by the following formula: re : In the formula, f m is the center frequency of the source wavelet, and the value range of factor a is [1,1.25]; In step S3, the forward continuation of the input source is performed, and in the step of saving the forward wave field corresponding to the first arrival travel time, the source is input into the second-order velocity-stress acoustic wave equation for forward simulation, and the formula is as follows: Where P s is the forward wave field data at any time, v is the velocity vector, x and z are the x and z directions in the 2D plane respectively, t is the time, f(t) is the earthquake source, x s and z s Respectively represent the x-direction and z-direction positions of the earthquake source, δ(xx s ,zz s )f(t) means setting the source at (x s ,z s ) In the process of forward wave field continuation, the first arrival travel time reference solution obtained in step S2 is used as the excitation time, and the forward wave field corresponding to the excitation time is saved as the excitation amplitude, and the formula is as follows: S max (x,z;T)=max{S(x,z;T)},T∈[T re -dt,T re +dt]; Where S(x,z;T) is the corresponding forward wave field within the range of T, max{} is the function that takes the maximum value, and S max (x, z; T) is the excitation amplitude; In step S4, input the reverse continuation of the seismic record, and save the reverse propagation wave field corresponding to the first arrival travel time. In the second-order velocity-stress acoustic wave equation, input the reverse propagation wave field continuation of the seismic record simulation, and the formula is as follows: Where P r is the back propagation wave field data at any time, u is the seismic record, as the boundary condition, and the back propagation wave field corresponding to the excitation time in step S3 is saved as R(x,z;T); In step S5, in the step of applying imaging conditions, the forward propagation wave field S at the excitation time saved in step S3 and step S4 is max (x, z; T) and the back propagation field R(x, z; T) are applied to the excitation amplitude imaging condition to obtain the migration image, the formula is as follows: I(x,z;s j )=R(x,z;T;s j )·S max (x,z;T;s j ); In the formula, s j represents the jth earthquake source, I(x,z;s j ) is the offset image of each shot, S max (x,z;T;s j ) is the forward wave field corresponding to the excitation time in each shot, R(x,z; T; s j ) is the back propagation wave field corresponding to the excitation moment in each shot.

Images