CN116660980A - A Microseismic Location Method Based on Improved Equation - Google Patents

Abstract

The invention discloses a microseism positioning method based on an improved equation of journey, which comprises the following steps: s1, dividing a delineated monitoring area into two-dimensional discrete grid points, and arranging detectors on the two-dimensional discrete grid points; s2, sequentially taking each detector as a new seismic source, and calculating the travel time of grid points near the new seismic source by adopting a spherical wave approximation algorithm; s3, calculating the travel time of the rest grid points by adopting a multi-template rapid travel algorithm according to the travel time of the grid points near the new seismic source obtained in the step S2; s4, calculating function values of travel time of all grid points by using a time residual function, and positioning the position of the seismic source according to the calculated function values; according to the invention, the spherical wave approximation algorithm is adopted to calculate the travel time of the grid points near the new seismic source, and the multi-template rapid travel algorithm is adopted to calculate the travel time of the rest grid points, so that the error of the travel time field in the diagonal direction and near the new seismic source point is reduced, the accuracy of the travel time field is improved, and the position of the microseismic seismic source can be more accurately positioned.

Description

A Microseismic Location Method Based on Improved Equation

technical field

The invention relates to the field of microseismic positioning, in particular to a microseismic positioning method based on an improved eigenfunction equation.

Background technique

For the exploration and development of unconventional oil and gas reservoirs, microseismic monitoring technology is an effective way to quantitatively analyze the development of fractures in the process of hydraulic fracturing reservoir stimulation.

Inversion location of microseismic events is the core issue in microseismic monitoring. According to different inversion theories, microseismic inversion problems can be roughly divided into two types: traveltime inversion and energy inversion. The earliest microseismic inversion positioning method borrowed from the positioning methods in natural earthquakes, such as Geiger positioning method and triangulation positioning method. This type of method belongs to the travel time inversion method, and its main principle is to use the difference between the theoretical travel time of the microseismic signal from the source to the receiver and the actual observed travel time as the objective function according to the space and time relationship of the propagation of the microseismic event in the formation. And solve the source location and the time of earthquake. Subsequently, microseismic source location methods based on energy superposition and migration-like methods appeared, such as grid search method, reverse time migration method, etc. This type of method avoids the first-arrival picking of microseismic events and finds the maximum energy after superposition point as the microseismic source.

The ray tracing method used in the traditional travel time inversion method has problems such as caustics and shadow areas when the geological model is strongly non-uniform, making it difficult to determine the ray path. The ray paths in inhomogeneous media can be obtained accurately by numerically solving the equations. At present, the method commonly used to solve the equation of the equation is the fast marching method, but the travel time field calculated by this method has a large error in the diagonal direction and near the source point.

Contents of the invention

Aiming at the above-mentioned deficiencies in the prior art, the present invention provides a microseismic positioning method based on the improved eigenfunction equation, which solves the shortcomings of the traditional ray tracing method in the case of strong and non-uniform geological conditions, and solves the problem with the fast marching algorithm to solve the problem. Compared with the function equation, the invention improves the accuracy of the travel time field, thereby realizing the precise positioning of microseismic signals.

In order to achieve the above-mentioned purpose of the invention, the technical scheme adopted in the present invention is:

A microseismic positioning method based on the improved eigenfunction equation, comprising the following steps:

S1. Divide the delineated monitoring area into two-dimensional discrete grids, and arrange geophones on the two-dimensional discrete grid points;

S2. Taking each geophone as a new seismic source in turn, using the spherical wave approximation algorithm to calculate the travel time of the grid points near the new seismic source;

S3. According to the travel time of the grid points near the new seismic source obtained in step S2, the travel time of the remaining grid points is calculated by using a multi-template fast-travel algorithm;

S4. Using the time residual function to calculate the travel time function values of all grid points, and locate the earthquake source according to the calculated function values.

Further, step S1 specifically includes:

According to the microseismic monitoring task, delineate the area to be monitored, establish a two-dimensional rectangular coordinate system, take the distance from the microseismic event plane projection to the origin as the abscissa, and take the depth of the microseismic event as the ordinate, divide the delineated monitoring area into Two-dimensional discrete grids, the size of each grid is d, and geophones are arranged in a vertical shaft on the two-dimensional discrete grid points.

Further, step S2 specifically includes:

S21. Using the geophone (x _k , z _k ) as a new seismic source;

_S22 . _{Insert a} grid point t ₀ ( _x , z ₀ ), use the Lagrange interpolation method to calculate the travel time of the new source to the grid point t ₀ (x, z ₀ ), and calculate the arrival time of the new source through the grid point t ₀ (x, z ₀ ) to the grid point t(x+ Δx, z ₂ ) travel time;

