CN105807316A - Surface observation microseism speed model correcting method based on amplitude stack - Google Patents

Abstract

The invention relates to a surface observation microseism speed model correcting method based on amplitude stack.Starting from the principle of a reverse time migration amplitude stack method, a very fast simulated annealing method is combined, seismic phase first arrival information does not need to be picked up, and a speed model is corrected by monitoring the energy focusing condition on a perforation point; perforation point relocation accuracy is used for judging whether the speed model can be used for follow-up microseism location.The problem that the maximal value and the maximum value are not separated in an existing algorithm can be effectively solved, the maximum value E of energy focusing on the perforation point can be found accurately, the perforation position is located accurately, and the speed module is effectively corrected.Noise interference can be effectively inhibited by adopting the reverse time migration amplitude stack algorithm, and a good location effect can still be obtained with a low signal-to-noise ratio.The more the speed model is close to an actual model, the more accurate perforation location becomes.Whether the speed model can be used for follow-up microseism location or not is judged according to perforation point relocation accuracy.

Description

Ground observation microseism velocity model corrections method based on amplitude superposition

Technical field

The present invention relates to a kind of oil field compression fracture microseism location, the monitoring materials being through microseism specifically carries out velocity model corrections, thus improving follow-up microseism positioning precision.

Technical background

Have become as the new focus of domestic and international oil and gas industry currently for the exploration of Low permeable oil and gas reservoirs and exploitation, and become a kind of popular tendency by hydraulic fracturing technology exploitation Low permeable oil and gas reservoirs.The article of " Improvedmicroseismicfracturemappingusingperforationtimin gmeasurementsforvelocitycalibration " by name that Warpinski in 2005 delivers on SPEJournal is pointed out in this technology implementation process, fracture extension causes that surrounding rock breaks, thus causing the micro-seismic event of a series of Observable record.The article of " Petroleumreservoircharacterizationusingdownholemicroseis micmonitoring " by name that the article of " Microseismicity-constrainedfracturemodelsforreservoirsim ulation " by name that Eisner in 2010 etc. deliver on TheLeadingEdge and Maxwell etc. deliver on Geophysics is pointed out by being accurately positioned micro-seismic event, may determine that fracture strike, evaluate fracturing effect and analyze inverting focal mechanism etc..Therefore, what the demand of raising micro-seismic event positioning precision became is more urgent.By perforation state event location to its actual value place, and setting up an effective rate pattern is the key reaching above-mentioned purpose.2010, the article of " Constrainedtomographyofrealisticvelocitymodelsinmicrosei smicmonitoringusingcalibrationshots " by name that Bardainne etc. deliver on GeophysicalProspecting points out that in seismic prospecting, the Technology of Seismic Tomography can be used for reference as microseism locating speed model correcting algorithm, but it is less in seismic source information amount, when geoceiver lazy weight and coverage are less, the Technology of Seismic Tomography is utilized to be difficult to obtain comparatively fine work area rate pattern.

Current existing micro-seismic monitoring velocity model corrections method is mostly set up fairly simple structure model of soil layer and is descended velocity structure to be described over the ground, and known perforating site, and perforating site is carried out inverting location, to reduce location risk.But Bardainne and Gaucher points out that this type of method needs to pick up P ripple or S ripple first break information from earthquake record, therefore it is required that have higher signal to noise ratio in earthquake record, such method is commonly used in well to observe data process.But for surface array formula observe for, particularly in reservoir deeper time, perforation record has low signal-to-noise ratio feature, utilizes the result that existing microseism velocity model corrections method obtains unsatisfactory.And Practical Project also usually there will be the situation of well-log information excalation, it is necessary to adopt the mathematical methods such as interpolation to fill up, affect the accuracy of rate pattern to a certain extent.

Inverse time amplitude excursion stacking method does not need first break picking information, by geological data is translated superposition, obtains Voice segment maximum of points and carries out perforation location, is that current microseism positioning field applies one of more method.But, when finding Voice segment maximum of points, it being commonly present maximum and situation that maximum is difficult to differentiate between, this will directly affect perforation positioning precision, cause velocity model corrections misalignment.

