CN110687598A - Method and device for accelerating microseismic numerical simulation - Google Patents

Method and device for accelerating microseismic numerical simulation Download PDF

Info

Publication number
CN110687598A
CN110687598A CN201911041762.5A CN201911041762A CN110687598A CN 110687598 A CN110687598 A CN 110687598A CN 201911041762 A CN201911041762 A CN 201911041762A CN 110687598 A CN110687598 A CN 110687598A
Authority
CN
China
Prior art keywords
seismic
point
travel time
namely
wave field
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
Application number
CN201911041762.5A
Other languages
Chinese (zh)
Inventor
王云宏
王保利
崔伟雄
赵朋朋
金丹
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Research Institute Co Ltd of CCTEG
Original Assignee
Xian Research Institute Co Ltd of CCTEG
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xian Research Institute Co Ltd of CCTEG filed Critical Xian Research Institute Co Ltd of CCTEG
Priority to CN201911041762.5A priority Critical patent/CN110687598A/en
Publication of CN110687598A publication Critical patent/CN110687598A/en
Pending legal-status Critical Current

Links

Images

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/282—Application of seismic models, synthetic seismograms
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Acoustics & Sound (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Business, Economics & Management (AREA)
  • Emergency Management (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention relates to a method and a device for simulating a microseismic numerical value, belongs to the field of microseismic monitoring, and particularly relates to a method and a device for accelerating the microseismic numerical simulation. The method combines the quantity characteristics of the seismic sources and the receiving detectors in the microseism monitoring, calculates the seismic wave field by taking the original seismic source points as the receiving points and the original receiving points as the seismic source points on the basis of the prior art according to the source detection interchange principle, can quickly obtain the wave field information of each receiving point in the microseism numerical simulation, and greatly improves the calculation efficiency compared with the prior art.

Description

Method and device for accelerating microseismic numerical simulation
Technical Field
The invention relates to a method and a device for simulating a microseismic numerical value, belongs to the field of microseismic monitoring, and particularly relates to a method and a device for accelerating the microseismic numerical simulation.
Background
The microseismic monitoring technology is a geophysical prospecting technology which is widely applied in the fields of petroleum and coal fields in recent years, and is widely applied to the fields of fracturing monitoring, rock burst monitoring and the like. The microseismic monitoring technology mainly relates to two technical problems of forward evolution and inversion. The microseismic monitoring forward modeling problem is a basic problem in microseismic monitoring and is also a basis for inversion, wherein the microseismic field (dynamics) and seismic wave travel time (kinematics) and other information are solved by knowing parameters such as a speed model, an observation system, a seismic source coordinate and the like, and a basis is provided for design of the microseismic monitoring observation system, correctness verification of a positioning method and the like. The microseismic monitoring inversion problem is a forward inverse process, namely a process of obtaining information such as a velocity model, seismic source parameters and the like by knowing information such as an observation system, a seismic wave field, travel time and the like.
The microseismic forward modeling method mainly comprises the following steps: finite difference method, finite element method, pseudo spectrum method, wave front fast propulsion method and the like based on wave equation theory, boundary integral method, boundary element method and the like based on integral equation, pilot method, bending method and the like based on ray tracing theory.
The basic flow of the conventional microseismic monitoring forward numerical simulation algorithm is as follows: 1. parameters such as a three-dimensional grid velocity model, a seismic source and a demodulator probe coordinate and the like of a simulation area are given; 2. calculating the seismic wave field (kinematics or dynamics) of each seismic source point in the three-dimensional model space; 3. information at each geophone point is extracted from the seismic wavefield. Since microseismic monitoring forward is a three-dimensional problem, the model space is large, and calculating the seismic wave field of the three-dimensional model space is a very time-consuming process, and the main influence factors of the process include the size of the model space, the grid division size and the like. If the model space is large and the grid division is fine, dozens of minutes or even hundreds of minutes are needed for calculating the three-dimensional seismic wave field of one seismic source point. Assuming that there are M microseismic sources, it takes N minutes to calculate 1 source, and the time to complete forward modeling of all the sources is approximately M × N minutes. Therefore, the step 2 is a main factor for limiting the microseismic forward modeling efficiency.
Disclosure of Invention
The invention mainly solves the technical problem of slow microseismic numerical simulation (forward operation) speed in the prior art, and provides a method and a device for accelerating microseismic numerical simulation. The method and the device are combined with the quantity characteristics of the seismic sources and the wave detection points in the microseismic monitoring, and the wave field information of each wave detection point in the microseismic numerical simulation can be quickly obtained according to the source detection interchange principle.
The technical problem of the invention is mainly solved by the following technical scheme:
a method for accelerating microseismic numerical simulation specifically comprises the following steps:
a forward modeling parameter setting step, namely establishing a grid speed model according to the geological condition of an object to be modeled, and specifying the coordinates of wave detection points;
and a wave field calculation step, namely determining a seismic wave field calculation mode according to the number Ns of the seismic source points and the number Nr of the wave detection points. If Ns is larger than Nr, performing source detection interchange and wave field calculation, namely taking each wave detection point as a seismic source point, and respectively calculating the point of each wave detection point as a seismic wave field of a seismic source; and if Ns is less than Nr, calculating the seismic wave field of each seismic source point according to a conventional method. The former is generally the case in microseismic monitoring.
And a wave field extraction step, namely extracting the seismic wave field of each seismic source point based on the three-dimensional seismic wave field obtained by calculation, and finishing forward modeling.
Preferably, the forward modeling parameter setting step performs the following substeps:
establishing an XYZ three-dimensional rectangular coordinate system according to a region to be simulated; meshing the area to be simulated based on the three-dimensional rectangular coordinate system; carrying out speed assignment on the grid points based on the grid division result, and establishing a speed model; the spatial location of each of the detection points is specified. Wherein, the grid velocity model, the space coordinates of the demodulator probe and the microseismic source point are respectively recorded as:
grid speed model: v. ofi,j,k(1≤i≤Nx,1≤j≤Ny,1≤k≤Nz)
Coordinate of wave detection point (xr)i,yri,zrri)(1≤i≤Nr)
Coordinates of the seismic source points: (xs)i,ysi,zsri)(1≤i≤Ns)
Description of variables: v. ofi,j,kVelocity, N, corresponding to mesh node (i, j, k)xIs the number of grids in the X direction of the model, NyIs the number of grids in the Y direction of the model, NzAs the number of grids in the Z direction of the model, xriIs the X coordinate of the i-th detection point, yriIs the Y coordinate of the ith demodulator probe, zriAs Z coordinate of the ith detection point, NrFor the number of detection points, xsiIs the X coordinate, ys, of the ith seismic source pointiIs the Y coordinate of the ith seismic source point, zsiIs the Z coordinate of the ith seismic source point, NsThe number of the seismic source points.
Preferably, the wavefield calculation module performs the sub-steps of:
firstly, determining a seismic wave field calculation mode according to the number Ns of seismic source points and the number Nr of wave detection points. And if Ns is less than Nr, traversing the seismic source points according to a conventional method, and if Ns is more than Nr, sequentially traversing each detection point, taking the detection point as the seismic source point, and marking the original seismic source point as the detection point.
Second, the seismic wavefield is calculated at the current seismic source point, under the known velocity model. Taking the fast-marching algorithm as an example, solving the equation of the equation under the following two-dimensional rectangular coordinate system (x-z):
Figure BDA0002253035820000031
where t is the seismic travel time, t (x, z) represents the seismic travel time at coordinate point (x, z), s is the slowness, i.e., the inverse of the velocity, s2(x, z) represents the square of the slowness value at coordinate point (x, z), which represents 2 coordinate axes of a two-dimensional rectangular coordinate system.
The windward difference format is adopted to replace a partial differential format in the equation of the equation, namely, the windward difference method is adopted to disperse the equation, and the dispersion formula can be expressed as follows:
Figure BDA0002253035820000032
d is a difference operator, i and j are grid node index numbers in the x direction and the z direction respectively, max is a maximum function, and min is a minimum function. The above equation can be simplified as:
Figure BDA0002253035820000041
Figure BDA0002253035820000042
represents the forward and backward difference operators of the travel time t in the x direction at the point (i, j),
Figure BDA0002253035820000043
represents the forward and backward difference operator of the travel time t in the z direction at point (i, j).
In order to ensure the calculation precision, a second-order upwind difference operator is adopted. Namely:
Figure BDA0002253035820000044
Figure BDA0002253035820000045
wherein Δ x and Δ z are respectively the grid size in the direction X, Z, ti,jRepresenting the value of the travel time at point (i, j).
Finally, the travel time values of the (i, j) point are calculated by 6 combination methods according to the travel time values of eight nodes (i, j ± 1), (i, j ± 2), (i ± 1, j), (i ± 2, j) around the target point (i, j), and the 6 combination methods have the following second-order difference calculation formulas:
Figure BDA0002253035820000046
in the formula, t is the seismic travel time, s is the slowness of the seismic, and h ═ Δ x ═ Δ z is the square grid size.
Calculating a plurality of travel times by using the formula, and selecting the minimum travel time point from the solved multi-value solution according to a narrow-band extension basic principle, wherein the minimum travel time point is the travel time of a target point, and the calculation formula is as follows:
ti,j＝min(t1 i,j,t2 i,j,t3 i,j,t4 i,j,t5 i,j,t6 i,j)
and traversing all the nodes to obtain the seismic wave time field.
Preferably, said wavefield extraction step performs the following sub-steps:
and (3) extracting the seismic wave field of each detection point, namely finding the time corresponding to the position of the detection point in the time field t (x, z) obtained in the previous step, namely the seismic travel time from the original seismic source point to the detection point, and finishing the forward modeling.
A device for accelerating microseismic numerical simulation specifically comprises the following steps:
and the forward modeling parameter setting device is used for establishing a grid speed model according to the geological condition of the object to be modeled and designating the coordinates of the wave detection points.
And the wave field calculation device determines the seismic wave field calculation mode according to the number Ns of the seismic source points and the number Nr of the wave detection points. If Ns is larger than Nr, performing source detection interchange and wave field calculation, namely taking each wave detection point as a seismic source point, and respectively calculating the point of each wave detection point as a seismic wave field of a seismic source; and if Ns is less than Nr, calculating the seismic wave field of each seismic source point according to a conventional method. The former is generally the case in microseismic monitoring.
And the wave field extraction device extracts the seismic wave field of each wave detection point based on the three-dimensional seismic wave field obtained by calculation, and the forward modeling is finished.
Preferably, the forward simulation parameter setting module performs the following sub-steps:
establishing an XYZ three-dimensional rectangular coordinate system according to a region to be simulated; carrying out meshing division on a region to be simulated based on the three-dimensional rectangular coordinate system, carrying out speed assignment on grid points based on the meshing division result, and establishing a speed model; and specifies the spatial location of each of the detection points. The grid velocity model, the spatial coordinates of the detection points and the potential seismic source points are respectively recorded as:
grid speed model: v. ofi,j,k(1≤i≤Nx,1≤j≤Ny,1≤k≤Nz)
Coordinate of wave detection point (xr)i,yri,zrri)(1≤i≤Nr)
Coordinates of the seismic source points: (xs)i,ysi,zsri)(1≤i≤Ns)
Description of variables: v. ofi,j,kVelocity, N, corresponding to mesh node (i, j, k)xIs the number of grids in the X direction of the model, NyIs the number of grids in the Y direction of the model, NzAs the number of grids in the Z direction of the model, xriIs the X coordinate of the i-th detection point, yriIs the Y coordinate of the ith demodulator probe, zriAs Z coordinate of the ith detection point, NrFor the number of detection points, xsiIs the X coordinate, ys, of the ith seismic source pointiIs the Y coordinate of the ith seismic source point, zsiIs the Z coordinate of the ith seismic source point, NsThe number of the seismic source points.
Preferably, the wavefield calculation module performs the sub-steps of:
firstly, determining a seismic wave field calculation mode according to the number Ns of seismic source points and the number Nr of wave detection points. And if Ns is larger than Nr, traversing the seismic source points according to a conventional method, and if Ns is smaller than Nr, sequentially traversing each detection point, taking the detection point as the seismic source point, and marking the original seismic source point as the detection point.
Second, the seismic wavefield is calculated at the current seismic source point, under the known velocity model. Taking the fast-marching algorithm as an example, solving the equation of the equation under the following two-dimensional rectangular coordinate system (x-z):
Figure BDA0002253035820000061
where t is the seismic travel time, s is the slowness, i.e., the reciprocal of the velocity, and x, z represent 2 coordinate axes of a two-dimensional rectangular coordinate system.
The windward difference format is adopted to replace a partial differential format in an equation function equation, namely, an windward difference method is adopted to disperse the above formula, and the dispersion formula can be expressed as follows:
Figure BDA0002253035820000062
d is a difference operator, i and j are grid node index numbers in the x direction and the z direction respectively, max is a maximum function, and min is a minimum function. The above equation can be simplified as:
Figure BDA0002253035820000063
Figure BDA0002253035820000064
represents the forward and backward difference operators of the travel time t in the x direction at the point (i, j),
Figure BDA0002253035820000065
represents the forward and backward difference operator of the travel time t in the z direction at point (i, j).
In order to ensure the calculation precision, a second-order upwind difference operator is adopted. Namely:
Figure BDA0002253035820000071
Figure BDA0002253035820000072
in the formula, Δ x and Δ z are grid sizes in the X, Z direction, respectively.
Finally, the travel time values of the (i, j) point are calculated by 6 combination methods according to the travel time values of eight nodes (i, j ± 1), (i, j ± 2), (i ± 1, j), (i ± 2, j) around the target point (i, j), and the 6 combination methods have the following second-order difference calculation formulas:
Figure BDA0002253035820000073
in the formula, t is the seismic travel time, s is the slowness of the seismic, and h ═ Δ x ═ Δ z is the square grid size.
Calculating a plurality of travel times by using the formula, and selecting the minimum travel time point from the solved multi-value solution according to a narrow-band extension basic principle, wherein the minimum travel time point is the travel time of a target point, and the calculation formula is as follows:
ti,j＝min(t1 i,j,t2 i,j,t3 i,j,t4 i,j,t5 i,j,t6 i,j)
and traversing all the nodes to obtain the seismic wave time field.
Preferably, the wavefield extraction device performs the following sub-steps:
and (3) extracting the seismic wave field of each detection point, namely finding the time corresponding to the position of the detection point in the time field t (x, z) obtained in the previous step, namely the seismic travel time from the original seismic source point to the detection point, and finishing the forward modeling.
Therefore, the invention has the following advantages:
according to the source detection interchange principle in the seismic wave field simulation, the original seismic source points are used as the detection points and the original detection points are used as the seismic source points on the basis of the prior art, the seismic wave field is calculated, and finally the seismic wave field information of each source detection pair is obtained by combining the characteristic that the number of the seismic source points in the microseismic monitoring is far larger than that of the detection points. After the source detection and the exchange, the microseismic numerical simulation efficiency is greatly improved.
Drawings
FIG. 1 is a flow chart of the present invention.
FIG. 2 is an embodiment velocity model (two-dimensional) and observation system.
Detailed Description
The technical scheme of the invention is further specifically described by the following embodiments and the accompanying drawings. The scope of the invention is not limited by the following examples.
According to the symmetric sampling principle proposed by Vermeer, the source point (or shot point) and the demodulator probe are interchanged in seismic exploration, and the same seismic signal can be obtained, namely, the signal is received at the excitation-demodulator probe at the source point and the signal is received at the excitation-source probe at the demodulator probe, and the signals received by the two are completely equivalent, and the principle is called as the source detection interchange principle (namely, the principle of shot detection interchange).
Based on the principle, firstly, a grid speed model is established according to the geological condition of an object to be simulated, and wave detection point coordinates are specified; then, determining a wave field calculation mode according to the number of the seismic source points and the number of the detection points, and calculating the seismic wave field of the model by sequentially taking each detection point as the seismic source point based on a source detection interchange principle when the number of the seismic source points is less than the number of the detection points; and finally, extracting the seismic wave field of each wave detection point (original seismic source point) based on the three-dimensional seismic wave field obtained by the calculation, and finishing forward modeling.
The method for simulating the acceleration microseismic value according to the present embodiment will be further described with reference to fig. 1.
As shown in fig. 1, the present invention adopts the following technical solutions:
firstly, forward modeling parameters are given, namely a grid speed model is established according to geological conditions of an object to be modeled, and wave detection point coordinates are designated.
Specifically, in the present embodiment, an XYZ three-dimensional rectangular coordinate system is established according to the region to be simulated; meshing the area to be simulated based on the three-dimensional rectangular coordinate system; carrying out speed assignment on the grid points based on the grid division result, and establishing a speed model; the spatial location of each of the detection points is specified. Wherein, the grid velocity model, the space coordinates of the demodulator probe and the microseismic source point are respectively recorded as:
grid speed model: v. ofi,j,k(1≤i≤Nx,1≤j≤Ny,1≤k≤Nz)
Coordinate of wave detection point (xr)i,yri,zrri)(1≤i≤Nr)
Coordinates of the seismic source points: (xs)i,ysi,zsri)(1≤i≤Ns)
Description of variables: v. ofi,j,kVelocity, N, corresponding to mesh node (i, j, k)xIs the number of grids in the X direction of the model, NyIs the number of grids in the Y direction of the model, NzAs the number of grids in the Z direction of the model, xriIs the X coordinate of the i-th detection point, yriIs the Y coordinate of the ith demodulator probe, zriAs Z coordinate of the ith detection point, NrFor the number of detection points, xsiIs the X coordinate, ys, of the ith seismic source pointiIs the Y coordinate of the ith seismic source point, zsiIs the Z coordinate of the ith seismic source point, NsThe number of the seismic source points.
Preferably, the wavefield calculation module performs the sub-steps of:
firstly, determining a seismic wave field calculation mode according to the number Ns of seismic source points and the number Nr of wave detection points. And if Ns is larger than Nr, traversing the seismic source points according to a conventional method, and if Ns is smaller than Nr, sequentially traversing each detection point, taking the detection point as the seismic source point, and marking the original seismic source point as the detection point. The former is generally the case in microseismic monitoring.
Second, the seismic wavefield is calculated at the current seismic source point, under the known velocity model. Taking the fast-marching algorithm as an example, solving the equation of the equation under the following two-dimensional rectangular coordinate system (x-z):
Figure BDA0002253035820000091
where t is the seismic travel time, s is the slowness, i.e., the reciprocal of the velocity, and x, z represent 2 coordinate axes of a two-dimensional rectangular coordinate system.
The windward difference format is adopted to replace a partial differential format in an equation function equation, namely, an windward difference method is adopted to disperse the above formula, and the dispersion formula can be expressed as follows:
Figure BDA0002253035820000101
d is a difference operator, i and j are grid node index numbers in the x direction and the z direction respectively, max is a maximum function, and min is a minimum function. The above equation can be simplified as:
Figure BDA0002253035820000102
represents the forward and backward difference operators of the travel time t in the x direction at the point (i, j),represents the forward and backward difference operator of the travel time t in the z direction at point (i, j).
In order to ensure the calculation precision, a second-order upwind difference operator is adopted. Namely:
Figure BDA0002253035820000105
Figure BDA0002253035820000106
in the formula, Δ x and Δ z are grid sizes in the X, Z direction, respectively.
Finally, the travel time values of the (i, j) point are calculated by 6 combination methods according to the travel time values of eight nodes (i, j ± 1), (i, j ± 2), (i ± 1, j), (i ± 2, j) around the target point (i, j), and the 6 combination methods have the following second-order difference calculation formulas:
Figure BDA0002253035820000107
in the formula, t is the seismic travel time, s is the slowness of the seismic, and h ═ Δ x ═ Δ z is the square grid size.
Calculating a plurality of travel times by using the formula, and selecting the minimum travel time point from the solved multi-value solution according to a narrow-band extension basic principle, wherein the minimum travel time point is the travel time of a target point, and the calculation formula is as follows:
ti,j＝min(t1 i,j,t2 i,j,t3 i,j,t4 i,j,t5 i,j,t6 i,j)
and traversing all the nodes to obtain the seismic wave time field.
Preferably, said wavefield extraction step performs the following sub-steps:
and then, extracting the seismic wave field of each detection point, namely finding the time corresponding to the position of the detection point in the time field t (x, z) obtained in the previous step, namely the seismic travel time from the original seismic source point to the detection point, and finishing forward evolution.
The method greatly improves the efficiency of the microseismic monitoring forward modeling. In microseismic monitoring, the number of seismic sources is far greater than the number of detection points (Ns > Nr), for example, in coal field horizontal well fracturing microseismic monitoring, the number of detection points (Nr) is generally dozens, while microseismic events (each event corresponds to one source point) (Ns) are often hundreds or thousands, and the difference between the detection points and the microseismic events is 1-2 orders of magnitude. When the forward modeling is performed on the process, if the seismic wave field of each source point is obtained according to the existing method, and the modeling time of one source point is assumed to be T, the time required for modeling the whole process is Ns x T. According to the method, after source detection is exchanged, only the seismic wave field of each detection point is required, and because the velocity model is unchanged, the simulation time of each detection point is also T, and the time required for the whole simulation process is Nr multiplied by T; the time difference between the two schemes is (Ns-Nr) multiplied by T, and it can be seen that the larger the difference between Ns and Nr, the larger the acceleration ratio.
In specific implementation, the difference between Ns and Nr needs to be compared, if Ns > Nr, each geophone point seismic wavefield is calculated, and if Ns < Nr, each seismic source point seismic wavefield is calculated, the process is shown in fig. 1, and microseismic monitoring forward modeling is generally suitable for the former.
The forward modeling method has the contrast effect with the conventional forward modeling method. The example is only illustrated in a two-dimensional case, and the specific physical problem is described as follows: 1 wave detection point is respectively placed at 4 vertexes of a square model with the side length of 100m, the total number of the wave detection points is 4, the model speed is 1000m/s at a constant speed, the grid size is 0.1m x 0.1m, the grid size is 1001 x 1001, 10 seismic sources are randomly arranged in the square, and the travel time of each source detection pair is calculated. The operation time of the method is 3058ms, the operation time of the conventional method is 7898ms, and the acceleration ratio is 2.58 times.
The embodiment simply illustrates the improvement of the invention on the calculation efficiency only by the two-dimensional condition, and in practical application, the invention is a three-dimensional model, the number of the seismic sources is far greater than that of the demodulator probes, and the acceleration ratio is larger.
The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.