S23. Minimize the travel time of the grid point t(x+Δx, z ₂ ), and calculate the effect of the travel time of the grid point t(x+Δx, z ₂ ) on the arrival of the new source at the grid point t ₀ (x, z ₀ ) Partial derivative at depth z ₀ is close to 0, the depth z ₀ at which the new source reaches the grid point t ₀ (x,z ₀ ) is obtained, and then when the travel time of the grid point t(x+Δx,z ₂ ) is the smallest, the grid point The travel time of t ₀ (x,z ₀ ) and grid point t(x+Δx,z ₂ ).

Further, the calculation formula of the travel time from the new source to the grid point t ₀ (x, z ₀ ) in step S22 is:

Among them, t ₀ represents the travel time of the grid point t ₀ (x, z ₀ ), t ₁ and t ₂ represent the travel time of the new source to the two grid points respectively, and z ₀ represents the travel time from the new source to the grid point t ₀ ( x, z ₀ ), z ₁ represents the depth from the new source to grid point t ₁ (x, z ₁ ), and z ₂ represents the depth from the new source to grid point t ₂ (x, z ₂ ).

Further, the formula for calculating the travel time of the new seismic source to the grid point t(x+Δx,z ₂ ) in step S22 is:

Among them, t represents the travel time of the grid point near the new source, S represents the slowness of the cell, t ₀ represents the travel time of the grid point t ₀ (x, z ₀ ), z ₀ represents the travel time from the new source to the grid point t ₀ (x , z ₀ ), z ₂ represents the depth from the new source to the grid point t ₂ (x, z ₂ ), and Δx represents the discrete grid spacing in the direction of the x coordinate axis.

Further, the calculation formula of the travel time of the minimized point t(x+Δx,z ₂ ) in step S23 is:

in, Indicates the partial derivative of the travel time of the grid point t(x+Δx,z ₂ ) to the depth z ₀ at which the new source arrives at the grid point t ₀ (x,z ₀ ), W represents the average slowness, and S represents the slowness of the cell degree, Δx represents the discrete grid spacing in the direction of the x-coordinate axis, t ₀ represents the travel time of grid point t ₀ (x, z ₀ ), z ₀ represents the travel time from the new source to grid point t ₀ (x, z ₀ ) Depth, z ₁ represents the depth from the new source to grid point t ₁ (x, z ₁ ), and z ₂ represents the depth from the new source to point t ₂ (x, z ₂ ).

Further, step S3 specifically includes:

S31. Mark the grid points near the new seismic source calculated in step S2 as near points, and mark the remaining grid points as far points;

S32. The multi-template rapid marching algorithm uses template 1 and template 2 to calculate respectively, and uses template 1 to calculate the grid points on the horizontal axis or vertical axis of the two-dimensional Cartesian coordinate system with a distance of one grid size from the near point using the rapid marching algorithm. The travel time of template 2 is calculated by using the fast marching algorithm on the diagonal of the two-dimensional Cartesian coordinate system, and the distance from the near point is The travel time of the grid points of the grid size is calculated according to the first-order difference formula and the second-order difference formula, and the calculated minimum value is used as the travel time of the grid points, and these points are marked as narrow-band points;

S33. Find the grid point with the minimum travel time among all the narrow-band points, and mark it as the near point;

S34, finding the far point adjacent to the near point in step S33 is marked as a narrowband point;

S35, using the first-order difference approximation and the second-order difference approximation formula obtained by the two templates to update the travel time of the narrow-band point adjacent to the newly added approach point found in step S34;

S36. Repeat step S33 to step S35 until the narrowband point is empty.

Further, the first-order difference approximation formula of template 1 is as follows:

Among them, i represents the number of discrete grid points along the direction of the x-coordinate axis, j represents the number of discrete grid points along the direction of the y-coordinate axis, Δx represents the discrete grid spacing along the direction of the x-coordinate axis, and Δy represents the y-coordinate axis The discrete grid spacing in the direction, T _{i, j} represents the arrival time of the grid point (i, j), T ₁ and T ₂ respectively represent the x-coordinate axis direction and the y-coordinate axis at the grid point (i, j) The minimum time forward and backward in the direction, T ₁ =min(T _i-1,j ,T _i+1,j ), T ₂ =min(T _i,j-1 ,T _i,j+1 ), S _i,j represents the slowness of grid point (i,j).

Further, the approximate formula of the second-order difference of template 1 is as follows:

Among them, T _i,j represents the arrival time of the grid point (i,j), Δx represents the discrete grid spacing in the direction of the x-coordinate axis, Δy represents the discrete grid spacing in the direction of the y-coordinate axis, S _i,j Indicates the slowness of the grid point (i, j), T ₁ and T ₂ respectively indicate the minimum time forward and backward in the direction of the x-coordinate axis and the y-coordinate axis at the grid point (i, j),