Therefore, for, in actual speed model trimming process, how improving perforation positioning precision, overcoming low signal-to-noise ratio, increase depth of reservoirs, it is thus achieved that rate pattern more accurately, thus improving micro-seismic event positioning precision is the problem that this area needs solution badly.

Summary of the invention

Present invention aims to above-mentioned problems of the prior art, based on the reverse-time migration amplitude of vibration addition method, and in conjunction with simulated annealing method, it is provided that a kind of ground observation microseism velocity model corrections method based on amplitude superposition.

Idea of the invention is that with the perforation position error that log data deviation and Gauusian noise jammer cause for principal element, using velocity model corrections as a kind of error compensation to improve perforation reorientation precision.Perforation positioning precision is more high, and gained rate pattern is more accurate.

Based on the ground observation microseism velocity model corrections method of amplitude superposition, comprise the following steps:

A, centered by shooting point, set up three-dimensional target region, select reference channel M；

B, arranging initial temperature T0, minimum temperature Tmin, dwell time T, simulated annealing calculates parameter；

D, set up initial velocity model, definition E (V)=0；

E, obtain initial velocity vector V, read perforation data；

F, calculate each road travel-time difference relative to library track；

G, by each track data translate superposition, obtain the Energy maximum value E under existing rate pattern and coordinate thereof；

H, the value of E is assigned to E (V)；

The coordinate of i, E (V) reaches T dead time；

J, energy are that the coordinate at (V) place is as the perforation elements of a fix；

K, whether meet perforation positioning precision, be；

L, end.

Simulated annealing described in step b calculates parameter setting, is the feature according to ground array microseism observation, arranges six simulated annealings and calculates parameter:

First parameter is velocity vector V, V=[V_p1,V_p2,V_p3,...,V_pn],V_pi

It is i-th layer of p wave interval velocity, just obtains initial velocity vector V according to well-log information；

Second parameter is object function E (V), in log data, selects, arbitrarily together as library track M, to adopt the method for ray tracing to ask for other each road theoretical travel-time difference relative to this library track:

Δt^cal=[t₁-t_M,t₂-t_M,...,t_N-t_M](1)

According to this theory travel-time difference, the perforation wave recording that each road cymoscope obtains being carried out reverse-time migration superposition, stack power (object function) mathematic(al) representation is as follows:

$E\left(V\right)=\underset{j=1}{\overset{N}{\Σ}}\[\underset{i=1}{\overset{M}{\Σ}}A(i,j)\]---\left(2\right)$

Wherein A is the i-th volume elements amplitude of vibration size in the jth moment, and the center point coordinate of i-th volume elements is (x_i,y_i,z_i), M is cymoscope number, and N is time window length；

3rd parameter is simulated annealing initial temperature T₀, the solution obtaining initial temperature is as follows, first give initial temperature one only small on the occasion of, be then constantly multiplied by the perseverance number more than 1, until the satisfied acceptance probability to any model is close to 1；

4th parameter is annealing cooling parameter, extremely fast simulated annealing cooling profiles such as following formula:

T_k=T₀exp(-ck^1/2N)(3)

Wherein k is iterations, T₀For initial annealing temperature, c is used for adjustment algorithm annealing temperature for given constant, allows c=0.5 here；N regulates the speed the number of plies for needs；

5th parameter is stochastic variable x, and be used for regulating the speed vector, specifically gives formula as follows:

$V_i^{k+1}=V_i^k+x*(V_i^{\max }-V_i^{\min })---\left(4\right)$

WithIt is the minimax boundary value of the i-th interval velocity, whereinX is stochastic variable, produces X expression formula as follows:

$x=\mathrm{sgn}(\mathrm{\μ}-0.5)T_k\[{(1+\frac{1}{T_k})}^{|2\mathrm{\μ}-1|}-1\]---\left(5\right)$

Wherein sgn is sign function, and x span is between [-1,1]；

