CN113885076A - A Method of Correcting Velocity Models for Microseismic Ground Monitoring - Google Patents

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

A Method of Correcting Velocity Models for Microseismic Ground Monitoring

technical field

The invention belongs to the technical field of microseismic monitoring, and in particular relates to a microseismic ground monitoring velocity model correction algorithm.

Background technique

Microseismic monitoring technology is to pick up the vibration signals generated by fracturing from the deployed microseismic detectors, and provide the spatial shape of the fractures through positioning processing technologies such as superimposed energy scanning analysis, and monitor the scale and shape of the artificial reservoir formed by hydraulic fracturing. and fracture network structure to guide the effective development of fracturing work and improve productivity. Microseismic positioning technology is the core of microseismic monitoring work. During the implementation of this technology, the extension of cracks causes the surrounding rocks to rupture, thereby triggering a series of microseismic events that can be observed and recorded. The velocity model within the monitoring area is the main factor affecting the accuracy of microseismic event location. Therefore, how to obtain an effective velocity model is a key issue in microseismic monitoring engineering.

At present, the correction of the microseismic ground monitoring velocity model is based on the horizontal layered stratigraphic model, and the velocity model is relatively simple, resulting in a certain error in the final perforation positioning results. Chinese patent CN107132578B discloses a microseismic ground monitoring velocity model correction algorithm. On the basis of the construction method of the microseismic velocity model based on the amplitude superposition, it is proposed to improve the accuracy of perforation relocation by expanding the solution space, that is, in the extremely fast simulated annealing method In the process, the uncertainty of the velocity of each layer in the formation and the uncertainty of the horizon of each layer are considered at the same time. However, the above method is based on the horizontal layered stratum model, and the velocity model is relatively simple, which leads to a certain error in the final perforation relocation result, especially when the actual stratum is complex, the perforation relocation error cannot meet the requirements. The velocity model required for the location of subsequent microseismic events cannot be obtained. In addition, the correction algorithm performs perforation relocation every time a new velocity model is obtained by adjusting the parameters, and the positioning process takes a long time. In addition, the positioning method is mesh division, and only one division is performed according to the set division accuracy. If the size of the mesh is small, the positioning accuracy is high, but the number of meshes will be very large, and the calculation time is very long. ; If the grid size is large, the calculation time will be short, but the positioning accuracy will be very low, and high calculation efficiency and positioning accuracy cannot be obtained at the same time.

Therefore, how to make up for the existing defects, improve the perforation relocation accuracy under complex geological conditions, obtain a more equivalent velocity model, and then improve the location accuracy of microseismic events is an urgent problem to be solved in this field.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide a microseismic ground monitoring velocity model correction algorithm to solve the problem of improving the perforation relocation accuracy under complex geological conditions, obtaining a more equivalent velocity model, and then improving the microseismic event positioning accuracy. . On the basis of the horizontal layered model, the velocity model is complicated, and the variation of the velocity value of each layer, the variation of the interface position of each layer, and the variation of the rotation angle of each layer are considered to further expand the solution space and improve the accuracy of perforation relocation.

The purpose of this invention is to realize through the following technical solutions:

A microseismic ground monitoring velocity model correction algorithm, comprising the following steps:

a. Arrange N geophones on the ground, and define a three-dimensional target area with the position of the perforation point as the center;

b. Establish the initial velocity model according to the logging curve and read the N-channel seismic wave data obtained by the N geophones;

c. Calculate the energy focus value E(V ₀ , H ₀ , θ ₀ ) of the perforation position by using the amplitude superposition method;

d. Taking the velocity value, horizon value and rotation angle value of the target layer as uncertain factors, the parameters of the target layer are adjusted by the extremely rapid simulated annealing method to obtain a new velocity model;

e. Recalculate the energy focus value E(V,H,θ) of the perforation position. If E(V,H,θ)<E(V ₀ ,H ₀ ,θ ₀ ), adjust the speed model by tempering; If E(V, H, θ)>E(V ₀ , H ₀ , θ ₀ ), the perforation is relocated by the grid successive division method, and the perforation relocation error is calculated;

f. If the positioning error does not meet the positioning requirements, perform simulated annealing and cooling treatment. If the temperature after cooling is higher than the set minimum temperature, perform tempering treatment to adjust the speed model. Otherwise, save the optimal speed model, perforation repositioning results and It ends after the error; if the positioning error meets the positioning requirements, it saves the optimal velocity model, the perforation relocation results and the error and ends.

Further, in step b, the concrete steps for establishing the initial velocity model are:

Obtain the initial velocity vector V, V=[V _p1 , V _p2 , V _p3 ,..., V _pn ] according to the logging curve, and the initial horizon vector H, H=[H ₁ , H ₂ , H ₃ ,..., H _{n -1} ], initial rotation angle vectors θ _x , θ _y , θ _x = [θ _x1 , θ _x x, θ _x3 , ..., θ _xn-1 ], θ _y = [θ _y1 , θ _y2 , θ _y3 , ... , θ _yn-1 ], where V _pi is the velocity of the P wave in the i-th layer, H _i is the depth coordinate of the i-th layer interface, θ _xi , θ _yi are the rotation of the horizontal plane around the x-axis and the y-axis to the i-th layer, respectively The angle value of the interface position of each layer, counterclockwise is positive, clockwise is negative.

Further, step c, the described amplitude superposition method, the specific steps are as follows:

c1, in the perforation data, select any one as the reference track M, and use the ray tracing method to obtain the theoretical travel time difference of the other tracks relative to the reference track:

Δt=[t ₁ -t _M , t ₂ -t _M , t ₃ -t _M , . . . , t _N -t _M ]

In the formula, t ₁ , t ₂ , t ₃ , ..., t _N is the theoretical travel time of each track except the reference track in the perforation data, and t _M is the theoretical travel time of the reference track;

c2, according to the theoretical travel time difference, perform reverse time migration and superposition of the perforation recording waveforms obtained by each detector. The mathematical expression of the objective function is as follows:

Among them, A(i, j) is the amplitude of the i-th channel data at the j-th time, N is the number of detectors, and W is the length of the time window.

Further, in step d, the specific steps of the extremely rapid simulated annealing method are as follows:

d1, calculate the initial temperature T ₀ of simulated annealing, and the solution to obtain the initial temperature is as follows:

First assign a small positive value to the initial temperature, and then keep multiplying it by a number that is always greater than 1 until the acceptance probability for any model is close to 1;

Among them, V ₀ , H ₀ , θ ₀ are the parameters of the initial velocity model, and V ₁ , H ₁ , θ ₁ are the model results of the first iteration;

d2, simulated annealing cooling calculation, the cooling formula used is as follows:

T _k =T ₀ exp(-ck ^1/2n )

Among them, k is the number of iterations, T ₀ is the initial temperature, c is a given constant, which determines the cooling process, where c=0.5, and n is the number of layers of the formation to be adjusted;

d3, in the process of simulated annealing calculation, the velocity vector V, horizon vector H and rotation angle vectors θ _x and θ _y in the formation model also need to be adjusted in each iteration. The adjustment formula is as follows:

为各层速度调整范围的的最大值和最小值，

为各层位置调整范围的的最大值和最小值，

为各层分界面与旋转角度调整范围的的最大值和最小值，

x₁、x₂、x₃、x₄为随机变量，取值范围为[-1，1]，产生x_{1，2，3，4}的表达式如下：in,

Adjust the maximum and minimum values of the range for the speed of each layer,

Adjust the maximum and minimum values of the range for each layer position,

The maximum and minimum values of the adjustment range for the interface and rotation angle of each layer,

x ₁ , x ₂ , x ₃ , and x ₄ are random variables with a value range of [-1, 1]. The expressions for generating x ₁ , 2, 3, and 4 are as follows:

where sgn is the sign function, μ∈[0, 1];

