Wave front rapid marching algorithm for anisotropic media

Abstract

The application discloses a wave front fast advancing algorithm suitable for anisotropic media, and is implemented according to the following steps: step 1, rewriting the eikonal equation, discretizing the eikonal equation by using the fast advancing algorithm, and initializing a source point as a known point; step 2, obtaining a function relation about an incident angle, an incident phase angle and a stratum dip angle by using the known point; step 3, obtaining a travel time of a to-be-solved point by using a phase angle and a phase velocity of the to-be-solved point, and calculating travel times of all grid nodes; and step 4, after the wave front time field is solved by using steps 1 and 3, determining a phase velocity field, calculating a group velocity field by using the phase velocity field, and finally completing ray tracing by using the group velocity. The application shapes the exit direction and exit velocity of the to-be-solved point as known quantities by using the snell equation, and then obtains high-precision travel times of all grid nodes, determines more accurate time field gradients of all grid nodes to determine group velocity vectors, and performs ray tracing.

Description

Technical Field

[0001] This invention belongs to the technical field of seismic wave propagation methods, specifically relating to a rapid wavefront propagation algorithm applicable to anisotropic media. Background Technology

[0002] Fast and accurate travel time calculation is crucial for the success of many seismic applications, including static and time-of-flight correction, travel time tomography, microseismic source location, and pre-stack migration. Therefore, travel time calculation algorithms must take into account the seismic anisotropy of high-resolution subsurface imaging.

[0003] Currently, the two most commonly used algorithms for calculating seismic wave travel times in anisotropic media are the Fast Progress Method (FMM) and the Fast Scan Method (FSM). FSM updates grid-point travel time information by locally solving equations and continuously scanning the model to stabilize the travel time field. The local solution is typically an upwind interpolation scheme. However, this method requires more iterations in complex models to obtain high-precision travel times, making it more time-consuming and computationally inefficient than FMM. FMM, on the other hand, utilizes the upwind difference scheme to determine the causal relationship of wavefront evolution while employing heap sorting techniques to reduce the computational complexity from O(N) of the wavefront extension method. 3 The computational efficiency is reduced to O(NlogN), where N is the number of nodes. However, in anisotropic media, the phase velocity is uncertain due to the change in velocity with the direction of propagation. Therefore, the original fast propagation algorithm cannot be reliably used to solve anisotropic equations. In recent years, many effective methods have been developed for solving anisotropic equations using the Fast Mine Model (FMM), but their computational efficiency and accuracy are lower than those in isotropic media. Summary of the Invention

[0004] The purpose of this invention is to provide a fast wavefront propagation algorithm suitable for anisotropic media, and to solve the problem of difficulty in solving the phase velocity in existing fast propagation algorithms in anisotropic media.

[0005] The technical solution adopted in this invention is a wavefront fast propagation algorithm applicable to anisotropic media, which is implemented according to the following steps:

[0006] Step 1: Rewrite the equations, discretize them using the fast propagation algorithm, and initialize the source points as known points;

[0007] Step 2: Use known points to obtain the functional relationship between the incident angle, the incident phase angle, and the dip angle of the formation;

[0008] Step 3: Using the phase angle and phase velocity of the point to be determined in Step 2, obtain the travel time of the point to be determined, and calculate the travel time of all grid nodes.

[0009] Step 4: After solving the wavefront time field using Step 1 and Step 3, determine its phase velocity field, then use the phase velocity field to calculate the group velocity field, and finally use the group velocity to complete ray tracing.

[0010] The invention is further characterized by:

[0011] In step 1, the functional equation is rewritten as:

[0012]

[0013] Where τ is the wavefront time and θ is the phase angle;

[0014] Equation (1) is discretized using a fast propagation algorithm, and the discretized equation (1) is further simplified to:

[0015]

[0016]

[0017] in

[0018] In the formula, It is the second-order difference operator τ in the x-direction during the travel of the point to be found, where Δx is the grid spacing in the x-direction. It is the second-order difference operator τ in the z-direction when the point to be determined travels, where Δz is the grid spacing in the z-direction.

[0019] In step 1, the fast-advance algorithm divides the entire computational domain into an upwind zone, a wavefront zone, and a downwind zone. The wave propagates from the upwind zone through the wavefront zone to the downwind zone. The upwind zone is the known area where the travel time has been solved, the wavefront zone is the area being solved, and the downwind zone is the unknown area that has not yet been solved. Within the narrow band of the wavefront zone, the minimum travel time point is found, and its attribute is changed to a known point while continuing to expand downwind until all grid points are in the upwind zone.

[0020] Step 2 is as follows:

[0021] During the propagation process, the travel time of point a is known. Using equation (3), the travel time of this point in the x and z directions can be obtained. and Using the following equation (4), we can know the incident direction β from point a to point b. a :

[0022]

[0023] In addition, β a Let θ be the incident direction at point b in the grid. inc Let b be the incident phase angle, and α be the incident phase angle. b Let the interface tilt angle be the factor, and the following relationship exists between the three:

[0024] θ inc =β a -α b (5)

[0025] By combining formulas (4), (5) and (2), we can obtain the velocity v from point a to point b. p (θ inc ).

[0026] The propagation of rays in TTI media follows the following equation (6):

[0027]

[0028] Where p is the ray parameter, θ inc θ b These are the incident and exit angles, v p (θ a ), v p (θ b θ represents the incident and exit velocities, respectively, and the exit direction θ at point b is... b and the propagation speed v p (θ b As shown in equations (7) and (8):

[0029]

[0030] v p (θ b ) = v n (1+ε b sin 4 θ b +δ b sin 2 θ b cos 2 θ b (8)

[0031] Among them, v n ε is the formation normal velocity, and δ and ε are Thomsen parameters.

[0032] The specific process in step 4 is as follows:

[0033] Using formulas (2) and (8), the travel time of all grid points is obtained, thus determining the wavefront time field. The entire phase velocity field can be determined using gradient calculations within the wavefront time field.

[0034]

[0035]

[0036] in, These represent the components of the wavefront gradient in the formation normal and the formation dip directions, respectively.

[0037] Group velocity field is obtained by utilizing the relationship between group velocity and phase velocity:

[0038]

[0039] in, and The unit vector representing the tangential and normal directions of the formation, v g and v p (θ) represents the group velocity and phase velocity, respectively. In TTI medium, the phase angle is obtained by coordinate rotation transformation of the wavefront gradient direction. and The group velocity vector can be expressed using the tilt angle, therefore equation (10) can be expressed as:

[0040]

[0041] Where α is the dip angle of the strata.

[0042] Finally, starting from the receiving point, the path is traced in the opposite direction of the group velocity vector.

[0043] The beneficial effects of this invention are as follows: This invention provides a fast wavefront advancement algorithm suitable for anisotropic media. Addressing the problem of difficulty in obtaining accurate phase slowness at the target point in conventional FMMs in anisotropic media, a fast advancement method based on Snell and FMM phase angle prediction (PAP-FMM) is proposed. By using the Snell equation to shape the exit direction and exit velocity of the target point into known quantities, high-precision travel time of each grid point is obtained, thereby determining a more accurate time field gradient for each grid point to determine the group velocity vector for ray tracing. This method has the advantages of high computational accuracy, speed, and stability. Attached Figure Description

[0044] Figure 1 This is a graphical description of the wavefront rapid propagation algorithm applicable to anisotropic media according to the present invention;

[0045] Figure 2 This is a schematic diagram illustrating the change in propagation direction during wavefront expansion in this invention;

[0046] Figure 3 This is a flowchart of the phase velocity estimation process of the present invention;

[0047] Figure 4 This is a schematic diagram of phase velocity and group velocity ray tracing in this invention;

[0048] Figure 5 This is a schematic diagram of anisotropic medium wavefront tracking according to the present invention;

[0049] Figure 6 These are schematic diagrams of ray tracing in this invention: (a) layered model, (b) TTI model, and (c) Marmousi model.

[0050] Figure 7 This is a comparison chart of the errors of the present invention with conventional FMM algorithms, test firing methods, and theoretical values. Detailed Implementation

[0051] The present invention will now be described in detail with reference to specific embodiments and accompanying drawings.

[0052] Fast wavefront propagation algorithms applicable to anisotropic media, such as Figure 1 As shown, please follow these steps:

[0053] Step 1: Rewrite the equations, discretize them using the fast propagation algorithm, and initialize the source points as known points;

[0054] Step 2: Use known points to obtain the functional relationship between the incident angle, the incident phase angle, and the dip angle of the formation;

[0055] Step 3: Using the phase angle and phase velocity of the point to be determined estimated in Step 2, substitute them into the discretization equation in Step 1 to obtain the travel time of the point to be determined. Then, select the minimum travel time point from the narrow band and change the attribute of the minimum point to the known point and place it in the upwind area. Then, determine whether the narrow band is empty. If the narrow band is empty, end the calculation. Otherwise, repeat Step 2 and Step 3 until the travel time of all grid nodes has been calculated.

[0056] Step 4: After solving the wavefront time field using Step 1 and Step 3, determine its phase velocity field, then use the phase velocity field to calculate the group velocity field, and finally use the group velocity to complete ray tracing.

[0057] Example 1

[0058] The specific process of step 1 is as follows:

[0059] The equation can be rewritten as:

[0060]

[0061] Where τ is the wavefront time, and θ is the angle between the wavefront normal direction and the formation normal direction, called the phase angle;

[0062] Then, the fast propagation algorithm is used to discretize equation (1), and the discretized equation (1) is further simplified to:

[0063]

[0064]

[0065] in

[0066] In the formula, It is the second-order difference operator τ in the x-direction during the travel of the point to be found, where Δx is the grid spacing in the x-direction. It is the second-order difference operator τ in the z-direction when the point to be determined travels, where Δz is the grid spacing in the z-direction.

[0067] like Figure 1 As shown, the fast-progression algorithm divides the entire computational domain into three regions: the upwind region, the wavefront region, and the downwind region. The wave propagates from the upwind region through the wavefront region to the downwind region. The upwind region is the known region where the departure time has been solved, the wavefront region is the region being solved, and the downwind region is the unknown region that has not yet been solved. Within the narrow band of the wavefront region, the minimum travel time point is found, and its attributes are changed to known points while continuing to expand downwind region until all grid points are in the upwind region.

[0068] Example 2

[0069] like Figure 2 As shown, the specific process of step 2 is as follows:

[0070] During the wavefront propagation process in step 1, since the velocities of the known point a and the point b to be determined are inconsistent, their directions will also change during propagation. Given the travel time of point a during propagation, equation (3) can be used to obtain the x-direction and z-direction of this point. and Using equation (4), we can determine the incident direction β from point a to point b. a :

[0071]

[0072] In addition, β a Let θ be the incident direction at point b in the grid. inc Let b be the incident phase angle, and α be the incident phase angle. b Let the interface tilt angle be the factor, and the following relationship exists between the three:

[0073] θ inc =β a -α b (5)

[0074] By combining formulas (4), (5) and (2), we can obtain the velocity v from point a to point b. p (θ inc ).

[0075] Example 3

[0076] The specific process of step 3 is as follows:

[0077] When the propagation speed in step 2 reaches the point b to be determined, it is unknown how its propagation direction and speed will change. If we still use the speed of the known point a for calculation, it will lead to a large error in the finite difference approximation solution. In addition, this error spreads to the entire computational domain, causing the accuracy of the travel time solution to be much lower than expected. A particularly attractive method to solve the above problem is to use the Snell equation to predict the propagation speed of point b based on the propagation speed of point a.

[0078] Rays propagate through TTI media according to formula (6):

[0079]

[0080] Where p is the ray parameter, θ inc θ b These are the incident and exit angles, v p (θ a ), v p (θ b The incident and exit velocities are respectively, and the direction θ from point a to point b can be determined by step 2. inc and velocity v p (θ inc The emission direction θ at point b b and the propagation speed v p (θ b As shown in equations (7) and (8):

[0081] For the outgoing phase angle θ b :

[0082]

[0083] Launch velocity v p (θ b ):

[0084] v p (θ b ) = v n (1+ε b sin 4 θ b +δ b sin 2 θ b cos 2 θ b (8)

[0085] Among them, v n ε is the formation normal velocity, and δ and ε are Thomsen parameters;

[0086] During the solution process, the output phase angle must satisfy 0°≤θ bWithin a range of ≤90°, the launch velocity v p (θ b Solving using equation (8) yields the result. Finally, to obtain a more accurate travel time solution, we need to continuously iterate the estimated launch velocity until it meets the set error range, and then make the next estimate, and so on. The specific process is as follows: Figure 3 As shown.

[0087] Example 4

[0088] The specific process of step 4 is as follows:

[0089] exist Figure 4 In the ray tracing diagram of phase velocity and group velocity shown, a slight difference in their deflection directions can be observed. This is because the directions of phase velocity and group velocity are not consistent in anisotropic media. Figure 5 As shown, the ray direction is consistent with the group velocity vector direction. Therefore, using formulas (2) and (8), the travel time of all grid points is obtained, thus determining the wavefront time field. The entire phase velocity field can be determined using gradient calculations within the wavefront time field.

[0090]

[0091] in, These represent the components of the wavefront gradient in the formation normal and the formation dip directions, respectively.

[0092] Group velocity field is obtained by utilizing the relationship between group velocity and phase velocity:

[0093]

[0094] in and The unit vector representing the tangential and normal directions of the formation, v g and v p (θ) represents the group velocity and phase velocity, respectively. In TTI medium, the phase angle is obtained by coordinate rotation transformation of the wavefront gradient direction. and The group velocity vector can be expressed using the tilt angle, therefore equation (10) can be expressed as:

[0095]

[0096] Where α is the dip angle of the strata;

[0097] Finally, starting from the receiving point, the path is traced against the direction of the group velocity vector.

[0098] This invention relates to a fast wavefront propagation algorithm for anisotropic media. Numerical simulations were performed in Fortran 2012, verifying the effectiveness of the method in layered models, TTI models, and Marmousi models. Figure 6 As shown, from Figure 6 As can be seen, the ray deflection position coincides with the interface position, and the density and spacing of the isochronous lines also correspond to the velocity change characteristics. Therefore, the improved fast-advancement algorithm is an effective, stable, and accurate ray tracing method.

[0099] Error analysis was performed on the results of this invention, such as... Figure 7 As shown, from Figure 7 As can be seen, the conventional FMM has a large error, with a maximum relative error of 0.559%. For the test firing method, the tracking results near the epicenter point show that its algorithm is extremely unstable, with a maximum relative error of nearly 6.7%. The improved fast propulsion algorithm has a smaller error, and the error is very stable with all errors less than 0.16%, which has the advantages of high calculation accuracy and fast stability.

Claims

1. A fast wavefront propagation algorithm applicable to anisotropic media, characterized in that, The specific steps are as follows: Step 1: Rewrite the equations, discretize them using the fast propagation algorithm, and initialize the source points as known points; Step 2: Use known points to obtain the functional relationship between the incident angle, the incident phase angle, and the dip angle of the formation; Step 3: Using the phase angle and phase velocity of the point to be determined in Step 2, obtain the travel time of the point to be determined, and calculate the travel time of all grid nodes. Step 4: After solving the wavefront time field using Step 1 and Step 3, determine its phase velocity field, then use the phase velocity field to calculate the group velocity field, and finally use the group velocity to complete ray tracing. The specific process of step 2 is as follows: During propagation, given the travel time of point a, the x- and z-direction coordinates of this point are... and , Given the wavefront time, the incident direction of the propagation from point a to point b can be determined using the following equation (4). : （4） also, The incident direction at point b in the grid, Let b be the incident phase angle. Let the interface tilt angle be the factor, and the following relationship exists between the three: （5） Get the speed at which the propagation from point a to point b is obtained. ; The propagation of rays in TTI media follows the following equation (6): (6) Where p is the ray parameter, , These are the incident and exit angles, respectively. , These are the incident and exit velocities, and the exit direction at point b. and the speed of transmission As shown in equations (7) and (8): （7） （8） Among them, v n ε is the velocity in the normal direction of the formation, δ and ε are Thomsen parameters, and x is the sine value corresponding to the incident angle; The specific process of step 4 is as follows: Using formula (8), the travel time of all grid points is obtained, thus determining the wavefront time field. The entire phase velocity field can be determined using gradient calculation within the wavefront time field. （9） in, , These represent the components of the wavefront gradient in the formation normal and the formation dip directions, respectively. Group velocity field is obtained by utilizing the relationship between group velocity and phase velocity: （10） in, and A unit vector representing the tangential and normal directions of the formation. and Let these represent the group velocity and phase velocity, respectively. In TTI media, the phase angle is obtained by transforming the wavefront gradient direction using coordinate rotation. and The group velocity vector can be expressed using the tilt angle, therefore equation (10) can be expressed as: (11） in, The dip angle of the strata; Finally, starting from the receiving point, the path is traced against the direction of the group velocity vector.

2. The wavefront fast propagation algorithm for anisotropic media according to claim 1, characterized in that, In step 1, the functional equation is rewritten as follows: （1） in, For wavefront time, Phase angle; Equation (1) is discretized using a fast propagation algorithm, and the discretized equation (1) is further simplified to: （2） in (3) In the formula, , Is it a travel spot to be sought? exist Second-order difference operators in direction, yes Grid spacing in the direction, , Is it a travel spot to be sought? Second-order difference operator in the z-direction, It is the grid spacing in the z-direction.

3. The wavefront fast propagation algorithm for anisotropic media according to claim 1, characterized in that, In step 1, the fast-advance algorithm divides the entire computational domain into an upwind zone, a wavefront zone, and a downwind zone. The wave propagates from the upwind zone through the wavefront zone to the downwind zone. The upwind zone is the known area where the departure time has been solved, the wavefront zone is the area being solved, and the downwind zone is the unknown area that has not yet been solved. Within the narrow band of the wavefront zone, the minimum travel time point is found, and its attribute is changed to a known point while continuing to expand to the downwind zone until all grid points are in the upwind zone.

Images