Further, the first-order difference approximation formula of template 2 is as follows:

Among them, T _{i, j} represents the arrival time of the grid point (i, j), T _v represents the minimum time forward and backward in the direction of the two diagonal lines at the grid point (i, j), S _{i, j} represents the slowness of the grid point (i, j), and d represents the grid size.

Further, the approximate formula of the second-order difference of template 2 is as follows:

Among them, T _{i, j} represents the arrival time of the grid point (i, j), T _v represents the minimum time forward and backward in the direction of the two diagonal lines at the grid point (i, j), S _{i, j} represents the slowness of the grid point (i, j), and d represents the grid size.

Further, step S4 specifically includes:

Find the minimum function value calculated by the time residual function, and the grid point corresponding to the found minimum function value is the location of the earthquake source.

Further, the time residual function in step S4 is specifically:

Among them, f _Ns (x, z) represents the time residual function value of the point (x, z), k represents the number of geophones, Ns represents the number of seismic events, T _i (x, z) represents the calculated value of the i-th geophone t _i (j) represents the arrival time of the jth microseismic event detected by the ith geophone, and τ ₀ (j) represents the occurrence of the jth microseismic event Shock time.

The present invention has the following beneficial effects:

The invention adopts spherical wave approximation algorithm to calculate the travel time of grid points near the new seismic source, and adopts multi-template fast marching algorithm to calculate the travel time of other grid points, which reduces the error of the travel time field in the diagonal direction and near the new seismic source point, and improves the The accuracy of the travel time field can more accurately locate the location of the microseismic source.

Description of drawings

Fig. 1 is a kind of schematic flow sheet of the microseismic positioning method based on the improved eigen equation of the present invention;

Fig. 2 is the schematic diagram of spherical wave approximation algorithm;

Fig. 3 is the schematic diagram of template 1 and template 2 used by multi-template fast marching algorithm;

Figure 4 is a schematic diagram of the error of the travel time field at the source position (51,51) under the T1 model during the test phase;

Fig. 5 is a schematic diagram of the error of the travel time field at the source position of (1,1) under the T1 model during the test phase;

Figure 6 is a schematic diagram of the travel time field error at the source position (51,51) under the T2 model during the test phase;

Fig. 7 is a schematic diagram of the error of the travel time field at the source position of (51,51) under the T3 model during the test stage;

Figure 8 is a schematic diagram of the layered velocity model in the testing phase;

Figure 9 is a schematic diagram of the two-dimensional layered observation system and the comparison of positioning effects in the testing stage.

Detailed ways

The specific embodiments of the present invention are described below so that those skilled in the art can understand the present invention, but it should be clear that the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art, as long as various changes Within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are included in the protection list.

As shown in Figure 1, a microseismic location method based on the improved equating equation includes the following steps S1-S4:

S1. Divide the delineated monitoring area into two-dimensional discrete grids, and arrange geophones on the two-dimensional discrete grid points.

In an optional embodiment of the present invention, in this embodiment, a ray propagation calculation method based on grid expansion is used to calculate the ray first arrival time.

Specifically, step S1 includes: delimiting the area to be monitored according to the microseismic monitoring task, establishing a two-dimensional rectangular coordinate system, taking the distance from the microseismic event plane projection to the origin as the abscissa, and taking the depth of the microseismic event as the ordinate, The delineated monitoring area is divided into two-dimensional discrete grids, the size of each grid is d, and the geophones are arranged in the form of shafts on the two-dimensional discrete grid points, and the number of geophones is k.

S2. Taking each geophone as a new seismic source in turn, and calculating the travel time of the grid points near the new seismic source using the spherical wave approximation algorithm.

In an optional embodiment of the present invention, in order to solve the disadvantage that the travel time calculated by the fast-moving algorithm has large errors near the source point, in this embodiment, according to the reciprocity principle of wave propagation, in the same model, the direct wave propagates from the source point The travel time to the geophone is equal to the travel time of the wave from the geophone to the source, and the geophone (x _k , z _k ) is used as the new source, and the spherical wave approximation algorithm is used to calculate the vicinity of the new source point (x _k -h: x _k +h:z _k -h:z _k +h) the travel time t of the grid point. And the seismic wave is calculated in the approximate form of spherical waves, which is more in line with the propagation characteristics of seismic waves. This method is suitable for any complex isotropic medium, and has high accuracy near the source point, but the calculation efficiency is low, so this embodiment only uses this method Calculate the travel time of some grid points near the new source point.

As shown in Figure 2, step S2 includes the following steps S21-S23:

S21. Use the geophone (x _k , z _k ) as a new seismic source.