Claims (5)

1. A method of accelerating microseismic numerical simulation comprising:
a forward parameter setting step, namely establishing a grid speed model according to the geological condition of an object to be simulated and setting the coordinates of a wave detector;
a time field calculation step, namely determining a seismic wave field calculation mode according to the number Ns of seismic source points and the number Nr of wave detection points; if Ns is larger than Nr, performing source detection interchange and wave field calculation, namely taking each wave detection point as a seismic source point, and respectively calculating the point of each wave detection point as a seismic wave field of a seismic source;
and a wave field extraction step, namely extracting the seismic wave field of each wave detection point, and finishing forward modeling.
2. The method according to claim 1, wherein the time field calculating step specifically comprises:
a seismic wave field determining sub-step, namely calculating a path function equation of the seismic wave field of the current seismic source point under a known velocity model;
a substep of equation calculation of an equation, namely discretizing the equation of the equation by adopting an upwind difference method;
the earthquake travel time operator step, which is to calculate a second-order difference calculation formula of the earthquake travel time by using travel time values at nodes in eight directions around a target point to obtain a multi-valued solution of the travel time of the target point;
and a target travel time determining sub-step, namely selecting the minimum travel time point from the multi-value solution of the calculated target travel time according to the narrow-band extension basic principle, namely the target travel time.
3. A method of accelerating a microseismic numerical simulation of the type defined in claim 2 wherein in the seismic wavefield determination substep the equation of the path function of the seismic wavefield is determined based on the equation:
Figure FDA0002253035810000011
where t is the seismic travel time, s is the slowness, i.e., the reciprocal of the velocity, and x, z represent 2 coordinate axes of a two-dimensional rectangular coordinate system.
4. The method according to claim 2, wherein in the earthquake travel time calculation operator step, the above formula is discretized by an upwind difference method, specifically:
Figure FDA0002253035810000021
represents the forward and backward difference operators of the travel time t in the x direction at the point (i, j),
Figure FDA0002253035810000023
a difference operator representing the travel time t forward and backward in the z direction at point (i, j);
wherein:
Figure FDA0002253035810000024
Figure FDA0002253035810000025
in the formula, Δ x and Δ z are grid sizes in the X, Z direction, respectively.
5. The method according to claim 2, wherein the seismic travel time operator step calculates a multi-valued solution of the travel time of the target point using travel time values at eight nodes around the target point based on the following formula:
Figure FDA0002253035810000026
where t is the seismic travel time, ti,jIs the seismic travel time at point (i, j), s is the slowness of the seismic, h ═ Δ x ═ Δ z is the square grid size,
calculating a plurality of travel times by using the formula, and selecting the minimum travel time point from the solved multi-value solution according to a narrow-band extension basic principle, wherein the minimum travel time point is the travel time of a target point, and the calculation formula is as follows:
ti,j＝min(t1 i,j,t2 i,j,t3 i,j,t4 i,j,t5 i,j,t6 i,j)。
CN201911041762.5A 2019-10-30 2019-10-30 Method and device for accelerating microseismic numerical simulation Pending CN110687598A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911041762.5A CN110687598A (en) 2019-10-30 2019-10-30 Method and device for accelerating microseismic numerical simulation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911041762.5A CN110687598A (en) 2019-10-30 2019-10-30 Method and device for accelerating microseismic numerical simulation

