CN110687598A - A method and device for accelerating microseismic numerical simulation - Google Patents

Abstract

The invention relates to a microseismic numerical simulation method and device, belonging to the field of microseismic monitoring, in particular to a method and device for accelerating microseismic numerical simulation. The invention combines the characteristics of the number of sources and receiving geophones in microseismic monitoring, according to the principle of source-detection interchange, and on the basis of the prior art, the original source point is used as the receiving point, the original receiving point is used as the source point, the seismic wave field is calculated, and the seismic wave field can be quickly obtained. Compared with the prior art, the wave field information of each receiving point in the microseismic numerical simulation greatly improves the calculation efficiency.

Description

A method and device for accelerating microseismic numerical simulation

technical field

The invention relates to a microseismic numerical simulation method and device, belonging to the field of microseismic monitoring, in particular to a method and device for accelerating microseismic numerical simulation.

Background technique

Microseismic monitoring technology is a geophysical technology that has been widely used in oil and coal fields in recent years, and is widely used in fracturing monitoring, rock burst monitoring and other fields. Microseismic monitoring technology mainly involves two major technical problems: forward modeling and inversion. Microseismic monitoring forward modeling problem, that is, given parameters such as velocity model, observation system, source coordinates, etc., to solve the information of microseismic wave field (dynamics) and seismic wave travel time (kinematics), and verify the correctness of microseismic monitoring and observation system design and positioning method It is the basic problem in microseismic monitoring and the basis of inversion. The inversion problem of microseismic monitoring is the inverse process of forward modeling, that is, the process of obtaining information such as velocity model, source parameters, etc., knowing the observation system, seismic wave field and travel time.

Microseismic forward modeling methods mainly include: finite difference method based on wave equation theory, finite element method, pseudospectral method, wavefront fast advancing method, etc., boundary integral method based on integral equation, boundary element method, etc. Shooting, bending, etc.

The basic process of the forward modeling algorithm for conventional microseismic monitoring is as follows: 1. Given the 3D grid velocity model of the simulation area, the coordinates of the source and receiver points, and other parameters; 2. Calculate the seismic wave field (motion) of each hypocenter in the 3D model space. 3. Extract the information at each detection point in the seismic wave field. Since the forward modeling of microseismic monitoring is a three-dimensional problem and the model space is large, it is a very time-consuming process to calculate the seismic wave field in the three-dimensional model space. The main influencing factors include the size of the model space and the size of grid division. If the model space is large and the grid is finely divided, it will take tens or even hundreds of minutes to calculate the three-dimensional seismic wave field of a hypocenter. Assuming that there are M microseismic sources, it takes N minutes to calculate one source, and the time to complete the forward modeling of all source points is roughly M×N minutes. Therefore, the second step above is the main factor restricting the efficiency of microseismic forward modeling.

SUMMARY OF THE INVENTION

The invention mainly solves the technical problem of slow microseismic numerical simulation (forward modeling) existing in the prior art, and provides a method and device for accelerating the microseismic numerical simulation. The method and device combine the characteristics of the number of sources and detection points in microseismic monitoring, and according to the principle of source-detection interchange, the wave field information of each detection point in the microseismic numerical simulation can be quickly obtained, and the calculation efficiency is greatly improved compared with the prior art.

The above-mentioned technical problems of the present invention are mainly solved by the following technical solutions:

A method for accelerating microseismic numerical simulation, which specifically includes the following steps:

In the given step of forward modeling simulation parameters, the grid velocity model is established according to the geological conditions of the object to be simulated, and the coordinates of the detection point are specified;

In the wave field calculation step, the seismic wave field calculation method is determined according to the number of source points Ns and the number of detection points Nr. If Ns>Nr, perform source detection interchange and wave field calculation, that is, take each receiver point as the source point, and calculate the seismic wavefield of each receiver point as the source point; if Ns<Nr, according to the conventional method, calculate each The seismic wave field of a hypocenter. In microseismic monitoring, it generally belongs to the former.