6th parameter is acceptance probability, and as E (V ') >=E (V), V ' substitutes V as currently most solution, as E (V ') < E (V), with probability

$P(V\&RightArrow;V^{\&prime;})=\exp \&lsqb;\frac{E\left(V\right)-E\left(V^{\&prime;}\right)}{T_k}\&rsqb;---\left(6\right)$

Substitute currently most solution, wherein T_kTemperature value during for kth time iteration.

Simulated annealing calculates parameter iteration end condition:

A, Simulated annealing are reduced to lowest set temperature T_min；

B, shooting point positioning precision reach setting value；

C, when shooting point place energy value reaches certain value, through long period T iterative computation still without being substituted, by adjust above-mentioned parameter, to reach the balance of computational accuracy and computational efficiency.

Beneficial effect: in existing velocity model corrections method, major part method is in accordance with perforation positioning precision to judge that whether rate pattern is accurate.The more high explanation rate pattern of perforation positioning precision is closer to the geographical rate pattern of reality, and rate pattern can position, for follow-up micro-seismic event, the foundation providing good accurately.The present invention adopts simulated annealing on the basis of inverse time amplitude excursion superposition algorithm, can effectively overcome maximum that other algorithms existing exist and maximum regardless of situation, find Voice segment maximum E accurately, be accurately positioned perforating site, effective correction rate model.Adopt inverse time amplitude excursion superposition algorithm can effectively suppress noise jamming, still can obtain better locating effect when signal to noise ratio is relatively low.The present invention need not pick up seismic phase first break information, by monitoring shooting point place Voice segment situation, carries out velocity model corrections.Because rate pattern is closer to actual geographic model, its perforation location will be more accurate.Can according to shooting point reorientation precision, it is determined that whether rate pattern can be used for follow-up microseism location.

Accompanying drawing illustrates:

Fig. 1 is based on the ground observation microseism velocity model corrections method flow diagram of amplitude superposition；

Fig. 2 cymoscope is structured the formation and shooting point forward simulation figure；

Fig. 3 (a) 96 road waveform synthetic data result figure；

Fig. 3 (b) utilizes initial velocity model shooting point place ripple stack result figure；

Shooting point place addition of waveforms result figure after Fig. 3 (c) algorithm correction herein；

Fig. 4 (a), (b) is for utilizing initial velocity model to perforation point location result figure；

Fig. 4 (c), perforation is positioned result figure by (d) after the inventive method corrects；

96 track data synthesis result figure after Fig. 5 superposition Gaussian noise；

Fig. 6 shooting point place addition of waveforms Comparative result figure；

Fig. 7 (a) is at the positioning result figure of X-Y plane；

Fig. 7 (b) is at the positioning result figure of X-Z plane.

Detailed description of the invention:

Below in conjunction with drawings and Examples, the present invention is described in further detail:

Based on the ground observation microseism velocity model corrections method of amplitude superposition, comprise the following steps:

A, centered by shooting point, set up three-dimensional target region, select reference channel M；

B, arranging initial temperature T0, minimum temperature Tmin, dwell time T, simulated annealing calculates parameter；

D, set up initial velocity model, definition E (V)=0；

E, obtain initial velocity vector V, read perforation data；

F, calculate each road travel-time difference relative to library track；

G, by each track data translate superposition, obtain the Energy maximum value E under existing rate pattern and coordinate thereof；

H, the value of E is assigned to E (V)；

The coordinate of i, E (V) reaches T dead time；

J, energy are that the coordinate at (V) place is as the perforation elements of a fix；

K, whether meet perforation positioning precision, be；

L, end.

Simulated annealing described in step b calculates parameter setting, is the feature according to ground array microseism observation, arranges six simulated annealings and calculates parameter:

First parameter is velocity vector V, V=[V_p1,V_p2,V_p3,...,V_pn],V_pi

It is i-th layer of p wave interval velocity, just obtains initial velocity vector V according to well-log information；

Second parameter is object function E (V), in log data, selects, arbitrarily together as library track M, to adopt the method for ray tracing to ask for other each road theoretical travel-time difference relative to this library track:

Δt^cal=[t₁-t_M,t₂-t_M,...,t_N-t_M](1)

According to this theory travel-time difference, the perforation wave recording that each road cymoscope obtains being carried out reverse-time migration superposition, stack power (object function) mathematic(al) representation is as follows:

$E\left(V\right)=\underset{j=1}{\overset{N}{\&Sigma;}}\&lsqb;\underset{i=1}{\overset{M}{\&Sigma;}}A(i,j)\&rsqb;---\left(2\right)$

Wherein A is the i-th volume elements amplitude of vibration size in the jth moment, and the center point coordinate of i-th volume elements is (x_i,y_i,z_i), M is cymoscope number, and N is time window length；

3rd parameter is simulated annealing initial temperature T₀, the solution obtaining initial temperature is as follows, first give initial temperature one only small on the occasion of, be then constantly multiplied by the perseverance number more than 1, until the satisfied acceptance probability to any model is close to 1；

4th parameter is annealing cooling parameter, extremely fast simulated annealing cooling profiles such as following formula:

T_k=T₀exp(-ck^1/2N)(3)

Wherein k is iterations, T₀For initial annealing temperature, c is used for adjustment algorithm annealing temperature for given constant, allows c=0.5 here；N regulates the speed the number of plies for needs；

5th parameter is stochastic variable x, and be used for regulating the speed vector, specifically gives formula as follows:

$V_i^{k+1}=V_i^k+x*(V_i^{\max }-V_i^{\min })---\left(4\right)$

WithIt is the minimax boundary value of the i-th interval velocity, whereinX is stochastic variable, produces X expression formula as follows:

$x=\mathrm{sgn}(\mathrm{\&mu;}-0.5)T_k\&lsqb;{(1+\frac{1}{T_k})}^{|2\mathrm{\&mu;}-1|}-1\&rsqb;---\left(5\right)$

Wherein sgn is sign function, and x span is between [-1,1]；

6th parameter is acceptance probability, and as E (V ') >=E (V), V ' substitutes V as currently most solution, as E (V ') < E (V), with probability

$P(V\&RightArrow;V^{\&prime;})=\exp \&lsqb;\frac{E\left(V\right)-E\left(V^{\&prime;}\right)}{T_k}\&rsqb;---\left(6\right)$

Substitute currently most solution, wherein T_kTemperature value during for kth time iteration.

Simulated annealing calculates parameter iteration end condition:

A, Simulated annealing are reduced to lowest set temperature T_min；

B, shooting point positioning precision reach setting value；

C, when shooting point place energy value reaches certain value, through long period T iterative computation still without being substituted, by adjust above-mentioned parameter, to reach the balance of computational accuracy and computational efficiency.

Embodiment 1:

Below in conjunction with accompanying drawing and simulated formation model, the present invention is carried out clear, complete description.It is embodied as step as follows:

A is as in figure 2 it is shown, stratigraphic model and perforation are just drilling process schematic, and ground acquisition station becomes star shapes arrangement, totally 6 surveys line, 16 cymoscopes of every survey line, amounts to 96 roads.First arrival time difference observation in model calculating is calculated acquisition by stratigraphic model actual value is just being drilled ray tracing, adopts the Ricker wavelet of 10Hz to describe waveform, and each road Ricker wavelet amplitude maximum is 1.The three-dimensional coordinate at definition shooting point place is X_s=236m, Y_s=-158m, Z_s=-947m, and centered by shooting point, define a sufficiently large three-dimensional target region.

B, according to location need to be divided into target area the volume elements that many sizes are identical.Target area is divided into volume elements that the length of side is 10m in this experiment, and each volume elements is considered as the potential site point that a micro-seismic event occurs.Initial temperature T is set₀, minimum temperature T_minWith maximum holding time T_s。

C, base area layer model simulate 96 road Wave datas, and Data Synthesis result is such as shown in Fig. 3 (a), and provides initial velocity model, and initial velocity vector is V=[800,1000,1700,2200,2700], and wherein the velocity amplitude of 1～5 layer is followed successively by；800m/s,1000m/s,1700m/s,2200m/s,2700m/s.

96 track datas select, arbitrarily together as library track, adopt the method for ray tracing according to formula (1)

Δt^cal=[t₁-t_M,t₂-t_M,...,t_N-t_M](1)

Ask for other each road theoretical travel-time difference Δ t relative to this library track^cal, all the other 95 track data times will all translate t under this rate pattern₁-t_M, row amplitude superposition of going forward side by side.Travel through all volume elements, it is thus achieved that this region self-energy focuses on maximum E and the coordinate of Voice segment maximum of points.

D, adopt extremely fast simulated annealing that destination layer is carried out speed adjustment.Initial annealing temperature T is set₀, formula (3) calculate temperature T after being annealed_k, regulate the speed model according to formula (4), (5), under new rate pattern, carry out shooting point reorientation:

By formula: T_k=T₀exp(-ck^1/2N)(3)

Calculating annealed after temperature T_k, arranging c=0.5, k is iterations, be 1, N is first the geologic structure number of plies, N=5.According to formula:

$x=\mathrm{sgn}(\mathrm{\&mu;}-0.5)T_k\&lsqb;{(1+\frac{1}{T_k})}^{|2\mathrm{\&mu;}-1|}-1\&rsqb;---\left(5\right)$

Calculate stochastic variable x, further according to following formula:

$V_i^{k+1}=V_i^k+x*(V_i^{\max }-V_i^{\min })---\left(4\right)$

Regulate the speed model, whereinWithIt is the minimax boundary value of the i-th interval velocity,For former rate pattern,For adjust after rate pattern, under this rate pattern, carry out shooting point reorientation；

E, travel through all volume elements, find Voice segment maximum E (V) under new rate pattern, and compare with E, if E (V) > E, then using new Energy maximum value E (V) as the object function in extremely fast simulated annealing method, if E (V) < E, then E (V) is with probability P alternative objective function.Wherein probability P computing formula is:

$P(V\&RightArrow;V^{\&prime;})=\exp \&lsqb;\frac{E\left(V\right)-E\left(V^{\&prime;}\right)}{T_k}\&rsqb;---\left(6\right)$

E (V) is former Energy maximum value, and E (V') is the energy value currently calculated.

After f, simulated annealing terminate, if perforation positioning precision meets the requirements, then show that rate pattern now is closer to actual geographic rate pattern, or Approximate Equivalent truly geographical rate pattern, it is possible to position for follow-up micro-seismic event.If perforation positioning precision is unsatisfactory for required precision, illustrate that the relatively actual geographic rate pattern error of rate pattern now is relatively big, then also need to change initial temperature T₀, carry out temper and regulate the speed model, until perforation position error is sufficiently small, i.e. the sufficiently close together actual geographic rate pattern of rate pattern.

Velocity model corrections flow process is as shown in Figure 1.

(b) and (c) in Fig. 3 is respectively adopted the result obtained after the inventive method calculates with first break pickup method, as can be seen from the figure adopt the rate pattern calculated addition of waveforms result that the inventive method builds better, Fig. 4 is that the rate pattern being respectively adopted after initial velocity model and the inventive method correction carries out the result of perforation location, it can be seen that the rate pattern after the inventive method corrects can obtain higher perforation positioning precision.

Method in order to verify the present invention still has good positioning precision when relatively low signal-to-noise ratio, adding signal to noise ratio in the data of Fig. 3 (a) is the Gaussian noise of S/N=0.1, add the generated data after noise as shown in Figure 5, Fig. 6 is respectively adopted the result obtained after the inventive method calculates with first break pickup method, as can be seen from the figure adopts the rate pattern calculated addition of waveforms result that the inventive method builds better.It is respectively adopted the rate pattern after the inventive method correction and initial velocity model carries out perforation location, positioning result is as shown in Figure 7, as can be seen from the figure the inventive method has good noise resisting ability, it is possible to obtain good positioning result when low signal-to-noise ratio.More than test result indicate that, adopting inventive algorithm can be effectively improved perforation positioning precision, it is more high that perforation is decided to be precision, illustrates that the rate pattern after correction is closer to the geographical rate pattern of reality.Good rate pattern is that follow-up microseism location provides foundation more accurately.

Claims (3)

1. the ground observation microseism velocity model corrections method based on amplitude superposition, it is characterised in that comprise the following steps:

A, centered by shooting point, set up three-dimensional target region, select reference channel M；

B, arranging initial temperature T0, minimum temperature Tmin, dwell time T, simulated annealing calculates parameter；

D, set up initial velocity model, definition E (V)=0；

E, obtain initial velocity vector V, read perforation data；

F, calculate each road travel-time difference relative to library track；

G, by each track data translate superposition, obtain the Energy maximum value E under existing rate pattern and coordinate thereof；

H, the value of E is assigned to E (V)；

The coordinate of i, E (V) reaches T dead time；

J, energy are that the coordinate at (V) place is as the perforation elements of a fix；

K, whether meet perforation positioning precision, be；

L, end.

2. the ground observation microseism velocity model corrections method based on amplitude superposition described in claim 1, it is characterized in that, simulated annealing described in step b calculates parameter setting, is the feature according to ground array microseism observation, arranges six simulated annealings and calculates parameter:

First parameter is velocity vector V, V=[V_p1,V_p2,V_p3,...,V_pn],V_pi

It is i-th layer of p wave interval velocity, just obtains initial velocity vector V according to well-log information；

Second parameter is object function E (V), in log data, selects, arbitrarily together as library track M, to adopt the method for ray tracing to ask for other each road theoretical travel-time difference relative to this library track:

Δt^cal=[t₁-t_M,t₂-t_M,...,t_N-t_M](1)

According to this theory travel-time difference, the perforation wave recording that each road cymoscope obtains being carried out reverse-time migration superposition, stack power (object function) mathematic(al) representation is as follows:

$E\left(V\right)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\&lsqb;\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}A(i,j)\&rsqb;---\left(2\right)$

Wherein A is the i-th volume elements amplitude of vibration size in the jth moment, and the center point coordinate of i-th volume elements is (x_i,y_i,z_i), M is cymoscope number, and N is time window length；

3rd parameter is simulated annealing initial temperature T₀, the solution obtaining initial temperature is as follows, first give initial temperature one only small on the occasion of, be then constantly multiplied by the perseverance number more than 1, until the satisfied acceptance probability to any model is close to 1；

4th parameter is annealing cooling parameter, extremely fast simulated annealing cooling profiles such as following formula:

T_k=T₀exp(-ck^1/2N)(3)

Wherein k is iterations, T₀For initial annealing temperature, c is used for adjustment algorithm annealing temperature for given constant, allows c=0.5 here；N regulates the speed the number of plies for needs；

5th parameter is stochastic variable x, and be used for regulating the speed vector, specifically gives formula as follows:

V_i ^k+1=V_i ^k+x*(V_i ^max-V_i ^min)(4)

V_i ^maxWith V_i ^minIt is the minimax boundary value of the i-th interval velocity, whereinX is stochastic variable, produces x expression formula as follows:

$x=\mathrm{sgn}(\mathrm{\&mu;}-0.5)T_k\&lsqb;{(1+\frac{1}{T_k})}^{\left|2\mathrm{\&mu;}-1\right|}-1\&rsqb;---\left(5\right)$

Wherein sgn is sign function, and x span is between [-1,1]；

6th parameter is acceptance probability, and as E (V ') >=E (V), V ' substitutes V as currently most solution, as E (V ') < E (V), with probability

$P(V\&RightArrow;V^{\&prime;})=\exp \&lsqb;\frac{E\left(V\right)-E\left(V^{\&prime;}\right)}{T_k}\&rsqb;---\left(6\right)$

Substitute currently most solution, wherein T_kTemperature value during for kth time iteration.

3. the ground observation microseism velocity model corrections method based on amplitude superposition described in claim 2, it is characterised in that simulated annealing calculates parameter iteration end condition and is:

A, Simulated annealing are reduced to lowest set temperature T_min；

B, shooting point positioning precision reach setting value；

C, when shooting point place energy value reaches certain value, through long period T iterative computation still without being substituted, by adjust above-mentioned parameter, to reach the balance of computational accuracy and computational efficiency.