_S22 . _{Insert a} grid point t ₀ ( _x , z ₀ ), use the Lagrange interpolation method to calculate the travel time of the new source to the grid point t ₀ (x, z ₀ ), and calculate the arrival time of the new source through the grid point t ₀ (x, z ₀ ) to the grid point t(x+ Δx, z ₂ ) travel time;

In an optional embodiment of the present invention, the spherical wave approximation algorithm in this embodiment proposes that the travel time of a grid point can be calculated from the travel time of the known grid points around, so this embodiment uses a new seismic source to reach the grid The travel time of point t ₁ (x, z ₁ ) and grid point t ₂ (x, z ₂ ) is used to calculate the travel time of the new source to grid point t ₀ (x, z ₀ ) and the travel time of the new source through grid point t ₀ ( x,z ₀ ) travel time to grid point t(x+Δx,z ₂ ).

Specifically, the calculation formula of the travel time from the new source to the grid point t ₀ (x, z ₀ ) in step S22 is:

Among them, t ₀ represents the travel time of the grid point t ₀ (x, z ₀ ), t ₁ and t ₂ represent the travel time of the new source to the two grid points respectively, and z ₀ represents the travel time from the new source to the grid point t ₀ ( x, z ₀ ), z ₁ represents the depth from the new source to grid point t ₁ (x, z ₁ ), and z ₂ represents the depth from the new source to grid point t ₂ (x, z ₂ ).

Specifically, the formula for calculating the travel time of the new seismic source to the grid point t(x+Δx,z ₂ ) in step S22 is:

Among them, t represents the travel time of the grid point near the new source, S represents the slowness of the cell, t ₀ represents the travel time of the grid point t ₀ (x, z ₀ ), z ₀ represents the travel time from the new source to the grid point t ₀ (x , z ₀ ), z ₂ represents the depth from the new source to the grid point t ₂ (x, z ₂ ), and Δx represents the discrete grid spacing in the direction of the x coordinate axis.

S23. Minimize the travel time of the grid point t(x+Δx, z ₂ ), and calculate the effect of the travel time of the grid point t(x+Δx, z ₂ ) on the arrival of the new source at the grid point t ₀ (x, z ₀ ) Partial derivative at depth z ₀ is close to 0, the depth z ₀ at which the new source reaches the grid point t ₀ (x,z ₀ ) is obtained, and then when the travel time of the grid point t(x+Δx,z ₂ ) is the smallest, the grid point The travel time of t ₀ (x,z ₀ ) and grid point t(x+Δx,z ₂ ).

Specifically, the calculation formula of the travel time of the minimized point t(x+Δx,z ₂ ) in step S23 is:

in, Indicates the partial derivative of the travel time of the grid point t(x+Δx,z ₂ ) to the depth z ₀ at which the new source arrives at the grid point t ₀ (x,z ₀ ), and W represents the average slowness, which is only related to the current surrounding coordinates It is related to the travel time, separated from the connection of the virtual source, S represents the slowness of the cell, Δx represents the discrete grid spacing in the direction of the x-coordinate axis, t ₀ represents the travel time of the grid point t ₀ (x,z ₀ ), z ₀ Denotes the depth from the new source to grid point t ₀ (x, z ₀ ), z ₁ represents the depth from the new source to grid point t ₁ (x, z ₁ ), z ₂ represents the depth from the new source to point t ₂ (x, z ₂ ) depth.

S3. According to the travel time of the grid points near the new seismic source obtained in step S2, the travel time of the remaining grid points is calculated using a multi-template fast-travel algorithm.

In an optional embodiment of the present invention, this embodiment uses the travel time t of the grid points near the new seismic source obtained in step S2 as the initial area, and brings it into the multi-template fast-travel algorithm to calculate the travel time of the remaining grid points to obtain the entire monitoring The travel time of grid points in the region, and the traditional multi-template fast-travel algorithm calculates the travel time of other points from the source point, which will lead to large errors near the source, while the multi-template fast-travel algorithm proposed in this embodiment On the basis of the fast marching algorithm, the travel time of the remaining grid points is calculated by combining multiple templates and directional derivatives. This calculation can not only make up for the low accuracy of the multi-template fast marching algorithm near the source, but also improve the calculation efficiency.

As shown in Figure 3, step S3 includes the following steps S31-S36:

S31. Mark the grid points near the new seismic source calculated in step S2 as near points, and mark the rest of the grid points as far points.