In the wave field extraction step, based on the three-dimensional seismic wave field obtained by the above calculation, the seismic wave field of each source point is extracted, and the forward modeling ends.

Preferably, the step of setting forward simulation parameters executes the following sub-steps:

Establish an XYZ three-dimensional Cartesian coordinate system according to the area to be simulated; perform grid division on the area to be simulated based on the three-dimensional Cartesian coordinate system; based on the grid division result, assign velocity values to the grid points to establish a velocity model; specify each detector The spatial location of the point. Among them, the grid velocity model, the spatial coordinates of the detection point and the microseismic source point are respectively recorded as:

Grid velocity model: vi _,j,k (1≤i≤N _x , 1≤j≤N _y , 1≤k≤N _z )

Detecting point coordinates: (xr _i ,yr _i ,zr _ri )(1≤i≤N _r )

The coordinates of the epicenter: (xs _i , ys _i , zs _ri ) (1≤i≤N _s )

Variable description: v _{i, j, k} is the velocity corresponding to the grid node (i, j, k), N _x is the number of grids in the X direction of the model, N _y is the number of grids in the Y direction of the model, and N _z is the Z direction of the model Number of grids, xri is the X coordinate of the _ith detection point, yri is the Y coordinate of the _ith detection point, zr _i is the Z coordinate of the ith detection point, N _r is the number of detection points, xs _i is the X coordinate of the ith hypocenter point, ys _i is the Y coordinate of the ith hypocenter point, zs _i is the Z coordinate of the ith hypocenter point, and _Ns is the number of hypocenter points.

Preferably, the wave field calculation module executes the following sub-steps:

First, according to the number of source points Ns and the number of receiver points Nr, determine the seismic wave field calculation method. If Ns<Nr, traverse the source points as usual; if Ns>Nr, traverse each receiver point in turn, take the receiver point as the source point, and mark the original source point as the receiver point.

Second, calculate the seismic wave field of the current hypocenter point under the known velocity model. Taking the fast forward algorithm as an example, solve the following equation of the program function in the two-dimensional rectangular coordinate system (x-z):

where t is the travel time of the seismic wave, t(x, z) represents the travel time of the seismic wave at the coordinate point (x, z), s is the slowness, the reciprocal of the velocity, and s ² (x, z) represents the travel time of the seismic wave at the coordinate point (x, z) ), the square of the slowness value, x and z represent the two coordinate axes of the two-dimensional rectangular coordinate system.

The upper wind difference method is used to replace the partial differential format in the functional equation of the above equation, that is, the upper wind difference method is used to discretize the above equation, and its discrete expression can be expressed as follows:

D is the difference operator, i and j are the grid node index numbers in the x and z directions, respectively, max is the maximum value function, and min is the minimum value function. The above formula can be simplified to:

表示走时t在点(i,j)处x方向上的向前、向后的差分算子，

表示走时t在点(i,j)处z方向上的向前、向后的差分算子。

represents the forward and backward difference operators of travel time t in the x direction at point (i, j),

Represents the forward and backward difference operators in the z direction at point (i, j) for travel time t.

In order to ensure the calculation accuracy, the second-order upwind difference operator is used. which is:

where Δx and Δz are the grid sizes in the X and Z directions, respectively, and t _i,j represents the travel time at point (i, j).

Finally, using the travel time values at the eight nodes (i,j±1), (i,j±2), (i±1,j), (i±2,j) around the target point (i,j), There are 6 combinations to calculate the travel time value of point (i, j), and the second-order difference calculation formula of the 6 combinations is as follows:

In the formula, t is the seismic wave travel time, s is the slowness of the seismic wave, and h=Δx=Δz is the square grid size.

Using the above formula to calculate multiple traveltimes, according to the basic principle of narrow-band expansion, select the minimum traveltime point in the obtained multi-valued solution, which is the target point traveltime. The calculation formula is as follows:

t _i,j = min(t ¹ _i,j ,t ² _i,j ,t ³ _i,j ,t ⁴ _i,j ,t ⁵ _i,j ,t ⁶ _i,j )

