CN116660980A - Microseism positioning method based on improved path function equation - Google Patents
Microseism positioning method based on improved path function equation Download PDFInfo
- Publication number
- CN116660980A CN116660980A CN202310730999.4A CN202310730999A CN116660980A CN 116660980 A CN116660980 A CN 116660980A CN 202310730999 A CN202310730999 A CN 202310730999A CN 116660980 A CN116660980 A CN 116660980A
- Authority
- CN
- China
- Prior art keywords
- grid
- point
- grid point
- representing
- travel
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 48
- 238000012544 monitoring process Methods 0.000 claims abstract description 13
- 238000004364 calculation method Methods 0.000 claims description 16
- 230000001351 cycling effect Effects 0.000 claims description 3
- 238000001514 detection method Methods 0.000 claims description 3
- 230000004807 localization Effects 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 25
- 238000012360 testing method Methods 0.000 description 6
- 230000007547 defect Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000000977 initiatory effect Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/288—Event detection in seismic signals, e.g. microseismics
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Emergency Management (AREA)
- Business, Economics & Management (AREA)
- Acoustics & Sound (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
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
Technical Field
The invention relates to the field of microseism positioning, in particular to a microseism positioning method based on an improved equation of journey.
Background
For the exploration and development of unconventional oil and gas reservoirs, the microseism monitoring technology is an effective method for quantitatively analyzing the crack development condition in the hydraulic fracturing oil reservoir reconstruction process.
Microseismic event inversion localization is a core problem in microseismic monitoring. Inversion problems of microseisms can be largely classified into two types, travel time inversion and energy inversion, according to different inversion theories. The earliest microseism inversion positioning method refers to the positioning method in natural earthquakes, such as Geiger positioning method, triangular positioning method and the like. The method belongs to a travel time inversion method, and the main principle is that according to the spatial and temporal relation of the propagation of microseism events in a stratum, the difference between the theoretical travel time of microseism signals from a seismic source to a detector and the actual observation travel time is used as an objective function to solve the position of the seismic source and the time of earthquake initiation. Then, a microseism focus positioning method based on an energy superposition and similar migration method, such as a grid search method, a reverse time migration method and the like, is adopted, the method avoids first arrival pickup of microseism events, and a superimposed energy maximum point is searched to serve as a microseism focus.
The ray tracing method used in the traditional travel time inversion method has the problems of scattered focus, shadow areas and the like when the geological model is strong and nonuniform, so that the ray path is difficult to determine. The ray paths in the non-uniform medium can be accurately obtained by a method of solving a program function equation through a numerical value. The current method for solving the equation of the path function is a fast travelling method, but the error of the travel time field calculated by the method is larger in the diagonal direction and near the source point.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a microseism positioning method based on an improved path function equation, solves the defect that a shadow area exists in the traditional ray tracing method under the condition of strong and nonuniform geology, and improves the accuracy of a travel time field compared with a fast travel algorithm solution path function equation, thereby realizing the accurate positioning of microseism signals.
In order to achieve the aim of the invention, the invention adopts the following technical scheme:
a microseism positioning method based on an improved equation of function 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 seismic source position according to the calculated function values.
Further, the step S1 specifically includes:
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.
Further, the step S2 specifically includes:
s21, detector (x) k ,z k ) As a new source;
s22, reaching grid point t according to the new seismic source 1 (x,z 1 ) And grid point t 2 (x,z 2 ) Is inserted between two grid points 0 (x,z 0 ) Calculating new seismic source arrival grid point t by Lagrange interpolation method 0 (x,z 0 ) And calculates the travel time of the new source past grid point t 0 (x,z 0 ) Reach grid point t (x+Δx, z 2 ) Is used for walking time;
s23, minimizing grid point t (x+Δx, z 2 ) Calculates the travel time of grid point t (x+Δx, z) 2 ) To the new source to reach grid point t 0 (x,z 0 ) Depth z of (2) 0 Is a deviator of (a)A point approaching 0, obtaining the arrival of the new seismic source at the grid point t 0 (x,z 0 ) Depth z of (2) 0 Further, the current grid point t (x+Δx, z) 2 ) At grid point t 0 (x,z 0 ) And grid point t (x+Δx, z 2 ) Is a running time size of (2).
Further, in step S22, the new source is moved to grid point t 0 (x,z 0 ) The calculation formula of the travel time of (2) is as follows:
wherein ,t0 Representing grid point t 0 (x,z 0 ) Time of travel, t 1 、t 2 Representing the travel time of the new source to reach the two grid points respectively, z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Depth, z 1 Representing arrival at grid point t from new source 1 (x,z 1 ) Depth, z 2 Representing arrival at grid point t from new source 2 (x,z 2 ) Is a depth of (c).
Further, in step S22, the new source reaches grid point t (x+Δx, z) 2 ) 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 0 Representing grid point t 0 (x,z 0 ) Is z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Depth, z 2 Representing arrival at grid point t from new source 2 (x,z 2 ) Δx represents the discrete grid spacing in the x coordinate axis direction.
Further, the minimum point t (x+Δx, z in step S23 2 ) The calculation formula of the travel time of (2) is as follows:
wherein ,representing the grid points t (x+Δx, z 2 ) To the new source to reach grid point t 0 (x,z 0 ) Depth z of (2) 0 W represents average slowness, S represents cell slowness, deltax represents discrete grid spacing in the x coordinate axis direction, t 0 Representing grid point t 0 (x,z 0 ) Is z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Depth, z 1 Representing arrival at grid point t from new source 1 (x,z 1 ) Depth, z 2 Representing the arrival at point t from the new source 2 (x,z 2 ) Is a depth of (c).
Further, the step S3 specifically includes:
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.
Further, the first order difference approximation formula of the 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 1 、T 2 Respectively representing the minimum time in the x-axis direction and the y-axis direction at grid point (i, j), T 1 ＝min(T i-1,j ,T i+1,j )，T 2 ＝min(T i,j-1 ,T i,j+1 )，S i,j Representing the slowness of grid points (i, j).
Further, the second order differential approximation formula of template 1 is as follows:
wherein ,Ti,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 1 、T 2 Representing the minimum time forward and backward in the x-axis direction and the y-axis direction at the grid point (i, j) respectively,
further, the first order difference approximation formula of the template 2 is as follows:
wherein ,Ti,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.
Further, the second order differential approximation formula of template 2 is as follows:
wherein ,Ti,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.
Further, the step S4 specifically includes:
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.
Further, in step S4, the time residual function is specifically:
wherein ,fNs (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, τ 0 (j) Representing the time of onset of the jth microseismic event.
The invention has the following beneficial effects:
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 source can be more accurately positioned.
Drawings
FIG. 1 is a schematic flow chart of a microseism locating method based on an improved equation of travel function according to the present invention;
FIG. 2 is a schematic diagram of a spherical wave approximation algorithm;
FIG. 3 is a schematic diagram of form 1 and form 2 used in the multi-form fast travel algorithm;
FIG. 4 is a schematic diagram of the time-to-field error for source location (51, 51) under the T1 model during the test phase;
FIG. 5 is a schematic diagram of the travel time field error for the source location (1, 1) under the T1 model during the test phase;
FIG. 6 is a schematic diagram of the time-to-field error for source location (51, 51) under the T2 model during the test phase;
FIG. 7 is a schematic diagram of the time-to-field error for source location (51, 51) under the T3 model during the test phase;
FIG. 8 is a schematic diagram of a layered velocity model during the test phase;
FIG. 9 is a schematic diagram showing a comparison of two-dimensional layered observation system and positioning effect during the test phase.
Detailed Description
The following description of the embodiments of the present invention is provided to facilitate understanding of the present invention by those skilled in the art, but it should be understood that the present invention is not limited to the scope of the embodiments, and all the inventions which make use of the inventive concept are protected by the spirit and scope of the present invention as defined and defined in the appended claims to those skilled in the art.
As shown in FIG. 1, a microseism locating method based on improved process function equations comprises the following steps S1-S4:
s1, dividing the delineated monitoring area into two-dimensional discrete grid points, and arranging detectors on the two-dimensional discrete grid points.
In an alternative embodiment of the present invention, a grid extension-based ray propagation calculation method is used in this embodiment to calculate the first arrival time of a ray.
Specifically, step S1 includes: 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, detectors are distributed on the two-dimensional discrete grid points in a vertical shaft mode, and the number of the distributed detectors is k.
S2, sequentially taking each detector as a new seismic source, and calculating the travel time of grid points nearby the new seismic source by adopting a spherical wave approximation algorithm.
In an alternative embodiment of the present invention, to solve the problem that the error in calculating the travel time in the vicinity of the source point is large in the fast travel algorithm, the travel time of the direct wave from the source to the detector is equal to the travel time of the wave from the detector to the source in the same model according to the reciprocity principle of wave propagation, the detector (x k ,z k ) As a new source, a spherical wave approximation algorithm is used to calculate the point (x k -h:x k +h:z k -h:z k +h) travel time t of grid points. The method is suitable for any complex isotropic medium, has high accuracy near the seismic source point, but has low calculation efficiency, so the embodiment only adopts the method to calculate the travel time of partial grid points near the new seismic source point.
As shown in fig. 2, step S2 includes the following steps S21 to S23:
s21, detector (x) k ,z k ) As a new source.
S22, reaching grid point t according to the new seismic source 1 (x,z 1 ) And grid point t 2 (x,z 2 ) Is inserted between two grid points 0 (x,z 0 ) Calculating new seismic source arrival grid by Lagrange interpolation methodPoint t 0 (x,z 0 ) And calculates the travel time of the new source past grid point t 0 (x,z 0 ) Reach grid point t (x+Δx, z 2 ) Is used for walking time;
in an alternative embodiment of the invention, the spherical wave approximation algorithm in this embodiment proposes a travel time for a grid point, which can be calculated from the surrounding known grid point travel, so this embodiment arrives at grid point t with the new source 1 (x,z 1 ) And grid point t 2 (x,z 2 ) Is calculated to reach grid point t 0 (x,z 0 ) Is passed through grid point t 0 (x,z 0 ) Reach grid point t (x+Δx, z 2 ) Is used for the travel time of the car.
Specifically, in step S22, the new source is set to grid point t 0 (x,z 0 ) The calculation formula of the travel time of (2) is as follows:
wherein ,t0 Representing grid point t 0 (x,z 0 ) Time of travel, t 1 、t 2 Representing the travel time of the new source to reach the two grid points respectively, z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Depth, z 1 Representing arrival at grid point t from new source 1 (x,z 1 ) Depth, z 2 Representing arrival at grid point t from new source 2 (x,z 2 ) Is a depth of (c).
Specifically, in step S22, the new source reaches grid point t (x+Δx, z 2 ) 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 0 Representing grid point t 0 (x,z 0 ) Is z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Is greater than the depth of (1)Degree, z 2 Representing arrival at grid point t from new source 2 (x,z 2 ) Δx represents the discrete grid spacing in the x coordinate axis direction.
S23, minimizing grid point t (x+Δx, z 2 ) Calculates the travel time of grid point t (x+Δx, z) 2 ) To the new source to reach grid point t 0 (x,z 0 ) Depth z of (2) 0 Is a deviator of (a)A point approaching 0, obtaining the arrival of the new seismic source at the grid point t 0 (x,z 0 ) Depth z of (2) 0 Further, the current grid point t (x+Δx, z) 2 ) At grid point t 0 (x,z 0 ) And grid point t (x+Δx, z 2 ) Is a running time size of (2).
Specifically, the minimum point t (x+Δx, z in step S23 2 ) The calculation formula of the travel time of (2) is as follows:
wherein ,representing the grid points t (x+Δx, z 2 ) To the new source to reach grid point t 0 (x,z 0 ) Depth z of (2) 0 W represents average slowness, is related to current surrounding coordinates and walk time only, is separated from the relation of virtual seismic sources, S represents cell slowness, deltax represents discrete grid spacing in the x coordinate axis direction, t 0 Representing grid point t 0 (x,z 0 ) Is z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Depth, z 1 Representing arrival at grid point t from new source 1 (x,z 1 ) Depth, z 2 Representing the arrival at point t from the new source 2 (x,z 2 ) Is a depth of (c).
And 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.
In an alternative embodiment of the present invention, the travel time t of the grid point near the new seismic source obtained in the step S2 is taken as an initial area, and is brought into the multi-template fast travel algorithm to calculate the travel times of the rest grid points, so as to obtain the grid point travel time of the whole monitoring area, while the conventional multi-template fast travel algorithm calculates the travel times of other points from one point of the seismic source, which results in larger error near the seismic source, and the multi-template fast travel algorithm provided in the present embodiment combines the multi-template and directional derivative to calculate the travel times of the rest grid points, so that the calculation can not only make up the defect that the precision of the multi-template fast travel algorithm is not high near the seismic source, but also can improve the calculation efficiency.
As shown in fig. 3, step S3 includes the following steps S31 to S36:
s31, marking the 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;
in an alternative embodiment of the present invention, since the fast traveling algorithm and the second order fast scanning algorithm both ignore the information of the grid points on the diagonal, a large numerical error is liable to occur in the diagonal direction, so that the present embodiment uses the template 1 to find the nearest neighbor covered and uses the template 2 to find the neighbor covered in the diagonal direction.
Specifically, the first-order differential approximation formula of the 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 1 、T 2 Respectively representing the minimum time in the x-axis direction and the y-axis direction at grid point (i, j), T 1 ＝min(T i-1,j ,T i+1,j )，T 2 ＝min(T i,j-1 ,T i,j+1 )，S i,j Representing the slowness of grid points (i, j).
Specifically, the second order differential approximation formula of the template 1 is as follows:
wherein ,Ti,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 1 、T 2 Representing the minimum time forward and backward in the x-axis direction and the y-axis direction at the grid point (i, j) respectively,
specifically, the first-order differential approximation formula of the template 2 is as follows:
wherein ,Ti,j Representing the arrival time of grid point (i, j), T v Expressed at grid points [ ]Minimum time between two diagonal directions of i, j), S i,j The slowness of grid point (i, j) is represented, and d represents the grid size.
Specifically, the second order differential approximation formula of template 2 is as follows:
wherein ,Ti,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.
And 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 the far point adjacent to the near point in the step S33, and marking the far point as a narrow-band point.
And S35, updating the narrow-band point travel time adjacent to the new nearest 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.
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.
In an alternative embodiment of the present invention, in this embodiment, according to the time of the grid points (x, z) in the travel time field calculated by all detectors, and according to the arrival time and the earthquake time corresponding to the event detected by the detectors, the function values of the travel times of the grid points are calculated by using a time residual function, the minimum function value calculated by the time residual function is found, and the grid point corresponding to the found minimum function value is the earthquake focus position.
Specifically, the time residual function in step S4 is specifically:
wherein ,fNs (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, τ 0 (j) Representing the time of onset of the jth microseismic event.
Specifically, in this embodiment, three speed models are adopted to verify the accuracy of the spherical wave approximation algorithm and the multi-template fast traveling algorithm MSFM_SW combined by the present invention, and compared with the fast traveling algorithm FMM and the multi-template fast traveling algorithm MSFM, the real travel time fields of the three speed models are calculated by the formulas given in Table 1, and the speed models can be calculated by the formulasThe results are shown in Table 1 below:
table 1 calculation formulas of three speed model real travel time fields
wherein ,x0 and y0 Representing the source coordinates, x and y representing any grid point within the grid area, T 1 (x) Representing the real travel time field under the T1 model, T 2 (x) Representing the real travel time field under the T2 model, T 3 (x) Representing the real travel time field under the T3 model.
In order to measure the travel time field and the real travel time field T of the three speed models obtained by calculation in the embodiment a (x) The error of the solution of (2) is normalized by the following error normalization formula:
l ∞ ＝mAx(|T-T a |)
where n represents the total number of grid points, T represents the calculated travel time field, T a Representing the real travel time field, L 1 、L 2 、L ∞ Three different error specifications are shown.
In this embodiment, table 2 and table 3 compare the spherical wave approximation algorithm and the multi-template fast traveling algorithm msfm_sw, the fast traveling algorithm FMM and the multi-template fast traveling algorithm MSFM, and the three algorithms are different in errors of the positions of the seismic sources under different speed models, and combine the three speed models of T1, T2 and T3 shown in fig. 4, 5, 6 and 7, and the travel time error diagrams of the three algorithms are shown when the seismic sources are at the (51, 51) and (1, 1) positions, so that the travel time error calculated by adopting the spherical wave approximation algorithm and the multi-template fast traveling algorithm msfm_sw combined by the invention is smaller than the travel time error calculated by adopting the fast traveling algorithm FMM and the multi-template fast traveling algorithm MSFM, and the accurate travel time is the premise of accurate positioning, so that the algorithm with smaller travel time error can make the positioning more accurate.
Table 2 in this example is shown below:
table 2 error specification of velocity model T1 calculated from different source points of a grid of size 101 x 101
Table 3 in this example is shown below:
table 3 error specification of velocity models T2, T3 calculated from grid center of size 101 x 101
Wherein, fig. 4 shows a schematic diagram of travel time field error of the source position (51, 51) under the T1 model, fig. 4 (a) shows a schematic diagram of travel time error of the spherical wave approximation algorithm and the multi-template fast travel algorithm msfm_sw combined in the present invention, fig. 4 (b) shows a schematic diagram of travel time error of the multi-template fast travel algorithm MSFM, and fig. 4 (c) shows a schematic diagram of travel time error of the fast travel algorithm FMM; FIG. 5 shows a schematic diagram of the travel time field error of the source position (1, 1) under the T1 model, wherein (a) in FIG. 5 shows a schematic diagram of the travel time error of the spherical wave approximation algorithm and the multi-template fast travel algorithm MSFM_SW combined by the present invention, (b) in FIG. 5 shows a schematic diagram of the travel time error of the multi-template fast travel algorithm MSFM, and (c) in FIG. 5 shows a schematic diagram of the travel time error of the fast travel algorithm FMM; FIG. 6 shows a schematic diagram of the travel time field error of the source position (51, 51) under the T2 model, wherein (a) in FIG. 6 shows a schematic diagram of the travel time error of the spherical wave approximation algorithm and the multi-template fast travel algorithm MSFM_SW combined by the present invention, (b) in FIG. 6 shows a schematic diagram of the travel time error of the multi-template fast travel algorithm MSFM, and (c) in FIG. 6 shows a schematic diagram of the travel time error of the fast travel algorithm FMM; fig. 7 shows a schematic diagram of a travel time field error of a seismic source position (51, 51) under a T3 model, fig. 7 (a) shows a schematic diagram of a travel time error of a spherical wave approximation algorithm and a multi-template fast travel algorithm msfm_sw combined in the present invention, fig. 7 (b) shows a schematic diagram of a travel time error of a multi-template fast travel algorithm MSFM, and fig. 7 (c) shows a schematic diagram of a travel time error of a fast travel algorithm FMM.
Specifically, in order to verify the travel time accuracy of the layered model of the present invention, in this embodiment, a layered velocity model as shown in fig. 8 is built, in fig. 8, an inverted triangle indicates the position of the seismic source, a black open circle indicates the travel time calculation point, and the position of the seismic source is set to be (400,321), 117 grid points are set in the observation system to detect the travel time, and the travel time calculated by the conventional ray tracing method is compared as the real travel time, so as to obtain a travel time error table 4 of three methods, where table 4 is as follows:
table 4 travel time error from 117 grid points under a 800 x 600 layered model
In the embodiment, the analysis table 4 can obtain that the accuracy of the method used by the invention is still high under the layered model, the positioning requirement can be met, and the travel time error obtained by adopting the spherical wave approximation algorithm and the multi-template fast travel algorithm MSFM_SW combined by the invention is smaller than the travel time error obtained by adopting the fast travel algorithm FMM and the multi-template fast travel algorithm MSFM.
Specifically, in order to verify the positioning accuracy of the present invention, a two-dimensional layered observation system as shown in fig. 9 (a) is established in this embodiment, in fig. 9 (a), a black inverted triangle indicates a station position, a black circle indicates a real position of a seismic event, 12 stations are utilized to position 12 events, the travel time calculated by using the conventional ray tracing method is used as the real travel time, the spherical wave approximation algorithm, the multi-template fast travel algorithm msfm_sw and the fast travel algorithm FMM combined by the present invention are adopted to position, the positioning effect is as shown in fig. 9 (b), fig. 9 (b) indicates a comparison graph of the positioning effects of the two algorithms, the left triangle indicates a fast travel algorithm FMM positioning position, and the black fork number indicates the spherical wave approximation algorithm and the multi-template fast travel algorithm msfm_sw positioning position combined by the present invention, so that the error of the spherical wave approximation algorithm and the multi-template fast travel algorithm msfm_sw combined by the present invention is smaller than the error of the positioning by adopting the fast travel algorithm FMM and is closer to the real position.
The principles and embodiments of the present invention have been described in detail with reference to specific examples, which are provided to facilitate understanding of the method and core ideas of the present invention; meanwhile, as those skilled in the art will have variations in the specific embodiments and application scope in accordance with the ideas of the present invention, the present description should not be construed as limiting the present invention in view of the above.
Those of ordinary skill in the art will recognize that the embodiments described herein are for the purpose of aiding the reader in understanding the principles of the present invention and should be understood that the scope of the invention is not limited to such specific statements and embodiments. Those of ordinary skill in the art can make various other specific modifications and combinations from the teachings of the present disclosure without departing from the spirit thereof, and such modifications and combinations remain within the scope of the present disclosure.
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 1 (x,z 1 ) Andgrid point t 2 (x,z 2 ) Is inserted between two grid points 0 (x,z 0 ) Calculating new seismic source arrival grid point t by Lagrange interpolation method 0 (x,z 0 ) And calculates the travel time of the new source past grid point t 0 (x,z 0 ) Reach grid point t (x+Δx, z 2 ) Is used for walking time;
s23, minimizing grid point t (x+Δx, z 2 ) Calculates the travel time of grid point t (x+Δx, z) 2 ) To the new source to reach grid point t 0 (x,z 0 ) Depth z of (2) 0 Is a deviator of (a)A point approaching 0, obtaining the arrival of the new seismic source at the grid point t 0 (x,z 0 ) Depth z of (2) 0 Further, the current grid point t (x+Δx, z) 2 ) At grid point t 0 (x,z 0 ) And grid point t (x+Δx, z 2 ) 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 0 (x,z 0 ) The calculation formula of the travel time of (2) is as follows:
wherein ,t0 Representing grid point t 0 (x,z 0 ) Time of travel, t 1 、t 2 Representing the travel time of the new source to reach the two grid points respectively, z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Depth, z 1 Representing arrival at grid point t from new source 1 (x,z 1 ) Depth, z 2 Representing arrival at grid point t from new source 2 (x,z 2 ) Is a depth of (2);
in step S22, the new source arrives at grid point t (x+Δx, z 2 ) 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 0 Representing grid point t 0 (x,z 0 ) Is z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Depth, z 2 Representing arrival at grid point t from new source 2 (x,z 2 ) Δx represents the discrete grid spacing in the x coordinate axis direction;
the minimization point t (x+Δx, z in step S23 2 ) The calculation formula of the travel time of (2) is as follows:
wherein ,representing the grid points t (x+Δx, z 2 ) To the new source to reach grid point t 0 (x,z 0 ) Depth z of (2) 0 W represents average slowness, S represents cell slowness, deltax represents discrete grid spacing in the x coordinate axis direction, t 0 Representing grid point t 0 (x,z 0 ) Is z 0 Representing arrival at grid point t from new source 0 (x,z 0 ) Depth, z 1 Representing arrival at grid point t from new source 1 (x,z 1 ) Depth, z 2 Representing the arrival at point t from the new source 2 (x,z 2 ) 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 1 、T 2 Respectively representing the minimum time in the x-axis direction and the y-axis direction at grid point (i, j), T 1 ＝min((T i-1,j ,T i+1,j )，T 2 ＝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 ,Ti,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 1 、T 2 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 ,Ti,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 ,Ti,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 ,fNs (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, τ 0 (j) Representing the time of onset of the jth microseismic event.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310730999.4A CN116660980A (en) | 2023-06-19 | 2023-06-19 | Microseism positioning method based on improved path function equation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310730999.4A CN116660980A (en) | 2023-06-19 | 2023-06-19 | Microseism positioning method based on improved path function equation |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116660980A true CN116660980A (en) | 2023-08-29 |
Family
ID=87727825
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310730999.4A Pending CN116660980A (en) | 2023-06-19 | 2023-06-19 | Microseism positioning method based on improved path function equation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116660980A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117607957A (en) * | 2024-01-24 | 2024-02-27 | 南方科技大学 | Seismic wave travel time solving method and system based on equivalent slowness rapid propulsion method |
-
2023
- 2023-06-19 CN CN202310730999.4A patent/CN116660980A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117607957A (en) * | 2024-01-24 | 2024-02-27 | 南方科技大学 | Seismic wave travel time solving method and system based on equivalent slowness rapid propulsion method |
| CN117607957B (en) * | 2024-01-24 | 2024-04-02 | 南方科技大学 | Seismic wave travel time solving method and system based on equivalent slowness rapid propulsion method |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106772577B (en) | Source inversion method based on microseism data and SPSA optimization algorithm | |
| CN106885576B (en) | AUV (autonomous Underwater vehicle) track deviation estimation method based on multipoint terrain matching positioning | |
| CN112328718B (en) | Road topology construction system and method based on vehicle dynamic trajectory tracking | |
| Hu et al. | Acoustic emission source location and experimental verification for two-dimensional irregular complex structure | |
| CN103136393B (en) | A kind of areal coverage computing method based on stress and strain model | |
| CN105241465B (en) | A kind of method of road renewal | |
| CN106597363A (en) | Pedestrian location method in indoor WLAN environment | |
| CN105093299A (en) | Observation system optimization method based on offset vector tile technology and apparatus thereof | |
| CN104807460A (en) | Indoor positioning method and system for unmanned aerial vehicle | |
| CN106249295B (en) | A kind of borehole microseismic P, S wave joint method for rapidly positioning and system | |
| CN108984818A (en) | Fixed-wing time domain aviation electromagnetic data intend restricted by three-dimensional space entirety inversion method | |
| CN105759311A (en) | Near-real time earthquake source position positioning method | |
| CN105246153A (en) | High-density rapid collection method for indoor fingerprint positioning database | |
| CN116660980A (en) | Microseism positioning method based on improved path function equation | |
| CN108845358B (en) | Tomography and the recognition methods of structural anomaly body and device | |
| CN106028446A (en) | Indoor parking lot location method | |
| CN104155694B (en) | A kind of residual static corrections reflecting converted shear wave common geophone stack section | |
| CN104507097A (en) | Semi-supervised training method based on WiFi (wireless fidelity) position fingerprints | |
| CN113960532B (en) | Microseism positioning method based on secondary positioning calculation of virtual source | |
| CN102706348B (en) | Gravimetric map fast matching method based on triangle | |
| CN109085642B (en) | Anisotropic medium microseism event positioning method | |
| CN109212594B (en) | Combined positioning method for longitudinal waves and transverse waves of anisotropic medium | |
| CN109387872B (en) | Surface multiple prediction method | |
| CN113568041B (en) | Repeatability analysis method and system for time-lapse seismic three-dimensional towing cable acquired data | |
| CN111189926B (en) | Method and system for identifying structure hole position based on global search |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |