Summary of the invention
The object of the present invention is to provide a kind of evaluation can implement the method for time shift earthquake.
Method provided by the invention may further comprise the steps:
Obtain the energy changing difference that oil deposit parameter changes the earthquake useful signal that causes; Obtain the noise level in the actual seismic data; More said energy changing difference and said noise level, if said energy changing difference greater than said noise level, explanation can be implemented the time shift earthquake;
Wherein: obtain the energy changing difference that oil deposit parameter changes the earthquake useful signal that causes and comprise following 1)-3) step:
1) carries out following step respectively at two time points and obtain geological data twice: the stratum basic parameter of measuring the target area; Make the oil reservoir section of said target area, set up the stratigraphic section module of said target area according to said oil reservoir section; Obtain the geological data 1a of said target area according to said stratum basic parameter and said stratigraphic section module) or 1b): 1a) two of velocity of longitudinal wave and density data volumes; 1b) three of velocity of longitudinal wave, shear wave velocity and density data volumes;
2) utilize ACOUSTIC WAVE EQUATION, twice geological data in the step 1) is scaled the seismic response record twice;
3) utilize following formula 5 with step 2) in twice seismic response record be scaled the energy changing difference that said oil deposit parameter changes the earthquake useful signal that causes;
In the formula 5, A
S1Be the primary seismic response record of simulation, A
S2Be secondary seismic response record of simulation, summation in window when ∑ is represented reservoir, K
oBe the energy changing difference that oil deposit parameter changes the earthquake useful signal that causes;
Wherein, the noise level that obtains in the actual seismic data comprises following step a)-c):
A) gather actual seismic data (earthquake wave amplitude)
B) data that obtains to step a) is handled and is obtained the two-dimension earthquake section;
C) through the SVD filtering method two-dimension earthquake section that step b) obtains is handled, calculated the useful signal section and the noise sections of this real seismic record, utilize formula 14 to calculate the noise level in the actual seismic data then;
In the formula 14, A
nBe isolated noise record from actual seismic section, A
rBe real seismic record, summation in window when ∑ is represented reservoir, K
nIt is the noise level of seismic data.
The relation of the geological data of above-mentioned steps actual seismic data and top step 1) a): be that step 1) is based on the subsurface deposit section and sets up the geological geophysical model; Obtain seismologic record through forward simulation, thus the contact between foundation and the actual seismic data.
Above-mentioned steps 1) in, said stratum basic parameter can be factor of porosity, pressure, WS, gas saturation and/or oil saturation.The application of method in reservoir monitoring that can above-mentioned evaluation implement the time shift earthquake also belongs within protection scope of the present invention.
Change the result's contrast with Analysis SNR through analog energy, can find out among the embodiment the gas field since the energy changing (13.3%) of the earthquake useful signal that the variation of degradation oil deposit parameter causes under the pressure far above the noise level (2.9%) of this area.Therefore think that this gas field can implement the time shift earthquake at present.The variation of time shift seismic energy is just thought and can be implemented the time shift earthquake greater than noise level.One of innovation of the present invention just is the time shift seismic energy change modeling based on reservoir model is come out; And the noise level of this energy changing and actual seismic data connected, the time shift earthquake feasibility analysis method that forms thus will provide important reference for everybody time shift earthquake feasibility analysis work from now on.
Embodiment
Below in conjunction with specific embodiment the present invention is described further, but the present invention is not limited to following examples.
Among the following embodiment,, be conventional method like no specified otherwise.
Embodiment 1, time shift earthquake feasibility analysis
One, based on the forward simulation analysis of model
1. the foundation of geological geophysical model
When carrying out the earthquake numerical simulation, must at first set up the geological geophysical model.Accurate geological-geophysical model is the basis of earthquake numerical simulation.
1) sets up the stratigraphic section module
Set up research block stratigraphic section module (Figure 1B) with femlab software according to the oil reservoir section (Figure 1A) of research block.Wherein the oil reservoir section is in the ECLIPSE numerical simulation software; To the target geologic model; Through numerical evaluation; Can export the grid property parameter (factor of porosity, permeability, saturation degree, pressure, the oil reservoir degree of depth etc.) of any development phase, arbitrary section, the sectional view that has the different attribute value is the oil reservoir sectional view.
2) input basic parameter
Measure the basic parameter (factor of porosity Por, pressure P p, WS Sw, gas saturation Sg, oil saturation So) on stratum to be measured; Be followed successively by each block of cells (like the block of setting up according to different colours among Figure 1A among Figure 1B) input basic parameter, set up in the good stratigraphic section module (Figure 1B) thereby be input to step 1).
Wherein the factor of porosity section is like Figure 1A (the factor of porosity section with this twice time in 2008 was the same in 2003), and Fig. 1 C is pressure traverse (last figure is 2003, and figure below is 2008), and Fig. 1 D is gas saturation section (last figure is 2003, and figure below is 2008)
3) well logging statistics
(1) according to log data, distinguish layer of sand and mud layer, count the relation of the shale index SH and the factor of porosity Por of layer of sand and mud layer respectively.
(2) for layer of sand, according to its fluid situations, applicating fluid replacement technology is the velocity of longitudinal wave Vp of sandstone and density p data conversion the velocity of longitudinal wave Vp of full water sandstone
SatAnd density p
Sat, count velocity of longitudinal wave Vp
SatAnd density p
SatRelation with factor of porosity Por; For mud layer, directly count the relation between its velocity of longitudinal wave and density and factor of porosity.
(3), use the shear wave data Vs that the mud stone formula obtains full water sandstone for layer of sand
SatWith velocity of longitudinal wave Vp
SatRelation; For mud layer, use the relation that the mud stone formula obtains rock shear wave data and velocity of longitudinal wave.
4) Model Calculation
For sandstone, the p-and s-wave velocity and the density (Vp of the shale index that step 3) is counted, full water sandstone
Sat, Vs
Sat, ρ
Sat) formula is input to and sets up in the good stratigraphic section module, applicating fluid replacement interpolation calculation obtain actual formation p-and s-wave velocity and density data (Vp, Vs, ρ); For mud stone, p-and s-wave velocity that step 3) counts and density are the actual formation data.(velocity of longitudinal wave such as Fig. 2 A, shear wave velocity such as Fig. 2 B, density such as Fig. 2 C).
Repeat the influence that above-mentioned steps considers that simultaneously pressure changes, can obtain the actual formation data (velocity of longitudinal wave, shear wave velocity and density) of (2008) for the second time.
2. the seismic wave field forward simulation calculates
Under the known situation of underground medium structural model and respective physical parameter,, and calculate on ground or the received seismic response record of underground each observation station based on the propagation law of ACOUSTIC WAVE EQUATION modeling effort seismic event in underground various media.
1) based on the pre-stack seismic analogy method of " grid method "
The correctly non-homogeneous influence with complex surfaces and interface shape of simulation medium of the entirely discrete method for numerical simulation of wave equation, so the value of particular importance is arranged in method of seismic prospecting is studied.The algorithm that the existing more rule-based grid of full discrete values analogy method is discrete is like limited method of difference, pseudo-spectrometry etc.Can use the discrete method of non-regular grid to comprise finite element method, spectral element method etc.Grid method is one type of new non-regular grid discrete method, and it is non-rule capable of using and unstructured grid meticulous depiction complex dielectrics interface and surface, with the discrete method of difference of regular grid suitable calculated amount is arranged again.
In conjunction with grid method; Absorbing boundary condition adopts the non-regular grid PML method that divides equation based on local coordinate system down; This absorbing boundary condition can promote based on the discrete seismic event analogy method of non-regular grid; Particularly the effect of grid method can pass through again to confirm artificial border neatly, reduces the zoning.
The core of sound wave grid method be set up under the local coordinate system the division ACOUSTIC WAVE EQUATION and based on the differential equation weak form of integral approach.
Non-homogeneous for treatment media better, the ACOUSTIC WAVE EQUATION that the sound wave grid method uses is following:
P is a pressure in the formula, and ρ and c are density of medium and speed.The core of grid method is under non-regular grid, provides the spatial spreading based on integral approach (weak) form of formula (1); For ease of the meticulous depiction complex structure, general non-rule, the unstructured grid that constitutes by triangle that adopt.With node k among Fig. 3 is example, at the k of dotted line neighbor domain of node formula (1) is made the area branch:
In the formula triangle number of m around node k, s
IlIt is the phantom line segments in each triangle.Utilize the lumped mass model in the dynamics calculation, promptly suppose 1/ ρ c along surface distributed
2Focused on each node, the area of formula (2) left end divides and can be approximately
Q
kBe ∫ ∫ in each triangle of
k vertex neighborhood
Ω1/ ρ c
2/ 3rd of a value summation of dxdz.Interpolating function in the numerical evaluation intermediate cam shape unit is always linear; Therefore
and
is constant in each triangular element, and available triangle difference operator is calculated based on the P value on the node.Based on above-mentioned 2 points, can obtain the spatial spreading based on integral approach (weak) form of formula (2):
Subscript 1 amount of showing is corresponding the 1st triangle around node k in the formula, b
kAnd a
kBe the geometric parameter of triangular element, its value be the corresponding length of side of this node on coordinate axis projection 1/2nd; With the triangular element ijk among Fig. 3 is example, b
k=(z
i-z
j)/2, a
k=(x
i-x
j)/2; D
xAnd D
zBe the triangle difference operator,, can be expressed as following simple form based on the geometric parameter of triangular element:
Subscript r represents leg-of-mutton three nodes in the formula, and A is leg-of-mutton area.
Given t is the P on each node constantly, can be tried to achieve the second time derivative of P on each node by (2) formula, based on time domain three central-difference formulas, can try to achieve the t+ Δ t P in the moment by t-Δ t and t P constantly, and this has just accomplished P renewal in time.Only utilized leg-of-mutton geometric parameter in the numerical evaluation, if calculated in advance and store these geometric parameters adopts rule and non-regular grid will not change calculated amount; And, only need store one group of geometric parameter, so required storage is also little to one group of triangle that equates.
2) based on the time shift seismic energy change calculations of simulating
The pressure information of underground wave field transfers amplitude information on ground to after wave detector receives.Because the energy propagated of ripple and wave amplitude square are directly proportional.So the quadratic sum of the amplitude of interior each sampling point of window has just been represented the relative size of the reflected energy on underground certain section stratum during seismic section.The difference that the time shift BEFORE AND AFTER EARTHQUAKE is gathered data for twice is caused by factors such as noise and changes in reservoir jointly.Therefore, need the difference of the seismologic record that changes in reservoir relatively causes and the relative size of noise energy.Poor according to the record of the seismic reflection before and after the actual production digital simulation changes in reservoir, this difference does not contain the noise influence, its energy (the following amplitude quadratic sum that all refers to) and the changes in reservoir of simulation preceding to the energy ratio K in the seasonable window
oReflected that changes in reservoir causes the relative size of energy changing.Concrete account form is following:
Wherein, A
S1Be the seismologic record of (2003) before the changes in reservoir of simulating, A
S2Be the record of (2008) after the changes in reservoir of simulating, summation in window when ∑ is represented reservoir.
The result obtains before and after the oil reservoir development twice pre-stack seismic single shot record (Fig. 4 A and Fig. 4 B), and the difference (Fig. 4 C) of the corresponding seismologic record that is caused by changes in reservoir.By the seismic signal capacity volume variance of objective interval, the number percent that can obtain energy changing is 13.3%.
Two, the Analysis SNR of actual seismic data
1) SVD (svd) denoising method:
SVD is a kind of new denoising method of introducing in recent years in the seismic data processing, and it has obtained widespread use in seismic data is handled.This method can keep the coherent signal of a certain direction (like horizontal direction), suppresses the ripple of irrelevant noise and other direction.Its advantage is not require that data are uniform samplings; Lose the dynamic characteristic of signal hardly.
Concrete realization principle is following:
If seismologic record has M road, there is N sampled point in each road, and the data of this seismologic record can be represented with a matrix X.
X=[x
ij],i=1,2,…,M;j=1,2,…,N (6)
Wherein i representes Taoist monastic name, and j representes the sampled point sequence number.According to matrix theory, any one M * N matrix can be made quadrature and decompose, that is:
X=U∑V
T (7)
U and V are respectively M * M and N * N rank orthogonal matrix, and ∑ is M * N matrix.If all row of matrix U are by covariance matrix XX
TEach proper vector (M dimension) form, V is by X
TEach proper vector of X (N dimension) is formed, and then ∑ is M * N diagonal matrix, and main diagonal element is by the singular value σ of X
i(XX
TOr X
TI the eigenvalue of X
iNon-negative square root) form, off-diagonal element is zero entirely, claims that at this moment (7) formula is the svd of matrix X.
(7) formula also can be write as
(8)
Here r is the order of matrix X, u
iBe XX
TCorresponding to σ
iProper vector, v
iBe X
TX is corresponding to σ
iProper vector.
is a M * N matrix, is called i the proper vector of X.Obviously, the svd formula of matrix X is not unique.For the convenience on using, the diagonal element of our regulation diagonal matrix ∑ is by the series arrangement of successively decreasing, that is:
σ
1>σ
2>…>σ
r>0
When the M road is separate, and r=min (M, N), generally speaking, r≤min (M, N).The Frobenius norm of matrix X is defined as:
Get F and equal 2, following formula becomes
Know that by (10) formula the F norm of matrix (when F=2) has in fact reflected the gross energy of M road earthquake record.Can prove,
(11) formula shows, the gross energy of M road earthquake record equals covariance matrix eigenwert (or square of singular values of the X) sum of X.
If M<N, when M seismologic record road linearity is irrelevant, σ
j(j=1,2 ..., M) all non-vanishing.So the complete reconstruct of X needs all proper vectors.When M seismologic record road of X linear when relevant, σ only then
1Non-vanishing, the complete reconstruct of X only needs first proper vector.The seismic trace correlativity is good more, and it is just few more that reconstruct should be write down needed proper vector.
If P proper vector approximate reconstruction X before only using, the seismologic record of reconstruction is X
P, then:
X
PGross energy be:
Through the SVD filtering method two-dimension earthquake section (Fig. 5 A) (seismic section is through seismic data acquisition, and processing obtains then) is handled.
2) Analysis SNR of actual seismic data
Through aforementioned calculation, can draw the useful signal section (Fig. 5 B) and the noise sections (Fig. 5 C) of this real seismic record.From actual seismic data, isolate noise, the ratio K of interior noise energy of window and earthquake raw readings energy in the time of just zone of interest can being calculated
n, also just reflected the noise level of this seismic data.Computing method are following:
A
nBe isolated noise record (noise sections of real seismic record just) from actual seismic section, A
rBe real seismic record (the useful signal section of real seismic record), summation in window when ∑ is represented reservoir.The noise level that can obtain actual seismic data objective interval according to said method is 2.9%.
Three, time shift earthquake feasibility analysis
Calculate when can obtain the objective interval place in the window useful signal percentage change that causes by changes in reservoir through forward simulation.Can obtain the noise level (number percent) of objective interval through Analysis SNR.
With objective interval should the time noise level in the window compare with the useful signal percentage change, if the useful signal percentage change is higher than noise level, just thinks so theoretically and can implement the time shift earthquake.Certainly, the variation of effective energy is compared with noise level, more remarkable just suitable more time shift earthquake.
Through result's contrast of analog energy variation and Analysis SNR, can find out that this gas field is owing to degradation oil deposit parameter under the pressure changes the noise level (2.9%) of the energy changing (13.3%) of the earthquake useful signal that causes far above this area in the present embodiment.Therefore think that this gas field can implement the time shift earthquake at present.