CN113885076A - Microseism ground monitoring speed model correction method - Google Patents

Microseism ground monitoring speed model correction method Download PDF

Info

Publication number
CN113885076A
CN113885076A CN202111159206.5A CN202111159206A CN113885076A CN 113885076 A CN113885076 A CN 113885076A CN 202111159206 A CN202111159206 A CN 202111159206A CN 113885076 A CN113885076 A CN 113885076A
Authority
CN
China
Prior art keywords
perforation
theta
value
model
positioning
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
CN202111159206.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.)
Jilin University
Original Assignee
Jilin University
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 Jilin University filed Critical Jilin University
Priority to CN202111159206.5A priority Critical patent/CN113885076A/en
Publication of CN113885076A publication Critical patent/CN113885076A/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/288—Event detection in seismic signals, e.g. microseismics
    • 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/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/10—Aspects of acoustic signal generation or detection
    • G01V2210/14—Signal detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/59—Other corrections
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/66—Subsurface modeling

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)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

The invention relates to a microseism ground monitoring speed model correction algorithm, which comprises the steps of arranging N detectors on the ground, defining a three-dimensional target area, establishing an initial speed model, reading N seismic wave data acquired by the N detectors, calculating an energy focus value of a perforation position by adopting an amplitude superposition method, adjusting parameters of a target layer by adopting an extremely-fast simulation annealing method to obtain a new speed model, recalculating the energy focus value of the perforation position, judging whether to perform tempering treatment to adjust the speed model or adopt a grid successive division method to relocate the perforation, calculating perforation relocation errors, judging whether to perform simulated annealing cooling treatment according to whether the positioning errors meet the positioning requirements, storing an optimal speed model, perforation relocation results and errors, and finishing. The correction algorithm can more accurately reposition the perforation event to the true value, and can effectively improve the positioning reliability of the microseism event near the perforation point.

Description