After traversing all nodes, the seismic wave time field can be obtained.

Preferably, the wave field extraction step executes the following sub-steps:

Extract the seismic wave field of each detection point, that is, find the time corresponding to the location of the detection point in the time field t(x, z) obtained in the previous step, which is the earthquake travel time from the original source point to the detection point, and the forward modeling ends.

A device for accelerating microseismic numerical simulation, specifically comprising the following steps:

The forward modeling parameter setting device establishes a grid velocity model according to the geological conditions of the object to be simulated, and specifies the coordinates of the detection point.

The wave field calculation device determines the seismic wave field calculation method according to the number of source points Ns and the number of detection points Nr. If Ns>Nr, perform source detection interchange and wave field calculation, that is, take each receiver point as the source point, and calculate the seismic wavefield of each receiver point as the source point; if Ns<Nr, according to the conventional method, calculate each The seismic wave field of a hypocenter. In microseismic monitoring, it generally belongs to the former.

The wave field extraction device extracts the seismic wave field of each detection point based on the three-dimensional seismic wave field obtained by the above calculation, and the forward modeling ends.

Preferably, the forward modeling parameter setting module executes the following sub-steps:

Establish an XYZ three-dimensional rectangular coordinate system according to the area to be simulated; perform grid division on the area to be simulated based on the three-dimensional rectangular coordinate system, and assign velocity values to grid points based on the results of the grid division to establish a velocity model; and specify each detector. The spatial location of the point. Among them, the grid velocity model, the spatial coordinates of the detection point and the potential source point are respectively recorded as:

Grid velocity model: vi _,j,k (1≤i≤N _x , 1≤j≤N _y , 1≤k≤N _z )

Detecting point coordinates: (xr _i ,yr _i ,zr _ri )(1≤i≤N _r )

The coordinates of the epicenter: (xs _i , ys _i , zs _ri ) (1≤i≤N _s )

Variable description: v _{i, j, k} is the velocity corresponding to the grid node (i, j, k), N _x is the number of grids in the X direction of the model, N _y is the number of grids in the Y direction of the model, and N _z is the Z direction of the model Number of grids, xri is the X coordinate of the _ith detection point, yri is the Y coordinate of the _ith detection point, zr _i is the Z coordinate of the ith detection point, N _r is the number of detection points, xs _i is the X coordinate of the ith hypocenter point, ys _i is the Y coordinate of the ith hypocenter point, zs _i is the Z coordinate of the ith hypocenter point, and _Ns is the number of hypocenter points.

Preferably, the wave field calculation module executes the following sub-steps:

First, according to the number of source points Ns and the number of receiver points Nr, determine the seismic wave field calculation method. If Ns>Nr, traverse the source point as usual; if Ns<Nr, traverse each receiver point in turn, take the receiver point as the source point, and mark the original source point as the receiver point.

Second, calculate the seismic wave field of the current hypocenter point under the known velocity model. Taking the fast forward algorithm as an example, solve the following equation of the program function in the two-dimensional rectangular coordinate system (x-z):

Among them, t is the travel time of the seismic wave, s is the slowness, that is, the reciprocal of the velocity, and x and z represent the two coordinate axes of the two-dimensional rectangular coordinate system.

The upper wind difference method is used to replace the partial differential format in the formula function equation, that is, the upper wind difference method is used to discretize the above equation, and its discrete expression can be expressed as follows:

D is the difference operator, i and j are the grid node index numbers in the x and z directions, respectively, max is the maximum value function, and min is the minimum value function. The above formula can be simplified to:

表示走时t在点(i,j)处x方向上的向前、向后的差分算子，

表示走时t在点(i,j)处z方向上的向前、向后的差分算子。

represents the forward and backward difference operators of travel time t in the x direction at point (i, j),

Represents the forward and backward difference operators in the z direction at point (i, j) for travel time t.

In order to ensure the calculation accuracy, the second-order upwind difference operator is used. which is:

where Δx and Δz are the grid sizes in the X and Z directions, respectively.

Finally, using the travel time values at the eight nodes (i,j±1), (i,j±2), (i±1,j), (i±2,j) around the target point (i,j), There are 6 combinations to calculate the travel time value of point (i, j), and the second-order difference calculation formula of the 6 combinations is as follows:

In the formula, t is the seismic wave travel time, s is the slowness of the seismic wave, and h=Δx=Δz is the square grid size.

Using the above formula to calculate multiple traveltimes, according to the basic principle of narrow-band expansion, select the minimum traveltime point in the obtained multi-valued solution, which is the target point traveltime. The calculation formula is as follows:

t _i,j = min(t ¹ _i,j ,t ² _i,j ,t ³ _i,j ,t ⁴ _i,j ,t ⁵ _i,j ,t ⁶ _i,j )

After traversing all nodes, the seismic wave time field can be obtained.

Preferably, the wavefield extraction device performs the following sub-steps:

Extract the seismic wave field of each detection point, that is, find the time corresponding to the location of the detection point in the time field t(x, z) obtained in the previous step, which is the earthquake travel time from the original source point to the detection point, and the forward modeling ends.

Therefore, the present invention has the following advantages:

Combined with the characteristic that the number of source points in microseismic monitoring is far greater than the number of receiver points, and based on the principle of source-detector interchange in seismic wavefield simulation, on the basis of the existing technology, the original source point is used as the receiver point, and the original receiver point is used as the source point. Seismic wavefield, and finally obtain the seismic wavefield information of each source-detection pair. After the exchange of source detection, the efficiency of microseismic numerical simulation is greatly improved.

Description of drawings

FIG. 1 is a flow chart of the present invention.

Fig. 2 is an embodiment velocity model (two-dimensional) and observation system.

Detailed ways

The technical solutions of the present invention will be further described in detail below through embodiments and in conjunction with the accompanying drawings. The protection scope of the present invention is not limited by the following examples.

According to the symmetric sampling principle proposed by Vermeer, the source point (or shot point) and the receiver point are exchanged in seismic exploration, and the same seismic signal can be obtained, that is, excitation at the source point - reception at the receiver point, excitation at the receiver point - source point The signals received by the two are completely equivalent, which is called the principle of source detection interchange (ie, the principle of shot detection interchange).

Based on this principle, the invention firstly establishes a grid velocity model according to the geological conditions of the object to be simulated, and specifies the coordinates of the detection points; then determines the wave field calculation method according to the number of source points and detection points, when the number of hypocenter points is less than the number of detection points , based on the principle of source-detection interchange, take each receiver point as the source point in turn, and calculate the seismic wave field of the model; finally, based on the three-dimensional seismic wave field obtained by the above calculation, extract the seismic wave field of each receiver point (original source point), and forward modeling Finish.

The method for numerical simulation of accelerated microseismic in this embodiment will be further described below with reference to FIG. 1 .

As shown in Figure 1, the present invention adopts following technical scheme:

First, the forward modeling parameters are given, that is, according to the geological conditions of the object to be simulated, a grid velocity model is established, and the coordinates of the detection point are specified.

Specifically, in this embodiment, an XYZ three-dimensional rectangular coordinate system is established according to the area to be simulated; the area to be simulated is divided into grids based on the three-dimensional rectangular coordinate system; based on the results of the grid division, velocity assignments are performed on the grid points, Build a velocity model; specify the spatial location of each pickup point. Among them, the grid velocity model, the spatial coordinates of the detection point and the microseismic source point are respectively recorded as:

Grid velocity model: vi _,j,k (1≤i≤N _x , 1≤j≤N _y , 1≤k≤N _z )

Detecting point coordinates: (xr _i ,yr _i ,zr _ri )(1≤i≤N _r )

The coordinates of the epicenter: (xs _i , ys _i , zs _ri ) (1≤i≤N _s )