S32. The multi-template rapid marching algorithm uses template 1 and template 2 to calculate respectively, and uses template 1 to calculate the grid points on the horizontal axis or vertical axis of the two-dimensional Cartesian coordinate system with a distance of one grid size from the near point using the rapid marching algorithm. The travel time of template 2 is calculated by using the fast marching algorithm on the diagonal of the two-dimensional Cartesian coordinate system, and the distance from the near point is The travel time of the grid points of the grid size is calculated according to the first-order difference formula and the second-order difference formula, and the calculated minimum value is used as the travel time of the grid points, and these points are marked as narrow-band points;

In an optional embodiment of the present invention, since both the fast marching algorithm and the second-order fast scanning algorithm ignore the information of the grid points on the diagonal, it is easy to generate a large numerical error in the direction of the diagonal, so, In this embodiment, template 1 is used to find the nearest neighbor point covered, and template 2 is used to find the neighbor point covered in the diagonal direction.

Specifically, the first-order difference approximation formula of template 1 is as follows:

Among them, i represents the number of discrete grid points along the direction of the x-coordinate axis, j represents the number of discrete grid points along the direction of the y-coordinate axis, Δx represents the discrete grid spacing along the direction of the x-coordinate axis, and Δy represents the y-coordinate axis The discrete grid spacing in the direction, T _{i, j} represents the arrival time of the grid point (i, j), T ₁ and T ₂ respectively represent the x-coordinate axis direction and the y-coordinate axis at the grid point (i, j) The minimum time forward and backward in the direction, T ₁ =min(T _i-1,j ,T _i+1,j ), T ₂ =min(T _i,j-1 ,T _i,j+1 ), S _i,j represents the slowness of grid point (i,j).

Specifically, the approximate formula of the second-order difference of template 1 is as follows:

Among them, T _i,j represents the arrival time of the grid point (i,j), Δx represents the discrete grid spacing in the direction of the x-coordinate axis, Δy represents the discrete grid spacing in the direction of the y-coordinate axis, S _i,j Indicates the slowness of the grid point (i, j), T ₁ and T ₂ respectively indicate the minimum time forward and backward in the direction of the x-coordinate axis and the y-coordinate axis at the grid point (i, j),

Specifically, the first-order difference approximation formula of template 2 is as follows:

Among them, T _{i, j} represents the arrival time of the grid point (i, j), T _v represents the minimum time forward and backward in the direction of the two diagonal lines at the grid point (i, j), S _{i, j} represents the slowness of the grid point (i, j), and d represents the grid size.

Specifically, the approximate formula of the second-order difference of template 2 is as follows:

Among them, T _{i, j} represents the arrival time of the grid point (i, j), T _v represents the minimum time forward and backward in the direction of the two diagonal lines at the grid point (i, j), S _{i, j} represents the slowness of the grid point (i, j), and d represents the grid size.

S33. Find the grid point with the minimum travel time among all the narrow-band points, and mark it as the near point.

S34. Search for the far points adjacent to the near point in step S33 to be marked as narrow-band points.

S35 , using the first-order difference approximation and the second-order difference approximation formula obtained from the two templates to update the travel time of the narrow-band point adjacent to the new approach point found in step S34 .

S36. Repeat step S33 to step S35 until the narrowband point is empty.

S4. Using the time residual function to calculate the travel time function values of all grid points, and locate the earthquake source according to the calculated function values.

In an optional embodiment of the present invention, in this embodiment, the time of the grid point (x, z) in the travel time field calculated according to all geophones, and the arrival time corresponding to the event detected by the geophones and the time of the earthquake, use the time residual function to calculate the travel time function values of these grid points, and find the minimum function value calculated by the time residual function, and the grid point corresponding to the minimum function value found is the location of the earthquake source.

Specifically, the time residual function in step S4 is specifically:

Among them, f _Ns (x, z) represents the time residual function value of the point (x, z), k represents the number of geophones, Ns represents the number of seismic events, T _i (x, z) represents the calculated value of the i-th geophone t _i (j) represents the arrival time of the jth microseismic event detected by the ith geophone, and τ ₀ (j) represents the occurrence of the jth microseismic event Shock time.

Specifically, in this embodiment, three speed models are used to verify the accuracy of the spherical wave approximation algorithm combined with the multi-template fast marching algorithm MSFM_SW in the present invention, and compared with the fast marching algorithm FMM and the multi-template fast marching algorithm MSFM, three The real traveltime field of this velocity model is calculated by the formula given in Table 1, and its velocity model can be calculated by the formula Obtained, Table 1 is as follows:

Table 1 Calculation formulas of the real traveltime fields of the three velocity models

Among them, x ₀ and y ₀ represent the source coordinates, x and y represent any grid point in the grid area, T ₁ (x) represents the real traveltime field under the T1 model, and t ₂ (x) represents the travel time field under the T2 model The real travel time field, T ₃ (x) represents the real travel time field under the T3 model.