Microseism ground monitoring speed model correction method
Technical Field
The invention belongs to the technical field of microseism monitoring, and particularly relates to a microseism ground monitoring speed model correction algorithm.
Background
The microseism monitoring technology is characterized in that a vibration signal generated by fracturing is picked up from a distributed microseism detector, and positioning processing technologies such as superposition energy scanning analysis and the like are used for providing the space form of a crack, monitoring the scale, the shape and the fracture network structure of an artificial reservoir formed by hydraulic fracturing, guiding the effective development of fracturing work and improving the productivity. Microseismic location techniques are central to microseismic monitoring efforts and are implemented such that crack propagation causes the surrounding rock to fracture, thereby triggering a series of observably recorded microseismic events. The velocity model in the monitoring work area range is a main factor influencing the positioning accuracy of the micro-seismic event, so that how to obtain an effective velocity model is a key problem in the micro-seismic monitoring project.
At present, the correction of the micro-seismic ground monitoring velocity model is based on a horizontal layered stratum model, and the velocity model is simple, so that the problem that a final perforation positioning result has certain errors is solved. Chinese patent CN107132578B discloses a microseism ground monitoring velocity model correction algorithm, which is based on an amplitude stacking microseism velocity model construction method and provides a method for improving perforation repositioning precision by expanding solution space, namely, the uncertainty of each layer velocity and the uncertainty of each layer position in the stratum are simultaneously considered in the process of an extremely-fast simulated annealing method. However, the method is based on a horizontal layered stratum model, and the velocity model is simple, so that a final perforation repositioning result has certain errors, and particularly when the actual stratum is complex, the repositioning error of the perforation cannot meet the requirements, and the velocity model required by positioning of a subsequent micro-seismic event cannot be obtained. In addition, the correction algorithm carries out perforation relocation every time a new speed model is obtained by adjusting parameters, and the positioning process is time-consuming. The positioning method is mesh division, only once division is carried out according to set division accuracy, if the size of the divided meshes is small, the positioning accuracy is high, the number of the meshes is very large, and the calculation time is long; if the grid size is large, the calculation time is short, but the positioning accuracy is low, so that higher calculation efficiency and positioning accuracy cannot be obtained at the same time.
Therefore, how to make up for the existing defects, improve the perforation repositioning accuracy under complex geological conditions, obtain a velocity model with stronger equivalence, and further improve the positioning accuracy of the microseism event is a problem to be solved urgently in the field.
Disclosure of Invention
The invention aims to provide a microseism ground monitoring velocity model correction algorithm to solve the problems of improving perforation repositioning precision under complex geological conditions, obtaining a velocity model with stronger equivalence and further improving microseism event positioning precision. On the basis of the horizontal layered model, the velocity model is complicated, and the change of velocity values of all layers, the change of interface positions of the layers and the change of rotation angles of all layers are considered, so that the solution space is further expanded, and the perforation repositioning precision is improved.
The purpose of the invention is realized by the following technical scheme:
a microseism ground monitoring velocity model correction algorithm comprises the following steps:
a. arranging N detectors on the ground, and defining a three-dimensional target area by taking the position of a perforation point as a center;
b. establishing an initial velocity model according to the logging curve and reading N seismic wave data acquired by the N detectors;
c. calculating the energy focus value E (V) of the perforation position by adopting an amplitude superposition method0,H0,θ0)；
d. Taking the speed value, the level value and the rotation angle value of the target layer as uncertain factors, and adjusting the parameters of the target layer by adopting a very fast simulated annealing method to obtain a new speed model;
e. recalculating the energy focus value E (V, H, theta) of the perforation position if E (V, H, theta)<E(V0,H0,θ0) Tempering to adjust the speed model; if E (V, H, theta) > E (V)0，H0，θ0) Repositioning the perforation by adopting a grid successive division method, and calculating perforation repositioning errors;
f. if the positioning error does not meet the positioning requirement, carrying out simulated annealing cooling treatment, if the temperature after cooling is higher than the set lowest temperature, carrying out tempering treatment to adjust the speed model, otherwise, storing the optimal speed model, the perforation repositioning result and the error, and then ending; and if the positioning error meets the positioning requirement, storing the optimal speed model, the perforation repositioning result and the error, and ending.
Further, in step b, the specific steps of establishing the initial velocity model are as follows:
obtaining an initial velocity vector V, V ═ V from the logp1，Vp2，Vp3，…，Vpn]Initial horizon vector H, H ═ H1，H2，H3，…，Hn-1]Initial rotation angle vector θx、θy，θx＝[θx1，θxx，θx3，…，θxn-1]，θy＝[θy1，θy2，θy3，…，θyn-1]Wherein V ispiVelocity of P-wave of i-th layer, HiIs the depth coordinate of the ith layer interface, θxi、θyiThe values of the angles of the horizontal plane rotating around the x axis and the y axis to the ith layer interface position are respectively positive anticlockwise and negative clockwise.
Further, in step c, the amplitude superposition method specifically includes the following steps:
c1, selecting any one channel in the perforation data as a reference channel M, and solving the theoretical travel time difference of other channels relative to the reference channel by adopting a ray tracing method:
Δt＝[t1-tM，t2-tM，t3-tM，…，tN-tM]
in the formula, t1，t2，t3，…，tNThe theoretical travel time t of each path except the reference path in the perforation dataMThe theoretical travel time of a reference track;
c2, performing reverse time migration superposition on the perforation record waveform obtained by each detector according to the theoretical travel time difference, wherein the mathematical expression of the objective function is as follows:
Figure BDA0003289480040000031
wherein A (i, j) is the amplitude of the ith channel data at the jth moment, N is the number of detectors, and W is the length of a time window.
Further, in step d, the extremely fast simulated annealing method specifically comprises the following steps:
d1, calculating the simulated annealing initial temperature T0The solution to obtain the initial temperature is as follows:
giving a small positive value to the initial temperature, and then continuously multiplying the initial temperature by a number which is constantly larger than 1 until the acceptance probability of any model is close to 1;
Figure BDA0003289480040000041
wherein, V0，H0，θ0As a parameter of the initial velocity model, V1，H1，θ1Is the model result of the first iteration;
d2, simulating annealing cooling calculation, wherein the cooling formula is as follows:
Tk＝T0exp(-ck1/2n)
where k is the number of iterations, T0C is a given constant for the initial temperature, and the cooling process is determined, wherein c is 0.5, and n is required to be adjustedThe number of layers of the whole stratum;
d3, in the course of simulated annealing calculation, each iteration needs to adjust the velocity vector V, horizon vector H and rotation angle vector theta in the stratum modelx、θyThe adjustment formula is as follows:
Figure BDA0003289480040000042
Figure BDA0003289480040000043
Figure BDA0003289480040000044
Figure BDA0003289480040000045
wherein,
Figure BDA0003289480040000046
the maximum and minimum values of the speed adjustment range for each layer,
Figure BDA0003289480040000047
Figure BDA0003289480040000048
Figure BDA0003289480040000049
the maximum and minimum values of the adjustment range for each layer position,
Figure BDA00032894800400000410
Figure BDA00032894800400000411
the maximum value and the minimum value of the adjustment range of the interfaces of the layers and the rotation angle are obtained,
Figure BDA00032894800400000412
x1、x2、x3、x4is a random variable with the value range of [ -1,1]Generating x1，2，3，4The expression of (a) is as follows:
Figure BDA0003289480040000051
wherein sgn is a sign function, μ ∈ [0,1 ];
d4, in each iteration, if E (V ', H', theta ') > E (V, H, theta), V', H ', theta' replaces V, H, theta as the current optimal solution, if E (V ', H', theta ') < E (V, H, theta), then the V', H ', theta' replaces V, H, theta as the current optimal solution with probability P (V, H, theta → V ', H', theta '), and the calculation formula of P (V, H, theta → V', H ', theta') is as follows:
Figure BDA0003289480040000052
wherein V ', H ', theta ' are velocity models in the iterative process, V, H, theta are the current optimal solution of the velocity models, and TkThe temperature value of the k iteration is the temperature value of the simulated annealing, and the iteration termination condition of the simulated annealing is that the positioning error meets the positioning requirement or the temperature TkBelow a set minimum temperature value.
Further, in step e, the mesh successive-subdivision positioning method specifically includes the following steps:
e1, setting the radius r of the target area, the center P of the target area as the perforation position, and the mesh division precision dinRequirement for positioning accuracy doutPositioning result Lout(0,0,0), energy focus maximum Eout＝0；
E2, performing mesh division on the target area, calculating the energy focus value of the central point of each mesh, and recording the maximum energy focus value as EmaxAnd the center point of the grid is Lmax；
E3, if Emax＞EoutThen order Eout＝Emax，Lout＝Lmax；
e4, if din＞doutIf r is equal to r/2, din＝din/2，P＝LoutRepeating the second step; if d isin＜doutThen store LoutAnd the final perforation positioning result.
Compared with the prior art, the invention has the beneficial effects that: on the basis of the horizontal layered model, the invention complicates the speed model, and considers the change of the layer interface position, the change of the rotation angle of each layer and the change of the speed value of each layer at the same time, so as to further expand the solution space and improve the perforation repositioning precision. The test result proves that compared with the prior method, the method can more accurately reposition the perforation event to the true value and can effectively improve the positioning reliability of the microseism event near the perforation point.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings needed to be used in the embodiments will be briefly described below, it should be understood that the following drawings only illustrate some embodiments of the present invention and therefore should not be considered as limiting the scope, and for those skilled in the art, other related drawings can be obtained according to the drawings without inventive efforts.
FIG. 1 is a flow chart of the method of the present invention;
FIG. 2a is a diagram of a three-dimensional layered relief structure stratigraphic model;
FIG. 2b a fifth layer formation interface;
FIG. 3 is a forward analog waveform diagram for each channel of the detector;
FIG. 4a is a ray tracing forward three-dimensional view;
FIG. 4b is a schematic normal vertical slice of the ray tracing;
5 a-5 b method perforation repositioning results before improvement;
fig. 5 c-5 d illustrate the perforation repositioning result of the method of the present invention.
Detailed Description
The present invention will be further described with reference to specific embodiments, which are implemented on the premise of the technology of the present invention, and detailed embodiments are given, but the scope of the present invention is not limited to the following examples.
As shown in FIG. 1, the microseism ground monitoring velocity model correction algorithm comprises the following steps:
a. in this embodiment, a 9-layer layered relief stratum model (as shown in fig. 2 a) is selected, and 8 curved surfaces are used as horizon interfaces (the curved surfaces of the interfaces are not parallel to each other, but the relief forms are all shown in fig. 2 b). In the local coordinate system, the designed perforation positions are (-280,360, -3550) (unit: m). The microseism ground detectors are arranged in a star shape, 6 measuring lines are provided, 10 detectors (60 lines) are arranged in each measuring line, the distance between the detectors is 50m (as shown in figure 2a, a blue cross is the position of the detector, and a red five-pointed star is the position of a perforation). The actual signal was simulated using the finite difference wave equation, and the waveform was described using a 20Hz Rake wavelet with a maximum amplitude of 1 for each Rake wavelet (as shown in FIG. 3). Subsequent velocity model correction schemes are evolving to employ ray tracing methods (as shown in fig. 4).
b. An initial velocity model is obtained from the logging information (shown as a solid black line in fig. 2 a) with an initial velocity vector V ═ 900, 1300, 1800, 2400, 2800, 3600, 4200, 4800,5400](unit: m/s), initial horizon vector H [ -200, -400, -700, -1000, -1400, -1900, -2500, -3200](unit: m), initial rotation angle vector θx＝θy＝[0，0，0，0，0，0，0，0](unit: °). In the process of correcting the speed model, the adjustment range of the speed of each layer is set to 600m/s, the adjustment range of the horizon value of each layer is set to 100m, the adjustment range of the rotation angle of each layer is set to 6 degrees, and the specific adjustment range of the parameters of each layer is shown in table 1.
c. Reading the analog signals of 60 detectors, and calculating the energy focus value E (V) of the perforation position under the condition of an initial velocity model by adopting an amplitude superposition method0，H0，θ0) 713.882, the specific steps of the amplitude superposition method are as follows:
c1, selecting any one channel in the perforation data as a reference channel M, and solving the theoretical travel time difference of other channels relative to the reference channel by adopting a ray tracing method:
Δt＝[t1-tM，t2-tM，t3-tM，…，tN-tM]
in the formula, t1，t2，t3，…，tNThe theoretical travel time t of each path except the reference path in the perforation dataMThe theoretical travel time of a reference track;
c2, performing reverse time migration superposition on the perforation record waveform obtained by each detector according to the theoretical travel time difference, wherein the mathematical expression of the objective function is as follows:
Figure BDA0003289480040000071
wherein A (i, j) is the amplitude of the ith data at the jth moment, N is the number of detectors, and W is the length of a time window;
d. and adjusting the speed value, the layer position value and the rotation angle value of the target layer by adopting an ultra-fast simulated annealing algorithm to obtain a new speed model. And (3) the stratum equivalent model is complicated into an inclined stratum model which is closer to the actual stratum and can be independently adjusted in each layer speed value, layer position value and interface rotation angle.
d1, calculating the simulated annealing initial temperature T0The solution to obtain the initial temperature is as follows:
giving a small positive value to the initial temperature, and then continuously multiplying the initial temperature by a number which is constantly larger than 1 until the acceptance probability of any model is close to 1;
Figure BDA0003289480040000081
wherein, V0，H0，θ0As a parameter of the initial velocity model, V1，H1，θ1As a result of the model of the first iteration(ii) a An approximate value of 200 ℃ was obtained for the initial temperature in this example.
d2, simulating annealing cooling calculation, wherein the cooling formula is as follows:
Tk＝T0exp(-ck1/2n)
where k is the number of iterations, T0C is a given constant, the cooling process is determined, n is the number of layers of the stratum to be adjusted, and k is 500, c is 0.5, and n is 9 in the calculation example;
d3, in the course of simulated annealing calculation, each iteration needs to adjust the velocity vector V, horizon vector H and rotation angle vector theta in the stratum modelx、θyThe adjustment formula is as follows:
Figure BDA0003289480040000082
Figure BDA0003289480040000083
Figure BDA0003289480040000084
Figure BDA0003289480040000085
wherein,
Figure BDA0003289480040000086
the maximum and minimum values of the speed adjustment range for each layer,
Figure BDA0003289480040000087
Figure BDA0003289480040000088
Figure BDA0003289480040000089
the maximum and minimum values of the adjustment range for each layer position,
Figure BDA0003289480040000091
Figure BDA0003289480040000092
the maximum value and the minimum value of the adjustment range of the interfaces of the layers and the rotation angle are obtained,
Figure BDA0003289480040000093
x1、x2、x3、x4is a random variable with the value range of [ -1,1]Generating x1，2，3，4The expression of (a) is as follows:
Figure BDA0003289480040000094
wherein sgn is a sign function, μ ∈ [0,1 ];
d4, in each iteration, if E (V ', H', theta ') > E (V, H, theta), V', H ', theta' replaces V, H, theta as the current optimal solution, if E (V ', H', theta ') < E (V, H, theta), then the V', H ', theta' replaces V, H, theta as the current optimal solution with probability P (V, H, theta → V ', H', theta '), and the calculation formula of P (V, H, theta → V', H ', theta') is as follows:
Figure BDA0003289480040000095
wherein V ', H ', theta ' are velocity models in the iterative process, V, H, theta are the current optimal solution of the velocity models, and TkThe temperature value of the k iteration is the temperature value of the simulated annealing, and the iteration termination condition of the simulated annealing is that the positioning error meets the positioning requirement or the temperature TkBelow a set minimum temperature value.
e. And calculating the energy focus value E (V, H, theta) of the perforation position under the condition of the new velocity model by adopting an amplitude superposition method (specifically, as described in c). If the ratio of E (V,H，θ)＜E(V0，H0，θ0) Carrying out simulated annealing cooling treatment; if E (V, H, theta) > E (V)0，H0，θ0) Repositioning the perforation by adopting a grid successive division method, and calculating perforation repositioning errors; the invention takes the energy focusing value E (V, H, theta) at the position of the perforation as the objective function of simulated annealing, and when the energy focusing value E (V, H, theta) > E (V, H, theta) is satisfied0，H0，θ0) Under the condition, the grid successive division algorithm is utilized to perform positioning calculation, so that the calculation efficiency is greatly improved.
The mesh successive subdivision specific steps are as follows:
e1, setting the radius r of the target area to be 400, setting the center P of the target area to be the perforation position, namely P (-280,360, -3550), and dividing the grid into grids with precision d in200, positioning accuracy requirement d out1, positioning result Lout(0,0,0), energy focus maximum Eout＝0；
E2, performing mesh division on the target area, calculating the energy focus value of the central point of each mesh, and recording the maximum energy focus value as EmaxAnd the center point of the grid is Lmax；
E3, if Emax＞EoutThen order Eout＝Emax，Lout＝Lmax；
e4, if din＞doutIf r is equal to r/2, din＝din/2，P＝LoutRepeating the second step; if d isin＜doutThen store LoutAnd the final perforation positioning result.
f. If the positioning error does not meet the positioning requirement (the positioning error is within 5 m), carrying out simulated annealing cooling treatment, if the temperature after cooling is higher than the set lowest temperature, carrying out tempering treatment to adjust the speed model, otherwise, storing the optimal speed model, the perforation repositioning result and the error, and then ending; and if the positioning error meets the positioning requirement, storing the optimal speed model, the perforation repositioning result and the error, and ending.
In the invention, the speed model is complicated from horizontal stratigraphy to inclined stratum, the consideration of stratum interface dip angle parameters is increased, the method is closer to the actual stratum, and a three-dimensional ray tracing algorithm which is more complicated than the horizontal stratigraphy is also needed.
The method is equivalent to the preliminary screening of a new speed model through the judgment condition of the grid successive subdivision, and only the speed model meeting the condition is subjected to perforation relocation, so that the calculation efficiency is greatly improved.
In the invention, a grid successive division method is adopted for positioning, so that the calculation precision and efficiency are greatly improved. For example, in a three-dimensional region with a division radius of 200m, the mesh division size is 1m, and if the mesh division method is adopted, the mesh division needs to be divided into 6.4x107If a mesh successive splitting method is adopted, the initial splitting size is set to be 100m, the splitting times are 8, the final mesh splitting size is 0.78125m < 1m, but 64 x 8 to 512 meshes are only needed to be calculated in total, and the calculation amount is reduced to 10 original meshes-5。
As shown in fig. 5 (the five-pointed star is the real perforation position, and the box is the perforation repositioning result), the perforation repositioning error obtained by adopting the method before improvement (the horizontal laminar velocity model) is 221.7m, and the perforation repositioning error obtained by adopting the method of the invention is within 5m, so that the method of the invention has obvious advantages under the conditions of large depth and complex stratum structure.
TABLE 1
Figure BDA0003289480040000111
It is to be noted that the foregoing is only illustrative of the preferred embodiments of the present invention and the technical principles employed. It will be understood by those skilled in the art that the present invention is not limited to the particular embodiments described herein, but is capable of various obvious changes, rearrangements and substitutions as will now become apparent to those skilled in the art without departing from the scope of the invention. Therefore, although the present invention has been described in greater detail by the above embodiments, the present invention is not limited to the above embodiments, and may include other equivalent embodiments without departing from the spirit of the present invention, and the scope of the present invention is determined by the scope of the appended claims.

Claims (5)

1. A microseism ground monitoring velocity model correction algorithm is characterized by comprising the following steps:
a. arranging N detectors on the ground, and defining a three-dimensional target area by taking the position of a perforation point as a center;
b. establishing an initial velocity model according to the logging curve and reading N seismic wave data acquired by the N detectors;
c. calculating the energy focus value E (V) of the perforation position by adopting an amplitude superposition method0,H0,θ0)；
d. Taking the speed value, the level value and the rotation angle value of the target layer as uncertain factors, and adjusting the parameters of the target layer by adopting a very fast simulated annealing method to obtain a new speed model;
e. recalculating the energy focus value E (V, H, theta) of the perforation position if E (V, H, theta)<E(V0,H0,θ0) Tempering to adjust the speed model; if E (V, H, theta)>E(V0,H0,θ0) Repositioning the perforation by adopting a grid successive division method, and calculating perforation repositioning errors;
f. if the positioning error does not meet the positioning requirement, carrying out simulated annealing cooling treatment, if the temperature after cooling is higher than the set lowest temperature, carrying out tempering treatment to adjust the speed model, otherwise, storing the optimal speed model, the perforation repositioning result and the error, and then ending; and if the positioning error meets the positioning requirement, storing the optimal speed model, the perforation repositioning result and the error, and ending.
2. The microseism ground monitoring velocity model correction algorithm of claim 1, wherein in the step b, the specific steps of establishing the initial velocity model are as follows:
obtaining an initial velocity vector V, V ═ V from the logp1,Vp2,Vp3,…,Vpn]The horizon vector H, H ═ H1,H2,H3,…,Hn-1]Turning roundCorner vector thetax、θy，θx＝[θx1,θx2,θx3,…,θxn-1]，θy＝[θy1,θy2,θy3,…,θyn-1]Wherein V ispiVelocity of P-wave of i-th layer, HiIs the depth coordinate of the ith layer interface, θxi、θyiThe values of the angles of the horizontal plane rotating around the x axis and the y axis to the ith layer interface position are respectively positive anticlockwise and negative clockwise.
3. The microseism ground monitoring velocity model correction algorithm of claim 1, wherein in the step c, the amplitude superposition method comprises the following specific steps:
c1, selecting any one channel in the perforation data as a reference channel M, and solving the theoretical travel time difference of other channels relative to the reference channel by adopting a ray tracing method:
Δt＝[t1-tM,t2-tM,t3-tM,…,tN-tM]
in the formula, t1,t2,t3,…,tNThe theoretical travel time t of each path except the reference path in the perforation dataMThe theoretical travel time of a reference track;
c2, performing reverse time migration superposition on the perforation record waveform obtained by each detector according to the theoretical travel time difference, wherein the mathematical expression of the objective function is as follows:
Figure FDA0003289480030000021
wherein A (i, j) is the amplitude of the ith channel data at the jth moment, N is the number of detectors, and W is the length of a time window.
4. The microseism ground monitoring velocity model correction algorithm of claim 1, wherein in step d, the extremely fast simulated annealing method comprises the following specific steps:
d1, calculating the simulated annealing initial temperature T0The solution to obtain the initial temperature is as follows:
giving a small positive value to the initial temperature, and then continuously multiplying the initial temperature by a number which is constantly larger than 1 until the acceptance probability of any model is close to 1;
Figure FDA0003289480030000022
wherein H0,V0,θ0As a parameter of the initial velocity model, H1,V1,θ1Is the model result of the first iteration;
d2, simulating annealing cooling calculation, wherein the cooling formula is as follows:
Tk＝T0exp(-ck1/2n)
where k is the number of iterations, T0C is an initial temperature, a given constant is given, and a cooling process is determined, wherein c is 0.5, and n is the number of layers of the stratum needing to be adjusted;
d3, in the course of simulated annealing calculation, each iteration needs to adjust the velocity vector V, horizon vector H and rotation angle vector theta in the stratum modelx、θyThe adjustment formula is as follows:
Figure FDA0003289480030000031
Figure FDA0003289480030000032
Figure FDA0003289480030000033
Figure FDA0003289480030000034
wherein,
Figure FDA0003289480030000035
the maximum and minimum values of the speed adjustment range for each layer,
Figure FDA0003289480030000036
Figure FDA0003289480030000037
the maximum and minimum values of the adjustment range are adjusted for each layer position,
Figure FDA0003289480030000038
Figure FDA0003289480030000039
Figure FDA00032894800300000310
adjusting the maximum value and the minimum value of the range for the interfaces of the layers and the rotating angle,
Figure FDA00032894800300000311
Figure FDA00032894800300000312
x1、x2、x3、x4is a random variable with the value range of [ -1,1]Generating x1,2,3,4The expression of (a) is as follows:
Figure FDA00032894800300000313
wherein sgn is a sign function, μ ∈ [0,1 ];
d4, in each iteration, if E (V ', H', theta ') > is equal to or more than E (V, H, theta), then V', H ', theta' replaces V, H, theta as the current optimal solution, if E (V ', H', theta ') is equal to or less than E (V, H, theta), then the probability P (V → V') is accepted to replace V, H, theta 'as the current optimal solution, and the calculation formula of P (V → V') is as follows:
Figure FDA00032894800300000314
wherein V ', H ', theta ' are velocity models in the iterative process, V, H, theta are the current optimal solution of the velocity models, and TkThe temperature value of the kth iteration is the iteration end condition of simulated annealing, and the iteration end condition is E (V, H, theta) > E (V)0,H0,θ0)。
5. The microseism ground monitoring velocity model correction algorithm according to claim 1, wherein in the step e, the grid successive subdivision positioning method comprises the following specific steps:
e1, setting the radius r of the target area, the center P of the target area as the perforation position, and the mesh division precision dinRequirement for positioning accuracy doutPositioning result Lout(0,0,0), energy focus maximum Eout＝0；
E2, performing mesh division on the target area, calculating the energy focus value of the central point of each mesh, and recording the maximum energy focus value as EmaxAnd the center point of the grid is Lmax；
E3, if Emax＞EoutThen order Eout＝Emax，Lout＝Lmax；
e4, if din＞doutIf r is equal to r/2, din＝din/2，P＝LoutRepeating the second step; if d isin＜doutThen store LoutAnd the final perforation positioning result.
CN202111159206.5A 2021-09-30 2021-09-30 Microseism ground monitoring speed model correction method Pending CN113885076A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111159206.5A CN113885076A (en) 2021-09-30 2021-09-30 Microseism ground monitoring speed model correction method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111159206.5A CN113885076A (en) 2021-09-30 2021-09-30 Microseism ground monitoring speed model correction method

Publications (1)

Publication Number Publication Date
CN113885076A true CN113885076A (en) 2022-01-04

Family

ID=79004627

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111159206.5A Pending CN113885076A (en) 2021-09-30 2021-09-30 Microseism ground monitoring speed model correction method

Country Status (1)

Country Link
CN (1) CN113885076A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114460635A (en) * 2022-02-09 2022-05-10 中国矿业大学(北京) Method and device for constructing microseism velocity model and electronic equipment

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105954795A (en) * 2016-04-25 2016-09-21 吉林大学 Grid successive dissection method used for microseismic positioning
CN106814391A (en) * 2015-11-27 2017-06-09 中国石油化工股份有限公司 Ground micro-seismic state event location method based on Fresnel zone tomographic inversion
CN107132578A (en) * 2017-04-06 2017-09-05 吉林大学 A kind of microseism ground monitoring velocity model corrections algorithm
CN108897035A (en) * 2018-05-14 2018-11-27 吉林大学 A kind of microseism weighting localization method based on wave detector weight
WO2019071504A1 (en) * 2017-10-12 2019-04-18 南方科技大学 Two-point ray tracing based seismic travel time tomography inversion method

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106814391A (en) * 2015-11-27 2017-06-09 中国石油化工股份有限公司 Ground micro-seismic state event location method based on Fresnel zone tomographic inversion
CN105954795A (en) * 2016-04-25 2016-09-21 吉林大学 Grid successive dissection method used for microseismic positioning
CN107132578A (en) * 2017-04-06 2017-09-05 吉林大学 A kind of microseism ground monitoring velocity model corrections algorithm
WO2019071504A1 (en) * 2017-10-12 2019-04-18 南方科技大学 Two-point ray tracing based seismic travel time tomography inversion method
CN108897035A (en) * 2018-05-14 2018-11-27 吉林大学 A kind of microseism weighting localization method based on wave detector weight

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114460635A (en) * 2022-02-09 2022-05-10 中国矿业大学(北京) Method and device for constructing microseism velocity model and electronic equipment
CN114460635B (en) * 2022-02-09 2022-07-29 中国矿业大学(北京) Method and device for constructing microseism velocity model and electronic equipment