Variable description: v _{i, j, k} is the velocity corresponding to the grid node (i, j, k), N _x is the number of grids in the X direction of the model, N _y is the number of grids in the Y direction of the model, and N _z is the Z direction of the model Number of grids, xri is the X coordinate of the _ith detection point, yri is the Y coordinate of the _ith detection point, zr _i is the Z coordinate of the ith detection point, N _r is the number of detection points, xs _i is the X coordinate of the ith hypocenter point, ys _i is the Y coordinate of the ith hypocenter point, zs _i is the Z coordinate of the ith hypocenter point, and _Ns is the number of hypocenter points.

Preferably, the wave field calculation module executes the following sub-steps:

First, according to the number of source points Ns and the number of receiver points Nr, determine the seismic wave field calculation method. If Ns>Nr, traverse the source point as usual; if Ns<Nr, traverse each receiver point in turn, take the receiver point as the source point, and mark the original source point as the receiver point. In microseismic monitoring, it generally belongs to the former.

Second, calculate the seismic wave field of the current hypocenter point under the known velocity model. Taking the fast forward algorithm as an example, solve the following equation of the program function in the two-dimensional rectangular coordinate system (x-z):

Among them, t is the travel time of the seismic wave, s is the slowness, that is, the reciprocal of the velocity, and x and z represent the two coordinate axes of the two-dimensional rectangular coordinate system.

The upper wind difference method is used to replace the partial differential format in the formula function equation, that is, the upper wind difference method is used to discretize the above equation, and its discrete expression can be expressed as follows:

D is the difference operator, i and j are the grid node index numbers in the x and z directions, respectively, max is the maximum value function, and min is the minimum value function. The above formula can be simplified to:

represents the forward and backward difference operators of travel time t in the x direction at point (i, j), Represents the forward and backward difference operators in the z direction at point (i, j) for travel time t.

In order to ensure the calculation accuracy, the second-order upwind difference operator is used. which is:

where Δx and Δz are the grid sizes in the X and Z directions, respectively.

Finally, using the travel time values at the eight nodes (i,j±1), (i,j±2), (i±1,j), (i±2,j) around the target point (i,j), There are 6 combinations to calculate the travel time value of point (i, j), and the second-order difference calculation formula of the 6 combinations is as follows:

In the formula, t is the seismic wave travel time, s is the slowness of the seismic wave, and h=Δx=Δz is the square grid size.

Using the above formula to calculate multiple traveltimes, according to the basic principle of narrow-band expansion, select the minimum traveltime point in the obtained multi-valued solution, which is the target point traveltime. The calculation formula is as follows:

t _i,j = min(t ¹ _i,j ,t ² _i,j ,t ³ _i,j ,t ⁴ _i,j ,t ⁵ _i,j ,t ⁶ _i,j )

After traversing all nodes, the seismic wave time field can be obtained.

Preferably, the wave field extraction step executes the following sub-steps:

Next, extract the seismic wave field of each detection point, that is, find the time corresponding to the position of the detection point in the time field t(x, z) obtained in the previous step, which is the earthquake travel time from the original source point to the detection point, forward modeling Finish.

This embodiment greatly improves the efficiency of microseismic monitoring forward modeling. In microseismic monitoring, the number of seismic sources is much larger than the number of detection points (Ns >> Nr). For example, in the microseismic monitoring of horizontal well fracturing in coal fields, the number of detection points (Nr) is generally more than a dozen, while the number of microseismic events (each event corresponds to a source Points) (Ns) are often hundreds or even thousands, and the difference between the two is 1 to 2 orders of magnitude. When simulating this process forward, if the seismic wave field of each source point is obtained according to the existing method, assuming that the simulation time of a source point is T, the time required to simulate the whole process is Ns×T. According to the method described in this patent, after the source detection is exchanged, only the seismic wave field of each detection point is required. Since the velocity model remains unchanged, the simulation time of each detection point is also T, and the time required to simulate the entire process is Nr ×T; the time difference between the two schemes is (Ns-Nr)×T. It can be seen that the greater the difference between Ns and Nr, the greater the speedup ratio.