In this embodiment, in order to measure the error of the solution of the calculated travel time field of the three speed models and the real travel time field T _a (x), the following error specification formula is used for specification:

l _∞ ＝mAx(|TT _a |)

Among them, n represents the total number of grid points, T represents the calculated travel time field, T _a represents the real travel time field, and L ₁ , L ₂ , L _∞ represent three different error specifications.

Table 2 and Table 3 in this embodiment compare the spherical wave approximation algorithm combined with the present invention and the multi-template fast marching algorithm MSFM_SW, the fast marching algorithm FMM and the multi-template fast marching algorithm MSFM. The errors at different times, combined with the three velocity models at T1, T2, and T3 given in Figure 4, Figure 5, Figure 6, and Figure 7, when the source is at (51,51) and (1,1), the above three The travel time error schematic diagram of this algorithm, from which it can be found that the travel time error calculated by the spherical wave approximation algorithm combined with the multi-template fast marching algorithm MSFM_SW is smaller than the travel time error of the fast marching algorithm FMM and the multi-template fast marching algorithm MSFM, and the accurate Travel time is the premise of accurate positioning, so the algorithm with smaller travel time error can make positioning more accurate.

Table 2 is as follows in the present embodiment:

Table 2 Error specifications for velocity model T1 calculated from different source points on a grid of size 101×101

Table 3 is as follows in the present embodiment:

Table 3. Error specifications of velocity models T2, T3 calculated from the grid center of size 101×101

Among them, Fig. 4 shows a schematic diagram of the traveltime field error with the focal position (51, 51) under the T1 model, and (a) in Fig. 4 shows a schematic diagram of the traveltime error of the spherical wave approximation algorithm combined with the multi-template fast marching algorithm MSFM_SW in the present invention , (b) in Figure 4 shows the schematic diagram of the travel time error of the multi-template fast-marching algorithm MSFM, and (c) in Figure 4 shows the schematic diagram of the travel-time error of the fast-marching algorithm FMM; Figure 5 shows that under the T1 model, the source position is (1,1 ) travel time field error schematic diagram, (a) in Fig. 5 shows the travel time error schematic diagram of the spherical wave approximation algorithm combined with the multi-template fast marching algorithm MSFM_SW in the present invention, and (b) in Fig. 5 shows the travel time error of the multi-template fast marching algorithm MSFM Schematic diagram, (c) in Fig. 5 represents the travel time error schematic diagram of fast marching algorithm FMM; Fig. 6 represents under the T2 model, and hypocenter position is the travel time field error schematic diagram of (51,51), in Fig. 6 (a) represents the present invention combines Schematic diagram of the traveltime error of the spherical wave approximation algorithm and the multi-template fast-marching algorithm MSFM_SW, Fig. 6 (b) shows the schematic diagram of the traveltime error of the multi-template fast-marching algorithm MSFM, and Fig. 6 (c) shows the schematic diagram of the traveltime error of the fast-marching algorithm FMM Fig. 7 shows under the T3 model, the hypocenter position is the travel time field error schematic diagram of (51,51), among Fig. 7 (a) represents the travel time error schematic diagram of the spherical wave approximation algorithm combined with the multi-template fast marching algorithm MSFM_SW of the present invention, (b) in Figure 7 shows a schematic diagram of the travel time error of the multi-template fast marching algorithm MSFM, and (c) in Figure 7 shows a schematic diagram of the travel time error of the fast marching algorithm FMM.

Specifically, in order to verify the travel time accuracy of the layered model of the present invention, the layered velocity model shown in Figure 8 is established in this embodiment, the inverted triangle in Figure 8 indicates the source position, and the black hollow circle indicates the travel time calculation point, and set The source location is (400,321), and 117 grid points are set in the observation system to detect the travel time, and the travel time calculated by the traditional ray tracing method is used as the real travel time for comparison, and the travel time errors of the three methods are obtained. Table 4, Table 4 is as follows Show:

Table 4 The travel time error of 117 grid points from the layered model with a size of 800×600

In this embodiment, by analyzing Table 4, it can be concluded that the accuracy of the method used in the present invention is still high under the layered model, which can meet the positioning requirements, and the spherical wave approximation algorithm combined with the present invention and the multi-template fast marching algorithm MSFM_SW are obtained. The travel time error is smaller than the travel time error obtained by using the fast marching algorithm FMM and the multi-template fast marching algorithm MSFM.