Similar Documents

Publication Publication Date Title
CN107132578B (en) A kind of microseism ground monitoring velocity model corrections algorithm
CN106353792B (en) Method suitable for positioning micro-seismic source of hydraulic fracturing
CN105093319B (en) Ground micro-seismic static correcting method based on 3D seismic data
CN104570106A (en) Near-surface tomographic velocity analysis method
CN110687602A (en) Shallow seismic multi-wave combined exploration method
CN108072900A (en) A kind of trace gather record processing method, device and computer storage media
CN103105622B (en) Based on the homotype ripple time difference positioning method of database technology
CN110389377B (en) Microseism offset imaging positioning method based on waveform cross-correlation coefficient multiplication
CN113885076A (en) Microseism ground monitoring speed model correction method
CN108897035A (en) A kind of microseism weighting localization method based on wave detector weight
CN113534259A (en) Vibroseis efficient acquisition real-time prestack time migration imaging method
CN110187386B (en) DTW seismic body attribute analysis method for automatically and rapidly identifying geological structure
CN104614762B (en) Loose sandstone gas reservoir boundary determining method and device
WO2024051834A1 (en) Full-waveform inversion method and device, and storage medium
CN105510958A (en) Three-dimensional VSP observation system designing method suitable for complex medium
CN110646840A (en) Angle gather extraction method and system
CN111474580A (en) Azimuth angle gather extraction method and system based on offset vector piece
CN113075636B (en) Parallel line coordinate transformation and weak target detection method for measuring points
CN111650645A (en) Variable offset VSP curved line correction processing method and device
CN111337975A (en) Automatic microseism event identification method based on variance fractal dimension
CN104502972A (en) Three-component earthquake wave integral migration method and device
CN113985493B (en) Intelligent modeling method for underground multi-information constrained isochronous stratum grillwork
CN105549075B (en) Method and device for solving shallow velocity field
CN110488350B (en) Seismic inversion big data generation method based on convolutional neural network
CN112379430B (en) Multi-component offset imaging method in angle domain

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: 20220104