In the specific implementation, it is necessary to compare the difference between Ns and Nr. If Ns>Nr, calculate the seismic wave field of each detection point; if Ns<Nr, calculate the seismic wave field of each source point. The process is shown in Figure 1. Monitoring forward modeling is generally suitable for the former.

The comparative effect of the forward modeling method described in the present invention and the conventional forward modeling method. This example is only described in two dimensions. The specific physical problem is described as follows: each of the four vertices of a square model with a side length of 100m is placed with a detector point, a total of 4 detector points, and the model speed is a constant speed of 1000m/s. The grid size is 0.1m*0.1m, the grid size is 1001*1001, 10 sources are randomly arranged in the square, and the travel time of each source pair is calculated. The operation time of the method of the present invention is 3058ms, the operation time of the conventional method is 7898ms, and the speed-up ratio is 2.58 times.

This embodiment only briefly illustrates the improvement of the computing efficiency of the present invention in a two-dimensional situation. In practical applications, all three-dimensional models are used, and the number of sources is far greater than the number of detection points, and the acceleration ratio will be larger.

The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which the present invention pertains can make various modifications or additions to the described specific embodiments or substitute in similar manners, but will not deviate from the spirit of the present invention or go beyond the definitions of the appended claims range.

Claims (5)

1. A method of accelerating microseismic numerical simulation comprising:

a forward parameter setting step, namely establishing a grid speed model according to the geological condition of an object to be simulated and setting the coordinates of a wave detector;

a time field calculation step, namely determining a seismic wave field calculation mode according to the number Ns of seismic source points and the number Nr of wave detection points; if Ns is larger than Nr, performing source detection interchange and wave field calculation, namely taking each wave detection point as a seismic source point, and respectively calculating the point of each wave detection point as a seismic wave field of a seismic source;

and a wave field extraction step, namely extracting the seismic wave field of each wave detection point, and finishing forward modeling.

2. The method according to claim 1, wherein the time field calculating step specifically comprises:

a seismic wave field determining sub-step, namely calculating a path function equation of the seismic wave field of the current seismic source point under a known velocity model;

a substep of equation calculation of an equation, namely discretizing the equation of the equation by adopting an upwind difference method;

the earthquake travel time operator step, which is to calculate a second-order difference calculation formula of the earthquake travel time by using travel time values at nodes in eight directions around a target point to obtain a multi-valued solution of the travel time of the target point;

and a target travel time determining sub-step, namely selecting the minimum travel time point from the multi-value solution of the calculated target travel time according to the narrow-band extension basic principle, namely the target travel time.

3. A method of accelerating a microseismic numerical simulation of the type defined in claim 2 wherein in the seismic wavefield determination substep the equation of the path function of the seismic wavefield is determined based on the equation:

where t is the seismic travel time, s is the slowness, i.e., the reciprocal of the velocity, and x, z represent 2 coordinate axes of a two-dimensional rectangular coordinate system.

4. The method according to claim 2, wherein in the earthquake travel time calculation operator step, the above formula is discretized by an upwind difference method, specifically:

represents the forward and backward difference operators of the travel time t in the x direction at the point (i, j),

a difference operator representing the travel time t forward and backward in the z direction at point (i, j);

wherein:

in the formula, Δ x and Δ z are grid sizes in the X, Z direction, respectively.

5. The method according to claim 2, wherein the seismic travel time operator step calculates a multi-valued solution of the travel time of the target point using travel time values at eight nodes around the target point based on the following formula:

where t is the seismic travel time, t_i,jIs the seismic travel time at point (i, j), s is the slowness of the seismic, h ═ Δ x ═ Δ z is the square grid size,

calculating a plurality of travel times by using the formula, and selecting the minimum travel time point from the solved multi-value solution according to a narrow-band extension basic principle, wherein the minimum travel time point is the travel time of a target point, and the calculation formula is as follows:

t_i,j＝min(t¹ _i,j,t² _i,j,t³ _i,j,t⁴ _i,j,t⁵ _i,j,t⁶ _i,j)。

Images