Publications (1)

Publication Number Publication Date
CN110687598A true CN110687598A (en) 2020-01-14

Family

ID=69114602

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911041762.5A Pending CN110687598A (en) 2019-10-30 2019-10-30 Method and device for accelerating microseismic numerical simulation

Country Status (1)

Country Link
CN (1) CN110687598A (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112462415A (en) * 2020-11-02 2021-03-09 中国电子科技集团公司第三研究所 Method and device for positioning multiple vibration sources
CN113640878A (en) * 2021-08-12 2021-11-12 西南石油大学 Method for constructing azimuth-apparent velocity radar map by using virtual seismic source scanning

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160306058A1 (en) * 2015-04-16 2016-10-20 The Board Of Trustees Of The Leland Stanford Junior Univerisity Reverse time migration based on geometric mean for imaging seismic sources
CN109100784A (en) * 2018-06-08 2018-12-28 恒泰艾普（北京）能源科技研究院有限公司 All-wave field imaging method is exchanged in the inspection of three-dimensional VSP source
CN110018516A (en) * 2019-05-07 2019-07-16 西安石油大学 A kind of decoupling wave field microseism inverse time interference localization method
CN110261900A (en) * 2019-06-10 2019-09-20 中北大学 A kind of underground shallow layer microseism positioning system based on velocity information

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160306058A1 (en) * 2015-04-16 2016-10-20 The Board Of Trustees Of The Leland Stanford Junior Univerisity Reverse time migration based on geometric mean for imaging seismic sources
CN109100784A (en) * 2018-06-08 2018-12-28 恒泰艾普（北京）能源科技研究院有限公司 All-wave field imaging method is exchanged in the inspection of three-dimensional VSP source
CN110018516A (en) * 2019-05-07 2019-07-16 西安石油大学 A kind of decoupling wave field microseism inverse time interference localization method
CN110261900A (en) * 2019-06-10 2019-09-20 中北大学 A kind of underground shallow layer microseism positioning system based on velocity information

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
刘洋: "基于快速推进算法的射线追踪研究", 《中国优秀硕士学位论文全文数据库•基础科学辑》 *
葛奇鑫: "基于逆时成像的被动源定位与识别方法研究", 《中国博士学位论文全文数据库•基础科学辑》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112462415A (en) * 2020-11-02 2021-03-09 中国电子科技集团公司第三研究所 Method and device for positioning multiple vibration sources
CN113640878A (en) * 2021-08-12 2021-11-12 西南石油大学 Method for constructing azimuth-apparent velocity radar map by using virtual seismic source scanning
CN113640878B (en) * 2021-08-12 2024-03-29 西南石油大学 Method for constructing azimuth-apparent velocity radar chart by utilizing virtual seismic source scanning

Similar Documents

Publication Publication Date Title
US8655632B2 (en) Gridless geological modeling
CN110031896B (en) Seismic random inversion method and device based on multi-point geostatistics prior information
EP2869096B1 (en) Systems and methods of multi-scale meshing for geologic time modeling
US20110213600A1 (en) Method and system for using multiple-point statistics simulation to model reservoir property trends
CN104504754B (en) Oil and gas reservoir multipoint statistical modeling method and device
GB2474740A (en) Gridless geological modeling of a structural framework
CN112231974B (en) Deep learning-based method and system for recovering seismic wave field characteristics of rock breaking seismic source of TBM (Tunnel boring machine)
WO2014099201A1 (en) Geophysical modeling of subsurface volumes based on horizon extraction
CN109459787B (en) coal mine underground structure imaging method and system based on seismic channel wave full-waveform inversion
CN110687598A (en) Method and device for accelerating microseismic numerical simulation
CN110827405A (en) Digital remote sensing geological mapping method and system
CN117852416B (en) Multimode grouting precontrolled analysis method and system based on digital geological model
CN110850469A (en) Imaging method for seismic channel wave depth migration based on kirchhoff product decomposition
CN108680968B (en) Evaluation method and device for seismic exploration data acquisition observation system in complex structural area
CN118015219A (en) Geological model generation method, device and equipment based on qualitative kriging interpolation
CN110610539A (en) Stratum curved surface construction method, device, equipment and storage medium
CN107817524A (en) The method and apparatus of three-dimensional seismic tomography
CN116524140A (en) Three-dimensional geological modeling system
CN108986217B (en) Multipoint geostatistical modeling method based on pattern vector distance
CN112946760B (en) Non-explosive three-dimensional imaging method, device and system based on regularization method
CN110532706B (en) Dam foundation rock mass blasting relaxation analysis method based on GOCAD
CN108537883B (en) Geological modeling method based on MapReduce framework
CN112346115A (en) Micro-seismic source positioning method under complex rock mass wave velocity environment with cavities in underground chamber group
CN110954956B (en) Method for evaluating acquisition trace of observation system and computer-readable storage medium
CN111815769A (en) Modeling method, computing device and storage medium for thrust-driven tectonic belt structure

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
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20200114