Specifically, in order to verify the positioning accuracy of the present invention, a two-dimensional layered observation system as shown in Figure 9(a) is established in this embodiment, the black inverted triangle in Figure 9(a) indicates the position of the station, and the black circle Indicate the real position of the seismic event, and use 12 stations to locate 12 events, use the travel time calculated by the traditional ray tracing method as the real travel time, adopt the spherical wave approximation algorithm combined with the present invention and the multi-template fast marching algorithm MSFM_SW and the fast marching Algorithm FMM performs positioning, and the positioning effect is shown in Figure 9 (b). Figure 9 (b) shows a comparison of the positioning effects of the two algorithms. The left triangle represents the positioning position of the fast-moving algorithm FMM, and the black cross represents the spherical surface combined by the present invention Wave approximation algorithm and multi-template fast marching algorithm MSFM_SW positioning position, from which it can be obtained that the error of adopting the spherical wave approximation algorithm combined by the present invention and multi-template fast marching algorithm MSFM_SW positioning will be smaller than the error of positioning using fast marching algorithm FMM. close to the real position.

In the present invention, specific examples have been applied to explain the principles and implementation methods of the present invention, and the descriptions of the above examples are only used to help understand the method of the present invention and its core idea; meanwhile, for those of ordinary skill in the art, according to this The idea of the invention will have changes in the specific implementation and scope of application. To sum up, the contents of this specification should not be construed as limiting the present invention.

Those skilled in the art will appreciate that the embodiments described here are to help readers understand the principles of the present invention, and it should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical revelations disclosed in the present invention without departing from the essence of the present invention, and these modifications and combinations are still within the protection scope of the present invention.

Claims (9)

1. The microseism positioning method based on the improved equation of the journey is characterized by comprising the following steps of:

s1, dividing a delineated monitoring area into two-dimensional discrete grid points, and arranging detectors on the two-dimensional discrete grid points;

s2, sequentially taking each detector as a new seismic source, and calculating the travel time of grid points near the new seismic source by adopting a spherical wave approximation algorithm;

s3, calculating the travel time of the rest grid points by adopting a multi-template rapid travel algorithm according to the travel time of the grid points near the new seismic source obtained in the step S2;

s4, calculating function values of travel time of all grid points by using a time residual function, and positioning the seismic source position according to the calculated function values.

2. The method for positioning a microseism based on an improved equation of travel function according to claim 1, wherein step S1 specifically comprises:

according to a microseism monitoring task, a region to be monitored is defined, a two-dimensional rectangular coordinate system is established, the distance from a microseism event plane to an origin is taken as an abscissa, the depth of a microseism event is taken as an ordinate, the defined monitoring region is divided into two-dimensional discrete grids, the size of each grid is d, and detectors are distributed on the two-dimensional discrete grid points in a vertical shaft mode.

3. The method for positioning a microseism based on an improved equation of travel function according to claim 1, wherein step S2 specifically comprises:

s21, detector (x) _k ,z _k ) As a new source;

s22, reaching grid point t according to the new seismic source ₁ (x,z ₁ ) Andgrid point t ₂ (x,z ₂ ) Is inserted between two grid points ₀ (x,z ₀ ) Calculating new seismic source arrival grid point t by Lagrange interpolation method ₀ (x,z ₀ ) And calculates the travel time of the new source past grid point t ₀ (x,z ₀ ) Reach grid point t (x+Δx, z ₂ ) Is used for walking time;

s23, minimizing grid point t (x+Δx, z ₂ ) Calculates the travel time of grid point t (x+Δx, z) ₂ ) To the new source to reach grid point t ₀ (x,z ₀ ) Depth z of (2) ₀ Is a deviator of (a)A point approaching 0, obtaining the arrival of the new seismic source at the grid point t ₀ (x,z ₀ ) Depth z of (2) ₀ Further, the current grid point t (x+Δx, z) ₂ ) At grid point t ₀ (x,z ₀ ) And grid point t (x+Δx, z ₂ ) Is a running time size of (2).

4. A method of microseismic localization based on improved process equations as claimed in claim 3, wherein in step S22 the new source is brought to grid point t ₀ (x,z ₀ ) The calculation formula of the travel time of (2) is as follows:

wherein ,t₀ Representing grid point t ₀ (x,z ₀ ) Time of travel, t ₁ 、t ₂ Representing the travel time of the new source to reach the two grid points respectively, z ₀ Representing arrival at grid point t from new source ₀ (x,z ₀ ) Depth, z ₁ Representing arrival at grid point t from new source ₁ (x,z ₁ ) Depth, z ₂ Representing arrival at grid point t from new source ₂ (x,z ₂ ) Is a depth of (2);

in step S22, the new source arrives at grid point t (x+Δx, z ₂ ) The calculation formula of the travel time of (2) is as follows:

where t represents the travel time of the grid point near the new source, S represents the cell slowness, t ₀ Representing grid point t ₀ (x,z ₀ ) Is z ₀ Representing arrival at grid point t from new source ₀ (x,z ₀ ) Depth, z ₂ Representing arrival at grid point t from new source ₂ (x,z ₂ ) Δx represents the discrete grid spacing in the x coordinate axis direction;