d4, in each iteration, if E(V', H', θ')>E(V, H, θ), then V', H', θ' replace V, H, θ as the current optimal solution, If E(V', H', θ')<E(V, H, θ), then with probability P(V, H, θ→V', H', θ'), accept V', H', θ′ replaces V, H, and θ as the current optimal solution. The calculation formula of P(V, H, θ→V′, H′, θ′) is as follows:

Among them, V′, H′, θ′ are the velocity model in the iterative process, V, H, θ are the current optimal solution of the velocity model, Tk is the temperature value at the _k -th iteration, the iteration termination condition of simulated annealing For the positioning error to meet the positioning requirements or the temperature _Tk is lower than the set minimum temperature value.

Further, in step e, the specific steps of the grid successive subdivision positioning method are as follows:

e1, set the radius r of the target area, the center P of the target area is the perforation position, the grid division accuracy d _in , the positioning accuracy requirement d _out , the positioning result L _out =(0, 0, 0), the maximum energy focus E _out = 0;

e2, carry out grid division to the target area, calculate the energy focus value of each grid center point, and record the maximum energy focus value as E _max and the grid center point as L _max ;

e3, if E _max >E _out , then let E _out =E _max , L _out =L _max ;

e4, if d _in >d _out , then set r=r/2, d _in =d _in /2, P=L _out , repeat the second step; if d _in <d _out , save L _out as the final perforation Positioning results.

Compared with the prior art, the beneficial effects of the present invention are: the present invention complicates the velocity model on the basis of the horizontal layered model, and simultaneously considers the variation of the layer interface position, the variation of the rotation angle of each layer, and the variation of the velocity value of each layer, In order to further expand the solution space and improve the accuracy of perforation relocation. The experimental results show that, compared with the existing method, the method can relocate the perforation event to its true value more accurately, and can effectively improve the positioning reliability of the microseismic event near the perforation point.

Description of drawings

In order to illustrate the technical solutions of the embodiments of the present invention more clearly, the following briefly introduces the accompanying drawings used in the embodiments. It should be understood that the following drawings only show some embodiments of the present invention, and therefore do not It should be regarded as a limitation of the scope, and for those of ordinary skill in the art, other related drawings can also be obtained according to these drawings without any creative effort.

Fig. 1 is a flow chart of the method of the present invention;

Figure 2a is a three-dimensional layered undulation structure stratigraphic model diagram;

Fig. 2b Stratigraphic interface of the fifth layer;

Figure 3 forward modeling waveform diagram of each channel of the detector;

Figure 4a ray tracing forward modeling 3D view;

Figure 4b ray tracing forward modeling vertical section view;

Fig. 5a-Fig. 5b perforation relocation results of the method before the improvement;

Figures 5c-5d are the results of perforation relocation according to the method of the present invention.

Detailed ways

The present invention will be further described below in conjunction with specific implementation cases. The implementation cases are implemented on the premise of the technology of the present invention, and provide detailed implementation modes, but the protection scope of the present invention is not limited to the following examples.

As shown in Figure 1, the microseismic ground monitoring velocity model correction algorithm of the present invention includes the following steps:

a. In this example, a 9-layer layered undulating formation model is selected (as shown in Figure 2a), and 8 curved surfaces are used as the horizon interface (the surfaces of each interface are not parallel to each other, but the undulation forms are shown in Figure 2b) . In the local coordinate system, the designed perforation position is (-280, 360, -3550) (unit: m). The microseismic ground detectors are arranged in a star shape, with 6 survey lines in total, each survey line has 10 detectors (60 channels), and the detector spacing is 50m (as shown in Figure 2a, the blue cross is the detector position, The red five-pointed star is the perforation location). The actual signal is simulated by the finite difference wave equation, and the waveform is described by the 20Hz rake wavelet, and the maximum amplitude of each rake wavelet is 1 (as shown in Figure 3). The subsequent forward modeling of the velocity model correction scheme adopts the ray tracing method (as shown in Fig. 4).

b. Obtain the initial velocity model according to the logging information (shown by the black solid line in Fig. 2a), and the initial velocity vector is 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 speed model calibration, the adjustment range of the speed of each layer is set to 600m/s, the adjustment range of the layer value of each layer is set to 100m, the adjustment range of the rotation angle of each layer is set to 6°, and the parameters of each layer are adjusted in detail. The range is shown in Table 1.

c. Read the analog signals of the 60-channel detectors, and use the amplitude superposition method to calculate the energy focus value E(V ₀ , H ₀ , θ ₀ )=713.882 at the perforation position under the initial velocity model condition. The specific steps of the amplitude superposition method are as follows :

c1, in the perforation data, select any one as the reference track M, and use the ray tracing method to obtain the theoretical travel time difference of the other tracks relative to the reference track:

Δt=[t ₁ -t _M , t ₂ -t _M , t ₃ -t _M , . . . , t _N -t _M ]

In the formula, t ₁ , t ₂ , t ₃ , ..., t _N is the theoretical travel time of each track except the reference track in the perforation data, and t _M is the theoretical travel time of the reference track;

c2, according to the theoretical travel time difference, perform reverse time migration and superposition of the perforation recording waveforms obtained by each detector. The mathematical expression of the objective function is as follows:

where A(i, j) is the amplitude of the i-th channel data at the j-th moment, N is the number of detectors, and W is the time window length;

d. Use the extremely fast simulated annealing algorithm to adjust the speed value, horizon value and rotation angle value of the target layer to obtain a new speed model. The equivalent stratigraphic model is complicated into a sloping stratigraphic model that is closer to the actual stratigraphy, and the velocity value, horizon value and interface rotation angle of each layer can be independently adjusted.

d1, calculate the initial temperature T ₀ of simulated annealing, and the solution to obtain the initial temperature is as follows:

First assign a small positive value to the initial temperature, and then keep multiplying it by a number that is always greater than 1 until the acceptance probability for any model is close to 1;

Among them, V ₀ , H ₀ , θ ₀ are the parameters of the initial velocity model, and V ₁ , H ₁ , θ ₁ are the model results of the first iteration; the approximate initial temperature obtained in this example is 200°C.

d2, simulated annealing cooling calculation, the cooling formula used is as follows:

T _k =T ₀ exp(-ck ^1/2n )

Among them, k is the number of iterations, T ₀ is the initial temperature, c is a given constant, which determines the cooling process, and n is the number of layers to be adjusted. In this example, k=500, c=0.5, n= 9;

d3, in the process of simulated annealing calculation, the velocity vector V, horizon vector H and rotation angle vectors θ _x and θ _y in the formation model also need to be adjusted in each iteration. The adjustment formula is as follows:

为各层速度调整范围的的最大值和最小值，

为各层位置调整范围的的最大值和最小值，

为各层分界面与旋转角度调整范围的的最大值和最小值，

x₁、x₂、x₃、x₄为随机变量，取值范围为[-1，1]，产生x_{1，2，3，4}的表达式如下：in,

Adjust the maximum and minimum values of the range for the speed of each layer,

Adjust the maximum and minimum values of the range for each layer position,

The maximum and minimum values of the adjustment range for the interface and rotation angle of each layer,

x ₁ , x ₂ , x ₃ , and x ₄ are random variables with a value range of [-1, 1]. The expressions for generating x ₁ , 2, 3, and 4 are as follows:

where sgn is the sign function, μ∈[0, 1];

d4, in each iteration, if E(V', H', θ')>E(V, H, θ), then V', H', θ' replace V, H, θ as the current optimal solution, If E(V', H', θ')<E(V, H, θ), then with probability P(V, H, θ→V', H', θ'), accept V', H', θ′ replaces V, H, and θ as the current optimal solution. The calculation formula of P(V, H, θ→V′, H′, θ′) is as follows:

Among them, V′, H′, θ′ are the velocity model in the iterative process, V, H, θ are the current optimal solution of the velocity model, Tk is the temperature value at the _k -th iteration, the iteration termination condition of simulated annealing For the positioning error to meet the positioning requirements or the temperature _Tk is lower than the set minimum temperature value.

e. Use the amplitude superposition method (specifically as described in c) to calculate the energy focus value E(V, H, θ) of the perforation position under the condition of the new velocity model. If E(V, H, θ)<E(V ₀ , H ₀ , θ ₀ ), perform simulated annealing cooling treatment; if E(V, H, θ)>E(V ₀ , H ₀ , θ ₀ ), Then, the grid successively dividing method is used to relocate the perforation, and the perforation relocation error is calculated; the present invention takes the energy focus value E (V, H, θ) at the perforation as the objective function of simulated annealing, and when it satisfies Only when E(V, H, θ)>E(V ₀ , H ₀ , θ ₀ ) is used, the grid successive division algorithm is used for positioning calculation, which greatly improves the calculation efficiency.

The specific steps of meshing successively are as follows:

e1, set the target area radius r=400, the center P of the target area is the perforation position, that is, P=(-280, 360, -3550), the grid division accuracy d _in = 200, the positioning accuracy requirement d _out = 1, the positioning result L _out =(0, 0, 0), the maximum energy focus E _out =0;

e2, carry out grid division to the target area, calculate the energy focus value of each grid center point, and record the maximum energy focus value as E _max and the grid center point as L _max ;

e3, if E _max >E _out , then let E _out =E _max , L _out =L _max ;

e4, if d _in >d _out , then set r=r/2, d _in =d _in /2, P=L _out , repeat the second step; if d _in <d _out , save L _out as the final perforation Positioning results.

f. If the positioning error does not meet the positioning requirements (the positioning error is within 5m), perform simulated annealing and cooling treatment. If the temperature after cooling is higher than the set minimum temperature, perform tempering treatment to adjust the speed model, otherwise, save the optimal speed model. , end after perforating relocation results and errors; if the positioning error meets the positioning requirements, save the optimal velocity model, perforation relocation results and errors and end.

In the present invention, the velocity model is complicated from the horizontal layer to the inclined layer, which increases the consideration of the dip angle parameter of the layer interface, and is closer to the actual layer, which also requires a more complex three-dimensional ray tracing algorithm than the horizontal layer.

In the present invention, the judgment condition of successive grid division is equivalent to preliminarily screening the new velocity model, and only the velocity model satisfying the condition is perforated and relocated, which greatly improves the calculation efficiency.

The positioning in the present invention adopts the grid successive subdivision method, which greatly improves the calculation accuracy and efficiency. For example, in a three-dimensional area with a division radius of 200m, the mesh size is 1m. If the mesh division method is used, the mesh division needs to be divided into 6.4x10 ⁷ meshes. If the mesh successive division method is used, The initial subdivision size is set to 100m and the subdivision times is 8, then the final mesh subdivision size is 0.78125m<1m, but only 64*8=512 grids need to be calculated in total, and the calculation amount is reduced to the original one. ^10-5 .

As shown in Figure 5 (the five-pointed star is the real position of the perforation, and the box is the result of the perforation relocation), the perforation relocation error obtained by the method before the improvement (horizontal layered velocity model) is 221.7m, and the method of the present invention is adopted. The perforation relocation error is within 5m, and it can be seen that the method of the present invention has obvious advantages under the conditions of relatively large depth and relatively complex formation structure.

Table 1

Note that the above are only preferred embodiments of the present invention and applied technical principles. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and various obvious changes, readjustments and substitutions can be made by those skilled in the art without departing from the protection scope of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments, and can also include more other equivalent embodiments without departing from the concept of the present invention. The scope is determined by the scope of the appended claims.

Claims (5)

1. a microseismic ground monitoring velocity model correction algorithm, is characterized in that, comprises the following steps: a. Arrange N geophones on the ground, and define a three-dimensional target area with the position of the perforation point as the center; b. Establish the initial velocity model according to the logging curve and read the N-channel seismic wave data obtained by the N geophones; c. Calculate the energy focus value E(V ₀ , H ₀ , θ ₀ ) of the perforation position by using the amplitude superposition method; d. Taking the velocity value, horizon value and rotation angle value of the target layer as uncertain factors, the parameters of the target layer are adjusted by the extremely rapid simulated annealing method to obtain a new velocity model; e. Recalculate the energy focus value E(V,H,θ) of the perforation position. If E(V,H,θ)<E(V ₀ ,H ₀ ,θ ₀ ), adjust the speed model by tempering; If E(V, H, θ)>E(V ₀ , H ₀ , θ ₀ ), the perforation is relocated by the grid successive division method, and the perforation relocation error is calculated; f. If the positioning error does not meet the positioning requirements, perform simulated annealing and cooling treatment. If the temperature after cooling is higher than the set minimum temperature, perform tempering treatment to adjust the speed model. Otherwise, save the optimal speed model, perforation repositioning results and It ends after the error; if the positioning error meets the positioning requirements, it saves the optimal velocity model, the perforation relocation results and the error and ends. 2. a kind of microseismic ground monitoring velocity model correction algorithm according to claim 1, is characterized in that, in step b, the described concrete steps of establishing initial velocity model are: Obtain the initial velocity vector V according to the logging curve, V=[V _p1 , V _p2 , V _p3 ,...,V _pn ], the horizon vector H, H=[H ₁ , H ₂ , H ₃ ,..., H _{n- 1} ], rotation angle vectors θ _x , θ _y , θ _x = [θ _x1 , θ _x2 , θ _x3 ,…, θ _xn-1 ], θ _y = [θ _y1 , θ _y2 , θ _y3 ,…, θ _{yn -1} ], where V _pi is the velocity of the P wave in the i-th layer, H _i is the depth coordinate of the i-th layer interface, θ _xi and θ _yi are the rotation of the horizontal plane around the x-axis and the y-axis to the i-th layer interface, respectively The angle value of the position, counterclockwise is positive, clockwise is negative. 3. a kind of microseismic ground monitoring velocity model correction algorithm according to claim 1, is characterized in that, step c, described amplitude superposition method, concrete steps are as follows: c1. In the perforation data, select any track as the reference track M, and use the ray tracing method to obtain the theoretical travel time difference of the other tracks relative to the reference track: Δt=[t ₁ -t _M , t ₂ -t _M , t ₃ -t _M ,...,t _N -t _M ] In the formula, t ₁ , t ₂ , t ₃ ,..., t _N is the theoretical travel time of each track except the reference track in the perforation data, and t _M is the theoretical travel time of the reference track; c2. According to the theoretical travel time difference, perform reverse time migration and superposition of the perforation recording waveforms obtained by each detector. The mathematical expression of the objective function is as follows:

Among them, A(i, j) is the amplitude of the i-th channel data at the j-th time, N is the number of detectors, and W is the time window length. 4. a kind of microseismic ground monitoring velocity model correction algorithm according to claim 1, is characterized in that, in step d, the concrete steps of described extremely fast simulated annealing method are as follows: d1, calculate the initial temperature T ₀ of simulated annealing, and the solution to obtain the initial temperature is as follows: First assign a small positive value to the initial temperature, and then keep multiplying it by a number that is always greater than 1 until the acceptance probability for any model is close to 1;

Among them, H ₀ , V ₀ , θ ₀ are the parameters of the initial velocity model, and H ₁ , V ₁ , θ ₁ are the model results of the first iteration; d2, simulated annealing cooling calculation, the cooling formula used is as follows: T _k =T ₀ exp(-ck ^1/2n ) Among them, k is the number of iterations, T ₀ is the initial temperature, c is a given constant, which determines the cooling process, where c=0.5, and n is the number of layers of the formation to be adjusted; d3, in the process of simulated annealing calculation, the velocity vector V, horizon vector H and rotation angle vectors θ _x and θ _y in the formation model also need to be adjusted in each iteration. The adjustment formula is as follows:

in,

Adjust the maximum and minimum values of the range for the speed of each layer,

Adjust the maximum and minimum values of the range for each layer position,

Adjust the maximum and minimum values for the interface and rotation angle of each layer,

x ₁ , x ₂ , x ₃ , and x ₄ are random variables whose value range is [-1, 1]. The expression to generate x ₁ , 2, 3, and 4 is as follows:

where sgn is the sign function, μ∈[0,1]; d4, in each iteration, if E(V',H',θ')≥E(V,H,θ), then V',H',θ' replace V,H,θ as the current optimal solution, If E(V',H',θ')≤E(V,H,θ), then accept V',H',θ' with probability P(V→V') instead of V,H,θ as the current maximum The optimal solution, the calculation formula of P(V→V') is as follows:

Among them, V', H', θ' is the velocity model in the iterative process, V, H, θ is the current optimal solution of the velocity model, T _k is the temperature value at the k-th iteration, the iteration termination condition of simulated annealing E(V, H, θ)>E(V ₀ , H ₀ , θ ₀ ). 5. a kind of microseismic ground monitoring velocity model correction algorithm according to claim 1, is characterized in that, in step e, the concrete steps of grid successively subdividing positioning method are as follows: e1, set the radius r of the target area, the center P of the target area is the perforation position, the grid division accuracy d _in , the positioning accuracy requirement d _out , the positioning result L _out =(0,0,0), the maximum energy focus E _out = 0; e2, carry out grid division to the target area, calculate the energy focus value of each grid center point, and record the maximum energy focus value as E _max and the grid center point as L _max ; e3, if E _max >E _out , then let E _out =E _max , L _out =L _max ; e4, if d _in >d _out , then set r=r/2, d _in =d _in /2, P=L _out , repeat the second step; if d _in <d _out , save L _out as the final perforation Positioning results.

Images