the minimization point t (x+Δx, z in step S23 ₂ ) The calculation formula of the travel time of (2) is as follows:

wherein ,representing the grid points t (x+Δx, z ₂ ) To the new source to reach grid point t ₀ (x,z ₀ ) Depth z of (2) ₀ W represents average slowness, S represents cell slowness, deltax represents discrete grid spacing in the x coordinate axis direction, t ₀ Representing grid point t ₀ (x,z ₀ ) Is z ₀ Representing arrival at grid point t from new source ₀ (x,z ₀ ) Depth, z ₁ Representing arrival at grid point t from new source ₁ (x,z ₁ ) Depth, z ₂ Representing the arrival at point t from the new source ₂ (x,z ₂ ) Is a depth of (c).

5. The method for positioning a microseism based on an improved equation of travel function according to claim 1, wherein step S3 specifically comprises:

s31, marking grid points nearby the new seismic source calculated in the step S2 as near points and marking the rest grid points as far points;

s32, calculating by using a template 1 and a template 2 respectively according to a multi-template rapid travel algorithm, calculating travel time of grid points with a grid size from a near point on a transverse axis or a longitudinal axis of a two-dimensional rectangular coordinate system by using the template 1 according to the rapid travel algorithm, calculating travel time of grid points with a grid size from a near point on a diagonal of the two-dimensional rectangular coordinate system by using the template 2 according to the rapid travel algorithm, and calculating travel time of the grid points with the same grid size from the near pointCalculating the travel time of grid points with the grid size according to a first-order difference formula and a second-order difference formula, taking the calculated minimum value as the travel time of the grid points, and marking the points as narrow-band points;

s33, searching a grid point with the smallest travel time in all the narrowband points, and marking the grid point as a near point;

s34, searching a far point adjacent to the near point in the step S33 to be marked as a narrow-band point;

s35, updating the narrow-band point travel time adjacent to the new near point found in the step S34 by using a first-order difference approximation formula and a second-order difference approximation formula obtained by the two templates;

s36, cycling the steps S33 to S35 until the narrowband point is empty.

6. The method for microseismic localization based on improved equation of elevation of claim 5 wherein the first order differential approximation formula for template 1 is as follows:

wherein i denotes a discrete grid point number along the x-axis direction, j denotes a discrete grid point number along the y-axis direction, Δx denotes a discrete grid pitch along the x-axis direction, Δy denotes a discrete grid pitch along the y-axis direction, and T _i,j Representing the arrival time of grid point (i, j), T ₁ 、T ₂ Respectively representing the minimum time in the x-axis direction and the y-axis direction at grid point (i, j), T ₁ ＝min((T _i-1,j ,T _i+1,j )，T ₂ ＝min(T _i,j-1 ,T _i,j+1 )，S _i,j Representing the slowness of grid points (i, j);

the second order differential approximation formula for template 1 is as follows:

wherein ,T_i,j When the grid points (i, j) arrive, Δx represents the discrete grid pitch in the x-coordinate axis direction, Δy represents the discrete grid pitch in the y-coordinate axis direction, and S _i,j Representing the slowness of grid point (i, j), T ₁ 、T ₂ Representing the minimum time forward and backward in the x-axis direction and the y-axis direction at the grid point (i, j) respectively,

7. the method of claim 5, wherein the first order differential approximation formula for template 2 is as follows:

wherein ,T_i,j Representing the arrival time of grid point (i, j), T _v Represents the minimum time between the two diagonal directions at grid point (i, j), S _i,j Representing the slowness of grid points (i, j), d representing the grid size;

the second order differential approximation formula for template 2 is as follows:

wherein ,T_i,j representing the arrival time of grid point (i, j), T _v Represents the minimum time between the two diagonal directions at grid point (i, j), S _i,j The slowness of grid point (i, j) is represented, and d represents the grid size.

8. The method for positioning a microseism based on an improved equation of travel function according to claim 1, wherein step S4 specifically comprises:

searching a minimum function value obtained by calculating a time residual function, wherein a grid point corresponding to the found minimum function value is the position of the seismic source.

9. The method for positioning microseism based on improved function equation according to claim 1, wherein the time residual function in step S4 is specifically:

wherein ,f_Ns (x, z) represents the time residual function value of the point (x, z), k represents the number of detectors, ns represents the number of seismic events, T _i (x, z) represents the time of the midpoint (x, z) of the travel-time field calculated by the ith detector, t _i (j) Indicating the arrival time of the ith receiver at the detection of the jth microseismic event, τ ₀ (j) Representing the time of onset of the jth microseismic event.